Three-dimensional stress analysis for laminated composite and sandwich structures

Accurate stress prediction in composite laminates is crucial for safe design under different loading conditions. Classical laminated theory, i.e. those based on the Euler-Bernoulli and Kirchhoff hypotheses, respectively for beams and plates/shells are inaccurate for relatively thick laminates as three-dimensional (3D) effects such as transverse shear and normal deformations are neglected. Therefore, 3D finite element models are often employed for accurate stress analysis. However, these models are computationally expensive when used for laminates with a large number of layers, in optimisation studies, or for non-linear analyses. To address this issue, a Unified Formulation approach is presented for the analysis of laminated, slender beam-like structures. To define the kinematic field over the beam's cross-section, a recently developed hierarchical set of expansion functions, based on Serendipity Lagrange expansions, are employed and adapted to the layer-wise approach. The present formulation, which has displacements as degrees of freedom, does not ensure continuous transverse stresses across layer interfaces. Thus, an extra post-processing step is required to capture these stresses accurately. The proposed model is benchmarked against a 3D closed-form solution, 3D finite elements, and results available in the literature by means of static analyses of highly heterogeneous, laminated composite and sandwich beams. A key advantage of the present model is its ability to predict accurate 3D stress fields efficiently, including boundary layer regions, i.e. towards clamped ends. As a result, global analyses (e.g. overall displacements, buckling, etc.) and local analyses (e.g. stress concentrations) are combined within a single, computationally efficient model. The performance of the proposed approach, in terms of computational cost and precision, is assessed. Significant computational efficiency gains over 3D finite elements are observed for similar levels of accuracy.


