New Shape Function for the Bending Analysis of Functionally Graded Plate

The bending analysis of thick and moderately thick functionally graded square and rectangular plates as well as plates on Winkler–Pasternak elastic foundation subjected to sinusoidal transverse load is presented in this paper. The plates are assumed to have isotropic, two-constituent material distribution through the thickness, and the modulus of elasticity of the plate is assumed to vary according to a power-law distribution in terms of the volume fractions of the constituents. This paper presents the methodology of the application of the high order shear deformation theory based on the shape functions. A new shape function has been developed and the obtained results are compared to the results obtained with 13 different shape functions presented in the literature. Also, the validity and accuracy of the developed theory was verified by comparing those results with the results obtained using the third order shear deformation theory and 3D theories. In order to determine the procedure for the analysis and the prediction of behavior of functionally graded plates, the new program code in the software package MATLAB has been developed based on the theories studied in this paper. The effects of transversal shear deformation, side-to-thickness ratio, and volume fraction distributions are studied and appropriate conclusions are given.


Introduction
Failure and delamination at the border between two layers are the biggest and the most frequently studied problem of the conventional composite laminates. Delamination of layers due to high local inter-laminar stresses causes a reduction of stiffness and a loss of structural integrity of a construction. In order to eliminate these problems, improved materials such as functionally graded materials (FGM), which are getting more and more popular, are used for innovative engineering constructions.
FGM is a composite material consisting of two or more constituents with the continuous change of properties in a certain direction. In other words, these materials can also be defined as materials which possess a gradient change of properties due to material heterogeneity. A gradient property can go in one or more directions and it can also be continuous or discontinuous from one surface to another depending on the production technique [1][2][3]. One of the most common uses of FGM materials is found in thermal barriers, one surface of which is in contact with high temperatures and is made of ceramic which can provide adequate thermal stability, low thermal conductivity, and fine antioxidant properties. The low-temperature side of the barrier is made of metal, which is superior in terms of mechanical in terms of mechanical strength, toughness, and high thermal conductivity. Functionally graded materials, which contain metal and ceramic constituents, improve thermo-mechanical properties between layers, because of which delamination of layers should be avoided due to continuous change between properties of the constituents. By varying the percentage of volume fraction content of the two or more materials, FGM can be formed so that it achieves a desired gradient property in specific directions. Figure 1 shows schematic of continuously graded microstructure with metal-ceramic constituents [4]. Depending on the nature of gradient, functionally graded materials may be grouped into fraction gradient type, shape gradient type, orientation gradient type and size gradient type ( Figure  2) [5]. With the expansion of the FGM material application area, it was necessary to improve fabrication methods for mentioned materials. Various fabrication methods have been developed for the preparation of bulk FGMs and graded thin films. The processing methods are commonly classified into four groups like powder technology methods (dry powder processing, slip vesting, tape casting, infiltration process or electrochemical gradation, powder injection molding and self-propagating high temperature synthesis, etc.), deposition methods (chemical vapor deposition, physical vapor deposition, electrophoretic deposition, slurry deposition, pulsed laser deposition, plasma spraying, etc.), in-situ processing methods (laser cladding, spray forming, sedimentation and solidification, centrifugal casting, etc.), and rapid prototyping processes (multiphase jet solidification, 3D-printing, laser printing, laser sintering, etc.) [6]. The basic difference between the mentioned production methods can be made according to whether the obtained materials have a stepwise or continuous Depending on the nature of gradient, functionally graded materials may be grouped into fraction gradient type, shape gradient type, orientation gradient type and size gradient type ( Figure 2) [5].
Materials 2018, 11, x FOR PEER REVIEW 2 of 28 in terms of mechanical strength, toughness, and high thermal conductivity. Functionally graded materials, which contain metal and ceramic constituents, improve thermo-mechanical properties between layers, because of which delamination of layers should be avoided due to continuous change between properties of the constituents. By varying the percentage of volume fraction content of the two or more materials, FGM can be formed so that it achieves a desired gradient property in specific directions. Figure 1 shows schematic of continuously graded microstructure with metal-ceramic constituents [4]. Depending on the nature of gradient, functionally graded materials may be grouped into fraction gradient type, shape gradient type, orientation gradient type and size gradient type ( Figure  2) [5]. With the expansion of the FGM material application area, it was necessary to improve fabrication methods for mentioned materials. Various fabrication methods have been developed for the preparation of bulk FGMs and graded thin films. The processing methods are commonly classified into four groups like powder technology methods (dry powder processing, slip vesting, tape casting, infiltration process or electrochemical gradation, powder injection molding and self-propagating high temperature synthesis, etc.), deposition methods (chemical vapor deposition, physical vapor deposition, electrophoretic deposition, slurry deposition, pulsed laser deposition, plasma spraying, etc.), in-situ processing methods (laser cladding, spray forming, sedimentation and solidification, centrifugal casting, etc.), and rapid prototyping processes (multiphase jet solidification, 3D-printing, laser printing, laser sintering, etc.) [6]. The basic difference between the mentioned production methods can be made according to whether the obtained materials have a stepwise or continuous On the other hand, continuum-based 3D elasticity theory could be used for the analysis of these plates. However, 3D solution methods are mathematically complex which consequently results in prolonged calculation time and the need for high performance hardware. Taking the aforementioned into consideration, developing and using 2D shear deformation plate theories, which consider the effects of previously mentioned shear and normal strains and provide the precision in the same way as 3D models do, represents a trend in the process of analysis of FGM plates. This paper presents, in detail, the methodology of the application of the HSDT theory based on the shape functions. A new shape function has been developed and the obtained results are compared to the results obtained with 13 different shape functions presented in the papers from the reference list. Also, the results have been verified through comparison with the results obtained with TSDT and 3D theories. In order to determine a procedure for the analysis and the prediction of behavior of FGM plates, the new program code in the software package MATLAB (MATrix LABoratory) has been developed based on theories studied in this paper.
Finally, the ultimate goal and the purpose of all the previously mentioned studies and analyses is the application of FGM in different areas of engineering and branches of industry. Although FGM were initially used as materials for thermal barrier in space shuttles, today they are becoming widely used in the field of medicine, dentistry, energy and nuclear sector, automotive industry, military, optoelectronics etc.

