A Size-Dependent Finite Element Method for the 3D Free Vibration Analysis of Functionally Graded Graphene Platelets-Reinforced Composite Cylindrical Microshells Based on the Consistent Couple Stress Theory

Within a framework of the consistent couple stress theory (CCST), a size-dependent finite element method (FEM) is developed. The three-dimensional (3D) free vibration characteristics of simply-supported, functionally graded (FG) graphene platelets (GPLs)-reinforced composite (GPLRC) cylindrical microshells are analyzed. In the formulation, the microshells are artificially divided into numerous finite microlayers. Fourier functions and Hermitian C2 polynomials are used to interpolate the in-surface and out-of-surface variations in the displacement components induced in each microlayer. As a result, the second-order derivative continuity conditions for the displacement components at each nodal surface are satisfied. Five distribution patterns of GPLs varying in the thickness direction are considered, including uniform distribution (UD) and FG A-type, O-type, V-type, and X-type distributions. The accuracy and convergence of the CCST-based FEM are validated by comparing the solutions it produces with the exact and approximate 3D solutions for FG cylindrical macroshells reported in the literature, for which the material length scale parameter is set at zero. Numerical results show that by increasing the weight fraction of GPLs by 1%, the natural frequency of FG-GPLRC cylindrical microshells can be increased to more than twice that of the homogeneous cylindrical microshells. In addition, the effects of the material length scale parameter, the GPL distribution patterns, and the length–to–thickness ratio of GPLs on natural frequencies of the FG-GPLRC cylindrical microshells are significant.


