大型风力机的几何非线性动态响应研究

曹九发， 王同光， 王　枭

(1. 扬州大学 水利与能源动力工程学院，江苏 扬州　225127; 2.南京航空航天大学 江苏省风力机设计高技术研究重点实验室,南京　210016)



大型风力机的几何非线性动态响应研究

曹九发1， 王同光2， 王枭2

(1. 扬州大学 水利与能源动力工程学院，江苏 扬州225127; 2.南京航空航天大学 江苏省风力机设计高技术研究重点实验室,南京210016)

摘要：随着风力机大型化、柔性化发展，风力机叶片的变形越来越大，几何非线性效应也越来越明显。采用自由涡尾迹方法提高风力机气动载荷计算精度，应用有限元方法建立了风力机叶片结构梁模型。应用伽辽金方法建立叶片的几何非线性刚度矩阵，同时构建了叶片的非线性动力学方程。采用Newmark直接积分法进行时间步长推进，并采用修正的Newton-Raphson方法进行增量迭代计算，实现了大型风力机叶片气动与非线性结构耦合计算模拟。最后，通过美国可再生能源实验室的 Phase VI风力机实验数据验证了计算模型的准确性和可靠性。同时以大型风力机NH1500为算例，计算了叶片的线性与非线性动态响应，分析了叶片几何非线性对动态响应和气动性能的影响。这对于提高大型风力机设计水平和各个部件载荷计算准确度具有重要意义。

关键词：风力机；自由涡尾迹方法；几何非线性；Newton-Raphson方法；动态响应

随着风力机大型化、柔性化发展，其叶片在气动力的作用下，挥舞和摆振方向变形进一步增大，这不仅会影响风轮气动性能，而且会改变叶片表面载荷分布，影响结构安全性，在一定情况下，甚至会引起结构破坏。因此，能够准确地计算出风力机的动态响应，对风力机的发展起着至关重要的作用。在叶片大尺度变形的情况下，将其简化为线性梁模型的方法不再适用，而采用几何非线性[1-2]假设能更准确地计算风力机动态响应。

叶素动量理论简单快捷，但准确度不高。CFD方法虽然能够较为准确的模拟叶片上的复杂流动，但计算效率低，消耗资源巨大。涡尾迹方法具有旋涡特性，计算成本较CFD更小，准确度较叶素动量理论更高，而且能够较好的模拟风力机尾流场。预定涡尾迹方法是根据大量的尾流场实验数据，建立尾迹形状的经验描述函数。相比之下，自由涡尾迹方法[3-5]不需要涡元位置的先验数据，而且尾迹涡元允许在当地速度场的影响下自由变形。考虑到叶片变形会导致风场及叶片气动外形变化，自由涡尾迹方法更适用于风力机气动结构耦合计算。

将风力机叶片简化为欧拉梁模型，其响应计算方法可分为工程梁方法，有限元方法[6-7]和多体动力学方法[8-10]。有限元方法求解动力学方程时，可采用振型叠加法[11]和直接积分法。振型叠加法对动力学方程进行模态分析，能够有效地减少结构模型的自由度数和矩阵规模，提高计算速度，但只适用于线性计算。直接积分法如Newmark方法，计算量比振型叠加法较大，但在节点数较少时，二者差别不明显，而且直接积分法适用于非线性计算。本文在Newmark直接积分法的基础上采用修正的Newton-Raphson增量计算法[12-13]，用来计算非线性动态响应。

本文采用自由涡尾迹方法计算了风力机的气动载荷，并与实验值进行对比，验证了气动计算方法的可靠性与准确性。然后，通过自由涡尾迹方法与有限元方法耦合计算，实现了大型风力机NH1500的动态响应求解。

1自由涡尾迹方法

1.1尾流场描述

自由涡尾迹方法中，对风力机流场作不可压和位流假设，气动模型可以简化为来流、叶片附着涡线和自由涡面的总和。如图1所示，叶片附着涡线位于1/4弦线处，并采用“arc-cosine”法离散，每段附着涡线代替每段叶素，叶素控制点位于3/4弦线处，从而叶片被模拟成一个Weissinger-L升力面模型。自由涡面是由叶片尾缘拖出涡线形成，可分为尾随涡线和脱体涡线，分别模拟附着环量在空间和时间上的变化。尾随涡强度定义为相邻叶素的附着环量之差：

(1)

