Formulation of a New Mixed Four-Node Quadrilateral Element for Static Bending Analysis of Variable Thickness Functionally Graded Material Plates

A new mixed four-node quadrilateral element (MiQ4) is established in this paper to investigate functionally graded material (FGM) plates with variable thickness. The proposed element is developed based on the ﬁrst-order shear deformation and mixed ﬁnite element technique, so the new element does not need any selective or reduced numerical integration. Numerous basic tests have been carried out to demonstrate the accuracy and convergence of the proposed element. Besides, the numerical examples show that the present element is free of shear locking and is insensitive to the mesh distortion, especially for the case of very thin plates. The present element can be applied to analyze plates with arbitrary geometries; it leads to reducing the computation cost. Several parameter studies are performed to show the roles of some parameters such as the power-law index, side-to-thickness ratio, boundary conditions (BCs), and variation of the plate thickness on the static bending behavior of the FGM plates.


Introduction
In the last decade of the 20th century and more than 20 years of the 21st century, functionally graded materials (FGMs) were extensively utilized in many fields of engineering as well as industry because of their advantages [1][2][3][4]. FGMs are usually produced by mixing two or more different ingredients with the volume varying continuously in the spatial direction. e most common type of FGM is ceramic-metal FGM, which is fabricated by ceramic and metal since the ceramic is compatible for working in high-temperature environments and the metal is suitable for working with the mechanical load.
at is the reason why numerous researchers focus on investigating the dynamic as well as static behaviors of the FGM structures, for example, beams [5][6][7][8][9][10], plates [11][12][13][14][15][16][17][18], shells [19][20][21], and other structures [22,23]. Many plate theories have been established and/or applied to investigate these structures, for example, classical plate theory (CPT), first-order shear deformation (FSDT), higherorder shear deformation theory (HSDT), and quasi-3D shear deformation theory. ai and Choi [24] considered the static bending and free vibration of FGM plates using a simple FSDT. A simple HSDT was developed and used to examine FGM plates by ai and Kim [25]. Zenkour [26] established a new shear deformation theory called generalized shear deformation theory (GSDT) to analyze the static bending of FGM plates. e influence of porosity on the bending of the FG plate had been studied by Hadj et al. [27] via refined plate theory. Although the FSDTs and HSDTs and their variations are suitable to predict the behavior of thin to moderate plates, it is necessary to develop quasi-3D theories to analyze thick and very thick plates. ese theories take into account the normal transverse stress which is important in these types of FGM structures. Neves et al. [28] established a quasi-3D sinusoidal shear deformation theory to study the bending and free vibration of FGM plates. ai and Kim [29] constructed a simple quasi-3D sinusoidal shear deformation theory to scrutinize FGM plates. More details on the development of the plate theory can read from these journals and their relevant references.
Although the analytical method can give exact solutions to many engineering problems, it is remarkable that it cannot be applied in complex problems such as analysis of complicated structures, arbitrary geometric plates, and variable thickness plates. erefore, large numerical methods have been developed to analyze these structures. Besides, the use of the FGM structures with variable thickness is becoming the new trend in advanced structures to reduce their weight and improve the art of these constructions [30,31], for example, the cover of the blades of the wind turbines, the integral shroud blades, and the cover of the spacecraft. ese structures usually have complex geometric and boundary conditions (BCs). In these cases, the application of an analytical method to analyze the complex geometric plate is very difficult. Some worthy numerical methods can be considered herein, including differential quadrature method (DQM), generalized differential quadrature method (GDQM), finite element method (FEM), smoothed finite element method (SFEM), extended finite element method (XFEM), meshless method (MM), and isogeometric analysis (IGA). Among them, FEM is widely employed in many complicated engineering problems because of its simplification, good accuracy, and high convergence rate. It also noticed that the development of plate elements based on the Mindlin-Reissner plate theory or the FSDT has been studied for a long time and achieved a lot of success. However, at the beginning of this development, the most and major difficulty of these types of plate elements was the shear-locking phenomenon in the case of thin to very thin plates. Numerous techniques and schemes have been introduced to address this challenge. Zienkiewicz et al. [32] introduced the reduced integration technique and Hughes et al. [33] presented selective reduced integration to analyze thin to very thin plates and shells using classical four-node quadrilateral plate element (Q4). Hinton and Huang [34] developed many quadrilateral plate elements using substitute shear strain fields. Bathe and Dvorkin [35] created some plate elements by introducing the mixed interpolated tensorial components (MITC) technique, which is widely used for the analysis of thin plates and shells. A discrete shear method was developed by Batoz and Lardeur [36] for both quadrilateral and triangular plate elements. Zienkiewicz et al. [37] developed the linked interpolation technique for the Mindlin-Reissner quadrilateral plate element; then it was extended to the higher-order linked interpolation quadrilateral thick plate element by Ribaric et al. [38]. Soh et al. [39,40] developed an improved shear strain interpolation technique based on Timoshenko's beam formulae. Castellazzi and Krysl [41] constructed a new displacement-based element by using the nodal integration method and assumed strain technique. Hansbo et al. [42,43] formulated two locking-free quadrilateral elements based on continuous/ discontinuous rotations. e mixed shear projected approach element (MiSP) was established by Ayad and his colleagues [44,45]. e combination of the discrete shear gap technique and alternative alpha FEM was used by Nguyen-anh et al. [46] to analyze isotropic plates. Liu et al. [47,48] developed the smoothed FEM, which integrates the strain smoothing technique into the standard FEM to construct many effective elements, such as cell-based SFEM [49][50][51], edge-based SFEM [52], and node-based SFEM [53]. e essential success of these elements is that they predict more accurate results, higher convergence rates, and insensitive to distorted mesh in comparison to conventional FEM. e generalized conforming elements were first developed by Long et al. [54][55][56][57] to analyze membrane, plate, and shell structures. Cen et al. [58] initiated a new element based on the quadrilateral area coordinate method for Mindlin-Reissner plates. Cen et al. [59] developed a hybrid-Trefftz stress element based on the hybrid displacement function (HDF) to analyze Mindlin-Reissner plates. A family of hybrid-Trefftz p-elements were established by Jirousek et al. [60] to analyze thick plates. Katili [61] used the discrete Kirchhoff-Mindlin method to establish a new plate element DKMQ to analyze thin to moderate plates. Among them, the classical Q4 element, MITC4 (four-node MITC element) element, and DKMQ element have been widely used to analyze composite and FGM plates, in which the classical Q4 element needs a reduced or selective integration and is very sensitive to the mesh distortion. e MITC4 is an efficient element with thick and moderate plates, but less accuracy with very thin plates.
e DKMQ element has been developed since 1993 and it is compatible to analyze thick to thin and very thin plate, but the convergence rate is slower than MITC4 [61].
Although a lot of plate elements were developed and applied successfully for the analysis of thin to moderate plates, there are some significant challenges that remain outstanding, such as the realizable results of the element with the distorted mesh, avoiding shear locking especially in very thin plates, and good precision for the stress of the element. is paper aims to establish a new mixed fournode quadrilateral plate element based on the mixed finite element formulation and the FSDT.
is is the main contribution of the current work. e displacement-based four-node element will be first used as the initiating point; then the formula of the mixed finite element will be applied to construct the proposed plate element named MiQ4 (mixed four-node quadrilateral element). Several numerical patch tests and verification studies will be considered to demonstrate the accuracy and efficiency of the proposed element MiQ4. e MiQ4 is a simple modification of the classical element Q4 so that the present element is simple in its formulation. In addition, the present element MiQ4 is free of shear locking without using reduced integration. Another advantage of the proposed element is that it is insensitive to mesh distortions, so it can be employed to analyze the FGM plates with complex geometry. en the proposed element will be employed to research the bending behavior of the variable thickness FGM plates. A comprehensive parameter study will be carried out to demonstrate the effects of some parameters on the bending behavior of the variable thickness FGM plates.

