电加热防/除冰过程中溢流水再冻结数值模拟

刘小暄，朱春玲，朱程香，曹岩; Liu Xiaoxuan; Zhu Chunling; Zhu Chengxiang; Cao Yan

2025年3月13日 23:17 星期四

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

电加热防/除冰过程中溢流水再冻结数值模拟 PDF

- ORCID：
刘小暄 ¹
- ORCID：
朱春玲 ^1,2
✉
- ORCID：
朱程香 ¹
- ORCID：
曹岩 ¹

1. 南京航空航天大学航空学院,南京 210016； 2. 南京航空航天大学航空航天结构力学及控制全国重点实验室,南京 210016

中图分类号： V244.15

最近更新：2024-04-30

DOI：10.16356/j.1005-2615.2024.02.010

摘要

当飞机飞行过程中遇到结冰工况时，为保证飞机的飞行安全以及气动性能，常使用热防护方法作为飞机防除冰的重要手段。本文通过使用共轭耦合方法建立了机翼内部多层导热与水膜流动的耦合计算模型，并针对电加热防除冰过程进行了一系列数值模拟计算。研究选取了可进行非稳态计算的Myers水膜模型，电加热组件导热模型则采用有限体积法隐式离散进行计算，水膜模型与导热模型通过交换边界值的共轭传热方法实现松耦合。研究发现在电加热防除冰过程中，由于积冰融化和水膜流动会产生溢流水再冻结现象。对电加热系统开启后期的积冰外形进行了流场计算与气动分析，发现机翼上下翼面溢流水冻结部位会对流场产生扰动。通过Q准则判断积冰后侧产生了涡结构，对比未发生结冰的干净翼型，机翼前缘位置溢流水冻结导致压力系数曲线发生较大震荡，因此溢流水再冻结问题影响了机翼的气动性能，未能使电加热防除冰系统达到理想的防护效果。

关键词

飞机防/除冰; 电加热防除冰; 耦合传热; 非稳态导热

当飞机在高空中飞行时，云层中含有的微尺度过冷水滴运动至机翼表面并发生撞击，水滴的聚集、流动以及结冰的过程会导致流经飞机表面的气流发生不利分离，进而严重地影响飞机的飞行稳定性以及气动性能。为解决由飞机结冰问题带来的安全隐患，需要使用不同的防/除冰系统以保障飞行安全。

自20世纪以来，国外已经对电加热防/除冰系统开展了一系列实验研究以及数值模拟研究。在电加热除冰的数值模拟方面，Wright等^［

1］使用焓法计算冰的相变，提出了求解二维电加热瞬态导热方程的方法，采用交替隐式方法求解了二维瞬态导热方程以及控制体的质量与能量守恒方程^{［参考文献 2

百度学术}2］，对二维瞬态的导热、冰增长以及冰脱落进行了数值模拟，通过假设节点相位模拟了冰与水的相变，并且将其与早期二维模拟代码LEWICE进行集成。Silva等^{［参考文献 3

百度学术}3］提出了能够预测温度和回流水流向分布的二维防冰数学模型，对水膜流动动量守恒方程以及热传导方程进行边界层积分以求解防冰表面的压力和温度梯度、耦合传热和传质效应^{［参考文献 4

百度学术}4］。Reid等^{［参考文献 5

百度学术}5］对飞行中非稳态电加热防/除冰的数值模拟使用了共轭热传导方法，在每一个时间步长内，对水膜、冰层和固体域之间的热力学平衡方程进行求解，进一步完善了FENSAP‑ICE中的三维计算部分^{［参考文献 6

百度学术}6］，其中的冰增长及水膜流动求解器基于Messinger模型^{［参考文献 7

百度学术}7］进行构造。Song等^{［参考文献 8

百度学术}8］在热力耦合周动力的框架下进行电加热除冰过程的模拟，建立了冰的本构模型并且提出了考虑温度效应的修正破坏准则。Esmaeilifar等^{［参考文献 9

百度学术}9］计算了非稳态的相变以及回流水，此研究中采用共轭传热将积冰与导热求解器紧密耦合，在每个时间步长都加入了子迭代，以确保界面处的温度和热通量收敛。

国内对电加热防除冰问题的数值模拟研究初期围绕工程应用展开，取得了丰富的研究成果。卜雪琴等^［

