跨音速风扇叶片高低转速下混合相积冰的对比研究

任虎虎，徐强仁，王立志，李广超，赵巍,4，赵庆军,4,5

（1. 沈阳航空航天大学 航空发动机学院，沈阳 110136；2. 中国科学院工程热物理研究所，北京 100190；3. 中国科学院轻型动力重点实验室，北京 100190；4. 中国科学院大学航空宇航学院，北京 100049；5. 中国科学院工程热物理研究所分布式冷热电联供系统北京市重点实验室，北京 100190）

飞机积冰分为地面积冰与飞行积冰，其中飞行积冰占比达到90%以上。在飞行过程中，飞机迎风面会发生积冰，包括机翼、旋转帽罩、旋翼、风扇叶片等部件[1-4]。目前，对机翼、旋转帽罩、旋翼等已开展了大量研究，对风扇叶片的研究还比较少，本文针对风扇叶片，展开积冰研究。相较于机翼等静止部件，风扇叶片积冰的形成及脱落均会受到离心力的作用。风扇叶片积冰会改变叶片型面，破坏流场，甚至使压气机流道堵塞，造成压气机喘振，风扇叶片上形成的积冰在脱落后可能进入压气机，与压气机叶片产生碰撞，威胁飞行安全。积冰研究有工程估算、实验研究和数值模拟等方法，由于数值模拟有成本低、周期短、流场信息丰富等优势而快速发展。

目前，旋转部件积冰研究主要集中在旋转帽罩及旋翼2 部件。Blez 等[5]通过一套观察发动机整流帽罩与风扇叶片的视频成像系统，给出了不同时刻的积冰图像。梁鹏[6]对旋转锥体进行了过冷水滴积冰研究，发现锥尖至底部局部水收集系数呈减小趋势，随着转速的提高，水滴受到的离心力增大，惯性力减小，碰撞到锥体的趋势减弱。Li 等[2]对带有不同形状帽罩的风扇叶片进行了积冰模拟，发现明冰条件的帽罩积冰与几何形状有较大关联，风扇叶片上的积冰相对不受其几何形状的影响。Dong 等[7]对旋转帽罩的研究表明，回流水对积冰有显著影响，转速对积冰的形状略有影响，固定圆锥上的积冰厚度大于旋转圆锥上的积冰厚度。Villedieu 等[8]对冰晶撞击过程中的粘附、破碎反弹等建立了相应的数学模型。Trontin 等[9-10]通过加入冰晶侵蚀的影响，改进了冰晶粘附效率模型，并对不同翼型进行了过冷大水滴、冰晶和混合相积冰仿真。Aliaga 等[11]以旋转的NACA0015 作为研究对象，将冰的增长与空气-液滴两相流耦合，使用动态缝合网格提高了积冰的预测精度。Liu 等[3]以某型无人机螺旋桨为研究对象，对其表面结冰过程进行了试验研究，发现在气动力和离心力的共同作用下，叶片表面未冻结的水从根部流向叶片尖端，形成了十分复杂的冰柱。陈宁立等[12]以NACA0015 作为研究对象，发展了一种适用于预测旋转表面积冰的数学模型。研究结果表明，离心力使流出控制体的水膜质量增加，导致驻点附近的积冰厚度减小。Baumert 等[13]使用NACA0012 翼型及圆柱模型，根据冰风洞实验数据完善了冰晶积冰计算软件，补充了冰晶积冰的实验与模拟数据。郭琪磊等[14]以NACA0012 翼型为对象，研究了混合相态的积冰模拟，发现混合相态积冰达到最大厚度，需要有充足的冰晶与液态水含量，增加飞行速度会使前缘驻点的积冰速率与积冰量增加。卜雪琴等[15]使用NACA0012 翼型计算了霜冰和明冰条件下的积冰情况，结果表明，冰晶黏附效应对混合相结冰量及冰形有很大影响。Farag 等[16]通过研制螺旋桨液滴撞击的BETAPROP 程序，发现部件的几何形状对水收集系数和撞击极限产生很大的影响。陈希等[17]通过对直升机旋翼的研究，证明了离心力在结冰数值模拟中的重要性。Chen 等[18]建立了旋翼的积冰数学模型，研究发现，转子积冰受到离心力和水膜运动的影响，考虑到离心力与水膜运动后，叶片上的积冰会减少。Chen 等[19]通过对旋翼的研究，验证了其建立的积冰预测模型。研究表明，转子转速增加，积冰厚度增加。

