不同浸潤性冷表面上水滴碰撞結(jié)冰的數(shù)值模擬

冷夢堯，常士楠，丁　亮

?

研究簡報

不同浸潤性冷表面上水滴碰撞結(jié)冰的數(shù)值模擬

冷夢堯，常士楠，丁亮

（ 北京航空航天大學(xué)航空科學(xué)與工程學(xué)院，北京 100191）

摘要:對冷水滴撞擊不同表面時的動力學(xué)行為和相變過程進(jìn)行了模擬。通過耦合 VOF和Level-set方法追蹤氣液自由界面，結(jié)合焓-孔隙度相變模型，模擬水滴撞擊冷表面的動力學(xué)行為及相變特征。選取親水(接觸角30°)、疏水(接觸角114°)和超疏水(接觸角163°) 3種典型浸潤性的表面，計算了多種壁溫條件下的水滴撞擊結(jié)冰過程。結(jié)果表明提高表面疏水性，將減小水滴與冷表面的接觸時間和接觸面積，降低水滴內(nèi)的相變速率，延緩水滴結(jié)冰的時間。在表面溫度高于-15℃時，超疏水表面可以避免冷水滴的凍結(jié)黏附，保持表面潔凈。將模擬得到的最大鋪展直徑、回縮速率以及凍結(jié)情況，與已有實驗結(jié)果進(jìn)行對比驗證，表明了模擬方法的有效性和準(zhǔn)確性。

關(guān)鍵詞:表面；液滴碰撞; 動力學(xué); 多相流；相變；數(shù)值模擬

引 言

水滴撞擊冷表面發(fā)生結(jié)冰的現(xiàn)象，屬于氣-液-固三相流的流動傳熱問題，包含水滴撞擊后的動力學(xué)行為以及水的凝固相變過程；而這二者之間的相互影響，使得水滴碰撞表面結(jié)冰往往會導(dǎo)致嚴(yán)重的后果且難以準(zhǔn)確預(yù)測[1-2]。近幾年研究者們通過實驗發(fā)現(xiàn)，提高表面疏水性能夠獲得減少結(jié)冰量[3]，降低防冰能耗[4]，削弱冰的黏附力[5]等效果。然而，傳統(tǒng)的結(jié)冰理論模型[6-7]主要從水滴撞擊時的各項熱流出發(fā)，通過熱力學(xué)平衡判斷水滴發(fā)生凍結(jié)與否，無法解釋表面特性對水滴碰撞后的相變行為的影響。其原因在于計算時假設(shè)相變過程的時間極短而將其忽略，缺少對水滴撞擊表面時其動力學(xué)行為與相變過程相互影響的考慮。

關(guān)于過冷水滴靜止?fàn)顟B(tài)下的結(jié)冰過程已有相關(guān)文獻(xiàn)研究[8-11]。實驗表明，水滴的完全凍結(jié)過程平均長達(dá)幾十秒[12]；在超疏水表面上，水滴凍結(jié)的時間會進(jìn)一步被延長[13]，而碰撞過程的時間則要短得多。對于冷水滴碰撞的動態(tài)過程，Mishchenko等[14]和 Li等[15]的研究均揭示出超疏水表面可以顯著改變結(jié)冰情況，有效防止水滴的凍結(jié)黏附。但是受實驗手段的限制，運動水滴內(nèi)部的相變詳情難以體現(xiàn)。作為實驗觀測的補充手段，研究者們通過數(shù)值模擬研究水滴動力學(xué)行為，如VOF和Level-set相結(jié)合的方法，獲得了與實驗結(jié)果符合良好的結(jié)果[16-17]。數(shù)值模擬能夠獲取更為全面的水滴內(nèi)部情況，通過分析水滴內(nèi)的壓力場[18]、速度矢量[19-20]和溫度分布[21-22]等，能夠進(jìn)一步解釋水滴動力學(xué)行為和傳熱機理。

本工作的目標(biāo)是建立能夠準(zhǔn)確地同時追蹤氣液界面和揭示液固相變的數(shù)值方法，獲得表面浸潤性對水滴碰撞凍結(jié)的影響。針對水滴碰撞不同浸潤性的冷表面，將多相流的界面追蹤方法（耦合Level-set和VOF方法）和相變模型（焓-孔隙度模型）相結(jié)合，模擬了冷水滴撞擊到不同表面上的鋪展-回縮行為，與實驗得到的濕潤因子對比驗證本方法的準(zhǔn)確性，并進(jìn)一步揭示水滴內(nèi)部相變區(qū)域的發(fā)展過程以及相變與回縮行為的相互作用，研究了表面特性對水滴內(nèi)部相變動態(tài)過程的影響，以期為表面結(jié)冰的界面現(xiàn)象和防冰/除冰機理研究提供參考。

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

