Size-Dependent Free Vibration and Buckling of Three-Dimensional Graphene Foam Microshells Based on Modified Couple Stress Theory

In this research, the vibration and buckling of three-dimensional graphene foam (3D-GrF) microshells are investigated for the first time. In the microshells, three-dimensional graphene foams can distribute uniformly or non-uniformly through the thickness direction. Based on Love’s thin shell theory and the modified couple stress theory (MCST), size-dependent governing equations and corresponding boundary conditions are established through Hamilton’s principle. Then, vibration and axial buckling of 3D-GrF microshells are analyzed by employing the Navier method and Galerkin method. Results show that the graphene foam distribution type, size effect, the foam coefficient, the radius-to-thickness ratio, and the length-to-radius ratio play important roles in the mechanical characteristics of 3D-GrF microshells.

In contrast to many experimental studies on the mechanical property of 3D-GrFs, few computational and theoretical investigations have been conducted to evaluate the relationship between the macro-mechanical characteristics and the intrinsic micro/nanostructures. Jinlong et al. [19] investigated and discussed the effects of graphene layer number and 3D-GrF defects on the electrical and mechanical properties of 3D-GrFs. Based on the coarse-grained molecular dynamics simulation, the constitutive relation between 3D-GrFs and microscopic deformation mechanisms were investigated by Wang et al. [20]. The mechanical properties in both tension and compression of 3D-GrFs at macro and nanoscales were evaluated by Nieto et al. [21]. Qin et al. [3] stated that 3D-GrFs have an exceptionally high tensile strength and they are 10 times as strong as mild steel. Nautiyal et al. [22] studied the dynamic mechanical properties of 3D-GrFs by nanoindentation technique.
As material size changes from macroscale to nano/microscale, size effect on mechanical properties of materials should be taken into account. It is reported that size effects in micro/nanostructures have been experimentally observed [23][24][25]. Due to the exclusion of size effect in the classical theory (CT), several non-classical theories were developed for analysis of micro/nanostructures. One of the non-classical theories incorporating size effect is the modified couple stress theory [26,27]. This theory has been employed in many aspects to interpret the size effect in micro/nanostructures [28][29][30][31][32][33][34].
Micro/nanoshells are an important and widely used form of micro/nanostructures in various engineering fields. Gholami et al. [35] analyzed the dynamic stability and axial buckling of functionally graded cylindrical microshells. Based on the strain gradient elasticity theory, Zhang et al. [36] developed a shear deformable functionally graded microshell model. Ghayesh and Farokhi [37] studied the nonlinear dynamical characteristics of doubly curved shallow microshells. SafarPour et al. [38] investigated the influences of various temperature distributions on the vibration of functionally graded rotating cylindrical microshells. Afterward, the buckling behavior of functionally graded sandwich microshells under axial loads was studied by Zeighampour and Shojaeian [39]. Wang et al. [40] studied nonlinear vibrations of cylindrical nanoshells conveying fluid in the framework of the surface stress elasticity theory.
To date, no work has been published on the mechanical properties of 3D-GrF microshells. In this paper, we aim to conduct size-dependent buckling and free vibration analysis of 3D-GrF cylindrical microshells. The couple stress theory and Love's thin shell theory are used to derive the governing equations. After that, natural frequencies and critical buckling loads are solved by adopting the Navier method and Galerkin method. Finally, the results are illustrated for different parameters.

Modified Couple Stress Theory
Based on the modified couple stress theory, the strain energy U s in a deformed elastic body occupying a volume Ω is [27] U s = 1 2 Ω (m : χ + σ : ε)dV (1) where m represents the deviatoric part of the couple stress tensor; χ is the symmetric curvature tensor; σ is the Cauchy stress tensor; ε is the strain tensor. These tensors are defined by [41]: where µ and λ represent Lamé constants; u represents the displacement vector; l represent a material length scale parameter; θ is the rotation vector given by The classical and higher-order strains are derived as [42,43] ε (i) Materials 2019, 12, 729 where η , and u (k) represent the physical components of higher-order displacement gradient η k i j , displacement gradient ε i j , and displacement vector u k , respectively; g ii and Γ i jk denote an individual diagonal covariant component of the Euclidean metric tensor and Christoffel symbols of the second kind, respectively. The underscores under the indices denote no summation over indices. In cylindrical coordinates, the components of metric tensor and Christoffel symbols are expressed as Substituting Equations (9) and (10) into Equation (8) gives

