曾飛雄，毛漢冬，章儒宸

（上海飛機(jī)設(shè)計(jì)研究院，上海 201210）

短艙防冰是指發(fā)動(dòng)機(jī)進(jìn)氣道唇口前緣的防冰方式，設(shè)計(jì)指標(biāo)包括干蒸發(fā)和濕蒸發(fā)。短艙防冰的方式主要分為電防冰、熱氣防冰。兩種防護(hù)方式在仿真建模時(shí)，主要在內(nèi)部熱流分析時(shí)存在差異。

國(guó)外從20 個(gè)世紀(jì)四五十年代就開(kāi)始采用數(shù)值計(jì)算的方法研究飛機(jī)結(jié)冰防護(hù)問(wèn)題。目前國(guó)際上通用的防冰傳熱耦合分析的基礎(chǔ)理論均基于Messinger 模型［1］，并以此來(lái)建立防冰表面的傳熱傳質(zhì)方程。至今數(shù)個(gè)發(fā)達(dá)國(guó)家都開(kāi)發(fā)了用于飛機(jī)結(jié)冰預(yù)測(cè)和防冰系統(tǒng)設(shè)計(jì)的計(jì)算軟件，如NASA和波音的LEWICE 軟件，空客公司的ONERA 軟件，龐巴迪公司的CANICE 軟件，羅羅公司的ICECREMO 軟件，ANSYS 公司的FENSAP?ICE軟件等。其中不少軟件被SAE?ARP5903 收錄為成熟軟件進(jìn)行推薦［2］。

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

1 性能分析關(guān)鍵指標(biāo)

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

唇口內(nèi)部冰脫落的影響通常大于外部冰脫落，因此主要性能考核指標(biāo)側(cè)重于進(jìn)氣道內(nèi)表面的結(jié)冰情況分析。由于冰型尺寸通常難以在試驗(yàn)中進(jìn)行復(fù)現(xiàn)驗(yàn)證，因此進(jìn)行短艙防冰性能分析時(shí)，通常也把表面溫度作為間接考核的判據(jù)，以此表明計(jì)算模型的準(zhǔn)確性。

進(jìn)氣道內(nèi)表面吞冰量通常轉(zhuǎn)換成某站位的結(jié)冰厚度進(jìn)行判定。通過(guò)仿真分析穩(wěn)態(tài)后流水量及持續(xù)時(shí)間得到溢流冰量，再參照FAR33.77［10］定義的發(fā)動(dòng)機(jī)名義結(jié)冰尺寸，得到某站位的結(jié)冰厚度。

短艙防冰的熱平衡與外界環(huán)境、發(fā)動(dòng)機(jī)抽吸量、飛行狀態(tài)、引氣狀態(tài)等相關(guān)。因此性能分析時(shí)，需要考量不同飛行工況作為計(jì)算狀態(tài)點(diǎn)。 通常均需考慮單發(fā)起飛、正常起飛、爬升、巡航、下降、待機(jī)、進(jìn)場(chǎng)以及著陸等飛行狀態(tài)。

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

2 總體計(jì)算模型

短艙防冰的性能計(jì)算模型重點(diǎn)是建立內(nèi)部熱流與外部熱流的平衡關(guān)系。計(jì)算時(shí)需要考慮4 個(gè)計(jì)算區(qū)域，包括外部流場(chǎng)區(qū)域（空氣?水滴兩相流）、表面水膜區(qū)域、蒙皮結(jié)構(gòu)區(qū)域以及防冰腔內(nèi)部流場(chǎng)區(qū)域。為了加快計(jì)算速度，同時(shí)保證足夠可靠的計(jì)算精度，本文熱流耦合計(jì)算模型以獨(dú)立的蒙皮網(wǎng)格為計(jì)算域，調(diào)用包含全部網(wǎng)格的完整內(nèi)外流場(chǎng)結(jié)果作為初始值，在后續(xù)平衡迭代計(jì)算中對(duì)外流場(chǎng)中離蒙皮一定距離以外的局部熱流場(chǎng)進(jìn)行更新，對(duì)于內(nèi)部流場(chǎng)進(jìn)行少量的更新。整個(gè)計(jì)算流程如圖1所示。

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

