EFFICIENT NUMERICAL MODELLING OF FUNCTIONALLY GRADED SHELL MECHANICAL BEHAVIOR

Numerical analysis of the static bending and free vibration mechanical behavior of FGM are performed using the UMAT-USDFLD subroutines in ABAQUS software. Different combinations of geometries, mechanical loading and boundary conditions are adopted. The material properties according to the coordinates of the integration points are defined in the developed numerical model. The First Order Deformation Theory is used for thin and moderately thick FG shells analysis. The accuracy and the robustness of the numerical model are illustrated through the solution of several non trivial structure problems. The proposed numerical procedure is significantly efficient from the computational point of view.


INTRODUCTION
Functionally graded material, known as FGMs, are non-homogeneous composite materials with mechanical properties that vary from one to three directions. They are generally composed of transition alloys from metal at the first surface to ceramic at the other opposite surface (Yanga & Shen, 2003;GhannadPour & Alinia, 2006;Draiche, Derras, Kaci & Tounsi, 2013). FGMs exhibit a smooth and continuous gradient in the composition and the material properties, which is the main difference with conventional composites which depicts a gradually properties variation with location. According to (Apetre, Sanka & Ambur, 2006), the closest similarities of FGMs are the laminated composites. However, their properties change abruptly across the interfaces. The improved properties of composite materials have led to their use in aerospace, automotive and biomedical applications. For this, new methodologies have to be developed to characterize their mechanical response when subjected to static and dynamic solicitations. These methodologies would be incorporated into available techniques in an optimized way.
Numerous research works are available in literature to analyze the linear mechanical behavior of FG shell structures (Simo, Fox & Rifai, 1989;Chung & Chen, 2007;Ghanned & Nejad, 2013;Shariyat & Alipour, 2014). The First Order Shear Deformation theory (FSDT) is widely used (Praveen & Reddy, 1997;Lin & Xiang, 2014;Thai & Kim, 2015). The use of FSDT is justified by the fact that the classical Kirchhoff theory (CPT) neglects the effects of transverse shear and normal strains of the structure (Wali, Hajlaoui, & Dammak, 2014). The Reissner-Mindlin theory provides a correct overall assessment. Nevertheless, equations of motion are more complicated to obtain (Frikha, Wali, Hajlaoui & Dammak, 2016a). Analytical works were conducted by Abrate (Abrate, 2006) on the free vibration and static deflections of FG square, circular, and skew plates with different combinations of boundary conditions on the basis of the FSDT.
In the objective to provide the numerical solution of FGM deformation and the effect of material inhomogeneity, the theoretical formulation is typically coupled with Finite Element method. In most studies, the solution procedure can be implemented into home codes (Abrate, 2006;Frikha, Wali, Hajlaoui & Dammak, 2016b;Frikha & Dammak, 2017); other authors obtained the numerical solution using the commercial FE ABAQUS commercial software (Nie & Zhong, 2007;Alipour & Shariyat, 2012). In this regard, variations of the material properties in the transverse direction are often modeled in ABAQUS through the following technique: each layer of the composite is divided into N slices to approximate the gradual material properties variations (Mao, Fu & Fang, 2013). This method is adopted as no such type of FGM element is available in the software element library. The main drawback of this method is non-continuous segmented distribution of material properties. This conventional numerical technique is also expensive in CPU time. Results of (Mao, Fu & Fang, 2013) showed an imperfect correlation between analytical model and numerical computation.
For robustness, an alternative method consists in using the user subroutine UMAT to define the material properties according to the coordinates of the integration points when considering one layer. To the best knowledge of the authors, there are no further accessible documents in literature on ABAQUS implementation of static and dynamic response of FG shells when taking into account the continuity of material point distribution. It is within this framework that this work is performed. The material properties of FGM shell are defined according to the coordinates of the integration points using both UMAT and USDFLD subroutines, in ABAQUS software. The accuracy of the developed model is illustrated through the solution of several FG structures problems.

BASIC CONCEPTS AND NUMERICAL IMPLEMENTATION 11 points break
A polynomial material law is adopted to control the heterogeneity of FGM properties as given by Zenkour (Zenkour, 2006) where PFGM(z), Pm, and Pc denote, respectively, the effective material property, the metal properties and the ceramic properties. EFGM(z) designs the Young modulus and ρFGM(z) is the density, h denotes the structure thickness and z is the coordinate measured along the thickness direction; n is the power-law index. , , φx and φy are the rotations of the transverse normal about the Cartesian axis x and y, respectively; u, v and w are the in-plane displacements and deflection of the mid-plane, respectively. The generalized displacement vector u is then u = [u,v,w,φx,φy] T . z is the thickness coordinate of the shell. The state of deformation can be decomposed in in-plane and transverse shear strains as: x y z e x y z x y where: εαβ, χαβ, χα3 are the membrane, bending and transverse shear strains, respectively. The strain vectors are presented in matrix notation: 11  11  13  22  22  23  12  12 , , The in-plane membrane and bending and transverse shear stresses resultants: The constitutive equations can then be written as a function of the in-plane and the out-of-plane linear elastic matrices depending on z: The constitutive relation between the generalized stress and strain is the following when applying the shear correction factor (Frikha, & Dammak, 2017): Numerical simulation is performed using the commercial FE software ABAQUS. To avoid stress discontinuity, the material properties according to the coordinates of the integration points are implemented in two interfaces: the first interface consists on the UMAT subroutine to implement the elastic mechanical behavior along thickness using the integration point number (KSPT). The second inter face consists on the USDFLD subroutine to predefine the density field variables at a material point. In ABAQUS the integration points through the thickness of the shell are numbered successively, starting from point 1, as described in Figure 2. To obtain accurate description of FGM structure response using shell elements, the number of the through-thickness integration points (n) was carefully fixed, since a small number of integration points leads to additional error of the numerical results. Considering Simpson's approach, the point (1) is located precisely on the bottom surface of the shell. Using Gauss quadrature, point (1) is placed close to the bottom surface.

