翼面結冰過程中的冰晶運動相變與黏附特性

馬乙楗，柴得林，王強，易賢，*，孔滿昭

1．中國空氣動力研究與發展中心 結冰與防除冰重點實驗室，綿陽 621000

2．中國空氣動力研究與發展中心 空氣動力學國家重點實驗室，綿陽 621000

3．航空工業第一飛機設計研究院，西安 710089

航空發動機攝入冰晶易發生結冰，進而導致發動機動力損失或損壞［1-2］。相比于傳統過冷水結冰，冰晶結冰存在冰晶運動、相變、黏附等物理現象，結冰機理復雜［3-5］。由于研究冰晶結冰的實驗設備復雜、成本昂貴［6］，有必要開展冰晶結冰數值計算方法研究［7］。在冰晶結冰計算中，冰晶的運動相變決定了撞擊壁面的冰晶的形態與撞擊位置，黏附特性決定了參與結冰相變的冰晶的質量通量，二者直接關系到壁面冰形的增長與發展，是準確模擬冰晶結冰現象的關鍵環節［8］。因此，為深刻認識冰晶結冰現象，需要開展冰晶運動相變與黏附計算方法研究，系統分析冰晶運動相變與黏附特性。

圍繞冰晶運動相變以及冰晶黏附過程，國內外研究者開展了一系列研究，研究表明冰晶形狀對冰晶運動以及相變過程有顯著影響［9-11］，Norde等［12］模擬了冰晶運動過程中的傳熱過程，指出粒子運動和傳熱強烈依賴于冰晶形狀，在冰晶大量融化之前冰晶形狀通常為非球狀；Nilamdeen和Habashz［13］假設冰晶為平板狀，開展了冰晶的動力學特性研究，指出冰晶在流場的流動特性接近于圓盤狀顆粒；Iuliano等［14］假設冰晶為圓柱或圓盤狀，開展了冰晶阻力模型和傳熱過程的數值模擬；Trontin等［10］提出了冰晶的傳熱傳質模型，并對自由流中的非球狀冰晶顆粒傳熱傳質進行了研究，通過實驗所得冰形完善了冰晶黏附模型。郭向東等［11］假設冰晶為六角板狀和六棱柱狀，采用Villedieu等［15］提出的熱力學模型研究了冰晶形狀對其運動過程中換熱和軌跡的影響。針對冰晶黏附效率，Zhang等［16］通過比較冰晶粒徑和結構表面水膜厚度來判斷冰晶是否黏附或破碎，該過程在Fluent軟件上進行，計算得到黏附冰晶質量通量，進而計算出結冰冰形。Villedieu等［15］提出了基于結構表面水膜厚度、粒子表面水膜厚度、粒子等效粒徑和法向動能恢復系數的冰晶黏附效率計算方法，該方法計算冰形與Trontin等［10］實驗結果的吻合效果并不理想。Trontin等［10］提出了一種基于冰晶粒子融化率、總液態水含量和總水含量的冰晶黏附效率的計算方法，并根據NASA-NRC混合氣象條件結冰的實驗冰形數據校核了該黏附效率計算模型。綜上，國內外學者已經圍繞冰晶運動相變過程以及冰晶黏附過程建立了一系列數值計算模型，但在相關研究中，僅根據少數結冰實驗、冰形數據確定了模型系數并做了簡單模型驗證，未開展冰晶形狀對冰晶運動相變特性影響的深入研究，也未開展來流參數和冰晶融化效應對冰晶黏附特性影響的系統性分析與研究。

因此，本文建立了拉格朗日框架下冰晶運動-傳熱傳質耦合的數值計算方法，分析基于Monte Carlo方法的冰晶撞擊與黏附收集系數計算方法；依托NNWICE平臺［17］開發相應的計算程序，實現冰晶運動-傳熱傳質耦合與冰晶黏附的數值計算。以NACA0012翼型為對象，完成翼面冰晶/混合相結冰過程中，冰晶運動相變與黏附過程的模擬，研究冰晶形狀對冰晶運動相變特性的影響，系統分析不同來流溫度和液態水含量條件下冰晶的黏附特性，研究冰晶運動相變與黏附特性的影響規律。

