10MW風力機葉片設計及其在隨機風載荷下的響應分析*

倪晨鋒，李德源，汪顯能，池志強

（廣東工業(yè)大學機電工程學院，廣州 510006）

10MW風力機葉片設計及其在隨機風載荷下的響應分析*

倪晨鋒，李德源?，汪顯能，池志強

（廣東工業(yè)大學機電工程學院，廣州 510006）

對大型風力機柔性葉片的設計方法及其在隨機風載荷作用下的動態(tài)響應與載荷特性進行了研究。根據風力機葉片空氣動力學和結構設計理論，將柔性葉片離散為多個剛體，形成一個多體系統(tǒng)。根據多體動力學的建模方法和葉片氣動模型，考慮兩者的相互作用，建立了柔性葉片的非線性耦合動力學方程并開發(fā)了相應的仿真程序。算例分析了葉片在隨機風載荷作用下的氣彈載荷與隨機振動響應，并對穩(wěn)定風速和紊流風速下的響應結果作了對比分析。

風力機葉片；超級單元；葉素動量理論；時域響應

0　 引 言

目前10 MW的巨型風電機組已成為風力機領域競相研發(fā)的對象［1］。由于葉片一端與輪轂相連，為類似懸臂梁的細長結構，隨著單機功率的增加，葉片愈加細長，柔性和變形越來越大，變形和氣動載荷的耦合以及所引起的風電機組氣動彈性問題將更加突出。

柔性葉片動力響應分析的困難在于，一方面葉片較大的彈性變形和機械振動會反作用于葉片的空氣動力載荷，引起氣彈耦合；另一方面風載荷是一個隨機變量。響應分析涉及到氣動載荷模型和結構動力模型。在氣動模型方面，葉素動量理論（blade element momentum， BEM）計算效率較高且結果比較符合實際，在風力機設計領域得到了廣泛應用［2］；而對于風力機葉片這類具有較大彈性變形和大范圍空間運動構件的動力學建模，應用傳統(tǒng)有限元數值進行葉片的氣彈耦合計算還有一定困難。ZHAO［3］及 HOLIERHOEK［4］等引入“超級單元”（super-element）將風力機葉片離散為若干個通過運動副、彈簧和阻尼器聯(lián)接的剛體，通過多體動力學理論可以用較少的自由度建立系統(tǒng)的非線性動力學方程。文獻［5］運用多體系統(tǒng)動力學的建模方法，建立了某5 MW風力機整機的剛-柔多體氣彈耦合方程，對整機的振動模態(tài)和氣彈耦合響應進行了數值模擬。

本文針對文獻［1］提出的10 MW概念型風力機葉片，應用“超級單元”模型，將葉片離散為6個超級單元，基于多體系統(tǒng)動力學理論，結合葉素動量理論，通過在翼型的入流角和攻角中引入葉片的變形與振動來修正，建立葉片氣彈耦合方程，應用Matlab編制動力學方程程序化的建模程序和數值求解程序，得出在穩(wěn)定風速和紊流風速下的氣彈耦合響應。為葉片結構和氣動性能的進一步優(yōu)化設計提供有效的分析手段。

1　 葉片的氣動外形設計和結構設計

文獻［1］風電機組是一臺三葉片、上風向的水平軸風電機組，采用變速變槳控制。額定風速取11.4 m/s，材料為玻璃鋼。

1.1葉片的氣動外形設計

風載荷作用在風輪葉片上，推動風輪轉動，將機械能轉化為電能。而葉片所受氣動載荷的大小主要決定于葉片的外形，所以葉片氣動設計主要是外形的設計。內容主要包括翼型的選擇，以及弦長、扭角分布的確定。

設計葉片外形時，氣動性能的分析通常采用葉素動量理論，同時考慮多種影響因素，如葉尖和輪轂損失修正、攻角的修正，并考慮氣流沿葉片軸向的流動，對軸向誘導因子的修正，建立風電機組氣動設計數學模型，綜合考慮機組的功率系數，設計完成的葉片弦長和扭角分布。

