# Analytical interlaminar stresses of composite laminated beams with orthotropic tapered layers

M.M.S. Vilar, P. Khaneh Masjedi, D.A. Hadjiloizi, Paul M Weaver

► To cite this version:

M.M.S. Vilar, P. Khaneh Masjedi, D.A. Hadjiloizi, Paul M Weaver. Analytical interlaminar stresses of composite laminated beams with orthotropic tapered layers. Composite Structures, 2023, 319, pp.117063. 10.1016/j.compstruct.2023.117063. hal-05339405

HAL Id: hal-05339405

https://hal.univ-brest.fr/hal-05339405v1

Submitted on 30 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Analytical interlaminar stresses of composite laminated beams with orthotropic tapered layers

M.M.S. Vilar<sup>a,c,d,*</sup>, P. Khaneh Masjedi<sup>b</sup>, D.A. Hadjiloizi<sup>a</sup>, Paul M. Weaver<sup>a</sup>

## ARTICLE INFO

### Keywords:

- Non-prismatic beams

- Tapered beams

- Composite laminated beams

- Analytical stresses

- Interlaminar stresses

## ABSTRACT

Innovative manufacturing methodologies offer opportunity to fabricate tapered composite laminates without layer terminations. However, the complex mechanical behaviour of tapered composite beams has yet to be studied. We introduce an accurate cost-efficient analytical method in line with Timoshenko theory to predict transverse stresses of composite laminated beams comprising orthotropic tapered layers under transverse loading, including pressure loads and transverse body forces. The transverse stress components are derived by introducing taper into the lamina constitutive relation followed by Cauchy stress equilibrium. The main conclusions of this study include: (i) transverse stresses are fully coupled with the internal forces and their derivatives; (ii) transverse stresses can be discontinuous through the thickness; (iii) Classical Laminate Theory underestimates transverse stress magnitudes of composite laminated beams. Results are validated with 2D and 3D solid-like finite element analyses inferring high levels of accuracy for the transverse stresses valid up to 8° of taper.

## 1. Introduction

Driven by the design requirements of modern technologies, the demand for lightweight slender structures has been increasing in several engineering fields, including energy harvesting of wind turbine blades, aircraft wings and building projects. In this context, composite structures are widely employed as a primary choice for structural elements because they can achieve large stiffness-to-mass ratios, especially when combined with variable cross-sectional area [1–13] and elastic properties spanwise [14–19]. However, the conventional manufacturing method to produce variable cross-section composite laminates remains the ply drop-off technique since the late 20th century [20], which consists of cutting constant thickness plies resulting in detrimental stress concentrations at associated resin pockets [21–24]. In turn, a new technology can manufacture pre-impregnated carbon-fibre laminates with gradual thickness and width variation by coupling a pressure device to a laser-assisted tape placement technique [25,26], which mitigates stress concentration effects. Modelling such structure in conventional finite element commercial software can be both computationally demanding and time-intensive for model set-up. Nevertheless, simplified

and accurate methodologies are yet to be developed. Hence, we propose an analytical stress recovery method to model straight fibre composite laminated beams where layer thicknesses vary linearly spanwise.

The shape of variable cross-section beams can be characterised in groups with distinct functionalities, each possessing particular intrinsic properties to be considered in a beam theory. However, no terminology convention exists in the literature. To avoid ambiguity, the expression “non-prismatic beams” concerns beam elements with generic cross-sectional variations, and “tapered beams” denotes the specific case of linear spanwise variation. The object of study of this investigation regards multiple tapered layers beams with orthotropic constitutive relations referred to as “tapered composite beams”. The tapering effect introduces complex mechanical behaviour to beam elements not appropriately described in classical 1D beam theories. For example, the shear and bending stiffnesses are all coupled with the curvature and generalised shear strains [1]. Specific stress relations governed by Cauchy traction-stress equilibrium arise on beam surfaces [1,3,7,8,10, 12,13,27–35] and layer interfaces [28,32] due to both being oblique to the reference system. Regardless that these challenges are well-known in the literature, no refinement attempt has been made to available 1D

<footnote><sup>a</sup> Bernal Institute, School of Engineering, University of Limerick, Limerick, Ireland</footnote>

<footnote><sup>b</sup> Department of Aerospace and Mechanical Engineering, South East Technological University, Kilkenny Road Campus, Carlow, Ireland</footnote>

<footnote><sup>c</sup> Institut de recherche Dupuy De Lôme, Université Bretagne Sud, Lorient, France</footnote>

<footnote><sup>d</sup> Laboratoire Roberval, Université de Technologie de Compiègne, Compiègne, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: vilarmot@univ-ubs.fr, vilarmms@gmail.com (M.M.S. Vilar), pedram.masjedi@setu.ie (P. Khaneh Masjedi), Demetra.Hadjiloizi@ul.ie (D.A. Hadjiloizi), Paul.Weaver@ul.ie (P.M. Weaver).</footnote>
## Nomenclature

### Alphanumeric characters

$a_{ij}$ Membrane compliance coefficients (N$^{-1}$m)

$\hat{A}$ Segment of cross-sectional area (m$^2$)

$b_{ij}$ Compliance coupling between membrane strains and moment resultants (N$^{-1}$)

$c_{ij}$ Compliance coupling between curvatures and in-plane force resultants (N$^{-1}$)

$d_{ij}$ Bending compliance coefficients (Nm$^{-1}$)

$C_{ij}^{(k)}$ Orthotropic elastic constants in local coordinates (Pa)

$e$ Discrepancy between results (%)

$E_{il}$ Young's modulus (Pa)

$f_k$ z-coordinate of the kth layer lower interface (m)

$F_i^{(k)}$ Body force acting on the kth layer (Nm$^{-3}$)

$G_{ij}$ Shear modulus (Pa)

$h$; $h^{(k)}$ Beam height; layer thickness (m)

$H_f$ Coefficient accounting for traction forces on interfaces oblique to the coordinate system (-)

$L$ Beam length (m)

$m$ Spanwise distributed bending moment (Nm/m)

$M$ Cross-sectional bending moment resultant (Nm)

$n$ Number of layers (-)

$Oxyz$ Global Cartesian coordinate system (-)

$O\bar{x}\bar{y}\bar{z}$ Local Cartesian coordinate system (-)

$Q_{ij}^{(k)}$ Stress-strain lamina coefficients in global coordinates (Pa)

$t_z^{(k)}$ Pressure load acting on the lower interface of the kth layer (Pa)

$T_z$ External load cross-sectional resultant (Nm$^{-1}$)

$V$ Shear force cross-sectional resultant (N)

$w$ Beam width (m)

### Greek characters & symbols

$\alpha_{fk}$ Angle formed by the x-axis the kth layer lower interface (°)

$\alpha^{(k)}$ Angle formed by the x-axis and the mean value of the upper and lower kth layer interfaces (°)

$\alpha_B$ Beam taper angle (°)

$\Phi_i^{(k)}$ Stress constant originating from indefinite integration (Pa)

$\Gamma_j$ Terms of the equilibrium of internal forces and their derivatives

$\Delta$ Discontinuity operator applied to interfaces

$\epsilon_{ij}$ Strain tensor (-)

$\hat{\epsilon}_0$ Membrane strains (-)

$\epsilon_i^{(k)}$ Stiffness coefficients related to the transverse direct stress

$\eta_i^{(k)}$ Stiffness coefficients related to the transverse shear stress

$\theta^{(k)}$ Fibre orientation (°)

$\kappa$ Curvatures (m$^{-1}$)

$\lambda_i^{(k)}$ Stiffness coefficients related to the transverse direct stress

$\nu$ Poisson's ratio (-)

$\rho_i^{(k)}$ Stiffness coefficients related to the transverse shear stress

$\sigma_{ij}$ Stress tensor (Pa)

$\psi$ Shear stress stiffness coefficients (m$^{-1}$)

beam formulations to include the tapering effects on composite tapered beam models, which is the primary objective of this study.

Analytical stress recovery methods of homogeneous non-prismatic beams with rectangular cross-sections have been under debate since the

