Thermal Buckling and Free Vibration Analysis of Functionally Graded Plate Resting on an Elastic Foundation According to High Order Shear Deformation Theory Based on New Shape Function

Abstract

Functionally graded square and rectangular plates of different thicknesses placed on the elastic foundation modeled according to the Winkler-Pasternak theory have been studied. The thermal and mechanical characteristics, apart from Poisson’s ratio, are considered to continuously differ through the thickness of the studied material as stated in a power-law distribution. A mathematical model of functionally graded plate which include interaction with elastic foundation is defined. The equilibrium and stability equations are derived using high order shear deformation theory that comprises various kinds of shape function and the von Karman nonlinearity. A new analytically integrable shape function has been introduced. Hamilton’s principle has been applied with the purpose of acquiring the equations of motion. An analytical method for identifying both natural frequencies and critical buckling temperature for cases of linear and nonlinear temperature change through the plate thickness has been established. In order to verify the derived theoretical results on numerical examples, an original program code has been implemented within software MATLAB. Critical buckling temperature and natural frequencies findings are shown below. Previous scientific research and papers confirms that presented both the theoretical formulation and the numerical results are accurate. The comparison has been made between newly established findings based on introduced shape function and the old findings that include 13 different shape functions available in previously published articles. The final part of the research provides analysis and conclusions related to the impact of the power-law index, foundation stiffness, and temperature gradient on critical buckling temperature and natural frequencies of the functionally graded plates.

1. Introduction

Due to a variety of organic and inorganic compounds, progress and growth has been made possible when it comes to present-day materials, advanced polymers, engineering alloys, structural ceramics, etc. [1]. The use of new materials happens as a consequence of current technology trends. Material properties are expected to adapt to current changes in technology and to have a spectrum of functions and characteristics that have not been introduced yet. As a result, materials are merged, and their advantages are preserved. Functionally graded materials (FGM) meet the needs of all the mentioned requests in technology.

Belonging to the family of engineering composite materials, FGM are modern materials that feature a continuous or discontinuous variation of the chemical composition through a defined direction. Detailed analysis and scientific experimentation have shown that FGM are able to constitute a gradient property, which is not the case with other homogeneous materials or composites. Present-day engineering faces a significant number of obstacles that could be overcome with these newly established materials with functionally graded composition. Mechanical characteristics like Poisson’s ratio, density of material, modulus of elasticity, shear modulus, and thermal expansion coefficient change through a defined direction, where a property gradient can be stepwise or continuous (linear, exponential, or parabolic) (Figure 1).

The existing materials and their qualities could be completely utilized by the FGM. The following factors are encompassed: reduction of transverse shear stress, the enhancement of mechanical and thermal characteristics as well as delamination prevention between the layers, which is one of the most crucial and the most commonly studied issues related to composite laminates [2] (Figure 2).

Metal/ceramics are the most frequently used FGM, metal being superior to ceramics regarding strength, toughness, and high thermal conductivity, while ceramics features a good temperature resistance, low thermal conductivity, and good antioxidant characteristics. FGM that contain both metal and ceramic constituents enhance thermal-mechanical characteristics between layers. As a consequence of continuous change of properties at the interface, FGM avoids delamination.

Functionally graded materials, being modern materials in the group composite materials, represent a popular topic discussed among numerous authors in recent years, as evidenced by a large number of publications in renowned journals in the field of composite materials. The actuality and the importance of the topic are addressed in numerous reviewed papers [3,4,5,6], which undoubtedly indicates the intention of the authors to illustrate the current state of research in this area and point to further research directions related to this very interesting area. Overall, FGM plates and shells under the impact of mechanical load or temperature can be studied using a 3D elastic theory or equivalent layer theories, which means classical plate theory (CPT), first-order shear deformation theory (FSDT), and higher-order shear deformation theory (HSDT). In order to eliminate the disadvantages of CPT in the analysis of moderately thick and thin plates [7], as well as to exclude the shear correction factors in FSDT [8], higher-order shear deformation theories (HSDT) were introduced. The most commonly used HSDT theory is third-order shear deformation theories (TSDT) developed by Reddy [9,10] for laminate composite materials, taking into account the effects of shear deformation and satisfying condition that the laminates upper and lower surface stress-free. Later, the aforementioned theory was applied to the analysis of the FGM plates [11]. Subsequently, a number of authors have used Reddy’s TSDT theory in the analysis of free vibrations and the dynamic stability [12,13] of FGM plates, with or without the interaction of the plate and the elastic foundation [14]. The impact of temperature, plate geometry, and material on free vibrations was studied in [15]. In addition to TSDT, a special HSDT group of theories, which has been developed in order to exclude the need to use correction factors, includes HSDT theories with shape functions. In general, there are different types of shape functions.

The initial idea of developing FGM was aimed at obtaining a material with high resistance to temperature gradient on one side and also good mechanical properties on the other side. For this reason, a number of authors have addressed the behavior problems of FG plates made of metal-ceramic constituents under mechanical and thermal, static, and dynamic loads, applying the theories mentioned above. The equilibrium and stability equations of thin, moderately thick, and thick FGM plates exposed to the impact of temperature have been considered in the area of linear [16,17,18] and nonlinear elasticity [19,20]. The effect of uniform, linear, and nonlinear temperature change in the direction of plate thickness has been analyzed in dynamic problems [21,22,23]. The plate/foundation interaction and the impact of the elasticity of the foundation, modeled by the Winkler-Pasternak model, were analyzed by the authors in [24,25,26,27]. The problem of constrained multi objective optimization performed for mass and material cost minimization as well as the minimization of stress failure criteria or maximization of natural frequency is studied in [28,29]. A recent trend of research in the area of FGM is quantifying uncertainty [30,31,32,33].

Ultimately, the final aim of all the previously mentioned research and studies is the application of FGM in various fields of engineering and industry. Although initially used as FGM material for thermal coating in spacecraft, due to their advantages over conventional materials, today FGM is increasingly being used in medicine [34], dentistry [35], the energy and nuclear sectors [36], the automotive industry [37], the military industry [38], optoelectronics [39], and others.

2. Mathematical Model of the Functionally Graded Plate Placed on Elastic Foundation

This paper deals with the FGM rectangular plate (a-length, b-width, and h-height) resting on an elastic foundation, where the z-axis has a direction of thickness h and the x-y plane represents the mid surface of the plate (Figure 3). Mathematical model of elastic foundation is defined by the use of Winkler-Pasternak type of two parameters elastic foundation:

• k₀ is stiffness of Winkler foundation,
• k₁ represents shear stiffness (Pasternak coefficient).

Power law distribution is used to define Young’s modulus of elasticity, thermal expansion coefficient, and temperature change through plate thickness [40]:

$$\begin{array}{c}\mathrm{E}\left(\mathrm{z}\right)={\mathrm{E}}_{\mathrm{m}}+{\mathrm{E}}_{\mathrm{cm}}{\left(\frac{1}{2}+\frac{\mathrm{z}}{\mathrm{h}}\right)}^{\mathrm{p}},\hspace{1em}\hspace{1em}\hspace{1em}{\mathrm{E}}_{\mathrm{cm}}={\mathrm{E}}_{\mathrm{c}}-{\mathrm{E}}_{\mathrm{m}},\\ \mathsf{\alpha }\left(\mathrm{z}\right)={\mathsf{\alpha }}_{\mathrm{m}}+{\mathsf{\alpha }}_{\mathrm{cm}}{\left(\frac{1}{2}+\frac{\mathrm{z}}{\mathrm{h}}\right)}^{\mathrm{p}},\hspace{1em}\hspace{1em}\hspace{1em}{\mathsf{\alpha }}_{\mathrm{cm}}={\mathsf{\alpha }}_{\mathrm{c}}-{\mathsf{\alpha }}_{\mathrm{m}},\\ \mathrm{T}\left(\mathrm{z}\right)={\mathrm{T}}_{\mathrm{m}}+\Delta {\mathrm{T}}_{\mathrm{cr}}{\left(\frac{1}{2}+\frac{\mathrm{z}}{\mathrm{h}}\right)}^{\mathrm{s}},\hspace{1em}\hspace{1em}\hspace{1em}\Delta {\mathrm{T}}_{\mathrm{cr}}={\mathrm{T}}_{\mathrm{c}}-{\mathrm{T}}_{\mathrm{m}}.\end{array}$$

(1)

The analytical procedure for determining natural frequencies as well as critical buckling temperature for both linear and nonlinear temperature change across the FG plate thickness is hereby developed.

3. Equilibrium and Stability Equations of FG Plate Placed on Elastic Foundation

By introducing HSDT with shape functions, numerous authors have eliminated the disadvantages of CPT and FSDT.

In order to produce good results for specific dynamical and static problems, many of these shape functions have been introduced. It is important to note that the shape functions proposed by various authors (

Table 1. Different type of shape functions.

Number of Function	Shape Function, f(z)
No. 1 [41]	$\left(\mathrm{z}/2\right)\left({\mathrm{h}}^2/4-{\mathrm{z}}^2/3\right)$
No. 2 [42]	$\left(5\mathrm{z}/4\right)\left(1-4{\mathrm{z}}^2/3{\mathrm{h}}^2\right)$
No. 3 [43]	$\left(\mathrm{h}/\mathsf{\pi }\right)\sin \left(\mathsf{\pi }\mathrm{z}/\mathrm{h}\right)$
No. 4 [44]	$\sin \left(\mathsf{\pi }\mathrm{z}/\mathrm{h}\right){\mathrm{e}}^{\frac{1}{2}\cos \left(\frac{\mathsf{\pi }}{\mathrm{h}}\mathrm{z}\right)}+\left(\mathsf{\pi }\mathrm{z}/2\mathrm{h}\right)$
No. 5–6 [45]	$\tan \left(\mathrm{mz}\right)-\mathrm{zm}{\sec }^2\left(\mathrm{mh}/2\right),\hspace{0.17em}\hspace{0.17em}\mathrm{m}=\left\{1/5\mathrm{h},\mathsf{\pi }/2\mathrm{h}\right\}$
No. 7 [46,47]	$\mathrm{z}\exp \left(-2{\left(\mathrm{z}/\mathrm{h}\right)}^2\right),\mathrm{z}\exp \left(\frac{-2{\left(\mathrm{z}/\mathrm{h}\right)}^2}{\ln \mathsf{\alpha }}\right), \forall \mathsf{\alpha }>0$
No. 8 [48]	$\mathrm{z}\cdot {2.85}^{-2{\left(\mathrm{z}/\mathrm{h}\right)}^2}+0.028\mathrm{z}$
No. 9 [49]	$\mathsf{\xi }\left[\left(\mathrm{h}/\mathsf{\pi }\right)\sin \left(\mathsf{\pi }\mathrm{z}/\mathrm{h}\right)-\mathrm{z}\right],\mathsf{\xi }=\left\{1,1/\cosh (\mathsf{\pi }/2)-1\right\}$
No. 10 [50]	$\mathrm{h}\sinh \left(\mathrm{z}/\mathrm{h}\right)-\mathrm{z}\cosh \left(1/2\right)$
No. 11 [51]	$\mathrm{z}\sec \mathrm{h}\left({\mathrm{z}}^2/{\mathrm{h}}^2\right)-\mathrm{z}\sec \mathrm{h}\left(\mathsf{\pi }/4\right)\left[1-\left(\mathsf{\pi }/2\right)\tanh \left(\mathsf{\pi }/4\right)\right]$
No. 12 [51]	$\left(3\mathsf{\pi }/2\right)\mathrm{h}\tanh \left(\mathrm{z}/\mathrm{h}\right)-\left(3\mathsf{\pi }/2\right)\mathrm{z}\sec {\mathrm{h}}^2\left(1/2\right)$
No. 13 [52]	$\frac{\mathrm{z}\cos \left(1/2\right)}{-1+\cos \left(1/2\right)}-\frac{\mathrm{h}\sin \left(\mathrm{z}/\mathrm{h}\right)}{-1+\cos \left(1/2\right)}$

) are not generally applicable to all types of problems.

This paper introduces a new shape function:

$$\mathrm{f}\left(\mathrm{z}\right)=\mathrm{z}\left(\cosh \left(\frac{\mathrm{z}}{\mathrm{h}}\right)-1.388\right)$$

(2)

The starting point for developing a new shape function was a comparative analysis of the advantages and disadvantages of existing shape functions given in Table 1. The initial conditions which had to be satisfied in developing a new shape function are: the function has to be an odd function of the thickness coordinate, the function has to satisfy zero stress conditions for out of plane shear stresses, the function has to be analytically integrable in order to additionally increases the precision of the results obtained. Table 1 shows that the newly developed shape function belongs to the category of simple mathematic functions. This makes the process of integration easier and it consequently reduces the calculation time significantly. Since the function is analytically integrable, it is not necessary to switch to numeric integration, and that fact additionally increases the precision of the results obtained.

Here, the assumed form of the displacement field [50] is:

$$\begin{array}{l}\mathrm{u}\left(\mathrm{x},\mathrm{y},\mathrm{z},\mathrm{t}\right)={\mathrm{u}}_0\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)-\mathrm{z}\frac{\partial \mathrm{w}}{\partial \mathrm{x}}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)+\mathrm{f}\left(\mathrm{z}\right){\mathsf{\theta }}_{\mathrm{x}},\\ \mathrm{v}\left(\mathrm{x},\mathrm{y},\mathrm{z},\mathrm{t}\right)={\mathrm{v}}_0\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)-\mathrm{z}\frac{\partial \mathrm{w}}{\partial \mathrm{y}}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)+\mathrm{f}\left(\mathrm{z}\right){\mathsf{\theta }}_{\mathrm{y}},\\ \mathrm{w}\left(\mathrm{x},\mathrm{y},\mathrm{z},\mathrm{t}\right)={\mathrm{w}}_0\left(\mathrm{x},\mathrm{y},\mathrm{t}\right).\end{array}$$

(3)

It is necessary to apply the relations between strains and displacements based on von Karman’s non-linear theory of elasticity so as to define the components of unit loads [53]. Considering the effect of the temperature change (1) and thermal expansion that cause a strain $\mathsf{\alpha }\Delta \mathrm{T}$, as well as using the generalized Hooke’s law, the following unit load components are obtained:

$$\begin{array}{l}\left\{\begin{array}{l}{\mathrm{N}}_{\mathrm{xx}}\\ {\mathrm{N}}_{\mathrm{yy}}\\ {\mathrm{N}}_{\mathrm{xy}}\end{array}\right\}={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\left\{\begin{array}{l}{\mathsf{\sigma }}_{\mathrm{xx}}\\ {\mathsf{\sigma }}_{\mathrm{yy}}\\ {\mathsf{\tau }}_{\mathrm{xy}}\end{array}\right\}\mathrm{dz}}={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\left[\begin{array}{ccc}{\mathrm{Q}}_{11}(\mathrm{z})& {\mathrm{Q}}_{12}(\mathrm{z})& 0\\ {\mathrm{Q}}_{12}(\mathrm{z})& {\mathrm{Q}}_{22}(\mathrm{z})& 0\\ 0& 0& {\mathrm{Q}}_{66}(\mathrm{z})\end{array}\right]\left\{\begin{array}{l}\frac{\partial {\mathrm{u}}_0}{\partial \mathrm{x}}+\frac{1}{2}{\left(\frac{\partial \mathrm{w}}{\partial \mathrm{x}}\right)}^2\\ \frac{\partial {\mathrm{v}}_0}{\partial \mathrm{y}}+\frac{1}{2}{\left(\frac{\partial \mathrm{w}}{\partial \mathrm{y}}\right)}^2\\ \frac{\partial {\mathrm{u}}_0}{\partial \mathrm{y}}+\frac{\partial {\mathrm{v}}_0}{\partial \mathrm{x}}+\frac{\partial \mathrm{w}}{\partial \mathrm{x}}\frac{\partial \mathrm{w}}{\partial \mathrm{y}}\end{array}\right\}\mathrm{dz}}\\ +{\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\left[\begin{array}{ccc}{\mathrm{Q}}_{11}(\mathrm{z})& {\mathrm{Q}}_{12}(\mathrm{z})& 0\\ {\mathrm{Q}}_{12}(\mathrm{z})& {\mathrm{Q}}_{22}(\mathrm{z})& 0\\ 0& 0& {\mathrm{Q}}_{66}(\mathrm{z})\end{array}\right]\left\{\begin{array}{l}-\frac{{\partial }^2{\mathrm{w}}_0}{\partial {\mathrm{x}}^2}\\ -\frac{{\partial }^2{\mathrm{w}}_0}{\partial {\mathrm{y}}^2}\\ -2\frac{{\partial }^2{\mathrm{w}}_0}{\partial \mathrm{x}\partial \mathrm{y}}\end{array}\right\}\mathrm{zdz}}+{\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\left[\begin{array}{ccc}{\mathrm{Q}}_{11}(\mathrm{z})& {\mathrm{Q}}_{12}(\mathrm{z})& 0\\ {\mathrm{Q}}_{12}(\mathrm{z})& {\mathrm{Q}}_{22}(\mathrm{z})& 0\\ 0& 0& {\mathrm{Q}}_{66}(\mathrm{z})\end{array}\right]\left\{\begin{array}{l}\frac{\partial {\mathsf{\theta }}_{\mathrm{x}}}{\partial \mathrm{x}}\\ \frac{\partial {\mathsf{\theta }}_{\mathrm{y}}}{\partial \mathrm{y}}\\ \frac{\partial {\mathsf{\theta }}_{\mathrm{x}}}{\partial \mathrm{y}}+\frac{\partial {\mathsf{\theta }}_{\mathrm{y}}}{\partial \mathrm{x}}\end{array}\right\}\mathrm{f}\left(\mathrm{z}\right)\mathrm{dz}}\\ -{\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\left[\begin{array}{c}{\mathrm{Q}}_{11}(\mathrm{z})+{\mathrm{Q}}_{12}(\mathrm{z})\\ {\mathrm{Q}}_{12}(\mathrm{z})+{\mathrm{Q}}_{22}(\mathrm{z})\\ 0\end{array}\right]\mathsf{\alpha }\left(\mathrm{z}\right)\mathrm{T}\left(\mathrm{z}\right)\mathrm{dz}},\\ \left\{\begin{array}{l}{\mathrm{M}}_{\mathrm{xx}}\\ {\mathrm{M}}_{\mathrm{yy}}\\ {\mathrm{M}}_{\mathrm{xy}}\end{array}\right\}={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\left\{\begin{array}{l}{\mathsf{\sigma }}_{\mathrm{xx}}\\ {\mathsf{\sigma }}_{\mathrm{yy}}\\ {\mathsf{\tau }}_{\mathrm{xy}}\end{array}\right\}\mathrm{zdz}}={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\left[\begin{array}{ccc}{\mathrm{Q}}_{11}(\mathrm{z})& {\mathrm{Q}}_{12}(\mathrm{z})& 0\\ {\mathrm{Q}}_{12}(\mathrm{z})& {\mathrm{Q}}_{22}(\mathrm{z})& 0\\ 0& 0& {\mathrm{Q}}_{66}(\mathrm{z})\end{array}\right]\left\{\begin{array}{l}\frac{\partial {\mathrm{u}}_0}{\partial \mathrm{x}}+\frac{1}{2}{\left(\frac{\partial \mathrm{w}}{\partial \mathrm{x}}\right)}^2\\ \frac{\partial {\mathrm{v}}_0}{\partial \mathrm{y}}+\frac{1}{2}{\left(\frac{\partial \mathrm{w}}{\partial \mathrm{y}}\right)}^2\\ \frac{\partial {\mathrm{u}}_0}{\partial \mathrm{y}}+\frac{\partial {\mathrm{v}}_0}{\partial \mathrm{x}}+\frac{\partial \mathrm{w}}{\partial \mathrm{x}}\frac{\partial \mathrm{w}}{\partial \mathrm{y}}\end{array}\right\}\mathrm{zdz}}\\ +{\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\left[\begin{array}{ccc}{\mathrm{Q}}_{11}(\mathrm{z})& {\mathrm{Q}}_{12}(\mathrm{z})& 0\\ {\mathrm{Q}}_{12}(\mathrm{z})& {\mathrm{Q}}_{22}(\mathrm{z})& 0\\ 0& 0& {\mathrm{Q}}_{66}(\mathrm{z})\end{array}\right]\left\{\begin{array}{l}-\frac{{\partial }^2{\mathrm{w}}_0}{\partial {\mathrm{x}}^2}\\ -\frac{{\partial }^2{\mathrm{w}}_0}{\partial {\mathrm{y}}^2}\\ -2\frac{{\partial }^2{\mathrm{w}}_0}{\partial \mathrm{x}\partial \mathrm{y}}\end{array}\right\}{\mathrm{z}}^2\mathrm{dz}}+{\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\left[\begin{array}{ccc}{\mathrm{Q}}_{11}(\mathrm{z})& {\mathrm{Q}}_{12}(\mathrm{z})& 0\\ {\mathrm{Q}}_{12}(\mathrm{z})& {\mathrm{Q}}_{22}(\mathrm{z})& 0\\ 0& 0& {\mathrm{Q}}_{66}(\mathrm{z})\end{array}\right]\left\{\begin{array}{l}\frac{\partial {\mathsf{\theta }}_{\mathrm{x}}}{\partial \mathrm{x}}\\ \frac{\partial {\mathsf{\theta }}_{\mathrm{y}}}{\partial \mathrm{y}}\\ \frac{\partial {\mathsf{\theta }}_{\mathrm{x}}}{\partial \mathrm{y}}+\frac{\partial {\mathsf{\theta }}_{\mathrm{y}}}{\partial \mathrm{x}}\end{array}\right\}\mathrm{zf}\left(\mathrm{z}\right)\mathrm{dz}}\\ -{\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\left[\begin{array}{c}{\mathrm{Q}}_{11}(\mathrm{z})+{\mathrm{Q}}_{12}(\mathrm{z})\\ {\mathrm{Q}}_{12}(\mathrm{z})+{\mathrm{Q}}_{22}(\mathrm{z})\\ 0\end{array}\right]\mathrm{z}\mathsf{\alpha }\left(\mathrm{z}\right)\mathrm{T}\left(\mathrm{z}\right)\mathrm{dz}},\\ \left\{\begin{array}{l}{\mathrm{P}}_{\mathrm{xx}}\\ {\mathrm{P}}_{\mathrm{yy}}\\ {\mathrm{P}}_{\mathrm{xy}}\end{array}\right\}={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\left\{\begin{array}{l}{\mathsf{\sigma }}_{\mathrm{xx}}\\ {\mathsf{\sigma }}_{\mathrm{yy}}\\ {\mathsf{\tau }}_{\mathrm{xy}}\end{array}\right\}\mathrm{f}\left(\mathrm{z}\right)\mathrm{dz}}={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\left[\begin{array}{ccc}{\mathrm{Q}}_{11}(\mathrm{z})& {\mathrm{Q}}_{12}(\mathrm{z})& 0\\ {\mathrm{Q}}_{12}(\mathrm{z})& {\mathrm{Q}}_{22}(\mathrm{z})& 0\\ 0& 0& {\mathrm{Q}}_{66}(\mathrm{z})\end{array}\right]\left\{\begin{array}{l}\frac{\partial {\mathrm{u}}_0}{\partial \mathrm{x}}+\frac{1}{2}{\left(\frac{\partial \mathrm{w}}{\partial \mathrm{x}}\right)}^2\\ \frac{\partial {\mathrm{v}}_0}{\partial \mathrm{y}}+\frac{1}{2}{\left(\frac{\partial \mathrm{w}}{\partial \mathrm{y}}\right)}^2\\ \frac{\partial {\mathrm{u}}_0}{\partial \mathrm{y}}+\frac{\partial {\mathrm{v}}_0}{\partial \mathrm{x}}+\frac{\partial \mathrm{w}}{\partial \mathrm{x}}\frac{\partial \mathrm{w}}{\partial \mathrm{y}}\end{array}\right\}\mathrm{f}\left(\mathrm{z}\right)\mathrm{dz}}\\ +{\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\left[\begin{array}{ccc}{\mathrm{Q}}_{11}(\mathrm{z})& {\mathrm{Q}}_{12}(\mathrm{z})& 0\\ {\mathrm{Q}}_{12}(\mathrm{z})& {\mathrm{Q}}_{22}(\mathrm{z})& 0\\ 0& 0& {\mathrm{Q}}_{66}(\mathrm{z})\end{array}\right]\left\{\begin{array}{l}-\frac{{\partial }^2{\mathrm{w}}_0}{\partial {\mathrm{x}}^2}\\ -\frac{{\partial }^2{\mathrm{w}}_0}{\partial {\mathrm{y}}^2}\\ -2\frac{{\partial }^2{\mathrm{w}}_0}{\partial \mathrm{x}\partial \mathrm{y}}\end{array}\right\}\mathrm{zf}\left(\mathrm{z}\right)\mathrm{dz}}+{\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\left[\begin{array}{ccc}{\mathrm{Q}}_{11}(\mathrm{z})& {\mathrm{Q}}_{12}(\mathrm{z})& 0\\ {\mathrm{Q}}_{12}(\mathrm{z})& {\mathrm{Q}}_{22}(\mathrm{z})& 0\\ 0& 0& {\mathrm{Q}}_{66}(\mathrm{z})\end{array}\right]\left\{\begin{array}{l}\frac{\partial {\mathsf{\theta }}_{\mathrm{x}}}{\partial \mathrm{x}}\\ \frac{\partial {\mathsf{\theta }}_{\mathrm{y}}}{\partial \mathrm{y}}\\ \frac{\partial {\mathsf{\theta }}_{\mathrm{x}}}{\partial \mathrm{y}}+\frac{\partial {\mathsf{\theta }}_{\mathrm{y}}}{\partial \mathrm{x}}\end{array}\right\}{\left(\mathrm{f}\left(\mathrm{z}\right)\right)}^2\mathrm{dz}}\\ -{\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\left[\begin{array}{c}{\mathrm{Q}}_{11}(\mathrm{z})+{\mathrm{Q}}_{12}(\mathrm{z})\\ {\mathrm{Q}}_{12}(\mathrm{z})+{\mathrm{Q}}_{22}(\mathrm{z})\\ 0\end{array}\right]\mathrm{f}\left(\mathrm{z}\right)\mathsf{\alpha }\left(\mathrm{z}\right)\mathrm{T}\left(\mathrm{z}\right)\mathrm{dz}},\\ \left\{\begin{array}{c}{\mathrm{R}}_{\mathrm{y}}\\ {\mathrm{R}}_{\mathrm{x}}\end{array}\right\}={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\left\{\begin{array}{c}{\mathsf{\tau }}_{\mathrm{xz}}\\ {\mathsf{\tau }}_{\mathrm{yz}}\end{array}\right\}}{\mathrm{f}}^{\prime }\left(\mathrm{z}\right)\mathrm{dz}={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\left[\begin{array}{cc}{\mathrm{Q}}_{44}(\mathrm{z})& 0\\ 0& {\mathrm{Q}}_{55}(\mathrm{z})\end{array}\right]\left\{\begin{array}{c}{\mathsf{\theta }}_{\mathrm{x}}\\ {\mathsf{\theta }}_{\mathrm{y}}\end{array}\right\}{\left({\mathrm{f}}^{\prime }\left(\mathrm{z}\right)\right)}^2\mathrm{dz}},\end{array}$$

