柯世堂，曹九發(fā)，王 瓏，王同光

（南京航空航天大學(xué) 江蘇省風力機設(shè)計高技術(shù)研究重點實驗室，江蘇 南京 210016）

風荷載是風力機系統(tǒng)設(shè)計最重要荷載之一[1].隨著兆瓦級大功率風力機的普遍應(yīng)用，葉片的旋翼直徑達到百米量級，從而滿足捕捉風能的要求，但同時也遭受了強大風推力，并通過機艙將風荷載傳遞給塔架，使得葉片和塔架之間的耦合效應(yīng)愈加明顯.以往針對風力機系統(tǒng)的抗風研究，關(guān)注風力機空氣動力學(xué)性能的研究者大多忽略塔架的影響，對葉片的流場特性[2－3]和氣彈響應(yīng)[4－5]進行了深入詳細的研究；而關(guān)注風力機結(jié)構(gòu)受力性能的研究者大多以塔架為研究對象，把葉片遭受的風載等效為擬靜力施加到塔頂進行風振計算，或是建模時考慮塔架－葉片的耦合模型，而風場模擬時沒有完全考慮葉片平穩(wěn)風修正、旋轉(zhuǎn)效應(yīng)和氣彈效應(yīng)等影響[6－8].

目前，對風力機全機結(jié)構(gòu)的抗風設(shè)計，基本上采用等效風荷載的擬靜力分析方法，少量文獻進行了頻域內(nèi)的風振分析，考慮到塔架－葉片耦合模型脈動風場模擬的復(fù)雜性，在時域內(nèi)對風力機塔架－葉片耦合模型進行風致響應(yīng)計算則更為復(fù)雜.然而時域分析可以克服頻域分析中基于線性化假設(shè)的不足，同時時域計算結(jié)果的直觀性有利于工程師了解風力機系統(tǒng)的受力性能.因此，基于風力機塔架－葉片耦合模型的脈動風荷載模擬結(jié)果，基于塔架－葉片耦合有限元模型進行風振響應(yīng)時域分析，對風力機系統(tǒng)抗風的精細化設(shè)計有重要的指導(dǎo)作用和實用價值.

鑒于此，本文基于時域方法分析風力機塔架－葉片耦合模型在順風向風荷載動力作用下的響應(yīng)特性，建立了風力機塔架－葉片耦合的計算模型，運用諧波合成法和葉素－動量理論進行數(shù)值模擬獲得了考慮平穩(wěn)風修正、葉片旋轉(zhuǎn)效應(yīng)和相干性的塔架－葉片耦合模型氣動載荷，在此基礎(chǔ)上利用有限元方法進行了風力機系統(tǒng)風致響應(yīng)時域計算，基于計算結(jié)果，分析了風力機塔架－葉片耦合模型的風致響應(yīng)特性.

1 風力機塔架－葉片耦合模型風場模擬

由于風力機塔架－葉片耦合系統(tǒng)結(jié)構(gòu)本身的復(fù)雜性及隨機脈動風作用的不確定性，風場模擬時需作如下簡化：①只考慮順風向風荷載數(shù)值模擬和作用；②脈動風被假定為零均值的平穩(wěn)高斯隨機過程；③不考慮葉片和氣流之間的氣動彈性作用.

1.1 平穩(wěn)風修正

風力機的平穩(wěn)風速由于受到風切變、塔影效應(yīng)和上游風機尾流的影響，風場模擬時必須要對平穩(wěn)風模型進行修正[9].其中風剪切主要采用指數(shù)模型，塔影效應(yīng)主要采用適用于葉片在塔架上風向運行的潛流模型，上游尾流影響主要采用包含由于尾流引起的附加湍流的渦流粘度模型.

風剪切的影響是指平穩(wěn)風速隨著高度而變化，其指數(shù)修正模型表達公式如下：

式中：V（h）為指高度h處的風速；V（h0）為輪轂的參考高度h0處的參考風速；h0為輪轂的位置；當不考慮風剪切的影響時，可以將α的值設(shè)為0，取值一般為0.1～0.25.

