短艙防冰系統三維內外流耦合計算方法

鄭梅，馮麗娟，秦娜，尹金鴿

中國航發商用航空發動機有限責任公司，上海 200241

當飛機穿越由過冷水滴組成的云層時，迎風部件表面會發生結冰現象［1-3］。發動機短艙進氣道表面產生的冰聚集會改變發動機進氣流場和進氣流量，可能導致發動機推力下降、喘振或失速［4］；當表面積冰發生破碎時，冰片可能被吸入發動機內部，撞擊風扇葉片，造成發動機機械損傷，嚴重危害飛行安全［5-7］。目前，在役的民用航空發動機普遍采用熱氣防冰結構，其中笛形管熱氣防冰系統由于其結構簡單、可靠性強、技術成熟度高等優點，被廣泛應用在發動機短艙進氣道上［8-9］。在短艙笛形管熱氣防冰系統設計過程中，防冰性能仿真分析是設計流程的一個重要環節，構建有效的防冰計算方法可以大幅提升設計效率，為防冰性能評估及防冰結構優化設計提供支撐。

笛形管熱氣防冰系統的性能仿真涉及復雜的內外流耦合傳熱［10-13］?；诜辣到y的運行環境與工作原理，防冰腔內、外流動換熱特性差異顯著［14-18］，為仿真計算帶來較大難度。文獻中關于笛形管防冰系統的耦合仿真大多針對機翼防冰腔構型［19-23］，所采用的耦合計算方法可分為強固耦合方法［24-25］和松散耦合方法［26-27］。強固耦合計算通常將防冰腔外流域、固體域及內流域進行一體化網格劃分，在計算過程中同時求解內外流場及固體域導熱［28-29］；而松散耦合方法則通常將計算域進行拆分并選取合適的迭代策略以實現內外流域的耦合計算。由于松散耦合并不要求內外表面網格節點完全一致，因此在減少計算網格和提升計算效率上具有較大優勢［30］。對于短艙進氣道笛形管防冰系統，由于全環模型尺寸大、笛形管上射流孔數多，且進氣道型面具有周向非對稱特性，無法參考機翼模型截取一段直翼防冰腔進行模擬。綜合計算網格與計算效率等因素，松散耦合方法更為適用。松散耦合迭代過程中涉及內外流計算域耦合面的數據交互。相較于干空氣條件，濕空氣條件下防冰表面熱質傳遞過程更為復雜，且進氣道表面還需考慮周向非對稱流場導致的三維溢流效應對防冰表面熱載荷的影響。

基于此，針對短艙全尺寸模型的笛形管熱氣防冰系統，采用松散耦合方法開展了三維防冰內外流仿真計算方法研究。將三維短艙進氣道計算模型中的外流場獨立為一個計算域，將蒙皮固體域與防冰腔內流場組合為另一個計算域，以防冰腔蒙皮外表面作為兩個計算域的耦合面，通過耦合面溫度和換熱系數（Heat Transfer Coeffi?cient， HTC）相互賦值為邊界條件來實現內外流計算域的數據交互。根據防冰計算流程，依次構建干、濕空氣條件下的防冰內外流松散耦合迭代策略，并在濕空氣條件下引入短艙周向非均勻性對防冰表面溢流流動換熱特性的影響。最終獲得迭代收斂后的短艙進氣道防冰腔內外表面溫度及換熱系數分布。

1 三維防冰內外流耦合計算方法

1.1 防冰計算流程

通常，防冰計算流程包括以下3個步驟：流場計算、水滴撞擊特性計算和防冰熱分析計算。

1.1.1 流場計算

采用ANSYS Fluent商用軟件進行流場計算，通過求解三維、定常、黏性流動的Navier-Stokes （N-S）方程，獲得防冰腔外部冷空氣流場、防冰腔內部熱空氣流場以及蒙皮固壁溫度場。此時，為干空氣條件下的流場計算，防冰腔內、外表面對流換熱系數均可通過CFD方法計算得到。其中，防冰腔外表面對流換熱系數以外流場環境溫度為參考溫度換算得到；而防冰腔內表面對流換熱系數則以熱氣總溫作為參考溫度進行提取。

1.1.2 水滴撞擊特性計算

水滴撞擊特性計算采用歐拉法［31］，引入水滴相體積分數，通過求解水滴相控制方程獲得水滴運動軌跡及防冰腔外表面局部水收集系數分布。

水滴相的連續性方程和動量方程可描述為