1.1 自由界面追蹤模型

本文采用VOF耦合Level-set方法[23]追蹤氣液自由界面，流場求解以 VOF的思想為主，在考慮界面效應(yīng)時，以Level-set方法為輔助。VOF方法能夠保持追蹤流體的“質(zhì)量”，并且能夠容易地追蹤拓?fù)浣Y(jié)構(gòu)變化的流體界面。在 VOF方法中，相與相之間的界面追蹤是通過計算次相的體積分?jǐn)?shù)連續(xù)方程完成的。由于本文的研究對象為兩相流，主相為空氣，次相為水，故而只需要求解次相體積分?jǐn)?shù)的連續(xù)方程

主相的體積分?jǐn)?shù)則通過總體積分?jǐn)?shù)之和為1計算得到，αP=1-αS。動量方程和能量方程中的物性參數(shù)（如密度、黏度等）則通過各相體積分?jǐn)?shù)進(jìn)行加權(quán)計算。

VOF方法通過幾何重構(gòu)法對界面進(jìn)行插值構(gòu)造，即基于部分填充單元內(nèi)的體積分?jǐn)?shù)及其梯度，用分段線性函數(shù)構(gòu)造相界面；然后根據(jù)通過界面的流量，更新單元內(nèi)體積分?jǐn)?shù)。但是分段線性界面不符合光滑和連續(xù)的要求，無法提供準(zhǔn)確的界面方向和曲率，故而在此引入 Level-set方法來解決這個問題。

Level-set方法是將氣液相界面用一個高階函數(shù)，即Level-set函數(shù)?的零值點來表示，由?的代數(shù)值來區(qū)分計算區(qū)域中的各相。Level-set函數(shù)?為一距離函數(shù)，其定義如下

但是由于數(shù)值計算過程的內(nèi)在效應(yīng)，Level-set函數(shù)在經(jīng)過幾個時層的計算后，將不再是距離函數(shù)，需要進(jìn)行距離函數(shù)的初始化運算，使得Level-set函數(shù)在不改變其0等值面位置的前提下重新構(gòu)成為距離函數(shù)。另外，Level-set方法的數(shù)值求解過程中，由于數(shù)值耗散和數(shù)值誤差等原因可能造成封閉界面內(nèi)流體質(zhì)量不守恒[24]。

耦合VOF和Level-set的自由界面追蹤方法將這兩者各自的優(yōu)勢結(jié)合起來:首先，在構(gòu)造相界面時，通過相體積分?jǐn)?shù)的值確定單元內(nèi)分界面位置，依靠Level-set函數(shù)?的梯度確定界面法向方向；然后采用函數(shù)?計算相界面的曲率及其相關(guān)量，提高相界面參數(shù)的求解精度；該時間步結(jié)束時，通過相體積分?jǐn)?shù)確定的界面位置更新流場中的Level-set函數(shù)，以克服Level-set方法求解對流輸運方程時的質(zhì)量不守恒問題。

1.2 流動計算

表面張力、壁面浸潤性以及相變過程對兩相流流動的影響，通過對基本動量方程添加源項來實現(xiàn)。式(4)為流場內(nèi)動量方程的表達(dá)式

等式左邊依次為非穩(wěn)態(tài)項以及對流項，右邊前3項分別為壓力項、黏性項和體積力項。等式右邊的-σκδ（?）??為基于連續(xù)表面力(CSF)模型添加的表面張力源項，以考慮表面張力σ和自由界面的曲率κ對流動的影響。SF為凝固引起的流動阻力。

壁面處控制體內(nèi)的自由界面法向可通過定義壁面接觸角來確定，即通過修正壁面附近的界面法向向量n，令其等于式(5)從而改變界面曲率及表面張力源項。式(5)中nw和tw分別表示壁面的法向和切向單位向量，而θw為定義的壁面接觸角

1.3 凝固-融化模型

實驗表明，水在過冷狀態(tài)下結(jié)冰與水從常溫下逐漸冷卻結(jié)冰的過程不同，其包含兩個過程[25]:①快速結(jié)冰過程，當(dāng)受到某種觸發(fā)作用（如降溫、振動、碰撞等）導(dǎo)致冰核成形，接著冰枝迅速生長并充滿水滴內(nèi)部，水滴相變部分變成疏松的海綿狀結(jié)構(gòu)，并升溫至冰-水混合平衡溫度；②緩慢凍結(jié)過程，若冷表面或冷空氣繼續(xù)從水滴中抽走熱量，則水滴會逐漸提高固體的比例，這一過程是等溫變化，過程的長短由傳熱速率決定。

