Thermoelastic stress analysis in composite cylindrical vessel with metallic liner using first-order shear deformation theory and differential quadrature method

Abstract

In this study, thermoelastic stress analysis in composite metal cylindrical vessel is investigated with the help of two theories, first-order shear deformation theory (FSDT) and layerwise theory, using differential quadrature method (DQM). Hamilton principle and FSDT are employed for derivation of the principle differential equations. DQM and Newmark’s methods are used for stress analysis of a composite metal cylindrical vessel subjected to internal dynamic pressure and thermal load achieved from experimental test. In order to more accurately estimate the thickness effect, the layerwise theory is applied to approximate the temperature field. In this solution method, DQM is used to discretize the resulted governing equations while Newmark’s method is implemented to solve the problem in the time domain. It was observed that the fast convergence rate of the DQM method and its accuracy is confirmed by comparing the results with those obtained using ABAQUS software and also with the exact solutions. Furthermore, it was achieved significantly less computational efforts of the proposed approaches with respect to the finite element method.

Introduction

Composite materials are becoming an important part of engineering materials because they present advantages as higher specific stiffness and strength, improved corrosion resistance, and better fatigue strength compared to other materials. The best structures in different industries such as chemical, petroleum industries, and so on are thick-walled cylinders that are used as pressure vessel and cylindrical shells. In advanced industries, the application of new materials such as composite and functionally graded materials (FGM) are increasing due to their high usage in controlling stresses due to mechanical and thermal loads. Ordinary pressure vessels made of single layer may very often satisfy the service requirement subject to high temperature, high pressure, and corrosive environments. Therefore, a multilayer composite pressure vessel is widely used to meet the requirements under various industrial service condition by using different layers of different material. Generally, the inner layer is made of high performance alloy material because the inner layer is subjected to high pressure and high temperature while the outer layer is made of steel or fiber. In aerospace structures such as multilayer composite plates and shells, the temperature changes are one of the most important factor for the stress fields that cause their failure. Composite materials of high efficiency combine a number of properties, including high specific stiffness and strength, and nearly zero coefficient of thermal expansion in the direction of fiber. These properties result in a growing use of composite materials in structures under extreme thermal environment, such as high temperatures, high gradients, and temperature cycling changes. Due to these implications, considerable works on laminated shells under mechanical/thermal loading are as follows. Noor et al. ¹ studied the effects of both mechanical loadings and high temperature have to be considered in the design process of such structures. Studies involving the thermoelastic behavior using classical or first-order theories are described by Kant et al. ² Khdeir et al. ³ described the thermoelastic behavior using classical or first-order theories. Barut et al. ⁴ analyze the nonlinear thermoelastic behavior of shells by means of the Finite Element Method, but the assigned temperature profile is linear. Miller et al. ⁵ investigated thermal stress analysis of layered cylindrical shells. Dumir et al. ⁶ remarked the importance of the zigzag form of displacements in the thermal analysis of composite shells. Lee ⁷ studied the analysis of thermal stresses within multilayered cylinder under axial symmetry periodic boundary conditions. Ootao et al. ⁸ investigated an angle-ply laminated cylindrical panel with simply supported edges due to a nonuniform heat supply in the circumferential direction. Shahani et al. ⁹ solved transiented thermoelasticity problem in an isotropic thick-walled cylinder analytically using the finite Hankel transform. Talor and Radu ¹⁰ investigated the sinusoidal transient temperature for an isotropic hollow cylinder in an axis symmetric condition for the analyses of the thermal fatigue in the pipelines. Hocine ¹¹ preformed on thermomechanical behavior of multilayer tubular composite in axis metric steady state conditions. His results express that thermal effect has the greatest influence on radial stresses and strains. Zamani Nejad et al. ¹² presented a complete and consistent 3-D set of field equations by tensor analysis to characterize the behavior of FGM thick shells of revolution with arbitrary curvature and variable thickness along the meridional direction. Bakaiyan et al. ¹³ investigated the stress and deformation of the filament-wound pipes under combined internal pressure and temperature variations with axisymmetric and steady-state consideration. Kumar et al. ¹⁴ studied an extensive thermostructural analysis of composite structures, including temperature-dependent properties, mechanical, thermal, and thermochemical loads, and presented a coupled thermostructural analysis for a solid rocket motor nozzle with different orthotropic and isotropic materials. In the study by Yoo et al., ¹⁵ a kick motor nozzle incorporating spatially made from reinforced composites was subjected to thermoelastic analysis following the development of a material model used to homogenize various 3-D reinforcement architectures. Wang et al. ¹⁶ investigated the thermoelastic dynamic solution of a multilayered orthotropic hollow cylinder in the state of axisymmetric plane strain. Atefi et al. ¹⁷ presented an analytical solution for obtaining thermal stresses in a pipe caused by periodic time varying of temperature of medium fluid. Jabbari et al. ¹⁸ studied a general analysis of one-dimensional steady-state thermal stresses in a hollow thick cylinder made of functionally graded material. Shao et al. ¹⁹ presented thermomechanical analysis of functionally graded hollow cylinder subjected to axisymmetric mechanical and transient thermal loads. Kant ^20,21 investigated a general theory for small deformations of a thick shell made up of a layered system of different orthotropic materials having planes of symmetry coincident with the orthogonal reference frame and subjected to mechanical and arbitrary temperature distribution. Khare et al. ²² presented a closed-form solution for the thermomechanical analysis of laminated orthotropic shells. Khdeir et al. ^23,24 studied the thermoelastic governing equations for laminated shells that are exactly solved and then a higher shear deformation theory is given. They assume a linear or constant temperature profile through the thickness. Birsan ²⁵ presented an interesting method to analyze the thermal stresses in cylindrical Cosserat elastic shells for two given temperature fields. Iesan ²⁶ with assuming of polynomial temperature variation in the axial coordinate investigated thermal stresses in heterogeneous anisotropic Cosserat elastic cylinders. Barut et al. ²⁷ studied nonlinear thermoelastic analysis of composite panels under nonuniform temperature distribution with linear temperature profile. Miller et al. ²⁸ investigated thermal stress analysis of layered cylindrical shells. Dumir et al. ²⁹ is noteworthy, they presented the importance of the zigzag form of displacements in the thermal analysis of composite shells. In the case of shells, further investigations were developed by Hsu et al. ³⁰ for both Finite Element method and closed form and by Ding ³¹ for a weak formulation study for thick open laminated shell including the boundary conditions. The importance of thermoelasticity of cylindrically anisotropic generally laminated cylinders in such investigations was remarked by Padovan. ³² Some interesting experimental results in deformation analysis of thermally loaded composite tubes can be found in Holstein et al. ³³ Carrera et al. ³⁴ presented some results on thermal stress of layered plates and shells by using unified formulation. Cinefra et al. ³⁵ investigated thermal stress analysis of laminated structures by a variable kinematic MITC9 shell element. Cheila et al. ³⁶ investigated dynamic mechanical and thermal behavior analysis of composites based on polypropylene recycled with vegetal leaves. Cinefra et al. ³⁷ presented heat conduction and thermal stress analysis of laminated composites by a variable kinematic MITC9 shell element. To the best of authors’ knowledge, there is no published paper on the thermoelastic analysis of laminated cylindrical shells by using the mixed layerwise and first-order shear deformation theories (FSDTs) with differential quadrature method (DQM). In this study, the thermoelastic behavior of laminated cylindrical shells subjected to internal dynamic pressure and thermal loading achieved from experimental test are investigated by DQM. In order to more accurately estimate the thickness effect, the layerwise theory is used to approximate the temperature field while FSDTs is used to estimate displacement field. DQM is implemented to discretize the resulted governing equations while Newmark’s method is applied to solve the problem in the time domain. The MATLAB programming code have been used to solve differential quadrature equations. The accuracy and convergence of the algorithm are proved via different references.

Thermal governing equations

Heat transfer formulation in composite cylindrical vessel

In this section, the basic concepts of heat transfer formulation in composite cylindrical vessel are presented briefly. Figure 1 shows the geometry of such composite cylinder.

Figure 1.

Direction of fibers in a cylindrical laminate.

In this figure, r, ∅, and z are elements of off-axis coordinate system (reference coordinate system). If L is tangent line on cylinder in direction of fibers and t is tangent line on cylinder in direction of ∅, hence; angle of fibers (θ) is angle between L and t. Axisymmetric unsteady heat transfer in longitudinal and radial directions (r, z) has been considered in this study. The fibers are wound around the cylinder and the direction of the fibers in each lamina can vary between layers. Fourier law for cylindrical multilayer composite laminate as shown schematically in Figure 1 system is given by ³⁸

$\left\{\begin{array}{c}{q}_r\\ q_{\varnothing }\\ q_z\end{array}\right\}=-\left[\begin{array}{c}{\overline{k}}_{11} {\overline{k}}_{12} {\overline{k}}_{13}\\ {\overline{k}}_{21} {\overline{k}}_{22} {\overline{k}}_{23}\\ {\overline{k}}_{31} {\overline{k}}_{32} {\overline{k}}_{33}\end{array}\right]\left\{\begin{array}{c}\frac{\partial T}{\partial r}\\ \frac{1}{r}\frac{\partial T}{\partial \varnothing }\\ \frac{\partial T}{\partial z}\end{array}\right\}$

