The Impacts of Random Distributed Vacancy Defects in Steady-State Thermal Conduction of Graphene

: The unavoidable vacancy defects dispersed throughout the entire pristine graphene tailor to the integrity of the lattice structure and thereby have complicated impacts on the mechanical and thermal properties of graphene. In order to analyze the inﬂuence of vacancy defects on the extraordinary thermal conductivity of graphene, three typical kinds of vacancy defects—namely center concentrated, periodic, and random distributed vacancy defects—are compared and discussed. In the steady-state thermal conduction, the ﬁnite element method (FEM) is performed to calculate the total thermal energy and temperature ﬁeld. The equivalent coe ﬃ cient of thermal conductivity is derived from thermal energy, amount of vacancy defects, and boundary condition. The chirality in graphene is discussed by the location of its heat source. Moreover, the Monte Carlo simulation is applied to propagate the uncertainty of random vacancy defects in the ﬁnite element model of pristine graphene. In this paper, we provide the robustness to defend the impacts of vacancy defects on thermal conduction and the ﬂuctuation and divergence caused by a certain number of random vacancy defects.


Introduction
The capability of heat conduction in materials is rooted in its atomic structure. Carbon materials form a variety of allotropes. The allotropes of carbon have completely different orders of magnitude in thermal conductivity, from 0.01 W/mK in amorphous carbon to 2000 W/mK in diamond or graphene in room temperature [1,2]. The mechanical exfoliation of graphene [3] and discovery of its extraordinary electrical conduction [4][5][6] provide strong evidence of thermal transportation in two-dimensional (2D) crystals.
Thermal conductivity is one of the essential issues in the exploration and development of optoelectronic and photonic devices. In the integrated circuits and three-dimensional electronics, thermoelectric energy conversion strongly suppressed thermal conductivity [7]. From the measurement of optothermal Raman in suspended graphene from mechanical exfoliation, the value of thermal conductivity is found to exceed 3000 W/mK [8][9][10]. The phonon mean-free path was estimated to be 775 nm near room temperature [11]. For high-quality chemical vapor deposited (CVD) graphene, the thermal conductivity is also beyond 2500 W/mK at 350 K and is as high as 1400 W/mK at 500 K [12]. When the environmental temperature is at 600 K, the thermal conductivity becomes 630 W/mK for where K x is the conductivity constant in the x direction and the lengths of the volume elements in y and z direction are represented as δy, δz respectively. The leaving quantity of heat can be expressed by Taylor's series as: Excluding the high order terms ε, the net heat inflow in x direction can be approximately computed as: where δV = δxδyδz is the volume. According to the heat inflow in x direction, the sum of the net amount of heat stored due to difference in conduction heat flow can be expressed as: K y , K z are the conductivity constants in the y and z direction, respectively. Besides, the heat generation in the element at a rate of per unit volume per unit time is introduced as q g : For homogeneous materials, K x = K y = K z . When there is no heat generation q g = 0, the above equation is simplified as: For the steady-state heat conduction, the derivative of the temperature with time is zero. Then, the heat conduction equation can be written as:

Thermal Conductivity
The thermal conductivity in matter is one of the natural properties. For graphene, optothermal Raman, and electrical measurements are performed to test the thermal conductivity of suspended graphene layers from exfoliation or chemical vapor deposition (CVD). Besides, MD numerical simulation and analytical methods based on the Boltzmann transport equation in valence force field are explored to determine the precise value of thermal conductivity. In Table 1, the ranges of thermal conductivity for graphene are listed from different experimental tests and simulation methods. Large deviation is demonstrated for thermal conductivity from 600 W/mK to 10000 W/Mk. However, the results of Nika [21,36], Savin [25] and Lindsay [31] are more approximated consistent to the measurements results of Balandin [8], Ghosh [10], Jaugeras [14], and Cai [12] than the result of MD [29], as shown in Figure 1. Therefore, the values of thermal conductivity for graphene are settled in the reported range of Balandin [8] and Nika [21].
The heat in graphene is mostly carried by acoustic phonons rather than by electrons. The heat flux along a graphene flake can be computed as [21]: Appl. Sci. 2019, 9,2363 4 of 17 The above equation can be used to calculate the thermal conductivity with the actual dependence on the phonon frequency ω s ( → q ), the phonon velocity → v (s, → q ), and the number of the phonon with the certain frequencies n( → q , ω s ).
→ v ω is the energy carried by one phonon and N(ω, → q ) is the number of phonons in the heat flux. From the macroscopic definition: where L x = d is the width of graphene flake, L y is the sample length, and h is supposed to be 0.35 nm as the thickness of graphene. The thermal conductivity tensor κ αβ in Equation (9) is expressed as: where τ tot is the combined phonon relaxation rate in graphene. Moreover, α, β are the subscripts that represent the elements in the matrix of thermal conductivity tensor. Similarly, the diagonal element of the thermal conductivity tensor is computed as: Table 1. Thermal conductivity of graphene in reference.

