考慮非平衡效應(yīng)的過冷水滴凝固特性

肖光明， 杜雁霞， 王橋， 郭龍， 王茂

中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點實驗室， 綿陽 621000

考慮非平衡效應(yīng)的過冷水滴凝固特性

肖光明， 杜雁霞*， 王橋， 郭龍， 王茂

中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點實驗室， 綿陽 621000

非平衡凝固是過冷條件下水滴凝固過程的重要現(xiàn)象。本文針對飛機(jī)結(jié)冰過程過冷水滴的非平衡凝固效應(yīng)，發(fā)展了改進(jìn)的凝固特性預(yù)測模型及數(shù)值計算方法，并自行搭建了實驗系統(tǒng)，開展了所建過冷水滴凝固模型與數(shù)值預(yù)測方法的實驗驗證。研究表明，所發(fā)展的改進(jìn)模型可有效表征水滴過冷階段的非平衡凝固效應(yīng)，因而對冷水滴凝固速率的預(yù)測有較好的改進(jìn)；當(dāng)過冷度為0 ℃時，過冷模型退化為傳統(tǒng)模型?；谒ǚ椒?，開展了過冷度及冷卻條件對水滴凝固特性的影響分析，獲得了不同條件下水滴凝固過程的溫度分布及相界面變化特征。研究表明，過冷度越大或水滴尺度越小，凝固速率相對越高；在考慮非平衡凝固效應(yīng)的條件下，過冷水滴凝固速率要高于不考慮非平衡凝固效應(yīng)的工況。相關(guān)研究可為結(jié)冰熱力學(xué)模型的改進(jìn)，以及結(jié)冰特性的精細(xì)化預(yù)測提供參考。

過冷水滴； 結(jié)冰； 非平衡效應(yīng)； 凝固； 改進(jìn)模型

飛機(jī)結(jié)冰是影響飛行安全的重要隱患之一，也是過冷水滴撞擊于飛機(jī)表面并發(fā)生凝固的一種特殊相變現(xiàn)象[1]。由于過冷條件的存在，因而飛機(jī)結(jié)冰具有顯著的非平衡凝固特征[2]。當(dāng)水滴低于凝固點溫度以液態(tài)形式存在時，往往會形成亞穩(wěn)平衡態(tài)[3-4]。此時，只要施加一個較小的擾動即可觸發(fā)凝固，并使其回到穩(wěn)定的平衡態(tài)[5]，而能量的波動、界面、雜質(zhì)、振動等均是觸發(fā)亞穩(wěn)態(tài)液體發(fā)生凝固的擾動源[5-6]。

過冷水滴的凝固通常分為2個典型階段[2,6]。第1階段由形核開始，水滴從熱力學(xué)非平衡態(tài)過渡到熱力學(xué)平衡態(tài)的階段，即枝晶形成階段，也有研究者稱其為部分凝固階段[7]。對于撞擊于飛機(jī)表面的過冷水滴，由于機(jī)體表面提供了異相形核的條件，因此在此階段形核過程從界面逐步發(fā)展到整個水滴，使水滴由液相變成冰水共存的模糊相，水滴溫度也由過冷態(tài)上升到凝固溫度所處的平衡態(tài)[6-7]。第2階段為完全凝固階段，即液固相界面由固相向液相推進(jìn)直至凝固完成的過程[7-8]。由于撞擊引起的異相形核作用，飛機(jī)結(jié)冰過程水滴凝固的第1階段顯著快于第2階段[2,6]。鑒于結(jié)冰物理過程的復(fù)雜性，目前大多關(guān)于過冷水滴結(jié)冰的研究均把重點放在凝固第2階段即相界面的推進(jìn)過程上。由于第2階段以平衡凝固為主，因此，多數(shù)研究者基于平衡凝固的相關(guān)理論與方法開展了結(jié)冰特性的預(yù)測研究, 如基于Enthalpy-Porosity法研究凝固過程的液/固相變行為[9-10]，但該方法主要針對平衡凝固過程，無法表征過冷水滴凝固的第1階段即非平衡凝固過程的影響特征。

