考虑遮蔽区影响的旋翼三维水滴撞击特性计算新方法

陈希， 招启军

南京航空航天大学 直升机旋翼动力学国家级重点实验室， 南京 210016

考虑遮蔽区影响的旋翼三维水滴撞击特性计算新方法

陈希， 招启军*

南京航空航天大学 直升机旋翼动力学国家级重点实验室， 南京 210016

针对直升机旋翼三维黏性流场特有的复杂环境，建立了一种基于欧拉法的旋翼三维水滴撞击特性计算的新方法。首先，在旋翼桨叶嵌套网格的基础上，发展了一套用于预测旋翼绕流流场的计算流体力学(CFD)模拟方法。然后，为克服传统直升机旋翼二维水滴撞击特性计算方法的不足，充分考虑旋翼流场的三维效应，在嵌套网格中基于欧拉法求解旋翼三维水滴撞击流场。其中，为解决尾流等区域的密度脉冲现象所引起的稳定性和收敛性问题，提出并建立了遮蔽区扩散模型。该模型通过判断遮蔽区变量，在计算域中动态生成遮蔽区域，并随迭代步数逐渐扩散。最后，通过与NACA0012翼型及国外UH-1H桨叶的试验和计算结果的对比，验证了旋翼三维水滴撞击特性计算新方法的可靠性，并进行了温度和水滴当量直径(MVD)对旋翼三维水滴撞击特性的影响分析。结果表明：遮蔽区扩散模型的加入，使二维情况的计算时间减少了22%，并增加了三维情况的计算稳定性，显著提高了旋翼三维水滴撞击特性的计算效率；沿着旋翼桨叶展向位置增大的方向，旋翼桨叶剖面水滴撞击范围有所增大，最大水滴局部收集系数呈先增加后减少再增加的变化趋势，其变化幅度接近50%；旋翼桨叶表面的水滴撞击区域和水滴局部收集系数随水滴当量直径的增加而增加。

旋翼； 结冰； 水滴撞击特性； 欧拉法； 遮蔽区； 计算流体力学

直升机在结冰环境下飞行时，空气中的过冷水滴会撞击到桨叶表面，造成旋翼结冰现象。旋翼结冰是一个危害直升机飞行安全的重要因素，它会改变旋翼绕流流场，破坏旋翼气动性能，进而影响直升机的操纵性，甚至导致空难的发生[1-3]。由于过冷水滴的运动轨迹决定了旋翼表面水滴局部收集系数，直接影响结冰情况，所以水滴撞击特性计算是直升机旋翼结冰数值模拟中的一个重要环节。

目前，水滴撞击特性的计算方法主要有拉格朗日法和欧拉法[4-8]。针对旋翼等复杂流动环境，拉格朗日法的实现过程较繁琐，计算效率和精度也需要提高。同时，欧拉法可以应用于非定常流场，固壁边界可以移动或转动，更适合于旋翼桨叶的水滴撞击特性计算[9-10]。

在进行旋翼水滴撞击特性计算时，研究者往往选择忽略旋翼的三维流动效应[11-14]。Narducci和Kreeger在旋翼结冰数值计算中，通过提取特征剖面的迎角和来流速度，对各剖面进行了二维的水滴撞击特性计算[11]。这种简化的方法给旋翼水滴撞击特性计算和结冰预测提供了便利，但是由于忽略了旋翼三维流场效应，计算得到的水滴撞击特性将会与真实值产生一定的偏差，正如文献[11]中的结果表明，桨叶中前段的结冰预测结果与试验数据的误差较大。因此，为获得更准确的旋翼水滴撞击特性，三维流场效应需要被充分考虑。

在欧拉法的水滴撞击特性求解过程中，由于控制方程的高度非线性，Bourgault等建立了基于体积分数的数值扩散项，以解决易出现的体积分数的数值振荡问题[15]。Tong和Luke将其归因于局部密度脉冲引起的高表观密度(Apparent Density)区域，并提出了对应的数值扩散系数[16]。杨胜华等进一步指出高表观密度区域易出现于流经机翼上下表面的两股气流交汇处，并提出自适应的数值扩散项，以解决计算不稳定、局部奇异解等数值问题[17]。在对固定翼的水滴撞击特性计算中，研究者往往将背风面的壁面密度设定为极小值[8,18]，且由于相对来流速度较快，局部密度脉冲现象引起的高密度区域可以很快被抹平。但是，旋翼三维水滴撞击流场具有桨尖涡干扰、相对桨叶的来流速度沿展向变化大、非定常流动等特征，气流交汇区域将更为频繁，局部密度脉冲现象所引起的高密度区域也将难以在迭代过程中被顺利抹平，给旋翼三维水滴撞击特性的计算增加了难度。