10］采用附面层积分法对机翼表面传热系数进行计算，将外部防冰载荷、蒙皮弦向导热和内部加热热流三者耦合，采用亚松弛系数稳定迭代过程，得到平衡表面温度、溢流水和溢流结冰的范围。钟国^{［参考文献 11

百度学术}11］根据Messinger模型，采用焓方法对机翼表面积冰的融冰过程进行了数值模拟。郑玉巧等^{［参考文献 12

百度学术}12］采用欧拉法及热传质原理构建叶片结冰过程的数学模型并进行数值模拟分析，并研究了结冰现象导致的气动特性变化。Zhu等^{［参考文献 13

百度学术}13］为优化防/除冰系统对过冷大水滴撞击风力涡轮机叶片进行了数值模拟，提出了针对大水滴改进的欧拉模型。辛苗等^{［参考文献 14

百度学术}14］提出基于Myers水膜流动模型与耦合传热的翼型防/除冰数值模拟方法，利用焓理论及有限体积法进行了复杂多层结构传热的数值计算，并且提出了一种基于焓理论的相变修正方法。Wang等^{［参考文献 15

百度学术}15］提出了超疏水/电加热协同的防除冰策略，并对此进行了冰风洞试验，与普通电加热系统相比降低了27.5%的能耗。刘宗辉等^{［参考文献 16

百度学术}16］则基于PHengLEI平台集成了非稳态电热除冰计算和非稳态导热计算功能，通过温度边界弱耦合的方式将导热和溢流相变计算耦合求解。Shen等^{［参考文献 17

百度学术}17］考虑了飞机结冰过程中撞击水滴的瞬态传热和传质过程，并且采用隐式格式建立非稳态热力学模型，模拟飞机表面水膜和冰层的动态发展。

学者在共轭耦合方法、多层非稳态导热模型等方面针对电加热防/除冰过程进行了大量数值模拟研究，然而仍然存在部分不足。目前的研究大多选择Messinger模型对飞机表面的结冰过程进行模拟，此种数值模拟方法为优化计算忽略控制体内水含量的变化，难以获得水膜内部温度分布等参数，不能准确地模拟水膜运动情况。本文采用的Myers模型由于引入了非定常水膜流动方程，与实际水膜流动的物理过程更加贴近，并且该模型通过计算水膜内部温度分布可以较好地模拟耦合传热过程，在现有使用Myers模型模拟电加热防除冰过程的研究中，学者们主要关注了水膜流动部分，而对溢流水再冻结现象的研究较少。

在飞机电加热防除冰过程后期，由于积冰融化等因素导致水膜流出热防护区域，会产生溢流水再冻结的情况，由于溢流水再冻结常发生在结冰后期，并且形成积冰的位置靠后，因此在电加热防除冰过程中常被忽视。本文研究使用Myers水膜流动模型进行结冰计算，建立了基于共轭传热法的水膜传热耦合模型，通过交换蒙皮处边界值分别对水膜模型以及多层导热模型进行计算。本文中的耦合计算方法对导热过程的描述更贴近实际的物理过程，并且由于Myers模型在非定常水膜计算方面的优势，可以对机翼表面溢流水再冻结问题进行更加深入的研究。

1 物理模型与方法

1.1　水膜流动模型

在电加热防冰过程中，由于电加热系统的开启，大量热量由加热元件传导至机翼表面，因此撞击到机翼表面的水滴将有少部分发生蒸发现象，其余的水滴在升温后会形成沿气流方向流动的水膜。因此为提高防除冰过程数值模拟的准确性，需要充分关注水膜流动的影响，本文采用Myers水膜流动模型^［

18］对电加热防冰过程中的水膜流动进行模拟计算。

在实际情况中，水膜流动由诸多因素影响，并且十分复杂，为了更好地对其进行数值模拟，在考虑空气压力梯度、重力分量以及空气剪切力等的影响下，对水膜流动进行以下简化处理^［

18］：

（1）由于水的动力黏度、比热容、表面张力等物性参数随温度改变较小，在水膜流动计算时取0 ℃时的值。

（2）水膜厚度很小，不考虑水膜沿壁面法线方向的流动速度。

（3）不考虑水膜破碎成溪流或者水珠的情况。

使用 $(η, ξ)$ 表示建立机翼表面的贴体坐标系，其中 $η$ 为与曲面相切方向， $ξ$ 代表垂直方向。根据润滑理论简化水膜流动的Navier‑Stokes控制方程，得到沿水膜流动 $η$ 方向的单位宽度体积流量 $Q_{η}$ ，即水膜的流动速度在高度上积分

Q_{η} = \int_{0}^{h_{w}} V_{η} d ξ = \frac{h_{w}^{2}}{2 μ_{w}} τ_{η} + \frac{h_{w}^{3}}{3 μ_{w}} (ν f_{η} - \frac{\partial p_{a}}{\partial η})

（1）

