納米氬液膜蒸發(fā)行為的分子動力學(xué)模擬研究

王寶和 ， 李瑞靜 ， 王 維

(1.大連理工大學(xué) 化工學(xué)院 ， 遼寧 大連 116024 ； 2.大連理工大學(xué) 化工機械與安全學(xué)院 ， 遼寧 大連 116024)

?開發(fā)與研究?

納米氬液膜蒸發(fā)行為的分子動力學(xué)模擬研究

王寶和1， 李瑞靜1， 王 維2

(1.大連理工大學(xué) 化工學(xué)院 ， 遼寧 大連 116024 ； 2.大連理工大學(xué) 化工機械與安全學(xué)院 ， 遼寧 大連 116024)

采用非平衡分子動力學(xué)模擬技術(shù)，探討模擬溫度、能量系數(shù)、粗糙度因子、相面積分數(shù)等對納米氬液膜蒸發(fā)行為的影響。模擬結(jié)果表明：納米氬液膜在光滑壁面或納米結(jié)構(gòu)粗糙壁面上的蒸發(fā)，存在恒速蒸發(fā)和降速蒸發(fā)兩個階段；在各個蒸發(fā)階段，蒸發(fā)通量相差不大；蒸發(fā)過程達到平衡后，納米欄柵形粗糙壁面吸附的氬原子數(shù)最多，納米方柱矩陣形的次之，光滑壁面的最少。對于光滑壁面上的納米氬液膜蒸發(fā)，隨著模擬溫度的提高，恒速蒸發(fā)階段的時間變短，蒸發(fā)通量迅速增大；在恒速蒸發(fā)階段，能量系數(shù)對蒸發(fā)通量幾乎沒有影響；在降速蒸發(fā)階段，隨著能量系數(shù)的增加，固體壁面吸附的氬流體分子數(shù)增大。對于納米欄柵形粗糙壁面上的納米氬液膜蒸發(fā)，隨著粗糙度因子或相面積分數(shù)的增大，恒速蒸發(fā)階段的時間增加，蒸發(fā)通量減??；當蒸發(fā)過程達到平衡后，固體壁面吸附的氬原子數(shù)增多。

納米氬液膜 ； 分子動力學(xué) ； 模擬 ； 蒸發(fā)

眾所周知，由于固體壁面上納米薄液膜的厚度處于納米尺度， 使得其蒸發(fā)過程的實驗研究工作難以開展，但隨著計算機技術(shù)的迅猛發(fā)展，采用分子動力學(xué)(MD)模擬技術(shù)，對這樣一個傳遞過程進行研究已經(jīng)成為可能。Yu等[1]利用分子動力學(xué)模擬技術(shù)，探討了脫離壓力對氬液膜蒸發(fā)—冷凝過程的影響，并獲得了納米薄液膜厚度的計算方法。Hens等[2]利用分子動力學(xué)模擬方法，研究了氬納米液膜在固體鉑壁面上的蒸發(fā)和沸騰過程，探究了溫度和固體表面潤濕性的影響。Shavik等[3]通過分子動力學(xué)模擬手段，研究了固—液之間的相互作用勢能以及固—液界面潤濕性(包括疏液性、親液性和中性)對固體鉑壁面上納米氬液膜的蒸發(fā)和爆炸沸騰過程的影響規(guī)律。Wang等[4]借助于分子動力學(xué)模擬方法，探索了納米氬液膜在鋁納米結(jié)構(gòu)粗糙表面上的蒸發(fā)和爆炸沸騰過程。模擬結(jié)果表明：與光滑表面相比，在一定范圍內(nèi)，表面的粗糙納米結(jié)構(gòu)有利于提高蒸發(fā)速率。本文擬通過非平衡分子動力學(xué)模擬方法，研究納米氬液膜在固體壁面上的蒸發(fā)特性，探討模擬溫度、能量系數(shù)、粗糙度因子、相面積分數(shù)等因素對納米氬液膜蒸發(fā)行為的影響。

1 模擬方法

