Three-dimensional stress analysis for beam-like structures using Serendipity Lagrange shape functions

Simple analytical and finite element models are widely employed by practising engineers for the stress analysis of beam structures, because of their simplicity and acceptable levels of accuracy. However, the validity of these models is limited by assumptions of material heterogeneity, geometric dimensions and slenderness, and by Saint-Venant’s Principle, i.e. they are only applicable to regions remote from boundary constraints, discontinuities and points of load application. To predict accurate stress fields in these locations, computationally expensive three-dimensional (3D) finite element analyses are routinely performed. Alternatively, displacement based high-order beam models are often employed to capture localised three-dimensional stress fields analytically. Herein, a novel approach for the analysis of beam-like structures is presented. The approach is based on the Unified Formulation by Carrera and co-workers, and is able to recover complex, 3D stress fields in a computationally efficient manner. As a novelty, purposely adapted, hierarchical polynomials are used to define cross-sectional displacements. Due to the nature of their properties with respect to computational nodes, these functions are known as Serendipity Lagrange polynomials. This new cross-sectional expansion model is benchmarked against traditional finite elements and other implementations of the Unified Formulation by means of static analyses of beams with different complex cross-sections. It is shown that Serendipity Lagrange elements solve some of the shortcomings of the most commonly used Unified Formulation beam models based on Taylor and Lagrange expansion functions. Furthermore, significant computational efficiency gains over 3D finite elements are achieved for similar levels of accuracy. © 2018 Published by Elsevier Ltd.


