Nonlinear Dynamics Analysis of FGM Shell Structures with a Higher Order Shear Strain Enhanced Solid-Shell Element 1

In this paper, non-linear dynamics analysis of functionally graded material (FGM) shell structures is investigated using the higher order solid-shell element based on the Enhanced Assumed Strain (EAS). With this element, a quadratic distribution of the shear stress through the thickness is considered in an enhancing part. Material properties of the shell structure are varied continuously in the thickness direction according to the general four-parameter power-law distribution in terms of the volume fractions of the constituents. Performance and accuracy of the present higher order solid-shell element are confirmed by comparing the numerical results obtained from finite element analyses with results from the literature.


Latin American Journal of Solids and
FGMs are varied continuously in one or more directions.Therefore, the stress distributions are smooth hence interface problems are eliminated when compared to laminate structure In order to avoid structural failure caused by dynamic loadings, linear dynamic characteristics of FGMs shell structures have been considered by many researchers.Finding analytic solutions for linear and non-linear dynamics analysis of functionally graded material is a very difficult task, and they only exist for particular cases.Meanwhile, the finite-element (FE) method has become the most widely used technique to model the processes of dynamics analysis.
There are several solution formulations for vibration and linear dynamic analysis of FGM shells structures.They can be classified into two groups: the 2D plate and shell models and the full 3D elasticity model.Zhi-Yuan and Hua-Ning (2007) studied free vibration characteristics of functionally graded cylindrical shells with holes using Kirchhoff Classical thin Plate Theory (CPT).On the other hand, Kadoli and Ganesan (2006) studied the buckling and free vibration analysis of functionally graded cylindrical shells subjected to a temperature-specified boundary condition.In their analysis, the finite element equations based on the First order Shear Deformation Theory (FSDT) and the Plane Stress Assumption (PSA) were formulated.Recently, Ansari et al. (2016) developed a nonclassical size-dependent plate model based on the modified strain gradient and the FSDT for the bending, buckling and free vibration analyses of microscale FG plates.
Based on the Third order Shear Deformation Theory (TSDT) and PSA, Reddy (2000) presented a theoretical formulation, Navier's solutions of rectangular plates, and finite element models to study the nonlinear dynamic response of FGM plates subjected to a suddenly applied uniform pressure.Yang and Shen (2002) analyzed free and forced vibration for initially stressed FGM plates.In this contribution, theoretical formulations are based on TSDT and PSA including the thermal effects.One-dimensional differential quadrature technique and Galerkin approach are used to determine the transient response of the plate subjected to dynamic loads.Gharooni and Ghannad (2015) investigated displacements and stresses in pressurized thick FGM cylinders with exponential variation of material properties based on TSDT.Furthermore, based on a TSDT of shells and PSA, Wali et al. (2015) Frikha et al. (2016) and Frikha et al. (2017) studied respectively the free vibration, linear dynamic response and fully geometrical nonlinear mechanical response of FGM shell structures using an efficient double directors shell element proposed in Wali et al. (2014).The shear stress boundary conditions on top and bottom faces are considered in a discrete form as given in (2005).
For full 3D elasticity formulations for vibration and linear dynamic analysis of FGM shells, one find the work of Vel and Batra (2004) based on three dimensional exact solutions for free and forced vibrations of simply supported FGM rectangular plates.In Asemi et al. (2014) the static and dynamic analyses of FGM skew plates are obtained based on the three-dimensional theory of elasticity.Graded elements, the principle of minimum energy and Rayleigh-Ritz energy method are used.Using 3D elasticity model, Nguyen and Nguyen-Xuan (2015) proposed an efficiently computational tool based on an isogeometric finite element formulation for static and dynamic response analysis of FGM plates.
Regarding the nonlinear dynamics analysis of FGM shells, three kinematics assumption are used, CPT, FSDT and TSDT.Structures 14 (2017) 72-91 Based on the CPT assumption and von-Karman geometrical nonlinearity, Woo et al. (2006) investigated the nonlinear free vibration behavior of FGM plates.Allahverdizadeh et al. (2008) evaluated the material properties of an FGM thin circular plate and investigated the nonlinear free and forced vibrations using the shooting technique.The formulation is based on CPT and von-Karman geometrical nonlinearity.Alijani et al. (2011) studied the nonlinear vibrations of FGM doubly curved shallow shells.They considered the thermal effect, the CPT assumption and Donnell nonlinearity.Using the CPT with an improved Donnell equations and PSA, Bich and Nguyen (2012) examined the nonlinear vibration of functionally graded circular cylindrical shells.Considering the CPT with PSA and von-Karman geometrical nonlinearity, Duc (2013) presented an analytic investigation on the nonlinear dynamic response of eccentrically stiffened functionally graded double curved shallow shells resting on elastic foundations and being subjected to axial compressive load and transverse load.Using the same considerations, Duc and Cong (2015) investigated the nonlinear dynamic response of imperfect symmetrical thin FGM plate on elastic foundation.

