納米結構表面液體沸騰傳熱的分子動力學模擬

張石重，陳占秀，劉峰瑞，龐潤宇，王清

（河北工業大學能源與環境工程學院，天津 300401）

隨著電子器件的尺寸不斷微型化，功能高度集成及能耗不斷增加，對其散熱性能要求日趨苛刻，高效散熱方式成為制約電子設備性能和可靠性的關鍵問題。池沸騰強化傳熱沒有運動部件，裝置結構簡單，池內所需液體量小，適合微型器件狹窄空間內的高熱流密度散熱，其相變潛熱比單相冷卻方案高1～3 個數量級，被視作一種重要的有效冷卻技術。

微型器件狹窄空間內發生的池沸騰包含復雜的傳熱與傳質過程，加之相變過程的不穩定性，增加了對其實驗測試和數值模擬的困難。付婷采用激光加工和燒結技術在金屬表面制作微槽和納米線等不同形狀的納米結構，分析了這些納米結構對沸騰傳熱的影響。Ahn 等實驗研究了液體飽和狀態下光滑及微結構表面對池沸騰的影響。由于目前加工及測試儀器在時間和空間尺度上的限制，采用實驗測試方法研究微納尺度下相變傳熱依然是具有挑戰的課題之一。

采用數值模擬成為微納尺度行之有效的方法之一。其中分子動力學模擬方法是對分子無原則運動不進行任何簡化，是研究原子和分子水平復雜系統的高級確定性方法之一。其優勢在于能夠以原子級精度描述流體的運動特性，確定每一個原子軌跡，有效得出微納尺寸下的傳熱特性，因而成為研究納米尺度系統沸騰現象的重要工具。

諸多學者利用分子動力學模擬方法對納米尺度下的沸騰傳熱過程進行研究，Hens等采用分子動力學模擬研究光滑基底上納米級厚度液態氬加熱沸騰過程，指出對于處在同一溫度下的不同浸潤性表面，沸騰僅發生在潤濕性強的表面。類似地，Mao等利用模擬技術研究了液態水的加熱過程，觀察到水分子也發生了與Hens 等的工作中類似的沸騰現象。王艷紅等從能量角度對沸騰現象進行分析，發現液膜厚度和潤濕性共同影響著相變形式：當薄膜厚度較大時，弱的固液間作用力會促進沸騰現象的發生，當厚度較小時，為達到沸騰現象，對壁面提出了更高的溫度和潤濕性要求；由于較低的Kapitza 熱阻，在同一壁面溫度下，潤濕性強的系統溫度梯度更大，在更短時間內有氣泡產生。

學者們也對固體壁面上增加微結構強化沸騰傳熱進行了研究。Chen 等采用分子動力學（molecular dynamics，MD）模擬了凹槽為“V 形”和“柱形”及潤濕性對池沸騰液膜上氣泡成核行為的影響，發現凹槽結構加速了初始氣泡成核時間。Diaz等研究了納米柱在水平銅基底上氬膜的蒸發及池沸騰傳熱過程的影響，發現隨著納米柱間距離的增加，熱流密度隨之增加，最小間距（2.17nm）的熱流密度是最大間距（10.66nm）的63%左右。Liu 等利用Weierstrass-Mandelbrot 函數生成了隨機粗糙度表面，發現表面粗糙結構提高了換熱效率并促進了成核。借助于表面納米結構，固體和流體原子之間的熱傳遞變得更加有效。Liu 等通過構造凹形半球納米結構表面，并進行液態氬薄膜在納米結構表面上的沸騰傳熱的MD模擬，發現納米結構可以有效地增強沸騰傳熱，可以通過增加凹形納米結構的親水性來改善傳熱性能。

本文總結了實驗及分子動力學模擬結果，將池沸騰表面微結構概況表述為凹、凸微結構表面，將液膜放置在設有凹、凸半球狀納米結構的粗糙表面進行加熱以引發沸騰過程，并與光滑表面進行比較分析；通過在親水壁面至超疏水壁面7種不同潤濕性的表面構建納米結構，揭示納米級凹、凸微結構表面和潤濕性對液膜池沸騰的綜合影響過程。