鉴于此，为获得更准确的旋翼水滴撞击特性，本文建立了适用于旋翼桨叶的三维水滴撞击特性计算新方法，并提出一种遮蔽区扩散模型以解决局部密度脉冲对计算稳定性的影响。为了验证本文方法的有效性，采用了美国直升机结冰飞行试验(The Helicopter Icing Flight Test)的旋翼结冰数据[19]。通过对比，遮蔽区扩散模型在能够保证计算精度的同时，可有效避免局部密度脉冲对计算稳定性的影响，并能够节省计算资源的消耗，提高了水滴撞击特性计算环节的效率。在此基础上，针对不同的环境温度和水滴当量直径的参数条件，分别对三维旋翼水滴撞击特性进行了计算和分析，获得了一些有意义的结果。

1 旋翼三维水滴流场的求解

在水滴撞击特性计算前，本文采用计算流体力学(CFD)方法对直升机旋翼三维流场进行了计算，其中，旋翼嵌套网格技术和旋翼CFD方法在文献[20]基础上进行了相应的改进。

1.1 水滴流场控制方程

在建立欧拉法求解水滴撞击流场前，先作如下假设[21]：① 水滴简化为直径为MVD(Median Volumetric Diameter)参数的球形，在运动过程中不会变形或被破坏；② 水滴之间没有碰撞或合并现象，水滴撞击在壁面上无飞溅现象；③ 水滴在空气流场中没有发生传热传质现象；④ 水滴上的作用力，只考虑重力和阻力，忽略其余非稳态力。

引入水滴体积分数(Water Volume Fraction)α的概念，它表示单位体积内水滴相所占有的体积，并由此构造水滴的表观密度ρd为

ρd=αρw

(1)

式中：ρw为水的密度。

通过水滴表观密度，建立水滴撞击流场的连续方程和动量方程。由于水滴与旋翼桨叶存在相对运动，将坐标系固连在旋转桨叶上。水滴撞击流场的积分形式控制方程为

(2)

源项由旋翼旋转引起的旋转通量和水滴上的作用力组成，它们的表达式为

(3)

式中：[uavawa]=qa为空气流场的绝对速度；Ki为惯性因子;ω为旋翼转速；g为重力加速度。

空气流场对水滴的作用力大小受惯性因子Ki的影响，惯性因子的表达式为

(4)

式中：μa为空气的动力黏度；dd为球形水滴的直径；ρa为当地空气密度；CDRe/24为Schiller-Naumann模型中定义的阻力函数，根据文献[17]的经验公式：

(5)

Re是相对雷诺数，其表达式为

(6)

在时间推进上，采用了五步显式Runge-Kutta迭代格式。为解决源项在控制方程中占优的问题，对时间步长进行了隐式处理。处理后的时间步长表达式为

(7)

式中：Δti为第i个单元的当地时间步长；m为五步Runge-Kutta法的迭代次数。

1.2 边界条件

水滴流场的边界条件设置如下：

1) 来流边界：qd=qa，α=LWC/ρw，LWC(Liquid Water Content)为液态水含量，表示单位体积内所含的液态水的质量。

2) 出流边界：所有变量法向一阶导数为0。

3) 壁面边界条件：与空气流场的固壁面不同，只允许水滴流出计算域，不允许水滴流入计算域。在不考虑水滴的碰撞飞溅现象后，水滴撞击在壁面时，会冻结成冰或以水膜形式吸附在壁面上，相当于水滴流出了计算域。

上述这种壁面边界条件可通过壁面法矢n与相对于桨叶的水滴速度矢量qd+qω的判定关系进行变换。若n·(qd+qω)≤0,则水滴撞击壁面，流出计算域；若n·(qd+qω)>0,则水滴未撞击壁面，以无滑移壁面边界条件处理。

1.3 局部收集系数的计算

水滴局部收集系数β为旋翼表面局部区域实际所收集的水量与该区域可能收集水量最大值的比值。该参数可以反映水滴的撞击范围和撞击区域内的收集水量分布，是影响最终结冰量的重要参数，其计算式为

(8)

式中：α∞为远场的水滴体积分数。

2 遮蔽区扩散模型

在旋翼三维水滴撞击流场求解过程中，某些局部单元会出现水滴表观密度的数值振荡，这可能会导致求解的不稳定。文献[16]将这种问题归源于欧拉法中对水滴速度值的单一假设。图1给出了沿桨叶表面运动的水滴路径图。

