垂直轴风力机叶片改进动态失速模型

张立军 赵昕辉 马东辰 米玉霞 王旱祥 姜 浩

中国石油大学(华东)机电工程学院，青岛，266580

0 引言

垂直轴风力机可以接收来自任何方向的风，其增速齿轮箱和发电机可以安装在地面，运行维修方便，但风力机运行时，其叶片常发生动态失速(低叶尖速比时尤为明显)[1]，因此准确计算动态失速下的叶片气动力系数是分析与设计垂直轴风力机的关键。动态失速是指叶片攻角发生周期性或非定常变化时，翼型的失速攻角比静态失速攻角要大得多，且翼型气动特性曲线(通常为法向力系数和切向力系数随攻角的变化曲线)要明显滞后于静态曲线的现象。翼型发生动态失速时测得的动态气动力系数与翼型静止时的静态气动力系数相差较大。

目前研究翼型动态失速的方法主要有三种：①以Navier-Stokes方程为基础的CFD数值方法[2-4]；②基于面元法和边界层理论的黏性与无黏耦合算法[5]；③基于实验数据建立的半经验动态失速模型方法[6-7]。其中，半经验动态失速模型计算效率高且通用性好，通过适当的修正便可用于垂直轴风力机翼型的动态气动系数计算。现有的半经验动态失速模型主要有以下几类[7-11]：基于对动态失速延迟数值修正的B-V模型和MIT模型；用三阶微分方程描述气动系数的ONERA模型；基于翼型绕流流动特性的B-L动态失速模型、丹麦Ris国家实验室建立的用于评价水平轴风力机气动特性的Ris模型。国内外一些科研单位对上述模型进行适当修正，并用于各自研究领域[12-14]。这些模型均是针对航空翼型或水平轴风力机翼型的，在计算垂直轴风力机叶片动态失速气动性能时，需要结合其实际工作特性对这些模型进行修正，但相关研究较少。

由于B-L模型和MIT更多地考虑了运动翼型的绕流物理特性，能更好地模拟翼型动态失速特性且实用性较强，因此本文首先介绍了目前常用的B-L模型和MIT模型的计算方法，并基于垂直轴风力机工作时的非定常气动特性对这两种模型进行修正。然后结合双致动盘多流管理论，用修正后模型分别求解Sandia实验室17 m垂直轴风力机叶片的动态切向力系数和动态法向力系数，将计算结果与实验数据进行对比，分析这两种修正模型对不同风区的预测精度。最终得到一种在整个风区都能保持较高计算精度的垂直轴风力机叶片动态气动力系数计算模型。

1 B-L模型修正

B-L模型[9]属于半经验模型，它虽依赖经验常数，但能部分预测动态失速现象。B-L模型将坐标系OXY固定在翼型上，如图1所示。翼型的升力系数CL和阻力系数CD向弦线方向(X方向)和弦线法向(Y方向)投影，可得切向力系数CT、法向力系数CN。图1中，W是作用在翼型上的合成风速，可通过双致动盘多流管理论求得；α为翼型攻角，是合成风速方向与翼型弦线方向的夹角。

图1 翼型的气动力系数Fig.1 Aerodynamic coefficients of airfoils

B-L模型主要是确定非定常附着流、分离流动和动态失速涡对翼型气动力系数的影响，在此基础上得出翼型在动态失速下的动态气动力系数。

1.1 非定常附着流

(1)

αE=φC(t)Δα+αn-1

φC(t)=A0-A1exp(-b1s)-A2exp(-b2s)

Δα=αn-αn-1

式中，CNα为对应翼型的静态法向力系数曲线线性段的斜率；αE为B-L模型中定义的有效攻角；φC(t)为指数形式的延迟函数；t为翼型的厚度；s为量纲一时间，s=2vt/c；v为平均来流风速；c为翼型弦长；αn-1为翼型在n-1时刻的攻角；常数A0=1.0，A1=0.30，A2=0.70，b1=0.35，b2=0.68[7]。

原始的B-L模型采用阶跃函数来描述翼型攻角随时间的变化规律。对于垂直轴风力机而言，其翼型攻角α的变化规律为类正弦曲线：

(2)

式中，θ为翼型当前所处的方位角；λ为风力机局部叶尖速比，λ=V/(ωR)；V为流体诱导速度；ω为风轮旋转角速度；R为风轮旋转半径。

本文使用式(2)代替原B-L模型的阶跃函数，既符合垂直轴风力机工作的实际情况，也使攻角变化的连续性更强。

脉冲分量反映了翼型运动引起的表面法向速度的气动响应，可以由活塞理论进行求解：

(3)

式中，Ma为马赫数。

