The correlation between graphene characteristic parameters and resonant frequencies by Monte Carlo based stochastic finite element model

The uncertainty and fluctuations in graphene characteristic parameters are inevitable issues in both of experimental measurements and numerical investigations. In this paper, the correlations between characteristic parameters (Young’s modulus, Poisson’s ratio and thickness of graphene) and resonant frequencies are analyzed by the Monte Carlo based stochastic finite element model. Based on the Monte Carlo stochastic sampling procedure, the uncertainty in the characteristic parameters are properly propagated and quantified. The displacements and rotation modes of graphene under the resonant vibration computed by the finite element method are verified. Furthermore, the result robustness of stochastic samples is discussed based on the statistic records and probability density distributions. In addition, both the Pearson and Spearman correlation coefficients of the corresponding characteristic parameters are calculated and compared. The work in this paper provides a feasible and highly efficient method for the characteristic parameter correlation discussion by taking uncertainty into consideration.

the atomic force microscope (AFM) based on nanoindentation method. In another study, the Young's modulus equals to 1.0 ± 0.1 TPa and the corresponding intrinsic stress is 130 ± 10 GPa at a strain of 0.25 10 . Therefore, it is necessary to take the uncertainty and fluctuations in graphene characteristic parameters into consideration and perform quantification research.
For the uncertainty quantification and propagation, the traditional FEM of graphene has limitations due to the deterministic values of characteristic parameters. Based on MC-SFEM, it is possible to extend the deterministic values into the interval ranges for the characteristic parameters of graphene. On the other hand, the complicated relationships between characteristic parameters and specific mechanical responses of graphene cannot be described by explicit expressions  . The correlation coefficients are the effective and feasible factors to analyze the sensitive impacts of characteristic parameters on mechanical responses. In terms of performing correlation coefficient computation, the MC-SFEM is an effective method to provide the sufficient stochastic samples in the broad interval ranges.
In this paper, the MC-SFEM is performed by the integration of stochastic sampling simulation with FEM computation for uncertainty quantification and propagation. The corresponding resonant frequencies of graphene are tracked and recorded. The basic deterministic FEM is verified by the comparison with the previous reported work. The robustness of results for the different stochastic sample space sizes is analyzed and discussed. The correlation coefficients of the Young's modulus, Poisson's ratio, thickness of graphene are calculated, respectively. A short conclusion is summarized in the last section.

Model formation
Characteristic parameters. According to the periodic honeycomb structure of graphene, the specific hexagonal cell in FEM is defined by geometrical and material parameters, including the length of carbon-carbon bonds, the thickness of the graphene (d), the Young's modulus (E), the Possion's ratio (v) and the physical density. The geometrical structure of graphene in FEM is presented in Fig. 1. The length of carbon-carbon covalent bonds represents the distance between two carbon atoms, which is also the side length of periodic hexagon cell. The interaction between two carbon atoms is represented by the beam connection. In addition, the thickness of the graphene is introduced as the diameter of the beam cross section as shown in Fig. 1b. Among the geometrical and material parameters of graphene, d, E, v are chosen as the characteristic parameters for the impact discussion and correlation coefficient computation, while the length of carbon-carbon bonds and the physical density are settled as the deterministic values as static references in the relative comparison.
In order to set precise interval ranges for the characteristic parameters, the related references are recorded in Table 1 and Fig. 2. Based on the statistic data in the references, the interval ranges are 0.05 to 3 TPa, 0.01 to 0.48, and 0.05 to 0.7 nm, for the Young's modulus, Poission's ratio, and the thickness of graphene, respectively.
Correlation coefficient computation. The estimated correlation matrix A is a symmetric matrix of the order N var and can be written as   www.nature.com/scientificreports/ where I is the identity matrix and L is the strictly lower triangular matrix with entries within the range �−1, 1� . There are N c pairwise correlations: Pearson correlation coefficient. The most well-known correlation measurement is the linear Pearson correlation coefficient (PCC) 74 . The PCC takes values between − 1 and + 1, and provides an effective factor to evaluate the strength of the linear relationship between two variables. The actual PCC between two variables, X i and X j , can be estimated using the sample correlation coefficient A ij : where x i,s is the actual data matrix, s = 1, 2, . . . , N sim is the data sequence number of each vector, and N sim is the data sequence length of each vector. i = 1, 2, . . . , N var is the variable index. In this study, N var equals to 3, with Young's modulus ( i = 1 ), Poisson's ratio ( i = 2 ) and thickness of graphene ( i = 3).s is directly related to the number of stochastic samples. N sim is settled as 10,000 in this work to have sufficiently huge sample spaces. j = 1, 2, . . . , N ord is corresponding with the mode of resonant vibration. In addition. X i is the mean of Young's modulus, Poisson's ratio and thickness of graphene; while X j is computed by the mean of the jth order resonant frequencies of graphene. When the actual data x i,s of each vector i = 1, 2, . . . , N var are standardized into z i,s that yield zero averages and unit sample variance estimates, the formula can be simplified to The mean of variables can also be written as X i = µ X i = E(X i ) and the variances are computed as σ X 2 . Therefore, the sample correlation coefficient can be expressed as Spearman correlation coefficient. The formula for Spearman correlation coefficient 75 estimation is identical to the one for Pearson linear correlation except that the values of random variables X i and X j are replaced by the ranks π i,s and π j,s , s = 1, 2, . . . , N sim . The ranks are permutation of numbers. It is convenient to transform the ranks into The rank correlation is then defined as, By noting that the sum of the first N sim squared integers is N sim (N sim +1)(2N sim +1) 6 , we find that , and the rank correlation is, Another formula exists for Spearman correlation suitable for data with no ties. The correlation coefficient between any two vectors that each consisting of permutations of integer ranks from 1 to N sim is, where D is the sum of square of d s = s − π s , which is the differences between the s th integer elements in the vectors, x j,s , www.nature.com/scientificreports/ Every mutual permutation of ranks can be achieved by permuting the ranks π s of the second variable against the identity permutation corresponding to that of the first variable. Therefore,

