Free Vibration Analysis of Functionally Graded Shells Using an Edge-Based Smoothed Finite Element Method

: An edge-based smoothed ﬁnite element method (ES-FEM) combined with the mixed interpolation of tensorial components technique for triangular shell element (MITC3), called ES-MITC3, for free vibration analysis of functionally graded shells is investigated in this work. In the formulation of the ES-MITC3, the sti ﬀ ness matrices are obtained by using the strain-smoothing technique over the smoothing domains that are formed by two adjacent MITC3 triangular shell elements sharing an edge. The strain-smoothing technique can improve signiﬁcantly the accuracy and convergence of the original MITC3. The material properties of functionally graded shells are assumed to vary through the thickness direction by a power–rule distribution of volume fractions of the constituents. The numerical examples demonstrated that the present ES-MITC3method is free of shear locking and achieves the high accuracy compared to the reference solutions in the literature.


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 stiffened functionally graded double curved shells on an

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 where P(z) represents the effective material properties: Young's modulus E, density ρ and Poisson ratio v; P c and P m denote the properties of the ceramic and metal, respectively; V c 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. variation of the volume fraction of ceramic and metal through the thickness via the volumefraction exponent n .

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 0 0 , u v , and 0 w be the displacements of the middle plane in the , x y , and z directions; , x y β β , 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 The linear strain-displacement relationship can be given as

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 0 0 , u v , and 0 w be the displacements of the middle plane in the , x y , and z directions; , x y β β , 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 The linear strain-displacement relationship can be given as

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; β x , β y , 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 The linear strain-displacement relationship can be given as From Hooke's law, the constitutive relations of FGM shells are expressed as: where A weak form of the free vibration analysis for FGM shells can be briefly given as:  In which 11 In Equation (11) A ij , B ij , D ij , and C ij are given by: In which In Equation (11) A ij , B ij , D ij , and C ij are given by: where k = 5/6 is transverse shear correction coefficient and m is the mass matrix containing ρ calculated as

Finite Element Formulation for Shell Analysis
Now, we discretize the bounded domain Ω into n e finite three-node triangular elements with n n nodes such that Ω ≈ n e e=1 Ω e and Ω i ∩ Ω j = ∅, i j. The displacement field u e = u e v e w e β e x β e y β e z T of the finite element solution can be expressed as where d e i = u e 0i v e 0i w e 0i β e xi β e yi β e zi T is the nodal displacement at the ith node; N i (x) is shape function for the ith node.
The approximation of the membrane, the bending and the shear strains of the triangular element can be written in matrix forms as follows where By substituting the discrete displacement field into Equation (9) the discretized equation for free vibration analysis can be written into matrix form such as where ω is the natural frequency, K and M are the global stiffness and mass matrices, respectively, with and in which T is the transformation matrix between the local coordinate system Oxyz and the global coordinate systemÔxŷẑ [14]. The problem of zero stiffness that appears with using the drilling degree of freedom β z , 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 −3 times the maximum diagonal value of the element stiffness matrix at the zero drilling degree of freedom to avoid the drill rotation locking [15].

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 where with 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 A e is the area of the triangular element.

The ES-MITC3 Formulation
In the ES-FEM, the strains are smoothed over local smoothing domains O x y z . In order to compute the edge-based smoothing strain k Ω for two non-planar adjacent elements, the virtual coordinate system Oxyz   is proposed as shown in Figure 6, whereas the -axis coinciding with the edge , the -axis with the average direction between the 1 z -axis and 2 z -axis, and the y  -axis is given by the cross-product of the unit vectors in the x  -axis and z  -axis.

The ES-MITC3 Formulation
In the ES-FEM, the strains are smoothed over local smoothing domains Ω 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 Ω = ∪ n k k=1 Ω k and Ω i ∩ Ω j = ∅ for i j, in which n k is the total number of edges of all the elements. On a curved geometry of shell models, an edge-based smoothing domain Ω k 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 Ω k for two non-planar adjacent elements, the virtual coordinate system O x y z 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ẑ 1 -axis andẑ 2 -axis, and the y-axis is given by the cross-product of the unit vectors in the x-axis and z-axis.

The ES-MITC3 Formulation
In the ES-FEM, the strains are smoothed over local smoothing domains O x y z . In order to compute the edge-based smoothing strain k Ω for two non-planar adjacent elements, the virtual coordinate system Oxyz   is proposed as shown in Figure 6, whereas the -axis coinciding with the edge , the -axis with the average direction between the 1 z -axis and 2 z -axis, and the y  -axis is given by the cross-product of the unit vectors in the x  -axis and z  -axis.
in which Λ k m1 , Λ k b1 , and Λ k s1 are strain transformation matrices between the global coordinate systemÔxŷẑ and the virtual coordinate system O x y z, respectively; Λ i m2 , Λ i b2 and Λ i s2 are the strain transformation matrices between the local coordinate system Oxyz of ith adjacent triangular elements and the virtual coordinate system O x y z, respectively; T i j is the transformation matrix between the local coordinate system Oxyz at the jth node of the ith adjacent triangular element and the global coordinate systemÔxŷẑ. More detailed information about these strain transformation matrices can be found in [14]. The area A k of the smoothing domain Ω k is computed by where n ek is the number of the adjacent triangular elements in the smoothing domain Ω k , and A i is the area of the ith 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

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][17][18][19]. For convenience to the numerical comparison, the non-dimensional frequency parameters ω * are expressed to the following equation as First, let us consider free vibration for analysis of clamped functionally graded cylindrical shell (R x = R, R y = ∞) with radius-to-length R/a = 100, a/h = 10. The functionally graded shell is made from Silicon nitride (Si 3 N 4 ) and Stainless steel (SUS304), which material properties are E c = 322.2715 GPa, v c = 0.24, ρ c = 2370 kg/m 3 , E m = 207.7877 GPa, v m = 0.31776, and ρ m = 8166 kg/m 3 . 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 (ED 555 ) 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 ED 555 [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 (ED 555 [19] is fixed).
Error (%) = 100×|Method−ED 555 [19]| ED 555 [19] .  Next, we investigate the first non-dimensional frequencies ω*. of functionally graded spherical (R x = R y = R) and cylindrical shells (R x = R, R y = ∞) 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 E m = 70, GPa, v m = 0.3, ρ m = 2707 kg/m 3 , E c = 380 GPa, v c = 0.3, and ρ c = 3000 kg/m 3 . 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 ED 555 [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.        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.

Conclusions
In this paper, the free vibration analysis of functionally graded shells is studied by using the ES-MITC3. Herein, the stiffness matrices obtained based on the strain-smoothing technique over the smoothing domains associated with edges of MITC3 shell elements. The present approach uses a triangular element and hence much easily generated automatically, even for complicated geometries. The numerical results showed that the ES-MITC3 is a good agreement with the reference solutions, which are a requirement of high computational costs, such as ED 555 [19] and a higher-order theory based on radial basis functions [17]. The ES-MITC3 presented herein is promising to be a simple and effective finite element method for analysis of functionally graded shells in practice.
The combination of finite element method (FEM) with an edge-based smoothed finite element method (ES-FEM) is very suitable for analyzing shell structures, especially for the complicated structures such as shell with reinforced stiffeners, reinforced nano grapheme, micro shell structures, nano shell structures, and so on. This combination allows us to calculate exactly plate and shell structures with thin thicknesses (h = a/10 8 ) due to overcoming the shear locking phenomenon.

Conflicts of Interest:
The authors declare no conflict of interest.