塔影是由于風力機塔架的存在影響了風場平穩(wěn)風速，其影響修正主要有3種模型：葉片在塔架上風向運行的潛流模型、葉片在塔架下風向運行的經(jīng)驗?zāi)Ｐ秃蛢烧叩慕M合模型.現(xiàn)有風力機基本都是上風向風機設(shè)計，用以避免葉片周期地通過塔架尾跡產(chǎn)生的附加噪聲和激振力，因此采用潛流模型修正：

式中：DT為開始考慮塔影影響的高度處的塔架直徑；F為塔架直徑修正因子；z為計算點到塔架中心的縱向距離；x為風矢量經(jīng)過時距離塔架中心橫向距離.

由于風力機很少會存在單個運營狀態(tài)，基本都是以群體形式建立.因此，在模擬風力機模型的風場時，當由于單個風力機系統(tǒng)轉(zhuǎn)子部分或者全部處于上游風力機的尾流中時，就必須考慮上游風力機尾流對風力機風場的影響.本文采用渦流粘度模型來考慮尾流引起的附加湍流影響，修正公式如下：

式中：r為局部速度到尾流中心的距離；w為尾流寬度；V0為未受擾動平均風速；Δ為速度逆差.

1.2 脈動風模擬方法

已有研究表明，風力機葉片的風場模擬可采用葉素－動量理論[10]，相較渦尾跡方法和 CFD方法，葉素－動量理論可以實現(xiàn)葉片旋轉(zhuǎn)風場的模擬，且計算量較小、耗時較短，滿足風力機日?？焖儆嬎愕囊?考慮到葉素－動量理論是針對葉片在某一特定來流風速和工況下對應(yīng)參考風速的求解，本文通過諧波合成法[11]模擬獲得考慮葉片和塔架相干性的來流風速時程，再對每個來流風速樣本采用改進的葉素－動量理論進行計算得到該來流風速下葉片的參考風速，從而考慮葉片的旋轉(zhuǎn)效應(yīng)和相干性，如此迭代循環(huán)獲得葉片的參考風速時程.風力機塔架－葉片耦合模型脈動風模擬流程示意圖見圖1.

圖1 風力機塔架－葉片耦合模型風場模擬流程示意圖Fig.1 The flow sketch figure of wind field simulation for wind turbine tower－blade coupled model

風譜模型采用改進的Von karman模型，因為相比Davenport模型和Kaimal模型，改進的Von karman風譜模型能夠更加準確地產(chǎn)生隨機風場樣本，該模型校正了基本模型在高度150m以下的缺陷，更加符合風力機自身風場特性.圖2為3種風譜模型曲線，在低頻區(qū)3種風譜數(shù)值接近，在峰值附近區(qū)域，Kaimal和Von karman頻譜分別要低于和高于改進Von karman頻譜；而在高頻區(qū)域，Kaimal和Von karman頻譜幾乎重合，但是改進Von karman頻譜會略高于另外兩個頻譜幅值.

圖2 3種標準化風譜曲線示意圖Fig.2 The different longitudinal normal powerspectrum

根據(jù)風力機風場的自身特性，選擇改進的Von karman風譜模型：

式中：Suu為隨機風速變化自頻譜；f為頻率；σu為隨機風的標準差；~nu為無量綱頻率參數(shù)＝fLux／U10，U10為10m高基本風速，Lux為隨機風縱向分量長度尺度.

采用Davenport相關(guān)系數(shù)經(jīng)驗公式來考慮葉片和塔架之間的相關(guān)性：

式中：Cx，Cy和Cz分別為葉片和塔架上任意2點橫向、順風向和垂直向的衰減系數(shù)；ω為脈動風頻率；v（H）為H處平均風速.僅考慮順風向相關(guān)性時，其余兩個方向衰減系數(shù)均取0.

定義風機上n個風速模擬節(jié)點，假定均為零均值的平穩(wěn)高斯過程，其風譜密度函數(shù)矩陣為：

