科技和产业

层数	材料	导热系数/ (W·m^-1·K^-1)	厚度/mm
1	防腐蚀罩	16.270	0.200 0
2	弹性体	0.256	0.280 0
3	加热元件	41.020	0.012 7
4	弹性体	0.256	0.280 0
5	玻璃纤维/ 环氧复合材料	0.294	0.890 0
6	绝缘泡沫层	0.121	3.430 0

层数

材料

导热系数/
(W·m^-1·K^-1)

厚度/mm

防腐蚀罩

16.270

0.200 0

弹性体

0.256

0.280 0

加热元件

41.020

0.012 7

弹性体

0.256

0.280 0

玻璃纤维/
环氧复合材料

0.294

0.890 0

绝缘泡沫层

0.121

3.430 0

层数	材料	导热系数/ (W·m^-1·K^-1)	厚度/mm
1	防腐蚀罩	16.270	0.200 0
2	弹性体	0.256	0.280 0
3	加热元件	41.020	0.012 7
4	弹性体	0.256	0.280 0
5	玻璃纤维/ 环氧复合材料	0.294	0.890 0
6	绝缘泡沫层	0.121	3.430 0

层数

材料

导热系数/
(W·m^-1·K^-1)

厚度/mm

防腐蚀罩

16.270

0.200 0

弹性体

0.256

0.280 0

加热元件

41.020

0.012 7

弹性体

0.256

0.280 0

玻璃纤维/
环氧复合材料

0.294

0.890 0

绝缘泡沫层

0.121

3.430 0

参数	工况1	工况2	工况3
来流温度T_∞/℃	-6.67	-6.67	-17.78
来流速度u_∞/(m·s^-1)	44.7	89.4	44.7
液态水含量LWC/(g·m^-3)	1.10	0.39	0.78
水滴直径MVD/μm	20	20	20

参数

工况1

工况2

工况3

来流温度T_∞/℃

-6.67

-17.78

来流速度u_∞/(m·s^-1)

44.7

89.4

44.7

液态水含量LWC/(g·m^-3)

1.10

0.39

0.78

水滴直径MVD/μm

参数	工况1	工况2	工况3
来流温度T_∞/℃	-6.67	-6.67	-17.78
来流速度u_∞/(m·s^-1)	44.7	89.4	44.7
液态水含量LWC/(g·m^-3)	1.10	0.39	0.78
水滴直径MVD/μm	20	20	20

参数

工况1

工况2

工况3

来流温度T_∞/℃

-6.67

-17.78

来流速度u_∞/(m·s^-1)

44.7

89.4

44.7

液态水含量LWC/(g·m^-3)

1.10

0.39

0.78

水滴直径MVD/μm

加热器	功率密度/(kW·m^-2)
F	2.635	4.030	4.805
D	3.100	5.425	5.890
B	4.030	5.270	6.975
A	4.960	6.510	8.990
C	3.565	5.270	5.890
E	3.565	5.890	6.200
G	2.635	3.720	4.495

加热器

功率密度/(kW·m^-2)

工况1

工况2

工况3

2.635

4.030

4.805

3.100

5.425

5.890

4.030

5.270

6.975

4.960

6.510

8.990

3.565

5.270

5.890

3.565

5.890

6.200

2.635

3.720

4.495

加热器	功率密度/(kW·m^-2)
F	2.635	4.030	4.805
D	3.100	5.425	5.890
B	4.030	5.270	6.975
A	4.960	6.510	8.990
C	3.565	5.270	5.890
E	3.565	5.890	6.200
G	2.635	3.720	4.495

加热器

功率密度/(kW·m^-2)

工况1

工况2

工况3

2.635

4.030

4.805

3.100

5.425

5.890

4.030

5.270

6.975

4.960

6.510

8.990

3.565

5.270

5.890

3.565

5.890

6.200

2.635

3.720

4.495

考虑水膜流动传热的电加热三维防冰数值研究

PDF下载

龚立锋 ¹ , 李世明 ² , 胡娅萍 ² , 刘蕾 ²

科技和产业 | 科技创新 2025,25(13): 79-88

收起

科技和产业 | 科技创新 2025, 25(13): 79-88

考虑水膜流动传热的电加热三维防冰数值研究

全屏

龚立锋¹, 李世明², 胡娅萍², 刘蕾²

作者信息

¹ 中国航发湖南动力机械研究所, 湖南株洲 412002

² 南京航空航天大学能源与动力学院, 南京 210016

龚立锋(1981—),男,河南光山人,硕士,副研究员,研究方向为适航;

