Modeling and design of thin bending wooden bilayers

In recent architectural research, thin wooden bilayer laminates capable of self-actuation in response to humidity changes have been proposed as sustainable, programmed, and fully autonomous elements for facades or roofs for shading and climate regulation. Switches, humidistats, or motor elements represent further promising applications. Proper wood-adapted prediction models for actuation, however, are still missing. Here, a simple model that can predict bending deformation as a function of moisture content change, wood material parameters, and geometry is presented. We consider material anisotropy and moisture-dependency of elastic mechanical parameters. The model is validated using experimental data collected on bilayers made out of European beech wood. Furthermore, we present essential design aspects in view of facilitated industrial applications. Layer thickness, thickness-ratio, and growth ring angle of the wood in single layers are assessed by their effect on curvature, stored elastic energy, and generated axial stress. A sensitivity analysis is conducted to identify primary curvature-impacting model input parameters.


Introduction
Wood, as a sustainable and natural-grown fiber composite material, has been used over many centuries as mankind's research on wood and wood-based products has enabled countless applications. More recently, novel and innovative concepts for application as self-actuationcapable, humidity-responsive structures have arisen in applied research for biomimetic architecture [1][2][3][4][5][6][7][8]. Bi-layered wooden structures, capable of complex and extensive shape changes in function of a given geometrical setup and change in ambient climate have been envisaged for diverse uses including shading elements [7], climate-adaptive facades [9], or motor elements [1]. The promising and sustainable principle found its inspirations in nature where anisotropic biological materials with inherent bi-layered and differential fiber structure use humidity changes to generate movement [10][11][12].
Transferring such innovative concepts to industrial production standards with essential economic benefits requires high reliability of the performance of the elements according to given design and specific application. Next to this, the accurate prediction of actuation as a function of geometry and climate change is crucial. A model was derived by Timoshenko [13] PLOS ONE | https://doi.org/10.1371/journal.pone.0205607 October 16, 2018 1 / 12 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 for predicting bending of thin bi-metal strips as a function of temperature change. This concept has recently been applied to the material wood as a linear-elastic prediction of thin bending wooden bilayers [1]. However, the consideration of the essential physical characteristics of the material wood for such simple models is missing. Wood, a biological, anisotropic and hygroscopic material shows a strong moisture-dependency of its material parameters that should be incorporated when modeling its mechanical behavior [14,15]. Furthermore, as a consequence of orthotropic material behavior, wood should not be treated as isotropic when reduced to two dimensions (2D) for modeling. We present a simple, wood-adapted, 2D linear-elastic model for predicting curvature for a given moisture content change of thin bending bilayers, based on Refs. [1,13]. The model is validated using experimental data collected on cross-grained bilayers made out of European beech wood. Using the model, we assess geometrical design aspects such as thicknesses, thickness-ratio, or growth ring angle of the wood in single layers in function of available climate and wood species. These insights help in avoiding potential secondary problems such as adhesive-bond delamination or increased deformation prediction errors. Furthermore, we analyze the effects of these parameters on resulting curvature, stored elastic energy, and maximum axial elastic stresses. Important dimensioning characteristics are suggested in view of optimized structures for industrial application.
High variability in material parameters, a natural characteristic inherent to wood [16], in combination with fabrication tolerances resulting from wood-processing, affecting bilayer geometry, may considerably influence resulting deformation magnitude and thus the precision of the model. To assess these uncertainties, we make use of total Sobol' indices [17] for a sensitivity analysis on each model input parameter. Through partial variances, the impact of variability in each input parameter on the variability of the model output, i.e., bilayer curvature, can be quantified. Critical parameters can be identified and separated from less important ones in view of efficient in-situ sample measurements and model parameter-updating for facilitated industrial processes.