在现实情况中，水滴在空间中分布稀疏，较难发生相互碰撞的现象。在尾流处交汇的时候，它们往往会互不干扰，彼此保持原有的状态继续运动。这正如拉格朗日法中，追踪的不同水滴经常会经过同一单元格，彼此间并无影响。在欧拉法中，来自不同路径的水滴交汇于同一单元格时，由于该单元只能以一个速度矢量来表示，这就使原本流向不同方向的水滴聚在一起运动，从而造成了一个突然增加的密度脉冲。这种局部密度脉冲现象在翼型或固定翼的流场中出现较少，它们引发的高密度区域往往会在迭代过程中被及时抹平。而在三维旋翼水滴撞击流场中，尾流、桨尖涡、桨毂附近等气流交汇频繁的区域极易发生密度脉冲现象，从而影响计算的稳定性，产生数值振荡。

为避免上述情况的发生，本文提出了遮蔽区扩散模型。在迭代过程中，通过使用扩散模型，计算域内会逐渐在背风面生成遮蔽区，并向下游扩散。同时，这些区域会被划出计算域，减少了计算资源的消耗，并有效避免了密度脉冲带来的数值问题。

图1 沿桨叶表面运动的水滴路径图Fig.1 Paths of droplet moving along blade surface

在遮蔽区扩散模型中，引入遮蔽区权重变量shadow，变量shadow有如下定义：

1) shadow的初始化。在水滴撞击流场求解前，需要对其进行初值处理，赋值为0。

2) shadow在桨叶网格壁面层的更新。在每一次计算迭代完后，都需要对壁面层的shadow进行重新赋值，并与前文壁面边界条件的判断同时进行。当n·(qd+qω)>0，表示没有水滴撞击到物面区域，即遮蔽区的底层，赋值为2；当n·(qd+qω)≤0，表示有水滴撞击到物面区域，即桨叶表面的撞击区，赋值为0。由此，随着水滴撞击流场速度的更新，桨叶表面上的撞击区域和遮蔽区也在逐渐变化，并最终趋于稳定。

3) shadow的传递。当第P个单元格的shadowP为0时，且邻近单元存在shadowNi等于2时，赋值为1。这表示该单元紧邻遮蔽区，有成为遮蔽区的可能性。

4) shadow的增值与减值处理。当单元格shadowP为1时，则根据单元格的水滴表观密度的增/减分别对该单元格进行减/增值处理，其计算式为

(9)

式中：η为收敛因子，控制shadow值的变化速度，影响遮蔽区生成的速度大小，本文取值为0.1；αP为单元P的水滴体积分数；上标n和n-1分别表示t和t-Δt时刻的值。由此，在单元格内，水滴表观密度的减少量越大，单元格成为遮蔽区的可能性就越大，在扩散模型中体现为shadow值的增加。

5) shadow值的范围限定。在计算迭代过程中，为减少不必要的计算量，shadow值范围为[0，2]，2表示单元格成为遮蔽区的可能性最大，0表示可能性最小。

6) 遮蔽区的扩散判定。当单元格shadow值为2时，且满足式(10)的条件，则该单元标记为遮蔽区。

(10)

式中：ε为比例因子，即定义遮蔽区内单元与来流值的表观密度之比，可以控制遮蔽区生成的范围大小，本文取值为0.01；下标inf表示表观密度的初始值。shadowcri为遮蔽区标准值，其范围为(3,4]，本文取值为3.5，该数值的意义为：邻近单元格中存在一个可能成为遮蔽区的单元时(shadowNi>1.5)，当前单元格才能确定为遮蔽区。

图2给出了旋翼桨尖附近的遮蔽区在迭代计算中的扩散过程。为显示清晰，遮蔽区的单元格均隔行隔列显示。图中可见：遮蔽区首先依附于桨叶表面，生成遮蔽区的大致范围；随着迭代步数的增加，遮蔽区底层向桨叶前缘推进，厚度也有所增加；在遮蔽区底层停止推进后，遮蔽区开始迅速增厚，并向后缘扩散。

图2 遮蔽区在三维旋翼桨叶尖部的扩散过程Fig.2 Dispersion process of shadow zone over 3D rotor blade tip

3 计算结果与分析

3.1 NACA0012二维翼型的水滴撞击特性计算

本文对NACA0012二维翼型进行了水滴撞击特性的计算。翼型弦长c=0.533 4 m，来流马赫数Ma=0.302，迎角AOA(Angle of Attack)分别选取了0°、4°、8°，水滴当量直径MVD=20 μm，液态水含量为1 g/m3。

