曹九發(fā)， 柯世堂， 王同光

(南京航空航天大學(xué) 江蘇省風(fēng)力機(jī)設(shè)計高技術(shù)研究重點實驗室，南京　210016)

?

復(fù)雜工況下的大型風(fēng)力機(jī)氣動彈性響應(yīng)和尾跡數(shù)值分析研究

曹九發(fā)， 柯世堂， 王同光

(南京航空航天大學(xué) 江蘇省風(fēng)力機(jī)設(shè)計高技術(shù)研究重點實驗室，南京210016)

摘要：在復(fù)雜工況下，大型風(fēng)力機(jī)受到載荷更加嚴(yán)重，導(dǎo)致風(fēng)力機(jī)氣動和結(jié)構(gòu)耦合響應(yīng)問題更加明顯。主要針對穩(wěn)態(tài)偏航、動態(tài)偏航、風(fēng)剪切和隨機(jī)風(fēng)速場等復(fù)雜工況，采用非定常自由渦尾跡方法計算尾跡形狀和氣動載荷，加入了復(fù)雜工況的模型，進(jìn)行了動態(tài)失速模型和三維旋轉(zhuǎn)效應(yīng)模型修正。在考慮氣動載荷、慣性載荷和重力載荷影響下，采用有限元法結(jié)合模態(tài)法建立起風(fēng)力機(jī)解耦動力學(xué)方程，并且通過Newmark方法進(jìn)行數(shù)值求解該方程。實現(xiàn)了復(fù)雜工況數(shù)值模擬計算，比較不同復(fù)雜工況的氣動彈性響應(yīng)結(jié)果。最后，得出大型風(fēng)力機(jī)在復(fù)雜工況下的氣動性能、載荷、動態(tài)響應(yīng)和尾跡葉尖渦線特性，并計算出風(fēng)力機(jī)在復(fù)雜工況下的遲滯時間。這為推進(jìn)自由渦尾跡方法應(yīng)用于大批工況載荷計算，以及提高大型風(fēng)力機(jī)載荷計算精度和設(shè)計水平等具有重要意義。

關(guān)鍵詞：自由渦尾跡；風(fēng)剪切；動態(tài)偏航；隨機(jī)風(fēng)速場；動態(tài)響應(yīng)

隨著風(fēng)能技術(shù)的發(fā)展，風(fēng)力機(jī)逐漸向著大型化發(fā)展。由于風(fēng)力機(jī)尺寸變大，非定常因素的影響變得越來越顯著，比如風(fēng)剪切、隨機(jī)風(fēng)速場和動態(tài)偏航等。復(fù)雜工況的氣動性能、載荷和氣動彈性響應(yīng)計算研究對于風(fēng)力機(jī)設(shè)計初級階段具有重要意義。

目前，對于風(fēng)力機(jī)氣動分析方法，主要有三類：葉素動量理論BEM(Blade Element Momentum)[1-2]、渦尾跡方法[3-4]和CFD(Computational Fluid Dynamics)方法[5-6]。BEM方法計算時間快，但是它是基于靜態(tài)平衡尾跡假設(shè)，其中，對于動態(tài)工況，如動態(tài)偏航和動態(tài)變槳等，必須加入動態(tài)修正的經(jīng)驗?zāi)Ｐ瓦M(jìn)行相應(yīng)工況修正，從而，不同風(fēng)力機(jī)和不同動態(tài)工況，該方法會出現(xiàn)局限性問題，可能導(dǎo)致風(fēng)力機(jī)氣動特性和氣動載荷的計算結(jié)果不準(zhǔn)確；而CFD方法能模擬三維非穩(wěn)態(tài)黏性流，但影響因素多且黏性耗散、邊界層轉(zhuǎn)捩等的處理有待解決[7]，而且網(wǎng)格量多，計算效率低，對于工程應(yīng)用具有局限性；渦尾跡方法是介于BEM與CFD之間，渦流理論采用升力線或者升力面對葉片進(jìn)行模擬，尾跡處理則有預(yù)定渦尾跡[4,8-9]和自由渦尾跡[10-12]，目前的渦尾跡方法還主要研究穩(wěn)態(tài)風(fēng)速來流工況和一些簡單的非定常工況，對于復(fù)雜非定常工況研究欠缺，特別是具有風(fēng)剪切的隨機(jī)風(fēng)速場風(fēng)力機(jī)氣動特性研究。渦尾跡方法適合動態(tài)復(fù)雜工況模擬計算，并且充分考慮了葉片與尾跡之間的相互干擾，非常適合于風(fēng)力機(jī)氣動特性分析。同時，渦尾跡方法計算效率比CFD方法高，具有與BEM方法批量計算載荷的優(yōu)點，對于工程實際應(yīng)用具有重要意義。

