Nonlinear Static Bending and Free Vibration Analysis of Bidirectional Functionally Graded Material Plates

Bidirectional functionally graded material (2D-FGM) plates have mechanical properties that vary continuously in both the thickness and one-edge directions; these plates are more and more widely used in design and engineering applications. When these structures are subjected to strong loads, they can be largely deformed; therefore, nonlinear calculations, in this case, are necessary. In this paper, nonlinear static bending and nonlinear free vibration behaviors of 2D-FGM plates are studied by using the finite element method based on the third-order shear deformation theory; the Newton-Raphson method is used to solve this problem. The accuracy of this approach is confirmed by comparing the results with respect to other papers. The effects of some numerical aspect ratios such as volume fraction index and thickness-to-length ratio on nonlinear static bending and free vibration of the plates are explored. This study shows that there is a big difference between the numerical results obtained from the nonlinear problem and those from the linear one.


Introduction
Functionally graded materials are a new smart type of composites that were introduced for the first time by Japanese researchers in 1984. The main difference between the new material and the classical composite layer is that the mechanical properties change continuously from one interface to the other. Thus, the stresses distribute smoothly in some directions, and the FGM structures become stronger as they mechanically deform. Nowadays, FGMs are applied widely in many engineering fields. For 1D-FGM structures (usually called FGM structures), they are made up of two different materials (ceramic and metal), in which the material properties vary smoothly from one surface to the other one by the thickness direction. For 2D-FGM structures, their formulations are based on three material components to take advantage of all three materials, where the mechanical properties change in two directions, which are the thickness direction and the longitudinal direction. Therefore, considering the mechanical behaviors of these types of structures is really important. It is very helpful in designing, manufacturing, and using them in engineering practice.
There are many works that have investigated the mechanical behaviors of one-direction [1][2][3][4][5][6] and twodirection FGM plates [7][8][9][10][11][12][13][14][15][16][17] (called FGM and 2D-FGM plates); the 2D-FGM plates with the volume fraction exponents are the power-law distributions in the longitudinal and thickness directions. Scientists also paid much attention to nonlinear static bending and free vibration analysis of FGM plates based on several plate theories. A nonpolynomial higher-order shear deformation theory with inverse hyperbolic shape function was adopted by Singh and Harsha [1] to research the free and forced nonlinear vibration characteristics of FGM plates. Ebrahimi and Rastgoo [2] investigated the nonlinear static problem of a circular FGM plate using an analytical method. Reddy [3] used the third-order shear deformation theory to solve the nonlinear static and dynamic bending problems of the FGM plate based on the finite element method and Navier's solutions. Bui et al. [4] investigated the static bending of FGM plates in the thermal environment based on the finite element solution and the third-order shear deformation theory. Lan et al. [5] explored the nonlinear bending behaviors of the FGM plates based on a four-node element within the context of the first-order shear deformation theory. Kien and his coworkers [6] researched the nonlinear static behaviors of planar beam and frame structures made of functionally graded material based on the finite element method and Bernoulli beam theory (classical beam theory). Civalek [18] used the discrete singular convolution method to study free vibration analysis of conical and cylindrical shells and annular plates made of composite laminated and functionally graded materials (FGMs). Arefi et al. [19] employed the higher-order shear deformation theory of Reddy to carry out the static bending analysis of FG graphene nanoplatelets (GNPs) reinforced composite microplates with porosity. Bendenia et al. [20] studied the static and free vibration behavior of FG-CNT reinforced sandwich plates resting on Pasternak elastic foundation based on the first shear deformation theory. In the works [21,22], the authors presented the static and dynamic responses of FGM plates by using a refined shear deformation theory. Rabhi et al. [23] investigated buckling and free vibration of exponentially graded sandwich plates resting on elastic foundations using a new innovative 3-unknown HSDT. Chikr and his coworkers [24] used Galerkin's approach and a refined trigonometric shear deformation theory to study the buckling response of FG sandwich plates resting on elastic foundations. Refrafi et al. [25] calculated the hygrothermal and mechanical buckling responses of a simply supported FG sandwich plate seated on the Winkler-Pasternak elastic foundation based on a novel shear deformation theory. Rahmani et al. [26] employed an original novel high order shear theory to examine the influence of boundary conditions on the bending and free vibration behavior of functionally graded sandwich plates resting on a two-parameter elastic foundation.
Besides, mechanical behaviors of 2D-FGM structures were also considered widely all over the world, and scientists obtained significant knowledge. Thom and his coworkers [7] studied static bending and buckling of 2D-FGM based on the finite element method and new thirdorder shear deformation theory. Nemat-Alla et al. [8] used a 3D finite element model to research the mechanical behavior of 2D-FGM plates in the thermal environment.
Simsek [9] studied free vibration and dynamic responses of the 2D-FGM Timoshenko beam based on the Lagrange equations and simple polynomials. The free vibration of the 2D-FGM Timoshenko beam was also studied by Deng and Cheng [10] using a new dynamic stiffness matrix solution. By using a semianalytical and seminumerical method, Wang et al. [11] focused on the natural frequencies of a 2D-FG Euler-Bernoulli beam. Rad [12] used a semianalytical method in order to research the static problem of bidirectional FGM auxetic porous circular plate resting on the Pasternak foundation. Esmaeilzadeh and Kadkhodayan [13] investigated dynamic responses of the bidirectional FGM porous plate structure; the solution was based on a dynamic relaxation method and a structure with moving loads. Mahinzare et al. [14] calculated the free vibration of a bidimensional functionally graded microplate by using a numerical method, while the structure rested on a Winkler-Pasternak foundation under the thermal load, in which temperature-dependent mechanical properties varied gradually in the thickness and radial direction of the plate. Lieu and Lee [15] considered taking the material optimization method into the free vibration analysis of multidirectional FGM plates based on the IGA method. The buckling and vibrational behaviors of the inplane bidirectional FGM plate were also explored by Lieu et al. [16]. Wu and Yu [17] used the combination of the prism method and Reissner's mixed variational method to investigate the free vibrational behavior of bidirectional FGM plates.
The functionally graded structures have been using widely in engineering practice; when they are subjected to strong loads, large deformations can be occupied. Therefore, linear problems cannot describe exactly the mechanical behavior of the structures. We need to use nonlinear problems to deal with these large deformations; this means that we have to consider the geometrical nonlinear factor in these cases. The explorations of mechanical responses of 2D-FGM plates taking into account nonlinear factors are still limited, especially problems related to nonlinear static bending and free vibration of structures in consideration of large deformation. These interesting issues require a higher complicated process than small deformation problems. In this paper, the nonlinear static bending and nonlinear free vibration of bidirection functionally graded material plates are investigated based on the finite element model and the third-order shear deformation theory; 4-node plate elements are used, in which each node with 5 degrees of freedom can reduce the working load much while the accuracy still remains. Numerical results of nonlinear mechanical behaviors of 2D-FGM plates are presented, then compared with linear results to show that the nonlinear results are much different from linear results.
The rest of this paper is structured as follows. Section 2 presents shortly bidirectional functionally graded plates used in this work. Finite element formulations for nonlinear static bending and free vibration analysis are introduced detail in Section 3. Section 4 presents the verification problems, numerical results, and discussions. Novel explorations are summed up in Section 5.

