**1. Introduction**

Functionally graded materials (FGM) are usually made from a mixture of metals and ceramics, whose material properties vary smoothly and continuously from one surface to the other of the structure according to volume fraction power–law distribution. It is well known that the ceramics are capable of resisting high temperature, while the metals provide structural strength and fracture toughness. They are therefore suitable to apply for aerospace structures, nuclear plants, and other applications. With the advantageous features of the FGM in many practical applications, the problem of static and free vibration behaviors of FGM shell structures are attractive to many researchers over the world. Woo and Meguid [1] used an analytical solution based on the von Karman theory to investigate nonlinear respond of FGM plates and shallow shells. Matsunaga [2] carried out the power series expansion of displacement component approach, which relied on higher-order shear deformation theory (HSDT) to analyze free vibration and buckling of FGM shells. Nguyen et al. [3] proposed an analytical solution using Reddy's HSDT to solve nonlinear dynamic and free vibration of piezoelectric FGM double curved shallow shells subjected to electrical, thermal, mechanical, and damping loads. Dao et al. [4] presented nonlinear vibration of sti ffened functionally graded double curved shells on an elastic foundation using the first order shear deformation theory (FSDT) and stress function. Due to

the complication of mathematics, it is generally difficult to use analytical methods for all problems. Thus, numerical methods have been devised to study the behavior of FGM structural components. Among these numerical approaches, the finite element method has become the most powerful, reliable, and simply tool to solve FGM shells. Arciniega and Reddy [5] presented a finite element formulation for nonlinear analysis of FGM shell based on the FSDT, consisting of seven parameters. Pradyumna and Bandyopadhyay [6] used the shell element, consisting of nine degrees of freedom per node, to investigate free vibration of FGM shells. Kordkheili and Naghdabadi [7] proposed a finite element model for geometrically nonlinear thermos-elastic analysis of FGM plates and shells.

In addition, in the recent trend of development of numerical methods, the edge-based smoothed finite element method (ES-FEM) combined with the mixed interpolation of tensorial components using triangular element (MITC3), named ES-MITC3, has been proposed to investigate plate and shell structures. For instance, Chau-Dinh et al. [8] proposed an ES-MITC3 to analyze static and free vibration of plates. Nguyen et al. [9] developed the ES-MITC3 for static and vibration analysis of isotropic and functionally graded plates. Pham et al. [10] examined the static and free vibration of composite shells using ES-MITC3 shell element. Pham et al. [11] used ES-MITC3 shell element to study geometrically nonlinear analysis of FGM shells based on FSDT. Pham-Tien et al. [12] investigated the dynamic response of composite shells based on the FSDT and ES-MITC3 element. Hoang-Nam Nguyen et al. [13] used FSDT to investigate dynamic composite shell with shear connectors. For shell structures, especially the shell with two curvatures, the employing of quadrilateral elements will not accurately describe the model due to the distorting during the bending process. In this case, the using of triangle elements is suitable because they can rotate freely around their three edges. However, the using of these elements can meet the shear locking phenomenon; thus, we propose the new method, in which the triangle element is combined with an edge-based smoothed finite element method (ES-MITC3) to analyze the shell structures. Its accuracy in comparison with other methods is shown in the numerical exploration.

This paper now further extends the ES-MITC3 method for free vibration analysis of functionally graded shell structures. The material properties of functionally graded shells are assumed to vary continuously and smoothly through the thickness based on a simple power–law of the volume fractions exponents. The formulation is based on the FSDT and flat shell theory due to the simplicity and computational efficiency. The accuracy and reliability of the present method are verified by comparing with those of others available numerical results.

### **2. Theoretical Formulation**

### *2.1. Functionally Grade Material*

FGM is formed by a mixture of ceramic and metal, as shown in Figure 1. The material properties change continuously from a surface to the other surface according to a power–law of volume fraction

$$P(z) = (P\_c - P\_m)V\_c + P\_m \tag{1}$$

$$V\_{\mathbb{C}}(z) = \left(\frac{1}{2} + \frac{z}{h}\right)^{n} \tag{2}$$

