Stress recovery of laminated non-prismatic beams under layerwise traction and body forces

Emerging manufacturing technologies, including 3D printing and additive layer manufacturing, offer scope for making slender heterogeneous structures with complex geometry. Modern applications include tapered sandwich beams employed in the aeronautical industry, wind turbine blades and concrete beams used in construction. It is noteworthy that state-of-the-art closed form solutions for stresses are often excessively simple to be representative of real laminated tapered beams. For example, centroidal variation with respect to the neutral axis is neglected, and the transverse direct stress component is disregarded. Also, non-classical terms arise due to interactions between stiffness and external load distributions. Another drawback is that the external load is assumed to react uniformly through the cross-section in classical beam formulations, which is an inaccurate assumption for slender structures loaded on only a sub-section of the entire cross-section. To address these limitations, a simple and efficient yet accurate analytical stress recovery method is presented for laminated non-prismatic beams with arbitrary cross-sectional shapes under layerwise body forces and traction loads. Moreover, closed-form solutions are deduced for rectangular cross-sections. The proposed method invokes Cauchy stress equilibrium followed by implementing appropriate interfacial boundary conditions. The main novelties comprise the 2D transverse stress field recovery considering centroidal variation with respect to the neutral axis, application of layerwise external loads, and consideration of effects where stiffness and external load distributions differ. A state of plane stress under small linear-elastic strains is assumed, for cases where beam thickness taper is restricted to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$15^{\circ }$$\end{document}15∘. The model is validated by comparison with finite element analysis and relevant analytical formulations.