李世明(1999—),男,河南郑州人,硕士研究生,研究方向为航空发动机结冰防冰模拟;

胡娅萍(1978—),女,湖北鄂州人,博士,副教授,研究方向为航空发动机结冰防冰模拟;

刘蕾(1997—),男,江西萍乡人,硕士研究生,研究方向为航空发动机结冰防冰模拟。

Numerical Study on Electrically Heated 3D Anti-icing Considering Water Film Flow and Heat Transfer

Lifeng GONG¹, Shiming LI², Yaping HU², Lei LIU²

Affiliations

¹ AECC Hunan Aviation Powerplant Research Institute, Zhuzhou 412002, Hunan, China

² College of Energy and Power Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China

出版时间: 2025-07-10

文章导航

摘要

收起

飞机结冰对飞行安全造成严重威胁,电加热是一种高效的防除冰方式。建立考虑水膜流动、传热的电加热防冰三维数学模型,并提出电加热防冰数值计算方法。以NACA0012翼型为研究对象,对持续性电加热条件下的稳态防冰过程进行数值模拟,通过将计算结果与现有实验数据进行对比,验证该计算方法的准确性。计算结果表明:加热功率较小时,水膜流出加热区域并在其下游生成溢流冰;同一来流条件下,加热功率越大,防冰表面温度越高,水膜覆盖的范围越小。此外,防冰表面温度计算结果与实验数据的误差在5 ℃以内,对于来流温度不低于6.67 ℃的工况,其误差小于3 ℃。

关键词

电加热 / 防冰 / 水膜流动 / 三维数值模拟

Abstract

收起

Icing on aircraft poses a serious threat to flight safety, and electric heating is an efficient method for anti-icing and de-icing. A three-dimensional mathematical model that considers water film flow and heat transfer for electric heating anti-icing was constructed. Additionally, a numerical calculation method for electric heating anti-icing was proposed. This method was applied to simulate the steady-state anti-icing process of the NACA0012 airfoil under continuous electric heating conditions. The accuracy of this calculation method was validated by comparison with existing experimental data and computational results. The results indicate that when the heating power is low, the water film flows out of the heated area and overflow ice forms downstream. Under the same inflow conditions, a higher heating power results in a higher anti-icing surface temperature and a smaller water film coverage area. Furthermore, the calculated anti-icing surface temperatures are within 5 ℃ error compared to experimental data, and for inflow temperatures no lower than 6.67 ℃, the error is less than 3 ℃.

Key words

electric heating / anti-icing / water film flow / 3D numerical simulation

引用本文

龚立锋, 李世明, 胡娅萍, 刘蕾. 考虑水膜流动传热的电加热三维防冰数值研究. 科技和产业, 2025 , 25 (13) : 79 -88 .

Lifeng GONG, Shiming LI, Yaping HU, Lei LIU. Numerical Study on Electrically Heated 3D Anti-icing Considering Water Film Flow and Heat Transfer[J]. Science Technology and Industry, 2025 , 25 (13) : 79 -88 .

正文

收起

云层中的过冷水滴撞击到飞机迎风表面时,水滴会发生破碎,覆盖在飞机表面铺展开并凝固成冰。结冰会使飞机升力降低、飞行阻力升高,严重影响飞行安全。因此,需要对飞机结冰部件采取防除冰措施。

电加热防除冰是现代飞机最广泛使用的防除冰方法。该方法将电加热片或其他加热元件通电产生的热量通过导热传递到防冰部件表面,使部件表面温度处于冰点以上,从而防止水凝固成冰,或将部件表面已凝结的冰除去,以达到防冰或除冰的目的。采用电加热进行防除冰也是飞机和风力涡轮机应用最有效的方法,因为该方法不需要从发动机引气,每个加热区域可以单独控制,避免过热或不足,提高了防除冰的效率。

在数值模拟方面,国外早期研究为针对多层复合材料组成的加热器固体进行的一维或二维导热建模与计算^[1-3]。之后,Wright^[4]开展了二维结冰、防除冰及冰脱落过程的数值研究。由此发展了非稳态电加热防冰/除冰的计算软件LEWICE/Thermal^[5]。Al-Khalil和Potapczuk^[6]、Dean等^[7]建立了一种沿流线的二维稳态防冰模型,其计算结果与冰风洞实验结果比较接近,并发展了稳态的气热和电加热防冰计算软件ANTICE。Silva等^[8-11]建立了一系列计算精度越来越高的二维防冰模型,包括溪流模型。Bennani等^[12]发展了考虑冰层融化和脱落的二维电加热防冰模型,基于焓法模型处理冰的融化。Remi等^[13]采用改进的Schwarz方法求解水膜质量、能量守恒方程和固体导热方程,该耦合算法可使求解过程快速收敛。Enache等^[14]发展了一维非等温相变模型,用于计算冰脱落时冰层底部融化形成的水膜厚度的临界值。在国内,常士楠等^[15-17]开发了模拟二维电热除冰过程的程序,采用焓法处理冰的相变,并与实验结果对比验证了算法的准确性。肖春华等^[18-19]研究了二维电加热除冰过程,计算了考虑冰脱落的除冰表面的温度变化。熊建军等^[20]研制了基于状态机的控制律电加热防除冰控制系统,并利用结冰风洞试验对不同控制律电加热防除冰的有效性进行了验证。

