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

第 11 卷 第 4 期 新 能 源 进 展 Vol.11 No.4 2023 年 8 月 ADVANCES IN NEW AND RENEWABLE ENERGY Aug.2023 *收稿日期：2023-03-08 修订日期：2023-04-12 基金项目：国家自然科学基金项目（51776044，51105079）；广东省自然科学基金项目（2020A1515010844）通信作者：夏鸿建，E-mail： 引用本文：曲珍壮,夏鸿建,李德源,等.弯扭耦合对风力机叶片气弹模态特性的影响J.新能源进展,2023,11(4):333-339.Citation：QU Zhenzhuang,XIA Hongjian,LI Deyuan,et al.Effect of bending-torsion coupling on aeroelastic modal characteristics of wind turbine bladesJ.Advances in new and renewable energy,2023,11(4):333-339.文章编号：2095-560X（2023）04-0333-07 弯扭耦合对风力机叶片气弹模态特性的影响*曲珍壮，夏鸿建，李德源，林伟豪（广东工业大学 机电工程学院，广州 510006）摘 要：弯扭耦合对风力机叶片气动弹性响应有重要影响。为研究弯扭耦合对风力机叶片气弹模态特性的影响，采用几何精确梁理论和叶素动量定理，建立风力机叶片的气弹耦合方程，在 Beddoes-Leishman 动态失速模型的基础上对气弹耦合方程线性化，构建特征平衡方程。由复特征值分析法计算 DTU 10 MW 叶片的气弹频率和气弹阻尼比。通过在截面刚度矩阵中引入耦合系数，考虑弯扭耦合对气弹模态的影响。挥舞顺桨耦合会增加挥舞模态气弹频率，降低挥舞模态气弹阻尼比。挥舞失速耦合会降低挥舞模态气弹频率，增大挥舞模态气弹阻尼比。挥舞顺桨耦合会降低风力机叶片的颤振风速。关键词：风力机叶片；几何精确梁；弯扭耦合；颤振分析 中图分类号：TK83；O313.7 文献标志码：A DOI：10.3969/j.issn.2095-560X.2023.04.005 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)Abstract: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 blades 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.Keywords:wind turbine blade;geometrically exact beam;bending-torsion coupling;flutter analysis 0 引 言 弯扭耦合对风力机叶片气动弹性响应有重要影响。叶片在气动力作用下弯曲后，会改变攻角的大小，进而改变叶片的气动力，形成反馈循环。弯扭耦合可以由叶片的几何形状（几何耦合）或各向异性叶片材料（材料耦合）引起。几何耦合是叶片几何形状弯曲的结果（预弯曲、后掠或负载偏移），从而引起额外的扭转。材料耦合是叶片表皮或梁帽的纤维方向偏离主方向，从而引起弯曲和扭转的耦合。弯扭耦合可以用来减轻风力机叶片的负载。根据 LOBITZ1、BOTTASSO2等的研究，挥舞扭转顺桨耦合会使风力机叶片的疲劳载荷降低10%20%。然而弯扭耦合除了降低叶片负载外，还会改变叶片334 新 能 源 进 展 第 11 卷 的气弹模态特性。LOBITZ3研究了弯扭耦合叶片与非弯扭耦合叶片发生颤振失稳临界风速，发现弯扭耦合叶片发生颤振失稳的临界风速会比非弯扭耦合叶片降低 12%左右。LARWOOD4研究后掠叶片的颤振特性，发现后掠叶片颤振失稳的模态是二阶挥舞与一阶扭转的耦合模态。随着风力机叶片长度和柔性的增加，叶片的非线性效应越加显著。KALLESE5采用非线性二阶 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 RWT10叶片为参考叶片，研究了挥舞扭转耦合对气弹模态特性的影响。1 几何精确梁变形的描述 几何精确梁采用平截面假设，梁的位形由梁轴线的位移和相应截面的转动决定。如图 1 有全局坐标系（不考虑叶片预弯预扭以及截面偏心等因素）、不变形位置坐标系（考虑叶片真实几何形状）和变形位置坐标系（考虑气弹耦合后叶片的变形）。图 1 梁的转动坐标系 Fig.1 The rotating coordinate system of the beam 在变形位置，梁上任意一点 p 的位置向量可以写为,ppt z x yt zx yuuR (1)式中：pu为剪切中心在全局坐标系下的坐标；z 为截面在梁轴线上的坐标；为截面转动矩阵；pR为截面坐标系下的位置向量。叶片坐标系的建立为，在叶片根部以各个截面的剪心为原点，建立局部坐标系，将局部坐标系平移和旋转后，得到叶片的真实几何形状。三维旋转矩阵 有 9 个元素，包含 6 个约束关系，难以直接作为姿态变量，欧拉四元数是描述截面转动非奇异的最小参数表达，本文采用欧拉四元数 q描述截面的转动11。T00123,qq q q qqq (2)式中：四元数之间满足222201231qqqq；T（）表示向量的转置。四元数 q和转动矩阵 的关系为 2T0021q2q2I+q+qq (3)2 结构动力学方程的建立 2.1 几何精确梁本构关系 几何精确梁理论基于 Reinssner 应变，根据参考文献12-13，定义应变 为平动应变向量，应变 K为转动应变向量，截面等效应力 N、M 和应变、K之间的关系为 NSMK (4)第 4 期 曲珍壮等：弯扭耦合对风力机叶片气弹模态特性的影响 335 式中：S 为 6 6 截面刚度矩阵。不考虑耦合时，S矩阵只包含对角线元素。通过在截面刚度矩阵 S 中设置耦合项46K可以引入挥舞扭转耦合。464466KK K (5)式中：为挥舞扭转耦合系数9；44K和66K分别为挥舞刚度和扭转刚度。应变与变形的关系为 T0u (6)T*12300,2K KKqqK (7)式中：0和0K为梁的初始位置和曲率参数；()为在参考构形中对截面在 z 轴上坐标的导数；为四元数乘法。2.2 动力学平衡方程 首先建立风力机叶片弱形式动力学平衡方程14，然后通过数值积分和微分方法离散求解。对内力沿单元轴线积分，得到内力虚功为15 eeTTint0dlZWNK M (8)对气动力沿轴线积分，重力进行体积积分，得到单元的外力虚功 eeTTTextaaaap0ddlVZVWu fmug (9)式中：au为轴线气动节点的虚位移；a 为轴线节点气动点的虚转动；af为轴线节点的气动力向量；am为轴线节点的气动弯矩向量；g 为重力加速度向量；为叶片体密度；el为单元长度；下标 e 表示叶片单元。各点的速度进行体积分，得到单元的惯性力虚功 eTinertdppVVWuu (10)根据虚功原理，单元内力、外力和惯性力虚功满足 eeeintextinert0WWW (11)通过 Gauss-Lobatto 积分和微分求积法14离散平衡方程，得到基于节点虚位移表达的内力、外力和惯性力虚功15。eintinteextexteinertinertWdRWdRWdReTeTeT (12)单元节点 k 虚位移表示为 kkkkxkykzkxyzd=(13)单元位移e12 Kdddd。把式（12）代入式（11），并在全局坐标系下进行组装，得到等效节点内力、外力和惯性力的离散平衡方程 intextinert0RRR (14)3 动力学方程线性化 为了构建叶片气弹线性平衡方程，形成特征值求解，需要对动力学平衡方程进行线性化处理。分为两个部分，第一部分是结构动力学方程线性化，第二部分是气动力线性化。3.1 结构动力学线性化 分别对节点内力、外力和惯性力矩阵线性化15得到 intintRKd (15)TaextextRKdBFB d (16)inertinertinertinertRKdCdMd (17)将式（15）式（17）代入平衡方程（14）可以得到 Ta0M dC dK dBFB d (18)式中：intextinertKKKK；inertMM；inertCC；为弱形式求积长度系数16；为Gauss-Lobatto积分系数；B为微分求积定位矩阵。3.2 气动力线性化 根据叶素动量定理，作用在叶片第k个截面的气动力可以写为 a2a12kkkkc WfC (19)式中：c为截面弦长；a为空气密度；W为相对来流风速；TLDM,CCCC，其中LC、DC、MC分别为升力系数、阻力系数和扭矩系数。336 新 能 源 进 展 第 11 卷 考虑叶片振动和入流角的影响，相对来流风速W和攻角可以写为 WVnxRy (20)1VnxtgRy (21)式中：Vn为垂直旋转平面的来流风速；为入流角；为桨距角；为叶片扭角；y 和x 分别为摆振和挥舞方向振动速度；R为翼型截面剪心到轮毂中心的距离；为叶片旋转速度，见图2。图 2 叶片翼型截面 Fig.2 Blade airfoil section 当来流风速和叶片旋转速度均已确定时，叶片在某一瞬时受到的气动力是振动速度的函数。为研究在某一工况下气动力对叶片振动特性的影响，首先将af进行坐标转换，得到全局坐标下的气动力aF，然后采用变分原理对气动力进行线性化处理。在线性化的过程中，为了考虑涡旋脱落和尾缘流动分离对气动特性的影响，引入基于状态变量Ta1234 x x x xx的B-L动态失速模型17，修正叶片各翼型截面处的升力系数及其斜率。气动载荷的线性化表示为8 aaaaafFMdCdKdKx (22)式中：aK、aM、aC分别是气动部分线性化得到的气动质量阵、刚度阵和阻尼阵；fK为状态变量系数矩阵。对状态变量线性化6得到 aadsssxAxMdCdKd (23)式中：sM、sK、sC分别为状态变量线性化引入的质量阵、刚度阵和阻尼阵；dA为系数矩阵。联立式（18）、式（22）和式（23），气弹线性平衡方程可以写为 aaaafaadsss0MMdCCdKKdKxxA xMdCdKd (24)4 特征值复模态计算 根据式（24）建立风力机颤振分析的特征值平衡方程 0yKM (25)设a(,)Tydd x，其解的形式为etyY；其中 aafssd 0 0 IKKKCCKKCA (26)as 0 00 00 IMMMMI (27)利用Matlab软件提供的eig函数计算特征值和特征向量，得到叶片的气弹频率、气弹阻尼比和振型，其中特征值的虚部为气弹频率，特征值的实部为气弹阻尼比。ImRe,22iiiiiff (28)式中：Im为复特征值的虚部；Re为复特征值的实部。5 计算结果 5.1 旋转梁算例 旋转梁作为能否处理梁大变形标准算例，受到广泛的关注。通过与参考文献18对比来验证动力学方程的正确性，梁的旋转角度随时间变化的函数为 22ssssssss2cos1 022()2TtttTTTtTttT(29)式中：Ts=15 s；s4 rad/s。与只有转动没有变形的刚性梁相比，柔性梁的自由端偏转位移如图3所示。图中可见，计算结果与文献18吻合，说明本文的动力学方程是正确的。第4期 曲珍壮等：弯扭耦合对风力机叶片气弹模态特性的影响 337 05101520-0.6-0.5-0.4-0.3-0.2-0.10.00.1 本文数据 参考文献 18 数据梁自由端横向位移/m时间/s(a)05101520-0.025-0.020-0.015-0.010-0.0050.0000.005 本文数据 参考文献 18 数据 梁自由端竖向位移/m时间/s(b)图3 梁自由端横向偏转（a）和竖向偏转（b）Fig.3 Longitudinal(a)and transverse(b)deflections at the tip of the rotating beam 5.2 运行工况下气弹阻尼和气弹频率 由于材料的各向异性，风力机叶片截面耦合特性比较复杂，为了单独研究挥舞扭转耦合对叶片颤振特性的影响，将丹麦技术大学10提供的10 MW风力机叶片（截面刚度不包含耦合项）设为参考叶片。在叶片中间部分某截面的刚度参数如式（30）。9cs0.514000000.37800007.57000103.3800sym.11.500.661K(30)截面耦合系数0.1 为挥舞顺桨耦合叶片，挥舞顺桨耦合会减小攻角。截面耦合系数0.1为挥舞失速耦合叶片，挥舞失速耦合会增大攻角。在运行工况下考察挥舞扭转耦合对叶片气弹模态特性的影响。图4显示了运行工况下，挥舞扭转耦合对一阶挥舞模态气弹频率和气弹阻尼的影响。在风力机运行工况下，三种叶片的一阶挥舞气弹频率和气弹阻尼比的变化趋势基本相同。具体为在风速4 7 m/s，一阶挥舞气弹阻尼比保持不变，在风速7 12 m/s，为上升趋势，风速12 m/s后保持平稳，此时风力机达到额定转速运行工况下，一阶挥舞气弹频率变化趋势与气弹阻尼比相似。相比于参考叶片，挥舞顺桨耦合增大了一阶挥舞气弹频率，降低了一阶挥舞气弹阻尼比，挥舞失速耦合降低了一阶挥舞气弹频率，增大了一阶挥舞气弹阻尼比。0.40.60.81.01.2510152025406080100 参考叶片 挥舞扭转顺桨 挥舞扭转失速气弹频率/Hz气弹阻尼比/%风速/(m/s)图4 运行工况下一阶挥舞气弹频率和气弹阻尼比 Fig.4 First flapwise aeroelastic frequency and damping under the operational wind speed 挥舞扭转耦合对一阶摆振气弹频率和气弹阻尼的影响如图5所示。在运行工况下，三种叶片的一阶摆振气弹频率基本重合，一阶气弹阻尼也比较接近，表明挥舞扭转耦合对一阶摆振模态气弹特性影响较小。在运行工况下，三种叶片的一阶摆振气弹频率变化较小，变化区间为0.945 0.979 Hz。一阶摆振气弹阻尼比变化趋势为，在风速为4 11 m/s时缓慢上升，然后短暂下降后又继续保持上升，一直到风速25 m/s时达到最大值，此时气弹阻尼比大约为4%，与风速11 m/s时的气弹阻尼接近。0.40.60.81.01.25101520250246 参考叶片 挥舞扭转顺浆 挥舞扭转失速气弹频率/Hz气弹阻尼比/%风速/(m/s)图5 运行工况下一阶摆振气弹频率和气弹阻尼比 Fig.5 First edgewise aeroelastic frequency and damping under the operational wind speed 338 新 能 源 进 展 第11卷 5.3 定叶尖速比工况下颤振风速 为了研究弯扭耦合对风力机叶片颤振风速的影响，将风力机叶片运行工况设为定叶尖速比10，即风力机叶尖速度为来流风速的10倍，并将叶片的桨距角设为0。结果如图6。6810120.60.81.01.21.41.61.82.02.2-8-404 参考叶片 挥舞扭转顺浆 挥舞扭转失速气弹频率/Hz气弹阻尼比/%转速/(rad/s)图6 定叶尖速比工况下一阶扭转气弹频率和气弹阻尼比 Fig.6 First torsion aeroelastic frequency and aeroelastic damping under constant tip speed ratio 挥舞顺桨耦合叶片一阶扭转气弹阻尼比在转速1.38 rad/s时出现负值，此时一阶扭转气弹频率为7.51 Hz，参考叶片大约在转速1.97 rad/s时出现负值，此时一阶扭转气弹频率为8.47 Hz，挥舞扭转失速耦合叶片没有出现负值。表明挥舞顺桨耦合会降低风力机叶片的颤振风速。这种现象的原因可能是，挥舞顺桨耦合叶片变形小，造成扭转刚度小于参考叶片，根据文献19，扭转刚度降低会使叶片的颤振风速降低。6 结 论 采用非线性几何精确梁模型和叶素动量定理建立风力机叶片的气弹模型，在Beddoes-Leishman动态失速模型的基础上对气弹模型线性化，构建特征平衡方程，通过特征值分析，计算风力机叶片的气弹频率和气弹阻尼比。通过在截面刚度矩阵中引入耦合系数，研究风力机叶片挥舞扭转耦合对气弹模态特性的影响。主要结论如下：（1）在运行工况下，挥舞扭转耦合主要影响挥舞模态，对摆振模态影响较小。挥舞顺桨耦合增大了一阶挥舞模态气弹频率，降低了一阶挥舞模态气弹阻尼比；挥舞失速耦合降低了一阶挥舞模态气弹频率，增大了一阶挥舞模态气弹阻尼比。（2）在定叶尖速比工况下，挥舞顺桨耦合风力机叶片一阶扭转气弹频率在风速1.38 rad/s出现负值，参考叶片一阶扭转气弹频率在风速1.97 rad/s出现负值，挥舞失速耦合一阶扭转气弹频率没有出现负值。截面耦合系数0.1 时，挥舞顺桨耦合会降低风力机叶片的颤振风速，约降低了29%。本文建立的气弹模型还可以用于分析风力机叶片其他的气弹问题。参考文献：1 LOBITZ D W,VEERS P S,LAINO D J.Performance of twist-coupled blades on variable speed rotorsR.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 alleviationJ.Wind energy,2013,16(8):1149-1166.DOI:10.1002/we.1543.3 LOBITZ D W.Aeroelastic stability predictions for a MW-sized bladeJ.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 bladeJ.Journal of solar energy engineering,2016,138(2):024504.DOI:10.1115/1.4032475.5 KALLESE B S.Effect of steady deflections on the aeroelastic stability of a turbine bladeJ.Wind energy,2011,14(2):209-224.DOI:10.1002/we.413.6 HANSEN M H.Aeroelastic properties of backward swept bladesC/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 bladesC/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 turbineR.Roskilde:DTU,2013.11 ZUPAN E,SAJE M,ZUPAN D.The quaternion-based three-dimensional beam theoryJ.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 beamsJ.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 IJ.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 theoryJ.Archive of applied mechanics,2013,83(9):1309-1325.DOI:10.1007/s00419-013-0748-3.第4期 曲珍壮等：弯扭耦合对风力机叶片气弹模态特性的影响 339 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 plateJ.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 formulationsR.Denmark:Ris National Laboratory,2004.18 FAN W,ZHU W D.An accurate singularity-free formulation of a three-dimensional curved EulerBernoulli beam for flexible multibody dynamic analysisJ.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 bladesJ.Composite structures,2018,192:317-326.DOI:10.1016/pstruct.2018.02.071.作者简介：曲珍壮（1996-），男，硕士研究生，主要从事基于几何精确梁模型的风力机叶片研究。夏鸿建（1978-），男，博士，副教授，硕士生导师，主要从事风力机气动与结构等研究。