叶型结冰计算及流场分析

2020-10-30任永鹏申连洋王忠义王艳华万雷王松

航空工程进展 2020年5期

任永鹏，申连洋，王忠义，王艳华，万雷，王松

(1. 哈尔滨工程大学动力与能源工程学院，哈尔滨 150001)

(2. 海装沈阳局驻哈尔滨地区第二军事代表室，哈尔滨 150001)

0 引言

飞机机翼、压气机进气部件结冰主要是由于过冷水滴(0 ℃以下仍为液态的水滴)撞击到固体表面冷凝，因此产生了结冰现象[1]。大气环境温度、液态水含量、过冷液滴直径和气流速度等因素都会对结冰的种类、结冰量和结冰区域产生影响[2-3]。翼型和叶型结冰都会改变其原有的气动外形，影响气动特性。因此，对于结冰机理的研究十分重要[4-5]。

20世纪60年代，科研人员开始对结冰问题开展数值模拟研究。21世纪随着计算机科学的迅猛发展，数值模拟已成为预测翼型结冰的重要手段之一。以计算流体力学和计算传热学为基础，流场、液滴撞击特性、结冰量的计算方法日趋成熟，数值模拟以其经济高效的优势，在结冰问题研究上得到了广泛的应用。

国外，美国代顿大学的C.MacArthur[6]建立了二维翼型上明冰和霜冰积累的数学模型，对含有液态水的结冰过程进行热力学分析，使液滴收集系数、传热和物质交换随着翼型表面积冰形状的改变而进行动态更新，从而对翼型周边的流场及液滴轨迹进行标准化的计算；美国通用电气全球研究中心的S.Saxena等[7]，模拟了液滴在动叶、静叶上的飞溅效应、叶片固结效应和叶片尾缘的蔓延效应，研究发现液滴与动、静叶间的相互作用会导致压气机喘振裕度的变化达到25%；S.Ozgen等[8]利用FORTRAN代码对二维NACA0012翼型积冰形状和二维Twin Otter翼型的水滴收集系数进行预测，计算结果与试验结果吻合较好。我国对于翼型结冰的研究虽起步较晚，但近些年来在结冰数值模拟方法上也积累了一定的科研成果，赵秋月[9]利用Fluent用户自定义设计，计算过冷水滴的运动轨迹、航空发动机三维进口支板的水收集系数、以及三维旋转航空发动机进口整流罩帽的水收集系数，并将数值计算结果与改进后的冰风洞运行相似准则结果进行对比，结果表明翼型的水收集系数和缩放后的翼型水收集系数吻合较好；张丽芬等[10]利用欧拉-拉格朗日法计算空气-液滴两相流，结合溢流水的溢流方向以及空气速度矢量，整合出一套计算三维翼型表面积冰的方法，利用该方法对NACA0012平直翼和截面为GLC-305的后掠翼翼型进行非稳态、非结构网格的数值模拟，其结果与试验值吻合较好。

目前对于结冰方面的研究主要集中在冰形预测上，对于结冰前后叶型气动特性的对比分析相对较少。虽然有少数研究人员分析了结冰前后翼型、叶型周围流场的变化，但缺少相关气动参数的分析[11-12]。

本文通过商业软件FENSAP-ICE模拟NACA0012翼型的结冰过程，对比NASA的结冰试验数据，验证结冰数值计算结果的正确性；对某型压气机进口导叶进行结冰计算，定量分析结冰前后叶片气动参数的变化情况，总结结冰对叶型气动特性的影响情况，以期为后续相关工作奠定基础。

1 结冰数值模拟方法

本文采用定常描述计算结冰，假定空气工质为理想气体，忽略质量力。利用微小颗粒偏微分方程计算液滴速度，采用欧拉法计算液滴撞击特性；利用Shallow-Water结冰模型，根据守恒定理对机翼表面的结冰量和结冰冰形进行求解。利用模块化思想，使用时间多步长法，将总的结冰过程分成N步进行计算，在每个步长内分别进行绕流流场、液滴撞击特性、热力学结冰的计算，单步完成后重新划分网格，继续计算直至总的结冰过程结束。

