王 磊，陳 柳，何玉林，劉 樺，金 鑫

（1.電子科技大學(xué) 能源科學(xué)與工程學(xué)院，成都 611731；2.重慶大學(xué) 機械學(xué)院，重慶 400044）

水平軸風力發(fā)電機系統(tǒng)作為具有強非線性流固耦合特性的周期時變動力學(xué)系統(tǒng)，其結(jié)構(gòu)特征及其運動較為復(fù)雜。風在風力發(fā)電機組部件上產(chǎn)生的氣動載荷與部件的變形之間相互作用，使風力發(fā)電機組系統(tǒng)具有氣彈耦合現(xiàn)象，系統(tǒng)柔性增加使各部件之間的運動和受力相互耦合更加顯著，解決這種柔性較大的系統(tǒng)的動力學(xué)問題，需要采用柔體動力學(xué)方法。近幾年來，國內(nèi)外許多學(xué)者對風力發(fā)電機組的動力學(xué)進行了研究：文獻［1］利用有限元法研究了槳葉、耦合轉(zhuǎn)子/機艙、塔筒系統(tǒng)的動力響應(yīng)。文獻［2］研究了風機氣動性能和結(jié)構(gòu)動力學(xué)特性，建立了槳葉半剛性模型。文獻［3］根據(jù)Hamilton原理推導(dǎo)了10自由度梁單元離散槳葉的運動方程。文獻［4］使用假設(shè)模態(tài)方法研究旋轉(zhuǎn)運動柔性梁動力特性，并通過與有限元離散化方法比較發(fā)現(xiàn)高速運動的系統(tǒng)，需要增加模態(tài)數(shù)目來降低誤差。

采用Kane方法［4-6］建立的風力機多剛體系統(tǒng)動力學(xué)模型，依據(jù)假設(shè)模態(tài)法進行離散化，并結(jié)合修正的BEM模型進行軸向誘導(dǎo)因子a過渡點修正，進行空氣動力載荷計算氣動力模型，建立風力發(fā)電機組整機系統(tǒng)動力學(xué)模型。以某2.5 MW風力發(fā)電機組為例，進行數(shù)值仿真計算，并對風力發(fā)電機組在隨機風場下的動力學(xué)運動規(guī)律以及動力學(xué)特性進行分析。

1 風力機多剛體系統(tǒng)動力學(xué)模型

文中采用Kane方法建立風力機多剛體系統(tǒng)動力學(xué)模型。首先確定以塔筒底部中心為原點中心，垂直于風輪平面的來流方向X方向，塔筒未變形時的中線方向為Z向的慣性坐標系E，如圖1所示。

圖1 三葉片風力機系統(tǒng)模型Fig.1 The model of three-bladed wind turbine

將風力發(fā)電機組視為具有N個自由度的一階線性多剛體系統(tǒng)，則風力機系統(tǒng)的運動可以由N個廣義坐標qi（i=1，2，…，N）描述，或可由N個廣義速率ur（r=1，2，…，N）描述。后者是從組成系統(tǒng)的質(zhì)點速度或剛體角速度的中任意選擇的N個獨立標量，表示為廣義坐標的導(dǎo)數(shù)的線性組合。

其中Yri、Zr為廣義坐標qi和時間t的函數(shù)，若ur為獨立變量，可以從方程（1）唯一的解出，得到：

當確定了每個剛體的偏速度和偏角速度，以及相應(yīng)的廣義主動力Fr和廣義慣性力之后，其Kane動力學(xué)方程為：

即每個廣義速率對應(yīng)的廣義主動力和廣義慣性力之和等于零。設(shè)風力發(fā)電機組系統(tǒng)由w個剛體組成。假設(shè)對于每一個剛體Ni，主動力施加在其質(zhì)心Xi和分別為作用在作用在每個剛體Ni上的主動力與力矩，則風力發(fā)電機組系統(tǒng)的廣義主動力為：

對應(yīng)的廣義慣性力為：

其中：

將式（6）、式（7）帶入式（5），風力發(fā)電機組系統(tǒng)動力學(xué)方程矩陣形式則為：