Finite Element Model for Graphene
FEM is applied for spatial discretization. In the weighted residual method, trial functions serve as weighting functions. In the finite element model, the conducting bars in graphene lattice is defined by two nodes, a cross-sectional area, and the material properties, as in Figure 2. Specific heat and density are ignored for steady-state solutions. The thermal conductivity is in the element longitudinal direction.  (11) Then, the 2D density of phonon states is taken into consideration to obtain the expression for the scalar thermal conductivity: Based on Equation (12), the dispersion of phonon and phonon polarization branches are taken into consideration using a valence force field. The room temperature is substituted into the theoretical expression of Equation (12). Then, the value of thermal conductivity in this study is set to be 3000 W/mK for steady-state heat conduction.

Finite Element Model for Graphene
FEM is applied for spatial discretization. In the weighted residual method, trial functions serve as weighting functions. In the finite element model, the conducting bars in graphene lattice is defined by two nodes, a cross-sectional area, and the material properties, as in Figure 2. Specific heat and density are ignored for steady-state solutions. The thermal conductivity is in the element longitudinal direction. The conductivity matrix of element can be expressed as: where A , K , and L is the area, thermal conductivity, and distance between nodes, respectively. The output is computed as: where q , i T , i T , and Q is the thermal flux, temperature at node i and j, and Q is the heat rate.
Energies are available in the solution for each element.   The conductivity matrix of element can be expressed as: where A, K, and L is the area, thermal conductivity, and distance between nodes, respectively. The output is computed as: where q, T i , T i , and Q is the thermal flux, temperature at node i and j, and Q is the heat rate. Energies are available in the solution for each element.
where N is the number of integration points, {σ} is the stress vector, ε ei is the elastic strain vector,  Figure 3 shows the finite element model based on honeycomb lattice structure of graphene. Its heat source is discussed in armchair and zigzag edges, respectively, as shown in Figure 3a, Figure 3 shows the finite element model based on honeycomb lattice structure of graphene. Its heat source is discussed in armchair and zigzag edges, respectively, as shown in Figure 3a,b. The colorful contour depicts the thermal conduction in the steady-state.  Vacancy defects are an unavoidable existence in graphene. However, the shape, size, and location in the entire graphene are unpredictable and uncertain. In order to discuss the impacts of vacancy defects in the steady-state thermal conduction of graphene, three different vacancy defects are compared and discussed in the following results. As shown in Figure 4, center concentrated vacancy defects, periodic vacancy defects, and random distributed vacancy defects are included as three typical vacancy defects in graphene.

Results and Discussion
In the following section, the three kinds of vacancy defects are sequentially presented.

Center Concentrated Vacancy Defects
The concentrated vacancy defects in the center of graphene are gradually amplified, as shown in Figure 5. According to the size amplification of the center concentrated vacancy defects, the thermal energy of entire graphene gradually decreases in both situations, no matter what its heat Vacancy defects are an unavoidable existence in graphene. However, the shape, size, and location in the entire graphene are unpredictable and uncertain. In order to discuss the impacts of vacancy defects in the steady-state thermal conduction of graphene, three different vacancy defects are compared and discussed in the following results. As shown in Figure 4, center concentrated vacancy defects, periodic vacancy defects, and random distributed vacancy defects are included as three typical vacancy defects in graphene.  Figure 3 shows the finite element model based on honeycomb lattice structure of graphene. Its heat source is discussed in armchair and zigzag edges, respectively, as shown in Figure 3a,b. The colorful contour depicts the thermal conduction in the steady-state. Vacancy defects are an unavoidable existence in graphene. However, the shape, size, and location in the entire graphene are unpredictable and uncertain. In order to discuss the impacts of vacancy defects in the steady-state thermal conduction of graphene, three different vacancy defects are compared and discussed in the following results. As shown in Figure 4, center concentrated vacancy defects, periodic vacancy defects, and random distributed vacancy defects are included as three typical vacancy defects in graphene.

Results and Discussion
In the following section, the three kinds of vacancy defects are sequentially presented.

Center Concentrated Vacancy Defects
The concentrated vacancy defects in the center of graphene are gradually amplified, as shown in Figure 5. According to the size amplification of the center concentrated vacancy defects, the thermal energy of entire graphene gradually decreases in both situations, no matter what its heat

Results and Discussion
In the following section, the three kinds of vacancy defects are sequentially presented.

