第一作者柯世堂男，博士，副教授，1982年生

偏航状态下风力机塔架-叶片耦合结构气弹响应分析

柯世堂，王同光

(南京航空航天大学江苏省风力机设计高技术研究重点实验室，南京210016)

摘要：提出一套快速预测偏航状态下风力机全机结构气弹响应的分析方法。以南京航空航天大学自主研发的5MW特大型概念风力机为例，建立风力机塔架-叶片耦合模型获取模态信息；采用谐波叠加法和改进的叶素-动量理论计算气动荷载，并考虑了偏航角对诱导速度的影响；再运用模态叠加法求解风力机耦合动力学方程，通过迭代循环更新叶片速度和气动力，对风力机塔架-叶片耦合结构进行气动载荷和气弹响应计算，并通过参数分析归纳出偏航角和气动弹性对风力机全机动态响应的影响规律。研究结论可为此类特大型风力机塔架-叶片耦合结构的抗风设计提供科学依据。

关键词：风力机塔架-叶片耦合结构；偏航角；气动载荷；叶素-动量理论；气弹响应

基金项目：国家重点基础研究计划(2014CB046200);中国博士后科学基金特别资助(2015T80551)

收稿日期：2014-06-24修改稿收到日期：2014-09-03

中图分类号:TK83；TP391.9文献标志码：A

Aero-elastic vibration analysis based on a tower-blade coupled model of wind turbine in yaw condition

KEShi-tang,WANGTong-guang(Jiangsu Key Laboratory of Hi-Tech Research for Wind Turbine Design, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China)

Abstract:A fast method to calculate aero-elastic responses of wind turbine based on a tower-blade coupled structure model was proposed. By taking the 5 MW wind turbine system designed by Nanjing University of Aeronautics and Astronautics as an example, a finite element model for investigating the wind turbine tower-blade coupled vibration was established to obtain the information of its dynamic characteristics. The harmonic superposition method and the modified blade element momentum theory were applied to calculate the aerodynamic load, considering the influence of yaw conditions. The mode superposition method was used to solve the kinetic equation of wind turbine system, the blade velocity and dynamic load were updated through iterative loop, and then the aero-elastic responses of wind turbine system were calculated. The influence of yaw angle on wind-induced responses was discussed. The research contributes a scientific basis to the wind-resistant structure design for the tower-blade system of large-scale wind turbines.

Key words:wind turbine tower-blade coupled system; yaw angle; aerodynamic load; blade element momentum theory; aero-elastic response

水平轴风力机系统结构(塔架和叶片)属于典型的风敏感结构，叶片的气动载荷和结构风效应是风力机抗风设计关注的重点内容[1]。伴随着全球风能开发规模的迅速增长，风力机结构也朝着特大柔性叶片和超高大塔架方向发展[2]。相比较传统风力机结构，特大型风力机结构的气动弹性效应愈加明显。且在日常运行中，由于风向连续变化和高湍流的影响，叶片转轴不能及时保持与风向平行，会出现某些情况下处于偏航状态，诱导速度将在方位角和径向上发生变化，将会改变叶片的气动载荷分布和引起塔架的气弹失稳或抖振疲劳等问题。

针对偏航状态下的风力机气弹响应问题，廖明夫等[3]对偏航状态下塔架的扭转振动进行理论研究，并给出了有效的减震措施；陈佳慧等[4]对偏航状态下风力机叶片的气弹响应进行研究，并探讨了离心力和气动弹性对叶片动态响应的影响；刘雄等[5]建立了包含偏航角参数的风力机叶片轴向诱导流动因子的计算模型，改进了叶片气动和结构的耦合求解模型；牛蔺楷等[6]基于坐标变换方法对风力机偏航状态下轴承气动载荷进行分析，提出了接触载荷的经验计算公式；查顾兵[7]建立了偏航状态下风力机叶片动态失速预测模型，探讨了偏航导致的动态失速发生的基本原理和特征。已有研究[8-11]大多忽略了风力机塔架和叶片耦合模型的气动弹性影响，随着百米量级塔架的增多，随机风作用下塔架变形不能忽略，其引起的叶片气动弹性问题和偏航工况共同影响风力机的风致动态响应。因此，准确计算风力机塔架-叶片耦合结构在偏航状态下的气动载荷和气弹响应，对提高风力机结构的抗风安全性具有重要意义。

