**1. Introduction**

Circular cylindrical shells have an essential role in various fields of engineering applications such as aircraft, pressure vessels, gas turbines, and many other industrial purposes because of their excellent performance. Due to the advancement of the knowledge and technology, in recent years a new category of materials with interesting properties, named as functionally graded materials (FGMs) has been successfully applied. Conventional types of these materials are made of two or more different constituent phases, namely, the ceramic and metal phases, which are distributed gradually according some fixed functions. Since FGMs have some extraordinary properties, namely, a high temperature and a corrosion resistance, as well as an improved residual stress distribution, they are widely studied in many field of the applied sciences and they are adopted as structural components in military, medical, or aerospace industries, as well as in power plants or vessels. Thus, due to their special privileges in comparison with traditional materials, most industries make effort to exert such materials in lieu of ordinary ones [1–4].

A large number of studies in literature has focused on the thermo-mechanical and buckling behavior of FGMs for shell and plate structures. In this context, only some research works associated with FG cylindrical shells will be reviewed here, in line with the perspective developed in the present work. Du et al. [5] investigated the nonlinear forced vibration response of FG cylindrical thin shells, and used the perturbation method and the numerical Poincaré maps to solve the governing equations of the problem. Rahimi et al. [6] studied the vibration of FG cylindrical shells with ring supports. It was found that symmetric and asymmetric boundary conditions affect significantly the vibrations of the structure, with a general increase or decrease, respectively. In a recent work, Ghasemi et al. [7] have studied the agglomeration effect of FG hybrid single-walled carbon nanotubes on the vibration of hybrid laminated cylindrical shell structures. Bich et al. [8] performed the nonlinear static and dynamic buckling of imperfect eccentrically stiffened FG thin circular cylindrical shells subjected to an axial compression. Beni et al. [9] presented a novel formulation based on a modified couple stress theory to study simply supported FG circular cylindrical shells in the framework of thin shell structures, whereby the vibration behavior based on a classical continuum was found to be quite unaffected by the length scale parameters. Da Silva et al. [10] studied the nonlinear vibrations of a simply supported fluid-filled FG cylindrical shell subjected to a lateral time-dependent load and an axial static preloading condition. Bich and Nguyen [11] applied the displacement method to study the nonlinear vibration of FG circular cylindrical shells subjected to an axial and transverse mechanical loading. Ghannad et al. [12] introduced an analytical solution for the deformation and stress response of axisymmetric clamped–clamped thick FG cylindrical shells with variable thickness, while applying the first-order shear deformation theory and the perturbation theory, based on the Donnell's nonlinear large deflection theory. To date, many analytical and numerical approaches have been proposed in literature to handle simple and coupled vibration problems of cylindrical shell structures, including thermo-elastic, piezoelectric, and thermo-piezoelectric multi-field problems (see refs. [13–27], among others). As far as FGMs are concerned, many recent studies about the free vibration and buckling response of conventional and bi-directional FG cylindrical shells have been recently performed in literature [28–33]. A key point of the static and dynamic response of FG shell structures is related to the presence of porosities, which can form during a fabrication process, with possible effects on the global structural response. Indeed, an increasing attention to this aspect has been devoted in the scientific community for a correct interpretation of the mechanical performances of FG materials and structures. For example, in a recent work, Kiran and Kattimani [34] assessed the possible effect of porosity on the vibration behavior and static response of FG magneto-electro-elastic plates, with a clear reduction of the natural frequencies for an increased porosity within the material. In another study, Kiran et al. [35] analyzed the effect of porosity on the structural behavior of skew FG magneto-electro-elastic plates. Barati et al. [36] performed the buckling analysis of higher order graded smart piezoelectric plates with porosities resting on an elastic foundation. It was found that the buckling behavior of piezoelectric plates is significantly influenced by the porosity distribution. In the further works by Wang et al. [37,38], the authors studied the vibration response of longitudinally traveling FG plates with porosities [37] while considering the thermo-mechanical coupled response in [38]. A similar free vibration problem was studied in [39,40] for FG cylindrical shells, by means of the sinusoidal shear deformation theory and the Rayleigh–Ritz method, accounting for the possible presence of defects and porosities. In the context of nanomaterials and nanostructures, some modified couple stress theories have been recently proposed as efficient theoretical tools to study their coupled thermomechanical vibration behavior, also in presence of different levels of porosity, see [41–47], among others.

Up to date, however, there is a general lack of works in literature focusing on the dynamic buckling response of bi-directional (BD)-FG cylindrical shell embedded in a Winkler–Pasternak foundation, including the simultaneous effect of porosity. To this end, we propose the third-order shear deformation theory (TSDT) to model the cylindrical shells with porosities, subjected to an axial compressive excitation. The Hamilton's principle will be employed to determine the governing partial equations of motion, whereby the generalized differential quadrature (GDQ) method is adopted to solve the problem together with the associated boundary conditions into a system of Mathieu–Hill equations. Afterward, the Bolotin method is employed to determine the boundaries of the dynamic instability region (DIR) of BD-FG cylindrical shells. A systematic study focuses on the sensitivity of the dynamic stability behavior to different dimensionless ratios, i.e., the thickness-to-radius ratio, or the length-to-radius ratio, as well as to different boundary conditions, transverse and longitudinal power indexes, porosity volume fractions and foundation constants. The paper is organized as follows. In Section 2 we determine the governing equations of the problem for porous BD-FG cylindrical shells, which are solved numerically according to the GDQ and Bolotin methods, as detailed in Section 3. Section 4 aims at validating the proposed approach and shows the main results from a broad numerical investigation, whereas the final remarks are discussed in Section 5.

#### **2. Governing Equations of the Problem**

Let us consider a BD-FG porous cylindrical shell embedded in an elastic foundation with thickness *h*, radius *R*, and length *L*, where two different porosity distributions of the constituent phases are accounted for the analysis, namely an even and an uneven distribution, see Figures 1 and 2. We assume a BD-FG material made of a metal (labeled as *m*) and a ceramic (labeled as *c*) in the inner and outer shell surfaces, respectively. While the material properties for a conventional FG model vary continuously along the thickness direction from a ceramic or metal to another one, the basic material properties selected herein, vary also along the shell length from the metal to the ceramic phase. To this end, the volume fractions of the ceramic and metal phases are defined as follows

$$V\_{\varepsilon}(\mathbf{x}, z) = \left(\frac{1}{2} + \frac{z}{h}\right)^{n\_{\varepsilon}} \left(\frac{\mathbf{x}}{L}\right)^{n\_{\varepsilon}},\\V\_{m}(\mathbf{x}, z) = 1 - V\_{\varepsilon}(\mathbf{x}, z),\tag{1}$$

where *nz* and *nx* refer to the non-negative volume fraction exponents defining the profile variation of the material properties along the shell thickness and length directions, respectively. In addition, *z* and *x* stand for the radial distance from the mid-plane and longitudinal distance from the origin of the BD-FG cylindrical shell, respectively. The effective material properties (i.e., Yong's modulus, density and Poisson's ratio) of the BD-FG porous cylindrical shell are assumed to change according to a modified power law model with a linear algebraic combination of volume fractions of two basic materials. Two types of BD-FG material models include both even and/or uneven porosity distributions, i.e.,

$$E(\mathbf{x}, z) = E\_m + (E\_c - E\_m) \left(\frac{1}{2} + \frac{z}{h}\right)^{n\_x} \left(\frac{\mathbf{x}}{L}\right)^{n\_x} - \frac{\zeta}{2} (E\_c + E\_m),\tag{2a}$$

$$
\rho(\mathbf{x}, z) = \rho\_{\rm m} + (\rho\_c - \rho\_{\rm m}) \left( \frac{1}{2} + \frac{z}{h} \right)^{n\_c} \left( \frac{\mathbf{x}}{L} \right)^{n\_x} - \frac{\zeta}{2} (\rho\_c + \rho\_{\rm m}), \tag{2b}
$$

$$\nu(\mathbf{x}, z) = \nu\_{\rm m} + (\nu\_{\rm c} - \nu\_{\rm m}) \left(\frac{1}{2} + \frac{z}{h}\right)^{n\_{\rm c}} \left(\frac{\mathbf{x}}{L}\right)^{n\_{\rm c}} - \frac{\zeta}{2} (\nu\_{\rm c} + \nu\_{\rm m})\_{\rm \epsilon} \tag{2c}$$

for even distributions, or

$$E(\mathbf{x}, z) = E\_m + (E\_c - E\_m) \left(\frac{1}{2} + \frac{z}{h}\right)^{n\_x} \left(\frac{x}{L}\right)^{n\_x} - \frac{\xi}{2} (E\_c + E\_m) \left(1 + \frac{2|z|}{h}\right) \tag{3a}$$

$$
\rho(\mathbf{x}, z) = \rho\_m + (\rho\_c - \rho\_m) \left(\frac{1}{2} + \frac{z}{h}\right)^{n\_z} \left(\frac{x}{L}\right)^{n\_x} - \frac{\xi}{2} (\rho\_c + \rho\_m) \left(1 + \frac{2|z|}{h}\right) \tag{3b}
$$

$$\nu(\mathbf{x}, z) = \nu\_m + (\nu\_c - \nu\_m) \left(\frac{1}{2} + \frac{z}{h}\right)^{n\_\varepsilon} \left(\frac{x}{L}\right)^{n\_x} - \frac{\xi}{2} (\nu\_c + \nu\_m) \left(1 + \frac{2|z|}{h}\right) \tag{3c}$$

for an uneven distribution. In the all the expression (2) and (3), ζ and ξ denote the volume fraction of an even or uneven porosity inside the phases, respectively. While an even model accounts for porosities evenly distributed across the radial direction, the porosities in an uneven model is mostly concentrated in the shell mid-plane. It is worth noting that the uneven porosity distribution is linearly reduced from a larger value at mid-plane to a smaller value at the top and bottom sides of the structure.

**Figure 1.** Geometrical scheme of a bi-directional (BD)-functionally graded (FG) porous cylindrical shell embedded in an elastic foundation.

**Figure 2.** Distribution of porosity: (**a**) even porosity, (**b**) uneven porosity.

In what follows, we apply the TSDT, such that the displacement field of an arbitrary point of the shell along the *x*, *y*, and *z* axes, is defined as

$$\begin{aligned} u(\mathbf{x}, y, z, t) &= u\_0(\mathbf{x}, y, t) + zq \mathbf{e}\_x(\mathbf{x}, y, t) - \tilde{c}z^3 \Big| q\_\mathbf{x}(\mathbf{x}, y, t) + \frac{\partial \mathbf{w}\_0(\mathbf{x}, y, t)}{\partial \mathbf{x}} \Big| \\ v(\mathbf{x}, y, z, t) &= v\_0(\mathbf{x}, y, t) + zq\_\mathbf{y}(\mathbf{x}, y, t) - \tilde{c}z^3 \Big| q\_\mathbf{y}(\mathbf{x}, y, t) + \frac{\partial \mathbf{w}\_0(\mathbf{x}, y, t)}{\partial \mathbf{y}} \Big| \\ w(\mathbf{x}, y, z, t) &= w\_0(\mathbf{x}, y, t), \end{aligned} \tag{4}$$

where *c* = <sup>4</sup> <sup>3</sup>*h*<sup>2</sup> , and *u*, *v*, *w* stand for the longitudinal, circumferential, and transverse (radial) displacement components, respectively. These are determined, in turn, by means of the kinematic quantities *u*0, *v*0, and *w*<sup>0</sup> at the middle surface, and the rotations ϕ*<sup>x</sup>* and ϕ*<sup>y</sup>* of a transverse normal section about the *x* and *y* axis, respectively.

According to the TSDT, the strain components of the cylindrical shell can be written as [48]

$$\varepsilon\_{\rm xx} = \frac{\partial u\_0(\mathbf{x}, y, t)}{\partial \mathbf{x}} + z \frac{\partial \varrho\_{\mathbf{x}}(\mathbf{x}, y, t)}{\partial \mathbf{x}} - \tilde{c}z^3 \Big( \frac{\partial \varrho\_{\mathbf{x}}(\mathbf{x}, y, t)}{\partial \mathbf{x}} + \frac{\partial^2 w\_0(\mathbf{x}, y, t)}{\partial \mathbf{x}^2} \Big) \tag{5a}$$

*Appl. Sci.* **2020**, *10*, 1345

$$\varepsilon\_{yy} = \frac{\partial v\_0(\mathbf{x}, y, t)}{\partial y} + z \frac{\partial \wp\_y(\mathbf{x}, y, t)}{\partial y} - \mathfrak{L}z^3 \left(\frac{\partial \wp\_y(\mathbf{x}, y, t)}{\partial y} + \frac{\partial^2 w\_0(\mathbf{x}, y, t)}{\partial y^2}\right) + \frac{w}{R},\tag{5b}$$