近年來，隨著飛機(jī)結(jié)冰預(yù)測精度要求的提高，結(jié)冰過程的非平衡凝固效應(yīng)引起了研究者越來越多的關(guān)注。Worster[11]、Ellen[12]等提出了平面生長理論，針對飛機(jī)結(jié)冰凝固過程的2個階段，采用了不同的預(yù)測方法。Feuillebois等[6]研究了非平衡凝固對結(jié)冰第2階段的初始物性參數(shù)的影響特性。Blake等[2]在考慮形核過程的基礎(chǔ)上，基于FLUENT二次開發(fā)發(fā)展了過冷水滴凝固特性的預(yù)測方法?？梢钥闯觯捎谶^冷條件下水滴的非平衡凝固效應(yīng)對后期結(jié)冰速率及結(jié)冰特征有著重要影響而受到了越來越多的關(guān)注[13-14]。但由于凝固過程的復(fù)雜性，非平衡凝固特性預(yù)測方法的相關(guān)研究目前還較為缺乏，預(yù)測精度也有待提高，因而非平衡凝固規(guī)律特征的相關(guān)研究也較為薄弱。本文針對過冷水滴結(jié)冰的非平衡凝固效應(yīng)，發(fā)展了能表征非平衡凝固效應(yīng)的過冷水滴凝固特性預(yù)測方法，并自行搭建了實驗系統(tǒng)開展所建方法的實驗驗證。相關(guān)研究可為結(jié)冰熱力學(xué)模型的改進(jìn)，以及結(jié)冰特性的精細(xì)化預(yù)測提供參考。

1 考慮非平衡效應(yīng)凝固模型的建立

對單個過冷水滴凝固特性的研究有助于深入揭示飛機(jī)結(jié)冰過程過冷水滴凝固的物理特性[3]。為便于觀測和實驗，本文以單個水滴為對象，研究過冷水滴結(jié)冰過程非平衡凝固現(xiàn)象的共性特征。凝固實驗研究表明，過冷水滴的凝固過程由枝晶形成和相界面推進(jìn)2個典型階段構(gòu)成。在枝晶形成階段，水滴由透明態(tài)轉(zhuǎn)變?yōu)槟：龖B(tài)，如圖1(a)所示；在相界面推動階段，固/液相界面由冷卻面即底面向頂部平行推進(jìn)直至凝固完成，如圖1(b)～圖1(d)所示。

為建立相應(yīng)的數(shù)理模型，圖2顯示了過冷水滴凝固過程的簡化示意圖。

在凝固初始階段即枝晶形成階段，過冷水滴溫度由過冷態(tài)上升至凝固點，并伴隨潛熱的部分釋放，是典型的非平衡凝固階段；凝固第2階段為由熱擴(kuò)散驅(qū)動的界面推進(jìn)階段，也是在等溫下進(jìn)行的平衡凝固階段。盡管研究表明凝固第1階段相對于第2階段的時間較短，但由于結(jié)冰條件的不同，使第1階段結(jié)束時形成了第2階段的不同初始條件，從而影響了凝固的后續(xù)特征。因此，過冷水滴整個凝固過程的有效預(yù)測應(yīng)綜合考慮第1階段的非平衡凝固效應(yīng)及第2階段的平衡凝固效應(yīng)。針對過冷水滴的凝固兩個典型階段的特點，如何在相界面推進(jìn)過程的預(yù)測中考慮第1階段的非平衡效應(yīng)，是有效預(yù)測結(jié)冰全過程凝固特性需解決的重要問題。

圖1 過冷水滴凝固的典型階段
Fig.1 Typical freezing stages of supercooled water droplet

圖2 過冷水滴凝固過程的簡化示意圖
Fig.2 Schematic of supercooled droplet freezing process

考慮到在凝固第1階段完成并形成混合態(tài)的過程中過冷水滴已經(jīng)有了潛熱的部分釋放，因此，在凝固第2階段即相界面推進(jìn)過程的預(yù)測中應(yīng)考慮凝固第1階段的影響。將無量綱過冷度ε表示為