图3给出了不同迎角下的水滴表观密度ρd的分布云图。图中可见，在翼型背风面会形成较大范围的遮蔽区域，该区域内水滴表观密度很低。随着迎角的增加，上翼面的遮蔽区范围逐渐增加，下翼面的遮蔽区范围逐渐减小。

图3 不同迎角下的水滴表观密度等值图Fig.3 Contours of droplet apparent density at different AOAs

图4给出了不同迎角下翼型前缘的水滴局部收集系数的分布，横坐标中的S为该单元中心点沿翼型表面到驻点处的距离。由于在二维翼型的水滴撞击特性计算和结冰预测研究中，LEWICE软件发展成熟，本文将计算结果与文献中LEWICE软件的计算数据[8]进行了对比。图中可见，水滴局部收集系数的计算结果与参考结果吻合良好，验证了本文方法在二维翼型水滴撞击特性计算中的有效性。通过不同迎角的水滴局部收集系数的结果对比，可以再次发现，随着迎角的增加，水滴撞击范围逐渐增加。同时，最大水滴局部收集系数不受迎角变化的影响，均保持在0.7附近，当迎角增加后，最大水滴收集处将偏离驻点，向上翼面移动。

图4 不同迎角下水滴局部收集系数图Fig.4 Local collection coefficients of droplets at different AOAs

3.2 UH-1H三维旋翼桨叶结冰环境下水滴撞击特性计算

本文对UH-1H型贝尔直升机旋翼在悬停状态下的水滴撞击特性及结冰情况进行了预测。其中，桨叶参数、气象条件、试验数据和文献预测数据均来自于文献[19]。桨叶参数如下：翼型为NACA0012，翼型弦长为0.533 4 m，桨叶长度为7.315 2 m，桨叶转速为33.9 rad/s，桨叶扭度为10.9°。气象条件如下：压强为101 300 Pa，液态水含量为0.7 g/m3，水滴平均直径为30 μm，来流温度为-19 ℃。

在旋翼流场和旋翼水滴撞击流场的计算过程中，本文采用了同一套旋翼嵌套网格系统，如图5所示。背景网格采用半圆柱型网格，K=1和K=108为对称面。桨叶网格采用C-O型网格，并在桨叶前缘进行了加密处理。其中：I方向是指网格单元从桨叶剖面的尾迹开始沿翼型表面以顺时针方向遍历一周；J方向是指网格单元从翼型最底层向最外层遍历；K方向是指网格单元从桨根向桨尖遍历。

图5 悬停状态下的旋翼嵌套网格系统Fig.5 Embedded grid system for rotor in hovering flight

图6 旋翼桨叶三维水滴撞击特性图Fig.6 Sketch of 3-D droplet impingement on rotor blade

图6给出了通过欧拉法计算得到的水滴撞击特性的示意图。在桨叶前缘形成了水滴撞击区域，其中，撞击区域主要存在于桨叶的下表面；遮蔽区范围从撞击区域边缘开始，贴附于桨叶表面，向后缘延伸，并逐渐增厚，如图中水滴体积分数α的分布所示。

图7给出了旋翼桨叶表面的水滴局部收集系数分布图和3个特征剖面的水滴表观密度的分布云图。从整体上看，水滴局部收集系数在桨叶前半段沿展向呈增加趋势。从桨根到r=0.35R剖面处，水滴撞击面积和局部收集系数迅速增加，变化显著；在桨叶中段附近，水滴局部收集系数增加趋势减缓；在r=0.7R附近，水滴局部收集系数开始减少，又在r=0.85R(邻近桨尖)附近时再次增加。

图7 桨叶表面的水滴局部收集系数分布图 Fig.7 Distribution of local collection coefficient of droplets on blade surface

图8 旋翼三维水滴局部收集系数计算结果Fig.8 Calculated results of local collection coefficients of 3-D droplets on rotor blade

为了进一步显示旋翼三维水滴撞击特性，图8 给出了水滴局部收集系数的数值曲线图。图8(a)给出了3个典型剖面的水滴局部收集系数的分布，可以看出最大水滴局部收集系数均分布在驻点附近，同时下翼面撞击区域较大，上翼面撞击区域较小；随着径向距离的增加，水滴撞击范围在弦向上有所增加。图8(b)给出了最大水滴局部收集系数沿展向的分布情况，从整体上看，最大水滴局部收集系数在展向上有近50%的变化幅度。其中，最大水滴局部收集系数在桨根附近增加速度较快，从0.6增加到了0.85；在0.35R到0.7R间，最大水滴局部收集系数增加量较少，仅增加了0.1；随后，最大水滴局部收集系数先于0.83R处减少到0.6，又在桨尖处回升至0.85附近。

