# A Spline Finite Point Method for Nonlinear Bending Analysis of FG Plates in Thermal Environments Based on a Locking-free Thin/Thick Plate Theory

Authors: Tian-Jing Mo; Jun Huang; Shuang-Bei Li; Hai Wu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

A spline finite point method (SFPM) based on a locking-free thin/thick plate theory, which is suitable for analysis of both thick and thin plates, is developed to study nonlinear bending behavior of functionally graded material (FGM) plates with different thickness in thermal environments. In the proposed method, one direction of the plate is discretized with a set of uniformly distributed spline nodes instead of meshes and the other direction is expressed with orthogonal functions determined by the boundary conditions. The displacements of the plate are constructed by the linear combination of orthogonal functions and cubic B-spline interpolation functions with high efficiency for modeling. The locking-free thin/thick plate theory used by the proposed method is based on the first-order shear deformation theory but takes the shear strains and displacements as basic unknowns. The material properties of the FG plate are assumed to vary along the thickness direction following the power function distribution. By comparing with several published research studies based on the finite element method (FEM), the correctness, efficiency, and generality of the new model are validated for rectangular plates. Moreover, uniform and nonlinear temperature rise conditions are discussed, respectively. The effect of the temperature distribution, in-plane temperature force, and elastic foundation on nonlinear bending under different parameters are discussed in detail.


---

## Body

## 1. Introduction

Functionally Graded Materials (FGMs), a new type of composite material developed in recent years, are obtained by mixing two (such as metal and ceramic) or several materials in a certain volume ratio [1]. Unlike traditional composite materials, the material composition of FGMs changes continuously along one direction without obvious delamination, so the physical properties mutation between different materials can be eliminated easily. FGMs are temperature dependent and usually used in high temperature environment and, moreover, can avoid stress concentration [2], reduce residual stress, and improve bond strength between material layers effectively. Therefore, FGMs have been widely applied in many fields nowadays, such as aerospace, medicine, mechanical, and construction engineering, which has high scientific research value and great potential for development [3].At present, there are a number of research studies on the deformation of FG plates and shells based on different theories, including the classical laminated plate theory (CLPT), the first-order shear deformation plate theory (FSDT), the higher-order shear deformation plate theory (HSDT), and the three-dimensional elasticity theory. BeikMohammadlou and EkhteraeiToussi [4] gave inclusive forms of governing equations for the elastoplastic buckling of FG plates based on the CLPT. Bagheri et al. [5] analyzed asymmetric thermal buckling behavior of elastically supported annular FG plates by the FSDT. Alshorbagy et al. [6] used the first-order shear deformation plate model to investigate uncoupled thermomechanical behavior of FG plates. Based on a novel integral first order shear deformation theory, Bousahla et al. [7] studied buckling and vibrational behavior of the composite beam armed with single-walled carbon nanotubes resting on Winkler-Pasternak elastic foundation. Using an nth-order shear deformation theory, thermomechanical flexural analysis of functionally graded material sandwich plates with P-FGM face sheets and E-FGM and symmetric S-FGM core is performed by Boussoula et al. [8]. Rahmani et al. [9] investigated bending and free vibration behavior of elastically supported functionally graded sandwich plates with different boundary conditions by an original novel high order shear theory. By employing a new type of quasi-3D hyperbolic shear deformation theory, Kaddari et al. [10] discussed statics and free vibration of functionally graded porous plates resting on elastic foundations, Zarga et al. [11] analyzed thermomechanical bending of FGM sandwich plates, and Mahmoudi et al. [12] studied the effect of the material gradation and elastic foundation on the deflections and stresses of functionally graded sandwich plate subjected to thermomechanical loading.It is worth noting that the above theories have their own advantages and limitations. The CLPT is simple and easy in mechanical analyzing, but will induce great errors when investigating thick plates because the transverse shear deformation of the plates is ignored, which is equivalent of over stiffness. Although the transverse shear deformation along the thickness direction is considered and the calculation accuracy is improved in the FSDT, it is necessary to use the shear correction factor in the calculation. The HSDT and three-dimensional elasticity theory have sufficient accuracy in the deformation analysis of plates and shells, but too many independent unknowns make the solution complicated and increase the computational cost. Moreover, the analysis of thin plate based on the first or higher-order plate theory has considered shear deformation may induce false shear strain, namely, shear-locking [13]. Therefore, some researchers devoted themselves to establish a general element and method to calculate both thick and thin plates [14–16], and a locking-free thin/thick plate theory was proposed. In this paper, this theory is based on the FSDT but takes the in-plane displacements, deflection, and shear strains as basic unknowns, which can effectively avoid shear locking phenomenon with simplicity and high generality.For the past decades, analytical, semianalytic, experimental, and numerical methods for FGM structures have been proposed [17–19]. The common numerical approaches, including the finite element methods (FEM) [20], differential quadrature methods (DQM) [21], and meshless methods, were quite popular in engineering because of their easy implementation with computers [22]. When using the traditional FEM to analyze nonlinear behaviors and large deformation problem, not only number of unknowns is large but also much effort needs to be expended in remeshing at each step of the problem solving [22]. On the contrary, when using the meshless method to analyze flat plate problem, the remeshing can be avoided naturally [23], and the boundary conditions processing is simple with less unknowns and high precision. Among those, the spline finite point method (SFPM) presented by Qin [24] is one of the spline meshless methods. As an active branch of modern function approximation, the B-spline function is widely used. Based on the B-spline function, orthogonal function, and minimum potential energy principle, the displacement modes of SFPM are constructed by the linear combination for the product of orthogonal functions and cubic B-spline interpolation functions. In this method, one direction of the flat structure is discretized by a set of uniformly distributed spline nodes, and the other one is expressed by orthogonal functions related to the corresponding boundary conditions. The equilibrium equation is established with the minimum potential energy principle. The SFPM has been proved more advantageous than the FEMs in analyzing the nonlinear behaviors of plates and beams. Li et al. [25] used a bidirectional B-spline QR method (BB-sQRM) to study the crack control of reinforced concrete beam, in which the plane model was discretized with a set of spline nodes along two directions. The results proved the high efficiency and low computational cost of this method. Liu et al. [26] investigated the free transverse vibration of FG Euler–Bernoulli beams through the SFPM and found that it is simple for modeling and boundary processing. Li et al. [27] derived the displacement and stress solutions of a piezoelectric laminated composite plate based on a bidirectional B-spline finite point method (BB-sFPM). The numerical results showed that the method was applied to the material parameter identification with a small loss of precision.Existing studies showed that, in complex temperature field, the heat conduction phenomenon and temperature dependence of FGMs can significantly affect the bending, buckling, and vibration of the plate [28], and the thermal load plays an important role in determining the mechanical behaviors of FGM structures. Besides, the FG plate generally deforms largely when subjected to the combined action of thermomechanical loads, and its mechanical behavior is usually nonlinear in this situation. For the nonlinear mechanical problems of FG plates in thermal environment, some relevant works were summarized in [29–32]. All these works proved that solving mechanical problems by simple linearization cannot meet the precision requirement in important structural calculations. The nonlinear analysis is very significant to the practical engineering application. Meanwhile, the temperature rising and lowing situations are varied, some papers have been published about different temperature distributions of FGM structures. Belbachir et al. [33] analyzed the response of antisymmetric cross-ply laminated plates subjected to a uniformly distributed nonlinear thermomechanical loading. Tounsi et al. [34] studied the static behavior of elastically supported advanced functionally graded ceramic-metal plates subjected to a nonlinear hygro-thermomechanical load. Dong and Li [35] proposed a linear temperature rise solution to analyze the nonlinear bending of FG rectangular plates. Farzad and Erfan [36] investigated the thermal effect on FG nanobeams subjected to two kinds of thermal loading, namely, linear and nonlinear temperature rise. Considering the uniform, linear, and nonlinear temperature distributions across the plate, Tebboune et al. [37] studied the thermoelastic buckling of thick FG plates. All the conclusions showed that the temperature distribution can affect the thermal deformation behavior of FGM structures.We found that there are less existing research studies for the nonlinear behaviors of FG plates using the locking-free thin/thick plate theory, and most of the models and analyses employing the traditional FEMs are based on the FSDT and HSDT. Thus, in order to enrich the theoretical system of FG plate and shell structures, increase the practicability of FGM and the diversity of research methods, and the SFPM is adopted in this paper to carry out the nonlinear bending analysis of FG plates on Winkler–Pasternak elastic foundations in thermal environments by the locking-free thin/thick plate theory. The correctness, efficiency, and generality of the model will be verified by comparing with the published studies. Meanwhile, several parameters, such as volume fraction index, thickness-length ratio, length-width ratio, temperature distribution, heating-up condition, and foundation stiffness, will be discussed about their influences on the in-plane temperature force and elastic foundation for the nonlinear bending of FG plates, which can supplement the test methods and provide further analytical means.


## 2. Spline Finite Point Method Formulation

Qin introduced the formulation of SFPM in his work [24]; accordingly, the n degree B-spine function can be defined by(1)φnx=∑k=0n+1−1kn+1kx−xk+nn!,where xk is the spline node, with(2)xk=k−n+12,n+1k=n+1!k!n+1−k!,x−xk+n=x−xkn,x−xk≥0,0,x−xk<0.The cubic B-spline function can be obtained withn = 3, i.e.,(3)φ3x=∑k=04−1k4kx−xk+33!=16x+2+3−4x+1+3+6x+3−4x−1+3+x−2+3.Suppose the length of a FG plate isa and the width is b, and its spline discretization model is shown in Figure 1.Figure 1
The spline finite point discretization model of a FGM square plate.For rectangular plates, there are no restrictions on the selection ofx-axis and y-axis by SFPM. As shown in Figure 1, the plate structure can be divided by the interval value ax in x direction, and the local coordinates of corresponding nodes are given by(4)xi=x0+iax,i=1,…,N,where N stands for the number of subdivisions in x direction, and(5)ax=aN,0=x0<x1<x2<⋯<xN−1<xN=a.The cubic B-spline functionφ3x is employed in this paper with the spline node xk=0 as its center, which can be expressed as(6)φ3x=16x+23,x∈−2,−1,x+23−4x+13,x∈−1,0,2−x3−41−x3,x∈0,1,2−x3x∈1,2,0,x≥2.The graphs of cubic B-spline function and its derivatives can by expressed as follows.Figure2 shows that the cubic B-spline function and its derivatives are piecewise, symmetrical, and compact.Figure 2
The cubic B-spline functions and its derivatives: (a) cubic B-spline functionϕ3x; (b) first-order derivative of the cubic B-spline function ϕ3′x; (c) second-order derivative of the cubic B-spline function ϕ3″x; (d) third-order derivative of the cubic B-spline function ϕ3″′x.
(a)(b)(c)(d)Thus, the displacement function inx direction of the plate can be represented by a linear combination of spline interpolation functions ϕix and generalized parameters ci. The displacement in y direction is expressed by the orthogonal function, which can meet the boundary conditions at both ends. For example, the transverse displacement can be obtained by(7)wx=∑i=−1N+1ciϕixZmy=ϕcZmy,in which(8)c=c−1c0c1⋯cN+1T,ϕ=ϕ−1ϕ0ϕ1⋯ϕN+1,where Zmy is the orthogonal function. For the two ends’ simply-supported plate in y direction, w=0 when y=0 and b, so the corresponding orthogonal function is suggested to be Zm=sinmπ/by. The spline interpolation functions ϕ can be constructed in the matrix form as follows:(9)ϕ=φ3iQs,(10)Qs=diagg,I,h,where I is a N−3×N−3 identity matrix, and(11)g=1−4101−12001,h=100−12101−41,φ3i=φ3xh+1φ3xhφ3xh−1⋯φ3xh−N−1.It is very simple to deal with the displacement boundary conditions using spline interpolation functions with the following properties:(12)ϕi0=0i≠−1,ϕia=0i≠N+1,ϕi′0=0i≠−1,0,ϕi′a=0i≠N,N+1.Therefore, the large number substitution method can be used to process the boundary conditions of the plate, that is, to substitute the corresponding row and column elements in the stiffness matrix with very large numbers in calculation. For instance, only elements corresponding toi=−1 and i=N+1 need to be substituted for the two ends’ simply supported plate in x direction.


## 3. Nonlinear Bending of FG Plates

### 3.1. Physical Properties of FGM

The properties of two materials are assumed to vary along the thickness of FG plates and follow the power function distribution; then, Young’s modulusEz, material density ρz, thermal expansion coefficient αz, and thermal conductivity κz are related to coordinate z in the thickness direction, which can be uniformly expressed as (Gupta and Talha [38])(13)Pz=Pt−Pb⋅h+2z2hk+Pb,where h denotes the thickness of the FG plate, k denotes the material volume fraction index, and t and b denote the upper and lower surface of the FGM, respectively.The temperature dependence of the FGM should be considered in a thermal environment. Since Young’s modulusE and thermal expansion coefficient α vary greatly with the temperature, while the Poisson ratio μ, material density ρ, and thermal conductivity κ generally are temperature independent, so only E and α in (13) are assumed to be related to the temperature and uniformly expressed by (Lal et al. [20] and Zhang and Zhou [32])(14)Pz,T=PtT−PbT⋅h+2z2hk+PbT,(15)PjT=P0P−1T−1+1+P1T1+P2T2+P3T3j=b,t,where T is the corresponding temperature and P−1, P0, P1, P2, and P3 are the temperature-dependent coefficients.


