0 引言

随着风电场大规模发展,风电场经济性和风力机能效改善逐步受到重视^[1]。涡流发生器(vortex generators,VGs)技术是一种基于控制边界层分离的增升减阻技术^[2]。VGs技术最早应用于飞机上,因不需要改造基本构型,使用简易方便,控制效果显著,经不断更新发展,成为现有飞机上应用最为广泛的湍流边界层分离的被动控制技术^[3]。后VGs技术也被用于风力机叶片表面的流动控制。随着风力发电机组容量增大,叶片变得越来越长。为确保在恶劣风况环境下风力机叶片能够安全运行,在叶片根部部位采用厚翼型设计,甚至为圆柱形。厚翼型在大攻角条件下容易发生气流分离的现象,使阻力增加,升力降低,影响叶片的气动性能^[4]。在这种情况下,采用VGs技术能使分离的气流重新附着在叶片表面,达到改善气动性能的效果^[5]。

VGs气动外形须根据流动分离形态进行设计,需要对控制的流动现象进行准确掌握,也需要对涡流发生器的形状、位置、角度及排列方式等结构参数进行优化,全面评估对叶片气动影响^[6]。以风洞试验为主开展上述研究需要耗费大量人力物力,周期长费用高。在先期开展相关数值模拟研究,给出流动控制方案,再进行风洞试验验证,是比较理想的方法。20世纪美国国家可再生能源实验室对安装有VGs的风力机进行了大量性能研究^[7]。Lin利用微型VGs抑制边界层分离,通过风洞试验和数值模拟对VGs叶片的流体动力学进行了详细描述^[8]。Johansen分别对安装了VGs的非旋转叶片段、旋转叶片段和非旋转风力机叶片进行了数值模拟,但是由于计算网格数量等因素的影响,没有得到理想的结果^[9]。宋娟娟采用CFD方法对在上表面增设环形扰流器的风力机叶片进行了数值模拟,结果表明扰流器有效减小叶片上表面的流动分离^[10]。Yang Ke采用CFD方法研究了涡流发生器对风力机专用翼型DU91-W2-250翼型和改进后的DU91-W2-250-6翼型气动特性的影响,研究表明在大攻角下,涡流发生器具有推迟流动分离,减小翼型阻力,增大翼型升力系数的作用^[11]。然而,现有数值计算研究基本采用湍流模型来模拟分析涡流发生器叶片的气动特性。

大型风力机叶片根部的雷诺数小于3×10⁶,这意味着转捩现象必然存在。对于光滑叶片而言,转捩现象对叶片气动特性具有重要影响。Yang于1994年对S809翼型采用SST模型进行了模拟研究,在研究结果指出：SST模型计算得到的升力系数与试验结果吻合良好,但阻力系数远大于实验值^[12]。Walter于1997年在指定了转捩点的方式下对S809翼型进行模拟研究,结果显示阻力系数误差明显减少^[13]。近几年不少学者考虑了转捩的影响,对翼型进行气动特性研究,结果都表明转捩效应的影响在叶片气动研究中不可忽略^[14-15]。

基于以上分析,在考虑转捩效应的前提下,采用γ-Re_θt转捩模型,对安装有三角翼的VGs的DU91-W2-250叶片段进行气动特性分析。

1 计算模型

1.1 转捩输运方程

γ-Re_θt 转捩模型是由langtry和Menter在2004年提出的一种耦合到SST k-ω湍流方程^[16-17]中的,完全基于流场当地变量的转捩模型^[18]。

γ-Re_θt 转捩模型考虑了自然转捩、跨越转捩和分离流转捩等多种转捩机制,可进行准确预测^[18-20]。转捩预测包括预测转捩起始位置和预测转捩区,模型求解当地转捩雷诺数Re_θt和间歇因子γ这2个变量的标准输运方程。其中,Re_θt构成预测转捩起始位置的判据,γ用来模拟转捩区域的流动。

间歇因子γ代表了在转捩区某一固定点流动处于有脉动和无脉动状态的时间比例,其输运方程为

$\frac{\partial (\rho \gamma )}{\partial t}+\frac{\partial (\rho {{U}_{j}}\gamma )}{\partial {{x}_{j}}}=\frac{\partial }{\partial {{x}_{j}}}[(\mu +\frac{{{\mu }_{t}}}{{{\sigma }_{\gamma }}})\frac{\partial \gamma }{\partial {{x}_{j}}}]+{{P}_{\gamma }}-{{E}_{\gamma }}$ (1)

${{P}_{\gamma }}={{F}_{\text{length}}}{{c}_{a1}}\rho S{{[\gamma {{F}_{\text{onset}}}]}^{0.5}}(1-\gamma )$ (2)