Description of the Problem
The subject of the analysis in this paper are FGM plate ( Figure 3a) and FGM plate on elastic foundation (Figure 3b). The plate (length a, width b and height h) is made of functionally graded material consisting of the two constituents, namely, metal and ceramics. On the other hand, continuum-based 3D elasticity theory could be used for the analysis of these plates. However, 3D solution methods are mathematically complex which consequently results in prolonged calculation time and the need for high performance hardware. Taking the aforementioned into consideration, developing and using 2D shear deformation plate theories, which consider the effects of previously mentioned shear and normal strains and provide the precision in the same way as 3D models do, represents a trend in the process of analysis of FGM plates.
This paper presents, in detail, the methodology of the application of the HSDT theory based on the shape functions. A new shape function has been developed and the obtained results are compared to the results obtained with 13 different shape functions presented in the papers from the reference list. Also, the results have been verified through comparison with the results obtained with TSDT and 3D theories. In order to determine a procedure for the analysis and the prediction of behavior of FGM plates, the new program code in the software package МATLAB (MATrix LABoratory) has been developed based on theories studied in this paper.
Finally, the ultimate goal and the purpose of all the previously mentioned studies and analyses is the application of FGM in different areas of engineering and branches of industry. Although FGM were initially used as materials for thermal barrier in space shuttles, today they are becoming widely used in the field of medicine, dentistry, energy and nuclear sector, automotive industry, military, optoelectronics etc.

Description of the Problem
The subject of the analysis in this paper are FGM plate ( Figure 3a) and FGM plate on elastic foundation (Figure 3b). The plate (length a, width b and height h) is made of functionally graded material consisting of the two constituents, namely, metal and ceramics. It is assumed that mechanical properties of the FGM in the thickness direction of the plate change according to the power law distribution (Figure 4a): This law defines the change of the mechanical properties as the function of the volume fraction of the FGM constituents in the thickness direction of the plate. It is assumed that mechanical properties of the FGM in the thickness direction of the plate change according to the power law distribution ( Figure 4a): In the Equation (1), h represents total thickness of the plate, and P(z) represents a material property in an arbitrary cross-section z, −h/2 < z < h/2. Pc represents the material property at the top of the plate z = h/2 − ceramic, and Pm represents the material property at the bottom of the plate z = −h/2 − metal. Index p is the exponent of the equation which defines the volume fraction of the constituents in FGM. Practically, by varying the index p, homogenous plates as well as FGM plate with precisely determined gradient structure could be obtained, as it is presented in Figure 4b: when p = 0 the plate is homogenous, made of ceramics, • when 0 < p < ∞ the plate has a gradient structure, • theoretically, when p = ∞ the plate becomes homogenous again, made of metal, although the plate can be considered homogenous even when p > 20.