Meanwhile, it can be obtained that
The lowest correlation is achieved for the reverse ordering of rank numbers and corresponds to the case when the sum D equals to . In turn, the maximum correlation is achieved for identical ranks and the sum equals zero.
Monte Carlo based stochastic finite element model. For the program implementation of MC-SFEM, two modules are compiled as presented in Fig. 3 with two different colors to remark. The process of finite element computation for the resonant frequencies of graphene is presented in the blue blocks in the flowchart, while the stochastic sampling procedure is presented in the red blocks. Finally, the input and output data are transferred into the correlation coefficient computation. To be clear, the integration of the Monte Carlo stochastic sampling procedure with the finite element computation is concluded in seven steps: Step 1: Define the initial configuration of graphene with the periodic hexagonal lattice, such as the thickness of graphene (diameter of the section in beam elements), which is stochastically sampled as a variable in the following process. The distance between the two carbon atoms and the number of hexagons in the longitude and latitude sides of graphene are defined with the deterministic values.
Step 2: Introduce material characteristic parameters, among which the Young's modulus and Poisson's ratio are also the variables in the Monte Carlo stochastic sampling procedure.
Step 3: Mesh the beam elements and define the boundary conditions, where the four edges of graphene are fixed, the six freedom degrees (displacements and rotations in x, y, z directions) of each node in the four edges of graphene are settled to be zeros.
Step 4: Perform the finite element computation and use the Block Lanczos eigenvalue extraction for the resonant frequency calculation. The validation of the deterministic FEM is the premise of the following MC-SFEM. www.nature.com/scientificreports/ Step 5: Apply the Monte Carlo stochastic sampling simulation to propagate the parameter uncertainty within specific intervals by the uniform probability distribution. Besides, the values of the characteristic parameters are recorded and stored as the input data for correlation analysis.
Step 6: Repeat the FEM computation and transfer the resonant frequencies of graphene as the output data.
Step 7: Based on the huge stochastic sample spaces of input and output data, the Pearson correlation coefficient and Spearman correlation coefficient are calculated.
It is important to be mentioned that for the Monte Carlo stochastic sampling method in Step 5, the Latin Hypercube sampling (LHS) method is performed to ensure the sampling efficiency and reduce computational costs. By dividing the range of each variable into disjoint intervals with equal probabilities, the samples are randomly selected from each interval. The LHS method effectively avoids the point clustering and duplication 76 , especially for the uniform probability density distribution in this study. Therefore, LHS procedure is used for random sample generation in the corresponding intervals of the characteristic parameters.
In addition, other probabilistic methods applied in the homogenization method are available alternatives as shown in literatures [77][78][79][80] . However, the uniform probability distribution is selected and performed in this study for several reasons. On the one hand, uncertainty and fluctuations in the characteristic parameters present broad interval ranges without deterministic probability density distributions. The uniform probability distribution is more comprehensive to fairly take each possible value in the certain intervals into consideration. On the other hand, the uniform probability distribution is more effective to provide stochastic samples for the correlation coefficient computation, since each possible value is unbiased sampled. Therefore, compared to the commonly used Gaussian, Weibull, Poisson, etc. probabilistic distributions, the uniform probability density distribution is more feasible to propagate the uncertainty for the characteristic parameters of graphene.