ε=cpsΔT/L

(1)

式中：ΔT=Tf-Tsupercooled，Tf為凝固溫度，Tsupercooled為過冷態(tài)溫度；L為相變潛熱；cps為固相定壓比熱容。

凝固第1階段結(jié)束時混合態(tài)中未發(fā)生相變的液相分?jǐn)?shù)可表示為

(2)

式中：cpl為液相定壓比熱容。

因此，混合態(tài)的液/固相變潛熱可表示為[6]

Lmix=Lf1

(3)

混合態(tài)的物性參數(shù)可表示為

ρmix=ρlf1+ρs1-f1

(4)

cpmix=cplf1+cps1-f1

(5)

λmix=λlf1+λs1-f1

(6)

式中：ρ、cp和λ分別為密度、定壓比熱容和熱導(dǎo)率；下標(biāo)mix、l和s則分別對應(yīng)混合態(tài)液相和固相物性參數(shù)。表1分別給出了液相水和固相冰對應(yīng)的物性參數(shù)[6]。

表1 冰和水的材料物性參數(shù)[6]Table 1 Physical properties of ice and water[6]

當(dāng)獲得凝固第1階段凝固過程的物性參數(shù)后，包含非平衡凝固的相變問題即可轉(zhuǎn)化為平衡凝固問題。為借鑒平衡凝固條件下液/固相變傳熱的預(yù)測方法，作以下幾點假設(shè)：

1) 由于撞擊形成的異相形核作用，凝固第1階段時間顯著小于第2階段。

2) 忽略水滴凝固過程因體積膨脹引起的水滴變形。

3) 界面推進(jìn)過程的液相區(qū)實質(zhì)為已在凝固第1階段發(fā)生了部分相變的混合態(tài)，因此凝固第2階段液相區(qū)物性參數(shù)由混合態(tài)參數(shù)代替。

基于上述假設(shè)并借鑒Enthalpy-Porosity方法，則凝固過程的傳熱控制方程可描述為

(7)

ρgβT-Tref+S

(8)

(9)

式中：u為液體流動速度矢量；p為液體壓力；h為熱焓；μ為液體黏性系數(shù)；g為重力加速度矢量；β為液體熱膨脹系數(shù)，且滿足Boussinesq近似；Tref為參考溫度；矢量S為凝固第2階段兩相區(qū)的源項，可描述為[10]

S=-C1-f2u

(10)

其中：C為界面推進(jìn)過程兩相共存區(qū)的特征參數(shù)；f2為結(jié)冰過程第2階段未凝固的液相分?jǐn)?shù)，可表示為

(11)

其中：Tl和Ts分別為液/固相變過程完全融化和完全凝固時的溫度。

2 過冷水滴凝固特性的數(shù)值分析

圖3顯示了直徑為4 mm、高度為2.5 mm，過冷度為-10 ℃條件下水滴凝固過程的溫度分布及相區(qū)分布特性預(yù)測結(jié)果，t為水滴凝固時間?？梢钥闯?，在相界面推進(jìn)的初始階段，在熱擴(kuò)散的作用下，相界面由固相區(qū)向液相區(qū)平行推進(jìn)。隨著凝固過程的發(fā)展，相界面與等溫線逐步向水滴頂部彎曲。

圖3 水滴凝固過程的溫度場及相界面變化特性
Fig.3 Temperature and interface evolution characteristics of droplet freezing process

圖4 過冷度對相界面變化速率的影響特性
Fig.4 Effect characteristics of supercooled temperature on freezing rate of phase interface

圖4顯示了過冷度分別為-15、-10、-5、0 ℃條件下相界面變化速率的比較(縱坐標(biāo)為相界面推進(jìn)高度ΔH與水滴半徑R的比值)。可以看出，在相同冷卻面溫度條件下，隨著水滴過冷度的增大，液/固相變的驅(qū)動力增加，相界面移動速度相應(yīng)增大；當(dāng)過冷度為0 ℃時，過冷模型退化為傳統(tǒng)Enthalpy-Porosity模型。

