A three-dimensional numerical analysis of moisture flow in wood and of the wood’s hygro-mechanical and visco-elastic behaviour

Abstract

A three-dimensional numerical model was employed in simulating nonlinear transient moisture flow in wood and the wood’s hygro-mechanical and visco-elastic behaviour under such conditions. The model was developed using the finite element software Abaqus FEA^®, while taking account of the fibre orientation of the wood. The purpose of the study was to assess the ability of the model to simulate the response of wood beams to bending and to the climate of northern Europe. Four-point bending tests of small and clear wood specimens exposed to a constant temperature and to systematic changes in relative humidity were conducted to calibrate the numerical model. A validation of the model was then performed on the basis of a four-point bending test of solid timber beams subjected to natural climatic conditions but sheltered from the direct effects of rain, wind and sunlight. The three-dimensional character of the model enabled a full analysis of the effects of changes in moisture content and in fibre orientation on stress developments in the wood. The results obtained showed a clear distinction between the effects of moisture on the stress developments caused by mechanical loads and the stress developments caused solely by changes in climate. The changes in moisture that occurred were found to have the strongest effect on the stress state that developed in areas in which the tangential direction of the material was aligned with the exchange surface of the beams. Such areas were found to be exposed to high-tension stress during drying and to stress reversal brought about by the uneven drying and shrinkage differences that developed between the outer surface and the inner sections of the beams.

Introduction

The renewable character and the low carbon footprint of wood as compared with for example concrete (Dodoo 2011), have made wood increasingly popular for use in structural applications (Malo et al. 2016). Due to the natural process of sorption that takes place, wood interacts very readily with the surrounding climate, trying to establish an equilibrium moisture content (EMC), even when the direct effects of rain, solar radiation or wind are lacking (Gustafsson et al. 1998). When wood is subjected to a combination of change in moisture content (MC) and state of stress brought on by external load (e.g., mechanical load) or internal constraint (e.g., through differential shrinkage or swelling), a continuous change in the level of stress and the occurrence of deformations can be the result. In bending, an unfavourable state of stress can lead to an exponential increase in deflection and, in consequence of this, to failure (Armstrong and Kingston 1962; Bodig and Jayne 1982) or to fracture in the longitudinal-radial plane (TL- and TR-oriented cracking) that can weaken the shearing capacity of beams (Larsen 2013; McMillen 1958). The techniques available for monitoring changes in MC and hygro-mechanical viscoelastic behaviour in-situ can only be employed in specific locations in or around the wood (Fortino et al. 2019a, b). With recent developments in three-dimensional modelling, predictions can be made both of moisture flow and of moisture-induced stress and distortion in three-dimensional space.

Many studies regarding the impact of moisture gradients on beam bending in which simulations have been involved have been carried out. These have concerned both controlled changes in climate (Honfi et al. 2014; Ranta-Maunus 2001) and changes in climate that occur naturally (Fortino et al. 2019b; Franke et al. 2019; Müller et al. 2007; Toratti 1992). A variety of constitutive models has been used in such studies, involving one-, two- or three-dimensional space considerations. Yet, the fibre orientations that can be present in three-dimensional space (such as annual ring patterns, conical shapes or spiral grain formations) have thus far never been fully incorporated into such studies, despite the strong effects these can have on the stress distribution and thus on the deformations that the wood can be subjected to (Angst-Nicollier 2012; Ormarsson et al. 1999). Spiral graining is a phenomenon in which the wood fibres spiral around the pith that lies at the centre of the log. The magnitude of the spiral graining at each point in the material can be expressed as an angle $\theta$ in the $lt$-plane (Ormarsson et al. 1998). Between the pith and the bark, $\theta$ can vary between a value of $+ 8$ and $- 5^\circ$ over a distance of some 200 mm (Dahlblom et al. 1999; Säll 2002). Spiral graining can strongly affect both the deflection to which the beams are subjected and the stress fields involved (Ormarsson 1999; Schniewind and Barrett 1972). Together with changes in moisture content, this can lead to complex stress states that can change continuously.

It was shown recently that use of a nonlinear single-Fickian model in connection with nonlinear Neumann boundary conditions is able to describe not only slow MC changes (Konopka et al. 2018), but also quick MC changes (Florisson et al. 2020). One challenge that brings this model about, is that there is not enough experimental data available to cover the moisture and temperature dependent diffusion coefficient (Avramidis and Siau1987; Rosenkilde 2002) and surface emission coefficient (Yeo et al. 2002, 2008; Yeo and Smith 2005). Most of the additional data available in this area concern only a specific temperature or moisture content, though often in a wide variety of wood species, only certain material directions, or only simply for sapwood and heartwood. An overview of various sets of data dealt with in the literature concerning the moisture dependency of the coefficients is presented in Fig. 1.

Fig. 1

Moisture- $w$ and temperature- $T$ dependent radial and tangential diffusion coefficient $D_{r,t}$, as well as surface emission coefficient $s,$ as reported in different literature sources. hwd refers to heartwood and swd to sapwood

For a moisture-sensitive and visco-elastic material such as wood, it is common to describe the total strain rate as a summation of the elastic, hygro-expansion, visco-elastic creep (in this paper regarded as creep) and mechano-sorptive strain rates (Ormarsson et al. 1998). The elastic component is defined according to Hook’s law and the hygro-expansion as the relationship between the hygro-expansion coefficient and the rate of change in moisture content. The most conventional way of describing mechano-sorption is with use of a model presented by Salin (1992). This is a model describing both the decreasing deflection rate phase and the constant deflection rate phase, two phases of mechano-sorption deflection that usually take place under moderate levels of stress, or in other words low-level bending (Florisson et al. 2021). The most commonly used model for describing creep is the hereditary model (Dahlblom 1987). Similar to the mechano-sorption model, this approach is used to describe not only the decreasing deflection rate and the constant deflection rate, but also the unloading response of the material.