式中：α為水滴相體積分數；ρd為水滴密度，由于水滴是不可壓縮的，因此式（1）中水滴密度項可約去；ud為水滴速度；CD為阻力系數；ρa為空氣密度；dd為水滴直徑；ua為空氣速度。

根據定義，表面局部水收集系數β的計算公式為

式中：n為表面撞擊位置的單位法向量；α∞與u∞分別為遠場區域位置水滴相體積分數和氣流速度。

將水滴相控制方程展開為Fluent用戶自定義標量方程形式，采用Fluent自帶的求解器，實現水滴相x、y、z方向速度與水滴相體積分數的迭代求解。

1.1.3 防冰熱分析計算

防冰熱分析計算采用改進的Messinger模型，考慮過冷水滴撞擊防冰腔外表面后未完全蒸發而形成的溢流流動影響，但暫不考慮防冰不足導致溢流冰形成的情況。此時防冰腔外表面的熱質傳遞如圖1［32］所示。

對于防冰腔外表面任意控制體積，在穩態情況下均滿足質量守恒與能量守恒，即

圖1 考慮溢流水流動的防冰腔外表面熱質傳遞［32］Fig.1 Heat and mass transfer on external surface of anti-icing chamber considering runback flow［32］

式中：為溢流流入質量流量；為過冷水滴撞擊質量流量；為溢流流出質量流量；為蒸發質量 流量；為防冰系 統外表 面加熱量；為水滴撞擊帶入的總焓值，包括水滴撞擊動能及水滴自身所攜帶的焓；為溢流流入帶入的滯止 焓；為表面 對 流 換 熱量；為 表 面蒸發帶走的總焓值；為溢流流出帶走的滯止焓。以上各項具體表達式詳見文獻［32］。

由于防冰表面溢流流動主要受氣動剪切力與表面壓力梯度驅動，在短艙進氣道周向非對稱流場的作用下，會呈現較強的三維特性。如圖2所示，圖中CVi為控制體積，Ai為面法向量，為流入當前控制體積的質量流量，為流出當前控制體積的質量流量，i=1,2,…,4，j=1,2，防冰表面任意控制體積與其相鄰的上下游控制體積之間可能存在多個溢流流入流出面，每個面上的流入流出質量流量可根據質量守恒并結合氣流流動特性計算得到。

圖2 三維防冰表面溢流分配示意圖Fig.2 Schematic of runback distribution on 3D antiicing surface

通過Fluent用戶自定義函數（User-Defined Functions ，UDFs），將水滴撞擊、溢流流動及表面蒸發等能量源項引入Fluent計算的能量方程中，根據能量守恒在給定的防冰腔外表面溫度邊界條件下求解表面加熱量Q?anti?icing，wet。

以外流環境溫度為參考溫度，根據防冰腔外表面加熱量換算得到濕空氣條件下防冰腔外表面等效換熱系數：

式中：Tw，wet為濕空氣條件下防冰腔外表面溫度；T∞為外流場環境溫度；A為換熱面積。

考慮固體域導熱，穩態情況下防冰腔內外表面達到熱平衡：

式中：為防冰腔內表面對流換熱量；為固體域導熱量。濕空氣條件下防冰腔內表面對流換熱系數計算方法與干空氣條件下保持一致，僅防冰腔外表面熱載荷發生變化，由此可迭代得到濕空氣條件下防冰腔固體域溫度場。

1.2 內外流松散耦合迭代策略

1.2.1 干空氣流場計算

圖3 干空氣流場計算的內外流松散耦合迭代策略Fig.3 Loose-coupling iterative strategy of internal and external flow for flow calculation under dry air condition

干空氣條件下內外流松散耦合迭代策略如圖3所示，每一迭代輪次中均包含外流計算、內流與固體域計算。除第1輪次（Run 1）為初場計算外，其余輪次（Run 2~N）內均先開展內流與固體域計算后再進行外流計算。同一迭代輪次內，內外流耦合的數據交互為：內流與固體域計算獲得的防冰腔外表面溫度分布賦值給外流計算作為第2類邊界條件。相鄰輪次耦合迭代計算的數據傳遞為：前一輪次（RunN?1）外流計算獲得的防冰腔外表面對流換熱系數分布作為后一輪次（RunN）內流與固體域計算的第3類邊界條件。計算過程中，不斷更新防冰腔外表面對流換熱系數分布與溫度分布，當相鄰迭代輪次的表面平均溫度偏差不超過±0.5 ℃時，可判定迭代收斂。

