Monte Carlo-Based Finite Element Method for the Study of Randomly Distributed Vacancy Defects in Graphene Sheets

This paper proposed an effective stochastic finite element method for the study of randomly distributed vacancy defects in graphene sheets. The honeycomb lattice of graphene is represented by beam finite elements. The simulation results of the pristine graphene are in accordance with literatures. The randomly dispersed vacancies are propagated and performed in graphene by integrating Monte Carlo simulation (MCS) with the beam finite element model (FEM). The results present that the natural frequencies of different vibration modes decrease with the augment of the vacancy defect amount. When the vacancy defect reaches 5%, the regularity and geometrical symmetry of displacement and rotation in vibration behavior are obviously damaged. In addition, with the raise of vacancy defects, the random dispersion position of vacancy defects increases the variance in natural frequencies. The probability density distributions of natural frequencies are close to the Gaussian and Weibull distributions.


Introduction
Graphene, discovered in 2004 [1], is a two-dimensional (2D) carbon nanomaterial composed of a hexagonal honeycomb lattice.Based on its superior properties, such as good flexibility and high thermal and electrical conductivity [2,3], graphene has a huge potential in the field of aerospace, microelectronics, microelectromechanical systems, and nanocomposites [4,5].For applications and development in different fields, the analysis of stress states, vibration behavior, and stability characteristics attract much attention of researchers and designers.
Among the material property corresponding parameters, the intrinsic strength and Young's modulus of graphene monolayers vary [2] in a large interval.It has been found that, owing to the appearance of defects in the microstructure of nanomaterials, the evident deviation and nonnegligible fluctuation occur in numerical simulations and physical experiments [25][26][27].Therefore, attempts and struggles are required for studying and simulating the defects in graphene to explain the fluctuation in simulations or experiments reasonably.
The Monte Carlo simulation (MCS) is a reliable sampling method with wide applications in engineering and research [28][29][30][31].With sufficiently large sampling space, a satisfactory accuracy can be achieved in the numerical simulation.The results of MCS are usually considered as an exact solution or comparison standard [32,33].This study utilizes the Monte Carlo-based finite element model (MC-FEM) to effectively discuss the influences of vacancy defects in the vibration behavior of graphene.MC-FEM has the advantages of convenient programming, high computational accuracy, and good convergence.In addition, MC-FEM can be used to overcome the difficulty of propagation of the unpredictable and stochastic vacancy defects in graphene.This study provides a fresh idea and perspective, which are meaningful and important to inspire researchers in nanomaterials and related fields.
The remaining parts of this article are structured as follows: in Section 2, there are a clear introduction about the lattice structure model of graphene, the mathematical foundation of the MCS, and the program design of the proposed MC-FEM.In Section 3, the beam finite element model (FEM) of pristine graphene sheets is verified by comparing with the results in literatures [34][35][36][37][38][39][40][41][42].Section 4 provides the discussion about the influence of vacancy defects in the vibration behavior of graphene.The statistical variances and probability density distribution produced by random dispersion of vacancy defects are analyzed and discussed.The last section gives a brief summary of this paper.

Material and Method
2.1.Graphene Sheets.Carbon atoms in graphene combine together with covalent bonds to form a honeycomb 2D lattice.The displacement of a single atom is constrained by a bond under external forces.As a result, the total deformation and displacement of graphene are the outcomes of the interactions between bonds.The elastic BEAM188 element in ANSYS is chosen to model C-C bonds as in the reported reference [43].BEAM188 is based on the Timoshenko beam theory with the first order shear deformation effect and is appropriate for analysis slender to medium short/thick beam structures.
A circular solid cross-sectional area is assumed in the finite element model.The wall thickness corresponds to the section diameter.It is desired to compute the natural frequencies of graphene in different vibration modes.Figure 1 illustrates a schematic of monolayer graphene in two dimensions.The Zigzag and Armchair types are both included.

