1.2MW风力机整机流场的数值模拟

李少华, 匡青峰, 吴殿文, 郭婷婷, 王梅丽

(1.东北电力大学 能源与动力工程学院,吉林 132012;2.大唐山东新能源有限公司,青岛 266061;3.北京国电龙源环保工程有限公司,北京 100052;4.云南龙源风力发电有限公司,曲靖 655600)

风能是一种洁净的可再生资源,其储能丰富,据有关研究资料表明:近地层的风能总量约为13×1012kW,如果其中的1%被有效利用,就可以满足人类对能源的需求[1-2].近年来,我国的风力发电产业持续快速增长.风力机的设计、气动性能以及质量是风力机组稳定运行的决定性因素.对于中、小型风力机,其叶片的气动性能可以通过风洞试验得到,但是随着风力机功率和塔架高度的增加,大型风力机叶片的风洞试验难以实现,因此只能采用CFD软件来计算.该计算软件具有简洁、周期短、通用性强及费用低等优势.

风力机整机的流场特性分析是研究风力机流场内空气流动的有效方法.在风力机叶轮气动性能的CFD软件计算中,N.N.Sorensen等[3]采用 EllipSys3D软件对叶片及风轮采用SST k-ε湍流模型在不同来流风速进行了稳态和非稳态流场计算研究;J.Johansen等[4]对NREL Phase V I叶片在非旋转情况下采用EllipSys 3D软件基于SST k-ω湍流模型进行了大涡模拟;E.P.N.Duque等[5]对由PhaseⅡ叶片组成的风力机使用Overflow软件进行了非稳态模拟,使用可压缩的N-S方程,湍流黏性采用Baldwin-Lomax模型;L.Makoto等[6]采用ATP软件求解可压缩N-S方程,利用多重网格系统和标准k-ε模型及Baldwin-Lomax模型,模拟风力机三维流场;杨瑞等[7]对NREL水平轴风力机在风速为 13m/s工况下,采用SST k-ε湍流模型以及分离DES 3种湍流模型进行了三维旋转流场数值模拟,计算结果与NASA非稳态气动力学风洞试验结果进行了对比和分析,得出SST k-ω模型能够很好地描述流场流动,比分离DES湍流模型的计算结果更接近试验结果;张玉良等[8-9]使用k-ε湍流模型,对整机进行了在旋转工况下的模拟,并对模拟的结果进行了尾涡分析;吴殿文[10]使用A nsys软件对风力机整机的受力及流场进行了分析,采用SST k-ε湍流模型对流场进行模拟,得出在风力机叶片振动中,挥舞振动和摆振振动的耦合振动将起到主要作用.

随着风力机功率的增大,叶片尾迹对塔架的相互作用不仅改变了塔架的绕流,而且对风机整机的输出功率产生重要影响.笔者对1.2 MW风力机整机流场进行了研究,并分析了叶片转动与塔架之间的相互干扰对风机输出功率的影响,从而为多机阵列的研究打下基础,同时对风力机系统的非定常特性以及气动弹性稳定性和噪声辐射的研究均有较重要的意义.

1 几何建模与网格划分

1.1 叶片设计

本文模拟用叶片在威尔森理论的基础上,考虑了风机叶尖损失及轮毂损失的工况下,采用Matlab编程设计求得,叶片设计的基本参数见表1.从表1可知:风机的叶轮直径为70 m,轮毂直径取2m,叶片长度取34 m;预弯发生在叶片距离轮毂1/3处背离塔架方向;翼型取NACA 44××系列.笔者采用Solidworks软件直接根据翼型的空间坐标生成翼型曲线放样成叶片几何模型(图1).

表1 叶片设计的基本参数Tab.1 Basic design parameters of b lade

图1 叶片的几何模型Fig.1 Geometric model of blade

1.2 整机建模

将叶片导入Gambit软件中,列出另外2个叶片,与原来叶片组成3叶片风轮.在画出导流罩后,将3叶片与导流罩进行布尔运算,合并为风轮实体,然后再将风轮与机舱、塔架合并,组成整机实体.机舱下底面半径为1.6 m,上底面半径为1.4m,圆台高为9 m,其中下底面与轮毂连接,塔架的高度选为80m,靠近地面的下底面半径为2m,与机舱连接部分的半径为1m.

1.3 流场的几何尺寸

风轮的转速为18.44 r/min,来流风速为11.26 m/s,入口距离风轮105 m,出口距离风轮350 m.图2为整机的流场区域和尺寸.流体旋转小区域前面距离风轮为8m,后面距离风轮为3 m,塔架距离风轮为5.2 m.

图2 整机的流场区域和尺寸Fig.2 Flow-field region and size of wind turbine

1.4 网格划分