Introduction
In recent decades, due to the successive discovery of nanoscale materials, namely, carbon nanotubes (CNTs) and graphene sheets (GSs), the development of material technology has undergone rapid changes [1,2]. Due to the superior physical, chemical, thermal, and electrical material properties of CNTs and GSs, the application of CNTs and GSs has gradually become popular in cutting-edge technology industries, such as micro-electromechanical systems and nano-electro-mechanical systems, sensors and actuators, aerospace, submarine, automobile, and biological technologies [3][4][5]. On the one hand, because of their extreme properties of high stiffness-to-weight ratio and high strength-to-weight ratio, CNTs and GPLs have been used as reinforcement materials and have been embedded in specific matrices to form nanocomposite structures, with some predesigned functions of CNT and GPL distribution patterns varying in the thickness direction [6][7][8]. On the other hand, various systems and components used in the above-mentioned cutting-edge technologies have been gradually miniaturized. Therefore, the relevant mechanical analyses of CNT-reinforced composite (CNTRC) and GPLRC microstructures have attracted considerable attention.
It has been well-known that due to the size-dependent effect, the classical continuum mechanics is no longer applicable to various mechanical analyses of microstructures [9][10][11]. In order to capture the size-dependent effect, some higher-order nonclassical continuum mechanics-based theories have been developed in the past several decades, including the nonlocal elasticity theory [12][13][14], the micropolar elasticity [15], doublet mechanics [16], the strain gradient theory [17], and the couple stress theory (CST) [18]. This article reviews the development of the CST, the CST-based classical shell theories, and the CST-based shear deformation shell theories.
Due to Eringen's criticism [19] of the uncertainty problem of the CST, Yang et al. [20] and Hadjesfandiari and Dargush [21,22] overcame these issues with the CST. They successively developed the modified CST (MCST) and the CCST based on a symmetric and a skew-symmetric couple stress tensor assumption, respectively. Afterward, the material length scale coefficients required for analyzing elastic isotropic solids were reduced from two to one, making the MCST and the CCST practically feasible. The MCST and CCST have been successfully applied to analyze the torsional behavior of a thin cylinder and the pure bending behavior of a flat plate with an infinite width [20,21].
Based on the MCST, some two-dimensional (2D) microstructure-dependent shell theories have been developed for analyzing various mechanical behaviors of FG cylindrical and conical microshells. Tang et al. [23] developed a microstructure-dependent Kirchhoff thin plate theory to study the microplate's static bending, buckling, and free vibration behaviors. Incorporating the kinematic model of Love's classical shell theory (LCST) into the MCST, Mehralian and Beni [24] developed a size-dependent LCST for the buckling analysis of single-walled CNTs (SWCNTs). The SWCNTs of interest were placed on a Winkler-Pasternak foundation and were subjected to axial compression, electric voltages, and constant temperature changes. Afterward, based on the MCST-based LCST, Zeighampour and Beni [25] examined the free vibration characteristics of a single-walled carbon conical nanoshell. The results showed that the material length scale parameter causes the nanoshell to become stiffer, increasing the nanoshell's lowest natural frequencies. Liu and Wang [26] presented size-dependent free vibration and buckling analyses of 3D graphene foam cylindrical microshells. In conjunction with the MCST and the kinematic model of the LCST, Soleimani et al. [27] developed the weak formulation of a two-node size-dependent axisymmetric shell element. Then, they applied it to carry out a buckling analysis of circular microplates, conical microshells, and cylindrical microshells. Combining the MCST with the kinematic model of the first-order shear deformation theory (FSDT), Gholami et al. [28] analyzed the axial buckling and dynamic instability behaviors of FG cylindrical microshells. The corresponding governing equations were derived using Hamilton's principle and were written as Mathieu-Hill equations. Bolotin's method was used to determine the instability regions. Afterward, with this MCST-based FSDT, Salehipour et al. [29] studied the size-dependent free vibration and static bending behaviors of FG porous microshells and nanoshells under clamped and simply-supported boundary conditions; Zeighampour and Beni [30] examined the free vibration characteristics of SWCNTs. Esfahani et al. [31] analyzed the vibration characteristics of a conical sandwich shell with a saturated FG porous core. Effective material properties of the porous shell were estimated using Biot's theory, and the deformations of the porous shell were described using the FSDT. Shameli et al. [32] studied nanorods' free torsional vibration characteristics with noncircular cross-sections based on the second-order strain gradient theory. Incorporating the displacement model of the higher-order shear deformation theory into the MCST, Lyu et al. [33] investigated the thermo-electro-mechanical free vibration and buckling behaviors of an FG piezoelectric porous cylindrical microshell. Based on the modified strain gradient theory, Thai et al. [34] presented an isogeometric analysis (IGA) for studying the postbuckling behavior of an FG microplate subjected to mechanical and thermal loads. In their formulation, the kinematic model of Reddy's refined shear deformation theory (RSDT) and the von Kármán geometrical nonlinearity (VKGN) were employed to describe the deformation of the microplate. The nonuniform rational B-splines basis functions were used to generate the corresponding Materials 2023, 16,2363 3 of 26 C 2 interpolation (shape) functions for each primary variable. The symbol C 2 represents the continuity conditions of the order of the derivatives for each primary variable up to and including the second order. Afterward, Thai et al. [35] employed this IGA approach to study the nonlinear static bending and the nonlinear transient vibration behaviors of FG microplates. Nguyen et al. [36] used this IGA approach to analyze the vibration characteristics of cracked FG microplates.
The kinematic models of the LCST, FSDT, RSDT, and sinusoidal shear deformation theory (SSDT) have also been incorporated into the MCST for analyzing various mechanical behaviors of FG-GPLRC microshells. Wang et al. [37] used the MCST-based LCST to analyze the size-dependent vibration of FG-GPLRC cylindrical microshells. In the analysis, they considered four GPL distribution patterns. They estimated the effective Young's modulus and Poisson's ratio using the modified Halpin-Tsai model [38] and the rule of mixtures [38], respectively. Baghbadorani and Kiani [39] presented a size-dependent analysis for the free vibration characteristics of FG-GPLRC cylindrical microshells using the MCST-based FSDT. With the MCST-based FSDT, Safarpour et al. [40] analyzed the thermal buckling and the free and forced vibration behavior of FG multilayered GPLRC nanostructures resting on a Winkler-Pasternak foundation; Salehi et al. [41] carried out the nonlinear vibration analysis of an imperfect FG-GPLRC porous cylindrical nanoshell. Wang et al. [42] conducted size-dependent research for the dynamic instability behavior of FG-GPLRC cylindrical micropanels using the MCST-based SSDT. Based on the MCST, Moayedi et al. [43] reformulated the classical continuum mechanics-based RSDT to study the size-dependent thermal buckling responses of FG-GPLRC micropanels.
After a close literature survey, we found that almost all relevant mechanical analyses of FG cylindrical microshells involved using size-dependent 2D theories, not the sizedependent 3D approach for elastic solids. In order to completely capture the 3D effects on the structural behavior of an FG microplate, including thickness stretching effects, shear deformation effects, 3D couple-stress tensor effects, and zig-zag deformation effects (for laminated microscale plates), Wu and Hsu [44] and Wu and Lu [45] developed the Lagrangian C 0 and the Hermitian C 1 FEMs based on the CCST to analyze FG elastic microplates. By specifying the material length scale parameter as zero, their results agree well with the exact and approximate 3D solutions of FG macroplates reported in the literature. The results also showed that the convergence rate of the Hermitian C 1 FEM was more rapid than that of the Lagrangian C 0 FEM. In order to further speed up the convergence rate of the Hermitian C 1 FEMs and extend their application from the analysis of microplates to the study of cylindrical shells, we aim to develop a Hermitian C 2 FEM based on the CCST to investigate the 3D free vibration characteristics of a simply-supported FG-GPLRC cylindrical microshell. A 3D weak formulation to address the current issue is first derived using Hamilton's principle by selecting the displacement components as the primary variables. Then, Fourier functions and Hermitian C 2 polynomials are used to interpolate the variations of the primary variables in the nodal surface and the thickness direction, respectively. As a result, the second-order derivative continuity conditions for the primary variables at each nodal surface are satisfied. In the numerical examples, five GPL distribution patterns varying in the thickness direction are considered, with a constant weight fraction of GPLs, including UD and FG A-type, O-type, V-type, and Xtype distributions. A parametric study related to some key effects on the lowest natural frequencies of the FG-GPLRC cylindrical microshell is conducted, including the impacts of the material length scale parameter, the GPL distribution patterns, the GPL weight fractions, the length-to-mid-surface radius and mid-surface radius-to-thickness ratios of the microshell, and the length-to-thickness ratio of GPLs.

Effective Material Properties
The schematic diagram of a simply-supported FG cylindrical microshell of interest in this study is shown in Figure 1, for which the dimensions of the microshell fall within 0 ≤ x ≤ L, 0 ≤ θ ≤ 2π, and a ≤ r ≤ b. The variables x, θ, and r are the cylindrical shell coordinates. The variables a and b denote the inner and outer radii of the microshell, respectively. The variables h, L, and R represent the thickness, the length, and the mid-surface radius of the microshell, respectively. Finally, the variable ζ denotes the global thickness coordinate of the microshell measured from the mid-surface. Some relationships between the above geometric variables are given as r = R + ζ, a = R − (h/2), and b = R + (h/2).

