An Efficient Quasi-3D Theory for Cylindrical Bending of Laminated Composite Plate with Thickness Stretching Effects

The composite laminates are widely adopted in engineering applications because of their high strength-to-weight ratio and ability to tailor stiffness through stacking sequences. Asymmetric configurations such as (0^°/90^°) are particularly attractive since they balance stiffness and enhance structural performance. Predicting bending and stress responses in such a laminate composite is challenging, especially under mechanical and thermal loads. The traditional theory often fails for thick laminates because it neglects transverse shear and normal effects, which has motivated the development of refined higher-order and quasi-3D formulations.

Several advanced theories have been proposed earlier to address these limitations. Early examined vibration and dynamic behaviour of composite plate [1, 2], while Sayyad and Ghugal [3] introduced sinusoidal shear and normal deformation plate theory (SSNPT) to incorporate transverse normal strain effects under mechanical loading. Quasi-3D and higher order models have been applied to static, vibration, and buckling problems [4-6] and thermomechanical analysis of shells and plates [7, 8]. An analytical model for bending and vibration of thick advanced plates using a four-variable quasi-3D higher-order theory has been developed by Khiloun et al. [9]. Nonlinear thermal effects on cylindrical bending have also been studied and presented by Kulkarni [10]. Review of shear deformation theories highlights the progression from classical to advanced quasi-3D formulations [11]. Higher-order displacement models have been extended to environmental loading scenarios [12], and Pagano [13] provided exact elasticity solutions that remain a benchmark for validating approximate theories. More recent quasi-3D approaches have been applied to functionally graded beams and sandwich structures [14]. Latest studies further extend quasi-3D and refined shear deformation theories to nonlinear thermomechanical bending, porous laminates and micro-scale plate structures [15-18]. A higher order formulation to evaluate the thermal response of an orthotropic plate subjected to cylindrical bending conditions with improved accuracy was developed by Kulkarni [19]. Among the foundational works in composite mechanics, Reddy [20] provides a comprehensive treatment of laminated plates and shells theories. The book systematically develops constitutive relations, bending behaviour and thermomechanical effects, making it a standard reference for researchers in composite structures. Sayyad et. al. [21] developed a quasi-3D higher-order model considering transverse normal and shear effects for composite plates; however, their work mainly addressed general flexural response and did not specifically focus on cylindrical bending. Gwak et al. [22] proposed an enhanced shear deformation finite element formulation for thermo-mechanical analysis, but their approach is primarily numerical and does not provide analytical benchmark solutions for validation. Junjie et al. [23] investigated the thermomechanical behavior of composite cylindrical shells considering temperature-dependent material characteristics under combined thermal and mechanical loading conditions.

Although many of these formulations include transverse shear and thickness stretching effects, most focus is noted on mechanical loading, vibration and buckling. Their application to thermal cylindrical bending of asymmetric laminates remains limited. In particular, the role of transverse normal effect under thermal loads has not been comprehensively addressed, leaving a gap in the literature.

The present work addresses this gap by applying a quasi-3D formulation, refined parabolic shear deformation theory (PSDT) and SSNPT to analyze the bending and stress response of an asymmetric laminate subjected to thermal loading. By explicitly incorporating transverse normal effects, the study aims to provide accurate predictions of displacements and stresses, establishing a reliable framework for thick laminates where conventional theories lose accuracy.

In the present study, the PSDT used is an existing refined theory and is not newly developed. The main contribution of this work is the application and detailed comparison of quasi-3D and PSDT formulations for thermal and mechanical analysis of a two-layer asymmetric laminated composite plate under cylindrical bending. A systematic comparison is carried out with available elasticity solutions and other higher theories. The study provides a consistent set of results for displacements and stress components under both mechanical and thermal loading. These results can serve as a useful reference for the validation of future models. Special attention is given to the role of transverse normal effects and their influence on the accuracy of prediction, particularly for thick laminates.