3 實驗驗證

為驗證所發(fā)展方法的有效性，本項目自行搭建了實驗系統(tǒng)并開展了相應(yīng)算法的實驗研究。實驗系統(tǒng)如圖5所示，由半導(dǎo)體制冷系統(tǒng)、Agilent 34970A多點數(shù)據(jù)采集儀、ANV TF100 PID溫度監(jiān)視器、溫度控制器、MotionXtra HG-100K高速攝像機(jī)、FLIR E60紅外熱像儀、LED無影光源系統(tǒng)、工控機(jī)及調(diào)壓電源組成。實驗時，采用滴管在冷表面產(chǎn)生不同尺度的水滴；開啟半導(dǎo)體制冷裝置，由PID溫度控制系統(tǒng)將制冷裝置冷卻面溫度維持在0 ℃以下實驗所需的溫度條件，使水滴緩慢冷卻至過冷態(tài)直至水滴完全凝固?；诮缑孀粉櫟姆椒?，采用高速攝像機(jī)記錄獲得水滴凝固過程相界面隨時間的變化特性；同時，采用紅外熱像儀記錄水滴凝固過程溫度的實時變化特性。每組工況重復(fù)3次，取每組平均值作最終實驗值。

在水滴凝固實驗中，在過冷期間的形核與枝晶生長階段，水滴逐步由透明態(tài)向模糊態(tài)過渡；當(dāng)水滴達(dá)到平衡溫度，固/液界面逐步出現(xiàn)并由液相區(qū)向固相區(qū)移動。如果以固/液界面出現(xiàn)的起始時刻作為凝固第2階段的起點，用高速相機(jī)即可記錄相界面隨時間的變化特征。圖6顯示了圖3相同條件下過冷水滴凝固第2階段相界面隨時間的變化特性?？梢钥闯觯嘟缑孀钕仍诶鋮s面產(chǎn)生，并逐步由固相區(qū)向液相區(qū)平行推進(jìn)，與計算相界面移動特性相似。

圖5 水滴凝固過程的實驗系統(tǒng)
Fig.5 Experimental setup for droplet freezing

圖6 實驗?zāi)踢^程相界面隨時間的變化特征
Fig.6 Experimental solidification evolution characteristics of phase interface vs time

圖7為本文發(fā)展的數(shù)值模型與傳統(tǒng)模型獲得的界面變化計算結(jié)果與實驗結(jié)果的比較。可以看出，在相同條件下，由于考慮了過冷水滴在過冷階段非平衡凝固效應(yīng)的影響，改進(jìn)的預(yù)測模型的計算與實驗結(jié)果的吻合程度優(yōu)于傳統(tǒng)Enthalpy-Porosity預(yù)測模型，且界面推進(jìn)速率要快于不考慮非平衡效應(yīng)的工況。改進(jìn)后的模型與實驗結(jié)果吻合得更好，表明了所建模型的有效性。

圖8顯示了半徑R分別為2 mm和3 mm兩種不同尺寸水滴，在過冷度為-10 ℃條件下的計算與實驗相界面移動特性的比較。從計算與實驗結(jié)果的比較可以看出，總體而言，計算與實驗相界面隨時間的移動特性吻合較好，但由于計算中忽略的凝固過程水滴體積膨脹效應(yīng)引起的變形及冒尖現(xiàn)象在凝固后期更為顯著，因而在凝固前期計算相界面與實驗相界面推進(jìn)速率的吻合程度要優(yōu)于后期。同時還可以看出，水滴尺度越小或冷卻面溫度越低，凝固速率相對越高。

為考查所建模型對凝固過程溫度變化特性預(yù)測的有效性，圖9顯示了半徑為2 mm的水滴在過冷度為-10 ℃條件下界面推進(jìn)過程的溫度變化特性紅外測試圖像及溫度曲線?？梢钥闯?，界面推進(jìn)階段，水滴由過冷態(tài)上升至平衡態(tài)，水滴溫度維持在0 ℃附近并保持相對恒定，當(dāng)凝固完成，水滴由平衡溫度逐步降低并趨于冷卻面溫度。