Effective Material Properties
The schematic diagram of a simply-supported FG cylindrical microshell of interest in this study is shown in Figure 1, for which the dimensions of the microshell fall within 0 x L ≤ ≤ , 0 2 θ π ≤ ≤ , and a r b ≤ ≤ . The variables x, θ , and r are the cylindrical shell coordinates. The variables a and b denote the inner and outer radii of the microshell, respectively. The variables h, L, and R represent the thickness, the length, and the mid-surface radius of the microshell, respectively. Finally, the variable ζ denotes the global thickness coordinate of the microshell measured from the mid-surface. Some relationships between the above geometric variables are given as r R ζ = + , a = R − (h/2), and b = R + (h/2). The microshell of interest is made of FG-GPLRC material, which is formed by mixing a specific matrix material and a reinforcement material GPL according to a particular function of GPL distribution pattern, varying in the thickness direction when the weight fraction ( * GPL W ) of GPLs is fixed. Five GPL distribution patterns varying in the thickness direction are considered in the following numerical examples: UD and FG A-type, O-type, Vtype, and X-type distributions. In addition, the Halpin-Tsai model [36] is used to estimate the effective Young's modulus, and the rule of mixtures [36] is used to estimate the effective Poisson's ratio and the effective mass density. Finally, the following equations are used to describe these material properties.
According to the Halpin-Tsai model, the effective Young's modulus of the FG-GPLRC cylindrical microshell eff E can be approximated with where L E and T E denote the longitudinal and transverse moduli, respectively, which are expressed following Yang et al. [46] as: The microshell of interest is made of FG-GPLRC material, which is formed by mixing a specific matrix material and a reinforcement material GPL according to a particular function of GPL distribution pattern, varying in the thickness direction when the weight fraction (W * GPL ) of GPLs is fixed. Five GPL distribution patterns varying in the thickness direction are considered in the following numerical examples: UD and FG A-type, O-type, V-type, and X-type distributions. In addition, the Halpin-Tsai model [36] is used to estimate the effective Young's modulus, and the rule of mixtures [36] is used to estimate the effective Poisson's ratio and the effective mass density. Finally, the following equations are used to describe these material properties.
According to the Halpin-Tsai model, the effective Young's modulus of the FG-GPLRC cylindrical microshell E e f f can be approximated with where E L and E T denote the longitudinal and transverse moduli, respectively, which are expressed following Yang et al. [46] as: where the variables E GPL and E matrix denote Young's modulus of GPLs and matrix material, respectively; the variable V GPL is the volume fraction of GPLs; the variables ξ L and ξ T are the parameters characterizing the geometrical dimensions of GPLs; the variables η L and η T are the parameters describing the geometrical dimensions of GPLs and Young's modulus ratio between GPLs and matrix material; the variables (L x ) GPL , L y GPL , and h GPL denote the length, width, and thickness of GPLs, respectively.
According to the rule of mixtures, the effective Poisson's ratio υ e f f and the effective mass density ρ e f f of the GPLRC material are given as where subscripts GPL, matrix, and eff denote GPLs, matrix material, and GPLRC material, respectively. The variables V GPL and V matrix are the volume fractions of GPLs and matrix material, respectively. In the numerical example, with a specific value of W * GPL , the weight fractions of five relevant distribution patterns of GPLs varying in the thickness direction of the microshell (W GPL (ζ)) are expressed as follows: where the weight fractions of the above five GPL distribution patterns varying in the thickness direction are shown in Figure 2. The total volume fraction of GPLs in each case in Equations (10)- (14) can be obtained using the relationship between V GPL and W GPL , which is

The 3D CCST for an Elastic Solid
Within a framework of the CCST for an elastic solid developed by Hadjesfandiari and Dargush [21,22], the force-stress tensor ( ij σ ) at a material point induced in a deformed elastic solid is assumed to be asymmetric, and the couple-stress tensor ( ij μ ) is skew-symmetric. Hadjesfandiari

The 3D CCST for an Elastic Solid
Within a framework of the CCST for an elastic solid developed by Hadjesfandiari and Dargush [21,22], the force-stress tensor (σ ij ) at a material point induced in a deformed elastic solid is assumed to be asymmetric, and the couple-stress tensor (µ ij ) is skew-symmetric. Hadjesfandiari and Dargush thus decomposed the force-stress tensor into the symmetric part (σ (ij) ) and the skew-symmetric part (σ [ij] ), with the former and the latter being distinguished using parentheses and brackets surrounding a pair of indices, respectively. In the formulation presented by Hadjesfandiari and Dargush, the skew-symmetric part of the force-stress tensor (σ [ij] ) was expressed in terms of the couple-stress tensor (µ ij ) as follows: where subscripts i, j, and k permute in a natural order; and µ k = µ ji = −µ ij . The strain energy density of an elastic microsolid is a function of the strain tensor (ε ij ) and the skew-symmetric part of the curvature tensor (κ ij ). The strain tensor is a symmetric tensor conjugated with the symmetric part of the force-stress tensor (σ (ij) ), and the skewsymmetric part of the curvature tensor (κ ij ) is conjugated with the couple-stress tensor (µ ij ), which is a skew-symmetric tensor. Therefore, the strain energy in a microscale elastic solid occupying a volume Ω can be written as follows: where σ (ij) = C ijkl ε kl , and C ijkl is the elastic coefficient; ε kl = (u k,l + u l,k )/2; , and θ kj denotes the rotation tensor, which is a skew-symmetric tensor; and µ ij = −8Gl 2 κ ij for an elastic solid, where G is the shear modulus and l denotes the material length scale parameter characterizing the size-dependent effects. The detailed expressions for the relationships between the above tensors and the displacement comments are given in the following section.