Model validation
Based on the program flowchart, the deterministic FEM of graphene is created by the ANSYS Parameter Design Language (Version 14.5, APDL, ANSYS, USA). The BEAM188 element is chosen to mesh the carbon-carbon bonds in graphene. Beam 188 element is suitable for analyzing slender to moderately stubby/thick beam structures. Each node of Beam 188 element has six degrees of freedom. In addition, Beam 188 element is based on the Timoshenko beam theory which is effective in resonant frequency computatioin.
There are totally 6226 beams with 16,664 nodes in the deterministic FEM of graphene. As presented in Fig. 4, the displacement vector sum and roation vector sum in different vibration modes have satisfied accuracy and convergence in FEM computaion. Both the resonant frequencies and the vibration modes reach a good aggreement with the previous reported references 2,4 . Therefore, the deterministic FEM is verified as the basic model for the following MC-SFEM.
In addition, the geometrical symmetry in the displacement and rotation contour results is evident in Fig. 4. However, the resonant frequencies of graphene in different vibration orders are discrepant. The periodic honeycomb lattice of graphene with symmetrical hexagons causes the geometrical symmetry in resonant vibration modes, but the resonant frequencies are distinct to observe the resonant vibration orders. Therefore, the resonant frequencies of graphene are convenient index to track discrepances in the different resonat vibration orders, while the vibration modes such as displacement and rotation contours are senstive to the microstruture of graphene.

Discussion and results
Result robustness. In order to record the result accuracy and convergency, the mean, maximum, minimum and standard variance values with the different stochastic sample sizes based on the MC-SFEM are recorded and compared in Fig. 5. The 10,000 groups of stochastic samples of FEM computation are supposed to be the precise results, and the compared sample space sizes are 1000 and 5000. It is evident to find that the discrepancies in the mean, maximum or standard variance values with different sample space sizes are not obvious. However, the deviations caused by the small stochastic space is observed in the minimum value of the 1000 sample groups, the statistic results of the 5000 sample groups are very approximated to that of the 10,000 sample groups. Even if the MC-SFEM has satisfied computational efficiency and accuracy, a sufficient number of stochastic samples is also a key factor. The larger sample space size requires heavier computation burden. Therefore, the trade-off between the result accuracy and computation cost exists. The discrepancies in the mean, maximum and standard variance are not as evident as that in the minimum for the size of stochastic sampling space.
Furthermore, the probability density distribution of the first ten resonant frequencies of graphene is also calculated according to the stochastic sampling space of MC-SFEM. The size of stochastic sampling space is 5000. Even if the characteristic parameters are distributed and sampled by the uniform distribution type, the probability density distributions of the resonant frequencies have sharp peak and long tails as shown in Fig. 6. The nonlinear relationship between the characteristic parameters and resonant frequencies are also presented in different types of probability density distributions. In addition, the right sides of the probability density distribution in Fig. 6 have wide interval ranges, which means that the resonant frequencies of graphene fluctuate within large interval ranges. Moreover, the shapes of probability density distributions of resonant frequenceis in Fig. 6 are not similar with the commonly used Gaussian, Weibull, etc., which well confirms the analytical methods by homogeneous chaos requires more further work. The explorations based on semi-analytical approaches by stochastic perturbation methods, polynomial chaos expansion, etc. are in prospect. More precisely, the change amplitudes of the Young's modulus are 1.18% and 1.24%, in the Pearson and Spearman correlation coefficients, respectively. From the first to the tenth vibration order, the largest reduction percentages are only 0.51% for the Pearson and 0.67% for Spearman correlation coefficients to the thickness of graphene. In terms of the Posson's ratio, the discrepancies with the first order and previous order resonant vibrations is presented in Fig. 8. Specifically, the smallest relative error is 0.39%, while the largest relative error equals to 6.27%. In short, the correlation coefficients of Young's modulus and thickness of graphene are settled in the narrow interval ranges with static tendancy, while the correlation coefficients of Poisson's ratio not only have negative values but also have evident flucutations lager than 6%.
Scatter sample results. In order to analyze the complicated impacts of Poisson's ratio on the resonant frequencies, the scatter sample results are presented and compared in Fig. 9. In Fig. 9a, the thickness and Young's modulus of graphene have the simple and linear influence on the first order resonant frequency. The results in Fig. 9b,c also prove that the Poisson's ratio of graphene has nonlinear impacts on the first order resonant frequency. Even the correlation coefficients of Poisson's ratio are evidently smaller than that of Young's modulus and thickness of graphene, the nonlinear and implict relationship between the Poisson's ratio and resonant fre- www.nature.com/scientificreports/ quencies are challenging. Therefore, the Young's modulus and thickness of graphene are more correlated with the resonant frequencies, but the Poisson's ratio has more implicit and nonlinear relationship with the resonant frequencies.