1 物理模型介紹及模擬細節

本文采用液氬作為傳熱工質，主要原因是氬以單原子形式存在，氬與氬僅有范德華力作用，沒有分子間不對稱相互作用、氫鍵作用及分子間纏繞等，分子振動及碰撞運動單一，反映分子間流動與傳熱過程簡單，很多學者也采用液氬作為流動與傳熱工質。池沸騰模型中固體壁面以鉑（Pt）原子構建，構建的池沸騰結構尺寸參考相關文獻[12-19]，計算的體系中凹、凸結構更小，適用于更小的納米結構對傳熱的影響，同時將固-液間的作用力范圍拓展為超親水到超疏水狀態。

1.1 物理模型

建立如圖1 所示的模擬區域，該模擬體系為L×L×L=7.056nm×7.056nm×50.176nm 的 長 方 體 盒子，由固、液和蒸氣區域組成。固壁采用Pt原子，依照面心立方結構（FCC）的排列方式構建池沸騰底部金屬基板，晶格常數為0.392nm。按照從下至上的順序，固體基板分別被設置為固定用外部原子、熱源和接觸墻3個部分，分別起著預防基底形變、提供熱源和進行能量交換的作用。固體壁面分為3種結構表面，分別為平面、凸半球納米結構表面及凹半球納米結構表面。對于光滑表面，其固體厚度為=3.136nm，液膜厚度為=5.096nm；對于含凸起半球納米結構的粗糙表面，固體平面的厚度設置為=2.744nm，液膜厚度=5.488nm，納米結構半徑=1.96nm；對于含凹陷半球納米結構的粗糙表面，固體平面厚度=3.528nm，液膜厚度=4.9nm，納米結構半徑=1.96nm；以上液膜厚度均指液體頂層到平面部分的距離。這些設置保證了氬原子數量在7609 上下，其差距不超過7.5‰。改變壁面厚度及液膜厚度使得氬原子數保持基本不變，由于固體壁面導熱層設置較厚且液膜改變較小，這些改動不會對模擬結果造成影響。

圖1 物理模型結構示意圖

1.2 力場選擇及參數設置

將彈簧力施加于固體鉑原子，使其被約束在初始位置。工質選用氬原子，其性質穩定，勢函數簡單，在沿模擬框高度方向的最高處設置反射墻以確保粒子不會丟失。將Lennard-Jones 勢能（12-6 勢能）模型應用于氬原子、鉑原子之間的相互作用力的計算，其中截斷半徑設置為1.19nm。在本文所有模擬中，基底厚度都超過了其截斷半徑距離，這保證了壁面厚度不會影響模擬結果。L-J（12-6）勢能函數如式(1)所示。

式中，r為計算粒子對之間的間隔距離。在計算氬-氬原子對之間的相互作用力時，兩個參數選擇為=0.0104eV，=0.3405nm；當計算固體-固體原子對之間的相互作用力時，選取=0.52eV，=0.2475nm。

固體和液體、氣體之間的L-J 參數由Lorentz-Berthelot 混合法則進行求解，并使用調節系數對其進行調節，從而獲得不同潤濕性壁面。其數學形式如式(2)、式(3)所示。

1.3 模擬過程

圖2 密度法計算接觸角示意圖

表1 不同潤濕性表面參數

整個模擬過程分為3 個階段進行。首先利用Nose-Hoover 熱浴法進行85K 的熱浴，這個階段選取系綜為正則系綜（NVT）并使用共軛梯度法對系統能量進行最小化處理，能量公差和力公差均設置為10，保證系統達到穩定狀態且粒子速度遵循麥克斯韋速度分布。積分時間步長選擇為5fs，進行了3ns 的計算使系統穩定；隨后撤去整體溫度控制，利用Langevin 控溫將基板溫度繼續控制在85K，同時應用微正則系綜（NVE）對全部氬原子進行更新，此過程將會持續1ns，在模擬過程中監測整個系統的溫度、壓力和能量變化，確保系統在該過程結束之前達到平衡狀態，此過程為弛豫過程，保證加熱開始前的初始狀態；最后一階段利用Langevin恒溫器將熱源原子在15ns時間內以43K/ns的升溫速率進行升溫，其余原子在NVE 環境下進行更新。