1.1 模擬體系的建立

模擬體系的初始構(gòu)型如圖1所示，模擬盒子尺寸為6.1 nm×6.1 nm×12.8 nm(包括2 nm的刪除區(qū))。3 136個氬原子隨機分布成鋪滿壁面的納米液膜，氬原子的初速度由隨機數(shù)發(fā)生器確定[5]。液膜上方為真空區(qū)域。壁面采用Pt原子參數(shù)[ε(Pt)=8.35×10-20J；σ(Pt)=2.475×10-10m，討論能量系數(shù)影響時除外]，以FCC的方式排布。壁面原子數(shù)量與粗糙度因子或相面積分數(shù)大小有關(guān)。

1.固體壁面 2.納米氬液膜 3.刪除區(qū)

1.2 勢能模型的選取

氬流體分子之間的作用勢能，采用Lennard-Jones(L-J)勢能函數(shù)進行描述，分子i和j之間的相互作用勢能函數(shù)U，如式(1)所示[6]。

(1)

式中：U為勢能函數(shù)，N為模擬分子個數(shù)，rij為分子i和j之間的距離，ε為能量參數(shù)，σ為尺寸參數(shù)；對于氬流體，ε=1.655×10-21J，σ=3.405×10-10m。

固體壁面原子與氬流體分子間的相互作用勢能，依然采用L-J的12-6勢能形式，如式(2)所示。相應(yīng)的尺度參數(shù)采用Lorentz-Berthelot混合規(guī)則來確定，如式(3)所示[7]，相應(yīng)的能量參數(shù)如式(4)所示[8]。

(2)

σLS=(σL+σS)/2

(3)

εLS=ErεL

(4)

式中：rLS為固體壁面原子與氬流體分子之間的距離，σS為固體壁面原子之間的尺度參數(shù)，σL為氬流體分子之間的尺度參數(shù)，σLS為固體壁面原子與氬流體分子之間的尺度參數(shù)，εL為氬流體分子之間的能量參數(shù)，εLS為固體壁面原子與氬流體分子之間的能量參數(shù)，Er為能量系數(shù)。由式(4)可見，可以通過調(diào)整能量系數(shù)Er的大小，來考察不同表面能材料對潤濕性的影響[8]。

1.3 模擬細節(jié)

模擬在x和y方向均采用周期性邊界條件，z方向采用固壁邊界條件。系統(tǒng)采用NVT系綜進行模擬。采用速度Verlet法積分運動方程，時間步長為1 fs。并用Woodcock變標度恒溫法進行恒溫控制。截斷半徑為10。模擬時間為3 ns，前1 ns為弛豫時間，確保1 ns后系統(tǒng)達到平衡狀態(tài)。采用刪除分子法，實現(xiàn)非平衡蒸發(fā)過程，即液相氬原子受熱氣化脫離液膜之后，一旦進入刪除區(qū)(如圖1所示)，模擬系統(tǒng)自動將其刪除。LAMMPS程序設(shè)置每隔0.002 ns，最多刪除20個氬原子[7]。

1.4 納米粗糙壁面的構(gòu)建

本文構(gòu)建的納米欄柵形和納米方柱矩陣形粗糙壁面的結(jié)構(gòu)參數(shù)，分別如圖2和圖3所示。其中，凸起部分的寬度為a，凸起高度為h，凹陷部分的寬度為b。

納米欄柵形和納米方柱矩陣形粗糙壁面的粗糙度因子r1、r2與結(jié)構(gòu)參數(shù)的關(guān)系，如式(5)和(6)所示；其相面積分數(shù)f1、f2與結(jié)構(gòu)參數(shù)的關(guān)系，如方程(7)和(8)所示[7]。

圖2 納米欄柵形粗糙壁面的結(jié)構(gòu)參數(shù)示意圖

圖3 納米方柱矩陣形粗糙壁面的結(jié)構(gòu)參數(shù)示意圖

(5)

(6)

(7)

(8)

1.5 蒸發(fā)通量的確定

