柴聰聰易 賢郭 磊王 俊
1.電子科技大學 計算機科學與工程學院,成都 611731;2.中國空氣動力研究與發展中心 結冰與防/除冰重點實驗室,四川 綿陽 621000
飛機在低于零度的環境中飛行,如果遇到含有過冷水滴的云層,且云層中的過冷水滴撞擊到飛機機翼表面,就會在碰撞區域及其附近發生結冰[1-2]。機翼結冰不僅會增加機翼部件的重量,還會改變機翼周圍的流場分布、破壞機翼的氣動性能[3],嚴重危害飛機的飛行安全。為確保飛機的飛行安全,防止墜機事故的發生,對翼型結冰特性進行快速、準確的預測非常有必要。
傳統的結冰預測方法主要有結冰試驗[4-5]和CFD計算[6-7]。其中,結冰試驗包括真實結冰和人工模擬結冰條件下的飛機飛行試驗以及冰風洞試驗。此類方法得到的冰形與真實條件相接近,但存在試驗周期長、試驗成本高等問題,甚至還伴隨一定的危險性[8]。相比之下,基于CFD計算的冰形預測成本低、安全性高,是科研人員常用的結冰預測方法。文獻[9]采用CFD計算對三段翼和MS-317后掠翼進行研究,結果顯示利用數值計算得到的三維積冰模型與真實數據相接近。雖然成本降低、安全性高,但CFD計算需要復雜的計算過程,預測精度易受模型的影響,無法實現翼型結冰的快速預測。
隨著人工神經網絡的出現和發展,人們開始利用人工神經網絡進行結冰預測[10-13]。文獻[10]采用BP神經網絡進行二維翼型結冰預測,預測結果表明:該方法可以快速預測冰形,具有較高的準確性。但是這種方法需要前緣近似轉換、保角映射等一系列復雜繁瑣的數學處理,而且前緣近似轉換只適用于對稱翼型,不適用非對稱翼型。針對此問題,文獻[13]提出利用翼面坐標轉換的方法進行冰形轉換,解決了前緣近似轉換無法用于非對稱翼型的問題。上述結冰預測研究均是對不同結冰條件下的翼型進行冰形定性預測,對冰形特征參數的預測較少。
冰形特征參數對臨界冰形確定、定量描述以及氣動特性具有重要影響。在冰形確定和定量描述方面,冰形特征參數中的冰角高度可以量化翼型結冰嚴重程度,結冰極限可直觀表現翼型表面結冰范圍等;在氣動特性方面,Kim[14]和Bragg教授[15]等進行了相關研究,結果表明:冰形特征參數的改變會對不同翼型氣動特性造成不同程度的影響,其中冰角高度及位置等冰形特征量影響最大。
本文利用BP神經網絡建立冰形特征參數預測模型,以飛行速度、結冰時間、液態水含量、環境溫度和平均水滴直徑作為輸入,以結冰極限、冰角高度和角度等冰形特征參數作為輸出。第一節介紹冰形特征參數,第二節描述樣本數據及數據相關處理,第三節介紹冰形特征參數預測網絡模型,第四節為仿真實驗和結果分析,第五節得出結論。
Ruff等[16]最早給出了結冰的幾何特征描述,主要包括駐點結冰厚度、最大結冰厚度、最大結冰寬度、撞擊極限寬度、冰角長度和冰角偏角等冰形特征參數(如圖1所示)。但是他們定義的冰形特征參數并沒有得到廣泛的應用。

圖1 結冰幾何特征描述[16]Fig.1 Description of icing geometric properties[16]
2012年美國汽車工程師協會(Society of Automotive Engineers,SAE)制定了有關結冰風洞試驗的一些標準[17],明確給出了用于定量描述二維翼型結冰幾何特性的冰形特征參數:結冰上極限Su、結冰下極限Sl、結冰面積Sice、駐點結冰厚度hsp、上冰角高度hu、下冰角高度hl、上冰角角度θu和下冰角角度θl,具體如圖2所示。其中:結冰上下極限為上下機翼表面結冰最遠位置到駐點的距離,上翼面為正,下翼面為負;駐點結冰高度為駐點沿法向的結冰高度;上(下)冰角為沿上(下)翼面法向結冰高度最大的冰角,對于飛機迎角不為零出現多個冰角的情況,規定上下冰角為靠近駐點且分別位于上下翼面結冰高度最大的冰角,對于只有一個冰角的霜冰,仍根據駐點位置劃分上下冰角(見圖3)。此外,上(下)冰角角度是指上(下)冰角頂點與翼型前緣之間連線的夾角;上(下)冰角高度為上(下)冰角頂點沿翼面法向到翼面的高度。

