H型垂直轴风力发电机气动性能分析

单 乐,贺 龙,马中武,赵玲霞

(1.河西学院物理与机电工程学院,甘肃 张掖 734000; 2.河西学院新能源研究所,甘肃 张掖 734000)

目前,大型水平轴风力发电机仍占据着主体地位,但垂直轴风力发电机由于起动风速低、结构简单、无需偏航对风、维护方便等优点,已逐步转变为关注和研究的热点[1-2]。应用于垂直轴风力发电机空气动力学性能的研究方法主要有3种:理论分析法、实验探究法、数值计算法。目前,数值计算方法已成为一种重要研究方法,它通过计算机模拟在实验或理论计算中很难得到的复杂流场,其适用性强。

风轮是H型垂直轴风力发电机捕获能量的关键部件,风轮外流场非常复杂,空气对于转动部件是定常流动,但实际流动为惯性系下的非定常流动,传统的理论分析方法不能准确描述其外部流场。计算流体力学(CFD,computational fluid dynamics)中的移动网格是一种非定常计算模型,分析流动随时间的变化情况,采用移动网格模型对H型垂直轴风力机的风轮流场进行仿真具有可行性[3]。

二维CFD计算数据能够描述风轮外流场的一般特性,但忽略了气体在叶展方向的流动。三维CFD计算数据能够更加准确地描述实际风轮的空气动力学特性,在风沙环境下,利用离散相模型(DPM,discrete phase model)计算数据可反应固体颗粒对风轮外流场的影响。通过CFD技术,为改善H型垂直轴风力发电机组风轮的气动性能提供理论支持,并为风轮的优化设计提供依据。

1 数值计算模型

1.1 研究对象

以某小型3叶片H型垂直轴风力发电机为研究对象,风轮直径1 200 mm,叶片高度1 000 mm,叶片翼型为NACA0012,叶片弦长100 mm,叶片连接点位于1/4弦长处。在对垂直轴风力发电机空气动力学性能进行CFD数值模拟的过程中,将转动轴、支撑杆以及连接法兰等对风力机气动转矩获取影响较小的部件去除,以简化计算模型,降低计算成本,提高计算效率。

1.2 计算区域建模

为空气绕H型垂直轴风力机叶片流动的CFD数值计算域模型如图1所示,计算区域划分如图2所示。该流动属于外流问题,对风轮气动性能进行数值模拟时,静止区域主要用于模拟外部流场,保持静止,如图2(b)所示。为了使风轮产生的尾迹得到充分的发展,将入口边界和出口边界分别置于风轮上、下风向10倍和15倍风轮直径处。

图1 计算域参数及边界条件

图2 计算区域划分

旋转区域模拟风轮所在处的区域,是1个包含3个叶片的圆形区域,如图2(a)所示。叶片作为风轮捕获风能的主要载体,需要在叶片附件进行网格加密,生成高质量的计算网格,所以在旋转区域内部叶片所在位置处,又划分出3个相互独立的叶片子区域,如图2(c)所示。

1.3 网格划分

用CFD前处理工具Gambit进行计算域的网格划分,为保证数值计算的准确性,对于高雷诺数湍流模型,叶片表面近壁面区域采用壁面函数法或增强壁面函数法,垂直轴风力机叶片表面有旋涡不断分离,而后又重新附着,故叶片边界层网格采用增强壁面函数法处理[4],取y+=1,叶尖速比在0.5～3.0范围内,并根据平板湍流边界层理论,计算叶片表面第1层网格高度,计算公式为

(1)

其中:ρ为流体密度;L为特征尺寸,取翼型弦长值;μ为动力黏度;W为来流相对于叶片的相对速度;Re为雷诺数;Cf为壁面摩擦系数,取平板湍流边界层指数规律壁面摩擦系数;τw为壁面切应力;V∞为来流速度;V*为摩擦速度;y为叶片表面第1层网格高度。

相对速度W计算公式为

(2)

其中:θ为风轮方位角,即叶片位置;λ为叶尖速比,当风速一定时,其值越大则风轮转速越大。

当风速为10 m/s时,雷诺数Re约为2.67×105,计算可得近壁面处第1层网格的高度为0.008 mm。计算区域网格如图3所示,从叶片表面到叶片子区域与旋转区域交接面的节点间距增长率为1.1,生成二维叶片子区域的网格,且相邻两交界面处的网格大小一致。为节省计算资源,可适当增大静止区域远场处的节点间距。综合考虑计算的准确性及计算机的可实现性,并验证网格无关性,最终将计算区域的网格总数定为435 468,其中转动区域网格总数为268 066,静止区域网格总数为167 402,动网格约占网格总数的61.5%。

图3 计算区域网格

2 求解设置

2.1 数值计算方法

用FLUENT 软件对垂直轴风力机流场进行瞬态计算,采用k-ωSST湍流模型[5-6],用SIMPLEC算法对速度和压力进行耦合,以加快收敛速度,动量方程、湍动能和湍流离散率采用二阶迎风格式,以保证计算精度。设定风轮每旋转0.005 s(即风轮每旋转1.2°)计算一次,则风轮旋转一周有300个时间步,在每个时间步内迭代20次直至在该时间步内迭代收敛,残差收敛标准为各个残差值小于10-6。当相邻2个旋转周期的转矩系数偏差小于1%时认为计算已经收敛[7-8]。

