Three-dimensional orthotropic nonlinear transient moisture simulation for wood: analysis on the effect of scanning curves and nonlinearity

Abstract

This paper introduces, with the development of user-subroutines in the finite-element software Abaqus FEA^®, a new practical analysis tool to simulate transient nonlinear moisture transport in wood. The tool is used to revisit the calibration of moisture simulations prior to the simulation of mechanical behaviour in bending subjected to climate change. Often, this calibration does not receive sufficient attention, since the properties and mechanical behaviour are strongly moisture dependent. The calibration of the moisture transport simulation is made with the average volumetric mass data experimentally obtained on a paired specimen of Norway spruce (Picea abies) with the dimensions $30\times 15\times 640\, {\mathrm{mm}}^{3}$. The data, from a 90-day period, were measured under a constant temperature of 60 °C and systematic relative humidity cycles between 40 and 80%. A practical method based on analytical expressions was used to incorporate hysteresis and scanning behaviour at the boundary surface. The simulation tool makes the single-Fickian model and Neumann boundary condition readily available and the simulations more flexible to different uses. It also allows for a smoother description of inhomogeneity of material. The analysis from the calibration showed that scanning curves associated with hysteresis cannot be neglected in the simulation. The nonlinearity of the analysis indicated that a coherent set of moisture dependent diffusion and surface emission coefficient is necessary for the correct description of moisture gradients and mass transport.

Introduction

Lack of standard test method leads to diversity

A popular method to determine the strength and stiffness of wood is the three- or four-point bending test. This method is not the most straightforward, though in bending the beam experiences both tension and compression simultaneously and therefore cannot result in valid material level properties (Morlier 1994; Muszyński 2006). However, the flexural test is still a common method to experimentally calibrate (numerical) models, especially long-term behaviour when subjected to dynamic changes in climate (Honfi et al. 2014; Ma et al. 2017; Mohager and Toratti 1993). Here, the definition of calibration is the iterative process of adjusting material properties and comparing the model to experimental data until good agreement is found.

There are no standard test methods available to study the long-term behaviour of timber beams. Therefore, much diversity is found between published experimental setups (Morlier 1994), but also between the calibrations from simulations based on experimental data. Often, only the total deflection measured between the bearing supports is used and no elaboration is given of the measurements made of mass change or hygro-expansion. The flow (Skaar 1988; Yeo and Smith 2005; Yeo et al. 2002) and mechanical properties (Dinwoodie 1981; Siau 1995) as well as the mechanical behaviour associated with long-term behaviour under climate variations are strongly dependent on moisture (Ormarsson 1999; Ranta-Maunus 1975). Therefore, the importance of performing a calibration of the moisture flow model (moisture model) prior to the calibration of the moisture induced mechanical behaviour (distortion model) (Angst-Nicollier 2012) is understood. However, in the context of the four-point bending test, this is not always sufficiently exercised.

The moisture model presented in this paper is part of a more thorough study that focuses on improving the calibration of models based on long-term experimental data obtained with the traditional four-point bending setup subjected to variations in climate. The overall experiment is based on the methodology, where the deflection curve is conceptually separated into an elastic, creep, hygro-expansion and mechano-sorptive component. A climate and load schedule is used to stimulate the dissection of the curve and, together with measurements of mass change, helps to isolate the effects of mechano-sorption. The mechano-sorption is a second order phenomenon, i.e. it cannot be measured independently, such as the elastic and creep component (Muszyński 2006). The current publication will focus on the calibration of the moisture model based on experimental data. The calibration of the distortion model will be treated in a future publication.

Single-Fickian model more suitable in practice

Moisture transport in wood at the microscopic level is observed below the fibre saturation point (FSP) as water vapour diffusion in the lumen, bound water diffusion in the cell walls and as sorption processes between these two phases (Frandsen 2007; Krabbenhøft 2003). The most traditional way to simulate moisture transport is with a single-Fickian model (Fick 1855–1995), which can describe moisture transport below FSP on the scale of annual rings. A more advanced model is the multi-Fickian model (Eitelberger and Hofstetter 2011; Eriksson et al. 2007; Fortino et al. 2013a; Frandsen et al. 2007a; Johannesson 2019; Konopka and Kaliske 2018) or the multi-phase moisture transport model (Alexandersson et al. 2016; Di Blasi 1997; Fortino et al. 2013b; Janssen et al. 2007; Krabbenhøft and Damkilde 2004; Pang et al. 1995; Perré and Turner 1999; Younsi et al. 2007). The multi-Fickian model uses separate equations to describe each moisture phase and coupling terms to describe the sorption rate between phases. The classical multiphase model assumes equilibrium between moisture phases at all times. The equations are not separated, and a total diffusion tensor or matrix is used to combine the phases (Krabbenhøft and Damkilde 2004).

Hence, the multiphase models are more accurate from a physics standpoint. The single-Fickian models tend to be suitable at lower relative humidity (RH) (below 65%) (Frandsen 2007; Konopka and Kaliske 2018; Krabbenhøft 2003). Here, the moisture transport is characterised by slow bound water transport and relatively fast sorption, and the ratio between surface emission and diffusion coefficient used to describe the transport is high (Konopka and Kaliske 2018). Compared to the multiphase model, the single-Fickian model is simpler to implement and a more efficient model to simulate engineering problems.

The importance of modelling sorption hysteresis

For wood, the sorption isotherm plays a fundamental role in the numerical simulations of moisture transport, both for the single-Fickian and multiphase models. The sorption isotherm, i.e. the relation between equilibrium moisture content (EMC) (0–FSP%) and ambient RH (0–100%), is characterized by an envelope created from the so-called border desorption (drying) and adsorption (wetting) isotherms, which are obtained at the first complete drying and wetting (Frandsen 2007; Time 1998). Traditionally, the sorption isotherm for wood is determined on wood fibres (Shi and Avramidis 2017a, c), but the effect is also present on large-scale wooden and timber elements (Svensson et al. 2011).

