湍流風與地震聯合作用下風力機塔架振動非線性特征研究

鄒錦華， 楊 陽， 李 春,2， 劉中勝， 袁全勇

(1. 上海理工大學 能源與動力工程學院， 上海 200093； 2. 上海市動力工程多相流動與傳熱重點實驗室， 上海 200093)

“十二五”規劃期間，我國風電年均增長34.6%，為第二大非化石能源，2016年新增裝機量為23 370 MW，總裝機量達168 732 MW。“十三五”時期計劃將風電占全國總發電量的比例提高至6%以實現非化石能源占一次能源消費比重達15%的重大目標[1-2]。至2020年，我國中東部和南方地區陸上風電新增并網裝機容量為42 000 MW以上，西北地區累計并網容量將達到48 500 MW[3]。上述區域分別處于環太平洋地震帶及亞歐地震帶附近，風電場極有可能受到地震影響。大型風力機塔架力學結構屬于典型的細長柔性體，加之塔架頂部的大質量機艙，在空間非均勻和時間非定常的湍流風作用下，塔架氣彈響應具有明顯的非線性和非平穩特征[4-5]，若同時受到高強度地震作用，可能發生屈曲及結構永久變形，影響結構安全[6]。因此，開展湍流風與地震聯合作用下風力機塔架動力學響應研究對風力機設計及安全評估具有重要意義。

關于地震作用對風力機塔架結構的影響，國內外學者開展了相關研究。Adam等[7]以1.5 MW風力機為研究對象，建立塔架有限元模型，分別研究了10種近斷層地震和遠場地震作用下塔架模態振型及在塑性變形時地震能量耗散情況，以5%的氣動阻尼代替湍流風作用，發現風力機主要受前兩階基礎振型的影響，且近斷層地震對塔架各段接縫處的破壞影響較大。季亮等[8]建立風力機單自由度和雙自由度簡化模型，運用剪力法計算8級地震作用下不同功率風力機的塔架底部彎矩與剪切力，因研究參考建筑結構設計方法簡化風力機模型，未能考慮風輪所受的非定常氣動載荷。Kim等[9]建立近海風力機樁-土體耦合模型，分析了土質對風力機地震動力學響應的影響，計算模型中忽略了氣動載荷的影響。Ikwulono等[10]采用故障樹方法分析風與地震分別獨立及聯合作用時風力機動力學響應，運用有限元軟件分析風載荷與地震載荷共同作用時塔基失效概率。李春等[11]基于開源風力機仿真軟件FAST建立湍流風-地震實時耦合模型，模擬多組湍流風況與地震共同作用風力機動態響應，發現地震激勵下風力機塔基載荷達到穩態風工況的2倍～15倍。

風力機的運行受地理環境、氣候條件及表面粗糙度等諸多因素的影響，加之湍流風與地震聯合激勵作用，導致風力機動態響應具有明顯非線性特征，因此，采用混沌理論這一非線性分析工具來更加準確揭示湍流風與地震聯合作用等復雜工況下塔架振動具有重要意義。

為研究風力機塔架結構在湍流風與地震聯合作用等不利情況下振動的非線性特性，以NREL 5 MW風力機為研究對象，根據Wolf方法[12]建立塔基與土體相互作用的耦合模型，采用模態截斷法建立風力機塔架模型，并基于混沌理論分析湍流風與地震聯合作用下塔架振動的非線性特征。對目標反應譜縮放后得150組不同強度地震運動?；陂_源軟件FAST[13]預留數據接口開發地震計算模塊，實現湍流風與不同強度地震共同作用仿真功能。得到塔架動態響應，為風力機塔架結構設計及安全評估提供參考。

1 研究對象及仿真模型

1.1 研究對象

以NREL 5 MW風力機為研究對象[14]，其主要結構及性能參數，如表1所示。

表1 NREL 5 MW風力機主要參數Tab.1 Configuration of NREL 5 MW wind turbine

1.2 土體-結構耦合模型

