A Cylindrical Superelement for Thermo-Mechanical Analysis of Thin Composite Vessels

A. Jafarzadeh, A. Taghvaeipour^* and M. R. Eslami

Department of Mechanical Engineering, Amirkabir University of Technology (Tehran Polytechnic), Hafez Ave., Tehran, Iran

E-mail: afshin.taghvaeipour@gmail.com

*Corresponding Author

Received 31 July 2020; Accepted 02 November 2020; Publication 04 January 2021

Abstract

In this study, a new cylindrical shell superelement with trigonometric shape functions is developed. This element is formulated based on the classical theory of shells, and it is especially designed for coupled-field analysis of thin cylindrical vessels or tubes made of composite materials. As a case study, a thermo-mechanical analysis of a thin composite cylinder is conducted. By invoking to the uniform and non-uniform meshing, the deformation and the stress results are calculated and compared with the analytical solutions. At the end, the efficiency and accuracy of the proposed superelement is also depicted via comparison of the corresponding results with the ones which are calculated by means of shell elements and via a commercial software package.

Keywords: Classical theory of shell, cylindrical vessels, trigonometric shape functions, shell superelement, composite cylinders.

1 Introduction

During the last decades, many researchers have studied mechanical and thermal behavior of the composite cylinders, which may be used in pressure vessels and tubes. In this regard, Xia et al. [1] used the three-dimensional anisotropic elasticity theory to investigate the influence of stacking sequence on the mechanical properties of a multi-layered filament-wound composite pipe with three different angle-ply and in the presence of an internal pressure. The authors also investigated the shear twisting coupling characteristic. The same authors used the classical laminated-plate theory and calculated the stress and deformation of laminated sandwich pipes with different winding angles and under combined internal pressure and temperature differences [2]. By resorting to the three-dimensional anisotropic elasticity theory, Bakayan et al. [3] conducted the thermal analysis of a composite pipe under the combined internal pressure and temperature gradient loadings. Xu et al. [4] predicted the failure strength of composite pipes by incorporating different failure criteria such as the maximum stress, Tsai-Hill, and etc. In fact, they used a three dimensional parametric FE model of the cylindrical part to explore the non-linear stress–strain relationship, and accordingly, the final failure pressure. Pagans et al. [5] used the classical lamination theory and generalized plane strain model to design and predict the behavior of fiber reinforced composite pressure vessels. The authors applied mechanical, thermal, and moisture loadings on the body, and finally, they obtained the optimum values for the winding angle, burst pressure, maximum axial force, and the maximum angular speed of the pressure vessel. Numerical and experimental studies on the filament wound composite tubes subjected to an internal pressure was also conducted by Martins et al. [6]. The numerical and experimental results were properly compatible for all fiber directions. While most of the reported researches were focused on the stress analysis of thick hollow cylinders, Sayman [7] developed a general stress analysis of thick and thin multi-layered composite cylinders under hygrothermal loadings. The solution was obtained for composite cylinders with symmetrical and antisymmetrical layers and for open and closed end conditions. The thermo-mechanical behavior of structures by analytical methods is, in general, complicated. However, with some simplifications and for a number of classical problems, the analytical solutions are available [8, 9]. Some researchers also developed analytical methods to cope with the thermo-mechanical analysis of cylindrical vessels made of functionally graded materials (FGM) [10, 11]. However, these works are mostly developed for thick vessels or tubes which are made of orthotropic materials.

The finite element method is proved to be a reliable and strong numerical approach in modeling and analysis of complex structures. In this method, a structure is discredited into many elements with simple geometries. Due to recent developments of the computational geometry, fairly advanced discretization methods have been proposed, and hence, with some effort, any complex geometry can be meshed. However, this usually leads to a bulk of nodes with many degrees-of-freedom (dof) which can be problematic in large systems. This can be even more severe when a nonlinear time-dependent analysis should be conducted. Accordingly, the accuracy of the method is mostly relied on the discritization quality. In the recent years, new type of large elements with trigonometric shape functions are introduced which are called superelements. These elements are especially designed for time-consuming iterative analyses of which less number of dof is demanded. These types of elements were first introduced by Koko et al., and called superelements [12]. Koko and Olson first used the super plate and beam finite elements to calculate the natural frequency and mode shapes of stiffened rectangular plates [13]. They validated their results with other numerical and experimental works. Jiang and Olson [14, 15] introduced three shell, curved beam, and straight beam superelements for the free vibration analysis of stiffened cylindrical shells. Moreover, the superelements were also incorporated in the static and dynamic analysis of orthogonally stiffened cylindrical shells. In this work, the stiffeners were assumed to be in form of curved and straight beams. In addition, the results including natural frequency and mode shapes were compared with the ones which were obtained by the conventional finite elements and experimental tests. As a result, high rate of convergence and accuracy were reported. In [16, 17] Ahmadian et al. employed a plate superelement to analyze the forced and free vibrations of an orthotropic rectangular plate with different boundary conditions. They considered bending and in-plane effects and assumed $C^0$ continuity for in-plane and $C^1$ continuity for out-of plane displacements. In one of the recent works, Ju and choo used a superelement to conduct static analysis of a cable passing through multiple pulleys [18]. In this paper deformations and stresses were calculated by superelements and compared with experimental tests; the results showed a good agreement. Ahmadian et al. [19] introduced a solid cylindrical superelement for the static and modal analyses of laminated hollow cylinders. Accordingly, in [20], Taghvaeipour et al. incorporated the foregoing element for the modal analysis of a thick hollow cylinder made of FGM. The reported results show a good accuracy with those reported in literature. Using the same concept, Ahmadian et al. [21, 22, 23] have defined the spherical and tapered super elements. With these new elements, they conducted structural analysis of spherical components such as biologic cells, nano bearings, pressure vessels, and fullerene. Recently, Pourhamid et al. [24] employed the cylindrical superelement to conduct the stress analysis on a cylinder-piston of a combustion engine which is made of FGM.

In this study, based on the classical theory of shell, a cylindrical superelement for thin cylinders, made of composite materials, is introduced. Comparing with other superelements in the literature, the introduced element is designed for shell structures in which one level of nodes exists through the thickness. Moreover, the radial deformation slope is considered as a continuous degree of freedom. In other words, this superelement is $C^1$ continuous in which the slop continuity condition results in a more accurate strain and stress calculation. The superelements introduced in the literature are often resorted on $C^0$ Continuous shape functions in which only the displacements are continuous, and thus, they demand higher number of elements for better accuracy in the calculation of strains and stresses. As a case study, the defined element is incorporated on a thin composite cylinder which is under thermo-mechanical loadings. The obtained results are first compared with the analytical solutions which were reported in the literature, and next with the ones obtained by shell elements in a commercial software package.