3D-GrF Circular Cylindrical Microshell
In Figure 1, a 3D-GrF microshell with the thickness h, the middle-plane radius R, and the length L is shown. u(x, θ, t), v(x, θ, t), and w(x, θ, t) represent the in-plane and transverse displacements of points at the middle plane; N 0 xx represents the axial load applied to the microshell.

3D-GrF Circular Cylindrical Microshell
In Figure 1, a 3D-GrF microshell with the thickness h, the middle-plane radius R, and the length L is shown. u(x, θ, t), v(x, θ, t), and w(x, θ, t) represent the in-plane and transverse displacements of points at the middle plane; 0 xx N represents the axial load applied to the microshell. Three types of foam distributions in the thickness direction were considered, as shown in Figure 2. Herein, Figure 2a,b depict non-uniform foam distribution while Figure 2c shows uniform foam distribution. They are denoted by 3D-GrF-I, 3D-GrF-II, and 3D-GrF-U, respectively. As shown in Figure 2, the largest foams are located on the mid-plane for 3D-GrF-I while on the top and bottom surfaces for 3D-GrF-II, leading to the variations of material properties given in Equations (12)-(14) Three types of foam distributions in the thickness direction were considered, as shown in Figure 2. Herein, Figure 2a,b depict non-uniform foam distribution while Figure 2c shows uniform foam distribution. They are denoted by 3D-GrF-I, 3D-GrF-II, and 3D-GrF-U, respectively. As shown in Figure 2, the largest foams are located on the mid-plane for 3D-GrF-I while on the top and bottom surfaces for 3D-GrF-II, leading to the variations of material properties given in Equations (12)- (14) for 3D-GrF-I, and Equations (15)-(17) for 3D-GrF-II. Material properties of 3D-GrF-U are described in Equations (18)- (20). They are given by [44][45][46]: where E(z), ρ(z), and G(z) are general Young's modules, mass density, and shear modules of the 3D-GrF microshell, respectively; E 1 , G 1 , and ρ 1 represent corresponding properties of solid graphenes without internal foams; κ 0 and κ m represent coefficients of foams and mass density for 3D-GrF-I, respectively; κ * 0 and κ * m are corresponding coefficients for 3D-GrF-II; ϑ and ϑ are corresponding coefficients for 3D-GrF-U. Thereinto, shear modulus G 1 is calculated by where ν represents Poisson's ratio. for 3D-GrF-I, and Equations (15)-(17) for 3D-GrF-II. Material properties of 3D-GrF-U are described in Equations (18)- (20). They are given by [44][45][46]: where E(z), ρ(z), and G(z) are general Young's modules, mass density, and shear modules of the 3D-GrF microshell, respectively; E1, G1, and ρ1 represent corresponding properties of solid graphenes without internal foams; 0 κ and m κ represent coefficients of foams and mass density for 3D-GrF-I, respectively; 0 κ * and m κ * are corresponding coefficients for 3D-GrF-II; ϑ and ' ϑ are corresponding coefficients for 3D-GrF-U. Thereinto, shear modulus G1 is calculated by ( ) where ν represents Poisson's ratio.
(a) The typical mechanical property of 3D-GrFs [3,20,21], shown in Equation (22), is employed to establish the relationships in Equation (23) between mass density coefficients and foam coefficients for different foam distributions: The masses of all 3D-GrF microshells with varying foams are set to be equivalent, namely which can be used to determine 0 κ * and ϑ with a given value of 0 κ , as tabulated in Table 1. It is seen that 0 κ * rises dramatically with the increase of 0 κ . When 0 κ reaches 0.65, 0 κ * is close to the upper limit ( 0 0.9976 ). Therefore, the selected range of 0 κ ∈ [0, 0.65] is applied in the following numerical calculations. The typical mechanical property of 3D-GrFs [3,20,21], shown in Equation (22), is employed to establish the relationships in Equation (23) between mass density coefficients and foam coefficients for different foam distributions: The masses of all 3D-GrF microshells with varying foams are set to be equivalent, namely which can be used to determine κ * 0 and ϑ with a given value of κ 0 , as tabulated in Table 1. It is seen that κ * 0 rises dramatically with the increase of κ 0 . When κ 0 reaches 0.65, κ * 0 is close to the upper limit (κ * 0 = 0.9976). Therefore, the selected range of κ 0 ∈ [0, 0.65] is applied in the following numerical calculations. The displacement fields, based on Love's thin shell theory, are expressed as [47] where t is time, and u 1 (x, θ, z, t), u 2 (x, θ, z, t), and u 3 (x, θ, z, t) are displacements of an arbitrary point of the microshell along x-, θ-, and z-axes, respectively. By substituting Equations (25)- (27) and Equation (10) into Equation (7), strain-displacement relations are expressed as: Nonzero components of χ can be obtained by substituting Equations (25)- (27) into Equation (11) and using Equation (5): According to Equation (1), the strain energy can be written as The non-classical and classical moments and forces are given by The kinetic energy of the microshell is The work performed by axial load N 0 xx applied on the middle surface of the 3D-GrF microshell is given by [48] Using Hamilton's principle [49][50][51][52][53][54] and applying Equations (37), (39) and (40) in Equation (41), the governing equations of motion for the 3D-GrF cylindrical microshell are −D 1,0 where Boundary conditions at the edges with x = constant are Boundary conditions at the edges with θ = constant are:

