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

2022-12-25曾飞雄毛汉冬章儒宸

南京航空航天大学学报 2022年6期

曾飞雄，毛汉冬，章儒宸

（上海飞机设计研究院，上海 201210）

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

国外从20 个世纪四五十年代就开始采用数值计算的方法研究飞机结冰防护问题。目前国际上通用的防冰传热耦合分析的基础理论均基于Messinger 模型［1］，并以此来建立防冰表面的传热传质方程。至今数个发达国家都开发了用于飞机结冰预测和防冰系统设计的计算软件，如NASA和波音的LEWICE 软件，空客公司的ONERA 软件，庞巴迪公司的CANICE 软件，罗罗公司的ICECREMO 软件，ANSYS 公司的FENSAP⁃ICE软件等。其中不少软件被SAE⁃ARP5903 收录为成熟软件进行推荐［2］。

目前国内针对机翼防冰仿真计算的研究较多，针对短舱防冰相对较少［3⁃6］。使用较多的商用软件为FENSAP，也有部分高校使用用户自定义函数（User defined function，UDF）自编程求解，通常算法是内外流场分别计算，再利用距离反比插值法进行内外流场数据交互耦合［3］，对流换热系数的计算则通常采用附面层积分法［7⁃8］或CFD 求解［9］。民机设计的工程应用重视计算效率和用户自定义，因此本团队基于基础软件工具自研开发了防冰计算软件，在经典算法上进行优化，从工程化角度出发，结合国内民机的设计研发和完整的适航取证经验积累，从条款提炼计算判据［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）所示。

式中：ṁimp为撞击水量，ṁin，n为微元体第n个边的流入水量，ṁice为结冰量，ṁevap为蒸发水量，ṁout，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个最近点的距离进行加权计算获得该节点数据。

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

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

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

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

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

4 对流换热系数计算

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

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

目前很多商业仿真软件可以进行复杂流动的CFD 求解，其获得的是包含了各种换热因素的总换热系数，但求解速度通常较慢。附面层积分法普遍用于结冰防冰求解领域，求得边界层外的速度、温度分布后，利用相应流态的积分方程得到对流换热系数的近似解［7］，但在求解积分方程时，通常建立的是沿流线的二维计算模型，对于转捩和流态的判定准确性仍需要进一步验证。

短舱唇口外形与圆管特征类似，考虑到外掠圆管的原理性试验已经有充分积累，本文基于外掠圆管的修正经验公式［15］，开展三维对流换热系数的求解分析。

式中：Nu为努赛尔数，C为试验修正系数，Re为雷诺数，m为雷诺数修正因子，Pr为普朗特数。修正因子见表1。

表1 外掠圆管流动修正因子Table 1 Flow correction factor of swept tube

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