基于降階模型的翼型結冰冰形預測方法

劉藤， 李棟， 黃冉冉， 張振輝

(西北工業大學航空學院, 西安 710072)

隨著航空業的不斷發展,飛機的應用也從早期的軍事用途逐漸擴展到飛行表演、緊急救援、商業運輸等多個領域。這就要求飛機可以適應各種復雜氣象條件完成預定飛行任務[1]。

飛機結冰現象主要發生在飛行中遇到含有過冷水滴的云層時，過冷水滴撞擊到機體迎風表面就有可能發生結冰現象[2]。通過研究發現，結冰對飛機的影響一方面是結冰后飛機覆冰會造成氣動外形破壞，形成的表面粗糙度會改變飛機的的原有流場狀態，從而造成阻力增加、升力減小、飛行穩定性降低等;另一方面還會造成儀表失靈，引發飛行員錯誤操作，嚴重時可能導致墜機[1]。

隨著計算機領域的飛速發展，研究及預測結冰問題已經可以利用數值方法來完成結冰模擬。然而利用數值手段模擬結冰問題需要分別得到機翼繞流流場和過冷水滴流動情況以及求解結冰熱力學方程[3]，因此結冰模擬也是比較復雜的,對于數分鐘的結冰過程通常需要較長時間才能得到模擬結果。

目前有2個因素需要開展結冰冰形快速預測方法研究。一方面，適航認證結冰試飛過程中，需要針對各個飛行狀態的結冰情況實時預測其對性能的影響，從而判斷防除冰系統的開啟[4]；另外即便現有飛機防冰系統已經比較完善,但除防冰系統并不能在各種狀態下都達到理想的效果,很多情況下依然需要帶冰飛行，考慮在設計階段就討論結冰情況的影響，在保證干凈翼型氣動性能的同時，提升結冰條件下翼型的氣動性能的容冰優化設計很有必要。 另一方面，考慮優化設計過程反復迭代時也需要發展快速預測不同外形的不同結冰條件下結冰冰形的方法。因此本文擬根據容冰優化設計需要，建立基于本征正交分解(POD)的結冰冰形快速預測降階模型。

考慮到飛機結冰問題的高度非線性，以及結冰數值模擬的時間成本，直接對CFD計算原始樣本數據進行建模構建代理模型需要較大的樣本數目來保證模型的精度，而降階模型在保證模擬精度同時可以快速預測結果。因此如果利用相關數據處理方法對樣本數據進行降階處理，在低階空間內對樣本數據進行建模會取得更為高效的效果。

降階模型是利用復雜的數學工具對大型復雜問題進行低階空間重構，將問題簡化為低階空間元素的線性組合，目前已經有Krylov子空間方法[5]、小波分析[6]、雙正交小波[7]和POD[8]等方法。POD方法基本原理[9]是利用已有樣本集提取一系列正交基向量，原樣本空間可以由這些向量的線性疊加來描述。然后便可以根據所需要的精度舍棄一些低能量模態，從而形成低階子空間，可以快速進行樣本空間重建，因此POD方法可以應用于快速預測飛機結冰冰形以及評估氣動性能損失。Nakakita等[10]首先將POD方法應用于飛機結冰冰形預測，并取得了較好的效果，但其采用的插值基于Akima插值法只能應用于兩參數情形。

本文介紹了建立基于POD方法和Kriging模型的翼型結冰冰形預測降階模型方法。首先，詳細介紹飛行迎角變化時的降階模型預測結冰冰形的實現步驟；然后，對結冰試飛認證框架下的五參數的冰形計算結果進行分析，將降階模型預測的冰形結果與CFD模擬得到的結果進行比較，分析驗證預測方法的可靠性；最后，進一步提出了基于更為先進的貝葉斯框架下的Blind-Kriging模型的改進預測模型的方法，建立了降階模型對結冰冰形進行預測，可以快速、準確地得到樣本空間內的任意非樣本點冰形。

1 基本方法

建立降階模型通常首先需要完成初始空間采樣，然后對樣本數據降階分析，最后完成數據插值等過程。具體流程如圖1所示。

圖1 降階模型建立流程Fig.1 Reduced order model building process

1.1 參數空間采樣

對于優化設計而言，合理的參數空間采樣是構建代理模型的關鍵，即采樣需要選擇合適的試驗設計方法。試驗設計概念最早是由英國生物統計學家Fisher于 20世紀20 年代所提出的[11]，目前已經發展了許多不同的方法和技術，比如基于拉丁超立方抽樣(LHS)方法、Halton序列方法、最小差異序列方法、質心Voronoi鑲嵌等。

