基于PHengLEI 的非穩態電熱除冰過程仿真

劉宗輝，卜雪琴，*，林貴平，李偉斌

(1.北京航空航天大學 航空科學與工程學院，北京 100083；2.中國空氣動力研究與發展中心，綿陽 621000)

0 引言

飛機在大氣環境中飛行時，其迎風表面可能出現冰層積聚，這一現象被稱為飛機結冰現象。發生在飛機機翼的結冰會顯著改變飛機的氣動外形，從而導致飛機氣動性能惡化、操縱性能下降、穩定性降低。升力表面的積冰會導致附近的空氣流場提前發生轉捩，增大飛機的摩擦和壓差阻力。飛機機翼表面結冰嚴重影響了飛行安全，飛機結冰防護研究對于飛機安全性具有重大意義。

目前，國內已有多種防除冰辦法，包括機械除冰、化學除冰、熱除冰等。其中周期性電加熱除冰系統是在水滴撞擊區和溢流區安裝電加熱片，采用電加熱的方式，使得機翼表面被周期性加熱：空氣中的過冷水滴撞擊到飛機機翼蒙皮表面時迅速結冰，通過電熱除冰系統對機翼表面進行加熱，在加熱過程中，表面溫度不斷升高，冰層開始融化；當系統不加熱時，表面溫度持續降低，可能繼續結冰；如此結冰-融冰周而復始，控制結冰惡化，達到結冰防護效果。因此合理安排電加熱片的位置分布，制定合適的電加熱控制率，對達到良好的除冰效果有重要作用。

國外較早對電熱除冰系統開展了數值計算和實驗研究。Chauvin 等[1]研究了結冰過程中機翼蒙皮的熱傳導，發現機翼熱傳導是不可忽略的。Fossati、Habashi等[2]提出了水滴流場計算的降階模型，可以快速得到水收集系數的低階近似，在此基礎上，Pourbagian、Habashi 等[3]針對電熱除冰系統，將降階模型（ROM）用于除冰耦合計算，進一步降低了除冰系統計算量，并基于能源和安全因素，提出了不同的目標函數作為除冰系統優化準則，提升了加熱控制率設計效率。Botura 等[4]提出脈沖電熱除冰系統，通過冰風洞實驗發現該系統在減少能量消耗和冰脊的形成上取得了良好的效果。同時，國外研究機構還開發了電熱除冰仿真軟件或模塊，例如加拿大NTI 公司開發了Fensap-Ice 軟件，具備電熱除冰計算模塊[5]。Bennani 等[6]的ONERA 軟件包含除冰計算功能，采用有限元方法分析了二維電熱除冰的過程，用破碎擴展模型來模擬冰層的脫落。

國內對于電熱除冰過程的研究內容大多集中在解決型號設計中面臨的實際問題。2005 年，艾劍波等[7]提出了針對直升機旋翼的電熱除冰系統。2007 年，常士楠[8-10]采用了焓方法分析了二維融冰過程。2010 年，卜雪琴等[11]提出了表面無結冰情況下的防冰模型及數值解法。肖春華等[12]采用多孔介質的方法，探討了加熱規律因素對表面溫度的影響，使用有限元方法計算了冰層內應力分布。2013 年，申曉斌等[13]基于Messinger 結冰模型，根據外流場空氣-水膜的剪切力邊界條件，建立了三維表面溢流水動力學模型，為三維電熱除冰模擬奠定了基礎。2014 年，傅見平等[14]基于Messinger 模型與改進的焓方法分析了二維電熱除冰過程，發現合理設計加熱片分布、加熱時間以及熱流密度能夠實現電熱除冰系統能源的高效利用。2016 年，楊詩雨等[15]研究了旋轉帽罩電加熱防冰過程，通過仿真分析，發現同樣的電加熱功率下，周期電加熱防冰更為節能。2018 年，穆作棟等[16]研究了三維電熱除冰模型，并對除冰過程中的溢流結冰、相變過程、固壁導熱進行了耦合計算。2020 年申曉斌等[17]利用NACA0012 翼型進行了非穩態電熱除冰模擬，得到了不同時刻表面溫度、溢流水和結冰厚度等結果。