对垂直轴风力机而言，Ma≪1，式(3)的计算结果接近于0，故可以忽略附着流的脉冲效应[12-13]，进而可对现有的B-L模型进行简化。因此非定常附着流的法向力系数为

(4)

切向力系数为[7]

(5)

式中，CLα为对应翼型的静态升力系数的曲线线性段的斜率。

1.2 分离流动

对分离流动的描述是B-L模型中至关重要的一部分，而描述分离流动最关键的是确定流动分离点的位置。通过Kirchhoff流动定理可以得到静态流动分离点与翼型在定常情况下的法向力系数、切向力系数：

(6)

(7)

式中，fN、fT分别为流动分离点在翼型表面法向和切向的位置，fN=x/c，fT=x/c；x为流动分离点距后缘的长度；α0为零升攻角。

将式(6)、式(7)重新整理可得对应的法向力系数和切向力系数的分离点位置：

(8)

(9)

在小攻角变化情况下，fN和fT的值非常接近，故在原B-L模型中，令流动分离点f=fN，其中，fN为式(8)的计算结果，并将f作为计算CT和CN共同的流动分离点参数，再通过函数拟合方法得到分离点f与攻角α的连续函数关系。这样对于任意的攻角α就都可以得到一个流动分离点f，进而通过式(6)、式(7)得到对应攻角下的静态法向力系数和切向力系数。对于攻角变化范围较大的垂直轴风力机翼型而言，通过上述方法得到的法向力系数和切向力系数与原始数据不能很好地吻合，文献[12]对这一现象做了相应的解释。为此，本文在将B-L模型用于垂直轴风力机翼型动态失速性能研究时，独立计算静态法向力分离点fN和静态切向力分离点fT，然后分别代入式(6)、式(7)计算法向力系数CN和切向力系数CT。

(10)

(11)

(12)

(13)

式中，Tf是B-L模型中定义的半经验时间常数。

(14)

(15)

1.3 动态失速涡

(16)

(17)

求解该式便可得到涡流法向力系数。

文献[13]指出，翼型在大攻角情况下，涡流对切向力的影响较大，需要加入B-L模型加以考虑。涡流对切向力系数的影响可表示为

(18)

式中，τV为量纲一涡时间。

(19)

(20)

2 MIT模型修正

对于垂直轴风力机，翼型攻角、气动力和风速关系如图2所示，V是诱导速度，W是诱导速度V和切向速度ωR的合成风速，FL、FD分别为叶素受到的升力和阻力，切向力FT、法向力FN分别为FL和FD在翼型弦线方向和法向上的分量。

图2 翼型气动参数分析示意图Fig.2 Sketch map of airfoil aerodynamic parameters analysis

MIT模型将攻角α的变化过程划分为4个阶段，分别用不同的公式计算当前攻角下翼型的动态切向力系数和动态法向力系数，具体计算方法如下：

(1)翼型攻角α从0开始增大且低于静态失速攻角αss时，动态升力系数和阻力系数均取静态值，这些值可通过查询翼型升/阻力系数数据库得到。数据库未给出的值，可通过线性插值的方法求得。

(2)当攻角继续增大，超过静态失速攻角αss但还未到达动态失速攻角αds时，动态阻力系数仍采用对应攻角的静态值，而升力系数为

(21)

式中，CLss为静态失速攻角αss对应的静态升力系数。

MIT模型用诱导速度V作为翼型的特征风速。对于垂直轴风力机而言，作用在叶片上的风速为诱导速度V和切向速度ωR的合成风速W(其大小为W)。动态失速攻角αds定义为

(22)

(23)

CDmax=CLmaxsinα

直到攻角开始减小。这种假设极大地简化了计算，但对模型的计算精度有一些影响。

式中，CLs、CDs分别为翼型在当前攻角下的升力系数和阻力系数的静态值，均以指数方式衰减至其静态值。

当攻角再次增大时，重复上述过程，便可得到各个攻角下的动态升阻力系数。通过下式可得到动态切向力系数和动态法向力系数：

(26)

3 修正模型中风速参数计算

B-L修正模型和MIT修正模型都需要求解各个方位下作用在叶片上的合成风速。本文中，合成风速采用双致动盘多流管理论计算得出。双致动盘多流管模型用2个串联的致动盘分别表示上风区的半个转子扫掠面和下风区的半个转子扫掠面，且假设各致动盘的诱导速度不同，如图3所示。

风速为v∞的无穷远来流经垂直轴风轮时会经历速度衰减过程，即流经上风区时，速度变成上风区诱导速度V；流经转轴处时，速度进一步降低至均衡诱导速度Ve；流经下风区时，速度为下风区诱导速度V′，三者速度之间的关系为

