民用飛機短艙防冰性能分析建模與仿真

曾飛雄，毛漢冬，章儒宸

（上海飛機設計研究院，上海 201210）

短艙防冰是指發動機進氣道唇口前緣的防冰方式，設計指標包括干蒸發和濕蒸發。短艙防冰的方式主要分為電防冰、熱氣防冰。兩種防護方式在仿真建模時，主要在內部熱流分析時存在差異。

國外從20 個世紀四五十年代就開始采用數值計算的方法研究飛機結冰防護問題。目前國際上通用的防冰傳熱耦合分析的基礎理論均基于Messinger 模型［1］，并以此來建立防冰表面的傳熱傳質方程。至今數個發達國家都開發了用于飛機結冰預測和防冰系統設計的計算軟件，如NASA和波音的LEWICE 軟件，空客公司的ONERA 軟件，龐巴迪公司的CANICE 軟件，羅羅公司的ICECREMO 軟件，ANSYS 公司的FENSAP?ICE軟件等。其中不少軟件被SAE?ARP5903 收錄為成熟軟件進行推薦［2］。

目前國內針對機翼防冰仿真計算的研究較多，針對短艙防冰相對較少［3?6］。使用較多的商用軟件為FENSAP，也有部分高校使用用戶自定義函數（User defined function，UDF）自編程求解，通常算法是內外流場分別計算，再利用距離反比插值法進行內外流場數據交互耦合［3］，對流換熱系數的計算則通常采用附面層積分法［7?8］或CFD 求解［9］。民機設計的工程應用重視計算效率和用戶自定義，因此本團隊基于基礎軟件工具自研開發了防冰計算軟件，在經典算法上進行優化，從工程化角度出發，結合國內民機的設計研發和完整的適航取證經驗積累，從條款提煉計算判據［10］，優化了對流換熱系數算法，并創新性地使用了KD 樹（K?dimensional tree）數據插值算法，在確保計算精度的同時大幅提升了計算效率。

1 性能分析關鍵指標

針對發動機和飛機對短艙防冰的功能需求，短艙防冰能力性能分析的考核判據主要包括：（1）各飛行工況下的唇口內表面吞冰量均小于發動機可接受吞冰量；（2）各飛行工況下的唇口外表面冰脫落量小于機體可承受的冰積聚尺寸。

唇口內部冰脫落的影響通常大于外部冰脫落，因此主要性能考核指標側重于進氣道內表面的結冰情況分析。由于冰型尺寸通常難以在試驗中進行復現驗證，因此進行短艙防冰性能分析時，通常也把表面溫度作為間接考核的判據，以此表明計算模型的準確性。

進氣道內表面吞冰量通常轉換成某站位的結冰厚度進行判定。通過仿真分析穩態后流水量及持續時間得到溢流冰量，再參照FAR33.77［10］定義的發動機名義結冰尺寸，得到某站位的結冰厚度。

短艙防冰的熱平衡與外界環境、發動機抽吸量、飛行狀態、引氣狀態等相關。因此性能分析時，需要考量不同飛行工況作為計算狀態點。 通常均需考慮單發起飛、正常起飛、爬升、巡航、下降、待機、進場以及著陸等飛行狀態。

計算吞冰量時需要用到的持續時間基于CCAR25 部附錄C 中的定義，通常考慮連續最大結冰條件下17.4 海里（1 海里=1.852 km）的水平穿云時間，間斷最大結冰條件下2.6 海里的水平穿云時間。通常較為嚴酷的工況是連續最大結冰條件下45 min 的待機飛行工況以及下降小推力工況。

2 總體計算模型

短艙防冰的性能計算模型重點是建立內部熱流與外部熱流的平衡關系。計算時需要考慮4 個計算區域，包括外部流場區域（空氣?水滴兩相流）、表面水膜區域、蒙皮結構區域以及防冰腔內部流場區域。為了加快計算速度，同時保證足夠可靠的計算精度，本文熱流耦合計算模型以獨立的蒙皮網格為計算域，調用包含全部網格的完整內外流場結果作為初始值，在后續平衡迭代計算中對外流場中離蒙皮一定距離以外的局部熱流場進行更新，對于內部流場進行少量的更新。整個計算流程如圖1所示。

圖1 防冰計算迭代流程圖Fig.1 Iteration process for anti-ice analysis

