Vibration and buckling analyses of functionally graded plates based on refined plate theory using airy stress function

Abstract

In this paper, an approach based on the refined plate theory and Airy stress function has been proposed to investigate the vibration and buckling behaviors of functionally graded (FG) rectangular plates. The significant feature of the proposed approach is considering only two unknowns in the displacement field in which the contribution due to shear and bending to the total transverse displacement are clarified. Using the extended Hamilton’s principle and defining an Airy stress function corresponding to the compatibility equation, the equations of motion, which do not explicitly include the in-plane displacements, are derived. The accuracy and effectiveness of the current model is shown by comparing the natural frequencies and buckling loads of various FG rectangular plates calculated by the proposed approach with even three-dimensional (3-D) and quasi-3D solutions. Besides, the exact dynamic response of a square FG plate due to a harmonic central force is investigated using modal analysis. This approach is capable of handling quasi-3D models, by selecting proper functions in the displacement field to consider the thickness stretching effect. By doing this, behavior of FG square plates with various coefficient functions are compared with each other. Therefore, through implementing the Airy stress function, the number of variables is reduced, in turn, simplifying the dynamic model. Consequently, it can be used for wide range study of static and dynamic behaviors of FG plates.

Introduction

In laminated composites, stepwise changes of elastic properties and thermal coefficients lead to high inter-laminar stresses causing delamination. FG materials are a class of composites having continuous material properties from one point to another, eliminating not only thermal and residual stresses, but also stress concentrations found in the usual composite materials. They are high thermal, wear and corrosion resistant materials widely used in many engineering applications such as mechanical, aerospace, nuclear, civil, fire retardants, and automotive industries.

Various theories have been presented through the history for investigating the mechanical behavior of plates. The simplest model is based on the classical plate theory (CPT), which ignores both transverse shear and normal deformation effects.

Numerous researchers have attempted to refine the CPT. The first order shear deformation theory (FSDT) developed by Mindlin ¹ accounts for the shear deformation effect with a linear variation of the in-plane displacements through the thickness.

Using FSDT with four variables, Rezaei et al. ² investigated natural frequencies of FG plates with porosities considering two boundary layer functions to decouple the governing equations.

Since in FSDT, a uniform shear strain is assumed through the plate thickness, a shear correction factor is needed. As a result, higher-order shear deformation theories (HSDTs) have been introduced in the literature. These theories do not need shear correction factor and lead to more accurate shear stress distribution.

HSDTs consider higher order variations of in-plane displacements or both in-plane and transverse displacements (i.e., quasi-three-dimensional (quasi-3D) theory) through the thickness. Based on an exact 3D method, Vel and Batra ³ studied vibration of rectangular FG plates. Using various shear deformation theories, Sarangan and Singh ⁴ investigated static, buckling, and free vibration behaviors of laminated composite and sandwich plates. Based on a non-polynomial HSDT, Mantari et al. ⁵ investigated free vibration behavior of FG plates resting on an elastic foundation. Ferreira et al. ⁶ studied free vibration of FG plates based on the FSTD and third-order shear deformation theory (TSDT) using the collocation method with radial basis functions.

Neves et al. ⁷ investigated bending deformation and free vibration of FG plates based on a quasi-3D hyperbolic sine shear deformation theory. Additionally, based on a quasi-3D higher-order shear deformation theory and a meshless technique, Neves et al. ⁸ studied static, free vibration and buckling of isotropic and sandwich FG plates. Thai et al. ⁹ investigated bending and free vibration of FG plates based on a quasi-3D hyperbolic shear deformation theory. Matsunaga ¹⁰ studied natural frequencies and buckling loads of FG plates considering the effect of transverse shear and normal deformation and rotary inertia.

Hellal et al. ¹¹ investigated vibration and buckling of FG sandwich plates supported by elastic foundation in hygro-thermal environment based on a HSDT. Free vibration, buckling, and static responses of graphene-reinforced sandwich plates under mechanical and thermal loadings were investigated by Yaghoobi and Taheri ¹² based on a refined shear deformation plate theory. Other variations of HSDTs to study buckling and vibration behaviors of FG plates under various conditions are also proposed.^13,14

The idea of partitioning the transverse displacement into the bending and shear parts was first proposed by Huffington ¹⁵ and later adopted by Murty ¹⁶ and Senthilnathan et al. ¹⁷ By implementing this idea, not only the number of unknowns is reduced but more importantly the contributions to the total displacement due to shear and bending are clarified. Shimpi ¹⁸ developed a refined plate theory (RPT) for isotropic plates by dividing the displacements into bending and shear components. Thai and Choi ¹⁹ studied bending, buckling, and free vibration behaviors of FG plates resting on elastic foundation using RPT based on the Navier and Levy approach. Free vibration behavior of FG plates based on a four-variable RPT was studied by Benachour et al. ²⁰ Hadji et al. ²¹ analyzed free vibration behavior of rectangular FG sandwich plates, applying four-variable RPT. El Meiche et al. ²² studied buckling and free vibration of thick FG sandwich plates through implementing a hyperbolic shear deformation theory. Thai and Choi ²³ investigated bending and free vibration behaviors of FG plates using various HSDTs, neglecting the effect of thickness stretching. Vu et al. ²⁴ analyzed static deflection, free vibration, and buckling of FG plates with a meshfree method by implementing a refined third-order shear deformation theory. Thermo-mechanical behavior of FG sandwich plates resting on elastic foundation was investigated by Mahmoudi et al. ²⁵ based on a quasi-3D shear deformation theory. Rouzegar et al. ²⁶ studied dynamic behavior of composite laminates integrated with a piezoelectric fiber-reinforced composite (PFRC) actuator layer, based on a four-variable RPT. It was concluded that their approach can predict the dynamic behavior of smart composite laminates accurately. Other kinds of refined theories were also proposed to consider through the thickness deformation effect.^27,28 It should be noted that wide usage of FG materials were studied more in detail for not only rectangular plates but also for other various mechanical models/components, for example, nano-plates, shells, cylinders, disks, beams, vessels, and spheres.^29–37

In this paper, a new method based on the RPT is proposed to investigate the vibration and buckling of rectangular FG plates. This approach does not require a shear correction factor while satisfying shear stress free surface conditions. Furthermore, by using Airy stress function, the in-plane variables are excluded and thus this approach includes only two state variables; transverse bending deformation and transverse shear deformation; compared to HSDTs containing a greater number of unknowns. Moreover, the contributions of shear and bending to the total displacement are completely clarified. The accuracy and effectiveness of the current model is shown by comparing natural frequencies and buckling loads with those computed using other approaches.

Furthermore, through implementing modal analysis, the exact responses of a square FG plate under a harmonic concentrated force are investigated. By using various coefficient functions in the displacement field to model quasi-3D thin and thick FG plates, the variations of transverse displacement and shear stress along the thickness direction are compared with each other by considering only two variables.

Displacement field

The geometry of an FG plate with top and bottom surfaces, being ceramic and metal rich, respectively, is shown in Figure 1. The in-plane coordinates are along x₁ and x₂ directions, while the out of plane coordinate is on the thickness direction z, with zero value in the middle of the plate. H is the plate thickness in z direction and $l_1$ and $l_2$ are the plate lengths in $x_1$ and $x_2$ directions, respectively. Figure 2 shows the deformed configuration of the FG plate with $u_0$ and $v_0$ , the in-plane displacements of the mid-plane in the $x_1$ and $x_2$ directions, respectively. $w_b$ and $w_s$ are transverse deflections of the mid-plane due to bending and shear, respectively.