本文基于谐波叠加法和改进的叶素动量理论，并考虑偏航对诱导速度的修正，计算了考虑风剪切、塔影、旋转和塔架-叶片相干效应的风力机气动载荷；再基于有限元方法建立塔架-叶片全机结构模型，运用模态叠加法求解风力机耦合动力学方程，通过迭代循环更新叶片速度和气动力，对不同偏航角下风力机全机结构进行气动载荷和气弹响应计算。

1计算方法描述

风力机风场可以分解为两部分：平稳风和湍流风。前者是宏观上的大气整体运动形成的，方向为水平纵向，数值与观测点高度相关；后者是局部的湍流运动，有纵向、横向和垂直向三个方向，而顺风向脉动风场能量相比横向和垂直向脉动风场来说是叶片和塔架主要承受的载荷，本文主要进行叶片和塔架的顺风向气动载荷和气弹响应分析。风力机气动载荷分析的三维坐标系见图1。

图1　风力机风场模型坐标系示意图 Fig.1 The wind field model coordinate system of the wind turbine

1.1气动载荷计算

采用修正的叶素动量理论计算风力机气动力，引入叶根和叶尖损失，在轴向诱导因子较大时使用Ct的经验模型，并加入动态入流和动态失速模型[12-13]。使用该方法，可以计算不同风速、转速、桨矩角的气动力。流经叶片的相对速度Vrel采用下式计算

(1)

式中：vrot为叶片旋转导致的线速度，W为诱导速度，Vbx、Vby为叶片振动速度，Vox、Voy为沿顺风向和横风向的来流风速，采用谐波合成法[14]计算，考虑叶片和塔架之间的相干性，计算公式如下

j=1,2,3，…,n

(2)

式中：风谱在频率范围内划分成N个相同部分；Δω=ω/N为频率增量；|Hjm(ωl)|为基于改进的Von Karman来流风谱矩阵进行Cholesky分解获得的下三角矩阵的模；θml为介于0和2π之间均匀分布的随机数，可采用Matlab的随机数生成函数，每次生成随机数后应恢复初始状态；ωl=l·Δω是频域的递增变量；ψjm(ωl)为两个不同作用点之间的相位角，是由Hjm(ωl)的虚部和实部的比值确定。其中改进的Von Karman风谱模型[15]和Davenport相干系数模型为

(3)

(4)

风力机的平稳风速由于受到风剪切和塔影效应的影响，风场模拟时必须要对平稳风进行修正[16]。其中风剪切采用指数模型，塔影效应采用适用于叶片在塔架上风向运行的潜流模型。指数模型如下

(5)

式中:V(h)为指高度h处的风速，V(h0)为指参考高度轮毂h0处的风速。当不考虑风剪切的影响时，可以将α的值设为0，取值一般范围为0.1～0.25。h0为轮毂的位置。潜流模型修正公式为：

V(x,z)=AV0

(6)

式中:

(7)

式中:DT为开始考虑塔影影响的高度处的塔架直径；F为塔架直径修正因子；z为计算点到塔架中心的纵向距离；x为风矢量经过时距离塔架中心横向距离。

这样，诱导速度W可由下式表示：

(8)

式中：B为叶片数，L为指升力，φ为入流角，ρ为空气密度，r为叶片截面的展向位置，n为推力方向的单位向量，F为普朗特叶尖损失因子，fg为Glauert修正。本文同时采用了动态入流模型和动态失速模型，修正叶片运转的非定常效应。

考虑风力机已经偏航(见图2)，诱导速度会产生一个偏角，此时指向上游的叶片比指向下游的叶片诱导速度要小，这是由于指向下游的叶片比指向上游的叶片更深入地进入尾流区域，经历了更大的风速，因此产生的气动载荷更大。此时由式(8)计算得到的诱导速度不再适用，本文采用格劳沃特提出的偏航模型对其进行修正。

图2　风力机叶片偏航示意图 Fig.2 Schematic diagram of the rotor in yaw condition

(9)