The difference in EMC between the desorption and adsorption isotherms for equal ambient climatic conditions is referred to as sorption hysteresis (Engelund et al. 2013; Fredriksson and Engelund 2018). This indicates that the EMC is not only dependent on RH, but also on moisture history. Between the border isotherms, the EMC follows a unique path that is dependent on the previous moisture content (MC) changes experienced, which creates so-called intermediate or scanning curves (Frandsen 2007). The paths occur under a slope that is not necessarily identical for desorption and adsorption curves. It is not uncommon to use the average sorption isotherm, without the effect of hysteresis, to describe the relation between EMC and RH when performing simulations (Salin 2011). However, such an assumption can lead to a significant deviation from the measured MC (Frandsen et al. 2007b; Time 2002b) and will affect the moisture gradients and rate of moisture change, important to describe the moisture dependent mechanical behaviour.

The most traditional way to model hysteresis and its associated scanning curves is with the independent domain model (Everett and Whitton 1952). Peralta (1995) first applied this model to wood. The challenges with this model are the lack of experimental data on scanning curves and the double integral that needs to be solved to obtain EMC (Frandsen 2007; Shi and Avramidis 2017b). A more recent hysteresis model was proposed by Frandsen et al. (2007b), which is largely used in the hygroscopic computational literature for wood, and other building materials (Fortino et al. 2019; Huč 2019; Johannesson and Janz 2009; Svensson et al. 2011). The model can be used in combination with the multi-Fickian model, where both RH and MC are driving forces behind moisture transport below FSP.

Novelty is of technical nature

The purpose of the current work is to present a new practical tool to simulate nonlinear transient moisture transport. The user-subroutine frameworks associated with this tool are provided by finite-element software Abaqus FEA^® (Dassault Systemes, Vélizy-Villacoublay, France). The work further develops the research previously published by Florisson et al. (2019). The research objectives related to the development of the simulation tool are: (1) the implementation of the single-Fickian approach and the moisture-dependent diffusion coefficient in user-subroutine UMATHT, which is usually used in the field of polymers, freezing soil and fire engineering; (2) the implementation of the moisture flux vector normal to the boundary surface in user-subroutine FILM, which is described by the moisture dependent material property named the surface emission coefficient (Yeo and Smith 2005; Yeo et al. 2002); (3) implementation of a practical method in user-subroutine FILM to describe sorption hysteresis and scanning behaviour at the boundary surface. Unique expressions are used for each adsorption and desorption curve to interpolate between RH and the experimentally obtained EMC at the end of each RH-period.

The simulation tool will be used to revisit the calibration of moisture transport model made prior to the simulation of moisture-dependent mechanical behaviour in bending. The calibration is based on experimentally obtained volumetric average moisture content. The experimental method is relatively simple compared to other methods, such as nuclear magnetic resonance (NMR) and computerised tomography (CT), which give a more detailed representation of the moisture gradient and change in wooden components (Rosenkilde 2002; Wiberg 1998). However, the test method is chosen to allow for the simultaneous measurement of mass and deflection in the same climatic environment for a long period of time. The research objectives related to the analysis of the calibration are: (1) to use the simulation tool to perform a fit on moisture content profiles obtained under similar climatic conditions and thus obtain a coherent set of moisture dependent diffusion and surface emission coefficients prior to the calibration of the moisture transport model, and (2) to visualise the effect of moisture gradients and mass transport by analysing sorption behaviour and nonlinearity of analysis.

Material and methods

The following section will introduce the governing equation for the nonlinear transient moisture transport together with the constitutive laws. A finite element formulation for the problem will be presented and the Jacobian matrix will be derived to iteratively solve the nonlinear equation system. The implementation of the equation system into Abaqus, using the user-subroutines, will be highlighted as well. Finally, the simulation model, the experiment used to calibrate this model, and the method used to evaluate the calibration will be discussed.

Moisture transport model

A single-Fickian approach is used to analyse the orthotropic nonlinear transient moisture flow in wood. The governing equation is derived from the mass balance law and expressed as

$$c \rho \frac{\partial w}{\partial t}+\varvec{\nabla} \boldsymbol{\cdot} \left({\varvec{q}}\right)=0$$

(1)

where $c$ is the moisture capacity in ${\text{m}}^{3}/{\text{kg}}$, $\rho$ is the density of liquid moisture in ${\text{kg}}/{\text{m}}^{3}$, $w$ is the unit-less moisture content, $t$ is time in $\mathrm{h}$, $\varvec{\nabla}$ is the divergence and ${\varvec{q}}$ is the moisture flux in ${\text{m/h}}$. The $(\boldsymbol{\cdot} )$ denotes the operation where the components of $\varvec{\nabla}$ are applied to the corresponding components of ${\varvec{q}}$ and summated. The constitutive law for the moisture flux is described by Fick’s law of diffusion defined as

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

(2)

where ${\varvec{D}}$ is a moisture and temperature dependent diffusion matrix in ${\text{m}}^{2}/{\text{h}}$ and $T$ is temperature in $^\circ{\mathrm{C}}$. The parameter with the over-line, i.e. $\left(\bar{\boldsymbol{\cdot}}\right)$ in Eq. (2) refers to the orthogonal local coordinate system. The constitutive law for the flux normal to the boundary ${q}_{n}$ in ${\text{m/h}}$ is defined as

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

(3)

where $s$ is a moisture and temperature dependent surface emission coefficient in ${\text{m/h}}$ and ${w}_{\infty }$ the unit-less moisture content in the ambient air. The geometry is assumed to be discretized into finite elements, each with an approximation of the moisture content, weight function, and their respective gradients as

$$\begin{aligned}w&={\varvec{N}}{\varvec{a}}, \quad v={{\varvec{c}}}^{{\mathrm{T}}}{{\varvec{N}}}^{{\mathrm{T}}}, \\ \quad \varvec{\nabla}w&={\varvec{B}}{\varvec{a}}, \quad \varvec{\nabla} v={{\varvec{c}}}^{{\mathrm{T}}}{{\varvec{B}}}^{{\mathrm{T}}}\end{aligned}$$

(4)

