The shear strain energy fluctuations caused by random porosities in graphene based on the stochastic finite element model

Strain-induced deformation is a promising strategy to modify and functionalize the material properties of graphene. However, the impacts of random porosities are inevitable and complicated in the microstructure. In order to quantify and analyze the effects of random porosities in graphene under shear stress, the energy fluctuations and the equivalent elastic modulus are computed and recorded based on the stochastic finite element model. The finite element computation is combined with the Monte Carlo stochastic sampling procedure to distribute and propagate the random porosities in pristine graphene. Two different boundary conditions are taken into consideration and compared. Furthermore, the probability statistics of shear strain energy and equivalent elastic modulus are provided based on the comparison with the results of pristine graphene. The inhomogeneous spatial randomness is founded in the statistic records of shear strain energy. The sensitivity to the graphene chirality and boundary conditions are also shown for the porous graphene. The work in this paper provides important references for strain-induced engineering and artificial functionalization through topological vacancy control in graphene.


Introduction
Strain-induced deformation has great application prospects because of the potentials in material property modification and functionalization. Shear strain was found to cause the formation of wrinkles and frictions in the multilayers of graphene [1]. Shear stress in graphene is influential, which leads to the reduction of infrared reflectivity, transformation into amorphous state and carbon onions, generation of a larger number of defects, fracture into nanometer clusters and fragmentations [2]. The shear modulus, fracture strength and related strain of both zigzag and armchair graphene vary according to different temperatures [3]. Under the elastic shear stress, the reversible phase transformation in the hexagonal lattice of graphene is implemented even if the shear strain is up to 35%, which is expected in the electromechanical systems [4]. Moreover, shear strain plays an essential role in the formation of ridge defects, which exhibit the associated Moiré superlattices and striations [5]. However, with the 30% hydrogen coverage and distribution, the shear strength and fracture strain of graphene are reduced by 50% [6]. Therefore, the sensitivity to shear stress provides feasibility for graphene property manipulation and functionalization.
The application of shear strain-induced deformation in graphene is investigated from the aspects of thermal conductivity, transport properties, flexible electronics, structural stability, etc. First, the thermal conductivity is feasible to be manipulated based on the shear strain-induced wrinkling deformation with the orientation dependence [7]. The shear strain in graphene nanoribbon influences the phonon density of states, whose higher frequency peaks (G band) shift to the low frequency and leads to a decrease in thermal conductivity [8]. In addition, a shear strain of up to 16.7% is applied to induce transport property modulation of graphene [9]. Moreover, in flexible electronics, the shear properties of graphene are controllable based on the Kirigami technique [10]. The combined effects of shear and atomic relaxation of the twisted bilayer graphene under the magic angle on the electronic structure are studied [11]. Furthermore, the interfacial strength between graphene sheets in the composite matrix is improved by shear strain-induced wrinkles [12]. In comparison, the stability of armchair graphene nanoribbons is improved with the increase of ribbon width under the shear strain [13]. Last but not least, the shear strain causes the rippling instability with strain-dependent direction and wavelength, and the large strain-induced shifts of the split components of the G optical phonon line, which may serve as a shear diagnosis [14].
For the reliable and successful implementation of shear strain engineering in graphene, more related research is conducted and reported. The pressurized microscale bubble loading devices are proposed to measure the interlayer shear stress of bilayer graphene [15]. In addition, the internal frictions and shear modulus of synthesized graphene films are also measured by the experiments [16]. The similarity in the shear modes of the layered materials provides a direct probe of interlayer interactions in the Raman spectroscopy [17]. For theoretical analysis and numerical simulation, the shear deformation and displacements are computed to discuss the dependence of wrinkling geometric patterns on the chirality of graphene [18]. The post-wrinkling behavior of square graphene nanoribbons under the in-plane shear deformation is discussed [19]. In addition, the dependences of the wrinkling amplitude, wavelength and out-of-plane displacements are observed [20]. Lattice shearing at the boundaries between grain domains with different orientations is also studied to analyze the impacts [21]. However, the random defects and porosities in graphene are implicit but essential issues requiring more concern.
The impacts of topological defects and random porosities on shear strain engineering are complicated. The thermal conductivity of graphene with low dislocation density can be modulated up to 23%, which is evidently higher than that of pristine graphene [22]. The unique topological line defects in zigzag graphene nanoribbons not only influence the local shear strain on the mechanical but also affect that on the electronic properties [23]. The stability and wrinkling of defective graphene sheets under shear deformation are not as explicit as that of pristine graphene [24]. The consideration of inhomogeneous spatial randomness in graphene is important for shear strain engineering [25].
In this paper, the random atomic vacancy defects are propagated and distributed in the honeycomb lattice of pristine graphene based on the proposed stochastic finite element model. The fluctuations and variations in the total strain energy of porous graphene are recorded by the large stochastic sample space. In addition, the impact sensitivity of the boundary conditions is taken into consideration by the chirality comparison. Both the statistic results and probability density distributions of the computed samples are provided to analyze the spatial randomness of atomic vacancy defects in graphene. A short summary is concluded to offer meaningful guidance and reference for strain-induced engineering under shear stress.