Latin American Journal of Solids and
Using the FSDT plate finite element method with PSA employing von-Karman nonlinearity, Praveen and Reddy (1998) investigated the nonlinear transient thermo-elastic response of FGM plates.Liew et al. (2006) studied the nonlinear vibration of a coating-FGM-substrate cylindrical panel subjected to a temperature gradient.The FSDT is considered with PSA and von-Karman geometrical nonlinearity.Also based on the FSDT with PSA and von-Karman geometrical nonlinearity, Zhang et al. (2012) analyzed the nonlinear dynamics of a clamped-clamped FGM circular cylindrical shell subjected to an external excitation and uniform temperature change.
Based on TSDT and PSA, Hao et al. (2008) presented an analysis of the nonlinear dynamics of a simply supported FGM rectangular plate subjected to the transversal and in-plane excitations in a thermal environment.The von-Karman geometrical nonlinearity assumption is used.Using the same von-Karman geometrical nonlinearity with TSDT and PSA, Duc et al. (2015) presented an analytical approach to investigate the nonlinear dynamic response and vibration of imperfect FGM thick circular cylindrical shells surrounded on elastic foundation.Based on TSDT with PSA and von-Karman type nonlinear kinematics, Liu et al. (2015) presented a nonlinear dynamic analysis of a slightly initial imperfect FGM circular cylindrical shell subjected to complex loads including aerodynamic pressure and thermal loading.
It can be seen from the previous literature that in most of the studies, full 3D elasticity formulations are only limited to linear dynamic and free vibration of FGM structures.Also, the von-Karman or Donell geometrical nonlinearity are the only kinematic that has been basically used for the nonlinear dynamic analysis of shell's structures where the PSA is largely considered.
However, in this paper we focus on the fully 3D non-linear dynamics analysis of FGM shell structures by using higher order shear strain enhanced solid-shell element.The present formulation constitute an extension of the higher order shear deformation solid-shell finite element, developed by Hajlaoui et al. (2016), to the full nonlinear dynamics.The present solid-shell element formulation is based on the partition of shear strain: one of the parts is independent of the thickness coordinate and formulated by the Assumed Natural Strain (ANS) method, which avoids the shear locking in thin limit structure; the other enhancing part ensures a quadratic distribution across the thickness.With this shear strain enhancement, the accuracy of transverse shear stresses will be improved and the shear correction factors will be avoided.Structures 14 (2017) 72-91 The remainder of this paper is organized as follows.Functionally graded materials are described in section two.After that, solid-shell finite element formulation and a transient analysis of the nonlinear formulation is described in section three and four respectively.Numerical results and discussions of the finite element model are investigated in detail in section five.Finally, some concluding remarks are analyzed and presented in section six.

FUNCTIONALLY GRADED MATERIALS
In this paper, we consider an FGM shell structures made from a mixture of metal and ceramics and the composition varies continuously in the thickness direction.In fact, the Young's modulus ( ) E z , density ( ) z r and Poisson's ratio ( ) are assumed to vary through the shell thicknesses according to a power-law distribution as in which the subscripts m and c refer to metal and ceramic components, respectively.In addition, the volume fraction c V follows two general four-parameter power-law distributions, Su et al. (2014).( ) ( ) where a, b and c are the parameters which determine the material variation profile through the FGM shell thickness and p is the power-law index.