### 3.2. Temperature Distribution

In view of that FGMs are often applied to high temperature environment, for simulating the stress state and deformation accurately, and it is necessary to determine the temperature distribution inside the materials and consider the influence of the thermal stress. We assume that the upper and lower surface of FG plates are in a steady temperature state, and the temperature distributes uniformly along any longitudinal section and only vary along the thickness, which is independent of coordinates (x, y). There are two temperature rise conditions to be considered in this paper.


#### 3.2.1. Uniform Temperature Rise Condition

In this situation, the temperature is uniformly raised from the lower surface to upper surface. The temperature variationΔT=Tz−T0 is a constant, in which Tz is the temperature distribution function along the thickness direction (Javaheri and Eslami [39]) and T0 = 300 K is the reference temperature without considering the thermal strain.


#### 3.2.2. Nonlinear Temperature Rise Condition

The steady-state heat transfer equation of FG plates can by expressed as follows (Lal et al. [20]):(16)−ddzκzdTzdz=0,where κz is given in detail in (13) and Tz is the temperature distribution function along the thickness direction. Both of them are the functions of z.The upper and lower surface temperature of the plate are assumed to beTc and Tm, respectively; then, the temperature distribution function inside FG plates can be determined as (Lal et al. [20])(17)Tz=Tm+Tc−Tmηz,where(18)ηz=1C2z+h2h−κcmn+1κm2z+h2hk+1+κcm22n+1κm22z+h2hk+1−κcm33n+1κm32z+h2h3k+1+κcm44n+1κm42z+h2h4k+1−κcm55n+1κm52z+h2h5k+1,C=1−κcmn+1κm+κcm22n+1κm2−κcm33n+1κm3+κcm44n+1κm4−κcm55n+1κm5,in which κcm=κc−κm, κc and κm are the material thermal conductivities of the upper and lower surfaces, respectively.


### 3.3. Nonlinear Geometric Equation

Using the locking-free thin/thick plate theory for thick and thin plates, which takes the in-plane displacementsu and v, deflection w, and shear strains γxz and γyz as the independent unknowns, the shear strain conditions of thick plate (γ≠0) and thin plate (γ=0) can be met simultaneously. The specific formula derivation and description is given in Appendix B. Thus, the displacement field can be obtained by(19)u=u0−zφx,v=v0−zφy,w=w0,γxz=∂w∂x−φx,γyz=∂w∂y−φy,where u, v, and w are the total displacements, u0, v0, and w0 are the corresponding displacements of a point on the neutral surface, and φx and φy are the rotations about the y-axis and x-axis, respectively. In order to establish the stiffness equation easily, the strain equation of the plate is expressed as(20)ε∗εxpεypγxypεxbεybγxybγxzγyz=εx0εy0γxy000000+000zεx1zεy1zγxy1γxz1γyz1+εxNLεyNLγxyNL00000,in which the linear strains are(21)ε0=εx0εy0γxy0=u0,xv0,yu0,y+v0,x,ε1=εx1εy1γxy1=−φx,x−φy,y−φx,y−φy,x,γ=γxz1γyz1=w0,x−φxw0,y−φy,and the nonlinear strain is(22)εNL=εxNLεyNLγxyNL=12w0,x212w0,y2w0,xw0,y.


### 3.4. Constitutive Equations in Thermal Environment

In a thermal environment, the constitutive equations of FG plates on Winkler–Pasternak elastic foundation can be defined by (Lal et al. [20])(23)σ=Qfε+Pf−QfαΔT,Pf=K1w+K2∇2w,ST=αTQfTε+pΔT,where Qf is the elastic modulus of FG plates considering temperature dependence, α is the thermal expansion coefficient array, Pf is the foundation reaction per unit area, K1 and K2 are the linear normal stiffness and shear stiffness of the elastic foundation, respectively, ∇2=∂2/∂x2+∂2/∂y2 is the second order Laplace differential operator, ST is the entropy density function, and p is the specific heat capacity of materials. The matrix and vector expressions for the aforementioned physical quantities are defined as follows:(24)Qf=Qfε00Qfγ,Qfε=Ez,T1−μ21μ0μ10001−μ2,Qfγ=5Ez,T121+μ005Ez,T121+μ,α=αz,Tαz,T000T.The following stress expression is employed to analyze the geometrically nonlinear problems in the thermal environment:(25)σ∗=σx0σy0τxy000000+000σx1σy1τxy1τxz1τyz1+σxNLσyNLτxyNL00000−NxθNyθNxyθMxθMyθMxyθ00.For the plate analysis, stresses are usually rewritten as the generalized forces as follows:(26)σ¯∗=NMFs=AB0BD000Cε0+εNLε1γ−NθMθ0,in which(27)N=NxNyNxyT,M=MxMyMxyT,Fs=FsxFsyT,(28)Nθ=∫−h/2h/2QfeαΔTdz,Mθ=∫−h/2h/2QfeαzΔTdz,where Nθ and Mθ are the thermal force and thermal moment caused by temperature, respectively. [A], [B], [D], and [C] are the stretching stiffness, stretching-bending coupling stiffness, bending stiffness, and shearing stiffness, respectively, and the specific expressions are given by(29)A,B,D=∫−h/2h/2Qfe1,z,z2dz,C=∫−h/2h/2Qfγdz.Therefore, (26) can be simplified as(30)σ¯∗=D¯∗ε∗−F¯∗,in which(31)D¯∗=AB0BD000C,F¯∗=NθMθ0,ε∗=ε∗L+ε∗NL=ε0ε1γ+εNL00.


### 3.5. Spline Discretization Modeling

The FG plate is discretized uniformly inx direction, and the displacement functions can be constructed by the cubic B-spline interpolation functions and orthogonal functions as follows:(32)u0=∑m=1rϕamXmy,v0=∑m=1rϕbmYmy,w0=∑m=1rϕcmZmy,γxz=∑m=1rϕ′fmFmy,γyz=∑m=1rϕgmGmy,where r is the convergent series, ϕ is the cubic B-spline interpolation function defined by (9), and Xmy, Ymy, Zmy, Fmy, and Gmy are the orthogonal functions which meet the boundary conditions in y direction. All FG plates are simple supported at all ends (SSSS) in this study, so the boundary conditions are given by(33)w0=u0=φy=0,atx=0,a,w0=v0=φx=0,aty=0,b.In order to satisfy the boundary conditions at y = 0 andb, the corresponding orthogonal functions are suggested to be(34)Xmy=cosmπby,Ymy=sinmπby,Zmy=sinmπby,Fmy=sinmπby,Gmy=cosmπby.Thus, (32) can be expressed in the matrix form as follows:(35)V=Nγδ,in which V=u0v0wγxzγyzT are the displacements and shear strains of a point on the neutral surface. However, in the calculation, the boundary conditions can be processed simply by substituting (19) into (35), i.e.,(36)V0=N¯δ,in which V0=u0v0wφxφyT, φx and φy are the rotations of normal to the neutral surface about the y-axis and x-axis respectively, δ is the generalized displacement parameter, and both Nγ and N¯ are the displacement shape functions of the plate, i.e.,(37)Nγ=Nγ1Nγ2Nγ3⋯Nγr,N¯=N¯1N¯2N¯3⋯N¯r,(38)Nγm=ϕXm00000ϕYm00000ϕZm00000ϕ′Fm00000ϕGm,N¯m=ϕXm00000ϕYm00000ϕZm0000ϕ′Zm−ϕ′Fm000ϕZm′0−ϕGm,(39)δ=δ1Tδ2Tδ3T⋯δrT,(40)δm=amTbmTcmTfmTgmT,m=1,2,3,…,r.So, the linear strain in (31) can be denoted by the generalized displacement parameter as follows:(41)ε∗L=BLpBLbBLγTδ,where BLp=BLp1BLp2⋯BLpr, which is similar to BLb and BLγ, and(42)BLpm=ϕ′Xm00000ϕYm′000ϕXm′ϕ′Ym000,BLbm=00−ϕ″Zmϕ″Fm000−ϕZm″0ϕGm′00−2ϕ′Zm′ϕ′Fm′ϕ′Gm,BLγm=000ϕ′Fm00000ϕGm.The nonlinear strain can be constructed as(43)ε∗NL=12Q¯Q,in which(44)Q¯=∂w∂x00∂w∂y∂w∂y∂w∂x,Q=∂w∂x∂w∂y.Substitute (37) into (44), i.e.,(45)ε∗NL=BNLδ,in which(46)BNL=BNLP00T=12Q¯00TJ,Jm=00ϕ′Zm0000ϕZm′00,Q=∑m=1rJmcm=Jδ.Thus, (31) can be rewritten as(47)ε∗=BL+BNLδ.


### 3.6. Nonlinear Equilibrium Equation

Based on the virtual work principle, the virtual work equation of a thermomechanical FG plate is given by(48)dδTΨ=∬Adε∗Tσ¯∗dA−dδTF=0,where Ψ is the vector sum of internal and external forces, {F} denotes the transverse external force, dδ and ε are the virtual displacement and virtual strain, respectively, and σ¯∗=NxNyNxyMxMyMxyFsxFsyT are the generalized forces and given in detail in (27).The nonlinear equilibrium equation of the FG plate on Winkler–Pasternak elastic foundation is given by(49)Ψ=KS+KF+Kθδ−F−FR=0,where KS stands for the nonlinear bending stiffness, KF denotes the elastic foundation stiffness, Kθ is the thermal geometric stiffness, and FR is the thermal load caused by the environment temperature, which can be expressed as follows, respectively:(50)KS=∬ABLTD¯∗BLdA+2∬ABNLTD¯∗BLdA+∬ABLTD¯∗BNLdA+2∬ABNLTD¯∗BNLdA,(51)KF=∬ARTKlRdA,(52)Kl=K1000K2000K2,R=00ϕZm0000ϕ′Zm0000ϕZm′00,(53)Kθ=12∬AJTN¯θJdA,N¯θ=NxθNxyθNxyθNyθ,(54)FR=∫VBpTQfαΔT+BbTQfαzΔTdV.The differential ofΨ can be written as(55)dΨ=K0+Kσ∗+KNLdδ=KT∗dδ,in which(56)K0=∬ABLTD¯∗BLdA,(57)Kσ∗=∬JTNfJdxdy+∬RTKlRdxdy+12∬JTN¯θJdxdy,(58)KNL=∬A2BNLTD¯∗BL+2BLTD¯∗BNL+4BNLTD¯∗BNLdA,where Nf=NxNxyNxyNy and the third term of (57) is the effect of the in-plane temperature force. It is caused by the temperature change and boundary constraint and regarded as an axial load in this paper. The nonlinear equilibrium equation (50) can be solved by the Newton–Raphson method.


## 3.1. Physical Properties of FGM

The properties of two materials are assumed to vary along the thickness of FG plates and follow the power function distribution; then, Young’s modulusEz, material density ρz, thermal expansion coefficient αz, and thermal conductivity κz are related to coordinate z in the thickness direction, which can be uniformly expressed as (Gupta and Talha [38])(13)Pz=Pt−Pb⋅h+2z2hk+Pb,where h denotes the thickness of the FG plate, k denotes the material volume fraction index, and t and b denote the upper and lower surface of the FGM, respectively.The temperature dependence of the FGM should be considered in a thermal environment. Since Young’s modulusE and thermal expansion coefficient α vary greatly with the temperature, while the Poisson ratio μ, material density ρ, and thermal conductivity κ generally are temperature independent, so only E and α in (13) are assumed to be related to the temperature and uniformly expressed by (Lal et al. [20] and Zhang and Zhou [32])(14)Pz,T=PtT−PbT⋅h+2z2hk+PbT,(15)PjT=P0P−1T−1+1+P1T1+P2T2+P3T3j=b,t,where T is the corresponding temperature and P−1, P0, P1, P2, and P3 are the temperature-dependent coefficients.


## 3.2. Temperature Distribution

In view of that FGMs are often applied to high temperature environment, for simulating the stress state and deformation accurately, and it is necessary to determine the temperature distribution inside the materials and consider the influence of the thermal stress. We assume that the upper and lower surface of FG plates are in a steady temperature state, and the temperature distributes uniformly along any longitudinal section and only vary along the thickness, which is independent of coordinates (x, y). There are two temperature rise conditions to be considered in this paper.


### 3.2.1. Uniform Temperature Rise Condition

In this situation, the temperature is uniformly raised from the lower surface to upper surface. The temperature variationΔT=Tz−T0 is a constant, in which Tz is the temperature distribution function along the thickness direction (Javaheri and Eslami [39]) and T0 = 300 K is the reference temperature without considering the thermal strain.


### 3.2.2. Nonlinear Temperature Rise Condition

