翼身融合布局與傳統布局飛機結冰特性對比分析

王灝東，桑為民，邱奧祥，李棟

西北工業大學，陜西 西安 710072

結冰是危害飛機飛行安全的重要因素，由結冰導致的飛行事故常常發生。機翼和發動機[1]等部位的結冰問題會對飛機的氣動特性產生劇烈影響，從而降低飛行質量，情況嚴重時甚至會導致機毀人亡等危害，因而受到人們的廣泛關注。目前國內外的研究人員致力于研究過冷大水滴和冰晶結冰[2]，為防/除冰方案設計提供了可靠依據。

對飛機結冰特性的分析研究方法主要分為數值模擬方法和試驗方法。數值模擬方法[3-4]具有計算效率和精度較高、周期短、容易操控計算條件、經濟性高等優點，是研究飛機結冰的重要手段。國內外研究人員在結冰數值模擬領域取得了眾多成果。V.H.Gray[5-6]研究了NACA 65A004翼型結冰過程，分析了帶冰翼型氣動性能，并提出了預測積冰產生的阻力的經驗公式。C.D.Macarthur[7]發展了一種數學模型，可用來計算二維翼型上的霜冰和光冰增長。M.B.Bragg[8]提出了一種能計算水滴軌跡的方法，并對該方法的改進提出了一些建議。C.E.Porter[9-10]研究了D8飛機上影響水滴撞擊特性的因素。譚燕[11]采用歐拉方法對某對稱楔形翼型結冰過程進行數值模擬，采用Spalart-Allmaras（SA）湍流模型獲得流場結果，應用歐拉方法獲得冰晶和液滴軌跡結果，并基于Messinger 模型來獲得冰形，之后通過NASA-NRC 第139 號試驗結果證實了該方法的可行性。張強等[12]利用歐拉法研究了ONERA M6 三維機翼表面的水滴收集系數，將結冰問題拓展到了三維。

本文研究工作對傳統布局飛機和新型翼身融合布局飛機的空氣流場計算結果與試驗結果進行了驗證，并將兩種布局飛機進行了對比研究分析?；谔岢龅慕Y冰數值模擬計算方法，對新型翼身融合布局構型的N2A翼身融合體和傳統布局的DLR-F4飛機在典型霜冰結冰條件下進行結冰數值模擬，分析結冰特性，并將兩者進行對比研究分析，初步探討新型翼身融合布局飛機和傳統布局飛機結冰特性的差異。

1 結冰數值模擬計算方法

1.1 空氣流場數值求解

空氣流場計算是研究結冰數值計算的基礎環節，也是研究水滴撞擊特性的關鍵。對結冰空氣流場的計算基于對計算域中Navier-Stokes（N-S）方程組的求解?？諝饬鲌鲈谶M行流動和傳熱時，可用控制方程對質量守恒定律、動量定律和能量守恒定律三個方程進行描述。N-S 方程積分形式為

式中，Q為守恒變量；F為對流通量；G為黏性通量；Ω為控制體。

1.2 水滴撞擊特性計算

水滴對飛機表面的撞擊特性是指水滴對飛機表面的撞擊區、撞擊量，以及水滴在撞擊區內的分布。過冷水滴運動方程的建立方法主要分為拉格朗日法和歐拉法[13-14]。本文基于歐拉法開展研究。

歐拉法的思想是將水滴與空氣看作兩相流建立控制方程。定義一個水滴容積參數α，建立水滴場的控制方程，其形式上與流場計算的控制方程可保持一致[15]。對于三維計算來說，采用歐拉法計算水滴撞擊特性不必在計算區域進行多次插值計算以及數值積分計算，并采用流場計算所用的網格，使水滴場的求解與流場求解在形式上統一，提高了代碼開發效率和計算效率[16]。歐拉法將云層中離散的過冷水滴當作連續相處理，并做出相關假設[17]。

歐拉法求解水滴運動特性的連續方程主要有以下幾個。

空氣連續方程為

水滴連續方程為

式中，ρ為過冷水滴的密度；α為體積因子。