LHS方法是由Mckay 等[12]在1979年提出的試驗設計方法，該方法的基本原理是：如果需要得到n個樣本點，則利用m個設計變量將參數空間均勻劃分為等概率的n個區間，這樣整個參數設計空間被分為nm個小區間，然后對其進行n次取樣。設計空間抽樣時一方面需要保證樣本點在每個小區間上都是隨機選取的；另一方面需要各個樣本點在任一維度上進行投影，每個小區間中只會有一個樣本[12]。通常取樣的算法為

(1)

(2)

其中:I為單位向量。抽樣具有隨機性，因此LHS 選取的樣本有的分布很好，有的分布很差。該方法有很多改進方法，包括中心化LHS、最優LHS[13]等。MATLAB 的LHS 已經實現了基于極大極小準則的最優LHS，因此本文采用該試驗設計方法，并將試驗點均放在小立方體的中心。

1.2 POD方法

POD方法是基于物理場原始數據快速準確地重構物理場，同時也可預測物理場在某些未知位置的數據的一種數值方法。POD方法的核心思想是利用特殊正交基函數的迭加來實現參數的近似[14]。它最初主要用于湍流擬序結構的分析研究中， POD如今已成功地應用于各類領域，如機器學習中的主成分分析(Principal Components Analysis,PCA)法，海洋學和氣象學中的經驗正交函數法，心理學和經濟學中常用的要素分析法等[15]。

POD的基本理論[16]為假設存在這樣的物理場y=y(x,t),在某個空間可以表示為如式(3)所示無窮級數:

(3)

式中：tn為參數變量，當前研究中可以為速度、壓強、升力、阻力、冰形坐標等；αk(tk)為經驗系數；φi(x) 為特征基函數。

利用相關的數學理論將式(3)求解轉換為求特征值問題：

(4)

式中：〈〉為內積運算符。對于此問題直接求解需要消耗巨大的資源和時間，所以實際應用中常采用由Sirovich提出的稱為快照(snapshot)[16]的方法。該方法可以將式(4)的高維矩陣求解特征值問題轉化為一個與快照個數N相等的特征值問題。若一個樣本包含M個數據，則取N個樣本組成一個樣本集U={U1,U2,…,UN} 這是一個M×N階的矩陣。該方法提出特征基函數可以由快照的線性組合構成，即

(5)

Aσ(n)=λnσ(n)n=1,2,…,N

(6)

其中：λn為特征值；A的元素定義如下：

(7)

其中特征值的大小可以表示該正交基函數在整個樣本空間中所占能量的比重，即第s個最佳正交基方向上的物場理能量Es占總能量Etotal的比重εs為[4]

(8)

式中：λs為第s個特征值。因此進行物理場的重建便可以根據能量的大小對基底進行選擇。并不一定需要選擇所有的基底來完成重建，假設M為所用正交基底的個數，有M≤N。一般M取得過大反而會受到舍入誤差影響影響重建效果。在進行矩陣分解得到最佳正交基之后，對選取的正交基底與樣本進行內積便可得到對應的系數向量：

(9)

本研究將POD理論用于翼型的二維冰形計算，主要考慮飛行迎角、飛行速度和結冰時間等多種物理參數對結冰冰形的影響。首先，利用不同狀態下冰形輪廓點坐標值來構建樣本矩陣；然后，利用前述的快照方法得到對應的快照矩陣A，對其進行矩陣特征分解得到對應的特征向量，也就得到了構建子空間的正交基及其系數向量。于是設計空間內的任意計算狀態下的結冰冰形可以由這些正交基線性疊加表示為

(10)

1.3 Kriging模型

假定已有N個快照的系數α定義了RK中N個超平面α(x)的值，其中K為自由參數的數目，x代表設計變量x=[x1,x2,…,xK]T，尋找α(x) 的問題也就成為了超平面中的多維插值問題，在大多數的降階模型應用中多采用多項式插值來決定系數，而在本研究中采用Kriging插值方法，Kriging插值方法可以接受任意數目的設計量，保證足夠的插值精度而不會對計算時間造成影響，可以應用于任意維度的插值過程。

Kriging模型[17]是由南非工程師克里金(Krige)于20世紀50年代提出的,現已廣泛應用于各個領域，本文以普通克里金插值方法為例來說明克里金插值的原理[18]：

(11)

只要得到了ω=[ω(1),ω(2),…,ω(n)]T的表達式，就可以得到參數空間中的任意設計點的估計值。Kriging模型的重要假設是將未知函數視為一個高斯隨機過程,該隨機過程可以表示為

