不同潤濕表面液體沸騰的分子動力學模擬

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

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

電子器件集成度的日趨增加，使散熱問題已經(jīng)對電子器件發(fā)展形成制約。相變過程被認為是解決高熱流散熱的重要途徑[1]。相變表面的潤濕性對相變過程有著重要影響，隨著微納米表面強化沸騰技術(shù)相關(guān)研究的深入進行，表面親疏水性對沸騰影響的研究逐步展開[2]。鄭曉歡等[3]對常壓下光表面、超親水和超疏水表面的池沸騰進行實驗研究，認為超親水表面親水疏氣，產(chǎn)生氣泡體積小但氣泡脫離頻率較快；而超疏水表面親氣疏水，使得激活活化點更為容易。馬強等[4]發(fā)現(xiàn)換熱系數(shù)與表面過熱度有關(guān)，當過熱度較高，親水和超親水表面換熱系數(shù)大于疏水和超疏水表面；但當過熱度較低時，疏水和超疏水表面換熱系數(shù)比親水和超親水表面大，這是由于疏水和超疏水表面有著較低的沸騰起點。在微納尺度條件下，由于尺寸限制實驗難以開展，利用基于經(jīng)典力學原理的分子動力學（molecular dynamics，MD）模擬可以深入了解原子的微觀動力學行為，將有助于理解在納米尺度上相變的機制。已有學者對潤濕性的影響進行了分子動力學模擬研究。Nagayama等[5]利用分子動力學研究了受限納米通道中的氣泡行為，發(fā)現(xiàn)親水表面為均相成核，疏水表面卻為異相成核。鄭文秀等[6]認為潤濕性越弱時，靠近壁面的分子向其他分子的能量傳遞越容易。相反的，張龍艷等[7]通過研究壁面浸潤性對氣泡初始核化和生長的影響，認為壁面浸潤性越強，能量傳遞越快，氣泡在固壁處越容易核化，氣泡生長也更為迅速。Peng等[8]認為界面強度在形成吸收層中具有重要作用，蒸發(fā)液膜的熱傳遞強烈依賴于吸收層。諸多學者對于納米尺寸下的相變過程進行了研究。Chen等[9]通過模擬在具有非均勻潤濕性納米結(jié)構(gòu)表面上的氣泡成核行為，發(fā)現(xiàn)較大面積的強親水區(qū)有利于氣泡核的形成和生長。Hens等[10]認為親水表面為氣泡成核和蒸汽薄膜的形成提供了有利條件。Mao 等[11]研究了在有限空間內(nèi)由熱銅板加熱的液態(tài)水膜的快速沸騰，發(fā)現(xiàn)靠近板的液態(tài)水分子立即過熱并快速相變。Seyf等[12]研究了壁面材料（鋁和銀）的種類對超薄液態(tài)氬膜在納米結(jié)構(gòu)上爆炸沸騰的影響，發(fā)現(xiàn)由于固液相互作用的降低，固體銀壁附近的原子的速度增加更多，導(dǎo)致液滴的上升速度增加，液滴與表面的分離發(fā)生得更快。與之不同，Hasan 等[13]在親水表面觀察到爆炸沸騰，但是對于疏水表面，觀察到擴散蒸發(fā)，且親水表面?zhèn)鳠岣哂谑杷砻娴膫鳠帷?/p>

綜上所述，壁面潤濕性對固壁表面上液膜的沸騰過程有著重要影響。盡管對于納米尺度下相變的研究已經(jīng)取得了諸多成果，但到目前為止，不同學者對壁面潤濕性影響傳熱過程機理并未達成一致，對于由于壁面潤濕性不同帶來的運動狀態(tài)不同的研究涉及較少。因此，需要進一步研究納米級相變傳熱和沸騰現(xiàn)象的行為和機理。本文利用分子動力學對不同潤濕表面液膜沸騰過程進行模擬，研究不同程度親疏水性壁面對沸騰的影響，重點分析不同潤濕壁面沸騰過程中能量轉(zhuǎn)換和液體運動中的差異。分析潤濕性對其影響，有助于更好地理解納米尺度下沸騰與核化的機理，以期望能夠填補實際應(yīng)用理論支撐的缺陷。本模擬采用開源分子動力學模擬軟件LAMMPS 實現(xiàn)，原子位型實現(xiàn)可視化采用OVITO軟件。

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