根據(jù)實驗觀測的結(jié)果，采用一種焓-孔隙度相變模型來模擬水滴內(nèi)部的凝固-融化相變過程[26]，計算求解的能量方程為

考慮相變后，能量方程中的焓值為顯熱及潛熱之和，H=cpT+ΔH。其中潛熱部分與水的凝固/融化相變潛熱的關(guān)系為

焓-孔隙度相變模型將冰水混合區(qū)域視為多孔結(jié)構(gòu)，孔隙率即為液相的體積分?jǐn)?shù)γ，取值區(qū)間為[0,1]，其與溫度的關(guān)系為

假設(shè)冰-水混合區(qū)內(nèi)的運動遵守多孔介質(zhì)內(nèi)流動的Darcy定律，且符合Carman-Koseny假設(shè)，則可推出動量方程源項為[26]

其中，ε為極小值（ε=0.001），以避免分母為零；Cmush為黏糊系數(shù)，與多孔介質(zhì)的形態(tài)相關(guān)，合理的Cmush常數(shù)取值，既可以在液相比例較高時保持流動性，又可以在液相比例降低后抑制流動速率。

從動量源項的表達(dá)式可以看出，在液相區(qū)，即γ=1，動量源項為0，動量方程與流體N-S方程形式一致；凝固開始變成冰-水混合物后，液相體積分?jǐn)?shù)γ降低，源項在動量方程中的比重逐漸超過瞬態(tài)項、對流項和擴(kuò)散項，直至完全轉(zhuǎn)變?yōu)楣滔嗪?，速度基本為零?/p>

圖1　計算區(qū)域Fig.1 Droplet impingement case setup

2　計算模型

通過模擬冷水滴撞擊不同表面上的動態(tài)相變過程，驗證本方法的可行性。本文的模擬區(qū)域為高4倍水滴直徑、寬3倍水滴直徑的二維軸旋轉(zhuǎn)對稱區(qū)域，計算模型及邊界條件如圖1所示，模型底部為定溫?zé)o滑移壁面邊界，頂部及側(cè)邊界均為大氣邊界條件。對網(wǎng)格進(jìn)行無關(guān)性檢驗后，決定以沿初始直徑分布200個網(wǎng)格為準(zhǔn)進(jìn)行劃分，最終計算區(qū)域由480000個正方形結(jié)構(gòu)化網(wǎng)格組成，能夠準(zhǔn)確地捕捉到氣液自由界面以及相界面的移動變化。計算中壓力速度耦合采用適合非定常流動問題的PISO方法，氣液自由界面重構(gòu)使用Geo-Reconstruct方法。除梯度和壓力分別使用Least Squares Cell Based和PRESTO！格式進(jìn)行離散外，剩余各項均使用二階迎風(fēng)格式。為節(jié)省計算資源并保證迭代結(jié)果的準(zhǔn)確性，計算過程采用變時間步長，控制Courant數(shù)始終小于0.25，時間離散格式為一階隱式格式，每個時間步內(nèi)迭代20次。

本文使用FLUENT 14.5軟件進(jìn)行模擬，計算工質(zhì)為空氣和水，空氣相為主相，水為第2相。為了同實驗結(jié)果進(jìn)行對比，選擇直徑3 mm的水滴，以初速度1.4 m·s-1自1倍直徑高度下落撞擊固體壁面，初始時刻，空氣溫度為-5℃，水滴溫度為相平衡溫度。計算了3種典型浸潤性表面，其接觸角分別為30°、114°和163°，壁面溫度為20、-15、-25℃共9種碰撞情形。計算過程中使用的水滴物性參數(shù)如表1所示。

表1　水滴物性參數(shù)Table 1 Material properties for water

3　結(jié)果與討論

3.1 表面溫度對水滴碰撞動力學(xué)行為的影響

圖2給出了3種典型浸潤性表面，在不同溫度條件下水滴碰撞后的形態(tài)演化過程。在Blake等[16]對純水結(jié)冰所做相關(guān)研究的基礎(chǔ)上，本文計算中使用的黏糊系數(shù) Cmush如表 2所示，從中可以看出當(dāng)表面溫度越低或浸潤性越強，即條件越有利于加快凝固時，則Cmush越大，這與Cmush象征的物理意義相符合。

圖2　不同溫度及不同浸潤性表面上液滴碰撞后的形態(tài)變化過程Fig.2 Dynamic behavior of droplet impinging on surfaces of different wettability and temperature