2.2 边界条件设置

设置进口为速度进口,方向与进口边界垂直;出口为压力出口,计示压强为0,即出口区域为未受扰动的标准大气。叶片表面为无滑移的固壁面(wall)边界条件,它和叶片子区域同步运动;静止区域和旋转区域、旋转区域和叶片子区域设置为交界面(interface)边界条件,保证转动网格与静止网格能量和动量的传递。静止区域上下两表面设置为对称面(symmetry)边界条件。

3 计算结果分析

3.1 二维计算域流场分析

入口风速为10 m/s,风轮转速为42 r/min时,二维计算域CFD速度分布云图如图4所示。以风轮最上方的叶片为研究对象,其初始方位角为0°,由于受上个周期尾流的影响,初始时刻该叶片的来流速度小,且攻角较小,故叶片上下表面速度分布差异不大,气动力较小;叶片方位角在0°～90°范围内变化时,叶片攻角逐渐增大,叶片上下表面速度分布差异增大,引起上下表面压力分布不均匀,气动力逐渐增大;叶片方位角在90°～180°范围内变化时,攻角先增大后减小,在攻角增大期间,叶片出现失速现象,气动力急剧减小,而攻角减小时,气动力随之逐渐减小;叶片方位角在180°～360°范围内变化时,该叶片处于上游叶片的尾流区域,叶片在前180°范围内产生并不断增强直至脱落的旋涡,在向下游发展的过程中,影响通过下游叶片的气流,限制下游叶片功率的捕获;方位角在180°～270°范围内变化时,叶片逐渐进入尾流区,攻角负向增大,气动力逐渐增大;叶片方位角在270°～360°范围内变化时,攻角先增大后减小,叶片完全进入尾流区,尾流影响增强,在攻角增大期间,叶片失速,气动力急剧减小,而攻角减小时,气动力随之逐渐减小。

3.2 三维计算域流场分析

洁净空气下的三维H型垂直轴风力发电机在入口风速为10 m/s、风轮转速为40 r/min时,采用与二维计算域相类似的计算方法得到的计算结果,分别用x=0,y=0,z=0 3个平面将计算域剖切开,得到截面处的速度场分布,如图5所示。

图5 三维计算域速度分布

三维垂直轴风力机绕流流场相对于二维流场流动更加复杂,表现出很强的非定常性。沿叶展方向(z向)速度分布对称,但由于气流在展向方向的相互影响,各叶片在不同水平面上(z=常数)速度并不相同,故平面上出现显著变化的压力分布,叶片上下表面沿叶展方向的压力变化趋势相同。就整个风轮而言,上游叶片的失速更为明显,上游叶片的尾流影响绕下游叶片的流动,导致下游叶片捕获的风能减小,所以整个风轮在每个旋转周期的风能利用系数不高,最大风能利用系数约为0.32。

洁净空气下的三维H型垂直轴风力发电机在入口风速为10 m/s、风轮转速为40 r/min时,分别用x=0,y=0,z=0 3个平面将三维计算域剖切开,得到截面处的涡量分布如图6所示。涡量即为速度的旋度,每个叶片在同一平面和叶展方向上的旋涡不断叠加,沿叶展方向,在叶片两端即叶尖处的涡量最大。就风轮而言,上游叶片处的旋涡强度高于下游叶片处的,上游叶片的旋涡相互叠加,在尾流中不断发展,使流经风轮的能量不断衰减,严重影响下游叶片风能的捕获。

图6 三维计算域涡量分布

3.3 风沙环境下流场分析

设置风沙来流质量流量为0.3 kg/s,颗粒物直径为0.15 mm,风沙来流速度为10 m/s,选择DPM模型进行仿真计算。分别用x=0,y=0,z=0 3个平面将三维计算域剖切开,得到风沙环境下截面处的涡量分布如图7所示。

图7 风沙环境下三维计算域涡量分布

相对于洁净空气,风沙环境下,在叶尖处的涡量增加,且涡量的最大值出现在上翼面的后半段,主要是由于边界层分离后,固体颗粒随气体沿翼面出现倒流,与随气流沿翼面顺流的固体颗粒在翼面后半段汇合,从而出现颗粒物堆积,改变翼型的几何形貌,形成“锋面”,使得该处的湍流强度增加,涡量增大。在风沙环境下,进行叶片气动外形设计时,最好选择对称的相对厚度较小的翼型(如NACA0012),减小翼型相对厚度,可减小边界层厚度增加趋势,延缓边界层分离。

4 结论

针对小型H型垂直轴风力发电机的气动性能,采用计算流体力学中的滑移网格技术,对叶片及风轮的流场进行计算,得到如下结论:

(1) 对于小型H型垂直轴风力发电机,由于其攻角变化范围大,流动复杂,其风轮的三维流场表现出具有显著非定常性的分离流动。

(2) 洁净空气下,风力机叶片在叶尖处的涡量最大,上游叶片处的失速更加明显,上游叶片的旋涡相互叠加,使尾流影响绕下游叶片的流动,减小风能的捕获,整个风轮的最大风能利用系数约为0.32。

(3) 风沙环境下,叶尖处的涡量进一步增大,且涡量最大处在上翼面的后半段,在考虑风沙作用进行叶片设计时,应选择对称的相对厚度不大的翼型,如NACA0012,从而延缓边界层分离,增加气动力。