以上學者的研究都說明機翼蒙皮電熱除冰效果除了受到飛行條件和氣象參數的影響外，還由電熱除冰系統的加熱功率、加熱控制率和電加熱片分布決定。在電熱除冰的優化設計上，尤其是電熱片控制率和分布的優化上有待研究。已有的防除冰算法和模型基本都是以國外成熟的CFD 軟件為基礎進行計算的，國內對于基于自主可控的CFD 軟件的防除冰功能開發有較大空白。風雷(PHengLEI)軟件作為國家數值風洞工程的大型通用CFD 計算軟件[18-19]，該軟件集成了高精度結構化網格、高超聲速計算、大規模并行計算等模塊，基于PHengLEI 進行電熱除冰計算功能的開發，有助于自主可控的國產化軟件工程的進一步推進。

本文以電熱除冰過程中的結冰-融冰耦合算法為基礎，開展了三維機翼模型的非穩態電熱除冰過程的數值計算方法研究，并對算法進行了驗證；基于國產PHengLEI 軟件平臺，進行了電熱除冰功能的開發，說明了基于PHengLEI 平臺開發的電熱除冰算法的可行性和正確性。在此基礎上，開展了不同電加熱控制率下的結冰-融冰耦合計算，結果說明不合理的加熱控制率不僅浪費了系統能源、導致過高的蒙皮溫度，對蒙皮材料產生影響，還可能在溢流區形成冰瘤，嚴重影響飛行安全，并據此提出了電熱除冰系統控制率的改進方案。本文研究，是國內首次將電熱除冰計算功能集成在我國自主可控CFD 軟件PHengLEI 上，對我國CFD 軟件進一步拓展多學科、多物理模型集成開發奠定基礎。

1 電熱除冰計算模型與方法

電熱除冰是一個傳熱、傳質耦合的物理過程，涉及外部水滴流動、表面溢流與相變以及材料導熱等過程，其物理現象如圖1 所示。在外部對流換熱、水滴撞擊以及結構加熱等因素的影響下，表面發生溢流及相變（凍結、融化、蒸發或升華）現象。由于加熱功率的周期性及加熱組件的空間布置，電熱除冰期間的溢流相變呈現非穩態特性與空間差異性。下面從水滴流場計算、表面溢流相變計算和導熱計算進行介紹。

圖1 電熱除冰涉及到的物理現象Fig.1 Physical phenomena involved in electrothermal deicing

1.1 水滴流場

空氣流場采用雷諾平均Navier-Stokes 方程（RANS）進行描述，其形式如下：

式中雷諾應力項由Boussinesq 定義：

式中，ρ為空氣密度，ui、uj為時均速度分量，p為壓強，μ為動力黏度，δij為 克羅內克張量分量，u′i、u′j為脈動速度分量，μt為湍動黏度，k為湍動能。

水滴流場的連續性方程和動量方程由下式給出：

式中，ρw為水滴的密度，α為水滴的容積分數，u為水滴的速度矢量，ua為空氣的速度矢量，F是作用在水滴上的除阻力外的外力，K是空氣-水滴動量交換系數，由下式定義：

其中，dp為水滴直徑大小，fa是阻力函數。

1.2 表面溢流與結冰相變

在除冰系統工作時，受到電加熱的作用，一般表面溫度高于273.15 K，撞擊的水升溫并少量蒸發，沒有蒸發的液態水沿著氣流方向溢流，在表面形成溢流水，溢流水流動方向由空氣流場剪切力方向決定；停止加熱階段，當溫度低于或等于273.15 K 的表面，則存在溢流結冰。建立如圖2 所示的表面控制體的質量和能量守恒方程，如下所述。