能量平衡方程主要考慮外部的導熱項、水滴顯熱項、水滴潛熱項、對流換熱項、氣動加熱項和內部的對流換熱項［1］。質量平衡主要考慮每個微元體中各項水量的平衡，如式（1）所示。

式中：m?imp為撞擊水量，m?in，n為微元體第n個邊的流入水量，m?ice為結冰量，m?evap為蒸發水量，m?out，n為微元體第n個邊的流出水量。防冰計算時，若表面溫度大于0 ℃，可將結冰量設置為0。

在進行計算域網格繪制時，常使用半模機體繪制外流場網格。需要將機身、機翼、吊掛、發動機包含在內。針對短艙防冰計算分析的特點，需要特別地對短艙區域進行網格加密，并將唇口耦合面單獨繪制網格，便于后續進行內外流場耦合時確定計算域。考慮到發動機的抽吸效應，需要在發動機建模時單獨設置壓力出口和壓力入口邊界，用于模擬發動機的抽吸氣量（圖2）。

圖2 發動機邊界設置圖Fig.2 Engine boundary condition setting

防冰腔體的內流場根據防冰形式主要分為基于笛形管的小孔射流和基于直流噴嘴的環狀流動。噴孔處需要映射加密以提高內流場的求解精度（圖3）。

圖3 防冰腔內流場加密效果示意Fig.3 Inner flow field mesh encryption

3 KD 插值算法

在防冰計算中涉及外部流場網格、蒙皮結構網格與內部流場網格3 套不同的網格，因此熱耦合計算時網格交界面間的數據需要通過插值計算進行映射傳遞。在一個耦合計算循環中需要進行兩次界面插值計算，分別為：

（1）蒙皮結構表面網格至內外流場壁面邊界網格的熱流和溫度插值。

（2）內外流場壁面邊界網格至蒙皮結構表面網格的熱流和溫度插值。

兩次插值計算在熱耦合計算流程中進行次序如圖4 所示。

圖4 防冰熱耦合計算流程Fig.4 Anti-ice thermal-coupling process

由于耦合計算的每一次單個迭代步內都需要進行兩次界面插值，計算至最終穩態結果需要進行多次循環迭代，并且邊界網格量級較大，因此，插值是整個計算中最耗資源的過程。另外，內外流場網格與蒙皮結構網格可能采用不同的拓撲，為使插值算法具有較好的適應性，插值計算應基于網格節點進行，從而使插值計算過程與網格單元類型無關。為在足夠的插值精度下盡可能提高插值計算速度，同時保證插值算法能夠適用于各類不同網格形狀（網格無關性），本文采用基于K 近鄰搜索的反距離加權插值算法對界面數據進行插值映射，該算法目前在計算機前沿的機器學習及大數據處理領域已有應用［11?14］。

耦合求解器基于K 近鄰搜索的反距離加權插值方法對界面數據進行插值映射，該算法基本思想如圖5 所示。

圖5 基于K 近鄰搜索的節點數據反距離加權插值Fig.5 K-nearest neighbor search for inverse distance weighted interpolation

圖5 中的網格1 為已獲得相應數據的表面網格，網格2 為待插值獲取數據的表面網格，為了得到網格2 中某一節點的數值，先對該節點周圍N個距離最近的節點進行查找，然后基于到這N個最近點的距離進行加權計算獲得該節點數據。

通常采用的反距離加權函數形式為

式中：d為距離，p為冪參數。當p=2 時，即為以空間距離（歐式距離）的倒數為權重。可以看出距離越近的臨近點對樣本點的影響越大。

不同于常見的窮舉搜索方法，K 近鄰搜索算法實現上采用KD 樹結構進行快速二分匹配，其算法時間復雜度為O（lg（2N）），與窮舉搜索時的O（N2）時間復雜度相比，其計算速度將極大提高，且數據量越大時優勢越明顯。相較于其他插值計算方法，基于K 近鄰搜索的反距離加權插值具有兩個優點：

（1）計算速度快。由于K 近鄰搜索具有較低的時間復雜度，且對于空間位置固定的表面網格節點拓撲，其KD 樹只需生成一次，因此在進行K 近鄰搜索和逆距離計算時計算開銷較小，當界面節點數大于100 時，使用K 近鄰插值的計算速度將比窮舉搜索插值快1 000 倍以上。

（2）適用于各類不同外形與網格間的數據映射。插值算法基于網格節點進行，與具體網格類型和拓撲無關，因此具備幾何無關性與網格無關性，適用于各類不同交界面網格間的數據映射。

