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

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. /e displacements of the plate are constructed by the linear combination of orthogonal functions and cubic B-spline interpolation functions with high efficiency for modeling. /e 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./e 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 newmodel are validated for rectangular plates. Moreover, uniform and nonlinear temperature rise conditions are discussed, respectively./e effect of the temperature distribution, in-plane temperature force, and elastic foundation on nonlinear bending under different parameters are discussed in detail.


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. erefore, 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. e 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. e 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]. erefore, some researchers devoted themselves to establish a general element and method to calculate both thick and thin plates [14][15][16], and a locking-free thin/thick plate theory was proposed. In this paper, this theory is based on the FSDT but takes the inplane 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][18][19]. e 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. e equilibrium equation is established with the minimum potential energy principle. e 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. e 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).
e 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][30][31][32]. All these works proved that solving mechanical problems by simple linearization cannot meet the precision requirement in important structural calculations. e 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. us, 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. e 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.

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 where x k is the spline node, with e cubic B-spline function can be obtained with n � 3, i.e., Suppose the length of a FG plate is a and the width is b, and its spline discretization model is shown in Figure 1.
For rectangular plates, there are no restrictions on the selection of x-axis and y-axis by SFPM. As shown in Figure 1, the plate structure can be divided by the interval value a x in x direction, and the local coordinates of corresponding nodes are given by where N stands for the number of subdivisions in x direction, and e cubic B-spline function φ 3 (x) is employed in this paper with the spline node x k � 0 as its center, which can be expressed as e graphs of cubic B-spline function and its derivatives can by expressed as follows. Figure 2 shows that the cubic B-spline function and its derivatives are piecewise, symmetrical, and compact.
us, the displacement function in x direction of the plate can be represented by a linear combination of spline interpolation functions ϕ i (x) and generalized parameters c i . e 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 where Z m (y) 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 Z m � sin((mπ/b)y). e spline interpolation functions [ϕ] can be constructed in the matrix form as follows: where [I] is a (N − 3) × (N − 3) identity matrix, and It is very simple to deal with the displacement boundary conditions using spline interpolation functions with the following properties:

Mathematical Problems in Engineering
erefore, 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 to i � − 1 and i � N + 1 need to be substituted for the two ends' simply supported plate in x direction.

Physical Properties of FGM.
e properties of two materials are assumed to vary along the thickness of FG plates and follow the power function distribution; then, Young's modulus E(z), 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]) 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. e temperature dependence of the FGM should be considered in a thermal environment. Since Young's modulus E 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]) where T is the corresponding temperature and P − 1 , P 0 , P 1 , P 2 , and P 3 are the temperature-dependent coefficients.

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). ere are two temperature rise conditions to be considered in this paper.

Uniform Temperature Rise Condition.
In this situation, the temperature is uniformly raised from the lower surface to upper surface. e temperature variation ΔT � T(z) − T 0 is a constant, in which T(z) is the temperature distribution function along the thickness direction (Javaheri and Eslami [39]) and T 0 � 300 K is the reference temperature without considering the thermal strain.

Nonlinear Temperature Rise Condition.
e steadystate heat transfer equation of FG plates can by expressed as follows (Lal et al. [20]): where κ(z) is given in detail in (13) and T(z) is the temperature distribution function along the thickness direction. Both of them are the functions of z. e upper and lower surface temperature of the plate are assumed to be T c and T m , respectively; then, the temperature distribution function inside FG plates can be determined as (Lal et al. [20]) where Mathematical Problems in Engineering 5 in which κ cm � κ c − κ m , κ c and κ m are the material thermal conductivities of the upper and lower surfaces, respectively.

Nonlinear Geometric Equation.
Using the locking-free thin/thick plate theory for thick and thin plates, which takes the in-plane displacements u and v, deflection w, and shear strains c xz and c yz as the independent unknowns, the shear strain conditions of thick plate (c ≠ 0) and thin plate (c � 0) can be met simultaneously. e specific formula derivation and description is given in Appendix B. us, the displacement field can be obtained by where u, v, and w are the total displacements, u 0 , v 0 , and w 0 are the corresponding displacements of a point on the neutral surface, and φ x and φ y are the rotations about the yaxis and x-axis, respectively. In order to establish the stiffness equation easily, the strain equation of the plate is expressed as in which the linear strains are and the nonlinear strain is

