民用飞机短舱防冰性能分析建模与仿真

曾飞雄，毛汉冬，章儒宸; ZENG Feixiong; MAO Handong; ZHANG Ruchen

2025年6月16日 15:15 星期一

网刊加载中。。。

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

确定继续浏览么?

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

摘要

基于民用飞机短舱防冰系统的特点，结合发动机结冰防护要求和适航条款的相关要求，阐述了短舱防冰性能分析的工程计算流程，并给出了提升计算效率的方法。在计算需求中重点分析归纳了短舱防冰的计算判据和常用计算工况，在计算模型中从外流场、内流场和耦合面3个域开展能量平衡和质量平衡分析，对计算网格提出了工程化的简化要求和方法；采用基于K近邻搜索的反距离加权插值算法，相比传统穷举法插值大幅提升了求解效率，且该方法具备网格无关性的特征；相较于经典的附面层积分法求解对流换热系数，本仿真模型在综合考虑三维效应和求解效率后，采用了外掠圆管的对流换热系数修正经验公式，并为后续模型修正提供了接口。基于自主搭建的民机防冰性能分析模型，开展了已有型号的短舱防冰性能分析，计算结果与试飞数据的表面温度对比表明，仿真中的优化满足工程分析需求且极大地提升了求解效率。

关键词

短舱防冰; 性能分析; KD树插值; 对流换热系数; 试飞对比

短舱防冰是指发动机进气道唇口前缘的防冰方式，设计指标包括干蒸发和湿蒸发。短舱防冰的方式主要分为电防冰、热气防冰。两种防护方式在仿真建模时，主要在内部热流分析时存在差异。

国外从20个世纪四五十年代就开始采用数值计算的方法研究飞机结冰防护问题。目前国际上通用的防冰传热耦合分析的基础理论均基于Messinger模型^［

1］，并以此来建立防冰表面的传热传质方程。至今数个发达国家都开发了用于飞机结冰预测和防冰系统设计的计算软件，如NASA和波音的LEWICE软件，空客公司的ONERA软件，庞巴迪公司的CANICE软件，罗罗公司的ICECREMO软件，ANSYS公司的FENSAP⁃ICE软件等。其中不少软件被SAE⁃ARP5903收录为成熟软件进行推荐^{［参考文献 2

百度学术}2］。

目前国内针对机翼防冰仿真计算的研究较多，针对短舱防冰相对较少^［

3⁃6］。使用较多的商用软件为FENSAP，也有部分高校使用用户自定义函数（User defined function，UDF）自编程求解，通常算法是内外流场分别计算，再利用距离反比插值法进行内外流场数据交互耦合^{［参考文献 3

百度学术}3］，对流换热系数的计算则通常采用附面层积分法^{［参考文献 7⁃8}7⁃8］或CFD求解^{［参考文献 9

百度学术}9］。民机设计的工程应用重视计算效率和用户自定义，因此本团队基于基础软件工具自研开发了防冰计算软件，在经典算法上进行优化，从工程化角度出发，结合国内民机的设计研发和完整的适航取证经验积累，从条款提炼计算判据^{［参考文献 10

百度学术}10］，优化了对流换热系数算法，并创新性地使用了KD树（K⁃dimensional tree）数据插值算法，在确保计算精度的同时大幅提升了计算效率。

1 性能分析关键指标

针对发动机和飞机对短舱防冰的功能需求，短舱防冰能力性能分析的考核判据主要包括：（1）各飞行工况下的唇口内表面吞冰量均小于发动机可接受吞冰量；（2）各飞行工况下的唇口外表面冰脱落量小于机体可承受的冰积聚尺寸。

唇口内部冰脱落的影响通常大于外部冰脱落，因此主要性能考核指标侧重于进气道内表面的结冰情况分析。由于冰型尺寸通常难以在试验中进行复现验证，因此进行短舱防冰性能分析时，通常也把表面温度作为间接考核的判据，以此表明计算模型的准确性。

