基于升力线理论的大型风力机气弹响应研究*

黄 维，李德源，吴俊霖

基于升力线理论的大型风力机气弹响应研究*

黄 维，李德源†，吴俊霖

（广东工业大学 机电工程学院，广州 510006）

针对风力机不断向大型化发展的趋势，导致结构柔度增加，气弹耦合特性和振动增强，研究了大型风力机高效精确的气弹响应分析方法。为了更准确模拟大型风力机气流沿叶片展向的三维流动现象，采用螺旋尾涡升力线模型代替传统叶素动量理论，建立了叶片气动载荷分析模型，进而结合风力机多体系统动力学模型，构建了机组的气弹耦合动力学方程和数值求解方法。以某10 MW风力机叶片为例，研究了稳态风况下不同风速的叶片气动性能，以及有效攻角、切向力等沿叶展方向的分布特点，并与采用修正叶素动量理论的气弹分析程序（HAWC）对比，结果表明，升力线理论无需引入经验修正模型即能获得叶素动量理论经修正后的分析精度。最后，通过非稳态风况下风力机的气弹响应分析，证明本文方法对大型风力机气弹耦合分析的有效性和准确性。

风力机；气动模型；升力线模型；气弹响应

0 引 言

为从风中吸取更多的能量，降低风力发电成本，风力机尺寸日益增大，风力机叶片也变得更细长，叶片结构柔度增加，气动载荷、重力、惯性载荷及操纵载荷同时随着叶片长度的增加急剧上升[1]，其带来的气弹问题也越发突出。气动弹性产生的振动将会显著影响风力机的运行寿命。因此，需要进行气动结构耦合分析来预测非稳态气动载荷与叶片弹性变形的复杂相互作用。

叶素动量（blade-element momentum, BEM）理论、涡尾迹方法和计算流体力学（computational fluid dynamics, CFD）是目前较常用的风力机气动分析方法[2]。BEM由于其简单易用、计算量小等优点被广泛使用。但由于传统的二维BEM模型忽略了气流沿叶片的展向流动，不能反映叶素单元间的相互影响，导致其不能很好地预测风轮的状态，需要人为引入修正模型，这在一定程度上降低了该理论的普适性和可靠性，尤其是对于气动流场存在三维展向流动以及非定常特性的情况[3]。CFD是一种高计算精度的叶片气动分析方法，通过求解N-S控制方程获得叶片详细的流场信息和气动载荷，可以准确地计算湍流、非定常等复杂流场结构对风力机气动性能的影响，但CFD不仅需要较高的计算成本，而且数值计算时需要与结构界面做特殊处理，这也使其在气弹分析方面应用较少[4]。与BEM和CFD这两种方法对比，涡尾迹方法从构造叶片表面的涡量出发，用不同涡量分布形式来模拟叶片与三维流场之间的相互作用。陈进格等[5]采用升力面配合自由尾迹预测法，计算尾流场对叶素的诱导作用，适用于大型柔性叶片三维气动分析。李德源等[6]采用升力线理论建立叶片的气动模型，计算5 MW叶片的气动性能并与引入修正模型的BEM对比，结果表明升力线气动模型具有较高的计算稳定性。GRECO等[7]采用一种自由涡尾迹和面元法全耦合的气动特性计算模型计算叶片在不发生严重的流量分离时的气动性能，不需BEM那样引入叶尖、叶根修正模型，不需CFD那样高的计算成本，便能得到精确的结果。

在气弹耦合分析中，赵鹰等[8]采用修正的BEM理论，结合基于有限元复合材料的叶片结构，建立了叶片的气弹耦合模型，结果表明叶片气动弹性效应会造成叶片功率输出下降，导致达不到设计功率。SAYED等[9]采用时间精确的CFD-计算结构力学（computational structural dynamics, CSD）弱耦合方法，研究了DTU10MW风力机的气弹响应，发现考虑叶片变形时风轮功率与推力提高了约1%和0.3%。钱晓航等[10]运用几何精确梁方法构建NREL5MW和IEA15MW风力机叶片，结合叶素动量理论建立了非线性叶片气弹耦合模型，采用线性分析与非线性分析方法计算了湍流风速时两个风力机组的动态响应，结果表明对于大型柔性叶片的风力机，应考虑几何非线性对风力机叶片气动弹性响应的影响，以准确评估风力机设计的安全性和稳定性。然而，目前综合考虑大型风力机整机气弹耦合特性的研究相对较少。