where *<sup>P</sup>*(*z*) represents the effective material properties: Young's modulus *E*, density ρ and Poisson ratio *v*; *Pc* and *Pm* denote the properties of the ceramic and metal, respectively; *Vc* is the volumefractions of the ceramic; *h* the thickness of structure; *n* ≥ 0 the volumefraction exponent; *z* ∈ [−*h*/2, *h*/2] is the thickness coordinate of the structure. Figure 2 illustrates the variation of the volume fraction of ceramic and metal through the thickness via the volumefraction exponent *n*.

**Figure 2.** Variation of the volume fraction versus the non-dimensional thickness.

1RQGLPHQVLRQDOWKLFNQHVV]K

 

### *2.2. The FGM Shell Model*

 

Consider an FGM shell element subjected to both in-plane forces and bending moments as shown in Figure 3. The middle (neutral) surface of the shell is chosen as the reference plane that occupies a domain Ω ∈ 3. Let *u*0, *v*0, and *w*0 be the displacements of the middle plane in the *x*, *y*, and *z* directions; β*<sup>x</sup>*, β*<sup>y</sup>*, and β*z* be the rotations of the middle surface of the shell around the *y*-axis, *x*-axis, and *z*-axis, respectively, as indicated in Figure 3. The unknown vector of an FGM shell including six independent variables at any point in the problem domain can be written as

$$\mathbf{u} = \begin{bmatrix} u\_0 & v\_0 & w\_0 & \beta\_x & \beta\_y & \beta\_z \end{bmatrix} \tag{3}$$

The linear strain–displacement relationship can be given as

$$\mathfrak{e} = \left\{ \begin{array}{c} \varepsilon\_x \\ \varepsilon\_y \\ \gamma\_{xy} \end{array} \right\} = \mathfrak{e}\_m + z\mathfrak{e} \tag{4}$$

$$\mathfrak{c}\_{m} = \begin{Bmatrix} \frac{\partial u\_{0}}{\partial x} \\ \frac{\partial v\_{0}}{\partial y} \\ \frac{\partial u\_{0}}{\partial y} + \frac{\partial v\_{0}}{\partial x} \end{Bmatrix}, \mathfrak{K} = \left\{ \begin{array}{c} \frac{\partial \mathfrak{k}\_{x}}{\partial x} \\ \frac{\partial \mathfrak{k}\_{y}}{\partial y} \\ \frac{\partial \mathfrak{k}\_{x}}{\partial y} + \frac{\partial \mathfrak{k}\_{y}}{\partial x} \end{array} \right\}; \tag{5}$$

$$\mathbf{Y} = \begin{Bmatrix} \boldsymbol{\mathcal{V}xx} \\ \boldsymbol{\mathcal{V}yz} \end{Bmatrix} = \begin{Bmatrix} \frac{\partial \boldsymbol{w}\_0}{\partial \boldsymbol{x}} + \boldsymbol{\beta}\_\mathbf{x} \\ \frac{\partial \boldsymbol{w}\_0}{\partial \boldsymbol{y}} + \boldsymbol{\beta}\_\mathbf{y} \end{Bmatrix} \tag{6}$$

From Hooke's law, the constitutive relations of FGM shells are expressed as:

⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎩

$$
\begin{pmatrix}
\sigma\_{xx} \\
\sigma\_{yy} \\
\tau\_{xy} \\
\tau\_{xz} \\
\tau\_{yz}
\end{pmatrix} = \begin{bmatrix}
Q\_{11} & Q\_{12} & 0 & 0 & 0 \\
Q\_{21} & Q\_{22} & 0 & 0 & 0 \\
0 & 0 & Q\_{66} & 0 & 0 \\
0 & 0 & 0 & Q\_{55} & 0 \\
0 & 0 & 0 & 0 & Q\_{44}
\end{bmatrix} \begin{pmatrix}
\varepsilon\_{xx} \\
\varepsilon\_{yy} \\
\gamma\_{xy} \\
\gamma\_{xz} \\
\gamma\_{yz}
\end{pmatrix} \tag{7}
$$

where

$$Q\_{11}(z) = \frac{E(z)}{1 - \nu(z)^2},\ Q\_{12}(z) = \nu(z)Q\_{11}(z),\ Q\_{22}(z) = Q\_{11}(z),\ Q\_{44}(z) = Q\_{55}(z) = Q\_{66}(z) = \frac{E(z)}{2(1 + \nu(z))}\tag{8}$$