在三维数值研究方面,Reid等^[21]研究了FENSAP-ICE软件中电加热除冰的数学模型和算法,并通过与实验数据对比验证了其计算精度。Esmaeilifar等^[22]提出了一种利用非均匀热边界条件对电热防冰进行数值模拟,根据瞬态壁温分布更新对流传热系数提高求解精度。刘小暄等^[23]使用共轭耦合方法建立机翼内部多层导热与水膜流动的耦合计算模型,发现积冰融化和水膜流动会产生溢流水再冻结现象。杨诗雨等^[24]对Fluent软件进行二次开发,计算得到旋转帽罩电加热防冰表面的温度变化。吴佩佩等^[25]使用FENSAP-ICE软件计算了进气道唇口电加热防冰系统在不同飞行高度、速度及攻角下的防冰特性,计算结果与实验数据吻合良好。

综上所述,国内外针对电加热防除冰的数值模拟大多基于二维数学模型,难以直接应用于三维数值研究。由于数学模型的限制,目前对于电加热部件表面水膜三维流动和传热特性的相关研究较少,而这些特性对于防冰效果至关重要。因此,为了更准确地模拟部件表面电加热防冰过程,本文建立考虑水膜流动、传热的电加热防冰三维数学模型,并提出电加热防冰数值计算方法。再以NACA0012翼型为研究对象,对其在持续性电加热工况下的稳态防冰过程进行数值模拟,并将计算结果与文献实验数据进行对比,验证本文计算方法的精度。

1 电加热稳态防冰三维数学模型

收起

1.1 空气-过冷水滴两相流

空气-过冷水滴两相流的控制方程与文献[26]中的一致,即两相流动采用欧拉-欧拉法求解,不考虑过冷水滴对空气流的影响。应用Fluent软件的UDS(User-Defined Scalar)开发功能,计算部件外两相流场。过冷水滴撞击部件表面的局部水收集系数β的计算式为

(1)

β = - α w ρ w u d, n o r m a l L W C u d, ∞

式中:

ρ w

为水滴的密度;

u d, n o r m a l

为水滴撞击到部件表面的法向速度;

α w

为水滴在撞击处的局部体积分数;

u d, ∞

为过冷水滴的来流速度;LWC为来流的液态水含量。

1.2 水膜流动与相变

考虑水膜流动、传热与相变的三维防冰计算的微元控制体如图1所示。当电加热提供的热量不足时,部件表面会有结冰,此时情况最复杂,控制体内自上而下依次为空气-过冷水滴两相流、水膜层和冰层。

根据部件防冰表面可能发生的情况,结合图1所示的控制体,建立描述部件表面水膜流动、传热与相变的连续方程、动量方程和能量方程。

1.2.1 水膜连续方程

考虑部件表面水膜流动过程中过冷水滴的撞击加入、水膜蒸发和结冰相变,建立其连续方程为

(2)

ρ w ∂ H w ∂ t + ∫ ρ w u - w n d A A s u b = m · i m p - m · e v a p - ρ i c e ∂ H i c e ∂ t

式中:t为水膜流动时间;

u - w

为水膜流在控制面的平均速度;

n

为控制面的法向单位矢量;A为控制面面积;

m · i m p

为单位面积上水滴撞击的质量流量;

m · e v a p

为单位面积上水膜蒸发的质量流量;

ρ i c e

为冰的密度。

1.2.2 水膜动量方程

水膜流动可近似为不可压层流流动,由此建立其三维流动的动量方程为

(3)

- ∂ p ∂ x + μ w ∂ 2 u w, x ∂ z 2 + F x = 0 - ∂ p ∂ y + μ w ∂ 2 u w, y ∂ z 2 + F y = 0 - ∂ p ∂ z + μ w ∂ 2 u w, z ∂ z 2 + F z = 0

式中:p为压力;

μ w

为水的动力黏度;

u w, x 、 u w, y 、 u w, z