Introduction
Multi-layered composite structures are widely used in engineering fields such as the automotive, aerospace, marine, sports and health sectors. The primary reasons are the high stiffness-and strength-toweight ratios of these materials. The increasing application of such structural members has stimulated interest in the development of tools for accurate stress predictions. Providing a robust and efficient tool, with advanced modelling and numerical techniques, is one of the major challenges in the field of computational mechanics, as the following issues must be addressed: 1. Severe transverse shear deformations due to high orthotropy ratio (E 11 /G 13 ), which increases the channelling of axial stresses [1], a phenomenon not captured by models with simple kinematics assumptions.
It is because of the aforementioned complexities, amongst others, that high-fidelity finite element methods (FEM) are often employed to obtain reliable three-dimensional (3D) stress analyses with the desired level of accuracy. However, these models are computationally expensive and require a vast amount of computer storage space. Thus, with the aim of developing computationally efficient, yet robust, design tools for the practicing engineer, there remains a need for efficient modelling techniques.
In this context, many efforts have been carried out over recent decades to accurately assess the response of laminated composites. The Euler-Bernoulli beam theory and Kirchhoff plate/shell models that support Classical Laminate Analysis (CLA) [4] are inaccurate for modelling moderately deep laminates with relatively low transverse shear modulus. The inaccuracy arising from transverse shear and normal strains across the laminate cross-section, as well as zig-zag effects in the displacement field approximation, being neglected. The First-Order Shear Deformation Theory (FSDT) [5,6] extends the kinematics of classical theories and captures the effect of transverse shear deformation in an average sense. This modification improves the prediction of the global structural response, such as bending displacements and lowfrequency buckling and vibrational modes, but does not improve the prediction of localised stress/strain. Furthermore, FSDT is limited by its uniform transverse shear strain assumption, and therefore, shear correction factors are needed to adjust the constant through-thickness strain profile. However, determining the magnitude of these shear correction factors is not a straightforward task, especially for highly heterogeneous, thick composite and sandwich laminates. Thus, to account for the actual higher-order distribution of transverse shear stresses through-thickness, and so as to guarantee that these vanish at the top and bottom surfaces when no shear tractions are applied, the so-called Higher-Order Shear Deformation Theories (HSDT) were introduced.
Over the years, several models based on HSDT have been proposed for the analysis of multi-layered composite beams. These models can be formally divided into two broad categories: global approximation theories based on Equivalent Single-Layer (ESL) models and discrete layer approximation theories based on Layer-Wise (LW) models.
The ESL theory condenses the laminate onto an equivalent single layer, such that the number of unknowns is independent of the number of layers. The major advantage of this theory is the significant reduction in the total number of mathematical variables and the required computational effort. Numerous theories based on the ESL concept have been proposed [7][8][9]. Reddy [10] proposed a third-order shear deformation theory for laminated plates, which provides a parabolic distribution of transverse shear stress through-thickness. Furthermore, to account for thickness stretching, i.e. transverse normal deformation, generalized higher-order theories have been developed [11]. For many applications, ESL theories provide an accurate description of the global laminate response. However, they are inadequate for capturing accurate three-dimensional ply-level stresses. This shortcoming is due to the displacement field approximation, which predicts continuous transverse strains across the interface of different material laminates. To overcome this deficiency, several attempts have been made to incorporate changes in the layerwise slopes of the in-plane displacements by employing zig-zag functions. A detailed historic review of zig-zag theories for multilayered structures can be found in Ref. [12]. Employing zig-zag functions within the ESL model yields fairly accurate global stress results. However, it fails to predict accurate ply level stress responses when employed for sandwich structures with large face-tocore stiffness ratios and thick laminates with general layups [13]. Furthermore, depending on the modelling technique chosen, zig-zag functions may encounter difficulties when dealing with variable-stiffness laminates [14].
With an aim to solve the aforementioned issues, many researchers [15][16][17] have adopted LW approaches that assumes separate displacement field expansions for each material layer. This assumption allows for a correct representation of the strain field and an accurate determination of 3D stresses at the layer level. Most of the LW theories available in the literature are displacement-based; meaning that the displacement components are the unknown variables, and all the strains and the stresses are derived from the displacement assumptions using the kinematic and constitutive relations, respectively. This, however, does not guarantee a priori the IC condition on transverse stresses. One way to overcome this limitation is to recover the transverse stresses by integration of the in-plane stresses in Cauchy's 3D stress equilibrium equations [18,19].
Another possible solution is to use a mixed formulation, which posits a simultaneous assumption of displacement and stress fields. Many authors have proposed mixed formulations based on the Hu-Washizu (HW) principle [20], Hellinger-Reissner (HR) principle [21] and Reissner's Mixed Variational Theory (RMVT) [22,23]. In this regard, Groh and Weaver [24] performed a detailed comparison between two mixed theories, namely HR and RMVT. The third-order refined zig-zag theory  derived from the Hellinger-Reissner mixed variational statement (HR3-RZT) is shown to predict accurate 3D stresses for arbitrary straight-fibre laminates. Moreover, being an equivalent single layer theory, it is computationally efficient, as the number of unknown variables is independent of the number of layers considered. Despite the high level of accuracy and efficiency, this model cannot be used as a general analysis tool for industrial applications, due to its inability to model complex geometries and boundary conditions. In addition, the mixed displacement/stress-based models have denser stiffness matrices, unlike the     Of relevance to the present work, a recent contribution to the field of structural mechanics is the Unified Formulation (UF) by Carrera and co-workers [25,26]. The available literature shows the capability of the UF to solve a wide range of structural mechanics problems in an efficient manner. The formulation supersedes classical theories by exploiting a compact, hierarchical notation that allows most classic and recent theories to be retrieved from one, hence unified, model. Importantly, and unlike many classic theories, the Unified Formulation applies to the partial differential equations governing three-dimensional elasticity. Full stress and strain fields are, therefore, recovered by its implementation. Although current implementations are found wanting in this respect, in a UF setting, complex geometries could easily be analysed. This is because the displacement field is expressed by means of classic 1D (beam-like case) and 2D (plate-and shell-like cases) finite element elements that need not be prismatic. Additional expansion functions are employed to approximate 3D kinematics over crosssections (beam-like case) and through-thickness (plate-and shell-like cases). Typical expansion functions include Taylor (TE) and Lagrange (LE) polynomials, exponential and trigonometric functions, or Chebyshev polynomials. Amongst these, TE and LE elements are the most widely adopted. However, TE models incur numerical instabilities when enriched to capture stresses near geometric discontinuities, such as corners; whilst LE models can have slow mesh convergence rates. A recently developed cross-sectional expansion model, namely the Serendipity Lagrange Expansion (SLE) model, overcomes the above limitations combining two of the main features of TE and LE elements, i.e. it is hierarchical and allows for numerically stable cross-sectional refinements via re-meshing. For the sake of brevity, the reader is referred to [27] for more detailed information. Another class of cross-sectional expansions based on Legendre polynomials, the so-called Hierarchical Legendre Expansion (HLE), and developed within the Unified Formulation framework shows similar advanced capabilities [28]. However, compared to HLE, SLE expansions are easier to implement as they are obtained in a straightforward manner from the product of linear two-dimensional Lagrange polynomials. Moreover, the numerical stability of HLE models is yet to be investigated, whereas, as demonstrated in Ref. [27], SLE models remain numerically stable increasing their order.
The aim of this paper is to exploit the capabilities of SLE-based finite elements within the Unified Formulation framework (UF-SLE) for the analysis of laminated composite and sandwich structures. Current focus is on prismatic, beam-like structures with one-dimensional extension. However, the models presented herein are of broader interest because they can be extended to complex, non-prismatic geometries via geometric mapping. A Layer-Wise approach is adopted, and together with the properties of SLE models, i.e. refinement by combined cross-sectional discretisation and hierarchical expansion, both local and global responses are obtained accurately and in a computationally efficient manner. In order to capture through-thickness transverse shear and normal stresses reliably, a post-processing step is employed where the transverse stresses are recovered by integrating the in-plane stresses in Cauchy's 3D indefinite equilibrium equations.
The remainder of the paper is organized as follows: Section 2 provides an overview of the displacement field approximation in a Unified Formulation approach. Then the geometrical and constitutive relations to compute the strains and the stresses for laminated composite structures are discussed. Section 3 presents a brief summary of the one-dimensional finite element formulation. In Section 4, numerical results obtained for various laminated composite and sandwich beams are compared to Pagano's exact 3D elasticity solutions [29] for a simplysupported bending load case. Furthermore, the performance of the UF-SLE model in capturing boundary layer effects near clamped ends is investigated and the results are compared to high-fidelity 3D finite element solutions and those given in the literature. Finally, the advantage of recovering the transverse normal stress, using 3D equilibrium equations, is highlighted. Conclusions are drawn in Section 5.

