On displacement-based and mixed-variational equivalent single layer theories for modelling highly heterogeneous laminated beams

Abstract The flexural response of laminated composite and sandwich beams is analysed using the notion of modelling the transverse shear mechanics with an analogous mechanical system of springs in series combined with a system of springs in parallel. In this manner a zig-zag function is derived that accounts for the geometric and constitutive heterogeneity of the multi-layered beam similar to the zig-zag function in the refined zig-zag theory (RZT) developed by Tessler et al. (2007). Based on this insight a new equivalent single layer formulation is developed using the principle of virtual displacements. The theory overcomes the problem in the RZT framework of modelling laminates with Externally Weak Layers but is restricted to laminates with zero B-matrix terms. Second, the RZT zig-zag function is implemented in a third-order theory based on the Hellinger–Reissner mixed variational framework. The advantage of the Hellinger–Reissner formulation is that both in-plane and transverse stress fields are captured to within 1% of Pagano’s 3D elasticity solution without the need for additional stress recovery steps, even for highly heterogeneous laminates. A variant of the Hellinger–Reissner formulation with Murakami’s zig-zag function increases the percentage error by an order of magnitude for highly heterogeneous laminates. Corresponding formulations using the Reissner Mixed Variational Theory (RMVT) show that the independent model assumptions for transverse shear stresses in this theory may be highly inaccurate when the number of layers exceeds three. As a result, the RMVT formulations require extra post-processing steps to accurately capture the transverse stresses. Finally, the relative influence of the zig-zag effect on different laminates is quantified using two non-dimensional parameters.