Monte Carlo Sampling
Method.Suppose ξ represents a random variable with the mathematical expectation of E ξ = I.The definition can be written as where the probability density function p x can be continuous or discrete.In order to approximate the variable I, the value of arithmetic mean is computed by The convergence of the arithmetic mean of the input random variables is predicted as ξ N → p I as N → ∞ The sequence of the random variables η 1 , η 2 , … , η N , … converges to the constant c, if for every h > 0, it follows that Thus, when N is large enough, ξ N ≈ I.
If the input random variable ξ exists a finite variance, the relative error of MCS can be expected as For the problem of the integral computation by MCS, where Ω represents an arbitrary domain and x ∈ Ω ⊂ R d is a vector in d-dimension and the density function p x is nonnegative with Ω p x dx = 1.ξ is a random point with probability density function p x .Introducing the random variable, With mathematical expectation equal to the value of integral I,

Journal of Nanomaterials
The input random points ξ 1 , ξ 2 , … , ξ N are independent, and according to the probability density function p x , an approximate value of I can be computed by Therefore, θ N will be convergent to I based on ξ N = 1/N∑ N i=1 ξ i .
2.3.Monte Carlo-Based Finite Element Method.With the purpose of combining the MCS with the beam FEM, the flowchart of MC-FEM is described in Figure 2, which is summed up in the 7 steps: (1) The original configuration of graphene lattice is defined, including the geometrical characteristic corresponding parameters of the bonds in the honeycomb microstructure, as well as the thickness or diameter of the cross section in the beam FEM.(2) Material parameters are fixed, including Young's modulus, physical density, and Poisson ratio.The precision of material parameters strongly influences the results of the finite element model of graphene sheets.
(3) The MCS is used to effectively propagate the vacancy defects randomly in graphene within a specific percentage.
(4) The amount of selected bonds representing the vacancy defects is generated by MCS.According to the stochastic numbers, the related bonds in graphene lattice are removed to create the vacancy defects.
(5) Meshing beam elements and definition of boundary conditions are the traditional steps of the finite element method and are accomplished in this step.
(6) The finite element computation is performed, and the natural frequencies of graphene sheets in each vibration mode are extracted and recorded.(7) The natural frequencies of the defected graphene are obtained, and the results are output.
(8) Considering the effect of the random distribution of vacancy defects, the Monte Carlo simulation repeats until sufficient times are completed to provide reliable results.

Validation of the Model
The original beam FEM of pristine graphene consists of 6226 beams, 18,678 elements, and 16,664 nodes.The parameters related to graphene material properties in the literature and the present model are demonstrated in Table 1.The results of the natural frequencies in the first four-order vibration modes are verified by the comparison with the reported results, as presented in Figure 3.
The beam finite element model of pristine graphene is created and performed; the results of natural frequencies are confirmed to be consistent with those of existing literatures.Specifically, the results of natural frequencies have good agreement with those of Lu, Wei, Reddy, Liu, Kudin, and Cadelano.Also, the results are relatively greater than those of Khatibi computed by the molecular dynamics, but smaller than those of Zhou and Gupta.Therefore, it is feasible to use the beam finite element model to represent the hexagonal structure of graphene and to simulate its modes of vibration.
Since the differences between the present results and the reported references are in an acceptable range, the specific vibration behavior of graphene is demonstrated in Figures 4  and 5. Figures 4 and 5 plot the contour results of the displacement vector sums and rotation vector sums for Zigzag graphene and Armchair graphene, respectively.For pristine graphene, the results of displacement vector sums and rotation vector sums are regular and symmetrically distributed according to the geometrical characteristics in both Zigzag and Armchair types.They are settled as criterions and will be compared with the following results of graphene sheets with vacancy defects.where D n is the amount of removed beams in vacancydefected graphene while A n is the total number of beams in pristine graphene.Figures 6 and 7 express the displacement vector sums and rotation sums of graphene with 1% vacancy defects in the Zigzag and Armchair types, respectively.By comparing with results of pristine graphene, it is easy to find that the existence of vacancy defects in graphene sheets makes the natural frequencies in different modes of vibration decrease.However, the effects of vacancy defects on contour results of displacement and rotation are not very evident either for the Zigzag or Armchair types.The contour results of displacement and rotation are still observed as regular and symmetrical according to the geometrical characteristics.6