風(fēng)力機(jī)氣動彈性問題的結(jié)構(gòu)分析方法分析過程中主要有三類模型被使用：基于工程梁理論或有限元方法的模型[13-14]；基于風(fēng)力機(jī)模態(tài)分析使用模態(tài)剛度和模態(tài)質(zhì)量的模型[15]；基于多體動力學(xué)方法的模型[16]。為了便于結(jié)構(gòu)動力學(xué)方程的解耦，本文采用有限元法結(jié)合模態(tài)法建立風(fēng)力機(jī)結(jié)構(gòu)模型。

本文首先，采用自由渦尾跡方法，進(jìn)行風(fēng)力機(jī)非定常工況的計算結(jié)果驗證分析，實現(xiàn)大型風(fēng)力機(jī)風(fēng)剪切、隨機(jī)風(fēng)速場和動態(tài)偏航數(shù)學(xué)模型，并且嵌入自由渦尾跡方法中；然后，考慮葉片速度和位移耦合，通過有限元法結(jié)合模態(tài)法解耦求解動力學(xué)方程，實現(xiàn)風(fēng)力機(jī)氣動和結(jié)構(gòu)的雙向耦合計算。最終，對穩(wěn)態(tài)風(fēng)速的動態(tài)偏航、風(fēng)剪切的動態(tài)偏航和隨機(jī)風(fēng)速場算例的氣動彈性響應(yīng)和尾跡進(jìn)行分析。

1自由渦尾跡方法

1.1尾流場描述

自由渦尾跡方法中，對風(fēng)力機(jī)流場作不可壓和位流假設(shè)，氣動模型可以簡化為來流、葉片附著渦線和自由渦面的總和，葉片附著渦線置于1/4弦線處，并采用“arc-cosine”法離散，每段附著渦線代替每段葉素，葉素控制點置于3/4弦線處，從而葉片被模擬成一個Weissinger-L升力面模型，如圖1所示。自由渦面是由葉片尾緣拖出渦線形成，可分為尾隨渦線和脫體渦線，分別模擬附著環(huán)量在空間和時間上的變化。尾隨渦強(qiáng)度定義為相鄰葉素的附著環(huán)量之差：

(Γt)i,j=

(1)

式中:j=1,2,…,NE，體現(xiàn)不同葉素。

圖1　風(fēng)力機(jī)尾跡離散描述示意圖Fig.1 Wake discrete description of wind turbine

相鄰方位角上葉素附著環(huán)量渦之差是脫體渦的強(qiáng)度，則第j個方位角下的脫體渦強(qiáng)度是：

(Γs)i,j=

(2)

式中:i=1,2,…,NT,體現(xiàn)不同方位角。

1.2渦線方程與求解

本文模型中每根渦線均在遠(yuǎn)場截斷，流場中渦線隨當(dāng)?shù)亓魉僖苿訒杂删砥?。渦線的偏微分控制方程可寫為：

(3)

推導(dǎo)尾跡控制方程中的對時間步微分方程的差分，令Δψ=Δζ，可得到控制方程離散格式[12]：

預(yù)估步：

(4)

校正步：

(5)

本文動態(tài)失速模型采用L-B模型，進(jìn)行風(fēng)力機(jī)的適當(dāng)動態(tài)失速模擬修正。根據(jù)動態(tài)失速L-B模型，可將二維翼型的非定常特性通過附著流、分離流和動態(tài)失速渦來模擬[3]。而對于三維旋轉(zhuǎn)效應(yīng)修正采用Du-Selig模型[3,17]。

2風(fēng)速場與動態(tài)偏航工況

2.1風(fēng)剪切模型

風(fēng)剪切就是指穩(wěn)態(tài)平均風(fēng)速隨高度的變化情況。在考慮風(fēng)剪切時，可選用以下模型。模型中V(h)是指高度h處的風(fēng)速，V(h0)是指參考高度h0處的風(fēng)速。常用風(fēng)剪切修正模型包括指數(shù)模型和對數(shù)模型，本文選用前者，其修正公式如下[2]：