Elastic properties of European beech
Considering Hooke's law in three dimensional Euclidean space σ = Cε, a strain-state (ε = ε kl e k � e l ) is mapped onto a stress state (σ = σ ij e i � e j ) by a 4 th -order stiffness tensor C (C = C ijkl e i � e j � e k � e l ). The compliance tensor S (S = C −1 ), for an orthotropic material such as wood, is expressed in terms of the engineering constants (w.r.t a basis e 1 , e 2 , e 3 ) and in Voigt-notation as where E i (Young's moduli), G ij (shear moduli), and ν ij (Poisson's ratios) can be collected from mechanical characterization experiments. In the case of wood, the local anatomical growth directions L, T, and R (Longitudinal or grain direction, tangential, and radial direction) can be assigned in arbitrary order to the basis e 1 , e 2 , e 3 . Hereby, the representative elementary volume is considered far away from the pith, i.e., the growth-ring-curvature is neglected. Elastic material properties have been collected for the hardwood species European beech (Fagus sylvatica L.), in L, R, and T directions as functions of wood moisture content ω in Ref [14]. A linear dependence on ω is assumed for the species beech as P i = b 0 + b 1 ω, where P i stands for the i th property P that represent the engineering constants, and b 0 and b 1 are coefficients found in Ref [15].

2D linear elastic model for wooden bilayer curvature
The 3D elastic material law ε = Sσ with S as described in Eq (1) is hereafter reduced to 2D by assuming (i) a plane strain state, and alternatively, (ii) a plane stress state in. A 2D schematic of a wooden bilayer is represented in Fig 2, where the undeformed state (with initial wood moisture content ω 0 ) is shown along a deformed state (curvature κ, wood moisture content ω 1 < ω 0 , for drying). Typically, the wood is oriented such that in layer 1, L aligns with e x and in Properties P i as function of moisture ω as found in Refs [14] and [15]. The moisture interval [8%, 18%] roughly corresponds to a relative humidity range of [35%, 85%] (adsorption at 85% and desorption at 35%).
https://doi.org/10.1371/journal.pone.0205607.g001 layer 2 R aligns with e x and T with e y . This results in a stiff and resisting layer 1 (called "passive" layer, very low swelling in L direction) to block axial deformation at the interface (an adhesive bond) coming from the driving layer 2 (the "active" layer, high swelling and shrinkage in R-or T-direction). For the resulting curvature κ of a 2D model, the stiffnesses Q i,x in direction of e x of each layer (i = 1, 2) are of interest. For cases (i) and (ii), and using the same basis as for Eq 1, they are derived as assuming a uni-axial stress state σ (2D) = (σ x , 0, 0) T to act on the bilayer. The curvature for a cross-grained wooden bilayer structure is a function of the applied difference in moisture content. Existing models assume a direct proportionality κ / Δω [1]. The presumed proportionality factor (referred to as K) is, in the case of wood, not constant but is itself a function of moisture (K = K(ω)). For any bilayer configuration, K contains the effects of geometry and material parameters and is obtained by formulating equilibrium, compatibility, and a deformation equilibrium at the interface of two bonded Euler-Bernoulli beams [13]. In terms of Eq 2 and considering the basis e x ,e y in Fig 2, the term is derived as where are the second moments of area per unit width of bilayerstrip. h 1 and h 2 denote layer thicknesses and α 1,x and α 2,x are the differential swelling coefficients in bilayer axial direction (e x ) such that α = ε ω Δω −1 where ε ω is the swelling strain. For beech wood, typical values used are α L = 0.0001, α R = 0.0019, and α T = 0.0040 in % −1 [15]. Here, we assume α as remaining constant over O, otherwise, the swelling strain would read ε ω = R O αdω. O stands for a moisture domain such that for beech wood O � [0, FSP(%)] contains all possible states of ω where FSP denotes the fiber saturation point. For any moisture increment Δω 2 O, the change in layer-thickness (direction e y ) of layer i can thus be calculated as Δh i = h i α i,y Δω, and added (swelling) or subtracted (shrinkage) to the actual layer thicknesses. Next to updating thicknesses and mechanical parameters over changes of ω 2 O, the growth ring inclination φ, considered only in layer 2 as φ 2 (Fig 2) affecting the orientation of the local wood-RT-plane, results in a rotation of the properties (analogous to an elementrotation in Mohr's circle). In layer 2, α is thus transformed as α 2,x = α 2,11 cos 2 (φ 2 ) + α 2,22 2D schematic of wooden bilayer. Initial flat (wet, ω 0 ) and resulting deformed state (dry, ω 1 ) after drying (ω 0 > ω 1 ). Global 2D coordinate system (basis e x , e y ) and local anatomical directions (R, T, L) for layers 1 and 2 with growth ring inclination φ 2 of layer 2. Passive layer (layer 1) with thickness h 1 and active layer (layer 2) with thickness h 2 . Stiffnesses Q i,x in e x -direction of layer i and axial stress distribution σ x over cross-section.
https://doi.org/10.1371/journal.pone.0205607.g002 sin 2 (φ 2 ), and equally, the transformed stiffness, where Q 2,11 is calculated as in Eq 2, reads Considering all moisture dependencies, thickness changes of layers and growth ring inclination, the resulting curvature κ over any O is calculated as The total work performed, or resilient elastic energy stored, W ([W] = Nm � m −2 ), to achieve κ over O can be expressed, in the 2D case, as The factor F i stands for contribution from layer i. For a bilayer with passive layer 1 and active layer 2, If, as in the calculation of Eq 3, it is assumed that the neutral axes situate at half the thicknesses of each layer and that stiffnesses Q i are constant over cross-section, the axial maximal and absolute stress values in layer i at layer-interface and layer-edge (in e x -direction) can be expressed by We note that h, I, Q, K, W, F, and σ are all functions of ω, and that in a strict sense, the integrals of Eqs 5-7 can be evaluated for any given Δω 2 O. However, for practical application, it is convenient to replace them using an incremental approach. For Eq 5, e.g. k À k 0 � P n i DK i Do i for Δω i = O/n, where in each increment, the geometry and engineering constants (all E, G, and ν) are updated for the actual moisture level.

Bilayer validation samples
Six different configurations of cross-grained bilayers were fabricated from beech wood (Fagus sylvatica L.) grown in the region of Zurich, Switzerland. Three samples were fabricated for each configuration. The wood was initially conditioned at a wet climate of 85% r.h. (relative humidity) and 20˚C. The single layers were cut as strips of 120 mm length and 20 mm width. For the six configurations, the thicknesses were varied according to the ratio r = h 1 (h 1 + h 2 ) −1 as shown in Table 1 where h 1 and h 2 denote thickness of layers 1 and 2 respectively (Fig 2). A growth ring angle φ 2 of zero was chosen for all experimental validation samples. The layers were glued using a one-component polyurethane glue (1cPUR, Purbond HB S709, Henkel & Cie AG, Switzerland) and manufacturer's standards were applied for the gluing procedure. The samples were placed inside a 35% r.h. and 20˚C climate room to dry and curvatures were measured after t = 1, 2, 4, 6, 8, 24, and 48 hours using image analysis (polynomial fits to edgethresholds). Simultaneously, wood moisture contents ω t at time t were measured using a gravimetric determination method, i.e. o t ¼ ðm t À m 0 Þm À 1 0 , where m t is the sample mass at time t and m 0 is the oven-dry sample mass. Finally, experimental data is verified against Eq 5.

Sensitivity analysis
Defining Eq 5 as the model with uncertain input, where each parameter is following a specific probability density function (PDF) with given parameters or moments, model output variability can be assessed using a Monte-Carlo (MC) sampling approach. This procedure allows for calculating partial variances of the model parameters and their combinations w.r.t. the model output variance. We make use of total Sobol' indices to conduct a global-type sensitivity analysis of the model [17,18]. The indices, normalized values between 0 and 1, for each parameter i, equal the contribution of variability of that parameter to the model output variability.
For the analysis, the Matlab-based uncertainty quantification tool UQLab was used [19]. The engineering constants E L , E R , and E T entering the model (for the plane strain case), the three swelling coefficients α L , α R , and α T , the two thicknesses h 1 , and h 2 , and the growth ring orientation φ 2 , were sampled using 10 7 MC samples. We assume a log-normal PDF for the material properties (fE L ; E R ; E T ; a L ; a R ; a T g � LN ), where a coefficient of variation (COV = σ/μ) of 10% is attributed to each property, and a Gaussian PDF for the geometrical properties (fh 1 ; h 2 ; φ 2 g � N ) where standard deviations (σ) of 0.1 mm, 0.1 mm and 1˚are attributed respectively. Exemplary, four beech bilayer configurations were analyzed (described along Results). The mean values μ i (first moments of the input PDFs of i) for the material properties are taken as in Fig 1. The described input parameters were chosen in view of direct applicability in industry. The Young's moduli, the differential swelling coefficients, and the geometry are deemed easier to record than other input parameters (i.e. shear moduli and Poisson ratios).

Model validation with experimental data
The plane strain (i) and the plane stress (ii) model formulations were validated against the obtained experimental results. Fig 3 presents the data of measured (κ d ) and predicted (κ m ) curvature versus wood moisture content (ω). The measured curvatures exhibit significant dependence on the ratio r. A maximum of κ d � 4.5 m −1 was reached for r = 0.17 and a minimum of κ d � 1.5 m −1 was reached for r = 0.49 after 48 hours in 35% r.h. climate, at ω � 9.5%. For the cases (i) and (ii), errors in prediction are displayed starting with a value of zero at initial conditions (ω � 17.5%). It can be seen, that the prediction errors are lower for model (i) than for model (ii). In fact, for low values of r (r = 0.12 and r = 0.17), both models slightly under-predict κ d , while for higher values (r = 0.43 and r = 0.49), κ d is over-predicted. At a value of r = 0.34, which corresponds approximately to h 1 : h 2 = 1: 2, both models predict κ d very well (R 2 i ¼ 0:986 and R 2 ii ¼ 0:992) with an absolute prediction error lower than 10%. Over all configurations, model (i) reached a mean goodness of fit (R 2 ) to experimental data of � This optimal value of r (maximizing κ, e.g. at r = 0.253 for φ 2 = 0˚) is, however, shifted in function of φ 2 . The higher φ 2 is chosen, the higher κ for a same value of h 1 + h 2 . High curvatures of  κ � 10 − 25 m −1 result for 2 mm thick bilayers (depending on φ 2 ). By choosing higher values of h 1 + h 2 , κ drastically drops in magnitude and seems to converge to a minimum value (at a given r) for large h 1 + h 2 according to the spacing of the curves in the semi-log plot for κ. It can be seen, that W tot is a linear function of total thickness h 1 + h 2 for any r (Fig 4b). Considering the contribution of layer 1 (W 1 ) and layer 2 (W 2 ) to W tot separately, it is seen that there are values of r for each layer i where W i is either minimized or maximized. As for κ, these r remain independent of h 1 + h 2 . W 1 is minimized where W 2 is maximized and vice versa, and W tot is minimized at r = 0.273. Surprisingly, the axial interface-stress does not depend on total bilayer thickness h 1 + h 2 but only on r (Fig 4c). Stresses in the passive layer are found to show higher dependence on r as those of the active layer. The minimum stress of the passive layer situates at r = 0.202, whereas for the active layer, stresses seem to remain approximately constant between r = 0.1 and r = 0.4.

Model input parameter sensitivity analysis
Results of the sensitivity analysis are displayed in Fig 5. Four different model configurations were considered for which r was set as the optimal value maximizing κ according to Fig 4 and O: 17.5% ! 9%. In the case of m h 1 ¼ 1 mm and m φ 2 ¼ 0 � , most of the model response variability can be attributed to α R . Input variabilities in thicknesses h 1 and h 2 also contribute to a small amount. In the case of m h 1 ¼ 2 mm, the contribution from the thicknesses is marginal, and all model variability seems to be exclusively attributed to the variability of α R . Setting m φ 2 ¼ 45 � , a major shift from α R to α T in total Sobol' indices can be observed. When m φ 2 ¼ 45 � , additionally, some contribution seems to be attributed to Young's moduli in bilayer axial direction e x of the passive layer (E L ) and the active layer (E R , but surprisingly not E T ). Regardless of the model configuration, no effects of input variability of the parameters φ 2 , E T , and α L are visible.

Model validity
It is seen that the proposed 2D linear elastic model can predict curvature of beam-like bilayer structures with narrow widths and thin layers. Model precision tends to get worse when high or low values of r (close to 1 or 0) are chosen. A plane strain state assumption predicted data better than a plane stress state assumption. This is in line with the known phenomena of strong out-of-plane effects of anisotropic materials when the out-of-plane material axis is stronger than the in-plane axis, as in the case of the active layer in a wooden bilayer. Besides, we note that curvature in such a bilayer is not a uni-axial phenomenon but, will happen about two perpendicular axes simultaneously. A wooden bilayer plate will, however, tend to minimize the change in Gaussian curvature. Curvature out of plane is thus expected to affect axial curvature to some extent. Besides modeling only a reduced geometry (2D), known beech-wood-inherent deformation mechanisms such as visco-elasticity [20,21], mechano-sorption [15], and plasticity [22] are here neglected for the sake of simplicity. Also, the presented model assumes steady-state moisture conditions, which is an ideal assumption, supposedly valid for thin layers. Time and moisture-diffusion effects, especially for increased layer thicknesses, are certainly of further interest when understanding wooden bilayer mechanical behavior. Overall, the relation κ / Δω appears to hold true despite the highly non-linear nature of the model and the complex physical effects in experimental data.

Design principles
Given the validity of the presented model, design aspects can be discussed using parametric studies as presented in Fig 4. For any arbitrary configuration defined by local orientations L, T, and R, material properties, and growth ring angle φ 2 , the optimal thickness ratio r opt , for which the curvature κ is maximized (over any O) can be calculated applying dκ/dr = 0 and solving for r. In Ref [23] an exemplary solution for the case of isotropic bimaterials with constant material properties is given as where E 1 and E 2 are the layer stiffnesses. For the presented model and wood, the optimal ratio r opt is, however, dependent on many input parameters. It was seen that κ can be maximized, or W tot and σ i can be minimized, and that in either case r opt is different. We state that in certain range of r, however, there is a tendency for the maximum curvature being reached while simultaneously, the system minimizes stored elastic energy and axial stresses. Furthermore, model prediction of validation data was better in the case of the samples with r close to r opt , tending to the conclusion that maximized κ and minimized W and σ influence model precision and justify the use of a 2D linear elastic model, as presented, to predict κ for such cases. Given the opportunity to vary h 1 + h 2 or Δω in order to set a target κ, we thus recommend to always set r to an optimal value r opt when designing wooden bilayers.

Increasing model precision
The model parameters of relevance were identified as the swelling coefficients in the active layer in all of the cases studied. Nearly all variability in κ was attributed to the variability of α R (or α T ). Other parameters appear to be of much lower relevance (vanishing total Sobol' indices), and it can be concluded that measuring and updating the model with those parameters is not necessary. By measuring the swelling coefficients and updating the model input to in situ values, most of the epistemic uncertainty in the model response can be reduced to nearly only including the aleatory uncertainty of that input (given the measurement is accurate). Aleatory uncertainty in model input parameters, however, can never be reduced. For wood, as a naturally grown material with high variability, a precision range of predicted model response compared to experimental validation data can realistically not be lower than the range of spread in measured input material parameters.

Remarks on cracking and delamination
Under mechanical, but mainly moisture loading, cross-laminated wooden structures are prone to cracking [24]. For thin wooden bilayers, two basic failure modes can be considered: Axial delamination of the adhesive bond at the layer interface, and transverse cracking in the active layer. Both modes supposedly originate from a combination of the crack separation modes I and II (mode-mixity). At the bilayer edges, the zero axial stress boundary condition results in a stress intensity factor K I at the interface. In combination with the interface-inherent stress intensity factor K II due to the shear transfer, the edge-near interface can be identified as a zone with high risk of failure. We refer to the work of Refs [25,26] where the failure of multilayers is analyzed, and, as general finding, the term h i s 2 i Q À 1 i G À 1 c is to be minimized to reduce failure risk. Γ c denotes the interface toughness and is a configuration-specific parameter to be experimentally determined. The dominant influence of the axial stress σ i on the term can directly be recognized, underlining the importance of choosing a ratio r, for which axial stress is minimized (for beech bilayers with φ 2 = 0, r opt,σ � 0.2). Layer thickness h i , stiffness Q i , and Γ c , only have secondary effects on failure risk. We state that failure was not observed for the experimental bilayer samples in this study, but that it can hypothetically occur for very high Δω in combination with unsuitable design.

Conclusions
We derive and propose a simple 2D linear elastic model for predicting curvature of thin wooden bilayers as a function of moisture change. The model takes into account moisturedependent orthotropic properties for wood and a physically accurate reduction to two dimensions. Experimental curvature data of beech bilayers is accurately modeled within a 10% precision range given an optimal ratio of layer thicknesses. In the design process, anatomical orientations of wood in each layer, the wood species, the growth ring angle in the active layer, the applied moisture difference, and the bilayer thickness combined with the thickness ratio of passive and active layers need to be considered. An optimal thickness ratio can be found for every configuration where either bilayer curvature is maximized or the passive layer axial stress at the interface and the edge is minimized. Axial stresses are found independent of total bilayer thickness and only depend on thickness ratio and applied moisture content difference. Further, they exert a considerable influence on delamination risk. The variability in axial swelling coefficient in the active layer was identified to have a decisive impact on bilayer curvature. Finally, we recommend considering effects such as moisture gradients or time-dependent mechanical behavior for modeling and design of wooden bilayers with increased layer thickness.