柔性風輪的動態入流效應研究

陳 嚴，沈 世，馬新穩，劉 雄

（汕頭大學能源研究所，廣東 汕頭 515063）

0 引言

隨著風力機技術的發展，風力機單機容量迅速增加，相應的風輪直徑從現代風力機技術發展早期的20m 增加到目前的160多m，為了充分利用風力資源和降低風電成本，單機大型化是現代風電機組發展的必然方案。隨著機組大型化，葉片向更大、更柔方向發展，葉片的柔性變形就無法忽略，基于小變形假設的氣彈模型已經不能適應風力機柔性增加所帶來的更多不確定因素的模擬。

考慮葉片柔性的非線性氣動彈性分析已成為研究熱點，并發展了一批有各自特色并在其適用范圍較成熟的分析方法，但現有方法都未能考慮機組柔性對于動態氣動模型的影響及與之相對應的氣動載荷變化對氣動彈性分析的反饋。隨著機組的大型化和近海風場環境的特殊要求，尋求包含風、波及機組結構流固耦合效應和機組非線性大柔性變形效應的整機全系統氣動彈性分析方法成為了重要的研究發展方向。本文基于BEM 修正模型和Pitt－Peter模型，綜合考慮了影響風力機氣動性能的各種因素，并且加入了柔性葉片變形對動態入流模型的反饋，主要考慮揮舞、擺振、扭轉三個自由度的反饋。準確描述了誘導速度隨各個參數改變的動態變化，實現了精確的柔性風輪動態入流效應分析。

1 基礎理論

1.1 BEM 修正模型

在BEM 葉素動量理論的基礎上考慮葉尖損失、輪轂損失、風剪切等因素的修正后，并建立風力機坐標系統［1］，并且考慮偏航，經坐標變換得到入流角和誘導因子計算公式分別為：

其中，θ為槳葉安裝角和截面扭角之和。β0 為錐角（葉片傾角初值），ψ為方位角，η輪轂的仰角。γ＊＝為偏航角。

氣流相對速度為：

由葉素理論得到葉素上的法向力和切向力分別為：

則葉素上的推力和轉矩為：

葉片截面的扭矩為：

作用于塔頂上的力為［2］：

最后求得偏航力矩和傾斜力矩分別為：

本文中的BEM 模型只用來計算誘導因子初始值，為了計算出推力系數、偏航力矩系數、傾斜力矩系數，為后面的模型迭代計算做準備。

1.2 Pitt－Peter修正模型

對Pitt－Peter模型做了相應修正，使其可以應用于風力機的空氣動力學，并推導出在非定長流下誘導因子的控制方程。Pitt－Peter動態入流模型是采用三個參數來描述誘導速度在風輪盤上的分布，風輪上誘導速度的分布和壓力的分布是有聯系的，其誘導因子求解公式為［3］：

其中a0、as、ac分別是風輪處的平均誘導速度、水平（橫向）誘導速度、垂直（縱向）誘導速度。Pitt－Peter給出了推力和力矩系數之間的關系，其合成矩陣的形式為（a）＝［L］（C），即：

其中CT、Cmx、Cmz分別是推力系數、傾斜力矩系數、偏航力矩系數。傾斜力矩系數和偏航力矩系數都可以用在固定偏航時的葉素動量理論求解，三者同時都是一個和時間有關的入流因數的函數。偏航力矩是繞偏航軸（垂直軸）產生的一個使風輪恢復到與風向一致的靜力矩。偏航力矩系數和傾斜力矩系數定義為：

矩陣［L］為流入增益陣，［L］中用尾流傾斜角代替偏航角，那么［L］修正為：

由于自然風不管在強度和方向上都不可能是穩定的，所以很少能夠滿足動量定理的條件。對于偏航中非定常流的情況，風輪盤上的非定常法向加速度分布與式中速度的線性變化具有相同的形式。用流量因數表示為：

此處的加速度與作用力系數的關系為：

與式（a）＝［L］（C）組合可得到完整的運動控制方程［4］為：

其中［M］是顯在質量矩陣，因為他導致誘導速度場延遲一段時間才達到穩定狀態，理解為速度變化引起的由于質量造成的時間延遲。

2 柔性結構反饋

本文重點是考慮柔性風輪，將柔性結構變形反饋到動態入流模型。反饋主要考慮揮舞、擺振和扭轉三個自由度的截面變形對葉片動態氣動力的影響。將葉片簡化成非均勻懸臂梁，采用二結點梁單元對葉片進行有限元離散。整個葉片的多自由度運動方程［5］為：

式中M為質量矩陣，C為阻尼矩陣，K為剛度矩陣，P（t）為載荷單元列陣，x（t）為揮舞／擺振的變形量，θ（t）為扭轉變形量。

對于揮舞、擺振方向［6］，

其中ρe為單元密度，l為單元長度，Ae為單元面積，Ee為彈性模量，Ie為截面慣性矩，Ne為單元所受軸向力。

對于扭轉方向［6］，

其中Je＝?r2dA是梁單元的極慣性矩，Ge為梁單元切變模量。

阻尼矩陣采用Rayleigh阻尼模型［7］：

