弯扭耦合对风力机叶片气弹模态特性的影响*

曲珍壮，夏鸿建，李德源，林伟豪

弯扭耦合对风力机叶片气弹模态特性的影响*

曲珍壮，夏鸿建†，李德源，林伟豪

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

弯扭耦合对风力机叶片气动弹性响应有重要影响。为研究弯扭耦合对风力机叶片气弹模态特性的影响，采用几何精确梁理论和叶素动量定理，建立风力机叶片的气弹耦合方程，在Beddoes-Leishman动态失速模型的基础上对气弹耦合方程线性化，构建特征平衡方程。由复特征值分析法计算DTU 10 MW叶片的气弹频率和气弹阻尼比。通过在截面刚度矩阵中引入耦合系数，考虑弯扭耦合对气弹模态的影响。挥舞顺桨耦合会增加挥舞模态气弹频率，降低挥舞模态气弹阻尼比。挥舞失速耦合会降低挥舞模态气弹频率，增大挥舞模态气弹阻尼比。挥舞顺桨耦合会降低风力机叶片的颤振风速。

风力机叶片；几何精确梁；弯扭耦合；颤振分析

0 引 言

弯扭耦合对风力机叶片气动弹性响应有重要影响。叶片在气动力作用下弯曲后，会改变攻角的大小，进而改变叶片的气动力，形成反馈循环。弯扭耦合可以由叶片的几何形状（几何耦合）或各向异性叶片材料（材料耦合）引起。几何耦合是叶片几何形状弯曲的结果（预弯曲、后掠或负载偏移），从而引起额外的扭转。材料耦合是叶片表皮或梁帽的纤维方向偏离主方向，从而引起弯曲和扭转的耦合。

弯扭耦合可以用来减轻风力机叶片的负载。根据LOBITZ[1]、BOTTASSO[2]等的研究，挥舞扭转顺桨耦合会使风力机叶片的疲劳载荷降低10% ～ 20%。然而弯扭耦合除了降低叶片负载外，还会改变叶片的气弹模态特性。LOBITZ[3]研究了弯扭耦合叶片与非弯扭耦合叶片发生颤振失稳临界风速，发现弯扭耦合叶片发生颤振失稳的临界风速会比非弯扭耦合叶片降低12%左右。LARWOOD[4]研究后掠叶片的颤振特性，发现后掠叶片颤振失稳的模态是二阶挥舞与一阶扭转的耦合模态。

随着风力机叶片长度和柔性的增加，叶片的非线性效应越加显著。KALLESØE[5]采用非线性二阶 Bernoulli-Euler梁理论，研究了5 MW叶片弯曲变形对气弹稳定性的影响，表明叶片的几何非线性对其气弹稳定性有明显的影响。传统采用欧拉梁结构模型[6]或模态法[7]构建叶片结构模型，难以表达柔性叶片的几何非线性变形。由于风力机叶片气弹耦合状态方程的建立比较复杂，基于频域的颤振分析研究相对较少，黄俊东等[8]采用多体动力学方法和叶素动量理论的方法建立叶片气弹耦合方程，研究了美国可再生能源实验室（National Renewable Energy Laboratory, NREL）发布的5 MW风力机气弹稳定性问题。和其他非线性模型相比，几何精确梁模型可以更为准确方便地计算叶片的变形位置以及叶片扭转的角度。通过在截面刚度矩阵中加入耦合项，研究由叶片材料引起的弯扭耦合问题。在结构动力学求解时，对结构动力学方程进行线性化处理，可以方便地与气动力线性化结合起来，进行特征值分析。

本文首先根据几何精确梁理论，结合叶素动量定理，建立风力机叶片气弹耦合方程，找出叶片在不同工况下的平衡位置。然后对气动力进行线性化，并引入Beddoes-Leishman（B-L）动态失速模型，构建特征平衡方程，采用特征值分析方法，研究风力机叶片的气弹模态特性。通过在截面刚度矩阵中引入耦合系数[9]，以丹麦技术大学（Technical University of Denmark, DTU）10 MW RWT[10]叶片为参考叶片，研究了挥舞扭转耦合对气弹模态特性的影响。

1 几何精确梁变形的描述

几何精确梁采用平截面假设，梁的位形由梁轴线的位移和相应截面的转动决定。如图1有全局坐标系（不考虑叶片预弯预扭以及截面偏心等因素）、不变形位置坐标系（考虑叶片真实几何形状）和变形位置坐标系（考虑气弹耦合后叶片的变形）。

