Free Vibration Analysis of Functionally Graded Porous Cylindrical Panels Reinforced with Graphene Platelets

The free vibration of functionally graded porous cylindrical shell panels reinforced with graphene platelets (GPLs) was numerically investigated. The free vibration problem was formulated using the first-order shear deformation shell theory in the framework of the 2-D natural element method (NEM). The effective material properties of the GPL-reinforced shell panel were evaluated by employing the Halpin–Tsai model and the rule of mixtures and were modified by considering the porosity distribution. The cylindrical shell surface was transformed into the 2-D planar NEM grid to avoid complex computation, and the concept of the MITC3+shell element was employed to suppress shear locking. The numerical method was validated through benchmark experiments, and the free vibration characteristics of FG-GPLRC porous cylindrical shell panels were investigated. The numerical results are presented for four GPL distribution patterns (FG-U, FG-X, FG-O, and FG-Λ) and three porosity distributions (center- and outer-biased and uniform). The effects of GPL weight, porosity amount, length–thickness and length–radius ratios, and the aspect ratio of the shell panel and boundary condition on the free vibration characteristics are discussed in detail. It is found from the numerical results that the proposed numerical method accurately predicts the natural frequencies of FG-GPLRC porous cylindrical shell panels. Moreover, the free vibration of FG-GPLRC porous cylindrical shell panels is significantly influenced by the distribution pattern as well as the amount of GPLs and the porosity.


