基于向量式有限元法的風力發(fā)電機組一體化仿真分析

陳俊嶺，陽榮昌，馬人樂

(1.同濟大學 建筑工程系，上海 200092; 2. 同濟大學 建筑設計研究院(集團)有限公司，上海 200092)

?

基于向量式有限元法的風力發(fā)電機組一體化仿真分析

陳俊嶺1?，陽榮昌2，馬人樂1

(1.同濟大學 建筑工程系，上海 200092; 2. 同濟大學 建筑設計研究院(集團)有限公司，上海 200092)

大型風電機組為強非線性剛柔耦合的周期時變多體系統(tǒng)，傳統(tǒng)有限元方法不能解決葉片轉動過程中由于剛體位移而導致的剛度矩陣奇異問題，向量式有限元法可有效考慮葉片的幾何非線性和大變形運動、塔架的彈性變形、氣動載荷等因素.基于MATLAB平臺編制空間梁系結構求解程序，選取懸臂梁受端部集中動荷載和歐拉梁繞定軸轉動兩個典型算例驗證程序的正確性.建立包含機艙、輪轂、葉片和塔架在內的風電機組一體化仿真分析模型，根據機組靜止狀態(tài)下自由振動響應，用模態(tài)參數識別方法提取了機組的自振頻率，計算結果與傳統(tǒng)有限元法吻合.采用諧波疊加法和本征正交分解法相結合的方式，運用B樣條曲面插值策略，模擬生成風電機組在正常運行狀態(tài)下的風速時程.采用向量式有限元法對正常運行狀態(tài)下的風電機組進行風振響應分析，較好地反映了重力對葉片內力周期性變化的影響以及葉片與塔架的共同作用.

風電機組；向量式有限元；仿真分析；運行狀態(tài)；非線性分析

風電機組運行時，作用在葉片上的空氣動力等外荷載使葉片和塔架產生振動，對風電機組關鍵部件造成疲勞損傷.因此建立包含塔架、葉片和機艙在內的風力發(fā)電機組一體化分析模型，對研究風力發(fā)電機組在不同運行狀態(tài)下的動力響應具有重要意義.Lobitz[1]提出了風力發(fā)電塔“質量-阻尼-彈簧”模型，通過連接矩陣實現塔體和葉片的耦合.王介龍等[2]綜合考慮葉片揮舞、擺振、扭轉和軸向位移4種變形運動與機艙剛體運動、塔架彈性變形的強非線性耦合作用，采用子空間迭代法計算在二維擬定常風荷載作用下系統(tǒng)的動力響應.柯世堂等[3]基于風力機塔架-葉片耦合模型，采用改進的葉素動量理論模擬了考慮平穩(wěn)風修正、葉片旋轉效應和空間相干性的風力機氣動載荷，并分析了氣動響應.Prowell等[4]將頂部葉片和機艙視為集中質量點建立簡化分析模型，通過與原型試驗結果進行對比，發(fā)現簡化分析模型對于整體結構的一階動力特性精度幾乎沒有影響.劉雄等[5]基于模態(tài)分析方法建立了風力機葉片和塔架的耦合動力學模型.李永建等[6]采用柔性多體動力學方法建立了風力機葉片-機艙-塔架耦合動力學方程.

風電機組為強非線性剛柔耦合的周期時變多體系統(tǒng)，采用傳統(tǒng)有限元方法不能解決葉片轉動過程中由于剛體位移而導致的剛度矩陣奇異問題，只能將葉片、輪轂和機艙簡化[7]；采用假設模態(tài)法不能考慮非線性問題，采用多體動力學方法進行建模無法體現結構的細部特征.而向量式有限元是基于向量力學與數值計算提出的一種新型數值計算方法[8-9]，可以有效處理連續(xù)體幾何變形、非線性與不連續(xù)的材料本構關系、多個連續(xù)體與剛體的運動及其相互耦合行為等各種復雜情形.本文采用向量式有限元法，基于MATLAB平臺編制空間梁系結構的求解程序，建了立包含機艙、輪轂、葉片和塔架在內的風電機組一體化仿真分析模型.

1 空間梁系結構向量式有限元法

