Vibration Analysis of Vacancy Defected Graphene Sheets by Monte Carlo Based Finite Element Method

The stochastic distributed placement of vacancy defects has evident effects on graphene mechanical property, which is a crucial and challenged issue in the field of nanomaterial. Different from the molecular dynamic theory and continuum mechanics theory, the Monte Carlo based finite element method (MC-FEM) was proposed and performed to simulate vibration behavior of vacancy defected graphene. Based on the Monte Carlo simulation, the difficulties in random distributed location of vacancy defects were well overcome. The beam element was chosen to represent the exact atomic lattice of the graphene. The results of MC-FEM have a satisfied agreement with that in the reported references. The natural frequencies in the certain vibration mode were captured to observe the mechanical property of vacancy defected graphene sheets. The discussion about the parameters corresponding with geometry and material property was accomplished by probability theory and mathematical statistics.


Introduction
For the parameters corresponding to material properties of defect-free graphene monolayer, Young's modulus and intrinsic strength are 1.0 TPa and 130 GPa [1]. Large deviation in simulations and experiments has been found, which is attributed to the presence and uncertainty of defects in the nanotube structure [2,3]. Various complicated defects exist in the graphene sheet, and some defects in the atomic structure can even deteriorate the mechanical and electronic properties of graphene materials [4]. The research of vacancy defects in graphene sheets is very necessary, by which the uncertainties and fluctuation in mechanical property of graphene sheets can be reasonably explained and comprehensively understood.
Graphene, as a two-dimensional structure, was introduced in 2004 [5]. The main tough task in the analysis of the mechanical behavior of graphene is its small size [6][7][8][9]. Obviously, it is difficult to conduct physical experiments in nanoscale. Analytical and numerical methods are available alternatives in nanomaterial research. These methods can be divided into two categories: atomistic-based methods and continuum mechanics theories. Atomistic-based methods include the molecular dynamics simulation (MD) [10], tight-binding molecular dynamics [11] and density function theory [12]. Some size dependent non-classical continuum theories have been studied in the research field of graphene sheets, which includes the strain gradient theory [13], modified couple stress theory [14], and nonlocal elasticity theory [15]. However, atomistic based methods are computationally expensive when the amount of atoms in the system is large. In continuum theories, the defects or uncertainties in the nanostructure are difficult to be taken into consideration. According to the overview of the existing literature, it can be found that the influence of the vacancy defects on the mechanical behavior of nanostructures has been studied in very few papers.
Attempts and struggles for a deeper understanding of the influence of defects in graphene sheets have been undertaken in this challenging field. For instance, the effect of vacancy defects on Young's modulus of graphene sheets was discussed by molecular dynamics [16]. The Stone-Wales defects in tensile behavior of graphene sheets and carbon nanotubes were accomplished by the same method [17]. In addition, molecular dynamics simulations are applied in studying mechanical property of defective single-layered graphene sheets [18]. The non-local elastic theory is also applied in vibration analysis of defective graphene sheets. However, the effects of structural defects on the vibration behavior of graphene sheets has rarely been investigated [19].
The motivation of this study is to effectively analyze vibration behavior of defected graphene sheets when the vacancy defects are randomly distributed. Since the physical experiments of vacancy defected graphene sheets are difficult and expensive, numerical simulation methods are another promising supplement that need to be developed. Besides, it is necessary to find an effective index to express and observe the effects of vacancy defects in vibration behavior of graphene sheets.
In the free vibration, natural frequencies of structures are the essential parameters, which are directly corresponding with the geometrical and material characteristics. The change of the vacancy distribution and location can cause fluctuation in natural frequencies of graphene sheets. The natural frequencies of vibration modes are very sensitive to the geometrical defects. Therefore, natural frequency is a good indicator to discuss the effects of vacancy defects in the graphene sheets.
Combining the Monte Carlo simulation with the finite element method to effectively simulate the random dispersed vacancy defects in graphene sheets is a perspective and feasible method. The Monte Carlo simulation as a sophisticated sampling method has been widely applied in the research of phase transition and magnetism of graphene sheets [20][21][22][23]. When the number of sampling is sufficient, it can reach a satisfied accuracy in numerical computations. It is common to take the results of the Monte Carlo simulation as the exact solution or comparison standard [24,25]. By combining the Monte Carlo method with the finite element model, an effective and feasible method is proposed in this paper. The difficulties in placement uncertainty, irregularity, and stochastic location are solved in the Monte Carlo simulation. The beam finite element model is well represented the specific hexagon microstructure of graphene sheets.
The remaining part of the paper is organized as follows: a clear introduction about the beam finite element model for representing graphene hexagon microstructure is presented. The Timoshenko beam theory for free vibration analysis is explicitly derived. Besides, the implementation of Monte Carlo based finite element method is expressed in computational program after the validation of deterministic beam finite element model. Good agreement is reached in the present model and reported literatures [26][27][28][29][30][31][32][33][34] for the vibration analysis of graphene sheets. The discussion extensively consists of the effect of amount of vacancy defects and geometrical and material parameters in the natural frequencies of vacancy defected graphene sheets. A short conclusion of this paper is provided.

