Parametric Instability of Functionally Graded Porous Cylindrical Panels under the Effect of Static and Time-Dependent Axial Loads

: This work presents an analytical study of the parametric instability of cylindrical panels containing functionally graded porous exposed to static and dynamic periodic axial loads under simply supported boundary conditions. Based on Hamilton’s principle, the governing equation of motion by using first-order shear deformation theory (FSDT) has been obtained. By applying the Galerkin technique, an excitation frequency expression is derived, which helps identify areas of instability of functionally graded porous cylindrical panels. Numerical simulations are used to validate the analytical results. Eventually, the impacts of the porosity coefficient, porosity distribution method, static and dynamic periodic axial loads, panel angle, circumferential wave number, and cylindrical panel characteristics on the region of instability are displayed in the section of results and discussions. The findings show that when the porosity is further from the surface, the more stable the structure is. Furthermore, a small angle of the cylindrical panels gives a better dynamic response than a large angle. In addition, increased static and dynamic loads lead to an expansion of areas of instability.


Introduction
The presence of porosity in a structure makes its strength-to-weight ratio very good as compared to the strength of homogeneous isotropic materials at the same weight, as well as being permeable and a thermal insulator [1]. Consequently, cylindrical panels are key structural elements present in a wide range of engineering fields, including aircraft, petrochemicals, nuclear reactors, the maritime industry, and mechanical engineering [2]. The stable performance of functionally graded porous (FGP) cylindrical panels is affected by a variety of mechanical loads (i.e., static axial load and time-dependent load). Many studies were performed on the free vibration of cylindrical shells made of homogeneous isotropic and functionally graded porous by existing authors [3][4][5][6][7], and for rotating cylindrical shells made of homogeneous isotropic, multilayered, functionally graded materials and CNTreinforced functionally graded materials, graded graphene, functionally graded porous materials, and stiffened functionally graded porous [8][9][10][11][12][13][14][15]. Authors have investigated linear and nonlinear free vibration of cylindrical panels made of homogeneous isotropic and functional graded porous materials, and analyses based on different shear deformation theories were performed in [16,17]. The effect of rotating on the cylindrical panel behavior made of an isotropic martial was examined in [18]. The study of a free vibration composite panel reinforced by CNT functionally graded panels based on the element-free KP-Ritz method under different boundary conditions was investigated in [19]. Hong et al. [20] examined the natural frequency of functionally graded cylindrical and parabolic panels based on the Fourier-Ritz approach under various boundary conditions. Dong et al. [21] studied the dynamic responses of the forward and backward traveling waves of functionally graded graphene-reinforced cylindrical shells under the effects of angular velocity and axial static load. Li et al. [22] investigated the combined impacts of subsonic airflow

Theoretical Approach
The parameters of functional graded porous cylindrical panels under study are symbolized as thickness: h; length: L; angle: θo; mean radius: R; square plan form: b; circumferential length: S, as shown in Figure 1. An orthogonal coordinate (X, θ, Z) is used to identify the structure's axes. The X and θ axes are the longitudinal and circumferential directions of the panels, respectively, and Z is the radial direction of the panels along the thickness.

The Governing Equations of the Distribution of Porosity
In this work, three types of porosity distributions were investigated, as shown in Figure 2 [17]. Type 1 is the uniform spread of porosity along the thickness direction and the properties of materials are constant (denoted by Type 1). The second type is a nonuniform symmetric distribution about the center. In this type, the midplane of the panel has the maximum porosity and the maximum stiffness is on the surfaces of panels (denoted by Type 2). The third type is a symmetric porosity about surfaces. This type has the maximum porosity on the top and bottom surfaces of the panel and the midplane of the panel has the maximum stiffness (denoted by Type 3).

The Governing Equations of the Distribution of Porosity
In this work, three types of porosity distributions were investigated, as shown in Figure 2 [17]. Type 1 is the uniform spread of porosity along the thickness direction and the properties of materials are constant (denoted by Type 1). The second type is a nonuniform symmetric distribution about the center. In this type, the midplane of the panel has the maximum porosity and the maximum stiffness is on the surfaces of panels (denoted by Type 2). The third type is a symmetric porosity about surfaces. This type has the maximum porosity on the top and bottom surfaces of the panel and the midplane of the panel has the maximum stiffness (denoted by Type 3).

Theoretical Approach
The parameters of functional graded porous cylindrical panels under study are symbolized as thickness: h; length: L; angle: θo; mean radius: R; square plan form: b; circumferential length: S, as shown in Figure 1. An orthogonal coordinate (X, θ, Z) is used to identify the structure's axes. The X and θ axes are the longitudinal and circumferential directions of the panels, respectively, and Z is the radial direction of the panels along the thickness.