由于土體的非絕對剛性，在地震載荷作用過程中，土體與結構體之間既有力的相互作用，也有變形的相互制約[15]，因此可通過在土體與結構體之間設置彈簧和阻尼器模擬土體與結構體之間的應力與應變。圖1為風力機基礎平臺與土體耦合模型示意圖，對于剛性圓柱型基礎，根據土體特性和塔基尺寸確定剛度K及阻尼系數C，公式如式(1)～(2)所示：

圖1 風力機基礎平臺與土體耦合模型Fig.1 Interaction model for wind turbine platform and soil

(1)

(2)

式中：下標ζ表示方向；當ζ為x和y時，κ取值為8，γ取值為4.6，υs=2-μs，當ζ為z時，κ取值為4，γ取值為3.4，υs=1-μs；Gs、ρs和μs分別為土體的切變模量、密度和泊松比，硬黏土所對應值分別為32 MPa、2 700 kg/m3和0.25，Rs為基礎平臺的半徑，值為9 m。由此計算得水平方向的剛度和阻尼分別為1 316.6 MN/m和62.58 MN·s/m，豎直方向彈簧阻尼器的剛度和阻尼分別為1 536 MN/m和107.93 MN·s/m。

1.3 塔架動力學模型

塔架可視作質量與剛度連續分布的倒置懸臂梁，x和y方向位移相互獨立[16]。采用模態疊加法建立有限自由度N的塔架模型?？紤]主振型數目，t時刻塔架高度h處的位移u(h,t)如式(3)所示：

(3)

式中：ci(t)為與振型階數相關的廣義坐標，表示振型對應的塔架撓度；φi(h)為模態函數。

塔架的模態函數由每個模態下的形函數確定，對形函數進行線性組合得：

(4)

式中：Ti,j為比例常數，由模態階數和形函數確定。

當塔架以v階模態振動時，對應的坐標函數為qv(t)=Qvsin(ωvt+ψv)，Qv為超級單元體的振幅，ωv和ψv分別為v階模態所對應的頻率和相位。則v階形函數對應的廣義坐標如式(5)所示：

ci(t)=qv(t)Tv,i，i=1,2,…,N

(5)

根據塔架x和y方向的前兩階模態振型，將塔架簡化為具有4個自由度的梁單元，采用拉格朗日方程描述系統振動守恒方程：

(6)

式中：mij和kij分別為廣義質量與廣義剛度，定義分別如式(7)～式(8)所示：

(7)

(8)

式中：mtop為塔架頂部質量，為風輪、機艙及偏航裝置質量之和；EIT(h)為塔架高度h處的剛度；ρ(h)為塔架高度h處的密度；H為塔架高度。

將式(5)代入式(6)中，矩陣形式為：

(-ω2[M]+[K])[C]=[0]

(9)

式中：ω為塔架固有頻率；[M]為廣義質量矩陣；[K]為廣義剛度矩陣；[C]是由各階模態相關的比例常數組成的特征向量。

式(9)是關于ω的特征方程，解之可得到塔架自振頻率，塔架結構振動方程如式(10)所示，采用Runge-Kutta法求解可獲得塔架結構的時域位移及速度：

(10)

2 湍流風場與地震運動

2.1 湍流風場

通過Turbsim模擬風力機運行環境，假設湍流風場三個方向速度分量相互獨立[17]。輪轂高度處平均風速為風力機額定風速11.4 m/s。選用美國可再生能源實驗室(National Renewable Energy Laboratory，NREL)基于實測風場數據和SMOOTH風譜模型建立的NWTCUP風譜模型模擬湍流風場[18]。湍流風場風速時域變化，如圖2所示。

圖2 湍流風場Fig.2 Turbulence wind field

2.2 地震運動

通過匹配不同的響應譜得到150種地震加速度時域歷程以表征不同強度的地震運動。地震加速度采用匹配目標響應譜的方法計算[19]。