式中：sii（ω）為節(jié)點脈動風自功率譜，采用式（2）中的風譜模型計算；sij（ω）為互功率譜，其表達式需要用到塔架和葉片以及塔架上任意兩點之間的相干性，計算公式為：

其中葉片需考慮旋轉(zhuǎn)平面內(nèi)的各點間的相干性，以及葉片和塔架之間的相干性影響.再將S（ω）進行Cholesky分解，此時風力發(fā)電塔架上的任何一個節(jié)點脈動風速時程可以由其功率譜決定，根據(jù)Shinozuka理論，模擬的風速時程可以表達為：

式中：風譜在頻率范圍內(nèi)劃分成N個相同部分，Δω＝ω／N為頻率增量；｜Hjm（ωl）｜為基于 Davenport來流風譜矩陣進行Cholesky分解獲得的下三角矩陣的模；θml為介于0和2π之間均勻分布的隨機數(shù)，可采用Matlab的隨機數(shù)生成函數(shù)，建議每次生成隨機數(shù)后應(yīng)恢復(fù)初始狀態(tài)；ωl＝lΔω為頻域的遞增變量；ψjm（ωl）為兩個不同作用點之間的相位角，它是由Hjm（ωl）的虛部和實部的比值確定.

采用修正的BEM 理論[12]，引入葉根損失和葉尖損失，在軸向誘導(dǎo)因子較大時使用Ct的經(jīng)驗?zāi)Ｐ停⒓尤雱討B(tài)入流和動態(tài)失速模型.使用該方法，可以計算風力機在不同風速、轉(zhuǎn)速、槳距角及偏航角情況下的動態(tài)載荷.進而獲取作用在葉片上的脈動風速時程.具體方法如下：

根據(jù)BEM理論，葉片上的相對風速Vrel采用下式計算：

式中：vox，voy分別為沿順風向和橫風向的來流脈動風速，采用式（9）諧波疊加法計算；vrot為葉片旋轉(zhuǎn)導(dǎo)致的線速度；W為誘導(dǎo)速度；vbx和vby為葉片振動速度.

誘導(dǎo)速度W可由下式表示：

式中：B為葉片數(shù)；L為指升力；φ為入流角；ρ為空氣密度；r為葉片截面的展向位置；n為推力方向的單位向量；F為普朗特葉尖損失因子；fg為Glauert修正.同時，本文還采用了動態(tài)入流模型和動態(tài)失速模型，修正葉片運轉(zhuǎn)的非定常效應(yīng).

根據(jù)下式計算葉片攻角α：

式中：β為槳矩角；θtwist為葉片剖面幾何扭角，計算公式為：

1.3 數(shù)值模擬結(jié)果

本文以某5MW三槳葉變槳距特大型風力機系統(tǒng)為例進行風場數(shù)值模擬和風致響應(yīng)時域分析.塔高124m，底徑4.8m，頂徑2.6m，塔體通長為變厚度結(jié)構(gòu)，底壁厚150mm，頂壁厚60mm，通長厚度由底部至頂部呈線性減小趨勢.機艙長12m，寬4.6m，高4.2m，總質(zhì)量140.2×103kg.各槳葉間成120°夾角，沿周向平均分布，風輪直徑為120m，寬度2.4m，厚度0.38m，長度60m，偏航角為0°，額定轉(zhuǎn)速為17r／min.

基于上述方法，采用Matlab語言編制相應(yīng)的數(shù)值模擬程序MBEM－HSM.在計算過程中脈動風上限頻率取為2π，脈動風頻率分割點數(shù)取為2 048，頻率增量Δω＝0.003 07Hz.圖3和4給出了10m高平均風速為24m／s模擬得到的葉片和塔架中部的風速時程曲線和對應(yīng)功率譜曲線，為對比顯示風力機塔架－葉片耦合模型風場的特殊性，風譜坐標系均為對數(shù)坐標.

圖3 葉片參考風速時程模擬結(jié)果Fig.3 Simulating result of fluctuating wind velocity of blades

圖4 塔架中部風速時程模擬結(jié)果Fig.4 Simulating result of fluctuating wind velocity of tower middle part