Modeling of Variable Thickness FGM Plate
In this study, a variable thickness FGM plate with the dimension of a × b and thickness of h as shown in Figure 1 is considered. e properties of the material vary continuously 2 Mathematical Problems in Engineering in the thickness direction of the plates by a power-law function as the following formula [1,14]: where P c and P m are, respectively, Young's modulus E or density ρ of ceramic and metal and p is the power-law index.
In this study, the thickness of the plates h depends on the variable x only h � h(x), and the four following types of thickness variation, which are uniform plate (U plate), linear plate (L plate), parabolic plate (P plate), and sinusoidal plate (S plate), are considered: Parabolic (P plate), where ζ is the varied parameter and h 0 is the thickness of the plate at x � 0, y � 0. In this study, the varied parameter ζ changes in the range of 0 ≤ ζ ≤ 0.5. It is noticed that the midplane of the variable thickness plate is still planar; it means that the plate is symmetric with respect to the Oxy plane.

Governing Equations
e displacement field of the plate based on FSDT is given as follows: u(x, y, z) � u(x, y) + zψ x (x, y), v(x, y, z) � v(x, y) + zψ y (x, y), w(x, y, z) � w(x, y), where u, v, w, ψ x , ψ y are five unknown displacement functions at the middle surface of the plate. e strain field of the plates is obtained as the following formulae: (4) e short form of equation (4) is where e vector of the transverse shear strain can be rewritten as e stress field of the plate is obtained as follows: Here, where E is Young's modulus of the material and k is the shear correction factor. e strain energy of the plate and the work done by external force can be calculated as the following formulae: Hamilton's principle is used to obtain the governing equation of motion of the plate Mathematical Problems in Engineering 3 e variation of the strain energy is Equation (12) can be rewritten in the following compact form: By integrating equation (13) across the thickness, one gets where P and Q are given by and δΩ � δε 0 Substituting equation (17) into equations (15) and (16), one gets P � HΩ, where

