風浪震及重力載荷聯(lián)合作用下大型風力機動力學響應及結構損傷分析

李志昊， 閆陽天， 李 春，2， 岳敏楠， 薛世成

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

為滿足日益增長的電力需求并減少化石能源引起的環(huán)境污染問題，可再生能源開發(fā)利用逐漸受到各國高度重視。海上風能作為可再生能源中最有潛力的能源之一，近年來發(fā)展迅速［1］。據(jù)Global Wind Energy Council（GWEC）預測，未來5年內全球海上風電新增裝機容量將超過469 GW［2］，未來10年超過40%的新增海上風力機將建于地震頻發(fā)區(qū)［3］，因此地震造成風力機結構破壞的風險不容忽視。單樁式風力機需置于海床內部以保證其穩(wěn)定性，但我國沿海地區(qū)地質條件較為復雜，土-結構耦合效應非線性明顯，因而由海底強烈地震引發(fā)的海床劇烈震蕩更易破壞原有基礎的穩(wěn)定性且易使具有靜定結構特性及高柔性的塔架振動加劇，甚至可能引發(fā)結構局部失效。因此，保證地震激勵下風力機結構安全極為重要。

關于地震載荷對風力機結構安全的影響，諸多學者展開了研究。Kim等［4］基于p-y曲線構建土-結構耦合效應模型并對地震作用下單樁式風力機進行了易損性分析，表明不同土層的地震相位差可能會導致更劇烈的地表運動但風力機動力學響應卻更小。Ma等［5］采用響應譜和時程分析法計算了地震作用下風力機動力學響應，表明土-結構耦合效應不可忽略。Yang等［6］以5 MW單樁式風力機模型為研究對象，對比分析了不同土-結構耦合效應模型對地震作用下風力機動力學響應的影響，表明將土-結構耦合效應簡化為剛性基礎可能會導致風力機動力學響應預估不準確。

此外，海上風力機安裝環(huán)境復雜，除地震載荷外，還受到湍流風載荷和波浪載荷沖擊。Zhao等［7］采用多體動力學方法計算風力機動力學響應，并用風力機輪轂處軸向推力代替氣動載荷，發(fā)現(xiàn)利用該簡化方法計算氣動載荷時可大幅提高計算效率，且該簡化方法下的動力學響應相差較小。Lesny等［8］基于p-y曲線構建土-結構耦合效應模型，研究了風浪載荷下風力機動力學響應，表明埋土樁基與土反力呈非線性關系。Cao等［9］基于非線性土-結構耦合效應模型對風浪載荷下風力機進行時頻分析，表明風浪載荷主要激發(fā)風力機一階模態(tài)。Banerjee等［10］采用有限元軟件計算了風浪聯(lián)合作用下5 MW單樁式風力機動力學響應，發(fā)現(xiàn)土-結構耦合效應將嚴重影響風浪載荷下的風力機動力學響應，此外，采用塔頂質量點既可保證一定的計算精度又可大幅節(jié)約計算資源。

上述研究僅考慮風浪或地震載荷，而關于風浪震載荷聯(lián)合作用下風力機動力學響應及穩(wěn)定性的研究較少。隨著海上風力機大型化發(fā)展，葉輪及支撐結構直徑大幅增加，湍流風和波浪載荷隨之增大，加之我國沿海地區(qū)地震多發(fā)，風浪震載荷聯(lián)合作用使得風力機塔架動力學響應更為復雜且失穩(wěn)風險進一步提高，這在評估風力機結構安全時不可忽略［11］。此外，以上研究均采用梁單元進行分析，難以獲取局部結構應力應變分布及結構失穩(wěn)現(xiàn)象，同時在研究結構響應時忽略了重力載荷。目前，為降低度電成本，風力機日趨大型化，重力載荷不斷增加，相較5 MW風力機，10 MW風力機葉片、機艙及輪轂等頂部結構和塔架的質量分別增加了326.7 t和852.5 t，塔架高度增加近30 m，柔性大幅增加，導致其更易受風浪震重力載荷而損壞［12］。因此，亟須開展風浪震重力載荷聯(lián)合作用下10 MW近海風力機動力學響應及穩(wěn)定性研究。因此，筆者選用丹麥技術大學（Danmarks Tekniske Universitet，DTU）10 MW單樁式風力機為研究對象，采用Kaimal風譜模型建立風場，通過P-M譜模擬波浪，選用實測地震位移數(shù)據(jù)為地震載荷，對不同環(huán)境載荷下風力機進行動力學響應及穩(wěn)定性研究，為風力機結構安全設計提供參考。

