On the Solution of Dynamic Stability Problem of Functionally Graded Viscoelastic Plates with Different Initial Conditions in Viscoelastic Media

Abstract

The widespread use of structural elements consisting of functionally graded (FG) materials in advanced technologies has led to extensive research. Due to the difficulties encountered during modeling and problem solving, the number of studies on the dynamic behavior of structural elements made of FG viscoelastic materials is quite limited compared to the number examining FG elastic materials. This study is one of the first attempts to solve the dynamical problem by the mathematical modeling of functionally graded viscoelastic plates (FG-VE-Ps) and viscoelastic media together with different initial conditions. FG-VE-Ps on viscoelastic foundations (VE-Fs) are assumed to be under compressive edge load in the longitudinal direction. The governing equations for FG-VE-Ps on VE-Fs are derived using Boltzmann and Volterra concepts. The problem is reduced to the solution of integro-differential equation system using the Galerkin method. Then, by performing Laplace transforms, new analytical expressions for the time-dependent deflection function and critical time at different initial conditions are found. The loss of stability of FG-VE-Ps on VE-Fs is modeled to cover three time-varying ranges: the first is the range in which the deflection function decreases; the second is the transition interval; the third is the increase range of deflection function, which leads to the loss of stability. The time corresponding to the sharp increase of the deflection function is defined as the critical time, and is determined both theoretically and numerically. The results are compared with the results obtained by various methods to confirm their accuracy. Finally, the effects of VE-Fs, VE material properties, and FG profiles on the critical time behavior of plates are studied numerically.

1. Introduction

The widespread use of composite materials leads to the need to consider the viscoelastic properties specific to composites with polymers. General questions of the theory of viscoelasticity are considered in the fundamental works [1,2,3], including the constructions of hereditary relations between stresses and strains, the reciprocity of such relations, methods of analyzing the kernels of creep and relaxation, and methods for solving static and dynamical problems. The applied problems of the theory of viscoelasticity, including the construction of continuous and discrete models of homogeneous viscoelastic plates, are considered in [4,5,6,7,8,9,10,11,12,13,14,15].

A new class of advanced heterogeneous materials, which are called FGMs, have great potential in many areas of technology. The study of the behavior of FGM structural elements under various conditions is very important for the development of mechanical engineering and aircraft construction. The field of missile technology is also interested in intensive research of this class of problems [16,17]. In most studies, although FGMs are used as elastic materials, their time-dependent behavior should be modeled to fully reveal their potential. In recent years, there have been some remarkable initiatives for the creation and application of the viscoelasticity theory of FGMs [18,19]. Some attempts have been made to examine the static and dynamic behavior of FG-VE structural elements in the past few years [20,21,22,23,24,25].

An important problem in many branches of technology is to increase the reliability and extend the service life of composite structures that interact with viscoelastic foundations. Fundamental results in the theory of viscoelastic foundations are described in [26,27]. The review of the literature shows that there are mainly studies of the behavior of homogeneous viscoelastic structural elements on VE-Fs. Pouresmaeeli et al. [28] analyzed the vibration behavior of VE orthotropic nanoplates on VE-Fs. Karlicic et al. [29] studied the transverse vibration of a nonlocal VE-orthotropic multi-nanoplate system embedded in VE medium. Zhang et al. [30] presented a vibration analysis of VE single-walled carbon nanotubes resting on VE-Fs. Zamani et al. [31] examined the free vibration of thick VE-Ps on visco-Pasternak foundations (V-PFs) using higher-order theory. Zeighampour et al. [32] examined wave propagation in a VE thin nanoshell resting on V-PFs based on nonlocal strain gradient theory. Zenkour and Sobhy [33] studied the nonlocal piezo-hygrothermal analysis of the vibration characteristics of a piezoelectric Kelvin–Voigt VE nanoplate embedded in VE medium. Sophy and Zenkour [34] examined the modified couple stress model for the bending of normal deformable VE nanobeams resting on visco-Pasternak foundations. Sobhy and Radwan [35] studied the influence of a 2D magnetic field on the hygrothermal bending of sandwich CNT-reinforced microplates with VE core embedded in VE medium. Frahlia et al. [36] examined the vibration response of a functionally graded plate resting on a viscoelastic foundation. Alazwari and Zenkour [37] developed a quasi-3D refined theory for the linear vibration of functionally graded plates resting on visco-Winkler-Pasternak foundations.