對比發(fā)現(xiàn)，考慮塔架－葉片耦合模型葉片脈動風速功率譜曲線在高頻處存在明顯的能量浮動，應(yīng)該是由于考慮改進的Von karman風譜模型、三維旋轉(zhuǎn)和干擾效應(yīng)而產(chǎn)生的現(xiàn)象，在風力機系統(tǒng)風場模擬中應(yīng)引起重視；塔架中部的風速功率譜和Davenport風譜吻合較好，這是由于隨著高度的降低，塔架受到葉片的干擾效應(yīng)減小，在高頻出現(xiàn)的能量浮動和數(shù)值逐漸減小.根據(jù)本文模擬過程和對比分析可以認為，采用本文的MBEM－HSM方法可以很好地模擬風力機塔架－葉片耦合模型風場.

2 風力機塔架－葉片耦合模型動力特性

基于ANSYS軟件平臺，建立了風力機塔架－葉片耦合模型.其中葉片和塔體采用SHELL91單元，機艙及其內(nèi)部結(jié)構(gòu)可作為整體采用梁單元BEAM189模擬.通過多點約束單元耦合命令將各部分連接在一起，形成整體的風力機塔架－葉片耦合模型.依據(jù)效率和精度均衡的原則，模型一共劃分了3 812個單元.動力特性分析時把葉片旋轉(zhuǎn)產(chǎn)生的離心力作為預(yù)應(yīng)力預(yù)先均勻施加在葉片上，后續(xù)計算的頻率和模態(tài)信息均考慮葉片轉(zhuǎn)動帶來的離心力作用.圖5給出了考慮／不考慮葉片離心力作用2種工況下系統(tǒng)前200階自振頻率的分布情況.

圖5 風力機塔架－葉片耦合系統(tǒng)固有頻率分布圖Fig.5 Scattergram of natural frequence for wind turbine tower－blade coupled system

從圖5可以看出：1）風力機塔架－葉片耦合模型的基頻很小，僅為0.24Hz，第50階自振頻率為9.2 Hz，各模態(tài)之間的間隔很小，其自振頻率數(shù)值遠遠小于相同高度和直徑的高聳結(jié)構(gòu)，說明葉片和塔架的耦合效應(yīng)明顯降低了風力機系統(tǒng)的基頻，使得系統(tǒng)結(jié)構(gòu)可能會對風荷載更為敏感.2）當考慮葉片轉(zhuǎn)動引起的離心力作用時，耦合系統(tǒng)的基頻（0.27Hz）要略大于不考慮離心力作用下的系統(tǒng)基頻，并且隨著模態(tài)數(shù)目的增加，離心力效應(yīng)帶來的頻率影響越來越大.這是由于風力機塔架－葉片耦合系統(tǒng)的高階頻率主要以葉片的左右擺動和前后振動為主，此時葉片旋轉(zhuǎn)引起的離心力效應(yīng)更為明顯，相應(yīng)增強了風力機系統(tǒng)的剛度.本文的后續(xù)計算均采用考慮葉片離心力作用的更加真實的模態(tài)參數(shù).

從圖6給出的風力機塔架－葉片系統(tǒng)模態(tài)振型圖中可以看出，塔架的彎曲變形與葉片的揮舞／擺振相互耦合.通過多階模型分析表明，系統(tǒng)振型主要以葉片的前后揮舞和左右擺動為主，高階模態(tài)下塔架本身也會出現(xiàn)結(jié)構(gòu)變形和失穩(wěn)形態(tài)，并且和葉片變形耦合作用.動力特性分析表明，風力機系統(tǒng)動力響應(yīng)計算必須要考慮葉片和塔架耦合作用.

圖6 風力發(fā)電塔－輪系統(tǒng)的典型模態(tài)振型示意圖Fig.6 Typical modes of vibration on wind rotors systems

3 風致響應(yīng)時域分析方法

根據(jù)達朗貝爾原理，風力機系統(tǒng)結(jié)構(gòu)風致動力響應(yīng)控制方程為：

