Refined Zigzag Theory: an appropriate tool for the analysis of CLT-plates and other shear-elastic timber structures

Cross laminated timber (CLT), as a structural plate-like timber product, has been established as a load bearing product for walls, floor and roof elements. In a bending situation due to the transverse shear flexibility of the crossing layers, the warping of the cross section follows a zigzag pattern which should be considered in the calculation model. The Refined Zigzag Theory (RZT) can fulfill this requirement in a very simple and efficient way. The RZT, founded in 2007 by A. Tessler (NASA Langley Research Center), M. Di Sciuva and M. Gherlone (Politecnico Torino) is a very robust and accurate analysis tool, which can handle the typical zigag warping of the cross section by introducing only one additional kinematic degree of freedom in case of plane beams and two more in case of biaxial bending of plates. Thus, the RZT-kinematics is able to reflect the specific and local stress behaviour near concentrated loads in combination with a warping constraint, while most other theories do not. A comparison is made with different methods of calculation, as the modified Gamma-method, the Shear Analogy method (SA) and the First Order Shear Deformation Theory (FSDT). For a test example of a two-span continuous beam, an error estimation concerning the maximum bending stress is presented depending on the slenderness L/h and the width of contact area at the intermediate support. A stability investigation shows that FSDT provides sufficiently accurate results if the ratio of bending and shear stiffness is in a range as stated in the test example. It is shown that by a simple modification in the determination of the zigzag function, the scope can be extended to beams with arbitrary non-rectangular cross section. This generalization step considerably improves the possibilities for the application of RZT. Furthermore, beam structures with interlayer slip can easily be treated. So the RZT is very well suited to analyze all kinds, of shear-elastic structural element like CLT-plate, timber-concrete composite structure or doweled beam in an accurate and unified way.