图1 梁的转动坐标系

在变形位置，梁上任意一点的位置向量可以写为

2 结构动力学方程的建立

2.1 几何精确梁本构关系

几何精确梁理论基于Reinssner应变，根据参考文献[12-13]，定义应变为平动应变向量，应变为转动应变向量，截面等效应力、和应变、之间的关系为

应变与变形的关系为

2.2 动力学平衡方程

首先建立风力机叶片弱形式动力学平衡方程[14]，然后通过数值积分和微分方法离散求解。

对内力沿单元轴线积分，得到内力虚功为[15]

对气动力沿轴线积分，重力进行体积积分，得到单元的外力虚功

各点的速度进行体积分，得到单元的惯性力虚功

根据虚功原理，单元内力、外力和惯性力虚功满足

通过Gauss-Lobatto积分和微分求积法[14]离散平衡方程，得到基于节点虚位移表达的内力、外力和惯性力虚功[15]。

单元节点虚位移表示为

把式（12）代入式（11），并在全局坐标系下进行组装，得到等效节点内力、外力和惯性力的离散平衡方程

3 动力学方程线性化

为了构建叶片气弹线性平衡方程，形成特征值求解，需要对动力学平衡方程进行线性化处理。分为两个部分，第一部分是结构动力学方程线性化，第二部分是气动力线性化。

3.1 结构动力学线性化

分别对节点内力、外力和惯性力矩阵线性化[15]得到

将式（15）～ 式（17）代入平衡方程（14）可以得到

3.2 气动力线性化

根据叶素动量定理，作用在叶片第个截面的气动力可以写为

图2 叶片翼型截面

气动载荷的线性化表示为[8]

对状态变量线性化[6]得到

联立式（18）、式（22）和式（23），气弹线性平衡方程可以写为

4 特征值复模态计算

根据式（24）建立风力机颤振分析的特征值平衡方程

利用Matlab软件提供的函数计算特征值和特征向量，得到叶片的气弹频率、气弹阻尼比和振型，其中特征值的虚部为气弹频率，特征值的实部为气弹阻尼比。

式中：Im为复特征值的虚部；Re为复特征值的实部。

5 计算结果

5.1 旋转梁算例

旋转梁作为能否处理梁大变形标准算例，受到广泛的关注。通过与参考文献[18]对比来验证动力学方程的正确性，梁的旋转角度随时间变化的函数为

与只有转动没有变形的刚性梁相比，柔性梁的自由端偏转位移如图3所示。图中可见，计算结果与文献[18]吻合，说明本文的动力学方程是正确的。

5.2 运行工况下气弹阻尼和气弹频率

由于材料的各向异性，风力机叶片截面耦合特性比较复杂，为了单独研究挥舞扭转耦合对叶片颤振特性的影响，将丹麦技术大学[10]提供的10 MW风力机叶片（截面刚度不包含耦合项）设为参考叶片。在叶片中间部分某截面的刚度参数如式（30）。

图4显示了运行工况下，挥舞扭转耦合对一阶挥舞模态气弹频率和气弹阻尼的影响。在风力机运行工况下，三种叶片的一阶挥舞气弹频率和气弹阻尼比的变化趋势基本相同。具体为在风速4 ～ 7 m/s，一阶挥舞气弹阻尼比保持不变，在风速7 ～ 12 m/s，为上升趋势，风速12 m/s后保持平稳，此时风力机达到额定转速运行工况下，一阶挥舞气弹频率变化趋势与气弹阻尼比相似。相比于参考叶片，挥舞顺桨耦合增大了一阶挥舞气弹频率，降低了一阶挥舞气弹阻尼比，挥舞失速耦合降低了一阶挥舞气弹频率，增大了一阶挥舞气弹阻尼比。

图4 运行工况下一阶挥舞气弹频率和气弹阻尼比

挥舞扭转耦合对一阶摆振气弹频率和气弹阻尼的影响如图5所示。在运行工况下，三种叶片的一阶摆振气弹频率基本重合，一阶气弹阻尼也比较接近，表明挥舞扭转耦合对一阶摆振模态气弹特性影响较小。在运行工况下，三种叶片的一阶摆振气弹频率变化较小，变化区间为0.945 ～ 0.979 Hz。一阶摆振气弹阻尼比变化趋势为，在风速为4 ～ 11 m/s时缓慢上升，然后短暂下降后又继续保持上升，一直到风速25 m/s时达到最大值，此时气弹阻尼比大约为4%，与风速11 m/s时的气弹阻尼接近。

