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

董琪琪 胡海豹2) 陳少?gòu)?qiáng) 何強(qiáng) 鮑路瑤

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

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

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

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

1 引 言

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

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

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

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

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

2.1 模擬系統(tǒng)

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

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

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

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

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

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

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

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

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

2.3 模擬算法

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

這里采用文萊特算法來(lái)迭代求解這些運(yùn)動(dòng)方程,

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

2.4 結(jié)冰判斷方法

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

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

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

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

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

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

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

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

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

圖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 壁面表面能屬性對(duì)水滴結(jié)冰的影響

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

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

4 結(jié) 論

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

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

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

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

另外,水滴晶格與表面晶格的匹配度同樣會(huì)影響水滴的結(jié)冰現(xiàn)象[26,27].因此,后續(xù)仍需要進(jìn)一步研究晶格結(jié)構(gòu)對(duì)結(jié)冰特性的影響,揭示其對(duì)結(jié)冰過(guò)程的影響機(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中國(guó)安全科學(xué)學(xué)報(bào)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航空動(dòng)力學(xué)報(bào)18 87]

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

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

[18]Sheng Q,Xing Y M,He C 2009Aeronaut.Comput.Tech.39 37(in Chinese)[盛強(qiáng),邢玉明,何超 2009航空計(jì)算技術(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)[胡海豹,何強(qiáng),余思瀟,張招柱,宋東2016物理學(xué)報(bào)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