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.

Keywords

refined plate theory airy stress function thickness stretching vibration buckling functionally graded plate

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} (x_{1}, x_{2}, z, t) = u_{0} (x_{1}, x_{2}, t) + g_{u} (z) \frac{\partial w_{b}}{\partial x_{1}} + f_{u} (z) \frac{\partial w_{s}}{\partial x_{1}} \\ u_{2} (x_{1}, x_{2}, z, t) = v_{0} (x_{1}, x_{2}, t) + g_{v} (z) \frac{\partial w_{b}}{\partial x_{2}} + f_{v} (z) \frac{\partial w_{s}}{\partial x_{2}} \\ u_{3} (x_{1}, x_{2}, z, t) = w_{b} (x_{1}, x_{2}, t) + s (z) w_{s} (x_{1}, x_{2}, t) \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} (z), f_{u} (z), g_{v} (z), f_{v} (z)$ , and $s (z)$ , 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

σ_{13} | \begin{matrix} z = \pm H / 2 \end{matrix} = 0, σ_{23} | \begin{matrix} z = \pm H / 2 \end{matrix} = 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_{e f f} (z) = P_{b} + (P_{t} - P_{b}) {(\frac{1}{2} + \frac{z}{H})}^{p} - \frac{H}{2} \leq z \leq \frac{H}{2}, 0 \leq p \leq \infty

(3)

where p is the power law index taking a positive value.

P_{e f f}

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

ε_{11} = \frac{\partial u_{0}}{\partial x_{1}} + g_{u} (z) \frac{\partial^{2} w_{b}}{\partial x_{1}^{2}} + f_{u} (z) \frac{\partial^{2} w_{s}}{\partial x_{1}^{2}}

(4)

ε_{22} = \frac{\partial v_{0}}{\partial x_{2}} + g_{v} (z) \frac{\partial^{2} w_{b}}{\partial x_{2}^{2}} + f_{v} (z) \frac{\partial^{2} w_{s}}{\partial x_{2}^{2}}

ε_{33} = \frac{d s (z)}{d z} w_{s}

γ_{12} = \frac{\partial u_{0}}{\partial x_{2}} + \frac{\partial v_{0}}{\partial x_{1}} + (g_{u} (z) + g_{v} (z)) \frac{\partial^{2} w_{b}}{\partial x_{1} \partial x_{2}} + (f_{u} (z) + f_{v} (z)) \frac{\partial^{2} w_{s}}{\partial x_{1} \partial x_{2}}

γ_{13} = (1 + \frac{d g_{u} (z)}{d z}) \frac{\partial w_{b}}{\partial x_{1}} + (s (z) + \frac{d f_{u} (z)}{d z}) \frac{\partial w_{s}}{\partial x_{1}}

γ_{23} = (1 + \frac{d g_{v} (z)}{d z}) \frac{\partial w_{b}}{\partial x_{2}} + (s (z) + \frac{d f_{v} (z)}{d z}) \frac{\partial w_{s}}{\partial x_{2}}

Constitutive equations

Three-dimensional stress–strain relations can be expressed as

{\begin{matrix} \begin{matrix} \begin{matrix} σ_{11} \\ σ_{22} \\ \begin{matrix} σ_{33} \\ σ_{12} \end{matrix} \end{matrix} \\ σ_{23} \end{matrix} \end{matrix}} = [\begin{matrix} \begin{matrix} \begin{matrix} C_{11} & C_{12} & C_{13} \\ C_{21} & C_{22} & C_{23} \\ C_{31} & C_{32} & C_{33} \end{matrix} \\ \begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} \end{matrix} & \begin{matrix} \begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} \\ \begin{matrix} C_{44} & 0 & 0 \\ 0 & C_{55} & 0 \\ 0 & 0 & C_{66} \end{matrix} \end{matrix} \end{matrix}] {\begin{matrix} \begin{matrix} \begin{matrix} ε_{11} \\ ε_{22} \\ \begin{matrix} ε_{33} \\ γ_{12} \end{matrix} \end{matrix} \\ γ_{13} \end{matrix} \\ γ_{23} \end{matrix}}

(5)

For elastic isotropic FGMs,

C_{i j}

are given in terms of Young’s modulus

E (z)

