黃鑫祥，郭小鋒，齊劍峰

（中原工學(xué)院 機(jī)電學(xué)院， 河南 鄭州 450007）

0 引言

水平軸風(fēng)力機(jī)葉片的氣彈變形包含揮舞、擺振和扭轉(zhuǎn)3 個(gè)方面。 常用的風(fēng)力機(jī)動(dòng)力學(xué)分析軟件（如GH-Bladed 和FAST）只考慮了葉片在揮舞和擺振兩個(gè)方向的自由度。 然而，隨著風(fēng)電機(jī)組容量的不斷擴(kuò)大，為降低載荷，大型低風(fēng)速及海上風(fēng)電葉片的柔韌度較早期同類型的風(fēng)電葉片有明顯的增大，如何準(zhǔn)確地分析葉片在實(shí)際工況下的氣彈扭轉(zhuǎn)變形，進(jìn)而分析葉片的動(dòng)態(tài)氣彈扭轉(zhuǎn)變形對風(fēng)力機(jī)氣動(dòng)功率、載荷及氣彈穩(wěn)定性的影響，是大型風(fēng)電葉片的氣彈特性評價(jià)以及氣彈剪裁設(shè)計(jì)中必須解決的問題。

郭小鋒[1]運(yùn)用修正的葉素-動(dòng)量理論和有限元方法，將風(fēng)力機(jī)的動(dòng)力學(xué)問題轉(zhuǎn)化為靜力學(xué)問題， 提出了考慮氣動(dòng)彈性的葉片性能分析方法，綜合分析了葉片扭轉(zhuǎn)變形對風(fēng)力機(jī)氣動(dòng)性能的影響。李德源[2]采用“超級單元”模型，應(yīng)用計(jì)算多體動(dòng)力學(xué)理論和風(fēng)力機(jī)氣動(dòng)模型，建立了受約束的風(fēng)力機(jī)系統(tǒng)氣彈耦合方程。 江術(shù)新[3]建立了整機(jī)在變槳工況下的氣-彈-控聯(lián)合仿真模型，時(shí)域分析了變槳過程中的槳距角變化、功率曲線變化和葉根處受力和力矩的變化。

現(xiàn)有的研究表明，靜氣彈性分析方法沒有反映柔性風(fēng)電機(jī)組的動(dòng)力學(xué)問題本質(zhì)。 風(fēng)力機(jī)的實(shí)際工況具有復(fù)雜多變性，對于風(fēng)力機(jī)葉片在實(shí)際復(fù)雜工況下的氣彈扭轉(zhuǎn)變形問題， 及其對風(fēng)力機(jī)氣動(dòng)性能及氣彈問題性的影響規(guī)律仍須進(jìn)一步研究。

本文基于葉片及塔筒部件的參數(shù)化表達(dá)方法， 采用多體動(dòng)力學(xué)仿真軟件SIMPACK 建立了風(fēng)力機(jī)整機(jī)剛?cè)狁詈夏Ｐ停?提出了一種葉片動(dòng)態(tài)氣彈扭轉(zhuǎn)變形分析的新方法。在幾種典型風(fēng)況下，對NREL5MW 風(fēng)力機(jī)進(jìn)行了氣彈動(dòng)力學(xué)分析，采用matlab 軟件對分析結(jié)果進(jìn)行處理，得到了葉片在典型風(fēng)況下的動(dòng)態(tài)氣彈扭轉(zhuǎn)變形數(shù)據(jù)[4]。 在此基礎(chǔ)上，對NREL5MW 風(fēng)力機(jī)在幾種典型工況下的氣動(dòng)功率、 載荷及氣彈穩(wěn)定性問題進(jìn)行了分析研究。

1 風(fēng)力機(jī)剛?cè)狁詈夏Ｐ偷臉?gòu)建

1.1 參數(shù)化建模

