A 3D multi-scale hygro-mechanical model of oak wood

Consider a composite material with a heterogeneous fine-scale structure, see Fig. 2b. The coarse-scale three-dimensional domain \(\Omega\) is constructed as a periodic repetition of the fine-scale unit cell Q depicted in Fig. 2a. Denoting \(\mathcal {L}\) and \(L=\eta \mathcal {L}\) as the characteristic lengths of the coarse- and fine-scale domains, the requirement of a strong separation of scales is represented by the condition \(\eta \ll 1\). The field quantities (displacement, strain, and stress) governing the hygro-mechanical behavior vary smoothly at the coarse scale, while they are periodic at the fine scale. Consequently, the field quantities may be assumed to explicitly depend on two variables: A slow variable \({\textbf {X}}\) that varies at the coarse-scale level and a fast variable \({\textbf {x}}={\textbf {X}}/\eta\) that varies at the fine-scale level. In the asymptotic homogenization procedure, which is based on the framework developed in Bosco et al. (2017a, 2017b, 2020a), the heterogeneous material at the fine scale is modeled as an equivalent homogeneous solid at the coarse scale. The procedure starts from the mechanical equilibrium equation, which, in the absence of body forces, is expressed bywhere \(\varvec{\sigma }\) is the Cauchy stress tensor that satisfies the hygro-elastic constitutive relation:Here, \({\textbf {u}} ({\textbf {X}})\) is the displacement field and \(\Delta \text {RH}({\textbf {X}})\) is the ambient, coarse-scale relative humidity variation, which is equal to the difference between the current relative humidity \(\text {RH}\) and the initial reference relative humidity \(\text {RH}_0\), i.e., \(\Delta \text {RH} = \text {RH} - \text {RH}_0\). The material properties at the fine-scale level determine the coarse-scale response and are represented by local constitutive tensors that are periodic functions of the fine-scale variable \({\textbf {x}}\), i.e., the elastic stiffness tensor \(^4{\textbf {C}}={^4{\textbf {C}}({\textbf {x}})}\) and the hygroscopic expansion tensor \(\varvec{\beta } = \varvec{\beta }({\textbf {x}})\). Although not explicitly indicated in Eq. (2), the local stiffness tensor \(^4{\textbf {C}}\) and the hygroscopic expansion tensor \(\varvec{\beta }\) can also depend on the current value of the relative humidity \(\text {RH}\). Note that the hygroscopic strain is commonly defined in terms of a moisture content variation, instead of a relative humidity variation as done in Eq. (2). The formulation in Eq. (2) is imposed by the homogenization procedure, which requires defining the effective properties as a function of a uniform field over the fine-scale domain. The moisture content depends on the water adsorption capacity of each of the constituent phases, and thus is non-uniform across the fine-scale configuration. Conversely, the relative humidity is uniform across the scales, and therefore is used as the reference field for expressing the hygroscopic strain in Eq. (2).

The displacement field is formulated as an asymptotic expansion in terms of \(\eta\):Expression (3) is inserted into the constitutive relation (2), after which the result is substituted into the equilibrium equation (1), leading to a set of boundary value problems—the so-called “cell problems”—defined at the fine scale, i.e.,where \({}^{4}{} {\textbf {I}}^{S}\) is the fourth-order symmetric identity tensor, given by \(I^S_{ijkl} = (\delta _{il}\delta _{jk} + \delta _{ik}\delta _{jl} )/2\), with \(\delta _{ij}\) the Kronecker delta symbol. The solution of these cell problems provides the fine-scale influence functions \(^3{\textbf {N}}_{1}\) and \({\textbf {b}}_1\), which represent the periodic fluctuations of the displacement field within the fine-scale domain, as caused by respective variations of the coarse-scale deformation and the relative humidity.

The cell problems in Eqs. (4) and (5) are completed by periodic boundary conditions expressed in terms of the unknown influence functions \(^3{\textbf {N}}_{1}\) and \({\textbf {b}}_1\). Periodic boundary conditions require that the values assumed by each influence function are equal at corresponding points on opposite boundaries (left and right, top and bottom). In addition, for warranting the uniqueness of the solution, the influence functions are prescribed to vanish on average across the dimensions of the unit cell. Finally, the cell problems are solved numerically by the finite element method, thereby considering a specific material configuration that represents the fine-scale domain. Once the influence functions \(^3{\textbf {N}}_{1}\) and \({\textbf {b}}_1\) have been obtained from Eqs. (4) and (5), the effective elastic tensor \({}^{4}\bar{{\textbf {C}}}\) and hygroscopic expansion tensor \(\bar{\varvec{\beta }}\) at the coarse scale can be computed asFurther details on the derivation of the effective properties, the solution of the cell problems, and the finite element implementation of the method can be found in Bosco et al. (2017a) and Yuan and Fish (2008).

