臺風過境下大型單樁式海上風機結構動力特性研究1)

秦夢飛 施 偉,?,2) 柴 威 付 興 李 昕

* (大連理工大學建設工程學部,遼寧大連 116024)

? (大連理工大學海岸與近海工程國家重點實驗室,遼寧大連 116024)

** (武漢理工大學船海與能源動力工程學院,武漢 430063)

引言

海上風電作為最具商業化前景的一種可再生能源,近年來得到世界各國的高度重視,呈現出跨越式發展.2020 年我國新增海上風電裝機規模306 萬千瓦,占據全球新增容量的50%[1].我國海岸線長度超過18 000 km,可開發利用的風能資源豐富,且我國經濟發達城市集中在沿海地區,發展海上風電可就近供給,緩解能源供應壓力并助力碳中和.然而,我國東南沿海處于臺風易發區域,以廣東為例,近60 年的統計資料表明對廣東省沿海有影響的年均臺風次數達到了5 次,其中超強臺風占到了年均總數的49%[2].臺風頻繁經過極大地影響了海上風電結構的安全性:一方面,大型風機的塔筒與葉片屬于高柔性結構,在臺風作用下易產生變形與振動,具有極強的風致敏感性;另一方面,臺風會引起巨大的臺風浪,對海上風機基礎造成災害[3].目前,全球并網的海上風機,80% 采用單樁式基礎.因此,開展臺風與臺風浪聯合作用下大型單樁海上風機響應研究具有重要的工程意義.

對臺風過境期間風場的實測數據表明,臺風與常態大風具有不同的湍流特性[4-7].Li等[8]依據臺風黑格比的實測數據,擬合出臺風不同區域的風譜,結果表明相對于常規風場,臺風區歸一化風譜的譜峰頻率向高頻轉移.美國國家可再生能源實驗室(NREL)[9]比較了風機分析標準中所用的多種湍流譜模型,發現即使平均風速與湍流度相同,采用不同風譜模型將導致結構響應產生較大差異.韓然等[10]學者將臺風風場劃分為渦旋區與大風區,比較了中國抗風設計指南中推薦的Simiu 譜與實測風譜下風力機動力響應結果,結果顯示在臺風大風區實測風場與規范推薦譜得到的風場脈動特性類似,但渦旋區差異則較大.Wang等[11]將臺風過境全程分為五個區域,對過境全程風力機載荷效應進行了數值模擬,結果表明臺風效應極大地增加了風力機的風載荷極值,且風載荷的增加呈現出非線性特性.李琪等[12]使用有限元方法分析了導管架與單樁基礎在臺風經過時的結構響應,得出了導管架基礎風力機動力響應較小的結論.Kim[13]使用氣象海浪耦合模式University of Miami Coupled Model (UMCM)模擬臺風路徑,得到風浪耦合下導管架基礎風力機動力響應,并在模擬臺風場時使用了IEC 規范中的Kaimal 譜.朱彬彬等[14]建立了臺風作用下運營期內海上風機基礎結構失效模式及評價指標.蘆直躍等[15]通過物理模型實驗,研究了臺風對海上風電單樁基礎累積變形的影響,結果指出可采用疊加法計算存在臺風時單樁基礎的累積轉角.

在風電基礎結構動力特性方面,姜貞強等[16]使用ABAQUS 有限元軟件,對單樁基礎在我國黃海50 年一遇的極端載荷作用下的動力響應進行了研究,結果表明增大樁徑及壁厚可有效增大單樁剛度,減小泥面線處動力響應.孫肖菲等[17]使用有限元對大直徑單樁基礎固有頻率影響因素進行了研究,研究結果表明基于梁單元及非線性彈簧單元體系計算得出的一二階固有頻率明顯小于實體單元所得結果.劉晨晨等[18]使用ABAQUS 研究了地震與波浪共同作用下樁的動力響應,結果指出相較于地震載荷,波浪載荷影響較小.Bergua等[19]基于國際能源署OC6 項目,采用不同數值模擬軟件,對10 MW級大型單樁結構的高精度樁土作用模型進行了對比驗證.Bachynski等[20]在SIMA 中使用宏元素考慮樁土相互作用,比較了等效風暴法中不同風暴面對結構響應的影響,其結果表明風暴中最大響應不一定出現在有義波高最大的區域.