The steady-state heat transfer equation of FG plates can by expressed as follows (Lal et al. [20]):(16)−ddzκzdTzdz=0,where κz is given in detail in (13) and Tz is the temperature distribution function along the thickness direction. Both of them are the functions of z.The upper and lower surface temperature of the plate are assumed to beTc and Tm, respectively; then, the temperature distribution function inside FG plates can be determined as (Lal et al. [20])(17)Tz=Tm+Tc−Tmηz,where(18)ηz=1C2z+h2h−κcmn+1κm2z+h2hk+1+κcm22n+1κm22z+h2hk+1−κcm33n+1κm32z+h2h3k+1+κcm44n+1κm42z+h2h4k+1−κcm55n+1κm52z+h2h5k+1,C=1−κcmn+1κm+κcm22n+1κm2−κcm33n+1κm3+κcm44n+1κm4−κcm55n+1κm5,in which κcm=κc−κm, κc and κm are the material thermal conductivities of the upper and lower surfaces, respectively.


## 3.2.1. Uniform Temperature Rise Condition

In this situation, the temperature is uniformly raised from the lower surface to upper surface. The temperature variationΔT=Tz−T0 is a constant, in which Tz is the temperature distribution function along the thickness direction (Javaheri and Eslami [39]) and T0 = 300 K is the reference temperature without considering the thermal strain.


## 3.2.2. Nonlinear Temperature Rise Condition

The steady-state heat transfer equation of FG plates can by expressed as follows (Lal et al. [20]):(16)−ddzκzdTzdz=0,where κz is given in detail in (13) and Tz is the temperature distribution function along the thickness direction. Both of them are the functions of z.The upper and lower surface temperature of the plate are assumed to beTc and Tm, respectively; then, the temperature distribution function inside FG plates can be determined as (Lal et al. [20])(17)Tz=Tm+Tc−Tmηz,where(18)ηz=1C2z+h2h−κcmn+1κm2z+h2hk+1+κcm22n+1κm22z+h2hk+1−κcm33n+1κm32z+h2h3k+1+κcm44n+1κm42z+h2h4k+1−κcm55n+1κm52z+h2h5k+1,C=1−κcmn+1κm+κcm22n+1κm2−κcm33n+1κm3+κcm44n+1κm4−κcm55n+1κm5,in which κcm=κc−κm, κc and κm are the material thermal conductivities of the upper and lower surfaces, respectively.


## 3.3. Nonlinear Geometric Equation

Using the locking-free thin/thick plate theory for thick and thin plates, which takes the in-plane displacementsu and v, deflection w, and shear strains γxz and γyz as the independent unknowns, the shear strain conditions of thick plate (γ≠0) and thin plate (γ=0) can be met simultaneously. The specific formula derivation and description is given in Appendix B. Thus, the displacement field can be obtained by(19)u=u0−zφx,v=v0−zφy,w=w0,γxz=∂w∂x−φx,γyz=∂w∂y−φy,where u, v, and w are the total displacements, u0, v0, and w0 are the corresponding displacements of a point on the neutral surface, and φx and φy are the rotations about the y-axis and x-axis, respectively. In order to establish the stiffness equation easily, the strain equation of the plate is expressed as(20)ε∗εxpεypγxypεxbεybγxybγxzγyz=εx0εy0γxy000000+000zεx1zεy1zγxy1γxz1γyz1+εxNLεyNLγxyNL00000,in which the linear strains are(21)ε0=εx0εy0γxy0=u0,xv0,yu0,y+v0,x,ε1=εx1εy1γxy1=−φx,x−φy,y−φx,y−φy,x,γ=γxz1γyz1=w0,x−φxw0,y−φy,and the nonlinear strain is(22)εNL=εxNLεyNLγxyNL=12w0,x212w0,y2w0,xw0,y.


## 3.4. Constitutive Equations in Thermal Environment

In a thermal environment, the constitutive equations of FG plates on Winkler–Pasternak elastic foundation can be defined by (Lal et al. [20])(23)σ=Qfε+Pf−QfαΔT,Pf=K1w+K2∇2w,ST=αTQfTε+pΔT,where Qf is the elastic modulus of FG plates considering temperature dependence, α is the thermal expansion coefficient array, Pf is the foundation reaction per unit area, K1 and K2 are the linear normal stiffness and shear stiffness of the elastic foundation, respectively, ∇2=∂2/∂x2+∂2/∂y2 is the second order Laplace differential operator, ST is the entropy density function, and p is the specific heat capacity of materials. The matrix and vector expressions for the aforementioned physical quantities are defined as follows:(24)Qf=Qfε00Qfγ,Qfε=Ez,T1−μ21μ0μ10001−μ2,Qfγ=5Ez,T121+μ005Ez,T121+μ,α=αz,Tαz,T000T.The following stress expression is employed to analyze the geometrically nonlinear problems in the thermal environment:(25)σ∗=σx0σy0τxy000000+000σx1σy1τxy1τxz1τyz1+σxNLσyNLτxyNL00000−NxθNyθNxyθMxθMyθMxyθ00.For the plate analysis, stresses are usually rewritten as the generalized forces as follows:(26)σ¯∗=NMFs=AB0BD000Cε0+εNLε1γ−NθMθ0,in which(27)N=NxNyNxyT,M=MxMyMxyT,Fs=FsxFsyT,(28)Nθ=∫−h/2h/2QfeαΔTdz,Mθ=∫−h/2h/2QfeαzΔTdz,where Nθ and Mθ are the thermal force and thermal moment caused by temperature, respectively. [A], [B], [D], and [C] are the stretching stiffness, stretching-bending coupling stiffness, bending stiffness, and shearing stiffness, respectively, and the specific expressions are given by(29)A,B,D=∫−h/2h/2Qfe1,z,z2dz,C=∫−h/2h/2Qfγdz.Therefore, (26) can be simplified as(30)σ¯∗=D¯∗ε∗−F¯∗,in which(31)D¯∗=AB0BD000C,F¯∗=NθMθ0,ε∗=ε∗L+ε∗NL=ε0ε1γ+εNL00.


## 3.5. Spline Discretization Modeling

The FG plate is discretized uniformly inx direction, and the displacement functions can be constructed by the cubic B-spline interpolation functions and orthogonal functions as follows:(32)u0=∑m=1rϕamXmy,v0=∑m=1rϕbmYmy,w0=∑m=1rϕcmZmy,γxz=∑m=1rϕ′fmFmy,γyz=∑m=1rϕgmGmy,where r is the convergent series, ϕ is the cubic B-spline interpolation function defined by (9), and Xmy, Ymy, Zmy, Fmy, and Gmy are the orthogonal functions which meet the boundary conditions in y direction. All FG plates are simple supported at all ends (SSSS) in this study, so the boundary conditions are given by(33)w0=u0=φy=0,atx=0,a,w0=v0=φx=0,aty=0,b.In order to satisfy the boundary conditions at y = 0 andb, the corresponding orthogonal functions are suggested to be(34)Xmy=cosmπby,Ymy=sinmπby,Zmy=sinmπby,Fmy=sinmπby,Gmy=cosmπby.Thus, (32) can be expressed in the matrix form as follows:(35)V=Nγδ,in which V=u0v0wγxzγyzT are the displacements and shear strains of a point on the neutral surface. However, in the calculation, the boundary conditions can be processed simply by substituting (19) into (35), i.e.,(36)V0=N¯δ,in which V0=u0v0wφxφyT, φx and φy are the rotations of normal to the neutral surface about the y-axis and x-axis respectively, δ is the generalized displacement parameter, and both Nγ and N¯ are the displacement shape functions of the plate, i.e.,(37)Nγ=Nγ1Nγ2Nγ3⋯Nγr,N¯=N¯1N¯2N¯3⋯N¯r,(38)Nγm=ϕXm00000ϕYm00000ϕZm00000ϕ′Fm00000ϕGm,N¯m=ϕXm00000ϕYm00000ϕZm0000ϕ′Zm−ϕ′Fm000ϕZm′0−ϕGm,(39)δ=δ1Tδ2Tδ3T⋯δrT,(40)δm=amTbmTcmTfmTgmT,m=1,2,3,…,r.So, the linear strain in (31) can be denoted by the generalized displacement parameter as follows:(41)ε∗L=BLpBLbBLγTδ,where BLp=BLp1BLp2⋯BLpr, which is similar to BLb and BLγ, and(42)BLpm=ϕ′Xm00000ϕYm′000ϕXm′ϕ′Ym000,BLbm=00−ϕ″Zmϕ″Fm000−ϕZm″0ϕGm′00−2ϕ′Zm′ϕ′Fm′ϕ′Gm,BLγm=000ϕ′Fm00000ϕGm.The nonlinear strain can be constructed as(43)ε∗NL=12Q¯Q,in which(44)Q¯=∂w∂x00∂w∂y∂w∂y∂w∂x,Q=∂w∂x∂w∂y.Substitute (37) into (44), i.e.,(45)ε∗NL=BNLδ,in which(46)BNL=BNLP00T=12Q¯00TJ,Jm=00ϕ′Zm0000ϕZm′00,Q=∑m=1rJmcm=Jδ.Thus, (31) can be rewritten as(47)ε∗=BL+BNLδ.


## 3.6. Nonlinear Equilibrium Equation

Based on the virtual work principle, the virtual work equation of a thermomechanical FG plate is given by(48)dδTΨ=∬Adε∗Tσ¯∗dA−dδTF=0,where Ψ is the vector sum of internal and external forces, {F} denotes the transverse external force, dδ and ε are the virtual displacement and virtual strain, respectively, and σ¯∗=NxNyNxyMxMyMxyFsxFsyT are the generalized forces and given in detail in (27).The nonlinear equilibrium equation of the FG plate on Winkler–Pasternak elastic foundation is given by(49)Ψ=KS+KF+Kθδ−F−FR=0,where KS stands for the nonlinear bending stiffness, KF denotes the elastic foundation stiffness, Kθ is the thermal geometric stiffness, and FR is the thermal load caused by the environment temperature, which can be expressed as follows, respectively:(50)KS=∬ABLTD¯∗BLdA+2∬ABNLTD¯∗BLdA+∬ABLTD¯∗BNLdA+2∬ABNLTD¯∗BNLdA,(51)KF=∬ARTKlRdA,(52)Kl=K1000K2000K2,R=00ϕZm0000ϕ′Zm0000ϕZm′00,(53)Kθ=12∬AJTN¯θJdA,N¯θ=NxθNxyθNxyθNyθ,(54)FR=∫VBpTQfαΔT+BbTQfαzΔTdV.The differential ofΨ can be written as(55)dΨ=K0+Kσ∗+KNLdδ=KT∗dδ,in which(56)K0=∬ABLTD¯∗BLdA,(57)Kσ∗=∬JTNfJdxdy+∬RTKlRdxdy+12∬JTN¯θJdxdy,(58)KNL=∬A2BNLTD¯∗BL+2BLTD¯∗BNL+4BNLTD¯∗BNLdA,where Nf=NxNxyNxyNy and the third term of (57) is the effect of the in-plane temperature force. It is caused by the temperature change and boundary constraint and regarded as an axial load in this paper. The nonlinear equilibrium equation (50) can be solved by the Newton–Raphson method.


## 4. Results and Discussion

### 4.1. Material Properties

For the common FG plate, the upper surface is usually ceramic rich and the lower surface is metal rich. In this paper, two kinds of FG plates are studied. The first kind is the combination of Al/ZrO2 and Al/Zr2O3, which is temperature independent. The typical values of Young’s modulus E, the thermal expansion coefficient α, thermal conductivity κ, and the Poisson ratio μ are listed in Table 1 from Valizadeh et al. [40]. The second kind is SUS304/Si3N4, which is temperature dependent. The temperature-dependent coefficients are listed in Table 2 from Fu et al. [41]. The thermal material properties of SUS304 and Si3N4 under different temperatures are calculated in Table 3 according to (15) in Section 3.Table 1
Material properties of the first kind of FGM (from Valizadeh et al. [40]).
MaterialE (GPa)μκW/mkα×10−6/KAl700.320423.0ZrO22000.32.0910Zr2O31510.32.0910Al2O33800.310.47.2Table 2
Temperature-dependent coefficients for the material properties of the second kind of FGM (from Fu et al. [41]).
MaterialPropertiesP0P−1P1P2P3Si3N4E (Pa)348.43e + 90−3.070e − 42.160e − 7−8.964e − 11μ0.280000α1/K5.8723e − 609.095e-400κW/mK12.040000SUS304E (Pa)201.04e + 903.079e − 4−6.534e − 70μ0.280000α1/K12.330e − 608.086e − 400κW/mK9.190000Table 3
Thermal material properties of Si3N4 and SUS304 under different temperatures (according to (15) in Section 3).
MaterialTemperatureE (Pa)μα1/KκW/mKSi3N4T = 300 K322.27e + 90.287.48e − 69.19T = 400 K315.69e + 90.288.01e − 69.19T = 600 K303.92e + 90.289.08e − 69.19SUS304T = 300 K207.79e + 90.2816.32e − 612.04T = 400 K204.78e + 90.2815.32e − 612.04T = 600 K190.89e + 90.2818.31e − 612.04


### 4.2. Model Comparisons and Validations

#### 4.2.1. Correctness and Generality Studies of SFPM