針對任意初始加速度時間序列a(t)和目標譜可通過圖3所示的方法進行匹配，圖4為地面加速度峰值(Peak Ground Acceleration，PGA)為0.3 g時來流方向加速度時域變化及目標譜匹配情況。tj時刻初始譜與目標譜之間的誤差為：

ΔRj=(Qj-Rj)

(11)

式中:Qj為目標譜值；Rj為初始譜值。

假設反應譜的峰值時間t不受小幅修正函數影響，則調整時間序列如式(12)所示：

(12)

式中：fi(t)為一組修正函數，定義如式(13)所示，bj為修正函數的幅值，N為需要匹配的樣本數。

fi(t)=hi(ti-t)

(13)

式中：hi(t)是單自由度的加速度脈沖響應函數，定義如式(14)所示。

(14)

加速度反應譜δRi為：

(15)

將式(12)代入式(15)中可得：

(16)

=[C]-1{δR}

(17)

計算式(17)可得修正函數幅值b，以此計算調整時間序列δa(t)。第一次匹配后，加速度時間序列為：

a1(t)=a(t)+γδa(t)

(18)

式中：γ為修正函數的松弛因子，γ∈(0,1)。

在第二次匹配中，a1(t)替代a(t)為新的時間序列，如此反復直至目標譜匹配達到目標精度為止。

圖3 目標譜匹配過程示意圖Fig.3 Flowchart of response spectrum match

PGA為0.3 g時所對應的偽譜加速度PSA(Pseudo-Spectral Acceleration)為0.127 g，匹配后PSA為0.128 g，誤差為0.8%，說明匹配后的加速度具有目標反應譜的特性[21]。

地震發生時，基礎平臺的目標加速度為地震加速度，此時基礎平臺地震載荷為：

(19)

3 塔架響應信號非線性分析方法

基于混沌理論對復雜工況下風力機動態響應進行非線性分析，采用相圖法及最大Lyapunov指數法對湍流風與地震聯合作用下的振動信號進行混沌特征識別。

(a) 地震加速度時域變化

(b) 地震加速度目標反應譜匹配情況圖4 地震來流方向加速度目標反應譜匹配情況及時域變化Fig.4 The acceleration of earthquake, target spectrum and matched spectrum in x direction

3.1 相圖法

相圖法是將風力機塔頂響應信號時間序列嵌入至三維相空間中進行相空間重構，相空間重構原理見文獻[22]。時間序列相圖可描述非線性系統隨時間的狀態變化，反映系統吸引子的空間結構。若系統相空間軌跡表現出在有限空間內不斷延伸和折疊的反復性永不相交的非周期運動狀態，即不是周期函數的重復性運動，又區別于毫無規律的隨機運動，則可定性判定系統具有混沌特征，且存在奇異吸引子[23-24]。即相圖法可依據風力機塔頂響應信號時間序列在三維相空間中表現出的特殊性，對塔頂響應信號混沌特征進行定性判定。

3.2 最大Lyapunov指數

非線性系統的最大Lyapunov指數是否大于0可作為系統振動信號是否具有混沌特征的定量指標[25]。其基本原理是：相空間重構后，相空間中相鄰兩條軌道隨時間逐漸發散或收斂，Lyapunov指數是其軌道發散或收斂的速率，用于刻畫系統動力學行為隨時間演化對初值的敏感程度。最大Lyapunov指數小于0表示時間序列具有隨機性或周期性。最大Lyapunov指數大于0，反應時間序列具有混沌特征，且值越大表明時間序列混沌特征越強。

4 結果分析

4.1 時域響應

圖5為湍流風與地震聯合作用及湍流風單獨作用時風力機塔頂位移動態響應。PGA為0.30 g和0.15 g分別對應設防烈度為7度和8度的工況。

圖5 塔頂位移動態響應對比Fig.5 Comparison for dynamic responses of tower top displacement on different conditions