蒸發(fā)通量是指單位時間內(nèi)離開固體表面吸附區(qū)域的原子數(shù)除以固體表面的面積。Amsterdam規(guī)則認為，對于氬原子i，若以r=1.5 為半徑的球體范圍內(nèi)存在至少4個氬原子，則可以認為氬原子i屬于液相[9]。在模擬計算過程中，對模擬得到的軌跡文件每隔0.4 ns統(tǒng)計汽、液相原子數(shù)，就可以計算出1 ns離開液膜的氬流體分子總數(shù)，再根據(jù)固體表面積，就可以得到蒸發(fā)通量的統(tǒng)計值。

2 結(jié)果與討論

2.1 納米氬液膜在光滑壁面上的蒸發(fā)

2.1.1 能量系數(shù)的影響

當氬流體分子數(shù)為3 136，模擬溫度為100 K，采用MD模擬計算得到的能量系數(shù)Er=0.4～10時，納米氬液膜在光滑壁面上的液相氬原子數(shù)與時間的關(guān)系，如圖4所示；蒸發(fā)通量與時間的關(guān)系如圖5所示。由圖4和圖5可見，納米氬液膜的蒸發(fā)過程存在恒速蒸發(fā)和降速蒸發(fā)兩個階段。在過程的初期，蒸發(fā)通量近似不變，即為恒速蒸發(fā)階段；然后，隨著過程的繼續(xù)進行，蒸發(fā)通量不斷下降，即為降速蒸發(fā)階段，直至蒸發(fā)通量趨近于零的平衡過程。在恒速蒸發(fā)階段，能量系數(shù)對蒸發(fā)通量幾乎沒有影響；在降速蒸發(fā)階段，隨著能量系數(shù)的增加，固體壁面吸附的氬流體分子數(shù)增多。

圖4 不同能量系數(shù)下的液相氬原子數(shù)隨時間的變化

圖5 不同能量系數(shù)下的蒸發(fā)通量隨時間的變化

2.1.2 溫度的影響

當氬流體分子數(shù)為3 136，能量系數(shù)Er=2時，采用分子動力學(xué)模擬方法，計算得到的不同模擬溫度下，納米氬液膜在光滑壁面上的液相氬原子數(shù)與時間的關(guān)系，見圖6；蒸發(fā)通量與時間的關(guān)系見圖7。由圖6和圖7可見，蒸發(fā)過程可以分為恒速蒸發(fā)和降速蒸發(fā)兩個階段；且隨著模擬溫度的提高，恒速蒸發(fā)階段的時間變短，蒸發(fā)通量迅速增大。

2.2 納米氬液膜在納米粗糙壁面上的蒸發(fā)

2.2.1 固體壁面形貌的影響

由式(5)和(6)可見，當取a=b=h=2個晶格長度時，納米欄柵形和納米方柱矩陣形兩種粗糙壁面的粗糙度因子相等(r=2)。

圖6 不同溫度下的液相氬原子數(shù)與時間的關(guān)系

圖7 不同溫度下的蒸發(fā)通量與時間的關(guān)系

當氬流體分子數(shù)為3 136，模擬溫度為100 K，能量系數(shù)Er=2，納米欄柵形和納米方柱矩陣形兩種粗糙壁面的粗糙度因子都等于2時，采用分子動力學(xué)模擬方法，計算得到的納米氬液膜在兩種形貌壁面上的液相氬原子數(shù)及蒸發(fā)通量隨時間的變化，分別見圖8和圖9。圖中還與光滑壁面的MD模擬結(jié)果進行了比較，可見，納米氬液膜在光滑壁面、納米欄柵形粗糙壁面及納米方柱矩陣形粗糙壁面上的蒸發(fā)過程，都存在恒速蒸發(fā)和降速蒸發(fā)兩個階段；在各個蒸發(fā)階段，蒸發(fā)通量相差不大；但蒸發(fā)過程達到平衡后，納米欄柵形粗糙壁面吸附的氬原子數(shù)最多，納米方柱矩陣形的次之，光滑壁面的最少。