式中：i=1,2,…,NE，体现不同叶素，Γt为尾随涡环量强度，Γb为叶素附着涡环量强度。

脱体涡强度为相邻方位角上叶素附着环量之差，则第j个方位角的脱体涡强度为：

(2)

式中：j=1,2,…,NT，体现不同方位角，Γs为脱体涡环量强度。

图1　风力机尾迹离散描述示意图Fig.1 Wake discrete description of wind turbine

1.2涡线控制方程

本文模型中每根涡线均在远场截断，流场中涡线随当地流速移动会自由卷起。涡线的偏微分控制方程可写为：

(3)

推导尾迹控制方程中的对时间步微分方程的差分，令Δψ=Δζ，可得到控制方程的离散格式：

预估步：

(4)

校正步：

2几何非线性动力学方程

风力机向大型化、柔性化方向发展，叶片变形进一步增大，尽管应变很小，没有超过弹性极限，但是位移较大，为了更准确的计算风力机动态响应，动力学方程应该建立在变形后的位形上。因此，本文考虑的几何非线性问题源于叶片的大位移、小应变问题。

风力机为旋转机械，动力学方程建立在随叶片转动的叶根坐标系上[14]。计算叶片变形时，将其根部固定，建立悬臂梁的动力学方程；进行非线性求解时，则在每个方位角时，采用完全拉格朗日格式(T.L.)进行增量分析。

2.1结构动力学方程

求解叶片动态响应，需要构造叶片的动力学方程：

(6)

式中：M为质量矩阵，C为阻尼矩阵，K为刚度矩阵，x为叶片变形量。

将风力机叶片简化为欧拉梁模型[15]，采用有限元方法进行线性与非线性分析[16]，构造动力学方程所需的质量矩阵，刚度矩阵和阻尼矩阵。采用伽辽金方法计算单元矩阵，其在力学范围内体现为虚功原理。质量矩阵在线性与非线性分析时，是没有区别的。线性与非线性分析的主要区别在于有没有考虑应变的二次项，即非线性项，其体现在刚度矩阵的求解过程中。阻尼矩阵表示为质量矩阵与刚度矩阵的线性组合，即C=c1M+c2K，c1和c2为比例常数，可通过实验确定。

忽略叶片扭转，梁模型节点有五个自由度：展向位移u，挥舞方向位移v和转角α，摆振方向位移w和转角β。由此，单元任一点位移表示为[u,v,w]T=N(x)δ。其中，N(x)为形函数矩阵，δ=[ui,vi,αi,wi,βi,uj,vj,αj,wj,βj]T为单元节点位移列阵。

单元动能为：

(7)

2.1.1线性分析

线性分析时，单元应变能为：

(8)式中：u′)为u对x方向一次导数，表示x方向应变；v″,w″为v,w对x方向二次导数，表示曲率；E为杨氏弹性模量。

外力做功为：

W=Fδ

(9)

由哈密顿原理得：

(10)

式中：K0为线性单元刚度矩阵。

得到单元质量矩阵和单元刚度矩阵之后，通过有限元方法，将各单元矩阵组合为一致质量矩阵和一致刚度矩阵，并求得一致阻尼矩阵，建立动力学方程。然后，采用Newmark方法[16]进行时间步长推进。

2.1.2非线性分析

非线性分析时，单元应变能为：

(11)

由哈密顿原理得：

(12)

式中：KN为非线性刚度矩阵，其中包含了与当前位形相关的元素，K0+KN为全量刚度矩阵。

由于每次位移和转角变化，都会导致非线性刚度矩阵变化，因此在时域分析时，宜采用增量计算的方法，需求解切线刚度矩阵：

(13)

式中：(K0+KT)为切线刚度矩阵。

时间步长推进时，在Newmark方法的基础上采用Newton-Raphson方法，从t时刻向t+Δt时刻递推公式如下：

(14)

式中变量可参考文献[13,16]。由于方程左端的切线刚度矩阵是非线性的，所以每一个时间步内都需要迭代求解。

2.2气动结构耦合计算实现

在一个时间步内，首先应用自由涡尾迹方法进行气动计算，然后将气动载荷转移到结构模型中，再考虑叶片所受重力与离心力，作为叶片受到的激励。然后计算叶片的刚度矩阵，质量矩阵和阻尼矩阵，建立动力学方程，求解叶片响应，并反馈到气动计算中，重新构建气动外形并进行下一个时间步的计算。