1.1 微小颗粒偏微分方程

FENSAP-ICE解决了微小颗粒偏微分方程对于颗粒速度和水浓度的问题，因此可以计算获得单次液滴喷射中的水浓度、液滴速度、液态水捕捉效率分布、撞击模式和整个区域的撞击极限，而不需要在喷射点处对液滴进行繁琐的迭代，从而加快了液滴撞击特性的计算效率。

1.2 Shallow-Water结冰模型

Shallow-Water结冰模型是一种基于表面水膜运动而建立的固体表面结冰模型，当液态水撞击到固体表面时形成一层薄膜，在气流的作用下，水膜会发生回流，加之热力学的作用，水膜会发生冻结、蒸发和升华。模型通过确定固体表面每个节点上的水膜厚度，对结冰时和液态水回流时的传热和传质问题进行计算。

Shallow-Water结冰模型假设水膜以液态水或积冰的形式附着在固体表面，在固体表面和法向上建立三维坐标系，并在该坐标系内构建水膜速度函数。同时，对水膜速度函数问题进行简化，假设水膜速度沿固体表面法线方向呈线性变化，在粘性壁面处速度为0。对水膜上的速度进行平均，得到水膜厚度与水膜平均速度间的函数关系。最终通过水膜运动偏微分方程，即水膜的质量守恒方程和能量守恒方程计算固体表面的结冰量和结冰冰形。

2 NACA0012翼型结冰数值模拟

2.1 模型及边界条件介绍

以结冰计算中常见的NACA0012翼型为计算模型，其翼型弦长0.533 4 m，在气流攻角4°，空气总压101 325 Pa，空气流速67 m/s，液态水含量1.0 g/m3，液滴颗粒直径20 μm，结冰时间6 min的情况下，给定翼型壁面为固体无滑移壁面，分别对翼型在环境温度为-2.22 ℃、-26.11 ℃时的结冰情况进行计算。其中，美国NASA的研究人员经过长期试验和计算研究，对比冰形发现以1 min为间隔进行时间多步长法数值计算最为合理[13]，因此使用时间多步长法计算的时间间隔也选取1 min。

2.2 网格划分

采用ICEM软件对模型进行网格划分，在翼型周围采用O型结构化网格，计算模型翼型前缘网格分布情况如图1所示，其中第一层网格高度0.003 mm，膨胀比1.1，保证网格y+值小于1。

图1 NACA0012翼型计算网格

在空气流速67 m/s，液态水含量1 g/m3，温度-2.22 ℃，结冰时间6 min，攻角4°，液滴颗粒直径20 μm的条件下，对计算模型进行网格无关性验证，网格数量对翼型结冰量的影响规律分布如图2所示，最终选择20万网格量进行计算。

图2 网格数量与翼型结冰量关系

2.3 数值计算结果分析

2.3.1 冰形计算结果

环境温度为-2.22 ℃、-26.11 ℃时翼型的结冰冰形计算结果如图3所示，其中c为弦长。

(a) -2.22 ℃

从图3可以看出：在其他条件一致的情况下，当环境温度为较高的-2.22 ℃时，由于撞击到翼型表面的液滴仅有一部分冻结为冰，另一部分受来流影响继续向翼型的上方或后方移动最终冻结，这一部分溢流水在翼型前缘形成了向上翘起的角冰，因此翼型在环境温度为-2.22 ℃时凝结成了典型的明冰；当环境温度为较低的-26.11 ℃时，由于来流液滴撞击到翼型表面后立即全部冻结，并无溢流水的形成，因此翼型在环境温度为-26.11 ℃时，冻结成为外形与翼型型线较为一致，形状较规则的毛冰。

2.3.2 计算结果误差分析