$$\begin{split} \gamma\_{xy} &= \frac{\partial u\_0(\mathbf{x}, \mathbf{y}, t)}{\partial y} + \frac{\partial v\_0(\mathbf{x}, \mathbf{y}, t)}{\partial \mathbf{x}} + z \Big( \frac{\partial p\_x(\mathbf{x}, \mathbf{y}, t)}{\partial y} + \frac{\partial p\_y(\mathbf{x}, \mathbf{y}, t)}{\partial \mathbf{x}} \Big) \\ &- \overline{c} z^3 \Big( \frac{\partial p\_x(\mathbf{x}, \mathbf{y}, t)}{\partial y} + \frac{\partial p\_y(\mathbf{x}, \mathbf{y}, t)}{\partial \mathbf{x}} + 2 \frac{\partial^2 w\_0(\mathbf{x}, \mathbf{y}, t)}{\partial \mathbf{y} \partial \mathbf{x}} \Big), \end{split} \tag{5c}$$

$$
\varphi\_{\rm xz} = \varphi\_x(\mathbf{x}, y, t) - 2\overline{\mathbf{c}}z^2 \Big( \varphi\_x(\mathbf{x}, y, t) + \frac{\partial w\_0(\mathbf{x}, y, t)}{\partial \mathbf{x}} \Big) + \frac{\partial w\_0(\mathbf{x}, y, t)}{\partial \mathbf{x}},\tag{5d}
$$

$$
\gamma\_{yz} = q\rho\_y(\mathbf{x}, y, t) - 2\tilde{\mathbf{z}}^2 \Big( q\rho\_y(\mathbf{x}, y, t) + \frac{\partial w\_0(\mathbf{x}, y, t)}{\partial y} \Big) + \frac{\partial w\_0(\mathbf{x}, y, t)}{\partial y},\tag{5e}
$$

which are related to the stress components as follows

⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎪⎪⎪⎩

$$
\begin{pmatrix}
\sigma\_{xx} \\
\sigma\_{yy} \\
\sigma\_{xz} \\
\sigma\_{yz} \\
\sigma\_{xy}
\end{pmatrix} = \begin{bmatrix}
Q\_{11} & Q\_{12} & 0 & 0 & 0 \\
Q\_{21} & Q\_{22} & 0 & 0 & 0 \\
0 & 0 & Q\_{44} & 0 & 0 \\
0 & 0 & 0 & Q\_{55} & 0 \\
0 & 0 & 0 & 0 & Q\_{66}
\end{bmatrix} \begin{pmatrix}
\varepsilon\_{xx} \\
\varepsilon\_{yy} \\
\varepsilon\_{xz} \\
\varepsilon\_{yz} \\
\varepsilon\_{xy}
\end{pmatrix} \tag{6}
$$

where the stiffness *Qij* is defined as

$$\begin{aligned} Q\_{11}(\mathbf{x}, \mathbf{z}) = Q\_{22}(\mathbf{x}, \mathbf{z}) &= \frac{\operatorname{E}(\mathbf{x}, \mathbf{z})}{1 - (\nu(\mathbf{x}, \mathbf{z}))^2}, \\ Q\_{44}(\mathbf{x}, \mathbf{z}) = Q\_{55}(\mathbf{x}, \mathbf{z}) &= Q\_{66}(\mathbf{x}, \mathbf{z}) = \frac{\operatorname{E}(\mathbf{x}, \mathbf{z})}{2(1 + \nu(\mathbf{x}, \mathbf{z}))}. \end{aligned} \tag{7}$$

The strain energy of the BD-FG porous cylindrical shell is expressed as follow

$$\Pi\_S = \frac{1}{2} \int\_0^L \int\_0^{2\pi R} \int\_0^{\frac{1}{2}} (\sigma\_{xx}\varepsilon\_{xx} + \sigma\_{yy}\varepsilon\_{yy} + \tau\_{xy}\gamma\_{xy} + \tau\_{xz}\gamma\_{xz} + \tau\_{yz}\gamma\_{yz}) dzdydx,\tag{8}$$

and the kinetic energy of the cylindrical shell reads

$$\Pi\_T = \frac{1}{2} \int\_0^L \int\_0^{2\pi R} \int\_{-\frac{1}{2}}^{\frac{1}{2}} \rho(\mathbf{x}, z) \left( \left( \frac{\partial u(\mathbf{x}, y, z, t)}{\partial t} \right)^2 + \left( \frac{\partial v(\mathbf{x}, y, z, t)}{\partial t} \right)^2 + \left( \frac{\partial w(\mathbf{x}, y, z, t)}{\partial t} \right)^2 \right) dz dy d\mathbf{x}. \tag{9}$$

The total potential energy corresponding to the axial compressive load *Fa*(*t*) together with the Winkler and/or Pasternak elastic foundation, can be written as follow

$$\Pi\_{\rm E} = \frac{1}{2} \int\_{0}^{L} \int\_{0}^{2\pi R} \int\_{0}^{\frac{\hbar}{2}} \left( k\_{w} w\_{0}^{\prime} + k\_{\rm g} \left( \left( \frac{\partial w\_{0}}{\partial \mathbf{x}} \right)^{2} + \left( \frac{\partial w\_{0}}{\partial y} \right)^{2} \right) + F\_{\rm a}(t) \left( \frac{\partial w\_{0}}{\partial \mathbf{x}} \right)^{2} \right) dz dydx,\tag{10}$$

where *kw* and *kg* refer to the Winkler foundation stiffness and shear layer stiffness of the elastic foundation, respectively.

Consistently with the Hamilton's principle, the following governing equations of motion for the BD-FG cylindrical shell are determined

$$\int\_{t\_1}^{t\_2} (\delta \Pi\_T - \delta \Pi\_S - \delta \Pi\_E) dt = 0,\tag{11}$$

where the symbol δ denotes the variation of the energy quantities. By combination of Equations (8)–(10) and Equation (11), after integration by parts, we get the following governing equations of motion, under the assumption of a null value for *u*0, *v*0, *w*0, ϕ*x*, and ϕ*y*.

$$\frac{\partial \mathcal{N}\_{\text{xx}}}{\partial \mathbf{x}} + \frac{\partial \mathcal{N}\_{\text{xy}}}{\partial \mathbf{y}} = l\_0(\mathbf{x}) \frac{\partial^2 u\_0(\mathbf{x}, \mathbf{y}, t)}{\partial t^2} + (l\_0(\mathbf{x}) - \overline{\epsilon} l\_3(\mathbf{x})) \frac{\partial^2 q\_x(\mathbf{x}, \mathbf{y}, t)}{\partial t^2} - \overline{\epsilon} l\_3(\mathbf{x}) \frac{\partial^3 w\_0}{\partial \mathbf{x} \partial t^2},\tag{12a}$$

$$\frac{\partial \mathcal{N}\_{\mathbf{x}\mathbf{y}}}{\partial \mathbf{y}} + \frac{\partial \mathcal{N}\_{\mathbf{y}}}{\partial \mathbf{y}} = l\_{\mathbf{0}}(\mathbf{x}) \frac{\partial^{2} v\_{\mathbf{0}}(\mathbf{x}, \mathbf{y}, \mathbf{t})}{\partial \mathbf{t}^{2}} + \left( l\_{\mathbf{0}}(\mathbf{x}) - \mathsf{T}l\_{\mathbf{0}}(\mathbf{x}) \right) \frac{\partial^{2} \rho\_{\mathbf{y}}(\mathbf{x}, \mathbf{y}, \mathbf{t})}{\partial \mathbf{t}^{2}} - \mathsf{T}l\_{\mathbf{0}}(\mathbf{x}) \frac{\partial^{3} w\_{\mathbf{0}}(\mathbf{x}, \mathbf{y}, \mathbf{t})}{\partial \mathbf{y} \partial \mathbf{t}^{2}},\tag{12b}$$

∂*Qxz* <sup>∂</sup>*<sup>x</sup>* + ∂*Qxy* <sup>∂</sup>*<sup>y</sup>* − 3*c* ∂*Rxz* <sup>∂</sup>*<sup>x</sup>* + ∂*Rxy* ∂*y* + *c* ∂2*Pxx* <sup>∂</sup>*x*<sup>2</sup> + <sup>2</sup> ∂2*Pxy* <sup>∂</sup>*x*∂*<sup>y</sup>* + ∂2*Pyy* ∂*y*<sup>2</sup> <sup>−</sup> *Nyy R* −*kww*0(*x*, *y*, *t*) + *kg* ∂2*w*0(*x*,*y*,*t*) <sup>∂</sup>*x*<sup>2</sup> <sup>+</sup> <sup>∂</sup>2*w*0(*x*,*y*,*t*) ∂*y*<sup>2</sup> + *Fa* ∂2*w*0(*x*,*y*,*t*) <sup>∂</sup>*x*<sup>2</sup> = *<sup>I</sup>*0(*x*) ∂2*w*0(*x*,*y*,*t*) ∂*t*<sup>2</sup> −*c* 2 *I*6(*x*) ∂4*w*0(*x*,*y*,*t*) <sup>∂</sup>*x*2∂*t*<sup>2</sup> <sup>+</sup> <sup>∂</sup>4*w*0(*x*,*y*,*t*) ∂*y*2∂*t*<sup>2</sup> + *cI*3(*x*) ∂3*u*0(*x*,*y*,*t*) <sup>∂</sup>*x*∂*t*<sup>2</sup> <sup>+</sup> <sup>∂</sup>3*v*0(*x*,*y*,*t*) ∂*y*∂*t*<sup>2</sup> <sup>+</sup>*c*(*I*4(*x*) <sup>−</sup> *cI*6(*x*)) ∂3ϕ*x*(*x*,*y*,*t*) <sup>∂</sup>*x*∂*t*<sup>2</sup> + ∂3ϕ*y*(*x*,*y*,*t*) ∂*y*∂*t*<sup>2</sup> , (12c) ∂*Mxx* ∂*Mxy* ∂*Pxx* ∂*Pxy* 

$$\begin{aligned} \frac{\partial M\_{\text{rx}}}{\partial x} + \frac{\partial M\_{\text{xy}}}{\partial y} - Q\_{\text{xz}} + 3\tilde{c}R\_{\text{xz}} - \tilde{c} \Big( \frac{\partial P\_{\text{rx}}}{\partial x} + \frac{\partial P\_{\text{xy}}}{\partial y} \Big) &= (I\_{0}(\mathbf{x}) - \tilde{c}I\_{3}(\mathbf{x})) \frac{\partial^{2} u\_{0}(\mathbf{x},y,t)}{\partial t^{2}} \\ + (I\_{2}(\mathbf{x}) - 2\tilde{c}I\_{4}(\mathbf{x}) + \tilde{c}^{2}I\_{6}(\mathbf{x})) \frac{\partial^{2} \varphi\_{\text{x}}(\mathbf{x},y,t)}{\partial t^{2}} \\ - \tilde{c}(I\_{4}(\mathbf{x}) - \tilde{c}I\_{6}(\mathbf{x})) \frac{\partial^{3} \varphi\_{\text{xy}}(\mathbf{x},y,t)}{\partial x \partial t^{2}} \Big) \end{aligned} \tag{12d}$$

$$\begin{cases} \frac{\partial M\_{yy}}{\partial y} + \frac{\partial M\_{xy}}{\partial x} - Q\_{yz} + 3\overline{\varepsilon}R\_{yz} - \overline{\varepsilon} \Big( \frac{\partial P\_{yy}}{\partial y} + \frac{\partial P\_{xy}}{\partial x} \Big) = (I\_0(\mathbf{x}) - \overline{\varepsilon}I\_3(\mathbf{x})) \frac{\partial^2 v\_0(\mathbf{x}, y, t)}{\partial t^2} \\ + (I\_2(\mathbf{x}) - 2\overline{\varepsilon}I\_4(\mathbf{x}) + \overline{\varepsilon}^2 I\_6(\mathbf{x})) \frac{\partial^2 \varphi\_y(\mathbf{x}, y, t)}{\partial t^2} \\ - \overline{\varepsilon}(I\_4(\mathbf{x}) - \overline{\varepsilon}I\_6(\mathbf{x})) \frac{\partial^3 w\_0(\mathbf{x}, y, t)}{\partial y \partial t^2} \end{cases} \tag{12e}$$

Note that the resultants in Equation (12) are computed by integration of the pertaining stress components along the shell structure, i.e.,

$$\begin{Bmatrix} N\_{xx} \\ N\_{yy} \\ N\_{xy} \end{Bmatrix} = \int\_{-\frac{1}{2}}^{\frac{1}{2}} \begin{Bmatrix} \sigma\_{xx} \\ \sigma\_{yy} \\ \tau\_{xy} \end{Bmatrix} dz,\tag{13}$$

$$
\begin{Bmatrix} M\_{xx} \\ M\_{yy} \\ M\_{xy} \end{Bmatrix} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \begin{Bmatrix} \sigma\_{xx} \\ \sigma\_{yy} \\ \tau\_{xy} \end{Bmatrix} zdz,\tag{14}
$$

$$
\begin{Bmatrix} P\_{xx} \\ P\_{yy} \\ P\_{xy} \end{Bmatrix} = \int\_{-\frac{k}{2}}^{\frac{k}{2}} \begin{Bmatrix} \sigma\_{xx} \\ \sigma\_{yy} \\ \tau\_{xy} \end{Bmatrix} z^3 dz,\tag{15}
$$

$$\mathbb{E}\left\{\begin{array}{c} Q\_{xz} \\ Q\_{yz} \end{array}\right\} = \int\_{-\frac{\mathbb{E}}{2}}^{\frac{\mathbb{E}}{2}} \begin{Bmatrix} \tau\_{xz} \\ \tau\_{yz} \end{Bmatrix} dz,\tag{16}$$