1 研究對象

以DTU 10 MW單樁式風力機為研究對象，模型圖及主要參數(shù)分別如圖1和表1［13］所示，其中X方向為前后向，Y方向為側向，Z方向為垂向。通過殼單元建立風力機有限元模型，離散模型如圖2所示，其中節(jié)點數(shù)量為1.927 8萬，網(wǎng)格數(shù)量為1.925萬。將葉片、輪轂及機艙等結構視為塔頂偏心質量點，并忽略焊縫、法蘭和油漆等因素，為避免風力機質量損失，將模型密度修正為8 500 kg／m3。

圖1 DTU 10 MW單樁式風力機示意圖Fig.1 Schematic diagram of DTU 10 MW monopile wind turbine

表1 風力機主要參數(shù)Tab.1 Main parameters of the wind turbine

圖2 風力機有限元離散模型網(wǎng)格Fig.2 Grid of finite element discrete model of the wind turbine

2 環(huán)境載荷

2.1 風載荷

通過TurbSim生成風力機輪轂處湍流風場，并假設其3個方向的速度分量相互獨立。采用國際電工委員會（IEC）Kaimal風譜模型，其功率譜密度Si（f）為：

式中：f為頻率；ˉuhub為輪轂處平均風速；M i和σi分別為積分尺度參數(shù)和標準差。

基于Kaimal風譜模型生成的湍流風場如圖3所示。

圖3 額定風速時輪轂處風場Fig.3 Wind field at hub height under rated wind speed

采用葉素動量理論求解風力機氣動載荷。將葉片沿展向分解為多段葉素，各葉素間相互獨立且互不影響，通過葉素速度三角形計算各葉素上的氣動力，并將其沿葉片展向積分，從而求解得到氣動載荷。圖4為葉素截面示意圖。其中，Ω為風輪轉速；R為風輪半徑；r為葉素展向半徑；U∞為來流風速；a和b分別為軸向和切向誘導因子；αi為攻角；βi為翼型扭角；φ為入流角。

圖4 第i段葉素截面及速度三角形示意圖Fig.4 Schematic diagram of blade element section i and velocity triangle

每個葉素的推力T和扭矩M［14］為：

式中：Cl為升力系數(shù)；Cd為阻力系數(shù)；ρa為空氣密度；c為翼型弦長；W為氣流相對速度。

選取風力機運行穩(wěn)定后100 s的湍流風載荷加載于風力機上，作用位置在風力機塔頂機艙處。

2.2 波浪載荷

波浪載荷可通過莫里森方程和繞射理論來求解。前者適用于計算小尺度結構的波浪載荷，并忽略了結構對入射波的影響，但隨著風力機大型化，風力機直徑不斷增大，其對入射波的影響不可忽略，因此選用繞射理論計算波浪載荷，作用位置在支撐結構處。選取可描述充分發(fā)展的波浪的P-M譜描述波浪分布，其功率譜密度［15］為：

式中：s（ω）為功率譜函數(shù)；ω為圓頻率；H s為有義波高，取3 m；T0為跨零周期，取6 s。

波浪表面由無限個余弦波疊加而成，其表面運動函數(shù)η（t）為：

式中：a k和ωk分別為組成波的振幅和圓頻率；εk為0～2π間均勻分布的初相位；k為波數(shù)。

基于線性波理論，結構波浪力F1（t）可表示為：

式中：ρw為海水密度；g為重力加速度；D為結構直徑；hw為水深；CM為水動力系數(shù)；z為笛卡爾坐標系空間坐標。

2.3 地震載荷

