0 引 言

随着航空科学与技术的发展，以飞机结冰气象条件为典型代表的气象因素成为威胁飞行安全的重要因素。过冷水滴撞击导致的飞机结冰现象可能破坏升力部件气动外形，影响发动机进气效率，甚至引起飞行事故。飞机结冰过程及其防护方法研究涉及复杂的传热传质耦合问题，近年来得到了极大的重视，国家重点基础研究发展计划、国家安全重大基础研究计划、国家自然基金等都将其纳入了研究计划。本文针对飞机结冰的重要环节，即水滴撞击特性分析，开展计算方法研究。

水滴撞击特性决定了结冰的范围与质量，其计算是飞机结冰与防除冰仿真分析的基础。国内外开展了大量的理论分析与计算方法研究，发展了拉格朗日法与欧拉法两种主要的水滴撞击计算方法，开发了一系列的计算分析软件^[1-4]，并将其应用于各种复杂的飞机部件表面，取得了较好的效果。目前，固定部件表面的水滴撞击计算方法已经基本成熟，其准确性得到了学者们的共识。

对于旋转部件，特别是具有非中心对称特点的旋转表面(如螺旋桨、发动机旋转叶片等)而言，由于部件旋转引入的离心力等惯性力影响，其水滴撞击过程与固定表面存在显著差异，针对旋转部件的水滴撞击研究分析较少。加拿大的FENSAP-ICE软件^[5-6]对旋转部件结冰问题的关注较早，分析了中心旋转体及复杂旋转表面的水滴撞击特性。国内，赵秋月、董威等^[7]用拉格朗日法对发动机旋转整流帽罩的水滴撞击特性进行了分析。王治国等^[8]应用CFX软件及其粒子输运模型计算了发动机旋转表面的水滴撞击特性。易贤等^[9]采用多参考坐标系方法下的空气流场结合欧拉法下的水滴计算，得到了风力机三维旋转表面水滴收集率。吴孟龙、常士楠等^[10]基于欧拉法模拟计算了旋转帽罩的水滴撞击特性。王健、胡娅萍等^[11]对旋转整流罩的结冰生长开展了试验研究。

本文采用单旋转坐标系模型来研究复杂表面稳定旋转时的空气流动与水滴运动，将其抽象为旋转坐标系下的定常流动，用欧拉法来计算水滴撞击特性，并分析水滴撞击计算方法的有效性及旋转速度对计算结果的影响规律。

1 旋转流数学模型

固定部件绕流流动模型通常是建立在惯性参考坐标系中，将其扩展到旋转体流动通常有四种处理方法：单旋转坐标系模型、多参考系坐标模型、混合平面模型及滑移网格模型。其中单旋转坐标系模型将坐标固定在旋转部件上，坐标随部件以相同的速度转动，因此旋转边界相对于坐标系是静止的，非定常的旋转流动转化为定常流动，从而简化计算过程。本文采用单旋转坐标系模型，建立旋转坐标系下的空气流动与水滴运动模型，发展旋转部件复杂表面水滴撞击计算方法。

1.1 空气流动模型

考虑部件匀角速度定轴旋转，建立单旋转参考坐标系。当牵连运动为定轴转动时，动点的绝对加速度等于它的牵连加速度、相对加速度和科里奥利加速度的矢量和。因此，单旋转坐标系下定常空气流动的控制方法变为^[12]：

(1)

(2)

式中，ρ_a为空气密度；p为压力； τ _r为基于相对速度的粘性应力； F 为外力矢量； ω 为旋转速度； r 为控制体矢径； v _ar为空气相对速度，由速度合成定理可得：

(3)

其中，v _a为空气绝对速度； u _r为牵连速度，由下式计算：

(4)

此外，式(2)中的2 ω×v _ar为科里奥利加速度； ω×(ω×r )为牵连加速度，即向心加速度。

1.2 水滴运动模型

在结冰气象条件下，过冷水滴随气流一起运动，迎风部件表面附近的水滴可能撞击到飞机壁面。由于水滴尺寸较小，空气中的液态水含量少，只考虑空气流动对水滴运动的影响，忽略水滴对气流的影响(单向耦合)。对于复杂表面的水滴撞击计算，欧拉法有显著的优势。欧拉法将水滴视为连续相，引入水滴容积分数来表示控制体内的水滴体积占总体积的大小，建立控制方程。用与空气流场同样的处理方法将水滴欧拉方程扩展到单旋转坐标系下，得：

(5)

(6)

式中，α是水滴容积分数；ρ为水滴的密度； v _r为水滴相对速度，由速度合成定理可知：

(7)

其中，v 为水滴绝对速度。

式(6)中的K为空气-水滴交换系数，由下式定义^[13]：

(8)

其中，μ_a表示空气的动力粘度；d_p为水滴直径；f是阻力函数，可由下式计算：

(9)

其中，阻力系数C_D由下式计算^[13]：

(10)

Re为相对雷诺数，由下式给出：

(11)

1.3 水滴收集系数计算模型

用欧拉法求解水滴控制方程后，可直接得到旋转部件表面的水滴撞击区域和各处的水滴撞击量，用于飞机结冰及防除冰计算分析。为便于分析比较，计算撞击表面水滴收集系数，其定义为表面上的实际水收集量与该表面上最大可能的水收集量之比。根据定义可推出旋转坐标系下的局部水滴收集系数β的表达式为^[9]：

(12)

式中， n为局部壁面的单位法向矢量；V_∞为远场的水滴速度；α_∞为远场的水滴容积分数。