Introduction
Graphene platelets (GPLs) have been widely used as advanced nanofiller materials to improve the physical properties of composites due to their extraordinary physical properties [1,2]. The elastic modulus of GPLs is much higher than that of polymers, so that the structural stiffness of polymer composites can be dramatically increased when only a small amount of graphene sheets are added [3]. Furthermore, it has been reported that graphene provides a higher aspect ratio, larger specific surface area, and lower production costs than carbon nanotubes (CNTs). Moreover, it has been found that the elastic modulus, tensile strength, and fracture toughness of graphene-platelet-reinforced composites (GPLRCs) are higher than those of CNT-reinforced nanocomposites [4]. However, the reinforcement with GPLs is subjected to the limitation of weight fraction owing to the high cost of FGPLs. Owing to this limitation, GPLs are generally blended into the polymer matrix according to the concept of functionally graded materials (FGMs) [5]. FGMs are characterized by the continuous and functional distributions of reinforcement materials through the thickness of composite structures [6]. The physical behavior of GPL-reinforced composites is influenced by this functional distribution pattern, and several purposeful distribution patterns of GPLs have been introduced and investigated: FG-U, FG-V, FG-X, FG-O, and FG-Λ [7,8].
More recently, the GPL-reinforced composite has been developed into a closed-cell porous reinforced foam due to rapidly developed nanotechnology [9][10][11]. As a typical example of porous materials, foams are characterized by high porosity, low density, and large specific surface. Thus, the porous GPL-reinforced composites can synergistically combine the excellent properties of porous foams and GPL nanocomposites [12]. The variation in microstructure and porosity within a porous material can also be purposefully tailored according to the concept of FGMs to enhance the designed material's performance. These so-called functionally graded porous materials have attracted great attention for the development of next-generation lightweight structures [13]. The porous composite materials reinforced with GPLs in which both internal pores and GPLs are functionally distributed through the structure's thickness are called functionally graded GPL-reinforced composite (FG-GPLRC) porous structures. In general, the functional distribution patterns are different to each other.
According to the literature survey on FG-GPLRC porous structures, Kitipornchai et al. [14] presented a micromechanics model of FG-GPLRC beams using Timoshenko beam theory and the Ritz method and parametrically investigated the free vibration and buckling behaviors. Gao et al. [15] performed a nonlinear primary resonance analysis of FG-GPLRC porous cylindrical shells under a uniformly distributed harmonic load. Akbas [16] investigated the free vibration and bending deformation of a simply supported functionally graded plate with the porosity effect using the first-order shear deformation plate theory. Barati and Zenkour [17] investigated the postbuckling behavior of geometrically imperfect porous beams reinforced with graphene platelets and resting on a nonlinear hardening foundation. Sahmani et al. [18] presented the size-dependent nonlinear bending of FG-GPLRC porous beams subjected to a uniform distributed load and an axial compressive load using the nonlocal strain gradient theory. Zhou et al. [19] performed a nonlinear buckling analysis of FG-GPLRC porous cylindrical shells under axial compressive load by considering the pre-buckling effect and in-plane constraint. Liu et al. [20] presented an analytical approach for nonlinear static responses and the stability of FG-GPLRC porous arches based on the Euler-Bernoulli hypothesis. Nguyen et al. [21] presented an efficient polygonal finite element method (PFEM) to numerically investigate the static and free vibrations of FG-GPLRC porous plates using Timoshenko's beam theory. Tao and Dai [22] investigated the postbuckling behavior of sandwich cylindrical shell panels with an FG-GPLRC porous core and two metallic face layers using a general higher-order shear deformation shell theory. Wang and Zhang [23] investigated the thermal buckling and postbuckling behaviors of FG-GPLRC porous beams by considering the temperature-dependent material properties. Wang et al. [24] investigated the forced vibration of FG-GPLRC porous cylindrical microshells subjected to time-dependent distributed impulsive loads. Zhang et al. [25] investigated the thermal buckling of FG-GPLRC porous cylindrical panels by considering the temperature-dependent material properties.
As revealed in the relevant literature survey, most of the studies on FG-GPLRC porous structures have been confined to beam and plate structures by analytical approaches based on the shear deformation theory or by numerical approaches using FEM. The studies on cylindrical structures have been poorly presented and mostly restricted to the buckling behavior. In this context, the purpose of this study was to thoroughly examine the free vibration of FG-CNTRC porous cylindrical shell panels by developing a 2-D effective and locking-free meshfree-based numerical method. The numerical method was developed in the framework of the 2-D planar natural element method (NEM), a previously introduced meshfree method characterized by high smooth Laplace interpolation functions [26,27], in which a geometry transformation and the MITC (mixed-interpolated tensorial components) approach [28] are integrated in order to relax the painstaking manipulation on the shell surface and the troublesome shear locking which occurs in the bending-dominated thin structures [29]. The displacement field of the cylindrical shell is expressed by the first-order shear deformation shell theory, and the physical shell surface and the 2-D rectangular planar NEM grid are correlated through the geometric transformation.
The numerical method was verified through experiments of benchmark examples, and the free vibration responses of FG-GPLRC porous cylindrical shell panels were investigated with respect to the GPL-and porosity-associated parameters, shell geometry dimensions, and boundary conditions. Following the introduction, the FG-GPLRC porous cylindrical shell panel and its displacement, strain and stress fields, and effective material properties are described in Section 2. The natural element approximation of the FG-GPLRC porous cylindrical panel is fully explained in Section 3. The benchmark and parametric numerical experiments are presented in Section 4 with a discussion on the numerical results. The final conclusion is made in Section 5. Figure 1a shows a cylindrical shell panel reinforced with graphene platelets (GPLs), where a coordinate system (x, y, z) is introduced on the mid-surface of the panel using the relation of x = Rθ. The geometric dimensions of the cylindrical panel are represented by radius R, length , sub-tended angle θ 0 , and uniform thickness h. Graphene platelets in this study were distributed according to a specific functionally graded pattern through the thickness. Figure 1b depicts the four patterns adopted for this study, where the GPLs are uniformly dispersed in FG-U, whereas they are rich at the mid-surface in FG-O, at the top surface in FG-X, and at the bottom surface in FG-Λ.