SOLID-SHELL FINITE ELEMENT FORMULATION
The developed solid-shell element is an eight nodes hexahedral element with three degrees of freedom per node.The transverse shear strain is composed of two parts.The first one is independent of the thickness coordinate and formulated by the assumed natural strain method (ANS).The second part is an enhancing part, which ensures a quadratic distribution through the thickness.
The enhanced assumed strain method consists in the enhancement of the compatible part of the Green Lagrange strain tensor, c E , with an enhanced part E  to have a total strain as follows

Variational Formulation
The point of departure is the well-known three-field variational functional in Lagrangean formulation.This variational functional is as follows Latin American Journal of Solids and Structures 14 (2017) 72-91 ( ) ( ) where y is the strain energy function and u, E  and S  are the independent tonsorial quantities which are: displacement, enhanced assumed strain and assumed stress fields respectively.Also in Eq. ( 5) appear the prescribed body force V F and surface traction S F .The orthogonality between the enhanced strain and stress fields leads to This orthogonality condition reduce the number of independent variables in the original functional to just two ( ) , u E  .The weak form of this modified reduced functional may be obtained with the direction derivative leading to where S is the Piola-Kirchhoff stress tensor

Finite Element Formulation
In each finite element domain, an eight-node hexahedral solid-shell element is considered.The position vectors in reference and current configurations respectively are where N is the tri-linear shape functions matrix given by x and n X are nodal coordinates.The displacement field, with the corresponding variation and increment, is interpolated in the same manner as follows where is the nodal displacements vector at the element level.The covariant base vectors obtained by partial derivative of the position vectors with respect to convective coordinate ( ) ( ) in reference and current configuration are given by Latin American Journal of Solids and Structures 14 (2017) 72-91 The covariant metric tensor at a material point ξ , in the reference and current configuration are defined by This leads to the following Green-Lagrangean strain tensor ( )

