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

西北工業(yè)大學(xué)，陜西 西安 710072

結(jié)冰是危害飛機(jī)飛行安全的重要因素，由結(jié)冰導(dǎo)致的飛行事故常常發(fā)生。機(jī)翼和發(fā)動(dòng)機(jī)[1]等部位的結(jié)冰問題會(huì)對飛機(jī)的氣動(dòng)特性產(chǎn)生劇烈影響，從而降低飛行質(zhì)量，情況嚴(yán)重時(shí)甚至?xí)?dǎo)致機(jī)毀人亡等危害，因而受到人們的廣泛關(guān)注。目前國內(nèi)外的研究人員致力于研究過冷大水滴和冰晶結(jié)冰[2]，為防/除冰方案設(shè)計(jì)提供了可靠依據(jù)。

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

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

1 結(jié)冰數(shù)值模擬計(jì)算方法

1.1 空氣流場數(shù)值求解

空氣流場計(jì)算是研究結(jié)冰數(shù)值計(jì)算的基礎(chǔ)環(huán)節(jié)，也是研究水滴撞擊特性的關(guān)鍵。對結(jié)冰空氣流場的計(jì)算基于對計(jì)算域中Navier-Stokes（N-S）方程組的求解。空氣流場在進(jìn)行流動(dòng)和傳熱時(shí)，可用控制方程對質(zhì)量守恒定律、動(dòng)量定律和能量守恒定律三個(gè)方程進(jìn)行描述。N-S 方程積分形式為

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

1.2 水滴撞擊特性計(jì)算

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

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

歐拉法求解水滴運(yùn)動(dòng)特性的連續(xù)方程主要有以下幾個(gè)。

空氣連續(xù)方程為

水滴連續(xù)方程為

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

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

式中，u和n分別為結(jié)冰表面當(dāng)?shù)氐乃嗡俣仁噶亢蛦挝环ㄊ噶浚籐WC 為空氣的液態(tài)水含量；U為自由流水滴速度大小；αN為水滴的無量綱（量綱一）體積分?jǐn)?shù)。

1.3 結(jié)冰熱力學(xué)模型

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

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

質(zhì)量守恒方程為

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

動(dòng)量方程為

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

能量方程為

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

聯(lián)立上述方程，可將其化簡為T和凍結(jié)系數(shù)n的方程，其中n可表示為

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

求解得到T和凍結(jié)系數(shù)n之后，再結(jié)合水滴撞擊求解得到的結(jié)果，積分可得整個(gè)翼面上的結(jié)冰質(zhì)量和厚度。

2 空氣流場及氣動(dòng)特性驗(yàn)證及對比分析研究

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

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

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

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

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

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

圖3 對比了不同來流迎角下，數(shù)值計(jì)算升力系數(shù)與試驗(yàn)數(shù)據(jù)[21]；圖4對比了不同迎角下所對應(yīng)的升力系數(shù)、阻力系數(shù)的值與試驗(yàn)數(shù)據(jù)[21]，計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)均對比良好，驗(yàn)證了空氣流場計(jì)算方法的準(zhǔn)確性。

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

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

從計(jì)算結(jié)果與試驗(yàn)結(jié)果的對比可以看出，本文所采用方法的計(jì)算結(jié)果與文獻(xiàn)數(shù)據(jù)均對比良好，誤差均在允許范圍內(nèi)，驗(yàn)證了本文計(jì)算方法的準(zhǔn)確性。

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

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

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

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

3 結(jié)冰數(shù)值模擬過程及結(jié)果對比分析研究

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

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

表1 結(jié)冰計(jì)算條件Table 1 Ice calculation conditions

在進(jìn)行N-S 方程求解過程中，對邊界層內(nèi)的流場進(jìn)行模擬時(shí)，結(jié)構(gòu)網(wǎng)格有著非結(jié)構(gòu)網(wǎng)格不可比擬的優(yōu)點(diǎn)，計(jì)算網(wǎng)格采用結(jié)構(gòu)化網(wǎng)格，與前文計(jì)算網(wǎng)格一樣，網(wǎng)格圖如圖7所示。網(wǎng)格總量約為870萬個(gè)，機(jī)身網(wǎng)格較粗，機(jī)翼前緣處為主要結(jié)冰區(qū)域，需要加密網(wǎng)格更好地捕捉流場信息，以良好地描述冰形，故加密此處網(wǎng)格。

圖7 DLR-F4計(jì)算網(wǎng)格Fig.7 DLR-F4 computing grid

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

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

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

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

3.2 N2A飛機(jī)表面冰形生成及其特征

選擇5.8%比例的N2A翼身融合體作為飛機(jī)計(jì)算模型。目前對N2A 翼身融合布局結(jié)冰數(shù)值模擬的公開文獻(xiàn)尚未出現(xiàn)，采用結(jié)冰軟件FENSAP-ICE來構(gòu)造N2A翼身融合布局飛機(jī)的冰形。

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

