常 林，劉廷瑞

（山東科技大學 機械電子工程學院，山東 青島 266590）

風能是一種具有較好發(fā)展前景的新型清潔能源[1]。風力機是捕獲、轉化風能的主要部分，傳統(tǒng)風力機形式粗大，葉片剛度和強度較高，但是近年來隨著材料科學和機械設計分析、制造水平的快速提升，風力機葉片的形式變得越來越細長，雖然風能捕獲量得到了很大提高，隨之而來的是剛度和強度較傳統(tǒng)葉片降低了很多。葉片在高風速或達到失速條件時發(fā)生的動態(tài)不穩(wěn)定現(xiàn)象稱為失速顫振。葉片一旦發(fā)生顫振，帶來的破壞是巨大的[2]。近年來發(fā)生在世界各地風場的風力機在大風速條件下的斷裂失效也印證了此問題的廣泛存在。2014年發(fā)生在湛江勇士風場的臺風事故中[3]，風場33臺風力機在臺風“威馬遜”的影響下幾乎被全部吹倒，有5臺風機完全損壞，其中葉片破壞是一種重要的破壞形式，且風力機葉片成本占風機總成本的15%～20%，所以葉片的保護至關重要。

風力機葉片所受載荷在氣動力、彈性力、重力影響下發(fā)生耦合，極易在大風速條件下產(chǎn)生氣動彈性不穩(wěn)定現(xiàn)象。失速顫振作為氣彈不穩(wěn)定的一種重要形式，是葉片破壞的重要原因之一。葉片一旦發(fā)生失速顫振，便會表現(xiàn)出振幅不衰減的高頻振動，而這種振動會在極短時間內(nèi)使葉片失效。水平軸風力機葉片顫振分析涉及葉片氣動彈性穩(wěn)定性，即葉片在氣動力作用下自身發(fā)生彈性變形，而這種變形反過來又影響周圍流體，因此葉片上氣動載荷在建模分析時應考慮其具有非線性[4]。此外，顫振分析還通常涉及彈性葉片的結構特征以及葉片附近非定常氣流場的描述以及相互耦合的機理分析。所以葉片的顫振研究十分復雜。

葉片常用分析模型有彈性鉸鏈模型、有限元模型、典型截面分析模型等。典型截面模型是一種在葉片氣彈分析中經(jīng)常采用的簡化型葉片結構模型，典型截面模型可以簡單、經(jīng)濟和快捷地獲得失速顫振特性的主要結構和氣動影響參數(shù)，是研究葉片揮舞/擺振氣彈穩(wěn)定性常用的一種有效結構模型，易于建模和模擬仿真。文獻[5]基于定值轉角下非線性氣動力模型分析了揮舞-扭轉以及揮舞-擺振的氣彈不穩(wěn)定問題，但過程未作線性化處理。文獻[6]采用2D截面分析模型研究了揮舞-擺振耦合運動，考慮了結構阻尼以及耦合結果的穩(wěn)定性。文獻[7]將葉片簡化為懸臂梁，考慮彎-彎耦合，在彈性梁建模的基礎上引入材料阻尼從而得到運動方程。但在非線性氣動載荷、大攻角下的失速顫振方面研究相對較少，且均未考慮風力機在實際工作過程中的塔影效應和風切變的影響。本文在課題組前期研究的基礎上[8]，考慮風速的空間分布，主要基于風切變和塔影效應的影響，模擬非線性氣動力的大攻角情形，進行葉片運動建模，提出MPC控制策略，抑制葉片非線性顫振。MPC控制策略與氣彈模型結合的方法具有較高的可操作性且該方式易于實現(xiàn)，可以得到很好的控制效果；基于風切變和塔影效應的影響的研究更具現(xiàn)實意義，分析過程佐證了該方法不失一般性。

1 風速建模

風力機在實際工作過程中，所受風載并非是各處相同的，在方位角不斷變化的情況下，葉片各點的作用風速也在變化，此處主要考慮兩點影響：風切變效應和塔影效應。風切變是指風在運動過程中，因地面摩擦和熱力因素的存在，使得作用在風力機上的風的速度隨著高度的升高而下降的現(xiàn)象。由于風切變的存在，葉輪在運轉過程中各處受載不均，甚至會出現(xiàn)葉片振動和功率波動，影響風力機的正常工作。圖1為三葉片水平軸風力機示意圖。

圖1 風力機風速作用示意圖

風切變在垂直高度方向的風速分布[9]，可以輪轂處中心風速為參考風速，得出某高度風速表達式

式中：V(H)是高度H處的風速，Vh為高度h處的參考風速，本文h取輪轂處高度。α為粗糙度指數(shù)，與地形相關。具體數(shù)值見表1。

表1 不同地形下粗糙度指數(shù)

在式（1）中以葉尖風速作為基準，考慮葉片的方位角θ以及葉輪半徑R，可以將含H的表達式寫為如下形式