圖8 不同形貌壁面上的液相氬原子數(shù)隨時間的變化

圖9 不同形貌壁面上的蒸發(fā)通量隨時間的變化

2.2.2 粗糙度因子的影響

對于納米欄柵形粗糙壁面，當取a=b時，由方程(7)計算得到，相面積分數(shù)f=0.5；在保持(a+b)=4個晶格長度的情況下，通過改變壁面突起高度h的大小分別為1、2個晶格長度，壁面粗糙度因子r分別為1.5、2.0。

對于納米欄柵形粗糙壁面，當氬流體分子數(shù)為3 136，模擬溫度為100 K，相面積分數(shù)f=0.50，能量系數(shù)Er=2時，采用分子動力學(xué)模擬技術(shù)，計算得到的納米氬液膜在粗糙度因子r分別為1.5、2.0時的液相氬原子數(shù)及蒸發(fā)通量隨時間的變化，分別如圖10和11所示?？梢姡S著粗糙度因子的增大，恒速蒸發(fā)階段的時間增加，蒸發(fā)通量減??；當蒸發(fā)過程達到平衡后，壁面吸附的氬原子數(shù)增多。

圖10 不同粗糙度因子壁面上的液相氬原子數(shù)隨時間的變化

圖11 不同粗糙度因子壁面上的蒸發(fā)通量隨時間的變化

2.2.3 相面積分數(shù)的影響

對于納米欄柵形粗糙壁面，當取h=2.5個晶格長度，(a+b)=5個晶格長度時，由式(5)計算得到，粗糙度因子r=2；在保持(a+b)=4個晶格不變的情況下，變化壁面凸起寬度a分別為1、1.5個晶格長度，使相面積分數(shù)f分別為0.25、0.5。

對于納米欄柵形粗糙壁面，當氬流體分子數(shù)為3 136，模擬溫度為100 K，粗糙度因子r1=2，能量系數(shù)Er=2 時，采用分子動力學(xué)模擬技術(shù)，計算得到的納米氬液膜在相面積分數(shù)f分別為0.25、0.5時的液相氬原子數(shù)及蒸發(fā)通量隨時間的變化，分別如圖12和13所示?？梢姡敶植诙纫蜃右欢〞r，隨著相面積分數(shù)的增大，恒速蒸發(fā)階段的時間增加，蒸發(fā)通量減??；當蒸發(fā)過程達到平衡后，壁面吸附的氬原子數(shù)增多。

圖12 不同相面積分數(shù)壁面上的液相氬原子數(shù)隨時間的變化

圖13 不同相面積分數(shù)壁面上的蒸發(fā)通量隨時間的變化

3 結(jié)論

本文采用非平衡分子動力學(xué)模擬方法，探討了模擬溫度、能量系數(shù)、粗糙度因子、相面積分數(shù)等因素對納米氬液膜蒸發(fā)行為的影響，得到如下結(jié)論：①納米氬液膜在光滑壁面或納米結(jié)構(gòu)粗糙壁面上的蒸發(fā)過程，都存在恒速蒸發(fā)和降速蒸發(fā)兩個階段；在各個蒸發(fā)階段，蒸發(fā)通量相差不大；蒸發(fā)過程達到平衡后，納米欄柵形粗糙壁面吸附的氬原子數(shù)最多，納米方柱矩陣形的次之，光滑壁面的最少。②對于光滑壁面上的納米氬液膜蒸發(fā)過程，隨著模擬溫度的提高，恒速蒸發(fā)階段的時間變短，蒸發(fā)通量迅速增大；在恒速蒸發(fā)階段，能量系數(shù)對蒸發(fā)通量幾乎沒有影響；在降速蒸發(fā)階段，隨著能量系數(shù)的增加，固體壁面吸附的氬流體分子數(shù)增大。③對于納米欄柵形粗糙壁面上的納米氬液膜蒸發(fā)，隨著粗糙度因子的增大，恒速蒸發(fā)階段的時間增加，蒸發(fā)通量減?。划斦舭l(fā)過程達到平衡后，固體壁面吸附的氬原子數(shù)增多。④對于納米欄柵形粗糙壁面上的納米氬液膜蒸發(fā)，隨著相面積分數(shù)的增大，恒速蒸發(fā)階段的時間增加，蒸發(fā)通量減?。划斦舭l(fā)過程達到平衡后，固體壁面吸附的氬原子數(shù)增多。