式中尾流斜交角x定义为尾流中的风速与叶片转轴的夹角(见图1)。W0为由方程求出的诱导速度。θ0为叶片指向尾流最深处时的角。斜交角x通过下式求得：

(10)

式中：n为一个法向量，指向旋转轴方向(见图1)。假设斜交角沿半径方向是常数，可以在接近r/R=0.7的径向位置计算。根据下式计算叶片攻角α

α=φ-(β+θtwist)

(11)

(12)

通过叶片翼型插值方法，可以得到升力系数Cl和阻力系数Cd，从而计算出升力L和阻力D

(13)

这样得到叶片的法向载荷Fn和切向载荷Ft

(14)

1.2风力机塔架-叶片耦合模型动力特性

基于南京航空航天大学自主研发的5 MW三叶片特大型概念风力机系统，该型号风力机是课题组依据以前开发的NH1500(1.5 MW)和NH3000(3 MW)风力机的结构概念延伸，主要是叶片和塔架的结构参数，没有涉及到发电机、齿轮箱和控制系统等内部结构。由于本文主要研究对象是风力机叶片-塔架耦合结构的气动问题，因此把机舱内的发电机和控制系统等效为一个刚性质量块，主要研究叶片和塔架结构的气动载荷和气弹响应。主要参数有：塔高124 m，底径4.8 m，顶径2.6 m，塔体通长为变厚度结构，底壁厚150 mm，顶壁厚60 mm，通长厚度由底部至顶部呈线性减小趋势。机舱长12 m，宽4.6 m，高4.2 m，总质量140.2×103kg。各叶片之间成120°夹角，沿周向平均分布，叶片直径为120 m，宽度2.4 m，厚度0.38 m，长度60 m，偏航角为0°，额定转速为17 r/min。

传统风力机的气弹响应研究大多假设塔架是完全刚性，仅考虑叶片的动态变形对气动载荷的影响。事实上，由于叶片载荷的存在使得超高柔塔架的变形越发明显，而塔架的变形又将影响叶片的动态响应，从而改变叶片的速度和气动力。因此，基于ANSYS软件平台，建立了风力机塔轮系统的“叶片-塔体”一体化有限元模型。其中风轮和塔体采用SHELL91单元，机舱及其内部结构可作为整体采用梁单元BEAM189模拟，圆形筏基基础采用SOLID65单元模拟，基础直径为10 m，高度为1.8 m，基础底部直接固结，不考虑土体和结构的相互作用。通过多点约束单元耦合命令将各部分连接在一起。模态分析时把叶片旋转产生的离心力作为预应力预先施加在叶片上，计算的频率和模态信息均考虑叶片转动带来的离心力效应。

有限元模型和频率分布曲线见图3。系统的第1阶模态频率很低，仅为0.27 Hz，第10阶模态频率为1.85 Hz，模态之间的间隔非常小。最后基于该风力机系统模型提取频率和振型信息，为后续的风致动力计算提供输入参数。从图4系统典型模态振型图可知，第5阶为三片风轮一起做前后舞动运动；第10阶为叶片和机舱带动塔架做前后舞动；第50阶为叶片复杂的前后舞动和左右摆振，并伴随塔架的弯曲变形。分析表明低阶模态下风力机叶片-塔架耦合模型主要以叶片变形为主，随着阵型频率的增加，叶片挥舞摆动则愈加显著，并且伴随着塔架的共同变形。

图3　风力机塔架-叶片耦合模型和频率分布 Fig.3 The finite element model and natural frequencies of wind turbine system

图4　风力机塔架-叶片耦合模型典型模态振型示意图 Fig.4 Typical modes of vibration of wind turbine systems

1.3风力机系统结构动态响应求解

根据达朗贝尔原理，风力机系统结构风致动力响应控制方程为

个人哲学是教师在教学情境中思考自身的方式，包含教师个人心智中的信念和价值观，作为一种无意识的经验假设支配着教师的行为。从个人角度来看，哲学是一种让人认识世界、了解世界的工具。由于每位教师的性格、成长环境、经验不同，其个人哲学也迥然不同，但相似的是，教师的个人哲学往往贯穿于生活和工作中，统领着他们的信念与价值观，指挥着他们的行为方式，使教师成为独特的个体。它通常包涵着个体的信念、价值观和行动原则，透过信念和价值观的外在表现，深入到经验中，成为教师实践性知识的个性化表征。