Compatible Strains
To treat the transverse shear locking and transverse normal locking problems, the ANS method is used.For the transverse shear strains, 13 c E and 23 c E , the ANS method, proposed by Bathe and Dvorkin (1985), is used and evaluated at four mid-points of the element edges A = (-1,0,0), B = (0,-1,0), C = (1,0,0), D = (0,1,0), Fig. 1.The transverse shear strains are given by ) ( )( ) 13 13 13 13 13 23 13 13 13 13 However, for the thickness strains, 33 c E , we adopt the ANS method as in Klinkel et al. (1999), Vu-Quoc and Tan (2003) and Hajlaoui et al. (2012Hajlaoui et al. ( , 2015Hajlaoui et al. ( , 2016) and evaluated at four collocation points defined in the reference surface ( ) ) ( ) For the detailed description of the ANS method, refer to Hajlaoui et al. (2012Hajlaoui et al. ( , 2015Hajlaoui et al. ( , 2016)).The compatible part is arranged in (6x1) column matrix as follows Using Eqs. ( 15)-( 17) and approximations Eqs. ( 9) and ( 11), the virtual and incremental compatible Green Lagrange strain tensor are then given by where B is the strain interpolation matrix, associated to node I ( ) Matrix T, in Eq. ( 19), is the transformation of the strain tensor from parametric coordinates to the local Cartesian coordinates, written as where .ij i j t = G T and j T (j = 1, 2, 3) are a local orthonormal base vectors.

Enhanced Strains
The enhanced Green-Lagrange strain part is related to the internal strain parameters vector α as: where E  , dE  and DE  are total, virtual and incremental enhanced Green Lagrange strain tensor respectively.The crucial assumption of the EAS method is the enforcement of the orthogonality conditions for the assumed stress field S  and the enhanced strain E  .This orthogonality conditions impose the following choice for the interpolation function matrix M  to be expressed as follows where the subscript '0' means the evaluation at the center of the element in the natural coordinates and G G G is the Jacobian matrix.Eq. ( 22) responds a minimum requirement that the enhancing strain be orthogonal to the constant stress state.The interpolation matrix xhz M , in Eq. ( 22), is expression in term of the parametric coordinates( In all these three matrices, the term ever, the three matrices are different in the plane and thickness strains enhancement.In a previous work, Hajlaoui et al. (2016), the analysis of stiffness matrix eigenvalue in the incompressible range demonstrate that only the 11 parameters solid shell element, Eq. (23.c), provided the correct eigenvalues.As FGM materials are far from incompressible conditions as in rubber or elastoplasticity, and to reduce the computational time, one can choose the nine parameters, C3D8C9 element, which will be the only considered in the numerical results.

Linearization
With the finite element approximation, Eq. ( 18) and ( 21), at hand, the continuum weak form, Eq. ( 7), became in a discrete form as ( ) where, int f , ext f and h , are given by the following expressions Equation ( 24) is a nonlinear equation that will be solved iteratively by the Newton-Raphson method.For this purpose one need the corresponding linearization ( ) , 0 where L , H and K are given by


is the 6 6  three dimensional material tangent moduli, or elasticity tensor, and D K is given by where the matrix B is given by Eq. ( 19) and G K is the geometric stiffness matrix Relative to a couple of nodes (I, J) the matrix G K is given by Latin American Journal of Solids and Structures 14 (2017) 72-91 The strain parameters D must be eliminated from Eq. ( 26) at the element level, which leads to following element tangent operator T K and residual vector R ( )

TRANSIENT ANALYSIS OF THE NON-LINEAR FORMULATION
In order to extend the variational formulation, Eq. ( 24), to accommodate transient analysis, the body force is replaced by Then, with this replacement the variational equation became The acceleration in Eq. ( 33) is computed from the isoparametric interpolation as in Eq. ( 11).Thus, the inertia term and residual vector may be written as where M is the consistent mass matrix.The Newmark method is a one-step method, which is used to advance the solution from time tn to tn+1.Firstly, given the initial displacement and velocity vectors, the initial acceleration is determined by solving Eq. ( 34).The Newmark formula to process the solution is given by with the updating formulas where , and are displacement, velocity and acceleration vectors given in the initial state by Latin American Journal of Solids and Structures 14 (2017) 72-91 The Newmark parameters b and g are chosen as 0.25 b = , 0.5 g = .

NUMERICAL RESULTS AND DISCUSSION
The finite element formulation of the present higher order shear strain enhanced solid-shell element for nonlinear dynamic analyses of functionally graded material shell structures as presented in the previous sections.In this section, a comprehensive investigation concerning the nonlinear dynamic response of functionally graded material shell structures with various boundary conditions is given to demonstrate the accuracy and reliability of the present higher order solid-shell element.

Pull-Out of an Open-Ended Isotropic Cylindrical Shell
The purpose of this example is to demonstrate the ability of the formulation to capture large displacement and illustrate the robustness of the numerical implementation in static analysis.A short cylinder, with two pinching vertical forces at the middle section is modeled using one octant Fig. 2a.´, 0.3125 n = . The concentrated load P=40000 is applied at point A. The automatic incremental/iterative Newton procedure is used with a total of 64 load steps.The final reel deformed configuration is given at Fig.

Shallow Spherical Cap Under a Concentrated Load
In this sub-section, we present numerical simulations in order to illustrate the good performance of the proposed formulation.The dynamic behavior of a clamped isotropic spherical shell under a concentrated apex load, F=100, is presented.This test is proposed in Duarte Filho and Awruch ( 2004).
Material and geometric properties for this test are

FGM Plate Under Constant Distributed Load
This example consists on a square functionally graded plate subjected to a constant uniformly distributed loading of Figs. 6 and 7 show the dynamic response of the Aluminum-Zirconia plates simply supported and clamped boundary conditions respectively.The accuracy of the present model is verified by comparing the computed non-dimensionalized center deflection with available results of Praveen and Reddy (1998).It is clear that the presently computed non-dimensionalized center deflection is in general in good agreement with Praveen and Reddy (1998) finite element solutions.However, a discrepancy exists for a higher magnitude of deflection in the simply supported plate, Fig. 6.This can be explained by the major difference between the present formulation and Praveen and Reddy (1998): the present formulation is full 3D enhanced Green-Lagrange formulation and Praveen and Reddy (1998)

Clamped Circular FGM Plate Under Distributed Transverse Load
This example, investigated in Nguyen and Nguyen-Xuan (2015), consists on clamped circular FGM plates under uniform pressure 10 6 N/m 2 .The radius is R=0.5 m and the thickness is h=0.1 m.The circular FGM plate is made of the same constituent materials as the previous example, Aluminum-Zirconia.The used law distribution of the FGM plate is the classical power-law which can be obtained by the general four-parameter power-law distribution Eq. ( 3): ( ) Fig. 9 illustrates the central deflection dynamic response of the clamped FGM circular plate.In this figure, ANSYS solutions, for full metal and full ceramic, using quadratic finite elements and approximately 140436 degrees of freedom, are given in Nguyen and Nguyen-Xuan (2015).The present solutions, using only 1302 degrees of freedom, match very well with the ones obtained using ANSYS software.In addition, Fig. 9 shows the results obtained using several values of the power low index, 0.2, 0.5, 1 and 2. The present results are in good agreement with the published ones given in Nguyen and Nguyen-Xuan (2015) using 2535 degrees of freedom.Fig. 13 illustrate the effect of the power law index p on the nonlinear dynamic response of ZrO2/Ti-6Al-4V FGM cylindrical shells with / 2 L R = , / 500 R h = and ( ) 1500 sin(600 ) q t t = .
As remarked in previous test, the amplitude of the nonlinear dynamic of the FGM shells increases when increasing the power low index p.Fig. 14 gives the effect of the thickness ration R/h on the nonlinear response of FGM cylindrical shell with / 2 L R = , 2 p = and ( ) 1500 sin(600 ) q t t = . As had been anticipated, the amplitude of the FGM cylindrical shells increase when the ration R/h increase.

CONCLUSION
This paper presents a numerical investigation on the nonlinear dynamic response of FGM shell structures.Material properties of the shell structure are assumed to vary continuously through the thickness according to the general four-parameter power-law distribution in terms of the volume fractions of the constituents.The finite element formulation is based on the higher order shear strain enhanced solid-shell element.To avoid locking in the thin shell cases, the ANS method, for the transverse shear strains and thickness strain, is used instead of the continuum components in the compatible part of strain tensor.In the enhanced strain part, parabolic functions in terms of natural thickness coordinate are considered for the transverse shear strains.This permit to remove shear correction factors and improves the accuracy of transverse shear stresses.In addition, the enhancements of membrane part and the transverse strain are considered to avoid the in-plane bending and volumetric locking.
In the present solid-shell element, the nonlinear dynamics of FGM shell is done by the linearization of the nonlinear dynamic weak form.This was performed for use in the iterative solution for the kinematics quantities via the Newton's method.A time stepping algorithm based on the Newmark's scheme is employed in our time discrete weak form.From a demonstrative numerical simulation of complex FGM structures, it is found that the proposed method shows accurate results with good convergence characteristics.
Latin American Journal of Solids and Structures 14 (2017) 72-91 Results of the present formulations are verified by comparing the present results with those available in literatures.The obtained results show the effects of material, boundary conditions, excitation and geometrical parameters on the dynamical response of FGM structures.Thus, it is obvious that the dynamic response of the considered FGM structures depends on many factors significantly: power low index, boundary conditions, excitation and geometrical parameters of the FGM structures.
In addition, the present higher order solid-shell model can be used in complex FGM or composite shell's structures analysis, unlike the analytic formulation that can be employed only for smooth differentiable geometry of shell structures.

Figure 1 :
Figure 1: Transverse shear strain and thickness strain interpolation points.

Fig. 1 .
Fig. 1.The thickness strains is given by the parabolic transverse shear distribution in thickness direction.This enhancement requires a numerical integration rule at least 2x2x3.How-Latin American Journal of Solids and Structures 14 (2017) 72-91

Figure 2 :
Figure 2: Pull-out of an open-ended cylindrical shell a/ Finite element model, b/ Reel deformed configuration at P=40000.

Figure 3 :
Figure 3: Pull-out force P vs. radial deflections at points A, B and C for the open-ended cylindrical shell.
, and h = 0.01576.Loading and finite element mesh are shown in Fig.4.
test has been addressed inPraveen and Reddy (1998) among others.The side length or the plate is a=b=0.2mand the thickness h = 0.01m.The plate is made of the constituent materials: aluminum (bottom surface) and zirconia (top surface).Young's modulus, Poisson's ratio and density for the aluminum are 70 used law distribution of the FGM plate is the classical power-law which can be obtained by the general four-parameter power-law distribution Eq.(3): noting the symmetry, only one quadrant of the plate is modeled with 8x8x1 elements.The time increment is

Figure 6 :
Figs.6 and 7show the dynamic response of the Aluminum-Zirconia plates simply supported and clamped boundary conditions respectively.The accuracy of the present model is verified by comparing the computed non-dimensionalized center deflection with available results ofPraveen and Reddy (1998).It is clear that the presently computed non-dimensionalized center deflection is in general in good agreement with Praveen and Reddy (1998) finite element solutions.However, a discrepancy exists for a higher magnitude of deflection in the simply supported plate, Fig.6.This can be explained by the major difference between the present formulation andPraveen and Reddy (1998): the present formulation is full 3D enhanced Green-Lagrange formulation and Praveen and Reddy (1998) is a 2D FSDT of plate with PSA and von-Karman nonlinearity.

Figure 7 :
Figure 7: Temporal evolution of center deflection of clamped FGM square plates under suddenly applied uniform load.

Figure 9 :
Figure 9: Temporal evolution of center deflection of clamped FGM circular plate under suddenly applied uniform load.

5. 5
Simply Supported FGM Plate Under Varying Distributed Load This numerical example consist on a simply supported square functionally graded plate of side length a=b=1m and thickness h = 0.01 m.The FGM plate is subjected to a uniformly distributed excited transverse load ( FGM plate is made of the constituent materials: aluminum (bottom surface) and alumina (top surface).Young's modulus, Poisson's ratio and density for the aluminum are 70 The used law distribution of the FGM plate is the classical power-law which can be obtained by the general four-parameter powerlaw distribution Eq.(3): numerical problem is investigated in Duc and Cong (2015) with a sigmoid-FGM material.By noting the symmetry, only one quadrant of the plate is modeled with 8 x 8 x 1 elements.A time step of 10 -3 s is used in this numerical test.Fig. 10.a shows the nonlinear dynamic response, represented by graph of maximum deflection (p=1, N=1500N/m², Ω= 500).In the same conditions Fig. 10.b shows the relation of maximum deflection and the velocity of maximum deflection.

Figure 10 :
Figure 10: Nonlinear dynamic response of the square FGM plate (a) maximum deflection versus time (b) velocity versus maximum deflection.

Figure 11 :
Figure 11: Nonlinear dynamic response of the square FGM plate under different loads.

Fig. 11
Fig.11 shows the nonlinear response of the FGM square plate with different intensity of loads: N=1500N/m² and N=2500N/m²; Ω= 500, p=1.From obtained results, one can see that amplitudes increase, when increasing the intensity of loads.The non-linear transient responses perform the phenomenon like periodic cycles.

Figure 12 :
Figure 12: Nonlinear dynamic response of the square FGM plate under different power-law index.

Figure 13 :
Figure 13: Effect of power law index p on the responses of FGM cylindrical shell.

Figure 14 :
Figure 14: Effect of ratio R/h on the responses of FGM cylindrical shell.

Journal of Solids and Structures 14 (2017) 72-91 Figure 8:
A quadrant of the clamped FGM circular plate modeled by 192 elements, 434 nodes.
Latin American