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

（河北工業(yè)大學(xué)能源與環(huán)境工程學(xué)院，天津300401）

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

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

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

1 物理模型及模擬細(xì)節(jié)

圖1 物理模型圖

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

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

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

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

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

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

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

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

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

2 模擬結(jié)果分析

2.1 氣泡生長核化過程

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

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

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

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

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

圖4 系統(tǒng)平衡后流體質(zhì)量密度沿Ly方向分布

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

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

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

圖6 勢能參數(shù)α=0.6、β=1流體質(zhì)量密度分布

2.2 氣泡核化機(jī)理

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

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