風力機塔架對氣流起到阻塞作用，改變風速的大小和方向，稱為塔影效應。塔影效應主要存在于葉輪下方以塔架中心線為軸線的兩側120o的范圍內(nèi)，在塔架正前方葉片位置影響最大。建模時增加一個速度影響矢量因子以修正風速模型

其中：q=1+α(α-1)R2(8h2)-1，c0為風機塔筒支架半徑，u0為分析截面微元距塔架中心軸線距離，稱懸垂距離。

因此在綜合考慮風切變和塔影效應之后，對葉片整體建立風速模型，以方位角和葉輪半徑作為自變量可得風速模型如下

以小型風力機參數(shù)為輸入?yún)?shù)，選取輪轂輸入風速為15 m/s，根據(jù)《風力等級》國家標準，此風速對應疾風，對小型風力機破壞較大。塔架半徑取為0.5 m，塔架高度取為6 m，風輪半徑取為2 m，粗糙度指數(shù)因子選取為0.2，懸垂距離取為0.5 m進行仿真，由風速模型經(jīng)仿真計算可得風速空間分布函數(shù)。將風速函數(shù)離散化為N位，取平均值作為葉片氣彈模型的輸入風速。

2 葉片結構模型及氣彈方程

近年來該領域?qū)W者多以大型風力機作為研究對象，而對小型風機研究相對較少。本文選取小型風機結構模型數(shù)值進行仿真分析?？紤]分析截面為葉片翼型典型對稱面，葉片氣動彈性中心與結構重心相距較近，因此彈性扭轉變形可忽略。分析翼型截面如圖3所示，其中U為輸入風速，此處選為經(jīng)風速模型處理后的平均風速；c為翼型截面弦長；v0為相對風速；α0為攻角；?為變槳角；ψ為相對風速角；翼型上作用的非線性氣動力包括非線性循環(huán)和非循環(huán)氣動升力TNc與TC1；DN1為非線性中的氣動阻力項。

圖2 翼型截面參數(shù)示意圖

如圖2所示的運動中，y方向表示葉片垂直于葉輪旋轉平面的揮舞運動；z0方向表示展長方向；z方向表示葉片在旋轉平面上垂直于揮舞的擺振運動；x方向表示葉片轉動。建模計算時綜合考慮系統(tǒng)動能和勢能，基于前期研究成果[10]，利用拉氏方程展開，忽略高階項，可以得到揮舞和擺振運動方向運動方程

根據(jù)前期研究成果，采用大攻角下修正后的ONERA氣動模型及非線性求解理論，利用傅氏變換提取氣動力中的諧波項，可以得到氣動力描述模型升力方程為

氣動阻力方程為

其中：KF1、KF2、KDF為升力和阻尼分別對應的循環(huán)結構變量，比例系數(shù)m=0.5，j=0.014，循環(huán)變量方程表達為

其中：v0是相對風速，與輸入風速和轉速相關。為翼型截面回轉半徑，b為翼型截面半弦長，ρ為空氣密度，SL、qL、ρL、LK分別為氣動項系數(shù)。ΔLKC為靜態(tài)氣動升力曲線延長線值與非線性升力曲線的差值。與之相對地，ΔLKD為阻力曲線差值。其余氣動參量，根據(jù)項目組前期研究成果，參考文獻[11]，取值如下

式中其余參數(shù)取值如下

3 氣彈方程求解仿真

上述氣彈方程的仿真需通過降階后化為系數(shù)矩陣狀態(tài)空間形式，將2階方程化為1階，使方程易于仿真。同時在輸出方程中取相應矩陣，通過觀察輸出值對應振動位移的振幅和頻率來判斷葉片相應方向的振動情形，從而判斷系統(tǒng)響應是否發(fā)生發(fā)散，顫振是否發(fā)生。本文取常值變槳角進行仿真運算，觀察風載作用下葉片揮舞-擺振方向的振動情形。

對氣彈方程組進行降階處理。定義Y=[XTX′T]T，將方程聯(lián)立 ，可以得到 ：可右乘輸入、控制矩陣，在時域內(nèi)表示為u(t)。其中

系統(tǒng)穩(wěn)定性的判定方法有許多種，諸如特征值判定法、伯德圖、內(nèi)奎斯特圖、李雅普諾夫判定法等。在實際應用中，應綜合分析系統(tǒng)響應狀況。本文所述葉片系統(tǒng)，應從葉片響應振幅大小、穩(wěn)定時間、靜差等多方面分析，且往往出現(xiàn)某一個方向穩(wěn)定而另外方向不穩(wěn)定的現(xiàn)象，所以應觀察各方向振動情形，且該種穩(wěn)定情形需符合物理和工作實際，否則穩(wěn)定性判別無現(xiàn)實意義。

使用4階5級龍格-庫塔方法進行計算仿真，模擬大攻角大風速條件下的系統(tǒng)響應。取恒葉尖速比系數(shù)為1.2，λ=LΩ/U，L為葉片長度，取L=0.9R，變槳角?=π/6，由以上給定條件可計算出攻角數(shù)值。由于風速和轉速的共同影響，攻角值或在大風速條件下表現(xiàn)出某一較大值，從而影響葉片氣彈振動。