Standard patch tests
The performance of the proposed numerical simulation is assessed with standard patch tests (Wali, Hajlaoui, & Dammak, 2014). Classically, the reference structure problems are threefold: (i) Bending of a simply supported rhombic plate, (ii) pinched hemispherical shell with 18 hole, (iii) pinched cylinder with end diaphragms. In the present study, shells are modeled with the standard quadrilateral 4-nodes element with three rotational and three translational degrees of freedom per node. Results, based on the FSDT of shell elements are obtained with the addition of an automatic computation of the shear correction factors as in (Frikha & Dammak, 2017).
FGM structure properties are: (Em, Ec, νm, νc)=(70 GPa, 380 GPa, 0.3, 0.3) for the metal and ceramic components, respectively. All material and geometrical properties are given in a coherent system of units in the UMAT subroutine. The power-law index is n = 6 for all cases.
The obtained numerical results of the normalized center-point deflection are gathered in table 1. Comparison with (Wali, Hajlaoui & Dammak, 2014) show a very good accuracy. In (Wali, Hajlaoui & Dammak, 2014), a 7DOF per node, 3d-shell model was applied based on a discrete double directors shell element which is expensive from a computational time point of view. So, one can conclude that the proposed technique exhibits high performance.

Free vibration of FG shear-diaphragm cylindrical shell
In this section, the numerical model is applied for the case of FGM cylindrical shell subjected to free vibration with shear-diaphragm boundary conditions. The cylindrical shell has a thickness to radius ratio fixed to 0.002 and length to radius ratio fixed to 20. FGM structure properties are: (Em, Ec,νm, νc, ρm, ρc) = (205.098 GPa, 207.788 GPa, 0.31, 0.317756, 8900 kg.m -3 ,8166 kg.m -3 ) for nickel and stainless steel components, respectively. For free vibration, the equations of motion take the form of a standard eigenvalue problem: where w is the eigen-frequencies of the FGM shell. The discretization of the cylinder is performed by means of S4 standard structural shell elements. The S4 shell element is widely used for industrial applications. It is suitable for thin to moderately thick shell structures. The accuracy of the proposed simulation is assessed with results in (Wali, Hentati, Jarraya & Dammak, 2015). The variation of the natural frequencies (Hz) against circumferential wave number is illustrated in Table 3. The longitudinal wave number is equal to 1. It is plausible to depict the effect of the power-law distribution choice on the frequency parameters. Some mode shapes for the shear-diaphragm FGM cylindrical shells are illustrated in Fig. 3.

Free vibration of FG conical panel
In this section, results of free vibration of moderate thick FG conical panel clamped on the bottom surface are presented. Results are compared with numerical findings of Tornabene (Tornabene, 2009). The four-parameter powerlaw distribution is: where, Vc denotes the ceramic volume fraction and the parameters a, b, c dictate the material variation profile through the FG shell thickness. Both constituents of FGM are the zirconia (ceramic) and aluminum (metal). Material properties and geometry for the zirconia and aluminum are detailed as considered in (Tornabene, 2009). The first ten frequencies for the FG conical panel as a function of the powerlaw exponent are gathered in Tables 3. The first six mode shapes are plotted in Figure 4. Results shows a good correlation with literature. The maximum relative error is about 1.2%. It is interesting to note that frequencies reach a minimum value for a shell made only of metal, due to the fact that aluminum has a much smaller modulus than zirconia. In fact, when increasing the power-law index, the frequencies decrease, until tending to the metal limit case. This is plausible as the ceramic content decreases by increasing the index n and the FGM shell approaches the case of the fully metal shell. This is in accordance with (Tornabene, 2009).
Mode shape 1 Mode shape 2 Mode shape 3 Mode shape 4 Mode shape 5 Mode shape 6

CONCLUSIONS 11 points break
In this research, a numerical FE model is implemented using the UMAT-USDFLD subroutines, in ABAQUS software in order to predict the FGMs material response to static and free vibration solicitations. The material properties according to the coordinates of the integration points are defined so that stress discontinuity at the interfaces are eliminated. To the best knowledge of the authors, there are no further accessible documents in literature on ABAQUS implementation of FG shells mechanical behavior when taking into account the continuity of material point distribution. This is the main contribution of the present research.
To check the performance of this technique, the present results are compared and validated with findings available in literature. It is also verified that the proposed solution procedure is significantly efficient from the computational point of view. The present study enables closed-form solutions for some fundamental solid mechanics problems, and will aid the development of finite element models for structures made of FGMs. Indeed, the present work can be applied in case of complex shells structures, contrarily to the analytical formulation which is limited to weak differentiable geometry of shell structures.