A weak form of the free vibration analysis for FGM shells can be briefly given as:

$$
\int\_{\Omega} \delta \overline{\boldsymbol{\varepsilon}}^T \mathbf{D} \overline{\boldsymbol{\varepsilon}} d\Omega + \int\_{\Omega} \delta \mathbf{y}^T \mathbf{C} \mathbf{y} d\Omega = \int\_{\Omega} \delta \mathbf{u}^T \mathbf{m} \ddot{\mathbf{u}} d\Omega \tag{9}
$$

where

$$
\overline{\boldsymbol{\varepsilon}} = \begin{bmatrix} \boldsymbol{\varepsilon}\_{\rm{W}} \\ \boldsymbol{\kappa} \end{bmatrix}, \; \mathbf{D} = \begin{bmatrix} \mathbf{D}^{\rm{W}} & \mathbf{B} \\ \mathbf{B} & \mathbf{D}^{\rm{b}} \end{bmatrix} \tag{10}
$$

**Figure 3.** Three-node triangular element.

In which

$$\mathbf{D}^{\rm m} = \begin{bmatrix} \mathbf{A}\_{11} & \mathbf{A}\_{12} & \mathbf{A}\_{16} \\ \mathbf{A}\_{12} & \mathbf{A}\_{22} & \mathbf{A}\_{26} \\ \mathbf{A}\_{16} & \mathbf{A}\_{26} & \mathbf{A}\_{66} \end{bmatrix}, \mathbf{B} = \begin{bmatrix} \mathbf{B}\_{11} & \mathbf{B}\_{12} & \mathbf{B}\_{16} \\ \mathbf{B}\_{12} & \mathbf{B}\_{22} & \mathbf{B}\_{26} \\ \mathbf{B}\_{16} & \mathbf{B}\_{26} & \mathbf{B}\_{66} \end{bmatrix}, \mathbf{D}^{b} = \begin{bmatrix} \mathbf{D}\_{11} & \mathbf{D}\_{12} & \mathbf{D}\_{16} \\ \mathbf{D}\_{12} & \mathbf{D}\_{22} & \mathbf{D}\_{26} \\ \mathbf{D}\_{16} & \mathbf{D}\_{26} & \mathbf{D}\_{66} \end{bmatrix} \tag{11}$$

In Equation (11) <sup>A</sup>*ij*, <sup>B</sup>*ij*, <sup>D</sup>*ij*, and C*ij*are given by:

$$\left(\mathbf{A}\_{i\circ}, \mathbf{B}\_{i\circ}, \mathbf{D}\_{i\circ}\right) = \int\_{-\hbar/2}^{\hbar/2} \mathbf{Q}\_{i\circ}(1, z, z^2) dz, \text{ i.e.} \ j = 1, 2, 6 \tag{12}$$

$$\mathbf{C}\_{ij} = k \int\_{-h/2}^{h/2} Q\_{ij} dz\_i \text{ i, j} = 4,5 \tag{13}$$

where *k* = 5/6 is transverse shear correction coefficient and **m** is the mass matrix containing ρ calculated as

$$\mathbf{m} = \begin{bmatrix} I\_0 & 0 & 0 & I\_1 & 0 & 0 \\ 0 & I\_0 & 0 & 0 & I\_1 & 0 \\ 0 & 0 & I\_0 & 0 & 0 & 0 \\ I\_1 & 0 & 0 & I\_2 & 0 & 0 \\ 0 & I\_1 & 0 & 0 & I\_2 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}; \begin{pmatrix} I\_0, I\_1, I\_2 \end{pmatrix} = \int\_{-h/2}^{h/2} \rho(1, z, z^2) dz \tag{14}$$

### *2.3. Finite Element Formulation for Shell Analysis*

Now, we discretize the bounded domain Ω into *n<sup>e</sup>* finite three-node triangular elements with *n<sup>n</sup>* nodes such that Ω ≈ *n<sup>e</sup> <sup>e</sup>*=1 Ω*e* and Ω*i* ∩ Ω*j* = ∅, *i* - *j*. The displacement field **u***<sup>e</sup>* = *u<sup>e</sup> v<sup>e</sup> w<sup>e</sup>* β*ex* β*ey* β*ez T* of the finite element solution can be expressed as

