電壓跌落對風力機葉片振動影響分析

尹堯杰， 廖明夫， 呂　品， 劉展翅, 王四季

(西北工業大學 動力與能源學院，西安　710029)

電壓跌落對風力機葉片振動影響分析

尹堯杰， 廖明夫， 呂品， 劉展翅, 王四季

(西北工業大學 動力與能源學院，西安710029)

摘要：風力機單機容量不斷增大，裝機容量不斷增長，風力機與電網之間的耦合效應越發突出。電網的干擾以及故障會對風力機機械部件產生沖擊。為了獲得葉片在電網電壓跌落時的動態響應，把葉片簡化為懸臂梁，用3D Timoshenko梁進行離散化建模，建立了葉片運動的動力學方程，基于某1.5 MW風力機低電壓穿越測試得到的轉速，變槳曲線，獲得葉片所受時變載荷的基礎上，運用Newmark法對葉片在電壓三相對稱跌落與不對稱跌落下的動態響應進行了計算。結果表明電壓的跌落會對葉片產生沖擊效應，在葉片拍打方向振動幅值比揮舞方向更大。葉片振動變化的總體特征隨著風機載荷和電壓跌落的幅值增大而增大的。同等工況和跌落幅度下，三相對稱跌落時振動的變化要大于不對稱跌落的振動變化。

關鍵詞：電壓跌落;風力機;葉片;Timoshenko梁模型;動態響應

隨著風力發電技術快速發展，風力機單機容量不斷增大，風電裝機容量高速增長，風能已經成為重要的能源供應選擇[1]。

圖1　風力機與電網之間的耦合[3]Fig.1 Interaction between wind turbine and grid

過去風力機的設計要求主要關注風力機自身的安全性能，并沒有考慮其對電網的影響[2]。然而，隨著風力機裝機容量的大規模增加，風力機與電網的耦合性增強，如圖1所示，風力機的動態特性影響并網電能質量，而電網的干擾以及故障也會對風力機電氣部件以及機械構件產生影響。

許多國家修訂與更新了風力機并網要求[4]，其中的一個重要方面就是要求風力機具備低電壓穿越能力，即當電網故障或擾動引起風電場并網點電壓跌落時，在一定的電壓跌落范圍內，風電機組能夠不間斷并網運行[5]。如圖2所示為中國國家電網對風力機低電壓穿越的要求，曲線上方區域的跌落不允許風機脫網。

圖2　風力機低電壓穿越要求[6]Fig.2 Wind turbine LVRT requirements

低電壓穿越能力在風力機的設計標準中也有體現。例如IEC 61400-1在風電場外部電網條件評估中規定為保證電氣設計條件的兼容性，外部電網條件中要考慮當地電網對低電壓穿越特性的要求[7]。

電網電壓跌落引發風力機機械、電氣功率的不平衡，使得機組載荷發生振蕩，嚴重時導致機組切除。電壓跌落導致風力發電機組脫網的原因不僅與電氣部件有關，也可能與機械振動或附屬系統的低電壓承受能力有關。因此在風力發電機組進行低電壓測試時，也要進行風電機組載荷相關測試[8]。在風力機的設計中需要評估電壓跌落對風機各個部件動力學的影響。

雙饋風力發電系統如圖3所示作為一種主流機型采用齒輪箱升速，通過控制轉子勵磁電流的頻率，使定子頻率恒定，實現變速恒頻發電[9]。雙饋風機定子側直接聯接電網，這種直接耦合使得電網電壓跌落直接反應在電機定子端電壓上，使得其對電網故障比較敏感，雙饋風力發電系統的低電壓穿越特性也成為了研究熱點。

圖3　雙饋風力發電系統Fig.3 Double-fed wind power generation system

風力機主要結構部件中，葉片是個展向長，弦向短的柔性體，很容易發生振動與變形[10-11]。葉片結構的變形和振動，不僅會對葉片自身的氣動性能產生影響，也會對風機電能質量產生影響。

