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

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

0 引言

飛機(jī)在大氣環(huán)境中飛行時(shí)，其迎風(fēng)表面可能出現(xiàn)冰層積聚，這一現(xiàn)象被稱為飛機(jī)結(jié)冰現(xiàn)象。發(fā)生在飛機(jī)機(jī)翼的結(jié)冰會(huì)顯著改變飛機(jī)的氣動(dòng)外形，從而導(dǎo)致飛機(jī)氣動(dòng)性能惡化、操縱性能下降、穩(wěn)定性降低。升力表面的積冰會(huì)導(dǎo)致附近的空氣流場(chǎng)提前發(fā)生轉(zhuǎn)捩，增大飛機(jī)的摩擦和壓差阻力。飛機(jī)機(jī)翼表面結(jié)冰嚴(yán)重影響了飛行安全，飛機(jī)結(jié)冰防護(hù)研究對(duì)于飛機(jī)安全性具有重大意義。

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

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

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

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

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

1 電熱除冰計(jì)算模型與方法

電熱除冰是一個(gè)傳熱、傳質(zhì)耦合的物理過程，涉及外部水滴流動(dòng)、表面溢流與相變以及材料導(dǎo)熱等過程，其物理現(xiàn)象如圖1 所示。在外部對(duì)流換熱、水滴撞擊以及結(jié)構(gòu)加熱等因素的影響下，表面發(fā)生溢流及相變（凍結(jié)、融化、蒸發(fā)或升華）現(xiàn)象。由于加熱功率的周期性及加熱組件的空間布置，電熱除冰期間的溢流相變呈現(xiàn)非穩(wěn)態(tài)特性與空間差異性。下面從水滴流場(chǎng)計(jì)算、表面溢流相變計(jì)算和導(dǎo)熱計(jì)算進(jìn)行介紹。

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

1.1 水滴流場(chǎng)

空氣流場(chǎng)采用雷諾平均Navier-Stokes 方程（RANS）進(jìn)行描述，其形式如下：

式中雷諾應(yīng)力項(xiàng)由Boussinesq 定義：

式中，ρ為空氣密度，ui、uj為時(shí)均速度分量，p為壓強(qiáng)，μ為動(dòng)力黏度，δij為 克羅內(nèi)克張量分量，u′i、u′j為脈動(dòng)速度分量，μt為湍動(dòng)黏度，k為湍動(dòng)能。

水滴流場(chǎng)的連續(xù)性方程和動(dòng)量方程由下式給出：

式中，ρw為水滴的密度，α為水滴的容積分?jǐn)?shù)，u為水滴的速度矢量，ua為空氣的速度矢量，F(xiàn)是作用在水滴上的除阻力外的外力，K是空氣-水滴動(dòng)量交換系數(shù)，由下式定義：

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

1.2 表面溢流與結(jié)冰相變

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

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

質(zhì)量守恒方程：

能量守恒方程：

因此根據(jù)表面相態(tài)可分為4 種情況考慮：

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

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

3）0 ＜f＜ 1，控制體液態(tài)水結(jié)冰，控制體處于冰水混合相。有：

4）f=1，控制體內(nèi)全部結(jié)冰。有：

式（11）與式（12）在形式上相同，但對(duì)應(yīng)熱流項(xiàng)有不同的含義。以上公式計(jì)算細(xì)節(jié)可查閱參考文獻(xiàn)[20]。

1.3 導(dǎo)熱計(jì)算模型

電熱除冰過程中蒙皮的導(dǎo)熱是非穩(wěn)態(tài)導(dǎo)熱過程，采用基于非結(jié)構(gòu)網(wǎng)格的數(shù)值計(jì)算方法。因?yàn)樵诒締栴}中，蒙皮溫度變化范圍不大，蒙皮材料物性認(rèn)為是常物性的，其非穩(wěn)態(tài)導(dǎo)熱過程可用如下微分方程描述：

式中，ρ為材料的密度，c為材料的比熱容，λ為材料的導(dǎo)熱系數(shù)，S為熱源項(xiàng)。

如圖3 所示的控制體單元，基于有限體積法，用隱式格式對(duì)中心網(wǎng)格單元的導(dǎo)熱方程進(jìn)行離散化處理：

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

圖3 導(dǎo)熱計(jì)算的控制體單元Fig.3 Control volume for the heat conduction calculation

整理成迭代形式：

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

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

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

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

1.4 計(jì)算方法