風(fēng)力機(jī)整機(jī)模型是典型的剛?cè)狁詈隙囿w系統(tǒng)。在工作狀態(tài)下，風(fēng)力機(jī)葉片存在明顯的彈性變形，所以將葉片作為柔性體處理。而風(fēng)力機(jī)的其他部件（如機(jī)艙、主軸及輪轂等部件）的剛性大，在工作時(shí)的彈性變形很小，可以看作剛體。為了縮短計(jì)算機(jī)的求解時(shí)間， 本文將風(fēng)力機(jī)模型進(jìn)行了合理簡化，只包括葉片、輪轂、塔筒、機(jī)艙等主要部件，而不考慮齒輪箱、發(fā)電機(jī)等部件。偏航法蘭與塔筒頂端固接，用于連接機(jī)艙，雖然機(jī)艙中不創(chuàng)建主機(jī)架、齒輪箱、發(fā)電機(jī)等模型，但是需要把忽略部件的質(zhì)量和轉(zhuǎn)動(dòng)慣量算入機(jī)艙中。 低速軸與機(jī)艙為轉(zhuǎn)動(dòng)連接，輪轂與低速軸進(jìn)行固接，葉片與輪轂由變槳法蘭連接。

采用EXCEL 表格對葉片的外形和鋪層結(jié)構(gòu)進(jìn)行參數(shù)化表達(dá)，采用matlab 編程調(diào)用PreComp軟件計(jì)算得到SIMPACK 葉片建模中所需各個(gè)截面的屬性參數(shù)（截面質(zhì)心、彈性中心、剪切中心、扭轉(zhuǎn)剛度、彈性模量等），并自動(dòng)生成Simpack 軟件Rotorblade Generation 計(jì)算模塊的輸入文件[5]，[6]。

風(fēng)力機(jī)的控制系統(tǒng)采用基本的變槳變速控制策略， 在SIMPACK 中設(shè)置風(fēng)輪轉(zhuǎn)速及槳距角在不同穩(wěn)態(tài)風(fēng)速下的參數(shù)，同時(shí)通過設(shè)置PID 參數(shù)實(shí)現(xiàn)對風(fēng)輪的動(dòng)態(tài)控制。

在SIMPACK 中， 剛?cè)狁詈隙囿w系統(tǒng)建模是在浮動(dòng)參考坐標(biāo)系法的基礎(chǔ)上建立多體動(dòng)力學(xué)方程， 并由SIMPACK 求解動(dòng)力學(xué)方程獲得風(fēng)力機(jī)的動(dòng)態(tài)響應(yīng)。在多體系統(tǒng)（MBS）理論中，柔性體看作有限元模型若干個(gè)節(jié)點(diǎn)的組合，其中每個(gè)節(jié)點(diǎn)的運(yùn)動(dòng)由跟隨整體坐標(biāo)系的整體平動(dòng)和轉(zhuǎn)動(dòng)以及在局部坐標(biāo)系下發(fā)生小的線性變形組成。 圖1 所示為柔性體的一般模型圖。

為描述節(jié)點(diǎn)的運(yùn)動(dòng)，一般采用慣性坐標(biāo)系下的笛卡爾坐標(biāo)和表示方位的歐拉角。

圖1 柔性體的一般模型Fig.1 General model of a flexible body

1.2 風(fēng)場模型

空氣動(dòng)力學(xué)軟件AeroDyn 是分析風(fēng)力機(jī)葉片風(fēng)載的有力工具， 為了實(shí)現(xiàn)與SIMPACK 軟件的聯(lián)合仿真，需要準(zhǔn)備AeroDyn 的輸入文件、風(fēng)場文件、翼型文件、塔筒文件等[7]。 本文的風(fēng)場文件由NREL 提供的IECWIND 軟件生成。