Figure 1.

Functionally graded plate and coordinate system.

Figure 2.

Deformed and undeformed configuration of the FG plate.

Considering the displacement field used by Senthilnathan et al. ¹⁷ and further improving the model by considering undefined functions, new displacement field is considered to be

$$\begin{array}{l}{u}_1\left(x_1,x_2,z,t\right)=u_0\left(x_1,x_2,t\right)+ g_u\left(z\right)\frac{\partial w_b}{\partial x_1}+f_u\left(z\right)\frac{\partial w_s}{\partial x_1}\\ u_2\left(x_1,x_2,z,t\right)=v_0\left(x_1,x_2,t\right)+g_v\left(z\right)\frac{\partial w_b}{\partial x_2}+f_v\left(z\right)\frac{\partial w_s}{\partial x_2}\\ u_3\left(x_1,x_2,z,t\right)=w_b\left(x_1,x_2,t\right)+s\left(z\right)w_s\left(x_1,x_2,t\right)\end{array}$$

(1)

where $u_1$ , $u_2$ , and $u_3$ , are the displacements of each point in the $x_1$ , $x_2$ , and $z$ directions, respectively, while t shows the time. The second terms in $u_1$ and $u_2$ , containing derivatives of $w_b$ , express the effect of bending on the in-plane displacements as used in the classical plate theory. The third terms, containing derivatives of $w_s$ , show the effect of shear. The transverse displacement $u_3$ consists of two parts, $w_b$ and $w_s$ representing the contribution of bending and shear, respectively.

Functions $g_u\left(z\right), f_u\left(z\right), g_v\left(z\right), f_v\left(z\right)$ , and $s\left(z\right)$ , called coefficient functions in this paper, used in the displacement field equation (1) should be selected in a manner for which the condition of zero transverse shear stress on the free surface is satisfied

$${\sigma }_{13}|\begin{array}{c}\\ z=\pm H/2\end{array}=0,{\sigma }_{23}|\begin{array}{c}\\ z=\pm H/2\end{array}=0$$

(2)

Material properties

Due to gradual change in the volumetric ratio of the constituent materials, functionally graded material (FGM) properties vary continuously along the z direction based on a simple power law distribution in the following form

