水滴撞擊結(jié)冰過程的分子動力學(xué)模擬?

董琪琪 胡海豹2) 陳少強 何強 鮑路瑤

1)(西北工業(yè)大學(xué)航海學(xué)院,西安 710072)

2)(西北工業(yè)大學(xué)深圳研究院,深圳 518057)

3)(中國船舶重工集團(tuán)公司第七〇五研究所,西安 710077)

(2017年10月6日收到;2017年12月11日收到修改稿)

1 引 言

水滴結(jié)冰是常見的自然現(xiàn)象,也給人類生活帶來無數(shù)的挑戰(zhàn).過度積冰會對基礎(chǔ)設(shè)施造成嚴(yán)重破壞,如建筑物、輸電線路、風(fēng)力渦輪機(jī)、太陽能電池板、船舶和氣動機(jī)翼結(jié)構(gòu)等[1?9].這些結(jié)冰現(xiàn)象實際為水滴撞擊到冷表面的動態(tài)結(jié)冰問題.因此,為真實反映壁面水滴結(jié)冰的規(guī)律,亟需深入研究水滴撞擊結(jié)冰過程.

在實驗研究方面,Quero等[10]在研究過冷水滴撞擊水膜表面的凍結(jié)過程時發(fā)現(xiàn),與靜態(tài)結(jié)冰相似,水滴撞擊冷表面的凍結(jié)過程受表面溫度和環(huán)境溫度的影響最大.Mishchenko等[11]通過水滴撞擊超疏水表面結(jié)冰實驗發(fā)現(xiàn),表面溫度在?25—?30°C之間的超疏水表面可使水滴在形核之前回彈,能抑制結(jié)冰.Jung等[12]發(fā)現(xiàn)材料的抑冰性能同時受接觸角和表面粗糙度的影響,同等粗糙度情況下,超疏水材料比親水材料具有更強的抑冰性能.Li等[13]研究了凍結(jié)對鋁板表面水滴撞擊物理過程的影響,結(jié)果表明,凍結(jié)不影響水滴的擴(kuò)散過程,但可以減小回縮過程.Yang等[14]對過冷水滴撞擊不同材質(zhì)金屬管表面的凍結(jié)機(jī)理進(jìn)行了研究,指出除環(huán)境溫度和冷表面溫度外,水滴過冷度、表面特性對結(jié)冰也有影響,且表面溫度越低撞擊速度對結(jié)冰的影響越弱.

在數(shù)值模擬方面,張大林等[15]模擬發(fā)現(xiàn),隨著過冷水滴平均直徑的增大,結(jié)冰區(qū)的面積和水收集系數(shù)越大.楊倩等[16]數(shù)值模擬給出了飛行高度、速度及水滴半徑對撞擊特性的影響規(guī)律.陳科和曹義華[17]則發(fā)現(xiàn)飛機(jī)發(fā)生結(jié)冰時,冷空氣會先在飛機(jī)表面凝華成霜,之后再誘導(dǎo)結(jié)冰.盛強等[18]對機(jī)翼流場進(jìn)行數(shù)值模擬,在計算機(jī)翼表面的結(jié)冰量的同時,分析了結(jié)冰對機(jī)翼阻力、升力和壓力系數(shù)的影響.

需要指出的是,現(xiàn)有研究報道主要關(guān)注了宏觀尺度下水滴撞擊冷表面的結(jié)冰過程,對微觀尺度下冷壁面上水滴撞擊結(jié)冰過程仍缺乏精細(xì)的探究手段.為此,本文采用三維分子動力學(xué)模擬方法,從納米尺度研究水滴撞擊冷壁面的結(jié)冰過程,并在分析水滴結(jié)冰判斷方法的基礎(chǔ)上,初步探究了溫度和表面能對水滴結(jié)冰的影響規(guī)律.

2 分子動力學(xué)模擬方法

2.1 模擬系統(tǒng)

建立了一個包含納米級水滴和固體壁面的三維分子動力學(xué)模型,如圖1所示.其中,固體壁面長度L=262.0875 ?,壁厚D=30.7 ?,z方向長度為289.2 ?.

壁面采用fcc晶格結(jié)構(gòu),晶格常數(shù)a=3.615 ?,整個壁面分為9層,共有196667個固體原子組成.其中,位于壁面最下端的3層原子為固定原子,作用是阻止模型中的壁面原子和水滴原子下移;中間4層原子為溫度控制原子組,用于調(diào)控整個模型,使其保持在一個穩(wěn)定的溫度;最上端的2層原子為自由原子組,該組原子的厚度為7.23 ?.