Journal of Nanomaterials
According to the increase of the vacancy defect amount in graphene, the natural frequencies in the first four-order vibration modes are gradually cut down.The vacancy defects deeply impact the vibration behavior of graphene.Especially for the fourth mode of vibration, as demonstrated in Figures 8 and 9, the displacement and rotation vector sums are largely different from those of pristine graphene.The regular and symmetrical properties are damaged owing to the arbitrarily distributed vacancy defects in both the Zigzag and Armchair graphene sheets.
In addition, by comparison between the above results of Zigzag and Armchair graphene sheets, the differences caused by graphene chirality in natural frequencies of the vibration modes are tiny and negligible.However, as a result of distinct microstructures in Zigzag and Armchair graphene sheets, the divergences happened in the results of displacement and   rotation, as presented in contour results of MC-FEM.To be more exact, the vibration behavior of Zigzag graphene sheets is close to 90 degrees of rotation of that in Armchair graphene sheets.This phenomenon is not only caused by the characteristics of the microstructures in Zigzag and Armchair graphene sheets but also due to the geometrical parameter setting in this study, and the length and width of graphene sheets are approximate to each other.

Statistical Results of MC-FEM.
The following study is focused on the natural frequencies of graphene sheets in the first four-order vibration modes.The difference of natural frequencies in Zigzag and Armchair graphene sheets is insignificant; therefore, Zigzag graphene sheets are chosen as the typical representation.Given that the vacancy defects are randomly dispersed in graphene sheets, one numerical simulation of the MC-FEM simulates one possibly existing situation.However, the uncertainty is unavoidable for the location of vacancy defects in graphene.It is crucial to repeat the MCS and perform the for sufficient times to prove the reliability of results.Table 2 presents the statistical results of natural frequencies in different modes of vibration by repeating the Monte Carlo sampling and performing the beam finite element method.The mean values of natural frequencies are more comprehensive and believable for the analysis of graphene vibration behavior.Obviously, with the raise of Per, the mean values of natural frequencies in the first four-order vibration modes identically diminish.When Per exceeds 1%, the decrease of natural frequencies is sudden.This is also proven in Figures 10 and 11.In Figure 11, the comparison between results of pristine graphene and graphene with vacancy defects can be defined as where M f is the mean value of natural frequencies in different vibration modes for the vacancy-defected graphene and P f is the natural frequencies in different vibration modes of pristine graphene.When Per is equal to 0.2%, 0.5%, or 1%, the relative errors are smaller than 0.05 THz, which are tiny and negligible, in the first four-order vibration modes.From the first to the fourth vibration mode, the relative error gradually increases with the raise of vacancy defects, as depicted in Figure 11.For Per is 3%, the reduction of natural frequencies compared with that of pristine graphene is less than 10%.When Per is equal to 5%, the reductions of natural frequencies are up to 53.1% in the fourth order vibration mode.It can be seen from the above results that the vacancy defects not only greatly change the natural frequencies and vibration behavior of graphene sheets but also affect the structure of graphene lattices.Specifically, the natural frequencies of the first four orders all decrease to some extent with the augment of vacancy defect amount.
Furthermore, statistical results are obtained by considering the uncertain dispersion of vacancy positions in graphene.With the increase of Per, the standard variances of natural frequency are amplified.As shown in Table 2 and Figure 12, the standard variances increase significantly in the first four-order vibration modes, especially if the vacancy defect percentage exceeds 1%. 9 Journal of Nanomaterials to 5%, the reduction of natural frequencies exceeds 40%.The characteristics of graphene lattice are strongly corresponding with the vibration behavior.When the amount of vacancy defects increases, the symmetrical and regular properties of vibration behavior in rotation vector sums and displacement vector sums are obviously damaged.
The randomly distributed placement of vacancy defects causes fluctuation and differences in the natural frequencies of vibration modes.The augment of the Per in the graphene sheets causes the reduction of natural frequencies in different vibration modes and also makes the variances and relative errors of natural frequencies become enlarged as in Figures 11 and 12. Therefore, the vacancy defects not only deeply influence the vibration behavior of graphene sheets but also amplify the uncertainty in each vibration mode.
Figures 13 and 14 are provided to analyze the probability density distribution of the natural frequencies for graphene with 0.2% and 1% vacancy defects, respectively.Both the mean and standard variance of graphene sheets with 0.2% vacancy defects demonstrate that its vibration behavior is close to pristine graphene.However, the probability density distribution in Figure 13 proves that the stochastic placement of vacancy defects in graphene causes fluctuation in the results of natural frequencies.As the same with that in Figure 13, when Per is 1%, the probability density distribution of the natural frequencies approaches to the Gaussian and Weibull distributions.When compared with other common probability density distributions, such as the gamma, exponential, and Rayleigh distribution, the Gaussian and Weibull distributions are much better to represent the distribution of the natural frequencies in all the first fourorder vibration modes.
In short, the analysis of the probability density distribution is a proper way to find the impacts of randomly distributed vacancy defects in the vibration behavior of graphene.The Gaussian and Weibull density distributions are the 10 Journal of Nanomaterials more approximated probability models to describe the effects of vacancy defects on uncertainty of natural frequencies.Besides, the results of the probability density distribution of vacancy-defected graphene sheets are meaningful for reliability and risk analysis.