and Poisson’s ratio

ν

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

(6)

If the thickness stretching effect is neglected (

ε_{33} = 0

), the constitutive relations of equation (5) can be rewritten as

{\begin{matrix} σ_{11} \\ σ_{22} \\ \begin{matrix} σ_{12} \\ \begin{matrix} σ_{13} \\ σ_{23} \end{matrix} \end{matrix} \end{matrix}} = [\begin{matrix} C_{11} & C_{12} & \begin{matrix} 0 & 0 & 0 \end{matrix} \\ C_{21} & C_{22} & \begin{matrix} 0 & 0 & 0 \end{matrix} \\ \begin{matrix} 0 \\ 0 \\ 0 \end{matrix} & \begin{matrix} 0 \\ 0 \\ 0 \end{matrix} & \begin{matrix} \begin{matrix} C_{44} \\ 0 \\ 0 \end{matrix} & \begin{matrix} 0 \\ C_{55} \\ 0 \end{matrix} & \begin{matrix} 0 \\ 0 \\ C_{66} \end{matrix} \end{matrix} \end{matrix}] {\begin{matrix} ε_{11} \\ ε_{22} \\ \begin{matrix} γ_{12} \\ \begin{matrix} γ_{13} \\ γ_{23} \end{matrix} \end{matrix} \end{matrix}}

(7)

where elastic constants are

\begin{array}{l} C_{11} (z) = C_{22} (z) = \frac{E (z)}{1 - ν^{2}} \\ C_{12} (z) = C_{21} (z) = \frac{ν E (z)}{1 - ν^{2}} C_{44} (z) = C_{55} (z) = C_{66} (z) = \frac{E (z)}{2 (1 + ν)} \end{array}

(8)

Governing equations

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

\int_{0}^{T} (δ U + δ W_{e} - δ K) d t = 0

(9)

where

δ U

is the variation of strain energy,

δ W_{e}

is the virtual work done by the applied external loads, and

δ K

is the variation of kinetic energy.

Variation of the plate strain energy can be expressed as

δ U = \int_{A} \int_{- H / 2}^{H / 2} (σ_{11} δ ε_{11} + σ_{22} δ ε_{22} + σ_{33} δ ε_{33} + σ_{12} δ γ_{12} + σ_{13} δ γ_{13} + σ_{23} δ γ_{23}) d A d z

(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} = \int_{A} - q (w_{b} + s (\frac{H}{2}) w_{s}) d A

(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} \int_{A} (P_{x} {(\frac{\partial u_{3}}{\partial x_{1}})}^{2} + P_{y} {(\frac{\partial u_{3}}{\partial x_{2}})}^{2}) d A

(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} \int_{A} \int_{- H / 2}^{H / 2} ρ (z) ({\dot{u_{1}}}^{2} + {\dot{u_{2}}}^{2} + {\dot{u_{3}}}^{2}) d A d z

(14)

where dot-superscript convention represents the partial differentiation w.r.t. time and

ρ (z)

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

δ 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)

δ 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} δ 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} (N_{12}^{g} + N_{21}^{g})}{\partial x_{1} \partial x_{2}} + P_{x} (\frac{\partial^{2} w_{b}}{\partial x_{1}^{2}} + \frac{\partial^{2} w_{s}}{\partial x_{1}^{2}}) \\ + P_{y} (\frac{\partial^{2} w_{b}}{\partial x_{2}^{2}} + \frac{\partial^{2} w_{s}}{\partial x_{2}^{2}}) - \frac{\partial}{\partial x_{1}} (Q_{1}^{00} + Q_{1}^{g}) - \frac{\partial}{\partial x_{2}} (Q_{2}^{00} + Q_{2}^{g}) - 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} δ 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} (N_{12}^{f} + N_{21}^{f})}{\partial x_{1} \partial x_{2}} + R^{s} + P_{x} (\frac{\partial^{2} w_{b}}{\partial x_{1}^{2}} + \frac{\partial^{2} w_{s}}{\partial x_{1}^{2}}) \\ + P_{y} (\frac{\partial^{2} w_{b}}{\partial x_{2}^{2}} + \frac{\partial^{2} w_{s}}{\partial x_{2}^{2}}) - \frac{\partial}{\partial x_{1}} (Q_{1}^{f} + Q_{1}^{s 0}) - \frac{\partial}{\partial x_{2}} (Q_{2}^{f} + Q_{2}^{s 0}) - s (\frac{H}{2}) 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} (I_{0}, I_{1}, I_{2}, I_{3}, I_{4}, I_{5}) = \int_{- \frac{H}{2}}^{\frac{H}{2}} ρ (z) (1, g_{u} (z), f_{u} (z), g_{v} (z), f_{v} (z), s (z)) d z \\ (J_{1}, J_{2}) = \int_{- \frac{H}{2}}^{\frac{H}{2}} ρ (z) (g_{u} (z) f_{u} (z), g_{v} (z) f_{v} (z)) d z (K_{1}, K_{2}, K_{3}, K_{4}, K_{5}) = \\ \int_{- H / 2}^{H / 2} ρ (z) ({(g_{u} (z))}^{2}, {(f_{u} (z))}^{2}, {(g_{v} (z))}^{2}, {(f_{v} (z))}^{2}, {(s (z))}^{2}) d z \end{array}