為準確描述地震發(fā)生時海床運動，從太平洋地震工程研究中心（PEERC）基于全球地震記錄建立的PEER NGA數(shù)據(jù)庫中選擇真實的地震運動數(shù)據(jù)。地震發(fā)生于Taiwan Chi-chi，震級為里氏7.62級，其地表位移時域曲線如圖5所示。

圖5 地震位移時域曲線Fig.5 Time-domain curve of seismic displacement

各有限元模型的環(huán)境載荷如表2所示。風浪震及重力載荷作用位置如圖6所示。

圖6 風浪震及重力載荷作用位置示意圖Fig.6 Location diagram of wind-wave-earthquake and gravity load

表2 不同環(huán)境載荷組合Tab.2 Different combinations of environmental loads

3 土-結構耦合效應

土體響應與結構運動間產生的相互作用稱為土-結構耦合效應。海床表面土質偏軟，土壤剛度對埋土樁基運動變化較為敏感，導致環(huán)境載荷作用下風力機土-結構耦合效應非線性十分明顯。為此，選取美國石油協(xié)會（API）規(guī)范推薦的p-y曲線及非線性彈簧單元建立土-結構耦合效應模型，該模型因應用方法簡便及可充分反映土層非線性而廣泛應用于近海工程［16］。埋土樁基總長度為30 m，間隔5 m設置垂直于樁基的彈簧以模擬樁周土水平抗力，其示意圖如圖7所示。

圖7 土-結構耦合效應模型Fig.7 Soil-structure coupling effect model

采用API推薦的雙曲正切模型計算p-y曲線，計算方法［17］如下：

式中：p為樁周土水平抗力；A為修正系數(shù)，取0.9；y為樁基水平方向位移（即樁基形變）；k1為地基反力系數(shù)；H為距泥面的距離（即土深）；pu為該處土體極限水平抗力；pus為淺層土壤極限土抗力；pud為深層土壤極限土抗力；D1為樁基直徑；γ為土體有效重度；C1、C2、C3為與砂礫內摩擦角的相關系數(shù)。

土壤參數(shù)如表3所示。不同土深處樁周土水平抗力計算結果見圖8。

圖8 不同土深處樁周土水平抗力Fig.8 Horizontal resistance of soil around piles in different depths

表3 土壤參數(shù)Tab.3 Soil parameters

4 有限元方法

4.1 模態(tài)分析

模態(tài)分析可確定結構的模態(tài)參數(shù)，即固有頻率及模態(tài)振型等，為結構振動特性分析提供論據(jù)且是瞬態(tài)動力學分析的基礎［18］。

無阻尼系統(tǒng)中，結構運動方程［19］為：

式中：M為質量矩陣；K為剛度矩陣；u··、u分別為節(jié)點加速度和位移向量。

4.2 瞬態(tài)動力學分析

瞬態(tài)動力學分析可通過求解結構運動方程來計算其在時變載荷下的形變和應力等動力學響應，運動方程［20］為：

式中：C為阻尼矩陣；u·為節(jié)點速度向量；F（t）為時變載荷。

4.3 屈曲分析

屈曲為工程設計中結構失效的一種形式，主要表現(xiàn)為當結構承受過大外載荷時，未完全發(fā)揮材料強度而失效。通過屈曲分析可獲得結構的屈曲模態(tài)、振型及臨界屈曲載荷。

結構平衡方程［21］為：

式中：Ke和Kg分別為彈性和初應力剛度矩陣；Δu為節(jié)點位移增量矢量；pj為結構承受的外載荷。

Kg可由載荷系數(shù)λ與單元幾何剛度矩陣ˉkg乘積求得：

將式（13）代入式（12）可得：

5 網(wǎng)格無關性驗證

有限元計算精度隨網(wǎng)格尺寸減小而提高，因此需要進行網(wǎng)格無關性驗證。根據(jù)結構內部最大應變能進行網(wǎng)格無關性驗證，不同網(wǎng)格數(shù)量對應的最大應變能如圖9所示。

圖9 不同網(wǎng)格數(shù)量對應的最大應變能Fig.9 The maximum strain energy under different grid numbers