Navier Solution
For simply supported-simply supported (SS-SS) 3D-GrF cylindrical microshells, the displacement functions using the Navier method can be expressed as: v(x, θ, t) = V mn sin mπx L sin(nθ)e iωt (56) w(x, θ, t) = W mn sin mπx L cos(nθ)e iωt (57) where U mn , V mn , and W mn represent the displacement amplitude components; m and n are mode numbers; and ω is natural circular frequency of the 3D-GrF microshell. Obviously, displacement functions in Navier solution procedure satisfy the SS-SS boundary condition. Substituting Equations (55)-(57) into Equations (42)- (44) and then removing trigonometric functions lead to where d = [U mn , V mn , W mn ] T ; K, K g , and M denote stiffness matrix, geometric stiffness matrix, and mass matrix, respectively. The non-zero elements in the above matrices are given in the Appendix A.
Neglecting axial load N 0 xx , Equation (58) is reduced to the following eigenvalue problem of the 3D-GrF microshell: If disregarding inertia terms and assuming N 0 xx = −P, Equation (58) is reduced to the equilibrium equations of a static buckling problem of the 3D-GrF microshell: where P represents buckling load. Therefore, the critical (minimum) buckling load P cr and vibration frequencies of the 3D-GrF microshell are obtained by solving the above eigenvalue problems [55][56][57].

