基于HHT方法的非穩(wěn)定工況風(fēng)力機(jī)結(jié)構(gòu)動(dòng)態(tài)響應(yīng)時(shí)頻特性分析

楊 陽， 李 春,2， 葉柯華， 繆維跑， 陽 君,2， 高 偉

(1. 上海理工大學(xué) 能源與動(dòng)力工程學(xué)院，上海 200093；2. 上海市動(dòng)力工程多相流動(dòng)與傳熱重點(diǎn)實(shí)驗(yàn)室，上海 200093； 3.大唐華創(chuàng)風(fēng)能有限公司，山東 青島 266000)

?

基于HHT方法的非穩(wěn)定工況風(fēng)力機(jī)結(jié)構(gòu)動(dòng)態(tài)響應(yīng)時(shí)頻特性分析

楊 陽1， 李 春1,2， 葉柯華1， 繆維跑1， 陽 君1,2， 高 偉3

(1. 上海理工大學(xué) 能源與動(dòng)力工程學(xué)院，上海 200093；2. 上海市動(dòng)力工程多相流動(dòng)與傳熱重點(diǎn)實(shí)驗(yàn)室，上海 200093； 3.大唐華創(chuàng)風(fēng)能有限公司，山東 青島 266000)

風(fēng)速變化劇烈的湍流風(fēng)場(chǎng)、開機(jī)啟動(dòng)、偏航以及緊急停機(jī)等典型非穩(wěn)定運(yùn)行工況均會(huì)增強(qiáng)風(fēng)力機(jī)非線性氣動(dòng)彈性響應(yīng)，時(shí)域和頻域的結(jié)構(gòu)動(dòng)力學(xué)響應(yīng)具有十分明顯的非平穩(wěn)特征。為此，基于湍流風(fēng)譜和相干結(jié)構(gòu)，建立了速度和方向均劇烈波動(dòng)的湍流風(fēng)，在氣動(dòng)-伺服-彈性仿真軟件FAST中計(jì)算了風(fēng)力機(jī)非穩(wěn)定工況下的動(dòng)力學(xué)特性，并與GH Bladed計(jì)算結(jié)果對(duì)比，驗(yàn)證了結(jié)果的有效性。使用HHT方法分析了塔架和葉片位移的時(shí)頻特性，結(jié)果表明：開機(jī)啟動(dòng)階段塔架和葉片位移均小幅振蕩約40 s后急劇增加，緊急停機(jī)均劇烈振蕩約20 s后恢復(fù)平穩(wěn)，偏航導(dǎo)致塔尖側(cè)向位移明顯上升。塔架位移響應(yīng)頻率主要集中于一階振動(dòng)頻率，偏航時(shí)幅值增大明顯。風(fēng)輪旋轉(zhuǎn)頻率為葉尖擺振的主要諧振動(dòng)頻率，葉片一階擺振頻率受到相干結(jié)構(gòu)影響，緊急停機(jī)時(shí)由于負(fù)氣動(dòng)阻尼影響而使得幅值增大，葉片設(shè)計(jì)時(shí)應(yīng)適當(dāng)增大阻尼以減小氣動(dòng)阻尼迅速降低帶來的振幅急速增加現(xiàn)象。

風(fēng)力機(jī)；湍流風(fēng)場(chǎng)；結(jié)構(gòu)動(dòng)力學(xué)；HHT方法

由于化石能源危機(jī)和環(huán)境污染日益嚴(yán)重，風(fēng)能因其資源廣泛、技術(shù)成熟和經(jīng)濟(jì)效益高等特點(diǎn)逐漸受到國家重視，2014年世界新增裝機(jī)容量達(dá)5 147萬kW，其中我國占45.1%[1]。近年來我國風(fēng)電場(chǎng)數(shù)量急劇增加，其中大部分選址于風(fēng)能資源相對(duì)豐富但自然環(huán)境十分復(fù)雜的“三北”和沿海地區(qū)。其風(fēng)速特點(diǎn)具有典型的時(shí)域非定常性和空間不均勻性，運(yùn)行于此類湍流風(fēng)的大型風(fēng)力機(jī)受到非穩(wěn)態(tài)氣動(dòng)載荷作用，氣動(dòng)彈性響應(yīng)和動(dòng)力學(xué)特性呈現(xiàn)強(qiáng)烈的非線性特點(diǎn)，在開機(jī)啟動(dòng)、緊急停機(jī)和偏航等非穩(wěn)定工況時(shí)愈加明顯，直接影響風(fēng)力機(jī)結(jié)構(gòu)安全和并網(wǎng)穩(wěn)定性。研究風(fēng)力機(jī)在此類典型非穩(wěn)定工況的動(dòng)力學(xué)特性已成為亟待解決的重要課題。

選擇一種合適的方法建立隨時(shí)間及空間變化的湍流風(fēng)是仿真首先需要解決的問題。較為常見的方法有：①基于經(jīng)典湍流譜模型如Kaimal譜或Von Karman譜，考慮空間相干性建立三維時(shí)變的風(fēng)場(chǎng)模型，該方法理論清晰容易實(shí)現(xiàn)，因此應(yīng)用最為廣泛[2-3]；②基于測(cè)風(fēng)塔采集的實(shí)際風(fēng)速數(shù)據(jù)，通過自回歸滑動(dòng)平均模型[4]、神經(jīng)網(wǎng)絡(luò)近似[5]或模糊邏輯預(yù)測(cè)[6]等方法得到小空間范圍的風(fēng)速分布規(guī)律，該方法較為準(zhǔn)確且適用性強(qiáng)，但對(duì)于大空間風(fēng)場(chǎng)模型需要極度豐富的實(shí)測(cè)數(shù)據(jù)作為近似樣本，成本較高；③考慮風(fēng)電場(chǎng)周圍實(shí)際地形和地表粗糙度等地貌條件，基于風(fēng)電場(chǎng)近期測(cè)量數(shù)據(jù)，通過氣象分析方法預(yù)測(cè)整個(gè)風(fēng)電場(chǎng)風(fēng)速分布[7]，該模型準(zhǔn)確度和可靠性高，但空間尺度過大容易導(dǎo)致氣動(dòng)載荷計(jì)算誤差偏大；④考慮大氣邊界層、地表粗糙度和科氏力等條件，通過大渦模擬的數(shù)值方法得到選定風(fēng)場(chǎng)域的三維風(fēng)場(chǎng)[8]，其空間及時(shí)間尺度均足以滿足仿真要求，但該方法需耗費(fèi)大量的計(jì)算資源。