Modeling of FG-GPLRC Porous Cylindrical Shell Panel
according to the linear rule of mixture.  Figure 2 represents three different porosity distributions: center-biased, outer-biased, and uniform, which are expressed as Letting V GPL (z) and V m (z) be the volume fractions of GPLs and the underlying matrix material, then both satisfy the physical relation given by In which the GPL volume fraction V GPL (z) is expressed by different functions in terms of the thickness coordinate z and the total GPL volume V * GPL such that for different GPL distribution patterns, where with g GPL being the GPL mass fraction, and ρ GPL and ρ m being the densities of the GPLs and matrix material. The GPLs are assumed to be uniformly dispersed within the matrix and act as an effective rectangular solid fiber with length l GPL , width w GPL , and thickness t GPL , and the graphene-reinforced composites are modeled as an isotropic material with the effective material properties. The effective elastic modulus E C of GPLRC is estimated by the Halpin-Tsai micromechanical model [30], which gives Here, E GPL and E m denote the elastic moduli of the GPLs and matrix material, and ξ L and ξ T are the geometric parameters given by Meanwhile, the effective mass density ρ C and Poisson's ratio ν C of the GPLRC are estimated as according to the linear rule of mixture. Figure 2 represents three different porosity distributions: center-biased, outer-biased, and uniform, which are expressed as with e 0 (0 ≤ e 0 ≤ 1) being the porosity coefficient. The porosity influences the elastic modulus E C , shear modulus G C = E C /2(1 + ν C ), and mass density ρ C of the GPLRC. Letting ℘(z) be the effective material properties (i.e., E, G, and ρ) when the porosity is considered, then it is calculated by from the material properties ℘ C (z) of the GPLRC. Here, for the effective mass density, the porosity coefficient e 0 should be modified using the relationship given by between the relative mass density and relative elastic property [31]. By letting e m be the porosity parameter for the mass density, it is modified as follows for Porosity 1, for example.
for Porosity 1, for example. By using the first-order shear deformation shell theory, the displacement field   By using the first-order shear deformation shell theory, the displacement field with d = u 0 , v 0 , w 0 , ϑ x , ϑ y T being the displacement components at the mid-surface of the shell panel. The strain-displacement relations are expressed as with r = R + z ≈ R. Here, ε xx , ε yy , ε xy and γ yz , γ zx are in-plane strains and transverse shear strains, and H and H s are the (3 × 5) and (2 × 5) partial derivative matrices defined by with H x = ∂/∂x and H y = ∂/∂y. Then, the constitutive relations are expressed as using the in-plane stresses σ xx , σ yy , σ xy and the transverse shear stresses τ yz , τ zx .

Analysis of Free Vibration Using 2-D NEM
For the free vibration analysis of the FG-GPLRC cylindrical shell panel by 2-D NEM, the mid-surface of the shell panel is discretized into a finite number of nodes and Delaunay triangles, as depicted in Figure 3. Then, the approximate displacement u h (x, y, z) is expressed as using Laplace interpolation functions ψ J (x, y) [26,32] and the nodal vector where the subscript J indicates the J-th node within the NEM grid C composed of N nodes.

Analysis of Free Vibration Using 2-D NEM
For the free vibration analysis of the FG-GPLRC cylindrical shell panel by 2-D NEM, the mid-surface ϖ of the shell panel is discretized into a finite number of nodes and Delaunay triangles, as depicted in Figure 3. Then, the approximate displacement [26,32] and the nodal vector    The definition of the Laplace interpolation function and its manipulation on the cylindrical surface are complex and painstaking. To relax this difficulty, a geometry transformation T C is introduced to correlate the physical NEM grid C = [0, Rθ 0 ] × [0, ] on the cylindrical surface and the computational NEM grid R = [0, θ 0 ] × [0, ] on the rectangular plane with coordinates ζ 1 and ζ 2 : Then, Laplace interpolation functions ψ J (x, s) are mapped to ϕ J (ζ 1 , ζ 2 ), and the relations of x = R · ζ 1 and y = ζ 2 lead to the inverse Jacobi matrix J −1 given by The partial derivatives H x and H y in Equations (18) and (19) on the cylindrical surface are switched to on the rectangular plane according to the chain rule. Introducing Equation (25) into Equations (18) and (19) results inĤ andĤ s in which H x and H y are replaced with H 1 and H 2 : Then, the NEM approximations of the in-plane strains ε in Equation (16) and the transverse shear strains γ in Equation (17) lead to The direct approximation (20) of transverse shear strain γ using the standard C 0 − interpolation functions can frequently suffer from a big approximation error caused by the shear locking [29,33,34]. To ensure numerical accuracy, the transverse shear strains are indirectly interpolated by adopting the three-node triangular MITC3+shell finite element represented in Figure 4 [28] with coordinates ξ and η. Each triangular element e in the physical NEM grid C shown in Figure 3 is mapped to its master triangular element . Letting N K (ξ, η) be the linear triangular FE shape functions [35], the approximate displacement field u h (x, y, z) in Equation (22)  The direct approximation (20) of transverse shear strain γ using the st interpolation functions can frequently suffer from a big approximatio caused by the shear locking [29,33,34]. To ensure numerical accuracy, the transvers strains are indirectly interpolated by adopting the three-node triangular MITC3+ nite element represented in Figure 4 [28] with coordinates ξ and η . Each tria element e ϖ in the physical NEM grid C ℑ shown in Figure 3 is mapped to its triangular element π . Letting  Then, referring to Figure 4, the element-wise transverse shear strains γˆ ar polated as: using the transverse shear strains at the tying points A, B, C, and D within the three-node triangular master element andĉ = γ yz . The analytic derivation of Equations (28) and (29) using Equation (17), the FE re-interpolation, and the chain rule between two coordinates (x, y) and (ξ, η) leads toγ e =B e d e with the (2 × 15) matricesB e in function of ξ, η, z and R and the (15 × 1) element-wise nodal vectors d e = {d e 1 , d e 2 , d e 3 }. Meanwhile, the standard Galerkin weak form for the free vibration analysis of FG-GPLRC cylindrical shell panels can be derived from the dynamic form of the energy principle [36] Here, m is a (5 × 5) symmetric matrix defined by: with the (3 × 3) identity matrix I and m 2 = diag z 2 , z 2 . Assuming the harmonic motion d = d · e jωt and plugging Equations (27)- (29) into Equation (30), together with the constitutive relations (20) and (21), one can derive the modal equation given by: to solve the natural frequencies {ω I } N I=1 and the natural modes d I N I=1 of the cylindrical panel which is discretized into M Delaunay triangles. Here, two stiffness matrices and the mass matrix are defined by: where (20), and D s given byD with the shear correction factor κ = 5/6, the longest side length L e of the Delaunay triangles, and a positive constant α(α > 0) called the shear stabilization parameter [34,37]. The value of α is chosen through the preliminary experiment, and this modification of the shear modulus matrix was proposed to stabilize the MITC3 element. Meanwhile, β is the porosity stabilization parameter which is dependent on the porosity distribution pattern.