Y(x)=F(β,x)+H(x)

(12)

式中：x為參數變量；F(β,x)稱為全局趨勢模型；H(x)為隨機分布的誤差。將隨機過程的方差記為σ2，可用協方差描述不同位置點這些隨機變量的相關性，其協方差可描述為

Cov(H(x),H(x′))=σ2Z(x,x′)

(13)

其中：Z(x,x′)為相關函數(只與空間距離相關)，相關性會隨距離增大而不斷減小。相關函數通常有指數函數、高斯函數、三次樣條函數等多種形式。

基于上述假設，Kriging模型尋找最優加權系數ω，使得均方差最小，并且滿足無偏差條件。

通常利用拉格朗日乘子法可以推導得出一個線性方程組，求解該方程組便可得到Kriging模型預估值為

(14)

式中：Vkrig只與已有樣本點相關；β0=(FTZ-1F)-1FTZ-1yk，yk為樣本數據，F為一維單位向量，Z和z則是前述的相關函數矩陣及其中的相關向量，故便可利用樣本數據得到模型的估計值。

構建Kriging模型時，還可以將對其模型參數(或超參數)也作為未知量來加入訓練，可以進一步增強代理模型的靈活性。利用已有樣本點和響應值可將問題轉化為一個非線性無約束優化問題，一般采用“最大似然估計”或“交叉驗證”得到超參數。對于最大似然估計而言，就是選擇超參數使式(15)函數值達到最大[19]：

(15)

式中：β0和σ2的最優值可以給出解析解為

(16)

將β0和σ2代入后，取對數忽略常數項，優化問題轉化為

(17)

利用數值優化算法如遺傳算法、單純形法、模式搜索法等訓練超參數，就可以構建最終的Kriging模型。

2 冰形模擬及樣本獲取

本文樣本的獲取依靠CFD的模擬獲得，結冰的CFD模擬一般包括以下步驟：①空氣流場計算；②過冷水滴或者冰晶、凍雨等流動情況求解及其撞擊特性計算；③求解結冰熱力學模型計算結冰量；④根據原有幾何外形加上結冰量，完成冰形計算及幾何形狀重構。

本文利用FENSAP-ICE軟件[20]來進行結冰數值模擬，首先對基于雷諾平均的Navier-Stokes方程來求解得到空氣流場，然后采用歐拉法計算模型表面的水滴撞擊特性，最后求解結冰熱力學模型計算結冰冰形。將結冰計算的步驟分為各個模塊，其中FENSAP模塊用于計算空氣外流場，DROP3D模塊用于計算水滴撞擊特性，ICE3D用于熱力學方程求解及冰形計算，模塊相互關系如圖2所示。

本文選用了NACA0012翼型為研究對象，根據相關文獻選取了指定計算狀態，用上述CFD方法先后調用各個模塊進行結冰模擬，模擬結果得到的冰形如圖3所示。

從圖3中可以看出，冰形的輪廓點可以描述完整的冰形。因此可以取冰形上各個輪廓點的二維坐標值做為樣本數據來進行后續計算。圖3翼型上表面網格點數為277，所以該算例一個樣本包含277個數據。

對于降階模型的建立而言，樣本數據的保真度對于模型的精度也是至關重要的，這里將上述模擬結果與實驗數據進行了對比驗證，對比結果如圖4所示，c為翼型弦長，表明了本文采用的數值模擬方法的可靠性。

圖3 FENSAP-ICE冰形模擬結果Fig.3 FENSAP-ICE ice shape simulation results

3 降階模型應用算例與結果分析

3.1 單參數冰形預測

本文先以單參數冰形預測為例，詳細描述構建降階模型的實施步驟。算例選擇NACA0012二維翼型為研究對象，結冰冰形樣本的獲得是通過上述CFD方法來完成的。考慮飛行過程中的迎角變化對結冰冰形的影響來建立POD預測模型。選取的計算條件如表1所示。

對于翼型而言，迎角的變化對氣動性能有著重要的影響，同樣它對結冰冰型也會產生很大的影響，分別選取了0°、1°、2°、3°和4°的結冰冰型作為樣本集進行研究(樣本個數為5)，CFD模擬后得到的樣本集冰形結果如圖5所示。

從圖5中可以看出，隨著迎角的增加，翼型上表面的結冰范圍越來越小，而下表面結冰范圍越來越大。這是由于隨著迎角的增大造成翼型下表面液滴撞擊區域擴大，因此結冰范圍向翼型下表面移動。本文對樣本空間中的結冰冰形樣本分別取出其冰形輪廓點的x、y坐標構成樣本矩陣，樣本矩陣的階數為277×5。