1

where q is heat flux, ${\overline{k}}_{ij}$ is conductive heat transfer coefficient and T is temperature. In order to study the heat conduction in composite laminates, two different coordinate systems should be defined. ^39,40 Here, (x₁, x₂, x₃) and (r, ∅, z) are considered as on-axis and off-axis coordinate systems, respectively. Off-axis coordinate system is defined to study the thermal properties in unique directions. The direction of on-axis coordinates depends on fiber orientation, in a way that x₁ is indirection of the fibers, x₂ is perpendicular to the fiber direction located in the plane of lamina, and x₃ is the third orthogonal direction. The fiber direction in each layer can be different from the other layers. Therefore, on-axis and off-axis systems have an angular deviation equal to θ in each lamina. In the off-axis coordinate system, by using the balance of energy for a cylindrical element, conductive heat transfer equation is obtained as follows relation

${\overline{k}}_{11}\frac{1}{r}\frac{\partial }{\partial r}\left(\frac{r\partial T}{\partial r}\right)+{\overline{k}}_{22}\frac{1}{r^2}\frac{{\partial }^{2}T}{\partial {\varnothing }^2}+{\overline{k}}_{33}\frac{{\partial }^{2}T}{\partial z^2}+\left({\overline{k}}_{12}+{\overline{k}}_{21}\right)\frac{1}{r}\frac{{\partial }^{2}T}{\partial \varnothing \partial r}+\left({\overline{k}}_{13}+{\overline{k}}_{31}\right)\frac{{\partial }^{2}T}{\partial r\partial z}+\frac{{\overline{k}}_{13}}{r}\frac{\partial T}{\partial z}+\left({\overline{k}}_{23}+{\overline{k}}_{32}\right)\frac{1}{r}\frac{{\partial }^{2}T}{\partial \varnothing \partial z}=\rho c_p\frac{\partial T}{\partial t}$

2

where ρ and c_p are density and specific heat capacity at constant pressure respectively and conductive heat coefficient ( $\overline{k}$ ) are defined as below ⁴¹

$\begin{array}{l}{\overline{k}}_{11}=k_{22}\\ {\overline{k}}_{22}=m_1^2{k}_{11}+n_1^2{k}_{22}\\ {\overline{k}}_{33}=n_1^2{k}_{11}+m_1^2{k}_{22}\\ {\overline{k}}_{12}={\overline{k}}_{21}=0\\ {\overline{k}}_{13}={\overline{k}}_{31}=0\\ {\overline{k}}_{23}={\overline{k}}_{32}=m_1{n}_1\left(k_{11}-k_{22}\right)\end{array}$

3

In these relations, k₁₁ is the conduction coefficient in parallel direction of fibers, and k₂₂ is the conduction coefficient in perpendicular direction of fibers, m₁ and n₁ represent cos θ and sin θ, respectively. With substituting equation (3) in equation (2), the heat transfer equation in cylindrical composite laminate is obtained

$k_{22}\frac{1}{r}\frac{\partial }{\partial r}\left(\frac{r\partial T}{\partial r}\right)+\left(m_1^2{k}_{11}+n_1^2{k}_{22}\right)\frac{1}{r^2}\frac{{\partial }^{2}T}{\partial {\varnothing }^2}+\left(n_1^2{k}_{11}+m_1^2{k}_{22}\right)\frac{{\partial }^{2}T}{\partial z^2}+2 m_1{n}_1\left(k_{11}-k_{22}\right)\frac{1}{r}\frac{{\partial }^2}{\partial \varnothing \partial z}=\rho c_p\frac{\partial T}{\partial t}$

4

With the assumption of axisymmetric unsteady heat transfer in longitudinal and radial directions (r, z), equation (4) can be written in a different way according to equation (5)

$\frac{{\partial }^{2}T}{\partial r^2}+\frac{1}{r}\frac{\partial T}{\partial r}+\frac{1}{{\mu }^2}\frac{{\partial }^{2}T}{\partial z^2}={\beta }^2\frac{\partial T}{\partial t}$

5

where

$\mu =\sqrt{\frac{k_{22}}{k_{11}{\text{sin}}^2\theta +k_{22}{\text{cos}}^2\theta }}$

6

$\beta =\sqrt{\frac{\rho c_p}{k_{22}}}$

7

Thermal governing equations in layerwise theory

Figure 2 shows schematic temperature distributions in linear layerwise theory. The temperature field in the form of layerwise theory for axisymmetric cylindrical shells is as follows ^42,43

$T\left(r,z,t\right)=\sum _{I=1}^{N_l}{T}_I\left(z,t\right){\phi }^I\left(r\right)N_l=N_e+1$

8

where N_e is the total number of mathematical layers, T_I is the values of the temperature at the bottom of the Ith mathematical layer in the r- and z-direction, respectively, and φ^I is interpolation approximation functions or a continuous function in terms of the thickness coordinate r. N^l represents the number of nodes in the r-direction of the shell, which depends on the number of considered mathematical layers N_e through the thickness.

Figure 2.

The through-thickness approximation of the temperature field.

By substituting equation (8) into equation (5) will be resulted in

$\sum _{I=1}^{N_l}{T}_I\frac{{\partial }^2\left({\phi }^I\right)}{\partial r^2}{k}_{22}+\sum _{I=1}^{N_l}\frac{T_I{k}_{22}}{r}\frac{\partial \left({\phi }^I\right)}{\partial r}+\sum _{I=1}^{N_l}\left(k_{22}{\text{cos}}^2\theta +k_{11}{\text{sin}}^2\theta \right){\phi }^I\frac{{\partial }^2\left(T_I\right)}{\partial z^2}=\sum _{I=1}^{N_l}{\rho }_0{c}_p\frac{\partial \left(T_I{\phi }^I\right)}{\partial t}$

9

The global linear interpolation function in the ith mathematical layer is defined as below ⁴⁴

$\begin{array}{l}{\phi }^1\left(r\right)=1-\frac{r-r_1}{r_2-r_1} r_1\le r\le r_2\\ {\phi }^I\left(r\right)=\{\begin{array}{c}\frac{r-r_{I-1}}{r_I-r_{I-1}} r_{I-1}\le r\le r_I\\ 1-\frac{r-r_I}{r_{I+1}-r_I} r_I\le r\le r_{I+1}\end{array}, I=2,3,\dots , N_e\\ {\phi }^{N_l}\left(r\right)=\frac{r-r_{Ne}}{r_{Ne+1}-r_{Ne}}, r_{Ne}\le r\le r_{Ne+1}\end{array}$

10

Form the first derivative of interpolation functions are as follows

$\begin{array}{l}\frac{\partial {\phi }^1}{\partial r}=-\frac{1}{r_2-r_1} r_1\le r\le r_2\\ \frac{\partial {\phi }^I}{\partial r}=\{\begin{array}{c}\frac{1}{r_I-r_{I-1}}, r_{I-1}\le r\le r_I\\ \frac{-1}{r_{I+1}-r_I}, r_{Ne}\le r\le r_{Ne+1}\end{array} I=2,3,\dots , N_e\\ \frac{\partial {\phi }^{N_l}}{\partial r}=\frac{1}{r_{Ne+1}-r_{Ne}}={\Gamma }^I r_{Ne}<r<r_{Ne+1}\end{array}$

11

Second-order derivatives interpolation functions will be equal to zero. Therefore, to compensate, this problem has been assisting from finite difference method

$\begin{array}{l}\frac{{\partial }^{2}T}{\partial r^2}=\sum _{I=1}^{N_l}\left(\frac{T_{I+1}+T_{I-1}-2T_I}{{\left(r_{I+1}-r_I\right)}^2}\right){\phi }^I:\text{ Inner layers}\\ \frac{{\partial }^{2}T}{\partial r^2}=\sum _{I=1}^{N_l}\left(\frac{2T_{Ne}-2T_{Ne+1}}{{\left(r_{Ne+1}-r_{Ne}\right)}^2}\right){\phi }^I:\text{ Top layer}\\ \frac{{\partial }^{2}T}{\partial r^2}=\sum _{I=1}^{N_l}\left(\frac{2T_2-2T_1}{{\left(r_2-r_1\right)}^2}\right){\phi }^I:\text{ Bottom layer}\end{array}$

12

So, with the help of derivatives interpolation function (equations (11) to (12)) and equation (9) will be resulted in

$k_{22}\frac{{\partial }^{2}T}{\partial r^2}+\sum _{I=1}^{N_l}\frac{T_I{k}_{22}}{r}{\Gamma }^I+\sum _{I=1}^{N_l}\left(k_{22}{\text{cos}}^2\theta +k_{11}{\text{sin}}^2\theta \right){\phi }^I\frac{{\partial }^2\left(T_I\right)}{\partial z^2}=\sum _{I=1}^{N_l}{\rho }_0{c}_p{\phi }^I\frac{\partial \left(T_I\right)}{\partial t}$

13

where Γ ^I is calculated from equation (11). The governing equations for the present layerwise theory can be derived using the variational principle ^{44 –48}

$\underset{0}{\overset{t}{\int }}\left(\delta E\right)\text{d}t=0$

14

where $\delta E=E\cdot \delta T$ is variational equation for heat transfer given by

$\delta E=\underset{\Omega }{\int }\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}\left(k_{22}\frac{{\partial }^{2}T}{\partial r^2}+\sum _{I=1}^{N_l}\frac{T_I{k}_{22}}{r}{\Gamma }^I+\sum _{I=1}^{N_l}\left(k_{22}{\text{cos}}^2\theta +k_{11}{\text{sin}}^2\theta \right){\phi }^I\frac{{\partial }^2\left(T_I\right)}{\partial z^2}-\sum _{I=1}^{N_l}{\rho }_0{c}_p{\phi }^I\frac{\partial \left(T_I\right)}{\partial t}\right).\left(\delta T\right)\text{d}r\text{d}A$

15

where h is thick laminated cylindrical shell. By substituting equations (12) and (8) into equation (15), obtain

$\delta E=\underset{\Omega }{\int }\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}\left(\sum _{I=1}^{N_l}{k}_{22}\left(\frac{T_{I+1}+T_{I-1}-2T_I}{{\left(r_{I+1}-r_I\right)}^2}\right){\phi }^I+\sum _{I=1}^{N_l}\frac{T_I}{r}{\Gamma }^I{k}_{22}+\sum _{I=1}^{N_l}\left(k_{22}{\text{cos}}^2\theta +k_{11}{\text{sin}}^2\theta \right){\phi }^I\frac{{\partial }^2\left(T_I\right)}{\partial z^2}-\sum _{I=1}^{N_l}{\rho }_0{c}_p{\phi }^I\frac{\partial \left(T_I\right)}{\partial t}\right)\left(\sum _{J=1}^{N_l}\delta T_J{\phi }^J\right)\text{d}r\text{d}A=\underset{\Omega }{\int }\left(\sum _{I=1}^{N_l}\sum _{J=1}^{N_l}\left(\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}{k}_{22}\left(\frac{T_{I+1}+T_{I-1}-2T_I}{{\left(r_{I+1}-r_I\right)}^2}\right){\phi }^I{\phi }^J+T_I\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}\frac{{\Gamma }^I{\phi }^J{k}_{22}}{r}+\frac{{\partial }^2\left(T_I\right)}{\partial z^2}\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}\left(k_{22}{\text{cos}}^2\theta +k_{11}{\text{sin}}^2\theta \right){\phi }^I{\phi }^J-\frac{\partial \left(T_I\right)}{\partial t}\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}{\rho }_0{C}_E{\phi }^I{\phi }^J\right)\delta T_J\right)\text{d}r\text{d}A$