其中的［C（q，t）］為系統(tǒng)加速度的系數(shù)矩陣，{f（，q，t）}為系統(tǒng)位移和速度相關(guān)的向量。求解時，在每個時間步，方程的數(shù)值解的第一步是采用四階Adams-Beshforth預(yù)測-修正算法的預(yù)測方法確定低階項的值，并以此以構(gòu)成方程的右邊項，然后采用Gauss消元法求解系統(tǒng)自由度的加速度，這些計算得到加速度值用于修正預(yù)測值，以提高預(yù)測精度。經(jīng)過幾次迭代后，采用四階Adams-Mounton預(yù)測-修正算法的修正方法確定加速度的值，并給出該時間步的最終解。由于該預(yù)測-修正算法不是自發(fā)的，前四個時間步的解需要用四階Runge-Kutta法確定。

2 基于假設(shè)模態(tài)法的風力機動力學(xué)模型離散化

將風力機葉片和塔筒視為質(zhì)量和剛度連續(xù)分布的柔性懸臂梁。根據(jù)連續(xù)力學(xué)理論，該結(jié)構(gòu)體具有無限多的自由度，因此需要無限多個坐標確定葉片和塔筒的每一個點的位置。但在假設(shè)模態(tài)法中，葉片和塔筒的連續(xù)體的變形被表示為一系列規(guī)格化的振動模態(tài)的振型線性疊加［7-8］。這樣可以將葉片和塔筒的自由度從無限大減小為N（N表示計算時，選取的假設(shè)模態(tài)數(shù)目）。因此，柔性懸臂梁（葉片和塔筒）在任意時間、任意位置的橫向變形u（z，t）是一系列規(guī)格化模態(tài)振型φa（z）的線性疊加，它們與廣義坐標qa（t）的關(guān)系為：

其中：φa（z）是為懸臂梁的第a個模態(tài)振型，它是懸臂梁縱向方向上的距離z的函數(shù)；而與模態(tài)a相關(guān)的廣義坐標qa（t）是時間的函數(shù)，它通常被取為該模態(tài)振型的懸臂梁自由端的變形。

采用保守的、不隨時間變化的系統(tǒng)的Lagrange方程，一個N個自由度系統(tǒng)的運動方程為：

其中的cb（t）是關(guān)于形函數(shù)φb（z）的廣義坐標；廣義質(zhì)量mij和廣義剛度kij分別由動能T和勢能V定義為：

當懸臂梁以特定的固有頻率ωa振動時，假設(shè)a=m，有：

其中的Qa是懸臂梁的自由端部在自由模態(tài)a時的變形大小。將qa（t）代入方程cb（t）=Cm，bqm（t），并將得到的結(jié)果代入方程（2），得到：

改寫成矩陣形式為：

其中，廣義質(zhì)量矩陣［M］和廣義剛度矩陣［K］都是N×N階矩陣，系數(shù)向量{C}為N×1階向量。對該矩陣方程進行特征值分析，得到特征值和特征向量{C}a。

2.1 塔筒的廣義質(zhì)量與廣義剛度

塔筒被模擬為一個倒立的懸臂梁，它的自由端固定著一個點質(zhì)量Mtop，表示基座盤、機艙、輪轂和葉片的總質(zhì)量。塔筒的廣義質(zhì)量表示為：

其中μT（h）是塔筒的質(zhì)量分布線密度。

塔筒的勢能由與梁的分布剛度有關(guān)的分量VBeam和與重力有關(guān)的分量VGravity組成：

其中與梁的分布剛度相關(guān)的勢能分量為：

其中的EIT（h）是塔筒的分布剛度，H是塔筒的高度。

與重力相關(guān)的勢能分量為：

其中負號表示重力將減小塔筒的廣義剛度，括號內(nèi)的第一項與塔頂質(zhì)量的重力有關(guān)，第二項與塔筒分布質(zhì)量的重力有關(guān)，并考慮了塔筒變形對重力勢能的影響。

與方程（4）比較并整理后，得到塔筒的廣義剛度為：

2.2 葉片的廣義質(zhì)量與剛度

風輪的每個葉片被模擬成自由端處固定一個點質(zhì)量MTip的旋轉(zhuǎn)懸臂梁，以角速度Ω繞垂直于風輪平面的軸旋轉(zhuǎn)。假設(shè)每個葉片的柔性部分獨立地向葉片的拍動方向（垂直于翼型弦線）和揮舞方向（平行于翼型弦線）運動。同時，葉片的變形也可以表示為面內(nèi)方向（平行于風輪旋轉(zhuǎn)平面）和面外方面（垂直于風輪旋轉(zhuǎn)平面）的變形運動，拍動-揮舞與面內(nèi)-面外兩個坐標系之間的關(guān)系如圖2所示。