where ${\varvec{N}}$ is a shape function vector, ${\varvec{a}}$ is a nodal moisture content vector, ${\varvec{c}}$ is an arbitrary vector and ${\varvec{B}}=\varvec{\nabla} {\varvec{N}}$. A backward Euler time integration scheme is used to account for the derivative in time and by following the steps suggested by Johannesson (2019) leads to the definition of the residual ${\varvec{\psi}}$ for each finite element as

$$\begin{aligned}{\varvec{\psi}}\left({{\varvec{a}}}_{n+1}\right)&=:\frac{1}{\Delta t}{\int }_{V}{{\varvec{N}}}^{{\mathrm{T}}}c\rho {\varvec{N}}\left({{\varvec{a}}}_{n+1}-{{\varvec{a}}}_{n}\right){\mathrm{d}}V\\& \quad +{\int }_{V}{{\varvec{B}}}^{{\mathrm{T}}}{\varvec{D}}{\varvec{B}}{{\varvec{a}}}_{n+1}{\mathrm{d}}V\\& \quad +{\int }_{A}{{\varvec{N}}}^{{\mathrm{T}}}s{\varvec{N}}{\varvec{a}}{\mathrm{d}}A-{\int }_{A}{{\varvec{N}}}^{{\mathrm{T}}}s{w}_{\infty }{\mathrm{d}}A\end{aligned}$$

(5)

where the subscript to a quantity, ${(\boldsymbol{\cdot})}_{n+1}$ denotes the value of the quantity at the current time step, and ${(\boldsymbol{\cdot})}_{n}$ the value at previous known time step at equilibrium, separated by a time increment $\Delta t$. The last two terms in Eq. (5) will only be nonzero for the elements at the boundary of the domain.

Numerical solution

The Newton–Raphson method is used to iteratively solve the nonlinear system of equations, see Eq. (5). Assuming an initially known iteration state $i-1$, with known residual ${{\varvec{\psi}}}^{i-1}={\varvec{\psi}}({{\varvec{a}}}_{n+1}^{i-1})$ and known ${{\varvec{a}}}_{n+1}^{i-1}$ and ${{\varvec{a}}}_{n}$, a truncated Taylor expansion is performed to achieve the residual ${{\varvec{\psi}}}^{i}={\varvec{\psi}}({{\varvec{a}}}_{n+1}^{i})$ at the next iteration $i$ as

$${{\varvec{\psi}}}^{i}={{\varvec{\psi}}}^{i-1}+\left(\frac{\partial {{\varvec{\psi}}}^{i-1}}{\partial {{\varvec{a}}}_{n+1}^{i-1}}\right)\left({{\varvec{a}}}_{n+1}^{i}-{{\varvec{a}}}_{n+1}^{i-1}\right)$$

(6)

For the iteration, $i$ to be at equilibrium, i.e. ${{\varvec{\psi}}}^{i}=0,$ Eq. (6) can be rewritten into a system of linearized equations used to determine ${{\varvec{a}}}_{n+1}^{i}$. More on how to derive and solve this linearized equation can be found in Johannesson (2019). The Jacobian (or tangent diffusion matrix), ${{\varvec{K}}}_{{\mathrm{t}}}$ is defined as

$${{\varvec{K}}}_{{\mathrm{t}}}=\left(\frac{\partial {{\varvec{\psi}}}^{i-1}}{\partial {{\varvec{a}}}_{n+1}^{i-1}}\right)=\frac{{{\varvec{C}}}^{i-1}}{\Delta t}+{\tilde{{\varvec{C}}}}^{i-1}+{{\varvec{K}}}^{i-1}+{\tilde{{\varvec{K}}}}^{i-1}+{{\varvec{K}}}_{{\mathrm{c}}}^{i-1}+{\tilde{{\varvec{K}}}}_{{\mathrm{c}}}^{i-1}-{{\varvec{F}}}_{{\mathrm{b}}}$$

(7)

where using Eq. (5), we have

$${\varvec{C}}={\int }_{V}{{\varvec{N}}}^{{\mathrm{T}}}c\rho {\varvec{N}}\,{\mathrm{d}}V$$

(8)

$$\stackrel{\sim }{{\varvec{C}}}={\int }_{V}{{\varvec{N}}}^{{\mathrm{T}}}\frac{\partial c}{\partial w}{\varvec{N}}\left({{\varvec{a}}}_{n+1}^{i-1}-{{\varvec{a}}}_{n}\right){\varvec{N}}{\mathrm{d}}V$$

(9)

$${\varvec{K}}={\int }_{V}{{\varvec{B}}}^{{\mathrm{T}}}{\varvec{D}}{\varvec{B}}{\mathrm{d}}V$$

(10)

$$\stackrel{\sim }{{\varvec{K}}}={\int }_{V}{{\varvec{B}}}^{{\mathrm{T}}}\frac{\partial {\varvec{D}}}{\partial w}{\varvec{B}}{{\varvec{a}}}_{n+1}^{i-1}{\varvec{N}}{\mathrm{d}}V$$

(11)

$${{\varvec{K}}}_{{\mathrm{c}}}={\int }_{A}{{\varvec{N}}}^{{\mathrm{T}}}s{\varvec{N}}{\mathrm{d}}A$$

(12)

$${\tilde{{\varvec{K}}}}_{{\mathrm{c}}}={\int }_{A}{{\varvec{N}}}^{{\mathrm{T}}}\frac{\partial s}{\partial w}{\varvec{N}}{{\varvec{a}}}^{i-1}{\varvec{N}}{\mathrm{d}}A$$

(13)

$${\tilde{{\varvec{F}}}}_{{\mathrm{b}}}={\int }_{A}{{\varvec{N}}}^{{\mathrm{T}}}\frac{\partial s}{\partial w}{w}_{\infty }{\varvec{N}}{\mathrm{d}}A$$

(14)

The Jacobian is updated in each Newton–Raphson iteration.

Implementation user-subroutines