The quasi-3D theory, as developed by Vo et al. [14], is an advanced structural modelling approach designed to overcome the limitations of traditional theory. The present quasi-3D formulation captures a realistic distribution of transverse shear strains and also permits thickness-direction stretching, such that $\left(\varepsilon_z \neq 0\right)$. This makes it effective for analyzing laminated composite structures, where the through-thickness effect plays a significant role. The theory introduces displacement fields that incorporate higher-order terms in the thickness coordinate, ensuring accurate representation of both bending and shear behaviour of the plate subjected to cylindrical bending conditions. From the above assumptions, the displacement field according to the quasi-3D formulation can be written as:

$u(x, z)=u_0-z \frac{\partial w_b(x)}{\partial x}-\frac{4}{3}\left(\frac{z^3}{h^2}\right) \frac{\partial w_s(x)}{\partial x}$       (1)

$w(x, z)=w_b(x)+w_s(x)+\left(1-4 \frac{z^2}{h^2}\right) w_z(x)$        (2)

In this formulation, $u(x, z)$ denotes the axial displacement along x axis, whereas $w(x, z)$ corresponds to the transverse displacement in the z-direction. Importantly, stretching in the thickness direction is permitted, reflecting the condition $\left(\varepsilon_z \neq 0\right)$.

3.1 Strain-displacement relationship

The strain-displacement relations in linear elasticity describe how the deformation of a body is related to the displacement of its points. According to this theory, the normal strain in the $x$ direction $\left(\varepsilon_x\right)$, the normal strain in the $z$ direction $\left(\varepsilon_z\right)$ and the shear strain $\left(\tau_{x z}\right)$ are obtained as follows.

$\varepsilon_x=\frac{\partial u(x, z)}{\partial x}$       (3)

$\varepsilon_z=\frac{\partial w(x, z)}{\partial z}$       (4)

$\tau_{x z}=\frac{\partial u(x, z)}{\partial z}+\frac{\partial w(x, z)}{\partial x}$       (5)

3.2 Constitutive relations

The constitutive relations define how stresses and strains are mathematically linked in a material under given loading and environmental conditions. The stress-strain relationship for a laminated composite plate in cylindrical bending is expressed by Eq. (6), which represents a reduced constitutive relation derived under the assumption of cylindrical bending of a laminated composite plate. In this case, only the stress components $\sigma_x, \sigma_z, \tau_{x z}$ and the corresponding strains are active, leading to a 3 × 3 constitutive matrix instead of the full 6 × 6 form. The reduced stiffness coefficients $\bar{Q}_{i j}$ are obtained from the transformed stiffness matrix of the laminate following classical laminate theory. Thermal expansion effects are included in the x and z directions. This reduction is consistent with the assumptions of plain strain in the x-z plane and neglect of stresses in the y direction as commonly adopted in laminated plate bending analysis [21].

$\left(\begin{array}{c}\sigma_x \\ \sigma_z \\ \tau_{x z}\end{array}\right)=\left[\begin{array}{ccc}\bar{Q}_{11} & \bar{Q}_{13} & 0 \\ \bar{Q}_{13} & \bar{Q}_{33} & 0 \\ 0 & 0 & \bar{Q}_{55}\end{array}\right]\left(\begin{array}{c}\varepsilon_x-\alpha_x \Delta T \\ \varepsilon_z-\alpha_z \Delta T \\ \gamma_{x_z}\end{array}\right)$                           (6)

The left-hand side of Eq. (6) represents stress components. The right-hand side of the equation consists of strain components with thermal expansion terms in x and z directions multiplied by reduced elastic constants $\left(\bar{Q}_{i j}\right)$. The reduced elastic constants are as given below.

$\begin{gathered}\bar{Q}_{11}=\frac{E_1\left(1-\mu_{23} \mu_{32}\right)}{\Delta} \\ \bar{Q}_{13}=\frac{E_1\left(\mu_{31}+\mu_{21} \mu_{32}\right)}{\Delta} \\ \bar{Q}_{33}=\frac{E_3\left(1-\mu_{12} \mu_{21}\right)}{\Delta} \\ \bar{Q}_{55}=G_{13}\end{gathered}$                       (7)