國內外針對臺風對海上風電結構的影響,已做了大量的研究工作.然而,上述研究僅計入了臺風影響,忽略了臺風浪的作用,或在模擬風速場時使用了規范推薦譜.本文結合臺風期間實測風譜,開展風浪聯合作用下10 MW 級大型單樁海上風機動力響應研究,重點關注整個臺風過程中,不同臺風階段海上風機的動力響應特性.研究發現,在臺風經過的不同階段,單樁海上風機結構表現出不同的動力特性.

1 計算理論

1.1 風載荷計算理論

風機葉片載荷計算多采用葉素動量理論.葉素動量理論[21]結合了葉素理論與動量理論,將葉片沿徑向分成多段,將每段視為獨立的二維翼型單元[22],假設各單元之間流場互不影響,每個單元中心處相對風速與雷諾數分別為

其中 ρair為空氣密度,c為葉片 截面弦長,νair為空 氣黏度,uwind為來流風速,ufoil為葉素單元速度,ur為相對速度,Vrx及Vry為相對速度在葉片坐標系下的分量.

基于雷諾數以及相對風速攻角,在翼形庫中選取升力系數CL、阻力系數CD、力矩系數CM,可計算出各單元所受風載荷

沿徑向積分即可得到葉片所受氣動載荷.

水面上的塔筒及樁基礎同樣受風載荷作用,單位長度風力為

對于圓柱形結構,阻力系數CD取1.2,D為結構外徑,V為相對風速大小.

1.2 波浪載荷計算理論

本文使用Morison 方程[23]計算樁基礎所受波浪載荷.假定樁基礎的存在不影響波浪運動,則樁基所受波浪力由慣性力及速度力構成

式中,ρ 為海水密度,取1 025 kg/m3,Cm和CD分別為無量綱質量力系數與拖曳力系數.ar為流體質點相對樁基礎的加速度,ur為流體質點與結構的相對速度

式中,下標w與s分別代表流體質點與結構點,流體質點速度與加速度由二階斯托克斯波浪理論計算得到.不規則波的一階波面升高及一階速度勢根據波浪譜由線性波疊加得到

式中,η與 ? 分別表示波面升高與速度勢,上標代表攝動展開階數,A為組成波幅,ε為組成波隨機相位.不規則波中的二階波面升高及二階速度勢在一階波面升高及一階速度勢的約束下得到,將一階速度勢代入二階自由表面條件中,求解拉普拉斯方程,可得二階速度勢

式中,上標 ± 表示分別取+和 -后求和,G與D分別為

進一步可求出二階波面升高,將一階結果與二階結果相加得到總波面升高及速度勢,進而由速度勢得到流體質點速度與加速度.

為保證二階波浪的有效性,二階規則波需滿足3 點:(1)二階項遠小于一階項;(2)波谷不會觸底;(3)波浪不發生破碎.對于二階不規則波浪,Hu 和Zhao[24]在大量數值模擬的基礎上,提出當Hs/Lz小于0.08時,可認為二階不規則波是有效的.其中Hs為有效波高,Lz為零交叉頻率對應的波長.本文據此準則選擇波譜的截止頻率,對最惡劣的工況FEWS 階段,當截止頻率為3 rad/s 時,Hs/Lz約為0.072,此時忽略了約0.13%的總能量.

1.3 臺風風場模擬

湍流風場可以看作不同時空尺度的漩渦疊加,湍流譜描述了風場中能量的分布,根據柯爾莫哥洛夫假設,慣性子區內的風譜為[8]