向量式有限元是用一組空間點近似描述結構幾何，以牛頓第二定律描述質點的運動，模型離散為空間點以及各點之間的物理關系，在一系列持續(xù)增加的時間點內描述結構變化[10].向量式有限元的分析對象為質點，因此，首先要將結構離散成一系列的空間質點，質點之間通過結構單元連接(圖1).對于兩結點空間梁單元，可將兩端結點取為質點，并將桿件的質量mi凝聚在質點處.質點上的荷載包括作用在其上的外力、單元對質點的內力以及阻尼力，三者一起影響質點的運動軌跡.質點數量越多，計算結果也越精確，但計算量會相應增加.

圖1 向量式有限元示意Fig.1 Scheme of vector form intrinsic finite element

在向量式有限元中，稱時間段ta≤t≤tb內滿足標準化控制方程的計算單元為途徑單元[11].途徑單元是更新單元內力的基礎，當結構具有復雜行為過程時，可以通過引入足夠多的途徑單元使得時間段足夠小，則在這個時段內可認為構件的幾何變形很小，可以用大變位和小變形理論來處理.若以ta時刻單元的狀態(tài)為參考，只要求得內力增量，則tb時刻單元的內力便可獲得，求解內力增量的關鍵是計算單元在該時間段內的純變形.為求得單元在時間段內的純變形，可令單元作虛擬的剛體平動和剛體轉動(即逆向運動)[11].取ta時刻梁的構形為參考構形，單元從ta時刻起始位置IaJa運動到tb時刻的位置IbJb，其逆向運動即為：首先將IbJb平移使得質點Ib和質點Ia重合，接著再以Ia為中心轉動至與IaJa在同一直線到達虛擬位置I’J’，則IaJa和I’J’之間的形態(tài)差異即為純變形(圖2).將初始的位移與角度定義為向量，逆向平移與旋轉均可以通過向量運算來完成，求得單元在虛擬位置上的單元內力之后，再進行與逆向運動相反的正向運動恢復至當前時刻的位置，而在正向運動的過程中，平移不改變內力大小和方向，旋轉只改變內力的方向不改變大小，這樣就求得當前時刻單元的內力.

圖2 空間梁單元逆向運動Fig.2 Reverse movement of space beam-element

2 計算程序驗證

基于向量式有限元理論，采用MATLAB編制了風電機組整體動力學分析程序，主要受力部件用空間梁單元模擬.選取懸臂梁受端部集中力作用(動載)和歐拉梁繞定軸轉動2個典型算例驗證程序的正確性.

2.1 懸臂梁受端部集中力作用(動載)

(a) 相關計算參數

時間/s (b) 端點加速度響應圖3 懸臂梁端部突然加載Fig.3 Impulsive loading at the cantilever beam end

2.2 歐拉梁繞定軸轉動

圖4(a)為無限剛歐拉梁繞定軸轉動示意圖，若長度為l，單位長度質量為ρ，橫截面面積為A，根據理論力學可求得梁端無量綱位移u/l的解析解(見式(1)).傳統(tǒng)有限元法無法求解剛體轉動問題，而向量式有限元則可有效處理該類問題.模擬時共劃分5個單元，阻尼系數取0，時間步長取1×10-4s，梁端無量綱位移u/l的時程曲線見圖4(b).從圖中可看出，向量有限元計算結果和解析解吻合.

(a) 相關計算參數

時間/s (b) 端點位移響應圖4 歐拉梁繞定軸轉動情況Fig.4 Rotation of Euler beam around a fixed axis

(1)

3 風電機組整體動力學模型

3.1 模型建立

風電機組由葉片、塔架、機艙等部件組成，機艙內部構造雖然非常復雜但并不是本文關注的重點，因此本文在整體建模時將機艙簡化為一根剛性梁(圖5).輪轂和機艙在各自質心位置凝聚成質點，機艙質點和輪轂質點分別與塔架頂部質點、各葉根質點與輪轂質點、各葉根質點之間通過剛性梁連接.當風電機組處于運行狀態(tài)時，釋放塔架頂部質點繞X軸的轉動自由度使得葉片能夠旋轉；若機組處于停機狀態(tài)，只需取消自由度釋放即可；若需要模擬機組的偏航過程，只需釋放塔架頂部節(jié)點繞塔架軸線的轉動自由度.文中以某1.5 MW機組為例，主要參數如下：輪轂高度65 m，風輪直徑75 m；塔架底部直徑4.2 m，頂部直徑2.74 m；底部壁厚25 mm，頂部壁厚20 mm，其余部分基本從25 mm線性減小至12 mm；機艙總質量66 t，輪轂質量16 t；塔架材料為Q345鋼材.組裝后的風電機組整體模型如圖6所示.

