基于升力线理论的风力机气动性能分析方法<sup>*</sup>

目前，在风力机叶片设计和机组气动性能分析中，普遍采用叶素

动量理论(BEM)或计算流体力学模型(CFD)[1]，如美国Sandia实验室开发的AeroDyn程序[2].BEM理论忽略了沿叶片的展向流动，不能反映叶素单元间的相互影响，无法预测复杂三维流动以及叶尖和轮毂损失等效应，需要引入大量的人为经验修正模型，这在一定程度上降低了该理论的普适性和可靠性，尤其是对于三维、非定常、偏航等工况[3-4].还有部分研究[5-6]采用CFD方法，该方法基于Navier-Stocks方程的数值求解，不仅可以给出详细的三维流场信息，而且可以准确地计算湍流等复杂流场结构对风力机气动性能的影响.但由于风力机三维流场的复杂性和多尺度效应，该方法存在计算量巨大的特性，目前，风电领域的应用集中在二维翼型的分析与优化设计，以及定常来流下的固定工况点的定常计算等方面[7-8]，还未能在风力机的叶片设计中得到广泛应用.

80年代Jeng提出升力线算法[9]，将叶片用一根环量连续变化的升力线来代替，根据霍姆亥兹涡量守恒定理，从升力线上拽出涡丝张成螺旋尾涡面，即由附着涡(升力线)和尾涡构成流动模型.诱导速度由Biot-Savart定理求解，由涡丝相互之间的诱导考虑有限叶片的效应，不需引入任何修正，是一种简洁有效的算法.

近年来，风力机组呈大型化发展趋势，气动载荷、重力、惯性载荷及操纵载荷随着叶片长度的增加急剧上升[10]，需要提高气动模型的三维计算能力.叶片近壁区域流场的展向流动和尾涡系统的三维运动等因素将使基于二维流动假设的BEM模型面临困难，需要对该模型引入人工修正模型，如叶根叶尖损失修正模型和三维无粘延迟失速模型等，才能得出符合实际的结果.而对升力线模型来说，虽然也需要根据二维翼型的升阻力曲线使其计算封闭，但通过Biot-Savart公式计算三维诱导速度在一定程度上提高了该模型的三维流场分析能力，尤其是尾涡模型的引入，大大地增加了模型对尾流分析的可靠性.

基于以上分析，本文采用螺旋尾涡升力线模型来计算叶片气动性能参数.通过对附着涡分布、控制点的诱导速度以及迭代法求解升力系数等问题进行研究和分析，建立了基于螺旋尾涡升力线模型的风力机气动性能数值分析算法，并在FORTRAN平台中创建了分析程序.算例分析了NREL实验室公布的5 MW风力机[11]的气动性能，并与NREL发布的结果进行了对比验证.

1 风力机叶片升力线模型

由伯努利的流体机械能守恒原理可知，由于流过翼型压力面和吸力面的流速不同导致作用在压力面和吸力面产生压力差，从而使翼型产生升力.根据Kutta-Joukowski定理，翼型的升力L可以由自由来流风速w与翼型环量Γ之间的相互作用来表示，即

式中：ρ为当地空气密度；dr为叶片沿展向方向的微段.根据这一思想，升力线模型利用一根涡量沿展向变化的线涡替代叶片，表征叶片与流场之间的相互作用，示意图如图1所示.

在升力线模型中，涡量沿展向分布的线涡代表了三维流场中真实三维叶片的附着环量分布，该线涡称为集中附着涡.此外，由Helmholtz第二定理可知，在无粘流体环境中，涡量并不能在流体内部终止，而只能延伸至流体边界或构成环.因此，附着涡沿展向的变化量并不能消失而必然脱落流向下游，形成尾流涡.升力线模型以一根集中涡表征叶片的空气动力作用，对于该集中涡来说，只能感受到力的作用而不能感受到力矩的作用，因此需要将该集中涡放置在翼型的某一点，使得该点的俯仰力矩系数不随攻角的改变而变化，而这样的点恰好为翼型的气动中心[12].

升力线模型通过计算附着涡线上参考点位置处的诱导速度，获取叶片在展向位置某一处的有效攻角，并通过插值获得该处翼型的升阻力系数，从而实现当前翼型的升阻力计算；再通过在整个叶片展向上进行积分，实现叶片的扭矩和功率等气动性能参数的计算.

1.1 诱导速度计算