${{E}_{\gamma }}={{c}_{a2}}\rho \Omega \gamma {{F}_{\text{turb}}}({{c}_{e2}}\gamma -1)$ (3)

式中：μ为动力黏度系数,Pa.s;S为应变率;F_length为转捩长度函数;F_onset、F_turb为转捩控制函数;Ω为涡量;c_a1、c_a2、c_e2、σ_γ为转捩常数。

动量厚度雷诺数\(\tilde{R}{{e}_{\theta t}}\)的输运方程为

\(\frac{\partial (\rho \tilde{R}{{e}_{\theta t}})}{\partial t}+\frac{\partial (\rho {{U}_{j}}\tilde{R}{{e}_{\theta t}})}{\partial {{x}_{j}}}=\frac{\partial }{\partial {{x}_{j}}}[{{\sigma }_{\theta t}}(\mu +{{\mu }_{t}})\frac{\partial \tilde{R}{{e}_{\theta t}}}{\partial {{x}_{j}}}]+{{P}_{\theta t}}\)(4)

${{P}_{\theta t}}={{c}_{\theta t}}\frac{\rho }{t}(R{{e}_{\theta t}}-\tilde{R}{{e}_{\theta t}})(1.0-{{F}_{\theta t}})$ (5)

$t=\frac{500\mu }{\rho {{U}^{2}}}$ (6)

式中：t为时间尺度;混合函数F_θt为边界层中关闭源项;c_θt、μ_t、σ_θt为转捩常数。

γ-Re_θt 转捩模型详细的描述,以及各参数的计算方法可参阅文献[14]。自由来流边界上γ 取1,

\(\tilde{R}{{e}_{\theta t}}\)依据来流湍流强度进行计算;固体边界上γ、\(\tilde{R}{{e}_{\theta t}}\)均为0。

1.2 修正的SST k-模型

间歇因子的大小代表流动脉动的强度,而湍动能方程控制流动脉动的生产和耗散,通过间歇因子对湍流动能方程产生项P_k和耗散项D_k的修正,实现了与SST k-ω湍流模型的耦合。修正后的生产和耗散项源项分别为式(7)和(8)^[15,20]。

\({{\bar{P}}_{k}}={{\gamma }_{\text{eff}}}{{P}_{k}}\) (7)

\({{\bar{D}}_{k}}=\min (max({{\gamma }_{\text{eff}}},0.1),1){{D}_{k}}\) (8)

式中下角标eff代表有效参数。

2 计算对象

荷兰Delft科技大学在低速风洞中,对装有VGs的采用风力机专用翼型DU91-W2-250的叶片段,进行了试验研究^[20],并有具体试验数据可以查询,这为很多学者数值模拟研究的验证提供了便利。DU91-W2-250翼型是由Delft科技大学设计的风力机专用翼型,其最大厚度为25%弦长,位于距前缘点32.2%弦长处,结构如图1所示。该试验的叶片段展长为0.175m,弦长为0.6m,在叶片弦长20%处,共安装了5组VGs。VGs几何尺寸如图2^[21]所示,VGs的叶片采用三角翼结构,在弦向有效长度为17mm、高5mm,三角翼的攻角为16.4°,1组内2个三角翼叶尖间距为10mm,相邻2组三角翼叶尖间距为25mm。为了便于与试验数据^[20]进行对比,本文对该试验模型进行了计算。

图1 几何模型 Fig. 1 Figuration of aerofoil and VGs

图2 VGs结构尺寸 Fig. 2 Parameters of structure of VGs

采用ANSYS ICEM生成结构化网格,总网格数为300万。对VGs处进行了网格加密处理,第1层网格高度为1×10^-5m,网格增长率为1.1,保证y+≤1.5。计算雷诺数Re=2.0×10⁶,攻角范围为0°到20°。定义叶片微段两侧断面为周期性边界条件。

3 模型验证

分别采用SST k-ω全湍流模型和γ-Re_θt 转捩模型对DU91-W2-250叶片段进行计算分析。图3为2种模型计算得到的升力系数曲线和Delft科技大学试验结果^[20]对比图。从图3可以看出,转捩模型计算结果与试验结果更加吻合。如在0°≤α≤18°攻角范围,转捩模型计算误差仅为5.6%;在α>18°攻角范围,转捩模型计算误差仅为7.9%。SST全湍流模型计算结果在0°≤α≤18°范围内与试验数据吻合良好,误差也较小;在α>18°后,升力系数曲线逐渐偏离试验数据,整体低于试验升力系数值,在α=20°时最大误差达23.1%。

