受限納米通道內氣泡核化的分子動力學模擬

苗瑞燦，張石重，陳占秀，楊歷

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

近年來，微納米科學技術在不同領域得到廣泛應用，由于工質核化相變可以有效提高流體的換熱特性，將工質核化相變原理應用到微納尺度中已經成為廣大學者研究的熱點問題。流體在相變換熱時會在通道內生成氣泡并形成兩相流。因此，研究微尺度下的流體相變及氣泡成核生長過程及其機理是提高流體強化換熱的一種重要方法。

固體壁面潤濕性對氣泡核化生長有重要影響，由經典核化理論可知，液體氣泡核化沸騰過程與固液接觸面粗糙度、液體過熱度、固液間的接觸角有一定的關系。固液壁面接觸角增加，氣泡核化沸騰所需要的過熱度越低并且氣泡發生沸騰所需時間越短。Bourdon等[1]利用實驗對壁面潤濕性影響大尺度沸騰初始核化情況進行研究，發現固體與壁面間潤濕性越強，液體初始核化溫度越高。Xu等[2]通過范德瓦爾斯理論研究沸騰核化在異質壁面形成過程，得到在異質情形下壁面潤濕性越弱，氣泡核化半徑越大，更容易形成膜態沸騰。Jo 等[3]研究發現在壁面浸潤性一致的情況下，異質核化相較于均質核化更容易發生沸騰現象。在常規尺度下，以流體為連續介質作為理論依據來研究固液之間傳熱傳質機理，而在微納尺度下流體以微觀粒子為依據來研究固液之間的傳熱傳質過程，固液間界面效應明顯，導致微尺度下流體相變過程較常規尺度下不同。在微尺度下，已有學者對核化沸騰問題進行研究。Mao 等[4]利用分子動力學研究氣泡在均質壁面核化過程，發現微尺度下的核化速率、極限過熱度、核化速率與經典核化理論結果相比有較大差異。Kinjo等[5]研究有限尺寸的納米通道內氣泡核化生長過程，結果發現強親水壁面形成均質核化過程，中等親水壁面形成氣泡異質核化過程以及疏水壁面形成氣膜過程。Nagayama等[6]利用分子動力學研究壁面不同浸潤性對氣泡核化位置的影響，結果發現固液間作用力較強時，流體在通道內發生均質核化，氣泡在通道中間的主流液體內生成，固液間相互作用力較弱時，流體在通道內發生異質核化，氣泡在固體表面處形成，固液間作用力為超疏水時，在固體表面會形成一層氣膜，而沒有氣泡生成。Hens等[7]發現超疏水壁面上很難形成氣泡，而親水性壁面由于較強的固液相互作用，在近壁面處容易形成類固體層，為氣泡核化提供了有利條件。Yamamoto等[8]認為流體與固體之間相互作用力大小影響界面間能量與熱量傳遞，壁面潤濕性較強時，熱量更容易由固體壁面傳到流體中，更有利于核化現象的產生。She 等[9]利用分子動力學方法模擬壁面帶有空腔的通道內氣泡核化生長過程，發現與光滑壁面相比，帶空腔壁面在近壁面處有較大的密度梯度、排斥力以及勢能梯度，從而更有利于氣泡核化生長。Chen等[10]利用分子動力學模擬壁面存在納米凸起對氣泡核化生長的影響，得到納米凸起結構可以強化氣泡核化，但存在最佳凸起高度。張龍艷等[11-12]利用分子動力學研究微納通道內氣泡核化生長過程，發現壁面潤濕性較弱時，氣泡在壁面凹槽內形成，壁面潤濕性較強時，氣泡在通道液體主流區形成。同時發現壁面潤濕性越強越有利于氣泡核化。此外，廣大學者對微納尺度下通道內流體傳熱傳質過程及機理也有不同觀點。鄭文秀等[13]研究壁面不同潤濕性對沸騰核化的影響，發現壁面潤濕性越強，氣泡形成位置離壁面越遠，此時界面熱阻較大不利于核化發生。Novak 等[14]利用分子動力學研究液體在固體壁面處的氣泡核化過程，發現壁面潤濕性越弱，氣泡在壁面處發生異質核化的時間越短。