圖7 水滴凝固過程相界面變化特性的實驗與計算對比
Fig.7 Comparison of experimental and simulated phase interface evolution characteristics of droplet freezing

圖8 不同水滴直徑凝固過程相界面變化特性比較
Fig.8 Phase interface evolution characteristics of droplet freezing with different diameters

圖9 水滴凝固過程溫度變化特性
Fig.9 Temperature evolution characteristics of droplet freezing process

圖10 水滴凝固過程溫度變化特性的實驗與計算對比
Fig.10 Comparison of experimental and simulated temperature evolution characteristics of droplet freezing

圖10則顯示了過冷度分別為-10 ℃和-14 ℃ 條件下，直徑為4 mm水滴凝固界面推進(jìn)過程中水滴表面計算溫度與實驗紅外溫度變化特性的比較?？梢钥闯?，由于測量誤差及計算模型忽略變形效應(yīng)假設(shè)帶來的誤差，計算與實驗溫度變化歷程在凝固后期存在一定程度的偏差，但總體變化趨勢基本相似。且在凝固期間，與一般的液/固相變過程相同，計算與實驗水滴溫度均維持在平衡溫度即0 ℃附近，較好反映了水滴相變過程的恒溫特性。

4 結(jié) 論

1) 基于Enthalpy-Porosity模型，發(fā)展了過冷水滴結(jié)冰特性預(yù)測模型及數(shù)值計算方法，并自行搭建了水滴凝固實驗臺，開展了數(shù)值方法的實驗驗證。所發(fā)展方法計算結(jié)果與實驗結(jié)果吻合較好，當(dāng)過冷度為0 ℃時，過冷模型退化為傳統(tǒng)Enthalpy-Porosity模型，表明所建方法的有效性。

2) 相對于傳統(tǒng)方法，所發(fā)展的過冷水滴結(jié)冰特性預(yù)測模型及數(shù)值方法能有效表征非平衡凝固效應(yīng)，可用于過冷水滴凝固特性的預(yù)測，從而將傳統(tǒng)基于平衡凝固的Enthalpy-Porosity模型拓展至非平衡凝固研究領(lǐng)域。

3) 基于所建方法，開展了水滴凝固特性的影響因素分析，獲得了不同過冷條件下水滴凝固過程的溫度分布及相界面變化特征。研究表明，過冷態(tài)溫度越大、水滴尺度越小或冷卻面溫度越低，水滴的相變速率越高；在考慮非平衡凝固效應(yīng)的條件下，過冷水滴凝固速率要高于不考慮非平衡凝固效應(yīng)的工況。

本研究將進(jìn)一步考慮水滴結(jié)冰過程的體積膨脹效應(yīng)，為結(jié)冰熱力學(xué)模型的改進(jìn)以及結(jié)冰特性的精細(xì)化預(yù)測提供參考。

[1] 王洪偉, 李先哲, 宋展. 通用飛機(jī)結(jié)冰適航驗證關(guān)鍵技術(shù)及工程應(yīng)用[J]. 航空學(xué)報, 2016, 37(1): 335-350.