线性计算时，刚度矩阵在整个时间流程中是不变的。非线性计算时，每次迭代都会更新叶片节点位移与转角，全量刚度矩阵和切线刚度矩阵也会不断更新，直到计算结果达到收敛条件。由于每次迭代都需计算切线刚度矩阵，计算量较大，因此采用修正的Newton-Raphson方法，即每个时间步中，全量刚度矩阵每次迭代都进行更新，而切线刚度矩阵只在第一个迭代步更新。这样虽然可能需要更多的迭代次数才能收敛，但是总的计算量减小了，而且此方法更稳定，不易发散。非线性动态响应计算流程图如图2所示。

图2　耦合计算流程图Fig.2 Coupling calculation flow chart

3计算结果及分析

3.1计算方法验证

通过美国可再生能源实验室进行的NREL Phase VI风力机实验数据[17]，验证本文计算模型的可靠性与准确性。NREL Phase VI风力机为2叶片风力机，风轮直径为10.058 m，轮毂中心高度为12.192 m，额定转速为72 r/min，翼型采用S809翼型，额定功率为19.8 kW，属于失速型功率控制。

图3是各个稳定风速工况下叶片上翼型的法向系数Cn和切向系数Ct，同时给出它们随着展向方向的变化曲线图。从图3可以看出，叶片展向位置各个截面的计算值与实验值具有较好的吻合度，特别是在低风速时。靠近叶尖处95%截面，整个风速区域Cn和Ct与实验值的吻合度比叶片中部和叶根要好，这是由于靠近叶尖处的诱导速度主要是叶尖涡在起作用，而且靠近叶尖处桨距角小，是整个叶片展向入流迎角最小的位置，失速程度较弱。但是，在高风速段，由于叶根处总是容易处于深失速的状态，失速程度较严重，因此，相对于叶片其他部位和低风速，靠近叶根位置预测计算精确度较弱。从整体上看，本文计算模型具有一定可靠性和准确性。

3.2大型风力机叶片非线性动态响应计算

首先计算了大型风力机叶片NH1500不同叶尖速比的响应，并与实验结果进行对比分析，如图4所示。NH1500风力机是变桨变速风力机，具体的风力机性能和几何参数见表1。

图3　稳态风速下叶片展向的Cn和CtFig.3 Variation of normal and tangential coefficient at different span with steady wind speed

本文计算结果与实验值比较吻合，但比实验值偏大。主要是因为本文气动计算方法只在涡核处考虑了黏性，而没有考虑全流场的黏性，而且实验采用缩比模型，雷诺数偏小，导致功率系数偏小。

然后，分别采用线性分析与非线性分析计算了风速10 m/s时NH1500动态响应，并将叶尖挥舞、摆振位移以及风力机气动性能进行对比分析。

表1　 NH1500风力机设计参数

图4　NH1500计算验证Fig.4 Validation of NH1500

对于叶尖挥舞位移，非线性与线性计算结果有明显差别。通过图5和表2可以看出，无论是平衡位置还是振幅，非线性计算结果都比线性计算结果要小。这是由于非线性分析过程中，会考虑到叶片由于挥舞和摆振位移而导致的展向长度减小和变形导致的弯曲刚度增强，其中后者为主要影响因素，因此叶尖的变形会比线性计算结果较小。

表2　 叶尖挥舞特性对比

图5　叶尖挥舞位移对比Fig.5 Comparison of the blade tip wave displacement

对于叶片摆振位移，非线性与线性计算结果差别较小,见图6。这是由于摆振方向叶片刚度较大，变形较小，而非线性刚度矩阵是与形变量相关的，在叶片摆振变形不明显，即形变量很小的情况下，非线性全量刚度矩阵与线性刚度矩阵差别非常小。因此，对于刚度较大，变形较小的摆振方向，几何非线性效应不明显。

图6　叶尖摆振位移对比Fig.6 Comparison of the blade tip shimmy displacement

从表3可以看出，对于风力机的气动性能，非线性与线性计算结果有一定差别。风力机气动性能主要与叶片控制点入流速度和展向位置有关。从图7(a)、图7(b)和表4可以看出，线性与非线性计算所得叶片挥舞速度有明显差别，摆振速度差别较小。并且非线性分析会考虑叶片由于弯曲变形导致的展向位置变化，显然叶片弯曲会导致叶片控制点距离叶根的展向长度减小，即气动力力臂减小，从而导致功率的减小，而风轮平面的减小则会导致风轮推力减小。对于兆瓦级以上的大型风力机，气动性能的稍许变化，就会对风能的捕获产生重要影响。因此，能够更准确的计算风力机气动性能，对于提高风力机效率及相关的控制策略，有着重要作用。

