陸 洋， 周桂林

（南京航空航天大學(xué) 旋翼動(dòng)力學(xué)國(guó)防科技重點(diǎn)實(shí)驗(yàn)室，江蘇 南京210016）

0 引 言

大型水平軸風(fēng)力機(jī)系統(tǒng)是強(qiáng)非線性流－剛?cè)狁詈系闹芷跁r(shí)變多體系統(tǒng)，結(jié)構(gòu)和運(yùn)動(dòng)關(guān)系非常復(fù)雜，對(duì)風(fēng)力機(jī)葉片進(jìn)行動(dòng)力學(xué)建模時(shí)必須考慮葉片的幾何非線性以及非定常氣動(dòng)載荷等因素，理論推導(dǎo)和數(shù)值計(jì)算都比較困難。另一方面，隨著風(fēng)力機(jī)額定功率不斷增加，葉片尺寸越來越大，細(xì)長(zhǎng)葉片在交變載荷下的揮舞變形更大。為防止葉尖與塔架碰撞，出現(xiàn)了預(yù)彎葉片設(shè)計(jì)概念，增加了風(fēng)力機(jī)葉片氣彈建模的難度。

風(fēng)力機(jī)葉片氣彈建模包括葉片結(jié)構(gòu)動(dòng)力學(xué)建模和空氣動(dòng)力學(xué)建模。文獻(xiàn)［1］采用模態(tài)分析法進(jìn)行葉片結(jié)構(gòu)動(dòng)力學(xué)建模，采用Greenberg's二維準(zhǔn)定常氣動(dòng)力模型計(jì)算氣動(dòng)力，研究了風(fēng)力機(jī)氣彈穩(wěn)定性；文獻(xiàn)［2］采用5節(jié)點(diǎn)18自由度梁?jiǎn)卧獙?duì)葉片進(jìn)行有限元離散，采用二維準(zhǔn)定常氣動(dòng)模型，利用擬線性法計(jì)算了水平軸風(fēng)機(jī)工作狀態(tài)下的氣動(dòng)彈性響應(yīng)；文獻(xiàn)［3］以Marc有限元軟件為基礎(chǔ)平臺(tái)，采用動(dòng)量葉素理論二維氣動(dòng)力模型，計(jì)算葉片的氣彈響應(yīng)；文獻(xiàn)［4］采用二節(jié)點(diǎn)梁?jiǎn)卧獙?duì)葉片進(jìn)行有限元離散，采用Leishman－Beddoes非定常模型計(jì)算氣動(dòng)力，分析了葉片的動(dòng)態(tài)響應(yīng)。

可見，為了保證葉片響應(yīng)的計(jì)算精度，有限元方法已逐漸替代模態(tài)法成為風(fēng)力機(jī)葉片結(jié)構(gòu)動(dòng)力分析的主流方法，氣動(dòng)力建模則仍以動(dòng)量葉素理論和二維準(zhǔn)定常氣動(dòng)力模型為主。

然而，現(xiàn)有的葉片氣彈建模方法仍存在以下局限性：1）基于有限元方法的風(fēng)力機(jī)動(dòng)力學(xué)方程通常建立在非慣性系下，即處理動(dòng)能部分時(shí)將其分為旋轉(zhuǎn)坐標(biāo)系的動(dòng)能和牽連運(yùn)動(dòng)引起的動(dòng)能兩部分，破壞了推導(dǎo)和編程的一致性；2）建模時(shí)，沒有考慮葉片預(yù)彎對(duì)葉片氣動(dòng)彈性動(dòng)力響應(yīng)的影響；3）氣動(dòng)力模型多采用二維準(zhǔn)定常空氣動(dòng)力計(jì)算模型，沒有考慮葉片動(dòng)態(tài)入流以及動(dòng)態(tài)失速影響。這些局限性，使得傳統(tǒng)建模方法的計(jì)算精度和適用性受到了制約。