16

Simplification of the above equation in equations (17) to (20) is expressed

$\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}\left(\frac{k_2}{{\rho }_0{c}_p}\right){\phi }^I{\phi }^J\text{d}r=\sum _{k=1}^{N_e}\underset{r_k}{\overset{r_{k+1}}{\int }}\left(\frac{k_2}{{\rho }_0{c}_p}\right){\phi }^I{\phi }^J\text{d}r=E^{IJ}$

17

$\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}\frac{{\Gamma }^I{\phi }^J\left(\frac{k_2}{{\rho }_0{c}_p}\right)}{r}\text{d}r=\sum _{k=1}^{N_e}\underset{r_k}{\overset{r_{k+1}}{\int }}\frac{{\Gamma }^I{\phi }^J\left(\frac{k_2}{{\rho }_0{c}_p}\right)}{r}\text{ d}r=B^{IJ}$

18

$\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}\frac{\left(k_{22}{\text{cos}}^2\theta +k_{11}{\text{sin}}^2\theta \right)}{{\rho }_0{c}_p}{\phi }^I{\phi }^J\text{d}r=\sum _{k=1}^{N_e}\underset{r_k}{\overset{r_{k+1}}{\int }}\frac{\left(k_{22}{\text{cos}}^2\theta +k_{11}{\text{sin}}^2\theta \right)}{{\rho }_0{c}_p}{\phi }^I{\phi }^J\text{ d}r=C^{IJ}$

19

$\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}{\phi }^I{\phi }^J\text{d}r=\sum _{k=1}^{N_e}\underset{r_k}{\overset{r_{k+1}}{\int }}{\phi }^I{\phi }^J \text{d}r=D^{IJ}$

20

where N_e is number of layers of envisaged in layerwise theory. By substituting equations (17) to (20) into equation (16), obtain

$\delta E=\underset{\Omega }{\int }\left(\sum _{I=1}^{N_l}\sum _{J=1}^{N_l}\left(\left(\frac{T_{I+1}+T_{I-1}-2T_I}{{\left(r_{I+1}-r_I\right)}^2}\right)E^{IJ}+T_I{B}^{IJ}+\frac{{\partial }^2\left(T_I\right)}{\partial z^2}{C}^{IJ}-\frac{\partial \left(T_I\right)}{\partial t}{D}^{IJ}\right)\delta T_J\right)\text{d}A$

21

By substituting equation (21) into equation (14) and by using the Euler-Lagrange equations are derived temperature governing equations in laminated cylindrical shells with layerwise theory as follows