where,

$\Delta=\left[1-\mu_{12} \mu_{21}-\mu_{23} \mu_{32}-\mu_{31} \mu_{13}-2 \mu_{12} \mu_{23} \mu_{21}\right]$           (8)

The Young's moduli along $x$ (direction 1) and $z$ (direction 3) are represented by $E_1$ and $E_3$ respectively. They measure resistance to normal deformation. The shear modulus in the $x$-$z$ plane is denoted by $G_{13}$. This represents resistance to shear deformation. The Poisson ratios describing how strain in one direction induces strain in another perpendicular direction are denoted by $\mu_{12}, \mu_{21}, \mu_{13}, \mu_{23}, \mu_{31}, \mu_{32}$.

3.3 Temperature distribution

For a laminated composite plate undergoing cylindrical bending, the temperature difference through the thickness is assumed to vary linearly and is given by Eq. (9).

$\Delta T=\frac{2 z}{h} T_1(x)$       (9)

where, $\Delta T$ represents the temperature change measured from a chosen reference state. The temperature varies with both the in-plane coordinate x and the thickness coordinate z. The term $T_1(x)$ denotes the thermal loading, which varies along the length of the plate, while the temperature distribution remains linear through the thickness of the laminated composite plate in cylindrical bending conditions. In this work, the temperature field is directly prescribed and is not obtained by solving a heat conduction equation or problem. This assumption is used to simplify the analysis and to compare the present results with earlier or available studies on laminated composite plates under cylindrical bending. The selected or chosen form represents a linear temperature change across the thickness, which varies along the length of the plate through $T_1(x)$. However, this does not correspond to any specific physical thermal boundary condition, such as applied surface temperature or heat transfer. Hence, the results of this study are valid only for the assumed temperature distribution. A more realistic analysis would require solving the heat conduction equation to find the actual temperature variation, which can be considered in future work. The linear variation of temperature across the thickness is shown in Figure 2. 

2.png

Figure 2. Linear variation of temperature across the thickness of a laminated composite plate subjected to cylindrical bending

Navier’s solution for simply supported laminated composite plates assumes trigonometric series expressions that satisfy boundary conditions along $(x=0)$ and $(x=a)$, the constraints are $w=0, M_x=0, N_x=0, M_x^s=0$, ensuring zero deflection, vanishing bending, in-plane force and shear moments. Here, $M_x^s$ denotes the moment resultant associated with the transverse shear deformation as defined in the adopted higher-order theory. Further, the problem is treated under a bending condition, where all field variables vary only along the x direction and remain constant along the y direction. This corresponds to an infinitely long plate with no variation in the transverse in-plane direction, thereby reducing the problem to a one-dimensional form.

$u_0(x)=\sum_{\infty}^{m=1} u_m \cos \left(\frac{m \pi x}{a}\right)$           (18)

$w_b(x)=\sum_{\infty}^{m=1} W_{b m} \sin \left(\frac{m \pi x}{a}\right)$        (19)

$w_s(x)=\sum_{\infty}^{m=1} W_{s m} \sin \left(\frac{m \pi x}{a}\right)$       (20)

$w_z(x)=\sum_{\infty}^{m=1} W_{z m} \sin \left(\frac{m \pi x}{a}\right)$     (21)

The thermal and mechanical loads are expanded in a single Fourier sine series as given below.

$T_1(x)=\sum_{\infty}^{m=1} T_{1 m} \sin \left(\frac{m \pi x}{a}\right)$               (22)

$q(x)=\sum_{m=1}^{\infty} q_m \sin \left(\frac{m \pi x}{a}\right)$          (23)

where, $q_m=q_0=1, T_{1 m}=T_0=1$ for sinusoidally distributed mechanical and thermal loads.

