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

Authors: Pham Van Vinh
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

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 first-order shear deformation and mixed finite 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.


---

## Body

## 1. 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–4]. FGMs are usually produced by mixing two or more different ingredients with the volume varying continuously in the spatial direction. The 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. That is the reason why numerous researchers focus on investigating the dynamic as well as static behaviors of the FGM structures, for example, beams [5–10], plates [11–18], shells [19–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), higher-order shear deformation theory (HSDT), and quasi-3D shear deformation theory. Thai 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 Thai 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. The 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. These 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. Thai 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. Therefore, 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. These 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. The mixed shear projected approach element (MiSP) was established by Ayad and his colleagues [44, 45]. The combination of the discrete shear gap technique and alternative alpha FEM was used by Nguyen-Thanh 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–51], edge-based SFEM [52], and node-based SFEM [53]. The 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. The generalized conforming elements were first developed by Long et al. [54–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. The MITC4 is an efficient element with thick and moderate plates, but less accuracy with very thin plates. The 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. This paper aims to establish a new mixed four-node quadrilateral plate element based on the mixed finite element formulation and the FSDT. This is the main contribution of the current work. The 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. The 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. Then 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.


## 2. Modeling of Variable Thickness FGM Plate

In this study, a variable thickness FGM plate with the dimension ofa×b and thickness of h as shown in Figure 1 is considered. The properties of the material vary continuously in the thickness direction of the plates by a power-law function as the following formula [1, 14]:(1)Pz=Pm+Pc−Pm2z+h2hp,where Pc and Pm 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=hx, 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:(2)h=hx=h0,UniformUplate,h01−ζxa,LinearLplate,h01−ζxa2,ParabolicPplate,h01−ζsinπx2a,SinusoidalSplate,where ζ is the varied parameter and h0 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.Figure 1
Modeling of variable thickness FGM plates.


## 3. Governing Equations

The displacement field of the plate based on FSDT is given as follows:(3)ux,y,z=ux,y+zψxx,y,vx,y,z=vx,y+zψyx,y,wx,y,z=wx,y,where u,v,w, ψx,ψy are five unknown displacement functions at the middle surface of the plate.The strain field of the plates is obtained as the following formulae:(4)εx=u,x+zψx,x,εy=v,y+zψy,y,γxy=v,x+u,y+zψx,y+ψy,x,γxz=w,x+ψx,γyz=w,y+ψy.The short form of equation (4) is(5)ε=ε0+zε1,where(6)ε0=u,xv,yu,y+v,x,ε1=ψx,xψy,yψx,y+ψy,x.The vector of the transverse shear strain can be rewritten as(7)γ=w,x+ψxw,y+ψy.The stress field of the plate is obtained as follows:(8)σ=Dε,τ=Dsγ.Here,(9)D=E1−ν21ν0ν10001−ν2,Ds=kE21+ν1001,where E is Young’s modulus of the material and k is the shear correction factor.The strain energy of the plate and the work done by external force can be calculated as the following formulae:(10)Πst=12∫VεT.σ+γT.τdV,Πext=∫SqwdS.Hamilton’s principle is used to obtain the governing equation of motion of the plate(11)δΠ=δΠst+δΠext=0.The variation of the strain energy is(12)δΠst=12δ∫VεTσ+δγTτdV=∫Vδε0+zδε1TDδε0+zδε1+δγTτdV.Equation (12) can be rewritten in the following compact form:(13)δΠst=∫Vδε0δε1DzDzDz2Dε0ε1+δγTτdV.By integrating equation (13) across the thickness, one gets(14)δΠst=∫SδΩTP+δγTQdS,where P and Q are given by(15)P=∫−hx/2hx/2DzDzDz2Dε0ε1dz,(16)Q=∫−hx/2hx/2τdz=∫−hx/2hx/2Dsγdz,and(17)δΩ=δε0δε1.Substituting equation (17) into equations (15) and (16), one gets(18)P=HΩ,Q=Hsγ,where(19)H=∫−hx/2hx/2DzDzDz2Ddz,(20)Hs=∫−hx/2hx/2Dsdz.


## 4. Finite Element Formulation

### 4.1. 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. The nodal displacement vector of thei‐th node is(21)qi=uiviwiψxiψyiT,i=1,4¯.The nodal displacement vector of the plate element,U, is defined as(22)U=q1Tq2Tq3Tq4TT.Then the displacement variables at any point of the element are approximated as the following formulae:(23)u=∑Niui,v=∑Nivi,w=∑Niwi,ψx=∑Niψxi,ψy=∑Niψyi,where Ni is the bilinear shape functions of the four-node quadrilateral plate element. The expressions of the shape functions in the local coordinate are(24)Ni=141+ξiξ1+ηiη,i=1,4¯.Equation (17) can be rewritten in short form as the following formula:(25)δΩ=BδU.Here,(26)B=B1B2B3B4,with(27)Bi=Ni,x00000Ni,y000Ni,yNi,x000000Ni,x00000Ni,y000Ni,yNi,x,i=1,4¯.The variation of the shear strain vector is obtained as follows:(28)δγ=δw,x+δψxδw,y+δψy.Equation (28) can be written in short form as(29)δγ=BsδU,where(30)Bs=Bs1Bs2Bs3Bs4,with(31)Bsi=00Ni,xNi000Ni,y0Ni,i=1,4¯.By substituting equations (25) and (29) into the expression of the variation of the strain energy, we get(32)δΠst=∫SδUTBTHBUdS+∫SδUTBsTHsBsUdS.The variation of the work done by an external force is calculated by(33)δΠext=−∫SqδwdS=−δUT∫SNwTqdS.Here,Nw is the matrix of the shape functions related to the transverse displacement; its expression is as the following formula:(34)Nw=Nw1Nw2Nw3Nw4,where(35)Nwi=00Ni00,i=1,4¯.Inserting equations (32) and (33) into equation (11) and using the trivial manner of classical FEM, one gets(36)KU=f.Here, the stiffness matrix and the nodal force vector of the plate element are computed by the following formulae:(37)K=Kb+Ks,where(38)Kb=∫SBTHBdS=∫−11∫−11BTHBJdξdη,Ks=∫SBsTHsBsdS=∫−11∫−11BsTHsBsJdξdη,f=∫SNwTqdS=∫−11∫−11NwTqJdξdη.The Gauss integration is used to estimate the stiffness matrixKb,Ks 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 Kb and f are calculated by using full Gauss integration 2×2, while the shear stiffness matrix Ks is estimated by using reduced Gauss integration 1×1 to avoid the well-known 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.