Graphene Sheets
Carbon atoms in graphene are bonded together with covalent bonds forming a hexagonal 2D lattice. The parameters to describe the bonds are the bond length and bond angle. The displacement of individual atoms under external forces is constrained by the bonds. In other words, the mechanical deformation of graphene sheets is the result of the interactions between the bonds. Graphene sheets are supposed as planar-frame structures, the bonds act as connecting load-carrying elements, and atoms are assumed as joints of the connecting elements. Then, classical structural mechanics methods such as the finite element method can be applied in the vibration analysis. Figure 1 depicts planar-frame model of graphene sheets by simulating the constitutional element carbon atoms as nodes and C-C action as linking beams. The entire graphene sheet's lattice is modeled by periodic and regular hexagons. C-C action as linking beams. The entire graphene sheet's lattice is modeled by periodic and regular hexagons. The elastic geometrical properties of the beam elements have been derived using an energy relationship between molecular mechanics and continuum mechanics developed in [35]. According to the computation, the diameter d , Young's modulus E , and shear modulus G of the beam elements representing the C-C bonds are derived from where , , r k k k   are the bond sketching, bond bending and torsional resistance force constants, respectively. In the finite element model of graphene sheets, the length of C-C bond is corresponding to the length of beam in planar-frame structure, wall thickness is related with the diameter of circular solid cross section of the beam elements, as in Figure 1.
The numerical simulation model of graphene sheets was created by the ANSYS parameter design language. For the modeling of the C-C bonds, the three-dimensional (3D) elastic BEAM188 element was used [20]. BEAM188 is suitable for analyzing slender to moderately stubby/thick beam structures. The element is based on the Timoshenko beam theory which includes first-order sheardeformation effects. The element is a linear, quadratic, and cubic two-node beam element in 3D. For each node, it has six degrees of freedom, which include translations in the x, y, and z directions and rotations around the x, y, and z axes. BEAM 188 is well-suited for linear, large rotation, and/or large strain nonlinear applications.

Finite Element Model in Vibration Analysis
The beam elements applied in graphene sheets are based on the Timoshenko beam theory, which is a first order shear deformation theory. In this theory, the transverse shear strain is constant through the cross section, which remains plane and undistorted after deformation.
The equation derived by Timoshenko that governs flexural vibrations of beams with constant cross section can be expressed as [36]  (1 ) The elastic geometrical properties of the beam elements have been derived using an energy relationship between molecular mechanics and continuum mechanics developed in [35]. According to the computation, the diameter d, Young's modulus E, and shear modulus G of the beam elements representing the C-C bonds are derived from where k r , k θ , k τ are the bond sketching, bond bending and torsional resistance force constants, respectively. In the finite element model of graphene sheets, the length of C-C bond is corresponding to the length of beam in planar-frame structure, wall thickness is related with the diameter of circular solid cross section of the beam elements, as in Figure 1.
The numerical simulation model of graphene sheets was created by the ANSYS parameter design language. For the modeling of the C-C bonds, the three-dimensional (3D) elastic BEAM188 element was used [20]. BEAM188 is suitable for analyzing slender to moderately stubby/thick beam structures. The element is based on the Timoshenko beam theory which includes first-order shear-deformation effects. The element is a linear, quadratic, and cubic two-node beam element in 3D. For each node, it has six degrees of freedom, which include translations in the x, y, and z directions and rotations around the x, y, and z axes. BEAM 188 is well-suited for linear, large rotation, and/or large strain nonlinear applications.

