基于假設模態法的風力機動力學分析

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

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

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

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

1 風力機多剛體系統動力學模型

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

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

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

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

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

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

對應的廣義慣性力為：

其中：

將式（6）、式（7）帶入式（5），風力發電機組系統動力學方程矩陣形式則為：

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

2 基于假設模態法的風力機動力學模型離散化

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

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

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

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

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

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

改寫成矩陣形式為：

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

2.1 塔筒的廣義質量與廣義剛度

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

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

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

其中與梁的分布剛度相關的勢能分量為：

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

與重力相關的勢能分量為：

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

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

2.2 葉片的廣義質量與剛度

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

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

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

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

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

與梁的分布剛度相關的葉片勢能分量為：

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

由風輪旋轉引起的勢能分量為：

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

經過整理簡化的葉片的廣義剛度為

2.3 基于假設模態表示的葉片和塔筒變形

根據葉片和塔筒的模態振型，可以表示葉片和塔筒的柔性變形位移。

3 修正的BEM空氣動力學模型

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

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

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

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

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

上面的修正關系如圖3所示。

圖3 推力系數修正Fig.3 Modified thrust Factor

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

4 算例與結果

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

表1 2MW風力發電機組主要結構參數Tab.1 Paremeters for 2MW wind turbine

圖4 塔筒模態Fig.4 Normalized mode of tower

圖5 葉片模態歸一化Fig.5 Normalized mode of blade

圖6 坎貝爾圖Fig.6 Campbell Diagram

圖7 動力學分析結果Fig.7 Results of Dynamic Analysis

本文以企業合作開發的2.5MW風力機進行數值仿真計算，并將數值計算結果與Bladed計算結果進行對比。

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

5 結論

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

［1］陳 彥.大型水平軸風力機結構動力響應與穩定性研究［D］.北京：清華大學，1999.

［2］竇秀榮.水平軸風力機氣動性能及結構動力學特性研究［D］.濟南：山東工業大學，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］蔡國平，洪嘉振.旋轉運動柔性梁的假設模態方法研究［J］.力學學報，2005，37（1）：48-56.

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

［6］盛國剛，李傳習，趙 冰.支承運動情況下旋轉梁的剛-柔耦合振動分析［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］劉延柱，陳文良，陳立群.振動力學［M］.北京：高等教育出版社，1998.