$q_m=\frac{4 q_0}{m \pi}$ and $T_{1 m}=\frac{4 T_0}{m \pi}$ for uniformly distributed mechanical and thermal loads. By substituting Eqs. (18)-(23) into the governing Eqs. (11)-(14), the system reduces to a set of coupled algebraic expressions. These can be compactly represented in matrix form, facilitating solution through standard linear algebra techniques and represented by the following Eq. (24).

$[K]\{X\}=\{f\}$       (24)

The stiffness matrix, unknowns and force vector are represented by $[K],\{X\}$ and $\{f\}$ respectively in the above Eq. (24). The stiffness matrix obtained and its elements are written as below.

$[K]=\left[\begin{array}{llll}k_{11} & k_{12} & k_{13} & k_{14} \\ k_{21} & k_{22} & k_{23} & k_{24} \\ k_{31} & k_{32} & k_{33} & k_{34} \\ k_{41} & k_{42} & k_{43} & k_{44}\end{array}\right]$            (25)

The elements of $[K]$ matrix are written below.

$\begin{gathered}k_{11}=A_{11} \alpha^2, k_{12}=-B_{11} \alpha^3, \\ k_{13}=-\frac{4}{3 h^2} B_{11} \alpha^3, k_{14}=\frac{8}{h^2} B_{13} \alpha, \\ k_{21}=k_{12}, k_{22}=D_{11} \alpha^4, k_{23}=\frac{4}{3 h^2} F_{11} \alpha^4, \\ k_{24}=-\frac{8}{h^2} D_{13} \alpha^2 \\ k_{31}=k_{13}, k_{32}=k_{23}, \\ k_{33}=\frac{16}{9 h^4} H_{11} \alpha^4+\left(A_{55}-\frac{8}{h^2} D_{55}+\frac{16}{h^4} F_{55}\right) \alpha^2 \\ k_{34}=\left(-\frac{32}{9 h^4} F_{13}+A_{55}-\frac{8}{h^2} D_{55}+\frac{16}{h^4} F_{55}\right) \alpha^2 \\ k_{41}=k_{14}, k_{42}=k_{24}, k_{43}=k_{34}, \\ k_{44}=\frac{64}{h^4} D_{33}+\left(A_{55}-\frac{8}{h^2} D_{55}+\frac{16}{h^4} F_{55}\right) \alpha^2  \end{gathered}$                 (26)

The unknowns to be determined and force matrix elements are as below.

$\{X\}=\left\{\begin{array}{c}u_m \\ w_{b m} \\ w_{s m} \\ w_{z m}\end{array}\right\}$           (27)

$\{f\}=\left\{\begin{array}{c}f_1=-\frac{2}{h} \alpha\left(T B_{11}+T B_{13}\right) T_{1 m} \\ f_2=\frac{2}{h} \alpha^2\left(T D_{11}+T D_{13}\_z\right) T_{1 m} \\ f_3=\frac{8}{3 h^3} \alpha^2\left(T F_{11}+T F_{13}\right) T_{1 m} \\ f_4=-\frac{16}{h^3}\left(T D_{13 \_x}+T D_{33}\right) T_{1 m}\end{array}\right\}$             (28)

In the above Eq. (26), $\alpha=\frac{m \pi}{a}$. The unkowns $\{X\}$ can be evaluated by solving above set of algebric equations. Further, displacements and stresses can be obtained under thermal and mechanical loads.

The transverse shear stress in elasticity can be obtained by applying the shear stress-free boundary conditions at the top and bottom surfaces and integrating the equilibrium equation for $\tau_{z x}$ with respect the thickness coordinate. This leads to the following expression.

$\tau_{x z}^{E E}=\int_{-h / 2}^{Z_{k+1}} \frac{\partial \sigma_{x x}}{\partial x} d z+C_1$         (29)

where, $\sigma_{xx}$ is normal stress in the $x$ direction, $h$ is the total thickness, $z_{k+1}$ represents the upper limit of integration along the thickness coordinate and $C_1$ is the constant of integration determined from the boundary conditions.