圖2 表面控制容積的質量與能量守恒示意圖Fig.2 Schematic diagram of the mass and energy conservation of the surface control volume

質量守恒方程：

能量守恒方程：

因此根據表面相態可分為4 種情況考慮：

1）f=0 且表面無積冰，即mice=0，控制體內全部處于液相。有：

2）f=0 且表面存在積冰，控制體積冰融化，處于冰水混合相。有：

3）0 ＜f＜ 1，控制體液態水結冰，控制體處于冰水混合相。有：

4）f=1，控制體內全部結冰。有：

式（11）與式（12）在形式上相同，但對應熱流項有不同的含義。以上公式計算細節可查閱參考文獻[20]。

1.3 導熱計算模型

電熱除冰過程中蒙皮的導熱是非穩態導熱過程，采用基于非結構網格的數值計算方法。因為在本問題中，蒙皮溫度變化范圍不大，蒙皮材料物性認為是常物性的，其非穩態導熱過程可用如下微分方程描述：

式中，ρ為材料的密度，c為材料的比熱容，λ為材料的導熱系數，S為熱源項。

如圖3 所示的控制體單元，基于有限體積法，用隱式格式對中心網格單元的導熱方程進行離散化處理：

其中：V為控制體體積；Tp1、Te1、Tw1、Tn1、Ts1、Tf1、Tb1為 控制體t+1 時刻的溫度，是未知量；Tp0是該控制體t時刻（上一時刻）溫度。Ae、Aw、An、As、Af、Ab為相鄰控制體接觸面積；Δde、Δdw、Δdn、Δds、Δdf、Δdb為相鄰控制體質心距離；其幾何關系示意圖見圖3。

圖3 導熱計算的控制體單元Fig.3 Control volume for the heat conduction calculation

整理成迭代形式：

蒙皮導熱計算時，蒙皮的外表面邊界采用溫度邊界條件（Dirichlet boundary condition），由表面水膜溢流相變的計算結果給出。相應的，導熱過程計算得到的外表面熱流結果，為水膜溢流相變的求解提供除冰加熱熱流，熱流由固體邊界的溫度梯度確定，如式（18）：

式中n為固體邊界的外法向量。具體水膜和蒙皮導熱計算的邊界條件交替如圖4。

圖4 邊界條件示意圖Fig.4 Schematic diagram of the boundary conditions

在實際物理過程中，加熱片非常薄，為了更接近物理過程，蒙皮的加熱面采用熱流邊界條件（Neumann boundary condition），由加熱片加熱功率直接給出。

1.4 計算方法

具體計算流程圖見圖5。考慮到本文重點在非穩態電熱除冰，因此水滴流場的計算在此不作贅述，其結果如剪切力、對流換熱系數和水收集系數等數據作為輸入條件供溢流相變-結構導熱計算使用，由于外流場網格與除冰蒙皮結構網格在交接界面處存在差異，利用距離反比插值方法將外部流場結果插值到蒙皮網格上。

圖5 電熱除冰的求解流程圖Fig.5 Flow chart of the electrothermal deicing solver

在電熱除冰過程中，空氣、水滴兩相流流場、表面溢流相變與固態導熱存在耦合關系。其中空氣流場、水滴流場的計算結果和蒙皮結構的導熱計算結果之間的耦合，通過表面溢流相變計算完成，因此需要一個合理的耦合方法完成復雜的耦合計算。

考慮到電熱除冰過程中，冰層厚度會在一個較小的范圍內，認為結冰不會影響到外部空氣流場和水滴流場，為了提高計算效率，本文忽略了除冰過程中外部流場和水滴流場的變化。基于穩態的外流場計算結果，利用歐拉方法計算水滴流場，得到的剪切力和水收集系數等結果用于溢流水相變模型，并通過溫度邊界弱耦合的方式將導熱和溢流相變計算耦合起來。