2 結果與討論

2.1 加熱相變過程

弛豫后壁面上液態氬層的分布狀態如圖3 所示。對于光滑平面和凸半球納米結構表面，所用工況下的形態近似，但是隨著接觸角增大，固壁與液氬的距離增大。對于凹半球納米結構表面，潤濕性較強的E-case1 至E-case4 工況下，其凹槽內充滿工質氬；但是隨著接觸角增大，即在潤濕性較弱的E-case5 至E-case7 工況下，其凹槽內已不再含有工質氬，液氬與固壁間形成部分空腔，如圖3(c)所示，說明不同潤濕性下固壁與液體層的初始狀態不同。

圖3 不同壁面工況下的初始狀態

對上述初始狀態的液氬層進行加熱，加熱過程液氬沸騰快照如圖4所示，對此池沸騰過程中起始沸騰時間進行統計，通過質心位置確定沸騰起始時間，即每隔固定時間取一次液體質心位置，取=+ ?，此過程上升距離為?=-，當時間周期為t～t時，h- h< 0.2nm；同時在時間周期為t～t時，h- h> 0.2nm，認為t時刻為沸騰現象發生的沸騰起始時間。統計起始沸騰時間如表2所示。

表2 平面、凹、凸結構表面各工況起始沸騰時間

對于接觸角小即潤濕性較強的E-case1、Ecase2、E-case3工況，液氬受熱后發生爆炸沸騰呈現大致相同的過程，如圖4(a)、(b)、(c)所示，在Ecase1的強親水潤濕性下，光滑平面起始沸騰發生在2.7ns左右，凸半球納米結構表面和凹半球納米結構表面起始沸騰時間分別是2.25ns和2.175ns；對比親水條件下逐漸減弱的E-case1、E-case2、E-case3工況，平面上液氬沸騰起始時間由2.7ns 左右逐漸延長為5.4ns 左右，凸半球納米結構表面起始沸騰時間由2.25ns左右逐漸延長為4.6ns左右，凹半球納米結構表面起始沸騰時間由2.175ns 左右逐漸延長為4.35ns左右，說明在親水條件下，相同結構的表面，潤濕性減弱延后了氣泡出現的起始時間；對于同一潤濕性下平面和含凸、凹半球納米結構的三種親水性表面，在E-case1親水工況下，凹半球納米結構表面發生沸騰所需時間最快，凸半球納米結構表面次之，平面所需時間最長，說明添加納米結構會促進氣泡提前發生，縮短起始沸騰時間。

疏水工況下E-case4、E-case5 如圖4(d)和(e)所示，含凸半球納米結構表面出現氣泡早于含凹半球納米結構表面，平面還是最晚，與親水E-case1、E-case2、E-case3 工況不同；E-case4 工況下三種表面均早于E-case5 工況下三種表面；超疏水Ecase6及E-case7，平面和含凹半球納米結構表面都在10.75ns 有氣泡產生，可以計算到起始沸騰時間點，而含凸半球結構表面氣泡產生時間點不明顯；對于更為疏水工況E-case7，加熱至計算完成，只有平面上有氣泡產生，約在12.8ns，含凹、凸半球納米結構表面均沒有明顯的氣泡產生，液膜離開固壁的距離也很小，主要表現為液膜上部氬原子蒸發過程。

圖4 不同工況壁面快照信息

2.2 相變過程中氣相分子數的變化

本文采用Maroo 等判定氣態原子的方法判斷加熱過程中氣相分子數量，以某原子為中心確定其周圍5.3nm半徑球體內氬原子數目，當該原子周圍有不超過7 個氬原子時，認為該原子處于氣體狀態。通過這種方法計算了上述工況下的氣態原子占比，如圖5所示，其數值上等于氣態原子數目與總的氬原子數目之比。

圖5 不同表面加熱過程氣態原子所占比例