V=uv∞

(27)

V′=u′(2u-1)v∞

(28)

式中，u、u′分别为上风区诱导因子和下风区诱导因子，u<1，u′<1。

对于垂直轴风力机， 采取分风区讨论的方法

图3 双致动盘多流管模型Fig.3 Model double-multiple stream-tubes

来研究叶片受力情况。如图4所示，方位角θ位于0°～180°的区域称为上风区；方位角θ位于180°～360°的区域称为下风区。

图4 叶片受力分析图Fig.4 Analysis of blade stress

叶片所受的切向力FT和法向力FN分别为

(29)

式中，ρ为空气密度；H为叶片高度。

当叶片处于上风区时，上风函数为

(30)

式中，N为叶片数目；n为双致动盘多流管理论中的流管数目。

通过fupu=1-u可迭代计算出上风区诱导因子u，进而可计算出上风区诱导速度、合成风速，以及各方位下的实际攻角值，最后得到叶片在上风区各个方位下的动态气动力系数。

同理，当叶片处于下风区时，计算方法与上风区类似，下风函数为：

(31)

将上风区迭代得到的各个流管的诱导因子作为下风区迭代公式fdwu′=1-u′的迭代初值进行计算，便可得到下风区各流管的诱导速度、合成风速，以及各方位下的实际攻角值，最终得到叶片在下风区各个方位下的动态气动力系数。

4 算例求解与分析

以美国Sandia实验室17 m垂直轴风力机[14-16]为算例进行计算，其基本参数如表1所示。

表1 风力机参数

将表1中的数据代入式(30)、式(31)，计算该垂直轴风力机的叶尖速比分别为2.33和3.09时，叶片的动态切向力系数和动态法向力系数。

对比图5可以发现，当叶片处于上风区且B-L修正模型和MIT修正模型在方位角0°～60°时，计算结果与实验数据吻合度较好；方位角60°～180°时，B-L模型的计算结果急剧减小，而MIT修正模型的计算结果基本反映实验数据的变化趋势。当叶片处于下风区时，B-L修正模型的计算结果与实验数据相差不大，而MIT修正模型的计算结果除在方位角210°～240°、330°～360°能吻合外，其他范围与实验结果的一致性较差，不推荐使用。

对比图6可以发现，这两种修正模型在计算动态法向力系数时，总体上可以反映实验结果的变化规律，但在一些方位角下，两种修正模型各有优缺点：对于上风区，B-L修正模型的总体趋势预测精度要高于MIT修正模型；对于下风区，B-L修正模型在方位角210°附近的预测值与实验结果相差较大，而MIT修正模型与实验结果的一致性较好。

(a)=2.33

(b)=3.09图5 两种叶尖速比的动态切向力系数Fig.5 Dynamic tangential force coefficientsof 2 tip speed ratios

(a)=2.33

(b)=3.09图6 两种叶尖速比的动态法向力系数Fig.6 Dynamic normal force coefficients of 2 tip speed ratios

基于上述对比分析，结合双致动盘多流管理论可以得到一种改进的垂直轴风力机叶片动态失速性能计算模型。该模型计算动态切向力系数时，上风区使用MIT修正模型，下风区使用B-L修正模型；计算动态法向力系数时，上风区使用B-L修正模型，下风区使用MIT修正模型。修正模型的计算流程如图7所示。

图7 垂直轴风力机叶片动态失速性能计算修正模型Fig.7 Correction model for dynamic stall performance calculation of vertical axis wind turbine blades

5 结论

(1)结合垂直轴风力机的非定常气动特性，对B-L和MIT动态失速模型进行了修正：简化了原B-L模型中的马赫数相关项及翼型非绕流升力部分，改进了流动分离点的计算方法，将原模型中的攻角变化规律修正为垂直轴风力机运行时攻角的实际变化规律；将MIT模型中的特征风速定义为垂直轴风力机工作时的实际风速，对动态失速攻角及最大升阻力系数的计算公式进行了修正，使该模型更加适合垂直轴风力机叶片的实际工况。

(2)鉴于垂直轴风轮流场的复杂性，采用分风区的方法研究叶片的动态受力。对比2种修正模型的计算结果和Sandia实验室17 m垂直轴风力机实验数据发现：MIT修正模型对上风区的切向力系数和下风区的法向力系数的预测精度较高；B-L修正模型对上风区法向力系数和下风区切向力系数的计算结果与实验数据吻合良好。

(3)该模型在计算动态切向力系数时，上风区使用MIT修正模型，下风区使用B-L修正模型；在计算动态法向力系数时，上风区使用B-L修正模型，下风区使用MIT修正模型。