图5 运行工况下一阶摆振气弹频率和气弹阻尼比

5.3 定叶尖速比工况下颤振风速

为了研究弯扭耦合对风力机叶片颤振风速的影响，将风力机叶片运行工况设为定叶尖速比10，即风力机叶尖速度为来流风速的10倍，并将叶片的桨距角设为0。结果如图6。

图6 定叶尖速比工况下一阶扭转气弹频率和气弹阻尼比

挥舞顺桨耦合叶片一阶扭转气弹阻尼比在转速1.38 rad/s时出现负值，此时一阶扭转气弹频率为7.51 Hz，参考叶片大约在转速1.97 rad/s时出现负值，此时一阶扭转气弹频率为8.47 Hz，挥舞扭转失速耦合叶片没有出现负值。表明挥舞顺桨耦合会降低风力机叶片的颤振风速。这种现象的原因可能是，挥舞顺桨耦合叶片变形小，造成扭转刚度小于参考叶片，根据文献[19]，扭转刚度降低会使叶片的颤振风速降低。

6 结 论

采用非线性几何精确梁模型和叶素动量定理建立风力机叶片的气弹模型，在Beddoes-Leishman动态失速模型的基础上对气弹模型线性化，构建特征平衡方程，通过特征值分析，计算风力机叶片的气弹频率和气弹阻尼比。通过在截面刚度矩阵中引入耦合系数，研究风力机叶片挥舞扭转耦合对气弹模态特性的影响。主要结论如下：

（1）在运行工况下，挥舞扭转耦合主要影响挥舞模态，对摆振模态影响较小。挥舞顺桨耦合增大了一阶挥舞模态气弹频率，降低了一阶挥舞模态气弹阻尼比；挥舞失速耦合降低了一阶挥舞模态气弹频率，增大了一阶挥舞模态气弹阻尼比。

本文建立的气弹模型还可以用于分析风力机叶片其他的气弹问题。

[1] LOBITZ D W, VEERS P S, LAINO D J. Performance of twist-coupled blades on variable speed rotors[R]. Albuquerque: Sandia National Labs, 1999.

[2] BOTTASSO C L, CAMPAGNOLO F, CROCE A, et al. Optimization-based study of bend-twist coupled rotor blades for passive and integrated passive/active load alleviation[J]. Wind energy, 2013, 16(8): 1149-1166. DOI: 10.1002/we.1543.

[3] LOBITZ D W. Aeroelastic stability predictions for a MW-sized blade[J]. Wind energy, 2004, 7(3): 211-224. DOI: 10.1002/we.120.

[4] LARWOOD S. Flutter of variations on a 5 MW swept wind turbine blade[J]. Journal of solar energy engineering, 2016, 138(2): 024504. DOI: 10.1115/1.4032475.

[5] KALLESØE B S. Effect of steady deflections on the aeroelastic stability of a turbine blade[J]. Wind energy, 2011, 14(2): 209-224. DOI: 10.1002/we.413.

[6] HANSEN M H. Aeroelastic properties of backward swept blades[C]//49th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition. Orlando: AIAA, 2011. DOI: 10.2514/6.2011-260.

[7] 刘雄, 李钢强, 陈严, 等. 水平轴风力机叶片动态响应分析[J]. 机械工程学报, 2010, 46(12): 128-134, 141. DOI: 10.3901/JME.2010.12.128.

[8] 黄俊东, 夏鸿建, 李德源, 等. 大型风力机柔性叶片非线性气弹模态分析[J]. 机械工程学报, 2020, 56(14): 180-187. DOI: 10.3901/JME.2020.14.180.

[9] LOBITZ D W, VEERS P S. Aeroelastic behavior of twist-coupled HAWT blades[C]//1998 ASME Energy Symposium. Reno: ASME, 1998: 75-83. DOI: 10.2514/6.1998-29.

[10] BAK C, ZAHLE F, BITSCHE R, et al. Description of the DTU 10 MW reference wind turbine[R]. Roskilde: DTU, 2013.

