两种二维碳同素异形体的热导率：分子动力学研究

介绍

碳材料，例如金刚石 [1]、碳纳米管 [2,3,4,5] 和石墨烯 [6,7,8,9,10,11,12]，由于其优异的热传输特性。特别是低维碳材料在传热方面表现出优异的性能。作为一维材料，单个碳纳米管的高热导率 (TC) 已通过实验 [2, 3] 和理论研究 [4, 5] 观察到。此外，作为一种单原子厚的扁平二维（2D）碳材料，石墨烯因其高 TC 被认为是下一代导热增强复合材料的革命性材料 [6,7,8,9,10 ,11,12]。也有报道称，石墨烯的TC可以达到石墨烯的40%，在热管理方面具有潜在的应用[13,14,15]。

受这些碳同素异形体迷人特性的启发，近年来研究人员对碳同素异形体及其衍生物进行了大量研究。已经采用实验和理论方法研究了新型二维碳同素异形体，例如 sp ² - 具有五元、六元和七元环的类碳层 [16]；具有四元环的二维无定形碳 [17]；平面碳五肽[18]；具有图案化缺陷的二维碳半导体 [19]；几个二维平面碳网络[20]；八石墨烯 [21]; T-石墨烯[22]；和 H-net [23]。确定这些二维碳同素异形体的独特性质对于未来电子、光子和热领域的纳米材料具有重要意义[16,17,18,19,20,21,22,23]。

随着人们对探索二维碳同素异形体的新结构越来越感兴趣，Su 等人。 [24]通过第一性原理计算提出了两种新的能量竞争性和动力学稳定的二维碳同素异形体，由八边形和五边形组成。这两种碳片的动力学稳定性通过计算它们的声子色散曲线得到证实。由于这两种碳同素异形体的结构可以看作是沿直线路径和锯齿形路径复制五五八元环（558）带，因此这两种碳同素异形体被命名为八角形和五边形石墨烯线（OPG-L）和八边形和五边形石墨烯之字形（OPG-Z），分别。这两种碳同素异形体的形成能分别为 0.31 eV/atom 和 0.34 eV/atom。这些值远低于先前合成的石墨炔的形成能，即 0.76 eV/原子 [25]。值得注意的是，OPG-Z 具有显着的电子结构各向异性，在电子器件中具有潜在的应用 [24]。因此，为了满足 OPG-L 和 OPG-Z 的电子应用要求，研究这两种新型结构的散热性能是必然和必要的。目前，这两种结构的热学性质仍不清楚。

在这项工作中，我们使用分子动力学模拟研究了两种新型二维碳同素异形体的热性质。探讨了尺寸、应变和温度对 TC 的影响。通过计算声子的振动态密度 (VDOS) 来分析结果。我们对这两种碳同素异形体热性能的研究表明了它们在热管理器件中的潜在应用。

模型和方法

OPG-L（图 1a）和 OPG-Z（图 1b）的结构包含由八边形和五边形组成的代表性细胞 [24]。为了区分结构的边缘类型，我们定义了扶手椅和锯齿形的手性，就像石墨烯一样（见图 1）。这两种结构可以由红色原子表示的代表性558条带使用沿绿色行的平移对称性形成。

a的示意图模型 OPG-L 和 b OPG-Z。黑色虚线框是 OPG-L 和 OPG-Z 的正交晶胞，其中 OA 和 OB 是晶格向量。 OPG-L的原胞以蓝色虚线框表示，OPG-Z的原胞与晶胞相同

所有 MD 模拟都是使用大规模原子/分子大规模并行模拟器 (LAMMPS) 包进行的 [26]。我们使用 Lindsay 和 Broido [27] 优化的 Tersoff 势，稍加修改，即修改后的优化 Tersoff 势，来描述碳原子之间的相互作用。与原始 Tersoff 势 [28] 相比，Lindsay 和 Brodio 优化了两个参数，一个用于平衡键角，另一个用于吸引力相互作用强度。根据此优化的 Tersoff 电位 [27]，石墨烯中的平衡键长为 1.4388 Å，大于实验值 1.42 Å [29]。因为 Tersoff 势中唯一与长度相关的参数是 λ 1 在排斥函数 (f ^R =A exp.(-λ 1r )) 和 λ 2 在吸引函数 (f ^A =B exp(-λ 2r ))，我们可以通过将这两个参数乘以因子 1.4388/1.42 来获得正确的键长。也就是说，我们改变λ 1 来自 3.4879 Å ⁻¹ 到 3.5333 Å ⁻¹ 并改变 λ 2 来自 2.2119 Å ⁻¹ 到 2.2407 Å ^{− 1} .这些修改仅以全局方式改变势能的长度尺度。基于这种改进的优化 Tersoff 势，MD 模拟中相应的平衡晶格参数如下：OA =3.63 Å，OPG-L 中的 OB =9.38 Å 和 OPG-L 中的 OA =6.78 Å，OB =5.04Z，O 中的 Å与 Su 等人先前的研究非常一致。 [24]，即 OPG-L 中的 OA =3.68 Å，OB =9.12 Å，OPG-Z 中的 OA =6.90 Å，OB =4.87 Å。