Results and Discussion
A Fortran program was coded according to the numerical formulae presented in Section 3 and integrated into the previously developed NEM program [38] which was developed for plate-like structures. The numerical integration of stiffness and mass matrices in Equations (33)-(35) was carried out triangle by triangle using seven Gauss integration points for K σ and M and one Gauss point for K e s . Referring to Figure 3, uniform 11 × 11 NEM grids were used for the numerical experiments, and a total of 15 modes were extracted using the Lanczos transformation and Jacobi iteration methods, unless otherwise stated. Meanwhile, three types of boundary conditions, simply supported (S), clamped (c), and free, were considered, where S and C were implemented as The present method was compared with the other methods with two isotropic, one FG-GPLRC, and one FG-GPLRC porous cylindrical panels.
The first example is a clamped isotropic cylindrical panel with the geometry dimensions given by R = = 0.762 m, s = 0.1016 m, h = 0.00033 m. The elastic modulus E, Poisson's ratio ν, and density ρ are 68.948 GPa, 0.33, and 2711.5 kg/m 3 , respectively. The maximum relative differences of the numerical results given in Table 1 are 5.819% compared with the experimental data [39] and 1.394% and 1.097% compared with references [40,41], respectively. Thus, it has been verified that the present method is in good agreement with the three reference results. The second example is the simply supported isotropic cylindrical panels with four different aspect ratios s/ . The relative geometry dimensions are /R = 0.1, /h = 10, the Poisson's ratio is ν = 0.3, and the fundamental frequencies are calibrated aŝ ω 1 = ω 1 ρ(1 − ν 2 )/E. Yang and Shen [41] and Kobayashi and Leissa [42] in Table 2 employed the higher-order and the first-order shear deformation shell theories, while Chen and Chao [42] adopted the 3-D elasticity theory. The maximum relative differences are 0.882% compared with reference [40] and 1.706% and 2.069% compared with references [42] and [43], respectively. It has been confirmed again that the present method is in excellent agreement with the existing reference solutions for different aspect ratios with the maximum relative difference equal to 2.069%. The third example is the functionally non-porous cylindrical panels reinforced with graphene platelets with different ratios of a/h and R/a for two different boundary conditions, SSSS and CCCC, where the combined four capital letters indicate a set of boundary conditions specified for the four sides 1 , 2 , 3 , and 4 of the cylindrical panel, as show in Figure 3. The shell radius R and a/b are 10 m and 1.0, and the GPL weight fraction g GPL is set to 1.0%. Epoxy is taken as the polymer matrix and its material properties are E m = 3.0 GPa, v m = 0.34, and ρ m = 1200 kg/m 3 . Referring to Yasmin and Daniel [44] and Rafiee et al. [4], the geometry dimensions of GPLs are l GPL = 2.5 µm, t GPL = 1.5 nm, and w GPL = 1.5 µm while the material properties are E GPL = 1.01 TPa, v GPL = 0.186, and ρ GPL = 1060 kg/m 3 . The fundamental frequencies are calibrated asω 1 = ω 1 2 /h ρ m /E m , and the computed normalized fundamental frequencies are compared in Table 3 with the numerical results obtained by Van Do and Lee [5] using the isogeometric analysis (IGA) method.
It can be seen in Table 3 that the fundamental frequencies obtained by the present method are as a whole smaller than those obtained by the IGA method, except for several exceptional cases. Regarding the GPL distribution pattern, the relative differences between the present method and the IGA method are shown to be relatively higher at FG-O. Meanwhile, the dependence of the relative difference inω 1 on the /h, R/ and the boundary condition is not apparent. The maximum relative difference inω 1 between the two methods is found to be 4.213% at the simply supported FG-O with /h = 50 and R/ = 50. Thus, it has been verified that the present method accurately predicts the fundamental frequencies of FG-GPLRC cylindrical panels for various GPL distribution patterns, geometry dimensions, and boundary conditions. Table 3. Comparison of the normalized fundamental frequenciesω 1 of functionally graded porous cylindrical panels reinforced with graphene platelets (R = 10 m, /b = 1.0, g GPL = 1.0%). Next, the normalized fundamental frequencyω 1 of non-porous FG-GPLRC cylindrical panels was parametrically investigated. Figure 5a,b represent the variations inω 1 with respect to the GPL mass fraction g GPL for different GPL distribution patterns and boundary conditions. It is observed thatω 1 uniformly increases with the increasing value of g GPL , regardless of the GPL distribution pattern and the boundary condition. This is because the panel stiffness increase due to the increase in the GPL amount is greater than the panel mass increase. Moreover, the normalized fundamental frequency is remarkably influenced by the GPL distribution pattern and the boundary condition. The order of the magnitude ofω 1 among the four GPL distribution patterns is FG-X > FG-U > FG-Λ> FG-O. This is because the thickness-wise GPL distribution pattern strongly affects the structural stiffness such that the structural stiffness increases as the GPL distribution becomes biased to the top and bottom surfaces of the panel. Regarding the effect of the boundary condition, the order of the magnitude ofω 1 is CCCC > CSCS > CFCF > SSSS, which is consistent with the order of the strength of the boundary constraint. Figure 6a represents the effect of the length-thickness ratio /h on the normalized fundamental frequency, where the panel length is kept unchanged at 1.0 m. It is seen that ω 1 increases in proportion to /h, but this owes entirely to the calibration with 2 /h. It was found from the numerical data that the non-calibrated absolute fundamental frequency ω 1 decreases with the increasing value of /h because the panel stiffness remarkably decreases as the panel becomes thinner. Figure 6b shows the effect of the radius-length R/ on the normalized fundamental frequency, where the panel length is also kept unchanged at 1.0 m. It is seen that the normalized fundamental frequency increases in reverse proportion to the shell radius R, because the structural stiffness increases while the total mass decreases as the shell radius becomes smaller. Figure 7a shows the variation inω 1 with respect to the length-thickness ratio /h for different values of g GPL . It is seen thatω 1 uniformly increases in proportion to /h, and the increase in the slopes is almost the same regardless of g GPL . The explanation for whŷ ω 1 increases with the increasing value of /h is the same as for Figure 6a. Figure 7b shows the variation inω 1 with the aspect ratio of the shell panel, where the shell length is kept unchanged at 1.0 m. It is observed thatω 1 dramatically increases in proportion to the value of /h, because the decrease in shell width b dramatically increases the panel stiffness but reduces the panel weight. Similar to the length-thickness ratio, the increase in the slope is almost insensitive to the GPL mass fraction g GPL . is remarkably influenced by the GPL distribution pattern and the boundary condition. The order of the magnitude of 1 ω among the four GPL distribution patterns is FG-X > FG-U > FG-Λ> FG-O. This is because the thickness-wise GPL distribution pattern strongly affects the structural stiffness such that the structural stiffness increases as the GPL distribution becomes biased to the top and bottom surfaces of the panel. Regarding the effect of the boundary condition, the order of the magnitude of 1 ω is CCCC > CSCS > CFCF > SSSS, which is consistent with the order of the strength of the boundary constraint.  . It is seen that the normalized fundamental frequency increases in reverse proportion to the shell radius R , because the structural stiffness increases while the total mass decreases as the shell radius becomes smaller.   , because the decrease in shell width b dramatically increases the panel stiffness but reduces the panel weight. Similar to the length-thickness ratio, the increase in the slope is almost insensitive to the GPL mass fraction GPL g . Next, the free vibration of the FG-GPLRC porous cylindrical panel was investigated by considering three porosity distributions shown in Figure 2. Table 4 compares the normalized fundamental frequencies of the simply supported FG-GPLRC porous cylindrical panel with the reference solutions of Zhou et al. [45]. The GPL distribution pattern is FG-U, and the geometry dimensions are given in the table caption. The reference solutions were obtained by employing Reddy's third-order shear deformation theory. The fundamental frequencies were normalized asω 1 = ω 1 R ρ m /E m , and the porosity stabilization parameter β included in Equation (36) was determined through the preliminary experiment: β = 1 for Porosity_1, β = 1/ 1 + 10e 2 0 2 for Porosity_2, and β = 1/ 1 − e 2 0 2 for Porosity_3. When compared with the reference solutions, the present results are higher at Porosity_1 but lower at Porosity_2 and 3. The maximum relative difference equal to 1.08% occurred at g GPL = 1.0 and e 0 = 0.3 for Porosity_3. Thus, the comparison verifies that the present method accurately predicts the fundamental frequency of FG-GPLRC porous cylindrical panels. , because the decrease in shell width b dramatically increases the panel stiffness but reduces the panel weight. Similar to the length-thickness ratio, the increase in the slope is almost insensitive to the GPL mass fraction GPL g .