图4为2模型得到的VGs叶片升阻比曲线与试验结果对比图。比较SST模型和γ-Re_θt 转捩模型计算结果可以看出,SST模型的升阻比曲线与试验结果存在较大误差,最大误差达51.8%。相对而言,γ-Re_θt转捩模型计算结果与试验数据吻合良好。

图3 升力系数曲线 Fig. 3 Curves of lift coefficient

图4 升阻比曲线 Fig. 4 Curves of ratio of lift to drag

从以上升力系数和升阻比曲线证明,转捩模型计算结果精度更高,其模拟结果可以更真实地反映VGs叶片的气动特性。

4 气动性能分析

结合图1和图2可以看出,在VGs叶片段上有2个典型位置,分别位于每组VGs中心线和相邻2组VGs的中心线处。下文主要对VGs中心线位置的气动特性参数计算结果进行了追踪和对比分析。

4.1 转捩点变化规律

定义摩阻系数C_f为

${{C}_{f}}=\frac{2{{\tau }_{x}}}{\rho {{v}^{2}}}$ (9)

式中：τ_x为沿流动方向的剪切应力,Pa;ρ为空气密度,kg/m³;v为来流风速,m/s。

图5是攻角为5°、12°、20°工况下叶片表面摩阻系数曲线。摩阻系数与叶片表面切向剪切应力成正比,在计算结果中转捩位置通常取壁面摩阻曲线层流段最低点^[14] ,从图5可以看出,转捩模型准确地捕捉了转捩现象。

在α=5°时,转捩模型摩阻系数的计算值明显小于SST模型的计算值,这是SST模型升阻比系数计算值误差较大的主要原因之一。进一步分析转捩模型计算结果可看出,上表面气流在0≤x/c≤0.2处摩阻系数剧烈波动,表明该区域流动参数变化剧烈。在x/c=0.2附近,摩阻系数骤升和骤降,这是由VGs引起的局部摩阻力突然大幅波动的结果。在VGs下游表面摩阻系数逐渐下降,并在x/c=0.53处摩阻系数产生最低点,表明叶片表面发生转捩。在转捩点前,转捩模型的摩阻系数远小于SST模型计算值,这是因为翼型表面流态基本为层流;在转捩点后x/c=0.6处,转捩模型和SST模型计算结果基本一致,这是因为翼型表面流体转捩完成,发展成为湍流。湍流模型计算结果显示,摩阻系数从叶片前缘点迅速增长为正值,该模拟结果表现为气流接触下表面后立即发生转捩,湍流气流附着在表面上。由此得出SST模型不能准确捕捉带VGs叶片的转捩现象。

图5 攻角为5°、12°、20°的摩阻系数曲线 Fig. 5 Friction coefficient curve of 5°, 12°, 20° attack angle

当攻角α=12°时,从转捩模型计算结果可以看出,上表面摩阻系数先出现负值随后上升为正值,表明上表面先在前缘出现分离泡,随后出现再层流化现象,因此导致摩阻系数上升。在VGs的影响下x/c=0.2摩阻系数骤升。在x/c=0.22附近,摩阻系数达到最低点,表明气流产生了由层流转变为湍流的转捩现象。SST模型计算值显示,在x/c=0.22转捩点之前,上表面摩阻系数整体比转捩模型的要大;在转捩完成x/c=0.32附近,摩阻系数基本相近。下表面摩阻系数2模型计算基本相近,都表现为下表面基本为层流状态。

当攻角α=20°时,转捩模型结果显示上表面摩阻系数首先出现负值,然后迅速增大较大正值,表明此时叶片表面在前缘点产生了分离泡,并在附近很快发生转捩,并迅速完成了转捩过程。在x/c=0.52处摩阻系数为负,出现气流分离。SST模型计算结果则显示,上表面在x/c=0.2之后就快速进入分离状态。x/c=0.52之后2模型流态基本一致。

从以上3个攻角下的摩阻系数分析可以得出,由于SST模型不能准确捕获转捩现象,因此其模拟的摩阻系数及其反映的叶片表面流体流态,与真实值均有较大差别。

4.2 表面速度型规律

结合图3和图4可以看出,VGs叶片段在5°攻角左右升阻比与试验值误差最大,为了进一步分析转捩给带VGs叶片气动特性造成的影响。图6给出了α=5°时,叶片表面不同位置的方向速度线型。图6中h为VGs高度,纵坐标y/h表示沿壁面法向的高度相对值,横坐标u/U为垂直于壁面的切向速度与自由来流速度的比值。