One difficulty in connection with three-dimensional modelling of moisture-induced distortion is the vast amount of experimental data needed to describe the material parameters that are involved. For Norway spruce, sufficient experimental data are available to cover both the moisture- and the temperature-dependent elastic and hygro-expansive behaviour that can occur (see Table 1). However, for mechano-sorption the data are limited to tension and compression in the principal material directions $l$, $r$ and $t$ (Aicher et al. 1998; Hunt 1984; Kangas and Ranta-Maunus 1989; Ranta-Maunus 1993; Toratti et al. 2000). To the author’s knowledge, not much has been reported on shear, although the work in this area by Hearmon and Paton (1964) can be cited. A small dataset concerned with mechano-sorption in the $rt$-plane was studied by Kangas and Ranta-Maunus (1989). Experimental data on creep are available, both in the principal material directions and in the $lr$- and $lt$-planes, yet it is scattered in the sense of only certain wood species being involved (Hayashi et al. 1993; Ozyhar et al. 2013; Schniewind and Barrett 1972; Svensson 1995, 1996), and often only short-term test periods $( < 1\; {\text{day}})$ being dealt with, so that there is a lack of predictors of long-term behaviour $( > 50 \;{\text{days}})$ when exponential expressions are used to describe creep deflection (Florisson et al. 2021).

Table 1 Inventory of literature on the elastic, hygro-mechanical, creep and mechano-sorption properties of Norway spruce, their moisture- and temperature-dependency, and the constitutive models used in the present paper to explain such properties

A three-dimensional numerical model was created in the FE-software Abaqus FEA® to simulate both the transient nonlinear moisture flow and the moisture-dependent distortion and stress that normally takes place in wood, account being taken of the fibre orientation. To determine whether the model was able to simulate in an adequate way the beam bending that occurs under northern European climatic conditions, the following steps were taken: (1) on the basis of experimental data available in the literature, a set of expressions was created to determine the moisture- and temperature-dependent diffusion coefficient and surface emission coefficient that appeared to be appropriate, (2) experimental results obtained for small beams tested under constant temperature and systematic RH-cycle (controlled climatic) conditions were used to calibrate an appropriate numerical model, account being taken of the spiral grain that applied and the annual ring curvature, and (3) test results for solid beams tested under natural climatic conditions were used to validate the numerical model, account being taken of the fibre orientation. Here, calibration involved the iterative process of changing material parameters until the simulation obtained agreed closely with the experimental data. This process indicated the order of range of the hygro-mechanical and the visco-elastic material parameters. Validation involved the assessment of whether the numerical model accurately described reality.

Materials and methods

Theory

Moisture flow model

The constitutive relationship used to describe the nonlinear transient moisture flow in wood was defined as

$$ \overline{\varvec{q}} = - \overline{\varvec{D}}\left( {w,T} \right)\overline{\varvec{\nabla }}w $$

(1)

where ${\varvec{q}}$ is the moisture flux vector, ${\varvec{D}}$ is a moisture- and temperature-dependent diffusion matrix, $w$ is the unit-less moisture content, $T$ is the temperature and ${\varvec{\nabla }}$ defines gradients. The over-line $\left( {\overline{\varvec{ \cdot } }} \right)$ in Eq. (1) refers to the orthogonal local coordinate system for the wood. The Neumann boundary condition that describes the flux that is normal for the exchange surface $q_{n}$ is defined as

$$ q_{n} = s\left( {w,T} \right)\left( {w - w_{\infty } } \right){ } $$

(2)

where $s$ is a moisture- and temperature-dependent surface emission coefficient and $w_{\infty }$ is the unit-less moisture content of the ambient air. Further information concerning derivation of the FE-formulation and how to iteratively solve the nonlinear system of equations here can be found in Florisson et al. (2020) and Johannesson (2019).

Stress model

The moisture-related total strain rate (column matrix) $\varvec{\dot{\bar{\varepsilon }}}$ is described by the constitutive equation presented in Eq. (3) as a summation of the elastic strain rate $\varvec{\dot{\overline{\varepsilon }}}_{{\text{e}}}$, the strain associated with hygro-expansion $\varvec{\dot{\overline{\varepsilon }}}_{{\text{h}}}$, the mechano-sorption strain rate $\varvec{\dot{\overline{\varepsilon }}}_{{{\text{ms}}}}$, and the creep strain rate $\varvec{\dot{\overline{\varepsilon }}}_{{\text{c}}}$. The notation ${\dot{(\varvec\cdot)}}$ refers to the time derivative $\partial /\partial t$ of the strains. The different matrix expressions used to calculate each strain rate in the constitutive model (see Eq. (3)) are presented in Eqs. (4) to (7)

$$ \varvec{\dot{\overline{\varepsilon }}} = \varvec{\dot{\overline{\varepsilon }}}_{{\text{e}}} + \varvec{\dot{\overline{\varepsilon }}}_{{\text{h}}} + \varvec{\dot{\overline{\varepsilon }}}_{{{\text{ms}}}} + \varvec{\dot{\overline{\varepsilon }}}_{{\text{c}}} $$

(3)

$$ \varvec{\dot{\bar{\varepsilon }}}_{e} = \varvec{\dot{\bar{\varepsilon }}}_{{e,1}} + \varvec{\dot{\bar{\varepsilon }}}_{{e,2}} = \varvec{\bar{C}}\left( {w,T} \right){\dot{\bar{\varvec\sigma }}} + \varvec{\dot{\bar{C}}}\left( {w,T} \right){\bar{\varvec\sigma }} $$

(4)

$$ {\mathbf{\dot{\overline{\varepsilon }}}}_{{\text{h}}} = {{\varvec{\upalpha}}}\dot{w}_{a} $$

(5)

$$ {\mathbf{\dot{\overline{\varepsilon }}}}_{{{\text{ms}}}} = \overline{\varvec{m}}\left( T \right)\left( {\overline{\varvec{\sigma }} - \varvec{\overline{n}\overline{\varepsilon }}_{{{\text{ms}}}} } \right)\left| {\dot{w}_{a} } \right| $$

(6)

$$ {\mathbf{\dot{\overline{\varepsilon }}}}_{{\text{c}}} = \mathop \sum \limits_{k = 1}^{K} \frac{1}{{\tau_{k} }}e^{{ - \frac{t}{{\tau_{k} }}}} \mathop \int \limits_{0}^{t} e^{{\frac{{t^{\prime}}}{{\tau_{k} }}}} \overline{\varvec{C}}_{{c_{k} }} \left( {w,T} \right)\frac{{{\text{d}}{\overline{\mathbf{\varvec\sigma }}}\left( {t^{\prime}} \right)}}{{{\text{d}}t^{\prime}}}{\text{d}}t^{\prime} $$

(7)