[11] ZUPAN E, SAJE M, ZUPAN D. The quaternion-based three-dimensional beam theory[J]. Computer methods in applied mechanics and engineering, 2009, 198(49/52): 3944-3956. DOI: 10.1016/j.cma.2009.09.002.

[12] REISSNER E. On fifinite deformation of space-curved beams[J]. Z. angew. Math Phys., 1981, 32: 734-744. https://doi.org/10.1007/BF00946983.

[13] SIMO J C. A finite strain beam formulation. The three-dimensional dynamic problem. Part I[J]. Computer methods in applied mechanics and engineering, 1985, 49(1): 55-70. DOI: 10.1016/0045-7825(85)90050-7.

[14] ZHANG R, ZHONG H Z. Weak form quadrature element analysis of planar slender beams based on geometrically exact beam theory[J]. Archive of applied mechanics, 2013, 83(9): 1309-1325. DOI: 10.1007/s00419-013-0748-3.

[15] 吕品, 廖明夫, 徐阳, 等. 基于几何精确梁理论的风力机叶片单元[J]. 太阳能学报, 2015, 36(10): 2422-2428. DOI: 10.3969/j.issn.0254-0096.2015.10.017.

[16] ZHONG H Z, YU T. Flexural vibration analysis of an eccentric annular Mindlin plate[J]. Archive of applied mechanics, 2007, 77(4): 185-195. DOI: 10.1007/s00419-006-0083-z.

[17] HANSEN M H, GAUNAA M, MADSEN H A. A Beddoes-Leishman type dynamic stall model in state-space and indicial formulations[R]. Denmark: Risφ National Laboratory, 2004.

[18] FAN W, ZHU W D. An accurate singularity-free formulation of a three-dimensional curved Euler–Bernoulli beam for flexible multibody dynamic analysis[J]. Journal of vibration and acoustics, 2016, 138(5): 051001. DOI: 10.1115/1.4033269.

[19] ZHOU X D, HUANG K F, LI Z. Effects of bend-twist coupling on flutter limits of composite wind turbine blades[J]. Composite structures, 2018, 192: 317-326. DOI: 10.1016/j.compstruct.2018.02.071.

Effect of Bending-Torsion Coupling on Aeroelastic Modal Characteristics of Wind Turbine Blades

QU Zhenzhuang, XIA Hongjian†, LI Deyuan, LIN Weihao

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

The bending-torsion coupling has an important influence on the aeroelastic response of wind turbine blades. To study the effect of bending-torsion coupling on the aeroelastic modal characteristics of wind turbine blades, the geometrically exact beam theory and blade element momentum theory are used to establish the aeroelastic coupling equation for wind turbine blades. Based on the Beddoes-Leishman dynamic stall model, the aeroelastic coupling equation is linearized and the characteristic equilibrium equation is constructed. The DTU 10 MW blade’s aeroelastic frequencies and damping ratio are calculated using complex eigenvalue analysis. The influence of bending-torsion coupling on the aeroelastic modes is considered by introducing coupling coefficients into the section stiffness matrix. The flap-twist to feather coupling increases the flapwise aeroelastic frequency and decreases the flapwise aeroelastic damping ratio, while the flap-twist to stall coupling decreases the flapwise aeroelastic frequency and increases the flapwise aeroelastic damping ratio. The flap-twist to feather can reduce the flutter speed of wind turbine blades.

wind turbine blade; geometrically exact beam; bending-torsion coupling; flutter analysis

2095-560X（2023）04-0333-07

TK83；O313.7

A

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

2023-03-08

2023-04-12

国家自然科学基金项目（51776044，51105079）；广东省自然科学基金项目（2020A1515010844）

夏鸿建，E-mail：hjxia@gdut.edu.cn

曲珍壮, 夏鸿建, 李德源, 等. 弯扭耦合对风力机叶片气弹模态特性的影响[J]. 新能源进展, 2023, 11(4): 333-339.

：QU Zhenzhuang, XIA Hongjian, LI Deyuan, et al. Effect of bending-torsion coupling on aeroelastic modal characteristics of wind turbine blades[J]. Advances in new and renewable energy, 2023, 11(4): 333-339.

曲珍壮（1996-），男，硕士研究生，主要从事基于几何精确梁模型的风力机叶片研究。

夏鸿建（1978-），男，博士，副教授，硕士生导师，主要从事风力机气动与结构等研究。