(19)

Stress resultants of the plate shown in equations (15–18) can be obtained by writing the expressions for

δ U

, using the strain–displacement relation equations equation (4)

(N_{1}^{0}, N_{1}^{g}, N_{1}^{f}) = \int_{- H / 2}^{H / 2} (1, g_{u} (z), f_{u} (z)) σ_{11} d z

(20)

(N_{2}^{0}, N_{2}^{g}, N_{2}^{f}) = \int_{- H / 2}^{H / 2} (1, g_{v} (z), f_{v} (z)) σ_{22} d z

R^{s} = \int_{- H / 2}^{H / 2} \frac{d s (z)}{d z} σ_{33} d z

(N_{12}^{0}, N_{12}^{g}, N_{12}^{f}, N_{21}^{g}, N_{21}^{f}) = \int_{- H / 2}^{H / 2} (1, g_{u} (z), f_{u} (z), g_{v} (z), f_{v} (z)) σ_{12} d z

(Q_{1}^{00}, Q_{1}^{s 0}, Q_{1}^{g}, Q_{1}^{f}) = \int_{- H / 2}^{H / 2} (1, s (z), \frac{d g_{u} (z)}{d z}, \frac{d f_{u} (z)}{d z}) σ_{13} d z

(Q_{2}^{00}, Q_{2}^{s 0}, Q_{2}^{g}, Q_{2}^{f}) = \int_{- H / 2}^{H / 2} (1, s (z), \frac{d g_{v} (z)}{d z}, \frac{d f_{v} (z)}{d z}) σ_{23} d z

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

(A_{11}^{0}, A_{11}^{g}, A_{11}^{f}) = \int_{- H / 2}^{H / 2} C_{11} (z) (1, g_{u} (z), f_{u} (z)) dz

(A_{22}^{0}, A_{22}^{g}, A_{22}^{f}) = \int_{- H / 2}^{H / 2} C_{22} (z) (1, g_{v} (z), f_{v} (z)) dz

[\begin{matrix} A_{12}^{0} & A_{12}^{g} & \begin{matrix} A_{12}^{f} & \begin{matrix} A_{12}^{g^{'}} & A_{12}^{f^{'}} \end{matrix} \end{matrix} \\ A_{44}^{0} & A_{44}^{g} & \begin{matrix} A_{44}^{f} & \begin{matrix} A_{44}^{g^{'}} & A_{44}^{f^{'}} \end{matrix} \end{matrix} \end{matrix}] = \int_{- H / 2}^{H / 2} {\begin{matrix} C_{12} (z) \\ C_{44} (z) \end{matrix}} (1, g_{u} (z), f_{u} (z), g_{v} (z), f_{v} (z)) dz

[\begin{matrix} D_{13}^{s} \\ D_{23}^{s} \end{matrix}] = \int_{- H / 2}^{H / 2} {\begin{matrix} C_{13} (z) \\ C_{23} (z) \end{matrix}} (\frac{ds (z)}{dz}) dz