$$\mathbf{u}^{\varepsilon} = \sum\_{i=1}^{n^{\varepsilon}} \begin{bmatrix} N\_i(\mathbf{x}) & 0 & 0 & 0 & 0 & 0 \\ 0 & N\_i(\mathbf{x}) & 0 & 0 & 0 & 0 \\ 0 & 0 & N\_i(\mathbf{x}) & 0 & 0 & 0 \\ 0 & 0 & 0 & N\_i(\mathbf{x}) & 0 & 0 \\ 0 & 0 & 0 & 0 & N\_i(\mathbf{x}) & 0 \\ 0 & 0 & 0 & 0 & 0 & N\_i(\mathbf{x}) \end{bmatrix} \mathbf{d}\_i^{\varepsilon} = \sum\_{i=1}^{n^{\varepsilon}} \mathbf{N}\_i \mathbf{d}\_i^{\varepsilon} \tag{15}$$

where **d***ei* = *ue*0*i ve*0*i <sup>w</sup><sup>e</sup>*0*i* β*exi* β*eyi* β*ezi T* is the nodal displacement at the *i*th node; *Ni*(**x**) is shape function for the *i*th node.

The approximation of the membrane, the bending and the shear strains of the triangular element can be written in matrix forms as follows

$$\mathbf{c}^{\varepsilon} = \begin{bmatrix} \mathbf{B}^{\varepsilon}\_{m1} & \mathbf{B}^{\varepsilon}\_{m2} & \mathbf{B}^{\varepsilon}\_{m3} \end{bmatrix} \mathbf{d}^{\varepsilon} = \mathbf{B}^{\varepsilon}\_{m} \mathbf{d}^{\varepsilon} \tag{16}$$

$$\mathbf{k}^{\varepsilon} = \begin{bmatrix} \mathbf{B}\_{b1}^{\varepsilon} & \mathbf{B}\_{b2}^{\varepsilon} & \mathbf{B}\_{b3}^{\varepsilon} \end{bmatrix} \mathbf{d}^{\varepsilon} = \mathbf{B}\_{b}^{\varepsilon} \mathbf{d}^{\varepsilon} \tag{17}$$

$$\mathbf{y}^{\varepsilon} = \begin{bmatrix} \mathbf{B}^{\varepsilon}\_{\varsigma 1} & \mathbf{B}^{\varepsilon}\_{\varsigma 2} & \mathbf{B}^{\varepsilon}\_{\varsigma 3} \end{bmatrix} \mathbf{d}^{\varepsilon} = \mathbf{B}^{\varepsilon}\_{\varsigma} \mathbf{d}^{\varepsilon} \tag{18}$$

where

$$\mathbf{B}\_{\rm mi}^{\varepsilon} = \begin{bmatrix} N\_{i,x} & 0 & 0 & 0 & 0 & 0 \\ 0 & N\_{i,y} & 0 & 0 & 0 & 0 \\ N\_{i,y} & N\_{i,x} & 0 & 0 & 0 & 0 \end{bmatrix} \tag{19}$$

$$\mathbf{B}\_{\text{li}}^{c} = \begin{bmatrix} 0 & 0 & 0 & N\_{i,\text{x}} & 0 & 0\\ 0 & 0 & 0 & 0 & N\_{i,\text{y}} & 0\\ 0 & 0 & 0 & N\_{i,\text{y}} & N\_{i,\text{x}} & 0 \end{bmatrix} \tag{20}$$

$$\mathbf{B}\_{si}^{\varepsilon} = \begin{bmatrix} 0 & 0 & N\_{i,\mathbf{x}} & N\_i & 0 & 0\\ 0 & 0 & N\_{i,\mathbf{y}} & 0 & N\_i & 0 \end{bmatrix} \tag{21}$$

By substituting the discrete displacement field into Equation (9) the discretized equation for free vibration analysis can be written into matrix form such as

$$(\mathbf{K} - \omega^2 \mathbf{M})\mathbf{\hat{d}} = 0,\tag{22}$$

where ω is the natural frequency, **K** and **M** are the global stiffness and mass matrices, respectively,