Constitutive Equations in ermal Environment.
In a thermal environment, the constitutive equations of FG plates on Winkler-Pasternak elastic foundation can be defined by (Lal et al. [20]) where [Q f ] is the elastic modulus of FG plates considering temperature dependence, α { } is the thermal expansion coefficient array, Pf is the foundation reaction per unit area, K 1 and K 2 are the linear normal stiffness and shear stiffness of the elastic foundation, respectively, ∇ 2 � (z 2 /zx 2 ) + (z 2 /zy 2 ) is the second order Laplace differential operator, S T is the entropy density function, and p is the specific heat capacity of materials. e matrix and vector expressions for the aforementioned physical quantities are defined as follows: e following stress expression is employed to analyze the geometrically nonlinear problems in the thermal environment: 6 Mathematical Problems in Engineering For the plate analysis, stresses are usually rewritten as the generalized forces as follows: in which where N θ and M θ are the thermal force and thermal moment caused by temperature, respectively.
, and [C] are the stretching stiffness, stretchingbending coupling stiffness, bending stiffness, and shearing stiffness, respectively, and the specific expressions are given by erefore, (26) can be simplified as in which

Spline Discretization Modeling.
e FG plate is discretized uniformly in x direction, and the displacement functions can be constructed by the cubic B-spline interpolation functions and orthogonal functions as follows: where r is the convergent series, [ϕ] is the cubic B-spline interpolation function defined by (9), and X m (y), Y m (y), Z m (y), F m (y), and G m (y) 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 In order to satisfy the boundary conditions at y � 0 and b, the corresponding orthogonal functions are suggested to be Mathematical Problems in Engineering us, (32) can be expressed in the matrix form as follows: 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., So, the linear strain in (31) can be denoted by the generalized displacement parameter as follows: e nonlinear strain can be constructed as in which Substitute (37) into (44), i.e., in which us, (31) can be rewritten as

Nonlinear Equilibrium Equation.
Based on the virtual work principle, the virtual work equation of a thermomechanical FG plate is given by  (27). e nonlinear equilibrium equation of the FG plate on Winkler-Pasternak elastic foundation is given by where K S stands for the nonlinear bending stiffness, K F denotes the elastic foundation stiffness, K θ is the thermal geometric stiffness, and F R is the thermal load caused by the environment temperature, which can be expressed as follows, respectively: 8 Mathematical Problems in Engineering e differential of Ψ can be written as in which where [N f ] � N x N xy N xy N y 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. e nonlinear equilibrium equation (50) can be solved by the Newton-Raphson method.

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. e first kind is the combination of Al/ZrO 2 and Al/Zr 2 O 3 , which is temperature independent.
e 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]. e second kind is SUS304/Si 3 N 4 , which is temperature dependent. e temperature-dependent coefficients are listed in Table 2 from Fu et al. [41]. e thermal material properties of SUS304 and Si 3 N 4 under different temperatures are calculated in Table 3 according to (15) in Section 3.  Table 4. Numerical results of CLPT and FSDT are also included. In this example, the spline discretization level N × r � 12 × 3. e FG plates are subjected to transverse sinusoidal load q � q 0 sin((π/a)x)sin((π/b)y), and q 0 is the intensity of the transverse load at the plate center. e nondimensional transverse displacement w � (10E c h 3 /q 0 a 4 )w, where E c � 380 GPa, and w is the center deflection.

Model Comparisons and
It can be noticed that, in Table 4, 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.

Convergence studies of SFPM.
e nonlinear bending of Al/ZrO 2 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]. e numerical results are listed in Table 5 with various volume fraction indexes. e 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 � 100wE m b 3 /12(1 − μ 2 )qa 4 , (a/h) � 5, (a/b) � 1, μ � 0.3, and E m � 70 GPa.
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.  Table 2: Temperature-dependent coefficients for the material properties of the second kind of FGM (from Fu et al. [41]).

Material
Properties   [20]. e FG plates subjected to uniformly distributed transverse loads q with a nonlinear temperature rise, the nondimensional pressure q * � (qb 4 /E m h 4 ), ΔT � T(z) − T 0 , and T 0 � 300 K. e values of E m with different temperature increments are listed in Table 3. e nondimensional linear normal and shear stiffness of the foundation are k 1 � (K 1 a 4 /E m h 3 ) and k 2 � (K 2 a 2 /E m h 3 ), respectively. e volume fraction index k � 2, (b/a) � 1, and (b/h) � 20. e model is shown in Figure 4 and the comparisons can be observed in Figure 5.
As shown in Figure 5, the results of this paper are in good agreement with the published results using different approaches.
e 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. is 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. is is because the thermal load F R increases with the increase of ΔT, which causes the nonlinear bending deformation to increase.

Parameter Studies.
In this section, the following dimensionless parameters, nondimensional uniformly distributed transverse load (q * ), and foundation stiffness (k 1 and k 2 ) are defined as (59) ere are three heating-up conditions to be considered (unless otherwise stated): C1: T m � 400 K, T c � 300 K; C2: T m � 300 K, T c � 400 K; and C3: T m � 400 K, T c � 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.