式中：p（t）為風荷載向量；M，C和K分別為結(jié)構(gòu)質(zhì)量、阻尼和剛度矩陣;（t），（t）和y（t）分別為節(jié)點加速度向量、速度向量和位移向量.

通過將節(jié)點位移向量從物理坐標轉(zhuǎn)換到廣義模態(tài)位移后，實現(xiàn)方程解耦，再引入模態(tài)阻尼，則各廣義模態(tài)對應(yīng)的運動方程為：

式中：qi（t）為第i階廣義模態(tài)位移響應(yīng)向量；mi，ci，ki和pi（t）分別為第i階模態(tài)的模態(tài)質(zhì)量、模態(tài)阻尼、模態(tài)剛度和模態(tài)力.這樣單自由度的運動方程可根據(jù)Duhamel積分原理，初始條件為零，數(shù)值解為：

式中：ωni＝（ki／mi）1／2為結(jié)構(gòu)的固有模態(tài)頻率；ζi＝bi／（2miωni）為模態(tài)阻尼比；ωdi＝ωni（1－ζ2）1／2為結(jié)構(gòu)阻尼振動頻率；Δτ為積分時間步長.本文選用采用Newmark－β逐步積分法和Newton－Raphson迭代理論，其逐步積分方法的基本假設(shè)為：

再將各模態(tài)下的廣義位移轉(zhuǎn)換為物理位移并進行疊加，可得到各節(jié)點風致響應(yīng)位移為：

式中：Φ為系統(tǒng)振型矩陣.

風力機塔架－葉片耦合模型風致響應(yīng)時域計算基于有限元軟件ANSYS平臺，將第2節(jié)數(shù)值模擬得到的風荷載時程作為外部激勵作用于風力機塔架－葉片耦合模型上，采用Newmark－β逐步積分法和Newton－Raphson迭代理論，其中各模態(tài)阻尼均為0.02，積分時間步長取為0.05s，響應(yīng)輸出時間步長取為0.025s，截取模態(tài)取系統(tǒng)的前50階.

4 風力機塔架－葉片耦合模型風致響應(yīng)

當平均風速過大時，風力發(fā)電機組會自動停止運行，這個臨界風速稱之為切出風速，對于兆瓦級風力機通常設(shè)定為25m／s左右.首先，進行切出風速下風力機塔架－葉片耦合模型風致響應(yīng)時域計算.限于篇幅，本文只給出了風力機耦合模型葉片尖部、機艙和塔架頂部位移響應(yīng)時程曲線，如圖7所示.

對比發(fā)現(xiàn)，基于風力機塔架－葉片耦合模型時域計算得到的葉片尖部位移響應(yīng)以低頻共振效應(yīng)為主，其脈動覆蓋范圍和數(shù)值最大，這是由于風力機遭受的風荷載卓越頻率與葉片的擺振／舞動頻率較為接近，從而引發(fā)了葉片的共振.機艙由于剛度較大，風荷載作用主要導(dǎo)致靜力或準靜力響應(yīng)，因此位移響應(yīng)時程曲線浮動范圍較小，均值約在0.295m.塔架頂部的位移響應(yīng)時程明顯小于葉片響應(yīng)數(shù)值，且高頻振型貢獻更加明顯，對應(yīng)的振型大多為葉片振動和塔架變形的耦合作用.分析表明，考慮了葉片的風場影響和結(jié)構(gòu)模態(tài)耦合作用后，塔架的風致響應(yīng)特性明顯呈現(xiàn)多振型響應(yīng)和多荷載形態(tài).

圖7 風力機塔架－葉片耦合模型典型部位位移響應(yīng)時程曲線Fig.7 The displacement time history on typical nodes of wind turbine tower－blade coupled model

