无人机螺旋桨三维积冰计算与积冰控制检测

栗 枢，曹广州

(1.南京航空航天大学 能源与动力学院，江苏 南京 210016；2.南京航空航天大学 无人机研究院，江苏 南京 210016)

当飞机穿过含有过冷水滴的云层时，空气中的过冷水滴会撞到发动机进口、机翼无人机螺旋桨等迎风部件表面，并产生结冰现象。这些迎风部件表面形成的结冰现象会破坏表面气动外形，危及飞机安全，甚至发生坠毁。因此研究飞机积冰问题，在工程上具有十分重要的意义。

能量有限、结构受限、对重量敏感等特点导致中小型无人机的防除冰设计相对困难，效费比太低，因此，在有人机中成熟使用的防除冰系统一般不能直接用于无人机上。目前来看，用户们普遍采取的应对方法是在环境条件较好的时间段进行无人机的飞行，比如：科索沃战争中美军为避免结冰危害，不得不在10月至4月期间禁飞无人机。这样就会大大影响无人机的出勤率。即便如此，仍有25%的架次遭遇结冰，导致不同程度的故障或坠毁。所以说，结冰是不可避免的，即使是环境条件相对好的情况下会产生结冰现象。

国外对有人机迎风表面积冰的数值模拟在20世纪40年代就已开展，1953年Messinger首次提出用于积冰数值计算的模型，此后该模型得到了广泛的应用。20世纪70年代，随着CFD的发展，积冰数值模拟技术快速发展，80年代之后陆续出现了一批基于Messinger模型开发的积冰模拟/防冰设计软件，如美国的LEWICE[1]、英国的 TRAJICE[2]、法国的 ONERA[3]、意大利的 CIRA[4]以及加拿大的 FENSAP-ICE[5]、Myers[6-7]等。但Messinger模型只是对霜冰的计算分析较为准确，认为当前控制体中未凝结的液态水全部流入下一个控制体中，没有考虑冰层表面薄水膜流动对积冰的影响[2]。

国内绝大多数的积冰模型都是以Messinger模型为基础开发的二维积冰模型，而曹广州[8-9]则在处理水膜流动的问题上通过量纲分析的方法简化了N-S方程，得到了3个水膜流动的方程，建立了国内第一个考虑水膜流动的三维积冰模型并开发其计算方法[10]。该模型认为水膜在重力和空气剪切力的共同作用下会分别向展向和弦向流动，并对国外典型的三维积冰进行了数值模拟研究，结果匹配良好。

对于旋转部件，陈宁立[11-12]在曹广州普通迎风部件表面三维积冰与薄水膜流动耦合的基础上，考虑了离心力对水膜流动的影响，发展了旋转部件表面薄水膜流动与三维积冰相变耦合的数学模型，并对直升机螺旋桨、发动机进口整流罩等典型旋翼结构进行了数值模拟研究，又通过与实验结果对比验证了考虑离心力作用下该三维积冰模型的准确性。梁鹏[13-14]自主搭建小冰风洞，通过可视化方法得到不同条件下的积冰形状，并通过图像分析了不同条件下的结冰规律。

本文考虑到螺旋桨翼型的三维扭转特性，开发了适合螺旋桨三维积冰的结构化-非结构化的混合网格划分方法。通过CFX流场计算，找到几组最佳速度-转速匹配，在考虑到飞行来流和旋转来流情况下使三维积冰模型更贴合螺旋桨气动外形。又考虑到旋转条件下水膜的甩脱，提出考虑飞行来流的旋转条件下的撞击特性和三维积冰模型的数学模型及其计算方法，利用Fortran开发了相应的积冰模拟程序，对整个积冰过程进行检测，并进行了不同速度/转速、温度、水滴含量和水滴直径下的积冰数值模拟和对比分析。

1 螺旋桨表面三维积冰模型与检测方法介绍

1.1 基本假设

对于旋转部件，将水膜受到离心力的作用而产生脱离水膜主体的现象称为“甩脱”。针对“甩脱”现象，做出以下假设。

① 旋转部件表面的薄水膜流动可以看成不可压缩的层流流动，且可以假设在小的时间微元内水膜的流动为稳态流动。