因此本文建立了3D Timoshenko梁的風力機葉片有限元模型，基于在某風場對某型1.5 MW雙饋風力發電機組低電壓穿越特性的測試數據[12]，計算與分析了對稱電壓跌落與不對稱電壓跌落對風力機葉片動態特性的影響。

1風力機葉片建模

1.1風力機坐標系統

圖4　葉片坐標系Fig.4 Coordinate of the blade

本文采用德國GL[13]認證所定義的坐標系作為葉片計算的坐標系，整個坐標系統如圖4所示，其中ZR沿葉片展長方向，XR垂直于ZR指向塔架方向，YR由右手規則確定，垂直于機艙主軸。

1.2葉片有限元模型

葉片展弦比很大，可以將葉片簡化為懸臂梁[14]，大型風力機的葉片越來越長，柔性增強，空腹特性突出，使得剪切變形不能忽略，因此采用考慮剪切變形的2節點3D Timoshenko梁單元進行離散化，離散化后的葉片如圖5所示。每個單元有2個節點，每個節點有6個自由度(3個位移，3個轉角)。

圖5　葉片離散化簡圖Fig.5 Discretization of the blade

1.3梁單元的剛度矩陣與質量矩陣

葉片離散化后，需要建立梁單元的剛度矩陣與質量矩陣?？紤]剪切變形影響時，梁的法向位移(撓度)可以表示為兩部分的疊加。

ω=ωb+ωs

(1)

式中，ωb是彎曲變形的法向位移；ωs是剪切變形的法向位移。

(2)

式中，Nb為Hermite插值函數；Ns為拉格朗日插值函數。

將式(2)代入梁彎曲最小勢能原理的泛函式(3)

(3)

從δΠ=0可推導出Timoshenko梁的剛度矩陣[15]，具體表達形式見附錄。

(4)

質量矩陣可以根據Hamilton原理推導，考慮梁的平移與旋轉具體推導過程可參考文獻[16]。

平移質量矩陣：

(5)

旋轉質量矩陣：

(6)

式中，l為單元長度；ρ為材料密度；A為截面面積；Ir為回轉半徑；Nv,Nθ為插值形函數；ξ為無量綱局部坐標。

單元剛度矩陣為：

M=[Mt]e+[Mr]e

(7)

單元質量矩陣為對稱矩陣，具體表達形式見附錄。

(8)

圖6　葉片截面圖Fig.6 Section of the blade

葉片抗彎剛度一般按截面主慣性軸(Xe-Ye)給出，需要分解到總體坐標系(Xs-Ys)上如圖6所示，公式按式(9)給出[17]。

EIx=EIyecos2β-EIxesin2β

EIy=EIyesin2β+EIxecos2β

(9)

式中：Ixe,Iye為主慣性軸上的截面慣性矩；β為葉片截面扭角；E為材料彈性模量。

葉片旋轉過程中產生的離心力會影響結構的剛度，產生離心剛化作用，影響葉片結構動力特性。離心剛化作用產生的單元應力剛度矩陣[18]為：

(10)

式中：

(11)

(12)

(13)

式中：F為作用于單元的軸向力；L為單元長度；Jz為極慣性矩。

將剛度矩陣與單元應力剛度矩陣相加得到單元組合剛度矩陣即：

Ke=K+Kσ

(14)

2葉片運動方程

組裝單元的剛度矩陣與質量矩陣可以得到葉片整體的剛度矩陣與質量矩陣。根據最小勢能原理可有導出葉片在時變載荷下的動力學運動方程為：

(15)

2.1空氣動力學載荷