具體計(jì)算流程圖見圖5?？紤]到本文重點(diǎn)在非穩(wěn)態(tài)電熱除冰，因此水滴流場(chǎng)的計(jì)算在此不作贅述，其結(jié)果如剪切力、對(duì)流換熱系數(shù)和水收集系數(shù)等數(shù)據(jù)作為輸入條件供溢流相變-結(jié)構(gòu)導(dǎo)熱計(jì)算使用，由于外流場(chǎng)網(wǎng)格與除冰蒙皮結(jié)構(gòu)網(wǎng)格在交接界面處存在差異，利用距離反比插值方法將外部流場(chǎng)結(jié)果插值到蒙皮網(wǎng)格上。

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

在電熱除冰過程中，空氣、水滴兩相流流場(chǎng)、表面溢流相變與固態(tài)導(dǎo)熱存在耦合關(guān)系。其中空氣流場(chǎng)、水滴流場(chǎng)的計(jì)算結(jié)果和蒙皮結(jié)構(gòu)的導(dǎo)熱計(jì)算結(jié)果之間的耦合，通過表面溢流相變計(jì)算完成，因此需要一個(gè)合理的耦合方法完成復(fù)雜的耦合計(jì)算。

考慮到電熱除冰過程中，冰層厚度會(huì)在一個(gè)較小的范圍內(nèi)，認(rèn)為結(jié)冰不會(huì)影響到外部空氣流場(chǎng)和水滴流場(chǎng)，為了提高計(jì)算效率，本文忽略了除冰過程中外部流場(chǎng)和水滴流場(chǎng)的變化?；诜€(wěn)態(tài)的外流場(chǎng)計(jì)算結(jié)果，利用歐拉方法計(jì)算水滴流場(chǎng)，得到的剪切力和水收集系數(shù)等結(jié)果用于溢流水相變模型，并通過溫度邊界弱耦合的方式將導(dǎo)熱和溢流相變計(jì)算耦合起來。

本文在風(fēng)雷框架下開發(fā)表面溢流相變求解器和導(dǎo)熱求解器，結(jié)構(gòu)導(dǎo)熱與溢流相變過程采用在邊界處耦合的方法進(jìn)行迭代計(jì)算。在一個(gè)時(shí)間步長(zhǎng)內(nèi)，依次求解非穩(wěn)態(tài)導(dǎo)熱方程與表面溢流相變方程，更新各自的邊界條件繼續(xù)迭代，直至計(jì)算結(jié)果滿足收斂判據(jù)（兩次迭代溫度殘差ΔT＜ 0.001℃）。在當(dāng)前時(shí)間步長(zhǎng)內(nèi)計(jì)算收斂后，進(jìn)行時(shí)間推進(jìn)，計(jì)算下一個(gè)時(shí)間步長(zhǎng)。耦合求解的具體步驟如下：

1）首先根據(jù)外流場(chǎng)剪切力計(jì)算溢流相變過程中的溢流水控制體計(jì)算順序，并將其存儲(chǔ)于數(shù)組中；

2）根據(jù)環(huán)境溫度對(duì)蒙皮進(jìn)行溫度的初始化；

3）根據(jù)電加熱控制率進(jìn)行加熱熱源的加載；

4）計(jì)算蒙皮導(dǎo)熱，得到蒙皮溫度分布，通過溫度梯度計(jì)算外表面導(dǎo)熱熱流密度qcond；

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

6）更新蒙皮外表面溫度，重復(fù)步驟4）和5），進(jìn)行一個(gè)時(shí)間步長(zhǎng)內(nèi)的迭代計(jì)算，直到前后兩次迭代中，所有溢流水控制體溫度變化量ΔTs滿足收斂要求；

7）更新結(jié)冰量等表面狀態(tài)，推進(jìn)到下一時(shí)間步長(zhǎng)，重復(fù)上述步驟3）至6），直至推進(jìn)到規(guī)定的計(jì)算時(shí)間，計(jì)算結(jié)束。

2 模型驗(yàn)證

2.1 導(dǎo)熱計(jì)算的驗(yàn)證

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

圖6 導(dǎo)熱驗(yàn)證蒙皮模型Fig.6 Skin model for the thermal conductivity verification

圖7 導(dǎo)熱驗(yàn)證蒙皮網(wǎng)格Fig.7 Skin mesh for the thermal conductivity verification