从x/c=0.1的速度型图可看出,在越靠近y=0的部位,u/U越小接近于零,但是在离开壁面后,u/U很快上升为1,这是因为此工况下该部位流体的流态为层流,其边界层厚度约为0.2h高度。此时SST模型和转捩模型模拟结果比较一致。

图6 5°攻角下速度型图 Fig. 6 Structure of velocity in 5° attack angle

从图6(b)可以看出,在VGs之后,叶片表面流体速度整体较x/c=0.1的要低,表明流体处于逆压流动而减速。在y/h=1高度的近壁面的速度产生最大速度,这是因为流体流过VGs时,在VGs顶端卷起叶尖涡,使这一高度流体的流动发生了变化,具有更高的能量,流速更高。但在距离壁面更远的地方,受VGs影响小,速度低于近壁面的流体流速。

在x/c=0.5处,转捩模型计算结果显示,近壁面流速比远壁面流速更大,这一差别现象比x/c=0.3的更加明显,这是因为在该位置基本处于转捩点位置,转捩现象也使得流速具有更高的能量,因此使得近壁面流速保持高能量流动;而离壁面远的流体受VGs和转捩效应影响较弱,随主体流体逆压流动,速度进一步减小。比较湍流模型和转捩模型的计算结果可以看出,湍流模型速度要小于转捩模型计算结果,这是因为湍流模型不能准确捕获转捩现象的原因。

在x/c=0.7附近,转捩效应结束和VGs的影响也逐渐减弱,使得壁面法向速度变化率变化,远壁面区气流与未扰动区域的流速保持一致,且边界层厚度增大,约为1.2h。

4.3 表面压力变化规律

图7为不同攻角下,叶片表面的压力系数分布图。上、下翼面的压力之差反映了叶片升力大小。从图7可知,在x/c=0.2,气流受VGs扰动,在叶片表面形成的压力稍有波动,波动幅度随攻角增大而减小。在α=5°时,在VGs之前,转捩模型的上表面压力系数稍大于SST模型,在VGs后两者计算结果差别不大;下表面压力系数基本重合。在α=12°时,VGs之前上表面压力系数差距加大,总体2模型值相差较小。在α=20°时,2模型的上翼面压力系数差别较大,但下翼面的压力系数仍保持一致。上述压力系数分布规律解释了,在小攻角下湍流模型和转捩模型升力系数模拟值差距较小,和在大攻角压力系数差异较大的原因;同时也给出,在大攻角下上翼面压力系数差异较大引起了升力系数的差异。

图7 叶片表面的压力系数分布图 Fig. 7 Pressure coefficient of blade surface

5 结论

1）转捩模型计算升力系数曲线和升阻比曲线与试验值更加吻合,转捩模型具有较高的精度。纯湍流模型的升力系数误差最大达23.1%;升阻比系数最大误差达51.8%。

2）转捩模型可以准确的捕捉到带VGs叶片段的上下表面的转捩点,随着攻角的增大,转捩点由VGs下游逐渐向VGs上游过渡。

3）纯湍流模型计算的摩阻系数远高于转捩模型。这是纯湍流模型模拟升阻比系数偏低的主要原因之一。

4）在VGs和转捩效应的影响下,近壁面流体具有高能量,使流速高于远壁面区域流速。

5）在小攻角下,在VGs之前,转捩模型的叶片表面压力系数高于湍流模型的计算值;在大攻角下,上表面压力系数差异加大,导致大攻角下转捩模型的升力系数大于湍流模型计算值。

参考文献

[1] Xing H,Cheng H,Zhang L.Demand response based and wind farm integrated economic dispatch[J]. CSEE Journal of Power and Energy Systems,2015,1(4):47-51.

[2] 黄红波,陆芳.涡流发生器应用发展进展[J].武汉理工大学学报:交通科学与工程版,2011,35(3):611-614,618. Huang Hongbo,Lu Fang.Research progress of vortex generator application[J].Journal of Wuhan University of Technology:Transportation Science & Engineering,2011,35(3):611-614,618(in Chinese).

[3] Reuss R L,Hoffmann M J,Gregorek G M.Effects of Surface Roughness and vortex generators on the NACA 4415 airfoil[R].Golden,Colorado:National Renewable Energy Laboratory,1995.

[4] 张磊,杨科,徐建中.涡流发生器对风力机专用翼型气动特性的影响[J].工程热物理学报,2010,31(5):749-752. Zhang Lei,Yang Ke,Xu Jianzhong.Effects on wind turbine airfoils by vortex generators[J].Journal of Engineering Thermophysics,2010,31(5):749-752(in Chinese).