風力機氣動性能計算方法采用葉素動量方法(BEM)，它綜合了動量理論和葉素理論，考慮了葉尖損失、葉柵效應、失速修正、間隙修正及偏航角等因數的影響，可以比較精確地計算風輪轉子的氣動性能，所以在風力機設計和氣動計算中得到廣泛的應用。另外，還考慮了風剪、偏航、風輪的結構參數和風力機安裝參數等[19]。在動態工況下，必須考慮動態入流與動態失速的影響對動量葉素-理論進行修正，具體的求解方法請參考文獻[20]。如圖7所示，葉片半徑r的位置，來流風速為U∞，葉片旋轉角速度為Ωr。軸向和切向誘導因子分別為a，a′。氣流軸向速度為U∞(1-a)，氣流的切向速度Ωr(1+a′)，合成的相對速度為W，由于風速方向的改變，將導致攻角α改變，定義來流角φ=α+β，β為安裝角和截面扭角之和。L和D是葉素的升力和阻力。

圖7　葉片的氣動力Fig.7 Aerodynamic loads of the blade

應用葉素動量理論，考慮葉尖與輪轂損失F，速度誘導因子求解公式為[21]，

如果a≤ac時，

(16)

如果a>ac時

(17)

式中，ac≈0.2

F為損失因子，表達式為，

(18)

σ為實度系數，表達式為，

(19)

式中，c為翼型弦長。

Cn為法向力系數，表達式為，

Cn=CLcosφ+Cdsinφ

(20)

式中，CL,Cd為截面升力系數與阻力系數。

(21)

經過迭代可以求得葉素軸向與切向誘導因子。根據葉素-動量理論，單位長度的葉片在推力方向(X)與切向方向(Y)的力為:

(22)

式中：CL,Cd為截面升力系數與阻力系數；W為氣流相對速度；φ為截面入流角；ρ為空氣密度；c為翼型弦長。

2.2重力載荷

葉片旋轉中，始終受到重力的作用，

重力在三個方向的分量表達式為：

(23)

2.3離心力載荷

葉片旋轉過程中產生的離心力沿葉片展向向外，離心力的表達式可以表示為：

(24)

式中，ω為葉片旋轉角速度。

根據求解得到的空氣動力學載荷、重力載荷、離心力載荷可以得到單位長度葉片所受的載荷，

(25)

3葉片固有屬性計算

采用建立的有限元模型對某1.5 MW的葉片動態特性進行計算，葉片R=35 m，其中輪轂半徑Rhub=1.75 m。計算得到的自振頻率與NREL軟件BModes[22]進行對比如表1所示，可以看出對于葉片常關注的前幾階自振頻率計算的誤差很小,驗證所建的Timoshenko梁模型的正確性。

表1　計算自振頻率與BModes對比

4葉片的動態響應

中國電科院在其低電壓穿越能力測試規程[23]規定，當風力發電有功功率輸出分別在以下范圍內時，測試風力發電機組在電網電壓跌落時的響應特性：① 大功率輸出0.9Pn之上;② 小功率輸出0.1Pn到0.3Pn之間。響應測試內容既包括電壓，有功功率，無功功率，電壓，電流等電學參數，還需要測試風力機的運行參數，例如轉速，變槳角的時間序列。規定的測試電網電壓跌落規格如表2所示。

表2　電網電壓跌落規格

對某型1.5 MW低電壓穿越測試得到的轉速與變槳響應曲線[12]，加載到所建立有葉片限元模型中，計算得到葉片所受的時變力，采用Newmark法積分可以得到葉片在低電壓穿越時的動態響應。

由于振動的變化具有一定的規律性，因此在進行不同機組工況振動分析時，只對同一機組工況的兩個極端電壓跌落情況(最大跌落和最小跌落)進行對比分析，如表3所示。

表3　低電壓穿越測試工況

4.1工況1與工況3對比分析

對比分析：工況1與工況3跌落幅值與跌落時間相同，只有機組的工況不同，工況1機組額定狀態，工況3機組偏載狀態。圖8與圖9可以看出，工況1機組工作在額定狀態，電壓的跌落引起了轉速的波動，轉速控制與變槳控制同時作用，而工況3機組工作在偏載狀態，轉速波動，并未變槳。

圖10與圖11可以看出平面外即葉片拍打方向振動的幅值比平面內即揮舞振動的幅值大。沖擊的幅值即電壓跌落時的振動到電壓恢復時的振動，工況1在拍打方向是工況3的2倍，揮舞方向為1.5倍，跌落影響逐漸變小后，工況1中葉片的振動仍然為工況3的1.5倍左右。