表1 算例計算條件Table 1 Calculation conditions of example

利用數值模擬得到的樣本矩陣采用快照方法進行簡化得到對應的快照矩陣。這里一共選取了5個樣本，故得到的快照矩陣為5×5階。

圖5 不同飛行迎角的冰形Fig.5 Ice shapes of different flying attack angles

通過上面的x坐標快照矩陣可以看出，上述矩陣的元素差別很小，這是由于不同的冰形樣本坐標大部分位置是相同的，只是在結冰位置略有差別，同時冰形的變化量相對于整個坐標來講是小量，所以計算過程需要保證一定的數據精度。

在得到快照矩陣后，分別計算x、y的特征值λx、λy與特征向量(見表2)，可以看到第1個特征值明顯大于其他，說明樣本集的數據相似度很高。利用獲得的特征向量便可以計算得到正交基，選取全部5個正交基進行POD冰形預測，可以得到每個樣本在不同正交基上的系數如下表所示，可以看到第1個正交基明顯系數大于其他，表明第1個正交基對POD計算結果的影響非常大。

通過進行上述計算，基本建立了POD模型，這里利用前述的Kriging模型進行插值。選取迎角為2.5°的非樣本點進行結果驗證，最終的POD預測x、y坐標與CFD的計算結果對比如圖6所示，從圖中可以看出POD的預測結果與CFD結果基本一致，只是在冰角以及結冰極限位置存在一些誤差，說明POD預測方法是可行有效的。

表2 快照矩陣特征值Table 2 Eigenvalue of snapshot matrix

圖6 單參數POD與CFD冰形對比Fig.6 Ice shape comparison between single-parameter POD and CFD

3.2 多參數冰形預測

將POD方法擴展到多維情形，這里選取通常影響結冰的5個參數迎角、速度、溫度、MVD、高度(選擇范圍如表3所示, 1°F=-17.2℃,1 ft=0.304 8 m)。結冰時間仍選擇為360 s，另外LWC參數的選取考慮到實際的工程應用，按照連續最大結冰條件(見圖7)[21](FAR-25附錄C)得到對應狀態的液態水含量。

考慮到多參數樣本選取的復雜性，為了保證一定的預測精度，本文選取了100個樣本(見圖8)，其樣本空間的采樣采用LHS方法。

LHS方法是一種“充滿空間的設計”，選取的樣本點可以均勻地散布于參數空間。

與3.1節的單參數POD方法類似，多參數情形同樣是利用上述樣本空間的結冰冰型結果，可得到對應的100組樣本二維坐標參數，然后通過POD方法得到對應的最佳正交基。由于多參數預測選取的樣本數較大也會得到較多的正交基，這里只需要從中選取一些所占能量較大的正交基進行樣本的重建即可。當選取的結冰冰形的預測正交基所占的能量總和如果大于99.9%，便可以認為它能完全反映樣本冰形的主要特征[4]。然后再利用Kriging模型對重建系數進行插值，得到POD系數及最終的冰形預測結果。

表3 多參數樣本選擇范圍Table 3 Selection range for multiparameter sample

圖7 連續最大結冰條件[21]Fig.7 Continuous maximum icing condition[21]

利用上述的方法建立預測模型，選取的算例參數如表4所示。

圖9為多參數的POD結果與CFD結果的比較，其中POD結果選取了能量總和占比達99.9%的最佳正交基計算得到，可以看出POD的預測結果與CFD結果吻合較好。

圖8 樣本空間分布Fig.8 Sample space distribution

算例狀態數值飛行迎角/(°)1.325飛行速度/(m·s-1)91.875MVD/μm25.5高度/ft12031結冰溫度/°F16結冰時間/s360

圖9 多參數POD與CFD冰形對比Fig.9 Ice shape comparison of multiparameter’ POD and CFD

4 插值模型

4.1 Blind-Kriging模型

傳統Kriging模型采用確定的趨勢模型，為了更好地改善近似精度, 文獻[22]研究了一種利用貝葉斯框架來識別趨勢模型的改進代理模型方法，稱為“Blind-Kriging”模型。該模型的趨勢模型不是確定的已知函數而是通過數據分析方法獲得，具體理論為

Y(x)=f(x)Tβm+H(x)

(18)

(19)