AeroDyn 基于BEM 理論建立迭代方程計(jì)算風(fēng)輪氣動(dòng)力和誘導(dǎo)速度， 通過迭代求解葉片單元的誘導(dǎo)速度、攻角、推力系數(shù)等參數(shù)，直到誘導(dǎo)因子和入流角的值收斂到最終值[8]。

AeroDyn 與SIMPACK 動(dòng)力學(xué)聯(lián)合仿真流程如圖2 所示。

圖2 動(dòng)力學(xué)仿真流程圖Fig.2 Diagram of simulation process

圖2 中，AeroDyn 與SIMPACK 兩 者 通 過SIMPACK 中241 號力元接口進(jìn)行數(shù)據(jù)交換。AeroDyn 同時(shí)從輸入文件和SIMPACK 求解模塊中收集t 時(shí)刻的葉片幾何形狀、工況條件、葉元速度、位置以及流入風(fēng)等信息，然后利用以上信息計(jì)算出葉片各截面的氣動(dòng)數(shù)據(jù)。 SIMPACK 的求解模塊利用241 號力元接口調(diào)用氣動(dòng)數(shù)據(jù)計(jì)算得到t+dt 時(shí)刻葉片各截面的氣彈變形和氣動(dòng)力， 計(jì)算得到的數(shù)據(jù)又通過241 號力元接口反饋給AeroDyn 用來計(jì)算t+dt 時(shí)刻各截面的氣動(dòng)數(shù)據(jù)，如此往復(fù)，直到仿真完成。而控制模塊根據(jù)輸出的反饋力矩和轉(zhuǎn)速對風(fēng)力機(jī)進(jìn)行變速變槳控制。

為了得到仿真結(jié)果的精確性， 本文選擇SODASRT_2 解算器對模型進(jìn)行離線積分求解，仿真時(shí)長設(shè)為120 s， 求解步長設(shè)為0.01 s。 由于SIMPACK 每個(gè)時(shí)間步長都會(huì)調(diào)用AeroDyn 計(jì)算氣動(dòng)力變化， 所以AeroDyn 求解時(shí)間步長也設(shè)為0.01 s。

為研究陣風(fēng)、 風(fēng)向角度對葉片扭轉(zhuǎn)變形的影響，選取如表1 所示風(fēng)場模型。考慮到海上風(fēng)況的復(fù)雜多變性， 本文在建立風(fēng)場模型時(shí)將方向變化的極端相干陣風(fēng)（EDC）和極端工作陣風(fēng)（EOGR）風(fēng)況下風(fēng)速和風(fēng)向的波動(dòng)設(shè)定為4 個(gè)變化周期，并且假定在EDC 和EOGR 風(fēng)況作用下控制系統(tǒng)故障， 研究風(fēng)力機(jī)在極端風(fēng)況下扭轉(zhuǎn)變形與整機(jī)性能的關(guān)系。

表1 極端風(fēng)況Table 1 Extreme wind conditions

1.3 模型驗(yàn)證

通過SIMPACK 中Eigenvalues 模塊對風(fēng)力機(jī)（NREL5MW）進(jìn)行模態(tài)分析，計(jì)算得到整機(jī)的各階自振頻率（表2）。 由表2 可知，本文的分析結(jié)果同NREL 報(bào)告中的數(shù)據(jù)誤差在5%以內(nèi)[4]，證實(shí)了SIMPACK 所建立風(fēng)力機(jī)整機(jī)模型的可靠性。

表2 5 MW 風(fēng)力機(jī)整機(jī)自振頻率對比Table 2 Comparison of natural frequency of 5 MW wind turbine

2 葉片氣彈扭轉(zhuǎn)變形分析

