混合潤濕性柱狀納米結構對銅板上納米氬膜沸騰傳熱的影響

劉峰瑞, 陳占秀, 李源華

(河北工業大學 能源與環境工程學院, 天津 300401)

1 引 言

近年來，隨著微電子技術的發展，系統設備的體積和面積越來越小，各種微納器件的熱管理成為亟待解決的核心問題，需要盡可能高效地強化熱交換. 在沸騰相變傳熱過程中，液體需要吸收大量的汽化潛熱，增加了傳熱速率. 因此，高效率、高傳熱率的沸騰相變熱迅速成為各類學者研發的重點[1-3]. 沸騰相變行為和傳熱效率受到許多因素的影響，例如流體性質、壁性質和外部條件. 其中，壁面潤濕性和壁面納米結構起主導作用.

對于沸騰傳熱，表面納米結構的幾何形狀和尺寸，作為刺激氣泡成核生成的構型[4]，對核態沸騰過程有重要影響. 研究人員通過修改加熱表面的拓撲結構，例如周期性正弦表面、納米通道結構表面和錐形納米結構表面，模擬了納米結構對爆炸沸騰的影響，發現與平面相比，它們顯著增強了從固體表面到液體的熱量傳遞，并促進了相變過程中的傳熱[5-7]. 此外，Wang和Zhang等人[8，9]通過分子動力學研究發現，納米結構對成核沸騰的初始成核時間和初始溫度有顯著影響. 與光滑表面相比，納米結構表面的臨界熱通量顯著增加，氣泡成核增強，初始成核時間縮短. Liu等[10，11]發現粗糙結構使沸騰起始時間提前，蒸發質量通量增大. 并提出隨著表面積的增加，界面吸附的水分子越多，單位面積的相互作用能就會增加，從而降低界面的熱阻. 在這項工作中，Liu等人[12]研究了內表面粗糙度對凹面納米結構沸騰傳熱性能的影響. 結果表明，傳熱性能強烈地依賴于粗糙元的結構. 低固體含量配置具有更好的傳熱性能，并可延遲沸騰過程的開始. 傅等人[13]進行了分子動力學模擬，以研究納米結構對被熱銅板突然加熱的水的快速沸騰的影響. 研究表明，納米結構對固體銅片到液態水的能量傳遞和液態水到蒸汽的相變過程有顯著影響. 仍有一些關于納米結構對固/液界面傳熱的影響的研究[14-18]，這也表明結構化表面可能是爆發沸騰的主導機制之一. 潤濕性對氣泡成核有很大影響. 周等[19]研究了液氬在不同潤濕性的圖形化表面上的氣泡成核. 發現疏水部分的面積比影響氣泡的成核溫度. 提出了一個最佳的面積分數，以盡量降低氣泡成核溫度. Wu等人[20-21]研究了在不同潤濕性條件下光滑和粗糙銅表面上水膜的快速沸騰. 研究發現，與純潤濕表面相比，混合潤濕表面通過調節蒸汽在加熱銅表面的擴散行為來強化沸騰傳熱，而納米通道高度的增加有利于成核點的形成和蒸發. Diaz等人[22]研究了疏水/親水納米圖案增強沸騰傳熱的功效. 結果表明，親水納米柱和親水性壁的組合表現出最佳的傳熱性能. Shahhardi等人[23]研究了表面潤濕性對沸騰引發、氣泡成核和生長的影響. 觀察到親水表面會加速氣泡的成核，但減緩氣泡的膨脹速度，從而延緩蒸汽膜的形成. Hens等人[24]發現氣泡成核和蒸汽膜的形成取決于固-液界面的潤濕性. 親水表面為氣泡成核和蒸汽膜形成提供了有利條件. 此外，也有研究者[25-27]通過設計混合潤濕性來研究薄液膜的爆炸沸騰過程. 結果表明，混合濕壁更容易引發沸騰現象，對氣泡起始溫度、膜態沸騰溫度和蒸發速率有顯著影響.

