Experimental investigation and beam-theory-based analytical model of cross-laminated timber panels buckling behavior

This paper investigates the buckling behavior of three-layered cross-laminated timber (CLT) panels, from both the experimental and analytical standpoints. Two different series of specimens are considered: the homogeneous ones, which are entirely made of beech, and the hybrid ones, whose inner layers are made of Corsican pine. The experimental tests aim to evaluate the failure limit loads of the specimens, when loaded by an increasing compression tip force. The analytical formulation is first obtained for a panel with a generic number of layers and after it is specialized for a three-layered panel. Timber layers are modeled as internally constrained planar Timoshenko beams linked together by adhesive layers, which are modeled as a continuous distribution of normal and tangential elastic springs. A closed-form solution of the buckling problem is obtained. The achieved Eulerian critical load of CLT panels depends on two parameters, which account for (1) the interaction between timber layers (due to the glue tangential stiffness) and (2) the rolling shear stiffness of the inner layer. Three different failure criteria are introduced to estimate the limit load. Finally, the analytical limit loads and the experimental ones are compared.


Introduction
Many types of timber structural elements (i.e., columns, trusses, and frame structures) undergo axial compression or combined axial compression and bending during their lifetime functions; the risk that they deflect laterally, under this load conditions, is called flexural buckling.
The load-bearing capacity of timber members is not easy to predict as it is influenced by several variables which can be gathered in two groups: the first one involves the geometrical features, the boundary conditions, the material and environmental properties (i.e., strength and service class) and the load duration; the second comprises the geometrical (i.e., initial curvature) and structural imperfections (i.e., knots, growth deviations, amount of moisture content). The effect of these latter variables was deeply investigated in Refs. [1][2][3][4][5] while the influence of the material properties on the load-bearing capacity was analyzed by Blaß [6][7][8].
Furthermore, two key aspects should be considered: (i) a geometric effect called P-delta which describes the nonlinear increasing of deformations due to the increasing eccentricity of the axial load, as well as (ii) the non-linear behavior of timber material under compression actions parallel to the grain. A large investigation on the P-delta effect on axial compressed timber elements was carried on by Tetmajer [9] and confirmed by test campaign by Larsen and Pedersen [10], the influence of non-linearity of timber on the interaction between moment and axial force was numerically defined by Buchanan [11,12].
The current National codes and Standards (i.e., Eurocode 5 [13]) provide two methods for the calculation of the loadbearing capacity of axially compressed columns with and without eccentricity; the most simple one is based on the effective length method (ELM), while, the other accounts for a second-order analysis of the structure. This latter method, by the way, considers timber as an elastic material so it is not able to take into account its material non-linearity.
The ELM considers a simple equivalent column hinged at both ends; the P-delta effect is taken implicitly into account through a buckling factor, k c . The buckling factor is used here to reduce the compressive strength of the timber column in the direction parallel to the grain and depends on the effective length of the structural element which can be expressed by the slenderness ratio λ.
Several experimental tests have been performed to define the influence of slenderness ratios on failure modes strain and ultimate bearing capacity of glulam columns made of different timber species: bamboo [14], larch [15], European beech [16]. Moreover, the effect of imperfections such as knots [17] and cracks [18], was also investigated as experimentally as performing FEM fittings and parametric analyses.
A novel strain-based model to analyze the load-bearing capacity with and without eccentricity is proposed in Ref. [19] for solid elements; in this work, it is highlighted how the non-linear behavior of timber when subjected to compression parallel to the grain considerably influences the load-bearing capacity. Based on this latter model and accounting for the Glos' stress-strain relation [20], in Ref. [21], Monte Carlo simulations were performed with the aim of investigating the influence of the varying material properties and slenderness.
Among the aforementioned solid and engineered wood products, an in-depth study of cross-laminated timber (CLT) column components has not yet been carried out. These latter structural elements result from the cut of prefabricated solid slabs made of adjacent planks layers glued together at a right angle.
In the last few decades, CLT has become popular in the mass timber construction field [22], especially for mid-and high-rise buildings, due to the numerous advantages [23] such as the high dimensional stability, elevated in-plane isotropic strength and stiffness with respect to sawn timber, high speed of installation [24], and the environmental benefits compared to other construction materials like steel, concrete, and masonry.
The growing interest on CLT for design purposes and its great mechanical properties has led several authors to analyze deepen its mechanical characteristics as bending, tensile, shear, and compressive strength [25][26][27][28][29][30][31]. Although, as previously highlighted, many studies have been conducted to clarify the in-plane and out-of-plane mechanical behavior of CLT panels, up to now, just limited research has been carried out to invesigate the stability bearing capacity of CLT members under axial compression loads. In Ref. [32], the different axial compression behaviors of cross-laminated timber columns (CLTCs) and control glued-laminated timber columns (GLTCs) were investigated. The GLT column specimens had a greater compressive loading capacity than CLT column specimens but the latter highlighted a better ductility and energy absorption than GLT. In Ref. [33], the main aims were to define the stability bearing capacity of the CLT members in axial compression and to propose a CLT stability coefficient calculation method by suitably modifying that implemented for GLT in the Chinese Standard. In Ref. [34], Ye et al. accounted for the CLT rolling shear sensitivity arising from the out-of-plane bending and compression load combination; an analytical model which takes rolling shear effects into account to predict the loadcarry capacity of CLT compression-bending members was provided.
Stability issues have to be addressed in view of challenging aims such as those analyzed in [35][36][37][38]. In the case of timber engineering, it is to realize increasingly high-rise timber buildings. Because of the low modulus of elasticity parallel to the grain of timber in general and of the even lower shear modulus which affects the ortogonal layers [39], CLT is more prone to undergo instability phenomena. Timber exhibits about one third of the modulus of elasticity of concrete in parallel to the grain, and the radial-tangential shear stiffness is roughly two order of magnitude lower than the longitudinal one [40], and as a consequence, timber structures may suffer from buckling.
Considering the technology of cross-laminated timber, the orthogonal layers are influenced by the radial-tangential shear stiffness, usually called rolling shear stiffness [41], and, therefore, it is expected that shear effects are important in this type of wood-based composite.
Some researches pointed out how the properties of the entire CLT panel are strongly influenced by the type of wood of which its layers are made [42,43], their thickness [44], and the type of glue adopted [45,46]. Other studies [47][48][49][50], have investigated panels made with wood species other than fir and spruce, which are currently the most used species for CLT industrial manufacturing. In this context, panels made of some native hardwoods species such as beech, yellow-poplar, and tropical hardwood (i.e., rubberwood, batai) have shown an excellent behavior. This evidence is relevant considering the impact that the use of hardwood, especially local, would have on the production of CLT panels. The development of a short supply chain of wood is a possible mean to economically enhance the woodlands and to achieve environmental and social benefits [51]. Wood can be effectively used to manufacture products for structural purposes, but it is necessary to investigate the potential of each local species to identify the most suitable technology to employ. An interesting research case was carried out in Italy, where maritime pine from Sardinia was employed to produce crosslaminated timber panels. Preliminary tests were required to evaluate the physical and mechanical properties so as to demonstrate the suitability of this locally grown species to produce CLT panels [52]. In addition, the activation of a short supply chain of timber could lead to a sustainable forest management and to a better exploitation of the raw material [53]. Nowadays in Europe, the attention is focusing on hardwoods especially for manufacturing wood-based composite due to their great mechanical properties [54]. Despite the fact that softwood has played a central rule in the market of glue-laminated timber and cross-laminated timber, the abundance of hardwoods in European forests has determined the need to analyze their performances [55].
The current work aims to deeply analyze the buckling and post-buckling behavior of CLT panels. The specimens'added value lies in (i) the short supply chain origin of both hardwood beech an softwood Corsican pine and in (ii) the homogeneous beech and hybrid beech-Corsican pine configuration, devised for defining a possible optimization of material. To reach the purposed aim, an experimental test campaign is first carried on. As a result, the ultimate load and the failure modes are accounted for. At a second stage, analytical models are formulated and fitted with the previous evidences and then discussed.