綜上所述，微納尺度界面效應對流體氣泡核化生長及其相變過程有重要影響，學者們在微通道流體傳熱傳質過程及機理分析中持有不同觀點，且較少學者涉及受限通道內流體核化生長研究，而生活中大部分應用工程都是有限通道內流體傳熱傳質相變過程，因此研究微尺度下受限通道內流體傳熱相變過程就顯得尤為重要。本文利用分子動力學方法對受限納米通道內流體氣泡核化過程進行模擬，研究不同程度親疏水性壁面對氣泡核化生長的影響，分析固體壁面浸潤性對氣泡核化生長的作用機制。本模擬采用開源分子動力學模擬軟件LAMMPS 實現，原子位型實現可視化采用VMD軟件。

1 物理模型及模擬細節

圖1 物理模型圖

采用分子動力學研究了受限通道內流體核化生長過程，模擬系統如圖1 所示，模擬體系尺寸為Lx×Ly×Lz=17.04nm×10.224nm×3.408nm，模擬體系在x、z方向上采用周期性邊界條件，y方向上采用固定邊界條件。固體壁面原子初始狀態按照面心立方晶格（FCC）結構排列，固體壁面共有13714 個原子。固體壁厚h為1.704nm，流體高度H為6.816nm。壁面最外兩層粒子固定不動，以維持模擬體系穩定，其余壁面粒子在初始位置做簡諧振動。流體原子初始狀態也按照照面心立方晶格（FCC）結構排列，系統初始溫度設為85K，其流體飽和密度ρsat為1.401×103kg/m3。設置流體初始密度為0.7ρsat，通道內共有5760個流體原子。

流體原子與流體原子間相互作用力采用Lennard-Jones勢能來描述，見式(1)。

式中，ε=0.01042eV；σ=0.3405nm；r為相互作用的兩個原子之間的距離；rc為截斷半徑，文中取rc=0.85nm；流體原子質量M=39.948g/mol。固體與流體之間采用修正的L-J勢能模型[15]，見式(2)。

固體原子間相互作用參數εs=0.4096eV，σs=0.23377nm，Ms=64g/mol，根據Lorentz-Berthelot公式計算能量參數εls，尺寸參數σls，計算如式(3)、(4)。

固體與固體原子之間相互作用力采用更加精確的EAM 嵌入原子勢來進行模擬，文中使用的銅EAM 勢能由Mishin 等[16]在2001 年提出，計算如式(5)、(6)。為了更好地觀察計算統計結果，在模擬過程中將通道高度平均分成m層，第t層（1≤t≤m）在第(JF-JE)時間段內流體溫度計算見式(7)。

式中，Fi為粒子i的嵌入能；φij(rij)為粒子i和粒子j間的對勢能；ρi為第i原子處的電子云密度；rij為粒子i和粒子j間的距離；kB為玻爾茲曼常數，數值為1.380649×10-23J/K；m為單個流體質量；i、j為第i、j個流體粒子；a為通道內流體x、y、z三個方向，在系統達到平衡狀態后uy、uz的值為0。

流體勢能通過將通道高度平均分成n份，通過計算每層的平均勢能，得到流體勢能在通道內的整體分布，具體計算如式(8)。

式中，rij為粒子i、j之間的距離；uij為兩粒子間的勢能。

在模擬過程中對系統首先采用正則系綜（NVT），維持系統溫度為85K，運行50萬步系統達到平衡，之后設置上壁面溫度為80K，下壁面溫度為140K，流體原子采用微正則系綜（NVE），系統模擬運行150萬步，觀察通道內流體分布情況。模擬過程中系統采用Velocity-Verlet 算法求解運動方程，時間步長為1fs。

