考虑转捩影响的风力机翼型气动特性计算研究

侯银珠 ,宋文萍 ,张 坤

(西北工业大学翼型、叶栅空气动力学国家级科技重点实验室,陕西 西安 710072)

0 引 言

风力机依靠风轮叶片获取风能,叶片的气动性能好坏直接影响着风力机的总体效率。构成叶片气动外形的翼型决定着叶片的性能,因此准确计算风力机翼型气动性能是进行叶片性能分析和气动设计的重要基础。传统的风力机通常选用NACA系列的航空翼型,这些翼型对前缘粗糙度非常敏感,前缘受到污染会导致翼型性能大幅度下降[1]。鉴于航空翼型不能很好地满足风力机及其特殊运行环境的要求,80年代后期风电发达国家开始对风力机翼型进行研究,并成功开发出风力机叶片专用翼型,比如美国NREL的S系列、丹麦 RIS Ø-A 系列、瑞典 FFA-W 系列和荷兰DU系列[2],目前DU翼型已被直径从29m到100m、最大功率从350KW到3.5MW的十多种不同型号的风电机组所采用[3]。

设计风力机翼型的技术基础大多是从60年代后期发展起来的,随着计算机硬件及软件技术的发展,先进的CFD技术已广泛应用于不同类型气动外形的计算。针对低雷诺数、高升力系数状态下的风力机运行条件,采用考虑粘性的NS控制方程分析叶片翼型的流场是非常必要的。由于风力机的工作环境和运行工况很复杂,来流方向不确定,风速大小多变及叶片表面存在的粗糙度,都会导致叶片上的流动更容易发生从层流向湍流的转捩。深刻认识翼型表面边界层内的流态,对于正确预估翼型升力和阻力、控制并减小流动分离以及翼型的优化设计有着重要意义。然而,由于对转捩机理认识的不足,多数对风力机翼型的绕流数值模拟都采用全湍流计算,难以准确地预估翼型表面边界层中转捩、分离等复杂现象。

为了更可靠地预测风力机翼型的气动性能,本文在求解RANS方程的过程中耦合Mark Drela[4]提出的简化eN转捩判断方法进行转捩点位置的自动判断,在翼型绕流数值计算中考虑流动转捩因素提高现有流场求解器的计算精度。

1 计算方法

在NS方程中耦合eN方法进行转捩点判断是一个迭代的过程,每一次转捩点位置的迭代需要调用三组程序。

首先是绕风力机翼型流动的计算程序,其功能是给层流边界层方程的求解提供物面压力分布。在连续介质假设下,忽略彻体力和热源的二维雷诺平均NS方程可表述为如下形式:

其次是二维不可压缩层流边界层方程求解程序,其功能是由RANS方程提供的物面压力分布确定外边界并求解出翼型表面层流边界层的解。本文采用二维不可压缩层流边界层方程,通过Falkner-Skan坐标变换转变原边界层方程的形式以降低边界层方程的求解难度,采用Keller提出的“BOX”格式离散求解边界层方程。

再者就是Mark Drela提出的的eN数据库转捩判断方法,这一方法将层流边界层的动量厚度、位移厚度及形状因子等参数与扰动的放大因子N联系起来,并给出了他们之间具体的关系,因此,利用边界层的解可以快速的得到边界层流动发生转捩的位置。eN数据库的具体形式为:

动量厚度雷诺数Reθ和空间扰动放大因子N的关系为:

临界动量厚度雷诺数Reθ可由形状因子来确定:

那么空间扰动放大因子N可由下式得到:

本文利用上述程序利用如图1所示流程实现在雷诺平均NS方程的迭代计算中耦合转捩点的自动判断,具体过程参见文献[7]。

2 算例及结果分析

根据上述方法本文对马赫数0.1,雷诺数1.0×106,攻角范围为-15.8°～ 28.1°的 DU91-W2-250翼型绕流分别进行了全湍流和耦合转捩判断的数值模拟计算,翼型几何形状如图2所示。计算网格为贴体的512×128 C型结构化网格,外边界距离翼型表面为15倍弦长,其中翼型上下表面共分布384个点,对翼型前后缘和表面处的网格进行了加密处理,图3所示为网格局部放大图。

图1 转捩点位置的自动判断过程Fig.1 Aotomatic transition prediction process

图2 翼型几何形状Fig.2 Airfoil profile

图3 计算网格放大图Fig.3 Computation grid with a local zoom