进气道内表面吞冰量通常转换成某站位的结冰厚度进行判定。通过仿真分析稳态后流水量及持续时间得到溢流冰量，再参照FAR33.77^［

10］定义的发动机名义结冰尺寸，得到某站位的结冰厚度。

短舱防冰的热平衡与外界环境、发动机抽吸量、飞行状态、引气状态等相关。因此性能分析时，需要考量不同飞行工况作为计算状态点。通常均需考虑单发起飞、正常起飞、爬升、巡航、下降、待机、进场以及着陆等飞行状态。

计算吞冰量时需要用到的持续时间基于CCAR25部附录C中的定义，通常考虑连续最大结冰条件下17.4海里（1海里=1.852 km）的水平穿云时间，间断最大结冰条件下2.6海里的水平穿云时间。通常较为严酷的工况是连续最大结冰条件下45 min的待机飞行工况以及下降小推力工况。

2 总体计算模型

短舱防冰的性能计算模型重点是建立内部热流与外部热流的平衡关系。计算时需要考虑4个计算区域，包括外部流场区域（空气‑水滴两相流）、表面水膜区域、蒙皮结构区域以及防冰腔内部流场区域。为了加快计算速度，同时保证足够可靠的计算精度，本文热流耦合计算模型以独立的蒙皮网格为计算域，调用包含全部网格的完整内外流场结果作为初始值，在后续平衡迭代计算中对外流场中离蒙皮一定距离以外的局部热流场进行更新，对于内部流场进行少量的更新。整个计算流程如图1所示。

图1 防冰计算迭代流程图

Fig.1 Iteration process for anti-ice analysis

能量平衡方程主要考虑外部的导热项、水滴显热项、水滴潜热项、对流换热项、气动加热项和内部的对流换热项^［

1］。质量平衡主要考虑每个微元体中各项水量的平衡，如式（1）所示。

{\dot{m}}_{i m p} + \sum_{n} {\dot{m}}_{i n, n} = {\dot{m}}_{i c e} + {\dot{m}}_{e v a p} + \sum_{n} {\dot{m}}_{o u t, n}

（1）

式中： ${\dot{m}}_{i m p}$ 为撞击水量， ${\dot{m}}_{i n, n}$ 为微元体第n个边的流入水量， ${\dot{m}}_{i c e}$ 为结冰量， ${\dot{m}}_{e v a p}$ 为蒸发水量， ${\dot{m}}_{o u t, n}$ 为微元体第n个边的流出水量。防冰计算时，若表面温度大于0 ℃，可将结冰量设置为0。

在进行计算域网格绘制时，常使用半模机体绘制外流场网格。需要将机身、机翼、吊挂、发动机包含在内。针对短舱防冰计算分析的特点，需要特别地对短舱区域进行网格加密，并将唇口耦合面单独绘制网格，便于后续进行内外流场耦合时确定计算域。考虑到发动机的抽吸效应，需要在发动机建模时单独设置压力出口和压力入口边界，用于模拟发动机的抽吸气量（图2）。

图2 发动机边界设置图

Fig.2 Engine boundary condition setting

防冰腔体的内流场根据防冰形式主要分为基于笛形管的小孔射流和基于直流喷嘴的环状流动。喷孔处需要映射加密以提高内流场的求解精度（图3）。

图3 防冰腔内流场加密效果示意

Fig.3 Inner flow field mesh encryption

3 KD插值算法

在防冰计算中涉及外部流场网格、蒙皮结构网格与内部流场网格3套不同的网格，因此热耦合计算时网格交界面间的数据需要通过插值计算进行映射传递。在一个耦合计算循环中需要进行两次界面插值计算，分别为：

（1）蒙皮结构表面网格至内外流场壁面边界网格的热流和温度插值。

（2）内外流场壁面边界网格至蒙皮结构表面网格的热流和温度插值。

两次插值计算在热耦合计算流程中进行次序如图4所示。

图4 防冰热耦合计算流程

Fig.4 Anti-ice thermal-coupling process