$\{\begin{array}{l}\sum _{I=1}^{N_l}\left(\left(\frac{T_{I+1}+T_{I-1}-2T_I}{{\left(r_{I+1}-r_I\right)}^2}\right)E^{I1}+T_I{B}^{I1}+\frac{{\partial }^2\left(T_I\right)}{\partial z^2}{C}^{I1}\right)=\sum _{I=1}^{N_l}\left(\frac{\partial \left(T_I\right)}{\partial t}{D}^{I1}\right)\\ \sum _{I=1}^{N_l}\left(\left(\frac{T_{I+1}+T_{I-1}-2T_I}{{\left(r_{I+1}-r_I\right)}^2}\right)E^{I2}+T_I{B}^{I2}+\frac{{\partial }^2\left(T_I\right)}{\partial z^2}{C}^{I2}\right)=\sum _{I=1}^{N_l}\left(\frac{\partial \left(T_I\right)}{\partial t}{D}^{I2}\right)\\ ⋮\\ ⋮\\ ⋮\\ ⋮\\ \sum _{I=1}^{N_l}\left(\left(\frac{T_{I+1}+T_{I-1}-2T_I}{{\left(r_{I+1}-r_I\right)}^2}\right)E^{I{N}_l}+T_I{B}^{I{N}_l}+\frac{{\partial }^2\left(T_I\right)}{\partial z^2}{C}^{I{N}_l}\right)=\sum _{I=1}^{N_l}\left(\frac{\partial \left(T_I\right)}{\partial t}{D}^{I{N}_l}\right)\end{array}$

22

where N_l is the number of nodes in layerwise theory. For the integration of the second derivative term in the above equations due to the equation (12) can be used the following equation

$\sum _{I=1}^{N_l}\sum _{J=1}^{N_l}\left(\left(\frac{T_{I+1}+T_{I-1}-2T_I}{{\left(r_{I+1}-r_I\right)}^2}\right)E^{IJ}\right)=\sum _{J=1}^{N_l}\left(\sum _{I=2}^{N_e}\left(\left(\frac{T_{I+1}+T_{I-1}-2T_I}{{\left(r_{I+1}-r_I\right)}^2}\right)E^{IJ}\right)+\left(\left(\frac{2T_2-2T_1}{{\left(r_2-r_1\right)}^2}\right)E^{IJ}\right)+\left(\frac{2T_{Ne}-2T_{N_l}}{{\left(r_{N_l}-r_{Ne}\right)}^2}\right)E^{IJ}\right)=\sum _{I=1}^{N_l}\sum _{J=1}^{N_l}\left(K_{\text{diff}}^{IJ}\right)\left\{\begin{array}{c}\left\{T_1\right\}\\ \left\{T_2\right\}\\ ⋮\\ \left\{T_{N_l}\right\}\end{array}\right\}$

23

where $\left\{T_1\right\}\dots .\left\{T_{N_l}\right\}$ are temperature in nodes.

DQM for thermal governing equations

At this stage, the governing differential equations are discretized into algebraic equations using DQM. In this method, the first- and second-order partial derivatives of a field variable at the ith discrete point in the computational domain are approximated by a weighted linear summation of the field variable at all discrete nodes located along the line that passes through that point in the domain. A differential quadrature approximation at the ith discrete point on a grid at the z-axis is given by ^49,50

$\frac{\partial T\left(z_i,r,t\right)}{\partial z}=T^{\left(1\right)}\left(z_i,r,t\right)=\sum _{j=1}^{N_z}{a}_{ij}T\left(z_i,r,t\right), \text{for} i=1,2,\dots ,N_z$

24

$a_{ij}=\frac{M^{\left(1\right)}\left(z_i\right)}{\left(z_i-z_k\right)\cdot M^{\left(1\right)}\left(z_j\right)},\text{ for}\hspace{0.17em}i\ne j$

25

$M^{\left(1\right)}\left(z_i\right)=\prod _{k=1,k\ne i}^{N_z}\left(z_i-z_k\right)$

26

$\sum _{j=1}^{N_z}{a}_{ij}=0or{a}_{ii}=-\sum _{j=1,j\ne i}^{N_z}{a}_{ij}$

27

where $T\left(z_i,r,t\right)$ represents the functional value at a grid point z_i, $T^{\left(1\right)}\left(z_i,r,t\right)$ indicates the first-order derivative of $T\left(z_i,r,t\right)$ at z_i, N_z is the number of grid points in the direction of z, and a_ij is the weighting coefficient of the first-order derivative. The generalized weighted coefficients in DQM are based on Lagrange interpolation method. Therefore, nondiagonal elements of the weight matrix and diagonal elements are obtained from equations (25) and (27), respectively. Second-order derivatives of the DQM can be obtained from the following equation.

$b_{ij}=2a_{ij}\left(a_{ii}-\frac{1}{z_i-z_j}\right),\text{ for}\hspace{0.17em}i\ne j$

28

$\sum _{j=1}^{N_z}{b}_{ij}=0or{b}_{ii}=-\sum _{j=1,j\ne i}^{N_z}{b}_{ij}$

29

Higher order derivatives of the DQM can be obtained as follows

$\begin{array}{l}{f}_{ij}{ }^{\left(m\right)}=m\left(f_{ii}{ }^{\left(m-1\right)}{f}_{ij}{ }^{\left(m-1\right)}-\frac{f_{ij}{ }^{\left(m-1\right)}}{z_i-z_j}\right)\\ for i\ne j \& m=2,\dots .N_z-1i, j=1,2\dots N_z\end{array}$

30

$\begin{array}{l}{f}_{ij}{ }^{\left(m\right)}=-\sum _{j=1, j\ne 1}^{N_z}{f}_{ij}{ }^{\left(m\right)}\\ for i=1,2,\dots ,N_z \& m=1,2,\dots .N_z-1\end{array}$

31

For variables r and t equations (24) to (31) exist that is omitted for brevity.

Newmark’s time integration method for thermal governing equations

As an alternative method and in order to obtain a comparative solution, Newmark’s time integration scheme is applied in the time domain to the DQM discretized form of the governing equations and the related boundaries. In this procedure, the dynamical system of equation (22) at time (t + Δt) is solved according to the equations (32) and (33) ⁴⁴

$\begin{array}{l}{\ddot{T}}_{t+\Delta t}=a_0\left(T_{t+\Delta t}-T_t\right)-a_2\dot{T_t}-a_3{\ddot{T}}_t\\ {\dot{T}}_{t+\Delta t}={\dot{T}}_t+a_6{\ddot{T}}_t+a_7{\ddot{T}}_t\end{array}$

32

$\begin{array}{l}{a}_0=\frac{1}{\alpha {\left(\Delta t\right)}^2}; a_1=\frac{\delta }{\alpha \Delta t};a_2=\frac{1}{\alpha \Delta t}; a_3=\frac{1}{2\alpha }-1\\ a_4=\frac{\delta }{a}-1; a_5=\frac{\Delta t}{2}\left(\frac{\delta }{\alpha }-2\right); a_6=\Delta t\left(1-\delta \right); a_7=\delta \Delta t;\hspace{0.17em}\Delta t=\frac{t_{\text{final}}}{N_t}\end{array}$

33

where Δt is time step, $\dot{T}$ is the first-order derivative of T at t, $\ddot{T}$ is the second-order derivative of T at t, N_t is the number of time steps, and α and δ are time integration parameters which are set equal to $\alpha =0.25{\left(1 + \gamma \right)}^2$ and δ = 0.5 + γ with γ = 0.005. For solution of problem, derived equations to layerwise theory (equation (22)) by using equations (24) to (33) can be written to solvable equations with the DQM as follows

$\begin{array}{l}\frac{\partial T_I}{\partial z}=\sum _{k=1}^{N_Z}{a}_{i,k}{T}_I{\hspace{0.17em}}_{k,j}=\left[a_{ij}\right]\left[T_I\right]\\ \frac{{\partial }^2{T}_I}{\partial z}=\sum _{k=1}^{N_Z}{b}_{i,k}{T}_I{\hspace{0.17em}}_{k,j}=\left[b_{ij}\right]\left[T_I\right]\\ \frac{\partial T_I}{\partial t}=\left[{\delta }_{\left(M_{\varnothing }{N}_z,M_{\varnothing }{N}_z\right)}\right]\left[{\dot{T}}_I\right]\end{array}$

34

where [T_I] is temperature in nodes; $\left[\dot{T}\right]$ is the first-order derivative of [T] at t; N_Z is the number of grid points in the direction of z ( $N_Z=N_z+1$ ); $\left[{\delta }_{\left(M_{\varnothing }{N}_z,M_{\varnothing }{N}_z\right)}\right]$ is the kronecker delta matrix(the unit or identity matrix); M_∅ and N_z the number of elements in the ∅ and z-direction; and [a_ij] and [b_ij] are the weighting coefficient of the first-order derivative and second-order derivatives, respectively. After substituting equations (34) into the temperature governing equations (equations (22)), result simplified in the following form

$\left[\begin{array}{cccc}{K}_{11}& K_{12}& \dots & K_{1N_l}\\ K_{21}& K_{22}& \dots & K_{2N_l}\\ ⋮& ⋮& ⋮& ⋮\\ K_{N_{l}1}& K_{N_{l}2}& \dots & K_{N_l{N}_l}\end{array}\right]\left\{\begin{array}{c}\left\{T_1\right\}\\ \left\{T_2\right\}\\ ⋮\\ \left\{T_{N_l}\right\}\end{array}\right\}-\left[\begin{array}{cccc}{\Lambda }_{11}& {\Lambda }_{12}& \dots & {\Lambda }_{1N_l}\\ {\Lambda }_{21}& {\Lambda }_{22}& \dots & {\Lambda }_{2N_l}\\ ⋮& ⋮& ⋮& ⋮\\ {\Lambda }_{N_{l}1}& {\Lambda }_{N_{l}2}& \dots & {\Lambda }_{N_l{N}_l}\end{array}\right]\left\{\begin{array}{c}\left\{{\dot{T}}_1\right\}\\ \left\{{\dot{T}}_2\right\}\\ ⋮\\ \left\{{\dot{T}}_{N_l}\right\}\end{array}\right\}=\left\{\begin{array}{c}\left\{0\right\}\\ \left\{0\right\}\\ ⋮\\ \left\{0\right\}\end{array}\right\}$

35

where B^JI, C^JI, and D^JI are obtained from equations (17) to (20); $K_{\text{diff}}^{IJ}$ is obtained from equation (23); $T_1\dots T_{N_l}$ and ${\dot{T}}_1\dots$ ${\dot{T}}_{N_l}$ are temperature and temperature changes versus time in nodes; N_l is the number of nodes in the r-direction of the vessel; I and J are the counter number of nodes integral equations in layerwise theory; and K_IJ and Λ _IJ are expressed by the following equations

$\begin{array}{l}{K}_{IJ}=\left(B^{JI}\right)\left[{\delta }_{\left(M_{\varnothing }{N}_z,M_{\varnothing }{N}_z\right)}\right]+C^{JI}\left[b_{ij}\right]+K_{diff}^{IJ}\left[{\delta }_{\left(M_{\varnothing }{N}_z,M_{\varnothing }{N}_z\right)}\right]J,I=1\dots N_l\\ {\Lambda }_{IJ}=D^{JI}\left[{\delta }_{\left(M_{\varnothing }{N}_z,M_{\varnothing }{N}_z\right)}\right]J,I=1\dots N_l\end{array}$

36

Basic equations in fiber-reinforced composites under thermoelastic load

In the present study, FSDT is employed to simulate the deformation of every layer of the cylinder in terms of deformation of mid-surface and rotation about outward axis of the mid-surface. The temperature field is calculated from layerwise theory and to be put in thermoelastic equations of FSDT. A circular cylindrical shell with radius of R, thickness h, and the length l is shown in Figure 3. The u₁, u₂, and u₃ are the displacement components in the axial, tangential, and radial directions, respectively, and that the deformations are assumed to be small. The displacement field of the first-order theory is of the form:

$\begin{array}{l}{u}_1\left(x_1,x_2,x_3,t\right)=u_0\left(x_1,x_2,t\right)+x_3{\phi }_1\left(x_1,x_2,t\right)\\ u_2\left(x_1,x_2,x_3,t\right)=v_0\left(x_1,x_2,t\right)+x_3{\phi }_2\left(x_1,x_2,t\right)\\ u_3\left(x_1,x_2,x_3,t\right)=w_0\left(x_1,x_2,t\right)\end{array}$

37

where ( $x_1,x_2,x_3$ ) are axis of the cylinder (see Figure 3), ( $u_0,v_0,w_0$ ) denote the displacements of a point on the plane x₃ = 0 (mid-plane), ( ${\phi }_1,{\phi }_2$ ) denote the rotation of a transverse normal about the x₁-axis and the x₂-axis, respectively, and t is an arbitrary time. Based on the FSDT, the equilibrium equations for a shell under axial loads are as follows ^44,51

$\begin{array}{l}\frac{\partial N_1}{\partial x_1}+\frac{\partial N_6}{\partial x_2}=I_0{\ddot{u}}_0+I_1{\ddot{\phi }}_1\\ \frac{\partial N_6}{\partial x_1}+\frac{\partial N_2}{\partial x_2}+\frac{Q_2}{R}=I_0{\ddot{v}}_0+I_1{\ddot{\phi }}_2\\ \frac{\partial Q_1}{\partial x_1}+\frac{\partial Q_2}{\partial x_2}-\frac{N_2}{R}+q-Ň\frac{{\partial }^2{w}_0}{\partial x_1^2}=I_0{\ddot{w}}_0\\ \frac{\partial M_1}{\partial x_1}+\frac{\partial M_6}{\partial x_2}-Q_1=I_1{\ddot{u}}_0+I_2{\ddot{\phi }}_1\\ \frac{\partial M_6}{\partial x_1}+\frac{\partial M_2}{\partial x_2}-Q_2=I_1{\ddot{v}}_0+I_2{\ddot{\phi }}_2\end{array}$

38

where Ň is the axial compressive load (positive in compression); I₀, I₁, and I₂ are the mass inertias and are defined in equation (39); $q=q_b+q_t$ , q_b, and q_t are internal and external dynamic pressures, respectively; and ( $N_1=N_{11}, N_2=N_{22}, N_6=N_{12})$ and ( $M_1=M_{11}, M_2=M_{22}, M_6=M_{12})$ are the in-plane force resultants and the moment resultants, respectively. Constitutive equations of composite shell based on classical laminate theory are defined by the following relations ⁴⁴

$\left(I_0,I_1,I_2\right)=\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}\left(1,x_3,x_3^2\right){\rho }_{k}d{x}_3$