圖2 拍動-揮舞坐標系和面內(nèi)-面外坐標系Fig.2 Relationship between the flapwise-edgewise frame and inplane-out-of-plane frame

在隨風輪旋轉(zhuǎn)的坐標系中，葉片的動能與塔筒的動能的形式是一致的，葉片的廣義質(zhì)量可以寫為如下的形式：

其中，μB（r）葉片的分布線密度，R是風輪半徑，RH是輪轂半徑。

忽略重力，葉片的勢能由與葉片分布剛度有關(guān)的分量VBeam和與葉片旋轉(zhuǎn)產(chǎn)生的離心力作用相關(guān)的分量VRot組成：

與梁的分布剛度相關(guān)的葉片勢能分量為：

其中的EIB（h）是葉片的拍打方向或揮舞方向的分布剛度。

由風輪旋轉(zhuǎn)引起的勢能分量為：

括號內(nèi)的第一項是與葉尖剎車部分質(zhì)量相關(guān)的勢能分量，第二項是與葉片分布質(zhì)量相關(guān)的勢能分量。兩者的離心勢能之和取正值，表示離心力將增加葉片的廣義剛度，即離心剛化現(xiàn)象。

經(jīng)過整理簡化的葉片的廣義剛度為

2.3 基于假設(shè)模態(tài)表示的葉片和塔筒變形

根據(jù)葉片和塔筒的模態(tài)振型，可以表示葉片和塔筒的柔性變形位移。

3 修正的BEM空氣動力學(xué)模型

由于建立BEM理論時所做的假設(shè)，經(jīng)典BEM理論僅適用于軸對稱穩(wěn)態(tài)氣流中的風輪氣動力計算，且沒有考慮三維效應(yīng)的影響，適用范圍有限。因此需要對經(jīng)典BEM理論有重載風輪的CT進行修正［9］，以滿足實際計算的需要。

當風輪處于高速重載的情況下，葉片上的軸向誘導(dǎo)因子a較大，若a＞0.5時，將出現(xiàn)尾流風速UW為負數(shù)的情況，即尾流出現(xiàn)了倒流，這顯然是不符合實際情況下。實際上，此時的尾流會從周圍空氣中獲得能量，使緩慢通過風輪的氣流重新獲得能量，風輪上下風向的壓力差將產(chǎn)生很大的推力，遠遠大于動量理論的預(yù)測值，因此需要對CT進行修正以反映實際的情況。修正CT值的方法有很多，它們之間的差別在于軸向誘導(dǎo)因子a過渡點和CT曲線的選擇。

GH的 Bladed中［9］，將a的過渡點的值取為0.353 9，當BEM理論計算得到a大于該值時，用拋物線經(jīng)驗公式：

NREL的Aerodyn中，將a的過渡點的值取為0.4，當計算得到的a大于該值時，用拋物線經(jīng)驗公式：

也可以將a的過渡點的值取為0.326 2，用線性經(jīng)驗公式：

上面的修正關(guān)系如圖3所示。

圖3 推力系數(shù)修正Fig.3 Modified thrust Factor

如圖3可見修正后的CT曲線在軸向誘導(dǎo)因子a大于過渡點的值之后，繼續(xù)上升，表示在重載風輪面存在很大的壓力差，相較于經(jīng)典BEM理論的動量理論［］計算得到的結(jié)果，更加符合實際情況。

4 算例與結(jié)果

在Fortran環(huán)境下編程進行風力機整機的系統(tǒng)動力學(xué)數(shù)學(xué)建模，將采用譜模擬法得到的湍流風風速時間歷程加載到修正的BEM空氣動力學(xué)模型中，進行風力機風輪氣動載荷計算。并將到風輪氣動載荷輸入基于Kane方法建立的風力機全剛體動力學(xué)模型中，得到風力機響應(yīng)曲線。再將所得的風力機響應(yīng)參數(shù)輸入到基于假設(shè)模態(tài)法離散的動力學(xué)模型中，獲取風力機各部位的變形等信息。同時通過獲取風力機葉片與塔架的質(zhì)量矩陣和剛度矩陣，進行模態(tài)分析。在風輪常用轉(zhuǎn)速下，避開風輪1P和3P頻率［12］，最終實現(xiàn)對風力發(fā)電機整機系統(tǒng)的動力學(xué)進行隨機響應(yīng)分析。

