寧 偉，孫文磊

（新疆大學(xué)機械工程學(xué)院，新疆 烏魯木齊 830049）

1 引言

為了在大規(guī)模利用風(fēng)能的同時降低生產(chǎn)成本，大型化是風(fēng)力發(fā)電機組目前的發(fā)展趨勢。于此同時，風(fēng)力機高度以及體積的持續(xù)增大，使其制造難度與成本不斷升高，其承受的載荷也更加難以控制并對風(fēng)力發(fā)電機組的安全運行造成影響，甚至產(chǎn)生破壞。因此，風(fēng)力發(fā)電機固有頻率、動態(tài)響應(yīng)等研究對于風(fēng)力發(fā)電機組的設(shè)計與優(yōu)化具有重要意義。

目前，在風(fēng)力機結(jié)構(gòu)力學(xué)研究方面，文獻[1-2]將葉片與塔架看成柔性懸臂梁，使用二節(jié)點梁單元進行離散，分別考慮了氣動阻尼、離心剛化的影響。在此基礎(chǔ)上，文獻[3]采用kane方法建立風(fēng)力機動力學(xué)模型，通過假設(shè)模態(tài)法離散進行柔性化。文獻[4]考慮材料非線性以及截面偏心的影響，建立一種葉片非線性單元，并應(yīng)用于整機分析。針對上述方法建模復(fù)雜、重復(fù)性操作過多的問題，文獻[5-6]采用剛體有限元法研究結(jié)構(gòu)阻尼對振動特性的影響。文獻[7]采用超級單元方法將柔性葉片通過力元彈簧及阻尼器進行離散推導(dǎo)出葉片多體動力學(xué)方程組。

根據(jù)以上調(diào)研，在柔性多體動力學(xué)理論基礎(chǔ)上，選取新穎、簡單、計算速度快的剛體有限元方法（Rigid Finite Element Method，RFEM），建立了系統(tǒng)動力學(xué)模型，編寫了相關(guān)程序。計算結(jié)果與GH-Bladed數(shù)據(jù)結(jié)果對比，驗證了剛體有限元方法和模型建立的準確性，適用于大型風(fēng)力發(fā)電機的動力學(xué)分析計算。

2 數(shù)學(xué)模型

2.1 剛體有限元方法

將塔架和葉片簡化為質(zhì)量與剛度連續(xù)分布的柔性懸臂梁，采用剛體有限元二次劃分方法，通過剛體單元與彈性阻尼單元分別對塔架和葉片進行離散。剛體有限元法，如圖1所示。

圖1 剛體有限元法Fig.1 Rigid Finite Element Method

剛體有限元法中，剛體單元存儲位移與速度，彈簧阻尼節(jié)點則存儲變形能?？紤]每個剛體單元有6個自由度，則任意剛體單元的自由度表示為：

在局部坐標系｛｝E中任意剛體單元i上某點A的位置矢量r′與該點在全局坐標系｛｝0下的位置矢量r有以下關(guān)系：

由此推導(dǎo)出剛體單元的動能Ei：

將系統(tǒng)所有剛體單元與彈性阻尼節(jié)點相關(guān)量整合，得出Lagrange方程建立的動力學(xué)模型[8]：

式中：Q—系統(tǒng)所受的廣義力；假設(shè)剛體單元的3個旋轉(zhuǎn)自由度方向的轉(zhuǎn)角較小，則各阻尼節(jié)點的勢能與阻尼耗散能可推導(dǎo)為：

式中：ki=diag（），ci=diag（）；彈簧剛度系數(shù)及阻尼系數(shù)基于Kelvin-Voigt模型推導(dǎo)。

各剛體單元動能的Lagrange方程變換形式為：

式中：（lj）=1…6—對該剛體段變換矩陣偏導(dǎo)的各自由度方向，Mi=diag（mi，mi，mi）。

將式（10）～式（11）寫為矩陣的形式：

同理，由式（8）可得系統(tǒng)中各彈性阻尼節(jié)點的阻尼矩陣為：

最終，推導(dǎo)出式（6）線性化后的統(tǒng)一形式的結(jié)構(gòu)動力學(xué)方程組[9]：

2.2 固有動力特性

結(jié)構(gòu)作自由振動時，系統(tǒng)的激勵為0，式（14）可寫為：

令節(jié)點位移 q（t）=φsin（ωt+θ），得特征矩陣方程：