2 模擬結果分析

2.1 氣泡生長核化過程

通過調節修正的L-J 勢能函數中的α、β值來控制壁面潤濕性，α、β的值越大，流體與固體間相互作用力越強，壁面潤濕性越好。圖2(a)～(e)為不同壁面潤濕性下氣泡核化生長過程。由圖可見壁面潤濕性不同，加熱情況下流體在通道內有不同的表現形式。當固液勢能參數α=1、β=1 時，流體與固體間作用力最強，流體被強作用力吸附在近壁面處，受熱初期小氣泡在液體薄層上生成，通道內發生均質核化，隨受熱時間增加，小氣泡在液體薄層上方逐漸融合形成一層氣膜，氣膜推動通道內液體向上運動。由于通道上壁面溫度為80K，且通道高度有限，系統平衡后氣膜會在近熱壁面處呈片條狀，如圖2(a)的2000ps快照圖所示。當固液勢能參數α=0.14、β=0.14 時，如圖2(e)所示，近壁面處流體與固體間排斥力最強，此時壁面呈超疏水性，流體從受熱初期到終了都不曾生成氣泡。當固液勢能參數α=0.6、β=1，α=0.14、β=1，加熱初期流體在壁面液膜層生成較小的氣泡，隨時間增加小氣泡逐漸融合形成氣膜，推動流體向上移動，隨著α的減小，氣膜尺寸逐漸減小近壁面處液膜厚度也在減薄。當固液勢能參數α=0.14、β=0.6 時，流體與固體間有較強排斥力，初期流體在固體壁面處生成小氣泡，流體在通道內發生異質核化，隨時間增加氣泡尺寸逐漸增大，最終氣泡大小基本保持不變。

圖2 不同壁面潤濕性下氣泡核化生長過程

圖3(a)、(b)中分別給出了勢能參數為α=1、β=1，α=0.14、β=0.14 初始快照圖，在近壁面處已將“類固體”層局部標出。“類固體”層是通道內流體由于壁面勢能的相互作用力在近壁面固體層處流體粒子形成和固體壁面粒子類似的排列結構，這種由流體粒子形成的具有和壁面固體類似的結構層，成為“類固體”層。對比兩個不同勢能參數下的“類固體”層快照圖，可以看出勢能參數為α=1、β=1時，近壁面處流體排列方式更加整齊規則，結構近似趨于壁面粒子排列方式。勢能參數為α=0.14、β=0.14時，在近壁面相同位置處流體排列方式與圖2 相比，整齊度和規則度不高，“類固體”層形成程度較小。這種在近壁面處流體形成的規則結構程度對流體核化生長有一定的促進作用。

圖3 不同壁面潤濕性下氣泡核化生長過程及初始快照圖

圖4是系統平衡后不同勢能參數條件下發現流體質量密度沿通道高度Ly方向分布圖，由圖可見在近壁面處流體與固體壁面存在一層液體層，隨固液勢能參數值減小液體層的密度逐漸降低，近熱壁面處液膜厚度也在逐漸減薄。在通道高度Ly方向上2.5～4nm 位置，除固液勢能參數α=0.14、β=0.14外，流體質量密度均處于最低值，此處為氣泡生成位置，由密度值可看出隨壁面潤濕性減弱氣泡尺寸在逐漸減小。當固液勢能參數為α=0.14、β=0.14時，流體與固體間排斥力強，近壁面處流體原子受較強排斥力作用，會在近壁面處形成一層流體原子極稀少的氣膜層，因此會出現近壁面流體原子質量密度幾乎為零的情況，由于固液間強排斥力作用流體原子均勻分布在通道3～7nm 位置處，很少存在于近壁面處。

圖4 系統平衡后流體質量密度沿Ly方向分布