首先對比常溫（20℃）條件無相變發(fā)生時，不同浸潤性表面上水滴的鋪展和回縮過程，結(jié)果表現(xiàn)出不同浸潤性表面上液滴碰撞的典型動態(tài)特征。在親水性表面（θ=30°）上，水滴可達(dá)到最大的鋪展面積，且回縮過程緩慢，最后水滴收縮至圓盤狀并停留在表面上。增大表面接觸角后，從水滴的邊緣形狀可以看出，水滴的鋪展?jié)櫇裥袨槭茏?，不僅最大鋪展直徑減小，達(dá)到最大鋪展的時間，即回縮過程開始的時間，也相應(yīng)提前。在疏水（θ=114°）及超疏水（θ=163°）表面上，水滴同表面之間摩擦損耗以及水滴內(nèi)部的黏性耗散減少，因此，水滴剩余的動能較大，其回縮過程加快，最終形成完全反彈離開表面。

表3　黏糊系數(shù)Cmush的取值Table 2 Value of mushy zone constant Cmush

對比相同浸潤性、不同溫度條件下水滴的動力學(xué)行為發(fā)現(xiàn)，水滴前半段的鋪展過程都十分相似，且在相同浸潤性的表面上，水滴最大鋪展形狀也十分相似。這說明水滴的撞擊熱以及液體流動時產(chǎn)生的黏性摩擦熱能夠一定程度上延緩相變的發(fā)生。由于在水滴達(dá)到最大鋪展直徑前后，也即是水滴動能最小的階段，水滴與表面之間換熱面積最大，且黏性及摩擦產(chǎn)生的加熱熱流最低，故而形成了最有利于結(jié)冰發(fā)生的條件。從水滴形態(tài)對比圖（圖 2）中可以看出，在常溫表面和遠(yuǎn)低于冰點的表面上，水滴碰撞后半段的回縮過程有明顯不同的動態(tài)行為。在低于冰點的表面上，由于相變是從水滴底部開始逐漸向上發(fā)展，故可以看到水滴回縮時分成了上下兩層流動，下層有明顯的黏滯感甚至停止了運動，而上層仍保持有液體的流動性。

從親水性表面（θ=30°）的水滴形態(tài)對比圖中可以看到，在壁面溫度為-25℃時，水滴撞擊鋪展后幾乎沒有回縮過程，僅表層薄水膜有輕微振蕩，最后水滴以最大鋪展的形狀凍結(jié)到了表面上，這一結(jié)果也同文獻(xiàn)[14]中實驗觀察到的現(xiàn)象一致。同樣降至-25℃時，水滴碰撞結(jié)果由反彈變?yōu)橥Ａ粼诒砻嫔?。比較疏水表面（θ=114°）上的兩個回縮過程發(fā)現(xiàn)，常溫壁面上，水滴收縮時呈現(xiàn)外緣凸起中心凹陷的形態(tài)，來自三相交界線處的表面張力是促進(jìn)水滴回縮的主要動力；而在-25℃的表面上，由于水滴底部被凍住，表面張力變成了阻止水滴彈離表面的作用力，水滴上層的液體部分僅依靠慣性回縮至中心區(qū)域，產(chǎn)生彈起的趨勢，隨后被拉回表面上。若進(jìn)一步提高表面的疏水性，在超疏水（θ=163°）表面上，由于最大鋪展直徑減小，水滴鋪展過程進(jìn)一步縮短，摩擦及黏性造成的損耗更小。水滴在凍結(jié)之前保有一定動能，能夠產(chǎn)生部分回縮，縮小凍結(jié)面積，且最終水滴的上層液體有足夠的慣性脫離凍結(jié)的底部，實現(xiàn)部分反彈。

此外，從不同浸潤性表面上水滴最大鋪展發(fā)生的時刻可以看出，接觸角越大，最大鋪展發(fā)生的時刻越靠前，表明水滴同表面之間換熱的概率也越低?？偟膩碚f增強表面疏水性對于延緩或防止水滴凍結(jié)有兩個方面的優(yōu)勢:一方面是通過阻止水滴潤濕表面，減小其鋪展面積，降低了散熱速率；另一方面是減少水滴碰撞過程中的動能損耗，加快水滴回縮過程，縮短水滴與表面接觸的時間。故而在表面溫度低于冰點時，超疏水表面更能保持表面的潔凈。

3.2 表面浸潤性對水滴濕潤系數(shù)的影響