$$P_{eff}\left(z\right)=P_b+\left(P_t-P_b\right){\left(\frac{1}{2}+\frac{z}{H}\right)}^p\hspace{0.17em}-\raisebox{1ex}{$H$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\le z\le \raisebox{1ex}{$H$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,\hspace{0.17em}0\le p\le \infty$$

(3)

where p is the power law index taking a positive value. $P_{eff}$ is the effective material property of the FG plate, and $P_t$ and $P_b$ are the properties of the top and bottom faces of the plate, respectively. z is the position in the thickness direction of the plate, with zero value in the middle of the plate thickness, while the full thickness of the plate is H.

Strain–displacement relations

Linear strains due to the displacement field equation (1) are

$${\epsilon }_{11}=\frac{\partial u_0}{\partial x_1}+g_u\left(z\right)\frac{{\partial }^2{w}_b}{\partial x_1^2}+f_u\left(z\right)\frac{{\partial }^2{w}_s}{\partial x_1^2}$$

(4)

$${\epsilon }_{22}=\frac{\partial v_0}{\partial x_2}+g_v\left(z\right)\frac{{\partial }^2{w}_b}{\partial x_2^2}+f_v\left(z\right)\frac{{\partial }^2{w}_s}{\partial x_2^2}$$

$${\epsilon }_{33}=\frac{ds\left(z\right)}{dz}{w}_s$$

$${\gamma }_{12}=\frac{\partial u_0}{\partial x_2}+\frac{\partial v_0}{\partial x_1}+\left(g_u\left(z\right)+g_v\left(z\right)\right)\frac{{\partial }^2{w}_b}{\partial x_1\partial x_2}+\left(f_u\left(z\right)+f_v\left(z\right)\right)\frac{{\partial }^2{w}_s}{\partial x_1\partial x_2}$$

$${\gamma }_{13}=\left(1+\frac{d{g}_u\left(z\right)}{dz}\right)\frac{\partial w_b}{\partial x_1}+\left(s\left(z\right)+\frac{d{f}_u\left(z\right)}{dz}\right)\frac{\partial w_s}{\partial x_1}$$

$${\gamma }_{23}=\left(1+\frac{d{g}_v\left(z\right)}{dz}\right)\frac{\partial w_b}{\partial x_2}+\left(s\left(z\right)+\frac{d{f}_v\left(z\right)}{dz}\right)\frac{\partial w_s}{\partial x_2}$$

Constitutive equations

Three-dimensional stress–strain relations can be expressed as

$$\left\{\begin{array}{c}\begin{array}{c}\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ \begin{array}{c}{\sigma }_{33}\\ {\sigma }_{12}\end{array}\end{array}\\ {\sigma }_{23}\end{array}\end{array}\right\}=\left[\begin{array}{cc}\begin{array}{c}\begin{array}{ccc}{C}_{11}& C_{12}& C_{13}\\ C_{21}& C_{22}& C_{23}\\ C_{31}& C_{32}& C_{33}\end{array}\\ \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\end{array}& \begin{array}{c}\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ \begin{array}{ccc}{C}_{44}& 0& 0\\ 0& C_{55}& 0\\ 0& 0& C_{66}\end{array}\end{array}\end{array}\right]\left\{\begin{array}{c}\begin{array}{c}\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{22}\\ \begin{array}{c}{\epsilon }_{33}\\ {\gamma }_{12}\end{array}\end{array}\\ {\gamma }_{13}\end{array}\\ {\gamma }_{23}\end{array}\right\}$$

(5)

For elastic isotropic FGMs, $C_{ij}$ are given in terms of Young’s modulus $E\left(z\right)$ and Poisson’s ratio $\nu$ as

$$\begin{array}{l}{C}_{11}\left(z\right)=C_{22}\left(z\right)=C_{33}\left(z\right)=\frac{\left(1-\nu \right)E\left(z\right)}{\left(1-2\nu \right)\left(1+\nu \right)}\\ C_{12}\left(z\right)=C_{21}\left(z\right)=C_{13}\left(z\right)=C_{31}\left(z\right)=C_{23}\left(z\right)=C_{32}\left(z\right)=\frac{\nu E\left(z\right)}{\left(1-2\nu \right)\left(1+\nu \right)}\\ C_{44}\left(z\right)=C_{55}\left(z\right)=C_{66}\left(z\right)=\frac{E\left(z\right)}{2\left(1+\nu \right)}\end{array}$$

(6)

If the thickness stretching effect is neglected ( ${\epsilon }_{33}=0$ ), the constitutive relations of equation (5) can be rewritten as

$$\left\{\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ \begin{array}{c}{\sigma }_{12}\\ \begin{array}{c}{\sigma }_{13}\\ {\sigma }_{23}\end{array}\end{array}\end{array}\right\}=\left[\begin{array}{ccc}{C}_{11}& C_{12}& \begin{array}{ccc}0& 0& 0\end{array}\\ C_{21}& C_{22}& \begin{array}{ccc}0& 0& 0\end{array}\\ \begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{ccc}\begin{array}{c}{C}_{44}\\ 0\\ 0\end{array}& \begin{array}{c}0\\ C_{55}\\ 0\end{array}& \begin{array}{c}0\\ 0\\ C_{66}\end{array}\end{array}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{22}\\ \begin{array}{c}{\gamma }_{12}\\ \begin{array}{c}{\gamma }_{13}\\ {\gamma }_{23}\end{array}\end{array}\end{array}\right\}$$

(7)

where elastic constants are

$$\begin{array}{l}{C}_{11}\left(z\right)=C_{22}\left(z\right)=\frac{E\left(z\right)}{1-{\nu }^2}\hspace{0.17em}\\ C_{12}\left(z\right)=C_{21}\left(z\right)=\frac{\nu E\left(z\right)}{1-{\nu }^2}{C}_{44}\left(z\right)=C_{55}\left(z\right)=C_{66}\left(z\right)=\frac{E\left(z\right)}{2\left(1+\nu \right)}\end{array}$$

(8)

Governing equations

The extended Hamilton’s principle can be used to obtain the equations of motion

$${\displaystyle {\int }_0^T\left(\delta U+\delta W_e-\delta K\right)dt=0}$$

(9)

where $\delta U$ is the variation of strain energy, $\delta W_e$ is the virtual work done by the applied external loads, and $\delta K$ is the variation of kinetic energy.

Variation of the plate strain energy can be expressed as

$$\delta U={\displaystyle {\int }_A{\displaystyle {\int }_{-H/2}^{H/2}\left({\sigma }_{11}\delta {\epsilon }_{11}+{\sigma }_{22}\delta {\epsilon }_{22}+{\sigma }_{33}\delta {\epsilon }_{33}+{\sigma }_{12}\delta {\gamma }_{12}+{\sigma }_{13}\delta {\gamma }_{13}+{\sigma }_{23}\delta {\gamma }_{23}\right)}dAdz}$$

(10)

The work done by the external transverse uniform distributed force q acting on the top surface of the plate ( $z=H/2$ ) is

$$W_{e,1}={\displaystyle {\int }_A-q\left(w_b+s\left(\frac{H}{2}\right)w_s\right)dA}$$

(11)

The work done by the uniform in-plane compressive distributed loads (i.e., $P_x$ and $P_y$ both negative with the dimension of force/length) applied along the sides at the edges $x_1=0,l_1$ and $x_2=0,l_2$ , respectively, which are due to the slope of the plate along $x_1$ and $x_2$ directions, can be written as

$$W_{e,2}=\frac{1}{2}{\displaystyle {\int }_A\left(P_x{\left(\frac{\partial u_3}{\partial x_1}\right)}^2+P_y{\left(\frac{\partial u_3}{\partial x_2}\right)}^2\right)dA}$$

(12)

Total work done by the external forces is the sum of equations (11) and (12), that is

$$W_e=W_{e,1}+W_{e,2}$$

(13)

Also, the Kinetic energy of the plate is given by

$$K=\frac{1}{2}{\displaystyle {\int }_A{\displaystyle {\int }_{-H/2}^{H/2}\rho \left(z\right)\left({\dot{u_1}}^2+{\dot{u_2}}^2+{\dot{u_3}}^2\right)dAdz}}$$

(14)

where dot-superscript convention represents the partial differentiation w.r.t. time and $\rho \left(z\right)$ is the effective mass density of the FG plate.

Based on the divergence theorem and collecting the coefficients of each variable, the equations of motion are derived as

$$\delta u_0: \frac{\partial N_1^0}{\partial x_1}+\frac{\partial N_{12}^0}{\partial x_2}=I_0{\ddot{u}}_0+I_1\frac{\partial {\ddot{w}}_b}{\partial x_1}+I_2\frac{\partial {\ddot{w}}_s}{\partial x_1}$$

(15)

$$\delta v_0: \frac{\partial N_2^0}{\partial x_2}+\frac{\partial N_{12}^0}{\partial x_1}=I_0{\ddot{v}}_0+I_3\frac{\partial {\ddot{w}}_b}{\partial x_2}+I_4\frac{\partial {\ddot{w}}_s}{\partial x_2}$$

(16)

$$\begin{array}{l}\delta w_b: \frac{{\partial }^2{N}_1^g}{\partial x_1^2}+\frac{{\partial }^2{N}_2^g}{\partial x_2^2}+\frac{{\partial }^2\left(N_{12}^g+N_{21}^g\right)}{\partial x_1\partial x_2}+P_x\left(\frac{{\partial }^2{w}_b}{\partial x_1^2}+\frac{{\partial }^2{w}_s}{\partial x_1^2}\right)\\ +P_y\left(\frac{{\partial }^2{w}_b}{\partial x_2^2}+\frac{{\partial }^2{w}_s}{\partial x_2^2}\right)-\frac{\partial }{\partial x_1}\left(Q_1^{00}+Q_1^g\right)-\frac{\partial }{\partial x_2}\left(Q_2^{00}+Q_2^g\right)-q =\\ I_1\frac{\partial {\ddot{u}}_0}{\partial x_1}+I_3\frac{\partial {\ddot{v}}_0}{\partial x_2}-I_0{\ddot{w}}_b+K_1\frac{{\partial }^2{\ddot{w}}_b}{\partial x_1^2}+K_3\frac{{\partial }^2{\ddot{w}}_b}{\partial x_2^2}-I_5{\ddot{w}}_s+J_1\frac{{\partial }^2{\ddot{w}}_s}{\partial x_1^2} +J_2\frac{{\partial }^2{\ddot{w}}_s}{\partial x_2^2}\end{array}$$

(17)

$$\begin{array}{l}\delta w_s: \frac{{\partial }^2{N}_1^f}{\partial x_1^2}+\frac{{\partial }^2{N}_2^f}{\partial x_2^2}+\frac{{\partial }^2\left(N_{12}^f+N_{21}^f\right)}{\partial x_1\partial x_2}+R^s+P_x\left(\frac{{\partial }^2{w}_b}{\partial x_1^2}+\frac{{\partial }^2{w}_s}{\partial x_1^2}\right)\\ +P_y\left(\frac{{\partial }^2{w}_b}{\partial x_2^2}+\frac{{\partial }^2{w}_s}{\partial x_2^2}\right)-\frac{\partial }{\partial x_1}\left(Q_1^f+Q_1^{s0}\right)-\frac{\partial }{\partial x_2}\left(Q_2^f+Q_2^{s0}\right) -s\left(\frac{H}{2}\right)q=\\ I_2\frac{\partial {\ddot{u}}_0}{\partial x_1}+I_4\frac{\partial {\ddot{v}}_0}{\partial x_2}-I_5{\ddot{w}}_b+J_1\frac{{\partial }^2{\ddot{w}}_b}{\partial x_1^2}+J_2\frac{{\partial }^2{\ddot{w}}_b}{\partial x_2^2}-K_5{\ddot{w}}_s +K_2\frac{{\partial }^2{\ddot{w}}_s}{\partial x_1^2}+K_4\frac{{\partial }^2{\ddot{w}}_s}{\partial x_2^2}\end{array}$$

(18)

where $I_i$ , $J_i$ , and $K_i$ are various inertias of the plate defined by

$$\begin{array}{l}\left(I_0,I_1,I_2,I_3,I_4,I_5\right)={\displaystyle {\int }_{-\frac{H}{2}}^{\frac{H}{2}}\rho \left(z\right)\left(1,g_u\left(z\right),f_u\left(z\right),g_v\left(z\right),f_v\left(z\right),s\left(z\right)\right)dz}\\ \left(J_1,J_2\right)={\displaystyle {\int }_{-\frac{H}{2}}^{\frac{H}{2}}\rho \left(z\right)\left(g_u\left(z\right)f_u\left(z\right),g_v\left(z\right)f_v\left(z\right)\right)dz\left(K_1,K_2,K_3,K_4,K_5\right)=}\\ {\displaystyle {\int }_{-H/2}^{H/2}}\rho \left(z\right)\left({\left(g_u\left(z\right)\right)}^2,{\left(f_u\left(z\right)\right)}^2,{\left(g_v\left(z\right)\right)}^2,{\left(f_v\left(z\right)\right)}^2,{\left(s\left(z\right)\right)}^2\right)dz\end{array}$$

(19)

Stress resultants of the plate shown in equations (15–18) can be obtained by writing the expressions for $\delta U$ , using the strain–displacement relation equations equation (4)

$$\left(N_1^0, N_1^g, N_1^f\right)={\displaystyle {\int }_{-H/2}^{H/2}\left(1, g_u\left(z\right), f_u\left(z\right)\right){\sigma }_{11}dz}$$

(20)

$$\left(N_2^0,N_2^g, N_2^f\right)={\displaystyle {\int }_{-H/2}^{H/2}\left(1,g_v\left(z\right),f_v\left(z\right)\right){\sigma }_{22}dz}$$

$$R^s= {\displaystyle {\int }_{-H/2}^{H/2}\frac{ds\left(z\right)}{dz}{\sigma }_{33}dz}$$

$$\left(N_{12}^0,N_{12}^g,N_{12}^f,N_{21}^g,N_{21}^f\right)= {\displaystyle {\int }_{-H/2}^{H/2}\left(1,g_u\left(z\right),f_u\left(z\right),g_v\left(z\right),f_v\left(z\right)\right){\sigma }_{12}dz}$$

$$\left(Q_1^{00},Q_1^{s0},Q_1^g,Q_1^f\right)={\displaystyle {\int }_{-H/2}^{H/2}\left(1,s\left(z\right),\frac{d{g}_u\left(z\right)}{dz},\frac{d{f}_u\left(z\right)}{dz}\right){\sigma }_{13}dz}$$

$$\left(Q_2^{00},Q_2^{s0},Q_2^g,Q_2^f\right)={\displaystyle {\int }_{-H/2}^{H/2}\left(1,s\left(z\right),\frac{d{g}_v\left(z\right)}{dz},\frac{d{f}_v\left(z\right)}{dz}\right){\sigma }_{23}dz}$$

By substituting the strain-displacement relations equation (4) and the stress–strain relations equation (5) in the stress resultants $N_1^0$ , $N_2^0$ , and $N_{12}^0$ described in equation (20), various coefficients used in equation (22), are obtained as

$$\left({\mathit{A}}_{11}^0,{\mathit{A}}_{11}^{\mathit{g}},{\mathit{A}}_{11}^{\mathit{f}}\right)={\displaystyle {\int }_{-\mathit{H}/2}^{\mathit{H}/2}{\mathit{C}}_{11}\left(\mathit{z}\right)\left(1,{\mathit{g}}_{\mathit{u}}\left(\mathit{z}\right),{\mathit{f}}_{\mathit{u}}\left(\mathit{z}\right)\right)\mathit{dz}}$$

$$\left({\mathit{A}}_{22}^0,{\mathit{A}}_{22}^{\mathit{g}},{\mathit{A}}_{22}^{\mathit{f}}\right)={\displaystyle {\int }_{-\mathit{H}/2}^{\mathit{H}/2}{\mathit{C}}_{22}\left(\mathit{z}\right)\left(1,{\mathit{g}}_{\mathit{v}}\left(\mathit{z}\right),{\mathit{f}}_{\mathit{v}}\left(\mathit{z}\right)\right)\mathit{dz}}$$

$$\left[\begin{array}{ccc}{\mathit{A}}_{12}^0& {\mathit{A}}_{12}^{\mathit{g}}& \begin{array}{cc}{\mathit{A}}_{12}^{\mathit{f}}& \begin{array}{cc}{\mathit{A}}_{12}^{{\mathit{g}}^{\mathit{\text{'}}}}& {\mathit{A}}_{12}^{{\mathit{f}}^{\mathit{\text{'}}}}\end{array}\end{array}\\ {\mathit{A}}_{44}^0& {\mathit{A}}_{44}^{\mathit{g}}& \begin{array}{cc}{\mathit{A}}_{44}^{\mathit{f}}& \begin{array}{cc}{\mathit{A}}_{44}^{{\mathit{g}}^{\mathit{\text{'}}}}& {\mathit{A}}_{44}^{{\mathit{f}}^{\mathit{\text{'}}}}\end{array}\end{array}\end{array}\right]={\displaystyle {\int }_{-\mathit{H}/2}^{\mathit{H}/2}\left\{\begin{array}{c}{\mathit{C}}_{12}\left(\mathit{z}\right)\\ {\mathit{C}}_{44}\left(\mathit{z}\right)\end{array}\right\}\left(1,{\mathit{g}}_{\mathit{u}}\left(\mathit{z}\right),{\mathit{f}}_{\mathit{u}}\left(\mathit{z}\right),{\mathit{g}}_{\mathit{v}}\left(\mathit{z}\right),{\mathit{f}}_{\mathit{v}}\left(\mathit{z}\right)\right)\mathit{dz}}$$

$$\left[\begin{array}{c}{\mathit{D}}_{13}^{\mathit{s}}\\ {\mathit{D}}_{23}^{\mathit{s}}\end{array}\right]={\displaystyle {\int }_{-\mathit{H}/2}^{\mathit{H}/2}\left\{\begin{array}{c}{\mathit{C}}_{13}\left(\mathit{z}\right)\\ {\mathit{C}}_{23}\left(\mathit{z}\right)\end{array}\right\}\left(\frac{\mathit{ds}\left(\mathit{z}\right)}{\mathit{dz}}\right)\mathit{dz}}\mathit{ }$$

(21)

Derivatives of the in-plane variables $u_0$ and $v_0$ can be obtained in terms of the stress resultants $N_1^0$ , $N_2^0$ , and $N_{12}^0$ described in equation (20) and the derivatives of variables. Substituting derivatives of the in-plane variables into the strain-displacement relations equation (4), normal strains along $x_1$ and $x_2$ directions and in-plane shear strain are obtained as

$${\epsilon }_{ij}={\mathrm{e}}_{i,1}{N}_1^0+{\mathrm{e}}_{i,2}{N}_2^0+{\mathrm{e}}_{i,3}\frac{{\mathrm{\partial }}^2{w}_b}{\mathrm{\partial }{x}_1^2}+{\mathrm{e}}_{i,4}\frac{{\mathrm{\partial }}^2{w}_b}{\mathrm{\partial }{x}_2^2}+{\mathrm{e}}_{i,5}\frac{{\mathrm{\partial }}^2{w}_s}{\mathrm{\partial }{x}_1^2}+{\mathrm{e}}_{i,6}\frac{{\mathrm{\partial }}^2{w}_s}{\mathrm{\partial }{x}_2^2}+{\mathrm{e}}_{i,7}{\mathit{\omega }}_s\left(i=j=\mathrm{1,2}\right)$$

$${\gamma }_{12}={\mathrm{\epsilon }}_1{N}_{12}^0+{\mathrm{\epsilon }}_2\frac{{\partial }^2{w}_b}{\partial x_1\partial x_2}+{\mathrm{\epsilon }}_3\frac{{\partial }^2{w}_s}{\partial x_1\partial x_2}$$

(22)

where coefficients ${\mathrm{e}}_{i,j}$ and ${\mathrm{\epsilon }}_i$ can be presented as

$$\begin{array}{l}{\mathrm{e}}_{\mathrm{1,1}}=\frac{A_{22}^0}{A_{11}^0{A}_{22}^0-{{\displaystyle A_{12}^0}}^2}, {\mathrm{e}}_{\mathrm{1,2}}=-\frac{A_{12}^0}{A_{11}^0{A}_{22}^0-{{\displaystyle A_{12}^0}}^2}\\ {\mathrm{e}}_{\mathrm{1,3}}=\frac{A_{12}^0{A}_{12}^g-A_{22}^0{A}_{11}^g}{A_{11}^0{A}_{22}^0-{{\displaystyle A_{12}^0}}^2}+g_u\left(\text{z}\right), {\mathrm{e}}_{\mathrm{1,4}}=\frac{A_{12}^0{A}_{22}^g-A_{22}^0{A}_{12}^{g^{\text{'}}}}{A_{11}^0{A}_{22}^0-{{\displaystyle A_{12}^0}}^2}{\mathrm{e}}_{\mathrm{1,5}}=\frac{A_{12}^0{A}_{12}^f-A_{22}^0{A}_{11}^f}{A_{11}^0{A}_{22}^0-{{\displaystyle A_{12}^0}}^2}+f_u\left(\text{z}\right),\\ {\mathrm{e}}_{\mathrm{1,6}}=\frac{A_{12}^0{A}_{22}^f-A_{22}^0{A}_{12}^{f^{\text{'}}}}{A_{11}^0{A}_{22}^0-{{\displaystyle A_{12}^0}}^2}{\mathrm{e}}_{\mathrm{1,7}}=\frac{A_{12}^0{D}_{23}^s-A_{22}^0{D}_{13}^s}{A_{11}^0{A}_{22}^0-{{\displaystyle A_{12}^0}}^2}, {\mathrm{e}}_{\mathrm{2,1}}=-\frac{A_{12}^0}{A_{11}^0{A}_{22}^0-{{\displaystyle A_{12}^0}}^2}{\mathrm{e}}_{\mathrm{2,2}}=\frac{A_{11}^0}{A_{11}^0{A}_{22}^0-{{\displaystyle A_{12}^0}}^2}, \\ {\mathrm{e}}_{\mathrm{2,3}}=\frac{A_{12}^0{A}_{11}^g-A_{11}^0{A}_{12}^g}{A_{11}^0{A}_{22}^0-{{\displaystyle A_{12}^0}}^2} {\mathrm{e}}_{\mathrm{2,4}}=\frac{A_{12}^0{A}_{12}^{g^{\text{'}}}-A_{11}^0{A}_{22}^g}{A_{11}^0{A}_{22}^0-{{\displaystyle A_{12}^0}}^2}+g_v\left(\text{z}\right), {\mathrm{e}}_{\mathrm{2,5}}=\frac{A_{12}^0{A}_{11}^f-A_{11}^0{A}_{12}^f}{A_{11}^0{A}_{22}^0-{{\displaystyle A_{12}^0}}^2} {\mathrm{e}}_{\mathrm{2,6}}=\frac{A_{12}^0{A}_{12}^{f^{\text{'}}}-A_{11}^0{A}_{22}^f}{A_{11}^0{A}_{22}^0-{{\displaystyle A_{12}^0}}^2}+\\ f_v\left(\text{z}\right), {\mathrm{e}}_{\mathrm{2,7}}=\frac{A_{12}^0{D}_{13}^s-A_{11}^0{D}_{23}^s}{A_{11}^0{A}_{22}^0-{{\displaystyle A_{12}^0}}^2}{\mathrm{\epsilon }}_1=\frac{1}{A_{44}^0}, {\mathrm{\epsilon }}_2=g_u\left(\text{z}\right)+g_v\left(\text{z}\right)-\frac{A_{44}^g+A_{44}^{g^{\prime }}}{A_{44}^0}{\mathrm{\epsilon }}_3=f_u\left(\text{z}\right)+f_v\left(\text{z}\right)-\frac{A_{44}^f+A_{44}^{f^{\prime }}}{A_{44}^0}\end{array}$$

(23)

Solution procedure

The compatibility equation of an FG plate can be expressed as

$$\frac{{\partial }^2{\text{ε}}_{11}}{\partial x_2^2}+\frac{{\partial }^2{\text{ε}}_{22}}{\partial x_1^2}-\frac{{\partial }^2{\text{γ}}_{12}}{\partial x_1\partial x_2}=0$$

(24)

$f$ , the stress function, should be selected in such a way that both the left side of equation (15) and the left side of equation (16) become zero. In other words, various in-plane stress resultants can be obtained as

$$N_1^0=\frac{{\partial }^{2}f}{\partial x_2^2}, N_2^0=\frac{{\partial }^{2}f}{\partial x_1^2}, N_{12}^0=- \frac{{\partial }^{2}f}{\partial x_1\partial x_2}$$

(25)

Therefore, the in-plane accelerations can be written as

$$\begin{array}{l}{\ddot{u}}_0=-\frac{I_1}{I_0}\frac{\partial {\ddot{w}}_b}{\partial x_1}-\frac{I_2}{I_0}\frac{\partial {\ddot{w}}_s}{\partial x_1}\\ {\ddot{v}}_0=-\frac{I_3}{I_0}\frac{\partial {\ddot{w}}_b}{\partial x_2}-\frac{I_4}{I_0}\frac{\partial {\ddot{w}}_s}{\partial x_2}\end{array}$$

(26)

By substituting equation (26) into equations (17) and (18), it can be easily seen that the equations of transverse motion are dependent of only $w_b$ and $w_s$ . Thus, the in-plane accelerations are taken into account in equations (17) and (18), while their related equations are omitted. Therefore, the total number of equations and variables is reduced without neglecting the effect of in-plane motions.

For movable simply supported boundary conditions, the moments at plate boundaries are zero and other conditions are

$$\begin{array}{l}{w}_b=w_s=0, N_1^0=N_{x0} at x_1=0,l_1\\ w_b=w_s=0, N_2^0=N_{y0} at x_2=0,l_2\end{array}$$

(27)

where $N_{x0}$ and $N_{y0}$ are compressive loads per unit length applied on the four sides of the plate.

The series to satisfy the mentioned boundary conditions can be selected as

$$\begin{array}{l}{w}_b\left(x_1,x_2,t\right)={\displaystyle \sum _{m_b=1}^{\infty }{\displaystyle \sum _{n_b=1}^{\infty }w{b}_{n_b{m}_b}\left(t\right)\sin \left(\frac{n_b\pi x_1}{l_1}\right)\sin \left(\frac{m_b\pi x_2}{l_2}\right)}}\\ w_s\left(x_1,x_2,t\right)={\displaystyle \sum _{m_s=1}^{\infty }{\displaystyle \sum _{n_s=1}^{\infty }w{s}_{n_s{m}_s}\left(t\right)\text{sin}\left(\frac{n_s\pi x_1}{l_1}\right)\text{sin}\left(\frac{m_s\pi x_2}{l_2}\right)}}\end{array}$$

(28)

where $w{b}_{n_b{m}_b}$ and $w{s}_{n_s{m}_s}$ are time dependent variables.

By using the expression of normal strains and in-plane shear strains in terms of $N_1^0$ , $N_2^0$ , and $N_{12}^0$ and variables of the problem described in equation (22) with the aid of definition of a proper stress function based on equation (25), the compatibility equation (24) gives the equation in terms of the stress function $f$ , bending variable $w_b$ , and shear variable $w_s$ . Using the series solution of variables (equation (28)), the compatibility equation gives the stress function f in the following form

$$f\left(x_1,x_2,t\right)={\displaystyle \sum }_{q=1}^{\infty }{\displaystyle \sum }_{p=1}^{\infty }{f}_{pq}\left(t\right)\sin \left(\frac{p\pi x_1}{l_1}\right)\sin \left(\frac{q\pi x_2}{l_2}\right)+\frac{1}{2}{C}_x{x}_1^2+\frac{1}{2}{C}_y{x}_2^2$$

(29)

where $f_{pq}$ is a function of $w{b}_{n_b{m}_b}$ and $w{s}_{n_s{m}_s}$ , obtained in the following form

$$\begin{array}{l}{f}_{pq}=-\frac{1}{{\mathrm{e}}_{\mathrm{2,2}}{\left(\frac{p\pi }{l_1}\right)}^4+\left({\mathrm{e}}_{\mathrm{1,2}}+{\mathrm{\epsilon }}_1+{\mathrm{e}}_{\mathrm{2,1}}\right){\left(\frac{p\pi }{l_1}\right)}^2{\left(\frac{q\pi }{l_2}\right)}^2+{\mathrm{e}}_{\mathrm{1,1}}{\left(\frac{q\pi }{l_2}\right)}^4}\\ \left\{\begin{array}{l}\left({\mathrm{e}}_{\mathrm{2,3}}{\left(\frac{n_b\pi }{l_1}\right)}^4+\left({\mathrm{e}}_{\mathrm{1,3}}-{\mathrm{\epsilon }}_2+{\mathrm{e}}_{\mathrm{2,4}}\right){\left(\frac{n_b\pi }{l_1}\right)}^2{\left(\frac{m_b\pi }{l_2}\right)}^2+{\mathrm{e}}_{\mathrm{1,4}}{\left(\frac{m_b\pi }{l_2}\right)}^4\right)\\ w{b}_{n_b{m}_b}\delta \left(n_b-p\right)\delta \left(m_b-q\right)+ \left({\mathrm{e}}_{\mathrm{2,5}}{\left(\frac{n_s\pi }{l_1}\right)}^4+\left({\mathrm{e}}_{\mathrm{1,5}}-{\mathrm{\epsilon }}_3+{\mathrm{e}}_{\mathrm{2,6}}\right){\left(\frac{n_s\pi }{l_1}\right)}^2{\left(\frac{m_s\pi }{l_2}\right)}^2+{\mathrm{e}}_{\mathrm{1,6}}{\left(\frac{m_s\pi }{l_2}\right)}^4-{\mathrm{e}}_{\mathrm{2,7}}{\left(\frac{n_s\pi }{l_1}\right)}^2-{\mathrm{e}}_{\mathrm{1,7}}{\left(\frac{m_s\pi }{l_2}\right)}^2\right)\\ w{s}_{n_s{m}_s}\delta \left(n_s-p\right)\delta \left(m_s-q\right)\end{array}\right\}\end{array}$$

(30)

In equation (3) $\delta$ is Dirac delta function. $C_x$ and $C_y$ shown in equation (29) are caused by uniform in-plane compressive loads $P_x$ and $P_y$ , applied at edges $x_1=0,l_1$ , and $x_2=0,l_2$ as

$$N_1^0=P_x at x_1=0,l_1 , N_2^0=P_y at x_2=0,l_2$$

(31)

Thus $C_x$ and $C_y$ can be obtained by averaging shown by the following equations

$$\begin{array}{l}{C}_x=\frac{1}{l_1}{\displaystyle {\int }_0^{l_1}\left(P_y\right)d{x}_1=P_y}\\ C_y=\frac{1}{l_2}{\displaystyle {\int }_0^{l_2}\left(P_x\right)d{x}_2=P_x}\end{array}$$

(32)

Based on the normal strains and in-plane shear strains relations equation (22) described in terms of $N_1^0$ , $N_2^0$ , $N_{12}^0$ , $w_b$ , and $w_s$ , various stress resultants (equation (20)) are obtained without including the in-plane variables, whereas $N_1^0$ , $N_2^0$ , and $N_{12}^0$ are related to the stress function $f$ described in equation (25). Based on equations (29) and (30), $f$ is a function of variables $w_b$ and $w_s$ . Hence, the stress resultants are obtained only in terms of $w_b$ and $w_s$ .

Furthermore, using proper forms of the coefficient functions used in the displacement field (equation (1)) and using the series solutions of variables described in equation (28) with any number of terms, equations (17) and (18), lead to the following matrix form of approximate equations of motion

$$\mathrm{M}\left\{\begin{array}{c}{\ddot{\mathit{w}b}}_{n_b{m}_b}\\ {\ddot{\mathit{w}s}}_{n_s{m}_s}\end{array}\right\}+\mathrm{K}\left\{\begin{array}{c}\mathit{w}{b}_{n_b{m}_b}\\ \mathit{w}{s}_{n_s{m}_s}\end{array}\right\}=\mathit{F}$$

(33)

where $\mathrm{M}$ , $\mathrm{K}$ , and $\mathit{F}$ are the mass matrix, the stiffness matrix, and the force vector, respectively. Therefore, by using only two continuous variables $w_b$ and $w_s$ , which can be each approximated by any arbitrary number of terms, the behavior of FG plates under various loads can be investigated efficiently and effectively.

Results and discussions

To demonstrate the validity and effectiveness of the proposed approach, two different FG materials are selected. Material properties of constituents of Al/ZrO₂ and Al/Al₂O₃ FG materials are shown in Table 1.

Table 1.

Material properties of FG plate constituents.

Material property	Material Set-1		Material Set-2
	Al, Metal layer	ZrO2, Ceramic layer	Al, Metal layer	Al2O3, Ceramic layer
E (GPa)	70	200	70	380
$\rho$ (Kg/m3)	2702	5700	2702	3800
$\nu$	0.3	0.3	0.3	0.3

In this study, functions used in the displacement field to satisfy equation ( 2 ) are chosen similar to those used by Senthilnathan et al. ¹⁷ as

$$g_u\left(z\right)=-z, f_u\left(z\right)=-\frac{4z^3}{3H^2}, g_v\left(z\right)=-z, f_v\left(z\right)=-\frac{4z^3}{3H^2}, s\left(z\right)=1$$

(34)

Hence ${\epsilon }_{33}$ described in equation (4) becomes zero and the constitutive equation (equation (7)) must be used to study the behavior of FG plates.

Free vibration analysis

In order to see how well the proposed approach works, non-dimensional fundamental frequency of a movable square simply supported Al/ZrO₂ FG plate (Material Set-1), with various thickness and power index, is obtained. The dimensionless frequency in the form of $\widehat{\omega }={\omega }_{1}H\sqrt{{\rho }_b/E_b}$ with ${\omega }_1$ as the first natural frequency of the FG plate is used.

In Table 2, the non-dimensional fundamental frequency obtained by the proposed approach is compared with those obtained by 3D solutions of Vel and Batra, ³ quasi-3D solutions,^7–9 TSDT solutions, ⁶ and HSDT solutions. ²³ As it is shown in the second column of Table 2, References 6 and 23, like our work, neglected the thickness-stretching effect. Number of variables used by each approach is shown in the third column. The numbers shown in the parenthesis present the percentage differences between the results of various studies with those of 3D solutions ³ as the most accurate results.

Table 2.

Dimensionless fundamental frequency $\widehat{\omega }={\omega }_{1}H\sqrt{{\rho }_b/E_b}$ of square Al/ZrO₂ FG plate.

Method	${\epsilon }_{33}$	No. of Variables	$p$ = 0		$p$ = 1			$l_1/H$ = 5
			$l_1/H$ = $\sqrt{10}$ .	$l_1/H$ = 10	$l_1/H$ = 5	$l_1/H$ = 10	$l_1/H$ = 20	$p$ = 2	$p$ = 3	$p$ = 5
3D exact 3	$\ne$ 0	—	0.4658	0.0578	0.2192	0.0596	0.0153	0.2197	0.2211	0.2225
Quasi-3D 7	$\ne$ 0	9	—	—	0.2193 (0.046%)	0.0596 (0.000%)	0.0153 (0.000%)	0.2201 (0.182%)	0.2216 (0.226%)	0.2230 (0.225%)
Quasi-3D 8	$\ne$ 0	9	—	—	0.2193 (0.046%)	—	—	0.2200 (0.136%)	0.2215 (0.181%)	0.2230 (0.225%)
Quasi-3D 9	$\ne$ 0	5	0.4661 (0.065%)	0.0578 (0.000%)	0.2192 (0.000%)	0.0597 (0.168%)	0.0153 (0.000%)	0.2201 (0.182%)	0.2214 (0.136%)	0.2225 (0.000%)
TSDT 6	=0	5	—	—	0.2188 (0.182%)	0.0592 (0.671%)	0.0147 (3.921%)	0.2188 (0.409%)	0.2202 (0.407%)	0.2215 (0.449%)
HSDT 23	=0	4	0.4622 (0.773%)	0.0577 (0.173%)	0.2169 (1.049%)	0.0592 (0.671%)	0.0152 (0.653%)	0.2178 (0.865%)	0.2193 (0.814%)	0.2206 (0.854%)
Present	=0	2	0.4623 (0.751%)	0.0577 (0.173%)	0.2169 (1.049%)	0.0592 (0.671%)	0.0152 (0.653%)	0.2178 (0.865%)	0.2193 (0.814%)	0.2206 (0.854%)

By softening the system through increasing the value of the power law index $p$ , the fundamental frequency decreases. As Table 2 shows, the natural frequencies of thin plates obtained by the proposed approach considering only two variables, namely transverse bending deformation and transverse shear deformation, are in good agreement even with the 3D solutions. ³ Showing that implementing Airy stress function simplifies the model in an effective way. The natural frequencies obtained by the proposed approach, just like those obtained by References 6 and 23, are slightly underestimated due to ignoring the thickness-stretching effect. Though the thinner the plate, the lower the difference percentages become.

Buckling analysis

The effect of transverse shear deformation variable $w_s$ on the dimensionless buckling behavior of square Al/Al₂O₃ FG plates (Material Set-2) is investigated in this part. The dimensionless buckling load along $x_1$ direction obtained by the proposed approach is compared in Table 3 with the results of Matsunaga ¹⁰ in which 29 variables were used to model the dynamic displacement field, including the effects of transverse shear and normal deformations and rotary inertia. The relation $\tilde{N}=\frac{P_{cr}}{H{E}_t}$ is used to obtain the results where $P_{cr}$ is the critical buckling load of the FG plate.

Table 3.

Comparison of dimensionless buckling load $\tilde{N}=P_{cr}/H{E}_t$ with and without considering transverse shear deformation for square Al/Al₂O₃ FG plate.

$\frac{l_1}{H}$ .	$p$	$\tilde{N}$
		Matsunaga, 10 29 Variables ${\epsilon }_{33}\ne 0$ .	Present, RPT ( $w_b,w_s$ ), 2 Variables ${\epsilon }_{33}=0$ .	Present, CPT ( $w_b$ ), 1 Variable ${\epsilon }_{33}=0$ .
2	0	0.3581	0.3346 (6.562%)	0.9038 (152.388%)
	0.5	0.2482	0.2420 (2.498%)	0.5871 (136.543%)
	1	0.1938	0.1933 (0.258%)	0.4505 (132.456%)
	4	0.1120	0.1034 (7.678%)	0.3059 (173.125%)
	10	0.08904	0.07725 (13.241%)	0.2708 (204.133%)
5	0	0.1140	0.1180 (3.509%)	0.1446 (26.842%)
	0.5	0.07571	0.07843 (3.593%)	0.09393 (24.066%)
	1	0.05826	0.06060 (4.016%)	0.07208 (23.721%)
	4	0.03721	0.03894 (4.649%)	0.04894 (31.524%)
	10	0.03183	0.03301 (3.707%)	0.04332 (36.098%)
10	0	0.03381	0.03422 (1.212%)	0.03615 (6.921%)
	0.5	0.02214	0.02238 (1.084%)	0.02348 (6.052%)
	1	0.01698	0.01720 (1.296%)	0.01802 (6.125%)
	4	0.01131	0.01149 (1.591%)	0.01223 (8.134%)
	10	0.009905	0.01004 (1.363%)	0.01083 (9.339%)

Numbers shown in the parenthesis demonstrate the percentage differences of the results obtained by the proposed approach (with and without considering shear deformation effect $w_s$ ) with Matsunaga. ¹⁰ It is evident that the results of the current study considering only two variables, are in good agreement with those of Matsunaga. ¹⁰

Considering only the bending deformation variable $w_b$ results in higher errors specifically in thick and even moderately thick plates. Thus, higher order theories must be used for modeling the aforementioned plates.

Harmonically excited vibration

Consider a harmonic concentrated force $q=f\sin \left(\omega t\right)\delta \left(x_1-\frac{l_1}{2}\right)\delta \left(x_2-\frac{l_2}{2}\right)$ with frequency of ω and magnitude of $f$ applied at the FG plate center with geometrical properties: $l_1=l_2=0.2 \left(m\right), H=0.01 \left(m\right), p=2$ . The mechanical properties are shown in Table 1 (Material Set_1). To decouple the equations of motion Eq. Error! Reference source not found., the modal matrix of the system (i.e., P) is used. Assuming X as the vector containing finite approximate degrees of freedom (i.e., X= $\left\{\begin{array}{c}w{b}_{n_b{m}_b}\\ w{s}_{n_s{m}_s}\end{array}\right\}$ ), through implementing the coordinate transformation X=PY, the system of equations is decoupled using orthogonality of the normal modes (modal analysis). Then, the exact solution of each degree of freedom can be written as

$$y_i\left(t\right)=y_i\left(0\right)\cos \left({\omega }_{i}t\right)+\frac{1}{{\omega }_i}\left(\dot{y_i}\left(0\right)-\frac{f_i}{k_{ii}}\frac{\omega }{1-{\left(\frac{\omega }{{\omega }_i}\right)}^2}\right)\sin \left({\omega }_{i}t\right)+\frac{f_i}{k_{ii}}\frac{\sin \left(\omega t\right)}{1-{\left(\frac{\omega }{{\omega }_i}\right)}^2}$$

(35)

where ${\omega }_i$ is the i-th natural frequency of the system, $y_i\left(0\right)$ and $\dot{y_i}\left(0\right)$ are the initial conditions of the i-th degree of freedom, k _ii is the i-th generalized stiffness, and f _i is the i-th generalized force.

Assuming zero initial velocity, the central deflection of the plate versus time for various frequencies and amplitudes of load is plotted in Figure 3, considering the first four symmetric terms for variables $w_b$ and $w_s$ . Initial condition of the plate considered here is $w_{b_{\mathrm{1,1}}}=0.1\times H$ while other degrees of freedom are assumed to be zero. First natural frequency of the plate is ${\omega }_1=2\pi \times 1281.53$ . A beating phenomenon can be seen occurring when $\omega =0.9\times {\omega }_1$ , because of the frequency of the applied load being close to the first natural frequency of the FG plate.

Figure 3.

The effect of frequency and amplitude of the concentrated harmonic force on the dynamic responses of the square Al/ZrO₂ FG plate.

Furthermore, for the frequencies of the applied force lower than the natural frequency of the system, as the load amplitude increases, the response of the system shifts to the left, while for frequencies bigger than the natural frequency of the system, it is vice versa.

Figure 4 shows the central deflection of the aforementioned plate with zero initial displacement and velocity (with various power law indices $p$ ) under a harmonic central force with magnitude of 15 (KN) and a frequency half of the first natural frequency of the system for each power law index. As can be seen the maximum response of the system is increased by increasing $p$ ; which means that the flexibility of the plate is enhanced subsequently.

Figure 4.

Dynamic responses of the square Al/ZrO₂ FG plate for various power law indices $\left(\omega =1/2{\omega }_1\right)$ .

Coefficient function analysis

According to the displacement field equation (1), coefficient functions can be selected in a way that not only the condition of zero transverse shear stress equation (2) is satisfied, but also ${\epsilon }_{33}\ne 0$ . This means that by considering only two variables, the effect of thickness stretching along the $z$ direction is also included, which leads to more accurate results for moderately thick and thick plates.

In this part, transverse displacement and transverse shear stress under the sinusoidal distributed load $q={\displaystyle \sum }_{n=1}^{\infty }{\displaystyle \sum }_{m=1}^{\infty }{q}_0\text{sin}\left(\frac{n\pi x_1}{l_1}\right)\text{sin}\left(\frac{m\pi x_2}{l_2}\right)$ , with $q_0$ as the maximum load, obtained by using various coefficient functions are compared with each other.

The coefficient functions shown in Table 4 are extracted from studies with and without considering the thickness stretching effect.^17,27–30 It is worth mentioning that $g_u\left(z\right)=g_v\left(z\right)=-z$ . The comparison is done for a square FG plate with Material Set-2 properties. Dimensionless transverse shear stress ${\overline{\sigma }}_{13}\left(z\right)=\frac{H}{q_0{l}_1}{\sigma }_{13}\left(0,l_2/2,z\right)$ and transverse displacement ${\overline{u}}_3\left(z\right)=\frac{10E_t{H}^3}{q_0{l}_1^4}{u}_3\left(l_1/2,l_2/2,z\right)$ are used to obtain the results.

Table 4.

Coefficient functions used in the displacement field.

$f_u\left(z\right)=f_v\left(z\right)= f\left(z\right)$ .	$s\left(z\right)$	${\epsilon }_{33}$
$-\frac{4z^3}{3H^2}$	1	$=0$
$t{g}^{-1}\left(\sin \left(\frac{\pi z}{H}\right)\right)-z$	1	$=0$
$\text{Hsinh}\left(\frac{z}{H}\right)-\frac{4z^3}{3H^2}\text{cosh}\left(\frac{1}{2}\right)$	$\frac{1}{12}f\prime \left(z\right)$	$\ne 0$
$\frac{\pi }{H}z-\frac{9\pi }{5H^3}{z}^3+\frac{28\pi }{25H^5}{z}^5$	$\frac{1}{8}f\prime \left(z\right)$	$\ne 0$
$-8z+\frac{10}{H^2}{z}^3+\frac{6}{5H^4}{z}^5+\frac{8}{7H^6}{z}^7$	$\frac{3}{20}f\prime \left(z\right)$	$\ne 0$

Figure 5 and Figure 6 show variations of ${\overline{u}}_3$ and ${\overline{\sigma }}_{13}$ along the thickness direction of the thick square FG plate (e.g., $l_1/H=2$ ) based on the five coefficient functions shown in Table 4 using various power law indices.

Figure 5.

Variation of ${\overline{u}}_3$ across the thickness of thick square Al/ZrO₂ FG plate under sinusoidal distributed load considering various power law indices ( $l_1/H=2$ ), (a) $p=0$ , (b) $p=1$ , (c) $p=4$ , (d) $p=10$ .

Figure 6.

Variation of ${\overline{\sigma }}_{13}$ across the thickness of thick square Al/ZrO₂ FG plate under sinusoidal distributed load considering various power law indices ( $l_1/H=2$ ), (a) $p=0$ , (b) $p=1$ , (c) $p=4$ , (d) $p=10$ .

As can be seen in Figure 5, for functions with $s\left(z\right)=1$ and thus ${\epsilon }_{33}=0$ , a straight line with no variation for ${\overline{u}}_3$ is observed. Although in quasi-3D models ( ${\epsilon }_{33}\ne 0$ ), curved variations with smaller transverse displacement values compared to simple RPTs can be seen. Furthermore, it is shown in Figure 6 that quasi-3D models predict smoother behavior with smaller transverse shear stress across the thickness. For the homogenous plate ( $p=0$ ), ${\overline{\sigma }}_{13}$ variation is parabolic along the thickness direction.

Figure 7 and Figure 8 shows variation of ${\overline{u}}_3$ and ${\overline{\sigma }}_{13}$ along the thickness direction for thin and thick square Al/ZrO₂ FG plates with $p=10$ and various thickness ratios based on the five previously mentioned coefficient functions (Table 4). For thick plates ( $l_1/H <2$ ), noticeable differences in ${\overline{u}}_3$ among the five curves are observed. It should be noted that the zero transverse shear stress condition on free surfaces of the plate (e.g., $z/H=\pm 1/2$ ) is satisfied for all models indifferent of whether the thickness stretching effect is considered or not. Moreover, as the plate becomes thinner, discrepancies among the models for both ${\overline{u}}_3$ and ${\overline{\sigma }}_{13}$ values are reduced, meaning that the stretching effect has the most influence on thick ones.

Figure 7.

Variation of ${\overline{u}}_3$ across the thickness of square Al/ZrO₂ FG plate ( $p=10$ ), under sinusoidal distributed load considering various plate thickness ratios ( $l_1/H=2, 5, 10, 20$ ).

Figure 8.

Variation of ${\overline{\sigma }}_{13}$ across the thickness of square Al/ZrO₂ FG plate ( $p=10$ ), under sinusoidal distributed load considering various plate thickness ratios ( $l_1/H=2, 5, 10, 20$ ).

Conclusion

An approach based on the refined plate theory was proposed for vibration and buckling analyses of FG plates with movable simply supported boundary conditions. By using the Airy stress function, the in-plane displacements and corresponding equations were omitted and only the transverse equation of motion containing two variables, transverse bending and transverse shear deformation, remained.

Validation and effectiveness of the proposed approach was demonstrated by comparing the obtained fundamental frequencies and buckling loads with those obtained using other studies containing a greater number of unknown variables, whether excluding or including the thickness stretching effect. It has been concluded that the refined plate theory containing only two unknowns is efficient enough to investigate the linear vibration and buckling of thin and thick FG plates.

Using modal analysis, the dynamic response of a square FG plate due to a harmonic central force was investigated. As the frequency of the applied load moves closer to the first natural frequency of the plate, a beating phenomenon can be observed.

The proposed approach is capable of handling wide variety of other theories with selecting proper coefficient functions defined in the displacement field.

In order to use quasi-3D modeling of thin and thick FG plates, various coefficient functions in the displacement field were used. Variations of transverse displacement and shear stress along the thickness direction, with and without considering the thickness stretching effect, were studied. It was shown that quasi-3D models lead to not only smaller transverse displacement and shear stress values, but also smoother stress behavior.

Since by using the proposed approach including only two continuous variables, FG plate modeling is simplified, this approach is a powerful tool for studying dynamical behavior of various FG plates.