Bidirectional Functionally Graded Plates
A plate made of three kinds of materials with length a, width b, and thickness h is considered as shown in Figure 1.
Young's modulus, Poisson's ratio, and mass density change continuously in the longitudinal and thickness directions, and these properties can be calculated as [7] where E i , ν i , ρ i , and V i (i = 1 -3) are Young's modulus, Poisson's ratio, mass density, and the volume fraction of i-th material, respectively; x and y are lines of inplane of 2 International Journal of Aerospace Engineering midsurface; and z defines the normal axis. The volume fractions of materials in this paper are given by [7] being everywhere V 1 + V 2 + V 3 = 1. The distributions of V 1 , V 2 , and V 3 through the 2D-FGM plate when n = 0:5 and q = 2 are also plotted in Figure 2. One can see that the volume fraction indexes of three materials change smoothly in both thickness and longitudinal directions. Therefore, the mechanical properties of the materials also vary gradually in the corresponding directions. This is the main difference from 1D-FGM structures, in which the material only changes in one thickness direction.

Finite Element Formulations for Nonlinear Bending of 2D-FGM Plates
Based on the Reddy third-order theory, the displacements at a point (x, y, z) in the plate can be written as [3] where u 0 ,v 0 , and w 0 are the displacements at a point (x, y, z = 0); φ x and φ y are the transverse normal rotations in xzand yz-surfaces.
The strains can be written clearly as According to Hooke's law, the relationship between stresses and strains is written by  International Journal of Aerospace Engineering Thus, the stress-deformation relation depends on both x and z coordinates; this is different from the case of 1D-FGM; therefore, this also makes the calculation process more complicated. The four-node element with five degrees of freedom in each node is used in this study, which is expressed as where N j is the Lagrangian interpolation function of the four-node element and ξ, η are natural coordinates, so that the strains relate with element displacement vector (q e ) as ,  2   6  6  6  6  6  6  6  6  6  4   3   7  7  7  7  7  7  7  7  7  5 ; After integrating from −h/2 to h/2, the normal force, bending moment, high order moment, shear force, and high order shear force can be expressed as where The total strain energy of the 2D-FGM plate can be given 5 International Journal of Aerospace Engineering whereQ is the surface loading, and Equation (13) is expressed in the matrix form as in which the linear matrix is calculated as follows: and the nonlinear matrix is expressed as follows: where F e is the element force vector. Hence, the governing equation for nonlinear bending analysis becomes The kinetic energy of a plate element is determined by the following expression: where N is the shape function and matrix L and the element mass matrix M e are expressed as follows: , For free vibration analysis, in order to find fundamental frequencies, we need to solve the following equation: Note that Equations (17) and (20) contain the element stiffness matrix, which depends on the displacement vector q e ; therefore, to solve this equation, the Newton-Raphson method is used. Equation (17) can be written in short form as follows: where KðqÞ is the stiffness matrix, which depends on nodal displacement q; q = ∑ e q e is an unknown displacement vector. Assuming that the value of the displacement is found in step i, denoted by q i , to find displacements in the next step (q i+1 ), the Taylor expansion to first-order is used as follows: in which K i T ðq i Þ = ð∂K/∂qÞ i is a Jacobian matrix, usually called the tangent stiffness matrix, andΔq i is a displacement increment. From Equation (22), the displacement increment Δq i can be found in the following equation: Now, Equation (23) will give displacement increment Δq i ; the displacement in the next step q i+1 is calculated as follows: The solution from Equation (24) does not satisfy the exact roots of nonlinear Equation (21). This time, there will exist an unbalanced force: This process is repeated until the unbalanced force is less than the given value; q i+1 will be accepted as solutions to nonlinear Equation (21). To determine this stationary time, a comparison of the convergence parameter is 6 International Journal of Aerospace Engineering carried out ε conv ≤ ξ 0 , in which ε conv is defined as and ξ 0 = 10 −4 can be acceptable.
To solve the nonlinear eigenvalue problems (20), an iterative procedure is used. Firstly, neglecting the nonlinear stiffness matrix component of Equation (20), the linear fundamental frequency is calculated, and then it is normalized. Next, the normalized vector is amplified/scaled up so that the maximum displacement is equal to the desired amplitude, supposing that w/h = 0:6 (w is the maximum lateral displacement, h is the thickness of the plate). This gives the initial vector, denoted by ζ. The iterative solution procedure for the nonlinear analysis starts with the initial vector, ζ. Based on this initial mode shape ( ζ), the nonlinear stiffness matrix, which depends on displacement, is formed, and subsequently, the updated eigenvalue and its corresponding eigenvector are obtained. This eigenvector is further normalized and scaled up by the same amplitude (w = h), and the iterative procedure adopted here continues till the frequency values and mode shapes evaluated from the subsequent two iterations satisfy the prescribed convergence criteria [27] as in which p, i, N, and r represent the mode number, degree of freedom of the finite element model, total degree of freedom, and iteration number, respectively.  [3] are shown in Figure 3; we can see that the results obtained by this approach are in a good agreement with the reference results of Reddy [3], which used the mesh with 8 × 8 elements.

