結冰數值模擬中網格收斂性驗證

陳 航, 謝 亮, 李佩哲, 王福新, 劉 洪

(上海交通大學航空航天學院，上海 200240)

飛機穿越含有過冷水滴的云層時，撞擊在機身表面的水滴會凝結成冰。積冰會改變飛機的氣動外形，惡化機翼的氣動特性，嚴重危及整架飛機的飛行安全。隨著航空工業的發展，飛機性能越來越先進，但由結冰引發的飛行事故仍然屢見不鮮，必須對飛機結冰問題予以高度重視[1]。

20世紀80年代，歐美等國家開始采用數值模擬方法研究飛機結冰問題。結冰數值模擬是隨著計算流體動力學的發展而興起的，具有高效率、低成本、適用范圍廣等優點。近年來，隨著對飛機結冰物理過程認識的深入和數值計算方法的發展，數值模擬方法在研究飛機結冰問題中的應用越來越廣泛。現已開發出數個具有實用價值的結冰預測軟件，其中最為著名的是美國航空航天局開發的LEWICE軟件[2]。LEWICE軟件以及法國航空研究中心的結冰計算軟件[3]被美國聯邦航空局和歐洲聯合航空局所采用，作為確認飛機是否滿足結冰適航要求的工具。中國的許多團隊在這一領域也開展了相應研究。空氣動力研究與發展中心易賢等[4]、南京航空航天大學孫志國等[5]、西北工業大學周志宏等[6]、北京航空航天大學楊勝華等[7]和常士楠等[8]均采用數值模擬的手段對飛機結冰問題開展了深入研究。

在飛機結冰數值模擬中，保證網格收斂性是冰形計算結果可信的必要條件。對此問題，中國學者們已經進行了探討。蔣勝矩等[9]將二維可壓黏流Navier-Stokes方程的解用于翼型結冰的模擬中，對霜冰和明冰冰形進行了計算，并使用三套疏密程度不同的網格對霜冰冰形進行模擬，分析了霜冰情況下的網格收斂性；朱程香等[10]也驗證了霜冰情況下的網格收斂情況。前人對網格收斂情況的討論主要集中在霜冰情況，得出目前算法可以做到網格收斂的結論，然而對于情況更為復雜的帶角狀冰的混合冰及明冰情況，當前算法是否可以得到網格收斂的解還未見系統的討論。

由于角狀冰對氣動特性影響更為顯著，因此研究目前結冰算法對于角狀冰的網格收斂性能力尤為重要。現基于目前普遍使用的算法開發仿真程序，并使用溫度從低到高的四種算例詳細定量地驗證當前結冰模擬算法的網格收斂性情況。

1 結冰數值模擬方法

結冰數值模擬流程如圖1所示。在結冰數值模擬中，首先需要計算空氣流場，之后根據流場狀況進行水滴撞擊特性計算，使用拉格朗日法計算所有水滴粒子的運動軌跡從而得出壁面的水滴收集系數，再使用考慮粗糙度影響的S-A(spalart-allmaras)湍流模型計算壁面熱流，結合結冰熱力學模型計算得出壁面網格的結冰量，以此進行邊界重構并通過網格變形的方法更新網格，重復這個過程直到達到指定的結冰時間為止[11]。

圖1 結冰數值模擬流程

1.1 空氣流場計算

根據質量、動量、能量守恒關系建立二維非定常Navier-Stokes方程，控制方程為

(1)

式(1)中：Q為守恒向量；?Ω為區域Ω的邊界；Fc表示對流矢通量；Fv表示黏性矢通量；n表示邊界的外法向量。

采用有限體積法離散，對流通量計算采用Roe格式，黏性通量計算采用中心格式，時間推進采用LU-SGS(lower-upper symmetric gauss-seidel)隱式時間推進方法，湍流模型采用S-A(spalart-allmaras)模型。

1.2 水滴撞擊特性計算

使用拉格朗日法建立水滴運動軌跡方程[12]，控制方程如式(2)所示：

(2)

式(2)中：μ為動力黏度；ud為水滴速度；ua為流體速度；ρd為水滴密度；ρa為空氣密度；g為重力加速度；D為水滴直徑；CD為阻力系數；Re為相對雷諾數，表達式為

(3)

阻力系數CD采用圓球阻力模型計算，表達式為

(4)

通過求解水滴運動軌跡方程，可以得出水滴撞擊到壁面的位置，從而求得壁面的水滴收集系數。

1.3 結冰量計算

1.3.1 結冰熱力學模型