4 對流換熱系數計算

對流換熱是影響防冰熱載荷的重要因素，基于工程經驗，外表面對流換熱在某些工況下能占外部總載荷的50%以上。

而對流換熱系數是數值求解對流換熱項的關鍵參數。從已有的研究來看，主要分為3 類求解方式：（1）計算流體力學（CFD）數值求解；（2）附面層積分法求解；（3）基于試驗修正的經驗公式求解。

目前很多商業仿真軟件可以進行復雜流動的CFD 求解，其獲得的是包含了各種換熱因素的總換熱系數，但求解速度通常較慢。附面層積分法普遍用于結冰防冰求解領域，求得邊界層外的速度、溫度分布后，利用相應流態的積分方程得到對流換熱系數的近似解［7］，但在求解積分方程時，通常建立的是沿流線的二維計算模型，對于轉捩和流態的判定準確性仍需要進一步驗證。

短艙唇口外形與圓管特征類似，考慮到外掠圓管的原理性試驗已經有充分積累，本文基于外掠圓管的修正經驗公式［15］，開展三維對流換熱系數的求解分析。

式中：Nu為努賽爾數，C為試驗修正系數，Re為雷諾數，m為雷諾數修正因子，Pr為普朗特數。修正因子見表1。

表1 外掠圓管流動修正因子Table 1 Flow correction factor of swept tube

式（4）中的相關無量綱數基于外流場計算結果中當地網格內的速度、密度、黏度等參數進行求解，雷諾數的結果會影響到試驗修正系數C和雷諾數修正因子m，由于不同機型短艙的前緣外形各有差別，會對特征尺寸有影響，因此這兩個修正系數也可基于風洞或干空氣試飛的溫度結果進行算法修正，極大地提升了結冰條件下仿真計算的準確性和一致性。

5 仿真求解與結果分析

本文使用的算例是參照某型號飛機的自然結冰試飛數據。計算網格模型基于飛機真實構型進行繪制。由于試飛過程只能采集唇口區域蒙皮溫度，而無法獲取結冰量等數據，因此本仿真求解主要對比計算與試驗的溫度差異情況。試飛狀態參數如表2 所示。

表2 計算算例工況點Table 2 Simulation cases

空氣外流場計算的主要邊界條件類型為：計算區域的外圍定為壓力遠場邊界；發動機進氣道入口為指定流量的壓力出口，其壓力值為達到目標流量時所匹配出的壓力；發動機尾噴出口為指定流量的壓力入口。

防冰腔內流場的主要邊界條件類型為：笛形管噴孔處設定為壓力入口，通過一維管路流動計算得到對應的入口總溫和壓力；笛形管表面設為恒溫壁面；排氣孔出口處設定為壓力出口，壓力和溫度取為大氣環境條件。

以發動機進氣道唇口前緣蒙皮為耦合面，通過插值獲取外流場和內流場的相關參數，基于蒙皮表面的離散單元進行能量和溢流水質量平衡求解。圖6、7 分別顯示的是算例1 和2 的唇口表面溫度分區云圖。從圖中可以看出，有顯著熱斑的區域集中在內表面，高溫區域出現在180°（下方）區域，這些均與設計目標相吻合。

圖6 算例1 表面溫度分布云圖Fig.6 Skin temperature contour in Case 1

圖7 算例2 表面溫度分布云圖Fig.7 Skin temperature contour in Case 2

基于試飛改裝安裝的表面溫度傳感器，對進氣道某截面的溫度數據進行對比驗證。驗證結果分別如圖8、9 所示。

圖8 算例1 某截面表面溫度對比圖Fig.8 Section temperature comparison in Case 1

圖9 算例2 某截面表面溫度對比圖Fig.9 Section temperature comparison in Case 2

仿真結果表明，截面表面溫度計算值與試飛實測值趨勢吻合較好，峰值位置較為匹配，偏差較大處出現在外唇口后緣，總體測點溫度偏差在10 ℃以內，在工程領域屬于可接受范圍。

6 結 論

（1）給出了基于工程應用的民用飛機短艙防冰系統性能分析的總體建模與仿真方法。

（2）采用KD 樹插值并使用基于修正經驗公式獲取了對流換熱系數，在保證計算準確性的同時大幅提升求解效率。

（3）對某型號的短艙防冰性能進行了分析，計算結果與試飛結果吻合較好，滿足工程應用需求，具有很好的工程實踐價值。