由式(12)与式(7)分析可知：旋转运动时的表面水滴收集系数不仅与远场空气速度有关系，还受旋转速度的影响；当旋转速度较大时，β值可能大于1^[9]，即旋转表面的水收集量比表面固定时最大可能的水收集量还要大。

2 计算条件与模型实现方法 2.1 几何模型与网格

以发动机进气道内部的旋转帽罩与叶片为研究对象，建立简化的几何模型，如图 1所示。图中的外部圆柱轮廓表示发动机腔的内壁面，中间位置用60°的三角圆锥来表示发动机整流帽罩，用三个互为120°的椭圆柱来表示旋转叶片。为保证计算的收敛性，将发动机腔壁向前后分别进行拉长处理，前后分别为速度进口与压力出口。

将几何模型导入ICEM软件中进行网格划分，选用非结构方法，并在内外壁面进行边界层网格的设置。由于水滴撞击的区域一般在迎风表面，在整流帽罩与椭圆柱表面进行网格加密处理，最终的网格结果如图 2所示。

2.2 旋转模型实现方法

以通用流体计算软件FLUENT为平台，通过二次开发来实现旋转部件复杂表面的水滴撞击特性计算。由于空气流动与水滴运动单向耦合，将空气与水滴进行分步计算。

发动机整流帽罩与叶片绕中心轴以恒定速度转动，以转轴为基础建立单旋转坐标系。采用三维粘性不可压流求解器，选用k-ε湍流模型来迭代求解旋转坐标系下的守恒方程式(1)、式(2)，得到旋转部件周围的空气绕流流场。需要注意的是，帽罩与叶片随旋转坐标系以相同的速度转动，其相对速度为0；而原本固定的发动机腔壁，在旋转坐标系下则是以相同的转速反向转动。

在空气流场的基础上，用欧拉法进行单旋转坐标系下的水滴流动。参考惯性系下的水滴求解方法，应用FLUENT软件提供的UDS(User Defined Scalar，用户自定义标量)输运方程来进行旋转水滴计算^[13-14]。分析旋转下的欧拉方程，将水滴容积分数、水滴相对速度沿坐标轴的三个分量分别设为自定义标量。同时将旋转坐标下产生的牵连项与科里奥利项放入水滴动量方程的源项中，建立UDS框架下的水滴输运方程。通过有限容积法对输运方程求解，便可以得到旋转作用下的水滴场分布。需要注意的是，进行水滴撞击壁面后的吸水处理时，考虑旋转作用的影响，用水滴的相对速度来进行。

2.3 计算状态

首先用基本成熟的固定部件水滴撞击计算方法对几何模型进行静止时的计算。之后设一个非常小的旋转速度，应用建立的单旋转坐标系模型进行求解，与传统的固定坐标方法进行比较，验证有效性。然后，取其他两个不同旋转速度进行计算，分析旋转速度对水滴撞击特性的影响。

模型计算时，速度进口条件取30m/s，出口压力取标准大气压，水滴直径取20μm，旋转计算时转速分别取0.01rad/s、30rad/s、60rad/s。

3 计算结果与分析 3.1 旋转模型验证

图 3所示为当叶片静止时，用传统固定坐标系模型^{[13, 15]}计算得到的局部水滴收集系数分布结果。从图 3中可以看出，水滴主要撞击区域为帽罩的前缘点以及叶片的前缘驻点附近。

图 4所示为旋转速度为0.01rad/s时，用单旋转坐标系模型计算得到的帽罩与叶片表面速度分布结果。从图 4中可以看到，旋转使壁面产生运动，离旋转轴越远速度越大。

图 5为通过旋转坐标下水滴撞击模型计算得到的转速为0.01rad/s时的水滴收集系数分布结果。从图 5中可以看到分布结果与静止时传统方法的结果类似。为更好地与传统固定参考系模型进行比较，取叶片的一个截面来进行显示，截面的位置为图 4中上边叶片较亮线处。

图 6为固定参考系模型计算得到的截面结果与旋转模型计算转速为0.01rad/s时水滴收集系数结果的比较。从图 6中可以看到水滴撞击范围与收集系数分布都非常一致，在一定程度上说明了所建立的单旋转坐标水滴撞击模型是有效的。

3.2 旋转速度的影响

图 7与图 8分别为旋转速度为30rad/s与60rad/s时(旋转速度方向为图中+x方向)，表面水滴收集系数分布云图。从图中可以看到旋转作用下的最大局部水滴收集系数比静止时的结果大，且当转速为60rad/s，水收集系数存在大于1的区域。图中叶片表面水滴撞击区域向右侧表面偏移，体现了旋转速度对水滴运动的影响，而帽罩由于中心对称，水滴收集系数变化较小。

着重分析旋转速度对叶片的影响，取上文所述截面进行比较显示，结果如图 9所示。从图 9中可以看到，截面上的局部水收集系数向-y轴方向移动，且旋转速度越大，区域移动越大。这是由于旋转作用使叶片表面的相对速度发生变化，从而改变驻点位置，且转速越快，合成速度产生的驻点变化越大。

4 结 论

采用单旋转坐标系建立了旋转表面水滴撞击模型，并对发动机旋转帽罩与叶片进行了计算分析。结果表明：所建立的方法能有效地计算旋转部件复杂表面的水滴撞击特性；旋转作用使水滴能撞击区域向相对迎风的叶片侧表面移动，且旋转速度越大，水滴撞击区域移动越大。所建立模型的计算准确性还需要经过试验来验证。