(4)

where:

$$\left\{\begin{array}{c}{\mathsf{\epsilon }}_{\mathrm{xx}}^{\left(0\right)}\\ {\mathsf{\epsilon }}_{\mathrm{yy}}^{\left(0\right)}\\ {\mathsf{\gamma }}_{\mathrm{xy}}^{\left(0\right)}\end{array}\right\}=\left\{\begin{array}{c}\frac{\partial {\mathrm{u}}_0}{\partial \mathrm{x}}+\frac{1}{2}{\left(\frac{\partial {\mathrm{w}}_0}{\partial \mathrm{x}}\right)}^2\\ \frac{\partial {\mathrm{v}}_0}{\partial \mathrm{y}}+\frac{1}{2}{\left(\frac{\partial {\mathrm{w}}_0}{\partial \mathrm{y}}\right)}^2\\ \frac{\partial {\mathrm{u}}_0}{\partial \mathrm{y}}+\frac{\partial {\mathrm{v}}_0}{\partial \mathrm{x}}\frac{\partial {\mathrm{w}}_0}{\partial \mathrm{x}}\frac{\partial {\mathrm{w}}_0}{\partial \mathrm{y}}\end{array}\right\}, \left\{\begin{array}{c}{\mathrm{k}}_{\mathrm{xx}}^{\left(0\right)}\\ {\mathrm{k}}_{\mathrm{yy}}^{\left(0\right)}\\ {\mathrm{k}}_{\mathrm{xy}}^{\left(0\right)}\end{array}\right\}=\left\{\begin{array}{c}-\frac{{\partial }^2{\mathrm{w}}_0}{\partial {\mathrm{x}}^2}\\ -\frac{{\partial }^2{\mathrm{w}}_0}{\partial {\mathrm{y}}^2}\\ -2\frac{{\partial }^2{\mathrm{w}}_0}{\partial \mathrm{x}\partial \mathrm{y}}\end{array}\right\}, \left\{\begin{array}{c}{\mathrm{k}}_{\mathrm{xx}}^{\left(1\right)}\\ {\mathrm{k}}_{\mathrm{yy}}^{\left(1\right)}\\ {\mathrm{k}}_{\mathrm{xy}}^{\left(1\right)}\end{array}\right\}=\left\{\begin{array}{c}\frac{\partial {\mathsf{\theta }}_{\mathrm{x}}}{\partial \mathrm{x}}\\ \frac{\partial {\mathsf{\theta }}_{\mathrm{y}}}{\partial \mathrm{y}}\\ \frac{\partial {\mathsf{\theta }}_{\mathrm{x}}}{\partial \mathrm{y}}+\frac{\partial {\mathsf{\theta }}_{\mathrm{y}}}{\partial \mathrm{x}}\end{array}\right\}, \left\{\begin{array}{c}{\mathrm{k}}_{\mathrm{xz}}^{\left(2\right)}\\ {\mathrm{k}}_{\mathrm{yz}}^{\left(2\right)}\end{array}\right\}=\left\{\begin{array}{c}{\mathsf{\theta }}_{\mathrm{x}}\\ {\mathsf{\theta }}_{\mathrm{y}}\end{array}\right\},$$

(5)

and ${\mathrm{f}}^{\prime }(\mathrm{z})=\frac{\mathrm{df}(\mathrm{z})}{\mathrm{dz}}$.

The coefficients of the constitutive elasticity tensor could be derived using engineering constants:

$${\mathrm{Q}}_{11}\left(\mathrm{z}\right)={\mathrm{Q}}_{22}\left(\mathrm{z}\right)=\frac{\mathrm{E}\left(\mathrm{z}\right)}{1-{\mathsf{\nu }}^2}, {\mathrm{Q}}_{44}\left(\mathrm{z}\right)={\mathrm{Q}}_{55}\left(\mathrm{z}\right)={\mathrm{Q}}_{66}\left(\mathrm{z}\right)=\frac{\mathrm{E}\left(\mathrm{z}\right)}{2\left(1+\mathsf{\nu }\right)}, {\mathrm{Q}}_{12}\left(\mathrm{z}\right)=\frac{\mathsf{\nu }\mathrm{E}\left(\mathrm{z}\right)}{1-{\mathsf{\nu }}^2}.$$

(6)

Based on the Equation (4), new matrices are defined:

$$\begin{array}{l}\left({\mathrm{A}}_{\mathrm{ij}},{\mathrm{B}}_{\mathrm{ij}},{ \mathrm{D}}_{\mathrm{ij}},{ \mathrm{E}}_{\mathrm{ij}}, {\mathrm{F}}_{\mathrm{ij}}, {\mathrm{G}}_{\mathrm{ij}}\right)={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}{\mathrm{Q}}_{\mathrm{ij}}\left(1,\mathrm{z},\mathrm{f}\left(\mathrm{z}\right),{\mathrm{z}}^2,\mathrm{zf}\left(\mathrm{z}\right),{\left(\mathrm{f}\left(\mathrm{z}\right)\right)}^2\right)\mathrm{dz}}\mathrm{i},\mathrm{j}=\left(1,2,6\right),\\ {\mathrm{H}}_{\mathrm{lr}}={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}{\mathrm{Q}}_{\mathrm{lr}}^{{}^{}}{\left({\mathrm{f}}^{\prime }\left(\mathrm{z}\right)\right)}^2\mathrm{dz}}, \left(\mathrm{l},\mathrm{r}\right)=\left(4,5\right).\end{array}$$

(7)

In order to use the principles of minimum potential energy, it is necessary to define strain energy U_s, the potential energy of the elastic foundation U_e and the total potential energy Π:

$${\mathrm{U}}_{\mathrm{s}}={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}{\displaystyle \underset{\mathrm{A}}{\int }\left({\mathsf{\sigma }}_{\mathrm{xx}}\left[{\mathsf{\epsilon }}_{\mathrm{xx}}-\mathsf{\alpha }\left(\mathrm{z}\right)\mathrm{T}\left(\mathrm{z}\right)\right]+{\mathsf{\sigma }}_{\mathrm{y}}\left[{\mathsf{\epsilon }}_{\mathrm{yy}}-\mathsf{\alpha }\left(\mathrm{z}\right)\mathrm{T}\left(\mathrm{z}\right)\right]+{\mathsf{\sigma }}_{\mathrm{zz}}{\mathsf{\epsilon }}_{\mathrm{zz}}+{\mathsf{\tau }}_{\mathrm{xy}}{\mathsf{\gamma }}_{\mathrm{xy}}+{\mathsf{\tau }}_{\mathrm{xz}}{\mathsf{\gamma }}_{\mathrm{xz}}+{\mathsf{\tau }}_{\mathrm{yz}}{\mathsf{\gamma }}_{\mathrm{yz}}\right)\mathrm{dAdz}}},$$

(8)

$${\mathrm{U}}_{\mathrm{e}}=\frac{1}{2}{\displaystyle \underset{\mathrm{A}}{\int }\left\{{\mathrm{k}}_0{\mathrm{w}}^2+{\mathrm{k}}_1\left[{\left(\frac{\partial \mathrm{w}}{\partial \mathrm{x}}\right)}^2+{\left(\frac{\partial \mathrm{w}}{\partial \mathrm{x}}\right)}^2\right]\right\}\mathrm{dA}},$$

(9)

$$\mathsf{П}={\mathrm{U}}_{\mathrm{s}}+{\mathrm{U}}_{\mathrm{e}}.$$

(10)

By applying the principles of minimum potential energy:

$$\begin{array}{l}\mathsf{\delta }\mathsf{П}=\mathsf{\delta }\left({\mathrm{U}}_{\mathrm{s}}+{\mathrm{U}}_{\mathrm{e}}\right)=\mathsf{\delta }{\mathrm{U}}_{\mathrm{s}}+\mathsf{\delta }{\mathrm{U}}_{\mathrm{e}}={\displaystyle \underset{\mathrm{A}}{\int }}({\mathrm{N}}_{\mathrm{xx}}\mathsf{\delta }{\mathsf{\epsilon }}_{\mathrm{xx}}^{\left(0\right)}+{\mathrm{N}}_{\mathrm{yy}}\mathsf{\delta }{\mathsf{\epsilon }}_{\mathrm{yy}}^{\left(0\right)}+{\mathrm{N}}_{\mathrm{xy}}\mathsf{\delta }{\mathsf{\gamma }}_{\mathrm{xy}}^{\left(0\right)}+{\mathrm{M}}_{\mathrm{xx}}\mathsf{\delta }{\mathrm{k}}_{\mathrm{xx}}^{\left(0\right)}+{\mathrm{M}}_{\mathrm{yy}}\mathsf{\delta }{\mathrm{k}}_{\mathrm{yy}}^{\left(0\right)}+{\mathrm{M}}_{\mathrm{xy}}\mathsf{\delta }{\mathrm{k}}_{\mathrm{xy}}^{\left(0\right)}\\ +{\mathrm{P}}_{\mathrm{xx}}\mathsf{\delta }{\mathrm{k}}_{\mathrm{xx}}^{\left(1\right)}+{\mathrm{P}}_{\mathrm{yy}}\mathsf{\delta }{\mathrm{k}}_{\mathrm{yy}}^{\left(1\right)}+{\mathrm{P}}_{\mathrm{xy}}\mathsf{\delta }{\mathrm{k}}_{\mathrm{xy}}^{\left(1\right)}+{\mathrm{R}}_{\mathrm{x}}\mathsf{\delta }{\mathrm{k}}_{\mathrm{xz}}^{\left(2\right)}+{\mathrm{R}}_{\mathrm{y}}\mathsf{\delta }{\mathrm{k}}_{\mathrm{yz}}^{\left(2\right)})\mathrm{dA}+{\displaystyle \underset{\mathrm{A}}{\int }\left\{{\mathrm{k}}_0\mathrm{w}\mathsf{\delta }\mathrm{w}+{\mathrm{k}}_1\left(\frac{\partial \mathrm{w}}{\partial \mathrm{x}}\frac{\partial \mathsf{\delta }\mathrm{w}}{\partial \mathrm{x}}+\frac{\partial \mathrm{w}}{\partial \mathrm{y}}\frac{\partial \mathsf{\delta }\mathrm{w}}{\partial \mathrm{y}}\right)\right\}\mathrm{dA}}=0,\end{array}$$

(11)

equilibrium equations become:

$$\begin{array}{l}\mathsf{\delta }{\mathrm{u}}_0: \frac{\partial {\mathrm{N}}_{\mathrm{xx}}}{\partial \mathrm{x}}+\frac{\partial {\mathrm{N}}_{\mathrm{xy}}}{\partial \mathrm{y}}=0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathsf{\delta }{\mathrm{v}}_0: \frac{\partial {\mathrm{N}}_{\mathrm{yy}}}{\partial \mathrm{y}}+\frac{\partial {\mathrm{N}}_{\mathrm{xy}}}{\partial \mathrm{x}}=0, \\ \mathsf{\delta }{\mathrm{w}}_0:\frac{{\partial }^2{\mathrm{M}}_{\mathrm{xx}}}{\partial {\mathrm{x}}^2}+2\frac{{\partial }^2{\mathrm{M}}_{\mathrm{xy}}}{\partial \mathrm{x}\partial \mathrm{y}}+\frac{{\partial }^2{\mathrm{M}}_{\mathrm{yy}}}{\partial {\mathrm{y}}^2}+{\mathrm{N}}_{\mathrm{xx}}\frac{{\partial }^2{\mathrm{w}}_0}{\partial {\mathrm{x}}^2}+2{\mathrm{N}}_{\mathrm{xy}}\frac{{\partial }^2{\mathrm{w}}_0}{\partial \mathrm{x}\partial \mathrm{y}}+{\mathrm{N}}_{\mathrm{yy}}\frac{{\partial }^2{\mathrm{w}}_0}{\partial {\mathrm{y}}^2}-{\mathrm{k}}_0{\mathrm{w}}_0+{\mathrm{k}}_1\left(\frac{{\partial }^2{\mathrm{w}}_0}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2{\mathrm{w}}_0}{\partial {\mathrm{y}}^2}\right)=0,\\ \mathsf{\delta }{\mathsf{\theta }}_{\mathrm{x}}: \frac{\partial {\mathrm{P}}_{\mathrm{xx}}}{\partial \mathrm{x}}+\frac{\partial {\mathrm{P}}_{\mathrm{xy}}}{\partial \mathrm{y}}-{\mathrm{R}}_{\mathrm{x}}=0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathsf{\delta }{\mathsf{\theta }}_{\mathrm{y}}: \frac{\partial {\mathrm{P}}_{\mathrm{xy}}}{\partial \mathrm{x}}+\frac{\partial {\mathrm{P}}_{\mathrm{yy}}}{\partial \mathrm{y}}-{\mathrm{R}}_{\mathrm{y}}=0.\end{array}$$

