转捩对风力机翼型和叶片失速特性影响的数值模拟

钟 伟，王同光

(南京航空航天大学江苏省风力机设计高技术研究重点实验室，江苏 南京 210016)

0 引言

风力机翼型的气动特性是风力机叶片设计的基础输入参数。要获得翼型在多个雷诺数和较大迎角范围内的完整风洞实验数据，所需费用较高。计算流体力学(CFD)数值模拟技术逐渐成为获取和研究翼型气动特性的一种更经济快捷的选择。在众多风力机翼型中，S809翼型［1］是最重要的风力机空气动力学研究对象之一，因为它不仅风洞实验数据充分，而且是美国国家可再生能源实验室(NREL)开展的水平轴风力机非定常空气动力学系列实验［2］的叶片翼型。S.L.Yang等人［3］于1994年最早对S809翼型进行了全湍流数值模拟。在附着流动阶段，升力系数与实验结果基本吻合，只在个别迎角下出现压力分布局部不符;但阻力系数远大于实验值。Walter P.Wolfe等人［4］于1997年使用指定转捩点的方式对S809翼型进行了数值模拟。附着流动阶段的阻力系数误差相对全湍流模拟显著下降，压力分布局部不符的问题也不再存在。他们的结果表明，转捩对翼型附着流动状态下的气动特性有影响，值得加以研究。

风力机叶片的气动特性与其基于的翼型的气动特性有紧密联系。在对S809翼型的数值模拟研究基础上，研究人员对基于该翼型的NREL Phase VI叶片开展了大量的数值模拟研究。Phase VI实验于2000年春季在世界上最大的风洞——NASA Ames研究中心24.4m×36.6m全尺寸风洞中进行，是迄今为止可信度最高测量数据最全面的风力机空气动力学实验。在实验结果公布之前，NREL即组织了世界上多名风力机空气动力学研究人员采用多种方法进行了计算盲比。然而盲比的结果出乎意料的不理想［5］。其中与实验数据最接近的结果来自CFD数值模拟［6］，但仍然存在较大误差。此后，众多学者采用多种解算器和计算模型对Phase VI叶片进行了CFD数值模拟［7-9］。这些结果较NREL盲比结果有所改进，但距离与实验数据的较好吻合仍然有相当的距离，且至今未有对造成数值模拟误差的原因形成合理解释。以上数值模拟都是全湍流模拟，关于转捩对叶片气动特性影响的研究很少。

R.B.Langtry等人［10］于2006 年使用 F.R.Menter等人发展的 Gamma-Theta转捩模式［11］对 S809翼型和Phase VI叶片进行了数值模拟。他们预测的翼型转捩点位置与实验吻合，且捕捉到了层流分离和湍流再附形成的分离泡，阐述了转捩对S809翼型在附着流动状态下气动特性的影响，但还缺乏对翼型失速特性影响的研究。他们对Phase VI叶片进行的转捩模拟获得的扭矩值在20m/s风速下相对全湍流模拟更接近实验值，在其它风速下则与全湍流模拟差别不明显。但其在20m/s风速下全湍流模拟得到的扭矩值本身高于实验值太多，误差明显大于其他学者的数值模拟结果，这削弱了其关于转捩影响的结论的说服力。

本文在以上学者的研究基础上，使用Gamma-Theta转捩模式对S809翼型和Phase VI叶片进行了数值模拟，着重研究了转捩对翼型和叶片失速特性的影响。

1 数值模拟方法

1.1 控制方程

求解的控制方程组为守恒型非稳态不可压缩雷诺平均的Navier-Stokes方程，其在静态坐标系下的表达式为:

式中ρ为空气密度，U为平均速度矢量，u为脉动速度矢量，p为静压，为分子应力张量，ρu⊗u为雷诺应力。

采用有限体积方法进行CFD定常求解。方程离散采用二阶迎风格式。

1.2 湍流模型与转捩模式

湍流模拟采用切应力输运(shear stress transport，SST)模型。该湍流模型是公认的较好的线性涡粘性模型，对边界层流动和分离流动都能有较好的模拟效果。采用的转捩模式为Gamma-Theta转捩模型。该湍流模型和转捩模型均由F.R.Menter提出，被结合使用。受篇幅限制，本文仅列出Gamma-Theta转捩模型的输运方程。关于该湍流模型和转捩模型的详尽描述及本文所列输运方程中各参数的含义可参考相关文献［11］。