1 計算方法

由于撞擊到翼型表面的冰晶不一定全部黏附，黏附冰晶量需在撞擊收集系數基礎上進一步計算冰晶黏附特性得到?；谏鲜龇治?，計算冰晶黏附特性的流程如圖1所示。首先，求解空氣流場，為冰晶運動相變、水滴運動計算提供背景流場；再通過冰晶運動-傳熱傳質計算得到冰晶撞擊特性與撞擊表面冰晶物態；最后，根據計算出的冰晶、水滴撞擊特性和冰晶物態完成冰晶黏附特性計算。采用NNW-FlowStar軟件完成空氣流場計算，空間離散采用一階迎風格式，通量格式采用Roe格式，湍流模型采用SST（Shear Stress Transport）k-ω模型。

圖1 冰晶黏附特性計算流程Fig. 1 Calculation process of ice crystal adhesion characteristics

1.1 冰晶運動-相變計算

1） 冰晶運動控制方程

根據Stokes定理、牛頓第二定律，在拉格朗日框架下的冰晶運動控制方程為

式中：up為冰晶運動速度；ua為空氣速度；CD為阻力系數；Rep為冰晶相對雷諾數；μa為空氣動力黏度；t為時間；ρa和ρp分別為空氣和冰晶粒子的密度；g為重力加速度；dp為冰晶等效體積直徑；系數

2） 冰晶形狀參數與阻力系數

采用以下形狀參數表征不同形狀冰晶［16］。

冰晶等效體積直徑（dp）表示具有和冰晶粒子相同體積的球狀冰晶的直徑，即

式中：Vp為冰晶粒子的體積。

冰晶球度（Φ）表示等效球狀冰晶表面積（Ap）與冰晶實際表面積（A）之比，即

冰晶橫向球度（Φ⊥）表示等效球狀冰晶的橫截面積（Ap，⊥）與冰晶流向投影的橫截面積（A⊥）之比，即

冰晶顆?？v橫比E（Aspect Ratio）表示冰晶軸向長度b（沿旋轉軸）與徑向長度a（垂直于旋轉軸）之比，即

通常冰晶形狀以板狀和柱狀為主［11］。為模擬對流云內的冰晶顆粒，表1給出了球狀和3種非球狀冰晶粒子的相關形狀參數計算公式［15］。

表1 冰晶形狀參數計算公式［15］Table 1 Calculation formulas for ice crystal shape parameters［15］

H?lzer和Sommerfeld［18］提出了一種考慮粒子方向并適用于所有形狀粒子的阻力系數表達式，本文采用式（7）計算式（1）中的冰晶阻力系數。

3） 冰晶相變模型

根據Trontin等［10］提出的相變模型，冰晶運動中的相變過程分為3個階段：① 過冷冰晶階段，冰晶溫度低于融化溫度，吸熱后溫度逐漸升高，僅發生升華；② 冰晶融化階段，冰晶溫度保持為融化溫度，吸熱后發生融化與蒸發，直至完全融化；③ 水滴蒸發階段，冰晶完全融化為水滴，溫度高于融化溫度，吸熱后溫度不斷升高，僅發生蒸發。

圖2展示了各階段冰晶狀態，其中Tp表示冰晶粒子溫度，Tm表示冰晶融化溫度。各階段傳熱傳質計算公式見文獻［10］。

4） 冰晶運動相變計算方法

冰晶運動方程與傳熱傳質方程采用一階歐拉向前時間離散格式進行求解，對于冰晶運動控制方程式（1），其離散形式為

圖2 各階段冰晶狀態示意圖Fig. 2 Schematic diagram of state of ice crystals at each stage

式中：n代表第n個時間步；Δt代表時間步長。而冰晶傳熱傳質方程根據冰晶所處狀態不同，具體公式見文獻［10］。