表3　气动性能对比

图7　叶尖速度对比Fig.7 Comparison of the blade tip velocities

然后对线性与非线性分析得到的叶根剪力进行对比，如图8，表5所示。挥舞方向叶根剪力的非线性计算结果明显小于线性计算结果，摆振方向则差别较小。这是由于摆振方向叶根剪力的来源主要是重力载荷，而挥舞方向叶根剪力主要为气动力，其受叶片变形影响较大，因此线性与非线性计算结果差别明显。非线性分析能够更准确的计算叶片变形及其对风力机气动载荷产生的影响，从而更准确的计算各部件载荷。因此，非线性分析对于风力机载荷设计与结构强度设计有重要意义。

表4　叶尖速度振幅对比

图8　叶根剪力对比Fig.8 Comparison of the blade root shear forces

线性非线性差别平衡位置/N71000705360.65%振幅/N28172821-0.14%

4结论

本文计算了Phase VI叶片与NH1500叶片的稳态性能，并与实验结果进行对比，验证了本文气动计算方法的可靠性与准确性。然后计算了NH1500在线性与非线性分析下的动态响应，并进行对比分析，得出了以下结论：

(1) 非线性分析所得叶片挥舞位移明显小于线性分析结果；非线性与线性分析所得叶片摆振位移差别较小，这是由于摆振方向刚度较大，叶片变形较小，几何非线性效应不明显造成的。

(2) 非线性与线性分析所得叶片挥舞速度差别明显，并且非线性分析能更准确地计算控制点位置变化，而风力机气动性能主要与叶片控制点位置及其入流速度有关，因此两种分析方法得到的气动性能差别明显。

(3) 非线性分析能够更准确的计算风力机部件载荷，对于风力机载荷设计与结构强度设计有重要意义，能够更好的保障风力机安全运行。

(4) 非线性分析能够更好的考虑由于叶片变形带来的几何非线性效应，更准确地计算风力机动态响应，更适用于大型风力机模拟计算。

风力机大型化、柔性化发展，随之叶片变形越来越大，线性分析势必会造成资源不合理分配和计算准确度低，采用非线性分析能够更准确地模拟叶片变形，对于更合理的分配资源，提高结构强度与寿命，以及准确模拟风力机气动性能，有着重要作用。

在考虑了几何非线性后，叶片挥舞位移计算结果明显减小，因此可以在保证安全系数的同时，进一步降低叶片的柔性或质量，从而降低叶片生产成本；功率系数的改变则需要考虑新的功率控制方案来改善风力机性能。对于更大型的风力机，或者柔性更强的叶片，非线性与线性响应的差别则会更大。因此，在风力机大型化、柔性化的趋势下，非线性分析对模拟风力机响应是必不可少的。

参 考 文 献

[1] Yuan W B. Geometric nonlinear bending response and buckling of angle-section beams under pure bending [J]. Journal of Constructional Steel Research, 2013, 85:73-80.

[2] Almeida Carlos A, Albino Juan C R, Menezes Ivan F M, et al. Geometric nonlinear analysis of functionally graded beams using a tailored Lagrangian formulation [J].Mechanics Research Communications, 2011, 38(8):553-559.

[3] Sebastian T, Lackner M A. Development of a free vortex wake method code for offshore floating wind turbines [J]. Renewable Energy, 2012, 46: 269-275.

[4] Jeong M S, Yoo S J, Lee I. Wind turbine aerodynamics prediction using free-wake method in axial flow[C]//4th International Symposium on Physics of Fluids. Lijiang, China, 2012.

[5] Zhou W P, Tang S L, Lu H. Computation on aerodynamic performance of horizontal axis wind turbine based on time-marching free vortex method [J]. Proceedings of the CSEE, 2011,31(29):124-130.

[6] Zheng Y Q, Zhao R Z, Liu H. Finite element modal analysis of large-scale composite wind turbine blade [J]. Advanced Materials Research, 2013,694/695/696/697:453-457.

[7] Zhu S F, Rustamov I. Structural design and finite element analysis of composite wind turbine blade [J]. Key Engineering Materials, 2013,525/526: 225-228.