建立圖1所示模擬體系，對固體壁面及置于系統(tǒng)內(nèi)氬原子進行分子動力學模擬。該模擬體系為Lx×Ly×Lz=34.57σ×227.15σ×6.91σ 的長方體盒子（σ為氬原子尺寸參數(shù)），由固體區(qū)、液體膜和蒸氣區(qū)域組成。系統(tǒng)底部為固體幾何模型，4680 個鉑原子置于模擬系統(tǒng)的底部，以面心立方結(jié)構(gòu)（FCC）排列方式構(gòu)建金屬固體，密度為ρσ3=2.613。固體幾何模型從下至上分別設(shè)置為3部分：將最底面的一層原子固定，防止模擬過程基板發(fā)生變形；中間兩層固體原子設(shè)置為“虛擬原子”，與外界恒溫加熱器相連作為熱源；上部金屬原子是真正的鉑原子，與氬原子進行熱量交換。對除固定原子外的固體鉑原子施加彈簧力，將其約束在初始位置，彈簧系數(shù)為3249.1εσ-2。選擇氬原子作為液體原子，7380 個液態(tài)氬原子被放置于固體鉑表面，密度為ρσ3=0.78，厚度為34.5σ。模擬系統(tǒng)其余部分充滿蒸氣氬原子。x 和z 方向上應(yīng)用周期性邊界條件，在y方向應(yīng)用固定邊界條件。模擬框y+方向頂部設(shè)置反射墻，當氬原子試圖穿過反射墻時，它將被反彈，該過程不會產(chǎn)生能量交換和損失。

圖1 物理模型圖

原子對之間的相互作用，本文采用Lennard-Jones（LJ）勢函數(shù)(1)描述。

式中，下角標i、j可取為l或s，l和s分別代表流體原子和固體原子；ε為能量參數(shù)，取值分別為εll=1（對應(yīng)的實際能量參數(shù)為1.6×10-21J）和εss=50εll（對應(yīng)的實際能量參數(shù)為8.36×10-20J）；σ是尺寸參數(shù)， σll=1 （對應(yīng)的實際尺寸參數(shù)為0.3405nm），σss=0.73σll（對應(yīng)的實際尺寸參數(shù)為0.2475nm）。上述參數(shù)均來自文獻[7]的報告。r是兩個原子之間的距離，取截斷半徑為3.5σAr-Ar，超過截斷半徑時忽略其原子間作用力。

根據(jù)Lorentz-Berthelot 公式計算能量參數(shù)εls、尺寸參數(shù)σls，計算公式如式(2)、式(3)。

為獲得不同潤濕性的壁面，使用參數(shù)α對固液間勢能進行調(diào)節(jié)，見式(4)。

表1 列出了模擬過程中不同粒子對之間的L-J參數(shù)和對應(yīng)的接觸角，本文中Case ε7.07表示固液作用為7.07的表面，其余兩種類推。為了確定不同表面的潤濕能力，如圖1(b)和圖1(c)所示，模擬了6nm 直徑氬液滴在不同的潤濕性基底上的擴散過程，以觀察接觸角。模擬過程在NVT系綜下進行，溫度控制為90K，模擬時間為2ns，以確保液滴充分在表面上蔓延。

采用velocity-Verlet 算法來求解整個粒子運動方程，時間步長為0.001τ，其中τ=(mσ2/ε)1/2為LJ約化時間常數(shù)。整個模擬過程分為3步：首先，使用共軛梯度（CG）算法實現(xiàn)能量最小化，以確保模型的初始狀態(tài)是穩(wěn)定的。采用正則系綜（NVT）控制系統(tǒng)在T=0.744εκB-1（相對應(yīng)實際溫度為90K）下運行50 萬步，在模擬過程中監(jiān)測整個系統(tǒng)的溫度、壓力和能量變化，確保系統(tǒng)在NVT 過程結(jié)束時達到平衡狀態(tài)。然后撤掉整體熱浴，僅對基板進行控溫，將基板的溫度提高到T=1.57εκB-1（相對應(yīng)實際溫度為190K），此溫度能夠保證液體沸騰的同時達不到固體熔點。當將固體壁面加熱至所研究的狀態(tài)后，將T=0.744εκB-1的氬與T=1.57εκB-1高溫壁面在微正則系綜（NVE）進行集合并采集所需數(shù)據(jù)，在此時從0時刻開始重新計時。