(6)

當(dāng)不考慮風(fēng)剪切的影響時，可以將α的值設(shè)為0，取值范圍和地表情況相對應(yīng)。h0是輪轂的位置，V(h0)是輪轂的參考風(fēng)速。

2.2動態(tài)偏航模型

當(dāng)風(fēng)輪出現(xiàn)動態(tài)偏航或者動態(tài)變槳時，風(fēng)輪尾跡的變化也是一個非定常過程，從而導(dǎo)致風(fēng)輪的氣動特性也發(fā)生改變。

通過風(fēng)向動態(tài)的改變來模擬風(fēng)輪動態(tài)偏航效果。參照風(fēng)力機(jī)設(shè)計認(rèn)證標(biāo)準(zhǔn)GL2003，采用極端風(fēng)向變化工況，其中風(fēng)向隨時間的函數(shù)為：

γ(t)=

(7)

式中:γe=30°，為風(fēng)向變化幅值；T=6 s，為風(fēng)速變化周期。

2.3隨機(jī)風(fēng)速模型

由于風(fēng)力機(jī)運行環(huán)境中的風(fēng)剪切和湍流的存在，導(dǎo)致來流風(fēng)速發(fā)生隨機(jī)性脈動變化，從而影響風(fēng)輪氣動特性，特別是大型風(fēng)力機(jī)受其影響更明顯。為了研究大型風(fēng)力機(jī)在隨機(jī)風(fēng)速場激勵下的氣動特性。把風(fēng)輪風(fēng)速場離散化成36個空間點，基于改進(jìn)Von karman隨機(jī)風(fēng)速頻譜函數(shù)，采用小波逆變換方法[18]，仿真出相互獨立的風(fēng)力機(jī)隨機(jī)風(fēng)場，再對其進(jìn)行空間相關(guān)性修正，建立起符合實際風(fēng)力機(jī)運行風(fēng)速場特性的隨機(jī)三維風(fēng)速場模型。風(fēng)力機(jī)風(fēng)速場嵌入自由渦尾跡方法中的流場圖如圖2。

圖2　自由渦尾跡方法的風(fēng)速場示意圖Fig.2 Wind field of the wake vortex method

3結(jié)構(gòu)動力學(xué)方程建立與求解

3.1氣動結(jié)構(gòu)耦合方式

本文采用自由渦尾跡方法進(jìn)行氣動載荷計算。但是，由于涉及到氣動與結(jié)構(gòu)的耦合計算，因此，根據(jù)自由渦尾跡方法的特點，在采用該方法進(jìn)行氣動載荷計算中，會出現(xiàn)兩個耦合項：葉片振動速度耦合項和葉片變形位移耦合項。其中，葉片振動速度耦合項主要是葉片揮舞速度和擺振速度，它們影響著每個葉素的相對速度大小，從而，每個葉素的相對速度計算可以表示如下：

(8)

式中：vrel_x和vrel_z分別是葉片擺振和揮舞方向的相對速度，v0_x和v0_z分別是葉片擺振和揮舞方向的初始來流速度，vrot是風(fēng)輪旋轉(zhuǎn)速度，vind_x和vind_z是流場中所有渦線對每個葉素上的擺振和揮舞方向的誘導(dǎo)速度，vb_x和vb_z是葉片的擺振和揮舞速度。

對于葉片變形位移耦合項是：揮舞方向位移、擺振方向位移和徑向位移。在本文的自由渦尾跡方法計算中，由于尾跡的形狀直接影響到流場中誘導(dǎo)速度的大小，因此，必須考慮葉片變形產(chǎn)生的位移對尾跡和誘導(dǎo)速度的影響。通過該項耦合可 以研究風(fēng)力機(jī)在考慮葉片變形后，尾跡和流場速度的變化情況，這對于采用BEM方法來說是很難實現(xiàn)。并且對提高風(fēng)力機(jī)載荷和流場速度計算精度，具有重要意義。尾跡位移耦合表達(dá)式可以表示如下式：

(9)

3.2有限元模型

通過模態(tài)疊加法對結(jié)構(gòu)動力學(xué)進(jìn)行解耦求解。系統(tǒng)自由振動方程的廣義特征值問題為：

Kφ-ω2Mφ=0

(10)