Abaqus uses the Newton Raphson method as a default numerical technique to solve the nonlinear equilibrium equations (Dassault Systèmes 2017a). The motivation for this choice is primarily the convergence rate obtained with this method. Abaqus FAE^® enters user-subroutines UMATHT and FILM at each iteration (Dassault Systèmes 2017b). The routines are used by Abaqus FAE^® to derive the Jacobian and define the necessary matrices and load vectors. Results are available at the beginning and end of each time increment. The subroutines are written in programming language FORTRAN (The Fortran Company, Colorado, USA). Routine UMATHT requires the user to define the following items:

• The moisture capacity $c$

• The derivative of the capacity towards moisture $\partial c/\partial w$

• The derivative of the capacity towards the spatial gradients of moisture $\partial c/\partial \bar{\varvec{\nabla} }w$

• The flux vector $\bar{{\varvec{q}}}$

• The derivative of the flux vector towards moisture $\partial \bar{{\varvec{q}}}/\partial w=\partial \bar{{\varvec{D}}}/\partial w\bar{\varvec{\nabla} }w$

• The derivative of the flux vector towards the spatial gradients of moisture $\partial \bar{{\varvec{q}}}/\partial \bar{\varvec{\nabla} }w=\bar{{\varvec{D}}}$

The routine also provides so-called state variables that can be used to store and reintroduce certain variables. Because the flux vector and moisture content are automatically reintroduced by the routines, and all other variables only need to be defined at the end of the time increment, state variables were not used in this analysis.

By assigning the subroutine FILM to the necessary surface, the routine can find the corresponding elements that need to be considered for the analysis. Important variables sent into the routine are the current moisture content values and corresponding time. The routine requires the user to define the following items:

• The ambient MC of the air ${w}_{\infty }$

• The surface emission coefficient $s$

• The derivative of the surface emission coefficient towards moisture $\partial s/\partial w$

For simplicity, it is assumed that Eq. (1) is divided by the density and capacity and that this inverse of density and capacity is implicitly taken care of by the diffusion coefficient later obtained from calibration. Further, the analysis is performed under constant temperature, simplifying the diffusion matrix. The orthotropic material orientation was defined in user-subroutine ORIENT as a local coordinate system (longitudinal, radial and tangential). The local system was described by direction cosiness, based on the global coordinate system $(x, \, y, \, \text{and} \, z)$. A more detailed description of theory and implementation can be found in Ormarsson (1999) and Ormarsson et al. (1998).

Abaqus model

A three-dimensional model was created to simulate the nonlinear transient moisture transport in a wood beam. The beam was obtained from the cross section of a log at a distance of $130\,{\text{mm}}$ from pith, see Fig. 1. The cut followed the conical slope of the tree. The wood beam had a constant width of $15\,{\text{mm}}$, a constant height of $30\,{\text{mm}}$ and a total length of $640\,{\text{mm}}$. The beam had a spiral grain angle $\theta$ of 4$^\circ$.

Fig. 1

Origin of the beam used for the moisture simulation

The moisture transport analysis was performed with a fixed time increment $\Delta t$, see “Numerical solution” section, of $0.1\,{\text{h}}$ over a total period of $2880\,{\text{h}}$$(120\,{\text{days}})$. The increment size, chosen based on the distortion simulation subsequent to the moisture transport simulation, was small enough to assure satisfactory results and large enough to get acceptable computation time. The simulation was performed with element type DC3D8 and an element size of $2.5\,{\text{mm}}$ over the cross section of the beam and $5\,{\text{mm}}$ over the length of the beam, see Fig. 2. The orthotropic material orientation of the beam as defined in user-subroutine ORIENT is also shown in Fig. 2.

Fig. 2

Local material orientation radial–tangential plane (front side beam including annual ring curvature) and longitudinal–tangential plane (right side beam including spiral grain angle $\theta$ of $4^\circ$) presented together with dimensions and mesh size

Figure 3 illustrates the temperature and RH schedules used to define the MC of the ambient air. A relatively high temperature and differential between subsequent RH were chosen to stimulate the development of mechano-sorption. This meant a constant temperature of 60 °C was used and systematic RH-cycles were initiated, covering 10 day-periods and consisting of two RH-phases; one at $40\%$ and one at $80\%$. Each phase resulted in EMC, which was used as the ambient moisture content ${w}_{\infty }$ as initiated by Eq. (3).

Fig. 3

Schedules used for temperature $\left(T\right)$, relative humidity $\left(\phi \right)$ and ambient moisture content $\left({w}_{\infty }\right)$, where $0.27$ is the fibre saturation point $\left({w}_{{\mathrm{f}}}\right)$

The RH-schedule (see Fig. 3) was defined in user-subroutine FILM. Each new phase was introduced under a slight slope that covered six hours (60 increments) to ensure numerical stability. The relation between RH $(\phi )$ and the experimentally obtained EMC $({w}_{{\mathrm{eq}}})$ was established with scanning curves. An approximation of these curves was made with Eq. (15) and (16) (Johannesson 2019). Equation (15) describes ${w}_{\infty }$ when in adsorption, and Eq. (16) when in desorption. The constants $a$ and $b$ are used to describe the shape of the curves, and ${w}_{f}$ is the moisture content at fibre saturation point as defined in Fig. 3. To increase the flexibility with which parametric studies and post-processing of results were done, scripts written in programming language Python^® (PSF, Delaware, United States) were used.

$${w}_{\infty }^{a}=\frac{a\boldsymbol{\cdot} {\phi }^{0.4}+b\boldsymbol{\cdot} {\phi }^{5}}{5}\boldsymbol{\cdot} {w}_{f}$$

(15)

$${w}_{\infty }^{d}=\left(\frac{a\boldsymbol{\cdot} {\phi }^{0.4}+b\boldsymbol{\cdot} {\phi }^{5}}{6}+0.08\boldsymbol{\cdot} \sin(\pi \boldsymbol{\cdot} \phi )\right)\boldsymbol{\cdot} {w}_{f}$$

(16)

Calibration of numerical model based on experiments