式中： $V_{η}$ 表示沿 $η$ 方向的水流速度， $h_{w}$ 表示水膜厚度， $μ_{w}$ 表示水的动力黏性系数， $τ_{η}$ 表示气流剪切力沿 $η$ 方向上的分量， $p_{a}$ 为空气压力， $ν f_{η}$ 表示体积力在 $η$ 方向上的分量。在飞机电加热防冰问题中，一般只考虑重力对机翼表面水膜流动的影响。

结冰过程中水膜流动质量平衡方程为

\frac{\partial h_{w}}{\partial t} + \frac{\partial Q_{η}}{\partial η} = \frac{L W C \cdot V_{0} \cdot β}{ρ_{w}} - \frac{m_{e v p}}{ρ_{w}} - \frac{ρ_{i}}{ρ_{w}} \frac{\partial h_{i}}{\partial t}

（2）

$式中 : L W C$ 为液态水含量， $V_{0}$ 为远场来流速度， $β$ 为局部水收集系数， $ρ_{i}$ 和 $ρ_{w}$ 分别为冰和水的密度， $m_{e v p}$ 为单位面积蒸发速率， $h_{i}$ 表示结冰厚度，水膜单位宽度体积流量 $Q_{η}$ 由式（1）计算得出。

在边界条件的处理部分，水膜与机翼界面采用第一类边界条件，由多层导热模型计算得到的机翼表面温度分布作为输入。

1.2　多层导热模型

电加热防冰系统中，机翼内部的电加热组件多为多层材料复合结构，为计算机翼表面以及加热组件内部的温度分布，建立了非稳态多层导热模型。其非稳态导热过程可表述为

ρ_{s} C_{p s} \frac{\partial T}{\partial t} = \nabla \cdot (k_{s} \nabla T) + S (t)

（3）

$式中 : ρ_{s}$ 、 $C_{p s}$ 分别为不同多层导热材料的密度和比热容，T为固体域温度， $k_{s}$ 为不同材料的导热率， $S (t)$ 为随时间变化的热源项。

在边界条件的处理部分，机翼蒙皮表面采用热流边界条件，由水膜流动模型计算得到的边界热流密度作为输入。

蒙皮表面热流示意图如图 1所示，其中 $q_{a}$ 、 $q_{k}$ 、 $q_{c}$ 、 $q_{d}$ 、 $q_{e}$ 分别表示来流气动加热、过冷水滴撞击、对流换热、入射过冷水滴的冷却作用及水膜蒸发带走的热流密度多层导热模型边界条件的热流密度计算公式为

q_{c} = κ_{w} {\frac{\partial T}{\partial z}|}_{z = 0^{+}} = κ_{w} \frac{ε_{0} + ε_{1} {T|}_{z = 0}}{1 - ε_{1} h}

（4）

其中

ε_{0} = \frac{q_{a} + q_{e} + ρ_{L W C} β V_{\infty} c_{p w} T_{d} + α_{c v} T_{\infty}}{κ_{w}}

（5）

ε_{1} = - \frac{α_{c v} + ρ_{L W C} β V_{\infty} c_{p w}}{κ_{w}}

（6）

式中： $q_{c}$ 代表流固域边界的热流密度， $κ_{w}$ 为水导热率， $V_{\infty}$ 和 $T_{\infty}$ 分别代表远场水流速度与温度， $α_{c v}$ 为对流换热系数， $c_{p w}$ 和 $T_{d}$ 分别表示液态水的定压比热容和过冷水滴的温度。

图1 蒙皮表面热流示意图

Fig.1 Schematic diagram of skin surface heat flow

对蒙皮内的多层导热计算区域进行结构化网格的构建，针对控制体单元，使用有限体积法，对中心网格的导热方程进行隐式格式的离散化

\frac{ρ C V}{Δ t} (T_{p}^{1} - T_{p}^{0}) = λ (\frac{T_{e}^{1} - T_{p}^{1}}{Δ d_{e}} A_{e} + \frac{T_{w}^{1} - T_{p}^{1}}{Δ d_{w}} A_{w} + \frac{T_{n}^{1} - T_{p}^{1}}{Δ d_{n}} A_{n} + \frac{T_{s}^{1} - T_{p}^{1}}{Δ d_{s}} A_{e}) + S

（7）

$式中 :$ C为材料比热容，V为控制体体积，下标p代表中心控制体，下标e、w、s、n分别代表东西南北4方向的控制体，0和1则分别代表t时刻以及t+1时刻的状态， $λ$ 为材料导热系数，A为相邻控制体接触面积， $Δ d$ 代表相邻控制体质心间的距离。