$$\mathbf{K} = \sum\_{e=1}^{n^e} \mathbf{T}^T \mathbf{K}^e \mathbf{T} \tag{23}$$

with

$$\mathbf{K}^{\varepsilon} = \int\_{\Omega\_{\varepsilon}} \left(\mathbf{B}^{\varepsilon}\right)^{T} \mathbf{D} \mathbf{B}^{\varepsilon} \mathrm{d}\Omega\_{\varepsilon} \tag{24}$$

and

$$\mathbf{B}^{\epsilon} = \begin{bmatrix} \mathbf{B}\_{\text{ff}}^{\epsilon} \\ \mathbf{B}\_{b}^{\epsilon} \\ \mathbf{B}\_{s}^{\epsilon} \end{bmatrix}, \mathbf{D} = \begin{bmatrix} \mathbf{D}^{\text{M}} & \mathbf{B} & 0 \\ \mathbf{B} & \mathbf{D}^{b} & 0 \\ 0 & 0 & \mathbf{C} \end{bmatrix} \tag{25}$$

$$\mathbf{M} = \sum\_{c=1}^{n^c} \mathbf{T}^T \mathbf{M}^c \mathbf{T},\tag{26}$$

$$\mathbf{M}^{\varepsilon} = \int\_{\Omega\_{\varepsilon}} \mathbf{N}^{T} \mathbf{m} \mathbf{N} d\Omega\_{\varepsilon} \tag{27}$$

in which **T** is the transformation matrix between the local coordinate system *Oxyz* and the global coordinate system *O* ˆ *x* ˆ *y* ˆ *z* ˆ [14].

The problem of zero stiffness that appears with using the drilling degree of freedom β*<sup>z</sup>*, which can cause a singularity in the global stiffness matrix when all the elements meeting at a node are coplanar. To deal with this issue, a simple modification coefficient is chosen to be 10−<sup>3</sup> times the maximum diagonal value of the element stiffness matrix at the zero drilling degree of freedom to avoid the drill rotation locking [15].

### **3. Formulation of ES-MITC3 Finite Element Method for FGM Shells**

### *3.1. Brief on the MITC3 Formulation*

In the linear triangular MITC3, the approximated displacement field **u** is simply interpolated using the linear basic functions for membrane, deflection, and rotation without adding any new variables. Herein, the membrane and bending strains of the standard finite elements are unchanged, while the transverse shearstrains, which are modified by the mixed interpolation of tensorial components [16].

As a result, the transverse shearstrain field [8,10] is being obtained as

$$\mathbf{y}\_{\text{MTTC},3}^{\varepsilon} = \begin{bmatrix} \mathbf{B}\_{s1}^{\varepsilon} & \mathbf{B}\_{s2}^{\varepsilon} & \mathbf{B}\_{s3}^{\varepsilon} \end{bmatrix} \mathbf{d}^{\varepsilon} = \mathbf{B}\_{s}^{\varepsilon} \mathbf{d}^{\varepsilon} \tag{28}$$

where

$$\mathbf{B}\_{s1}^{\epsilon} = J^{-1} \begin{bmatrix} 0 & 0 & -1 & \frac{6}{3} + \frac{4}{6} & \frac{6}{3} + \frac{6}{6} & 0\\ 0 & 0 & -1 & \frac{4}{3} + \frac{4}{6} & \frac{6}{3} + \frac{6}{6} & 0 \end{bmatrix} . \tag{29}$$

$$\mathbf{B}\_{s2}^{\varepsilon} = \boldsymbol{J}^{-1} \begin{bmatrix} 0 & 0 & 1 & \frac{a}{2} - \frac{d}{6} & \frac{b}{2} - \frac{c}{6} & 0 \\ 0 & 0 & 0 & \frac{d}{6} & \frac{c}{6} & 0 \end{bmatrix} \tag{30}$$

$$\mathbf{B}\_{s3}^{c} = J^{-1} \begin{bmatrix} 0 & 0 & 0 & \frac{a}{6} & \frac{b}{6} & 0\\ 0 & 0 & 1 & \frac{d}{2} - \frac{a}{6} & \frac{c}{2} - \frac{b}{6} & 0 \end{bmatrix} \tag{31}$$