Consider a three-dimensional laminate \(\Omega\) composed of a finite number of layers. Each layer k is characterized by uniform material properties, defined through the elastic stiffness tensor \(^4{{\textbf {C}}}_{k}\) and the hygroscopic expansion tensor \(\varvec{\beta }_k\). The effective properties of the laminate are obtained from a Voigt averaging procedure, which assumes full kinematic compatibility between the composing layers. Hence, the strains in the individual layers are supposed to be equal, and thereby are equal to the effective strain in the laminate. Furthermore, the effective stress in the laminate is consistently determined by the average of the stresses in the individual layers weighted by the layer thicknesses. The effective elastic stiffness tensor \(^4{\bar{{\textbf {C}}}}\) and the effective hygroscopic expansion tensor \(\varvec{\bar{\beta }}\) of the multi-layer laminate then follow as (Gao and Oterkus 2019)with \(t_k\) the relative thickness of layer k, as determined with respect to the total thickness of the laminate.

The meso-scale model is based on high-resolution microscopy images of an oak wood sample. For obtaining the images, a cubic oak sample was cut from a larger specimen along the three anatomical (L, T, R) planes, which refer to the longitudinal (L), transverse (T), and radial (R) material directions of wood. Using a disposable blade and a sliding Reichert microtome, three sections (i.e., one per principal material direction) were cut and subsequently prepared into high-quality permanent slides containing several growth rings. The surface of the slides was equal to \(14\times 8.5\) \(\hbox {mm}^2\), and the thickness ranged between 15 and \(25\,\upmu\)m. A solution of 1 g of Safranin powder in 100 ml of distilled water was prepared for staining the thin slides. After staining, the slides were washed with water to remove the residual stain, until the running water was clear. Next, the slides were dehydrated with an ethanol concentration that was gradually increased in a stepwise fashion from \(50\%\) to \(75\%\), and finally to \(96\%\), in order to prevent damage of cell walls under rapid shrinkage (Yeung et al. 2015). The slides were embedded in Canada balsam under a cover glass and were dried for a period of 12 to 24 hours at a temperature of \(60\,^\circ\)C (Yeung et al. 2015). Subsequently, the permanent slides were scanned using an automated digital microscope \(100\times\) to create high-quality images of the complete slides.

Figure 8a illustrates the transverse (R–T) section of the oak sample. The image reveals the complex cellular structure of oak, in which libriform fibers, fiber tracheids, axial parenchyma, vessels, and rays are present. Further, three different growth rings A, B, and C with width W and height H have been selected for the construction of meso-scale unit cells. The variations in morphological features within the growth ring are clearly visible. Specifically, a pattern of large earlywood vessels can be identified across the width \(W_e\) of the earlywood region of growth ring B. In contrast, much smaller vessels are present across the width \(W-W_e\) of the latewood region, resulting in a relatively dense material structure. Moreover, parallel regions of multi-seriate and uni-seriate rays between the axial cells can be detected. Figure 8b, c show the radial (R–L) and tangential (T–L) sections of the oak sample, respectively. Here, the large voids with an elongated oval shape correspond to vessel sections. Additionally, the relatively dark regions indicate fiber sections, while the light regions correspond to parenchyma cells. It can be further observed that the cross section of rays in the R-T plane is substantially different from that in the T-L plane, in correspondence with respective cuts made along the longitudinal and transverse axes of the ray cells. The resolution of the original images is very high so that cell walls can be accurately identified when zooming in at higher magnification. The high resolution enables an accurate image processing required for the construction of three-dimensional finite element configurations of oak growth rings at the meso-scale, as described in Sect. “Image–based meso–scale unit cell” below.