表1 不同潤濕性表面參數(shù)

2 模擬結(jié)果分析

2.1 液膜沸騰過程密度變化

疏水性和親水性通常與液體為水時的潤濕性現(xiàn)象有關(guān)，本文中雖然液體不是水，但遵循相同的慣例，仍用疏水性和親水性表示，它們實際上意味著疏液性和親液性。為明確顯示液體薄膜在加熱過程中的動態(tài)行為，圖2給出了不同潤濕表面上模擬系統(tǒng)運動過程快照。與前人的工作結(jié)果[10-13]一致：在最初階段，液態(tài)氬在固體鉑板表面被加熱；隨著加熱的進行，氬原子開始蒸發(fā)；在加熱進行一定程度時，近固體表面處的液體開始成核出現(xiàn)氣泡并逐漸擴展成密度較小的氣體層［圖2(a)］，而距壁面較遠的液氬分子仍為液體狀態(tài)，近壁面處氣體層推動未蒸發(fā)的液氬上升，出現(xiàn)固液分離現(xiàn)象。在液體區(qū)域呈現(xiàn)出3個部分：近壁面氣體層，其上的液體區(qū)域以及系統(tǒng)最上面的蒸氣區(qū)域。ε=7.07的表面，在150τ 左右時出現(xiàn)了分離現(xiàn)象，液體區(qū)域向上運動至頂部與反射墻碰撞被彈回，由于近固體表面蒸氣產(chǎn)生的壓力導(dǎo)致液體薄膜重新向上運動，在運動過程中液體層幾乎保持不變，而ε=0.50 表面液體層在上升過程出現(xiàn)了明顯縮小過程。隨著潤濕性的減弱，固液分離的發(fā)生時間延后，ε=1.00 表面與ε=0.50表面分別在900τ～1000τ和4800τ～4900τ之間出現(xiàn)了固液分離現(xiàn)象?；诒疚呐c文獻[10]在模型以及加熱過程中出現(xiàn)的現(xiàn)象上的相似性，故認為本文結(jié)果是可靠的。

圖2 不同潤濕性壁面上沸騰過程

為了更進一步了解氬在y方向上運動情況和密度分布情況，將模擬體系沿著y 方向平均分成200個切片，計算在每個切片的原子數(shù)以獲得每個切片的平均數(shù)密度。不同時刻y方向的數(shù)密度變化顯示在圖3中，縱坐標為原子的量綱為1數(shù)密度ρσ3，橫坐標為y方向上的層數(shù)。隨著加熱時間的增加，近壁面處形成了低密度區(qū)，距離基板較遠處氬密度增加。系統(tǒng)中液體氬的密度峰值逐漸減小，體積也逐漸縮小，縮小速率隨著潤濕性減弱而增大。Case ε7.07 液膜體積縮小速率極為緩慢，液體厚度和密度峰值基本不變，Case ε0.50 峰值密度降低速率更明顯，意味著液體蒸發(fā)更快。Case ε7.07 在1300τ時數(shù)密度大于900τ，是因為液體運動到頂部受到擠壓使得密度增大。同時應(yīng)當指出，初始階段壁面附近液氬數(shù)密度略高于液體整體密度，這是由于固液間相對作用力較強時，壁面附近會吸附氬原子形成“類固體”[14-15]。

圖3 不同模型沿y方向氬數(shù)密度分布

圖4 不同壁面潤濕性下近壁面氬數(shù)密度分布

整個過程中壁面附近3.5σAr-Ar內(nèi)氬原子數(shù)密度如圖4所示，可以發(fā)現(xiàn)壁面附近氬原子數(shù)密度隨著加熱的進行呈顯出先降低后穩(wěn)定的趨勢，潤濕性強的壁面附近氬原子數(shù)密度下降較快，達到近似穩(wěn)定狀態(tài)也較早；隨著潤濕性減弱，下降速率也隨之降低，達到近似穩(wěn)定需要時間更長。從圖2(a)可以看出，壁面附近始終存在氬原子，根據(jù)潤濕性的不同，3 種表面上氬原子數(shù)密度分別穩(wěn)定在0.453、0.136和0.102附近。隨著親水性的減弱，固液之間作用力降低，壁面附近氬原子數(shù)密度逐漸降低。