1.2.2 濕空氣防冰熱分析計算

圖4 濕空氣流場計算的內外流松散耦合迭代策略Fig.4 Loose-coupling iterative strategy of internal and external flow for flow calculation under wet air condition

濕空氣條件下，內外流松散耦合迭代策略如圖4所示。以干空氣條件下內外流耦合計算的收斂結果作為初場，引入水滴撞擊對防冰表面換熱特性的影響，外流流場計算的同時還需開展防冰表面熱分析計算。同一迭代輪次中，先開展外流防冰熱分析計算，提取防冰腔外表面等效換熱系數作為第3類邊界條件賦值給內流與固體域計算。相鄰迭代輪次中，傳遞的數據則為內流與固體域計算獲得的防冰腔外表面溫度分布。在反復多次耦合迭代后，不斷更新防冰腔外表面參數，迭代收斂依據仍為相鄰迭代輪次防冰腔外表面平均溫度偏差不超過±0.5 ℃。

2 計算模型與網格劃分

2.1 計算模型與計算域

短艙笛形管熱氣防冰系統三維全尺寸模型如圖5所示。將短艙所處外流域獨立為一個計算域，防冰腔固體域和內流域合并為另一個計算域。以防冰腔蒙皮外表面作為兩個計算域的耦合面，通過耦合面數據交互實現內外流域的耦合迭代。根據計算域拆分方式，依次完成外流計算域和內流計算域的網格劃分。

圖5 計算模型與內外流計算域Fig.5 Computational model and computational do?mains of internal and external flow

2.2 外流計算域網格與邊界條件

采用ANSYS ICEM軟件對計算域進行網格劃分。外流計算域網格如圖6所示，總網格數約444萬。為了準確模擬附面層內的復雜流動特性，短艙及其延長段的近壁面網格均進行了局部加密處理，附面層第1層網格高度約為0.01 mm，保證壁面y+≈1，網格尺度與選用的剪切應力輸運（Shear Stress Transfer， SST）κ-ω湍流模型相匹配；另外，在進氣道唇口附近等流場參數變化較為劇烈的區域也適當進行了網格加密處理。

圖6 外流計算域網格Fig.6 Mesh for external computational domain

外流計算域邊界條件設置如下：外流邊界為壓力遠場；風扇進口為壓力出口；進氣道防冰腔外表面為內外流松散耦合迭代時的數據交互面，盡管該表面邊界類型始終為壁面邊界條件，但賦值的溫度分布會隨迭代輪次不斷更新；除上述邊界以外，其余均為絕熱壁面。

外流計算選用密度基求解器，流域中空氣物性參數選用理想氣體模型進行計算。松散耦合迭代過程中每一輪次的外流計算均同時監控流場殘差、防冰腔外表面平均換熱系數及風扇進口流量，計算收斂后流場結果才可用于后續迭代。計算的收斂標準為各殘差指標下降3個數量級或殘差曲線不再下降，并且外表面平均換熱系數和風扇進口流量保持穩定。

2.3 內流計算域網格與邊界條件

內流與固體域組成的計算網格如圖7所示。為精確模擬熱氣由射流孔噴射沖擊至防冰腔內表面的過程，需對每個小孔的射流方向及其沿程周向進行網格局部加密。同時，由于防冰腔內表面射流沖擊流動換熱特性較為復雜，為提高計算精度，需對內表面近壁面網格也進行局部加密，第1層網格高度約為0.01 mm，確保壁面y+與SSTκ-ω湍流模型相匹配。三維全環模型的內流與固體域計算網格總數約為1 822萬。

圖7 內流計算域網格Fig.7 Mesh for internal computational domain

內流與固體域計算邊界條件如下：防冰腔和前隔框的內表面均為流固耦合面；前隔框外表面為絕熱壁面；笛形管管壁為定壁溫邊界條件；管壁上均布的射流孔為壓力入口邊界條件；防冰腔外表面仍為內外流松散耦合計算的數據交互面，壁面邊界條件賦值的換熱系數也會隨耦合迭代過程不斷更新。

內流與固體域計算時，同樣選用密度基求解器，流域中熱空氣物性參數選用理想氣體模型進行計算。松散耦合迭代過程中每一輪次內流與固體域計算均同時監控流場殘差和防冰腔外表面平均溫度。計算的收斂標準同樣要求流場殘差與監控的表面參數收斂曲線保持不變。