Results and Discussion
In order to verify the correctness of the derivation in this paper, we chose an SS-SS isotropic homogeneous shell to make a comparison without considering the small-scale effect (l = 0). The system parameters used are as follows: ρ = 2300 kg/m 3 , ν = 0.3, E = 1.06 TPa, R = 2.32 nm, L/R = 5. As can be seen in Table 3, very good agreement was achieved.
Considering the scale effect, dimensionless natural frequencies of an isotropic homogeneous cylindrical nanoshell were calculated, and are compared with the literature in Table 4. Here the parameters used are the same as those in Table 3. It was found that the present results match those given by Beni et al. [42] well, bespeaking the validity of the present study. To make further comparisons, we considered a homogeneous cylindrical shell with the C-C boundary condition. The comparison result is listed in Table 5. It can be seen that the present results agree well with those in the literature [64].  For all six of the boundary conditions, one can see the natural frequencies first decreased and then increased as circumferential number n increased. The minimum frequency occurred at n = 2 for the C-F boundary condition and occurred at n = 3 for the other boundary conditions. In addition, the natural frequency of the C-F microshell was the lowest while that of the F-F one was the highest. When the circumferential wave number was greater than 5, natural frequencies under various boundary conditions tended toward the same value. This result shows that the boundary condition effect was closely associated with the circumferential wave number.
Hereinafter, we conducted free vibration and buckling analyses of 3D-GrF cylindrical microshells. If not specified, the geometrical and material properties of the 3D-GrF microshell are: boundary conditions, one can see the natural frequencies first decreased and then increased as circumferential number n increased. The minimum frequency occurred at n = 2 for the C-F boundary condition and occurred at n = 3 for the other boundary conditions. In addition, the natural frequency of the C-F microshell was the lowest while that of the F-F one was the highest. When the circumferential wave number was greater than 5, natural frequencies under various boundary conditions tended toward the same value. This result shows that the boundary condition effect was closely associated with the circumferential wave number.       Based on the modified couple stress theory and classical continuum theory, Figure 6 shows the natural frequency versus dimensionless length scale parameter for the SS-SS 3D-GrF microshell, where n = 3 and κ0 = 0.2. One can see that the natural frequency from the modified couple stress theory increased while that from the classical theory did not change with length scale parameter. This is because size effect tends to increase the stiffness of the 3D-GrF microshell. Nevertheless, the classical shell theory fails to incorporate this effect and thus produces inaccurate results. Based on the modified couple stress theory and classical continuum theory, Figure 6 shows the natural frequency versus dimensionless length scale parameter for the SS-SS 3D-GrF microshell, where n = 3 and κ 0 = 0.2. One can see that the natural frequency from the modified couple stress theory increased while that from the classical theory did not change with length scale parameter. This is because size effect tends to increase the stiffness of the 3D-GrF microshell. Nevertheless, the classical shell theory fails to incorporate this effect and thus produces inaccurate results.  Figure 7 gives the variation of natural frequency against the foam coefficient for different foam distributions. The results show that the increase in the foam coefficient led to a decrease in the natural frequencies of 3D-GrF cylindrical microshells. When the foam coefficient was small, the 3D-GrF-I microshell had the highest natural frequency while the 3D-GrF-II one had the lowest natural frequency. With the increase in the foam coefficient, however, the natural frequency of the  Figure 7 gives the variation of natural frequency against the foam coefficient for different foam distributions. The results show that the increase in the foam coefficient led to a decrease in the natural frequencies of 3D-GrF cylindrical microshells. When the foam coefficient was small, the 3D-GrF-I microshell had the highest natural frequency while the 3D-GrF-II one had the lowest natural frequency. With the increase in the foam coefficient, however, the natural frequency of the 3D-GrF-U microshell dropped faster than that of the 3D-GrF-II one. At last, the 3D-GrF-U microshell had the lowest natural frequency among the three types of foam distribution, except in the C-F boundary condition.  Figure 7 gives the variation of natural frequency against the foam coefficient for different foam distributions. The results show that the increase in the foam coefficient led to a decrease in the natural frequencies of 3D-GrF cylindrical microshells. When the foam coefficient was small, the 3D-GrF-I microshell had the highest natural frequency while the 3D-GrF-II one had the lowest natural frequency. With the increase in the foam coefficient, however, the natural frequency of the 3D-GrF-U microshell dropped faster than that of the 3D-GrF-II one. At last, the 3D-GrF-U microshell had the lowest natural frequency among the three types of foam distribution, except in the C-F boundary condition. The variations of natural frequency against the length-to-radius ratio for different foam distributions are illustrated in Figures 8-10, where n = 3 and κ0 = 0.2. It was found that the natural frequencies kept decreasing with the length-to-radius ratio. It can also be seen that the natural frequencies of 3D-GrF microshells under various boundary conditions tended toward the same value as the length-to-radius ratio increased. The variations of natural frequency against the length-to-radius ratio for different foam distributions are illustrated in Figures 8-10, where n = 3 and κ 0 = 0.2. It was found that the natural frequencies kept decreasing with the length-to-radius ratio. It can also be seen that the natural frequencies of 3D-GrF microshells under various boundary conditions tended toward the same value as the length-to-radius ratio increased. The variations of natural frequency against the length-to-radius ratio for different foam distributions are illustrated in Figures 8-10, where n = 3 and κ0 = 0.2. It was found that the natural frequencies kept decreasing with the length-to-radius ratio. It can also be seen that the natural frequencies of 3D-GrF microshells under various boundary conditions tended toward the same value as the length-to-radius ratio increased.        Table 6. Natural frequency ω (MHz) of 3D-GrF microshell with different radius-to-thickness ratios  Tables 6-8 show the natural frequency versus radius-to-thickness ratio for different foam distributions, where n = 3 and κ 0 = 0.2. It is clear that natural frequencies decreased as the radius-to-thickness ratio rose for all the boundary conditions. Table 6. Natural frequency ω (MHz) of 3D-GrF microshell with different radius-to-thickness ratios (3D-GrF-I).  Table 7. Natural frequency ω (MHz) of 3D-GrF microshell with different radius-to-thickness ratios (3D-GrF-II).  Table 8. Natural frequency ω (MHz) of 3D-GrF microshell with different radius-to-thickness ratios (3D-GrF-U).  Table 9 lists the variation of buckling load against circumferential wave number for the SS-SS 3D-GrF microshell, where L/R = 3 and κ 0 = 0.2. Results show the buckling loads initially decreased and then increased with the circumferential wave number. It can be seen that the critical buckling load P cr occurred at n = 3. Therefore, mode (m = 1, n = 3) is chosen as a representative mode in the next buckling studies.