风扇叶片的混合相积冰研究主要为冰晶轨迹的研究。吴涛[4]以Rotor-67 跨音速转子为研究对象，研究其设计转速积冰，发现过冷水滴的撞击区域主要为叶片吸力面，随着过冷水滴粒径的增大，过冷水滴受到的惯性变大，运动方向不易改变，易沿着原始速度方向进行直线运动。Das 等[20]对高旁通涡扇发动机的压气机转子进行了数值模拟，研究表明，在较低的转子转速下，积冰从叶根到叶尖厚度逐渐减小，转子的转速增加，叶根到叶尖积冰的厚度变化减小，同一截面吸力面积冰逐渐消失，压力面积冰向前缘移动。Norde 等[21]针对涡扇发动机的冰晶结冰，考虑冰晶传热和相变对轨迹及撞击结冰表面的影响，开发了一种积冰方法，使用NACA0012 翼型进行了验证，结果表明，该方法能较准确地预测冰厚与冰形。姜飞飞等[22]对冰晶相变传热传质方程进行了离散化，计算了冰晶在涡扇发动机内涵通道内运动过程中的粒子半径、冰晶温度与速度、冰水混合粒子中的液态水质量分数等，获得了冰晶在低压压气机内涵道的运动轨迹和叶片碰撞特性。Zhang 等[23]采用拉格朗日法测定了冰晶和液滴的运动轨迹，考虑了冰晶的破碎和反弹，以及薄膜的飞溅，建立了冰晶撞击模型，提出了混合相积冰热力学模型。通过不同混合相条件下NACA0012翼型的积冰模拟，验证了混合相结冰方法的有效性和合理性。综上所述，积冰时考虑离心力得到的积冰厚度减小，冰晶的相变换热及侵蚀影响积冰冰形，考虑冰晶的粘附效应将使积冰厚度减小，但在风扇叶片混合相积冰中，转速对风扇叶片过冷水滴与冰晶积冰的变化规律考虑较少。

本文对风扇叶片进行混合相积冰，分析过冷水滴与冰晶在流道内的变化规律，对比风扇叶片不同转速时过冷水滴与冰晶对叶片的撞击角度，获得风扇叶片高低转速时积冰差异的形成原因，为预测风扇叶片在高低转速工作状态下的积冰提供理论支持。

1 研究方法与验证

积冰的数值模拟一般主要由3 部分组成。第一部分是对部件周围的空气-粒子（过冷水滴与冰晶）流场进行模拟，此部分为部件表面的积冰模拟提供空气剪切力。第二部分是计算过冷水滴与冰晶的撞击特性，获得粒子在流场中的速度、蒸汽与过冷水滴和冰晶之间的质量和能量传输、局部水收集系数和局部冰晶收集系数。第三部分是对已获得的变量采用合适的模型计算积冰，通过迭代得到对应时间部件表面的积冰量。本文采用文献[24]所发展基于牛顿剪切定律描述水膜流动的SWIM（Shallow Water Icing Model）积冰模型，该模型对叶轮机械的研究中，增加了离心力与科氏力对积冰的影响，积冰求解时的交互关系如图1 所示。