圖5所示為加熱過程中氣態原子生成分數，曲線斜率表示氣態原子生成速率。由圖5發現氣態原子變化可以分為兩個階段，曲線拐點位置在起始沸騰時間點附近。第一階段是開始加熱過程，氣態氬原子上升速度很快，氣態氬逐漸擺脫液氬層進入氣體空間，說明在此階段液氬層原子在固壁附近吸熱，固壁的熱量通過液氬層傳遞到其上層氬原子，上層氬原子吸熱后以蒸發方式進入氣體空間；第二階段是沸騰加熱過程，氣體分子數也在增加，但是增速減緩，直至全部液體原子加熱后轉變為氣態原子，之后氣態原子不再增加，曲線不再上升，Ecase1至E-case4工況下在計算時間內液氬全部轉變為氣態氬；曲線隨著潤濕性減弱，拐點位置逐漸下移，說明潤濕性減弱，液態原子轉換為氣態原子的速度也隨之降低，E-case5 至E-case7 工況下在計算時間內液氬未全部轉變為氣態氬，且對于超疏水E-case6 和E-case7，液氬轉變為氣態氬的數量未超過30%。對比三種表面上氣態氬生成速度，Ecase1至E-case4工況下，含凹半球微結構表面略大于凸半球微結構表面，兩者都大于平面；而對于E-case5 工況，三種表面上氣態生成速度相差不大，對于E-case6 和E-case7，平面速度最快，含凹半球微結構表面次之，凸半球微結構表面最慢。

2.3 不同工況下氬膜層溫度變化

溫度變化曲線如圖6～圖8所示，溫度曲線的斜率反映了液氬層升溫速率。圖6 是親水工況Ecase1 至E-case3 下液氬層溫度隨時間的變化趨勢，溫度曲線變化過程中也有拐點，此拐點位置也在起始沸騰時間點附近。相同親水工況下均顯示出含凹半球納米結構表面上液氬溫度最高，升溫速率最快，含凸半球納米結構表面次之，平面工況中液氬溫度最小和升溫速率最慢的順序，說明相同親水情況下納米結構的存在能夠強化換熱，但潤濕性占據主要影響因素，這與前人的研究結果一致。

圖6 親水壁面溫度隨時間變化情況

溫度曲線也可分為兩段，曲線轉折發生也在起始沸騰時間點附近（約2.7ns），強親水工況Ecase1 與E-case2 和E-case3 比較，溫度升高轉折變化比較明顯，反映了起始沸騰前液氬在固壁附近的快速升溫過程；發生起始沸騰后，液氬層離開壁面，是升溫速率減慢的過程。E-case1 至E-case3工況，隨著潤濕性減弱，升溫速度減緩，凸半球微結構表面升溫速率最快，凹半球微結構表面次之，平面最弱。

圖7 顯示了疏水E-case4 和E-case5 工況下液氬層溫度隨時間變化趨勢，兩種工況下溫度曲線均有轉折點，也在起始沸騰時間點附近。E-case4 工況下三種表面上兩個階段的變化趨勢比較相近，凹、凸半球納米結構表面有溫度交叉，在4.5ns之前顯示出＞＞的溫度排序；而在4.5ns之后，如插圖P2-Concave-E-case4所示情況，由于納米凹槽初始狀態有液氬原子，三種表面溫度變化相差不大。而對于接觸角為132.72°的E-Case5工況，開始階段凸半球納米結構表面溫度升高領先，而由于凹、凸槽內液氬原子蒸發逃逸，使得平面上液氬層的溫度升高變快，超過凹、凸槽上液氬溫度。

圖7 疏水壁面溫度隨時間變化情況

圖8顯示了接觸角在150°～180°的超疏水E-case6和E-case7 工況下液氬層溫度隨時間的變化趨勢，三種表面上溫度曲線轉折也在起始沸騰時間點附近，與前面快照圖4相符。這種超疏水情況，整體溫升情況都比較差，加上初始時刻凹、凸半球納米結構表面便存在空隙（氣泡），直至計算完成，最高溫升在115K 范圍內。超疏水表面出現的氣泡（氣體層）初始時刻便存在，不是加熱過程形成的，在強疏水表面上固-液間作用力非常小，液氬難以完全浸潤壁面，該氣體層的出現是固-液弱相互作用力導致，壁面上能量傳遞效率很低，加熱產生氣泡出現較晚。這與“壁面越親水傳熱效果越好氣泡出現越早”是不同的。比較親水、輸水表面上的溫升最大相差45K，并且平面上溫升大于含微結構表面上的溫升，呈現出>>。加熱過程結束，液體溫度只是從85K提高到100～110K，再一次說明了超疏水表面不適用于強化傳熱過程。