收集系數β為微元表面上水滴的實際收集率與該微元表面上最大可能的收集率之比，是表征結冰表面法向水滴流率的無量綱參數，在歐拉法中其定義式為

式中，u和n分別為結冰表面當地的水滴速度矢量和單位法矢量；LWC 為空氣的液態水含量；U為自由流水滴速度大?。沪罭為水滴的無量綱（量綱一）體積分數。

1.3 結冰熱力學模型

基于空氣和水滴流場、水滴撞擊特性的計算結果，獲得機翼表面局部水收集系數，通過對機翼建立結冰模型并求解，得到結冰后的結冰形狀。本文采用了考慮初始水膜流動的Shallow-Water結冰熱力學模型。

Shallow-Water 結冰模型是基于表面水膜運動而建立的機翼表面結冰模型[18]。在結冰預測研究中，使用Shallow-Water結冰模型進行結冰模擬，該模型水膜流動受三維偏微分方程（PDE）控制，并在結構化和非結構化表面網格上采用穩定的有限體積格式離散。Shallow-Water 模型由質量和能量平衡的兩個偏微分方程與動量系統的代數方程組成。模型方程概述主要有以下幾個[19]。

質量守恒方程為

式中，ρw為水的密度；hf為水膜高度；u為水膜速度矢量；m?β、m?ice、m?evap為液滴撞擊、凍結、蒸發的質量通量。

動量方程為

式中，τ為表面切應力矢量；μw為水動力黏度；a為水膜受到的累積加速度（離心力、科里奧利力、重力）。

能量方程為

式中，T為取決于霜冰或光冰狀態的水膜或冰的溫度；Cw為水的比熱；Cice為冰的比熱；Td為初始撞擊的水滴的溫度；Lf為水的凍結潛熱；Q?evap、Q?h、Q?cond為分別來自蒸發、對流、傳導的熱通量。

聯立上述方程，可將其化簡為T和凍結系數n的方程，其中n可表示為

式中，?為水滴能量傳遞參數；ε為空氣能量傳遞參數；b為相對熱因子。

求解得到T和凍結系數n之后，再結合水滴撞擊求解得到的結果，積分可得整個翼面上的結冰質量和厚度。

2 空氣流場及氣動特性驗證及對比分析研究

首先，驗證翼身組合體布局飛機的空氣流場計算。選擇美國航空航天學會（AIAA）第二屆阻力預測研討會模型DLR-F4 翼身組合體作為本文傳統布局飛機計算模型，并將數值計算研究與AIAA 第二屆阻力預測會議文獻（DPW2）進行對比驗證。計算條件為：馬赫數Ma=0.6，雷諾數Re=3×106，參考弦長c=0.1412m。計算模型為對稱模型，在對稱面上采用對稱邊界條件，網格總量約為870萬個，與參考文獻[20]中量級相當。

在不同來流迎角下，數值計算得到的升力系數與試驗數據[20]對比如圖1所示；在不同來流迎角下，數值計算得到的阻力系數與試驗數據[20]對比如圖2所示。計算結果與試驗數據較相符，驗證了空氣流場計算方法的準確性。

圖1 升力特性曲線Fig.1 Lift characteristic curve

圖2 阻力特性曲線Fig.2 Resistance characteristic curve

然后驗證翼身融合飛機空氣流場。本節計算模型選擇5.8%縮比的N2A翼身融合體，并與NASA蘭利亞聲速風洞試驗結果[21]進行對比驗證。參考文獻試驗條件為：馬赫數Ma=0.2，雷諾數Re=6.6×106，參考弦長c=1.538m。

計算模型為對稱模型，采用半模計算，在對稱面上采用對稱邊界條件來減少計算量。在0°～10°范圍內不同來流迎角下，對比數值計算升力系數與試驗數據[20]的誤差。為保證計算誤差符合要求，采用S-A 湍流模型進行計算，本文N2A 翼身構型計算網格總量約為2800 萬個，與參考文獻[22]網格量級相當。