上述研究揭示了納米尺度下池沸騰的部分機理. 在傳熱過程中，納米結構和表面潤濕性都起到促進作用，但很少有研究將上述兩個影響因素結合起來進行分析. 本文采用分子動力學方法在納米尺度上模擬了液態氬在銅表面的快速沸騰，構建了具有混合潤濕性的圓柱納米結構壁，然后改變納米結構的高度和疏水部分的比例，觀察其對沸騰傳熱的影響.

2 模擬細節與方法

2.1 物理模型

創建了一個尺寸為91 ?(X)×91 ?(Y)×543 ?(Z)的模擬框，模擬系統的初始配置如圖1所示，系統沿x軸和y軸采用周期性邊界條件，z軸采用固定邊界，模擬體系由固體區域、液體區域和氣體區域三部分構成. 此外，假設模擬系統的頂部邊界是反射的，當氬原子試圖穿過反射墻時，它將被反彈，該過程不會產生能量交換和損失. 模擬工質為銅和氬. 工質選用氬原子是因為氬的性質穩定，勢函數簡單，可以很大程度上節約計算時間. 模擬過程中求解L-J勢能時采用球形截斷處理，截斷半徑選取為10 ?. 氬的晶格常數為5.72 ?. 銅原子呈面心立方 (FCC)晶格結構，晶格常數為3.62 ?. 對于所有壁面，固定層原子保持靜止，由2500個銅原子組成，防止壁面受熱變形. 緊接著是作為系統中的熱源的幻影原子“phantom atoms”[28]，由6250個銅原子組成，其溫度通過恒溫器來維持. 銅原子受到大小為-Kr的彈簧力的約束，其中K是彈性常數，為0.293 eV/nm[29]，r是原子從當前位置到其初始位置所在位置的距離，以便銅原子在整個模擬過程中交換能量時會圍繞其原始位置振動. 將Lennard-Jones(12-6)勢能模型應用于氬原子、銅原子之間的相互作用力的計算.

圖1 模擬系統初始配置示意圖：(a)總體視圖；(b)壁面正視圖；(c)壁面俯視圖Fig. 1 Schematic diagram of the initial configuration of the simulation system：(a)overall view；(b)front view of the wall；(c)top view of the wall

在本文所有模擬中，基底厚度都超過了其截斷半徑距離，這保證了壁面厚度不會影響模擬結果. L-J(12-6)勢能函數表達式：

(1)

其中，rij為計算粒子對之間的間隔距離，在計算氬-氬原子對之間相互作用力時，能量參數ε和尺寸參數σ選擇為：ε=0.0104 eV，σ=0.3400 nm. 當計算銅-銅原子對之間相互作用力[30]時，能量參數ε和尺寸參數σ選取為ε=0.2047 eV，σ=1.9297 nm. 為了描述氬原子和銅原子之間的相互作用，采用Lorentz-Berthelot混合法則[31]進行求解，特征長度由下式計算：

(2)

正如之前的研究[32]所提到的，可以通過改變描述壁面和液態氬原子之間相互作用的 LJ 勢的能量參數來調整潤濕性. 當εAr-Pt=0.5εAr-Ar(0.0052 eV)，獲得親水材料，當εAr-Pt=0.25εAr-Ar(0.0026 eV)，獲得疏水材料. 如圖2顯示了潤濕性占比不同的納米結構壁面，底面黃色和圓柱上深藍色代表親水部分，淺藍色代表疏水部分. 在這項工作中，圓柱納米結構柱寬為1.6 nm，有兩個不同的高度hc，分別為1.5 nm和2.5 nm，其中疏水部分占圓柱比β分別為0、33%、53%和73%.

圖2 潤濕性占比不同的納米結構壁面案例：(a)β=0 (b)β=33% (c)β=53% (d)β=73%Fig. 2 Cases of nanostructured walls with different wettability ratios：(a)β=0 (b)β=33% (c)β=53% (d)β=73%

2.2 模擬方法