由于葉片截面翼型直接決定了氣動載荷的大小，因此翼型的選擇至關重要?？紤]到葉根部分與輪轂連接的結構要求，以及葉片氣動載荷主要集中在中部到葉尖部分，所以，選擇了具有一定尾緣厚度的 FFA-W 系列翼型，分別是 FFA-W3-241、FFA-W3-301、 FFA-W3-360、 FFA-W3-480 和FFA-W3-600，最大相對厚度分別為24.1%、30.1%、36%、48%和60%，靠近葉尖部分選擇了最大相對厚度為24.1%。最終得到葉片全長為86.366 m，葉片截面弦長、扭角和厚度隨r（r表示葉片的截面到葉根距離）的分布情況分別如圖1、圖2和圖3所示［1］。

圖1　弦長隨r分布情況Fig. 1 The distribution of the chord along with r

圖2　扭角隨r分布情況Fig. 2 The distribution of the twist along with r

圖3　厚度隨r分布情況Fig. 3 The distribution of the thickness along with r

1.2復合材料葉片的內部結構設計

復合材料葉片的結構從作用上可以分為蒙皮和主梁兩個部分。蒙皮形成葉片外形，同時承受部分氣動與機械載荷（如重力等）。主梁是主要的承載構件，通常由梁帽和腹板組成，梁帽承擔氣動和機械載荷引起的彎曲載荷，腹板置于梁帽空腔內，支撐梁帽，以提高剛度，防止結構屈曲，梁帽一般采用單向布制作，而腹板采用夾層結構制作。葉片的大部分重量集中在主梁上。

文獻［1］的葉片鋪層中，梁的結構采用了矩形梁形式，復合材料采用了玻璃鋼材料。葉片總質量為41 716 kg，結構阻尼比為3%。葉片質量和截面剛度隨r的分布情況分別如圖4、圖5所示。

圖4　剛體質量隨r分布情況Fig. 4 The mass of the rigid body along with r

圖5　剛體截面剛度隨r分布情況Fig. 5 The sectional stiffness of the rigid body along with r

2　 柔性葉片多體動力學模型的建立

將柔性葉片離散為若干具有一定質量與形狀的體，剛體之間用萬向節(jié)、轉動副和約束力元連接，形成一個多體系統(tǒng)。

2.1葉片的多體模型

將葉片離散為6個超級單元來建立柔性葉片的拓撲模型。每個超級單元包含4個剛體，前后兩個超級單元相連的兩個剛體剛性連接（即將前一個超級單元的最后一個剛體與后一個超級單元的第一個剛體相連接），合并為一個剛體，形成19個剛體，總共31個自由度的多體系統(tǒng)［5］。

2.2模型中剛體的質量分布與彈簧剛度系數的確定

根據葉片質量分布圖（圖4），結合超級單元模型中剛體長度，可以提取劃分的各剛體質量，如表1所示。彈簧的剛度系數可以通過把葉片簡化成懸臂梁模型來確定，如果一個超級單元的長度為L，則剛體B1和B4長為kL，k（k ＜ 1）為長度系數，剛體B2和B3長為（1 - 2kL）/2。根據彈性梁的彎曲及扭轉理論，計算出彈簧的剛度系數［3］，計算表達式如下：

式中：E為彈性模量，G為剪切模量，I為截面慣性矩，1xC 為扭轉彈簧剛度系數，其余為彎曲彈簧剛度系數。由于長度系數k會對固有頻率的近似解產生影響，根據RAUH［6］的結論：長度系數為1/5 ≤ k ≤ 1/4時，則可用較少的超級單元來逼近固有頻率精確解，本文取k為0.217 6。根據式（1）～ 式（5）和葉片的截面數據（表 1）得到葉片上各剛體之間的彈簧剛度系數如表2。

表1　葉片截面屬性Table 1 Data of the blade section