除此之外需要注意的是，風(fēng)力機(jī)不僅受到非穩(wěn)定的風(fēng)載荷，開機(jī)啟動(dòng)、緊急停機(jī)或偏航等導(dǎo)致的突變載荷亦為瞬時(shí)響應(yīng)主要載荷[9-10]。風(fēng)力機(jī)運(yùn)行過程中受到非線性非平穩(wěn)的隨機(jī)載荷作用，其結(jié)構(gòu)動(dòng)力學(xué)時(shí)域特性是隨機(jī)載荷與結(jié)構(gòu)耦合作用的時(shí)程響應(yīng)結(jié)果，對(duì)時(shí)域結(jié)果進(jìn)行分析具有一定的參考意義，但根本在于頻率響應(yīng)。而簡(jiǎn)單的FFT方法的頻率響應(yīng)是基于全局時(shí)域特征，不能區(qū)分引起非平穩(wěn)響應(yīng)的真正頻率。盡管小波變換可同時(shí)從時(shí)域和頻域?qū)憫?yīng)信號(hào)進(jìn)行多尺度聯(lián)合分析，其局部細(xì)化分析能力可有效區(qū)分突變信號(hào)頻率和非平穩(wěn)白噪聲頻率，但結(jié)果準(zhǔn)確性嚴(yán)重依賴于小波基函數(shù)的選擇，如何針對(duì)特定問題準(zhǔn)確的選擇基函數(shù)十分困難。

針對(duì)上述兩個(gè)問題，首先采用NREL提出的NWTCUP湍流譜模型，該模型的特殊之處在于結(jié)合了SMOOTH[11-12]風(fēng)譜模型和San Gorgonio風(fēng)電場(chǎng)實(shí)測(cè)風(fēng)速。其中SMOOTH風(fēng)譜模型適合于地表粗糙度較低的地形，而San Gorgonio風(fēng)電場(chǎng)與我國東海附近風(fēng)電場(chǎng)緯度、氣候和地形均十分接近，具有地表粗糙度低、能量密度高和湍流度強(qiáng)等特點(diǎn)。該風(fēng)譜模型一定程度上可以代表我國東部近海湍流風(fēng)。同時(shí)，為表示時(shí)有發(fā)生的風(fēng)速變化突然加劇的風(fēng)況，在基礎(chǔ)湍流風(fēng)上加入相干結(jié)構(gòu)，以增強(qiáng)風(fēng)速的擾動(dòng)程度。通過該方法建立相應(yīng)的湍流風(fēng)，在專用于風(fēng)力機(jī)氣動(dòng)伺服彈性仿真的FAST[13]軟件中計(jì)算風(fēng)力機(jī)在開機(jī)啟動(dòng)、緊急停機(jī)和偏航等非穩(wěn)定工況的結(jié)構(gòu)動(dòng)力學(xué)響應(yīng)。為克服FFT方法和小波分析可能存在的不足，采用HHT(Hilbert-Huang Transform)方法，同時(shí)從時(shí)域和頻域分析風(fēng)力機(jī)動(dòng)態(tài)特性。該方法通過瞬時(shí)相位獲得瞬時(shí)響應(yīng)頻率，不受Heisenberg測(cè)不準(zhǔn)原則約束，可同時(shí)在時(shí)域和頻域達(dá)到較高精度[14-15]。通過對(duì)NREL 5 MW樣機(jī)進(jìn)行非穩(wěn)定工況的結(jié)構(gòu)動(dòng)力學(xué)分析，以期為風(fēng)力機(jī)非穩(wěn)定工況安全運(yùn)行提供一定的參考價(jià)值。

1 湍流風(fēng)場(chǎng)建模

1.1 風(fēng)譜模型和相干結(jié)構(gòu)

1.1.1 風(fēng)譜模型

NWTC/LIST項(xiàng)目[16]是由NWTC和GE Wind 合作，于洛杉磯附近的San Gorgonio風(fēng)電場(chǎng)開展的一項(xiàng)低海拔風(fēng)速測(cè)量項(xiàng)目，NWTC基于實(shí)測(cè)風(fēng)場(chǎng)數(shù)據(jù)和SMOOTH風(fēng)譜模型建立了NWTCUP風(fēng)譜模型，通過比例縮放SMOOTH風(fēng)模型功率譜密度函數(shù)，具體為[17]：

(1)

式中:f為周期頻率，K表示風(fēng)向，K=u,v,w；縮放系數(shù)Pi,K和Fi,K(i=1,2)為經(jīng)驗(yàn)總結(jié)，其函數(shù)形式十分復(fù)雜，具體見文獻(xiàn)[18]。SK(f)為SMOOTH風(fēng)譜模型，具體為：

(2)

φe=(1+2.5ZL0.6)3/2

φm=1+4.7ZL