### 4.2. 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:(39)ψx=∑Niψxi+Nmψxm,ψy=∑Niψyi+Nmψym,where(40)Nm=1161−ξ21−η2.Substituting equation (39) into equation (25), the expression of δε1 can be obtained as follows:(41)δε1=∑Ni,xδψxi∑Ni,yδψyi∑Ni,yδψxi+∑Ni,xδψyi+Nm,xδψxmNm,yδψymNm,yδψxm+Nm,xδψym,or(42)δε1=∑Ni,xδψxi∑Ni,yδψyi∑Ni,yδψxi+∑Ni,xδψyi+Nm,x00Nm,yNm,yNm,xδψxmδψym.The expression of the strain vector is obtained as follows:(43)δΩ=BδU+Bmδψm,where(44)Bm=0Nm,x00Nm,yNm,yNm,x,δψm=δψxmδψym.The shear strainδγ is expressed as the following formula:(45)δγ=∑Ni,xδwi+∑Niδψxi+Nmδψxm∑Ni,yδwi+∑Niδψyi+Nmδψym,or(46)δγ=∑Ni,xδwi+∑Niδψxi∑Ni,yδwi+∑Niδψyi+Nm00Nmδψxmδψym.We rewrite equation (46) in the matrix form as follows:(47)δγ=BsδU+Mmδψm,where(48)Mm=Nm00Nm,δψm=δψxmδψym.We introduce a constant shear strain fieldγ¯=Hs−1Q0, and following the procedure of the mixed FEM, the expression of the potential energy is evaluated as [62](49)Π=12∫SΩTHΩ+γTHsγ¯+γ¯THsγ−γ¯THsγ¯dS+Πext+Πbc,where Πbc is the energy of the constrain to consider the BC effects.The variations of equation (49) are(50)∫SδΩTHΩ+δγTHsγ¯dS−δΠext=0,(51)∫Sδγ¯THsγ−γ¯dS=0.Inserting equations (43) and (47) into equations (50) and (51), one gets(52)∫SBδU+BmδψmTHBU+Bmψm+BsδU+MmδψmTHsγ¯dS+δΠext=0,(53)∫Sδγ¯THsBsU+Mmψm−γ¯dS=0.Insertingγ¯=Hs−1Q0 into equations (52) and (53), one gets(54)∫SBδU+BmδψmTHBU+Bmψm+BsδU+MmδψmTQ0dS+δΠext=0.(55)∫SδT0TBsU+Mmψm−Hs−1Q0dS=0,or(56)∫SδUTBTHBU+δUTBTHBmψm+δψmTBmTHBU+δψmTBmTHBmψm+δUTBsTQ0+δψmTMmTQ0dS−δUT∫SNwTqdS=0,(57)∫SδQ0TBsU+δQ0TMmψm+δQ0T−Hs−1Q0dS=0.By combining equations (53) and (57), one gets(58)δUTδQ0TδψmTKuuKuTKuψKuTTKTTKTψKuψTKTψTKψψUQ0ψm−fu00=0,where(59)Kuu=∫SBTHBdS=∫−11∫−11BTHBJdξdη,(60)Kuψ=∫SBTHBmdS=∫−11∫−11BTHBmJdξdη,(61)Kψψ=∫SBmTHBmdS=∫−11∫−11BmTHBmJdξdη.(62)KuT=∫SBsTdS=∫−11∫−11BsTJdξdη,(63)KTψ=∫SMmdS=∫−11∫−11MmJdξdη,(64)KTT=−∫SHs−1dS=−∫−11∫−11Hs−1Jdξdη,(65)fu=∫SNwTqdS=∫−11∫−11NwTqJdξdη,where J is the determinant of the Jacobian matrix.By eliminatingψm and Q0 from equation (58) at the element level, one gets(66)KU=f,in which(67)K=K11−K12K22−1K21,(68)f=fu,where(69)K11=Kuu−KuψKψψ−1KuψT,K12=KuT−KuψKψψ−1KTψT,K21=KuTT−KTψKψψ−1KuψT,K22=KTT−KTψKψψ−1KTψT.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. They are substantial benefits of the proposed plate element in comparison with the classical four-node quadrilateral element Q4 and other four-node quadrilateral elements.


## 4.1. 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. The nodal displacement vector of thei‐th node is(21)qi=uiviwiψxiψyiT,i=1,4¯.The nodal displacement vector of the plate element,U, is defined as(22)U=q1Tq2Tq3Tq4TT.Then the displacement variables at any point of the element are approximated as the following formulae:(23)u=∑Niui,v=∑Nivi,w=∑Niwi,ψx=∑Niψxi,ψy=∑Niψyi,where Ni is the bilinear shape functions of the four-node quadrilateral plate element. The expressions of the shape functions in the local coordinate are(24)Ni=141+ξiξ1+ηiη,i=1,4¯.Equation (17) can be rewritten in short form as the following formula:(25)δΩ=BδU.Here,(26)B=B1B2B3B4,with(27)Bi=Ni,x00000Ni,y000Ni,yNi,x000000Ni,x00000Ni,y000Ni,yNi,x,i=1,4¯.The variation of the shear strain vector is obtained as follows:(28)δγ=δw,x+δψxδw,y+δψy.Equation (28) can be written in short form as(29)δγ=BsδU,where(30)Bs=Bs1Bs2Bs3Bs4,with(31)Bsi=00Ni,xNi000Ni,y0Ni,i=1,4¯.By substituting equations (25) and (29) into the expression of the variation of the strain energy, we get(32)δΠst=∫SδUTBTHBUdS+∫SδUTBsTHsBsUdS.The variation of the work done by an external force is calculated by(33)δΠext=−∫SqδwdS=−δUT∫SNwTqdS.Here,Nw is the matrix of the shape functions related to the transverse displacement; its expression is as the following formula:(34)Nw=Nw1Nw2Nw3Nw4,where(35)Nwi=00Ni00,i=1,4¯.Inserting equations (32) and (33) into equation (11) and using the trivial manner of classical FEM, one gets(36)KU=f.Here, the stiffness matrix and the nodal force vector of the plate element are computed by the following formulae:(37)K=Kb+Ks,where(38)Kb=∫SBTHBdS=∫−11∫−11BTHBJdξdη,Ks=∫SBsTHsBsdS=∫−11∫−11BsTHsBsJdξdη,f=∫SNwTqdS=∫−11∫−11NwTqJdξdη.The Gauss integration is used to estimate the stiffness matrixKb,Ks 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 Kb and f are calculated by using full Gauss integration 2×2, while the shear stiffness matrix Ks is estimated by using reduced Gauss integration 1×1 to avoid the well-known 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.


## 4.2. 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:(39)ψx=∑Niψxi+Nmψxm,ψy=∑Niψyi+Nmψym,where(40)Nm=1161−ξ21−η2.Substituting equation (39) into equation (25), the expression of δε1 can be obtained as follows:(41)δε1=∑Ni,xδψxi∑Ni,yδψyi∑Ni,yδψxi+∑Ni,xδψyi+Nm,xδψxmNm,yδψymNm,yδψxm+Nm,xδψym,or(42)δε1=∑Ni,xδψxi∑Ni,yδψyi∑Ni,yδψxi+∑Ni,xδψyi+Nm,x00Nm,yNm,yNm,xδψxmδψym.The expression of the strain vector is obtained as follows:(43)δΩ=BδU+Bmδψm,where(44)Bm=0Nm,x00Nm,yNm,yNm,x,δψm=δψxmδψym.The shear strainδγ is expressed as the following formula:(45)δγ=∑Ni,xδwi+∑Niδψxi+Nmδψxm∑Ni,yδwi+∑Niδψyi+Nmδψym,or(46)δγ=∑Ni,xδwi+∑Niδψxi∑Ni,yδwi+∑Niδψyi+Nm00Nmδψxmδψym.We rewrite equation (46) in the matrix form as follows:(47)δγ=BsδU+Mmδψm,where(48)Mm=Nm00Nm,δψm=δψxmδψym.We introduce a constant shear strain fieldγ¯=Hs−1Q0, and following the procedure of the mixed FEM, the expression of the potential energy is evaluated as [62](49)Π=12∫SΩTHΩ+γTHsγ¯+γ¯THsγ−γ¯THsγ¯dS+Πext+Πbc,where Πbc is the energy of the constrain to consider the BC effects.The variations of equation (49) are(50)∫SδΩTHΩ+δγTHsγ¯dS−δΠext=0,(51)∫Sδγ¯THsγ−γ¯dS=0.Inserting equations (43) and (47) into equations (50) and (51), one gets(52)∫SBδU+BmδψmTHBU+Bmψm+BsδU+MmδψmTHsγ¯dS+δΠext=0,(53)∫Sδγ¯THsBsU+Mmψm−γ¯dS=0.Insertingγ¯=Hs−1Q0 into equations (52) and (53), one gets(54)∫SBδU+BmδψmTHBU+Bmψm+BsδU+MmδψmTQ0dS+δΠext=0.(55)∫SδT0TBsU+Mmψm−Hs−1Q0dS=0,or(56)∫SδUTBTHBU+δUTBTHBmψm+δψmTBmTHBU+δψmTBmTHBmψm+δUTBsTQ0+δψmTMmTQ0dS−δUT∫SNwTqdS=0,(57)∫SδQ0TBsU+δQ0TMmψm+δQ0T−Hs−1Q0dS=0.By combining equations (53) and (57), one gets(58)δUTδQ0TδψmTKuuKuTKuψKuTTKTTKTψKuψTKTψTKψψUQ0ψm−fu00=0,where(59)Kuu=∫SBTHBdS=∫−11∫−11BTHBJdξdη,(60)Kuψ=∫SBTHBmdS=∫−11∫−11BTHBmJdξdη,(61)Kψψ=∫SBmTHBmdS=∫−11∫−11BmTHBmJdξdη.(62)KuT=∫SBsTdS=∫−11∫−11BsTJdξdη,(63)KTψ=∫SMmdS=∫−11∫−11MmJdξdη,(64)KTT=−∫SHs−1dS=−∫−11∫−11Hs−1Jdξdη,(65)fu=∫SNwTqdS=∫−11∫−11NwTqJdξdη,where J is the determinant of the Jacobian matrix.By eliminatingψm and Q0 from equation (58) at the element level, one gets(66)KU=f,in which(67)K=K11−K12K22−1K21,(68)f=fu,where(69)K11=Kuu−KuψKψψ−1KuψT,K12=KuT−KuψKψψ−1KTψT,K21=KuTT−KTψKψψ−1KuψT,K22=KTT−KTψKψψ−1KTψT.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. They are substantial benefits of the proposed plate element in comparison with the classical four-node quadrilateral element Q4 and other four-node quadrilateral elements.


