146 第二十六届（2010）全国直升机年会论文 桨叶结冰对旋翼气动特性影响的计算 胡立芃 刘国强 唐正飞 (南京航空航天大学直升机旋翼动力学重点实验室，南京，210016)
摘 要：用CFD、动量-叶素等方法建立一套桨叶结冰后旋翼气动特性分析方法。包括：二维翼型结冰的建模；二维翼型结冰后气动特性的计算方法；旋翼结冰对旋翼悬停性能的影响分析。本文使用UH-1的结冰实验数据验证了算法，并对国产直八直升机进行结冰前后的对比计算。计算结果表明积冰状态下桨叶翼型的升力系数降低，阻力系数增大。最终结果表明：积冰严重影响着旋翼的性能，这与由相关实验得出的结论一致。 关键词：结冰；直升机；旋翼
1引言 直升机结冰是一个老问题，也是公认的危及飞行安全的严重问题，特别是直升机旋翼系统结冰比固定翼飞机更加敏感，由于其自身具有可用功率有限，操纵面较小等特点，积冰更易使直升机造成危险。美国在1982～2000年统计，因结冰引起583起飞行事故，造成800多人死亡[1]。为了实行直升机全天候安全飞行的目标，各国针对直升机结冰问题做了大量研究。目前国内外针对旋翼桨叶结冰的研究方法主要有：在真实结冰气象条件下进行飞行试验；在人工气候实验室制造的模拟云中进行飞行测试；数值计算和风洞试验。如西科斯公司已经从事结冰研究以及防冰系统的研发超过57年，进行了直升机真实结冰气象的飞行试验，旋翼和翼型全尺寸模型以及缩比例模型的风洞试验。19世纪70年代，西科斯公司研制的电加热除冰系统，安装在2400架黑鹰直升机，现在仍然是一种有效的防冰系统。NASA也承担了结冰对旋翼翼型影响，直升机旋翼桨叶防/ 除冰系统等问题的研究，建立了实用的旋翼翼型结冰数据库,同时发展了LEWIC结冰计算程序。LEWIC可以计算直升机在结冰条件下的流场、水滴撞击特性，模拟结冰过程。国内也对结冰模型[2]、流场计算[3]、水滴撞击特性等方面做了大量研究，分析了结冰对翼型气动性能的影响，以及旋翼的防/ 除冰系统的防护范围[4]等。本文将采用数值计算方法，通过建立二维翼型结冰的数学模型，预测结冰后翼型形状，进而分析翼型结冰后的气动特性，采用动量-叶素结合方法分析结冰对旋翼气动特性影响。
2结冰后的外形确定 结冰的数值模拟主要分为：网格生成、空气流场计算、求解水滴运动方程、冰生长模型。由于过冷水滴在流场中的体积含量很小，不能够影响空气的流动,因此空气控制方程同水滴控制方程可以单独求解。流程如图1，首先计算空气流场，在求解出水滴粒子的撞击特性和运动轨迹后，建立单位控制体积内的质量和能量平衡冰生长模型，模拟结冰过程.当翼型形状改变后需要重新生成网格计算，直到所需的结冰时间为止。
tT结冰外形计算结果 t=0
总时间达到水滴收集系数
指定时间
初始外形生成计算网格空气流场计算结冰外形结 冰 模 型水滴流场求解
图 1 结冰外形计算流程图 2.1 空气流场计算
本文将采用FLUENT计算空气绕机翼运动的外流场，由于本文假设的结冰部分发生在靠近桨根 147
部分，空气速度较低，使用用定常、不可压纳维－斯托克斯（N-S）方程，积分形式如下： 0uuuxyz （2-1）
2222()()uuuupuvuxyxyx （2-2）
2222()()vvvvpuvuxyxyx （2-3）
EEQWxy
（2-4）
其中为流体的瞬时密度，u、v为流场速度在两个坐标方向的分量，E是系统的总能，Q为是通过系统界面以热传导形式传递给系统的热量。 压力-速度耦合采用SIMPLE算法[5]，在中心有限体积法基础上，采用中心格式离散N-S方程。湍流模型使用Spalart-Allmaras模型。流场计算结果再作为已知条件，用于水滴运动方程的求解。 2.2 水滴运动方程求解 采用欧拉参考系下的水滴连续性方程和动量方程[6]分别为：
（2-5） （2-6） 其中，水滴的容积分数，u水滴的速度矢量，ou空气的速度矢量，K空气-水滴交换系数:
218opfKd （2-7）
公式中：o空气动力粘度，pd水滴直径，f阻力函数,采用Schiller – Naumann模型时： （2-8） 其中，水滴阻力系数DC：
（2-9） ||oopeouudR （2-10）
o为空气密度，通过求解水滴连续方程与动量方程得到水滴的运动方程与水滴容积分数，
以及欧拉两相流法中,翼型表面局部水滴收集率：
()0tuu
()()otuuuKuu
24DeCRf
1000eR1000eR0.68724(10.15)/0.44C{DEeRR 148
ssuLWCu
（2-11）
其中LWC与u为来流的液态水含量和水滴速度，su与s为翼型表面的水滴法向速度与水滴的密度。 2.3 结冰模型 根据霜冰的形成过程，假设所有的水滴在碰撞后就完全凝结, 并且冰沿着与翼型表面法向一致的方向增长，可以只考虑质量守恒。根据水滴流场求解结果，一个时间步长内翼型表面控制体的水滴收集质量为：
TMLWCuAT （2-12）
其中A为该控制体的水滴撞击面积，T为时间步长。 当翼型表面结冰产生新的外形后，需要重新的结冰翼型进行流畅计算，然后求解水滴运动方程、计算结冰厚度，直到需要的时间为止。 2.4 算法验证 取翼型NACA0012初始计算条件为： 翼型弦长：C = 0. 533m；气流速度: U∞ = 119m / s； 空气液态水含量：LWC = 0. 75g/m3 ； 结冰时间：4.5min;迎角：6°； 水滴当量直径：deq = 20μm； 图 2 本文计算结冰外形 环境温度: T∞ = 262K 计算值与实验值的比较如图2~3： 图3为NASA结冰实验所得翼型结冰后外形，实验所用翼型以及结冰条件与本文相同。实验所获积冰形状与计算结果趋势一致，结冰厚度基本相同，从而验证了本文计算翼型结冰形状的准确性。 本文以直八直升机旋翼为例，研究结冰后旋翼气动特性。直八直升机旋翼转速：207 rpm；桨叶长度R=9.45m；桨叶片数：b=6；桨叶翼型：NACA0012；翼型弦长：C = 0. 54m；桨叶负扭转5°50 图3 NASA结冰实验结果 根据ＮＡＣＡ结冰实验以及常士楠等人直升机旋翼桨叶防/ 除冰系统防 护范围研究[4]，将桨叶分为5段：0~0.25R、0.25~0.35R、0.35~0..45R、 0.45~0.6R、0.6~1R,其中，0.6~1R段未发生结冰，取各段中间状态：分别 是沿桨叶展向２０％３０％４０％52.5％等处（依次定义为1#2#3#4#截面），计算环境条件为： 空气液态水含量：LWC = 0. 75g/m3； 结冰时间：4.5min;水滴当量直径：deq = 20μm； 环境温度：T∞ = 262K； 通过方程：4-1与4-2计算各个界面处的来流速度与迎角。 计算结果如图4~7所示：
图 4 桨叶0. 2R处结冰外形 149
图 5 桨叶0.3R处结冰外形 图 6桨叶0.4R处结冰外形 图7桨叶0.525R处结冰外形 3 结冰后翼型气动特性计算 结冰翼型的气动特性计算是旋翼气动特性研究的基础，建立翼型的气动特性的计算方法，对准确计算旋翼气动特性具有重要意义。本章将利用计算流体力学(CFD)方法计算翼型气动特性，将结冰后翼型与光滑翼型气动特性计算结果进行对比验证。 流场计算模型与本文2.1相同，计算得到的结冰前后升阻力系数变化如下：
图8为结冰前后阻力系数随迎角的变化曲线 图9为结冰前后升力系数随迎角的变化曲线 如图8、9所示，结冰对翼型的气动性能有明显的影响，这是由于结冰改变了翼型的外气动外形，气流在结冰翼型后方形成的旋涡对外流场产生了剧烈扰动。对于升力系数，在结冰翼型后方形成的低压区使得升力减小，随着迎角的增大影响越来越明显。由于翼型的摩擦阻力和压差阻力均变大，阻力系数也明显增大，且随着迎角的增大影响也越来越显著；各段翼型气动特性的降低，必然会导致整个桨叶气动特性的降低。
4 旋翼气动特性计算与分析 研究旋翼气动特性的分析方法主要包括：动量理论、叶素理论、涡流理论以及旋翼流场计算的CFD 方法。本文利用动量-叶素理论相结合的方法计算桨叶不同半径处的诱导速度以及翼型迎角，进而计算翼型气动特性，利用叶素理论得到旋翼的气动力、力矩和功率。 4.1动量-叶素理论计算模型
根据动量理论，距离桨盘中心为r，宽度为dr的环带桨盘的拉力增量为：214Tvrr，根据
叶素理论 21()()2vTbracrr