風(fēng)力機(jī)在啟動(dòng)階段的各種變化均不穩(wěn)定，為更好地研究各風(fēng)況對葉片氣彈變形的影響， 本文假設(shè)各風(fēng)況發(fā)生變化時(shí)風(fēng)力機(jī)處于穩(wěn)定階段。 假定各風(fēng)況在風(fēng)力機(jī)啟動(dòng)后40 s 時(shí)刻發(fā)生變化，并于80 s 之后逐漸恢復(fù)為12 m/s 的穩(wěn)態(tài)風(fēng)。 對風(fēng)力機(jī)整機(jī)進(jìn)行了仿真分析， 得到4 種工況下葉尖揮舞、擺振的動(dòng)態(tài)響應(yīng)（圖3）。

圖3 各風(fēng)況下的葉尖位移Fig.3 Tip deformation under various wind conditions

為分析葉尖部位翼型的氣彈扭轉(zhuǎn)變形數(shù)值，本文首先采用SIMPACK 軟件得到葉尖某截面各節(jié)點(diǎn)的初始坐標(biāo)以及其相對于初始坐標(biāo)的變形量， 將各節(jié)點(diǎn)的初始坐標(biāo)與對應(yīng)的位移量相加得到節(jié)點(diǎn)在各個(gè)時(shí)刻的動(dòng)態(tài)坐標(biāo)值。 將各節(jié)點(diǎn)的動(dòng)態(tài)坐標(biāo)值投影到初始時(shí)刻的截面上， 得到截面扭轉(zhuǎn)變形后的節(jié)點(diǎn)坐標(biāo)值。 對于NREL5MW 葉片，距離葉根60 m 處截面的初始形狀和在48 s 時(shí)的投影形狀對比如圖4 所示。

圖4 在ECD 風(fēng)況下，0 s 和48 s 時(shí)葉尖翼型扭轉(zhuǎn)形變的對比Fig.4 Comparison of tip airfoil torsional deformation at 0 s and 48 s under the ECD wind condition

在風(fēng)力機(jī)的剛?cè)狁詈夏Ｐ椭校?假設(shè)葉片截面在扭轉(zhuǎn)變形時(shí)自身形狀不發(fā)生變化， 葉片截面扭轉(zhuǎn)變形前后兩翼型弦長線間的夾角即為該截面在ECD 風(fēng)況48 s 時(shí)刻的氣彈扭轉(zhuǎn)變形量。

為分析葉片的動(dòng)態(tài)氣彈扭轉(zhuǎn)時(shí)域變形趨勢，由上述方法得到距離葉根60 m 處截面各個(gè)時(shí)刻的動(dòng)態(tài)氣彈扭轉(zhuǎn)變形曲線（圖5）。

圖5 各風(fēng)況下的葉尖扭轉(zhuǎn)變形Fig.5 Tip torsional deformation under various wind conditions

由圖5 可知，當(dāng)風(fēng)力機(jī)運(yùn)轉(zhuǎn)穩(wěn)定后，各風(fēng)況下的葉尖動(dòng)態(tài)扭轉(zhuǎn)角度在2°上下浮動(dòng)，說明葉片具有的扭轉(zhuǎn)剛度使得葉片在一定風(fēng)速下的扭轉(zhuǎn)角度會(huì)在某一個(gè)值上下變動(dòng)。 而葉片各個(gè)截面的扭轉(zhuǎn)剛度不同，導(dǎo)致了葉片各截面扭轉(zhuǎn)角度不同，計(jì)算葉片實(shí)際工作中的扭轉(zhuǎn)角度對葉片的優(yōu)化設(shè)計(jì)具有指導(dǎo)作用。

對比EOGR，EDC 風(fēng)況作用下葉尖扭轉(zhuǎn)變化曲線可知， 陣風(fēng)變化對葉尖扭轉(zhuǎn)變形的影響要大于風(fēng)向角度變化， 其中風(fēng)向角度變化使得扭轉(zhuǎn)變化頻率增加。 在ECD 風(fēng)況作用下，葉尖扭轉(zhuǎn)變形量最為明顯，葉尖最大扭轉(zhuǎn)為4.78 °，最小扭轉(zhuǎn)為0.71 °，陣風(fēng)風(fēng)速越大，葉尖扭轉(zhuǎn)變形量越大， 而風(fēng)速和風(fēng)向角度的共同作用增加了葉片的扭轉(zhuǎn)頻率。