針對(duì)前述葉片建模方法的局限性，本文開展了如下工作：1）引入了Hartenberg－Denavit增廣轉(zhuǎn)換矩陣［5］，通過遞歸計(jì)算的方式得到慣性坐標(biāo)系下的位置矢量，從而在慣性系下建立了風(fēng)力機(jī)葉片動(dòng)力學(xué)方程，將動(dòng)能部分合為一體，避免了傳統(tǒng)建模中的繁瑣；2）引入上反坐標(biāo)系方法描述風(fēng)力機(jī)葉片的預(yù)彎變形，考慮了預(yù)彎變形影響；3）入流模型采用Suzuki提出的針對(duì)風(fēng)機(jī)葉片的廣義動(dòng)態(tài)尾跡入流模型（GDW）［6］，并采用Pierce修正的風(fēng)機(jī)Leishman－Beddoes模型［7］計(jì)算翼型非定常氣動(dòng)力，綜合考慮了葉片動(dòng)態(tài)入流以及動(dòng)態(tài)失速影響。利用Newmark數(shù)值積分方法求解葉片動(dòng)力學(xué)方程獲得穩(wěn)態(tài)周期解。最后，分別以美國(guó)可再生能源實(shí)驗(yàn)室（NREL）Phase VI非定常空氣動(dòng)力學(xué)實(shí)驗(yàn)風(fēng)力機(jī)［8－9］及其公開的1.5MW風(fēng)力機(jī)葉片［10］為算例計(jì)算了有／無預(yù)彎葉片的氣彈響應(yīng)，驗(yàn)證了本文所建模型的正確性和有效性。

1 系統(tǒng)坐標(biāo)系

1.1 增廣轉(zhuǎn)換矩陣

為描述風(fēng)力機(jī)葉片復(fù)雜的結(jié)構(gòu)形式及運(yùn)動(dòng)特點(diǎn)，引入多個(gè)坐標(biāo)系，表1給出了建立風(fēng)力機(jī)系統(tǒng)的坐標(biāo)系定義和描述。圖1為風(fēng)力機(jī)系統(tǒng)坐標(biāo)系示意圖。對(duì)于開環(huán)拓?fù)浣Y(jié)構(gòu)，如果要描述梁剖面上的任意一點(diǎn)P，在變形后的坐標(biāo)系下，可以將位置矢量表示為：

將位置坐標(biāo)行向量S增廣為：

在未變形坐標(biāo)下，增廣后P點(diǎn)的位置增廣向量及增廣轉(zhuǎn)換矩陣為：

根據(jù)式（3）類推，在n個(gè)坐標(biāo)系中第i個(gè)坐標(biāo)系下，P點(diǎn)的位置增廣向量及增廣轉(zhuǎn)換矩陣為：

依次類推，在慣性坐標(biāo)系下，葉片上任意點(diǎn)P的位置增廣向量表示為：

由上面的遞歸表達(dá)式（4）可以看出，只要確定了各坐標(biāo)系原點(diǎn)的位置和轉(zhuǎn)換關(guān)系，就可以通過遞歸計(jì)算的方式得到慣性坐標(biāo)系下的位置矢量，從而簡(jiǎn)化了剛?cè)狁詈嫌?jì)算，使推導(dǎo)和編程也得到簡(jiǎn)化，提高了算法對(duì)風(fēng)力機(jī)結(jié)構(gòu)和運(yùn)動(dòng)形式的適用性。

圖1 風(fēng)力機(jī)系統(tǒng)坐標(biāo)系示意圖Fig.1 Schematic of coordinate system of wind turbine

表1 風(fēng)力機(jī)坐標(biāo)系統(tǒng)Table 1 Coordinate system of wind turbine

1.2 上反坐標(biāo)系

葉片設(shè)計(jì)時(shí)預(yù)先向偏離塔架方向彎曲，可以增大葉尖與塔架之間的凈空間隙，減少二者發(fā)生碰撞的可能性。大型風(fēng)力機(jī)葉片長(zhǎng)度通常超過30m，葉片一般從1／3長(zhǎng)度位置處開始預(yù)彎。對(duì)于長(zhǎng)度為35m左右的1.5MW風(fēng)機(jī)葉片，預(yù)彎高度一般取值為600～750mm。由于預(yù)彎值遠(yuǎn)小于葉片長(zhǎng)度，為了描述預(yù)彎部分葉片的位置矢量，本文提出一種簡(jiǎn)化處理辦法，如圖2所示，用直線葉片BC′來替代連續(xù)預(yù)彎葉片部分BB′C′，直線葉片BC′的上反角度為ζ0。因此，為了描述整個(gè)葉片的位置矢量，需要將葉片劃分為AB和BC′兩段進(jìn)行處理，在BC′段引入上反坐標(biāo)系Rs：osxsyszs。

圖2 預(yù)彎葉片示意圖Fig.2 Sketch of pre－curved rotor blade

上反坐標(biāo)系與變距坐標(biāo)系之間的轉(zhuǎn)換關(guān)系Ts為：

