Thermo-hygro-mechanically coupled modelling of wood including two-phase moisture diffusion for transient simulation of wooden structures at mechanical and climatic loads

A model to predict the thermo-hygro-mechanical behaviour of wood is introduced. The description of the transport processes of moisture and heat are combined with a model for the mechanical response. Moisture transport is represented by a two-phase multi-Fickian approach, considering bound water and water vapour. For the mechanical response, a moisture- and temperature-dependent, orthotropic, elastic material formulation is used. The theoretical basis of the model and the numerical implementation of the monolithic solution into a Finite Element framework are discussed as well as its verification and validation. With this model at hand, arbitrary wooden structures can be simulated in a transient manner subjected to climatic and mechanical loads. In the contribution, the approach is applied to the analysis of a panel painting by L. Cranach the Elder.


Introduction
Wood is a strongly anisotropic material with moisture-and temperature-dependent properties. Comprehensive descriptions on the physical characteristics of wood can be found in e.g. [1] and [2]. To predict the performance, damage, durability or bearing capacity of wooden structures, models are needed, which are capable to describe the multiphysical transport processes as well as the elastic, plastic, viscous and fracture behaviour. Beside other application fields, especially wooden artwork shall be protected from even small cracks and damage. Damage of wooden artwork is often caused by large changes of climatic conditions, e.g. inappropriately conditioned exhibition or storage locations and during transport. The climatic conditions of the environment have large impact on the mechanical features. Environmental moisture and temperature changes affect the state of moisture and temperature in the material and lead to transport processes in the wood continuum. Moisture is found in wood at three states of aggregation: fluid water and water vapour in the lumens as well as bound water in the cell walls. The phases exhibit different diffusion properties. In most cases, wood is not exposed to fluid water, therefore, many descriptions of water transport do not consider this phase. Inside wood, exchange processes take place between these phases of water. Bound water is the phase that directly influences the structural deformation and mechanical behaviour. Adsorption and desorption of moisture cause swelling and shrinking. Furthermore, moisture affects the mechanical behaviour, e.g. by moisture-dependent elasticity and strength characteristics. In addition, mechano-sorptive creep takes place. For complex structures, transient moisture transport yields extensive deformations and stresses due to moisture gradients and, thereby, different degrees of swelling and shrinking and moisture-dependent mechanical properties in the same sample. Additionally, temperature influences the mechanical behaviour of wood, but mostly indirectly by affecting moisture transport. During sorption, heat is absorbed or emitted. The saturation vapour pressure as well as the equilibrium of bound water and water vapour pressure are temperaturedependent. Also moisture diffusion is influenced by temperature. Beside the indirect influence of temperature on the mechanical behaviour, deformation due to thermal expansion and temperature-dependent elasticity and strength properties is observed, but the influence is relatively small. In complex structures, transient heat transport processes yield temperature gradients, different degrees of thermal expansion in a structure and, thus, stresses and deformations. To simulate the coupling of mechanical behaviour, moisture and heat transport, thermo-hygro-mechanical modelling is required.
For the non-linear numerical analysis of complex structures, the Finite Element Method (FEM) evolved in the last decades. By the spatial discretisation of a structure with an assemblage of elements with simple geometries, approximate solutions can be achieved for the whole complex structure. This method can be applied with respect to the mechanical behaviour as well as to transport processes. It also can be used for coupled analyses, computing the multiphysical solution simultaneously. This is reasonable, since the different physical processes depend on each other.
Up to now, modelling of the mechanical behaviour and transport processes in wood and the simulation by the FEM have been subject of intensive research. Approaches are available considering the coupling of mechanical behaviour and moisture transport of wood, e. g. [ 3,4,5,6]. Many of them include moisture as one single-phase, but [7] shows the inadequacy of single-phase moisture description whenever no instant equilibrium of water vapour and bound water is present. Therefore, a two-phase moisture transport model is discussed by Frandsen et al. in [11] including approaches for phase transition and sorption hysteresis. Based on the model presented by Frandsen et al. a coupled two-phase hygro-mechanical formulation is presented in [11] and applied to the analysis of a panel painting in [12]. In the publications [13] and [13], hygro-thermally coupled moisture and heat transport are described by two moisture phases and temperature as third physical field. In the publication at hand, a new thermo-hygro-mechanically coupled threedimensional model with degrees of freedom for displacements, two moisture phases, i.e. bound water and water vapour, and temperature is presented as well as the implementation into a monolithic solution algorithm within the FEM. It combines the advantages of the hygro-mechanical descriptions in [14] and the multi-physical moisture-and heat-transport formulations in [14] and [8][9][10]. For the mechanical part, elastic behaviour without time-dependent mechanical effects is considered. By coupling of models for two moisture phases and temperature to the mechanical formulation, the presented description is expected to improve the numerical analysis of wood at transient changes of moisture and temperature significantly. Moisture and heat transport are expected to influence the stress-and strain-state significantly so that a thermo-hygro-mechanical simulation is necessary for the damage analysis of climatically loaded wood.
In the present article, modelling of the mechanical behaviour and moisture-and heat-transport is described. At first, the physical background and the mathematical description are introduced and the required model parameters are identified. Subsequently, the algorithmic implementation into the in-house FE-software WoodFEM is introduced. The model is verified and validated by comparing simulation results to several experimental data. Furthermore, the material model is applied to analyse the structure of a panel painting by L. cranach the Elder.
2 Constitutive thermo-hygro-mechanical model for wood