本文在風雷框架下開發表面溢流相變求解器和導熱求解器，結構導熱與溢流相變過程采用在邊界處耦合的方法進行迭代計算。在一個時間步長內，依次求解非穩態導熱方程與表面溢流相變方程，更新各自的邊界條件繼續迭代，直至計算結果滿足收斂判據（兩次迭代溫度殘差ΔT＜ 0.001℃）。在當前時間步長內計算收斂后，進行時間推進，計算下一個時間步長。耦合求解的具體步驟如下：

1）首先根據外流場剪切力計算溢流相變過程中的溢流水控制體計算順序，并將其存儲于數組中；

2）根據環境溫度對蒙皮進行溫度的初始化；

3）根據電加熱控制率進行加熱熱源的加載；

4）計算蒙皮導熱，得到蒙皮溫度分布，通過溫度梯度計算外表面導熱熱流密度qcond；

5）假定溢流水控制單元溫度為273.15 K，按照步驟1）中得到的控制體計算順序，求解溢流水溢流相變模型，得到所有溢流水控制體結冰速率或融冰速率，同時得到新的表面溫度Tsnew；

6）更新蒙皮外表面溫度，重復步驟4）和5），進行一個時間步長內的迭代計算，直到前后兩次迭代中，所有溢流水控制體溫度變化量ΔTs滿足收斂要求；

7）更新結冰量等表面狀態，推進到下一時間步長，重復上述步驟3）至6），直至推進到規定的計算時間，計算結束。

2 模型驗證

2.1 導熱計算的驗證

為了驗證本文導熱模型和求解器的準確性，將在風雷框架下自編導熱求解器的計算結果與Fluent 軟件計算結果進行了對比驗證。導熱計算驗證選用某機翼的一部分，其幾何模型如圖6。機翼蒙皮厚8 mm，密度為7 800 kg/m3，導熱系數為16.7 W/(m·℃)，比熱容為502 J/(kg·℃)。蒙皮采用六面體結構化網格，網格如圖7 所示，網格總數為4 300 個。計算初始溫度為263 K，設置外側蒙皮為熱流邊界條件，其熱流密度為4 000 W/m2，其余面全部為絕熱邊界，計算時間為50 s。

圖6 導熱驗證蒙皮模型Fig.6 Skin model for the thermal conductivity verification

圖7 導熱驗證蒙皮網格Fig.7 Skin mesh for the thermal conductivity verification

對自編三維非穩態導熱求解器進行了網格無關性（網格量4 300～ 6 780）和時間步長（1～ 0.5 s）的無關性驗證，50 s 時刻外表面網格的溫度結果見圖8，不同網格量和時間步長下的導熱計算結果一致，說明導熱計算結果與網格數量、時間步長無關。基于風雷自編程序的導熱計算結果與Fluent 計算結果的誤差控制在0.05%以內。誤差產生的原因是本文導熱計算采用了一階精度的離散格式，對于該模型中部區域，網格不是嚴格正交，在輸運方程中，體現在面法向和網格質心連線不平行，該誤差可通過調整網格劃分來控制。

圖8 導熱計算結果及對比Fig.8 Heat conduction calculation results and comparison

2.2 電熱除冰計算的驗證

公開文獻中的電熱除冰實驗記載甚少，選取Al-Khalil 等利用NASA Lewis 冰風洞（IRT）開展的非穩態電熱除冰實驗[21]，以驗證本文模型及算法；該實驗結果同樣被FENSAP 軟件作為驗證參考。計算模型如圖9 所示，其實驗環境條件和材料參數見表1 和表2，加熱控制率見圖10，其中序號1～7 分別對應圖9 中的不同位置加熱片編號。圖11 給出了加熱區1溫度的模擬結果與文獻測量結果對比，數據吻合良好，溫度變化曲線具有顯著相似性，溫度誤差在5 K以內，說明本文模型和算法的正確性。

圖9 驗證模型示意圖Fig.9 Schematic diagram of the validation model