Kinematics Assumptions
Based on the CCST for an elastic solid, we develop a weak formulation for analyzing the free vibration characteristics of a simply-supported FG-GPLRC cylindrical microshell. In the formulation, the microshell is artificially divided into n l cylindrical microlayers, for which the domain of each microlayer falls within 0 ≤ x ≤ L x , 0 ≤ θ ≤ 2π, (−h m /2) ≤ z m ≤ (h m /2), and m = 1 − n l . The variable z m represents the local thickness coordinate measured from the midsurface of each microlayer. The elastic displacement components of each microlayer of the FG-GPLRC cylindrical microshell for each Hermitian C 2 element are given as follows: where t denotes the time variable; n d represents the number of nodal surfaces in each microlayer; the superscript m indicates the mth-microlayer; the variables u  (i = 1, 2, . . . , 3n d ) are the shape (or interpolation) functions that consist of Hermitian C 2 polynomial functions and satisfy the continuity conditions for the second-order derivatives of the primary variables at each nodal surface; and ψ For each microlayer, the linear constitutive equations, which are valid for the orthotropic elastic materials, following Hadjesfandiari and Dargush [21,22], are given as follows: where the variables σ represents the shear modulus related to the i-j surface for the mth-microlayer The variable l (m) i is the material length scale parameter related to the kj-surface for the mth-microlayer. When the isotropic material is considered, the above-mentioned b (m) for which l denotes the material length scale parameter. It should be noticed that the value of the material length scale parameter (l) in the CCST is one-half of the value of the material length scale parameter (l) in the MCST. The above relation is because the relationship between the couple stress tensor and the skew-symmetric part of the curvature tensor is µ = −8Gl 2 κ in the CCST, while the relationship between the couple stress tensor and the symmetric part of the curvature tensor is m = 2Gl 2 χ. For example, the value ofl of the epoxy material isl = 17.6 × 10 −6 m, and that of l is l = 8.8 × 10 −6 m.
The strain-displacement relationships for each microlayer are given as follows: where m = 1, 2, . . . , n l ; the commas denote partial differentiation with respect to the suffix variables; and Dψ The skewsymmetric parts of the curvatures-to-the-displacements relations for each microlayer are given by:

Hamilton's Principle
The Euler-Lagrange equations to address the issue that we are discussing now are derived using Hamilton's principle, and its corresponding energy functional is expressed as follows: where T (m) and U (m) s denote the kinetic and strain energy of a typical mth-microlayer, respectively, and they are given as follows: where Ω denotes the domain of the cylindrical microshell on the x − θ surface. As mentioned above, we take the elastic displacement components as the primary variables subject to variation. Therefore, we first perform the first-order variation of the kinetic and strain energy based on the generalized kinematics assumptions given in Equations (3)-(5), and then employ the integration by parts technique, which results in the following equations: where the superscript T denotes the transposition of the matrices or vectors, and where i = 1, 2, . . . , n d , and m = 1, 2, . . . , n l .

Hermitian C 2 Finite Element Equations
In this section, we develop weak form-based Hermitian C 2 FEMs for analyzing the free vibration characteristics of simply-supported FG-GPLRC cylindrical microshells.
The boundary conditions of each microlayer at two edges of the cylindrical microshell are taken to be completely simple supports with traction loads prescribed to be zero, and they are specified as follows: The primary variables of each microlayer given in Equations (17)- (19) are further expanded as a double Fourier series in the spatial domain and as harmonic functions in the time domain such that the boundary conditions of the simply-supported edges are exactly satisfied, and they are expressed as follows: where ω denotes the natural frequency of the FG-GPLRC cylindrical microshell; m =mπ/L x , andm andn are the half-wave and full-wave numbers in x and θ directions, respectively, andm is positive integers andn is either zero or positive integers. The element equations for the free vibration problems of the FG-GPLRC cylindrical microshell can be obtained by introducing Equations (33)- (35) into Equation (32) and then using Hamilton's principle (i.e., δI = 0), which leads to the following equation: where K Assembling the element stiffness matrix and the element mass matrix for each microlayer leads to the structure stiffness matrix and the structure mass matrix as follows: Equation (42) represents the system equations for analyzing the free vibration behavior of a simply-supported FG-GPLRC cylindrical microshell. A nontrivial solution of Equation (42) exists if and only if the determinant of its coefficient matrix vanishes. Finally, the natural frequencies of the FG-GPLRC cylindrical microshell for fixed values of the wave number pair (m,n) can be obtained by solving the following equation as follows:

Numerical Examples
Based on the CCST, the two-node and three-node Hermitian C 2 FEMs are developed as discussed above. Then, they are applied to investigate the free vibration characteristics of simply-supported FG-GPLRC cylindrical microshells. In addition, an h-convergence technique is adopted for this work to achieve the convergent solutions of each Hermitian C 2 FEM.