将此离散后的控制方程整理为迭代形式

a_{p}^{1} T_{p}^{1} = a_{e}^{1} T_{e}^{1} + a_{w}^{1} T_{w}^{1} + a_{n}^{1} T_{n}^{1} + a_{s}^{1} T_{s}^{1} + a_{p}^{0} T_{p}^{0}

（8）

其中

a_{e}^{1} = \frac{λ A_{e}}{Δ d_{e}} a_{w}^{1} = \frac{λ A_{w}}{Δ d_{w}} a_{n}^{1} = \frac{λ A_{n}}{Δ d_{n}} a_{s}^{1} = \frac{λ A_{s}}{Δ d_{s}}

（9）

a_{p}^{1} = a_{p}^{0} = \frac{ρ C V}{Δ t} a_{p}^{1} = a_{e}^{1} + a_{w}^{1} + a_{n}^{1} + a_{s}^{1} + a_{p}^{0}

（10）

为简便计算，使用中间变量 $a$ 表示不同时间及不同方位下的导热系数与控制体关系。

采用追赶法求解三对角方程组的方法对离散后的微分控制方程进行迭代求解。

1.3　水膜流动与传热耦合方法

为更准确地对飞机电加热防冰过程进行数值模拟，本研究采用共轭传热方法对机翼表面水膜以及蒙皮内部多层加热元件的导热过程进行耦合计算，可以在保证程序稳定性的前提下，有效地提高数值计算的效率。

共轭传热法^［

19］是指流体和固体之间热相互作用的耦合分析，本方法对流体域与固体域分别进行求解，再通过交换界面处边界值，即流体‑固体界面处的温度和热通量分布进行耦合计算。

在耦合计算过程中，采用松耦合方法加速计算收敛，在使用Myers模型对飞机表面水膜进行求解后，将计算获得的边界热流密度作为输入条件代入进多层传热模型进行计算，随后通过传热模型的计算结果更新机翼表面的温度分布，并将其作为水膜流动新的边界条件进行计算。在耦合计算过程中，CHT方法使水膜流动计算过程与导热部分计算可以采用不同时间步长，在一定程度上简便了计算，并提高了稳定性。具体的计算流程如图 2所示。

图2 耦合计算流程图

Fig.2 Flowchart of cupling calculation

本文飞机电加热防冰过程的计算流程如下：

（1）对流场和水滴场进行求解，输出剪切力、对流换热参数、水收集系数等参数。

（2）对计算开始时的水膜厚度、冰层高度、蒙皮表面以及固体域温度进行初始化。

（3）由初始化给定的水膜厚度以及蒙皮表面温度作为水膜求解的边界条件进行计算，得到耦合界面处的热流密度q_c。

（4）对机翼蒙皮内部的多层导热部分进行计算，并更新耦合界面边界的温度分布作为下一循环计算中水膜模型的边界条件。

（5）根据蒙皮表面温度判断当前控制体的结冰类型，计算结冰增长率，并更新冰型。

（6）将冰增长率作为源项加入水膜流动模型并重复上述步骤（3~6）直到迭代至设定的计算时间，计算结束。

2 电加热防除冰计算

2.1　几何模型与计算条件

为验证本研究的飞机电加热防除冰计算方法，采用Al‑Khalil等在NASA Lewis冰风洞所进行的非稳态电加热防除冰实验^［

20］进行对比验证。该实验采用NACA 0012翼型，弦长为0.914 4 m，此翼型采用了多层复合材料的蒙皮结构，多层结构与加热片A~G的分布如图 3所示。

图3 多层复合蒙皮结构与电加热分区示意图

Fig.3 Schematic diagram of multi‑layer composite skin structure and electric heating zoning

蒙皮结构中的多层材料物性参数与实验工况条件分别如表 1、2所示。

表1 复合蒙皮各层材料物性参数

Table 1 Material property parameters of each layer of composite skin

材料	厚度/ mm	导热率/ （W∙（m•K）^-1）	密度/ （kg•m^-3）	定压比热容/ （J• （kg•K）^-1）
耐磨表层	0.203	16.3	80 225.25	502.32
加热片	0.013	41.018	8 906.26	385.112
弹性体绝缘层	0.279	0.256	1 383.99	1 255.88
玻璃纤维/ 环氧树脂复材	0.889	0.294	1 794.07	1 569.75
硅胶泡棉绝热层	3.429	0.121	648.75	1 130.22

表2 实验条件

Table 2 Experimental condition

参数	数值
环境温度/K	266.5
来流速度/（m•s^-1）	44.7
迎角/（°）	0
液态水含量/（g•m^-3）	0.78
过冷水滴等效直径/μｍ	20
机翼表面初始温度/K	264.2