式中,au為常數,取0.5;K為波數;ε為能量傳輸率,在表面層中為

其中,k為卡爾曼常數,u*為摩擦速度,U為平均速度,?ε是風切變的無量綱莫寧-奧布霍夫函數,在中性大氣條件下取1.

結合式(12)與式(13),湍流風譜為

式中,Au為無量綱常數,通常等于0.27,n為頻率,f為經高度z與平均速度u無量綱化后的頻率,f=nz/u.現存的風譜模型大都基于式(14) 得到,結合Wang等[11]提出的方法,可將臺風風場分為五個區域:臺風外圍渦旋前緣FOVS (front outer vortex stage)、前眼壁強風區FEWS (front eye wall stage)、臺風眼TES (typhoon eye stage)、后眼壁強風BEWS(back eye wall stage)、外圍渦旋后緣(back outer vortex stage).將臺風經過各階段的縱向風功率譜考慮如式(15),其中對風速平緩的臺風眼區域,使用了IEC 規范中推薦的Kaimal 譜,其余區域風譜則使用Li等[8]基于2008 年14 號臺風黑格比在南海的實測數據提出的湍流譜

式中,σ 為湍流風速標準差,Vhub為輪轂處風速,L為紊流尺度,現有規范中的風譜模型忽略了縱向風譜與橫向風譜的相關性,縱向風譜與橫向風譜各自基于實測數據由式(14)擬合得到[25-28].馮·卡門基于各向同性假設,提出橫向風譜與縱向風譜應有以下關系

式中,系數0.5 是由各向同性假設所決定的,實際臺風過程中湍流風場并非各向同性的,因此,需對式(16)進行修改.Tao等[29]提出可用參數km代替系數0.5 以考慮湍流場的各向異性

湍流風速標準差與湍流譜有如下關系

式中,下標i=u,v表示縱向、橫向.將式(14)及式(18)代入式(17),可得

對臺風風場的大量實測數據表明σv/σu在0.7至0.8 之間[30],本文中取 σv/σu=0.75,即km=0.28 .

湍流風模擬主要采用線性濾波法(AR 法)或諧波疊加法.線性濾波法[31]法將零均值隨機序列通過濾波器變換成具有指定頻譜特征的隨機過程.諧波疊加法[30,32-33]則將風速視為一系列正弦波的疊加,理論簡易且模擬結果精度較高,但較為耗時.由于線性濾波法算法繁瑣,因此,本文采用諧波疊加法在Matlab 軟件中模擬得到臺風各階段風場.模擬過程中采用Davenport 公式考慮各點風速之間的相關性,如式(20)所示,模擬過程中橫向衰減系數Cy為7,垂向衰減系數Cz為10,u為平均風速,i,j表示空間中任意兩點.圖1 所示為模擬得到的前眼壁強風(FEWS)階段輪轂處風速時程圖及由風速時程曲線得到的擬合譜與目標譜的對比,FEWS 階段輪轂處樣本平均風速49.05 m/s,輪轂處標準差5.985,標準差與目標值相差僅3%,可見模擬結果較為準確

圖1 FEWS 階段輪轂處風速曲線及譜對比Fig.1 Wind speed at hub height and the comparison between target spectrum and fitted spectrum

1.4 臺風浪模擬

海浪譜描述了海浪內部能量隨頻率的分布,目前海洋工程中常使用PM 譜與JONSWAP 譜模擬隨機海浪[34].PM 譜由無限風距測得,用于描述充分發展的海浪,JONSWAP 譜是在有限風距情況下測得的,適用于不同成長階段的海浪.現有研究表明,JONSWAP 譜能夠較好的描述臺風浪,本文使用Young[35]提出的經驗公式由十米高度處平均風速計算有義波高,結合JONSWAP 譜中譜峰頻率及有義波高的關系計算譜峰頻率