結冰計算基于經典的Messinger模型[13]，將結冰表面劃分為若干個控制體，對每一個控制體建立質量和能量守恒方程以求解結冰量。質量守恒方程和能量守恒方程分別如式(5)、式(6)所示[14]。

根據質量守恒，控制體內結冰的水量等于進入該控制體的水量減去離開該控制體的水量，即：

mice=mim+min-mout-mev

(5)

式(5)中：mice為結冰的水量；mim為撞擊到該控制體內的水量；min為從上一控制體流入的水量；mout為流出到下一控制體的水量；mev為蒸發或升華損失的水量。

根據能量守恒

eim+ein+qice=eout+eev+qc

(6)

式(6)中：eim為撞擊到控制體內的水帶來的能量；ein為從上一控制體流入的水帶來的能量；qice為結冰釋放的潛熱；eout為流出到下一控制體的水帶走的能量；eev為蒸發或升華的水帶走的能量；qc為對流換熱帶走的熱量。

1.3.2 對流換熱系數計算

在能量守恒方程中，最重要的熱流項是對流換熱項，因此，對流換熱系數的計算是結冰模擬中的關鍵問題。使用CFD (computational fluid dynamics)方法求解對流換熱系數通過湍流模型來完成，其中，Spalart等[15]提出的S-A模型具有計算簡單、與試驗結果吻合較好等優點。

S-A湍流模型輸運方程為

(7)

湍動黏度μt由式(8)求得：

(8)

式(8)中：fv1為黏性阻尼系數，定義為

(9)

式(9)中：Cv1為常數；χ為脈動速度與分子黏度的比值，即：

(10)

雷諾應力由式(11)計算：

(11)

式(11)中：ρ為流場當地密度；Sij為速度散度，ui、uj分別代表i、j方向的速度。

(12)

(13)

式中：d為到最近壁面的距離；S為渦度，定義為

(14)

式(14)中：Ωij為時均旋轉率張量，定義為

(15)

輸運方程中fw函數為

(16)

式(16)中：Cw3為常數；g的表達式如式(17)所示：

g=r+Cw2(r6-r)

(17)

式(17)中：r為脈動速度與到最近壁面距離關系的函數，其表達式如式(18)所示：

(18)

壁面邊界條件為

(19)

上述S-A湍流模型適合于光滑表面的湍流模擬。在結冰計算中，翼面由于積冰的影響會變得粗糙。由于表面粗糙組織對氣流的擾動作用，流動轉捩的位置將大大前移，壁面與氣流間的對流換熱也將被大大強化，因此，在結冰計算中必須考慮粗糙度對壁面對流換熱的影響。

Aupoix等[16]基于等效沙粒粗糙高度模型對S-A模型進行了擴展，得到了考慮粗糙度影響的S-A湍流模型的擴展模型。粗糙壁面的擴展模型在原始湍流模型的基礎上作如下改動。

向到最近壁面的距離d中加入了粗糙度的影響，d變為

d=dmin+0.03ks

(20)

式(20)中：ks為等效沙粒粗糙高度。

(21)

式(21)中：n表示邊界的外法向量。

調整χ的值以改變黏性阻尼函數fv1，χ改為

(22)

式(22)中：常數Cr1=0.5。

fv2函數變為

(23)

使用CFD方法求解流場后，在壁面附近，根據傅里葉定律可以求得壁面對流換熱量q：

(24)

根據對流換熱系數定義即可求出粗糙表面對流換熱系數：

(25)

式(25)中：Twall為壁面溫度；T∞為環境溫度。

1.3.3 結冰表面粗糙度計算

結冰表面粗糙度的產生機理非常復雜，與水膜的流動穩定性、水滴及氣流的擾動等很多因素都有關系，是一個耦合的多尺度問題，很難用解析方法模擬其特征，因此在結冰數值模擬中一般不對其做數值預測，而是利用試驗得到的經驗公式進行計算。采用Shin等[17]提出的等效沙粒粗糙高度模型計算粗糙度，將等效沙粒粗糙高度ks與結冰氣象條件中的液態水含量(liquid water content，LWC)、環境溫度T∞和水滴粒徑(median volume diameter，MVD)關聯起來，其中：

(26)

式(26)中：kLWC、kT∞和kMVD分別為表征含水量、環境溫度、水滴粒徑影響的函數，其表達式分別為

kLWC=0.571 4+0.245 7LWC+1.257 1LWC2

(27)

(28)

(29)

(30)

式(30)中：c為弦長。