2 Element Definition

2.1 Geometry and Coordinate System of Element

In general, shell structures are divided into two types: thin- and thick-walled. Here, we consider the superelement for a thin shell with a cylindrical geometry. The element has a mid-surface radius of $R$, thickness of $t$, and the length of $2L$. The word ”thin-walled” means, that the ratio $t/R$ is larger than 1/1000 and less than 1/20 [25]. The thickness $t$ is considered constant all over the element. To derive the corresponding equations, the cylindrical coordinate system is used which is located at the center of the element. The nodes are distributed uniformly along two edges of the element and at the mid-surface. Geometrical properties, coordinates, and nodes are shown in Figure 1. The cylindrical global coordinate system is selected with $r$, $\alpha$, and $z$ components, namely [12],

$$\{\begin{array}{cc}0\le \alpha \le 2\pi \to -1\le \gamma \le 1\hfill & \\ -L\le z\le L\to -1\le \zeta \le 1\hfill & \end{array}$$

(1)

where

$$\gamma =\frac{\alpha }{\pi }-1,\zeta =\frac{z}{L}$$

(2)

Figure 1 Cylindrical Shell Superelement.

2.2 Element Degrees of Freedom and Shape Functions

As it was mentioned, the element is defined based on the classical shell theory in which the radial displacement component is independent of the radial coordinate [25, 26]. In this regard, the displacements in an arbitrary point of the shell are defined as [27],

	$w=w^0,$	(3)
	$\nu ={\nu }^0+(r-R){\beta }_{\alpha }$
	$u=u^0+(r-R){\beta }_z$

where $\nu$, $u$, and $w$ are the tangential, axial, and radial displacements, respectively, and ${\nu }^0$, $u^0$, and $w^0$ are the corresponding displacements at the mid-surface of the shell. Also, ${\beta }_{\alpha }$ and ${\beta }_z$ are rotations in the tangential and axial directions, respectively. Equation 3 shows that $\nu$ and $u$ depend on the radial coordinate $r$, where $w$ is independent of it. The rotations are also defined as [27]:

$${\beta }_{\alpha }=\frac{{\nu }^0}{R}-\frac{1}{R}\frac{\partial w}{\partial \alpha },{\beta }_z=-\frac{\partial w}{\partial z}$$

(4)

Since the rotations depend on the first derivative of $w$, it should be approximated by $C^1$-continuous shape functions. However, other displacements need $C^0$-continuous shape functions [28, 29]. Therefore, the displacement vector at a point at the mid-surface is defined as

$$𝒖={\left[\begin{array}{cccccc}\hfill w\hfill & \hfill \frac{\partial w}{\partial z}\hfill & \hfill \frac{\partial w}{\partial \alpha }\hfill & \hfill \frac{{\partial }^{2}w}{\partial z\partial \alpha }\hfill & \hfill {\nu }^0\hfill & \hfill u^0\hfill \end{array}\right]}^T$$

(5)

where the first four degrees of freedom need to be approximated by $C^1$-continuous shape functions. Hence, the corresponding shape function will have $32$ constants (four times the number of nodes). Therefore, the radial displacement is approximated as [28, 29]:

$w(z,\alpha )$	$=$	$N_1{w}_1+{N_1^{\prime }{\displaystyle \frac{\partial w}{\partial z}}|}_1+{N_1^{\prime \prime }{\displaystyle \frac{\partial w}{\partial \alpha }}|}_1+{N_1^{\prime \prime \prime }{\displaystyle \frac{{\partial }^{2}w}{\partial \alpha \partial z}}|}_1+\mathrm{\cdots }+N_8{w}_8$	(6)
		$+{N_8^{\prime }{\displaystyle \frac{\partial w}{\partial z}}|}_8+{N_8^{\prime \prime }{\displaystyle \frac{\partial w}{\partial \alpha }}|}_8+{N_8^{\prime \prime \prime }{\displaystyle \frac{{\partial }^{2}w}{\partial \alpha \partial z}}|}_8$

where $N_i$, $N_i^{\prime }$, $N_i^{\prime \prime }$, and $N_i^{\prime \prime \prime }$ are $C^1$-continuous shape functions. Accordingly, the approximation function is assumed as a combination of trigonometry and polynomial functions, namely,

$w(z,\alpha )=$	$\left(a_1+a_{2}z+a_3{z}^2+a_4{z}^3\right)$	(7)
	$\times (b_1+b_2\sin \alpha +b_3\cos \alpha +b_4\sin 2\alpha +b_5\cos 2\alpha +b_6\sin 3\alpha$
	$+b_7\cos 3\alpha +b_8\cos 4\alpha )$

By invoking the necessary conditions on the approximation function, the constants are calculated. For the sake of brevity, the shape functions are presented in the appendix A. Moreover, the tangential and axial displacements at the mid surface can be approximated with $C^0$-continuous identical functions, namely,

$$\nu =u=\left(a_1+a_{2}z\right)\left(b_1+b_2\cos \alpha +b_3\sin \alpha +b_4\cos 2\alpha \right)$$

(8)

Likewise, the eight constants can be calculated by applying the necessary conditions on the approximation function. As a result, the displacement function can be written in terms of nodal displacements and shape functions as

${\nu }^0(z,\alpha )$	$=M_1{\nu }_1^0+M_2{\nu }_2^0+\mathrm{\dots }+M_8{\nu }_8^0$	(9)
$u^0(z,\alpha )$	$=M_1{u}_1^0+M_2{u}_2^0+\mathrm{\dots }+M_8{u}_8^0$

in which ${\nu }_i^0$ and $u_i^0$ are tangential and axial nodal displacements, respectively. The shape functions $M_i$ are given in the appendix A. In total, this element has $32+8=40$ shape functions. By resorting to the FEM notations, the vector $𝐮$, which is the vector of displacement components (Equation 5), can be obtained as follows

$$𝐮={𝐍}_{𝐮}𝐪$$

(10)

where $𝐪={\left[\begin{array}{ccc}\hfill {𝐪}_1\hfill & \hfill \mathrm{\dots }\hfill & \hfill {𝐪}_8\hfill \end{array}\right]}^T$ is the nodal displacement vector in which

$${𝐪}_i=\left[\begin{array}{cccccc}\hfill w_i\hfill & \hfill \frac{\partial w_i}{\partial z}\hfill & \hfill \frac{\partial w_i}{\partial \alpha }\hfill & \hfill \frac{\partial w_i}{\partial \alpha \partial z}\hfill & \hfill {\nu }_i^0\hfill & \hfill u_i^0\hfill \end{array}\right]$$

(11)

Also, ${𝐍}_u$ is the matrix of shape functions which is defined as

$${𝐍}_u=\left[\begin{array}{ccccc}\hfill {𝐍}_1\hfill & \hfill 𝐎\hfill & \hfill \mathrm{\dots }\hfill & \hfill {𝐍}_8\hfill & \hfill 𝐎\hfill \\ \hfill {𝐎}^T\hfill & \hfill {𝐌}_1\hfill & \hfill \mathrm{\dots }\hfill & \hfill {𝐎}^T\hfill & \hfill {𝐌}_8\hfill \end{array}\right]$$

(12)

where $𝐎$ is a $4\times 2$ zero matrix and

$${𝐍}_i=\left[\begin{array}{cccc}\hfill N_i\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill N_i^{\prime }\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill N_i^{\prime \prime }\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill N_i^{\prime \prime \prime }\hfill \end{array}\right],{𝐌}_i=\left[\begin{array}{cc}\hfill M_i\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill M_i\hfill \end{array}\right]$$

(13)

3 The FEM Formulations

3.1 Kinematic Relations

In the classical shell theory, the shear strains ${\gamma }_{r\alpha }$ and ${\gamma }_{rz}$ and the radial strain ${ϵ}_{rr}$ are neglected. The remained strains are functions of strains at the mid-surface and the curvatures, namely [27]:

${ϵ}_z$	$={ϵ}_z^0+\left(r-R\right){\chi }_z$	(14)
${ϵ}_{\alpha }$	$={ϵ}_{\alpha }^0+\left(r-R\right){\chi }_{\alpha }$
${\gamma }_{\alpha z}$	$={\gamma }_{\alpha z}^0+2\left(r-R\right){\chi }_{\alpha z}$

where ${ϵ}_z^0$, ${ϵ}_{\alpha }^0$, and ${\gamma }_{\alpha z}^0$ are strains at the mid-surface which are functions of displacements [27]:

$${ϵ}_z^0=\frac{\partial u^0}{\partial z},{ϵ}_{\alpha }^0=\frac{1}{R}\frac{\partial {\nu }^0}{\partial \alpha }+\frac{w}{R},{\gamma }_{\alpha z}^0=\frac{1}{R}\frac{\partial u^0}{\partial \alpha }+\frac{\partial {\nu }^0}{\partial z}$$

(15)

in which ${\chi }_z$ , ${\chi }_{\alpha }$, and ${\chi }_{\alpha z}$ are components of the shell curvatures that are defined below [27]:

${\chi }_z$	$={\displaystyle \frac{\partial {\beta }_z}{\partial z}}=-{\displaystyle \frac{{\partial }^{2}w}{\partial z^2}}$	(16)
${\chi }_{\alpha }$	$={\displaystyle \frac{1}{R}}{\displaystyle \frac{\partial {\beta }_{\alpha }}{\partial \alpha }}={\displaystyle \frac{1}{R}}{\displaystyle \frac{\partial }{\partial \alpha }}\left({\displaystyle \frac{{\nu }^0}{R}}-{\displaystyle \frac{1}{R}}{\displaystyle \frac{\partial w}{\partial \alpha }}\right)$
${\chi }_{\alpha z}$	$={\displaystyle \frac{\partial {\beta }_{\alpha }}{\partial z}}+{\displaystyle \frac{1}{R}}{\displaystyle \frac{\partial {\beta }_z}{\partial \alpha }}={\displaystyle \frac{1}{R}}{\displaystyle \frac{\partial }{\partial z}}\left({\nu }^0-2{\displaystyle \frac{\partial w}{\partial \alpha }}\right)$

Finally, the strains and shell curvatures can be presented in terms of nodal displacements, namely,

$${\left[\begin{array}{cccccc}\hfill {ϵ}_{zz}^0\hfill & \hfill {ϵ}_{\alpha \alpha }^0\hfill & \hfill {\gamma }_{z\alpha }^0\hfill & \hfill {\chi }_{zz}\hfill & \hfill {\chi }_{\alpha \alpha }\hfill & \hfill {\chi }_{z\alpha }\hfill \end{array}\right]}^T={\mathrm{𝐋𝐍}}_u𝐪=\mathrm{𝐁𝐪}$$

(17)

in which $𝐋$ is an $6\times 6$ operator matrix which is defined as

$$𝐋=\left[\begin{array}{cccccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial }{\partial z}\hfill \\ \hfill \frac{1}{R}\hfill & \hfill \frac{1}{R}\hfill & \hfill \frac{1}{R}\hfill & \hfill \frac{1}{R}\hfill & \hfill \frac{1}{R}\frac{\partial }{\partial \alpha }\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial }{\partial z}\hfill & \hfill \frac{1}{R}\frac{\partial }{\partial \alpha }\hfill \\ \hfill -\frac{{\partial }^2}{\partial z^2}\hfill & \hfill -\frac{{\partial }^2}{\partial z^2}\hfill & \hfill -\frac{{\partial }^2}{\partial z^2}\hfill & \hfill -\frac{{\partial }^2}{\partial z^2}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -\frac{1}{R^2}\frac{{\partial }^2}{\partial {\alpha }^2}\hfill & \hfill -\frac{1}{R^2}\frac{{\partial }^2}{\partial {\alpha }^2}\hfill & \hfill -\frac{1}{R^2}\frac{{\partial }^2}{\partial {\alpha }^2}\hfill & \hfill -\frac{1}{R^2}\frac{{\partial }^2}{\partial {\alpha }^2}\hfill & \hfill \frac{1}{R^2}\frac{{\partial }^2}{\partial {\alpha }^2}\hfill & \hfill 0\hfill \\ \hfill -\frac{1}{R^2}\frac{{\partial }^2}{\partial \alpha \partial z}\hfill & \hfill -\frac{1}{R^2}\frac{{\partial }^2}{\partial \alpha \partial z}\hfill & \hfill -\frac{1}{R^2}\frac{{\partial }^2}{\partial \alpha \partial z}\hfill & \hfill -\frac{1}{R^2}\frac{{\partial }^2}{\partial \alpha \partial z}\hfill & \hfill \frac{1}{2R}\frac{\partial }{\partial z}\hfill & \hfill 0\hfill \end{array}\right]$$