Finite Element Model in Vibration Analysis
The beam elements applied in graphene sheets are based on the Timoshenko beam theory, which is a first order shear deformation theory. In this theory, the transverse shear strain is constant through the cross section, which remains plane and undistorted after deformation.
The equation derived by Timoshenko that governs flexural vibrations of beams with constant cross section can be expressed as [36] EI ρA ∂ 4 ξ ∂z 4 − where ξ = ξ(z, t) is the transversal displacement along the x-axis at point z and time t, and E is the Young's modulus, I is the inertia moment, G is the shear modulus, ρ is the mass density, and A is the cross section area. In this theory, the Timoshenko shear coefficient κ is a free parameter. Besides ξ, an angular variable θ is introduced. During flexural motion, cross sections are supposed to remain flat and perpendicular to the deflected neutral axis at any point of this axis. The angle θ between the z-axis and a vector orthogonal to the cross section is equal to the angle between the neutral axis tangent line and the z-axis. Note that θ equals the slope of the deflected neutral axis, that is In a normal mode, ξ(z, t) varies harmonically with time as where A, B, C, and ϕ are the corresponding constants to be determined.w = 2π f is the angular frequency and χ(z) is a function that determines the normal mode amplitude. Substituting this form of ξ into Equation (1), It is well known that solutions of the above equation behave differently according to w 2 − w 2 c . The general solution can be written as where Usually, coefficients A i , B i , C i , D i are different from zero, and the solutions of equations include functions depending on both K 1 and K 2 . Where K 1 and K 2 are defined as positive square roots.
For free vibration analysis for Timoshenko beam, based on principal of virtual work, the weak form of equation can be written as [36] As defined above, ξ is the transversal displacement in Timoshenko beam, where θ is the transversal rotation, while .. θ are the transverse and rotary accelerations, respectively, L is the length of the beam and δ denotes that the terms are virtual.
For free vibration where K and M are the global stiffness and mass matrices which contain contributions from element stiffness and mass matrices.
A general solution can be expressed as Substituting into Equation (8) yields where φ k is a set of displacement-type amplitude at the control points otherwise known as the model vector, w k is the natural frequency associated with the kth mode. This is an eigenvalue problem and for non-zero solutions the determinant of the equation must be zero The discretization of governing equation can be done by unidirectional linear finite element with two nodes. In each axes direction, there are two degrees of freedom at each node, and six degrees of freedom for each node in total. The approximated solution in the displacement field, transversal displacement and rotation, can be written as [37]  Consider the natural coordinate ς = [−1, 1] in the element domain, In matrix form, it can be presented by (15) where [H] is the matrix of shape functions and {u} is vector of nodal displacements.
[H] = By substituting to the weak form of equation, the stiffness and mass matrix in element domain are written as (17) where [B] is the strain displacement matrix, [D] is the constitutive matrix, and |J| is the Jacobian determinant.
For large symmetric eigenvalue problems, the Block Lanczos eigenvalue extraction method is popular and available [38]. A fast convergence rate can be reached by this method. In modal analysis, Block Lanczos method is combined with Sturm sequence checks, the number of eigenvalues requested is extracted by an automated shift strategy. In addition, the Sturm sequence check also makes sure that the requested number of eigen-frequencies beyond the user supplied shift frequency is found without missing any modes [39].