Because the hygro mechanical behaviour is strongly dependent on moisture gradient and mass transport, calibration of the numerical model used to perform the moisture simulation is essential. In this section, the experiment is presented that was needed for this calibration. The timber log of Norway spruce used in the experiment came from a private stand named Engaholms skogar, located near Lyngsåsa in southern Sweden. The specimens were obtained with the procedure previously illustrated in Fig. 1. Figure 4 illustrates the method to identify the beams. The first object in the identifier indicated the tree, the second was the log number and the third pointed towards the primary wind direction. Additionally, a fourth object was used to show whether the beam was tested in the long-term (for calibration) experiment or the static experiment (future publication). The beams used for the moisture simulation were identified as T1–2-WW-s-LT.

Fig. 4

Adopted method for identification of beams

The specimens used for the long-term experiment were all provided with a silicone sealant to prevent moisture transport in the lengthwise direction. The limited space inside the climate chamber allowed for one specimen to be monitored for mass change. This was justified by the fact that between the specimens in the sample set of four, only a $10\%$ difference in density was found, and therefore, a similar thickness of cell wall and lumen size can be assumed. Further, the specimens had similar MC, and a similar amount of latewood on visual inspection, and therefore, a comparable moisture transport.

The temperature and RH-schedules used in the experiment are displayed in Fig. 3. Each RH-phase resulted in a mass equilibrium state. This state was reached when the change in mass was below $1\%$ between consecutive days. The recorded average volumetric mass change was converted to average volumetric MC $({w}_{{\mathrm{a}}}^{{\mathrm{exp}}})$ by means of the oven-dry method (EN 13183-1 2003), see Eq. (17). Here, ${m}_{w}$ indicates the wet mass and ${m}_{d}$ the oven dry mass of the beam. The oven dry mass was reached when a change in mass of less than $1\%$ was determined between two consecutive hours. The experimental results were compared with the simulated volumetric average MC ${w}_{{\mathrm{a}}}^{{\mathrm{sim}}}$, see Eq. (18),

$${w}_{{\mathrm{a}}}^{{\mathrm{exp}}}=\left(\frac{{m}_{w}-{m}_{d}}{{m}_{d}}\right)\times 100$$

(17)

$${w}_{{\mathrm{a}}}^{{\mathrm{sim}}}=\left(\frac{1}{V}\sum_{i=1}^{m}{w}_{i}{V}_{i}\right)\times 100$$

(18)

where $V$ is the total volume, and ${w}_{i}$ is the MC and ${V}_{i}$ is the effected volume for each nodal point $i$ corresponding to the finite element mesh. The calibration of the moisture transport model based on the experimentally obtained average volumetric MC resulted in a moisture-dependent diffusion and surface emission coefficient. The expressions that were used to describe the coefficients mathematically are presented by Eq. (19) and (20), where $D$ indicates the diffusion coefficient, $s$ the surface emission coefficient, $w$ the MC and indexes $a$ to $c$ reference the constant parameters needed to describe each coefficient. The constants were established by calibration. Here, index $k$ refers to the orthotropic material directions $l$, $r$ or $t$.

$${D}_{k}={D}_{{\mathrm{a}}}\boldsymbol{\cdot} {w}^{{D}_{{\mathrm{b}}}}+{D}_{{\mathrm{c}}}$$

(19)

$$s={s}_{{\mathrm{a}}}\boldsymbol{\cdot} {w}^{{s}_{{\mathrm{b}}}}+{s}_{{\mathrm{c}}}.$$

(20)

Accuracy of calibration

The calibration was evaluated by determining the relative difference (RD) between experimental data and simulation data, see Eq. (21). The RD is based on the absolute difference between the simulation ${y}_{m}^{{\mathrm{sim}}}$ and the experimental ${y}_{m}^{{\mathrm{exp}}}$ value of each data point $m$. The most optimal value for the RD is zero per cent. Such value indicates accurate calibration, though it does not indicate possible experimental measuring errors or simulation errors due to geometric simplification, discretisation, numerical imperfections, or the fact that the finite element method is based on approximations.

$${\text{RD}}= \frac{{\text{abs}}\left({y}_{m}^{{\mathrm{sim}}}-{y}_{m}^{{\mathrm{exp}}}\right)}{{y}_{m}^{{\mathrm{exp}}}}\times 100$$

(21)

$${\text{RD}}= \frac{{\text{abs}}\left({s}_{m}^{x}-{s}_{m}^{1}\right)}{{s}_{m}^{1}}\times 100$$

(22)

The method of RD was also used to analyse the outcome of the parametric studies presented in the results and discussion chapter, see Eq. (22). Here, the reference value is ${s}_{m}^{1}$, indicating the data corresponding to the calibrated numerical model. The variable ${s}_{m}^{x}$ refers to specific data obtained from the parametric study and $x$ is a number between 2 and 8, referring to a specific simulation.

Results and discussion

In the first “Moisture content profiles” section, the numerical model will be used to simulate MC-profiles associated with flow in radial, tangential and longitudinal directions separately. The profiles come from Rosenkilde (2002). To validate the routines, the moisture-dependent diffusion and surface emission coefficient obtained from this analysis are compared to the coefficients obtained in a previous publication by Florisson et al. (2019). In the second “Hysteresis” section, the experimentally observed hysteretic behaviour, together with its implementation into the routine that covers the boundary condition, is discussed. The flow properties found in the first section are used as a starting point for the calibration presented in the third “Calibration of numerical model based on experiment” section. The calibrated model is used to analyse the effect of scanning curves and nonlinearity in the final “Moisture change and gradients” and “Nonlinearity analysis” sections of this chapter. These items have significant effect on moisture flow and gradients, which are fundamental to describe the hygro mechanical behaviour of wood.

Moisture content profiles

To assess numerical models for moisture transport, it is recommended to validate the models with moisture profiles. Rosenkilde (2002) provides MC-profiles in radial, tangential and longitudinal directions for Scots pine dried from green state to EMC at $60\,^\circ{\mathrm{C}}$ and $59\%$ RH. The profiles were obtained on specimens with dimensions $29\times 29\times 29\,{\text{mm}}^3$. The results in radial direction were obtained with a combination of the so-called slicing technique and oven-dry method, and in longitudinal direction with CT-scanning. More about these techniques and the experiment can be read in Rosenkilde (2002). The fits of the MC-profiles in the radial, tangential and longitudinal directions below FSP using the numerical moisture transport model are presented in Fig. 5.