Conclusion
In this paper, detailed studies on the vibration behavior of vacancy-defected graphene are carried out and performed by MC-FEM.The random dispersion of vacancy defects in graphene is taken into consideration, and the effects of the amount and randomly distributed placement of vacancy defects are discussed.The beam FEM has a good accuracy for computing natural frequencies of different vibration modes.With the augment of vacancy defect amount, the regular and symmetrical characteristics in the results of displacement and rotation are obviously damaged.Large amount of vacancy defects also amplifies the uncertainties in vibration behavior of graphene.The Weibull and Gaussian distributions are more appropriate to represent the effects of randomly dispersed placement of vacancy defects.Therefore, MC-FEM is a powerful model to simulate the mechanical behavior of vacancy-defected graphene sheets.

Figure 6 :
Figure 6: Results of MC-FEM for Zigzag graphene sheets with 1% vacancy defects ((a), (b), (c), and (d) are the displacement vector sums in the first four-order vibration modes, resp., and (e), (f), (g), and (h) are the rotation vector sums in the first four-order vibration modes, resp.).

Figure 7 :
Figure 7: Results of MC-FEM for Armchair graphene sheets with 1% vacancy defects ((a), (b), (c), and (d) are the displacement vector sums in the first four-order vibration modes, resp., and (e), (f), (g), and (h) are the rotation vector sums in the first four-order vibration modes, resp.).

Figure 8 :
Figure 8: Results of MC-FEM for Zigzag graphene sheets with 5% vacancy defects ((a-d) are the displacement vector sums in the first four-order vibration modes, respectively, and (e-h) are the rotation vector sums in the first four-order vibration modes, resp.).

Figure 9 :
Figure 9: Results of MC-FEM for Armchair graphene sheets with 5% vacancy defects ((a-d) are the displacement vector sums in the first four-order vibration modes, resp., and (e-h) are the rotation vector sums in the first four-order vibration modes, resp.).

Figure 10 :
Figure 10: Mean values of natural frequencies in different vibration modes.

Figure 11 :
Figure 11: Relative error of natural frequencies in different vibration modes.

Figure 12 :
Figure 12: Standard variance of natural frequency.

Figure 14 :
Figure 14: Probability density distribution of natural frequencies of graphene with 1% vacancy defects ((a), (b), (c), and (d) are the distributions of the first four-order natural frequencies, resp.).

Table 2 :
Statistical results of natural frequencies by MC-FEM.