倪晨鋒，李德源，汪顯能，池志強(qiáng)

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

10MW風(fēng)力機(jī)葉片設(shè)計及其在隨機(jī)風(fēng)載荷下的響應(yīng)分析*

倪晨鋒，李德源?，汪顯能，池志強(qiáng)

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

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

風(fēng)力機(jī)葉片；超級單元；葉素動量理論；時域響應(yīng)

0　 引 言

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

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

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

1　 葉片的氣動外形設(shè)計和結(jié)構(gòu)設(shè)計

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

1.1葉片的氣動外形設(shè)計

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

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

由于葉片截面翼型直接決定了氣動載荷的大小，因此翼型的選擇至關(guān)重要?？紤]到葉根部分與輪轂連接的結(jié)構(gòu)要求，以及葉片氣動載荷主要集中在中部到葉尖部分，所以，選擇了具有一定尾緣厚度的 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復(fù)合材料葉片的內(nèi)部結(jié)構(gòu)設(shè)計

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

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

圖4　剛體質(zhì)量隨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　 柔性葉片多體動力學(xué)模型的建立

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

2.1葉片的多體模型

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

2.2模型中剛體的質(zhì)量分布與彈簧剛度系數(shù)的確定

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

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

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

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

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

葉片在極端載荷作用下的強(qiáng)度分析工作已經(jīng)在文獻(xiàn)［1］中完成，驗(yàn)證了該葉片靜強(qiáng)度和剛度。

2.3帶約束的葉片多體動力學(xué)方程及其數(shù)值求解

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

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

3　 葉片氣動載荷分析

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

葉素動量理論

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

葉素理論沒有考慮氣流的展向流動，將構(gòu)成葉片展向的多個葉素看成相互獨(dú)立的二維翼型。風(fēng)輪繞轉(zhuǎn)動中心旋轉(zhuǎn)一周受到的推力T及轉(zhuǎn)矩Q為：

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

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

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

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

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

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

4　 葉片的隨機(jī)風(fēng)載荷響應(yīng)分析

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

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

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

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

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

圖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)定和紊流風(fēng)速下某一截面攻角隨時間的變化情況?？梢钥吹剑谌~片運(yùn)行穩(wěn)定后，穩(wěn)定風(fēng)速條件下攻角的波動比紊流風(fēng)速下攻角的波動要小，對氣動力的影響也會更小。

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

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

表3　兩種風(fēng)速下的數(shù)據(jù)對比Table 3 The contrast of the data in the two kinds of wind

5　 結(jié) 論

（1）在綜合考慮風(fēng)力機(jī)的氣動性能與結(jié)構(gòu)強(qiáng)度基礎(chǔ)上，設(shè)計了10 MW風(fēng)力機(jī)葉片的氣動外形和內(nèi)部結(jié)構(gòu)。建立了柔性葉片的多剛體動力學(xué)模型，結(jié)合多體動力學(xué)建模方法和風(fēng)力機(jī)氣動分析模型，構(gòu)建了柔性葉片非線性流固耦合方程，并開發(fā)了相應(yīng)的數(shù)值分析軟件，能夠分析在隨機(jī)風(fēng)速作用下，實(shí)現(xiàn)在葉片實(shí)時變形位置分析葉片的氣動載荷，并獲得葉片的動力學(xué)響應(yīng)情況。

（2）分析兩種風(fēng)況條件下葉片的振動響應(yīng)和氣動載荷，并進(jìn)行了對比分析，結(jié)果表明所設(shè)計的葉片滿足設(shè)計要求。采用所提出的葉片氣彈耦合分析方法，能夠較為有效地展現(xiàn)葉片的振動與氣動載荷之間的耦合關(guān)系。

［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］ 李德源， 莫文威， 嚴(yán)修紅， 等. 基于多體模型的水平軸風(fēng)力機(jī)氣彈耦合分析［J］. 機(jī)械工程學(xué)報， 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］ 李德源， 莫文威， 夏鴻建， 等. 水平軸風(fēng)力機(jī)柔性葉片氣彈耦合分析［J］. 太陽能學(xué)報， 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)動力學(xué)［M］. 北京: 科學(xué)出版社，2008.

［13］ 洪嘉振. 計算多體系統(tǒng)動力學(xué)［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-），男，教授，碩士生導(dǎo)師，主要從事大型風(fēng)力機(jī)氣動與結(jié)構(gòu)分析、風(fēng)力機(jī)系統(tǒng)測試和計算機(jī)軟件的開發(fā)與應(yīng)用等方面的科研工作。

TK83

A

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

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

2016-03-03

2016-04-17

國家自然科學(xué)基金項(xiàng)目（51276043）

倪晨鋒（1991-），男，碩士研究生，主要從事風(fēng)力機(jī)氣彈耦合研究及風(fēng)力機(jī)結(jié)構(gòu)設(shè)計。