圖5 是模擬時間為600ps 時流體質量密度沿Ly方向的分布圖。圖中給出了不同壁面潤濕性下同一時刻的流體質量密度分布，由圖可知系統模擬相同時間后固液勢能參數α=1、β=1 時，在氣泡核化位置處流體密度最小，固液勢能參數α=0.14、β=0.6時，在氣泡核化位置處流體密度最大，隨著固液勢能參數α、β增加，通道內氣泡核化尺寸逐漸增大。固液間較大勢能參數使流體在近壁面處有較強的流固作用力，壁面上吸附的流體粒子數越多質量密度越大，形成“類固體”層就越厚，熱量更容易通過“類固體”層傳給通道內流體，氣泡越容易核化。因此，表面潤濕性強的固體壁面更有利于流體發生氣泡核化。

圖5 時間為600ps時流體質量密度沿Ly方向的分布

圖6是固液勢能參數α=0.6、β=1不同時刻流體質量密度沿Ly方向的分布圖。由圖可知初始時刻通道內流體密度維持在1.0kg/m3左右，近壁面處流體密度相對較大，流體在近壁面處形成液膜薄層，隨著壁面受熱時間增加，液膜厚度逐漸減薄，在2.5～4nm 位置流體密度顯著下降，流體在此位置發生均質核化生成氣泡如圖3(a)所示。此外，流體發生均質核化過程主要在受熱初期階段，500～730ps 是流體在通道內核化的主要階段，氣泡在此時間段內生成，之后隨受熱時間增加氣泡尺寸不在發生變化，小氣泡會逐漸融合形成氣膜，氣膜推動流體向上運動。

圖6 勢能參數α=0.6、β=1流體質量密度分布

2.2 氣泡核化機理

圖7 不同壁面潤濕性下流體溫度分布

通過分析不同勢能參數下流體溫度變化。近一步研究壁面潤濕性對流體傳熱過程的影響。將系統模型沿Ly方向平均分成50 層，計算每層流體粒子的平均溫度，只統計流體區域溫度變化。圖7給出了系統平衡后不同壁面潤濕性下流體溫度分布情況。由圖可見，固液勢能參數值越大近壁面處流體溫度越高，即壁面潤濕性越強，固液傳熱效率越好，因此，流體在潤濕性強的壁面處獲得的能量多，更容易達到核化條件形成氣泡。勢能參數α=1、β=1時，在Ly為2.8nm位置有溫度突變，由于系統穩定后氣泡融合形成氣膜[如圖2(a)2000ps時的快照圖]，氣膜層處含有極少數流體原子，溫度從固體壁面傳到氣膜層再經氣膜層傳給通道上層流體，由于氣膜層傳熱效率低，極少熱量通過氣膜層傳給上層流體，因此溫度會在有氣膜層的地方出現突變。

圖8 給出了勢能參數為α=0.6、β=1，α=0.14、β=1不同時刻流體溫度沿Ly方向變化。由圖可知壁面溫度經過固液接觸面傳給通道內流體，靠近熱壁面處的流體溫度率先升高，隨通道高度Ly方向增加流體溫度逐漸降低。此外，流體溫升主要發生在1000ps 之前，之后隨時間增加溫度幾乎保持不變，說明在受熱初期流體在近熱壁面處發生核化后，氣泡核化尺寸逐漸趨于平穩，隨加熱時間增長溫度也達到平衡。由圖8(a)、(b)對比發現，在800ps 時，勢能參數為α=0.6、β=1，近壁面溫度已經接近140K，隨時間增加近壁面處溫度保持不變。勢能參數為α=0.14、β=1的情況下，近壁面處溫度隨時間增加逐漸升高，最終近壁面處溫度接近132K，說明在相同時間內，勢能參數值越大，熱量越容易傳遞給通道內流體，即潤濕性強的壁面界面熱阻小，固液間傳熱效率高，氣泡核化時間短。

圖8 不同勢能參數、不同時刻流體溫度沿Ly方向分布