通常冰晶運動伴隨傳熱傳質發生，本文按以下計算方法求解冰晶運動-相變耦合模型：① 初始化冰晶參數（冰晶初始位置、粒徑、形狀、初始溫度）；② 采用徑向基函數插值方法由冰晶空氣流場變量插值得到冰晶所在位置處物理量；③ 求解冰晶相變方程更新冰晶質量、形態、球度、溫度等參數與變量；④ 將與求解運動方程相關的參數變量代入運動方程，計算更新冰晶位置、速度等變量；⑤ 重復步驟2~步驟4直至冰晶撞擊壁面/飛出指定區域或完全蒸發，完成軌跡求解。

1.2 冰晶黏附特性與收集系數計算

1） 冰晶撞擊收集系數計算

采用基于統計原理的Monte Carlo方法［19-20］計算冰晶撞擊收集系數，該方法既能適應復雜三維外形，又能克服軌跡交叉、冰晶運動過程中質量損失等問題。

在Monte Carlo方法下，冰晶撞擊收集系數（βs，MC）的計算式為

式 中：Ms為 網 格 單 元s從τ0到τtotal時 刻 的 質 量 流率積分；nτ為τ時刻撞擊到s上的冰晶數量；mˉp為單個粒子的質量；T為總時間步數；As為表面單元s的面積；u∞表示來流速度；IWC表示來流冰水含量。

2） 冰晶融化率與黏附效應模型

定義mp，i為冰核質量即冰晶中冰的質量，則融化率ηm可定義為

式中：m表示冰晶總質量。

文獻［21］提出了將黏附效率系數εs作為冰晶黏附效應的數學描述，其定義為黏附的冰晶質量與撞擊冰晶質量之比，并給出了計算冰晶黏附效率的一種計算方式。首先定義當空氣中只有冰晶以及融化冰晶時的黏附效率系數εs，c，其計算式為

式中：Kc為比例常數，通常取Kc=2.5。

定義空氣中同時含有冰晶和水滴時的黏附效率系數εs，d，計算式為

式中：Kd為比例常數，其計算式為

式中：Kd*=0.6；Tw為壁面處溫度；Tf為冰晶融化溫度，根據水的三相圖可知在來流入口壓力為98 000 Pa時，Tf=273.152 75 K；T*為冰晶完全無法黏附的閾值溫度，取文獻［21］中的建議值273.15 K。

φc表示撞擊冰晶的質量分數，φd表示撞擊水滴的質量分數，其計算式為

式中：βimp，c、βd分別表示冰晶、水滴撞擊局部收集系數； LWC表示來流液態水含量。

對式（12）~式（15）進行推導，可得

根據εs，c和εs，d，給出冰晶黏附效率εs和冰晶黏附收集系數βcatch的計算公式為

此外，式（13）表明壁面溫度也會對冰晶黏附率產生影響。為探究融化率與來流液態水含量對冰晶黏附收集系數的影響，本文假設翼面溫度為恒定值Tw=273.5 K，收集過程均在Tf

2 冰晶運動相變特性

2.1 冰晶形狀對運動軌跡的影響

通過球度對冰晶形狀特征進行數學描述，本節將討論不同來流條件下，冰晶球度對冰晶運動軌跡的影響。

以二維NACA0012翼型作為結冰對象，翼型弦長為1 m，參照文獻［13］，按照表2設置來流參數并開展計算與分析，表中MVD表示冰晶粒子的平均容積直徑。相同最大尺度（軸向長度與徑向長度的較大者）的不同形狀冰晶在運動過程中的沿程球度與運動軌跡結果如圖3~圖6所示。

當來流總溫為?5 ℃時，由圖3可看出，冰晶球度在運動過程中保持恒定。這是由于冰晶未發生融化。圖4展示了同一位置釋放的4種形狀冰晶的運動軌跡，3種非球狀冰晶運動軌跡較為接近且與球狀差距較大。這是由于球狀冰晶的質量與球度都明顯大于非球狀冰晶。圖5給出了來流總溫35 ℃冰晶沿程球度的變化。由圖5可知，非球狀冰晶融化前球度不變，開始融化后球度沿程增大，長橢球狀冰晶球度將逐漸增大至1（融化為水滴）；六角平板狀以及扁橢球狀冰晶球度接近且變化趨勢基本一致。圖6展示了來流總溫35 ℃下同一位置釋放4種冰晶的運動軌跡，與圖4相似，非球狀冰晶的運動軌跡接近且與球狀冰晶的差距明顯。算例1、算例2的數值實驗結果證明，冰晶形狀將對冰晶運動軌跡產生顯著影響。