## 5. Numerical Results and Discussions

### 5.1. Convergence Study

In this subsection, the convergence studies are considered for square isotropic plates subjected to uniform load as shown in Figure2. The convergences of the proposed element MiQ4 with full Gauss integration are compared with those of classical Q4 with reduced Gauss integration. The nondimensional central deflections of reference wref.100D/qa4 in this subsection are calculated using the formulation given by Jirousek et al. [60]. The dimension and material properties of the plate are a=10,E=10.92, and ν=0.3 and the thickness-to-side ratio h/a of the plates is 0.1, 0.01, and 0.001, respectively. The plates are divided by a regular rectangular mesh of N×N.Figure 2
Model of square isotropic plates.
(a)(b)The comparisons of nondimensional central deflectionwc.100D/qa4 of the plates using MiQ4 and Q4 and reference results are presented in Tables 1–3, in which the error of the numerical results is computed in comparison to the solutions of Jirousek et al. [60]. The 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.Table 1
Nondimensional central deflectionwc.100D/qa4 of the square plate with h/a=0.1.
MeshSSSSCCCCQ4Error (%)MiQ4Error (%)Q4Error (%)MiQ4Error (%)2 × 20.354582−17.010.4633338.440.035714−76.260.22860951.944 × 40.420696−1.540.4345011.690.145856−3.060.17140813.926 × 60.424538−0.640.4300720.650.148648−1.200.1591825.808 × 80.425739−0.360.4287720.350.149469−0.660.1553453.2510 × 100.426292−0.230.4282120.220.149828−0.420.1535642.0612 × 120.426592−0.160.4279180.150.150021−0.290.1526071.4314 × 140.426773−0.120.4277440.110.150137−0.210.1520331.0516 × 160.426890−0.090.4276320.080.150212−0.160.1516620.80[60]0.4272840.150460Table 2
Central deflectionwc.100D/qa4 of the square plate with h/a=0.01.
MeshSSSSCCCCQ4Error (%)MiQ4Error (%)Q4Error (%)MiQ4Error (%)2 × 20.319225−21.460.4279765.300.000357−99.720.19325353.384 × 40.397125−2.290.4109661.110.121379−3.670.14604515.916 × 60.402650−0.930.4081920.430.124209−1.420.1351317.258 × 80.404346−0.520.4073810.230.125326−0.530.1314964.3610 × 100.405109−0.330.4070300.140.125850−0.120.1297853.0012 × 120.405518−0.230.4068440.100.1261340.110.1288652.2714 × 140.405763−0.170.4067340.070.1263050.240.1283101.8316 × 160.405921−0.130.4066630.050.1264170.330.1279511.55[60]0.4064460.12600Table 3
Nondimensional central deflectionwc.100D/qa4 of the square plate with h/a=0.001.
MeshSSSSCCCCQ4Error (%)MiQ4Error (%)Q4Error (%)MiQ4Error (%)2 × 20.318872−21.510.4276225.260.000004−100.000.19289953.094 × 40.396889−2.300.4107311.110.121123−3.870.14578315.706 × 60.402431−0.940.4079740.430.123950−1.630.1348837.058 × 80.404132−0.520.4071680.230.125068−0.740.1312464.1610 × 100.404897−0.330.4068180.140.125593−0.320.1295352.8112 × 120.405307−0.230.4066340.100.125879−0.100.1286152.0814 × 140.405553−0.170.4065240.070.1260510.040.1280601.6316 × 160.405711−0.130.4064530.050.1261630.130.1277011.35[60]0.4062370.12600Figure 3
Comparison of central deflectionwc.100D/qa4 of the square plate: (a) h/a=0.1, (b) h/a=0.01, and (c) h/a=0.001.
(a)(b)(c)Additionally, a benchmark test about the shear locking free of the proposed plate element MiQ4 is presented. The 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). The 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 correctly when the thickness-to-side ratio is smaller than h/a=10−8.Table 4
Nondimensional central deflectionwc.100D/qa4 of the square plate with different values of h/a.
BCsh/a10–310–410–610–810–1010–30SSSSQ40.4057110.4057090.4057010.3321420.1379.1020.2206.10−41MiQ40.4064530.4064510.4064510.4064510.4064510.406451Cen et al. [59]0.406240.406240.406240.406240.406240.40624Analytical [59]0.40620.40620.40620.40620.40620.4062CCCCQ40.1261630.1261610.1261590.1175940.4408470.9374.10−39MiQ40.1277010.1276980.1276980.1276980.1276980.127698Cen et al. [59]0.126530.126530.126530.126530.126530.12653Analytical [59]0.12650.12650.12650.12650.12650.1265


### 5.2. Distorted Study

#### 5.2.1. Morley’s Acute Skew Plate

Continuously, we consider a skew plate (see Figure4) subjected to a uniform load which is carried out by Morley [63]. The skew angle of the plate is β=30o and the BCs of the plate are soft simply supported (w=0) at all edges. The geometry and material properties of the plate are a=100,E=10.92, and ν=0.3 and the side-to-thickness ratio a/h is 100 and 1000.Figure 4
The model of Morley’s acute plate.The comparisons of the nondimensional central deflectionw∗=wc103D/qa4 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]. The 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. The 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.Figure 5
Comparison of the central deflection of Morley’s acute plate subjected to uniform load: (a)a/h=100 and (b) a/h=1000.
(a)(b)Table 5
Central deflection of Morley’s acute plate subjected to uniform load witha/h=100.
MeshDKMQ [61]Error (%)MITC4 [61]Error (%)MiQ4Error (%)2 × 20.757085.540.3590−12.010.3830−6.138 × 80.504023.530.3570−12.500.457312.0716 × 160.44108.090.3830−6.130.43556.7432 × 320.42303.680.4040−0.980.42865.04Morley [63]0.408Table 6
Central deflection of Morley’s acute plate subjected to uniform load witha/h=1000.
MeshDKMQ [61]Error (%)MITC4 [61]Error (%)MiQ4Error (%)2 × 20.760086.270.3580−12.250.3821−6.368 × 80.507024.260.3430−15.930.456211.8216 × 160.44308.580.3430−15.930.43396.3532 × 320.42504.170.3620−11.270.42584.37Morley [63]0.408


#### 5.2.2. 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. The dimension and material properties of the plate area=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. The 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. This 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.Figure 6
Four types of distorted mesh of square plates.
(a)(b)(c)(d)Table 7
Central deflectionwc.100D/qa4 of the square SSSS plates with h/a=0.1.
MeshType AType BType CType DQ4MiQ4Q4MiQ4Q4MiQ4Q4MiQ44 × 40.3915760.4391430.3789810.4212750.4003120.4262770.3555840.3936558 × 80.4195340.4294460.4178780.4258530.4144450.4236220.4032640.41437616 × 160.4252050.4276160.4250210.4269120.4230670.4256910.4179270.42214932 × 320.4267490.4273520.4267150.4271830.4260740.4267820.4273970.427632[60]0.427284Table 8
Central deflectionwc.100D/qa4 of the square SSSS plates with h/a=0.01.
MeshType AType BType CType DQ4MiQ4Q4MiQ4Q4MiQ4Q4MiQ44 × 40.2396790.4097800.2273670.3932470.1169950.4008560.0880910.3630068 × 80.3807370.4071290.3858810.4036340.1869030.3956330.1974000.38135616 × 160.4034290.4064970.4033690.4058340.2640440.3946330.2879500.38836132 × 320.4058120.4064500.4058040.4062990.3400330.3960930.3560150.395804[60]0.406446Table 9
Central deflectionwc.100D/qa4 of the square plates with h/a=0.001.
MeshType AType BType CType DQ4MiQ4Q4MiQ4Q4MiQ4Q4MiQ44 × 40.1146210.4094540.1421550.3929420.0345040.4005890.0269150.3626348 × 80.3232620.4068920.3634500.4033850.0940290.3952050.1011940.38071816 × 160.3913830.4062800.3989100.4056070.1416530.3934400.1775680.38677532 × 320.4049630.4062390.4051680.4060840.1727640.3924160.2161360.391490[60]0.406237Table 10
Central deflectionwc.100D/qa4 of the square SSSS plates with h/a=10−6.
MeshType AType BType CType DQ4MiQ4Q4MiQ4Q4MiQ4Q4MiQ44 × 40.1180060.4094510.1107970.3929390.0240060.4005870.0202520.3626308 × 80.3183350.4068890.3706210.4033820.0581590.3952000.1137830.38071216 × 160.0000000.4062780.7232660.4056050.1322560.3934270.1217220.38675832 × 320.3824670.4062370.0000000.4060820.3280340.3923660.0000000.391430[59]0.40624Figure 7
Comparison of central deflection of the square SSSS plates: (a)h/a=0.1, (b) h/a=0.01, (c) h/a=0.001, (d) h/a=10−6.
(a)(b)(c)(d)