Fig. 5

Comparison of simulated (solid line) and experimentally obtained moisture content profiles (scattered data) below fibre saturation point in a radial/tangential, b longitudinal direction (Rosenkilde 2002)

The distribution of the moisture-dependent flow parameters obtained, with the fit on the data by Rosenkilde (2002), is presented in Fig. 6. The diffusion coefficient is the same in both the radial and tangential directions. Figure 6a presents these coefficients together with the values obtained in Florisson et al. (2019). Great similarities are found between the simulation made with the standard analysis (without routines) provided by Abaqus FAE^® in combination with tabulated flow properties, and the simulation made using the user-subroutines in combination with the expressions provided by Eq. (19) and (20). Figure 6b presents the flow properties associated with the longitudinal direction. The diffusion coefficient is constant with respect to moisture change, which was also seen in Rosenkilde (2002). The distribution of the surface emission coefficient is similar to that of the radial/tangential directions, which have an RD of maximum $21\%$ and a minimum of $19\%$ in the MC-domain between $10$ and $35\%$. The variables that describe the flow properties by using Eq. (19) and (20) are collected in Table 1.

Fig. 6

Moisture-dependent diffusion and surface emission coefficient obtained in fit (solid line) on moisture content profiles covering a radial and tangential direction, together with data (scattered) obtained in Florisson et al. (2019), b longitudinal direction

Table 1 Variables needed to describe the diffusion coefficient $\left(D\right)$ in the orthotropic directions $\left(l,r,t\right)$ and corresponding surface emission coefficient $\left(s\right)$ used to describe the flux vector normal to the boundary. Values are obtained with fit made of moisture content profiles by Rosenkilde (2002) in the MC-range $8$–$30\%$ at 60 °C $\left(a\right)$ and $\left(b\right)$ relate to Fig. 6

Hysteresis

The analysis of hysteresis and the scanning behaviour experienced due to the cyclic movement between $40$ and $80\%$ RH are illustrated in Fig. 7. Figure 7a presents the experimentally obtained values for EMC at $40$ and $80\%$ RH, together with a plot of a sorption isotherm of Norway spruce obtained at a temperature of 20 °C by Bratasz et al. (2012). The sorption isotherm was obtained with the Guggenheim, Anderson and de Boer (GAB) method, presented by Eq. (23). Here, ${V}_{m}$ is the monolayer capacity, ${c}_{{\mathrm{e}}}$ is the energy constant related to the difference of free enthalpy of water molecules in the upper sorption layers and monolayers, and $k$ is the free enthalpy of water molecules in the pure liquid and upper sorption layers. The method assumes water molecules to be present in a monolayer or upper layer in the cell wall, as well as capillary condensation in the lumen, dependent on the RH. The water molecules in the monolayer have direct contact with the surface layer, whereas those in the upper layer do not. In this context, enthalpy is the thermodynamic quantity, which is equivalent to the total heat content of a system.

Fig. 7

Interpretation of sorption behaviour at a constant temperature of 60 °C used in calibration of numerical model based on long-term experiment, a plot of sorption isotherm of spruce at 20 °C (solid line) (Bratasz et al. 2012) including an adjustment of this plot to enclose the data points at $40$ and $80\%$ RH obtained in long-term experiment (dotted line), b interpolation data points between $40$ and $80\%$ RH at first adsorption, c interpolation data points between $40$ and $80\%$ RH at first desorption, d approximation of scanning curves for all data points in long-term experiment

$${w}_{{\mathrm{eq}}}=\frac{{V}_{m}\boldsymbol{\cdot} {c}_{{\mathrm{e}}}\boldsymbol{\cdot} k\boldsymbol{\cdot} \phi }{(1-k\boldsymbol{\cdot} \phi )(1+({c}_{{\mathrm{e}}}-1)\boldsymbol{\cdot} k\boldsymbol{\cdot} \phi )}$$

(23)

Cycling the wood between two RH-values will result in closed hysteresis loops (scanning curves) bound by the adsorption and desorption isotherm. To get an indication of the two boundary curves at 60 °C, the variables used to describe the sorption isotherm at 20 °C were adjusted to enclose the experimentally obtained datapoints. The results show that the sorption isotherm lowers for higher temperatures (Engelund et al. 2013; Shi and Avramidis 2017a; Skaar 1988), as seen when the average EMC of $12.0$ and $7.7\%$ ($80$ and $40\%$ RH) obtained in the present study is compared with the average of $16.4$ and $8.5\%$ MC found for the sorption isotherm by Bratasz et al. (2012) at the same RH-levels. The variables needed to describe both sorption isotherms are collected in Table 2.

Table 2 Three parameters needed to describe the sorption isotherm using the GAB method presented in Eq. (23) and illustrated in Fig. 7

The scanning curves generally do not develop in a straight manner (Salin 2011; Shi and Avramidis 2017a; Time 2002a), but in a similar trend to the boundary curves. With Eq. (15) and (16), an approximation is made of what these scanning curves could look like. In Fig. 7b, c, the interpolation between the experimentally obtained data points at first desorption and adsorption are presented. In Fig. 7d, all scanning curves are collected. It can be seen that the hysteretic behaviour results in larger differences in EMC at $80 \%$ RH $(\Delta 1.1 \%)$ than at $40 \%$ RH $(\Delta 0.4 \%)$. It is difficult to compare the experimental results with other publications (Fredriksson and Engelund 2018; Shi and Avramidis 2017a), where the scanning curves are limited and often start at adsorption. However, it is seen in Fredriksson and Engelund (2018) that the distance between desorption and adsorption boundary curve is much smaller at $40\%$ RH than at $80\%$ RH, also suggesting a similar scanning range. Additionally, this difference is seen in Johannesson and Janz (2009), as an example of the hysteretic behaviour simulated with the hysteresis model by Frandsen et al. (2007b) starting from the adsorption curve and scanning between $74$ and $88\%$ RH. All variables needed to make the scanning curves illustrated in Fig. 7 are collected in Table 3.