39

$\begin{array}{l}\left[\begin{array}{c}{N}_1\\ N_2\\ N_6\end{array}\right]=\left[\begin{array}{ccc}{A}_{11}& A_{12}& A_{16}\\ A_{21}& A_{22}& A_{26}\\ A_{61}& A_{62}& A_{66}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_1{\hspace{0.17em}}^{\left(0\right)}\\ {\epsilon }_2{\hspace{0.17em}}^{\left(0\right)}\\ {\epsilon }_6{\hspace{0.17em}}^{\left(0\right)}\end{array}\right\}+\left[\begin{array}{ccc}{B}_{11}& B_{12}& B_{16}\\ B_{21}& B_{22}& B_{26}\\ B_{61}& B_{62}& B_{66}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_1{\hspace{0.17em}}^{\left(1\right)}\\ {\epsilon }_2{\hspace{0.17em}}^{\left(1\right)}\\ {\epsilon }_6{\hspace{0.17em}}^{\left(1\right)}\end{array}\right\}-\left[\begin{array}{c}{N}_1^T\\ N_2^T\\ N_6^T\end{array}\right]\\ \left[\begin{array}{c}{M}_1\\ M_2\\ M_6\end{array}\right]=\left[\begin{array}{ccc}{B}_{11}& B_{12}& B_{16}\\ B_{21}& B_{22}& B_{26}\\ B_{61}& B_{62}& B_{66}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_1{\hspace{0.17em}}^{\left(0\right)}\\ {\epsilon }_2{\hspace{0.17em}}^{\left(0\right)}\\ {\epsilon }_6{\hspace{0.17em}}^{\left(0\right)}\end{array}\right\}+\left[\begin{array}{ccc}{D}_{11}& D_{12}& D_{16}\\ D_{21}& D_{22}& D_{26}\\ D_{61}& D_{62}& D_{66}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_1{\hspace{0.17em}}^{\left(1\right)}\\ {\epsilon }_2{\hspace{0.17em}}^{\left(1\right)}\\ {\epsilon }_6{\hspace{0.17em}}^{\left(1\right)}\end{array}\right\}-\left[\begin{array}{c}{M}_1^T\\ M_2^T\\ M_6^T\end{array}\right]\\ \left[\begin{array}{c}{Q}_1\\ Q_2\end{array}\right]=K_S\left[\begin{array}{cc}{A}_{55}& A_{45}\\ A_{45}& A_{44}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_4{\hspace{0.17em}}^{\left(0\right)}\\ {\epsilon }_5{\hspace{0.17em}}^{\left(0\right)}\end{array}\right\}\end{array}$

40

where ρ_k is the density for each layer; K_S is the shear correction factor introduced by Mindlin and is equal to $\frac{{\pi }^2}{12}$ ⁵¹ ; ( $N_1^T$ , ${\text{N}}_2^{\text{T}}$ , $N_6^T$ , $M_1^T$ , $M_2^T$ , $M_6^T$ ) denote thermal resultants are defined in equation (41); A_ij, B_ij and D_ij are the extensional, coupling, and bending stiffness matrices, respectively, and they are defined in equation (42) ⁴⁴ ; and the strains ( ${\epsilon }_1^0={\epsilon }_{11}^0$ , ${\epsilon }_2^0={\epsilon }_{22}^0$ , ${\epsilon }_6^0={\gamma }_{12}^0,{\epsilon }_4^0={\gamma }_{23}^0,{\epsilon }_5^0={\gamma }_{13}^0$ ) are the mid-surface strains and ( ${\epsilon }_1^1={\epsilon }_{11}^1$ , ${\epsilon }_2^1={\epsilon }_{22}^1$ , ${\epsilon }_6^1={\gamma }_{12}^1$ ) are linear strains through the laminate thickness, which are defined in equation (43) ⁴⁴

$\begin{array}{l}\left\{N^T\right\}=\sum _{k=1}^n\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}\left({\left[\overline{Q}\right]}^{\left(k\right)}{\left\{\overline{\alpha }\right\}}^{\left(k\right)}\Delta T\right)\text{d}{x}_3\\ \left\{M^T\right\}=\sum _{k=1}^n\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}\left({\left[\overline{Q}\right]}^{\left(k\right)}{\left\{\overline{\alpha }\right\}}^{\left(k\right)}\Delta T\right)x{}_3\text{d}{x}_3\end{array}$

41

where ΔT = T − T₀; T₀ is the reference temperature; ${\left[\overline{Q}\right]}^{\left(k\right)}={\overline{Q}}_{ij}{\hspace{0.17em}}^{\left(k\right)}$ is the transformed elastic stiffness matrix ⁴⁴ ; n is the number of lamina; and ${\left\{\overline{\alpha }\right\}}^{\left(k\right)}={\overline{\alpha }}_{ij}^{\left(k\right)}$ is the transformed thermal coefficients of expansion matrix ⁴⁴

$(A_{ij},B_{ij},D_{ij})=\underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}\left(1,x_3,x_3^2\right){\overline{\text{Q}}}_{\text{ij}}{\hspace{0.17em}}^{\left(\text{k}\right)}d{x}_3\left(i,j=1\dots 6\right)$

42

$\begin{array}{l}\left\{\begin{array}{c}{\epsilon }_1{\hspace{0.17em}}^{\left(0\right)}\\ {\epsilon }_2{\hspace{0.17em}}^{\left(0\right)}\\ {\epsilon }_6{\hspace{0.17em}}^{\left(0\right)}\end{array}\right\}=\left\{\begin{array}{c}\frac{\partial u_0}{\partial x_1}+\frac{1}{2}{\left(\frac{\partial w_0}{\partial x_1}\right)}^2\\ \frac{\partial v_0}{\partial x_2}+\frac{1}{2}{\left(\frac{\partial w_0}{\partial x_2}\right)}^2\\ \frac{\partial u_0}{\partial x_2}+\frac{\partial v_0}{\partial x_1}+\frac{\partial w_0}{\partial x_1}\frac{\partial w_0}{\partial x_2}\end{array}\right\} \left\{\begin{array}{c}{\epsilon }_1{\hspace{0.17em}}^{\left(1\right)}\\ {\epsilon }_2{\hspace{0.17em}}^{\left(1\right)}\\ {\epsilon }_6{\hspace{0.17em}}^{\left(1\right)}\end{array}\right\}=\left\{\begin{array}{c}\frac{\partial {\varphi }_1}{\partial x_1}\\ \frac{\partial {\varphi }_2}{\partial x_2}\\ \frac{\partial {\varphi }_1}{\partial x_2}+\frac{\partial {\varphi }_2}{\partial x_1}\end{array}\right\}\\ \left\{\begin{array}{c}{\epsilon }_4{\hspace{0.17em}}^{\left(0\right)}\\ {\epsilon }_5{\hspace{0.17em}}^{\left(0\right)}\end{array}\right\}=\left\{\begin{array}{c}\frac{\partial w_0}{\partial x_2}+{\varphi }_2\\ \frac{\partial w_0}{\partial x_1}+{\varphi }_1\end{array}\right\}\end{array}$

43

Figure 3.

Laminated shell geometry in cylindrical coordinate system and stress resultants on a shell element.

At this stage, the governing differential equations are discretized into algebraic equations using DQM. Newmark’s time integration scheme is applied in the time domain to the DQM discretized form of the governing equations. DQM and Newmark’s method are presented in “Differential quadrature method for thermal governing equations” and “Newmark’s time integration method for thermal governing equations” sections for temperature field. In the following, this relationship is used for displacement field as follows

