石墨烯-碳納米管復(fù)合結(jié)構(gòu)熱導(dǎo)率的非平衡分子動力學(xué)模擬

摘要：為制備高性能的熱界面材料，采用非平衡分子動力學(xué)方法計(jì)算了三維石墨烯–碳納米管復(fù)合結(jié)構(gòu)的熱學(xué)特性，研究了復(fù)合結(jié)構(gòu)的尺寸及系統(tǒng)溫度對于整體結(jié)構(gòu)熱導(dǎo)率的影響。研究表明，復(fù)合結(jié)構(gòu)z 方向長度的增加和溫度的升高都能使復(fù)合結(jié)構(gòu)的熱導(dǎo)率在一定程度上增加并逐漸達(dá)到一個(gè)飽和值。為進(jìn)一步探究整體及局部結(jié)構(gòu)中聲子能量傳遞的情況，計(jì)算了復(fù)合整體和連接處界面兩側(cè)的聲子振動態(tài)密度，結(jié)果表明，碳納米管與石墨烯界面處形變是阻礙復(fù)合結(jié)構(gòu)熱導(dǎo)率提高的主要因素。

關(guān)鍵詞：熱界面材料；非平衡分子動力學(xué)；熱導(dǎo)率；振動態(tài)密度

中圖分類號：TK 12 文獻(xiàn)標(biāo)志碼：A

隨著半導(dǎo)體技術(shù)的發(fā)展，芯片的特征尺寸不斷縮小，高性能微處理器的功率密度急劇上升[1]。由此產(chǎn)生的大量熱量需要良好的散熱措施，以保證電子器件的性能和使用壽命。幾乎所有的大功率電子器件都會配合銅或鋁制散熱器將熱量導(dǎo)出[2]，但由于絕對粗糙度的存在，使得熱源與散熱器接觸的界面存在空氣間隙。而空氣的導(dǎo)熱系數(shù)比金屬低約4 個(gè)數(shù)量級，導(dǎo)致熱量會在這些區(qū)域內(nèi)聚集，無法快速散出[3]。這時(shí)，在界面間的空隙中填充具有高導(dǎo)熱性能且形狀易加工的熱界面材料是提高熱源與散熱器之間傳熱效率的常用方法[4]，因此，尋找并制備高性能的熱界面材料也逐漸成為當(dāng)前的研究熱點(diǎn)之一。

相比于傳統(tǒng)熱界面材料，碳材料本身具有優(yōu)異的導(dǎo)熱性能，并在改進(jìn)制備高性能熱界面材料方面展示出良好的應(yīng)用前景[5]。常見的碳材料主要有石墨烯[6]、碳納米管[7]、氧化石墨烯等[8]，對于石墨烯和碳納米管這兩種材料，它們的導(dǎo)熱性都具有明顯的各向異性，在實(shí)際應(yīng)用中存在一定的限制。2008 年，Dimitrakakis 等[9] 首次引入了三維石墨烯–碳納米管復(fù)合結(jié)構(gòu)（GCHs），并開始嘗試在氫能源儲存領(lǐng)域應(yīng)用它。2022 年，Jyoti 等[10] 進(jìn)一步探究了GCHs 這一三維結(jié)構(gòu)在能量轉(zhuǎn)換和儲存方面的應(yīng)用潛力。Park 等[11] 對GCHs 的熱傳導(dǎo)特性進(jìn)行了研究，發(fā)現(xiàn)這種復(fù)雜結(jié)構(gòu)不僅大幅度提高了石墨烯和碳納米管在熱傳導(dǎo)方向的機(jī)械特性，同時(shí)也顯示出卓越的熱傳導(dǎo)能力。Loh 等[12]利用分子動力學(xué)的方法，對GCHs 中的機(jī)械應(yīng)變對聲子動力學(xué)界面產(chǎn)生的影響進(jìn)行了深入的研究。Zhang 等[13] 發(fā)現(xiàn)錐形界面的熱傳導(dǎo)率優(yōu)于垂直界面。近年來，GCHs 的熱傳導(dǎo)性質(zhì)引起了許多研究者的關(guān)注[14]，但對它在微觀水平的熱傳導(dǎo)機(jī)制的理解仍是有限的，其復(fù)合結(jié)構(gòu)中各個(gè)部分例如碳納米管的長度、石墨烯片層的面積以及溫度等對熱傳導(dǎo)特性的影響也還沒有得到明確，對各種相關(guān)理論的完善和深入探討還需要更多的研究。本文利用郎之萬的非平衡分子動力學(xué)方法模擬了GCHs 的導(dǎo)熱率，運(yùn)用聲子振動態(tài)密度等解析方法深入研究了GCHs 內(nèi)部的熱傳遞機(jī)理，并考察了結(jié)構(gòu)規(guī)模和系統(tǒng)溫度對總體熱傳遞的影響。該研究成果對于GCHs 在熱界面材料中的實(shí)際應(yīng)用有豐富的益處，并能為高效熱界面材料的創(chuàng)新提供理論支持。

