一種基于水滴收集量代理模型的結(jié)冰試飛空域確定方法

牛俊杰，桑為民，2，*，李棟，郝蓮，王澤林

1．西北工業(yè)大學(xué) 航空學(xué)院，西安 710072

2．中國(guó)空氣動(dòng)力研究與發(fā)展中心 結(jié)冰與防除冰重點(diǎn)實(shí)驗(yàn)室，綿陽(yáng) 621000

3．中國(guó)商用飛機(jī)有限責(zé)任公司 上海飛機(jī)設(shè)計(jì)研究院，上海 201210

4．中國(guó)飛行試驗(yàn)研究院，西安 710089

飛機(jī)在含有過冷水滴的云層中飛行時(shí)，會(huì)在其表面產(chǎn)生結(jié)冰。飛機(jī)結(jié)冰會(huì)嚴(yán)重威脅飛行安全［1-2］。據(jù) 統(tǒng) 計(jì)，僅2003—2008年 間 就有380余起與結(jié)冰相關(guān)的飛行事故［3］。中國(guó)自2006年6·3空難以來，飛機(jī)結(jié)冰以及由此引起的飛行事故越來越受到關(guān)注，特別是大飛機(jī)項(xiàng)目的啟動(dòng)，使得與飛機(jī)結(jié)冰相關(guān)的研究成為研究熱點(diǎn)［4］。適航規(guī)章規(guī)定飛機(jī)必須經(jīng)過結(jié)冰適航審定才能夠容許在結(jié)冰環(huán)境下飛行［5-6］。飛行試驗(yàn)驗(yàn)證的目的是驗(yàn)證飛機(jī)在結(jié)冰條件下飛機(jī)的飛行性能、飛行品質(zhì)能夠滿足適航審定要求，防除冰系統(tǒng)能夠履行其功能。與美國(guó)、歐洲等航空事業(yè)較發(fā)達(dá)的國(guó)家相比，中國(guó)飛機(jī)結(jié)冰適航驗(yàn)證工作起步較晚。中國(guó)雖然完成了Y12-II、Y7-200A、ARJ21-700等飛機(jī)結(jié)冰適航驗(yàn)證，但在適航驗(yàn)證方法和許多關(guān)鍵技術(shù)方面與國(guó)外差距較大，尚未形成系統(tǒng)的飛機(jī)結(jié)冰適航驗(yàn)證和試飛技術(shù)［7］。通常飛行試驗(yàn)驗(yàn)證分為干空氣飛行試驗(yàn)、帶模擬冰型的飛行試驗(yàn)和自然結(jié)冰飛行試驗(yàn)3類。自然結(jié)冰試飛試驗(yàn)是在自然結(jié)冰狀態(tài)下，驗(yàn)證飛機(jī)在暴露于適航條例規(guī)定的連續(xù)和間斷最大結(jié)冰狀態(tài)下能夠安全運(yùn)行，其系統(tǒng)能夠達(dá)到預(yù)定的功能和性能。自然結(jié)冰試飛試驗(yàn)的關(guān)鍵是尋找適合的結(jié)冰氣象條件，在一定的結(jié)冰氣象條件下進(jìn)行試飛試驗(yàn)，驗(yàn)證飛機(jī)的各項(xiàng)功能滿足適航要求。中國(guó)基于自然科學(xué)研究的結(jié)冰探測(cè)試驗(yàn)過程“屈指可數(shù)”，完整的結(jié)冰探測(cè)試驗(yàn)過程的具體數(shù)量尚沒有權(quán)威的統(tǒng)計(jì)，主要有：1986—1987年Y-12II在新疆完成了5次結(jié)冰試飛［7］；1997年Y7-200A飛機(jī)自然結(jié)冰飛行試驗(yàn)選擇新疆天山地區(qū)為試驗(yàn)區(qū)，在滿足飛機(jī)自然結(jié)冰條件時(shí)進(jìn)行了2次自然結(jié)冰試飛工作［8］；2018年中國(guó)生產(chǎn)的某型直升機(jī)在新疆五家渠自然結(jié)冰試飛［9］以及從2011年開始ARJ21結(jié)冰試飛［10］。ARJ21首先在新疆進(jìn)行結(jié)冰試飛，最終在加拿大完成結(jié)冰試飛全部科目［10］。從中可以發(fā)現(xiàn)中國(guó)之前的結(jié)冰試飛主要在新疆以及加拿大完成。ARJ21前期連續(xù)3年在新疆找結(jié)冰環(huán)境，真正飛成功的架次較少，之后去了北美五大湖1.5個(gè)月完成結(jié)冰試飛所有科目。從實(shí)際業(yè)務(wù)發(fā)展需求角度出發(fā)，自然結(jié)冰試飛氣象預(yù)測(cè)已成為嚴(yán)重制約國(guó)產(chǎn)大飛機(jī)的研制和交付使用的一個(gè)重要方面［10］。