$\begin{array}{l}\frac{\partial u_0}{\partial x_1}=\left[a_{ij}\right]\left[u_0\right] \frac{\partial v_0}{\partial x_1}=\left[a_{ij}\right]\left[v_0\right] \frac{\partial w_0}{\partial x_1}=\left[a_{ij}\right]\left[w_0\right] \frac{\partial {\phi }_1}{\partial x_1}=\left[a_{ij}\right]\left[{\phi }_1\right] \frac{\partial {\phi }_2}{\partial x_1}=\left[a_{ij}\right]\left[{\phi }_2\right]\\ \frac{\partial u_0}{\partial x_2}=\left[c_{ij}\right]\left[u_0\right] \frac{\partial v_0}{\partial x_2}=\left[c_{ij}\right]\left[v_0\right] \frac{\partial w_0}{\partial x_2}=\left[c_{ij}\right]\left[w_0\right] \frac{\partial {\phi }_1}{\partial x_2}=\left[c_{ij}\right]\left[{\phi }_1\right] \frac{\partial {\phi }_2}{\partial x_2}=\left[c_{ij}\right]\left[{\phi }_2\right]\\ \frac{{\partial }^2{u}_0}{\partial x_1^2}=\left[b_{ij}\right]\left[u_0\right] \frac{{\partial }^2{v}_0}{\partial x_1^2}=\left[b_{ij}\right]\left[v_0\right] \frac{{\partial }^2{w}_0}{\partial x_1^2}=\left[b_{ij}\right]\left[w_0\right] \frac{{\partial }^2{\phi }_1}{\partial x_1^2}=\left[b_{ij}\right]\left[{\phi }_1\right] \frac{{\partial }^2{\phi }_2}{\partial x_1^2}=\left[b_{ij}\right]\left[{\phi }_2\right]\\ \frac{{\partial }^2{u}_0}{\partial x_2^2}=\left[d_{ij}\right]\left[u_0\right] \frac{{\partial }^2{v}_0}{\partial x_2^2}=\left[d_{ij}\right]\left[v_0\right] \frac{{\partial }^2{w}_0}{\partial x_2^2}=\left[d_{ij}\right]\left[w_0\right] \frac{{\partial }^2{\phi }_1}{\partial x_2^2}=\left[d_{ij}\right]\left[{\phi }_1\right] \frac{{\partial }^2{\phi }_2}{\partial x_2^2}=\left[d_{ij}\right]\left[{\phi }_2\right]\\ \frac{{\partial }^2{u}_0}{\partial x_1\partial x_2}=\left[e_{ij}\right]\left[c_{ij}\right]\left[u_0\right] \frac{{\partial }^2{v}_0}{\partial x_1\partial x_2}=\left[e_{ij}\right]\left[c_{ij}\right]\left[v_0\right] \frac{{\partial }^2{w}_0}{\partial x_1\partial x_2}=\left[e_{ij}\right]\left[c_{ij}\right]\left[w_0\right] \frac{{\partial }^2{\phi }_1}{\partial x_1\partial x_2}=\left[e_{ij}\right]\left[c_{ij}\right]\left[{\phi }_1\right] \frac{{\partial }^2{\phi }_2}{\partial x_1\partial x_2}=\left[e_{ij}\right]\left[c_{ij}\right]\left[{\phi }_2\right]\\ \frac{\partial T^{*}}{\partial x_1}=\left[a_{ij}\right]\left[T^{*}\right] \frac{\partial T^{*}}{\partial x_2}=\left[c_{ij}\right]\left[T^{*}\right] \frac{{\partial }^2{T}^{*}}{\partial x_1^2}=\left[b_{ij}\right]\left[T^{*}\right] \frac{{\partial }^2{T}^{*}}{\partial x_2^2}=\left[d_{ij}\right]\left[T^{*}\right]\end{array}$

44

where a_ij, b_ij, c_ij, d_ij, and e_ij are the weighting coefficient matrices in DQM method; $\left(u_0,v_0,w_0\right)$ are the displacements of a point ( $x_1,x_2,0)$ on the mid-surface of the shell and $\left({\phi }_1,{\phi }_2\right)$ are the rotations of a normal to the reference surface that are defined in equations (45); and $T^{*}=T-T_0$ , T₀ is the reference temperature and is defined in equation (46)

$\left\{u_0\right\}=\left\{\begin{array}{c}{u}_0{\hspace{0.17em}}_{11}\\ u_0{\hspace{0.17em}}_{12}\\ ⋮\\ u_0{\hspace{0.17em}}_{1M}\\ ⋮\\ u_0{\hspace{0.17em}}_{NM}\end{array}\right\},\left\{v_0\right\}=\left\{\begin{array}{c}{v}_0{\hspace{0.17em}}_{11}\\ v_0{\hspace{0.17em}}_{12}\\ ⋮\\ v_0{\hspace{0.17em}}_{1M}\\ ⋮\\ v_0{\hspace{0.17em}}_{NM}\end{array}\right\}\dots .\left\{{\phi }_2\right\}=\left\{\begin{array}{c}{\phi }_2{\hspace{0.17em}}_{11}\\ {\phi }_2{\hspace{0.17em}}_{12}\\ ⋮\\ {\phi }_2{\hspace{0.17em}}_{1M}\\ ⋮\\ {\phi }_2{\hspace{0.17em}}_{NM}\end{array}\right\}$

45

$\left\{T^{*}\left(x_1,x_2,{\xi }_k,t\right)\right\}=\left\{\begin{array}{c}{T}^{*}{\left({\xi }_k,t\right)}_{11}\\ T^{*}{\left({\xi }_k,t\right)}_{12}\\ ⋮\\ T^{*}{\left({\xi }_k,t\right)}_{1M}\\ ⋮\\ T^{*}{\left({\xi }_k,t\right)}_{NM}\end{array}\right\}$

46

where M and N are the number of elements in the x₁- and x₂-direction, $\xi =x_3$ , the superscript k denotes the layer number, and t is an arbitrary time. After substituting equations (40) and (44) into the equations of motion (38) and with omitted the small strain (nonlinear terms in the strain), result simplified in the following form

$\left[\begin{array}{ccccc}{K}_{11}& K_{12}& K_{13}& K_{14}& K_{15}\\ K_{21}& K_{22}& K_{23}& K_{24}& K_{25}\\ K_{31}& K_{32}& K_{33}& K_{34}& K_{35}\\ K_{41}& K_{42}& K_{43}& K_{44}& K_{45}\\ K_{51}& K_{52}& K_{53}& K_{54}& K_{55}\end{array}\right]\left\{\begin{array}{c}\left\{u_0\right\}\\ \left\{v_0\right\}\\ \left\{w_0\right\}\\ \left\{{\phi }_1\right\}\\ \left\{{\phi }_2\right\}\end{array}\right\}+\left[\begin{array}{ccccc}\left[I_0\right]& 0& 0& \left[I_1\right]& 0\\ 0& \left[I_0\right]& 0& 0& \left[I_1\right]\\ 0& 0& \left[I_0\right]& 0& 0\\ \left[I_1\right]& 0& 0& \left[I_2\right]& 0\\ 0& \left[I_1\right]& 0& 0& \left[I_2\right]\end{array}\right]\left\{\begin{array}{c}\left\{{\ddot{u}}_0\right\}\\ \left\{{\ddot{v}}_0\right\}\\ \left\{{\ddot{w}}_0\right\}\\ \left\{{\ddot{\phi }}_1\right\}\\ \left\{{\ddot{\phi }}_2\right\}\end{array}\right\}=\left\{\begin{array}{c}\left\{F_1^T\right\}\\ \left\{F_2^T\right\}\\ \left[q{\delta }_{ij}\right]+\left\{F_3^T\right\}\\ \left\{F_4^T\right\}\\ \left\{F_5^T\right\}\end{array}\right\}$

47

where $K_{11}\dots K_{55}$ are the element stiffness matrices and are shown in Appendix 1; $\left\{F_1^T\right\}\dots \left\{F_5^T\right\}$ denote thermal resultants are shown in Appendix 1; $\left\{u_0\right\}\dots$ $\left\{{\phi }_2\right\}$ are displacement field vectors on the mid-surface that were defined in equations (45); $\left\{{\ddot{u}}_0\right\}\dots$ $\left\{{\ddot{\phi }}_2\right\}$ are the second derivative of $\left\{u_0\right\}\dots$ $\left\{{\phi }_2\right\}$ ; I₀, I₁, and I₂ are the mass inertias that were defined in equations (39); and δ_ij is the kronecker delta matrix.

Results and discussions of numerical analysis

In order to test the accuracy and efficiency of the developed algorithm, here we present numerical results for a number of sample problems. Then, the convergence study and further analysis are performed.

Two-layer $\left[-45°/45°\right]$ angle-ply clamped cylindrical shell under internal pressure

The following geometric parameters and lamina material properties were used in the study

$\begin{array}{l}{E}_1=7.5\times {10}^6\text{ psi};\hspace{0.17em}{E}_2=2\times {10}^6\text{ psi};\hspace{0.17em}{G}_{12}=G_{23}=G_{31}=1.25\times {10}^6\text{ psi};\\ {\vartheta }_{12}={\vartheta }_{12}=0.25;\hspace{0.17em}\rho =1\frac{{\text{lbs}}^2}{{\text{in}}^4};\hspace{0.17em}l=20\text{ in};\hspace{0.17em}h=1\text{ in};\hspace{0.17em}R=20\text{ in}\end{array}$

and subjected to an internal uniform step pulse (i.e. applied at time t = 0 and kept indefinitely) of intensity q₀ = 5000 psi is considered. The results of the present study are compared in Figure 4 with those obtained by Reddy and Chandrashekhara. ⁵² The results have good agreement with the DQM solution. This shows the correctness and reliability of the present formulation.