(18)

3.2 Stress-strain relations

The stress-strain relation for the $k^{th}$ layer of a composite cylindrical shell with the thermo-mechanical effects is defined in [30] as:

$${\left[\begin{array}{c}\hfill {\sigma }_{\alpha \alpha }\hfill \\ \hfill {\sigma }_{zz}\hfill \\ \hfill {\tau }_{r\alpha }\hfill \\ \hfill {\tau }_{zr}\hfill \\ \hfill {\tau }_{\alpha z}\hfill \end{array}\right]}_k={\overline{𝐐}}_k\left[\begin{array}{c}\hfill {ϵ}_{\alpha \alpha }-{\alpha }_{\alpha \alpha }\mathrm{\Delta }T-{\beta }_{\alpha \alpha }c\hfill \\ \hfill {ϵ}_{zz}-{\alpha }_{zz}\mathrm{\Delta }T-{\beta }_{zz}c\hfill \\ \hfill {\gamma }_{r\alpha }\hfill \\ \hfill {\gamma }_{zr}\hfill \\ \hfill {\gamma }_{\alpha z}-2{\alpha }_{\alpha z}\mathrm{\Delta }T-2{\beta }_{\alpha z}c\hfill \end{array}\right]$$

(19)

where ${\sigma }_{ii}$ are the normal stresses, ${\tau }_{ij}$ are the shear stresses, ${\alpha }_{ij}$ are the heat expansion coefficients, $\mathrm{\Delta }T$ is the temperature gradient, and ${\beta }_{ij}$ and $c$ are the hygroscopic expansion and the moisture concentration coefficients, respectively. The mechanical properties matrix is also expressed as [30]:

$${\overline{𝐐}}_k={\left[\begin{array}{ccccc}\hfill {\overline{Q}}_{11}\hfill & \hfill {\overline{Q}}_{12}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {\overline{Q}}_{16}\hfill \\ \hfill {\overline{Q}}_{12}\hfill & \hfill {\overline{Q}}_{22}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {\overline{Q}}_{26}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\overline{Q}}_{44}\hfill & \hfill {\overline{Q}}_{45}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\overline{Q}}_{45}\hfill & \hfill {\overline{Q}}_{55}\hfill & \hfill 0\hfill \\ \hfill {\overline{Q}}_{16}\hfill & \hfill {\overline{Q}}_{26}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {\overline{Q}}_{66}\hfill \end{array}\right]}_k$$

(20)

The components of matrix ${\overline{𝐐}}_k$ for each layer and with the specific angle ply are given in [30].

3.3 Stiffness matrix and force vectors

Finally, a system of linear equations is obtained as follows

$${\left[{𝐊}_{uu}\right]}^e{𝐮}^e={𝐅}^e$$

(21)

where ${\left[{𝐊}_{uu}\right]}^e$ yield from the membrane strain energy of shell and can be obtained from the relation:

$${\left[{𝐊}_{uu}\right]}^e=RL\pi {\int }_1^{-1}{\int }_1^{-1}{\left[{𝐁}_u(\zeta ,\gamma )\right]}^T\left[𝐃\right]\left[{𝐁}_u(\zeta ,\gamma )\right]𝑑\zeta 𝑑\gamma$$

(22)

in which ${𝐁}_u(\zeta ,\gamma )$ is the nodal strain matrix-displacement matrix which is equal to:

$${𝐁}_u=𝐋\times {𝐍}_u$$

(23)

The operator matrix $𝐋$ was defined in Equation 18. Also, ${𝐃}_c$ is the mechanical properties matrix [30], namely,:

$${𝐃}_c=\left[\begin{array}{cccccc}\hfill A_{11}\hfill & \hfill A_{12}\hfill & \hfill A_{16}\hfill & \hfill B_{11}\hfill & \hfill B_{12}\hfill & \hfill 2B_{16}\hfill \\ \hfill A_{12}\hfill & \hfill A_{22}\hfill & \hfill A_{26}\hfill & \hfill B_{12}\hfill & \hfill B_{22}\hfill & \hfill 2B_{26}\hfill \\ \hfill A_{16}\hfill & \hfill A_{26}\hfill & \hfill A_{66}\hfill & \hfill B_{16}\hfill & \hfill B_{26}\hfill & \hfill 2B_{66}\hfill \\ \hfill B_{11}\hfill & \hfill B_{12}\hfill & \hfill B_{16}\hfill & \hfill D_{11}\hfill & \hfill D_{12}\hfill & \hfill 2D_{16}\hfill \\ \hfill B_{12}\hfill & \hfill B_{22}\hfill & \hfill B_{26}\hfill & \hfill D_{12}\hfill & \hfill D_{22}\hfill & \hfill 2D_{26}\hfill \\ \hfill B_{16}\hfill & \hfill B_{26}\hfill & \hfill B_{66}\hfill & \hfill D_{16}\hfill & \hfill D_{26}\hfill & \hfill 2D_{66}\hfill \end{array}\right]$$

(24)

where

$A_{ij}$	$={\displaystyle \sum _{k=1}^{N_l}}{\left({\overline{Q}}_{ij}\right)}_k\left(h_k-h_{k+1}\right)\mathit{\hspace{1em}}i,j=1,2,6$	(25)
$B_{ij}$	$={\displaystyle \frac{1}{2}}{\displaystyle \sum _{k=1}^{N_l}}{\left({\overline{Q}}_{ij}\right)}_k\left(h_k^2-h_{k+1}^2\right)\mathit{\hspace{1em}}i,j=1,2,6$
$D_{ij}$	$={\displaystyle \frac{1}{3}}{\displaystyle \sum _{k=1}^{N_l}}{\left({\overline{Q}}_{ij}\right)}_k\left(h_k^3-h_{k+1}^3\right)\mathit{\hspace{1em}}i,j=1,2,6$

In Eqs. (18), (19), and (20), $h_i$ is the $i^{t}h$ layer thickness, which is depicted in Figure 2. The force vectors, including the thermal and moisture effects, yield as [30]:

$F_T^e$	$=RL\pi \mathrm{\Delta }T{\displaystyle {\int }_1^{-1}}{\displaystyle {\int }_1^{-1}}{\left[{𝐁}_u(\zeta ,\gamma )\right]}^T\left[\begin{array}{c}\hfill {\overline{𝐍}}^T\hfill \\ \hfill {\overline{𝐌}}^T\hfill \end{array}\right]𝑑\zeta 𝑑\gamma$	(26)
$F_m^e$	$=RL\pi \mathrm{\Delta }m{\displaystyle {\int }_1^{-1}}{\displaystyle {\int }_1^{-1}}{\left[{𝐁}_u(\zeta ,\gamma )\right]}^T\left[\begin{array}{c}\hfill {\overline{𝐍}}^m\hfill \\ \hfill {\overline{𝐌}}^m\hfill \end{array}\right]𝑑\zeta 𝑑\gamma$