Gamma-Theta转捩模型中间歇因子γ和转捩动量厚度雷诺数～Reθt的输运方程分别为:

该转捩模型的求解顺序为:根据上一个时间步的平均流场及 γ值，通过的输运方程求出，然后通过间歇因子γ的输运方程求解出当前时间步的γ，最后通过有效粘性系数γut作用于平均流场。

1.3 多参考系模型(MFR)

在对风力机叶片的数值模拟中，叶片相对地面坐标系是旋转的，流场相对地面坐标系是非稳态流动。如果直接以地面坐标系为参考系进行求解，需要对叶片的实际运动进行非稳态模拟，且涉及网格的运动问题，计算量和难度都较静止物体的求解要大。如果定义参考系跟随叶片一起旋转，则在此旋转坐标系下流场转化为定常状态。因此，考虑采用旋转坐标系的方法来模拟叶片的转动。但为了正确地设置边界条件，又需要远场边界在计算坐标系下是静止的。为了解决模拟叶片旋转与正确设置远场边界条件这一矛盾，采用多参考系模型(MFR)，将计算域分为内外两个区域，包围叶片的内区采用旋转坐标系，靠近远场的外区则使用地面坐标系。

在实际解算过程中，处于旋转坐标系区域内的流体被赋予了一个与坐标系旋转速度相反的旋转速度。另外，还需要考虑离心力和科氏力的影响，在动量方程中添加相关源项。离心力和科氏力的源项分别为:

1.4 计算域、网格与边界条件

本文S809翼型的数值模拟不考虑翼型实验的风洞洞壁影响(用作参考的OSU实验报告［12］的数据已经过洞壁干扰修正)，计算域外边界为开口大气。计算域内网格拓扑结构为C型，翼型表面共430个网格点，首层网格高度保证y+≤1。计算域入口边界距离翼型前缘10倍弦长，出口边界距离翼型后缘15倍弦长，上下边界距离翼型10倍弦长。入口边界指定来流的速度矢量和湍流强度，出口边界指定出口静压。

由于工艺原因，风洞实验模型的尾缘不可能是完全尖的，而是稍有钝化的。因此，本文用于模拟的S809翼型尾缘也按照OSU实验报告描述的方法进行了微小尺度的钝化。

Phase VI叶片数值模拟的计算域如图1所示，分为叶片附近的旋转坐标系区域和外层的地面坐标系区域，也没有考虑风洞洞壁的影响。计算域的外边界分为迎风的入口边界、背风的出口边界和平行于来流的开口边界。入口边界指定来流的速度矢量和湍流强度，出口边界指定出口静压，开口边界允许气流流入或流出。应用了旋转周期边界，以利用风轮的轴对称特点节省网格数。对于两叶片的Phase VI风轮，周期边界的应用使得计算域只有实际区域的一半，从而使网格数和计算量减半。计算域入口边界与叶片的距离为8L(L为叶片展长)，出口边界距离叶片15L，开口边界距离叶片10L。叶片表面首层网格高度保证y+≤1。对不同网格密度的模型进行了网格无关性测试，发现网格单元总数大于450万时达到网格无关性要求，因此选用网格单元总数为450万的模型用于数值模拟。

NREL的报告中没有描述Phase VI叶片在加工过程中尾缘的实际钝化情况。考虑到工艺上可达到的尾缘尖锐程度应当是相近的，本文用于模拟的Phase VI叶片模型的尾缘按照与S809翼型相同的方法进行了微小尺度的钝化。

图1 Phase VI叶片数值模拟的计算域Fig.1 Simulation domain for blade Phase VI

2 S809翼型数值模拟结果

考虑到Phase VI实验中叶片各剖面的工作雷诺数大部分处于1×106附近，为使S809翼型的数值模拟结果与Phase VI叶片的数值模拟结果互相具有参考性，本文对S809翼型数值模拟的雷诺数选择为1×106。