從圖5看出，湍流風單獨作用時，塔頂位移在429.99 s達到峰值0.45 m，PGA為0.30 g和0.15 g時，地震與湍流風聯合作用，塔頂位移均在431.80 s達到峰值，分別為0.52 m和0.60 m，較湍流風單獨作用時分別增大15.6%和33.3%。在地震結束，塔頂位移可較快恢復到湍流風單獨作用狀態。

圖6為湍流風與地震聯合作用及湍流風單獨作用時風力機塔頂加速度動態響應。

圖6 塔頂加速度動態響應對比Fig.6 Comparison for dynamic responses of tower top acceleration on different conditions

從圖6可以看出，湍流風單獨作用時，塔頂加速度在438.37 s達到峰值0.31 m/s2。PGA為0.15 g和0.30 g的地震與湍流風聯合作用時，塔頂加速度分別在434.98 s和433.24 s時達到峰值0.66 m/s2和1.23 m/s2，分別增大了113%和297%。此外，PGA為0.30 g時，塔頂在地震激勵結束后仍有較大的加速度響應。表明地震作用對塔架加速度響應影響較大。

4.2 塔頂動態響應分析

不同強度地震下，風力機不同方向塔頂位移如圖7所示。

從圖7可以發現，來流方向塔頂存在約0.3 m初始位移，這是因為來流方向受湍流風影響較大。在415 s，塔架在不同方向均出現明顯的位移，且隨地震強度的增大而增大。425 s時，塔頂位移出現峰值，來流方向和側向的最大位移相當，說明地震強度較大時，地震載荷是塔架的主要載荷。同時，隨著地震過程的結束，來流方向的塔頂位移逐漸變小，而側向的塔頂位移峰值衰減程度較小，出現強度較大時間較長的振動過程，顯示來流方向塔架振動受氣動阻尼影響明顯，能量耗散較快。此外，來流方向的塔頂位移更易受到地震激勵的影響，在PGA為0.1 g時，來流方向的塔頂位移出現明顯變化，而側向變化較小，說明在小強度地震作用下，地震與湍流風耦合使來流方向塔頂處于不穩定狀態。

(a) 來流方向

(b) y方向圖7 不同強度地震下塔頂位移響應Fig.7 Tower top displacement at different seismic intensity in x and y directions

功率譜將原先對時間域的振動描述轉為頻率域的振動描述，從而得到各個頻率的振動能量分布。不同強度地震下，塔頂位移功率譜如圖8所示。

從圖8可以看出，塔頂來流方向和側向振動頻率特征較為相似，均在一階固有頻率0.32 Hz發生強迫共振。由于來流方向氣動阻尼較大，能量耗散較快，塔頂來流方向位移的功率譜密度較小。同時，來流方向塔頂位移在低頻能量較大，而側向塔頂位移能量主要集中在一階固有頻率附近，說明側向振動頻率主要為風力機結構一階固有頻率。

(a) 來流方向

(b) 側向圖8 不同強度地震下塔頂位移功率譜密度Fig.8 Tower top displacement spectrum power density at different seismic intensity in x and y directions

不同強度地震下，風力機不同方向塔頂加速度如圖9所示。

(a) 來流方向

(b) 側向圖9 不同強度地震下塔頂加速度Fig.9 Tower top acceleration at different seismic intensity in x and y directions

從圖9可以看出，不同方向的塔頂加速度受湍流風影響較小，說明塔頂振動主要受地震影響。428 s時，塔頂加速度出現峰值。在地震加速度衰減過程中，來流方向加速度逐漸減小，而側向加速度仍保持一定幅度。這說明塔頂來流方向振動受湍流風影響較大，地震激勵引起的振動能量可通過氣動阻尼耗散，而塔頂側向振動能量只能通過結構阻尼耗散。

不同強度地震下，塔頂加速度功率譜如圖10所示。

從圖10可以看出，不同方向的塔頂加速度能量密度峰值均出現在一階固有頻率0.32 Hz附近，來流方向峰值相比側向較小。來流方向能量主要分布在0.15～0.5 Hz區間內，在風輪轉動頻率0.2 Hz處有小峰值，而側向能量集中在0.25～0.4 Hz區間內，一階固有頻率對應的峰值遠大于其他頻率下的能量密度。說明來流方向塔頂振動受風輪轉動影響，而塔頂側向振動特性主要由風力機自身結構確定。

