風力機塔架葉片耦合模型風致響應時域分析
2014-09-27 14:17:18柯世堂曹九發王瓏王同光
湖南大學學報·自然科學版 2014年4期
柯世堂+曹九發+王瓏+王同光
文章編號：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.