模擬初始情況下，將包含15200個氬原子的液膜放置在基板上方. 為了減少模擬計算成本，在每次模擬開始時使用共軛梯度 (CG)算法最小化整個系統的能量，時間步長設置為5fs. 隨后將NVT(恒定原子數、體積和溫度)系綜應用于整個系統，使其在90 K下運行1ns達到平衡狀態. 達到穩定后，改為微正則(NVE，即恒定原子數、體積和能量)系綜. 最后，通過恒溫器將熱源層銅原子的溫度瞬間升至160K，以誘導液膜的爆炸沸騰現象. 所有模擬均使用開源經典MD 代碼大規模原子/分子并行模擬器 (LAMMPS)[33]進行，并使用開放可視化工具 (OVITO)[34]進行可視化.

2.3 模型驗證

模型驗證采用Yin 等人[35]的工作，模擬了液氬在光滑銅表面上的池沸騰過程. 計算時采用文獻[35]中的模擬參數. 計算如圖3所示，將液氬溫度隨時間的變化情況和液氬發生沸騰脫離銅表面的時間與文獻中的結構進行了比較. 本文計算結果與文獻中的結果相比，液氬溫度的最大誤差為 5.9%，發生沸騰時間的誤差為5.0%，兩者計算非常相近證明本文計算結果的正確性.

圖3 驗證結果：(a)液氬溫度隨時間變化；(b)液氬發生沸騰脫離銅表面的時間.文獻[35](左)和當前模擬(右)Fig. 3 Verification results：(a)temperature of liquid argon as a function of time；(b)time for liquid argon to boil off the copper surface. Reference [35](left)and current simulation (right)

3 結果與分析

3.1 固壁上液氬的池沸騰過程

本文給出了兩種圓柱高度情況下不同混合潤濕性表面液氬池沸騰的過程，快照如圖4所示.在這里，β=0表示疏水占比為0，即固體壁面為全親水. 模擬初期，液膜均與底面相接觸，并未因增加疏水占比而浮于圓柱上方. 首先，如圖4(a)當hc=1.5 nm，β=0時，在加熱初期，系統處于純蒸發階段，僅在液-氣界面處發生液膜向蒸汽氬原子的轉變直至13.50 ns時，底部液體層出現氣泡或縫隙，在14.40 ns氣泡聚結形成蒸汽膜，將液體與壁面分開并推動其向上運動，發生沸騰；之后壁面與液膜間的氬原子及液膜中的氬原子仍然繼續被加熱形成氣態氬原子. 對比圖4(b)圓柱高度為2.5 nm，初始時表面蒸發形成的氣態氬原子較多，固液之間形成蒸氣膜的時間提前到了8.30 ns. 其余六個增加疏水占比的表面，它們都有著相似的沸騰過程. 然而，就沸騰發生時間而言，發現相同柱高情況下，β越大發生的越快，這表明傳熱效率更高. 具體來講，Case7到達沸騰起始時間為10.25 ns，而Case5需要的時間相對較長，為11.75 ns. 另一方面，當柱高增加到2.5 nm，β=33%、β=53%和β=73%時發生沸騰的時間分別比柱高為1.5 nm時提前了4.55 ns、5.1 ns和4.4 ns. 我們注意到，增加疏水占比后氣泡并不再只出現于圓柱底部，而是出現在親疏水交界處或疏水部分頂部. 其原因可以推斷為與疏水部分接觸的Ar原子受到其影響，在接下來的非平衡分子動力學過程中，這部分Ar原子不斷吸收熱量并率先突破液體內部分子間作用力，脫離兩原子之間的相互束縛而運動，且距離不斷增大，使氣泡時間提前，促進氣泡形核，加速了爆炸沸騰的過程. 總結液氬沸騰的起始溫度和時間如圖5所示. 從圖中可以看出，對于兩種柱高，隨著疏水占比β的增大，都會降低爆炸沸騰的起始溫度和起始時間. 因此我們得出本文提出的親疏水性混合的圓柱納米結構壁面對爆炸沸騰起著促進作用.

3.2 液氬沿Z方向的密度分布