Numerical Results of Nonlinear Static
Bending for 2D-FGM Plate. In this section, the nonlinear static mechanical displacements of 2D-FGM plates are analyzed using the proposed method. The fully simply supported plate is subjected to a uniform load P 0 and made from three materials with properties [7] E 1 = 151 GPa, E 2 = 205:1 GPa, E 3 = 70 GPa, and ν 1 = ν 2 = ν 3 = 0:3. In this analysis, the deflection and load are normalized by the following formulas: (1) Effect of Volume Fraction. To study the effect of the volume fraction on the nonlinear static behavior of 2D-FGM (a/b = 1, a/h = 100), the values of n and q are changed from 0 to 10, maximum deflections w * of this plate depending on n and q are plotted in Figure 4. The deflection changes differently with the increase in n and q when n increases and deflection increases; however, it decreases when q increases. The linear and nonlinear maximum deflection in two values of q (0.5 and 10) are shown in Figures 5 and 6; these figures show that when increasing the applied load, the nondimensional deflection in the case of the linear problem changes much than that in the case of the nonlinear one, and the nonlinear curve of w * changes sharply when P * increases in a range of 0 to 40; for the other range of 40-100, the curve of w * changes gradually. And one can see that when the load applied to the structure has small intensity (P * <5), the numerical results of the linear problem are very close to those of the nonlinear problem. However, when P * >5, the results between the linear problem are much different from those of the nonlinear problem. This demonstrates that nonlinear results are much different from linear results when the load applied to the structure has great intensity, so that nonlinear analysis is more useful than the linear one. (2) Effect of the Thickness-to-Length Ratio a/h. In order to investigate the effect of the thickness-to-length ratio a/h on the nonlinear mechanical bending response of 2D-FGM plate (a/b = 1), new normalization for this specified analysis is used as follows: w * * = W max /h 0 , h 0 = a/100, P * = P 0 a 4 0 /E 3 h 4 0 = 50. Thicknesses are considered as different values of the thickness-to-length ratio a/h = 40-100. The numerical results of maximum deflections in two cases of q are plotted in Figure 6. One can see that by reducing the plate thickness, the higher nonlinear mechanical deflections of 2D-FGM plates are obtained. It is also explained that when the plate thickness is reduced, its stiffness will decrease, and the displacements will increase; therefore, it can be understood that the ratio a/h has a great influence on the nonlinear bending responses of this plate. At the same time, the change of w * is sharper when n increases from 0 to 2, and w * varies gradually when n > 2.   Figure 7. It can be seen that the higher the aspect ratio b/a is taken, the higher the nondimensional deflections of 2D-FGM plates are obtained. The reason is that when b/a increases, the plate becomes softer; thus, the nondimensional deflection also increases. The material properties, as given in Shen [30], are E m = 105 GPa, ν m = 0:3, ρ m = 4429 kg/m 3 for Ti-6Al-4V and E c = 168 GPa, ν c = 0:3, and ρ c = 3000 kg/m 3 for ZrO 2 . The nonlinear to linear frequency ratios ω NL /ω L presented in Table 2 are com-pared with Shen [30] (exact results) and Behjat and Khoshravan [31] (finite element model with mesh of 6 × 6 elements); good agreement can be seen from this table when the number of elements is large enough. Besides, it also can be observed that the mesh 14 × 14 ensures the accuracy; therefore, for all succeeding investigations, this mesh will be employed.    Tables 3, 4, and 5 and Figures 8, 9, and 10. It can be seen that for the case of q = 0:5 and 1, when n increases, the ratio ω NL /ω L also increases and obtains the peak point around the value of n = 1, then it goes down. For the case of q = 10, when n increases, the ratio ω NL /ω L increases, too. Thus, depending on the specific case of material distribution, the case that corresponds to the maximum value of the ratio ω NL /ω L can be found. In addition, when the ratio w c /h = 0:2, then the difference between ω NL and ω L is small (maximum 3.05% for n = 1). However, when w c /h increases, the difference between ω NL and ω L also increases; the largest one can be up to 44.84%. This also shows that when the deformation is small, the effect of the nonlinear factor is not much. The larger the deformation is, the larger the component K N e in Equations (16) and (17) gets. The stiffness of the structure is increased, and as a result, the nonlinear specific frequency and the difference between ω NL and ω L are also increased. Figure 11 shows the first five fundamental vibration mode shapes of the 2D-FGM plate in the case of n = 10, q = 1, and w c /h = 0:8.
(2) Effect of the Plate Thickness. Now in the last subsection, in order to explore the effect of the plate thickness on the free vibration of the 2D-FGM plate, the thickness changes in a range of a/10 to a/100, where other parameters are n = 0-10, q = 0:5, and w c /h = 0:8. The ratio ω NL /ω L is shown in Figure 12. Table 3: Nonlinear-to-linear frequency ratio ω NL /ω L of 2D-FGM square plate, q = 0:5, a/h = 10.     It can be obtained that the thinner the plate thickness is, the higher the ratio ω NL /ω L gets. This phenomenon means that when the plate becomes thinner, the difference between the nonlinear frequency and the linear frequency is higher. Furthermore, when increasing the value of the volume fraction index n, the nonlinear frequency will reach the maxi-mum value then drop down gradually. This means that there is one value of n (around n = 1) in order that the nonlinear frequency obtains the maximum value. It also means that the difference between ω NL and ω L reaches its maximum value when n is around 1; the difference between ω NL and ω L can be as close as 33%.