最大水滴局部收集系数在r=0.8R附近的下降现象主要源于旋翼桨尖涡的影响，图9(a)给出了旋翼流场的桨尖涡和水滴撞击流场的低表观密度区域(3.5×10-7kg/m3)。桨尖涡经过桨叶的尾流区，这导致了桨尖涡区域的水滴表观密度较低。同时，在桨尖涡的诱导作用下，大部分的水滴不经过桨尖涡区域，如图中剖面所示。由于旋翼做旋转运动，这个低水滴表观密度的桨尖涡区域会对后一片桨叶产生很大影响。图9(b)所示，在r=0.8R剖面的前缘附近，水滴表观密度较低，这导致了最大水滴局部收集系数在0.83R处的波谷。

图9 旋翼桨尖涡对水滴撞击流场表观密度分布的影响Fig.9 Influence of blade tip vortex on distribution of droplet apparent density in rotor droplet flowfield

图10给出了y=0平面的水滴表观密度的局部分布云图，图中ABCD区域内的黑色箭头表示水滴运动轨迹在y=0平面上的投影。在桨根附近(图中A所示)，水滴运动方向变化显著，具有明显的三维特征，如果采用二维翼型剖面的简化计算方法，可能会与真实值产生较大误差。在r=0.45R附近(z=3.2 m，图中B所示)，水滴运动轨迹仍具有明显的方向变化，后文选择了这里的剖面冰形进行了对比验证。在靠近桨尖附近(图中C、D所示)，由于相对来流速度的增加、展向速度的影响减小，水滴运动轨迹的方向主要沿x轴方向，此时，与二维翼型的水滴运动轨迹特点一致。

由于该试验中并未直接记录桨叶表面的水滴局部收集系数的试验值，这里通过计算冰形与试验值的对比来间接验证三维水滴撞击特性计算新方法的可靠性。在结冰计算环节中，采用了文献[22]的结冰预测方法。图11给出了展向最大结冰厚度的对比(两种测量方法的试验值来自文献[19]，结冰时间为180 s)。图12给出了桨叶剖面的结冰预测结果(结冰时间为180 s，计算冰形来自文献[11]的计算结果，试验冰形来自文献[19]中的Tracings方法)。由图11可见，在径向最大结冰厚度的对比上，计算值沿径向的变化趋势与试验值一致，且在各剖面上，预测冰形与试验值均吻合良好。由于旋翼桨叶的气动加热，桨叶的表面温度沿展向逐渐上升。剖面r=0.45R处的温度较低，所结冰型接近于霜冰类型，结冰量可以在一定程度上反映水滴收集系数的大小。因此，该剖面冰形可以用来验证水滴局部收集系数的计算结果。由图12(a)可见，在结冰范围和结冰量上，计算结果均与文献[19]的试验值一致，并优于文献[11]的计算结果。需要指出的是，文献[11]的计算结果是基于二维剖面的水滴撞击特性计算的方法，而r=0.45R剖面处于图10的区域B，具有较明显的三维效应。因此，通过与文献[11]计算结果在r=0.45R的剖面冰形对比，进一步说明了旋翼三维效应在水滴撞击特性计算中的重要性。

图10 y=0剖面上的水滴表观密度分布图Fig.10 Distribution of droplet apparent density in y=0 plane

图11 旋翼最大结冰厚度分布图Fig.11 Maximum depth of ice accretion on rotor

图12 不同径向位置的剖面冰形图 Fig.12 Ice shape on different sections for different radial position

3.3 遮蔽区扩散模型在水滴撞击特性计算过程中的影响

为了讨论遮蔽区扩散模型对水滴撞击特性计算的影响，分别给出了二维翼型和三维旋翼桨叶情况下的典型算例。

图13给出了前文AOA为0°时二维翼型算例的计算情况。在加入遮蔽区模型后，残值在580步附近即达到10-3量级，比未加入模型情况提前了约150步，节省了近22%的计算时间。同时，通过水滴局部收集系数的对比，可以看出，遮蔽区模型的引入，并不会影响二维水滴撞击特性的数值计算结果。

图14(a)给出了遮蔽区单元的动态生成过程。图中可见，在计算初期，遮蔽区单元生成速度较快，随着迭代步数的推进，遮蔽区单元增加速度逐渐减少，在约500步附近，达到最大值。之后，遮蔽区单元个数将保持稳定，遮蔽区也不再扩散。图中也分别给出了第200步、第300步和第500步的翼型尾缘的局部遮蔽区域分布情况。图14(b)给出了遮蔽区域整体的动态扩散过程。在第200步时，遮蔽区主要依附于翼型表面，此时，遮蔽区在整个计算域中所占比例较小；在第300步时，翼型表面的遮蔽区域逐渐增厚，并有向后方扩散的趋势；在第500步时，遮蔽区单元进一步向尾流方向扩散，同时，翼型表面遮蔽区的厚度增加速度减缓。