由圖9可知，當網(wǎng)格數(shù)量從6 000增至19 250時，塔架最大應變能逐漸減小最后趨于不變，因此本文中19 250網(wǎng)格數(shù)量可保證較高計算精度。

6 結果及分析

6.1 模態(tài)

為避免塔架與葉片共振，可通過模態(tài)分析獲取塔架振動特性，如固有頻率及模態(tài)振型。通過計算精度高且收斂速度快的Block Lanczos方法提取塔架的前兩階固有頻率及模態(tài)振型。葉片及輪轂簡化前后的固有頻率對比如表4所示，模態(tài)振型如圖10所示。

由表4可知，將輪轂及葉片簡化成塔頂質量點對風力機模態(tài)特性的影響在允許范圍內。由圖10可知，塔架一、二階前后及側向模態(tài)振型分別為前后及側向彎曲振動，一階模態(tài)振型相對位移峰值位于塔頂，二階模態(tài)振型相對位移峰值位于塔架中部。

圖10 風力機前兩階模態(tài)振型Fig.10 The first two modes of the wind turbine

表4 風力機固有頻率Tab.4 Natural frequency of the wind turbine

6.2 時域

6.2.1 位移

塔頂位移為反映外載荷下風力機動力學響應的特征響應之一。為觀察不同環(huán)境載荷下風力機塔架動力學響應差異，輸出如圖11和圖12所示的塔頂前后向及側向位移時域響應曲線。

圖11 塔頂前后向位移時域曲線Fig.11 Time-domain curve of tower top front and rear displacement

圖12 塔頂側向位移時域曲線Fig.12 Time-domain curve of tower top lateral displacement

由圖11和圖12可知，不同環(huán)境載荷下風力機塔頂前后向及側向位移時域曲線波動形式有較大不同。圖11中，由模型1～模型3對比可知，地震未發(fā)生時風浪載荷為引起塔頂前后向位移的主要載荷。地震發(fā)生后，海床產生劇烈振動，導致原本僅受風浪載荷作用的風力機塔頂前后向位移急劇增大直至達到峰值且波動形式發(fā)生較大改變。模型1～模型3的塔頂前后向位移分別在25.75 s、7.86 s及38.28 s達到峰值，峰值分別為0.862 m、2.35 m及2.71 m，其中模型2的塔頂前后向位移在地震發(fā)生前達到峰值。相較于模型1和模型2，模型3塔頂前后向位移峰值增大22.3倍和18.3%。與模型3相比，模型4的塔頂前后向位移在地震發(fā)生前達到峰值，且地震發(fā)生后，塔頂前后向位移略微減小且其波動周期有明顯滯后性。此外，模型1～模型4的塔頂前后向位移標準差分別為0.401 m、0.526 m、0.664 m及0.849 m，相較于模型1和模型2，模型3的塔頂前后向位移標準差分別增大65.6%及26.2%，模型4的塔頂前后向位移標準差較模型3增大27.9%。因此，風浪震載荷共同作用時，地震發(fā)生后導致風力機塔頂前后向位移峰值及波動程度相較僅受風浪載荷大幅增大。此外，考慮重力載荷后，塔頂前后向位移波動周期明顯推遲且其波動程度進一步加劇。

由圖12可知，模型1～模型4的塔頂側向位移峰值分別為1.42 m、0.018 6 m、1.41 m及1.63 m。相較于模型1和模型2，模型3的塔頂側向位移峰值分別減小0.71%及增大74.8倍，而相較于模型3，模型4的塔頂側向位移峰值增大15.6%，模型1～模型4的塔頂側向位移標準差分別為0.605 m、0.005 99 m、0.603 m及0.704 m，相較于模型1和模型2，模型3的塔頂側向位移標準差分別減小0.331%及增大99.7倍，模型4的塔頂側向位移標準差較模型3增大16.7%。因此，塔頂側向位移主要由地震載荷引起，風浪載荷因其阻尼作用略微抑制塔頂側向位移峰值及波動程度，而重力載荷則略微增大了塔頂側向位移峰值及波動程度，且波動周期有明顯滯后性。