4.2工況2與工況4的對比

對比分析：工況2與工況4跌落幅值與跌落時間相同，只有機組的工況不同，工況2機組額定狀態，工況4機組偏載狀態。圖12與圖13可以看出，工況2機組工作在額定狀態，電壓的跌落引起了轉速的波動，轉速控制與變槳控制同時作用，而工況4機組工作在偏載狀態，轉速波動，并未變槳。

圖14和圖15可以看出葉片的沖擊效應很明顯，葉片拍打方向的振動幅值比揮舞方向的振動大。沖擊的幅值工況2在拍打方向是工況4的3.5倍，揮舞方向為4倍，跌落影響逐漸變小后，振動仍然為工況4的2倍。

4.3工況2與工況6的對比

對比分析：工況2與工況6機組工況，跌落時間相同，電壓跌落的形式不同，工況2是電壓對稱跌落，工況6是不對稱稱跌落。從圖12和圖16可以看出兩組工況電壓的跌落使得轉速控制與變槳控制同時作用。

從圖17與圖18可以看出，兩組工況振動規律相似，都有明顯的沖擊現象，同樣的葉片拍打方向的振動要比揮舞方向大。工況2比工況6在平面外方向的沖擊幅值高出約12%，平面內方向高出約13%。

4.4工況1與工況2的對比

對比分析：工況1與工況2機組工況相同，電壓跌落的幅值與跌落的時間不同。從圖8和圖12中可以看出，兩組工況轉速與變槳發生振蕩。

從圖19和圖20可以看出， 沖擊幅值在平面外方向工況2是工況1的1.8倍，平面內方向為2倍。

所有工況中工況2的電壓跌落沖擊效應最明顯。圖21，圖22為葉尖振動的速度與加速度，可以看出振動變化很大。

圖8　工況1的電壓跌落曲線、轉速與變槳角響應曲線Fig.8 Curve of voltage drop,speed and pich of test 1

圖9　工況3的電壓跌落曲線轉速與變槳角響應曲線Fig.9 Curve of voltage drop,speed and pich of test 3

圖10　葉尖平面外方向響應Fig.10 Dynamic response of blade tip in flap wise

圖11　葉尖平面內方向響應Fig.11 Dynamic response of blade tip in edgewise

4.5低穿過程中葉片振動的機理解釋

當電網電壓發生跌落時，發電機端的電磁轉矩會瞬間變小，使得氣動轉矩與發電機轉矩產生不平衡，葉輪轉速增加，如果機組工作在偏載工況下，風力機組只有轉速控制，如果機組工作在滿載工況下，當轉速瞬間發生較大的波動時，轉速控制與變槳控制同時作用。

如上文葉片載荷分析，可知作用在葉片上的力，尤其是氣動力與轉速以及變槳角密切相關，并且是強烈的非線性關系，轉速與變槳角波動，引起作用于葉片上的氣動力發生波動，進而激起葉片發生振動。

另外從結構特性分析，當變槳角發生改變時候，葉片揮舞與拍打方向的剛度分布會發生改變，使得葉片的自振頻率發生改變，如圖23所示，也是引起葉片振動的原因。

圖12　工況2的電壓跌落曲線轉速與變槳角響應曲線Fig.12 Curve of voltage drop,speed and pich of test 2

圖13　工況4的電壓跌落曲線轉速與變槳角響應曲線Fig.13 Curve of voltage drop,speed and pich of test 2

圖14　葉尖平面外方向響應Fig.14 Dynamic response of blade tip in flap wise

圖15　葉尖平面內方向響應Fig.15 Dynamic response of blade tip in edgewise

圖16　工況6的電壓跌落曲線 轉速與變槳角響應曲線Fig.16 Curve of voltage drop,speed and pich of test 2