水滴與壁面之間的接觸面積決定了換熱量的大小，故而監(jiān)測水滴碰撞過程中潤濕面積的變化能夠反映表面浸潤性及表面溫度對水滴碰撞結(jié)冰過程的影響。首先，圖3比較了不同浸潤性和不同溫度下能夠達(dá)到的最大鋪展系數(shù)，其定義為水滴濕潤面積同初始時刻水滴在表面上的投影面積之比。由于相變過程是從水滴的最大鋪展?fàn)顟B(tài)前后開始，所以在表面溫度相同時，水滴的最大鋪展系數(shù)越大，水滴凍結(jié)的可能性就越高。從圖3中可以看出最大鋪展系數(shù)主要決定于表面浸潤性的大小，接觸角相同的表面上，最大鋪展系數(shù)的變化量在 5%以內(nèi)。超疏水表面上最大鋪展系數(shù)為親水性表面上的 1/2，即水滴的最大換熱面積可減少 1/2。此外，從圖 3中還可看出，在親水性表面（θ=30°）上，隨著壁面溫度降低，最大鋪展系數(shù)減小，表明水滴在完全鋪展之前就凍結(jié)在了表面上。而在超疏水表面（θ=163°）上，水滴最大鋪展系數(shù)隨著壁面溫度降低略有提高，表明在此情況下結(jié)冰發(fā)生在完全鋪展之后，由于水滴底部的相變影響了三相接觸線附近表面張力的作用，使得水滴收縮過程有所推遲。

圖4和圖5定量地比較了各個碰撞表面上水滴濕潤面積的動態(tài)變化過程。其中，不同過程中的濕潤面積以各自的最大濕潤面積進(jìn)行量綱1化，結(jié)果表明與文獻(xiàn)[14]實驗測得的量綱1濕潤面積變化趨勢一致，最終時刻的濕潤面積比相符。從圖4中可看到，水滴的鋪展過程隨著表面浸潤性的提高而延長。在θ=30°的親水性表面上，水滴持續(xù)鋪展至10 ms才開始回縮，且回縮過程緩慢，而疏水表面和超疏水表面上達(dá)到最大鋪展面積的時間依次為 6 ms和4 ms。當(dāng)表面降低至-15℃后，親水表面上水滴達(dá)到最大鋪展的時間略提前，且回縮過程消失；同樣地，疏水表面上的水滴最終產(chǎn)生部分凍結(jié)。但是超疏水表面上水滴的濕潤面積同常溫下的變化過程基本一致，這說明在-15℃時，超疏水表面能夠?qū)崿F(xiàn)防止水滴碰撞結(jié)冰。

圖3　不同溫度和浸潤性壁面上水滴最大鋪展系數(shù)的對比Fig.3 Comparison of maximum wetting factor with different temperatures and contact angles

圖4　-15℃壁面溫度下液滴潤濕因子的變化Fig.4 Variation of wetting area ratio on different wettability surfaces at -15℃

圖5　-25℃壁面溫度下液滴潤濕因子的變化Fig.5 Variation of wetting area ratio on different wettability surfaces at -25℃

當(dāng)表面溫度進(jìn)一步降低至-25℃時，如圖 5所示，水滴碰撞后的鋪展過程與常溫壁面下仍然保持一致，并且達(dá)到最大鋪展的時間也基本不隨表面溫度的變化而改變。但是降低表面溫度后，各個表面上的凍結(jié)面積均有所增加，在-25℃時，最終疏水表面上的凍結(jié)面積比為0.8左右，而超疏水表面上的凍結(jié)面積比大約為0.5。比較兩者的回縮過程發(fā)現(xiàn)，超疏水表面能夠延長回縮的時間，說明提高表面疏水性能夠減少水滴同表面間的換熱，且一定程度上使水滴保持了更高的動能以克服相變帶來的流動阻力。

值得注意的是，疏水表面上數(shù)值模擬的水滴收縮過程均比實驗結(jié)果偏快，這是由于本文研究中使用靜態(tài)接觸角模型表征表面浸潤性，忽略了水滴滯后角的影響。實驗中使用的疏水表面約有30°的滯后角，而最終模擬的結(jié)果表明，在滯后角較小（Δθ=5°）的超疏水表面上與實驗結(jié)果符合得更好。

3.3 疏水性表面的減冰/防冰機理

同親水性表面及弱疏水性表面相比，超疏水表面能夠進(jìn)一步減少水滴與壁面之間的接觸面積，降低碰撞過程中的黏性耗散，從而防止或者延緩水滴凍結(jié)在表面上。但是水滴碰撞冷表面并完全反彈，并不意味著在此過程中水滴內(nèi)部沒有相變發(fā)生。分析水滴內(nèi)的相變云圖對探究表面特性對水滴碰撞后動態(tài)行為的影響有重要意義。以-15℃的溫度下，3種不同浸潤性表面上的水滴碰撞相變過程為例，模擬結(jié)果表明隨著表面疏水性的提高，水滴碰撞后的結(jié)果逐漸由完全鋪展變成部分反彈，進(jìn)而是完全反彈，如圖6～圖8所示。