在ECD 風(fēng)況下，由于風(fēng)速和風(fēng)向角度的共同作用造成了葉尖扭轉(zhuǎn)變形的突然增大， 這種短時(shí)間內(nèi)葉尖扭轉(zhuǎn)變化的突然加劇容易對葉片造成致命的損壞。 而風(fēng)向角度多變的風(fēng)況使得葉尖扭轉(zhuǎn)頻率增大，這種風(fēng)況下容易造成葉片的疲勞破壞。

3 氣彈扭轉(zhuǎn)變形對整機(jī)性能的影響

3.1 氣動(dòng)性能

根據(jù)BEM 理論計(jì)算可得，在11.4 m/s 的穩(wěn)態(tài)風(fēng)下， 風(fēng)力機(jī)功率為5 MW。 而在實(shí)際工作過程中，由于塔筒變形、葉片的氣彈變形等影響，風(fēng)力機(jī)功率并不是一個(gè)恒定的值。 不同風(fēng)況的作用對葉片氣彈變形的影響很大， 導(dǎo)致各風(fēng)況下的風(fēng)力機(jī)功率存在明顯波動(dòng)。

風(fēng)力機(jī)控制模塊根據(jù)檢測到的實(shí)時(shí)轉(zhuǎn)速做出變槳響應(yīng)，使其穩(wěn)定在額定功率附近。 仿真時(shí)，在40～80 s，ECD，EOGR 風(fēng)況的風(fēng)速和風(fēng)向角急劇變化，導(dǎo)致風(fēng)力機(jī)變槳速率、葉片翼型攻角變化幅度增大， 葉片氣動(dòng)載荷增加。 隨著葉片變槳速率增大，葉片在風(fēng)速上升期卸載越快，對功率的影響也就越明顯。 另外，當(dāng)控制系統(tǒng)故障時(shí)，風(fēng)力機(jī)變槳控制明顯滯后于實(shí)際風(fēng)況， 在風(fēng)速、 風(fēng)向變化劇烈，但未到達(dá)切出風(fēng)速的極端風(fēng)況下，控制系統(tǒng)無法實(shí)時(shí)調(diào)整到對應(yīng)風(fēng)速的槳距角， 從而導(dǎo)致控制系統(tǒng)暫時(shí)失效， 造成功率波動(dòng)異常。 圖6 為在ECD，EOGR 風(fēng)況下，風(fēng)力機(jī)功率的變化曲線。 由圖6 可知，風(fēng)力機(jī)功率出現(xiàn)了大幅度的波動(dòng)，峰值功率甚至達(dá)到了13 MW。

圖6 各風(fēng)況下的風(fēng)力機(jī)功率變化Fig.6 Variation of wind turbine power under different wind conditions

風(fēng)力機(jī)的功率波動(dòng)主要受陣風(fēng)的影響， 陣風(fēng)變化越大，葉片的氣彈變形越大，風(fēng)力機(jī)功率變化也就越大，在極端風(fēng)況下，功率和載荷波動(dòng)異常，在實(shí)際工作中會(huì)造成風(fēng)力機(jī)的損壞。

3.2 載荷波動(dòng)

圖7 為各風(fēng)況下輪轂處載荷的變化曲線。 由圖7 可知： 相比于EDC 風(fēng)況，ECD 和EOGR 風(fēng)況有非常明顯的載荷波動(dòng)， 最大輪轂推力值分別為1 280 kN 和1 175 kN；ECD 風(fēng)況作用下的輪轂處最大載荷要小于EOGR 風(fēng)況，葉片的扭轉(zhuǎn)變形在一定程度上能削減風(fēng)力機(jī)上的極限載荷。

圖7 各風(fēng)況下輪轂處載荷的變化Fig.7 Load changes at the hub under different wind conditions