表2　葉片彈簧剛度系數 / （N·m/rad）Table 2 Spring stiffness coefficients of the blade / （N·m/rad）

根據以上參數，在動力學分析軟件adams中建立葉片多體模型，在 adams/vibration振動分析模塊中進行模態(tài)計算，得到前六階振型及其對應的固有頻率，其結果與DTU實驗室［1］提供的結果基本吻合，驗證了超級單元方法建模的正確性。

葉片在極端載荷作用下的強度分析工作已經在文獻［1］中完成，驗證了該葉片靜強度和剛度。

2.3帶約束的葉片多體動力學方程及其數值求解

由上所述，離散的葉片多體系統(tǒng)由多個剛體通過運動副連接，取各運動副的相對運動為廣義坐標θ（t），則由此多剛體系統(tǒng)的樹結構形式可得系統(tǒng)的廣義坐標矩陣為：

根據多體系統(tǒng)動力學普遍方程，可以導出葉片系統(tǒng)的非線性動力學微分方程組：θ為加速度列陣，推導過程與表達形式可參閱文獻［13］。葉片在額定轉速下轉動，確定的轉速相當于運動約束，表達為一代數方程，引入拉格朗日乘子λ，結合方程（7），構成系統(tǒng)封閉的微分-代數方程組。采用違約穩(wěn)定法處理，采用數值積分可以求出廣義加速度

3　 葉片氣動載荷分析

作用在葉片上的空氣動力載荷分析，一般根據葉素動量理論，但是，由于葉片的變形和振動（揮舞、擺振和扭轉振動）反過來會影響氣動載荷，考慮流固耦合，本文對傳統(tǒng)的氣動理論分析迭代過程進行了修正。

葉素動量理論

傳統(tǒng)的葉素動量理論通過引入軸向和切向誘導因子，綜合考慮了氣流的三維流動效應，比較準確地表達了作用在葉片上的氣動載荷，在風力機領域應用廣泛。計算流固耦合，本文通過在翼型的入流角和攻角中引入葉片的變形與振動來修正模型，使模型更加完善。

葉素理論沒有考慮氣流的展向流動，將構成葉片展向的多個葉素看成相互獨立的二維翼型。風輪繞轉動中心旋轉一周受到的推力T及轉矩Q為：

其中：B為葉片數；c是翼型弦長；ρ為空氣密度；W是微段dr的轉速rω（ω為葉片轉動角速度）與來流風速∞U 的合速度；LC 、DC 分別是升力、阻力，隨攻角α的變化而變；L為剛體沿葉片徑向的長度；r為翼型截面距離輪轂中心的距離；φ為入流角。

而動量理論通過引入軸向誘導因子a和切向誘導因子a′來模擬氣流的展向流動，同時考慮葉尖與輪轂渦處復雜的流動，推力和轉矩分別表達為：

其中，F為綜合考慮葉尖損失和輪轂損失的修正因子。各氣動參數對葉素產生氣動載荷的作用機理如下圖所示。由前面葉素理論的公式可知，確定入流角即可得到葉素上的氣動力。

考慮氣彈耦合，葉片變形和振動對氣動載荷的影響體現在對入流角φ及攻角α的影響上，修正的計算公式為：

其中，θ為各剛體的扭轉角；opeV-為翼型垂直于旋轉平面（揮舞）速度；ipeV-為翼型平行于旋轉平面（擺振）速度。由上述式（8）～式（13），參考DTU實驗室［1］提供的二維翼型升阻數據，通過迭代即可算出各氣動參數。

圖6　剛體質心截面上的氣動參數及氣動中心連體基Fig. 6 Aerodynamic variables of center-mass section and body-fixed coordinate system of aerodynamic center of rigid body

4　 葉片的隨機風載荷響應分析