圖3 對比了不同來流迎角下，數值計算升力系數與試驗數據[21]；圖4對比了不同迎角下所對應的升力系數、阻力系數的值與試驗數據[21]，計算結果與試驗數據均對比良好，驗證了空氣流場計算方法的準確性。

圖3 升力特性曲線Fig.3 Lift characteristic curve

圖4 升阻極曲線Fig.4 Lifting resistor curve

從計算結果與試驗結果的對比可以看出，本文所采用方法的計算結果與文獻數據均對比良好，誤差均在允許范圍內，驗證了本文計算方法的準確性。

圖5 所示的是DLR-F4 飛機在來流迎角為5.384°時的壓力分布云圖。從圖5中可以看出，傳統布局飛機DLR-F4飛機機翼前后緣處上下表面壓力差較大，機身處上下表面壓力系數較為相近，可知傳統布局飛機DLR-F4 升力幾乎都由機翼提供，機身提供的升力可忽略。

圖5 DLR-F4飛機上下表面壓力分布云圖Fig.5 DLR-F4 aircraft upper and lower surface pressure distribution cloud chart

圖6所示的是N2A翼身融合布局飛機在來流迎角為5°時的壓力分布云圖。從圖6 中可以看出，與傳統布局飛機不同的是，翼身融合飛機機身、機翼和翼身融合處都是升力面。從壓力系數分布云圖中可知，翼身融合布局飛機全機都是升力面，機翼前緣和機身頭部位置上下表面的壓力系數差距都十分明顯。

圖6 N2A飛機上下表面壓力分布云圖Fig.6 N2A aircraft upper and lower surface pressure distribution cloud chart

3 結冰數值模擬過程及結果對比分析研究

3.1 DLR-F4飛機表面冰形生成及其特征

選擇DLR-F4 作為本節飛機計算模型，但DLR-F4 帶冰形飛機模型目前還未公開，采用結冰軟件FENSAP-ICE來構造DLR-F6 翼身組合體的冰形。計算條件為：大氣壓力為39000Pa，雷諾數Re=3×106。第2節通過與試驗結果對比，驗證了空氣流場計算方法，保證了計算結果的準確性。由于結冰通常都是發生在飛機低速起降和爬升等的穿云飛行過程中，對應的馬赫數較低，所以結冰計算工況選擇典型的霜冰工況，見表1。本文采用單步法對DLR-F4翼身組合體進行結冰數值模擬。

表1 結冰計算條件Table 1 Ice calculation conditions

在進行N-S 方程求解過程中，對邊界層內的流場進行模擬時，結構網格有著非結構網格不可比擬的優點，計算網格采用結構化網格，與前文計算網格一樣，網格圖如圖7所示。網格總量約為870萬個，機身網格較粗，機翼前緣處為主要結冰區域，需要加密網格更好地捕捉流場信息，以良好地描述冰形，故加密此處網格。

圖7 DLR-F4計算網格Fig.7 DLR-F4 computing grid

圖8所示為DLR-F4飛機表面生成的冰形圖，雖然結冰時間只有360s，但是在機翼前緣和機頭部位仍形成較大尺度的結冰。

圖8 DLR-F4冰形生成圖Fig.8 DLR-F4 ice formation diagram

圖9 為機翼沿展向在z/b=20%，z/b=50%，z/b=80%三個不同截面的冰形圖，其中z為展向長度，b為半展長，各截面冰形呈現流線型，且相對翼型尺寸由翼根至翼尖逐漸變大。

圖9 DLR-F4飛機不同截面冰形圖Fig.9 Ice diagram of different sections of DLR-F4 aircraft

3.2 N2A飛機表面冰形生成及其特征

選擇5.8%比例的N2A翼身融合體作為飛機計算模型。目前對N2A 翼身融合布局結冰數值模擬的公開文獻尚未出現，采用結冰軟件FENSAP-ICE來構造N2A翼身融合布局飛機的冰形。

本節結冰計算工況選擇典型的霜冰工況，見表2。采用單步法對N2A翼身組合體進行結冰數值模擬。

表2 結冰計算條件Table 2 Ice calculation conditions