为水膜速度在3个坐标方向的分量,z坐标轴正方向为贴体坐标系下垂直于部件表面向外的法向;

F x 、 F y 、 F z

为水膜所受体积力源项在3个坐标方向的分量。水膜上表面受空气剪切力作用,水膜下表面与冰面或部件表面接触,为无滑移边界,则动量方程的边界条件为

(4)

μ w ∂ u w ∂ z = τ, z = H w u w = 0, z = 0

式中:

τ

为空气的剪切力矢量。

1.2.3 水膜能量方程

水膜流动的能量守恒方程为

(5)

ρ w ∂ (H w c w T - w) ∂ t A s u b + ρ w ∫ (c w T - w u w n) d A s u b = (q · e v a p + q · c o n d + q · c o n v - q · i m p - q · i c e) A s u b

式中:

c w

为水的比热容;

T - w

为水膜温度沿其厚度方向的平均值;

q · e v a p

为水膜蒸发的热流密度;

q · c o n d

为水膜与冰面或部件表面传递的热流密度;

q · c o n v

为水膜与空气间的对流换热热流密度;

q · i m p

为撞击过冷水滴带入的能量;

q · i c e

为结冰放出的潜热。

1.2.4 水膜厚度方向的温度分布方程

因防冰表面的水膜厚度通常很小,故假设水膜温度沿其厚度方向呈线性变化,则有

(6)

∂ T w ∂ z = c o n s t (常 量)

水膜的上表面与外部空气-过冷水滴两相流接触,采用第2类边界条件:

(7)

- λ w ∂ T w ∂ z = - q · i m p + q · e v a p + q · c o n v, z = H w

式中:

λ w

为水膜的导热系数。水膜下表面若紧邻冰层(当电加热量不足时部件表面仍结冰),则水膜下表面温度与冰点温度T_f相等,其边界条件为

(8)

T w = T f, z = 0

若水膜下表面与被加热部件表面直接接触(电加热量较大时部件表面不结冰),其边界条件为

(9)

- λ w ∂ T w ∂ z = q · c o n d, z = 0

1.3 固体导热

电加热元件产生的热量通过部件固体传导至其外表面。部件固体内的导热方程为

(10)

∂ ∂ t (ρ h) + ∇ · (u i ρ h) = ∇ · (λ s ∇ · T)

式中:h为固体的焓;λ_s为固体的导热系数;T为固体的温度;

∇ · (u i ρ h)

为因固体平移或旋转运动传输的热量。

2 电加热稳态防冰数值计算方法

收起

采用Fluent软件及二次开发的程序求解电加热稳态防冰特性,计算方法及整体计算流程如图2所示,主要包括4个计算模块:电加热条件下部件外空气流动传热与部件固体导热耦合计算、过冷水滴流动与撞击特性计算、水膜流动传热与相变计算、防冰表面温度计算。具体计算流程如下。

(1)电加热条件下外部空气流动传热与部件固体导热耦合计算。在计算外部空气流场时,同时启用电加热器,通过热流边界的方式将加热功率添加到固体内壁面。利用Fluent软件的单相流和固体导热求解器计算得到外部空气的压力、速度、部件外壁面空气剪切力以及加热热流分布等。

(2)过冷水滴流动与撞击特性计算。在完成部件外空气流场计算后,使用Fluent软件编译自主开发的UDS程序,将过冷水滴的体积分数和水滴速度的3个坐标分量作为4个标量添加到软件中,并通过求解这4个标量的输运方程获得过冷水滴的流场分布。然后,利用自主开发的UDF程序中的宏命令计算表面局部水收集系数和过冷水滴的撞击质量流量。

(3)水膜流动、传热与相变计算。在Fluent中执行自主开发的UDF程序,基于(1)和(2)计算得到的两相流结果。首先,求解水膜的温度分布方程和水膜能量方程,得到结冰量和冰层厚度;然后求解动量方程,获得水膜流动速度;最后,迭代求解水膜连续方程,得到水膜厚度。若当前时间步的水膜厚度计算收敛,判断防冰计算时间是否达到;若没有达到,则继续进行下一时间步的计算。本模块计算完成后,更新由于过冷水滴撞击、表面水膜流动传热和相变引起的能量项

q · i m p + q · i c e + q · e v a p

。

(4)防冰表面温度计算。将(3)计算得到的能量项

q · i m p + q · i c e + q · e v a p

作为部件固体的热源条件,以UDF源项的形式添加到部件表面。然后,再次进行外部空气流动传热与固体导热的耦合计算,即可获得防冰表面的温度分布。