[1] Yu J,Wang H.Thin-film phase change driven by disjoining pressure difference[J].Heat and Mass Transfer,2012,48(7):1135-1140.

[2] Hens A,Agarwal R,Biswas G.Nanoscale study of boiling and evaporation in a liquid Ar film on a Pt heater using molecular dynamics simulation[J].International Journal of Heat and Mass Transfer,2014,71:303-312.

[3] Shavik S M,Hasan M N,Morshed A K M M,et al.Molecular dynamics study of effect of different wetting conditions on evaporation and rapid boiling of ultra-thin argon layer over platinum surface[J].Procedia Engineering,2015,105:446-451.

[4] Wang W D,Zhang H Y,Tian C H,et al.Numerical experiments on evaporation and explosive boiling of ultra-thin liquid argon film on aluminum nanostructure substrate[J].Nanoscale Research Letters,2015,10(1):1-14.

[5] Heermann D W.Computer-simulation methods[M].Springer Berlin Heidelberg,1990.

[6] Long L N,Micci M M,Wong B C.Molecular dynamics simulations of evaporation[J].Computer Physics Communications,1996,96:167-172.

[7] 李瑞靜.納米液膜傳遞特性的分子動力學(xué)模擬[D].大連:大連理工大學(xué),2016.

[8] 蔣大林.固液界面潤濕的分子動力學(xué)研究及實驗[D].蘇州:江蘇大學(xué),2007.

[9] Sumardiono S,Fischer J.Molecular simulations of droplet evaporation by heat transfer[J].Microfluidics and Nanofluidics,2007,3(2):127-140.

Study on Molecular Dynamics Simulation of Evaporation Behavior of Argon Nanofilm

WANG Bao-he1, LI Rui-jing1, WANG Wei2

(1.School of Chemical Engineering , Dalian University of Technology , Dalian 116024 , China ; 2.School of Chemical Machinery and Safety , Dalian University of Technology , Dalian 116024 , China)

Non-equilibrium molecular dynamics simulation technology is used to research influence of factors such as simulated temperature,energy factor,roughness factor,phase area fraction on evaporation behavior of argon nanofilm on solid surfaces.Simulation results show that there are two stages of a constant rate evaporating period and a falling rate evaporating period for evaporation process of argon nanofilm on smooth or nanotextured rough surfaces;evaporation flux has little change during the constant rate evaporating period or the falling rate evaporating period,but there is difference in number of argon molecules absorbed on the solid surfaces and the ranking (large to small) is as follows:fencing form surface,square column matrix form surface,smooth surface after evaporation process reaches equilibrium state.For evaporation of argon nanofilm on the smooth surface,the constant rate evaporating period becomes short and evaporation flux increases quickly with the increasing of simulated temperature,and with the increase of energy factor,evaporation flux changes hardly at the constant rate evaporating stage;however,number of argon molecules absorbed on the solid surface increases at the falling rate evaporating stage.For evaporation argon nanofilm on the fencing form surface,with increase of roughness factor or phase area fraction,the constant rate evaporating period becomes long and evaporation flux decreases;but number of argon molecules absorbed on the solid surface increases after evaporation process reaches equilibrium state.

argon nanofilm ; molecular dynamics ; simulation ; evaporation

2016-10-28

國家自然科學(xué)基金(21676042)

王寶和(1959-)，男，副教授，從事不同形貌微納結(jié)構(gòu)的制備、干燥及分子動力學(xué)模擬研究，電話：0411-84986167，E-mail：wbaohe@163.com。

TB383；O561

A

1003-3467(2017)01-0016-05