图13 遮蔽区扩散模型对二维翼型计算结果的影响Fig.13 Influence of shadow zone model on calculated results of 2-D airfoil

图14 遮蔽区单元扩散过程Fig.14 Dispersion process of cells among shadow zone

图15是三维旋翼桨叶的水滴撞击流场残值对数收敛曲线。在未加入遮蔽区模型时，残值在第55迭代步开始发散，不能继续进行计算，这正是因为尾流附近发生了密度脉冲。在加入遮蔽区扩散模型后，由于三维桨叶遮蔽区域的底层生成数量较多，残值在经过初期发生的小规模振荡后，迅速收敛。

图15 遮蔽区扩散模型对旋翼三维算例结果的影响Fig.15 Influence of shadow zone model on calculated results of 3-D rotor

图16 未加入遮蔽区扩散模型的第55迭代步的计算结果Fig.16 Calculated results at the 55th iterative step without present model

考虑到未加入遮蔽区扩散模型的算例在第55迭代步后发散，图16给出了当前迭代步的水滴流场数据。由图可见，由于桨叶上下来流的交汇，在桨叶尾流后呈现大面积的高密度区域。在r=0.15R处剖面的水滴表观密度示意图中，尾流出现呈带状的高密度区域，其中，某些单元的表观密度高达来流值的40倍。

加入遮蔽区扩散模型后，密度脉冲引起的高密度区域面积减小，如图17所示。第55步时，与未加入模型的结果相比，高密度区域虽然存在，但面积和数值均有所降低。同时，因遮蔽区域的生成，从桨叶尾缘上脱出的低密度区域面积有所增加，可以有效抑制密度脉冲现象的发生，满足数值计算的稳定性要求。当迭代步数为100步时，由于遮蔽区域的进一步扩散，高密度区域消失。此时，在尾流区将不再出现因密度脉冲现象引起的高密度区域，增强了计算的稳定性。

图17 加入扩散模型后桨叶剖面表观密度及遮蔽区的计算结果(r=0.15R)Fig.17 Calculated results of apparent density and shadow zone of blade with present model (r=0.15R)

3.4 不同气象参数下的旋翼桨叶水滴撞击特性计算

为了分析气象参数对旋翼水滴撞击特性的影响，本文在不同的水滴平均当量直径和环境温度的情况下(其余气象参数与前文一致)，对悬停状态下的三维旋翼进行了计算。

表1给出了水滴当量直径(MVD)对水滴撞击特性的影响。从表1中看出，MVD为15 μm时，水滴撞击面积占整个桨叶表面积的2.6%，当MVD增加到30 μm时，撞击面积比例增加到6.75%，单位时间的水滴收集量也随之增加。可见，MVD参数对旋翼水滴撞击特性的影响很大。

表1 不同MVD下的水滴撞击特性Table 1 Droplet impingement property at different MVDs

图18给出了在不同MVD参数下的三维旋翼表面水滴局部收集系数的分布云图。图中可见，随着MVD参数的增加，水滴的惯性(保持原有运动方向的能力)随之增强，这使得水滴更容易撞击到旋翼桨叶表面，导致了水滴撞击区域扩大、水滴局部收集系数增加。

表2给出了环境温度对水滴撞击特性的影响。在易发生结冰的环境温度下，随温度的下降，水滴收集量和撞击面积都有小幅度的增加。总体上，水滴撞击特性受温度的影响较小。

图18 不同MVD参数下的水滴局部收集系数分布图 Fig.18 Distribution of local collection coefficient at different MVDs

表2不同温度下的水滴撞击特性

Table2Dropletimpingementpropertyatdifferenttemperatures

CaseT/℃Amountofimpingingwater/(10-3g·s-1)Proportionofimpactareatosurfacearea/%1-47．1826．692-87．1906．703-147．2466．754-197．2956．835-247．3296．89

4 结 论

1) 建立了一种基于欧拉法的三维旋翼水滴撞击特性计算新方法。通过对NACA0012翼型和UH-1H旋翼的验证计算，表明了该方法的有效性。同时，提出的遮蔽区扩散模型有效解决了密度脉冲所引起的稳定性和收敛性问题，并大幅减少了水滴撞击特性环节的计算时间。