Kinematic Displacement-Strain Relations and Constitutive Equation of Elasticity for FGM
According to HSDT based on the shape functions, displacements could be presented in the following way: In the reference literature there are many shape functions which can be polynomial, trigonometric, exponential, hyperbolic. Some examples of the shape functions are given in Table 1. This law defines the change of the mechanical properties as the function of the volume fraction of the FGM constituents in the thickness direction of the plate.
In the Equation (1), h represents total thickness of the plate, and P(z) represents a material property in an arbitrary cross-section z, −h/2 < z < h/2. P c represents the material property at the top of the plate z = h/2 − ceramic, and P m represents the material property at the bottom of the plate z = −h/2 − metal. Index p is the exponent of the equation which defines the volume fraction of the constituents in FGM. Practically, by varying the index p, homogenous plates as well as FGM plate with precisely determined gradient structure could be obtained, as it is presented in Figure 4b: when p = 0 the plate is homogenous, made of ceramics, • when 0 < p < ∞ the plate has a gradient structure, • theoretically, when p = ∞ the plate becomes homogenous again, made of metal, although the plate can be considered homogenous even when p > 20.

Kinematic Displacement-Strain Relations and Constitutive Equation of Elasticity for FGM
According to HSDT based on the shape functions, displacements could be presented in the following way: u(x, y, z, t) = u 0 (x, y, t) − z ∂w 0 (x,y,t) where: u 0 , v 0 , w 0 are displacement components in the middle plane of the plate, ∂w 0 ∂x , ∂w 0 ∂y are rotation angles of transverse normal in relation to x and y axes, respectively, θ x , θ y are rotations of the transverse normal due to transverse shear and f (z) is the shape function.
In the reference literature there are many shape functions which can be polynomial, trigonometric, exponential, hyperbolic. Some examples of the shape functions are given in Table 1.  [54] (h/π) sin(πz/h) SF 4 Mantari et al. [55] sin(πz/h)e cos (πz/h)/2 + (πz/2h) SF 5-6 Mantari et al. [45] tan(mz) − zm sec 2 (mh/2), m = {1/5h, π/2h} SF 7 Karama et al. [56], Aydogdu [44] z exp −2(z/h) 2 , z exp −2(z/h) 2 / ln α , ∀α> 0 SF 8 Mantari et al. [46] z · 2.85 −2(z/h) 2 + 0.028z SF 9 El Meiche et al. [47] ξ[(h/π) sin(πz/h) − z], ξ = {1, 1/ cosh(π/2) − 1} SF 10 Soldatos [43] hsinh(z/h) − z cosh(1/2) SF 11 Akavci and Tanrikulu [49] z sec h z 2 /h 2 − z sec h(π/4)[1 − (π/2)tanh(π/4)] SF 12 Akavci and Tanrikulu [49] (3π/2)htanh(z/h) − (3π/2)z sec h 2 (1/2) SF 13 Mechab et al. [48] z cos(1/2) This paper proposes a new shape function as follows: The introduced shape function is an odd function of the thickness coordinate z and satisfies zero stress conditions for out of plane shear stresses. Observing the shape functions in the Table 1, may see that the proposed function belongs to the group of simple mathematic functions. This fact makes the integration process easier and thus reduces considerably the calculation time. Having in mind that the function is analytically integrable, there is no need to switch to numeric integration, which additionally increases the precision of the obtained results. The verification of the above claims is shown in the comparative diagrams ( Figure 5) of the newly introduced shape function and the shape functions given in the Table 1. These shape functions' diagrams can be categorized into two groups of functions. In both cases it can be seen in the diagram that, in the case of the ratio z/h = 0.5, all shape functions have extreme values, which are different (Figure 5a For small displacements and moderate rotations of a transverse normal in relation to x axis and y axis, normal and shear strain components are obtained by well-known relations in linear elasticity between displacements and strains: where: For small displacements and moderate rotations of a transverse normal in relation to x axis and y axis, normal and shear strain components are obtained by well-known relations in linear elasticity between displacements and strains: where: where f (z) = d f (z) dz is the first derivative of the shape function in the thickness direction of the plate. The elastic constitutive relations for FGM are given as follows: where the coefficients of the constitutive elasticity tensor could be defined through engineering constants: Due to the gradient change of the plate structure in the direction of the z coordinate, based on (1), the modulus of elasticity could be defined as: while Poisson's ratio ν is considered constant due to a small value variation in the thickness direction of the plate, ν = const.
As it could be seen, the coefficients of the constitutive tensor are functionally dependent on the z coordinate which practically means that for p = 0 there is a finite number of planes parallel to the middle plane, where each of these planes has different values of the constitutive tensor C ij .