试验采用了A~G七个不同的加热区域，在实验过程中控制不同加热片的热流密度以及加热时间，加热区域的位置以及各加热片的控制率如表 3、4所示。试验中电加热的控制率单循环时间为120 s，其中加热片A在整个实验过程中保持开启状态。

表3 加热片位置

Table 3 Position of heating elements

加热片	s/c
加热片	左侧	右侧
G	-0.102	-0.061
E	-0.061	-0.033
C	-0.033	-0.005
A	-0.005	0.016
B	0.016	0.043
D	0.043	0.071
F	0.071	0.113

表4 电加热系统控制率

Table 4 Control rate of electric heating system

加热片	热流密度/（W•m^-2）	开启时间/s
A	7 750	0~120
B、C	15 500	100~110
D、E、F、G	12 400	110~120

2.2　电加热防除冰计算验证

根据实验条件以及电加热控制率，使用NUAAICE进行计算，将计算结果与试验结果进行对比验证如图 4所示。通过对比加热片A、B、D的平均温度随时间变化曲线，可以看出计算结果中加热区域D的升温速率较快，在模拟过程的20 s左右便基本达到了稳定温度，而实验结果则在40 s左右达到稳定温度，此外在温度峰值方面，加热片B以及加热片D的计算结果与实验结果均相近，加热片A的计算结果相对偏小。可以看出，本研究的电加热防除冰计算模型所得到的加热区域温度响应曲线与实验结果间的偏差在合理范围之内。

图4 加热片A、B、D的温度响应曲线计算结果对比

Fig.4 Comparison of calculated temperature response curves of heating plates A, B, and D

在与FENSAP计算结果进行对比时可以发现，本研究提出基于Myers模型的计算方法由于引入了非定常水膜流动方程，对水膜内部温度分布进行计算并与加热组件导热方程进行耦合，对电加热系统温度响应过程的模拟更加贴近实验结果。

图 5为不同时刻机翼表面的温度分布的计算结果与实验结果对比图。为研究飞机电加热防除冰过程中的耦合传热特性，选取加热时间分别为100，105和110 s时的温度曲线进行对比分析，其中100 s时加热片B、C开始工作，110 s时停止工作。从图中可以看到NUAAICE的计算结果与实验结果温度分布曲线基本吻合，在温度分布范围方面，计算结果略小于实验结果，差异在允许的误差范围内。在100 s时刻，由于加热片B、C的开启，机翼表面存在积冰发生相变融化为水膜，因此机翼表面的温度停留在冰点，当积冰完全融化为水膜后，机翼表面温度则继续升高。

图5 不同时刻机翼表面温度分布曲线与实验结果对比

Fig.5 Comparison of wing surface temperature distribution curves and experimental results at different time points

不同时间机翼内部温度分布云图如图 6所示，可以明显看出加热片A的持续开启在电加热防除冰系统中起到热刀的结冰防护作用，尽管其加热功率较低，但随着时间增加，加热片温度可以稳定在285 K，保证机翼蒙皮温度高于冰点，实现防冰目的。其余的加热片以较高的功率按照控制率设定的顺序依次开启，能在短时间内有效地除去其防护范围表面的积冰。

图6 不同时刻蒙皮内部电加热组件温度云图

Fig.6 Temperature cloud map of internal electric heating components of the skin at different time

本研究使用共轭耦合法对机翼内部电加热组件导热模型与水膜流动模型进行耦合计算，此种耦合方法通过交换界面处的边界值实现，图 7为电加热系统开启前50 s的边界热流密度。从图中可以看出在由于加热区域A的表面温度达到冰点以上，因此入射到此区域的水滴不会冻结成冰，而是形成流动的水膜，此区域的边界热流密度较小，而由于A防护区两侧表面的产生了积冰，阻碍水膜流动，导致此处的热流密度较大。

图7 不同时刻水膜与蒙皮交换的热流密度计算结果

Fig.7 Calculation results of heat flux density between water film and skin at different time

图8为电加热系统开启后的机翼表面积冰示意图，从图中可以看出由于加热片A的持续作用，机翼前缘区域在整个模拟过程中均未发生结冰现象，随着加热片B、C在100 s时刻的开启，此区域的积冰在短时间内发生融化，而加热片D、E、F、G于120 s开启后也使其防护区域的积冰融化。由此可以看出此电加热防除冰系统在防护范围内可以保证前缘防冰，并快速地除去表面积冰，有效地避免了大量积冰对飞行安全造成的影响。然而由于防冰过程的水膜流动以及除冰过程冰融化所产生的大量水将发生溢流水再冻结情况，于机翼上下表面的溢流水产生的积冰同样将影响翼型的气动性能。