表2 結(jié)冰計(jì)算條件Table 2 Ice calculation conditions

計(jì)算網(wǎng)格采用結(jié)構(gòu)化網(wǎng)格，與前文計(jì)算網(wǎng)格一樣，網(wǎng)格圖如圖10 所示。而新型翼身融合布局構(gòu)型將機(jī)身融合成機(jī)翼的一部分，使飛機(jī)的升阻比顯著提高，但這一部分機(jī)身迎風(fēng)面也變成飛機(jī)主要結(jié)冰區(qū)域，故進(jìn)行機(jī)身頭部加密。

圖10 N2A計(jì)算網(wǎng)格Fig.10 N2A computing grid

圖11所示為N2A飛機(jī)表面結(jié)冰圖。根據(jù)計(jì)算條件，雖然結(jié)冰時(shí)間只有360s，但是在機(jī)翼前緣仍形成較大尺度的冰。與傳統(tǒng)翼身組合體構(gòu)型由機(jī)翼和機(jī)身兩個(gè)部件結(jié)合而成，機(jī)翼前緣明顯是主要結(jié)冰區(qū)域不同，該構(gòu)型飛機(jī)機(jī)體成為一個(gè)完整的升力面，升力面上的結(jié)冰區(qū)域比較模糊，從大約沿展向在z/b=10%位置處發(fā)生結(jié)冰現(xiàn)象。

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

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

3.3 DLR-F4飛機(jī)與N2A飛機(jī)結(jié)冰特性對比分析

圖13 和圖14 所示分別為DLR-F4 和N2A 飛機(jī)表面的冰生長質(zhì)量流量圖。從圖13與圖14中可以看出，翼身組合體飛機(jī)的機(jī)翼前緣與機(jī)頭是結(jié)冰速率最快的區(qū)域，但N2A翼身融合布局飛機(jī)結(jié)冰速率最快的區(qū)域?yàn)闄C(jī)頭、翼身融合段與機(jī)翼的前緣。

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

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

圖15是DLR-F4翼身組合體結(jié)冰后飛機(jī)上下表面的冰形圖。傳統(tǒng)翼身構(gòu)型飛機(jī)由機(jī)翼和圓筒形機(jī)身兩個(gè)部件結(jié)合而成，飛機(jī)的機(jī)身和機(jī)翼區(qū)別明顯。如圖15 所示，飛機(jī)的迎風(fēng)面發(fā)生了結(jié)冰現(xiàn)象，機(jī)翼作為飛機(jī)主要升力面生成的冰形最厚，機(jī)翼前緣是主要結(jié)冰區(qū)域。

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

圖16 是N2A 翼身融合體結(jié)冰后前后表面的冰形圖。新型翼身融合布局構(gòu)型為了減小翼身之間的干擾阻力和誘導(dǎo)阻力，減小飛機(jī)的總阻力，將機(jī)翼與翼身融合，使得整個(gè)飛機(jī)機(jī)體成為一個(gè)大的升力面。如圖16所示，飛機(jī)的迎風(fēng)面發(fā)生了結(jié)冰現(xiàn)象，由于整個(gè)飛機(jī)機(jī)體為飛機(jī)提供升力，翼身融合體整機(jī)前緣是主要結(jié)冰區(qū)域。

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

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

4 結(jié)論

本文針對兩種布局飛機(jī)結(jié)冰，首先構(gòu)建了結(jié)冰數(shù)值模擬方法，驗(yàn)證了兩種布局飛機(jī)的空氣流場，之后通過數(shù)值計(jì)算方法針對新型翼身融合布局構(gòu)型和傳統(tǒng)布局飛機(jī)進(jìn)行霜冰結(jié)冰數(shù)值模擬，對比分析其結(jié)冰特性，并得出以下結(jié)論：

（1）驗(yàn)證了DLR-F4與N2A空氣流場計(jì)算方法，并對比分析氣動(dòng)特性，得出傳統(tǒng)布局飛機(jī)升力由機(jī)翼提供，翼身融合布局飛機(jī)升力由機(jī)翼和機(jī)身共同提供，為后續(xù)結(jié)冰特性分析研究奠定良好基礎(chǔ)。

（2）對翼身融合布局構(gòu)型N2A在典型霜冰結(jié)冰條件下開展結(jié)冰數(shù)值模擬，分析其結(jié)冰區(qū)域及表面冰生長質(zhì)量流量，并將N2A 結(jié)冰特性與傳統(tǒng)翼身組合體DLR-F4 結(jié)冰特性進(jìn)行比較，發(fā)現(xiàn)需要考慮機(jī)翼與機(jī)身融合處的結(jié)冰及防/除冰問題，翼身融合布局飛機(jī)機(jī)翼部分的結(jié)冰與傳統(tǒng)布局飛機(jī)較為相似。