在自然結(jié)冰試飛過程中，常常使用結(jié)冰指數(shù)來表征可能結(jié)冰的區(qū)域以及結(jié)冰強(qiáng)度，如：基于Schultz和Politovich研 究 成 果 的Icing Indices［11］，該指數(shù)基于溫度和相對(duì)濕度對(duì)結(jié)冰概率進(jìn)行預(yù)測(cè)，結(jié)冰概率為溫度與濕度的函數(shù)；由Belo-Pereira提出一種新的結(jié)冰指數(shù)SFIP（Simplified Forecast Icing Potential），該指數(shù)基于數(shù)值氣象預(yù)測(cè)的結(jié)果，認(rèn)為結(jié)冰指數(shù)是溫度、相對(duì)濕度以及液態(tài) 水 含 量（Liquid Water Content， LWC）的 函數(shù)［12］；Politovich和Sand提出的一種結(jié)冰嚴(yán)重指數(shù)［13］，該指數(shù)考慮影響結(jié)冰的3個(gè)主要參數(shù)：大氣環(huán)境溫度、云中過冷水含量及云滴的中位數(shù)體積直徑；Thompson等提出結(jié)冰方案，該方案根據(jù)每個(gè)探空層上的溫度（T）、露點(diǎn)（Td）以及在本層與上層之間的溫度遞減率，確定該層的結(jié)冰強(qiáng)度和類型，結(jié)冰強(qiáng)度劃分為8個(gè)等級(jí)［14］；Thompson等提出的RAP（Research Applications Program）算法根據(jù)溫度、相對(duì)濕度定義了4種結(jié)冰類型：層狀、凍雨、不穩(wěn)定、普通［15］；NCAR （National Center for Atmospheric Research）提出模糊邏輯預(yù)報(bào)方法，該方法是將和結(jié)冰有關(guān)的溫度、濕度、云量、云水等與結(jié)冰的可能性相聯(lián)系，做成曲線，找出各種情況下的結(jié)冰可能性和結(jié)冰強(qiáng)度，該方法又被稱為當(dāng)前結(jié)冰潛勢(shì)分析（Current Icing Potential，CIP）［16］。李佰平等［10］對(duì)自然結(jié)冰潛勢(shì)算法進(jìn)行了改進(jìn)，該算法直接基于大氣溫濕層結(jié)給出結(jié)冰潛勢(shì)。

然而使用結(jié)冰指數(shù)大多只能給出結(jié)冰的可能性或者對(duì)結(jié)冰強(qiáng)度進(jìn)行預(yù)測(cè)。飛機(jī)結(jié)冰受到多種因素的影響，包括結(jié)冰氣象因素（如溫度、液態(tài)水含量、水滴中值體積直徑以及高度等）和飛行狀態(tài)因素（如速度、前緣直徑、攻角以及結(jié)冰時(shí)間等）。一般而言，溫度主要影響冰的形態(tài)，液態(tài)水含量以及速度影響結(jié)冰量，水滴中值體積直徑影響水滴的撞擊極限位置，高度影響溫度。然而，溫度會(huì)對(duì)液態(tài)水含量產(chǎn)生影響，溫度越低，空氣中的液態(tài)水含量越低。這是由于空氣中的過冷液態(tài)水會(huì)隨著溫度的降低逐漸向冰晶轉(zhuǎn)化，進(jìn)而導(dǎo)致液態(tài)水含量的降低。當(dāng)溫度低于?40 ℃時(shí)，被認(rèn)為很少有過冷液態(tài)水存在［5］。結(jié)冰氣象因素屬于不可控因素，需要特定的天氣條件，而飛行狀態(tài)因素中速度和時(shí)間是可控因素。因此在一定的結(jié)冰氣象條件下，可以通過調(diào)節(jié)飛行速度或者時(shí)間達(dá)到使結(jié)冰氣象條件達(dá)到滿足結(jié)冰試飛的條件。

因此，提出一種基于水滴收集量代理模型結(jié)合WRF（Weather Research and Forecasting Model）氣象預(yù)測(cè)的結(jié)冰試飛空域確定方法。首先對(duì)美國(guó)聯(lián)邦航空條例（FAR）25部附錄C進(jìn)行采樣，對(duì)采樣點(diǎn)進(jìn)行空氣流場(chǎng)和水滴流場(chǎng)數(shù)值模擬，獲得不同采樣點(diǎn)下翼型的水滴收集量，利用本征正交分解（Proper Orthogonal Decomposition，POD）方法及Kriging插值建立氣象因素（溫度、高度、液態(tài)水含量、水滴中值體積直徑）、飛行狀態(tài)（速度、高度）與水滴收集量之間的代理模型。基于WRF獲得的氣象條件，結(jié)合水滴收集量代理模型，則可以獲得目標(biāo)區(qū)域內(nèi)水滴收集量的空間分布及隨時(shí)間的變化，以中度結(jié)冰強(qiáng)度對(duì)目標(biāo)區(qū)域的水滴收集量進(jìn)行劃分，進(jìn)而獲得可用于結(jié)冰試飛的空域并且對(duì)結(jié)冰強(qiáng)度進(jìn)行預(yù)測(cè)。

1 計(jì)算方法及控制方程