Preliminaries
Consider a laminated beam of length L, rectangular cross-section of width b and thickness h, composed of N layers. The material properties and the thickness of each layer may be entirely different. The beam is referred to a Cartesian coordinate system (x,y,z), where the y-direction is defined to be along the principle beam axis, while the z-axis is in the transverse stacking direction as shown in Fig. 1. Let θ denote the fibre orientation angle and the subscript k be used to refer to layer k.

Displacement field
The three-dimensional displacement field is given as A typical way to overcome the limitations of classical beam models is to enrich the kinematics of the approximated displacement field. The use of Taylor expansions, for instance, is common to many theories where higher-order terms are included. However, an accurate analysis of multi-layered beams requires the use of more sophisticated models. The present work employs the UF framework where a 3D structure is discretised with a finite number of transverse planes that runs along the structure's longitudinal axis as shown in Fig. 2. For simplicity, the structure's longitudinal axis can be thought of as a beam and the   transverse plane as its cross-section. The beam is discretised with traditional 1D finite elements, as depicted in Fig. 3, and the cross-sectional deformations are approximated using SLE functions. Adopting this expansion model, cross-sections are further discretised using four-noded Lagrange sub-domains (SLE nodes). The displacement field within subdomains can be enriched by increasing the order of the local Serendipity Lagrange expansion. It is important to point out that the crosssectional mesh captures the warping of the cross-section with one set of 2D shape functions, and the axial behaviour is captured by a separate 1D mesh. This approach differentiates the method from classic 3D FEM, where 3D shape functions are used to discretise and model a volumetric brick of the structure.
The cross-sectional displacement field at the i th beam node is expressed as where m is the number of terms which depends on the order of expansion and u iτ are generalized displacement vectors.
This model allows a layer-wise approach to be implemented directly where each layer can be modelled as one sub-domain and the kinematics within each layer (or sub-domain) can be varied hierarchically as depicted in Fig. 4 (where the shading denotes hierarchical functions spanning the sub-domain). The reader is referred to [27] for more detailed implementation and treatment of SLE models.

