陳 嚴(yán)，沈 世，馬新穩(wěn)，劉 雄

（汕頭大學(xué)能源研究所，廣東 汕頭 515063）

0 引言

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

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

1 基礎(chǔ)理論

1.1 BEM 修正模型

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

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

氣流相對速度為：

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

則葉素上的推力和轉(zhuǎn)矩為：

葉片截面的扭矩為：

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

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

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

1.2 Pitt－Peter修正模型

對Pitt－Peter模型做了相應(yīng)修正，使其可以應(yīng)用于風(fēng)力機的空氣動力學(xué)，并推導(dǎo)出在非定長流下誘導(dǎo)因子的控制方程。Pitt－Peter動態(tài)入流模型是采用三個參數(shù)來描述誘導(dǎo)速度在風(fēng)輪盤上的分布，風(fēng)輪上誘導(dǎo)速度的分布和壓力的分布是有聯(lián)系的，其誘導(dǎo)因子求解公式為［3］：

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

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

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

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

此處的加速度與作用力系數(shù)的關(guān)系為：

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

其中［M］是顯在質(zhì)量矩陣，因為他導(dǎo)致誘導(dǎo)速度場延遲一段時間才達(dá)到穩(wěn)定狀態(tài)，理解為速度變化引起的由于質(zhì)量造成的時間延遲。

2 柔性結(jié)構(gòu)反饋

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

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

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

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

對于扭轉(zhuǎn)方向［6］，

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

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

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

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

其中Fc＝mrω2cosβ0 為離心力，eem是截面剛心與質(zhì)心之間的有效距離，F(xiàn)g為重力。

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

3 算 例

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

3.1 風(fēng) 剪

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

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

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

3.2 風(fēng)湍流

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

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

圖2 風(fēng)速隨時間的變化Fig.2 Wind speed step on the time

圖3 半徑32.25m 處的截面在風(fēng)湍流情況下軸向誘導(dǎo)因子的變化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 偏 航

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

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

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

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

3.4 變槳距

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

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

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

3.5 變轉(zhuǎn)速

快速變轉(zhuǎn)速過程，在40s時由額定轉(zhuǎn)速1.26rad／s變?yōu)?.8rad／s，保持到80s，然后轉(zhuǎn)速又快速變回1.26rad／s，仿真時間為120s，運行風(fēng)速為11.4m／s不變。分別用修正BEM 模型、動態(tài)入流模型和加柔性結(jié)構(gòu)反饋的動態(tài)入流模型計算半徑36.35m 處軸向誘導(dǎo)因子的值。

從圖8可以看出，變轉(zhuǎn)速時三種模型計算的誘導(dǎo)樣子值和變槳距時大體趨勢是相同的，也反映出了加柔性結(jié)構(gòu)反饋的動態(tài)入流模型能更好的描述誘導(dǎo)因子的動態(tài)變化。

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

4 結(jié)論

以BEM 模型和Pitt－Peter模型為理論基礎(chǔ)，建立了考慮柔性風(fēng)輪結(jié)構(gòu)變形反饋的動態(tài)入流模型，分別計算了風(fēng)力機在風(fēng)剪、風(fēng)湍流、偏航、葉片槳距角和風(fēng)輪轉(zhuǎn)速變化過程中的誘導(dǎo)速度流場的漸變機理，通過對比不同模型，得出了以下結(jié)論：

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

風(fēng)力機在自然環(huán)境下工作，綜合考慮風(fēng)剪、湍流風(fēng)、偏航和機組本身大范圍的快速運動，如變槳距和變轉(zhuǎn)速，是精確預(yù)測風(fēng)力機動態(tài)載荷的基礎(chǔ)。

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

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

［2］劉吉輝.風(fēng)力機全系統(tǒng)載荷分析及優(yōu)化設(shè)計軟件的集成開發(fā)［D］.汕頭大學(xué)，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［美］.結(jié)構(gòu)動力學(xué)［M］.彭津，王光遠(yuǎn)，等，譯.北京：高等教育出版社，2006.

［6］徐榮橋.結(jié)構(gòu)分析的有限元法和MATLAB 程序設(shè)計［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］漢森［丹］.風(fēng)力機空氣動力學(xué)［M］.肖勁松，譯.北京：中國電力出版社，2009.