相變過程的變化也可以用壁面上部的瞬態空間密度變化來描述. 當氣泡核在液膜中產生時，近壁成核區的密度會降低. 為獲得密度剖面，將仿真系統沿z方向分為50個高度為10.84 ?的切片，每1 ns輸出一次每個切片的密度. 液氬在不同柱高混合潤濕性表面上密度剖面的變化趨勢如圖6所示. 由圖可知，柱高為1.5 nm的不同混合潤濕性壁表面經歷9 ns的蒸發階段后，近壁面處液氬數密度大概在0.016 ρσ3，柱高為2.5 nm的不同混合潤濕性壁表面經歷4 ns的蒸發階段后，近壁區液氬數密度大概在0.018 ρσ3. 隨著加熱時間的增加，液體中會產生蒸汽，導致密度下降，近壁面處形成低密度區，距離壁面較遠處氬密度增加同時，產生的蒸汽將一層厚厚的液膜提升到z方向較高的位置，并且由于系統溫度的升高，該液膜的密度峰值也逐漸減小. 我們發現，由于case1沸騰起始時間較晚，長時間的蒸發使其近壁面區氬密度降低到了0.01ρσ3左右，遠壁面處氬密度才開始增加；然而case2、case4、case6和case8隨著β的增大，發生沸騰的時間不斷加快，高密度區出現在系統越來越高的地方. 另一方面，在50 ?左右的位置總是有一個密度讀數，這是由于表面附近吸附氬原子而形成的液態氬薄層.

圖6 液氬在不同柱高混合潤濕性表面上沿z方向的密度分布Fig. 6 Density distribution of liquid argon along the Z-direction on mixed wetted surfaces with different column heights

3.3 液態氬原子蒸發速率

凈蒸發原子數的變化能夠很好的顯示出加熱過程中的相變速率. 本文計算凈蒸發原子數的方法是通過在給定截止距離0.53 nm內的每個氬原子的相鄰原子數，當指定氬原子周圍的相鄰原子數小于7，則認為它是蒸氣原子. 否則，它被認為是液態氬原子. 該方法已在文獻[36，37]中廣泛使用. 觀察圖7，隨著液膜不斷吸收熱量，不同柱高混合潤濕性表面上液態原子均被蒸發，其數量首先在蒸發階段快速減少. 當液體蒸發一段時間后發生沸騰，液態原子隨時間減少趨勢開始變緩，這種趨勢隨著β的增大出現的越來越早，柱高為1.5 nm時，β=0、β=33%、β=53%和β=73%分別在14.5 ns、12 ns、11.8 ns和10 ns附近減少的速率放緩，這正好大致對應于液氬沸騰的起始時間. 正如我們所預想的，在相同的疏水占比情況下，柱高為1.5 nm的壁面上液體原子數均小于柱高為2.5 nm的，如β=53%時，因為柱高增加后的沸騰溫度128 K要低于1.5 nm時的135 K，且在更短時間內發生沸騰. 由于較高的沸騰溫度和更長的沸騰起始時間，使得蒸發過程持續時間更久，氣態氬原子更多，壁面上液體原子數較少.

圖7 不同柱高混合潤濕性表面上液體原子數隨時間變化Fig. 7 Variation of the number of liquid atoms with time on the mixed wetted surfaces with different column heights

3.4 液氬溫度及能量變化

圖8(a)給出了不同柱高混合潤濕性表面液氬溫度隨時間變化情況，觀察到，兩種柱高的傳熱速率和最終溫度呈現出Tβ=0>Tβ=33%>Tβ=53%>Tβ=73%的規律.β=0%時，液氬升溫速率最快，但因其為純親水結構，對液氬具有很強的吸引力，需要更高的溫度才能發生爆炸沸騰現象. 但是當圓柱增加了73%疏水占比后，最先達到沸騰溫度點，發生固-液分離現象并使溫度變化出現拐點. 溫度變化體現了液氬吸收能量的大小，故圖8(b)給出了不同柱高混合潤濕性表面上液氬的動能與勢能變化. 我們知道，體系吸收的動能與勢能之和為體系吸收的總能量，總能量顯示出與溫度幾乎一致的變化規律. 如β=53%時，系統總能量增量先迅速增大，隨后由于固-液分離，吸收能量減緩，并保持較低的增長率直至模擬時間結束，表明液氬仍在與壁面進行能量交換. 此外，柱高相同時，隨著β的增大，能量增量曲線出現拐點的時間提前，即系統達到沸騰所需要的能量隨之減少. 仔細觀察動能勢能變化ΔE*，可以得知吸收的能量主要轉換成了勢能，使其在某刻沖破勢能壁壘，即導致氣泡核的形成和生長.