依据文献[25]可知，Flow Solver ALE（FENSAP）模块用于计算气流流场，可求解定常和非定常可压缩三维Navier-Stokes 方程。流体可以是无黏性的，也可以是黏性的，流动可以是层流或湍流，湍流由单方程或双方程模型模拟。通过对固体表面上能量方程再求解，可以以二阶精度直接计算出对明冰结冰至关重要的壁面热通量。由文献[25]可知，气流流场计算也可使用Fluent、CFX 等软件进行计算。DROP3D 是FENSAP-ICE 系统过冷水滴或冰晶三维欧拉单次撞击模块。DROP3D 可以接受FENSAP、Fluent 和CFX流动求解数据，它可以求解粒子速度和液态水浓度的偏微分方程，处理外部和内部流的冲击。因此，DROP3D 可以在一次粒子注入中提供整个区域的液态水浓度、液滴速度矢量、液态水捕获效率分布、撞击方式和撞击极限，而无需在注入点上进行复杂的迭代过程。进入流道的过冷水滴与冰晶逐渐升温，其升温速度与气流升温速度不同，通过粒子热平衡方程，可以计算不同位置处过冷水滴与冰晶的温度，过冷水滴和冰晶与气流之间的能量转移（包括气流对流、辐射能的吸收与发射、质量耦合效应引起的蒸发、凝结、冻结以及融化的增加与损失）。SWIM 是ICE3D 模块使用的积冰模型，是基于积冰形成的复杂热力学偏微分方程，它可以在复杂三维表面上求解冰形、水膜厚度和表面温度。在计算中，使用NTI 弹跳模型判断晶体粘附。该模型中，冰晶的粘附数量由冰晶撞击速度、晶体尺寸和薄膜高度确定。

本文通过对风扇叶片进行冰风洞试验验证积冰位置的准确性。冰风洞试验在东北农业大学冰风洞展开，该冰风洞出口截面为1 m×1 m，最大来流速度为19 m/s，液滴粒径为15～50 μm，液态水质量浓度为0.3～3 g/m3，来流温度不可调节，需依靠来流温度进行调整。试验总温为268.15 K，液滴粒径为20 μm，液态水质量浓度为0.45 g/m3，积冰时间为1 350 s，详细试验条件见表1。

表1 积冰试验条件Tab.1 Icing test conditions

试验积冰位置结果如图2a 所示，FENSAP-ICE计算积冰位置如2b 所示。依据积冰在弦向与沿叶高方向的积冰位置变化进行初步验证。积冰位置主要在叶片前缘及叶身中部与叶根处，在叶尖部位积冰生成较少，甚至无积冰生成，在叶片尾缘中部附近生成较薄的积冰。试验积冰与模拟积冰位置基本一致，验证了本文模拟方法计算积冰位置的可靠性。

图2 积冰位置验证Fig.2 Verification of icing position: a) test icing position; b) simulated icing position

使用文献[26]的NACA0012 机翼的计算条件，对本文积冰厚度的准确性进行检验。计算条件：过冷水滴质量浓度为0.48 g/m3，液滴粒径为27.97 μm，来流速度74.97 m/s，具体验证条件见表2。本文计算方法得出的局部水收集系数和冰形与文献[26]中的局部水收集系数和冰形的对比如图3 所示。由图3 可知，采用本文计算方法得出的积冰极限位置与厚度结果与文献[26]的结果基本一致，从而验证了本文模拟计算方法的可靠性。

表2 验证计算条件Tab.2 Calculation conditions for validation

图3 收集系数和冰形验证Fig.3 Verification of collection efficiency and ice shape: a) collection efficiency; b) ice shape

2 结果和分析

2.1 研究对象与边界条件

本文的研究对象是叶片数为27 的跨音速风扇转子，该风扇转子进口轮毂比为0.6，外径为0.206 2 m，为增加计算效率，在模拟中使用加周期边界的单通道模型。

本文对飞行器处于高度为4 600 m 时，过冷水滴与冰晶混合相在风扇转子叶片积冰进行数值模拟，计算域的进口条件设置为总压进口边界，进口总压为65 kPa，进口总温为264.15 K，出口条件设置为平均静压出口，壁面为无滑移绝热边界。积冰采用间断最大结冰条件计算，该高度下水的总质量浓度为2.826 g/m3，液态水的质量浓度选择为0.5 g/m3，液滴与冰晶粒径均为20 μm，液滴与冰晶的初始温度为264.15 K。积冰计算条件见表3。

表3 积冰计算条件Tab.3 Icing computational condition

2.2 网格划分

本文使用商业软件Numeca 对该模型进行计算域网格划分，使用网格为六面体网格，网格如图4 所示。为保证计算结果的可靠性，对网格无关性的验证依据文献[12]，使用叶片壁面上收集系数β作为判断变量。分别采用了37 万、141 万、336 万和553 万等4 套网格进行计算，通过水收集系数与冰晶收集系数共同判断网格对计算结果的影响。经过计算，发现网格数大于141 万后，随着网格数增加，水收集系数与冰晶收集系数变化较小（如图5 所示），所以本文最终采用141 万网格进行计算。