3 电加热稳态防冰数值模拟

收起

3.1 计算模型

3.1.1 物理模型

本文数值计算结果将与文献[27]中NASA(美国航空航天局)的实验数据进行对比验证,因此本文按该文献研究的翼型建立计算模型。图3为所研究的电加热NACA0012翼型的物理模型,其展向长度为1.828 8 m、弦长为0.914 4 m。翼型前缘由多层复合材料构成,分布有A、B、C、D、E、F、G 7个独立可控的电加热器,上下表面各布置3个加热器,前缘驻点附近布置一个加热器A,如图3所示。前缘的多层结构从外向内分别由防腐蚀不锈钢罩、薄加热元件和厚绝缘层等组成。表1列出了各层材料的导热系数和厚度。

文献[27]中实验数据涉及翼型防冰表面的温度分布、表面水膜流动和溢流冰的分布。与翼型固体内部的物理变化相比,机翼表面发生的物理现象更值得关注。因此,本文参考文献[28],将多层材料组成的NACA0012翼型固体域等效简化为单一材料的固体域。简化后的固体域厚度为5 mm,等效导热系数为0.918 2 W/(m·K),电加热元件分布在固体内壁面。

3.1.2 计算域及边界条件

图4为NACA0012翼型外部两相流与防冰计算的计算域。翼型上游流场取半径为翼型弦长3倍的半圆柱域,下游流场取长度为4倍弦长、宽度为6倍弦长的长方体域。

将计算域的进口边界设为速度进口,出口边界设为压力出口,压力为1个标准大气压(1标准大气压=101 325 Pa),计算域的上下两个端面设为对称边界,机翼的内表面电加热区域设置为热流边界,热流大小取决于每个工况的电加热功率,翼型外表面设置为无滑移耦合传热壁面。本文开展3个工况下的稳态防冰数值计算,防冰时间均为600 s,且600 s内电加热器一直处于开启态。其中,每个工况的来流参数及其加热功率分别见表2^[27]和表3^[28]。

3.1.3 网格划分及独立性验证

采用ICEM-CFD软件对计算域进行网格划分,如图5所示。整个计算域划分为六面体结构化网格,分别设置流体区域和固体区域,并对翼型附近的网格进行加密,以保证数值计算的精度,如图6所示,其中灰色区域表示流体域,橙色区域表示固体域。外部空气流场求解采用可实现的k-ε湍流模型,结合Scalable型壁面函数;根据湍流计算对翼型近壁面网格的要求,翼型壁面外第1层网格高度设为0.9 mm,对应的反映壁面附近流体黏性影响的无量纲参数y⁺约为30。

为了消除网格疏密度对计算结果的影响,在工况1的来流条件下,分别采用总单元数为34.5万、38.7万、47.1万、53万和59.8万的5套网格进行空气-过冷水滴两相流场的计算,通过监测并对比翼型驻点附近上游某处的水滴速度进行网格独立性的验证。当网格总数大于53万时,监测位置处的水滴速度大小基本保持不变,如图7所示,最终采用单元总数为53万的网格进行计算。

3.2 两相流场计算结果

对表2和表3中的3个防冰工况进行考虑水膜流动传热的电加热三维防冰数值计算。因3个工况的来流条件相近,故对于两相流场,以工况1的计算结果为代表进行分析。为方便分析,沿着垂直于翼型展向在计算域中部截取一个中截面,如图8所示,接下来主要分析该中截面上的计算结果。

工况1中截面上空气的相对压力分布如图9所示,为了显示更多细节,仅展示出翼型周围区域的流场,可以看出,翼型驻点附近的压力最大,相对压力值约为1.33 kPa,沿着流向往下游,翼型的上下表面压力有所降低,小于标准大气压,即翼型上下表面附近出现相对负压。NACA0012翼型为上下对称结构,因此翼型上下表面两侧压力呈对称分布。为了防冰内部加热单元的设计,机翼尾缘处被切成一个台阶状,尾缘附近压力有所升高。

图10为工况1中截面上翼型附近空气速度分布,可以看出,翼型驻点附近形成滞止区域,空气速度很小。由于流道减小,翼型上下表面两侧区域的气流速度增大,最大速度约52.3 m/s,这与图10中两侧的负压区域对应。在机翼尾缘处,被机翼分开的上下两股气体又重新汇合,速度大小发生变化。

图11给出了工况1中截面上过冷水滴在翼型附近的速度分布。可以看出,水滴相和空气相的速度分布基本相同,翼型驻点附近速度比较小,上下表面两侧区域呈现局部加速、尾部有涡的脱落。水滴的惯性大,其最大速度稍大于空气相的最大速度,尾缘处流道发生变化,导致水滴速度有所减小。

