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]
(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。
基于上述讨论可以得到翼型的静态流动分离点，翼型处于动态失速状态时，需要修正静态流动分离点fN和fT，得到与前缘压力相应滞后有关的后缘分离点
和
修正过程需要引入有效分离攻角α′：
(10)
(11)
式中，
为滞后法向力系数，可描述法向力系数CN对应的前缘压力的滞后效应；Tp为B-L模型中定义的时间常数。
将有效分离攻角α′代替式(8)、式(9)中的α，便可得到和
然后通过对和
的一阶延迟修正来得出最终的非定常尾迹的动态分离点：
(12)
(13)
式中，Tf是B-L模型中定义的半经验时间常数。
最后，将式(12)、式(13)得到的
和
代入式(6)、式(7)，便可得到流动分离非线性法向力系数
和切向力系数
(14)
(15)
1.3 动态失速涡
翼型动态失速的显著特点是在翼型前缘产生分离涡，该分离涡形成后，将在翼型上表面沿着弦线向后缘方向移动，这一过程导致了翼型的气动力系数发生非定常变化。B-L模型引入当量值Cv来描述涡流的产生、发展和分离对气动力系数的影响。Cv是非定常绕流法向力系数
和流动分离非线性法向力系数的差值：
(16)
总的涡法向力系数
随时间按指数规律减小，但同时也因新涡流的出现而增大，因此涡流法向力系数的变化规律可表示[12]为
(17)
求解该式便可得到涡流法向力系数。
文献[13]指出，翼型在大攻角情况下，涡流对切向力的影响较大，需要加入B-L模型加以考虑。涡流对切向力系数的影响可表示为
(18)
式中，τV为量纲一涡时间。
B-L模型在确定翼型整体的动态切向力系数和动态法向力系数时，还要考虑该翼型在零升力攻角的阻力系数CD0。翼型发生动态失速时，总的动态法向力系数
和切向力系数
可由式(16)～式(18)中的分量叠加得到：
(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)
式中，γ为攻角变化曲线的斜率；
为翼型攻角变化率；
是攻角关于时间的导数。
(3)攻角超出动态失速攻角αds。翼型发生动态失速时，翼型的前缘附近形成失速涡，失速涡沿翼型的上表面运动并在后缘脱落，这一过程伴随有翼型升力系数的迅速增大。MIT模型假设这一过程中，翼型的动态升力系数
和动态阻力系数
始终保持最大值CLmax和CDmax：
(23)
CDmax=CLmaxsinα
直到攻角开始减小。这种假设极大地简化了计算，但对模型的计算精度有一些影响。
(4)翼型攻角从αds开始下降时，MIT模型假设动态升力系数
和动态阻力系数
分别为
式中，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为叶片高度。
原双致动盘多流管模型(式(29))中，切向力系数CT和法向力系数CN的值是通过查询翼型升阻力系数数据库得到的。垂直轴风力机运行时，叶片攻角随方位角不断变化，导致叶片处于失速状态，静态的CT和CN已经不能反映叶片真实的受力状态，此处应代入动态切向力系数
和动态法向力系数
文中，
通过B-L修正模型和MIT修正模型求解，因此可以建立2种计及动态失速效应的双致动盘多流管模型。
当叶片处于上风区时，上风函数为
(30)
式中，N为叶片数目；n为双致动盘多流管理论中的流管数目。
通过fupu=1-u可迭代计算出上风区诱导因子u，进而可计算出上风区诱导速度、合成风速，以及各方位下的实际攻角值，最后得到叶片在上风区各个方位下的动态气动力系数。
同理，当叶片处于下风区时，计算方法与上风区类似，下风函数为：
(31)
将上风区迭代得到的各个流管的诱导因子作为下风区迭代公式fdwu′=1-u′的迭代初值进行计算，便可得到下风区各流管的诱导速度、合成风速，以及各方位下的实际攻角值，最终得到叶片在下风区各个方位下的动态气动力系数。
4 算例求解与分析
以美国Sandia实验室17 m垂直轴风力机[14-16]为算例进行计算，其基本参数如表1所示。
表1 风力机参数
Tab.1 Parameters of wind turbine
名称取值翼型NACA0015弦长(m)0.53风轮高度(m)17风轮直径(m)16.7叶片数2空气密度ρ(kg/m3)1.001 2空气运动黏度υ(m2/s)1.79×10-5
将表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修正模型。
[1] PARASCHIVOU I．Wind Turbine Design with Emphasis on Darrieus Concept[M]．Montreal: Ecole Polytechnique de Montreal, 2002．
[2] 杨从新, 巫发明, 张玉良. 基于滑移网格的垂直轴风力机非定常数值模拟[J]. 农业机械学报, 2009, 40(6):98-102.
YANG Congxin, WU Faming, ZHANG Yuliang. Numerical Simulation on Unsteady Rotated Flow of a Vertical Axis Wind Turbine Based on Moving Meshes[J]. Transactions of the Chinese Society for Agricultural Machinery, 2009, 40(6): 98-102.
[3] BRIAN H, GER K, ANDREW C. Numerical Simulation of a Vertical Axis Wind Turbine Airfoil experiencing Dynamic Stall at High Reynolds Numbers[J]. Computers and Fluids, 2017, 149: 12-30.
[4] 李银然, 李仁年, 王秀勇,等. 计算模型维数对风力机翼型气动性能预测的影响[J]. 农业机械学报, 2011, 42(2): 115-119.
LI Yinran, LI Rennian, WANG Xiuyong, et al. Effects of the Calculation Models with Different Dimension on the Aerodynamic Performance Prediction for Wind Turbine Airfoil[J].Transactions of the Chinese Society for Agricultural Machinery, 2011, 42(2): 115-119.
[5] VASILIS A R, SPYROS G V．Dynamic Stall Modeling on Airfoils Based on Strong Viscous-inviscid Interaction Coupling[J]. International Journal for Numerical Methods, 2007, 30(1): 153-162．
[6] LARSEN J W, NIELSEN S R．Dynamic Stall Model for Wind Turbine Airfoils[J]. Journal of Fluids and Structures, 2007, 23(7): 959-982．
[7] 刘占芳,颜世军,张凯,等. 立轴风力机叶片动态失速特性与气动性能分析[J].太阳能学报, 2012, 33(2): 204-209.
LIU Zhanfang，YAN Shijnu, ZHANG Kai, et al. Dynamic Stall Characteristic and Aerodynamic Performance Prediction for Blade of VAWT[J]. Acta Energiae Solaris Sinica, 2012, 33(2): 204-209.
[8] TRAN C T, PETOT D. Semi-empirical Model for the Dynamic Stall of Airfoils in View of the Application to the Calculation of Response of a Helicopter Blade in Forward Flight[J]. Vertica, 1981, 5(1): 35-53.
[9] LEISHMAN J G，BEDDOES T S．A Semi-empirical Model of Dynamic Stall[J]. Journal of American Helicopter Society, 1989, 34(3): 3-17．
[10] SANDEEP G, LEISHMAN J G．Dynamic Stall Modelling of the S809 Airfoil and Comparison with Experiments[J]. Wind Energy, 2006, 9(1) : 521-547.
[11] HANSEN M H, GAUNAA M, MADSEN H A. A Beddoes-Leishman Type Dynamic Stall Model in State-space and Indicial Formulations [Risφ-r-1354][R]. Copenhagen, Denmark: Risφ National Laboratory, 2004.
[12] 刘雄,梁湿,陈严,等. 风力机翼型动态失速气动特性仿真[J].工程力学, 2015, 32(3): 203-211.
LIU Xiong, LIANG Shi, CHEN Yan, et al. Dynamic Stall Simulation of Wind Turbine Airfoils[J]. Engineering Mechanics, 2015, 32(3): 203-211.
[13] 伍艳,谢华,王同光.风力机叶片的非定常气动特性计算方法的改进[J]. 工程力学, 2008, 25(10): 54-59.
WU Yan, XIE Hua, WANG Tongguang. Modification of Calculating Unsteady Aerodynamic Characteristics of Wind Turbine Blades[J]. Engineering Mechanics, 2008, 25(10): 54-59.
[14] PIERCE K G. Wind Turbine Load Prediction Using the Beddoes-Leishman Model for Unsteady Aerodynamics and Dynamic Stall[D]. Salt Lake: The University of Utah, 1996.
[15] ALLET A，PARASCHIVOIU I．Viscous Flow and Dynamic Stall Effects on Vertical-axis Wind Turbines[J]. International Journal of Rotating Machinery, 1995, 2(1): 1-14．
[16] ROBERT E S，PAUL C K．Aerodynamic Characteristics of Seven Symmetrical Airfoil Sections through 180-degree Angle of Attack for Use in Aerodynamic Analysis of Vertical Axis Wind Turbines，SAND-80-2114[R]. Albuquerque，1981.