4.3 塔頂振動信號非線性分析

通過C-C算法[26]計算得不同塔架高度加速度的最佳延遲時間τ，塔頂響應信號嵌入至三維相空間中。PGA分別為0、0.15 g和0.30 g時重構后塔頂位移三維相空間如圖10所示。

從圖11可以發現，在不同工況下，風力機塔頂位移相圖軌跡均呈現出圍繞一個或多個中心反復纏繞的現象，表明存在奇異吸引子，說明其塔頂位移信號具有混沌特征。湍流風單獨作用時，大部分塔架位移軌跡被限制在一定區域內，少數形成發散的尾跡，說明塔頂在湍流風的作用下在0.3 m附近小幅振動。而湍流風與地震聯合作用時，塔頂位移軌跡均形成圍繞一個中心點的圓環狀，混沌特征最為明顯，說明塔架處于極度不穩定的運動狀態中。同時，相圖也反應出塔頂動力學特性，隨著地震強度的增大，塔頂位移也會隨之增大。

(a) 來流方向

(b) 側向圖10 不同強度地震下塔頂加速度功率譜Fig.10 Tower top acceleration spectrum power density at different seismic intensity in x and y directions

(a) PGA = 0

(b) PGA = 0.15 g

(c) PGA = 0.30 g圖11 不同地震強度下的塔頂位移時間序列三維相圖Fig.11 The 3-D phase space of tower top displacement time series at different seismic intensity

為了更準確地研究塔頂位移信號的混沌特征，計算了不同工況下塔頂位移的最大Lyapunov指數，其結果如表2所示。由表可知，3種不同工況時的塔頂位移響應信號時間序列的最大Lyapunov指數各不相同，但均大于0，體現了塔頂位移信號均表現出不同程度的混沌特征。其中，PGA為0.15 g的地震與湍流風聯合作用時塔頂位移響應信號時間序列最大Lyapunov指數最大，為0.021 4，表明其混沌特征最明顯；而湍流風單獨作用時的最大Lyapunov指數最小，為0.001 2，說明其混沌特征最弱。

表2 塔頂位移信號最大Lyapunov指數Tab.2 The max Lyapunov index of tower top displacement signal at different seismic intensity

5 結 論

以NREL 5 MW風力機為研究對象，通過模態截斷法理論建立塔架模型，根據Wolf理論考慮塔基結構與土體之間的相互作用。計算了150組不同強度地震與湍流風聯合作用下風力機動力學響應，并基于混沌理論，采用相圖法及最大Lyapunov指數法分析塔頂響應非線性特征，主要得到以下結論：

(1) 地震激勵對塔頂加速度影響較大，PGA為0.15 g和0.30 g時，塔頂位移較無地震工況分別增大了15.6%和33.3%，塔頂加速度較無地震工況分別增大了113%和297%。

(2) 由于湍流風作用，來流方向塔頂存在約0.3 m初始位移。隨著地震加速度的衰減，來流方向的塔頂位移和加速度逐漸變小，而側向的塔頂位移和加速度仍保持一定幅值。塔頂來流方向振動受湍流風影響較大，地震激勵引起的振動能量可通過氣動阻尼耗散，而塔頂側向振動能量只能通過結構阻尼耗散。來流方向塔頂振動受風輪轉動影響，而側向塔頂振動頻率主要為風力機結構一階固有頻率。

(3) 湍流風單獨作用或與地震聯合作用時，塔頂位移信號均具有混沌特性。湍流風與地震耦合會加劇塔架振動的復雜程度和不穩定性，塔頂位移的最大Lyapunov指數較大，混沌特征較為明顯，而湍流風單獨作用時塔頂位移的最大Lyapunov指數較小，混沌特征相對較弱。