1 模型與模擬方法

1.1 模型的構(gòu)建

首先，構(gòu)造出面積為3.01 nm×4.88 nm、層面間的間距為2.4 nm 的多層石墨烯片狀組織，接下來在每塊石墨烯片中挖出稍微大一些的孔，孔的尺寸略超過碳納米管的直徑。接著，將長度達(dá)到2.35 nm 的碳納米管移動填充到對應(yīng)的孔洞位置，重新連接構(gòu)成三維石墨烯–碳納米管單晶體結(jié)構(gòu)。最后，在z 軸擴(kuò)大晶胞以構(gòu)建三維GCHs 的多晶胞結(jié)構(gòu)，如圖1 所示。

1.2 模擬方法

分子動力學(xué)計(jì)算熱導(dǎo)率的方法主要有兩種：平衡態(tài)方法和非平衡態(tài)方法[15]。本研究選取非平衡態(tài)的分子動力學(xué)方法對GCHs 熱導(dǎo)率進(jìn)行模擬運(yùn)算。這種非平衡態(tài)的分子動力學(xué)方式中，熱能傳播通過對系統(tǒng)施加擾動來執(zhí)行，然后通過運(yùn)算獲取熱流和溫度梯度，熱流僅沿z 軸傳遞。根據(jù)傅里葉的熱傳導(dǎo)法則，分別引出模擬體的熱傳導(dǎo)率和界面熱阻：

k= -Jz/ （dT/dz）（1）

式中：k為熱導(dǎo)率；Jz為沿z軸的熱流密度，表示 單位時(shí)間、單位截面積上的熱量流通；dT/dz為z軸方向的溫度梯度；負(fù)號表示熱量向溫度較低處進(jìn) 行傳導(dǎo)。

R =AΔT/Jz（2）

式中：R為界面熱阻；ΔT為界面處的溫度階躍。

關(guān)于勢函數(shù)的決定，考慮本文的研究系統(tǒng)完全由碳元素構(gòu)成，選擇使用經(jīng)過優(yōu)化的Tersoff 勢函數(shù)[16] 來描繪碳–碳原子間的作用，Tersoff 勢函數(shù)的表述方式如下：

式中：rij表示原子i和j之間的距離；фR和фA分別 用于描述原子間的排斥力和吸引力，下標(biāo)R和 A分別對應(yīng)排斥項(xiàng)和吸引項(xiàng)；Bij為鍵角作用力的 度量參數(shù)，由于Bij并非常量，因此Tersoff勢函數(shù) 并非對勢。

在非平衡分子動力學(xué)方法中，產(chǎn)生溫差的方法包括恒熱流法、動量交換法以及基于溫控器的恒溫法。在恒熱流法和動量交換法中，溫差都不能精確控制。因此，本文使用恒溫法來產(chǎn)生溫差。同時(shí)，鑒于郎之萬（Langevin）溫控器的隨機(jī)性和局部性，相比于胡佛（Nosé-Hoover）（鏈）溫控器的全局控溫，其性能更優(yōu)越。所以采用郎之萬溫控器[17] 來產(chǎn)生溫差。在郎之萬方法中，恒溫區(qū)粒子的運(yùn)動方程為：