(3)

式中：ZL一般取值為0.4。

1.1.2 加入相干結(jié)構(gòu)

考慮到風(fēng)是一種典型的湍流運(yùn)動(dòng)，相干結(jié)構(gòu)(Coherent Structure)是湍流運(yùn)動(dòng)的基本特征，因此可以通過加入相干結(jié)構(gòu)的方法，建立風(fēng)速突然劇烈變化的湍流風(fēng)。KHB(Kelvin-Helmholtz Billow)流動(dòng)是一種典型的相干結(jié)構(gòu)，描述自然界連續(xù)流場(chǎng)中由于速度梯度而導(dǎo)致渦迅速卷起并破裂的較規(guī)則周期性的流動(dòng)現(xiàn)象，如圖1所示[19]。

圖1 Kelvin-Helmholtz Billow流動(dòng)Fig.1 Flow of Kelvin-Helmholtz Billow

因此，通過在時(shí)間和空間維度將KHB的數(shù)值模擬的速度場(chǎng)信息無量綱化，并加入到基礎(chǔ)湍流風(fēng)中以增強(qiáng)風(fēng)速大小和方向的脈動(dòng)，表示風(fēng)速劇烈變化的風(fēng)況。

根據(jù)所選研究對(duì)象風(fēng)力機(jī)結(jié)構(gòu)參數(shù)，設(shè)置以輪轂點(diǎn)為中心，200 m寬175 m高的風(fēng)場(chǎng)計(jì)算域，網(wǎng)格節(jié)點(diǎn)分別為17和19，如圖2所示。考慮空間相干性，通過湍流風(fēng)譜模型每一個(gè)節(jié)點(diǎn)的風(fēng)速分布，并對(duì)圖中選定陰影部分區(qū)域加入相干結(jié)構(gòu)。

圖2 風(fēng)場(chǎng)計(jì)算域及網(wǎng)格節(jié)點(diǎn)Fig.2 Wind field of simulation and the grid

首先，將相干結(jié)構(gòu)量綱化，假設(shè)選定區(qū)域共有m×n個(gè)節(jié)點(diǎn)，其中第i行的第j個(gè)節(jié)點(diǎn)基礎(chǔ)風(fēng)速為ui,j(t)，相干結(jié)構(gòu)無量綱速度為CohUi,j(t)，則加入相干結(jié)構(gòu)后的風(fēng)速Ui,j(t)為：

(4)

1.2 風(fēng)場(chǎng)模擬結(jié)果

模擬風(fēng)場(chǎng)平均風(fēng)速為11.4 m/s，時(shí)間步長(zhǎng)為0.01 s，總時(shí)長(zhǎng)為600 s，加入相干結(jié)構(gòu)的湍流風(fēng)在100 s、200 s、300 s、400 s、500 s和600 s的總風(fēng)速分布如圖3所示。

由圖3可知，不同時(shí)刻的風(fēng)速均呈現(xiàn)地面風(fēng)速小高空風(fēng)速大的分布規(guī)律，且具有一定的無序性。其中，200 s、300 s和400 s均加入了相干結(jié)構(gòu)，在風(fēng)輪旋轉(zhuǎn)區(qū)域風(fēng)速空間分布差異愈加明顯，而且存在較為明顯的漩渦，整體風(fēng)速更紊亂，說明了加入相干結(jié)構(gòu)的有效性。

圖4為輪轂點(diǎn)三個(gè)方向風(fēng)速分量時(shí)域變化曲線。

圖3 各時(shí)刻風(fēng)場(chǎng)計(jì)算域的總風(fēng)速分布Fig.3 Wind speeddistributions of the field in various moment

圖4 輪轂處風(fēng)速分布Fig.4 Wind speed in time domain at hub point

由圖4可知，相干結(jié)構(gòu)加劇了三個(gè)方向的風(fēng)速變化，各個(gè)分量的大小和方向均發(fā)生較大變化，其中垂直方向速度w劇烈程度增大最為明顯，風(fēng)速脈動(dòng)劇烈。

圖5為加入相干結(jié)構(gòu)前后風(fēng)速頻譜，由圖5可知，加入相干結(jié)構(gòu)后增大了風(fēng)在高于0.1 Hz頻率的能量。

圖5 輪轂點(diǎn)處風(fēng)速頻譜響應(yīng)Fig.5 Frequency spectrum of wind speed at hub point

為了表明相干結(jié)構(gòu)加入的可行性，通過計(jì)盒維數(shù)法計(jì)算輪轂點(diǎn)處總風(fēng)速的分形維數(shù)，具體算法參見文獻(xiàn)[20]。計(jì)算結(jié)果為：未加入相干結(jié)構(gòu)的分形維數(shù)為1.520 6，加入相干結(jié)構(gòu)后為1.548 0，二者均在1.5左右，且相差很小，表明了所建立的湍流風(fēng)具有非常明顯的混沌特征和自然屬性。

2 研究對(duì)象及研究方法

2.1 風(fēng)力機(jī)模型

大型風(fēng)力機(jī)動(dòng)態(tài)響應(yīng)非線性更加明顯，因此選用由NREL開發(fā)的額定功率為5 MW的風(fēng)力機(jī)作為研究對(duì)象[21]，其額定風(fēng)速和額定風(fēng)輪轉(zhuǎn)速分別為11.4 m/s和12.1 r/min，主要結(jié)構(gòu)參數(shù)和風(fēng)力機(jī)實(shí)體如圖6所示。

圖6 NREL 5 MW風(fēng)力機(jī)結(jié)構(gòu)參數(shù)示意圖Fig.6 The structural parameters of the NREL 5 MW wind turbine

2.2 研究方法