3.3 稳态防冰计算结果与分析

对3个持续性电加热工况的翼型防冰表面的局部水收集系数、温度、水膜和结冰厚度等参数的分布进行分析,并将防冰表面温度与实验数据和国外仿真软件的计算结果进行对比,验证自主开发的程序的计算精度。

3.3.1 工况1

采用穿透边界计算表面的局部水收集系数。为了更直观地分析,截取图8所示的中截面与翼型表面相交的型线,定义该型线在驻点处的弧长坐标S为0,翼型上表面该曲线的弧长坐标S>0,向下游方向S逐渐增大;翼型下表面对应S<0,向下游方向S的绝对值逐渐增大。

工况1翼型表面局部水收集系数沿弧长的分布如图12所示,将本文计算结果与ANTICET和ONERA2D软件模拟结果进行对比。可以看出,在驻点处,本文程序计算的最大局部水收集系数与ONERA2D的基本相同,都大于ANTICE的预测结果;局部水收集系数沿弧长的分布范围与ONERA2D软件的基本相同,都小于ANTICE预测的撞击范围。

提取翼型中截面型线上的防冰温度,并与NASA的实验数据及ANTICE软件的计算结果进行对比,如图13所示。可见,本文程序计算的表面温度与实验数据非常接近,最大误差不超过3 ℃,而且与ANTICE预测的温度变化趋势基本相同。在驻点附近水滴撞击区,表面温度变化平缓,在加热区域以外,表面温度陡降。相比于实验温度值,在驻点和上下游的最后一个加热器附近区域,本文程序模拟的结果优于ANTICE软件,证明了本文自主开发的程序具有较好的防冰预测效果。

工况1翼型中截面型线上水膜和溢流冰分布分别如图14、图15所示,其中图右侧的纵坐标Q代表单位面积上的加热功率。可以看出,水膜不但覆盖了整个加热区域,而且溢流到了加热区域以外。因翼型上表面的加热区域稍大于下表面,故而水膜在机翼上表面覆盖范围稍大于下表面。另外,加热区域表面温度较高,而加热区域以外表面温度低于0 ℃,上表面水膜覆盖范围较广,进而导致溢流冰分布较多,最大冰厚约4.3 mm。

定义来流方向为x正方向,垂直于翼型弦长的方向为y方向,展向为z方向。工况1翼型中截面型线上水膜速度沿x、y、z 3个坐标方向上的速度分布如图16所示,可以看出水膜流动的速度非常小,翼型上表面水膜流动的范围略大于下表面,与图14的水膜分布符合。驻点处水膜流动的速度为0,水膜从驻点向翼型表面下游流动的速度先明显增大,然后以相对较大的速度在加热区域流动,流出加热区域后速度迅速减小为0。沿x方向水膜流动的最大速度约为0.041 m/s,y方向最大速度约0.025 m/s。由于攻角为零,水膜在展向(z方向)的运动速度为0。

3.3.2 工况2

工况2的翼型表面局部水收集系数沿弧长的分布如图17所示。由于文献没有记载该工况下CANICE软件计算的局部水收集系数,因此将本文计算结果与ANTICE软件及法国防冰软件ONERA2D的计算结果进行对比。可以看出,在驻点处,本文计算的局部水收集系数最大,计算的撞击极限介于ANTICE软件与ONERA2D软件的计算结果之间,三者计算结果在驻点下游附近的区域基本相同。

另外,由于工况2的来流速度是工况1的两倍,两种工况的水滴平均直径相同,因此对比图17和图12可见,来流速度较大时翼型表面的局部水收集系数也越大。

图18为翼型中截面型线上防冰温度分布与实验数据和ANTICE软件的计算结果对比图。工况2来流速度较大,除了表面部分水膜蒸发带走一部分热量,对流换热也带走一部分热量,比工况1的更显著,故防冰表面水滴撞击区域的最高温度仅2 ℃,整个加热区域最高温度为4.86 ℃,加热区域外,表面温度迅速下降。本文计算的表面温度与实验数据非常接近,最大误差不超过3 ℃,证明了本文自主开发的程序具有较高的计算精度,并且与ANTICE软件预测的温度变化趋势相同。

图19和图20分别显示了工况2翼型中截面型线上的水膜和溢流冰的厚度及其分布。可以看出,工况2的表面水膜较薄,最厚不超过6 μm,加热区域以外的水膜厚度迅速减小为0,加热区域没有结冰,水膜流到加热区域以外冻结形成溢流冰,最大厚度达到了7 mm。