5.1 Numerical implementation

The governing equations are solved using Navier’s method by expressing the displacement variables in the form of a trigonometric series that satisfy the simply supported conditions. For the considered loading condition, the solution is dominated by the fundamental term $(m=1)$, which satisfied the simply supported boundary conditions and represented the primary deformation mode. It was observed that the inclusion of higher order terms does not produce any significant change in the results. Therefore, only the first is retained in the present computation.

The constant of integration $\left(C_1\right)$ in Eq. (29) is determined by enforcing the transverse shear stress-free boundary condition at the plate surfaces, i.e., $\tau_{x z}=0$ at $z= \pm \frac{h}{2}$. This condition is applied after performing the thickness-wise integration of the equilibrium equation, ensuring that the computed shear stress distribution satisfies the physical boundary requirements.

For laminated plates, the integration through the thickness is carried out layer-wise. The stress components are first evaluated within each layer using the corresponding material properties and continuity of stresses is maintained at the layer interfaces. The total response is then obtained by assembling the contributions from all layers.

The complete formulation and solution procedure are implemented in MATLAB software. Standard numerical routines are used to evaluate derivatives and perform integrations involved in the formulation.

7.1 Response of the laminate under mechanical loading

The results of stresses and displacements under sinusoidal and uniform mechanical load obtained using quasi-3D and PSDT are shown in Table 2.

Table 2. Non-dimensional displacements and stress components of a two-layer (0°/90°) asymmetric laminated plate subjected to cylindrical bending under mechanical loading

Note: PSDT = parabolic shear deformation theory

The results evaluated by quasi-3D and PSDT are very close to each other for both thick and thin laminates in cylindrical bending.

These results are also close to the exact elasticity solution and the results evaluated by SSNPT available in the literature. At low aspect ratio (a/h = 4), these models predict displacement and stresses that nearly match the exact elasticity solution, confirming their reliability. When the aspect ratio (a/h) changes to 10, the refined theories continue to align well with the exact elasticity solution.

The variation of the normalized in-plane or axial displacement $\bar{u}$, normal stress $\bar{\sigma}_{xx}$ and transverse shear stress $\bar{\tau}_{x z}$ across the thickness of the plate under sinusoidal and uniform mechanical load is shown in Figures 3-5, respectively. The variation of axial displacement through the thickness shows distinct differences between sinusoidal and uniform loadings. Under sinusoidal load, the displacement distribution is smoother and more gradual, reflecting the periodic nature of the applied force. Whereas a uniform load produces a more consistent displacement profile with noticeable change near the mid-plane (Figure 3).

3.png

Figure 3. Axial displacement ($\bar{u}$) under sinusoidal and uniform mechanical load for an aspect ratio of 4

This highlights how the type of mechanical load strongly influences the deformation pattern of the laminate in cylindrical bending. The normal stresses (Figure 4) vary significantly across the thickness. Figure 4 indicates that sinusoidal loading introduces localized stress variation, while uniform loading distributes stresses more evenly across the laminates.

4.png

Figure 4. Normal stress ($\bar{\sigma}_{xx}$) under sinusoidal and uniform mechanical load for an aspect ratio of 4

5.png

Figure 5. Transverse shear stress ($\bar{\tau}_{x z}$) under sinusoidal and uniform mechanical load for an aspect ratio of 4

The transverse shear stress shows strong dependence on the type of loading (Figure 5). In the case of sinusoidal loading, the shear stress curve is more oscillatory, whereas uniform loading produces a smoother distribution with maximum shear stress occurring at the mid-plane. This demonstrates that sinusoidal loads amplify shear variations across the thickness of the plate in cylindrical bending, whereas uniform loads give a more predictable shear response.

In addition, the quasi-3D formulation includes transverse normal strain ($\varepsilon_z$), which introduces thickness stretching and coupling with transverse normal stress ($\sigma_z$). Although this effect is small under mechanical loading in cylindrical bending, it slightly improves the accuracy of stress prediction compared to PSDT, which neglects this contribution.

