- Email: [email protected]
A molecular dynamics simulation study on thermal conductivity of functionalized bilayer graphene sheet
Chemical Physics Letters 622 (2015) 104–108
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
A molecular dynamics simulation study on thermal conductivity of functionalized bilayer graphene sheet Y.Y. Zhang a,∗ , Q.X. Pei b , X.Q. He c , Y.-W. Mai d a
School of Computing, Engineering and Mathematics, University of Western Sydney, Penrith, NSW 2751, Australia Institute of High Performance Computing, A*STAR, 1 Fusionopolis Way, Singapore 138632, Singapore Department of Civil and Architectural Engineering, City University of Hong Kong, Tat Chee Avenue, Hong Kong d Centre for Advanced Materials Technology, School of Aerospace, Mechanical and Mechatronic Engineering, The University of Sydney, Sydney, NSW 2006, Australia b c
a r t i c l e
i n f o
Article history: Received 5 November 2014 In final form 20 January 2015 Available online 23 January 2015
a b s t r a c t We investigate the thermal conductivity of chemically functionalized bilayer graphene sheet by using molecular dynamics simulations, focusing on the effects of different functional groups and tensile strain. Our simulation results show that chemical functionalization leads to the reduction in the thermal conductivity of bilayer graphene sheet and the reduction is proportional to the coverage of functional groups. Unlike the pristine graphene sheet, the tensile strain can enhance the thermal conductivity of functionalized bilayer graphene sheet. The underlying mechanisms for the thermal conductivity change are elaborated via the vibrational density of states of the graphene sheets. © 2015 Elsevier B.V. All rights reserved.
1. Introduction Graphene is a one-atom-thick planar sheet consisting of sp2 bonded carbon atoms in a honeycomb crystal lattice [1–4]. The strong sp2 -bond has endowed graphene the extremely high in-plane thermal conductivity of 2000–5000 W/mK at room temperature [2]. Therefore, graphene holds great potential as the ideal heat transfer material to be used in the next-generation nanodevices. Enormous research efforts have been devoted to the exploration of the thermal properties of graphene and the effective ways to modulate its thermal conductivity. But compared with the extensive and intensive researches on graphene [5–14], relatively few researches have been conducted on the graphene derivatives, such as graphene oxide (GO). GO is the oxygenated derivative of graphene. It is characteristic of hydroxyl ( OH) and epoxy ( O ) groups on the basal plane while carboxyl ( COOH) and carbonyl ( COH) groups in sheet edges [15,16]. Most recently, for the first time, Lin and Buehler [17] investigated the thermal transport in monolayer GO with randomized surface epoxy and hydroxyl groups at various degrees of oxidation using non-equilibrium molecular dynamics (MD). They showed that the thermal conductivity of GO drops dramatically with respect to the oxidation degree due to the phonon-defect scattering. The suppressed thermal conductivity of
∗ Corresponding author. E-mail address: [email protected] (Y.Y. Zhang). http://dx.doi.org/10.1016/j.cplett.2015.01.034 0009-2614/© 2015 Elsevier B.V. All rights reserved.
the GO suggested that the GO is not a very good heat conductor as compared with the pristine graphene. Almost at the same time, Mu et al. [18] studied the role of oxygen adatoms on the thermal transport in graphene sheet. They found that graphene experiences highly ballistic and amorphous thermal transport with increasing number of oxygen adatoms. The presence of the oxygen adatoms causes dramatic structural deformation and induces significantly enhanced phonon scattering, and thus leads to significant reduction in thermal conductivity. These two available works focus on the monolayer graphene oxide. It is well known that at macroscopic level, GO can form GO paper which is free-standing film consisting of aligned GO nanosheets. These ordered structures have possessed excellent mechanical properties (stiffness up to 40 GPa and tensile strength up to 125 MPa) [19] while remaining highly flexible and ductile due to the interlayer interaction. Among the GO nanosheets, bilayer GO has been shown to have great potential to be used in nano and bio-devices. In addition to the mechanical properties, the thermal behaviour of bilayer GO is also an important property, which remains unexplored. In this work, we carry out large-scale MD simulations to investigate the thermal properties of functionalized bilayer graphene sheet (GS). The effectiveness in suppressing the thermal conductivity of bilayer GS is explored by using different functional groups, such as epoxy, hydroxyl ( OH), hydrogen and the mixture of the first two groups at the coverage up to 20%. The strain engineering is also employed to tune the thermal conductivity of functionalized bilayer graphene sheet.
Y.Y. Zhang et al. / Chemical Physics Letters 622 (2015) 104–108
Figure 1. Schematics of the reverse non-equilibrium molecular dynamics simulation.
2. Simulation model The MD simulations are performed using LAMMPS [20] in which the C/H/O bond interaction is described by ReaxFF [21,22]. The first-principles-based ReaxFF was previously developed for carbon–carbon interactions and hydrocarbon oxidation. It is able to treat large systems with thousands of atoms with near quantum-chemical accuracy. It has been shown to accurately characterize the chemical and mechanical behaviours of hydrocarbons, graphite, diamond, carbon nanotubes and other carbon nanostructures. The thermal conductivity is obtained by using the reverse non-equilibrium molecular dynamics (RNEMD) based on MullerPlathe’s approach [23]. The essence of the RNEMD is to apply a heat flux to the model and determine the temperature gradient induced by this flux. As shown in Figure 1 the heat source and sink slabs are located at the middle and the two ends of the model, respectively. Periodic boundary condition is applied in both the x and y directions. A constant heat flux J is induced along the x direction by continuously exchanging the kinetic energies between the hottest atom in the heat sink slab and the coldest atom in the heat source slab in a microcanonical NVE ensemble (i.e., number of atoms, volume and energy are constant in the system). As the total momentum and the energy of the system are conserved, the heat flux J due to the exchange of atoms is then given. Simultaneously, the temperature profile T(x) with respect to the length of the model is obtained after a time averaging over a prescribed time interval. The thermal conductivity is then calculated by using the Fourier law as = J/2A(∂T/∂x), where ∂T/∂x is the temperature gradient along the heat flux direction, and A is the cross sectional area that the heat flux passes through. The factor 2 in the denominator is used to account for the fact that the heat current flows in two directions as shown in Figure 1. In the present simulation, the initial equilibrium configuration of the functionalized bilayer GS is achieved by the conjugate gradient minimization method. Then the system is equilibrated at 300 K in the NPT ensemble (constant number of atoms, pressure and temperature) for 1 × 106 time steps with each time step of 0.1 fs. The small time step is used to ensure the stability of the simulations and capture the relatively high vibrational frequency of the hydrogen atoms. The environmental temperature is maintained by using the Nose–Hoover thermostat. Thereafter, the heat flux is applied and the system is switched to NVE ensemble. Before the temperature profile is computed, the system is evolved freely for 1 × 107 MD time steps during which the kinetic energies of the hottest atom in the hot slab and the coldest atom in the cold slab are exchanged at every 200 time steps. The temperature gradient along x direction is then obtained by averaging over 5 × 107 time steps. All the MD simulations are performed on the bilayer GSs with a length of 170.4 A˚ and a width of 24.46 A˚ with or without functional groups. The functional groups are randomly distributed on both sides of the GS (see Figure 2) with different coverages varying from 0% to 20%. The coverage of the functional groups is defined as the number of the functional groups divided by the total number of carbon atoms. Each of the H and OH groups bonded with one carbon atom, while the O group bonded with two carbon atoms.
105
Figure 2. Molecular models of bilayer GS with epoxy and hydroxyl groups. Grey balls are C atoms, red balls are O atom, and blue balls are H atoms. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Throughout the Letter, the bilayer GS with hydrogen is denoted as GH, with epoxy as GO, hydroxyl group as GOH, and a mix of 1:1 ratio of epoxy and hydroxyl groups as GMix for easy reference. The effects of different functional groups on the thermal conductivity are assessed by using the pristine bilayer GS as the benchmark. 3. Results and discussions First of all, to validate the accuracy of the simulation methods, we calculate the thermal conductivity of monolayer GSs with ˚ By assuming the thickness of length of 170.4 A˚ and width of 24.6 A. the GS as 0.142 nm [5,24], the thermal conductivity of the monolayer graphene is 135.01 ± 13.38 W/mK and 129.58 ± 8.45 W/mK based on the adaptive intermolecular reactive empirical bond order (AIREBO) [25] potential and ReaxFF force field, respectively. These results are in good agreement with 140 W/mK (for a 20 nm long GS using AIREBO potential) [5] and 114 W/mK (for a 25 nm long GS using ReaxFF force field) [17]. It is clearly shown that the AIREBO and ReaxFF result in close thermal conductivity for the pristine GSs. However, ReaxFF is more computationally expensive but versatile than AIREBO since it is capable to simulate oxygen atoms. In the present work, ReaxFF is employed to simulate the bilayer GS with and without functional groups. It is worth noting that the magnitude of the thermal conductivity is dependent on the selection of the thickness of the GS. In the present work, we limit our attention to the effect of the functional groups and the external strain, rather than the absolute value of the thermal conductivity. Therefore, in the following sections, we use relative thermal conductivity instead, which is the ratio between the thermal conductivities of the functionalized and pristine GSs to avoid the effect of the different thickness assumptions. In addition, the thermal conductivity of monolayer GS (129.58 ± 8.45 W/mK) is higher than that of its bilayer counterpart (120.67 ± 10.89 W/mK) based on ReaxFF force field. It is indicated that the layer number exert negative effect on the thermal conductivity, in consistent with the experimental findings in Ref. [2]. 3.1. Functional groups The relative thermal conductivity of the functionalized bilayer GSs is illustrated in Figure 3. For all the functionalized GS, the relative thermal conductivity is smaller than one, indicating that their thermal conductivity is smaller than that of their pristine counterpart. The presence of the various functional groups leads to the reduction of the thermal conductivity. At the same coverage, the GH has the highest thermal conductivity, followed by the GOH, GMix and GO. The reduction in thermal conductivity induced by the functional groups is attributed to the mass effect as well as the structural deformation. The presence of the functional groups transforms the local sp2 bonding into sp3 bonding, and the associated bond deformations make the GS lose its flat structure. The curvy geometry leads to a change in the local vibrational characteristics
106
Y.Y. Zhang et al. / Chemical Physics Letters 622 (2015) 104–108
Figure 4. VDOS of the GO with different coverages of epoxy groups. Figure 3. Relative thermal conductivity of the bilayer GS with different functional groups with respect to coverage.
and thus an acoustic mismatch that scatters phonon and reduces the phonon mean free path (MFP). As a result, the thermal conductivity is reduced. For the GH and GOH, each hydrogen atom and hydroxyl group are connected to one carbon atom in the GS. However, the hydroxyl group is heavier than hydrogen atom, leading to more mass effects on the phonon and thus a smaller thermal conductivity. The mass effect can also be explained by using the Klemens model [26]. Based on the Klemens model, the phonon scattering rate 1/ due to point defect is given as 1/ ∝ (M/M)2 , where M is the average atomic mass and M is the mass difference of a single substitution defect in a crystal. In the case of C13 isotope defect, M/M = 1/12. Similarly, the M/M of the GH and GOH can be estimated as 1/12 and 17/12. The GOH has larger phonon scattering rate and thus smaller thermal conductivity. In addition, from the initial equilibrium structures, it is found that the presence of hydroxyl group in GOH induces larger out-of-plane deformation at local GS than hydrogen atoms in GH. In combination of these two factors (mass effect and local structural deformation), GOH possesses lower thermal conductivity than GH. However, when half of the hydroxyl groups are replaced by epoxy, GMix possesses even lower thermal conductivity than GOH. This is because the epoxy group is bonded with two carbon atoms in the GS and thus it has more sp3 bonds and more local deformation. Therefore, the presence of the epoxy groups reduces the thermal conductivity further. Likewise, for the GO decorated with epoxy groups only, at the same coverage level, the GO has sp3 bonds twice as many as GOH and 25% more than the GMix, thus the GO has the lowest thermal conductivity. In addition, Figure 3 shows the thermal conductivity of the functionalized bilayer GS decreases with increasing coverage. The functional groups can be regarded as defects. At higher coverage, more phonon-defect scattering occurs, which reduces the overall phonon MFP of the model, and thus leads to the degradation of the thermal conductivity. To gain more understanding on the suppression of thermal conductivity by functionalization, we calculated the phonon spectra or the vibrational density of states (VDOS) for both the pristine and functionalized GS. In the present work, the VDOS of the GS with different functional groups is computed using the autocorrelation function of the atomic velocities of all the atoms. As shown in Figure 4, the presence of epoxy groups in GO has induced a significant reduction in the VDOS of the high frequency phonon modes, and the reduction increases with increasing coverage. This reduction in the density of phonon modes reduces the specific heat of
the phonon modes based on the Debye model for the specific heat of acoustic phonons [27] (specific heat Ci ∝ D(ω) where D(ω) is the VDOS for phonon mode i with frequency ω), and thus reduces the thermal conductivity. In addition, the reduction in the thermal conductivity is also dependent on the coverage, as more functional groups introduce more phonon scattering, which in turn reduce the effective MFP. The calculated VDOS of GSs with different functional groups at the coverage of 5% are shown in Figure 5. As stated in the previous section, the presence of functional groups activates more phonon scattering and reduces the VDOS of C C bonds (see zoomed part in Figure 5b), which lead to the reduction of the effective phonon MFP and the specific heat, respectively. Therefore, pristine GS possesses the highest thermal conductivity. For the VDOS of the different functional groups in Figure 5a, there are three distinct peaks at 40 THz, 57 THz and 100 THz, reflecting the C H (in GH), C C and O H (in GOH and GMix) bond interactions, respectively. Based on the Debye model [27], the much higher VDOS at 40 THz and slightly higher magnitude of VDOS at 57 THz for GH lead to higher specific heat for GH. Although GOH and GMix have more phonon modes at 100 THz due to the O H bond interactions, the VDOS are much smaller compared to C H phonon mode at 40 THz. Therefore, the thermal conductivity of the functionalized GS exhibits the trend of GH > GOH > GMix > GO. In order to further elucidate the effect of the functional groups on the suppression of the thermal conductivity, we plot the VDOS of top and bottom layers in the bi-layer pristine GS and GO in Figure 6a and b, respectively. In the pristine GS, the VDOS of the two layers match well with each other. The good overlap in VDOS means that phonons can easily pass through the system, resulting in a high thermal conductivity. But in the presence of the functional groups and the associated formation of out-of-plane sp3 bonds in GO, the VDOS of the two layers exhibit obvious mismatch in the low (<40 TPa) and high (>80 TPa) frequency regions, which consequently lowers the thermal conductivity of GO. 3.2. Tensile strain Strain engineering is widely used in the manipulation of the thermal properties of nanomaterials [28,29]. In the following, we will investigate how the thermal conductivity responds to the external tensile strain by performing simulations on GMix models with 1% coverage for illustration purpose. Figure 7 shows the relative thermal conductivity of the GMix with respect to the tensile strain up to 0.1. All the relative thermal conductivity is higher than 1, meaning that the GMix has larger thermal conductivity under
Y.Y. Zhang et al. / Chemical Physics Letters 622 (2015) 104–108
107
Figure 5. (a) VDOS of the bi-layer GS with different functional groups at 5% coverage; (b) Zoom of (a) in the frequency range of 50–65 THz.
elongation. The thermal conductivity increases proportionally with increasing tensile strain and it reaches to a peak value at a strain of 0.04, thereafter, the thermal conductivity decreases with increasing tensile strain up to 0.1. It was well known that the thermal conductivity of pure graphene decreases monotonously with increasing tensile strain [5,6]. Therefore, this increasing-and-then-decreasing behaviour of thermal conductivity of GMix is in contrast with that of the pure graphene. In order to explore the underlying mechanism for the strain effect on the thermal conductivity of GMix, the VDOS of the GMix under different strain was calculated and is shown in Figure 8. In order to have a clearer view of the in-plane phonon mode, Figure 8 shows the VDOS in the frequency range from 40 to 80 THz. It is seen that the tensile strain leads to very obvious phonon softening (left shift of phonon modes) when the strain is larger than 0.04. This phonon softening reduces the phonon group velocity and thus reduces the thermal conductivity [5]. Therefore, the decreasing trend of thermal conductivity of GMix for the tensile strain beyond 0.04 can be attributed to the phonon softening-induced lower thermal conductivity. For the increasing trend of thermal conductivity during the small strain region from 0 to 0.04 shown in Figure 7, it can be explained as follows. The functional groups on graphene lead to the local C C bonding transition from in-plane sp2 bonds to outof-plane sp3 bonds, which results in wrinkles in the graphene sheet [30]. These wrinkles increase the phonon scattering and reduce the thermal conductivity. At the small strain region, the increase of the
Figure 6. VDOS at top and bottom layers for (a) pristine bilayer GS and (b) GO with 5% coverage.
Figure 7. Relative thermal conductivity of GMix with 1% coverage with respect to tensile strain.
108
Y.Y. Zhang et al. / Chemical Physics Letters 622 (2015) 104–108
pristine counterpart. We have elaborated the underlying mechanism for the thermal conductivity change induced by the functional groups and the tensile stain by the vibrational density of states. Acknowledgement The first author is grateful for the computational support provided by Intersect Australia Ltd. References
Figure 8. VDOS of the GMix under different tensile strain.
tensile strain reduces the wrinkles and so increases the thermal conductivity. 4. Conclusions We performed non-equilibrium MD simulations based on ReaxFF force field to investigate the thermal conductivity of bilayer GS functionalized with various functional groups. Our simulation results reveal that the thermal conductivity of bilayer GS can be tuned by chemical functionalization and tensile strain to meet different thermal applications. The thermal conductivity of the GS can be reduced significantly by using different functional groups, which is beneficial for its thermoelectronic applications requirements. In particular, the GS functionalized by the pure epoxy (GO) exhibits the lowest thermal conductivity and thus is the most suitable candidature in these applications. For the thermal applications that require efficient heat conduction, the thermal conductivity of the functionalized GS can be enhanced via strain engineering unlike its
[1] K.S. Novoselov, et al., Science 306 (2004) 666. [2] S. Ghosh, W.Z. Bao, D.L. Nika, S. Subrina, E.P. Pokatilov, C.N. Lau, A.A. Balandin, Nat. Mater. 9 (2010) 555. [3] A.K. Geim, K.S. Novoselov, Nat. Mater. 6 (2007) 183. [4] C. Lee, X. Wei, J.W. Kysar, J. Hone, Science 321 (2008) 385. [5] N. Wei, L. Xu, H.Q. Wang, J.C. Zheng, Nanotechnology 22 (2011) 105705. [6] L. Lindsay, W. Li, J. Carrete, N. Mingo, D.A. Broido, T.L. Reinecke, Phys. Rev. B 89 (2014) 155426. [7] J. Huang, C.H. Wong, Comput. Mater. Sci. 92 (2014) 192. [8] Q.X. Pei, Z.D. Sha, Y.W. Zhang, Carbon 49 (2011) 4752. [9] W.X. Huang, Q.X. Pei, Z.S. Liu, Y.W. Zhang, Chem. Phys. Lett. 552 (2012) 97. [10] Y.Y. Zhang, Q.X. Pei, C.M. Wang, Comput. Mater. Sci. 65 (2012) 406. [11] D. Konatham, D.V. Papavassiliou, A. Striolo, Chem. Phys. Lett. 527 (2012) 47. [12] T. Ouyang, Y. Chen, Y. Xie, G.M. Stocks, J. Zhong, Appl. Phys. Lett. 99 (2011) 233101. [13] T.Y. Ng, J.J. Yeo, Z.S. Liu, Carbon 50 (2012) 4887. [14] U. Ray, G. Balasubramanian, Chem. Phys. Lett. 599 (2014) 154. [15] W. Gao, L.B. Alemany, L. Ci, P.M. Ajayan, Nat. Chem. 1 (2009) 403. [16] A. Bagri, C. Mattevi, A. Acik, Y.J. Chabal, M. Chhowalla, V.B. Shenoy, Nat. Chem. 2 (2010) 581. [17] S.C. Lin, M.J. Buehler, Carbon 77 (2014) 351. [18] X. Mu, X.F. Wu, T. Zhang, D.B. Go, T.F. Luo, Sci. Rep. 4 (2014) 3909. [19] O.C. Compton, S.W. Cranford, K.W. Putz, Z. An, L.C. Brinson, M.J. Buehler, S.T. Nguyen, ACS Nano 6 (2012) 2008. [20] S. Plimpton, J. Comput. Phys. 117 (1995) 1. [21] K. Chenoweth, A.C.T. van Duin, W.A. Goddard, J. Phys. Chem. A 112 (2008) 1040. [22] H.M. Aktulga, J.C. Fogarty, S.A. Pandit, A.Y. Grama, Parallel Comput. 38 (2012) 245. [23] F. Muller-Plathe, J. Chem. Phys. 106 (1997) 6082. [24] Z. Guo, D. Zhang, X.G. Gong, Appl. Phys. Lett. 95 (2009) 163103. [25] S.J. Stuart, A.B. Tutein, J.A. Harrison, J. Chem. Phys. 112 (2000) 6472. [26] C.A. Ratsifaritana, P.G. Klemens, Int. J. Thermophys. 8 (1987) 737. [27] G. Chen, Nanoscale Energy Transport and Conversion: A Parallel Treatment of Electrons, Molecules, Phonons, and Photons, Oxford University Press, Oxford, New York, 2005. [28] Q.X. Pei, Y.W. Zhang, Z.D. Sha, V.B. Shenoy, J. Appl. Phys. 114 (2013) 033526. [29] M. Hu, X. Zhang, D. Poulikakos, Phys. Rev. B 87 (2013) 195417. [30] Q.X. Pei, Y.W. Zhang, V.B. Shenoy, Carbon 48 (2010) 898.
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
A molecular dynamics simulation study on thermal conductivity of functionalized bilayer graphene sheet Y.Y. Zhang a,∗ , Q.X. Pei b , X.Q. He c , Y.-W. Mai d a
School of Computing, Engineering and Mathematics, University of Western Sydney, Penrith, NSW 2751, Australia Institute of High Performance Computing, A*STAR, 1 Fusionopolis Way, Singapore 138632, Singapore Department of Civil and Architectural Engineering, City University of Hong Kong, Tat Chee Avenue, Hong Kong d Centre for Advanced Materials Technology, School of Aerospace, Mechanical and Mechatronic Engineering, The University of Sydney, Sydney, NSW 2006, Australia b c
a r t i c l e
i n f o
Article history: Received 5 November 2014 In final form 20 January 2015 Available online 23 January 2015
a b s t r a c t We investigate the thermal conductivity of chemically functionalized bilayer graphene sheet by using molecular dynamics simulations, focusing on the effects of different functional groups and tensile strain. Our simulation results show that chemical functionalization leads to the reduction in the thermal conductivity of bilayer graphene sheet and the reduction is proportional to the coverage of functional groups. Unlike the pristine graphene sheet, the tensile strain can enhance the thermal conductivity of functionalized bilayer graphene sheet. The underlying mechanisms for the thermal conductivity change are elaborated via the vibrational density of states of the graphene sheets. © 2015 Elsevier B.V. All rights reserved.
1. Introduction Graphene is a one-atom-thick planar sheet consisting of sp2 bonded carbon atoms in a honeycomb crystal lattice [1–4]. The strong sp2 -bond has endowed graphene the extremely high in-plane thermal conductivity of 2000–5000 W/mK at room temperature [2]. Therefore, graphene holds great potential as the ideal heat transfer material to be used in the next-generation nanodevices. Enormous research efforts have been devoted to the exploration of the thermal properties of graphene and the effective ways to modulate its thermal conductivity. But compared with the extensive and intensive researches on graphene [5–14], relatively few researches have been conducted on the graphene derivatives, such as graphene oxide (GO). GO is the oxygenated derivative of graphene. It is characteristic of hydroxyl ( OH) and epoxy ( O ) groups on the basal plane while carboxyl ( COOH) and carbonyl ( COH) groups in sheet edges [15,16]. Most recently, for the first time, Lin and Buehler [17] investigated the thermal transport in monolayer GO with randomized surface epoxy and hydroxyl groups at various degrees of oxidation using non-equilibrium molecular dynamics (MD). They showed that the thermal conductivity of GO drops dramatically with respect to the oxidation degree due to the phonon-defect scattering. The suppressed thermal conductivity of
∗ Corresponding author. E-mail address: [email protected] (Y.Y. Zhang). http://dx.doi.org/10.1016/j.cplett.2015.01.034 0009-2614/© 2015 Elsevier B.V. All rights reserved.
the GO suggested that the GO is not a very good heat conductor as compared with the pristine graphene. Almost at the same time, Mu et al. [18] studied the role of oxygen adatoms on the thermal transport in graphene sheet. They found that graphene experiences highly ballistic and amorphous thermal transport with increasing number of oxygen adatoms. The presence of the oxygen adatoms causes dramatic structural deformation and induces significantly enhanced phonon scattering, and thus leads to significant reduction in thermal conductivity. These two available works focus on the monolayer graphene oxide. It is well known that at macroscopic level, GO can form GO paper which is free-standing film consisting of aligned GO nanosheets. These ordered structures have possessed excellent mechanical properties (stiffness up to 40 GPa and tensile strength up to 125 MPa) [19] while remaining highly flexible and ductile due to the interlayer interaction. Among the GO nanosheets, bilayer GO has been shown to have great potential to be used in nano and bio-devices. In addition to the mechanical properties, the thermal behaviour of bilayer GO is also an important property, which remains unexplored. In this work, we carry out large-scale MD simulations to investigate the thermal properties of functionalized bilayer graphene sheet (GS). The effectiveness in suppressing the thermal conductivity of bilayer GS is explored by using different functional groups, such as epoxy, hydroxyl ( OH), hydrogen and the mixture of the first two groups at the coverage up to 20%. The strain engineering is also employed to tune the thermal conductivity of functionalized bilayer graphene sheet.
Y.Y. Zhang et al. / Chemical Physics Letters 622 (2015) 104–108
Figure 1. Schematics of the reverse non-equilibrium molecular dynamics simulation.
2. Simulation model The MD simulations are performed using LAMMPS [20] in which the C/H/O bond interaction is described by ReaxFF [21,22]. The first-principles-based ReaxFF was previously developed for carbon–carbon interactions and hydrocarbon oxidation. It is able to treat large systems with thousands of atoms with near quantum-chemical accuracy. It has been shown to accurately characterize the chemical and mechanical behaviours of hydrocarbons, graphite, diamond, carbon nanotubes and other carbon nanostructures. The thermal conductivity is obtained by using the reverse non-equilibrium molecular dynamics (RNEMD) based on MullerPlathe’s approach [23]. The essence of the RNEMD is to apply a heat flux to the model and determine the temperature gradient induced by this flux. As shown in Figure 1 the heat source and sink slabs are located at the middle and the two ends of the model, respectively. Periodic boundary condition is applied in both the x and y directions. A constant heat flux J is induced along the x direction by continuously exchanging the kinetic energies between the hottest atom in the heat sink slab and the coldest atom in the heat source slab in a microcanonical NVE ensemble (i.e., number of atoms, volume and energy are constant in the system). As the total momentum and the energy of the system are conserved, the heat flux J due to the exchange of atoms is then given. Simultaneously, the temperature profile T(x) with respect to the length of the model is obtained after a time averaging over a prescribed time interval. The thermal conductivity is then calculated by using the Fourier law as = J/2A(∂T/∂x), where ∂T/∂x is the temperature gradient along the heat flux direction, and A is the cross sectional area that the heat flux passes through. The factor 2 in the denominator is used to account for the fact that the heat current flows in two directions as shown in Figure 1. In the present simulation, the initial equilibrium configuration of the functionalized bilayer GS is achieved by the conjugate gradient minimization method. Then the system is equilibrated at 300 K in the NPT ensemble (constant number of atoms, pressure and temperature) for 1 × 106 time steps with each time step of 0.1 fs. The small time step is used to ensure the stability of the simulations and capture the relatively high vibrational frequency of the hydrogen atoms. The environmental temperature is maintained by using the Nose–Hoover thermostat. Thereafter, the heat flux is applied and the system is switched to NVE ensemble. Before the temperature profile is computed, the system is evolved freely for 1 × 107 MD time steps during which the kinetic energies of the hottest atom in the hot slab and the coldest atom in the cold slab are exchanged at every 200 time steps. The temperature gradient along x direction is then obtained by averaging over 5 × 107 time steps. All the MD simulations are performed on the bilayer GSs with a length of 170.4 A˚ and a width of 24.46 A˚ with or without functional groups. The functional groups are randomly distributed on both sides of the GS (see Figure 2) with different coverages varying from 0% to 20%. The coverage of the functional groups is defined as the number of the functional groups divided by the total number of carbon atoms. Each of the H and OH groups bonded with one carbon atom, while the O group bonded with two carbon atoms.
105
Figure 2. Molecular models of bilayer GS with epoxy and hydroxyl groups. Grey balls are C atoms, red balls are O atom, and blue balls are H atoms. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Throughout the Letter, the bilayer GS with hydrogen is denoted as GH, with epoxy as GO, hydroxyl group as GOH, and a mix of 1:1 ratio of epoxy and hydroxyl groups as GMix for easy reference. The effects of different functional groups on the thermal conductivity are assessed by using the pristine bilayer GS as the benchmark. 3. Results and discussions First of all, to validate the accuracy of the simulation methods, we calculate the thermal conductivity of monolayer GSs with ˚ By assuming the thickness of length of 170.4 A˚ and width of 24.6 A. the GS as 0.142 nm [5,24], the thermal conductivity of the monolayer graphene is 135.01 ± 13.38 W/mK and 129.58 ± 8.45 W/mK based on the adaptive intermolecular reactive empirical bond order (AIREBO) [25] potential and ReaxFF force field, respectively. These results are in good agreement with 140 W/mK (for a 20 nm long GS using AIREBO potential) [5] and 114 W/mK (for a 25 nm long GS using ReaxFF force field) [17]. It is clearly shown that the AIREBO and ReaxFF result in close thermal conductivity for the pristine GSs. However, ReaxFF is more computationally expensive but versatile than AIREBO since it is capable to simulate oxygen atoms. In the present work, ReaxFF is employed to simulate the bilayer GS with and without functional groups. It is worth noting that the magnitude of the thermal conductivity is dependent on the selection of the thickness of the GS. In the present work, we limit our attention to the effect of the functional groups and the external strain, rather than the absolute value of the thermal conductivity. Therefore, in the following sections, we use relative thermal conductivity instead, which is the ratio between the thermal conductivities of the functionalized and pristine GSs to avoid the effect of the different thickness assumptions. In addition, the thermal conductivity of monolayer GS (129.58 ± 8.45 W/mK) is higher than that of its bilayer counterpart (120.67 ± 10.89 W/mK) based on ReaxFF force field. It is indicated that the layer number exert negative effect on the thermal conductivity, in consistent with the experimental findings in Ref. [2]. 3.1. Functional groups The relative thermal conductivity of the functionalized bilayer GSs is illustrated in Figure 3. For all the functionalized GS, the relative thermal conductivity is smaller than one, indicating that their thermal conductivity is smaller than that of their pristine counterpart. The presence of the various functional groups leads to the reduction of the thermal conductivity. At the same coverage, the GH has the highest thermal conductivity, followed by the GOH, GMix and GO. The reduction in thermal conductivity induced by the functional groups is attributed to the mass effect as well as the structural deformation. The presence of the functional groups transforms the local sp2 bonding into sp3 bonding, and the associated bond deformations make the GS lose its flat structure. The curvy geometry leads to a change in the local vibrational characteristics
106
Y.Y. Zhang et al. / Chemical Physics Letters 622 (2015) 104–108
Figure 4. VDOS of the GO with different coverages of epoxy groups. Figure 3. Relative thermal conductivity of the bilayer GS with different functional groups with respect to coverage.
and thus an acoustic mismatch that scatters phonon and reduces the phonon mean free path (MFP). As a result, the thermal conductivity is reduced. For the GH and GOH, each hydrogen atom and hydroxyl group are connected to one carbon atom in the GS. However, the hydroxyl group is heavier than hydrogen atom, leading to more mass effects on the phonon and thus a smaller thermal conductivity. The mass effect can also be explained by using the Klemens model [26]. Based on the Klemens model, the phonon scattering rate 1/ due to point defect is given as 1/ ∝ (M/M)2 , where M is the average atomic mass and M is the mass difference of a single substitution defect in a crystal. In the case of C13 isotope defect, M/M = 1/12. Similarly, the M/M of the GH and GOH can be estimated as 1/12 and 17/12. The GOH has larger phonon scattering rate and thus smaller thermal conductivity. In addition, from the initial equilibrium structures, it is found that the presence of hydroxyl group in GOH induces larger out-of-plane deformation at local GS than hydrogen atoms in GH. In combination of these two factors (mass effect and local structural deformation), GOH possesses lower thermal conductivity than GH. However, when half of the hydroxyl groups are replaced by epoxy, GMix possesses even lower thermal conductivity than GOH. This is because the epoxy group is bonded with two carbon atoms in the GS and thus it has more sp3 bonds and more local deformation. Therefore, the presence of the epoxy groups reduces the thermal conductivity further. Likewise, for the GO decorated with epoxy groups only, at the same coverage level, the GO has sp3 bonds twice as many as GOH and 25% more than the GMix, thus the GO has the lowest thermal conductivity. In addition, Figure 3 shows the thermal conductivity of the functionalized bilayer GS decreases with increasing coverage. The functional groups can be regarded as defects. At higher coverage, more phonon-defect scattering occurs, which reduces the overall phonon MFP of the model, and thus leads to the degradation of the thermal conductivity. To gain more understanding on the suppression of thermal conductivity by functionalization, we calculated the phonon spectra or the vibrational density of states (VDOS) for both the pristine and functionalized GS. In the present work, the VDOS of the GS with different functional groups is computed using the autocorrelation function of the atomic velocities of all the atoms. As shown in Figure 4, the presence of epoxy groups in GO has induced a significant reduction in the VDOS of the high frequency phonon modes, and the reduction increases with increasing coverage. This reduction in the density of phonon modes reduces the specific heat of
the phonon modes based on the Debye model for the specific heat of acoustic phonons [27] (specific heat Ci ∝ D(ω) where D(ω) is the VDOS for phonon mode i with frequency ω), and thus reduces the thermal conductivity. In addition, the reduction in the thermal conductivity is also dependent on the coverage, as more functional groups introduce more phonon scattering, which in turn reduce the effective MFP. The calculated VDOS of GSs with different functional groups at the coverage of 5% are shown in Figure 5. As stated in the previous section, the presence of functional groups activates more phonon scattering and reduces the VDOS of C C bonds (see zoomed part in Figure 5b), which lead to the reduction of the effective phonon MFP and the specific heat, respectively. Therefore, pristine GS possesses the highest thermal conductivity. For the VDOS of the different functional groups in Figure 5a, there are three distinct peaks at 40 THz, 57 THz and 100 THz, reflecting the C H (in GH), C C and O H (in GOH and GMix) bond interactions, respectively. Based on the Debye model [27], the much higher VDOS at 40 THz and slightly higher magnitude of VDOS at 57 THz for GH lead to higher specific heat for GH. Although GOH and GMix have more phonon modes at 100 THz due to the O H bond interactions, the VDOS are much smaller compared to C H phonon mode at 40 THz. Therefore, the thermal conductivity of the functionalized GS exhibits the trend of GH > GOH > GMix > GO. In order to further elucidate the effect of the functional groups on the suppression of the thermal conductivity, we plot the VDOS of top and bottom layers in the bi-layer pristine GS and GO in Figure 6a and b, respectively. In the pristine GS, the VDOS of the two layers match well with each other. The good overlap in VDOS means that phonons can easily pass through the system, resulting in a high thermal conductivity. But in the presence of the functional groups and the associated formation of out-of-plane sp3 bonds in GO, the VDOS of the two layers exhibit obvious mismatch in the low (<40 TPa) and high (>80 TPa) frequency regions, which consequently lowers the thermal conductivity of GO. 3.2. Tensile strain Strain engineering is widely used in the manipulation of the thermal properties of nanomaterials [28,29]. In the following, we will investigate how the thermal conductivity responds to the external tensile strain by performing simulations on GMix models with 1% coverage for illustration purpose. Figure 7 shows the relative thermal conductivity of the GMix with respect to the tensile strain up to 0.1. All the relative thermal conductivity is higher than 1, meaning that the GMix has larger thermal conductivity under
Y.Y. Zhang et al. / Chemical Physics Letters 622 (2015) 104–108
107
Figure 5. (a) VDOS of the bi-layer GS with different functional groups at 5% coverage; (b) Zoom of (a) in the frequency range of 50–65 THz.
elongation. The thermal conductivity increases proportionally with increasing tensile strain and it reaches to a peak value at a strain of 0.04, thereafter, the thermal conductivity decreases with increasing tensile strain up to 0.1. It was well known that the thermal conductivity of pure graphene decreases monotonously with increasing tensile strain [5,6]. Therefore, this increasing-and-then-decreasing behaviour of thermal conductivity of GMix is in contrast with that of the pure graphene. In order to explore the underlying mechanism for the strain effect on the thermal conductivity of GMix, the VDOS of the GMix under different strain was calculated and is shown in Figure 8. In order to have a clearer view of the in-plane phonon mode, Figure 8 shows the VDOS in the frequency range from 40 to 80 THz. It is seen that the tensile strain leads to very obvious phonon softening (left shift of phonon modes) when the strain is larger than 0.04. This phonon softening reduces the phonon group velocity and thus reduces the thermal conductivity [5]. Therefore, the decreasing trend of thermal conductivity of GMix for the tensile strain beyond 0.04 can be attributed to the phonon softening-induced lower thermal conductivity. For the increasing trend of thermal conductivity during the small strain region from 0 to 0.04 shown in Figure 7, it can be explained as follows. The functional groups on graphene lead to the local C C bonding transition from in-plane sp2 bonds to outof-plane sp3 bonds, which results in wrinkles in the graphene sheet [30]. These wrinkles increase the phonon scattering and reduce the thermal conductivity. At the small strain region, the increase of the
Figure 6. VDOS at top and bottom layers for (a) pristine bilayer GS and (b) GO with 5% coverage.
Figure 7. Relative thermal conductivity of GMix with 1% coverage with respect to tensile strain.
108
Y.Y. Zhang et al. / Chemical Physics Letters 622 (2015) 104–108
pristine counterpart. We have elaborated the underlying mechanism for the thermal conductivity change induced by the functional groups and the tensile stain by the vibrational density of states. Acknowledgement The first author is grateful for the computational support provided by Intersect Australia Ltd. References
Figure 8. VDOS of the GMix under different tensile strain.
tensile strain reduces the wrinkles and so increases the thermal conductivity. 4. Conclusions We performed non-equilibrium MD simulations based on ReaxFF force field to investigate the thermal conductivity of bilayer GS functionalized with various functional groups. Our simulation results reveal that the thermal conductivity of bilayer GS can be tuned by chemical functionalization and tensile strain to meet different thermal applications. The thermal conductivity of the GS can be reduced significantly by using different functional groups, which is beneficial for its thermoelectronic applications requirements. In particular, the GS functionalized by the pure epoxy (GO) exhibits the lowest thermal conductivity and thus is the most suitable candidature in these applications. For the thermal applications that require efficient heat conduction, the thermal conductivity of the functionalized GS can be enhanced via strain engineering unlike its
[1] K.S. Novoselov, et al., Science 306 (2004) 666. [2] S. Ghosh, W.Z. Bao, D.L. Nika, S. Subrina, E.P. Pokatilov, C.N. Lau, A.A. Balandin, Nat. Mater. 9 (2010) 555. [3] A.K. Geim, K.S. Novoselov, Nat. Mater. 6 (2007) 183. [4] C. Lee, X. Wei, J.W. Kysar, J. Hone, Science 321 (2008) 385. [5] N. Wei, L. Xu, H.Q. Wang, J.C. Zheng, Nanotechnology 22 (2011) 105705. [6] L. Lindsay, W. Li, J. Carrete, N. Mingo, D.A. Broido, T.L. Reinecke, Phys. Rev. B 89 (2014) 155426. [7] J. Huang, C.H. Wong, Comput. Mater. Sci. 92 (2014) 192. [8] Q.X. Pei, Z.D. Sha, Y.W. Zhang, Carbon 49 (2011) 4752. [9] W.X. Huang, Q.X. Pei, Z.S. Liu, Y.W. Zhang, Chem. Phys. Lett. 552 (2012) 97. [10] Y.Y. Zhang, Q.X. Pei, C.M. Wang, Comput. Mater. Sci. 65 (2012) 406. [11] D. Konatham, D.V. Papavassiliou, A. Striolo, Chem. Phys. Lett. 527 (2012) 47. [12] T. Ouyang, Y. Chen, Y. Xie, G.M. Stocks, J. Zhong, Appl. Phys. Lett. 99 (2011) 233101. [13] T.Y. Ng, J.J. Yeo, Z.S. Liu, Carbon 50 (2012) 4887. [14] U. Ray, G. Balasubramanian, Chem. Phys. Lett. 599 (2014) 154. [15] W. Gao, L.B. Alemany, L. Ci, P.M. Ajayan, Nat. Chem. 1 (2009) 403. [16] A. Bagri, C. Mattevi, A. Acik, Y.J. Chabal, M. Chhowalla, V.B. Shenoy, Nat. Chem. 2 (2010) 581. [17] S.C. Lin, M.J. Buehler, Carbon 77 (2014) 351. [18] X. Mu, X.F. Wu, T. Zhang, D.B. Go, T.F. Luo, Sci. Rep. 4 (2014) 3909. [19] O.C. Compton, S.W. Cranford, K.W. Putz, Z. An, L.C. Brinson, M.J. Buehler, S.T. Nguyen, ACS Nano 6 (2012) 2008. [20] S. Plimpton, J. Comput. Phys. 117 (1995) 1. [21] K. Chenoweth, A.C.T. van Duin, W.A. Goddard, J. Phys. Chem. A 112 (2008) 1040. [22] H.M. Aktulga, J.C. Fogarty, S.A. Pandit, A.Y. Grama, Parallel Comput. 38 (2012) 245. [23] F. Muller-Plathe, J. Chem. Phys. 106 (1997) 6082. [24] Z. Guo, D. Zhang, X.G. Gong, Appl. Phys. Lett. 95 (2009) 163103. [25] S.J. Stuart, A.B. Tutein, J.A. Harrison, J. Chem. Phys. 112 (2000) 6472. [26] C.A. Ratsifaritana, P.G. Klemens, Int. J. Thermophys. 8 (1987) 737. [27] G. Chen, Nanoscale Energy Transport and Conversion: A Parallel Treatment of Electrons, Molecules, Phonons, and Photons, Oxford University Press, Oxford, New York, 2005. [28] Q.X. Pei, Y.W. Zhang, Z.D. Sha, V.B. Shenoy, J. Appl. Phys. 114 (2013) 033526. [29] M. Hu, X. Zhang, D. Poulikakos, Phys. Rev. B 87 (2013) 195417. [30] Q.X. Pei, Y.W. Zhang, V.B. Shenoy, Carbon 48 (2010) 898.