Specimens description
To carry on this experimental campaign, two of the most interesting Italian hardwood and softwood timber species, beech and Corsican pine, were chosen. The raw material for both species came from Calabrian forests, South Italy. All the steps that led to obtaining the final product were carried out with short supply chain procedures. An amount of 365 beech boards and 90 Corsican pine boards were specifically classified as D40 and C20 [50] and planed to a thickness of h = 18 mm to guarantee the adequate planar surface necessary for the face gluing. Two different series of three-layered cross-laminated timber (CLT) panels, made, respectively, of beech and Corsican pine and only beech, were realized to carry out an extended experimental campaign to assess their mechanical properties prior to the industrial production. The hybrid configuration is composed by two beech outer layers and one Corsican pine inner layer. The layering of the boards and the gluing process were performed by hand.
Timber boards were tiled, overlapped, and glued together by a melamine urea formaldehyde (MUF) adhesive [Grip-Pro™ Design-AkzoNobel, that was obtained by mixing with the ratio 2-1 the two components: hardener (H002) and liquid flexible resin (A002)]. The narrow board edges were not bonded. The choice of employing the melamine adhesive was due to the usage of beech hardwood boards in both homogeneous and hybrid beam configurations with Corsican pine boards. In particular, the GripPro™ adhesive is specifically designed for softwood as well as for hardwood timber, such as beech, birch, oak, and chestnut. The average of glue grammage was 190 g/m 2 . The adhesive was applied in a factory condition of 17 °C, with an environmental moisture of 50%. Formed beams were then cold-pressed for 2 h, with a pressure of 1.4 N/mm 2 . No finger joints were forecasted in cross-laminated panels of this experimental campaign.
From this production, an amount of seven specimens for each series was accounted for these tests, their features are summarized in Table 1.