1.4 網格變形

隨著冰的增長，外部流場和水滴收集系數等參數會隨著結冰外形變化而變化，為保證計算結果的準確性，需要對計算網格進行更新。網格變形技術實際上是將初始狀態下壁面的位移量通過插值的方法均勻地分布到流場內部的節點中，達到網格更新的效果。采用網格變形技術能夠最大限度地保留初始網格特征，在網格更新速度上比重新生成網格要快得多，而且，利用上一次結冰計算的流場值作為新網格流場計算的初始值，可以大幅減少網格更新后的流場計算時間[18]。

在網格變形開始之前需要判斷冰的生長方向，選取的冰生長方向如圖2所示[19]。

圖2 結冰生長方向[19]

節點A、B、C、D分別為原始機翼表面的網格節點，A＇、B＇、C＇、D＇分別為結冰后新的網格節點。圖2中的兩條虛線BB＇和CC＇分別為∠ABC和∠BCD的平分線，以虛線BB＇和CC＇ 為冰生長方向進行模擬。

之后根據Messinger結冰模型計算出每個網格節點的結冰高度，在下一時間步的結冰計算開始之前，采用徑向基函數插值網格變形技術來更新結冰翼型網格[20]，結果如圖3～圖6所示。

圖4～圖6是在干凈翼型網格的基礎上，根據冰生長方向和結冰高度，通過徑向基函數網格變形技術得到的變形之后的結冰網格，可以看出采用網格變形技術能夠最大限度地保留原始翼型網格的特征，在變形之后，網格質量仍然較好。

圖3 繞翼型網格

圖4 霜冰冰形網格

圖5 混合冰冰形網格

圖6 明冰冰形網格

2 網格收斂性驗證算例

對NACA0012翼型使用5套疏密程度不同的結構網格進行結冰計算，壁面首層網格的法向距離y+取1，不同網格對應的網格數量如表1所示。驗證算例來自文獻[21]，各算例計算條件如表2所示。

表1 網格編號與網格單元數量

表2 驗證算例計算條件

對一個霜冰算例、兩個混合冰算例以及一個明冰算例進行計算，環境溫度由250.37～265.37 K依次升高，計算結果如圖7～圖10所示，各結冰條件下采用不同網格計算得出的結冰冰形對比情況如表3～表6所示。

對于具有一個冰角的霜冰和明冰冰形，量化標準為最大結冰高度和與水平線的夾角，最大結冰高度為網格變形后移動距離最大的網格節點的移動距離，與水平線的夾角為網格變形后移動距離最大的網格節點與機翼前緣內切圓圓心的連線和水平線之間的夾角。

對于具有兩個冰角的混合冰冰形，量化標準為上冰角最大結冰高度、下冰角最大結冰高度以及冰形張角，冰形張角為上冰角最大高度的節點和機翼前緣內切圓圓心的連線與下冰角最大高度的節點和機翼前緣內切圓圓心的連線之間的夾角。

圖7 霜冰冰形計算結果(T=250.37 K)

圖8 混合冰冰形計算結果(T=262.04 K)

圖9 混合冰冰形計算結果(T=263.71 K)

圖10 明冰冰形計算結果(T=265.37 K)

表3 霜冰冰形對比情況(T=250.37 K)

表4 混合冰冰形對比情況(T=262.04 K)

表5 混合冰冰形對比情況(T=263.71 K)

由霜冰算例可知，計算結果與試驗結果吻合良好。由圖7、表3中的定量數據來看，隨著網格的加密，霜冰冰形逐漸收斂，且收斂的規律性良好。

由帶角狀冰的三個算例可知，計算結果與試驗結果大致吻合；從圖8～圖10及表4～表6的定量數據來看，隨著網格的加密，冰形基本逐漸收斂；但從定量數據來看，收斂的規律性不如霜冰情況。

表6 明冰冰形對比情況(T=265.37 K)

3 結論

(1)對于霜冰情況，目前采用的結冰模擬算法可以得到網格收斂的結果，這與其他文獻中得出的結論相符。

(2)對于混合冰和明冰情況，從結冰冰形來看，隨著網格的加密，其形狀趨于收斂；從冰角高度和張角角度來看，也可以得知混合冰和明冰是逐漸趨于收斂的，但收斂的規律性不如霜冰情況。

(3)目前普遍使用的結冰模擬算法具備得到網格收斂結果的能力，可用于工程實際，但對于帶有角狀冰的情況，在實踐中需要詳細謹慎地驗證結果的網格收斂情況。