Bending of FGM Plates and FGM Plates on Elastic Foundation
It is assumed that the plate is loaded with an arbitrary transverse load q(x, y). Work under external load is defined as: where: is the sinusoidal transverse load with an amplitude q 0 . Plate strain energy is defined as: where force, moments and higher order moments vectors are obtained in the following form: Matrices in the developed form could be presented in the following way: In the Equation (12) by grouping the terms with the elements of constitutive tensor, new matrices with the following components could be defined: Therefore, load vectors could now be defined in the following form: By exchanging plate strain energy (11) and work under external load (9) into the equation which defines the minimum total potential energy principle: The following form is obtained: Materials 2018, 11, 2381 9 of 24 By exchanging the strain components (5) and by applying the calculus of variations, the following equilibrium equations are obtained: δu 0 : N xx,x + N xy,y = 0, δv 0 : N yy,y + N xy,x = 0, δw 0 : M xx,xx + 2M xy,xy + M yy,yy + q = 0, δθ x : P xx,x + P xy,y − R x = 0, δθ y : P xy,x + P yy,y − R y = 0. (18) which could be further solved through analytical and numerical methods.
In the case of a plate on elastic foundation, in the Equation (16) deformation energy of the elastic foundation should be taken into consideration, which is defined using Winkler-Pasternak model in the following way: Using the previously mentioned the minimum total potential energy principle, the equilibrium equations of the plate on elastic foundation are the following:

Analytical Solution of the Equilibrium Equations
Although analytical solution methods are limited to simple geometrical problems, boundary conditions and loads, they can provide a clear understanding of the physical aspect of the problem and its solutions are very precise. Since analytical solutions are extremely important for developing new theoretical models, primarily due to their understanding of the physical aspects of the problem, and considering that a new HSDT theory based on a new shape function has been developed in this paper, the analytical solution of the equilibrium equations for a rectangular plate will be presented in the following part of the paper. For complex engineering calculations, which include solving the system of a large number of equations, it is necessary to use numerical methods which provide approximate, but satisfactory results.
For a simply supported rectangular FGM plate, boundary conditions are defined based on [57] as: In order to satisfy these kinematic boundary conditions, assumed forms of Navier's solutions are introduced: Materials 2018, 11, 2381 10 of 24 The equilibrium equation is further developed into: or: Through the matrix multiplication of the Equation (24) with L −1 , the following is obtained: The Equation (25) fully defines the amplitudes of the assumed displacement components. The displacement components are obtained if the displacement amplitude matrix is multiplied with the vector from trigonometric functions which depend on x and y.

Numerical Results
In order to apply the previously obtained theoretical results to a simulation of real problems, a new program code for static analysis of FGM plates has been developed within the software package MATLAB. Material properties of the used materials are shown in Table 2 [58]. Table 2. Material properties of FGM constituents.