從而可將系統(tǒng)的動力學(xué)方程寫為式(12)。風(fēng)力機(jī)簡化模型如圖3。通過梁單元進(jìn)行葉片的模擬，輪轂和塔架采用殼單元模擬，機(jī)艙看成0D質(zhì)量點。部件之間都采用剛性連接。

圖3　風(fēng)力機(jī)簡化模型圖Fig.3 The simplified model of Wind turbine diagram

本文以NH1500風(fēng)力機(jī)為算例，風(fēng)輪旋轉(zhuǎn)速度為17.2 r/min。采用商業(yè)軟件Patran/nastran，分別計算了風(fēng)輪和風(fēng)力機(jī)整機(jī)模態(tài)，相應(yīng)的固有頻率如表1所示，并且與商業(yè)軟件GH Bladed計算結(jié)果對比。從表格中可以看出， 頻率計算結(jié)果和Bladed軟件計算結(jié)果很接近，塔架模態(tài)頻率最低，風(fēng)輪存在揮舞和擺振模態(tài)，并且每一階揮舞和擺振都具有兩個反對稱和一個對稱振型，相比固有頻率而言，揮舞頻率較低些。通過只有風(fēng)輪模型和風(fēng)力機(jī)整機(jī)模型計算值對比，可以發(fā)現(xiàn)整機(jī)大部分模態(tài)比只有風(fēng)輪模型的固有頻率低，這是由于考慮輪轂、機(jī)艙和塔架后，使得整個風(fēng)力機(jī)模型剛度變小導(dǎo)致；并且在旋轉(zhuǎn)作用下，風(fēng)力機(jī)模態(tài)固有頻率都被提高，可見離心力的“剛化作用”分量大于“柔化作用”分量。

表1　風(fēng)力機(jī)固有頻率對比

圖4　考慮旋轉(zhuǎn)的整機(jī)前11階模態(tài)振型圖Fig.4 11 Order modes before of the wind turbine in rotational condition

3.3結(jié)構(gòu)動力學(xué)方程與求解

風(fēng)力機(jī)結(jié)構(gòu)系統(tǒng)結(jié)合模態(tài)動力方程可表示為：

(12)

式中:M是質(zhì)量矩陣，C是比例阻尼矩陣，K是剛度矩陣，F(xiàn)(t)是廣義力，x是廣義位移。其中，風(fēng)力機(jī)載荷包括重力，慣性力，氣動力。結(jié)構(gòu)阻尼本文都用0.005來進(jìn)行計算。對其微分方程的求解采用Newmark方法離散求解，具體形式如下：

(13)

整個氣動結(jié)構(gòu)耦合求解實現(xiàn)過程見圖5。

圖5　程序求解過程示意圖Fig.5 Procedure solving process diagram

4算例計算與結(jié)果分析

4.1計算程序驗證與分析

為了驗證本文非定常氣動模型的可靠性，通過美國可再生能源實驗室進(jìn)行的NREL Phase VI風(fēng)力機(jī)的非定?？諝鈩恿W(xué)系列實驗[19]，進(jìn)行驗證計算與分析。

圖6是穩(wěn)態(tài)風(fēng)速10m/s下偏航30°工況下的法向力系數(shù)Cn和切向力系數(shù)Ct。從圖中可以看出，“2D計算值”和“3D計算值”葉根的氣動載荷預(yù)測會比葉片中部和葉尖差，并且主要是體現(xiàn)在0°方位角附近(“2D計算值”是程序計算中沒有加三維旋轉(zhuǎn)效應(yīng)修正和動態(tài)失速模型修正；“3D計算值”是自由渦尾跡方法中嵌入了Du-Selig模型三維旋轉(zhuǎn)效應(yīng)修正模型和LB動態(tài)失速修正模型)。通過偏航時迎角的變化情況，可得知風(fēng)輪在0°方位角附近，葉片的迎角處于大迎角，并且變化幅度較大，從而存在迎角大、變化頻率快的動態(tài)失速修正困難的問題，從而導(dǎo)致“2D計算值”在葉根處的Cn和Ct都很難與實驗值匹配，然而，通過三維旋轉(zhuǎn)效應(yīng)和動態(tài)失速模型修正的“3D計算值”會使得數(shù)值計算值精度大大提高，但相對于葉片中部和葉尖位置的計算誤差仍然較大，這需要進(jìn)一步研究更好的修正模型來提高葉根的氣動載荷計算精度。而對于“2D計算值”和“3D計算值”的計算精度比較，不管是葉片的展向位置來看，還是方位角來看，“3D計算值”大部分都是要比“2D計算值”好很多。因此，本文氣動計算模塊具有一定的可靠性。