Balance equations
The 3D mechanical sub-model includes three degrees of freedom, i.e. the displacement vector u = u x , u y , u z T . For moisture transport, as mentioned before, two phases are considered, which are represented in the model by the concentration of bound water c b and the water vapour pressure p v . Water vapour pressure correlates with the concentration of water vapour It has to be noted, that Eqs. (1) to (4) do not consider any inertia phenomena and strain rate dependencies, which is reasonable for the analysis of climatically loaded wood.

Elasticity
The total strain tensor = u∕ x is in case of wood typically split additively into with the components of mechanical strain mech , swelling and shrinking m and thermal expansion . The elastic stress-strain-relation is described by hooke's law with the material tensor C . In the paper at hand, only elastic material features are considered for simplicity. In this case, mech equals the elastic strain el . The approach can easily be expanded to plastic and viscous behaviour. The consideration of plastic strain is described in [15], viscous and mechano-sorptive creep are accounted for in [6].
Swelling and shrinking are calculated by with the shrinking tensor and the moisture increment Δm b . Thermal expansion is calculated by with the thermal expansion tensor and the temperature increment ΔT.

Moisture
Based on Fick's first law, moisture flux is given as for bound water diffusion and for water vapour diffusion. In Eqs. (9) and (10), the first term describes the diffusion driven by the gradient of concentration according to Fick's first law. The second term describes the soret-effect, i.e. moisture diffusion driven by the gradient of temperature. This formulation is derived in [10].

Heat
Heat flux is described by Fourier's law as The conduction process with the conduction tensor K is driven by the gradient of the temperature.

Boundary conditions
At the faces of a specimen, boundary conditions allow for moisture and heat exchange. For a realistic simulation of the exchange process, the boundary layer theory (BLT) is used considering the flow properties of the ambient air. A model for mass-and heat-flux is given in [16]. Moisture flux and heat flux through the surface are driven by the difference of water vapour concentration or temperature between the surface and the ambient air. The emission parameters consider the properties of the air by use of the diffusion and conduction in air, D a and K a , the reynoLd-number Re, the schmidt-number Sc and the PrandtL-number Pr, as well as the stream properties of the boundary layer by use of Re.

Algorithmic implementation
To solve the non-linear partial differential balance equations numerically, weak forms of these equations are required. Employing the principle of virtual displacements, applying the divergence theorem and spatial and time discretization leads in combination with the constitutive equations to (16) Fig. 1 Sorption curves for a sequence of moistening and drying (0% RH → 95% RH → 40% RH → 80% RH), comparison of the sorption models of (a) haiLwood-horrobin [14,21], (b) avramidis [6,22], (c) Patera [23], (d) Frandsen [9, 10, 14] In contrast to previous approaches, in the publication at hand, the solution is computed by the monolithic solution of a fully coupled system of all six equations, i.e. three mechanical Eq. 16, two mass balance Eqs. 17 and 18, and one heat balance Eq. 19. For the solution of the non-linear equations, the newton-raPhson-Method is used. Therefore, the derivatives of Eqs. (16) to (19) with respect to the degrees of freedom have to be computed. Linearization yields the overall equation system where the entries of the tangential matrix on the left hand side are given as The term F apl on the right-hand side marks the external mechanical loads. For time integration, a generalised trapezoidal method is used. The index n indicates the values of the current timestep. The index n − 1 stands for the last time-step. By the parameter , the influence of the current and the last timestep on the computation of the increment is described.
For the modelling of wood samples, solid elements with 20 nodes and, for the surface modelling, surface elements with 8 nodes, both based on quadratic shape-functions, are used. By Gaussian quadrature, the equations are integrated numerically over the volume using 27 Gauss points or 9 in the case of surface elements.