在保证精度的前提下，为降低模型的求解规模，提高大型风力机的气弹分析效率。本文采用“超级单元”方法将DTU10MW风力机[11]离散成若干个由运动副和力元联接的多体系统，应用Roberson-Wittenburg递推建模[12]方法，结合刚性螺旋尾涡升力线模型建立叶片非线性气弹耦合模型，将升力线模型应用于大型风力机整机的气弹响应分析，并利用双循环隐式积分更稳定地求解风力机整机的动力学微分代数方程。通过对DTU10MW风力机在稳态风工况下响应的对比分析，证明大型柔性叶片采用升力线模型的有效性和正确性，并模拟风力机在非稳态风速下的动态响应。

1 风力机多体动力学建模

1.1 整机机组离散模型

大型风力机的叶片为细长的柔性体，翼型截面形状与大小沿展向变化，非线性变形明显。为此，采用6个“超级单元”模型将叶片离散为由刚体、弹簧、阻尼器和运动副组成叶片的多刚体系统[13-14]，其中超级单元模型方法是通过将需要分析的构件分成若干单元，每个单元由4个刚体组成，相邻刚体之间通过运动副连接，并用弹簧和阻尼器约束刚体间的相对运动，采用相关计算方法获取能准确模拟构建运动的弹簧刚度系数。同样用3个超级单元离散塔架，用一个刚体简化机舱，塔顶与机舱用万向节连接，用来模拟机舱的偏航与俯仰运动，将主轴划分为两个刚体，中间用转动铰连接，其中一个刚体与机舱用转动铰连接，用来模拟主轴的转动，另一个刚体与轮毂固定连接，用一个刚体简化轮毂，三个叶片的第一个刚体分别与轮毂以转动铰连接，用来模拟叶片的变桨，从而将整机离散为71个刚体、116个自由度，由弹簧阻尼约束的刚柔混合多体系统，如图1所示。经过单元相关性数值分析，风力机系统可以用相对较少的自由度，准确地描述整机和叶片的结构动态特性，满足时域气弹响应分析要求。

图1 DTU10MW整机混合多体模型

1.2 风力机多体动力学

约束方程对时间求二阶导数，得到

联立式（3）、式（4）可将动力学方程缩并至仅包含独立自由度：

2 基于升力线叶片的气动载荷计算

大型风力机叶片弯扭耦合现象明显，结构变形使得气流沿展向流动显著，而升力线模型考虑了风力机叶片流场的三维流动，能够准确地模拟大型风力机叶片的气动性能。在升力线模型中，通过一根涡量沿展向变化的集中线涡，也称作附着涡，来模拟风力机叶片三维流场的涡量分布，用来表示叶片与流场之间的相互作用[16]。根据Helmholtz第二定理，在无黏性流体环境中，涡量不能在流体的内部终止，只能延伸至流体边界或构成环，变化的涡量会从附着涡上脱落形成尾涡[12]，如图2所示。

图2 涡方法模型示意图

从附着涡上脱落的螺旋尾涡流向叶片下游至无限远处，形成诱导速度场，本文采用刚性螺旋尾涡模型模拟风力机叶片流场[12]，尾流坐标系如图3所示。

图3 螺旋尾涡示意图

线涡对空间中的目标点位置处产生的诱导速度可通过Biot-Savart公式描述：

由Kutta-Joukowski定理可知，环量与升力系数的关系为

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

综合式（6）～ 式（10）

3 气弹耦合程序设计