基于水滴收集量代理模型及WRF氣象預(yù)測(cè)的結(jié)冰試飛空域確定包括2部分內(nèi)容：一部分為水滴收集量代理模型構(gòu)建，另一部分為WRF氣象條件預(yù)測(cè)。在水滴收集量代理模型構(gòu)建中，采用優(yōu)化拉丁超立方采樣對(duì)附錄C進(jìn)行采樣，獲得連續(xù)最大結(jié)冰狀態(tài)下的40個(gè)采樣點(diǎn)，對(duì)采樣點(diǎn)進(jìn)行空氣流場(chǎng)和水滴撞擊特性計(jì)算，獲得不同采樣點(diǎn)工況下的水滴收集量。采用POD方法和Krig?ing插值構(gòu)建水滴收集量的代理模型。采用WRF進(jìn)行氣象預(yù)測(cè)，獲得目標(biāo)區(qū)域內(nèi)的氣象參數(shù)，然后基于氣象參數(shù)和水滴收集量代理模型，對(duì)目標(biāo)區(qū)域進(jìn)行水滴收集量快速預(yù)測(cè)，獲得目標(biāo)區(qū)域內(nèi)水滴收集量的空間分布及隨時(shí)間的變化，以中度結(jié)冰強(qiáng)度對(duì)目標(biāo)區(qū)域的水滴收集量進(jìn)行劃分，進(jìn)而獲得可用于結(jié)冰試飛的空域并且對(duì)結(jié)冰強(qiáng)度進(jìn)行預(yù)測(cè)。結(jié)冰試飛區(qū)域確定流程如圖1所示。在本節(jié)中對(duì)水滴收集量的代理模型構(gòu)建進(jìn)行說明，WRF氣象預(yù)測(cè)的說明見2.2節(jié)。

圖1 結(jié)冰試飛區(qū)域確定流程Fig. 1 Flow chart of icing test flight area determination

1.1 水滴撞擊特性求解

結(jié)冰數(shù)值模擬過程通常包括空氣流場(chǎng)求解、水滴撞擊特性求解以及結(jié)冰傳熱傳質(zhì)模擬3部分內(nèi)容。本文旨在求解壁面處水滴收集量分布，因此只針對(duì)上述計(jì)算過程中前2個(gè)步驟進(jìn)行數(shù)值模擬，即空氣流場(chǎng)的求解和水滴撞擊特性求解。在空氣流場(chǎng)求解的基礎(chǔ)上，采用拉格朗日方法建立過冷水滴運(yùn)動(dòng)方程，并進(jìn)行求解，獲得過冷水滴的撞擊特性及壁面水滴收集量分布。

空氣流場(chǎng)通過求解Navier-Stokes方程組獲得，方程組的張量形式為

式中：ρ為空氣密度；t為時(shí)間；u為空氣速度，下標(biāo)i、j為張量指標(biāo)；x表示方向；p為壓力；Sij為應(yīng)力張量；μ為空氣動(dòng)力黏性系數(shù)；cp為空氣定壓比熱容；T為溫度；κ為空氣導(dǎo)熱系數(shù)分別為湍流應(yīng)力張量和湍流熱流通量。湍流模型采 用SST（Shear Stress Transport）k-ω湍 流 模型［17］，空氣流場(chǎng)求解使用SIMPLE算法，壓力、動(dòng)量、湍流動(dòng)能和能量項(xiàng)以二階方式離散。壓力和動(dòng)量松弛因子為0.4，收斂精度設(shè)定為10?6。在獲得空氣流場(chǎng)分布基礎(chǔ)上，采用拉格朗日法獲得過冷水滴的運(yùn)動(dòng)軌跡。拉格朗日法對(duì)水滴軌跡進(jìn)行求解時(shí)，對(duì)每個(gè)水滴進(jìn)行跟蹤，其中假設(shè)［18-19］：水滴在運(yùn)動(dòng)過程中既不相互碰撞也不分解，水滴的密度、溫度等在運(yùn)動(dòng)中保持不變；水滴的初速度與自由來流相同，水滴流場(chǎng)不會(huì)對(duì)空氣流場(chǎng)產(chǎn)生影響。水滴的運(yùn)動(dòng)軌跡及撞擊位置如圖2所示。

圖2 水滴軌跡示意圖Fig. 2 Schematic diagram of water droplet trajectory

水滴的運(yùn)動(dòng)過程可以表示為

式中：uw為水滴速度；ua為空氣速度；ρa(bǔ)為空氣密度；ρw為水滴密度；g為重力加速度；Cd為水滴曳力系數(shù)；deq為水滴直徑；μa為空氣黏性系數(shù)；Rew為相對(duì)雷諾數(shù)，定義為

在計(jì)算得到空氣流場(chǎng)的結(jié)果后，采用四階龍格-庫(kù)塔法對(duì)式（4）進(jìn)行求解，獲得了水滴的運(yùn)動(dòng)軌跡及撞擊位置。由水滴的初始位置和碰撞位置，得到局部水滴收集系數(shù)β：

式中：dy為初始位置時(shí)相鄰水滴的距離；ds為相鄰水滴翼面碰撞位置的距離，水滴位置如圖2所示。壁面水滴收集量為

其中：LWC為液態(tài)水含量；V∞為遠(yuǎn)場(chǎng)速度；β為水滴收集系數(shù)。

1.2 POD方法

本征正交分解（Proper Orthogonal Decom?position，POD）是一種數(shù)值方法，可以降低計(jì)算流體力學(xué)和結(jié)構(gòu)分析等計(jì)算模擬的復(fù)雜性。通常在流體動(dòng)力學(xué)和湍流分析中，其用更簡(jiǎn)單的模型來代替Navier-Stokes方程［20-21］。假設(shè)存在這樣的物理場(chǎng)y=y（x， tn），該物理場(chǎng)可以表示為無(wú)窮級(jí)數(shù)的形式：

