非均質土中海上風電單樁基礎動力響應特性

孔德森，劉 一，鄧美旭，李亞洲

1) 山東科技大學山東省土木工程防災減災重點實驗室，青島 266590 2) 山東科技大學大學土木工程與建筑學院，青島 266590 3) 中鐵建工集團山東有限公司，青島 266100

風能源已經逐漸成為取代傳統能源發電的一種主要形式，具有無污染、可再生等特點，是一種取之不盡、用之不竭的新能源[1]. 我國具有漫長的海岸線，近海資源豐富，空間區域廣闊，海上風電場逐漸成為風力發電的開發重點[2?3]. 大直徑單樁基礎由于其制作與安裝方便快捷、相對于群樁而言受力條件明確、造價經濟等優勢，在海上發電領域得到了廣泛應用[4]. 大直徑單樁基礎在服役期間內受到來自風、洋流、波浪等水平荷載的作用，其水平循環受荷特性顯著[5]. 研究大直徑單樁基礎在水平循環荷載下的動力響應特性，對于消除安裝工程中的隱患，保證服役期間正常使用均具有重要意義.

關于水平循環荷載下單樁基礎動力響應特性，國內外學者進行了相關研究，并且取得了一定的成果. Leblanc等[6]和Peng等[7]利用離心試驗研究循環次數、樁徑、荷載大小以及砂土密實度對模型樁側向位移的影響. 郭鵬飛等[8]基于非飽和土的動力控制方程，考慮橫向慣性效應，建立了三相非飽和介質中樁的豎向動力響應連續介質模型，得到了樁側土體剪應力及豎向振動位移的表達式. 張光建[9]對ABAQUS進行二次開發建立土體剛度衰減模型，重點研究了不同影響因素下大直徑單樁基礎水平位移的發展趨勢. Basack與Banerjee[10]基于數值模型研究了層狀土中樁基礎在橫向荷載作用下的響應，同時利用現場試驗數據對模型進行了驗證，并進行了參數化研究.Bhattacharya與Adhikari[11]對樁?土相互作用進行了實驗研究，結果表明，考慮樁?土相互作用時設計出的基礎結構使用壽命更長. 朱斌等[12]研究了砂土中大直徑單樁離心模型試驗結果，定義了臨界循環應力比，給出了其與循環折減系數的關系.Kuo等[13]和Achmus等[14]建立了砂土的循環剛度衰減模型，得到了樁?土體系在N次循環后的承載變形特性. 但是對于水平循環荷載作用下非均質土中海上風電單樁基礎動力響應特性的研究比較少.

本文以海上風電單樁基礎為研究對象，采用有限元分析軟件ABAQUS建立了非均質土中海上風電單樁基礎數值計算模型進行研究，在模型中將樁基礎受到的波浪、洋流、風荷載等效成雙向對稱循環荷載，對水平循環荷載作用下樁身水平位移、樁身剪力、樁身彎矩、樁內外壁土抗力沿埋深變化特性進行了研究，并對不同循環次數下樁身水平位移進行了對比分析.

1 數值計算模型的建立

1.1 模型概況

以我國東海某近海海域風電場為例，建立了非均質土中海上風電單樁基礎數值計算模型，由于荷載和結構的對稱性，選擇半個物理模型的樁土體系進行建模，數值模型為半圓柱體. 為了忽略邊界條件對樁土體系的影響，在幾何模型上，用大尺寸來模擬半無限空間體，土體直徑取20D (D為樁徑)，土體高度取2hem(hem為樁的嵌固深度)，該尺度可滿足最小邊界尺寸和計算精度的要求[15?16].