图4 计算域网格Fig.4 Mesh used for calculation

图5 网格独立性验证Fig.5 Verification of mesh independency: a) water collection efficiency; b) ice crystal collection efficiency

2.3 数值模拟结果分析

2.3.1 转速对风扇叶片表面收集系数的影响

风扇叶片在18 000 r/min 时，不同截面的马赫数云图如图6 所示。由图6 可知，在90%叶高入口处存在“λ”激波，在50%叶高入口处存在斜激波，在10%叶高入口处为高亚音速流场。这说明该风扇叶片扭转较大，跨音速流影响过冷水滴与冰晶的速度，改变撞击角度，进而影响积冰。

图6 不同截面马赫数云图Fig.6 Mach number contours at different spans: a) 90% span from hub; b) mid-span; c) 10% span from hub

低转速与高转速下叶片表面的水收集系数如图7所示。低转速时，过冷水滴的高收集区域在叶根中部，随着转速的增加，高收集区域向叶尖尾缘方向发展；低收集区域主要存在于叶尖，随着转速的增加，叶尖的低收集区域减小，在叶根尾缘处开始出现低收集区并扩大。高转速时，过冷水滴的高收集区域为叶根中部与距轮毂3/5 尾缘线性连接区域，随着转速的增加，高收集区域无明显变化；低收集区域为叶尖中部与叶根尾缘，转速增加，低收集区域变化不明显。低转速时高收集区域的变化是由于随转速增加，气流曳力改变过冷水滴运动状态的能力减弱，过冷水滴撞击到叶片表面的数量增加导致的。叶根叶型弯角较大（如图6c 所示），增大了过冷水滴和冰晶的撞击角度，提高了风扇叶片对应位置的收集系数，这也同时导致了叶根尾缘过冷水滴与冰晶含量减少，出现低收集区域。高转速高收集区域无明显变化是由于过冷水滴在流道中速度较大，运动状态不易被气流改变，在惯性下运动造成的。叶尖叶型弯角较小（如图6a 所示），过冷水滴和冰晶的撞击角度减小，叶尖的收集系数减小。叶尖低收集区域低转速时，收集系数随转速增加而减小，高转速时无明显变化。这是因为过冷水滴低转速时容易受到气流影响，转速增加使撞击到叶尖的过冷水滴数量增加，高转速时气流对过冷水滴的影响较小，收集数量变化较小。

图7 叶片表面水收集系数Fig.7 Water collection efficiency on the blade surface: a) low velocity; b) high velocity

低转速与高转速下叶片表面的冰晶收集系数如图8 所示。冰晶低收集区域占据叶身较大区域，随着转速的增加，低收集区域减小。高转速时，冰晶高收集区域主要在叶根中部，低收集区域主要集中在叶尖。随着转速的增加，叶片表面的收集系数变化较小，变化主要集中在尾缘。低转速时，冰晶收集系数的变化主要是由于过冷水滴的收集变化引起的，过冷水滴收集增加，叶片表面未凝结的液态水增加，表面未融化的冰晶捕获增加；高转速时，冰晶收集系数变化较小是由于对应转速时叶片表面的液态水充足，撞击到水膜上的冰晶未捕获的极少，尾缘处的变化则是冰晶融化，表面形成水膜后撞击到叶片表面被捕获的。

与生理盐水组比较,各剂量组雌性大鼠的肝、肾、脑、胸腺、卵巢和子宫的脏/体比无显著性差异。各剂量组雌性大鼠的肾上腺和脾的脏/体比部分有显著性差异,但是综合考虑脏器的重量和大鼠的终期空腹体重,认为其差异无实际生物学意义。

图8 叶片表面冰晶收集系数Fig.8 Ice crystal collection efficiency on the blade surface: a) low velocity; b) high velocity