模態(tài)截?cái)喾ㄅc多體動(dòng)力學(xué)相結(jié)合的計(jì)算方法是目前實(shí)現(xiàn)風(fēng)力機(jī)結(jié)構(gòu)動(dòng)力學(xué)仿真的主流方式之一。該方法假設(shè)風(fēng)力機(jī)為有限個(gè)剛性體和柔性體的組合系統(tǒng)，通過形函數(shù)描述風(fēng)力機(jī)柔性體的模態(tài)振型，考慮的模態(tài)數(shù)量即為風(fēng)力機(jī)系統(tǒng)的自由度數(shù)目，從而加快氣動(dòng)結(jié)構(gòu)耦合計(jì)算速度。其中，由NWTC(National Wind Technology Center)針對(duì)水平軸風(fēng)力機(jī)研發(fā)的開源軟件FAST應(yīng)用最為廣泛[22-23]。

FAST是耦合氣動(dòng)-伺服-彈性在時(shí)域求解風(fēng)力機(jī)結(jié)構(gòu)動(dòng)力學(xué)響應(yīng)的CAE軟件，主要包含3個(gè)模塊：氣動(dòng)模塊(AeroDyn)、彈性模塊(ElastoDyn)和伺服控制模塊(ServoDyn)。其中，氣動(dòng)模塊采用Pitt-Peters加速度勢(shì)動(dòng)態(tài)入流理論，求解風(fēng)輪平面誘導(dǎo)速度；考慮Prandtl葉尖損失及葉輪損失，通過葉素動(dòng)量理論結(jié)合翼型靜態(tài)氣動(dòng)力特性求解風(fēng)輪氣動(dòng)力，若需要翼型動(dòng)態(tài)氣動(dòng)特性則通過Beddoes-Leishman動(dòng)態(tài)失速模型修正。在彈性模塊中，通過Kane方法建立多體動(dòng)力學(xué)模型，將風(fēng)力機(jī)視為由葉片、低速軸和塔架等柔性體及輪轂、變速箱、高速軸、發(fā)電機(jī)和機(jī)艙等剛性體組成的多結(jié)構(gòu)體系統(tǒng)。采用模態(tài)截?cái)喾枋鋈~片和塔架等柔性連續(xù)體彈性變形，假設(shè)其結(jié)構(gòu)變形為一系列振動(dòng)模態(tài)的線性疊加，以氣動(dòng)模塊求解的風(fēng)輪氣動(dòng)力作為輸入激勵(lì)，得到該時(shí)間步的結(jié)構(gòu)動(dòng)力學(xué)及運(yùn)動(dòng)學(xué)響應(yīng)并反饋至伺服模塊和氣動(dòng)模塊。伺服模塊則根據(jù)彈性模塊反饋信息作出相應(yīng)的控制指令，主要包括調(diào)節(jié)葉片槳距角、風(fēng)輪轉(zhuǎn)速和高速軸轉(zhuǎn)速等。具體仿真流程如圖7所示，其中Tmax為仿真時(shí)間，dt為時(shí)間步長(zhǎng)。

圖7 動(dòng)力學(xué)仿真流程圖Fig.7 The flowchart of dynamic simulation

3 HHT數(shù)據(jù)分析方法

HHT是美國NASA研究中心的Norden Huang首次提出的一種分析非平穩(wěn)數(shù)據(jù)的方法，具有自適應(yīng)性局部時(shí)頻特性分析的能力，對(duì)復(fù)雜的風(fēng)力機(jī)非穩(wěn)定的非線性氣動(dòng)彈性響應(yīng)和結(jié)構(gòu)動(dòng)力學(xué)特性具有更高的時(shí)頻特性分辨率。

HHT方法主要包含兩部分：經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)和Hilbert 變換[14-15]。首先通過EMD方法提取時(shí)間序列的固有模態(tài)函數(shù)(Intrinsic Mode Function，IMF)，對(duì)每組IMF進(jìn)行Hilbert 變換得到對(duì)應(yīng)的Hilbert譜，整合所有的Hilbert譜即可得到對(duì)應(yīng)時(shí)間序列的時(shí)頻特性。

3.1 EMD方法

EMD方法將時(shí)間序列數(shù)據(jù)分解成一系列具有不同特征尺度的數(shù)據(jù)系列，即IMFs。對(duì)于任意具體的時(shí)間序列數(shù)據(jù)X(t)，獲取IMFs的具體步驟如下：

(1) 分別使用樣條曲線擬合所有的極大值和極小值，得到該數(shù)據(jù)曲線的上包絡(luò)線Eup(t)、下包絡(luò)線Elow(t)，包絡(luò)線均值Emean(t)=[Eup(t)+Elow(t)]/2。

(2) 原始數(shù)據(jù)X(t)減去包絡(luò)線均值Emean(t)得到第一個(gè)IMF初始值h1(t)。

(3) 將h1(t)視為初始數(shù)據(jù)，重復(fù)步驟(1)～步驟(2)k次直至滿足IMF條件，得到第一個(gè)IMF：h1k(t)，判斷準(zhǔn)則為方差SD值在0.2～0.3之間,SD通過下式計(jì)算：

(5)

(4) 記第一個(gè)IMF為c1(t)，余項(xiàng)r1(t)為：

r1(t)=X(t)-c1(t)

(6)

將余項(xiàng)r1(t)作為原始數(shù)據(jù)，重復(fù)步驟(1)～步驟(4)，直至余項(xiàng)rn(t)成為單調(diào)函數(shù)或定值(即極值點(diǎn)小于2個(gè))，則分解結(jié)束。原始數(shù)據(jù)X(t)可表示為：

(7)

3.2 Hilbert 變換