相變云圖能夠解釋這種變化產(chǎn)生的原因。從圖6中可以看到，在親水表面上，相變發(fā)生在水滴達(dá)到最大鋪展(10 ms)之前。由之前的研究可知，在水滴的邊緣區(qū)域，由水滴表面張力及界面曲率引起的內(nèi)部高壓是促使水滴回縮的主要動力。但在親水性表面上，由于接觸角很小，在三相接觸線附近沒有足夠的動力克服相變帶來的阻力，最后水滴剩余的動能被黏性振蕩消耗，直至靜止。而在疏水表面上，如圖 7所示，雖然相變也從較早的時刻t=5.0327 ms 發(fā)生，但是水滴邊緣區(qū)域提供了較強的表面張力促進(jìn)液體向中心區(qū)域回流。從12.992 ms 和27.835 ms的相變云圖中能夠看到，水滴底部雖然已經(jīng)完全變成了固體，紅色區(qū)域的液體體積分?jǐn)?shù)為 0，但其邊緣處被稍微拉離了壁面。若進(jìn)一步提高表面疏水性，則水滴回縮的動力會更強，從圖 8 中7.9915 ms時刻的相變云圖可以看到，水滴回縮過程中就可以將底部中間的部分相變區(qū)域卷起。在水滴離開壁面后，仍未見有完全固化的區(qū)域，表明在超疏水表面上，水滴與壁面接觸時間十分短暫，足以令水滴相變過程的時間長度不可忽略。由于同超疏水表面的接觸時間很短，水滴底部發(fā)生的相變并不完全，加上表面張力的作用較強，最終水滴剩余的動能足以將發(fā)生相變的部位裹挾起，共同反彈離開表面。

圖6　-15℃親水表面上水滴碰撞結(jié)冰過程Fig.6 Dynamic behavior and phase transition of droplet impinging on hydrophilic surface at -15℃

圖7　-15℃疏水表面上水滴碰撞結(jié)冰過程Fig.7 Dynamic behavior and phase transition of droplet impinging on hydrophobic surface at -15℃

圖8　-15℃超疏水表面上水滴碰撞結(jié)冰過程Fig.8 Dynamic behavior and phase transition of droplet impinging on super-hydrophobic surface at -15℃

4　結(jié) 論

（1）基于耦合VOF和Level-set的方法，結(jié)合焓-孔隙度模型，模擬了冷水滴碰撞不同溫度、不同浸潤性的表面之后的結(jié)冰過程，準(zhǔn)確捕捉到水滴碰撞過程中的動態(tài)行為，與實驗結(jié)果符合較好。且本方法能反映水滴內(nèi)部的相變細(xì)節(jié)，為研究表面浸潤性對水滴碰撞結(jié)冰的影響提供依據(jù)。

（2）從水滴內(nèi)部相變的發(fā)展過程分析不同浸潤性表面上水滴凍結(jié)情況，結(jié)果表明提高表面疏水性不僅能夠減小水滴碰撞后的鋪展?jié)駶櫭娣e，而且能夠縮短水滴同壁面的接觸時間，從而降低水滴內(nèi)部的相變速率。

（3）在水滴碰撞的回縮過程中，疏水或超疏水表面能夠產(chǎn)生較大的表面張力克服部分相變引起的黏滯阻力，從而起到減小凍結(jié)面積，降低結(jié)冰危害程度的效果。在高于某一溫度的范圍內(nèi)(-15℃)，超疏水表面能使水滴完全彈離表面。

（4）為提高數(shù)值模擬的準(zhǔn)確性，下一步工作可以關(guān)注在模型中考慮液體運動及相變引起的接觸角和物性參數(shù)的變化。

符 號 說 明

Cmush——黏糊系數(shù)

cp——比熱容，kJ·(kg·K)-1

d ——流場內(nèi)控制體到自由界面的距離

H ——焓，kJ·K-1

k ——熱擴(kuò)散系數(shù)，kJ·(K·m·s)-1

L ——潛熱，kJ·K-1

n,t ——分別為面積法向、面積切向

P ——壓力，Pa

T ——溫度，℃

Tsolid,Tliquid——分別為凝固溫度、融化溫度，℃

t ——時間，s

v ——速度，m·s-1

α ——體積分?jǐn)?shù)

γ ——液相體積分?jǐn)?shù)

θ ——接觸角，(°)

κ ——曲率，m-1

μ ——黏度，kg·(m·s)-1

ρ ——密度，kg·m3

σ ——表面張力，N·m-1

下角標(biāo)

P ——兩相流中的主相

S ——兩相流中的次相

w ——壁面

References