Displacement Finite Element Formulation.
In this subsection, a brief review of a classical four-node quadrilateral element based on displacements with five degrees of freedom (DOFs) is considered as the starting point. e nodal displacement vector of the i-th node is e nodal displacement vector of the plate element, U, is defined as en the displacement variables at any point of the element are approximated as the following formulae: where N i is the bilinear shape functions of the four-node quadrilateral plate element. e expressions of the shape functions in the local coordinate are Equation (17) can be rewritten in short form as the following formula: Here, e variation of the shear strain vector is obtained as follows: Equation (28) can be written in short form as where with By substituting equations (25) and (29) into the expression of the variation of the strain energy, we get e variation of the work done by an external force is calculated by Here, N w is the matrix of the shape functions related to the transverse displacement; its expression is as the following formula: where Inserting equations (32) and (33) into equation (11) and using the trivial manner of classical FEM, one gets Here, the stiffness matrix and the nodal force vector of the plate element are computed by the following formulae: where e Gauss integration is used to estimate the stiffness matrix K b , K s as well as the nodal force vector f. For each point of Gauss integration, the integration through the thickness direction of equations (19) and (20) will be evaluated via Newton-Cotes quadrature rules. It is noticed that K b and f are calculated by using full Gauss integration (2 × 2), while the shear stiffness matrix K s is estimated by using reduced Gauss integration (1 × 1) to avoid the wellknown shear-locking phenomenon. In this study, the acronym Q4 is used to denote the classical four-node quadrilateral plate element based on the displacement formula with reduced Gauss integration.