整机气弹耦合方程是一个典型的非线性代数微分方程，求解方程组可以采用数值积分方法[14]。整机的气弹耦合程序的非线性代数微分方程的建立以及数值积分流程如图4所示，首先程序将会读取系统的基本数据和仿真数据，如系统的拓扑构型参数、刚体的惯性参数、铰的参数等，接着计算不随时间变化的中间常值数据，如通路阵等。完成数据的输入与前处理后，程序进入数值积分过程，将上一时刻气动模块计算的气动载荷以外力的形式加入整机的多体动力学模块中，通过求解动力学方程得到这个时刻的广义坐标、广义速度，进而计算出各刚体气动中心的速度、扭角等，并与风速一起作为这个时刻升力线模块的输入参数，通过升力线模块计算有效攻角、入流角等参数，结合升阻力系数表，计算出各刚体上的气动载荷，并作为下一时刻整机的多体动力学模块的外力，如此循环直至仿真结束时间end。本文采用MATALB中的ode15s求解器来求解系统的常微分方程组，积分精度设为1 × 10−5，输出步长为0.01 s。

图4 气弹耦合分析流程图

4 计算结果与分析

4.1 稳态风工况

为证明所建立的提气弹耦合模型的有效性与可行性，选用DTU10MW风力机为示例，在不考虑重力、主轴仰角、塔架与主轴变形、风剪切，入流风为稳定常值并垂直于风轮平面的条件下，分析叶片在来流风速为4 ～ 25 m/s时的稳态气动性能，并与采用修正叶素动量理论的HAWC程序计算结果[18]对比。HAWC程序是一个气动弹性程序，由丹麦Risø实验室开发，其结构部分采用Timoshenko梁单元构建，用于计算风力机的时域响应。各风速下风力机对应的转速、变桨距角可查阅文献[11]。

如图5所示，当风速低于12 m/s时，升力线模型计算的风轮功率略高于HAWC的结果，而推力基本一致，仅在风速11 m/s下相差较大，当风速大于12 m/s时，升力线模型计算的风轮功率小于HAWC的结果，推力基本一致，仅在高风速下略小于HAWC的推力。

图5 各风速功率、推力

为进一步了解两组模型的异同，以稳态风速10 m/s为例，研究两组模型气动参数沿叶片展向分布的异同。由图6可知两组模型有效攻角和相对速度沿叶展方向的趋势基本一致，有效攻角的大小略有偏差，在叶根与叶尖位置有效攻角相差约7°，在叶展30 ～ 80 m处，升力线模型的有效攻角约比HAWC程序结果大1.2°；由于有效攻角的偏差，导致升力系数与阻力系数也存在偏差，如图7所示，进而导致了翼型截面气动载荷的偏差，如图8所示。

综上可知，叶素动量理论通过引入诱导速度修正和叶尖损失修正来考虑尾流场对叶素的诱导影响，而升力线理论未考虑任何修正，直接计算尾流场对叶素的诱导速度，两组模型有效攻角和气动载荷沿叶展方向的趋势基本一致，表明升力线理论结合超级单元结构的气弹耦合模型是有效的、可行的。

图6 有效攻角、相对风速沿展向分布

图7 升力、阻力系数沿展向分布

图8 切向力、轴向力沿展向分布

4.2 非稳态风况气弹响应

由于大气边界层的空气运动是一种随机的湍流运动，瞬时风速随时间变化不规则。采用自回归模型模拟风速时间历程，具体过程可参阅文献[19]，以轮毂处平均风速10 m/s为例，获得轮毂处随机风速的时间历程曲线如图9所示。

图9 轮毂处风速随时间变化曲线

考虑DTU10MW风力机风轮直径达到178 m，轮毂高度达119 m，叶片在转动过程中，叶素所处高度变化大，因此引入风速随高度变化的指数规律分布加以修正：

图10为考虑随机风时风力机功率、推力的时域响应曲线。在考虑随机风时，功率与推力的变化趋势与随机风的变化趋势一致，功率最大可达11.53 MW，超过额定功率10 MW，最小4.17 MW，变化幅度为7.36 MW，功率推力最大可达1.49 MN，最小0.97 MN，变化幅度为0.52 MN，在设计大型风力机时，需要考虑功率的大幅度变化对电网电能质量的影响以及推力的大幅度变化对风力机安全运行的影响。

图10 功率、推力随时间变化曲线