with

$$J^{-1} = \frac{1}{2A^{\epsilon}} \begin{bmatrix} c & -b \\ -d & a \end{bmatrix} \tag{32}$$

in which *a* = *x*2 − *x*1, *b* = *y*2 − *y*1, *c* = *y*3 − *y*1, and *d* = *x*3 − *x*1, as pointed out in Figure 4 and *Ae* is the area of the triangular element.

**Figure 4.** Three-node triangular element and local coordinates.

### *3.2. The ES-MITC3 Formulation*

In the ES-FEM, the strains are smoothed over local smoothing domains <sup>Ω</sup>*k*, the computation for stiffness matrix is no longer based on elements, but on these smoothing domains. These smoothing domains are formed based on edge of elements such as Ω = ∪*<sup>n</sup><sup>k</sup> <sup>k</sup>*=1Ω*<sup>k</sup>* and Ω*<sup>i</sup>* ∩ Ω*j* = ∅ for *i* - *j*, in which *nk* is the total number of edges of all the elements. On a curved geometry of shell models, an edge-based smoothing domain Ω*<sup>k</sup>* associated with the inner edge *k* is created two sub-domains of two non-planar adjacent MITC3 triangular elements as shown in Figure 5. These triangular elements are defined by two local coordinate systems *O*1*x*1 *y*1*z*1 and *O*2*x*2 *y*2*z*2. In order to compute the edge-based smoothing strain Ω*<sup>k</sup>* for two non-planar adjacent elements, the virtual coordinate system *<sup>O</sup>xyz* is proposed as shown in Figure 6, whereas the *x*-axis coinciding with the edge *k*, the *z*-axis with the average direction between the *z*ˆ1-axis and *z*ˆ2-axis, and the *y*-axis is given by the cross-product of the unit vectors in the *x*-axis and *z*-axis.

**Figure 5.** The smoothing domain; Ω*<sup>k</sup>* is formed by triangular elements.

Hence, a smoothed membrane strain ε*k*, a smoothed bending strain κ*k*, a smoothed shearstrain γ*k* of the smoothing domain **Ω***<sup>k</sup>* in the global coordinate system *O*ˆ *x*<sup>ˆ</sup>*y*<sup>ˆ</sup>*z*<sup>ˆ</sup> can be derived as

$$
\overline{\boldsymbol{\mathfrak{e}}}^{k} = \sum\_{j=1}^{n^{nk}} \overline{\mathbf{B}}\_{mj}^{k} \mathbf{d}\_{j}^{k}; \\
\overline{\boldsymbol{\mathfrak{k}}}^{k} = \sum\_{j=1}^{n^{nk}} \overline{\mathbf{B}}\_{bj}^{k} \mathbf{d}\_{j}^{k}; \\
\overline{\mathbf{y}}^{k} = \sum\_{j=1}^{n^{nk}} \overline{\mathbf{B}}\_{bj}^{k} \mathbf{d}\_{j}^{k} \tag{33}
$$

where *nnk* is the number of the neighboring nodes of edge *k*. **d***kj* is the nodal degrees of freedom at the *j*th node of the smoothing domain Ω*<sup>k</sup>* in *O*ˆ *x*<sup>ˆ</sup>*y*<sup>ˆ</sup>*z*ˆ. **<sup>B</sup>***kmj*, **B***kbj*, and **B***ksj* are the membrane, the bending and the MITC3 shear smoothed gradient matrices at the *j*th node of the smoothing domain Ω*<sup>k</sup>* in the global coordinate system *O* ˆ *x* ˆ *y* ˆ *z* ˆ, respectively. The **B** *k mj*, **B** *k bj* and **B** *k sj* can be computed by

$$\overleftarrow{\mathbf{B}}\_{mj}^{k} = \frac{1}{A^{k}} \sum\_{i=1}^{n^{k}} \frac{1}{3} A^{i} \boldsymbol{\Lambda}\_{m1}^{k} \boldsymbol{\Lambda}\_{m2}^{i} \mathbf{B}\_{mj}^{i} \mathbf{T}\_{j}^{i} \tag{34}$$