表2 工況參數Table 2 Running parameters

圖3 來流總溫?5 ℃下冰晶沿程球度Fig. 3 Sphericity of ice crystals along trajectory at total temperature of ?5 ℃

圖4 來流總溫?5 ℃下冰晶運動軌跡Fig. 4 Trajectories of ice crystals at total temperature of ? 5 ℃

圖5 來流總溫35 ℃下冰晶沿程球度Fig. 5 Sphericity of ice crystals along trajectory at total temperature of 35 ℃

圖6 來流總溫35 ℃下冰晶運動軌跡Fig. 6 Trajectories of ice crystals at total temperature of 35 ℃

2.2 冰晶形狀對冰晶融化率的影響

采用表2所示的工況參數，僅改變來流總溫為15 ℃、35 ℃設置算例3和算例4。計算得到相同最大尺度的不同形狀冰晶的沿程溫度及融化率變化曲線如圖7、圖8所示。

圖7 來流總溫15 ℃下冰晶粒子溫度及融化率變化曲線Fig. 7 Variation of particle temperature and melting ratio at total temperature of 15 ℃

圖8 來流總溫35 ℃下冰晶粒子溫度及融化率變化曲線Fig. 8 Variation of particle temperature and melting ratio at total temperature of 35 ℃

圖7中不同形狀的冰晶運動過程中溫度與融化率變化曲線差異明顯，長橢球狀冰晶溫升與融化最快，其次是扁橢球狀和六角平板狀冰晶，球狀最慢；所有冰晶撞擊到壁面時均未完全融化，而球狀冰晶剛開始融化。這是由于球狀冰晶的質量明顯大于其余冰晶，升溫所需外界能量更多。由于扁橢球狀與六角平板狀冰晶形狀和球度類似，故其溫升與融化速率差異較小。

圖8給出了來流總溫35 ℃時4種冰晶粒子溫度及冰晶融化率變化曲線，冰晶溫度與融化率變化趨勢與算例3基本相同。由于來流總溫更高，冰晶粒子換熱更劇烈，升溫與融化速率顯著增加；相同冰晶撞擊到壁面時，算例4中融化率較算例3增大；長橢球狀冰晶已完全融化。

算例3和算例4的數值實驗結果證明，冰晶形狀對運動過程中冰晶溫度與融化率具有顯著影響，這將影響冰晶黏附特性及后續冰晶結冰相變過程中撞擊到結構表面的冰晶粒子的顯熱、物態進而影響冰形計算；此外，在最大尺度相同的條件下，相較于扁橢球狀和六角平板狀冰晶，長橢球狀冰晶在運動過程中換熱、融化現象更顯著。

3 冰晶黏附特性

云霧場通常是指包含由液態水滴和固態冰晶形成的云霧的空間區域，云霧場性能主要由以下云霧參數衡量：粒子的平均容積直徑（MVD）、液態水含量（LWC）、冰水含量（IWC）、總水含量（TWC）等，上述參數是影響翼面結冰的重要因素。根據Currie等［22］的實驗結果，冰晶黏附效率取決于撞擊粒子的液態水含量與總水含量的比值。通常LWC來自于2部分：冰晶融化與混合相來流中的過冷水滴。為分析冰晶融化以及過冷水滴含量對冰晶收集系數的影響，將通過改變來流溫度以及云霧參數開展數值計算研究。

3.1 來流溫度對黏附特性的影響

采用150 μm六角平板狀冰晶，參照表2設置其余工況參數，并改變來流總溫分別為?5、5、15、35 ℃，對算例5~算例8進行數值實驗。