(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

ε_{i j} = e_{i, 1} N_{1}^{0} + e_{i, 2} N_{2}^{0} + e_{i, 3} \frac{\partial^{2} w_{b}}{\partial x_{1}^{2}} + e_{i, 4} \frac{\partial^{2} w_{b}}{\partial x_{2}^{2}} + e_{i, 5} \frac{\partial^{2} w_{s}}{\partial x_{1}^{2}} + e_{i, 6} \frac{\partial^{2} w_{s}}{\partial x_{2}^{2}} + e_{i, 7} ω_{s} (i = j = 1,2)

γ_{12} = ε_{1} N_{12}^{0} + ε_{2} \frac{\partial^{2} w_{b}}{\partial x_{1} \partial x_{2}} + ε_{3} \frac{\partial^{2} w_{s}}{\partial x_{1} \partial x_{2}}

(22)

where coefficients

e_{i, j}

and

ε_{i}

can be presented as

\begin{array}{l} e_{1,1} = \frac{A_{22}^{0}}{A_{11}^{0} A_{22}^{0} - {A_{12}^{0}}^{2}}, e_{1,2} = - \frac{A_{12}^{0}}{A_{11}^{0} A_{22}^{0} - {A_{12}^{0}}^{2}} \\ e_{1,3} = \frac{A_{12}^{0} A_{12}^{g} - A_{22}^{0} A_{11}^{g}}{A_{11}^{0} A_{22}^{0} - {A_{12}^{0}}^{2}} + g_{u} (z), e_{1,4} = \frac{A_{12}^{0} A_{22}^{g} - A_{22}^{0} A_{12}^{g^{'}}}{A_{11}^{0} A_{22}^{0} - {A_{12}^{0}}^{2}} e_{1,5} = \frac{A_{12}^{0} A_{12}^{f} - A_{22}^{0} A_{11}^{f}}{A_{11}^{0} A_{22}^{0} - {A_{12}^{0}}^{2}} + f_{u} (z), \\ e_{1,6} = \frac{A_{12}^{0} A_{22}^{f} - A_{22}^{0} A_{12}^{f^{'}}}{A_{11}^{0} A_{22}^{0} - {A_{12}^{0}}^{2}} e_{1,7} = \frac{A_{12}^{0} D_{23}^{s} - A_{22}^{0} D_{13}^{s}}{A_{11}^{0} A_{22}^{0} - {A_{12}^{0}}^{2}}, e_{2,1} = - \frac{A_{12}^{0}}{A_{11}^{0} A_{22}^{0} - {A_{12}^{0}}^{2}} e_{2,2} = \frac{A_{11}^{0}}{A_{11}^{0} A_{22}^{0} - {A_{12}^{0}}^{2}}, \\ e_{2,3} = \frac{A_{12}^{0} A_{11}^{g} - A_{11}^{0} A_{12}^{g}}{A_{11}^{0} A_{22}^{0} - {A_{12}^{0}}^{2}} e_{2,4} = \frac{A_{12}^{0} A_{12}^{g^{'}} - A_{11}^{0} A_{22}^{g}}{A_{11}^{0} A_{22}^{0} - {A_{12}^{0}}^{2}} + g_{v} (z), e_{2,5} = \frac{A_{12}^{0} A_{11}^{f} - A_{11}^{0} A_{12}^{f}}{A_{11}^{0} A_{22}^{0} - {A_{12}^{0}}^{2}} e_{2,6} = \frac{A_{12}^{0} A_{12}^{f^{'}} - A_{11}^{0} A_{22}^{f}}{A_{11}^{0} A_{22}^{0} - {A_{12}^{0}}^{2}} + \\ f_{v} (z), e_{2,7} = \frac{A_{12}^{0} D_{13}^{s} - A_{11}^{0} D_{23}^{s}}{A_{11}^{0} A_{22}^{0} - {A_{12}^{0}}^{2}} ε_{1} = \frac{1}{A_{44}^{0}}, ε_{2} = g_{u} (z) + g_{v} (z) - \frac{A_{44}^{g} + A_{44}^{g^{'}}}{A_{44}^{0}} ε_{3} = f_{u} (z) + f_{v} (z) - \frac{A_{44}^{f} + A_{44}^{f^{'}}}{A_{44}^{0}} \end{array}

(23)

Solution procedure

The compatibility equation of an FG plate can be expressed as

\frac{\partial^{2} ε_{11}}{\partial x_{2}^{2}} + \frac{\partial^{2} ε_{22}}{\partial x_{1}^{2}} - \frac{\partial^{2} γ_{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_{x 0} a t x_{1} = 0, l_{1} \\ w_{b} = w_{s} = 0, N_{2}^{0} = N_{y 0} a t x_{2} = 0, l_{2} \end{array}

(27)

where

N_{x 0}

and

N_{y 0}

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} (x_{1}, x_{2}, t) = \sum_{m_{b} = 1}^{\infty} \sum_{n_{b} = 1}^{\infty} w b_{n_{b} m_{b}} (t) \sin (\frac{n_{b} π x_{1}}{l_{1}}) \sin (\frac{m_{b} π x_{2}}{l_{2}}) \\ w_{s} (x_{1}, x_{2}, t) = \sum_{m_{s} = 1}^{\infty} \sum_{n_{s} = 1}^{\infty} w s_{n_{s} m_{s}} (t) sin (\frac{n_{s} π x_{1}}{l_{1}}) sin (\frac{m_{s} π x_{2}}{l_{2}}) \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 (x_{1}, x_{2}, t) = \sum_{q = 1}^{\infty} \sum_{p = 1}^{\infty} f_{p q} (t) \sin (\frac{p π x_{1}}{l_{1}}) \sin (\frac{q π x_{2}}{l_{2}}) + \frac{1}{2} C_{x} x_{1}^{2} + \frac{1}{2} C_{y} x_{2}^{2}

(29)

where

f_{p q}

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_{p q} = - \frac{1}{e_{2,2} {(\frac{p π}{l_{1}})}^{4} + (e_{1,2} + ε_{1} + e_{2,1}) {(\frac{p π}{l_{1}})}^{2} {(\frac{q π}{l_{2}})}^{2} + e_{1,1} {(\frac{q π}{l_{2}})}^{4}} \\ {\begin{cases} (e_{2,3} {(\frac{n_{b} π}{l_{1}})}^{4} + (e_{1,3} - ε_{2} + e_{2,4}) {(\frac{n_{b} π}{l_{1}})}^{2} {(\frac{m_{b} π}{l_{2}})}^{2} + e_{1,4} {(\frac{m_{b} π}{l_{2}})}^{4}) \\ w b_{n_{b} m_{b}} δ (n_{b} - p) δ (m_{b} - q) + (e_{2,5} {(\frac{n_{s} π}{l_{1}})}^{4} + (e_{1,5} - ε_{3} + e_{2,6}) {(\frac{n_{s} π}{l_{1}})}^{2} {(\frac{m_{s} π}{l_{2}})}^{2} + e_{1,6} {(\frac{m_{s} π}{l_{2}})}^{4} - e_{2,7} {(\frac{n_{s} π}{l_{1}})}^{2} - e_{1,7} {(\frac{m_{s} π}{l_{2}})}^{2}) \\ w s_{n_{s} m_{s}} δ (n_{s} - p) δ (m_{s} - q) \end{cases}} \end{array}

(30)

In equation (3)

δ

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}

N_{1}^{0} = P_{x} a t x_{1} = 0, l_{1}, N_{2}^{0} = P_{y} a t 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}} \int_{0}^{l_{1}} (P_{y}) d x_{1} = P_{y} \\ C_{y} = \frac{1}{l_{2}} \int_{0}^{l_{2}} (P_{x}) 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