图2显示了雷诺数1×106下S809翼型的实验值、SST模型转捩模拟结果和全湍流模拟结果。升力系数的实验值在迎角0°～8°之间基本呈线性增长，8°～16°之间变得平缓，16°～20°急剧下降，20°以后再次上升。结合OSU实验报告中S809翼型在各迎角下的压力分布可以判断，以上四个阶段的流场基本特征分别为:附着流动状态、尾缘分离的发生和发展、尾缘分离向完全分离的急剧转变、完全分离流动状态。其中除附着流动状态以外的其它三个状态可以分别被称为轻失速状态、轻失速向深失速的转变、深失速状态。

图2 S809翼型的升力系数(Re=1×106)Fig.2 Lift coefficient of airfoil S809(Re=1×106)

全湍流模拟的升力系数在附着流动阶段与实验结果吻合，在轻失速阶段则明显高于实验值。这与其他学者采用包括SST模型在内的各种线性涡粘性模型开展数值模拟得到的结果是一致的。随着迎角的继续增大，全湍流模拟由轻失速向深失速的转变较实验结果显著延迟，当实验升力系数已因深失速而急剧下降时，全湍流模拟的升力系数仍处于高位，导致在迎角20°左右出现很大的升力系数计算误差(高达约80%)。进入深失速以后，全湍流模拟的升力系数低于实验值，这可能是两方面原因导致的:一方面，大迎角下翼型本身对风洞的阻塞度增大，同时深失速状态下翼型背风面的漩涡的尺度比翼型本身还大，洞壁干扰的影响大增导致实验数据的洞壁干扰修正不足;另一方面，深失速状态的实际流场为非定常，且分离漩涡内湍流的发展对流场结构有较大影响，当前的定常数值模拟只是实际流场的平均近似，湍流模型也不足以准确描述漩涡内的湍流发展。

相对全湍流模拟结果，转捩模拟的升力系数在附着流动阶段略微高一些。关于转捩在该阶段的影响，Wolfe等人［4］和 R.B.Langtry 等人［10］已有相关论述，本文不再展开分析。

当迎角增大至轻失速状态后，转捩模拟的升力系数逐渐低于全湍流模拟结果，但仍然明显高于实验值。这说明转捩对S809翼型轻失速状态的气动特性有一定影响，但并不是造成数值模拟升力系数偏高的主要原因。图3显示了迎角16°下翼型周围的流线。由图可见，转捩模拟和全湍流模拟的流场结构基本相同，只在流场细节上有两点差异:一是转捩模拟结果在前缘附近存在层流分离泡，而全湍流模拟结果不存在;二是转捩模拟的流动分离区域稍大于全湍流模拟结果。

图3 迎角16°下S809翼型周围流线(上:转捩，下:全湍流)Fig.3 Streamlines around airfoil S809 at AOA 16°(up:transitional，down:full turbulence)

迎角16°～22°是转捩模拟和全湍流模拟的升力系数差别最大的区域。转捩模拟的升力系数与实验值同步急剧下降，全湍流模拟的升力系数则在更大的迎角下才开始急剧下降。图4显示了迎角20°下翼型周围的流线。由图可见，转捩模拟和全湍流模拟的流场结构有明显不同:前者发生了从前缘到尾缘的完全分离，漩涡区域覆盖了翼型的整个背风面;后者在前缘附近仍然为附着流动，漩涡尺度较前者明显小。结合16°迎角下转捩模拟中观察到的前缘分离泡推测，转捩模拟比全湍流模拟更早进入深失速的原因可能是前缘层流分离转捩后因逆压梯度较大不能再附，与已经推进至接近前缘的尾缘分离涡一同促成了完全分离的发生。

图4 迎角20°下S809翼型周围流线(上:转捩，下:全湍流)Fig.4 Streamlines around airfoil S809 at AOA 20°(up:transitional，down:full turbulence)

流动分离是附面层不能克服逆压梯度而在其底层发生“返流”导致的，附面层内速度剖面的丰满程度决定了其可以抵御的逆压梯度的强弱，而速度剖面的形状又与附面层的发展历程有关。因此，翼型的尾缘分离是受从前缘驻点开始的附面层发展历程影响的。在失速阶段，转捩模拟与全湍流模拟获得的附面层发展历程中最显著的区别在于前缘分离泡的存在。为了进一步分析前缘分离泡的影响，作者开展了指定转捩点的数值模拟，将十分接近前缘的地方(即将发生层流分离处)设定为强制转捩点，以避免层流分离泡的形成。在失速范围内，强制转捩的数值模拟结果与全湍流模拟结果相近。这说明，前缘分离泡的存在与否对翼型的失速特性有着决定性的影响。因此，转捩模拟与全湍流模拟在S809翼型失速特性上的差异，与其说是转捩的影响，不如说是前缘层流分离泡的影响更确切。