式中：ri表示原子i的位置；pi表示原子i的動 能；Fi表示原子i受到的外力；mi為原子i的質(zhì) 量；τ是一個(gè)時(shí)間參數(shù)，作為恒溫器和系統(tǒng)之間的 耦合；fi是一個(gè)隨機(jī)力，是用來恢復(fù)正則分布的波 動-耗散關(guān)系決定的。

1.3 模擬過程

在研究過程中，本文運(yùn)用開源軟件LAMMPS[18]執(zhí)行所有的模擬計(jì)算。為了預(yù)防系統(tǒng)過高的溫度變動，設(shè)置z 軸方向最外層的兩層原子保持不變，剩余的部分原子分別設(shè)為冷源和熱源，這樣便可使熱量從熱源沿z 軸方向輸送至冷源。在這種設(shè)定下，將z 軸方向作為熱傳遞路徑并設(shè)為固定邊界，而非熱傳遞路徑設(shè)為周期性的邊界。在模擬的整個(gè)過程中，先在NVT （體積固定、溫度穩(wěn)定）系綜環(huán)境下開始進(jìn)行系統(tǒng)的能量最小化，以此來達(dá)到優(yōu)化后的穩(wěn)定狀態(tài)結(jié)構(gòu)。接下來，在NVT （體積固定、能量穩(wěn)定）系綜環(huán)境下進(jìn)行模擬，設(shè)定的溫度為300 K，時(shí)間步長為0.5 fs，而模擬步數(shù)達(dá)到300 萬步，總共模擬了1.5 ns。系統(tǒng)達(dá)到穩(wěn)定狀態(tài)后，使用Langevin 熱控器來引入熱流，構(gòu)成溫度梯度，并同時(shí)采用NVE 系統(tǒng)和Langevin 溫控器相結(jié)合的控制方式來調(diào)整系統(tǒng)的壓力和溫度。為了保障計(jì)算結(jié)果的精確度，對每一模擬流程作了3 次重復(fù)，并取其平均值作為最終結(jié)果。

2 結(jié)果與分析

2.1 熱導(dǎo)率的計(jì)算和驗(yàn)證

本文把模擬系統(tǒng)分成3 個(gè)部分：冷源、熱源和熱傳輸區(qū)。把冷源和熱源的溫度設(shè)定在270 K和320 K。當(dāng)系統(tǒng)平衡，得到GCHs 在平衡狀態(tài)下的溫度分布。將模擬的結(jié)果和文獻(xiàn)[11] 中的數(shù)據(jù)進(jìn)行了比較，在圖2 中能看出兩者高度一致。通過3 次重復(fù)模擬計(jì)算得到GCHs 的熱導(dǎo)率分別為1.45 ，1.46 ，1.50 W/（m·K），取其平均值為1.47 W/（m·K），與文獻(xiàn)中熱導(dǎo)率的值（1.42±0.05） W/（m·K） 相符合。