### 5.3. Validation Study

#### 5.3.1. Bending of Isotropic Homogeneous Plates

In this subsection, the static bending of an isotropic homogeneous plate under uniform load is considered to prove the validity of the new plate element MiQ4. The geometrical and material properties of the quadratic plate area=b=1,E=1,q=1, and ν=0.3, while the thickness of the plate is considered by three values of h=0.01,h=0.03, and h=0.1. The comparison of the numerical results of the plates using the proposed element MiQ4 with other published results is presented in 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.Table 11
The central deflection of isotropic homogeneous plates under uniform load.
hClassical [26]RS-FSDT [18]SSDT [26]3D [26]PresentDifferent with 3D (%)0.0144360.944385.4144383.8444384.744388.77−0.0090.031643.001651.1691650.6461650.941650.85−0.0050.144.360946.8127146.6548146.744346.6678−0.164


#### 5.3.2. Bending of Isotropic FGM Plates

Continuously, an FGM square plate of Al/Al2O3 subjected to uniform load and sinusoidal load is considered to exhibit the accuracy of the proposed element MiQ4. Young’s modulus is Em=70GPa for Al and Ec=380GPa for Al2O3, while Poisson’s ratios ν are constant and equal 0.3, and the side-to-thickness ratio is a/h=10. The material properties of the plate vary through the thickness direction by the power-law function. The comparison of nondimensional displacement and stresses using MiQ4 and those of Zenkour [26] is shown in Table 12. The comparison shows that the numerical results of the plate using MiQ4 are identical to those of Zenkour [26] for both uniform and sinusoidal loads. The maximum error of the deflection is 1.35%, while the maximum error of the stresses is 1.83%. Therefore, it can conclude that the proposed element is compatible to analyze the FGM structure. The nondimensional quantities are obtained as the following formulae:(70)w¯=10Ech3q0a4wa2,b2,σ¯x=hq0aσxa2,b2,h2,σ¯y=hq0aσya2,b2,h3,τ¯xy=hq0aτxy0,0,−h3.Table 12
Nondimensional deflection and stresses of the square SSSS FGM plates.
pMethodUniform loadSinusoidal loadw¯σ¯xσ¯yτ¯xyw¯σ¯xσ¯yτ¯xyCeramicZenkour [26]0.46652.89321.91031.28500.29601.99551.31210.7065Present0.46672.86961.91311.28710.29601.96971.31310.7071Error (%)0.04−0.820.150.16−0.00−1.290.080.081Zenkour [26]0.92874.47452.16921.11430.58893.08701.48940.6110Present0.92894.43512.17381.11200.58873.04421.49210.6107Error (%)0.02−0.880.21−0.21−0.03−1.390.18−0.052Zenkour [26]1.19405.22962.03380.99070.75733.60941.39540.5441Present1.19105.17862.04140.99110.75483.55461.40120.5443Error (%)−0.25−0.970.370.04−0.32−1.520.410.043Zenkour [26]1.32005.61081.85931.00470.83773.87421.27480.5525Present1.31245.55041.86931.00710.83203.80981.28310.5531Error (%)−0.57−1.080.540.24−0.68−1.660.650.114Zenkour [26]1.38905.89151.71971.02980.88194.06931.17830.5667Present1.37715.82401.73141.03360.87323.99751.18840.5676Error (%)−0.86−1.150.680.37−0.99−1.760.860.175Zenkour [26]1.43566.15041.61041.04510.91184.24881.10290.5755Present1.42066.07781.62301.04980.90104.17171.11400.5766Error (%)−1.04−1.180.780.45−1.18−1.811.000.186Zenkour [26]1.47276.40431.52141.05360.93564.42441.04170.5803Present1.45576.32811.53431.05890.92344.34361.05310.5816Error (%)−1.16−1.190.850.51−1.30−1.831.100.227Zenkour [26]1.50496.65471.44671.05890.95624.59710.99030.5834Present1.48696.57611.45951.06460.94334.51381.00180.5847Error (%)−1.20−1.180.880.54−1.35−1.811.160.228Zenkour [26]1.53436.89991.38291.06280.97504.76610.94660.5856Present1.51596.81971.39531.06870.96194.68100.95770.5870Error (%)−1.20−1.160.900.56−1.35−1.791.180.239Zenkour [26]1.56177.13831.32831.06620.99254.93030.90920.5875Present1.54357.05701.34041.07230.97944.84380.92000.5889Error (%)−1.17−1.140.910.57−1.32−1.751.190.2410Zenkour [26]1.58767.36891.28201.06941.00895.08900.87750.5894Present1.56987.28671.29351.07570.99625.00150.88790.5908Error (%)−1.12−1.120.900.59−1.26−1.721.180.24MetalZenkour [26]2.53272.89321.91031.28501.60701.99551.31210.7065Present2.53342.86961.91311.28711.60661.96971.31310.7071Error (%)0.03−0.820.150.16−0.02−1.290.080.08


### 5.4. Parameter Study

The proposed plate element MiQ4 is now applied to examine the static bending behavior of variable thickness FGM plates subjected to uniform load. The ingredients of the FGM plates are Al as metal and Al2O3 as ceramic. Young’s modulus is Em=70GPa for Al and Ec=380GPa for Al2O3, while Poisson’s ratios ν are constant and equal to 0.3, and the side-to-thickness ratio is a/h0=10. The material properties of the plate vary through the thickness direction by the power-law function. The nondimensional quantities are computed as the following formulae:(71)wmax∗=10Ech03q0a4wmax,σx∗=h0q0aσxa2,b2,z,τxy∗=h0q0aτxya,0,z.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 ofN×N=32×32 is handled in the rest of this study.


#### 5.4.1. The Influence of the Power-Law Index

Firstly, we study the effects of the power-law indexp on the deflection of the variable thickness FGM plates. A square FGM plate with the side-to-thickness ratio of a/h0=10 and the varied parameter ζ=0.5 is considered. The 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 power-law 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.Table 13
The maximum nondimensional deflection of the FGM plates as a function ofp.
BCsThickness variationspCeramic0.51510MetalCCCCU plate0.16460.24960.32210.50230.56270.8935L plate0.40100.61240.79271.22111.35672.1762P plate0.28940.44090.57000.88150.98221.5707S plate0.55010.84161.09011.67421.85632.9854SCSCU plate0.24140.36720.47440.73630.82221.3104L plate0.59120.90441.17161.79971.99553.2085P plate0.44320.67690.87611.34901.49872.4054S plate0.77921.19331.54662.37072.62524.2286SSSSU plate0.46670.71600.92891.42061.56982.5334L plate1.16361.78952.32403.53983.90076.3151P plate0.82031.26041.63612.49532.75264.4518S plate1.62562.50123.24874.94305.44388.8222CFCFU plate0.36310.55290.71471.10681.23401.9708L plate1.64082.51813.26674.99115.51318.9048P plate1.20601.84942.39833.66924.05706.5452S plate2.21833.40624.41956.74557.446612.0390


#### 5.4.2. The Eeffects of the Side-to-Thickness Ratio