The recent emergence of new-generation functionally graded viscoelastic materials and the introduction of their structural elements, which are often used in viscoelastic media, require the study of their dynamical problems. Zenkour et al. [38] presented a bending analysis of FG-VE beams resting on Pasternak elastic foundations. Shariyat and Alipor [39] analyzed the damping effect on the behavior of FG-VE-Ps on elastic foundations with variable thickness. Liu et al. [40] investigated the vibration of FG-VE magneto-electro porous nanobeams on V-PFs. Sofiyev et al. [41] examined the stability behavior of heterogeneous viscoelastic plates resting on two-parameter elastic foundations. Zenkour and El-Shahrany [42] presented the active control of a sandwich plate with reinforced magnetostrictive faces and a viscoelastic core resting on an elastic foundation. Li et al. [43] investigated the vibration behavior of a fiber-reinforced polymer composite plate with viscoelastic damping boundary conditions resting on a viscoelastic foundation. Yuan et al. [44] presented a magneto-hygro-thermal vibration analysis of viscoelastic nanobeams reinforced with carbon nanotubes resting on Kerr elastic foundations.

To the author’s knowledge, the dynamic buckling and vibrational behavior of functionally graded viscoelastic plates on a viscoelastic foundation have not yet been analytically reported due to difficulties in modeling and problem solving. The novelty of the present work is the proposal of an analytical solution for the dynamical problem under various initial conditions of a functionally graded plate and soil having viscoelastic property in contact with each other. The rest of this paper is structured as follows: in Section 2, the basic relations and equations are derived after modelling of the material and soil properties; in Section 3, the solution of the problem under various initial conditions is presented by applying the Galerkin and Laplace methods; Section 4 consists of comparisons and an original numerical analysis; Section 5 summarizes conclusions.

2. Mathematical Modeling of the Problem

2.1. Basic Relationships

An FG-VE plate with length $a$, width $b$ and thickness $h$ on a VE-F and the coordinate system in the midplane is shown in Figure 1. It is subjected to compressive edge load, $N$, in the $x$ direction and is resting on the VE-F. Displacement $w$ is in the positive direction of the $z$ axis, and is assumed to be quite small compared to the thickness. It is also assumed that the normal stresses acting on the midplane are negligible compared to other forces.

The effective FGM properties change exponentially depending on the thickness coordinate [16]:

$$E^f=E^m{e}^{\left(\frac{z+0.5h}{h}\right)\mathrm{l}\mathrm{n}\left(\frac{E^c}{E^m}\right)},\hspace{1em} {\nu }^f={\nu }^m{e}^{\left(\frac{z+0.5h}{h}\right)\mathrm{l}\mathrm{n}\left(\frac{{\nu }^c}{{\nu }^m}\right)},\hspace{1em} {\rho }^f={\rho }^m{e}^{\left(\frac{z+0.5h}{h}\right)\mathrm{l}\mathrm{n}\left(\frac{{\rho }^c}{{\rho }^m}\right)}$$

(1)

where the Young moduli of metal and ceramic surfaces of FGM are shown with $E^m$ and $E^c$, Poisson ratios with ${\nu }^m$ and ${\nu }^c$, and the densities with ${\rho }^m$ and ${\rho }^c$.

The viscoelastic foundations are modelled as [27,35]:

$$K(w,t)=K_{0}w+K_1\frac{\partial w}{\partial t}+\dots +K_{n-1}\frac{{\partial }^{n-1}w}{\partial t^{n-1}}+K_n\frac{{\partial }^{n}w}{\partial t^n}$$

(2)

in which $K_i,\hspace{1em} (i=0,1,\dots ,n)$ are constants for the VE-Fs and can be obtained from the experiments.

If the relaxation core of the integral operator is denoted as $R(t-\tau )$, the constitutive equations of FG-VE-Ps in the framework of the Boltzmann–Volterra principle are modeled as follows [1,5,41]:

$$\begin{array}{l}{\tau }_{11}=\frac{E^f}{1-{\left({\nu }^f\right)}^2}\left[{\epsilon }_{11}+{\nu }^f{\epsilon }_{22}-{\displaystyle \underset{0}{\overset{t}{\int }}R(t-\tau )\left({\epsilon }_{11}+{\nu }^f{\epsilon }_{22}\right)d\tau }\right],\\ {\tau }_{22}=\frac{E^f}{1-{\left({\nu }^f\right)}^2}\left[{\epsilon }_{22}+{\nu }^f{\epsilon }_{11}-{\displaystyle \underset{0}{\overset{t}{\int }}R(t-\tau )\left({\epsilon }_{22}+{\nu }^f{\epsilon }_{11}\right)d\tau }\right],\\ {\tau }_{12}=\frac{E^f}{2\left(1+{\nu }^f\right)}\left[{\gamma }_{12}-{\displaystyle \underset{0}{\overset{t}{\int }}R(t-\tau ){\gamma }_{12}d\tau }\right]\end{array}$$

(3)

The forces and moments are found from the following integrals [1]:

$$\left[\left(N_{11},N_{12},N_{22}\right),\left(M_{11},M_{12},M_{22}\right)\right]={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}(1,z)\left({\tau }_{11},{\tau }_{12},{\tau }_{22}\right)dz}$$

(4)