Homogeneous Isotropic SWCNTs
No benchmark solutions for the free vibration problem of simply-supported FG cylindrical microshells are reported in the literature. Thus, we compare the solutions of natural frequencies obtained using the Hermitian C 2 FEMs with the accurate solutions of the natural frequencies for homogeneous isotropic SWCNTs reported in the literature, following Zeighampour and Beni [47], by setting the material property gradient index at κ p = 0 and the material length scale parameter atl/h = 0 and 1. Table 1 shows the convergence and validation studies of the two-node Lagrangian C 0 [44], Hermitian C 1 [45], and Hermitian C 2 FEMs for the natural frequencies of simplysupported SWCNTs. The material properties of the SWCNT are E = 1.06 Tpa, υ = 0.3, and ρ = 2300 kg/m 3 . The geometric parameters of the SWCNT are R = 2.32 nm, L/R = 5, and h/R = 0.05. The wave number pairs are (m,n) = (1, 1), (2, 2), (3,3), (4,4), and (5,5). In addition, a dimensionless frequency parameter is defined as ω = ωR ρ/E.   It can be seen in Table 1 that the convergence rate for various FEMs is arranged in descending order: Hermitian C 2 FEM > Hermitian C 1 FEM > Lagrangian C 0 FEM. The convergent solutions obtained using various FEMs are identical to one another. In the cases of macroshells (l/h = 0), the relative errors between the solutions obtained using the current Hermitian C 2 FEM and Zeighampour and Beni's solutions obtained using the MCST-based LCST are 0.2%, 1.6%, 3.0%, 3.0%, and 6.9% for the wave number pairs (m,n) = (1, 1), (2, 2), (3,3), (4,4), and (5, 5), respectively. These deviations between them are mainly due to shear deformation effects and thickness stretching effects. The relative errors for high-frequency vibration modes are greater than those for low-frequency vibration modes, which indicates that these effects on the natural frequency are more significant for high-frequency vibration modes than for low-frequency vibration modes. In the cases of microshells (l/h = 1), the relative errors increase to 7.7%, 10.5%, 5.0%, 4.1%, and 10.0% for the wave number pairs (m,n) = (1, 1), (2, 2), (3,3), (4,4), and (5, 5), respectively. These deviations are mainly due to 3D couple-stress tensor effects, in addition to shear deformation effects and thickness stretching effects.

FG Cylindrical Macroshells
In this section, we compare the solutions of natural frequencies obtained using the twonode and three-node Hermitian C 2 FEMs with the exact and approximate 3D solutions of the natural frequencies for FG cylindrical macroshells reported in the literature by assigning the material length scale parameter a value of zero. Table 2 shows the convergence and validation studies of the two-node and threenode Hermitian C 2 FEMs for the natural frequencies of simply-supported FG cylindrical macroshells. The FG cylindrical macroshell of interest is considered to be composed of nickle and stainless steel. The inner surface of the macroshell is completely composed of nickel, and the outer surface of the macroshell is completely composed of stainless steel. The material properties of the macroshell are assumed to obey a power-law distribution of the volume fractions of the constituents varying in the thickness direction. The effective material properties are estimated using a rule of mixtures which are expressed as follows: P e f f = P n (1 − Γ ss ) + P ss Γ ss (44) where P represents some specific material properties, including Young's modulus, Poisson's ratio, and mass density; the subscripts eff, n, and ss denote the FG material, nickel, and stainless steel, respectively; Γ ss is the volume fraction function of stainless steel varying in the thickness direction, and Γ ss = [(r − a)/h] κ p , where κ p is defined as the material property gradient index; and the material properties of the constituents (i.e., nickel and stainless steel), following Loy et al. [48], are listed as: E ss = 2.07788 × 10 11 Pa, υ ss = 0.317756, ρ ss = 8166 kg/m 3 ; E n = 2.05098 × 10 11 Pa, υ n = 0.3100, ρ n = 8900 kg/m 3 . Table 2 shows the convergence and validation studies for the natural frequencies [ω/(2π)] (Hz) of a simply-supported FG cylindrical macroshell obtained using the two-node and threenode Hermitian C 2 FEMs. The relevant parameters in this analysis are κ p = 0.5, 1, and 30; the wave number pairs arem = 1 andn = 0-10; and L/a = 20, and h/a = 0.002.
It can be seen in Table 2 that the convergent solutions of natural frequencies are obtained when n l = 4 is used. The results also show that the convergence rate of the three-node Hermitian C 2 FEM is faster than that of the two-node Hermitian C 2 FEM. The convergent solutions of natural frequencies of the FG macroshell are in excellent agreement with the solutions obtained by Liu et al. [49] with the state space method and by Loy et al. [48] with the Rayleigh-Ritz method based on the LCST. It is noted that the natural frequencies of the FG macroshell decrease when the value of κ p becomes greater. The above phenomenon observed is because an increase in the value of κ p indicates both a decrease in the volume fraction of stainless steel and an increase in the volume fraction of nickel, leading to the overall stiffness of the FG macroshell decreasing and the total mass of the FG macroshell increasing, which in turn decreases the natural frequencies of the macroshell. In addition, the fundamental frequency of the very thin FG macroshell (R/h = 500.5) occurs when the wave number pair is (m,n)= (1, 3). Because of the excellent performance of the three-node Hermitian C 2 FEM shown above, it is selected for the following free vibration analyses of laminated GPLRC cylindrical microshells and FG-GPLRC cylindrical microshells.