The Governing Equations of the Distribution of Porosity
In this work, three types of porosity distributions were investigated, as shown in Figure 2 [17]. Type 1 is the uniform spread of porosity along the thickness direction and the properties of materials are constant (denoted by Type 1). The second type is a nonuniform symmetric distribution about the center. In this type, the midplane of the panel has the maximum porosity and the maximum stiffness is on the surfaces of panels (denoted by Type 2). The third type is a symmetric porosity about surfaces. This type has the maximum porosity on the top and bottom surfaces of the panel and the midplane of the panel has the maximum stiffness (denoted by Type 3).  The governing equations of the three types of distribution can be defined as [15]: where e 0 is the porosity coefficient, its value (0 < e 0 < 1), which can be calculated as [17] e 0 = 1 − E min /E 0 e m * is the porosity coefficient of the mass density of Type 1 distributions, which can be calculated as [17] e m * = (1 − √ 1 − e 0 λ)/λ , E 0 , ρ 0 , E min , ρ min . λ in Equation . e m is the porosity coefficient of mass density for Type 2 and Type 3 distributions, which can be calculated as e m = 1 − √ 1 − e 0 or e m = 1 − ρ min /ρ 0 . E 0 , E min , ρ 0 , ρ min are the maximum and minimum values of elastic moduli and density, respectively.

Formulation of Dynamic Stability
The displacement field equation can be defined based on Donell's shell theory using first-order shear deformation theory (FSDT) [41]: where (u, v, w) symbolize the displacement of a point in x, θ, z directions, respectively, at z = 0. ϕ x symbolizes the rotation of the normal about the θ axis, ϕ θ denotes the rotation of the normal about the x axis. The strain-displacement relation can be calculated as: To derive the governing partial differential equations of motion, Hamilton's principle is used in the following form: where δ is the variation operator, t 1 , t 2 are two arbitrary times. U denotes potential energy. T denotes kinetic energy, V denotes work done generated from static and dynamic axial forces.
The stresses in Equation (7) are given based on Hook's law, as follows: µ denotes Poisson's ratio and is assumed a constant value compared with other mechanical characteristics of the FGP cylindrical panel [17]. k denotes the shear correction factor (k = 5/6) [30].
Vibration 2022, 5 is the uniformly distributed axial periodic load per unit length, N 0 is the static axial load, N t is the time-dependent axial load, Ω refers to the excitation frequency in radian per unit time, and t denotes the time variable. Substituting Equations (7)-(10) into Equation (6), we obtain the governing partial differential equations of motion: where (N x , N θ , N xθ , Q x , Q θ ) are the force resultants and (M x , M θ , M xθ ) are the moment resultants, which are expressed as follows [40].
Now, we substitute Equation (5) into Equation (8), and then into Equation (13), and rewrite the equations of motion in the following form: The effects of the inertia forces of rotation and midplane displacement have been neglected because they are smaller than the inertia force from transverse displacement (w) [40], where L ij (i, j = 1, 2, 3, 4, 5) is the operator of the partial differential equation. The cylindrical panels are assumed to be simply supported at both ends at x = 0, L and corresponding to the boundary condition as The approximate solution that satisfies boundary conditions can be assumed to solve Equation (14) of the functionally graded porous cylindrical panel is as follows [21]: where (n, m) are the half-circumferential and axial wave numbers, respectively, β = nπ/θ 0 , α = mπ/L. After substituting Equation (15) into Equation (14) and converting the partial differential equation to an ordinary differential equation by using superposition theory, the ordinary differential equations of cylindrical panels are obtained as: where Γ ij (i, j = 1, 2, 3, 4, 5) are defined in Appendix A. So, from Equations (16) and (17) we can derive the relation between (u i , v i , ϕ xi , ϕ θi ) and (w i ), (i = 1, 2), defined as: The constants D, K, t, H are shown in Appendix B, and then, substituting Equations (19) and (20) into Equation (18), after some mathematical manipulation, we obtain the following relationship: The critical buckling periodic loads are obtained from Equation (21) by neglecting the mass term [30], so (N cr = F/α 2 ).
If we need natural frequency without periodic axial loads obtained from Equation (21), In the present study, the parametric instability is bounded by two periods, T and 2T, with T = 2π [39]. The solution to Equation (21) corresponds to 2T because the width of instability at period 2T is usually larger than T, and is assumed as [38]: where (a, b) are constants and p is the excitation frequency. Substituting Equation (23) into Equation (21) and applying Galerkin's method with periodic 2T to obtain an expression of Vibration 2022, 5 576 the excitation frequency of (FGP) a cylindrical panel with the effect of the periodic load, we can define the expression as an eigenvalue problem: Now, we rewrite Equation (24) in a new form after mathematical manipulation, as follows: N 01 is static load factor = N 0 /N cr , N t1 is dynamic load factor = N t /N cr , the negative signal gives Ω 1 , and the positive signal gives Ω 2 . The dimensionless excitation frequency is defined as follows:

Validation
The instability analysis has been validated, so three numerical examples are examined for a cylindrical panel, a cylindrical shell made of homogeneous isotropic material, an FG porous material, and sandwich shells containing an FG core material, respectively. Example 1. The validation of the natural frequency of the isotropic cylindrical panel is shown in Table 1 and against results taken from Soldatos and Hadjigeorgiou [16].  (22) is used to compare with Soldatos and Hadjigeorgiou's smallest first four frequency parameters [16]. The reliability of the results is good between the present result and the result of [16], as shown in Table 1.  [39] for sandwich cylindrical shells containing an FG core. The top layer is made of a metal-rich material, i.e., nickel (Ni), the bottom layer is made of a ceramic-rich material, i.e., silicon nitride (Si3N4), and the core of the shell is functionally graded through thickness from Si3N4 to Ni and the equations of elastic modulus and density are:  Table 2 for different wave numbers (m,n) = (1,1) and (1,3), and for a different volume friction index, d = 0.5, 2. The results have been obtained based on the first-order shear deformation theory. The results show good agreement between the present study and the results of SofIiyev and Kuruoglu [39]. Example 3. We compare the results of Yuewu and Dafang [7] for functional graded porous materials for the cylindrical shells with different values of the porosity coefficient and circumferential wave number based on Type 2 distributions, as shown in Table 3. The numerical results were performed based on Equation (22)

Results of the Present Study
In the present study, the dynamic stability of FG porous cylindrical panels is investigated under the impact of the porosity coefficient, porosity distribution, wave numbers, geometrical dimensions, static load factor, and dynamic load factor. The dimensionless excitation frequencies are obtained by applying Equation (26). The mechanical properties are [17] Figure 3 shows the effect of porosity distribution versus dynamic load on the instability of the FGP cylindrical panel, e 0 = 0.4, static load factor 0.5, and wave number (m,n) is (1,2). The results show that Type 2 distributions give the largest values of dimensionless excitation frequencies compared to other types. The results appear to show that Type 2 distributions provide the highest values of dimensionless excitation frequencies compared to other types, which is due to the surfaces' resistance to loads and deformation due to the lack of porosity at the surfaces, resulting in higher density and stiffness at the surface panel than other types. On the other hand, to enhance the stability of the structure, the porosity of the surface must be prevented. The closer the porosity is to the center, the stronger the structure will be than if the porosity was closer to the surface.

The Effect of the Panel's Angle
The effect of the angle of the panels on the instability of the structure is shown in Figure 4, with the values of N 01 = 0.2, e 0 = 0.4, and m = n = 1. The results reveal that the width of the instability region of the small angle, such as θ 0 = 30 • , increases greater than the large angle of the panel, and the origin point of instability is greater than the large angle. The reason is that when the angle increases, the flexural rigidity of the panel drops, and the decreasing amplitude of the natural frequency becomes larger with the rise of the panel angle.

The Effect of the Porosity Coefficient
The examined coefficient of porosity for Type 2 distributions against the dynamic load factor on the dimensionless excitation frequencies is presented in Figure 5. The values of the static load and wave number are N 01 = 0.2, m = n = 1, respectively. The results show the instability of the structure decreases and is wider as the porosity coefficient increases. The reason for that is the decrease in both the stiffness and mass density of the panel.    The examined coefficient of porosity for Type 2 distributions against the dynamic load factor on the dimensionless excitation frequencies is presented in Figure 5.    2). It is detected that the value of the origin point of the boundary instability region shifts forward to a high value when the static load factor decreases, and the width of the instability region is wider when the static load factor and dynamic load factor increase. The existence of dynamic load leads to the instability region appearing and the boundary of instability beginning to shift to the left. That means the structure's resistance and response dynamics decrease when the static load increases. The structure approaches failure whenever the static and dynamic load factors equal the critical value, and this means that the structure does not weaken only when the porosity increases but also weakens when the static load factor increases.   2). It is detected that the value of the origin point of the boundary instability region shifts forward to a high value when the static load factor decreases, and the width of the instability region is wider when the static load factor and dynamic load factor increase. The existence of dynamic load leads to the instability region appearing and the boundary of instability beginning to shift to the left. That means the structure's resistance and response dynamics decrease when the static load increases. The structure approaches failure whenever the static and dynamic load factors equal the critical value, and this means that the structure does not weaken only when the porosity increases but also weakens when the static load factor increases.