Introduction
The application of multi-layered composite materials in loadbearing structures is finding widespread application particularly in the aeronautical, marine and renewable energy industries. Reasons include their high specific strength and stiffness, good fatigue resistance and enhanced design freedom on a micro-and macromechanical level. Furthermore, as material costs reduce the use of carbon fibre composites in large-scale automotive applications is expected to grow considerably in the coming years (Lucintel, 2013).
The design of primary load-bearing structures requires accurate tools for stress analysis. When used around areas of stress concentration or in fail-safe design frameworks composite laminates are often designed to have thicker cross-sections. Under these circumstances non-classical effects, such as transverse shear and normal deformation become important factors in the failure event. These considerations mean that Euler-Bernoulli beam and Kirchhoff plate/shell models that underpin Classical Laminate Analysis (CLA) (Jones, 1998) inaccurately predict global and local deformations. Transverse shear deformations are particularly pronounced in composite materials because the ratio of longitudinal to shear modulus is approximately one order of magnitude larger than for isotropic materials (E iso =G iso ¼ 2:6; E 11 =G xz % 140=5 ¼ 28). The analysis of layered composites is also more cumbersome due to transverse anisotropy, and interlaminar continuity (IC) conditions on displacement, transverse shear and transverse normal stress fields.
Most notably, transverse anisotropy, i.e. the difference in layerwise transverse shear and normal moduli, leads to sudden changes in slope of the three displacement fields u x ; u y ; u z at layer interfaces. This is known as the zig-zag (ZZ) phenomenon. In fact, Carrera (2001) points out that ''compatibility and equilibrium, i.e., ZZ and IC, are strongly connected to each other.'' Thus, while Demasi (2012) showed that the ZZ form of the in-plane displacements u x ; u y and u z can be derived directly from s xz ; s yz and r zz continuity, respectively. Therefore, an accurate model for multilayered composite and sandwich structures should ideally address the modelling issues named C 0 z -requirements by Carrera (2002Carrera ( , 2003b: 1. Through-thickness z-continuous displacements and transverse stresses i.e. the IC condition. 2. Discontinuous first derivatives of displacements between layers with different mechanical properties i.e. the zig-zag effect. For this purpose high-fidelity 3D Finite Element (FE) models are often employed for accurate structural analysis. However, these models can become computationally prohibitive when employed for laminates with large number of layers, in optimisation studies, for non-linear problems that require iterative solution techniques or for progressive failure analyses. Thus, with the aim of developing rapid, yet robust design tools for industrial purposes there remains a need for further efficient modelling techniques. In this regard particular focus is required on equivalent single layer (ESL) theories because the number of unknowns in these formulations is independent of the number of layers.
One of the earliest examples of ESL theories including nonclassical effects is the First Order Shear Deformation Theory (FSDT) (Timoshenko, 1934;Mindlin, 1951;Yang et al., 1966). However, Whitney and Pagano (1970) demonstrated that FSDT only yields improvements on CLA for global structural phenomenon but does not improve in-plane strain and stress predictions for highly heterogeneous and thick laminates. Furthermore, FSDT produces piecewise constant transverse shear stresses that violate IC and traction-free conditions at the top and bottom surfaces.
To overcome these shortcomings the so-called Higher-Order Shear Deformation Theories (HOT) were introduced. In general, the cross-section is allowed to deform in any form by including higher-order terms in the axiomatic expansions of the in-plane displacements u x and u y . Vlasov (1957) refined Mindlin's theory by guaranteeing that transverse shear strains and stresses disappear at the top and bottom surfaces in the absence of shear tractions. Taking Vlasov's condition into consideration, Reddy (1983) presented a higher-order shear deformation theory by expanding the in-plane displacement field to a third order polynomial in z. A large number of different shear shape functions have been published in the past ranging from polynomial (Ambartsumyan, 1958a;Reissner, 1975;Reddy, 1986) to trigonometric (Levy, 1877;Stein, 1986;Touratier, 1991;Karama et al., 1998;Ferreira et al., 2005), hyperbolic (Soldatos, 1992;Neves et al., 2013) and exponential (Karama et al., 2003;Mantari et al., 2011).
Transverse normal strains may be incorporated by extending the expansion of the out-of-plane displacement u z to yield a class of theories denoted as Advanced Higher-Order Theories (AHOT). Here, the class of theory is often denoted by fa; bg where a refers to the order of expansion of the in-plane displacements u x and u y , and b to the order of the transverse displacement u z . Examples of such theories are given by Tessler (1993), Cook and Tessler (1998) and Barut et al. (2001) but these theories generally only provide improvements that are worth their additional computational effort for sandwich panels with compliant, thick cores (Demasi, 2012) or when one face laminate is considerably stiffer than the other (Gherlone, 2013).
All of the previously discussed theories are based on displacement formulations where the displacements u x ; u y ; u z are treated as the unknown variables, and the strains and stresses are derived using kinematic and constitutive equations, respectively. In this case, the governing field and boundary equations may be derived using the principle of virtual displacements (PVD). Being formulated on a displacement-based assumption the transverse shear stresses typically do not guarantee the C 0 z -requirements a priori. More accurate transverse stresses are recovered a posteriori by integration of the in-plane stresses in Cauchy's 3D indefinite equilibrium equations (Whitney, 1972).
This post-processing operation can be precluded if some form of stress assumption is made. One class of model is based on applying the Hellinger-Reissner mixed variational principle. Here the strain energy is written in complementary form in terms of in-plane and transverse stresses, and the transverse equilibrium equation is introduced as a constraint condition using a Lagrange multiplier (Reissner, 1944(Reissner, , 1945.  and  used the Hellinger-Reissner mixed variational theorem to develop a higher-order theory for studying vibrations and plane waves in piezoelectric and anisotropic plates, accounting for both transverse shear and transverse normal deformations with all functional unknowns expanded in the thickness direction using orthonormal Legendre polynomials. The researchers showed that the major advantage of the Hellinger-Reissner theory is that by enforcing stresses to satisfy the natural boundary conditions at the top and bottom surfaces, and deriving transverse stresses from the plate equations directly, the stress fields are closer to 3D elasticity solutions than a pure displacement-based equivalent that relies on Hooke's law to derive the stress fields. In particular this means that boundary layers near clamped and free edges, and asymmetric stress profiles due to surface tractions on one surface only can be captured accurately. Cosentino and Weaver (2010) applied the Hellinger-Reissner principle to symmetrically laminated straight-fibre composites to develop a single sixth-order differential equation in just two variables: transverse deflection w and stress function X. The formulation of this theory is less general than the one proposed by  as its aims are to realise accurate 3D dimensional stress predictions for practical composite laminates at minimum computational cost.
Later, Reissner (1984) had the insight that when considering multi-layered structures, it is sufficient to restrict the stress assumptions to the transverse stresses because only these have to be specified independently to guarantee the IC requirements. This variational statement is known as Reissner's Mixed Variational Theory (RMVT), which makes model assumptions on the three displacements u x ; u y ; u z and independent assumptions on the transverse stresses s xz ; s yz ; r zz . Compatibility of the transverse strains derived from kinematic and constitutive equations is enforced by means of Lagrange multipliers.
HOT and AHOT provide considerable improvements in terms of transverse stress profiles and accurate modelling of global structural effects. However, these theories are not capable of explicitly capturing ZZ effects as the in-plane variables u x ; u y are defined to be at least C 1 z -continuous. In this regard ESL theories that incorporate ZZ kinematics present a good compromise between local, layerwise accuracy and computational cost. Based on an historical review of the topic by Carrera (2003a) the ZZ theories can generally be divided into three groups: 1. Lekhnitskii Multilayered Theory (LMT). 2. Ambartsumyan Multilayered Theory (AMT). 3. Reissner Multilayered Theory (RMT). Lekhnitskii (1935) appears to be the first author to propose a ZZ theory originally formulated for multilayered beams. This was later extended to the analysis of plates by Ren (1986a,b). Ambartsumyan (1958a,b) developed a ZZ-theory for symmetric, specially orthotropic laminates by making a parabolic assumption for the transverse shear stresses. The corresponding displacement field is derived by integrating the shear strain kinematic equation through the thickness and solving for the integration constants by enforcing displacement-IC. Whitney (1969) later extended the analysis to symmetric laminates with off-axis plies and noted that the theory provides excellent results for gross laminate behaviour when compared to the 3D elasticity solutions of Pagano (1969Pagano ( , 1970a. However, Whitney also noticed that the theory fails to accurately capture the slope discontinuity at layer interfaces for the transverse shear stresses. Di Sciuva (1984Sciuva ( , 1985 introduced a displacement based ZZ theory where piece-wise linear, ZZ contributions in the thickness direction enhance a FSDT expansion for u x and u y . The slopes of the layerwise ZZ functions are obtained by enforcing the same transverse shear stress for all layers and by defining the ZZ function to vanish across the bottom layer. As a result, the transverse shear stress in all layers is identical to that of the bottom layer, causing a bias towards the transverse shear stiffness of this layer. To overcome this counterintuitive property Averill (1994) and Averill and Yip (1996) introduced a penalty term in the variational principle that enforces continuity of the transverse shear stresses as the penalty term becomes large. Tessler et al. (2007) note that the formulations based on Di Sciuva's early works present two major issues: 1. The in-plane strains are functions of the second derivative of transverse deflection w. This fact means less attractive C 1 continuous shape functions of w are required for implementation in FE codes. 2. The physical shear forces derived from the first derivatives of the bending moments are different from the shear forces derived by integrating the transverse shear stresses over the cross-section.
To remedy these drawbacks Tessler et al. developed a refined zig-zag theory (RZT) (Tessler et al., 2007(Tessler et al., , 2009(Tessler et al., , 2010a. The kinematics of RZT are essentially those of FSDT enhanced by a zig-zag field w a ðx; yÞ multiplied by a piecewise continuous transverse function / k a , u ðkÞ a ðx; y; zÞ ¼ u a þ zh a þ / ðkÞ a ðzÞw a for a ¼ x; y; ð1aÞ u z ðx; y; zÞ ¼ wðx; yÞ: In this theory the ZZ slopes b x ¼ @/ ðkÞ x =@z and b y ¼ @/ ðkÞ y =@z for u x and u y , respectively, are defined by the difference between the transverse shear rigidity of a layer G az , and the effective transverse shear rigidity G of the entire layup where t k and t are the thickness of the layer k and total laminate thickness, respectively. RZT has shown excellent results compared to the 3D elasticity solutions by Pagano (1969Pagano ( , 1970a for both general composite laminations and sandwich constructions. Recently RZT has also been expanded to include transverse normal stretching and higher-order displacements for a ZZ theory of order {2, 2} (Barut et al., 2012). Similarly, Murakami (1986) enhanced the axiomatic FSDT expansion by including a zig-zag function, herein denoted as Murakami's ZZ function (MZZF), that alternatively takes the values of +1 or À1 at layer interfaces. Therefore the slope purely depends on geometric differences between plies and is not based on continuity of transverse shear stresses. In addition Murakami made independent, piecewise parabolic assumptions for the transverse shear stresses and applied RMVT to obtain the governing field equations. In recent years the MZZF has been applied to functionally graded materials (Neves et al., 2013), sandwich structures (Brischetto et al., 2009a,b) and in the framework of the Carrera Unified Formulation (Carrera et al., 2013a,b) for static and dynamic analyses. Carrera (2004) investigated the effect of including the MZZF in first-order and higher-order displacementbased and mixed-variational theories, showing that superior representation of displacements and stresses, combined with less computational cost can be achieved by including a single ZZ compared to a higher-order continuous term. On the other hand, Gherlone (2013) showed that the MZZF leads to inferior results than RZT for sandwiches with large face-to-core stiffness ratios and laminates with general layups. Thus, an accurate choice of the ZZ function seems to be of paramount importance.
A multiscale approach for modelling the multifaceted structural behaviour of composite laminates in one unified model has been proposed by Williams (2005). The theory uses a general framework with non-linear von Kármán displacement fields and additional temperature and solute diffusion variables on two lengthscales, global and local levels, with the transverse basis functions of the two lengthscales enforced to be independent. This results in two sets of variationally consistent governing equations such that the theory is capable of capturing, in a coupled fashion, the thermomechanical-diffusional phenomena of laminates at the micro, meso and macro levels simultaneously. The use of interfacial constitutive models allows the theory to model delamination initiation and growth, as well as non-linear elastic or inelastic interfacial constitutive relations in a unified form. Williams (2001) has shown that multi-lengthscale theories can be more computationally efficient than pure layerwise models as the order of theory can be increased on both the local and global level. The displacement-based theory features 3ðV g þ V l NÞ unknowns where V g and V l are the global and local number of variables and N the number of layers. In general V g ¼ V l ¼ 3 is sufficient for accurate 3D stress field predictions as derived from constitutive relations. In later work Williams (2008) developed an improved formulation by deriving the governing equations from the method of moments over different length scales and enforcing the interfacial continuity of transverse stresses in a strong sense.
The aim of this paper is to provide additional insight into the fundamental mechanics of the ZZ phenomenon and develop a computationally efficient analysis tool for industrial applications. The aim is not to develop a unified general theory as presented by Carrera (2003b), Demasi (2008),  or Williams (2005) but to provide evidence that non-classical effects due to highly heterogeneous laminations and their associated 3D stress fields can be captured adequately and efficiently by global third-order moments combined with a local ZZ moment. Section 2 outlines a physical explanation for the source of ZZ-displacement fields based on an analogous system of mechanical springs, which to the authors' knowledge, has previously not been given in this manner. In Section 3 these insights are used to combine Reddy's transverse shear function with an Ambartsumyan-type ZZ theory to derive a displacement-based model (MRZZ). In Section 4, the Hellinger-Reissner mixed variational principle is used to derive a ZZ theory that may be used alongside the RZT ZZ function (HR-RZT) or the MZZF (HR-MZZF). The theory is different from general theories in that in-plane and transverse stress fields share the same variables thereby greatly reducing the number of unknowns. In Section 5 the MRZZ, HR-RZT and HR-MZZF theories are compared to Pagano's exact 3D elasticity solutions (Pagano, 1969) in a simply-supported bending load case of a thick beam. Furthermore, the performance of the Hellinger-Reissner principle and the Reissner Mixed Variational Theory are compared. Finally, the influence of transverse shear, transverse normal and ZZ effects on bending deformations are analysed.

Origin of zig-zag displacements
In the following discussion our analysis is restricted to 1D composite beams with negligible transverse normal strain and stress effects. Thus, consider an N-layer composite beam of arbitrary constitutive properties as depicted in Fig. 1. The beam may be of entirely anisotropic or sandwich construction and is subjected to bending moments or shear forces causing it to deflect transversely to the stacking direction. In all cases the x-direction is defined to be along the principle beam axis (0°fibre-direction) while the z-axis is in the transverse stacking direction. Individual layers are prevented from sliding such that the IC conditions for the displacement field u x and the transverse shear stress s xz are satisfied: s ðkÞ where superscripts and subscripts ðkÞ indicate layerwise and interfacial quantities, respectively. If the composite beam comprises layers with different transverse shear modulii then the IC condition on transverse stress inherently results in discontinuous transverse shear strains across ply interfaces. Assuming linear geometric deformation, the kinematic relation for the transverse shear strain is given by where the comma notation denotes differentiation. As transverse normal strain is assumed to be negligible, i.e. u z is constant for all layers, discontinuous transverse shear strains require a change in u x;z across ply interfaces. Thus the slope of the displacement field u x in the thickness direction must change at ply interfaces giving rise to the so-called ''zig-zag'' displacement field. This effect is depicted graphically by the in-plane displacement (u x ), and transverse shear stress (s xz ) and strain (c xz ) plots through the thickness of a [90/0/90/0/90] laminate in Figs. 2(a) and (b), respectively. Here, the transverse shear modulus G xz of the 90°layers are 2.5 times less than the value of the 0°oriented layers causing a step change in transverse shear strain at the ply interfaces. Fig. 2 also shows an example of a laminate with ''Externally Weak Layers'' (EWL). As discussed by Gherlone (2013) these laminates have external layers (k ¼ 1 or k ¼ N) with transverse shear modulii lower than the adjacent internal layers (k ¼ 2 and k ¼ N À 1, respectively), and do not appear to have a ZZ displacement at these interfaces. Gherlone attributed this phenomenon to the stiffer inner layers dominating the more compliant external layers. Furthermore, these phenomena cannot be captured by RZT such that the layer material properties need to be altered artificially as follows: EWL Implementation in RZT.
However, there is, in fact, a slope discontinuity at the interfaces of the EWL. This discontinuity is considerably less pronounced than at the interface between the internal 0°and 90°layers and is consequently not noticed as easily. This phenomenon may be explained by observing the general shape of the transverse shear stress and shear strain profile of a [90/0/ 90/0/90] laminate as shown in Fig. 2(b). The transverse shear stress at the interface between the outer layers is an order of magnitude less than the transverse shear stress at the inner interfaces. Therefore the discontinuity in transverse shear strain is much larger for inner layers than for outer layers, making it appear that there is no ZZ effect for the EWL. Even though the ratio of shear strains at the outer and inner [90/0] interface remains the same, the difference in magnitude is considerably larger for the inner layers. It is this difference in transverse shear strains, rather than the ratio that drives the slope discontinuity of the displacement field.
This also means that discontinuities in transverse shear strain for EWL laminates such as The difficulty in accurately modelling the ZZ phenomenon is that the displacement and transverse shear stress fields are interdependent. As shown in Fig. 2 the layerwise slopes of the ZZdisplacement field u x depend on the transverse shear stress distribution. At the same time the transverse shear stress is a function of the kinematic equations. The difficulty in axiomatic, displacementbased theories is that the ad hoc assumptions for u x and u z need to derive accurate transverse shear stresses, if the IC on s xz is to be used to define layerwise ZZ slopes. Similarly, Ambartsumyan-type models (Ambartsumyan, 1958a) need to include all pertinent variables that influence the distribution of the transverse shear stress to derive an accurate through-thickness distribution for u x .

Spring model for zig-zag displacements
The IC requirements on in-plane displacements and transverse shear stresses are mechanically similar to a combined system of ''springs-in-series'' and ''springs-in-parallel'' (see Fig. 3). For example, a set of springs in series acted upon by a constant force extends the springs by different amounts. By analogy, a constant transverse shear stress acting on a laminate with layers of different shear modulii results in different shear strains in the layers. This represents a smeared, average value of the actual piecewise, parabolic transverse shear distribution. At the same time a system of springs in parallel elongated by a common displacement develops different reaction forces in the springs. This case may be interpreted as layerwise transverse shear stresses following the path of highest stiffness. Conceptually, these two spring systems combine to capture the interplay between transverse shear stress and strain as influenced by the IC conditions.
The average transverse shear stress condition of the ''springsin-series'' model is expressed via Hooke's law as an effective shear modulus G multiplied by an average shear strain c xz The effective shear modulus G is found using the stiffness equation of a set of springs in series Note that the shear modulus of each layer is normalised by the layer thickness fraction to guarantee that G ¼ G xz for a laminate with layers of equal shear modulii. It is worth emphasising that the effective stiffness G is the same as the expression for G a in Eq. (2) found by Tessler et al. in RZT. The change in displacement slope u x;z at layer interfaces depends on the differences in transverse shear strain at interfaces. By inserting Eq. (5) into the transverse shear constitutive equation, we see that the transverse shear strain is a function of the layerwise stiffness ratio g ðkÞ ¼ G=G ðkÞ xz . This ratio is used to capture the differences in layerwise displacement slopes. Fig. 2(b) shows that the shear stress profile of a multi-layered beam differs from a single layer beam in that the z-direction curvature of the transverse shear stress profile in the stiffer 0°plies is increased whereas the curvature in the more compliant 90°is reduced. Integrating the axial stress r x derived from CLA for a zero B-matrix laminate in Cauchy's equilibrium equations, shows that the magnitude of the quadratic term z 2 is influenced by the transformed layer stiffness Q ðkÞ , noting that j x is the flexural curvature. The ''springs-in-series'' analogy is now used to define an effective in-plane stiffness E, The change in layerwise z-direction curvature of the transverse shear stress profile is a function of the relative magnitude of Q ðkÞ to the equivalent laminate stiffness. Therefore, a layerwise in-plane stiffness ratio e ðkÞ ¼ Q ðkÞ =E is defined to quantify the change in transverse shear stress curvature of each layer.

Modified Reddy zig-zag theory
In this section the laminate stiffness ratios g ðkÞ and e ðkÞ are used to derive a new Ambartsumyan-type ZZ theory by modifying Reddy's polynomial shear function (Reddy, 1983) to account for the ZZ effect. The use of the shear function guarantees that transverse shear stresses disappear at the top and bottom surfaces and that transverse shear stresses are continuous at layer interfaces. The model also resolves the issue of modelling ''Externally Weak Layers'' that was addressed by Gherlone (2013) and provides accurate a priori transverse shear predictions for composite laminates with zero B-matrix terms. Being a displacement-based method,  the governing equations are derived using the principle of virtual displacements.

Transverse shear stress and displacement fields
As a starting point Hooke's Law of Eq. (5) is enhanced by the continuous, parabolic shear stress profile proposed by Reddy (1983), A numerical investigation of various composite and sandwich laminates using Pagano's 3D elasticity solution was performed.
This showed that the layerwise z-direction curvatures of s xz can be quantified empirically by the modification factor m ðkÞ ¼ e ðkÞ ðg ðkÞ þ 1=g ðkÞ À 1Þ. This modification factor reduces to m ðkÞ ¼ 1 for a homogeneous laminate. Eq. (9) is rewritten to account for the differences in z-direction profile curvature for different layers, where the layerwise constants A ðkÞ are found by enforcing transverse shear stress continuity at layer interfaces. Details of this derivation are shown in Appendix A.
From the constitutive relation the layerwise shear strain c ðkÞ xz is given by Next, the displacement field u x ðx; zÞ is derived by integrating the kinematic equation (4) in the thickness z-direction while assuming that the transverse displacement u z ðx; zÞ is constant through the thickness The layerwise constants c ðkÞ are found by enforcing continuity of displacements at layer interfaces and the condition that u x ðx; 0Þ ¼ 0 due to the midplane symmetry of the beam. The derivation of the layerwise constants c ðkÞ is provided in Appendix A.

Derivation of the governing equations
Consider a beam as represented in Fig. 4 undergoing static deformations under a specific set of externally applied loads and boundary conditions. The static behaviour of this structure is analysed using the ZZ displacement fields derived in Eq. (12) by means of the two kinematic unknowns wðxÞ and c xz ðxÞ.
The principle of virtual displacements states that a body is in equilibrium if the virtual work done by the equilibrium forces, when the body is perturbed by a virtual amount dũ from the true configurationũ is zero. With regard to the elastic body depicted in Fig. 4 the virtual work done by the virtual displacement dũ is wherer x andŝ xz are the prescribed stresses at the boundary points x A and x B . The superscript G indicates that the strains are calculated via the geometric strain-displacement relations: where the shear function s ðkÞ ðzÞ ¼ A ðkÞ À 4 t 2 m ðkÞ z 2 and the displacement function f ðkÞ ðzÞ ¼ g ðkÞ ðA ðkÞ z À 4 3t 2 m ðkÞ z 3 Þ þ c ðkÞ . The superscript H indicates that the stresses are calculated via the material constitutive equations (Hooke's Law) r HðkÞ where Q ðkÞ ¼ E ðkÞ x for a beam in plane-stress in the width-direction, and Q ðkÞ ¼ E ðkÞ Furthermore, minimisation of P gives the essential and natural boundary conditions at ends x A and x B : The stress resultants are defined as follows: Àt=2 g ðkÞ s ðkÞ ðzÞs HðkÞ where the beam stiffness constants are calculated using the following integrals: Gg ðkÞ ðs ðkÞ ðzÞÞ 2 dz: M x is the bending moment as defined in CLA whereas the stress resultant L x is a higher-order moment that captures the ZZ behaviour. Similarly,M x andL x are the prescribed moments on the boundary. Q x represents the transverse shear force in the field equations whereasV x is the shear force on the boundary. It is worth mentioning that Eq. (17b) and (17c) show an inconsistency between the two shear forces Q x andV x that is discussed by Groh and Weaver, 2015.

Hellinger-Reissner zig-zag theory
Previous studies (Tessler et al., 2007(Tessler et al., , 2009Whitney, 1972) have shown that accurate transverse stress fields can be obtained in a post-processing step by integrating the axial stresses in Cauchy's indefinite equilibrium equations. It is expedient to perform this step a priori based on an accurate assumption of the axial stresses, and then derive new sets of governing equations using the Hellinger-Reissner mixed variational principle (Reissner, 1945). This approach was recently applied to straight-fibre and variable stiffness composites (Cosentino and Weaver, 2010;Groh et al., 2013) but the authors did not include cubic axial stresses, ZZ effects and the possibility of modelling laminates with non-zero B-matrix terms in their models. The cubic behaviour is independent of the ZZ effect and is driven by the material orthotropy ratio E x =G xz and aspect ratio t=L that causes a ''stress-channelling'' effect towards the surfaces of the beam (Everstine and Pipkin, 1971;Groh and Weaver, 2015). These previous works are extended here to remedy these shortcomings.

Higher-order zig-zag theory
We assume a cubic in-plane displacement field of the form, where u 0 is the reference surface axial displacement, h is the rotation of the beam cross-section, f and n are higher-order rotations, w is the ZZ rotation and / ðkÞ is a pertinent ZZ function. The row vector f ðkÞ / describes the through-thickness displacement variation and U is the vector of in-plane variables, As outlined in the work by Tessler et al. (2009) where G is the equivalent ''springs-in-series'' stiffness defined in Eq. (6). Murakami's ZZ function (MZZF) (Murakami, 1986) is given by where z k m is the midplane co-ordinate of layer k. The layerwise definitions in Eqs. (22) and (23) are re-written in the following general form where g ðkÞ ZZF and c ðkÞ ZZF are the z-coefficient and constant term of either the RZT or MZZF definitions. In Section 5 the accuracy of the RZT ZZ function and the MZZF are compared for a number of laminates and these two implementations are denoted by HR-RZT and HR-MZZF, respectively.
The following in-plane stress resultants can be derived from the displacement field of Eq. (19) where S is the constitutive stiffness matrix between stress resultants F and strains . The beam stiffness constants are calculated using the following integrals, Here ðO; PÞ and L are higher-order moments that capture the ''stress-channelling'' and ZZ effects, respectively. The axial stress field of this higher-order theory, written in terms of the stress resultants F , can be used alongside the Hellinger-Reissner mixed variational principle (HR) to develop a new higher-order ZZ theory that enforces Cauchy equilibrium equations a priori.

Derivation of transverse shear and transverse normal stresses
The axial strain corresponding to the displacement field of Eq. (19) is given by, The constitutive relation between stress resultants F and strains Applying the stress-strain Eq. (14a) in combination with Eqs. (27) and (28) the axial stress is, An expression for the transverse shear stress is found by integrating the axial stress of Eq. (29) in Cauchy's indefinite equilibrium equation, where g ðkÞ / ðzÞ captures the quartic variation of s ðkÞ xz through each ply k of the laminate, The N layerwise constants a ðkÞ are found by enforcing N À 1 interfa- where by definition Q 0 ¼ 0. In the derivation of Eq. (32) the surface traction on the top surface is not enforced explicitly. However, this condition is automatically satisfied if equilibrium of the axial stress field Eq. (29) and transverse shear stress Eq. (30) is guaranteed. As we are dealing with an equivalent single layer the equilibrium equation is integrated through the thickness z-direction, An expression for N ;x is easily derived from Eq. (29), Now the only undefined quantity in Eq. (33) is s ðNÞ xz ðz N Þ and an expression for this is sought using Eq. (30), (32) and (34), such that by substituting back into Eq. (33) we have and as s ð1Þ xz ðz 0 Þ ¼T b the expression in Eq. (36) is satisfied. Thus, as long as Eq. (33) is enforced in the theory the top surface shear traction is automatically recovered.
Next, an expression for the transverse normal stress is derived in a similar fashion. Integrating Cauchy's transverse equilibrium equation yields r ðkÞ where the transverse normal matrix h where by definition where Q is the transverse shear force. An expression for Q ;x is derived by integrating Eq. (30) and substituting Eq. (32) for a ðkÞ , where t ðkÞ is the thickness of the kth layer. An expression for r ðNÞ z ðz N Þ is defined using Eq. (37), (39) and (41), such that by substituting back into Eq. (40) we have r ðkÞ

Hellinger-Reissner mixed variational principle
A new set of equilibrium equations is derived by means of minimising the potential energy functional P defined in Castigliano's Theorem of Least Work. In this case, P is a functional of the stress resultants F that define the internal strain energy of the beam and the work done by external tractions. As shown in the previous section, equilibrium equations (33) and (40) should be satisfied to guarantee that the transverse stresses are recovered accurately. First, the transverse shear force Q is eliminated from Eq. (40) using the moment equilibrium condition, such that equilibrium equation (40) becomes For equilibrium of the system the first variation of the potential energy functional P must vanish in such a manner that equilibrium equations (33) and (46) are satisfied over the whole beam domain x 2 ½x A ; x B . Following the rules of the calculus of variations this condition is enforced using Lagrange multipliers k 1 ðxÞ and k 2 ðxÞ, respectively, and adding these to the variation of functional P.
The transverse normal strain ðkÞ z is derived from Hooke's Law, written in terms of the full compliance matrix S ij in a state of plane strain in y, as this is the condition assumed in Section 5. Thus, The new set of governing equations is derived by substituting all stress and strain expressions Eqs. (29), (44a), (44b), (48) and (49) into Eq. (47) and setting the first variation to zero. The corresponding Euler-Lagrange field equations in terms of the functional unknowns k 1 ðxÞ; k 2 ðxÞ and F are, where the superscript T denotes the transpose of a matrix. The pertinent essential and natural boundary conditions are given by, The governing equations related to dF are written in matrix notation, with each row defining a separate equation. In total there are seven field equations that can be solved by defining the ten boundary equations and four surface tractions. Eq. (50c) are enhanced versions of the CLA constitutive equation M ¼ Dw ;xx , taking into account transverse shearing and transverse normal effects, and the influence of the higher-order moments. The members of matrix g are transverse shear correction factors, while members of q and x are correction factors related to the transverse normal stresses and Poisson's effect of transverse normal stresses, respectively. The members of row vectors v; q t ; q p ; x t and x p are correction factors enforcing transverse shear and transverse normal behaviour related to the surface tractions. Column vectors K only include the Lagrange multipliers k 1 ; k 2 and their derivatives. The full derivation of the governing equations including details of all coefficients are given in Appendix B. Finally, the expressions found in Eq. (51a) can be used to determine the deformation vector U of the reference surface, whereas the second row of Eq. (51b) can be used to find an expression for the bending deflection w.

Results and model validation
The governing equations for the Modified Reddy zig-zag theory (MRZZ) and the Hellinger-Reissner zig-zag theory (HR-RZT and HR-MZZF) were derived for laminated beams in plane-strain because this allows the results to be compared against Pagano's 3D elasticity solution. However, both theories are not restricted to beams and may readily be extended to the analysis of laminated plates.
Consider a multilayered, laminated beam comprising N orthotropic composite layers as illustrated in Fig. 4 with the midplane and normal to the beam aligned with the cartesian x-and z-axes. The layers may be arranged in any general fashion with different ply thicknesses, material properties or material orientations. The beam is assumed to be simply-supported at the two ends x A ¼ 0 and x B ¼ a as shown in Fig. 5, and is considered to undergo static deformation in plane strain under the applied sinusoidal distributed load equally divided between the top and bottom surfaceŝ This boundary value problem is analysed using the governing equations derived in Eq. (15) for the Modified Reddy displacement-based theory and Eq. (50) for the Hellinger-Reissner mixed theory.

Modified Reddy zig-zag theory
A closed form solution for ðwðxÞ; c xz ðxÞÞ is sought that satisfies the boundary conditions at the two ends, The expressions, satisfy the above boundary conditions exactly, where W and C are the bending deflection and shear rotation amplitudes, respectively. Substituting Eq. (53) into the constitutive relations (17) and then into the governing differential Eqs. (15) results in the following algebraic expressions for W and C, The coefficients ðW; CÞ are now used to calculate displacements, strains and stresses using the pertinent formulae in Section 3.2.

Hellinger-Reissner zig-zag theory
For the Hellinger-Reissner zig-zag theories HR-RZT and HR-MZZF variable assumptions that satisfy the boundary conditions, are given by The boundary condition N ¼ 0 at x ¼ a; b in Eq. (55) combined with the absence of surface shear tractionsT b ¼T t ¼ 0 means that the membrane force N vanishes over the whole beam domain. Therefore the membrane force amplitude N 0 ¼ 0 in Eq. (56a) and equilibrium Eq. (50a) need not be considered. Substituting Eq. (56) into the governing differential equations (50b) and (50c) results in six simultaneous algebraic equations in six unknowns x ¼ ðM 0 O 0 P 0 L 0 k 1 0 k 2 0 Þ T . These equations are readily solved by matrix inversion, where the stiffness matrix K is comprised of the coefficients of the F ; k 1 and k 2 terms in Eqs. (50b) and (50c) and the column load vector q is comprised of the terms associated withT b ;T t ;P b andT t .

Numerical results
The MRZZ and HR-RZT theories introduced within are compared against the 3D elasticity solution by Pagano (1969) for various symmetrically and arbitrarily laminated composite and sandwich beams. Even though the 3D elasticity solution by Pagano was developed for cylindrical bending of an infinitely wide plate, the solution is equally applicable to beams under plane strain. In order to emphasise the effects of transverse shear and ZZ deformability, relatively deep beams of length-to-thickness ratios a=t ¼ 8 are considered for all stacking sequences. The material properties and stacking sequences are shown in Tables 1 and 2, respectively. Material p was orginally defined by Pagano (1969) and is representative of a carbon-fibre reinforced plastic whereas material m features increased transverse stiffness and is based on the work by Toledano and Murakami (1987). Material pvc is a closed-cell polyvinyl chloride foam modelled as an isotropic material. The honeycomb core h is modelled as transversely isotropic and features significantly lower transverse shear stiffness than material p to exacerbate the ZZ effect. As all results are normalised data, the Young's and shear modulii in Table 1 are presented in nondimensionalised form, in this case with respect to the shear modulus G 12 of material h.
The stacking sequences in Table 2 are split into a group of zero B-matrix laminates A-H and general laminates I-M. Laminates A-D are symmetric cross-ply composite laminates with 0 and 90 layers progressively more dispersed through the thickness. Even though thick blocks of 0 and 90 plies (as in Laminate A) are not commonly used in industry due to transverse cracking issues, this sequence maximises the ZZ effect for validation purposes. Laminates E-G are symmetric thick-core sandwich beams with unidirectional or cross-ply outer skins. Laminate G may be considered as a challenging test case in that the sandwich construction maximises the ZZ effect, the stacking sequence is a combination of three distinct materials and the 90 outer plies act as Externally Weak Layers (EWL). Laminate J is an example of an anti-symmetrically laminated beam with zero B-matrix terms. As Pagano's 3D elasticity solution does not include modelling of off-axis, anisotropic layers the AE45 plies were modelled with effective orthotropic material properties using the transformed axial modulii Q ðkÞ x and transverse shear moduli G ðkÞ xz . Laminates I and J are non-symmetric counterparts to the cross-ply laminates A-D mentioned above. Finally, laminates K-M are highly heterogeneous laminates with general laminations in terms of ply orientations, ply thicknesses and ply material properties, and present a challenging test case for any ESL theory.
Normalised quantities of the bending deflection w, axial stress r x , transverse shear stress s xz and transverse normal stress r z are used as metrics to assess the accuracy of the present theories. These normalised quantities are defined as follows, w ¼  Table 3. In each case errors greater than 3% have been underlined to indicate an error outside the acceptable accuracy margin. For comparison, the table also includes the results of a third-order RMVT implementation using the cubic in-plane displacement assumption of Lo et al. (1977) enhanced with a ZZ variable, combined with the piecewise, parabolic transverse shear stress assumption of Murakami (1986). Both the HR and RMVT formulations have been implemented with the RZT ZZ function and Murakami's ZZ function (MZZF) (Murakami, 1986) to compare their performances. Similarly, the results for the general laminates I-M are shown in Table 4. This table does not include the results for MRZZ as this theory is based on the assumption of zero B-matrix lamination. Finally, to qualitatively compare the stress fields through the full laminate thicknesses, the normalised axial stresses r x and transverse shear stresses s xz are plotted in Figs. 6-18.   Laminates A-E the accuracy of the maximum transverse shear stress s max xz is also within 2.3%. Furthermore, the corresponding Figs. 6-10 show that the axial stress and transverse shear stress for these laminates is accurately captured through the entire laminate thickness. The error in maximum transverse shear stress for the anti-symmetric Laminate H is slightly greater at 4.13%. However, the transverse shear profile in Fig. 13(b) shows that this error is on the conservative side and that the overall profile is captured well. Considering that MRZZ only requires two degrees of freedom, the spring modification factors g ðkÞ and e ðkÞ appear to be efficient methods of achieving both accurate axial and transverse shear stress solutions for zero B-matrix laminates. This result is especially noteworthy as most displacement-based models with the principle of virtual displacements do not obtain accurate transverse shear stresses a priori and require a post-processing stress recovery step (Carrera, 2003a).
However, the errors in s max xz for the honeycomb sandwich beams F and G are excessive. The corresponding Figs. 11(b) and 12(b) show that the MRZZ shear function is unable to capture the reversal of the transverse shear stress field in the stiffer face layers. This effect happens because the shear stress assumption Eq. (10) does not allow the magnitude of transverse shear stress to decrease between the outer surface and the laminate mid-plane. This behaviour only occurs for extreme cases of transverse orthotropy, when the transverse shear rigidity of an inner layer is insufficient to support the peak transverse shear stress of the adjacent outer layer. In essence, it is a load redistribution effect that occurs because the shear force must remain constant for any laminate under the same external load. Thus, as the transverse shear orthotropy between different layers is increased the transverse shear stress is increasingly shifted towards the stiffer layers. The extreme case of transverse orthotropy occurs when stiffer outer layers are essentially bending independently with fully reversed transverse shear profiles (i.e. the inner layer carries no transverse shear). To capture this effect an extra z-term with an appropriate coefficient would have to be added to the MRZZ shear stress assumption of Eq. (10). Figs. 11(b) and 12(b) also show that the shear force, being the integral through the thickness, is much less for the MRZZ formulation than for Pagano's solution. This occurs because the transverse shear stress profile is dominated by the smallest value of G ðkÞ xz in the laminate as a result of the reciprocal sum definition of G in Eq. (6).
On the other hand, the corresponding MRZZ axial stress fields for Laminates F and G in Figs. 11(a) and 12(a) remain accurate throughout the entire thickness. Thus, more accurate transverse shear stress profiles for laminates F and G could be obtained in a post-processing stress recovery step although these would not satisfy the original beam equilibrium equations. Finally, the axial stress plots in Figs. 8(a), 12(a) and 13(a) show that the slope of the externally weak 90 layers is rigorously and accurately captured using MRZZ without the need for the RZT implementation rule described in Section 2.1.
For all analysed laminates the accuracy of HR-RZT is within 1% for all three metrics w; r max x and s max xz . The corresponding through-thickness plots in Figs. 6-18 show that both axial stress and transverse shear stress profiles are closely matched to Pagano's 3D elasticity solution for any type of laminate. Most importantly, the transverse shear stress profile is captured accurately from the a priori model assumption. Moreover, the throughthickness plots of r z in Figs. 19 and 20 show that the transverse normal stress field is also captured accurately. Thus, the HR-RZT formulation is shown to be an ESL theory with seven unknowns that provides 3D stress field predictions to within nominal errors of Pagano's 3D elasticity solution for arbitrarily laminated, thick (thickness-to-length ratio 1:8), anisotropic, composite and sandwich beams without the need for post-processing stress recovery steps.
The accuracy of the HR-MZZF formulation is within the same range as HR-RZT for most laminates, and for cross-ply Laminates A,B,D and I the results are identical. Small discrepancies exist for cross-ply laminates Laminates C and J because of the presence of EWLs which are not taken account of in the MZZF. For laminates with at least three unique plies, arising either due to different material properties or varying ply orientations, the HR-MZZF formulation generally gives less accurate results for all three metrics metrics w; r max  Tables 3 and 4 suggest that HR-MZZF captures the maximum value of the transverse shear stress accurately for all laminates. However, the throughthickness shear stress profiles indicate that this is generally not the case. For example, in Fig. 12(b) the transverse shear stress is accurately captured in the outer 0 p-layers of Laminate G whereas there are visible discrepancies in the stress profile for all other layers. The authors want to emphasise that while the results of HR-MZZF are not as accurate as the HR-RZT formulation, overall the results are satisfactory given the challenging laminates considered here.
The HR a priori model assumption of transverse shear stress provides superior results to the model assumption in the RMVT formulation. The transverse shear stress profiles for laminates with a small number of layers, such as A, B, I and M, follow Pagano's solution closely for both RMVT-RZT and RMVT-MZZF. As the number of layers is increased the transverse shear profiles of both formulations oscillate around the 3D elasticity solution, most clearly shown in Fig. 9(b). In Laminates E, F, G, H and L the oscillations in the RMVT-MZZF solution significantly increases the maximum value of the transverse shear stress s max xz as indicated by the 2700% error for Laminate F in Table 3. In RMVT two independent assumptions are made for the displacement and transverse shear stress fields which are enforced to be kinematically compatible in the variational statement. However, as was shown in Section 2, the ZZ effect in the displacement field is directly related to the presence of C 1 discontinuous transverse shear strains that result from continuity requirements on transverse shear stresses at layer interfaces, and as such, the independence of the two fields is not absolute.
Whereas the minimisation of the strain energy in RMVT guarantees that the geometric and assumed model shear strains are compatible there is no such condition on equilibrium between the axial stress and the transverse shear stress. In the HR principle the situation is reversed in that compatibility of geometric and assumed model strains is not guaranteed whereas equilibrium of stresses is enforced. In terms of deriving accurate stress fields, which are the critical measures for failure analyses, it seems that the enforcement of equilibrium is more critical than displacement compatibility, and hence HR seems to be a better formulation than RMVT.
For most laminates the through-thickness profiles of fields may be integrated in the equilibrium equations to compute more accurate transverse stresses, a step which was advised by Toledano & Murakami in their original papers on RMVT (Toledano and Murakami, 1987) and later reinforced by Carrera (2000). While the HR formulation introduced here features seven functional unknowns, the RMVT formulation only features six variables. Thus, the overall computational efficiency of the RMVT formulation with respect to the HR theory depends on the computational effort involved in the extra post-processing step.
It is also observed that for all EWL laminates the RMVT-RZT formulation considerably improves the accuracy of s max xz compared to RMVT-MZZF. As mentioned above, this occurs because the effects of EWLs are not taken into consideration within the MZZF. Furthermore, Figs. 11(a), 12(a) and 17(a) show significant discrepancies between the RMVT-MZZF through-thickness profiles of r x compared to Pagano's solution. Combined with the greater accuracy of HR-RZT compared to HR-MZZF, this corroborates the findings of Gherlone (2013) that the MZZF may lead to inferior results compared to the RZT ZZ functions for general laminates. In fact Toledano & Murakami point out that the ''inclusion of the zig-zag shaped C 0 function was motivated by the displacement microstructure of periodic laminated composites'' and that ''for general laminate configurations, this periodicity is destroyed'' such that the ''theory should be expected to break down in these particular cases'' (Toledano and Murakami, 1987). The present authors want to emphasise that the MZZF results in nominal errors for most practical laminates when employed in a third-order theory coupled with the Hellinger-Reissner mixed variational statement. However, for sandwich beams with very flexible cores or highly heterogeneous laminates the constitutive independence of the MZZF may lead to larger errors.
In conclusion, the HR-RZT formulation is the most accurate of the formulations investigated here for predicting bending deflections, axial bending stresses and transverse bending stresses from a priori model assumptions as a result of guaranteeing that stresses equilibrate. The RMVT-RZT theory provides similar accuracy for axial stresses but requires a posteriori stress recovery step for accurate transverse shear stresses. In terms of computational efficiency, there is a tradeoff between the extra degree of freedom in the HR-RZT formulation and the post-processing step in the RMVT-RZT formulation.

Assessment of transverse shear, transverse normal and zig-zag effects
Previous authors (Murakami, 1986;Tessler et al., 2009;Carrera, 2004) have shown that ignoring the ZZ effect may lead to significant underestimations of the peak axial and transverse stresses. The inaccuracies are especially pronounced for sandwich beams because of the large degree of transverse orthotropy between the flexible core and stiff face layers. Even though a relatively large thickness-length ratio of 1:8 was analysed in this study the ZZ effect is important for longer beams as well. Furthermore, a major aim of sandwich construction is to separate the stiff face layers as far as possible for maximum bending stiffness. This means that sandwich panels are often of larger thickness-to-length ratios than most other laminated structures.
However, the ZZ effects for practical composite laminates may not be as significant. Many industrial lamination guidelines prevent the use of thick blocks of same orientation plies to prevent problems associated with transverse cracking. As such Laminates A-C are not representative of typical laminates used in practice and the results for Laminate D in Fig. 9 show that dispersing the 0 and 90 layers through the thickness greatly reduces the ZZ effect.
The relative effects of transverse shear deformation, transverse normal deformation and the ZZ effect may be analysed numerically using the bending displacement magnitude of the Hellinger-Reissner theory. For simplicity, we only examine symmetric laminates (N ¼ O ¼ 0), and ignore the effect of the higher-order moment P. Also the entire surface traction acts on the top surface such thatP b ¼ 0 andP t ¼ q 0 sin px=a. To start, the effect of the ZZ deformation is also ignored such that the relative importance of transverse shear and transverse normal deformation can be compared. In this case, the governing Eqs. (50) of the HR theory reduce to M ;xx þP t ¼ 0 ð59aÞ and boundary condition Eq. (51b) is required to calculate the bending magnitude w 0 from M and k 2 . Using the solution assumption in Eq. (56a) results in where s CLA ¼ 1=D CLA , and g; x and q are equal to g 22 ; x 22 and q 22 respectively, when calculated without the influence of terms related to N; O; P and L. Note, that the normal shear correction factor q does not influence the bending deflection w 0 result in Eq. (61). The bending deflection can be non-dimensionalised into three separate entities by dividing by the factor q 0 s CLA a 4 =p 4 , where w TS and w TN refer to transverse shear and transverse normal deflection, respectively. These three factors are plotted against the thickness to length ratio (t=a) in Fig. 21 for Laminate D. Furthermore, this plot shows a comparison between the total normalised deflection w 0 and Pagano's normalised solution w pag . Laminate D is chosen here to minimise the ZZ effect on Pagano's 3D elasticity solution and allow a fair comparison with w 0 . compressibility effects are to be included their effects have to be factored into the initial assumption of the displacement field in Eq. (19) that underlies the derivation of the theory. Finally, the comparison between w 0 and w pag shows that the HR model maintains accurate predictions of the bending deflection up to very deep aspect ratios of t=a ¼ 0:4. Fig. 21 shows that the effect of w TN is insignificant compared to the magnitude of w TS for a ½0=ð90=0Þ 25 laminate. Although the transverse normal correction factor x varies with layup, a numerical study showed that for practical carbon-and glass-reinforced composite and sandwich panels the magnitude of this factor is always negligible compared to the transverse shear factor g.
Therefore, the transverse normal component is ignored in the assessment of the influence of the ZZ effect on the bending behaviour. From the boundary condition Eq. (51b) it is seen that Lagrange multiplier k 2 ¼ w when all transverse normal correction factors vanish. As such, when ZZ effects are included the bending deflection equals where s ij and g ij are equal to s ij and g ij respectively, when calculated without the influence of terms related to N; O, and P. The bending deflection is non-dimensionalised using the factor q 0 s CLA a 4 =p 4 , thus Note that the bending flexibility s CLA of CLA and s 22 of RZT are generally not equal. The terms are related by a constant of proportionality that is independent of the laminate thickness t. The two components w ZZ CLA and w ZZ TS in Eq. (64) can be compared to the corresponding bending components that ignore the ZZ effect in Eq. (62). Thus the total deflection ratio r w ¼ w ZZ 0 = w 0 and shear deflection ratio r TS ¼ w ZZ TS = w TS are used as metrics to assess the influence of the ZZ effect on the bending displacement. Fig. 22(a) shows that the ratio of transverse shear components is invariant with t=a but may vary with the stacking sequence. Furthermore, the ZZ effect always reduces the magnitude of transverse shear deformation. Fig. 22(b) shows that the ZZ effect reduces the overall bending deflection of all analysed laminates and that this effect is greatest for the two honeycomb core sandwich beams F and G, i.e. stiffening is most for laminates with the greatest ZZ effect. Furthermore, the reduction in bending displacement is highly non-linear in t=a and converges to a constant value as the thickness of the beam approaches the length. The authors believe that this stiffening effect occurs partially due to the lower sensitivity to transverse shear deformation ( Fig. 22(a)) and due to the action of the higher-order moment L. In general, the sum of the two non-dimensionalised coefficients w ZZ CLA and w ZZ TS may be used to gauge the ZZ effect on a specific combined laminate and thickness to length ratio case.

Conclusions
In this paper, the fundamental mechanics driving the zig-zag (ZZ) effect in multilayered structures was analysed. The ZZ effect was attributed to differences in transverse shear strains at layer interfaces that require changes in slope of the deformation field in order to satisfy the kinematic equations. The dual requirement of transverse shear stress and displacement continuity at layer interfaces led to the notion of modelling the transverse shear mechanics of a multi-layered structure using a combination of ''springs-in-series'' and ''springs-in-parallel'' systems. By defining equivalent spring stiffnesses for the transverse shear modulii and axial modulii, Reddy's parabolic shear function was modified using a layerwise curvature modification factor. This modifed Reddy zigzag theory (MRZZ), guarantees that the shear stress assumption satisfies the layer interface continuity conditions, traction free surface conditions and that laminates with Externally Weak Layers (EWL) can be modelled in an appropriate manner. However, the MRZZ theory is only applicable to laminates with zero B-matrix terms. Furthermore, the notion that accurate transverse stresses can be obtained by integrating the axial stress in Cauchy's equilibrium equations led to the development of a second theory using the Hellinger-Reissner (HR) variational principle. In this theory the ZZ function of the refined zig-zag theory (RZT) introduced by Tessler et al. (2007) and Murakami's ZZ function (MZZF) were used to develop two alternative theories HR-RZT and HR-MZZF, respectively. The governing equations for the MRZZ and HR formulations were derived for laminated beams under the plane-strain condition and the theories were validated with respect to Pagano's 3D elasticity solution.
The results for different laminated composite and sandwich beams with zero B-matrix terms show that the bending deflection and axial stress field is captured to within 1.1% percent by the MRZZ theory. MRZZ is also capable of making accurate predictions of the transverse shear stress field from constitutive equations when the profile does not cause local peaks at z-wise positions other than the laminate mid-plane. This latter effect occurs for extreme cases of transverse orthotropy, i.e. when the transverse shear rigidity of an inner layer is insufficient to support the peak transverse shear stress of the adjacent outer layer. In these cases an additional posteriori stress recovery step should be performed for accurate transverse stress prediction. Nevertheless, for practical, symmetrically laminated composite laminates MRZZ provides accurate transverse shear stress fields that satisfy interfacial continuity conditions directly from the constitutive equations at the same computational cost as the Timoshenko theory.
The HR-RZT is the best performing theory considered, predicting the bending deflection and three stress fields r x ; s xz and r z to within 1% of Pagano's solution even for highly heterogeneous laminates with arbitrary thickness ratios, ply material orientations and layer material properties. This result is noteworthy because the  transverse shear and transverse normal stresses are captured accurately without requiring a post-processing stress recovery step. In the HR principle the equilibrium of in-plane and transverse stress fields is enforced in a weak sense such that boundary layers near clamped and free edges, and accurate transverse stresses can be derived directly from the plate equations. Furthermore, by using the same variables for in-plane and transverse stress fields the computational effort is kept low.
The results of the HR-MZZF theory showed that even though the MZZF is not based on the constitutive heterogeneity of a laminate it captures the three-dimensional stress field to similar accuracy of HR-RZT. For some laminates the errors in HR-MZZF compared to Pagano's solution are as great as 10% whereas the HR-RZT formulation remains to within 1% (Laminate G). The performance of the HR formulations was also compared to corresponding theories developed using the Reissner Mixed Variational Theory (RMVT). Whereas the RMVT-RZT and RMVT-MZZF give accurate predictions for the bending deflection and axial stress, the model assumptions for transverse shear stress may be highly inaccurate when the number of layers exceeds three. As a result, the RMVT formulations require extra post-processing steps to guarantee accurate transverse stress results. However, compared to the HR formulation the RMVT formulation reduces the variable count by one. Thus, the overall computational efficiency of the RMVT formulation with respect to the HR theory depends on the effort involved in the extra post-processing step.
Finally, for practical non-sandwich beams used in industry the ZZ effect does not seem to be of great importance as many lamination codes prohibit thick blocks of 0 and 90 plies. In these cases higher-order effects such as ''stress-channelling'' are more important. Furthermore, two non-dimensional factors have been identified that allow the influence of the ZZ effect on the classical bending deflection and transverse shear behaviour to be quantified. The results show that including the ZZ effect in the model reduces the effect of transverse shear deformation and generally acts to stiffen the structure.
Similarly, if the reference surface z ¼ 0 lies on the interface between two layers, defined as k þ and k À , the constants are given by Performing the variations on the functionals in Eqs. (B.1a)-(B.1e) following the rules of the calculus of variations results in the following expressions. For the potential of axial stress we have, For the potential of transverse shear stress, where g is a 5 Â 5 matrix of shear coefficients that automatically includes pertinent shear correction factors, and v is a 1 Â 5 row vector of correction factors that enforce transverse shearing effects of the surface shear tractions. Matrix g and row vector v are defined by For the potential of transverse normal stress we expand the parantheses in Eq. (B.1c), define the following transverse normal correction factors where x and q are 5 Â 5 matrixes and x t ; x p ; q t and q p are 1 Â 5 row vectors, and integrate by parts to give

ðB:7Þ
Finally, the potential of the Lagrange multipliers and the potential of contour loads are given by, x A

ðB:9Þ
The integral expressions in Eqs. (B.2), (B.5), (B.7), (B.8) and (B.9) combine to form the governing field Eqs. (50), while the terms evaluated at x ¼ x A and x ¼ x B combine to form the governing boundary Eq. (51). These equations feature three column vectors K eq ; K bc1 ; K bc2 that include the Lagrange multipliers k 1 ; k 2 and their derivatives. These are given by,