整机由风轮、机舱和塔架组成.叶片表面扭曲复杂,需要在风力机周围区域进行网格加密处理.为了获得较高的网格质量,将整个流场区域划分为几块:旋转小区域包含风轮、轮毂及半截机舱;剩余机舱及塔架另行建立小区域,分别对这2个小区域及边界近壁面采用size function函数进行加密处理,并使用Tgrid对该区域进行网格划分,其余流场分块采用六面体网格.在网格中导入Fluent软件后,进行网格检查和尺寸转换.图3为流场区域的网格划分.从图 3可知:包含风轮的旋转小区域网格数为1 732 448,包含半个机舱及塔架的小区域网格数为976 714,其余区域为六面体网格,网格数为858 144,总网格数约为356万.为验证网格无关性,笔者在进行计算时,比较了3种网格数312万、356万和391万的不同模拟结果,并分析了其不同截面上湍动能的变化;以穿过机舱中心线上的湍动能变化为参考,发现粗网格的计算结果与2种较细网格的计算结果差别较大,而后2种较细网格的计算结果吻合较好.因此,为了节省计算时间,笔者选用

图3 流场区域的网格划分Fig.3 G rid division of the flow field

2 数学模型与计算方法

2.1 数学模型

在本文中,假设叶片为刚体,且模拟过程不考虑叶片表面的变形.笔者基于稳态不可压缩流动三维定常雷诺时均N-S方程(RANS)进行了数值模拟,采用Segregated隐式求解器三维稳态算法,紊流模型使用SST k-ω模型,压力-速度耦合采用Simp le算法,对流项差分格式采用二阶迎风格式[10].

控制方程通用形式:

式中:φ为通用变量;Γ为广义扩散系数;S为广义源项,为张力的常量.

SST k-ω紊流模型方程:

k方程

式中:F1、F2分别为混合函数,F2用来修正F1在自由剪切流中的误差.

这就是以k-ω模型为基础的SST(Shear Stress Transport)模型,式中参数见文献[11].

2.2 数值计算方法与边界条件

采用SST k-ω湍流模型,通用控制方程的离散采用有限容积法,对流项差分格式采用二阶迎风格式,壁面无滑移,其边界条件如下.

速度入口:u=u∞=11.26m/s,v=w=0 m/s;

速度出口:选为自由出流,设定为outflow;

叶片及轮毂:选w all边界,以 y轴为旋转轴绝对旋转,转速 ω=18.44 r/m in,壁面无滑移;

机舱和塔架:选w all为边界,绝对静止,壁面无滑移;

包含叶片的流场区域:以y轴为旋转轴旋转,转速 ω=18.44 r/min;

地面:绝对静止,绝热无滑移.

3 结果与分析

整机的模拟收敛后,对流场中某些特定的区域,建立切片并进行性能的定性分析.图4为流场区域的切片截面.在图4中,在y方向建立4个截面:y=-2m、y=2 m、y=8 m和y=20 m;x方向建立4个截面:x=-33 m、x=-20m、x=20m和 x=33 m;z方向建立2个面:z=-33m和z=33m,然后分析这些截面上的结果.

图4 流场区域的切片截面(单位:m)Fig.4 Sliced cross sections of the flow field(unit:m)

3.1 风力机的输出功率

表2为整机和轮毂在各坐标轴的转矩.从表2可知:3个叶片和轮毂对 y轴的总转矩为577 438.920 N◦m,经计算得到整机的输出功率为1 114 487.91W.表3为叶片和轮毂在各坐标轴的转矩.通过表2和表3中的转矩可得到输出功率,风轮模拟的输出功率为1 170 999.8W,而整机的输出功率为1 114 487.91W,从而可以得到输出功率的相对误差,其中风轮的相对误差为2.42%,整机的相对误差为7.13%.

从表2和表3可以看出:在模拟条件均一致的情况下,整机比风轮多出机舱和塔架部分,但其输出功率却比风轮减少了56 511.89W.可见,由于叶片与塔架间的相互干扰以及塔架对流场的影响等原因,导致整机的输出功率减少.

表2 整机和轮毂在各坐标轴的转矩Tab.2 Torques obtained through turbine and hub N◦m

表3 叶片和轮毂在各坐标轴的转矩Tab.3 Torques obtained through wheel and hub N◦m

3.2 y方向不同截面上的压力

图5为y方向不同截面上的静压力.

从图5可以看到气流通过整机时气流压力的变化:图5(a)为风机前方2 m截面上的静压力,叶尖处压力最大,并向轮毂逐渐减小,塔架的存在使得塔架前方因速度减小而出现压力增大区;图5(b)为风力机后2 m的静压力,此截面位于风轮之后塔架之前,由于风轮的作用使得风轮后方形成类似圆形的低压区,叶片对后方压力分布的影响减小,而由于塔架的存在,使得叶尖部位压力比其余区域大很多;从图5(c)可以看到:叶片对压力的影响较小,因此可以忽略叶片的作用,看作仅在风轮的作用下形成了类似圆形的低压区,由于塔架的存在,气流绕流后在塔架后方形成高压区;从图5(d)可以看到:随着距离增大,风轮和塔架对流场压力的影响进一步减弱,但塔架对压力的影响范围相对增大.

由此可见,塔架的存在使得风机后方的压力场分布发生了较大变化,而且叶片和塔架间的干涉也对压力的分布产生了影响.但从总体看,塔架对风力机后方远场的压力分布影响有限.