Dependencies of physical properties
The coupling of the degrees of freedom, represented in Eq. 20 by the terms K i,j with i ≠ j , arises on the one hand from the coupling terms of the soret-effect, heat convection by moisture transport and the sorption rate in the heat balance equation, and, on the other hand, from material properties, which depend on the bound water concentration, water vapour pressure and temperature. An overview of material parameters and dependencies is given in Table 1.
For the temperature-and history-dependency of the moisture equilibrium and different modelling approaches, see Fig. 1.

Validation of moisture transport
For the validation of moisture transport, an experiment of wadsö [24], see also [8], is simulated. Figure 2 shows the specimen with material directions, dimensions and boundary conditions. Four sides are coated to yield one-dimensional moisture flow. Moisture and heat transfer are possible at the two opposite faces at x = 0 mm and x = 16.2 mm. The transfer is described by surface elements or by prescribing water vapour pressure and temperature at the surface nodes directly. In Table 2, dry density as well as thermal and hygrical expansion coefficients are given.
In the original experiment in [24], the fractional weight increase due to moisture uptake is measured.
Three different moisture steps are conducted in [24]. Since the first step is just done for initial conditioning, only the second step, step B, from 54% RH to 75% RH, and the third step, step C, from 75% RH to 84% RH, are evaluated. Table 3 shows the experimental conditions.
For the simulations, the haiLwood-horrobin-Isotherm [21] is used with the parameters taken from [8].
In the experiment in [24], the fractional weight increase due to moisture uptake is measured. Figure 3 compares the experimental results to the simulation results of the presented model. A good agreement can be seen for both steps, step B and step C. As described in [24], a single phase moisture transport model shows larger discrepancies in step c.
In [8], a moisture diffusion model is validated by comparing the simulation results for the weight increase to the experimental results of [24]. With this model, the moisture distribution in the direction of diffusion is calculated and given for different time steps. Figure 4 compares the simulation results of the moisture distribution [8] and the model used in the publication at hand. In contrast to [8], the climatic conditions are considered by a boundary layer at air velocity v = 3 m s −1 (see [8]) and are not prescribed directly for the boundary nodes. The porosity is supposed to be constant in [8], but dependent on the moisture content in here. As supposed in [11], the slower wetting might result from a decreasing porosity with increasing moisture content. The results at the surfaces are close to the results in [8]. The compared results show differences in the middle of the specimen, which are a consequence of the faster diffusion in [8]. Figure 4b shows the temperature distribution of the experiment. At the beginning, the temperature is rising due to enthalpy gain induced by wetting. When sorption slows down, temperature returns to the temperature of ambient air.

Validation of heat transport
For the validation of heat transport, an experiment described in [25] is simulated. A beech wood specimen is heated and dried by increasing the temperature and decreasing the relative humidity in the ambient air. The specimen with the material directions and boundary conditions is shown in Fig. 5. Dry density as well as thermal and hygric expansion coefficients are given in Table 4.
Solid-elements are used to represent the beech wood specimen and surface-elements to represent the boundary layer. In this simulation, the sorption isotherms of avramidis are used, because the temperature dependency of the sorption is considered. The beech specimen is heated from 20 • C to 120 • C within 60 seconds. Then, this temperature is kept for 45 minutes. The initial moisture content is measured as 8 %, the final relative humidity is 0.1 %. This procedure results in a final moisture content of 0.01 % and an initial relative humidity of 75 %. The climatic conditions are also given in Table 5. For the ambient air, the flow characteristics are not given, but the convection parameters k m = 1 ⋅ 10 −10 kg mm −2 s −1 and k = 15kg s −3 K −1 are documented in [25]. The parameter k v is calculated as During the experiment, the temperature is measured at five locations. The positions of the measuring points can be seen in Fig. 6. To compare the simulation results to the experimental results, the simulation is also evaluated at the five measuring points.
In Fig. 7, the simulations are compared to the experimental results of the temperature at the five indicated measuring points.
Since the simulation yields underestimated temperature using the calculation of the conduction tensor as in [18] and the convection parameters given in [25], these parameters are varied. It should be noted that the coefficients for the conductivity presented in [18] are valid for spruce, which is a softwood, while beech is a hardwood. Additionally, the regression polynomial is valid for temperatures just up to 80 %. In [20], the thermal conductivities are given for beech wood as functions of temperature. With this conductivity parameters, the simulation results are fitting well for point M01 and M03 at the surfaces in longitudinal direction. In the middle of the specimen, the temperature is underestimated. Therefore, another simulation with the same conduction parameters in every material direction is carried out. Here, the simulation results are in good agreement for the middle part of the specimen, but the temperatures for M01 and M03 are overestimated. For the surfaces 0.015 20 ⋅ 10 −6 20 ⋅ 10 −6 3.5 ⋅ 10 −6 Fig. 6 Position of the evaluation points for the experiment of FLeischhauer et al. [25] 0.01 75 0.1 1.76 199 with maximal and minimal y-coordinate, lower temperature is calculated than measured. Due to this fact, it is supposed, that realistic convection parameters have to be larger than given in [25]. In a third simulation, the thermal conductivities from [20] are used, and the convection parameters are fitted. This parameter set yields acceptable results, which are in the range of the results of the previous two simulations.
Since the conduction parameter in longitudinal direction affects mainly points M01 and M03, but rarely the points in the middle, it seems reasonable to assume a smaller conduction parameter in longitudinal direction and a larger convection parameter in radial and tangential direction. It can be seen, that the model is able to predict the temperature distribution versus time well for transient climatic conditions. It is evident, that the identification of material parameters and climatic conditions is important and has large impact on the prediction accuracy. In case of material parameters and boundary conditions, that cannot be clearly identified, data uncertainty has to be taken into account.