圖8 超疏水壁面溫度隨時間變化情況

2.4 不同工況下氬原子熱流密度變化

熱流密度參考了Liu 等的做法，使用氬原子總能量的變化與液態氬薄膜的橫截面積的商來獲得熱流密度，即= ?(?·)。式中，和分別表示整個系統氬原子能量增量和系統加熱方向的截面面積，選取的截面積均為底面面積。將熱流密度隨時間的變化曲線記錄在圖9～圖11 中。對于光滑的強潤濕性表面，利用本文模擬程序，按照前人的成果中的相應參數進行模擬驗證，得到的熱流密度趨勢與上述文獻一致，且熱流密度峰值約為0.00305eV/(ps·nm)和0.00312eV/(ps·nm)，與 兩 種熱流密度峰值0.003eV/( ps·nm)和0.00306 eV/( ps·nm)相差分別為1.67%和1.96%，驗證了本文模擬程序的準確性。

圖9 所示為親水E-case1、E-case2、E-case3的熱流密度變化曲線，熱流密度變化分為三部分，曲線有兩次轉折變化，統計發現第一個拐點均發生在起始沸騰時間點之前，第二個拐點在時間方面沒有總結出明顯特征，但是都在液膜離開固壁2nm以上出現，這個問題本文不做深入討論。

圖9 親水工況熱流隨時間變化

親水E-case1至E-case3工況下，含凹半球納米結構表面上熱流密度峰值出現最早，但與凸半球納米結構表面工況中峰值出現時間相差不大，平面工況中熱流密度峰值出現最晚，這說明存在納米結構能夠強化換熱，強化效果與納米結構有關；親水E-case1、E-case2、E-case3工況下隨潤濕性減弱，最大熱流密度降低，并且需要花費更多時間到達最高熱流密度點，相應的下降熱流密度也表現出相似性質。

疏水工況的熱流變化如圖10 所示，疏水Ecase4 工況下熱流密度上升與下降的拐點也發生在起始沸騰時間點之前，發生在4ns 附近，直至固-液完全分離后，起始沸騰時間點后，熱流密度才陡然下降。說明在疏水E-case4 下，即使產生氣泡，發生固-液完全分離，液氬層距離固壁位置也不遠，不像親水表面那樣產生較大能量推動液膜遠離壁面。而工況E-case5初始時刻，液膜與凹陷處就已分離，6.0～7.0ns期間由于底部加熱才發生固-液分離現象，這歸因于固-液間弱作用力不能完全被液氬浸潤，疏水表面在加熱初始時，凹陷納米結構存在的氣泡核將會阻礙熱量傳遞，導致疏水凹半球納米結構平面的固-液分離時間大幅度延遲。

對于超疏水E-case6和E-case7工況(圖11)，熱流密度小于疏水工況E-case5，也出現上升和下降過程，含凸半球微結構表面表現有異常，在5.875ns 時便出現了氣泡波動，但直至11.435ns 時才發生固-液分離現象，根據氣泡質心計算，沒有計算得到氣泡起始沸騰點。對于超疏水E-case6和E-case7 工況，由于強疏水性導致壁面對氬原子束縛力很小，升高的溫度對其稍加擾動，便會產生氣泡，出現的氣體核反而阻礙了熱量傳遞，導致液膜上層溫度較低，無法產生足夠的啟動壓力，導致固-液分離現象出現較晚。

圖11 超疏水工況熱流隨時間變化

2.5 影響傳熱因素及其機理分析