The Effect of the h/R Ratio
In Figure 7, the influence of the thickness-to-radius ratio on the behavior of the instability of a structure has been presented. The cylindrical panel parameters are 01 N = 0.1, and the various values of the radius-to-thickness ratio are 0.01, 0.1, and 0.2. Figure 7 shows that the dimensionless excitation frequencies increase by about 28.2% when the thickness increases from 0.01 to 0.2, and the instability regions are significantly wider when the thickness-to-radius ratio increases. To explain this behavior, we can say that the stiffness of the structure increases with the increase of the thickness.

The Effect of the h/R Ratio
In Figure 7, the influence of the thickness-to-radius ratio on the behavior of the instability of a structure has been presented. The cylindrical panel parameters are N 01 = 0.1, and the various values of the radius-to-thickness ratio are 0.01, 0.1, and 0.2. Figure 7 shows that the dimensionless excitation frequencies increase by about 28.2% when the thickness increases from 0.01 to 0.2, and the instability regions are significantly wider when the thickness-to-radius ratio increases. To explain this behavior, we can say that the stiffness of the structure increases with the increase of the thickness.   The results show that when the size of the length increases, the origin point of dimensionless excitation frequency decreases, the influence of axial dynamic load on the structure decreases significantly, and the width of the instability region decreases. It is observed that the excitation frequency is more sensitive to a change in length. When the panel is shorter, the onset of the instability region shifts to the right. 3.2.7. The Effect of the Circumferential Wave Number Table 4 shows the effect of the circumferential wave number with an angle on the instability of the cylindrical panel. The excitation frequency increases as the circumferential wave number (n) increases for the smallest angle of the panel, (i.e., 30°), but for the largest angle of the panel, such as 0 100 ,180 , the excitation frequency decreases and then increases when the circumferential wave number increases. It is observable that the lower point of the instability region increases with the increase of porosity after a specific value of porosity coefficient, in fact, the reduction in stiffness and weight remains with the porosity increase, but at this value, this phenomenon happens because the average reduction in stiffness is significantly smaller than the cross-section inertia.  3.2.7. The Effect of the Circumferential Wave Number Table 4 shows the effect of the circumferential wave number with an angle on the instability of the cylindrical panel. The excitation frequency increases as the circumferential wave number (n) increases for the smallest angle of the panel, (i.e., 30 • ), but for the largest angle of the panel, such as θ 0 = 100 • , 180 • , the excitation frequency decreases and then increases when the circumferential wave number increases. It is observable that the lower point of the instability region increases with the increase of porosity after a specific value of porosity coefficient, in fact, the reduction in stiffness and weight remains with the porosity increase, but at this value, this phenomenon happens because the average reduction in stiffness is significantly smaller than the cross-section inertia.

Conclusions
This study analyzed the parametric instability of simply supported functional graded porous cylindrical panels subjected to combined static and time-dependent periodic axial loads. The problem has been solved based on first-order shear deformation theory (FSDT) with Hamilton's principle. The general equations are converted to ordinary differential equations by using superposition theory. Herein, the influences of different parameters on the instability of cylindrical panels have been examined, and the results are fully discussed. The theoretical approach was validated compared to other studies.

•
According to the results shown, it may be concluded that the symmetric porosity distribution around the midplane (i.e., Type 2) has more stiffness and results in the excitation frequency shifting forward to a high value, so it is preferred over the rest of the types. • A small angle (i.e., θ 0 = 30 • ) of the panel results in a large excitation frequency, but it is more influenced by the dynamic load. This behavior is fully reversed when the angle of the panel is large.

•
By increasing the circumferential wave number, the excitation frequencies decrease and then increase when the angle of the panel is large (i.e., θ 0 = 100 • , 180 • ), but when the angle of the panel is small (i.e., θ 0 = 30 • ) the excitation frequencies increase by increasing the circumferential wave number.

•
The excitation frequencies increase when the static load factor decreases. Additionally, the width of instability increases as the static load factor and dynamic load factor increase.

•
The structure is less stable when it has porosity, and the weakness of the structure increases with an increase in the porosity.

•
For design purposes, care should be taken to decide the values of the static load factor and the porosity coefficient because a wrong selection leads to instability and then an early failure in the structure.

•
The dimensionless excitation frequencies and width of the instability region have a large value when the thickness is large and when the length is small.