由于耦合计算的每一次单个迭代步内都需要进行两次界面插值，计算至最终稳态结果需要进行多次循环迭代，并且边界网格量级较大，因此，插值是整个计算中最耗资源的过程。另外，内外流场网格与蒙皮结构网格可能采用不同的拓扑，为使插值算法具有较好的适应性，插值计算应基于网格节点进行，从而使插值计算过程与网格单元类型无关。为在足够的插值精度下尽可能提高插值计算速度，同时保证插值算法能够适用于各类不同网格形状（网格无关性），本文采用基于K近邻搜索的反距离加权插值算法对界面数据进行插值映射，该算法目前在计算机前沿的机器学习及大数据处理领域已有应用^［

11⁃14］。

耦合求解器基于K近邻搜索的反距离加权插值方法对界面数据进行插值映射，该算法基本思想如图5所示。

图5 基于K近邻搜索的节点数据反距离加权插值

Fig.5 K-nearest neighbor search for inverse distance weighted interpolation

图5中的网格1为已获得相应数据的表面网格，网格2为待插值获取数据的表面网格，为了得到网格2中某一节点的数值，先对该节点周围N个距离最近的节点进行查找，然后基于到这N个最近点的距离进行加权计算获得该节点数据。

通常采用的反距离加权函数形式为

u (x) = \sum_{i = 0}^{N} \frac{w_{i} (x) u_{i}}{\sum_{j = 0}^{N} w_{j} (x)}

（2）

式中：u为样本点数据，N为临近点数目，w为权重量。第i个临近点的权重w表达式为

w_{i} (x) = \frac{1}{d (x, x_{i})^{p}}

（3）

式中：d为距离，p为幂参数。当p=2时，即为以空间距离（欧式距离）的倒数为权重。可以看出距离越近的临近点对样本点的影响越大。

不同于常见的穷举搜索方法，K近邻搜索算法实现上采用KD树结构进行快速二分匹配，其算法时间复杂度为O（lg（2N）），与穷举搜索时的O（N²）时间复杂度相比，其计算速度将极大提高，且数据量越大时优势越明显。相较于其他插值计算方法，基于K近邻搜索的反距离加权插值具有两个优点：

（1）计算速度快。由于K近邻搜索具有较低的时间复杂度，且对于空间位置固定的表面网格节点拓扑，其KD树只需生成一次，因此在进行K近邻搜索和逆距离计算时计算开销较小，当界面节点数大于100时，使用K近邻插值的计算速度将比穷举搜索插值快1 000倍以上。

（2）适用于各类不同外形与网格间的数据映射。插值算法基于网格节点进行，与具体网格类型和拓扑无关，因此具备几何无关性与网格无关性，适用于各类不同交界面网格间的数据映射。

4 对流换热系数计算

对流换热是影响防冰热载荷的重要因素，基于工程经验，外表面对流换热在某些工况下能占外部总载荷的50%以上。

而对流换热系数是数值求解对流换热项的关键参数。从已有的研究来看，主要分为3类求解方式：（1）计算流体力学（CFD）数值求解；（2）附面层积分法求解；（3）基于试验修正的经验公式求解。

目前很多商业仿真软件可以进行复杂流动的CFD求解，其获得的是包含了各种换热因素的总换热系数，但求解速度通常较慢。附面层积分法普遍用于结冰防冰求解领域，求得边界层外的速度、温度分布后，利用相应流态的积分方程得到对流换热系数的近似解^［

7］，但在求解积分方程时，通常建立的是沿流线的二维计算模型，对于转捩和流态的判定准确性仍需要进一步验证。

短舱唇口外形与圆管特征类似，考虑到外掠圆管的原理性试验已经有充分积累，本文基于外掠圆管的修正经验公式^［

15］，开展三维对流换热系数的求解分析。

N u = C R e^{m} P r^{1 / 3}

（4）

式中： $N u$ 为努赛尔数， $C$ 为试验修正系数， $R e$ 为雷诺数， $m$ 为雷诺数修正因子， $P r$ 为普朗特数。修正因子见表1。