表1 2MW風力發(fā)電機組主要結(jié)構(gòu)參數(shù)Tab.1 Paremeters for 2MW wind turbine

圖4 塔筒模態(tài)Fig.4 Normalized mode of tower

圖5 葉片模態(tài)歸一化Fig.5 Normalized mode of blade

圖6 坎貝爾圖Fig.6 Campbell Diagram

圖7 動力學(xué)分析結(jié)果Fig.7 Results of Dynamic Analysis

本文以企業(yè)合作開發(fā)的2.5MW風力機進行數(shù)值仿真計算，并將數(shù)值計算結(jié)果與Bladed計算結(jié)果進行對比。

圖4～圖5，分別為塔架和葉片的模態(tài)振型。由圖6可知，在風輪常用轉(zhuǎn)速7～14 r/min內(nèi)，塔架與風輪的固有頻率均能避開風輪自激頻率1P和3P。由圖7所示，數(shù)值模擬得到動力學(xué)結(jié)果與Bladed分析結(jié)果比較，使用風力機葉片與塔筒的動力學(xué)響應(yīng)結(jié)果吻合較好，表明該方法能夠較為準確進行風力機動力學(xué)建模與分析。

5 結(jié)論

本文利用Kane方法建立了風力機系統(tǒng)動力學(xué)方程，并根據(jù)假設(shè)模態(tài)法進行離散化，推導(dǎo)出基于柔性梁變形位移場一階完備的動力學(xué)模型。將修正BEM空氣動力學(xué)模型得到風載荷加載到該風力機結(jié)構(gòu)動力學(xué)模型上，進行風力機整機動力學(xué)分析。通過某型2.5MW風力機的數(shù)值計算，并與同權(quán)威計算軟件GH Bladed結(jié)果進行比較，表明使用Kane方法和假設(shè)模態(tài)法建立的風力機動力學(xué)模型充分考慮了風力發(fā)電機組整機結(jié)構(gòu)運動與風速的耦合效應(yīng)，所建立的整機系統(tǒng)動力學(xué)模型更為準確，為開發(fā)具備我國自主知識產(chǎn)權(quán)的水平軸風力發(fā)電機組設(shè)計軟件奠定了基礎(chǔ)。

［1］陳 彥.大型水平軸風力機結(jié)構(gòu)動力響應(yīng)與穩(wěn)定性研究［D］.北京：清華大學(xué)，1999.

［2］竇秀榮.水平軸風力機氣動性能及結(jié)構(gòu)動力學(xué)特性研究［D］.濟南：山東工業(yè)大學(xué)，1995.

［3］Buhl Jr L.A new empirical relationship between thrust coefficient and induction factor for the turbulent windmill state［R］.TechnicalReportNREL/TP-500-36834-August，2005.

［4］蔡國平，洪嘉振.旋轉(zhuǎn)運動柔性梁的假設(shè)模態(tài)方法研究［J］.力學(xué)學(xué)報，2005，37（1）：48-56.

［5］蔣麗忠，趙躍宇.作大范圍運動柔性結(jié)構(gòu)的耦合動力學(xué)［M］.北京：科學(xué)出版社，2007，5.

［6］盛國剛，李傳習，趙 冰.支承運動情況下旋轉(zhuǎn)梁的剛-柔耦合振動分析［J］.振動與沖擊，2006，25（6）：117-120.

［7］ Peter D A.Fast floquet theory and trim for multi-blade in rotorcraft［J］.J.American Helicopter Society，1994，39（4）：82-89.

［8］Jonkman J M.Modeling of the UAE wind turbine for refinement of FAST_AD［R］.NREL/TP-500-34755，National Renewable Energy Laboratory，2003.

［9］ Patrick J，Moriarty A，Hansen C.Aerodyn theory mannal［M］.National Renewable Energy Laboratory，2005.

［10］ Bossanyi E A.Bladed G H.Version 3.67 user manual［M］.Garrad Hassan ＆Partners Limited，2005.

［11］ Burton T，Sharpe D，Jenkins N.Wind energy handbook［M］.England John Wiley ＆ Sons，Ltd，2001.

［12］劉延柱，陳文良，陳立群.振動力學(xué)［M］.北京：高等教育出版社，1998.