Central axial compressive test
The experimental tests were performed at the Laboratory of Structures, Tests and Materials of the University of L'Aquila (Italy), Department of Civil, Architecture and Building and Environmental Engineering. The typical duration of buckling experiments was in the range of 30 min long and an appropriate setup was developed, to provide the desired loading and boundary conditions to the tested CLT panel. Each specimen was positioned between a couple of one-way knife hinges whose shape has been designed to ensure the right positioning in the underlying and overlying wedges (Fig. 1a, c). The distance between the axes of end hinges is l = 1064 mm . The wedges were connected to the machine test support devices by bolts.
During the test, a total of four laser displacement transducers were used. Their position (Fig. 1b) was on both sides of the span area of the specimen; the transducers 1 and 3  were used to measure the middle deflections of the specimen, while the transducers 2 and 4 could measure some torsional effects. The in-plane compressive load F, continuously monitored by means of a load cell, was applied by means of an Instron-Schenck hydraulic servo-controlled jack along the width B.
In the process of specimen loading, a preload F 0 was applied to the specimen to eliminate the gap between the ends of the specimen and the supports. In this experimental test, the preloaded load value F 0 was taken as 30% of the estimated ultimate load value of the specimen.
The tests were carried out under displacement control; in particular, a limit of 15 mm was fixed for the vertical lowering of the machine crossbeam. Two different speed rates have been set: the first, until the critical load was reached, was equal to 0.44 mm/s, and once the critical load was reached, the speed rate has been reduced to 0.33 mm/s.

Experimental results
The ultimate loads F b and the failure modes are summarized in Table 2. Some examples of failure modes are defined in Fig. 2. Regardless of configuration, panels reached similar average limit load levels equal to 200 kN for the homogeneous configuration and 192 kN for the hybrid configuration. Panels having an homogeneous configuration were characterized by a bending due to bucking (tension and compression) failure while the main failure mode of hybrid specimens was rolling shear. The standard deviation s y and the 5th percentile m k are given in Table 2.
The axial displacement versus compressive load F curves are shown in Fig. 3, both for the seven homogeneous specimens ( Fig. 3a) and for the seven hybrid specimens (Fig. 3b). The curves of the homogeneous and hybrid specimens reach the final slope after an initial settlement of the specimens. It is worth observing that after such initial settlement, the slope of the linear increasing branches of the curve is almost the same both for the homogeneous and hybrid panels.
The recorded mid-span displacement versus compressive load F curves are shown in Fig. 4, both for the homogeneous (Fig. 4a) and hybrid (Fig. 4b) specimens. They all show the existence of some problems of the experimental apparatus. Specifically, due to the unavoidable imperfections of the experimental setup, it was necessary to re-elaborate the acquired data to extract the displacement component linked to the deformation of the specimen only. The imperfections of the setup that cannot be fully eliminated include: the misalignment of the two hinges with respect to the longitudinal axis of the specimen, the imperfect contact of the end sections of the specimen with base and top steel plates, and the presence of a gap between the lateral steel plates of the end restraints and the specimen. The accurate procedure to clean up the experimental curves from the imperfections is described in the Appendix. Finally, the re-elaborated midspan displacement versus compressive load F curves are shown in Fig. 4c, d for the homogeneous and hybrid specimens, respectively.

General glued multilayer beam model
To understand the results of the experimental campaign and in the attempt to extend such results for practical use, an analytical model is here developed. A general model is first introduced with the aim to approach the buckling of a glued multilayer beam subject to a concentrated compression load, after which it will be adapted to the three layers CLT panel object of the present paper. As the model is intended to interpret the experimental results described before, it will be made reference to an initial imperfection of the panels as in Fig. 5. It is worth observing that the length l is taken equal to the distance between the two end hinges of the experimental apparatus (l = 1064 mm), even though the panel length is of 800 mm (see Table 1), since it is assumed that the entire system between the hinges undergoes bending-buckling deformation. The general mechanical model is made of N overlapped planar Timoshenko beams, with rectangular cross-sections of different height h i (i = 1…, N) and unique width b; the beams are connected through N − 1 glue lines, represented by a continuous distribution of linear elastic springs that work along the tangential and normal directions, as in Refs. [56,57]. Since the model aims to analyze its buckling behavior, it is developed according to the classical linearized theory approach, which consider the hypothesis of small displacements and small strains, though imposing the equilibrium in a deformed shape adjacent to the reference equilibrium configuration (trivial solution).
The mathematical model is defined by three sets of three equations (kinematic, static, and constitutive) for each layer, which involve kinematic and static quantities related to the glue lines.