where ${\varvec{C}}$ is the compliance matrix, ${\varvec{\sigma}}$ is the stress vector, ${\varvec{\alpha}}$ is a vector that contains the hygro-expansion coefficients, $\dot{w}_{a}$ is the moisture content rate below FSP, ${\varvec{m}}$ is the mechano-sorption property matrix, ${\varvec{n}}$ is the retardation matrix, $\tau_{k}$ is a value that indicates in which time range the different creep terms have an effect on the change of sum, $t^{\prime}$ is the time corresponding to a recent stress change $\Delta {\varvec{\sigma}}$, and ${\varvec{C}}_{{c_{k} }}$ is the creep compliance matrix. The creep compliance matrix consists of ${\varvec{C}}$ scaled with use of the creep parameter $\phi_{{\sigma_{l} }}^{k} \left( - \right)$. The notation $\left| \varvec\cdot \right|$ indicates an absolute value. The time integral that is part of Eq. (7) is evaluated from time $t$ to time $\Delta t$ $\left( {t + \Delta t} \right)$, since the value of the integral is known at time $t$ due to its storage in a so-called state variable (Dassault Systèmes 2017). The creep model is able to account for both recoverable and irrecoverable creep. More about the process of deriving the FE-formulation and of how to iteratively solve this set of equations can be found, for example, in Ormarsson (1999) and Ottosen and Ristinmaa (2005).

Numerical model

The three-dimensional numerical model employed used a sequential solution scheme. First, the transient moisture flow problem was solved, which was then used as input (in terms of moisture content nodal values) to a moisture-induced distortion analysis (Fig. 2). Both analyses used the same geometry and fibre orientation. The moisture flow model was supported by the material routine UMATHT and the boundary condition routine FILM (Florisson et al. 2020). The constitutive theory used for the distortion model was implemented in the user-subroutine UMAT. This routine was used very recently in Florisson et al. (2019) to simulate the drying behaviour of timber boards, and is extended here by the mechano-sorption model of Salin (1992) and the hereditary creep model of Dahlblom (1987). The orthotropic material orientation was determined by the user-subroutine ORIENT (Ormarsson et al. 1998), which describes the local coordinate system $\left( {l,r,t} \right)$ in relation to the global coordinate system $\left( {x,y,z} \right)$ being employed. In the present study, the numerical model was used in two different applications: (1) calibration of four-point bending behaviour under controlled climate conditions, and (2) validation of four-point bending behaviour under natural climate conditions. Both applications use experimental data from self-performed tests.

Fig. 2

Principal components of the numerical model developed in FE-software Abaqus FAE® and its applications

Experimental program

Application 1

Three beams (beams A1-SS/NN/EE) were tested in terms of four-point bending inside a climate chamber of type Vötsch VCR 4033/S (Vötsch Industrietchnik GmbH, Balingen, Germany). The beams were 30 mm in height, 15 mm in width, and 640 mm in length when testing began (Fig. 3). A $5$ mm thick steel plate was placed beneath the loads and above the supports, each with an area of $16$ by $55$ mm, which resulted in beam spans of $595$ mm. A detailed description of the experimental methodology, the setup and the results obtained can be found in Florisson et al. (2021). The experiment provided data on mass change and hygro-expansion over the height of the beam, information regarding both being obtained from an unloaded specimen. In addition, deflection from the constant moment area $u_{y,1}$ was collected; see location LVDT in Fig. 3.

Fig. 3

Details concerning the four-point bending test setup for beam A1, carried out under controlled climate conditions: a close-up of the custom-made roller support, and b close-up of the custom-made pinned support

To show the effect of $\theta$ on calibration in the three long-term bending tests, two simulations were carried out for each specimen and for each of the deflection components (simulations of elastic creep and of mechano-sorption). To obtain the initial values, the first simulation was carried out with $\theta$ equal to $0^\circ$. The only contributing material parameters were those related to the $l$-direction. The second simulation was carried out using a positive $\theta$ of $4^\circ$. The material parameters belonging to the $l$- and to the $t$-direction and to the $lt$-plane affected the fit and required an adjustment once again to meet $u_{y,1}$. The procedure is explained in greater detail in the section on numerical simulations.

Application 2

Three solid timber beams (each of them an A2 beam) were tested in terms of four-point bending and were subjected to natural climate conditions inside a research facility located in the Asa Research Park in Southern Sweden. The facility was not isolated, but it protected the test material against the direct effects of rain, wind and sunlight. The beams were 195 mm in height, 45 mm in width, and 3100 mm in length, with a ±${12}\;\%$ moisture content (Fig. 4). The steel contact surface between the load points and the beam, or between the supports and the beam, had an area 45 mm by 100 mm in size. The beam spanned 3000 mm, which was divided by the location of the load into 3 equal segments, each 1000 mm in length.

Fig. 4

Details of the test setup for the four-point bending test, performed under natural climate conditions, in the case of beam A2, where a is a close-up of the roller support used for both the supports and for the loading points, and b is a close-up of the potentiometer, located on the bearing frame (not shown)

In a room adjacent to the setup, the temperature, $T$ and the humidity, RH, $\phi$ were monitored by use of a pin-type SHT75 temperature (± 0.3 °C) and humidity (${1.8}\;\%$ RH) sensor (Sensirion, Stäfa, Switzerland). The deflection $u_{y,2}$ of the wood beams was measured by an S18FLPA100 potentiometer (Sakae, Kawasaki, Japan) that had a measuring range of 100 mm $\left( {{ 0.3 \,\text{mm}}}\right)$ at mid-span 20 mm from the lower surface. The SAAB group (SAAB AB, Stockholm, Sweden) in Växjö developed the data collection system. $T$, $\phi$ and $u_{y,2}$ were logged every 10 $\text{min}$.

The $T$ and $\phi$ values that were recorded were used as input to the numerical model, whereas the average $u_{y,2}$ value of the three beams that were tested was used to carry out validation of beam A2. The effect of beam size on moisture transport and mechano-sorption was analysed by performing the simulation under natural climate conditions, also on the smaller geometry of beam A1, hereinafter referred to as beam A3. Beam A3 was assigned the same material parameters, the same expression for determining $\theta$ (the beam being located 130 mm from the pith), and the same bending stress and climate as beam A2. The deflection $u_{y,2}$ was analysed at mid-span at a level 3 mm above the lower exchange surface, at a position corresponding to the location at which the data for beam A2 was analysed.

Numerical simulations