2) 在二维翼型上，最大水滴局部收集系数受来流迎角的变化影响较小，主要受来流马赫数影响。在三维旋翼桨叶上，随着径向距离的增加(相对来流马赫数变大)，旋翼桨叶剖面撞击范围有所增加，最大水滴局部收集系数会呈现先增加后减少再增加的变化趋势。同时，在桨尖附近，水滴运动轨迹三维特征较小；在桨叶中段至桨根附近，水滴运动轨迹具有较明显的三维特征，用二维剖面拟合的方法在该部位可能会偏离真实值。

3) 随着MVD参数的增加，旋翼桨叶表面的水滴撞击区域扩大，水滴局部收集系数有所增加，并可能导致后续结冰量的增加。

[1] FLEMMING R J. The past twenty years of icing research and development at Sikorsky aircraft: AIAA-2002-0238[R]. Reston: AIAA, 2002.

[2] FLEMMING R J, ALLDRIDGE P, DOEPPNER R. Artificial icing tests of the S-92 helicopter in the McKinley Climatic Laboratory: AIAA-2004-0737[R]. Reston: AIAA, 2004.

[3] BRITTON R K, BOND T H, FLEMMING R J. An overview of a model rotor icing test in the NASA Lewis icing research tunnel: AIAA-1994-0716[R]. Reston: AIAA, 1994.

[4] WRIGHT W B, POTAPCZUK M G, LEVINSON L H. Comparison of LEWICE and GlennICE in the SLD regime: AIAA-2008-0439[R]. Reston: AIAA, 2008.

[5] IULIANO E, BRANDI V, MINGIONE G, et al. Water impingement prediction on multi-element airfoils by means of Eulerian and Lagrangian approach with viscous and inviscid air flow: AIAA-2006-1270[R]. Reston: AIAA, 2006.

[6] 胡剑平, 刘振侠, 张丽芬. 发动机整流支板大尺寸过冷水滴撞击特性[J]. 航空学报, 2011, 32(10): 1778-1785.

HU J P, LIU Z X, ZHANG L F.Supercooled large droplet impact behaviors on an aero-engine strut[J]. Acta Aeronautica et Astronautica Sinica, 2011, 32(10): 1778-1785 (in Chinese).

[7] 易贤, 王开春, 桂业伟, 等. 结冰面水滴收集率欧拉计算方法研究及应用[J]. 空气动力学学报, 2010, 28(5): 596-601.

YI X, WANG K C, GUI Y W, et al. Study on Eulerian method for icing collection efficiency computation and its application[J]. Acta Aerodynamica Sinica, 2010, 28(5): 596-601 (in Chinese).

[8] KIM J W, DENNIS P G, SANKAR L N, et al. Ice accretion modeling using an Eulerian approach for droplet impingement: AIAA-2013-0246[R]. Reston: AIAA, 2013.

[9] KINZEL M P, SAROFEEN C M, NOACK R W, et al. A finite-volume approach to modeling ice accretion: AIAA-2010-4230[R]. Reston: AIAA, 2010.

[10] BAIN J, SANKAR L N,GARZA D. A methodology for the prediction of rotor blade ice formation and shedding[C]//Proceedings of the SAE 2011 Aircraft and Engine icing and Ground Testing Conference. Warrendale, PA: SAE International, 2011.

[11] NARDUCCI R, KREEGER R E.Analysis of a hovering rotor in icing conditions: NASA/TM-2012-217126[R]. Washington, D.C.: NASA, 2012.

[12] CAO Y H, LI G Z, ZHONG G. Tandem helicopter trim and flight characteristics in the icing condition[J]. Journal of Aircraft, 2010, 47(5): 1559-1569.

[13] RAJMOHAN N, BAIN J, NUCCI M, et al. Icing studies for the UH-60A rotor in forward flight[C]//AHS Aeromechanics Specialists’ Conference. Fairfax, VA: American Helicopter Society International, 2010.

[14] 张威, 林永峰, 陈平剑. 基于欧拉方法的旋翼翼型结冰数值模拟及参数影响[J]. 南京航空航天大学学报, 2011, 43(3): 375-380.