Stress components
The laminates considered in the present study are assumed to be homogeneous and operate in the linear elastic range. The stress-strain relation for an orthotropic laminate takes the form as given below where C is the transformed material stiffness matrix that depends on the mechanical properties of the laminate material and fibre orientation angle. The coefficients C ij are the transformed elastic coefficients referred to the (x y z , , ) coordinate system, which are related to the elastic coefficients in the material coordinates C ij by the transformation matrix T as given below The coefficients C ij and the matrix T are not included here for sake of brevity, but can be found in Ref. [30].
It is common practice to compute stresses using the constitutive relation as given by Eq. (4). However, this may lead to discontinuities of stresses at the interface of two adjacent layers (particularly in a displacement-based approach) and thus violates traction continuity. Accurate modelling of a laminated structure requires a description of interlaminar continuous transverse stresses (shear and normal components). In order to improve the 3D stress fields predicted by displacement-based models, transverse stresses can be recovered by employing the indefinite equilibrium equations of 3D elasticity and integrating in-plane stresses in the thickness direction. The 3D stress equilibrium equations for the static case, and in the absence of volume forces, are In-plane stresses, σ xx , σ yy and τ xy , and their derivatives are computed conventionally using constitutive relations. Transverse shear stresses, τ xz and τ yz , are recovered from Eqs. (7) and (8) and the transverse normal stress σ zz is calculated afterwards from Eq. (9) as given by where σ z ( ) zz k is the stress value in the k th -layer and σ zz k b is the stress value at the bottom of the k th -layer.
Furthermore, it is noted that, in order to recover the transverse stresses accurately from the stress equilibrium equations, exact derivatives of the in-plane stresses are required. With the hierarchical nature of the SLE model, such accuracy can be achieved by including higher-order terms in the displacement field approximation. This level of accuracy is not possible for conventional 3D FE elements, as linear or quadratic elements are usually employed, the derivatives of the in-plane stresses are obtained by using numerical schemes, such as finite differences, which may not be sufficiently accurate.

Finite element formulation
The UF relies on a displacement-based formulation of the finite element method. The advantage of a finite element discretisation is that arbitrary geometries and boundary conditions can readily be modelled. In this setting, the volume is discretised into a series of N e -noded subdomains (the elements), so that the displacement field can be approximated element-wise by means of local shape functions N i , and generalised nodal displacements, u i , such that In the present work, 1D Lagrange polynomials with cubic shape functions are used to define the finite element as given in Ref. [26]. Finally, by introducing the cross-sectional approximation of Eq. (1) into the FE discretisation along the beam axis of Eq. (11), the displacement field reads Fig. 11. Through-thickness distribution of the normalized transverse normal stress σ zz at 5%, 10%, 15% and 20% from the clamped end A, for laminate 1.
where m denotes the order of the Serendipity Lagrange expansion. With the current methodology, we overcome the limitation on the aspect ratio of a 3D brick element in FE analysis by decoupling the shape functions, along the longitudinal axis (beam) and across the transverse planes (cross-section). Elastic equilibrium is enforced via the Principle of Virtual Displacements, which, in a quasi-static setting, states that where δ denotes virtual variation with respect to displacements, and W int and W ext denote the internal and external work, respectively. By definition, the internal work is the work done by the internal stresses over the corresponding internal virtual strains and is equivalent to the elastic strain energy. Noting that W W int int e e = ∑ and letting l e be the length of the generic beam element and A be the cross-section area, In the UF notation, the internal work can be re-written as The term K τsij e is referred to as the element Fundamental Nucleus. Its explicit form can be found in Refs. [25,26]. Fundamental nuclei may be assembled into a global stiffness matrix following the standard finite element procedure. For the sake of brevity, the derivation of the fundamental nucleus of the loading vector from the virtual variation of the external work is not reported here, but can be found in Ref. [25].