Secondly, the effects of the side-to-thickness ratioa/h0 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. The 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 side-to-thickness ratio increases, the maximum nondimensional deflection of the plate decreases. In addition, the BCs affect strongly the deflection of the plates. The maximum deflection of the CCCC plate is smallest, while the maximum deflection of the CFCF plate is greatest.Table 14
The maximum nondimensional deflection of the FGM plates as a function ofa/h0.
BCsThickness variationsa/h051020304050CCCCU plate0.57450.41340.37070.36250.35960.3582L plate1.24451.01650.95650.94510.94100.9392P plate0.93390.73110.67750.66740.66380.6621S plate1.65321.39761.33041.31761.31311.3110SCSCU plate0.80540.60870.55640.54650.54300.5414L plate1.78061.50231.42961.41591.41101.4088P plate1.36911.12331.05891.04671.04241.0404S plate2.29571.98261.90091.88551.88011.8776SSSSU plate1.34241.19101.15311.14611.14371.1425L plate3.19092.97902.92592.91602.91262.9110P plate2.28642.09722.04982.04102.03792.0365S plate4.40334.16354.10314.09194.08804.0862CFCFU plate1.18830.91690.84510.83070.82530.8227L plate4.65064.18744.05954.03314.02324.0184P plate3.47503.07442.96662.94492.93702.9331S plate6.21055.66425.50845.47515.46245.4562


#### 5.4.3. The Influence of the Thickness 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 ofp=1, the varying parameter of ζ=0.5, and the side-to-thickness ratio of a/h0=10 is examined. The dependents of the deflected shapes on the BCs and the type of thickness variation are exhibited in Figures 8–11 . The 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. The 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.Figure 8
Deflected shapes of the CCCC FGM plates depend on the thickness variations.Figure 9
Deflected shapes of the SCSC FGM plates depend on the thickness variations.Figure 10
Deflected shapes of the SSSS FGM plates depend on the thickness variations.Figure 11
Deflected shapes of the CFCF FGM plates depend on the thickness variations.Additionally, we investigate the dependence of the stresses of the variable thickness FGM plate on the type of thickness variation. The 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. The stresses are strongly affected by the variation of the thickness. The 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. The stress distribution also depends on the law of thickness variation. The 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.Figure 12
The distribution of normal and in-plane shear stresses through the thickness of the FGM plates.
(a)(b)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/h0=10 is considered. The 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. The 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.Figure 13
The influence of the thickness variation on the nondimensional maximum deflections of the FGM plates with different BCs: (a) CCCC plate, (b) SCSC plate, (c) SSSS plate, and (d) CFCF plate.
(a)(b)(c)(d)


## 5.1. Convergence Study

In this subsection, the convergence studies are considered for square isotropic plates subjected to uniform load as shown in Figure2. The convergences of the proposed element MiQ4 with full Gauss integration are compared with those of classical Q4 with reduced Gauss integration. The nondimensional central deflections of reference wref.100D/qa4 in this subsection are calculated using the formulation given by Jirousek et al. [60]. The dimension and material properties of the plate are a=10,E=10.92, and ν=0.3 and the thickness-to-side ratio h/a of the plates is 0.1, 0.01, and 0.001, respectively. The plates are divided by a regular rectangular mesh of N×N.Figure 2
Model of square isotropic plates.
(a)(b)The comparisons of nondimensional central deflectionwc.100D/qa4 of the plates using MiQ4 and Q4 and reference results are presented in Tables 1–3, in which the error of the numerical results is computed in comparison to the solutions of Jirousek et al. [60]. The 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.Table 1
Nondimensional central deflectionwc.100D/qa4 of the square plate with h/a=0.1.
MeshSSSSCCCCQ4Error (%)MiQ4Error (%)Q4Error (%)MiQ4Error (%)2 × 20.354582−17.010.4633338.440.035714−76.260.22860951.944 × 40.420696−1.540.4345011.690.145856−3.060.17140813.926 × 60.424538−0.640.4300720.650.148648−1.200.1591825.808 × 80.425739−0.360.4287720.350.149469−0.660.1553453.2510 × 100.426292−0.230.4282120.220.149828−0.420.1535642.0612 × 120.426592−0.160.4279180.150.150021−0.290.1526071.4314 × 140.426773−0.120.4277440.110.150137−0.210.1520331.0516 × 160.426890−0.090.4276320.080.150212−0.160.1516620.80[60]0.4272840.150460Table 2
Central deflectionwc.100D/qa4 of the square plate with h/a=0.01.
MeshSSSSCCCCQ4Error (%)MiQ4Error (%)Q4Error (%)MiQ4Error (%)2 × 20.319225−21.460.4279765.300.000357−99.720.19325353.384 × 40.397125−2.290.4109661.110.121379−3.670.14604515.916 × 60.402650−0.930.4081920.430.124209−1.420.1351317.258 × 80.404346−0.520.4073810.230.125326−0.530.1314964.3610 × 100.405109−0.330.4070300.140.125850−0.120.1297853.0012 × 120.405518−0.230.4068440.100.1261340.110.1288652.2714 × 140.405763−0.170.4067340.070.1263050.240.1283101.8316 × 160.405921−0.130.4066630.050.1264170.330.1279511.55[60]0.4064460.12600Table 3
Nondimensional central deflectionwc.100D/qa4 of the square plate with h/a=0.001.
MeshSSSSCCCCQ4Error (%)MiQ4Error (%)Q4Error (%)MiQ4Error (%)2 × 20.318872−21.510.4276225.260.000004−100.000.19289953.094 × 40.396889−2.300.4107311.110.121123−3.870.14578315.706 × 60.402431−0.940.4079740.430.123950−1.630.1348837.058 × 80.404132−0.520.4071680.230.125068−0.740.1312464.1610 × 100.404897−0.330.4068180.140.125593−0.320.1295352.8112 × 120.405307−0.230.4066340.100.125879−0.100.1286152.0814 × 140.405553−0.170.4065240.070.1260510.040.1280601.6316 × 160.405711−0.130.4064530.050.1261630.130.1277011.35[60]0.4062370.12600Figure 3
Comparison of central deflectionwc.100D/qa4 of the square plate: (a) h/a=0.1, (b) h/a=0.01, and (c) h/a=0.001.
(a)(b)(c)Additionally, a benchmark test about the shear locking free of the proposed plate element MiQ4 is presented. The 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). The 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 correctly when the thickness-to-side ratio is smaller than h/a=10−8.Table 4
Nondimensional central deflectionwc.100D/qa4 of the square plate with different values of h/a.
BCsh/a10–310–410–610–810–1010–30SSSSQ40.4057110.4057090.4057010.3321420.1379.1020.2206.10−41MiQ40.4064530.4064510.4064510.4064510.4064510.406451Cen et al. [59]0.406240.406240.406240.406240.406240.40624Analytical [59]0.40620.40620.40620.40620.40620.4062CCCCQ40.1261630.1261610.1261590.1175940.4408470.9374.10−39MiQ40.1277010.1276980.1276980.1276980.1276980.127698Cen et al. [59]0.126530.126530.126530.126530.126530.12653Analytical [59]0.12650.12650.12650.12650.12650.1265


## 5.2. Distorted Study

### 5.2.1. Morley’s Acute Skew Plate

Continuously, we consider a skew plate (see Figure4) subjected to a uniform load which is carried out by Morley [63]. The skew angle of the plate is β=30o and the BCs of the plate are soft simply supported (w=0) at all edges. The geometry and material properties of the plate are a=100,E=10.92, and ν=0.3 and the side-to-thickness ratio a/h is 100 and 1000.Figure 4
The model of Morley’s acute plate.The comparisons of the nondimensional central deflectionw∗=wc103D/qa4 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]. The 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. The 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.Figure 5
Comparison of the central deflection of Morley’s acute plate subjected to uniform load: (a)a/h=100 and (b) a/h=1000.
(a)(b)Table 5
Central deflection of Morley’s acute plate subjected to uniform load witha/h=100.
MeshDKMQ [61]Error (%)MITC4 [61]Error (%)MiQ4Error (%)2 × 20.757085.540.3590−12.010.3830−6.138 × 80.504023.530.3570−12.500.457312.0716 × 160.44108.090.3830−6.130.43556.7432 × 320.42303.680.4040−0.980.42865.04Morley [63]0.408Table 6
Central deflection of Morley’s acute plate subjected to uniform load witha/h=1000.
MeshDKMQ [61]Error (%)MITC4 [61]Error (%)MiQ4Error (%)2 × 20.760086.270.3580−12.250.3821−6.368 × 80.507024.260.3430−15.930.456211.8216 × 160.44308.580.3430−15.930.43396.3532 × 320.42504.170.3620−11.270.42584.37Morley [63]0.408