圖17 葉尖平面外方向響應Fig.17Dynamicresponseofbladetipinflapwise圖18 葉尖平面內方向響應Fig.18Dynamicresponseofbladetipinedgewise圖19 葉尖平面外方向響應Fig.19Dynamicresponseofbladetipinflapwise

圖20 葉尖平面內方向響應Fig.20Dynamicresponseofbladetipinedgewise圖21 葉尖振動速度Fig.21Vibrationofbladetip圖22 葉尖振動加速度Fig.22Vibrationaccelerationofthebladetip

圖23　變槳角對葉片自振頻率的影響Fig.23 The impact of pitch angle on the blade natural frequency

葉片的振動會與氣動力形成耦合，使得葉片的振動超標引起機組的振動報警，影響機組的安全穩定運行。

5結論

通過以上的建模計算，結果的對比與分析可以得出以下結論：

(1) 風力機在低電壓穿越過程中，由于載荷發生振蕩，葉片的振動具有明顯的沖擊。其中葉片拍打方向振動幅值更明顯，葉片拍打方向的載荷變化會引起風力機整機軸向載荷的變化，影響機組的安全穩定運行。

(2) 葉片振動變化的總體特征是隨著電壓跌落的幅值增大而增大的。同等工況下葉片的振動在電壓跌落90%工況下，振動變化不是很明顯，在電壓跌落20%工況下，出現較大的沖擊振動，其沖擊振動是90%跌落的2倍。

(3) 葉片振動變化隨著風機載荷的增大而增大。跌落90%時額定工況是偏載工況的2倍，跌落20%時額定工況是偏載工況的3.5倍.電壓跌落影響逐漸消失，振動仍然是分別跌落前的1.5倍與2倍。

(4) 同等工況和跌落幅度下，三相對稱跌落時振動的變化要比不對稱跌落的振動大10%左右。

參 考 文 獻

[1] Gasch R,Twele J.Windkraftanlagen: Grundlagen, Entwurf, Planung und Betrieb[M].Berlin：Springer Vieweg，2013.

[2] Fadaeinedjad R,Moschopoulos G,Moallem M.Investigation of voltage sag impact on wind turbine tower vibrations[J].Wind Energy,2008，11(4):351-375.

[3] Burton T,Jenkins N,Sharpe D,et al. Wind energy handbook[M].2nd ed.Chichester:John Wiley&Sons，2011.

[4] Hansen A D,Cutululis N A,Markou H,et al.Grid fault and design-basis for wind turbines—Final report[R]. Risφ-R-1714, Roskilde:RisφNational Laboratory for Sustainable Energy Technical University of Denmark，2010.

[5] 廖明夫,Gasch R,Twele J.風力發電技術[M].西安：西北工業大學出版社，2009.

[6] 國家電網公司風電場接入電網技術規定[S].北京：國家電網中國電力科學研究院，2009.

[7] IEC 61400-1.Wind turbine-part1:design requirements[S]. Geneva: International Electrotechnical Commission,2005.

[8] IEC-61400-21-2008. Part 21: Measurement and assessment of power quality characteristics of grid connected wind turbines[S]. Geneva: International Electrotechnical Commission,2008.

[9] Heier S.Grid integration of wind energy conversion systems[M]. 2nd ed. Chichester:John Wiley&Sons, 2006.

[10] 劉雄,李剛強,陳嚴,等.水平軸風力機葉片動態響應分析[J].機械工程學報，2010，46(12)：128-134.

LIU Xiong, LI Gang-qiang,CHEN Yan,et al. Dynamic response analysis of the blade of horizontal axis wind turbines[J]. Journal of Mechanical Engineering, 2010,46(12)：128-134.

[11] 黃巍.大型水平軸風力機整機動力學研究[D].西安：西北工業大學，2012.

[12] 廖明夫，金路，劉展翅.某1.5MW風力機低電壓穿越測試報告[R].西安：旋轉機械與風能裝置測控研究所，2010.

[13] Guideline for the Certification of Wind Turbine [S].Hamburg:Germanischer Lioyd,2004.