氣彈耦合仿真首先在求解動力學微分方程的過程中，調入氣動模塊算出該時刻的氣動載荷，將其反饋到動力學模塊中，分析得出下一時刻剛體的運動參數并再代入氣動模塊中，從而可計算出下一時間步的氣動載荷，如此循環(huán)下去，最后達到設定的仿真時間后結束。氣彈耦合分析流程圖參考文獻［7］。

分析中，葉片以恒定的轉速繞輪轂中心轉動，這樣建立約束方程限制葉片轉速：

其中，廣義坐標q1為葉根剛體相對于風輪轉動中心線的轉動；ωe為額定轉速，9.6 r/min［1］。取α=β=10，仿真時間為50 s，積分精度為0.000 1，輸出積分步長0.001 s。

以仿真分析了設計的10 MW葉片在穩(wěn)定風速和紊流風速下（這里風速為11.4 m/s）葉片扭轉變形和振動對氣動載荷的影響。對比了兩種風速下，葉尖的位移和葉根的位移以及各剛體的攻角變化對氣動載荷的影響。

圖7是在穩(wěn)定風速和紊流風速下葉片第1個剛體（葉根位置）位移隨時間的變化曲線。圖8是在穩(wěn)定和紊流風速下第19個剛體（葉尖位置）力矩隨時間的變化情況。從葉尖位移變化圖可以看出，在葉片進入穩(wěn)定運行后，穩(wěn)定風速下葉尖揮舞位移按照一定的周期變化，位移范圍在6.3 ～ 7.3 m；而紊流風速下葉尖揮舞位移變化沒有規(guī)律波動較大，變化范圍4.9 ～ 8.1 m。這兩種風速下葉尖的最大位移都沒有超過葉片與塔架的距離18.26 m［1］（防止葉片在運行過程中與塔架相撞）。

圖7　第1個剛體揮舞與擺振力矩Fig. 7 The flapwise and edgewise bending moment of the first rigid body

圖8 第19個剛體揮舞與擺振方向位移Fig. 8 The flapwise and edgewise displacement of the 19th rigid body

圖9是在穩(wěn)定和紊流風速下某一截面攻角隨時間的變化情況?？梢钥吹剑谌~片運行穩(wěn)定后，穩(wěn)定風速條件下攻角的波動比紊流風速下攻角的波動要小，對氣動力的影響也會更小。

圖9　某一截面攻角隨時間的變化情況Fig. 9 The curves of time versus varying angle of attack of a cross-section

從葉根力矩圖（圖 7）同樣可以看到，在穩(wěn)定風速下的葉根力矩的波動比紊流風速下的葉根力矩要小，具體的對比結果如下表 3。由于在兩種不同的風速下，葉片的扭轉速度和扭轉變形不同，導致攻角有很大的差異，使得氣動載荷也不相同，當然最終得到的葉根力矩也會有很大的差別。因此，對于大型風力機，需要考慮葉片的振動和扭轉對氣動載荷的耦合。

表3　兩種風速下的數據對比Table 3 The contrast of the data in the two kinds of wind

5　 結 論

（1）在綜合考慮風力機的氣動性能與結構強度基礎上，設計了10 MW風力機葉片的氣動外形和內部結構。建立了柔性葉片的多剛體動力學模型，結合多體動力學建模方法和風力機氣動分析模型，構建了柔性葉片非線性流固耦合方程，并開發(fā)了相應的數值分析軟件，能夠分析在隨機風速作用下，實現在葉片實時變形位置分析葉片的氣動載荷，并獲得葉片的動力學響應情況。

（2）分析兩種風況條件下葉片的振動響應和氣動載荷，并進行了對比分析，結果表明所設計的葉片滿足設計要求。采用所提出的葉片氣彈耦合分析方法，能夠較為有效地展現葉片的振動與氣動載荷之間的耦合關系。

［1］ BAK C， ZAHLE F， BITSCHE R， et al. Description of the DTU 10MW reference wind turbine［R］. DTU Wind Energy Report-I-0092， 2013.