(12)

Based on the equilibrium Equation (12) and using the displacement components ${\mathrm{u}}_0,\hspace{0.17em}{\mathrm{v}}_0,\hspace{0.17em}{\mathrm{w}}_0,\hspace{0.17em}{\mathsf{\theta }}_{\mathrm{x}0}\hspace{0.17em}\mathrm{and}\hspace{0.17em}{\mathsf{\theta }}_{\mathrm{y}0}$, the stability equation could be defined. The displacement components of the next stable configuration are:

$$\begin{array}{c}\mathrm{u}={\mathrm{u}}_0+{\mathrm{u}}_1,\hspace{1em}\hspace{1em}\mathrm{v}={\mathrm{v}}_0+{\mathrm{v}}_1,\hspace{1em}\hspace{1em}\mathrm{w}={\mathrm{w}}_0+{\mathrm{w}}_1,\\ {\mathsf{\theta }}_{\mathrm{x}}={\mathsf{\theta }}_{\mathrm{x}0}+{\mathsf{\theta }}_{\mathrm{x}1}, {\mathsf{\theta }}_{\mathrm{y}}={\mathsf{\theta }}_{\mathrm{y}0}+{\mathsf{\theta }}_{\mathrm{y}1},\end{array}$$

(13)

where ${\mathrm{u}}_1,\hspace{0.17em}{\mathrm{v}}_1,\hspace{0.17em}{\mathrm{w}}_1,\hspace{0.17em}{\mathsf{\theta }}_{\mathrm{x}1}\hspace{0.17em}\mathrm{and}\hspace{0.17em}{\mathsf{\theta }}_{\mathrm{y}1}$ represent the displacement components of arbitrarily small deviation from the stable configuration. If it is assumed that the temperature is constant in xy-plane of the FG plate and that it changes only in the direction of z-axis, the stability equation can be obtained by substituting the Equation (13) into the Equation (12):

$$\begin{array}{l}\mathsf{\delta }{\mathrm{u}}_0: \frac{\partial {\mathrm{N}}_{\mathrm{xx}}^1}{\partial \mathrm{x}}+\frac{\partial {\mathrm{N}}_{\mathrm{xy}}^1}{\partial \mathrm{y}}=0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathsf{\delta }{\mathrm{v}}_0: \frac{\partial {\mathrm{N}}_{\mathrm{yy}}^1}{\partial \mathrm{y}}+\frac{\partial {\mathrm{N}}_{\mathrm{xy}}^1}{\partial \mathrm{x}}=0, \\ \mathsf{\delta }{\mathrm{w}}_0:\frac{{\partial }^2{\mathrm{M}}_{\mathrm{xx}}^1}{\partial {\mathrm{x}}^2}+2\frac{{\partial }^2{\mathrm{M}}_{\mathrm{xy}}^1}{\partial \mathrm{x}\partial \mathrm{y}}+\frac{{\partial }^2{\mathrm{M}}_{\mathrm{yy}}^1}{\partial {\mathrm{y}}^2}+{\mathrm{N}}_{\mathrm{xx}}^0\frac{{\partial }^2{\mathrm{w}}_0}{\partial {\mathrm{x}}^2}+2{\mathrm{N}}_{\mathrm{xy}}^0\frac{{\partial }^2{\mathrm{w}}_0}{\partial \mathrm{x}\partial \mathrm{y}}+{\mathrm{N}}_{\mathrm{yy}}^0\frac{{\partial }^2{\mathrm{w}}_0}{\partial {\mathrm{y}}^2}-{\mathrm{k}}_0{\mathrm{w}}_1+{\mathrm{k}}_1\left(\frac{{\partial }^2{\mathrm{w}}_1}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2{\mathrm{w}}_1}{\partial {\mathrm{y}}^2}\right)=0,\\ \mathsf{\delta }{\mathsf{\theta }}_{\mathrm{x}}: \frac{\partial {\mathrm{P}}_{\mathrm{xx}}^1}{\partial \mathrm{x}}+\frac{\partial {\mathrm{P}}_{\mathrm{xy}}^1}{\partial \mathrm{y}}-{\mathrm{R}}_{\mathrm{x}}^1=0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{1em}\mathsf{\delta }{\mathsf{\theta }}_{\mathrm{y}}: \frac{\partial {\mathrm{P}}_{\mathrm{xy}}^1}{\partial \mathrm{x}}+\frac{\partial {\mathrm{P}}_{\mathrm{yy}}^1}{\partial \mathrm{y}}-{\mathrm{R}}_{\mathrm{y}}^1=0.\end{array}$$

(14)

where ${\mathrm{N}}_{\mathrm{xx}}^0,{\mathrm{N}}_{\mathrm{yy}}^0\hspace{0.17em}\mathrm{and}\hspace{0.17em}{\mathrm{N}}_{\mathrm{xy}}^0$ are the resultants of pre-buckling forces:

$${\mathrm{N}}_{\mathrm{xx}}^0={\mathrm{N}}_{\mathrm{yy}}^0=-{\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\frac{\mathrm{E}\left(\mathrm{z}\right)\mathsf{\alpha }\left(\mathrm{z}\right)\mathrm{T}\left(\mathrm{z}\right)}{1-\mathsf{\upsilon }}}\mathrm{dz},\hspace{1em}\hspace{1em}{\mathrm{N}}_{\mathrm{xy}}^0=0.$$

(15)

Analytical solutions are obtained by using assumed solution forms and boundary conditions in accordance with Navier’s methods [54]. Boundary conditions along the edges of the rigidly fixed-simply supported plate rectangular plate are the following:

$$\begin{array}{l}{\mathrm{u}}_1={\mathrm{v}}_1={\mathrm{w}}_1={\mathsf{\theta }}_{\mathrm{y}1}={\mathrm{N}}_{\mathrm{xx}}^1={\mathrm{M}}_{\mathrm{xx}}^1={\mathrm{P}}_{\mathrm{xx}}^1=0\hspace{1em}\hspace{1em}\mathrm{along} \mathrm{sides} \mathrm{x}=0 \mathrm{and} \mathrm{x}=\mathrm{a},\\ {\mathrm{u}}_1={\mathrm{v}}_1={\mathrm{w}}_1={\mathsf{\theta }}_{\mathrm{x}1}={\mathrm{N}}_{\mathrm{yy}}^1={\mathrm{M}}_{\mathrm{yy}}^1={\mathrm{P}}_{\mathrm{yy}}^1=0\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{along} \mathrm{sides} \mathrm{y}=0 \mathrm{and} \mathrm{y}=\mathrm{b}.\end{array}$$

(16)

In order to satisfy the previously defined boundary conditions, the following Navier’s solution is assumed [54]:

$$\begin{array}{l}{\mathrm{u}}_1\left(\mathrm{x},\mathrm{y}\right)={\displaystyle \sum _{\mathrm{m}=1}^{\infty }{\displaystyle \sum _{\mathrm{n}=1}^{\infty }{\mathrm{U}}_{\mathrm{mn}}^1}}\cos \frac{\mathrm{m}\mathsf{\pi }\mathrm{x}}{\mathrm{a}}\sin \frac{\mathrm{n}{\pi }\mathrm{y}}{\mathrm{b}},\hspace{1em}\hspace{1em}\hspace{1em}{\mathrm{v}}_1\left(\mathrm{x},\mathrm{y}\right)={\displaystyle \sum _{\mathrm{m}=1}^{\infty }{\displaystyle \sum _{\mathrm{n}=1}^{\infty }{\mathrm{V}}_{\mathrm{mn}}^1}}\sin \frac{\mathrm{m}\mathsf{\pi }\mathrm{x}}{\mathrm{a}}\cos \frac{\mathrm{n}\mathsf{\pi }\mathrm{y}}{\mathrm{b}},\\ {\mathrm{w}}_1\left(\mathrm{x},\mathrm{y}\right)={\displaystyle \sum _{\mathrm{m}=1}^{\infty }{\displaystyle \sum _{\mathrm{n}=1}^{\infty }{\mathrm{W}}_{\mathrm{mn}}^1}}\sin \frac{\mathrm{m}\mathsf{\pi }\mathrm{x}}{\mathrm{a}}\sin \frac{\mathrm{n}\mathsf{\pi }\mathrm{y}}{\mathrm{b}},\\ {\mathsf{\theta }}_{\mathrm{x}1}\left(\mathrm{x},\mathrm{y}\right)={\displaystyle \sum _{\mathrm{m}=1}^{\infty }{\displaystyle \sum _{\mathrm{n}=1}^{\infty }{\mathrm{T}}_{\mathrm{xmn}}^1}}\cos \frac{\mathrm{m}\mathsf{\pi }\mathrm{x}}{\mathrm{a}}\sin \frac{\mathrm{n}\mathsf{\pi }\mathrm{y}}{\mathrm{b}},\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}{\mathsf{\theta }}_{\mathrm{y}1}\left(\mathrm{x},\mathrm{y}\right)={\displaystyle \sum _{\mathrm{m}=1}^{\infty }{\displaystyle \sum _{\mathrm{n}=1}^{\infty }{\mathrm{T}}_{\mathrm{ymn}}^1}}\sin \frac{\mathrm{m}\mathsf{\pi }\mathrm{x}}{\mathrm{a}}\cos \frac{\mathrm{n}\mathsf{\pi }\mathrm{y}}{\mathrm{b}}.\end{array}$$

(17)

where ${\mathrm{U}}_{\mathrm{mn}}^1,{\mathrm{V}}_{\mathrm{mn}}^1,{\mathrm{W}}_{\mathrm{mn}}^1,{\mathrm{T}}_{\mathrm{xmn}}^1,{\mathrm{T}}_{\mathrm{ymn}}^1$ are parameters which are to be determined.

Based on Navier’s solution, the stability equation becomes:

$$\left[L-\mathsf{\Omega }I\right]U=0,$$

(18)

where $U={\left\{\begin{array}{ccccc}{\mathrm{U}}_{\mathrm{mn}}^1& {\mathrm{V}}_{\mathrm{mn}}^1& {\mathrm{W}}_{\mathrm{mn}}^1& {\mathrm{T}}_{\mathrm{xmn}}^1& {\mathrm{T}}_{\mathrm{ymn}}^1\end{array}\right\}}^{\mathrm{T}}$ and $\mathsf{\Omega }$ buckling parameter. Coefficients ${\mathrm{L}}_{\mathrm{ij}},$ $\left(\mathrm{i},\mathrm{j}=1÷5\right)$ are defined in the following way:

$$\begin{array}{l}{\mathrm{L}}_{11}={\mathsf{\alpha }}^2{\mathrm{A}}_{11}+{\mathsf{\beta }}^2{\mathrm{A}}_{66}, {\mathrm{L}}_{12}=\mathsf{\alpha }\mathsf{\beta }\left({\mathrm{A}}_{12}+{\mathrm{A}}_{66}\right), {\mathrm{L}}_{13}=0, {\mathrm{L}}_{14}={\mathsf{\alpha }}^2{\mathrm{D}}_{11}+{\mathsf{\beta }}^2{\mathrm{D}}_{66},\\ {\mathrm{L}}_{15}=\mathsf{\alpha }\mathsf{\beta }\left({\mathrm{D}}_{12}+{\mathrm{D}}_{16}\right), {\mathrm{L}}_{22}={\mathsf{\alpha }}^2{\mathrm{A}}_{66}+{\mathsf{\beta }}^2{\mathrm{A}}_{22}, {\mathrm{L}}_{23}=0, {\mathrm{L}}_{24}=\mathsf{\alpha }\mathsf{\beta }\left({\mathrm{D}}_{12}+{\mathrm{D}}_{16}\right),\\ {\mathrm{L}}_{25}={\mathsf{\alpha }}^2{\mathrm{C}}_{66}+{\mathsf{\beta }}^2{\mathrm{C}}_{22}, {\mathrm{L}}_{33}={\mathsf{\alpha }}^4{\mathrm{E}}_{11}+2{\mathsf{\alpha }}^2{\mathsf{\beta }}^2{\mathrm{E}}_{12}+4{\mathsf{\alpha }}^2{\mathsf{\beta }}^2{\mathrm{E}}_{66}+{\mathsf{\beta }}^4{\mathrm{E}}_{22},\\ {\mathrm{L}}_{34}=-{\mathsf{\alpha }}^3{\mathrm{F}}_{11}-\mathsf{\alpha }{\mathsf{\beta }}^2{\mathrm{F}}_{12}-2\mathsf{\alpha }{\mathsf{\beta }}^2{\mathrm{E}}_{66}, {\mathrm{L}}_{35}=-{\mathsf{\alpha }}^2\mathsf{\beta }{\mathrm{F}}_{12}-2{\mathsf{\alpha }}^2\mathsf{\beta }{\mathrm{F}}_{66}-{\mathsf{\beta }}^3{\mathrm{F}}_{22},\\ {\mathrm{L}}_{44}={\mathrm{H}}_{44}+{\mathsf{\alpha }}^2{\mathrm{G}}_{11}+{\mathsf{\alpha }}^2{\mathrm{G}}_{66}, {\mathrm{L}}_{45}=\mathsf{\alpha }\mathsf{\beta }\left({\mathrm{G}}_{12}+{\mathrm{G}}_{66}\right), {\mathrm{L}}_{55}={\mathrm{H}}_{55}+{\mathsf{\alpha }}^2{\mathrm{G}}_{66}+{\mathsf{\alpha }}^2{\mathrm{G}}_{22},\end{array}$$