The numerical model was used in two different applications that differed in terms of geometry, fibre orientation, element mesh, boundary conditions and load configurations. These different aspects of the situation involved will be discussed here, together with an assessment of the instantaneous elastic stress state. For both applications, element type DC3D8 was used to create the mesh for the moisture flow analysis, and element type C3D8 was used for the distortion analysis. A comparison with a quadratic stress element showed acceptable agreement between results.

Application 1

Application 1 concerned calibration of the numerical model on the basis of the experimentally obtained bending behaviour of three beams (A1-SS/NN/EE) that were tested under controlled climate conditions. The geometry of the model, as well as the element mesh, the boundary conditions and the load configuration, are shown in Fig. 5. The input parameters for the moisture flow analysis and its calibration are to be found in Florisson et al. (2020). The centre of the beams was located 137.5 mm from the pith, this resulting in a slight curvature of the annual rings in cross section. The beams also had a positive spiral grain angle of $4^\circ$. In performing the calibration, the differences in the deflection data between nodes $p_{1}$ and $p_{2}$ $\left( {p_{1} - p_{2} }\right)$ were determined and were compared with the experimentally obtained $u_{y,1}$ values (Fig. 5). Here, $p_{1}$ and $p_{2}$ indicate the position of the head of the LVDT and the right-hand side of the yoke, respectively, as shown in Fig. 3. The numerical analysis was carried out over a period of ${90\,\text{days}}$.

Fig. 5

Geometry, mesh, boundary conditions and material orientation used in the simulation of beam A1, $F$ indicating load and $A$ the area of the steel plate

The simulation began after a $5$-day period (when the value of $t$ was $0$) of constant climate conditions (60 °C and ${40}\;\%$ RH), after which conditions of pure elastic and creep deflection were established (and maintained for a period of $5$ days) (Fig. 6b, c), followed by $85$ days of testing at a constant $T$ of 60 °C and with systematic $\phi$-phases of between $40$ and ${80}\;\%$ for generating deflection caused by mechano-sorption. The experimentally obtained EMC at the end of each $\phi$-cycle was used to define $w_{\infty }$ in the subroutine FILM, which is compared in Fig. 6d with values determined by use of Simpson’s formula (Simpson 1971). This latter expression concerns the relationship between EMC and $\phi$ as well as $T,$ without hysteresis being taken into account.

Fig. 6

Schedules used for a load ${ }F$, b relative humidity ${ }\phi$, c temperature $T$ and d the moisture content of the ambient air $w_{\infty }$, in which d, the experimentally obtained $w_{\infty }$ (black), is used as a boundary condition in the moisture analysis and as compared with values obtained with use of Simpson’s formula (grey)

