随着航空科技的发展，以飞机结冰为典型代表的相关问题，对飞行安全构成了严重威胁。因受到过冷水滴撞击，发动机的气动外形可能会出现结冰现象，使发动机进气效率受到影响，甚至可能造成飞行事故。涉及传热传质耦合复杂问题的飞机结冰过程及其防护措施，近年来备受关注。
国内外学者进行了大量的理论分析及方法研究，发展了基于拉格朗日方程和欧拉方程两种主要数值计算方法，并将其应用于相关计算分析软件系列程序中，在分析飞机机翼表面［1-3］及发动机表面［4］的水滴撞击特性方面取得了良好的进展。目前，用于计算固定表面水滴撞击特性的方法已经基本完善和成熟，其精度已经得到了许多实验的验证。
对于旋转部件，特别是轴对称旋转表面，例如旋转整流帽罩和航空发动机叶片，由于受切向速度引入的离心力影响，水滴撞击与固定表面存在明显差异。关于旋转部件水滴撞击的研究相对较少。Reid等［5］、Switchenko等［6］较早地关注了旋转部件的结冰问题，利用FENSAP-ICE软件分析了轴对称旋转表面和其他复杂旋转表面的水滴撞击特征。赵秋月等［7］采用拉格朗日法对航空发动机旋转整流帽罩的水滴撞击特性进行了分析。王治国等［8］利用CFX模型和粒子输运模型对航空发动机旋转面水滴撞击进行了计算。易贤等［9］采用多参考坐标系计算了空气流场，并利用欧拉法对风力机三维旋转叶片的集水效率进行了计算。吴孟龙等［10］基于欧拉法模拟了水滴在旋转整流帽罩上的撞击。王健等［11］对旋转整流罩的积冰现象进行了实验研究。
本文针对飞机结冰问题的部分重要环节，即整流帽罩的水滴撞击特性，开展了模拟仿真计算。通过引入惯性力将非稳态周期性旋转边界条件简化为稳态边界条件。用欧拉法将空气流场和水滴流场描述为单向耦合形式。在非惯性参考系中引入科氏加速度以修正欧拉控制方程。使用有限体积求解器求解空气、水滴流动的控制方程，将惯性力定义为动量方程的源项，获得空气、水滴流动的速度、体积分数及水滴收集系数。采用上述方法，模拟计算并分析了航空发动机整流帽罩在静止和不同转速情况下的水滴撞击特性。
1 旋转流场的数学模型
惯性坐标系下通常建立固定部件的绕流流动模型。将非旋转条件扩展为旋转条件的方法有4种，分别是单旋转坐标模型、多参考坐标系模型、混合平面模型及滑动网格模型。其中，单旋转坐标模型是假定坐标固定于旋转元件上，坐标旋转时的切向速度与旋转元件的切向速度一致。因此，本文将旋转边界转换为静止参考坐标，将非定常流动转换为定常流动，从而简化模型的处理和计算过程。
1.1 空气流场模型
首先，考虑部件匀角速度定轴旋转，建立单旋转参考坐标系。当牵连运动为定轴运动时，绝对加速度等于牵连加速度、相对加速度和科氏加速度的总和。因此，在单旋转坐标系下的空气流场控制方程为［11］
式中： 为空气密度； 为压力； 为以相对速度为基础的黏性应力； 为外力矢量； 为转速； 为控制体积的失径；2 v ar为科氏加速度； （ ）为牵连加速度。 可由式（3） 计算得到。
式中： 为空气的绝对速度； 为牵连速度，可由式（4） 计算得到。
1.2 水滴撞击模型
在结冰的天气条件下，过冷水滴随着气流一起运动，可能会撞击到飞机的迎风表面。由于水滴直径小、空气中水滴的含量少，因此忽略水滴对气流的影响，仅考虑空气流动对水滴运动的影响（单向耦合）。对于复杂表面，欧拉方法具有明显的优势，它将水滴视为连续相，通过引入水滴体积分数来表示控制体内的水滴体积占总体积的大小，从而建立控制方程。水滴的欧拉方程在单旋转坐标系中表示为［12］
▽ρα v r=0
▽（ρα v r v r）+ρ［2 v r + ］
=ραK（ v ar- v r）
式中： 为水滴体积分数；ρ为水的密度； v r为水滴的相对速度，通过式（7） 计算。
v r = v - u r
K=
式中： 为空气的动力黏度； 为水滴直径； 为阻力函数，通过式（9） 计算。
f =
式中： 为球体的阻力系数，可由式（10） 计算得到。
C D=
为相对雷诺数，定义为
1.3 水滴收集系数模型
求解了水滴的欧拉方程之后，可直接得到旋转表面的水滴撞击区域及水滴撞击量，进一步计算出积冰和除冰热。将实际水滴收集值与局部表面可能的最大值之比定义为水滴的收集系数，使计算更为简单化。根据定义，旋转坐标下的水滴收集系数 可以描述为
式中： 为局部法向量； 为远场水滴速度； 为水滴在远场中的体积分数。
由式（7） 和式（12） 可知，旋转表面的水滴收集系数既与远场空气流速有关，也受旋转速度的影响。当转速足够大时， 可能大于1，说明旋转表面的收集系数可能大于固定表面的最大收集系数。
2 计算模型
2.1 几何和网格
本文采用ICEM软件对三维模型进行网格划分。由于六面体结构化网格在数值模拟中具有较强的优势，因此本文采用六面体结构化网格对流场进行拓扑划分。为了更准确地计算表面水滴撞击特性，对整流帽罩的表面网格进行局部加密。在下游和远场位置，适当增加网格尺寸以加快计算速度，网格细节如图2所示。
2.2 仿真模型计算方法
由于空气运动和水滴流动是单向耦合的，因此分两个单独步骤并利用Fluent软件完成计算。
发动机整流帽罩沿中心轴线以恒定速度旋转，建立以中心轴为基准的单旋转坐标系。采用三维可压缩流体求解器和 湍流模型来迭代求解方程（1）和方程（2），得到旋转部件周围的空气流场，由于整流帽罩与坐标系同步旋转，整流帽罩相对速度为零。实际的固定远场则表现为相反的转速。
2.3 计算条件
首先，计算在非旋转条件下不同来流速度的影响；其次，引入旋转边界条件进行对比。计算条件如表1所示。
表1 计算条件 |
| 参数 | 条件1 | 条件2 | 条件3 | 条件4 | 条件5 |
|---|---|---|---|---|---|
| Mach | 0.15 | 0.30 | 0.60 | 0.15 | 0.15 |
| RPM | 0 | 0 | 0 | 3 120 | 5 200 |
所采用的欧拉模型只能求解单直径的水滴。但在真实的大气中，水滴的直径通常具有某种分布类型。此外，单个水滴直径不能反映由旋转引入的惯性影响。因此，本文采用L-D分布，其平均容积直径（mean volume cliameter，MVD）为11.5 ，由7种不同直径的水滴组成。表2给出了直径分布和相应的体积分数。计算完7个直径的水滴收集系数后，最终的总收集系数可由式（13） 计算。
式中： 和 分别为各直径的体积分数和收集系数。
表2 直径分布和相应的体积分数 |
| di /μm | 3.6 | 6.0 | 8.2 | 11.5 | 15.8 | 20.0 | 25.5 |
|---|---|---|---|---|---|---|---|
| p i | 0.05 | 0.10 | 0.20 | 0.30 | 0.20 | 0.10 | 0.05 |
3 仿真结果与分析
3.1 非旋转情况下来流速度的影响
水滴直径的L-D分布也与单一直径的计算结果存在很大的差异。从图4可以看出，L-D分布会导致高撞击区域特别是驻点处的收集系数值变小，但会扩大撞击区域并增加下游位置的收集系数值。因此，由直径分布所导致的水滴撞击特性的改点是显著的，在分析结冰或设计除冰系统时必须计算这些特性。
3.2 旋转的影响
本文计算了不同转速下空气流速大小和水滴体积分数，如图5所示。
旋转通过壁面剪切应力对空气流场的边界层产生影响，因此，不同的转速有着不同的壁面绝对速度。由于整流帽罩形状的特殊性，其边界层较薄，所以在不同的转速下，大部分空气流场都是相同的。由图5还可以看出，不同条件下水滴体积分数云图基本相同。这是因为旋转只在水滴到达边界层并非接近表面时才会影响水滴的运动。由于水滴的惯性，水滴受到剪切应力的影响较小，同时也没有足够的时间来改变运动。
不同转速下的水滴最大收集系数和碰撞极限如表3所示。最大收集系数为驻点对应的值，碰撞极限是沿曲线的长度。结果表明，随着转速的增大，最大收集系数降低，碰撞极限增大。但是差值非常小，可以忽略不计。因此，在大多数情况下，旋转对水滴撞击特性的影响可以忽略不计。
表3 不同转速下的水滴最大收集系数和碰撞极限 |
| 转速/（rad·m·s-1） | 最大收集系数 | 碰撞极限 |
|---|---|---|
| 0 | 0.313 591 08 | 0.069 422 00 |
| 326.73 | 0.313 569 31 | 0.069 837 80 |
| 544.54 | 0.313 533 21 | 0.070 212 83 |
4 结论
本文利用单一旋转坐标系和欧拉方法，建立了水滴在旋转表面撞击特性的三维撞击模型。并利用该方法对航空发动机旋转整流帽罩进行了数值计算和分析。结果表明，随着马赫数的增加，水滴收集系数增加，撞击区域增大。水滴直径分布对撞击特性也有明显的影响，其结果比单一直径更准确全面。而整流帽罩的旋转对水滴撞击特性的影响则可以忽略不计。