圖6　穩(wěn)態(tài)風(fēng)速下葉片展向的Cn和CtFig.6 Variation of normal and tangential coefficient at different span with wind speed

4.2大型風(fēng)力機(jī)算例結(jié)果與分析

NH1500風(fēng)力機(jī)是變槳變速風(fēng)力機(jī)，具體的風(fēng)力機(jī)性能和幾何參數(shù)見表2。針對該風(fēng)力機(jī)主要進(jìn)行兩個算例數(shù)值計算，第一個是考慮風(fēng)剪切的偏航和動態(tài)偏航對比算例；第二個是考慮風(fēng)剪切的三維隨機(jī)風(fēng)場工況算例。采用風(fēng)力機(jī)轉(zhuǎn)速為17.2 r/min。其中，風(fēng)速場模型的主要參數(shù)為：參考風(fēng)速為輪轂處風(fēng)速即10 m/s，粗糙度為0.01，采樣時間為0.1 s，采樣點數(shù)為8 192，風(fēng)剪切指數(shù)模型的系數(shù)0.2。

表2　NH1500風(fēng)力機(jī)參數(shù)

基于本文的風(fēng)力機(jī)整機(jī)氣動和結(jié)構(gòu)耦合計算算法，模擬計算的穩(wěn)態(tài)偏航和動態(tài)偏航工況都是仿真80 s，并且兩個工況都考慮風(fēng)剪切。

圖7是動態(tài)偏航在第40 s時，風(fēng)力機(jī)葉尖渦線形狀圖，從圖中可以看出，風(fēng)力機(jī)在動態(tài)偏航過程中，尾跡變化情況。其中，葉尖渦線出現(xiàn)了嚴(yán)重的不對稱性，隨著偏航的運動，葉尖渦線也開始緩慢地偏轉(zhuǎn)，并且YZ平面可以看出風(fēng)力機(jī)受到風(fēng)剪切的影響整個尾跡出現(xiàn)“上稀疏下密集”分布特點。

圖7　動態(tài)偏航40 s時的葉尖渦線圖Fig.7 Tip vortex lines at 40s in dynamic yaw condition

由于在揮舞方向上氣動和結(jié)構(gòu)的相互影響比較明顯，因此，給出了穩(wěn)態(tài)偏航和動態(tài)偏航的揮舞方向的葉尖位移和葉根載荷，如圖8和圖9所示。(在動態(tài)偏航工況前，增加20 s進(jìn)行收斂計算)。從圖中可以看出，大概20 s時各個工況都已經(jīng)計算收斂了，收斂之后穩(wěn)態(tài)偏航工況表現(xiàn)出比無偏航工況振蕩幅度大，并且出現(xiàn)響應(yīng)周期曲線不對稱的現(xiàn)象；動態(tài)偏航在25 s處開始發(fā)生葉尖位移和葉根載荷的降低，這是由于25 s開始發(fā)生偏航動作，并且其變化趨勢是：逐漸減小，接著緩慢回到平衡位置，最后與穩(wěn)態(tài)偏航的效應(yīng)曲線重合。由此可見，在動態(tài)偏航過程中，葉尖揮舞位移和葉根揮舞載荷響應(yīng)變化明顯，雖然總體載荷會出現(xiàn)減小，但會有高頻響應(yīng)出現(xiàn)，這對于風(fēng)力機(jī)疲勞具有重要影響。

圖10和圖11是風(fēng)輪的扭矩和推力系數(shù)變化曲線，給出無偏航、穩(wěn)態(tài)偏航、動態(tài)偏航和剛體風(fēng)力機(jī)動態(tài)偏航的數(shù)據(jù)。從圖中可以觀察出，偏航、風(fēng)剪切和風(fēng)力機(jī)結(jié)構(gòu)振動因素給風(fēng)輪動態(tài)響應(yīng)曲線帶來的不同變化振蕩幅度和非對稱周期響應(yīng)特點：無偏航風(fēng)剪切響應(yīng)曲線振蕩幅度最大，而有偏航風(fēng)剪切的振蕩幅度減弱了，但是出現(xiàn)鋸齒形狀的非對稱周期響應(yīng)曲線。

