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

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

?

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

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

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

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

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

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

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

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

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

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

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

1.1.1 風(fēng)譜模型

NWTC/LIST項目[16]是由NWTC和GE Wind 合作，于洛杉磯附近的San Gorgonio風(fēng)電場開展的一項低海拔風(fēng)速測量項目，NWTC基于實測風(fē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)驗總結(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)是一種典型的湍流運動，相干結(jié)構(gòu)(Coherent Structure)是湍流運動的基本特征，因此可以通過加入相干結(jié)構(gòu)的方法，建立風(fēng)速突然劇烈變化的湍流風(fēng)。KHB(Kelvin-Helmholtz Billow)流動是一種典型的相干結(jié)構(gòu)，描述自然界連續(xù)流場中由于速度梯度而導(dǎo)致渦迅速卷起并破裂的較規(guī)則周期性的流動現(xiàn)象，如圖1所示[19]。

圖1 Kelvin-Helmholtz Billow流動Fig.1 Flow of Kelvin-Helmholtz Billow

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

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

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

首先，將相干結(jié)構(gòu)量綱化，假設(shè)選定區(qū)域共有m×n個節(jié)點，其中第i行的第j個節(jié)點基礎(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)場模擬結(jié)果

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

由圖3可知，不同時刻的風(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為輪轂點三個方向風(fēng)速分量時域變化曲線。

圖3 各時刻風(fēng)場計算域的總風(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)加劇了三個方向的風(fēng)速變化，各個分量的大小和方向均發(fā)生較大變化，其中垂直方向速度w劇烈程度增大最為明顯，風(fēng)速脈動劇烈。

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

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

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

2 研究對象及研究方法

2.1 風(fēng)力機模型

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

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

2.2 研究方法

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

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

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

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

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

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

3.1 EMD方法

EMD方法將時間序列數(shù)據(jù)分解成一系列具有不同特征尺度的數(shù)據(jù)系列，即IMFs。對于任意具體的時間序列數(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)得到第一個IMF初始值h1(t)。

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

(5)

(4) 記第一個IMF為c1(t)，余項r1(t)為：

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

(6)

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

(7)

3.2 Hilbert 變換

(8)

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

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

(9)

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

(10)

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

4 結(jié)果分析

4.1 有效性驗證

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

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

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

4.2 時域結(jié)果分析

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

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

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

圖9 風(fēng)輪轉(zhuǎn)速、發(fā)電機功率以及塔架和葉片位移時域響應(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)增強了風(fēng)力機動態(tài)響應(yīng)的非線性和非穩(wěn)定特征，其中，相干結(jié)構(gòu)對塔尖位移的影響很大，對葉尖位移的影響相對較小。塔尖側(cè)向位移在正常運行狀態(tài)下受到相干結(jié)構(gòu)的極大影響，變化范圍增大8倍。偏航后，由于風(fēng)輪旋轉(zhuǎn)30°，受到風(fēng)力機重心偏移的影響，相干結(jié)構(gòu)對其影響相對較小。相干結(jié)構(gòu)對葉尖擺振影響最小，但在停機后，前期的微弱影響極大改變了擺振的動態(tài)響應(yīng)，兩種風(fēng)況出現(xiàn)巨大反差。

4.3 時頻特性分析

風(fēng)輪轉(zhuǎn)速和功率的整體時域變化規(guī)律驗證了三種非穩(wěn)定工況的特征響應(yīng)，一定程度地反映了仿真設(shè)置的有效性。塔尖位移和葉尖位移的動態(tài)響應(yīng)說明了不同運行狀態(tài)下塔架和葉片結(jié)構(gòu)動力學(xué)響應(yīng)的非線性和非平穩(wěn)特征，加入相干結(jié)構(gòu)后的湍流風(fēng)場則進(jìn)一步增大了結(jié)構(gòu)振動，對于此類非線性非平穩(wěn)的振動分析，通過FFT方法無法準(zhǔn)確甄別引起非平穩(wěn)時域響應(yīng)的特征頻率。因此，通過HHT方法對風(fēng)力機動態(tài)響應(yīng)進(jìn)行分析，以準(zhǔn)確獲得其時頻響應(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處理后的余項均為0，表明其動態(tài)響應(yīng)由若干組IMF代替。對于塔尖前后位移，其中第一組IMF(c1)主要表示湍流風(fēng)和偏航的影響，開機和停機的振動主要體現(xiàn)在c3～c5，更高階的IMF表示全局低頻振動。湍流風(fēng)和偏航對葉尖擺振的影響主要體現(xiàn)在c1～c3，高階IMF頻率較低幅值較大。

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

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

5 結(jié) 論

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

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

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

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

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

[2] 岳一松，蔡旭. 風(fēng)場與風(fēng)力機模擬系統(tǒng)的設(shè)計與實現(xiàn)[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)場與風(fēng)力機尾流模型研究[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)力機柔性葉片隨機響應(yīng)分析[J].振動與沖擊，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)力機塔架-葉片耦合結(jié)構(gòu)氣彈響應(yīng)分析[J].振動與沖擊，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)力機不同風(fēng)況的動力學(xué)響應(yīng)研究[J].中國電機工程學(xué)報，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)力機動力學(xué)分析[J].振動與沖擊，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é)基金資助項目(51176129;51506126)；上海市教育委員會科研創(chuàng)新(重點)項目(13ZZ120;13YZ066)；教育部高等學(xué)校博 士學(xué)科點專項科研基金(博導(dǎo)類)項目(20123120110008)；上海市科委項目資助(13DZ2260900);上海市研究生創(chuàng)新基金項目(JWCXSL1402)

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

楊陽 男，博士生，1992年生

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

TK83

A

10.13465/j.cnki.jvs.2016.21.004