式中,Hs為有義波高,U10為十米高度處平均風速,g為重力加速度,E為波浪能量,可由有義波高得到,ε為無因次波浪能量,υ為無因次譜峰頻率,fp為譜峰頻率.JONSWAP 譜形式是在深水區域測得,在淺水區域,受海底影響,波譜產生變形.為考慮淺水效應,對得到的J 譜進行修正.將原有J 譜乘以一個修正因子Φ,該因子與水深h及波數k有關,如式(24)所示

波數k可借助色散關系由頻率ω計算,由于色散方程需通過迭代方法求解,為加快計算速度,本文采用下列6 位精度擬合多項式求解

式中,y為無因次波浪頻率,系數Cn取值見文獻[23].

2 數值計算模型

2.1 模型簡介

本文采用一體化仿真軟件SIMA (simulation workbench for marine application)開展耦合動力分析.在SIMA 中建立臺風經過全程,風浪聯合作用下數值模型,塔筒、樁基礎、葉片均采用歐拉梁單元建模.單樁基礎外徑9 m,厚度0.11 m,塔筒為多段等截面梁連接而成,每段梁沿長度方向外徑與厚度不變,塔筒底部外徑8.16 m,厚度0.07 m,塔筒頂部外徑5.65 m,厚度0.03 m.考慮到實際塔筒結構中螺栓、焊接、油漆、法蘭及電纜等各種設備的影響,將海床上方塔筒及樁基結構鋼的密度考慮為8500 kg/m3,略大于實際中鋼的密度7500 kg/m3,彈性模量取210 GPa,剪切模量取80.8 GPa.風力機使用由丹麥技術大學提出的DTU-10 MW 風機,該型風機具有較好的代表性,其主要參數見表1.

表1 DTU 10 MW 風力機參數Table 1 Main parameters of DTU 10 MW wind turbine

本文使用等效樁法考慮樁土作用,模型采用OC6 Phase2[36]中的單樁基礎.插入海床中的樁基礎總長為28.15 m,外徑9 m,內徑8.78 m.為模擬泥面線處轉角與位移,將插入海床下方的樁基分為兩部分,分別具有不同的長度與彈性模量,如圖2 所示,圖中Ll為23.15 m,L2為5 m,E1為821 GPa,E2為110.89 GPa.

圖2 樁土模型示意圖Fig.2 Illustration of the soil-structure interaction model

時域模擬采用Newmark-beta 方法求解式(27)所示動力方程,γ為1/2,β為1/4,式中Mt為質量矩陣,Kt為剛度矩陣,由結構構型在時刻t計算得出,分別為位移增量、速度增量、加速度增量.結構阻尼Ct采用了瑞利阻尼,對應一階模態1%的臨界阻尼

2.2 工況選取

本文重點探究臺風經過全程風浪聯合作用下大型單樁風機響應,基于文獻[8]中實測10 min 平均風速最大值及式(21)～式(23),經淺水修正后計算出臺風經過前后眼壁強風區及渦旋前后緣大風區的波浪條件.需要指出,以上公式是基于中國海域大量臺風期間海浪實測數據基礎上提出的半經驗方法,采用國際驗證的DTU 10 MW 風機模型,具有一定的局限性.工況如表2 所示.單次模擬時長4000 s,每個工況采用不同種子數后模擬十次,其中,去除400 s 的瞬態反應.坐標系原點位于靜水面上塔筒中心處,坐標軸如圖3 所示,模擬過程中風浪同向,均沿x 軸正向,實際海況下受海底地形及涌浪影響,風浪不一定同向.風機葉片在模擬開始后變槳至90°,之后一直處于順槳狀態.葉片位置1 豎直向上,位置2、位置3 依次順時針旋轉120°.

表2 工況定義Table 2 Definition of load cases

圖3 海上風機模型示意圖(單位:m)Fig.3 Illustration of the offshore wind turbine model (unit:m)