ABAQUS中通過樁?土表面定義接觸屬性以模擬樁與土之間的剪力傳遞和相對位移，采用主?從接觸算法，選擇剛度大的樁體為主控面，土體表面為從屬面，樁?土法向行為采用硬接觸，切向行為采用摩爾?庫倫摩擦罰函數形式，界面滑動摩擦系數選取為土體內摩擦角)[17]，接觸對采用面對面接觸與有限滑移. 邊界條件是約束斷面處y方向位移，約束模型側面x和y方向位移，模型底端為固定約束. 樁體和土體都采用8節點6面體線性減縮積分三維實體單元(C3D8R)，采用掃掠的方式對網格進行劃分，從而使網格更加合理. 為了減小計算誤差，同時也為了縮短計算時間，采用為邊布種的方式，在樁土接觸面附近單元網格劃分得較細，而在遠離接觸面的土體，網格劃分的相對稀疏. 在實際狀況中，假定海床泥面處位移為零，土體內部是有應力存在的，因此在施加水平荷載前必須進行初始地應力平衡[18?19]. 本文采用ODB導入法進行初始地應力的平衡. 非均質土中海上風電單樁基礎數值計算模型如圖1所示.

圖1 非均質土中海上風電單樁基礎數值計算模型Fig.1 Numerical model of offshore wind power monopile foundation in heterogeneous soil

為了模擬樁?土之間的非線性，根據土層性質采用基于Mohr-Coulomb破壞準側的理想彈塑性本構模型來模擬，孔位學等[20]認為在非關聯流動法則條件下采用剪脹角 ψ =φ/2所得到的滑移線場與Prandtl理論一致. 樁周土及樁端土參數列于表1.樁體采用線彈性模型來模擬，樁的物理力學參數列于表2，其中15 m位于水中，1 m位于水面之上.

表1 樁周土及樁端土參數Table 1 Parameters of soil around pile and soil at pile end

表2 樁的物理力學參數Table 2 Physical and mechanical parameters of pile

1.2 荷載的確定與施加

為了建模方便和獲得較強的規律性，采用的分析步為靜力通用分析步，將海上風電單樁基礎在工作期間受到的波浪、洋流、風等荷載等效成雙向對稱循環荷載[21?23]的形式來模擬非均質土中海上風電單樁基礎動力響應特性，在模型中沿x軸方向不斷循環加載. ABAQUS中采用周期型幅值曲線來定義水平循環荷載[24]，周期型幅值曲線用傅里葉（Fourier）級數表示.

t≥t0時，幅值表達式為：

式中： ω 為圓頻率，其值為 ω =2πf，f為 頻率，t0為起始時間，A0為初始幅值，An為cos項的系數（n=1,2,3 , ···,N），Bn為 sin 項的系數.

1.3 樁側土抗力分析

樁側土抗力分布如圖2所示. 由圖（a）可知，水平荷載施加之前，初始土壓力沿樁周均勻分布，初始有效壓力值為 σ0. 荷載施加之后，土壓力沿樁周分布形式出現了變化. 由圖2（b）可知，施加水平荷載之后，樁周順載側有效土壓力增大到 σmax，樁周逆載側有效土壓力減小至 σmin. 陰影部分Ⅰ區為樁周逆載側減小的土抗力，陰影部分Ⅱ區為樁周順載側增加的土抗力. 施加水平荷載之后，作用在樁身的凈土抗力Pnet為Ⅰ區和Ⅱ區之和[25?26].

2 數值計算結果分析

在模型中定義參考點，將參考點與水面處樁身橫截面建立分布耦合約束，在參考點上施加水平循環荷載. 由于不同循環時樁身位移、樁身剪力、樁身彎矩和樁側土抗力沿埋深隨著時間的變化規律一致，循環的次數增加必然造成計算時間的增大，且對不同時間點的規律變化幾乎無影響，循環次數較大時意義不大，為了節約計算時間，且不影響其規律分析，故本文對第20次循環時水平循環荷載作用下樁身位移、剪力、彎矩、樁內外壁土抗力沿埋深變化特性進行了研究，并對不同循環次數下樁身位移進行了對比分析.

2.1 水平極限承載力的確定

由于海上風機單樁基礎樁徑較大且采用鋼管樁，樁身強度非常大[2]，且樁的水平荷載?位移曲線為緩變型，該曲線沒有明顯的拐點，可以認為是漸進式破壞，所以海上風機單樁基礎的水平極限承載力主要由樁體的水平變形控制. 采用位移控制法，對樁頂施加0.3 m的水平位移，在有限元軟件ABAQUS后處理中提取水平支反力和樁身水平位移的相關數據，然后利用繪圖軟件Origin繪制水平支反力與水平位移之間的關系曲線，得到樁頂荷載位移曲線，如圖3所示.