Center Concentrated Vacancy Defects
The concentrated vacancy defects in the center of graphene are gradually amplified, as shown in Figure 5. According to the size amplification of the center concentrated vacancy defects, the thermal energy of entire graphene gradually decreases in both situations, no matter what its heat source is in the armchair or zigzag edges. The thermal energy of graphene with its heat source in armchair edges is larger than that of graphene with its heat source in zigzag edges, as shown in Figure 6a. The location of its heat source, as one of the boundary conditions, contributes to the difference of thermal energy in the steady-state thermal conduction of graphene. Furthermore, the definition of equivalent coefficient of thermal conductivity for vacancy defected graphene can be expressed as: armchair edges is larger than that of graphene with its heat source in zigzag edges, as shown in Figure  6a. The location of its heat source, as one of the boundary conditions, contributes to the difference of thermal energy in the steady-state thermal conduction of graphene. Furthermore, the definition of equivalent coefficient of thermal conductivity for vacancy defected graphene can be expressed as: In addition, the percentage of vacancy defects in the entire graphene can be counted as: where n and m are the number of elements in pristine graphene and vacancy defected graphene.
The total thermal energy is computed by the sum of the thermal energy of each element and T pi E in pristine graphene and vacancy defected graphene, independently. Based on Equation (18), the amount of vacancy defects is introduced in an equivalent coefficient of the thermal conduction. In Figure 6b, the reduction of an equivalent coefficient of thermal conductivity becomes sharper with the size amplification of a concentrated center. Under both boundary conditions, the impact of vacancy defects in a steady-state thermal conduction is non-linear. Further, the gradient of curve in the equivalent coefficient is augmented with the enlargement of center concentrated vacancy defects. Moreover, the vacancy defected graphene sheets with its heat source in armchair and zigzag edges have consistent results in the thermal condition's equivalent coefficient, especially when the size of the center concentrated vacancy defects is not large. With the size enlargement of the center concentrated vacancy defects, the difference between two boundary conditions becomes evident. The graphene with its heat source in 6a. The location of its heat source, as one of the boundary conditions, contributes to the difference of thermal energy in the steady-state thermal conduction of graphene. Furthermore, the definition of equivalent coefficient of thermal conductivity for vacancy defected graphene can be expressed as: In addition, the percentage of vacancy defects in the entire graphene can be counted as: where n and m are the number of elements in pristine graphene and vacancy defected graphene.
The total thermal energy is computed by the sum of the thermal energy of each element and T pi E in pristine graphene and vacancy defected graphene, independently. Based on Equation (18), the amount of vacancy defects is introduced in an equivalent coefficient of the thermal conduction. In Figure 6b, the reduction of an equivalent coefficient of thermal conductivity becomes sharper with the size amplification of a concentrated center. Under both boundary conditions, the impact of vacancy defects in a steady-state thermal conduction is non-linear. Further, the gradient of curve in the equivalent coefficient is augmented with the enlargement of center concentrated vacancy defects. Moreover, the vacancy defected graphene sheets with its heat source in armchair and zigzag edges have consistent results in the thermal condition's equivalent coefficient, especially when the size of the center concentrated vacancy defects is not large. With the size enlargement of the center concentrated vacancy defects, the difference between two boundary conditions becomes evident. The graphene with its heat source in In addition, the percentage of vacancy defects in the entire graphene can be counted as: where n and m are the number of elements in pristine graphene and vacancy defected graphene. The total thermal energy is computed by the sum of the thermal energy of each element and E T pi in pristine graphene and vacancy defected graphene, independently.
Based on Equation (18), the amount of vacancy defects is introduced in an equivalent coefficient of the thermal conduction. In Figure 6b, the reduction of an equivalent coefficient of thermal conductivity becomes sharper with the size amplification of a concentrated center. Under both boundary conditions, the impact of vacancy defects in a steady-state thermal conduction is non-linear. Further, the gradient of curve in the equivalent coefficient is augmented with the enlargement of center concentrated vacancy defects. Moreover, the vacancy defected graphene sheets with its heat source in armchair and zigzag edges have consistent results in the thermal condition's equivalent coefficient, especially when the size of the center concentrated vacancy defects is not large. With the size enlargement of the center concentrated vacancy defects, the difference between two boundary conditions becomes evident. The graphene with its heat source in armchair edges has higher equivalent coefficient of thermal conductivity than that with its heat source in zigzag edges.
In addition, the contour of temperature in graphene is revealed in Figures 7 and 8. In order to track the exact impacts of center concentrated vacancy defects in the steady-state heat conduction of graphene, the local region around the center concentrated vacancy defects is compared. When the size of the vacancy defects is small, the effects in temperature are not obvious. Along with the size amplification of center concentrated vacancy, the influenced region in graphene become larger and the regularity and band characteristic of temperature results in pristine graphene is destroyed under both boundary conditions. The regions close to the center concentrated vacancy defects tend to have a higher temperature than the pristine graphene. Therefore, the center concentrated vacancy defects in graphene is a crucial factor to change the temperature field in the steady-state thermal conduction. Heat accumulation near the center concentrated vacancy defects in graphene leads to local higher temperature and can cause thermal damages in devices.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 8 of 17 armchair edges has higher equivalent coefficient of thermal conductivity than that with its heat source in zigzag edges. In addition, the contour of temperature in graphene is revealed in Figures 7 and 8. In order to track the exact impacts of center concentrated vacancy defects in the steady-state heat conduction of graphene, the local region around the center concentrated vacancy defects is compared. When the size of the vacancy defects is small, the effects in temperature are not obvious. Along with the size amplification of center concentrated vacancy, the influenced region in graphene become larger and the regularity and band characteristic of temperature results in pristine graphene is destroyed under both boundary conditions. The regions close to the center concentrated vacancy defects tend to have a higher temperature than the pristine graphene. Therefore, the center concentrated vacancy defects in graphene is a crucial factor to change the temperature field in the steady-state thermal conduction. Heat accumulation near the center concentrated vacancy defects in graphene leads to local higher temperature and can cause thermal damages in devices.