$$\left\{ \begin{array}{c} R\_{xz} \\ R\_{yz} \end{array} \right\} = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \left\{ \begin{array}{c} \tau\_{xz} \\ \tau\_{yz} \end{array} \right\} z^2 dz. \tag{17}$$

The generalized inertia moments are defined as

$$\rho\left(I\_0(\mathbf{x}), I\_1(\mathbf{x}), I\_2(\mathbf{x}), I\_4(\mathbf{x}), I\_5(\mathbf{x}), I\_6(\mathbf{x})\right) = \bigcap\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} \{1, z, z^2, z^3, z^4, z^6\} \rho(\mathbf{x}, z) dz. \tag{18}$$

Three types of boundary conditions are considered along the shell edges, namely

➢ **Simply-Simply** (S-S) supports

$$\mathbf{x} = \mathbf{0}, \mathbf{L} \implies \mathbf{v}\_0 = w\_0 = q\_\mathbf{y} = \mathbf{M}\_\mathbf{x} = \mathbf{N}\_\mathbf{x} = \mathbf{0},\tag{19}$$

➢ **Clamped-Clamped** (C-C) supports

$$\mathfrak{x} = 0, L \implies \mathfrak{u}\_0 = \mathfrak{v}\_0 = \mathfrak{w}\_0 = \mathfrak{q}\_x = \mathfrak{q}\_y = 0,\tag{20}$$

### ➢ **Clamped-simply** (C-S) supports

$$
\mathbf{x} = \mathbf{0} \Rightarrow \mathbf{u}\_0 = \mathbf{v}\_0 = \mathbf{w}\_0 = q\_\mathbf{x} = q\_\mathbf{y} = \mathbf{0},\tag{21a}
$$

$$\mathbf{x} = L \Rightarrow \mathbf{v}\_0 = \mathbf{z} \\ \mathbf{v}\_0 = \boldsymbol{\varphi}\_y = \mathbf{M}\_x = \mathbf{N}\_x = \mathbf{0},\tag{21b}$$

#### **3. Solution Procedure**

In this section, we want to determine the dynamic stability of BD-FG porous cylindrical shells, where the governing equations of motion are expressed through the following expansion for the kinematic quantities

$$u\_0(\mathbf{x}, y, t) = \mathcal{U}(\mathbf{x}) \sin(n \frac{y}{R}) \overline{\mathcal{U}}(t) \tag{22a}$$

$$w\_0(\mathbf{x}, y, t) = V(\mathbf{x}) \cos(n \frac{y}{R}) \overline{V}(t) \tag{22b}$$

$$m\_0(\mathbf{x}, \mathbf{y}, t) = \mathcal{W}(\mathbf{x}) \sin(n \frac{\mathbf{y}}{R}) \overline{\mathcal{W}}(t) \tag{22c}$$

$$
\varphi\_x(\mathbf{x}, \mathbf{y}, t) = \Phi\_x(\mathbf{x}) \sin(n \frac{\mathcal{Y}}{R}) \overline{\Phi}\_x(t) \tag{22d}
$$

$$\Phi\_y(\mathbf{x}, y, t) = \Phi\_y(\mathbf{x}) \cos(n \frac{y}{R}) \overline{\Phi}\_y(t) \tag{22e}$$

where *n* is the circumferential half wave number, *U*(*t*), *V*(*t*), *W*(*t*), Φ*x*(*t*) and Φ*y*(*t*) are the time functions. The admissible displacement functions in Equation (22) satisfy both the equations of motion and their boundary conditions. Afterward, the governing equations of the problem are discretized according to the GDQ method.

Upon substitution of Equation (22a–e) and Equation (5a–e) into Equation (12), after a proper manipulation, we obtain the equations of motion in their final form, as detailed in Appendix A.

The above-mentioned equations of motion are solved numerically in a strong form by means of the GDQ method, as largely discussed in [49] and in a review paper [50] in terms of accuracy, stability and reliability, and successfully applied for many numerical applications, namely, the buckling, free vibration, or dynamic problems of composite structures [51–55], as well as the fracture mechanics problems [56,57] or non-linear transient problems [58,59]. In addition, the Bolotin method [60] is proposed herein to determine the DIRs for the differential equations system, known as Mathieu–Hill system of equations. More details about the basics of the proposed numerical tools are recalled in what follows.

### *3.1. The GDQ Method*

The GDQ method approximates the fundamental system of differential equations, by discretizing the derivatives of a function *J*(*x*) respect to a spatial variable at a given discrete grid distribution, by means of the weighting coefficients. For a one-dimensional problem where the whole domain 0 ≤ *x* ≤ *L* is discretized in *N* grids points, the approximation of the *n*th-order derivatives of *J* function respect to *x* variable can be expressed as [49]

$$\frac{d^n J(\mathbf{x}\_p)}{dx^n} = \sum\_{r=1}^N \chi\_{pr}^{(n)} J(\mathbf{x}\_r) \quad n = 1, 2, \dots, N - 1,\tag{29}$$

χ(*n*) *pr* being the weighting coefficients, defined as follows [49]

$$\chi\_{pq}^{(1)} = \frac{\mathbf{Y}(\mathbf{x}\_p)}{\left(\mathbf{x}\_p - \mathbf{x}\_q\right)\mathbf{Y}(\mathbf{x}\_q)}, \ p \neq q; \ p, q = 1, 2, \dots, N \tag{30}$$

and Υ *xp* is the Lagrangian operator expressed as [49]

$$\Upsilon(\mathbf{x}\_p) = \prod\_{\substack{p \neq q. \ q = 1}}^N (\mathbf{x}\_p - \mathbf{x}\_q). \tag{31}$$

For higher order derivatives of the weighting coefficients it is [49]

$$\chi\_{pq}^{(n)} = n \left( \chi\_{pp}^{(n-1)} \chi\_{pq}^{(1)} - \frac{\chi\_{pq}^{(n-1)}}{\left(\chi\_p - \chi\_q\right)} \right) \tag{32}$$

It is well known in literature that the type of grid distribution within the domain can affect significantly the accuracy of the proposed method [50]. In what follows we apply a Chebyshev–Gauss–Lobatto non-uniform pattern, due to its great performances, as verified by Shu [49]

$$\mathbf{x}\_p = \frac{\mathbf{L}}{2} \mathbf{\hat{z}} \left( \mathbf{1} - \cos \left( \frac{\mathbf{x}\_p - 1}{N - 1} \right) \pi \right) \ y = \mathbf{1}, 2, \dots, N. \tag{33}$$

Thus, the governing differential equations of motion and boundary conditions are discretized according to the GDQ approach, as detailed in Appendix B. Let us denote the periodic axial compressive load as

$$F\_a(t) = aF\_{cr} + \beta F\_{cr} \cos(\omega t) \tag{34}$$

where α and β refer to the static and dynamic load factors, respectively. Furthermore, *Fcr* denotes the critical static load and ω stands for the excitation frequency. By substituting Equation (34) into the third Equation (A3) from the Appendix A, and by combining the discretized equations of motion along with the associated boundary conditions, the problem can be redefined in the following matrix form

$$
\begin{pmatrix} M\_{bb} & M\_{bd} \\ M\_{db} & M\_{dd} \end{pmatrix} \begin{pmatrix} \ddot{\Gamma}\_{b} \\ \ddot{\Gamma}\_{d} \end{pmatrix} + \begin{pmatrix} \begin{pmatrix} K\_{bb} & K\_{bd} \\ K\_{db} & K\_{dd} \end{pmatrix} + F\_{\mathcal{F}l}(\alpha + \beta \cos(\omega t)) \begin{pmatrix} K\_{bb}^{G} & K\_{bd}^{G} \\ K\_{db}^{G} & K\_{dd}^{G} \end{pmatrix} \end{pmatrix} \begin{pmatrix} \Gamma\_{b} \\ \Gamma\_{d} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} \tag{35}
$$

where *M*, *K*, and *K<sup>G</sup>* are the mass, stiffness, and geometric stiffness matrixes, respectively, and Γ = & *U*, *V*, *W*, Φ*x*, Φ*<sup>y</sup>* '*<sup>T</sup>* denotes the unknown dynamic displacement vector. In addition, indexes *<sup>b</sup>* and *d* indicate the boundary points and domain points, respectively.

#### *3.2. Bolotin Method*

The second order system of differential Equation (35) is known in literature as Mathieu–Hill system of equations due to presence of the periodic coefficient, accordingly. In the present study we propose the Bolotin method [60] to define the boundaries associated to the DIR of the BD-FG porous cylindrical shell. Based on this method, the dynamic displacement vector Γ can be defined in a Fourier series as follows [60]

$$\langle \Gamma \rangle = \sum\_{s=1,3,\dots}^{\infty} \left( \langle \vartheta\_s \rangle \sin \left( \frac{s\omega t}{2} \right) + \langle \upsilon\_s \rangle \cos \left( \frac{s\omega t}{2} \right) \right) \tag{36}$$

where ϑ*<sup>s</sup>* and υ*<sup>s</sup>* denote the arbitrary time invariant vectors. It should be mentioned that the first DIR with period 2*T* is generally much meaningful and wider than the secondary one with period *T*. For this reason, in this work we consider the solutions with period 2*T*. By substitution of Equation (36) into Equation (35) and by mathematical manipulation, we get the following first order equation

$$\left| \left[ K \right] - F\_{cr} \left( \alpha \pm \frac{\beta}{2} \right) \left[ K^G \right] - \frac{\alpha^2}{4} \left[ M \right] \right| = 0,\tag{37}$$

which represents a classical eigenvalue problem. The critical buckling load can be computed from Equation (35) by neglecting the inertia terms and by setting α + β cos(ω*t*) = 1. Then, solving Equation (37) for some fixed values of α and *Fcr*, the variation of the excitation frequency ω in regards to β can be drawn as DIR for the BD-FG structure.

#### **4. Numerical Investigation**

In this section some illustrative example are shown, starting with a preliminary validation of the proposed method with respect to the available literature, and continuing with a parametric investigation of the problem, whose results are evaluated comparatively in order to evaluate the sensitivity of the mechanical response.

#### *4.1. Validation*

Due to the general lack of works in the literature on the dynamic buckling behavior of BD-FG porous cylindrical shells, the proposed model is validated, herein, for an axial buckling problem of a simply supported conventional FG cylindrical shell based on two different theories. Thus, for comparative purposes, we select the same material properties and shell geometry as reported in [48,61], and neglect the inertia terms, foundation parameters and porosity effects, while assuming a null value for *nx*. In Table 1 we summarize the results in terms of critical axial buckling load *Fcr* for a FG cylindrical shell with *h* = 0.001 m, *L*/*R* = 0.5, *Ec* = 380 GPa, *Em* = 70 GPa, ν = 0.3, *c* = 0, with a clear excellent agreement between our results and predictions in Ref. [61] based on the first order shear deformation theory. This first numerical example could be considered as limit case, where the TSDT reverts to the FSDT, since it refers to a thin shell structure, just for validation purposes. More

accurate results, however, are always expected under a TSDT assumption for increased values of the shell thickness, as done in the next parametric investigation.


**Table 1.** Comparative evaluation of the critical axial buckling load (MN) for a FG cylindrical shell with *h* = 0.001 m, *L*/*R* = 0.5, *Ec* = 380 GPa, *Em* = 70 GPa, ν = 0.3, *nx* = 0.

As further comparative study, Table 2 compares the dimensionless critical buckling load (*Pcr* = *FcrL*2/π2*Dm*; *Dm* = *Emh*3/12 <sup>1</sup> <sup>−</sup> <sup>ν</sup><sup>2</sup> *m* ) for a FG cylindrical shell with *h* = 0.001 m, *Ec* = 380 GPa, *Em* = 70 GPa, ν = 0.3, *c* = 4/3*h*2. Based on Table 2, it is worth noticing the high precision between our results and predictions by Bagherizadeh et al. [48], which confirms the accuracy of the GDQ method. This method is thus proposed in the following parametric study, as efficient numerical tool to solve the problem.


**Table 2.** Comparative evaluation of the dimensionless critical axial buckling load for a FG cylindrical shell with *h* = 0.001 m, *Em* = 70 GPa, ν = 0.3, *nz* = 2.

### *4.2. Parametric Study*

We refer to a BD-FG cylindrical structure with constituent phases of properties listed in Table 3, where the following dimensionless parameters are considered to compute the dimensionless structural excitation frequencies

$$
\Omega = \omega R \sqrt{\frac{\rho\_m}{E\_m}},\\ K\_\% = \frac{k\_\% R^2}{E\_m h^3},\\ K\_\text{\textmu} = \frac{k\_\text{\textmu} R^4}{E\_m h^3}.\tag{38}
$$

**Table 3.** Material properties of the BD-FG cylindrical shell.


We determine the DIR, and highlight the effects of different parameters such as the thickness-to-radius ratio (*h*/*R*), the length-to-radius ratio (*L*/*R*), the static load factor, the boundary conditions, the power law indexes (*nx*,*nz*), the type and volume fraction of porosity, and the foundation parameters, on the dynamic buckling behavior of the BD-FG cylindrical shell.