图21为翼型中截面上水膜速度沿x、y、z 3个坐标方向上的速度分布。可以看出,由于工况2来流速度较大,导致水膜流动的速度也比工况1的大,沿x方向水膜流动的最大速度约0.082 m/s,y方向最大速度约0.047 m/s,流出加热区域后速度很快减小为0。由于该工况的来流速度是工况1的两倍,并结合工况1的水膜流动速度(图16)对比可发现,该工况水膜流动速度最大值也大约是工况1的两倍。

3.3.3 工况3

工况3与工况1的来流速度和水滴平均直径均相同,计算得到的防冰表面水滴撞击特性也相同,故不再赘述。

工况3的翼型中截面型线上的温度与实验数据对比如图22所示。由图可见,本文计算得到的表面温度与实验数据的误差不超过5 ℃,该工况的来流温度较低,加热区域以外,表面温度很快降低到0 ℃以下。

工况3翼型中截面型线上的水膜和溢流冰的厚度及其分布如图23和图24所示。可以看出,沿机翼的上下表面布置的电加热器的加热功率几乎呈对称分布,水膜沿机翼的上下表面也基本上呈对称分布,加热区域的水膜厚度维持在10 μm左右,加热区域以外的水膜厚度迅速减小为0。整个加热区域没有结冰,由于翼型外部环境温度很低,一部分水膜流到加热区域以外冻结形成溢流冰,溢流冰沿机翼表面的覆盖范围不大,最大冰厚约为9 mm。

4 结论

收起

针对NACA0012翼型建立三维计算模型,采用商业软件Fluent及自主开发的程序,对3个持续性电加热工况的防冰过程开展表面水滴撞击特性、防冰表面温度分布验证、表面水膜和溢流冰厚度、表面水膜速度分布等研究。主要研究结论如下。

(1)建立了考虑水膜流动和传热的电加热防冰三维数学模型和计算方法,突破了以往主要基于二维模型的研究局限。

(2)在研究的参数范围内,来流水滴直径均相等,当来流速度增大时,翼型表面局部水收集系数随之增大,来流速度相同时,翼型表面局部水收集系数基本不变。

(3)翼型防冰表面温度的计算结果与实验数据的误差基本都在5 ℃以内,对于来流温度不低于-6.67 ℃的工况,计算误差不超过3 ℃,验证了自主开发的程序具有较高的计算精度。

(4)当加热功率密度小于10 kW/m²时,防冰表面加热区域的温度较低,但都高于冰点,表面加热区域不结冰,水膜向下游流动,在表面加热区域的下游生成溢流冰。

本文的研究成果可应用于飞机部件、航空发动机进气道与进气部件、无人机等的热防冰模拟、防冰设计及优化,有望提高飞机在结冰气象条件下的飞行安全性,降低结冰风险,同时优化加热功率分配,提高能源利用效率。

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

STALLABRASS

J R

. Thermal aspects of deicer design[C]// International Helicopter Icing Conference. Ottawa: SAE: 1972: 410057.

[2]

MARANO

J J

. Numerical simulation of an elect-rothermal deicer pad[D]. Toledo: University of Toledo, 1983.

[3]

MASIULANIEC

K C

. A numerical simulation of the full two-dimensionalelectrothermal deicer pad (Heat transfer)[D]. Toledo: University of Toledo, 1987.

[4]

WRIGHT

W B

. A comparison of numerical methods for the prediction of two-dimensional heat transfer in anelectrothermal deicer pad[M]. Washington: NASA, 1988.

[5]

AL-KHALIL

K M

. Validation of NASA thermal ice protection computer codes. part 2-LEWICE/ THERMAL[C]// 35th Aerospace Sciences Meeting and Exhibit. Nevada: AIAA: 1997: 970050.

[6]

AL-KHALIL

, POTAPCZUK

. Numrical modeling of anti-icing systems and comparison to test results on a NACA 0012 airfoil[C]// 31st Aerospace Sciences Meeting. Nevada: AIAA, 1993: 930170.

[7]

DEAN

, THOMAS

, DAVID

, et al. Validation of NASA thermal ice protection computer codes.Part 3-Validation of ANTICE[C]// 35th Aerospace Sciences Meeting and Exhibit. Nevada: AIAA, 1997: 970051.

[8]

SILVA

, SILVARES

, ZERBINI

. Airfoil anti-ice system modeling and simulation[C]// 41st Aerospace Sciences Meeting and Exhibit. Nevada: AIAA, 2003: 030734.

[9]

SILVA

, SILVARES

, ZERBINI

