尺寸效應對微通道內固液界面溫度邊界的影響*

張龍艷 徐進良 雷俊鵬

(華北電力大學, 低品位能源多相流與傳熱北京市重點實驗室, 北京 102206)

(2018 年10 月19日收到; 2018 年11 月19日收到修改稿)

采用非平衡分子動力學方法模擬不同浸潤性微通道內液體的傳熱過程, 分析了尺寸效應對固液界面熱阻及溫度階躍的影響. 研究結果表明, 界面熱阻隨微通道尺寸的變化可分為兩個階段, 即小尺寸微通道的單調遞增階段和大尺寸微通道的恒定值階段. 隨著微通道尺寸的增加, 近壁區液體原子受對側固體原子的約束程度降低, 微通道中央的液體原子自由移動, 固液原子振動態密度近似不變, 使得尺寸效應的影響忽略不計.上述兩種階段的微通道尺寸過渡閾值受固液作用強度與壁面溫度的共同作用: 減弱壁面浸潤性, 過渡閾值向大尺寸區域遷移; 相較于低溫壁面, 高溫壁面處的過渡閾值更大. 增加微通道尺寸, 固液界面溫度階躍呈單調遞減趨勢, 致使壁面溫度邊界和宏觀尺度下逐漸符合. 探討尺寸效應有助于深刻理解固液界面能量輸運及傳遞機制.

1 引 言