[1] LAMRAOUI F, FORTIN G, BENOIT R, et al. Atmospheric icing impact on wind turbine production[J]. Cold Regions Science and Technology, 2014, 100: 36-49. DOI: 10.1016/j.coldregions.2013. 12.008.

[2] BRAGG M B, BROEREN A P, BLUMENTHAL L A. Iced-airfoil aerodynamics[J]. Progress in Aerospace Sciences, 2005, 41(5): 323-362. DOI:10.1016/j.paerosci.2005.07.001.

[3] MORITA K, OKAMOTO K, AOKI A, et al. Hydrophobic coating study for anti-icing aircraft[R]. SAE Technical Paper, 2011.

[4] ANTONINI C, INNOCENTI M, HORN T, et al. Understanding the effect of superhydrophobic coatings on energy reduction in anti-icing systems[J]. Cold Regions Science and Technology, 2011, 67(1/2): 58-67. DOI: 10.1016/j.coldregions.2011.02.006.

[5] KIMURA S, YAMAGISHI Y, SAKABE A, et al. A new surface coating for prevention of icing on airfoils[R]. SAE Technical Paper, 2007.

[6] MESSINGER B L. Equilibrium temperature of an unheated icing surface as a function of air speed[J]. Journal of the Aeronautical Sciences, 1953, 20(1): 29-42. DOI:10.2514 /8.2520.

[7] BRAKEL T W, CHARPIN J P F, MYERS T G. One-dimensional ice growth due to incoming supercooled droplets impacting on a thin conducting substrate[J]. International Journal of Heat and Mass Transfer, 2007, 50(9/10): 1694-1705. DOI:10.1016/j.ijheatmasstransfer. 2006.10.014.

[8] JUNG S, DORRESTIJN M, RAPS D, et al. Are super-hydrophobic surfaces best for icephobicity?[J]. Langmuir, 2011, 27(6): 3059-3066. DOI:10.1021/la104762g.

[9] FEUILLEBOIS F, LASEK A, CREISMEAS P, et al. Freezing of a subcooled liquid droplet[J]. Journal of Colloid and Interface Science, 1995, 169(1): 90-102. DOI:10.1006 /jcis.1995.1010.

[10] 王皆騰, 劉中良, 勾昱君, 等. 冷表面上水滴凍結(jié)過程的研究[J].工程熱物理學(xué)報, 2007, 28(6): 989-991. DOI:10.3321/j.issn: 0253-231X.2007.06.028.

WANG J T, LIU Z L, GOU Y J, et al. Study of freezing process of water droplet on cold surface[J]. Journal of Engineering Thermophysics, 2007, 28(6): 989-991. DOI:10.3321/j.issn: 0253-231X. 2007.06.028.

[11] CHAUDHARY G, LI R. Freezing of water droplets on solid surfaces: an experimental and numerical study[J]. Experimental Thermal and Fluid Science, 2014, 57: 86-93. DOI:10.1016/j.expthermflusci. 2014.04.007.

[12] HINDMARSH J P, RUSSELL A B, CHEN X D. Experimental and numerical analysis of the temperature transition of a suspended freezing water droplet[J]. International Journal of Heat and Mass Transfer, 2003, 46(7): 1199-1213. DOI:10.1016/S0017-9310(02)00399-X.

[13] WANG H, HE G, TIAN Q. Effects of nano-fluorocarbon coating on icing[J]. Applied Surface Science, 2012, 258(18): 7219-7224. DOI:10.1016/j.apsusc.2012.04.043.

[14] MISHCHENKO L, HATTON B, BAHADUR V, et al. Design of ice-free nanostructured surfaces based on repulsion of impacting water droplets[J]. ACS Nano, 2010, 4(12): 7699-7707. DOI:10.1021/ nn102557p.

[15] LI H, ROISMAN I V, TROPEA C. Influence of solidification on the impact of supercooled water drops onto cold surfaces[J]. Experiments in Fluids, 2015, 56(6) DOI:10.1007/s00348-015-1999-2.

[16] BLAKE J, THOMPSON D, RAPS D, et al. Simulating the freezing of supercooled water droplets impacting a cooled substrate[J]. AIAA Journal, 2015, 53(7): 1725-1739. DOI:10.2514/1.J053391.

[17] 梁超, 王宏, 朱恂, 等. 液滴撞擊不同浸潤性壁面動態(tài)過程的數(shù)值模擬[J]. 化工學(xué)報, 2013, 64(8): 2745-2751. DOI: 10.3969/j.issn. 0438-1157.2013.08.006

LIANG C, WANG H, ZHU X, et al. Numerical simulation of droplet impact on surfaces with different wettabilities[J]. CIESC Journal, 2013, 64(8): 2745-2751. DOI: 10.3969/j.issn. 0438-1157.2013.08.006