in which

${\overline{𝐍}}^T$	$={\displaystyle \sum _{k=1}^{N_l}}{\displaystyle {\int }_{R-h_k}^{R-h_{k+1}}}{\overline{𝐐}}_k\left[\begin{array}{c}\hfill {\alpha }_z\hfill \\ \hfill {\alpha }_{\alpha }\hfill \\ \hfill 2{\alpha }_{\alpha z}\hfill \end{array}\right]𝑑r,$	(27)
${\overline{𝐌}}^T$	$={\displaystyle \sum _{k=1}^{N_l}}{\displaystyle {\int }_{h_{k-1}}^{R-h_k}}{\overline{𝐐}}_k\left[\begin{array}{c}\hfill {\alpha }_z\hfill \\ \hfill {\alpha }_{\alpha }\hfill \\ \hfill 2{\alpha }_{\alpha z}\hfill \end{array}\right]\left(r-R\right)𝑑r$
${\overline{𝐍}}^m$	$={\displaystyle \sum _{k=1}^{N_l}}{\displaystyle {\int }_{R-h_k}^{R-h_{k+1}}}{\overline{𝐐}}_k\left[\begin{array}{c}\hfill {\beta }_z\hfill \\ \hfill {\beta }_{\alpha }\hfill \\ \hfill 2{\beta }_{\alpha z}\hfill \end{array}\right]𝑑r,$
${\overline{𝐌}}^m$	$={\displaystyle \sum _{k=1}^{N_l}}{\displaystyle {\int }_{h_{k-1}}^{R-h_k}}{\overline{𝐐}}_k\left[\begin{array}{c}\hfill {\beta }_z\hfill \\ \hfill {\beta }_{\alpha }\hfill \\ \hfill 2{\beta }_{\alpha z}\hfill \end{array}\right]\left(r-R\right)𝑑r$

Figure 2 Layers in the laminate.

Figure 3 Schematic of the cylindrical composite shell.

4 A Case Study

In this case study, a thin cylindrical vessel made of a laminated composite material with a lay-up configuration of ${[0/\pm 45/90]}_s$ is considered, as shown in (Figure 3). Each layer is a unidirectional carbon/epoxy (T300/934) material with a thickness of 0.5 mm; the corresponding mechanical properties are given in Table 1. The cylindrical shell is assumed to have a inner radius of 0.2 m and a length of 0.8 m, and both ends are assumed to be closed with rigid caps. The internal pressure of 250 psi is applied inside the cylinder while the inside and outside temperatures are ${120}^{\circ }C$ and ${23}^{\circ }C$, respectively. The caps rigidly hold both ends of the cylinder, therefore, except the axial translation at one end, the other translational displacements are considered fixed at both boundaries, while the rotational degrees-of-freedom are allowed. The cylinder with the applied conditions is depicted in Figure 3. In this problem, three types of force vectors are considered:

a) the radial force vector due to the internal pressure,

$${𝐟}_1^e=\int {\int }_A{𝐍}_u^T{𝐏}_i𝑑A$$

(28)

in which ${𝐏}_i$ is the vector of internal pressure, namely,

$${𝐏}_i={\left[\begin{array}{cccccc}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]}^T\mathit{\hspace{1em}}(MPa)$$

(29)

b) The axial force vector due to the internal pressure whose total magnitude equals to

$$F_s=P_i\left(\pi R_i^2\right)$$

(30)

Because there are four nodes at each end, the above magnitude is divided by four and applied on each node. In this problem, one end is free to move axially, hence, the axial forces are only applied to the nodes at the fixed end (as it is shown in Figure 3). Therefore, the axial force vector is defined as

$${𝐟}_2^n={\left[\begin{array}{cccccccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{F_s}{4}\hfill & \hfill \frac{F_s}{4}\hfill & \hfill \frac{F_s}{4}\hfill & \hfill \frac{F_s}{4}\hfill \end{array}\right]}^T$$

(31)

It is noteworthy to mention that this axial force vector only exists for the $n^{th}$ element (the last one), and for the remaining elements it equals to zero.

c) The last force vector is related to the thermal load, namely,

$${𝐟}_3^e=RL\pi \mathrm{\Delta }T{\int }_{-1}^1{\int }_{-1}^1{𝐁}_u^T(\zeta ,\gamma )\left[\begin{array}{c}\hfill {\overline{𝐍}}^T\hfill \\ \hfill {\overline{𝐌}}^T\hfill \end{array}\right]𝑑\zeta 𝑑\gamma$$

(32)

The total force vector that is applied on the structure equals

$$𝐟={𝐟}_1+{𝐟}_2+{𝐟}_3$$

(33)

where ${𝐟}_1$, ${𝐟}_2$, and ${𝐟}_3$ are the assembled radial, axial, and thermal force vectors, respectively.

Table 1 Mechanical properties of the unidirectional carbon/epoxy (T300/934) [2]

$E_1(GPa)$	$E_2(GPa)$	$G_{12},G_{13}(GPa)$	$G_{23}(GPa)$
$141.6$	$10.7$	$3.88$	$3.58$
${\nu }_{12},{\nu }_{13}$	${\nu }_{23}$	${\alpha }_1(1{/}^{o}C)$	${\alpha }_2(1{/}^{o}C)$
$0.268$	$0.495$	$-0.006\times {10}^{-6}$	$30.04\times {10}^{-6}$

Figure 4 (a) Uniform superelements, (b) Non-uniform superelements.

Figure 5 (a)–(d) The radial displacements obtained by 10 and 20 uniform and non-uniform elements.

Figure 6 (a)–(d) The axial displacements obtained by 10 and 20 uniform and non-uniform elements.

4.1 Results

In the first attempt, the cylindrical shell is meshed longitudinally by 10 uniform superelements. Figure 5 shows the radial displacements obtained by the superelements and an analytical method which was presented in [27]. As it is apparent, at places where the bending is large, the superelement result deviated from the analytical solution. The bending boundary layer of a shell structure $L_b$ can be calculated via the following formula [27]

$$L_b=4\sqrt{Rt{(D_{11}/D_{22})}^{1/2}}$$

(34)

