鄒錦華, 楊 陽, 李 春,2, 劉中勝, 袁全勇
(1. 上海理工大學(xué) 能源與動力工程學(xué)院, 上海 200093; 2. 上海市動力工程多相流動與傳熱重點實驗室, 上海 200093)
“十二五”規(guī)劃期間,我國風(fēng)電年均增長34.6%,為第二大非化石能源,2016年新增裝機量為23 370 MW,總裝機量達(dá)168 732 MW?!笆濉睍r期計劃將風(fēng)電占全國總發(fā)電量的比例提高至6%以實現(xiàn)非化石能源占一次能源消費比重達(dá)15%的重大目標(biāo)[1-2]。至2020年,我國中東部和南方地區(qū)陸上風(fēng)電新增并網(wǎng)裝機容量為42 000 MW以上,西北地區(qū)累計并網(wǎng)容量將達(dá)到48 500 MW[3]。上述區(qū)域分別處于環(huán)太平洋地震帶及亞歐地震帶附近,風(fēng)電場極有可能受到地震影響。大型風(fēng)力機塔架力學(xué)結(jié)構(gòu)屬于典型的細(xì)長柔性體,加之塔架頂部的大質(zhì)量機艙,在空間非均勻和時間非定常的湍流風(fēng)作用下,塔架氣彈響應(yīng)具有明顯的非線性和非平穩(wěn)特征[4-5],若同時受到高強度地震作用,可能發(fā)生屈曲及結(jié)構(gòu)永久變形,影響結(jié)構(gòu)安全[6]。因此,開展湍流風(fēng)與地震聯(lián)合作用下風(fēng)力機塔架動力學(xué)響應(yīng)研究對風(fēng)力機設(shè)計及安全評估具有重要意義。
關(guān)于地震作用對風(fēng)力機塔架結(jié)構(gòu)的影響,國內(nèi)外學(xué)者開展了相關(guān)研究。Adam等[7]以1.5 MW風(fēng)力機為研究對象,建立塔架有限元模型,分別研究了10種近斷層地震和遠(yuǎn)場地震作用下塔架模態(tài)振型及在塑性變形時地震能量耗散情況,以5%的氣動阻尼代替湍流風(fēng)作用,發(fā)現(xiàn)風(fēng)力機主要受前兩階基礎(chǔ)振型的影響,且近斷層地震對塔架各段接縫處的破壞影響較大。季亮等[8]建立風(fēng)力機單自由度和雙自由度簡化模型,運用剪力法計算8級地震作用下不同功率風(fēng)力機的塔架底部彎矩與剪切力,因研究參考建筑結(jié)構(gòu)設(shè)計方法簡化風(fēng)力機模型,未能考慮風(fēng)輪所受的非定常氣動載荷。Kim等[9]建立近海風(fēng)力機樁-土體耦合模型,分析了土質(zhì)對風(fēng)力機地震動力學(xué)響應(yīng)的影響,計算模型中忽略了氣動載荷的影響。Ikwulono等[10]采用故障樹方法分析風(fēng)與地震分別獨立及聯(lián)合作用時風(fēng)力機動力學(xué)響應(yīng),運用有限元軟件分析風(fēng)載荷與地震載荷共同作用時塔基失效概率。李春等[11]基于開源風(fēng)力機仿真軟件FAST建立湍流風(fēng)-地震實時耦合模型,模擬多組湍流風(fēng)況與地震共同作用風(fēng)力機動態(tài)響應(yīng),發(fā)現(xiàn)地震激勵下風(fēng)力機塔基載荷達(dá)到穩(wěn)態(tài)風(fēng)工況的2倍~15倍。
風(fēng)力機的運行受地理環(huán)境、氣候條件及表面粗糙度等諸多因素的影響,加之湍流風(fēng)與地震聯(lián)合激勵作用,導(dǎo)致風(fēng)力機動態(tài)響應(yīng)具有明顯非線性特征,因此,采用混沌理論這一非線性分析工具來更加準(zhǔn)確揭示湍流風(fēng)與地震聯(lián)合作用等復(fù)雜工況下塔架振動具有重要意義。
為研究風(fēng)力機塔架結(jié)構(gòu)在湍流風(fēng)與地震聯(lián)合作用等不利情況下振動的非線性特性,以NREL 5 MW風(fēng)力機為研究對象,根據(jù)Wolf方法[12]建立塔基與土體相互作用的耦合模型,采用模態(tài)截斷法建立風(fēng)力機塔架模型,并基于混沌理論分析湍流風(fēng)與地震聯(lián)合作用下塔架振動的非線性特征。對目標(biāo)反應(yīng)譜縮放后得150組不同強度地震運動?;陂_源軟件FAST[13]預(yù)留數(shù)據(jù)接口開發(fā)地震計算模塊,實現(xiàn)湍流風(fēng)與不同強度地震共同作用仿真功能。得到塔架動態(tài)響應(yīng),為風(fēng)力機塔架結(jié)構(gòu)設(shè)計及安全評估提供參考。
以NREL 5 MW風(fēng)力機為研究對象[14],其主要結(jié)構(gòu)及性能參數(shù),如表1所示。
表1 NREL 5 MW風(fēng)力機主要參數(shù)Tab.1 Configuration of NREL 5 MW wind turbine
由于土體的非絕對剛性,在地震載荷作用過程中,土體與結(jié)構(gòu)體之間既有力的相互作用,也有變形的相互制約[15],因此可通過在土體與結(jié)構(gòu)體之間設(shè)置彈簧和阻尼器模擬土體與結(jié)構(gòu)體之間的應(yīng)力與應(yīng)變。圖1為風(fēng)力機基礎(chǔ)平臺與土體耦合模型示意圖,對于剛性圓柱型基礎(chǔ),根據(jù)土體特性和塔基尺寸確定剛度K及阻尼系數(shù)C,公式如式(1)~(2)所示:
圖1 風(fēng)力機基礎(chǔ)平臺與土體耦合模型Fig.1 Interaction model for wind turbine platform and soil
(1)
(2)
式中:下標(biāo)ζ表示方向;當(dāng)ζ為x和y時,κ取值為8,γ取值為4.6,υs=2-μs,當(dāng)ζ為z時,κ取值為4,γ取值為3.4,υs=1-μs;Gs、ρs和μs分別為土體的切變模量、密度和泊松比,硬黏土所對應(yīng)值分別為32 MPa、2 700 kg/m3和0.25,Rs為基礎(chǔ)平臺的半徑,值為9 m。由此計算得水平方向的剛度和阻尼分別為1 316.6 MN/m和62.58 MN·s/m,豎直方向彈簧阻尼器的剛度和阻尼分別為1 536 MN/m和107.93 MN·s/m。
塔架可視作質(zhì)量與剛度連續(xù)分布的倒置懸臂梁,x和y方向位移相互獨立[16]。采用模態(tài)疊加法建立有限自由度N的塔架模型??紤]主振型數(shù)目,t時刻塔架高度h處的位移u(h,t)如式(3)所示:
(3)
式中:ci(t)為與振型階數(shù)相關(guān)的廣義坐標(biāo),表示振型對應(yīng)的塔架撓度;φi(h)為模態(tài)函數(shù)。
塔架的模態(tài)函數(shù)由每個模態(tài)下的形函數(shù)確定,對形函數(shù)進行線性組合得:
(4)
式中:Ti,j為比例常數(shù),由模態(tài)階數(shù)和形函數(shù)確定。
當(dāng)塔架以v階模態(tài)振動時,對應(yīng)的坐標(biāo)函數(shù)為qv(t)=Qvsin(ωvt+ψv),Qv為超級單元體的振幅,ωv和ψv分別為v階模態(tài)所對應(yīng)的頻率和相位。則v階形函數(shù)對應(yīng)的廣義坐標(biāo)如式(5)所示:
ci(t)=qv(t)Tv,i,i=1,2,…,N
(5)
根據(jù)塔架x和y方向的前兩階模態(tài)振型,將塔架簡化為具有4個自由度的梁單元,采用拉格朗日方程描述系統(tǒng)振動守恒方程:
(6)
式中:mij和kij分別為廣義質(zhì)量與廣義剛度,定義分別如式(7)~式(8)所示:
(7)
(8)
式中:mtop為塔架頂部質(zhì)量,為風(fēng)輪、機艙及偏航裝置質(zhì)量之和;EIT(h)為塔架高度h處的剛度;ρ(h)為塔架高度h處的密度;H為塔架高度。
將式(5)代入式(6)中,矩陣形式為:
(-ω2[M]+[K])[C]=[0]
(9)
式中:ω為塔架固有頻率;[M]為廣義質(zhì)量矩陣;[K]為廣義剛度矩陣;[C]是由各階模態(tài)相關(guān)的比例常數(shù)組成的特征向量。
式(9)是關(guān)于ω的特征方程,解之可得到塔架自振頻率,塔架結(jié)構(gòu)振動方程如式(10)所示,采用Runge-Kutta法求解可獲得塔架結(jié)構(gòu)的時域位移及速度:
(10)
通過Turbsim模擬風(fēng)力機運行環(huán)境,假設(shè)湍流風(fēng)場三個方向速度分量相互獨立[17]。輪轂高度處平均風(fēng)速為風(fēng)力機額定風(fēng)速11.4 m/s。選用美國可再生能源實驗室(National Renewable Energy Laboratory,NREL)基于實測風(fēng)場數(shù)據(jù)和SMOOTH風(fēng)譜模型建立的NWTCUP風(fēng)譜模型模擬湍流風(fēng)場[18]。湍流風(fēng)場風(fēng)速時域變化,如圖2所示。
圖2 湍流風(fēng)場Fig.2 Turbulence wind field
通過匹配不同的響應(yīng)譜得到150種地震加速度時域歷程以表征不同強度的地震運動。地震加速度采用匹配目標(biāo)響應(yīng)譜的方法計算[19]。
針對任意初始加速度時間序列a(t)和目標(biāo)譜可通過圖3所示的方法進行匹配,圖4為地面加速度峰值(Peak Ground Acceleration,PGA)為0.3 g時來流方向加速度時域變化及目標(biāo)譜匹配情況。tj時刻初始譜與目標(biāo)譜之間的誤差為:
ΔRj=(Qj-Rj)
(11)
式中:Qj為目標(biāo)譜值;Rj為初始譜值。
假設(shè)反應(yīng)譜的峰值時間t不受小幅修正函數(shù)影響,則調(diào)整時間序列如式(12)所示:
(12)
式中:fi(t)為一組修正函數(shù),定義如式(13)所示,bj為修正函數(shù)的幅值,N為需要匹配的樣本數(shù)。
fi(t)=hi(ti-t)
(13)
式中:hi(t)是單自由度的加速度脈沖響應(yīng)函數(shù),定義如式(14)所示。
(14)
加速度反應(yīng)譜δRi為:
(15)
將式(12)代入式(15)中可得:
(16)
=[C]-1{δR}
(17)
計算式(17)可得修正函數(shù)幅值b,以此計算調(diào)整時間序列δa(t)。第一次匹配后,加速度時間序列為:
a1(t)=a(t)+γδa(t)
(18)
式中:γ為修正函數(shù)的松弛因子,γ∈(0,1)。
在第二次匹配中,a1(t)替代a(t)為新的時間序列,如此反復(fù)直至目標(biāo)譜匹配達(dá)到目標(biāo)精度為止。
圖3 目標(biāo)譜匹配過程示意圖Fig.3 Flowchart of response spectrum match
PGA為0.3 g時所對應(yīng)的偽譜加速度PSA(Pseudo-Spectral Acceleration)為0.127 g,匹配后PSA為0.128 g,誤差為0.8%,說明匹配后的加速度具有目標(biāo)反應(yīng)譜的特性[21]。
地震發(fā)生時,基礎(chǔ)平臺的目標(biāo)加速度為地震加速度,此時基礎(chǔ)平臺地震載荷為:
(19)
基于混沌理論對復(fù)雜工況下風(fēng)力機動態(tài)響應(yīng)進行非線性分析,采用相圖法及最大Lyapunov指數(shù)法對湍流風(fēng)與地震聯(lián)合作用下的振動信號進行混沌特征識別。
(a) 地震加速度時域變化
(b) 地震加速度目標(biāo)反應(yīng)譜匹配情況圖4 地震來流方向加速度目標(biāo)反應(yīng)譜匹配情況及時域變化Fig.4 The acceleration of earthquake, target spectrum and matched spectrum in x direction
相圖法是將風(fēng)力機塔頂響應(yīng)信號時間序列嵌入至三維相空間中進行相空間重構(gòu),相空間重構(gòu)原理見文獻(xiàn)[22]。時間序列相圖可描述非線性系統(tǒng)隨時間的狀態(tài)變化,反映系統(tǒng)吸引子的空間結(jié)構(gòu)。若系統(tǒng)相空間軌跡表現(xiàn)出在有限空間內(nèi)不斷延伸和折疊的反復(fù)性永不相交的非周期運動狀態(tài),即不是周期函數(shù)的重復(fù)性運動,又區(qū)別于毫無規(guī)律的隨機運動,則可定性判定系統(tǒng)具有混沌特征,且存在奇異吸引子[23-24]。即相圖法可依據(jù)風(fēng)力機塔頂響應(yīng)信號時間序列在三維相空間中表現(xiàn)出的特殊性,對塔頂響應(yīng)信號混沌特征進行定性判定。
非線性系統(tǒng)的最大Lyapunov指數(shù)是否大于0可作為系統(tǒng)振動信號是否具有混沌特征的定量指標(biāo)[25]。其基本原理是:相空間重構(gòu)后,相空間中相鄰兩條軌道隨時間逐漸發(fā)散或收斂,Lyapunov指數(shù)是其軌道發(fā)散或收斂的速率,用于刻畫系統(tǒng)動力學(xué)行為隨時間演化對初值的敏感程度。最大Lyapunov指數(shù)小于0表示時間序列具有隨機性或周期性。最大Lyapunov指數(shù)大于0,反應(yīng)時間序列具有混沌特征,且值越大表明時間序列混沌特征越強。
圖5為湍流風(fēng)與地震聯(lián)合作用及湍流風(fēng)單獨作用時風(fēng)力機塔頂位移動態(tài)響應(yīng)。PGA為0.30 g和0.15 g分別對應(yīng)設(shè)防烈度為7度和8度的工況。
圖5 塔頂位移動態(tài)響應(yīng)對比Fig.5 Comparison for dynamic responses of tower top displacement on different conditions
從圖5看出,湍流風(fēng)單獨作用時,塔頂位移在429.99 s達(dá)到峰值0.45 m,PGA為0.30 g和0.15 g時,地震與湍流風(fēng)聯(lián)合作用,塔頂位移均在431.80 s達(dá)到峰值,分別為0.52 m和0.60 m,較湍流風(fēng)單獨作用時分別增大15.6%和33.3%。在地震結(jié)束,塔頂位移可較快恢復(fù)到湍流風(fēng)單獨作用狀態(tài)。
圖6為湍流風(fēng)與地震聯(lián)合作用及湍流風(fēng)單獨作用時風(fēng)力機塔頂加速度動態(tài)響應(yīng)。
圖6 塔頂加速度動態(tài)響應(yīng)對比Fig.6 Comparison for dynamic responses of tower top acceleration on different conditions
從圖6可以看出,湍流風(fēng)單獨作用時,塔頂加速度在438.37 s達(dá)到峰值0.31 m/s2。PGA為0.15 g和0.30 g的地震與湍流風(fēng)聯(lián)合作用時,塔頂加速度分別在434.98 s和433.24 s時達(dá)到峰值0.66 m/s2和1.23 m/s2,分別增大了113%和297%。此外,PGA為0.30 g時,塔頂在地震激勵結(jié)束后仍有較大的加速度響應(yīng)。表明地震作用對塔架加速度響應(yīng)影響較大。
不同強度地震下,風(fēng)力機不同方向塔頂位移如圖7所示。
從圖7可以發(fā)現(xiàn),來流方向塔頂存在約0.3 m初始位移,這是因為來流方向受湍流風(fēng)影響較大。在415 s,塔架在不同方向均出現(xiàn)明顯的位移,且隨地震強度的增大而增大。425 s時,塔頂位移出現(xiàn)峰值,來流方向和側(cè)向的最大位移相當(dāng),說明地震強度較大時,地震載荷是塔架的主要載荷。同時,隨著地震過程的結(jié)束,來流方向的塔頂位移逐漸變小,而側(cè)向的塔頂位移峰值衰減程度較小,出現(xiàn)強度較大時間較長的振動過程,顯示來流方向塔架振動受氣動阻尼影響明顯,能量耗散較快。此外,來流方向的塔頂位移更易受到地震激勵的影響,在PGA為0.1 g時,來流方向的塔頂位移出現(xiàn)明顯變化,而側(cè)向變化較小,說明在小強度地震作用下,地震與湍流風(fēng)耦合使來流方向塔頂處于不穩(wěn)定狀態(tài)。
(a) 來流方向
(b) y方向圖7 不同強度地震下塔頂位移響應(yīng)Fig.7 Tower top displacement at different seismic intensity in x and y directions
功率譜將原先對時間域的振動描述轉(zhuǎn)為頻率域的振動描述,從而得到各個頻率的振動能量分布。不同強度地震下,塔頂位移功率譜如圖8所示。
從圖8可以看出,塔頂來流方向和側(cè)向振動頻率特征較為相似,均在一階固有頻率0.32 Hz發(fā)生強迫共振。由于來流方向氣動阻尼較大,能量耗散較快,塔頂來流方向位移的功率譜密度較小。同時,來流方向塔頂位移在低頻能量較大,而側(cè)向塔頂位移能量主要集中在一階固有頻率附近,說明側(cè)向振動頻率主要為風(fēng)力機結(jié)構(gòu)一階固有頻率。
(a) 來流方向
(b) 側(cè)向圖8 不同強度地震下塔頂位移功率譜密度Fig.8 Tower top displacement spectrum power density at different seismic intensity in x and y directions
不同強度地震下,風(fēng)力機不同方向塔頂加速度如圖9所示。
(a) 來流方向
(b) 側(cè)向圖9 不同強度地震下塔頂加速度Fig.9 Tower top acceleration at different seismic intensity in x and y directions
從圖9可以看出,不同方向的塔頂加速度受湍流風(fēng)影響較小,說明塔頂振動主要受地震影響。428 s時,塔頂加速度出現(xiàn)峰值。在地震加速度衰減過程中,來流方向加速度逐漸減小,而側(cè)向加速度仍保持一定幅度。這說明塔頂來流方向振動受湍流風(fēng)影響較大,地震激勵引起的振動能量可通過氣動阻尼耗散,而塔頂側(cè)向振動能量只能通過結(jié)構(gòu)阻尼耗散。
不同強度地震下,塔頂加速度功率譜如圖10所示。
從圖10可以看出,不同方向的塔頂加速度能量密度峰值均出現(xiàn)在一階固有頻率0.32 Hz附近,來流方向峰值相比側(cè)向較小。來流方向能量主要分布在0.15~0.5 Hz區(qū)間內(nèi),在風(fēng)輪轉(zhuǎn)動頻率0.2 Hz處有小峰值,而側(cè)向能量集中在0.25~0.4 Hz區(qū)間內(nèi),一階固有頻率對應(yīng)的峰值遠(yuǎn)大于其他頻率下的能量密度。說明來流方向塔頂振動受風(fēng)輪轉(zhuǎn)動影響,而塔頂側(cè)向振動特性主要由風(fēng)力機自身結(jié)構(gòu)確定。
通過C-C算法[26]計算得不同塔架高度加速度的最佳延遲時間τ,塔頂響應(yīng)信號嵌入至三維相空間中。PGA分別為0、0.15 g和0.30 g時重構(gòu)后塔頂位移三維相空間如圖10所示。
從圖11可以發(fā)現(xiàn),在不同工況下,風(fēng)力機塔頂位移相圖軌跡均呈現(xiàn)出圍繞一個或多個中心反復(fù)纏繞的現(xiàn)象,表明存在奇異吸引子,說明其塔頂位移信號具有混沌特征。湍流風(fēng)單獨作用時,大部分塔架位移軌跡被限制在一定區(qū)域內(nèi),少數(shù)形成發(fā)散的尾跡,說明塔頂在湍流風(fēng)的作用下在0.3 m附近小幅振動。而湍流風(fēng)與地震聯(lián)合作用時,塔頂位移軌跡均形成圍繞一個中心點的圓環(huán)狀,混沌特征最為明顯,說明塔架處于極度不穩(wěn)定的運動狀態(tài)中。同時,相圖也反應(yīng)出塔頂動力學(xué)特性,隨著地震強度的增大,塔頂位移也會隨之增大。
(a) 來流方向
(b) 側(cè)向圖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
為了更準(zhǔn)確地研究塔頂位移信號的混沌特征,計算了不同工況下塔頂位移的最大Lyapunov指數(shù),其結(jié)果如表2所示。由表可知,3種不同工況時的塔頂位移響應(yīng)信號時間序列的最大Lyapunov指數(shù)各不相同,但均大于0,體現(xiàn)了塔頂位移信號均表現(xiàn)出不同程度的混沌特征。其中,PGA為0.15 g的地震與湍流風(fēng)聯(lián)合作用時塔頂位移響應(yīng)信號時間序列最大Lyapunov指數(shù)最大,為0.021 4,表明其混沌特征最明顯;而湍流風(fēng)單獨作用時的最大Lyapunov指數(shù)最小,為0.001 2,說明其混沌特征最弱。
表2 塔頂位移信號最大Lyapunov指數(shù)Tab.2 The max Lyapunov index of tower top displacement signal at different seismic intensity
以NREL 5 MW風(fēng)力機為研究對象,通過模態(tài)截斷法理論建立塔架模型,根據(jù)Wolf理論考慮塔基結(jié)構(gòu)與土體之間的相互作用。計算了150組不同強度地震與湍流風(fēng)聯(lián)合作用下風(fēng)力機動力學(xué)響應(yīng),并基于混沌理論,采用相圖法及最大Lyapunov指數(shù)法分析塔頂響應(yīng)非線性特征,主要得到以下結(jié)論:
(1) 地震激勵對塔頂加速度影響較大,PGA為0.15 g和0.30 g時,塔頂位移較無地震工況分別增大了15.6%和33.3%,塔頂加速度較無地震工況分別增大了113%和297%。
(2) 由于湍流風(fēng)作用,來流方向塔頂存在約0.3 m初始位移。隨著地震加速度的衰減,來流方向的塔頂位移和加速度逐漸變小,而側(cè)向的塔頂位移和加速度仍保持一定幅值。塔頂來流方向振動受湍流風(fēng)影響較大,地震激勵引起的振動能量可通過氣動阻尼耗散,而塔頂側(cè)向振動能量只能通過結(jié)構(gòu)阻尼耗散。來流方向塔頂振動受風(fēng)輪轉(zhuǎn)動影響,而側(cè)向塔頂振動頻率主要為風(fēng)力機結(jié)構(gòu)一階固有頻率。
(3) 湍流風(fēng)單獨作用或與地震聯(lián)合作用時,塔頂位移信號均具有混沌特性。湍流風(fēng)與地震耦合會加劇塔架振動的復(fù)雜程度和不穩(wěn)定性,塔頂位移的最大Lyapunov指數(shù)較大,混沌特征較為明顯,而湍流風(fēng)單獨作用時塔頂位移的最大Lyapunov指數(shù)較小,混沌特征相對較弱。