Introduction
Cross laminated timber (CLT), as a structural plate-like timber product, has been established as a load bearing product for walls, floor and roof elements and is in operation since about 30 years. Due to its high resistance, stiffness and adaptability to in-plane and out-of-plane load situations it is becoming increasingly popular among architects and engineers. Despite of its important role and relevance in the timber construction sector it is not included in Eurocode 5, the European design standard for timber structures apart from National Application Documents. The current research activities are summarized in the state-of-the-art reports by Brandner et al. (2016Brandner et al. ( , 2018. For the analysis of other types of shear-elastic composite members, several analytical and numerical methods can be found in the literature. According to its origin, the main focus of application of the RZT has so far mainly been to the engineering divisions of aviation and shipping as well as the automotive industry, where different composite materials have to be tailored to optimize the structural and thermal behaviour of such components.
When wooden parts are assembled into composite structural elements, their connections are usually characterized by a relatively high degree of deformability arising from shear strains. With CLT, this results from the orthotropic properties and the arrangement of the wooden boards. In other cases, this is due to the flexibility of discrete shear connectors. It is the explicit aim of this paper to show that RZT is a suitable and reliable means to analyze such structural elements in a uniform manner, which is demonstrated by examples given in Sect. 3 and 4.

Existing analysis methods
In the following, mechanics of cross laminated timber plates under uniaxial bending caused by out-of-plane loads only is discussed. When dealing with CLT-plates, the transverse shear flexibility of the layers has to be considered. The large difference of the shear modulus of the neighbouring board layers leads to a non-planar warping of the cross section. To take into account this important fact, different kinematic models are developed. Well-known procedures are the Gamma-method (EN 1995(EN -1-1 2015, originally used for composite beams with interlayer slip by Möhler (1956), which can be modified to apply to CLT and the Shear Analogy (SA) method introduced by Kreuzinger (1999) and worked out by Scholz (2002Scholz ( , 2004. For beams, the First Order Shear Deformation Theory (FSDT) by Reddy (2004), and Timoshenko's beam theory can be used when an appropriate shear correction factor is installed. The Gammamethod and the FSDT are not able to tackle the characteristic normal stress peaks that occur near intermediate support locations, nor the variability of shear stress pattern along such a beam. For the SA model, Scholz (2004) gave a series of analytical solutions for a single-span beam. For more complex beams, a numerical model with growing effort has to be established. As reported (Thiel 2013), the SA method overestimates the shear stress peaks close to point loads and inner supports, while the normal stress is quite accurate. In addition to the above mentioned simple beam models, there exist a large group of single layer (Abrate and Di Sciuva 2017) and layer-wise theories (Reddy 2004;Guggenberger and Moosbrugger 2006). The latter can be characterized as very accurate but very elaborate, because the number of kinematic variables is dependent on the number of layers.

The Refined Zigzag Theory
Among a further group of theories, named as Zigzag theories (Carrera 2003) the so-called Refined Zigzag Theory (RZT), developed by Tessler et al. (2009Tessler et al. ( , 2010 and Tessler (2015) has proven to be one of the most promising approaches. This is due to its simple kinematics and superior capacity for the prediction of static, dynamic and buckling behaviour. The RZT takes the FSDT-kinematics as a base line and enhances it by a zigzag rotation (x) , which controls a warping function (z) that can be estimated in advance from the sequence of layers and its material data. The kinematic relations for a beam in the x-z-plane are given in the fundamental paper (Tessler et al. 2009), where the axial displacements u (k) in x-direction and the transversal deflection w (k) in z-direction of the k-th layer were defined (Fig. 1).
The zigzag function is defined in terms of its layer-interface values (i) ( i = 0, 1, … , N ). Using a local dimensionless coordinate (k) ∈ [−1, 1] , the linear function within each layer is defined by Fig. 1 Laminate of three layers with coordinate system and additive composition of axial displacements Using the comma notation to denote partial differentiation, the corresponding strains from linear elasticity theory yield: The function (k) is the first derivative of the zigzag function (k) after the variable z and is layer-wise constant (see Fig. 2). Thus, the shear strain and shear stress are layerwise constant. It represents the mean value and violates the demand of continuity of shear stresses beyond the layer interfaces. The key factor of RZT is the determination of the (k) -function in advance depending only on the stacking sequence and the transverse shear stiffness of each layer, represented by the shear modulus G (k) xz . In the fundamental paper by Tessler et al. (2009), the given procedure is related to rectangular cross sections, which means that all layers have the same width b (k) , as shown in Fig. 2. Therefore, for a more general representation, the calculations are modified by considering a surrogate rectangular cross section (terms noted by a bar, see Fig. 3). With the widely used approximate formula of Jourawski (see Gere 2001) for the shear flow t (k) in each layer, the mean shear strains of the two cross sections are equated, and hence, the transformed k-th shear modulus Ḡ (k) xz is obtained (Fig. 3).
With these values, the calculation for the overall modulus Ḡ of the whole laminate and the modified terms for the (k) -function are performed. For details see Tessler (2015).
And finally the modified (k) -function reads With the assumptions that each layer is linearly elastic and orthotropic corresponding to the Cartesian coordinates (x, z), the constitutive relations take the form For a plate in plain strain in the width direction xz are the Young's modulus and shear modulus of the k-th layer, while (k) xy , (k) yx denote the major and minor Poisson's ratios (Reddy 2004). The resultant stress vectors, ̃ x and ̃ xz , and the generalized constitutive matrix are obtained with Eqs. (4) and (5) by integration over the cross section area (each layer has thickness h (k) and can exhibit a width b (k) ).
Combining Eqs. (14a), (14b) leads to The variation of strain energy can be expressed in terms of generalized stresses and strains It contains derivatives of the kinematic quantities of first order only. Consequently, for a finite element model a C 0 -approximation is sufficient. For further details, the reader is referred to Wimmer and Gherlone (2017). The classical RZT is verified by comparison with exact analytical solutions given by Pagano (1969) and high fidelity 2D-FE calculations. Such assessments are performed by Iurlaro et al. (2013Iurlaro et al. ( , 2016). An extended version of RZT, called RZT (m) is derived from Reissner's Mixed Variational Theorem by Tessler (2015). It leads to very accurate shear stress distributions that satisfy the equilibrium conditions along all layer interfaces and at the outer surfaces as well.

Implementations of RZT-beams
For a finite element approximation only C 0 continuous kinematic approximations are required. A first element was developed by Gherlone et al. (2011), who have introduced an anisoparametric shape function for transverse deflection w(x), which guarantees that no shear locking will occur. Oñate et al. (2012) presented an isoparametric element with selective reduced integration to prevent shear locking. Wimmer and Gherlone (2017) gave explicit expressions for the linear stiffness, the geometric stiffness, mass matrices as well as for the consistent load vector based on the anisoparametric approach. Nallim et al. (2017) presented a hierarchical p-type finite element, while Iurlaro et al. (2015) developed a beam element based on a higher order zigzag theory, where thickness stretching is also taken into account.
Recently, applications of the mixed version of RZT are given by Groh and Tessler (2017) and Kefal et al. (2018). Wimmer and Nachbagauer (2019) treated beams on two-parameter elastic foundation. Additional aspects concerning the zigzag function are presented by Gherlone (2013) and Flores et al. (2018). A different numerical approach is presented by Wimmer and Nachbagauer (2018), where the governing first order differential equation system of the mechanical problem is formulated and solved by the transfer matrix method. For the static case, an exact expression of the transfer matrix is developed by similarity transformation, which provides analytically exact results for the kinematic variables and their first derivatives as well. This approach has the advantage that no error-prone smoothing process is involved to get the beam strains and stresses. Moreover, the transfer matrix allows for the derivation of a further finite beam element, which provides exact deformation results with coarse discretization. For the dynamic case, a series solution can be invoked to get the transfer matrix.

Composite timber structures with interlayer slip
The RZT is well suited for the analysis of composite beams, where special shear connectors (shear studs, nails etc.) are responsible for the interaction of the subcomponents. The behaviour of such members is quite complex and requires the consideration of the compliance of the connection, which leads to an interlayer slip between the different members.
There are numerous proposals for the analysis of mostly twolayer composite members based on analytical or numerical models (Girhammar and Gopu 1993;Pan 1993, 2007;Adam et al. 1997 Herein c s denotes the shear stiffness and b (i) the width of interfacing layer. Otherwise, we get the shear slip by Hooke's law (Fig. 4) Equating (17) and (18) leads to the fictitious shear modulus G (i) xz of the interfacing layer With the help of the proposed model, it is also possible to simulate the effect of thin adhesive layers in any kind of sandwich and composite beam (Huang 2003).

Examples
As a first reference example, a double span beam "T2" with equal length of L = 4.80 m is chosen (Bogensperger et al. 2012). The cross section consists of five layers with 32 mm thickness each, the width is fixed with 1000 mm. Material data are as follows: Material M1: E 0 = 11600 N/mm 2 , G 0 = 720 N/mm 2 , Material M2: E 90 = 0 N/mm 2 , G 90 = G R =72 N/mm 2 . Layer-Sequence: M1-M2-M1-M2-M1. Distributed load p z = 5 kN/m 2 . For details of the 2D-FEM models see Bogensperger et al. (2012). Table 1 shows the differences of various calculation methods. The result of RZT is obtained via 200 anisoparametric 2-noded beam elements in each span. Further comparisons were made based on a more realistic support condition, whereas the reaction force is uniformly distributed over a length a = 192 mm. Figures 5, 6, 7 and 8 show the distribution of normal stress, mean shear 1 3 stresses and first derivative ′ of the zigzag rotation. With the latter we are able to localize the specific sector, where the normal stresses deviate from the FSDT type. Using the latter, the stress peak at the inner support is considerably underestimated. Figure 9 reports the deviation of FSDT-based results from the more accurate values of RZT concerning the maximum bending stress at the intermediate support for two cases of contact. The difference decreases with rising slenderness L/h and contact width a of the reaction force. A noticeable difference between the two approaches is shown in Fig. 7, where the FSDT considerably underestimates the mean shear stress in the weak layer (rolling shear). As a second example a stability investigation is made with a comparison between FSDT and RZT. The same beam as before was chosen but with a span L = 3.00 m. The initial external normal force is set to F x = 100 kN at the right end. In Table 2, the first three eigenvalues of buckling analysis are shown for pinned supports and for clamped support on the left end. In Fig. 10, the corresponding buckling modes are presented.

Fig. 4 Interlayer slip and shear gap
It should be noticed that the current approach is able to distinguish additional types of boundary conditions in an The third example, related to a non-rectangular cross section, is the widely referenced composite beam of Girhammar and Gopu (1993) and Girhammar and Pan (1993). This member is analyzed by introducing a fictitious interface layer of 1mm thickness (Fig. 11). Following Eq.     Table 3 shows the excellent convergence of the anisoparametric finite elements (Wimmer and Gherlone 2017) and a comparison with existing solutions for a wide range of slenderness.
For the dynamic case, refer to Xu and Wu (2007b) who gave analytically exact solutions for a composite Timoshenko beam including the effect of shear deformation and rotary inertia. The comparison reflects three different boundary conditions (S...simply supported, C...clamped, F... free). In the clamped situations, all four degrees of freedom are set to zero. Table 4 shows very good agreement of RZT values with the solutions given in the above mentioned paper. It should be noticed, that the referenced solution is dependent on the choice of the shear correctors for which different approaches exist. The fourth example deals with a composite timber beam of non-rectangular cross section connected with discrete shear dowels (Kneidl and Hartmann 1995), see Fig. 12. The material data are given by .90 ⋅ 10 2 MPa. As in the example before, the problem is solved by Fig. 11 Composite beam with two materials and partial interaction b (2) = 75⋅0.001 0.12 = 0.625 MPa. Kneidl and Hartmann (1995) solved this problem by a two-part Bernoulli beam model, each part lying in the center-line of the partial cross section and connected by a truss ensemble which generates the shear flexibility of the connection. Figure 13 reports the distribution of the shear flow along the beam. In Fig. 14, representative cross sectional stress distributions are shown. The values in brackets and the dotted lines originate from the paper of Kneidl and Hartmann (1995). The stresses in the cantilever part form a self-equilibrating system, which is recorded very accurately.

Conclusion
This paper sets out an alternative, reliable and robust approach for the analysis of cross laminated timber elements under uniaxial bending in statics, stability and dynamics. The Refined Zigzag Theory, assessed as a competitive and accurate model for laminated composite structures is applied to CLT-beams. Comparisons with available beam solutions were made and show very good results. Hence, RZT-based finite elements, as implemented in the author's program ZZBEAM, can serve as an appropriate way to analyse the mechanical behaviour of CLT-structures. The layerwise distortion of the laminate's cross section is captured while the effort is independent of the number of layers. The method does not use any kind of shear correction factor and shows no locking behaviour. The RZT-kinematics is able to reflect the specific stress situation near concentrated loads in combination with a warping constraint as it can be usually detected only by 2D finite element calculations with high computational cost. The order of magnitude between FSDT-and RZT-solutions concerning the maximum axial bending stress are presented for a two-span CLT-beam example. Stability investigations suggest that FSDT-based models are sufficiently accurate if the ratio of bending and shear stiffness is in a range as stated in the above example. By some modifications in the calculation of the zigzag function, non-rectangular cross sections can also be analyzed. Combined with the introduction of fictitious layers to simulate shear-compatible connectors, the method is applicable to a variety of composite wooden structural members, as demonstrated by a number of examples under static and dynamic loads.