式中：?i(x)為特征基函數(shù)；αk(tk)為經(jīng)驗(yàn)系數(shù)；tn為參數(shù)變量。利用相關(guān)的數(shù)學(xué)理論可將式（8）轉(zhuǎn)換為求解特征值的問題：

其中： 為內(nèi)積符號(hào)；λ為特征值。在實(shí)際使用中，可應(yīng)用快照方法將式（9）轉(zhuǎn)換為一個(gè)與快照個(gè)數(shù)N相等的特征值問題。若一個(gè)樣本包含M個(gè)數(shù)據(jù)，則取N個(gè)樣本組成一個(gè)樣本集U={U1，U2，…，UN}，特征基函數(shù)可以由快照的線性組合構(gòu)成，即

式中：λn為特征值。A中各元素為

其中：矩陣A的特征值大小代表該正交基函數(shù)在整個(gè)樣本空間中所占能量的比重。因此進(jìn)行物理場(chǎng)重建時(shí)便可以根據(jù)能量的大小對(duì)基函數(shù)進(jìn)行篩選，而并不需要選擇所有的基函數(shù)。假設(shè)選用了M個(gè)正交基函數(shù)，則有M≤N。在進(jìn)行矩陣分解得到最佳正交基后，對(duì)選取的正交基與樣本進(jìn)行內(nèi)積即可得到對(duì)應(yīng)的系數(shù)：

因此，重構(gòu)的物理場(chǎng)可以使用正交基的線性疊加表示為

1.3 Kriging插值

Kriging插值方法是20世紀(jì)50年代由南非工程師Krige所提出的［22］，已廣泛應(yīng)用于各個(gè)領(lǐng)域。在飛機(jī)結(jié)冰和防/除冰領(lǐng)域，Pellissier等［23］基于POD、Kriging方法以及遺傳算法，以噴射角、噴射間距及前緣距離為控制變量對(duì)防冰腔布局進(jìn)行了優(yōu)化設(shè)計(jì)；劉藤等［24］建立了基于POD方法和Kriging模型的冰形快速預(yù)測(cè)方法，從而實(shí)現(xiàn)了多參數(shù)的冰形快速預(yù)測(cè)。

假定N個(gè)快照的系數(shù)α定義了Rk中N個(gè)超平面α（x），其中k為自由參數(shù)的數(shù)目，x為設(shè)計(jì)變量，x=[x1，xx，…，xk]T，尋找α（x）的問題也就轉(zhuǎn)變?yōu)槌矫嬷械亩嗑S插值問題。Kriging插值具體原理為

Kriging插值將未知函數(shù)視為一個(gè)高斯隨機(jī)過程，該隨機(jī)過程可以表示為

式中：F(β，x)為全局趨勢(shì)模型；H(x)為隨機(jī)分布的誤差；x為參數(shù)變量。

使用協(xié)方差來表示不同相關(guān)點(diǎn)處隨機(jī)變量的相關(guān)性，即