The relationships between Airy stress function ($F$) and forces are formed as follows [1,23,24]:

$$\left(N_{11},N_{12},N_{22}\right)=\left(hF{,}_{yy},\hspace{1em}-hF{,}_{xy},\hspace{1em}hF{,}_{xx}\right)$$

(5)

2.2. Governing Equations

Using the relations of (3)–(5), the force and moment components are expressed by the functions of $w$ and $F$, then these expressions and relation (2) are taken into account in the basic equations, the following system of integro-differential equations for dynamic stability and compatibility equations of FG-VE-Ps on VE-Fs is obtained:

$$\begin{array}{l}{L}_{11}F(x,y,t)+\left(L_{12}+L_{13}\right)w(x,y,t)-{\displaystyle \underset{0}{\overset{t}{\int }}\left[L_{11}F(x,y,\tau )+L_{12}w(x,y,\tau )\right]R(t-\tau )d\tau }=0\\ L_{21}F(x,y,t)+L_{22}w(x,y,t)-{\displaystyle \underset{0}{\overset{t}{\int }}\left[L_{21}F(x,y,\tau )+L_{22}w(x,y,\tau )\right]R(t-\tau )d\tau }=0\end{array}$$

(6)

where $L_{ij}(i=1,2,j=1,2,3)$ are expressed as:

$$\begin{array}{l}{L}_{11}=h{q}_2\frac{{\partial }^4}{\partial x^4}+2h\left(q_1-q_5\right) \frac{{\partial }^4}{\partial x^2\partial y^2}+h{q}_2\frac{{\partial }^4}{\partial y^4}\\ L_{12}=-q_3\frac{{\partial }^4}{\partial x^4}-2\left(q_4+q_6\right) \frac{{\partial }^4}{\partial x^2\partial y^2}-q_4 \frac{{\partial }^4}{\partial y^4}\\ L_{13}=\hspace{1em} -N\frac{{\partial }^2}{\partial x^2}-K_0-K_1\frac{\partial }{\partial t}-K_2\frac{{\partial }^2}{\partial t^2}-{\overline{\rho }}^f\frac{{\partial }^2}{\partial t^2}\\ L_{21}=h{p}_1\frac{{\partial }^4}{\partial x^4}+h\left(2p_2+p_5\right) \frac{{\partial }^4}{\partial x^2\partial y^2}+h{p}_1\frac{{\partial }^4}{\partial y^4}\\ L_{22}=-p_4\frac{{\partial }^4}{\partial x^4}-\left(2p_3-p_6\right) \frac{{\partial }^4}{\partial x^2\partial y^2}-p_4 \frac{{\partial }^4}{\partial y^4}\end{array}$$

(7)

in which

$$\begin{array}{l}{q}_1=e_1^1{p}_1+e_2^1{p}_2,\hspace{1em} q_2=e_1^1{p}_2+e_2^1{p}_1,\hspace{1em} q_3=e_1^1{p}_3+e_2^1{p}_4+e_1^2,\hspace{1em} q_4=e_1^1{p}_4+e_2^1{p}_3+e_2^2,\hspace{1em}{q}_5=e_6^1{p}_5,\hspace{1em} q_6=e_6^1{p}_6+2e_6^2,\hspace{1em}\\ p_1=\frac{e_1^0}{\eta },\hspace{1em} p_2=-\frac{e_2^0}{\eta },\hspace{1em} p_3=\frac{e_2^0{e}_1^1-e_1^1{e}_1^0}{\eta },\hspace{1em}{p}_4=\frac{e_2^0{e}_1^1-e_2^1{e}_1^0}{\eta },\hspace{1em} p_5=\frac{1}{e_6^0},p_6=-\frac{2e_6^1}{e_6^0},\hspace{1em} \eta ={\left(e_1^0\right)}^2-{\left(e_2^0\right)}^2,\hspace{1em}\\ e_1^k=\hspace{1em} {\displaystyle \underset{-h/2}{\overset{h/2}{\int }} }\frac{E^f}{1-{\left({\nu }^f\right)}^2}{z}^{k}dz,\hspace{1em} e_2^k=\hspace{1em}{\displaystyle \underset{-h/2}{\overset{h/2}{\int }} }\frac{{\nu }^f{E}^f}{1-{\left({\nu }^f\right)}^2}{z}^{k}dz,\hspace{1em}{e}_6^k=\hspace{1em} {\displaystyle \underset{-h/2}{\overset{h/2}{\int }} }\frac{E^f}{2\left(1+{\nu }^f\right)}{z}^{k}dz,\hspace{1em} k=0,\hspace{1em}1,\hspace{1em}2.\end{array}$$

(8)

It should be emphasized that in the viscoelastic foundation model, higher-order derivatives than the 2nd-order derivative have little effect on the behavior of the plate, so they can be neglected and those terms are not taken into account in the next steps. Therefore, $K(w,t)=K_{0}w+K_1\frac{\partial w}{\partial t}+K_2\frac{{\partial }^{2}w}{\partial t^2}$ is used in the basic equations instead of Equation (2).