通过将节点位移向量从物理坐标转换到广义模态位移后，实现方程解耦，再引入模态阻尼，则各广义模态对应的运动方程为

i=1,2,3,…,n

(16)

式中:qi(t)为第i阶广义模态位移响应向量；mi、ci、ki、pi(t)分别为第i阶模态的模态质量、模态阻尼、模态刚度和模态力。这样单自由度的运动方程可根据Duhamel积分原理，初始条件为零，数值解为

i=1,2,3,…,n

(17)

式中:ωni=(ki/mi)1/2为结构的固有模态频率，ζi=bi/(2miωni)为模态阻尼比，ωdi=ωni(1-ζ2)1/2为结构阻尼振动频率，△τ为积分时间步长。

再将各模态下的广义位移转换为物理位移并进行叠加，可得到各节点动态响应位移

(18)

式中:Φ为系统振型矩阵。

风力机塔架-叶片耦合模型风致响应时域计算基于有限元软件ANSYS平台，将数值模拟得到的风荷载作为外部激励作用于风力机塔架-叶片耦合模型上，采用Newmark-β逐步积分法和Newton-Raphson迭代理论，其中各模态阻尼均为0.02，积分时间步长取为0.05 s，响应输出时间步长取为0.025 s，截取模态取系统的前50阶。

1.4考虑偏航角风力机塔架-叶片耦合结构气弹响应计算方法

综上所述，图5整理出偏航状态下考虑风剪切、塔影、旋转和相干性的风力机塔架-叶片耦合结构气弹响应计算流程。

图5　偏航状态下风力机结构气弹响应计算流程图 Fig.5 Flow chart of calculating aero-elastic responses for wind turbines under yaw condition

(1)基于改进的Von Karman风谱模型和Davenport相干函数模拟考虑塔架-叶片相干性的来流风速，将其作为后续改进的BEM方法计算的输入样本，并对其进行风剪切和塔影效应修正；

(2)采用改进的叶素-动量理论计算考虑旋转和相干效应的风力机诱导速度，并对其进行偏航状态的修正，进而获得风力机系统的气动力；

(3)建立风力机塔架-叶片耦合模型，求解对应气动力作用下的结构动态响应，并将结构响应返回引入气动力，更新叶片诱导速度和气动力，进行迭代循环计算获得风力机塔架-叶片耦合结构的气弹响应时程，直到满足预先设定的截止时间或样本数目为止。

2算例分析

为验证本文方法模拟风力机气动载荷的正确性，通过自编程序进行风场模拟和气动力计算，获得了不同风速偏航状态下叶片的功率系数和平面内/外气动载荷(是气动升力和阻力的组合载荷)的分布特性，并与GH Bladed软件的功率系数计算结果进行对比(见图6和图7)。

对比结果显示本文方法的计算结果比GH Bladed软件的计算结果要稍低，在11 m/s风速对应的最大功率系数值误差约5%，基本上两者的结果吻合较好，验证了本文方法模拟气动力的正确性。图7为偏航状态下计算获得的叶片平面外气动载荷明显小于无偏航的计算结果，而偏航状态下叶片平面内气动载荷与不考虑偏航的计算曲线比较吻合，说明其对平面内的气动载荷影响较小。

图6　偏航状态下风力机功率系数计算结果对比图 Fig.6 Power coefficient contrast of wind turbine under yaw condition

图7　考虑/不考虑偏航工况叶片气动载荷分布曲线 Fig.7 Distribution curves of aerodynamic loads of blades with or without yaw condition

为了深入分析偏航角对风力机塔架-叶片耦合结构动态响应的影响机理，图8和图9分别给出了不同偏航角下叶尖舞动和摆振两个方向气弹响应时程曲线。可以看出：

(1)不同偏航角下的叶尖舞动和摆振方向动态气弹响应都以长周期成分为主，随着偏航角的增大长周期成分越加明显，且位移时程曲线中存在间歇性的高幅响应数值；