Conclusions
Nonlinear static bending and free vibration responses of 2D-FGM plates are investigated. The plate is made from three materials, which change continuously in longitudinal and thickness directions. The study is based on the finite element method and third-order shear deformation theory. The results obtained by this method are compared with those of other methods to show the accuracy of the proposed theory and mathematical model. The paper explored the effect of geometric and material parameters on nonlinear static behaviors of the 2D-FGM plates. It shows that there is much difference between linear and nonlinear deflections. When the plate has a large deformation, some highlighted results in specific problems are as follows: (i) For the nonlinear static bending: when the load applied to the structure has small intensity (P * <5), the numerical results of the linear problem are very close to those of the nonlinear problem. However, when P * >5, the maximum deflection w max between the linear problem is much different from that of the nonlinear problem (ii) For the nonlinear vibration: when the plate has a small deformation (w c /h = 0:2), the influence of strain w c on the nonlinear stiffness matrix composition is also small. Therefore, the difference between the nonlinear oscillation frequency and the linear oscillation frequency is also small (3.05%). As w c /h increases, the difference between the nonlinear oscillation  14 International Journal of Aerospace Engineering frequency and the linear one increases as well, and the maximum can reach up to 44.84% when w c /h = 1 These results will give good information for the design and application of these structures in practice.

Data Availability
Data used to support the findings of this study are included within the article.

Conflicts of Interest
The author declares that there is no conflict of interest regarding the publication of this paper.