Introduction
In engineering design, long slender structures are typically analysed using axiomatic beam models. These models are valid under the premise that the longitudinal dimension of a structure is at least one order of magnitude larger than representative crosssectional dimensions. This geometric feature allows the governing elasticity equations to be reduced from three to one dimension, (with the reference axis coinciding with the beam axis), and in so doing, brings about significant physical insight and computational benefits. The aim of this current work is to build a model capable of capturing three-dimensional stress fields for beam-like structures in a computationally efficient manner. In this regard, many tudinal surface, correction factors are commonly employed in a TB setting ( Timoshenko and Goodier, 1970;Cowper, 1966;Sokolnikoff, 1956 ).
To account for effects that are not captured by classical axiomatic theories, several refined finite element (FE) models have been developed. However, geometric complexities and accurate approximations of the displacement field can lead to computationally expensive models, where a large number of unknown variables is required.
An approach developed by Ladevèze and Simmonds (1996) reduces a three-dimensional (3D) model to a beam-like structure thereby simplifying the 3D elasticity equations. Using this method a beam model can be constructed as the sum of a Saint-Venant part and a residual, higher-order part. In a following work, Ladevèze and Simmonds (1998) used linear shape functions on beams with general cross-section and developed an exact beam theory for calculating 3D displacements and stresses. However, the theory only works if one neglects localised effects that occur at extremities and geometric discontinuities. Surana and Nguyen (1990) developed a two-dimensional (2D) curved beam element using Lagrange interpolating polynomials along the crosssection to describe transverse shear stress distributions in composite structures. Although accurate, the computational cost of Surana's model grows rapidly as the number of composite layers increases. Recently, Groh and Weaver (2016) used highorder Equivalent Single Layer theories to study the stretching and bending of multilayered variable stiffness, anisotropic flat plates. Kameswara et al. (2001) used Taylor series to include displacement components along the cross-section and proposed an analytical solution based on a mixed formulation, whereby, to ensure continuity, transverse stresses are invoked as degrees of freedom (DOFs). Kameswara's model has been employed for static and dynamic analyses of laminated plates and beams.
Another powerful tool to develop structural models is the asymptotic method. In the beam model scenario, the works by Berdichevsky (1976) and Berdichevsky et al. (1992) are among the earliest contributions that exploited the Variational Asymptotic Method (VAM). More recently, Yu et al. (2002Yu et al. ( , 2012 have developed the so called variational asymptotic beam sectional analysis (VABS) which uses VAM to split the 3D elastic problem into a 2D linear problem in the cross-section and a 1D beam problem in longitudinal direction.
Classical approaches have also been enhanced by the Generalized Beam Theory (GBT) for thin-walled structures, as given by Silvestre and Camotim (2002) , where transverse cross-sectional displacements are obtained from the axial ones. In GBT, in order to obtain a displacement representation compatible with classical beam theories, each component of displacement is expressed as a product of two single-variable functions-one depending on the longitudinal position along the reference axis and the other on cross-sectional coordinates. However, since thin plate assumptions are adopted ( Silvestre and Camotim, 2002 ), through-thickness strains are set to be zero and full 3D stress fields cannot be captured. Following on from early implementations of the GBT, many other high-order theories, based on enriched cross-section displacement fields, have been developed in order to describe effects that classical models cannot capture. A complete account of the literature, however, goes beyond the scope of this paper. The keen reader is referred to Carrera et al. (2015) for further details.
Of relevance to the present work, is one of the most recent contributions to the development of refined beam theories and that is the Unified Formulation by  . The formulation provides one-dimensional (beam)  and two-dimensional (plate and shell) ( Carrera, 2002 ) models that extend the classical approximations by exploiting a compact, hierarchical notation that allows most classic and recent formula-tions to be retrieved from one, hence unified , model. The displacement field is expressed over the cross-section (beam case) and through the thickness (plate and shell cases) by employing various expansion functions including Taylor polynomials , Lagrange polynomials , exponential and trigonometric functions ( Carrera et al., 2013 ), Chebyshev ( Filippi et al., 2015 ) and Legendre polynomials ( Pagani et al., 2016 ). Amongst these, Taylor (TE) and Lagrange expansion (LE) models are most widely adopted. TE models are hierarchical and the degree of accuracy with which kinematic variables are captured is enriched by increasing the order of the cross-sectional expansion. On the other hand, LE models are based on cross-sectional discretisations using Lagrange elements of given kinematic order and refinement is obtained by increasing the mesh density, i.e. by increasing the number of Lagrange elements in the cross-section. Both models are found to be accurate and computationally efficient, but have limitations. Namely, 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. Another known limitation of Carrera's Unified Formulation (CUF) is the oscillation of shear stresses along the beam axis that appears if the mesh along the beam length is not sufficiently fine.
In this work, we propose a new approach for the analysis of beam-like structures that overcome all of the above limitations. The approach is based on CUF and, as a novelty, hierarchical Lagrange polynomials are used to define cross-sectional displacement fields. This new element class, called Serendipity Lagrange (SL), is based on the Trunk (or Serendipity) Space which is a polynomial space spanned by the set of monomials ξ where N is the order of the polynomial ( Szabó and Babuška, 2011 ). SL expansions combine two of the main features of TE and LE models, i.e. they are hierarchical and facilitate numerically stable crosssectional refinements via remeshing. The advantages of using SL elements for beams of general cross-section compared with finite elements, Taylor and Lagrange type models are discussed in the following sections. In addition we investigate the effect of collocating beam nodes towards the boundaries using Chebyshev biased grids, which reduce problematic oscillations in numerical solutions (the Runge effect) ( Kreyszig, 2011 ).
The remainder of the paper is structured as follows. Section 2 summarises the derivation of the governing elasticity equations in weak form. In Section 3 , the equations are cast in a form suitable for the Unified Formulation. Section 4 provides an overview of the Taylor and Lagrange expansion models and details of the derivation of the new Serendipity Lagrange expansions. In Section 5 , Chebyshev biased meshes are described. Numerical results obtained for various beams are found in Section 6 , where accuracy, computational efficiency and numerical stability of the proposed model are discussed. Finally, conclusions are drawn in Section 7 .

Finite element formulation
CUF relies on a displacement-based formulation of the finite element method. Fundamental equations are summarised here for completeness and clarity of exposition.
Let us consider an elastic continuum of volume V , embedded in R 3 . In a finite element setting, the volume is discretised into a series of N -noded subdomains (the elements), so that displacement fields of the form  JID: SAS [m5G;February 28, 2018;7:59 ] can be approximated element-wise by means of local shape functions, N i , and generalised nodal displacements, u i , such that

ARTICLE IN PRESS
In the previous expression and throughout remainder of the paper, the Einstein summation convention is implied over repeated indices.
As per the classical theory of elasticity, the stress and strain tensors can be expressed by six-term vectors as These tensors are related through the material's stiffness matrix C by Hooke's law, stating that σ = Cε .
(2) For the sake of brevity, the reader is referred to Reddy (2003) for an explicit definition of the coefficients in C .
Using Eq. (1) , the strain-displacement relationship in its linear form may be recast as Elastic equilibrium is enforced via the Principle of Virtual Displacements, which, in a quasi-static setting, states that where W int and W ext are the internal and external works, respectively, and δ denotes virtual variation with respect to displacements.
By definition, the internal work is the work done by stresses over corresponding virtual strains. Noting that W int = e W (e) int and letting V (e) be the volume of the generic element Substituting (2) and (3) If we now denote body forces per unit volume as g , surface forces per unit area as p , line forces per unit length as q and concentrated forces acting on Q as P , the external work is which is a statement of elastic equilibrium in weak form, where K (e) i j and f (e) are, respectively, the structural stiffness matrix and the generalised load vector of the generic element.

Unified formulation
A typical way to overcome the limitations of classical beam models and to refine the structural analyses that employ them is to enrich the kinematics of the approximated displacement field. The use of Taylor expansions, for instance, is common to many theories where high-order terms are included to enrich the kinematic approximation. In general, the higher the order, the higher the computational effort required. One of the advantages of CUF is that, owing to the notation adopted, beam models of increasing kinematic refinement are readily developed. Let us consider a beamlike structure as shown in Fig. 1 , where the beam extends along the y -axis and cross-sections lie in the xz -plane. In CUF, the beam is discretised along the length with traditional 1D finite elements. Cross-sectional deformations can be approximated using different expansions as explained in Sections 4.1 and 4.2 . Mathematically, this means that the displacement field and its virtual variations may be written as a product of two functions: cross-sectional expansion functions, F ( x, z ), and 1D Lagrange shape functions, N ( y ), along the beam axis. In principle, these functions can have as many terms as desired. The more terms there are, the richer the kinematics. With reference to Eq. (1) , selectfont where M is the number of terms in the cross-sectional expansion depending on the order; N e is the number of Lagrange nodes within each element along the beam; and u i τ and u js are generalized displacement vectors. Substituting (10) into Eq. (5) and following the steps as shown in (7) gives and substituting (10) into Eq. (8) gives Finally, equating internal and external work which is a statement of elastic equilibrium in weak form in CUF notation. The term k i jτ s (e) is referred to as the element Fundamental Nucleus . Its explicit form can be found in Carrera et al. (2014Carrera et al. ( , 2011 . Fundamental nuclei may be assembled in a global stiffness matrix following the standard finite element procedure resulting in The latter equation is then solved to find the generalised displacements. Please    In CUF, the choice of F τ and M is arbitrary. In the literature, different kinds of expansion functions have been used, including Taylor, Lagrange, Legendre, exponential trigonometric and Chebyshev polynomials Carrera et al., 2013;Filippi et al., 2015;Pagani et al., 2016 ). In this work, as a novelty, we introduce and adopt Serendipity Lagrange expansions, which are described in the following sub-sections, together with more traditional models for comparison.

Cross-sectional expansion models
As mentioned in the previous section, in the Unified formulation, cross-sectional expansion functions can be chosen arbitrarily. That said, the most widely adopted expansions are based on Taylor and Lagrange polynomials. These two types of functions are used in fundamentally different ways, with profound implications on some computational and numerical aspects of the implementation.

Taylor expansion model
In Taylor expansion models, the cross-sectional displacement field at the i th Lagrange beam node is expressed with complete, Taylor polynomials containing terms of the form F τ = x n z m . For example, a second-order expansion ( N = 2 ) has constant, linear and quadratic terms as follows   . . . .   els, i.e. models with high-order kinematics, can be obtained by enriching the polynomial expansion. The reader is referred to Carrera et al. (2014) for a more detailed treatment of TE models.

Lagrange expansion model
In Lagrange expansion models, beam elements are further discretised by dividing cross-sections into a number of local subdomains as shown in Fig. 2 b. Two-dimensional Lagrange polyno-mials are used as expansion functions over the sub-domains. The order of the Lagrange polynomials spanning each sub-domain depends on the number of computational nodes therein. For instance, a 9-noded Lagrange type element ( L9 ) is spanned by quadratic expansions so that the displacement field at the i th beam node becomes   where the expansion functions, F τ = L τ , are 2D Lagrange polynomials and the terms u i τ on the right hand side are unknown variables. Unlike in TE models, these global unknowns are pure displacement components at the computational nodes defined across the sub-domains. Refined model accuracy is obtained by increasing the number of sub-domains or the number of nodes therein, or in other words, by increasing the cross-sectional mesh density. A more detailed description of LE models can be found in  .

Numerical integration over CUF elements
Sections 4.1 and 4.2 highlight one of the fundamental differences between TE models and LE models. Taylor expansions are defined to span cross-sections starting from the origin of xz -planes along the y -axis. Lagrange expansions are defined on quadrilateral sub-domains. This difference is illustrated graphically in Fig. 2 a and b.
In practical terms, the fact that Lagrange expansions are defined on cross-sectional sub-domains implies that an additional mapping is required for integrations over V (e) . To clarify, like traditional beam elements, N -noded CUF elements based on Taylor ex- , i.e. the element position in global coordinates. An identical operation is required, for elements based on Lagrange expansions, to integrate V (e) B T j C B i dV, however an additional mapping is required to link physical subdomains in cross-sectional xz -planes to the master computational 1] . A visual representation of this two-dimensional map is given in Fig. 3 .
Throughout this paper, cross-sectional sub-domains defined in ( x, z ) are mapped and interpolated using linear Lagrange polynomials L k where ξ is the vector of mapped coordinates and x k are the physical positions of the nodes of the generic quadrilateral sub-domain. As customary, by using (17) one can compute the Jacobian of the transformation, which is required for integrals oven the master domain.

Serendipity Lagrange expansion model
In TE models, it is straightforward to enrich the displacement field by choosing higher order expansions. On the other hand, in LE models, the displacement field is enriched by increasing the number of nodes in the beam cross-section. In choosing TEs over LEs one trades-off numerical stability for ease of refinement, i.e. no need for remeshing. We now introduce alternative expansion functions, based on hierarchical Serendipity Lagrange polynomials, that eliminate this duality. Adopting this expansion model, crosssections are discretised using four-noded Lagrange sub-domains. In addition, and as a novelty, the displacement field within subdomains can be enriched by increasing the order of the local Serendipity Lagrange expansion as depicted in Fig. 2 c, where the shading indicates enrichment hierarchy. The proposed expansion model is based on the hierarchical finite element shape functions as derived from Trunk (or Serendipity) polynomial spaces in Szabó and Babuška (2011) .
In als are combined and used as expansion functions for the displacement field within the computational sub-domains. Enrichment of the model kinematics can then be achieved by increasing the expansion order and/or the number of nodes in the cross-section, which will be shown to be tantamount to combining the benefits of TE and LE models, whilst also eliminating their limitations.

1D Lagrange-type polynomials
In this section, we introduce the 1D polynomials used to build the 2D SL expansions.
Let us consider the set 1D = { ξ ∈ R : −1 ≤ ξ ≤ 1 } and let N ≥ 2 be the number of equally spaced points ξ i within 1D . 2 Starting at ξ = −1 , An N th -order polynomial, p N ( ξ ), can be found such that The explicit form of this polynomial is such that, for instance, p 2 (ξ ) = (ξ + 1)(ξ − 1) , 2 By construction N will also be the order of the polynomial.
Traditional Lagrange polynomials can readily be derived from (20)    We note that the property of vanishing values on the boundary of 1D is essential to ensure continuity of the displacement field at the interfaces between cross-sectional sub-domains, which in turn allows for the formulation of hierarchical shape functions.

2D Lagrange-type polynomials
Polynomials of the family p N ( ξ ) can be used to define their N thorder 2D counterparts in 2D = { (ξ , η) ∈ R 2 : −1 ≤ ξ ≤ 1 , −1 ≤ η ≤ 1 } . These 2D polynomials are to be employed as hierarchical Lagrange-type shape functions. With this aim in mind, we need three different sets of functions, each with specific requirements: 1. A set of four first-order Lagrange polynomials. These are bilinear polynomials that take value 1 at each of the four nodes and 0 on the others. These are named polynomials of type I. 2. A set of N th -order polynomials that vanish along three sides of 2D in order to satisfy the continuity of displacements across cross-sectional sub-domains. These are named polynomials of type IIA and IIB. 3. A set of N th -order polynomials defined in the interior subset of 2D that vanish along its four sides. These are named as polynomials of type III.
Letting r = 1 , . . . , N, and s = 1 , 2 , 3 , 4 , the Serendipity expansion functions are indicated by where the subscript τ is an index defined as (1 + ξ s ξ )(1 + η s η) , where δ ij is the Kronecker delta and the argument of p r (−ξ ) and p r (−η) is negative to ensure that all L (IIA , IIB) τ polynomials of odd order are identical and separated by a 90 degree rotation; a property of shape functions required to ensure uniqueness and completeness. And finally, with n, m = 2 , 3 , . . . N, constrained by n + m = r and n + m ≤ N. Fig. 4 shows the first few polynomials L (t) τ , sorted by order, type and index τ . Henceforth, N th -order Serendipity Lagrange models are implicitly assumed to include all of the shape functions of orders 1 to N , as opposite to just order N . As an example, a model of order N = 5 contains: 1. Four bi-linear Lagrange polynomials (type I). Subscripts 1 to 4; 2. Four second-order polynomials (type II). Subscripts 5 to 8; 3. Four third-order polynomials (type II). Subscripts 9 to 12; 4. Five fourth-order polynomials (4 type II, 1 type III). Subscripts 13 to 17; 5. Six fifth-order polynomials (4 type II, 2 type III). Subscripts 18 to 23.
Similarly, cross-sectional displacements of order N = 2 , at the i th Lagrange beam node, take the form (using the notation F τ = L (t) τ ): In conclusion, the SL expansion model is beneficial in that it has characteristics of both TE and LE models, because: (a) Serendipity polynomials have the same hierarchical nature of TEs; (b) as in LE models, they are defined on sub-domains thus enabling local refinement and enhanced numerical stability via cross-sectional discretisation. A schematic representation of the trade-offs between the three expansion models is shown in Fig. 5 .

Chebyshev node distribution
CUF models can lead to inaccurate prediction of shear stresses near the boundaries. For instance, considering a cantilever beam, shear stress oscillations may be observed along the axis near the fixed support. One way to overcome this problem is to increase the number of beam elements or to use a high-order expansion in the longitudinal direction. Both choices can significantly increase the number of unknowns required for convergence. These meshdependent stress oscillations are of numerical nature and are detrimental to the objective we set out to achieve, i.e. performing detailed, yet inexpensive, stress analyses around localised features in beam-like structures. To solve this issue, we propose a simple, yet  effective, approach to redistribute the nodes with a bias towards boundaries and features. Namely, the nodes are distributed using a Chebyshev biased mesh. Chebyshev polynomials are known to give better convergence criteria and minimise Runge phenomena ( Boyd and Petschek, 2014 ). These polynomials of the first kind of order n , denoted as T n ( y ) ∈ [ −1 , 1] , are a set of orthogonal functions defined as the solutions to the Chebyshev differential equation. They are related to Legendre and Jacobi Polynomials ( Rivlin and Wayne, 1969;Arfken et al., 2013 ) and may be defined using a series expansion: Chebyshev meshes are defined using the set of zeros of (29) in [ −1 , 1] , i.e.
which can be mapped in [0, L ] as follows: As seen in Fig. 6 , we use y k to place n nodes along the length L of the beam. Consequently, the nodes are positioned more compactly towards the boundaries.

Numerical examples and discussion
Capturing 3D stress fields accurately using displacement-based weak formulations can be challenging. Since stresses and strains are obtained by differentiating the displacement field components, the stress equilibrium equations are satistied in a weak sense and not necessarily point-wise. In all numerical cases assessed here, displacements and stresses are computed at various locations along the beam. Results are compared with 3D finite element solution.

Comparison of Chebyshev and uniform node distribution
This section draws a comparison between the convergence behaviour of stress fields obtained using Chebyshev and uniform beam meshes. For this purpose, a clamped-free, square crosssection beam of length L = 1 m , height h = 0 . 1 m and width b = 0 . 1 m is considered. A load P z = −10 N is applied at the end (y = L ) , on the neutral axis, as shown in Fig. 7 . The constituent material is isotropic with Young's modulus E = 75 GPa and Poisson's ratio ν = 0 . 33 . A 3D FE analysis, performed with ANSYS , is used as a reference for validation, where the beam is discretised using 40 0 0 0 SOLID186 (3D 20-noded) elements to yield converged results.
One-dimensional CUF models, based on Taylor expansions, are used for the analyses presented in this section, as they are known to perform well with beams of square cross-section. The analyses are carried out with expansion order N = 5 and different meshes of 4-noded ( B4 ) elements with uniform and Chebyshev distributions. Ensuing nodes and respective degrees of freedom are shown in Fig. 8 , where it can be seen that the Chebyshev and uniform meshes, with 10 B4 elements, have almost half the DOFs of the uniform mesh with 20 B4 elements.
Normal stress ( σ yy ) values along the beam, at x = 0 , z = h/ 2 , are plotted in Fig. 9 a, showing that results match the ANSYS model throughout the length, except for the region near the clamped end. For further clarity, Fig. 9  As expected, results show clearly that a Chebyshev grid of 10 elements provides enhanced accuracy near the boundary than uniform meshes of 10 and 20 elements. This conclusion confirms that biased CUF meshes, refined towards regions of high stress gradients, can improve accuracy with no need for increasing the total number of nodes (and DOFs). For this reason, Chebyshev meshes are adopted for longitudinal discretisations in all of the following analyses.

Comparison between TE, LE and SL models
In this section, we compare SL expansion models with the traditional TE and LE models. First, a cantilevered, square crosssection beam is considered, as in Section 6.1 . Ten B4 elements, with a Chebyshev type distribution, are employed for the mesh in the longitudinal direction. 3D FE results are used as a reference. Analytical results, obtained with classical theories such as Euler-Bernoulli (EB) and Timoshenko beam (TB), are provided for comparison. In addition, results are also compared to Timoshenko's enhanced analytical (TB-EN) solution obtained using Airy's stress function ( Timoshenko and Goodier, 1970 ). This enhanced formulation predicts accurate transverse shear stress distribution. In chapter 11 of reference Timoshenko and Goodier (1970) , the formulation is termed as "exact". However, it is derived by enforcing certain stress components to be zero and assumes that the bending stress varies linearly along the thickness coordinate. As such, strictly speaking, the formulation is not exact, as these conditions hold true when measuring the stress distribution remote from boundary constraints. In contrast, the present formulation accounts for all stress components without any of the abovementioned assumptions and is expected to predict the stress re- sponse accurately in all regions within the structure. The following analytical expressions are employed to calculate deflection and stresses ( Timoshenko and Goodier, 1970;Filippi et al., 2015 ): (−1) n n 2 cosh (nπ ) where, as is customary, G is the shear modulus, I is the second moment of area with respect to the x axis and A is the area of the cross-section.  Transverse displacement, normal and shear stresses are evaluated at various locations as shown in Table 1 . The throughthickness variation of shear stress at the beam's midspan is plotted in Fig. 12 for SL ( N = 5 ), Taylor ( N = 5 ) and three Lagrange models with different cross-sectional meshes. Plots of the percentage error of displacement, normal and shear stress (with respect to 3D FE solution) versus DOFs are shown in Fig. 13 .
Results show that the SL model with one cross-sectional element of order N = 1 provides identical results to the Lagrange model with one L4 element. This result is expected because the models have identical kinematical descriptions. The benefits of us-ing the SL elements can be seen for expansions of order greater than one ( N > 1). SL, Taylor and Lagrange models perform similarly in terms of convergence of displacement and normal stress. Turning our attention to shear stresses, SL and Taylor expansions achieve convergence at around 20 0 0 DOFs. Conversely, as shown in Fig. 13  which upon differentiation can only provide piecewise constant or linear stress variations respectively. To demonstrate the capabilities of the proposed SL model in predicting the local variation of 3D stresses towards the clamped edges, relevant stress components are measured at several locations along the beam. In the present example, in order to capture 3D stress fields accurately, the beam's cross-section is divided into a 2 × 2 mesh of SL domains of order N = 8 . Fig. 15 shows the through-thickness variation of shear ( τ yz ) and transverse normal stress ( σ zz ) at different locations from the clamped support. In the latter region, significant localised changes in σ zz occur, which can be characterized by the presence of an inflection point. Moving away from the clamped end, boundary layer effects are less evident. Our calculations are in good agreement with 3D FE results at a significantly reduced computational cost ( ≈ 1/10 of DOFs). Similar analyses, carried out with a TE model of order N = 8 , are found to produce similar results, with some differences. For instance, Fig. 15 b shows σ zz to match the reference solution almost everywhere, except in a small region near the free surfaces, where ∂ σ zz / ∂ z is expected to vanish. Unlike the SL model, the TE expansion fails to capture this feature. This discrepancy can be explained by the fact that SLs allow not only the order of expansion to be increased, but also to discretise the cross-section. Owing to these ca-  JID: SAS [m5G;February 28, 2018;7:59 ] pabilities, boundary effects in the stress profiles can be more readily captured. In a TE setting, the only way to improve the prediction of transverse normal stresses along the beam's free surface is to increase the expansion order. However, this leads to numerical instabilities, which may be measured by computing the conditioning number ( r c ) of the ensuing stiffness matrix ( Kreyszig, 2011 ). Fig. 16 a is a plot of 1/ r c , reciprocal of the conditioning number, versus, N , the expansion order of SL and TE models with one cross-sectional element. From the figure, we observe that, for increasing N , the stiffness matrix of TE models becomes ill-conditioned (i.e. r c diverges). Conversely, the conditioning properties of SL models are almost independent from the expansion order. This is shown to be the case also for LE models, proving that cross-sectional discretisation improves numerical stability.

T-section beam
In order to show the enhanced capabilities of SL expansions, in comparison with TEs and LEs, an additional beam of more complex geometry is examined. Specifically, we consider the T-section beam shown in Fig. 17 . Material properties are the same as in the previous example. The beam is clamped at one end and loaded with a concentrated force, P z = −10 N , at the other end. The analysis is performed with Taylor, Lagrange and SL models. Converged 3D finite element results from ANSYS , computed by discretising the structure with 554,036 SOLID186 elements, are taken as a reference for comparison. Displacement fields, as well as, normal and shear stresses are evaluated at several locations and a convergence analysis is performed by varying the order of Taylor and SL expansion and by refining the cross-sectional mesh for LE. For an accurate estimation of the stress field at the intersection between flange and web, several cross-sectional LE and SL meshes have been trialled. Resulting discretisations are shown in Fig. 18 , where it can be seen that local refinement is required in the regions with high stress gradients. For the LE mesh, convergence is achieved with 488 L9 elements; In comparison, the SL model necessitates some 66 SL8 elements. Fig. 16 b confirms that, also in this case, TE models lose numerical stability for increasing N , which limits our analyses to order 9. In contrast, LE and SL are found to be numerically stable again.
Elastic field results are reported in Table 2 . As expected, Taylor models produce accurate and converged displacement and normal stress values, but fail to represent shear stresses to an acceptable degree of precision. Lagrange and SL models are numerically stable, as such they are able to capture the response of the structure better than TEs, particularly localised stresses concentrations. The reason for this difference is that Lagrange and SL expansions rely on local discretisations at cross-sectional level, whereas Taylor expansions are constructed with displacement shape functions spanning the entire cross-section from the beam reference axis, which affects the conditioning number negatively, thus preventing indefinite refinement.
In the remainder of this section, particular attention is given to τ yz , which, as indicated by Table 2 , is the most problematic field variable to be modelled accurately. Figs. 19 and 20 show the variations of shear stress at the beam's mid span, respectively, along the flange and through the web at x = 0 . 07 and x = 0 . 06 . In addition, the models are interrogated throughout the beam's length.
Values of τ yz through z , at 2%, 5% and 50% of the span from the clamped end, are reported in Figs. 21 and 22 . The latter, shows the shear stress distribution along the T-section's web. At y = 50% L, such distributions can be calculated analytically using Jourawski's formula ( Gere and Goodno, 2011 ). This is done to highlight an example of the intrinsic limitations that may affect sim plified models. Specifically, it is observed that the formula deviates from the numerical results, proceeding from the top of the section towards the flange. This result is as expected due to the assumptions in Jourawski's model. In summary, shear stresses from the LE, SL and 3D finite element solutions match almost exactly and can capture localised features in the 3D stress field.
Finally, for further appraisal of SL discretisations, 3D stress profiles across full cross-sections are compared to the reference AN-SYS solution through contour plots of transverse shear and normal stresses at various span-wise locations. These positions are shown in Figs. 23-25 . Overall agreement is excellent, except at the corner between the flange and web, which theoretically is a singular point. No model is accurate in capturing stresses exactly in this location.
In conclusion, from the results presented in this section it is evident that the SL models are capable of accurate stress predictions with considerably less DOFs than 3D FE, which is a proxy for computational cost. From a numerical standpoint, SL and Lagrange models behave identically. This result allows either of the two models to be used with confidence. SL meshes, however, give an extra advantage, because, unlike LE meshes, they facilitate element-wise hierarchical refinement thereby reducing the need for cross-sectional remeshing.

Conclusions
We have aimed to capture 3D stress fields accurately using 1D models with greater computational efficiency than 3D finite element analyses. The Serendipity Lagrange expansion model is introduced within the framework of Carrera's Unified Formulation. Static analyses of square and T-section beams have been carried out to both challenge and exemplify the merits of the proposed approach. The model is benchmarked against traditional Taylor and Lagrange expansions, 3D finite element solutions as well as analytical formulae (where available). The main findings from the present study can be summarised as follows: 1. The effect of collocating beam nodes using a Chebyshev biased mesh has been studied. The mesh was refined in the regions where stress fields are expected to change rapidly. It has been observed that, by employing this node distribution, accurate results can be obtained near constraints, without the need to increase the total number of beam nodes. This type of discretisation also precludes spurious oscillations in the solutions, previously observed in CUF models. 2. For the numerical cases assessed, the Serendipity Lagrange expansion model retains benefits of both the Lagrange model (cross-sectional discretization) and the Taylor model (hierarchical approximations), eliminating their disadvantages, as described in the following points. 3. In order to capture the response of beam-like structures accurately, high-order models may be required. For Taylor models, as the order of expansion increases, the conditioning number of the stiffness matrix decreases exponentially. This problem makes the system ill-conditioned and numerically unstable. Serendipity Lagrange expansions overcome this limitation and are therefore suitable for analysing beams with complex cross-sections. 4. Similarly to Lagrange expansion models, the Serendipity Lagrange ones allow for cross-sectional discretisation. This feature, together with the hierarchical nature of the local expansions, makes Serendipity Lagrange elements particularly suited for capturing localised stress fields near boundaries, discontinuities and points of load application, unlike the Taylor expansion model.  JID: SAS [m5G;February 28, 2018;7:59 ] 5. Cross-sections are also discretised in the Lagrange model, however model building is cumbersome because remeshing is the only way to improve accuracy.

ARTICLE IN PRESS
The proposed Serendipity Lagrange expansion models proved to be an efficient and effective means for computing 3D stress fields for solid and thin-walled isotropic beam structures. Future work will focus on composite structures and a comparison with other Unified Formulation models, such as Legendre based models ( Pagani et al., 2016 ), that show similar advanced capability.