(2)偏航角对叶片前后舞动方向位移影响显著，随着偏航角度的增大，叶片挥舞方向动态位移的振动幅度变大，但均值明显减小。这是由于偏航角的增大使得叶片平面法向的气动载荷逐渐变小，所以导致叶片的动态位移均值逐渐减小；然后由于诱导速度指向下游比指向上游更快的进入尾流区，以致叶片指向下游的诱导速度要更大，从而引起气动载荷的变化，并且随着偏航角的增大气动载荷变化也越大，因此引起的叶片挥舞方向振动幅度变大；

(3)偏航角对叶片左右摆振方向的位移影响微弱，随着偏航角度的变化，叶片摆振方向的动态响应曲线并无明显变化特征。主要是因为叶片摆振方向主要包含气动力和重力载荷，偏航状态下叶片所受载荷周期与旋转周期接近，重力载荷也随叶片旋转呈周期性变化，此时偏航角的改变对于叶片摆振方向动态响应的影响较弱。

图8　不同偏航角下叶片挥舞方向气弹响应时程 Fig.8 The aero-dynamic responses in back-forth waving direction under different yaw angles

图9　不同偏航角下叶片摆振方向气弹响应时程 Fig.9 The aero-dynamic responses in right-left swing direction under different yaw angles

表1给出了不同偏航角工况下风力机塔架-叶片耦合结构的典型目标气弹响应的均值、根方差和极值对比。分析发现，偏航角的增大会减小塔顶顺风向位移均值，但明显增大其根方差，计算得到的极值随着偏航角的增大显著变大，当偏航角为30°时，增幅达到20%。偏航角对塔底弯矩或叶片剪力的影响明显没有塔顶位移显著，特别是对基底弯矩极值的影响，偏航角为30°时引起的极值增幅仅为3.6%，但变化规律与塔顶位移相同。这就说明偏航角下风力机塔架-叶片耦合结构气弹响应分析时叶尖和塔顶的气弹位移变化最为明显、其次是叶根剪力、而基底弯矩影响最弱。

表1　不同偏航角下风力机系统典型目标气弹响应参数对比

为考虑气动弹性对风力机系统结构响应的影响，图10分别计算了引入迭代循环前后塔架顶部位移随时间变化的曲线。对比可知，考虑气动弹性会减小风力机系统的动态响应根方差，特别是在高风速作用下的计算误差最为明显。因此，气动弹性对准确计算风力机气动载荷和动态响应影响较大，在抗风计算中应加以考虑。

图10　不同风速下风力机塔架顶部顺风向位移响应均方差 Fig.10 RMS of along-wind displacement responses of tower top under different wind speed

3结论

本文分析了风力机塔架-叶片耦合结构的动力特性，计算了偏航状态下风力机系统气动载荷，并采用时域模态叠加法求解了风力机全机结构的气弹响应。主要得到以下结论：

(1)本文提出的偏航状态下风力机整机结构气弹响应计算方法可以准确地计算系统气动载荷和气弹响应，与商用软件相比，本文方法在系统气动载荷模拟和结构模态分析等方面有一定的改进；

(2)偏航角的存在会明显减小叶片平面外气动载荷，但对平面内的气动载荷分布影响微弱；

(3)偏航角对风力机塔架-叶片耦合结构的气弹响应以叶尖和塔顶位移最为显著、其次是叶根剪力、最弱的是基底弯矩。影响规律均是随着偏航角的增大，气弹响应均值逐渐减小，但根方差明显变大，导致其极值响应均不同程度的变大；

(4)气动弹性会减小风力机系统的动态响应根方差，尤其是高风速作用下的计算误差不能忽略。

参考文献

[1]Ronold K O, Larsen G C. Optimization of a design code for wind turbine rotor blades in fatigue[J]. Engineering Structure, 2001, 23: 993-1002.

[2]任勇生, 张明辉. 水平轴风力机叶片的弯扭耦合气弹稳定性研究[J]. 振动与冲击, 2010, 29(7): 196-200.

REN Yong-sheng, ZHANG Ming-hui. Aeroelastic stability study on coupled flutter for horizontal axis wind turbine blades[J]. Journal of Vibration and Shock, 2010, 29(7): 196-200.