R/h SS-SS C-SS C-C C-F F-F F-SS
The buckling load against radius-to-thickness ratio for the SS-SS 3D-GrF microshell is listed in Table 10, where different foam distributions are taken into account. It is shown that the larger radius-to-thickness ratio led to the lower buckling load. One can also see that the buckling loads evaluated by the modified couple stress theory were higher than those evaluated by the classical theory, showing that the classical theory underestimated the buckling loads of 3D-GrF microshells.  Figure 11 examines the critical buckling load against the foam coefficient for SS-SS 3D-GrF microshells with different foam distributions. One can see that the higher foam coefficient led to the lower critical buckling load. Additionally, the 3D-GrF-I microshell had the highest critical buckling load. This shows that compact inner and outer surfaces and sparse mid-plane can lead to higher structural stiffness of 3D-GrF microshells.   Figure 11 examines the critical buckling load against the foam coefficient for SS-SS 3D-GrF microshells with different foam distributions. One can see that the higher foam coefficient led to the lower critical buckling load. Additionally, the 3D-GrF-I microshell had the highest critical buckling load. This shows that compact inner and outer surfaces and sparse mid-plane can lead to higher structural stiffness of 3D-GrF microshells.  Figure 12 depicts the buckling load against length-to-radius ratio for SS-SS 3D-GrF microshells with different foam distributions. With the increase of length-to-radius ratio, buckling load decreased initially and then increased. The lowest value of buckling loads occurred at nearby L/R = 3.5 for all types of foam distribution.  Figure 12 depicts the buckling load against length-to-radius ratio for SS-SS 3D-GrF microshells with different foam distributions. With the increase of length-to-radius ratio, buckling load decreased initially and then increased. The lowest value of buckling loads occurred at nearby L/R = 3.5 for all types of foam distribution.

Conclusions
In this study, size-dependent free vibration and buckling of cylindrical 3D-GrF microshells were investigated. Formulation of the 3D-GrF microshell was extracted using the modified couple stress theory and the thin shell theory. The Navier method and the Galerkin method were utilized to get critical buckling loads and natural frequencies of 3D-GrF microshells. Results indicate that the foam distribution and foam coefficient play important roles in the buckling and vibration behavior of 3D-GrF microshells. An increase in the foam coefficient led to the smaller natural frequency and critical buckling load. When the foam coefficient was small, the 3D-GrF-I microshell had the highest while its 3D-GrF-II counterpart had the lowest natural frequency and critical buckling load; when the foam coefficient was large, the 3D-GrF-U microshell had the lowest natural frequency, except in the C-F boundary condition. Moreover, compact inner and outer surfaces and sparse mid-plane can the enhance structural stiffness of 3D-GrF microshells.

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

Appendix A
The non-zero elements in the matrices of Equation (58)