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

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

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

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

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

1 引 言

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

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

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

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

2 分子動力學模擬方法

2.1 模擬系統(tǒng)

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

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

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

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

2.2 數(shù)學模型

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

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

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

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

2.3 模擬算法

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

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

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

2.4 結冰判斷方法

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

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

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

3 模擬結果與分析

3.1 溫度對水滴結冰的影響

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

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

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

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

圖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 壁面表面能屬性對水滴結冰的影響

固/液間作用勢能可以影響壁面的表面能屬性,作用勢越大,壁面越親水,反之越疏水.本文模擬了不同表面能壁面上水滴撞擊結冰過程,其中,水滴起始溫度為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時,水滴完全結冰時間分別為133,136,198 ps,水滴完全降至壁面溫度所需時間分別為462,560,755 ps;即隨著親水性的減弱,水滴完全結冰和降至壁面溫度的時間都逐漸增多.這是由于隨著壁面親水性降低,水滴內(nèi)部的熱量傳遞的速率會減慢(特別是冷壁面與近壁面的水分子層間的熱傳遞),延長了結冰時間.

4 結 論

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

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

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

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

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

[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中國安全科學學報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航空動力學報18 87]

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

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

[18]Sheng Q,Xing Y M,He C 2009Aeronaut.Comput.Tech.39 37(in Chinese)[盛強,邢玉明,何超 2009航空計算技術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物理學報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