M {\begin{matrix} {\ddot{w b}}_{n_{b} m_{b}} \\ {\ddot{w s}}_{n_{s} m_{s}} \end{matrix}} + K {\begin{matrix} w b_{n_{b} m_{b}} \\ w s_{n_{s} m_{s}} \end{matrix}} = F

(33)

where

M

K

, and

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
Material property	Al, Metal layer	ZrO₂, Ceramic layer	Al, Metal layer	Al₂O₃, Ceramic layer
E (GPa)	70	200	70	380
$ρ$ (Kg/m³)	2702	5700	2702	3800
$ν$	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} (z) = - z, f_{u} (z) = - \frac{4 z^{3}}{3 H^{2}}, g_{v} (z) = - z, f_{v} (z) = - \frac{4 z^{3}}{3 H^{2}}, s (z) = 1

(34)

Hence

ε_{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 $\hat{ω} = ω_{1} H \sqrt{ρ_{b} / E_{b}}$ with $ω_{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 $\hat{ω} = ω_{1} H \sqrt{ρ_{b} / E_{b}}$ of square Al/ZrO₂ FG plate.

Method	$ε_{33}$	No. of Variables	$p$ = 0		$p$ = 1			$l_{1} / H$ = 5
Method	$ε_{33}$	No. of Variables	$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³	$\neq$ 0	—	0.4658	0.0578	0.2192	0.0596	0.0153	0.2197	0.2211	0.2225
Quasi-3D⁷	$\neq$ 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⁸	$\neq$ 0	9	—	—	0.2193 (0.046%)	—	—	0.2200 (0.136%)	0.2215 (0.181%)	0.2230 (0.225%)
Quasi-3D⁹	$\neq$ 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⁶	=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²³	=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_{c r}}{H E_{t}}

is used to obtain the results where

P_{c r}

is the critical buckling load of the FG plate.

Table 3.

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

$\frac{l_{1}}{H}$ .	$p$	$\tilde{N}$
$\frac{l_{1}}{H}$ .	$p$	Matsunaga,¹⁰ 29 Variables $ε_{33} \neq 0$ .	Present, RPT ( $w_{b}, w_{s}$ ), 2 Variables $ε_{33} = 0$ .	Present, CPT ( $w_{b}$ ), 1 Variable $ε_{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 (ω t) δ (x_{1} - \frac{l_{1}}{2}) δ (x_{2} - \frac{l_{2}}{2})$ with frequency of ω and magnitude of $f$ applied at the FG plate center with geometrical properties: $l_{1} = l_{2} = 0.2 (m), H = 0.01 (m), 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= ${\begin{matrix} w b_{n_{b} m_{b}} \\ w s_{n_{s} m_{s}} \end{matrix}}$ ), 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} (t) = y_{i} (0) \cos (ω_{i} t) + \frac{1}{ω_{i}} (\dot{y_{i}} (0) - \frac{f_{i}}{k_{i i}} \frac{ω}{1 - {(\frac{ω}{ω_{i}})}^{2}}) \sin (ω_{i} t) + \frac{f_{i}}{k_{i i}} \frac{\sin (ω t)}{1 - {(\frac{ω}{ω_{i}})}^{2}}