式中：σ2為方差；Z(x，x')為相關(guān)函數(shù)，只與空間距離有關(guān)。Kriging模型尋找最優(yōu)加權(quán)系數(shù)ω，可使得均方差最小，并且滿足無(wú)偏差條件。通常使用式（18）獲得Kriging模型的預(yù)估值：

其 中 ：Vkrig只 與 樣 本 點(diǎn) 有 關(guān) ；β0=(FTZ?1F)?1FTZ?1yk，yk為樣本數(shù)據(jù)，F(xiàn)為一維單位向量，Z和z分別為前述相關(guān)函數(shù)矩陣及其中的相關(guān)向量，因此可以利用樣本數(shù)據(jù)得到模型的估計(jì)值。關(guān)于POD和Kriging的詳細(xì)介紹可以參考文獻(xiàn)［24］。

2 數(shù)值模擬驗(yàn)證

對(duì)NACA0012翼型進(jìn)行數(shù)值模擬，獲得壁面水滴收集系數(shù)分布，與文獻(xiàn)［25］中結(jié)果進(jìn)行對(duì)比，以驗(yàn)證水滴撞擊特性數(shù)值模擬的準(zhǔn)確性。NACA 0012水滴撞擊特性驗(yàn)證工況如表1所示。

表1 水滴撞擊特性驗(yàn)證算例工況Table 1 Conditions of water droplets impingement characteristics validation cases

2.1 NACA0012水滴撞擊特性驗(yàn)證

圖3和圖4分別為NACA0012翼型周圍網(wǎng)格示意圖和水滴軌跡示意圖。圖5為NACA0012翼型表面水滴收集系數(shù)分布對(duì)比，其中實(shí)線為本文模擬結(jié)果，三角和黑色點(diǎn)劃線為文獻(xiàn)［25］中結(jié)果。橫坐標(biāo)（S?Ssp）/C表示壁面點(diǎn)與滯止點(diǎn)之間的距離與弦長(zhǎng)的比值。圖中可以看出本文模擬結(jié)果與文獻(xiàn)結(jié)果符合良好。由于來流攻角為4°，水滴碰撞到上翼面范圍較小、下翼面范圍較大，相應(yīng)的水滴收集系數(shù)在駐點(diǎn)附近達(dá)到最大值，下翼面處水滴收集系數(shù)變化相較于上翼面更為緩慢，水滴撞擊極限也相對(duì)較大。

圖3 NACA0012網(wǎng)格示意圖Fig. 3 Schematic diagram of NACA0012 mesh

圖4 NACA0012水滴軌跡示意圖Fig. 4 Schematic diagram of droplet trajectory around NACA0012

圖5 NACA0012表面水滴收集系數(shù)分布Fig. 5 Distribution of droplet collection coefficient on NACA0012 surface

2.2 WRF氣象模擬驗(yàn)證

WRF是中尺度數(shù)值天氣預(yù)報(bào)系統(tǒng)，適用于高分辨率大氣模擬。WRF模型由許多不同的物理方案組成，包括對(duì)微物理、積云參數(shù)化、表面物理、行星邊界層以及大氣輻射的不同描述。各種物理方案的介紹可以參考文獻(xiàn)［26］。在本文研究中微物理方案采用Morrison等［27］方案，該方案適用于液滴、云冰、雨、雪的數(shù)值模擬。積云參數(shù)化 方 案 采 用Kain-Fritsch方 案［28］。行 星 邊 界 層（Planet Boundary Layer， PBL）方案采用延世大學(xué)（YSU）PBL方案［29］。氣象預(yù)測(cè)的初始邊界條件采用美國(guó)國(guó)家環(huán)境預(yù)測(cè)中心（National Centers for Environmental Prediction， NCEP）的FNL再分析數(shù)據(jù)（數(shù)據(jù)精度為1°×1°）。

WRF可以獲得諸如溫度、壓力、液態(tài)水含量等參數(shù)。在WRF中可以直接輸出溫度、壓力以及云水混合比等參數(shù)，液態(tài)水含量可以通過云水混合比（Qcloud）獲得，具體方式為

在本文中，假設(shè)當(dāng)溫度低于273.15 K時(shí)，則認(rèn)為液態(tài)水為過冷態(tài)，并且液態(tài)水溫度與空氣溫度一致。云層中的液滴在與空氣進(jìn)行長(zhǎng)時(shí)間的熱量交換后，可以認(rèn)為液滴與空氣處于平衡態(tài)，液滴溫度與空氣溫度一致。Cober等［30］指出，CASP （Canadian Atlantic Storms Program） II曾開展了31架次共119 h的結(jié)冰飛行試驗(yàn)，對(duì)冷鋒、暖鋒、低壓區(qū)域以及低空層云中的中值體積直徑MVD （Median Volume Diameter） 進(jìn) 行 統(tǒng) 計(jì)，發(fā)現(xiàn)冷鋒、暖鋒、低壓區(qū)域云中的MVD平均值為20 μm，低空層云中的MVD平均值為16 μm。因此在本文中MVD則統(tǒng)一假定為較大值20 μm。

由于結(jié)冰試飛氣象條件比較苛刻且文獻(xiàn)中可用的數(shù)據(jù)較少，因此本文選用2013年4月1日華盛頓山（44.267N，71.3W，海拔1 917 m）的結(jié)冰事件進(jìn)行數(shù)值模擬，獲得的溫度和液態(tài)水含量LWC的結(jié)果與實(shí)際觀測(cè)數(shù)據(jù)和文獻(xiàn)［31］中的結(jié)果進(jìn)行比較。采用3層區(qū)域嵌套，內(nèi)外層網(wǎng)格分辨率分別為9、3、1 km，最內(nèi)層區(qū)域內(nèi)采用225×225×26個(gè)網(wǎng)格點(diǎn)進(jìn)行數(shù)值模擬。圖6為嵌套區(qū)域的示意圖。

圖7為華盛頓山溫度隨時(shí)間的變化與文獻(xiàn)［31］結(jié)果的對(duì)比。其中黑色原點(diǎn)為觀察值，實(shí)線為本文模擬結(jié)果，虛線為文獻(xiàn)［31］中的模擬結(jié)果。

圖8為華盛頓山LWC隨時(shí)間變化與文獻(xiàn)［31］結(jié)果的對(duì)比。與圖7類似，黑色原點(diǎn)為觀察值，實(shí)線是本文模擬結(jié)果，虛線是文獻(xiàn)［31］中的模擬結(jié)果。

圖6 嵌套區(qū)域示意圖Fig. 6 Schematic diagram of nested domains

圖7 溫度隨時(shí)間變化對(duì)比Fig. 7 Comparison of temperature variation with time

圖8 LWC隨時(shí)間變化對(duì)比Fig. 8 Comparison of LWC variation with time

分別以平均絕對(duì)誤差（Mean Absolute Er?ror， MAE）、均 方 根 誤 差（Root Mean Squared Error， RMSE）、平均絕對(duì)百分比誤差（Mean Ab?solute Percentage Error， MAPE）對(duì)預(yù)測(cè)結(jié)果與試驗(yàn)結(jié)果以及文獻(xiàn)［31］中的結(jié)果進(jìn)行比較，其具體定義見式（20）~式（22）。

式中：ymeas和ypred分別為結(jié)冰氣象參數(shù)觀測(cè)值與模擬結(jié)果。

表2為結(jié)冰氣象參數(shù)溫度、液態(tài)水含量預(yù)測(cè)結(jié)果與觀測(cè)結(jié)果以及文獻(xiàn)［31］預(yù)測(cè)結(jié)果與觀測(cè)結(jié)果的誤差對(duì)比?？梢园l(fā)現(xiàn)無(wú)論溫度還是液態(tài)水含量，本文預(yù)測(cè)結(jié)果均要好于文獻(xiàn)中結(jié)果，平均絕對(duì)誤差、均方根以及平均絕對(duì)百分比誤差均較小。溫度T的平均絕對(duì)誤差為1.45 K，平均絕對(duì)百分比誤差為0.5%；液態(tài)水含量的平均絕對(duì)誤差為0.07 g/m3，平均絕對(duì)百分比誤差為25.6%，好于文獻(xiàn)［31］中的結(jié)果28.7%。

表2 結(jié)冰氣象預(yù)測(cè)誤差分析Table 2 Error analysis of icing weather prediction

3 結(jié)果與討論

3.1 水滴收集量代理模型

通過對(duì)FAR 25部附錄C中連續(xù)最大結(jié)冰條件進(jìn)行優(yōu)化拉丁超立方采樣，獲得40個(gè)工況點(diǎn)，采樣點(diǎn)分布如圖9所示。以NACA0012翼型為模型，分別對(duì)40個(gè)采樣點(diǎn)進(jìn)行數(shù)值模擬，獲得不同工況下的水滴收集量分布，其中部分采樣點(diǎn)工況如表3所示。NACA0012網(wǎng)格如圖3所示。

圖10為表3中工況下水滴收集量分布云圖。圖11為表3中采樣點(diǎn)水滴收集量在翼面分布的對(duì)比。其中橫坐標(biāo)Y/C表示翼面位置Y坐標(biāo)與弦長(zhǎng)C的比值，其中正值表示上翼面，負(fù)值表示下翼面。NACA0012位置示意圖可參考圖3。從圖10和圖11可以看出，在翼型前緣處水滴收集量最大：受攻角的影響，水滴收集量的最大值出現(xiàn)在翼型前緣下沿；而在前緣點(diǎn)之后沿著翼型表面水滴收集量逐漸減小，下翼面水滴收集量變化相較于上翼面更為緩慢，水滴撞擊極限也相對(duì)更大。與此同時(shí)，水滴收集量極值以及撞擊極限受LWC、MVD、溫度、高度、來流速度等因素共同影響。水滴收集系數(shù)極值受LWC和來流速度的影響較大，如Case 1與Case 4。而由于直徑（或質(zhì)量）較大、速度較大的液滴具有較大的慣性，在撞擊機(jī)翼表面時(shí)液滴軌跡不易發(fā)生改變，因此水滴撞擊極限范圍通常隨MVD和來流速度增大而增加，如Case 4與Case 5。

圖9 采樣點(diǎn)分布Fig. 9 Distribution of samples

圖12為水滴的撞擊范圍，橫坐標(biāo)表示翼面位置X軸坐標(biāo)與弦長(zhǎng)C的比值，NACA0012位置示意圖可參考圖3。撞擊范圍內(nèi)一共139個(gè)點(diǎn)，40組工況則可以構(gòu)成139×40的矩陣，對(duì)該矩陣進(jìn)行POD降階，獲得對(duì)應(yīng)的特征向量及系數(shù)。

圖13給出了壁面處水滴收集量最大誤差隨積累特征值的變化。圖中可以看出隨著積累特征值的增加，最大誤差迅速減小，當(dāng)積累特征值達(dá)到20時(shí)，水滴收集量最大誤差減小至接近于0。說明采用POD方法可利用少量特征基函數(shù)對(duì)原物理場(chǎng)進(jìn)行重構(gòu)。因此在本文中使用前20個(gè)特征基函數(shù)即可滿足精度要求。對(duì)于非采樣點(diǎn)（即預(yù)測(cè)點(diǎn)）處的系數(shù)，需要通過Kriging多維插值獲得。即針對(duì)非樣本點(diǎn)，需要對(duì)前20個(gè)系數(shù)分別進(jìn)行Kriging多維插值，獲得20個(gè)對(duì)應(yīng)的系數(shù)。本文的Kriging插值使用MATLAB的OODACE工具箱完成。最終，壁面的水滴收集量分布可以由20個(gè)特征基函數(shù)及對(duì)應(yīng)的系數(shù)相乘后獲得。

表3 部分采樣點(diǎn)工況Table 3 Conditions of part samples

圖10 Case 1~Case 5水滴收集量云圖Fig. 10 Water droplet collection amount contours of Case 1-Case 5

圖11 Case 1~Case 5水滴收集量分布Fig. 11 Water droplet collection amount distributions of Case 1-Case 5

圖12 水滴撞擊范圍示意圖Fig. 12 Schematic diagram of water droplet impinge?ment extent

圖13 水滴收集量最大誤差隨積累特征值變化Fig. 13 Maximum error variation of water droplet col?lection amount with cumulative eigenvalues

選用9組工況進(jìn)行驗(yàn)證，具體驗(yàn)證算例工況如表4所示。圖14為數(shù)值模擬與代理模型水滴收集量對(duì)比。圖中可以看出無(wú)論是水滴收集量的極值還是水滴收集量的分布范圍，快速預(yù)測(cè)結(jié)果均與數(shù)值模擬結(jié)果符合良好，9組算例水滴收集量最大值的平均絕對(duì)誤差MAE為0.002、均方根誤差RMSE為0.003，平均絕對(duì)百分比誤差MAPE為6.2%，說明所構(gòu)建的代理模型可以準(zhǔn)確地預(yù)測(cè)翼型表面的水滴收集量。

表4 驗(yàn)證工況Table 4 Validation conditions

圖14 數(shù)值模擬與代理模型水滴收集量對(duì)比Fig. 14 Comparison of water droplet collection amount between simulations and surrogated models

3.2 結(jié)冰氣象預(yù)測(cè)

用WRF模擬了華盛頓山（44.267N，71.3W，海拔1 917 m）2013年4月1日結(jié)冰過程。WRF的詳細(xì)信息可參考2.2節(jié)。圖7和圖8顯示了中心點(diǎn)的溫度和LWC的比較，結(jié)果顯示模擬獲得的溫度與LWC與觀測(cè)值良好符合。圖15為模擬區(qū)域和不同高度層的示意圖。圖16和圖17為4月1日6：00時(shí)高度層eta=2，4，6的LWC和溫度云圖?？梢钥闯?，LWC和溫度在不同層的變化很大。圖18和圖19分別為6：00、12：00、18：00和24：00時(shí)eta=6的LWC云圖和溫度云圖。圖中顯示LWC和溫度在1 d中均發(fā)生了很大變化。

圖15 模擬區(qū)域及高度層示意圖Fig. 15 Schematic diagram of simulated area and height layers

圖16 6：00時(shí)高度層eta=2，4，6的LWC云圖Fig. 16 LWC contours of layers eta=2，4，6 at 6：00

圖17 6：00時(shí)各高度層eta=2，4，6的溫度云圖Fig. 17 Temperature contours of layers eta=2，4，6 at 6：00

圖18 不同時(shí)刻下高度層eta=6時(shí)LWC云圖Fig. 18 LWC contours of layer eta=6 at different time

3.3 結(jié)冰試飛空域確定

3.3.1 結(jié)冰試飛空域

以WRF氣象模擬獲得的華盛頓山此次結(jié)冰事件中4月1日0：00—24：00時(shí)的氣象數(shù)據(jù)為基礎(chǔ)，假定飛行速度為50 m/s，利用水滴收集量代理模型對(duì)目標(biāo)區(qū)域內(nèi)的水滴收集量進(jìn)行快速預(yù)測(cè)，獲得了目標(biāo)區(qū)域內(nèi)在6：00、12：00、18：00以及24：00時(shí)刻下水滴收集量分布。以中度結(jié)冰強(qiáng)度對(duì)應(yīng)的水滴收集量為閾值對(duì)各個(gè)時(shí)刻下的空間點(diǎn)進(jìn)行識(shí)別。結(jié)冰強(qiáng)度指單位時(shí)間內(nèi)機(jī)體表面所形成的冰層的厚度，美國(guó)聯(lián)邦航空條例和航空信息手冊(cè)將結(jié)冰強(qiáng)度劃分為4個(gè)等級(jí)，用以說明結(jié)冰情況的嚴(yán)重性。將水滴收集量除以冰的密度即可獲得可能的最大增長(zhǎng)速率。由于本文主要集中于自然結(jié)冰試飛空域的確定，并未考慮結(jié)冰過程以及結(jié)冰后的形態(tài)，在后續(xù)的研究中將加入溫度對(duì)冰增長(zhǎng)速率的影響。在本文的研究中，冰的密度取917 kg/m3，且認(rèn)為冰的密度不變。結(jié)冰強(qiáng)度等級(jí)見表5［4］。

圖19 不同時(shí)刻下高度層eta=6時(shí)溫度云圖Fig. 19 Temperature contours of layer eta=6 at differ?ent time

表5 結(jié)冰強(qiáng)度等級(jí)［4］Table 5 Icing intensity rating［4］

取中度結(jié)冰時(shí)水滴收集量0.016 8 kg/（s·m2）作為劃分結(jié)冰試飛空域的判定標(biāo)準(zhǔn)，認(rèn)為最大水滴收集量超過該值時(shí)，該空域具有結(jié)冰試飛氣象條件，低于該值時(shí)則認(rèn)為不具有結(jié)冰試飛氣象條件。對(duì)目標(biāo)區(qū)域內(nèi)的空間點(diǎn)進(jìn)行篩選，獲得不同高度層在6：00、12：00、18：00以及24：00滿足條件的點(diǎn)的總數(shù)與高度層的曲線圖，如圖20所示。圖中可以看出在6：00、12：00和18：00時(shí)刻下均存在結(jié)冰試飛氣象條件，而在24：00時(shí)沒有結(jié)冰試飛氣象條件存在，曲線呈一條直線；不同時(shí)刻下滿足結(jié)冰試飛氣象條件的空間點(diǎn)的數(shù)量也不相同，滿足條件的點(diǎn)數(shù)隨著高度的增加先增加后減少，在大約高度為2~3 km的范圍內(nèi)達(dá)到峰值。圖中指出在12：00時(shí)刻下eta為5~7時(shí)滿足結(jié)冰試飛氣象條件的點(diǎn)數(shù)較高。圖21和圖22為12：00時(shí)刻下eta為5~7時(shí)的水滴收集量云圖和結(jié)冰強(qiáng)度云圖。圖中顯示在eta為5~7時(shí)存在一個(gè)水滴收集量較大的核心區(qū)域，如黑色線框所示。在核心區(qū)域外，水滴收集量或者結(jié)冰強(qiáng)度較低。不同高度的核心區(qū)坐標(biāo)平均高度及平均結(jié)冰強(qiáng)度見表6。

圖20 在6：00、12：00、18：00和24：00時(shí)滿足自然結(jié)冰試飛條件的空間點(diǎn)總數(shù)隨高度層的變化Fig. 20 Variation of total number of points meeting natural icing flight test conditions with height layers at 6：00、12：00、18：00 and 24：00

圖21 12：00時(shí)滿足結(jié)冰試飛條件的高度層水滴收集量云圖Fig. 21 Water droplet collection amount contours of lay?ers meeting conditions of icing flight test at 12：00

圖22 12：00時(shí)滿足結(jié)冰試飛條件的高度層結(jié)冰強(qiáng)度云圖Fig. 22 Icing intensity contours of layers meeting con?ditions of icing flight test at 12：00

3.3.2 飛行速度對(duì)自然結(jié)冰試飛空域的影響

飛機(jī)的結(jié)冰受到多種因素的影響，包括氣象因素和飛行條件。氣象因素中溫度、液態(tài)水含量以及水滴中值體積直徑是不可控因素，而飛行條件中的飛行速度和飛行時(shí)間是可以根據(jù)需要進(jìn)行調(diào)節(jié)的。本文通過對(duì)華盛頓山4月1日6：00~24：00時(shí)飛行速度為50 m/s和80 m/s時(shí)的水滴收集量的變化探究飛行速度對(duì)結(jié)冰試飛空域選擇的影響。

表6 核心區(qū)坐標(biāo)及平均結(jié)冰強(qiáng)度Table 6 Core area coordinates and average icing intensity

圖23為4月1日6：00、12：00、18：00以及24：00時(shí)刻不同高度層滿足結(jié)冰氣象條件的空間點(diǎn)總數(shù)對(duì)比。圖中可以看出對(duì)于6：00、12：00、18：00飛行速度為80 m/s時(shí)滿足結(jié)冰氣象條件的空間點(diǎn)數(shù)相比飛行速度為50 m/s時(shí)均有所增加；而在24：00時(shí)刻時(shí)滿足結(jié)冰氣象條件的空間點(diǎn)數(shù)在各個(gè)高度層上仍然均為0，說明在24：00時(shí)刻時(shí)目標(biāo)區(qū)域不具有滿足結(jié)冰試飛的氣象條件。圖中也可以看出，在6：00、12：00、18：00時(shí)刻下在eta為6時(shí)滿足結(jié)冰試飛氣象條件的空間點(diǎn)總數(shù)均達(dá)到最大值。

圖23 不同速度下滿足自然結(jié)冰試飛條件的點(diǎn)總數(shù)對(duì)比Fig. 23 Comparison of total number of points meeting natu?ral icing flight test conditions at different speeds

圖24 不同速度下高度層eta=6時(shí)水滴收集量云圖對(duì)比Fig. 24 Comparison of water droplet collection amount contours at different speeds of layer eta=6

圖24為eta為6時(shí)6：00、12：00、18：00下速度為50 m/s和80 m/s時(shí)滿足結(jié)冰試飛氣象條件的空間點(diǎn)分布。圖中可以看出由于速度的增加，目標(biāo)區(qū)域中滿足結(jié)冰試飛氣象條件的空間點(diǎn)均有所增加，在速度為50 m/s情況下目標(biāo)區(qū)域不具有結(jié)冰氣象條件的點(diǎn)，在速度為80 m/s的情況下出現(xiàn)了滿足結(jié)冰試飛氣象條件的點(diǎn)，如圖24（b）所示。此外，由于速度的增加，導(dǎo)致水滴收集量在一些區(qū)域內(nèi)過大，如圖24（b）、圖24（d）以及圖24（f）中的深色區(qū)域，這些區(qū)域內(nèi)的結(jié)冰強(qiáng)度可能達(dá)到嚴(yán)重結(jié)冰，若在這些區(qū)域內(nèi)進(jìn)行結(jié)冰試飛，可能導(dǎo)致事故發(fā)生。

4 結(jié) 論

1） 通過POD降階模型和Kriging插值方法構(gòu)建了水滴收集量的代理模型，代理模型可以快速預(yù)測(cè)結(jié)冰氣象參數(shù)（溫度、LWC、MVD）和飛行參數(shù)（速度、高度）對(duì)水滴收集量分布的影響；將水滴收集量等效轉(zhuǎn)換為結(jié)冰強(qiáng)度，作為試飛空域確定的劃分依據(jù)。

2） 通過WRF模式對(duì)2013年4月1日華盛頓山處一次結(jié)冰事件的氣象預(yù)測(cè)，獲得目標(biāo)區(qū)域內(nèi)結(jié)冰氣象條件；基于水滴收集量的代理模型快速獲得目標(biāo)區(qū)域內(nèi)水滴收集量分布，利用試飛空域劃分依據(jù)獲得了適合結(jié)冰試飛的空域、時(shí)間以及結(jié)冰強(qiáng)度。

3） 對(duì)速度對(duì)自然結(jié)冰試飛區(qū)域的影響進(jìn)行研究，發(fā)現(xiàn)速度的增加使得水滴收集量增加，進(jìn)而使得原本沒有結(jié)冰試飛氣象條件的位置出現(xiàn)了滿足結(jié)冰試飛的條件，也可能導(dǎo)致原本滿足結(jié)冰試飛條件的位置水滴收集量過大從而威脅飛行安全。