其中α、β為模型常量，ωm通常取多自由度體系的基頻，而ωn在對動力反應有顯著貢獻的高階振型中選取。通常假設ξm＝ξn＝ξ。

考慮柔性葉片變形，葉素上揮舞、擺振作用力及扭矩分別為：

其中Fc＝mrω2cosβ0 為離心力，eem是截面剛心與質心之間的有效距離，Fg為重力。

由式（34）～式（36）計算葉片上的載荷，代入式（26）多自由度運動方程，計算相應變形量，包括截面扭角變化Δθ、葉片揮舞傾角變化Δβ、擺振傾角σ，其中Δθ反饋回θ，Δβ反饋回β0。經計算擺振方向的反饋較小，故忽略擺振方向的影響。前面的動態入流模型加上結構反饋構成完整的計算模型。

3 算 例

本文考慮大型風力機葉片變形情況下誘導速度的漸變機理，以某5MW 風力機參數作為計算參數［8］，在考慮柔性結構反饋的情況下，分別研究了在風剪、風湍流、偏航、葉片槳距角和風輪轉速變化過程中誘導速度的變化。風力機一般都在野外自然環境下工作，所以考慮以上因素的影響是必要的。由于切向誘導因子對攻角的影響較小，且切向誘導因子可以由軸向誘導因子算出，所以本文中主要研究軸向誘導因子在各個情況下的動態變化。在文獻［1］中已經比較過修正后的BEM 模型和Bladed的計算的各項氣動參數結果，比較一致，故本文中不再和Bladed計算結果進行比較。

3.1 風 剪

在風速恒定為11.4m／s 的情況下，計算半徑32.25m 處軸向誘導因子的值，仿真時間為20s，分別計算有無考慮風剪情況下軸向誘導因子的值，程序計算結果如圖1所示。

圖1 半徑32.25m 處的截面在有無風剪情況下軸向誘導因子的變化Fig.1 Axial induced factor at r＝32.25m for step on the wind shear

從圖1中可以看出兩種情況下軸向誘導因子隨時間都是呈現正玄波周期變化，主要是由于葉片傾角和輪轂仰角所引起的，而且在考慮風剪的情況下軸向誘導因子的波動范圍明顯變大，小型風力機有無考慮風剪影響相對較小。但隨著風力機單機容量迅速增加，相應的風輪直徑從現代風力機技術發展早期的20m 增加到目前的160多m，風剪切的影響也越來越大，所以說對于大型風力機載荷計算時考慮風剪切是非常必要的。在后面的幾種情況中都有考慮風剪切的影響。

3.2 風湍流

在建立三維湍流風的基礎上，風速隨時間變化如圖2所示。分別用修正BEM 模型、動態入流模型和加了柔性結構反饋的動態入流模型進行編程計算，同樣計算半徑32.25m 處軸向誘導因子的值，仿真時間是120s，程序計算結果如圖3 所示。平均風速為11.4m／s，縱向湍流強度為14.87%。

從圖3的計算結果可以看出，用修正BEM 模型計算時，軸向誘導因子是隨湍流風速變化而時刻變化的，且波動劇烈。而用動態入流模型計算時，軸向誘導因子值隨湍流風速的波動較小，這就可以體現出修正BEM 模型和動態入流模型的區別，本文中所建立的動態入流模型所求的誘導因子是和上一時刻所求的誘導因子值有關的，修正BEM 模型只是對單一時刻的誘導因子就行求解。加了結構反饋的動態入流模型計算出的誘導因子的值比原來稍大，整體變化趨勢相同。其中反饋包括截面扭角變化Δθ、葉片揮舞傾角變化Δβ，其中Δθ反饋回θ，Δβ反饋回β0，實際情況中扭矩的方向與葉素扭轉中心、重心和氣動中心的位置有關，本文假設扭矩使翼型順時針方向旋轉，即使翼型抬頭，使截面上扭角減?。?］，扭角減小致使攻角增大，攻角增大反過來又增加扭矩，在扭轉剛度較好的情況下，則不會產生發散的現象，Δβ使得槳距角變小，同樣導致攻角變大，在反復迭代計算后使得軸向誘導因子略大于原來的值，正如圖3中計算結果所示。

圖2 風速隨時間的變化Fig.2 Wind speed step on the time

圖3 半徑32.25m 處的截面在風湍流情況下軸向誘導因子的變化Fig.3 Axial induced factor at r＝32.25m for step on the turbulence

圖4 偏航角隨時間的變化Fig.4 Yaw angel step on the time

3.3 偏 航

風速維持在11.4m／s不變，用加柔性結構反饋的動態入流模型分別計算偏航0°、15°、30°、45°四種情況下半徑43.45m 處軸向誘導因子的值，仿真時間為60s，結果如圖5所示。當偏航角隨時間變化如圖4所示時，用加柔性結構反饋的動態入流模型計算半徑43.45m 處軸向誘導因子的值，仿真時間為60s，結果如圖6所示。