The meso-scale unit cells are generated based on the microscopy images presented in Fig. 8. Accordingly, three regions of dimensions \(W\times H\) are selected in the R-T plane of the meso-scale material structure, as indicated by the black rectangular boxes A, B, and C in Fig. 8a. Note that the three regions correspond to three distinct growth rings, where each growth ring is characterized by a specific proportion and distribution of the various cell types, and thus has a specific local density, denoted by \(\rho _{\text {A}}, \rho _{\text {B}}\), and \(\rho _{\text {C}}\). Under the assumption of a periodic meso-structure, see Fig. 2, the three selected regions may be considered as representative of three distinct oak samples that are characterized by the overall densities \(\rho _{\text {A}}, \rho _\text {B}\), and \(\rho _{\text {C}}\). Each of the three two-dimensional regions extracted from the R-T plane image is next converted into a binary image with an appropriate pixel intensity threshold. The threshold is fine-tuned such that the black pixels of the binary image represent cell wall material, as shown in Fig. 9a. An extra homogeneous layer with a thickness of three black pixels is added at the rectangular outer boundary of each of the three regions to impose geometrical periodicity on the domain. Subsequently, the cell wall edges are identified using the image processing toolbox for boundary detection of MATLAB, which contains an edge detection scheme based on localizing the variations of gray-level images. Since the resolution of the obtained images is insufficient to provide smooth edges for the cell walls, jagged-shaped edges for the black pixel regions are retrieved instead, as illustrated in Fig. 9b. From the jagged geometry, the smooth edges of the cell walls are reconstructed using a second-order B-spline interpolation function, for which the corner points along the interface between the black and white pixels are selected as the control points of the B-spline function, as indicated by the solid circles in Fig. 9b. The smooth cell edges identified with this procedure are provided as input for the generation of a finite element mesh. Figure 9c shows a detail of the two-dimensional finite element mesh based on triangular elements. From the discretized two-dimensional domain, the three-dimensional FEM mesh can be subsequently constructed. In the literature, three-dimensional FEM meshes for wood are generally established through the extrusion of the two-dimensional mesh in the out-of-plane direction (Malek and Gibson 2017; Mendez et al. 2019). However, when extruding the current 2D image in the (longitudinal) out-of-plane direction, the section of rays becomes fully “open” along the longitudinal direction, see the local geometry within the dashed rectangle shown in Fig. 9e, which may lead to an underprediction of the effective, macro-scale Young’s moduli, especially the effective Young’s modulus \(\bar{E}_T\) in the tangential direction (Malek and Gibson 2017). To minimize this discrepancy, the following strategy is proposed for generating the three-dimensional meso-scale domain. First, the two-dimensional image is extruded along the longitudinal (L) direction, see Fig. 9e, where the thickness of the extruded volume is set equal to the average thickness of rays in the longitudinal direction, as determined from the image in the T-L plane presented in Fig. 8b. Next, in the regions containing rays, the domain is closed in the longitudinal direction by adding elements with ray cell wall properties in the two outer element layers, see Fig. 9f. The three-dimensional unit cell obtained is depicted in Fig. 9g, and thus reflects the presence of single ray cells in the longitudinal direction. Upon extrusion, the triangular finite elements initially defining the 2D mesh in the R-T plane are replaced by linear, six-node wedge elements equipped with a two-point Gauss integration scheme to construct the 3D mesh. Depending on the material density, the 3D finite element configurations contain between \(3.2 \times 10^6\) and \(4.0 \times 10^6\) elements. A mesh-refinement study not presented here has shown that the constructed meshes are sufficiently fine for providing mesh-independent simulation results.

Table 5 summarizes the characteristics of the generated three-dimensional meso-scale unit cells. The width W and height H of the unit cell and the width \(W_e\) of the earlywood region are measured from the microscopy image in the transverse plane, as indicated in Fig. 8a for sample B. Further, the porosity \(\phi\), the ray volume fraction \(v_\text {ray}\) and the vessel volume fraction \(v_{\text {vess,EW}}\) in the earlywood (EW) of samples A, B, and C are determined by measuring the corresponding regions in the discretized three-dimensional meso-scale unit cell. The porosity values range between \([0.45-0.60]\); assuming a cell wall density in the dry state of \(\rho _{cw} = 1440\) kg/\(\hbox {m}^3\) (Kellogg and Wangaard 1969), the densities \(\rho = (1-\phi )\rho _{cw}\) of the generated meso-scale unit cells follow as \(\rho _\text {A} = 582\) kg/\(\hbox {m}^3\), \(\rho _\text {B} = 685\) kg/\(\hbox {m}^3\), and \(\rho _\text {C} = 798\) kg/\(\hbox {m}^3\). These densities are in good agreement with the densities reported in the literature for Pedunculate oak, which vary in the range \([461-874]\) kg/\(\hbox {m}^3\) (Badel and Perré 2006). The solid volume fraction of ray cells is defined as \(v_\text {ray} = V_\text {ray}/(V_\text {ax}+V_\text {ray})\), with \(V_\text {ray}\) and \(V_\text {ax}\) being the volumes of axial cells and ray cells. From the value of \(v_\text {ray}\), the solid volume fraction of axial cells, \(v_\text {ax} = V_\text {ax}/(V_\text {ax}+V_\text {ray})\), is simply obtained as \(v_\text {ax}=1-v_\text {ray}\).