The total load $F$ of 413 N was applied instantaneously and remained constant during the entire analysis (Fig. 6a). Simple hand calculations based on elastic beam theory indicated that a maximum normal stress $\sigma_{x}$ of $18.6\; {\text{MPa}}$ could be expected for the outer fibres at mid-span, and a maximum shear stress $\tau_{xz}$ of 0.69 MPa and a stress perpendicular to the beam direction $\sigma_{z}$ beneath the loads and supports of 0.25 MPa could be expected. In Fig. 7, these values are compared with results of an elastic analysis performed in the global coordinate system. Rigid shells were used to represent the steel plates above the supports and below the loads. The interaction between the steel plates and the wood was defined as involving a combination of normal (hard) friction and tangential friction (isotropic static friction with a friction coefficient of ${0.3} \;\left( - \right)$ behaviour.

Fig. 7

a Stress variation of $\sigma_{x}$, $\sigma_{z}$, and $\tau_{xz}$ directly after the loading of beam A1-SS, and b variation in $\sigma_{x}$ and $\tau_{xz}$ over the height of the beam $h$ at a quarter and half of the span in bending $l_{s}$. Note that the magnification of deflection in the z-direction is $10\times .$

Figure 7 indicates differences between the stress values obtained by hand calculations and those obtained by simulation. The load $F/2A$ applied to the top of the steel plates, where $A$ is the contact area between the load and the steel, was found to result in a compression stress $\sigma_{z}$ beneath the plates. This stress was found to vary between quite significant beneath the loading points to being zero on the bottom surface of the beam, resulting in an asymmetric distribution of $\sigma_{x}$ over the height of the beam. A similar stress situation was found to exist above the supports, but due to the asymmetrical division of the steel plate by the location of the line supports, the stress is concentrated at the ends of the beam. Simulation of the long-term behaviour was carried out without reference to the steel plates. The boundary conditions, in contrast, were directly associated with the wood surface.

Material parameters

The moisture flow analysis was performed with use of the $w$- and $T$-dependent ${\varvec{D}}$ and $s$ found in Florisson et al. (2020). Plots of $D_{r}$, $D_{t}$ and $s$ are shown in Fig. 1. The dependency of the elastic moduli $E$ and the shear moduli $G$ on $T$ and $w$, and of the mechano-sorption moduli $m$ on $T$ was described using expressions presented by Ormarsson et al. (1998). The reference values needed for these expressions are indicated by $\left( \varvec\cdot \right)_{0}$ (Table 2). The reference values designated by $\left( \times \right)$ in Table 2 were calibrated by use of the experimentally obtained $u_{y,1}$ value.

Table 2 Hygro-mechanical (reference) material parameters used in the stress analysis, and the parameters to be determined by calibration ( ×)

The expressions describing the remaining elastic and shear moduli ($r$, $t$, $lr$, $lt$, $rt$) were calibrated using experimental data from Keunecke et al. (2007). The hygro-expansion coefficients are from Bengtsson (2001), except for $\alpha_{t}$, which is from Florisson et al. (2021). Both $E_{l}$ and the hygro-expansion coefficient $\alpha_{l}$ are assumed to be constant in the radial direction. In view of the lack of creep deflection data reported in the literature, it was assumed that the creep parameters that pertained to the $t$-direction and $lt$-plane were consistent with the relationship obtained by use of Eq. (7) and data reported by Schniewind and Barrett (1972). The mechano-sorption parameters stem from Ormarsson (1999), with the exception of the $n_{lt}$ parameter, which was made consistent with the relationships of the other parameters.

Application 2

The second application concerned validation of the numerical model used to simulate the bending behaviour of beam A2 under natural climate conditions. The geometry of the model, the element mesh, the boundary conditions and the load configuration are shown in Fig. 8. At both ends of beam A2, the pith of the tree was located in the centre of the beam’s height $h$ on the right exchange surface (see Fig. 8). The conical shape of the tree was taken to be $- 0.5^\circ$ and the spiral grain angle to be $\theta = 4 - 40r$, $r$ being the radial distance (in ${\text{m}}$) between the pith and the material point (Ormarsson 1999). The numerical analysis concerned a period of ${556\,\text \quad {days}}$ altogether.

Fig. 8

Geometry, mesh, boundary conditions and material orientation used in the simulation of beam A2, $F$ indicating load, $A$ being the area of the steel plate, $h$ the height of the beam and $r$ the distance between pith and analyzed material point

Figure 9 presents the load and the climate schedules used in simulating beam A2. Figure 9b and c presents the experimentally obtained raw $\phi$ and $T$ data, together with a 2-week average used as input to the numerical model. Subroutine FILM uses Simpson’s formula to determine $w_{\infty }$ (Fig. 9d) by use of the averaged $\phi$ and $T$ values.

Fig. 9

Schedules used for a the load ${ }F$, b the relative humidity ${ }\phi$, c the temperature $T$, and d the moisture content of the ambient air $w_{\infty }$, where for $\phi$, $T$ and $w_{\infty }$ both the raw data (grey) and the averaged data (black) are plotted over a $2$-week time period

The beam carried a total mechanical load $F$ of 6900 N (see Fig. 9a), which resulted in a bending stress of 12 MPa in the outer fibres, a nominal stress of 0.77 MPa perpendicular to the steel plates beneath the loading points and the supports, and a shear stress of 0.59 MPa. In Fig. 10, these results can be compared with the results of an elastic analysis carried out directly after the loading. Similar arguments apply to explaining the differences that can be seen in the stress values obtained for beam A1 (see Fig. 7). Similar simplifications were also made regarding the boundary conditions in the long-term analysis that was carried out.

Fig. 10

a Stress variations of $\sigma_{x}$, $\sigma_{z}$, and $\tau_{xz}$ directly after loading in the case of beam A2 and b variations in $\sigma_{x}$ and $\tau_{xz}$ over the height of the beam $h$ at a quarter and half of the span in bending $l_{s}$. Note that the magnification of the deflection in the $z$-direction is $10 \times$

Material parameters

To describe the moisture flow under natural circumstances, exponential expressions were created for ${\varvec{D}}$ and $s$ with regards to $w$ and $T$ with the help of experimental data found in the literature (see Fig. 11). The experimental data used to determine ${\varvec{D}}$ were from Avramidis and Siau (1987) (in the r- and $t$-directions), from Rosenkilde (2002) (in the $l$-direction), and for ${{s}}$ from Yeo et al. (2008). The expressions are to be found in Eqs. (8) to (10), where $D_{l}$ is constant with regards to $w$, and—due to a lack of data—is assigned the same temperature dependency as $D_{r,t}$.

Fig. 11

Moisture- $\left( w \right)$ and temperature- ($T)$ dependent diffusion coefficient $D$ as related to $l$, $t$ and $r$ material directions, and the surface emission coefficient as presented in Eqs. (8) to (10)

The hygro-mechanical parameters used in the distortion analysis of beam A2 are presented in Table 3. Here too, similar to what was found for beam A1, the expressions of Ormarsson (1999) are adopted to describe the $w$ and $T$ dependency of $E$ and $G$, and the $T$ dependency of $m$. The reference value $E_{l0}$ and the values for $m_{l0}$ and $n_{l}$ were determined on the basis of $u_{y,2}$.

$$ D_{l} \left( T \right) = 2.7 \times 10^{ - 8} \cdot \left( {8e^{0.07T} } \right) $$

(8)

$$ D_{r,t} \left( {w,T} \right) = \left( {1.7 \times 10^{ - 9} e^{11w} + 9.3 \times 10^{ - 13} } \right) \cdot \left( {8e^{0.07T} } \right) $$

(9)

$$ s\left( {w,T} \right) = \left( {2.1 \times 10^{ - 4} e^{ - 13w} + 4.5 \times 10^{ - 6} } \right) \cdot \left( {0.9e^{0.04T} } \right) $$

(10)

Table 3 Hygro-mechanical (reference) material parameters employed in connection with the stress analysis

The variation of $\alpha_{l}$ over the cross section of the log is described by the expression presented in Table 3, with $r$ being the distance between the pith and the bark (Ormarsson et al. 1998). The values for $\alpha_{r}$ and $\alpha_{t}$ are from Bengtsson (2001). All the other parameters were those adopted for beam A1. The creep parameters related to the $r$-direction, and the $lr$- and $rt$-plane are assumed to involve the same relationships as those applying to $\phi_{{\sigma_{l} }}^{ }$, $\phi_{{\sigma_{t} }}^{ }$ and $\phi_{{\sigma_{lt} }}^{ }$. The mechano-sorption moduli and the retardation parameters for the same material direction and planes are from Ormarsson (1999).

Results and discussion

Application 1

Deflection

The experimentally obtained $u_{y,1}$ values (the grey line) presented in Fig. 12a show the usual characteristics observed in initially dry wood tested at low-level bending (Armstrong and Christensen 1961; Armstrong and Kingston 1962). The first adsorption (uptake of water molecules) and the subsequent desorption (release of adsorbed water molecules) were found to result in a substantial increase in $u_{y,1}$. When the waves observed between high and low $\phi$-cycles are neglected, the expected decreasing deflection rate phase and its subsequent constant deflection rate phase were observed. The results obtained using the numerical model (black line) show a good fit with these two phases. The model is unable, however, to predict the characterising waves.

Fig. 12

Results of calibration, where a is the fit with the experimentally obtained $u_{y,1}$ value during $90{\text{ days}}$ of testing and b is the fit with the experimentally obtained $u_{y,1}$ value during the first $5{\text{ days}}$ of testing. Note that b starts where the instantaneous elastic deflection $u_{e,i}$ ends

Figure 12b shows that the series of Kelvin modules and the creep parameters from Tables 2 and 4 used to describe $u_{c}$ result in a good fit with the experimental results in the given time range of $5$ days. The figure also shows the fit used to obtain the material parameters (see Table 3) needed for the second application. Due to the numerical limitations placed on the exponent by the FE-program, the value selected for $\tau$ must be large enough for long-term analyses.

Table 4 Material (reference) parameters obtained in calibrations made in simulations in which the spiral grain angle $\theta$ is and is not $\left( {\theta = 0} \right)$ included, as regards beams A1-SS, NN and EE

The effects of spiral grain on material parameters

Table 4 presents the material parameters obtained on the basis of calibrations made for beams A1-SS, NN and EE. The rotation of the longitudinal-tangential plane caused by the spiral grain results in an increase in the longitudinal elastic modulus $E_{l}$ as compared with the values obtained for the simulations that include $\theta$ as $0^\circ$. The values obtained for $E_{l}$ when $\theta$ is $4^\circ$ show relative differences (RD) of $16.2$, $19.7$ and ${41.4}\;\%$ (EE, NN and SS) as compared with $E_{l}$ when $\theta$ is $0^\circ$. Since the other elastic and shear moduli do not change during the calibration, there is a lesser difference for stiffer specimens.

Since creep is described in terms of creep compliance ${\varvec{C}}_{{c_{k} }}$, calibration of the curve when $\theta$ is $4^\circ$ is considered to be made with $E_{l}$ rather than with $E_{l,\theta = 0}$. The approach adopted results, except in the case of beam A1-SS, in a slight increase in $\phi_{{\sigma_{l} }}$, due to a spiral grain in all of the beams. In the case of mechano-sorption, a spiral grain leads to a decrease in $m_{l}$ and an increase in $n_{l}$. The values obtained for $m_{l}$ when $\theta$ is $4^\circ$ show RD values of $30.4$, $20.0$ and ${44.4}\;\%$ (EE, NN and SS) with respect to $m_{l}$ when $\theta$ is $0^\circ$. When $n_{l}$ is considered in combination with an $\theta$ of $4^\circ$, relative differences of $8.3$, $22.2$ and ${51.4}\;\%$ (EE, NN and SS) were found with respect to $n_{l}$ when $\theta$ is $0^\circ$.

Stress distribution

Under appropriate circumstances, a beam subjected both to bending and to climate change is prone to drying and to mechanically driven cracks. During the period of $90$ days, no cracks were found to have occurred. The colour plots in Fig. 13 show the stress distributions of $\sigma_{l}$ and $\tau_{lt}$ for beam A1-SS directly after loading, as well as the effects of $\theta$. The graphs that apply indicate the effects of climate on the various stress types. A spiral grain does not appear to have any particularly strong effect on $\sigma_{l}$ for the mesh element adjacent to the lower surface at mid-span $\left( \times \right)$. However, the change in moisture content after each change in RH-phase results in a spike in the stress, which is slightly greater when $\theta$ is equal to $4^\circ$. These spikes are the result of rapid changes in moisture content and of high moisture gradients close to the surface. With regard to the initial state, a significant increase, respectively, of ${33.4}\;\%$ in $\sigma_{l}$ can be observed at the first desorption when $\theta$ equals $4^\circ$. For a stress state close to the ultimate limit state, such a change may possibly lead to failure.

Fig. 13

Stress development a $\sigma_{l}$ and b $\tau_{lt}$ for beam A1-SS when $\theta$ is $0$ or $4^\circ$. Note that the graphs are taken from specific mesh elements (×), that the stress field is obtained directly after loading and is projected onto half of the respected beam, and that the deflection is magnified by a factor of 10

The shear stress distribution $\tau_{lt}$ is strongly affected by a spiral grain, which enters the so-called constant moment area when $\theta$ is $4^\circ$. The maximum positive $\tau_{lt}$ moves then from the area close to the support (see $\times$ for $\theta = 0^\circ$ in Fig. 13) to the lower surface at mid-span. This is caused by the translation of a much larger normal stress $\sigma_{x}$ to $\tau_{lt}$. This translation also results in a vast increase in the extreme stress state from a level of about 0.69 MPa to about 1.3 Mpa, and in a much stronger effect of moisture on $\tau_{lt}$ (a ${14.7}\;\%$ increase at the first desorption compared to the initial state). In EN 1995-1-1 (2004), the effect of moisture is incorporated into the design of timber elements as a reduction in strength by use of a modification coefficient, $k_{{{\text{mod}}}}$. Since the stress is a local effect and the peaks are affected by $\theta$, the question remains of whether such a coefficient is able to adequately describe these effects, and of whether the climate effects should not instead be regarded as a load.

In addition, after each change in $\phi$, a steep moisture gradient develops quickly at the exchange surface, and within a short period of time a $\sigma_{t}$ stress field involving either compression or tension appears, depending upon whether it is swelling or shrinkage that is involved (Fig. 14). As a counterreaction, the inner section of the beam develops a stress state opposite to that at the surface. Such moisture-induced stress states can result in surface or internal checks when the tensile stress exceeds the tangential strength, and under appropriate conditions, cracks or so-called splits can develop. For Norway spruce, the tangential strength at 60 °C of different moisture content levels varies between 1 and 3 MPa (Gustafsson 2003; Hanhijärvi 1998; Larsen 2013). For simulations in which $\theta$ either is (0.91 MPa) or is not (1.68 MPa) accounted for, the strength is not exceeded (see Fig. 14). This can be confirmed by a visual inspection of specimens that are tested.

Fig. 14

Stress development $\sigma_{t}$ for beam A1-SS when $\theta$ is $0$ or ${ }4^\circ$. Note that the graphs are taken from specific mesh elements ( ×), that the stress field is either obtained after 10.25 or 15.25 days and is projected onto half of the respective beam, and that the deflection is magnified by a factor of 10

Due to the continuous change in RH, a continuous change in stress can be noted in the graphs presented in Fig. 14 for the locations $\left( \times \right)$ that are selected, the surface and the centre alternating between a state of tension and of compression, depending upon whether the surface undergoes shrinkage or swelling. In following the stress peaks from the first to the last adsorption and desorption, the effects of hysteresis can be seen in the values for extreme stress moving closer and closer to zero. In addition, it was found, as shown in Fig. 14, that $\theta$ lowers the maximum $\sigma_{t}$ level but does lead to a more permanent state of stress. A matter not dealt with here, but interesting to consider, is that, due to drying and the fact of there being a spiral grain, possible mixed failure modes that can lead to brittle failure become accessible. Figure 14 shows, for example, for a situation where $\theta$ is $0^\circ ,$ a relatively high value for $\sigma_{t}$ because of drying occurring at the same location and time as $\tau_{lt,}$ due to shear (see Fig. 13), and in Fig. 13 due to a value for $\theta$ of $4^\circ$, $\sigma_{l}$ and $\tau_{lt}$ on the upper and lower surfaces coinciding.

Application 2

Moisture data

The moisture content of the ambient air, $w_{\infty } ,$ was determined on the basis of Simpson’s formula and the raw $T$ and $\phi$-data recorded in Asa, showing a minimum value of $7$ and a maximum value of ${29}\;\%$ (Fig. 15). The three service classes used by engineers in designing timber structures involve the following moisture content ranges: ${< 12}\;\%$, $12$–${20}\;\%$, and $({>20}\;\%)$ (EN 1995-1-1 2004). Since in the wintertime $\phi$ exceeds ${85}\;\%$ for a period of several weeks, the climate can be categorized as being of class 3 ($({>20}\;\%)$). However, after a two-week averaging, the extremes change to $10$ and ${23}\;\%$, respectively. Finally, after averaging done by the FE-software concerning the time increment, the extremes change to $11$ and $20\,{\%}$, respectively (solid black line in Fig. 15). At this latter stage, the climate is classified as belonging to class 2 ($12$–${20}\,\%$). This averaging process eliminates spikes in $w_{\infty }$, which affect mainly the moisture content of the surface of the wood and thus the drying stress.

Fig. 15

Changes in moisture content $w$ at the centre and the surface of a beam A3 and b beam A2, plotted together with the original and the averaged $w_{\infty }$ values, determined by use of Simpson’s formula and moisture content profiles over the width $b$ of c beam A3, and d beam A2, after ${185\, {\text {days}}}$ of testing

Two different beam sizes are displayed in Fig. 15. As expected, beam A3 (Fig. 15a), which has a smaller cross section, reacts much quicker to changes in climate than beam A2 (see Fig. 15b, c), with its larger cross section. To quantify this, $\Delta t_{1}$ and $\Delta t_{2}$ were introduced for analysing the shift in moisture content peaks between surface and centre of the beam that occurred around March (after $120$ days) and in May (after $185$ days) (Fig. 15a, b). For beam A3, the first shift takes place during a period of $10.6$ days and results in a difference in moisture content of ${0.2}\,{\%}$. The second shift takes place during a period of $13.6$ days and leads to a difference in moisture content of ${0.77}\,{\%}$. For beam A2, the shifts $\Delta t_{1}$ and $\Delta t_{2}$ extend over a much longer time frame than for beam A3. Periods of $51.1$ and $103.2$ days are involved, and result in a difference in moisture content of $1.62$ and ${3.24}\;\%$, respectively. Comparing beam A3 in Fig. 15c with beam A2 in Fig. 15d, over the same period of time in both cases, shows there is a much quicker drying process for beam A2 than for beam A3, the latter being much closer to having a ${10}\;\%$ moisture content and a larger moisture gradient between centre and surface.

Deflection

The experimentally obtained deflection curve $u_{y,2}$ (grey) was found to be smooth, with a slight decrease during the period between the months of April (about ${150\,{\text {days}}}$ after the testing had begun) and November (about ${400\,{\text {days}}}$ after it had begun) (Fig. 16). During these months, a relatively low w ${\text{value}}$ and high $T$ value could be noted (with extremes of 22.3 °C and 10 % MC) as compared with the winter months (extremes of 4.7 °C and ${23}\quad{\%}\quad{\text{MC}}$. Figure 16 shows there to be a relatively close agreement between the simulation results and the experimental results. A discrepancy can be noted, however, between the experimental data for the first ${200\,{\text {days}}}$ of testing and the simulation data (see $\Delta u_{y,2}$ in Fig. 16), which can be attributed to the choice of creep parameter, that will be discussed in the next section, In addition, the simulated deflection curve showed clear fluctuations, caused by changes in moisture content.

Fig. 16

Validation of the numerical model on the basis of experimentally obtained deflection $u_{y,2}$, obtained with use of the four-point bending test setup exposed to natural climate conditions

Strain terms

The different terms that serve to define the strain component $\varepsilon_{l}$ as it applies to beams A2 and A3 can be found in Fig. 17. Values of the strain component for these two beams were measured in a mesh element located on the lower surface of the beams at mid-span. Note that the material points related to the elements used for beam A2 and those related to the elements used for beam A3 are not located at the same distance from the surface. The total strain $\varepsilon_{l}$ for both beams develops in a similar way, and in the same range order for both. The second elastic strain term, $\varepsilon_{e,2}$ (see Eq. (4)), is the most dominant one, being followed by the mechano-sorptive strain, $\varepsilon_{ms}$, which increases continuously over time. Hygro-expansion, $\varepsilon_{h}$ and the first elastic strain term, $\varepsilon_{e,1}$ are the two measures most strongly affected by the climate, their values during the summer also differing markedly from those during the winter. Regarding the creep parameters (see Fig. 12b), the creep strain variable $\varepsilon_{c}$ moves in the direction of becoming an asymptote almost directly after loading. Since the different parameters here were obtained at a ${7.35}\;\%$ moisture content level of the wood and are highly moisture-content dependent (Engelund and Salmen 2012), the creep strain is very likely to be underestimated.

Fig. 17

Different strain-related terms that serve to define the strain component $\varepsilon_{l}$, which are taken from a specific mesh element located on the lower surface of a beam A2 and b beam A3 at mid-span (×). Note that h is the height of the beam, and b the width of the beam

Comparing the strain terms for beam A2 with those for beam A3, one can note that the $\varepsilon_{h}$ and $\varepsilon_{ms}$ values are much higher for beam A3 than for beam A2. This is due to the fact that the interaction between the beam’s cross section and the surrounding climate is much greater in the A3 than in the A2 beam (Fig. 15 and Table 5), and in the case of $\varepsilon_{h}$ also due to a radial variation in $\alpha_{l}$, which results in a larger hygro-expansion coefficient for beam A2 in the location that is analysed. For more details, see the next section ‘Stress distribution’.

Table 5 Magnitude of different strain terms that make up $\varepsilon_{l} ,$ in the case of beam A2 and beam A3, respectively, after $556{\text{ days}}$ of testing

Stress distribution

In Fig. 18, the effects of changes in the moisture level of the wood on the stress states $\sigma_{l}$ and $\tau_{lt} ,$ being brought about by mechanical loads, are shown for beams A2 and A3 when exposed to natural climate conditions. For beam A2, the pith was located within the cross section of the beam ($y = 97.5 \;{\text{mm}}, z = 45\;{\text{mm}}$), whereas for beam A3, it was located outside the cross section ($y = 15\; {\text{mm}},{ }z = 145\;{\text{ mm}}$). The graphs in the figure displaying the $\sigma_{l}$ stress-level changes for beams A2 and A3, as obtained from a mesh element located at mid-span on the lower exchange surface, greater fluctuations for beam A3, and more seasonal variation for beam A2. The spikes seen in the stress-level changes in beam A2 show a ${17.1}\;\%$ increase as compared with the initial stress level. The differences that can be seen between beams A2 and A3 can be explained by the much quicker change in moisture content in beam A3. Both beams have a hygro-expansion coefficient $\alpha_{l}$ that changes in a radial direction. At the location analysed in Fig. 18a, this results in a much lower value of $0.0019 \left( - \right)$ for beam A3 and a much higher value for beam A2 of $0.0035 \left( - \right)$. However, when a constant and much higher value of $\alpha_{l}$ for beam A3 was initiated ($0.012 \left( - \right)$ (Bengtsson 2001)), a ${6.4}\;\%$ increase is seen for beam A2 and a $1.9\;\%$ increase for beam A3 after around ${185\,{\text{days}}}$ of testing. The results concerning the stress $\tau_{lt} { }$ variable were similar. This can be explained by the fact that the stress state here is a product of the translation of the much larger normal stress value $\sigma_{x}$.

Fig. 18

Stress changes in $\sigma_{l}$ and $\tau_{lt}$ for beams A2 and A3. Note that $h$ is the height of the beam, $b$ is the width of the beam, the stress field is projected onto half of the beam, and the deflection that occurs is magnified by a factor $10.$

The changes that occur under natural climate conditions in beams A2 and A3 lead to phenomena similar to those that occur in beam A1 (Figs. 14, 19). In the areas in which the tangential material directions align with the exchange surface, high concentrations of $\sigma_{t}$ develop. Here, steep moisture gradients located directly below the exchange surface result in shrinkage or swelling, which is counteracted by the inner section of the beam, which is still unaffected by the climate changes that have occurred (Fig. 19).

Fig. 19

Development of $\sigma_{t}$ over time of beams A2 and A3 as obtained from the mesh element located on the exchange surface or at the centre ${ }\left( \times \right)$. Note that $h$ is the height of the beam, $b$ is the width of the beam, the stress field is projected onto half of the beam, and the deflection is magnified by a factor of 10

Both beams are subjected to stress overturn. However, the smaller geometry of beam A3 results in changes occurring more quickly and in the stress values being lower than those in the larger beam, A2. The maximum tensile stress for $\sigma_{t} ,$ obtained in the month of May after about ${200\,{\text{days}}}$ of testing, is much lower for beam A3 $\left( { \approx 0.68 \;{\text{MPa}}} \right)$ than for beam A2 $\left( { \approx 1.4 \;{\text{MPa}}} \right)$. For Norway spruce, the tangential strength of the wood at different moisture content levels and at a temperature of 20 °C is about $2\; {\text{MPa}}$ (Hanhijärvi 1998). Since $\phi$ and $T$ are averaged before being inserted into the analysis and averaged too in the analysis itself with each time increment, these stress peaks are readily underestimated.

Conclusion

The three-dimensional numerical model that was employed was found to be able to simulate long-term beam bending under both controlled and natural climate conditions. The model, implemented in a powerful FE-software, appeared to be ideal for computing and visualizing complicated stress fields characterized by a combination of differing material orientations, as well as marked climate changes, long-term effects, and mechanical loads. On the basis of simulations made here for both applications, the following conclusions could be drawn:

• The small clear-wood beams analysed as part of the second application show the strong effect that spiral grain and climate have on deflection, calibrated material parameters and normative stress states.

• The results from the third application show that the larger beam exposed to mechanical load and natural climate experience a slower change in moisture, smaller moisture gradients, more seasonal fluctuation in longitudinal stress, tangential stress and longitudinal-tangential shear stress, and higher drying stress in tension compared to the smaller beam.

• Analysis of the material orientation contributes to an understanding of why bending tests are not a particularly suitable method of determining how reliable various material parameters are here. The results obtained showed spiral grain in particular to have a strong effect on the stress distribution and deflection that is characteristic of wooden beams generally and is difficult to be eliminated from slender beams.

• Climate had the strongest effect on the stress states concerning the $t$-direction of the material orientation. For each beam configuration and each type of climatic condition that was involved, stress change between tension and compression was found to occur, this leading to an alternating of the stress states involved between the surface and the inner sections of the beams. This created an environment in which both surface and inner checks could develop and could result in splits or cracks, weakening the shear capacity of the beams.

• Calibration of the numerical model benefitted from the experimental methodology that was adopted, whereas it provided a framework in which not only moisture content but also elastic, creep, hygro-expansive and mechano-sorptive behaviour could be fitted to the experimental data. Although the experiment to validate the model was sufficient, it would have benefitted from constant recordings of changes in mass.

• A major challenge in regard to three-dimensional models is that of obtaining a sufficient amount of reliable experimental data for gaining an adequate understanding of the many material parameters needed to describe the moisture flow, and the hygro-mechanical and viscoelastic behaviour of Norway spruce, or of whatever type of wood is of interest.

• The model generally employed for describing creep trends towards an asymptote, making it difficult to predict creep outside of the time frame in which the model is calibrated.

Future studies

The adjustment of the material parameters of the wood that was studied to differing climatic conditions was used here to simulate moisture flow and material stress in the wood by use of linear or exponential mathematical expressions. Most of the expressions were calibrated on the basis of different sample sets that consisted of a variety of wood species, specimen sizes, test methods and climatic conditions. To overcome its limitations, an experimental investigation could be conducted, aimed at identification of the moisture and temperature dependency of material parameters related to specific models of relevance that deal with elastic, hygro-expansive, creep and mechano-sorption behaviour. The investigation can be conducted at relevant levels of temperature and moisture content of the wood, loaded in primary material directions using small clear-wood specimens of Norway spruce. By means of matched specimens, the sample set or sets involved could be used to simulate and validate the long-term behaviour of both medium and large sized beams.