风力机塔架叶片耦合模型风致响应时域分析
2014-09-27柯世堂曹九发王珑王同光
柯世堂+曹九发+王珑+王同光
文章编号：16742974(2014)04008707
收稿日期：20130728
基金项目：国家重点基础研究计划资助项目(2014CB046200)；国家自然科学基金资助项目(51208254)；江苏省自然科学基金资助项目(BK2012390)；江苏高校优势学科建设工程资助项目
作者简介：柯世堂（1982-），男，安徽池州人，南京航空航天大学讲师，博士
通讯联系人，E-mail:keshitang@163.com
摘要：基于风力机塔架叶片耦合模型，采用改进的叶素动量理论模拟了考虑平稳风修正、叶片旋转效应和空间相干性的风力机气动载荷，并基于有限元方法对该耦合模型进行了动力特性分析和风致响应时域计算．基于目标响应时程探讨了风力机塔架叶片耦合系统在随机风荷载作用下的动力响应特性，并与不考虑叶片影响的风力机塔架风致响应进行对比分析.研究表明,在进行风力机的抗风设计时，应该考虑塔架叶片的耦合作用.
关键词：风力机塔架叶片耦合模型；风场模拟；风致响应；时域分析
中图分类号：TK83；TP391.9 文献标识码：A
TimedomainAnalysisoftheWindinducedResponsesofthe
CoupledModelofWindTurbineTowerbladeCoupledSystem
KEShitang,CAOJiufa，WANGLong,WANGTongguang
(JiangsuKeyLaboratoryofHiTechResearchforWindTurbineDesign,NanjingUniv
ofAeronauticsandAstronautics,Nanjing，Jiangsu210016,China)
Abstract:Basedonthewindturbinetowerbladecoupledmodel,themodifiedbladeelementmomentumwasusedtosimulatetheaerodynamicloadsofwindturbinemodel,whichconsiderstherotationaleffectcausedbybladesandbladetowerinteraction.Thewindturbinetowerbladecoupledmodelconsideringthecentrifugalforcefromrotationalbladeswasestablishedinfiniteelementmethod,andthetimedomainanalysisofwindinducedresponsesofwindturbinetowerbladecoupledmodelwascarriedout.Then,thedynamicresponsecharacteristicsofthewindturbinesystemunderstochasticwindloadswerediscussedonthebasisofthecalculatedresultsandthroughthecontrastwiththewindinducedresponsesofatowerwithouttheblades.Itisproposedthatthewindturbinetowerbladecoupledmodelshouldbeadoptedinwindresistingdesign.
Keywords:windturbinetowerrotorcoupledsystem;fluctuatingwindfieldsimulation;windinducedresponses;timedomainanalysis
风荷载是风力机系统设计最重要荷载之一[1].随着兆瓦级大功率风力机的普遍应用，叶片的旋翼直径达到百米量级，从而满足捕捉风能的要求，但同时也遭受了强大风推力，并通过机舱将风荷载传递给塔架，使得叶片和塔架之间的耦合效应愈加明显.以往针对风力机系统的抗风研究，关注风力机空气动力学性能的研究者大多忽略塔架的影响，对叶片的流场特性[2-3]和气弹响应[4-5]进行了深入详细的研究；而关注风力机结构受力性能的研究者大多以塔架为研究对象，把叶片遭受的风载等效为拟静力施加到塔顶进行风振计算，或是建模时考虑塔架叶片的耦合模型，而风场模拟时没有完全考虑叶片平稳风修正、旋转效应和气弹效应等影响[6-8].
目前，对风力机全机结构的抗风设计，基本上采用等效风荷载的拟静力分析方法，少量文献进行了频域内的风振分析，考虑到塔架叶片耦合模型脉动风场模拟的复杂性，在时域内对风力机塔架叶片耦合模型进行风致响应计算则更为复杂.然而时域分析可以克服频域分析中基于线性化假设的不足，同时时域计算结果的直观性有利于工程师了解风力机系统的受力性能.因此，基于风力机塔架叶片耦合模型的脉动风荷载模拟结果，基于塔架叶片耦合有限元模型进行风振响应时域分析，对风力机系统抗风的精细化设计有重要的指导作用和实用价值.
鉴于此，本文基于时域方法分析风力机塔架叶片耦合模型在顺风向风荷载动力作用下的响应特性，建立了风力机塔架叶片耦合的计算模型，运用谐波合成法和叶素动量理论进行数值模拟获得了考虑平稳风修正、叶片旋转效应和相干性的塔架叶片耦合模型气动载荷，在此基础上利用有限元方法进行了风力机系统风致响应时域计算，基于计算结果，分析了风力机塔架叶片耦合模型的风致响应特性.
1风力机塔架叶片耦合模型风场模拟
由于风力机塔架叶片耦合系统结构本身的复杂性及随机脉动风作用的不确定性，风场模拟时需作如下简化：①只考虑顺风向风荷载数值模拟和作用；②脉动风被假定为零均值的平稳高斯随机过程；③不考虑叶片和气流之间的气动弹性作用.
1.1平稳风修正
风力机的平稳风速由于受到风切变、塔影效应和上游风机尾流的影响，风场模拟时必须要对平稳风模型进行修正[9].其中风剪切主要采用指数模型，塔影效应主要采用适用于叶片在塔架上风向运行的潜流模型，上游尾流影响主要采用包含由于尾
流引起的附加湍流的涡流粘度模型.
风剪切的影响是指平稳风速随着高度而变化，其指数修正模型表达公式如下：
V(h)=V(h0)(hh0)α.(1)
式中:V(h)为指高度h处的风速；V(h0)为轮毂的参考高度h0处的参考风速；h0为轮毂的位置；当不考虑风剪切的影响时，可以将α的值设为0，取值一般为0.1~0.25.
塔影是由于风力机塔架的存在影响了风场平稳风速，其影响修正主要有3种模型：叶片在塔架上风向运行的潜流模型、叶片在塔架下风向运行的经验模型和两者的组合模型.现有风力机基本都是上风向风机设计，用以避免叶片周期地通过塔架尾迹产生的附加噪声和激振力，因此采用潜流模型修正：
V(x,z)=1+(FDT2)2(x2－z2)(x2+z2)2V0.(2)
式中:DT为开始考虑塔影影响的高度处的塔架直径；F为塔架直径修正因子；z为计算点到塔架中心的纵向距离；x为风矢量经过时距离塔架中心横向距离.
由于风力机很少会存在单个运营状态，基本都是以群体形式建立.因此，在模拟风力机模型的风场时，当由于单个风力机系统转子部分或者全部处于上游风力机的尾流中时，就必须考虑上游风力机尾流对风力机风场的影响.本文采用涡流粘度模型来考虑尾流引起的附加湍流影响，修正公式如下：
V=V0(1－Δe－r22w2).(3)
式中：r为局部速度到尾流中心的距离；w为尾流宽度；V0为未受扰动平均风速；Δ为速度逆差.
1.2脉动风模拟方法
已有研究表明，风力机叶片的风场模拟可采用叶素动量理论[10]，相较涡尾迹方法和CFD方法，叶素动量理论可以实现叶片旋转风场的模拟，且计算量较小、耗时较短，满足风力机日常快速计算的要求.考虑到叶素动量理论是针对叶片在某一特定来流风速和工况下对应参考风速的求解，本文通过谐波合成法[11]模拟获得考虑叶片和塔架相干性的来流风速时程，再对每个来流风速样本采用改进的叶素动量理论进行计算得到该来流风速下叶片的参考风速，从而考虑叶片的旋转效应和相干性，如此迭代循环获得叶片的参考风速时程.风力机塔架叶片耦合模型脉动风模拟流程示意图见图1.
图1风力机塔架叶片耦合模型
风场模拟流程示意图
Fig.1Theflowsketchfigureofwindfieldsimulation
forwindturbinetowerbladecoupledmodel
风谱模型采用改进的Vonkarman模型，因为相比Davenport模型和Kaimal模型，改进的Vonkarman风谱模型能够更加准确地产生随机风场样本，该模型校正了基本模型在高度150m以下的缺陷，更加符合风力机自身风场特性.图2为3种风谱模型曲线，在低频区3种风谱数值接近，在峰值附近区域，Kaimal和Vonkarman频谱分别要低于和高于改进Vonkarman频谱；而在高频区域，Kaimal和Vonkarman频谱几乎重合，但是改进Vonkarman频谱会略高于另外两个频谱幅值.
频率/Hz
图23种标准化风谱曲线示意图
Fig.2Thedifferentlongitudinalnormalpowerspectrum
根据风力机风场的自身特性，选择改进的Vonkarman风谱模型：
fSuu(f)σ2u=β12.987u/a［1+(2πu/a)2］5/6+β2×
1.294u/a［1+（πu/a)2］5/6F1.(4)
式中：Suu为随机风速变化自频谱；f为频率；σu为随机风的标准差；u为无量纲频率参数，u=fLux/U10，U10为10m高基本风速，Lux为随机风纵向分量长度尺度.
采用Davenport相关系数经验公式来考虑叶片和塔架之间的相关性：
Cij=exp×
－ωCxxi－xj2+Cyyi－yj2+Czzi－zj22πvH.(5)
式中：Cx,Cy和Cz分别为叶片和塔架上任意2点横向、顺风向和垂直向的衰减系数；ω为脉动风频率；v(H)为H处平均风速.仅考虑顺风向相关性时，其余两个方向衰减系数均取0.
定义风机上n个风速模拟节点，假定均为零均值的平稳高斯过程，其风谱密度函数矩阵为：
S(ω)=s11(ω)s12(ω)...s1n(ω)s21(ω)s22(ω)...s2n(ω)
sn1(ω)sn2(ω)...snn(ω).(6)
式中:sii(ω)为节点脉动风自功率谱，采用式(2)中的风谱模型计算；sij(ω)为互功率谱，其表达式需要用到塔架和叶片以及塔架上任意两点之间的相干性，计算公式为:
sij(ω)=sii(ω)sjj(ω)Cij.(7)
其中叶片需考虑旋转平面内的各点间的相干性，以及叶片和塔架之间的相干性影响.再将S(ω)进行Cholesky分解，此时风力发电塔架上的任何一个节点脉动风速时程可以由其功率谱决定，根据Shinozuka理论，模拟的风速时程可以表达为:
vj(t)=∑jm=1∑Nl=1Hjm(ωl)2Δωcos［ωlt+
ψjm(ωl)+θml］,j=1，2，3，…，n．(8)
式中：风谱在频率范围内划分成N个相同部分，Δω=ω/N为频率增量；|Hjm(ωl)|为基于Davenport来流风谱矩阵进行Cholesky分解获得的下三角矩阵的模；θml为介于0和2π之间均匀分布的随机数，可采用Matlab的随机数生成函数，建议每次生成随机数后应恢复初始状态；ωl=lΔω为频域的递增变量；ψjm(ωl)为两个不同作用点之间的相位角，它是由Hjm(ωl)的虚部和实部的比值确定.
采用修正的BEM理论[12]，引入叶根损失和叶尖损失，在轴向诱导因子较大时使用Ct的经验模型，并加入动态入流和动态失速模型.使用该方法，可以计算风力机在不同风速、转速、桨距角及偏航角情况下的动态载荷.进而获取作用在叶片上的脉动风速时程.具体方法如下：
根据BEM理论，叶片上的相对风速Vrel采用下式计算：
vrel,xvrel,y=voxvoy+0vrot+WxWy－vbxvby．(9)
式中：vox，voy分别为沿顺风向和横风向的来流脉动风速，采用式(9)谐波叠加法计算；vrot为叶片旋转导致的线速度；W为诱导速度；vbx和vby为叶片振动速度.
诱导速度W可由下式表示：
Wx=－BLcosφ4ρπrFv0+fgn(n•W),
Wy=－BLsinφ4ρπrFv0+fgn(n•W)．(10)
式中：B为叶片数；L为指升力；φ为入流角；ρ为空气密度；r为叶片截面的展向位置；n为推力方向的单位向量；F为普朗特叶尖损失因子；fg为Glauert修正.同时，本文还采用了动态入流模型和动态失速模型，修正叶片运转的非定常效应.
根据下式计算叶片攻角α：
α=φ－(β+θtwist).(11)
式中:β为桨矩角；θtwist为叶片剖面几何扭角，计算公式为:
tanφ=vrel,xvrel,y.(12)
1.3数值模拟结果
本文以某5MW三桨叶变桨距特大型风力机系统为例进行风场数值模拟和风致响应时域分析.塔高124m，底径4.8m，顶径2.6m，塔体通长为变厚度结构，底壁厚150mm，顶壁厚60mm，通长厚度由底部至顶部呈线性减小趋势.机舱长12m，宽4.6m，高4.2m，总质量140.2×103kg.各桨叶间成120°夹角，沿周向平均分布，风轮直径为120m，宽度2.4m，厚度0.38m，长度60m，偏航角为0°，额定转速为17r/min.
基于上述方法，采用Matlab语言编制相应的数值模拟程序MBEMHSM.在计算过程中脉动风上限频率取为2π，脉动风频率分割点数取为2048，频率增量Δω=0.00307Hz.图3和4给出了10m高平均风速为24m/s模拟得到的叶片和塔架中部的风速时程曲线和对应功率谱曲线，为对比显示风力机塔架叶片耦合模型风场的特殊性，风谱坐标系均为对数坐标.
时间/s(a)脉动风速时程曲线
f/Hz(b)脉动风速计算谱校核曲线
图3叶片参考风速时程模拟结果
Fig.3Simulatingresultoffluctuatingwindvelocityofblades
时间/s(a)脉动风速时程曲线
f/Hz(b)脉动风速计算谱校核曲线
图4塔架中部风速时程模拟结果
Fig.4Simulatingresultoffluctuatingwindvelocityoftowermiddlepart
对比发现，考虑塔架叶片耦合模型叶片脉动风速功率谱曲线在高频处存在明显的能量浮动，应该是由于考虑改进的Vonkarman风谱模型、三维旋转和干扰效应而产生的现象，在风力机系统风场模拟中应引起重视；塔架中部的风速功率谱和Davenport风谱吻合较好，这是由于随着高度的降低，塔架受到叶片的干扰效应减小，在高频出现的能量浮动和数值逐渐减小.根据本文模拟过程和对比分析可以认为，采用本文的MBEMHSM方法可以很好地模拟风力机塔架叶片耦合模型风场.
2风力机塔架叶片耦合模型动力特性
基于ANSYS软件平台，建立了风力机塔架叶片耦合模型.其中叶片和塔体采用SHELL91单元，机舱及其内部结构可作为整体采用梁单元BEAM189模拟.通过多点约束单元耦合命令将各部分连接在一起，形成整体的风力机塔架叶片耦合模型.依据效率和精度均衡的原则，模型一共划分了3812个单元.动力特性分析时把叶片旋转产生的离心力作为预应力预先均匀施加在叶片上，后续计算的频率和模态信息均考虑叶片转动带来的离心力作用.图5给出了考虑/不考虑叶片离心力作用2种工况下系统前200阶自振频率的分布情况.
模态阶数
图5风力机塔架叶片耦合系统固有频率分布图
Fig.5Scattergramofnaturalfrequenceforwind
turbinetowerbladecoupledsystem
从图5可以看出：1)风力机塔架叶片耦合模型的基频很小，仅为0.24Hz，第50阶自振频率为9.2Hz，各模态之间的间隔很小，其自振频率数值远远小于相同高度和直径的高耸结构，说明叶片和塔架的耦合效应明显降低了风力机系统的基频，使得系统结构可能会对风荷载更为敏感.2)当考虑叶片转动引起的离心力作用时，耦合系统的基频(0.27Hz)要略大于不考虑离心力作用下的系统基频，并且随着模态数目的增加，离心力效应带来的频率影响越来越大.这是由于风力机塔架叶片耦合系统的高阶频率主要以叶片的左右摆动和前后振动为主，此时叶片旋转引起的离心力效应更为明显，相应增强了风力机系统的刚度.本文的后续计算均采用考虑叶片离心力作用的更加真实的模态参数.
从图6给出的风力机塔架叶片系统模态振型图中可以看出，塔架的弯曲变形与叶片的挥舞/摆振相互耦合.通过多阶模型分析表明,系统振型主要以叶片的前后挥舞和左右摆动为主，高阶模态下塔架本身也会出现结构变形和失稳形态，并且和叶片变形耦合作用.动力特性分析表明,风力机系统动力响应计算必须要考虑叶片和塔架耦合作用.
图6风力发电塔轮系统的典型模态振型示意图
Fig.6Typicalmodesofvibrationonwindrotorssystems
3风致响应时域分析方法
根据达朗贝尔原理，风力机系统结构风致动力响应控制方程为:
My¨(t)+C(t)+
Ky(t)=p(t)．(13)
式中:p(t)为风荷载向量；M,C和K分别为结构质量、阻尼和刚度矩阵；y¨(t),(t)和y(t)分别为节点加速度向量、速度向量和位移向量.
通过将节点位移向量从物理坐标转换到广义模态位移后，实现方程解耦，再引入模态阻尼，则各广义模态对应的运动方程为：
miq¨i(t)+cii(t)+kiqi(t)=pi(t),i=1,2,3,…,n．(14)
式中：qi(t)为第i阶广义模态位移响应向量；mi，ci，ki和pi(t)分别为第i阶模态的模态质量、模态阻尼、模态刚度和模态力.这样单自由度的运动方程可根据Duhamel积分原理，初始条件为零，数值解为：
qi(t)=1miωdi∑tτ=0pi(τ)exp－ξiωni(t－τ)×sinωdi(t－τ)Δτ,i=1,2,3,…,n．(15)
式中：ωni=(ki/mi)1/2为结构的固有模态频率；ζi=bi/(2miωni)为模态阻尼比；ωdi=ωni(1-ζ2)1/2为结构阻尼振动频率；Δτ为积分时间步长.本文选用采用Newmarkβ逐步积分法和NewtonRaphson迭代理论，其逐步积分方法的基本假设为：
(t+Δt)=(t)+(1－δ)U¨(t)+δU¨(t+Δt)Δt,(16)
U(t+Δt)=U(t)+(t)Δt+(12－α)U¨(t)+α(t+Δt)Δt2．(17)
再将各模态下的广义位移转换为物理位移并进行叠加，可得到各节点风致响应位移为：
y(t)=Φq(t)．(18)
式中：Φ为系统振型矩阵.
风力机塔架叶片耦合模型风致响应时域计算基于有限元软件ANSYS平台，将第2节数值模拟得到的风荷载时程作为外部激励作用于风力机塔架叶片耦合模型上，采用Newmarkβ逐步积分法和NewtonRaphson迭代理论，其中各模态阻尼均为0.02，积分时间步长取为0.05s，响应输出时间步长取为0.025s，截取模态取系统的前50阶.
4风力机塔架叶片耦合模型风致响应
当平均风速过大时，风力发电机组会自动停止运行，这个临界风速称之为切出风速，对于兆瓦级风力机通常设定为25m/s左右.首先,进行切出风速下风力机塔架叶片耦合模型风致响应时域计算.限于篇幅，本文只给出了风力机耦合模型叶片尖部、机舱和塔架顶部位移响应时程曲线，如图7所示.
对比发现，基于风力机塔架叶片耦合模型时域计算得到的叶片尖部位移响应以低频共振效应为主，其脉动覆盖范围和数值最大，这是由于风力机遭受的风荷载卓越频率与叶片的摆振/舞动频率较为接近，从而引发了叶片的共振.机舱由于刚度较大，风荷载作用主要导致静力或准静力响应，因此位移响应时程曲线浮动范围较小，均值约在0.295m.塔架顶部的位移响应时程明显小于叶片响应数值，且高频振型贡献更加明显，对应的振型大多为叶片振动和塔架变形的耦合作用.分析表明，考虑了叶片的风场影响和结构模态耦合作用后，塔架的风致响应特性明显呈现多振型响应和多荷载形态.
时间/s(a)叶片尖部位移响应时程曲线
时间/s(b)机舱位移响应时程曲线
时间/s(c)塔架顶部位移响应时程曲线
图7风力机塔架叶片耦合模型典型部位位移响应时程曲线
Fig.7Thedisplacementtimehistoryontypicalnodes
ofwindturbinetowerbladecoupledmodel
为研究叶片对风力机塔架风致响应特性的影响，再对不同风速下考虑和不考虑叶片影响塔架风致响应进行分析.塔顶顺风向位移和塔底弯矩脉动响应均方根随风速的变化曲线如图8和图9所示，由图可见，在低风速(≤10m/s)下是否考虑叶片影响对风力机塔架风致响应影响不大，而随着风速的逐渐增大，考虑叶片影响计算得到的塔架风致响应要远大于不考虑叶片作用时的计算结果，其中考虑叶片影响时塔架顶部位移较未考虑的要大将近2倍，对于塔架底部弯矩来说要大1.5倍.
风速/（m•s-1）
图8不同风速下风力机塔架顶部顺风向位移响应均方差
Fig.8RMSofalongwinddisplacementresponses
oftowertopunderdifferentwindspeed
风速/（m•s-1）
图9不同风速下风力机塔底弯矩响应均方差
Fig.9RMSofbendingmomentresponsesoftower
bottomunderdifferentwindspeed
5结束语
基于风力机塔架叶片耦合模型，实现了给定风环境的系统风场模拟和考虑离心力效应的动力特性分析，并在此基础上进行了系统风致响应的时域分析.
数值算例分析表明，风力机叶片风场存在高频能量，塔架风谱更接近自然风谱；塔架叶片耦合模型以塔架变形伴随着叶片振动耦合振型为主，且叶片旋转引起的离心力作用不能忽略；叶片风致响应以低阶模态的共振效应为主，塔架的风致响应中高阶模态的贡献逐渐显著，忽略叶片影响带来的误差非常显著.综上所述，风力机的抗风设计应当考虑塔架-叶片的耦合效应，本文仿真算法可为风力机系统风场模拟和风致响应时域分析提供参考.
参考文献［1］ BUTTERFIELD S, MUSIAL W, SCOTT G. Definition of a 5MW reference wind turbine for offshore system development[M]. Colorado: National Renewable Energy Laboratory, 2009: 129-142.［2］ BAZILEVS Y, HSU M C, KIENDL J, et al. A computational procedure for prebending of wind turbine blades[J]. International Journal for Numerical Methods in Engineering, 2012, 89(3): 323-336.［3］ HOOGEDOORN E, JACOBS G B, BEYENE A. Aeroelastic behavior of a flexible blade for wind turbine application: a 2D computational study[J]. Energy, 2010, 35(2): 778-785.［4］ 金志昊, 范宣华, 苏先樾, 等. 风力机叶片顺风向风致振动研究[J]. 南京航空航天大学学报, 2011, 43(5): 677-681.JIN Zhihao, FAN Xuanhua, SU Xianyu，et al. Research on blade windinduced vibration in wind direction [J]. Journal of Nanjing University of Aeronautics & Astronautics, 2011, 43(5): 677-681.（In Chinese)［5］ 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.［6］ 李德源, 刘胜祥. 风波联合作用下的风力机塔架疲劳特性分析[J]. 太阳能学报, 2009,30(10):1250-1256．LI Deyuan, LIU Shengxiang. Fatigue characteristics analysis of wind turbine tower under windwave combined effect[J]. Acta Energiae Solaris Sinica, 2009,30(10):1250-1256．（In Chinese)［7］ NAESS A, GAIDAI O, HAVER S. Efficient estimation of extreme response of dragdominated offshore structures by monte carlo simulation [J]. Ocean Engineering, 2007, 34: 2188-2197．［8］ 柯世堂, 王同光, 赵林, 等. 风力机风振背景、共振响应特性及耦合项分析[J]. 中国电机工程学报, 2013, 33(26): 101-108．KE Shitang, WANG Tongguang, ZHAO Lin, et al. Background, resonant components and coupled effect of windinduced responses on wind turbine systems[J]. Proceeding of the CSEE, 2013,33(26): 101-108.（In Chinese)［9］ 廖明夫. 风力发电技术[M]. 西安： 西北工业大学出版社, 2008：60-72.LIAO Mingfu. Wind turbine technologies[M]. Xian:The Northwest Industry University Press, 2008：60-72.（In Chinese)［10］伍艳, 谢华, 王同光. 风力机叶片的非定常气动特性计算方法的改进[J]. 工程力学，2008, 25(10): 54-60．WU Yan, XIE Hua, WANG Tongguang. Modification of calculating unsteady aerodynamic characteristics of wind turbine blades[J]. Engineering Mechanics,2008, 25(10): 54-60.（In Chinese)［11］KAREEM A．Numerical simulation of wind effects：a probabilistic perspective[J]．Journal of Wind Engineering and Industrial Aerodynamics，2008(96)：1472-1497.［12］WANG T G, COTON F N. Prediction of the unsteady aerodynamic characteristics of horizontalaxis wind turbines including threedimensional effects[J]. Proceedings of the Institution of Mechanical Engineers Part A: Journal of Power and Energy, 2000, 214 (A5):385-400.
4风力机塔架叶片耦合模型风致响应
当平均风速过大时，风力发电机组会自动停止运行，这个临界风速称之为切出风速，对于兆瓦级风力机通常设定为25m/s左右.首先,进行切出风速下风力机塔架叶片耦合模型风致响应时域计算.限于篇幅，本文只给出了风力机耦合模型叶片尖部、机舱和塔架顶部位移响应时程曲线，如图7所示.
对比发现，基于风力机塔架叶片耦合模型时域计算得到的叶片尖部位移响应以低频共振效应为主，其脉动覆盖范围和数值最大，这是由于风力机遭受的风荷载卓越频率与叶片的摆振/舞动频率较为接近，从而引发了叶片的共振.机舱由于刚度较大，风荷载作用主要导致静力或准静力响应，因此位移响应时程曲线浮动范围较小，均值约在0.295m.塔架顶部的位移响应时程明显小于叶片响应数值，且高频振型贡献更加明显，对应的振型大多为叶片振动和塔架变形的耦合作用.分析表明，考虑了叶片的风场影响和结构模态耦合作用后，塔架的风致响应特性明显呈现多振型响应和多荷载形态.
时间/s(a)叶片尖部位移响应时程曲线
时间/s(b)机舱位移响应时程曲线
时间/s(c)塔架顶部位移响应时程曲线
图7风力机塔架叶片耦合模型典型部位位移响应时程曲线
Fig.7Thedisplacementtimehistoryontypicalnodes
ofwindturbinetowerbladecoupledmodel
为研究叶片对风力机塔架风致响应特性的影响，再对不同风速下考虑和不考虑叶片影响塔架风致响应进行分析.塔顶顺风向位移和塔底弯矩脉动响应均方根随风速的变化曲线如图8和图9所示，由图可见，在低风速(≤10m/s)下是否考虑叶片影响对风力机塔架风致响应影响不大，而随着风速的逐渐增大，考虑叶片影响计算得到的塔架风致响应要远大于不考虑叶片作用时的计算结果，其中考虑叶片影响时塔架顶部位移较未考虑的要大将近2倍，对于塔架底部弯矩来说要大1.5倍.
风速/（m•s-1）
图8不同风速下风力机塔架顶部顺风向位移响应均方差
Fig.8RMSofalongwinddisplacementresponses
oftowertopunderdifferentwindspeed
风速/（m•s-1）
图9不同风速下风力机塔底弯矩响应均方差
Fig.9RMSofbendingmomentresponsesoftower
bottomunderdifferentwindspeed
5结束语
基于风力机塔架叶片耦合模型，实现了给定风环境的系统风场模拟和考虑离心力效应的动力特性分析，并在此基础上进行了系统风致响应的时域分析.
数值算例分析表明，风力机叶片风场存在高频能量，塔架风谱更接近自然风谱；塔架叶片耦合模型以塔架变形伴随着叶片振动耦合振型为主，且叶片旋转引起的离心力作用不能忽略；叶片风致响应以低阶模态的共振效应为主，塔架的风致响应中高阶模态的贡献逐渐显著，忽略叶片影响带来的误差非常显著.综上所述，风力机的抗风设计应当考虑塔架-叶片的耦合效应，本文仿真算法可为风力机系统风场模拟和风致响应时域分析提供参考.
参考文献［1］ BUTTERFIELD S, MUSIAL W, SCOTT G. Definition of a 5MW reference wind turbine for offshore system development[M]. Colorado: National Renewable Energy Laboratory, 2009: 129-142.［2］ BAZILEVS Y, HSU M C, KIENDL J, et al. A computational procedure for prebending of wind turbine blades[J]. International Journal for Numerical Methods in Engineering, 2012, 89(3): 323-336.［3］ HOOGEDOORN E, JACOBS G B, BEYENE A. Aeroelastic behavior of a flexible blade for wind turbine application: a 2D computational study[J]. Energy, 2010, 35(2): 778-785.［4］ 金志昊, 范宣华, 苏先樾, 等. 风力机叶片顺风向风致振动研究[J]. 南京航空航天大学学报, 2011, 43(5): 677-681.JIN Zhihao, FAN Xuanhua, SU Xianyu，et al. Research on blade windinduced vibration in wind direction [J]. Journal of Nanjing University of Aeronautics & Astronautics, 2011, 43(5): 677-681.（In Chinese)［5］ 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.［6］ 李德源, 刘胜祥. 风波联合作用下的风力机塔架疲劳特性分析[J]. 太阳能学报, 2009,30(10):1250-1256．LI Deyuan, LIU Shengxiang. Fatigue characteristics analysis of wind turbine tower under windwave combined effect[J]. Acta Energiae Solaris Sinica, 2009,30(10):1250-1256．（In Chinese)［7］ NAESS A, GAIDAI O, HAVER S. Efficient estimation of extreme response of dragdominated offshore structures by monte carlo simulation [J]. Ocean Engineering, 2007, 34: 2188-2197．［8］ 柯世堂, 王同光, 赵林, 等. 风力机风振背景、共振响应特性及耦合项分析[J]. 中国电机工程学报, 2013, 33(26): 101-108．KE Shitang, WANG Tongguang, ZHAO Lin, et al. Background, resonant components and coupled effect of windinduced responses on wind turbine systems[J]. Proceeding of the CSEE, 2013,33(26): 101-108.（In Chinese)［9］ 廖明夫. 风力发电技术[M]. 西安： 西北工业大学出版社, 2008：60-72.LIAO Mingfu. Wind turbine technologies[M]. Xian:The Northwest Industry University Press, 2008：60-72.（In Chinese)［10］伍艳, 谢华, 王同光. 风力机叶片的非定常气动特性计算方法的改进[J]. 工程力学，2008, 25(10): 54-60．WU Yan, XIE Hua, WANG Tongguang. Modification of calculating unsteady aerodynamic characteristics of wind turbine blades[J]. Engineering Mechanics,2008, 25(10): 54-60.（In Chinese)［11］KAREEM A．Numerical simulation of wind effects：a probabilistic perspective[J]．Journal of Wind Engineering and Industrial Aerodynamics，2008(96)：1472-1497.［12］WANG T G, COTON F N. Prediction of the unsteady aerodynamic characteristics of horizontalaxis wind turbines including threedimensional effects[J]. Proceedings of the Institution of Mechanical Engineers Part A: Journal of Power and Energy, 2000, 214 (A5):385-400.
4风力机塔架叶片耦合模型风致响应
当平均风速过大时，风力发电机组会自动停止运行，这个临界风速称之为切出风速，对于兆瓦级风力机通常设定为25m/s左右.首先,进行切出风速下风力机塔架叶片耦合模型风致响应时域计算.限于篇幅，本文只给出了风力机耦合模型叶片尖部、机舱和塔架顶部位移响应时程曲线，如图7所示.
对比发现，基于风力机塔架叶片耦合模型时域计算得到的叶片尖部位移响应以低频共振效应为主，其脉动覆盖范围和数值最大，这是由于风力机遭受的风荷载卓越频率与叶片的摆振/舞动频率较为接近，从而引发了叶片的共振.机舱由于刚度较大，风荷载作用主要导致静力或准静力响应，因此位移响应时程曲线浮动范围较小，均值约在0.295m.塔架顶部的位移响应时程明显小于叶片响应数值，且高频振型贡献更加明显，对应的振型大多为叶片振动和塔架变形的耦合作用.分析表明，考虑了叶片的风场影响和结构模态耦合作用后，塔架的风致响应特性明显呈现多振型响应和多荷载形态.
时间/s(a)叶片尖部位移响应时程曲线
时间/s(b)机舱位移响应时程曲线
时间/s(c)塔架顶部位移响应时程曲线
图7风力机塔架叶片耦合模型典型部位位移响应时程曲线
Fig.7Thedisplacementtimehistoryontypicalnodes
ofwindturbinetowerbladecoupledmodel
为研究叶片对风力机塔架风致响应特性的影响，再对不同风速下考虑和不考虑叶片影响塔架风致响应进行分析.塔顶顺风向位移和塔底弯矩脉动响应均方根随风速的变化曲线如图8和图9所示，由图可见，在低风速(≤10m/s)下是否考虑叶片影响对风力机塔架风致响应影响不大，而随着风速的逐渐增大，考虑叶片影响计算得到的塔架风致响应要远大于不考虑叶片作用时的计算结果，其中考虑叶片影响时塔架顶部位移较未考虑的要大将近2倍，对于塔架底部弯矩来说要大1.5倍.
风速/（m•s-1）
图8不同风速下风力机塔架顶部顺风向位移响应均方差
Fig.8RMSofalongwinddisplacementresponses
oftowertopunderdifferentwindspeed
风速/（m•s-1）
图9不同风速下风力机塔底弯矩响应均方差
Fig.9RMSofbendingmomentresponsesoftower
bottomunderdifferentwindspeed
5结束语
基于风力机塔架叶片耦合模型，实现了给定风环境的系统风场模拟和考虑离心力效应的动力特性分析，并在此基础上进行了系统风致响应的时域分析.
数值算例分析表明，风力机叶片风场存在高频能量，塔架风谱更接近自然风谱；塔架叶片耦合模型以塔架变形伴随着叶片振动耦合振型为主，且叶片旋转引起的离心力作用不能忽略；叶片风致响应以低阶模态的共振效应为主，塔架的风致响应中高阶模态的贡献逐渐显著，忽略叶片影响带来的误差非常显著.综上所述，风力机的抗风设计应当考虑塔架-叶片的耦合效应，本文仿真算法可为风力机系统风场模拟和风致响应时域分析提供参考.
参考文献［1］ BUTTERFIELD S, MUSIAL W, SCOTT G. Definition of a 5MW reference wind turbine for offshore system development[M]. Colorado: National Renewable Energy Laboratory, 2009: 129-142.［2］ BAZILEVS Y, HSU M C, KIENDL J, et al. A computational procedure for prebending of wind turbine blades[J]. International Journal for Numerical Methods in Engineering, 2012, 89(3): 323-336.［3］ HOOGEDOORN E, JACOBS G B, BEYENE A. Aeroelastic behavior of a flexible blade for wind turbine application: a 2D computational study[J]. Energy, 2010, 35(2): 778-785.［4］ 金志昊, 范宣华, 苏先樾, 等. 风力机叶片顺风向风致振动研究[J]. 南京航空航天大学学报, 2011, 43(5): 677-681.JIN Zhihao, FAN Xuanhua, SU Xianyu，et al. Research on blade windinduced vibration in wind direction [J]. Journal of Nanjing University of Aeronautics & Astronautics, 2011, 43(5): 677-681.（In Chinese)［5］ 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.［6］ 李德源, 刘胜祥. 风波联合作用下的风力机塔架疲劳特性分析[J]. 太阳能学报, 2009,30(10):1250-1256．LI Deyuan, LIU Shengxiang. Fatigue characteristics analysis of wind turbine tower under windwave combined effect[J]. Acta Energiae Solaris Sinica, 2009,30(10):1250-1256．（In Chinese)［7］ NAESS A, GAIDAI O, HAVER S. Efficient estimation of extreme response of dragdominated offshore structures by monte carlo simulation [J]. Ocean Engineering, 2007, 34: 2188-2197．［8］ 柯世堂, 王同光, 赵林, 等. 风力机风振背景、共振响应特性及耦合项分析[J]. 中国电机工程学报, 2013, 33(26): 101-108．KE Shitang, WANG Tongguang, ZHAO Lin, et al. Background, resonant components and coupled effect of windinduced responses on wind turbine systems[J]. Proceeding of the CSEE, 2013,33(26): 101-108.（In Chinese)［9］ 廖明夫. 风力发电技术[M]. 西安： 西北工业大学出版社, 2008：60-72.LIAO Mingfu. Wind turbine technologies[M]. Xian:The Northwest Industry University Press, 2008：60-72.（In Chinese)［10］伍艳, 谢华, 王同光. 风力机叶片的非定常气动特性计算方法的改进[J]. 工程力学，2008, 25(10): 54-60．WU Yan, XIE Hua, WANG Tongguang. Modification of calculating unsteady aerodynamic characteristics of wind turbine blades[J]. Engineering Mechanics,2008, 25(10): 54-60.（In Chinese)［11］KAREEM A．Numerical simulation of wind effects：a probabilistic perspective[J]．Journal of Wind Engineering and Industrial Aerodynamics，2008(96)：1472-1497.［12］WANG T G, COTON F N. Prediction of the unsteady aerodynamic characteristics of horizontalaxis wind turbines including threedimensional effects[J]. Proceedings of the Institution of Mechanical Engineers Part A: Journal of Power and Energy, 2000, 214 (A5):385-400.