Table 3 Variables used to describe the scanning curves at adsorption $(\mathrm{ad})$ and desorption $(\mathrm{de})$ with a fibre saturation point of $0.27\, ({-})$

Calibration of numerical model based on experiment

Figure 8 presents the results of the calibration made from the moisture transport model based on the average MC. The calibration started with the flow properties obtained in “Moisture content profiles” section, see Table 1. Figure 8a presents these results, together with the experimental data. This figure also shows the RD between simulated and experimental data, which was determined based on Eq. (21). The next step was to adjust the flow properties until a better RD distribution was found. Figure 8b illustrates the final fit made in the calibration process. The flow properties related to this fit are collected in Table 4. To emphasize the EMC after each RH-phase and the effect hysteresis has on the EMC, Fig. 9 is presented.

Fig. 8

Average moisture content $({w}_{{\mathrm{a}}})$ and relative difference $({\text{RD}})$ between experiment $(\mathrm{e})$ and simulation $(\mathrm{s})$ data for a simulation made with flow properties obtained with fit of moisture content profiles presented by (Rosenkilde 2002), b calibration made of moisture model based on long-term experiment

Table 4 Variables needed to describe the diffusion coefficient $\left(D\right)$ in orthotropic directions $\left(l,r,t\right)$ and surface emission coefficient $\left(s\right)$ used to describe the flux vector normal to the boundary. Values are obtained with calibration made of numerical model based on experimental data between $7{-}13\%$ MC and $60\,^\circ{\mathrm{C}}$

Fig. 9

Simulated average volumetric moisture content $({w}_{{\mathrm{a}}})$ including colour plots emphasising the effect of hysteresis on adsorption and desorption side for specific points in time (${t}_{1}$ to ${t}_{13})$

The obtained RD-datasets were statistically analysed with a built in function of MATLAB^®. The analysis showed that both sets were not normally distributed. The RD-dataset presented in Fig. 8a had a mean value of $2.7\%$. By adjusting the diffusion and surface emission coefficient, a better fit between simulation and experiment was obtained in Fig. 8b. The mean value of this dataset is $2.3\%$.

Although the mean RD-values did not reduce significantly, a visual inspection of the graphs in Fig. 8a, b shows considerable differences. In the RD-plot, three locations are highlighted. The peaks, indicated by the first point, are largely caused by the subtle difference in an initial rapid increase or decrease in MC after a change in RH. These points have a significant influence on the mean RD, but are difficult to eliminate. However, the RH is initiated under a time span in the simulation to avoid numerical instability issues. The second and third points highlight the beginning of a desorption and adsorption phase. Despite the higher peaks seen in adsorption, the general fit results in similar values for RD. Figure 9 shows that with hysteresis included in the numerical model, the EMC has considerably different values after every RH-phase.

To finalise the discussion about the calibration of the model, the flow properties of Table 1 are compared with those collected in Table 4. For the diffusion coefficient corresponding to the radial/tangential direction, a maximum RD is found of $82\%$ around $10\%$ MC and a minimum RD of $18\%$ around $35\%$ MC. The RD for the surface emission coefficient finds a maximum of $15.7\%$ at $35\%$ MC and a minimum of $0\%$ at $15\%$ MC. The diffusion coefficient in the longitudinal direction was adopted from Table 1. However, despite the spiral grain, no large reduction in the RD presented in Fig. 8 could be established by changing this property in value.

The radial/tangential diffusion coefficient and surface emission coefficient presented in Table 4 are shown in Fig. 10. The variation of these properties with respect to MC at various times and at different locations across the beam-section is shown. Figure 10a, c illustrate, in particular, the opposing behaviour of diffusion and surface emission coefficient to moisture, which has been observed by many researchers, such as Siau and Avramidis (1996), Yeo and Smith (2005) and Yeo et al. (2002). Especially interesting, shown by Fig. 10b, d, is the effect of hysteresis on the flow properties, and, in Fig. 10e, the variation of the diffusion coefficient over the cross section with respect to time seen for the first two RH-phases.

Fig. 10

Calibrated flow properties, a moisture-dependent diffusion coefficient in radial/tangential direction $(D)$, b variation of the same $D$ over time for a nodal point at the surface and the centre of the cross section, c moisture dependent surface emission coefficient $(s)$, d variation of $s$ with respect to time for a nodal point at the surface, e colour plots of $D$ corresponding to the points illustrated in b (colour figure online)

Moisture change and gradients

The average sorption isotherm is often used to simulate the moisture transport needed to analyse the moisture dependent mechanical behaviour. These choices influence the rate with which moisture changes, and the steepness of the gradients are experienced inside the material, which can have a significant effect on strain development in wood. In this section, the simulation from “Calibration of numerical model based on experiment” section (simulation 1) will be compared with a simulation where hysteresis at the boundary is neglected (simulation 2). The latter simulation uses the final values found in the calibration experiment, which is a constant EMC of $11.65\%$ at $80\%$ RH and an EMC of $7.6\%$ at $40\%$ RH as boundary condition.

In Fig. 11, the development of the average volumetric MC with respect to time is presented for both simulations. Eight points in time are selected for further study, which are strategically chosen $25\,{\text{h}}$ after an increase in RH, when a moisture gradient is clearly present. The points were chosen in adsorption, when the hysteretic effects on the sorption isotherm are most significant. Colour plots were provided to visualise the MC-distribution over the cross section. The last plot in the sequence of colour plots shows the point where simulations 1 and 2 coincide.

Fig. 11

Simulated average volumetric moisture content ${(w}_{{\mathrm{a}}})$ with respect to time $(t)$ for a simulation with (solid line) and without (dotted line) the scanning curves included at the boundary, together with colour plots that show the cross sectional variation of moisture at eight specific points in time $({p}_{1}$ to ${p}_{8}$) for the simulation with and one colour plot for the simulation without $({p}_{9}$) sorption effects(colour figure online)

For the eight points illustrated in Fig. 11, path plots that run from the upper exchange surface to the lower exchange surface for both simulations 1 and 2 were created, see Fig. 12. As a reference, the volumetric average MC belonging to the eight points are also plotted. The MC during adsorption shows to be systematically higher for simulation 1 than for simulation 2.