圖2 SAE規定的冰形特征參數示意圖Fig.2 Icing characteristic parameters provided by SAE

圖3 霜冰的冰形特征參數Fig.3 Icing characteristic parameters of glaze ice
結合Bragg教授的研究結果,本文主要對結冰上、下極限,上、下冰角高度和上、下冰角角度等這6個冰形特征參數開展預測研究。
本文采用二維NACA0012翼型結冰數據,主要利用文獻[18]中前緣結冰數值計算方法得到:首先利用SIMPLE方法求解低速黏流時均N-S方程進行流場計算,之后進行部件表面水滴運動及撞擊特性計算,最后基于改進的Messinger結冰熱力學模型開展結冰計算。圖4為計算與結冰風洞實驗得到的冰形對比圖,由圖可知該方法計算得到的冰形與實驗結果較為一致,誤差在合理范圍內,說明該方法可行。

圖4 數值計算與實驗冰形對比Fig.4 Comparison between the numerical calculation method and the experiment
基于以上數值方法,本文一共產生了900個二維翼型結冰原始數據,結冰條件有關設置如表1所示,其中飛行迎角保持不變。

表1 氣象和飛行條件參數設置Table 1 Meteorological and flight condition parameters setting
利用文獻[13]中的翼面坐標轉換方法處理上文得到的翼型結冰原始數據,提取冰形特征參數,并與對應的結冰條件參數組成訓練網絡用的有效樣本。其中每組樣本數據都由5個結冰條件的輸入向量和6個冰形特征參數的目標向量組成。
由于輸入向量的結冰條件量綱不同,不同條件參數之間的數量級相差較大,為加快網絡收斂速度,防止神經元出現輸出飽和現象,對輸入向量進行Zscore標準化處理:

式中,xij_std為標準化后的輸入向量第j個條件的第i個樣本值;xij為原始輸入向量第j個條件的第i個樣本值;xj_mean為第j個輸入條件的均值;xj_std為第j個輸入條件的標準差;n為樣本個數,n=900;m為輸入向量條件個數,m=5。通過上式將輸入向量標準化為均值為0、標準差為1的數據。
BP神經網絡是一種按照誤差逆向傳播算法(error back-propagation)訓練的前饋神經網絡,由輸入層、隱藏層和輸出層組成,具有較強的自適應性和良好的容錯特性,對于未訓練的輸入數據也會得到較為合適的預測結果[19]。BP神經網絡的訓練過程主要由信號前向傳播和誤差反向傳播兩部分組成:信號從輸入層輸入,經過隱藏層,最后到達輸出層;而誤差則從輸出層到隱藏層,最后到輸入層,依次調節隱藏層到輸出層的權重和偏置、輸入層到隱藏層的權重和偏置,多次訓練直至達到誤差閾值或訓練上限。
在網絡結構方面,本文采用如圖5所示的單層BP神經網絡作為預測模型,以影響翼型結冰的液態水含量(Liquid Water Content,LWC)、平均水滴直徑(Median Droplet Diameter,MVD)、結冰時間t、飛行速度v、環境溫度T等作為網絡的輸入,以冰形特征參數作為預測網絡模型的輸出。

圖5 單層BP神經網絡模型Fig.5 Single-layer BP network
在網絡訓練方面,本文采用適用于中等規模以下且訓練速度最快的LM(Levenberg-Marquardt)學習算法,以均方誤差(Mean Square Error,MSE)作為誤差性能函數。相比Sigmoid激活函數的收斂速度緩慢,Tanh函數的0均值更加有利于提高訓練效率,因此本文選擇Tanh函數作為網絡的激活函數,函數表達式為:

除此之外,隱藏層神經元個數對網絡的泛化性能有很大的影響,個數過多會導致過擬合、訓練時間過長;個數過少則會導致欠擬合,因此需要選擇出泛化性能較優模型對應的神經元個數。本文采用k折交叉驗證方法選擇隱藏層神經元的個數。
k折交叉驗證是將原數據集D劃分為k個大小相同的子集,每個子集通過原數據集分層采樣得到,之后每次用k-1個子集作為訓練集,余下的子集作為測試集,進行k次訓練和測試,最后返回k個測試結果的平均誤差值作為該模型的性能評估(如圖6所示)。

圖6 k折交叉驗證Fig.6 k-fold cross validation
本文采用10折交叉驗證進行網絡模型選擇和評估,隱藏層神經元個數選擇范圍為:4、6、8、10、12、14、16、18、20、22。圖7給出了計算得到的各個模型平均測試誤差。由圖可知,神經元個數為22的模型誤差最小,但神經元個數過多將導致過擬合,所以本文選擇隱藏層神經元個數為14的預測模型,并利用全部的樣本數據訓練模型,為后續的預測做準備。

圖7 不同神經元模型的10折交叉驗證平均測試誤差Fig.7 Average test error of 10 fold cross validation for different neuron models
本文選擇樣本數據之外的28組算例數據進行網絡預測,其中霜冰9組,明冰12組,混合冰5組。表2為選取的最具代表性的8組數據,其中算例1、2和3是霜冰,算例4、5和6是明冰,算例7和8是混合冰。圖8展示了這8個算例二維翼型結冰外形。

表2 算例的氣象和飛行條件參數Table 2 Meteorological and flight parameters of cases

圖8 各算例的翼型結冰結果Fig.8 Airfoil icing results of cases
冰形特征參數的預測誤差利用絕對百分比誤差(Absolute Percentage Error,APE)和平均絕對百分比誤差(Mean Absolute Percentage Error,MAPE)進行評估,具體公式如下:

其中,ypre為網絡預測的冰形特征參數,ycal為數值計算得到的冰形特征參數,N為算例的個數,本文取N=8。
表3為算例冰形特征參數數值計算結果C和網絡預測結果P,表4為冰形特征參數預測的絕對百分比誤差APE和平均絕對百分比誤差MAPE。

表3 冰形特征參數計算與預測結果Table 3 Calculation and prediction results of ice shape characteristic parameters
由表4可知,冰形特征參數平均絕對百分比預測誤差的排列順序為:下冰角高度>結冰上極限>結冰下極限>上冰角高度>上冰角角度>下冰角角度,其中冰形特征參數中的冰角角度平均絕對百分比預測誤差最小,MAPE在1.50%以下,結冰上、下極限次之,其MAPE分別為4.68%和3.32%,而下冰角高度最大,超過其他冰形特征參數之和,MAPE為18.28%。整體而言,該神經網絡模型的預測效果與數值模擬計算結果較為接近,說明了該方法的可行性。

表4 冰形特征參數預測誤差Table 4 Prediction error of ice shape characteristic parameters
本文開展了冰形特征參數預測研究,建立了神經網絡預測模型,并利用算例進行了預測分析,得到以下結論:
1)本文提供了一種冰形特征參數快速預測方法,以結冰條件作為網絡輸入,冰形特征參數作為輸出,可以實現不同結冰條件下的冰形幾何特性的快速預測,簡單且有效。
2)對明冰、霾冰、混合冰而言,冰形特征參數的冰角角度、結冰范圍(結冰上下極限)和上冰角高度的預測結果與數值模擬計算結果相接近,平均絕對百分比誤差MAPE在5%以下;而下冰角高度預測結果不佳,后續需要對其進行單獨預測研究。
需要注意的是,該方法預測的范圍易受到訓練樣本的氣象和飛行條件輸入向量的限制,僅在一定范圍內才能保證預測結果可靠,因此后續的研究需要獲取更多的冰形數據,加入結冰試驗的數據進一步改善網絡,提高預測精度。
致謝:感謝中國空氣動力研究與發展中心結冰與防/除冰重點實驗室易賢研究員、王強和劉宇工程師,感謝他們在翼型結冰研究方面提供的幫助和文章撰寫方面提出的建設性意見。