3. Solution Procedure

Assuming that the edges of an FG-VE plate on a VE-F are simply supported, the solution of the basic equations is sought as follows [1,5,41]:

$$w=w_1(t)\sin \left(\overline{m}x\right)\cos \left(\overline{n}y\right),\hspace{1em} F=F_1(t)\sin \left(\overline{m}x\right)\cos \left(\overline{n}y\right)$$

(9)

where $w_1(t)$ and $F_1(t)$ are the functions depending on the time, $\left(\overline{m},\hspace{1em}\overline{n}\right)=\left(\frac{m}{a},\frac{n}{b}\right)\pi$, where $m$ and $n$ are the modes.

Substituting (9) into the set of Equation (6), performing the Galerkin procedure and removing $F_1(t)$ from the system, the following integro-differential equation is obtained:

$$\left({\overline{\rho }}^f+K_2\right)\frac{d^2{w}_1(t)}{d{t}^2}+K_1\frac{d{w}_1(t)}{dt}+\left[Q+K_0-N{\overline{m}}^2\right]w_1(t)-\hspace{1em}Q\hspace{1em}{\displaystyle \underset{0}{\overset{t}{\int }}R(t-\tau )w_1(\tau )d\tau }=0$$

(10)

where

$$\begin{array}{l}Q=q_3{\overline{m}}^4+2\left(q_4+q_6\right){\overline{m}}^2{\overline{n}}^2+q_3{\overline{n}}^4\hspace{1em}-\frac{\left[q_2{\overline{m}}^4+2\left(q_1-q_5\right){\overline{m}}^2{\overline{n}}^2+q_2{\overline{n}}^4\right]}{p_1{\overline{m}}^4+\left(2p_2+p_5\right){\overline{m}}^2{\overline{n}}^2+p_1{\overline{n}}^4}\left[p_4{\overline{m}}^4+\left(2p_3-p_6\right){\overline{m}}^2{\overline{n}}^2+p_4{\overline{n}}^4\right]\\ {\overline{\rho }}^f={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}{\rho }^{f}dz}\end{array}$$

(11)

The following initial conditions (ICs) are taken into account:

$$w_1(0)\ne 0,\hspace{1em}{w^{\prime }}_1(0)\ne 0,\hspace{1em} \mathrm{at}\hspace{1em} t=0$$

(12)

By using the Laplace transform in Equation (10) for the ICs in (12), the following equation is obtained in the images:

$$W_1(s_1)=\frac{(1+s_1)\left[(s_1+{\tilde{K}}_{11})w_{10}+w_{11}\right]}{s_1^3+(1+{\tilde{K}}_{11})s_1^2+r_1(1-p^{VF})s_1-r_1(p^{VF}-q^{VF})}$$

(13)

where $R(t-\tau )=\gamma \hspace{1em}{e}^{-\gamma (t-\tau )}\hspace{1em}(0<\gamma <1),\hspace{1em} W_1(s)={\displaystyle \underset{0}{\overset{\infty }{\int }}{w}_1(t)e^{-st}dt,\hspace{1em} }s=s_1\gamma \hspace{1em},$ in which $W_1(s)$ is the Laplace transformation of the original function $w_1(t)$, $\gamma$ denotes the viscoelasticity parameter, $s$ is the variable, and the following symbols are used [1]:

$$\begin{array}{l}{w}_{10}=\frac{w_1(0)}{\gamma },\hspace{1em} w_{11}=\frac{{w^{\prime }}_1(0)}{{\gamma }^2},\hspace{1em} r_1=\frac{r}{{\gamma }^2},\hspace{1em}{\tilde{K}}_{11}=\frac{{\tilde{K}}_1}{\gamma },\hspace{1em}{p}^{VF}=\frac{N}{N_{cr}^{VF}},\hspace{1em}{q}^{VF}=\frac{K_0}{{\overline{m}}^2{N}_{cr}^{VF}},\hspace{1em}\\ r=\frac{{\overline{m}}^2{N}_{cr}^{VF}}{{\overline{\rho }}^f+K_2}={\left({\omega }^{VF}\right)}^2,\hspace{1em}{\tilde{K}}_1=\frac{K_1}{K_2+{\overline{\rho }}^f},\hspace{1em}{N}_{cr}^{VF}=\frac{Q+K_0+\gamma K_1}{{\overline{m}}^2}\end{array}$$

(14)

Since $p^{VF}>q^{VF}$ and $r_1>>1$, one of the positive values that converts the denominator of Equation (13) to zero is indicated by $u_0^{VF}$, and the other two roots are easily found as follows:

$$u_{2,3}=-\frac{1+{\tilde{K}}_{11}+u_0^{VF}}{2}\pm i\Lambda$$