為研究葉片對風力機塔架風致響應(yīng)特性的影響，再對不同風速下考慮和不考慮葉片影響塔架風致響應(yīng)進行分析.塔頂順風向位移和塔底彎矩脈動響應(yīng)均方根隨風速的變化曲線如圖8和圖9所示，由圖可見，在低風速（≤10m／s）下是否考慮葉片影響對風力機塔架風致響應(yīng)影響不大，而隨著風速的逐漸增大，考慮葉片影響計算得到的塔架風致響應(yīng)要遠大于不考慮葉片作用時的計算結(jié)果，其中考慮葉片影響時塔架頂部位移較未考慮的要大將近2倍，對于塔架底部彎矩來說要大1.5倍.

圖8 不同風速下風力機塔架頂部順風向位移響應(yīng)均方差Fig.8 RMS of along－wind displacement responses of tower top under different wind speed

圖9 不同風速下風力機塔底彎矩響應(yīng)均方差Fig.9 RMS of bending moment responses of tower bottom under different wind speed

5 結(jié)束語

基于風力機塔架－葉片耦合模型，實現(xiàn)了給定風環(huán)境的系統(tǒng)風場模擬和考慮離心力效應(yīng)的動力特性分析，并在此基礎(chǔ)上進行了系統(tǒng)風致響應(yīng)的時域分析.

數(shù)值算例分析表明，風力機葉片風場存在高頻能量，塔架風譜更接近自然風譜；塔架－葉片耦合模型以塔架變形伴隨著葉片振動耦合振型為主，且葉片旋轉(zhuǎn)引起的離心力作用不能忽略；葉片風致響應(yīng)以低階模態(tài)的共振效應(yīng)為主，塔架的風致響應(yīng)中高階模態(tài)的貢獻逐漸顯著，忽略葉片影響帶來的誤差非常顯著.綜上所述，風力機的抗風設(shè)計應(yīng)當考慮塔架－葉片的耦合效應(yīng)，本文仿真算法可為風力機系統(tǒng)風場模擬和風致響應(yīng)時域分析提供參考.

[1] BUTTERFIELD S，MUSIAL W，SCOTT G.Definition of a 5－MW reference wind turbine for offshore system development[M].Colorado：National Renewable Energy Laboratory，2009：129－142.

[2] BAZILEVS Y，HSU M C，KIENDL J，etal.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.Aero－elastic behavior of a flexible blade for wind turbine application：a 2D computational study[J].Energy，2010，35（2）：778－785.

[4] 金志昊，范宣華，蘇先樾，等.風力機葉片順風向風致振動研究[J].南京航空航天大學(xué)學(xué)報，2011，43（5）：677－681.JIN Zhi－h(huán)ao，F(xiàn)AN Xuan－h(huán)ua，SU Xian－yu，etal.Research on blade wind－induced 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] 李德源，劉勝祥.風波聯(lián)合作用下的風力機塔架疲勞特性分析[J].太陽能學(xué)報，2009，30（10）：1250－1256.LI De－yuan，LIU Sheng－xiang.Fatigue characteristics analysis of wind turbine tower under wind－wave 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 drag－dominated offshore structures by monte carlo simulation [J].Ocean Engineering，2007，34：2188－2197.

[8] 柯世堂，王同光，趙林，等.風力機風振背景、共振響應(yīng)特性及耦合項分析[J].中國電機工程學(xué)報，2013，33（26）：101－108.KE Shi－tang，WANG Tong－guang，ZHAO Lin，etal.Background，resonant components and coupled effect of wind－induced responses on wind turbine systems[J].Proceeding of the CSEE，2013，33（26）：101－108.（In Chinese）

[9] 廖明夫.風力發(fā)電技術(shù)[M].西安：西北工業(yè)大學(xué)出版社，2008：60－72.LIAO Ming－fu.Wind turbine technologies[M].Xi’an：The Northwest Industry University Press，2008：60－72.（In Chinese）

[10] 伍艷，謝華，王同光.風力機葉片的非定常氣動特性計算方法的改進[J].工程力學(xué)，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.（In Chinese）

[11] KAREEM A.Numerical simulation of wind effects：aprobabilistic 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 horizontal－axis wind turbines including three－dimensional effects[J].Proceedings of the Institution of Mechanical Engineers Part A：Journal of Power and Energy，2000，214（A5）：385－400.