(19)

while the matrix ${\mathrm{I}}_{\mathrm{ij}}, \left(\mathrm{i},\mathrm{j}=1÷5\right)$ is determined as:

$$I=\{\begin{array}{l}0,\\ {\mathsf{\alpha }}^2{\mathrm{N}}_{\mathrm{x}}^0+{\mathsf{\beta }}^2{\mathrm{N}}_{\mathrm{y}}^0+{\mathrm{k}}_0+{\mathrm{k}}_1\left({\mathsf{\alpha }}^2+{\mathsf{\beta }}^2\right), \left(\mathrm{i},\mathrm{j}=3\right), \mathsf{\alpha }=\mathrm{m}\mathsf{\pi }/\mathrm{a}, \mathsf{\beta }=\mathrm{n}\mathsf{\pi }/\mathrm{b} \end{array}$$

(20)

The determinant in (18) must be equal to zero value in order to get nontrivial solutions:

$$\left|L-\mathsf{\Omega }I\right|=0.$$

(21)

4. Equations of Motion of FG Plate Placed on Elastic Foundation

As the subject of this chapter is linear dynamic analysis, the kinematic relations of displacements and deformations are defined under assumptions of small deformations. Since the total potential energy $\mathsf{П}$ is represented as the sum of the strain energy of the plate and the potential energy of the elastic foundation (10), for the application of the Hamilton’s principle it is necessary to further define the kinetic energy:

$$\mathrm{K}=\frac{1}{2}{\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}{\displaystyle \underset{\mathrm{A}}{\int }\mathsf{\rho }(\mathrm{z})\left[{\left(\frac{\partial \mathrm{u}}{\partial \mathrm{t}}\right)}^2+{\left(\frac{\partial \mathrm{v}}{\partial \mathrm{t}}\right)}^2+{\left(\frac{\partial \mathrm{w}}{\partial \mathrm{t}}\right)}^2\right]}}\mathrm{dAdz},$$

(22)

wherein:

• $\mathsf{\rho }(\mathrm{z})$-material density in an arbitrary cross-section z,

As can be seen, as a consequence of the gradient structure of the plate material, the material density represents a function of the z coordinate. The change of density in the direction of z-axis is defined in accordance to the power-law distribution as:

$$\mathsf{\rho }\left(\mathrm{z}\right)={\mathsf{\rho }}_{\mathrm{m}}+{\mathsf{\rho }}_{\mathrm{cm}}{\left(\frac{1}{2}+\frac{\mathrm{z}}{\mathrm{h}}\right)}^{\mathrm{p}}, {\mathsf{\rho }}_{\mathrm{cm}}={\mathsf{\rho }}_{\mathrm{c}}-{\mathsf{\rho }}_{\mathrm{m}},$$

(23)

By substituting the strain energy of the plate (8), the potential energy of the elastic foundation (9) and the kinetic energy (22) into the Hamilton’s principle:

$$\begin{array}{l}\mathsf{\delta }{\displaystyle \underset{{\mathrm{t}}_1}{\overset{{\mathrm{t}}_1}{\int }}\left[\mathrm{K}-\left({\mathrm{U}}_{\mathrm{s}}+{\mathrm{U}}_{\mathrm{e}}\right)\right]}\mathrm{dt}={\displaystyle \underset{{\mathrm{t}}_1}{\overset{{\mathrm{t}}_2}{\int }}[-{\displaystyle \underset{\mathrm{A}}{\int }({\mathrm{N}}_{\mathrm{xx}}{\mathsf{\epsilon }}_{\mathrm{xx}}^{\left(0\right)}+{\mathrm{N}}_{\mathrm{yy}}{\mathsf{\epsilon }}_{\mathrm{yy}}^{\left(0\right)}+{\mathrm{N}}_{\mathrm{xy}}{\mathsf{\gamma }}_{\mathrm{xy}}^{\left(0\right)}+{\mathrm{M}}_{\mathrm{xx}}{\mathrm{k}}_{\mathrm{xx}}^{\left(0\right)}+{\mathrm{M}}_{\mathrm{yy}}{\mathrm{k}}_{\mathrm{yy}}^{\left(0\right)}+{\mathrm{M}}_{\mathrm{xy}}{\mathrm{k}}_{\mathrm{xy}}^{\left(0\right)}}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{1em}\hspace{1em}\hspace{1em}+{\mathrm{P}}_{\mathrm{xx}}{\mathrm{k}}_{\mathrm{xx}}^{\left(1\right)}+{\mathrm{P}}_{\mathrm{yy}}{\mathrm{k}}_{\mathrm{yy}}^{\left(1\right)}+{\mathrm{P}}_{\mathrm{xy}}{\mathrm{k}}_{\mathrm{xy}}^{\left(1\right)}+{\mathrm{R}}_{\mathrm{x}}{\mathrm{k}}_{\mathrm{xz}}^{\left(2\right)}+{\mathrm{R}}_{\mathrm{y}}{\mathrm{k}}_{\mathrm{yz}}^{\left(2\right)})\mathrm{dA}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}-{\displaystyle \underset{\mathrm{A}}{\int }\left\{{\mathrm{k}}_0\mathrm{w}\mathsf{\delta }\mathrm{w}+{\mathrm{k}}_1\left(\frac{\partial \mathrm{w}}{\partial \mathrm{x}}\frac{\partial \mathsf{\delta }\mathrm{w}}{\partial \mathrm{x}}+\frac{\partial \mathrm{w}}{\partial \mathrm{y}}\frac{\partial \mathsf{\delta }\mathrm{w}}{\partial \mathrm{y}}\right)\right\}\mathrm{dA}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}{\displaystyle \underset{\mathrm{A}}{\int }\mathsf{\rho }(\mathrm{z})\left(\frac{\partial \mathrm{u}}{\partial \mathrm{t}}\frac{\partial \mathsf{\delta }\mathrm{u}}{\partial \mathrm{t}}+\frac{\partial \mathrm{v}}{\partial \mathrm{t}}\frac{\partial \mathsf{\delta }\mathrm{v}}{\partial \mathrm{t}}+\frac{\partial \mathrm{w}}{\partial \mathrm{t}}\frac{\partial \mathsf{\delta }\mathrm{w}}{\partial \mathrm{t}}\right)}}\mathrm{dAdz}]\mathrm{dt}=0.\end{array}$$

(24)

By substituting the strain components expressed by assumed displacement forms as well as by applying the calculus of variations and group the members with the $\mathsf{\delta }{\mathrm{u}}_0,\hspace{0.17em}\mathsf{\delta }{\mathrm{v}}_0,\hspace{0.17em}\mathsf{\delta }{\mathrm{w}}_0,\hspace{0.17em}\mathsf{\delta }{\mathsf{\theta }}_{\mathrm{x}}\hspace{0.17em}\mathrm{and}\hspace{0.17em}\mathsf{\delta }{\mathsf{\theta }}_{\mathrm{y}}$, the equations of motion are obtained:

$$\begin{array}{l}\mathsf{\delta }{\mathrm{u}}_0: \frac{\partial {\mathrm{N}}_{\mathrm{xx}}}{\partial \mathrm{x}}+\frac{\partial {\mathrm{N}}_{\mathrm{xy}}}{\partial \mathrm{y}}={\mathrm{I}}_1{\ddot{\mathrm{u}}}_0-{\mathrm{I}}_2\frac{\partial {\ddot{\mathrm{w}}}_0}{\partial \mathrm{x}}+{\mathrm{I}}_4{\ddot{\mathsf{\theta }}}_{\mathrm{x}},\\ \mathsf{\delta }{\mathrm{v}}_0: \frac{\partial {\mathrm{N}}_{\mathrm{yy}}}{\partial \mathrm{y}}+\frac{\partial {\mathrm{N}}_{\mathrm{xy}}}{\partial \mathrm{x}}={\mathrm{I}}_1{\ddot{\mathrm{v}}}_0-{\mathrm{I}}_2\frac{\partial {\ddot{\mathrm{w}}}_0}{\partial \mathrm{y}}+{\mathrm{I}}_4{\ddot{\mathsf{\theta }}}_{\mathrm{y}}, \\ \mathsf{\delta }{\mathrm{w}}_0: \frac{{\partial }^2{\mathrm{M}}_{\mathrm{xx}}}{\partial {\mathrm{x}}^2}+2\frac{{\partial }^2{\mathrm{M}}_{\mathrm{xy}}}{\partial \mathrm{x}\partial \mathrm{y}}+\frac{{\partial }^2{\mathrm{M}}_{\mathrm{yy}}}{\partial {\mathrm{y}}^2}-{\mathrm{k}}_0{\mathrm{w}}_0+{\mathrm{k}}_1\left(\frac{{\partial }^2{\mathrm{w}}_0}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2{\mathrm{w}}_0}{\partial {\mathrm{y}}^2}\right)={\mathrm{I}}_1{\ddot{\mathrm{w}}}_0+{\mathrm{I}}_2\left(\frac{\partial {\ddot{\mathrm{u}}}_0}{\partial \mathrm{x}}+\frac{\partial {\ddot{\mathrm{v}}}_0}{\partial \mathrm{y}}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}-{\mathrm{I}}_3\left(\frac{{\partial }^2{\ddot{\mathrm{w}}}_0}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2{\ddot{\mathrm{w}}}_0}{\partial {\mathrm{y}}^2}\right)+{\mathrm{I}}_5\left(\frac{\partial {\ddot{\mathsf{\theta }}}_{\mathrm{x}}}{\partial \mathrm{x}}+\frac{\partial {\ddot{\mathsf{\theta }}}_{\mathrm{y}}}{\partial \mathrm{y}}\right),\\ \mathsf{\delta }{\mathsf{\theta }}_{\mathrm{x}}: \frac{\partial {\mathrm{P}}_{\mathrm{xx}}}{\partial \mathrm{x}}+\frac{\partial {\mathrm{P}}_{\mathrm{xy}}}{\partial \mathrm{y}}-{\mathrm{R}}_{\mathrm{x}}={\mathrm{I}}_4{\ddot{\mathrm{u}}}_0-{\mathrm{I}}_5\frac{\partial {\ddot{\mathrm{w}}}_0}{\partial \mathrm{x}}+{\mathrm{I}}_6{\ddot{\mathsf{\theta }}}_{\mathrm{x}},\\ \mathsf{\delta }{\mathsf{\theta }}_{\mathrm{y}}: \frac{\partial {\mathrm{P}}_{\mathrm{xy}}}{\partial \mathrm{x}}+\frac{\partial {\mathrm{P}}_{\mathrm{yy}}}{\partial \mathrm{y}}-{\mathrm{R}}_y={\mathrm{I}}_4{\ddot{\mathrm{v}}}_0-{\mathrm{I}}_5\frac{\partial {\ddot{\mathrm{w}}}_0}{\partial \mathrm{y}}+{\mathrm{I}}_6{\ddot{\mathsf{\theta }}}_{\mathrm{y}}.\end{array}$$

(25)

where ${\mathrm{I}}_{\mathrm{i}}\hspace{0.17em}(\mathrm{i}=1,2,3,4,5\hspace{0.17em}и\hspace{0.17em}6)$ are terms due to inertia defined as:

$$\begin{array}{l}{\mathrm{I}}_1={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\mathsf{\rho }(\mathrm{z})\mathrm{dz}},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}{\mathrm{I}}_2={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\mathsf{\rho }(\mathrm{z})\mathrm{zdz}},\\ {\mathrm{I}}_3={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\mathsf{\rho }(\mathrm{z})\mathrm{f}\left(\mathrm{z}\right)\mathrm{dz}},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathrm{I}}_4={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\mathsf{\rho }(z){\mathrm{z}}^2\mathrm{dz}},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\\ {\mathrm{I}}_5={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\mathsf{\rho }(\mathrm{z})\mathrm{zf}\left(\mathrm{z}\right)\mathrm{dz}},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{1em}{\mathrm{I}}_6={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}\mathsf{\rho }(\mathrm{z}){\left(\mathrm{f}\left(\mathrm{z}\right)\right)}^2\mathrm{dz}}.\end{array}$$

(26)

Analytical solutions will be derived for the simply supported rectangular FGM plate, wherein the boundary conditions are defined according to [54] as:

$$\begin{array}{l}{\mathrm{v}}_0={\mathrm{w}}_0={\mathsf{\theta }}_{\mathrm{y}}={\mathrm{N}}_{\mathrm{xx}}={\mathrm{M}}_{\mathrm{xx}}={\mathrm{P}}_{\mathrm{xx}}={\mathrm{R}}_{\mathrm{y}}=0, \mathrm{on} \mathrm{the} \mathrm{sides} \mathrm{where} \mathrm{x}=0 \mathrm{or} \mathrm{x}=\mathrm{a},\\ {\mathrm{u}}_0={\mathrm{w}}_0={\mathsf{\theta }}_{\mathrm{x}}={\mathrm{N}}_{\mathrm{yy}}={\mathrm{M}}_{\mathrm{yy}}={\mathrm{P}}_{\mathrm{yy}}={\mathrm{R}}_{\mathrm{x}}=0, \mathrm{on} \mathrm{the} \mathrm{sides} \mathrm{where} \mathrm{y}=0 \mathrm{or} \mathrm{y}=\mathrm{b}.\end{array}$$

(27)

In order to satisfy the previously defined boundary conditions, the following Navier’s solution is assumed:

$$\begin{array}{l}{\mathrm{u}}_0\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)={\displaystyle \sum _{\mathrm{m}=1}^{\infty }{\displaystyle \sum _{\mathrm{n}=1}^{\infty }{\mathrm{U}}_{\mathrm{mn}}}}\cos \frac{\mathrm{m}\mathsf{\pi }\mathrm{x}}{\mathrm{a}}\sin \frac{\mathrm{n}\mathsf{\pi }\mathrm{y}}{\mathrm{b}}{\mathrm{e}}^{\mathrm{i}\mathsf{\omega }\mathrm{t}}, \hspace{0.17em}{\mathrm{v}}_0\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)={\displaystyle \sum _{\mathrm{m}=1}^{\infty }{\displaystyle \sum _{\mathrm{n}=1}^{\infty }{\mathrm{V}}_{\mathrm{mn}}}}\sin \frac{\mathrm{m}\mathsf{\pi }\mathrm{x}}{\mathrm{a}}\cos \frac{\mathrm{n}\mathsf{\pi }\mathrm{y}}{\mathrm{b}}{\mathrm{e}}^{\mathrm{i}\mathsf{\omega }\mathrm{t}},\\ {\mathrm{w}}_0\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)={\displaystyle \sum _{\mathrm{m}=1}^{\infty }{\displaystyle \sum _{\mathrm{n}=1}^{\infty }{\mathrm{W}}_{\mathrm{mn}}}}\sin \frac{\mathrm{m}\mathsf{\pi }\mathrm{x}}{\mathrm{a}}\sin \frac{\mathrm{n}\mathsf{\pi }\mathrm{y}}{\mathrm{b}}{\mathrm{e}}^{\mathrm{i}\mathsf{\omega }\mathrm{t}},\\ {\mathsf{\theta }}_{\mathrm{x}}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)={\displaystyle \sum _{\mathrm{m}=1}^{\infty }{\displaystyle \sum _{\mathrm{n}=1}^{\infty }{\mathrm{T}}_{\mathrm{xmn}}}}\cos \frac{\mathrm{m}\mathsf{\pi }\mathrm{x}}{\mathrm{a}}\sin \frac{\mathrm{n}\mathsf{\pi }\mathrm{y}}{\mathrm{b}}{\mathrm{e}}^{\mathrm{i}\mathsf{\omega }\mathrm{t}}, \hspace{0.17em}{\mathsf{\theta }}_{\mathrm{y}}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)={\displaystyle \sum _{\mathrm{m}=1}^{\infty }{\displaystyle \sum _{\mathrm{n}=1}^{\infty }{\mathrm{T}}_{\mathrm{ymn}}}}\sin \frac{\mathrm{m}\mathsf{\pi }\mathrm{x}}{\mathrm{a}}\cos \frac{\mathrm{n}\mathsf{\pi }\mathrm{y}}{\mathrm{b}}{\mathrm{e}}^{\mathrm{i}\mathsf{\omega }\mathrm{t}}.\end{array}$$

(28)

Comparing the assumed form of Navier’s solutions (17) and (28), it can be observed that the only difference between these forms is in the terms ${\mathrm{e}}^{\mathrm{i}\mathsf{\omega }\mathrm{t}}$, where $\mathsf{\omega }$ is the natural frequency of the system and ${\mathrm{U}}_{\mathrm{mn}},{\mathrm{V}}_{\mathrm{mn}},{\mathrm{W}}_{\mathrm{mn}},{\mathrm{T}}_{\mathrm{xmn}},{\mathrm{T}}_{\mathrm{ymn}}$ are parameters to be determined. By substituting (28) into (25) the following equation is obtained:

$$\left[K-{\mathsf{\omega }}^{2}I\right]\overline{U}=0,$$

(29)

$$K=\left[\begin{array}{ccccc}{\mathrm{K}}_{11}& {\mathrm{K}}_{12}& {\mathrm{K}}_{13}& {\mathrm{K}}_{14}& {\mathrm{K}}_{15}\\ & {\mathrm{K}}_{22}& {\mathrm{K}}_{23}& {\mathrm{K}}_{24}& {\mathrm{K}}_{25}\\ & & {\mathrm{K}}_{33}& {\mathrm{K}}_{34}& {\mathrm{K}}_{35}\\ & \mathrm{sym}& & {\mathrm{K}}_{44}& {\mathrm{K}}_{45}\\ & & & & {\mathrm{K}}_{55}\end{array}\right],\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\overline{U}=\left\{\begin{array}{c}{\mathrm{U}}_{\mathrm{mn}}\\ {\mathrm{V}}_{\mathrm{mn}}\\ {\mathrm{W}}_{\mathrm{mn}}\\ {\mathrm{T}}_{\mathrm{xmn}}\\ {\mathrm{T}}_{\mathrm{ymn}}\end{array}\right\}.$$

(30)

Coefficients ${\mathrm{K}}_{\mathrm{ij}},$$\left(\mathrm{i},\mathrm{j}=1÷5\right)$ are defined as:

$$\begin{array}{l}{\mathrm{K}}_{11}={\mathsf{\alpha }}^2{\mathrm{A}}_{11}+{\mathsf{\beta }}^2{\mathrm{A}}_{66}, {\mathrm{K}}_{12}=\mathsf{\alpha }\mathsf{\beta }\left({\mathrm{A}}_{12}+{\mathrm{A}}_{66}\right), {\mathrm{K}}_{13}=-3{\mathrm{B}}_{16}{\mathsf{\alpha }}^2\mathsf{\beta }-{\mathrm{B}}_{26}{\mathsf{\beta }}^3,\\ {\mathrm{K}}_{14}=2{\mathrm{D}}_{16}\mathsf{\alpha }\mathsf{\beta }, {\mathrm{K}}_{15}={\mathsf{\alpha }}^2{\mathrm{D}}_{16}+{\mathsf{\beta }}^2{\mathrm{D}}_{26}, {\mathrm{K}}_{22}={\mathsf{\alpha }}^2{\mathrm{A}}_{66}+{\mathsf{\beta }}^2{\mathrm{A}}_{22},\\ {\mathrm{K}}_{23}=-{\mathrm{B}}_{16}{\mathsf{\alpha }}^3-3{\mathrm{B}}_{26}\mathsf{\alpha }{\mathsf{\beta }}^2, {\mathrm{K}}_{24}={\mathsf{\alpha }}^2{\mathrm{E}}_{16}+{\mathsf{\beta }}^2{\mathrm{E}}_{26}, {\mathrm{K}}_{25}=2\mathsf{\alpha }\mathsf{\beta }{\mathrm{E}}_{26},\\ {\mathrm{K}}_{33}={\mathsf{\alpha }}^4{\mathrm{E}}_{11}+2{\mathsf{\alpha }}^2{\mathsf{\beta }}^2{\mathrm{E}}_{12}+4{\mathsf{\alpha }}^2{\mathsf{\beta }}^2{\mathrm{E}}_{66}+{\mathsf{\beta }}^4{\mathrm{E}}_{22}, {\mathrm{K}}_{34}=-{\mathsf{\alpha }}^3{\mathrm{F}}_{11}-\mathsf{\alpha }{\mathsf{\beta }}^2{\mathrm{F}}_{12}-2\mathsf{\alpha }{\mathsf{\beta }}^2{\mathrm{F}}_{66},\\ {\mathrm{K}}_{35}=-{\mathsf{\alpha }}^2\mathsf{\beta }{\mathrm{F}}_{12}-2{\mathsf{\alpha }}^2\mathsf{\beta }{\mathrm{F}}_{66}-{\mathsf{\beta }}^3{\mathrm{F}}_{22}, {\mathrm{K}}_{44}={\mathrm{H}}_{44}+{\mathsf{\alpha }}^2{\mathrm{G}}_{11}+{\mathsf{\beta }}^2{\mathrm{G}}_{66},\\ {\mathrm{K}}_{45}=\mathsf{\alpha }\mathsf{\beta }\left({\mathrm{G}}_{12}+{\mathrm{G}}_{66}\right), {\mathrm{K}}_{55}={\mathrm{H}}_{55}+{\mathsf{\alpha }}^2{\mathrm{G}}_{66}+{\mathsf{\beta }}^2{\mathrm{G}}_{22}.\end{array}$$

(31)

while ${\mathrm{I}}_{\mathrm{ij}}, \left(\mathrm{i},\mathrm{j}=1÷5\right)$ is defined as:

$$I=\left[\begin{array}{ccccc}{\mathrm{I}}_1& 0& -\mathsf{\alpha }{\mathrm{I}}_2& {\mathrm{I}}_4& 0\\ 0& {\mathrm{I}}_1& -\mathsf{\beta }{\mathrm{I}}_2& 0& {\mathrm{I}}_4\\ -\mathsf{\alpha }{\mathrm{I}}_2& -\mathsf{\beta }{\mathrm{I}}_2& {\mathrm{I}}_3\left({\mathsf{\alpha }}^2+{\mathsf{\beta }}^2\right)+{\mathrm{I}}_1& -\mathsf{\alpha }{\mathrm{I}}_5& -\mathsf{\beta }{\mathrm{I}}_5\\ {\mathrm{I}}_4& 0& -\mathsf{\alpha }{\mathrm{I}}_5& {\mathrm{I}}_6& 0\\ 0& {\mathrm{I}}_4& -\mathsf{\beta }{\mathrm{I}}_5& 0& {\mathrm{I}}_6\end{array}\right], \mathsf{\alpha }=\frac{\mathrm{m}\mathsf{\pi }}{\mathrm{a}}, \mathsf{\beta }=\frac{\mathrm{n}\mathsf{\pi }}{\mathrm{b}}$$

(32)

The determinant in (29) must be equal to zero value in order to get nontrivial solutions:

$$\left|K-{\mathsf{\omega }}^{2}I\right|=0.$$

(33)

5. Numerical Examples and Results

In order to verify the derived theoretical results on numerical examples, an original program code for determination of critical buckling temperature as well as for determination natural frequencies of FGM plates has been implemented within software MATLAB. The main goal of this chapter is to check the exactness and the effectiveness of the proposed theory and new shape function. Different numerical examples were done in order to compare obtained results to results based on other shear deformation theory. Analysis were performed for FG plate with metal and ceramic constituents, the mechanical and thermal characteristics of which are given in

Table 2. Material properties of functionally graded materials (FGM) constituents.

Material	Material Properties
	Elasticity Modulus $\mathrm{E}[\mathrm{GPa}]$	Poisson’s Ratio $\mathsf{\nu }$	Thermal Expansion Coefficient $\mathsf{\alpha }[°{\mathrm{C}}^{-1}]$	Density $\mathsf{\rho }\hspace{0.17em}[\mathrm{kg}/{\mathrm{m}}^3]$
Aluminum ($\mathrm{Al}$)	${\mathrm{E}}_{\mathrm{m}}=70$	$\mathsf{\nu }=0.3$	${\mathsf{\alpha }}_{\mathrm{M}}=23\cdot {10}^{-6}[°{\mathrm{C}}^{-1}]$	${\mathsf{\rho }}_{\mathrm{m}}=2702$
Alumina ($\mathrm{A}{\mathrm{l}}_2{\mathrm{O}}_3$)	${\mathrm{E}}_{\mathrm{c}}=380$	$\mathsf{\nu }=0.3$	${\mathsf{\alpha }}_{\mathrm{C}}=7.4\cdot {10}^{-6}[°{\mathrm{C}}^{-1}]$	${\mathsf{\rho }}_{\mathrm{c}}=3800$

.

5.1. Thermal Buckling Analysis

This section provides results of the research based on comparative analysis of all shape functions (Table 1) and new proposed function. The obtained results for critical buckling temperatures of FG square and rectangular plates placed on an elastic foundation are totally in accordance with the results by authors in [55] which is applied trigonometric shear deformation plate theory and [25] which is applied HSDT based on the just one shape function.

In

Table 3. Critical buckling temperatures (∆t_cr) of FGM plates placed on an elastic foundation for case of linear temperature change across plate thickness (a/b = 1, m = n = 1, and T_m = 5°).

p	Source	∆tcr
		k0 = 0, k1 = 0		k0 = 10, k1 = 0		k0 = 10, k1 = 10
		a/h = 10	a/h = 20	a/h = 10	a/h = 20	a/h = 10	a/h = 20