(15)

where

$$u_0^{VF}=\frac{r_1(p^{VF}-q^{VF})}{r_1(1-p^{VF})+u_0^{VF}(1+{\tilde{K}}_{11}+u_0^{VF})},\hspace{1em} \Lambda =\sqrt{r_2^{VF}-0.25{\left(1+{\tilde{K}}_{11}-u_0^{VF}\right)}^2}$$

(16)

in which

$$r_2^{VF}={\left(u_0^{VF}\right)}^2+r_1(1-p^{VF})\hspace{1em} \mathrm{and} r_2^{VF}-0.25{\left(1+{\tilde{K}}_{11}-u_0^{VF}\right)}^2>0$$

When the following symbols are replaced in Equation (13),

$$B_1=\frac{1+{\tilde{K}}_{11}+u_0^{VF}}{2},\hspace{1em} B_2=\frac{1+{\tilde{K}}_{11}-u_0^{VF}}{2},\hspace{1em}{B}_3={\tilde{K}}_{11},A=\frac{1+{\tilde{K}}_{11}+3u_0^{VF}}{2}\hspace{1em},s_2=s_1+\frac{1+{\tilde{K}}_{11}+u_0^{VF}}{2}$$

(17)

after some mathematical operations, it takes the following standard form:

$$\begin{array}{ll}{W}_1(s_2)& =\frac{\left[A(B_2-B_1)+A^2-B_1{B}_2+B_3\right]w_{10}+(A+B_2)w_{11}}{A^2+{\Lambda }^2}\frac{1}{s_2-A}\\ & +\frac{\left[{\Lambda }^2-A(B_2-B_1)+B_1{B}_2+B_3\right]w_{10}-(A+B_2)w_{11}}{A^2+{\Lambda }^2}\frac{s_2}{s_2^2+{\Lambda }^2}\\ & +\frac{\left[A+(B_2-B_1)\right]{\Lambda }^2{w}_{10}+A\left(B_1{B}_2-B_3\right)w_{10}-A({\Lambda }^2-A{B}_2)w_{11}}{A^2+{\Lambda }^2}\frac{\Lambda }{s_2^2+{\Lambda }^2}\frac{1}{\Lambda }\end{array}$$

(18)

By using the inverse Laplace transform, the original function of $W_1(s_2)$ from Equation (18) turns into the following form:

$${\overline{w}}_1(t_1)=\left(d_1+d_4\right)e^{u_0^{VF}{t}_1}+\left[(d_2+d_5)\cos \left(\Lambda t_1\right)+(d_3+d_6)\sin \left(\Lambda t_1\right)\right]e^{-\frac{1+{\tilde{K}}_{11}+u_0^{VF}}{2}{t}_1}$$

(19)

where