直徑為幾個納米的水滴就可以準(zhǔn)確反映水滴特性[19?21],因此,為減少計算量,選取水滴直徑為9 nm.水分子共15625個,由15625個氧原子和31250個氫原子組成.模型中的氣相成分則是由水分子在相應(yīng)模擬條件下自由擴(kuò)散形成,且所有氣相的密度、原子質(zhì)量以及勢能參數(shù)均接近真實的氣相環(huán)境.

圖1 模擬系統(tǒng)示意圖Fig.1.Schematic of simulation system.

2.2 數(shù)學(xué)模型

這里使用TIP4P/ice模型來模擬氧原子、氫原子和水分子,具體勢能函數(shù)形式如下:

其中roo為兩氧原子之間的距離,ε和σ是LJ勢能的特征能量和特征長度,e為質(zhì)子電荷,ε0為介電常數(shù),a和b為分子的帶電位置.水分子的鍵長和鍵角則使用SHAKE算法來進(jìn)行固定.

溫度低于熔點時,水滴處于固態(tài),水分子會表現(xiàn)出近程和長程均有序的特性;溫度高于熔點時,水滴處于液態(tài),水分子會表現(xiàn)出近程有序、長程無序的特性.利用TIP4P/ice模型,能夠自動捕捉到水滴結(jié)冰和融化過程.在其他特性不受影響的情況下,該模型預(yù)測的六邊形冰的最優(yōu)熔點為?0.95°C/1 bar(1 bar=0.1 MPa)[22].因此,該模型非常適合研究水滴撞擊結(jié)冰.值得注意的是水滴的表面張力(大量水分子間實際存在的相互作用的統(tǒng)計結(jié)果)已經(jīng)包含在了TIP4P/ice模型中.

另外參考文獻(xiàn)[23],采用LJ/126模型模擬水分子和固體原子之間的相互作用,通過調(diào)控勢能參數(shù)來模擬不同潤濕性壁面.

2.3 模擬算法

模擬中,統(tǒng)計系統(tǒng)均采用微正則系綜(NVE),溫度控制策略采用速度定標(biāo)法.而在速度定標(biāo)法中,采用場基溫控法(profile-unbiased thermostat)[24]來計算溫度.系統(tǒng)中所有水滴分子運動均滿足牛頓運動方程:

這里采用文萊特算法來迭代求解這些運動方程,

模擬時間步設(shè)定為0.002 ps,運行總步數(shù)為400000步.

2.4 結(jié)冰判斷方法

相比于固體分子,液體分子之間空隙較大,分子不會長時間地在同一個固定位置上振動,而是振動一段時間后,轉(zhuǎn)到另一個平衡位置上振動,即液體分子可以在液體內(nèi)移動.因此,若水分子中氧原子位置保持不變,則可判定其所在位置已從液態(tài)轉(zhuǎn)變?yōu)楣虘B(tài).圖2為采用上述模型模擬出的?40°C冷壁面上10260號氧原子的運動軌跡,可以看出,約75 ps之后氧原子三個方向的坐標(biāo)基本保持不變(僅在該位置振動),即此時該氧原子所屬水分子被凍結(jié).表明通過此方法理論上可以分析出每個水分子的結(jié)冰時間.

圖2 10260號氧原子的軌跡變化坐標(biāo)圖Fig.2.The coordinate graph of trajectory variation of oxygen atom No.10260.

但由于水滴中含水分子數(shù)眾多,對不同時刻下各分子位置的統(tǒng)計過于復(fù)雜,致使計算量在現(xiàn)有硬件條件下無法承受.為此,這里利用TIP4P/ice模型中水的凝固溫度恒為?0.95°C的特點,提出了一種通過統(tǒng)計不同時刻水滴內(nèi)水分子的溫度來判斷水滴結(jié)冰的方法.考慮到水滴結(jié)冰自下而上的規(guī)律,這里采用垂直分層統(tǒng)計溫度的方法來簡化計算量.當(dāng)水滴頂部溫度也達(dá)到結(jié)冰點(?0.95°C)時,判定水滴完全結(jié)冰.

3 模擬結(jié)果與分析

3.1 溫度對水滴結(jié)冰的影響

模擬4種壁面溫度下的水滴撞擊結(jié)冰過程.模擬參數(shù):固-液/液-液勢能比ε=11.5,水滴起始溫度20°C,撞擊速度為300 m/s,壁面溫度與空間溫度(環(huán)境溫度)相同,分別為20,?40,?50,?100°C.