圖3 樁頂荷載位移曲線圖Fig.3 Load?displacement curve of pile top

根據允許變形法[9]，將樁頂水平位移達到0.02D時對應的水平荷載，確定為海上風機單樁基礎的水平極限承載力. 根據荷載位移曲線，當位移為0.1 m時對應荷載1.31 MN，為了反映單樁基礎所受的荷載水平，單樁上施加的水平循環荷載幅值的大小為其水平極限荷載Fu的一定比值，為了避免荷載幅值過小或過大對分析結果造成影響，故取水平荷載幅值為0.6Fu，根據工程資料，荷載頻率取為0.1 Hz進行研究.

2.2 樁身位移

對第20次循環時前5 s和后5 s樁身位移沿埋深變化規律進行分析，第20次循環時樁身水平位移變化曲線如圖4所示.

由圖 4（a）可知，191～195 s泥面處樁身水平位移分別為 13.7、25.3、26.4、17.2和1.83 mm，隨著時間的增加，位移先增加后減小，且位移增加的幅度越來越小，減小的幅度越來越大，由于荷載呈對數形式增加—荷載逐漸增大且增大的幅度逐漸減小，造成樁身位移增加幅度變小，193 s時樁身泥面處位移最大，屈曲現象最明顯，由于193 s之后荷載呈指數形式減小—荷載逐漸減小且減小的幅度逐漸增大，故樁身位移減小幅度增加. 193 s時與192 s時荷載值相同，但泥面處位移值比192 s時增加了 1.1 mm，194 s時與 191 s時荷載值相同，泥面處位移值比191 s時增加了3.5 mm，正向卸載與正向加載相同荷載值下樁身位移增大，說明循環荷載作用下樁身位移是不斷累積的，不同時間點的樁身位移零點位于泥面下32～34 m，樁繞樁身軸線上某一點轉動，表現出剛柔樁的性質.

由圖 4（b）可知，196～200 s泥面處樁身水平位移分別為?13.7、?25.7、?26.8、?17.5和?1.82 mm，由于荷載和結構的對稱性，負向荷載與正向荷載作用時樁身位移規律一致，關于y軸大致呈對稱分布；196～200 s時樁身位移零點主要位于泥面下34～36 m范圍內，樁身位移零點相比正向加載時發生了下移，這可能是由于循環過程中樁周土體軟化引起的.

圖4 第 20 次循環時樁身水平位移變化曲線. （a）前 5 s；（b）后 5 sFig.4 Horizontal displacement variation curves along the pile shaft during the 20th cycle: (a) first 5 s; (b) next 5 s

2.3 樁身剪力

對第20次循環時前5 s和后5 s樁身剪力沿埋深變化規律進行分析，第20次循環時樁身剪力變化曲線如圖5所示.

圖5 第 20 次循環時樁身剪力變化曲線. （a）前 5 s；（b）后 5 sFig.5 Shear force variation curves along the pile shaft during the 20th cycle: (a) first 5 s; (b) next 5 s

由圖5（a）可知，隨著時間的增加，樁身剪力先增加后減小，且剪力增加的幅度越來越小，減小的幅度越來越大，193 s時泥面處樁身剪力最大，正向卸載與正向加載相比相同荷載值下剪力是增加的. 由于樁周土體的抵抗作用出現反向的剪力值，不同時間點的樁身剪力反彎點位于泥面下7～9 m范圍內. 191～192 s為正向加載過程，樁身剪力零點下移，193～195 s為正向卸載過程，樁身剪力零點逐漸上移. 不同時間點的剪力沿x軸正方向最大值均位于泥面處，沿x軸負方向最大值均位于泥面以下31～32 m范圍內，此范圍內應適當的增加壁厚，提高抗剪強度.

由圖5（b）可知，負向加載時樁身剪力曲線與正向加載時樁身剪力曲線關于y軸大致呈對稱分布，規律類似.

2.4 樁身彎矩

對第20次循環時前5 s和后5 s樁身彎矩沿埋深變化規律進行分析，第20次循環時樁身彎矩變化曲線如圖6所示.

圖6 第 20 次循環時樁身彎矩變化曲線. （a）前 5 s；（b）后 5 sFig.6 Bending moment variation curves along the pile shaft during the 20th cycle: (a) first 5 s; (b) next 5 s

由圖6（a）可知，樁身彎矩由泥面處開始沿埋深先增大后減小，隨著時間的增加，彎矩先增大后減小，且彎矩增加的幅度越來越小，減小的幅度越來越大. 不同時間點下樁身彎矩均沒有出現反彎點，最大彎矩均位于泥面下7～9 m范圍內，說明此范圍內最容易發生破壞，在實際的工程設計時，此處樁身的壁厚或強度應適當增加，以減小樁體在該位置的轉角. 隨著時間的增加樁身最大彎矩對應的樁的埋深位置先增大后減小，與剪力零點表現出來的規律一致.

由圖 6（b）可知，后 5 s樁身彎矩曲線與前 5 s關于y軸大致呈對稱分布，規律類似.

2.5 樁側土抗力

對第20次循環時前5 s和后5 s樁外壁土抗力沿埋深變化規律進行分析，第20次循環時樁外壁土抗力沿埋深分布曲線如圖7所示.

圖7 第 20次循環時樁身外壁土抗力埋深分布曲線. （a）前 5 s順載側；（b）后 5 s逆載側Fig.7 Soil resistance distribution curves of pile outer wall along the buried depth during the 20th cycle: (a) forward side of the first 5 s; (b) reverse load side after 5 s

由圖7（a）可知，不同時間點樁順載側土抗力曲線隨時間的變化存在分界點. 分界點以上樁側土抗力隨著時間的增加先增大后減小，這種變化主要集中在埋深20 m范圍內，這是因為此范圍內樁身位移為正，樁與外部土體相互擠壓，樁外壁土壓力為被動土壓力；在193 s時達到最大，說明此范圍內隨時間的增加樁身承擔荷載變化比較明顯. 分界點以下樁側土抗力隨著時間的增加先減小后增大，這種變化主要集中在距樁底端4 m范圍內，這是由于此范圍內樁身位移為負，樁與外部土體產生分離的趨勢，樁外壁土壓力為主動土壓力；在193 s時達到最小，說明此范圍內隨著時間的增加對樁側土抗力影響較大. 在淤泥和粉砂土分界面處土抗力增加顯著，究其原因，該處由于土層之間性質差異較大發生了明顯的應力集中現象.

軟黏土中樁外壁土抗力沿埋深逐漸增大，這是因為軟黏土水平位移沿樁埋深逐漸增大，樁身需要提供較大的支撐力去約束較大的位移；而砂土中樁外壁土抗力沿埋深逐漸減小，究其原因，一方面是由于淤泥土層與粉砂土層相比土層性質較差，另一方面是由于淤泥土與粉砂土相比水平位移較大，且由于靠近樁端，粉砂土水平位移沿樁埋深逐漸減小，故樁身提供的支撐力沿埋深逐漸減小，表現為樁側土抗力沿埋深逐漸降低. 樁外壁土抗力沿埋深變化不一致，說明樁順載側外壁土抗力的發揮不僅僅取決于樁土相對位移，還與埋深和樁側土層性質有關.

由圖7（b）可知，不同時間點樁側土抗力發揮沿埋深不同，198 s時土抗力發揮位置埋深最大，這主要是因為此時樁身負向位移最大，與逆載側土體脫離程度比較大，導致樁側土抗力發揮沿埋深滯后. 不同時間點樁逆載側土抗力沿埋深存在分界點，樁逆載側土抗力隨時間變化規律與順載側土抗力正好相反. 樁逆載側土抗力沿埋深變化與樁身順載側外壁土抗力變化大致相同. 上部土層土抗力較小，下部土層土抗力較大，承擔的荷載較大，可作為持力層.

對第20次循環時前5 s和后5 s樁內壁土抗力沿埋深變化規律進行分析，第20次循環時樁內壁土抗力沿埋深分布曲線如圖8所示.

圖8 第 20次循環時樁身內壁土抗力埋深分布曲線. （a）前 5 s順載側；（b）后 5 s逆載側Fig.8 Soil resistance distribution curves of pile inner wall along the buried depth during the 20th cycle: (a) forward side of the first 5 s; (b) reverse load side after 5 s

由圖 8（a）可知，191～195 s樁側土抗力曲線沿埋深幾乎重合，承擔的荷載基本不變，在淤泥質黏土層中土抗力呈指數型增加，在淤泥土層中土抗力呈對數型增加. 粉砂土中土抗力先增大后減小，在距樁底端1/3粉砂土層層厚處土抗力值達到最大，說明樁內部土體繞此處發生旋轉，使得此處以上樁側土抗力增加，而以下由于樁與土脫離，而造成樁側土抗力減小.

由圖8（b）可知，樁逆載側內壁土抗力沿埋深分布規律與順載側內壁土抗力基本一致，僅在距樁底端2 m范圍內有所不同.

2.6 不同循環次數下樁身位移

為了研究水平位移隨時間的變化規律，對不同循環時泥面處和樁底端水平位移時程曲線進行了對比分析，第N次循環時泥面處樁身和樁底端水平位移變化曲線如圖9所示.

圖9 第 N 次循環時泥面處樁身和樁底端水平位移. （a）泥面處；（b）樁底端Fig.9 Horizontal displacement of pile shaft at mud surface and pile bottom at the Nth cycle: (a) at the mud level; (b) bottom end of pile

由圖9（a）可知，前3 s泥面處樁身位移隨著循環次數的增加沿y軸逐漸下移，位移絕對值先減小后增加，3～5 s位移全部為正，而 5～8 s泥面處樁身位移隨著循環次數的增加沿y軸逐漸上移，后2 s位移全部為負，由于結構和荷載的對稱性，正向荷載作用時和負向荷載作用時樁身水平位移規律類似. 不同循環時泥面處樁身水平位移均呈現出周期性的變化.

從圖9（a）還可以看出，第20次和第100次循環時位移最大正值和最小負值分別發生在3和8 s，在第200次和第300次循環時在4和9 s分別產生最大正值和最小負值，當循環次數為400和500次時，5 s位移最大，第 10 s位移最小，說明隨著循環次數的增加，產生最大水平位移的時間點滯后. 不同循環次數下正向位移最大值分別為26.41、27.51、26.28、27.93、26.8和 28.32 mm，負向位移最大值分別為?26.79、?27.85、?26.46、?28.02、?26.72和?28.14 mm，同一時刻下發生的樁身最大位移值隨循環次數的增加而增大，說明樁身最大位移不斷累積. 由圖9（b）可知，樁底端水平位移與泥面處樁身水平位移規律類似.

第N次循環結束時泥面處樁身和樁底端水平位移變化曲線如圖10所示.

圖10 第 N 次循環結束時泥面處樁身和樁底端水平位移. （a）泥面處；（b）樁底端Fig.10 Horizontal displacement of pile body at the mud surface and pile bottom at the end of the Nth cycle: (a) at the mud level; (b) bottom end of pile

由圖10（a）可知，不同循環結束時泥面處樁身水平位移均為負，分別為?1.82、?8.32、?16.31、?22.71、?26.72和?28.14 mm，隨著循環次數的增加，泥面處樁身水平位移絕對值逐漸增大，且增加的幅度越來越小，說明樁周土體強度隨循環次數的增加而逐漸弱化，引起樁及樁周土體變形較大，且隨著循環次數的增加變形逐漸趨于穩定.

由圖10（b）可知，不同循環結束時樁底端位移均為正，分別為 0.13、0.52、0.97、1.33、1.57和1.65 mm，隨著循環次數的增加，樁底端水平位移逐漸增大，且增加的幅度逐漸變小.

3 數值計算結果對比分析

由于模型試驗中樁的直徑相對較小，無法直接與本文數值模型中的大直徑樁進行對比驗證，故與已有的有限元分析結果進行對比分析. 文獻[27]中陳新奎采用ABAQUS對大直徑樁的水平受荷特性進行了模擬，并將數值模擬結果與風電場現場試樁結果進行了對比分析，驗證了ABAQUS數值模擬的可靠性. 取文獻[27]中4 MN（此荷載與本文荷載幅值最接近）水平力作用下樁身水平位移、剪力、彎矩與193 s時樁身水平位移、剪力、彎矩（此時樁身水平位移、剪力、彎矩最大）進行對比分析，如圖11～13所示.

從圖 11～13可以看出，圖（a）與圖（b）中泥面以下范圍內樁身水平位移、剪力、彎矩沿埋深變化規律類似. 由圖11可知，樁身水平位移曲線沿埋深均出現零點，均呈現出非線性的變化，樁均繞樁身軸線上某一點轉動，表現出剛柔樁的性質. 由圖12可知，泥面以下范圍內樁身剪力均出現反彎點，圖12（a）中樁身剪力變化幅度較大，究其原因，樁身水平位移較大，引起樁周土體壓縮范圍較大，能夠提供較大的土抗力[28?29]，導致剪力變化程度較大. 由圖13可知，泥面以下范圍內樁身彎矩沿埋深均沒有出現反彎點，均沿埋深先變大后變小，最大值均發生在淺層土體，樁底端彎矩均接近于零.

圖11 樁身水平位移對比分析. （a）4 MN 水平力作用下樁身位移[27]；（b）193 s 時樁身位移Fig.11 Comparative analysis of horizontal displacement along the pile shaft: (a) pile displacement under 4-MN horizontal force[27]; (b) pile displacement at 193 s

圖12 樁身剪力對比分析. （a）4 MN 水平力作用下樁身剪力[27]；（b）193 s時樁身剪力圖Fig.12 Comparative analysis of shear force along the pile shaft: (a) pile displacement under 4-MN horizontal force[27]; (b) pile displacement at 193 s

圖13 樁身彎矩對比分析. （a）4 MN 水平力作用下樁身彎矩[27]；（b）193 s 時樁身彎矩圖Fig.13 Comparative analysis of bending moment along the pile shaft: (a) pile displacement under 4-MN horizontal force[27]; (b) pile displacement at 193 s

綜上可知，本文運用ABAQUS建立的數值模型是可靠的，且數值模擬結果是準確的.

4 結論

（1）樁繞樁身軸線上某一點轉動，表現出剛柔樁的性質；不同循環時泥面處樁身水平位移均呈現出周期性的變化；隨著循環次數的增加，泥面處樁身最大位移發生的時間點滯后，樁身水平位移絕對值逐漸增大，且增加的幅度越來越小，逐步趨于穩定，樁身水平位移隨時間變化逐漸累積.

（2）由于樁周土體的抵抗作用出現反向剪力值，不同時間點的樁身剪力反彎點位于埋深7～9 m范圍內，隨樁身剪力絕對值的增大而沿埋深有所下降；不同時間點剪力沿x軸正方向最大值均位于泥面處，沿x軸負方向最大值均位于埋深31～32 m 范圍內.

（3）不同時間點樁身彎矩均沒有出現反彎點，彎矩最大值均位于泥面下7～9 m范圍內，與樁身剪力零點位置相對應，發生在淺層土體；不同時間點樁底端的彎矩接近于零，說明樁底部土體對樁的嵌固作用明顯.

（4）淤泥土和粉砂土分界面處由于土層之間性質差異較大，導致樁外壁土抗力增加顯著；樁身外壁土抗力沿埋深曲線隨時間的變化出現分界點，分界點上下范圍內土抗力隨時間變化規律正好相反；樁身外壁土抗力的發揮不僅僅取決于樁土相對位移，還與埋深和樁側土層性質有關；不同時間點樁身內壁土抗力曲線沿埋深幾乎重合，承擔的荷載基本不變.