Graphene with random porosities
The special periodic honeycomb lattice of graphene leads to its extraordinary properties in thermal, electronic, optical and mechanical aspects [26]. The characteristic parameters corresponding with the geometrical properties include the thickness, length of the hexagon edge, and the length and width of graphene sheets [27]. According to the measurements and computation in the reported literature [28,29], the mentioned geometrical parameters are settled with certain values for parallel comparison. Other existing types of defects (topological defects, random porosities, grain boundaries, dislocations, etc) are really important to affect the mechanical, thermal, electrical and other properties of graphene. This study focuses on the vacancy defects for the feasibility in artificial manipulation and functionalization. In order to propagate the random porosities in the pristine graphene, the selected atoms and the related three carbon covalent bonds are supposed to be disappeared for the implementation of atomic vacancy defects, as shown in figure 1. All the carbon atoms in the lattice of graphene sheets are sequentially numbered, while the random numbers are generated by the Monte Carlo stochastic sampling procedure. The generated random numbers reflect the graphene lattice to select certain atoms for the implementation of vacancy defects. Figure 1 presents an example of graphene sheets with randomly distributed porosities.

Finite element computation
The geometrical configurations of porous graphene are figured out using the above-mentioned method. The repeated stochastic sampling loops are performed to simulate all the possible geometrical configurations of porous graphene. In addition, the element selection for geometrical model meshing, loading of the boundary conditions and computation of the mechanical response are essential steps in the finite element model. In order to take both the bending and shear stresses into consideration, the Beam 188 (Timoshenko beam) is chosen in the ANSYS Parameter design language (Version 14.5).
For the displacement and deformation computation in the finite element model of the Timoshenko beam, the deflection ω and rotation θ have the independent interpolation, where n is the number of element nodes, N i is the Lagrangian interpolation polynomials. According to the principle of minimum potential energy, where,  2 , and Q, respectively. P 1 , P 2 , P 3 , P 4 , and P 5 represent 0.1%, 0.5%, 1%, 3%, and 5%, respectively).
In the above equations, E and G are the elastic modulus and shear modulus, respectively; I is the moment of inertia, while A and l are the area of cross-section and length of the beam, respectively. Besides, There are two nodes in each beam element, n = 2, therefore, where, Then, according to the Gaussian integration, the K e s matrix can be written as, Based on the finite element method, the strain energy and vector sum of displacement for each node in the honeycomb lattice of graphene are computed and recorded as S j and D j , where j = 1, 2, . . . , N n , and N n is the number of nodes in the finite element model of graphene. Therefore, the total the strain energy and vector sum of displacement can be expressed as: (2.13) In order to analyze the fluctuations and variations of the strain energy caused by the random atomic vacancy defects, the finite element computation is performed in all the stochastic samples of porous graphene. For each stochastic sample, the corresponding S and D 2 are feasible to be computed and tracked, and S and D 2 are expanded as S k and D 2 k , where k = 1, 2, . . . , N m , and N m is the number of Monte Carlo stochastic samples for porous graphene. The results of pristine graphene are recorded as S 0 and D 2 0 for further comparison. Therefore, the following relative results are computed as, (2.14) Furthermore, the quotient of total strain energy and the sum of displacement squares are also introduced for pristine graphene and stochastic samples of porous graphene Therefore, the dimensionless results S k , D 2 k and Q k are the important index in the result comparison and discussion to observe the impacts of random porosities in graphene.