(8)

式中:P為柯西主值，一般取值1。

則c(t)的解析信號(hào)ca(t)為：

(9)

對(duì)X(t)的每組IMF進(jìn)行Hilbert變換可得：

(10)

將時(shí)間t、頻率ω(t)/2π和幅值a(t)表示在同一圖中即可得到X(t)的時(shí)頻響應(yīng)。

4 結(jié)果分析

4.1 有效性驗(yàn)證

為驗(yàn)證計(jì)算結(jié)果的有效性，與風(fēng)力機(jī)動(dòng)態(tài)載荷專業(yè)計(jì)算軟件GH Bladed比較，圖8為葉尖擺振加速度和揮舞加速度動(dòng)態(tài)響應(yīng)對(duì)比。

圖8 葉尖加速度比較Fig.8 Comparison for blade tip accelerations

從圖8中可以看出，本文計(jì)算結(jié)果與GH Bladed結(jié)果趨勢(shì)相同，均能明顯反映出開機(jī)啟動(dòng)、偏航和緊急停機(jī)的非平穩(wěn)動(dòng)態(tài)響應(yīng)過程，且響應(yīng)幅值大小基本一致，驗(yàn)證了本文方法的有效性。

4.2 時(shí)域結(jié)果分析

在FAST中分別計(jì)算NREL 5 MW風(fēng)力機(jī)在加入相干結(jié)構(gòu)前后的湍流風(fēng)場(chǎng)下的動(dòng)態(tài)響應(yīng)，仿真總時(shí)間為600 s，時(shí)間步長(zhǎng)為0.005 s。為反映非穩(wěn)定工況的動(dòng)態(tài)特性，設(shè)定起始風(fēng)輪轉(zhuǎn)速為0，仿真風(fēng)力機(jī)自啟動(dòng)工況。在200 s時(shí)開始偏航運(yùn)動(dòng)，偏航速度為3°/s，持續(xù)10 s，并在400 s時(shí)緊急順槳停機(jī)，風(fēng)輪和發(fā)電機(jī)轉(zhuǎn)速均降為0。

圖9為風(fēng)輪轉(zhuǎn)速、發(fā)電機(jī)功率、塔尖位移和葉尖位移的時(shí)域仿真結(jié)果。由圖9可知，在開機(jī)啟動(dòng)階段，風(fēng)輪轉(zhuǎn)速緩慢上升，尖速比逐漸增大，風(fēng)力機(jī)氣動(dòng)性能逐漸提高，大約40 s后風(fēng)力機(jī)因受力較大，風(fēng)輪轉(zhuǎn)速迅速上升至額定轉(zhuǎn)速，從55 s開始在額定轉(zhuǎn)速附近微小波動(dòng)。表明已開機(jī)成功，此時(shí)處于正常運(yùn)行工況。200 s時(shí)執(zhí)行偏航運(yùn)動(dòng)并持續(xù)10 s，由于并未停機(jī)，風(fēng)力機(jī)仍盡可能輸出功率，但由于處于偏航狀態(tài)，葉片實(shí)際攻角范圍內(nèi)的氣動(dòng)性能相對(duì)較差，風(fēng)輪功率略有降低。400 s時(shí)緊急停機(jī)，風(fēng)輪轉(zhuǎn)速急劇降低，風(fēng)輪功率出現(xiàn)劇烈波動(dòng)，約20 s后穩(wěn)定至0，說明此時(shí)停機(jī)成功。

由塔尖位移和葉尖位移動(dòng)態(tài)響應(yīng)結(jié)果可知，不同運(yùn)行狀態(tài)其振動(dòng)形態(tài)差異較大。開機(jī)啟動(dòng)階段，結(jié)構(gòu)位移均為0左右，啟動(dòng)成功后出現(xiàn)微小振蕩，約為最大位移的1/3左右。正常運(yùn)行工況時(shí)塔架和葉片位移均在相對(duì)穩(wěn)定的范圍內(nèi)變化，偏航時(shí)塔尖位移出現(xiàn)突然變化，但逐漸回穩(wěn)。在緊急停機(jī)時(shí)，各結(jié)構(gòu)動(dòng)態(tài)響應(yīng)除葉尖擺振外，均急劇降低后逐漸恢復(fù)至0值附近。

圖9 風(fēng)輪轉(zhuǎn)速、發(fā)電機(jī)功率以及塔架和葉片位移時(shí)域響應(yīng)Fig.9 Response of rotor speed, electric power, tower tip and blade tip deflections in time domain

各響應(yīng)結(jié)果說明加入相干結(jié)構(gòu)增強(qiáng)了風(fēng)力機(jī)動(dòng)態(tài)響應(yīng)的非線性和非穩(wěn)定特征，其中，相干結(jié)構(gòu)對(duì)塔尖位移的影響很大，對(duì)葉尖位移的影響相對(duì)較小。塔尖側(cè)向位移在正常運(yùn)行狀態(tài)下受到相干結(jié)構(gòu)的極大影響，變化范圍增大8倍。偏航后，由于風(fēng)輪旋轉(zhuǎn)30°，受到風(fēng)力機(jī)重心偏移的影響，相干結(jié)構(gòu)對(duì)其影響相對(duì)較小。相干結(jié)構(gòu)對(duì)葉尖擺振影響最小，但在停機(jī)后，前期的微弱影響極大改變了擺振的動(dòng)態(tài)響應(yīng)，兩種風(fēng)況出現(xiàn)巨大反差。

4.3 時(shí)頻特性分析