Numerical examples and discussion
The aim of this section is to assess the accuracy and robustness of the Unified Formulation finite element model, based on the Serendipity Lagrange cross-sectional expansions, in analysing laminated composite structures. To do so, a number of benchmark tests are performed. In Section 4.1, results obtained using the present approach are validated against Pagano's closed-form 3D elasticity solution [29] for simplysupported, laminated composite and sandwich beams. As the exact solution is available for an infinitely wide plate subject to cylindrical bending [14], in the present approach, a plane strain condition is enforced by removing some terms from the material stiffness matrix as shown in Appendix A. Section 4.2 highlights the ability of the proposed Fig. 12. Through-thickness distribution of the normalized axial normal stress σ yy at 5%, 10%, 15% and 20% from the clamped end A, for laminate 2. model in predicting the structure's response near clamped ends, i.e. where 3D effects become relevant and computationally inexpensive classic theories are not applicable. In this case the results are compared with high-fidelity, yet computationally expensive, finite element solutions. Finally, Section 4.3 presents through-thickness plots of the transverse normal stress computed for various laminated beams, by employing the stress recovery scheme as described in Section 2.4.
Throughout, axial stress σ yy , transverse shear stress τ yz and transverse normal stress σ zz are normalised as follows

Model validation
The validation is carried out for a relatively thick square crosssection beam of length-to-thickness ratio, L t / 8 = . The beam is aligned with the Cartesian y-axis and the cross-section is in the xz-plane. The layers are arranged in a general fashion with different ply thickness, material properties and material orientations. The beam is simplysupported at the two ends y 0 = and y L = , and loaded by a sinusoidal distributed load, equally divided between the top and the bottom surface, P P q/2 z z , as shown in Fig. 5. It is to be noted that compared to Pagano's original benchmark [29], Groh and Weaver [24] split the sinusoidal load between the top and bottom surfaces to minimise through-thickness normal stretching. This loading condition allowed for a fairer comparison with their equivalent single layer (HR3-RZT) model and demonstrated that it could accurately capture the tractions on the top and bottom surfaces without a priori assumptions. As the HR3-RZT model is also used as a reference solution herein, we use the benchmark with split sinusoidal tractions (although this is not strictly necessary for the present UF-SLE model). The material properties and stacking sequences adopted are shown in Tables 1 and 2, respectively. This combination of materials, ply thickness and stacking sequence are taken from a recent paper by Groh and Weaver [24] on modelling highly heterogeneous laminated beams using the Hellinger-Reissner mixed formulation. The wide range of laminates considered, from simple to challenging, allows the full capabilities of the current formulation to be tested and validated. Material p represents a carbonfibre reinforced plastic, material m a reinforced plastic with increased   transverse stiffness, pvc is a poly-vinyl chloride foam modelled as an isotropic material and h represents a honey-comb core modelled as a transversely isotropic material. The plies made of these materials are stacked together in different combinations to form laminates as presented in Table 2. Laminates A-D are symmetric; I-J are non-symmetric cross-ply composites. Although these are not widely used in industry due to transverse cracking issues, it is a good test case for model validation as the 0°/90°sequence maximises the zig-zag effect. Laminates E-G are symmetric thick-core sandwich construction. In laminate F, the low transverse shear stiffness of material h compared to that of p exacerbates the zig-zag effect. Laminate G is a challenging sandwich construction with a combination of three distinct materials. Finally, laminates K-M represent highly heterogeneous laminated plates.
The UF-SLE model is used for the analyses presented in this section. The structure is discretised with 30 B4 (four-noded) elements along the length; the cross-section is divided into sub-domains (one per layer). Within each sub-domain (Serendipity Lagrange element) a fifth-order expansion is employed. The number of beam elements and the order of expansion in the cross-section are decided by performing a convergence analysis. For the sake of brevity, only the converged results for all the cases are presented in the paper. Table 3 shows the maximum through-thickness normalised axial stress σ yy max and transverse shear stress τ yz max at y L/2 = and y 0 = , respectively. The results obtained are validated against Pagano's 3D elasticity solution and are also compared with those given by Groh and Weaver [24] using the Hellinger-Reissner third-order refined zig-zag theory (HR3-RZT). For all the cases assessed, the accuracy of results obtained with the proposed model is within 0.01%. Out of all the laminates, F and G are challenging constructions, and therefore, are considered as a particularly important test cases for model validation. The normalised axial normal stress σ yy (at y L/2 = ) and transverse shear stress τ yz (at y 0 = ) are plotted in Figs. 6 and 7. From the plots, it is clearly shown that the commonly used Third-order Shear Deformation (TSDT), employed for laminated composites, is incapable of capturing the extreme cases of transverse orthotropy in laminates F and G, where a reversal of the transverse shear stress in the stiffer layers is observed. This stress distribution is due to the low transverse shear stiffness of the inner layer which makes it insufficient to support the peak transverse shear stress of the adjacent outer layer. Moreover, this behaviour cannot be predicted by a Reissner's Mixed-Variational Theory (RMVT) model implemented with zig-zag functions as the stress assumptions used in the variational statement are not inherently equilibrated [24]. However, the present model is able to capture the effect accurately and the results are in excellent agreement with 3D elasticity solution and those obtained by employing the HR3-RZT model. For more detailed comparisons of the 3D stress fields for laminates A-E and H-M, the reader is referred to Appendix B. It is to be noted that all the results presented in this section and in Appendix B are based on a plane strain assumption (to model an infinitely wide plate). This assumption simplifies the problem as there is no effect of the Poisson's coupling (C 31 ) term and the transverse normal-in-plane shear coupling (C 36 ) term, as shown in Appendix A. Thus, the current kinematic fidelity of the model is sufficient to naturally satisfy the stress equilibrium equations and to accurately capture the transverse stresses without employing the stress recovery post-processing step.