直线线涡对空间内目标点位置(即控制点)处产生的诱导速度可用Biot-Savart公式描述，即

式中：wt为叶片控制点处的诱导速度；γ为流场中分布的涡量；r为叶片半径向量；ds为线涡微段矢量.

由于本文采用的是刚性尾涡模型，该涡线为螺旋状，螺旋尾涡向叶片下游延伸至无限远处，形成诱导速度场，坐标系如图2所示.

将螺旋尾涡离散为若干段直线线涡，由式(2)可推出由螺旋尾涡引起的叶片展向位置处的诱导速度，其表达式为

式中：s为轮毂中心指向微端dη的向量；wi为叶片展向位置处的诱导速度；r′为轮毂中心指向目标位置点的向量；dη为第k个叶片展向位置r处发散出来的螺旋线涡微段矢量；Γ为附着涡的环量.将dwi沿着坐标轴分解成分量形式，可得

根据刚性尾涡理论可知，dwx=0.若将诱导速度wi沿着垂直和平行于相对风速w′方向进行分解(如图3所示)，则有

式中：φ为入流角；V0为来流风速；Ω为叶片转速.

诱导速度总量可以通过式(5)、(6)从叶片轮毂到叶尖处积分得到.Moriya证明了在叶片上任意位置均有wt=0，设叶片上展向无量纲位置参数为

式中，R为叶片半径.叶片展向位置任意位置ξ′处，由所有尾涡引起的总诱导速度为

式中：λ0=V0/RΩ为尖速比；ψ为叶片锥角；RL= width=704,height=111,dpi=110

为叶片轮毂半径；θk为从第k个叶片展向位置r处发散出来的螺旋线涡的方位角，若存在N个叶片，则有

其中，K为从控制点所在叶片的下一个叶片(任意方向均可，但方向需要保持一致)开始算起，对应的第K个叶片(如图2所示).

当θ=θk=0时，积分在点ξ=ξ′处的值变为无穷大，产生奇点，而且奇点附近的积分值也无法确定，为此，本文引进一个额外的参数I，即诱导因子，其定义为

式中：wn2和wn1分别为涡量强度相同的螺旋线涡和直线涡产生的诱导速度；I为连续函数.

由文献[9]可知，诱导因子I是关于ξ、ξ′和λ0的函数，即

当θ=θk=0时，积分在点ξ=ξ′处的值同样会变为无穷大，实际上可以近似地认为在这个点处dwn2=dwn1，即I=1，且由于诱导因子为连续函数，可以通过在奇点附近进行样条插值来确定奇点以及其附近的诱导因子的值.

当叶片上各点处的诱导因子确定之后，其诱导速度也可确定，即

诱导速度与各攻角的关系如图4所示，其中，图4中各变量的表达式为

根据图4可以得到

设环量为傅里叶正弦级数，并且满足边界条件Γ(ξhub)=Γ(1)=0，则环量函数为

式中，无穷级数求和可近似用有限项求和来计算，项数m越多，环量Γ(ξ)的计算越精确，一般m取为计算控制点的数目.

根据Kutta-Joukovski定理可知，环量与升力系数之间的关系为

式中：c为翼型弦长；CL为升力系数.对于每一个有效攻角，均可以通过二维翼型数据表并通过线性插值得到一个与之相对应的升力系数，即CL=CL(αe).

对于初始计算，合成风速w即为来流风速V0，根据式(19)可得

在小攻角范围内，升力系数与有效攻角之间近似为线性关系，即

式中，a0、b0为常数.综合式(12)、(17)、(18)、(21)可得

将方程(22)应用于升力线各计算点ξ′处，即可得到M个方程，求解该方程组可以得到Am的值，将Am代入式(18)即可得到环量Γ的值；将式(21)代入式(19)，可得有效攻角与环量之间的关系，即

由环量值即可求出有效攻角αe.

1.2 迭代法求解有效攻角

一般采用迭代求解有效攻角αe，其迭代步骤为：1)利用得到的有效攻角值，通过查询二维翼型数据表用线性插值法求出对应的升力系数CL；2)由式(19)计算出环量；3)求出环量之后，由式(18)计算出Am；4)根据Am的值便可通过式(22)得到诱导攻角αi；5)根据诱导攻角、有效攻角与几何攻角之间的关系，得到新的有效攻角αe；6)通过查询二维翼型数据表用线性插值法求出新的升力系数CL；7)将新的升力系数与前一次迭代所得到的升力系数进行比较，若满足一定的收敛条件则迭代结束，否则