表1 實驗環境條件Table 1 Experimental conditions

表2 實驗材料物性參數Table 2 Physical parameters of the experimental materials

圖10 NASA 實驗加熱規律（功率單位：W/m2）Fig.10 Heating law of the NASA experiment（power in W/m2）

圖11 加熱區1 溫度曲線Fig.11 Temperature curve of heating zone 1

3 電熱除冰計算與分析

3.1 幾何模型及計算條件

機翼蒙皮模型如圖6 所示，蒙皮厚8 mm，展向長20.9 cm，密度為7 800 kg/m3，導熱系數為16.7 W/(m·℃)，比熱容為502 J/(kg·℃)。加熱組件采用5 個獨立的電加熱片，布置在蒙皮內表面。加熱分區示意圖見圖12，計算條件如表3。電加熱控制率周期為12 s，其電加熱初步控制率見表4。

表3 計算條件Table 3 Calculation conditions

表4 電加熱初步控制率Table 4 Preliminary control law of the electric heating

圖12 電加熱分區示意圖Fig.12 Schematic diagram of the electric heating zone

3.2 基準條件下的計算結果

圖13、圖14 是一個加熱周期內，各區域蒙皮表面溫度隨時間變化云圖和曲線圖，從圖中可以看出，在該加熱控制率下，整個外蒙皮溫升較大。最高溫度達到了350 K，比273 K 高出了約70 K，對蒙皮材料的耐溫性能造成了一定影響。其中在整個控制率中加熱時間較長，加熱功率較大的Heater3 和Heater4 區域，整體溫度較高，而加熱時間短的Heater1 和Heater5 區域溫度較低。

圖13 除冰周期部分時刻的溫度分布云圖Fig.13 Temperature contours at different time instances of the deicing cycle

圖14 加熱周期各區域溫度曲線Fig.14 Temperature curve of each heater during the heating cycle

圖15、圖16 分別是加熱周期內幾個時刻的溢流水量、冰層厚度云圖和結冰厚度曲線圖。從整個溢流水計算表面的結冰量看來，計算開始時，結冰區域主要集中在1 區和2 區，1 區、2 區的加熱很快將積冰融化，形成溢流水區域，隨著加熱的進行，撞擊區上部的積冰也不斷融化，約4 s 時完全融化，融化的積冰分別往上游和下游溢流。溢流到加熱區4 和5 的水，受限于加熱時間，溢流水在對流換熱和與蒙皮的熱傳導過程中，溫度很快降低，重新結冰，在y坐標0.03 m位置處形成冰瘤。

圖15 除冰周期各時刻的溫度分布云圖Fig.15 Temperature contours at different time instances of the deicing cycle

圖16 結冰厚度變化曲線Fig.16 Ice thickness variation curve

由此可見，該加熱方案很快將撞擊區積冰融化，融化后的水向后溢流，溢流到下端的水在加熱區外重復結冰，溢流到上端的水在上端溢流區重復結冰，但沒能融化整個防冰區域的積冰和重復結冰的溢流水。另外，由于加熱區2、3、4 配置的加熱時間過長，造成區域溫度過高。過高的區域溫度不僅不滿足蒙皮材料性能要求，而且造成了系統能源的浪費。

3.3 改進方案及結果分析

對于該工況下的電熱除冰情況進行分析可知，撞擊區的積冰在加熱融化后，部分未蒸發的液態水會向后緣溢流，至未加熱區，因此需要合理配置各加熱區位置、電加熱控制率，使各加熱區按照一定規則加熱，以保證電熱除冰的效果。在原加熱方案中，對于撞擊區對應的2 和3 區域，配置加熱時間過長，加熱功率過大，而對于溢流區域加熱區1 和5 的加熱時間不足，導致加熱區5 溢流水重復結冰。