[14] Gasch R,Twele J.Wind power plants—fundamental,design,construction and operation[M].Berlin:Solarpraxis AG,2002.

[15] 王勖成.有限單元法[M].北京：清華大學出版社，2003.

[16] Yokoyama T.Vibrations of a hanging Timoshenko beam under gravity[J].Journal of Sound and Vibration,1990,141(2):245-258.

[17] 李立本，宋憲耕，賀德馨，等.風力機結構動力學[M].北京：北京航空航天大學出版社，1999.

[18] 王新敏，李義強，許宏偉.ANSYS結構分析單元與應用[M].北京：人民交通出版社，2011.

[19] 金鑫,何玉林,劉樺,等.風力發電機耦合振動分析[J].振動與沖擊，2007，26(8)：144-148.

JIN Xin，HE Yu-lin，LIU Hua，et al.Coupled vibration analysis of wind turbine[J].Journal of Vibration and Shock, 2007，26(8)：144-148.

[20] 董禮.水平軸變速變槳風力機氣動設計與控制技術研究[D].西安：西北工業大學，2011.

[21] Hansen M O L.Aerodynamics of wind turbines[M].London:Earthscan,2008.

[22] Bir G S.User’s guide to bmodes[R].Golden:National Renewable Energy Laboratory,2007.

[23] NB/T 31051-2014. 風電機組低電壓穿越能力測試規程[S].北京：中國電科院，2014.

附錄：

剛度矩陣具體形式

其中：A為截面面積；E為彈性模量；L為單元長度；G為剪切模量；J為扭轉慣性矩。

ay=a(Iy,φx)

ax=a(Ix,φy)

by,bx,cy,cx,dy,dx,ey,ex,fy,fx形式同ay,ax

其中：

其中：Ix為x方向的截面慣性矩；Iy為y方向的截面慣性矩；kx為x方向的剪切系數；ky為y方向的剪切系數。

質量矩陣具體形式

其中：

Ay=A(ry,φx)

Ax=A(rx,φy)

By,Bx,Cy,Cx,Dy,Dx,Ey,Ex,Fy,Fx形式同Ay,Ax

其中：

F(r,φ)=

其中：Jz為極慣性矩；rx為回轉半徑；ry為回轉半徑。

Investigation of voltage sag impact on wind turbine blade vibrations

YIN Yao-jie, LIAO Ming-fu, Lü Pin, LIU Zhan-chi,WANG Si-ji

(School of Power and Energy, Northwestern Polytechnical University, Xi’an 710129, China)

Abstract:Recently the size and the capacity of wind turbines increase rapidly, the interactions between the wind turbines and the electric grid becomes more obvious. The disturbances from the grid and the faults of the grid will have impact on mechanical components of wind turbines. In order to obtain the dynamic responses of the blade under voltage sag, the blade was modeled as a cantilever and discretized by applying 3D Timoshenko beam element, considering the centrifugal stiffening of the rotating blade. The structural dynamic equation of the blade was established and the time-varying loads acting on the blade were calculated based on the data from a 1.5 MW wind turbine LVRT(low voltage ride through) test. The dynamic responses of the blade under the three-phase symmetrical voltage dip and three-phase asymmetry dip were calculated by applying Newmark method. The results show that, the voltage dip produces the transient loads and excites the vibrations of the blade. The amplitude of vibration in the flap wise is bigger than that in the edgewise. The vibration increases with the increase of the load of wind turbine and the magnitude of voltage drop. The dynamic responses of the blade under the three-phase symmetrical voltage dip are larger than those under three-phase asymmetry dip.

Key words:voltage sag; wind turbine; blade; Timoshenko beam ; dynamic response

基金項目：中央高?；究蒲袠I務專項資金資助(310201513J(Ⅱ)CG015)

收稿日期：2015-07-28修改稿收到日期：2015-10-13

通信作者廖明夫 男，博士,教授，博士生導師，1960年2月生

中圖分類號:TK83

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2016.10.031

第一作者 尹堯杰 男，博士生，1987年11月生