Dynamic Stability of Bi-Directional Functionally Graded Porous Cylindrical Shells Embedded in an Elastic Foundation

This paper investigates the dynamic buckling of bi-directional (BD) functionally graded (FG) porous cylindrical shells for various boundary conditions, where the FG material is modeled by means of power law functions with even and uneven porosity distributions of ceramic and metal phases. The third-order shear deformation theory (TSDT) is adopted to derive the governing equations of the problem via the Hamilton’s principle. The generalized differential quadrature (GDQ) method is applied together with the Bolotin scheme as numerical strategy to solve the problem, and to draw the dynamic instability region (DIR) of the structure. A large parametric study examines the effect of different boundary conditions at the extremities of the cylindrical shell, as well as the sensitivity of the dynamic stability to different thickness-to-radius ratios, length-to-radius ratios, transverse and longitudinal power indexes, porosity volume fractions, and elastic foundation constants. Based on results, the dynamic stability of BD-FG cylindrical shells can be controlled efficiently by selecting appropriate power indexes along the desired directions. Furthermore, the DIR is highly sensitive to the porosity distribution and to the extent of transverse and longitudinal power indexes. The numerical results could be of great interest for many practical applications, as civil, mechanical or aerospace engineering, as well as for energy devices or biomedical systems.


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][2][3][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][14][15][16][17][18][19][20][21][22][23][24][25][26][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][29][30][31][32][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][42][43][44][45][46][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-toradius 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.

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 where z n and x n 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., for even distributions, or 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.
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 )   0  3  0   ,,  , , , , , , , , , , xx w x y t u x y z t u x y t z x y t cz x y t x  According to the TSDT, the strain components of the cylindrical shell can be written as [48] ( where the stiffness ij Q is defined as 21 11 , , , , , , The strain energy of the BD-FG porous cylindrical shell is expressed as follow xx xx yy yy xy xy xz xz yz yz h dzdydx (8) and the kinetic energy of the cylindrical shell reads The total potential energy corresponding to the axial compressive load where w k and g k refer to the Winkler foundation stiffness and shear layer stiffness of the elastic foundation, respectively.

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 where n is the circumferential half wave number, 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][52][53][54][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.

The GDQ Method
The GDQ method approximates the fundamental system of differential equations, by discretizing the derivatives of a function ( ) Jx 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 xL  is discretized in N grids points, the approximation of the nth-order derivatives of J function respect to x variable can be expressed as [49] ( )  being the weighting coefficients, defined as follows [49] and ( ) p x  is the Lagrangian operator expressed as [49] ( ) 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] 1 1 cos 1, 2,..., . 21 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 where  and  refer to the static and dynamic load factors, respectively. Furthermore, cr F 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 where M , K , and

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]       where s  and s  denote the arbitrary time invariant vectors. It should be mentioned that the first DIR with period 2T 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 2T . By substitution of Equation (36) into Equation (35) and by mathematical manipulation, we get the following first order equation (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 ( ) Equation (37) for some fixed values of  and cr F , the variation of the excitation frequency  in regards to  can be drawn as DIR for the BD-FG structure.

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.

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 x n . In Table 1 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.
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.

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 We determine the DIR, and highlight the effects of different parameters such as the thicknessto-radius ratio ( / hR ), the length-to-radius ratio ( / LR ), the static load factor, the boundary conditions, the power law indexes ( x n , z n ), 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 ( / hR ), where a clear shift of the DIR is observed for increasing / hR 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 / hR 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 / LR ratios, while keeping fixed / 0.01 hR = . In detail, for /1 LR = the structure has a wider DIR in comparison with the other values, whereby an increasing / LR 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 / LR ranging between 1 and 10. When the / LR ratio features higher magnitudes, the bending resistance gradually reduces and yields an increased bending deformation.
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. 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.
A further investigation is devoted to study the influence of the power law index along the length x n 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.   sensitive to the dynamic instability, due to their higher stiffness. The contrary occurs for higher values of x n and z n 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. 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   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 x n , z n , for different values of  . Figure 14 is devoted to check for the influence of  on the DIR of the structure for 0.15 xz nn == . 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 1.5 xz nn ==, 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.
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.

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 thicknessto-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

•
An increased thickness-to-radius ratio causes a general shift of the DIR origin towards higher excitation frequencies. Moreover, the DIR gets wider at a certain value of the dynamic load factor.

•
An increased length-to-radius dimensionless ratio moves the DIR origin towards lower excitation frequencies, whereas the DIR gets smaller. • A simultaneous increase of longitudinal and transverse power indexes yields an overall decrease in the excitation frequencies associated with the DIR origin.

•
The control of the dynamic instability for a BD-FG cylindrical shell, is convenient for an appropriate selection of the power indexes. • BD-FG cylindrical shells with a simply support at both ends are more unstable than the clampedclamped or clamped-simply supported structures, because of their higher deformability.

•
The effect of coefficients and type of porosity on the structural DIR depend on the extent of the longitudinal and transverse power law indexes. There exists a certain value for these indexes, for which the excitation frequencies corresponding to the DIR can invert their behavior.
• A general increase in the elastic foundation coefficients yields higher excitation structural frequencies especially when a Pasternak foundation is assumed instead of a Winkler foundation. Anyway, the presence of an elastic foundation makes the structure stiffer and more stable. 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

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