② 普通积冰部件表面未凝结的水会形成连续的薄水膜，该水膜厚度的量级比较小[15]，一般在10-4m量级，旋转部件考虑到离心力的作用，并考虑到转速的不固定性，水膜厚度的量级一般在10-7，转速相对较高的无人机螺旋桨可达10-8。

③ 暂认为水膜主体的能量、动量的不受水膜甩脱的影响。

④ 暂认为周围的气液稀疏两相流场不受水膜甩脱的影响。

⑤ 暂不考虑水滴撞击所带来的溅洒。

⑥ 旋转积冰部件表面的薄水膜的温度一般很低，因此暂不考虑辐射换热所带来的影响。

1.2 螺旋桨表面水膜流动与三维积冰模型

积冰的形成首先是水滴撞击到飞机的迎风部件表面形成水膜，之后在空气-过冷水滴两相流场中对流换热逐渐形成冰层。所采用的积冰控制体如图1所示，该控制体由3个部分组成，最上方的是空气-过冷水滴两相流场，中间是未冻结水膜，最下方则是冰层。其中，Hw、Tw分别为水膜厚度和温度；Hi、Ti分别为冰层厚度和温度；mimp为水滴撞击的质量；mice为积冰的质量；min为流入该控制体的质量；mout为流向下一个控制体的质量。

图1 积冰控制体

建立积冰与水膜流动相耦合的数学模型，具体如下：

(1)

(2)

(3)

(4)

(5)

式中，ρi、ρw分别为冰和水的密度；rx、ry、rz分别为r在x、y、z方向上的分量，r为该点到旋转轴的距离向量的反向;Lf为积冰的相变潜热，本文取334400 J/kg；λi、λw分别为冰和水的导热系数。

其中，式(1)～式(3)分别为冰层表面水膜流动的连续方程、动量方程和能量方程，式(4)为冰层中的能量方程，式(5)为水膜的积冰相变方程。

1.3 局部水收集系数的选取

局部水收集系数β是积冰模拟过程中用于表征水滴撞击特性的无量纲参数，反映出结冰的范围与程度。局部水收集系数是指当地微元面上的实际水滴撞击量与相同面积上按自由来流条件获得的水滴撞击量之比。采用欧拉-欧拉法计算两相流场时，可用式(6)和式(7)来计算局部水收集系数。

(6)

(7)

式中，αw为壁面上水滴的体积含量；Unw为壁面上水滴的法向撞击速度，取壁面外法线为正方向，因此Unw均为负值，式(6)中的负号是为了保证局部水收集系数β的值均为正；U0为螺旋桨表面受到合来流速度，数值上等于飞行来流U∞和旋转来流Ur的合速度，相比于普通翼型的水收集系数计算多了一项旋转来流。

1.4 积冰表面水膜甩脱的依据

当气流曳力做的功和水膜由于受到离心力增加的潜能大于水膜和螺旋桨表面的黏附功时，水膜会发生甩脱。具体得出：

We>2

(8)

得到临界水膜厚度为

(9)

式中，We为韦伯数，是衡量冰层表面水膜动力学特性的无量纲参数；ρa为气流密度；Va为当地气流速度；σiv为表面张力系数。

因此，对螺旋桨水膜甩脱的判据均是基于韦伯数We来判断的，认为当韦伯数We>2时水膜就会发生甩脱。

1.5 计算方法

本文的计算方法步骤如下：

① 积冰从未积冰的光滑壁面开始，t=0时，Hw=0，Hi=0，开始进行三维建模、网格划分、两相流场计算，并通过自开发软件提取流场数据。

② 计算得到积冰厚度、水膜厚度等参数。

③ 计算We>2是否成立，若满足，计算求解新的水膜厚度。

④ 时间步推进，继续下一时间步的求解。

2 螺旋桨表面三维积冰过程的检测内容与积冰控制

2.1 模型介绍

选取2714型螺旋桨进行三维积冰模拟，如图2所示，桨的直径27 in，约为0.6858 m(1 in=0.0254 m)。由于前方机身的遮挡，螺旋桨转轴附近受流场影响较弱，故做简化处理，去除螺旋桨转轴附近的桨叶部分，由于桨尖效率低且结构复杂，计算模型半径缩减至为0.31 m。

图2 扭转性极强的2714型螺旋桨计算模型

由于中小型无人机一般采用两叶螺旋桨，故只选取一半进行研究，所建立的计算域为图3所示的半圆柱结构，图3(a)中红色区域为整个计算域，蓝色区域为模拟的机身部分，黄色区域为内部结构化网格部分，包含螺旋桨，同时考虑前面机身遮挡的影响而挖去一部分。螺旋桨的弦长在各个截面处都不相同，范围为0.03～0.05 m之间，内部计算域为提高计算精度，螺旋桨所受到的旋转来流的上游和下游取单倍弦长，厚度取桨长。而在外部非结构化区域，要考虑部分机身的长度，飞行来流的上游和下游为2∶1的15倍弦长，圆柱的半径为1.5倍的桨长。半圆柱的两个侧面采用周期性边界条件，图3(b)为内部结构化网格的放大图。利用ANSYS-ICEM对螺旋桨进行网格划分，为了减少网格量，同时确保计算精度，采用了结构-非结构的混合网格，如图3所示，对螺旋桨的近壁区采用质量较好、密度较大的结构化网格。而在外部区域采用较为粗糙的非结构化网格。采用表1的参数进行流场设置和数据提取。

图3 2714型螺旋桨的计算域及网格

表1 2714型螺旋桨的计算条件

图4为选取不同半径位置螺旋桨的截面图，以螺旋桨各个截面上最左端为驻点位置，上表面即压力面的弧长坐标值为正值，下表面即吸力面坐标值为负值。由图4可见，螺旋桨在不同的半径位置的截面在形状、角度甚至几何尺度上都各不相同，体现了其强扭转特性。由1.3节局部水收集系数公式可知，其本身与气流速度有关，这正与要研究的流场速度-转速匹配相关联，其次本文最终旨在研究积冰的影响，局部水收集系数作为积冰的重要参数，影响着整个积冰过程和结果，故选用该参数作为一些算例展示及验证分析。

螺旋桨表面的局部水收集系数分布曲线如图4所示，其中横坐标为螺旋桨的弧长坐标。由图4可知，压力面即S>0方向的水滴撞击极限要远远大于吸力面即S<0方向上的水滴撞击极限，因此局部水收集系数呈不对称分布。且在S>0方向，旋转半径越大的位置，水滴撞击极限越广。这是由于随着旋转半径的增大，使得线速度增大，来流水滴撞击速度增大，水滴惯性也随之增大，水滴就更容易撞击到螺旋桨表面，其局部水收集系数也相应增大。

图4 光滑螺旋桨不同半径位置截面形状图

2.2 网格的划分和网格独立性验证

在算例1的条件下，通过网格单元数分别为65万、85万、100万、120万、140万进行流场计算，并通过数据提取得到相应的局部水收集系数数据，如图5所示，工程上螺旋桨总半径的0.7倍的位置一般为螺旋桨的最初设计点。故取螺旋桨表面在r/R=0.7截面上的局部水收集系数，由图5可知，当网格单元数量达到100万时，局部水收集系数基本上不再变化，可以判断在所采用的网格疏密范围内，计算结果受网格量的影响可以忽略不计。最终采用网格单元数为100万的网格进行计算分析。采取网格单元的弦长随着展向位置的不同而不断变化，故本文在不加说明的情况下，均取螺旋桨总长度的0.7倍位置处的积冰参数进行研究。

图5 网格独立性验证(r/R=0.7)

2.3 不同温度对螺旋桨积冰影响的检测

选取温度为-5 ℃、-10 ℃、-15 ℃、-20 ℃，研究不同温度对螺旋桨积冰的影响。图6为T=-10 ℃下r/R=0.7截面处的局部水收集系数分布。由局部水收集系数公式可知，局部水收集系数仅与水滴密度、体积分数、水滴含量、速度有关，与环境温度的高低无关，而且温度的改变不会影响水滴撞击量的改变，所以4种不同环境温度下的局部水收集系数分布相同。