[3]廖明夫, 黄巍, 董礼, 等. 风力机偏航引起的失稳振动[J]. 太阳能学报, 2009, 30(4): 488-492.

LIAO Ming-fu, Huang Wei, Dong Li, et al. Fatigue characteristics analysis of wind turbine tower under wind-wave combined effect[J]. Acta Energiae Solaris Sinica,2009,30(4): 488-492.

[4]陈佳慧, 王同光. 偏航状态下的风力机叶片气弹响应计算[J]. 南京航空航天大学学报,2011, 43(5): 629-635.

CHEN Jia-hui, WANG Tong-guang.Aeroelastic responses calculation of wind turbine blade in yaw condition[J]. Journal of Nanjing University of Aeronautics and Astronuatics,2011, 43(5): 629-635.

[5]刘雄, 陈严, 叶枝全. 风力机气动与结构 CAD 软件[J]. 太阳能学报, 2001, 22(3): 346-350.

LIU Xiong, CHEN Yan, YE Zhi-quan. Wind turbine aerodynamic performance and structure CAD sortware[J]. Acta Energiae Solaris Sinica, 2001, 22(3): 346-350.

[6]牛蔺楷, 杨洁明, 高俊云. 基于坐标变换方法的风力机偏航轴承载荷分布的分析与计算[J]. 工程力学, 2012, 29(10): 282-287.

NIU Lin-kai,YANG Jie-ming,GAO Jun-yun. Determination of load distribution in yaw bearing of wind turbine using coordinate transformation method[J]. Engineering Mechanics, 2012, 29(10): 282-287.

[7]查顾兵, 竺晓程, 沈昕, 等. 水平轴风力机在偏航情况下动态失速模型分析[J]. 太阳能学报,2009, 30(9): 1297-1300.

ZHA Gu-bin, ZHU Xiao-cheng, SHEN Xi. Dynamic stallmodeling of horizontal axis wind turbine in yaw considion[J]. Acta Energiae Solaris Sinica,2009, 30(9): 1297-1300.

[8]Bazilevs Y, Hsu M C, Kiendl J, et al. 3D simulation of wind turbine rotors at full scale. Part II: Fluid-structure interaction modeling with composite blades[J]. International Journal for Numerical Methods in Fluids. 2011, 65(1):236-253.

[9]Tempel J V D. Design of support structures for offshore windturbines[D]. Netherlands: Delft University of Technology, 2006.

[10]Wang T G, Wang L, Zhong W, et al. Large-scale wind turbine blade design and aerodynamic analysis[J]. Chinese Science Bulletin,2012, 57(5):466-472.

[11]陈小波, 陈建云. 海上风力发电塔脉动风速时程数值模拟[J]. 中国电机工程学报,2008, 28(32): 111-116.

CHEN Xiao-bo, CHEN Jian-yun. Numericalsimulation of fluctuating wind velocity time series of offshore wind turbine[J]. Proceedings of the CSEE,2008, 28(32): 111-116.

[12]伍艳, 谢华, 王同光. 风力机叶片的非定常气动特性计算方法的改进[J]. 工程力学,2008, 25(10): 54-60.

WU Yan, XIE Hua, WANG Tong-guang. Modification of calculating unsteady aerodynamic characteristics of wind turbine blades[J]. Engineering Mechanics,2008, 25(10): 54-60.

[13]柯世堂, 王同光, 赵林,等. 风力机风振背景、共振响应特性及耦合项分析[J]. 中国电机工程学报, 2013, 33(26): 101-108.

KE Shi-tang, WANG Tong-guang, ZHAO Lin, et al. Background, resonant components and coupled effect of wind-induced responses on wind turbine systems[J]. Proceeding of the CSEE, 2013,33(26): 101-108.

[14]Kareem A．Numerical simulation of wind effects：a probabilistic perspective[J]．Journal of Wind Engineering and Industrial Aerodynamics，2008(96)：1472-1497.

[15]廖明夫. 风力发电技术[M]. 西安:西北工业大学出版社, 2008.

[16]Burton T, Sharpe D, Jenkins N, et al. Wind energy handbook [M]. Chichester, John Wiley&Sons,2001.