### 5.2.2. 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. The dimension and material properties of the plate area=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. The 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. This 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.Figure 6
Four types of distorted mesh of square plates.
(a)(b)(c)(d)Table 7
Central deflectionwc.100D/qa4 of the square SSSS plates with h/a=0.1.
MeshType AType BType CType DQ4MiQ4Q4MiQ4Q4MiQ4Q4MiQ44 × 40.3915760.4391430.3789810.4212750.4003120.4262770.3555840.3936558 × 80.4195340.4294460.4178780.4258530.4144450.4236220.4032640.41437616 × 160.4252050.4276160.4250210.4269120.4230670.4256910.4179270.42214932 × 320.4267490.4273520.4267150.4271830.4260740.4267820.4273970.427632[60]0.427284Table 8
Central deflectionwc.100D/qa4 of the square SSSS plates with h/a=0.01.
MeshType AType BType CType DQ4MiQ4Q4MiQ4Q4MiQ4Q4MiQ44 × 40.2396790.4097800.2273670.3932470.1169950.4008560.0880910.3630068 × 80.3807370.4071290.3858810.4036340.1869030.3956330.1974000.38135616 × 160.4034290.4064970.4033690.4058340.2640440.3946330.2879500.38836132 × 320.4058120.4064500.4058040.4062990.3400330.3960930.3560150.395804[60]0.406446Table 9
Central deflectionwc.100D/qa4 of the square plates with h/a=0.001.
MeshType AType BType CType DQ4MiQ4Q4MiQ4Q4MiQ4Q4MiQ44 × 40.1146210.4094540.1421550.3929420.0345040.4005890.0269150.3626348 × 80.3232620.4068920.3634500.4033850.0940290.3952050.1011940.38071816 × 160.3913830.4062800.3989100.4056070.1416530.3934400.1775680.38677532 × 320.4049630.4062390.4051680.4060840.1727640.3924160.2161360.391490[60]0.406237Table 10
Central deflectionwc.100D/qa4 of the square SSSS plates with h/a=10−6.
MeshType AType BType CType DQ4MiQ4Q4MiQ4Q4MiQ4Q4MiQ44 × 40.1180060.4094510.1107970.3929390.0240060.4005870.0202520.3626308 × 80.3183350.4068890.3706210.4033820.0581590.3952000.1137830.38071216 × 160.0000000.4062780.7232660.4056050.1322560.3934270.1217220.38675832 × 320.3824670.4062370.0000000.4060820.3280340.3923660.0000000.391430[59]0.40624Figure 7
Comparison of central deflection of the square SSSS plates: (a)h/a=0.1, (b) h/a=0.01, (c) h/a=0.001, (d) h/a=10−6.
(a)(b)(c)(d)


## 5.2.1. Morley’s Acute Skew Plate

Continuously, we consider a skew plate (see Figure4) subjected to a uniform load which is carried out by Morley [63]. The skew angle of the plate is β=30o and the BCs of the plate are soft simply supported (w=0) at all edges. The geometry and material properties of the plate are a=100,E=10.92, and ν=0.3 and the side-to-thickness ratio a/h is 100 and 1000.Figure 4
The model of Morley’s acute plate.The comparisons of the nondimensional central deflectionw∗=wc103D/qa4 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]. The 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. The 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.Figure 5
Comparison of the central deflection of Morley’s acute plate subjected to uniform load: (a)a/h=100 and (b) a/h=1000.
(a)(b)Table 5
Central deflection of Morley’s acute plate subjected to uniform load witha/h=100.
MeshDKMQ [61]Error (%)MITC4 [61]Error (%)MiQ4Error (%)2 × 20.757085.540.3590−12.010.3830−6.138 × 80.504023.530.3570−12.500.457312.0716 × 160.44108.090.3830−6.130.43556.7432 × 320.42303.680.4040−0.980.42865.04Morley [63]0.408Table 6
Central deflection of Morley’s acute plate subjected to uniform load witha/h=1000.
MeshDKMQ [61]Error (%)MITC4 [61]Error (%)MiQ4Error (%)2 × 20.760086.270.3580−12.250.3821−6.368 × 80.507024.260.3430−15.930.456211.8216 × 160.44308.580.3430−15.930.43396.3532 × 320.42504.170.3620−11.270.42584.37Morley [63]0.408


## 5.2.2. 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. The dimension and material properties of the plate area=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. The 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. This 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.Figure 6
Four types of distorted mesh of square plates.
(a)(b)(c)(d)Table 7
Central deflectionwc.100D/qa4 of the square SSSS plates with h/a=0.1.
MeshType AType BType CType DQ4MiQ4Q4MiQ4Q4MiQ4Q4MiQ44 × 40.3915760.4391430.3789810.4212750.4003120.4262770.3555840.3936558 × 80.4195340.4294460.4178780.4258530.4144450.4236220.4032640.41437616 × 160.4252050.4276160.4250210.4269120.4230670.4256910.4179270.42214932 × 320.4267490.4273520.4267150.4271830.4260740.4267820.4273970.427632[60]0.427284Table 8
Central deflectionwc.100D/qa4 of the square SSSS plates with h/a=0.01.
MeshType AType BType CType DQ4MiQ4Q4MiQ4Q4MiQ4Q4MiQ44 × 40.2396790.4097800.2273670.3932470.1169950.4008560.0880910.3630068 × 80.3807370.4071290.3858810.4036340.1869030.3956330.1974000.38135616 × 160.4034290.4064970.4033690.4058340.2640440.3946330.2879500.38836132 × 320.4058120.4064500.4058040.4062990.3400330.3960930.3560150.395804[60]0.406446Table 9
Central deflectionwc.100D/qa4 of the square plates with h/a=0.001.
MeshType AType BType CType DQ4MiQ4Q4MiQ4Q4MiQ4Q4MiQ44 × 40.1146210.4094540.1421550.3929420.0345040.4005890.0269150.3626348 × 80.3232620.4068920.3634500.4033850.0940290.3952050.1011940.38071816 × 160.3913830.4062800.3989100.4056070.1416530.3934400.1775680.38677532 × 320.4049630.4062390.4051680.4060840.1727640.3924160.2161360.391490[60]0.406237Table 10
Central deflectionwc.100D/qa4 of the square SSSS plates with h/a=10−6.
MeshType AType BType CType DQ4MiQ4Q4MiQ4Q4MiQ4Q4MiQ44 × 40.1180060.4094510.1107970.3929390.0240060.4005870.0202520.3626308 × 80.3183350.4068890.3706210.4033820.0581590.3952000.1137830.38071216 × 160.0000000.4062780.7232660.4056050.1322560.3934270.1217220.38675832 × 320.3824670.4062370.0000000.4060820.3280340.3923660.0000000.391430[59]0.40624Figure 7
Comparison of central deflection of the square SSSS plates: (a)h/a=0.1, (b) h/a=0.01, (c) h/a=0.001, (d) h/a=10−6.
(a)(b)(c)(d)


## 5.3. Validation Study

### 5.3.1. Bending of Isotropic Homogeneous Plates

In this subsection, the static bending of an isotropic homogeneous plate under uniform load is considered to prove the validity of the new plate element MiQ4. The geometrical and material properties of the quadratic plate area=b=1,E=1,q=1, and ν=0.3, while the thickness of the plate is considered by three values of h=0.01,h=0.03, and h=0.1. The comparison of the numerical results of the plates using the proposed element MiQ4 with other published results is presented in 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.Table 11
The central deflection of isotropic homogeneous plates under uniform load.
hClassical [26]RS-FSDT [18]SSDT [26]3D [26]PresentDifferent with 3D (%)0.0144360.944385.4144383.8444384.744388.77−0.0090.031643.001651.1691650.6461650.941650.85−0.0050.144.360946.8127146.6548146.744346.6678−0.164


### 5.3.2. Bending of Isotropic FGM Plates