2.2 沸騰過程中能量變化

圖5 不同壁面潤濕性下流體溫度變化

圖5(a)給出了不同壁面潤濕性下流體溫度隨時間變化規(guī)律，3種潤濕表面氬溫度均呈現(xiàn)先迅速增加后緩慢增加趨勢，Case ε7.07 潤濕性最強，溫度快速上升最快，升溫速率出現(xiàn)拐點時間最早，最終溫度最低，隨著潤濕性減弱，升溫速率隨之減弱，溫度上升變緩的時刻也隨之延后，但最終的平均溫度更高。這說明潤濕性較強Case ε7.07在模擬開始傳熱更好，溫度上升更快，發(fā)生固液分離較早；發(fā)生固液分離后氣體層的存在阻礙了熱量傳遞，使得溫度上升速率減緩。隨著潤濕性減弱，初始階段傳熱效果變差，溫度上升減緩，發(fā)生固液分離時間延后。

圖5(b)給出了加熱開始120τ時不同壁面潤濕性下流體溫度沿厚度方向變化規(guī)律，此時氣體層尚未出現(xiàn)?？梢园l(fā)現(xiàn)，靠近壁面處溫度并不一致，潤濕性最強的Case ε7.07近壁面溫度也最高，潤濕性降低，溫度隨之降低；液氬層遠離固體壁一側(cè)的溫度近似一致。通過圖5(b)可以看出，潤濕性越強，壁面與流體間的溫降越小，傳熱效果越好；潤濕性強，液體中的溫度梯度更大，近壁面氣體更容易達到沸騰條件產(chǎn)生氣體層，出現(xiàn)固液分離現(xiàn)象。

圖6 不同壁面潤濕性下流體能量變化

溫度變化體現(xiàn)出氬原子吸收能量的差異。圖6(a)給出了不同壁面潤濕性下流體能量隨時間變化規(guī)律，能量增量指每個時間步長上從固體基板流入流體的能量，通過統(tǒng)計所有氬原子的每原子勢能和動能來獲得能量增量?？梢钥闯?，能量增量與潤濕性的關(guān)系和溫度類似，先迅速增加，然后由于固液分離，能量交換變慢。此外，在模擬后期流體能量以較低的增長率保持增長，表明基板和液態(tài)氬之間仍有能量交換，這解釋了為什么圖5中溫度能夠繼續(xù)上升。Case ε7.07 潤濕性最強，能量上升最快，上升變緩時間最早，最終能量最低，隨著潤濕性減弱，能量增加速率隨之減弱，但最終能量增量更大。結(jié)合圖6(b)，可以得知吸收的能量主要轉(zhuǎn)換為了勢能。

不同潤濕性壁面相變過程在時間上有很大差異，為探究差異的原因，計算了熱通量的變化。通過氬原子總能量的變化除以液態(tài)氬薄膜的橫截面積來獲得熱通量，即q=(1/A)×(?EAr/?T)，式中，EAr和A分別表示整個系統(tǒng)氬原子能量增量和系統(tǒng)y方向的截面面積。Cao 等[16]同樣采用了這種方法計算了潤濕性對傳熱的影響。不同壁面的熱通量變化展示在圖7中，不同表面的平均熱通量在0τ達到最大值，然后隨時間逐漸減小。這是因為初始時刻固液之間溫差是最大的，隨著傳熱過程的進行，固液溫差減小，使得熱通量減小，在150τ 左右的時候，ε=7.07表面熱通量開始小于ε=1.00表面熱通量，在約1000τ 的時候，ε=1.00 表面熱通量開始小于ε=0.50表面熱通量。初始階段親水表面的平均熱通量明顯高于疏水表面，這說明固液間相互作用力越強，傳熱速率越快。這種傳熱速率上的差異導(dǎo)致了潤濕性強的表面發(fā)生爆炸沸騰的時間要早于潤濕性弱的表面。