Mixed Finite Element Formulation.
In this subsection, we introduce a quadratic interpolation for ψ x , ψ y , and then the expressions of ψ x , ψ y are obtained as follows: Substituting equation (39) into equation (25), the expression of δε 1 can be obtained as follows: e expression of the strain vector is obtained as follows: where e shear strain δc is expressed as the following formula: We rewrite equation (46) in the matrix form as follows: where Mathematical Problems in Engineering We introduce a constant shear strain field c � H − 1 s Q 0 , and following the procedure of the mixed FEM, the expression of the potential energy is evaluated as [62] where Π bc is the energy of the constrain to consider the BC effects. e variations of equation (49) are Inserting equations (43) and (47) into equations (50) and (51), one gets Inserting c � H − 1 s Q 0 into equations (52) and (53), one gets By combining equations (53) and (57), one gets where 6 Mathematical Problems in Engineering where |J| is the determinant of the Jacobian matrix. By eliminating ψ m and Q 0 from equation (58) at the element level, one gets in which where For convenience, the acronym MiQ4 is used to denote the proposed mixed four-node quadrilateral plate element. It is noticed that all integrations in equations (59)-(65) relate to the element stiffness matrix, the nodal force vector, and they are calculated using full Gauss integration, but the shear-locking phenomenon does not appear as we will appreciate and discuss in the next section. Besides, the proposed plate element MiQ4 is insensitive to mesh distortion, so it can be employed to analyze the plates without worry about the complicated geometries. Moreover, the shape functions of MiQ4 are bilinear and the number of DOFs of MiQ4 is equal to that of Q4, so the computation cost is reduced. ey are substantial benefits of the proposed plate element in comparison with the classical four-node quadrilateral element Q4 and other four-node quadrilateral elements.  [60]. e errors of the numerical results using MiQ4 and Q4 and reference results with different meshes and BCs are displayed in Figure 3. According to Tables 1-3 and Figure 3, it can be seen that the proposed element MiQ4 is an excellent convergence speed for the case of fully simple supported (SSSS) plates and good convergence speed for the case of fully clamped (CCCC) ones. In comparison with Q4, although the MiQ4 element is using full Gauss integration for all its integral formulae of stiffness matrix and nodal force vector, the convergence rate of MiQ4 is better than the Q4 element, especially for coarse meshes and very thin plates. Besides, the proposed element MiQ4 seems stiffer than the classical element Q4. In the case of SSSS plates, the convergence of the MiQ4 element is better than that of the Q4 element for thick to thin plates. In the case of CCCC plates, the convergence of the MiQ4 element is worse than that of the Q4 element; this is a disadvantage of the proposed element. However, the element MiQ4 is much insensitive with mesh distortion which will be demonstrated in the following subsections.

Numerical Results and Discussions
Additionally, a benchmark test about the shear locking free of the proposed plate element MiQ4 is presented. e geometrical and material properties of the plate are retained while the thickness-to-side ratio changes from the thin case (h/a � 10 − 3 ) to the very thin case (h/a � 10 − 30 ). e numerical results of the classical Q4 element with reduced integration, the proposed element MiQ4 with full integration, HDF element HDF − P4 − 11β of Cen et al. [59], and analytical solution [59] are given in Table 4. It is obvious that the proposed element is free of shear locking even in very thin plates, while the classical Q4 element cannot work     Figure 4) subjected to a uniform load which is carried out by Morley [63]. e skew angle of the plate is β � 30 o and the BCs of the plate are soft simply supported (w � 0) at all edges. e geometry and material properties of the plate are a � 100, E � 10.92, and ] � 0.3 and the side-tothickness ratio a/h is 100 and 1000. e comparisons of the nondimensional central deflection (w * � w c 10 3 D/qa 4 ) of the acute skew plates using the proposed element MiQ4, the numerical results using DKMQ [61], and those of MITC4 [61] element are presented in Tables 5 and 6 and Figure 5, in which the error of the numerical results is calculated in comparison with the solutions of Morley [63]. e comparison shows that the results of MiQ4 are in good agreement with the reference solution of Morley [63]. According to Tables 5 and 6, it can be seen that the convergence of the proposed element MiQ4 element for both thin and very thin plate is similar. e convergence of the proposed element MiQ4 is better than that of DKMQ and equable with MITC4 for thin plates (a/h � 100). For the case of very thin plates (a/h � 1000), the convergence of the present element MiQ4 is better than that of both DKMQ and MITC4 elements. Hence, the MiQ4 element can be applied to analyze both thin and very thin plates.