Kinematics
With reference to Fig. 6a, b, the kinematics on the entire system is described by a 3N displacement components vector , expressed as a function on a unique variable x, lying along the beam length, in which T encompasses the displacements of the i-th layer. By neglecting the dependence on the variable x for simplicity, a generalized 5N − 2 strain vector follows as = 1 , g1,2 .., gi−1,i , i , gi,i+1 .., N T in which the classical beam strains for the i-th layer are defined, adopting from now on the Lagrangian notation, as  It should be observed that, since the thickness of a real glue line may be of difficult estimation or control, it is worthy to express the components of gi,i+1 independently from the thickness t i,i+1 ; therefore, such components represent displacements, and thus they are dimensionally equal to a length.

Statics
The static quantities of the general model, (represented in Figs. 7, 8 as positive), are the generalized internal forces N i (x), T i (x) and M i (x) (i = 1,..,N), which represent the axial force, the shear force, and the bending moment, respectively, acting on the i-th layer, whereas quantities i,i+1 (x) and i,i+1 (x) are, respectively, the normal and tangential stresses (per unit of length) acting along the contact interface between layers i and i + 1. These can be collected into a (5N − 2) × 1 force vector T with (neglecting again the dependence on the variable x for With reference to Fig. 8, the equilibrium equations, for each layer, should be written with respect to the adjacent deformed configuration defined by the .   (i = 1,..,N), which may include some initial imperfection and, in general, should be different for each layer as the glue lines are allowed to deform out of their plane. Since the forces N i and T i (i = 1,..,N) are not parallel and the sections stepwise, it may be convenient to express the equilibrium condition with reference to an infinitesimal portion of the deformed beam, cut by two vertical sections, spaced dx, and projecting the forces N i and T i in the horizontal and vertical direction. In Fig. 9, the forces involved in the equilibrium of the single layer are portrayed and the apex (⋅) + indicates the small increment derived by a series expansion up to the first order, so that (x) + = (x) + � (x)dx , (x) being one of the kinematic or static entities shown in Fig. 9.
In the same figure, the rotations of the axis line are split in two parts, having distinguished the initial rotation i0 = v � 0i , which may be present in the unloaded imperfect beam (Fig. 5), from the relative rotation i = v � i , which rise from the initial imperfect shape, as an effect of the applied loads. It should be recalled that while the bending response of the layer is related to the (relative) rotation i , the bending effect induced by the horizontal force H + i depends on the (total) With reference to Fig. 9, it should be also observed that strictly the glue stress and vary in an unknown way along the layer's border; however, for simplicity, they are supposed uniform in dx and acting under the rotation it .
For the i-th layer, it can be then demonstrated that the equilibrium equations read: However, recalling that H = ∑ i H i = −F for the reference system (Fig. 5), the equilibrium equations can be simplified supposing that it is small enough for each layer, such ; similarly the terms t and t can be overlooked, thus determining the equilibrium equations in the form:

Constitutive laws
Finally, the constitutive laws are introduced, based on a linear elastic uncoupled behavior for all the materials, according to the following expressions: (as the layers present a rectangular section), I i , E i , and G i are, respectively, the cross-section area, the shear area, the moment of inertia, the Young modulus, and the shear modulus of the i-th layer. For the glue line between layers i and i + 1, the stiffness of the two continuous bed of springs in the normal and tangential direction is defined, respectively ( Fig. 7): having considered the glue as an elastic continuum with Young modulus E gi,i+1 and shear modulus G gi,i+1 . In the next numerical calculations, it is always assumed t i,i+1 = 0.1 mm.
It can be noted that globally, the general mechanical model is defined by 13N − 4 kinematic and static unknowns against 13N − 4 equations; therefore, the algebraic-differential problem associated to the proposed model is balanced. (2)

Three-layer CLT beam model
The base mechanical system described above can be adapted to fit the specific three-layer CLT panels (N = 3) used in the experimental campaign described in the previous sections, resulting in a global simplification of the mathematical formulation. Since the three layers of the specimens have the same thickness h and the same width b, a unique cross-section area A i = A = bh , shear (rectangular) area A it = A t = (5∕6) bh , and moment of inertia I i = I = bh 3 ∕12 is considered for each layer ( ∀i = 1, 2, 3).
Due to the mechanical symmetry of the section of the CLT panels, as the outer layers are made of the same kind of wood, a single value for the elastic modulus E 1 = E 3 and for the shear modulus G 1 = G 3 are taken into account for the external layers.
Finally, for both the two glue lines, the same thickness t, Young modulus E g , and shear modulus G g are supposed; therefore, a single normal stiffness k = k 1,2 = k 2,3 as well as a single shear stiffness g = g 1,2 = g 2,3 are considered.
Equilibrium Eq. (3) now became: The CLT beam model can be further simplified by introducing some internal constraints, based on suitable assumptions on the mechanical behavior of the system. These constraints are expressed as follows: By means of Eq. (6) 1,2 , the impossibility of the layers to detach from each other in the normal direction is introduced, thus resulting in a single displacement component v 1 = v 2 = v 3 = v orthogonal to the layers' axes. Likewise it is supposed that the initial imperfection v i0 is equal for each layer, thus assuming v i0 = v 0 ∀i = 1, 2, 3 (Fig. 5). Based on the low thickness to length ratio of each layer, Eq. (6) 3,4 enforces the shear non-deformability of the external timber layers (Euler-Bernoulli beam theory); such a constraint makes the rotation angles 1 and 3 slave quantities of the orthogonal displacement v, thus assuming the same value . Finally, Eq. (6) 5 is based on the fact that often in CLT panels, as in the case of the ones tested experimentally, the cross layers are made of planks whose narrow borders are not glued or not even in contact; it is then supposed that the central layer is unable to bear any axial force. Nevertheless, as the inner planks hold some cross bending stiffness, which can be involved by the deflection of the outer layers, it is supposed that the central layer is able to bear some bending moment (M 2 ≠ 0). From Eqs. (6) 5 and (5) 2 , it can now be written 1,2 = 2,3 = , which means that the two glue lines are subject to the same shear stress and the same shear strain ∕g = t = t1,2 = t2,3 while the global compressive force acting on the panels (Fig. 5) reads: N = N 1 + N 3 = −F keeping itself constant along the panel.
Due to the internal constraints, some internal forces ( T 1 , T 3 , 1,2 , 2,3 ) become reactive quantities, as they cannot be expressed by the constitutive laws. Such reactive quantities can be then obtained by solving Eq. (5) 7,9 with respect to T 1 and T 3 , and Eq. (5) 4,6 with respect to 1,2 and 2,3 . Together with Eq. (5) 1,3 , the condensed equilibrium equations are, hence, obtained by substituting the relationships that provide the reactive forces 1,2 and 2,3 into Eq. (5) 5 , and finally substituting, in the same equation, the derivatives of T 1 , T 2 , and T 3 expressed from Eq. (5) 4,5,6 . The condensed equilibrium equations, thus, read: in which it has been considered, from Eq. (6) 3,4 , that . By substituting the constitutive laws (Sect. 3.1.3), and the kinematics equations (Sect. 3.1.1) (modified according to Eq. (6)) into Eq. (7), the field equations are obtained as: which represent a system of four coupled linear differential equations facing five kinematic unknowns ( u 1 , u 3 , t ,v, 2 ); therefore, since the system is indeterminate, a new constraint is needed. However, before doing this, further considerations have to be made.
It should be now taken into account the relationships that exist among the various displacement components (Fig. 10), according to Eq. (1), from which the two axial displacements for the external layers are obtained in the form: Summing both sides of Eq. (9) and deriving it, it follows: 1 + 3 = 2 2 which means that the axial displacement and strain of the central layer identify, respectively, the mean displacement and strain of the outer layers. Therefore, it can be written −F = E 1 A 1 + E 1 A 3 = E 1 A( 1 + 3 ) = 2E 1 A 2 and then 2 = −F∕ (2E 1 A) which implies that the central layer axial strain is constant along the panel. This means that deriving twice Eq. (9) and substituting them into the field equations, Eq. (8) 1,2 become equals, thus resulting in a differential system twice indeterminate (three independent equations against five unknowns). Therefore, two more constraints are needed.

Coupling parameters
To eliminate the indetermination of the problem, a simple but consistent attempt, from a mechanic standpoint, is to introduce the following two non-dimensional parameters: Although as well as should be, generally, functions of x, it is now supposed that they are both constant showing later the conditions under which such a hypothesis turns to be true. Then introducing Eq. (10) into Eq. (9), after derivation, the outer layers axial strains become: From Eqs. (10) and (11) 1 , the field Eqs. (8), which have now become an algebraic-differential system of three equations facing three unknowns ( v , , ), read: in which the elastic modulus ratio E = E 2 ∕ E 1 has been introduced. Equation (12) 2 can be then written in a compact form: which, in absence of an initial imperfection (i.e..v 0 = 0 ), clearly recalls the classic elastic curve equation for the buckling of a Euler-Bernoulli's beam, having introduced the equivalent moment of inertia: The two parameters and can now be determined from Eq. (12) 1 , 3 , thus obtaining: These last equations effectively express two constant quantities provided that the differential Eq. (13) admits a solution v in the form v(x) = − e sin ( x∕ l) . As known, such a circumstance may stem out assuming v 0 (x) = − e 0 sin ( x∕ l) for the imperfect beam problem of Fig. 5, which, from now on, will be considered as the reference system for the analytical model, as it can be directly associated to the experimental set up.
From the abovementioned sine solution, it can be determined: which, according to Eq. (14), determines the degree of coupling of the external layers. In this context, the parameter represents the coupling flexibility of the glue lines whereas the parameter represents the coupling efficiency of the central layer. Among all the parameters that appear in Eq. (16), those, which a specific attention is given to in this paper, are the glue line shear stiffness g and the shear modulus G 2 , the latter representing the rolling shear modulus of the central layer. With reference to Fig. 11, it is possible to define some limit conditions associated to the extreme values that g and G 2 may achieve.
No slip in the glue lines, with the global section remaining plane while the quantity Ah 2 (1 + − ) (Steiner's term) reaches its maximum value as 1 + − = 2.

Coupled beam
No slip in the glue lines, with the global section that exhibits a broken line profile.

Uncoupled beam
With no tangential stress in the glue line, with a stepwise global section and with the Steiner's term that vanishes, the uncoupled panel behaves as if it was made of three stacked independent beams.
On the basis of the geometrical dimensions exposed in Table 1 and the mechanical properties shown in Table 3 (Sect. 4) for beech wood and glue, the variation of the coupling parameters and with g and G 2 is portrayed in Fig. 12. As can be observed especially from the right graphs of the figure, whatever value of the rolling shear stiffness G 2 , both the coupling parameters and depend very little on the stiffness of the glue g, since over a very small value of g > 10,000 N/mm 2 , the curve is substantially horizontal.

Euler's buckling load and the elastic solution for the imperfect three-layer CLT beam
Recalling Eq. (13), the Euler's critical load for the threelayer CLT panels follows the classical form: over which, according to Eq. (14), it is possible to evaluate the influence of g and G 2 . With reference to the geometrical (Table 1) and mechanical (  Fig. 13, the variation of F cr with g   and G 2 is portrayed, by means of a buckling load efficiency F cr = F cr ∕F cr,FC , being F cr,FC the critical load for the full composite beam.
Clearly, the influence of glue stiffness g is negligible whereas the rolling shear modulus plays a crucial role.
From the sine initial imperfection v 0 (x) , the solution of Eq. (13), assuming the following boundary conditions v(0) = v(l) = v �� (0) = v �� (l) = 0 (Fig. 5), reads: in which the applied force F and the initial eccentricity e 0 are assumed known. Specifically, the initial imperfection e 0 can be estimated by comparing the experimental midspan displacement-compressive load F curves of Fig. 4c, d and the correspondent analytical ( v(l∕2) versus F ) curves obtained from the solution of Eq. (21). Figure 14 shows the experimental mid-span displacement-compressive load F of the homogeneous specimens HO-1, HO-4, HO-5, and HO-7, already shown in Fig. 4c (thin lines), and the analytical curves (thick lines) obtained for different initial imperfection e 0 . As can be observed, the best fitting among experimental and analytical curves occur for e 0 = 0.25 mm (≈ L/4000). As a final clarification, the curves in Fig. 14 are shown in a zoomed window with respect to those shown in Fig. 4c. Moreover, Fig. 14 shows the absolute value of the mid-span displacement. In the next section, the value of e 0 = 0.25 mm will be always used.
From Eq. (21), it is now possible to evaluate all the displacements, strains, and forces of the system. In particular, focusing on the static quantities, introducing the single i-th (i = 1, 2, 3) layer's critical load as F cri = E i I 2 l 2 (F cr1 = F cr3 ) and recalling the kinematic and constitutive equations of Sect. 0, it can be written from Eq. (11) 1 : from Eq. (10) 2 : from Eq. (5) 7,9 : e cos l x from Eq. (10) 1 : eventually from Eq. (5) 4,6 : On the basis of the geometrical dimensions exposed in Table 1 and the mechanical properties shown in Table 3 (Sect. 4) for beech wood and glue, the elastic response of the homogeneous three-layer CLT panel is then described by the functions portrayed in Fig. 15.

Stability verification
Recalling the theoretical nature of the Euler's buckling load and the fact that, in principle, it could never be reached, at least within the boundaries of the linearized theory, to interpret the results of the experimental campaign, it is now necessary to define some specific failure criteria, as in Ref. [58].
Basically, as it was described in Sect. 2.3, three failure modes have been observed in the tests: a bending-buckling mode, occurring in the middle of the panel, a rolling shear mode, detected at the ends, and a delamination mode, involving the glue at the extremities as well. To each of these modes, it is possible to associate an ultimate compression load which account for the achievement of the relevant failure stress in a specific part of the structural element. The ultimate load of the panel will be then the lowest of the three, for which an assessment is here presented.
Bending buckling mode The flexural failure mode stems from the minimum (absolute maximum) compression stress reached in the middle of the panel at the very external surface of layer 1. Recalling the classical Navier's formula, such a stress can be calculated as: and therefore, from Eqs. (22) and (23) Given f cu the compression failure stress of the layer's wood, it determines a limit compression load, for the straight panel, equal to F cu = 2 A f cu . The search for the ultimate load F = F cb , which makes x,min = f cu according to Eq. (31), follows the Ayrton-Perry criterion. After introducing the following non-dimensional parameter: from Eq. (31) it can be written: It is now introduced the relative slenderness in the form:  where F rb represents the limit compression load which determines the achievement of the rolling shear strength in the inner layer. This is, hence, defined as:

Delamination mode
The delamination of the glue at the ends of the panel should likely derive from the combined action of the shear stress and the normal traction (Fig. 15); however,in this paper, for the sake of simplicity, it will be considered solely the first one.
The maximum shear stress (per unit of length) is evaluated from Eq. (25) for x = l , thus obtaining: Given the failure shear stress of the adhesive, which determines the delamination, its associated value u = u b stems from an ultimate load F = F gb which makes u = max according to Eq. (43). Therefore, defining: it can be written: and finally: Finally, the limit compression load for the CLT panel reads:

Comparison between analytical and experimental results
The numerical value of the mechanical properties of the wood species considered in this analysis is reported in Table 3. The mean elastic moduli E 1 and E 2 were identified in Ref. [59], where the same CLT panels were mechanically characterized from bending tests. For beech wood, the rolling shear modulus and strength have been taken from [60], whereas the compression strength has been approximately defined from EN338, considering the mean experimental MOR shown in Ref. [59].
For Corsican pine wood, the rolling shear modulus has been calculated from EN338, as a fraction of 10% of the shear modulus normal to grain for a C20 graded soft wood, while the rolling shear strength derives from [50].
On the basis of a shear modulus of the glue G g = 642 N/mm 2 , [59] and of an alleged glue line thickness t = 0.1 mm, it can be calculated from Eq. (4), the shear stiffness of the glue: g = 924.480 N/mm 2 . However, regarding this last value, it should be considered the uncertainty related to the thickness of the glue line; therefore, it seems reasonable to consider g in the range 10 5 ÷ 10 6 N/mm 2 .
The adhesive shear strength has been taken from [61].
The limit failure loads of the CLT panel are reported in Table 4. Specifically, the theoretical critical load for both the homogeneous and hybrid specimens is shown in the first line of Table 4. Such a critical load F cr is obtained from Eq. (20). The limit loads that refer to the three different failure modes of the CLT panel are reported in the second, third, and fourth rows of Table 4. They are obtained from Eqs. (37), (42), and (46), respectively. Finally, the limit compression load F b for the CLT panel, identified with the minimum value of the previous three limit loads, is reported in the last row of Table 4.
A comparison among the experimental and the analytical limit loads, here obtained, is finally performed. The analytical limit loads of the homogeneous and hybrid configuration of the CLT panels are displayed in the first row of Table 5. The mean experimental limit loads are reported in the second row of the same table. As can be observed, the percentage differences among analytical and experimental values are quite small. Specifically, the analytical model seems to evaluate better the limit failure behavior of the homogeneous CLT panels.
The minor ability of the analytical model to describe the failure limit behavior of the hybrid CLT panels seems to be related to the statistically scattered values provided by the experimental test. From a deeper observation of the numerical data in Table 2, the limit loads of the specimens HB-1 and HB-4 are higher than the highest limit load provided from the homogeneous CLT panels, which are supposed to be of higher resistance as the inner layer holds better mechanical properties. Since such a result appeared conceptually implausible, and thus to be traced back exclusively to the material uncertainty, it was deemed of interest to consider, beside the full range of test results for hybrid panels, also a reduced range in which the most unlikely results were excluded. This fact suggests to remove the results of the specimens HB-1 and HB-4 from the calculus of the experimental mean value of the limit load. The updated experimental limit load of the hybrid CLT panels, obtained from the HB-2, HB-3, HB-5, HB-6, and HB-7 specimens, is reported in the second row of Table 6. Considering this updated value, the relative difference between the analytical and experimental values decreases and becomes of the same order of magnitude of that of the homogeneous CLT panels.
Clearly, the comparison presented before relies on the specific value of the initial imperfection calibrated on the basis of the experimental results (Fig. 14). However, for designing purposes, the practical application of the analytical model would require the adoption of a general value of the initial imperfection e 0 , consistent with the standard production's tolerance. Obviously higher imperfections would result in lower limit compression loads for the CLT panel. For example, applying a value e 0 = 1.00 mm (≈ L/1000) to the reference system described before, it can be calculated a limit compression load F b = 161 kN for HO specimens and F b = 140 kN for HB specimens, with a relative reduction, with respect to the use of e 0 = 0.25 mm, respectively, of almost 13% and 11%.

Conclusion
In this paper, the buckling behavior of three-layered CLT panels was investigated, first by means of 14 experimental tests and then by developing an analytical formulation calibrated on the specific test setup.
The experimental campaign was conducted on two different groups of seven specimens, each one related to a specific configuration characterized by the panels crosssection's arrangement. The first configuration (homogeneous, labeled HO) presented three layers made of the same timber specie (beech), whereas, in the second one (hybrid, labeled HB), the inner layer was made of a different timber specie (Corsican pine). The experimental tests aimed to evaluate the failure limit loads of the specimens, under an increasing tip compression force. The homogeneous panels showed almost exclusively a mid-span bending failure mode while for the hybrid panels, different failure mechanisms, such as rolling shear failure, as well as delamination failure, were observed at the ends, singularly or in conjunction. Such a circumstance is likely to follow from the lower rolling shear strength of the inner layer in the hybrid panels. Nevertheless, the two groups showed quite similar average limit loads; however, the hybrid specimens presented results with a higher standard deviation with respect to those related to the homogeneous ones. An analytical model was then developed in the attempt to describe theoretically the system object of the experimentation. The moderately high length to width ratio of the specimens, the absence of side restraints in the test setup and, above all, the intent to develop a simplified approach to the problem at hand, led straight to a beam-based theoretical formulation. This was first introduced for a panel with a generic number of layers, modeled as a stack of planar Timoshenko beams connected each other by continuous distributions of normal and tangential elastic springs representing the glue lines. The general model was then adapted to a three-layer system, similar to the panels of the experimental campaign, by means of some internal constraints, and considering a side displacement field in the form of a sine, as a result of an initial sine shape for the imperfect panel's longitudinal axis line. From such a model, a closed-fform solution for the elastic problem, together with the related Eulerian critical load, was determined; this was performed by introducing two non-dimensional parameters, and which express the interaction among the wooden layers due, respectively, to the tangential stiffness of the glue and the rolling shear tangential stiffness of the inner layer. It should be noted that the abovementioned closed-form solution stems from the hypothesis on the initial imperfect axis line's sine shape, which introduces some strict similarities with the -process [62]. Nevertheless, the proposed model includes the influence of the elastic response of the glue lines, which represents the novelty of the presented theoretical approach with respect to others [62] which generally overlook the presence of the adhesive layers. However, although it was recognized that the influence of the glue stiffness on the critical load is negligible, it should be highlighted that the modeling of the adhesive layers, beyond allowing a richer description of the response of the three-layer CLT panel, permitted the investigation of the delamination failure mode observed during the laboratory tests. Such a delamination mode enlarges the range of possible failure modes for a CLT panel, which includes the classical bending mode and the shear mode [58].
On the basis of the achieved closed-form solution and recalling the different failure mechanisms [63][64][65] detected during the experimental campaign, three different failure criteria were introduced, leading to a unique theoretic compression limit load. The first one accounts for the bending-buckling failure, the second refers to the rolling shear failure mechanism of the inner layer, and finally the last takes into account the delamination phenomena.
A comparison among the experimental and the theoretic limit loads was performed after the calibration of the analytical model on the basis of some of the experimental results as well as on some parameters derived from the literature. The comparison showed that the analytical formulation describes sufficiently well the real behavior of the CLT panels, both with a homogeneous or hybrid configuration, with relative differences on the limit loads lower than 20%.
As a beam-based formulation, the proposed model, besides being adopted for a comparison with the specific test described above, seems to configure a potential preliminary design tool for stability verification, since it considers one of the most demanding configuration for a three-layer panel wall, though substantially unusual in the construction practice, with no vertical sides' restraints and hinges on both the horizontal ends. However, it should be recognized that the model is clearly unable to describe the two-dimensional buckling behavior of a panel wall in the presence of restraints along its vertical border and/or of a width to height ratio higher than unity. displacement, the displacement component linked to the deformation of the specimen alone can be rewritten: Specifically, Fig. 16, in abscissa, shows the mid-span displacement evaluated by Eq. (48).
Funding Open access funding provided by Università degli Studi dell'Aquila within the CRUI-CARE Agreement. This study was funded by Italian Ministry of University and Research-PRIN 2015, Grant number 2015YW8JWA_002.

Conflict of interest
The authors declare no conflicts of interest.
Ethical approval All procedures performed in studies do not involve human/animal participants.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.