對(duì)于預(yù)彎量較大的葉片，可以采用類似方法，引入多個(gè)上反坐標(biāo)系。引入的上反坐標(biāo)系越多，則預(yù)彎變形情況描述越精確，綜合考慮計(jì)算效率，可選擇適當(dāng)數(shù)量的上反坐標(biāo)系進(jìn)行計(jì)算。

2 葉片動(dòng)力學(xué)建模

2.1 葉片單元彈性勢(shì)能貢獻(xiàn)項(xiàng)

采用中等變形梁理論得到Green應(yīng)變表達(dá)式［11］，帶入單元彈性勢(shì)能表達(dá)式并對(duì)其求變分，得到彈性廣義力及彈性勢(shì)能貢獻(xiàn)項(xiàng)：

2.2 葉片單元?jiǎng)幽茇暙I(xiàn)項(xiàng)

單片葉片上的動(dòng)能變分［12］可以表示為：

經(jīng)推導(dǎo)得動(dòng)能產(chǎn)生的第i個(gè)廣義力為：

動(dòng)能所產(chǎn)生的切線質(zhì)量、阻尼和剛度陣為：

為了計(jì)算質(zhì)量、阻尼和剛度陣，還需求出增廣位置矢量的相關(guān)導(dǎo)數(shù)和偏導(dǎo)數(shù)，利用遞推算法后，可使推導(dǎo)通用性大大提高，避免了傳統(tǒng)建模中的繁瑣。

2.3 重力功貢獻(xiàn)項(xiàng)

在計(jì)算風(fēng)力機(jī)風(fēng)輪葉片氣彈響應(yīng)時(shí)，必須考慮葉片的重力影響，葉片重力引起的外力功和重力廣義力分別為：

2.4 氣動(dòng)力貢獻(xiàn)項(xiàng)

誘導(dǎo)速度計(jì)算采用Suzuki廣義動(dòng)態(tài)尾跡模型［6］。用任意諧波次數(shù)和任意階次徑向型函數(shù)的級(jí)數(shù)形式來描述風(fēng)輪平面垂向誘導(dǎo)速度：

其中，Δ為環(huán)量項(xiàng)，Δ為非環(huán)量項(xiàng)。利用Kirchhoff理論處理翼型后緣分離的非線性問題，翼型法向力系數(shù)和弦向力系數(shù)與分離點(diǎn)的關(guān)系為：

式中、ψ分別是無量綱半徑、方位角和無量綱時(shí)間。是徑向函數(shù)，是基于時(shí)間的狀態(tài)量。

翼型剖面非定常氣動(dòng)力計(jì)算采用Pierce修正后應(yīng)用于風(fēng)力機(jī)計(jì)算的Leishman－Beddoes模型［7］。該模型通過階躍響應(yīng)的疊加模擬附著流條件下的非定常效應(yīng)：

2.5 葉片動(dòng)力學(xué)方程

應(yīng)用Hamilton原理建立系統(tǒng)運(yùn)動(dòng)方程：

其中，CN是法向力系數(shù)，CC是弦向力系數(shù)，f是分離點(diǎn)位置。根據(jù)動(dòng)態(tài)失速的特性，模擬動(dòng)態(tài)失速渦沿翼型上表面的移動(dòng)及因而引起的壓力中心的移動(dòng)，從而計(jì)算動(dòng)態(tài)失速渦誘導(dǎo)產(chǎn)生的升力和俯仰力矩［9］?？諝鈩?dòng)力FA及力矩MA所做虛功的變分可以表示為：

對(duì)式（16）進(jìn)行變分運(yùn)算和有限元離散，組集動(dòng)能、應(yīng)變能及空氣動(dòng)力等產(chǎn)生的虛功，得到基于廣義力形式的葉片非線性動(dòng)力學(xué)隱式微分方程：

式中右上標(biāo)T、E、A、G分別表示動(dòng)能、應(yīng)變能、氣動(dòng)力及重力引起的廣義力，q、、為廣義位移、速度、加速度向量。

由于建模過程中考慮了葉片整體的剛性運(yùn)動(dòng)與葉片自身彈性運(yùn)動(dòng)間的剛?cè)狁詈咸攸c(diǎn)及非定常和動(dòng)態(tài)失速等非線性因素，最終采用適應(yīng)性較強(qiáng)的Newmark數(shù)值積分方法進(jìn)行運(yùn)動(dòng)方程求解。

3 算例驗(yàn)證