WANG H W, LI X Z, SONG Z. Key airworthiness validation technologies for icing of general aviation aircraft and their engineering application[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(1): 335-350 (in Chinese).

[2] 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.

[3] FUMOTO K, KAWANAMI T. Study on freezing characteristics of supercooled water droplets impacting on solid surfaces[J]. Journal of Adhesion Science and Technology, 2012, 26(4-5): 463-472.

[4] KING W D. Freezing rates of water droplets[J]. Journal of the Atmospheric Sciences, 1975, 32(2): 403-408.

[5] TABAKOVA S, FEUILLEBOIS F, RADEV S. Freezing of a suspended supercooled droplet with a heat transfer mixed condition on its outer surface[C]//1st International Conference on Applications of Mathematics in Technical and Natural Sciences, 2009, 1186(1): 240-247.

[6] 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.

[7] BURTNETT E. Volume of fluid simulations for droplet impact on dry and wetted hydrophobic and superhydrophobic surfaces[D]. Mississippi: Mississippi State University, 2012.

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

[9] VOLLER V, 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.

[10] CRANK J. Free and moving boundary problems[M]. Oxford: Oxford University Press, 1987.

[11] WORSTER M G. Solidification of fluids[M]. Cambridge: Cambridge University Press, 2000.

[12] ELLEN N, JACCO M H, EDWIN W, et al. Aircraft icing in flight: Effects of impact of supercooled large droplets[C]//29th Congress of the Aeronautical Sciences, 2014.

[13] TABAKOVA S, FEUILLEBOIS F. On the solidification of a supercooled liquid Droplet lying on a surface[J]. Journal of Colloid and Interface Science, 2004, 272(1): 225-234.

[14] SZIMMAT J. Numerical simulation of solidification process in enclosures[J]. Heat Mass Transfer, 2002, 38(4): 279.

(責(zé)任編輯： 李明敏)

URL：www.cnki.net/kcms/detail/11.1929.V.20161205.1640.002.html

Freezingcharacteristicsofsupercooledwaterdropletinconsiderationofnon-equilibriumeffect

XIAOGuangming,DUYanxia*,WANGQiao,GUOLong,WANGMao

StateKeyLaboratoryofAerodynamics,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China

Non-equilibriumeffectisanimportantphenomenoninfreezingofsupercooledwaterdropletinaircrafticingprocess.Basedontheenthalpy-porositymodel,anumericalpredictionmethodforfreezingofsupercooledwaterdropletisdeveloped.Theexperimentalsystemfordropletfreezingisbuiltupandseveralexperimentsareperformedtovalidatethenumericalmethodproposed.Theresultsindicatethatthedevelopedmodelisvalidandcanbeusedtopredictthefreezingcharacteristicsofsupercooledwaterdroplet.Basedontheimprovedfreezingmodel,theinfluenceofthedegreeofsupercooledandcoolingconditionsonthecharacteristicsofsupercooleddropletareanalyzed.Whenthedegreeofsupercooledisdecreasedtozero,thedevelopedmodeldegeneratestothetraditionalmodel.Thegreaterorsmallerthedegreeofsupercooledorsmallerthedropletsis,therelativelyhigherfreezingrateis.Inconsiderationoftheeffectofnon-equilibriumconditions,thefreezingandmovingrateofinterfaceishigherthanthetraditionalmodel.Relatedresearchcanprovideimportantreferenceforimprovingicingthermodynamicmodelandrefiningthepredictionmethodforicingaccretion.

supercooledwaterdroplet;icing;non-equilibriumeffect;solidification;improvedmodel

2016-08-24；Revised2016-10-25；Accepted2016-11-23；Publishedonline2016-12-051640

s：NationalNaturalScienceFoundationofChina(51308531,11672322);NationalBasicResearchProgramofChina(2015CB755800)

.E-mailyanxiadu@163.com

2016-08-24；退修日期2016-10-25；錄用日期2016-11-23； < class="emphasis_bold">網(wǎng)絡(luò)出版時間

時間：2016-12-051640

www.cnki.net/kcms/detail/11.1929.V.20161205.1640.002.html

國家自然科學(xué)基金 (51308531，11672322)； 國家“973”計劃 (2015CB755800)

.E-mailyanxiadu@163.com

肖光明， 杜雁霞， 王橋， 等. 考慮非平衡效應(yīng)的過冷水滴凝固特性J． 航空學(xué)報,2017,38(2):520703.XIAOGM,DUYX,WANGQ,etal.Freezingcharacteristicsofsupercooledwaterdropletinconsiderationofnon-equilibriumeffectJ.ActaAeronauticaetAstronauticaSinica,2017,38(2):520703.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0309

V211.3

A

1000-6893(2017)02-520703-07