Figure 4.

Transient response of two-layer [−45°/45°] angle-ply clamped cylindrical shells under internal pressure (q₀ = 5000 psi, Δt = 0.005 s).

Two-layer laminated cylindrical shell with heat transfer analysis

To validate heat transfer analysis in the present study (heat transfer analysis with layerwise theory), comparisons with the exact solution of references. ^53,54 This example includes a two-layered cylinder with the following characteristics

$r_0=1, r_1=2, r_2=4, k_1=5, k_2=1, {\rho }_1={\rho }_2=c_{p1}=c_{p2}=1, T\left(r_0, t\right)=1, T\left(r_2, t\right)=0$

The temperature distribution versus layers position in direction of cylinder radius is shown in Figure 5 at t = 0.1 for the number of nodes across the thickness is considered to be N_l= 2, 4, 8, 10, 12, and 14. It is seen that the temperature distribution converges already with 16 grid points through the thickness of the cylinder.

Figure 5.

Temperature distribution in two-layered cylinder.

Laminated shell with metallic liner under thermoelastic load

In this example, a sample of metal composite cylindrical vessel is considered. A schematic of the vessel is shown in Figure 6. The vessel consists of five layers. The first layer is aluminum that its properties are listed in Table 1. This layer with thickness of 2 mm is used for vessel liner. Four other layers are used from the composite T700 that their properties are listed in Table 2. Thickness and angles of composite layers are (1.5, 0.5, 0.5, 2) millimeters and (90, 75, 75, 90) degrees, respectively. This vessel is used to a specific user. The cylinder is subjected to internal dynamic pressure and thermal load. Dynamic pressure has been obtained from an experimental test. For this work has been considered a metallic vessel. Some points from the vessel have been punctured for installation of pressure sensor adaptors. Then, with reinforcing connector parts, welding has been performed upon the shell by argon weld. The vessel dimensions are listed in Table 3. For loading or create pressure in the vessel has been utilized 800 g double-base propellant.

Figure 6.

Schematic of metal composite cylindrical vessel.

Table 1.

Thermal properties of Aluminum. ^55,56

Property	Unite	Value
Density	$\rho \left({\text{kg/m}}^{\text{3}}\right)$	2700
Specific heat capacity	$c_p\left(\frac{\text{J}}{\text{kg}\cdot °\text{k}}\right)$	900
Thermal expansion coefficient	$\alpha \left(\frac{1}{°\text{k}}{.10}^{-6}\right)$	24
Thermal conductivity	$k\left(\frac{\text{W}}{\text{m}\cdot °\text{k}}\right)$	210

Table 2.

Composite properties T700 with fiber volume fraction of 50%. ^{57 –60}

Property	Unite	Value
Density	$\rho \left(\text{kg}/{\text{m}}^3\right)$	1600
Specific heat capacity	$c_p\left(\frac{\text{J}}{\text{kg}\cdot °\text{k}}\right)$	800
Thermal expansion coefficient—longitudinal	${\alpha }_{11}\left(\frac{1}{°k}\times {10}^{-6}\right)$	3.1
Thermal expansion coefficient—transverse	${\alpha }_{22}\left(\frac{1}{°k}\times {10}^{-6}\right)$	24.7
Thermal conductivity—longitudinal	$k_{11}\left(\frac{W}{\text{m}\cdot °\text{k}}\right)$	7.46
Thermal conductivity—transverse	$k_{22}\left(\frac{W}{\text{m}\cdot °\text{k}}\right)$	0.56

Table 3.

Size of test vessel.

The inner radius(mm)	Thickness (mm)	Length (mm)
65(mm)	5.5	1200

The schematic arrangement of the sensors is shown in Figure 7. Piezoelectric sensors (kistler model 601H) have been utilized for dynamic pressure measurement. After preparation, experimental test was done that a sample of pressure–time diagram obtained from test results has been shown in Figure 8.

Figure 7.

The schematic arrangement of the sensors (dimensions in mm).

Figure 8.

Pressure–time diagram obtained from experimental test.

Vessel inside temperature is calculated from internal ballistic analysis of double-base propellant combustion that is shown in Figure 9 and vessel outside temperature is ambient temperature. Figure 10 shows a schematic of two-dimensional from the vessel. For calculation of free heat convection coefficient, value 10 $\frac{\text{W}}{\text{m}\cdot °\text{k}}$ is considered to ambient air and for hot gases of vessel inside is shown in Figure 11 that is calculated from references. ^{61 –64}

Figure 9.

Change of vessel inside temperature versus time.

Figure 10.

Schematic of two-dimensional from the cylindrical vessel and arrangement of layers.

Figure 11.

Change of convective heat transfer coefficient in heat gases versus time.

Boundary conditions are expressed as follows:

Boundary conditions for heat transfer analysis with layerwise theory: this type of boundary condition is as follows:

$\begin{array}{l}r=r_{Ne}:K_{Ne}\frac{\partial T}{\partial r}+h_{\text{out}}\left(T_{N_l}-T_{\infty }^{\text{out}}\right)=0\\ r=r_0:-K_1\frac{\partial T}{\partial r}+h_{\text{in}}\left(T_1-T_{\infty }^{\text{in}}\right)=0\end{array}$

48

The continuity conditions at the interfaces ( $r=r_{i, i=1,2,3,\dots .Ne})$ are

$\begin{array}{l}{T}_{i+1}\left(r_i,t\right)=T_i\left(r_i,t\right)\\ K_{i+1}\frac{\partial T_{i+1}\left(r_i,t\right)}{\partial r}=K_i\frac{\partial T_i\left(r_i,t\right)}{\partial r}\end{array}$

49

where h_out and h_in are the heat convection coefficient for ambient air and heat gases (vessel outside and vessel inside), respectively; K_Ne and K₁ are the heat conduction coefficient for up and down layers in cylindrical vessel; and $T_{\infty }^{\text{out}}$ and $T_{\infty }^{\text{in}}$ are temperature of ambient air and heat gases, respectively. Heat transfer by radiation is very small and it may be ignored. By substituting the boundary conditions (equations (48) to (49)) into the thermal governing equations (equations (35)), a system of algebraic equations with $(M_{\varnothing }+1)(N_z+1)(N_e+1)N_t$ equations and unknowns is derived. Then, the time histories of the temperature at each grid point are determined by solving the system of equations.

Boundary conditions for thrmoelastic analysis:

Clamped boundary condition: this type of boundary condition is as follows:

$u_0=v_0=w_0={\phi }_1={\phi }_2=0$

50

Simply supported boundary condition: this type of boundary condition is as follows:

$w_0=M_1=N_6=N_1={\phi }_2=0$

51

By substituting the boundary conditions (equations (50) to (51)) into the thermoelastic governing equations (equations (47)), a system of algebraic equations with $5\left(M+1\right)\left(N+1\right)N_t$ equations and unknowns is derived. Then, the time histories of the displacement field at each grid point are determined by solving the system of equations. On the other hand, in this example, a metal composite cylindrical vessel is modeled using the commercial finite element software of ABAQUS. The results of finite element software are compared with the DQM solution. Therefore, in this stage, the eight types of analysis is done: (1) thermoelastic load with simple support and DQM (TELSS-DQM), (2) thermoelastic load with clamp support and DQM (TELCS-DQM), (3) elastic load with simple support and DQM (ELSS-DQM), (4) elastic load with clamp support and DQM (ELCS-DQM), (5) thermoelastic load with simple support and finite element method by ABAQUS (TELSS-ABAQUS), (6) thermoelastic load with clamp support and finite element method by ABAQUS (TELCS-ABAQUS), (7) elastic load with simple support and finite element method by ABAQUS (ELSS-ABAQUS), and (8) elastic load with clamp support and finite element method by ABAQUS (ELCS-ABAQUS). The convergence behavior of the temperature field is examined in Figure 12. The grid points in the axial direction are changed from 5 to 20 while the number of nodes is considered along the thickness direction N_e= 6, 8, 10, 12, 14, and 16. It is observed that the temperature is converges with 12 grid points through the thickness of the vessel. Therefore, the density of grid points in the DQM analysis is suitable for the temperature components with N_e = 12, $N_{x_1}=N_z$ = 15. As well as, the mentioned analysis is performed to the displacement and stress components that the convergence behavior is obtained similar to the above results but it is not repeated here. The convergence behavior of the time step is examined and time step equal to 0.01 ms is selected. The grid points in the ∅-direction similar to z-direction is achieved in M_∅ = 12. For summarizing result, these results are not shown here.

Figure 12.

The convergence behavior of temperature field for DQM code at t=13 ms.

Figure 13 shows time histories of the radial displacement at the mid-surface of the cylinder with l = 300 mm for the two described boundary conditions. The results achieved by the DQM method are compared with those obtained by the ABAQUS software. It is observed that the distribution of the radial displacement near the ends exhibits a very low dependency on the type of boundary condition. As shown in figure, the comparison shows good agreement between the DQM algorithm and ABAQUS software. It is also seen that the radial displacement in the simply supported is about 20% more than from the clamped boundary condition in the maximum value.

Figure 13.

