曹 莉， 孫文磊， 周建星

(新疆大學(xué) 機(jī)械工程學(xué)院，烏魯木齊 830047)

由于全球性能源短缺和環(huán)境污染等問題，風(fēng)能作為一種清潔的可再生能源，變得越來越重要[1]。風(fēng)力發(fā)電機(jī)組常年工作在變速、變載荷的工作環(huán)境中，在隨機(jī)風(fēng)的作用下會(huì)受到頻繁的擾動(dòng)和動(dòng)載激勵(lì)，對(duì)其工作性能和使用壽命有很大的影響[2]。隨著風(fēng)電機(jī)組不斷向大型化、大功率方向發(fā)展，風(fēng)力發(fā)電機(jī)組的動(dòng)態(tài)特性更加復(fù)雜。因此，在風(fēng)力機(jī)的設(shè)計(jì)過程中，需要考慮風(fēng)力機(jī)在典型工況下的動(dòng)態(tài)特性。

當(dāng)前，在風(fēng)力機(jī)動(dòng)力學(xué)分析領(lǐng)域，常用的方法有多體動(dòng)力學(xué)方法MBS (Multibody Systems)、有限元方法FES(Finite Element System)、模態(tài)分析方法及連續(xù)系統(tǒng)COS(Continuous Systems)等[3]，MBS多體動(dòng)力學(xué)方法將實(shí)際的機(jī)械構(gòu)件視為剛體來建立系統(tǒng)動(dòng)力學(xué)微分方程(組)，該方法模擬的精度有限[4]；COS連續(xù)系統(tǒng)方法建立的偏微分方程組僅在特殊、簡單的幾何結(jié)構(gòu)及載荷下才可以求解。在風(fēng)力機(jī)整機(jī)動(dòng)力學(xué)方面主要采用有限元方法和模態(tài)分析方法[5]。但有限元方法具有較多的自由度，計(jì)算和分析成本較高[6]；模態(tài)疊加法通過提取結(jié)構(gòu)的主模態(tài)集(主振型陣)形成坐標(biāo)變換并實(shí)現(xiàn)降價(jià)來縮短計(jì)算時(shí)間，相對(duì)于其它方法可使用較少的自由度對(duì)風(fēng)力機(jī)動(dòng)態(tài)行為進(jìn)行可靠分析，具有計(jì)算速度較快、效率高等優(yōu)點(diǎn)。

近年來，國內(nèi)外學(xué)者對(duì)風(fēng)力發(fā)電機(jī)組的動(dòng)態(tài)特性問題進(jìn)行了研究： Asareh等[7]以FAST軟件為平臺(tái)，結(jié)合氣動(dòng)力分析了5 MW風(fēng)力機(jī)在不同工況下的動(dòng)態(tài)響應(yīng)；劉雄等[8]研究了風(fēng)力機(jī)在湍流風(fēng)作用下的動(dòng)態(tài)特性，并揭示了葉片的離心剛化和氣動(dòng)阻尼對(duì)分析結(jié)果的影響較大；呂計(jì)男等[9]對(duì)葉片和塔架進(jìn)行簡化并建立了風(fēng)力機(jī)整機(jī)模型，計(jì)算了風(fēng)力機(jī)在啟動(dòng)工況下的動(dòng)態(tài)響應(yīng)及葉片幾何非線性對(duì)結(jié)果的影響。李明等[10]采用車載法模擬自然風(fēng)場(chǎng)，對(duì)風(fēng)力機(jī)葉片的相似模型進(jìn)行測(cè)試，得到的葉片各截面的載荷數(shù)據(jù)，與BLADED仿真結(jié)果較為一致。研究中風(fēng)力機(jī)所處工況多為單一工況，關(guān)于風(fēng)力機(jī)在不同工況下的動(dòng)態(tài)特性研究較少；研究對(duì)象多為葉片、塔架、輪轂等單個(gè)部件，鮮有文獻(xiàn)對(duì)風(fēng)力機(jī)整機(jī)的動(dòng)態(tài)特性研究。

因此，筆者針對(duì)復(fù)雜工況下的大型風(fēng)力機(jī)，采用了基于模態(tài)疊加法的動(dòng)態(tài)特性分析方法，建立了風(fēng)力機(jī)動(dòng)態(tài)特性計(jì)算模型，采用該方法對(duì)2 MW大型風(fēng)力機(jī)進(jìn)行分析，同時(shí)總結(jié)了不同工況下的風(fēng)力機(jī)的動(dòng)態(tài)特性。

