不同浸潤性冷表面上水滴碰撞結冰的數值模擬

冷夢堯，常士楠，丁　亮

?

研究簡報

不同浸潤性冷表面上水滴碰撞結冰的數值模擬

冷夢堯，常士楠，丁亮

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

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

關鍵詞:表面；液滴碰撞; 動力學; 多相流；相變；數值模擬

引 言

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

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

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

1　數學模型

1.1 自由界面追蹤模型

本文采用VOF耦合Level-set方法[23]追蹤氣液自由界面，流場求解以 VOF的思想為主，在考慮界面效應時，以Level-set方法為輔助。VOF方法能夠保持追蹤流體的“質量”，并且能夠容易地追蹤拓撲結構變化的流體界面。在 VOF方法中，相與相之間的界面追蹤是通過計算次相的體積分數連續方程完成的。由于本文的研究對象為兩相流，主相為空氣，次相為水，故而只需要求解次相體積分數的連續方程

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

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

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

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

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

1.2 流動計算

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

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

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

1.3 凝固-融化模型

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

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

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

焓-孔隙度相變模型將冰水混合區域視為多孔結構，孔隙率即為液相的體積分數γ，取值區間為[0,1]，其與溫度的關系為

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

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

從動量源項的表達式可以看出，在液相區，即γ=1，動量源項為0，動量方程與流體N-S方程形式一致；凝固開始變成冰-水混合物后，液相體積分數γ降低，源項在動量方程中的比重逐漸超過瞬態項、對流項和擴散項，直至完全轉變為固相后，速度基本為零。

圖1　計算區域Fig.1 Droplet impingement case setup

2　計算模型

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

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

表1　水滴物性參數Table 1 Material properties for water

3　結果與討論

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

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

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

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

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

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

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

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

3.2 表面浸潤性對水滴濕潤系數的影響

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

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

圖3　不同溫度和浸潤性壁面上水滴最大鋪展系數的對比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℃

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

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

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

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

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

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

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

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

4　結 論

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

（2）從水滴內部相變的發展過程分析不同浸潤性表面上水滴凍結情況，結果表明提高表面疏水性不僅能夠減小水滴碰撞后的鋪展濕潤面積，而且能夠縮短水滴同壁面的接觸時間，從而降低水滴內部的相變速率。

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

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

符 號 說 明

Cmush——黏糊系數

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

d ——流場內控制體到自由界面的距離

H ——焓，kJ·K-1

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

L ——潛熱，kJ·K-1

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

P ——壓力，Pa

T ——溫度，℃

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

t ——時間，s

v ——速度，m·s-1

α ——體積分數

γ ——液相體積分數

θ ——接觸角，(°)

κ ——曲率，m-1

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

ρ ——密度，kg·m3

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

下角標

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] 王皆騰, 劉中良, 勾昱君, 等. 冷表面上水滴凍結過程的研究[J].工程熱物理學報, 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] 梁超, 王宏, 朱恂, 等. 液滴撞擊不同浸潤性壁面動態過程的數值模擬[J]. 化工學報, 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] 鄭志偉, 李大樹, 仇性啟, 等. 液滴碰撞球形凹曲面復合 level set-VOF法的數值分析[J]. 化工學報, 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].航空動力學報, 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] 車得福, 李會雄. 多相流及其應用[M]. 西安: 西安交通大學出版社, 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收到修改稿。

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

Received date: 2015-12-14.

中圖分類號:TK 01; O 359.1

文獻標志碼:A

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

DOI:10.11949/j.issn.0438-1157.20151888

基金項目:國家重點基礎研究發展計劃項目（2015CB755803）；國家自然科學基金項目（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