圖9為勢能參數α=0.14、β=1時x、y、z方向平均受力分布圖，由圖可知，在x、z方向上流體粒子受力幾乎為零，說明流體粒子在x、z方向沒有優先運動，因為x、z方向設置為周期性邊界條件。在近熱壁面處觀察到y方向上力先從0.0058下降到零，之后又從零增加到0.003，最后從0.003一直減小到零左右，力的數值為負代表力的作用方向為y軸負方向，在通道中心力的大小維持在零左右，近冷壁面處力從0 增加到0.003 左右。觀察發現在近壁面處有力的梯度出現，在通道中心處幾乎沒有力的梯度，并且熱壁面處力的梯度要大于冷壁面處力的梯度。近壁面處流體所受勢能大小相同，但力的梯度卻不同，說明高溫導致流體在近壁面處表現出較大的力度梯度，從而流體更容易出現核化現象。圖9 是不同壁面潤濕性下y方向受力分布圖，根據圖9可知流體粒子在x、z方向上受力幾乎為零，所以圖10只給出不同壁面潤濕性情況下y方向上的受力情況。由圖可知，不同壁面潤濕性下，流體粒子在近熱壁面處力的梯度都要大于冷壁面附近力的梯度，在近熱壁面處，不同潤濕性壁面力的梯度雖然不同，但相差很小。固液勢能參數越大，流體在近壁面處越容易出現“類固體”現象，更容易束縛在近壁面不易進行自由運動，因此近壁面處會出現隨勢能參數值下降，流體粒子在y方向受力增加的現象。

圖9 勢能參數為α=0.14、β=1時x、y、z方向受力分布

圖10 不同壁面潤濕性下y方向受力情況

圖11 不同壁面潤濕性下流體勢能分布

圖11 是不同壁面潤濕性下流體在通道內勢能分布情況，由圖可以看出，流體在通道中心勢能幾乎都保持在0.05左右，在靠近壁面處，流體勢能出現較大差異，當勢能參數α=1、β=1 時，近壁面處流體勢能最大，約為0.25。當勢能參數α=0.14、β=0.6 時，近壁面處流體勢能最小，幾乎接近于零。近壁面處勢能大小隨α、β值的增加而增加，流體在近壁面處表現的勢能強度導致壁面不同程度的潤濕性，近壁面處流體勢能強度越大，壁面原子對流體原子吸附力越強，壁面潤濕性越強。由于壁面勢能參數不同，流體在近壁面處表現出不同的現象，當勢能參數最強時，流體會在壁面處形成一層液膜，流體排列方式出現“類固體”形式，壁面通過“類固體”將熱量傳遞給通道內流體，傳熱效果較好。當勢能參數較小時，在近壁面處由于壁面粒子與流體粒子間有較強排斥力，導致在近壁面處出現一層氣膜，壁面熱量透過氣膜將熱量傳遞給流體，傳熱效果相對較差。因此，當勢能參數較小時，熱量很難傳遞到流體中，通道內流體將不會成核形成氣泡。

3 結論

采用分子動力學方法模擬了流體在受限通道內沸騰核化生長過程，研究了壁面潤濕性對流體成核過程的影響機理，結果如下。

（1）固體壁面潤濕性不同，會對通道內流體成核過程產生很大影響。壁面潤濕性強，固體與液體間勢能相互作用力大，固液間有強大的吸附力，流體原子在近壁面處會出現“類固體”排列形式，熱量由壁面通過“類固體”層進一步傳遞給通道內流體，流體傳熱效果最好，在通道內流體受熱容易出現均質核化現象。

（2）壁面潤濕性相對較弱時，流體與固體間作用強度相對較小，在近壁面處不會形成“類固體”層，熱量由固體壁面直接傳給通道內流體，熱量傳遞速度相對較慢，傳熱效果較好，流體在通道內會出現異質成核現象。

（3）當壁面潤濕性最弱時，流體與固體間有較大的排斥力，在近壁面附近會出現一層氣膜，熱量由壁面透過氣膜層傳遞給流體，傳熱效果較差，流體在通道內不會出現成核現象。