3.3 不同截面的湍动能

图6为不同截面上的湍动能.其中,图6(a)为以旋转轴为轴线、半径为20 m圆柱侧面上的湍动能,可以看到叶片的旋转效应非常明显.在该圆柱面上,气流经过风轮后湍动能在叶片周围小区域内;由于附着涡的存在,使得叶片周围的湍动能较大,在叶片后方,随着附着涡及叶尖旋转尾涡的脱落,在风轮后方形成了湍动能较大区域,在两叶片之间,前一叶片的影响较小,因而湍动能较小.

图5 y方向不同截面上的静压力Fig.5 Static pressure distribution on various sections in y direction

在图6(b)中,由于整机中机舱和塔架的存在,使得气流经过风力机后发生了复杂的变化,在风轮直径范围内,湍动能变化剧烈,中心区域的湍动能未必小于外侧的湍动能.由此可见,由于机舱和塔架的存在使得风力机后方的湍动能发生明显变化,且在风力机后方区域内,风力机对流场的影响复杂.

图6 不同截面上的湍动能Fig.6 Turbulent kinetic energy on different sections

在图6(a)和图6(b)的截面上,从迎风面看位于风力机的左部,而图6(c)和图6(d)的截面为风轮的右部,风轮为顺时针旋转,由此可以看到:风轮下半部分湍动能较大,主要是因为此处存在叶片的缘故,而风轮上半部分湍动能较小,是因为此区域位于2叶片之间.前一叶片旋转的影响消除,使得风轮上半部分的后方湍动能较大,下一叶片的旋转还未对此区域产生影响,因此出现湍动能较小区域.而图6(c)和图6(d)则相反,在叶片的影响下,风轮下半部分湍动能较大,而上半部分湍动能较小.图6(d)的深色区域主要是由于叶片旋转影响及附着涡和螺旋尾涡的脱落造成的.

由于风力机整机包括机舱和塔架,对湍动能有较大的影响,直接改变了风轮后方的湍动能分布.因此,在研究多机阵列之前,对整机进行较为详细的模拟研究十分必要.

4 结 论

(1)叶片旋转过程中与塔架相互干扰,减小了风力机通过叶片获得的转矩,导致风力机的输出功率降低.

(2)塔架对风力机风场产生重要影响,使得塔架后方的湍动能变大,但叶片中心涡、附着涡对风力机流场区域的影响较小.

(3)风轮和整机对风力机后方几百米区域的流场产生影响,但流场受整机影响的区域小于受风轮影响的区域.

[1] 张义华,李隆键.水平轴风力机空气动力学数值模拟[D].重庆大学动力工程学院,2007.

[2] KENISARIN M K,VEDATM A A,MEHM ET A A.Wind power engineering in the world and perspectives of its development in Turkey[J].Renewable and Sustainable Energy Reviews,2006,10(4):341-369.

[3] SORENSEN N N,M ICHELSEN J A.Aerodynamic predictions for the unsteady aerodynamics experiment Phase II rotor at the national renew able energy laboratory[J].AIAA,2000(37):161-165.

[4] JOHANSEN J,SORENSEN N N,M ICHELSEN J A,et al.Detached-eddy simulation o f flow around the NREL phase V I:blade[J].wind Energy,2002,5(2/3):185-197.

[5] DUQUE E P N,VANDAM C P,HUGHES S C.Navier-Stokes simulation o f the NREL combined experiment,phase II:rotor[J].AIAA,1999(27):143-148.

[6] MAKOTO L,CHUICHI A.Three dimensional Navier-Stokes flow-field computations through horizontal axis wind turbine blade[J].AIAA,2001(58):112-115.

[7] 杨瑞,李仁年,张士昂,等.水平轴风力机CFD计算湍流模型研究[J].甘肃科学学报,2008,28(4):90-93.YANG Rui,LI Rennian,ZHANG Shiang,et al.A study on the turbulence models for the CFD calculation of the horizontal axis wind turbine[J].Journal of Gansu Sciences,2008,28(4):90-93.

[8] 张玉良.水平轴大功率高速风力机风轮空气动力学计算[D].兰州理工大学能源与动力工程学院,2006.

[9] 张玉良,李仁年,杨从新.水平轴风力机的设计与流场特性数值预测[J].兰州理工大学学报,2007,33(2):54-57.ZHANG Yu liang,LI Rennian,YANG Congxin.Design and numerical prediction o f flow field characteristics of horizontal-axis wind turbine[J].Journal of Lanzhou University of Technology,2007,33(2):54-57.

[10] 吴殿文.风力发电机叶片预弯设计及其数值研究[J].动力工程学报,2010,30(6):450-455.WU Dianwen.Pre-bend design and numerical simulation of wind turbo-generator[J].Journal of Chinese Society of Power Engineering,2010,30(6):450-455.

[11] 王福军.计算流体动力学分析——CFD软件原理与应用[M].北京:清华大学出版社,2004:125-129.