Continuously, an FGM square plate of Al/Al2O3 subjected to uniform load and sinusoidal load is considered to exhibit the accuracy of the proposed element MiQ4. Young’s modulus is Em=70GPa for Al and Ec=380GPa for Al2O3, while Poisson’s ratios ν are constant and equal 0.3, and the side-to-thickness ratio is a/h=10. The material properties of the plate vary through the thickness direction by the power-law function. The comparison of nondimensional displacement and stresses using MiQ4 and those of Zenkour [26] is shown in Table 12. The comparison shows that the numerical results of the plate using MiQ4 are identical to those of Zenkour [26] for both uniform and sinusoidal loads. The maximum error of the deflection is 1.35%, while the maximum error of the stresses is 1.83%. Therefore, it can conclude that the proposed element is compatible to analyze the FGM structure. The nondimensional quantities are obtained as the following formulae:(70)w¯=10Ech3q0a4wa2,b2,σ¯x=hq0aσxa2,b2,h2,σ¯y=hq0aσya2,b2,h3,τ¯xy=hq0aτxy0,0,−h3.Table 12
Nondimensional deflection and stresses of the square SSSS FGM plates.
pMethodUniform loadSinusoidal loadw¯σ¯xσ¯yτ¯xyw¯σ¯xσ¯yτ¯xyCeramicZenkour [26]0.46652.89321.91031.28500.29601.99551.31210.7065Present0.46672.86961.91311.28710.29601.96971.31310.7071Error (%)0.04−0.820.150.16−0.00−1.290.080.081Zenkour [26]0.92874.47452.16921.11430.58893.08701.48940.6110Present0.92894.43512.17381.11200.58873.04421.49210.6107Error (%)0.02−0.880.21−0.21−0.03−1.390.18−0.052Zenkour [26]1.19405.22962.03380.99070.75733.60941.39540.5441Present1.19105.17862.04140.99110.75483.55461.40120.5443Error (%)−0.25−0.970.370.04−0.32−1.520.410.043Zenkour [26]1.32005.61081.85931.00470.83773.87421.27480.5525Present1.31245.55041.86931.00710.83203.80981.28310.5531Error (%)−0.57−1.080.540.24−0.68−1.660.650.114Zenkour [26]1.38905.89151.71971.02980.88194.06931.17830.5667Present1.37715.82401.73141.03360.87323.99751.18840.5676Error (%)−0.86−1.150.680.37−0.99−1.760.860.175Zenkour [26]1.43566.15041.61041.04510.91184.24881.10290.5755Present1.42066.07781.62301.04980.90104.17171.11400.5766Error (%)−1.04−1.180.780.45−1.18−1.811.000.186Zenkour [26]1.47276.40431.52141.05360.93564.42441.04170.5803Present1.45576.32811.53431.05890.92344.34361.05310.5816Error (%)−1.16−1.190.850.51−1.30−1.831.100.227Zenkour [26]1.50496.65471.44671.05890.95624.59710.99030.5834Present1.48696.57611.45951.06460.94334.51381.00180.5847Error (%)−1.20−1.180.880.54−1.35−1.811.160.228Zenkour [26]1.53436.89991.38291.06280.97504.76610.94660.5856Present1.51596.81971.39531.06870.96194.68100.95770.5870Error (%)−1.20−1.160.900.56−1.35−1.791.180.239Zenkour [26]1.56177.13831.32831.06620.99254.93030.90920.5875Present1.54357.05701.34041.07230.97944.84380.92000.5889Error (%)−1.17−1.140.910.57−1.32−1.751.190.2410Zenkour [26]1.58767.36891.28201.06941.00895.08900.87750.5894Present1.56987.28671.29351.07570.99625.00150.88790.5908Error (%)−1.12−1.120.900.59−1.26−1.721.180.24MetalZenkour [26]2.53272.89321.91031.28501.60701.99551.31210.7065Present2.53342.86961.91311.28711.60661.96971.31310.7071Error (%)0.03−0.820.150.16−0.02−1.290.080.08


## 5.3.1. Bending of Isotropic Homogeneous Plates

In this subsection, the static bending of an isotropic homogeneous plate under uniform load is considered to prove the validity of the new plate element MiQ4. The geometrical and material properties of the quadratic plate area=b=1,E=1,q=1, and ν=0.3, while the thickness of the plate is considered by three values of h=0.01,h=0.03, and h=0.1. The comparison of the numerical results of the plates using the proposed element MiQ4 with other published results is presented in 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.Table 11
The central deflection of isotropic homogeneous plates under uniform load.
hClassical [26]RS-FSDT [18]SSDT [26]3D [26]PresentDifferent with 3D (%)0.0144360.944385.4144383.8444384.744388.77−0.0090.031643.001651.1691650.6461650.941650.85−0.0050.144.360946.8127146.6548146.744346.6678−0.164


## 5.3.2. Bending of Isotropic FGM Plates

Continuously, an FGM square plate of Al/Al2O3 subjected to uniform load and sinusoidal load is considered to exhibit the accuracy of the proposed element MiQ4. Young’s modulus is Em=70GPa for Al and Ec=380GPa for Al2O3, while Poisson’s ratios ν are constant and equal 0.3, and the side-to-thickness ratio is a/h=10. The material properties of the plate vary through the thickness direction by the power-law function. The comparison of nondimensional displacement and stresses using MiQ4 and those of Zenkour [26] is shown in Table 12. The comparison shows that the numerical results of the plate using MiQ4 are identical to those of Zenkour [26] for both uniform and sinusoidal loads. The maximum error of the deflection is 1.35%, while the maximum error of the stresses is 1.83%. Therefore, it can conclude that the proposed element is compatible to analyze the FGM structure. The nondimensional quantities are obtained as the following formulae:(70)w¯=10Ech3q0a4wa2,b2,σ¯x=hq0aσxa2,b2,h2,σ¯y=hq0aσya2,b2,h3,τ¯xy=hq0aτxy0,0,−h3.Table 12
Nondimensional deflection and stresses of the square SSSS FGM plates.
pMethodUniform loadSinusoidal loadw¯σ¯xσ¯yτ¯xyw¯σ¯xσ¯yτ¯xyCeramicZenkour [26]0.46652.89321.91031.28500.29601.99551.31210.7065Present0.46672.86961.91311.28710.29601.96971.31310.7071Error (%)0.04−0.820.150.16−0.00−1.290.080.081Zenkour [26]0.92874.47452.16921.11430.58893.08701.48940.6110Present0.92894.43512.17381.11200.58873.04421.49210.6107Error (%)0.02−0.880.21−0.21−0.03−1.390.18−0.052Zenkour [26]1.19405.22962.03380.99070.75733.60941.39540.5441Present1.19105.17862.04140.99110.75483.55461.40120.5443Error (%)−0.25−0.970.370.04−0.32−1.520.410.043Zenkour [26]1.32005.61081.85931.00470.83773.87421.27480.5525Present1.31245.55041.86931.00710.83203.80981.28310.5531Error (%)−0.57−1.080.540.24−0.68−1.660.650.114Zenkour [26]1.38905.89151.71971.02980.88194.06931.17830.5667Present1.37715.82401.73141.03360.87323.99751.18840.5676Error (%)−0.86−1.150.680.37−0.99−1.760.860.175Zenkour [26]1.43566.15041.61041.04510.91184.24881.10290.5755Present1.42066.07781.62301.04980.90104.17171.11400.5766Error (%)−1.04−1.180.780.45−1.18−1.811.000.186Zenkour [26]1.47276.40431.52141.05360.93564.42441.04170.5803Present1.45576.32811.53431.05890.92344.34361.05310.5816Error (%)−1.16−1.190.850.51−1.30−1.831.100.227Zenkour [26]1.50496.65471.44671.05890.95624.59710.99030.5834Present1.48696.57611.45951.06460.94334.51381.00180.5847Error (%)−1.20−1.180.880.54−1.35−1.811.160.228Zenkour [26]1.53436.89991.38291.06280.97504.76610.94660.5856Present1.51596.81971.39531.06870.96194.68100.95770.5870Error (%)−1.20−1.160.900.56−1.35−1.791.180.239Zenkour [26]1.56177.13831.32831.06620.99254.93030.90920.5875Present1.54357.05701.34041.07230.97944.84380.92000.5889Error (%)−1.17−1.140.910.57−1.32−1.751.190.2410Zenkour [26]1.58767.36891.28201.06941.00895.08900.87750.5894Present1.56987.28671.29351.07570.99625.00150.88790.5908Error (%)−1.12−1.120.900.59−1.26−1.721.180.24MetalZenkour [26]2.53272.89321.91031.28501.60701.99551.31210.7065Present2.53342.86961.91311.28711.60661.96971.31310.7071Error (%)0.03−0.820.150.16−0.02−1.290.080.08


## 5.4. Parameter Study

The proposed plate element MiQ4 is now applied to examine the static bending behavior of variable thickness FGM plates subjected to uniform load. The ingredients of the FGM plates are Al as metal and Al2O3 as ceramic. Young’s modulus is Em=70GPa for Al and Ec=380GPa for Al2O3, while Poisson’s ratios ν are constant and equal to 0.3, and the side-to-thickness ratio is a/h0=10. The material properties of the plate vary through the thickness direction by the power-law function. The nondimensional quantities are computed as the following formulae:(71)wmax∗=10Ech03q0a4wmax,σx∗=h0q0aσxa2,b2,z,τxy∗=h0q0aτxya,0,z.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 ofN×N=32×32 is handled in the rest of this study.