The macroscopic elastic response of oak following from the three-level homogenization procedure is approximately orthotropic and thus is represented by the nine independent elastic parameters shown in Figs. 13, 14, and 15. In these figures, solid, dashed, and dotted lines refer to the meso-scale unit cells A, B, and C, which are respectively characterized by a relatively low density (\(\rho _\text {A}=582\) kg/\(\hbox {m}^3\)), a medium density (\(\rho _\text {B}=685\) kg/\(\hbox {m}^3\)), and a high density (\(\rho _\text {C}=798\) kg/\(\hbox {m}^3\)), see also Table 5. Further, the results are depicted for MFAs of \(10^\circ\) and \(20^\circ\) in the S2 layer of axial cells.

Figure 13 illustrates the effective Young’s moduli of oak as a function of the macro-scale moisture content \(\bar{m}\), together with the experimental data reported in Ozyhar et al. (2016), Güntekın et al. (2016a), Luimes et al. (2018b) and Güntekın et al. (2016b). The Young’s moduli in the radial direction \(\bar{E}_R\), tangential direction \(\bar{E}_T,\) and longitudinal direction \(\bar{E}_L\) are presented in Fig. 13(a), Fig. 13(b), and Fig. 13(c), respectively, and clearly show that a higher material density generally leads to a higher stiffness value. Specifically, the Young’s modulus \(\bar{E}_L\) increases almost proportionally with the material density \(\rho\), as also reported in other modeling works on hardwood (Malek and Gibson 2017). Additionally, all three Young’s moduli show a slightly increasing trend at relatively small moisture content, followed by a monotonic decrease for higher moisture content values. The initial increase in the curves is caused by the moisture-dependent behavior of lignin at the nano-scale, see Fig. 4(c). Despite that the experimental data in Fig. 13 represent some scatter, they lie close to the model predictions and show a comparable decreasing trend under a growing moisture content. It can be further observed that the Young’s modulus \(\bar{E}_L\) in the longitudinal direction has the largest value, and is less influenced by moisture content variations than the Young’s moduli \(\bar{E}_R\) and \(\bar{E}_T\) in the radial and tangential directions. This feature can be ascribed to the predominantly longitudinal orientation of the stiff, hydrophobic cellulose micro-fibrils in the S2 layer of axial cells, i.e., MFA = \(10^\circ , 20^\circ\), see Table 3. Note hereby that an increase in the MFA orientation from \(10^\circ\) to \(20^\circ\) in the S2 layer of axial cells indeed leads to a substantial drop in the Young’s modulus \(\bar{E}_L\) in the longitudinal direction, and has a negligible influence on the Young’s moduli \(\bar{E}_T\) and \(\bar{E}_R\) in the transverse and radial directions. Figures 13(a) and 13(b) also illustrate that the radial Young’s modulus \(\bar{E}_R\) generally has a higher value than the tangential Young’s modulus \(\bar{E}_T\). This feature can be ascribed to the prevalent alignment and the denser configuration of fiber cells in the radial direction compared to the transverse direction and the stiffening effect of the rays in the radial direction (Perré and Badel 2003). At the highest moisture content considered, \(\bar{m} = 21\%\), the ratio \(\bar{E}_R/\bar{E}_T\) ranges between 2.16 and 2.32, with the exact value depending on the material density. This result is comparable to the experimental value \(\bar{E}_R/\bar{E}_T \approx 2\) reported in Burgert et al. (2001) for a fully saturated Pedunculate oak wood sample with a dry density of 530 kg/\(\hbox {m}^3\).

Figure 14 shows the effective Poisson’s ratios of oak as a function of the macro-scale moisture content \(\bar{m}\), together with the experimental data reported in Ozyhar et al. (2016), Güntekın et al. (2016a), Luimes et al. (2018b) and Güntekın et al. (2016b). The influence of the material density \(\rho\) on the Poisson’s ratios \(\bar{\nu }_{LT}\) (Fig. 14(a), tangential plane), and \(\bar{\nu }_{RT}\) (Fig. 14(c), transverse plane) is mild, and on the Poisson’s ratio \(\bar{\nu }_{LR}\) (Fig. 14(b), radial plane) is negligible. Further, in correspondence with the experimental results presented in Mendez et al. (2019), the Poisson’s ratios only slightly increase under an increasing moisture content, with the exception being the somewhat more substantial increase for \(\bar{\nu }_{LT}\) in case the S2 layer of axial cells is characterized by the larger micro-fiber angle of \(20^\circ\), see Fig. 14(a). The effect on the Poisson’s ratio by a change in MFA from 10\(^\circ\) to 20\(^{\circ }\) is the largest for \(\bar{\nu }_{LT}\), followed by \(\bar{\nu }_{LR}\) and \(\bar{\nu }_{RT}\). Despite their spread, the experimental data for the Poisson’s ratios show a consistent, approximately constant trend under a varying moisture content, for which the values in most cases are comparable to the modeling results obtained for an MFA of \(20 ^\circ\) in S2 layer of axial cells.