迎角大于22°以后，转捩模拟与全湍流模拟获得的升力系数一致。在这样大的迎角下，紧邻翼型前缘的逆压梯度非常大，以至于无论层流边界层还是湍流边界层都无法抵御，在几乎相同的位置发生分离，因此翼型表现出来的气动特性也几乎相同。

3 Phase VI叶片数值模拟结果

风速7m/s～25m/s下Phase VI风轮扭矩的实验值［5］、转捩模拟和全湍流模拟的数值模拟结果如图5所示。在风速为7m/s和9m/s时，转捩模拟和全湍流模拟结果相近，与实验值基本吻合;在风速9m/s～20m/s之间，全湍流模拟结果明显高于实验值，转捩模拟结果则与实验值接近。当风速继续增大到20m/s以上，转捩模拟和全湍流模拟的结果一致。

图5 Phase VI风轮的扭矩Fig.5 Low speed shaft torque of Phase VI rotor

叶片各个剖面的当地迎角是随着来流风速的增大而增大的。叶片气动力随风速的变化与S809翼型气动力随迎角的变化必然存在一定的关联。图5中同时画出了S809翼型的升力系数。通过对比可见，转捩模拟与全湍流模拟结果的差异在Phase VI风轮的扭矩与S809翼型的升力系数之间呈现出很强的相似性:小风速(小迎角)时，转捩模拟与全湍流模拟的风轮扭矩(翼型升力系数)相近;随着风速(迎角)的增大，转捩模拟结果明显低于全湍流模拟结果;当风速(迎角)继续增大至某一值以后，转捩模拟与全湍流模拟的结果再次一致。

图6显示了Phase VI叶片在风速9m/s、13 m/s和20 m/s时吸力面的摩擦力线，这几个风速下的流场结构代表了叶片从轻失速到深失速状态的典型流态。

风速9m/s时，叶片吸力面存在局部尾缘分离，越靠近叶根分离越严重，叶尖附近则为附着流动状态(说明叶片剖面迎角由叶根至叶尖逐渐减小)。转捩模拟和全湍流模拟的分离区域大小相近，但转捩模拟描述了层流分离泡的存在。这与S809翼型模拟中轻失速阶段的情况是相符的。

图6 Phase VI叶片吸力面的表面摩擦力线Fig.6 Friction lines on suction surface of blade Phase VI

风速13m/s时，转捩模拟的叶片吸力面除叶尖局部外的大部分区域已处于完全分离流动状态(深失速);全湍流模拟结果则在叶尖附近存在大一些的尾缘分离流动(轻失速)区域，完全分离(深失速)的区域相对转捩模拟结果小一些。这与S809翼型模拟中转捩模拟更早进入深失速的情况是相符的，也是转捩模拟获得的风轮扭矩在中等风速下低于全湍流模拟结果并与实验值更接近的主要原因。

风速20m/s时，转捩模拟和全湍流模拟的叶片吸力面均已处于完全分离流动状态(深失速)，两者获得的流场细节未见明显差别。这与S809翼型模拟中转捩模拟和全湍流模拟结果在深失速以后变得一致的情况是相符的。

以上关于各风速下流场形态的分析进一步将叶片模拟结果与翼型模拟结果联系起来。说明转捩对叶片失速特性和翼型失速特性产生影响的方式是相似的，都是通过前缘层流分离泡的作用体现出来。不过需要说明的是，本文转捩模拟获得的Phase VI风轮扭矩与实验数据符合较好，除了考虑转捩使叶片进入深失速的时机与实验更接近外，还因为轻失速剖面的正误差和深失速剖面的负误差相互有所抵消，因此不能认为转捩的引入完全解决了Phase VI叶片气动力预测不准确的难题。

4 结论