通过粒子撞击在叶片上的角度云图（图9）可知，低转速与高转速时，粒子在叶尖的撞击角度在5°以下，叶根的撞击角度最大，低转速时为19°，高转速时不小于30°。随着转速的增加，低转速时叶身尾缘的粒子撞击角度变大，高转速时粒子撞击叶片的角度基本不变。由此可见，叶型弯角影响粒子撞击角度，进而影响收集系数。由叶片压力面撞击角度的变化规律可知，相较于高转速，低转速时粒子易受到气流的影响而改变运动状态，影响收集系数。

图9 粒子撞击叶片角度Fig.9 Angle at which the particle hit the blade: a) low velocity; b) high velocity

考虑过冷水滴与冰晶的结冰机理，尤其是冰晶被捕获时其表面融化的液态水质量分数需达到一定值。在混合相积冰中，过冷水滴的捕获将在叶片表面生成水膜，有利于冰晶捕获。因此，为确定过冷水滴与冰晶的捕获原因，对叶片表面收集的过冷水滴温度与冰晶温度进行进一步分析。

低转速与高转速叶片表面收集到的过冷水滴温度如图10 所示。低转速时，撞击到叶片表面的过冷水滴温度低于冰点，且整个叶片表面过冷水滴温度无明显变化，过冷水滴撞击到叶片表面后易凝结形成积冰。高转速时，撞击到叶片表面的过冷水滴温度也低于冰点，但从前缘到尾缘的温度升高，且变化明显，这使得前缘未凝结的液态水向尾缘移动，并逐渐凝结，影响最终的积冰。

图10 叶片表面过冷水滴温度Fig.10 Temperature of supercooled droplets on the blade surface: a) low velocity; b) high velocity

低转速与高转速叶片表面收集到的冰晶的温度如图11 所示。低转速时，收集到的冰晶温度低于冰点，依据冰晶结冰机理，此时冰晶不易被捕获，但由于捕获的过冷水滴在叶片表面生成水膜，使表面未融化的冰晶可被捕获。高转速时，被捕获的冰晶从前缘至尾缘温度升高，并在尾缘处高于冰点，这使得除了表面未融化的冰晶被叶片表面过冷水滴形成的水膜捕获外，在不存在水膜的位置（如叶根尾缘）或水膜较少的区域，冰晶可以通过其自身表面融化形成的水膜被叶片捕获。

2.3.2 转速对截面冰形的影响

为进一步明确低转速与高转速积冰的异同及原因，依据不同转速叶片表面的积冰情况，选择距轮毂10%的截面进行具体分析。风扇叶片距轮毂10%截面处在不同转速下的积冰冰形如图12 所示。低转速时，叶片表面积冰厚度除前缘外，基本无明显变化，且积冰整体表面光滑。随着转速的增加，风扇叶片压力面及前缘积冰厚度增加。这是因为低转速时，过冷水滴与冰晶易受到气流的影响发生偏转，转速增加，气流偏转减小，撞击至风扇叶片表面的过冷水滴与冰晶增加，即对应位置处的过冷水滴与冰晶含量增加，增加了收集系数。另外，低转速时过冷水滴与冰晶的温度不易超过冰点，过冷水滴撞击后易凝结，水膜流动较少，故积冰厚度增加且表面光滑。高转速时，叶片表面积冰厚度变化较大，自前缘至中后部先增加、后减小，积冰表面粗糙，前缘生成明显的角状冰（如图13 所示），尾缘处几乎无积冰生成。随着转速的增加，前缘积冰厚度减小，积冰末端厚度增加。这是因为高转速时，转速增加，叶片表面未凝结的水膜受到的气流剪切力增加，水膜运动程度加剧，在前缘积冰量减少，前缘积冰厚度减小，水膜在流动中的冻结量减小，积冰极限位置向尾缘移动。

图12 不同转速冰形Fig.12 Ice shape at different velocities

图13 叶片前缘冰形Fig.13 Ice shape on the leading edge of the blade

由图14 可知，水膜的径向流动主要在前缘驻点及压力面叶身上半区，风扇转子前缘驻点附近水膜径向速度较大。在水膜流动过程中，该部分生成的积冰相较于其两侧较少，积冰在前缘位置生成角状冰。随着转速降低，驻点附近的水膜径向速度减小，积冰增多，前缘角状冰特征减弱。