對(duì)自編三維非穩(wěn)態(tài)導(dǎo)熱求解器進(jìn)行了網(wǎng)格無關(guān)性（網(wǎng)格量4 300～ 6 780）和時(shí)間步長(zhǎng)（1～ 0.5 s）的無關(guān)性驗(yàn)證，50 s 時(shí)刻外表面網(wǎng)格的溫度結(jié)果見圖8，不同網(wǎng)格量和時(shí)間步長(zhǎng)下的導(dǎo)熱計(jì)算結(jié)果一致，說明導(dǎo)熱計(jì)算結(jié)果與網(wǎng)格數(shù)量、時(shí)間步長(zhǎng)無關(guān)?；陲L(fēng)雷自編程序的導(dǎo)熱計(jì)算結(jié)果與Fluent 計(jì)算結(jié)果的誤差控制在0.05%以內(nèi)。誤差產(chǎn)生的原因是本文導(dǎo)熱計(jì)算采用了一階精度的離散格式，對(duì)于該模型中部區(qū)域，網(wǎng)格不是嚴(yán)格正交，在輸運(yùn)方程中，體現(xiàn)在面法向和網(wǎng)格質(zhì)心連線不平行，該誤差可通過調(diào)整網(wǎng)格劃分來控制。

圖8 導(dǎo)熱計(jì)算結(jié)果及對(duì)比Fig.8 Heat conduction calculation results and comparison

2.2 電熱除冰計(jì)算的驗(yàn)證

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

圖9 驗(yàn)證模型示意圖Fig.9 Schematic diagram of the validation model

表1 實(shí)驗(yàn)環(huán)境條件Table 1 Experimental conditions

表2 實(shí)驗(yàn)材料物性參數(shù)Table 2 Physical parameters of the experimental materials

圖10 NASA 實(shí)驗(yàn)加熱規(guī)律（功率單位：W/m2）Fig.10 Heating law of the NASA experiment（power in W/m2）

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

3 電熱除冰計(jì)算與分析

3.1 幾何模型及計(jì)算條件

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

表3 計(jì)算條件Table 3 Calculation conditions

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

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

3.2 基準(zhǔn)條件下的計(jì)算結(jié)果

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

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

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

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

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

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

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

3.3 改進(jìn)方案及結(jié)果分析

對(duì)于該工況下的電熱除冰情況進(jìn)行分析可知，撞擊區(qū)的積冰在加熱融化后，部分未蒸發(fā)的液態(tài)水會(huì)向后緣溢流，至未加熱區(qū)，因此需要合理配置各加熱區(qū)位置、電加熱控制率，使各加熱區(qū)按照一定規(guī)則加熱，以保證電熱除冰的效果。在原加熱方案中，對(duì)于撞擊區(qū)對(duì)應(yīng)的2 和3 區(qū)域，配置加熱時(shí)間過長(zhǎng)，加熱功率過大，而對(duì)于溢流區(qū)域加熱區(qū)1 和5 的加熱時(shí)間不足，導(dǎo)致加熱區(qū)5 溢流水重復(fù)結(jié)冰。

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

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

表5 改進(jìn)的電加熱控制率Table 5 Improved control law for the electric heating

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

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

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

圖20 結(jié)冰厚度隨時(shí)間變化曲線Fig.20 Variation curve of the icing thickness with time

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

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

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

引入每個(gè)周期內(nèi)的能量消耗Ecycle來定量描述除冰系統(tǒng)的性能指標(biāo)。其計(jì)算公式如下：

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

4 結(jié)論

本文在電熱除冰模型和非穩(wěn)態(tài)導(dǎo)熱模型研究的基礎(chǔ)上，將表面溢流相變和蒙皮導(dǎo)熱耦合，基于風(fēng)雷平臺(tái)，構(gòu)建了電熱除冰計(jì)算功能模塊。導(dǎo)熱計(jì)算結(jié)果與Fluent 計(jì)算結(jié)果對(duì)比發(fā)現(xiàn)，導(dǎo)熱的表面溫度誤差控制在0.1%以內(nèi)；電熱除冰仿真計(jì)算與實(shí)驗(yàn)數(shù)據(jù)對(duì)比發(fā)現(xiàn)，表面溫度變化顯著相似，誤差控制在3%以內(nèi)。說明基于風(fēng)雷軟件的非穩(wěn)態(tài)電熱除冰二次開發(fā)能夠很好地進(jìn)行電熱除冰耦合計(jì)算。在計(jì)算中發(fā)現(xiàn)，不合理的加熱片布局和加熱控制率會(huì)導(dǎo)致大量溢流水向后方溢流，在溢流區(qū)結(jié)冰后形成瘤狀冰，不能達(dá)到預(yù)期的除冰效果，需結(jié)合具體工況條件下的電熱除冰情況，通過重新布置加熱片、調(diào)整加熱控制率等方式調(diào)整電熱除冰方案。經(jīng)過方案修正，除冰效果得到了顯著改善，溢流區(qū)沒有瘤狀冰形成，蒙皮表面溫度大幅降低，系統(tǒng)能耗降低20%以上。