early 20th century. For instance, Bleich [31] adapted the Jourawski formula [36], i.e., the parabolic shear stress profile of prismatic beams to symmetric tapered shapes. Concurrently, several studies using elasticity theory derived the shear stress of infinite wedge elements [37–39], followed by truncated wedge elements of symmetric [40,41] as well as I-shape cross-sections [30,42,43]. Boley [35] tested the limits of Navier's hypothesis in non-prismatic beams under pure bending via Airy's stress potential and concluded that the error associated with such simplification escalates with the rate that the cross-section varies. Romano [44] predicted the shear stress of tapered beams under Timoshenko's kinematic assumptions, and Zhou et al. [33] of asymmetric tapered beams. Variational asymptotic methods have been developed to model symmetric tapered beams subject to extension, bending, and shear loads [12,13]. Baluzzi et al. [1] reconstructed the shear stress of non-prismatic beams by superposing Jourawski's formulation to a linear function that satisfies traction-free boundary equilibrium. Beltempo et al. [10] proposed a variational formulation based on the Hellinger-Reissner Principle to predict the mechanical behaviour of non-prismatic beams, but despite being a mixed-energy method, the transverse direct stress is linearly interpolated through the beam's depth. Later, Mercuri et al. [9] adapted the cubic transverse stress profile obtained from Cauchy-stress equilibrium considerations to non-prismatic beams under concentrated loading. Subsequently, Vilar et al. [27,29] extended the former theory to account for beams subject to general spanwise distributed loads, uniformly and partially distributed through the cross-section, respectively. Migliaccio et al. [45] relaxed the Navier's longitudinal stress assumption to predict the shear stress profile considering deformable cross-sections with large displacements. Analytical solutions for thin-walled cross-sections have also gained increasing attention in recent years. Relevant research investigated the stress field of beams with box girder cross-sections symmetrically tapered in one [6,7] and two directions [8], as well as being asymmetrically tapered [34]. Additionally, Chockalingam et al. recovered the shear stress of I-shaped cross-section tapered beams [46,47].

Research considering stress fields of beams comprising tapered layers has been increasing due to the interest on heterogeneous lightweight-stiff structures. Nevertheless, state-of-the-art analytical methods compromise applicability for the sake of modelling simplicity. For example, Gayen and Roy [48] investigated the stress distribution of circular constant-thickness tapered composite beams under hydrothermal loading. Vel et al. [49] proposed an interlaminar stress recovery procedure to sandwich tapered panels with an isotropic soft core and anisotropic constant-thickness skins. Baluzzi et al. [32] and Vilar et al. [28] derived the interlaminar stresses of multilayered isotropic non-prismatic beams. Ai and Weaver [50] utilised the Ritz method and a layerwise displacement field to predict the mechanical response of symmetrically tapered sandwich beams with a functionally graded core. Within the scope of concrete-steel box girder cross-sections, Zhou et al. [34,51] derived the shear stress distribution of asymmetric tapered and non-prismatic beams, respectively. Later, Zhou et al. [52] refined the former investigation to include the Resal effect, i.e., an effective shear force reduction due to vertical compressive force component acting on tapered flanges. Marín and Graciani [53] developed a strength of material theory to determine normal stress flow acting on composite wingbox cross-sections subject to a torsional load. As noted,
<figure><img src="image_3.png" /><figcaption>Fig. 1. (a) Side view of a tapered composite beam; (b) arbitrary cross-section.</figcaption></figure>

all previous studies restrict geometric features or material character-
istics for simplification. To the authors' knowledge, no analytical 1D
beam theory has been appropriately refined to model the complex stress
field of composite beams with orthotropic tapered layers.

To partly address the state-of-the-art limitations, we introduce a
cost-efficient yet accurate analytical methodology in line with the
Timoshenko beam theory to predict interlaminar stresses of compos-
ite beams comprising orthotropic tapered layers with a rectangular
cross-section and small thickness variation. The beam deforms lin-
ear elastically under plane stress in conjunction with small strains
and displacements. The adopted methodology invokes Cauchy stress
equilibrium to recover the transverse stress field from Classical Lam-
inate Theory (CLT) in-plane stress expressions considering the taper
effect, followed by enforcing appropriate interfacial boundary con-
ditions. Body forces and pressure loads in the transverse direction
are prescribed in a layerwise fashion, allowing consideration of non-
uniform loads through the cross-section in a 1D beam theory. To
validate the proposed method, numerical examples of tapered com-
posite laminated beams with symmetric and non-symmetric stacking
sequences are considered for taper angles up to 8°. Results show good
agreement with 2D and 3D solid-like finite element analyses, while
the CLT does not predict the transverse direct stresses with reasonable
accuracy.

The outline of this work is as follows: Section 2 introduces the hypotheses and limitations as well as the equilibrium and constitutive relations. Section 3 recovers the transverse stress components. Section 4 compares results of the developed formulation with finite element analyses. Finally, Section 5 summarises the novelties of this study and suggests future work.

## 2. Problem idealisation

Relevant remarks and definitions regarding the tapered beam's
mechanical response, as well as for the equilibrium of force resultants,
constitutive relations, and interfacial conditions, are introduced in this
section.

Consider the tapered beam illustrated in Fig. 1 comprising *n* homo-
geneous orthotropic layers of length *L* and width *w* measured along the
*x*- and *y*-axes, respectively. A total of (*n* + 1) known interface functions
*f*<sub>*k*</sub>(*x*) are introduced corresponding to the *z*-coordinate each *k*th layer's
lower interface, noting that *f*<sub>1</sub> refers to the bottom surface and *f*<sub>*n*+1</sub> to
the top surface. The *beam height* is calculated by *h*(*x*) = *f*<sub>*n*+1</sub> − *f*<sub>1</sub>.

The taper angle is a metric commonly used to quantify the spanwise
variation of a non-prismatic beam's cross-sectional area. In multilay-
ered media, the rate of change in the thickness of a specific layer affects
the mechanical response of the overall laminate. To clarify the adopted
terminology, the *interface taper angle* α<sub>f_k</sub> is delimited by f_k and the *x*-axis, evaluated by α<sub>f_k</sub> = arctan f'<sub>k</sub>. Secondly, the *lamina taper angle* α<sup>(k)</sup>

refers to the mean value between two adjacent interface taper angles.
Finally, the *beam taper angle* is calculated by α<sub>B</sub> = α<sub>f<sub>n+1</sub></sub> − α<sub>f<sub>1</sub></sub> (positive
counter-clockwise).

It is convenient to adopt a local orthonormal coordinate system Ox̄ȳz̄ to describe the fundamentals of an arbitrary lamina. The x̄ȳ-plane corresponds to a given lamina's midplane, where the x̄- and ȳ-axes are oriented to the material principal directions. The angle between the x̄-axis projection to the xz-plane and the x-axis matches α<sup>(k)</sup> (Fig. 2(a)), and the fibre orientation θ<sup>(k)</sup> (Fig. 2(b)) is defined by the angle between the x-axis projection to the x̄ȳ-plane and the x̄-axis.

It is assumed that the lamina thickness is relatively thin under a
state of plane stress that undergoes linear strains together with small
displacements. After deformation, the cross-section remains plane but
not necessarily perpendicular to the neutral axis, hence is compatible
with Timoshenko beam theory.

The adopted hypotheses are only valid for cases where the lamina taper angle is “sufficiently small”, and whose reference magnitude has not been validated for tapered composite beams. With regards to isotropic materials, previous related research reported that assuming Navier’s hypothesis to a 10°-taper results in a 7.5% error in the estimation of the longitudinal stress [35]. In many practical applications, taper angles in beam-like structures are relatively small. For example, a 10 MW wind turbine blade appears to have a rapidly varying cross-section near the root, but it has a maximum beam taper angle of 5.96° [6].

Further restrictions are imposed to ensure the problem is well
formulated:

(1) the cross-sectional dimensions are small compared to the beam
length to meet the slenderness requirement of beam theory;

(2) Layer interfaces are perfectly bonded and their projections to the cross-section are parallel to the xy-plane;

(3) cut-outs are disregarded to avoid stress concentration effects;

(4) force and bending moment resultants act in the xz-plane, not-
ing that only the couplings between curvature and in-plane
force resultants, and membrane-bending moment resultants, are
considered.

Although no limitation on stacking sequence is directly imposed,
remark (4) suggests the accuracy of the proposed method reduces
with an increase in the laminate anisotropy due to non-symmetry.
Therefore, the proposed formulation is recommended for beam bending
in the thickness direction, with limited cross-sectional warping and
bending in the width direction. This restriction is further investigated in
Section 4.2, where numerical examples with laminates of symmetrical
and non-symmetrical layups are solved and their accuracy compared.
Furthermore, the cross-section is defined with respect to the global
Oxyz coordinate system, hence it is non-orthogonal to the beam's neu-
tral axis, noting that related research on non-prismatic beams adopted
the same cross-sectional definition [1–3,9,27–29,32,54].
<figure><img src="image_4.png" /><figcaption>Fig. 2. Perspectives of tapered lamina: (a) side view; (b) top view.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Equilibrium of force resultants.</figcaption></figure>

## 2.1. Force resultant equilibrium

Consider the infinitesimal tapered beam segment depicted in Fig. 3 subject to body loads $F_i(x, z)$ and surface forces $t_i(x, z)$. The shear force cross-sectional resultant $V(x)$ and the resultant bending moment $M(x)$ are known variables, not necessarily evaluated from statically determinate cases, given by the following classical expressions

$$V(x) = w \int_{f_1}^{f_{n+1}} \tau_{xz} dz \qquad (2.1a)$$

$$M(x) = w \int_{f_1}^{f_{n+1}} z \sigma_{xx} dz \qquad (2.1b)$$