Laminated GPLRC Cylindrical Microshells
In this section, we apply the three-node Hermitian C 2 FEM to investigate the free vibration characteristics of simply-supported, laminated GPLRC cylindrical macroshells (l = 0) and microshells (l = 0). Epoxy is used as the matrix material, and GPLs are used as the reinforcement material, which is randomly dispersed in the thickness direction. Five layer-wise constant functions of the GPL distribution patterns are considered: UD and FG V-type, A-type, X-type, and O-type distributions. The relationship between the GPL volume fraction and its weight fraction for the mth-layer, in the light of Song et al. [50,51], is expressed in the following forms: where the superscript m denotes the mth-layer, and m = 1 − n l , with n l denoting the total number of the layers constituting the laminated cylindrical shell. In this work, following Liu et al. [49], the value of n l is assigned as n l = 20. Further, W (m) GPL represents the GPL weight fraction for the mth-layer.
When a specific value of the total weight fraction (W * GPL ) of the laminated GPLRC cylindrical shell is given, the GPL weight fractions of the mth-layer for different GPL distribution patterns are expressed as follows: where m = 1 − n l , with m being counted from the bottom. The effective Young's modulus is estimated using the Halpin-Tsai model and is given in Equations (1)- (7), and the effective Poisson's ratio and the effective mass density are estimated using the rule of mixtures and are given in Equations (8) and (9), respectively. Unless otherwise stated, the material properties and geometric parameters of GPL, following Liu et al. [48], are given as follows: E GPL = 1.01 TPa, ρ GPL = 1.06 g/cm 3 = 1060 kg/m 3 , υ GPL = 0.186; (49) E epoxy = 3.0 GPa, ρ epoxy = 1.2 g/cm 3 = 1200 kg/m 3 , υ epoxy = 0.34; (50) (L x ) GPL = 2.5 µm, L y GPL = 1.5 µm,h GPL = 1.5 nm, W * GPL = 1.5%.
A nondimensional natural frequency parameter ω is defined as the same form as that used in Liu et al. [49] and is given as: (52) Table 3 shows the solutions obtained using the three-node Hermitian C 2 FEM for the first three dimensionless natural frequencies of the laminated GPLRC cylindrical macroshells l /h = 0 and microshells l /h = 0 for different values of the half-wave number,m, withm = 2, 5, and 10. The geometrical parameters of the laminated GPLRC cylindrical shell under observation are L/R = 2π/3; R/h = 3/2. It can be seen in Table 3 that in the cases of macroshells l /h = 0 , the solutions of dimensionless natural frequencies obtained using the Hermitian C 2 FEM are in excellent agreement with the solutions obtained by Liu et al. [49] with the state space method. The relative errors for the cases considered in Table 2 are less than 0.05%. The results also show that the dimensionless natural frequency increases when thel/h ratio increases. The above results are because an increase in the material length scale parameter causes the GPLRC microshell to become stiffer, in turn increasing its natural frequency. The results also show that in the cases of the total weight fraction W * GPL = 0.015, the dimensionless natural frequency of laminated GPLRC shells is much greater than that of homogeneous epoxy shells. For a particular case ofl/h = 0.4, and (m,n) = (2, 1), the dimensionless fundamental frequencies of laminated GPLRC microshells for different GPL distribution patterns are arranged following the descending order as FG V-type > FG X-type > UD > FG O-type > FG A-type distributions. The dimensionless fundamental frequencies of the laminated GPLRC shells for various GPL distribution patterns: FG V-type, FG X-type, UD, FG O-type, and FG A-type distributions are 2.464, 2.457, 2.452, 2.362, and 2.154 times the dimensionless fundamental frequencies of homogeneous epoxy shells, respectively.

FG-GPLRC Cylindrical Microshells
This section applies the three-node Hermitian C 2 FEM to examine the free vibration characteristics of simply-supported FG-GPLRC cylindrical microshells. Again, epoxy is used as the matrix material, and GPLs are used as the reinforcement material. Five GPL distribution patterns varying in the thickness direction are considered: UD and FG Vtype, A-type, X-type, and O-type distributions. The relationship between the GPL volume fraction (V GPL (ζ)) and its weight fraction (W GPL (ζ)) for each surface in the thickness direction is given in Equation (47).
When a specific value of the total weight fraction (W * GPL ) of the FG cylindrical shell is given, five GPL weight fractions varying in the thickness direction for different GPL distribution patterns are expressed in Equation (10)- (14). A nondimensional natural frequency parameter ω is defined as follows: The results in Table 4 show that the convergent solutions obtained using the three-node Hermitian  Table 4 that for a moderately thick microshell (R/h = 10), the dimensionless fundamental frequencies occur at wave number pair (m,n) = (1, 1), and the effect of the material length scale parameter on the dimensionless fundamental frequencies ((m,n) = (1, 1)) is more unsignificant than that on the dimensionless natural frequencies for other wave number pairs ((m,n) =(1, 1)). In those cases, for the V-type GPL distribution pattern, the magnitude ratios of the lowest natural frequency between the FG-GPLRC microshell (l/h = 1.0) and the FG-GPLRC macroshell (l/h = 0) are 1.0013, 2.3696, 1.8649, 2.4326, 2.4018, and 2.3036 for the wave number pairs (m,n) = (1, 1), (1, 2), (2, 2), (1, 3), (2, 3), and (3, 3), respectively. The results also show that the fundamental natural frequencies ((m,n) = (1, 1)) of various GPL distribution patterns are arranged in descending order as: FG V-type > FG X-type > UD > FG O-type > FG A-type distributions. This is because the more the GPL is dispersed away from the center of the GPLRC microshell, the stiffer the GPLRC microshell is, which in turn increases its dimensionless natural frequency. Table 3. Comparison study for the first three lowest frequency parameters ω k of laminated GPLRC cylindrical macroshells (l/h = 0) and microshells (l/h = 0) obtained using the CCST-based three-node Hermitian C 2 FCLM, for whichm = 2, 5, and 10; n l = 20; and

