Luiz Felipe C. Pereira*a,
Bohayra Mortazavib,
Meysam Makaremic and
Timon Rabczukde
aDepartamento de Física Teórica e Experimental, Universidade Federal do Rio Grande do Norte, 59078-970 Natal, Brazil. E-mail: pereira@fisica.ufrn.br
bInstitute of Structural Mechanics, Bauhaus-Universität Weimar, Marienstr. 15, D-99423 Weimar, Germany. E-mail: bohayra.mortazavi@gmail.com
cChemical Engineering Department, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA
dCollege of Civil Engineering, Tongji University, 1239 Siping Road, Shanghai 20092, PR China
eChair of Computational Mechanics, Bauhaus-Universität Weimar, Marienstraße 15, 99423 Weimar, Germany
First published on 8th June 2016
Phagraphene is a novel 2D carbon allotrope with interesting electronic properties which has been recently theoretically proposed. Phagraphene is similar to a defective graphene structure with an arrangement of pentagonal, heptagonal and hexagonal rings. In this study we investigate the thermal conductivity and mechanical properties of phagraphene using molecular dynamics simulations. Using the non-equilibrium molecular dynamics method, we found the thermal conductivity of phagraphene to be anisotropic, with room temperature values of 218 ± 20 W m−1 K−1 along the armchair direction and 285 ± 29 W m−1 K−1 along the zigzag direction. Both values are one order of magnitude smaller than for pristine graphene. The analysis of the phonon group velocities also shows a significant reduction in this quantity for phagraphene in comparison to graphene. By performing uniaxial tensile simulations, we studied the deformation process and mechanical response of phagraphene. We found that phagraphene exhibits a remarkable high tensile strength around 85 ± 2 GPa, whereas its elastic modulus is also anisotropic along the in-plane directions, with values of 870 ± 15 GPa and 800 ± 14 GPa for the armchair and zigzag directions, respectively. The lower thermal conductivity of phagraphene along with its predicted electronic properties suggests that it could be a better candidate than graphene in future carbon-based thermoelectric devices.
Graphene presents an extraordinary and unique combination of superior physical properties, such as ballistic electronic transport,9,10 high thermal conductivity11–13 and remarkable mechanical strength.14 Among the novel carbon-based materials there is a wealth of physical properties. For example, nitrogenated holey graphene is a semiconductor with a band gap ≈1.96 eV,5 while penta-graphene is predicted to have a 3.25 eV band gap.7 Similar to graphene, phagraphene also presents peculiar electronic properties due to the prediction of distorted Dirac cones. Nevertheless, to the best of our knowledge, the thermal conduction and mechanical properties of phagraphene have not yet been studied theoretically or experimentally. The objective of this work is therefore to evaluate the thermal transport properties and mechanical strength of phagraphene via computer simulations in order to guide future experimental investigations.
In the present work we employ classical molecular dynamics simulations to calculate the thermal conductivity of phagraphene at room temperature, as well as to determine its elastic modulus and mechanical strength. The in-plane thermal conductivity of phagraphene was obtained via non-equilibrium molecular dynamics simulations (NEMD), which include a finite temperature gradient and resemble most experimental setups.12 The mechanical response of phagraphene was investigated by performing uniaxial tensile stress simulations, from which we determine the elastic modulus and the breaking strength of the material. We find that the thermal conductivity of phagraphene is one order of magnitude smaller than that of graphene, while its elastic modulus and tensile strain are comparable to graphene, albeit a little smaller. Remarkably, our simulations show that the thermal conductivity and the elastic modulus of phagraphene are direction dependent, presenting an in-plane anisotropy not observed in the case of graphene and other two-dimensional carbon-based materials.
![]() | ||
Fig. 1 A periodic supercell in the atomic structure of phagraphene. A 20-atom unit cell is also shown. In this study the properties are investigated along the armchair and zigzag directions. |
In order to verify the accuracy of the optimized Tersoff potential for describing the atomic bonding structure of phagraphene we calculated the phonon dispersions via the lattice dynamics software GULP.26,27 Fig. 2 presents the phonon dispersions along high symmetry points of the Brillouin zone, obtained from the 20-atom unit cell shown in Fig. 1. The absence of phonon modes with negative (imaginary) frequencies in the dispersion indicates that the crystal structure of phagraphene is stable when modeled by the chosen potential. In particular, of the three acoustic phonon modes, both in-plane modes present linear dispersions while the flexural mode presents a parabolic dispersion around the Γ point. Furthermore, the presence of high frequency modes (∼50 THz) such as in graphene is related to the high strength of the sp2 chemical bonds between carbon atoms in phagraphene.
Next we performed room temperature molecular dynamics simulations with LAMMPS.28 All simulations employed periodic boundary conditions along both in-plane directions and free boundary conditions in perpendicular directions, such that our simulations would correspond to suspended phagraphene samples. It is important to notice that the presence of 5- and 7-membered rings in the structure of phagraphene is similar to some of the topological defects found in graphene and may result in technical difficulties such as instabilities during simulations at finite temperature.19,21 To avoid this problem, in all molecular dynamics simulations the equations of motion were integrated with a relatively small time increment of 0.25 fs.
NEMD simulations were performed to investigate the thermal conductivity of phagraphene along both in-plane directions, labeled armchair and zigzag in Fig. 1, independently. For each simulation, the system was divided into 22 slabs along the sample length (the direction of heat transport), and atoms at both ends where fixed. The whole structure, excluding fixed atoms, was then thermalized for 100 ps at room temperature (300 K) with a Nosé–Hoover thermostat (NVT).29,30 In our simulations, a temperature gradient is imposed in the system via independent thermostats while the heat flux is calculated from the energy exchanged with the thermostats. Therefore, the original thermostat was turned off, and two independent thermostats were coupled to the “hot” slab (first slab, 310 K) and “cold” slab (22nd slab, 290 K), giving rise to a temperature gradient along the system length. All the atoms between the cold and hot regions evolve freely from thermostats, and the temperature of each slab can be calculated from the average kinetic energy of its atoms via the equipartition theorem. The temperature of slab n is given by:
![]() | (1) |
![]() | (2) |
![]() | (3) |
For the evaluation of phagraphene’s mechanical properties we also applied periodic boundary conditions along both in-plane directions, and the structures were equilibrated at 300 K and zero-pressure via a Nosé–Hoover thermostat and barostat. The loading of the structure was performed via deformation of the simulation box, which was increased along the loading direction by a constant engineering strain rate of 108 s−1 at every simulation time step. In order to guarantee uniaxial stress conditions, a barostat was applied to the direction perpendicular to the loading and set for zero-stress. Virial stresses were calculated at each strain level to obtain the stress–strain curves presented in the next section, and the elastic modulus was obtained directly from the slope of the stress–strain curves.
![]() | ||
Fig. 3 Inverse thermal conductivity versus inverse length. The lines are best fits to the data from which we extract anisotropic thermal conductivities and the effective phonon MFP for phagraphene. |
In general, the size-dependent thermal conductivity is related to the intrinsic thermal conductivity of the material, via:31,35
![]() | (4) |
In order to better understand the difference between thermal conductivities of graphene and phagraphene, we have also calculated the phonon group velocities from the phonon dispersions of each material. The group velocities are shown in Fig. 4 as a function of phonon frequencies, obtained from unit cells with approximately the same dimensions (24 atoms for graphene and 20 atoms for phagraphene). Although the difference in group velocities for the acoustic modes is not so pronounced (notice the modes at 0 THz), for higher frequency modes the group velocities are consistently larger for graphene, as evidenced by the dashed horizontal line. Considering all the phonon frequencies in Fig. 4 we can calculate an average group velocity of 780 m s−1 for graphene, while the equivalent average for phagraphene is only 125 m s−1. This reduction in phonon group velocities is not the only reason for the lower thermal conductivity of phagraphene relative to graphene, but it provides further insight into the origin of the conductivity reduction. Interestingly, the lower thermal conductivity of phagraphene along with its predicted electronic properties suggests that phagraphene could be a better candidate than graphene in future carbon-based thermoelectric devices.8,36,37
The stress–strain response curves of phagraphene along armchair and zigzag directions are presented in Fig. 6, where an anisotropy in the stress–strain curves along the zigzag and armchair directions is apparent. Our simulations predict an elastic modulus of 870 ± 15 GPa along the armchair direction and 800 ± 14 GPa along the zigzag direction. We also predict a remarkable tensile rigidity of 85 ± 2 GPa for phagraphene, which is 64% of the value for pristine graphene and larger than the corresponding values for graphynes.41 Interestingly, our simulations show that the tensile strength of phagraphene is direction independent, as indicated by points a3 and z3 in Fig. 6.
For pristine graphene, the cutoff-modified optimized Tersoff potential predicts that the sample extends uniformly during uniaxial tensile loading and that it remains free of defects up to its tensile strength. At that point, the first debonding occurs between two adjacent carbon atoms. This results in the formation of a crack that rapidly grows and eventually leads to sample rupture. Therefore, for pristine graphene the initial void formation and subsequent failure occur at almost the same strain level, which suggests a brittle failure mechanism.
In Fig. 7, the deformation and breakage process of phagraphene sheets along zigzag and armchair directions is depicted. For phagraphene under tensile loading along the zigzag direction, we observe the initial bond ruptures at strain levels of 0.07 (Fig. 7z1). These debonding events do not happen in pentagon–heptagon pairs but occur in bonds that connect two hexagons, consequently forming ten-membered rings from two carbon hexagons (Fig. 7z1). We observe that at the initial stages of loading, when the sample is stretched along the zigzag direction there are contraction and tension stresses in the sample. In this case we found that the stresses around the pentagon and the bonds connecting the two hexagons are contraction and tension stress, respectively. These diverse stress conditions result in a non-uniform stress distribution and cause stress concentration in specific parts of the samples. Consequently the high tensile stress around the bonds connecting two hexagons leads to local bond ruptures (Fig. 7z1 and z2) forming a more uniform stress distribution. These void initiations can be identified in stress–strain curves at the point in which the slope of the stress–strain curve decreases (point z1 in Fig. 6). Increasing the strain level along the zigzag direction leads to more bond breakage and ten-membered rings extend through the entire sample (Fig. 7z2). After this point, the defect concentration along the sample remains almost unchanged and the formation of new ten-membered rings slows down. In this case, the density of defects is approximately uniform throughout the structure and the sample starts to extend uniformly. This process is distinguishable from the stress–strain response through an increase in the slope of the stress–strain curve (Fig. 7z2 and z2 point in Fig. 6). The ultimate tensile strength happens at the point when the defect coalescence happens (Fig. 7z3) forming a crack that grows and leads to the sample rupture (Fig. 7z4). Meanwhile, our simulations reveal that when phagraphene is under tensile loading along the armchair direction, it remains defect-free up to the high stages of loading (Fig. 7a1). In this case, bonds connecting two hexagons are almost perpendicular to the loading direction and by increasing the strain level their bond length does not change considerably, so the stress distribution is much more uniform in this case in comparison to the zigzag direction. A more uniform stress distribution therefore postpones the first bond breakage in the structure. Nevertheless, at strain levels higher than 0.1, ten-membered rings start to form along the sample, which gradually extend throughout the structure (Fig. 7a2). Similarly to our results along the zigzag direction, the coalescence of ten-membered rings results in the formation of a crack and the subsequent rupture of the sample (Fig. 7a4).
![]() | ||
Fig. 7 Uniaxial deformation process of single-layer phagraphene at various engineering strain levels along the zigzag (a) and armchair (z) loading directions, respectively. The strain and stress levels (z1–z4 and a1–a4) are depicted in Fig. 6. Images generated with OVITO.42 |
Unlike pristine graphene, in phagraphene the initial debonding and subsequent rupture happen at distinctly different strain levels, which suggests ductile failure mechanics for phagraphene. It is worth noticing that, based on our simulations, phagraphene presents a higher ductile failure mechanism along the zigzag direction when compared to its armchair direction. In this case we found that during the tensile deformation and up to sample rupture, phagraphene stretched along the zigzag direction can absorb around 14% more energy than when it is stretched along the armchair direction. In addition, along the zigzag direction phagraphene’s failure strain is found to be 0.15 while this is around 0.13 for the armchair direction.
Table 1 summarizes the predicted thermal transport and mechanical properties of phagraphene, along with the values obtained for pristine graphene with the same empirical potential. While graphene presents anisotropic (in-plane) thermal transport and mechanical properties, phagraphene presents a clear anisotropy for the same properties. Interestingly, while the thermal conductivity and the effective phonon MFP are predicted to be larger along the zigzag direction, the elastic modulus is predicted to be higher along the armchair direction. Nonetheless, while κ along the zigzag direction is 25% larger than along the armchair direction, the elastic modulus along the armchair direction is less than 10% larger than along the zigzag direction. Furthermore, while we are very confident in our predicted anisotropic thermal conductivity, it is possible that the anisotropy observed in tensile strength is an artefact of our simulations, particularly due to the chosen empirical potential which was not explicitly fitted to structures such as phagraphene.
Graphene | Phagraphene | ||
---|---|---|---|
Armchair | Zigzag | ||
κ (W m−1 K−1) | 3050 ± 100 | 218 ± 20 | 285 ± 29 |
Λeff (nm) | 287.3 ± 5.2 | 74.9 ± 5.6 | 94.3 ± 6.1 |
Elast. Mod. (GPa) | 960 ± 20 | 870 ± 15 | 800 ± 14 |
Tens. strength (GPa) | 132 ± 3 | 85 ± 2 | 85 ± 2 |
This journal is © The Royal Society of Chemistry 2016 |