圖9給出了不同溫度條件下的冰晶撞擊收集系數和冰晶黏附效率（εs）。圖中S表示翼型上下翼面的點到翼型前緣駐點的距離，L表示特征長度，此處為翼型的弦長。從圖9（a）可看出來流總溫為?5 ℃和5 ℃下的冰晶撞擊收集系數（βimp，c）曲線接近，隨著來流總溫升高冰晶的撞擊極限和撞擊收集系數峰值均增大。?5 ℃與5 ℃的冰晶均未進入融化階段，質量損失均很小，撞擊收集系數曲線相似。隨著來流總溫升高至冰晶融化，冰晶在靠近翼型表面繞流時隨流性變弱，冰晶撞擊收集系數峰值、撞擊極限增大。

圖9 冰晶撞擊收集系數與黏附效率曲線Fig. 9 Ice crystal impingement coefficient and sticking efficiency

由圖9（b）可知算例5與算例6的黏附效率曲線接近且呈鐘形。算例7與算例8冰晶黏附效率的范圍以及峰值較算例5與算例6顯著增大且維持在某一定值附近。對于算例5與算例6，由于冰晶并未融化，黏附效率完全取決于水滴收集量，曲線呈鐘形。隨著來流總溫的升高冰晶開始融化，這時黏附效率大小主要取決于融化率（ηm）和水滴收集系數（βd），其值取水滴對應黏附效率（εs，d）與 冰 晶 融 化 對 應 黏 附 效 率（εs，c）中 較 大 值。故算例7與算例8中冰晶黏附效率范圍與冰晶撞擊收集系數范圍一致。此外，算例7中靠近前緣附近水滴導致的黏附效應相較于冰晶融化導致的黏附效應更明顯，該區域的黏附效率明顯增高。

圖10給出了?5 ℃、15 ℃和35 ℃下冰晶、水滴收集系數和冰晶黏附收集系數（βcatch）曲線，圖10（d）給出了3種來流總溫條件下冰晶黏附收集系數曲線。選取收集系數曲線主要特征量——收集系數峰值和撞擊上極限弧面坐標值，并進行對比，表3給出了相應的數據，表中Tt∞為來流總溫，峰值降幅指βcatch，峰值相對于βimp，c峰值的降幅，上極限降幅指黏附上極限弧面坐標相對于撞擊上極限弧面坐標的降幅。

圖10 混合相條件下冰晶和水滴收集系數曲線Fig. 10 Ice crystal and droplet collection coefficient under mixed phase conditions

表3 不同來流總溫收集系數曲線主要特征對比Table 3 Comparison of main characteristics of collection coefficient at different total temperatures

如圖10（a）所示，在來流總溫?5 ℃下冰晶不融化，冰晶僅在液滴撞擊區域內發生黏附，黏附收集系數極限與水滴的撞擊極限相同。從圖10（b）可看出來流總溫15 ℃冰晶融化，冰晶撞擊收集系數峰值與黏附效率均增大，故其黏附收集系數的峰值增大；冰晶黏附收集系數極限與冰晶撞擊極限一致。如圖10（c）所示，來流總溫為35 ℃時，冰晶進一步融化，此時融化導致的黏附效應占主導，黏附效率值趨于1，故冰晶撞擊收集系數曲線與冰晶黏附收集系數曲線接近。由圖10（d）可知，冰晶黏附收集系數峰值與極限隨著來流總溫的升高呈增大趨勢。分析表明冰晶融化對冰晶黏附有顯著影響。

3.2 來流IWC和LWC對黏附特性的影響

針對混合相條件下，來流云霧參數對冰晶黏附效率和黏附收集系數產生的影響展開研究與討論。

由于不同來流總溫條件下，冰晶融化程度存在差異，這對冰晶黏附收集系數結果產生影響。為去除該影響因素，參照表2中算例1設置工況參數，來流總溫取?5 ℃保持不變，但改變來流LWC和IWC值 設 置 算 例9：LWC = 0.3 g/m3，IWC = 0.3 g/m3；算 例10：LWC = 0.3 g/m3，IWC = 0.7 g/m3；算 例11：LWC = 0.7 g/m3，IWC = 0.3 g/m3；算 例12：LWC = 0.7 g/m3，IWC = 0.7 g/m3；并進行數值實驗。

4種工況下計算結果中冰晶撞擊上極限弧面坐標值均為0.052 68。冰晶黏附上極限弧面坐標值均為0.037 2，這是由于冰晶均未融化，冰晶黏附極限均與水滴撞擊極限一致。表4給出撞擊與黏附收集系數峰值的對比結果。