GPL Distribution
Patterns  In Figures 3-9, we present a parametric study related to some key effects on the lowest dimensionless natural frequencies of simply-supported FG-GPLRC cylindrical microshells, including the effects of the ratios ofl/h, R/h, and L/R, the GPL weight fraction, and the length-to-thickness ratio of the GPL, and different GPL distribution patterns. The material properties of epoxy and GPLs and the geometric dimensions of GPLs are given in Equations (49)- (51), the value of the material length scale parameter is set atl = 2l = 17.6 × 10 −6 m, and the dimensionless natural frequency ω is defined as the same form in Equation (45).
for different GPL distribution patterns can be arranged in descending order as FG V-type > FG X-type > UD > FG O-type > FG A-type distributions. Likewise, they are in descending order as FG X-type > (FG A-type, UD) > (FG O-type, FG V-type) distributions for the case of ( ), m n = (1, 2) in Figure 5b. The results also show that the effect of the material length scale parameter on the lowest dimensionless natural frequencies for the vibration modes with ( ), m n = (1, 2), is more significant than that for those with ( ),            Figure 8b. It can be seen in Figure 8a,b that the dimensionless fundamental frequency increases when the weight fraction of the GPL increases. This observation indicates that an increase in the weight fraction of GPL leads to an increase in the overall stiffness of the microshell and a decrease in the total mass of the microshell, increasing its dimensionless fundamental frequency. Figure 9 shows the variations in the dimensionless fundamental frequency of an FG-GPLRC cylindrical microshell with the ratios of the length-to-thickness of GPLs, for dif- = 0.01. Figure 9 shows that the dimensionless fundamental frequency increases when the length-to-thickness ratio increases. This observation indicates that when the values of the length, the width, and the GPL weight fraction are fixed, an increase in the length-tothickness ratio leads to an increase in the total contact surface area between GPL reinforcements and epoxy matrix and an increase in the overall stiffness of the microshell, which in turn increases the dimensionless fundamental frequency.

Concluding Remarks
Within a framework of the CCST, we developed a 3D weak formulation of the CCSTbased Hermitian C 2 FEM for analyzing the free vibration characteristics of simply-supported FG-GPLRC cylindrical microshells using Hamilton's principle. In addition, a parametric study related to some key effects on the lowest natural frequencies of simply-supported FG-GPLRC cylindrical microshells was carried out, including the material length scale parameter, the GPL distribution patterns, the GPL weight fractions, the length-tomid-surface radius and mid-surface radius-to-thickness ratio of the microshell, and the length-to-thickness ratio of GPLs. Some conclusions drawn from the parametric study were summarized as follows: 1. The CCST-based Hermitian C 2 FEM was validated to be accurate and to converge rapidly by comparing the solutions it produced with the 3D exact and the quasi-3D solutions for FG cylindrical macroshell reported in the literature by assigning the length scale parameter to a value of zero. 2. In numerical examples, the results showed that: (a) An increase in the material length scale parameter caused the FG cylindrical microshell to become stiffer, increasing the lowest natural frequency of the microshell. The effect of the material length scale parameter on the lowest natural frequency for the high-frequency vibration modes was more significant than that for the fundamental vibration mode. (b) For moderately thick (R/h = 10) and thick (R/h = 5) cylindrical microshells, the  Figures 3 and 4 that for a fixed value ofm, when the value ofn increases from zero to seven, the dimensionless natural frequency first decreases to a relative minimum value, and then it gradually increases to a relative maximum value. The fundamental natural frequency for different GPL distribution patterns always occurs when the wave number pair is (m,n) = (1, 1) for the thick microshell (R/h = 5) and is (m,n) = (1, 2) for the thin microshell (R/h = 20).
By comparing the dimensionless natural frequency of an FG-GPLRC cylindrical microshell in Figures 3a-e and 4a-e with those of a homogeneous epoxy cylindrical microshell in Figures 3f and 4f, respectively, it is shown that the dimensionless natural frequency of a homogeneous epoxy cylindrical microshell increases significantly by adding a small amount of GPLs. For example: for the V-type GPL distribution pattern, by increasing the weight fraction of GPLs by 1%, the natural frequency of FG-GPLRC cylindrical microshells can be increased to 2.13 and 2.05 times that of the homogeneous epoxy cylindrical microshells for the mid-surface radius-to-thickness ratios R/h = 5 and R/h = 20, respectively. Figure 5a,b show that the variations in the lowest dimensionless natural frequency of an FG-GPLRC cylindrical microshell with the ratio ofl/h for different GPL distribution patterns. Various parameters are considered asl/h = 0-1; (m,n) = (1, 1) and (1,2) in Figure 5a,b, respectively; R/h = 5, and L/R = 10; and W * GPL = 0.01. It can be seen in Figure 5a that for the case of (m,n) = (1, 1), the values of the dimensionless natural frequencies for different GPL distribution patterns can be arranged in descending order as FG V-type > FG X-type > UD > FG O-type > FG A-type distributions. Likewise, they are in descending order as FG X-type > (FG A-type, UD) > (FG O-type, FG V-type) distributions for the case of (m,n) = (1, 2) in Figure 5b. The results also show that the effect of the material length scale parameter on the lowest dimensionless natural frequencies for the vibration modes with (m,n) = (1, 2), is more significant than that for those with (m,n) = (1, 1). Figures 6 and 7 show the variations in the dimensionless fundamental frequency of an FG-GPLRC cylindrical microshell with the ratios of R/h and L/R, respectively, for different GPL distribution patterns. The relevant parameters are (m,n) = (1, 1), R/h = 5-15, L/R = 10, l/h = 1.0, and W * GPL = 0.01 given in Figure 6, and (m,n) = (1, 1), R/h = 5, L/R = 2-10, l/h = 1.0, and W * GPL = 0.01 given in Figure 7. It can be seen in Figures 6 and 7 that the dimensionless fundamental frequency decreases when the microshell becomes thinner and longer. This observation indicated that an increase in the ratios of R/h and L/R causes the microshell to become softer, decreasing its dimensionless fundamental frequencies. The results also show that in the ranges of R/h = 5-15 and L/R = 10 considered in Figure 6 and of R/h = 5 and L/R = 2-10 considered in Figure 7, the fundamental frequency occurs when the wave number pair is (m,n) = (1, 1). Again, the values of the dimensionless natural frequencies for different GPL patterns can be arranged in descending order as FG V-type > FG X-type > UD > FG O-type > FG A-type distributions. Figure 8a,b show the variations in the dimensionless fundamental frequency of an FG-GPLRC thick and an FG-GPLRC thin cylindrical microshell, respectively, with the weight fraction of GPL for different GPL distribution patterns.
The relevant parameters arel/h = 1.0, R/h = 5, L/R = 5, and (m,n) = (1, 1) in Figure 8a, andl/h = 1.0, R/h = 20, L/R = 5, and (m,n) = (1, 2) in Figure 8b. It can be seen in Figure 8a,b that the dimensionless fundamental frequency increases when the weight fraction of the GPL increases. This observation indicates that an increase in the weight fraction of GPL leads to an increase in the overall stiffness of the microshell and a decrease in the total mass of the microshell, increasing its dimensionless fundamental frequency. Figure 9 shows the variations in the dimensionless fundamental frequency of an FG-GPLRC cylindrical microshell with the ratios of the length-to-thickness of GPLs, for different GPL distribution patterns. The relevant parameters are R/h = 20, L/R = 5, (m,n) = (1, 2), and W * GPL = 0.01. Figure 9 shows that the dimensionless fundamental frequency increases when the length-to-thickness ratio increases. This observation indicates that when the values of the length, the width, and the GPL weight fraction are fixed, an increase in the length-tothickness ratio leads to an increase in the total contact surface area between GPL reinforcements and epoxy matrix and an increase in the overall stiffness of the microshell, which in turn increases the dimensionless fundamental frequency.