圖5 風機各部件間的連接Fig.5 Connections among components of wind turbine

圖6 風電機組整體模型Fig.6 Integated model of wind turbine

3.2 關鍵參數取值

3.2.1 時間步長

運用中心差分法在求解結構非線性動力響應時的穩(wěn)定性受時間步長影響，時間步長Δt必須小于臨界值Δtcr(見式(2)).

(2)

式中：Tn為系統(tǒng)的最小自振周期，可用最小尺寸單元的最小自振周期代替[12].因此，在將模型空間離散時需要選擇合適的單元尺寸，則Δtcr根據單元的最小尺寸L按式(3)近似求得.

(3)

式中：E和ρ分別為材料的彈性模量和密度.當不考慮剛性梁、機艙和輪轂集中質量影響時，根據式(4)估算的Δtcr的量級為10-4s，最終選取時間步長Δt為4×10-5s.

3.2.2 阻尼系數

傳統(tǒng)有限元分析中采用瑞利阻尼假定 (見式(4))，但向量式有限元中無剛度矩陣的概念，因此需要將振型阻尼比等效成阻尼系數.

C=αM+βK

(4)

其中，C為阻尼矩陣；M和K分別為質量矩陣和剛度矩陣；α和β分別為質量阻尼系數和剛度阻尼系數.瑞利阻尼模型下第n階阻尼比ξn，可表述為式(5)[13].

(5)

當給出任意兩階振型的阻尼比時，α和β可根據式(5)得到的兩個聯(lián)立方程求得.根據文獻[14]的取值，葉片和塔架所有模態(tài)的阻尼比取1%.在計算葉片阻尼系數時ωm,ωn分別取葉片前兩階彎曲振動模態(tài)對應圓頻率，在計算塔架阻尼系數時ωm,ωn分別取塔架前兩階順風向模態(tài)對應圓頻率.由于在向量式有限元中不存在剛度矩陣，且β主要對高頻起作用，所以可通過適當提高α的方式考慮β的影響，具體取值可通過自由衰減數值試驗確定.

3.3 停機狀態(tài)下模態(tài)分析

采用傳統(tǒng)有限元法建立塔架的整體模型進行模態(tài)分析，前10階模態(tài)列于表1.

表1 傳統(tǒng)有限元和向量式有限元計算所得自振頻率

Tab.1 Natural frequencies by traditional finite element and vector form intrinsic finite element Hz

將有限元法的結點轉換成質點，保持單元結點的連接關系、材料及截面屬性不變，將模型轉換成向量式有限元模型，基于結構的響應通過峰值法識別結構自振頻率[15].通過在0.1 s內對輪轂處突然施加200 kN水平荷載，使機組在停機狀態(tài)自由振動，阻尼系數取0，提取典型質點的加速度時程用于模態(tài)識別.由于突加荷載使得機組產生前后方向的振動，只能識別前后振動對應的那部分模態(tài)的頻率，識別結果同樣列于表1.從表中可以看出，基于向量式有限元響應所識別的頻率和有限單元法計算結果也十分接近，若以有限單元結果為基準，最大誤差為1.6%，由此驗證了向量式有限元模型的正確性.

4 風振響應分析

風電機組在不同的運行狀態(tài)下因風輪的轉動所受風荷載并不相同，正常運行狀態(tài)下葉片的旋轉使得葉片上平均風隨時間變化，如果荷載不經處理而直接施加會造成突然加載的現象，使得機組的動力響應失真.為避免產生突然加載的現象，在荷載時程前續(xù)一段持時20 s的從0緩慢增加至實際值的時程.由于風荷載模擬的時間間隔和中心差分時間步長不一致，中間時刻的風速通過線性插值計算.風振響應分析時外部荷載條件為：輪轂高度處平均風速為額定風速11 m/s，葉片轉速為額定轉速17.4 r/min，考慮重力荷載.

