刘宗辉，卜雪琴，*，林贵平，李伟斌
(1.北京航空航天大学 航空科学与工程学院，北京 100083；
2.中国空气动力研究与发展中心，绵阳 621000)
飞机在大气环境中飞行时，其迎风表面可能出现冰层积聚，这一现象被称为飞机结冰现象。发生在飞机机翼的结冰会显著改变飞机的气动外形，从而导致飞机气动性能恶化、操纵性能下降、稳定性降低。升力表面的积冰会导致附近的空气流场提前发生转捩，增大飞机的摩擦和压差阻力。飞机机翼表面结冰严重影响了飞行安全，飞机结冰防护研究对于飞机安全性具有重大意义。
目前，国内已有多种防除冰办法，包括机械除冰、化学除冰、热除冰等。其中周期性电加热除冰系统是在水滴撞击区和溢流区安装电加热片，采用电加热的方式，使得机翼表面被周期性加热：空气中的过冷水滴撞击到飞机机翼蒙皮表面时迅速结冰，通过电热除冰系统对机翼表面进行加热，在加热过程中，表面温度不断升高，冰层开始融化；
当系统不加热时，表面温度持续降低，可能继续结冰；
如此结冰-融冰周而复始，控制结冰恶化，达到结冰防护效果。因此合理安排电加热片的位置分布，制定合适的电加热控制率，对达到良好的除冰效果有重要作用。
国外较早对电热除冰系统开展了数值计算和实验研究。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 电热除冰涉及到的物理现象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.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.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%。
本文在电热除冰模型和非稳态导热模型研究的基础上，将表面溢流相变和蒙皮导热耦合，基于风雷平台，构建了电热除冰计算功能模块。导热计算结果与Fluent 计算结果对比发现，导热的表面温度误差控制在0.1%以内；
电热除冰仿真计算与实验数据对比发现，表面温度变化显著相似，误差控制在3%以内。说明基于风雷软件的非稳态电热除冰二次开发能够很好地进行电热除冰耦合计算。在计算中发现，不合理的加热片布局和加热控制率会导致大量溢流水向后方溢流，在溢流区结冰后形成瘤状冰，不能达到预期的除冰效果，需结合具体工况条件下的电热除冰情况，通过重新布置加热片、调整加热控制率等方式调整电热除冰方案。经过方案修正，除冰效果得到了显著改善，溢流区没有瘤状冰形成，蒙皮表面温度大幅降低，系统能耗降低20%以上。