水滴撞擊結(jié)冰過程如圖3所示,可以看出,冷壁面上水滴最大鋪展直徑明顯小于常溫,且水滴最大鋪展直徑隨著壁面溫度的降低而減小.同時,圖3中水滴撞擊20°C壁面時,壁面上固/液接觸線在740 ps趨于穩(wěn)定;而撞擊?40,?50,?100°C壁面時,接觸線分別在95,81和60 ps后基本保持不變,且壁面溫度越低,水滴接觸線振蕩時間越短.胡海豹等[25]的實驗結(jié)果表明,同種材料的壁面上,溫度T越低,最大鋪展系數(shù)βmax越小,而且,低溫壁面上固/液接觸線在振蕩過程中始終保持不變.這與圖3中的仿真結(jié)果相一致,證明此模型是可行的.

圖3 光滑壁面上水滴撞擊結(jié)冰時序圖Fig.3.Sequence diagram of freezing behavior of water droplet on the smooth wall.

水滴撞擊?50°C和?40°C壁面時,垂直方向各層水分子的溫度變化情況如圖4所示.分析?40°C壁面水滴撞擊結(jié)冰過程(圖4(a))可知,120 ps時水滴底部(靠近固體層)的水分子溫度已低于結(jié)冰點,而水滴頂部的水分子溫度仍高于結(jié)冰點,這說明在水滴由底部到頂端的結(jié)冰過程中存在著較大的溫度梯度136 ps時水滴頂部的溫度達(dá)到結(jié)冰點,水滴完全結(jié)冰;330 ps時水滴完全降至壁面溫度.而水滴撞擊?50°C壁面結(jié)冰過程中(圖4(b)),132 ps時水滴頂部的溫度就已達(dá)到結(jié)冰點,471 ps后水滴才完全降至壁面溫度.說明壁面溫度越低,水滴完全結(jié)冰所需的時間越少,但水滴降至壁面溫度所需的時間卻越多.這是由于隨著壁面溫度的降低,固/液間熱流密度增大,單位時間內(nèi)水滴被吸收的熱量增多,從而使水滴整體溫度快速降低,減少了完全結(jié)冰的時間;但水滴需要降低較多的溫度和進(jìn)行較多的熱量傳遞才能降至較低的壁面溫度,并且越接近壁面溫度時,熱量傳遞的速率就會越慢,所以壁面溫度越低,水滴降至壁面溫度所需的時間越多.

圖4 不同溫度壁面水滴垂直方向各層水分子的溫度分布(a)T=?40°C;(b)T=?50°CFig.4. Thetemperaturedistributionofwater molecules in the vertical direction of different temperature wall:(a)T= ?40°C;(b)T=?50°C.

3.2 壁面表面能屬性對水滴結(jié)冰的影響

固/液間作用勢能可以影響壁面的表面能屬性,作用勢越大,壁面越親水,反之越疏水.本文模擬了不同表面能壁面上水滴撞擊結(jié)冰過程,其中,水滴起始溫度為20°C,撞擊速度為300 m/s,壁面溫度和環(huán)境溫度為?40°C,固-液/液-液作用勢能比(ε)依次為8,6,2(壁面親水性越來越弱).

圖5 不同表面能壁面水滴垂直方向各層水分子的溫度分布 (a) ε =8;(b) ε =6;(c) ε =2Fig.5. Thetemperaturedistributionofwater molecules in the vertical direction of different surface energy wall:(a) ε =8;(b) ε =6;(c) ε =2.

圖5為不同表面能壁面水滴垂直方向各層水分子的溫度分布.可以明顯看出,同一時刻下,ε=2時的溫度分布明顯高于ε=8,6時的溫度分布;且ε=2時液滴各分子層間溫差較小,內(nèi)部溫度相對均勻.ε=8,6,2時,水滴完全結(jié)冰時間分別為133,136,198 ps,水滴完全降至壁面溫度所需時間分別為462,560,755 ps;即隨著親水性的減弱,水滴完全結(jié)冰和降至壁面溫度的時間都逐漸增多.這是由于隨著壁面親水性降低,水滴內(nèi)部的熱量傳遞的速率會減慢(特別是冷壁面與近壁面的水分子層間的熱傳遞),延長了結(jié)冰時間.

4 結(jié) 論

采用分子動力學(xué)方法模擬研究了水滴撞擊冷壁面的結(jié)冰過程,結(jié)果發(fā)現(xiàn):