Fig. 12

Moisture content plots $({w}_{{\mathrm{p}}})$ over the height of the cross section at unique nodes for eight specific points in time previously illustrated in Fig. 11 for the simulation where hysteresis was and was not included as a boundary condition

To increase the understanding of the path plots presented in Fig. 12, values of the average MC, the MC at the centre of the path plot and the moisture gradient at a $2.5\,{\text{mm}}$ distance from the upper and lower exchange surface are collected in Table 5. Between the two simulations, a significant difference is found in the average MC and MC at the centre of the path plot, where both find their maximum at around $10\%\, {\text{RD}}$. The moisture gradient is unaffected by spiral grain and is higher for simulation 2 than for simulation 1, with a maximum RD of $10\%$ after $385\,{\text{h}}$ of testing. For simulation 1, the gradient gradually increases with increasing test time.

Table 5 Data from the eight points in time treated in Figs. 11 and 12 for a simulation with $(h)$ and without $(nh)$ the effect of hysteresis included at the boundary

Nonlinearity analysis

When simulating moisture transport, it is generally known to incorporate the moisture dependency of flow properties, where they can have a strong effect on the speed with which MC changes and the formation of moisture gradients (Angst-Nicollier 2012; Florisson et al. 2019; Toratti 1992). To show the effect regarding the properties obtained in this paper, simulation 1 from the previous section will be compared to simulations where either constant flow properties are assumed or a mixture between constant and moisture dependent properties.

In total, six additional simulations were made, numbered 3–8. Simulations 3 and 4 were performed with a moisture-dependent diffusion coefficient and a constant surface emission coefficient. Simulations 5 and 6 were performed with moisture-dependent surface emission coefficient, but with a constant diffusion coefficient, and simulations 7 and 8 both ran with constant coefficients and therefore did not experience any nonlinearity. All simulations used the constant diffusion coefficient in longitudinal direction presented in Table 4.

For the moisture-dependent properties, the curves in Fig. 10a, c are used. The constant values are based on the minimum and maximum values presented in these curves. The constant surface emission coefficient in simulation 3 had a value of $0.0153\, {\text{m/h}}$ found at $7.35\%$ MC. For simulation 4, the value of $0.007\, {\text{m/h}}$ corresponded to $12.75\%$ MC. The constant diffusion coefficient used in simulation 5 was the value of $1.07\times {10}^{-6}\, ({\text{mm}}^{2}/{\text{h}})$, corresponding to $7.35 \%$ MC. The value in simulation 6 was a value of $2.15\times {10}^{-6}\, ({\text{mm}}^{2}/{\text{h}})$ at $12.75\%$ MC. Finally, simulation 7 runs with the constant values of simulations 3 and 5, where simulation 8 runs with the constant values of simulation 4 and 6.

To show the general difference between constant and moisture-dependent properties, the volumetric average MC with respect to time for the first two RH-phases is presented in Fig. 13 for simulations 1, 7 and 8. The experimental data presented in Fig. 8 showed a difference between the curvatures of the average MC-curve in adsorption and desorption. Figure 13 shows that this behaviour is only visible in simulation 1, where the equilibrium moisture state of the volume is reached sooner in adsorption than in desorption due to the moisture dependency of the properties. This behaviour is not present when constant properties are used as in simulations 7 and 8. The RD-distribution in Fig. 13 shows that in adsorption simulation 8 comes closer to simulation 1, where in desorption simulation 7 comes closer to simulation 1.

Fig. 13

Average volumetric moisture content $\left({w}_{{\mathrm{a}}}\right)$ change in first two RH-phases for simulations 1, 7 and 8 together with a display of the relative difference $({\text{RD}})$

A brief impression of the effect of nonlinearity of the analysis on moisture profiles is given in Fig. 14. The profiles are taken $145\,{\text{h}}$ into the analysis. The results of simulation 1 presented in Fig. 14 agree with the first point illustrated for simulation 1 in Fig. 11. Compared with the surface emission coefficient (simulations 5 and 6), the moisture dependency of the diffusion coefficient (simulations 3 and 4) has the most effect on the flow speed and formation of gradients, with a maximum RD of $8\%$ compared to $0.59\%$.

Fig. 14

Moisture content plots over the width of the cross section at unique nodes at time period ${p}_{1}$$(145\,\mathrm{h})$ illustrated in Fig. 12 for the different simulations in the nonlinearity analysis including estimation of the relative difference $({\text{RD}})$

Conclusion

The simulation tool makes the nonlinear orthotropic single-Fickian model and Neumann boundary condition readily available and the simulations more flexible to different uses. It allowed for the practical incorporation of the scanning curves associated with hysteresis at the boundary and a smoother definition of the moisture-dependent flow properties. The analysis of the calibration made with this tool on volumetric average moisture content data obtained in a bending test subjected to systematic climate changes showed that: (1) the continuous recordings of average moisture content is beneficial due to its simplicity and should be a standard measurement when performing bending tests subjected to climate; (2) the calibration benefitted from the fit made of MC-profiles, where it gives a better representation of moisture content gradient; (3) it is necessary to incorporate hysteresis into the moisture simulation to obtain moisture profiles that best describe the moisture-dependent mechanical behaviour; and (4) for the same reasons, a coherent set of moisture-dependent diffusion and surface emission coefficients are needed to perform a moisture simulation.

The experimental methodology presented in this paper will be used in future publications to analyse mechano-sorption analytically and with simulations. The results obtained in this paper will be used as input to the moisture-induced distortion simulation. The routines can easily be extended by including inhomogeneity of material properties and effects due to temperature on diffusion and surface emission coefficient, making it suitable in conjunction with kiln drying operations. The simulation tool could benefit from a simple hysteresis model that is able to describe hysteresis as a boundary condition in combination with a single-Fickian model based on theoretical equations that are derived from underlying theory. To better incorporate the scanning behaviour and the temperature and moisture dependency of the flow properties, more experiments are needed with different RH-ranges and temperature levels.