Periodic Vacancy Defects
Different with the pattern of center concentrated vacancy defects, periodic vacancy defects are gradually dispersed in one quarter of a graphene sheet, as presented in Figure 9. The thermal energy and equivalent coefficient of thermal conductivity are computed by the finite element model in graphene with its heat source in armchair and zigzag edges.  armchair edges has higher equivalent coefficient of thermal conductivity than that with its heat source in zigzag edges. In addition, the contour of temperature in graphene is revealed in Figures 7 and 8. In order to track the exact impacts of center concentrated vacancy defects in the steady-state heat conduction of graphene, the local region around the center concentrated vacancy defects is compared. When the size of the vacancy defects is small, the effects in temperature are not obvious. Along with the size amplification of center concentrated vacancy, the influenced region in graphene become larger and the regularity and band characteristic of temperature results in pristine graphene is destroyed under both boundary conditions. The regions close to the center concentrated vacancy defects tend to have a higher temperature than the pristine graphene. Therefore, the center concentrated vacancy defects in graphene is a crucial factor to change the temperature field in the steady-state thermal conduction. Heat accumulation near the center concentrated vacancy defects in graphene leads to local higher temperature and can cause thermal damages in devices.

Periodic Vacancy Defects
Different with the pattern of center concentrated vacancy defects, periodic vacancy defects are gradually dispersed in one quarter of a graphene sheet, as presented in Figure 9. The thermal energy and equivalent coefficient of thermal conductivity are computed by the finite element model in graphene with its heat source in armchair and zigzag edges.

Periodic Vacancy Defects
Different with the pattern of center concentrated vacancy defects, periodic vacancy defects are gradually dispersed in one quarter of a graphene sheet, as presented in Figure 9. The thermal energy and equivalent coefficient of thermal conductivity are computed by the finite element model in graphene with its heat source in armchair and zigzag edges.
With the increase in the amount of periodic vacancy defects, the thermal energy in graphene with its heat source in the armchair and zigzag edges generally decrease. The thermal energy in the graphene with heat source in the armchair edges are absolutely larger than that in the graphene with its heat source in zigzag edges. However, the results of the equivalent coefficient of thermal conductivity is more complicated than that of the center concentrated vacancy defects, as shown in Figure 10. The equivalent coefficient of the thermal conductivity in graphene with its heat source in armchair edges has a strengthening effect, owing to periodic vacancy defects. Furthermore, the equivalent coefficient of thermal conductivity in graphene with its heat source in zigzag edges does not monotonically decrease with the development of periodic vacancy defects. The fluctuation exists in the second mode of periodic vacancy defects.