图11为考虑随机风时叶片叶尖相对于轮毂中心的挥舞、摆振位移时域响应曲线。由图可知，考虑随机风时，叶尖挥舞位移最大可达7.50 m，最小4.75 m，变化幅度为2.75 m，叶尖摆振位移最大为0.441 m，最小为0.161 m，变化幅值为0.28 m，总的变化趋势与随机风一致。在设计大型风力机叶片时，需要考虑叶片叶尖的大位移可能会导致叶片与塔架的碰撞，如设计预弯叶片增大叶尖与塔架的间距，或设计自适应降载的后掠叶片减小叶尖位移。相较于叶尖挥舞位移，叶尖摆振位移随风速的快速变化振动得更剧烈。

图11 叶尖挥舞、摆振位移随时间变化曲线

图12为考虑随机风时塔顶的横向、纵向位移的时域响应曲线。可以看出，塔顶横向、纵向位移的变化趋势与随机风的变化趋势一致，随机风对塔顶位移的影响很大，在大型风力机的设计过程中，需要校核随机风况下塔架的强度，并需要开发更能反映实际风场特性的湍流风模型。

图12 塔顶横向、纵向位移随时间变化曲线

5 结 论

基于升力线理论与超级单元模型耦合建立风力机整机的气弹耦合模型，并与HAWC模型的稳态风况对比验证，然后用建立的气弹耦合模型进行随机风的数值仿真，结果表明：

（1）基于升力线理论，可以在不考虑修正的情况下计算出与修正后的叶素动量理论相近的风力机气动性能和响应。

（2）所建立的气弹模型能够很好地完成气动和结构的耦合仿真，可以很好地模拟风力机在随机风况下的气动载荷与结构位移的动态响应，该结果对风力机叶片设计具有一定的参考价值。

（3）相对于需要引入修正模型的叶素动量理论而言，升力线理论考虑了风力机三维尾流场对叶片的诱导速度，更适合分析积叠线弯曲、三维流动显著的叶片，如后掠、预弯叶片。

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

[2] 王同光, 田琳琳, 钟伟, 等. 风能利用中的空气动力学研究进展Ⅰ: 风力机气动特性[J/OL]. 空气动力学学报, 2021: 1-24. [2022-03-09].http://kns.cnki.net/kcms/detail/51.1192.TK.20211127.1024.002.html.

[3] KIM Y, MADSEN H A, APARICIO-SANCHEZ M, et al. Assessment of blade element momentum codes under varying turbulence levels by comparing with blade resolved computational fluid dynamics[J]. Renewable energy, 2020, 160: 788-802. DOI: 10.1016/j.renene.2020. 06.006.

[4] CHAVIAROPOULOS P K, NIKOLAOU I G, AGGELIS K A, et al. Viscous and aeroelastic effects on wind turbine blades. The VISCEL project. Part I: 3D Navier–Stokes rotor simulations[J]. Wind energy, 2003, 6(4): 365-385. DOI: 10.1002/we.100.

[5] 陈进格, 沈昕, 王广, 等. 基于升力面和大变形梁的风力机叶片气弹模型[J]. 工程热物理学报, 2018, 39(7): 1469-1475.

[6] 李德源, 何正举, 郭坤翔. 基于升力线理论的风力机气动性能分析方法[J]. 沈阳工业大学学报, 2020, 42(2): 148-154. DOI: 10.7688/j.issn.1000-1646.2020.02.06.

[7] GRECO L, TESTA C. Wind turbine unsteady aerodynamicsand performance by a free-wake panel method[J]. Renewable energy, 2021, 164: 444-459. DOI: 10.1016/j.renene.2020. 08.002.

[8] 赵鹰, 廖猜猜, 秦志文, 等. 基于高精度有限元模型的叶片气弹耦合分析[J]. 工程热物理学报, 2017, 38(4): 748-753.

[9] SAYED M, LUTZ T, KRÄMER E, et al. Aeroelastic analysis of 10 MW wind turbine using CFD–CSD explicit FSI-coupling approach[J]. Journal of fluids and structures, 2019, 87: 354-377. DOI: 10.1016/j.jfluidstructs. 2019.03.023.