(35)

where

ω_{i}

is the i-th natural frequency of the system,

y_{i} (0)

and

\dot{y_{i}} (0)

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_{1,1}} = 0.1 \times H$ while other degrees of freedom are assumed to be zero. First natural frequency of the plate is $ω_{1} = 2 π \times 1281.53$ . A beating phenomenon can be seen occurring when $ω = 0.9 \times ω_{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 $(ω = 1 / 2 ω_{1})$ .

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 $ε_{33} \neq 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 = \sum_{n = 1}^{\infty} \sum_{m = 1}^{\infty} q_{0} sin (\frac{n π x_{1}}{l_{1}}) sin (\frac{m π x_{2}}{l_{2}})$ , 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} (z) = g_{v} (z) = - z

. The comparison is done for a square FG plate with Material Set-2 properties. Dimensionless transverse shear stress

{\bar{σ}}_{13} (z) = \frac{H}{q_{0} l_{1}} σ_{13} (0, l_{2} / 2, z)

and transverse displacement

{\bar{u}}_{3} (z) = \frac{10 E_{t} H^{3}}{q_{0} l_{1}^{4}} u_{3} (l_{1} / 2, l_{2} / 2, z)

are used to obtain the results.

Table 4.

Coefficient functions used in the displacement field.

$f_{u} (z) = f_{v} (z) = f (z)$ .	$s (z)$	$ε_{33}$
$- \frac{4 z^{3}}{3 H^{2}}$	1	$= 0$
$t g^{- 1} (\sin (\frac{π z}{H})) - z$	1	$= 0$
$Hsinh (\frac{z}{H}) - \frac{4 z^{3}}{3 H^{2}} cosh (\frac{1}{2})$	$\frac{1}{12} f' (z)$	$\neq 0$
$\frac{π}{H} z - \frac{9 π}{5 H^{3}} z^{3} + \frac{28 π}{25 H^{5}} z^{5}$	$\frac{1}{8} f' (z)$	$\neq 0$
$- 8 z + \frac{10}{H^{2}} z^{3} + \frac{6}{5 H^{4}} z^{5} + \frac{8}{7 H^{6}} z^{7}$	$\frac{3}{20} f' (z)$	$\neq 0$