1 風(fēng)力機(jī)動(dòng)態(tài)特性計(jì)算理論

1.1 風(fēng)力機(jī)葉片載荷計(jì)算理論

風(fēng)力機(jī)在運(yùn)行時(shí)主要通過葉片捕獲風(fēng)能，捕獲的風(fēng)能帶動(dòng)風(fēng)輪轉(zhuǎn)動(dòng)，再通過傳動(dòng)系統(tǒng)將風(fēng)能轉(zhuǎn)化為機(jī)械能，因此葉片是主要的受力部件[11]。葉片所受的力主要包括空氣動(dòng)力、離心力、和重力。葉片的受力和變形坐標(biāo)系如圖1所示，其中 為迎風(fēng)向， 為葉片展長方向， 由右手規(guī)則確定。

圖1 葉片受力和變形坐標(biāo)系Fig.1 Coordinate system for blade loads and deflections

對(duì)葉片的載荷進(jìn)行綜合可以計(jì)算出葉片單位長度所受的載荷[12]：

(1)

式中：FXg，F(xiàn)Ya為葉片單位長度上的氣動(dòng)載荷；FXg，F(xiàn)Yg，F(xiàn)Zg為葉片單位長度上的重力；FYc，F(xiàn)Zc為葉片單位長度上的離心力。由于筆者所考慮的風(fēng)力機(jī)風(fēng)輪是靜止的，氣動(dòng)載荷和離心力可忽略不計(jì)，因此迎風(fēng)向載荷FXa是主要載荷。

1.2 基于模態(tài)疊加法的風(fēng)力機(jī)動(dòng)態(tài)特性計(jì)算理論模型

筆者將有限元法和模態(tài)疊加法這兩種方法進(jìn)行了結(jié)合，其基本思想是以有限元為平臺(tái)，反復(fù)利用瞬態(tài)分析求解風(fēng)力機(jī)在周期性載荷作用下的穩(wěn)態(tài)響應(yīng)。瞬態(tài)分析包括完全法、縮減法和模態(tài)疊加法。完全法計(jì)算量大且計(jì)算時(shí)間較長；縮減法不能在時(shí)間單元上添加載荷且所有荷載必須加載用戶定義的主自由度上；模態(tài)疊加法相對(duì)于其它方法計(jì)算速度較快、效率高。因此本文選用模態(tài)疊加法對(duì)風(fēng)力機(jī)整機(jī)動(dòng)態(tài)特性進(jìn)行求解，得到風(fēng)力機(jī)在時(shí)變載荷作用下的多自由度運(yùn)動(dòng)程[13]：

(2)

[ψ]為系統(tǒng)的振型矩陣，作如下坐標(biāo)變換：

u=[ψ]η

(3)

式中：η=[u1u2…un]T。

則系統(tǒng)的強(qiáng)迫振動(dòng)方程為：

(4)

式(4)乘以特征振型[ψ]T可得：

(5)

由于[ψ]具有正交性，則主坐標(biāo)下的方程為：

(6)

由主振型的正交性可得：

(7)

將式(7)代入式(6)，得到二階微分方程：

(8)

求出每個(gè)振型坐標(biāo)上的位移分量后，采用模態(tài)疊加法可以得到各個(gè)節(jié)點(diǎn)自由度上的位移響應(yīng)：

(9)

2 風(fēng)力機(jī)動(dòng)態(tài)特性有限元計(jì)算方法

2.1 風(fēng)力機(jī)有限元模型

選取陸地某2 MW大型風(fēng)力機(jī)為研究對(duì)象，葉輪直徑為40 m，輪轂直徑為2.5 m，塔架高75 m。葉片材料采用玻璃纖維，彈性模量為1.76×1010Pa，泊松比為0.17；塔架為Q345E合金鋼，彈性模量為2.06×1011Pa，泊松比為0.28。葉片結(jié)構(gòu)采用梁單元，塔頂質(zhì)量采用點(diǎn)質(zhì)量單元?jiǎng)澐志W(wǎng)格。塔架底部截面采用全約束，約束該截面節(jié)點(diǎn)的所有自由度。風(fēng)力機(jī)有限元模型如圖2所示，單元節(jié)點(diǎn)數(shù)為93 389，單元數(shù)為462 396。