where in the current case $L_b=133.4mm$. This means that the maximum bending deformation occurs approximately 134 mm from each end. With 10 superelements, each element has a length of 80 mm, and hence, two elements cannot accurately report the aggressive change of bending deformation. Thus, the number of elements should be elevated via a) decreasing the size of elements uniformly, b) using non-uniform superelements. In this case, 10 uniform and non-uniform superelements are incorporated. Figures 5(a) and 5(b) show the meshing with 10 uniform and non-uniform superelements, respectively. In Figures 5(a) and 5(b), the radial displacements which were calculated by 10 uniform and non-uniform elements are compared with the analytical results which were presented in [27]. The corresponding axial displacements are also depicted in Figures 6(a) and 6(b). In order to conduct a mesh independence analysis, the cylindrical shell is meshed with 20 and 30 uniform and non-uniform elements. The radial and axial displacements which are calculated by 20 uniform and non-uniform elements are presented in Figures 5(c), 5(d), 6(c), and 6(d), respectively. As it is apparent, the displacements converged rapidly as the number of elements is increased. Moreover, by comparing the results depicted in Figures 5 and 6, it is seen that proper accuracy can be obtained if the superelements are properly concentrated at places with abrupt changes. The maximum relative error (MRE) for the displacement results calculated with 10, 20, and 30 uniform and non-uniform elements are presented in Table 2. The last column of Table 2 also shows the computation run time corresponding to each number of elements. As it is expected, the maximum error occurs at places with sudden bending deformation changes. According to Table 2, non-uniform meshing leads to better results.

Table 2 Maximum Relative error (MRE) for displacement results obtained by different numbers of elements and the corresponding run time (CPU: 1.6 GHz Core i7, RAM: 8Gb.)

Number of Elements	MRE $(u_r\%)$	MRE $(u_z\%)$	Run time (s)
$10$	$5.76$	$12.86$	$0.55$
$10$ N	$2.41$	$6.41$	$0.43$
$20$	$0.23$	$1.87$	$1.28$
$20$ N	$0.09$	$0.64$	$1.28$
$30$	$0.20$	$0.38$	$3.27$
$30$ N	$0.07$	$0.12$	$3.27$

4.1.1 Stress Results

After calculation of the nodal displacements, the stress results can be derived at each element. Here, the stresses are calculated via the stress recovery method which is extensively explained in [31]. For the case study at hand, the stresses are calculated at each layer and at two different sections A and B, which are shown in Figure 7. The stresses are derived from the displacements. Hence, in order to keep accurate results, it is needed to use smaller elements, especially at places with drastic changes. Figure 8 shows the axial, tangential, and shear stresses at the plane A for 10 uniform and non-uniform and 20 non-uniform elements. The corresponding error for 10, 20, and 30 uniform and non-uniform elements are presented in Table 3. However, as it is depicted in Figure 9, at the plane B even with 10 elements the result is accurate enough. This is because the stress is constant over the central part of the composite cylinder.

Figure 7 Sections A and B.

Figure 8 (a), (b) and (c) The tangential stresses, (c), (d), and (e) the axial stresses, (e), (f), and (g) the shear stresses calculated with 10 uniform, 10 non-uniform and 20 non-uniform elements.

Figure 9 (a)–(c) The tangential, axial and shear stresses at section B calculated with 10 uniform elements.

The maximum relative error (MRE) of the stress results at the section A that are obtained by 10, 20, and 10 non-uniform elements are presented in Table 3. As the stress results are recovered from the nodal displacements, the errors are elevated, which is apparent from Table 1. However, in stress analysis the non-uniform meshing still leads to better results.

Table 3 Maximum Relative error (MRE) for stress results obtained by different number of elements.

Number of Elements	MRE $({\sigma }_{\alpha \alpha }\%)$	MRE $({\sigma }_{zz}\%)$	MRE $({\tau }_{zt}\%)$
$10$	$34.13$	$84.1$	$22.94$
$10$ N	$20.51$	$39.42$	$14.25$
$20$	$4.47$	$10.32$	$2.85$
$20$ N	$0.69$	$6.89$	$0.48$
$30$	$1.96$	$5.26$	$1.25$
$30$ N	$0.35$	$2.78$	$0.29$

Figure 10 The cylinder which is meshed by shell elements.

Figure 11 (a) The radial and (b) the axial displacement plots of the cylinder obtained by a commercial software package.

5 Comparison with the results of shell elements

In this section, the results of the mentioned case study, which were obtained by means of uniform superelements, will be compared with the ones which are calculated by shell elements and via a commercial software package. Figure 10 shows the cylinder understudy which is meshed by shell elements. The boundary conditions, the applied internal pressure, temperature and the caps effect are defined same as to what was mentioned in the case study explanation. The obtained radial and axial displacement plots are depicted in Figures 11(a) and 11(b), respectively. In this regard, Figure 12 compares the radial and axial displacements which are obtained by superelements, shell elements and an analytical method. In Figure 12(a), the radial displacements which are estimated by 400 shell elements are compared with the ones resulted from only 10 superelements. Apparently, the results with 400 shell elements are compatible with the ones obtained by 10 superelements, however, both are deviated from the analytical solution close to the boundaries where the cylinder deform significantly. This deviation is resolved when the number of shell elements are decreased to 1632 and the superelements to 20, as it is depicted in Figure 12(b). Likewise, the corresponding axial displacements are also compared in Figures 12(c) and 12(d). In order to quantitatively verify the results of superelements, the maximum value of radial and axial displacements which are obtained by 1632 shell elements and 20 superelements are compared with the analytical ones in Table 4. After the comparison of displacement plots, now, the tangential, axial and shear stresses at sections A and B of cylinder (Figure 7) which are obtained by all three approaches are compared and depicted in Figures 13(a)–13(b). As it is apparent, because the section A is closer to one of the boundaries, the finite element results are deviated from the analytical ones; this can be resolved by adding more elements in that region as previously discussed. Table 5 shows the relative error of the stresses at section A in the layers where the maximum deviations happen between the ones obtained by the superelements and the analytical approach. As it can be seen, few cylindrical superelements can properly estimate the displacements and stresses in the model, comparing with the ones obtained by the analytical method and many shell elements.

Figure 12 (a) The radial and (b) the axial displacement plots of the cylinder obtained by a commercial software package.

Table 4 The maximum radial and axial displacements obtained by superelements and shell elements, and the corresponding relative errors with the analytical results.

	Analytical	Superelements	Error%	Shell Elements	Error%
Radial Disp.	0.336	0.335	0.3	0.344	2.3
Axial Disp.	-0.421	-0.429	1.9	-0.433	2.8

Figure 13 The tangential, axial and shear stresses calculated with 20 uniform superelements, 1632 shell elements and the analytical method.

Table 5 The tangential, axial and shear stresses which are obtained at section A by 20 superelements and 1632 shell elements, and the corresponding relative errors with the analytical results.