The bending of Al/Al2O3 FGM square plate with different thicknesses and volume fraction indexes k are calculated and compared with some published research studies in Table 4. Numerical results of CLPT and FSDT are also included. In this example, the spline discretization level N × r = 12 × 3. The FG plates are subjected to transverse sinusoidal load q¯=q0sinπ/axsinπ/by, and q0 is the intensity of the transverse load at the plate center. The nondimensional transverse displacement w¯=10Ech3/q0a4w, where Ec = 380 GPa, and w is the center deflection.Table 4
Comparisons of the nondimensional center transverse displacement for Al/Al2O3 FGM square plates with different thicknesses and volume fraction indexes.
Volume fraction indexkSourcew¯a/h = 4a/h = 10a/h = 1001CLPT0.56230.56230.5623FSDT0.72910.58890.5625Carrera et al. [42]0.71710.58750.5625Neves et al. [43]0.69970.58450.5624Zenkour [44]0.70360.58480.5624Present0.70770.58600.56254CLPT0.82810.82810.8281FSDT1.11250.87360.8281Carrera et al. [42]1.15850.88210.8286Neves et al. [43]1.11780.87500.8286Zenkour [44]1.11140.87370.8286Present1.11660.87730.828510CLPT0.93540.93540.9354FSDT1.31780.99660.9360Carrera et al. [42]1.37451.00720.9361Neves et al. [43]1.34900.87500.8286Zenkour [44]1.32110.99770.9361Present1.32260.99640.9399It can be noticed that, in Table4, for the bending of FGM thin and thick plates, the present approach is very close to the result of FSDT and other published literature methods. It illustrates that the proposed SFPM based on the locking-free thin/thick plate theory is applicable to the plates of various thickness, and the model established in this paper is correct and general.


#### 4.2.2. Convergence studies of SFPM

The nonlinear bending of Al/ZrO2 FGM square plates with different discretization levels are calculated and compared with the NS-FEM results that are given by Nguyen-Xuan et al. [45]. The numerical results are listed in Table 5 with various volume fraction indexes. The deflection errors obtained by comparing the SFPM and NS-FEM results under different discretization levels are shown in Figure 3. In this example, the FG plates are subjected to a transverse uniformly pressure, the nondimensional center transverse displacement w¯=100wEmb3/121−μ2qa4, a/h=5, a/b=1, μ=0.3, and Em = 70 GPa.Table 5
Comparisons of the nondimensional center transverse displacement for Al/ZrO2 FGM square plates with different volume fraction indexes and discretization levels.
MethodDiscretization levelVolume fraction indexk0.512SFPM (N × r)4 × 30.23690.28350.33006 × 30.23520.27990.32378 × 30.23450.27670.319710 × 30.23400.27480.316412 × 30.23330.27310.313214 × 30.23280.27230.311916 × 30.23270.27200.311312 × 10.23730.27850.321012 × 50.23290.27260.312312 × 70.23280.27210.3114NS-FEM [45] (mesh)8 × 80.23320.27230.311612 × 120.23260.27160.310716 × 160.23240.27130.310420 × 200.23220.27110.310132 × 320.23210.27090.3099N stands for the number of spline subdivisions in x direction, and r is the convergent series.Figure 3
Convergence efficiency of the SFPM results compared with the NS-FEM results: (a) under different spline subdivisionsN; (b) under different convergent series r.
(a)(b)It was proved by Nguyen-Xuan et al. [45] that the NS-FEM can gain super-accurate and super-convergent stress solutions, and it can produce the stress solution that can be much more accurate than that of the FEM. As shown in Table 5 and Figure 3, with the increase of the spline subdivisions in x direction N and convergent series r, the results of SFPM converge to the NS-FEM results gradually. When the discretization level N × r is 12 × 3, 16 × 3, 20 × 3, 12 × 5, and 12 × 7, the maximum error is 0.80%, 0.29%, 0.23%, 0.52%, and 0.21%, respectively, which are all less than 1%. It means that the spline discretization level N × r = 12 × 3 can reach high precision. In order to improve the computational efficiency and reduce the calculation cost, the spline discretization level N × r = 12 × 3 is used in all of the following analyses, in which the number of unknowns is (12 + 3) × 3 × 5 = 225, and the results are able to meet the precision requirement. In Figure 3, we can see that the SFPM can nearly reach the same accurate and convergent as NS-FEM, and with less unknowns.


#### 4.2.3. Validation of the Elastically Supported FG Plate in Thermal Environments

The center deflection-load curves for SUS304/Si3N4 FGM square plates on Winkler–Pasternak elastic foundations in thermal environments are calculated and compared with the results from Lal et al. [20]. The FG plates subjected to uniformly distributed transverse loads q with a nonlinear temperature rise, the nondimensional pressure q∗=qb4/Emh4, ΔT=Tz−T0, and T0=300K. The values of Em with different temperature increments are listed in Table 3. The nondimensional linear normal and shear stiffness of the foundation are k1=K1a4/Emh3 and k2=K2a2/Emh3, respectively. The volume fraction index k = 2, b/a=1, and b/h=20. The model is shown in Figure 4 and the comparisons can be observed in Figure 5.Figure 4
The model of FGM square plates on Winkler–Pasternak elastic foundations.Figure 5
The nondimensional transverse deflection-load curves for SUS304/Si3N4 FG plates on Winkler–Pasternak elastic foundations subjected to uniform pressure and nonlinear temperature rise.As shown in Figure5, the results of this paper are in good agreement with the published results using different approaches. The model for FG plates on Winkler–Pasternak elastic foundations in the thermal environment with a nonlinear temperature rise is correct. Since the uniform temperature rise solution can be obtained by simplifying the nonlinear temperature rise function, the correctness of the model for FG plates with a uniform temperature rise also can be verified. As the number of the foundation parameters increases, the nonlinear center deflection decreases. This is because the foundation parameters make the stiffness of the plate increase. Meanwhile, the larger the temperature increment ΔT is, the larger the nonlinear center deflection is. This is because the thermal load FR increases with the increase of ΔT, which causes the nonlinear bending deformation to increase.


### 4.3. Parameter Studies

In this section, the following dimensionless parameters, nondimensional uniformly distributed transverse load (q∗), and foundation stiffness (k1 and k2) are defined as(59)q∗=qb4Emh4,k1=K1a4Emh3,k2=K2a2Emh3.There are three heating-up conditions to be considered (unless otherwise stated):C1:Tm = 400 K, Tc = 300 K; C2: Tm = 300 K, Tc = 400 K; and C3: Tm = 400 K, Tc = 400 K.And two kinds of situations for the in-plane temperature force:S1: the in-plane temperature force is considered and S2: the in-plane temperature force is neglected. S1 is the default situation unless otherwise stated.


#### 4.3.1. Discussion of Temperature Distribution

(1) Effect of Volume Fraction Index and Length-Thickness Ratio. In order to compare the nonlinear center deflection of FG plates under uniform and nonlinear temperature rise conditions, the uniformly distributed transverse load and elastic foundation are not considered temporarily. Taking the first kind of FGM as the research object, the center deflection for Al/ZrO2 square plates a=b=1m under two kinds of temperature rise conditions with different length-thickness ratios and volume fraction indexes are shown in Table 6.Table 6
The center deflection for Al/ZrO2 FGM square plates under two kinds of temperature rise conditions and C2 with different volume fraction index curves and length-thickness ratios.
Volume fraction indexa/hUniform temperature risew (10−3m)Nonlinear temperature risew (10−3m)Relative deviation (%)Ceramic (k = 0)101.12560.929217.45404.51283.857914.5110011.04789.906110.33k = 1101.23890.949823.33404.73443.859918.4710012.00379.970016.94k = 5101.66561.078135.27406.34684.397730.7110015.318511.451025.25Metal (k ⟶ +∞)101.99291.589020.27407.74686.511915.9410021.313117.127111.32Relative deviation is evaluated between the values in columns 3 and 4.Table6 presents that the center deflections with the uniform temperature rise are larger than that with the nonlinear temperature rise. It illustrates that the effect of the in-plane temperature force is weaker under the nonlinear temperature rise condition with C2. This is because the overall temperature of the plates is higher with the uniform temperature rise, so the influence of the thermal geometric stiffness Kθ is greater than the nonlinear bending stiffness Ks in this case. With the increase of volume fraction index k and thickness h of the plate, the relative deviation between two kinds of temperature rise conditions becomes bigger and bigger until the FGM tends to homogeneous the material again (k ⟶ +∞). It shows that the type of temperature rise condition has a greater influence on FGM thick plates.(2) Effect of FGM Combination and Width-Length Ratio. The nonlinear bending of two kinds of FG plates with different width-length ratios is calculated and compared in Table 7. In this example, the volume fraction index k = 2, b = 0.2 m, heating-up condition is C2, and two temperature rise conditions are considered.Table 7
The center deflection for Al/Zr2O3 and Si3N4/SUS304 FG plates under two kinds of temperature rise conditions and C2 with different width-length ratios h/b=1/20.
b/aUniform temperature risew (10−3m)Nonlinear temperature risew (10−3m)Relative deviation (%)Al/Zr2O310.42950.378111.972.50.75050.619017.5250.88620.645027.22Si3N4/SUS30410.20970.19168.622.50.34590.298413.7450.41620.311325.21Relative deviation is evaluated between the values in columns 3 and 4.Table7 shows that the relative deviation between two kinds of temperature rise conditions increases with width-length ratio b/a of the plate. It illustrates that the type of temperature distribution function has a greater influence on FGM rectangle plates. Besides, the relative deviation of Si3N4/SUS304 is smaller which presents that the effect of temperature rise condition is related to material combination. The temperature dependent FGM is less affected because the thermal material properties decrease with the increase of temperature, and the variation of nonlinear bending stiffness Ks is different under two kinds of temperature rise conditions. The following analyses only consider the nonlinear temperature rise condition.


#### 4.3.2. Discussion of In-Plane Temperature Force

(1) Effect of Volume Fraction Index and Heating-Up Condition. To compare the nonlinear center deflection of FG plates under in-plane temperature forces, the uniformly distributed transverse load and elastic foundation are still not considered. The Al/Zr2O3 square plates are under different temperature increments, a=b=1m and h/a=1/20. And the center deflection-volume fraction index curves are shown in Figure 6.Figure 6
The center deflection-volume fraction index curves for Al/Zr2O3 FGM square plates under in-plane temperature force with different temperature increments: (a) heating the lower surface of the plate separately; (b) heating the upper surface of the plate separately; (c) heating the upper and lower surfaces of the plate simultaneously.
(a)(b)(c)As shown in Figure6(a), when the lower surface of the plate is heated separately, the center deflection is negative and changes more greatly with the volume fraction index k = 0∼2, and it is basically unchanged after k > 2. This is because the FGM properties change greatly with k = 0∼2, after that the FGM tends to be homogeneous with k increasing. In Figure 6(b), the upper surface of the plate is heated separately. The bending deformation of the FG plate under the in-plane temperature force is positive and similar to that of the plate subjected to transverse mechanical loads in normal temperature. The center deflection reaches the smallest when k = 0 and increases with the increase of k. This is because the thermal expansion coefficient α and thermal conductivity κ of ceramic are smaller than those of metal. Therefore, the temperature change of the upper surface has little influence on the deflection. Figure 6(c) shows that the homogeneous plate (k = 0) does not deform and the center deflection is 0 when the upper and lower surfaces of the plate are heated simultaneously. This is because the temperature stays the same along the thickness and the thermal moment Mθ=0. However, for FG plates, the temperature distribution function and thermal expansion coefficient are not only the functions of thickness z but also related to k, so the deformation is similar to Figure 6(a).It can be observed from Figure6 that the center deflection increases with the temperature. This is because the in-plane temperature force is pressure when the plate is heated up, which makes the geometric stiffness matrix in Kσ∗ negative and the total stiffness decrease. Thus, the in-plane temperature force can soften the structure. It has a greater effect on FG plates than on homogeneous plates when the bending deformation is small, and the effect is related to the material properties and environment temperature.(2) Coupling Effect with Transverse Loads. The uniformly distributed transverse load is added into the consideration on the basis of (1). The nondimensional center deflection-load curves for FGM square plates with different volume fraction indexes and heating-up conditions are shown in Figure 7.Figure 7
Comparisons of the nondimensional center deflection-load curves for Al/Zr2O3 FGM square plates with different volume fraction indexes and heating-up conditions: (a) C1; (b) C2; (c) C3.
(a)(b)(c)It can be observed that the nonlinear center deflection increases with the volume fraction indexk in Figure 7. This is because the larger k is, the more metal and the less ceramic the plate has. The composition changes from pure ceramic to ceramic-metal composite, which makes the elastic modulus and stiffness of the FG plate decrease. Besides, under different heating-up conditions, the center deflection values with the same volume fraction index are much different from one another when q∗ is small, and they tend to be close with the increase of q∗. For example, for the FG plate with k = 2 under three heating-up conditions, the nondimensional center deflections are −0.112, 0.002, and -0.004, respectively, and the maximum deviation is 0.114 while q∗ = 0. However, when q∗ = 50, the center deflections are 0.948, 1.013, and 0.980, respectively, and the maximum deviation becomes 0.065. It illustrates that the influence of the environment temperature and in-plane temperature force decreases with the increase of nonlinear bending deformation of the plate. The large deflection is mainly controlled by the transverse distributed loads.(3) Effect of Length-Width Ratio and Thickness-Length Ratio. The nondimensional center-load curves for Al/Zr2O3 FGM square plates with various length-width ratios are shown in Figure 8. The nonlinear bending of two kinds of FG plates with different length-width ratios and thickness-length ratios are calculated and compared in Table 8. In this example, the volume fraction index k = 2.Figure 8
The nondimensional center deflection-load curves for Al/Zr2O3 FG plates with C2 and various length-width ratios.Table 8
The nondimensional center deflection for two kinds of FG plates under C1andq∗ = 25 with different length-width ratios and thickness-length ratios (k1 = k2 = 0).
a/bh/aAl/Zr2O3SUS304/Si3N4Al/ZrO2w/h (S1)w/h (S2)Relative deviation (%)w/h (S1)w/h (S2)Relative deviation (%)w/h (S1)w/h (S2)Relative deviation (%)0.51/80−0.6564−0.2506161.93−0.8646−0.3371156.48−0.6424−0.2378170.141/40−0.1436−0.0703104.27−0.2491−0.150665.41−0.1716−0.0795115.851/20−0.0209−0.010991.740.04960.041918.38−0.0329−0.016797.0111/80−0.5611−0.2446129.39−0.02810.1177123.87−0.5760−0.2291151.421/400.52890.407429.820.73490.625517.490.46890.349034.361/200.60980.501321.640.57210.52199.620.71380.553828.891.51/800.65300.488233.760.90580.720225.770.65420.467539.941/401.12810.948019.001.09270.978011.731.03510.850621.691/201.00850.914310.300.99970.94206.130.95140.813716.9221/801.02400.832023.081.31841.197610.090.66090.513028.831/401.13821.018611.741.28921.22625.141.15041.042310.371/201.17091.08437.991.26411.24411.611.07340.99378.02Relative deviation is evaluated between the values in S1 and S2.As shown in Figure8, the nonlinear center deflection values increase with the increase of length-width ratio a/b of the plate. The larger the uniformly distributed transverse load q∗ is, the more significant the influence of the length-width ratio is. Table 8 shows that, with the increase of length-width ratio a/b and thickness-length ratio h/a, the relative deviations between situation S1 and S2 become smaller. It illustrates that the thermal environment effect decreases gradually, and the in-plane temperature force has a greater influence on the thin plate with the small length-width ratio. Meanwhile, the calculation method employed in this paper is suitable for the analysis of both thick and thin plates, and the locking-free thin/thick plate theory has a wide applicability.