### 5.4.1. The Influence of the Power-Law Index

Firstly, we study the effects of the power-law indexp on the deflection of the variable thickness FGM plates. A square FGM plate with the side-to-thickness ratio of a/h0=10 and the varied parameter ζ=0.5 is considered. The 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 power-law 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.Table 13
The maximum nondimensional deflection of the FGM plates as a function ofp.
BCsThickness variationspCeramic0.51510MetalCCCCU plate0.16460.24960.32210.50230.56270.8935L plate0.40100.61240.79271.22111.35672.1762P plate0.28940.44090.57000.88150.98221.5707S plate0.55010.84161.09011.67421.85632.9854SCSCU plate0.24140.36720.47440.73630.82221.3104L plate0.59120.90441.17161.79971.99553.2085P plate0.44320.67690.87611.34901.49872.4054S plate0.77921.19331.54662.37072.62524.2286SSSSU plate0.46670.71600.92891.42061.56982.5334L plate1.16361.78952.32403.53983.90076.3151P plate0.82031.26041.63612.49532.75264.4518S plate1.62562.50123.24874.94305.44388.8222CFCFU plate0.36310.55290.71471.10681.23401.9708L plate1.64082.51813.26674.99115.51318.9048P plate1.20601.84942.39833.66924.05706.5452S plate2.21833.40624.41956.74557.446612.0390


### 5.4.2. The Eeffects of the Side-to-Thickness Ratio

Secondly, the effects of the side-to-thickness ratioa/h0 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. The 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 side-to-thickness ratio increases, the maximum nondimensional deflection of the plate decreases. In addition, the BCs affect strongly the deflection of the plates. The maximum deflection of the CCCC plate is smallest, while the maximum deflection of the CFCF plate is greatest.Table 14
The maximum nondimensional deflection of the FGM plates as a function ofa/h0.
BCsThickness variationsa/h051020304050CCCCU plate0.57450.41340.37070.36250.35960.3582L plate1.24451.01650.95650.94510.94100.9392P plate0.93390.73110.67750.66740.66380.6621S plate1.65321.39761.33041.31761.31311.3110SCSCU plate0.80540.60870.55640.54650.54300.5414L plate1.78061.50231.42961.41591.41101.4088P plate1.36911.12331.05891.04671.04241.0404S plate2.29571.98261.90091.88551.88011.8776SSSSU plate1.34241.19101.15311.14611.14371.1425L plate3.19092.97902.92592.91602.91262.9110P plate2.28642.09722.04982.04102.03792.0365S plate4.40334.16354.10314.09194.08804.0862CFCFU plate1.18830.91690.84510.83070.82530.8227L plate4.65064.18744.05954.03314.02324.0184P plate3.47503.07442.96662.94492.93702.9331S plate6.21055.66425.50845.47515.46245.4562


### 5.4.3. The Influence of the Thickness 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 ofp=1, the varying parameter of ζ=0.5, and the side-to-thickness ratio of a/h0=10 is examined. The dependents of the deflected shapes on the BCs and the type of thickness variation are exhibited in Figures 8–11 . The 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. The 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.Figure 8
Deflected shapes of the CCCC FGM plates depend on the thickness variations.Figure 9
Deflected shapes of the SCSC FGM plates depend on the thickness variations.Figure 10
Deflected shapes of the SSSS FGM plates depend on the thickness variations.Figure 11
Deflected shapes of the CFCF FGM plates depend on the thickness variations.Additionally, we investigate the dependence of the stresses of the variable thickness FGM plate on the type of thickness variation. The 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. The stresses are strongly affected by the variation of the thickness. The 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. The stress distribution also depends on the law of thickness variation. The 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.Figure 12
The distribution of normal and in-plane shear stresses through the thickness of the FGM plates.
(a)(b)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/h0=10 is considered. The 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. The 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.Figure 13
The influence of the thickness variation on the nondimensional maximum deflections of the FGM plates with different BCs: (a) CCCC plate, (b) SCSC plate, (c) SSSS plate, and (d) CFCF plate.
(a)(b)(c)(d)


## 5.4.1. The Influence of the Power-Law Index

Firstly, we study the effects of the power-law indexp on the deflection of the variable thickness FGM plates. A square FGM plate with the side-to-thickness ratio of a/h0=10 and the varied parameter ζ=0.5 is considered. The 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 power-law 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.Table 13
The maximum nondimensional deflection of the FGM plates as a function ofp.
BCsThickness variationspCeramic0.51510MetalCCCCU plate0.16460.24960.32210.50230.56270.8935L plate0.40100.61240.79271.22111.35672.1762P plate0.28940.44090.57000.88150.98221.5707S plate0.55010.84161.09011.67421.85632.9854SCSCU plate0.24140.36720.47440.73630.82221.3104L plate0.59120.90441.17161.79971.99553.2085P plate0.44320.67690.87611.34901.49872.4054S plate0.77921.19331.54662.37072.62524.2286SSSSU plate0.46670.71600.92891.42061.56982.5334L plate1.16361.78952.32403.53983.90076.3151P plate0.82031.26041.63612.49532.75264.4518S plate1.62562.50123.24874.94305.44388.8222CFCFU plate0.36310.55290.71471.10681.23401.9708L plate1.64082.51813.26674.99115.51318.9048P plate1.20601.84942.39833.66924.05706.5452S plate2.21833.40624.41956.74557.446612.0390


## 5.4.2. The Eeffects of the Side-to-Thickness Ratio

Secondly, the effects of the side-to-thickness ratioa/h0 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. The 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 side-to-thickness ratio increases, the maximum nondimensional deflection of the plate decreases. In addition, the BCs affect strongly the deflection of the plates. The maximum deflection of the CCCC plate is smallest, while the maximum deflection of the CFCF plate is greatest.Table 14
The maximum nondimensional deflection of the FGM plates as a function ofa/h0.
BCsThickness variationsa/h051020304050CCCCU plate0.57450.41340.37070.36250.35960.3582L plate1.24451.01650.95650.94510.94100.9392P plate0.93390.73110.67750.66740.66380.6621S plate1.65321.39761.33041.31761.31311.3110SCSCU plate0.80540.60870.55640.54650.54300.5414L plate1.78061.50231.42961.41591.41101.4088P plate1.36911.12331.05891.04671.04241.0404S plate2.29571.98261.90091.88551.88011.8776SSSSU plate1.34241.19101.15311.14611.14371.1425L plate3.19092.97902.92592.91602.91262.9110P plate2.28642.09722.04982.04102.03792.0365S plate4.40334.16354.10314.09194.08804.0862CFCFU plate1.18830.91690.84510.83070.82530.8227L plate4.65064.18744.05954.03314.02324.0184P plate3.47503.07442.96662.94492.93702.9331S plate6.21055.66425.50845.47515.46245.4562


## 5.4.3. The Influence of the Thickness 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 ofp=1, the varying parameter of ζ=0.5, and the side-to-thickness ratio of a/h0=10 is examined. The dependents of the deflected shapes on the BCs and the type of thickness variation are exhibited in Figures 8–11 . The 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. The 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.Figure 8
Deflected shapes of the CCCC FGM plates depend on the thickness variations.Figure 9
Deflected shapes of the SCSC FGM plates depend on the thickness variations.Figure 10
Deflected shapes of the SSSS FGM plates depend on the thickness variations.Figure 11
Deflected shapes of the CFCF FGM plates depend on the thickness variations.Additionally, we investigate the dependence of the stresses of the variable thickness FGM plate on the type of thickness variation. The 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. The stresses are strongly affected by the variation of the thickness. The 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. The stress distribution also depends on the law of thickness variation. The 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.Figure 12
The distribution of normal and in-plane shear stresses through the thickness of the FGM plates.
(a)(b)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/h0=10 is considered. The 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. The 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.Figure 13
The influence of the thickness variation on the nondimensional maximum deflections of the FGM plates with different BCs: (a) CCCC plate, (b) SCSC plate, (c) SSSS plate, and (d) CFCF plate.
(a)(b)(c)(d)


## 6. 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. Then 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:The proposed element MiQ4 utilizes the bilinear shape function which is simple and similar to that of the well-known classical Q4 element and allows us to calculate easily the numerical integration.The 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.The 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.The 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.The 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.



---

*Source: 6653350-2021-03-23.xml*