為了對比碳納米管對復(fù)合結(jié)構(gòu)導(dǎo)熱性能的作用，本文也計(jì)算得到了多層石墨烯結(jié)構(gòu)的熱導(dǎo)率，結(jié)果為0.289 W/（m·K）?？梢园l(fā)現(xiàn)三維GCHs的熱導(dǎo)率與之相比提高了一個(gè)量級，說明在石墨烯片層之間添加碳納米管可以大幅度提高復(fù)合結(jié)構(gòu)的法向?qū)崧剩?這是因?yàn)樘技{米管可以在GCHs 法向方向上構(gòu)建更高效的熱傳遞通道。

2.2 振動態(tài)密度分析

對于非金屬晶體的納米材料，其熱導(dǎo)率與晶格的振動有關(guān)，將晶格振動的能量量子化，其最小單位通過引入一種準(zhǔn)粒子——聲子來定義，通俗地講，聲子就是原子（或分子）晶格振動的集體激發(fā)。因此， 通過計(jì)算聲子的振動態(tài)密度[19]（vibrational density of states，VDOS），可以獲取晶格單位頻率間隔內(nèi)的振動狀態(tài)和能量在晶體結(jié)構(gòu)中的分布?；谶@些信息，可以剖析哪些原子或結(jié)構(gòu)在整體結(jié)構(gòu)中限制了能量的傳輸。聲子振動態(tài)密度是通過對原子速度的自相關(guān)函數(shù)（velocityautocorrelation function，VACF）進(jìn)行傅里葉變換而得來的，其計(jì)算公式如下：

本文借由計(jì)算出的GCHs 和多層石墨烯的VDOS 來研究熱導(dǎo)率與原子及其構(gòu)造的聯(lián)系。首先，計(jì)算出GCHs 的標(biāo)準(zhǔn)化速度自相關(guān)函數(shù)，如圖3 所示，原子的速度改變是無定向的，其變動特性不會在時(shí)間序列中復(fù)現(xiàn)。然后通過對速度自相關(guān)函數(shù)進(jìn)行傅里葉變換得到GCHs 的聲子振動態(tài)密度， 如圖4 所示。圖4 顯示的結(jié)果可以看出，GCHs 的VDOS 在低頻區(qū)相比于多層石墨烯的下降，但在中高頻區(qū)，GCHs 展現(xiàn)出多個(gè)高峰。由此可見，在GCHs 中，激發(fā)并參與振動的聲子數(shù)量更多，從而顯著提升了其熱傳輸效率。同時(shí)，對于GCHs，VDOS 的峰值主要集中于中低頻率段，且峰值從高頻率段向中低頻率段移動，而通過后續(xù)計(jì)算可知，中低頻率段的聲子對整體的熱導(dǎo)率貢獻(xiàn)更大。

為深化對GCHs 熱傳導(dǎo)特性的認(rèn)識，本文詳細(xì)研究了聲子振動頻率如何影響熱導(dǎo)率。首先，對圖4 中數(shù)據(jù)進(jìn)行了以20 THz 頻率間隔的反傅里葉轉(zhuǎn)化，然后計(jì)算了特定頻率范圍內(nèi)聲子對熱導(dǎo)率的影響，如圖5 所示。結(jié)果顯示，盡管在所有頻率段，聲子都由于不同程度的激活參加了熱能傳遞，但主要是0～20 THz，20～40 THz 的低頻聲子起了主要作用，對熱導(dǎo)率貢獻(xiàn)最大。這主要是因?yàn)?，在聲子從碳管傳遞到石墨烯再到碳管的過程中，振動的不匹配會影響聲子的能量傳輸，導(dǎo)致能量重新分配。這個(gè)時(shí)候，界面處會產(chǎn)生較強(qiáng)的聲子散射，而低頻聲子的平均自由程比高頻聲子大，對熱導(dǎo)率的貢獻(xiàn)更大，散射越強(qiáng)，熱導(dǎo)率越小。

觀察圖2 中的溫度分布圖還可以發(fā)現(xiàn)，GCHs的溫度躍升主要呈現(xiàn)在碳納米管和石墨烯結(jié)合的部分，這表明在GCHs 的構(gòu)造中，石墨烯片層與碳納米管接合區(qū)的形變是垂直熱傳遞的最大阻礙。為了進(jìn)一步分析連接界面處聲子能量傳遞的情況，計(jì)算連接界面處兩側(cè)原子的聲子振動態(tài)密度分布如圖6 所示。觀察二者的振動態(tài)密度分布圖可以發(fā)現(xiàn)，二者的振動態(tài)密度峰值分布的頻率基本不一致，這說明在此連接界面處兩側(cè)原子的聲子振動模式耦合效果不理想，從而導(dǎo)致兩側(cè)聲子振動的不匹配，對復(fù)合結(jié)構(gòu)在法向方向的熱量傳遞帶來極大的限制。

2.3 不同碳納米管長度對體系熱導(dǎo)率的影響