風(fēng)輪轉(zhuǎn)速和功率的整體時(shí)域變化規(guī)律驗(yàn)證了三種非穩(wěn)定工況的特征響應(yīng)，一定程度地反映了仿真設(shè)置的有效性。塔尖位移和葉尖位移的動(dòng)態(tài)響應(yīng)說明了不同運(yùn)行狀態(tài)下塔架和葉片結(jié)構(gòu)動(dòng)力學(xué)響應(yīng)的非線性和非平穩(wěn)特征，加入相干結(jié)構(gòu)后的湍流風(fēng)場(chǎng)則進(jìn)一步增大了結(jié)構(gòu)振動(dòng)，對(duì)于此類非線性非平穩(wěn)的振動(dòng)分析，通過FFT方法無法準(zhǔn)確甄別引起非平穩(wěn)時(shí)域響應(yīng)的特征頻率。因此，通過HHT方法對(duì)風(fēng)力機(jī)動(dòng)態(tài)響應(yīng)進(jìn)行分析，以準(zhǔn)確獲得其時(shí)頻響應(yīng)。圖10為加入相干結(jié)構(gòu)工況的塔尖前后位移和葉尖擺振EMD分解結(jié)果。

(a) 塔尖前后位移

(b) 葉尖擺振圖10 塔尖前后位移和葉尖擺振EMD結(jié)果Fig.10 EMD results of tower tip deflection at fore-after direction and blade tip deflection in plane

從圖10中可以看出，塔尖前后位移和葉尖擺振通過EMD處理后的余項(xiàng)均為0，表明其動(dòng)態(tài)響應(yīng)由若干組IMF代替。對(duì)于塔尖前后位移，其中第一組IMF(c1)主要表示湍流風(fēng)和偏航的影響，開機(jī)和停機(jī)的振動(dòng)主要體現(xiàn)在c3～c5，更高階的IMF表示全局低頻振動(dòng)。湍流風(fēng)和偏航對(duì)葉尖擺振的影響主要體現(xiàn)在c1～c3，高階IMF頻率較低幅值較大。

將IMFs分別進(jìn)行Hilbert變換，得到的時(shí)頻響應(yīng)Hilbert譜，如圖10所示。

圖11 塔尖前后位移和葉尖擺振Hilbert譜Fig.11 Hilbert spectrums of tower tip deflection at fore-after direction and blade tip deflection in plane

由圖11(a)可知，塔尖前后位移在塔架一階振動(dòng)頻率0.32 Hz附近有明顯的響應(yīng)，且為典型的非穩(wěn)定響應(yīng)，幅值隨著時(shí)間不斷變化，其中在開機(jī)啟動(dòng)階段，響應(yīng)頻率相對(duì)集中但幅值較小，44 s風(fēng)力機(jī)成功啟動(dòng)后，響應(yīng)頻率出現(xiàn)劇烈波動(dòng)，這是因?yàn)闅鈩?dòng)載荷阻尼增大使得結(jié)構(gòu)阻尼作用降低。200 s時(shí)進(jìn)行主動(dòng)偏航運(yùn)動(dòng)，產(chǎn)生的載荷突變使得結(jié)構(gòu)阻尼作用增大，因此其幅值增大且響應(yīng)頻率集中。同時(shí)由于偏航后一直處于相干湍流風(fēng)作用，0.32 Hz頻率附近風(fēng)的能量更大，諧振影響增強(qiáng)，相比于150 s之前，該階段的響應(yīng)幅值更大。

圖11(b)說明葉尖擺振的諧振頻率明顯集中于葉輪旋轉(zhuǎn)頻率，在開機(jī)啟動(dòng)階段，風(fēng)輪轉(zhuǎn)速逐漸增大，其響應(yīng)頻率隨之增大，在44 s開機(jī)成功后，風(fēng)輪轉(zhuǎn)速比較穩(wěn)定，響應(yīng)頻率也在0.2 Hz左右波動(dòng)。在400 s后，風(fēng)輪轉(zhuǎn)速急劇降低至0，該頻段的響應(yīng)消失。說明葉尖擺振的刺激頻率主要為葉輪旋轉(zhuǎn)頻率，因此風(fēng)速大小和變化對(duì)其影響有限，從側(cè)面佐證了相干結(jié)構(gòu)對(duì)葉尖擺振時(shí)域響應(yīng)影響很小。除此之外，葉尖擺振在一階擺振頻率1.1 Hz附近具有明顯的非平穩(wěn)動(dòng)態(tài)響應(yīng)，隨時(shí)間變化其響應(yīng)頻率在1.1 Hz上下劇烈波動(dòng)，并在緊急停機(jī)時(shí)頻率中心下降至0.7 Hz。因?yàn)榧尤胂喔山Y(jié)構(gòu)后增大了固有頻率激勵(lì)載荷，使得葉尖擺振在該頻段響應(yīng)幅值增大，而緊急停機(jī)時(shí)風(fēng)輪轉(zhuǎn)速突然降低且葉片順槳，導(dǎo)致在升力方向的葉片運(yùn)動(dòng)產(chǎn)生負(fù)氣動(dòng)阻尼，而此時(shí)葉片正處于較明顯的一階擺振模態(tài)，由于順槳產(chǎn)生的負(fù)氣動(dòng)阻尼大于其結(jié)構(gòu)阻尼，由于之前的擾動(dòng)導(dǎo)致葉尖擺振產(chǎn)生增幅振蕩，從側(cè)面佐證了時(shí)域幅值響應(yīng)增大的現(xiàn)象。這一現(xiàn)象也說明在葉片設(shè)計(jì)時(shí)，應(yīng)當(dāng)考慮緊急順槳停機(jī)使得葉片擺振方向位移幅值增大的現(xiàn)象，適當(dāng)增加葉片阻尼大小，以緩解氣動(dòng)阻尼減小帶來的影響。