4.1 風速時程模擬

采用諧波疊加法和本征正交分解法相結合的方式，運用B樣條曲面插值策略，模擬生成風電機組在正常運行狀態(tài)下的風速時程，主要步驟如下[16]：首先，確定風輪平面內和塔架上的采樣點分布形式及風速模擬的截止頻率、地貌等相關參數；然后，由諧波疊加法模擬采樣點處的脈動風速時程，并用本征正交分解法將風輪平面內采樣點的脈動風速時程進行分解；接著，不斷更新葉片位置，用B樣條曲面插值法求得葉片上插值點的本征模態(tài)值并根據本征正交分解法重構當前時刻的脈動風速，直到模擬時間結束；最后，將各個時刻的脈動風速匯總并和平均風疊加形成總的風速時程.風速模擬參數見表2，生成的輪轂質點和葉尖質點1處的風速時程見圖7.葉片上風荷載按照葉素動量理論方法計算，單管塔和機艙風荷載體型系數分別取0.6與1.2.

表2 風速模擬參數Tab.2 Parameters for wind simulation

4.2 風振響應

將外荷載施加于相應質點，采用向量式有限元法對機組正常運行狀態(tài)下的風振響應進行分析，提取葉片轉速達到目標值以后200 s的響應，塔架頂部加速度響應和塔架底部風輪平面外彎矩響應如圖8所示.從圖中可以看出，在正常運行階段，塔架頂部位移最大值為280 mm，約為塔架高度的1/232.根據VDI3834[17]規(guī)定：風電機組正常運行的機艙加速度均方根限值為0.3 m/s2，計算值為0.28 m/s2，主機廠要求正常運行階段的加速度最大值不超過0.1g，計算所得最大幅值為0.89 m/s2.

時間/s (a) 輪轂質點風速時程

時間/s (b) 葉尖1質點風速時程圖7 典型質點風速時程Fig.7 Time history at typical particles

時間/s (a)塔頂加速度時程

時間/s (b) 塔底彎矩圖8 塔頂加速度和塔底彎矩Fig.8 Top acceleration and base moment of the tower

以葉片1為例給出葉片根部的軸力、風輪平面外彎矩和剪力響應如圖9所示.從圖中可以看出，葉根的軸力和平面內剪力主要呈正弦式波動，這主要是由于葉片受重力荷載作用并周期性轉動所致.雖然面內剪力比面外幅值小，但是對葉片的疲勞壽命也有重要的影響，因此在進行風電機組風振響應分析時必須考慮重力的影響.

時間/s (a) 葉根軸力

時間/s (b) 葉根風輪平面內剪力

時間/s (c) 葉根風輪平面外剪力

時間/s (d) 葉根風輪平面外彎矩圖9 葉根內力Fig.9 Internal force at the root of blade

將典型響應變換至頻域以考察其響應譜特征(見圖10).圖10(a)為塔頂x向加速度響應譜，從圖中可以看出，塔頂加速度以共振響應為主且受到多個振型影響，主要以整體一階和葉片揮舞一階為主，高頻振型影響相對較小.圖10(b)為塔底彎矩響應譜，從圖中可以看出，塔底彎矩的共振響應主要以整體一階振型為主，葉片與塔架的相互作用也有一定影響.圖10(c)為葉尖x向位移響應，背景響應同樣占很大部分，而共振響應主要葉片揮舞一階為主，在1P和2P處有峰值.圖10(d)為葉根平面外彎矩響應譜，從圖中可以看出，葉根平面外彎矩的能量分布和葉尖x向位移響應譜相似，共振響應受到多個振型影響，1P和2P處也有較明顯的峰值，但主要以葉片揮舞一階為主.

頻率/Hz (a) 塔頂x向加速度響應譜

頻率/Hz (b) 塔底彎矩響應譜

頻率/Hz (c) 葉尖x向位移響應譜

頻率/Hz (d) 葉根平面外彎矩響應譜圖10 典型響應譜Fig.10 Spectra of typical dynamic responses

5 結 論

本文采用向量式有限元法建立了風力發(fā)電機組一體化仿真分析模型，分析了風電機組在靜止狀態(tài)下的自振特性和正常運行狀態(tài)下的風振響應.主要結論如下：