為了探究不同的碳納米管長度對體系熱導(dǎo)率的影響，設(shè)置三維GCHs 復(fù)合結(jié)構(gòu)在z 方向的尺寸分別為5，10，15，20 nm，計(jì)算其熱導(dǎo)率，結(jié)果如圖7 所示。從圖中可以看到GCHs 體系的熱導(dǎo)率隨著碳納米管長度的增加而增加，這是因?yàn)閱螌邮┰诜ㄏ蚍较蛏系膶?dǎo)熱能力極差，對于整個(gè)復(fù)合結(jié)構(gòu)的體系來說，法向方向上的導(dǎo)熱能力主要由碳納米管的軸向?qū)崮芰Q定，碳納米管越長導(dǎo)熱越強(qiáng)。從圖中還可以觀察到，當(dāng)GCHs 體系在z 方向尺寸較小時(shí)，熱導(dǎo)率隨碳納米管長度的增加大幅度增加，而后當(dāng)GCHs 體系在z 方向的尺寸逐漸增大后，體系的熱導(dǎo)率值隨碳納米管長度的增加逐漸放緩。這主要受兩方面因素的影響：一方面，隨著碳納米管長度的增加，接近其平均自由程，其內(nèi)部的導(dǎo)熱方式逐漸由彈道輸運(yùn)轉(zhuǎn)變?yōu)閺椀括C擴(kuò)散輸運(yùn)；而另一方面，隨著石墨烯片層之間的距離越來越遠(yuǎn)，相互之間的范德華力也隨之變小，其輔助導(dǎo)熱的能力也隨之下降。在這兩方面的共同作用下，雖然碳納米管的長度增加，復(fù)合結(jié)構(gòu)的熱導(dǎo)率仍是增長的，但增長速度是逐步變緩的。

2.4 系統(tǒng)溫度對體系熱導(dǎo)率的影響

為了探究GCHs 作為熱界面材料在不同溫度下的工作情況， 分別設(shè)置了GCHs 在平衡溫度為100～700 K 情況下的導(dǎo)熱情況，計(jì)算得到的熱導(dǎo)率如圖8 所示。從圖中可以發(fā)現(xiàn)：隨著體系溫度的上升，復(fù)合結(jié)構(gòu)的熱導(dǎo)率開始是隨之近似線性增加的，這是因?yàn)榉肿拥臒徇\(yùn)動隨體系溫度的上升而增大，提升了熱傳遞的效率；而當(dāng)體系溫度達(dá)到500 K 后，由于溫度較高時(shí)，原子振動會加劇，從而導(dǎo)致不同的碳納米管之間會發(fā)生相互作用，進(jìn)而影響復(fù)合結(jié)構(gòu)法向方向熱量的輸運(yùn)過程，熱導(dǎo)率因此基本不再增加。

3 結(jié) 論

本文采用郎之萬控溫策略的非平衡分子動力學(xué)進(jìn)行了三維石墨烯–碳納米管復(fù)合結(jié)構(gòu)熱學(xué)特征的模擬，并考察了碳納米管的長度及體系溫度對復(fù)合結(jié)構(gòu)熱導(dǎo)率的影響，從而得出以下結(jié)論：

a．由于碳納米管在法向方向上構(gòu)建了更高效的熱傳遞通道，相對于具有0.29 W/（m·K） 熱導(dǎo)率的多層石墨烯結(jié)構(gòu)，GCHs 的熱導(dǎo)率得到了顯著增加，提升到了1.47 W/（m·K）。

b．碳納米管長度越長，GCHs 復(fù)合結(jié)構(gòu)的熱導(dǎo)率越大。并且當(dāng)復(fù)合結(jié)構(gòu)在z 方向上的尺寸較小時(shí)，熱導(dǎo)率增加的幅度較大，而后隨著復(fù)合結(jié)構(gòu)尺寸的增大，熱導(dǎo)率增加的趨勢逐漸變緩慢。

c．隨著體系溫度的上升，復(fù)合結(jié)構(gòu)的熱導(dǎo)率開始是隨之近似線性增加，而后當(dāng)體系達(dá)到一定溫度時(shí)，熱導(dǎo)率逐漸趨于穩(wěn)定。

參考文獻(xiàn)：

[1] GARIMELLA S V， PERSOONS T， WEIBEL J A， et al.Electronics thermal management in information and communications technologies： challenges and future"directions[J]. IEEE Transactions on Components， Packaging and Manufacturing Technology， 2017， 7（8）： 1191-1205.