5 結(jié) 論

基于NWTCUP湍流風(fēng)譜模型建立了基礎(chǔ)湍流風(fēng)場(chǎng)，無量綱化Kelvin-Helmholtz Billow的數(shù)值模擬結(jié)構(gòu)，以相干結(jié)構(gòu)的形式加入到基礎(chǔ)流中，構(gòu)建了湍動(dòng)劇烈的湍流風(fēng)場(chǎng)作為風(fēng)力機(jī)動(dòng)力學(xué)仿真環(huán)境，在FAST軟件中計(jì)算了風(fēng)力機(jī)開機(jī)啟動(dòng)、偏航和緊急停機(jī)三種特征非穩(wěn)定過程，并與GH Bladed計(jì)算結(jié)果比較以驗(yàn)證研究方法的有效性，通過HHT方法對(duì)塔尖前后方向的位移非平穩(wěn)響應(yīng)進(jìn)行分析，主要得到以下結(jié)論：

(1) 加入相干結(jié)構(gòu)后的風(fēng)場(chǎng)紊亂度更強(qiáng)，風(fēng)速變化更劇烈，高于0.1 Hz頻率的能量更大，其分形維數(shù)為1.548 0，具有非常明顯的分形特征和混沌屬性，表明模型建立的正確性。

(2) 開機(jī)啟動(dòng)階段塔架和葉片位移均在0值附近緩慢變化，大約40 s后響應(yīng)幅值急劇增加。在偏航運(yùn)動(dòng)時(shí)，塔尖側(cè)向位移出現(xiàn)明顯的階躍響應(yīng)。緊急停機(jī)時(shí)，塔架和葉片位移(除葉尖擺振外)均突然降低并在短時(shí)間內(nèi)持續(xù)振蕩，約20 s停機(jī)成功后逐漸恢復(fù)至0值附近小幅度變化。

(3) HHT分析結(jié)果表示，塔架一階振動(dòng)頻率0.32 Hz下的響應(yīng)幅值變化為明顯的非平穩(wěn)過程，在開機(jī)啟動(dòng)和停機(jī)后響應(yīng)頻率尤為集中，偏航時(shí)結(jié)構(gòu)阻尼作用增強(qiáng)，響應(yīng)幅值明顯增大。葉尖擺振主要刺激頻率為風(fēng)輪旋轉(zhuǎn)頻率，在一階擺振頻率下的非平穩(wěn)響應(yīng)主要受到相干結(jié)構(gòu)影響，在緊急停機(jī)時(shí)由于葉片順槳產(chǎn)生的負(fù)氣動(dòng)阻尼作用而導(dǎo)致響應(yīng)幅值增大，在葉片設(shè)計(jì)時(shí)應(yīng)當(dāng)考慮這一現(xiàn)象并適當(dāng)增大其阻尼值，以減小氣動(dòng)阻尼迅速降低帶來的振幅急速增大的現(xiàn)象。

[1] STEVE S, KLAUS R. Global wind report-annual market update 2014[R]. Brussels: GWEC, 2015.

[2] 岳一松，蔡旭. 風(fēng)場(chǎng)與風(fēng)力機(jī)模擬系統(tǒng)的設(shè)計(jì)與實(shí)現(xiàn)[J]. 電機(jī)與控制應(yīng)用，2008，35(4)：17-21. YUE Yisong，CAI Xu．Design and actualization of wind farm and wind turbine imitation system[J]．Electric Machines & Control Application, 2008，35(4)：17-21．

[3] 陳曉明．風(fēng)場(chǎng)與風(fēng)力機(jī)尾流模型研究[D]．蘭州：蘭州理工大學(xué)，2010．

[4] BOSSANYI E. Short-term stochastic wind prediction and possible control applications[C]//Proceedings of the Delphi workshop on wind energy applications, 1985.

[5] BILGILI M, SAHIN B, YASAR A. Application of artificial neural networks for the wind speed prediction of target station using reference stations data[J]. Renewable Energy, 2007, 32(14): 2350-2360.

[6] DAMOUSIS I G, ALEXIADIS M C, THEOCHARIS J B, et al. A fuzzy model for wind speed prediction and power generation in wind parks using spatial correlation[J]. Energy Conversion, IEEE Transactions on, 2004, 19(2): 352-361.

[7] RATHMANN O, MORTENSEN N G, LANDBERG L. The numerical wind atlas-the KAMM/WAsP method[R]. Riso National Laboratory, 2001.

[8] FLEMING P, GEBRAAD P, VAN WINGERDEN J W, et al. The SOWFA super-controller: a high-fidelity tool for evaluating wind plant control approaches[C]//Proceedings of the EWEA Annual Meeting, Vienna, Austria. 2013.

[9] 徐磊，李德源，莫文威，等. 基于非線性氣彈耦合模型的風(fēng)力機(jī)柔性葉片隨機(jī)響應(yīng)分析[J].振動(dòng)與沖擊，2015，34(10)：20-27. XU Lei, LI Deyuan, MO Wenwei, et al. Random response analysis for flexible blade of a wind turbine based on nonlinear aero-elastic coupled model[J]. Journal of Vibration and Shock, 2015,34(10): 20-27.

[10] 柯世堂，王同光. 偏航狀態(tài)下風(fēng)力機(jī)塔架-葉片耦合結(jié)構(gòu)氣彈響應(yīng)分析[J].振動(dòng)與沖擊，2015，34(18)：33-38. KE Shitang, WANG Tongguang. Aero-elastic vibration analysis based on a tower-blade coupled model of wind turbine in yaw motion[J]. Journal of Vibration and Shock, 2015, 34(18): 33-38.