. Simulation of an airfoil electro-thermal anti-ice system operating in running wet regime[C]// 43rd Aerospace Sciences Meeting and Exhibit. Nevada: AIAA, 2005: 051374.

[10]

SILVA

, SILVARES

, ZERBINI

. Water film breakdown and rivulets formation effects on thermal anti-ice operation simulation[C]// 9th AIAA/ASME Joint Thermophysics and Heat Transfer Conference. California: AIAA/ASME, 2006: 3785.

[11]

SILVA

G A

, SILVARES

, ZERBINI

. Numerical simulation of airfoil thermal anti-ice operation, part 1: mathematical modelling[J]. Journal of Aircraft, 2007, 44(2): 627-633.

[12]

BENNANI

, VILLEDIEU

, SALAUN

. Two dimensional model of an electro-thermal ice protection system[C]// 5th AIAA Atmospheric and Space Environments Conference. San Diego: AIAA, 2013: 2936.

[13]

REMI

, VILLEDIEU

, TRONTIN

, et al. A robust coupling algorithm applied to thermal ice protection system unsteady modeling[C]// 6th AIAA Atmospheric and Space Environments Conference. Atlanta: AIAA, 2014: 2061.

[14]

ENACHE

, BERNAY

, GLABEKE

, et al. Ice shedding phenomenon: an experimental and numerical investigation[C]// AIAA Aviation 2020 Forum. San Diego: AIAA, 2020: 2817.

[15]

常士楠, 候雅琴, 袁修干. 周期电加热控制律对除冰表面温度的影响[J]. 航空动力学报, 2007, 22(8): 1247-1251.

[16]

常士楠, 艾素霄, 霍西恒, 等. 改进的电热除冰系统仿真[J]. 航空动力学报, 2008, 23(10): 1753-1758.

[17]

CHANG

S N

. Numerical simulation of electro-thermal deicing process and the effects of icing conditions[C]// International Conference on Aircraft and Engine Icing and Ground Deicing. Chicago: SAE, 2011: 380096.

[18]

肖春华. 飞机电热除冰过程的传热特性及其影响研究[D]. 四川: 中国空气动力研究与发展中心, 2010.

[19]

肖春华, 林贵平, 桂业伟, 等. 冰脱落对电热除冰传热特性的影响研究[J]. 空气动力学学报, 2012, 30(4): 551-556.

[20]

熊建军, 刘锡, 冉林, 等. 基于控制律的电加热防除冰系统设计与验证[J]. 测控技术, 2021, 40(2): 130-134, 139.

[21]

REID

, BARUZZI

G S

, HABASHI

W G

. FENSAP-ICE: unsteady conjugate heat transfer simulation ofelectrothermal de-icing[J]. Journal of Aircraft, 2012, 49(4): 1101-1109.

[22]

ESMAEILIFAR

, SENGUPTA

, RAJ

L P

, et al. In-flight anti-icing simulation of electrothermal ice protection systems with inhomogeneous thermal boundary condition[J]. Aerospace Science and Technology, 2024, 150: 109210.

[23]

刘小暄, 朱春玲, 朱程香, 等. 电加热防/除冰过程中溢流水再冻结数值模拟[J]. 南京航空航天大学学报, 2024, 56(2): 281-290.

[24]

杨诗雨, 常士楠, 高艳欣, 等. 旋转帽罩电加热防冰瞬态过程研究[J]. 空气动力学学报, 2016, 34(3): 289-294, 307.

[25]

吴佩佩, 晏涛, 马赛强, 等. 发动机唇口电热防冰系统性能仿真[J]. 航空动力学报, 2020, 35(10): 2056-2063.

[26]

FLUENT

. Ansys fluent theory guide[R]. Pittsburgh, Pennsylvania: ANSYS Inc., USA, 2011.

[27]

MILLER

, BOND

, SHELDON

, et al. Validation of NASA thermal ice protection computer codes. Part 1-Program overview[C]// 35th Aerospace Sciences Meeting and Exhibit. Nevada: AIAA, 1997: 970049.

[28]

MORENCY

, TEZOK

, PARASCHIVOIU

. Anti-icing system simulation using CANICE[J]. Journal of Aircraft, 1999, 36(6): 999-1006.

2025年第25卷第13期

PDF下载

238

124

引用本文

BibTeX

文章信息

接收时间：2025-01-16
首发时间：2025-12-17
出版时间：2025-07-10

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2025-01-16

基金

作者信息

¹ 中国航发湖南动力机械研究所, 湖南株洲 412002

² 南京航空航天大学能源与动力学院, 南京 210016

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/kjhcy/CN/1208054447612208096

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT