直升机旋翼结冰后三维桨叶气动特性数值分析

王正之,杜孟华,朱春玲,刘森云

(1.南京工程学院能源与动力工程学院，南京，211167；2.南京航空航天大学航空学院，南京，210016；3.中国空气动力研究与发展中心结冰与防除冰重点实验室，四川绵阳，621000)

直升机相较于固定翼飞机而言飞行高度低，前飞的飞行速度慢，非常容易遭遇结冰气象环境。直升机旋翼桨叶结冰后，改变了桨叶翼型剖面的形状，旋翼的空气动力学特性被破坏，使得升力降低，阻力增加。同时，结冰后，旋翼在高速旋转，每个截面提供的升力被改变，使升力分布不再对称，加剧了直升机的振动，危及直升机的安全。报道表明我国现役的直升机中，诸如米-8、超黄蜂等在冬季飞行时，极易遭遇结冰气象条件，导致结冰事故[1-2]。因此结冰严重威胁到了直升机的飞行安全，需要引起高度重视。

直升机旋翼结冰过程相对于固定翼结冰而言非常复杂。随着计算流体力学技术的逐渐发展，相关理论的完善和成熟，数值模拟方法模拟旋翼桨叶结冰问题并分析其对气动特性的影响成为了可能。在结冰对直升机气动特性影响数值模拟方面，美国NASA于2012年开发的LEWICE软件增加了直升机旋翼结冰模拟模块[3]，首先模拟三维旋翼流场，然后截取其中的二维翼型流场结果，在此基础上采用固定翼翼型结冰数值计算方法预测冰形，最后将多个截面的计算结果组合并拓展至三维，利用三维结冰外形求解结冰后的气动特性。目前该方法已经应用在UH-1H直升机悬停的试验程序中。在国内，胡立芃[4]计算了旋翼桨叶不同位置翼型的结冰后气动特性变化，基于动量-叶素理论计算了旋翼气动特性，分析了结冰对直升机飞行性能的影响。胡国才等[5]分析了桨叶结冰对尾桨气动特性的影响，根据某直升机尾桨桨叶结冰实验结果，分析了桨叶翼型结冰前后气动特性，利用动量叶素理论计算了尾桨结冰前后的气动特性，结果表明结冰使得升力下降阻力上升，拉力系数下降而扭矩上升，气动性能恶化。李国知、曹义华等[6-7]以旋翼翼型冰风洞试验结果为基础，考虑了桨叶结冰脱落，研究了结冰对翼型气动特性的影响。 基于动量-叶素理论和涡流理论，以UH-60A为模型建立了旋翼结冰后的直升机飞行动力学模型，研究了直升机结冰后的平衡特性和稳定性的影响，分析了结冰时间、温度、LWC、MVD对旋翼结冰后的直升机操纵特性、姿态敏捷性、轴间耦合特性以及垂直轴操纵功效的影响。杨虎[8]以旋翼翼型冰风洞的试验数据作为依据，研究了结冰对飞行力学的影响，对直升机在旋翼结冰后的平衡性、稳定性进行了综合性的分析。刘国强等[9]分析了旋翼结冰对黑鹰直升机飞行性能影响，基于西科斯基公司翼型冰风洞结冰试验推导的工程公式，计算了翼型结冰后的气动特性；在此基础上建立了黑鹰直升机悬停和前飞计算模型。同时分析了桨叶结冰随飞行速度、温度和LWC变化规律，沿桨叶展向结冰厚度的变化情况，以及桨叶不同位置处空气动力对结冰的抑制作用[10]。可以看出，目前直升机旋翼桨叶结冰后气动特性数值计算主要是以二维翼型计算为主，由于目前数值方法较难直接获得桨叶三维结冰外形，较少针对三维桨叶直接开展结冰后气动特性分析。

本文在前期三维旋翼结冰数值模拟研究的基础之上，采用CFD方法建立了三维旋翼桨叶结冰数值模拟方法，针对结冰后的三维几何计算了结冰前后桨叶的气动特性，分析了结冰对直升机旋翼桨叶升阻特性以及悬停效率等的影响。

1 计算模型

对于旋翼桨叶流场的计算，由于桨叶相对于地面坐标系是旋转运动的，如果直接以地面坐标系为参考系进行计算，那么桨叶的运动是一个非稳态的运动，这不仅涉及到动网格的问题，而且计算量和难度都很大。

为了解决旋翼桨叶旋转问题，本文采用多参考系模型(MRF)对旋翼流场进行求解。MRF模型将计算区域几何划分为多个子区域，每个子区域有自己的运动方式，求解区域之间信息通过交界面进行传递。将整个区域划分为2个部分，将包围桨叶的内场区域定义为旋转坐标系区域，远离叶片的远场区域为地面坐标系。将坐标系定义与桨叶一同旋转则桨叶周围的流场将是一个定常状态。同时考虑到远场区域的流场是静止状态，所以将远场区域定义为地面坐标系。整个计算域示意如图 1所示。

图1 旋翼桨叶流场计算域

2 结冰数值计算方法

三维旋翼结冰数值模拟流程包括网格划分、空气流场计算、水滴撞击特性计算、结冰生成计算和几何模型重构。基于得到的三维几何外形就可以计算结冰前后旋翼桨叶的气动特性。

2.1 空气流场计算

采用旋转坐标系下以绝对速度表示的N-S方程计算旋翼空气流场，方程可以表示为以下形式：

(1)

(2)

(3)

2.2 水滴撞击特性计算

基于空气流场计算的结果，本文采用欧拉法计算水滴在旋转过程中的运动轨迹，对水滴列出质量，动量方程，可以得到如下方程：

(4)

Kαρ(ua-uw)+αρF

(5)

式(4)～(5)中：ρw为水滴密度；uw为水滴速度；F为所受外力，这里表示重力加速度；α为水滴容积分数；K为惯性因子。

水收集系数β是表征水滴撞击特性的重要指标，其定义为：微元上实际收集到的水量与可能收集到的最大水量之比，基于欧拉方法的水滴收集系数定义为：

(6)

式中：αn表示当地的水滴容积分数；α表示来流的水滴容积分数；V和n分别表示当地的水滴速度矢量和表面法相矢量；V表示水滴来流速度，对于旋翼结冰计算来说，水滴是相对旋翼桨叶运动进而撞击在桨叶表面，其来流速度应由转速来表征，因此可以采用桨尖速度表示水滴来流速度。

2.3 结冰生成计算

采用Messinger控制容积方法，将表面划分为多个控制体，对表面每个控制体列出质量和能量守恒方程：

(7)

(8)

2.4 几何模型重构

在数值模拟中网格划分过程会将原本的几何模型离散成众多的网格点，在此基础之上开展计算。在结冰数值模拟中，如何将结冰后的外形重构为完整几何成为计算结冰后气动特性的重要一环。

根据结冰后三维桨叶表面坐标及其拓扑关系，通过网格相邻单元排布将网格节点按照弦向方向和展向方向排列，采用双三次非均匀B样条曲面重构算法，重新构建新的三维桨叶几何外形[11]。随后利用得到的几何模型计算结冰后旋翼桨叶的气动特性。

本文研究重点在于分析结冰前后旋翼桨叶的气动特性变化情况，对旋翼结冰数值模拟方法的详细介绍和模拟方法正确性验证可以参考作者前期研究成果[12-13]。

3 计算域和网格划分

计算采用的旋翼桨叶为Caradonna-Tung(C-T)桨叶，桨叶的展长为1.143 m，弦长为0.190 5 m，展弦比为6，桨叶剖面形状为NACA0012翼型。流场的边界需要对桨叶的流场没有影响，所以计算域的出入口要与桨叶保持一定的距离。将流场计算域入口定义在距桨叶前缘20倍旋翼桨叶展长处，而流场计算域出口边界定义在后缘15倍展长处，开口边界距翼型为20倍展长。在模拟过程中，定义入口边界的气流速度，出口边界为压力出口，其余边界可使气流自由的进出。考虑到C-T旋翼有2片桨叶，为了减少网格量以及计算量，采用了周期性边界条件。

采用结构网格对桨叶进行划分，生成的计算网格如图 2所示。为了更好地计算桨叶的气动参数，对近壁面网格进行了加密。为了验证网格无关性，以未结冰旋翼桨叶为例，采用不同数量的网格开展流场计算。总共选取了3 240 000，5 520 000和6 240 000这3种数量的网格，从表1可以看出，随着网格量的逐渐增多，总扭矩逐渐降低，这说明网格量对计算结果有很大影响。5 520 000和6 240 000网格计算出的旋翼桨叶的总扭矩分别为47.8 N·m和47.6 N·m，总体结果相差不大。为了减少计算量，同时保证计算结果的准确性，选择5 520 000网格量来绘制网格。保持结冰桨叶和未结冰桨叶网格节点分布保持一致，因此在进行结冰前后桨叶的气动性能对比时可以消除网格量的影响。

表1 网格无关性验证

图2 旋翼桨叶网格划分

4 计算结果及分析

本节利用数值模拟方法得到结冰后三维旋翼几何外形，计算了结冰前后旋翼空气流场，对比分析了结冰前后旋翼气动性能。旋翼转速Ω=1 500 r/m，桨尖马赫数Matip=0.520，桨距角为8°。

整体气动性能对比如下。图3给出了悬停状态下结冰前后旋翼桨叶上表面压力分布云图。从图中可以看出，结冰后桨叶前缘出现了较为明显的负压区域，这是由于旋翼结冰破坏了桨叶原有的气动外形，结冰外形遮挡改变了气流流动情况，导致冰形后方出现了较为明显的负压区域，负压区域的出现会对旋翼升阻力产生较大影响。

图3 结冰前后桨叶表面压力分布的影响对比

表2中给出了结冰前后旋翼的气动特征参数。从表中可以明显看出结冰后旋翼的升力减小了8.8%，扭矩增加了11.1%，升力减小明显，扭矩急剧增大。衡量旋翼系统气动效率的重要参数悬停效率从0.670 8降低到0.525 7，降幅达21.6%，说明结冰对旋翼气动性能产生较大的影响，会导致直升机整体性能大幅下降。本文在计算过程中假设结冰前后旋翼实度不变，如果考虑结冰后旋翼实度增加，那么旋翼气动特性下降将更为剧烈。

表2 结冰对旋翼悬停性能的影响

图4为旋翼桨叶各截面拉力系数沿展向长度的分布情况，可以看出，相比于干净旋翼，结冰后桨叶的拉力系数出现了下降，特别是在桨叶尖部下降幅度较大，说明结冰对各个截面的拉力都产生了影响。

图4 不同位置桨叶拉力系数沿展向分布的对比

图5给出了不同截面处桨叶压力系数分布曲线，截面位置分别为r/R=0.70,0.80,0.91,0.96。从图中可以看出，结冰前后压力系数分布发生较明显的变化，结冰后压力系数的最大值虽然增加，但是在压力系数峰值后方压力迅速减小，导致此处出现了明显的压差阻力，使得旋翼阻力大幅增加，而整个压力系数曲线包围的面积也有所减小，说明升力有所减小。这是由于结冰改变了桨叶外形，导致冰形后方出现了明显的负压区，使得气动特性恶化。从图 6和图 7中的压力分布也可以看出结冰前桨叶表面的形状尚未被破坏，此时流线紧贴壁面。结冰后由于桨叶表面的形状被破坏，在前缘冰形后方出现一个负压区域，流线中也可以看出此处存在涡，这导致桨叶阻力增加，升力减小，且这种趋势随着展向位置逐渐增强。

图5 不同截面处压力系数曲线对比

图7 结冰后各个截面压力分布情况

5 结语

本文对直升机旋翼桨叶结冰后气动特性进行了研究，首先给出了一套三维旋翼桨叶结冰数值模拟方法，对其中的流场计算、水滴运动、结冰生成和几何重构进行介绍。随后计算了结冰前后旋翼桨叶的气动特性，分析了结冰对三维桨叶气动性能的影响。计算结果表明：

1)本文方法可以针对三维旋翼桨叶定量分析结冰对气动特性的影响。结冰会使桨叶整体升力下降，阻力上升，悬停效率下降超过20%。

2)结冰后旋翼桨叶各个截面的升力系数均发生降低，前缘处压力系数波动比较大，压力系数会有一个突变的过程，阻力明显增加。

3)结冰后桨叶原有气动外形被破坏，气流流线发生分离，冰形后方甚至产生了涡，导致了气动特性的恶化。