1)基于向量式有限元方法通過自由振動響應識別的機組頻率和基于傳統(tǒng)有限元法計算所得結果吻合較好，驗證了向量式有限元模型的準確性.

2)采用向量式有限元方法建模，可以較好反映重力對葉片內力周期性變化的影響，這種周期性變化對葉片的疲勞壽命不利，但對于塔架的影響相對較小.

3)正常運行狀態(tài)下塔頂位移和塔底彎矩的共振響應主要以整體一階振型為主，但葉片與塔架的相互作用也有一定影響；塔頂加速度響應則以共振響應為主且受到多個振型影響，其中整體一階和葉片揮舞一階貢獻較大.葉尖平面外位移和彎矩響應中背景響應都占很大部分，而共振響應主要以葉片揮舞一階為主；面內位移則主要受定軸轉動影響，在葉片轉動頻率處存在尖峰.

[1] LOBITZ D W. A nastran-based computer program for structural dynamic analysis of horizontal axis wind turbines[C]//Proceedings of the Horizontal Axis Wind Turbine Technology Workshop. Cleveland:Department of Energy and NASA-Lewis, 1984.

[2] 王介龍, 陳彥, 薛克宗. 風力發(fā)電機耦合轉子/機艙/塔架的氣彈響應[J] .清華大學學報:自然科學版, 2002, 2(2):211-215.

WANG Jie-long, CHEN Yan, XUE Ke-zong.Aeroelastic response analysis of coupled rotor/nacelle/tower for a horizontal axis wind turbine[J]. Journal of Tsinghua University :Science and Technology, 2002,2(2): 211-215. (in Chinese)

[3] 柯世堂,曹九發(fā),王瓏,等. 風力機塔架-葉片耦合模型風致響應時域分析[J].湖南大學學報:自然科學版, 2014, 41(4):87-93.

KE Shi-tang, CAO Jiu-fa, WANG Long,etal. Time-domain analysis of the wind-induced responses of the coupled model of wind turbine tower-blade coupled system[J]. Journal of Hunan University:Natural Sciences, 2014, 41(4):87-93.(In Chinese)

[4] PROWELL I, VELETZOSM, ELGAMALA,etal. Experimental and numerical seismic response of a 65 kW wind turbine[J]. Journal of Earthquake Engineering, 2009, 13:1172-1190.

[5] 劉雄, 張憲民, 陳嚴,等. 水平軸風力機結構動力響應分析[J]. 太陽能學報, 2009, 30(6):804-809.

LIU Xiong, ZHANG Xian-min, CHEN Yan,etal. Structure dynamic response analysis of horizontal axis wind turbines[J]. ACTA Energiae Solaris Sinica, 2009, 30(6):804-809.(In Chinese)

[6] 李永建, 殷玉楓,丁健剛,等. 隨機風載荷下大型風力機葉片/機艙/塔架耦合動力學分析[J]. 可再生能源, 2014, 32(7):992-997.

LI Yong-jian,YIN Yu-feng,DING Jian-gang,etal. Coupling dynamics analysis on blades/nacelle/tower of large-scalewind turbine with stochastic wind load[J]. Renewable Energy Resources, 2014, 32(7):992-997.(In Chinese)

[7] 陳俊嶺，陽榮昌，馬人樂. 大型風電機組組合式塔架結構優(yōu)化設計. 湖南大學學報:自然科學版, 2015, 42(5):29-35.

CHEN Jun-ling, YANG Rong-chang, MA Ren-le. Design optimization of wind turbine tower with composite structure using particle swarm approach[J]. Journal of Hunan University:Natural Science, 2015, 42(5):29-35. (In Chinese)

[8] TING E C, SHI H C, WANG Y K. Fundamentals of a vector form intrinsic finite element: Part 1.Basic procedure and a plane frame element [J].Journal of Mechanics, 2004, 20(2):113-122.

[9] TING E C, SHI H C, WANG Y K. Fundamentals of a vector form intrinsic finite element: Part 2.Plane solid element [J].Journal of Mechanics,2004,20(2):123-132.

[10]朱明亮，董石麟. 基于向量式有限元的弦支穹頂失效分析[J]. 浙江大學學報:工學版，2012，46(9):1611-1618, 1632.