执行反向非平衡分子动力学 (rNEMD) [30] 模拟以计算 TC。在 x 和 y 维度上采用周期性边界条件。 OPG-L 和 OPG-Z 的结构最初通过共轭梯度算法的 Polak-Ribiered 版本进行优化 [31]，然后使用 0.25-ns Nosé-Hoover 热浴 [32, 33] 以确保系统达到300 K 时的平衡状态（时间步长为 0.25 fs）。接近平衡状态后，模型沿传热方向分为50块板。如图 2a 所示，第 1 块板被指定为散热器，而第 26 块板（样品的中间板）是热源，热通量从热源（热区）转移到散热器（寒冷地区）。热通量传输方向定义为长度方向 (L)，而横向定义为宽度 (W) 方向。通过在热沉板中具有最高动能的最热原子和在热源中具有最低动能的最冷原子之间交换动能，在这两个板之间释放/注入热通量 J平板。根据下式计算散热器与热源板之间的动能交换量即可得到热通量J。

一 rNEMD 方法的示意图。热通量从热源（热区）传递到散热器（冷区）。热通量传输方向定义为长度方向 (L)，而横向定义为宽度 (W) 方向。 b 平均温度随板坯的分布

其中 t swap 是交换动能的总时间，N swap 表示交换原子对的数量，m 是原子的质量，v h 和 v c 分别代表交换原子的速度（散热片中动能最高的最热原子和热源片中动能最低的最冷原子）。当系统达到非平衡稳态时（1.5 ns 后），每个板的温度被收集并平均超过 3.0 ns，以获得温度分布。 TC (κ ) 然后使用傅立叶定律计算为

其中 A 是传热的横截面积 (A 乘以模型的宽度和厚度得到），和∂T /∂L 表示系统达到非平衡稳态后的温度梯度（见图 2b）。因子 2 表示热通量在远离热源的两个方向上传输的事实。假设模型厚度为石墨烯的层间平衡间距（0.34 nm）[8, 10, 34, 35]。

结果和讨论

我们首先检查系统尺寸对两种碳同素异形体的 TC 的影响。生成的模拟样本宽度相同，为 3 纳米，但长度不同，从 50 到 1000 纳米不等。需要注意的是，本文中提到的样本长度的所有值都是有效长度（L 效果 ) 的传热。即有效样本长度为样本长度的一半（L )，即 L 效果 =L /2，这归因于rNEMD方法中从样品中间（热源）到两端（散热片）的热通量。特别是，我们通过计算固定长度为 50 nm 但宽度分别为 3 nm、6 nm、9 nm 和 12 nm 的样品的热导率，确认了 TC 不依赖于样品宽度，如图所示图 3 OPG-L 沿锯齿形和扶手椅方向的 TC 命名为 κ OPG-LZ 和 κ OPG-LA，分别。类似地，κ OPG-ZZ 和 κ OPG-ZA 用于表示 OPG-Z 沿锯齿形和扶手椅方向的 TC。模拟结果表明，OPG-L 和 OPG-Z 在两个手性方向上的 TC 随着样品长度从 50 到 1000 nm 变化而单调增加。这归因于在长样本中，波长较长的声子参与传热[9, 36]。 50 nm 和 1000 nm 长的 OPG-L 和 OPG-Z 沿锯齿形方向的 TC 分别为 κ OPG-LZ50 =125 W/mK，κ OPG-LZ1000 =296 W/mK，κ OPG-ZZ50 =94 W/mK，并且 κ OPG-ZZ1000 =236 W/mK。沿扶手椅方向，OPG-L和OPG-Z的TC分别为κ OPG-LA50 =105 W/mK，κ OPG-LA1000 =316 W/mK，κ OPG-ZA50 =93 W/mK，和 κ OPG-ZA1000 =214 W/mK。