$$\overleftarrow{\mathbf{B}}\_{b\circ j}^{k} = \frac{1}{A^{k}} \sum\_{i=1}^{n^{k}} \frac{1}{\mathfrak{Z}} \mathcal{A}^{i} \mathbf{A}\_{b1}^{k} \mathbf{A}\_{b2}^{i} \mathbf{B}\_{bj}^{i} \mathbf{T}\_{j}^{i} \tag{35}$$

$$\overline{\mathbf{B}}\_{sj}^{k} = \frac{1}{A^k} \sum\_{i=1}^{n^k} \frac{1}{3} A^i \mathbf{A}\_{s1}^k \mathbf{A}\_{s2}^i \mathbf{B}\_{sj}^i \mathbf{T}\_j^i \tag{36}$$

in which **<sup>Λ</sup>***km*1, **<sup>Λ</sup>***kb*1, and **<sup>Λ</sup>***ks*1 are strain transformation matrices between the global coordinate system *O* ˆ *x* ˆ *y* ˆ *z* ˆ and the virtual coordinate system *O x y z*, respectively; **<sup>Λ</sup>***im*2, **<sup>Λ</sup>***ib*2 and **<sup>Λ</sup>***is*2 are the strain transformation matrices between the local coordinate system *Oxyz* of *i*th adjacent triangular elements and the virtual coordinate system *O x y z*, respectively; **T***ij* is the transformation matrix between the local coordinate system *Oxyz* at the *j*th node of the *i*th adjacent triangular element and the global coordinate system *O* ˆ *x* ˆ *y* ˆ *z* ˆ. More detailed information about these strain transformation matrices can be found in [14]. The area *A<sup>k</sup>* of the smoothing domain Ω*<sup>k</sup>* is computed by

$$A^k = \int\_{\Omega^k} \mathrm{d}\Omega = \frac{1}{3} \sum\_{i=1}^{n^k} A^i \tag{37}$$

where *nek* is the number of the adjacent triangular elements in the smoothing domain <sup>Ω</sup>*k*, and *A<sup>i</sup>* is the area of the *i*th triangular element around the edge *k*.

**Figure 6.** Global, local, and virtual coordinates.

As a result, the global stiffness matrix of the FGM shell in Equation (22) is rewritten as

$$(\overline{\mathbf{K}} - \omega^2 \mathbf{M})\hat{\mathbf{d}} = 0,\tag{38}$$

$$
\widetilde{\mathbf{K}} = \sum\_{i=1}^{n^{\text{v}}} \widetilde{\mathbf{K}}^{k} \tag{39}
$$

where

$$\widetilde{\mathbf{K}}^k = \int\_{\Omega\_k} \left(\overline{\mathbf{B}}^k\right)^T \mathbf{D} \overline{\mathbf{B}}^k d\Omega\_k \tag{40}$$

with

$$
\overleftarrow{\mathbf{B}}^{k} = \begin{bmatrix} \overleftarrow{\mathbf{B}}\_{m}^{k} & \overleftarrow{\mathbf{B}}\_{b}^{k} & \overleftarrow{\mathbf{B}}\_{s}^{k} \end{bmatrix}^{T} \tag{41}
$$

### **4. Numerical Results**

In the section, several numerical examples are provided to show the performance of the ES-MITC3 element for free vibration analysis of FGM shell and results obtained are compared to those published [6, 16–19]. For convenience to the numerical comparison, the non-dimensional frequency parameters ω<sup>∗</sup> are expressed to the following equation as

$$
\omega^\* = \omega \mathfrak{a}^2 \sqrt{\rho\_m \mathfrak{h} / D\_{m\nu}^\*} \ D\_m^\* = E\_m \mathfrak{h}^3 / 12 (1 - v\_m^3) \tag{42}
$$