綜上所述，地震發(fā)生前，風浪載荷為導致風力機前后向位移主要載荷，地震發(fā)生后塔頂前后向和側向位移峰值及波動幅度均有不同程度增大，其中塔頂側向位移增幅更大，導致塔頂側向位移的載荷為地震載荷。此外，風浪載荷顯著增大塔頂前后向位移而對塔頂側向位移有略微抑制作用，重力載荷致使塔頂前后向及側向位移響應波動更加劇烈，其中塔頂前后向位移波動更劇烈，且使波動周期有一定滯后性。因此，日趨大型化且安裝于地震頻發(fā)區(qū)的風力機必須考慮風浪震重力載荷聯(lián)合作用，否則將導致對塔頂動力學響應預估不足且難以準確預估風力機塔頂位移達到峰值的時間。

6.2.2 等效應力

海上風力機支撐結構設計及建造成本花費巨大，約占總成本50%左右［22］，此外，復雜環(huán)境載荷也為支撐結構安全帶來巨大挑戰(zhàn)，易導致支撐結構等效應力集中，若等效應力過大可能致使支撐結構發(fā)生破壞。因此，對風浪震重力載荷聯(lián)合作用下支撐結構進行等效應力分析十分重要。

圖13給出了支撐結構等效應力時域響應曲線。由圖13可知，不同環(huán)境載荷導致支撐結構等效應力響應峰值及波動產生巨大差異。模型1～模型4支撐結構等效應力分別在26.76 s、7.98 s、38.3 s及23.75 s達到峰值，峰值分別為155.3 MPa、186.5 MPa、224.8 MPa及244.4 MPa，標準差分別為41.1 MPa、42.0 MPa、55.7 MPa及61.8 MPa，相較于模型1和模型2，模型3等效應力峰值分別增大44.8%和20.5%，其標準差增大35.5%和32.6%。相較于模型3，模型4等效應力峰值增大8.72%，其標準差增大11.0%。因此，相較于地震及風浪載荷，風浪震載荷聯(lián)合作用下風力機塔架等效應力峰值及波動程度都大幅增加，且達到峰值時間也略有延后，此外，考慮重力載荷后，塔架等效應力峰值及波動程度進一步增加。

圖13 支撐結構等效應力時域響應曲線Fig.13 Time-domain curve of equivalent stress of the support structure

為進一步觀察支撐結構等效應力分布，輸出如圖14所示的不同環(huán)境載荷下支撐結構等效應力達到峰值時的等效應力響應云圖。由圖14可知，模型1～模型4支撐結構迎風面等效應力峰值均位于支撐結構頂部，而云圖中支撐結構中部等效應力灰度值較低，等效應力集聚區(qū)域面積相較支撐結構兩端明顯更小，且云圖中相較于背風面，迎風面等效應力灰度值更高，高應力區(qū)域面積更大。此外，相較于模型1和模型2，模型3支撐結構處等效應力集聚嚴重，模型4等效應力集聚更甚。因此，地震載荷將會導致支撐結構等效應力響應更為劇烈，考慮重力載荷后，支撐結構等效應力響應進一步加劇。此外，必須考慮風浪震重力載荷聯(lián)合作用，否則將對風力機支撐結構等效應力響應預估不準確。

圖14 支撐結構等效應力響應云圖Fig.14 Cloud image of equivalent stress response of the support structure

6.3 頻域

由于塔頂位移時域響應為風載荷、波浪載荷、地震載荷及重力載荷等隨機載荷聯(lián)合作用結果，具有非線性及非平穩(wěn)特征，需通過頻域分析更深入地分析塔架振動特性。通過快速傅里葉變換（FFT）將塔頂位移時域響應轉換為頻域數(shù)據(jù)，如圖15所示。