3 結果與討論

3.1 臺風經過全程運動響應特性分析

對于單樁式海上風力發電機,塔頂、塔基及泥面線運動較能體現出安全性能,因此,本文對以上結果進行對比分析,鑒于篇幅有限,僅給出塔頂運動全過程時程曲線.

圖4(a)～圖4(c)所示為臺風過境全過程塔頂在前后向(FA)運動的時程曲線,圖5(a)為塔頂運動的功率譜密度(PSD)圖.從時程圖可知五階段塔頂運動具有相似的波動特性,響應譜圖顯示五階段運動頻率均存在兩個峰值,控制頻率約為1.7 rad/s,接近塔筒FA 向固有頻率,另一峰值較小,接近各階段波浪頻率.同時在風速較高的FEWS 階段及BEWS 階段,塔頂運動在低頻區域分布有能量.表3 為塔頂運動統計值,表中可見塔頂FA 向運動均值隨平均風速增加而增大,且表現出非線性性質:BOVS 階段至FOVS 階段、BEWS 階段至FEWS 階段平均風速均增加5 m/s,但塔頂運動均值分別增加了0.028 m 與0.06 m,增加了約2 倍,PSD 圖也體現出這一點,FEWS 階段PSD 約為BEWS 階段兩倍.這對設計載荷的準確性提出了較高的要求,高風速下較小的風速差值可能導致較大的響應差別.

圖4 臺風經過全過程塔頂FA 向運動響應Fig.4 Tower top-aft (FA) motion during typhoon event

圖5 運動響應功率譜密度Fig.5 PSD of motion response

表3 塔頂運動時程統計Table 3 Statistics of tower-top motion

圖5(b)和圖5(c)分別為塔基、泥面線FA 向運動響應譜,圖5(b)中可見臺風經過各階段塔基處運動頻率表現出相似特性:各階段控制頻率均為塔筒一階頻率,同時受到波浪頻率影響.值得注意的是盡管BEWS 階段風速及波高均大于BOVS 階段,一階頻率處BOVS 階段塔基及泥面線PSD 均大于BEWS 階段,這可能是由于BEWS 階段結構的氣動阻尼較大.有義波高較低時,波浪激發了塔筒一階頻率,隨著有義波高增加,響應譜能量向波頻轉移.在臺風眼區域外,樁基礎泥面線處FA 向運動在一階頻率的功率譜密度隨風速增高略有增加,但相差不大.

3.2 臺風經過全程載荷特性分析

本節將重點分析塔頂、塔基、泥面線處FA 向彎矩與剪力響應,限于篇幅,僅給出FEWS 階段塔頂動力響應時域曲線.

圖6 為臺風全過程塔頂剪力時程曲線,圖7 為臺風經過全過程塔頂、塔基、泥面線處FA 向彎矩及剪力的時程統計圖,規定使塔筒向x正向彎曲的彎矩為正,圖8 為三個位置處的彎矩及剪力響應譜.

圖6 塔頂剪力Fig.6 Tower top shear force

圖7 載荷響應統計圖Fig.7 Statistics of load response