where $\sigma_{ij}$ is the stress tensor and $\tau_{xz}$ is the shear stress component acting on the $xz$-plane.

The external load is prescribed in a layerwise fashion such that the body forces $F_z$ and the pressure loads $t_z$ are uniformly distributed within the layer boundaries as expressed by

$$\begin{aligned} F_z(x, z) &= F_z^{(k)} | \{f_k \le z < f_{k+1}\} \\ t_z(x, z) &= t_z^{(k)} | \{z = f_k\} \end{aligned} \qquad (2.2)$$

At a given cross-section, the resultant of the external load $T_z(x)$ is given by

$$T_z(x) = w \int_{f_1}^{f_{n+1}} (F_z + t_z) dz \qquad (2.3)$$

Relations between stress resultants and external loads are derived for tapered beams by considering the equilibrium of forces and moments. In particular, the internal force resultants, or generalised stresses, vary linearly through an infinitesimal segment. It is assumed

that the Timoshenko equilibrium expressions are valid as given in the following

$$V' = -T_z \qquad (2.4a)$$

$$M' + V = -m \qquad (2.4b)$$

It is noteworthy that the equilibrium of force resultants introduced in this section is valid for general non-prismatic beam shapes [1,2,9,27–29,32,55]. For isotropic beams, Eq. (2.4) simplifies to the equilibrium equations of homogeneous non-prismatic beams [1,2,27,55].

From now on, the resultant bending moment $M$ is referred to as the “bending moment” while “bending load” denotes the spanwise distributed bending moment $m$.

## 2.2. Constitutive relations

Relevant considerations for the lamina and laminate stress/strain relations are discussed in this section. In particular, the plane stress condition imposed on a traction-free arbitrary lamina is expressed as follows

$$\left\{ \begin{array}{c} \sigma_{\bar{x}\bar{x}} \\ \sigma_{\bar{y}\bar{y}} \\ 0 \\ 0 \\ \tau_{\bar{x}\bar{y}} \end{array} \right\}^{(k)} = \begin{bmatrix} C_{11} & C_{12} & C_{13} & 0 & 0 & 0 \\ C_{12} & C_{22} & C_{23} & 0 & 0 & 0 \\ C_{13} & C_{23} & C_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{66} \end{bmatrix}^{(k)} \begin{bmatrix} \epsilon_{\bar{x}\bar{x}} \\ \epsilon_{\bar{y}\bar{y}} \\ \epsilon_{\bar{x}\bar{y}} \\ \gamma_{\bar{y}\bar{z}} \\ \gamma_{\bar{x}\bar{z}} \\ \gamma_{\bar{x}\bar{y}} \end{bmatrix} \qquad (2.5)$$

where $\epsilon_{ij}$ is the strain tensor, $C_{ij}^{(k)}$ corresponds to the orthotropic material elastic constants that determines the state of stress for a given state of strain in the $O\bar{x}\bar{y}\bar{z}$ system, whose expressions are derived from the Young’s modulus $E_{il}$, the shear modulus $G_{ij}$ and the Poisson’s ratio $\nu_{ij}$ [56–58].

Because of the $\bar{x}\bar{y}$-plane stress hypothesis, the transverse strains ($\epsilon_{zz}$, $\gamma_{xz}$ and $\gamma_{yz}$) can be expressed as a linear combination of the in-plane strains. Therefore, the reduced stress–strain lamina relation in global coordinates yields

$$\{\hat{\sigma}\} = [Q(\alpha^{(k)}, \theta^{(k)})] \{\hat{\epsilon}\} \qquad (2.6)$$

where $\{\hat{\sigma}\} = \{(\cdot)_{xx}, (\cdot)_{yy}, (\cdot)_{xy}\}^T$, and $Q_{ij}^{(k)}$ is the reduced stiffness matrix transformed from local to global coordinates for the $k$th layer.

The relation between the force resultants and generalised strains is described by the well-known ABD constitutive matrices [56–58] and the following compliance expressions

$$\begin{bmatrix} \hat{\epsilon}_0 \\ \hat{k} \end{bmatrix} = \begin{bmatrix} a(x) & b(x) \\ c(x) & d(x) \end{bmatrix} \begin{bmatrix} 0 \\ \hat{M} \end{bmatrix} \qquad (2.7)$$

where $\hat{\epsilon}_0$ and $\hat{k}$ are the membrane strains and curvatures, respectively, such that $\hat{\epsilon} = \hat{\epsilon}_0 + z\hat{k}$. The compliance coefficients are evaluated
<figure><img src="image_6.png" /><figcaption>Fig. 4. Interface equilibrium in: (a) x-direction; (b) y-direction; (c) z-direction.</figcaption></figure>

from the inversion of the ABD matrices: $a_{ij}(x)$ represents membrane compliance; $b_{ij}(x)$ membrane strain-moment coupling; $c_{ij}(x)$ curvature-in-plane force coupling and $d_{ij}(x)$ bending compliance.

In general, the $\bar{x}\bar{y}$-plane is not under a state of plane stress because it is oblique to both lamina surfaces, given the arbitrariness of $f_k$ and $f_{k+1}$. Unless the lamina surfaces are parallel, the stress components transverse to the $\bar{x}\bar{y}$-plane do not vanish on both interface surfaces simultaneously. As such, assuming the $\bar{x}\bar{y}$-plane is under plane stress undermines the physics of the proposed problem, resulting in an error that escalates with $a^{(k)}$ (in line with the small taper angle assumption). Further discussion on the repercussions of this limitation can be found in Section 4.2, where numerical examples are studied for multiple values of $a^{(k)}$.

## 2.3. Interfacial boundary equilibrium

The equilibrium of force resultants invoked on an infinitesimal interfacial surface point results in specific stress relations between adjacent layers. For traction-free prismatic laminates, a consistent stress field requires that transverse stress components are continuous through the thickness [59]. Conversely, recent investigations on multilayered tapered beams have reported that transverse stresses expressed in global coordinates can be discontinuous at interfaces [28,32].

It is possible to derive the interlayer boundary condition by invoking the equilibrium of force resultants on an infinitesimal representative unit cell located at layer interfaces, as illustrated in Fig. 4.

In the $xy$-plane, the infinitesimal unit cell has independent dimensions $dx$ and $dy$, with dimension $f'_k$dx in the $z$-direction. Therefore,

considering the equilibrium of force resultants in the proposed unit cell results in the following expressions

$$ \Delta\tau_{xz}^{(k)} = f'_{k}\Delta\sigma_{xx}^{(k)} \qquad (2.8a) $$

$$ \Delta\tau_{yz}^{(k)} = f'_{k}\Delta\tau_{xy}^{(k)} \qquad (2.8b) $$

$$ \Delta\sigma_{zz}^{(k)} = f'_{k}\Delta\tau_{xz}^{(k)} - t_{z}^{(k)}H_{k} \qquad (2.8c) $$

where $H_k$ is the non-dimensional length of the oblique line shown in Fig. 4, and the $\Delta(.)$ operator indicates a discontinuity of its argument evaluated at $f_k$

$$ \Delta(.) = (.)^{(k)} - (.)^{(k-1)} \qquad (2.9a) $$