Stochastic finite element model
Different from the conventional stochastic finite element model for mechanical response computation, the proposed stochastic finite element model for the prediction of strain energy fluctuations in the porous graphene involves two layers of stochastic sampling procedures: one is the random distribution of porosities in the pristine graphene, and the other is the repeated stochastic sampling of the porous graphene, as presented in the flowchart.
In the flowchart, the inner layer stochastic sampling procedure is for the geometrical configurations of porous graphene, which is for the propagation and distribution of the random atomic vacancy defects in the honeycomb lattice of pristine graphene. In comparison, the outer layer stochastic sampling procedure is the repeated finite element computation for a large number of samples of porous graphene sheets. The implementation of the stochastic finite element model for the porous graphene is expressed as 7 steps.
Step 1. Define the characteristic geometrical parameters, and figure out the periodic honeycomb lattice of pristine graphene with the geometrical parameters.
Step 2. Number all the carbon atoms in the pristine graphene in sequence. Determine the percentage of porosities and the corresponding number of the atomic vacancy defects in the lattice of graphene.
Step 3. Generate the random numbers according to the percentage of the atomic vacancy defects, and reflect the random numbers in the lattice of graphene to select the random locations of the atomic vacancy defects. Step 4. Perform the disappearance of carbon atoms and the related carbon covalent bonds for all the selected random atoms, then form the geometrical configurations of porous graphene.
Step 5. Mesh the porous graphene with Beam 188, and define the boundary conditions for the shear stress.
Step 6. Calculate the displacements and strain energy of porous graphene by the finite element model.
Step 7. Repeat the finite element computation for a large number of stochastic samples of porous graphene, and record the sum of displacements and the total strain energy of nodes in the finite element model for each porous graphene (figure 2).

Results and discussion
Based on the stochastic finite element model of porous graphene, numerous random samples are performed to track and observe the strain energy fluctuations. In order to discuss the effects of chirality of graphene, two boundary conditions are taken into consideration, and the computed results are compared.

Boundary condition (1)
The boundary condition (1) is defined as presented in figure 3, the freedom degrees of nodes on the left edge of the graphene sheet are fixed, that is, the displacements and rotations of the X, Y, Z axial directions of nodes on the left edge of graphene sheet are settled as zero. In comparison, the nodes on the right edge are loaded by the force of 100 N in the vertical upward direction. With the increase of the atomic vacancy defects, the displacements of porous graphene become larger. The atomic vacancy defects change the local displacements, especially the places near the right edge of the porous graphene. In order to track and record the energy fluctuations of porous graphene, 5000 samples are performed under the first boundary condition. In figure 4, the magnitudes of Y axis directly present the discrepancies between porous graphene and pristine graphene under the same boundary condition. By comparing the results of pristine graphene, the influences caused by the magnitudes of loading forces on the right edge are neglectable. As defined above, the total strain energy of all nodes compared with that of pristine graphene, as well as the relative sum of displacement squares of each porous graphene sample are calculated and recorded as S and D 2 . In addition, the quotient Q of the total strain energy with the sum of displacement squares are computed by the comparison with that of pristine graphene, which is an important index to evaluate the stiffness of the graphene sheets.
According to the results in figure 4, both the total strain energy and sum of displacement squares are amplified with the increase of the number of atomic vacancy defects. The mean of total strain energy of P 1 -P 5 reaches 1.0183, 1.0435, 1.0877, 1.2978, 1.5854 times of the results of pristine graphene, respectively; while the mean of the sum of displacement squares of P 1 -P 5 presents 1.0320, 1.0767, 1.1577, 1.5765, and 2.2408 times of that of pristine graphene. In addition, with the increase of the number of atomic vacancy defects in porous graphene, the interval ranges of both S and D 2 also rise, which means the strain energy and sum of displacement squares vary in the wider intervals. However, the quotient Q decreases sharply with the increase of atomic vacancy defects. Therefore, the increase of atomic vacancy defects leads to the amplification of the total strain energy and sum of displacement squares in the wider variation intervals, but the stiffness of porous graphene reduces.
In figure 5, the probability density distributions of porous graphene samples are presented. For P 1 , P 2 and P 3 , the histograms of S, D 2 , and Q are clustered closely. Furthermore, probability density distributions of Q in figure 5(c) are closer to the Gaussian distribution compared with that in figures 5(a) and (b). Furthermore, the skewness in figures 5(a) and (b) for the probability density distributions is more evident than that of figure 5(c). It is evident that the randomness of the total strain energy and sum of displacement squares is neither symmetrical nor homogeneous. In order to further analyze the randomness caused by the stochastically distributed porosities, the statistic computation is performed based on the data of the stochastic finite element model.