2.4 計算工況

計算工況如表1所示，包括飛行狀態參數、結冰氣象條件以及射流孔出流參數。

表1 計算工況Table 1 Computational conditions

3 計算結果與分析

3.1 算例驗證

截取短艙進氣道12點鐘方向局部模型，依次開展強固耦合計算與松散耦合計算，通過對比分析以驗證本文松散耦合計算方法及耦合迭代策略的有效性。圖8為驗證模型的幾何及計算網格。

圖8 計算方法驗證模型及網格Fig.8 Model and mesh for computational method validation

圖9 干空氣條件下強固耦合與松散耦合計算結果對比 （z=0）Fig.9 Comparison between tight-coupling and loosecoupling methods under dry air condition （z=0）

圖10 濕空氣條件下強固耦合與松散耦合計算結果對比（z=0）Fig.10 Comparison between tight-coupling and loosecoupling methods under wet air condition （z=0）

圖9和圖10分別為干、濕空氣條件下強固耦合計算方法與松散耦合計算方法的對比結果，其中h0為對流換熱系數基準值，T0為溫度基準值，圖9所示干空氣條件下松散耦合計算結果為迭代3個輪次后的收斂結果，圖10所示濕空氣條件下計算結果為迭代4個輪次后的收斂結果?？梢钥闯鰞煞N計算方法所得到的防冰腔內外表面各參數基本保持一致。由此可知，所采用的松散耦合計算方法及迭代策略可以達到與強固耦合計算方法基本一致的仿真精度。

3.2 干空氣流場

圖11 干空氣條件下短艙進氣道防冰腔外表面對流換熱系數分布Fig.11 Convective heat transfer coefficient distributions on external surface of nacelle inlet anti-icing chamber under dry air condition

圖12 干空氣條件下短艙進氣道防冰腔外表面溫度分布Fig.12 Temperature distributions on external surface of nacelle inlet anti-icing chamber under dry air condition

圖13 干空氣條件下短艙進氣道防冰腔內表面對流換熱系數分布Fig.13 Convective heat transfer coefficient distributions on internal surface of nacelle inlet anti-icing chamber under dry air condition

圖11~圖13為干空氣條件下短艙進氣道笛形管防冰腔外表面對流換熱系數、溫度以及內表面對流換熱系數在內外流松散耦合不同迭代輪次下的計算結果。在經過3個輪次的內外流耦合迭代計算后，防冰腔內外表面各參數分布逐漸趨于收斂。其中，由于干空氣條件下第1輪次為初場計算，防冰腔外表面溫度設置為定壁溫，因此圖12中并未給出該輪次下的溫度分布，但增加了第4輪次內流與固體域計算結果。對比耦合迭代的第3輪次和第4輪次防冰腔外表面溫度分布可以看出，兩者基本達到一致。

圖14 干空氣條件下短艙進氣道防冰腔內外表面參數分布曲線（12點鐘方向）Fig.14 Parameter distribution curves on internal and ex?ternal surfaces of nacelle inlet anti-icing chamber under dry air condition （12 o’clock direction）

如圖14所示，取短艙進氣道防冰腔在12點鐘方向上的參數分布曲線，可以更加準確地分析干空氣條件下內外流松散耦合迭代計算的收斂過程?？梢钥闯觯?輪次防冰腔內外表面對流換熱系數分布均與第2輪次計算結果有較為明顯的差距，但第2輪次和第3輪次各參數分布曲線基本重合；第3、4輪次的表面溫度曲線吻合良好，可達到迭代收斂要求；防冰腔內表面對流換熱系數分布曲線在迭代過程中變化不大，僅射流駐點區域隨迭代輪次略有下降。綜上，干空氣條件下防冰腔內、外表面平均對流換熱系數在耦合迭代過程中均呈現先減后增再逼近收斂結果的變化規律，而防冰腔外表面平均溫度則先增后減至迭代收斂。

3.3 水滴撞擊特性

圖15 短艙進氣道防冰腔外表面局部水收集系數分布Fig.15 Local collection efficiency distributions on ex?ternal surface of nacelle inlet anti-icing chamber

圖15為短艙進氣道防冰腔外表面局部水收集系數分布，圖中β0為局部水收集系數基準值?？梢钥闯?，水滴撞擊區域主要集中在唇口前緣駐點位置附近，且唇口內表面撞擊極限要明顯大于唇口外表面，沿唇口周向各截面上局部水收集系數分布并不均勻，12點鐘方向具有最大的局部水收集系數峰值和撞擊范圍，而6點鐘方向撞擊范圍明顯減小，3點鐘和9點鐘方向左右對稱，具有相同的局部水收集系數分布。由于水滴撞擊特性計算結果是后續防冰熱分析計算的輸入，唇口表面局部水收集系數的周向分布不均勻性很大程度上會影響濕空氣條件下防冰腔表面的換熱特性及溫度分布。

3.4 濕空氣防冰熱分析

圖16~圖18為濕空氣條件下短艙進氣道笛形管防冰腔外表面等效換熱系數、溫度以及內表面對流換熱系數在內外流松散耦合不同迭代輪次下的計算結果。相比于干空氣條件，濕空氣條件下耦合迭代輪次明顯增加，這是由于該工況條件下表面撞擊水量不能完全蒸發并在防冰腔表面形成周向不均勻分布的溢流水流動，如圖19所示，其中m0為溢流水質量流量基準值。溢流水量分布與表面換熱特性相互影響，導致表面溫度分布收斂速度減緩。防冰腔外表面由于過冷水滴撞擊引入的熱質傳遞源項，濕空氣條件下外表面等效換熱系數要明顯大于干空氣條件；而對于防冰腔內表面對流換熱系數，干、濕空氣條件下數值變化并不大。此外，相較于干空氣條件，濕空氣條件下進氣道防冰腔外表面溫度明顯下降，且表面溫度均勻性有所提升。

圖17 濕空氣條件下短艙進氣道防冰腔外表面溫度分布Fig.17 Temperature distributions on external surface of nacelle inlet anti-icing chamber under wet air condition

圖18 濕空氣條件下短艙進氣道防冰腔內表面對流換熱系數分布Fig.18 Convective heat transfer coefficient distributions on internal surface of nacelle inlet anti-icing chamber under wet air condition

圖19 濕空氣條件下短艙進氣道防冰腔外表面溢流水質量流量分布（第8輪次迭代）Fig.19 Runback mass flux on external surface of na?celle inlet anti-icing chamber under wet air con?dition （Run 8）

圖20 濕空氣條件下短艙進氣道防冰腔內外表面參數分布曲線（12點鐘方向）Fig.20 Parameter distributions curves on internal and external surfaces of nacelle inlet anti-icing chamber under wet air condition （12 o’clock direction）

結合圖20所示的短艙進氣道防冰腔內外表面各參數在12點鐘方向上的分布曲線可以看出，在撞擊區域內松散耦合迭代的第1輪次到第4輪次防冰腔外表面溫度下降明顯，且外表面等效換熱系數分布也發生較大的變化，但從第5輪次開始到第8輪次計算結果都較為接近，尤其在唇口內表面各參數曲線吻合良好，但在唇口外表面撞擊極限之后的蒙皮外表面等效換熱系數及溫度分布還略有差異?？偟膩碚f，濕空氣條件下防冰腔內外表面參數在耦合迭代過程中所呈現的變化規律與干空氣條件恰好相反。防冰腔外表面平均等效換熱系數與內表面平均對流換熱系數在迭代過程中均先增加后減小再逼近收斂，而防冰腔外表面平均溫度則先降低后升高直至迭代收斂。

4 結 論

針對三維全尺寸短艙進氣道笛形管熱氣防冰系統，開展了內外流耦合計算方法研究，結合進氣道防冰計算流程，確定了干、濕空氣條件下短艙防冰系統內外流松散耦合迭代策略及計算收斂所需的迭代輪次，得到的主要結論如下：

1）干、濕空氣條件下內外流松散耦合迭代均以防冰腔外表面作為內、外計算域的耦合面，并通過耦合面溫度和換熱系數相互賦值為邊界條件來實現內外計算域的數據交互，其中濕空氣條件下考慮過冷水滴撞擊引入的熱質傳遞影響及三維溢流效應，采用防冰腔外表面等效換熱系數作為耦合面邊界條件。

2）干空氣條件下短艙進氣道防冰腔內外流松散耦合僅需3個輪次的迭代計算即可收斂。

3）濕空氣條件下進氣道防冰腔內表面在松散耦合迭代4個輪次后趨于收斂，但受表面溢流流動換熱影響，防冰腔外表面溫度分布在耦合迭代過程中仍有一定波動。后續將針對該問題，對內外流耦合迭代策略作進一步優化。