回轉半徑取在葉根處，模擬危險截面振動。一旦葉片發(fā)生顫振，就會表現(xiàn)出振幅不自衰減的振動，不斷吸收能量維持發(fā)散振動。葉片發(fā)生顫振時極易破壞失效，且顫振的發(fā)生很難察覺，除已有報道的顫振發(fā)生前的“蜂鳴聲”外幾乎無征兆可循。所以顫振的抑制是風電技術發(fā)展的重要方面。模型結構參數(shù)選取匹配于小型風力機的結構數(shù)值。

在確定模型輸入風速的前提下，代入小型風力機結構數(shù)值模擬大攻角狀態(tài)下系統(tǒng)響應，通過4階5級龍格-庫塔法計算出風力機葉片振動位移。通過圖3可以明顯看出，在大攻角大風速條件下，揮舞和擺振兩個方向都出現(xiàn)了振幅不自衰減的高頻振動，也由此可判斷出葉片此條件下發(fā)生顫振，需加以抑制。

表2 結構參數(shù)取值

4 振動控制策略

模型預測控制（MPC）是上世紀80年代發(fā)展起來的一種計算機控制算法。該種控制策略一經(jīng)出現(xiàn)便運用在實際工作過程中，并且在工作實踐中不斷完善和發(fā)展。該種算法采用步長可給定的滾動優(yōu)化和反饋校正等控制策略。MPC的控制思路是：根據(jù)模型未來的預測行為，基于誤差校正和滾動優(yōu)化選擇對受控目標最合適的控制動作進行控制[12]。選取葉片氣彈狀態(tài)空間模型為受控模型，將該模型離散化處理可得

其中兩式式尾新增兩項為葉片氣彈模型誤差修正項。令ΔY(k)=Y(k)-Y(k-1)，

所以可以將式（14）化為增量形式

輸出可寫為增量形ΔY0(k+1)=Y0(k+1)-Y0(k)，所以式（15）可以寫為

圖3 未控制時揮舞、擺振方向振動情形

式中W1=010×10，W2=I10×10。定義新狀態(tài)變量y(k)=[ΔYT(k)YT(K)]T，所以式（16）可以寫為

式中各系數(shù)矩陣均可由式（16）推算出，此處不贅述?；谏鲜剿邢到y(tǒng)狀態(tài)空間模型，可以得到系統(tǒng)在離散化后的某時刻系統(tǒng)軌跡。按照引入的系統(tǒng)模型誤差以及預設目標值進行離散化的預測時域內(nèi)的滾動優(yōu)化和誤差控制，按照MPC控制策略的一般規(guī)則，選取預測時域長度為11，選取控制時域長度為10，采樣步長和采樣次數(shù)分別為0.5、100，選取無約束條件下輸入、輸出，基于系統(tǒng)狀態(tài)空間的模型預測控制策略運用MATLAB軟件進行仿真，可以得出葉片控制后振動位移。

通過圖4可以明顯看出，經(jīng)MPC策略控制后，揮舞和擺振方向振動最終穩(wěn)定，很好抑制了葉片顫振問題。且揮舞和擺振方向的穩(wěn)定位移減少為控制前的3%、5.83%，穩(wěn)定時間少，靜差可接受。

在實際工程中，有多種方法來實現(xiàn)該控制策略：運用貝加萊PLC可以實現(xiàn)從MATLAB/SIMULINK代碼模型直接轉化為PLC控制程序；運用三菱PLC可以實現(xiàn)從MATLAB代碼-C語言-PLC語言的轉換；運用西門子TIA技術可以直接實現(xiàn)從MATLAB語言到PLC程序的轉換[13]?；蛲ㄟ^對控制器編程實現(xiàn)MPC控制算法與傳感系統(tǒng)和執(zhí)行裝置的結合，組成閉環(huán)振動控制系統(tǒng)。執(zhí)行裝置可選擇變槳激勵器、蒙皮內(nèi)置作動器（如形狀記憶合金電加熱、壓電材料電流控制等）進行物理實現(xiàn)，限于篇幅不再贅述。

5 結語

本文基于實際應用意義提出了可以實踐于葉片氣彈穩(wěn)定性控制的控制策略和系統(tǒng)模型。經(jīng)過仿真對比，MPC控制策略可以很好抑制葉片的氣彈顫振：葉片兩個方向上的位移發(fā)散得到根本抑制；振動的振幅、頻率大大減小，且具有較短的穩(wěn)定時間和極小的靜差。并且該方法對受控模型精確度要求不高，具有控制效果好、魯棒性強等特點，在某些行業(yè)已取得了一定的發(fā)展，但在風力機葉片振動抑制方面應用較少，本文在葉片氣彈振動抑制技術和實踐控制理論相結合方面做出了一點探索。

圖4 經(jīng)MPC控制后揮舞、擺振方向振動情形