對系統(tǒng)固有頻率及振型的求解就歸結(jié)到式（16）的特征值與特征向量的求解，上式有非零解的條件為系數(shù)行列式為0，即：

得出其特征值與特征向量：

3 動態(tài)響應(yīng)計算

由于模態(tài)疊加法相比于其它動力學(xué)分析方法具有高效、計算速度快的優(yōu)點，動態(tài)響應(yīng)計算采用了模態(tài)疊加法，通過在模態(tài)計算中所獲得的固有頻率與振型來描述結(jié)構(gòu)的動態(tài)響應(yīng)。

3.1 模態(tài)疊加原理

將式（20）代入式（14），則系統(tǒng)的強迫振動方程：

由于振型具有正交性，假設(shè)阻尼矩陣也符合正交性。則可將式（22）離散化為非耦合的n個單自由度方程：

式中：Mi—第 i階模態(tài)質(zhì)量；Ci—第 i階模態(tài)阻尼；Ki—第 i階模態(tài)剛度，在式（23）兩邊同時除以模態(tài)質(zhì)量Mi：

經(jīng)坐標變換，物理坐標q轉(zhuǎn)換為模態(tài)坐標x：

式中：ξi—第 i階阻尼比；ωi=—第i階振動頻率。

在求解每個振型下的位移分量后帶入式（20）最終通過模態(tài)疊加法得到總位移：

3.2運動方程求解

由Newmark法可知各剛體段在每個時間步長下的位移、速度表達式：

式中：h—步長，α=0.25，β=0.5，將上式帶入式（23），得到加速度表達式：

求解式（26）～式（28），取時間步長 0.01s，得出位移、速度、加速度的動態(tài)響應(yīng)。

總結(jié)動態(tài)響應(yīng)的求解過程：

（1）計算固有頻率及振型。

（2）由振型的正交特性，分別計算各階模態(tài)質(zhì)量、模態(tài)阻尼與模態(tài)剛度。

（3）由每個時間步長下所受的載荷，計算得出各階模態(tài)載荷。

（4）離散方程為非耦合的n個單自由度方程。

（5）分別給定位移、速度和加速度的初始值，通過Newmark積分法求解各個單自由度方程。

（6）由式（25）得到動態(tài)響應(yīng)。

4 風(fēng)輪與塔架耦合

由于受風(fēng)輪載荷的影響，塔架將產(chǎn)生前后、左右方向的擺動，同時系統(tǒng)集合形狀也在不斷變化，若用標準的模態(tài)分析來分析風(fēng)力機的整體動態(tài)特性會變得越來越復(fù)雜。因此，需要將模態(tài)振型與各階頻率分開考慮，塔架模態(tài)中將風(fēng)輪視為剛性的，而風(fēng)輪模態(tài)中風(fēng)輪類似于懸掛在剛性主軸上，如圖2所示。

圖2 風(fēng)輪和塔架的基本形狀Fig.2 The Basic Shape of Wind Wheel and Tower

由于風(fēng)輪與塔架模態(tài)不成直角，因此葉片與塔架的軌跡方程彼此關(guān)聯(lián)，則由圖2可知在某時刻葉片J的偏轉(zhuǎn)可寫為[10]：

式中：φTJ—相對于剛性塔架的葉片模態(tài)振型；

φTJ—葉片J的位置偏差；

φTJ=1cosψ，ψ—葉片 J的方位角；

φT—塔架模態(tài)振型。

將式（29）帶入式（14），則葉片的運動方程為：

其中，塔架模態(tài)質(zhì)量與塔架、機艙、風(fēng)輪的質(zhì)量及風(fēng)輪的轉(zhuǎn)動慣量相關(guān)，Mi=MTφT+MN+MR+IR/L2。

式（30）和式（31）的動態(tài)分析過程為：

（1）帶入初始時刻的位移、速度、加速度和空氣動力載荷。

（2）假設(shè)方程右邊耦合項為不隨時間改變的常數(shù)，以便在解耦后的方程中消去同類項。

（3）解方程算出位移與速度的增量，從而得出葉片J的葉尖位移與速度增量。

（4）解運式（30）、式（31），得出最后一個時間步長下的加速度。

（5）改變方程右邊的耦合項，重新求解運動增量方程，得到位移和速度增量的修正值。

5 算例分析