表4 不同來流云霧條件下收集系數曲線主要特征對比（算例9~算例12）Table 4 Comparison of main characteristics of collection coefficient under different conditions of cloud field （Case9-Case12）

?5 ℃來流總溫下冰晶黏附效率和冰晶黏附收集系數計算結果如圖11所示。圖11顯示，隨著LWC/TWC比值增大，冰晶黏附效率和黏附收集系數峰值逐漸增大。當LWC/TWC值一定時（算例9與算例12），在來流總溫一致的情況下，冰晶將有著同樣的黏附效率和黏附收集系數。這是因為冰晶不融化時，黏附效率取決于撞擊壁面處水滴質量分數大小。

為進一步研究冰晶融化時LWC和IWC對冰晶黏附效率與冰晶黏附收集系數的影響，來流其余參數設置參照算例9~算例12，改變來流總溫為15 ℃，分別對算例13~算例16進行數值實驗。

算例13~算例16計算結果中冰晶撞擊上極限與冰晶黏附上極限的弧面坐標值均為0.073 93，這是由于4種工況下冰晶均融化，冰晶黏附極限均與冰晶撞擊極限一致。表5給出了收集系數曲線特征對比。

圖12給出了15 ℃條件下不同IWC與LWC下冰晶黏附效率和黏附收集系數的計算結果。如圖12（a）所示，除算例15外，其余算例的黏附效率曲線差距很小，整體保持在0.4左右，算例15中靠近翼型前緣位置冰晶黏附效率相較于其他算例明顯增高，由圖12（b）可得相似結果；此外，圖中各算例的冰晶黏附收集系數曲線極限均與冰晶撞擊收集系數極限一致。

圖11 來流總溫?5 ℃不同云霧參數下黏附效率與黏附收集系數Fig. 11 Sticking efficiency and catching coefficient with different cloud field parameters of total temperature of incoming flow ?5 ℃

表5 不同來流云霧條件下收集系數曲線主要特征對比（算例13~算例16）Table 5 Comparison of main characteristics of collection coefficient under different conditions of cloud field （Case13-Case16）

圖12 來流總溫15 ℃不同云霧參數下黏附效率與黏附收集系數Fig. 12 Sticking efficiency and catching coefficient with different cloud field parameters of total temperature of incoming flow 15 ℃

由于來流總溫為15 ℃時冰晶融化，冰晶黏附收集系數極限范圍與冰晶撞擊收集系數一致。對于算例15外的算例，撞擊水滴的質量分數較小，冰晶融化導致的黏附效應占主導，故具有相同的黏附效率曲線和黏附收集系數曲線；而算例15的來流LWC/TWC值較大，在翼型前緣位置附近收集到的液態水滴質量分數較高，水滴對應 黏 附 效 率（εs，d）大 于 冰 晶 融 化 對 應 黏 附 效 率（εs，c），故在該區域冰晶黏附效率和黏附收集系數相較于其余算例有明顯增高。

4 結 論

基于NNWICE平臺，對冰晶/混合相結冰過程中冰晶運動相變與黏附的計算方法與影響因素開展了研究，得到以下結論：

1） 建立了拉格朗日框架下冰晶運動-傳熱傳質耦合的數值計算方法，分析了冰晶幾何形狀對冰晶運動相變特性的影響，發現同粒徑條件下類似柱狀的長橢球狀冰晶相較于扁橢球狀冰晶和六角平板狀冰晶在運動過程中換熱、融化現象更顯著。

2） 分析了基于Monte Carlo方法的冰晶/混合相結冰黏附收集系數計算方法。該方法可以在考慮冰晶運動過程中質量損失的情況下，有效模擬冰晶收集系數。

3） 系統分析了來流溫度等云霧參數對冰晶黏附特性的影響。來流總溫升高，冰晶融化率增加，撞擊表面微元液態水質量分數增大，冰晶更易黏附；來流LWC/TWC增大，撞擊表面微元液態水質量分數增大，冰晶更易黏附。