图4给出了采用全湍流和耦合转捩判断两种方法计算得到的DU91-W2-250翼型升力、阻力特性曲线并与实验值的比较结果,图中实验数据记为“Experimental”,全湍流计算结果记为“Cal-Fully Turbulent”,加转捩判断的计算结果记为“Cal-Transitional”。从图中可以明显看出,在风力机翼型的使用工况攻角变化范围(-10°～10°)内,加转捩判断的计算结果与实验数据基本吻合,而全湍流计算的升力系数在-5°～10°攻角范围内比实验值要低0.1左右,阻力系数在攻角-10°～10°范围内也略高于实验值,这说明在考虑流动转捩因素后的雷诺平均NS方程计算对翼型的升阻力特性的预测精度都有了一定提高。同样可以看到当攻角小于-10°或大于10°后,计算结果与实验值变化趋势基本一致,但计算的升力系数明显偏高而阻力系数偏小,这是由于数值计算所采用的湍流模型与真实情况还存在差异,高估了升力,低估了阻力,难以准确模拟大分离流动。

图4 翼型升、阻力系数随攻角变化曲线Fig.4 Curves of lift and drag coefficient versus angle of attack

图 5分别给出了翼型攻角为7.686°、9.742°时采用全湍流和耦合转捩判断计算的表面压力系数分布与实验数据的对比。从图可以看出,计算结果与实验测量结果基本一致,尤其是在攻角为7.686°、9.742°时,加转捩判断的压力分布与实验值符合很好,而全湍流计算的上翼面前半段压力分布与前两者相差较大,表明考虑转捩影响后的雷诺平均NS方程求解能更准确地模拟翼型绕流流场。

图5 翼型表面压力分布计算与实验结果比较Fig.5 Comparison of the pressure distribution between the calculated and the experimental data

3 小 结

本文以风力机翼型DU91-W2-250为研究对象,通过求解雷诺平均NS方程对其绕流进行数值模拟,计算采用SA湍流模型,并耦合eN数据库方法进行流动转捩的自动判断。采用该方法计算了雷诺数1.0×106时,攻角在-15.8°～28.1°范围内该翼型的升力和阻力特性曲线及表面压力分布,并与全湍流计算结果和实验数据比较,对比结果显示:

(1)在正常工况攻角范围内,耦合转捩判断的计算结果与实验结果相比符合很好,与全湍流假设得到的计算结果相比与实验值更吻合,说明要准确模拟风力机翼型使用工况下的气动特性,必须考虑转捩带来的影响;

(2)随着攻角不断增大,转捩点位置提前到翼型前缘,导致翼型的整个上表面几乎全是湍流,此时再考虑流动转捩因素进行转捩判断的计算与全湍流计算的效果相同。

[1]叶枝全,黄继雄,陈严.适用于风力机的新翼型气动性能的实验研究[J].太阳能学报,2003,24(4):548-554.(YE Z Q,HUANG J X,CHEN Y.Experimental study on aerodynamic performances of new airfoils for windturbine[J].Acta Energiae Solaris Sinicas,2003,24(4):548-554.)

[2]BERTAGNOLIO F,SORENSEN N,JOHANSEN J,FUGLSANG P.Wind turbine airfoil catalogue[R].RisØ National Laboratory,Roskilde,Denmark,August 2001.

[3]TIMMER W A,ROOIJ R P J O M.Summary of the delft university wind turbine dedicated airfoils[R].AIAA-2003-0352.

[4]DRELA M,GILES M B.Viscous-inviscid analysis of transonic and low reynold number airfoils[J].AIAA Journal,1986,25(10):1347-1355.

[5]韩忠华.旋翼绕流的高效数值计算方法及主动流动控制研究[D].西安:西北工业大学,2007.(HAN Z H.Efficient method for simulation of viscous flows past helicopter rotors and active flow control[D].Xi'an:Northwestern Polytechnical University,2007.)

[6]SPALART P R,ALLMARAS S R.A one-equation turbulence model for aerodynamic flows[R].AIAA-92-0439.

[7]张坤,宋文萍.基于转捩自动判断的翼型阻力精确计算方法研究[A].第六届西北地区计算物理学术会议论文集[C].2008:274-279.(ZHANG K,SONG W P.Accurate airfoils drag calculations based on the automatic transition prediction[A].In:6th computational physics conference papers of northwest in China[C].2008,pp.274-279.)