#### 4.3.3. Discussion of Elastic Foundation

(1) Effect of Transverse Distributed Load, Length-Width Ratio, and Heating-Up Condition. The nondimensional center deflection-load curves for Al/Zr2O3 FG plates are shown in Figure 9. The nonlinear bending of Al/ZrO2 plates under different length-width ratios and heating-up conditions are calculated and compared in Table 9. In this example, the plates are rested on Winkler–Pasternak elastic foundations with a=b=1m and volume fraction index k = 2.Figure 9
The nondimensional center deflection-load curves for elastically supported Al/Zr2O3 FG plates with (a) C1; (b) C2; (c) C3.
(a)(b)(c)Table 9
The center deflection for elastically supported Al/ZrO2 FG plates with different length-width ratios and heating-up conditions under q∗ = 25 and h/a=1/40.
a/bHeating-up conditions (Tm/Tc)k1=0k2=0w (10−3m)k1=10k2=10w (10−3m)k1=100k2=10w (10−3m)k1=100k2=100w (10−3m)Relative deviation I (%)Relative deviation II (%)Relative deviation III (%)0.5C110.740−1.297−0.524−0.097112.0859.6081.49C28.9881.0540.5320.09688.2749.5381.95C31.096−0.674−0.220−0.041161.5067.3681.361C111.7232.4161.8210.28579.3924.6384.35C220.0794.6693.2670.49076.7530.0385.00C310.1693.5502.4820.35165.0930.0885.861.5C122.1209.7307.9801.17056.0117.9985.34C223.92511.4029.3611.38352.3417.9085.23C322.04911.5299.2531.25047.7119.7486.492C128.76118.49414.5062.52735.7021.5682.58C229.19819.50417.3742.73733.2010.9284.25C328.90317.55615.5772.62439.2611.2783.15Relative deviation I is evaluated between the values in columns 3 and 4; relative deviation II is evaluated between the values in columns 4 and 5; Relative deviation III is evaluated between the values in columns 5 and 6.Figure9 shows that the stiffness of the FG plates increases and the center deflection decreases obviously after considering the effect of the elastic foundations, which indicates that the elastic foundations can stiffen the structure and weaken the bending behavior effectively. The center deflection decreases continuously with the increase of foundation parameter values. It can be observed that the shear stiffness K2 has a more obvious stiffening effect on the plates than the linear normal stiffness K1. As the uniformly distributed transverse load q∗ increases and the nonlinear deformation becomes larger and larger, the stiffening effect caused by K2 becomes more and more obvious.As shown in Table9, the relative deviations III are close to one another under different length-width ratios, which indicates that the influence of K2 is unrelated to the shape and size of the plate. However, the relative deviations I and II decrease with the increase of the length-width ratio. It illustrates that the shape and size can affect the stiffening effect caused by K1. The smaller the length-width ratio is, the larger the stiffening effect is. Besides, the relative deviations I under the three heating-up conditions become closer to one another as the length-width ratio increases. It shows that the influence of thermal environment decreases gradually, but the stiffening effect of the elastic foundation is unrelated to the heating-up condition.(2) Effect of Volume Fraction Index and Thickness-Length Ratio. The Al/Zr2O3 square plates a=b=1m are elastically supported. The center deflection with different volume fraction indexes and thickness-length ratios are shown in Table 10.Table 10
The center deflection for elastically supported Al/Zr2O3 FG plates with different volume fraction indexes and thickness-length ratios under q∗ = 25 and C2.
h/aVolume fraction indexk1=0k2=0w (10−3m)k1=10k2=10w (10−3m)k1=100k2=10w (10−3m)k1=100k2=100w (10−3m)Relative deviation I (%)Relative deviation II (%)Relative deviation III (%)1/80k = 0.56.2834.1602.6200.32933.7937.0287.44k = 28.7964.3962.7030.33150.0238.5187.75k = 39.0234.5052.7420.33250.0739.1387.891/50K = 0.517.5293.8262.6850.40078.1729.8285.10k = 219.4744.1682.8460.41378.6031.7285.49k = 319.9494.2102.8660.41478.9031.9285.551/25K = 0.525.8246.3024.5050.70875.6028.5184.28k = 226.6436.4344.5900.71775.8528.6684.38k = 327.2766.5504.6100.71975.9929.6284.401/10k = 0.593.08414.74810.5471.69884.1628.4983.90k = 295.22014.95310.6551.71084.3028.7483.95k = 399.58515.54011.0041.71684.4029.1984.41Relative deviation I is evaluated between the values in columns 3 and 4; Relative deviation II is evaluated between the values in columns 4 and 5; Relative deviation III is evaluated between the values in columns 5 and 6.Table10 shows that, for the thin plates (h/a = 1/80), the relative deviations I are smaller than that of the thick plates (h/a ≥ 1/50). The relative deviations II are larger than and the relative deviations III are close to those of the thick plates. It illustrates that the linear normal stiffness K1 has a greater influence on the thin plates, while the influence of shear stiffness K2 and the stiffening effect of the elastic foundation are unrelated to the thickness-length ratio. By comparing the three relative deviations with the same thickness-length ratio and different volume fraction indexes, the stiffening effect of the elastic foundation on FG plates is found to increase with k.


## 4.1. Material Properties

For the common FG plate, the upper surface is usually ceramic rich and the lower surface is metal rich. In this paper, two kinds of FG plates are studied. The first kind is the combination of Al/ZrO2 and Al/Zr2O3, which is temperature independent. The typical values of Young’s modulus E, the thermal expansion coefficient α, thermal conductivity κ, and the Poisson ratio μ are listed in Table 1 from Valizadeh et al. [40]. The second kind is SUS304/Si3N4, which is temperature dependent. The temperature-dependent coefficients are listed in Table 2 from Fu et al. [41]. The thermal material properties of SUS304 and Si3N4 under different temperatures are calculated in Table 3 according to (15) in Section 3.Table 1
Material properties of the first kind of FGM (from Valizadeh et al. [40]).
MaterialE (GPa)μκW/mkα×10−6/KAl700.320423.0ZrO22000.32.0910Zr2O31510.32.0910Al2O33800.310.47.2Table 2
Temperature-dependent coefficients for the material properties of the second kind of FGM (from Fu et al. [41]).
MaterialPropertiesP0P−1P1P2P3Si3N4E (Pa)348.43e + 90−3.070e − 42.160e − 7−8.964e − 11μ0.280000α1/K5.8723e − 609.095e-400κW/mK12.040000SUS304E (Pa)201.04e + 903.079e − 4−6.534e − 70μ0.280000α1/K12.330e − 608.086e − 400κW/mK9.190000Table 3
Thermal material properties of Si3N4 and SUS304 under different temperatures (according to (15) in Section 3).
MaterialTemperatureE (Pa)μα1/KκW/mKSi3N4T = 300 K322.27e + 90.287.48e − 69.19T = 400 K315.69e + 90.288.01e − 69.19T = 600 K303.92e + 90.289.08e − 69.19SUS304T = 300 K207.79e + 90.2816.32e − 612.04T = 400 K204.78e + 90.2815.32e − 612.04T = 600 K190.89e + 90.2818.31e − 612.04


## 4.2. Model Comparisons and Validations

### 4.2.1. Correctness and Generality Studies of SFPM

The bending of Al/Al2O3 FGM square plate with different thicknesses and volume fraction indexes k are calculated and compared with some published research studies in Table 4. Numerical results of CLPT and FSDT are also included. In this example, the spline discretization level N × r = 12 × 3. The FG plates are subjected to transverse sinusoidal load q¯=q0sinπ/axsinπ/by, and q0 is the intensity of the transverse load at the plate center. The nondimensional transverse displacement w¯=10Ech3/q0a4w, where Ec = 380 GPa, and w is the center deflection.Table 4
Comparisons of the nondimensional center transverse displacement for Al/Al2O3 FGM square plates with different thicknesses and volume fraction indexes.
Volume fraction indexkSourcew¯a/h = 4a/h = 10a/h = 1001CLPT0.56230.56230.5623FSDT0.72910.58890.5625Carrera et al. [42]0.71710.58750.5625Neves et al. [43]0.69970.58450.5624Zenkour [44]0.70360.58480.5624Present0.70770.58600.56254CLPT0.82810.82810.8281FSDT1.11250.87360.8281Carrera et al. [42]1.15850.88210.8286Neves et al. [43]1.11780.87500.8286Zenkour [44]1.11140.87370.8286Present1.11660.87730.828510CLPT0.93540.93540.9354FSDT1.31780.99660.9360Carrera et al. [42]1.37451.00720.9361Neves et al. [43]1.34900.87500.8286Zenkour [44]1.32110.99770.9361Present1.32260.99640.9399It can be noticed that, in Table4, for the bending of FGM thin and thick plates, the present approach is very close to the result of FSDT and other published literature methods. It illustrates that the proposed SFPM based on the locking-free thin/thick plate theory is applicable to the plates of various thickness, and the model established in this paper is correct and general.


### 4.2.2. Convergence studies of SFPM

The nonlinear bending of Al/ZrO2 FGM square plates with different discretization levels are calculated and compared with the NS-FEM results that are given by Nguyen-Xuan et al. [45]. The numerical results are listed in Table 5 with various volume fraction indexes. The deflection errors obtained by comparing the SFPM and NS-FEM results under different discretization levels are shown in Figure 3. In this example, the FG plates are subjected to a transverse uniformly pressure, the nondimensional center transverse displacement w¯=100wEmb3/121−μ2qa4, a/h=5, a/b=1, μ=0.3, and Em = 70 GPa.Table 5
Comparisons of the nondimensional center transverse displacement for Al/ZrO2 FGM square plates with different volume fraction indexes and discretization levels.
MethodDiscretization levelVolume fraction indexk0.512SFPM (N × r)4 × 30.23690.28350.33006 × 30.23520.27990.32378 × 30.23450.27670.319710 × 30.23400.27480.316412 × 30.23330.27310.313214 × 30.23280.27230.311916 × 30.23270.27200.311312 × 10.23730.27850.321012 × 50.23290.27260.312312 × 70.23280.27210.3114NS-FEM [45] (mesh)8 × 80.23320.27230.311612 × 120.23260.27160.310716 × 160.23240.27130.310420 × 200.23220.27110.310132 × 320.23210.27090.3099N stands for the number of spline subdivisions in x direction, and r is the convergent series.Figure 3
Convergence efficiency of the SFPM results compared with the NS-FEM results: (a) under different spline subdivisionsN; (b) under different convergent series r.
(a)(b)It was proved by Nguyen-Xuan et al. [45] that the NS-FEM can gain super-accurate and super-convergent stress solutions, and it can produce the stress solution that can be much more accurate than that of the FEM. As shown in Table 5 and Figure 3, with the increase of the spline subdivisions in x direction N and convergent series r, the results of SFPM converge to the NS-FEM results gradually. When the discretization level N × r is 12 × 3, 16 × 3, 20 × 3, 12 × 5, and 12 × 7, the maximum error is 0.80%, 0.29%, 0.23%, 0.52%, and 0.21%, respectively, which are all less than 1%. It means that the spline discretization level N × r = 12 × 3 can reach high precision. In order to improve the computational efficiency and reduce the calculation cost, the spline discretization level N × r = 12 × 3 is used in all of the following analyses, in which the number of unknowns is (12 + 3) × 3 × 5 = 225, and the results are able to meet the precision requirement. In Figure 3, we can see that the SFPM can nearly reach the same accurate and convergent as NS-FEM, and with less unknowns.