Concluding Remarks
Within a framework of the CCST, we developed a 3D weak formulation of the CCST-based Hermitian C 2 FEM for analyzing the free vibration characteristics of simplysupported FG-GPLRC cylindrical microshells using Hamilton's principle. In addition, a parametric study related to some key effects on the lowest natural frequencies of simplysupported FG-GPLRC cylindrical microshells was carried out, including the material length scale parameter, the GPL distribution patterns, the GPL weight fractions, the length-to-mid-surface radius and mid-surface radius-to-thickness ratio of the microshell, and the length-to-thickness ratio of GPLs. Some conclusions drawn from the parametric study were summarized as follows: 1.
The CCST-based Hermitian C 2 FEM was validated to be accurate and to converge rapidly by comparing the solutions it produced with the 3D exact and the quasi-3D solutions for FG cylindrical macroshell reported in the literature by assigning the length scale parameter to a value of zero.

2.
In numerical examples, the results showed that: (a) An increase in the material length scale parameter caused the FG cylindrical microshell to become stiffer, increasing the lowest natural frequency of the microshell. The effect of the material length scale parameter on the lowest natural frequency for the high-frequency vibration modes was more significant than that for the fundamental vibration mode.
For moderately thick (R/h = 10) and thick (R/h = 5) cylindrical microshells, the fundamental mode of vibration always occurred when the wave number pair was (m,n) = (1, 1). In that case, the values of the dimensionless natural frequencies for different GPL distribution patterns can be arranged in descending order as follows: FG V-type > FG X-type > UD > FG O-type > FG A-type distributions. (c) For thin (R/h = 20) microshells, the fundamental mode of vibration always occurred when the wave number pair was (m,n) = (1, 2). In that case, the values of the dimensionless natural frequencies for different GPL distribution patterns can be arranged in descending order as follows: FG X-type > (FG A-type, UD) > (FG O-type, FG V-type) distributions. (d) By increasing the weight fraction of GPLs by 1%, the natural frequency of FG-GPLRC cylindrical microshells can be increased to more than twice that of the homogeneous cylindrical microshells. (e) The dimensionless fundamental frequency decreased when the microshells became thinner and longer. (f) The dimensionless fundamental frequency increased when the weight fraction of the GPL increased and when the length-to-thickness ratio of the GPL increased.
The 3D solutions for the free vibration characteristics of FG-GPLRC cylindrical microshells are relatively limited in the literature. Therefore, the current solutions can be used as a reference to evaluate the theoretical accuracy of various 2D size-dependent shear deformation theories of FG cylindrical microshells. Furthermore, due to the excellent performance of the CCST-based Hermitian C 2 FEM, it is recommended to analyze various mechanical behaviors of FG elastic/piezoelectric cylindrical microshells and various FG elastic/piezoelectric microshells of revolution.