[11] H?JSTRUP J. Velocity spectra in the unstable planetary boundary layer[J]. Journal of the Atmospheric Sciences, 1982, 39(10): 2239-2248.

[12] OLESEN H R, LARSEN S E, H?JSTRUP J. Modelling velocity spectra in the lower part of the planetary boundary layer[J].Boundary-Layer Meteorology,1984,29(3): 285-312.

[13] JONKMAN J M, BUHI JR M L. FAST user’s guide[R]. National Renewable Energy Laboratory, Golden, CO, NREL/EL-500-38230, 2005.

[14] HUANG N E, SHEN Z, LONG S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[C]//Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. The Royal Society, 1998, 454(1971): 903-995.

[15] HUANG N E, WU Z. A review on hilbert-huang transform: method and its applications to geophysical studies[J]. Reviews of Geophysics, 2008, 46(2).

[16] SHIRAZI M, JAGER D, WILDE S, et al. Lamar low-level jet project interim report[R]. Golden, CO, USA: National Renewable Energy Laboratory, 2004.

[17] KELLEY N, HAND M, LARWOOD S, et al. The NREL large-scale turbine inflow and response experiment: preliminary results[C]//ASME 2002 Wind Energy Symposium. American Society of Mechanical Engineers, 2002: 412-426.

[18] SHIRAZI M, JAGER D, WILDE S, et al. Lamar low-level jet project interim report[R]. Golden, CO, USA: National Renewable Energy Laboratory, 2004.

[19] WIKIPERIA. Kelvin-Helmholtz instability[EB/OL]. [2015-07-30]. https://en.wikipedia.org/wiki/Kelvin%E2%80%93Helmholtz_instability.

[20] CHANG T P, KO H H, LIU F J, et al. Fractal dimension of wind speed time series[J]. Applied Energy,2012,93: 742-749.

[21] JONKMAN J, BUTTERFIELD S, MUSIAL W, et al. Definition of a 5-MW reference wind turbine for offshore system development [R]. National Renewable Energy Laboratory, 2009, NREL/TP 500-38060.

[22] 吳攀，李春，李志敏，等. 風(fēng)力機(jī)不同風(fēng)況的動(dòng)力學(xué)響應(yīng)研究[J].中國電機(jī)工程學(xué)報(bào)，2014，34(26)：4539-4545. WU Pan, LI Chun, LI Zhimin, et al. Research on dynamic characteristics simulation for wind turbine with different wind[J].Proceedings of the CSEE, 2014, 34(26): 4539-4545.

[23] 王磊，陳柳，何玉林，等. 基于假設(shè)模態(tài)法的風(fēng)力機(jī)動(dòng)力學(xué)分析[J].振動(dòng)與沖擊，2012，31(11)：122-126. WANG Lei, CHEN Liu, HE Yulin, et al. Dynamic analysis of a wind turbine base on assumed mode method[J]. Journal of Vibration and Shock, 2012, 31(11): 122-126.

Structural dynamic response characteristics of a wind turbine in time-frequency domain under non-stationary operating conditions based on HHT method

YANG Yang1, LI Chun1,2, YE Kehua1, MIAO Weipao1, YANG Jun1,2, GAO Wei3

(1. School of Energy and Power Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China;2. Shanghai Key Laboratory of Multiphase Flow and Heat Transfer in Power Engineering, Shanghai 200093, China;3. China Creative Wind Energy Co., Ltd., Qingdao 266000, China)

The typical non-stationary operating conditions, such as, turbulent inflow of wind speed changing dramatically, starting up, yawing motion, and emergency shutdown enhance nonlinear aero-elastic responses of a wind turbine, its structural dynamic responses in time domain and frequency domain both have obvious non-stationary characteristics. The turbulent inflow with velocity and direction both varying intensely was built based on its spectral model and coherent structure. The dynamic characteristics of the wind turbine under non-stationary operating conditions were computed with the aerodynamic-servo-elastic software called FAST. The numerical calculation results were compared with those of GH Bladed to verify the results’ validity. The dynamic characteristics in time-frequency domain for blades and tower deflections were analyzed based on Hilbert-Huang Transformation (HHT) method. The results showed that the deflections of tower and blade tip fluctuate within a narrow range of about 40 s and then increase rapidly in starting up, but they fluctuate wildly about 20 s and return a steady state after emergency shutdown; yawing motion leads to a clearly rise of tower tip deflection in lateral direction; the fluctuation of tower tip deflection is concentrated at the first order vibration frequency and its amplitude increases in yawing motion; the fluctuation of blade tip deflection in plane is mainly caused by the rotor rotating and blades’ the first order shimmy frequency is affected by coherent structure; the negative aerodynamic damping makes amplitudes rising during emergency shutdown, enlarging structural damping should be considered in the design of blades to depress their amplitude rising due to aerodynamic damping rapidly dropping.

wind turbine; turbulent inflow; structural dynamic characteristics; HHT method

國家自然科學(xué)基金資助項(xiàng)目(51176129;51506126)；上海市教育委員會(huì)科研創(chuàng)新(重點(diǎn))項(xiàng)目(13ZZ120;13YZ066)；教育部高等學(xué)校博 士學(xué)科點(diǎn)專項(xiàng)科研基金(博導(dǎo)類)項(xiàng)目(20123120110008)；上海市科委項(xiàng)目資助(13DZ2260900);上海市研究生創(chuàng)新基金項(xiàng)目(JWCXSL1402)

2015-09-07 修改稿收到日期：2015-10-15

楊陽 男，博士生，1992年生

李春 男，教授,博導(dǎo)，1963年生

TK83

A

10.13465/j.cnki.jvs.2016.21.004