?？〗?，桑為民，2，*，李棟，郝蓮，王澤林

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

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

3．中國商用飛機有限責任公司 上海飛機設計研究院，上海 201210

4．中國飛行試驗研究院，西安 710089

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

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

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

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

1 計算方法及控制方程

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

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

1.1 水滴撞擊特性求解

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

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

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

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

水滴的運動過程可以表示為

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

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

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

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

1.2 POD方法

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

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

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

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

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

因此，重構的物理場可以使用正交基的線性疊加表示為

1.3 Kriging插值

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

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

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

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

使用協(xié)方差來表示不同相關點處隨機變量的相關性，即

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

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

2 數(shù)值模擬驗證

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

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

2.1 NACA0012水滴撞擊特性驗證

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

圖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氣象模擬驗證

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

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

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

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

圖7為華盛頓山溫度隨時間的變化與文獻［31］結果的對比。其中黑色原點為觀察值，實線為本文模擬結果，虛線為文獻［31］中的模擬結果。

圖8為華盛頓山LWC隨時間變化與文獻［31］結果的對比。與圖7類似，黑色原點為觀察值，實線是本文模擬結果，虛線是文獻［31］中的模擬結果。

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

圖7 溫度隨時間變化對比Fig. 7 Comparison of temperature variation with time

圖8 LWC隨時間變化對比Fig. 8 Comparison of LWC variation with time

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

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

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

表2 結冰氣象預測誤差分析Table 2 Error analysis of icing weather prediction

3 結果與討論

3.1 水滴收集量代理模型

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

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

圖9 采樣點分布Fig. 9 Distribution of samples

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

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

表3 部分采樣點工況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組工況進行驗證，具體驗證算例工況如表4所示。圖14為數(shù)值模擬與代理模型水滴收集量對比。圖中可以看出無論是水滴收集量的極值還是水滴收集量的分布范圍，快速預測結果均與數(shù)值模擬結果符合良好，9組算例水滴收集量最大值的平均絕對誤差MAE為0.002、均方根誤差RMSE為0.003，平均絕對百分比誤差MAPE為6.2%，說明所構建的代理模型可以準確地預測翼型表面的水滴收集量。

表4 驗證工況Table 4 Validation conditions

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

3.2 結冰氣象預測

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

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

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

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

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

3.3 結冰試飛空域確定

3.3.1 結冰試飛空域

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

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

表5 結冰強度等級［4］Table 5 Icing intensity rating［4］

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

圖20 在6：00、12：00、18：00和24：00時滿足自然結冰試飛條件的空間點總數(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時滿足結冰試飛條件的高度層水滴收集量云圖Fig. 21 Water droplet collection amount contours of lay?ers meeting conditions of icing flight test at 12：00

圖22 12：00時滿足結冰試飛條件的高度層結冰強度云圖Fig. 22 Icing intensity contours of layers meeting con?ditions of icing flight test at 12：00

3.3.2 飛行速度對自然結冰試飛空域的影響

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

表6 核心區(qū)坐標及平均結冰強度Table 6 Core area coordinates and average icing intensity

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

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

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

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

4 結 論

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

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

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