圖2 風(fēng)力機(jī)有限元模型Fig.2 Finite element model of the wind turbine

2.2 風(fēng)力機(jī)動(dòng)態(tài)響應(yīng)有限元計(jì)算方法

針對(duì)上述建立的基于模態(tài)疊加法的風(fēng)力機(jī)動(dòng)態(tài)特性計(jì)算模型，采用有限元法對(duì)其動(dòng)態(tài)特性進(jìn)行計(jì)算，進(jìn)行求解的具體計(jì)算流程，如圖3所示。

圖3 風(fēng)力機(jī)整機(jī)動(dòng)態(tài)特性有限元法計(jì)算流程Fig.3 Dynamic characteristics computing process forwind turbines

圖3中的外部程序1是采用葉素動(dòng)量定理計(jì)算葉片上的氣動(dòng)力[14]，主要是葉片氣動(dòng)性能計(jì)算和速度誘導(dǎo)因子迭代，Wilson等[15]詳細(xì)論述了該計(jì)算過程，本文不再贅述；主程序主要是通過二次開發(fā)語言APDL將葉片迎風(fēng)向載荷讀入ANSYS，并將其離散成沖擊載荷，在運(yùn)用模態(tài)疊加法求解風(fēng)力機(jī)動(dòng)態(tài)響應(yīng)；外部程序2主要判斷風(fēng)力機(jī)在求解動(dòng)響應(yīng)時(shí)是否達(dá)到穩(wěn)定狀態(tài)，周期求解動(dòng)響應(yīng)所允許的誤差為ε，若滿足式(10)，則表明動(dòng)響應(yīng)求解已進(jìn)入穩(wěn)定狀態(tài)，則可停止計(jì)算；通過主程序及外部程序2求得風(fēng)力機(jī)動(dòng)態(tài)響應(yīng)后，根據(jù)外部程序3進(jìn)行擴(kuò)展可進(jìn)一步得到風(fēng)力機(jī)的Von Mises動(dòng)應(yīng)力。

(10)

2.3 風(fēng)力機(jī)模態(tài)分析

模態(tài)分析是動(dòng)響應(yīng)分析的基礎(chǔ)，采用Block Lanczons對(duì)風(fēng)力機(jī)進(jìn)行模態(tài)分析并提取前10階固有頻率，如表1所示。由于篇幅有限。僅列出風(fēng)力機(jī)前五階模態(tài)振型，如圖4。

表1 風(fēng)力機(jī)模態(tài)計(jì)算結(jié)果

圖4 風(fēng)力機(jī)模態(tài)振型(1～5階)Fig.4 Modal vibration mode (1st to 5st order) of thewind turbine

由圖5可知：一階振型主要時(shí)風(fēng)力機(jī)上面兩葉片沿X軸前后彎曲振動(dòng)；二階振型主要風(fēng)力機(jī)整機(jī)沿X軸前后彎曲振動(dòng)；三階、四階、五階振型主要是葉片一階揮舞振動(dòng)。對(duì)于三階振型，其下葉片保持不動(dòng)，其余兩葉片沿X軸異向振動(dòng)；四階振型是下葉片沿X軸反向振動(dòng)；五階振型為三葉片沿Y軸同向振動(dòng)。

3 典型工況下風(fēng)力機(jī)動(dòng)態(tài)特性分析

3.1 穩(wěn)態(tài)風(fēng)況下風(fēng)力機(jī)動(dòng)態(tài)特性分析

計(jì)算可得風(fēng)力機(jī)輪轂在恒定風(fēng)速(即額定風(fēng)速10 m/s)下的動(dòng)載激勵(lì)，如圖5所示。根據(jù)風(fēng)力機(jī)動(dòng)態(tài)特性計(jì)算流程可得到風(fēng)力機(jī)輪轂的位移動(dòng)響應(yīng)，如圖6。由圖5、6可知，輪轂載荷的在266.25～274.74 kN之間呈周期性變化，載荷頻率由1/3次諧波、0.9 Hz主頻成分及二倍頻組成；輪轂位移動(dòng)響應(yīng)在0.88～0.80 m之間波動(dòng)，動(dòng)響應(yīng)主頻率成分為0.30 Hz，與輪轂1/3次諧波頻率(葉片主頻率)一致。