$$ H_k = \sqrt{1 + (f'_k)^2} \qquad (2.9b) $$

noting that the second term in the right-hand side of Eq. (2.9a) vanishes for the first layer.

The interface conditions expressed in Eq. (2.8) are valid for points remote from the lateral surfaces in the $xz$-plane located at the interfaces of tapered non-prismatic beams. Thus, Eq. (2.8) does not satisfy the equilibrium conditions on the $xz$-lateral surfaces since the respective state of stress in this region for an actual structure is more complex than the plane stress assumption adopted here. Additionally, no small taper angle approximation has been assumed to derive Eq. (2.8); hence the introduced interface conditions are valid regardless of the magnitude of $\alpha_{f_k}$.

It should be noted that the local stresses at interfaces are continuous due to equilibrium conditions, and the local stress transformation to the global coordinate system results in equivalent relations as those
<figure><img src="image_7.png" /><figcaption>Fig. 5. Local and global stresses in equilibrium on layer interfaces.</figcaption></figure>

expressed in Eq. (2.8). To confirm this information, consider the equilibrium of force resultants between local and global stresses of the infinitesimal element at the traction-free layer interface illustrated in Fig. 5.

Invoking the equilibrium of forces in the x- and z-directions of the lower infinitesimal elements of Fig. 5 results in the following expressions

$$-\sigma_{xx}^{(k)} \tan \alpha_k + \tau_{xz}^{(k)} + (\sigma_{\bar{z}\bar{z}}^{(k)} + \tau_{\bar{x}\bar{z}}^{(k)}) \sec \alpha = 0 \quad (2.10)$$

$$\sigma_{xx}^{(k+1)} \tan \alpha_k - \tau_{xz}^{(k+1)} - (\sigma_{\bar{z}\bar{z}}^{(k+1)} + \tau_{\bar{x}\bar{z}}^{(k+1)}) \sec \alpha = 0 \quad (2.11)$$

Analogously, the equilibrium of forces in the x- and z-directions of the upper infinitesimal elements of Fig. 5 yields

$$\tau_{xz}^{(k)} \tan \alpha_k - \sigma_{\bar{z}\bar{z}}^{(k)} + (\sigma_{\bar{z}\bar{z}}^{(k)} - \tau_{\bar{x}\bar{z}}^{(k)}) \sec \alpha = 0 \quad (2.12)$$

$$-\tau_{xz}^{(k+1)} \tan \alpha_k + \sigma_{\bar{z}\bar{z}}^{(k+1)} - (\sigma_{\bar{z}\bar{z}}^{(k+1)} - \tau_{\bar{x}\bar{z}}^{(k+1)}) \sec \alpha = 0 \quad (2.13)$$

According to Newton's Third Law, a force acting on the hypotenuse of the lower infinitesimal element is equal in magnitude and opposite in direction to that acting on the upper element, leading to the continuity of local stresses $\tau_{\bar{x}\bar{z}}^{(k)} = \tau_{\bar{x}\bar{z}}^{(k+1)}$ and $\sigma_{\bar{z}\bar{z}}^{(k)} = \sigma_{\bar{z}\bar{z}}^{(k+1)}$. Substituting Newton's Third Law to Eqs. (2.10) to (2.13) results in direct relations between global stresses acting on a given interface obtained in Eq. (2.8), which analytically proves that the local stresses are continuous through layer interfaces. This subject is further discussed in the numerical example given in Section 4.2.

As such, a discontinuity in transverse stresses at layer interfaces exists in global coordinates if two conditions are met: (i) the interface surface is oblique to the reference system ($f'_k \neq 0$); (ii) the elastic constants or the material orientation changes from one layer to the next. Furthermore, body forces do not contribute to the interface force equilibrium due to being negligible compared to the stress and traction resultants, as body loads are applied to infinitesimal volumes while the latter are applied to infinitesimal areas.

### 3. Transverse stresses recovery

This section introduces analytical expressions for the in-plane and transverse stress fields of tapered composite beams. The proposed formulation recovers the transverse stress field from Cauchy stress equilibrium expressions, starting with the shear components, followed by the imposition of interface stress conditions discussed in Section 2.3.

#### 3.1. In-plane stresses

Substituting Eq. (2.7) into Eq. (2.6) yields the well-known CLT in-plane stress expressions

$$\{\hat{\sigma}\} = 1/w [Q]^{(k)} ([b] + [d]z\{\hat{M}\}) \quad (3.1)$$

To simplify further operations, the in-plane stresses are expressed as a linear combination by

$$\hat{\sigma}_i = 1/w \sum_{j=1}^{3} Q_{ij}^{(k)} (b_{j1} + d_{j1}z) M \quad (3.2)$$

noting that because the proposed method does not consider moments and forces in the width direction, only the first column of $a_{ij}$, $b_{ij}$, $c_{ij}$ and $d_{ij}$ are necessary to define the in-plane stresses, thus the notation of these variables are simplified to $a_j$, $b_j$, $c_j$ and $d_j$ from now on [56].

#### 3.2. Transverse shear stresses

The stress recovery methodology is similar for both transverse shear stresses; thus, $\tau_{xz}$ is obtained, and $\tau_{yz}$ is analogously derived subsequently.

Cauchy stress equilibrium in the longitudinal direction applied to an interior point of an arbitrary layer yields

$$\tau_{xz}^{(k)}(x, z) = - \iint \frac{\partial \sigma_{xx}^{(k)}}{\partial x} dA \quad (3.3)$$

where $dA = w dz$

To simplify operations on Eq. (3.3), the vector $\Gamma$ containing the terms of the equilibrium relations and their derivatives is introduced

$$\Gamma = \{V m M V' m'\}^T \quad (3.4)$$

The longitudinal Cauchy stress integrated in the thickness direction can be stated as a linear combination as follows

$$\tau_{xz}^{(k)}(x, z) = \sum_{i=1}^{3} (\eta_{i}^{(k)}(x, z)\Gamma_{i}(x)) + C_{1}^{(k)}(x) \quad (3.5)$$

where $C_{1}^{(k)}$ is the constant originating from indefinite integration.

Variables $\eta_{i}^{(k)}$ are stiffness coefficients that relate force resultants to $\tau_{xz}$ (e.g. $\eta_{M}^{(k)}$ relates bending moment to shear stress). As such, $\eta_{i}^{(k)}$
can be understood as layerwise functions that recover the shear stress profile terms from equilibrium relations and their derivatives. The shear stress profile can be deduced for different laminated materials from Eq. (3.5) by defining appropriate compliance coefficients discussed in Section 2.2. In a previous investigation, $\eta_i^{(k)}$ coefficients were obtained for multilayered beams [28] under longitudinal and transverse loading. For orthotropic tapered laminated beams under transverse loading, expressions for $\eta_i^{(k)}$ are given by

$$ \eta_V^{(k)} = \eta_m^{(k)} = Q_{1j}^{(k)} \left[ \frac{1}{2} d_j z^2 + b_j z \right] \quad (3.6a) $$

$$ \eta_M^{(k)} = -Q_{1j}^{(k)} \left[ \frac{1}{2} d_j' z^2 + b_j' z \right] \quad (3.6b) $$

where the dummy index indicating the first column of compliance coefficients is dropped to avoid repetition (e.g. $a_{j1}$ is simplified to $a_j$), and summation from $j=1$ to $j=3$ is implicit.

The transverse shear stress dependency on all internal forces has been referred to in the literature as the non-triviality of stress field [30]. Additionally, $\tau_{xz}$ can be understood as a superposition of the classical parabolic shear stress of prismatic beam theory ($\eta_V^{(k)}$ and $\eta_M^{(k)}$) and the contribution of the bending moment originating from variable cross-section ($\eta_M$).

The interfacial equilibrium is imposed by evaluating $C_1^{(k)}$ after substituting $z = f_k$ and Eq. (2.8a) in Eq. (3.5).

$$ C_1^{(k)} = - \sum_{i=1}^{3} \left( \eta_i^{(k)}(x, f_k) \Gamma_i(x) + f_k' \Delta \sigma_{xz}^{(k)} + \tau_{xz}^{(k-1)}(x, f_k) \right) \quad (3.7) $$

For $k=1$, in cases where no traction is prescribed, $C_1^{(k)}$ reduces to the traction-free boundary condition for homogeneous beams well-established in previous works [1,13,27,29–32,60]. Furthermore, $C_1^{(k)}$ requires shear stress evaluation of the lower adjacent layer at $f_k$. For example, the $\tau_{xz}^{(2)}$ profile can be determined only if $C_1^{(1)}$ and $\tau_{xz}^{(1)}$ have already been evaluated. As such, after performing successive substitutions from the $k$th to the bottom layer, the shear stress throughout-thickness distribution is expressed by

$$ \tau_{xz}^{(k)} = \sum_{i=1}^{k} \left( \sum_{j=1}^{3} \left( \left[ \eta_i^{(k)} - \Delta \eta_i^{(l)} \right] \Gamma_i + \Delta \tau_{xz}^{(l)} \right) \right) \quad (3.8) $$

dropping the $\Delta$ operator dependency on $f_k$.

The methodology to derive $\tau_{yz}$ is similar to that of $\tau_{xz}$. Cauchy stress equilibrium in the y-direction is integrated through a segment of cross-sectional area

$$ \tau_{yz}^{(k)}(x, z) = - \iint \frac{\partial \tau_{xy}^{(k)}}{\partial x} dA \quad (3.9) $$

Subsequently, the interface boundary condition given by Eq. (2.8b) are imposed analogously to Eq. (3.7), resulting in

$$ \tau_{yz}^{(k)} = \sum_{i=1}^{k} \left( \sum_{j=1}^{3} \left( \left[ \rho_i^{(k)} - \Delta \rho_i^{(l)} \right] \Gamma_i + \Delta \tau_{xy}^{(l)} \right) \right) \quad (3.10) $$

where $\rho_i$ coefficients hold the same expressions as $\eta_i$ but $Q_{1j}^{(k)}$ is replaced by $Q_{6j}^{(k)}$ in the latter.

### 3.3. Transverse direct stress

Consider Cauchy stress equilibrium in the transverse direction

$$ \sigma_{zz}^{(k)} = - \iint \left( \frac{\partial \tau_{xz}^{(k)}}{\partial x} + F_z^{(k)} \right) dA \quad (3.11) $$

Expressing the transverse direct stress as a linear combination and substituting $\tau_{xz}$ given by Eq. (3.8) into Eq. (3.11) yield

$$ \sigma_{zz}^{(k)} = \sum_{i=1}^{5} \left( \xi_i^{(k)} \Gamma_i \right) + \left( \frac{F_z^{(k)}}{T_z} - C_1^{(k)} \right) z + C_2^{(k)} \quad (3.12) $$

where $C_2^{(k)}$ is the constant from the indefinite integration of Eq. (3.11). Note that the definitions of external load resultants and equilibrium of force resultants, i.e., Eqs. (2.3) and (2.4) were invoked to derive Eq. (3.12). Furthermore, the summation operator varies from $i=1$ to $i=5$ to account for all terms of $\Gamma$ given by Eq. (3.4).

Likewise Eq. (3.5), the transverse direct stress given by Eq. (3.12) is undefined because it does not assume specific materials properties. A previous study deduced expressions for $\xi_i^{(k)}$ for the case of multilayered non-prismatic beams [28]. For composite beams comprising orthotropic tapered layers, $\xi_i^{(k)}$ are expressed as

$$ \begin{aligned} \xi_i^{(k)}(x, z) &= - \iint \lambda_i^{(k)}(x, z) dA \\ \lambda_V^{(k)} &= \lambda_m^{(k)} = \frac{\partial \eta_V^{(k)}}{\partial x} - \eta_M^{(k)} \quad (3.13) \\ \lambda_M^{(k)} &= \frac{\partial \eta_M^{(k)}}{\partial x} \\ \lambda_{V'}^{(k)} &= \lambda_{m'}^{(k)} = \eta_V^{(k)} \end{aligned} $$

noting that closed form solutions for $\xi_i^{(k)}$ and $\lambda_i^{(k)}$ can be found in Appendix.

Analogous to $C_1^{(k)}$, the constant $C_2^{(k)}$ is deduced by letting $z = f_k$ and invoking the interlaminar boundary condition in the z-direction expressed in Eq. (2.8c)

$$ C_2^{(k)} = - \sum_{i=1}^{5} \left( \xi_i^{(k)} \Gamma_i + C_1^{(k)} \bar{A} + f_k' \Delta \tau_{xz}^{(k)} + \sigma_{zz}^{(k-1)} \right) \quad (3.14) $$

Following an analogous method of successive substitutions explained in Eq. (3.8), an analytical expression for the transverse direct stress is given by

$$ \sigma_{zz}^{(k)} = \sum_{i=1}^{k} \left( \sum_{j=1}^{5} \left( \left[ \xi_j^{(k)} - \Delta \xi_j^{(l)} + \Delta \lambda_j^{(l)} \hat{A}^{(l)} \right] \Gamma_i + \left( \frac{F_z^{(k)}}{T_z} - \Delta \tau_{xz}^{(l)} \right) \hat{A}^{(l)} + \Delta \sigma_{zz}^{(l)} \right) \right) \quad (3.15) $$

The introduced variable $\hat{A}^{(l)}$ represents a segment of cross-sectional area described as follows

$$ \hat{A}^{(l)} = (z - f^{(l)}) | z > f_l \quad (3.16) $$

### 4. Numerical results

Numerical examples are presented in this section to give insight into the merits of the current methodology. Results are compared to CLT as well as to 2D and 3D solid-like finite element analysis solutions, referred to as “CLT”, “FEA 2D” and “FEA 3D”, respectively, while the current formulation is denoted as “meth”. The following beam cases are proposed:

(1) a symmetric sandwich beam with different beam taper angles ($\alpha_1 = -4^\circ$ and $\alpha_2 = -8^\circ$) subject to a concentrated transverse force at the tip;

(2) an asymmetric 10-layer tapered beam case with different stacking sequences, i.e., unidirectional [45]$_{10}$, crossply [0/90]$_5$, and non-symmetric [0$_5$/90$_5$], subject to a pressure load linearly distributed over the lower surface.

The terminology “Case (.)” refers to the beam cases (1) and (2) solved in this section.

A mesh refinement study was performed to ensure finite element analyses had converged, where the number of elements between mesh refinements was multiplied by a factor of 1.5 to 2.5. Relative errors defined as $e_i = (\sigma_{zz}^{(l)} - \sigma_{zz}^{(l-1)})/\sigma_{zz}^{(l)}$ were calculated between mesh refinements at specific locations (index $i$ refers to the $i$th mesh refinement), and the stopping criteria was $e_i \le 0.1\%$. Finally, the relative error of the proposed methodology in relation to the converged finite element
<figure><img src="image_8.png" /><figcaption>Fig. 6. Case (1) boundary conditions.</figcaption></figure>

<table><caption>**Table 1**
Orthotropic elastic properties.</caption><thead><tr><th>E<sub>1</sub></th><th>E<sub>2</sub> = E<sub>3</sub></th><th>G<sub>23</sub></th><th>G<sub>12</sub> = G<sub>13</sub></th><th>ν<sub>23</sub></th><th>ν<sub>12</sub> = ν<sub>13</sub></th></tr><tr><th>GPa</th><th></th><th>GPa</th><th></th><th>(-)</th><th></th></tr></thead><tbody><tr><td>150</td><td>12</td><td>3.2</td><td>4.4</td><td>0.45</td><td>0.25</td></tr></tbody></table>

mesh is given by $e_{meth} = (\sigma_{zz}^{(meth)} - \sigma_{zz}^{(i)})/\sigma_{zz}^{(meth)}$. It is relevant to highlight that $e_{meth}$ should not be interpreted as a measure of the accuracy of the proposed method because it is evaluated at a specific location in the beam.

To extract results from the numerical solutions, a straight path is created from the beam's lower surface ($f_1$) to its upper surface ($f_{n+1}$) at selected cross-sections, noting that the selected path intersects the width midpoints for the case of the 3D finite element analyses. More information on the mesh structure and element type is given in the respective beam case discussion.

## 4.1. Symmetrically tapered sandwich

Consider the symmetric sandwich cantilever beam illustrated in Fig. 6, of length $L = 100$ mm and width $w = 0.2$ mm.

The core, made of isotropic material ($E = 150$ MPa and $ν = 0.25$), separates the stiffer orthotropic skins of stacking sequence [0/90/0] made of 0.2 mm-thickness layers with the elastic properties given in Table 1. Two magnitudes of beam taper are considered: $α_1 = -4°$ and $α_2 = -8°$, resulting in the beam height varying from $h^{(4)}(0) = 15.2$ mm at the fixed end to $h^{(4)}(L) = 8.2$ mm, for the case of $α_1$, and to $h^{(4)}(L) = 1.2$ mm for $α_2$.

The following interface surfaces for Case(1) are given by (in mm)

$$ \left\{ \begin{array}{l} f_k = \tan(\alpha_i)x + 0.2(k-1) - 7.6 \\ f_{(2+n-k)} = -f_k |k = (1, 2, 3, 4), n = 7 \end{array} \right. \qquad (4.1) $$

The shear force acting at $x = 0.25L$ is $V = -1$ kN and the bending moment is $M = 75$ Nm.

A side view of a typical finite element mesh for Case (1) is depicted in Fig. 7(a) and (b) for models $α_1$ and $α_2$, respectively. Regarding the 2D model, 8-node quadrilateral plane stress elements (CPS8R in Abaqus terminology) were used for the skins and the core of model $α_1$, while the 6-node modified quadratic plane stress triangle (CPS6M) was used for the core of $α_2$ due to its triangular shape. Concerning the 3D finite element analysis, the 20-node quadratic brick element type (C3D20R) was assigned to the elements in the skin region and to the core of the $α_1$ model while the 15-node quadratic triangular prism (C3D15) was used for the core of the $α_2$ numerical example.

The finite element studies reported in Tables 2 and 3 were performed to assure the convergence of the numerical method, where values for $σ_{zz}$ have been evaluated at the lower surface (converged mesh information is highlighted blue in the respective table), noting that the proposed methodology resulted in $σ_{zz} = -89.482$ kPa ($α_1 = -4°$)

<table><caption>**Table 2**
Case (1) mesh convergence analysis ($α_1 = -4°$).</caption><thead><tr><th>FEA model</th><th>Number of elements</th><th>σ<sub>zz</sub> (kPa)</th><th>e<sub>i</sub> (%)</th><th>e<sub>meth</sub> (%)</th></tr></thead><tbody><tr><td rowspan="2">2D</td><td>871</td><td>-90.71</td><td>-</td><td>-1.38</td></tr><tr><td>1,330</td><td>-90.695</td><td>-0.016</td><td>-1.37</td></tr><tr><td rowspan="2">3D</td><td>1,500</td><td>-90.676</td><td>-</td><td>-1.34</td></tr><tr><td>3,325</td><td>-90.640</td><td>-0.04</td><td>-1.30</td></tr></tbody></table>

<table><caption>**Table 3**
Case (1) mesh convergence analysis ($α_2 = -8°$).</caption><thead><tr><th>FEA model</th><th>Number of elements</th><th>σ<sub>zz</sub> (kPa)</th><th>e<sub>i</sub> (%)</th><th>e<sub>meth</sub></th></tr></thead><tbody><tr><td rowspan="2">2D</td><td>2,096</td><td>-415.63</td><td>-</td><td>0.39</td></tr><tr><td>4,235</td><td>-415.68</td><td>0.014</td><td>0.38</td></tr><tr><td rowspan="2">3D</td><td>4,233</td><td>-415.63</td><td>-</td><td>0.48</td></tr><tr><td>9,546</td><td>-415.42</td><td>0.036</td><td>0.44</td></tr></tbody></table>

and $σ_{zz} = -417.258$ kPa ($α_2 = -8°$). In the skin region, the finite element solution converged with three elements per lamina.

Stresses on the xz-planes are given in Figs. 8 and 9, respectively, noting a state of plane stress has been invoked in the y-direction, reflecting the small dimension that exists in this coordinate. As expected, the $σ_{xx}$ magnitude is larger in the stiffer, rather than the more compliant 90° layers in the longitudinal direction.

The shear stress results of the developed formulation match with the finite element solution as depicted in Fig. 9(a). Conversely, the CLT methodology underestimates the shear stress magnitude in the skin region and it overestimates it in the core. Remarkably, tapering sandwich beams drastically changes the core functionality since the skin region experiences a larger transverse shear stress magnitude than the core. This observation is explained by the satisfaction of the stress boundary requirements discussed in Section 2.3.

Fig. 9(b) illustrates the results for $σ_{zz}$. Notably, the proposed method agrees well with the finite element solution. To better visualise the discrepancy between results, inset illustrations showing $σ_{zz}$ distributions are provided for the first layer. A small error is observed at the lower boundary between the 3D finite element analysis and the current formulation ($e_{meth} = -1.3\%$ for $α_1$ and $e_{meth} = -0.44\%$ for $α_2$). Furthermore, CLT does not capture the transverse direct stress distribution because $V' = 0$.

## 4.2. Asymmetrically tapered laminate

Consider the 10-layer tapered composite laminate beam depicted in Fig. 10 with one end fixed under a transverse pressure $t_z^{(1)} = -10x$ (N/mm) linearly distributed over the bottom surface as depicted in Fig. 10. The beam length is $L = 100$ mm and width $w = 1$ mm, and the lamina thickness varies from $h^{(k)}(0) = 1$ mm to $h^{(k)}(L) = 0.5$ mm, resulting in the beam height varying from $h(0) = 10$ mm to $h(L) = 5$ mm
<figure><img src="image_9.png" /><figcaption>Fig. 7. Case (1) finite element mesh: (a) $\alpha_1 = -4^\circ$; (b) $\alpha_2 = -8^\circ$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Case (1) $\sigma_{xx}$ results.</figcaption></figure>

<table><caption>**Table 4**
Case (2) force resultants and their derivatives for all layers.</caption><thead><tr><th>V<br>kN</th><th>M<br>kNm</th><th>V'<br>N/m</th></tr></thead><tbody><tr><td>-468.75</td><td>21.615</td><td>2.5</td></tr></tbody></table>

($\alpha_B = -2.86^\circ$). The elastic properties chosen for Case (2) are the same as in the former example given by Table 1. Results are provided for three different stacking sequences: (i) unidirectional [45]$_{10}$; (ii) crossply [0/90]$_5$; (iii) non-symmetric [0$_5$/90$_5$].

The interface surfaces for Case (2) are given by (in mm)

$$f_k = -0.05(k-1)x + (k-1) \quad (4.2)$$

Values for $\Gamma$ at $x = 0.25L$ are reported in Tables 4.

Element types used for the finite element meshes were CPS8R and C3D20R for the 2D and 3D finite element analyses, respectively. A side view of a typical finite element mesh for Case (2) is illustrated in Fig. 11. The convergence study is reported in Tables 5 to 7, where values for $\hat{\sigma}_{zz}$, i.e., the direct stress magnitude normalised by the magnitude of the pressure load at the tip $t_z^{(1)}(L) = 10$ N/mm, have been evaluated at the lower surface, noting that the proposed method resulted in $\hat{\sigma}_{zz} = 0.178$ ([0$_5$/90$_5$]), $\hat{\sigma}_{zz} = 0.0586$ ([0/90]$_5$) and $\hat{\sigma}_{zz} = 0.4235$ ([45]$_{10}$). The finite element solution converged with three elements per lamina.

Results for the asymmetric tapered beam under pressure load are discussed in this section. First, the longitudinal stress component is shown in Fig. 12. Subsequently, the transverse stresses $\tau_{xz}$ and $\sigma_{zz}$ are given in Figs. 13 and 14, respectively.

<table><caption>**Table 5**
Case (2) mesh convergence analysis for [0$_5$/90$_5$].</caption><thead><tr><th rowspan="2">Number of<br>elements<br>(×10<sup>3</sup>)</th><th colspan="3">FEA 2D</th><th colspan="3">FEA 3D</th></tr><tr><th>σ̄<sub>zz</sub><br>(-)</th><th>e<sub>i</sub><br>(%)</th><th>e<sub>meth</sub><br>(%)</th><th>σ̄<sub>zz</sub><br>(-)</th><th>e<sub>i</sub><br>(%)</th><th>e<sub>meth</sub><br>(%)</th></tr></thead><tbody><tr><td>1</td><td>0.1756</td><td>-</td><td>-1.51</td><td>0.4168</td><td>-</td><td>-1.56</td></tr><tr><td>4</td><td>0.1756</td><td>-0.53</td><td>-2.03</td><td>0.4146</td><td>-0.54</td><td>-2.08</td></tr><tr><td>12</td><td>0.1755</td><td>-0.09</td><td>-2.12</td><td>0.4143</td><td>-0.09</td><td>-2.17</td></tr></tbody></table>

<table><caption>**Table 6**
Case (2) mesh convergence analysis for [0/90]$_5$.</caption><thead><tr><th rowspan="2">Number of<br>elements<br>(×10<sup>3</sup>)</th><th colspan="3">FEA 2D</th><th colspan="3">FEA 3D</th></tr><tr><th>σ̄<sub>zz</sub><br>(-)</th><th>e<sub>i</sub><br>(%)</th><th>e<sub>meth</sub><br>(%)</th><th>σ̄<sub>zz</sub><br>(-)</th><th>e<sub>i</sub><br>(%)</th><th>e<sub>meth</sub><br>(%)</th></tr></thead><tbody><tr><td>1</td><td>0.0591</td><td>-</td><td>2.41</td><td>0.0591</td><td>-</td><td>2.40</td></tr><tr><td>4</td><td>0.0596</td><td>-0.57</td><td>1.83</td><td>0.0596</td><td>-0.56</td><td>1.83</td></tr><tr><td>12</td><td>0.0596</td><td>-0.09</td><td>1.73</td><td>0.0596</td><td>-0.10</td><td>1.72</td></tr></tbody></table>

<table><caption>**Table 7**
Case (2) mesh convergence analysis for [45]$_{10}$.</caption><thead><tr><th rowspan="2">Number of<br>elements<br>(×10<sup>3</sup>)</th><th colspan="3">FEA 2D</th><th colspan="3">FEA 3D</th></tr><tr><th>σ̄<sub>zz</sub><br>(-)</th><th>e<sub>i</sub><br>(%)</th><th>e<sub>meth</sub><br>(%)</th><th>σ̄<sub>zz</sub><br>(-)</th><th>e<sub>i</sub><br>(%)</th><th>e<sub>meth</sub><br>(%)</th></tr></thead><tbody><tr><td>1</td><td>0.4171</td><td>-</td><td>-0.83</td><td>0.4169</td><td>-</td><td>-0.84</td></tr><tr><td>4</td><td>0.4149</td><td>-0.50</td><td>-1.32</td><td>0.4146</td><td>-0.48</td><td>-1.32</td></tr><tr><td>12</td><td>0.4145</td><td>-0.08</td><td>-1.40</td><td>0.4142</td><td>-0.08</td><td>-1.40</td></tr></tbody></table>

The longitudinal stress predicted by the current method depicted in Fig. 12(a) agrees well with the finite element analyses. A small error
<figure><img src="image_11.png" /><figcaption>Fig. 9. Case (1) results: (a) $\tau_{xz}$; (b) $\sigma_{zz}$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Case (2) boundary conditions.</figcaption></figure>

between the proposed formulation and the finite element solutions is observed at the lower boundary, corresponding to 2.94% for the non-symmetric, 2.36% for the crossply, and 2.48% for the unidirectional laminates, which are in line with the expected error associated with the small taper angle assumption.

The results for the transverse shear stress $\tau_{xz}$ illustrated in Fig. 13 match well with the finite element analyses for all layups. The

maximum discrepancy obtained was 1.52% and 3.62% for the non-symmetrical and crossply laminates, respectively, located at the mid interface, and 2.30% for the unidirectional layup at the upper boundary. Remarkably, the discontinuity in the transverse shear stress was accurately captured by the developed theory. This observation confirms that the transverse stress field expressed in global coordinates is discontinuous through the thickness in tapered composite beams, as
<figure><img src="image_13.png" /><figcaption>Fig. 11. Case (2) finite element mesh.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Case (2) $\sigma_{xx}$ stress distribution.</figcaption></figure>

deduced in Eq. (2.8). Furthermore, the traction-free condition at the upper boundary was captured by the developed methodology but not by CLT, noting that the latter could not predict an appropriate shear stress profile.

The proposed method shows good levels of accuracy for $\sigma_{zz}$ when compared to the finite element analyses, with similar error magnitude of the transverse shear stress component. More precisely, the maximum discrepancy observed between finite elements and the developed theory results was 1.28% and 3.64% for the non-symmetrical and crossply laminates, respectively, located at the mid interface, and 2.23% for the unidirectional layup at the upper boundary. Overall, the [0/90]$_5$ profile was less accurate than the other stacking sequences, which indicates the accuracy of the developed theory is sensitive to laminate anisotropy due to non-symmetry. Notably, all methodologies captured the pressure load at the lower boundary, and the developed method matched the traction-free condition at the top boundary predicted by the finite element analyses. CLT obtained a different stress distribution than the benchmarked solution, underestimating $\sigma_{zz}$ magnitude for all stacking sequences with high levels of discrepancy in relation to finite element analysis results.

To verify the continuity of local stresses discussed in Section 2.3, the state of stress at $f_6$ (i.e. $\bar{z} = 0$) for layers $k = 5$ and $k = 6$ of [0$_5$/90$_5$] are rotated by $\alpha_{f_6} = -1.432^\circ$ as given in Table 8, which obtains the state of stress associated with a coordinate system parallel to $f_6$. Note that the stress components transverse to $f_6$, highlighted in light blue colour, have the same magnitudes for both layers, hence confirming the local transverse stress continuity (values for $\hat{\tau}_{xy}$ were calculated from the CLT and $\hat{\sigma}_{yy} = 0$).

## 5. Conclusions

An analytical methodology has been introduced to predict the stress field of untwisted composite laminated beams comprising orthotropic layers of linearly varying thickness. A state of plane stress with small and linearly varying strains, known internal forces, and

<table><caption>Table 8
Verification of local transverse stress continuity.</caption><thead><tr><th rowspan="2">Stress components</th><th colspan="2">k = 5</th><th colspan="2">k = 6</th></tr><tr><th>Global</th><th>Local</th><th>Global</th><th>Local</th></tr></thead><tbody><tr><td>σ̂<sub>xx</sub></td><td>32.712</td><td>32.749</td><td>2.183</td><td>2.201</td></tr><tr><td>τ̂<sub>xy</sub></td><td>0.653</td><td>0.653</td><td>0.651</td><td>0.651</td></tr><tr><td>τ̂<sub>xz</sub></td><td>-1.1429</td><td>-0.3252</td><td>-0.3796</td><td>-0.3252</td></tr><tr><td>τ̂<sub>yz</sub></td><td>0</td><td>0.0163</td><td>0</td><td>0.0163</td></tr><tr><td>σ̂<sub>zz</sub></td><td>0.0428</td><td>0.0061</td><td>0.0237</td><td>0.0060</td></tr></tbody></table>

perfectly bonded layers has been assumed. Additionally, body loads and pressure forces are prescribed in the transverse direction in a layerwise fashion. Within the scope of tapered composite beams, the main conclusions of this investigation include the following:

(1) transverse stresses are discontinuous at oblique layer interfaces whose elastic properties of adjacent layers are different;

(2) the maximum transverse stress magnitudes and their cross-sectional locations are not intuitive and require stress analysis to be determined;

(3) bending moment induces shear stress due to the taper effect;

(4) Classical Laminate Theory does not predict appropriate interfacial and surface equilibrium conditions of composite laminate tapered beams and underestimates transverse stress magnitudes.

The results suggest good agreement between the current formulation and 2D- and 3D- solid finite element analyses for non-symmetrical and symmetrical laminates. Regarding the error associated with the small angle restriction, satisfactory results have been obtained for up to 8° taper angle, noting that more significant cross-sectional variations may not be practicable for beam modelling due to the slenderness requirement of beam theory. Furthermore, the developed theory shows a substantial improvement in transverse stress field prediction when compared to Classical Laminate Theory.
<figure><img src="image_15.png" /><figcaption>Fig. 13. Case (2) results for $\tau_{xz}$: (a) $[0_{5}/90_{5}]$ and $[45]_{10}$; (b) $[0/90]_{5}$.</figcaption></figure>

The proposed theory is efficient from a computational standpoint as it does not require modelling the structure in commercial software, numerical convergence studies, and post processing of results. Practitioners can benefit from the low computational cost of the developed formulation and make fast and accurate early design stage decisions. However, it is advised that engineers assess the limitations for practical applications since the error of the introduced formulation escalates with the taper angle, material anisotropy, and mechanical response in the width direction. Potential further developments involve extending this work to thin-walled cross-sections, functionally graded materials and variable angle tow laminates.

## Data availability

No data was used for the research described in the article

## Acknowledgments

The authors would like to thank Science Foundation Ireland (SFI) for funding Spatially and Temporally VARIable COMPositive Structures (VARICOMP) Grant No. (15/RP/2773) under its Research Professor programme.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Appendix

Full expressions for the stiffness coefficients introduced in Eq. (3.13) are given in this section. The kth layer reference index is dropped for simplicity reasons ($Q_{ij}^{(k)} = Q_{ij}$).

### $\zeta$ Coefficients

$$ \zeta_V = Q_{1j} \left( \frac{1}{3} d'_j z^3 + b'_j z^2 \right) $$

$$ \zeta_m = -\zeta_V $$

$$ \zeta_M = Q_{1j} \left( \frac{1}{6} d''_j z^3 + \frac{1}{2} b''_j z^2 \right) $$

$$ \zeta_{V'} = \eta_{m'} = -Q_{1j} \left( \frac{1}{6} d_j z^3 + \frac{1}{2} b_j z^2 \right) $$

(A.1)
<figure><img src="image_16.png" /><figcaption>Fig. 14. Case (2) results for $\sigma_{zz}$: (a) [0$_{5}$/90$_{5}$] and [45$_{10}$]; (b) [0/90]$_{5}$.</figcaption></figure>

## λ Coefficients

$$ 
\begin{aligned}
\lambda_V &= -Q_{1j} (f^2 d'_j + f (d_j f' + 2b'_j) + f' b_j) \\
\lambda_m &= -\lambda_V \\
\lambda_M &= -Q_{1j} \left[ \frac{f^2}{2} d''_j + f (b''_j + d'_j f') + f' b'_j \right] \\
\lambda_{V'} &= \eta_{m'} = Q_{1j} \frac{f}{2} (2b_j + f d_j)
\end{aligned}
\qquad (A.2) 
$$

## References

[1] Baluzzi G, Aminbaghai M, Sacco E, Füssl J, Eberhardsteiner J, Auricchio F. Non-prismatic beams: a simple and effective Timoshenko-like model. Int J Solids Struct 2016;90:236–50.

[2] Baluzzi G, Morganti S, Auricchio F, Reali A. Non-prismatic Timoshenko-like beam model: Numerical solution via isogeometric collocation. Comput Math Appl 2017;74(7):1531–41.

[3] Baluzzi G, Hochreiner G, Füssl J. Stress recovery from one dimensional models for tapered bi-symmetric thin-walled I beams: Deficiencies in modern engineering tools and procedures. Thin-Walled Struct 2017;119:934–45.

[4] Baluzzi G, Morganti S, Füssl J, Aminbaghai M, Reali A, Auricchio F. Modeling the non-trivial behavior of anisotropic beams: a simple Timoshenko beam with enhanced stress recovery and constitutive relations. Compos Struct 2019;229:111265.

[5] Ojo SO, Weaver PM. Efficient strong unified formulation for stress analysis of non-prismatic beam structures. Compos Struct 2021;272:114190.

[6] Taglialegne L. Analytical study of stress fields in wind turbine blades [Ph.D. thesis], Architecture, Civil Engineering and Environmental Engineering. Universities of Florence, Perugia and Pisa - TU Braunschweig; 2018, p. 1–181.

[7] Bertolini P, Eder M, Taglialegne L, Valvo P. Stresses in constant tapered beams with thin-walled rectangular and circular cross sections. Thin-Walled Struct 2019;137:527–40.

[8] Bertolini P, Taglialegne L. Analytical solution of the stresses in doubly tapered box girders. Eur J Mech A Solids 2020;81:103969.

[9] Mercuri V, Baluzzi G, Asprone D, Auricchio F. Structural analysis of non-prismatic beams: Critical issues, accurate stress recovery, and analytical definition of the finite element (FE) stiffness matrix. Eng Struct 2020;213:110252.

[10] Beltempo A, Baluzzi G, Alfano G, Auricchio F. Analytical derivation of a general 2D non-prismatic beam model based on the Hellinger-Reissner principle. Eng Struct 2015;101:88–98.

[11] Auricchio F, Baluzzi G, Lovadina C. A new modeling approach for planar beams: Finite-element solutions based on mixed variational derivations. J Mech Mater Struct 2010;5(5):771–94.

[12] Hodges DH, Ho J, Yu W. The effect of taper on section constants for in-plane deformation of an isotropic strip. J Mech Mater Struct 2008;3(3):425–40.

[13] Hodges DH, Rajagopal A, Ho J, Yu W. Stress and strain recovery for the in-plane deformation of an isotropic tapered strip-beam. J Mech Mater Struct 2011;5(6):963–75.

[14] Masjedi PK, Maheri A, Weaver PM. Large deflection of functionally graded porous beams based on a geometrically exact theory with a fully intrinsic formulation. Appl Math Model 2019;76:938–57.

[15] Masjedi PK, Weaver PM. Analytical solution for the fully coupled static response of variable stiffness composite beams. Appl Math Model 2020;81:16–36.

[16] Masjedi PK, Weaver PM. Variable stiffness composite beams subject to non-uniformly distributed loads: An analytical solution. Compos Struct 2020;112975.

[17] Doeva O, Masjedi PK, Weaver PM. Static analysis of composite beams on variable stiffness elastic foundations by the homotopy analysis method. Acta Mech 2021;232:4169–88.

[18] Viglietti A, Zappino E, Carrera E. Analysis of variable angle tow composites structures using variable kinematic models. Composites B 2019;171:272–83.
[19] Viglietti A, Zappino E, Carrera E. Free vibration analysis of variable angle-tow composite wing structures. Aerosp Sci Technol 2019;92:114-25.

[20] Kemp B, Johnson E. Response and failure analysis of a graphite-epoxy laminate containing terminating internal plies. In: 26th structures, structural dynamics, and materials conference. 1985, p. 608.

[21] Curry JM, Johnson ER, Starnes Jr JH. Effect of dropped plies on the strength of graphite-epoxy laminates. AIAA J 1992;30(2):449-56.

[22] Mukherjee A, Varughese B. Design guidelines for ply drop-off in laminated composite structures. Composites B 2001;32(2):153-64.

[23] Varughese B, Mukherjee A. A ply drop-off element for analysis of tapered laminated composites. Compos Struct 1997;39(1-2):123-44.

[24] Edwin Sudhagar P, Ananda Babu A, Rajamohan V, Jeyaraj P. Structural optimization of rotating tapered laminated thick composite plates with ply drop-offs. Int J Mech Mater Des 2017;13(1):85-124.

[25] Clancy G, Peeters D, O'Higgins RM, Weaver PM. In-line variable spreading of carbon fibre/thermoplastic pre-preg tapes for application in automatic tape placement. Mater Des 2020;194:108967.

[26] Clancy GJ. Manufacture of variable angle tow laminates using carbon fibre/thermoplastic pre-preg tapes [Ph.D. thesis], University of Limerick; 2021.

[27] Vilar M, Hadjiloizi D, Masjedi PK, Weaver PM. Stress analysis of generally asymmetric non-prismatic beams subject to arbitrary loads. Eur J Mech A Solids 2021;104284.

[28] Vilar M, Hadjiloizi D, Khaneh Masjedi P, Weaver P. Stress recovery of laminated non-prismatic beams under layerwise traction and body forces. Int J Mech Mater Des 2022;18(3):719-41.

[29] Vilar M, Masjedi PK, Hadjiloizi DA, Weaver PM. Analytical plane-stress recovery of non-prismatic beams under partial cross-sectional loads and surface forces. Eng Struct 2022;252:113169.

[30] Vu-Quoc L, Léger P. Efficient evaluation of the flexibility of tapered I-beams accounting for shear deformations. Internat J Numer Methods Engng 1992;33(3):553-66.

[31] Bleich F. Stahlhochbauten: Ihre theorie, berechnung und bauliche gestaltung (1. Band). Springer-Verlag; 1932.

[32] Balduzzi G, Aminbaghai M, Auricchio F, Füssl J. Planar Timoshenko-like model for multilayer non-prismatic beams. Int J Mech Mater Des 2018;14(1):51-70.

[33] Zhou M, Fu H, An L. Distribution and properties of shear stress in elastic beams with variable cross section: Theoretical analysis and finite element modelling. KSCE J Civ Eng 2020;1-15.

[34] Zhou M, Zhang J, Zhong J, Zhao Y. Shear stress calculation and distribution in variable cross sections of box girders with corrugated steel webs. J Struct Eng 2016;142(6):1-10.

[35] Boley BA. On the accuracy of the Bernoulli-Euler theory for beams of variable section. J. Appl Mech 1963;373-8.

[36] Jourawski D. Sur la résistance d'un corps prismatique et d'une pièce composée en bois ou en tôle de fer à une force perpendiculaire à leur longueur. In: Annales des Ponts et Chaussées, vol. 12, 1856, p. 328-51.

[37] Michell J. The stress in an elotropic elastic solid with an infinite plane boundary. Proc Lond Math Soc 1900;1(1):247-57.

[38] Carothers S. XXVI.—Plane strain in a wedge, with applications to masonry dams. Proc Roy Soc Edinburgh 1914;33:292-306.

[39] Timoshenko SP, Goodier JN. Theory of elasticity. McGraw-Hill; 1951.

[40] Krahula JL. Shear formula for beams of variable cross section. AIAA J 1975;13(10):1390-1.

[41] Knops R, Villaggio P. Recovery of stresses in a beam from those in a cone. J Elasticity 1998;53(1):65-75.

[42] Blodgett OW. Design of welded structures. 1966, Cleveland: James F. Lincoln Arc Welding Foundation.

[43] Trahair N, Ansourian P. In-plane behaviour of web-tapered beams. Eng Struct 2016;108:47-52.

[44] Romano F. Deflections of Timoshenko beam with varying cross-section. Int J Mech Sci 1996;38(8-9):1017-35.

[45] Migliaccio G, Ruta G, Barsotti R, Bennati S. A new shear formula for tapered beamlike solids undergoing large displacements. Meccanica 2022;1-22.

[46] Chockalingam SN, Nithyadharan M, Pandurangan V. Shear stress distribution in tapered I-beams: Analytical expression and finite element validation. Thin-Walled Struct 2020;157:107152.

[47] Chockalingam SN, Pandurangan V, Nithyadharan M. Timoshenko beam formulation for in-plane behaviour of tapered monosymmetric I-beams: Analytical solution and exact stiffness matrix. Thin-Walled Struct 2021;162:107604.

[48] Gayen D, Roy T. Hygro-thermal effects on stress analysis of tapered laminated composite beam. Int J Compos Mater 2013;3(3):46-55.

[49] Vel SS, Caccese V, Zhao H. Elastic coupling effects in tapered sandwich panels with laminated anisotropic composite facings. J Compos Mater 2005;39(24):2161-83.

[50] Ai Q, Weaver PM. Simplified analytical model for tapered sandwich beams using variable stiffness materials. J Sandw Struct Mater 2017;19(1):3-25.

[51] Zhou C, Li L, Wang J. Modified bar simulation method for shear lag analysis of non-prismatic composite box girders with corrugated steel webs. Thin-Walled Struct 2020;155:106957.

[52] Zhou M, Liao J, Zhong J, An L, Wang H. Unified calculation formula for predicting the shear stresses in prismatic and non-prismatic beams with corrugated steel webs. In: Structures, vol. 29. Elsevier; 2021, p. 507-18.

[53] Marín J, Graciani E. Normal stress flow evaluation in composite aircraft wing sections by strength of material models. Compos Struct 2022;282:115088.

[54] Rajagopal A, Hodges DH. Asymptotic approach to oblique cross-sectional analysis of beams. J Appl Mech 2014;81(3).

[55] Balduzzi G, Hochreiner G, Füssl J, Auricchio F. Serviceability analysis of non-prismatic timber beams: derivation and validation of new and effective straightforward formulas. Open J Civ Eng 2017;7(1):32-62.

[56] Reddy JN. Mechanics of laminated composite plates and shells: Theory and analysis. CRC Press; 2003.

[57] Jones RM. Mechanics of composite materials. CRC Press; 2018.

[58] Daniel IM, Ishai O, Daniel IM, Daniel I. Engineering mechanics of composite materials, vol. 1994. Oxford university press New York; 2006.

[59] Pagano N. Stress fields in composite laminates. Int J Solids Struct 1978;14(5):385-400.

[60] Bruhns O. Advanced mechanics of solids. Springer; 2003.