本文采用Gamma-Theta转捩模型对S809翼型和NREL Phase VI叶片进行了气动力数值模拟，分析了转捩对该翼型和叶片失速特性的影响。S809翼型的数值模拟结果表明:在轻失速阶段，转捩模拟得到的升力系数稍低于全湍流模拟结果，仍明显高于实验值;在轻失速向深失速的转变阶段，转捩模拟的升力系数与实验值同步急剧下降，而全湍流模拟结果进入深失速的时机明显滞后;进入深失速以后，转捩模拟与全湍流模拟获得的升力系数一致。结合典型迎角下的流场细节，发现前缘层流分离泡的存在影响了翼型的失速特性，导致深失速更早发生。

Phase VI叶片的数值模拟与S809翼型的数值模拟结果之间存在较显著的关联性。小风速时转捩模拟与全湍流模拟结果相近;中等风速时转捩模拟得到的风轮扭矩显著低于全湍流模拟结果，更接近实验值;大风速时转捩模拟与全湍流模拟的结果一致。这与翼型数值模拟中随着迎角的增大，转捩模拟与全湍流模拟得到的升力系数的变化规律是相符的。结合对典型风速下叶片表面流动状态的分析，认为转捩对叶片失速特性和翼型失速特性产生影响的方式是相似的，都是通过前缘层流分离泡的作用体现出来。

综上所述，本文认为与层流和转捩相关的前缘分离泡的存在会导致S809翼型和Phase VI叶片失速特性的变化，特别是会导致深失速的更早发生。考虑到前缘分离泡对边界层流动的作用机制并不因翼型或叶片的不同而不同，因此认为以上结论对其它翼型和叶片同样适用。但另一方面，由于S809翼型和Phase VI叶片流动显示实验结果的缺乏，导致本文数值模拟中捕捉到的前缘层流分离泡暂不能得到实验验证。不过这并不影响本文关于前缘层流分离泡对翼型和叶片失速特性影响的结论在性质上的正确性。

［1］DAN M，SOMERS.Design and experimental results for the S809 airfoil［R］.Colorado:National Renewable Energy Laboratory，1997.

［2］HAND M M，SIMMS D A，FINGERSH L J，etal.Unsteady aerodynamics experiment phase VI:wind tunnel test configurations and available data campaigns［R］.Colorado:National Renewable Energy Laboratory，2001.

［3］YANG S L，CHANG Y L，ARICI O.Incompressible Navier-Stokes computation of the NREL airfoils using a symmetric total variational diminishing scheme［J］.Journal of solar energy engineering，1994，116:174-182.

［4］WALTER P W，STUART S O.CFD calculations of S809 aerodynamic characteristics［R］.AIAA 97-0973，1997.

［5］SIMMS D，SCHRECK S，HAND M，et al.NREL unsteady aerodynamics experiment in the NASA-Ames wind tunnel:a comparison of predictions to measurements［R］.Colorado:National Renewable Energy Laboratory，2001.

［6］SØRENSEN N N，MICHELSEN J A，SCHRECK S.Navier-Stokes predictions of the NREL Phase VI rotor in the NASA Ames 80-by-120 wind tunnel［J］.Wind Energy，2002，5:151-169.

［7］BENJANIRAT S，SANKAR L N.Evaluation of turbulence models for the prediction of wind turbine aerodynamics［R］.AIAA 2003-0517，2003.

［8］EARL P N.DUQUE，BURKLUND M D.Navier-Stokes and comprehensive analysis performance predictions of the NREL Phase VI experiment［R］.AIAA 2003-0355，2003.

［9］SCHMITZ S，CHATTOT Jean-Jacques.Wind turbine blade aerodynamics of the NREL phase VI rotor near peak power［R］.AIAA 2005-4850，2005.

［10］LANGTRY R B，GOLA J，MENTER F R.Predicting 2D airfoil and 3D wind turbine rotor performance using a transition model for general CFD codes［R］.AIAA-2006-395，2006.

［11］MENTER F R，LANGTRY R，VOLKER S.Transition modelling for general purpose CFD codes［J］.Flow Turbulence Combust，2006，77:277-303.

［12］RAMSAY R R，HOFFMANN M J，GREGOREK G M，et al.Effects of grit roughness and pitch oscillations on the S809 airfoil［R］.Colorado:National Renewable Energy Laboratory，1995.