Blind-Kriging最重要的步驟是識別出未知函數fi，這是通過函數(變量)選擇技術在簡單函數構成的候選函數集中選擇得到，然后就可以得到預測模型的第一部分f(x)Tβm。

構建Blind-Kriging模型可以被看作是1個兩階段的過程。在第1階段，建立一個普通Kriging模型，即具有常數回歸函數的Kriging模型，完成模型構建和超參數θ估計。在第2階段擴展該初始Kriging模型的回歸函數，根據估計的α系數選擇新的特征，生成一系列中間Kriging模型。當交叉驗證錯誤滿足要求時，中間Kriging停止迭代，然后選擇最佳的特征組來構建最終Blind-Kriging模型，重新估計超參數θ。

比利時根特大學的研究人員已經開發了OODACE工具箱可以實現Blind-Kriging模型構建，本文的工作也是基于此工具箱開展的，其工具箱偽代碼[23]如下。

X←samples

b1,…,bt{Candidate features}

b=Cte{Selected features}

θ0=maxθlikelihood(X,θ,b)

M0=construct(X,θ,b){Construct ordinary Kriging model}

α0=evaluateMeasure(M0){Assess accuracy}

i=0

while improvement(α){Accuracy improves?}

i=i+1

β=rank(b1,…,bt)

j=maxj(|βj|)

b=b∪bj

θi=maxθlikelihood(X,θ,b){Optional}

Mi=update(Mi-1,θi,b){Intermediate Kriging model}

αi=evaluateMeasure(Mi){Assess accuracy}

endwhile

θfinal=maxθlikelihood(X,θ)

Mfinal=update(Mi,θfinal,b){Final Blind Kriging model}

4.2 模型應用

冰形預測仍選擇為NACA0012翼型，考慮方法研究的目的，這里只選取參數為迎角與飛行速度，進行兩參數的POD冰型預測，該方法同樣可以擴展到多更多參數的應用中。分別選取了迎角為0°、1°、2°、3°和4°和飛行速度為90、100、110、120和130 m/s的結冰冰型作為樣本集進行研究(M=25)。該算例選取的計算條件見表5。

利用本文1.2節介紹的法進行了POD分析，并分別采用傳統Kriging與Blind-Kriging模型進行預測。選取迎角為2.5°，飛行速度為115 m/s的非樣本點進行模型對比，最終的POD預測與CFD的計算結果對比如圖10所示。從圖中可以看出，POD的預測結果與CFD結果基本一致，不同插值模型的預測結果也比較類似，采用Blind-Kriging模型對于結冰冰型的粗糙度(非均勻性)預測未有很明顯改善，為了改善預測精度，更好地發揮Blind-Kriging模型特點，選擇在趨勢模型中引入更高階項，建立Blind-Kriging高階模型可以發現對于冰型預測相比于低階的Blind-Kriging模型預測結果精度有了較大改善，尤其對于結冰極限位置的預測相較于其他模型精度明顯提高，表明了Blind-Kriging模型對于復雜冰型預測有著較好的適應性，尤其采用高階模型能更好的發揮其先進性。

表5 計算狀態Table 5 Calculation conditions

圖10 冰形預測結果對比Fig.10 Comparison of ice shape prediction results

5 結 論

本文建立了一種基于POD和Kriging插值方法的降階模型，利用CFD數值模擬的冰形結果作為樣本來源，考慮了飛行迎角、飛行速度以及結冰溫度的等多種物理參數對冰形影響，先后完成了單參數、兩參數以及五參數情形下的翼型結冰冰形快速預測，得到：

1) 通過對POD預測結果與CFD結果進行對比可知：降階模型可以得到較好的冰形預測效果，只是在冰角附近及結冰極限位置預測結果有細微的差別，表明該方法應用于翼型結冰冰形的預測是準確、可行的。

2) 利用貝葉斯框架下的Blind-Kriging模型進行數據插值雖然需要更多的時間完成模型建立，但可以改善降階模型的精度，可以取得更為理想的冰形預測效果。

3) 雖然樣本的取得需要一定的計算時間，但是降階模型建立完成之后，便能夠快速、準確地預測一定范圍內任意狀態的冰形。

4) 本文方法對涉及大參數研究的應用極具價值，如果預先計算足夠的快照，則通過構建POD/Kriging逼近而不是使用全部的CFD解決方案來進行參數研究可以很容易節省大量時間。這為結冰試飛認證等工程應用提供了一種有效可行的方法，同時該方法可以推廣應用至翼型設計階段，對于容冰設計中的翼型結冰影響分析以及容冰優化設計是一種可靠準確的預測手段。