#### 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.

#### 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 *n*th-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 bywhere *x _{k}* is the spline node, with

The 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 bywhere *N* stands for the number of subdivisions in *x* direction, and

The cubic B-spline function is employed in this paper with the spline node as its center, which can be expressed as

The 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.

**(a)**

**(b)**

**(c)**

**(d)**

Thus, the displacement function in *x* direction of the plate can be represented by a linear combination of spline interpolation functions and generalized parameters . 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 byin whichwhere is the orthogonal function. For the two ends’ simply-supported plate in *y* direction, when and *b*, so the corresponding orthogonal function is suggested to be . The spline interpolation functions can be constructed in the matrix form as follows:where is a identity matrix, and

It is very simple to deal with the displacement boundary conditions using spline interpolation functions with the following properties:

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 to and 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 modulus , material density , thermal expansion coefficient , and thermal conductivity 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.

The 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.

##### 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 is a constant, in which 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.

###### 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]):where is given in detail in (13) and 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 be *T*_{c} and *T*_{m}, respectively; then, the temperature distribution function inside FG plates can be determined as (Lal et al. [20])wherein which , and 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 displacements *u* and , deflection , and shear strains and as the independent unknowns, the shear strain conditions of thick plate () and thin plate () can be met simultaneously. The specific formula derivation and description is given in Appendix B. Thus, the displacement field can be obtained bywhere *u*, , and are the total displacements, , , and are the corresponding displacements of a point on the neutral surface, and and 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 asin which the linear strains areand the nonlinear strain is

##### 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])where 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, is the second order Laplace differential operator, 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:

The following stress expression is employed to analyze the geometrically nonlinear problems in the thermal environment:

For the plate analysis, stresses are usually rewritten as the generalized forces as follows:in whichwhere and 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

Therefore, (26) can be simplified asin which

##### 3.5. Spline Discretization Modeling

The 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 , , , , and 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

Thus, (32) can be expressed in the matrix form as follows:in which 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.,in which , and are the rotations of normal to the neutral surface about the *y*-axis and *x*-axis respectively, is the generalized displacement parameter, and both and are the displacement shape functions of the plate, i.e.,

So, the linear strain in (31) can be denoted by the generalized displacement parameter as follows:where , which is similar to and , and

The nonlinear strain can be constructed asin which

Substitute (37) into (44), i.e.,in which

Thus, (31) can be rewritten as

##### 3.6. Nonlinear Equilibrium Equation

Based on the virtual work principle, the virtual work equation of a thermomechanical FG plate is given bywhere is the vector sum of internal and external forces, {*F*} denotes the transverse external force, and are the virtual displacement and virtual strain, respectively, and 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 bywhere stands for the nonlinear bending stiffness, denotes the elastic foundation stiffness, is the thermal geometric stiffness, and is the thermal load caused by the environment temperature, which can be expressed as follows, respectively:

The differential of can be written asin whichwhere 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/ZrO_{2} and Al/Zr_{2}O_{3}, 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/Si_{3}N_{4}, 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 Si_{3}N_{4} under different temperatures are calculated in Table 3 according to (15) in Section 3.

##### 4.2. Model Comparisons and Validations

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

The bending of Al/Al_{2}O_{3} 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 , and *q*_{0} is the intensity of the transverse load at the plate center. The nondimensional transverse displacement , where *E*_{c} = 380 GPa, and is the center deflection.

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.

###### 4.2.2. Convergence studies of SFPM

The 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]. 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 , , , , and *E*_{m} = 70 GPa.

**(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/Si_{3}N_{4} 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 , , and . The values of with different temperature increments are listed in Table 3. The nondimensional linear normal and shear stiffness of the foundation are and , respectively. The volume fraction index *k* = 2, , and . The 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. 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 is, the larger the nonlinear center deflection is. This is because the thermal load increases with the increase of , which causes the nonlinear bending deformation to increase.

##### 4.3. Parameter Studies

In this section, the following dimensionless parameters, nondimensional uniformly distributed transverse load (), and foundation stiffness (*k*_{1} and *k*_{2}) are defined as

There 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.

###### 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/ZrO_{2} square plates 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. This is because the overall temperature of the plates is higher with the uniform temperature rise, so the influence of the thermal geometric stiffness is greater than the nonlinear bending stiffness 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 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 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/Zr_{2}O_{3} square plates are under different temperature increments, and _{.} And the center deflection-volume fraction index curves are shown in Figure 6.

**(a)**

**(b)**

**(c)**

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. 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 . 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. This is because the in-plane temperature force is pressure when the plate is heated up, which makes the geometric stiffness matrix in 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.

**(a)**

**(b)**

**(c)**

It can be observed that the nonlinear center deflection increases with the volume fraction index *k* 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 is small, and they tend to be close with the increase of *.* 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 = 0. However, when = 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/Zr_{2}O_{3} 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.

As shown in Figure 8, 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 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/Zr_{2}O_{3} FG plates are shown in Figure 9. The nonlinear bending of Al/ZrO_{2} 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 and volume fraction index *k* = 2.

**(a)**

**(b)**

**(c)**

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. The 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 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}. 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/Zr_{2}O_{3} square plates are elastically supported. The center deflection with different volume fraction indexes and thickness-length ratios are shown in Table 10.

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). 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 *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*.

#### 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 stiffness *K*_{1} has a greater influence on thin plates. The 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. 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.

#### Appendix

#### A. Calculation of the Nonlinear Spline Integral

The integral of product of three B-spline interpolation functions is needed when using SFPM for nonlinear calculation. Taking the calculation of as an example,

According to the matrix form of spline interpolation functions equations (32) and (33), equation (A.1) can be rewritten asgiven

Substituting equation (A.3) into equation (A.2) givesin which

Forgiven , and , so equation (A.6) can be rewritten aswhere .

#### B. The Specific Formula Derivation and Description of Locking-Free Thin/Thick Plate Theory

Based on the locking-free thin/thick plate theory for thick and thin plates, the displacement field can be expressed as follows:

The bending shear strain energy of thick plates is given byin whichgiven

Substituting equation (B.4) into equation (B.2) gives

For , i.e.,in which

For the thin plates, , , i.e., . So, equation (B.5) degenerate to the bending shear strain energy of thin plates, i.e.,

Equation (B.2) degenerates to

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. The locking-free thin/thick plate theory is suitable for the analysis of different thickness plates without shear locking.

#### Data Availability

The 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

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This research was supported by the grants awarded by the National Natural Science Foundation of China (no. 11562001), Science and Technology Major Project of Guangxi Province (no. AA18118029), and Innovation Project of Guangxi Graduate Education (no. YCSW2018047).