### 4.2.3. Validation of the Elastically Supported FG Plate in Thermal Environments

The center deflection-load curves for SUS304/Si3N4 FGM square plates on Winkler–Pasternak elastic foundations in thermal environments are calculated and compared with the results from Lal et al. [20]. The FG plates subjected to uniformly distributed transverse loads q with a nonlinear temperature rise, the nondimensional pressure q∗=qb4/Emh4, ΔT=Tz−T0, and T0=300K. The values of Em with different temperature increments are listed in Table 3. The nondimensional linear normal and shear stiffness of the foundation are k1=K1a4/Emh3 and k2=K2a2/Emh3, respectively. The volume fraction index k = 2, b/a=1, and b/h=20. The model is shown in Figure 4 and the comparisons can be observed in Figure 5.Figure 4
The model of FGM square plates on Winkler–Pasternak elastic foundations.Figure 5
The nondimensional transverse deflection-load curves for SUS304/Si3N4 FG plates on Winkler–Pasternak elastic foundations subjected to uniform pressure and nonlinear temperature rise.As shown in Figure5, the results of this paper are in good agreement with the published results using different approaches. The model for FG plates on Winkler–Pasternak elastic foundations in the thermal environment with a nonlinear temperature rise is correct. Since the uniform temperature rise solution can be obtained by simplifying the nonlinear temperature rise function, the correctness of the model for FG plates with a uniform temperature rise also can be verified. As the number of the foundation parameters increases, the nonlinear center deflection decreases. This is because the foundation parameters make the stiffness of the plate increase. Meanwhile, the larger the temperature increment ΔT is, the larger the nonlinear center deflection is. This is because the thermal load FR increases with the increase of ΔT, which causes the nonlinear bending deformation to increase.


## 4.2.1. Correctness and Generality Studies of SFPM

The bending of Al/Al2O3 FGM square plate with different thicknesses and volume fraction indexes k are calculated and compared with some published research studies in Table 4. Numerical results of CLPT and FSDT are also included. In this example, the spline discretization level N × r = 12 × 3. The FG plates are subjected to transverse sinusoidal load q¯=q0sinπ/axsinπ/by, and q0 is the intensity of the transverse load at the plate center. The nondimensional transverse displacement w¯=10Ech3/q0a4w, where Ec = 380 GPa, and w is the center deflection.Table 4
Comparisons of the nondimensional center transverse displacement for Al/Al2O3 FGM square plates with different thicknesses and volume fraction indexes.
Volume fraction indexkSourcew¯a/h = 4a/h = 10a/h = 1001CLPT0.56230.56230.5623FSDT0.72910.58890.5625Carrera et al. [42]0.71710.58750.5625Neves et al. [43]0.69970.58450.5624Zenkour [44]0.70360.58480.5624Present0.70770.58600.56254CLPT0.82810.82810.8281FSDT1.11250.87360.8281Carrera et al. [42]1.15850.88210.8286Neves et al. [43]1.11780.87500.8286Zenkour [44]1.11140.87370.8286Present1.11660.87730.828510CLPT0.93540.93540.9354FSDT1.31780.99660.9360Carrera et al. [42]1.37451.00720.9361Neves et al. [43]1.34900.87500.8286Zenkour [44]1.32110.99770.9361Present1.32260.99640.9399It can be noticed that, in Table4, for the bending of FGM thin and thick plates, the present approach is very close to the result of FSDT and other published literature methods. It illustrates that the proposed SFPM based on the locking-free thin/thick plate theory is applicable to the plates of various thickness, and the model established in this paper is correct and general.


## 4.2.2. Convergence studies of SFPM

The nonlinear bending of Al/ZrO2 FGM square plates with different discretization levels are calculated and compared with the NS-FEM results that are given by Nguyen-Xuan et al. [45]. The numerical results are listed in Table 5 with various volume fraction indexes. The deflection errors obtained by comparing the SFPM and NS-FEM results under different discretization levels are shown in Figure 3. In this example, the FG plates are subjected to a transverse uniformly pressure, the nondimensional center transverse displacement w¯=100wEmb3/121−μ2qa4, a/h=5, a/b=1, μ=0.3, and Em = 70 GPa.Table 5
Comparisons of the nondimensional center transverse displacement for Al/ZrO2 FGM square plates with different volume fraction indexes and discretization levels.
MethodDiscretization levelVolume fraction indexk0.512SFPM (N × r)4 × 30.23690.28350.33006 × 30.23520.27990.32378 × 30.23450.27670.319710 × 30.23400.27480.316412 × 30.23330.27310.313214 × 30.23280.27230.311916 × 30.23270.27200.311312 × 10.23730.27850.321012 × 50.23290.27260.312312 × 70.23280.27210.3114NS-FEM [45] (mesh)8 × 80.23320.27230.311612 × 120.23260.27160.310716 × 160.23240.27130.310420 × 200.23220.27110.310132 × 320.23210.27090.3099N stands for the number of spline subdivisions in x direction, and r is the convergent series.Figure 3
Convergence efficiency of the SFPM results compared with the NS-FEM results: (a) under different spline subdivisionsN; (b) under different convergent series r.
(a)(b)It was proved by Nguyen-Xuan et al. [45] that the NS-FEM can gain super-accurate and super-convergent stress solutions, and it can produce the stress solution that can be much more accurate than that of the FEM. As shown in Table 5 and Figure 3, with the increase of the spline subdivisions in x direction N and convergent series r, the results of SFPM converge to the NS-FEM results gradually. When the discretization level N × r is 12 × 3, 16 × 3, 20 × 3, 12 × 5, and 12 × 7, the maximum error is 0.80%, 0.29%, 0.23%, 0.52%, and 0.21%, respectively, which are all less than 1%. It means that the spline discretization level N × r = 12 × 3 can reach high precision. In order to improve the computational efficiency and reduce the calculation cost, the spline discretization level N × r = 12 × 3 is used in all of the following analyses, in which the number of unknowns is (12 + 3) × 3 × 5 = 225, and the results are able to meet the precision requirement. In Figure 3, we can see that the SFPM can nearly reach the same accurate and convergent as NS-FEM, and with less unknowns.


## 4.2.3. Validation of the Elastically Supported FG Plate in Thermal Environments

The center deflection-load curves for SUS304/Si3N4 FGM square plates on Winkler–Pasternak elastic foundations in thermal environments are calculated and compared with the results from Lal et al. [20]. The FG plates subjected to uniformly distributed transverse loads q with a nonlinear temperature rise, the nondimensional pressure q∗=qb4/Emh4, ΔT=Tz−T0, and T0=300K. The values of Em with different temperature increments are listed in Table 3. The nondimensional linear normal and shear stiffness of the foundation are k1=K1a4/Emh3 and k2=K2a2/Emh3, respectively. The volume fraction index k = 2, b/a=1, and b/h=20. The model is shown in Figure 4 and the comparisons can be observed in Figure 5.Figure 4
The model of FGM square plates on Winkler–Pasternak elastic foundations.Figure 5
The nondimensional transverse deflection-load curves for SUS304/Si3N4 FG plates on Winkler–Pasternak elastic foundations subjected to uniform pressure and nonlinear temperature rise.As shown in Figure5, the results of this paper are in good agreement with the published results using different approaches. The model for FG plates on Winkler–Pasternak elastic foundations in the thermal environment with a nonlinear temperature rise is correct. Since the uniform temperature rise solution can be obtained by simplifying the nonlinear temperature rise function, the correctness of the model for FG plates with a uniform temperature rise also can be verified. As the number of the foundation parameters increases, the nonlinear center deflection decreases. This is because the foundation parameters make the stiffness of the plate increase. Meanwhile, the larger the temperature increment ΔT is, the larger the nonlinear center deflection is. This is because the thermal load FR increases with the increase of ΔT, which causes the nonlinear bending deformation to increase.


## 4.3. Parameter Studies

In this section, the following dimensionless parameters, nondimensional uniformly distributed transverse load (q∗), and foundation stiffness (k1 and k2) are defined as(59)q∗=qb4Emh4,k1=K1a4Emh3,k2=K2a2Emh3.There are three heating-up conditions to be considered (unless otherwise stated):C1:Tm = 400 K, Tc = 300 K; C2: Tm = 300 K, Tc = 400 K; and C3: Tm = 400 K, Tc = 400 K.And two kinds of situations for the in-plane temperature force:S1: the in-plane temperature force is considered and S2: the in-plane temperature force is neglected. S1 is the default situation unless otherwise stated.


### 4.3.1. Discussion of Temperature Distribution

(1) Effect of Volume Fraction Index and Length-Thickness Ratio. In order to compare the nonlinear center deflection of FG plates under uniform and nonlinear temperature rise conditions, the uniformly distributed transverse load and elastic foundation are not considered temporarily. Taking the first kind of FGM as the research object, the center deflection for Al/ZrO2 square plates a=b=1m under two kinds of temperature rise conditions with different length-thickness ratios and volume fraction indexes are shown in Table 6.Table 6
The center deflection for Al/ZrO2 FGM square plates under two kinds of temperature rise conditions and C2 with different volume fraction index curves and length-thickness ratios.
Volume fraction indexa/hUniform temperature risew (10−3m)Nonlinear temperature risew (10−3m)Relative deviation (%)Ceramic (k = 0)101.12560.929217.45404.51283.857914.5110011.04789.906110.33k = 1101.23890.949823.33404.73443.859918.4710012.00379.970016.94k = 5101.66561.078135.27406.34684.397730.7110015.318511.451025.25Metal (k ⟶ +∞)101.99291.589020.27407.74686.511915.9410021.313117.127111.32Relative deviation is evaluated between the values in columns 3 and 4.Table6 presents that the center deflections with the uniform temperature rise are larger than that with the nonlinear temperature rise. It illustrates that the effect of the in-plane temperature force is weaker under the nonlinear temperature rise condition with C2. This is because the overall temperature of the plates is higher with the uniform temperature rise, so the influence of the thermal geometric stiffness Kθ is greater than the nonlinear bending stiffness Ks in this case. With the increase of volume fraction index k and thickness h of the plate, the relative deviation between two kinds of temperature rise conditions becomes bigger and bigger until the FGM tends to homogeneous the material again (k ⟶ +∞). It shows that the type of temperature rise condition has a greater influence on FGM thick plates.(2) Effect of FGM Combination and Width-Length Ratio. The nonlinear bending of two kinds of FG plates with different width-length ratios is calculated and compared in Table 7. In this example, the volume fraction index k = 2, b = 0.2 m, heating-up condition is C2, and two temperature rise conditions are considered.Table 7
The center deflection for Al/Zr2O3 and Si3N4/SUS304 FG plates under two kinds of temperature rise conditions and C2 with different width-length ratios h/b=1/20.
b/aUniform temperature risew (10−3m)Nonlinear temperature risew (10−3m)Relative deviation (%)Al/Zr2O310.42950.378111.972.50.75050.619017.5250.88620.645027.22Si3N4/SUS30410.20970.19168.622.50.34590.298413.7450.41620.311325.21Relative deviation is evaluated between the values in columns 3 and 4.Table7 shows that the relative deviation between two kinds of temperature rise conditions increases with width-length ratio b/a of the plate. It illustrates that the type of temperature distribution function has a greater influence on FGM rectangle plates. Besides, the relative deviation of Si3N4/SUS304 is smaller which presents that the effect of temperature rise condition is related to material combination. The temperature dependent FGM is less affected because the thermal material properties decrease with the increase of temperature, and the variation of nonlinear bending stiffness Ks is different under two kinds of temperature rise conditions. The following analyses only consider the nonlinear temperature rise condition.


### 4.3.2. Discussion of In-Plane Temperature Force