由圖15可知，模型1～模型4前后及側向模態(tài)的一階固有頻率處均存在明顯位移峰值，且相較于模型1和模型2，模型3的位移峰值較大，而模型4的位移峰值最大。但只有模型3一階前后模態(tài)固有頻率的二倍頻處存在明顯位移峰值。因此，地震及風浪載荷僅激發(fā)塔架一階模態(tài)，風浪震載荷聯(lián)合作用下塔架一階前后模態(tài)二倍頻處將出現(xiàn)明顯位移峰值，而考慮重力載荷后，一階前后模態(tài)二倍頻處被抑制。此外，相較地震及風浪載荷，風浪震載荷聯(lián)合作用下塔架一階前后及側向模態(tài)處塔頂位移峰值更大，而風浪震重力載荷聯(lián)合作用下塔頂位移峰值最大。

圖15 塔頂位移頻域曲線Fig.15 Frequency-domain curve of tower top displacement

6.4 屈曲

隨著風力機日趨大型化，塔架柔性逐漸變大，且其結構為偏心受壓薄壁結構，這決定了局部屈曲為塔架主要破壞形式之一［23］。因此，對不同環(huán)境載荷下塔架進行屈曲分析，表5和圖16給出了不同環(huán)境載荷下塔架一階屈曲因子及模態(tài)。

表5 不同環(huán)境載荷下塔架一階屈曲因子Tab.5 First-order buckling factor of tower under different environmental loads

圖16 不同環(huán)境載荷下塔架一階屈曲模態(tài)Fig.16 First-order buckling mode of tower under different environmental loads

由表5可知，相較于模型1和模型2，模型3一階屈曲因子大幅降低，模型4一階屈曲因子進一步降低。因此，相較地震及風浪載荷，風浪震載荷聯(lián)合作用下塔架屈曲風險更大，風浪震重力載荷聯(lián)合作用下塔架屈曲風險進一步增大，必須考慮風浪震重力載荷聯(lián)合作用，否則將低估屈曲風險。

由圖16可知，僅模型2塔架屈曲區(qū)域為塔架底端，而其余模型屈曲區(qū)域均位于塔頂，這是由于塔頂壁厚較薄，相較其余地方剛度較小，更易發(fā)生屈曲。故工程中需要重點關注塔頂并附加防屈曲裝置，以降低風力機在風浪震重力載荷聯(lián)合作用下產生過大振動而發(fā)生屈曲的風險。

7 結 論

（1）風浪載荷為導致風力機發(fā)生塔頂前后向位移響應的主要載荷，地震載荷為引起塔頂側向位移的主要載荷。地震發(fā)生后塔頂前后向和側向位移響應均有所增大，相較于塔頂前后向位移，塔頂側向位移增幅更大。風浪載荷因其阻尼作用將略微抑制塔頂側向位移。重力載荷致使塔頂前后向及側向位移波動更為劇烈，且波動周期有一定滯后性，也改變了塔頂位移達到峰值的時間。

（2）相較于地震及風浪載荷，風浪震載荷聯(lián)合作用下支撐結構等效應力響應更大且波動更為劇烈，考慮重力載荷后等效應力響應進一步增大。不同環(huán)境載荷下塔架等效應力均集中于支撐結構上部。相較背風面，支撐結構迎風面等效應力集聚更嚴重。

（3）風力機結構設計時必須考慮風浪震重力載荷聯(lián)合作用，否則將導致風力機動力學響應預估不準確。

（4）不同環(huán)境載荷下塔架前后及側向一階固有頻率處均存在明顯位移峰值，即一階模態(tài)被環(huán)境載荷誘發(fā)。此外，風浪震載荷聯(lián)合作用下塔架一階前后模態(tài)頻率在二倍頻處出現(xiàn)明顯峰值。考慮重力載荷后，一階前后模態(tài)二倍頻處被抑制。

（5）不同環(huán)境載荷下，塔架屈曲因子略微不同。相較于地震及風浪載荷，風浪震載荷聯(lián)合作用下一階屈曲因子大幅降低，更易發(fā)生局部屈曲，考慮重力載荷后，一階屈曲因子進一步降低。此外，塔頂為易發(fā)生屈曲區(qū)域，在工程設計時應重點關注此處，可通過附加防屈曲結構等方式增大塔頂剛度。