相對穩(wěn)態(tài)偏航，動態(tài)偏航的風(fēng)輪扭矩和推力都發(fā)生了遲滯現(xiàn)象，并且遲滯現(xiàn)象會比葉尖位移和葉根載荷明顯，柔性和剛性風(fēng)力機(jī)扭矩遲滯時間分別為13.7 s和17.9 s。通過對比剛體風(fēng)力機(jī)和柔性風(fēng)力機(jī)響應(yīng)曲線，出現(xiàn)柔性風(fēng)力機(jī)有剛體3倍的響應(yīng)周期，這是由于風(fēng)力機(jī)葉片重力和結(jié)構(gòu)振動導(dǎo)致了周期改變，另外使得扭矩和推力系數(shù)均值都變大了。

圖8 穩(wěn)態(tài)和動態(tài)偏航葉尖揮舞位移Fig.8Displacementofflapwisedirectionofbladetipinsteadyanddynamicyawcondition圖9 穩(wěn)態(tài)和動態(tài)偏航葉根揮舞方向剪力Fig.9Shearingforceofflapwisedirectionofbladerootinsteadyanddynamicyawcondition圖10 穩(wěn)態(tài)和動態(tài)偏航的風(fēng)輪扭矩曲線Fig.10Torqueofrotorinsteadyanddynamicyawcondition

本文隨機(jī)風(fēng)速場輸入工況仿真計算，該隨機(jī)風(fēng)速場包括三個方向：縱向風(fēng)速、橫向風(fēng)速和垂直向風(fēng)速，并且風(fēng)速場包含有風(fēng)剪切效應(yīng)。風(fēng)力機(jī)氣動和結(jié)構(gòu)全耦合模型進(jìn)行計算，其中風(fēng)力機(jī)模型為整機(jī)模型(全耦合是指速度和位移都耦合的雙向動態(tài)響應(yīng)計算，無耦合計算是指速度和位移不進(jìn)行耦合的單向動態(tài)響應(yīng)計算)。

圖12為風(fēng)力機(jī)葉尖位移的動態(tài)響應(yīng)時間曲線圖。從圖中可以看出，無耦合計算的葉尖位移波動幅值大于全耦合計算結(jié)果，因此，速度和位移耦合起到了氣動阻尼作用，減弱了振動幅度。圖13為風(fēng)輪的扭矩隨時間變化曲線。從圖中可以看出，在受到隨機(jī)風(fēng)速場的干擾后，風(fēng)輪扭矩也出現(xiàn)隨機(jī)變化趨勢，并且柔性風(fēng)力機(jī)扭矩和推力系數(shù)均值都比剛性風(fēng)力機(jī)大，趨勢與動態(tài)偏航分析的結(jié)果保持一致。由于扭矩和Cp曲線是對應(yīng)關(guān)系，同時可以得出風(fēng)力機(jī)氣動性能在變?nèi)嵝院笤谀承┕r可以提高風(fēng)能利用系數(shù)。

圖11 穩(wěn)態(tài)和動態(tài)偏航的風(fēng)輪推力系數(shù)曲線Fig.11Thrustcoefficientofrotorinsteadyanddynamicyaw圖12 隨機(jī)風(fēng)速場工況的葉尖揮舞位移Fig.12Displacementofflapwisedirectionofbladetipinstochasticwindspeedfieldcondition圖13 隨機(jī)風(fēng)速場工況的風(fēng)輪扭矩Fig.13Torqueofrotorinstochasticwindspeedfield

5結(jié)論

本文基于自由渦尾跡方法，針對大型風(fēng)力機(jī)在復(fù)雜工況下，完成了具有風(fēng)剪切效應(yīng)的穩(wěn)態(tài)偏航、動態(tài)偏航和隨機(jī)風(fēng)速場模塊的數(shù)值求解程序，并且其中加入了動態(tài)失速模型和三維旋轉(zhuǎn)效應(yīng)模型修正。同時，通過PhaseVI風(fēng)力機(jī)非定常實驗數(shù)據(jù)進(jìn)行驗證，最后，進(jìn)行了響應(yīng)的數(shù)值模擬計算以及氣動彈性響應(yīng)、載荷和尾跡葉尖渦線分析，得出以下結(jié)論：