Elasticity Modulus, E[GPa]
Poisson's Ratio, ν Normalized values of a vertical displacement w (deflection), normal stresses σ xx and σ yy , shear stress τ xy , and transverse shear stresses τ xz and τ yz are given by using HSDT theory based on the new shape function. Normalization of the aforementioned values has been conducted according to (26) as: Table 3 shows comparative results of the normalized values of displacement and stresses of square plate for two different ratios of length and thickness of the plate (a/h = 5 and a/h = 10) and for different values of the index p. Verification of the results obtained in this paper has been conducted by comparing them to the results from the reference papers when a/h = 10. Based on that, the results when a/h = 5 are provided for different values of the index p, i.e., different volume fraction of the constituents in FGM. Using HSDT theory with the new shape function, the obtained results are compared to the results obtained using 13 different shape functions as well as to the results obtained using quasi 3D theory of elasticity [59] and TSDT theory [58]. The results based on the CPT theory are also presented [60] in order to find certain disadvantages of the theory. Based on the comparative results of displacement and stresses, which are provided in this paper and in previously mentioned theories, it could be seen that there is a match with both TSDT theory and quasi 3D theory of elasticity. On the other hand, it is clearly seen that there are some significant differences in the results obtained by CPT theory, especially related to the stress σ xx which shows that CPT theory does not provide satisfying results in the analysis of thick and moderately thick FGM plates. A comparative review of these results with the results obtained using 13 different shape functions shows that the newly given shape function provides almost identical results. However, since these results are given for the plane on a certain height z (for example, stress σ xx on the height of h/3 etc.), a real insight into the values obtained by varying the new function could be offered by presenting stress distribution across the thickness of the plate, which is done through appropriate diagrams. Figure 6 shows the distribution of normal stresses σ xx and σ yy across the thickness of the plate for different values of the index p. By analyzing the diagrams, it could be noticed that the curves representing both stresses are identical. Also, the basic property of FGM could be noticed, namely, the shift of a neutral plane in relation to the plane z/h = 0. It can also be seen that for the planes at the height z/h = 0.    - [59] 0.  --  Figure 7 shows the distribution of the shear stress τ xy across the thickness of the plate for different values of the index p ( Figure 7a) and for different shape functions (Figure 7b), but for the unchanging values of a/h = 10 and a/b = 1. While analyzing the diagrams, it should be considered that when p = 0 the plate is homogenous made of ceramics, when p = 20 the plate is homogenous made of metal, and when 0 < p < 20 the plate is made of FGM. By analyzing the diagram in the Figure 7a, it could be noticed that for all values of the index p, the stress τ xy achieves the maximum value on the upper edge of the plate. Ceramic plate has the lowest maximum value. Therefore, with an increase of the metal volume fraction when p = 1, maximum stress value also increases and the highest value is achieved when the plate is homogenous made of metal. Moreover, apart from affecting the maximum stress values, the variation of the index p value also affects the shape of the τ xy stress distribution curve across the thickness of the plate. In order to conduct a comparative analysis of the results for different shape functions and to estimate the application of the new shape function to the given problems, Figure 7b shows the distribution of the shear stress τ xy by using newly developed shape function and the shape functions given in Table 1. It is clearly seen that all the previously mentioned shape functions give identical results to the results obtained with the new shape function. Figure 8 shows the distribution of transverse shear stresses τ xz and τ yz across the thickness of the plate for different values of the index p and for different shape functions. By analyzing transverse shear stresses in Figure 8a,c, a basic distinction between homogenous and FGM plates can be noticed. When plates are made of ceramics (p = 0) or metal (p = 20), it can be noticed that both stresses achieve maximum values in the plane at the height z/h = 0, due to the homogeneity of the material. On the other hand, when FGM plates are considered, there is an asymmetry in relation to the plane z/h = 0, therefore, when p = 1 stresses achieve maximum values in the plane z/h = 0.15, and when p = 5 stresses achieve maximum values when z/h = 0.3. In contrast to the homogenous ceramic plate, where stress distribution curve is a parabola with the maximum value in the plane z/h = 0, plates with the larger volume fraction of metal (p = 10) also achieve the maximum value of stress when z/h = 0, but the distribution curve is not a parabola. With the further increase of the metal volume fraction (p = 20), and although the plate can be practically considered homogenous, the diagram still shows the curve which is not a parabola. Generally, due to insignificant but still present ceramic fractions in the upper part of the plate, there is a slight deformity of the curve.
By conducting comparative analysis of the stresses τ xz and τ yz for different shape functions, and with fixed values of a/h = 10, a/b = 1 and p = 5, it could be seen in Figure 8b and Figure 8d   In order to conduct a comparative analysis of the results for different shape functions and to estimate the application of the new shape function to the given problems, Figure 7b shows the distribution of the shear stress τ xy by using newly developed shape function and the shape functions given in Table 1. It is clearly seen that all the previously mentioned shape functions give identical results to the results obtained with the new shape function. Figure 8 shows the distribution of transverse shear stresses τ xz and τ yz across the thickness of the plate for different values of the index p and for different shape functions. By analyzing transverse shear stresses in Figure 8a,c, a basic distinction between homogenous and FGM plates can be noticed. When plates are made of ceramics (p = 0) or metal (p = 20), it can be noticed that both stresses achieve maximum values in the plane at the height z/h = 0, due to the homogeneity of the material. On the other hand, when FGM plates are considered, there is an asymmetry in relation to the plane z/h = 0, therefore, when p = 1 stresses achieve maximum values in the plane z/h = 0.15, and when p = 5 stresses achieve maximum values when z/h = 0.3. In contrast to the homogenous ceramic plate, where stress distribution curve is a parabola with the maximum value in the plane z/h = 0, plates with the larger volume fraction of metal (p = 10) also achieve the maximum value of stress when z/h = 0, but the distribution curve is not a parabola. With the further increase of the metal volume fraction (p = 20), and although the plate can be practically considered homogenous, the diagram still shows the curve which is not a parabola. Generally, due to insignificant but still present ceramic fractions in the upper part of the plate, there is a slight deformity of the curve. In order to understand the effects of increasing the index p as well as the effect of the thickness and geometry, Figure 9 shows the diagram of the normalized values of the displacement w for different a/h and a/b ratios and values of the index p. By analyzing Figure 9a and Figure 9b, it could be noticed that the displacement values w for the metal plate (p = 20) are the highest, for the ceramic plate they are the lowest, and for the FGM plate they are somewhere in between. Moreover, by varying the volume fraction of metal or ceramics, a desired bending rigidity of the plate could be achieved. In Figure 9a, it could be seen that the curves gradually become closer when a/b > 4. In contrast to that, Figure 9b shows that with an increase of the ratio a/h the curves do not become closer, namely, the difference of the displacement ratio remains constant regardless of the index p change. This conclusion comes from the fact that in thin plates it is less possible to vary the volume fraction of the FGM constituents in the thickness direction of the plate and, thus, the index p has no effects. By conducting comparative analysis of the stresses τ xz and τ yz for different shape functions, and with fixed values of a/h = 10, a/b = 1 and p = 5, it could be seen in Figure 8b,d that, unlike the stress τ xy , the results do not match for all the shape functions. The most significant deviation could be noticed in the results for the El Meiche's and Karama's shape functions. The Akavci's function also shows a slight deviation and it achieves maximum stress value at the height z/h = 0.25, while the results for all the other shape functions are almost identical, achieving the maximum stress value in the plane z/h = 0.25.
In order to understand the effects of increasing the index p as well as the effect of the thickness and geometry, Figure 9 shows the diagram of the normalized values of the displacement w for different a/h and a/b ratios and values of the index p. By analyzing Figure 9a,b, it could be noticed that the displacement values w for the metal plate (p = 20) are the highest, for the ceramic plate they are the lowest, and for the FGM plate they are somewhere in between. Moreover, by varying the volume fraction of metal or ceramics, a desired bending rigidity of the plate could be achieved. In Figure 9a, it could be seen that the curves gradually become closer when a/b > 4. In contrast to that, Figure 9b shows that with an increase of the ratio a/h the curves do not become closer, namely, the difference of the displacement ratio remains constant regardless of the index p change. This conclusion comes from the fact that in thin plates it is less possible to vary the volume fraction of the FGM constituents in the thickness direction of the plate and, thus, the index p has no effects. In order to determine the effect of the elastic foundation on the displacements and stresses of the FGM plate, the results of different combinations of the FGM constituents have been presented, as well as different combinations of the Winkler (k0) and Pasternak (k1) coefficient of the elastic foundation. Apart from the normalization given in Error! Reference source not found., it is necessary to apply the normalization of the coefficients k0 and k1, in the following form: where the bending stiffness of the plate is ν The Tables 4 and 5 show the results of the normalized values of displacements and stresses of the square plate on elastic foundation for p = 5, and p = 10, different values of k0 and k1 coefficients, as well as for two different ratios length/thickness of the plate (a/h = 10 and a/h = 5). In order to determine the effect of the elastic foundation on the displacements and stresses of the plate, the values of displacements and stresses for k0 = 0 and k1 = 0 are first shown, which practically matches the case of the plate without the elastic foundation. Afterwards, the values of the given coefficients are varied in order to conclude which of the two has greater influence. Based on the results, it is concluded that the introduction of the coefficient k0 has less influence on the change of the displacements and stresses values then when only k1 coefficient is introduced. By introducing k0 and k1 coefficients, bending stiffness of the plate increases, i.e., displacement and stresses values decrease and the influence of the Winkler coefficient is smaller than the influence of the Pasternak coefficient. This phenomenon is especially noticeable in the diagram dependency which is to be shown later. In order to determine the effect of the elastic foundation on the displacements and stresses of the FGM plate, the results of different combinations of the FGM constituents have been presented, as well as different combinations of the Winkler (k 0 ) and Pasternak (k 1 ) coefficient of the elastic foundation. Apart from the normalization given in (26), it is necessary to apply the normalization of the coefficients k 0 and k 1 , in the following form: where the bending stiffness of the plate is D = E c h 3 12(1−ν 2 ) . The Tables 4 and 5 show the results of the normalized values of displacements and stresses of the square plate on elastic foundation for p = 5, and p = 10, different values of k 0 and k 1 coefficients, as well as for two different ratios length/thickness of the plate (a/h = 10 and a/h = 5). In order to determine the effect of the elastic foundation on the displacements and stresses of the plate, the values of displacements and stresses for k 0 = 0 and k 1 = 0 are first shown, which practically matches the case of the plate without the elastic foundation. Afterwards, the values of the given coefficients are varied in order to conclude which of the two has greater influence. Based on the results, it is concluded that the introduction of the coefficient k 0 has less influence on the change of the displacements and stresses values then when only k 1 coefficient is introduced. By introducing k 0 and k 1 coefficients, bending stiffness of the plate increases, i.e., displacement and stresses values decrease and the influence of the Winkler coefficient is smaller than the influence of the Pasternak coefficient. This phenomenon is especially noticeable in the diagram dependency which is to be shown later. Table 4. Normalized values of displacement and stresses of square plate on elastic foundation for p = 5, different values of the k 0 и k 1 and the ratio a/h (a/b = 1).    Figure 10 shows the effect of the Winkler coefficient k 0 on the distribution of the normal stress σ xx , shear stress τ xy and transversal shear stresses τ xz and τ yz across the thickness of the plate on the elastic foundation. By analyzing the diagram, it can be seen that the value of the stresses σ xx and τ xy equals zero for z/h = 0.15. On the other hand, the maximum values of τ xz and τ yz stresses are at z/h = 0.2 when the new proposed shape function is applied, while the maximum values of mentioned stresses is respectively at z/h = 0.15 i.e., z/h = 0.25 for Karama's shape function.    Figure 11 shows a comparative review of shear transversal stresses τ xz and τ yz distribution across the thickness of the plate on elastic foundation for different shape functions. As in the case of bending the plate without the elastic foundation, the shape functions do not give the same results. Therefore, it can be seen that for the Mantari's and Akavci's shape functions, stresses achieve their maximum values in the plane z/h = 0.25, and for El Meiche's function in the plane z/h = 0.15, while for all the other shape functions as well as new proposed function, maximum values of the stresses are in the plane z/h = 0.2. In order to get a clear insight on the effect of Winkler and Pasternak coefficients of the elastic foundation, Figure 12 shows the diagram of the normalized values of the displacement w plate on the elastic foundation for different values of the index p and coefficients k0 and k1. By comparing the two diagrams, it can be seen that the change of the displacement value w is higher with the increase of the coefficient k1 value than with the increase of the coefficient k0. For example, for the FGM plate when p = 5, and the increase of the coefficient from k0 = 0 to k0 = 100, the value of deflection changes twice its value. In the other case, with the change of the coefficient from k1 = 0 to k1 = 100, the value of deflection changes 8 times its value.