(1) Effect of Volume Fraction Index and Heating-Up Condition. To compare the nonlinear center deflection of FG plates under in-plane temperature forces, the uniformly distributed transverse load and elastic foundation are still not considered. The Al/Zr2O3 square plates are under different temperature increments, a=b=1m and h/a=1/20. And the center deflection-volume fraction index curves are shown in Figure 6.Figure 6
The center deflection-volume fraction index curves for Al/Zr2O3 FGM square plates under in-plane temperature force with different temperature increments: (a) heating the lower surface of the plate separately; (b) heating the upper surface of the plate separately; (c) heating the upper and lower surfaces of the plate simultaneously.
(a)(b)(c)As shown in Figure6(a), when the lower surface of the plate is heated separately, the center deflection is negative and changes more greatly with the volume fraction index k = 0∼2, and it is basically unchanged after k > 2. This is because the FGM properties change greatly with k = 0∼2, after that the FGM tends to be homogeneous with k increasing. In Figure 6(b), the upper surface of the plate is heated separately. The bending deformation of the FG plate under the in-plane temperature force is positive and similar to that of the plate subjected to transverse mechanical loads in normal temperature. The center deflection reaches the smallest when k = 0 and increases with the increase of k. This is because the thermal expansion coefficient α and thermal conductivity κ of ceramic are smaller than those of metal. Therefore, the temperature change of the upper surface has little influence on the deflection. Figure 6(c) shows that the homogeneous plate (k = 0) does not deform and the center deflection is 0 when the upper and lower surfaces of the plate are heated simultaneously. This is because the temperature stays the same along the thickness and the thermal moment Mθ=0. However, for FG plates, the temperature distribution function and thermal expansion coefficient are not only the functions of thickness z but also related to k, so the deformation is similar to Figure 6(a).It can be observed from Figure6 that the center deflection increases with the temperature. This is because the in-plane temperature force is pressure when the plate is heated up, which makes the geometric stiffness matrix in Kσ∗ negative and the total stiffness decrease. Thus, the in-plane temperature force can soften the structure. It has a greater effect on FG plates than on homogeneous plates when the bending deformation is small, and the effect is related to the material properties and environment temperature.(2) Coupling Effect with Transverse Loads. The uniformly distributed transverse load is added into the consideration on the basis of (1). The nondimensional center deflection-load curves for FGM square plates with different volume fraction indexes and heating-up conditions are shown in Figure 7.Figure 7
Comparisons of the nondimensional center deflection-load curves for Al/Zr2O3 FGM square plates with different volume fraction indexes and heating-up conditions: (a) C1; (b) C2; (c) C3.
(a)(b)(c)It can be observed that the nonlinear center deflection increases with the volume fraction indexk in Figure 7. This is because the larger k is, the more metal and the less ceramic the plate has. The composition changes from pure ceramic to ceramic-metal composite, which makes the elastic modulus and stiffness of the FG plate decrease. Besides, under different heating-up conditions, the center deflection values with the same volume fraction index are much different from one another when q∗ is small, and they tend to be close with the increase of q∗. For example, for the FG plate with k = 2 under three heating-up conditions, the nondimensional center deflections are −0.112, 0.002, and -0.004, respectively, and the maximum deviation is 0.114 while q∗ = 0. However, when q∗ = 50, the center deflections are 0.948, 1.013, and 0.980, respectively, and the maximum deviation becomes 0.065. It illustrates that the influence of the environment temperature and in-plane temperature force decreases with the increase of nonlinear bending deformation of the plate. The large deflection is mainly controlled by the transverse distributed loads.(3) Effect of Length-Width Ratio and Thickness-Length Ratio. The nondimensional center-load curves for Al/Zr2O3 FGM square plates with various length-width ratios are shown in Figure 8. The nonlinear bending of two kinds of FG plates with different length-width ratios and thickness-length ratios are calculated and compared in Table 8. In this example, the volume fraction index k = 2.Figure 8
The nondimensional center deflection-load curves for Al/Zr2O3 FG plates with C2 and various length-width ratios.Table 8
The nondimensional center deflection for two kinds of FG plates under C1andq∗ = 25 with different length-width ratios and thickness-length ratios (k1 = k2 = 0).
a/bh/aAl/Zr2O3SUS304/Si3N4Al/ZrO2w/h (S1)w/h (S2)Relative deviation (%)w/h (S1)w/h (S2)Relative deviation (%)w/h (S1)w/h (S2)Relative deviation (%)0.51/80−0.6564−0.2506161.93−0.8646−0.3371156.48−0.6424−0.2378170.141/40−0.1436−0.0703104.27−0.2491−0.150665.41−0.1716−0.0795115.851/20−0.0209−0.010991.740.04960.041918.38−0.0329−0.016797.0111/80−0.5611−0.2446129.39−0.02810.1177123.87−0.5760−0.2291151.421/400.52890.407429.820.73490.625517.490.46890.349034.361/200.60980.501321.640.57210.52199.620.71380.553828.891.51/800.65300.488233.760.90580.720225.770.65420.467539.941/401.12810.948019.001.09270.978011.731.03510.850621.691/201.00850.914310.300.99970.94206.130.95140.813716.9221/801.02400.832023.081.31841.197610.090.66090.513028.831/401.13821.018611.741.28921.22625.141.15041.042310.371/201.17091.08437.991.26411.24411.611.07340.99378.02Relative deviation is evaluated between the values in S1 and S2.As shown in Figure8, the nonlinear center deflection values increase with the increase of length-width ratio a/b of the plate. The larger the uniformly distributed transverse load q∗ is, the more significant the influence of the length-width ratio is. Table 8 shows that, with the increase of length-width ratio a/b and thickness-length ratio h/a, the relative deviations between situation S1 and S2 become smaller. It illustrates that the thermal environment effect decreases gradually, and the in-plane temperature force has a greater influence on the thin plate with the small length-width ratio. Meanwhile, the calculation method employed in this paper is suitable for the analysis of both thick and thin plates, and the locking-free thin/thick plate theory has a wide applicability.


### 4.3.3. Discussion of Elastic Foundation

(1) Effect of Transverse Distributed Load, Length-Width Ratio, and Heating-Up Condition. The nondimensional center deflection-load curves for Al/Zr2O3 FG plates are shown in Figure 9. The nonlinear bending of Al/ZrO2 plates under different length-width ratios and heating-up conditions are calculated and compared in Table 9. In this example, the plates are rested on Winkler–Pasternak elastic foundations with a=b=1m and volume fraction index k = 2.Figure 9
The nondimensional center deflection-load curves for elastically supported Al/Zr2O3 FG plates with (a) C1; (b) C2; (c) C3.
(a)(b)(c)Table 9
The center deflection for elastically supported Al/ZrO2 FG plates with different length-width ratios and heating-up conditions under q∗ = 25 and h/a=1/40.
a/bHeating-up conditions (Tm/Tc)k1=0k2=0w (10−3m)k1=10k2=10w (10−3m)k1=100k2=10w (10−3m)k1=100k2=100w (10−3m)Relative deviation I (%)Relative deviation II (%)Relative deviation III (%)0.5C110.740−1.297−0.524−0.097112.0859.6081.49C28.9881.0540.5320.09688.2749.5381.95C31.096−0.674−0.220−0.041161.5067.3681.361C111.7232.4161.8210.28579.3924.6384.35C220.0794.6693.2670.49076.7530.0385.00C310.1693.5502.4820.35165.0930.0885.861.5C122.1209.7307.9801.17056.0117.9985.34C223.92511.4029.3611.38352.3417.9085.23C322.04911.5299.2531.25047.7119.7486.492C128.76118.49414.5062.52735.7021.5682.58C229.19819.50417.3742.73733.2010.9284.25C328.90317.55615.5772.62439.2611.2783.15Relative deviation I is evaluated between the values in columns 3 and 4; relative deviation II is evaluated between the values in columns 4 and 5; Relative deviation III is evaluated between the values in columns 5 and 6.Figure9 shows that the stiffness of the FG plates increases and the center deflection decreases obviously after considering the effect of the elastic foundations, which indicates that the elastic foundations can stiffen the structure and weaken the bending behavior effectively. The center deflection decreases continuously with the increase of foundation parameter values. It can be observed that the shear stiffness K2 has a more obvious stiffening effect on the plates than the linear normal stiffness K1. As the uniformly distributed transverse load q∗ increases and the nonlinear deformation becomes larger and larger, the stiffening effect caused by K2 becomes more and more obvious.As shown in Table9, the relative deviations III are close to one another under different length-width ratios, which indicates that the influence of K2 is unrelated to the shape and size of the plate. However, the relative deviations I and II decrease with the increase of the length-width ratio. It illustrates that the shape and size can affect the stiffening effect caused by K1. The smaller the length-width ratio is, the larger the stiffening effect is. Besides, the relative deviations I under the three heating-up conditions become closer to one another as the length-width ratio increases. It shows that the influence of thermal environment decreases gradually, but the stiffening effect of the elastic foundation is unrelated to the heating-up condition.(2) Effect of Volume Fraction Index and Thickness-Length Ratio. The Al/Zr2O3 square plates a=b=1m are elastically supported. The center deflection with different volume fraction indexes and thickness-length ratios are shown in Table 10.Table 10
The center deflection for elastically supported Al/Zr2O3 FG plates with different volume fraction indexes and thickness-length ratios under q∗ = 25 and C2.
h/aVolume fraction indexk1=0k2=0w (10−3m)k1=10k2=10w (10−3m)k1=100k2=10w (10−3m)k1=100k2=100w (10−3m)Relative deviation I (%)Relative deviation II (%)Relative deviation III (%)1/80k = 0.56.2834.1602.6200.32933.7937.0287.44k = 28.7964.3962.7030.33150.0238.5187.75k = 39.0234.5052.7420.33250.0739.1387.891/50K = 0.517.5293.8262.6850.40078.1729.8285.10k = 219.4744.1682.8460.41378.6031.7285.49k = 319.9494.2102.8660.41478.9031.9285.551/25K = 0.525.8246.3024.5050.70875.6028.5184.28k = 226.6436.4344.5900.71775.8528.6684.38k = 327.2766.5504.6100.71975.9929.6284.401/10k = 0.593.08414.74810.5471.69884.1628.4983.90k = 295.22014.95310.6551.71084.3028.7483.95k = 399.58515.54011.0041.71684.4029.1984.41Relative deviation I is evaluated between the values in columns 3 and 4; Relative deviation II is evaluated between the values in columns 4 and 5; Relative deviation III is evaluated between the values in columns 5 and 6.Table10 shows that, for the thin plates (h/a = 1/80), the relative deviations I are smaller than that of the thick plates (h/a ≥ 1/50). The relative deviations II are larger than and the relative deviations III are close to those of the thick plates. It illustrates that the linear normal stiffness K1 has a greater influence on the thin plates, while the influence of shear stiffness K2 and the stiffening effect of the elastic foundation are unrelated to the thickness-length ratio. By comparing the three relative deviations with the same thickness-length ratio and different volume fraction indexes, the stiffening effect of the elastic foundation on FG plates is found to increase with k.


## 4.3.1. Discussion of Temperature Distribution

(1) Effect of Volume Fraction Index and Length-Thickness Ratio. In order to compare the nonlinear center deflection of FG plates under uniform and nonlinear temperature rise conditions, the uniformly distributed transverse load and elastic foundation are not considered temporarily. Taking the first kind of FGM as the research object, the center deflection for Al/ZrO2 square plates a=b=1m under two kinds of temperature rise conditions with different length-thickness ratios and volume fraction indexes are shown in Table 6.Table 6
The center deflection for Al/ZrO2 FGM square plates under two kinds of temperature rise conditions and C2 with different volume fraction index curves and length-thickness ratios.
Volume fraction indexa/hUniform temperature risew (10−3m)Nonlinear temperature risew (10−3m)Relative deviation (%)Ceramic (k = 0)101.12560.929217.45404.51283.857914.5110011.04789.906110.33k = 1101.23890.949823.33404.73443.859918.4710012.00379.970016.94k = 5101.66561.078135.27406.34684.397730.7110015.318511.451025.25Metal (k ⟶ +∞)101.99291.589020.27407.74686.511915.9410021.313117.127111.32Relative deviation is evaluated between the values in columns 3 and 4.Table6 presents that the center deflections with the uniform temperature rise are larger than that with the nonlinear temperature rise. It illustrates that the effect of the in-plane temperature force is weaker under the nonlinear temperature rise condition with C2. This is because the overall temperature of the plates is higher with the uniform temperature rise, so the influence of the thermal geometric stiffness Kθ is greater than the nonlinear bending stiffness Ks in this case. With the increase of volume fraction index k and thickness h of the plate, the relative deviation between two kinds of temperature rise conditions becomes bigger and bigger until the FGM tends to homogeneous the material again (k ⟶ +∞). It shows that the type of temperature rise condition has a greater influence on FGM thick plates.(2) Effect of FGM Combination and Width-Length Ratio. The nonlinear bending of two kinds of FG plates with different width-length ratios is calculated and compared in Table 7. In this example, the volume fraction index k = 2, b = 0.2 m, heating-up condition is C2, and two temperature rise conditions are considered.Table 7
The center deflection for Al/Zr2O3 and Si3N4/SUS304 FG plates under two kinds of temperature rise conditions and C2 with different width-length ratios h/b=1/20.
b/aUniform temperature risew (10−3m)Nonlinear temperature risew (10−3m)Relative deviation (%)Al/Zr2O310.42950.378111.972.50.75050.619017.5250.88620.645027.22Si3N4/SUS30410.20970.19168.622.50.34590.298413.7450.41620.311325.21Relative deviation is evaluated between the values in columns 3 and 4.Table7 shows that the relative deviation between two kinds of temperature rise conditions increases with width-length ratio b/a of the plate. It illustrates that the type of temperature distribution function has a greater influence on FGM rectangle plates. Besides, the relative deviation of Si3N4/SUS304 is smaller which presents that the effect of temperature rise condition is related to material combination. The temperature dependent FGM is less affected because the thermal material properties decrease with the increase of temperature, and the variation of nonlinear bending stiffness Ks is different under two kinds of temperature rise conditions. The following analyses only consider the nonlinear temperature rise condition.


## 4.3.2. Discussion of In-Plane Temperature Force