Boundary condition (2)
For the second boundary condition, the freedom degrees of the nodes on the lower edge are all fixed, and the displacements and rotation of the fixed nodes in X, Y, and Z axial directions are supposed to be zero. The forces are loaded to the nodes on the upper edge of the graphene sheets, and the value of the shear force is also 100 N in the horizontal right direction. The contour results of displacements are different from those under the first boundary condition. The sensitivity to the boundary conditions of the shear displacements and deformation is also proved in the reported literature [28] (figure 6).
Furthermore, both the statistic results and probability density distributions under the second boundary condition are provided, as shown in figures 7 and 8. The results in figure 7 reach agreements with that in figure 4. The increase of the atomic vacancy defects leads to the augmentations of the total strain energy, sum of displacement squares, and fluctuations in the wider interval ranges. However, the quotient Q decreases evidently when the percentage of atomic vacancy defects is more than 1%. The mean of G in figure 7 decreases from 0.9852 to 0.6856 for graphene sheets with 0.1% and 5% of atomic vacancy defects.
In addition, the probability density distributions of S and D 2 are presented in figures 8(a) and (b). The differences in figures 8 and 5 reveal the effects of the boundary conditions. Both the S and D 2 are sensitive to the boundary conditions of the porous graphene. Moreover, the probability density distributions of S and D 2 also present evident skewness as shown in figures 5(a) and (b). The effects of random atomic vacancy defects in porous graphene are neither homogeneous nor symmetrical on S and D 2 . However, the probability density distributions of Q are more symmetrical and approximated to the Gaussian distribution, which will be compared based on the skewness and kurtosis computation.

Comparison and discussion
The variations and fluctuations in the total strain energy and sum of displacement squares caused by the location of random atomic vacancy defects are presented with the corresponding probability density distributions and related interval ranges. In order to further explore the statistical characteristics of the random sample space, the kurtosis, skewness, mean and variance of S, D 2 , and Q are computed according to the random samples of the stochastic finite element model. In addition, since S, D 2 , and Q are computed by comparison with that of pristine graphene, the presented results are relative values, and the mean values of S, D 2 , and Q vary around 1, as shown in tables 1 and 2.
According to tables 1, 2 and figure 9, the characteristic statistic results demonstrate that the porous graphene sheets under the first boundary condition provide smaller S and D 2 than that under the second boundary condition when the number of atomic vacancy defects is equal. It means that the shear stress loading on the armchair edge of porous graphene cause lower total strain energy and displacements in the statistical term. In addition, the variances of S and D 2 in figures 9(c) and (d) also prove that the porous graphene with the shear stress on the armchair edge is more robust than that on the zigzag edge. Therefore, loading the shear stress on the zigzag edge of porous graphene presents larger strain energy storage, while loading the shear stress on the armchair edge is more robust to defense displacements and variance.
In the term of quotient of the total strain energy and sum of displacement squares, the more important information is presented for strain-induced engineering. Firstly, the mean of Q under the first boundary condition has a larger value than that under the second boundary condition. With the increase of atomic vacancy defects, the reduction of Q of the porous graphene with shear stress on the armchair edge is less than that on the zigzag edge, that reaches good agreements with the points mentioned above.
Moreover, the skewness in figure 10(d) reveals that the probability density distribution of Q under the first boundary condition is more symmetrical than that under the second boundary condition, since the values of the red curve are more approaching zero than the values of the blue curve. In addition, according to the statistic theory, when the Kurtosis of sample space is closer to 3, it means that the probability density distribution is more similar to the Gaussian distribution. Then, as shown in figure 10(c), the probability density distribution of Q under the second boundary condition are closer to the Gaussian distribution, especially when the percentage of atomic vacancy defects is higher than 1%.

Conclusion
Based on the proposed stochastic finite element model, the random porosities are successfully distributed in the periodic honeycomb lattice of graphene. Under two different boundary conditions, the total strain energy and other corresponding results are compared with that of pristine graphene. According to the statistical results of the porous graphene, the following points can be concluded for strain-induced engineering.
• The total strain energy and sum of displacement squares vary in wide intervals with the increase of the atomic vacancy defects. • The effects of the random distributed atomic vacancy on the total strain energy and the sum of displacement squares are neither homogeneous nor symmetrical. • The probability density distributions of the quotient of the total strain energy and sum of displacement squares are approximated to the Gaussian distribution, especially when the percentage of atomic vacancy defects is larger than 1%. • The chirality of porous graphene is an important factor to be taken into consideration. The fluctuations and variations in the total strain energy and sum of displacement squares are sensitive to the boundary conditions. • Loading the shear stress on the armchair edge is more robust to defense displacements and variance, while loading the shear stress on the zigzag edge has greater potential for strain energy storage. The work in this study provides a meaningful reference for strain-induced engineering and artificial functionalization through topological vacancy control in graphene.