以美國國家新能源實驗室的某NREL1.5WM風(fēng)力發(fā)電機[11]為研究對象，對其進行模態(tài)與動力學(xué)分析，主要參數(shù)，如表1所示。

表1 風(fēng)力機結(jié)構(gòu)參數(shù)Tab.1 Structural Parameters of Wind Turbine

5.1 模態(tài)分析

由于在風(fēng)力發(fā)電機模態(tài)分析中，低頻振動要比高頻振動危險，高階模態(tài)對響應(yīng)的影響較小，同時結(jié)構(gòu)阻尼的作用使響應(yīng)中高階部分衰減很快，因此忽略了高階模態(tài)的影響[12]。首先計算葉片平面外與平面內(nèi)前三階固有頻率，葉片固有頻率和振型的計算結(jié)果同GH-Bladed仿真結(jié)果對比，如表2、圖3所示。

表2 葉片模態(tài)Tab.2 Natural Frequencies of Blade

圖3 葉片固有振型Fig.3 Mode Shapesof Blade

計算塔架前后與左右方向的前兩階彎曲固有頻率，其彎曲固有頻率和振型結(jié)果的對比，如表3、圖4所示。

表3 塔架模態(tài)分析Tab.3 Natural Frequencies of Tower

圖4 塔架固有振型Fig.4 Mode Shapesof Tower

由于塔架的對稱性，一階前后向與左右向、二階前后向與左右向的振型相同。從塔架振型圖可以看出：一階振型塔架頂端的振幅最大，二階振型塔架頂端的振幅較小，中間的振幅較大。

通過模態(tài)分析對比，初步驗證了模型建立的正確性，同時為動態(tài)響應(yīng)的計算提供了計算參數(shù)。

5.2 動態(tài)響應(yīng)分析

風(fēng)力發(fā)電機組運行時，對于風(fēng)輪所產(chǎn)生的時變動態(tài)載荷，通過GH-Bladed來獲取，輪轂高度處的縱向風(fēng)速，如圖5所示。湍流風(fēng)采用von-karman模型，平均風(fēng)速為11.8m/s，仿真時間70s，步長0.05s。

圖5 縱向風(fēng)速Fig.5 Longitudinal Wind Speed

由于風(fēng)力機在來流方向上的變形遠遠大于其左右方向變形，因此先考慮了葉尖平面外和塔架前后向的變形情況，計算對葉片的前三階模態(tài)與塔架的前兩階模態(tài)進行了疊加。

葉片與塔架的動態(tài)響應(yīng)計算結(jié)果分別，如圖6、圖7所示。從分析結(jié)果可以看出實際工作狀況下，葉尖在風(fēng)輪平面外的偏移量較大，最大值為2.59 m，由于數(shù)值計算時初始位移、速度、加速度的值均從零開始，因此在起始時刻葉片與塔架的位移與速度均出現(xiàn)劇烈變化7s后趨于穩(wěn)定，計算同GH-Bladed仿真結(jié)果保持一致。

圖6 葉片尖端動態(tài)響應(yīng)Fig.6 The Dynamic Analysis of Bladed Tip

圖7 塔架頂部動態(tài)響應(yīng)Fig.7 The Dynamic Analysis of Tower Top

葉片與塔架在來流方向較大的變形量不僅影響氣動載荷的變化，更對風(fēng)力發(fā)電機組的運行安全造成較大影響，甚至可能會產(chǎn)生破壞。因此，設(shè)計過程中對其來流方向的動態(tài)響應(yīng)需要充分考慮。

仿真結(jié)果表明：剛體有限元法在兼顧較少的建模計算時間的情況下，能夠較為準確地進行風(fēng)力發(fā)電機動力學(xué)分析。

6 結(jié)論

（1）針對風(fēng)力發(fā)電機組動力學(xué)問題，利用剛體有限元法進行了離散化建模，建立了葉片與塔架耦合的結(jié)構(gòu)動力學(xué)模型，給出了剛體有限元法的模態(tài)分析與動態(tài)響應(yīng)求解方法與步驟。

（2）利用模型對1.5MW風(fēng)力機的模態(tài)與動態(tài)響應(yīng)進行計算。計算結(jié)果同GH-Bladed進行了對比驗證，保持了良好的一致性，證明建立的模型簡化合理、計算準確，能較好預(yù)測風(fēng)力機組對于時變載荷的動態(tài)響應(yīng)，具有良好的實用價值，適用于大型風(fēng)力機的動力學(xué)分析計算。