Monte Carlo Based Finite Element Method
The Monte Carlo approach is a method developed to compute integrals by random number generation. Suppose a complex integral like Decompose h(x) into the production of a function f (x) and a probability density function p(x) defined in the interval (a, b), then the integral can be written as If a large number of random variables from the probability density p(x) are sampling, then Consider the integral I(y) = f (y|x )p(x)dx, it can be approximated by Monte Carlo simulation asÎ Therefore, the estimated Monte Carlo standard error can be expressed as Based on Equation (23), when n is amplified, the Monte Carlo standard error will become tiny, which means that the Monte Carlo simulation can reach a satisfied accuracy when the number of sampling is large enough. However, the huge sampling sets for the finite element model cost expensive computation. There is a trade-off between the results' accuracy and computational time costs. Depending on the work of CHU in 2015 [24], the number of Monte Carlo simulation is chosen as 500 in this study.
In order to implement the Monte Carlo based finite element method for vibration analysis of graphene sheets, the flowchart is present in Figure 2. The blue boxes represent the program of the finite element method, the red boxes indicate the Monte Carlo simulation. In order to combine the Monte Carlo simulation with the beam finite element method, the exact design of the computational simulation in the flowchart can be concluded in seven steps: Step 1: Define the initial configuration of graphene sheets, which includes corresponding parameters of bonds' length, height, and width of a hexagonal 2D lattice, and also thickness or diameter of the section.
Step 2: Fix material property parameters, which consist of the Young's modulus, Poisson ratio, physical density. The accuracy of the finite element model for graphene sheets heavily depends on exact parameters of material property.
Step 3: Mesh the beam elements and apply the boundary conditions, which is a traditional step of the finite element method.
Step 4: Perform the finite element method and use the Block Lanczos eigenvalue extraction for natural frequencies in each mode vibration of graphene sheets. If the deterministic finite element model is confirmed to be validated, the Monte Carlo simulation will commence, otherwise, return back to the step of initial configuration.
Step 5: Apply the Monte Carlo simulation method to randomly disperse vacancy defects in the graphene sheets in the certain percentage. The corresponding number of chosen bonds is evident and can be used to create the vacancy defects in the graphene.
Step 6: Mesh the beam elements and apply the boundary conditions, and perform the finite element model computation.
Step 7: Capture the natural frequencies of graphene with vacancy defects in the vibration mode and output the results.
Taking consideration of the effects of stochastic distributed vacancy, the loop needs to repeat to provide reliable and accurate output results. The loop continues until the Monte Carlo simulation is completed. In order to combine the Monte Carlo simulation with the beam finite element method, the exact design of the computational simulation in the flowchart can be concluded in seven steps: Step 1: Define the initial configuration of graphene sheets, which includes corresponding parameters of bonds' length, height, and width of a hexagonal 2D lattice, and also thickness or diameter of the section.
Step 2: Fix material property parameters, which consist of the Young's modulus, Poisson ratio, physical density. The accuracy of the finite element model for graphene sheets heavily depends on exact parameters of material property.
Step 3: Mesh the beam elements and apply the boundary conditions, which is a traditional step of the finite element method.
Step 4: Perform the finite element method and use the Block Lanczos eigenvalue extraction for natural frequencies in each mode vibration of graphene sheets. If the deterministic finite element model is confirmed to be validated, the Monte Carlo simulation will commence, otherwise, return back to the step of initial configuration.
Step 5: Apply the Monte Carlo simulation method to randomly disperse vacancy defects in the graphene sheets in the certain percentage. The corresponding number of chosen bonds is evident and can be used to create the vacancy defects in the graphene.
Step 6: Mesh the beam elements and apply the boundary conditions, and perform the finite element model computation.
Step 7: Capture the natural frequencies of graphene with vacancy defects in the vibration mode and output the results.
Taking consideration of the effects of stochastic distributed vacancy, the loop needs to repeat to provide reliable and accurate output results. The loop continues until the Monte Carlo simulation is completed.

Validation of the Model
By commercial software of ANSYS (Version 14.5, APDL, Cannonsburg, PA, USA), the typical graphene sheet lattice structure is created, which is in accord with the certain sizes in the literatures. In the truss finite element model, there are 6226 trusses, 4212 points and 16,664 nodes created. Performing the beam finite element model for graphene without vacancy defects, a good agreement between the present results and that of relevant references is achieved. Table 1 lists the corresponding parameters of material property of graphene in literatures, and in proposed model, the results of natural frequencies for the first four modes of vibration.
From Table 1 and Figure 3, it is obvious that the present results of the beam finite element model have satisfied agreement with the corresponding results in the reported references. To be more exact, the natural frequencies of graphene in this study are close to the results of Reddy, Lu, Wei, Kudin, Liu, and Cadelano; slightly higher than the results of the molecular dynamics used by Khatibi; and lower than the results of Zhou and Gupta. The present beam finite element model is feasible to be applied in the subsequent study of vacancy defected graphene sheets.

Validation of the Model
By commercial software of ANSYS (Version 14.5, APDL, Cannonsburg, PA, USA), the typical graphene sheet lattice structure is created, which is in accord with the certain sizes in the literatures. In the truss finite element model, there are 6226 trusses, 4212 points and 16,664 nodes created. Performing the beam finite element model for graphene without vacancy defects, a good agreement between the present results and that of relevant references is achieved. Table 1 lists the corresponding parameters of material property of graphene in literatures, and in proposed model, the results of natural frequencies for the first four modes of vibration.
From Table 1 and Figure 3, it is obvious that the present results of the beam finite element model have satisfied agreement with the corresponding results in the reported references. To be more exact, the natural frequencies of graphene in this study are close to the results of Reddy, Lu, Wei, Kudin, Liu, and Cadelano; slightly higher than the results of the molecular dynamics used by Khatibi; and lower than the results of Zhou and Gupta. The present beam finite element model is feasible to be applied in the subsequent study of vacancy defected graphene sheets.   The vacancy defects are quantified according to defect density, which is defined as where D n is the amount of vacancy defects and A n is the total quantity of beams in graphene.
The relative error is used to calculate the difference of the mean of the natural frequencies with that of graphene without defects. Error can be written as where M f is the mean of natural frequencies of vacancy defected graphene, P f is the natural frequency of pristine graphene.

Results and Discussion
In the following part, the explicit discussion about the amount of vacancy defects and parameters corresponding to geometrical and material characteristics of graphene sheets are sequentially presented. The results are computed by the Monte Carlo based finite element method and demonstrated by the probability theory and mathematical statistics.