Figure 5 and Figure 6 show variations of ${\bar{u}}_{3}$ and ${\bar{σ}}_{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 ${\bar{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 ${\bar{σ}}_{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 (z) = 1$ and thus $ε_{33} = 0$ , a straight line with no variation for ${\bar{u}}_{3}$ is observed. Although in quasi-3D models ( $ε_{33} \neq 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$ ), ${\bar{σ}}_{13}$ variation is parabolic along the thickness direction.

Figure 7 and Figure 8 shows variation of ${\bar{u}}_{3}$ and ${\bar{σ}}_{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 ${\bar{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 ${\bar{u}}_{3}$ and ${\bar{σ}}_{13}$ values are reduced, meaning that the stretching effect has the most influence on thick ones.

Figure 7.

Variation of ${\bar{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 ${\bar{σ}}_{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.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iDs

Leila Monajati

Navid Farid

Appendix

References

Mindlin

. Influence of rotary inertia and shear on flexural motions of isotropic, elastic plates. J Appl Mech 1951; 18: 31–38.

Rezaei

Saidi

Abrishamdari

, et al. Natural frequencies of functionally graded plates with porosities via a simple four variable plate theory: an analytical approach. Thin-Walled Struct 2017; 120: 366–377.

Vel

Batra

. Three-dimensional exact solution for the vibration of functionally graded rectangular plates. J Sound Vib 2004; 272: 703–730.

Sarangan

Singh

. Higher-order closed-form solution for the analysis of laminated composite and sandwich plates based on new shear deformation theories. Compos Struct 2016; 138: 391–403.

Mantari

Granados

Soares

, et al. Vibrational analysis of advanced composite plates resting on elastic foundation. Composites Part B 2014; 66: 407–419.

Ferreira

Batra

Roque

, et al. Natural frequencies of functionally graded plates by a meshless method. Compos Struct 2006; 75: 593–600.

Neves

Ferreira

Carrera

, et al. A quasi-3D hyperbolic shear deformation theory for the static and free vibration analysis of functionally graded plates. Compos Struct 2012; 94: 1814–1825.

Neves

Ferreira

Carrera

, et al. Static, free vibration and buckling analysis of isotropic and sandwich functionally graded plates using a quasi-3D higher-order shear deformation theory and a meshless technique. Composites Part B 2013; 44: 657–674.

Thai

Bui

, et al. A quasi-3D hyperbolic shear deformation theory for functionally graded plates. Acta Mech 2014; 225: 951–964.

10.

Matsunaga

. Free vibration and stability of functionally graded plates according to a 2-D higher-order deformation theory. Compos Struct 2008; 82: 499–512.

11.

Hellal

Bourada

Hebali

, et al. Dynamic and stability analysis of functionally graded material sandwich plates in hygro-thermal environment using a simple higher shear deformation theory. J Sandwich Struct Mater 2019. DOI: 10.1177/1099636219845841

12.

Yaghoobi

Taheri

. Characterization of the vibration, stability and static responses of graphene-reinforced sandwich plates under mechanical and thermal loadings using the refined shear deformation plate theory. J Sandwich Struct Mater 2021. DOI: 10.1177/1099636221993865

13.

Tahir

Chikh

Tounsi

, et al. Wave propagation analysis of a ceramic-metal functionally graded sandwich plate with different porosity distributions in a hygro-thermal environment. Compos Struct 2021; 269: 114030.

14.

Zaitoun

Chikh

Tounsi

, et al. Influence of the visco-Pasternak foundation parameters on the buckling behavior of a sandwich functional graded ceramic–metal plate in a hygrothermal environment. Thin-Walled Struct 2022; 170: 108549.

15.

Huffington

. Response of elastic columns to axial pulse loading. AIAA J 1963; 1: 2099–2104.

16.

Murty

. Flexure of composite plates. Compos Struct 1987; 7: 161–177.

17.

Senthilnathan

Lim

Lee

, et al. Buckling of shear-deformable plates. AIAA J 1987; 25: 1268–1271.

18.

Shimpi

. Refined plate theory and its variants. AIAA J 2002; 40(1):137–146.

19.

Thai

Choi

. A refined plate theory for functionally graded plates resting on elastic foundation. Compos Sci Technol 2011; 71: 1850–1858.

20.

Benachour

Tahar

Tounsi

, et al. A four variable refined plate theory for free vibrations of functionally graded plates with arbitrary gradient. Composites Part B 2011; 42: 1386–1394.

21.

Hadji

Atmane

Tounsi

, et al. Free vibration of functionally graded sandwich plates using four-variable refined plate theory. Appl Math Mech 2011; 32: 925–942.

22.

El Meiche

Tounsi

Ziane

, et al. A new hyperbolic shear deformation theory for buckling and vibration of functionally graded sandwich plate. Int J Mech Sci 2011; 53: 237–247.

23.

Thai

Choi

. Efficient higher-order shear deformation theories for bending and free vibration analyses of functionally graded plates. Arch Appl Mech 2013; 83: 1755–1771.

24.

Khosravifard

Hematiyan

, et al. A new refined simple TSDT-based effective meshfree method for analysis of through-thickness FG plates. Appl Math Modell 2018; 57: 514–534.

25.

Mahmoudi

Benyoucef

Tounsi

, et al. A refined quasi-3D shear deformation theory for thermo-mechanical behavior of functionally graded sandwich plates on elastic foundations. J Sandwich Struct Mater 2019; 21: 1906–1929.

26.

Rouzegar

Koohpeima

Abad

. Dynamic analysis of laminated composite plate integrated with a piezoelectric actuator using four-variable refined plate theory. Iran J Sci Technol Trans Mech Eng 2020; 44: 557–570.

27.

Nguyen

Abdel-Wahab

, et al. A refined quasi-3D isogeometric analysis for functionally graded microplates based on the modified couple stress theory. Comput Methods Appl Mech Eng 2017; 313: 904–940.

28.

Zenkour

. A simple four-unknown refined theory for bending analysis of functionally graded plates. Appl Math Model 2013; 37: 9041–9051.

29.

Nguyen

Hui

Lee

, et al. An efficient computational approach for size-dependent analysis of functionally graded nanoplates. Comput Methods Appl Mech Eng 2015; 297: 191–218.

30.

Nguyen-Xuan

Tran

Kulasegaram

, et al. Isogeometric analysis of functionally graded plates using a refined plate theory. Composites Part B 2014; 64: 222–234.

31.

Khorsand

. Dynamic analysis of a functionally graded piezoelectric spherical shell under mechanical and thermal shocks. Proc Inst Mech Eng Part C: J Mech Eng Sci 2014; 228: 632–645.

32.

Abdalla

Casagrande

Moro

, et al. Thermo-mechanical analysis and optimization of functionally graded rotating disks. J Strain Anal Eng Des 2020; 55: 159–171.

33.

Kurpa

Shmatko

. Buckling and free vibration analysis of functionally graded sandwich plates and shallow shells by the Ritz method and the R-functions theory. Proc Inst Mech Eng Part C: J Mech Eng Sci 2021; 235: 4582–4593.

34.

Huang

. Accurate solution for functionally graded beams with arbitrarily varying thicknesses resting on a two-parameter elastic foundation. J Strain Anal Eng Des 2020; 55: 222–236.

35.

Salehi Kolahi

Rahmani

Moeinkhah

, et al. Mechanical analysis of shrink-fitted thick FG cylinders based on first order shear deformation theory and FE simulation. Proc Inst Mech Eng Part C: J Mech Eng Sci 2021. DOI: 10.1177/09544062211012729

36.

Karimiasl

Ebrahimi

Akgöz

, et al. Buckling and post-buckling responses of smart doubly curved composite shallow shells embedded in SMA fiber under hygro-thermal loading. Compos Struct 2019; 223: 110988.

37.

Moallemi

Baghani

Almasi

, et al. Large deformation and stability analysis of functionally graded pressure vessels: an analytical and numerical study. Proc Inst Mech Eng Part C: J Mech Eng Sci 2018; 232: 3300–3314.