Discussion
For a full validation of the presented model, validation experiments would be necessary accounting for changes in moisture and temperature measuring moisture-distribution, temperature-distribution, stress and strain. Up to now, adequate experimental data are missing in the literature. However, a brief example is presented to show the necessity of considering heat transport in the hygro-mechanical model. Figure 8 shows the specimen of the simulated experiment, in Tables 6 and 7 material parameters and boundary conditions of the simulation are given. A 200 mm long piece of wood is loaded climatically by changes of relative humidity and temperature in the environment. Moisture and heat exchange take place at two faces, so that transport processes are one-dimensional in longitudinal direction. A drop of relative humidity from 85% RH to 20% RH and a temperature increase from 10°C to 40°C are simulated. Simulation results of the presented model and a two-phase hygro-mechanical model are depicted in Fig. 9. In the hygro-mechanical simulation, the temperature change is also considered, but at constant temperature distribution in the specimen.
It can be seen, that different moisture distributions result from the two approaches. The hygro-mechanical model predicts faster transport, the thermo-hygro-mechanical model results in higher moisture gradients. As result of the climatic loading, also mechanical stresses occur, see Fig. 9b. The maximum stresses in longitudinal direction along the x-axis are presented. As can be seen, the discrepancy of both approaches is at some points about 70%. It must be noted, that the discrepancy might be higher in more complex structures consisting of several parts with different material properties and grain orientations, since moisture gradients exhibit a larger effect.