Square Plates with Severe Mesh Distortion.
To exhibit the advantage of the MiQ4 element in comparison with the classical Q4 elements, a fully simply supported square plate is investigated with four types of distorted mesh. e      [63] 0.408 dimension and material properties of the plate are a � 10, E � 10.92, and ] � 0.3 and the thickness-to-side ratio h/a is 0.1, 0.01, and 0.001. Four types of distorted mesh are displaced in Figure 6 for the case of 64 elements. e nondimensional central deflection of the plates using MiQ4 and Q4 and the reference solutions [60] are presented in Tables 7-10 and Figure 7, where the error is calculated in comparison with the reference solutions [60] and those of Cen [59]. In comparison with Q4, MiQ4 is survival for all cases of distorted mesh, especially in the case of very thin plates, while the classical element Q4 cannot work accurately. is is a significant benefit of MiQ4 elements because the process of meshing always takes much time and computation cost in finite element analysis. Besides, the proposed element MiQ4 is suitable for the analysis of plates with arbitrary and complex geometries.  Table 11. It is obvious that the numerical results of the plates using the proposed element agree very well with those of CPT [26], the results of Nguyen et al. [18] using a refined simple FSDT, the solutions of Zenkour [26] using GSDT, and the 3D exact solutions [26]. In comparison with 3D solutions, the presented numerical results are minor difference.  [26] is shown in Table 12.

Bending of
e comparison shows that the numerical results of the plate using MiQ4 are identical to those of Zenkour [26] for both uniform and sinusoidal loads. e maximum error of the deflection is 1.35%, while the maximum error of the stresses is 1.83%. erefore, it can conclude that the proposed element is compatible to analyze the FGM structure. e nondimensional quantities are obtained as the following formulae: (70)

Parameter Study.
e proposed plate element MiQ4 is now applied to examine the static bending behavior of variable thickness FGM plates subjected to uniform load. e ingredients of the FGM plates are Al as metal and Al 2 O 3 as ceramic. Young's modulus is E m � 70 GPa for Al and E c � 380 GPa for Al 2 O 3 , while Poisson's ratios ] are constant and equal to 0.3, and the side-to-thickness ratio is a/h 0 � 10. e material properties of the plate vary through the thickness direction by the power-law function. e nondimensional quantities are computed as the following formulae: xy (a, 0, z).

(71)
It is noticed that the proposed element, MiQ4, is of good accuracy and convergence with coarse meshes, but the plates are divided with a fine mesh because of its thickness variation. Accordingly, a regular fine mesh of N × N � 32 × 32 is handled in the rest of this study.

e Influence of the Power-Law
Index. Firstly, we study the effects of the power-law index p on the deflection of the variable thickness FGM plates. A square FGM plate with the Type D (d) Figure 6: Four types of distorted mesh of square plates.          side-to-thickness ratio of a/h 0 � 10 and the varied parameter ζ � 0.5 is considered. e maximum nondimensional deflections of the plate for four cases of the BCs including fully clamped plates (CCCC), two opposite sides are simply supported and two other sides are clamped plates (SCSC), fully simple supported plates (SSSS), and two opposite sides are clamped and two other sides are free plates (CFCF) with different values of the power-law index are presented in Table 13. According to Table 13, it can be seen clearly that the increasing power-law index leads to increasing the maximum deflection of the plates for all cases of BCs and four types of thickness variation. It is due to the fact that when the powerlaw index increases, the plate becomes a metal-rich plate. In a comparison of the deflection between four cases of thickness variation of the plates, it can be concluded that the maximum deflection of the S plate is greatest, while the maximum deflection of the U plate is smallest.

5.4.2.
e Eeffects of the Side-to-ickness Ratio. Secondly, the effects of the side-to-thickness ratio a/h 0 on the deflection of the variable thickness FGM plates are investigated in this subsection. In this investigation, a square FGM plate with p � 2 and ζ � 0.5 is considered. e maximum nondimensional deflections of the plate for four cases of the BCs and thickness variations with different values of the side-to-thickness ratio are shown in Table 14. It seems that the change of the side-to-thickness ratio has strong effects on the deflection of the plate. When the sideto-thickness ratio increases, the maximum nondimensional deflection of the plate decreases. In addition, the BCs affect strongly the deflection of the plates. e maximum deflection of the CCCC plate is smallest, while the maximum deflection of the CFCF plate is greatest.