3.3 氣彈穩(wěn)定性

在復(fù)雜工況下，風(fēng)速、風(fēng)向的劇烈變化引起了葉片揮舞、擺振和扭轉(zhuǎn)方向的急劇變化，從而造成葉片氣彈失穩(wěn)。 特別是葉片的扭轉(zhuǎn)變形直接影響攻角變化，從而造成升力系數(shù)的變化，嚴(yán)重時(shí)會(huì)發(fā)生失速現(xiàn)象。 本文通過AeroDyn 對加載風(fēng)載的葉片進(jìn)行計(jì)算分析， 選取50～60 s 時(shí)刻各風(fēng)況下距離葉根60 m 處截面翼型的氣動(dòng)數(shù)據(jù)進(jìn)行處理，得到攻角-升力系數(shù)曲線，并與靜態(tài)翼型攻角-升力系數(shù)曲線進(jìn)行比較（圖8）。

圖8 各風(fēng)況下葉尖截面攻角-升力系數(shù)曲線Fig.8 Angle of attack-lift coefficient curve of blade tip section under various wind conditions

由圖8 可知：NWP，EDC 風(fēng)況下的升力系數(shù)隨攻角呈線性變化，且該時(shí)間段內(nèi)EDC 風(fēng)況下的攻角大于NWP 風(fēng)況下的攻角，但沒達(dá)到翼型的靜態(tài)失速攻角；EOGR 風(fēng)況下的攻角進(jìn)一步增大，接近翼型靜態(tài)失速攻角，仍呈線性變化；ECD 風(fēng)況下的攻角超過靜態(tài)失速攻角并繼續(xù)增大， 升力系數(shù)也隨著攻角的增大總體呈上升趨勢，達(dá)到峰值后攻角逐漸減小，升力系數(shù)也隨之減小，曲線開始回歸線性區(qū)，這是典型的動(dòng)態(tài)失速現(xiàn)象[9]；在ECD 風(fēng)況作用下，葉尖截面存在明顯的動(dòng)態(tài)失速現(xiàn)象，而其他風(fēng)況作用下，葉尖截面均未出現(xiàn)失速現(xiàn)象，所以在風(fēng)速、風(fēng)向變化幅度大的風(fēng)況下，葉片容易出現(xiàn)動(dòng)態(tài)失速現(xiàn)象，從而造成葉片氣彈失穩(wěn)， 嚴(yán)重時(shí)葉片會(huì)發(fā)生顫振甚至斷裂。

4 結(jié)論

本文基于SIMPACK 軟件建立了剛?cè)狁詈系娘L(fēng)力機(jī)整機(jī)模型，對NREL5MW 風(fēng)力機(jī)進(jìn)行了動(dòng)力學(xué)分析研究，得到以下結(jié)論。

①提出了一種分析葉片動(dòng)態(tài)氣彈扭轉(zhuǎn)變形的新方法， 該方法可得到葉片在復(fù)雜風(fēng)況下的動(dòng)態(tài)氣彈扭轉(zhuǎn)變形數(shù)據(jù)。

②極端復(fù)雜工況下， 風(fēng)速和風(fēng)向的變化會(huì)導(dǎo)致葉片產(chǎn)生較大的氣彈變形，進(jìn)而對風(fēng)力機(jī)氣動(dòng)、載荷及氣彈穩(wěn)定性產(chǎn)生較大的影響。

③在ECD 工況下， 葉片氣彈扭轉(zhuǎn)變形急劇增大，會(huì)導(dǎo)致葉片處于短暫失速狀態(tài)，葉片容易出現(xiàn)失速現(xiàn)象。 風(fēng)力機(jī)功率和載荷波動(dòng)增大，從而造成葉片氣彈失穩(wěn)，嚴(yán)重時(shí)葉片會(huì)發(fā)生顫振甚至斷裂。