[5] Velte C M,Hansen M O L,Meyer K E,et al.Evaluation of the Performance of Vortex Generators on the DU 91-W2-250 Profile using Stereoscopic PIV[C]// Proceedings of the 14th International Conference on Information Systems Analysis and Synthesis.Orlando,United States:IEIC,2008.

[6] 郝礼书,乔志德,宋文萍.涡流发生器布局方式对翼型失速流动控制效果影响的实验研究[J].西北工业大学学报,2011,29(4):524-528. Hao Lishu,Qiao Zhide,Song Wenping.Experimentally studying effects of different layouts of vortex generator on controlling stall flow over airfoil[J].Journal of Northwestern Polytechnical University,2011,29(4):524-528(in Chinese).

[7] Griffin D A.Investigation of vortex generators for augmentation of wind turbine power performance [R].Golden,Colorado:National Renewable Energy Laboratory,1996.

[8] Lin J C.Review of research on low-profile vortex generators to control boundary-layer separation[J].Progress in Aerospace Sciences,2002,38(4-5):389-420.

[9] Johansen J,Sφrensen N N,Reck M,et al.Rotor blade computations with 3D vortex generators[R].Rosklide Denmark:Risø National Laboratory,2005.

[10] 宋娟娟,徐宇,陈洪胜,等.叶片扰流器对风力机性能影响的数值研究[J].工程热物理学报,2012,33(3):405-407. Song Juanjuan,Xu Yu,Chen Hongsheng,et al.Numerical simulation of blade spoiler effect on wind turbine performance[J].Journal of Engineering Thermophysics,2012,33(3):405-407(in Chinese).

[11] Yang Ke,Zhang Lei,Xu Jianzhong.Simulation of aerodynamic performance affected by vortex generators on blunt trailing-edge airfoils[J].Science in China Series E:Technological Sciences,2010,53(1):1-7.

[12] 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(4):174-182.

[13] Wolfe W P,Ochs S S.CFD calculations of S809 aerodynamic characteristics[R].Office of Scientific & Technical Information Technical Reports,AIAA 97-0973,1997.

[14] Menter F R,Langtry R B,Likki S R,et al.A correlation- based transition model using local variables-part I:model formulation[C]//Proceedings of 2004 ASME Turbo Expo.Vienna,Austria:ASME,2004.

[15] 陈奕,高正红.Gamma-Theta转捩模型在绕翼型流动问题中的应用[J].空气动力学学报,2009,27(4):411-418. Chen Yi,Gao Zhenghong.Application of gamma-theta transition model to flows around airfoils[J].Acta Aerodynamica Sinica,2009,27(4):411-418(in Chinese).

[16] 张骏,袁奇,吴聪,等.大型风力机叶片表面粗糙度效应数值研究[J].中国电机工程学报,2014,34(20):3384-3391. Zhang Jun,Yuan Qi,Wu Cong,et al.Numerical simulation on the effect of surface roughness for large wind turbine blades[J].Proceedings of the CSEE,2014,34(20):3384-3391(in Chinese).

[17] 徐浩然,杨华,刘超,等.不同湍流模型对MEXICO风力机气动性能预测精度的研究[J].中国电机工程学报,2013,33(35):95-101. Xu Haoran,Yang Hua,Liu Chao,et al.Research on the effect of different turbulence models on the aerodynamic performance prediction accuracy of MEXICO wind turbines[J].Proceedings of the CSEE,2013,33(35):95-101(in Chinese).

[18] 侯银珠,宋文萍,张坤.考虑转捩影响的风力机翼型气动特性计算研究[J].空气动力学学报,2010,28(2):234-237. Hou Yinzhu,Song Wenping,Zhang Kun.Calculation of aerodynamic performance of wind turbine airfoil incorporating transition prediction[J].Acta Aerodynamica Sinica,2010,28(2):234-237(in Chinese).

[19] 钟伟,王同光.转捩对风力机翼型和叶片失速特性影响的数值模拟[J].空气动力学学报,2011,29(3):385-390. Zhong Wei,Wang Tongguang.Numerical analysis of transition effect on stall performance of wind turbine airfoils and blades[J].Acta Aerodynamica Sinica,2011,29(3):385-390(in Chinese).

[20] 王强,赵宁,王同光,等.考虑转捩的风力机翼型动态失速数值模拟[J].太阳能学报,2012,33(1):113-119. Wang Qiang,Zhao Ning,Wang Tongguang,et al.Numerical simulation of wind turbine airfoil dynamic stall with transition modeling[J].Acta Energiae Solaris Sinica,2012,33(1):113-119(in Chinese).

[21] Timmer W A,van Rooij R P J O M.Summary of the Delft university wind turbine dedicated airfoils[J].Journal of Solar Energy Engineering,2003,125(4):488-496.