0	[55]	3.2276	0.833	3.3154	0.855	5.0479	1.2881
	[25]	3.2273	0.833	3.3151	0.855	5.0476	1.2881
	Present study	3.2274	0.8331	3.3151	0.855	5.0476	1.2881
	No. 1	3.2273	0.833	3.3151	0.855	5.0476	1.2881
	No. 2	3.2273	0.833	3.3151	0.855	5.0476	1.2881
	No. 3	3.2276	0.833	3.3154	0.855	5.0479	1.2881
	No. 4	3.2333	0.8334	3.3211	0.8554	5.0536	1.2885
	No. 5	3.2273	0.833	3.3151	0.855	5.0476	1.2881
	No. 6	3.2282	0.8331	3.316	0.855	5.0485	1.2882
	No. 7	3.2284	0.8331	3.3162	0.855	5.0487	1.2882
	No. 8	3.2285	0.8331	3.3163	0.855	5.0488	1.2882
	No. 9	3.2285	0.8331	3.3163	0.855	5.0488	1.2882
	No. 10	3.2273	0.833	3.3151	0.855	5.0476	1.2881
	No. 11	3.2288	0.8331	3.3166	0.8551	5.0491	1.2882
	No. 12	3.2275	0.833	3.3152	0.855	5.0477	1.2881
	No. 13	3.2273	0.833	3.3151	0.855	5.0476	1.2881
1	[55]	1.413	0.3587	1.4897	0.3778	3.004	0.7564
	[25]	1.4129	0.3587	1.4896	0.3778	3.0039	0.7564
	Present study	1.413	0.3587	1.4897	0.3779	3.0039	0.7564
	No. 1	1.4129	0.3587	1.4896	0.3778	3.0039	0.7564
	No. 2	1.4129	0.3587	1.4896	0.3778	3.0039	0.7564
	No. 3	1.413	0.3587	1.4897	0.3778	3.004	0.7564
	No. 4	1.4152	0.3588	1.4919	0.378	3.0061	0.7565
	No. 5	1.4129	0.3587	1.4896	0.3778	3.0039	0.7564
	No. 6	1.4133	0.3587	1.49	0.3779	3.0042	0.7564
	No. 7	1.4133	0.3587	1.49	0.3779	3.0043	0.7564
	No. 8	1.4134	0.3587	1.4901	0.3779	3.0043	0.7564
	No. 9	1.4134	0.3587	1.4901	0.3779	3.0043	0.7564
	No. 10	1.4129	0.3587	1.4896	0.3778	3.0039	0.7564
	No. 11	1.4135	0.3587	1.4902	0.3779	3.0044	0.7564
	No. 12	1.413	0.3587	1.4897	0.3778	3.0039	0.7564
	No. 13	1.4129	0.3587	1.4896	0.3778	3.0039	0.7564
5	[55]	1.16	0.2986	1.2576	0.323	3.1839	0.8046
	[25]	1.1606	0.2987	1.2582	0.3231	3.1845	0.8046
	Present study	1.1608	0.2987	1.2584	0.3231	3.1846	0.8047
	No. 1	1.1606	0.2987	1.2582	0.3231	3.1845	0.8046
	No. 2	1.1606	0.2987	1.2582	0.3231	3.1845	0.8046
	No. 3	1.16	0.2986	1.2576	0.323	3.1839	0.8046
	No. 4	1.1604	0.2986	1.2579	0.323	3.1842	0.8046
	No. 5	1.1607	0.2987	1.2582	0.3231	3.1845	0.8046
	No. 6	1.1626	0.2988	1.2602	0.3232	3.1865	0.8048
	No. 7	1.1597	0.2986	1.2573	0.323	3.1835	0.8045
	No. 8	1.1597	0.2986	1.2572	0.323	3.1835	0.8045
	No. 9	1.1597	0.2986	1.2572	0.323	3.1835	0.8045
	No. 10	1.1607	0.2987	1.2583	0.3231	3.1846	0.8046
	No. 11	1.1631	0.2988	1.2607	0.3232	3.187	0.8048
	No. 12	1.1602	0.2986	1.2578	0.323	3.184	0.8046
	No. 13	1.1606	0.2987	1.2582	0.323	3.1844	0.8046
10	[55]	1.2183	0.3156	1.3317	0.344	3.5699	0.9035
	[25]	1.2186	0.3156	1.332	0.344	3.5701	0.9035
	Present study	1.2187	0.3157	1.3321	0.344	3.5703	0.9036
	No. 1	1.2186	0.3156	1.332	0.344	3.5701	0.9035
	No. 2	1.2186	0.3156	1.332	0.344	3.5701	0.9035
	No. 3	1.2183	0.3156	1.3317	0.344	3.5699	0.9035
	No. 4	1.2206	0.3158	1.3339	0.3441	3.5721	0.9037
	No. 5	1.2186	0.3156	1.332	0.344	3.5702	0.9035
	No. 6	1.2201	0.3157	1.3335	0.3441	3.5717	0.9036
	No. 7	1.2184	0.3156	1.3318	0.344	3.57	0.9035
	No. 8	1.2185	0.3156	1.3318	0.344	3.57	0.9035
	No. 9	1.2185	0.3156	1.3318	0.344	3.57	0.9035
	No. 10	1.2186	0.3156	1.332	0.344	3.5702	0.9035
	No. 11	1.2203	0.3158	1.3337	0.3441	3.5719	0.9036
	No. 12	1.2184	0.3156	1.3318	0.344	3.5699	0.9035
	No. 13	1.2186	0.3156	1.3319	0.344	3.5701	0.9035

, critical buckling temperatures (∆t_cr = ∆Tcr·10⁻³) of FGM plates placed on an elastic foundation for case a linear temperature change through the plate thickness are presented. Values such as the power law index p, Winkler and Pasternak coefficient k₀, k₁ and thickness ratios a/h were varied during the analysis of the influence on the critical buckling temperature. By drawing an analogy between these findings where 13 different shape functions were applied, it can be clearly seen that the newly established shape function demonstrates quite similar results. From the findings, we can infer that critical buckling temperatures drops off with the rise of power law index p and ratio a/h. Additionally, Pasternak coefficient k₁ has a bigger effect than Winkler coefficient k₀ on the critical buckling temperatures.

Table 4. Critical buckling temperatures (∆t_cr) of FGM plates placed on an elastic foundation for case of nonlinear temperature change across plate thickness (a/b = 1, m = n = 1, s = 3, and T_m = 5°).

p	Source	∆tcr
		k0 = 0, k1 = 0		k0 = 10, k1 = 0		k0 = 10, k1 = 10
		a/h = 10	a/h = 20	a/h = 10	a/h = 20	a/h = 10	a/h = 20
0	[55]	6.4552	1.6661	6.6308	1.71	10.0958	2.5763
	Present study	6.4547	1.6661	6.6303	1.71	10.0953	2.5763
	No. 1	6.4547	1.6661	6.6302	1.71	10.0953	2.5762
	No. 2	6.4547	1.6661	6.6302	1.71	10.0953	2.5762
	No. 3	6.4552	1.6661	6.6308	1.71	10.0958	2.5763
	No. 4	6.4667	1.6669	6.6422	1.7108	10.1073	2.577
	No. 5	6.4547	1.6661	6.6302	1.7182	10.0953	2.5762
	No. 6	6.4564	1.6662	6.632	1.71	10.097	2.5764
	No. 7	6.4569	1.6662	6.6324	1.7101	10.0974	2.5764
	No. 8	6.4571	1.6663	6.6326	1.7101	10.0977	2.5764
	No. 9	6.4571	1.6663	6.6326	1.7101	10.0977	2.5764
	No. 10	6.4547	1.6661	6.6302	1.7101	10.0953	2.5762
	No. 11	6.4577	1.6663	6.6332	1.7139	10.0983	2.5764
	No. 12	6.455	1.6661	6.6305	1.7102	10.0956	2.5762
	No. 13	6.4547	1.6661	6.6302	1.71	10.0953	2.5762
1	[55]	2.8269	0.7176	2.9804	0.756	6.0097	1.5133
	Present study	2.8268	0.7176	2.9802	0.756	6.0096	1.5133
	No. 1	2.8267	0.7176	2.9802	0.7559	6.0095	1.5133
	No. 2	2.8267	0.7176	2.9802	0.7559	6.0095	1.5133
	No. 3	2.8269	0.7176	2.9804	0.756	6.0097	1.5133
	No. 4	2.8312	0.7179	2.9847	0.7562	6.014	1.5136
	No. 5	2.8267	0.7176	2.9802	0.7559	6.0095	1.5133
	No. 6	2.8274	0.7176	2.9808	0.756	6.0102	1.5133
	No. 7	2.8275	0.7176	2.981	0.756	6.0103	1.5133
	No. 8	2.8276	0.7176	2.9811	0.756	6.0104	1.5133
	No. 9	2.8276	0.7176	2.9811	0.756	6.0104	1.5133
	No. 10	2.8267	0.7176	2.9802	0.7559	6.0095	1.5133
	No. 11	2.8278	0.7177	2.9813	0.756	6.0106	1.5134
	No. 12	2.8268	0.7176	2.9803	0.756	6.0096	1.5133
	No. 13	2.8267	0.7176	2.9802	0.7559	6.0095	1.5133
5	[55]	2.0152	0.5188	2.1847	0.5612	5.5309	1.3977
	Present study	2.0165	0.5189	2.186	0.5613	5.5322	1.3978
	No. 1	2.0162	0.5188	2.1858	0.5612	5.5319	1.3978
	No. 2	2.0162	0.5188	2.1858	0.5612	5.5319	1.3978
	No. 3	2.0152	0.5188	2.1847	0.5611	5.5309	1.3977
	No. 4	2.0157	0.5188	2.1853	0.5612	5.5315	1.3977
	No. 5	2.0163	0.5188	2.1858	0.5612	5.532	1.3978
	No. 6	2.0197	0.5191	2.1892	0.5615	5.5354	1.398
	No. 7	2.0146	0.5187	2.1841	0.5611	5.5303	1.3977
	No. 8	2.0145	0.5187	2.1841	0.5611	5.5302	1.3977
	No. 9	2.0145	0.5187	2.1841	0.5611	5.5302	1.3977
	No. 10	2.0164	0.5189	2.1859	0.5612	5.5321	1.3978
	No. 11	2.0206	0.5191	2.1901	0.5615	5.5363	1.3981
	No. 12	2.0154	0.5188	2.1849	0.5612	5.5311	1.3977
	No. 13	2.0161	0.5188	2.1856	0.5612	5.5318	1.3978
10	[55]	2.0971	0.5433	2.2923	0.5921	6.1448	1.5552
	Present study	2.0978	0.5434	2.2929	0.5922	6.1455	1.5553
	No. 1	2.0976	0.5433	2.2928	0.5921	6.1453	1.5553
	No. 2	2.0976	0.5433	2.2928	0.5921	6.1453	1.5553
	No. 3	2.0971	0.5433	2.2923	0.5921	6.1448	1.5552
	No. 4	2.101	0.5436	2.2961	0.5924	6.1487	1.5555
	No. 5	2.0976	0.5433	2.2928	0.5921	6.1453	1.5553
	No. 6	2.1002	0.5435	2.2954	0.5923	6.1479	1.5554
	No. 7	2.0973	0.5433	2.2925	0.5921	6.145	1.5552
	No. 8	2.0973	0.5433	2.2925	0.5921	6.145	1.5552
	No. 9	2.0973	0.5433	2.2925	0.5921	6.145	1.5552
	No. 10	2.0977	0.5433	2.2928	0.5921	6.1454	1.5553
	No. 11	2.1005	0.5435	2.2957	0.5923	6.1482	1.5555
	No. 12	2.0972	0.5433	2.2924	0.5921	6.1449	1.5552
	No. 13	2.0975	0.5433	2.2927	0.5921	6.1452	1.5553

illustrates comparative findings of critical buckling temperatures of FGM plates placed on an elastic foundation for case a nonlinear temperature change through the plate thickness. By drawing an analogy between these findings, the ones where 13 different shape functions were applied, and the ones included in the sources [55], it can be clearly seen that the newly established shape function demonstrates quite similar results.

Figure 4a and Figure 5a clarify the impact of ratio a/b on the critical buckling temperature for case of linear and nonlinear temperature change through the plate thickness. It should be stressed that the rise of a/b ratio leads to the rise of critical temperature, e.g., curves show the harsh rise of a/b ratio functions. Furthermore, the critical buckling temperatures in the same instance are higher for nonlinear than the linear temperature change through the plate thickness.

Figure 4b and Figure 5b illustrate the impact of ratio a/h on the critical buckling temperature for case of linear and nonlinear temperature change through the plate thickness. This diagram exhibits the values of p = 0, 1, 5, 10 and it is evident that the most prominent curve is the one whose value is p = 0, while the remaining curves lie over each other when the value of ratio is a/h > 15.

The impact of the elastic foundation on the critical buckling temperature for case of linear and nonlinear temperature change through the plate thickness is shown in Figure 6a,b. It should be underlined that results are presented for the shape function is No.10 because of its analytical integration. Numerical integration does not have to be employed.

The numerical calculations were done to figure out the influence of every elastic foundation parameter (k₀ and k₁). This was achieved by varying one of the parameters and setting the other one as a constant. Similarly to previous cases, the analysis is conducted for a linear and nonlinear rise in temperature and the findings are illustrated in Figure 7a and Figure 8a (fixed value k₁ and variation of k₀) and also in Figure 7b and Figure 8b (fixed value k₀ and variation of k₁). The Figures provide information about the curves that represent critical buckling temperature changes and their rapid rise that occurs because of the change of coefficient k₁, rather than when coefficient k₀ changes. Another important observation is that the curve with the most rapid rise has the value p = 10.