A(x)
Arbitrary cross-sectional area with undefined boundaries ( m 2 ) A(x) Arbitrary cross-sectional area with defined lower boundary ( m 2 ) A(x) Generalization of the cross-sectional area multiplied by Young's modulus ( N) c (x) z-direction distance between the reference system and centroidal position ( m) E Young's modulus (Pa) f (k) (x) Functions describing the lower surface of the kth layer ( m) F (k) i (x, z) Body force acting on the kth layer ( Coefficient accounting for traction forces oblique to the coordinate system ( m)

I(x)
Generalization of the second moment of area multiplied by Young's modulus ( Nm 2 ) L Beam length ( m) m(x) Bending moment due to eccentric axial loading ( Nm∕ m)

M(x)
Bending moment resultant ( Nm) n

S(x)
Generalization of the first moment of area multiplied by Young's modulus ( Nm) t (k) i (x) Surface force acting on the lower interface of the kth layer ( N∕ m) External load resultant ( N∕ m) Beam width ( m) Greek characters Γ j (x) Terms of the equilibrium of internal forces and their derivatives Δ(.) Discontinuity operator applied to interfaces Discrepancy between results ( %) (k) i (x, z) Stiffness-related coefficients necessary for zz (k) i (x, z) Stiffness-related coefficients necessary for xz (k) i (x) Constants variables originated from undefined integration ( Pa) (k) i (x, z) Stiffness-related coefficients necessary for zz Poisson's ratio (-) ij Stress tensor ( Pa) (x) Ratio between first moment of area in relation to the centroid position and Ĩ ( m∕ N) Ω (k) Domain of the kth layer ( − )

Introduction
Laminated non-prismatic, or multilayered tapered, beams are variable cross-section slender structures comprising multiple materials stacked in layers in a specific order (Balduzzi et al. 2017a(Balduzzi et al. , 2016. To manufacture tapered beams, modern technologies such as 3D printing and automated welding can shape beam-like structures of complex geometry. As a result, more cost-efficient structural elements can be produced by optimizing stiffness distributions (Timoshenko and Young 1945;Bai et al. 2018;Vo et al. 2021) while meeting design criteria. Regarding the manufacture of tapered composite laminates, two methods are employed: ply drop-off and variable prepreg tape spreading. The ply drop-off technique has gained attention during the late twentieth century, and involves terminating constant thickness plies with associated resin pocket tapering (Curry et al. 1992; Mukherjee and Varughese 2001;Varughese and Mukherjee 1997;Sudhagar et al. 2017). However, eliminating layers results in stress concentrations due to cutting of load bearing fibers (Varughese and Mukherjee 1997;Her 2002). Conversely, the variable tape spreading technique developed by Clancy et al. (2020) gradually introduces tapering to laminates composed of carbon fiber pre-impregnated with thermoplastic resin. This technique utilizes laserassisted tape placement (LATP) in conjunction with a spreading device before consolidation takes place, thus avoiding stress concentrations by doing so. As technologies emerge to manufacture heterogeneous tapered beams, sophisticated analytical methods are required to comprehend the mechanical behavior of innovative structures. Hence this investigation focuses on laminated non-prismatic beams with variable thickness layers. In addition to being lighter, non-prismatic beams can be manufactured to meet specific geometrical shapes driven by design requirements. Well-known examples include bridges, helicopter and wind turbine blades and aircraft wings (Hodges 1990;Balduzzi et al. 2017d;Vinod et al. 2007;Vo et al. 2021). Nevertheless, predicting tapered beam behavior through 1D beam formulations poses some challenges. For example, the state of stress on beam surfaces is governed by Cauchy's traction-stress relation, resulting in non-vanishing transverse stresses (Balduzzi et al. 2016(Balduzzi et al. , 2017a(Balduzzi et al. , b, c, d, 2019Mercuri et al. 2020;Hodges et al. 2008Hodges et al. , 2011Vilar et al. 2021a, b). For heterogeneous media, Cauchy's traction-stress equation must also be satisfied on interface surfaces, which causes transverse stress discontinuities (Balduzzi et al. 2017a). Furthermore, beam eccentricity, which describes the offset of centroidal position from the neutral axis, plays an essential role in deformation and stress distribution (Balduzzi et al. 2016), but is often neglected in non-prismatic beam formulations. Another drawback is that the Euler-Bernoulli beam theory assumes the external load is reacted uniformly through the cross-section, which is referred to in this work as generalizing the external load. This hypothesis, however, is not appropriate to cases where the external load is applied to a portion of the cross-sectional area. Practical examples include aircraft wings, wind turbine and helicopter rotor blades under fluid pressure and shear forces; arched beams in bridges under vehicular traction forces and tensile stresses in tendons of prestressed concrete (Vilar et al. 2021b). As such, an accurate and effective analytical method for stresses is necessary to exploit the full potential of laminated non-prismatic beams.
Various analytical solutions have been proposed to predict the stress field of homogeneous non-prismatic beams. Bleich (1932) derived the shear stress of linearly tapered beams based on a generalization of the parabolic shear stress distribution of prismatic beams, also known as Jourawski's formula (Jourawski 1856). Subsequently, many investigations proposed shear stress recovery procedures to infinite wedge elements via Theory of Elasticity approaches (Michell 1900;Carothers 1914;Timoshenko and Goodier 1951), and later to truncated linearly tapered beams (Krahula 1975), including that of I-sections (Blodgett 1966;Vu-Quoc and Léger 1992;Trahair and Ansourian 2016;Balduzzi et al. 2017b). Boley (1963), using Airy's stress function, examined the limits of a linear longitudinal stress distribution in non-prismatic beams under pure bending, and concluded that the error associated with Navier's hypothesis escalates with the taper angle, such that a 10 • -taper results in 7.5% of error in predicting the longitudinal stress. Romano (1996) offered a solution to shear stress for linearly tapered beams using Timoshenko beam theory. Hodges et al. (2008Hodges et al. ( , 2011 used the variational asymptotic method to model linearly tapered beams under extension, bending, and transverse loads. Furthermore, Beltempo et al. (2015) employed an energy-based approach in conjunction with the Hellinger-Reissner Principle to predict the mechanical behavior of non-prismatic beams. Later, Mercuri et al. (2020) used the Hellinger-Reissner Principle to propose solutions for the full 2D stress field. Balduzzi et al. (2016) reconstructed the shear stress of tapered beams by superposing Jourawski's formula to a linear function that satisfies traction-free boundary equilibrium on oblique surfaces. (Zhou et al. 2016(Zhou et al. , 2020b predicted the shear stress of asymmetric linearly tapered beams with box girder and rectangular crosssections, respectively. Moreover, Taglialegne (2018) and Bertolini et al. (2019) derived the stress field of symmetric linearly tapered thin-walled beams. Subsequently, Bertolini and Taglialegne (2020) proposed a similar methodology to include taper in the width direction. Chockalingam et al. (2020) derived the shear stress of I-beams tapered in the thickness and width directions. Also, Vilar et al. (2021aVilar et al. ( , 2021b recovered the 2D stress field of non-prismatic beams subject to generalized and partial cross-sectional loads, respectively. In contrast to homogeneous media, the fewer investigations concerning laminated tapered beams leave scope for important limitations to be addressed. For example, Balduzzi et al. (2017a) proposed a solution for shear stress and displacements in laminated tapered beams under plane stress. However, beam eccentricity and the spanwise distributed bending moment arising from eccentric axial loading are disregarded in the stress recovery procedure, although both are included in the equilibrium relations of generalized forces and moments. Ai and Weaver (2017) modeled tapered sandwich beams with a functionally graded core based on a layerwise displacement field approximated via the Ritz method, yet the beam shape is limited to being symmetric with respect to the longitudinal axis. Zhou et al. (2020aZhou et al. ( , 2021 proposed solutions for shear stress in non-prismatic beams with a corrugated steel web considering the Resal effect, i.e., the reduction of the effective shear force due to compressive forces acting on the tapered flange, and shear lag analysis, respectively. However, in both investigations, cross-sections were restricted to thinwalled members. In general, classical beam models generalize the external load through the entire crosssection. This hypothesis, although not best representing surface forces and partial cross-sectional loads, is widely accepted in beam modeling approaches (Masjedi et al. 2019;Weaver 2020a, b, 2022;Masjedi et al. 2021;Doeva et al. 2020aDoeva et al. , b, 2021Doeva et al. , 2022Balduzzi et al. 2019;Bertolini et al. 2019;Bertolini and Taglialegne 2020;Chockalingam et al. 2020Chockalingam et al. , 2021Vilar et al. 2021a). Another issue not addressed in analytical laminated beam formulations is for cases where the stiffness distribution is different from the transverse and axial load distributions, which gives rise to non-classical transverse stresses. More specifically, derivatives of the axial force, i.e., the generalized longitudinal stress, are required to predict the transverse stresses: the first derivative for the shear component and up to the second for the transverse direct stress. Additionally, the first derivative of the shear force is necessary to define the transverse direct stress, noting that the contribution of this term differs from that of prismatic beams.
To remedy the state-of-the-art limitations, we provide an analytical stress recovery procedure to laminated non-prismatic beams subject to layerwise body forces and traction loads. Beam eccentricity is accounted for in the equilibrium equations and stressfield derivation. Non-classical terms arise from crosssectional interactions between stiffness distribution, and the transverse and longitudinal external load distributions. The cross-sectional shape is arbitrary but has simplifications, including cross-sections with no cut-outs and interfacial surface projections are parallel to the global reference system. To give insight into the potential of the developed formulation, closedform solutions are deduced for the case of a rectangular cross-section.
The proposed method can be applied to untwisted laminated beams tapered in the thickness direction with perfectly bonded interfaces. A state of plane stress is assumed, coupled with linear-elastic material undergoing small strains. A linear longitudinal stress distribution with discontinuities at interfaces is adopted. Subsequently, the transverse stress field is recovered following Jourawski's formulation adapted to heterogeneous non-prismatic beams. In the context of analytical methods for laminated non-prismatic beams, the main novelties of the proposed approach include: (1) the recovery of the shear stress distribution considering beam eccentricity; (2) consideration of non-classical terms originating from the mismatch between stiffness and external load cross-sectional distributions; (3) derivation of the transverse direct stress component and (4) application of layerwise loading and traction forces.
The outline of this work is as follows: Sect. 2 provides the hypotheses and limitations as well as the equilibrium relations. Section 3 recovers the transverse stress field. Section 4 compares results of the developed theory with relevant analytical formulations and finite element analyses. Finally, Sect. 5 summarizes the novelties of this study and suggests further developments.

Problem idealization
This section introduces the underlying principles of the proposed theory. Herein the Cartesian coordinate system is adopted such that the x-axis is parallel to the beam axis and the y-and z-axes are parallel to the beam width and height (or thickness), respectively.
The object of study is a laminated non-prismatic beam comprising n layers as depicted in Fig. 1. The xz-plane is under plane stress undergoing small displacements and linear strains. Also, the crosssection of the beam remains plane after deformation but not necessarily perpendicular to the neutral axis, thus consistent with the First Order Shear Deformation Theory. It is assumed that the beam length L is significantly greater than the depth h(x) and width w(x, z) to meet the slenderness requirements of beam theory.
To be consistent with the hypotheses assumed, the following limitations are imposed to the cross-sectional shape: (1) symmetry with respect to the z-axis is adopted to avoid cross-section warping effects; (2) the ratio w/h should be sufficiently small to reduce shear stress in the yz-plane (Timoshenko and Goodier 1951), noting that an optimized ratio depends on specific cross-sectional shapes; (3) w∕ x ≈ 0 to be consistent with the plane stress hypothesis; (4) interfacial boundaries are perfectly bonded and their projections to the cross-section are parallel to the xy-plane; (5) cut-outs are disregarded to avoid stress concentrations.
Introducing (n + 1) interface surfaces f (k) , that represent the z-coordinate of the kth layer's lower interface, note that f (n+1) and f (1) correspond to the upper and lower beam surfaces, respectively (dependency on the x-coordinate is omitted). The taper angle, i.e., the angle formed by the x-axis and the plane tangent to the beam surfaces, is smaller than 15 • throughout the longitudinal axis for reasons explained in Sect. 3. The adopted geometrical description implies that cross-sections are non-orthogonal to the beam's neutral axis. Instead, orthogonality is preserved in the Cartesian coordinate system, noting similar methodology is adopted in related works (Balduzzi et al. 2016;Gimena et al. 2008;Rajagopal and Hodges 2014;Balduzzi et al. 2017a, c, d;Mercuri et al. 2020;Vilar et al. 2021a). Consequently, the cross-sectional layer area A (k) and total cross-sectional area A are given by The domain of each layer Ω (k) is then defined by

Material and mechanical properties
Laminated non-prismatic beams comprise multiple homogeneous layers stacked in a specific order. As such, the Young's modulus of the entire beam "E(x, z)" can be expressed as a piecewise function (2.1) Since the material properties vary through the crosssection, it is appropriate to include Young's modulus in the generalization of stiffness properties. For example, the longitudinal stiffness is given by Noting that Eq. (2.4) agrees with the classical definition of longitudinal stiffness of homogeneous beams for the case where n = 1 . Analogous to Ã (x) , the first moment of stiffness S (x) is defined as Consequently, the neutral axis c(x) reads Clearly, the stacking sequence plays an important role in defining the neutral axis position. Finally, the bending stiffness in relation to the neutral axis yields The expressions introduced in this section are necessary to adapt essential concepts of homogeneous beam theories to heterogeneous media. However, Eqs. (2.4), (2.5) and (2.7) should not be interpreted as work conjugate of generalized stresses and strains but as stiffness-related variables that simplify the stress recovery procedure. It is noteworthy that curvature, axial, and shear deformations are coupled with all generalized stresses (Vu-Quoc and Léger 1992;Balduzzi et al. 2016Balduzzi et al. , 2017aMercuri et al. 2020); hence constitutive relations between generalized strains and stresses are not readily extracted from application of compatibility relations to internal force definitions, as is common practice in prismatic beam theories. To overcome this impasse, previous works used the Complementary Virtual Work Principle (Vu-Quoc and Léger 1992), the Hellinger-Reissner Principle (Beltempo et al. 2015;Auricchio et al. 2015;Mercuri et al. 2020), and stress potential approaches (Balduzzi et al. 2016(Balduzzi et al. , 2017b. Since the proposed methodology is based on Newton's Third Law, it only requires equilibrium relations (derived in the next section) and a definition for the longitudinal stress to recover the transverse stress components, which are expressed in terms of Eqs. (2.4)-(2.7).

Equilibrium relations
The beam element is subject to layerwise body forces , uniformly distributed through a layer and variable spanwise, and to traction forces t (k) i applied to f (k) , as Fig. 2 depicts. The body force resultants in the longitudinal and transverse directions, namely T x and T z , as well as the bending moment m originating from eccentric axial loading, are given by The stress resultants are known variables but not necessarily evaluated from statically determinate models. The generalized longitudinal and shear stresses N(x) and V(x), respectively, and the generalized bending moment M(x) are given by the following expressions where xx and xz are the longitudinal and shear stress components, respectively. From now on, the generalized bending moment M(x) is referred to as bending moment while the terminology bending load is adopted to refer to the bending moment m(x) caused by eccentric axial forces. Note that, even if the external load is reacted uniformly distributed through the cross-section, and no traction is prescribed, m(x) ≠ 0 because the neutral axis position depends on the stiffness distribution, which is irrespective of the external load idealization.
Consider a segment of a laminated non-prismatic beam in equilibrium as illustrated in Fig. 2, where internal forces vary linearly along the infinitesimal length.
The Timoshenko equilibrium expressions are adapted to non-prismatic beam shapes by accounting for the neutral axis variation as follows The cross-sectional stiffness distribution does not affect the equilibrium of forces, but it influences the equilibrium of bending moments. For single layered beams, Eq. (2.10) simplifies to the equilibrium equations of homogeneous tapered beams (Vilar et al. 2021a;Balduzzi et al. 2016Balduzzi et al. , 2017c.

Interfacial boundary equilibrium
Consideration of equilibrium of forces requires specific relations between stress components at interfaces. One method to derive the interlayer boundary condition is to adapt the Cauchy stress equilibrium by selecting an appropriate representative unit cell, such as that illustrated in Fig. 3.
According to classic Cauchy stress equilibrium, an interior point of a 2D structure is represented by a rectangular shape of two independent infinitesimal dimensions. However, to best represent the interfacial boundary, the adopted unit cell dimensions are related by f �(k) . Also, classic Cauchy stress equilibrium assumes the stress field varies linearly within an infinitesimal point, but this hypothesis does not hold for regions of discontinuous material properties. Instead, the stresses acting on the unit cell faces are independent variables. Invoking the equilibrium of forces in the longitudinal and transverse directions shown in Fig. 3 yields where zz is the transverse direct stress component and H (k) f the hypotenuse of Fig. 3, and the Δ(.) (k) oper- ator indicates a discontinuity of its argument evaluated at interface surface f (k) The second term of the right-hand side of Eq. (2.12a) vanishes for the first layer. It is noteworthy that the body force contribution to the interlayer force equilibrium is infinitesimal compared to the stress and traction resultants, and that one can substitute Eq. (2.11a) into Eq. (2.11b) to obtain the transverse interface boundary condition in terms of longitudinal stresses. Equation (2.11) suggests that the transverse stresses are continuously distributed through the depth of the beam for traction-free prismatic beams, as reported in related studies (Bareisis 2006;Balduzzi et al. 2017a;Reddy 2003;Pagano 1978). Conversely, the transverse stress components are discontinuously distributed through the cross-section in non-prismatic beams due to interface surfaces being oblique to the global reference system (Balduzzi et al. 2017a).
It is relevant to highlight that the equilibrium derived in Eq. (2.11) is valid for an interior point of the crosssection over the interface boundary. Hence the tractionfree boundary condition on the cross-section contour is not satisfied since the plane stress hypothesis is not valid in this region.

Stress field recovery
This section introduces analytical expressions for the transverse stress field of laminated non-prismatic beams of general cross-sectional shapes. Closed-form solutions can be obtained from the analytical formulation by deriving the stiffness and geometric properties for specific cross-sectional shapes, which are given in the Appendix for the specific case of rectangular cross-sections.
The proposed method integrates the Cauchy stress equilibrium of an interior point in an arbitrary lamina, followed by imposition of interfacial boundary conditions derived in Sect. 2.3. As a result of the linear strain assumption coupled with a state of plane stress and material idealization, the longitudinal stress yields the well-known expression (Zenkert 1997) Noting that the taper angle is limited to 15 • to reduce error of Eq. (3.1) (Boley 1963;Balduzzi et al. 2016;Vilar et al. 2021a).

Shear stress
Cauchy stress equilibrium applied to an interior point of an arbitrary layer yields where Ā is an arbitrary slice of a cross-sectional area with an arbitrary z value for the upper boundary and undefined lower boundary. Information on the lower boundary coordinate is recovered by imposition of interface boundary conditions. To simplify operations on Eq. (3.2), the vector Γ containing terms of the equilibrium relations and their derivatives is introduced where () � denotes differentiation.
Substituting for (k) xx from Eq. (3.1) and the generalized external load and internal forces defined in Eqs.
(2.8a) and (2.10) leads to the shear stress expressed as a linear combination as follows 1 is a measure of stress, constant in the yzplane, as a result of the indefinite integration performed in Eq. (3.2). Expressions for (k) are given by Ā noting that the last term of the right hand side of (3.2) is expressed in terms of generalized longitudinal stress in the expressions for (k) N � . The introduced variable is defined as i are defined as longitudinal Cauchy stress coefficients. These stiffness coefficients relate the internal forces to the shear stress distribution (e.g.

(k)
M relates bending moment to shear stress). As such, can be understood as intra-layer functions that recover the shear stress distribution from terms of the equilibrium relations and their derivatives.
The xz dependency on all internal forces is referred to as the non-triviality of stress field in previous works (Vu-Quoc and Léger 1992;Hodges et al. 2011;Balduzzi et al. 2016Balduzzi et al. , 2017aBruhns 2003;Vilar et al. 2021a). In addition to the internal forces and the bending load, N � (x) contributes to the shear stress profile of laminated beams. In homogeneous media, N ′ vanishes if the external load is generalized through the cross-section (Taglialegne 2018;Bertolini et al. 2019;Bertolini and Taglialegne 2020;Vilar et al. 2021a Conversely, this balance does not occur in heterogeneous beams since the ratio E (k) ∕Ã could take any value from zero to one. Hence, in general, (k) N � does not vanish for laminated beams, and its effects should be considered for practical applications.
Apart from (k) V , all longitudinal Cauchy stress coefficients carry information on cross-sectional variation and beam eccentricity. For single layered prismatic beams, Eq. (3.4) reduces to Jourawski's formulation. This observation suggests that the classic solution for prismatic beams is part of the solution for tapered shapes (Vilar et al. 2021a, b).
The constant (k) 1 is evaluated by setting z = f (k) in Eq. (3.4), followed by imposition of the longitudinal interfacial condition stated in Eq. (2.11a) The expression for (k) 1 shown in Eq. (3.7) can be separated into two different groups: (1) the first term balances the Eq. (3.2) residual originating from integration of Cauchy stress equilibrium; (2) the remaining variables enforce longitudinal interfacial equilibrium.
The constant (k) 1 requires information on the lower adjacent layer to define the shear stress. For instance, xz profile can be determined only if (1) 1 and (1) xz have already been evaluated. As such, after performing successive substitutions from the kth to the bottom layer, the shear stress profile is expressed as follows Note that the successive substitution process is considered by the summation operator from dummy index l varying from one to the kth layer.
The shear stress given by Eq. (3.8) is novel and more accurate than related analytical expressions reported in up-to-date literature (Balduzzi et al. 2017a). The key novelties regarding xz derived include consideration of: the longitudinal traction and body forces, both defined for each layer; the beam eccentricity; the first derivative of the generalized longitudinal stress.

Transverse direct stress
Analogous to the methodology derived in Sect. 3.1, Cauchy stress equilibrium is invoked in the transverse direction Substituting the expression for (k) xz given by Eq. (3.4) and the generalized external load and internal forces defined in Eqs. (2.8a) and (2.10), yields where (k) 2 is a constant value of stress in the crosssection as a result of indefinite integration of Eq. (3.9).
Expressions for (k) j are given by Note that the dummy index in the summation operator of Eq. (3.10) varies from j = (1, 2, … , 7) , representing a range of values that differs from that in Eq. (3.4) because the transverse direct stress requires all terms of Γ (see Eq. (3.3)). As such, the transverse direct stress depends on all internal forces, the bending load, and their derivatives (up to the second derivative for the case of the generalized longitudinal stress). The constant (k) 2 is deduced by invoking the transverse interlaminar boundary condition expressed in Eq. (2.11b) Note that z = f (k) has been omitted for Ā and j . Analogous to (k) 1 , variable (k) 2 can be characterized by different contributions: (3.10) (1) the first term balances the Eq. (3.9) residual that arises due to indefinite integration; (2) the last two terms ensure transverse interlayer equilibrium; (3) the term involving (k) 1 arises from mathematical operations on Eq. (3.9) regarding longitudinal interlayer equilibrium.
Following an analogous procedure of successive substitutions explained in Sect. 3.1, an analytical expression for the transverse direct stress is obtained where the dependence of (k) j on (x, z) and the Δ(.) operators on (x, f (k) ) have been omitted. The introduced variable Â (l) is defined by The transverse direct stress expressed by Eq. (3.13) is new and applies to laminated tapered beams under layerwise body loads and traction forces, considering the simplifications discussed in Sect. 2.

Numerical results
This section solves numerical examples to discuss and validate the proposed methodology. To benchmark the developed theory, results are compared to finite element analysis solutions, which are referred to as "FEA" in the next section. The current formulation is labeled as "mod", whose longitudinal and shear stress distributions are compared with Balduzzi et al.'s (2017a) solution, denoted as "Bal". The differences between Balduzzi et al.'s (2017a) shear stress to that derived in Eq. (3.8) concern the neglect of beam eccentricity and bending load, as well as the assumption that the external load is generalized through the cross-section. To the best of the authors' knowledge, no other analytical solution for zz has been derived for laminated tapered beams.
The following beam models are solved considering a rectangular cross-section: (1) a symmetrically tapered sandwich cantilever subject to a transverse concentrated force at the tip; (2) a simply supported double tapered beam under a transverse body force applied to a thin layer; (3) an asymmetric cantilever beam with the upper surface flat under: (3.1) a longitudinal uniform body force; (3.2) constant traction and pressure loads applied to the upper surface.
The terminology "case (.)" is used to refer to the beam models in this section. A mesh refinement study was performed to ensure the finite element analyses converged, noting that the number of elements between mesh refinements is multiplied by a factor of 1.5-2.5. Relative errors defined as i = ( (i) zz − (i−1) zz )∕ (i) zz were calculated between mesh refinements at a specific point of the beam domain (index i refers to the ith mesh refinement). Finite element analysis results were considered to be converged when i ≤ 0.02% . Finally, the relative error of the proposed methodology in relation to the converged finite element mesh is given by (mod concerns the developed model). Note that mod should not be interpreted as a measure of accuracy of the proposed model because it is evaluated at a specific point. The 8-node quadrilateral plane stress element (CPS8R in Abaqus terminology) was used for cases (1) and (3) while the 6-node triangular plane stress element (CPS6R) was used for case (2). More information on the mesh structure is given in the respective model discussion. In each case, the processing time has been calculated to compare performance, noting that the proposed formulation was coded in Matlab2019a. The computer used for both analyses was a model i7-6800 CPU @ 3.40GHz with 32 GB RAM installed.

Linearly symmetric
Consider the laminated cantilever symmetric beam illustrated in Fig. 4, of length L = 10 m and width w = 1 m . Due to its symmetry in the x-direction, this model has zero beam eccentricity. As such, the aim of case (1) is to test the proposed formulation for the specific case of c � (x) = 0.
This beam has outer layers of Young's modulus E (1) = E (3) = 800 GPa separated by a softer core E (2) = 50 GPa . The following interface surfaces are set for this model (in mm) The thickness of the core is double that of the outer layers, and the beam depth at the fixed end is h(0) = 1.25 m , while at h(L) = 0.625 m , resulting in a constant 1.79 • taper angle.
The beam is subject to a concentrated force at the tip such that T z (L) = −1 kN , resulting in the internal forces at cross-sections x = 0.25L and x = 0.90L given in Table 1. A typical mesh layout is shown in Fig. 5; the number of elements through the thickness is constant, resulting in an average aspect ratio of 1.07 for the converged mesh. Values for zz are evaluated at the top surface of the cross-section at x = 0.25L for the mesh convergence analysis given in Table 2. The adopted Poisson's ratio was = 0.25 for all layers. The processing time of the finite element analysis was 5.6 s using Abaqus software against 0.053 s from the proposed method. (4.1) x − 375  Table 1 Case (1) internal forces The longitudinal stress field is shown in Fig. 6. Notably, all methodologies result in the same distribution. The magnitude of xx at the outer layers is significantly greater than at the core due to their larger Young's modulus. Moreover, the discontinuity at f (2) and f (3) represent the sudden change in elastic properties on the outer-layer/core interfaces. Figure 7 depicts the results for xz . Remarkably, all solutions provide the same parabolic distribution for both cross-sections. A larger discrepancy is observed between the finite element analysis and the current formulation for x = 0.90L due to localized effects closer to the applied load. The maximum absolute shear stress magnitude is located at the beam surfaces for x = 0.25L and the outer-layer/core interfaces for x = 0.90L . This observation contradicts Jourawski's solution for prismatic beams and highlights the nontriviality of the shear stress distributions. Additionally, the shear stress discontinuity at interfaces confirms the interface boundary equilibrium discussed in Sect. 2.3.
The transverse direct stress is illustrated in Fig. 8, which shows that the solutions from the current formulation and that of finite element analysis match each other for both cross-sections. A small discrepancy between the two formulations is observed at x = 0.90L due to localized Saint Venant's effects. Note that, the maximum zz occurs at the same location as maximum xz , i.e. at the beam surfaces   for x = 0.25L and at the outer-layer/core interfaces for x = 0.90L . It is relevant to highlight that the discontinuities at the interfacial surfaces for the case of zz are smaller than that for xz because Δ zz = f � Δ xz and the small taper angle assumption drives f � (.) ≤ 0.268 . Consequently, more significant stress discontinuities are expected in zz rather than in xz for the case of f � (.) > 1 , which corresponds to a local taper angle larger than 45 • , noting that this magnitude is an unlikely circumstance for most beam applications.

Double-taper
The section solves the beam model illustrated in Fig. 9. Case (2) consists of a simply supported double-tapered beam comprising two different materials: a 30 mm stiffer layer of Young's modulus E (2) = 210 GPa embedded in softer layers with E (1) = E (3) = 25 GPa . Moreover, the beam has L = 10 m and w = 1 m. The depth of the beam varies linearly over its span such that h(0) = h(L) = 0.5 m and h(0.5L) = 1 m . Consequently, the taper angle is 9.09 • for the left half of the beam and −9.09 • for the right half. The interfacial surfaces of case (2) are given by the following expressions (in mm) This numerical example highlights a practical problem in engineering employed in building construction. Figure 9 could represent an application of prestressed concrete, which consists of highly-tensioned robust steel bars with a curved shape (Gilbert and Mickleborough 1990). According to concrete design codes (ACI Committee and International Organization for Standardization 2008; EN 2004), steel bar stiffening effects can be disregarded for design purposes. To give an insight into the error associated with this simplification, results for the specific case where the Young's modulus E = 25 GPa is constant through the thickness are provided (referred to as "mod (E (2) = E (1) )").
The beam is subject to a distributed body force F (2) z = 1.111(10 4 kN∕m 3 ) applied to the stiffest layer such that the transverse resultant over a cross-section is constant spanwise T z = 333.33 kN∕m . Table 3 reports the internal forces and their derivatives for x = 0.25L and x = 0.60L . It is relevant to mention that the current formulation is not able to predict the stress field at the midspan due to stress-channeling effects (Everstine and Pipkin 1971). Instead a higherorder theory that allows non-linear distribution of axial stress through-the thickness would be required, as for example shown by Groh and Weaver (2015) for Fig. 9 Case (2) boundary conditions the case of prismatic beams, and will be considered in future work.
In the finite element analysis, each layer of the beam has a different mesh arrangement, which follows a structured pattern only for the second layer, as illustrated in Fig. 10. Consequently, the mesh density " k " varies through the thickness but is approximately uniform spanwise (where "k" indicates the layer index). More specifically, the ratio between mesh densities corresponds to 1 ∕ 2 = 1.5 and 1 ∕ 3 = 3 . The mesh refinement study is given in Table 4, with values of zz calculated at the upper surface of x = 0.25L . The adopted Poisson's ratio was (1) = (3) = 0.3 and (2) = 0.2 . It is noteworthy that the choice for triangular elements arose from the curved geometry of the second layer. The processing time of the finite element analysis was 25.5s against 0.041s from the proposed method.
The xx distribution is shown in Fig. 11. For both cross-sections, the current methodology matches Balduzzi et al.'s formulation (Balduzzi et al. 2017a) and the finite element analysis. However, a greater discrepancy is observed between the proposed theory and "FEA" results for x = 0.25L between f (2) and f (3) . More specifically, the locations of maximum stress do not agree but their magnitudes are identical. The simplification suggested by design codes (ACI Committee and International Organization for Standardization 2008; EN 2004) fails to capture the correct neutral axis position and stress discontinuities at material interfaces, resulting in significant underestimation of xx for the stiffer layer. Figure 12 depicts the shear stress distribution. Notably, the proposed theory matches well with the finite element analysis, but discrepancies were observed at the second layer for both cross-sections. This inconsistency is explained by the neglect of Poisson's ratio effects in the current method. State-ofthe-art formulations have not captured the shear stress distribution accurately due to the neglect of beam eccentricity ("Bal") and stiffer layer properties ("mod (E (2) = E (1) )"). Regarding the traction-free boundary requirement, all models successfully predict the vanishing of shear stress on the flat surface, but only   (2): the developed formulation matched the finite element analysis at the top surface. Finally, the transverse direct stress field is shown in Fig. 13. The developed formulation matches well with the finite element analysis solution. Overall, a greater discrepancy is observed for x = 0.60L . Furthermore, the simplification considered by assuming E (1) = E (2) does not capture the zz distribution well and overestimates the transverse direct stress at the top surface. Similar to the previous example, the maximum absolute zz magnitude occurs at the same position as the maximum shear stress, i.e., at the top surface. The cross-sections x = 0.25L and x = 0.75L have been selected to compare results due to their remote distance from localized effects occurring near boundaries in the finite element analysis. A (4.3) typical finite element mesh for case (3) is illustrated in Fig. 14; note that the number of elements in the thickness direction is constant spanwise. The convergence study is reported in Table 5, with zz values evaluated at the lower surface of cross-section x = 0.25L for the refinement analysis (Table 5 zz values refer to case (3.2)). The processing time of the finite element analyses were 19 s and 19.6 s against 0.067 s and 0.063 s from the proposed method for cases (3.1) and (3.2), respectively.

Case (3.1)
Consider the linearly asymmetric tapered beam with a longitudinal body force F x = 1 kN∕m 3 applied to all layers, as shown in Fig. 15. The body force resultant is a linear function of the x-coordinate given by T x = 1 − 0.05x (kN/m). The motivation to solve this load case is to validate the proposed formulation when N ′′ ≠ 0 , a variable necessary to define the transverse direct stress field. Furthermore, a longitudinal body force applied to a slender structure is a typical load case of beams under centrifugal forces, which occurs in wind turbine blades and helicopter rotor blades.
Approximate results for the internal forces of case (3.1) and their derivatives are reported in Table 6.
The longitudinal stress results suggest that all solutions match the same piecewise linear distribution depicted in Fig. 16. As expected, the gradient of xx decreases through the thickness as a result of Young's modulus reduction from the bottom to the top, noting that a similar trend is anticipated for the other stress components. Due to this particular material stiffness distribution, the centroid is located closer to the lower surface. Figure 17 shows that the proposed model matches the shear stress distribution of the finite element analysis. Notably, the curvature of xz is more pronounced for the cross-section at x = 0.75L rather than at x = 0.25L , which occurs because the contributions of the parabolic terms, i.e., N(x), M(x) and m(x), when added together, are more significant at x = 0.75L compared to x = 0.25L . In turn, the contribution of N � (x) , a linear function through the thickness, is accentuated at x = 0.25L . The solution given by Balduzzi et al. (2017a) does not capture the appropriate shear distribution because of the neglect of beam eccentricity, bending load,    and N � (x) . It is noteworthy that the maximum absolute shear stress occurs at the lower surface for x = 0.25L and at the interface between the second and the third laminae for x = 0.75L . Moreover, the current formulation successfully captures the traction-free boundary conditions at both surfaces. Figure 18 depicts the excellent agreement between the proposed methodology and the finite element analysis for the transverse direct stress. The largest discrepancy observed was 0.1% , located at the lower surface, as reported in Table 5. Notably, the maximum absolute zz occurs at the lower surface for x = 0.25L . In contrast, the maximum zz magnitude at x = 0.75L resides within the first layer, which is a different position from that of maximum shear stress. Additionally, the current formulation successfully captures the traction-free boundary conditions at the lower surface and the vanishing shear stress at the horizontal surface.

Case (3.2)
The loading case (3.2) illustrated in Fig. 19 is the linearly tapered beam of the former example subject to x = 1 kN∕m and an outwards pressure load t (5) z = 1 kN∕m , both applied to the horizontal surface, resulting in the internal forces given in Table 7.
The solutions for the longitudinal stress field shown in Fig. 20 result in the same piecewise linear distribution for all methodologies. Similar to the previous load case, the gradient of xx is proportional to the stiffness distribution, which is more significant for the lower layers. Note that the neutral axis position is closer to the centroid for cross-section x = 0.25L rather than at x = 0.75L because the bending moment is more pronounced in relation to the generalized longitudinal stress for x = 0.25L.
The shear stress prediction of the current methodology matches that of the finite element analysis, as depicted in Fig. 21. Notably, the analytical solution "Bal" has not captured the correct xz distribution due to three different reasons: the neglect of beam eccentricity and bending load, and the assumption of external load being generalized to the entire cross-section. The appropriate traction requirement has been successfully captured by the developed methods ( xz = 1 MPa at z = f (5) for both cross-sections). With regards to the xz profile, the   analyzed cross-sections exhibit different trends. The curvature of xz is more pronounced for x = 0.75L rather than x = 0.25L . Furthermore, the maximum absolute xz magnitude is located at the lower surface for x = 0.25L whereas at the interface between the first and second layers for x = 0.75L. The transverse direct stress component of the current model and the finite element analysis match the third-order piecewise function depicted in Fig. 22, including the correct magnitude of the pressure load at the horizontal surface. Furthermore, the position of the maximum zz magnitude is located at the interface between layers k = 1 and k = 2 for x = 0.25L whereas it is at the upper surface for x = 0.75L , which represent different positions of maximum shear. Hence, the maximum absolute shear and transverse direct stresses may not occur at the same coordinate.

Conclusion
An analytical methodology to predict the stress field of untwisted tapered beams comprising multiple isotropic layers has been derived. A state of plane stress with linear and small strains, known internal forces, and perfectly bonded layers, has been assumed. In the scope of laminated tapered beams, new developments include: (1) consideration of beam eccentricity and the first derivative of the longitudinal internal force to define the shear stress profile; (2) recovery of the transverse direct stress component; (3) layerwise body loads and layerwise traction forces; (4) closed-form solutions for the stress field of rectangular cross-sections.
Non-classical terms arise in the transverse stress field for heterogeneous non-prismatic beams due to the arbitrary stiffness distribution through the crosssection. More precisely, the longitudinal external load is typically applied eccentrically to the central axis because the neutral axis position depends on the stiffness distribution through the cross-section. As such, a bending load should be included in the equilibrium relations and considered in the derivation of stress field. Moreover, the term involving the first derivative of the generalized longitudinal stress does not vanish in the shear stress profile of laminated tapered beams, as it would be the case for homogeneous media under generalized external load. In conclusion, the position of maximum shear and transverse direct stresses can not be easily anticipated and may differ from each other. Also, the need to appropriately account for beam eccentricity is required to define the transverse stress field. The proposed theory is efficient from a computational standpoint as well as processing time being less than finite element analysis, not to mention the analytical formulation does not require mesh convergence studies and the design of each layer as in computer-aided design software. Furthermore, the results suggest excellent agreement with 2D-solid finite element analysis. Practitioners could benefit from the stress field derived in this work to design laminated non-prismatic beams under more complex loading conditions including partial crosssectional loads, discontinuous loads through the cross-section, pressure and traction forces. Potential further developments involve extending this work to composite laminates, anisotropic media, and functionally graded materials.