[2] ZHUANG D W， YANG Y F， DING G L， et al.Optimization of microchannel heat sink with rhombus fractal-like units for electronic chip cooling [J].International Journal of Refrigeration， 2020， 116： 108-118.

[3] YU W H， FRANCE D M， ROUTBORT J L， et al. Review and comparison of nanofluid thermal conductivity and heat transfer enhancements[J]. Heat Transfer Engineering， 2008， 29（5）： 432-460.

[4] XING W K， XU Y， SONG C Y， et al. Recent advances in thermal interface materials for thermal management of high-power electronics[J]. Nanomaterials， 2022， 12（19）： 3365.

[5] GUO X X， CHENG S J， CAI W W， et al. A review of carbon-based thermal interface materials： mechanism， thermal measurements and thermal properties[J]. Materials amp; Design， 2021，209： 109936.

[6] 茹越，戚桂村，王湘，等.聚合物基石墨烯復(fù)合材料導(dǎo)熱性能的研究進(jìn)展[J].石油化工，2022， 51（7）： 844-849.

[7] FABRIS D， ROSSHIRT M， CARDENAS C， et al. Application of carbon nanotubes to thermal interface materials[J]. Journal of Electronic Packaging， 2011， 133（2）： 020902.

[8] 陳利濤，劉妮，劉庭崧，等.碳納米粒子促進(jìn)水合物生成 的研究進(jìn)展[J].能源研究與信息，2022， 38（2）： 109-116.

[9] DIMITRAKAKIS G K， TYLIANAKIS E， FROUDAKIS G E. Pillared graphene： a new 3-D network nanostructure for enhanced hydrogen storage[J]. Nano Letters， 2008， 8（10）： 3166-3170.

[10] JYOTI J， GUPTA T K， SINGH B P， et al. Recent advancement in three dimensional graphene-carbon nanotubes hybrid materials for energy storage and conversion applications[J]. Journal of Energy Storage， 2022， 50： 104235.

[11] PARK J， PRAKASH V. Thermal transport in 3D pillared SWCNT-graphene nanostructures[J]. Journal of Materials Research， 2013， 28（7）： 940-951.

[12] LOH G C， TEO E H T， TAY B K. Tuning the Kapitza resistance in pillared-graphene nanostructures[J]. Journal of Applied Physics， 2012， 111（1）： 013515.

[13] ZHANG Z A， KUTANA A， ROY A， et al. Nanochimneys： topology and thermal conductance of 3D nanotube-graphene cone junctions[J]. The Journal of Physical Chemistry C， 2017， 121（2）： 1257-1262.

[14] KHANNA V， KUMAR V， BANSAL S A. Mechanical"properties of aluminium-graphene/carbon nanotubes （CNTs） metal matrix composites： advancement，opportunities and perspective[J]. Materials Research Bulletin， 2021， 138： 111224.

[15] LI Z， XIONG S Y， SIEVERS C， et al. Influence of thermostatting on nonequilibrium molecular dynamics simulations of heat conduction in solids[J]. The Journal of Chemical Physics， 2019， 151（23）： 234105.

[16] TERSOFF J. Modeling solid-state chemistry： interatomic potentials for multicomponent systems[J]. Physical Review B， 1989，39（8）： 5566-5568.

[17] HU Y， FENG T L， GU X K， et al. Unification of nonequilibrium molecular dynamics and the mode-resolved phonon Boltzmann equation for thermal transport simulations[J]. Physical Review B， 2020， 101（15）： 155308.

[18] PLIMPTON S. Fast parallel algorithms for short-range molecular dynamics[J]. Journal of Computational Physics， 1995， 117（1）： 1-19.

[19] CUI L， FENG Y H， ZHANG X X. Enhancement of heat conduction in carbon nanotubes filled with fullerene molecules[J]. Physical Chemistry Chemical Physics， 2015， 17（41）： 27520-27526.

（編輯：丁紅藝）

基金項(xiàng)目：國家自然科學(xué)基金資助項(xiàng)目（51476102）