$$\begin{array}{l}{d}_1=\frac{\left[A(B_2-B_1)+A^2-B_1{B}_2+B_3\right]w_{10}}{A^2+{\Lambda }^2},\hspace{1em} d_2=\frac{\left[{\Lambda }^2-A(B_2-B_1)+B_1{B}_2+B_3\right]w_{10}}{A^2+{\Lambda }^2},\\ d_3=\frac{\left[A{\Lambda }^2+{\Lambda }^2(B_2-B_1)+A(B_1{B}_2-B_3\right]w_{10}}{A^2+{\Lambda }^2}\frac{1}{\Lambda },\hspace{1em} d_4=\frac{A+B_2}{A^2+{\Lambda }^2}{w}_{11},\\ \hspace{1em}{d}_5=-\frac{A+B_2}{A^2+{\Lambda }^2}{w}_{11},\hspace{1em} d_6=\frac{{\Lambda }^2-A{B}_2}{A^2+{\Lambda }^2}\frac{w_{11}}{\Lambda }\end{array}$$

(20)

The loss of stability of FG-VE-Ps is a time-dependent process that develops over time. This process covers three ranges: the first is the range of decreasing the deflection function; the second is the transition interval; the third is the range of the deflection function, which leads to the loss of stability. The time corresponding to a sharp increase of the deflection function is defined as the critical time, and, in this study, it is determined both theoretically and numerically.

(a)

In the particular case, as $w_1(0)\ne 0,\hspace{1em}{w^{\prime }}_1(0)=0$ or the first initial condition (FIC) and $r_2^{VF}>>1$ are satisfied, $d_4=d_5=d_6=0,\hspace{1em} d_3<<d_2$, (19) is transformed as:

$${\overline{w}}_{11}(t_1)=d_{11}{e}^{u_0^{VF}{t}_1}+d_{12}{e}^{-\frac{1+{\tilde{K}}_{11}+u_0^{VF}}{2}{t}_1}$$

(21)

where

$$d_{11}=\frac{A(B_2-B_1)+A^2-B_1{B}_2+B_3}{A^2+{\Lambda }^2},\hspace{1em} d_{12}=\frac{{\Lambda }^2-A(B_2-B_1)+B_1{B}_2+B_3}{A^2+{\Lambda }^2},\hspace{1em} {\overline{w}}_{11}(t_1)=\frac{w_1(t_1)}{w_{10}}$$

(22)

The critical time of FG-VE-Ps on VE-Fs is found from the solution of ${\overline{w^{\prime }}}_{11}(t_1)=0$ as follows:

$$t_{11cr}^{VF}=\frac{2}{1+{\tilde{K}}_{11}+3u_0^{VF}}\ln \left[\frac{1+{\tilde{K}}_{11}+u_0^{VF}}{2u_0^{VF}}\frac{d_{12}}{d_{11}}\right]$$

(23)

(b)

In the particular case, as $w_1(0)=0,\hspace{1em}{w^{\prime }}_1(0)\ne 0$ or second initial conditions (SICs) and $r_2^{VF}>>1$ are satisfied, $d_1=d_2=d_3=0,\hspace{1em} d_5<<d_6$, Equation (19) turns into the following form:

$${\overline{w}}_{12}(t_1)=d_{14}{e}^{u_0^{VF}{t}_1}+d_{16}{e}^{-\frac{1+u_0^{VF}}{2}{t}_1}$$

(24)

where

$$d_{14}=\frac{A+B_2}{A^2+{\Lambda }^2},\hspace{1em} d_{16}=\frac{{\Lambda }^2-A{B}_2}{A^2+{\Lambda }^2}\frac{1}{\Lambda },\hspace{1em} {\overline{w}}_{12}(t_1)=\frac{w_1(t_1)}{w_{11}}$$

(25)

The critical time of FG-VE-Ps on VE-Fs is defined as:

$$t_{12cr}^{VF}=\frac{2}{1+{\tilde{K}}_{11}+3u_0^{VF}}ln\left[\frac{1+{\tilde{K}}_{11}+u_0^{VF}}{2u_0^{VF}}\frac{d_{16}}{d_{14}}\right]$$

(26)

In particular, at $N=0$, Equation (19) is converted to the equation for the deflection function of the vibration of FG-VE-Ps on VE-Fs. Equations (23) and (26) are converted to the formulas for unconstrained FG-VE-Ps, with $K_i=0\hspace{1em} (i=0,1,2)$.

4. Numerical Analysis

To demonstrate the accuracy of the analytical expressions obtained in this study, our results are compared with the results of [4], which solves the free vibration problem of homogeneous isotropic VE shells without VE-Fs using the average method. In the current study, at $N=0$, the first equation of the system (10) becomes the motion equation, and Equation (19) becomes the vibration equation of the FG-VE-Ps. The comparison is made for the following initial conditions: $w_1(0)\ne 0,\hspace{1em}{w^{\prime }}_1(0)=0$ or FIC. For comparison, the obtained Equation (26) for the vibration of homogenous isotropic VE shells is converted to the formula for the vibration of homogeneous isotropic VE-Ps (taking into account that $b=\pi R$). In our study, the influences of the compressive load and VE-Fs are ignored, and calculations are made by converting FGM into the homogeneous material ($V_c=0$). The following data are used in the comparison:$a=40h,\hspace{1em} b=0.4a,\hspace{1em} E^m=2.05098\times {10}^5(\mathrm{Mpa}),\hspace{1em}{\rho }^m=8.9\times {10}^3\hspace{1em}\mathrm{kg}/{\mathrm{m}}^3,\hspace{1em}{\nu }^m=0.3262,\hspace{1em}n=4,w_{10}={10}^{-3}$ and $\gamma =0.05$. As can be seen from Figure 2, the results of our study are in harmony with the results of [4].

In the next step, the specific numerical analyses are presented for the CTs of FG-VE-Ps resting on VE-Fs. The FGM consists of a mixture of silicon nitride (Si₃N₄) and Ni (Ni), and its properties are as follows [16]: $E^c=3.22271\times {10}^5(\mathrm{Mpa}),\hspace{1em}{E}^m=2.05098\times {10}^5(\mathrm{Mpa}),$${\nu }^c=0.24,\hspace{1em} {\nu }^m=0.31,\hspace{1em} {\rho }^c=2370\hspace{1em}\mathrm{kg}/{\mathrm{m}}^3,\hspace{1em}{\rho }^m=8900\hspace{1em} \mathrm{kg}/{\mathrm{m}}^3$, where the above properties are calculated at room temperature, that is, for $T=300\hspace{1em}\mathrm{K}$ [16]. The viscoelastic parameter is $\gamma =1/2$, and initial conditions are $w_1(0)\ne 0,\hspace{1em} {w^{\prime }}_1(0)=0$ and $w_1(0)=0,\hspace{1em} {w^{\prime }}_1(0)\ne 0$. In the tables and graphics created for the analysis, the width-to-length and the height-to-thickness ratios of VE plates, as well as the soil coefficients, are varied. The numbers of transverse and longitudinal modes corresponding to the critical times are obtained as $\left(m, n\right)=\left(1, 1\right)$ in all tables and graphs.

In Figure 3 and Figure 4, the time-dependent change of ${\overline{w}}_1(t_1)$ is calculated numerically using the Maple program for Equation (19) to verify the accuracy of Equations (23) and (26), where the critical time values are obtained. In numerical calculations, the following data are used: $a=50h,\hspace{1em} a=2b,\hspace{1em} \gamma =0.5,\hspace{1em}{p}^{VF}=0.7,\hspace{1em}{K}_0=4.2\times {10}^7(\mathrm{N}/{\mathrm{m}}^3),$ $K_1=10\hspace{1em}(\mathrm{N}\times \sec /{\mathrm{m}}^3)$ and $K_2=5\hspace{1em}(\mathrm{N}\times {\sec }^2/{\mathrm{m}}^3)$ [1,3,4,5,27,33,41]. The units for the VE-Fs are also valid for the subsequent tables and figures. These calculations use FG-VE-Ps consisting of the mixture of materials presented above. The time-dependent variation of ${\overline{w}}_1(t_1)$ of FG-VE-Ps on VE-Fs according to Equation (19), and the critical time corresponding to the moment of increasing ${\overline{w}}_1(t_1)$ for the initial conditions $w_1(0)\ne 0,\hspace{1em} {w^{\prime }}_1(0)=0$ and $w_1(0)=0,\hspace{1em} {w^{\prime }}_1(0)\ne 0$, respectively, are presented in Figure 3 and Figure 4. The moment when ${\overline{w}}_1(t_1)$ begins to increase is defined as the critical time. When the abovementioned data are used in Equations (23) and (26), the following values are obtained for the critical time: $t_{11cr}^{VF}=3.6733$ and $t_{12cr}^{VF}=1.755$, respectively. In Figure 3 and Figure 4, the critical times corresponding to an increase in ${\overline{w}}_1(t_1)$ maintain the accuracy of the critical times obtained by Equations (23) and (26).

The changes of the critical times of FG-VE-Ps on VE-Fs depending on the viscoelasticity parameter $\gamma$ with FIC and SIC are tabulated in

Table 1. Changes of critical times of FG-VE-Ps on VE-Fs depending on the viscoelastic parameter $\gamma$ with FIC and SIC ($h/a=0.02,\hspace{1em}b=0.5a,\hspace{1em}{p}^{VF}=0.7$).

	${\mathit{t}}_{11\mathit{c}\mathit{r}}^{\mathit{V}\mathit{F}}$			${\mathit{t}}_{12\mathit{c}\mathit{r}}^{\mathit{V}\mathit{F}}$
	Ni	FG-exp.	Si3N4	Ni	FG-exp.	Si3N4
$\gamma$	$K_0=\hspace{1em}{K}_1=K_2=0$
0.1	4.236	4.434	4.669	2.031	2.13	2.248
0.3	3.686	3.884	4.12	1.757	1.856	1.973
0.5	3.431	3.629	3.864	1.629	1.728	1.845
0.7	3.263	3.461	3.696	1.545	1.644	1.761
$\gamma$	$K_0=4.2\times {10}^7,\hspace{1em} K_1=K_2=0$
0.1	4.6	4.742	4.93	2.207	2.279	2.374
0.3	4.011	4.16	4.354	1.912	1.988	2.086
0.5	3.737	3.889	4.086	1.775	1.853	1.952
0.7	3.556	3.711	3.91	1.685	1.764	1.864
$\gamma$	$K_0=4.2\times {10}^7,\hspace{1em} K_1=10,\hspace{1em} K_2=0$
0.1	3.971	3.705	3.136	1.918	1.797	1.529
0.3	3.807	3.8	3.648	1.821	1.825	1.762
0.5	3.62	3.679	3.659	1.724	1.759	1.758
0.7	3.476	3.565	3.607	1.65	1.699	1.728
$\gamma$	$K_0=4.2\times {10}^7,\hspace{1em} K_1=10,\hspace{1em} K_2=5$
0.1	3.987	3.759	3.313	1.925	1.822	1.613
0.3	3.803	3.806	3.71	1.819	1.827	1.789
0.5	3.611	3.673	3.68	1.719	1.755	1.766
0.7	3.465	3.553	3.609	1.644	1.692	1.727

. Additional data required for the numerical analysis are presented in Table 1. As can be seen, the critical time values decrease with and without VE-Fs with increasing $\gamma$. It is noted that the critical time has the largest value when $K_0\ne 0$, while it has the smallest value when $K_2=0$. In addition, if the case of $K_i\ne 0(i=0,1,2)$ is compared with the case of $K_2=0$, although the critical time values increase, these values are less than those of unconstrained FG-VE-Ps for the small values of $\gamma$. The greatest influences of VE-Fs occur when $K_2=0$, and these effects are −16.4% and −15.22% for FIC and SIC, respectively, at $\gamma =0.1$. When the VE-Fs’ effects are compared to the critical times of FG-VE-Ps for the first and second initial conditions, the effect of the soil is great in the SIC for $K_0\ne 0$ only, while the other effects are higher in the FCI. When the CTs of an FG-VE-P and Ni-VE-P are compared, the influence of the FG-exp. profile on the critical time parameter increases in both initial conditions due to the increase of $\gamma$ in the groundless state. When $K_0\ne 0$ and other VE-F coefficients are zero, although the influence of the FG-exp. profile decreases to the critical time, the increase continues due to the increase of $\gamma$. In both boundary conditions, the influence of the FG profile on the critical times for $K_0\ne 0$ and $K_1\ne 0$ varies irregularly depending on $\gamma$, while it is considerable. Compared with FG-VE-P and Si₃N₄-VE-P, the greatest effect of the FG profile occurs at $K_2=0$, and is approximately +18% for both initial conditions.

The changes of the critical times of FG-VE-Ps on VE-Fs depending on $K_1$ with FIC and SIC for different $K_2$ are tabulated in Figure 5. The data required for numerical analysis are: $a=50h, a=2b, p^{VF}=0.7 ,\gamma =0.5, K={10}^8(\mathrm{N}/{\mathrm{m}}^3)$. In both initial conditions, the critical time values increase weakly with increasing $K_1$, while the critical time values decrease weakly as $K_1$ increases. In addition, it is observed that the influences of the VE-F on the critical time increase as $K_1$ and $K_2$ increase, the effect of the FG profile on the critical time increases weakly as $K_1$ increases, and this effect decreases with the increase of $K_2$.

The changes of critical times of FG-VE-Ps on VE-Fs (with various $K_i(i=0,1,2)$) against $a/h$ for FIC and SIC are shown in Figure 6. The critical times of FG-VE-Ps for various FIC and SIC decrease as $a/h$ increases. The data for the numerical analysis are: $a/b=2, p^{VF}=0.7, \gamma =0.5$ and $K_i(i=0,1,2)$, as shown in Figure 6. With increasing $a/h$ ratio, the influence of the VE-F (for $K_0=4.5\times {10}^7, K_1=K_2=0$) on the values of the critical time increases, while the effect of the VE-F (for $K_0=4.5\times {10}^7, K_1=40, K_2=0$) on the critical time is significantly greater than the effect of the previous VE-F type, although the influence of the VE-F is quickly decreasing. In addition, it is found that the influence of the VE-F (for $K_0=4.5\times {10}^7, K_1=40 ,K_2=35$) on the critical time is smaller than the previous VE-F, and this effect also decreases with increasing $a/h$. Note that the effects of the FG profile on the critical time decrease for all soil types as $a/h$ increases.

The distributions of critical times of FG-VE-Ps on VE-Fs versus $a/b$ for FIC and SIC are plotted in Figure 7. The data required for the numerical analysis are: $a/h=50, p^{VF}=0.7, \gamma =0.5$ and $K_0=2.5\times {10}^6, K_1=40, K_2=30$. The influence of VE-F (for $K_0=2.5\times {10}^6, K_1=K_2=0$) on the critical time decreases with increasing $a/b$, whereas the influence of the VE-F (for $K_0=2.5\times {10}^6, K_1=40, K_2=0$) on the critical time is significantly greater than the influence of the previous VE-F type, and the influence of the VE-F increases as $a/b$ increases. In addition, it is found that the influence of the VE-F (for $K_0=2.5\times {10}^6, K_1=40, K_2=30$) on the critical time is small compared to the previous VE-F, and this effect increases with increasing $a/b$. The influences of the FG profile on the critical time decrease for $K_0=2.5\times {10}^6, K_1=K_2=0$, while they increase for $K_0=2.5\times {10}^6, K_1=40, K_2=0$ and $K_0=2.5\times {10}^6,$ $K_1=40, K_2=30$ as $a/b$ increases.

5. Conclusions

The dynamic stability of FG-VE-Ps with different initial conditions under compressive edge load resting on VE-Fs was solved. The basic equations for FG-VE-Ps on VE-Fs were derived using the concepts of Boltzmann and Volterra. The problem was reduced to the solution of an integro-differential equation system using the Galerkin method. By using Laplace transforms the new analytical expressions were obtained for deflection function and critical time with two initial conditions. The results were compared with the results obtained by various methods, and the accuracy of the results was confirmed. Finally, the effects of the VE-Fs, VE properties and the FG exponential profile on the dynamic characteristics of plates with various initial conditions were investigated numerically.