計算過程中發現，勢能分布與固-液界面的潤濕性相關，故提取超親水和超疏水兩個工況下初始時刻的勢能部分展示在圖12中，為了對比明顯，將標度統一到-0.12～0eV；統計壁面附近原子的數密度并將其記錄在圖13 中。可以看到潤濕性較強的壁面的勢能顯著大于潤濕性弱的壁面。同樣壁面親水性強會吸附更多的原子；凹半球納米結構近壁面處的氬原子數密度更大，界面處更多的流體分子參與熱傳遞會降低Kapitza 阻值，因此，親水性凹槽里原子排布要密集，傳熱性質更優。在親水表面上的凹半球納米結構表面傳熱性能要優于其余兩個表面；凸半球納米結構表面的勢能絕對值最小，但其接觸面積相對平面更大，近壁面的數密度相對平面仍占據優勢，所以凸半球納米結構表面的傳熱性能也超過了光滑平面。而對于疏水壁面，因其勢能較小，壁面附近原子更容易掙脫束縛，導致可進行熱傳遞的原子數量減少，因而疏水表面熱阻大，傳熱效果差。

圖12 初始階段勢能分布

圖13 初始階段密度分布

圖14 展示了出現氣泡時的溫度隨著接觸角的變化趨勢。圖14(a)中溫度指的是出現氣泡時的壁面溫度，可以看到，隨著浸潤性的增強，在接觸角小于150°內氣泡產生的溫度呈上升趨勢，說明接觸角小即潤濕性強的固-液表面氣泡產生時間要短，這與模擬文獻是一致的，這個微觀結論與宏觀上“疏水表面容易產生氣泡”的結論正好相反，宏觀上連續介質理論不能簡單應用于微尺度系統，在微尺度系統中，宏觀上可以忽略的界面熱阻，在微納尺度變得不可忽略，因此使用溫度階躍對氣泡出現的溫度進行了修正，當考慮固-液界面熱阻時，氣泡產生的溫度與接觸角的關系表示在圖14(b)中，其展現出與圖14(a)相反的趨勢，即疏水壁面出現氣體層需要的過熱度反而更低，這與前人的實驗結果相符。也就是說，在考慮溫度階躍和界面熱阻的情況下，微尺度模擬結果與宏觀實驗本質上是一樣的，疏水表面形成氣泡并發生沸騰現象需要的溫度更低。親水表面固-液相互作用較強，拖延氣泡出現，但是其較強的傳熱性能，傳遞給液體較多的熱量，加速了氣泡產生，而疏水表面固-液相互作用較弱，可促進氣泡產生，傳熱性能卻較差，氣泡產生是兩種作用相互博弈的結果，因微納尺度上界面熱阻不容忽略，導致與宏觀現象的不同。

圖14 出現氣泡時的溫度隨著接觸角的變化

3 結論

采用分子動力學方法模擬計算了不同潤濕性平面和含凹、凸半球納米結構三種表面上液氬層的池沸騰傳熱過程，分析了氣相原子數、氬膜層溫度、熱流密度等參數的變化，通過固-液勢能和接觸角隨溫度變化分析親疏表面產生氣泡的機理，得到以下結論。

（1）平面和含凹、凸半球納米結構三種表面在親水工況下，液氬層的池沸騰傳熱過程可以明顯分為兩個階段，即起始沸騰時間前后，表現為固-液分離前的加熱和固-液分離后的加熱沸騰過程，含凹半球納米結構表面上起始沸騰時間、氣相原子數、溫升最高點及熱流密度均大于相同潤濕性條件的含凸半球納米結構表面，兩種微結構表面均大于平面情況。

（2）對疏水表面上微結構的影響不同于親水表面，固-液間潤濕性差，微結構的存在反而影響傳熱，起始沸騰時間延后，氣相原子數減少，溫升和熱流密度降低，甚至在超疏水工況出現平面上液膜溫升及熱流密度均大于微結構表面的情況。

（3）從固-液作用勢能和界面熱阻討論了微納尺度親疏水表面的微結構沸騰氣泡產生的機理，親水表面上固-液相互作用強，但是固-液間界面熱阻小，而疏水及超疏水工況由于固-液相互作用小，固-液界面有空腔，造成界面熱阻大，傳熱性能差，產生氣泡晚于親水表面。