图8 不同时刻机翼表面结冰情况

Fig.8 Ice formation on wing surface at different time

3 溢流水再冻结及气动分析

3.1　溢流水结冰计算结果

在飞机电加热防/除冰过程中，随着机翼表面温度不断升高，存在积冰的部分开始融化成水向后流动，即溢流水，但由于加热片铺设范围有限，导致热防护区域外的蒙皮温度仍处于冰点以下，此时向后流动的溢流水则会在此处发生结冰。由于溢流水产生的积冰，导致防除冰系统无法达到理想的防护效果，与此同时积冰会破坏机翼的气动外形，将在一定程度上影响飞行安全。

图 9为使用FENSAP计算得出的机翼表面冰增长速率与本研究的计算结果对比图。本文分别选取100~110 s以及110~120 s的冰增长速率曲线对溢流水再冻结问题进行分析，FENSAP软件的水膜与结冰模块基于Messinger模型进行计算，而Myers模型由于引入了非定常水膜流动方程，可以更加准确地对水膜在流动过程中的冻结进行模拟。从图 9（a）中可以看出，在100 s时本研究计算结果的结冰范围相对较大，而图 9（b）中120 s时的结冰速率显示本计算模拟出溢流水流动至防护范围外并发生冻结的情况。从计算结果对比中可以看出两模型的计算结果总体趋势较为相似，Myers模型可以有效地增加非定常结冰模拟的精度，并且对溢流水冻结现象进行预测。

图9 FENSAP与本研究冰增长速率对比

Fig.9 Comparison of ice growth rates between FENSAP and this paper

其次在结冰情况方面，在电加热循环中由于加热片A持续开启，因此机翼前缘并未产生结冰现象，循环前期积冰的产生范围集中在A加热区域两侧，随着加热时间增加，加热片B、C于100 s开启，表面温度逐渐升高，加热片上方的积冰融化并产生水膜向后流动。图 10则详细展示了电加热系统循环后期的溢流水再冻结过程，由图中可以看出在110 s时前缘的积冰已经由于壁面导热融化，其融化后产生的水膜则部分向后流动，并在现有的积冰后侧发生冻结，此处的结冰位置已经处于热防护区域以外，而在后侧加热片开启后，成功地将其防护范围的积冰除去，却无法融化防护范围外的溢流水再冻结。此时的积冰高度较低，但随着结冰过程的继续，流过此处的水膜继续发生冻结，最终导致此处的积冰高度增加，并破坏原有的气动外形而影响飞行安全。

图10 结冰后期机翼表面积冰变化

Fig.10 Ice variation of wing surface area in the later stage of icing

3.2　溢流水结冰气动分析

由3.1节中的计算结果可以看出当电加热片开启时间在110 s后，热防护区域的积冰已被基本去除，此时产生的溢流水发生向后流动并结冰的情况，改变了机翼原有的气动外形，下文将从飞机的升阻力系数等方面对产生溢流水再冻结现象的机翼进行气动分析。

当飞机电加热防除冰系统开启160 s时机翼表面的积冰情况如图 11所示，已经产生了溢流水再冻结的情况，并且结冰位置的冰层高度在随着时间变化不断增加。

图11 溢流水再冻结过程冰型示意图

Fig.11 Schematic diagram of ice shape during the re‑freezing process of overflow water

从结冰前后的马赫数云图（图12）中可以看出在溢流水再冻结位置，由于气动外形的改变，导致此处流场中的气流发生扰动，进而使马赫数分布产生变化，从图中可以发现在结冰前后，马赫数总体分布较为相似，而图11（b）中发生溢流水再冻结后，马赫数的极值出现在积冰的上方。

图12 结冰前后流场马赫数云图

Fig.12 Mach number cloud map of flow field before and after icing

在进行气动分析中，常使用涡识别方法中的Q准则研究流场中湍流对飞机气动性能的影响。Q准则根据速度梯度矩阵的特征值判断流场中是否存在涡，其中Q的负值表示应变速率或黏性应力占主导的区域，而正值表示流场中涡占主导的区域。

因此由Q准则云图（图13）中可以看出在溢流水再冻结导致流场中出现了明显的涡结构，随后气流在向下流运动的过程中，涡结构逐渐变大并趋向于充分发展。

图13 结冰前后Q准则云图及涡结构局部放大图

Fig.13 Q‑criterion cloud map and local enlarged view of vortex structure before and after icing