OPG-L 和 OPG-Z 的 TC 作为宽度的函数

为了提取无限长样本的 TC，采用逆拟合程序。逆长与逆TC的关系表示为[37,38,39]：

其中 κ ∞ 是无限样本的外推 TC，l 是声子平均自由程，L 效果 是传热的有效长度。等式(3)表明反长度和反TC之间的关系应该是线性的。如图 4 所示，观察到反长度和反 TC 之间的线性关系。通过外推到 L ⁻¹ =0，无限样本的TC，即κ OPG-LZ =310 W/mK，κ OPG-LA =332 W/mK，κ OPG-ZZ =247 W/mK，并且 κ 得到OPG-ZA =228 W/mK。

a 的逆 TC OPG-L 和 b OPG-Z 作为 300 K 时样品反长度的函数。空心蓝色菱形和红色点分别代表沿锯齿形和扶手椅方向的 TC

此外，我们还在平衡分子动力学 (EMD) 方法中通过建立相同长度和宽度为 20 nm 的样品来表达运行 TC（此模拟样品尺寸已被测试为足够大以消除有限尺寸效应） .根据 Fan 等人的工作。 [39, 40]，EMD方法中的TC计算基于Green-Kubo公式[41, 42]，其中沿x方向运行的TC可表示为：

其中 κ B 是玻尔兹曼常数，V 是系统的体积，T 是系统的绝对温度，<J x (0)J x (t ^' )〉是热通量自相关函数，t 是相关时间，J x 是 x 方向的热通量。符号〈〉代表 EMD 模拟中的时间平均值。最大相关时间为 2 ns，经测试足够大。如图 5 所示，OPG-L 和 OPG-Z 在 300 K 的两个手性方向上的运行 TC 是通过对不同初始速度的 100 次独立模拟的结果求平均值来表示的。我们可以通过对 1.0 到 2.0 ns 的相关时间中运行的 TC 进行平均来进一步获得无限样本的 TC。即 OPG-LZ、OPG-LA、OPG-ZZ 和 OPG-ZA 的收敛 TC 分别为 313 W/mK、344 W/mK、261 W/mK 和 233 W/mK，分别为与rNEMD方法外推的结果合理吻合。

a TC 的演变 OPG-LZ，b OPG-LA，c OPG-ZZ 和 d OPG-ZA 在 300 K 时作为相关时间的函数。细线代表 100 次独立模拟的结果，粗实线和虚线代表它们的平均值和误差范围。 κ ∞ 是无限样本的 TC，它是通过在相关时间从 1.0 到 2.0 ns 中对运行的 TC 进行平均而获得的

发现这两种碳同素异形体的 TC 远低于石墨烯 (3000–5000 W/mK) [7, 43]。 为了解释这种现象并探索物理洞察力，我们计算了三个重要参数，即 C v , v g , 和 l , 基于经典的晶格热传输方程：

其中Cv为热容，vg为有效声子群速度，l为声子平均自由程。

采用长度和宽度均为 20 nm 的样品来研究 300 K 时的热容量。 热容量是按照 McGaughey 和 Kaviany [44] 的方法计算的，该方法已用于平衡分子动力学的方法模拟 [45]。我们计算总能量 E 在 T 的温度下 =290 K、295 K、300 K、305 K、310 K 在规范系综中，结果在具有不同初始速度的 10 次独立模拟的 60 ps 上取平均值。如图6所示，能量-温度曲线线性拟合的斜率为热容。

a 的能量随温度的变化 OPG-L 和 b OPG-Z。能量-温度曲线的斜率表示热容量。每个原子对应的热容分别为 4.163 E-23 J/K 和 4.126 E-23 J/K

需要注意的是，我们这里计算的声子群速度是有效声子群速度v g 而不是平均声子群速度 v .如图 7 所示，可以通过比较 rNEMD 和 EMD 模拟的结果来获得有效声子群速度。即有效系统长度L 效果 可以在 EMD 方法中通过乘以相关时间的上限 t 来定义 在 Green-Kubo 公式中。 (4) 由有效声子群速度 v g , L 效果 ≈ v g t .正在运行的 TC κ (t ) 的 EMD 方法也可以看作是系统长度 κ 的函数 (L 效果 ）。与平均声子群速度相比，有效声子群速度是粗略估计，但它已广泛用于研究低维晶格模型中的热传输 [46]，也已用于石墨烯 [40] 和同素异形体司[39]。

a 的 TC OPG-LZ，b OPG-LA，c OPG-ZZ 和 (d ) OPG-ZA 作为来自 EMD 和 rNEMD 模拟的有效样本长度的函数。有效声子群速度V g 结合EMD和rNEMD模拟得到

基于方程。 (3)、声子平均自由程可以通过rNEMD方法外推得到。为了比较这两种碳同素异形体的 TC 与石墨烯的 TC，我们还介绍了石墨烯的这三个参数。石墨烯的热容通过上述方法计算，而有效声子群速度和声子平均自由程在其他工作中获得[7, 40]。可以发现，这两种碳同素异形体的热容接近于石墨烯；然而，有效声子群速度和声子平均自由程远低于石墨烯，导致两种材料的 TC 较低（见表 1）。

此外，我们探索了 TC 对温度的依赖性，如图 8 所示。200 K 到 300 K 的温度区域是我们关注的主要范围。生成的模拟样本宽度相同，为 3 纳米，但长度不同，分别为 50 纳米、75 纳米、100 纳米、150 纳米和 200 纳米。如图 8a、b 所示，我们给出了 OPG-LZ 和 OPG-LA 在不同温度下的逆 TC，作为逆样本长度的函数。与 300 K 时的尺寸效应外推类似，通过外推程序提取无限样品在不同温度下的热导率。如图 8c、d 所示，所有收敛的热导率均由 300 K (κ 0).

a 的逆 TC OPG-LZ，b OPG-LA 在不同温度下作为逆样本长度的函数，以及相对 TC (κ /κ 0) 的 c OPG-L 和 d OPG-Z 作为温度的函数。 κ 0 是 300 K 时的 TC，对于 κ，其为 310 W/mK、332 W/mK、247 W/mK 和 227 W/mK OPG-LZ，κ OPG-LA，κ OPG-ZZ 和 κ OPG-ZA，分别

图 8 表明，沿着之字形和扶手椅方向，OPG-L 和 OPG-Z 的 TC 都随着温度的升高而降低。 TC 随温度（从 200 到 500 K）变化的趋势与之前石墨烯 TC 研究的趋势非常一致 [8, 36, 47]。这种现象源于 Umklapp 散射过程的增强，该过程在热传输中起着关键作用 [8, 36, 47]。此外，当温度从 300 到 500 K 变化时，κ OPG-LZ，κ OPG-LA，κ OPG-ZZ 和 κ OPG-ZA 分别下降了 42%、40%、36% 和 37%。这两种碳同素异形体的TC对温度的依赖性表明，在实际应用中必须考虑温度的影响。

二维材料的热性能，例如石墨烯 [48, 49]、硅烯 [34, 50, 51] 和磷烯 [37]，对应变工程很敏感。据报道，小尺寸石墨烯的TC随着拉伸应变的增加而降低[48]，并且当样品大于500μm时，TC也可以通过增加应变来增强[49]。 TC 对样本大小和应变的异常依赖归因于边界散射和声子-声子散射之间的竞争。此外，由于面内模式中的声子软化和面外模式中的声子硬化之间的竞争，发现硅烯的 TC 在小拉伸应变下增加，但在大应变下降低 [34, 50, 51]。因此，研究OPG-L和OPG-Z结构的TC行为与拉伸应变之间的关系具有重要意义和必要性。

我们首先研究了这两种碳同素异形体的机械性能。样本大小约为 5 纳米长和 5 纳米宽。为了避免任何虚假的高键合力和非物理应变硬化 [52, 53]，截止距离固定在 (R =S =1.95 Å）。修改后的优化 Tersoff 电位中的截止距离也与之前用于模拟 C-C 键的 Tersoff 电位 (1.8–2.1 Å) [28, 53,54,55] 中的截止距离一致。所有模拟都是通过将结构的原子配置松弛到最小势能状态来启动的。施加单轴拉伸应变，应变率为 0.0002 ps ^-1 .需要注意的是，石墨烯的层间平衡间距（3.4 Å）用于表示两种结构的层间平衡距离。表 2 中列出了这两种碳同素异形体的机械性能，并比较了石墨炔和石墨烯 [56]。 z 的上标特征 和 a 分别代表锯齿形和扶手椅床单。

从表 2 中可以看出，沿锯齿形方向，OPG-L 和 OPG-Z 的杨氏模量分别为 538 GPa 和 492 GPa，沿扶手椅方向，杨氏模量分别为 648 GPa 和 550 GPa。这表明 OPG-L 和 OPG-Z 的杨氏模量接近于石墨炔 (503.1 ^z 和 525.0 ^a ) 但低于石墨烯 (856.4 ^z 和 964.0 ^a ）。沿锯齿形和扶手椅方向的两种碳同素异形体的应力应变关系如图 9 所示。根据这两种碳同素异形体的断裂行为，我们进一步得到了这两种碳同素异形体的极限应变（张力）。分别沿Z字形方向，OPG-L和OPG-Z的极限应变（拉力）分别为17.2%和10.9%，沿扶手椅方向，极限应变（拉力）分别为8.7%和7.9%。我们发现 OPG-L 的结构在锯齿形方向的拉伸应变下具有更高的强度。然而，与石墨炔和石墨烯相比，两种碳同素异形体的极限应变（张力）较低。

两种碳同素异形体沿锯齿形和扶手椅方向的应力-应变关系

然后，我们通过沿传热方向施加单轴拉伸应变来研究这两种碳同素异形体对 TC 的应变效应。模拟样品具有相同的 3 纳米宽度但不同的长度分别为 50 纳米、75 纳米、100 纳米、150 纳米和 200 纳米。通过外推程序提取无限样品在各种应变下的热导率（见图 10a、b）。如图 10c 和 d 所示，所有收敛的热导率都通过 300 K (κ 0)，我们进一步给出相对 TC (κ /κ 0) 作为各种单轴应变的函数的两个碳同素异形体。图 10 清楚地表明，OPG-L 和 OPG-Z 的 TC 都随着拉伸应变的增加而单调降低，这与之前对石墨烯的研究一致 [34, 48]，但与硅烯 [34, 50, 51] 和磷烯[37]。如图 10 所示，κ 的最大减少 OPG-LZ，κ OPG-LA 和 κ OPG-ZZ，κ OPG-ZA 分别为 49%、44%、37% 和 31%。特别是OPG-L沿Z字形方向的TC可以通过应变在大范围内进行调节。

a 的逆 TC OPG-LZ，b OPG-LA 在各种单轴应变下作为逆样本长度的函数，以及相对 TC (κ /κ 0) 的 c OPG-L 和 d OPG-Z 作为应变的函数。 κ 0 是 300 K 时无应力的 TC，对于 κ，分别为 310 W/mK、332 W/mK、247 W/mK 和 227 W/mK OPG-LZ，κ OPG-LA，κ OPG-ZZ 和 κ OPG-ZA，分别

为了进一步阐明应变对 OPG-L 和 OPG-Z 热传输特性的影响，我们计算了 OPG-LZ 在典型应变下的声子的 VDOS。 VDOS 是通过原子速度的自相关函数的傅立叶变换来计算的。函数定义如下：

如图 11 所示，观察到面内和面外方向的声子软化（红移）。这种现象与之前在拉伸应变下对石墨烯的研究非常吻合 [34, 48]。特别是，与面外方向的VDOS相比，面内方向的声子软化是明显的。这表明OPG-L和OPG-Z的TC下降主要是由于面内方向的应变诱导声子软化。

平面内a 和平面外b OPG-L 的 VDOS 与沿锯齿方向的单轴拉伸应变

结论

总之，已经进行了 EMD 和 rNEMD 模拟来研究由八边形和五边形组成的两种新型二维碳同素异形体的热性能。获得了尺寸、温度和应变对 TC 的影响。我们的结果表明，TC 随着尺寸的增加而单调增加。无限大的热导率是通过反长度和反 TC 的线性关系获得的。发现逆非平衡分子动力学方法通过外推得到的收敛TC与平衡分子动力学方法具有合理的一致性。与石墨烯相比，TC 低得多，这归因于较低的声子群速度和声子平均自由程。我们的研究结果为尺寸、温度和应变对 OPG-L 和 OPG-Z 热传输特性的影响提供了重要见解，并表明其在微/纳米电子领域的热管理器件中的潜在应用。

缩写

558：

五五八元环

OPG-L：

八边形和五边形石墨烯线

OPG-Z：

八边形和五边形石墨烯-锯齿形

rNEMD：

逆非平衡分子动力学

TC：

Thermal conductivity

VDOS:

Vibrational density of states