能量平衡方程主要考慮外部的導(dǎo)熱項(xiàng)、水滴顯熱項(xiàng)、水滴潛熱項(xiàng)、對(duì)流換熱項(xiàng)、氣動(dòng)加熱項(xiàng)和內(nèi)部的對(duì)流換熱項(xiàng)［1］。質(zhì)量平衡主要考慮每個(gè)微元體中各項(xiàng)水量的平衡，如式（1）所示。

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

在進(jìn)行計(jì)算域網(wǎng)格繪制時(shí)，常使用半模機(jī)體繪制外流場(chǎng)網(wǎng)格。需要將機(jī)身、機(jī)翼、吊掛、發(fā)動(dòng)機(jī)包含在內(nèi)。針對(duì)短艙防冰計(jì)算分析的特點(diǎn)，需要特別地對(duì)短艙區(qū)域進(jìn)行網(wǎng)格加密，并將唇口耦合面單獨(dú)繪制網(wǎng)格，便于后續(xù)進(jìn)行內(nèi)外流場(chǎng)耦合時(shí)確定計(jì)算域。考慮到發(fā)動(dòng)機(jī)的抽吸效應(yīng)，需要在發(fā)動(dòng)機(jī)建模時(shí)單獨(dú)設(shè)置壓力出口和壓力入口邊界，用于模擬發(fā)動(dòng)機(jī)的抽吸氣量（圖2）。

圖2 發(fā)動(dòng)機(jī)邊界設(shè)置圖Fig.2 Engine boundary condition setting

防冰腔體的內(nèi)流場(chǎng)根據(jù)防冰形式主要分為基于笛形管的小孔射流和基于直流噴嘴的環(huán)狀流動(dòng)。噴孔處需要映射加密以提高內(nèi)流場(chǎng)的求解精度（圖3）。

圖3 防冰腔內(nèi)流場(chǎng)加密效果示意Fig.3 Inner flow field mesh encryption

3 KD 插值算法

在防冰計(jì)算中涉及外部流場(chǎng)網(wǎng)格、蒙皮結(jié)構(gòu)網(wǎng)格與內(nèi)部流場(chǎng)網(wǎng)格3 套不同的網(wǎng)格，因此熱耦合計(jì)算時(shí)網(wǎng)格交界面間的數(shù)據(jù)需要通過(guò)插值計(jì)算進(jìn)行映射傳遞。在一個(gè)耦合計(jì)算循環(huán)中需要進(jìn)行兩次界面插值計(jì)算，分別為：

（1）蒙皮結(jié)構(gòu)表面網(wǎng)格至內(nèi)外流場(chǎng)壁面邊界網(wǎng)格的熱流和溫度插值。

（2）內(nèi)外流場(chǎng)壁面邊界網(wǎng)格至蒙皮結(jié)構(gòu)表面網(wǎng)格的熱流和溫度插值。

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

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

由于耦合計(jì)算的每一次單個(gè)迭代步內(nèi)都需要進(jìn)行兩次界面插值，計(jì)算至最終穩(wěn)態(tài)結(jié)果需要進(jìn)行多次循環(huán)迭代，并且邊界網(wǎng)格量級(jí)較大，因此，插值是整個(gè)計(jì)算中最耗資源的過(guò)程。另外，內(nèi)外流場(chǎng)網(wǎng)格與蒙皮結(jié)構(gòu)網(wǎng)格可能采用不同的拓?fù)?，為使插值算法具有較好的適應(yīng)性，插值計(jì)算應(yīng)基于網(wǎng)格節(jié)點(diǎn)進(jìn)行，從而使插值計(jì)算過(guò)程與網(wǎng)格單元類型無(wú)關(guān)。為在足夠的插值精度下盡可能提高插值計(jì)算速度，同時(shí)保證插值算法能夠適用于各類不同網(wǎng)格形狀（網(wǎng)格無(wú)關(guān)性），本文采用基于K 近鄰搜索的反距離加權(quán)插值算法對(duì)界面數(shù)據(jù)進(jìn)行插值映射，該算法目前在計(jì)算機(jī)前沿的機(jī)器學(xué)習(xí)及大數(shù)據(jù)處理領(lǐng)域已有應(yīng)用［11?14］。

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