为明确溢流水再冻结带来的影响，分别对不同迎角下的干净机翼与发生溢流水再冻结后的机翼进行流场计算，对计算获得的阻力系数与升阻比进行分析研究。图 14为发生溢流水再冻结现象前后的机翼气动特性对比，从图（a）阻力系数曲线对比中可以看出，溢流水结冰后机翼的阻力系数整体都略高于干净翼型，机翼迎角大于0°时的阻力系数升高了15%~23%，并且在0°迎角时，由于溢流水结冰破坏了原有的气动外形，导致阻力系数发生了较大的偏差。与此同时，图14（b）为不同迎角的升力系数与阻力系数之比得出的升阻比曲线，从计算结果中可以看出在负迎角状态时，溢流水再冻结对机翼气动特性的影响较小，而当迎角大于0°时，溢流水结冰对周围流场会产生较大的影响，使气流发生不利分离，导致升阻比的下降。当迎角为6°时，升阻比最下降了约21%，对飞机的气动特性造成较大影响。

图14 溢流水再冻结前后气动特性对比图

Fig.14 Comparison of aerodynamic characteristics before and after re‑freezing of overflow water

经过对电加热防除冰过程的数值模拟，可以发现在开启电加热系统时，由于湿防冰的水膜流动或除冰过程中冰融化产生的溢流水都会在机翼后侧产生溢流水再冻结的现象，而溢流水再冻结将在一定程度上对飞机的气动性能产生影响。从计算结果可以看出，随着结冰时间的增加，溢流水再冻结的冰层厚度也随之增加，将更大程度地破坏飞机的气动外形。由于电加热系统设计的防护范围往往主要参考水滴的入射范围，而较少的考虑到溢流水再冻结的问题，因此在电加热防除冰系统的设计过程中应充分考虑水膜流动，改变控制率或热防护范围以解决溢流水再冻结对飞机气动性能的影响。

4 结论

本文提出了基于共轭传热法耦合传热与水膜的计算模型以模拟飞机电加热防除冰过程，对机翼内部加热元件导热过程与机翼外部水膜流动模型进行分别计算，通过交换界面处的热流密度和温度边界条件以实现两模型的耦合，所得结论如下：

（1）本文中所提出的耦合计算模型计算得出的电加热防除冰过程中表面温度分布与实验结果误差较小，因此使用Myers水膜模型与多层非稳态导热模型耦合计算可对电加热防除冰过程进行较为准测的模拟。

（2）湿防冰过程中水膜流动超出热防护区域，除冰过程中积冰融化产生的溢流水在防护范围外的机翼表面将发生溢流水再冻结现象。

（3）通过涡识别以及流场分析可知溢流水再冻结后对流场产生较大扰动，并使机翼前缘处压力系数曲线发生震荡。溢流水再冻结现象不仅会导致电加热系统的工作效果下降，同时还会对飞行安全造成危害。

参考文献

WRIGHT W B. A comparison of numerical methods for the prediction of two-dimensional heat transfer in an electrothermal deicer pad［M］. Washington， USA：National Aeronautics and Space Administration， 1988. [百度学术]

WRIGHT W， DEWITT K， KEITH JR T. Numerical simulation of icing， deicing， and shedding［C］// Proceedings of the 29th Aerospace Sciences Meeting. Reno， Nevada， USA： American Institute of Aeronautics and Astronautics， 1991. [百度学术]

SILVA G， SILVARES O， ZERBINI E. Airfoil anti-ice system modeling and simulation［C］//Proceedings of the 41st Aerospace Sciences Meeting and Exhibit. Reno， Nevada， USA： American Institute of Aeronautics and Astronautics， 2003. [百度学术]

DA SILVA G A L， DE MATTOS SILVARES O， DE JESUS ZERBINI E J G. Numerical simulation of airfoil thermal anti-ice operation， Part 1： Mathematical modelling［J］. Journal of Aircraft， 2007， 44（2）： 627-633. [百度学术]

REID T， BARUZZI G， ALIAGA C. FENSAP-ICE： Application of unsteady CHT to de-icing simulations on a wing with inter-cycle ice formation［C］//Proceedings of the AIAA Atmospheric and Space Environments Conference. Toronto， Ontario， Canada： American Institute of Aeronautics and Astronautics， 2010. [百度学术]

WRIGHT W B， KEITH T G， DE WITT K J. Two-dimensional simulation of electrothermal deicing of aircraft components［J］. Journal of Aircraft， 1989，26（6）： 554-562. [百度学术]

MESSINGER B L. Equilibrium temperature of an unheated icing surface as a function of air speed［J］. Journal of the Aeronautical Sciences， 1953， 20（1）： 29-42. [百度学术]