图6 T=-10 ℃下局部水收集系数分布

不同环境温度下r/R=0.7截面位置冰层厚度分布情况如图7所示。由图7可见，随着环境温度的降低，相同截面上的冰形高度越大，冰形的覆盖范围也越大。

图7 不同环境温度下的冰层厚度分布

不同环境温度下r/R=0.7截面处的最终冰型轮廓图如图8所示。由图8可知，随着温度的降低，结冰位置基本上不受影响，因为未能从气动上改变整个流场，只是温度更低了，水滴更易凝结成冰，在相同的位置上，积冰量更大，冰型更厚，逐渐形成更大的角状冰。当水滴撞击量相同时，温度越低，温度的变化量增大，水膜表面与两相流场间的换热量增大，相同质量流量下水滴吸收的热量增大，根据热力学平衡定律，螺旋桨表面的积冰量就变大。同时环境温度越低，空气能吸收的热量也越多，由热力学平衡定律可知，相应的积冰放热量增大，螺旋桨表面的积冰量增大。

2.4 不同水滴含量对螺旋桨积冰影响的检测

选取水滴含量分别为0.5 g/m3、1.0 g/m3、1.5 g/m3、2.0 g/m3,研究不同水滴含量对螺旋桨积冰的影响。不同水滴含量下r/R=0.7截面位置冰层厚度和最终冰型轮廓图分布情况如图9、图10所示。由图9可见，随着水滴含量的增加，结冰位置逐渐向下游发展，这是因为水滴含量增加了，更多的水膜流向下游，边流动边冻结成冰，即所谓的明冰。所以在上游的冰型几乎不变，而在下游处，水滴含量较大的条件使冰型凸起，积冰量变大，逐渐形成更大的角状冰。

图9 不同水滴含量的冰层厚度分布

图10 不同水滴含量的冰型轮廓图

2.5 不同水滴直径对螺旋桨积冰影响的检测

选取水滴直径分别为15 μm、20 μm、25 μm、35 μm，研究不同水滴直径对螺旋桨表面积冰的影响，不同水滴直径下r/R=0.7截面的局部水收集系数分布情况如图11所示。由图11可知，随着水滴直径的增大，局部水收集系数有所增加。因为水滴直径越大，水滴的惯性就越大，从而使水滴更容易撞击到螺旋桨表面。

图11 不同水滴直径的局部水收集系数分布

不同水滴直径下r/R=0.7截面的冰层厚度分布情况如图12所示。由图12可知，随着水滴直径的增加，积冰的范围和厚度均有所扩大和增加。

图12 不同水滴直径的冰层厚度分布

不同水滴直径下r/R=0.7截面的最终冰型轮廓图如图13所示。由图13可知，上游的冰型轮廓基本保持一致，在下游处，随着水滴直径的增加，冰型有所扩大。

图13 不同水滴直径的冰型轮廓图

3 结论

本文将适用于旋转部件表面积冰模拟的模型应用到扭转性极强的2714型螺旋桨中，介绍了基于2714型螺旋桨强扭转特性下的模型，并进行了积冰的检测，通过研究环境温度分别为-5 ℃、-10 ℃、-15 ℃、-20 ℃，水滴含量分别为0.5 g/m3、1.0 g/m3、1.5 g/m3、2.0 g/m3以及水滴直径分别为15 μm、20 μm、25 μm、35 μm对螺旋桨表面积冰规律的影响，分析得出以下结论：

① 随着环境温度的降低，局部水收集系数和对流换热系数基本不受影响，检测出空气中的水滴更容易结冰，冰层厚度增大，水膜厚度增大，冰型轮廓逐渐扩大。

② 随着水滴含量的增大，检测出冰层厚度增大，水膜厚度增大，上游的冰型轮廓几乎不变，更多的水滴流向下游，从而形成明显的角状冰。

③ 随着水滴直径的增大，局部水收集系数增大，对流换热系数几乎不受影响，检测出积冰量增大，水膜厚度略有增加，上游积冰形态基本上不变，下游的积冰量越来越大，逐渐形成明显的角状冰。