Time histories of the radial displacement at the mid-surface of the cylinder with l = 300 mm obtained by TELSS-DQM, TELCS-DQM, TELCS-ABAQUS, and TELSS-ABAQUS.

Figure 14 illustrates time histories of the von Mises stress at the liner with l = 300 mm for the simply supported boundary conditions and two types of loading. The results achieved by the DQM method are compared with those obtained by the ABAQUS software. The current solutions show a very good agreement with the ABAQUS results. It is also observed that the stress in the TELSS loading is between 0% and 50% more than from the ELSS loading.

Figure 14.

Time histories of the Von Mises stress at the liner with l = 300 mm obtained by TELSS-DQM, ELCS-DQM, TELSS-ABAQUS, and ELSS-ABAQUS.

Figure 15 demonstrates time histories of the fiber stress in the circumferential direction with l = 300 mm for the simply supported boundary conditions and two types of loading. It can be seen that the normal stress in the TELSS loading is between 0% and 34% more than from the ELSS loading. It is also observed that maximum stress for both types of loading appears at θ = 90°.

Figure 15.

Time histories of the fiber stress in the circumferential direction with l = 300 mm obtained by TELSS-DQM and ELSS-DQM.

Figure 16 depicts time histories of the longitudinal stress at the liner with l = 0, 300, and 1200 mm for the simply supported boundary conditions and two types of loading. It can be seen that the stress in the TELSS loading is between 0% and 85% more than from the ELSS loading. It is also seen that the longitudinal stress is about half the amount of von Mises stress in the maximum value.

Figure 16.

Time histories of the longitudinal stress at the liner with l = 0, 300, and 1200 mm obtained by TELSS-DQM and ELSS-DQM.

Time histories of the tangential stresses at the liner for the two described boundary conditions with = 0, 300, and 1200 mm are shown in Figure 17. The results predicted by the DQM algorithm are compared with those obtained by the ABAQUS software. It is observed that the distribution of the tangential stress near the ends exhibits a very low dependency on the type of boundary condition. As shown in Figure 17, the comparison shows good agreement between the DQM algorithm and ABAQUS software. It is also seen that the tangential stress in the clamped boundary conditions is nearly twice more than from the simply supported boundary conditions in the maximum value.

Figure 17.

Time histories of the tangential stress in liner with (a) l = 300 mm and (b) l = 0 and 1200 mm obtained by TELCS-DQM, TELSS-DQM, TELSS-ABAQUS, and TELCS-ABAQUS.

The tangential and shear stresses with l = 300 mm in the thickness direction are shown in Figure 18 by TELSS and ELSS. The results predicted by the DQM algorithm are compared with those obtained by the ABAQUS software. It is observed that the distribution of the circumferential stress in 90° layers is ascending but in 75° layers is descending to mid-plane (the reference plane). Amount of shear stresses at the top and bottom of the laminate equal to zero. Note that interlaminar stresses exhibit apparent singular behavior near the intersection of the layers interfaces. The stress distributions in Figure 18 show good qualitative agreement between the DQM algorithm with results ABAQUS software.

Figure 18.

The stress with l = 300 mm in the thickness direction (a) tangential stress and (b) shear stress obtained by TELSS-DQM, ELSS-DQM, TELSS-ABAQUS, and ELSS-ABAQUS.

In Table 4, the stress components as well as the radial displacement for liner at t = 7.5 ms are obtained and compared with ABAQUS results for various types of boundary conditions and loading. The solutions presented show a good agreement between the DQM method and ABAQUS results. In Tables 5 to 8, the stress components for the 1th to 4th fiber-reinforced layer at t = 7.5 ms and l = 300 mm are presented for various types of boundary conditions and loading. The results of the DQM code are compared with those obtained by ABAQUS.

Table 4.

Amounts of the stress components as well as the radial displacement for liner at t = 7.5 ms and l = 300 mm obtained by DQM and ABAQUS.

	TELCS-DQM	TELCS-ABAQUS	TELSS-DQM	TELSS-ABAQUS	ELCS-DQM	ELCS-ABAQUS	ELSS-DQM	ELSS-ABAQUS
$S_{\text{von}}\left(\text{MPa}\right)$	141.28	139.12	139.36	135	71.4	69.36	63.08	61.68
${\text{S}}_l\left(\text{MPa}\right)$	337.88	334.24	329.8	324.8	221.24	217.04	216.04	213.68
${\text{S}}_{\text{t}}\left(\text{MPa}\right)$	406.28	390.76	379.24	368.2	284.08	279.56	263.24	252.2
$w\hspace{0.17em}\left(\text{mm}\right)$	0.1896	0.18552	0.24	0.23472	0.1334	0.13152	0.1723	0.1701

DQM: differential quadrature method.

Table 5.

Amounts of the stress components for the 1th fiber-reinforced layer at t = 7.5 ms and l = 300 mm obtained by DQM and ABAQUS.

	TELCS-DQM	TELCS-ABAQUS	TELSS-DQM	TELSS-ABAQUS	ELCS-DQM	ELCS-ABAQUS	ELSS-DQM	ELSS-ABAQUS
$S_{11}\left(\text{MPa}\right)$	666.16	655.8	622.8	613.8	521	517.36	485.44	479.24
$S_{22}\left(\text{MPa }\right)$	8.32	7.92	7.72	7.4	6.68	6.44	6.32	6.08
$S_{12}\left(\text{MPa }\right)$	0.392	0.384	0.38	0.364	0.156	0.14	0.128	0.116

DQM: differential quadrature method.

Table 6.

Amounts of the stress components for the 2th fiber-reinforced layer at t = 7.5 ms and l = 300 mm obtained by DQM and ABAQUS.

	TELCS-DQM	TELCS-ABAQUS	TELSS-DQM	TELSS-ABAQUS	ELCS-DQM	ELCS-ABAQUS	ELSS-DQM	ELSS-ABAQUS
$S_{11}\left(\text{MPa}\right)$	609.72	604.52	577.72	574.6	501.68	479.92	452.68	447.68
$S_{22}\left(\text{MPa }\right)$	11.24	11.04	8.72	8.64	10.84	10.76	7.16	6.96
$S_{12}\left(\text{MPa }\right)$	10.24	10.04	8.32	8.12	6.72	6.6	6.48	6.4

DQM: differential quadrature method.

Table 7.

Amounts of the stress components for the 3th fiber-reinforced layer at t = 7.5 ms and l = 300 mm obtained by DQM and ABAQUS.

	TELCS-DQM	TELCS-ABAQUS	TELSS-DQM	TELSS-ABAQUS	ELCS-DQM	ELCS-ABAQUS	ELSS-DQM	ELSS-ABAQUS
$S_{11}\left(\text{MPa}\right)$	573.8	569.72	532.12	529.36	453.68	439.56	400.48	399.92
$S_{22}\left(\text{MPa }\right)$	11.12	10.6	8.68	8.6	10.76	10.28	7.12	6.92
$S_{12}\left(\text{MPa }\right)$	−10.16	−9.72	−8.2	−8.04	−6.76	−6.64	−6.56	−6.48

DQM: differential quadrature method.

Table 8.

Amounts of the stress components for the 4th fiber-reinforced layer at t = 7.5 ms and l = 300 mm obtained by DQM and ABAQUS.

	TELCS-DQM	TELCS-ABAQUS	TELSS-DQM	TELSS-ABAQUS	ELCS-DQM	ELCS-ABAQUS	ELSS-DQM	ELSS-ABAQUS
$S_{11}\left(\text{MPa}\right)$	666.12	655.76	622.84	613.8	524.04	517.4	485.4	479.32
$S_{22}\left(\text{MPa }\right)$	8.36	7.96	7.72	7.44	6.72	6.48	6.32	6.12
$S_{12}\left(\text{MPa }\right)$	0.388	0.384	0.38	0.368	0.152	0.136	0.128	0.116

DQM: differential quadrature method.

All solutions show close agreements. The total CPU time to thermoelastic analysis of considered laminated shell with metallic liner is about 80 s using the DQM method. On the other hand, the results of the finite element solutions using the ABAQUS software use 36 times more CPU time than the DQM code for the same problem. Therefore, in comparison with the simulation of finite element method that require very fine meshes, the DQM method is an efficient way for numerical analysis due to less computational cost.

Conclusion

Comparison of thermoelastic stress with elastic stress in metal composite vessel is investigated under experimental dynamic load. Elastic analysis of laminated cylindrical shells is done by using FSDT. Thermoelastic analysis of laminated cylindrical shells is done by using the mixed layerwise and FSDTs with DQM. DQM is used to discretize the resulted governing equations while Newmark’s method is implemented to solve the problem in the time domain. The MATLAB programming code has been used to solve differential quadrature equations. It was observed that the fast convergence rate of the DQM method and its accuracy is confirmed by comparing the results with those obtained using ABAQUS software and also with the exact solutions. It was achieved significantly less computational efforts of the proposed approaches with respect to the finite element method. Also, calculation of temperature field with layerwise theory provides the ability to model accurately of laminated composite shells. Furthermore, comparative studies are done between two methods of elastic and thermoelastic stress, it is found that about 48% in the maximum amount higher stress in the thermoelastic state for the metal liner and 32% for the composite layers. Since the results of this research are derived from the experimental load can be used as criterion solutions for future studies.