(1) 對幾種風(fēng)力機(jī)簡化模型模態(tài)進(jìn)行了分析，總結(jié)了無旋轉(zhuǎn)風(fēng)輪、有旋轉(zhuǎn)風(fēng)輪模型、無旋轉(zhuǎn)整機(jī)模型和有旋轉(zhuǎn)整機(jī)模型的固有頻率特點，旋轉(zhuǎn)帶來的“剛化作用”分量大于“柔化作用”分量。

(2) 在動態(tài)偏航的工況尾跡分析中，動態(tài)偏航尾跡不僅體現(xiàn)出“上稀疏下密集”的分布風(fēng)剪切特點，還捕捉到動態(tài)偏航的尾跡變化過程。這些尾跡特點的捕捉對于保證風(fēng)力機(jī)載荷計算精度有重要意義。

(3) 對比了考慮風(fēng)剪切的無偏航、穩(wěn)態(tài)偏航和動態(tài)偏航工況的氣動彈性響應(yīng)特點，偏航使得響應(yīng)均值的減小，并且動態(tài)偏航在完成偏航動作的過程中高頻響應(yīng)多于穩(wěn)態(tài)偏航工況。從響應(yīng)周期和響應(yīng)幅度方面看，柔性風(fēng)力機(jī)偏航下的風(fēng)輪扭矩和推力都與剛性風(fēng)力機(jī)有著明顯的差別，并且柔性和剛性風(fēng)力機(jī)扭矩遲滯時間分別為13.7 s和17.9 s。因此，對于大型風(fēng)力機(jī)的柔性問題必須被考慮。

(4) 在隨機(jī)風(fēng)速場模擬計算中，葉片位移和速度的全耦合起到氣動阻尼的效果，減緩了葉片氣動彈性響應(yīng)，并且柔性風(fēng)力機(jī)的扭矩會比剛性風(fēng)力機(jī)更大些。

綜上所述，本文基于自由渦尾跡方法的復(fù)雜工況計算結(jié)果具有一定可靠性，并且隨機(jī)風(fēng)速場的嵌入，從計算時間和準(zhǔn)確性，使得自由渦尾跡方法應(yīng)用于工程的大批量工況載荷計算具有可行性。對于大型風(fēng)力機(jī)的氣動彈性響應(yīng)分析、載荷計算、設(shè)計以及優(yōu)化等具有重要意義。

參 考 文 獻(xiàn)

[1] Rijs R P P, Jacobs P, Smulders P T. Parameter study of the performance of slow running rotors[J]. Wind Energy and Industrial Aerodyn,1992, 39(1/2/3):95-103.

[2] Burton T,Sharpe D,Jenkins N,et al.Wind Energy Handbook[M].2005.

[3] Wang T G. Unsteady aerodynamic modeling of horizontal axis wind turbine performance[D]. Glasgow: University of Glasgow, 1999.

[4] 王芳，王同光. 基于渦尾跡方法的風(fēng)力機(jī)非定常氣動特性計算[J].太陽能學(xué)報,2009,30(9):1286-1291.

WANG Fang, WANG Tong-guang. Wind turbine unsteady aerodynamic performance prediction based on the vortex wake method[J]. ACTA Energlae Solaris Sinica, 2009, 30(9): 1286-1291.

[5] Chaviaropoulos P K, Hansen M O L. Investigating three-dimensional and rotational effects on wind turbine blades by means of a quasi-3D Navier-Stokes solver [J]. Journal of Fluids Engineering,2000, 122(2): 518-548.

[6] Sandersen B, Pijl S P, Koren B. Review of computational fluid dynamics for wind turbine wake aerodynamics [J]. Wind Energy, 2011, 14(7): 799-819.

[7] Vermeer L J，Sorensen J N,Crespo A．Wind turbine aerodynamics[J] ．Progress in Aerospace Sciences,2003(39)：467-510.

[8] Coton F N, Wang T. The prediction of horizontal axis wind turbine performance in yawed flow using an unsteady prescribed wake model [J]. Proceedings of the Institution of Mechanical Engineers, Part A, Journal of Power and Energy,1999, 213: 33-43.

[9] Wang T, Coton F N. An unsteady aerodynamic model for HAWT performance including tower shadow effects [J]. Wind Engineering,1999,23: 255-268.