为考察数值模拟方法的准确性，本文利用商业软件FENSAP-ICE计算得到的翼型结冰冰形与现有的试验冰形、结冰代码预测冰形进行对比。试验冰形来源于NASA Lewis研究中心的结冰风洞试验台[14-15]。Lewis研究中心利用自行研发的2D LEWICE/IBL代码对NACA0012翼型结冰情况进行预测，计算数据丰富。为了精准地对翼型结冰计算冰形与试验冰形进行定量分析，获得结冰数值模拟方法的计算精度，引入三个无因次结冰特性量(无单位)对翼型结冰数值计算结果进行误差分析，如表1～表2所示。对于NACA0012这类对称翼型特征量[16]有三个。

表1 -2.22 ℃时翼型结冰特征量计算结果与误差分析

表2 -26.11 ℃时翼型结冰特征量计算结果与误差分析

无因次结冰面积ηs：

(1)

无因次结冰上极限ηLu：

(2)

无因次结冰下极限ηLd：

(3)

式中：A为翼型的面积，m2；c为翼型的弦长，m；S为翼型结冰面积，m2；Lu为结冰上极限，m；Ld为结冰下极限，m。

翼型结冰特征量示意图如图4所示。

图4 翼型结冰特征量

从表1～表2可以看出：在环境温度为-2.22 ℃和-26.11 ℃时，利用时间多步长法计算获得的冰形在无因次结冰面积、无因次结冰上极限、无因次结冰下极限这三个结冰特征量的误差整体小于由2D LEWICE/IBL代码获得的翼型结冰冰形误差，主要是由于FENSAP-ICE将气流粘性考虑在流场计算中，并利用Shallow-Water结冰模型进行冰形计算；结合翼型结冰冰形和结冰特征量误差分析来看，利用FENSAP-ICE时间多步长数值模拟方法可以获得与试验冰形吻合情况较好的明冰、毛冰冰形，能够清晰的反映出冰形的增长趋势和变化情况，计算结果处于研究所允许的误差范围内。因此，本文采用的数值模拟方法获得的明冰和毛冰冰形可用于后续研究。

3 压气机进口导叶结冰数值模拟

3.1 模型及边界条件

以某型号压气机进口导叶为模型，不考虑叶片根部、冠部区域以及相邻叶片间的影响，计算压气机进口导叶50%叶高处的结冰情况。其边界条件为：气流攻角0°，空气总压101 325 Pa，空气流速100 m/s，液态水含量0.25 g/m3，液滴颗粒直径10 μm，结冰时间60 s，多步长计算时间间隔20 s，叶型壁面为固体无滑移壁面；环境总温分别为265 K和276 K。

3.2 网格划分

使用ICEM软件对模型进行网格划分，在叶型周围采用O型结构化网格，进口导叶叶型前缘、尾缘网格分布情况如图5所示。

(a) 叶型前缘网格分布

在叶型计算网格中，第一层网格高度0.001 mm，膨胀比1.1。保证网格y+值小于1，在空气流速100 m/s，液态水含量0.25 g/m3，温度276 K，结冰时间60 s，攻角0°，液滴颗粒直径10 μm的条件下，对计算模型进行网格无关性验证，网格数量对叶型结冰量的影响规律分布如图6所示，本文选择20万网格量进行计算。

图6 网格数量与叶型结冰量关系

3.3 数值计算结果分析

3.3.1 冰形计算结果

环境温度为265、276 K时叶型结冰冰形的计算结果如图7所示。

(a) 265 K叶型结冰冰形

从图7可以看出：当环境温度为较低的265 K时，叶片前缘形成了与叶片形线较为一致的毛冰；当环境温度为较高的276 K时，叶片前缘形成了向上翘起的明冰。这与NACA0012对称翼型的结冰情况相同，对于具有非对称结构的进口导叶来说，除了叶片前缘产生积冰外，在叶片压力面上也发生了结冰现象。但压力面上的积冰厚度不大，且冰形与叶片压力面型线基本保持一致。

3.3.2 结冰前后流场分析