表 1 外掠圆管流动修正因子

Table 1 Flow correction factor of swept tube

Re	C	m
0.4~4	0.989	0.330
4~40	0.911	0.385
40~4 000	0.683	0.466
4 000~40 000	0.193	0.618
40 000~400 000	0.027	0.805

式（4）中的相关无量纲数基于外流场计算结果中当地网格内的速度、密度、黏度等参数进行求解，雷诺数的结果会影响到试验修正系数 $C$ 和雷诺数修正因子 $m$ ，由于不同机型短舱的前缘外形各有差别，会对特征尺寸有影响，因此这两个修正系数也可基于风洞或干空气试飞的温度结果进行算法修正，极大地提升了结冰条件下仿真计算的准确性和一致性。

5 仿真求解与结果分析

本文使用的算例是参照某型号飞机的自然结冰试飞数据。计算网格模型基于飞机真实构型进行绘制。由于试飞过程只能采集唇口区域蒙皮温度，而无法获取结冰量等数据，因此本仿真求解主要对比计算与试验的温度差异情况。试飞状态参数如表2所示。

表2 计算算例工况点

Table 2 Simulation cases

序

号

海拔/

环境

温度/

℃

马赫数

液态

水含量/

（g·m^-3）

平均体积直径/

μm

进气量/

（kg·s^-1）

2 400

-6.80

0.35

0.83

25.4

113

1 430

-11.24

0.36

0.71

16.8

236

空气外流场计算的主要边界条件类型为：计算区域的外围定为压力远场边界；发动机进气道入口为指定流量的压力出口，其压力值为达到目标流量时所匹配出的压力；发动机尾喷出口为指定流量的压力入口。

防冰腔内流场的主要边界条件类型为：笛形管喷孔处设定为压力入口，通过一维管路流动计算得到对应的入口总温和压力；笛形管表面设为恒温壁面；排气孔出口处设定为压力出口，压力和温度取为大气环境条件。

以发动机进气道唇口前缘蒙皮为耦合面，通过插值获取外流场和内流场的相关参数，基于蒙皮表面的离散单元进行能量和溢流水质量平衡求解。图6、7分别显示的是算例1和2的唇口表面温度分区云图。从图中可以看出，有显著热斑的区域集中在内表面，高温区域出现在180°（下方）区域，这些均与设计目标相吻合。

图6 算例1表面温度分布云图

Fig.6 Skin temperature contour in Case 1

图7 算例2表面温度分布云图

Fig.7 Skin temperature contour in Case 2

基于试飞改装安装的表面温度传感器，对进气道某截面的温度数据进行对比验证。验证结果分别如图8、9所示。

图8 算例1某截面表面温度对比图

Fig.8 Section temperature comparison in Case 1

图9 算例2某截面表面温度对比图

Fig.9 Section temperature comparison in Case 2

仿真结果表明，截面表面温度计算值与试飞实测值趋势吻合较好，峰值位置较为匹配，偏差较大处出现在外唇口后缘，总体测点温度偏差在10 ℃以内，在工程领域属于可接受范围。

6 结　　论

（1）给出了基于工程应用的民用飞机短舱防冰系统性能分析的总体建模与仿真方法。

（2）采用KD树插值并使用基于修正经验公式获取了对流换热系数，在保证计算准确性的同时大幅提升求解效率。

（3）对某型号的短舱防冰性能进行了分析，计算结果与试飞结果吻合较好，满足工程应用需求，具有很好的工程实践价值。

参考文献

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. [百度学术]

ARP5903. Droplet impingement and ice accretion computer codes［S］. ［S.l.］： SAE， 2003. [百度学术]

郁嘉，赵柏阳，卜雪琴，等. 某型飞机发动机短舱热气防冰系统性能数值模拟［J］. 空气动力学学报， 2016，34（3）： 302-307. [百度学术]