[10] Sebastian T, Lackner M A. Development of a free vortex wake method code for offshore floating wind turbines [J]. Renewable Energy, 2012, 46: 269-275.

[11] Zhou W P, Tang S L, Lü H, Computation on aerodynamic performance of horizontal axis wind turbine based on time-marching free vortex method [J]. Chin.Soc.for Elec.Eng., 2011, 31(29):124-130.

[12] 許波峰.基于渦尾跡方法的風(fēng)力機(jī)氣動特性研究[D]. 南京:南京航空航天大學(xué)，2013.

[13] 柯世堂, 王同光, 趙林, 等. 風(fēng)力機(jī)風(fēng)振背景、共振響應(yīng)特性及耦合項分析[J]. 中國電機(jī)工程學(xué)報, 2013,33(26):101-108.

KE Shi-tang, WANG Tong-guang, ZHAO Lin, et al.Background, resonant components and coupled effect of wind-induced responses on wind turbine systems[J]. Proceeding of the CSEE, 2013,33(26): 101-108.

[14] 任勇生, 張明輝. 水平軸風(fēng)力機(jī)葉片的彎扭耦合氣彈穩(wěn)定性研究[J]. 振動與沖擊, 2010, 29(7): 196-200.

REN Yong-sheng, ZHANG Ming-hui. Aeroelastic stability study on coupled flutter for horizontal axis wind turbine blades[J]. Journal of Vibration and Shock, 2010, 29(7): 196-200.

[15] Hansen M O L. Aerodynamics of wind turbines[M].Second Edition,Earthscan,2008.

[16] Schiehlen W. Multibody system dynamics: Roots and perspectives[J]. Multibody System Dynamics, 1997, 1: 149-188.

[17] Du Z, Selig M S. A 3-D stall delay model for horizontal axis wind turbine performance prediction[R]. AIAA-98-0021, 1998.

[18] Kitagawa T, Nomura T. A wavelet-bassed method to generate artificial wind fluctuation data[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2003, 91:943-964.

[19] Hand M M, Simms D A, Fingersh L J, et al. Unsteady aerodynamics experiment Phases Ⅵ: Wind tunnel test configurations and available data campaigns[R]. Golden: National Renewable Energy Laboratory, NREL/TP-500-29955, 2001.

基金項目：國家973計劃項目(2014CB046200)大型風(fēng)力機(jī)的關(guān)鍵力學(xué)問題研究及設(shè)計實現(xiàn)；國家自然科學(xué)基金(51208254)復(fù)雜環(huán)境下超大型冷卻塔風(fēng)振機(jī)理與等效靜風(fēng)荷載研究；江蘇高校優(yōu)勢學(xué)科建設(shè)工程資助項目

收稿日期：2014-09-23修改稿收到日期：2014-12-03

通信作者王同光 男，教授，博士生導(dǎo)師，1962年生

中圖分類號:O357; TK89

文獻(xiàn)標(biāo)志碼：A

DOI:10.13465/j.cnki.jvs.2016.01.009

Numerical analysis for aero-elastic responses and wake of a large scale wind turbine under complicated conditions

CAO Jiu-fa, KE Shi-tang, WANG Tong-guang

(Jiangsu Provincial Key Laboratory of Hi-Tech Research for Wind Turbine Design, Nanjing University of Aeronautics & Astronautics, Nanjing 210016, China)

Abstract:Large scale wind turbines suffer serious unsteady loads under complicated conditions, it leads to their obvious aero-elastic coupled responses. For steady yaw, dynamic yaw, wind shear and stochastic wind field, the free vortex method was used to calculate their aerodynamic loads and wake shapes. The dynamic stall model and the three-dimension stall delay model were taken into account. At the same time, considering the aerodynamic load, inertial load and gravity load, the finite element method was combined with the modal method to build the decoupled dynamic equations of a wind turbine, these equations were solved numerically with Newmark method. The aero-elastic dynamic responses under different complicated conditions were compared. Finally, the aerodynamic performance, load, dynamic responses and tip vortex line characteristics of wind turbines were deduced under complicated conditions. The results were significant for applying the free vortex method in load calculations of wind turbines, and improving the load calculations accuracy and design levels of large scale wind turbines.

Key words:free wake method; wind shear; dynamic yaw; stochastic wind field; dynamic response

第一作者 曹九發(fā) 男，博士生，1986年生

郵箱: tgwang@nuaa.edu.cn