石墨烯-碳納米管復合結構熱導率的非平衡分子動力學模擬

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

關鍵詞：熱界面材料；非平衡分子動力學；熱導率；振動態密度

中圖分類號：TK 12 文獻標志碼：A

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

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

1 模型與模擬方法

1.1 模型的構建

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

1.2 模擬方法

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

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

式中：k為熱導率；Jz為沿z軸的熱流密度，表示 單位時間、單位截面積上的熱量流通；dT/dz為z軸方向的溫度梯度；負號表示熱量向溫度較低處進 行傳導。

R =AΔT/Jz（2）

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

關于勢函數的決定，考慮本文的研究系統完全由碳元素構成，選擇使用經過優化的Tersoff 勢函數[16] 來描繪碳–碳原子間的作用，Tersoff 勢函數的表述方式如下：

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

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

式中：ri表示原子i的位置；pi表示原子i的動 能；Fi表示原子i受到的外力；mi為原子i的質 量；τ是一個時間參數，作為恒溫器和系統之間的 耦合；fi是一個隨機力，是用來恢復正則分布的波 動-耗散關系決定的。

1.3 模擬過程

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

2 結果與分析

2.1 熱導率的計算和驗證

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

為了對比碳納米管對復合結構導熱性能的作用，本文也計算得到了多層石墨烯結構的熱導率，結果為0.289 W/（m·K）。可以發現三維GCHs的熱導率與之相比提高了一個量級，說明在石墨烯片層之間添加碳納米管可以大幅度提高復合結構的法向導熱率， 這是因為碳納米管可以在GCHs 法向方向上構建更高效的熱傳遞通道。

2.2 振動態密度分析

對于非金屬晶體的納米材料，其熱導率與晶格的振動有關，將晶格振動的能量量子化，其最小單位通過引入一種準粒子——聲子來定義，通俗地講，聲子就是原子（或分子）晶格振動的集體激發。因此， 通過計算聲子的振動態密度[19]（vibrational density of states，VDOS），可以獲取晶格單位頻率間隔內的振動狀態和能量在晶體結構中的分布。基于這些信息，可以剖析哪些原子或結構在整體結構中限制了能量的傳輸。聲子振動態密度是通過對原子速度的自相關函數（velocityautocorrelation function，VACF）進行傅里葉變換而得來的，其計算公式如下：

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

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

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

2.3 不同碳納米管長度對體系熱導率的影響

為了探究不同的碳納米管長度對體系熱導率的影響，設置三維GCHs 復合結構在z 方向的尺寸分別為5，10，15，20 nm，計算其熱導率，結果如圖7 所示。從圖中可以看到GCHs 體系的熱導率隨著碳納米管長度的增加而增加，這是因為單層石墨烯在法向方向上的導熱能力極差，對于整個復合結構的體系來說，法向方向上的導熱能力主要由碳納米管的軸向導熱能力決定，碳納米管越長導熱越強。從圖中還可以觀察到，當GCHs 體系在z 方向尺寸較小時，熱導率隨碳納米管長度的增加大幅度增加，而后當GCHs 體系在z 方向的尺寸逐漸增大后，體系的熱導率值隨碳納米管長度的增加逐漸放緩。這主要受兩方面因素的影響：一方面，隨著碳納米管長度的增加，接近其平均自由程，其內部的導熱方式逐漸由彈道輸運轉變為彈道–擴散輸運；而另一方面，隨著石墨烯片層之間的距離越來越遠，相互之間的范德華力也隨之變小，其輔助導熱的能力也隨之下降。在這兩方面的共同作用下，雖然碳納米管的長度增加，復合結構的熱導率仍是增長的，但增長速度是逐步變緩的。

2.4 系統溫度對體系熱導率的影響

為了探究GCHs 作為熱界面材料在不同溫度下的工作情況， 分別設置了GCHs 在平衡溫度為100～700 K 情況下的導熱情況，計算得到的熱導率如圖8 所示。從圖中可以發現：隨著體系溫度的上升，復合結構的熱導率開始是隨之近似線性增加的，這是因為分子的熱運動隨體系溫度的上升而增大，提升了熱傳遞的效率；而當體系溫度達到500 K 后，由于溫度較高時，原子振動會加劇，從而導致不同的碳納米管之間會發生相互作用，進而影響復合結構法向方向熱量的輸運過程，熱導率因此基本不再增加。

3 結 論

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

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

b．碳納米管長度越長，GCHs 復合結構的熱導率越大。并且當復合結構在z 方向上的尺寸較小時，熱導率增加的幅度較大，而后隨著復合結構尺寸的增大，熱導率增加的趨勢逐漸變緩慢。

c．隨著體系溫度的上升，復合結構的熱導率開始是隨之近似線性增加，而后當體系達到一定溫度時，熱導率逐漸趨于穩定。

參考文獻：

[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] 茹越，戚桂村，王湘，等.聚合物基石墨烯復合材料導熱性能的研究進展[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].能源研究與信息，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.

（編輯：丁紅藝）

基金項目：國家自然科學基金資助項目（51476102）