跳转至步骤2)，继续迭代过程.

迭代结束之后，可以计算出叶片的轴向力F、力矩Q以及功率P等叶片气动性能参数，其表达式为

式中，CD为阻力系数，对于每一个有效攻角，都可以通过二维翼型数据表并通过线性插值得到一个与之相对应的阻力系数CD.本文中ψ取为0°.

2 算例分析

根据以上理论分析与数值求解分析方法，在FORTRAN平台上建立了完整的气动性能数值计算程序，实现输入已知参数即可求得多叶片风轮的气动性能参数.

为了验证建立的升力线模型的准确性，以NREL实验室公布的5 MW风力机叶片数据为研究案例，对其气动性能进行计算，并将计算结果与NREL实验室提供的结果(基于叶素动量理论进行的计算)进行对比.

为了保证计算效率和计算结果的可靠性，离散控制点数和位置的选取，一般可根据叶素动量理论计算控制点的选取方法，理论上控制点越多，相应的计算精度越高，但一般不宜取得太多，否则增大了计算量且对精度没有太大提高[9].本文计算中，根据周传捷等[13]的研究可知，将叶片离散为13段，既能保证计算效率，又能得到较为满意的结果.由于叶片靠近轮毂中心的前两段为圆筒状，其升阻力系数在任意攻角下均为恒定值，故其升阻力系数与攻角的函数曲线斜率为0，在进行式(12)的计算时分母为0，会导致结果变为无穷，使得计算无法进行，因此，可以将其剥离出来单独讨论，由于升阻力系数均为已知，故可通过式(24)、(25)直接计算出前两段圆筒的轴向力以及力矩.

基于以上分析，假定前两段圆筒形的环量均为0，则可将叶片第三段起始截面假想为轮毂半径ξhub，因此，叶片的扭角和弦长分布如图5所示.

根据NREL公布的5 MW风力机叶片二维翼型数据表[11]，各翼型的升力系数曲线的斜率以及截距如表1所示.

将上述已知参数作为输入文件导入用FORTRAN编写的升力线模型程序，计算出该5 MW风力机在不同风速下的气动力特性.图6为不同风速下叶片展向位置的环量分布图.

根据文献[11]可知，当风速低于额定风速时，无需对桨距角进行调整，而在更高的风速下，为了使风力机产生稳定的机械功率，在恒定的转速12.1 r/min下需对叶片的桨距角进行相应调整，其相应风速下的变桨距角为文献[11]通过风洞试验所得.表2为不同风速下变桨距后利用升力线模型所求得的风力机轴向力、功率以及扭矩.

将不同风速下采用升力线理论计算所得结果中的推力、功率和扭矩，与NREL实验室公布的5 MW风力机采用引进修正模型的叶素动量理论计算所得的相应参数[11]进行对比，结果如图7所示.

由图7可以看出，用升力线理论计算得到的风力机的功率、力矩、轴向力与文献[11]基于修正的叶素动量理论计算结果比较接近.通过NREL实验室公布的5 MW风力机叶片计算案例可以看出，风力机叶片气动性能参数的升力线模型具有较高的准确性.

3 结论

本文采用基于螺旋尾涡的升力线模型来研究风力机三维流场的空气动力学性能，研究了升力线模型中附着涡分布和诱导速度的计算，并通过对升力系数进行迭代收敛来修正环量，从而确定准确的攻角和升阻力系数，进而计算出风力机的各项气动性能参数.值得注意的是，诱导速度的求解考虑了风力机各个叶片之间的相互诱导效应，因此，相比于叶素动量理论而言，升力线理论可以极大提高风力机的三维计算能力.

基于升力线模型的理论基础，本文建立了风力机叶片的升力线数值计算模型，编译了适用于多个叶片风轮的气动性能参数分析程序，并与NREL实验室公布的基于叶素动量理论的5 MW风力机气动性能参数进行对比，结果表明，基于螺旋刚性尾涡模型的升力线算法通过结合翼型的升阻力特性可以有效计算风力机的气动性能参数，而且无需引入应用叶素动量理论时需要的大量人为修正模型即可以满足工程要求的精度.对于叶素动量理论而言，人为引入模型的不同可能会导致计算结果存在一定偏差，而升力线模型由于无需引入各种修正模型，减少了繁琐的模型引入，相对于叶素动量理论和CFD理论而言，其计算效率也会提高，同时也可以在一定程度上保证计算的稳定性.由于升力线理论考虑到了叶片展向之间的流动，其不仅适用于直叶片，对于积叠线弯曲的、三维流动更强的叶片(如后掠叶片)设计与气动性能计算同样有效，并能与叶片结构模型结合进行高效而准确的气弹耦合分析.