［2］ HANSEN M O L， S?RENSEN J N， VOUTSINAS S， et al. State of the art in wind turbine aerodynamics and aeroelasticity［J］. Progress in aerospace sciences， 2006，42（4）: 285-330. DOI: 10.1016/j.paerosci.2006.10.002.

［3］ ZHAO X Y， MAI?ER P， WU J Y. A new multibody modelling methodology for wind turbine structures using a cardanic joint beam element［J］. Renewable energy，2007， 32（3）: 532-546. DOI: 10.1016/j.renene.2006.04.010.

［4］ HOLIERHOEK J G. Aeroelasticity of large wind turbines［D］. Delft: Delft University of Technology， 2008.

［5］ 李德源， 莫文威， 嚴修紅， 等. 基于多體模型的水平軸風力機氣彈耦合分析［J］. 機械工程學報， 2014， 50（12）: 140-150. DOI: 10.3901/JME.2014.12.140.

［6］ RAUH J， SCHIEHLEN W. Various approaches for the modeling of flexible robot arms［C］. Proceedings of the Euromech-Colloquium 219 on Refined Dynamical Theories of Beams， Plates and Shells and Their Applications. Berlin Heidelberg: Springer， 1987: 420-429. DOI: 10.1007/978-3-642-83040-2_36.

［7］ 李德源， 莫文威， 夏鴻建， 等. 水平軸風力機柔性葉片氣彈耦合分析［J］. 太陽能學報， 2015， 36（3）: 734-742.

［8］ HANSEN M O L. Aerodynamics of wind turbines［M］. 2nd ed. Denmark: Routledge， 2007.

［9］ HOLIERHOEK J G. Rigid body modeling of wind turbines for aeroelastic stability investigations［C］//45th AIAA Aerospace Sciences Meeting and Exhibit. Reno，Nevada: American Institute of Aeronautics and Astronautics Inc， 2007: 2485-2497. DOI: 10.2514/6.2007-210.

［10］ LI L， LI Y H， LIU Q K， et al. Flapwise non-linear dynamics of wind turbine blades with both external and internal resonances［J］. International journal of non-linear mechanics， 2014， 61: 1-14. DOI: 10.1016/j.ijnonlinmec. 2014.01.006.

［11］ COX K， ECHTERMEYER A. Structural design and analysis of a 10MW wind turbine blade［J］. Energy procedia， 2012， 24: 194-201. DOI: 10.1016/j.egypro. 2012.06.101.

［12］ 齊朝暉. 多體系統(tǒng)動力學［M］. 北京: 科學出版社，2008.

［13］ 洪嘉振. 計算多體系統(tǒng)動力學［M］. 北京: 高等教育出版社， 1999.

Design of the10 MW Wind Turbine Blade and Its Response Analysis under Random Wind Loads

NI Chen-feng， LI De-yuan， WANG Xian-neng， CHI Zhi-qiang
（School of Electromechanical Engineering， Guangdong University of Technology， Guangzhou 510006， China）

This paper presented a design method of the flexible blade and its dynamic response under random wind. The flexible blade was discretized into multiple rigid bodies， and the multi-body system was established according to the theory of aerodynamics and the structural design of the blade. The nonlinear aeroelastic coupling equations of a constrained flexible blade are derived by the dynamics of multi-body system theory and the blade aerodynamics model， and the simulation program was developed. The examples analyzed time-domain responses of the blade under stable and turbulent wind speed， and the results were compared under these two kinds of wind speed.

wind turbine blade； super-element； blade element momentum； time-domain response

李德源（1965-），男，教授，碩士生導師，主要從事大型風力機氣動與結構分析、風力機系統(tǒng)測試和計算機軟件的開發(fā)與應用等方面的科研工作。

TK83

A

10.3969/j.issn.2095-560X.2016.03.007

2095-560X（2016）03-0206-07

2016-03-03

2016-04-17

國家自然科學基金項目（51276043）

倪晨鋒（1991-），男，碩士研究生，主要從事風力機氣彈耦合研究及風力機結構設計。