Amount of Vacancy Defects
The amount of vacancy defects accounted in this study is by the absence of beam in the stochastic placement. Table 2 provides the statistical results of the first four order natural frequencies of vibration modes for graphene sheets with different number of vacancy defects. By taking the uncertainty of randomly distributed location of vacancy defects into consideration, the probability and statistical results of natural frequencies are more reliable and suitable for the study of vibration behavior of graphene sheets. It is clear that, with the increase of vacancy defects, the mean of the first four order natural frequencies are diminished. When Per is 0.5% and 1%, the changes of natural frequencies caused by randomly distributed vacancy defects are not evident. However, the reduction of natural frequencies are abrupt when the percentage of vacancy defects exceeds 1%.
The natural frequencies of the beam finite element model as derived in Equation (9) are corresponding to the stiffness and mass matrices. The fluctuation in natural frequencies in different vibration modes is a comprehensive effects of change in Young's modulus and mass reduction by vacancy defects. According to the relative reference [40], the elastic modulus of graphene can increase by well-controlled defect creation. The Young's modulus reaches at the peak when the percentage of vacancy defects of carbon atoms is 0.2%. It is reasonable to find that when Per equals to 0.5% and 1%, the changes of natural frequencies ( Error ) in Table 2 and Figure 5c are not evident.

Results and Discussion
In the following part, the explicit discussion about the amount of vacancy defects and parameters corresponding to geometrical and material characteristics of graphene sheets are sequentially presented. The results are computed by the Monte Carlo based finite element method and demonstrated by the probability theory and mathematical statistics.

Amount of Vacancy Defects
The amount of vacancy defects accounted in this study is by the absence of beam in the stochastic placement. Table 2 provides the statistical results of the first four order natural frequencies of vibration modes for graphene sheets with different number of vacancy defects. By taking the uncertainty of randomly distributed location of vacancy defects into consideration, the probability and statistical results of natural frequencies are more reliable and suitable for the study of vibration behavior of graphene sheets. It is clear that, with the increase of vacancy defects, the mean of the first four order natural frequencies are diminished. When Per is 0.5% and 1%, the changes of natural frequencies caused by randomly distributed vacancy defects are not evident. However, the reduction of natural frequencies are abrupt when the percentage of vacancy defects exceeds 1%.
The natural frequencies of the beam finite element model as derived in Equation (9) are corresponding to the stiffness and mass matrices. The fluctuation in natural frequencies in different vibration modes is a comprehensive effects of change in Young's modulus and mass reduction by vacancy defects. According to the relative reference [40], the elastic modulus of graphene can increase by well-controlled defect creation. The Young's modulus reaches at the peak when the percentage of vacancy defects of carbon atoms is 0.2%. It is reasonable to find that when Per equals to 0.5% and 1%, the changes of natural frequencies (Error) in Table 2 and Figure 5c are not evident.  When Per is 3%, the reduction of the first four order natural frequencies compared with that of pristine graphene are 7.1%, 8.7%, 6.04%, and 8.7%, respectively. Besides, the reduction reaches 32.1%, 46.7%, 41.3%, and 53.1% for first four order vibration mode when Per equals to 5%. As shown in Table 2 and Figure 5, with the increase of vacancy defects, the first four natural frequencies all have different degrees of reduction. Besides, the standard variance of natural frequencies becomes larger When Per is 3%, the reduction of the first four order natural frequencies compared with that of pristine graphene are 7.1%, 8.7%, 6.04%, and 8.7%, respectively. Besides, the reduction reaches 32.1%, 46.7%, 41.3%, and 53.1% for first four order vibration mode when Per equals to 5%. As shown in Table 2 and Figure 5, with the increase of vacancy defects, the first four natural frequencies all have different degrees of reduction. Besides, the standard variance of natural frequencies becomes larger along with the increase of Per. Furthermore, in the first four vibration modes, the standard variances increase sharply when Per is above 1%. The vacancy defects greatly affect the vibration behavior of graphene sheets. The existence of vacancy defects changes the structure of graphene lattice and has a profound impact on natural frequencies.
As mentioned above, natural frequencies are the ratios of stiffness and mass matrices. The appearance of vacancy defects influences both the stiffness and mass matrices. The mass of graphene sheets decreases according to the increase of the amount of vacancy defects. The influence of vacancy defects in Young's modulus of graphene sheets is not monotonous and non-linear [31,41]. When the number of vacancy defects is small, by effectively controlling the distribution of vacancy defects, the augmentation of the amount of vacancy defects in graphene sheets can increase the Young's modulus [40]. However, when the amount of vacancy defects reaches a certain value, the Young's modulus dramatically reduces. The augmentation of vacancy defects leads to not only cutting down in Young's modulus, tensile strength, and failure strain, but also fracture propagation and crack formation [42]. Vibration analysis is an appropriate and feasible way to observe the influence of the amount of vacancy defects in graphene sheets.