圖7 不同壁面潤濕性下熱通量隨時間變化規(guī)律

2.3 沸騰過程中液膜運動

氬原子平均質(zhì)心位置隨時間變化趨勢展示在圖8中，以便更為全面地了解整個過程運動情況。在ε=7.07 表面，在0～150τ 時間內(nèi)質(zhì)心在y 方向以一個較緩慢速度上升，在150τ 左右上升速度突然增加，直至上升到215σ，此時液氬層已到達模擬體系頂部。這與圖3 中液體氬原子先緩慢受熱膨脹、然后突然爆炸沸騰的現(xiàn)象具有一致性。由于頂部反射墻的存在，液氬原子運動方向在碰撞后發(fā)生了改變，開始向下運動，在向底部運動一定距離后，再次向上運動，做往返運動。對于ε=1.00 表面，氬原子質(zhì)心上升到1580σ開始第1次下降，此時液氬層尚未達到頂部，而ε=0.50 表面上，在計算時間內(nèi)并未發(fā)現(xiàn)質(zhì)心有明顯下降。此外，對于不同潤濕性ε=7.07 與ε=0.50，質(zhì)心能夠到達的位置分別為215σ 和115σ。液體質(zhì)心隨時間變化的斜率可反映其上升速度，潤濕性強的ε=7.07 表面上升速度要高于潤濕性較弱的表面。

圖8 不同壁面潤濕性下氬質(zhì)心變化情況

圖9給出了整個過程上部氣體層和下部氣體層壓力隨時間的變化情況，P 和PL 分別用來指代液體層上方和下方的氣體的壓力。由于氣體層出現(xiàn)時間不同，下層壓力起始點并非0τ?？梢钥闯觯瑢τ讦?7.07 表面，其下部壓力遠大于且一直大于上部氣體壓力，這可以解釋液氬的迅速上升行為以及碰撞后液體層能夠再次上升的原因；隨著潤濕性減小，在初始上升階段的壓差也隨之減小。

ε=7.07 和ε=1.00 兩種表面從固液分離到液氬首次改變運動方向時間段內(nèi)壓力變化如圖10所示，ε=7.07表面下部壓力要遠大于上部壓力，隨著液體向上移動，下部壓力在一定時間內(nèi)呈下降趨勢，上部壓力增大；ε=1.00 表面上壓力出現(xiàn)了交叉，在1200τ～1600τ 時間區(qū)間內(nèi)，下部氣體壓力要低于上部氣體，推測是由于蒸發(fā)使得上部壓力升高而初始階段下部壓力較小。這種不同可以用于解釋ε=7.07表面上的液氬有著很高的上升速度且到達頂部發(fā)生碰撞才改變方向，而ε=1.00 表面首次改變方向時并未到達頂部。

圖9 不同壁面潤濕性下壓力變化情況

圖10 不同壁面潤濕性下流體壓力

3 結(jié)論

通過分子動力學模擬研究固體基板上氬的相變過程，分析了壁面潤濕性對沸騰過程中能量傳遞和液體運動情況的影響，得出了以下結(jié)論。

（1）在近壁面附近存在未蒸發(fā)液體原子，其數(shù)密度隨著潤濕性增加而增加。3種表面均發(fā)生了固液分離現(xiàn)象，出現(xiàn)固液分離的時間隨潤濕性增加而減小。固液分離后液體氬的密度、體積逐漸減小，減小速率隨潤濕性增加而降低。

（2）潤濕性影響基板與液體間傳熱過程。潤濕性強，能量交換速率快，熱通量大，升溫迅速且液體中溫度梯度大。氣膜的產(chǎn)生會阻礙傳熱，模擬結(jié)束時壁面潤濕性越強，系統(tǒng)內(nèi)流體溫度越低。氬原子吸收的能量主要轉(zhuǎn)換為勢能。

（3）出現(xiàn)固液分離現(xiàn)象后，下部氣體壓力大于上部氣體壓力，推動液體層上升。潤濕性越強，上下壓差越大。在ε=1.00 表面，會在一段時間發(fā)現(xiàn)上部氣體壓力較大的情況。液體質(zhì)心能夠達到高度隨著潤濕性增強而增加。