ZHU Ming-liang, DONG Shi-lin. Failure analysis of suspen-dome by vector form intrinsic finite element method[J]. Journal of Zhejiang University:Engineering Science, 2012，46(9):1611-1618, 1632.(In Chinese)

[11]喻瑩,羅堯治.基于有限質點法的結構屈曲行為分析[J].工程力學,2009,26(10):23-29.

YU Ying, LUO Yao-zhi. Buckling analysis of structures by the finite particle method[J]. Engineering Mechanics,2009,26(10):23-29.(In Chinese)

[12]王勖成.有限單元法[M].北京:清華大學出版社,2003:617-662.

WANG Xu-cheng. Finite element method[M]. Beijing: Tsinghua University Press,2003:617-662. (In Chinese)

[13]CLOUGH R W,PENZIEN J. Dynamics of structures[M]. Berkeley, Calif:Computers and Structures Inc, 2004: 234-245.

[14]WOUDEV C, NARASIMHAN S. A study on vibration isolation for wind turbine structures[J]. Engineering Structures, 2014, 60: 223-234.

[15]馬人樂, 馬躍強, 劉慧群, 等. 風電機組塔筒模態(tài)的環(huán)境脈動實測與數值模擬研究[J]. 振動與沖擊, 2011, 30(5):152-155.

MA Ren-le,MA Yue-qiang,LIU Hui-qun，etal. Ambient vibration test and numerical simulation for modes of wind turbine towers[J]. Journal of Vibration and Shock, 2011, 30(5):152-155.

[16]馬人樂，陽榮昌，陳俊嶺. 本征正交分解法在風電機組風場模擬中的應用[J]. 太陽能學報, 2014,35(9): 1764-1770.

MA Ren-le, YANG Rong-chang, CHEN Jun-ling. Application of proper orthogonal decomposition method in wind field simulation for wind turbine system[J]. ACTA Energiae Solaris Sinica, 2014, 35(9):1764-1770. (In Chinese)

[17]Measurement and evaluation of the mechanical vibration of wind energy turbines and their components VDI 3834[S]. Düsseldorf: Verein Deutscher Ingenieure, 2009:8-12.

Integrated Simulation of Wind Turbine Based on Vector Form Intrinsic Finite Element

CHEN Jun-ling1?, YANG Rong-chang2, MA Ren-le1

(1. Dept of Structural Engineering, Tongji Univ, Shanghai 200092, China; 2.Tongji Architectural Design (Group) Co Ltd, Shanghai 200092, China)

Large wind turbine system is a periodic time-varying system with rigid-flexible coupling multi-bodies. The traditional finite element method cannot solve the singular stiffness matrix produced by the rotation of blades. However, the vector form intrinsic finite element method can effectively solve the geometric deformation of elastic continuum, the nonlinear or discrete constitutive model, the coupling motion of continuum and rigid body, and so on. In this study, a solver program of space beam elements was developed by the vector form intrinsic finite element method with MATLAB code, and verified by two typical examples where one is a cantilever method with a dynamic force acting at the end, and the other is an Euler beam rotating around a fixed axis. The integrated simulation of wind turbine system that consists of tower, rotor blades, and nacelle was then established, and its dynamic response of free vibration under parking was analyzed. The natural frequencies of the turbine system were obtained by modal parameter identification, and they agree well with the results obtained by traditional finite element method. Moreover, the weighted amplitude wave superposition method and proper orthogonal decomposition method as well as B-spline surface interpolation were employed to obtain the wind time series of wind turbine under the normal operation condition. The wind-induced dynamic response of wind turbine system was also calculated by vector form intrinsic finite element method. The results reflect the periodic influence of gravity on the internal forces of blades and the interaction between blades and the tower.

wind turbine; vector form intrinsic finite element; numerical simulation; operation condition; nonlinear analysis

1674-2974(2016)11-0141-08

2015-11-24

國家自然科學基金資助項目(51378381), National Natural Science Foundation of China(51378381)

陳俊嶺(1974-)，女，河北景縣人，同濟大學副教授，工學博士?通訊聯(lián)系人，E-mail:chenjl@#edu.cn

TK83；O322

A