計算網格采用結構化網格，與前文計算網格一樣，網格圖如圖10 所示。而新型翼身融合布局構型將機身融合成機翼的一部分，使飛機的升阻比顯著提高，但這一部分機身迎風面也變成飛機主要結冰區域，故進行機身頭部加密。

圖10 N2A計算網格Fig.10 N2A computing grid

圖11所示為N2A飛機表面結冰圖。根據計算條件，雖然結冰時間只有360s，但是在機翼前緣仍形成較大尺度的冰。與傳統翼身組合體構型由機翼和機身兩個部件結合而成，機翼前緣明顯是主要結冰區域不同，該構型飛機機體成為一個完整的升力面，升力面上的結冰區域比較模糊，從大約沿展向在z/b=10%位置處發生結冰現象。

圖11 N2A冰形生成圖Fig.11 N2A ice formation diagram

圖12 為機翼沿展向在z/b=20%，z/b=50%，z/b=80%三個不同截面的冰形圖，從圖12 中能看出，各截面冰形附著在翼型表面，為流向冰，相對翼型尺寸由翼根至翼尖逐漸變大。

3.3 DLR-F4飛機與N2A飛機結冰特性對比分析

圖13 和圖14 所示分別為DLR-F4 和N2A 飛機表面的冰生長質量流量圖。從圖13與圖14中可以看出，翼身組合體飛機的機翼前緣與機頭是結冰速率最快的區域，但N2A翼身融合布局飛機結冰速率最快的區域為機頭、翼身融合段與機翼的前緣。

圖13 DLR-F4表面冰生長質量流量Fig.13 DLR-F4 surface ice growth mass flow

圖14 N2A表面冰生長質量流量Fig.14 N2A surface ice growth mass flow

圖15是DLR-F4翼身組合體結冰后飛機上下表面的冰形圖。傳統翼身構型飛機由機翼和圓筒形機身兩個部件結合而成，飛機的機身和機翼區別明顯。如圖15 所示，飛機的迎風面發生了結冰現象，機翼作為飛機主要升力面生成的冰形最厚，機翼前緣是主要結冰區域。

圖15 DLR-F4翼身組合體結冰后上下表面的冰形圖Fig.15 Ice diagram of the front and rear surfaces of the DLR-F4 wing-body combination after freezing

圖16 是N2A 翼身融合體結冰后前后表面的冰形圖。新型翼身融合布局構型為了減小翼身之間的干擾阻力和誘導阻力，減小飛機的總阻力，將機翼與翼身融合，使得整個飛機機體成為一個大的升力面。如圖16所示，飛機的迎風面發生了結冰現象，由于整個飛機機體為飛機提供升力，翼身融合體整機前緣是主要結冰區域。

圖16 N2A翼身融合體結冰后上下表面的冰形圖Fig.16 Ice diagram of the front and rear surfaces of the N2A blended-wing-body after freezing

從計算結果中可以看出，相較于DLR-F4翼身組合體，N2A 翼身融合體在機身、翼身融合段和機翼前緣均發生結冰現象。在后續的結冰研究和防/除冰研究中，需要考慮翼身融合布局飛機全翼面上的結冰問題。

4 結論

本文針對兩種布局飛機結冰，首先構建了結冰數值模擬方法，驗證了兩種布局飛機的空氣流場，之后通過數值計算方法針對新型翼身融合布局構型和傳統布局飛機進行霜冰結冰數值模擬，對比分析其結冰特性，并得出以下結論：

（1）驗證了DLR-F4與N2A空氣流場計算方法，并對比分析氣動特性，得出傳統布局飛機升力由機翼提供，翼身融合布局飛機升力由機翼和機身共同提供，為后續結冰特性分析研究奠定良好基礎。

（2）對翼身融合布局構型N2A在典型霜冰結冰條件下開展結冰數值模擬，分析其結冰區域及表面冰生長質量流量，并將N2A 結冰特性與傳統翼身組合體DLR-F4 結冰特性進行比較，發現需要考慮機翼與機身融合處的結冰及防/除冰問題，翼身融合布局飛機機翼部分的結冰與傳統布局飛機較為相似。