由圖5的計算結果可知，固定偏航下，偏航角在0°至45°之間時，軸向誘導因子隨著偏航角的增加而變大。由圖6動態偏航的計算結果也可以看出軸向誘導因子的變化，且精確計算出了動態態偏航情況下軸向誘導因子的隨偏航角變化而變化的值。圖5和圖6中波動變化主要是考慮了風剪的原因。

圖5 半徑43.45m 處的截面在幾種固定偏航角下軸向誘導因子的變化Fig.5 Axial induced factor at r＝43.45m for step on four different yaw angels

圖6 半徑43.45m 處的截面在動態偏航情況下軸向誘導因子的變化Fig.6 Axial induced factor at r＝43.45m for step on the dynamic yaw

3.4 變槳距

快速變漿距過程，在40s時槳距角快速由0°變成6°，保持到80s，然后又快速變回0°，仿真時間為120s，運行風速為11.4m／s。分別用修正BEM 模型、動態入流模型和加柔性結構反饋的動態入流模型計算半徑32.25m 處軸向誘導因子的值。

從圖7 中可以看出，使用修正BEM 模型計算時，軸向誘導速度是即時變化的，而采用動態入流模型和加柔性結構反饋的動態入流模型時，軸向誘導因子表現出了明顯的滯后，經過一個過渡的過程才達到階躍突變后槳距角所對應的狀態。而在槳距角無變化的時間里誘導因子的值大致相同，只是隨方位角變化幅度不一樣。實際情況中誘導因子也不會出現突變的情況，這也說明動態入流模型更合理的反映了誘導因子的動態變化過程。

圖7 半徑32.25m 處的截面在快速變漿距時軸向誘導因子的變化Fig.7 Axial induced factor at r＝32.25m for step on the pitch angel

3.5 變轉速

快速變轉速過程，在40s時由額定轉速1.26rad／s變為0.8rad／s，保持到80s，然后轉速又快速變回1.26rad／s，仿真時間為120s，運行風速為11.4m／s不變。分別用修正BEM 模型、動態入流模型和加柔性結構反饋的動態入流模型計算半徑36.35m 處軸向誘導因子的值。

從圖8可以看出，變轉速時三種模型計算的誘導樣子值和變槳距時大體趨勢是相同的，也反映出了加柔性結構反饋的動態入流模型能更好的描述誘導因子的動態變化。

圖8 半徑36.35m 處的截面在快速變轉速時軸向誘導因子的變化Fig.8 Axial induced factor at r＝36.35m for step on the rotor speed

4 結論

以BEM 模型和Pitt－Peter模型為理論基礎，建立了考慮柔性風輪結構變形反饋的動態入流模型，分別計算了風力機在風剪、風湍流、偏航、葉片槳距角和風輪轉速變化過程中的誘導速度流場的漸變機理，通過對比不同模型，得出了以下結論：

1）隨著風力機單機容量迅速增加，風剪的影響越來越大，所以對于大型風力機計算載荷時考慮風剪是十分必要的。2）湍流風情況下動態入流模型計算的誘導因子值更符合實際情況，不會隨著風速的突變而即刻變化，加了結構反饋的動態入流模型比原模型計算值稍大，主要原因是結構反饋使得扭角減小，攻角增大。3）固定偏航下，偏航角在0°至45°之間時，軸向誘導因子隨著偏航角的增加而變大。動態偏航的情況下，精確計算出了軸向誘導因子的隨偏航角變化而變化的值。4）在變槳距和變轉速的情況下，采用動態入流模型的情況下，誘導因子都表現出了明顯的滯后，經過一個過渡的過程才達到階躍突變后槳距角和轉速所對應的狀態。

風力機在自然環境下工作，綜合考慮風剪、湍流風、偏航和機組本身大范圍的快速運動，如變槳距和變轉速，是精確預測風力機動態載荷的基礎。

本文的研究內容為大型機組整機氣彈耦合分析做了很好的準備工作，從而對大型柔性風力機進行精確的氣動載荷計算。

［1］劉雄，陳嚴，葉枝全.水平軸風力機氣動性能計算模型［J］.太陽能學報，2005，26（6）：792－800.

［2］劉吉輝.風力機全系統載荷分析及優化設計軟件的集成開發［D］.汕頭大學，2006.

［3］BURTON T，SHARPE D，JENKINS N，et al.Wind energy handbook［M］.John Wiley ＆Sons Ltd，2001.

［4］SUZUKI A.Application of dynamic inflow theory to wind turbine rotors［D］.Salt Lake City：University of Utah，2000.

［5］克拉夫R J［美］.結構動力學［M］.彭津，王光遠，等，譯.北京：高等教育出版社，2006.

［6］徐榮橋.結構分析的有限元法和MATLAB 程序設計［M］.北京：人民交通出版社，2001.

［7］CLOUGH R W，JOSEPH PENZIEN.Dynamics of structures［M］.New York：McGraw Hill，1993.

［8］JONKMAN J，BUTTERFIELD S，MUSIAL W，et al.Definition of a 5MW reference wind turbine for offshore system development［R］.NREL／TP 500－38060.National Renewable Energy Laboratory，Golden，CO，2009.

［9］漢森［丹］.風力機空氣動力學［M］.肖勁松，譯.北京：中國電力出版社，2009.