近年來, 隨著微納技術的高速發展, 微納設備中流體傳熱傳質機理的研究日益受到學術界的廣泛關注[1-3]. 當熱流流過固液界面時, 界面熱阻的存在導致界面處發生溫度不連續現象. 界面熱阻的影響程度可由熱阻長度LK來衡量(即Kapitza長度其中為固液界面處的階躍溫度,為界面處流體側的溫度梯度. 實驗測量和分子模擬的研究結果都表明, 在液固相互作用較弱的界面, 熱阻長度可達數十納米[4-7]. 在常規尺度下, 固液界面處的溫度邊界條件直接采用固壁邊界, 不考慮溫度階躍. 但是在微納結構中,微尺度效應的存在使得能量輸運及傳遞規律明顯區別于常規尺度[8,9], 導致各種物理參量與系統尺寸有著緊密聯系. 因此, 研究微納系統中界面熱阻及界面溫度階躍隨尺寸的變化關系具有重要的物理意義與實際應用價值.

現有研究表明, 界面熱阻和溫度階躍受到眾多因素的影響, 包括固液相互作用強度[10-13]、固體結構[14]、表面粗糙度[15]、壁面溫度[16,17]、流體剪切率[18]等. 在數值模擬方面, 分子動力學方法作為探測微納尺度現象的有效手段, 被廣泛應用于界面熱阻的研究. Liang和Tsai[19]探究了液膜厚度小于1.3 nm時的界面傳熱過程, 發現固液界面熱阻隨液膜厚度的增加而急劇增大. Wang等[20]模擬體系中狹縫尺寸為0.361—1.62 nm, 得到與上述類似的規律. 而Kim等[21]卻發現當固壁間距為3.24—12.96 nm時, 固液界面熱阻并未發生明顯的變化. 盡管關于尺寸效應對界面熱阻的影響已有一定數量的研究報道, 但是受計算資源的限制, 目前的模擬體系尺寸多數限制在10 nm的數量級. 在實驗研究方面, 原子力顯微鏡等微納測量技術的發展, 促使界面熱阻的實驗測量更加精細化. 目前比較先進的實驗方法包括TTR (transient thermo reflectance)技術和3ω法[22], 其測試材料的最小尺寸在100 nm數量級.

綜上所述, 現有數值模擬研究體系與實驗測量之間存在較大的尺寸跨度, 該范圍內的界面熱阻變化情況很少涉及, 關于溫度邊界的研究仍處于空白. 此外, 對于小于10 nm的微通道內, 已有模擬結果中, 界面熱阻隨微通道尺寸的變化趨勢存有爭議, 觀點尚未統一. 因此, 針對尺寸效應如何影響固液界面熱阻及溫度階躍有待進一步探究. 本文利用非平衡分子動力學方法模擬在給定溫差下不同浸潤性微通道內液體傳熱的過程, 分析了尺寸效應對固液界面熱阻及溫度階躍現象的影響, 并深入剖析界面熱阻作用機制. 研究結果為從物理意義角度揭示溫度邊界條件的“宏觀”與“微觀”尺度的內在機理提供了有益啟發. 本模擬采用開源分子動力學模擬軟件LAMMPS實現.

2 物理模型及模擬細節

圖1為物理模型和界面熱阻長度定義示意圖.系統沿著x和y方向均采用周期性邊界條件,z方向為固定邊界條件. 模擬體系尺寸為(h為變化量,σ為氬原子的尺寸參數), 研究表明, 在x和y方向進一步增大尺寸對計算結果沒有影響. 固體壁面原子Pt按照面心立方(FCC)晶格排列, 初始密度為對應的晶格常數為1.15σ, 其(100)晶面與液體原子接觸, 共有4800個固體原子. 上下壁面每一側的固體原子共有13層, 靠近流體區域的10層固體原子構成熱墻, 與液體區域進行熱量傳輸. 固體板之間充滿液體氬. 初始液氬原子也按照FCC晶格排列方式布置, 初始密度為對應晶格常數為1.72σ, 共有2366—10829 (隨著h發生變化)個液體原子.

液體氬原子之間的相互作用采用Lennard-Jones (LJ)勢能模型, 表達式為

式中r為原子間的距離; 液體氬原子之間的尺寸參數為能量參數為原子質量固體鉑原子之間的相互作用同樣采用LJ勢能模型, 尺寸參數為能量參數為其中下標s表示固體. 液體原子與固體原子的勢能參數均來自Nagayama等[23]的研究報告.

圖1 (a) 物理模型圖; (b) 界面熱阻長度定義Fig.1. (a) Physical model of system; (b) definition of thermal resistance length at interface.

固液之間的勢能作用通過對LJ勢能模型進行修正[24]:

采用Velocity-Verlet算法求解運動方程, 時

為了控制固體墻的溫度, 在熱墻的最外層增加3層虛擬原子, 以模擬其無限大勢場. 虛擬原子類似熱墻的外延擴展. 最外兩層虛擬原子為固定原子, 運行過程中使其固定在晶格位置處保持原子靜止不動, 以便維持系統的穩定. 靠近熱墻的虛擬原子用來模擬與熱墻鏈接的虛擬熱源, 熱浴的實現采用Langevin方法[25], 對其額外施加一個隨機力和阻尼力, 并使其滿足三維郎之萬方程. 這兩個力共同確保實現固體墻及此層固體原子的溫度恒定. 其動量方程表示為

式中pi為第 i個原子的動量矢量;為阻尼常數; f(t)為原子之間的相互作用力; F(t)為隨機作用力. 可以從具有零均值和標準偏差的高斯分布中進行簡化, 其標準偏差為在本文中, 考慮固體原子Pt的簡諧振動模型, Pt原子的相互作用力為簡諧力, 認為固體原子均在其晶格固定位置附近做簡諧振動. 彈簧系數為k = 46.8 N/m. 在模擬中實際固體原子與虛擬原子之間靠虛擬彈簧相連, 在水平方向彈簧系數為0.5k, 在垂直方向為2k. 虛擬原子與固定原子之間的彈簧系數在水平方向為2k, 垂直方向為3.5k.

3 結果與討論

3.1 溫度和密度分布

微通道內液體分布狀態影響固液界面能量傳遞, 引起固液界面溫度邊界條件的變化. 為探索微通道尺寸效應對溫度邊界的作用機制, 首先需要分析微通道內液體的溫度及密度分布特征. 本文為了觀察統計液體的密度和溫度分布, 將液體區域沿著z方向劃分為n層(n隨微通道尺寸H發生變化), 第 l個切片層 (1≤ l ≤ n)在 JStart—JEnd時間段內的無量綱粒子密度

式中Nk為第l層的液體原子數目為各液體分層的厚度; A為液體計算區域在xy平面的面積,

第 l個切片層 (1 ≤ k ≤ n)在 JStart—JEnd時間段內的溫度為

圖2 固液勢能參數 α =0.14 時, 流體沿著z方向的溫度分 布 和 密 度 分 布 (a) H =8.07σ ; (b) H =57.65σ ;(c)H=103.77σFig.2. Temperature and density profile of liquid along z direction when solid-liquid potential energy parameter α=0.14: (a) H =8.07σ ; (b) H =57.65σ ; (c)H=103.77σ.

觀察圖2液體的密度分布特征, 分別從近壁區和通道中央兩個方面進行分析. 一方面, 近壁區液體原子受固體壁面作用較強而呈現類固體排列方式, 液體密度表現出明顯的分層現象. 當微通道尺寸H=8.07σ時, 近熱壁區內液體無量綱原子密度分布的第一個波峰值約為冷壁區約為當微通道尺寸增大到H=57.65σ時,液體原子密度波峰值明顯降低, 熱壁區約為冷壁區約為近壁區液體原子振蕩區間寬度約為 6.34σ. 對于微通道H=103.77σ, 熱壁區液體無量綱密度波峰值2.358, 冷壁區近壁區密度振蕩區間寬度約為6.27σ. 比較上述兩種大尺寸微通道內流體近壁區分布微觀狀態, 發現改變微通道尺寸, 對近壁處振蕩液體區的范圍及密度波動幅度的影響幾乎忽略不計. 另一方面, 對于小尺寸微通道H=8.07σ, 通道中央的液體原子排布仍呈現出微弱的振蕩現象, 這是由于兩側固壁原子的共同作用引起的; 而對于大尺寸微通道H=57.65σ及 1 03.77σ,固體壁面原子對微通道中央區域液體原子的約束程度微弱. 液體原子在微通道內可以自由移動, 原子均勻排布, 其平均密度均約為0.79σ3.

綜上所述, 微通道尺寸的變化引起內部流體溫度和密度分布的改變. 為了分析固液界面傳熱與微通道尺寸的關系是否受其他因素的影響, 進一步探究了壁面浸潤性的作用機制. 圖3為微通道尺寸為H=46.12σ時, 不同固液勢能參數下液體無量綱密度和溫度分布. 圖中黑色虛線表示固液界面無溫度階躍時液體的理論溫度分布. 固液勢能參數決定了固體原子與液體原子的相互作用強度, 它不但改變近壁區液體原子的分布狀態, 還會影響固液界面處的能量傳遞. 當固液勢能參數為α=1.0 時, 固液原子相互作用最強; 當α=0.14 時, 固液原子相互作用最弱. 由圖3(a)的密度分布可以看出, 在近壁區流體排布結構呈現出類固體的層狀排列方式, 當勢能參數α=0.14 時, 其熱壁面處的四層流體數密度峰值依次為1.413 , 1.036 和 0.876 ,根據已有文獻資料[26]可知, 當溫度為100 K (無量綱溫度為 0.826 )時, 晶體氬的密度為1.48 g/cm3(無量綱數密度為ρ?=0.881 ), 對比發現, 近壁處的第一層氬原子密度遠高于結晶密度, 甚至第四層也與晶體密度相近, 因此, 近壁區的氬液體原子表現出類固體的結構特征. 增強壁面浸潤性,近壁區氬原子無量綱數密度波動幅度增大, 即其類固體排列特征更加顯著. 然而, 密度振蕩區域的范圍近似不變, 均為 Δz? ≈6.34 , 壁面浸潤性的影響甚微. 觀察圖3(b)液體無量綱溫度分布可知, 當α=1.0 時, 固液原子傳遞熱量效率最高, 界面處溫度階躍最低, 尤其在熱壁面處幾乎與固體壁面溫度相等.液體內部溫度分布斜率也接近無溫度階躍情況下的溫度梯度. 當α=0.14 時, 液體溫度曲線變平緩,其明顯偏離無階躍溫度分布狀態. 在熱壁面處無量綱溫度階躍約為0.0796, 是α=1.0 時的4.14倍.

圖3 微通道尺寸 H=46.12σ時流體沿著z方向的(a) 密度分布和(b)溫度分布Fig.3. (a) Density profile and (b) temperature profile of liquid along zdirection when microchannel height H= 46.12σ.

3.2 溫度階躍

為了直觀地描述微通道尺寸效應對固液界面溫度階躍的影響, 圖4給出了固液勢能參數α=0.14, 2.0時, 固液界面溫度階躍與微通道尺寸的變化規律. 由圖4可知, 增加微通道尺寸, 界面溫度階躍隨之減小, 微觀尺度下固液界面處的階躍溫度邊界逐漸向宏觀尺度的無階躍溫度邊界靠攏.尤其, 強浸潤性壁面α=0.2 在微通道尺寸H=103.77σ時, 熱壁面處無量綱溫度階躍為0.01175, 界面處液體溫度與壁面溫度的比值兩者幾乎接近. 此外, 界面溫度階躍還與壁面浸潤性和壁面溫度有關. 在壁面溫度相同的情況下, 液勢能參數α=0.2 的相互作用力更強, 其溫度階躍總小于α=0.14 的值. 當壁面浸潤性相同時, 熱壁面處的溫度階躍總小于冷壁面的值. 這是因為溫度升高, 近壁區的液體原子具有更高的能量掙脫固體原子的束縛, 從而與固體壁面碰撞更加頻繁, 能量交換能力增強[27]. 因此, 在高溫壁面處界面熱阻較小, 固液界面溫度階躍相應的較低.

圖4 固液界面溫度階躍隨微通道尺寸的變化Fig.4. Variation of temperature jump at solid-liquid interface with microchannel height.

根據液體的溫度分布, 給出液體實際溫度梯度與微通道尺寸H和冷熱壁面處界面熱阻長度LK,c,LK,h的理論表達式為

此式表明在界面熱阻長度LK不隨通道尺寸變化的情況下, 微通道尺寸H越大, 實際溫度梯度越接近理論溫度梯度, 即在無溫度階躍的情況下的液體溫度梯度. 圖5展示了實際溫度與理論溫度的接近程度隨微通道尺寸的變化關系, 其中表示溫度梯度, 圖中縱坐標表示液體實際理論溫度梯度與理論溫度梯度的差值比理論溫度梯度值. 由圖5可知, 在小尺寸微通道H=5.19σ時, 實際溫度梯度比理論溫度梯度幾乎小1/2. 增大微通道尺寸, 液體實際溫度分布與接近理論狀態之間的區別逐漸減小. 綜上所述, 隨著微通道尺寸的增大, 固液界面溫度階躍呈單調遞減趨勢, 內部流體實際溫度梯度與理論梯度無限逼近, 致使壁面溫度邊界和宏觀尺度下逐漸吻合.

圖5 液體實際溫度與理論溫度的接近程度隨微通道尺寸的變化關系Fig.5. Degree of closeness of actual and theoretical temperature of liquid versus microchannel height.

3.3 界面熱阻

在宏觀尺度下, 無階躍溫度邊界被認為是理所當然的; 而微尺度下, 由于界面熱阻的存在, 固液界面處發生溫度不連續現象, 使得微尺度下溫度邊界條件明顯區別于常規尺度. 如果從物理意義的角度區分宏觀與微觀的溫度邊界條件界限, 那么就需要考慮兩者之間是否具有內在的聯系, 并剖析界面熱阻對溫度邊界的作用機理. 探討尺寸效應對固液邊界熱阻的影響, 將有助于加對強溫度邊界條件物理本質的理解.

圖6為不同浸潤性下, 冷熱壁面處界面熱阻長度隨微通道尺寸的變化. 由圖6(a)可知, 對于較強浸潤性壁面, 即固液勢能參數α＞0.3 , 當微通道尺寸小于10.95σ的尺寸范圍時, 增加微通道尺寸,界面熱阻將隨之略微增大. 對于較弱浸潤性壁面,即勢能參數α=0.2 , 0.14的情況, 界面熱阻隨微通道尺寸的變化更加顯著, 其在小尺寸范圍內迅速上升. 直至微通道尺寸分別達到12.68σ, 23.06σ時,界面熱阻的變化與微通道尺寸近似無關, 尺寸效應可以忽略不計. 總結不同浸潤性下微通道尺寸效應對固液界面熱阻的影響規律, 可將其歸納為兩個部分: 小尺寸范圍的單調遞增趨勢與大尺寸范圍的恒定值規律. 上述兩種規律的尺寸過渡閾值隨壁面浸潤性的減弱向大尺寸范圍遷移. 觀察圖6(b), 可以發現冷壁面處的界面熱阻變化規律與圖6(a)中熱壁面的一致, 區別在于微通道尺寸過渡閾值對壁面浸潤性的依賴程度減弱, 尤其對于弱固液相互作用, 即勢能參數α=0.14 的情形, 其微通道尺寸過渡閾值降低為17.29σ, 相比熱壁面, 約小于5.77σ,產生明顯的差異. 此外, 保持微通道尺寸不變, 界面熱阻隨壁面浸潤性增強而減小; 而且對于相同浸潤性及微通道尺寸的情況, 熱壁面處界面熱阻總是低于冷壁面處界面熱阻值, 該結論與文獻[28]結果一致.

圖6 固液界面熱阻隨微通道尺寸的變化 (a) 熱壁面; (b) 冷壁面Fig.6. Variation of thermal resistance at solid-liquid interface with microchannel height: (a) Hot wall; (b) cold wall.

固液界面熱阻主要受固體與近壁區流體之間傳熱的影響[19-21], 因此對其機理分析需從近壁區液體原子與固體結構特征出發. 由微通道中流體密度分布特征可以看出, 在近壁區液體原子尤其是第一層液體原子按類固體方式排列, 由固體壁面向氬薄膜傳遞熱量的過程可以認為是載熱子的傳遞過程. 在小尺度下, 微通道中的氬薄膜受兩側固體壁面的影響, 呈現出顯著的類固體排列結構特征. 如圖2(a)所示, 隨著通道尺寸的增大, 近壁區氬原子受對側壁面的影響減小, 導致氬原子排布的無序程度增加, 使得其類固體結構特征減弱. 因此, 在熱量傳遞過程中, 增加納米通道的高度, 載熱子的邊界散射增強, 通過彈道輸運的載熱子相應地減少.此外, 增加微通道尺寸, 引起短程載熱子不能彈道輸運穿過界面, 從而增強了散射程度, 使得彈道效應進一步減弱[29,30]. 由此推斷得到, 界面熱阻隨微通道尺寸的增加而增大. 當微通道尺寸增加到一定程度, 近壁區氬原子受對側固體原子的影響忽略不計, 其類固體排列結構幾乎保持不變. 該結論可由圖2(b)和圖2(c)密度分布得到驗證, 在大尺寸通道中, 近壁區氬原子密度分布的振蕩幅度及振蕩范圍變化甚微, 原子排列的有序程度幾乎保持不變.載熱子穿過界面時的散射程度不再隨尺寸變化, 界面熱阻維持不變. 改變壁面浸潤性, 如圖3(a)所示, 固液勢能作用越強, 則近壁區氬原子排列越有序, 從而削弱了載熱子的界面散射, 即增強了彈道效應. 因此, 當微通道尺寸相同時, 界面熱阻隨壁面浸潤性的增強而減小.

綜上所述, 界面熱阻由載熱子的散射程度決定. 目前, 對于界面熱阻的微觀機理解釋主要有聲失配理論(AMM)和散射失配理論(DMM)等. 根據上述理論, 界面熱阻是由相鄰兩種材料的原子之間振動耦合, 即振動態密度(VDOS)失配程度決定的. 固液原子之間的聲子VDOS失配程度越低,更多的能量將以簡諧振動的方式傳遞, 界面熱阻將越小[31]. 聲子VDOS可由速度自相關函數經過傅里葉變換得到:

式中v為原子速度;〈〉表示對不同時間點的統計平均;w為振動頻率. 因此固體與近壁區液體原子之間的VDOS失配程度可表示為

圖7 固液勢能參數為 α=0.14 時固液原子VDOS分布Fig.7. VDOS profile of solid atoms and liquid atoms located at interfaces when solid-liquid potential energy parameter α=0.14 .

計算得到不同微通道尺寸下的固體與近壁區液體原子VDOS分布如圖7所示. 由圖7可知, 固體壁面原子的VDOS的波峰分布在高頻區域, 而近壁區液體原子的VDOS向低頻區域遷移. 這是因為相較于固體原子而言, 液體原子之間的結合力更弱. 對于小尺寸微通道, 內部液體受兩側固壁原子的約束作用顯著, 通道中液體的密度振蕩特征明顯, 彈道聲子傳輸占主要地位, 界面處以簡諧振動的方式傳遞能量, 固體與液膜之間的VDOS失配程度降低, 導致固液熱阻減小. 增大微通道尺寸,處于微通道中央的液體原子幾乎不受固體原子影響, 能夠在微通道內自由移動, 近壁區液體原子VDOS受對側壁面固體原子的束縛甚微, 固液原子振動耦合程度保持不變. 當微通道尺寸為H=34.59σ與H=103.77σ時液體原子VDOS分布幾乎沒有差異, Δ VDOS 未發生變化; 而H=8.07σ時,固體原子VDOS未發生改變, 液體原子的VDOS波峰值明顯下降, 使得固液之間 Δ VDOS 減小.因此, 界面熱阻長度隨微通道尺寸先增大后保持不變.

4 結 論

采用非平衡分子動力學方法模擬不同浸潤性微通道內液體的傳熱過程, 對液體密度和溫度分布展開研究, 分析了尺寸效應對固液界面熱阻及溫度階躍的作用機制, 得到如下結論.

1)界面熱阻隨微通道尺寸的變化規律可劃分為兩個階段: 小尺寸微通道的單調遞增階段, 大尺寸微通道的恒定值階段. 對于小尺寸微通道, 液體原子受固體原子影響顯著, 導致固液原子VDOS失配程度較低, 可以有效減小固液界面熱阻. 對于大尺寸微通道, 只有近壁區液體原子受同側固體原子的約束明顯, 遠離固體壁面的液體原子則自由移動, 固液原子VDOS幾乎保持不變, 導致固液界面熱阻不再隨著微通道尺寸的變化而發生改變.

2)雖然界面熱阻隨微通道尺寸的變化趨勢與壁面浸潤性及壁面溫度無關, 但是上述兩種階段的微通道尺寸過渡閾值卻受二者的影響: 減弱壁面浸潤性, 過渡閾值向大尺寸區域遷移; 相較于高溫壁面, 低溫壁面處的過渡閾值更小.

3)固液界面溫度階躍隨微通道尺寸的增加呈減小趨勢, 使其逐漸接近宏觀尺度下的無階躍溫度邊界條件. 探討這種尺寸效應有助于深刻理解微觀和宏觀尺度下的溫度邊界條件的作用機制, 為微納設備的制造提供理論支撐.