Figure 15 finally shows the effective shear moduli of oak as a function of the macro-scale moisture content \(\bar{m}\), together with experimental data taken from the literature (Ozyhar et al. 2016; Güntekın et al. 2016b; Luimes et al. 2018b). The shear moduli associated with the tangential plane \(\bar{G}_{LT}\), the radial plane \(\bar{G}_{LR}\), and the transverse plane \(\bar{G}_{RT}\) are presented in Fig. 15(a), Fig. 15(b), and Fig. 15(c), respectively. Similar to the Young’s moduli, the shear moduli increase with increasing material density \(\rho\). Further, all shear moduli show a comparable trend under a growing moisture content, which is characterized by a small, initial increase, followed by a substantial, monotonic decrease at larger moisture content. Additionally, the MFA orientation in the S2 layer of axial cells does not influence the effective shear modulus \(\bar{G}_{RT}\) in the transverse plane, while the moduli \(\bar{G}_{LT}\) and \(\bar{G}_{LR}\) in the tangential and radial planes become larger for a larger micro-fiber angle. This result has also been reported in Astley et al. (1998). Observe further that the (limited) experimental data points lie relatively close to the modeling results, especially to those related to the larger micro-fiber angle of \(20^\circ\).

The effective hygroscopic coefficients of oak following from the three-level homogenization procedure are depicted in Fig. 16 as a function of the macro-scale moisture content \(\bar{m}\), where the solid, dashed, and dotted lines refer to the meso-scale unit cells A, B, and C with densities \(\rho _\text {A}=582\) kg/\(\hbox {m}^3\), \(\rho _\text {B}=685\) kg/\(\hbox {m}^3\), and \(\rho _\text {C}=798\) kg/\(\hbox {m}^3\), respectively. The results correspond to an MFA of 10\(^\circ\) and 20\(^\circ\) in the S2 layer of axial cells. The hygroscopic behavior of oak is approximately orthotropic and thus is essentially defined by the three coefficients \(\bar{\beta }_\text {R}\), \(\bar{\beta }_\text {T}\) and \(\bar{\beta }_\text {L}\) depicted in Fig. 16(a), Fig. 16(b) and Fig. 16(c), respectively. These coefficients represent the diagonal elements of the hygroscopic expansion tensor; it has been indeed confirmed that the computed off-diagonal elements are negligible. The experimental results presented in Fig. 16 are based on the test data presented in Ozyhar et al. (2016), Badel and Perré (2006) and Kránitz (2014); since the values of the moisture content at which the hygro-expansion coefficients were measured are not reported in these references, only the spread in experimental data has been displayed by means of a gray region of constant bandwidth. Note that the experimental bandwidths for the three hygroscopic coefficients indeed mostly fall within the ranges of the modeling results. It is further seen that the hygro-expansion coefficient is hardly affected by both the material density and the MFA, with the latter result being confirmed by the experimental-numerical study presented in Yamamoto et al. (2001). In addition, the hygro-expansion coefficients typically decrease as a function of the moisture content; only for the longitudinal component \(\bar{\beta }_\text {L}\) a slight initial increase caused by the typical moisture-dependent behavior of lignin can be observed. Note further that the hygro-expansion coefficient in the transverse direction, \(\bar{\beta }_{T}\), is larger than the hygro-expansion coefficient in the radial direction, \(\bar{\beta }_{R}\), for which, depending on the material density, the ratio \(\bar{\beta }_{T}\)/\(\bar{\beta }_{R}\) lies between 1.69 to 1.90. This result corresponds well with experimental results presented in Badel and Perré (2006) for oak samples with an average density of 700 kg/\(\hbox {m}^3\), for which \(\bar{\beta }_{T}\)/\(\bar{\beta }_{R} \approx 1.78\). Finally, the hygro-expansion coefficient in the longitudinal direction is much smaller than in the two other directions, which is a typical characteristic of (oak) wood (Ozyhar et al. 2016; Kránitz 2014).