Localised stress fields towards clamped ends
To assess the capability of the present formulation in capturing boundary layer effects and localised stress gradients towards boundaries, the second validation example is carried out for a square section laminated beam of length-to-thickness ratio L t / 10 = , clamped at both the ends. The beam is subjected to a uniformly distributed load, equally divided between the top and the bottom surface, P P q/2 z z Fig. 8. A plane strain condition is enforced as described in Appendix A. Two laminates as shown in Table 4 are considered, where laminates 1 and 2 are non-symmetric, composite and sandwich beams, respectively, comprised of materials p and pvc as defined in Table 1.
In the present approach, the beam is discretised using 40B4 elements and a fifth-order expansion is employed within each cross-section element (one element per layer). As Pagano's closed-form solutions are acceptable only for simply-supported beams, 3D FE results from the commercial code, ABAQUS as given in [14], are used as the reference solution. The 3D model, 1 m long, 0.1 m thick and 0.001 m wide, is meshed with 96,000 C3D8R brick elements. To model a plane strain condition, the lateral faces are restrained from expanding and one element is used in the width direction.
Through-thickness distribution of the stress fields, σ yy , τ yz and σ zz , at four locations 5%, 10%, 15% and 20% from the clamped end (y 0 = ) are plotted in Figs. 6 and 7. The results obtained are also compared with those given in Ref. [14] using the HR3-RZT model. The boundary layer effect is clearly shown in these plots as there is a clear change in the stress profiles at different locations from the clamped support, for all three stress fields. In addition to the boundary layer effect, the high orthotropy ratio in a composite laminate causes channelling of the axial     stress towards the surface [1,2]. This effect requires the non-classical complexity of a higher-order model. Fig. 9 shows the through-thickness distribution of the normalised axial stress σ yy for laminate 1. The stresschannelling effect can be clearly observed in the 0°laminates, first and third layer from the bottom, with an orthotropy ratio E G / y yz = 50. However, near the clamped support, the linear behaviour of the zig-zag effect reduces the relative magnitude of these higher-order throughthickness variations. This effect can be observed by looking at the variation of σ yy between 20% in Fig. 9d and the 5% in Fig. 9a. To accurately capture this stress distribution, at least a fifth-order expansion function is required, as employed in the present formulation. In the case of the HR3-RZT model, based on a third-order theory, discrepancies with 3D FE results are observed. In contrast, the present results are in an excellent agreement with the 3D FE solutions.
The through-thickness profiles of τ yz and σ zz for laminate 1 are plotted in Fig. 10. The effect of the boundary layer induced by the clamped support is observed from the transverse shear and normal stress distributions at 5% location as shown in Figs. 10a and 11a, respectively. The clamped boundary condition exacerbates the zig-zag deformations within the laminate which results in the redistribution of the transverse shear and normal stresses across the section. This effect reduces as we move away from the clamped end. The plot for the 20% location in Figs. 10d and 11d presents the converged solution free from boundary layer effects. Similarly, Figs. 12-14 show through-thickness distributions of the three stress fields for the sandwich beam (laminate 2). The flexible core and the stiff face layers with clamped supports make this a challenging test case to analyse. The increasing effect of the zig-zag deformations towards the clamped ends is shown from the stress profiles from the 20%-5% locations.
From the results presented in this section it is evident that the UF-SLE model is capable of accurate stress predictions compared to the HR3-RZT. However, this comparison is incomplete without highlighting the computational cost incurred by the models. Therefore, we compare the degrees of freedom (the number of unknown variables) required to solve the system, which gives an estimate of computational efficiency. The FE model requires 582,498 DOFs (for laminates 1 and 2), the UF-SLE model requires 26,862 DOFs (for laminate 1) and 33,033 DOFs (for laminate 2), and the HR3-RZT model employs only 217 DOFs (for laminates 1 and 2). Clearly, the HR3-RZT model, based on an equivalent single layer approach, is more computationally efficient than the UF-SLE model, followed by the 3D FE approach. However, the inability of the HR3-RZT model to analyse large and complex structures, subject to a variety of loads and boundary conditions, makes it unfit as a design tool for industrial applications. This requirement of solving complex structural problems is rather important and therefore, analysts often use alternative approaches, for example FE techniques. However, the present formulation can be a good compromise between the two numerical models discussed, when problem (or geometrical) complexity and computational efforts are of concern, as depicted in Fig. 15.