(1) Effect of Volume Fraction Index and Heating-Up Condition. To compare the nonlinear center deflection of FG plates under in-plane temperature forces, the uniformly distributed transverse load and elastic foundation are still not considered. The Al/Zr2O3 square plates are under different temperature increments, a=b=1m and h/a=1/20. And the center deflection-volume fraction index curves are shown in Figure 6.Figure 6
The center deflection-volume fraction index curves for Al/Zr2O3 FGM square plates under in-plane temperature force with different temperature increments: (a) heating the lower surface of the plate separately; (b) heating the upper surface of the plate separately; (c) heating the upper and lower surfaces of the plate simultaneously.
(a)(b)(c)As shown in Figure6(a), when the lower surface of the plate is heated separately, the center deflection is negative and changes more greatly with the volume fraction index k = 0∼2, and it is basically unchanged after k > 2. This is because the FGM properties change greatly with k = 0∼2, after that the FGM tends to be homogeneous with k increasing. In Figure 6(b), the upper surface of the plate is heated separately. The bending deformation of the FG plate under the in-plane temperature force is positive and similar to that of the plate subjected to transverse mechanical loads in normal temperature. The center deflection reaches the smallest when k = 0 and increases with the increase of k. This is because the thermal expansion coefficient α and thermal conductivity κ of ceramic are smaller than those of metal. Therefore, the temperature change of the upper surface has little influence on the deflection. Figure 6(c) shows that the homogeneous plate (k = 0) does not deform and the center deflection is 0 when the upper and lower surfaces of the plate are heated simultaneously. This is because the temperature stays the same along the thickness and the thermal moment Mθ=0. However, for FG plates, the temperature distribution function and thermal expansion coefficient are not only the functions of thickness z but also related to k, so the deformation is similar to Figure 6(a).It can be observed from Figure6 that the center deflection increases with the temperature. This is because the in-plane temperature force is pressure when the plate is heated up, which makes the geometric stiffness matrix in Kσ∗ negative and the total stiffness decrease. Thus, the in-plane temperature force can soften the structure. It has a greater effect on FG plates than on homogeneous plates when the bending deformation is small, and the effect is related to the material properties and environment temperature.(2) Coupling Effect with Transverse Loads. The uniformly distributed transverse load is added into the consideration on the basis of (1). The nondimensional center deflection-load curves for FGM square plates with different volume fraction indexes and heating-up conditions are shown in Figure 7.Figure 7
Comparisons of the nondimensional center deflection-load curves for Al/Zr2O3 FGM square plates with different volume fraction indexes and heating-up conditions: (a) C1; (b) C2; (c) C3.
(a)(b)(c)It can be observed that the nonlinear center deflection increases with the volume fraction indexk in Figure 7. This is because the larger k is, the more metal and the less ceramic the plate has. The composition changes from pure ceramic to ceramic-metal composite, which makes the elastic modulus and stiffness of the FG plate decrease. Besides, under different heating-up conditions, the center deflection values with the same volume fraction index are much different from one another when q∗ is small, and they tend to be close with the increase of q∗. For example, for the FG plate with k = 2 under three heating-up conditions, the nondimensional center deflections are −0.112, 0.002, and -0.004, respectively, and the maximum deviation is 0.114 while q∗ = 0. However, when q∗ = 50, the center deflections are 0.948, 1.013, and 0.980, respectively, and the maximum deviation becomes 0.065. It illustrates that the influence of the environment temperature and in-plane temperature force decreases with the increase of nonlinear bending deformation of the plate. The large deflection is mainly controlled by the transverse distributed loads.(3) Effect of Length-Width Ratio and Thickness-Length Ratio. The nondimensional center-load curves for Al/Zr2O3 FGM square plates with various length-width ratios are shown in Figure 8. The nonlinear bending of two kinds of FG plates with different length-width ratios and thickness-length ratios are calculated and compared in Table 8. In this example, the volume fraction index k = 2.Figure 8
The nondimensional center deflection-load curves for Al/Zr2O3 FG plates with C2 and various length-width ratios.Table 8
The nondimensional center deflection for two kinds of FG plates under C1andq∗ = 25 with different length-width ratios and thickness-length ratios (k1 = k2 = 0).
a/bh/aAl/Zr2O3SUS304/Si3N4Al/ZrO2w/h (S1)w/h (S2)Relative deviation (%)w/h (S1)w/h (S2)Relative deviation (%)w/h (S1)w/h (S2)Relative deviation (%)0.51/80−0.6564−0.2506161.93−0.8646−0.3371156.48−0.6424−0.2378170.141/40−0.1436−0.0703104.27−0.2491−0.150665.41−0.1716−0.0795115.851/20−0.0209−0.010991.740.04960.041918.38−0.0329−0.016797.0111/80−0.5611−0.2446129.39−0.02810.1177123.87−0.5760−0.2291151.421/400.52890.407429.820.73490.625517.490.46890.349034.361/200.60980.501321.640.57210.52199.620.71380.553828.891.51/800.65300.488233.760.90580.720225.770.65420.467539.941/401.12810.948019.001.09270.978011.731.03510.850621.691/201.00850.914310.300.99970.94206.130.95140.813716.9221/801.02400.832023.081.31841.197610.090.66090.513028.831/401.13821.018611.741.28921.22625.141.15041.042310.371/201.17091.08437.991.26411.24411.611.07340.99378.02Relative deviation is evaluated between the values in S1 and S2.As shown in Figure8, the nonlinear center deflection values increase with the increase of length-width ratio a/b of the plate. The larger the uniformly distributed transverse load q∗ is, the more significant the influence of the length-width ratio is. Table 8 shows that, with the increase of length-width ratio a/b and thickness-length ratio h/a, the relative deviations between situation S1 and S2 become smaller. It illustrates that the thermal environment effect decreases gradually, and the in-plane temperature force has a greater influence on the thin plate with the small length-width ratio. Meanwhile, the calculation method employed in this paper is suitable for the analysis of both thick and thin plates, and the locking-free thin/thick plate theory has a wide applicability.


## 4.3.3. Discussion of Elastic Foundation

(1) Effect of Transverse Distributed Load, Length-Width Ratio, and Heating-Up Condition. The nondimensional center deflection-load curves for Al/Zr2O3 FG plates are shown in Figure 9. The nonlinear bending of Al/ZrO2 plates under different length-width ratios and heating-up conditions are calculated and compared in Table 9. In this example, the plates are rested on Winkler–Pasternak elastic foundations with a=b=1m and volume fraction index k = 2.Figure 9
The nondimensional center deflection-load curves for elastically supported Al/Zr2O3 FG plates with (a) C1; (b) C2; (c) C3.
(a)(b)(c)Table 9
The center deflection for elastically supported Al/ZrO2 FG plates with different length-width ratios and heating-up conditions under q∗ = 25 and h/a=1/40.
a/bHeating-up conditions (Tm/Tc)k1=0k2=0w (10−3m)k1=10k2=10w (10−3m)k1=100k2=10w (10−3m)k1=100k2=100w (10−3m)Relative deviation I (%)Relative deviation II (%)Relative deviation III (%)0.5C110.740−1.297−0.524−0.097112.0859.6081.49C28.9881.0540.5320.09688.2749.5381.95C31.096−0.674−0.220−0.041161.5067.3681.361C111.7232.4161.8210.28579.3924.6384.35C220.0794.6693.2670.49076.7530.0385.00C310.1693.5502.4820.35165.0930.0885.861.5C122.1209.7307.9801.17056.0117.9985.34C223.92511.4029.3611.38352.3417.9085.23C322.04911.5299.2531.25047.7119.7486.492C128.76118.49414.5062.52735.7021.5682.58C229.19819.50417.3742.73733.2010.9284.25C328.90317.55615.5772.62439.2611.2783.15Relative deviation I is evaluated between the values in columns 3 and 4; relative deviation II is evaluated between the values in columns 4 and 5; Relative deviation III is evaluated between the values in columns 5 and 6.Figure9 shows that the stiffness of the FG plates increases and the center deflection decreases obviously after considering the effect of the elastic foundations, which indicates that the elastic foundations can stiffen the structure and weaken the bending behavior effectively. The center deflection decreases continuously with the increase of foundation parameter values. It can be observed that the shear stiffness K2 has a more obvious stiffening effect on the plates than the linear normal stiffness K1. As the uniformly distributed transverse load q∗ increases and the nonlinear deformation becomes larger and larger, the stiffening effect caused by K2 becomes more and more obvious.As shown in Table9, the relative deviations III are close to one another under different length-width ratios, which indicates that the influence of K2 is unrelated to the shape and size of the plate. However, the relative deviations I and II decrease with the increase of the length-width ratio. It illustrates that the shape and size can affect the stiffening effect caused by K1. The smaller the length-width ratio is, the larger the stiffening effect is. Besides, the relative deviations I under the three heating-up conditions become closer to one another as the length-width ratio increases. It shows that the influence of thermal environment decreases gradually, but the stiffening effect of the elastic foundation is unrelated to the heating-up condition.(2) Effect of Volume Fraction Index and Thickness-Length Ratio. The Al/Zr2O3 square plates a=b=1m are elastically supported. The center deflection with different volume fraction indexes and thickness-length ratios are shown in Table 10.Table 10
The center deflection for elastically supported Al/Zr2O3 FG plates with different volume fraction indexes and thickness-length ratios under q∗ = 25 and C2.
h/aVolume fraction indexk1=0k2=0w (10−3m)k1=10k2=10w (10−3m)k1=100k2=10w (10−3m)k1=100k2=100w (10−3m)Relative deviation I (%)Relative deviation II (%)Relative deviation III (%)1/80k = 0.56.2834.1602.6200.32933.7937.0287.44k = 28.7964.3962.7030.33150.0238.5187.75k = 39.0234.5052.7420.33250.0739.1387.891/50K = 0.517.5293.8262.6850.40078.1729.8285.10k = 219.4744.1682.8460.41378.6031.7285.49k = 319.9494.2102.8660.41478.9031.9285.551/25K = 0.525.8246.3024.5050.70875.6028.5184.28k = 226.6436.4344.5900.71775.8528.6684.38k = 327.2766.5504.6100.71975.9929.6284.401/10k = 0.593.08414.74810.5471.69884.1628.4983.90k = 295.22014.95310.6551.71084.3028.7483.95k = 399.58515.54011.0041.71684.4029.1984.41Relative deviation I is evaluated between the values in columns 3 and 4; Relative deviation II is evaluated between the values in columns 4 and 5; Relative deviation III is evaluated between the values in columns 5 and 6.Table10 shows that, for the thin plates (h/a = 1/80), the relative deviations I are smaller than that of the thick plates (h/a ≥ 1/50). The relative deviations II are larger than and the relative deviations III are close to those of the thick plates. It illustrates that the linear normal stiffness K1 has a greater influence on the thin plates, while the influence of shear stiffness K2 and the stiffening effect of the elastic foundation are unrelated to the thickness-length ratio. By comparing the three relative deviations with the same thickness-length ratio and different volume fraction indexes, the stiffening effect of the elastic foundation on FG plates is found to increase with k.


## 5. Conclusions

This paper proposes a SFPM to study the thick and thin FG plates in thermal environments based on a locking-free thin/thick plate theory. The material parameters and temperature are assumed to vary along the thickness direction of the plates. By comparing with the published models, the new model is validated to be correct with high precision and good convergence for rectangular plates. The numerical results show that the proposed method performs well in analyzing the nonlinear bending of elastically supported FG plates without shear locking. The effect of the temperature rise condition, in-plane temperature force and elastic foundation on the nonlinear bending under different parameters, including the volume fraction index, thickness-length ratio, length-width ratio, heating-up condition, and foundation stiffness, are discussed through several examples, and the conclusions are summarized as follows:(1)
For the nonlinear bending analysis of FG plates, the effect of thermal environment is related to material combination, and the temperature independent FGM is more affected. And thermal environment has greater influence on the plates with larger length-thickness ratios and with larger width-length ratios. Two kinds of temperature rise conditions are discussed, as the length-thickness ratios decrease and width-length ratios increase, the stiffnesses of the plates increase, so the relative deviations between two kinds of temperature rise conditions increase. And the uniform temperature rise condition has more effect on bending behavior than the nonlinear one with C2.(2)
The in-plane temperature force is generated when the cold shrinkage and thermal expansion deformation of the material is constrained by the boundary. It is related to the environment temperature and the material properties and can soften the structure and make the bending deformation increase with temperature rises because the in-plane temperature force produces bending moments that deform the plate. The softening effect is obvious on thin plates with small length-width ratio, and greater on FGMs than homogeneous materials. The phenomenon conforms to the general law of mechanics, the smaller the stiffness is, the larger the deformation is.(3)
As the nonlinear deformations of FG plates increase, the influence of the in-plane temperature and thermal environment decreases gradually. The large deflection is mainly controlled by external mechanical loads and obviously effected by the length-width ratio. Due to the small thickness of plate, the bending moment generated by the in-plane temperature force is limited, so the influence of the in-plane temperature and thermal environments is limited.(4)
The proposed method is applicable not only to theoretical research but also to some engineering application, such as the elastic foundation analysis. The analysis found that the elastic foundations can stiffen the structure and weaken the bending behavior effectively, and the nonlinear bending decreases with the increase of the foundation parameters. The linear normal stiffnessK1 has a greater influence on thin plates. The shear stiffness K2 has a more obvious stiffening effect than K1, and the larger the mechanical load is, the more significant the effect is. The stiffening effect is unrelated to the thickness-length ratio and the heating-up condition, but it can be affected by the length-width ratio a/b and volume fraction index k. The smaller the value of a/b and the larger the k, the stronger the effect.



---

*Source: 2943705-2020-07-10.xml*