Periodic Vacancy Defects
Different with the pattern of center concentrated vacancy defects, periodic vacancy defects are gradually dispersed in one quarter of a graphene sheet, as presented in Figure 9. The thermal energy and equivalent coefficient of thermal conductivity are computed by the finite element model in graphene with its heat source in armchair and zigzag edges.  With the increase in the amount of periodic vacancy defects, the thermal energy in graphene with its heat source in the armchair and zigzag edges generally decrease. The thermal energy in the graphene with heat source in the armchair edges are absolutely larger than that in the graphene with its heat source in zigzag edges. However, the results of the equivalent coefficient of thermal conductivity is more complicated than that of the center concentrated vacancy defects, as shown in Figure 10. The equivalent coefficient of the thermal conductivity in graphene with its heat source in armchair edges has a strengthening effect, owing to periodic vacancy defects. Furthermore, the equivalent coefficient of thermal conductivity in graphene with its heat source in zigzag edges does not monotonically decrease with the development of periodic vacancy defects. The fluctuation exists in the second mode of periodic vacancy defects. Contrary to graphene with its heat source in armchair edges, graphene with its heat source in zigzag edges has an apparent reduced equivalent coefficient of thermal conductivity in the first mode of periodic vacancy defects. Although thermal energy of graphene with its heat source in armchair edges is larger than that of graphene with its heat source in zigzag edges in all modes of periodic vacancy defects, graphene with its heat source in zigzag edges has greater equivalent coefficient of thermal conductivity after the first mode. Thus, graphene with its heat source in armchair edges has a strengthening effect caused by the first mode of periodic vacancy defects. However, it decreases abruptly with the development of periodic vacancy defects, while graphene with its heat source in zigzag edges is more robust and defends against the effects of periodic vacancy defects with fluctuation.
Moreover, the influence of periodic vacancy defects in the temperature field is not as evident as that in the center concentrated vacancy. In Figures 11 and 12, the changes in the temperature field of the local region near the periodic vacancy defects in graphene with its heat source in armchair and zigzag edges are not obvious when compared with the results in the original pristine graphene. Contrary to graphene with its heat source in armchair edges, graphene with its heat source in zigzag edges has an apparent reduced equivalent coefficient of thermal conductivity in the first mode of periodic vacancy defects. Although thermal energy of graphene with its heat source in armchair edges is larger than that of graphene with its heat source in zigzag edges in all modes of periodic vacancy defects, graphene with its heat source in zigzag edges has greater equivalent coefficient of thermal conductivity after the first mode. Thus, graphene with its heat source in armchair edges has a strengthening effect caused by the first mode of periodic vacancy defects. However, it decreases abruptly with the development of periodic vacancy defects, while graphene with its heat source in zigzag edges is more robust and defends against the effects of periodic vacancy defects with fluctuation.
Moreover, the influence of periodic vacancy defects in the temperature field is not as evident as that in the center concentrated vacancy. In Figures 11 and 12, the changes in the temperature field of the local region near the periodic vacancy defects in graphene with its heat source in armchair and zigzag edges are not obvious when compared with the results in the original pristine graphene.
with its heat source in zigzag edges is more robust and defends against the effects of periodic vacancy defects with fluctuation.
Moreover, the influence of periodic vacancy defects in the temperature field is not as evident as that in the center concentrated vacancy. In Figures 11 and 12, the changes in the temperature field of the local region near the periodic vacancy defects in graphene with its heat source in armchair and zigzag edges are not obvious when compared with the results in the original pristine graphene.

Random Distributed Vacancy Defects
The center concentrated and periodic vacancy defects are two typical vacancy defects in graphene. Actually, the vacancy defects are stochastically distributed; thus, it is difficult to predict the placement in the entire graphene. By combing the finite element model with the Monte Carlo simulation [46], the vacancy defects are randomly dispersed in graphene with a specific amount. The statistical results of thermal energy for graphene with heat source in armchair and zigzag edges are listed in Table 2. In the procedure of the Monte Carlo based finite element method, the minimum and maximum values of the thermal energy are captured from many possible situations. The mean and variance values are computed from series of thermal energy from samples of the Monte Carlo based finite element method. Different than the center concentrated and periodic vacancy defects, the minimum, maximum, and mean values of thermal energy in graphene with random distributed vacancy defects have linear declined tendency with the increase of the amount of vacancy defects. The minimum, maximum, and mean values of graphene's thermal energy with its heat source in armchair edges are larger than when its heat source is in zigzag edges, but the linear results are parallel to each