SONG Y， LI S， ZHANG S. Peridynamic modeling and simulation of thermo-mechanical de-icing process with modified ice failure criterion［J］. Defence Technology， 2021， 17（1）： 15-35. [百度学术]

ESMAEILIFAR E， RAJ L P， MYONG R S. Computational simulation of aircraft electrothermal de-icing using an unsteady formulation of phase change and runback water in a uniﬁed framework［J］. Aerospace Science and Technology， 2022（130）： 1-19. [百度学术]

卜雪琴，林贵平，郁嘉. 机翼电加热防冰表面内外传热的耦合计算［J］. 航空动力学报， 2010（7）： 1491-1495. [百度学术]

BU Xueqin， LIN Guiping， YU Jia. Coupled heat transfer calculation on an airfoil electrothermal anti-icing surface［J］. Journal of Aerospace Power， 2010（7）： 1491-1495. [百度学术]

钟国. 翼型电热防/除冰系统的数值模拟［J］. 航空制造技术， 2011（4）： 75-79. [百度学术]

ZHONG Guo. Simulation of airfoil electro-thermal anti-ice/de-ice system［J］. Aeronautical Manufacturing Technology， 2011（4）： 75-79. [百度学术]

郑玉巧，潘永祥，魏剑峰. 叶片翼型结冰形态及其气动特性［J］. 南京航空航天大学学报，2020， 52（4）： 632-638. [百度学术]

ZHENG Yuqiao， PAN Yongxiang， WEI Jianfeng. Icing morphology and aerodynamic characteristics of blade airfoil［J］. Journal of Nanjing University of Aeronautics & Astronautics， 2020， 52（4）： 632⁃638. [百度学术]

ZHU Chengxiang， ZHU Chunling， FU Bin. Simulation of SLD impingement on wind turbine blade airfoil［J］. Journal of Nanjing University of Aeronautics & Astronautics， 2016，33（1）： 112-120. [百度学术]

辛苗，钟国，曹义华. 基于水膜流动与耦合传热的翼型防/除冰计算［J］. 航空动力学报， 2021， 36（4）：783-794. [百度学术]

XIN Miao， ZHONG Guo， CAO Yihua. Calculation of airfoil ant-icing/deicing characteristics based on water film flow and coupleoheat transfer［J］.Journal of Aerospace Power，2021，36（4）： 783-794. [百度学术]

WANG Jingxin， YU Dachuan， YANG Zaili. Experimental investigation of super‑hydrophobic/electro‑thermal synergistically anti‑icing/de‑icing strategy in ice wind tunnel［J］. Transactions of Nanjing University of Aeronautics and Astronautics， 2023，40（2）：193‑204. [百度学术]

刘宗辉，卜雪琴，林贵平，等. 基于PHengLEI的非稳态电热除冰过程仿真［J］. 空气动力学学报， 2022， 41（2）： 53-64. [百度学术]

LIU Zonghui， BU Xueqin， LIN Guiping， et al. Simulation of unsteady electrothermal deicing process based on PHengLEI［J］. Acta Aerodynamica Sinica， 2022， 41（2）： 53-64. [百度学术]

SHEN X B， QI Z C， ZHAO W Z， et al. Numerical simulation of aircraft icing with an unsteady thermodynamic model considering the development of water film and ice layer［J］. International Journal of Aerospace Engineering， 2022（10）： 1-18. [百度学术]

MYERS T G， CHARPIN J P F. A mathematical model for atmospheric ice accretion and water flow on a cold surface［J］. International Journal of Heat and Mass Transfer， 2004， 47（25）： 5483-5500. [百度学术]

YAU L C. Conjugate heat transfer with the multiphysics coupling library preCICE［D］. Munich， Germany：Technische Universität München， 2016. [百度学术]

WRIGHT W， AL-KHALIL K， MILLER D. Validation of NASA thermal ice protection computer codes II—LEWICE/Thermal［C］//Proceedings of the 35th Aerospace Sciences Meeting and Exhibit. Reno， NV， USA： American Institute of Aeronautics and Astronautics， 1997. [百度学术]

请使用 Firefox、Chrome、IE10、IE11、360极速模式、搜狗极速模式、QQ极速模式等浏览器，其他浏览器不建议使用!

电加热防/除冰过程中溢流水再冻结数值模拟 PDF

摘要

关键词

1 物理模型与方法

1.1 水膜流动模型

1.2 多层导热模型

1.3 水膜流动与传热耦合方法

2 电加热防除冰计算

2.1 几何模型与计算条件

2.2 电加热防除冰计算验证