算例1：對(duì)NREL Phase VI非定?？諝鈩?dòng)力學(xué)實(shí)驗(yàn)中葉片長(zhǎng)度為10.058m的小型水平軸風(fēng)力機(jī)進(jìn)行計(jì)算并與Bladed軟件以及實(shí)驗(yàn)結(jié)果比較。本文分別計(jì)算了來流速度分別為 7m／s、10m／s、15m／s、20m／s、25m／s，偏航角為0°時(shí)，下風(fēng)向風(fēng)機(jī)徑向位置為63%截面的法向力系數(shù)Cn，切向力系數(shù)Ct，俯仰力矩系數(shù)Cm以及葉片根部的揮舞力矩。圖3～圖6為本文計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)及Bladed計(jì)算結(jié)果的對(duì)比曲線。

從圖3～圖5可以看出，相比較Bladed計(jì)算結(jié)果，本文計(jì)算得到的截面氣動(dòng)力系數(shù)與實(shí)驗(yàn)結(jié)果吻合程度更佳，表明了本文葉片氣彈建模中空氣動(dòng)力學(xué)建模方法的正確性和精確性。圖6中計(jì)算結(jié)果與實(shí)驗(yàn)值的葉片根部揮舞彎矩吻合程度較高，而葉片根部揮舞力矩主要由氣動(dòng)力貢獻(xiàn)，再次證明了本文氣動(dòng)力模型的正確性。

算例2：以NREL 1.5MW下風(fēng)向風(fēng)機(jī)葉片為算例進(jìn)行氣彈響應(yīng)計(jì)算。該風(fēng)機(jī)基本參數(shù)如下：風(fēng)輪仰角η＝0°，風(fēng)輪錐角χ＝0°，葉片長(zhǎng)度33.25m，輪轂半徑1.75m，來流速度10m／s，風(fēng)輪轉(zhuǎn)速20rpm。葉片無預(yù)彎，剖面特性具體參數(shù)見文獻(xiàn)［8］。為了分析葉片預(yù)彎對(duì)葉片氣彈響應(yīng)的影響，假設(shè)該葉片在1／3位置處開始連續(xù)預(yù)彎，至葉尖處預(yù)彎高度為3.5m，即上反角10°，計(jì)算該預(yù)彎葉片的氣彈響應(yīng)。將上述計(jì)算結(jié)果與NREL公開文獻(xiàn)的計(jì)算結(jié)果一并進(jìn)行對(duì)比分析，以驗(yàn)證本文算法的正確性。圖7～圖12給出了葉片旋轉(zhuǎn)系下單片葉片作用于輪轂三個(gè)方向的載荷及力矩沿方位角變化對(duì)比曲線。

圖7 葉片根部徑向載荷Fig.7 Blade axial force at blade root

對(duì)比無預(yù)彎情況下葉片在旋轉(zhuǎn)系下的載荷和力矩，本文計(jì)算結(jié)果與NREL公開文獻(xiàn)的計(jì)算結(jié)果的變化趨勢(shì)以及幅值符合較好，其中葉片徑向載荷的相位存在約50°的差異。由于葉片徑向載荷主要由風(fēng)輪旋轉(zhuǎn)產(chǎn)生的離心力和周期變化的重力載荷分量構(gòu)成，考慮重力載荷在葉片徑向周期變化，可知本文計(jì)算結(jié)果更加合理。通過無預(yù)彎情況的算例對(duì)比，驗(yàn)證了本文建模算法的正確性。

圖12 葉片根部擺振力矩Fig12 Blade in－plane moment at blade root

對(duì)比有／無預(yù)彎兩種情況下基于本文模型計(jì)算所得的葉片載荷和力矩，可以看出對(duì)于下風(fēng)向風(fēng)機(jī)預(yù)彎對(duì)葉片的徑向載荷以及擺振方向的載荷影響不大（圖7、圖8），這是因?yàn)槿~片擺振方向的載荷主要是周期變化的重力載荷，預(yù)彎對(duì)重力載荷基本沒有影響。同理，由于葉片擺振力矩主要為擺振載荷貢獻(xiàn)，預(yù)彎對(duì)擺振力矩影響也很?。▓D12）。由圖9和圖11可知，預(yù)彎使得葉片揮舞方向載荷、力矩的穩(wěn)態(tài)值變小。揮舞方向載荷的穩(wěn)態(tài)值主要為氣動(dòng)力，表明預(yù)彎降低了葉片氣動(dòng)力載荷。由圖10葉片變距力矩對(duì)比曲線可以看出，預(yù)彎使得葉片變距力矩的絕對(duì)值增大。影響葉片變距力矩的因素很多，其中最重要的是由葉片自身質(zhì)量引起的離心力作用而產(chǎn)生的慣性力矩和變距結(jié)構(gòu)各種摩擦副產(chǎn)生的摩擦阻力距。下風(fēng)向風(fēng)機(jī)葉片由于預(yù)彎增加了葉片離心力作用產(chǎn)生的慣性力矩，從而使葉片變距力矩顯著增加。通過上述分析可知，本文采用上反坐標(biāo)系處理葉片預(yù)彎變形的方法所得到的預(yù)彎影響與理論分析相符，表明該處理方法是合理的。