5.2. Free Vibration Analysis

In order to properly evaluate the behaviour of FGM plates, in addition to the results of the static analysis, an analysis of the behaviour of plates in a dynamic environment is also required. This section presents the results of free vibrations of FGM plates placed on an elastic foundation for different values of Winkler coefficient k₀, Pasternak coefficient k₁, and power law index p. As with the results of the thermal analysis described above, the verification of the developed and implemented theoretical results based on the newly introduced shape function was performed through a tabular representation of the obtained results, in comparison with the results from the literature. The procedure verified in this way was used to obtain other results for plates of different gradient structure. Based on the results of the dynamic analysis, appropriate interpretations and the comments were provided, and certain conclusions made. In order to display the numerical values of natural frequencies, for rectangular and square FGM plates, it is necessary to normalize the obtained values according to:

$$\tilde{\mathsf{\omega }}=\mathsf{\omega }\mathrm{h}\sqrt{\frac{{\mathsf{\rho }}_{\mathrm{m}}}{{\mathrm{E}}_{\mathrm{m}}}}.$$

(34)

Table 5. Non-dimensional natural frequencies $(\tilde{\mathsf{\omega }})$ of rectangular FGM plates placed on an elastic foundation for different values of the Winkler coefficient k₀, Pasternak coefficient k₁, and power law index p (a/b = 0.5, a/h = 5, m = 1, n = 1).

a/b	k0	k1	Theory	$\tilde{\mathsf{\omega }}$
				a/h = 5
				p = 0	p = 1	p = 5	p = 10
0.5	0	0	Present study	6.761	5.2016	4.3761	4.206
			No. 1	6.7609	5.2015	4.3757	4.2058
			No. 2	6.7609	5.2015	4.3757	4.2058
			No. 3	6.7616	5.202	4.3733	4.205
			No. 4	6.775	5.2108	4.3753	4.2136
			No. 5	6.7609	5.2015	4.3757	4.2058
			No. 6	6.7628	5.2027	4.3832	4.211
			No. 7	6.7636	5.2033	4.3722	4.2055
			No. 8	6.7638	5.2034	4.3721	4.2056
			No. 9	6.7638	5.2034	4.3721	4.2056
			No. 10	6.7609	5.2015	4.3759	4.2059
			No. 11	6.7642	5.2036	4.3852	4.2117
			No. 12	6.8031	5.2291	4.4434	4.2771
			No. 13	6.7609	5.2015	4.3754	4.2056
	100	0	Present study	7.2125	5.8654	5.2358	5.1214
			No. 1	7.2125	5.8653	5.2354	5.1211
			No. 2	7.2125	5.8653	5.2354	5.1211
			No. 3	7.2132	5.8657	5.2336	5.1205
			No. 4	7.2256	5.8734	5.2353	5.1274
			No. 5	7.2125	5.8653	5.2355	5.1212
			No. 6	7.2142	5.8664	5.2415	5.1252
			No. 7	7.215	5.8668	5.2327	5.1209
			No. 8	7.2152	5.867	5.2326	5.121
			No. 9	7.2152	5.867	5.2326	5.121
			No. 10	7.2125	5.8653	5.2357	5.1212
			No. 11	7.2155	5.8672	5.2431	5.1258
			No. 12	7.2517	5.8893	5.2902	5.1777
			No. 13	7.2125	5.8653	5.2352	5.121
	0	100	Present study	11.115	10.845	10.992	11.079
			No. 1	11.115	10.845	10.9919	11.0793
			No. 2	11.115	10.845	10.9919	11.0793
			No. 3	11.1154	10.8452	10.9914	11.0791
			No. 4	11.1226	10.8484	10.9922	11.0814
			No. 5	11.115	10.845	10.9919	11.0793
			No. 6	11.116	10.8455	10.9936	11.0803
			No. 7	11.1164	10.8457	10.9912	11.0793
			No. 8	11.1166	10.8457	10.9912	11.0794
			No. 9	11.1166	10.8457	10.9912	11.0794
			No. 10	11.115	10.845	10.992	11.0793
			No. 11	11.1168	10.8458	10.994	11.0804
			No. 12	11.138	10.8552	11.0077	11.0949
			No. 13	11.115	10.845	10.9918	11.0792
	100	100	Present study	11.395	11.178	11.3593	11.4558
			No. 1	11.3952	11.178	11.3593	11.4558
			No. 2	11.3952	11.178	11.3593	11.4558
			No. 3	11.3956	11.1782	11.3588	11.4557
			No. 4	11.4026	11.1812	11.3596	11.4578
			No. 5	11.3952	11.178	11.3593	11.4558
			No. 6	11.3962	11.1784	11.3608	11.4567
			No. 7	11.3966	11.1786	11.3587	11.4559
			No. 8	11.3967	11.1787	11.3587	11.4559
			No. 9	11.3967	11.1787	11.3587	11.4559
			No. 10	11.3952	11.178	11.3593	11.4558
			No. 11	11.3969	11.1787	11.3612	11.4568
			No. 12	11.4174	11.1876	11.3737	11.47
			No. 13	11.3952	11.178	11.3592	11.4557

and

Table 6. Non-dimensional natural frequencies $(\tilde{\mathsf{\omega }})$ of square FGM plates placed on an elastic foundation for different values of the Winkler coefficient k₀, Pasternak coefficient k₁, and power law index p (a/b = 1, a/h = 5, m = 1, n = 1).

a/b	k0	k1	Theory	$\tilde{\mathsf{\omega }}$
				a/h = 5
				p = 0	p = 1	p = 5	p = 0
1	0	0	Present study	10.3761	8.0121	6.6687	6.3883
			No. 1	10.3761	8.0121	6.6677	6.3879
			No. 2	10.3761	8.0121	6.6677	6.3879
			No. 3	10.3779	8.0133	6.663	6.3864
			No. 4	10.4086	8.0336	6.6684	6.4062
			No. 5	10.3761	8.0121	6.6679	6.3879
			No. 6	10.38	8.0147	6.6838	6.3987
			No. 7	10.3824	8.0163	6.6607	6.3877
			No. 8	10.3831	8.0167	6.6606	6.388
			No. 9	10.3831	8.0167	6.6606	6.388
			No. 10	10.3761	8.0121	6.6683	6.3881
			No. 11	10.383	8.0166	6.6881	6.4001
			No. 12	10.4698	8.0739	6.8155	6.5423
			No. 13	10.3762	8.0122	6.6672	6.3876
	100	0	Present study	10.6722	8.4517	7.2542	7.0179
			No. 1	10.6723	8.4517	7.2534	7.0175
			No. 2	10.6723	8.4517	7.2534	7.0175
			No. 3	10.674	8.4528	7.2491	7.0162
			No. 4	10.7037	8.4718	7.2541	7.0339
			No. 5	10.6722	8.4517	7.2535	7.0176
			No. 6	10.676	8.4541	7.2678	7.0271
			No. 7	10.6783	8.4556	7.2471	7.0174
			No. 8	10.679	8.456	7.247	7.0177
			No. 9	10.679	8.456	7.247	7.0177
			No. 10	10.6722	8.4517	7.2539	7.0177
			No. 11	10.6789	8.4559	7.2717	7.0284
			No. 12	10.7629	8.5096	7.3865	7.1553
			No. 13	10.6723	8.4517	7.2529	7.0173
	0	100	Present study	15.1867	14.3818	14.3054	14.376
			No. 1	15.1867	14.3818	14.3052	14.3759
			No. 2	15.1867	14.3818	14.3052	14.3759
			No. 3	15.1878	14.3823	14.304	14.3757
			No. 4	15.2066	14.391	14.3064	14.3818
			No. 5	15.1867	14.3818	14.3052	14.376
			No. 6	15.1891	14.3829	14.3094	14.3785
			No. 7	15.1906	14.3836	14.3036	14.3763
			No. 8	15.191	14.3838	14.3036	14.3764
			No. 9	15.191	14.3838	14.3036	14.3764
			No. 10	15.1867	14.3818	14.3053	14.376
			No. 11	15.191	14.3838	14.3105	14.3788
			No. 12	15.2444	14.4086	14.3463	14.4167
			No. 13	15.1868	14.3818	14.305	14.3759
	100	100	Present study	15.3904	14.6304	14.5846	14.6636
			No. 1	15.3904	14.6305	14.5843	14.6636
			No. 2	15.3904	14.6305	14.5843	14.6636
			No. 3	15.3914	14.6309	14.5833	14.6634
			No. 4	15.4099	14.6394	14.5856	14.6692
			No. 5	15.3904	14.6305	14.5844	14.6636
			No. 6	15.3927	14.6315	14.5883	14.666
			No. 7	15.3941	14.6322	14.5829	14.6639
			No. 8	15.3945	14.6323	14.5829	14.664
			No. 9	15.3945	14.6323	14.5829	14.664
			No. 10	15.3904	14.6305	14.5845	14.6636
			No. 11	15.3945	14.6324	14.5894	14.6663
			No. 12	15.447	14.6564	14.6233	14.7021
			No. 13	15.3904	14.6305	14.5842	14.6635

show the non-dimensional natural frequencies $\tilde{\mathsf{\omega }}$ of rectangular (a/b = 0.5) and square (a/b = 1) plates resting on an elastic foundation for different values of Winkler coefficient (k₀), Pasternak coefficient (k₁) and index p. In order to observe the impact of the elastic foundation, the values of k₀ = 0 and k₁ = 0 were first taken, introducing one coefficient at the time in order to determine which of the two coefficients has a greater impact. By analyzing the results, it is evident that the impact of the coefficient k₁ on the $\tilde{\mathsf{\omega }}$ is far greater than the impact of the coefficient k₀.

Figure 9 shows the diagrams of the non-dimensional natural frequencies $\tilde{\mathsf{\omega }}$ of the plates resting on an elastic foundation for different ratios a/h, a/b, and the values of the coefficients k₀, k₁, and index p. Figure 9a shows the impact of the Winkler coefficient (k₀) and Pasternak coefficient (k₁) on the values of natural frequencies for the first oscillation mode. It can be clearly observed that the introduction of the coefficient k₀ leads to very small changes in the $\tilde{\mathsf{\omega }}$ values in comparison to the case of absence of an elastic foundation (k₀ = 0 and k₁ = 0). On the other hand, when you introduce the coefficient k₁, there is an evident change in the value of natural frequencies $\tilde{\mathsf{\omega }}$. The effect of the geometry changes of the plate (a/b ratio) and the index p is shown in the diagram 9b. With the increase of a/b ratio the curves become further away from each other, i.e., the fastest change in the $\tilde{\mathsf{\omega }}$ value occurs in a ceramic plate, the slowest change occurs in a metal plate, and in FGM the rate of change depends on the ratio of constituents. Figure 9c shows the effect of different a/h ratios on the values $\tilde{\mathsf{\omega }}.$ For smaller values of a/h ratios, the changes in values $\tilde{\mathsf{\omega }}$ are greater, while at values of a/h > 5, the effect on the values $\tilde{\mathsf{\omega }}$ is decreasing.

Figure 10 shows 3D diagrams of the non-dimensional natural frequencies $\tilde{\mathsf{\omega }}$ of the plates resting on an elastic foundation for different ratios a/h, a/b, and values k₀, k₁, and the index p. This visualization provides a more transparent insight into the previously described effects of certain parameters on the natural frequency values and the conclusions reached. For example, at Figure 10c it can be clearly seen that the impact of the coefficient k₁ on the frequency values is far greater than the impact of the coefficient k₀.

6. Conclusions

The obtained results presented in the already published articles have been the foundation for developing and introducing the new shape function. The results obtained put an emphasis on the significance and topicality of the research in the field of functionally graded materials. A comprehensive and detailed investigation and systematization of the literature according to the topic have been guided related to the type of problems that authors tried to solve during the analysis of the functionally graded materials. Special focus and attention has been paid to different shear deformation theories that authors have used during the research. The new introduced shape function has been compared to 13 other shape functions that were originally proposed by different authors for the purpose of analysing composite laminates. However, this article implemented the previously mentioned shape 13 shape functions as well as new proposed shape function in order to analyse FGM plates. By comparing obtained results related to the static and dynamic analysis of moderately thick and thick plates, it is possible to conclude that the newly presented shape function can be applied during the analysis of FGM plates. Generally, based on the above research as well as obtained results, the following conclusions could be emphasized:

• Decreasing the volume fraction of ceramics and increasing the volume fraction of metal in the FGM (the value of p index increases) decreases the value of the critical buckling temperature for both linear and nonlinear cases of temperature distribution through plate thickness
• Comparative analysis of the results for the linear and nonlinear distribution of the temperature across the plate thickness, and for other fixed parameters of the plate, it can be concluded that higher critical buckling temperatures are obtained for nonlinear distribution
• The elastic foundation effect shows that critical buckling temperature rapid rise because of the change of Pasternak coefficient k₁, rather than when Winkler coefficient k₀ changes
• Based on the analysis of the impact of the Winkler-Pasternak elastic foundation model parameters, similar to the thermal analysis, it was pointed out that the Pasternak coefficient k₁ has a far greater influence on natural frequencies than the Winkler coefficient k₀