YU Jia， ZHAO Boyang， BU Xueqin， et al. Numerical simulation of the performance of an engine nacelle hot-air anti-icing system［J］. Acta Aerodynamica Sinica， 2016， 34（3）： 302-307. [百度学术]

常士楠，艾素霄. 飞机发动机进气道防冰系统的设计计算［J］.北京航空航天大学学报，2007，33（6）：649-652. [百度学术]

CHANG Shinan， AI Suxiao. Design and calculation for the anti-icing system of an aircraft engine inlet［J］. Journal of Beijing University of Aeronautics and Astronautics，2007，33（6）： 649-652. [百度学术]

朱永峰，方玉峰，封文春. 某型飞机发动机短舱防冰系统设计计算［J］.航空动力学报，2012，27（6）： 1326-1331. [百度学术]

ZHU Yongfeng， FANG Yufeng， FENG Wenchun. Design and calculation of aircraft nacelle anti-icing system［J］. Journal of Aerospace Power， 2012， 27（6）：1326-1331. [百度学术]

陈景松.发动机进气道前缘防冰腔性能研究［D］.南京：南京航空航天大学，2008. [百度学术]

CHEN Jingsong. Research on engine inlet anti-icing performance［D］. Nanjing：Nanjing University of Aeronautics and Astronautics， 2008. [百度学术]

宋馨，林贵平，卜雪琴. 防冰表面的对流换热计算分析［J］. 北京航空航天大学学报， 2011， 37（1）： 49-53. [百度学术]

SONG Xin，LIN Guiping，BU Xueqin. Convective heat transfer analysis of anti-icing surface［J］.Journal of Beijing University of Aeronautics and Astronautics，2011， 37（1）： 49-53. [百度学术]

刘重阳. 飞机电热防/除冰系统数值计算及实验研究［D］. 南京：南京航空航天大学，2020. [百度学术]

LIU Chongyang. Simulation and test research on aircraft electrical anti-icing/de-icing system［D］. Nanjing：Nanjing University of Aeronautics and Astronautics，2020. [百度学术]

徐佳佳，史献林，王向转. 民用飞机风挡表面对流换热系数的校核［J］. 民用飞机设计与研究， 2015（1）：87-90. [百度学术]

XU Jiajia， SHI Xianlin， WANG Xiangzhuan. Modification method of windshield HTC for civil aircraft［J］. Civil Aircraft Design & Research，2015（1）： 87-90. [百度学术]

FAA. Airworthiness standards： Aircraft engines： 14 CFR PART 33［S］. Washington DC： FAA，2022. [百度学术]

毋雪雁，王水花，张煜东.K最近邻算法理论与应用综述［J］. 计算机工程与应用， 2017， 53（21）： 1-7. [百度学术]

WU Xueyan， WANG Shuihua， ZHANG Yudong. Survey on theory and application of K-nearest neighbor algorithm［J］. Computer Engineering and Applications，2017， 53（21）： 1-7. [百度学术]

HASSANAT A B A. Furthest-pair-based binary search tree for speeding big data classification using K-nearest neighbors［J］.Big Data，2018，6（3）： 225-235. [百度学术]

赵俭辉，龙成江，丁乙华，等.一种基于立方体小栅格的K邻域快速搜索算法［J］，武汉大学学报（信息科学版），2009， 34（5）： 615-618. [百度学术]

ZHAO Jianhui， LONG Chengjiang， DING Yihua， et al. A new K-nearest neighbors search algorithm based on 3D cell grids［J］. Geomatics and Information Science of Wuhan University， 2009， 34（5）： 615-618. [百度学术]

Zhang Haida，Xing Zhihao，Chen Lu， et al. Efficient metric all-K-nearest-neighbor search on datasets without any index［J］. Journal of Computer Science and Technology， 2016， 31（6）： 1194-1211. [百度学术]

THEODORE L. Bergman， fundamentals of heat and mass transfer［M］. 7th Edition. ［S.l.］： John Wiley & Sons， 2011： 458. [百度学术]

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

民用飞机短舱防冰性能分析建模与仿真 PDF

摘要

关键词