基于三維地震作用下樁-土耦合效應近海超大型風力機動力學響應研究

李志昊， 岳敏楠， 李 春,2， 閆陽天， 楊 陽

(1.上海理工大學 能源與動力工程學院，上海 200093； 2.上海市動力工程多相流動與傳熱重點實驗室，上海 200093；3.Department of Maritime and Mechanical Engineering, Liverpool John Moores University, Byrom Street, Liverpool L3 3AF, UK)

海上風能因風速穩定、可開發區域廣及對人類影響小等優勢逐漸成為最具增長潛力的可再生能源之一[1]。至2019年，全球海上風電新增裝機容量27.2 GW，另據國際能源署(International Energy Agency)預測，自2018—2040年，全球海上風電裝機容量將增加15倍[2-3]。

為降低度電成本，海上風力機呈大型化發展趨勢，其中因技術成熟、安裝簡單及成本低等優點，多采用固定式支撐結構風力機。其支撐結構可分為單樁式、重力式、三樁式及導管架式。因單樁式支撐結構具有安裝簡單和成本低等優勢，成為近海風力機中應用最廣的支撐結構[4]。然而，隨單樁式風力機逐漸大型化，其各組成結構尺寸顯著增加，多采用降低單個部件重量以減小安裝難度，此雖增加了結構柔性但也更易受地震激勵影響[5]。由于單樁式風力機多安裝于地震頻發且海床土質松軟的東南沿海地區風場，其地質特征導致樁-土耦合效應明顯，樁基極有可能因地震時的土體劇烈震動而遭破壞，導致支撐結構發生局部失穩及屬高聳結構的風力機塔架振動加劇。因此，研究準確的樁-土耦合效應模型及地震作用下超大型單樁式風力機動力學響應尤為重要。

單樁式海上風力機樁基直徑達9 m，屬于超大直徑受荷樁，其樁-土耦合效應廣受重視。徐海濱等[6]基于莫爾-庫倫模型建立土體模型，并模擬了樁徑為6 m的大直徑樁的樁-土耦合效應。Yang等[7]基于有限元法及現場水平載荷試驗進行對比分析。Min等[8]等采用有限元軟件建立砂土土體有限元模型，在此基礎上研究了樁的受力特點。Byrne等[9]以大長徑比的樁為研究對象，提出一種新的樁-土耦合效應模型。Abdullahi等[10]以REL 5 MW單樁式風力機為研究對象，基于分布式彈簧和p-y曲線建立樁-土耦合效應模型，并將其與四彈簧樁-土耦合效應模型對比，發現二者雖有一定誤差，但誤差較小。

以上針對樁-土耦合效應模型的研究均未考慮地震載荷，由上文可知，地震對近海風力機造成較大威脅。為此，國內外諸多學者開展了地震作用下風力機動力學響應研究。陳法波[11]通過時域法和頻域法研究地震作用下3.6 MW風力機動力學響應，發現頻域法因無法考慮結構非線性，致其計算結果誤差較大。季亮等[12]基于時域法和頻域法分析了NREL 5 MW風力機在地震作用下的塔底應力、彎矩及剪力，結果表明：頻域法計算結果較時域法偏低。因此，采用時域法分析地震作用下風力機動力學響應十分重要。

地震激勵作用下風力機結構響應是準確計算結構承載力和施工成本的重要依據，而樁-土耦合效應是評估地震作用下結構響應的重要因素[13]。劉中勝等[14]基于Wolf方法建立樁-土耦合效應模型，并通過風力機多體動力學仿真軟件FAST模擬考慮樁-土耦合效應時風力機在不同地震強度下結構響應，時程結果表明：樁-土耦合效應對風力機在地震作用下的結構響應影響顯著，不可忽略。Ma等[15]對比分析地震下是否考慮樁-土耦合效應的3 MW海上單樁式風力機瞬態動力學響應研究發現：樁-土耦合效應明顯改變塔架一階固有頻率及塔架振動響應。因此，對于地震作用下風力機動力學響應，必須考慮樁-土耦合效應。

基于p-y曲線法建立的樁-土耦合效應模型因簡便、具有一定精度且可較好地描述水平地震作用下樁-土耦合效應，在學術界得到廣泛應用[16]。龐文彥[17]以1.5 MW風力機為研究對象，對比分析現有計算樁基水平承載力方法，結果表明：p-y曲線可較好地描述樁基與土體間相互作用。王明超等[18]基于線性p-y曲線確定三種描述樁-土耦合效應的彈簧模型剛度，并研究其對單樁式風力機前兩階固有頻率影響，但線性p-y曲線尚無法體現土反力與樁基撓度間非線性關系。Lesny等[19]基于p-y曲線研究5 MW風力機在風、浪載荷作用下的結構動力學響應表明：當樁基受載荷作用產生撓曲時，樁周土反力與樁基撓度呈非線性關系。Reese[20]通過隨土深變化的非線性p-y曲線模擬風力機樁基埋于砂土中的樁-土耦合效應，然未考慮地震作用影響。Ju等[21]以非線性p-y曲線建立樁-土耦合效應模型，并以此為基礎研究NREL 5 MW 桁架式風力機在風、浪和地震作用下的動力學響應，結果表明：地震作用于風力機時，必須考慮樁-土耦合效應；此外，地震加速度峰值大于0.52g時，較風浪載荷，地震為風力機控制載荷。

顯然，地震對結構的空間影響是三維的，地震不僅引起水平方向土體運動，也將引發垂向土體運動。但，上述研究未考慮地震作用影響或僅分析水平而忽略垂向地震作用下樁-土耦合效應，無法完整反映三維地震作用下樁-土耦合效應。Kj?rlaug等[22]對比分析是否考慮垂向地震作用下樁-土耦合效應的陸上風力機動力學響應，發現考慮垂向地震激勵下樁-土耦合效應使塔架加速度響應放大2倍～5倍，但未對海上風力機垂向地震作用下樁-土耦合效應展開研究。此外，隨風力機日趨大型化，而目前針對超大型10 MW單樁式風力機研究較少，且上述文獻均僅以地震加速度數據為地震作用于風力機，其未能真實反映地表運動，而能準確反映地表運動的地震位移數據作為地震作用于風力機卻鮮有研究。

為此，選取實測地震位移數據為地震作用于風力機，基于p-y、t-z及q-z曲線建立樁-土耦合效應模型，并以超大型DTU 10 MW單樁式近海風力機為研究對象，對比分析不同樁-土耦合效應下動力學響應，為風力機結構抗震設計提供參考。

1 地震作用

為準確描述地震發生時地表運動，選取Pacific Earthquake Engineering Research(PEER)[23]提供的1999年發生于臺灣里氏7.62級地震位移數據，地震位移峰值為0.146 m，地震位移時域曲線如圖1所示，圖中X為流向，Y為橫向，Z為垂向。

圖1 地震位移時域曲線Fig.1 Seismic displacement time domain curve

2 模型建立

2.1 超大型單樁式風力機

選取近海超大型DTU 10 MW單樁式風力機為研究對象，如圖2，主要參數見表1[24]。分析模型中，采用殼單元建模并將葉片、輪轂及機艙等按實際質量等效為塔頂偏心質量點，該模型采用結構化網格，網格尺寸為0.5 m，網格數為25 074，節點數為25 423，邊界條件為流向、橫向及垂向地震發生時地表位移數據，此外，鑒于油漆、焊接及法蘭等因素，計算模型有效密度取8 500 kg/m3。以瑞利阻尼反映結構阻尼，瑞利阻尼系數按下式確定

(1)

(2)

式中：ζ為阻尼比，1%；ωi、ωj為自振頻率，分別為0.2及0.21 Hz；α為0.002，β為0.05[25]。

圖2 風力機及樁-土耦合效應示意圖Fig.2 Schematic diagram of wind turbine and pile-soilcoupling effect

表1 DTU 10 MW風力機主要參數Tab.1 Main parameters of DTU 10 MW wind turbine

2.2 樁-土耦合效應模型

地震發生時，水平和垂直方向地震波通過土層傳遞并作用于結構，結構由于振動所產生水平及垂直方向慣性力反作用于土層，稱為三維地震作用下樁-土耦合效應[26]。本文采用美國石油協會(American Petroleum Institute，API)推薦的p-y、t-z及q-z曲線彈簧模型建立三維地震作用下樁-土耦合效應模型，分別模擬水平及垂直方向地震作用下的樁-土耦合效應，即樁側土水平抗力與樁撓度間關系、樁周土摩擦力及樁周土摩阻力與樁沉降距離間關系，該模型因綜合考慮土體非線性、不同土體種類及各土層特性等因素廣泛應用于近海工程[27]。樁基長30 m，以1 m及5 m為間隔布置沿樁身方向和垂直于樁的彈簧模擬樁側土水平抗力和樁周土摩擦力，樁端則通過垂向彈簧模擬樁周土摩阻力，選取非線性彈簧單元(Combin 39)模擬彈簧約束，見圖2。該彈簧單元在一、二及三維應用中具有拉伸和扭轉能力且可模擬力與變形間非線性關系軸向單元，其需定義力-變形(F-D)曲線以確定彈簧受力，如圖3所示。

圖3 非線性彈簧單元力-變形曲線Fig.3 Nonlinear spring element force-deformation curve

如圖3所示，非線性彈簧單元將根據輸入的離散點以擬合成力-變形曲線，從而實現彈簧非線性變形。

為觀察樁周土摩阻力及三維地震作用下樁-土耦合效應對風力機的影響，以單樁式風力機為基礎建立三種樁-土耦合效應模型，詳情如表2所示。

表2 不同樁-土耦合效應模型

2.2.1p-y曲線

按API推薦，砂性土p-y曲線表達式為

p=Aputanh[kHy/(Apu)]

(3)

pus=(C1H+C2D)γH

(4)

pud=C3DγH

(5)

pu=min{pus,pud}

(6)

式中:pus為淺層土壤極限土抗力；pud為深層土壤極限土抗力；C1、C2、C3為砂礫摩擦角函數；D為樁徑；γ為土體有效容重；H為土體深度；A為經驗參數，循環載荷時取0.9，定載荷時取(3～0.8H/D)≥0.9；pu為深度為H時土體極限抗力；y為樁基撓度；k為初始土體切變模量。土壤參數如表3所示。

表3 土壤參數

p-y曲線計算結果如圖4所示。

2.2.2t-z曲線

砂性土t-z曲線可由式(7)、(8)和表4得到，其表達式為

fs=KP0tanδ

(7)

tmax=fs·AS

(8)

式中:K為無因次土層側壓力系數；P0為計算點處有效復土壓力，P0=γH，γ為土體有效容重；δ為樁土間摩擦角；AS為樁身面積；fs為單位面積極限樁周土摩擦力；tmax為極限樁周土摩擦力。

圖4 不同土深處p-y曲線Fig.4 p-y curves at different depths

表4 砂性土t-z曲線參數Tab.4 Parameter of t-z curve of sandy soil

表4中:z為樁沉降距離；D為樁直徑；t為樁周土摩擦力。t-z曲線計算結果如圖5所示。

圖5 不同土深處t-z曲線Fig.5 t-z curves at different depths

2.2.3q-z曲線

將式(9)、(10)及表4聯立可得q-z曲線，計算式為：

qp=NqP0

(9)

qmax=qp·Ap

(10)

式中：Nq為無量綱承載力系數；qp為樁端單位面積極限阻力；Ap為樁端橫截面積；qmax為極限樁周土摩阻力。

表5 砂性土q-z曲線參數

計算q-z曲線結果由圖6所示。

圖6 q-z曲線Fig.6 q-z curves

2.2.4 三維實體土體模型

為衡量三種樁-土耦合效應模型的準確性，采用基于粘彈性理論的等效線性模型作為土體的本構模型并建立三維實體土體模型以模擬樁-土耦合效應，土體模型如圖7所示。

3 驗 證

原有整機在進行葉片及輪轂的簡化后，其將造成一定誤差，為此對不考慮樁-土耦合效應的整機模型與簡化模型分別進行模態計算，以此判斷簡化模型與整機結構動力學特性差異，前四階固有頻率如表6所示。

圖7 三維實體土體模型Fig.7 3-D solid soil model

表6 簡化模型與整機結構前四階固有頻率Tab.6 Simplified model and the first four natural frequencies of the machine structure

由表6可知，簡化模型與整機結構固有頻率誤差范圍較小，均在5%以內。因此，可判斷簡化模型與整機結構動力學特性在誤差允許范圍內及二者一致性較好。

4 結果及分析

4.1 模 態

為研究垂向樁-土耦合效應對塔架動力學特性的影響，采用Block Lanczos方法分別對模型1、2和3模態分析，將得到的模態振型和固有頻率對比分析。模型1、2和3前四階模態振型如圖8所示，前六階模態固有頻率如表7所示。

(a) 模型1

表7 模型1、2和3前六階模態固有頻率

由圖8可知，模型1和2前四階模態振型相似，模型3前兩階模態振型與模型1及2類似，三階及四階模態振型發生較大改變。模型1和2塔架前兩階模態振型發生前后彎曲振動，三階及四階模態振型為側向彎曲振動，而模型3塔架前四階模態振型均為前后彎曲振動，且模型1和2三階及四階模態振型相對位移峰值位于塔架中部，模型3三階及四階模態振型相對位移峰值位于塔頂。此外，愈加復雜的樁-土耦合效應下，塔架模態振型振動相對位移峰值愈大。因此，樁周土摩阻力未對塔架前四階模態振型產生較大影響，僅增加塔架前四階模態振型相對位移，而垂向樁-土耦合效應改變塔架三階及四階模態振型及相對位移峰值位置，同時亦加劇前四階模態相對位移，且其增幅更大。

由表7可知，較之模型1，模型2和3塔架前四階固有頻率均有不同程度變化，其中模型3變化幅度更大，而三種模型的五階及六階固有頻率差異較小。故前四階固有頻率受樁周土摩阻力和垂向樁-土耦合效應影響顯著，受后者影響更甚，而高階模態固有頻率受二者影響較小。由于未考慮門洞和爬梯等設備，塔架為幾乎完全對稱結構，因此相鄰兩階固有頻率較為接近，進一步說明計算結果的準確性。此外，工程實際中，常以塔架一階固有頻率是否避開風輪旋轉頻率及葉片通過頻率(1P：0.230～0.288 Hz及3P：0.779～0.865 Hz)判斷塔架是否將與風輪及葉產生共振，通常認為1P～3P為安全區，為更直觀地觀察塔架一階固有頻率是否在安全區內，將其與安全區范圍繪于圖9中比較[28]。

由圖可知，模型1、2和3塔架均位于安全區內，但較之模型1和2，模型3塔架固有頻率安全裕度更大。因此，考慮垂向樁-土耦合效應后，塔架固有頻率遠離共振頻率，更不易與上部結構發生共振。

4.2 動力學響應

4.2.1 塔 架

風力機大型化導致塔架柔性大幅增加，使其在極具破壞力的地震作用下易發生過大振動造成風力機損壞，為監測地震過程中不同樁-土耦合效應下的塔架振動響應，采用求解精度高的New-mark積分法對其進行動力學響應分析[29]，其中塔頂加速度為反映塔架振動的重要特征參數。為判斷及分析地震作用下三種樁-土耦合效應模型與三維實體土體模型塔架的動力學響應情況，輸出如圖10所示的塔架加速度時域曲線。

圖9 模型與塔架固有頻率安全區對比示意圖Fig.9 Comparison diagram of the model and the naturalfrequency safe zone of the tower

圖10 塔頂加速度時域曲線Fig.10 Time domain curve of tower top acceleration

由圖10可知，模型3與三維實體土體模型的塔頂加速度響應最為接近，因此模型3與實際情況最為接近。此外，不同樁-土耦合效應模型塔頂加速度時域曲線差別較大。模型1、2和3塔頂加速度分別于29.7 s、16.175 s及15.5 s達到峰值，后因結構阻尼和地震能量衰減，塔頂加速度響應逐漸減小，峰值分別為2.82 m/s2、4.21 m/s2及6.31 m/s2，標準差分別為0.54 m/s2、0.55 m/s2及0.82 m/s2。較模型1，模型2及3位移響應峰值分別增大49%和124%，標準差分別增加2%和52%。因此，較之僅考慮水平樁-土耦合效應，樁周土摩阻力和三維地震作用下樁-土耦合效應將導致塔頂加速度響應更早達到峰值，且峰值和波動均增加，其中三維地震作用下樁-土耦合效應增幅更甚。綜上所述，若僅采用p-y曲線及p-y、q-z曲線模擬樁-土耦合效應，將導致對地震作用下塔架動力學響應評估不足。

圖11、12及13為三種樁-土耦合效應塔頂流、橫和垂向加速度響應，由圖可知樁周土摩阻力及三維樁-土耦合效應不可忽略。較之橫和垂向，流向塔頂加速度響應峰值更大，并流和橫向塔頂加速度響應均以0為中線上下振蕩。為詳細分析流、橫及垂向加速度波動情況，輸出模型1、2及3的流、橫和垂向加速度變化范圍如表8所示。

圖11 流向塔頂加速度時域曲線Fig.11 Time domain curve of acceleration to the top of flow tower

圖12 橫向塔頂加速度時域曲線Fig.12 Time domain curve of horizontal top acceleration

圖13 垂向塔頂加速度時域曲線Fig.13 Vertical tower top acceleration time domain curve

較之模型1，模型2和3三個方向塔頂加速度響應波動均有所增大，其中橫向塔頂加速響應變化范圍增幅最大，分別增大209%和317%，且三維地震作用下樁-土耦合效應相較于垂向地震作用中僅考慮樁周土摩阻力，流、橫和垂向塔頂加速度響應變化范圍均有增大，其中流向塔頂加速度響應變化范圍增幅最大，增加169%。

4.2.2 支撐結構

單樁式風力機塔架、葉片、輪轂和機艙等設備完全由支撐結構支撐，而支撐結構在地震作用下易發生應力集中現象，若應力過大將導致其局部屈服甚至結構破壞，使風力機存在倒塌風險，因此支撐結構處的等效應力分析十分重要，輸出如圖14所示的不同樁-土耦合效應時最大等效應力響應云圖。

表8 模型1、2和3流、橫和垂向加速度變化范圍

圖14 模型1、2和3支撐結構最大等效應力響應云圖Fig.14 Models 1, 2 and 3 support structures respond tomaximum equivalent stress clouds

由圖14可知，樁周土摩阻力和三維地震作用下樁-土耦合效應未改變最大等效應力位置，其為過渡段和入水段連接處。對比不同模型最大等效應力可得，較模型1，模型2和3最大等效應力分別減少283%和566%，且模型2及3應力集中區出現明顯縮減，模型3應力集中區最小。

此外，由結構鋼的屈服極限為250 MPa，強度極限為460 MPa可知，模型1超過材料強度極限460 MPa，結構發生破壞，模型2和3均未超過材料強度極限，其中模型2等效應力大于屈服極限250 MPa，發生局部屈服。因此，樁周土摩阻力降低了支撐結構最大等效應力，三維地震作用下樁-土耦合效應降幅更大。

4.2.3 樁 基

樁基為樁-土耦合效應直接作用對象，因其深埋土體中對其修復異常困難且花費巨大。但對于樁基分析，尤其是考慮垂向地震作用樁-土耦合效應下的樁基分析常被忽略，為分析垂向地震作用樁-土耦合效應對樁基的影響，對模型1、2及3的樁基進行動力學分析。

由圖15可知，模型2和3較模型1樁基位移響應明顯增大，模型3位移響應最劇烈。此外，模型1和2地震開始作用后，樁基位移響應迅速增大，達到峰值后，樁基位移響應迅速減小，而模型3樁基位移響應激增至峰值后，樁基位移開始衰減但衰減速度遠小于模型1和模型2。由此可見，樁周土摩阻力和三維地震作用下樁-土耦合效應均使樁基位移響應顯著增加，三維地震作用下樁-土耦合效應下樁基位移響應最劇烈，同時，考慮三維地震作用下樁-土耦合效應地震后期樁基位移響應衰減速度大幅下降。

圖15 樁基位移時域曲線Fig.15 Time domain curve of pile foundation displacement

由圖16可知，較之模型1，模型2和3樁基各高度處位移響應均有不同程度增大，模型3增幅更大。同時，模型1、2及3樁基位移響應峰值均位于樁基頂部，且模型1和3樁基位移響應隨樁基高度增大而減小。綜上所述，樁周土摩阻力及三維地震作用下樁-土耦合效應未改變樁基位移峰值位置，均位于樁頂，但其增大了樁基各高度處位移響應，三維地震作用下樁-土耦合效應增幅最大。同時，樁基位移響應隨樁基高度增加而減少。

圖16 不同高度處樁基位移響應Fig.16 Pile foundation displacement response at different heights

綜上所述，三維地震作用下樁-土耦合效應與僅考慮水平樁-土耦合效應風力機動力學響應有很大不同，若僅以p-y曲線及p-y、q-z曲線模擬樁-土耦合效應，將導致對地震作用下風力機動力學響應的評估嚴重不足。

5 結 論

基于p-y、t-z及q-z曲線描述樁-土耦合效應，以超大型DTU 10 MW單樁式風力機為研究對象，分析不同樁-土耦合效應下風力機動力學響應，結論如下：

(1) 采用集中質量點模擬葉片等上部結構的簡化模型與風力機整機結構動力學特性差異較小；樁周土摩阻力和垂向樁-土耦合效應使前四模態相對位移增加，垂向樁-土耦合效應增幅更大，且樁周土摩阻力未明顯改變塔架前四階模態振型及位移峰值位置，考慮垂向樁-土耦合效應后，三階及四階模態振型由側向彎曲振動變為前后彎曲振動且位移峰值位置由塔架中部移至頂部，此外，塔架固有頻率安全裕度更大。

(2) 樁周土摩阻力和三維地震作用下樁-土耦合效應增大了塔頂加速度響應波動及峰值，峰值分別增加49%和124%，標準差分別增大2%和52%，并使得流、橫及垂向加速度響應變化范圍激增，其中橫向增幅最大，分別增加209%和318%，且較樁周土摩阻力，三維地震作用下樁-土耦合效應使流向塔頂加速度響應增幅最大，增大169%；若僅以p-y曲線及p-y、q-z曲線模擬樁-土耦合效應，將導致對地震作用下風力機動力學響應的評估嚴重不足。

(3) 樁周土摩阻力和三維地震作用下樁-土耦合效應未改變支撐結構最大等效應力位置，位于過渡段和入水段連接處，且其降低了最大等效應力，峰值減小了283%及566%，并減小了應力集中區。

(4) 樁周土摩阻力和三維地震作用下樁-土耦合效應未改變樁基位移峰值位置，位于樁頂；此外，較之僅考慮水平載荷樁-土耦合效應，樁周土摩阻力和三維地震作用下樁-土耦合效應下樁基及其各高度處位移響應均有不同程度增大，其中三維地震作用下樁-土耦合效應下樁基位移增幅更大，且使地震后期樁基位移響應衰減速度大幅下降。