Random Distributed Vacancy Defects
The center concentrated and periodic vacancy defects are two typical vacancy defects in graphene. Actually, the vacancy defects are stochastically distributed; thus, it is difficult to predict the placement in the entire graphene. By combing the finite element model with the Monte Carlo simulation [46], the vacancy defects are randomly dispersed in graphene with a specific amount. The statistical results of thermal energy for graphene with heat source in armchair and zigzag edges are listed in Table 2. In the procedure of the Monte Carlo based finite element method, the minimum and maximum values of the thermal energy are captured from many possible situations. The mean and variance values are computed from series of thermal energy from samples of the Monte Carlo based finite element method. Different than the center concentrated and periodic vacancy defects, the minimum, maximum, and mean values of thermal energy in graphene with random distributed vacancy defects have linear declined tendency with the increase of the amount of vacancy defects. The minimum, maximum, and mean values of graphene's thermal energy with its heat source in armchair edges are larger than when its heat source is in zigzag edges, but the linear results are parallel to each other in the two different boundary conditions, as presented in Figure 13. In addition, the divergences between the minimum and maximum values of the thermal energy with different amounts of vacancy defects are non-negligible (Table 2 and Figure 13). The unit of the thermal energy is 10 5 J. With the increase of the number of random distributed vacancy defects in graphene, the differences between the minimum and maximum values of the thermal energy are amplified. The exact computation results are provided and compared in the following discussion. Furthermore, the variance of thermal energy is amplified with the enlargement in the amount of stochastically distributed vacancy defects. However, the parallel characteristic between the results of the two boundary conditions disappears in the variance of thermal energy. The variances of graphene with its heat source in armchair edges have linear growth with an increase of the number of stochastic vacancy defects. However, the graphene with its heat source in zigzag edges is tougher to shield than the variation and fluctuation in thermal energy. The gradient of variance in Figure 13d declines when Per equals 3%. The difference in variance of thermal energy in two boundary conditions becomes large.
To be more exact, the probability density distribution of thermal energy in graphene under two boundary conditions are illustrated in Figures 14 and 15. With the variance of thermal energy, the probability density is distributed in a narrow range when the number of random vacancy defects is small in both boundary conditions. When the amount of vacancy defects is amplified, the probability density covers a wider range with larger variances. When Per is 0.2% for graphene with its heat Figure 13. The statistical results of graphene with random distributed vacancy defects (the unit for thermal energy is 10 5 J; (a-d) are the minimum, maximum, mean, and variance values of the thermal energy, respectively. Furthermore, the variance of thermal energy is amplified with the enlargement in the amount of stochastically distributed vacancy defects. However, the parallel characteristic between the results of the two boundary conditions disappears in the variance of thermal energy. The variances of graphene with its heat source in armchair edges have linear growth with an increase of the number of stochastic vacancy defects. However, the graphene with its heat source in zigzag edges is tougher to shield than the variation and fluctuation in thermal energy. The gradient of variance in Figure 13d declines when Per equals 3%. The difference in variance of thermal energy in two boundary conditions becomes large.
To be more exact, the probability density distribution of thermal energy in graphene under two boundary conditions are illustrated in Figures 14 and 15. With the variance of thermal energy, the probability density is distributed in a narrow range when the number of random vacancy defects is small in both boundary conditions. When the amount of vacancy defects is amplified, the probability density covers a wider range with larger variances. When Per is 0.2% for graphene with its heat source in armchair edges, the range of probability density distribution ranges from 114.9-115.7 kJ; the range becomes 95.71-100.6 kJ when Per equals 5%. The same situation happens in graphene with its heat source in zigzag edges. The range of probability density distribution is at 92.3-93.0 kJ when Per is 0.2% and 77.18-80.50 kJ when Per equals to 5%. When the percentage of random vacancy defects is petite, the probability density is more condensed as a peak, while the probability density distribution is gentle with a longer drag when Per is large. Importantly, even though vacancy defects are stochastically dispersed in the entire graphene, the range of thermal energy cannot cover the whole result domain and may be limited in certain intervals when Per is determined. The range of thermal energy in graphene with two boundary conditions is not continuous but discrete in the specific intervals according to Per. The discontinuity (between Per is 1%, 3%, and 5%) in thermal energy is useful to distinguish the amount of vacancy defects and also can be designed as control switch in precise devices and instruments.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 12 of 17 Figure 14. The histogram of probability density distribution of thermal energy in graphene with its heat source in armchair edges (P1-P6 represents the percentage of vacancy as 0.2%, 0.5%, 0.8%, 1%, 3%, and 5%, respectively, and the unit for thermal energy is J). Figure 15. The histogram of probability density distribution of thermal energy in graphene with its heat source in zigzag edges (P1-P6 represents the percentage of vacancy as 0.2%, 0.5%, 0.8%, 1%, 3%, and 5%, respectively, and the unit for thermal energy is J).
In addition, the ratio of thermal energy reduction in graphene with two boundary conditions is compared in Figure 16. In the graphene with its heat source in armchair edges, the minimum value of thermal energy has 1%, 2.37%, 3.16%, 3.85%, 10.86%, and 17.56% of reduction compared with the thermal energy in original pristine graphene, when Per equals to 0.2%, 0.5%, 0.8%, 1%, 3%, and 5%. The decrease in the maximum value of thermal energy is smaller when Per equals 0.28%, 0.98%, 1.65%, 2.25%, 7.39%, and 13.32%, respectively. In graphene with its heat source in zigzag edges, the degrees of reduction in minimum and maximum values of thermal energy are similar to that in armchair edges. Precisely, the minimum values of thermal energy have 1.05%, 2.09%, 3.19%, 4.09%, 11.04%, and 17.28% of declines, while the maximum values of thermal energy have 0.31%, 1.03%, Figure 14. The histogram of probability density distribution of thermal energy in graphene with its heat source in armchair edges (P1-P6 represents the percentage of vacancy as 0.2%, 0.5%, 0.8%, 1%, 3%, and 5%, respectively, and the unit for thermal energy is J).
Appl. Sci. 2019, 9, x FOR PEER REVIEW 12 of 17 Figure 14. The histogram of probability density distribution of thermal energy in graphene with its heat source in armchair edges (P1-P6 represents the percentage of vacancy as 0.2%, 0.5%, 0.8%, 1%, 3%, and 5%, respectively, and the unit for thermal energy is J). Figure 15. The histogram of probability density distribution of thermal energy in graphene with its heat source in zigzag edges (P1-P6 represents the percentage of vacancy as 0.2%, 0.5%, 0.8%, 1%, 3%, and 5%, respectively, and the unit for thermal energy is J).
In addition, the ratio of thermal energy reduction in graphene with two boundary conditions is compared in Figure 16. In the graphene with its heat source in armchair edges, the minimum value of thermal energy has 1%, 2.37%, 3.16%, 3.85%, 10.86%, and 17.56% of reduction compared with the thermal energy in original pristine graphene, when Per equals to 0.2%, 0.5%, 0.8%, 1%, 3%, and 5%. The decrease in the maximum value of thermal energy is smaller when Per equals 0.28%, 0.98%, 1.65%, 2.25%, 7.39%, and 13.32%, respectively. In graphene with its heat source in zigzag edges, the degrees of reduction in minimum and maximum values of thermal energy are similar to that in armchair edges. Precisely, the minimum values of thermal energy have 1.05%, 2.09%, 3.19%, 4.09%, 11.04%, and 17.28% of declines, while the maximum values of thermal energy have 0.31%, 1.03%, Figure 15. The histogram of probability density distribution of thermal energy in graphene with its heat source in zigzag edges (P1-P6 represents the percentage of vacancy as 0.2%, 0.5%, 0.8%, 1%, 3%, and 5%, respectively, and the unit for thermal energy is J).
In addition, the ratio of thermal energy reduction in graphene with two boundary conditions is compared in Figure 16. In the graphene with its heat source in armchair edges, the minimum value of thermal energy has 1%, 2.37%, 3.16%, 3.85%, 10.86%, and 17.56% of reduction compared with the thermal energy in original pristine graphene, when Per equals to 0.2%, 0.5%, 0.8%, 1%, 3%, and 5%. The decrease in the maximum value of thermal energy is smaller when Per equals 0.28%, 0.98%, 1.65%, 2.25%, 7.39%, and 13.32%, respectively. In graphene with its heat source in zigzag edges, the degrees of reduction in minimum and maximum values of thermal energy are similar to that in armchair edges. Precisely, the minimum values of thermal energy have 1.05%, 2.09%, 3.19%, 4.09%, 11.04%, and 17.28% of declines, while the maximum values of thermal energy have 0.31%, 1.03%, 1.81%, 2.14%, 7.94%, and 13.72% of declines, respectively, when Per equals to 0.2%, 0.5%, 0.8%, 1%, 3% and 5%. Therefore, the difference caused by boundary conditions is not as obvious as that attributed to the random distribution of vacancy defects. The stochastic dispersion of vacancy defects plays a more important role in the deviation of thermal energy.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 13 of 17 Figure 16. Ratio of thermal energy reduction (A-1 and A-2 are for minimum and maximum thermal energy in graphene with its heat source in armchair edges, respectively; Z-1 and Z-2 for minimum and maximum thermal energy in graphene with its heat source in zigzag edges, respectively).
In order to analyze the robustness of graphene with random distributed vacancy defects, the index β is introduced as: where M is the mean value of thermal energy and V is the variance of thermal energy. The divergence due to the randomness of vacancy defects in graphene is measured using the variance of thermal energy; the quadratic square root of variance is written as a denominator compared with the mean value of thermal energy. In Figure 17a, the equivalent coefficient of thermal conductivity is nearly identical in the two boundary conditions. With the increment of the random vacancy defects, the equivalent coefficient of mean thermal conductivity is linearly declined; the β gradient becomes smaller with the increase of random vacancy defects. When Per is less than 1%, β decreases sharply, but when Per equals 3% and 5%, β experiences a moderate reduction. Similar with the equivalent coefficient of thermal conductivity, the results of β in the two boundary conditions are consistent with small deviation. β equals to 2.47, 1.54, 1.21, 1.12, 0.58, and 0.40 in graphene with its heat source in armchair edges, and 2.41,1.60, 1.22, 1.11, 0.54, and 0.41 when Per is 0.2%, 0.5%, 0.8%, 1%, 3%, and 5%, respectively. Thus, the robustness of graphene to resist the fluctuation is due to the randomness of vacancy defects, which sharply cuts down when Per is less than 1%. Figure 16. Ratio of thermal energy reduction (A-1 and A-2 are for minimum and maximum thermal energy in graphene with its heat source in armchair edges, respectively; Z-1 and Z-2 for minimum and maximum thermal energy in graphene with its heat source in zigzag edges, respectively).
In order to analyze the robustness of graphene with random distributed vacancy defects, the index β is introduced as: where M is the mean value of thermal energy and V is the variance of thermal energy. The divergence due to the randomness of vacancy defects in graphene is measured using the variance of thermal energy; the quadratic square root of variance is written as a denominator compared with the mean value of thermal energy. In Figure 17a, the equivalent coefficient of thermal conductivity is nearly identical in the two boundary conditions. With the increment of the random vacancy defects, the equivalent coefficient of mean thermal conductivity is linearly declined; the β gradient becomes smaller with the increase of random vacancy defects. When Per is less than 1%, β decreases sharply, but when Per equals 3% and 5%, β experiences a moderate reduction. Similar with the equivalent coefficient of thermal conductivity, the results of β in the two boundary conditions are consistent with small deviation. β equals to 2.47, 1.54, 1.21, 1.12, 0.58, and 0.40 in graphene with its heat source in armchair edges, and 2.41,1.60, 1.22, 1.11, 0.54, and 0.41 when Per is 0.2%, 0.5%, 0.8%, 1%, 3%, and 5%, respectively. Thus, the robustness of graphene to resist the fluctuation is due to the randomness of vacancy defects, which sharply cuts down when Per is less than 1%.
deviation. β equals to 2.47, 1.54, 1.21, 1.12, 0.58, and 0.40 in graphene with its heat source in armchair edges, and 2.41,1.60, 1.22, 1.11, 0.54, and 0.41 when Per is 0.2%, 0.5%, 0.8%, 1%, 3%, and 5%, respectively. Thus, the robustness of graphene to resist the fluctuation is due to the randomness of vacancy defects, which sharply cuts down when Per is less than 1%.  In Figures 18 and 19, the temperature results of the samples of random vacancy defected graphene are presented. The local temperature field is not convenient when distinguishing the impacts of the random distributed vacancy defects. The discrete and stochastic vacancy defects in the usual samples cannot present the comprehensive influence in thermal conductivity. However, the thermal energy of the entire graphene is more appropriate and reliable to observe the impacts of vacancy defects. In Figures 18 and 19, the temperature results of the samples of random vacancy defected graphene are presented. The local temperature field is not convenient when distinguishing the impacts of the random distributed vacancy defects. The discrete and stochastic vacancy defects in the usual samples cannot present the comprehensive influence in thermal conductivity. However, the thermal energy of the entire graphene is more appropriate and reliable to observe the impacts of vacancy defects.

Conclusions
In this paper, we analyzed and discussed the impacts of vacancy defects in graphene for steady-state thermal conduction. Two boundary conditions as chirality in graphene are considered. The center concentrated, periodic, and random distributed vacancy defects in graphene are compared using thermal energy, the equivalent coefficient of thermal conductivity, and thermal field. The results show that: Figure 18. Temperature results of random vacancy defected graphene with its heat source in armchair edges (the percentages of vacancy defects in (a-d) are 0.2%, 0.5%, 1%, and 3%, respectively). In Figures 18 and 19, the temperature results of the samples of random vacancy defected graphene are presented. The local temperature field is not convenient when distinguishing the impacts of the random distributed vacancy defects. The discrete and stochastic vacancy defects in the usual samples cannot present the comprehensive influence in thermal conductivity. However, the thermal energy of the entire graphene is more appropriate and reliable to observe the impacts of vacancy defects.

Conclusions
In this paper, we analyzed and discussed the impacts of vacancy defects in graphene for steady-state thermal conduction. Two boundary conditions as chirality in graphene are considered. The center concentrated, periodic, and random distributed vacancy defects in graphene are compared using thermal energy, the equivalent coefficient of thermal conductivity, and thermal field. The results show that: 1. The center concentrated vacancy defects evidently destroyed the regularity and band Figure 19. Temperature results of random vacancy defected graphene with its heat source in zigzag edges (the percentages of vacancy defects in (a-d) are 0.2%, 0.5%, 1%, and 3%, respectively).

Conclusions
In this paper, we analyzed and discussed the impacts of vacancy defects in graphene for steady-state thermal conduction. Two boundary conditions as chirality in graphene are considered. The center concentrated, periodic, and random distributed vacancy defects in graphene are compared using thermal energy, the equivalent coefficient of thermal conductivity, and thermal field. The results show that: 1.
The center concentrated vacancy defects evidently destroyed the regularity and band characteristic of the temperature field in the steady-state thermal conduction. Moreover, the impacts of periodic and random distributed vacancy defects in temperature field of graphene is not convenient to track. 2.
The graphene with periodic vacancy defects is sensitive to boundary conditions. The graphene with its heat source in armchair edges has a strengthening effect in the equivalent coefficient of thermal conductivity when periodic vacancy defects is distributed in the first mode. However, graphene with its heat source in zigzag edges is more robust to defend the reduction of the equivalent coefficient of thermal conductivity in the development of periodic vacancy defects. On the contrary, center concentrated and randomly distributed vacancy defects have a consistent equivalent coefficient of thermal conductivity under two boundary conditions. 3.
Furthermore, the amount and randomness of the vacancy defects are more important than the chirality in stochastically vacancy defected graphene. When Per equals to 0.2%, the reduction of thermal energy fluctuates from 0.3-1%, and from 13-17% when Per is 5%. The thermal energy has discrete range for certain quantity of random vacancy defects. The range of thermal energy is useful when reflecting on the amount of vacancy defects designed as a control switch in devices or instruments. 4.
The robustness of graphene to resist the fluctuation owing to the randomness of vacancy defects declines sharply when Per is less than 1%. The variance caused by random dispersion of vacancy defects is amplified, according to the increment of Per.

5.
Different than the center concentrated and periodic vacancy defects, the minimum, maximum, and mean values of thermal energy in graphene with random vacancy defects have a linear decreased tendency with the increase of vacancy amounts.

Data Availability Statement
The original data used to support the findings of this study are available from the corresponding author upon request.