Method
(a) (b) Next, the free vibration of the FG-GPLRC porous cylindrical panel was investigated by considering three porosity distributions shown in Figure 2. Table 4 compares the normalized fundamental frequencies of the simply supported FG-GPLRC porous cylindrical panel with the reference solutions of Zhou et al. [45]. The GPL distribution pattern is FG-U, and the geometry dimensions are given in the table caption. The reference solu-   Figure 8a represents the variation inω 1 with respect to the porosity parameter e 0 , where theω 1 uniformly decreases in proportion to e 0 . This is because the decrease in panel stiffness owing to the increase in the porosity amount is greater than the decrease in the panel mass. Meanwhile, the decreasing trend is affected by the porosity distribution such that the decrease in the slope is highest at Porosity_2 and lowest at Porosity_1. This is because the decrease in the structural stiffness becomes larger as the porosity distribution becomes biased to the top and bottom surfaces of the panel. Figure 8b represents the effect of the GPL distribution pattern on the variation inω 1 with respect to the porosity parameter e 0 . It is found that the GPL distribution pattern affects the order of the magnitude ofω 1 , but its effect on the decrease in the slope ofω 1 with respect to the porosity parameter e 0 is not remarkable. This trend is different from the effect of the GPL distribution pattern on the variation inω 1 with respect to the GPL mass fraction shown in Figure 5a. This is because the decrease in the slope itself of the panel stiffness and the panel mass along the porosity parameter e 0 is almost insensitive to the GPL distribution pattern.    Figure 9a represents the effect of GPL mass fraction g GPL on the decrease inω 1 with respect to the porosity parameter e 0 . It is seen that the g GPL affects the magnitude ofω 1 , but its effect on the decrease in the slope ofω 1 along the porosity parameter is not significant. This trend is also seen in Figure 9b which represents the effect of the boundary condition on the decrease inω 1 with the porosity parameter e 0 . It is found that the magnitude ofω 1 is apparently influenced by the boundary condition, but the dependence of its decreasing slope with respect to the porosity parameter on the boundary condition is not shown to be remarkable. This is because the decrease in the slope itself of the panel stiffness and the panel mass with respect to e 0 is not remarkably influenced by the GPL mass fraction g GPL and the boundary condition.  ω , but its effect on the decrease in the slope of 1 ω along the porosity parameter is not significant. This trend is also seen in Figure 9b which represents the effect of the boundary condition on the decrease in 1 ω with the porosity parameter 0 e . It is found that the magnitude of 1 ω is apparently influenced by the boundary condition, but the dependence of its decreasing slope with respect to the porosity parameter on the boundary condition is not shown to be remarkable. This is because the decrease in the slope itself of the panel stiffness and the panel mass with respect to 0 e is not remarkably influenced by the GPL mass fraction GPL g and the boundary condition.   Figure 10a represents the effect of porosity distribution on the increase inω 1 with respect to the GPL mass fraction g GPL . It is found that the order of the magnitude ofω 1 is apparently affected by the porosity distribution, but the increase in the slope ofω 1 with the GPL mass fraction is not influenced by the porosity distribution. This trend is different from the effect of porosity distribution on the decrease inω 1 with respect to the porosity parameter shown in Figure 8a. Meanwhile, Figure 10b represents the variation inω 1 with respect to the GPL mass fraction g GPL for different porosity parameters. It is observed that the porosity parameter significantly affects the magnitude ofω 1 , but its effect on the increase in the slope ofω 1 with respect to g GPL is not shown to be remarkable. This is because the increase in the slope itself of the panel stiffness and the panel mass with respect to g GPL is not remarkably influenced by the GPL distribution and the porosity distribution. with the GPL mass fraction is not influenced by the porosity distribution. This trend is different from the effect of porosity distribution on the decrease in 1 ω with respect to the porosity parameter shown in Figure 8a. Meanwhile, Figure 10b represents the variation in 1 ω with respect to the GPL mass fraction GPL g for different porosity parameters. It is observed that the porosity parameter significantly affects the magnitude of 1 ω , but its effect on the increase in the slope of 1 ω with respect to GPL g is not shown to be remarkable. This is because the increase in the slope itself of the panel stiffness and the panel mass with respect to GPL g is not remarkably influenced by the GPL distribution and the porosity distribution.