Stress (MPa)	Analytical	Superelements	Error%	Shell Elements	Error%
Tangential (Layer 1)	-12.69	-12.86	1.34	-12.38	2.44
Axial (Layer 8)	69.94	71.23	1.84	72.67	3.9
Shear (Layer 2)	-82.93	-82.45	0.24	-82.65	0.33

6 Conclusions

In this paper, a cylindrical thin superelement is formulated based on the classical theory of shell and the trigonometric shape functions. The element targets the thermo-mechanical analysis of cylindrical thin structures made of composite materials; however, its application can be extended to other analyses such as coupled field or nonlinear analyses in which a large number of degrees of freedom, and accordingly, large number of equations exist. As a case study, a thermo-mechanical analysis of a thin cylindrical composite vessel, under an internal pressure and a temperature gradient, was conducted. In this example, a mesh independence analysis was performed by incorporating 10, 15, and 20 uniform elements. It is shown that as the number of superelements is increased the deformations are rapidly converged to the ones which were obtained in the literature analytically. Moreover, it is concluded that the same number of non-uniform elements can lead to more accurate results. After the calculation of displacements, the stress results are also derived, and accordingly, compared with the analytical ones. Finally, for the same case study, the displacement and stress plots which were obtained by superelements were compared with the ones calculated by shell elements. The comparison revealed that the incorporation of the superelements in the mechanical analysis of cylindrical structures is accurate and efficient.

Appendix
The Shape Functions

The shape functions of the cylindrical thin superelement corresponding to the $C^1$-continuous approximation function of radial displacement are obtained as

$N_{i,j}$	$=f_i(\zeta )g_j(\gamma ),N_{i,j}^{\prime }=F_i(\zeta )g_j(\gamma )$
$N_{i,j}^{\prime \prime }$	$=f_i(\zeta )G_j(\gamma ),N_{i,j}^{\prime \prime \prime }=F_i(\zeta )G_j(\gamma )$
	$i=1,2,j=1,\mathrm{\dots },4$

in which

$f_1(\zeta )$	$={\displaystyle \frac{1}{4}}\left({\zeta }^3-3\zeta +2\right),f_2(\zeta )=-{\displaystyle \frac{1}{4}}\left({\zeta }^3-3\zeta -2\right)$
$F_1(\zeta )$	$={\displaystyle \frac{1}{4}}\left({\zeta }^3-{\zeta }^2-\zeta +1\right),F_2(\zeta )={\displaystyle \frac{1}{4}}\left({\zeta }^3+{\zeta }^2-\zeta -1\right)$
$g_1(\gamma )$	$={\displaystyle \frac{1}{8}}\left(-3\cos \pi \gamma +2\cos 2\pi \gamma -\cos 3\pi \gamma +2\right)$
$g_2(\gamma )$	$={\displaystyle \frac{1}{8}}\left(-3\sin \pi \gamma -2\cos 2\pi \gamma +\sin 3\pi \gamma +2\right)$
$g_3(\gamma )$	$={\displaystyle \frac{1}{8}}\left(3\cos \pi \gamma +2\cos 2\pi \gamma -\cos 3\pi \gamma +2\right)$
$g_4(\gamma )$	$={\displaystyle \frac{1}{8}}\left(3\sin \pi \gamma -2\cos 2\pi \gamma -\sin 3\pi \gamma +2\right)$
$G_1(\gamma )$	$={\displaystyle \frac{1}{8\pi }}\left(-\sin \pi \gamma +\sin 2\pi \gamma -\sin 3\pi \gamma +0.5\sin 4\pi \gamma \right)$
$G_2(\gamma )$	$={\displaystyle \frac{1}{8\pi }}\left(-\cos \pi \gamma -\sin 2\pi \gamma -\cos 3\pi \gamma +0.5\sin 4\pi \gamma \right)$
$G_3(\gamma )$	$={\displaystyle \frac{1}{8\pi }}\left(\sin \pi \gamma +\sin 2\pi \gamma +\sin 3\pi \gamma +0.5\sin 4\pi \gamma \right)$
$G_4(\gamma )$	$={\displaystyle \frac{1}{8\pi }}\left(-\cos \pi \gamma -\sin 2\pi \gamma +\cos 3\pi \gamma +0.5\sin 4\pi \gamma \right)$

Also, the shape functions corresponding to the $C^0$-continuous approximation functions of tangential and axial displacements are presented as

$M_{i,j}$	$=h_i(\zeta )I_j(\gamma ),i=1,2\mathit{\hspace{1em}}j=1,\mathrm{\dots },4$
$h_1(\zeta )$	$={\displaystyle \frac{1}{2}}\left(1-\zeta \right),h_2(\zeta )={\displaystyle \frac{1}{2}}\left(1+\zeta \right)$

where,

$I_1(\gamma )$	$={\displaystyle \frac{1}{4}}\left(1-2\cos \pi \gamma +\cos 2\pi \gamma \right),I_2(\gamma )={\displaystyle \frac{1}{4}}\left(1-2\sin \pi \gamma -\cos 2\pi \gamma \right)$
$I_3(\gamma )$	$={\displaystyle \frac{1}{4}}\left(1+2\cos \pi \gamma +\cos 2\pi \gamma \right),I_2(\gamma )={\displaystyle \frac{1}{4}}\left(1+2\sin \pi \gamma -\cos 2\pi \gamma \right)$

References

[1] M. Xia, H. Takayanagi, K. Kemmochi, ‘Analysis of multi-layered filament-wound composite pipes under internal pressure’, Composite Structures 53 (2001) 483–491.

[2] M. Xia, H. Takayanagi, K. Kemmochi, ‘Analysis of lament-wound ber-reinforced sandwich pipe under combined internal pressure and thermo mechanical loading’, Composite Structures 51 (2001) 473–283.

[3] H. Bakaiyan, H. Hosseini, E. Ameri, ‘Analysis of multi-layered filament-wound composite pipes under combined internal pressure and thermomechanical loading with thermal variations’, Composite Structures 88 (2009) 532–541.

[4] P. Xu, J. Y. Zheng, ‘Finite element analysis of burst pressure of composite hydrogen storage vessels’, Composite Structures 30 (2009) 2295–2301.

[5] L. Parnas, N. Katirci, ‘Design of fiber-reinforced composite pressure vessels under various loading conditions’, Composite Structures 58 (2002) 83–95.

[6] L. A. L. Martins, F. L. Bastian, T. O. Netto, ‘Structural and functional failure pressure of filament wound composite tubes’, Materials and Design 36 (2012) 779–787.