ZHANG W, LIN Y F, CHEN P J. Numerical simulation of ice accretion and parameter effects based on Eulerian droplet model[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2011, 43(3): 375-380 (in Chinese).

[15] BOURGAULT Y, HABASHI W G, DOMPIERRE J, et al. A Eulerian approach to supercooled droplets impingement calculations: AIAA-1997-0176[R]. Reston: AIAA, 1997.

[16] TONG X L, LUKE E A. Eulerian simulations of icing collection efficiency using a singularity diffusion model: AIAA-2005-1246[R]. Reston: AIAA, 2005.

[17] 杨胜华, 林贵平, 申晓斌. 三维复杂表面水滴撞击特性计算[J]. 航空动力学报, 2010, 25(2): 284-290.

YANG S H, LIN G P, SHEN X B. Water droplet impingement prediction for three-dimensional complex surfaces[J]. Journal of Aerospace Power, 2010, 25(2): 284-290 (in Chinese).

[18] JUNG K Y, JUNG S K, MYONG R S, et al. A three-dimensional unstructured finite volume method for droplet impingement in aircraft icing: AIAA-2013-2576[R]. Reston: AIAA, 2013.

[19] LEE J D, HARDING R, PALKO R L. Documentation of ice shapes on the main rotor of a UH-1H helicopter in hover: NASA-CR-168332[R]. Washington, D.C.: NASA, 1984.

[20] 招启军, 徐国华. 新型桨尖旋翼悬停气动性能试验及数值研究[J]. 航空学报, 2009, 30(3): 422-429.

ZHAO Q J, XU G H. Aerodynamic performance of rotor with new type blade-tip in hover based upon test and numerical investigations[J]. Acta Aeronautica et Astronautica Sinica, 2009, 30(3): 422-429 (in Chinese).

[21] CAO Y H, MA C, ZHANG Q, et al. Numerical simulation of ice accretions on an aircraft wing[J]. Aerospace Science and Technology, 2012, 23(1): 296-304.

[22] 陈希, 招启军, 赵国庆. 计入离心力影响的直升机旋翼翼型结冰数值模拟[J]. 航空动力学报, 2014, 29(9): 2158-2165.

CHEN X, ZHAO Q J, ZHAO G Q. Numerical simulation for ice accretion on rotor airfoil of helicopter considering the effects of centrifugal force[J]. Journal of Aerospace Power, 2014, 29(9): 2158-2165 (in Chinese).

(责任编辑： 鲍亚平，徐晓)

AeronauticsandAstronautics,Nanjing210016,China

New method for predicting3-D water droplet impingement on rotorconsidering influence of shadow zone

CHENXi,ZHAOQijun*

NationalKeyLaboratoryofScienceandTechnologyonRotorcraftAeromechanics,NanjingUniversityof

For the complex 3-D viscous flowfield of the helicopter rotor, a new 3-D Eulerian method is established for calculating the droplet impingement on the rotor. The Computional Fluid Dynamics (CFD) simulation method for predicting the rotor flowfield is developed based on the embedded grid method. Considering the 3-D effect of the rotor fully, the 3-D flowfield of droplets on the same embedded grids is solved by the Eulerian method to overcome the defects of the traditional 2-D calculation method. To overcome the problem of numerical stability and convergence by the density impulse in the wake area, a shadow zone dispersion model is proposed. In the model, a shadow variable is used to control the generation and dispersion of the shadow cells. The new Eulerian method is validated by comparing the calculation and experiment results of NACA0012 airfoil and UH-1H rotor in hover. In addition, the effects of the atmospheric temperature and Median Volumetric Diameter (MVD ) of the droplet on the droplet impingement on the rotor are calculated and analyzed. Results show that the computation time of the 2-D droplet impingement property is reduced by 22% and the calculation stability of 3-D droplet impingement property is improved by using the shadow zone dispersion model. The droplet impingement area increases with the spanwise distance on the rotor blade. The maximum droplet local collection efficiency presents the change tendency of ‘increase-decrease-increase’, and the variation amplitude is close to 50%. The droplet impingement area and the collection efficiency increase with the increase of the MVD.

rotor; ice accretion; water droplet impingement property; Eulerian method; shadow zone; CFD

2016-09-01；Revised2016-10-11；Accepted2016-11-25；Publishedonline2016-12-051640

NationalNaturalScienceFoundationofChina(11272150)

2016-09-01；退修日期2016-10-11；录用日期2016-11-25； < class="emphasis_bold">网络出版时间

时间：2016-12-051640

www.cnki.net/kcms/detail/11.1929.V.20161205.1640.004.html

国家自然科学基金 (11272150)

*

.E-mailzhaoqijun@nuaa.edu.cn

陈希， 招启军． 考虑遮蔽区影响的旋翼三维水滴撞击特性计算新方法J． 航空学报,2017,38(6):120745.CHENX,ZHAOQJ.Newmethodforpredicting3-DwaterdropletimpingementonrotorconsideringinfluenceofshadowzoneJ.ActaAeronauticaetAstronauticaSinica,2017,38(6):120745.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0313

V211.3

A

1000-6893(2017)06-120745-12

URL：www.cnki.net/kcms/detail/11.1929.V.20161205.1640.004.html

*Correspondingauthor.E-mailzhaoqijun@nuaa.edu.cn