圖7為風(fēng)力機(jī)整機(jī)的Von Mises應(yīng)力云圖，可以看出塔架底部應(yīng)力最大，平均應(yīng)力達(dá)152×106Pa，且塔架底部至塔架頂部應(yīng)力逐漸減小；下葉片的葉跟至葉尖的動(dòng)載激勵(lì)和Von Mises應(yīng)力如圖8所示，可以看出葉根應(yīng)力最大，最大值為2.19×107Pa；從葉根至葉片11 m處時(shí)Von Mises應(yīng)力呈拋物線型減小，葉片11 m處至葉尖，Von Mise應(yīng)力呈減小趨勢(shì)但變化不大。

圖5 10 m/s恒定風(fēng)速下風(fēng)力機(jī)輪轂動(dòng)載激勵(lì)Fig.5 Dynamic loads of hub under a constant wind speedof 10 m/s

圖6 10m/s恒定風(fēng)速下風(fēng)力機(jī)輪轂位移動(dòng)響應(yīng)Fig.6 Displacement dynamic response of hub under a constantwind speed of 10m/s

圖7 10 m/s恒定風(fēng)速下風(fēng)力機(jī)應(yīng)力云圖Fig.7 The stress nephogram of the wind turbine under a constantwind speed of 10 m/s

圖8 葉片根部至葉尖動(dòng)載激勵(lì)和Von Mises應(yīng)力Fig.8 The loads and stress of blade from root to tip

3.2 不同風(fēng)況下風(fēng)力機(jī)動(dòng)態(tài)特性分析

根據(jù)風(fēng)力機(jī)整機(jī)動(dòng)態(tài)特性計(jì)算方法計(jì)算了風(fēng)力機(jī)在不同風(fēng)況下的動(dòng)態(tài)特性，總結(jié)不同工況下風(fēng)力機(jī)位移動(dòng)響應(yīng)、Von Mises應(yīng)力最大值及其相應(yīng)位置，如表2所示。

表2 不同工況下風(fēng)力機(jī)動(dòng)態(tài)特性

由表2可知，風(fēng)力機(jī)在不同風(fēng)況下的最大位移動(dòng)響應(yīng)位置在都葉尖處；額定風(fēng)速下風(fēng)力機(jī)的位移響應(yīng)最大值在所有工況中最小，為2.59 m；正常停機(jī)工況下的位移響應(yīng)最大值為2.89 m，由于停機(jī)時(shí)產(chǎn)生了強(qiáng)烈的振動(dòng)使得葉尖位移響應(yīng)值瞬間減小到了-1.86 m；極端運(yùn)行陣風(fēng)作用下的位移動(dòng)響應(yīng)最大值在所有工況中最大，為7.19 m，較額定風(fēng)速下的位移動(dòng)響應(yīng)最大值增加可177.22%。風(fēng)力機(jī)在啟動(dòng)工況下的Von Mises應(yīng)力最大在在所有工況中最小，為0.44×108Pa；極端運(yùn)行陣風(fēng)工況下的Von Mises應(yīng)力最大值在所有工況中最大，為4.05×108Pa，較額定風(fēng)速下的Von Mises應(yīng)力最大值增加了820.45%；且極端運(yùn)行工況下的Von Mises應(yīng)力最大值的位置在塔架頂部，其余工況都在塔架底部。在所有工況中，極端運(yùn)行陣風(fēng)的位移響應(yīng)最大值和Von Mises應(yīng)力最大值最大，在風(fēng)力機(jī)設(shè)計(jì)中應(yīng)著重考慮。

3.3 極端運(yùn)行陣風(fēng)下風(fēng)力機(jī)動(dòng)態(tài)特性分析

風(fēng)力機(jī)運(yùn)行環(huán)境十分復(fù)雜，由于丘陵、建筑物和森林等障礙物的阻擋空氣的流動(dòng)會(huì)造成許多不規(guī)則的渦旋，渦旋的流動(dòng)方向與空氣流動(dòng)方向一致時(shí)，會(huì)產(chǎn)生瞬時(shí)極大風(fēng)速，相反，會(huì)產(chǎn)生瞬時(shí)極小風(fēng)速。極端運(yùn)行陣風(fēng)是指風(fēng)速突然下降，接著陡然上升，然后再突然下降，最后又上升到初始值的過程。通過風(fēng)力機(jī)動(dòng)態(tài)特性計(jì)算方法可以求得風(fēng)力機(jī)輪轂、葉尖的位移響應(yīng)(圖9)，輪轂、塔架頂部的Von Mises應(yīng)力(圖10)。