4 結(jié) 論

本文基于Hamilton原理，建立了水平軸風(fēng)力機(jī)葉片的氣彈模型。通過數(shù)學(xué)建模和算例分析，可以得到以下結(jié)論：1）引入增廣轉(zhuǎn)換矩陣，在慣性系下建立葉片動(dòng)力學(xué)方程，可以簡(jiǎn)化剛?cè)狁詈嫌?jì)算，避免傳統(tǒng)建模中的繁瑣；2）氣動(dòng)力建模時(shí)，采用Suzuki廣義動(dòng)態(tài)尾跡方法計(jì)算誘導(dǎo)速度，采用Pierce修正的風(fēng)機(jī)Leishman－Beddoes模型計(jì)算翼型氣動(dòng)力，可提高了氣動(dòng)力計(jì)算精度，比Bladed計(jì)算結(jié)果更符合實(shí)驗(yàn)值；3）本文所建模型預(yù)測(cè)的無預(yù)彎葉片的載荷與NREL計(jì)算結(jié)果符合較好；4）采用上反坐標(biāo)系處理葉片預(yù)彎變形，計(jì)算所得的預(yù)彎影響與理論分析相符，驗(yàn)證了本文所建立的可考慮葉片預(yù)彎的氣彈動(dòng)力學(xué)模型的合理性。

［1］李本立，安玉華.風(fēng)力機(jī)氣動(dòng)彈性穩(wěn)定性研究［J］.太陽(yáng)能學(xué)報(bào)，1996，17（4）：314－320.

［2］王介龍，陳彥，薛克宗.風(fēng)力發(fā)電機(jī)耦合轉(zhuǎn)子／機(jī)艙／塔架的氣彈響應(yīng)［J］.清華大學(xué)學(xué)報(bào)（自然科學(xué)版），2002，42（2）：211－215.

［3］傅程，王延榮.風(fēng)力發(fā)電機(jī)葉片氣動(dòng)彈性響應(yīng)分析［J］.機(jī)械設(shè)計(jì)與研究，2009，25（1）：68－71.

［4］LIU X，ZHANG X M.Dynamic response analysis of the rotating blade of horizontal axis wind turbine［J］.WindEngineering，2010，34（5）：543－559.

［5］GERADIN M，CARDONA A.Kinematics and dynamics of rigid and flexible mechanisms using finite elements and quaternion algebra［J］.ComputationalMechanics，1989，18（4）：651－659.

［6］SUZUKI.Application of dynamic inflow theory to wind turbine rotors［D］.［Doctoral Dissertation］.Salt Lake City：Department of Mechanical Engineering，University of Utah，2000.

［7］KIRK GEE PIERCE.Wind turbine load prediction using the Beddoes－Leishman model for unsteady aerodynamics and dynamic stall［D］.［Doctoral Dissertation］.Salt Lake City：Department of Mechanical Engineering，University of Utah，1996.

［8］HAND M M，SIMMS D A.Unsteady aerodynamics experiment phase VI：wind tunnel test configurations and available data campaigns［R］.NREL／TP－500－29955，2001.

［9］HAND M M，SIMMS D A.NREL unsteady aerodynamics experiment in the NASA－Ames Wind Tunnel：A comparison of predictions to measurements［R］.NREL／TP－500－29494，2001.

［10］JASON M.FAST user's guide［R］.Technical Report NREL／EL－500－38230，2005.

［11］王浩文，高正，鄭兆昌.前飛狀態(tài)下直升機(jī)旋翼系統(tǒng)氣彈響應(yīng)及穩(wěn)定性分析［J］.振動(dòng)工程學(xué)報(bào)，1999，12（4）：521－528.

［12］王益鋒.直升機(jī)旋翼槳葉動(dòng)力學(xué)建模與彈性碰撞動(dòng)響應(yīng)分析［D］.［博士學(xué)位論文］.南京：南京航空航天大學(xué)，2009.