[18] 鄭志偉, 李大樹, 仇性啟, 等. 液滴碰撞球形凹曲面復(fù)合 level set-VOF法的數(shù)值分析[J]. 化工學(xué)報, 2015, 66(5): 1667-1675. DOI: 10.11949/j.issn.0438-1157.20141116.

ZHENG Z W, LI D S, QIU X Q, et al. Numerical analysis of coupledlevel set-VOF method on droplet impact on spherical concave surface[J]. CIESC Journal, 2015, 66(5): 1667-1675. DOI: 10.11949/j.issn.0438-1157.20141116.

[19] 徐爽, 趙寧, 王春武, 等. 低速下水滴撞擊固體表面運動模擬[J].航空動力學(xué)報, 2013, 28(3): 695-700. DOI: 10.13224/j.cnki. jasp.2013.03.030.

XU S, ZHAO N, WANG C W, et al. Simulation of droplets impacting on solid surface at low velocity[J]. Journal of Aerospace Power,2013, 28(3): 695-700. DOI: 10.13224/j.cnki.jasp.2013.03.030.

[20] BRIONES A M, ERVIN J S, PUTNAM S A, et al. Micrometer-sized water droplet impingement dynamics and evaporation on a flat dry surface[J]. Langmuir, 2010, 26(16): 13272-13286. DOI:10.1021/ la101557p.

[21] STROTOS G, ALEKSIS G, GAVAISES M, et al. Nondimensionalisation parameters for predicting the cooling effectiveness of droplets impinging on moderate temperature solid surfaces[J]. International Journal of Thermal Sciences, 2011, 50(5): 698-711. DOI:10.1016/j.ijthermalsci.2010. 11.021.

[22] KUMAR A. Solidification of impinging molten metal droplet on a cold substrate[J]. International Journal of Mechanical Engineering and Robotics Research, 2014, 3(2): 86-497.

[23] SUSSMAN M, PUCKETT E G. A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows[J]. Journal of Computational Physics, 2000, 162(2): 301-337. DOI: 10.1006/jcph.2000.6537.

[24] 車得福, 李會雄. 多相流及其應(yīng)用[M]. 西安: 西安交通大學(xué)出版社, 2007: 260.

CHE D F, LI H X. Multiphase Flow and Its Applications[M]. Xi’an: Xi’an Jiaotong University Press, 2007: 260.

[25] JUNG S, TIWARI M K, DOAN N V, et al. Mechanism of supercooled droplet freezing on surfaces[J]. Nature Communications, 2012, 3: 615. DOI:10.1038/ncomms1630.

[26] VOLLER V R, PRAKASH C. A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems[J]. International Journal of Heat and Mass Transfer, 1987, 30(8): 1709-1719. DOI:10.1016/0017-9310(87)90317-6.

2015-12-14收到初稿，2016-04-26收到修改稿。

聯(lián)系人:常士楠。第一作者:冷夢堯（1991—），女，博士研究生。

Received date: 2015-12-14.

中圖分類號:TK 01; O 359.1

文獻(xiàn)標(biāo)志碼:A

文章編號:0438—1157（2016）07—2784—09

DOI:10.11949/j.issn.0438-1157.20151888

基金項目:國家重點基礎(chǔ)研究發(fā)展計劃項目（2015CB755803）；國家自然科學(xué)基金項目（11372026）。

Corresponding author:Prof. CHANG Shinan, sn_chang@buaa.edu.cn supported by the National Basic Research Program of China (2015CB755803) and the National Natural Science Foundation of China (11372026).

Numerical simulation of droplet impinging and freezing on cold surfaces with different wettability

LENG Mengyao, CHANG Shinan, DING Liang
(School of Aeronautics Science and Engineering, Beihang University, Beijing 100191, China)

Abstract:A strategy is presented to simulate the impact and solidification of the water droplet on different substrates. Simulations were performed using a coupled volume-of-fluid and level-set method to tracking the air-water interface and an enthalpy-porosity method to capture the phase transition. Three surface types were investigated: hydrophilic surface (contact angle 30°), hydrophobic surface (contact angle 114°) and superhydrophobic surface (contact angle 163°). The results showed that decreasing the wettability would reduce the contact time and area with the cold surface, and delay the freezing of the droplet. When temperature of the surface is higher than -15℃, the superhydrophobic surface can remain entirely ice-free under the simulated conditions. By comparing the maximum droplet spread, retraction response and time for solidification with the experimental results, the effectiveness and precision of the simulation strategy were demonstrated.

Key words:surface; droplet impact; dynamic behavior; multiphase flow; phase change; numerical simulation