常 林，劉廷瑞

（山東科技大學(xué) 機(jī)械電子工程學(xué)院，山東 青島 266590）

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

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

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

1 風(fēng)速建模

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

圖1 風(fēng)力機(jī)風(fēng)速作用示意圖

風(fēng)切變?cè)诖怪备叨确较虻娘L(fēng)速分布[9]，可以輪轂處中心風(fēng)速為參考風(fēng)速，得出某高度風(fēng)速表達(dá)式

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

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

在式（1）中以葉尖風(fēng)速作為基準(zhǔn)，考慮葉片的方位角θ以及葉輪半徑R，可以將含H的表達(dá)式寫為如下形式

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

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

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

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

2 葉片結(jié)構(gòu)模型及氣彈方程

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

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

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

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

氣動(dòng)阻力方程為

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

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

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

3 氣彈方程求解仿真

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

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

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

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

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

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

表2 結(jié)構(gòu)參數(shù)取值

4 振動(dòng)控制策略

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

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

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

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

圖3 未控制時(shí)揮舞、擺振方向振動(dòng)情形

式中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)在離散化后的某時(shí)刻系統(tǒng)軌跡。按照引入的系統(tǒng)模型誤差以及預(yù)設(shè)目標(biāo)值進(jìn)行離散化的預(yù)測(cè)時(shí)域內(nèi)的滾動(dòng)優(yōu)化和誤差控制，按照MPC控制策略的一般規(guī)則，選取預(yù)測(cè)時(shí)域長(zhǎng)度為11，選取控制時(shí)域長(zhǎng)度為10，采樣步長(zhǎng)和采樣次數(shù)分別為0.5、100，選取無約束條件下輸入、輸出，基于系統(tǒng)狀態(tài)空間的模型預(yù)測(cè)控制策略運(yùn)用MATLAB軟件進(jìn)行仿真，可以得出葉片控制后振動(dòng)位移。

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

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

5 結(jié)語

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

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