通过水收集系数与冰晶收集系数对积冰冰形进行分析。风扇叶片在不同转速下的水收集系数和冰晶收集系数如图15 所示。前缘过冷水滴与冰晶收集系数差距较小，这是由于在进口处粒子的速度含量差别较小，流道中粒子收集系数差别较大是由于粒子在流道中的含量及速度差异较大。高转速时，叶身水收集系数在最大值点前后变化均比较明显。在最大值点前，水收集系数自前缘向尾缘方向增加；过最大值点后，水收集系数减少，至尾缘附近后水收集系数基本为0，最大值点后的变化比最大值点前的变化更加显著。低转速时，叶身水收集系数最大值点到前缘方向变化相对平缓，到尾缘方向变化相对明显。于冰晶而言，压力面的收集系数均比较平缓，只有高转速时前缘与尾缘有较为明显的波动。由撞击角度与冰晶温度可知，前缘收集系数较低是因为冰晶碰撞角度小，冰晶易弹跳而使收集系数减小，尾缘处冰晶表面融化，增加了收集系数。因此，水收集系数与冰晶收集系数较为完整地体现了积冰冰形的变化。

图15 不同转速下的水收集系数和冰晶收集系数Fig.15 Collection efficiency of water and ice crystal at different velocities: a) water collection efficiency;b) ice crystal collection efficiency

收集系数与对应部件位置的粒子含量相关，通过流道中的液态水含量与冰晶含量（如图16 所示）可知，过冷水滴与冰晶在流道中聚集。随着转速的增加，流道中过冷水滴与冰晶的聚集位置在低转速时向前缘移动，高转速时变化较小。此外，过冷水滴与冰晶在低转速时的聚集位置比在高转速时偏向出口，这与水收集系数和冰晶收集系数在叶身上的峰值和分布一致。

图16 流场中的液态水含量和冰晶含量Fig.16 Liquid water content and ice crystal content in the flow field: a) liquid water content; b) ice crystal content

过冷水滴与冰晶的聚集边界在低转速时呈现明显的弧形，高转速时过冷水滴边界近乎直线。通过过冷水滴与冰晶在流道中轴向速度（如图17 所示）可知，低转速时，过冷水滴与冰晶的轴向速度小，受到的惯性小，过冷水滴与冰晶随气流穿过流道。随着转速的增加，流道中过冷水滴与冰晶速度增加，运动状态不易被改变，故聚集位置向前缘移动。高转速时，过冷水滴与冰晶的轴向速度较高，过冷水滴与冰晶不易受到气流曳力的影响，因此聚集边界无明显变化。

图17 流场中过冷水滴与冰晶轴向速度Fig.17 Axial velocity of supercooled droplets and ice crystals in the flow field: a) supercooled droplets; b) ice crystals

流道中的冰晶温度如图18 所示。低转速时，流道中的冰晶表面温度未超过冰点，且无明显变化。高转速时，前缘附近冰晶表面温度低于冰点，自前缘向尾缘，冰晶表面温度升高，在尾缘附近，冰晶表面温度高于冰点，这是由于高转速时流道中的气流温度较高所造成的。此外，流道中的冰晶表面温度与叶片表面对应截面收集到的冰晶表面温度一致。

图18 流场中不同转速的冰晶温度Fig.18 Temperature of ice crystals at different velocities in the flow field

3 结论

通过对风扇叶片高转速与低转速积冰的研究，得出以下结论：

1）低转速撞击到叶片表面的粒子的运动方向与叶片的夹角小于高转速时形成的夹角，使得低转速时冰晶容易弹跳再次进入流道，叶表捕获量相对减少。

2）低转速时风扇叶片表面的温度均低于冰点，且整个叶身温度差较小，过冷水滴被捕获后未冻结的水膜极少，表面生成的积冰较为光滑。高转速时叶表温度的变化大，且存在高于冰点的位置，水膜流动较多，积冰表面粗糙。

3）低转速时流道内气流温度未升高至冰点以上，而高转速时气流温度在尾缘附近超过冰点，使得高转速时冰晶在不存在水膜的情况下可以被捕获，而低转速时的冰晶捕获只能依靠过冷水滴形成的水膜。