Assessment of transverse normal stress via stress recovery
In the previous sections, results for laminated composite and sandwich beams are computed and compared with the analytical and various numerical solutions available in the literature. All the analyses performed were based on the plane strain assumption in the beam's width direction (x-direction). This assumption forces the normal and shear strains with x-components (ε xx , γ xz and γ xz ) to be zero. In order to assess the performance of the present model in predicting the full 3D stress response of the structure, the analyses performed in Section 4.1 are repeated without any plane strain assumption.
The 3D finite element analysis is performed using the commercial code, ANSYS and the results are used as the reference solution. The 3D model is meshed with the SOLID186 (20-noded brick) element and a mesh convergence analysis is performed to define the optimal mesh size for each laminate considered. In the present SLE model, the beam is dicretised with 30B4 elements and a fifth-order expansion function is used within each Serendipity Lagrange element in the cross-section (one element per layer). Figs. 16-19 present the converged solution for through-thickness transverse normal stress obtained from the 3D FEA and UF-SLE models.
From these figures, it can be clearly seen that like other displacement-based weak-form formulations, the SLE model based on the unified formulation approach, is unable to capture the transverse stresses accurately. The zig-zag effect due to the transverse anisotropy and the C 1 -discontinuous displacements field make the transverse normal stress profile discontinuous at the laminar interfaces. This issue can be addressed either by increasing the fidelity of the model, which is a computationally expensive solution, or by employing the stress-recovery scheme, as used in the present case, where the stress equilibrium equations are integrated along the thickness direction as described in Section 2.4. This feature of recovering the transverse stresses from Cauchy's equilibrium equations creates a stronger condition than simply post-processing the stresses from the displacement unknowns in the kinematic and constitutive relations. The stress distribution profiles obtained are denoted by UF-SLE-SR in the plots, where SR denotes stress recovery. Results show an excellent agreement with the 3D FE solutions.
For all of the laminates considered, the number of unknowns required in the UF-SLE model is less than those required in 3D FEA. However, it is believed that comparing models based on DOFs only is not a fair assessment of computational efficiency. Instead, computational time must be the criterion for comparison. Because it is tricky to compare in-house codes with a commercial software, we compare other parameters which directly relate to computational time and memory requirements. For instance, to solve a linear static analysis, the most time consuming steps are the stiffness matrix inversion and multiplication, which further depends upon the solver type. The first choice employs a sparse direct solver as based on the direct elimination of equations (usually the Gaussian Elimination algorithm) and the solution obtained is stable without being affected by the numerical characteristics of the matrix. However, the direct solver demands a significant memory space and a large amount of calculations for large problems, in which case an iterative solver requiring less memory is more desirable (e.g. the Conjugate Gradient method). To give a detailed mathematical insight into these algorithms is beyond the scope of this paper. The reader is referred to [31,32] for more details.
For both cases, the time and space complexities are measured, which quantifies the amount of time and storage taken by an algorithm [32,33]. The time complexity is estimated by counting the number of elementary operations performed and the space complexity is measured by the input size. Both are commonly expressed using a big O notation [34]. These quantities are calculated for a few laminates as shown in Table 5 for the UF-SLE and the 3D FE model. Despite the large number of degrees of freedom in the 3D FE model, the time required for matrix inversion in both models is the same. This result is due to the fact that in the 3D FE model, the stiffness matrix is more sparse than the UF-SLE model. However, the advantage of the present approach becomes clear when memory requirements are considered. Due to the huge number of degrees of freedom in 3D FE, the memory required is 10 times more than the case of the UF-SLE model. Moreover, the direct solver uses computer's RAM for storing the matrix and for performing other operations. If sufficient RAM is not available, the solver must be changed to iterative, which in turn makes the computation more expensive in case of 3D FE compared to the UF-SLE model.