采用FENSAP-ICE计算叶型结冰，该软件会自动生成由于结冰而发生位移后的叶型计算网格，如图8所示。

(a) 265 K叶型前缘 (b) 265 K叶型尾缘

利用商业软件Fluent，设置与初始流场计算相同的边界条件，对结冰后叶型进行流场计算。结冰前后叶型的流场情况如图9所示，可以看出：即使气流以0°攻角流向叶型，但叶型的滞止区域仍发生在叶型的前缘偏下部，主要是由于气流流过叶片吸力面的速度要大于气流流过叶片压力面的速度，因此在叶型前缘处生成了一个下大上小的压力梯度，导致气流以一定角度流向叶型，在结冰前后叶型的尾部上侧都存在一个低速区，这个低速区内存在一个顺时针的大面积分离涡和一个逆时针的小面积分离涡，其中顺时针的分离涡主要是由于吸力面上发生的气流分离造成的，而较小面积的逆时针分离涡是由于顺时针分离涡与主流发生掺混而形成的；当环境温度为265 K时，结冰后叶型仅在前缘上部形成了一片较小的气流分离区，但并没有影响叶片周围流场的分布情况，结冰前后叶型的气流分离都发生在叶片吸力面侧，起始于弦长的50%处，毛冰的形成并没有对叶片周围流场产生较大影响；当环境温度在276 K时，由于向上翘起的角冰，导致结冰后叶型的气流分离加剧，分离点提前至叶片弦长的1%处，产生了一个面积约为原来2倍的低速分离区，低速区内的两个分离涡面积也随之增大，导致叶型的气动性能急剧恶化。

(a) 265 K结冰前

结冰前后叶型出口处，即x=0.06 mm处气流速度大小分布如图10所示。当环境温度为265 K时，结冰前后速度值变化不大。此时，低速区主要集中在-0.01～0.0 025 mm，这与图9(a)，(b)所展现的马赫云图分布一致。但当环境温度为276 K时，叶型出口处的最低速度值下降了约80%，低速区范围扩大了40%。同时，当环境温度为276 K时，结冰后叶型尾部低速区两侧的气流速度明显高于结冰前，可能是由于结冰导致的气流分离加剧，涡流流速增加，掺混主流后造成的。此外，结冰后尾流低速区中部存在一个速度值先增大后减小的过程，而结冰前却没有这个过程，这主要是叶型前缘结冰，使得尾部分离区沿弦向扩大，分离区中两个方向不同的分离涡在出口处共同作用而导致的气流速度值增大。综上，结冰会明显增大叶型出口速度不均匀性分析。

(a) 265 K结冰前后叶型出口处速度分布

3.3.3 结冰前后总压损失系数分析

叶型总压损失系数ω的计算公式为

(4)

式中:PT1为叶型进口，即x=-0.02 mm处气流总压；PT2为叶型进口气流动压；PD1为叶型出口，即x=0.06 mm处气流总压。

不同环境温度下叶型总压损失系数分布如图11所示。

(a) 265 K结冰前后叶型总压损失系数分布

(b) 276 K结冰前后叶型总压损失系数分布

从图11(a)可以看出：环境温度为265 K时，在y=-0.00 6 mm处总压损失系数最大，整体来看，结冰后叶型的总压损失系数有所增大，但增量极小。因此，叶型在265 K冻结形成的毛冰并不会额外产生较大的总压损失。从图11(b)可以看出：当环境温度为276 K时，结冰后总压损失系数明显升高，在y=-0.006 mm处总压损失系数达到最大值，较结冰前升高了约22%，随后总压损失系数呈现先下降，后上升再下降的趋势，主要是由于出口处速度先上升后下降再上升而导致的。还可以看出：无论环境温度高低，结冰后叶型压力损失系数的增大主要集中在尾流分离区的上部，这表明，结冰后叶型的压力损失大小主要由分离区上部的顺时针分离涡影响。叶型前缘结冰上翘角度越大，叶型尾部分离区域面积越大，分离涡越强烈，导致结冰后叶型压力损失越大。