5.4.3.
e Influence of the ickness Variation and BCs. Continuously, the effects of the thickness variation and BCs on the deflection and stresses of the variable thickness FGM plate are studied in this subsection. A square FGM plate with the power-law index of p � 1, the varying parameter of ζ � 0.5, and the side-to-thickness ratio of a/h 0 � 10 is examined. e dependents of the deflected shapes on the BCs and the type of thickness variation are exhibited in Figures [8][9][10][11] . e change of the BCs leads to the change of the deflected shape of the plates. Besides, the strong effects of the thickness variation on the deflection of the plates are that the maximum deflection of the varied plate will shift to the thin side of the plate except for U plates. e deflected shapes are symmetrical in the case of U plates, while they are asymmetrical in other cases. In the case of CFCF plates, the maximum deflection appears at the midpoint of the thin edge of the plates. Additionally, we investigate the dependence of the stresses of the variable thickness FGM plate on the type of thickness variation. e nondimensional normal stress σ * x and in-plane shear stress τ * xy of the plate for four cases of thickness variation are demonstrated in Figure 12 for the case of square SSSS FGM plate and p � 1, ζ � 0.5. e stresses are strongly affected by the variation of the thickness. e maximum tensile normal stress and in-plane shear stress occur at the ceramic-rich surface of the plates. Furthermore, the maximum stresses of the S plates are greatest in comparison to the other cases, so it can lead to the stress concentration in the variable thickness structures. e stress distribution also depends on the law of thickness variation. e stresses at the middle FGM plate does not equal zero as in isotropic homogeneous; it moves toward the ceramic-rich surface of the FGM plates; especially, in the case of P plates, the normal stress equals zero at two points of the vertical fiber; it differs with other cases.
Finally, the effects of the varied parameter ζ on the maximum deflection of the plates are studied. A square FGM plate with p � 1 and a/h 0 � 10 is considered. e varied parameter varies in the range of 0÷0.5. It can see from Figure 13 that when the varied parameter increases, the maximum deflection of the plates increases rapidly for all cases of thickness variation except for U plate. e speed of the increase depends on the type of the thickness variation in the order of P plate, L plate, and S plate. It also seems that the  deflections of the CFCF plates are greatest while the CCCC plates are smallest. Additionally, the increasing rate of deflection of the CFCF plate is fastest.

Conclusion
In this work, a new mixed four-node quadrilateral plate element named MiQ4 with only five DOFs per node has been developed based on the FSDT and the mixed finite element formulation. Several examples were given to demonstrate the accuracy, high convergence rate, and efficiency of the proposed element. en the proposed element MiQ4 is applied to research the static bending behavior of the FGM plates with variable thickness. Based on the results of the presented work, some remarkable conclusions can be clarified as follows: e proposed element MiQ4 utilizes the bilinear shape function which is simple and similar to that of the wellknown classical Q4 element and allows us to calculate easily the numerical integration.
e new plate element is of good accuracy, high convergence rate, and reliability, especially for thin and very thin plates, in competition with other available elements.
e new element is insensitive to the mesh distortion, so it can be used to examine arbitrary and complex geometrical plates; as a consequence, the computation cost of the meshing process is reduced. e proposed element is free of shear locking without employing selective or reduced integrations. Some numerical examples with different power-law index, side-to-thickness ratio, BCs, and thickness variation on the bending behavior of FGM plates are investigated.
e results of this study suggest some future works such as investigating the complex geometrical plates, studying the free vibration and dynamic response of the FGM plates subjected to complex loading in the thermal environment, and improving the convergent rate of the presented element.

Data Availability
No data were used to support this study.

Conflicts of Interest
e author declares that there are no conflicts of interest.