Geometrical Parameters
On the basis of molecular mechanics by the energy approach, the tube diameter and length strongly effect the material properties, such as Young's modulus, shear modulus, bend-angle and strain energy of single-and multi-walled nanotubes [43]. The discussion about geometrical parameters, namely bond length L, diameter of cross section D, and the amount of hexagons in y axis N, is important and necessary in vibration analysis of graphene sheets.
The results are based on the statistical results by consideration of the random dispersion of vacancy location in the graphene sheets. In Figure 6, with the amplification of L, the mean of first four order natural frequencies is distinctly cut down. The standard variance and error in Figure 6b,c also follow this rule. When the amount of the vacancy defects is 1% and L is 0.15, with the augmentation of length of each bond in hexagon of graphene sheets, the increase of the first four order natural frequencies compared with that of graphene sheets in Table 2 are 223.98%, 224.01%, 224.00%, and 223.90%, respectively. When L is equal to 0.2, the increases are 82.26%, 82.27%, 82.25%, and 82.22%, respectively. When L is 0.32, the reductions of the first four order natural frequencies of vibration modes are 28.82%, 28.79%, 28.82%, and 28.82%, respectively. Therefore, the length of bond in hexagon of graphene sheets is a crucial factor, which shows a good agreement of size effects in physical experiments and numerical simulation. According to the results above, the geometrical parameter L has a negative effect on first four order natural frequencies. The larger L is, the smaller the natural frequencies are. Different from the effects of L in vibration behavior of vacancy defected graphene sheets, with the increase of D, the first four order natural frequencies are amplified. Figure 7 well proves this point. The standard variance in Figure 7b also follows this rule, but the results of errors in Figure 7c are different. Compared with the results in Table 2 when Per equals to 1% and D is 0.02, the reduction of the first four order natural frequencies of vibration modes are 37.51%, 37.49%, 37.49%, and 37.50%, respectively. When D is equal to 0.026, the reduction of the first four order natural frequencies of vibration modes compared with that of graphene in Table 2 are 18.73%, 18.74%, 18.74%, and 18.75%, respectively. When L is 0.04, the increases of the first four order natural frequencies of vibration modes are 24.97%, 25.02%, 24.97%, and 24.99%, respectively. In Equations (2) and (8), both the cross section and inertia moment are related with the diameter of cross section D, which has more complicated effects in natural frequencies of vibration modes. Even though the change of extent in natural frequencies is not as large as the effect of the bond length L, more emphasis should be put on the diameter of cross section D. According to the results above, the geometrical parameter L has a negative effect on first four order natural frequencies. The larger L is, the smaller the natural frequencies are. Different from the effects of L in vibration behavior of vacancy defected graphene sheets, with the increase of D, the first four order natural frequencies are amplified. Figure 7 well proves this point. The standard variance in Figure 7b also follows this rule, but the results of errors in Figure 7c are different. Compared with the results in Table 2 when Per equals to 1% and D is 0.02, the reduction of the first four order natural frequencies of vibration modes are 37.51%, 37.49%, 37.49%, and 37.50%, respectively. When D is equal to 0.026, the reduction of the first four order natural frequencies of vibration modes compared with that of graphene in Table 2 are 18.73%, 18.74%, 18.74%, and 18.75%, respectively. When L is 0.04, the increases of the first four order natural frequencies of vibration modes are 24.97%, 25.02%, 24.97%, and 24.99%, respectively. In Equations (2) and (8), both the cross section and inertia moment are related with the diameter of cross section D, which has more complicated effects in natural frequencies of vibration modes. Even though the change of extent in natural frequencies is not as large as the effect of the bond length L, more emphasis should be put on the diameter of cross section D. There is another geometrical parameter, N, the amount of hexagons in y axes, which is not only related to the width of graphene sheets, but also corresponding with the amount of vacancy defects if Per is determined as a certain value. Compared with the results in Table 2 and Figure 8, when Per equals to 1% and N is 30, the natural frequencies of the first four orders of vibration mode are 47.90%, 20.83%, 66.55%, and 22.85% larger, respectively. When N is equal to 35, the increase of the natural frequencies of the first four vibration modes compared with that of graphene in Table 2 are 18.32%, 7.81%, 26.04%, and 15.68%, respectively. When N is 45, the reductions of the natural frequencies of the first four order vibration modes are 11.84%, 6.34%, 16.43%, and 11.76%, respectively. There is another geometrical parameter, N, the amount of hexagons in y axes, which is not only related to the width of graphene sheets, but also corresponding with the amount of vacancy defects if Per is determined as a certain value. Compared with the results in Table 2 and Figure 8, when Per equals to 1% and N is 30, the natural frequencies of the first four orders of vibration mode are 47.90%, 20.83%, 66.55%, and 22.85% larger, respectively. When N is equal to 35, the increase of the natural frequencies of the first four vibration modes compared with that of graphene in Table 2 are 18.32%, 7.81%, 26.04%, and 15.68%, respectively. When N is 45, the reductions of the natural frequencies of the first four order vibration modes are 11.84%, 6.34%, 16.43%, and 11.76%, respectively. In short, different with the effects of bond length L and the diameter of section cross D, the influence of N in natural frequencies of vacancy defected graphene sheets is more intricate and disproportionate, as in Figure 8b,c. In general, more emphasis is put on the bond length and the diameter of the cross section, the total width and length of graphene sheets are controlled and limited by experimental operation and chemical methods. However, the number of hexagons in the axes directions is another sensitive factor related to natural frequencies of vacancy defected graphene sheets.