1)通過統(tǒng)計水滴垂直方向水分子的溫度判定水滴是否結(jié)冰比通過統(tǒng)計微觀原子的位置坐標(biāo)來判定更為簡潔;

2)壁面溫度是水滴結(jié)冰現(xiàn)象的重要影響因素之一,壁面溫度越低,水滴完全結(jié)冰所需的時間越少,但降至壁面溫度所需的時間卻越大;

3)壁面親水性是水滴結(jié)冰現(xiàn)象的另一重要影響因素,隨著壁面親水性降低,水滴內(nèi)部的熱傳遞速度減慢(尤其是冷壁面與水滴底端分子層間),內(nèi)部溫度趨于均勻,延長了水滴結(jié)冰時間.

另外,水滴晶格與表面晶格的匹配度同樣會影響水滴的結(jié)冰現(xiàn)象[26,27].因此,后續(xù)仍需要進(jìn)一步研究晶格結(jié)構(gòu)對結(jié)冰特性的影響,揭示其對結(jié)冰過程的影響機(jī)制.

[1]Jung S,Tiwari M K,Doan N V,Poulikakos D 2012Nat.Commun.3 615

[2]Jin Z,Wang Z,Sui D 2016Int.J.Heat.Mass.Trans.97 211

[3]Wang Y,Orol D,Owens J,Simpson K,Lee H J 2013Mater.Sci.Appl.04 347

[4]Dalili N,Edrisy A,Carriveau R 2009Renew Sust.Energ.Rev.13 428

[5]Zou L,Xu H J,Gong S K,Li D W 2010China Safety Sci.J.20 105(in Chinese)[周莉,徐浩軍,龔勝科,李大偉2010中國安全科學(xué)學(xué)報20 105]

[6]Xiao S,He J,Zhang Z 2017Acta Mech.Solida Sin.30 224

[7]Yao Y,Li C,Zhang H,Yang R 2017Appl.Surf.Sci.419 52

[8]Zou M,Beckford S,Wei R,Ellis C,Hattonc G,Millerb M A 2011Appl.Surf.Sci.257 3786

[9]Zhang C,Liu H 2016Phys.Fluids28 260

[10]Quero M,Hammond D W,Purvis R,Smith F T 2006AIAA466

[11]Mishchenko L,Hatton B,Bahadur V,Taylor J A,Krupenkin T,Aizenberg J 2010ACS Nano4 7699

[12]Jung S,Dorrestijn M,Raps D,Das A,Megaridis C M,Poulikakos M 2011Langmuir27 3059

[13]Li H,Roisman I V,Tropea C 2011Proceeding of the Sixth International Conference on Fluid Mechanics1376 451

[14]Yang G,Guo K,Li N 2011Int.J.Refrig.34 2007

[15]Zhang D L,Yang X,Ang H S 2003J.Propul.Power18 87(in Chinese)[張大林,楊曦,昂海松2003航空動力學(xué)報18 87]

[16]Yang Q,Chang S N,Yuan X G 2002Acta Aeronaut.Astronaut.Sin.23 173(in Chinese)[楊倩,常士楠,袁修干2002航空學(xué)報23 173]

[17]Chen K,Cao Y H 2008Aeronaut.Comput.Tech.38 36(in Chinese)[陳科,曹義華2008航空計算技術(shù)38 36]

[18]Sheng Q,Xing Y M,He C 2009Aeronaut.Comput.Tech.39 37(in Chinese)[盛強,邢玉明,何超 2009航空計算技術(shù)39 37]

[19]Yuan Q Z,Zhao Y P 2010Phys.Rev.Lett.104 246101

[20]Xiao S,He J Y,Zhang Z X 2016Nanoscale8 14625

[21]Bi Y,Cao B,Li T 2017Nat.Commun.8 15372

[22]Abascal J L,Sanz E,García F R,Vega C 2005J.Chem.Phys.122 234511

[23]Hong S D,Ha M Y,Balachandar S 2009J.Colloid Interf.Sci.339 187

[24]Evans D J,Morriss G P 1986Phys.Rev.Lett.56 2172

[25]Hu H B,He Q,Yu S X,Zhang Z Z,Song D 2016Acta Phys.Sin.65 104703(in Chinese)[胡海豹,何強,余思瀟,張招柱,宋東2016物理學(xué)報65 104703]

[26]Fitzner M,Sosso G C,Cox S J,Michaelides A 2015J.Am.Chem.Soc.137 13658

[27]Liu K,Wang C,Ma J,Shi G,Yao X,Fang H,Song Y L,Wang J J 2016Proc.Natl.Acad.Sci.USA113 14739