[8] Wu A, Sun W L, Wang P P. Large wind turbine dynamic modeling and simulation [J].Advanced Materials Research, 2010, 34/35:1098-1103.

[9] Riziotis V A, Politis E S, Voutsinas S G, et al. Stability analysis of pitch-regulated, variable-speed wind turbines in closed loop operation using a linear eigenvalue approach [J]. Wind Energy, 2008, 11(5):517-535.

[10] 莫文威，李德源，夏鸿建，等. 水平轴风力机柔性叶片多体动力学建模与动力特性分析[J]. 振动与冲击,2013,32(22):99-105.

MO Wen-wei, LI De-yuan, XIA Hong-jian, et al. Multibody dynamic modeling and dynamic characteristics analysis of flexible blades for a horizontal axis wind turbine [J]. Journal of Vibration and Shock, 2013, 32(22):99-105.

[11] 马剑龙，汪建文，董波，等. 风力机风轮低频振动特性的实验模态研究[J]. 振动与冲击,2013,32(16):164-170.

MA Jian-long, WANG Jian-wen, DONG Bo, et al.Experimental modal analysis on low-frequency vibration characteristics of wind turbine [J]. Journal of Vibration and Shock, 2013, 32(16):164-170.

[12] 巴斯 K J. 工程分析中的有限元法[M]. 北京：机械工业出版社,1991.

[13] 王勖成. 有限单元法[M]. 北京：清华大学出版社,2003.

[14] Martin O L H. Aerodynamics of wind turbines [M]. Oxford:Taylor & Francis Group,2008.

[15] 程危危，曹登庆，初世明. 风力机叶片动力学建模研究[J]. 动力学与控制学报, 2011, 9(4):342-347.

CHENG Wei-wei, CAO Deng-qing, CHU Shi-ming. Dynamics modeling for the blade of a wind turbine [J]. Journal of Dynamics and Control,2011, 9(4):342-347.

[16] 彭细荣，杨庆生，孙卓. 有限单元法及其应用[M].北京:清华大学出版社,2012:369-373.

[17] Hand M M, Simms D A, Fingersh L J, et al.Unsteady aerodynamics experiment phases Ⅵ: Wind tunnel test configurations and available data campaigns[R]. Golden: National Renewable Energy Laboratory, NREL/TP-500-29955, 2001.

基金项目：国家973计划项目(2014CB046200)大型风力机的关键力学问题研究及设计实现；江苏省自然科学基金(BK20140059)；江苏高校优势学科建设工程资助项目

收稿日期：2015-04-28修改稿收到日期：2015-08-06

通信作者王同光 男，博士，研究员，1962年生

中图分类号:O355;TK89

文献标志码：A

DOI:10.13465/j.cnki.jvs.2016.14.030

Geometric nonlinear dynamics response of large-scale wind turbine

CAO Jiu-fa1, WANG Tong-guang2, WANG Xiao2

(1. School of Hydraulic, Energy and Power Engineering, Yangzhou University, Yangzhou 225127, China;2. Jiangsu Key Laboratory of Hi-tech Research for Wind Turbine Design, Nanjing University of Aeronautics & Astronautics, Nanjing 210016, China)

Abstract:With the development of large-scale flexible wind turbine, the deflection of the blades gets larger, which makes the geometric nonlinear effect more obvious. A free vortex wake model was presented to improve the accuracy of the aerodynamics load. The blade was assumed to be a beam model to analyze the structural response with the finite element method. Then a nonlinear kinetic equation was built, in which the geometric nonlinear stiffness matrix was derived by the Galerkin method. The equation was solved with the Newmark immediate integration method and the corrected Newton-Raphson method for time marching. Thus, the coupling simulation between the aerodynamics performance and the geometric nonlinear structural response of the large-scale wind turbine blade was realized. The computation results of the NREL Phase VI rotor agree well with the experimental data, validating the prediction accuracy of the calculation model. The linear and nonlinear dynamics responses of the large-scale wind turbine NH1500 were then calculated. The results of the geometric nonlinear analysis indicate the influence on the structural response and aerodynamics performance. The model is helpful to improve the design standard of the large-scale wind turbine and the accuracy of the load calculation.

Key words:wind turbine; free vortex wake method; geometric nonlinearity; Newton-Raphson method; dynamics response

第一作者 曹九发 男，博士生，1986年生