In Figure 3 we plot the variation of the DIR for different thickness-to-radius ratios (*h*/*R*), where a clear shift of the DIR is observed for increasing *h*/*R* ratios. This means that the DIR becomes wider for a certain value of the dynamic load factor, and occurs with a sort of delay. An increased *h*/*R* ratio from 0.01 to 0.1 yields a global shift of the DIR origin point towards high excitation frequencies.

Figure 4 shows the sensitivity of the DIR for varying *L*/*R* ratios, while keeping fixed *h*/*R* = 0.01. In detail, for *L*/*R* = 1 the structure has a wider DIR in comparison with the other values, whereby an increasing *L*/*R* ratio yields the DIR to take place at lower excitation frequencies. Based on the plots in Figure 4, we can observe a reduction of about 41.96% in the excitation frequencies corresponding to the origin of the instability region, for a *L*/*R* ranging between 1 and 10. When the *L*/*R* ratio features higher magnitudes, the bending resistance gradually reduces and yields an increased bending deformation.

**Figure 3.** Effect of the thickness-to-radius ratio on the dynamic instability region (DIR) for a BD-FG cylindrical shell with *R* = 0.5 m, *L*/*R* = 1, *nx* = *nz* = 1, α = 0.3.

**Figure 4.** Effect of the length-to-radius ratio on the DIR for a BD-FG cylindrical shell with *R* = 0.5 m, *h*/*R* = 0.02, *nx* = *nz* = 1, α = 0.3.

In Figure 5, we evaluate the effect of the static load factor on the instability region of the BD-FG cylindrical shell. As expected, in absence of a static load on the structure, the width of DIR gets smaller, whereby for an increased static load factor, it becomes gradually greater for a fixed value of dynamic load factor (i.e., β = 1), and its origin tends to move on the left side. This proves the sensitivity of the structural instability to the static load factor.

**Figure 5.** Effect of the static load factor on the DIR for a BD-FG cylindrical shell with *R* = 0.5 m, *h*/*R* = 0.02, *L*/*R* = 1, *nx* = *nz* = 1, α = 0.3.

As also visible in Figure 6, we evaluate the impact of different boundary conditions on the DIR of the cylindrical shell. Here, we consider three different boundary conditions, namely, S-S, C-S, and C-C boundary conditions. Based on the plots of Figure 7, it is worth noting that a C-C boundary conditions yields higher values of the dimensionless excitation frequencies than those ones provided by a S-S or C-S supports, due to an increased stiffness of the structure. Furthermore, the origin of the instability region tends to get away from the origin. Once the dynamic load factor β reaches the unit value, the width of the DIR for S-S boundary condition becomes smaller, compared to the other boundary conditions. This means that, for lower values of dimensionless excitation frequency, a BD-FG cylindrical shell with S-S supports tends to become more unstable compared to the other boundary conditions.

**Figure 6.** Effect of boundary conditions on the DIR for a BD-FG cylindrical shell with *R* = 0.5 m, *h*/*R* = 0.02, *L*/*R* = 1, *nx* = *nz* = 1, α = 0.3.

**Figure 7.** Effect of the longitudinal power index on the DIR for a FG cylindrical shell with *R* = 0.5 m, *h*/*R* = 0.02, *L*/*R* = 1, *nz* = 0, α = 0.3.

A further investigation is devoted to study the influence of the power law index along the length *nx* on the dynamic buckling behavior of one-directional FG cylindrical shell, as depicted Figure 7. In such a case, the DIR takes place at lower frequencies owning to an increased magnitude of the power law index. The effect of an increased dimensionless excitation frequency related to the origin of the DIR is meaningful within the range 0.2 ≤ *nx* ≤ 5. For greater values of *nx*, the variation in frequency corresponding to the origin of DIR becomes less remarkable. A one-directional FG cylindrical shell with *nx* = 10 or *nx* = 8 is more sensitive to the dynamic instability for lower excitation frequencies compared with those ones with *nx* ≤ 5.

For a conventional FG cylindrical shell, we also investigate the effect of the transverse power index *nz* on the DIR, as plotted in Figure 8. It can be mentioned that for a constant value of the dynamic load factor, the enhancement of *nz* yields a reduction width of the dynamic instability region, especially for lower excitation frequencies. In addition, the origin DIR moves to a lower dimensionless excitation frequency. Comparing the results from Figures 7 and 8, it can be concluded that, a double increase of both *nx*, and *nz*, leads to a reduction in the excitation frequency. Nevertheless, *nx* plays an important role in the reduction of the excitation frequency and in the increase of the structural instability. For β = 0, for example, we note a reduction of the excitation frequency equal to 49.3% and 44.49%, for and increasing value of *nx* and *nz* from 0.2 up to 10, respectively.

Figure 9 shows the effect of the transverse and longitudinal volume fraction indexes on the DIR. In detail, for an increase of these two parameters, the origin of DIR moves to higher dimensionless excitation frequency, and the DIR width declines. In conclusion, the double increase of *nx* and *nz* leads a metal phase reinforcement, with a overall decrease of the structural stiffness and an increase in the structural instability. BD-FG cylindrical shells with lower values of *nz* and *nx*, are less sensitive to the dynamic instability, due to their higher stiffness. The contrary occurs for higher values of *nx* and *nz* due to an increased deformability of the structure. The importance of applying BD-FG materials is highlighted when the variation of material properties is considered in a single direction or more directions simultaneously. This issue can be beneficial for the fabrication and design purposes of modern FG structures. Subsequently, the dynamic stability of BD-FG cylindrical shells can be controlled selecting appropriate power indexes corresponding to the desired direction.

**Figure 8.** Effect of the transverse power index on the DIR for a BD-FG cylindrical shell with *R* = 0.5 m, *h*/*R* = 0.02, *L*/*R* = 1, *nx* = 0, α = 0.3.

**Figure 9.** Effect of the longitudinal and transverse power indexes on the DIR for a FG cylindrical shell with *R* = 0.5 m, *h*/*R* = 0.02, *L*/*R* = 1, α = 0.3.

Moreover, Figure 10 shows the effect of an even porosity volume fraction on the dimensionless excitation frequency of a BD-FG porous cylindrical shell. Note that the dimensionless excitation frequency decreases for increasing values of *nx* and *nz*. In addition, for a fixed value of *nx*, *nz*, the excitation frequency tends to converge to a common point. As expected, at the intersection point, the effect of an even porosity volume fraction, ζ, on the dimensionless excitation frequency is almost negligible. However, for different values of *nx*, *nz*, the effect of even porosity between the ceramic and metal phases can change significantly, such that before the intersection point, an increased even porosity volume fraction increases the dimensionless excitation frequency, and the contrary occurs

after the intersection point, with a gradual decrease in the excitation frequency for an enhanced ζ. The additional Figures 11 and 12 show the different response for different values of ζ, while assuming the same value for *nx* and *nz*. In detail, under the assumption *nx* = *nz* = 0.15 (Figure 11), it is visible that an increased porosity ζ moves the DIR towards higher excitation frequencies. A reversed behavior occurs in Figure 12 under the assumption *nx* = *nz* = 1.5, since an increased value of ζ causes a general shift of the DIR to lower excitation frequencies. This confirms the effect of either the even porosity volume fraction ζ and the power indexes *nx* and *nz* on the stability response of the structure.

**Figure 10.** Effect of the even porosity volume fraction ζ on the dimensionless excitation frequency versus the longitudinal and transverse power indexes for a BD-FG porous cylindrical shell with *R* = 0.5 m, *h*/*R* = 0.02, *L*/*R* = 1, α = 0.3, β = 0.

**Figure 11.** Effect of the even porosity volume fraction ζ on the DIR for a BD-FG porous cylindrical shell with *R* = 0.5 m, *h*/*R* = 0.02, *L*/*R* = 1, α = 0.3, *nx* = *nz* = 0.15.

**Figure 12.** Effect of the even porosity volume fraction ζ on the DIR for a BD-FG porous cylindrical shell with *R* = 0.5 m, *h*/*R* = 0.02, *L*/*R* = 1, α = 0.3, *nx* = *nz* = 0.15.

In Figures 13–15 we repeat the parametric analysis to evaluate the effect of an uneven porosity between two phases of the second ceramic and second metal. In detail, Figure 13 shows the variation of the dimensionless excitation frequency versus *nx*, *nz*, for different values of ξ. Figure 14 is devoted to check for the influence of ξ on the DIR of the structure for *nx* = *nz* = 0.15. Additionally, in this case, we can observe as the DIR moves to the right side by increasing ξ and it takes place at higher excitation frequencies. Nevertheless, by assuming *nx* = *nz* = 1.5, a different trend is noticed for the DIR in Figure 15, since an increased value of ξ yields the DIR to occur at lower excitation frequencies and its width gets smaller.

**Figure 13.** Effect of the uneven porosity volume fraction ξ on the dimensionless excitation frequency versus the longitudinal and transverse power indexes for a BD-FG porous cylindrical shell with *R* = 0.5 m, *h*/*R* = 0.02, *L*/*R* = 1, α = 0.3, β = 0.

**Figure 14.** Effect of the uneven porosity volume fraction ξ on the DIR for a BD-FG porous cylindrical shell with *R* = 0.5 m, *h*/*R* = 0.02, *L*/*R* = 1, α = 0.3, *nx* = *nz* = 0.15.

**Figure 15.** Effect of the uneven porosity volume fraction ξ on the DIR for a BD-FG porous cylindrical shell with *R* = 0.5 m, *h*/*R* = 0.02, *L*/*R* = 1, α = 0.3, *nx* = *nz* = 0.15.

The last parametric analysis considers the possible sensitivity of the response to the elastic foundation. For this reason, Figures 16 and 17 plot the variation of the DIR with the Winkler or the Pasternak foundation coefficients, respectively. A noteworthy increase in stiffness emerges from both figures, where the origin of the DIR moves towards higher values of frequency. According to a comparative evaluation of the results, it seems that the best dynamic behavior of the cylindrical shell is reached for a structure surrounded by a Pasternak elastic foundation. Hence, the effect of the Pasternak elastic coefficient is more remarkable than the Winkler-based one, where a BD-FG cylindrical shell becomes more stable if embedded in a Pasternak foundation.

**Figure 16.** Effect of the Winkler coefficient on the DIR for a BD-FG cylindrical shell with *R* = 0.5 m, *h*/*R* = 0.02, *L*/*R* = 1, α = 0.3, *Kg* = 0, *nx* = *nz* = 1.5.

**Figure 17.** Effect of the Pasternak coefficient on the DIR for a BD-FG cylindrical shell with *R* = 0.5 m, *h*/*R* = 0.02, *L*/*R* = 1, α = 0.3, *Kw* = 0, *nx* = *nz* = 1.

#### **5. Conclusions**

This work investigates the dynamic stability of BD-FG cylindrical shells embedded in an elastic foundation, including possible effects related to porosity. The material properties of BD-FG porous cylindrical shells are computed according to a modified BD power law model. Using the Hamilton's principle, we determine the governing equations of the problem, under the classical TSDT assumptions. The aforementioned equations are rewritten into a system of Mathieu–Hill equations, according to a GDQ approach. The work is also devoted to determine the DIR of the structure while applying the Bolotin method. After a preliminary validation of the proposed formulation, with respect to the

available literature, we perform a large numerical investigation to check for the sensitivity of the response both in terms of excitation frequencies and DIRs, for different thickness-to-radius ratios, length-to-radius ratios, boundary conditions, transverse and longitudinal power law indexes, even and uneven porosity volume fractions, and foundation parameters. Based on the systematic numerical investigation, the main conclusions can be summarized as follows


**Author Contributions:** Conceptualization, F.A., H.T., R.D. and F.T.; Formal analysis, F.A., H.T., R.D. and F.T.; Investigation, F.A., H.T. and F.T.; Validation, H.T., R.D. and F.T.; Writing—Original Draft, F.A., H.T., R.D. and F.T.; Writing—Review & Editing, R.D. and F.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

Here below are the equations of motion in their final form

 d*A*11(*x*) d*x* d*U*(*x*) d*x* + *A*11(*x*) d2*U*(*x*) d*x*<sup>2</sup> <sup>−</sup> *<sup>n</sup>*2*A*66(*x*)*U*(*x*) *R*2 −*nV*(*x*) *R* d*A*12(*x*) d*x* <sup>−</sup> *nA*12(*x*) *R* d*V*(*x*) d*x* <sup>−</sup> *nA*66(*x*) *R* d*V*(*x*) d*x* + *<sup>W</sup>*(*x*) *R* d*A*12(*x*) d*x* + *<sup>A</sup>*12(*x*) *R* d*W*(*x*) d*x* − *c* d*D*11(*x*) d*x* d2*W*(*x*) d*x*<sup>2</sup> −*cD*11(*x*) d3*W*(*x*) d*x*<sup>3</sup> <sup>+</sup> *cn*2*W*(*x*) *R*2 d*D*12(*x*) d*x* <sup>+</sup> <sup>2</sup>*cn*2*D*66(*x*) *R*2 d*W*(*x*) d*x* <sup>+</sup>*cn*2*D*12(*x*) *R*2 d*W*(*x*) d*x* + d*B*11(*x*) d*x* dΦ*x*(*x*) d*x* + *B*11(*x*) d2Φ*x*(*x*) d*x*<sup>2</sup> −*cD*11(*x*) d2Φ*x*(*x*) d*x*<sup>2</sup> − *c* d*D*11(*x*) d*x* dΦ*x*(*x*) d*x* <sup>−</sup> *<sup>n</sup>*2*B*66(*x*)Φ*x*(*x*) *<sup>R</sup>*<sup>2</sup> <sup>+</sup> *cn*2*D*66(*x*)Φ*x*(*x*) *R*2 −*n*Φ*y*(*x*) *R* d*B*12(*x*) d*x* <sup>−</sup> *nB*12(*x*) *R* dΦ*y*(*x*) d*x* <sup>−</sup> *nB*66(*x*) *R* dΦ*y*(*x*) d*x* +*cnD*12(*x*) *R* dΦ*y*(*x*) d*x* <sup>+</sup> *cn*Φ*y*(*x*) *a* d*D*12(*x*) d*x* + *cnD*66(*x*) *a* dΦ*y*(*x*) d*x* = <sup>1</sup> *U*(*t*) *I*0(*x*)*U*(*x*) *d*2*U*(*t*) *dt*<sup>2</sup> + (*I*0(*x*) <sup>−</sup> *cI*3(*x*))Φ*x*(*x*) *d*2Φ*x*(*t*) *dt*<sup>2</sup> −*cI*3(*x*) d*W*(*x*) d*x d*2*W*(*t*) *dt*<sup>2</sup> (A1)

*nU*(*x*) *R* d*A*66(*x*) d*x* + *nA*66(*x*) *R* d*U*(*x*) d*x* + *nA*12(*x*) *R* d*U*(*x*) d*x* + d*A*66(*x*) d*x* d*V*(*x*) d*x* + *A*66(*x*) d2*V*(*x*) d*x*<sup>2</sup> <sup>−</sup> *<sup>n</sup>*2*A*11(*x*)*V*(*x*) *R*2 <sup>−</sup>2*cn R* d*D*66(*x*) d*x* d*W*(*x*) d*x* <sup>−</sup> <sup>2</sup>*cnD*66(*x*) *R* d2*W*(*x*) d*x*<sup>2</sup> + *nA*11(*x*)*W*(*x*) *R*2 <sup>+</sup>*cn*3*D*11(*x*)*W*(*x*) *<sup>R</sup>*<sup>3</sup> <sup>−</sup> *cnD*12(*x*) *R* d2*W*(*x*) d*x*<sup>2</sup> +*n*Φ*x*(*x*) *R* d*B*66(*x*) d*x* + *nB*66(*x*) *R* dΦ*x*(*x*) d*x* <sup>−</sup> *cn*Φ*x*(*x*) *R* d*D*66(*x*) d*x* −*cnD*66(*x*) *R* dΦ*x*(*x*) d*x* + *nB*12(*x*) *R* dΦ*x*(*x*) d*x* <sup>−</sup> *cnD*12(*x*) *R* dΦ*x*(*x*) d*x* + d*B*66(*x*) d*x* dΦ*y*(*x*) d*x* + *B*66(*x*) d2Φ*y*(*x*) d*x*<sup>2</sup> − *c* d*D*66(*x*) d*x* dΦ*y*(*x*) d*x* −*cD*66(*x*) d2Φ*y*(*x*) d*x*<sup>2</sup> <sup>−</sup> *<sup>n</sup>*2*B*11(*x*)Φ*y*(*x*) *<sup>R</sup>*<sup>2</sup> <sup>+</sup> *cn*2*D*11(*x*)Φ*y*(*x*) *R*2 = <sup>1</sup> *V*(*t*) *I*0(*x*)*V*(*x*) *d*2*V*(*t*) *dt*<sup>2</sup> + (*I*0(*x*) <sup>−</sup> *cI*3(*x*))Φ*y*(*x*) *d*2Φ*y*(*t*) *dt*<sup>2</sup> <sup>−</sup>*cn <sup>R</sup> I*3(*x*)*W*(*x*) *d*2*W*(*t*) *dt*<sup>2</sup> (A2) −2*cn*2*U*(*x*) *R*2 d*D*66(*x*) d*x* <sup>−</sup> <sup>2</sup>*cn*2*D*66(*x*) *R*2 d*U*(*x*) d*x* −*cn*2*D*12(*x*) *R*2 d*U*(*x*) d*x* <sup>−</sup> *<sup>A</sup>*12(*x*) *R* d*U*(*x*) d*x* +*c* d2*D*11(*x*) d*x*<sup>2</sup> d*U*(*x*) d*x* + 2*c* d*D*11(*x*) d*x* d2*U*(*x*) d*x*<sup>2</sup> + *cD*11(*x*) d3*U*(*x*) d*x*<sup>3</sup> −*cnV*(*x*) *R* d2*D*12(*x*) d*x*<sup>2</sup> <sup>−</sup> <sup>2</sup>*cn R* d*D*12(*x*) d*x* d*V*(*x*) d*x* <sup>−</sup> *cnD*12(*x*) *R* d2*V*(*x*) d*x*<sup>2</sup> <sup>−</sup>2*cn R* d*D*66(*x*) d*x* d*V*(*x*) d*x* <sup>−</sup> <sup>2</sup>*cnD*66(*x*) *R* d2*V*(*x*) d*x*<sup>2</sup> <sup>+</sup> *cn*3*D*11(*x*)*V*(*x*) *<sup>R</sup>*<sup>3</sup> <sup>+</sup> *nA*11(*x*)*V*(*x*) *R*2 +2*<sup>c</sup>* 2 *n*2 *R*2 d*G*12(*x*) d*x* d*W*(*x*) d*x* + <sup>4</sup>*<sup>c</sup>* 2 *n*2 *R*2 d*G*66(*x*) d*x* d*W*(*x*) d*x* − *c* 2 *n*4*G*11(*x*)*W*(*x*) *R*4 +*c* 2 *n*2*W*(*x*) *R*2 d2*G*12(*x*) d*x*<sup>2</sup> + <sup>2</sup>*<sup>c</sup>* 2 *n*2*G*12(*x*) *R*2 d2*W*(*x*) d*x*<sup>2</sup> <sup>−</sup> <sup>9</sup>*<sup>c</sup>* 2 *n*2*F*66(*x*)*W*(*x*) *R*2 +4*<sup>c</sup>* 2 *n*2*G*66(*x*) *R*2 d2*W*(*x*) d*x*<sup>2</sup> <sup>+</sup> <sup>6</sup>*cn*2*C*66(*x*)*W*(*x*) *<sup>R</sup>*<sup>2</sup> <sup>−</sup> <sup>2</sup>*cn*2*D*11(*x*)*W*(*x*) *R*3 −*n*2*A*66(*x*)*W*(*x*) *<sup>R</sup>*<sup>2</sup> <sup>+</sup> *cW*(*x*) *R* d2*D*12(*x*) d*x*<sup>2</sup> + <sup>2</sup>*<sup>c</sup> R* d*D*12(*x*) d*x* d*W*(*x*) d*x* +2*cD*12(*x*) *R* d2*W*(*x*) d*x*<sup>2</sup> + *A*66(*x*) d2*W*(*x*) d*x*<sup>2</sup> + d*A*66(*x*) d*x* d*W*(*x*) d*x* −2*c* 2 d*G*11(*x*) d*x* d3*W*(*x*) d*x*<sup>3</sup> <sup>−</sup> *<sup>c</sup>*2*G*11(*x*) d4*W*(*x*) d*x*<sup>4</sup> <sup>−</sup> *<sup>A</sup>*11(*x*)*W*(*x*) *R*2 −6*c* d*C*66(*x*) d*x* d*W*(*x*) d*x* − 6*cC*66(*x*) d2*W*(*x*) d*x*<sup>2</sup> +2*<sup>c</sup>* 2 *n*2*G*66(*x*) *R*2 dΦ*x*(*x*) d*x* + <sup>2</sup>*<sup>c</sup>* 2 *n*2Φ*x*(*x*) *R*2 d*G*66(*x*) d*x* <sup>−</sup> *cn*2*F*12(*x*) *R*2 dΦ*x*(*x*) d*x* −2*cn*2Φ*x*(*x*) *R*2 d*F*66(*x*) d*x* <sup>−</sup> <sup>2</sup>*cn*2*F*66(*x*) *R*2 dΦ*x*(*x*) d*x* + *<sup>c</sup>* 2 *n*2*G*12(*x*) *R*2 dΦ*x*(*x*) d*x* +*cD*12(*x*) *R* dΦ*x*(*x*) d*x* + d*A*66(*x*) d*x* Φ*x*(*x*) + *A*66(*x*) dΦ*x*(*x*) d*x* −*c* 2 *G*11(*x*) d3Φ*x*(*x*) d*x*<sup>3</sup> <sup>−</sup> *<sup>B</sup>*12(*x*) *R* dΦ*x*(*x*) d*x* − 6*c* d*C*66(*x*) d*x* Φ*x*(*x*) −6*cC*66(*x*) dΦ*x*(*x*) d*x* + 9*c* 2 d*F*66(*x*) d*x* Φ*x*(*x*) + 9*c* 2 *F*66(*x*) dΦ*x*(*x*) d*x* +*cF*11(*x*) d3Φ*x*(*x*) d*x*<sup>3</sup> + *c* d2*F*11(*x*) d*x*<sup>2</sup> dΦ*x*(*x*) d*x* + 2*c* d*F*11(*x*) d*x* d2Φ*x*(*x*) d*x*<sup>2</sup> −*c* 2 d2*G*11(*x*) d*x*<sup>2</sup> dΦ*x*(*x*) d*x* − 2*c* 2 d*G*11(*x*) d*x* d2Φ*x*(*x*) d*x*<sup>2</sup> +2*<sup>c</sup>* 2 *n R* d*G*12(*x*) d*x* dΦ*y*(*x*) d*x* − *c* 2 *n*3*G*11(*x*)Φ*y*(*x*) *<sup>R</sup>*<sup>3</sup> <sup>−</sup> <sup>9</sup>*<sup>c</sup>* 2 *nF*66(*x*)Φ*y*(*x*) *R* +2*<sup>c</sup>* 2 *n R* d*G*66(*x*) d*x* dΦ*y*(*x*) d*x* + *<sup>c</sup>* 2 *nG*12(*x*) *R* d2Φ*y*(*x*) d*x*<sup>2</sup> <sup>+</sup> *cn*3*F*11(*x*)Φ*y*(*x*) *R*3 (A3)

−*cnD*11(*x*)Φ*y*(*x*) *<sup>R</sup>*<sup>2</sup> <sup>+</sup> <sup>2</sup>*<sup>c</sup>* 2 *nG*66(*x*) *R* d2Φ*y*(*x*) d*x*<sup>2</sup> <sup>+</sup> <sup>6</sup>*cnC*66(*x*)Φ*y*(*x*) *R* <sup>−</sup>2*cn R* d*F*66(*x*) d*x* dΦ*y*(*x*) d*x* <sup>−</sup> <sup>2</sup>*cnF*66(*x*) *R* d2Φ*y*(*x*) d*x*<sup>2</sup> <sup>−</sup> *cn*Φ*y*(*x*) *R* d2*F*12(*x*) d*x*<sup>2</sup> <sup>−</sup>2*cn R* d*F*12(*x*) d*x* dΦ*y*(*x*) d*x* <sup>−</sup> *cnF*12(*x*) *R* d2Φ*y*(*x*) d*x*<sup>2</sup> +*c* 2 *n*Φ*y*(*x*) *R* d2*G*12(*x*) d*x*<sup>2</sup> <sup>+</sup> *nB*11(*x*)Φ*y*(*x*) *<sup>R</sup>*<sup>2</sup> <sup>−</sup> *nA*66(*x*)Φ*y*(*x*) *R* −*kwW*(*x*) + *kg d*2*W*(*x*) *dx*<sup>2</sup> <sup>−</sup> *<sup>n</sup>*<sup>2</sup> *<sup>R</sup>*<sup>2</sup> *<sup>W</sup>*(*x*) − *Fa d*2*W*(*x*) *dx*<sup>2</sup> <sup>=</sup> <sup>1</sup> *W*(*t*) *I*0(*x*)*W*(*x*) *d*2*W*(*t*) *dt*<sup>2</sup> −*c* 2 *I*6(*x*) *d*2*W*(*x*) *dx*<sup>2</sup> *d*2*W*(*t*) *dt*<sup>2</sup> <sup>−</sup> *<sup>n</sup>*<sup>2</sup> *<sup>R</sup>*<sup>2</sup> *<sup>W</sup>*(*x*) *d*2*W*(*t*) *dt*<sup>2</sup> + *cI*3(*x*) *dU*(*x*) *dx d*2*U*(*t*) *dt*<sup>2</sup> −*n <sup>R</sup>V*(*x*) *d*2*V*(*t*) *dt*<sup>2</sup> <sup>+</sup> *<sup>c</sup>*(*I*4(*x*) <sup>−</sup> *cI*6(*x*)) *d*Φ*x*(*x*) *dx d*2Φ*x*(*t*) *dt*<sup>2</sup> <sup>−</sup> *<sup>n</sup> <sup>R</sup>*Φ*y*(*x*) *d*2Φ*y*(*t*) *dt*<sup>2</sup> 

*cn*2*D*66(*x*)*U*(*x*) *<sup>R</sup>*<sup>2</sup> <sup>−</sup> *<sup>n</sup>*2*B*66(*x*)*U*(*x*) *<sup>R</sup>*<sup>2</sup> + d*B*11(*x*) d*x* d*U*(*x*) d*x* + *B*11(*x*) d2*U*(*x*) d*x*<sup>2</sup> −*c* d*D*11(*x*) d*x* d*U*(*x*) d*x* − *cD*11(*x*) d2*U*(*x*) d*x*<sup>2</sup> −*nB*66(*x*) *R* d*V*(*x*) d*x* <sup>−</sup> *nV*(*x*) *R* d*B*12(*x*) d*x* <sup>−</sup> *nB*12(*x*) *R* d*V*(*x*) d*x* +*cnD*12(*x*) *R* d*V*(*x*) d*x* + *cnD*66(*x*) *R* d*V*(*x*) d*x* + *cnV*(*x*) *R* d*D*12(*x*) d*x* −*cW*(*x*) *R* d*D*12(*x*) d*x* <sup>−</sup> *cD*12(*x*) *R* d*W*(*x*) d*x* <sup>−</sup> *<sup>c</sup>*2*n*2*G*12(*x*) *R*2 d*W*(*x*) d*x* <sup>+</sup>*cn*2*F*12(*x*) *R*2 d*W*(*x*) d*x* <sup>+</sup> *cn*2*W*(*x*) *R*2 d*F*12(*x*) d*x* −2*c* 2 *n*2*G*66(*x*) *R*2 d*W*(*x*) d*x* <sup>+</sup> <sup>2</sup>*cn*2*F*66(*x*) *R*2 d*W*(*x*) d*x* −*c*2*n*2*W*(*x*) *R*2 d*G*12(*x*) d*x* − *A*66(*x*) d*W*(*x*) d*x* − 9*c* 2 *F*66(*x*) d*W*(*x*) d*x* +*c* 2 d*G*11(*x*) d*x* d2*W*(*x*) d*x*<sup>2</sup> + *c* 2 *G*11(*x*) d3*W*(*x*) d*x*<sup>3</sup> + *<sup>W</sup>*(*x*) *R* d*B*12(*x*) d*x* + *<sup>B</sup>*12(*x*) *R* d*W*(*x*) d*x* − *c* d*F*11(*x*) d*x* d2*W*(*x*) d*x*<sup>2</sup> −*cF*11(*x*) d3*W*(*x*) d*x*<sup>3</sup> + 6*cC*66(*x*) d*W*(*x*) d*x* −*n*2*C*66(*x*)Φ*x*(*x*) *<sup>R</sup>*<sup>2</sup> <sup>+</sup> <sup>2</sup>*cn*2*F*66(*x*)Φ*x*(*x*) *<sup>R</sup>*<sup>2</sup> <sup>−</sup> *<sup>c</sup>* 2 *n*2*G*66(*x*)Φ*x*(*x*) *<sup>R</sup>*<sup>2</sup> + d*C*11(*x*) d*x* dΦ*x*(*x*) d*x* +*C*11(*x*) d2Φ*x*(*x*) d*x*<sup>2</sup> − *A*66(*x*)Φ*x*(*x*) − 9*c* 2 *F*66(*x*)Φ*x*(*x*) + *c* 2 d*G*11(*x*) d*x* dΦ*x*(*x*) d*x* +*c* 2 *G*11(*x*) d2Φ*x*(*x*) d*x*<sup>2</sup> − 2*c* d*F*11(*x*) d*x* dΦ*x*(*x*) d*x* −2*cF*11(*x*) d2Φ*x*(*x*) d*x*<sup>2</sup> <sup>+</sup> <sup>6</sup>*cC*66(*x*)Φ*x*(*x*) <sup>−</sup> *nC*66(*x*) *R* dΦ*y*(*x*) d*x* −*n*Φ*y*(*x*) *R* d*C*12(*x*) d*x* <sup>−</sup> *nC*12(*x*) *R* dΦ*y*(*x*) d*x* − *c* 2 *nG*66(*x*) *a* dΦ*y*(*x*) d*x* +2*cnF*66(*x*) *R* dΦ*y*(*x*) d*x* + <sup>2</sup>*cnF*12(*x*) *R* dΦ*y*(*x*) d*x* − *c* 2 *n*Φ*y*(*x*) *R* d*G*12(*x*) d*x* <sup>+</sup>2*cn*Φ*y*(*x*) *R* d*F*12(*x*) d*x* − *c* 2 *nG*12(*x*) *R* dΦ*y*(*x*) d*x* = <sup>1</sup> Φ*x*(*t*) (*I*0(*x*) − *cI*3(*x*))*U*(*x*) *d*2*U*(*t*) *dt*<sup>2</sup> <sup>+</sup> *I*2(*x*) − 2*cI*4(*x*) + *c* 2 *I*6(*x*) Φ*x*(*x*) *d*2Φ*x*(*t*) *dt*<sup>2</sup> <sup>−</sup>*c*(*I*4(*x*) <sup>−</sup> *cI*6(*x*)) *dW*(*x*) *dx d*2*W*(*t*) *dt*<sup>2</sup> , (A4)

*nB*12(*x*) *R* d*U*(*x*) d*x* + *nU*(*x*) *R* d*B*66(*x*) d*x* + *nB*66(*x*) *R* d*U*(*x*) d*x* −*cnD*12(*x*) *R* d*U*(*x*) d*x* <sup>−</sup> *cnU*(*x*) *R* d*D*66(*x*) d*x* <sup>−</sup> *cnD*66(*x*) *R* d*U*(*x*) d*x* + d*B*66(*x*) d*x* d*V*(*x*) d*x* + *B*66(*x*) d2*V*(*x*) d*x*<sup>2</sup> <sup>−</sup> *<sup>n</sup>*2*B*11(*x*)*V*(*x*) *R*2 <sup>+</sup>*cn*2*D*11(*x*)*V*(*x*) *<sup>R</sup>*<sup>2</sup> − *c* d*D*66(*x*) d*x* d*V*(*x*) d*x* − *cD*66(*x*) d2*V*(*x*) d*x*<sup>2</sup> +*nB*11(*x*)*W*(*x*) *<sup>R</sup>*<sup>2</sup> <sup>−</sup> *nA*66(*x*)*W*(*x*) *<sup>R</sup>* <sup>−</sup> *cnF*12(*x*) *R* d2*W*(*x*) d*x*<sup>2</sup> <sup>−</sup> <sup>9</sup>*<sup>c</sup>* 2 *nF*66(*x*)*W*(*x*) *R* −2*cnF*66(*x*) *R* d2*W*(*x*) d*x*<sup>2</sup> <sup>+</sup> *cn*3*F*11(*x*)*W*(*x*) *<sup>R</sup>*<sup>3</sup> <sup>−</sup> *<sup>c</sup>* 2 *n*3*G*11(*x*)*W*(*x*) *R*3 <sup>−</sup>2*cn R* d*F*66(*x*) d*x* d*W*(*x*) d*x* + <sup>2</sup>*<sup>c</sup>* 2 *nG*66(*x*) *R* d2*W*(*x*) d*x*<sup>2</sup> +6*cnC*66(*x*)*W*(*x*) *<sup>R</sup>* <sup>+</sup> <sup>2</sup>*<sup>c</sup>* 2 *n a* d*G*66(*x*) d*x* d*W*(*x*) d*x* + *<sup>c</sup>* 2 *nG*12(*x*) *R* d2*W*(*x*) d*x*<sup>2</sup> −*cnD*11(*x*)*W*(*x*) *<sup>R</sup>*<sup>2</sup> <sup>+</sup> *nC*12(*x*) *R* dΦ*x*(*x*) d*x* +*n*Φ*x*(*x*) *R* d*C*66(*x*) d*x* + *nC*66(*x*) *R* dΦ*x*(*x*) d*x* + *<sup>c</sup>* 2 *n*Φ*x*(*x*) *R* d*G*66(*x*) d*x* −2*cnF*12(*x*) *R* dΦ*x*(*x*) d*x* <sup>−</sup> <sup>2</sup>*cnF*66(*x*) *R* dΦ*x*(*x*) d*x* + *<sup>c</sup>* 2 *nG*66(*x*) *R* dΦ*x*(*x*) d*x* +*c* 2 *nG*12(*x*) *R* dΦ*x*(*x*) d*x* <sup>−</sup> <sup>2</sup>*cn*Φ*x*(*x*) *R* d*F*66(*x*) d*x* + d*C*66(*x*) d*x* dΦ*y*(*x*) d*x* + *C*66(*x*) d2Φ*y*(*x*) d*x*<sup>2</sup> <sup>−</sup> *<sup>A</sup>*66(*x*)Φ*y*(*x*) <sup>−</sup> *<sup>n</sup>*2*C*11(*x*)Φ*y*(*x*) *R*2 <sup>+</sup>2*cn*2*F*11(*x*)Φ*y*(*x*) *<sup>R</sup>*<sup>2</sup> <sup>−</sup> *<sup>c</sup>* 2 *n*2*G*11(*x*)Φ*y*(*x*) *<sup>R</sup>*<sup>2</sup> − 2*c* d*F*66(*x*) d*x* dΦ*y*(*x*) d*x* −2*cF*66(*x*) d2Φ*y*(*x*) d*x*<sup>2</sup> + 6*cC*66(*x*)Φ*y*(*x*) − 9*c* 2 *F*66(*x*)Φ*y*(*x*) +*c* 2 d*G*66(*x*) d*x* dΦ*y*(*x*) d*x* + *c* 2 *G*66(*x*) d2Φ*y*(*x*) d*x*<sup>2</sup> = <sup>1</sup> Φ*y*(*t*) (*I*0(*x*) − *cI*3(*x*))*V*(*x*) *d*2*V*(*t*) *dt*<sup>2</sup> <sup>+</sup> *I*2(*x*) − 2*cI*4(*x*) + *c* 2 *I*6(*x*) Φ*y*(*x*) *d*2Φ*y*(*t*) *dt*<sup>2</sup> <sup>−</sup>*c*(*I*4(*x*) <sup>−</sup> *cI*6(*x*)) *<sup>n</sup> R d*2*W*(*t*) *dt*<sup>2</sup> , (A5)

where

$$\left(A\_{i\bar{j}}(\mathbf{x}), B\_{i\bar{j}}(\mathbf{x}), \mathbb{C}\_{i\bar{j}}(\mathbf{x}), D\_{i\bar{j}}(\mathbf{x}), F\_{i\bar{j}}(\mathbf{x}), \mathbb{G}\_{i\bar{j}}(\mathbf{x})\right) = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} Q\_{\bar{i}\bar{j}}(\mathbf{x}, z) \{1, z, z^2, z^3, z^4, z^6\} dz \tag{A6}$$

#### **Appendix B**

In what follows we rewrite the equations of motion (A1)–(A5) in a discretized form, according to the GDQ method.


+*B*11 *xp* - *N r*=1 χ(2) *pr* Φ*x*(*xr*) <sup>−</sup> *cD*11 *xp* - *N r*=1 χ(2) *pr* Φ*x*(*xr*) −*c* - *N r*=1 χ(1) *pr D*11(*xr*) - *N r*=1 χ(1) *pr* Φ*x*(*xr*) <sup>−</sup> *<sup>n</sup>*2*B*66(*xp*)Φ*x*(*xp*) *<sup>R</sup>*<sup>2</sup> <sup>+</sup> *cn*2*D*66(*xp*)Φ*x*(*xp*) *R*2 <sup>−</sup>*n*Φ*y*(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr B*12(*xr*) <sup>−</sup> *nB*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr* Φ*y*(*xr*) <sup>−</sup> *nB*66(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr* Φ*y*(*xr*) <sup>+</sup>*cnD*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr* Φ*y*(*xr*) <sup>+</sup> *cn*Φ*y*(*xp*) *<sup>a</sup>* - *N r*=1 χ(1) *pr D*12(*xr*) <sup>+</sup> *cnD*66(*xp*) *<sup>a</sup>* - *N r*=1 χ(1) *pr* Φ*y*(*xr*) = <sup>1</sup> *U*(*t*) *I*0 *xp U xp <sup>d</sup>*2*U*(*t*) *dt*<sup>2</sup> <sup>+</sup> *I*0 *xp* − *cI*<sup>3</sup> *xp* Φ*<sup>x</sup> xp <sup>d</sup>*2Φ*x*(*t*) *dt*<sup>2</sup> −*cI*<sup>3</sup> *xp* - *N r*=1 χ(1) *pr W*(*xr*) *d*2*W*(*t*) *dt*<sup>2</sup> *nU*(*xp*) *R* - *N r*=1 χ(1) *pr A*66(*xr*) <sup>+</sup> *nA*66(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr U*(*xr*) <sup>+</sup> *nA*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr U*(*xr*) + - *N r*=1 χ(1) *pr A*66(*xr*) - *N r*=1 χ(1) *pr V*(*xr*) + *A*66 *xp* - *N r*=1 χ(2) *pr V*(*xr*) <sup>−</sup> *<sup>n</sup>*2*A*11(*xp*)*V*(*xp*) *R*2 <sup>−</sup>2*cn R* - *N r*=1 χ(1) *pr D*66(*xr*) - *N r*=1 χ(1) *pr W*(*xr*) <sup>−</sup> <sup>2</sup>*cnD*66(*xp*) *<sup>R</sup>* - *N r*=1 χ(2) *pr W*(*xr*) <sup>+</sup> *nA*11(*xp*)*W*(*xp*) *R*2 <sup>+</sup>*cn*3*D*11(*xp*)*W*(*xp*) *<sup>R</sup>*<sup>3</sup> <sup>−</sup> *cnD*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(2) *pr W*(*xr*) <sup>+</sup>*n*Φ*x*(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr B*66(*xr*) <sup>+</sup> *nB*66(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr* Φ*x*(*xr*) <sup>−</sup> *cn*Φ*x*(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr D*66(*xr*) <sup>−</sup>*cnD*66(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr* Φ*x*(*xr*) <sup>+</sup> *nB*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr* Φ*x*(*xr*) <sup>−</sup>*cnD*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr* Φ*x*(*xr*) + - *N r*=1 χ(1) *pr B*66(*xr*) - *N r*=1 χ(1) *pr* Φ*y*(*xr*) +*B*66 *xp* - *N r*=1 χ(2) *pr* Φ*y*(*xr*) − *c* - *N r*=1 χ(1) *pr D*66(*xr*) - *N r*=1 χ(1) *pr* Φ*y*(*xr*) <sup>−</sup>*cD*66 *xp* - *N r*=1 χ(2) *pr* Φ*y*(*xr*) <sup>−</sup> *<sup>n</sup>*2*B*11(*xp*)Φ*y*(*xp*) *<sup>R</sup>*<sup>2</sup> <sup>+</sup> *cn*2*D*11(*xp*)Φ*y*(*xp*) *R*2 = <sup>1</sup> *V*(*t*) *I*0(*x*)*V xp <sup>d</sup>*2*V*(*t*) *dt*<sup>2</sup> <sup>+</sup> *I*0 *xp* − *cI*<sup>3</sup> *xp* Φ*<sup>y</sup> xp <sup>d</sup>*2Φ*y*(*t*) *dt*<sup>2</sup> <sup>−</sup>*cn <sup>R</sup> I*<sup>3</sup> *xp W xp <sup>d</sup>*2*W*(*t*) *dt*<sup>2</sup> (A8) −2*cn*2*U*(*xp*) *R*2 - *N r*=1 χ(1) *pr D*66(*xr*) <sup>−</sup> <sup>2</sup>*cn*2*D*66(*xp*) *R*2 - *N r*=1 χ(1) *pr U*(*xr*) −*cn*2*D*12(*xp*) *R*2 - *N r*=1 χ(1) *pr U*(*xr*) <sup>−</sup> *<sup>A</sup>*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr U*(*xr*) +*c* - *N r*=1 χ(2) *pr D*11(*xr*) - *N r*=1 χ(1) *pr U*(*xr*) + 2*c* - *N r*=1 χ(1) *pr D*11(*xr*) - *N r*=1 χ(2) *pr U*(*xr*) +*cD*11 *xp* - *N r*=1 χ(3) *pr U*(*xr*) <sup>−</sup> *cnV*(*xp*) *<sup>R</sup>* - *N r*=1 χ(2) *pr D*12(*xr*) <sup>−</sup>2*cn R* - *N r*=1 χ(1) *pr D*12(*xr*) - *N r*=1 χ(1) *pr V*(*xr*) <sup>−</sup> *cnD*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(2) *pr V*(*xr*) - - - 

$$\begin{split} & -\frac{2\tilde{c}n}{\mathbb{R}} \Big( \sum\_{r=1}^{N} \chi\_{pr}^{(1)} D\_{\theta\theta}(\mathbf{x}\_{r}) \Big) \Big( \sum\_{r=1}^{N} \chi\_{pr}^{(1)} V(\mathbf{x}\_{r}) \Big) - \frac{2\tilde{c}n D\_{\theta\theta} \Big( \mathbf{x}\_{p} \big)}{\mathbb{R}} \Big( \sum\_{r=1}^{N} \chi\_{pr}^{(2)} V(\mathbf{x}\_{r}) \Big) \\ & + \frac{2\mathbf{n}^{3} D\_{11} \big( \mathbf{x}\_{p} \big) V(\mathbf{x}\_{p})}{R^{3}} + \frac{n A\_{11} \big( \mathbf{x}\_{p} \big) V(\mathbf{x}\_{p})}{R^{2}} + \frac{2\tilde{c}^{2} n^{2}}{R^{2}} \Big( \sum\_{r=1}^{N} \chi\_{pr}^{(1)} G\_{12}(\mathbf{x}\_{r}) \Big) \Big( \sum\_{r=1}^{N} \chi\_{pr}^{(1)} W(\mathbf{x}\_{r}) \Big) \Big) \end{split}$$

+4*<sup>c</sup>* 2 *n*2 *R*2 - *N r*=1 χ(1) *pr G*66(*xr*) - *N r*=1 χ(1) *pr W*(*xr*) − *c* 2 *n*4*G*11(*xp*)*W*(*xp*) *R*4 +*c* 2 *n*2*W*(*xp*) *R*2 - *N r*=1 χ(2) *pr G*12(*xr*) + <sup>2</sup>*<sup>c</sup>* 2 *n*2*G*12(*xp*) *R*2 - *N r*=1 χ(2) *pr W*(*xr*) − 9*c* 2 *n*2*F*66(*xp*)*W*(*xp*) *R*2 +4*<sup>c</sup>* 2 *n*2*G*66(*xp*) *R*2 - *N r*=1 χ(2) *pr W*(*xr*) <sup>+</sup> <sup>6</sup>*cn*2*C*66(*xp*)*W*(*xp*) *<sup>R</sup>*<sup>2</sup> <sup>−</sup> <sup>2</sup>*cn*2*D*11(*xp*)*W*(*xp*) *R*3 −*n*2*A*66(*xp*)*W*(*xp*) *<sup>R</sup>*<sup>2</sup> <sup>+</sup> *cW*(*xp*) *<sup>R</sup>* - *N r*=1 χ(2) *pr D*12(*xr*) + <sup>2</sup>*<sup>c</sup> R* - *N r*=1 χ(1) *pr D*12(*xr*) - *N r*=1 χ(1) *pr W*(*xr*) <sup>+</sup>2*cD*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(2) *pr W*(*xr*) + *A*66 *xp* - *N r*=1 χ(2) *pr W*(*xr*) + - *N r*=1 χ(1) *pr A*66(*xr*) - *N r*=1 χ(1) *pr W*(*xr*) −2*c* 2 - *N r*=1 χ(1) *pr G*11(*xr*) - *N r*=1 χ(3) *pr W*(*xr*) <sup>−</sup> *<sup>c</sup>*2*G*11 *xp* - *N r*=1 χ(4) *pr W*(*xr*) <sup>−</sup> *<sup>A</sup>*11(*xp*)*W*(*xp*) *R*2 −6*c* - *N r*=1 χ(1) *pr C*66(*xr*) - *N r*=1 χ(1) *pr W*(*xr*) <sup>−</sup> <sup>6</sup>*cC*66 *xp* - *N r*=1 χ(2) *pr W*(*xr*) +2*<sup>c</sup>* 2 *n*2*G*66(*xp*) *R*2 - *N r*=1 χ(1) *pr* Φ*x*(*xr*) + <sup>2</sup>*<sup>c</sup>* 2 *n*2Φ*x*(*xp*) *R*2 - *N r*=1 χ(1) *pr G*66(*xr*) −*cn*2*F*12(*xp*) *R*2 - *N r*=1 χ(1) *pr* Φ*x*(*xr*) <sup>−</sup> <sup>2</sup>*cn*2Φ*x*(*xp*) *R*2 - *N r*=1 χ(1) *pr F*66(*xr*) −2*cn*2*F*66(*xp*) *R*2 - *N r*=1 χ(1) *pr* Φ*x*(*xr*) + *c* 2 *n*2*G*12(*xp*) *R*2 - *N r*=1 χ(1) *pr* Φ*x*(*xr*) <sup>+</sup>*cD*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr* Φ*<sup>x</sup> xp* + - *N r*=1 χ(1) *pr A*66(*xr*) Φ*<sup>x</sup> xp* + *A*66 *xp* - *N r*=1 χ(1) *pr* Φ*x*(*xr*) −*c* 2 *G*11 *xp* - *N r*=1 χ(3) *pr* Φ*x*(*xr*) <sup>−</sup> *<sup>B</sup>*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr* Φ*x*(*xr*) − 6*c* - *N r*=1 χ(1) *pr C*66(*xr*) Φ*<sup>x</sup> xp* <sup>−</sup>6*cC*66 *xp* - *N r*=1 χ(1) *pr* Φ*x*(*xr*) + 9*c* 2 - *N r*=1 χ(1) *pr F*66(*xr*) Φ*<sup>x</sup> xp* +9*c* 2 *F*66 *xp* - *N r*=1 χ(1) *pr* Φ*x*(*xr*) + *cF*11(*x*) - *N r*=1 χ(3) *pr* Φ*x*(*xr*) +*c* - *N r*=1 χ(2) *pr F*11(*xr*) - *N r*=1 χ(1) *pr* Φ*x*(*xr*) + 2*c* - *N r*=1 χ(1) *pr F*11(*xr*) - *N r*=1 χ(2) *pr* Φ*x*(*xr*) −*c* 2 - *N r*=1 χ(2) *pr G*11(*xr*) - *N r*=1 χ(1) *pr* Φ*x*(*xr*) − 2*c* 2 - *N r*=1 χ(1) *pr G*11(*xr*) - *N r*=1 χ(2) *pr* Φ*x*(*xr*) +2*<sup>c</sup>* 2 *n R* - *N r*=1 χ(1) *pr G*12(*xr*) - *N r*=1 χ(1) *pr* Φ*y*(*xr*) − *c* 2 *n*3*G*11(*xp*)Φ*y*(*xp*) *R*3 −9*c* 2 *nF*66(*xp*)Φ*y*(*xp*) *<sup>R</sup>* <sup>+</sup> <sup>2</sup>*<sup>c</sup>* 2 *n R* - *N r*=1 χ(1) *pr G*66(*xr*) - *N r*=1 χ(1) *pr* Φ*y*(*xr*) +*c* 2 *nG*12(*xp*) *R* - *N r*=1 χ(2) *pr* Φ*y*(*xr*) <sup>+</sup> *cn*3*F*11(*xp*)Φ*y*(*xp*) *R*3 −*cnD*11(*xp*)Φ*y*(*xp*) *<sup>R</sup>*<sup>2</sup> <sup>+</sup> <sup>2</sup>*<sup>c</sup>* 2 *nG*66(*xp*) *R* - *N r*=1 χ(2) *pr* Φ*y*(*xr*) <sup>+</sup>6*cnC*66(*xp*)Φ*y*(*xp*) *<sup>R</sup>* <sup>−</sup> <sup>2</sup>*cn R* - *N r*=1 χ(1) *pr F*66(*xr*) - *N r*=1 χ(1) *pr* Φ*y*(*xr*) <sup>−</sup>2*cnF*66(*xp*) *<sup>R</sup>* - *N r*=1 χ(2) *pr* Φ*y*(*xr*) <sup>−</sup> *cn*Φ*y*(*xp*) *<sup>R</sup>* - *N r*=1 χ(2) *pr F*12(*xr*) <sup>−</sup>2*cn R* - *N r*=1 χ(1) *pr F*12(*xr*) - *N r*=1 χ(1) *pr* Φ*y*(*xr*) <sup>−</sup> *cnF*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(2) *pr* Φ*y*(*xr*) +*c* 2 *n*Φ*y*(*xp*) *R* - *N r*=1 χ(2) *pr G*12(*xr*) <sup>+</sup> *nB*11(*xp*)Φ*y*(*xp*) *<sup>R</sup>*<sup>2</sup> <sup>−</sup> *nA*66(*xp*)Φ*y*(*xp*) *<sup>R</sup>* <sup>−</sup> *kwW xp* +*kg* - *N r*=1 χ(2) *pr <sup>W</sup>*(*xr*) <sup>−</sup> *<sup>n</sup>*<sup>2</sup> *<sup>R</sup>*<sup>2</sup> *<sup>W</sup> xp* − *Fa N r*=1 χ(2) *pr W*(*xr*) = <sup>1</sup> *W*(*t*) *I*0 *xp W xp <sup>d</sup>*2*W*(*t*) *dt*<sup>2</sup> (A9)

$$\begin{split} -\tilde{\mathsf{c}}^{2}I\_{6}\Big(\mathbf{x}\_{p}\Big)\Big(\sum\_{r=1}^{N}\chi\_{pr}^{(2)}\big)\mathcal{W}(\mathbf{x}\_{r})\frac{d^{2}\overline{\mathsf{W}}(t)}{dt^{2}} - \frac{\mathsf{n}^{2}}{R^{2}}\mathcal{W}\Big(\mathbf{x}\_{p}\Big)\frac{d^{2}\overline{\mathsf{W}}(t)}{dt^{2}}\Big) + \tilde{\mathsf{c}}I\_{3}\Big(\mathbf{x}\_{p}\Big)\Big(\sum\_{r=1}^{N}\chi\_{pr}^{(1)}\big)\mathcal{U}(\mathbf{x}\_{r})\frac{d^{2}\overline{\mathsf{W}}(t)}{dt^{2}} \\ - \frac{\mathsf{n}}{R}\mathcal{V}\Big(\mathbf{x}\_{p}\Big)\frac{d^{2}\overline{\mathcal{V}}(t)}{dt^{2}}\Big) + \tilde{\mathsf{c}}\Big(I\_{4}\big(\mathbf{x}\_{p}\big) - \mathsf{\tilde{c}}I\_{6}\big(\mathbf{x}\_{p}\big)\Big)\Big(\sum\_{r=1}^{N}\chi\_{pr}^{(1)}\Phi\_{\mathbf{x}}\big)\frac{d^{2}\overline{\Phi}\_{\mathbf{x}}(t)}{dt^{2}} \\ - \frac{\mathsf{n}}{R}\Phi\_{\mathbf{y}}\Big(\mathbf{x}\_{p}\Big)\frac{d^{2}\overline{\Phi}\_{\mathbf{y}}(t)}{dt^{2}}\Big) \end{split}$$

−*n*2*B*66(*xp*)*U*(*xp*) *<sup>R</sup>*<sup>2</sup> <sup>+</sup> *cn*2*D*66(*xp*)*U*(*xp*) *<sup>R</sup>*<sup>2</sup> + - *N r*=1 χ(1) *pr B*11(*xr*) - *N r*=1 χ(1) *pr U*(*xr*) +*B*11 *xp* - *N r*=1 χ(2) *pr U*(*xr*) − *c* - *N r*=1 χ(1) *pr D*11(*xr*) - *N r*=1 χ(1) *pr U*(*xr*) <sup>−</sup> *cD*11 *xp* - *N r*=1 χ(2) *pr U*(*xr*) <sup>−</sup>*nB*66(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr V*(*xr*) <sup>−</sup> *nV*(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr B*12(*xr*) <sup>−</sup> *nB*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr V*(*xr*) <sup>+</sup>*cnD*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr V*(*xr*) <sup>+</sup> *cnD*66(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr V*(*xr*) <sup>+</sup> *cnV*(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr D*12(*xr*) <sup>−</sup>*cW*(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr D*12(*xr*) <sup>−</sup> *cD*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr W*(*xr*) <sup>−</sup> *<sup>c</sup>*2*n*2*G*12(*xp*) *R*2 - *N r*=1 χ(1) *pr W*(*xr*) <sup>+</sup>*cn*2*F*12(*xp*) *R*2 - *N r*=1 χ(1) *pr W*(*xr*) <sup>+</sup> *cn*2*W*(*xp*) *R*2 - *N r*=1 χ(1) *pr F*12(*xr*) −2*c* 2 *n*2*G*66(*xp*) *R*2 - *N r*=1 χ(1) *pr W*(*xr*) <sup>+</sup> <sup>2</sup>*cn*2*F*66(*xp*) *R*2 - *N r*=1 χ(1) *pr W*(*xr*) −*c*2*n*2*W*(*xp*) *R*2 - *N r*=1 χ(1) *pr G*12(*xr*) <sup>−</sup> *<sup>A</sup>*66 *xp* - *N r*=1 χ(1) *pr W*(*xr*) − 9*c* 2 *F*66 *xp* - *N r*=1 χ(1) *pr W*(*xr*) +*c* 2 - *N r*=1 χ(1) *pr G*11(*xr*) - *N r*=1 χ(2) *pr W*(*xr*) + *c* 2 *G*11 *xp* - *N r*=1 χ(3) *pr W*(*xr*) <sup>+</sup> *<sup>W</sup>*(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr B*12(*xr*) <sup>+</sup> *<sup>B</sup>*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr W*(*xr*) − *c* - *N r*=1 χ(1) *pr F*11(*xr*) - *N r*=1 χ(2) *pr W*(*xr*) <sup>−</sup>*cF*11 *xp* - *N r*=1 χ(3) *pr W*(*xr*) + 6*cC*66 *xp* - *N r*=1 χ(1) *pr W*(*xr*) −*n*2*C*66(*xp*)Φ*x*(*xp*) *<sup>R</sup>*<sup>2</sup> <sup>+</sup> <sup>2</sup>*cn*2*F*66(*xp*)Φ*x*(*xp*) *<sup>R</sup>*<sup>2</sup> <sup>−</sup> *<sup>c</sup>* 2 *n*2*G*66(*xp*)Φ*x*(*xp*) *R*2 + - *N r*=1 χ(1) *pr C*11(*xr*) - *N r*=1 χ(1) *pr* Φ*x*(*xr*) + *C*11 *xp* - *N r*=1 χ(2) *pr* Φ*x*(*xr*) <sup>−</sup> *<sup>A</sup>*66 *xp* Φ*<sup>x</sup> xp* −9*c* 2 *F*66 *xp* Φ*<sup>x</sup> xp* + *c* 2 - *N r*=1 χ(1) *pr G*11(*xr*) - *N r*=1 χ(1) *pr* Φ*x*(*xr*) +*c* 2 *G*11 *xp* - *N r*=1 χ(2) *pr* Φ*x*(*xr*) − 2*c* - *N r*=1 χ(1) *pr F*11(*xr*) - *N r*=1 χ(1) *pr* Φ*x*(*xr*) <sup>−</sup>2*cF*11 *xp* - *N r*=1 χ(2) *pr* Φ*x*(*xr*) + 6*cC*66 *xp* Φ*<sup>x</sup> xp* <sup>−</sup> *nC*66(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr* Φ*y*(*xr*) <sup>−</sup>*n*Φ*y*(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr C*12(*xr*) <sup>−</sup> *nC*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr* Φ*y*(*xr*) (A10)

$$\llcorner$$

<sup>+</sup> <sup>2</sup>*cnF*66(*xp*) *<sup>R</sup>*



*N r*=1 χ(1) *pr* Φ*y*(*xr*) 

*N r*=1 χ(1) *pr G*12(*xr*)

*N r*=1 χ(1) *pr* Φ*y*(*xr*) 

*d*2*W*(*t*) *dt*<sup>2</sup> ,

−*c* 2 *nG*66(*xp*) *a*

<sup>+</sup>2*cnF*12(*xp*) *<sup>R</sup>*

<sup>+</sup>2*cn*Φ*y*(*xp*) *<sup>R</sup>*

= <sup>1</sup> Φ*x*(*t*) *I*0 *xp* − *cI*<sup>3</sup> *xp U xp <sup>d</sup>*2*U*(*t*) *dt*<sup>2</sup> <sup>+</sup> *I*2 *xp* − 2*cI*<sup>4</sup> *xp* 

+*c* 2 *I*6 *xp* Φ*<sup>x</sup> xp <sup>d</sup>*2Φ*x*(*t*) *dt*<sup>2</sup> − *c I*4 *xp* − *cI*<sup>6</sup> *xp <sup>N</sup> r*=1 χ(1) *pr W*(*xr*)




*N r*=1 χ(1) *pr* Φ*y*(*xr*) 

*N r*=1 χ(1) *pr* Φ*y*(*xr*) − *c* 2 *n*Φ*y*(*xp*) *R* -

*N r*=1 χ(1) *pr F*12(*xr*) − *c* 2 *nG*12(*xp*) *R*

*nB*12(*xp*) *R* - *N r*=1 χ(1) *pr U*(*xr*) <sup>+</sup> *nU*(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr B*66(*xr*) <sup>+</sup> *nB*66(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr U*(*xr*) <sup>−</sup>*cnD*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr U*(*xr*) <sup>−</sup> *cnU*(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr D*66(*xr*) <sup>−</sup> *cnD*66(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr U*(*xr*) + - *N r*=1 χ(1) *pr B*66(*xr*) - *N r*=1 χ(1) *pr V*(*xr*) + *B*66 *xp* - *N r*=1 χ(2) *pr V*(*xr*) <sup>−</sup> *<sup>n</sup>*2*B*11(*xp*)*V*(*xp*) *R*2 <sup>+</sup>*cn*2*D*11(*xp*)*V*(*xp*) *<sup>R</sup>*<sup>2</sup> − *c* - *N r*=1 χ(1) *pr D*66(*xr*) - *N r*=1 χ(1) *pr V*(*xr*) <sup>−</sup> *cD*66 *xp* - *N r*=1 χ(2) *pr V*(*xr*) <sup>+</sup>*nB*11(*xp*)*W*(*xp*) *<sup>R</sup>*<sup>2</sup> <sup>−</sup> *nA*66(*xp*)*W*(*xp*) *<sup>R</sup>* <sup>−</sup> *cnF*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(2) *pr W*(*xr*) − 9*c* 2 *nF*66(*xp*)*W*(*xp*) *R* <sup>−</sup>2*cnF*66(*xp*) *<sup>R</sup>* - *N r*=1 χ(2) *pr W*(*xr*) <sup>+</sup> *cn*3*F*11(*xp*)*W*(*xp*) *<sup>R</sup>*<sup>3</sup> <sup>−</sup> *<sup>c</sup>* 2 *n*3*G*11(*xp*)*W*(*xp*) *R*3 <sup>−</sup>2*cn R* - *N r*=1 χ(1) *pr F*66(*xr*) - *N r*=1 χ(1) *pr W*(*xr*) + <sup>2</sup>*<sup>c</sup>* 2 *nG*66(*xp*) *R* - *N r*=1 χ(2) *pr W*(*xr*) <sup>+</sup>6*cnC*66(*xp*)*W*(*xp*) *<sup>R</sup>* <sup>+</sup> <sup>2</sup>*<sup>c</sup>* 2 *n a* - *N r*=1 χ(1) *pr G*66(*xr*) - *N r*=1 χ(1) *pr W*(*xr*) +*c* 2 *nG*12(*xp*) *R* - *N r*=1 χ(2) *pr W*(*xr*) <sup>−</sup> *cnD*11(*xp*)*W*(*xp*) *<sup>R</sup>*<sup>2</sup> <sup>+</sup> *nC*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr* Φ*x*(*xr*) <sup>+</sup>*n*Φ*x*(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr C*66(*xr*) <sup>+</sup> *nC*66(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr* Φ*x*(*xr*) +*c* 2 *n*Φ*x*(*xp*) *R* - *N r*=1 χ(1) *pr G*66(*xr*) <sup>−</sup> <sup>2</sup>*cnF*12(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr* Φ*x*(*xr*) <sup>−</sup>2*cnF*66(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr* Φ*x*(*xr*) + *c* 2 *nG*66(*xp*) *R* - *N r*=1 χ(1) *pr* Φ*x*(*xr*) +*c* 2 *nG*12(*xp*) *R* - *N r*=1 χ(1) *pr* Φ*x*(*xr*) <sup>−</sup> <sup>2</sup>*cn*Φ*x*(*xp*) *<sup>R</sup>* - *N r*=1 χ(1) *pr F*66(*xr*) + - *N r*=1 χ(1) *pr C*66(*xr*) - *N r*=1 χ(1) *pr* Φ*y*(*xr*) + *C*66 *xp* - *N r*=1 χ(2) *pr* Φ*y*(*xr*) <sup>−</sup>*A*66 *xp* Φ*<sup>y</sup> xp* <sup>−</sup> *<sup>n</sup>*2*C*11(*xp*)Φ*y*(*xp*) *<sup>R</sup>*<sup>2</sup> <sup>+</sup> <sup>2</sup>*cn*2*F*11(*xp*)Φ*y*(*xp*) *R*2 −*c* 2 *n*2*G*11(*xp*)Φ*y*(*xp*) *<sup>R</sup>*<sup>2</sup> − 2*c* - *N r*=1 χ(1) *pr F*66(*xr*) - *N r*=1 χ(1) *pr* Φ*y*(*xr*) <sup>−</sup>2*cF*66 *xp* - *N r*=1 χ(2) *pr* Φ*y*(*xr*) + 6*cC*66 *xp* Φ*<sup>y</sup> xp* − 9*c* 2 *F*66 *xp* Φ*<sup>y</sup> xp* +*c* 2 - *N r*=1 χ(1) *pr G*66(*xr*) - *N r*=1 χ(1) *pr* Φ*y*(*xr*) + *c* 2 *G*66 *xp* - *N r*=1 χ(2) *pr* Φ*y*(*xr*) = <sup>1</sup> Φ*y*(*t*) *I*0 *xp* − *cI*<sup>3</sup> *xp V xp <sup>d</sup>*2*V*(*t*) *dt*<sup>2</sup> <sup>+</sup> *I*2 *xp* − 2*cI*<sup>4</sup> *xp* + *c* 2 *I*6 *xp* Φ*<sup>y</sup> xp <sup>d</sup>*2Φ*y*(*t*) *dt*<sup>2</sup> −*c I*4 *xp* − *cI*<sup>6</sup> *xp n R d*2*W*(*t*) *dt*<sup>2</sup> , (A11)

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