[1]Capuzzi M，Pirrera A，Weaver P M.A novel adaptive blade concept for large-scale wind turbines，part II：structural design and power performance [J].2014，73：15-24.

[2]Ashwill T D.Sweep-twist adaptive rotor blade：final project report [R].United States：Department of Energy，Sandia National Laboratories，2010.

[3]Pourazarm P，Modarres-Sadeghi Y，Lackner M.A pa-rametric study of coupled-mode flutter for MW-size wind turbine blades [J].Wind Energy，2016，19：497-514.

[4]石亚丽，左红梅，杨华，等.偏航角对风力机气动性能的影响 [J].农业工程学报，2015，31(16)：78-85.

(SHI Ya-li，ZUO Hong-mei，YANG Hua，et al.Aerodynamic performance of wind turbine under different yaw angles [J].Transactions of the Chinese Society of Agricultural Engineering，2015，31(16)：78-85.)

[5]Amano R S，Malloy R J.CFD analysis on aerodyna-mic design optimization of wind turbine rotor blades [J].World Academy of Science，Engineering and Technology，2009，60：71-75.

[6]张旭，刘安宇.风力机复合材料叶片流固耦合分析 [J].沈阳工业大学学报，2019，41(2)：20-24.

(ZHANG Xu，LIU An-yu.Fluid-structure coupling analysis for composite blade of wind turbine [J].Journal of Shenyang University of Technology，2019，41(2)：20-24.)

[7]李媛，康顺，樊小莉，等.后掠风力机叶片气动性能数值模拟 [J].工程热物理学报，2011，32(9)：1497-1500.

(LI Yuan，KANG Shun，FAN Xiao-li，et al.Numerical investigation of aerodynamic performance of a back-swept wind turbine blade [J].Journal of Engineering Thermophysics，2011，32(9)：1497-1500.)

[8]Jimenez A，Crespo A，Migoya E.Application of a LES technique to characterize the wake deflection of a wind turbine in yaw [J].Wind Energy，2010，13(6)：559-572.

[9]Yung C N，De Witt K J，Keith T G.Three-dimensional steady flow through a bifurcation [J].Journal of Biomechanical Engineering，1990，112(2)：189-197.

[10]李俊.大型风力发电机组整机及关键部件仿真分析与优化设计研究 [D].重庆：重庆大学，2011.

(LI Jun.Study on the simulation and optimization of the large wind turbine and its key components [D].Chongqing：Chongqing University，2011.)

[11]Jonkman J，Butterfield S，Musial W，et al.Definition of a 5-MW reference wind turbine for offshore system development [R].Golden：National Renewable Energy Laboratory，2009.

[12]王强.水平轴风力机三维空气动力学计算模型研究 [D].北京：中国科学院工程热物理研究所，2014.

(WANG Qiang.Study on 3D aerodynamic computational models of HAWT [D].Beijing：Institute of Engineering Thermophysics，Chinese Academy of Sciences，2014.)

[13]周传捷，李德源，赵世林.风力机柔性叶片的混合多体模型研究 [J].太阳能学报，2010，31(8)：1011-1017.

(ZHOU Chuan-jie，LI De-yuan，ZHAO Shi-lin.Research for the hybrid multibody systems model of the flexible blade of horizontal axis wind turbines [J].Acta Energiae Solaris Sinica，2010，31(8)：1011-1017.)

基于升力线理论的风力机气动性能分析方法*

1 风力机叶片升力线模型

1.1 诱导速度计算

1.2 迭代法求解有效攻角

2 算例分析

3 结论

Wind turbine aerodynamic performance analysis method based on lift line theory

基于升力线理论的风力机气动性能分析方法*

1 风力机叶片升力线模型

1.1 诱导速度计算

1.2 迭代法求解有效攻角

2 算例分析

3 结 论

Wind turbine aerodynamic performance analysis method based on lift line theory

3 结论