Conclusions
This paper presents a free vibration analysis of FG-GPLRC porous cylindrical shell panels using a NEM-based 2-D numerical method. In the framework of 2-D planar NEM, the numerical method was developed by integrating a geometry transformation between the shell surface and the 2-D planar NEM grid and the MITC3+shell element. Benchmark and parametric experiments were carried out to validate the proposed numerical method and to thoroughly investigate the free vibration characteristics of FG-GPLRC porous cylindrical panels with respect to the associated parameters. The numerical results reveal the following main observations:

•
The numerical method accurately analyzes the free vibration of FG-GPLRC porous cylindrical shell panels, without causing shear locking, with the maximum relative difference of 4.213% even for coarse and 2-D planar NEM grids.

Conclusions
This paper presents a free vibration analysis of FG-GPLRC porous cylindrical shell panels using a NEM-based 2-D numerical method. In the framework of 2-D planar NEM, the numerical method was developed by integrating a geometry transformation between the shell surface and the 2-D planar NEM grid and the MITC3+shell element. Benchmark and parametric experiments were carried out to validate the proposed numerical method and to thoroughly investigate the free vibration characteristics of FG-GPLRC porous cylindrical panels with respect to the associated parameters. The numerical results reveal the following main observations:

•
The numerical method accurately analyzes the free vibration of FG-GPLRC porous cylindrical shell panels, without causing shear locking, with the maximum relative difference of 4.213% even for coarse and 2-D planar NEM grids.

•
The normalized natural frequencyω 1 uniformly increases in proportion to the GPL mass fraction g GPL while it uniformly decreases with the increasing value of porosity parameter e 0 , and it uniformly increases in proportion to the length-thickness ratio /h, the length-radius ratio /R, and the aspect ratio /b of the shell panel.

•
The distribution patterns of both the GPL and porosity significantly affect the variations inω 1 with respect to the values of g GPL and e 0 such that the order of the magnitude ofω 1 among the four GPL distribution patterns is FG-X > FG-U > FG-Λ> FG-O while that among the three porosity distributions is Porosity_1 > Porosity_3 > Porosity_2.

•
The increase in the slope ofω 1 with respect to the GPL mass fraction is influenced by the GPL distribution pattern, /h, /R, and /b, but it is independent of the magnitude and distribution of the porosity. Meanwhile, the decrease in the slope ofω 1 with respect to the porosity parameter is influenced by the porosity distribution, but it is independent of the mass fraction and distribution of the GPL and the boundary condition.