First, let us consider free vibration for analysis of clamped functionally graded cylindrical shell (*Rx* = *R*, *Ry* = ∞) with radius-to-length *R*/*a* = 100, *a*/*h* = 10. The functionally graded shell is made from Silicon nitride (Si3N4) and Stainless steel (SUS304), which material properties are *Ec* = 322.2715 GPa, *vc* = 0.24, ρ*c* = 2370 kg/m3, *Em* = 207.7877 GPa, *vm* = 0.31776, and ρ*m* = 8166 kg/m3. The first four non-dimensional frequency of the present method list in Table 1 are compared with MITC3 [16], a higher-order theory based on radial basis functions collocation including transverse normal deformation (HSDT RBFC-1) and discarding transverse normal deformation (HSDT RBFC-2) by Neves et al. [17], a higher-order theory and finite element formalation (HSDT FEM) by Pradyumna and Bandyopadhyay [6], a higher-order theory and semi-analytical method relied on Galerkin (HSDT SAG) by Yang and Shen [18], and Quasi-3D Ritz model (ED555) by Fazzolari and Carrera [19]. From Table 1 we can see that this proposed method (ES-MITC3) is more accurate than other methods, such as MITC3, HSDT RBFC-1, HSDT RBFC-2, HSDT FEM and HSDT SAG. The errors are less than 3% in comparison with the exact solution ED555 [19]. Figure 7 shows non-dimensional frequency parameter for four first modes of clamped functionally graded cylindrical shell using various methods.


**Table 1.** Non-dimensional frequency parameter for clamped cylindrical functionally graded materials (FGM) shell with *R*/*a* = 100, and relative error between methods (ED555 [19] is fixed). Error (%) = <sup>100</sup>×|Method−ED555[19]| ED555[19].

**Figure 7.** Non-dimensional frequency parameter for four first modes. (**a**) Mode 1; (**b**) mode 2; (**c**) mode 3; (**d**) mode 4.

Next, we investigate the first non-dimensional frequencies ω\*. of functionally graded spherical (*Rx* = *Ry* = *R*) and cylindrical shells (*Rx* = *R*, *Ry* = ∞) with geometric data: radius to edge *R*/*a* and *a*/*h* are varied from 5 to 50 and 10, respectively. The functionally graded shells in these studies are made from aluminum, and alumina whose material properties are *Em* = 70, GPa, *vm* = 0.3, ρ*m* = 2707 kg/m3, *Ec* = 380 GPa, *vc* = 0.3, and ρ*c* = 3000 kg/m3. Again, it is seen from Tables 2–5 that the results of the present approach are very close to an HSDT RBFC-1, HSDT RBFC-2 [17], and ED555 [19]. Figures 8–11 show non-dimensional frequency parameter for the first mode of cylindrical FGM shell and spherical FGM shell with different *n*, respectively. The first six mode shapes of simply supported spherical FGM shell are illustrated in Figure 12.

**Figure 8.** Non-dimensional frequency parameter for the first mode of clamped cylindrical FGM shell with different *n*. (**a**) R/a = 5; (**b**) R/a = 10; (**c**) R/a = 50.

**Table 2.** Non-dimensional frequency parameter for clamped cylindrical FGM shell with *a*/*h* = 10 and different *R*/*a* ratios.



**Table 3.** Non-dimensional frequency parameter for simply supported cylindrical FGM shell with *a*/*h* = 10 and different *R*/*a* ratios.

**Figure 9.** Non-dimensional frequency parameter for the first mode of simple supported cylindrical FGM shell with different *n*. (**a**) R/a = 5; (**b**) R/a = 10; (**c**) R/a = 50.


**Table 4.** Non-dimensional frequency parameter for simply supported spherical FGM shell with *a*/*h* = 10 and different *R*/*a* ratios.

**Figure 10.** Non-dimensional frequency parameter for the first mode of simply supported spherical FGM shell with different *n*. (**a**) R/a = 5; (**b**) R/a = 10; (**c**) R/a = 50.


**Table 5.** Non-dimensional frequency parameter for clamped spherical FGM shell with different *R*/*a* ratios.

**Figure 11.** Non-dimensional frequency parameter for the first mode of clamped spherical FGM shell with different *n*. (**a**) R/a = 5; (**b**) R/a = 10; (**c**) R/a = 50.

**Figure 12.** First six mode shapes of simply supported spherical FGM shell (R/a = 10, a/h = 10, *n* = 0.2).

For the fully clamped spherical shell, the second vibration mode shape and the third vibration mode shape are similar to each other (their natural frequencies are equal), they just have different views. Adding, the value of non-dimensional frequency of the fifth and the sixth vibration mode shape are approximatively each other. This thing is consistent with the actual symmetrical shell structures with the same boundary conditions.