Conclusions
The purpose of the present work is to capture 3D stress fields accurately and with greater computational efficiency than 3D finite element analysis. The Serendipity Lagrange expansion (SLE) model is introduced within the Unified Formulation (UF) framework for the first time to analyse laminated composite and sandwich beam structures. The model is benchmarked against a 3D elasticity solution, 3D finite element solutions and a mixed formulation based on a Hellinger-Reissner third-order refined zig-zag model (HR3-RZT). The findings from the present study can be summarised as follows: 1. The UF-SLE model is sufficient to obtain a Layer-Wise model and therefore captures the zig-zag effect. The beam's cross-section is discretised such that each layer represents a four-node Lagrange element and the precision of the solution is tuned by varying the polynomial order. In contrast, 3D FE models require a large number of elements per layer and furthermore, the condition on the aspect ratio of a 3D brick element increases overall mesh density. 2. For all of the laminates considered, the present model predicts 3D stress fields accurately and the results are in excellent agreement with Pagano's 3D elasticity solution. In most cases, the results are more accurate than those obtained by the mixed beam benchmark problem. 3. The UF-SLE model is a displacement-based layer-wise approach and the HR3-RZT is a mixed-variational equivalent single layer theory. Both models provide similar levels of accuracy and the HR3-RZT is shown to be computationally more efficient. Despite this relative inefficiency, the present approach has significant benefits as it is more general in terms of the variety of structural mechanics problems that can be solved. 4. As the present approach is displacement-based, i.e. the equilibrium of stresses is guaranteed in a weak sense, the inter-laminar continuity condition of transverse stresses is not satisfied. To ensure that the transverse stresses are captured accurately, a posteriori stress recovery step is employed. 5. The proposed model accurately predicts the boundary layer effects that arise due to local variations in the 3D stresses towards clamped ends. The boundary layer intensifies the through-thickness transverse shear and normal stresses. These stresses play an important role in delamination initiation, thus robust numerical models that capture these effects are essential. 6. With the UF-SLE formulation, global stiffness and buckling, as well as detailed localised stress analyses can be performed in a single model. As such, the need for running low-fidelity models for global response, and high-fidelity models for accurate stress predictions is removed. Potentially, the modelling approach for structural analysts in industry could be simplified. 7. All of the above mentioned points are valid for a 3D finite element model. However, the computational efficiency gain obtained with the proposed model in comparison with finite elements is significant. Thus, the combination of accuracy and computational expense makes the Unified Formulation, based on Serendipity Lagrange expansion model, an attractive method for industrial design tools.
The three zero strain entries in the strain vector indicate that we can ignore their associated columns in the stiffness matrix (i.e. columns 1, 5, and 6). If we also ignore the rows associated with the stress components with x-subscripts, the stiffness matrix reduces to a simple 3×3 matrix, as given below In order to model the plane strain behaviour within the UF framework, we use the following material stiffness matrix