Material Parameters
A large deviation between simulations and experiments in material properties has been found because of the presence of defects in the nanotube structure [2,44]. The mechanical properties of defect-free carbon nanomaterials have not been exactly measured yet. For discussion, material parameters (Young's modulus and Poisson ratio of graphene sheets) are supposed as specific values.
From Figure 9, it is distinct to find that, with the increase of E, namely Young's modulus, the mean of the first four natural frequencies evidently increased as well. The standard variance and error in Figure 9b also follow this rule. However, the results in Figure 8c are more complicated. When the amount of the vacancy defects is fixed to 1%, and E is 800 GPa, the reduction of the natural frequencies of the first four vibration modes compared with that of graphene sheets in Table 2 are In short, different with the effects of bond length L and the diameter of section cross D, the influence of N in natural frequencies of vacancy defected graphene sheets is more intricate and disproportionate, as in Figure 8b,c. In general, more emphasis is put on the bond length and the diameter of the cross section, the total width and length of graphene sheets are controlled and limited by experimental operation and chemical methods. However, the number of hexagons in the axes directions is another sensitive factor related to natural frequencies of vacancy defected graphene sheets.

Material Parameters
A large deviation between simulations and experiments in material properties has been found because of the presence of defects in the nanotube structure [2,44]. The mechanical properties of defect-free carbon nanomaterials have not been exactly measured yet. For discussion, material parameters (Young's modulus and Poisson ratio of graphene sheets) are supposed as specific values.
From Figure 9, it is distinct to find that, with the increase of E, namely Young's modulus, the mean of the first four natural frequencies evidently increased as well. The standard variance and error in Figure 9b also follow this rule. However, the results in Figure 8c are more complicated. When the amount of the vacancy defects is fixed to 1%, and E is 800 GPa, the reduction of the natural frequencies of the first four vibration modes compared with that of graphene sheets in Table 2 are 18.38%, 18.35%, 18.34%, and 18.36%, respectively. When E is equal to 1500 GPa, the increase of the natural frequencies of the first four vibration modes compared with that of pristine graphene are 11.79%, 11.82%, 11.80%, and 11.80%, respectively. When E is 2000 GPa, the increments of the natural frequencies of the first four vibration modes are 29.11%, 29.13%, 29.11%, and 29.09%, respectively.
With the augment of Young's modulus, the values in stiffness matrix of vacancy defected graphene sheets are amplified, while the stiffness matrix is acted as the numerator, which is reasonable to find the enlargement of natural frequencies in vibration modes. Furthermore, as mentioned in the above results, the first four natural frequencies of graphene sheets all have the close magnitude of change. The effects of Young's modulus in natural frequencies are convergent to the specific percentage in the different vibration modes.  18.38%, 18.35%, 18.34%, and 18.36%, respectively. When E is equal to 1500 GPa, the increase of the natural frequencies of the first four vibration modes compared with that of pristine graphene are 11.79%, 11.82%, 11.80%, and 11.80%, respectively. When E is 2000 GPa, the increments of the natural frequencies of the first four vibration modes are 29.11%, 29.13%, 29.11%, and 29.09%, respectively. With the augment of Young's modulus, the values in stiffness matrix of vacancy defected graphene sheets are amplified, while the stiffness matrix is acted as the numerator, which is reasonable to find the enlargement of natural frequencies in vibration modes. Furthermore, as mentioned in the above results, the first four natural frequencies of graphene sheets all have the close magnitude of change. The effects of Young's modulus in natural frequencies are convergent to the specific percentage in the different vibration modes. Poisson ratio as one of the parameters corresponding with material properties is discussed in Figure 10. The change of Poisson ratio in vacancy defected graphene sheets does not play significant role in the natural frequencies of graphene sheets. As in Figure 10a, with the increase of V, the results of natural frequencies are very approaching. In other words, the Poisson ratio can be simplified as a neglected factor in the analysis of natural frequencies for vacancy defected graphene sheets. Poisson ratio as one of the parameters corresponding with material properties is discussed in Figure 10. The change of Poisson ratio in vacancy defected graphene sheets does not play significant role in the natural frequencies of graphene sheets. As in Figure 10a, with the increase of V, the results of natural frequencies are very approaching. In other words, the Poisson ratio can be simplified as a neglected factor in the analysis of natural frequencies for vacancy defected graphene sheets.

Graphene Sheets with Vacancy Defects
In order to observe the vibration behavior of graphene with vacancy defects, Figures 11 and 12 present the displacement vector sum and rotation vector sum for graphene with 5% vacancy defects, respectively. Along with the augment of vacancy defects in the graphene sheet, the natural frequencies of vibration mode are evidently cut down. As demonstrated in Figures 11 and 12, the regular and symmetrical property is also destroyed because of the randomly distributed vacancy defects. Therefore, besides the change of the first four order natural frequencies in the vibration behavior of graphene sheets, the vacancy defects can also deeply influence the displacement and rotation vector of graphene sheets.

Graphene Sheets with Vacancy Defects
In order to observe the vibration behavior of graphene with vacancy defects, Figures 11 and 12 present the displacement vector sum and rotation vector sum for graphene with 5% vacancy defects, respectively. Along with the augment of vacancy defects in the graphene sheet, the natural frequencies of vibration mode are evidently cut down. As demonstrated in Figures 11 and 12, the regular and symmetrical property is also destroyed because of the randomly distributed vacancy defects. Therefore, besides the change of the first four order natural frequencies in the vibration behavior of graphene sheets, the vacancy defects can also deeply influence the displacement and rotation vector of graphene sheets. (c) (d) Figure 11. Displacement vector sum for graphene with 5% vacancy defects (a-d show the first order, second order, third order, and fourth order of vibration mode, respectively). Figure 11. Displacement vector sum for graphene with 5% vacancy defects (a-d show the first order, second order, third order, and fourth order of vibration mode, respectively).
(a) (b) (c) (d) Figure 12. Rotation vector sum for graphene with 5% vacancy defects (a-d represent the first order, second order, third order, and fourth order of vibration mode, respectively).
In addition, besides the effects of vacancy defects, impurities in graphene sheets are another unavoidable existence and need more attention [45]. Natural frequencies of vacancy defected graphene sheets are a sensitive and appropriate index to comprehensively discuss and study the influence of vacancy defects in stiffness and mass matrices of entire graphene sheets. In addition to mechanical properties, thermal and electronic transport properties are promising in perspective for vacancy defected graphene sheets [46]. In addition, besides the effects of vacancy defects, impurities in graphene sheets are another unavoidable existence and need more attention [45]. Natural frequencies of vacancy defected graphene sheets are a sensitive and appropriate index to comprehensively discuss and study the influence of vacancy defects in stiffness and mass matrices of entire graphene sheets. In addition to mechanical properties, thermal and electronic transport properties are promising in perspective for vacancy defected graphene sheets [46].

Conclusions
In this paper, an effective Monte Carlo based finite element method to analyze the vibration behavior of vacancy defected graphene sheets is proposed. By combining the Monte Carlo method with the beam FE model, the stochastic dispersed placement of vacancies in graphene is successfully propagated and performed. The results show that, 1.
The natural frequencies of vibration decrease dramatically with the increase of the vacancy defects when Per is larger than 1%. 2.
The regular and symmetrical properties of vibration behavior in displacement and rotation are evidently destroyed when Per reaches 5%.

3.
Furthermore, with the augment of the vacancy defects quantity, the variances caused by randomly dispersion location of vacancy defects become large.

4.
From the discussion about geometrical and material parameters, it can be concluded that the increment of length of each bond in hexagons of graphene sheets distinctly reduces the first four order natural frequencies.

5.
For vacancy defected graphene sheets, the amount of hexagons in the long side of the graphene sheet has a similar effect on natural frequencies with length of each bond in hexagon, but it is more complicated. 6.
With the increase of thickness of the graphene sheets and Young's modulus, the natural frequencies of vacancy defected graphene become smaller.