Conclusions
The results obtained in the previously published papers have been a starting point for developing and applying the new shape function. They have emphasized the importance and topicality of the research on the application of the functionally graded materials. A thorough and comprehensive systematization and investigation of the literature on the matter have been conducted according to the problem type which authors tried to solve during FGM plate analysis. Special attention and focus have been given to different deformation theories which authors had used in their analyses. The new shape function has been presented along with the comparative review of it with In order to get a clear insight on the effect of Winkler and Pasternak coefficients of the elastic foundation, Figure 12 shows the diagram of the normalized values of the displacement w plate on the elastic foundation for different values of the index p and coefficients k 0 and k 1. By comparing the two diagrams, it can be seen that the change of the displacement value w is higher with the increase of the coefficient k 1 value than with the increase of the coefficient k 0 . For example, for the FGM plate when p = 5, and the increase of the coefficient from k 0 = 0 to k 0 = 100, the value of deflection changes twice its value. In the other case, with the change of the coefficient from k 1 = 0 to k 1 = 100, the value of deflection changes 8 times its value. In order to get a clear insight on the effect of Winkler and Pasternak coefficients of the elastic foundation, Figure 12 shows the diagram of the normalized values of the displacement w plate on the elastic foundation for different values of the index p and coefficients k0 and k1. By comparing the two diagrams, it can be seen that the change of the displacement value w is higher with the increase of the coefficient k1 value than with the increase of the coefficient k0. For example, for the FGM plate when p = 5, and the increase of the coefficient from k0 = 0 to k0 = 100, the value of deflection changes twice its value. In the other case, with the change of the coefficient from k1 = 0 to k1 = 100, the value of deflection changes 8 times its value.