圖6 所示時程曲線在模擬時使用了相同的隨機波浪種子及不同的隨機風速種子,結果顯示在風速較低的TES,FOVS,BOVS 區域塔頂剪力的相位保持一致,而FEWS 階段及BEWS 階段則不一致.結合圖5(a)可知塔頂剪力波動主要受波浪激發的塔筒運動所引起的慣性力影響.風速較低時,塔筒頂部運動受一階頻率控制,且響應譜較窄,因此同一種子數下塔頂剪力相位一致.風速較高時,塔筒頂部運動PSD 在一階頻率處較寬且在低頻及波頻處分布有較多能量,因此相位出現離散.從圖7(a) 中可見TES 階段塔頂彎矩均值約-5600 kN·m,FEWS 階段塔頂彎矩均值約-2800 kN·m,可見風速的增加有助于抵消風力機重力造成的彎矩,但風速增大極大地增加了彎矩的標準差,FEWS 階段達到了3600 kN·m.塔頂剪力均值反應了上部風機葉片所受風力,圖7(b)為塔頂及塔基處剪力均值隨風速的變化,由于各風速下剪力均為負值,圖中取絕對值.圖中可見塔頂剪力均值的最大值出現在FEWS 階段,約為200 kN,遠小于10 MW 風機額定推力1500 kN.可見順槳有效降低了葉片所受風載荷.且塔頂剪力隨風速變化較小,在各階段變化率接近,而塔基剪力隨風速變化較大,高風速下剪力增長率更高.塔基剪力均值來自塔筒及風力機葉片所受風載荷,對比可知變槳停機狀態下風力機塔筒所受風載荷遠大于葉片.由圖7(c)可知泥面線處剪力均值與塔基處接近,但方差遠大于泥面線處.

圖8 為各截面處剪力和彎矩的PSD 圖.由圖可知,塔基剪力控制頻率為FA 向一階頻率,同時在低頻處有風頻反應.而泥面線處剪力控制頻率則為波浪頻率,波浪載荷極大地增加了泥面線處剪力的方差.塔基FA 向彎矩PSD 在各階段均只存在一個峰值,控制頻率為一階自振頻率.泥面線處彎矩PSD 存在兩個峰值,各階段控制頻率均為一階自振頻率.

圖8 載荷響應功率譜密度Fig.8 PSD of load response

圖9 為葉根FA 向彎矩及剪力時程統計圖,三個葉片的載荷響應差異明顯:圖9(a)為葉根彎矩響應統計圖,從圖中可見各階段葉片處于位置1 時葉根彎矩均值基本相同,而方差異明顯,FEWS 階段方差接近1500 kN·m.三個位置處葉根彎矩響應方差最大值均出現在FEWS 階段.圖9(b)為葉根剪力統計圖,圖中可見葉片處于位置2 時葉根剪力均值在三個葉片中最大,且其對風速變化不敏感.綜合可知,葉片位于位置2 時處于較危險位置.

圖9 葉片載荷響應統計圖Fig.9 Statistics of blade load

4 結論

本文基于DTU 10 MW 風機,應用諧波疊加法模擬臺風風場,開展了臺風經過全過程風浪聯合作用下大型單樁海上風機的耦合動力特性分析,可得出如下結論.

(1) 臺風經過全程中,FEWS 階段風機處于最危險位置,此時塔基及泥面線彎矩均值及標準差最大.

(2) 在風速45 m/s 下時,塔筒運動波動主要受波浪激發一階頻率控制,此時塔基上方各截面處結構動力載荷主要為慣性載荷,同樣也受一階頻率控制.而塔基下方泥面線處剪力及彎矩則表現出不同特性:剪力受波浪頻率控制,彎矩則受一階頻率控制.同時可看到一方面波高增加時引起的結構固有頻率處的響應能量增長有限,響應能量向波頻轉移.另一方面BEWS 階段及FEWS 階段風速相差不大,而塔基及泥面線處低頻響應差異較大.因此風速繼續增大,波高增加時響應特性可能出現較大變化.

(3)風力機葉片變槳有效降低了葉片所受風載荷,臺風經過各階段中,塔筒所受風力遠大于風力機葉片.三個葉片中,葉片處于位置2 時較危險,在臺風過后可對處于位置2 的葉片進行針對性檢修[37],減少維護成本,同時可開發針對波浪載荷的減振裝置,減小臺風過境下大型風力機的振動響應.

本文的研究工作僅考慮了單次臺風這一特例,由實測臺風數據推測海浪要素,重點關注臺風過境時風力機響應特性.下一步將基于實測的臺風和海浪數據,對風機的長期影響和安全評估進行深入研究.