Investigation of a panel painting at climate changes
A panel painting of L. cranach the Elder is simulated employing the presented models. The transient thermohygro-mechanically coupled analysis enables to evaluate the influence of climatic loads on complex structures. The chosen panel is the inner side of the left wing of the Katharinenaltar from 1506 (see [12]), shown in Fig. 10.
The panel painting consists of nine bonded basswood (Tilia) panels, coated with a primer, several oil paints and a varnish layer. The on both sides originally coated altar wing has been sawn in an inner and an outer part, so that the thickness is now about 2.6 to 4 mm and the back side is uncoated. On the back, a support construction made of pine has been installed to stabilize the panel and minimise bending deformations. Due to constraining forces, caused by the massive supporting framework, different used wood types and fibre directions, the relatively rigid support construction led to damage of the panel and the paint layers. Therefore, the supporting framework has been removed in 2015 during restoration measures.
For the simulation, solid-elements with 20 nodes and a linear-elastic orthotropic material formulation are used. Since only few properties for basswood are published in literature, the properties of cherrywood, which can be found in [6], are used instead. For the support frame, the properties of spruce are taken instead of pine. This modification is appropriate, because the dry densities are similar and the properties are mainly dependent on the dry density. The coating is included in the front surface layer, modelled by surface-elements with 8 nodes. The finite element model is shown in Fig. 11, where the panel is depicted in grey and the support construction is shown in red and blue.
For the thermal convection parameter, the boundary layer theory is used. The properties of shellac are assumed for the hygric surface resistance at the front surface. The thickness-dependent permeability is given as P t = 3.08989 ⋅ 10 −15 kg m −1 s −1 Pa −1 (see [6]), the supposed thickness of t = 50 μm leads to a permeability of P p = 6.1798 ⋅ 10 −11 kg m −2 s −1 Pa −1 . The hygric convection parameter is calculated using the ideal gas equation as 0.008 23.9 ⋅ 10 −6 34.6 ⋅ 10 −6 3.5 ⋅ 10 −6  The climatic boundary conditions are shown in Fig. 12. Three climatic scenarios have been simulated, one cyclic humidity and temperature change and two single humidity changes at cyclic temperature change. For climate 1, the relative humidity, beginning at RH = 65 %, changes alternating between 50 % and 80 %, simulating day night fluctuation. Temperature starts at 20 % and alternates with an amplitude of 4 %, while the temperature is overall increasing by 1 K per day. Climate 2 is a single change from RH = 65 % to RH = 80 % . In climate 3, relative humidity also starts at RH = 65 % and changes to RH = 50 %. Figure 12c and d show the results for the average concentration of bound water and the deflection of the middle of the panel painting. The change of the bound water concentration and the deflection follow the changes of the climatic conditions. It can be seen, that the moisture content does not reach equilibrium within 12 hours. The graph of the moisture just flattens when approaching equilibrium. Due to the overall increase of temperature, moisture is decreasing overall. The deflection results from the different material directions of the panel and the support construction. The x-direction equals to the tangential direction in the panel, but to the longitudinal direction in the support construction. Due to the different hygro-expansion coefficients of the material directions, much larger swelling and shrinking arises in the panel compared to the support beams in x-direction. Figure 12 shows also the deflection due to the single humidity changes of climate 2 and climate 3. It reveals, that also with constant relative humidity in the ambient air, oscillations result from the temperature changes, due to the temperature-dependent sorption equilibrium. Furthermore, constant relative humidity together with changing temperature yield changing water vapour pressure. Figure 13 shows the deformed panel painting as well as temperature and moisture distribution after 118 h (dashed line in Fig. 12a and b), when a change from dry, warm to wet, cool climate takes place It can be seen, that the dependency of the temperature and moisture in the structure on the overall thickness of the panel is represented by the introduced approach. At the areas, where no support frame is connected to the panel, the cross-section is the thinnest, therefore, cooling is fastest. Most heat is buffered at the thickest parts of the panel. The back of the panel cools down faster than the front, due to the larger surface. Also, the moisture distribution depends on the thickness. Here, the overall drying shows the effect, that the points at the thickest cross-sections of the panel painting on the coated side are the wettest, due to the small permeability of the surface and the long diffusion distance through the material. The overall drying has to be considered instead of local change, because In conclusion, it can be stated, that temperature influences moisture transport significantly and, thus, also the mechanical response of the panel painting. Sorption and diffusion are temperature-dependent processes. Simultaneously, moisture exchange affects temperature. Therefore, thermo-hygromechanical modelling seems to be reasonable and relevant.

Summary and conclusion
A new thermo-hygro-mechanically coupled monolithic FEM algorithm with multi-Fickian moisture transport is presented, which links models for two phase moisture transport, heat transport and mechanical behaviour of wood. The specific approaches are introduced and combined by an overall equation system for the monolithic computation of mechanical, thermal and hygrical states. The model is validated by one moisture transport experiment and one heat transport experiment. Both simulations show a good agreement with the experimental results. Therefore, it can be stated, that the model is able to well predict the transport processes individually. For the validation of coupling of transport processes and the mechanical response, experimental data are missing. But it can be seen, that results of thermo-hygro-mechanical simulations differ from results of hygro-mechanical simulations at combined moisture-and temperature load. Heat transport affects the moisture distribution and by this results in higher climate-induced stresses. Since the thermo-hygromechanical model considers in reality appearing heat transport, it is supposed to yield more realistic results than the hygro-mechanical model.
To illustrate a possible application of the introduced model, a panel painting of L. cranach the Elder is simulated. Daily cycles of relative humidity and temperature are computed as climatic load. Moisture and temperature distribution as well as the deformation response of the panel painting are presented. The application of the model shows, that it is able to analyse also relatively complex structures at changing climatic conditions. The simulation results are reasonable. It is possible to estimate critical climatic conditions by evaluation of the mechanical response.
A huge need for experimental investigation of material properties can be seen due to large scattering of wood properties depending on the wood species, location and conditions of growth, moisture state, temperature, age and history. Furthermore, the uncertainty of wood properties (see e.g. [26]) has to be taken into account numerically as well. Moreover, the formulation can be extended by more phases, e.g. fluid water as in [27]. Finally, the application of further mechanical descriptions, like thermo-hygromechanically coupled fracture, plasticity (e.g. [15]), viscous creep and mechano-sorption (e.g. [28]) approaches has to be investigated.