Conclusions
The results obtained in the previously published papers have been a starting point for developing and applying the new shape function. They have emphasized the importance and topicality of the research on the application of the functionally graded materials. A thorough and comprehensive systematization and investigation of the literature on the matter have been conducted according to the problem type which authors tried to solve during FGM plate analysis. Special attention and focus have been given to different deformation theories which authors had used in their analyses. The new shape function has been presented along with the comparative review of it with

Conclusions
The results obtained in the previously published papers have been a starting point for developing and applying the new shape function. They have emphasized the importance and topicality of the research on the application of the functionally graded materials. A thorough and comprehensive systematization and investigation of the literature on the matter have been conducted according to the problem type which authors tried to solve during FGM plate analysis. Special attention and focus have been given to different deformation theories which authors had used in their analyses. The new shape function has been presented along with the comparative review of it with 13 different shape functions which were primarily developed by different authors for the analysis of composite laminates but, in this paper, they have been adjusted and implemented in appropriate relations for the analysis of FGM plates. Based on the obtained results of the static analysis of moderately thick and thick plates, it can be concluded that the newly developed shape function could be applied in the analysis of FGM plates.
By analyzing the obtained results, the following could be concluded: • the values of the vertical displacement w (deflection) and the corresponding stresses, which were obtained in this paper by using HSDT theory based on the new shape function, match the results of the same values obtained in the reference papers by using TSDT theory [58], quasi 3D theory of elasticity [59] and HSDT theories based on 13 different shape functions. However, in contrast to that, there are significant deviations of the results obtained for the values of the vertical displacement, especially for stresses σ xx , from the results obtained by CPT theory from the reference papers [60]. • the diagram of the distribution of transverse shear stresses τ xz and τ yz across the thickness of the plate shows the difference in behavior between a homogenous, ceramic or metal, plate and FGM plate. A basic property of FGM can be clearly seen, and that is the asymmetry of the stress distribution in relation to the middle plane of the plate (z = 0). The maximum values of stresses, depending on the volume fraction of certain constituents, are shifted in relation to the plane z = 0, which represents a neutral plane in homogenous plates. • the highest values of the displacement w are obtained in a metal plate, the lowest in a ceramic plate and in an FGM plate, the values are somewhere in between and they depend on the volume fraction of the constituents. Based on that, it can be concluded that by varying the volume fraction of metal and ceramic, a desired bending rigidity of the plate can be achieved. • a comparative analysis of the change of transverse shear stresses τ xz and τ yz across the thickness of the plate shows that, unlike the stress τ xy , their values do not match for all the shape functions. Funding: This research received no external funding.