7.2 Response of the laminate under thermal loading

For the thermal loading case, results from the exact elasticity solution reported by Pagano [13] are not included in Table 3, as thermal benchmark data for the (0^°/90^°) configuration are not available in that reference. Therefore, validation of the present results is carried out through comparison with other refined plate theories, such as PSDT and SSNPT, which are well established in the literature. The results are obtained by using quasi-3D, PSDT and SSNPT theories and compared with each other for validation purposes.

Table 3. Non-dimensional displacements and stress components of a two-layer (0°/90°) asymmetric laminated plate subjected to cylindrical bending under thermal loading

The results of thermal stresses and displacement under sinusoidal thermal load evaluated using quasi-3D, parabolic (PSDT) and sinusoidal shear and normal plate theory (SSNPT) are shown in Table 3. In the case of a lower aspect ratio (a/h = 4), the quasi-3D and PSDT models give identical predictions for displacement and stresses, showing that both theories capture the thermal response very well.

A slight difference is noted in the values evaluated by the SSNPT model. The results of all three refined theories remain close to each other, confirming their accuracy in representing thermal effects in laminated plates in cylindrical bending. The transverse shear stresses drop significantly when the aspect ratio changes to 10. This indicates the reduced influence of shear effect as the plate becomes thinner.

The variation of the normalized in-plane or axial displacement $\bar{u}$, normal stress $\bar{\sigma}_{xx}$ and transverse shear stress $\bar{\tau}_{x z}$ through the thickness of the plate in cylindrical bending under sinusoidal thermal load is shown in Figures 6-8, respectively. The axial displacement (Figure 6) occurs near the outer surfaces, while the mid-plane experiences a smaller value. The normal stress changes significantly through the thickness (Figure 7).

6.png

Figure 6. Axial displacement ($\bar{u}$) under sinusoidal thermal load for an aspect ratio of 4

7.png

Figure 7. Normal stresses ($\bar{\sigma}_{xx}$) under sinusoidal thermal load for an aspect ratio of 4

8.png

Figure 8. Transverse shear stress ($\bar{\tau}_{x z}$) under sinusoidal thermal load for an aspect ratio of 4

Under sinusoidal thermal load, the normal stress distribution is asymmetric with higher values near the outer layers. This reflects the strong effect of thermal gradients on the laminates under cylindrical bending conditions. The transverse shear stress profile shows peaks closer to the mid-plane and decreasing towards the surfaces.

The identical results obtained using the quasi-3D and PSDT theories are due to the present cylindrical bending condition and the assumed thermal loading; both theories reduce to a similar form in terms of displacement field representation. In particular, the contribution of thickness stretching and higher order terms in the quasi-3D model becomes negligible for the selected laminate configuration and loading case. As a result, both models predict nearly the same response for in-plane displacement, transverse displacement and normal stress. However, a difference between the two theories is observed in the transverse shear values, especially under uniform thermal loading. This is because the quasi-3D theory accounts for transverse normal deformation. Further, it is observed that transverse shear stress is more sensitive to the type of thermal loading and plate thickness. For thicker plate (lower a/h), shear stresses are high due to stronger through-thickness interaction. As the plate becomes thinner (a/h = 10), the shear stress reduces significantly, indicating that bending behaviour dominates and shear effects become less.

The close agreement between quasi-3D and PSDT results is limited to the present loading and boundary conditions and may not hold for cases involving strong thickness stretching effects or different thermal boundary conditions.

Under thermal loading, the effect of transverse normal strain ($\varepsilon_z$) becomes more relevant due to temperature variation across the thickness. The quasi-3D model captures this through thickness stretching while PSDT neglects it. However, for present loading, this contribution remains small for displacement and normal stress, leading to identical results. Its influence is more noticeable in transverse shear stress due to coupling between thermal gradients and through-thickness deformation.