圖8 不同柱高混合潤濕性表面上液氬溫度及能量變化隨時間變化情況Fig. 8 Variation of liquid argon temperature and energy with time on the mixed wetted surfaces with different column heights

為了進一步研究β對沸騰傳熱的影響，本文計算了不同柱高混合潤濕性表面上的熱通量隨時間變化情況，如圖11所示. 固體和流體原子之間的熱通量由下式計算：

(3)

其中Asub是x-y平面上的截面積. 由于在t=0 ns時表面溫度突然從90 K升高到160 K，而液膜的溫度仍然平衡狀態時的90 K，此時液膜與壁面的溫差達到最大. 因此，不同柱高混合潤濕性表面上的熱通量在初始時刻達到最大值. 隨著傳熱過程的進行，固-液溫差不斷減小，使得熱通量下降. 觀察圖9，β相同時，柱高為2.5 nm的混合潤濕性表面上熱通量曲線出現拐點的時間比柱高為1.5 nm時提前，且熱通量下降時的斜率更大，傳熱更劇烈. 我們發現，出現拐點原因是在發生沸騰時底部出現蒸汽膜，使得傳熱性能惡化，柱高相同時，隨著β的增大，拐點提前出現. 這說明親疏水性混合的圓柱納米結構壁面可以強化換熱速率，從而導致更快的爆炸沸騰.

傳熱系數是反映傳熱速率的一個重要參數. 在本研究中，固/液界面處的傳熱系數HTC(Heat transfer coefficient)計算如下：

(4)

其中Ts是固體表面的溫度，Tl是液態氬的溫度. 混合潤濕性表面上液氬的HTC如圖10所示. 對于hc=1.5 nm，液氬的沸騰時間較晚，由于較大的分子間作用力，純親水表面的HTC在蒸發階段略有增加. 隨著β的增加，疏水部分增加，成核時間變短，液層中出現更多的蒸汽，降低了傳熱性能，液氬的HTC越來越低.hc=2.5 nm時，HTC的變化趨勢與hc=1.5 nm時相似，不同的是，液氬的HTC先略有上升，然后緩慢下降，這是由于沸騰時間加快，固液之間產生蒸汽膜，降低了傳熱效果.

4 結 論

在本文中，采用分子動力學方法構建了疏水-親水混合柱納米結構壁對池沸騰的影響，討論了疏水占比為0、33%、53%和73%時，兩種圓柱高度下爆發沸騰的起始溫度和時間. 此外，還深入研究了具有混合親水性和疏水性的圓柱形納米結構對爆發沸騰的影響. 得出的結論總結如下：

(1)與純親水的納米結構壁面相比，增大疏水部分占比和圓柱高度后可以降低爆炸沸騰的起始溫度并使發生沸騰的時間提前.

(2)通過觀察不同配置壁面的沸騰過程快照，發現增加疏水占比將改變氣泡出現的位置，不再只是從底部產生，而是轉變為親水疏水交接處或者圓柱頂端. 證明混合潤濕性納米結構有利于氣泡成核，促進爆炸沸騰.

(3)當疏水占比相同時，柱高為2.5 nm的混合潤濕壁面的爆炸沸騰起始溫度低于柱高為1.5 nm的混合潤濕壁面. 這是由于增大了與液氬之間的接觸面積，在保持親疏水比例不變的情況下，進一步放大了混合潤濕性納米結構的影響.

(4)隨著疏水性在混合潤濕性納米結構上所占比例的增加，體系爆發沸騰發生的速度加快，使得液氬在相變過程中需要從壁面吸收的能量減少，平均熱流降低，液氬的HTC降低，導致沸騰后傳熱效果不佳.