[7] O. Sayman, ‘Analysis of multi-layered composite cylinders under hygrothermal loading’, Composites: Part A 36 (2005) 923–933.

[8] Q. Zhang, Z. Wang, C. Tang, D. Hu, P. Liu, L. Xia, ‘Analytical solution of the thermo-mechanical streeses in a multilayered composite pressure vessel considering the influence of the closed ends’, International journal of pressure vessels and piping 98 (2012) 102–110.

[9] K. Vedeld, H. Sollund, ‘Stresses in Heated pressurized multi-layer cylinders in generalized plane strain conditions’, International journal of pressure vessels and piping 120-121 (2014) 27–35.

[10] M. Jabbari, S. Sohrabpour, M. Eslami, ‘Mechanical and thermal stresses in a functionally graded hollow cylinder due to a radially symmetric loads’, International journal of pressure vessels and piping 79 (2002) 493–497.

[11] Z. Shao, ‘Mechanical and thermal stresses of a functionally graded circular hollow cylinder with finite length’, International journal of pressure vessels and piping 82 (2005) 155–163.

[12] T. S. Koko, ‘Super finite elements for nonlinear static and dynamic analysis of stiffened plate structures’, Ph.D. Thesis.

[13] T. Koko, M. D. Olson, ‘Vibration analysis of stiffened plates by super elements’, J. Sound Vib. 158 (1991) 149–167.

[14] J. Jiang, M. D. Olson, ‘Vibration Analysis of Orthogonally Stiffened Cylindrical Shells Using Super Finite Elements’, J. Sound Vib. 173 (1991) 73–83.

[15] J. Jiang, M. D. Olson, ‘Nonlinear analysis of orthogonally stiffened cylindrical shells by a super element approach’, Finite Elem. Anal. Des. 18 (1994) 99–110.

[16] M. T. Ahmadian, M. Zanganeh, ‘Nonlinear analysis of orthogonally stiffened cylindrical shells by a super element approach’, Comput. Methods Appl. Mech. Eng. 191 (2002) 2069–2075.

[17] M. T. Ahmadian, M. Zanganeh, ‘Application of super elements to free vibration analysis of laminated stiffened plates’, J. Sound Vib. 259 (2003) 1243–1252.

[18] F. Ju, Y. S. Choo, ‘Super element approach to cable passing through multiple pulleys’, Int. J. Solids Struct. 42 (2005) 3533–3547.

[19] M. T. Ahmadian, M. Bonakdar, ‘A new cylindrical element formulation and its application to structural analysis of laminated hollow cylinders’, Finite Element in Analysis and Design 44 (2008) 617–630.

[20] A. Taghvaeipour, M. Bonakdar, M. T. Ahamadian, ‘Application of a new cylindrical element formulation in finite element structural analysis of FGM hollow cylinders’, Finite Element in Analysis and Design 50 (2012) 1–7.

[21] M. Nasiri Sarvi, M. T. Ahmadian, ‘Design and implementation of a new spherical super element in structural analysis’, Finite Element in Analysis and Design 218 (14) (2012) 7546–7561.

[22] M. Nasiri Sarvi, M. T. Ahmadian, Deformation, stress and natural frequency analysis of the fullerene by finite super element method, International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, August 28-31, Washington, DC, USA, 2011.

[23] M. T. Ahmadian, M. R. Movahhedy, M. M. Rezaei, ‘Design and application of a new tapered superelement for analysis of revolving geometries’, Finite Element in Analysis and Design 218 (47) (2011) 1242–1252.

[24] R. Pourhamid, M. T. Ahmadian, H. Mahdavy Moghaddam, A. R. Mohammadzadeh, ‘Mechanical analysis of a functionally graded cylinder-piston under internal pressure due to a combustion engine using a cylindrical super element and considering thermal loading’, Scientia Iranica 22 (2) (2015) 493–5–3.

[25] A. E. H. Love, A Treatise On The Mathematical Theory of Elasticity, Dover, New York, 1944.

[26] S. Timoshenko, S. Woinowski-Krieger, Theory of plate and shells, McGraw-Hill, New York, 1959.

[27] E. Venstel, T. Krauthammer, Thin plates and shells, Marcel Dekker, New York, 2001.

[28] M. R. Eslami, Finite Elements Methods in Mechanics, Springer International Publishing, Switzerland, 2014.

[29] S. S. Rao, The Finite Element Method in Engineering, Butterworth-Heinemann, United Kingdom, 2011.

[30] J. R. Vinson, R. L. Sierakowski, The Behavior of Structures Composed of Composite Materials, Kluwer Academic Publishers, New York, 2002.

[31] C. Felippa, Introduction to Finite Element Methods (ASEN 5007), http://www.colorado.edu/engineering/CAS/courses.d/IFEM.d.

Biographies

A. Jafarzadeh received his B.Sc. degree in mechanical engineering from University of Shahroud, Iran, in 2014 and the M.Sc. degree from the Amirkabir University of Technology (AUT), Iran in 2017. His research interests include finite element methods, numerical methods in engineering, creep and fatigue Life assessment in the hot section parts of gas turbines, composite structures design, designing machines and modulus for assembly and production process.

A. Taghvaiepour is an assistant professor at Amirkabir University of Technology (AUT), Iran. He received his B.Sc. from AUT in 2006, his M.Sc. from Sharif University of Technology (SUT) in 2008, and his Ph.D. form McGill University in 2012. His research interests are Structural and Dynamic analysis of Mechanical and Robotic Systems and Finite Element Analysis.

M. R. Eslami is an emeritus professor of mechanical engineering at Amirkabir University of Technology (AUT), Iran. He received his B.Sc. degree in mechanical engineering from Tehran Polytechnic, Iran, in 1968, the M.Sc. and Ph.D. degrees from the Louisiana State University, Baton Rouge La., USA, in 1970 and 1973, respectively. He is an ASME fellow. His main research fields are thermoelasticity and finite element methods (FEM).

European Journal of Computational Mechanics, Vol. 29_2-3, 173–198.
doi: 10.13052/ejcm2642-2085.29231
© 2020 River Publishers

Abstract

1 Introduction

2 Element Definition

2.1 Geometry and Coordinate System of Element

2.2 Element Degrees of Freedom and Shape Functions

3 The FEM Formulations

3.1 Kinematic Relations

3.2 Stress-strain relations

3.3 Stiffness matrix and force vectors

4 A Case Study

4.1 Results

4.1.1 Stress Results

5 Comparison with the results of shell elements

6 Conclusions

Appendix
The Shape Functions

References

Biographies