圖5 基于K 近鄰搜索的節(jié)點(diǎn)數(shù)據(jù)反距離加權(quán)插值Fig.5 K-nearest neighbor search for inverse distance weighted interpolation

圖5 中的網(wǎng)格1 為已獲得相應(yīng)數(shù)據(jù)的表面網(wǎng)格，網(wǎng)格2 為待插值獲取數(shù)據(jù)的表面網(wǎng)格，為了得到網(wǎng)格2 中某一節(jié)點(diǎn)的數(shù)值，先對(duì)該節(jié)點(diǎn)周圍N個(gè)距離最近的節(jié)點(diǎn)進(jìn)行查找，然后基于到這N個(gè)最近點(diǎn)的距離進(jìn)行加權(quán)計(jì)算獲得該節(jié)點(diǎn)數(shù)據(jù)。

通常采用的反距離加權(quán)函數(shù)形式為

式中：d為距離，p為冪參數(shù)。當(dāng)p=2 時(shí)，即為以空間距離（歐式距離）的倒數(shù)為權(quán)重?？梢钥闯鼍嚯x越近的臨近點(diǎn)對(duì)樣本點(diǎn)的影響越大。

不同于常見(jiàn)的窮舉搜索方法，K 近鄰搜索算法實(shí)現(xiàn)上采用KD 樹(shù)結(jié)構(gòu)進(jìn)行快速二分匹配，其算法時(shí)間復(fù)雜度為O（lg（2N）），與窮舉搜索時(shí)的O（N2）時(shí)間復(fù)雜度相比，其計(jì)算速度將極大提高，且數(shù)據(jù)量越大時(shí)優(yōu)勢(shì)越明顯。相較于其他插值計(jì)算方法，基于K 近鄰搜索的反距離加權(quán)插值具有兩個(gè)優(yōu)點(diǎn)：

（1）計(jì)算速度快。由于K 近鄰搜索具有較低的時(shí)間復(fù)雜度，且對(duì)于空間位置固定的表面網(wǎng)格節(jié)點(diǎn)拓?fù)?，其KD 樹(shù)只需生成一次，因此在進(jìn)行K 近鄰搜索和逆距離計(jì)算時(shí)計(jì)算開(kāi)銷較小，當(dāng)界面節(jié)點(diǎn)數(shù)大于100 時(shí)，使用K 近鄰插值的計(jì)算速度將比窮舉搜索插值快1 000 倍以上。

（2）適用于各類不同外形與網(wǎng)格間的數(shù)據(jù)映射。插值算法基于網(wǎng)格節(jié)點(diǎn)進(jìn)行，與具體網(wǎng)格類型和拓?fù)錈o(wú)關(guān)，因此具備幾何無(wú)關(guān)性與網(wǎng)格無(wú)關(guān)性，適用于各類不同交界面網(wǎng)格間的數(shù)據(jù)映射。

4 對(duì)流換熱系數(shù)計(jì)算

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

而對(duì)流換熱系數(shù)是數(shù)值求解對(duì)流換熱項(xiàng)的關(guān)鍵參數(shù)。從已有的研究來(lái)看，主要分為3 類求解方式：（1）計(jì)算流體力學(xué)（CFD）數(shù)值求解；（2）附面層積分法求解；（3）基于試驗(yàn)修正的經(jīng)驗(yàn)公式求解。

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

短艙唇口外形與圓管特征類似，考慮到外掠圓管的原理性試驗(yàn)已經(jīng)有充分積累，本文基于外掠圓管的修正經(jīng)驗(yàn)公式［15］，開(kāi)展三維對(duì)流換熱系數(shù)的求解分析。

式中：Nu為努賽爾數(shù)，C為試驗(yàn)修正系數(shù)，Re為雷諾數(shù)，m為雷諾數(shù)修正因子，Pr為普朗特?cái)?shù)。修正因子見(jiàn)表1。

表1 外掠圓管流動(dòng)修正因子Table 1 Flow correction factor of swept tube

式（4）中的相關(guān)無(wú)量綱數(shù)基于外流場(chǎng)計(jì)算結(jié)果中當(dāng)?shù)鼐W(wǎng)格內(nèi)的速度、密度、黏度等參數(shù)進(jìn)行求解，雷諾數(shù)的結(jié)果會(huì)影響到試驗(yàn)修正系數(shù)C和雷諾數(shù)修正因子m，由于不同機(jī)型短艙的前緣外形各有差別，會(huì)對(duì)特征尺寸有影響，因此這兩個(gè)修正系數(shù)也可基于風(fēng)洞或干空氣試飛的溫度結(jié)果進(jìn)行算法修正，極大地提升了結(jié)冰條件下仿真計(jì)算的準(zhǔn)確性和一致性。

5 仿真求解與結(jié)果分析

本文使用的算例是參照某型號(hào)飛機(jī)的自然結(jié)冰試飛數(shù)據(jù)。計(jì)算網(wǎng)格模型基于飛機(jī)真實(shí)構(gòu)型進(jìn)行繪制。由于試飛過(guò)程只能采集唇口區(qū)域蒙皮溫度，而無(wú)法獲取結(jié)冰量等數(shù)據(jù)，因此本仿真求解主要對(duì)比計(jì)算與試驗(yàn)的溫度差異情況。試飛狀態(tài)參數(shù)如表2 所示。

表2 計(jì)算算例工況點(diǎn)Table 2 Simulation cases

空氣外流場(chǎng)計(jì)算的主要邊界條件類型為：計(jì)算區(qū)域的外圍定為壓力遠(yuǎn)場(chǎng)邊界；發(fā)動(dòng)機(jī)進(jìn)氣道入口為指定流量的壓力出口，其壓力值為達(dá)到目標(biāo)流量時(shí)所匹配出的壓力；發(fā)動(dòng)機(jī)尾噴出口為指定流量的壓力入口。

防冰腔內(nèi)流場(chǎng)的主要邊界條件類型為：笛形管噴孔處設(shè)定為壓力入口，通過(guò)一維管路流動(dòng)計(jì)算得到對(duì)應(yīng)的入口總溫和壓力；笛形管表面設(shè)為恒溫壁面；排氣孔出口處設(shè)定為壓力出口，壓力和溫度取為大氣環(huán)境條件。

以發(fā)動(dòng)機(jī)進(jìn)氣道唇口前緣蒙皮為耦合面，通過(guò)插值獲取外流場(chǎng)和內(nèi)流場(chǎng)的相關(guān)參數(shù)，基于蒙皮表面的離散單元進(jìn)行能量和溢流水質(zhì)量平衡求解。圖6、7 分別顯示的是算例1 和2 的唇口表面溫度分區(qū)云圖。從圖中可以看出，有顯著熱斑的區(qū)域集中在內(nèi)表面，高溫區(qū)域出現(xiàn)在180°（下方）區(qū)域，這些均與設(shè)計(jì)目標(biāo)相吻合。

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

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

基于試飛改裝安裝的表面溫度傳感器，對(duì)進(jìn)氣道某截面的溫度數(shù)據(jù)進(jìn)行對(duì)比驗(yàn)證。驗(yàn)證結(jié)果分別如圖8、9 所示。

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

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

仿真結(jié)果表明，截面表面溫度計(jì)算值與試飛實(shí)測(cè)值趨勢(shì)吻合較好，峰值位置較為匹配，偏差較大處出現(xiàn)在外唇口后緣，總體測(cè)點(diǎn)溫度偏差在10 ℃以內(nèi)，在工程領(lǐng)域?qū)儆诳山邮芊秶?/p>

6 結(jié) 論

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

（2）采用KD 樹(shù)插值并使用基于修正經(jīng)驗(yàn)公式獲取了對(duì)流換熱系數(shù)，在保證計(jì)算準(zhǔn)確性的同時(shí)大幅提升求解效率。

（3）對(duì)某型號(hào)的短艙防冰性能進(jìn)行了分析，計(jì)算結(jié)果與試飛結(jié)果吻合較好，滿足工程應(yīng)用需求，具有很好的工程實(shí)踐價(jià)值。