基于上述計算和分析結果，針對該工況條件，對加熱片位置和加熱功率進行改進，如圖17 和表5。由前文計算結果可知，撞擊區會反復結冰，為了避免撞擊區結冰和撞擊區溫度過高的問題，此區域改為低功率恒加熱區域。同時采用逐級加熱的辦法，在更低的能耗下，實現融冰效果。

圖17 加熱片的新布局Fig.17 New layout of the heating zone

表5 改進的電加熱控制率Table 5 Improved control law for the electric heating

對于該工況下的電熱除冰，做了網格無關性驗證（網格量6 780～22 302）和時間步長無關性驗證（1～0.5 s），如圖18 是計算時間為2 s 時刻的表面溫度曲線，說明計算結果與網格數量和時間步長無關。圖19 是加熱片局部溫度隨時間變化的曲線圖，可以看出在改進后的加熱片布局和控制率下，整體溫度明顯降低，且運行8 s 之后，全局溫度高于273.15 K，整個除冰區域已經沒有積冰。機翼結冰量隨時間變化曲線如圖20 所示。可見改進后的機翼結冰情況得到了很大的改善，撞擊區由于恒加熱，溫度一直高于結冰溫度，Heater1 和Heater6 區域在0～6 s 有一些結冰，該區域在開啟加熱后，冰層很快融化。Heater5 區域在6 s 后關閉加熱后，12 s 時刻溫度再次降低到結冰溫度，出現了結冰，在下個加熱周期開始時，冰層融化。

圖18 2 s 時刻表面溫度曲線Fig.18 Surface temperature curves at 2 s

圖19 優化后加熱周期各區域溫度曲線Fig.19 Temperature curve of each region in the optimized heating cycle

圖20 結冰厚度隨時間變化曲線Fig.20 Variation curve of the icing thickness with time

在修改后的加熱區位置和加熱控制率作用下，第11 s 和第15 s 的溢流水量和結冰厚度云圖如圖21、圖22 所示。結合不同時刻的結冰厚度曲線圖，可以看出冰層在加熱片作用下很快融化，融化后的溢流主要往下端溢流，溢流到上端的水再次結冰，隨后在相應區域加熱作用下融化，完成該工況的除冰過程。對比之前的方案，現方案沒有產生溢流水重復結冰形成的瘤狀冰，蒙皮溫度有較大的下降，系統整體能耗大幅下降，針對該工況條件，改進后方案的除冰性能明顯改善。

圖21 11 s 時刻溢流水量和結冰厚度云圖Fig.21 Contours of the overflow water volume and the icing thickness at 11 s

圖22 15 s 時刻溢流水量和結冰厚度云圖Fig.22 Contours of the overflow water volume and the icing thickness at 15 s

引入每個周期內的能量消耗Ecycle來定量描述除冰系統的性能指標。其計算公式如下：

計算得到修改前的除冰方案：Ecycle=0.109 J，修改后的方案：Ecycle=0.087 J，改進后的除冰方案系統整體能耗下降20.1%。

4 結論

本文在電熱除冰模型和非穩態導熱模型研究的基礎上，將表面溢流相變和蒙皮導熱耦合，基于風雷平臺，構建了電熱除冰計算功能模塊。導熱計算結果與Fluent 計算結果對比發現，導熱的表面溫度誤差控制在0.1%以內；電熱除冰仿真計算與實驗數據對比發現，表面溫度變化顯著相似，誤差控制在3%以內。說明基于風雷軟件的非穩態電熱除冰二次開發能夠很好地進行電熱除冰耦合計算。在計算中發現，不合理的加熱片布局和加熱控制率會導致大量溢流水向后方溢流，在溢流區結冰后形成瘤狀冰，不能達到預期的除冰效果，需結合具體工況條件下的電熱除冰情況，通過重新布置加熱片、調整加熱控制率等方式調整電熱除冰方案。經過方案修正，除冰效果得到了顯著改善，溢流區沒有瘤狀冰形成，蒙皮表面溫度大幅降低，系統能耗降低20%以上。