[10] 钱晓航, 郜志腾, 王同光, 等. 百米级大柔性风电叶片非线性气弹响应分析[J/OL]. 空气动力学学报, 2021: 1-11. (2021-12-21). http://kns.cnki.net/kcms/detail/51.1192.TK.20211217.1932.002.html.

[11] BAK C, ZAHLE F, BITSCHE R, et al. The DTU 10-MW reference wind turbine[R]. Fredericia: DTU, 2013. https://orbit.dtu.dk/en/publications/the-dtu-10-mw-reference-wind-turbine.

[12] 李德源, 刘佳宇, 夏鸿建, 等. 基于涡旋升力线模型的风机后掠叶片性能优化[J]. 沈阳工业大学学报, 2021, 43(4): 384-390. DOI: 10.7688/j.issn.1000-1646.2021.04.05.

[13] 李德源, 莫文威, 严修红, 等. 基于多体模型的水平轴风力机气弹耦合分析[J]. 机械工程学报, 2014, 50(12): 140-150. DOI: 10.3901/JME.2014.12.140.

[14] 李德源, 汪显能, 莫文威, 等. 非定常条件下风力机柔性叶片气弹耦合分析[J]. 太阳能学报, 2017, 38(4): 966-975.

[15] 张乐, 章定国. 基于向后差分法求解多体系统动力学微分-代数方程组的双循环隐式积分方法[J]. 机械工程学报, 2016, 52(7): 79-87. DOI: 10.3901/JME.2016.07.079.

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

[17] JENG D R, KEITH T G, ALIAKBARKHANAFJEH A. Aerodynamic analysis of a horizontal axis wind turbine by use of helical vortex theory. Volume 1: theory[R]. Washington: NASA, 1982.

[18] VERELST D R, HANSEN M H, PIRRUNG G. Steady state comparisons HAWC2 v12.2 vs HAWCStab2 v2.12[R]. DTU Wind Energy, 2016.https://orbit.dtu.dk/en/publications/steady-state-comparisons-hawc2-v122-vs-hawcstab2-v212#.

[19] 徐磊. 基于非线性气弹耦合模型的风力机随机响应分析[D]. 广州: 广东工业大学, 2014.

Aeroelastic Response Model of Large Wind Turbine Blade Based on Lift Line Theory

HUANG Wei, LI De-yuan, WU Jun-lin

(School of Electromechanical Engineering, Guangdong University of Technology, Guangzhou 510006, China)

In view of the continuous development trend of large-scale wind turbine, which leads to the increase of structural flexibility, aeroelastic coupling characteristics and vibration, an efficient and accurate aeroelastic response analysis method of large-scale wind turbine was studied. In order to more accurately simulate the three-dimensional flow phenomenon of large-scale wind turbine airflow along the spanwise direction of the blade, the spiral tail vortex lift line model was used to replace the traditional blade element momentum theory, and the blade aerodynamic load analysis model was established. Then, combined with the multi-body system dynamic model of wind turbine, the aeroelastic coupling dynamic equation and numerical solution method of wind turbine were constructed. Taking a 10 MW wind turbine blade as an example, the aerodynamic performance of the blade with different wind speeds under steady-state wind conditions, as well as the distribution characteristics of effective angle of attack and tangential force along the blade span direction were studied. Compared with the aeroelastic analysis program (HAWC) using modified blade element momentum theory, the results showed that the modified analytical accuracy of blade element momentum theory can be obtained by lift line model without introducing empirical correction model. Finally, through the aeroelastic response analysis of wind turbine under unsteady wind conditions, the effectiveness and accuracy of this method for aeroelastic coupling analysis of large wind turbine were verified.

wind turbine; aerodynamic model; lift line model; aeroelastic response

2095-560X（2022）03-0195-08

TK89

A

10.3969/j.issn.2095-560X.2022.03.002

2022-03-01

2022-03-31

李德源，E-mail：lidey@gdut.edu.cn

黄 维（1996-），男，硕士研究生，主要从事风力机气动与结构方面的研究。

李德源（1965-），男，教授，硕士生导师，主要从事大型风力机气动与结构分析、风力机系统测试和计算机软件的开发与应用等方面的研究。