Discussion of Temperature Distribution (1) Effect of Volume Fraction Index and Length-ickness
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/ZrO 2 square plates (a � b � 1 m) under two kinds of temperature rise conditions with different length-thickness ratios and volume fraction indexes are shown in Table 6. Table 6 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.
is 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 K s 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. e 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 shows that the relative deviation between two kinds of temperature rise conditions increases with widthlength 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. e temperature dependent FGM is less affected because the thermal material properties decrease with the increase of temperature, and the variation of nonlinear bending stiffness K s is different under two kinds of temperature rise conditions. e following analyses only consider the nonlinear temperature rise condition.

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. e Al/Zr 2 O 3 square plates are under different temperature increments, a � b � 1 m and (h/a) � (1/20) . And the center deflection-volume fraction index curves are shown in Figure 6.
As shown in Figure 6(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.
is 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. e 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. e center deflection reaches the smallest when k � 0 and increases with the increase of k. is is because the thermal expansion coefficient α and thermal conductivity κ of ceramic are smaller than those of metal. erefore, 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. is 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 Figure 6 that the center deflection increases with the temperature. is 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. us, the inplane 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. e uniformly distributed transverse load is added into the consideration on the basis of (1). e nondimensional center deflectionload curves for FGM square plates with different volume  fraction indexes and heating-up conditions are shown in Figure 7. It can be observed that the nonlinear center deflection increases with the volume fraction index k in Figure 7. is is because the larger k is, the more metal and the less ceramic the plate has. e 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. e large deflection is mainly controlled by the transverse distributed loads.

(3) Effect of Length-Width Ratio and ickness-Length Ratio.
e nondimensional center-load curves for Al/Zr 2 O 3 FGM square plates with various length-width ratios are shown in Figure 8. e 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. As shown in Figure 8, the nonlinear center deflection values increase with the increase of length-width ratio a/b of the plate. e larger the uniformly distributed transverse load q * is, the more significant the influence of the lengthwidth ratio is. Table 8 shows that, with the increase of lengthwidth 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.

Discussion of Elastic Foundation. (1) Effect of Transverse Distributed Load, Length-Width Ratio, and Heating-Up
Condition.
e nondimensional center deflection-load curves for Al/Zr 2 O 3 FG plates are shown in Figure 9. e nonlinear bending of Al/ZrO 2 plates under different lengthwidth 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 � 1 m and volume fraction index k � 2.  Figure 9 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. e center deflection decreases continuously with the increase of foundation parameter values. It can be observed that the shear stiffness K 2 has a more obvious stiffening effect on the plates than the linear normal stiffness K 1 . As the uniformly distributed transverse load q * increases and the nonlinear deformation becomes larger and larger, the stiffening effect caused by K 2 becomes more and more obvious.
As shown in Table 9, the relative deviations III are close to one another under different length-width ratios, which indicates that the influence of K 2 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 K 1 . e 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 ickness-Length Ratio. e Al/Zr 2 O 3 square plates (a � b � 1 m) are elastically supported. e center deflection with different volume fraction indexes and thickness-length ratios are shown in Table 10.  Relative deviation is evaluated between the values in S1 and S2. 16 Mathematical Problems in Engineering

a/b
Heating-up conditions (T m /T c )  Table 10 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). e 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 K 1 has a greater influence on the thin plates, while the influence of shear stiffness K 2 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.

Conclusions
is paper proposes a SFPM to study the thick and thin FG plates in thermal environments based on a locking-free thin/ thick plate theory. e 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. e numerical results show that the proposed method performs well in analyzing the nonlinear bending of elastically supported FG plates without shear locking. e 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 lengththickness 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) e 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. e softening effect is obvious on thin plates with small length-width ratio, and greater on FGMs than homogeneous materials. e 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. e large deflection is mainly controlled by external mechanical loads and obviously effected by the lengthwidth 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 inplane temperature and thermal environments is limited.
(4) e proposed method is applicable not only to theoretical research but also to some engineering application, such as the elastic foundation analysis. e 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. e linear normal stiffness K 1 has a greater influence on thin plates. e shear stiffness K 2 has a more obvious stiffening effect than K 1 , and the larger the mechanical load is, the more significant the effect is. e 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. e smaller the value of a/b and the larger the k, the stronger the effect.

(B.7)
For the thin plates, h ⟶ 0, λ ⟶ 0, i.e., V c � 0. So, equation (B.5) degenerate to the bending shear strain energy of thin plates, i.e., Equations (B.9) and (B.2) show that there is no false shear strain for thin plates and the shear deformation is considered for thick plates. e locking-free thin/thick plate theory is suitable for the analysis of different thickness plates without shear locking.

Data Availability
e results in the paper are all implemented by MATLAB, and the data used to support the findings of this study are cited within the article at references and are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.