由圖9可知在穩(wěn)態(tài)風(fēng)(0～15 s)作用下，輪轂位移響應(yīng)在0.55～0.80 m之間波動(dòng)，葉尖位移響應(yīng)在0.76～2.92 m之間波動(dòng)。在陣風(fēng)(15～40 s)作用下，風(fēng)力機(jī)發(fā)生了劇烈的振動(dòng)，輪轂振幅最大值達(dá)2.08 m，較穩(wěn)態(tài)振幅均值增加了208.15%，最小值達(dá)在-1.18 m，較穩(wěn)態(tài)振幅均值反向增加了274.18%；葉尖振動(dòng)幅值最大為7.19 m，較葉尖穩(wěn)態(tài)振動(dòng)均值增加了290.76%，最小為-4.78 m，較穩(wěn)態(tài)均值反向增加了359.78%。陣風(fēng)消失后，風(fēng)力機(jī)逐漸回歸到穩(wěn)定狀態(tài)。

圖9 陣風(fēng)作用下輪轂、葉尖位移動(dòng)響應(yīng)Fig.9 Displacement dynamic responses of hub and blade tipunder transient forcing-flurry

由圖10可知在穩(wěn)態(tài)風(fēng)(0～15 s)作用下，風(fēng)力機(jī)輪轂的Von Mises應(yīng)力范圍為0.30×106～2.19×106Pa，塔架頂部的Von Mises應(yīng)力范圍為0.08×108～1.69×108Pa；在陣風(fēng)(15～40 s)作用下，輪轂的Von Mises應(yīng)力最大值為5.3×106Pa，相對(duì)于其穩(wěn)態(tài)應(yīng)力均值增加了325.70%；塔架頂部的Von Mises應(yīng)力最大值為4.05×108Pa，相對(duì)于其穩(wěn)態(tài)應(yīng)力均值增加了357.63%。陣風(fēng)結(jié)束后，輪轂及塔架頂部的Von Mises應(yīng)力減小至最小值，然后逐漸增加，最終呈周期性波動(dòng)。

圖10 陣風(fēng)作用下風(fēng)力機(jī)輪轂塔架頂部的Von Mises應(yīng)力Fig.10 The stress of hub and the top of tower under transientforcing-flurry

4 結(jié) 論

針對(duì)典型工況下風(fēng)力機(jī)整機(jī)動(dòng)態(tài)特性問題，提出了基于模態(tài)疊加法的風(fēng)力機(jī)動(dòng)態(tài)響應(yīng)計(jì)算模型，進(jìn)而采用有限元方法對(duì)其動(dòng)態(tài)響應(yīng)進(jìn)行計(jì)算。對(duì)2 MW風(fēng)力機(jī)在不同工況下進(jìn)行動(dòng)態(tài)特性分析，得到以下結(jié)論：

(1)穩(wěn)態(tài)風(fēng)作用下風(fēng)力機(jī)塔架底部應(yīng)力最大，平均應(yīng)力達(dá)152×106Pa，且塔架底部至塔架頂部應(yīng)力逐漸減??；葉片的葉跟處應(yīng)力最大，最大值達(dá)2.19×107Pa；其Von Mises應(yīng)力從葉根至尖處時(shí)呈拋物線型減小。

(2)計(jì)算并總結(jié)了不同工況下的風(fēng)力機(jī)動(dòng)態(tài)特性，發(fā)現(xiàn)風(fēng)力機(jī)的動(dòng)響應(yīng)最大的位置都在葉尖處，極端運(yùn)行陣風(fēng)下的Von Mises應(yīng)力最大位置為塔架頂部，其余工況都在塔架底部；且極端運(yùn)行陣風(fēng)作用下的最大位移響應(yīng)值、最大Von Mises應(yīng)力值在所有工況中均為最大，在風(fēng)力機(jī)設(shè)計(jì)時(shí)應(yīng)著重考慮。

(3)極端運(yùn)行陣風(fēng)條件下，葉尖的位移動(dòng)響應(yīng)最大值較其穩(wěn)態(tài)響應(yīng)均值增加了290.76%，最小值較其穩(wěn)態(tài)響應(yīng)均值反向增加了359.78%；塔架頂部的Von Mises應(yīng)力最大值較其穩(wěn)態(tài)應(yīng)力均值增加了357.63%。