長期循環荷載下海上單樁基礎變形特性的數值模擬研究

楊鈺榮，孫毅龍，許成順，翟恩地

（1.北京工業大學城市與工程安全減災教育部重點實驗室，北京 100124；2.溫州大學建筑工程學院，浙江溫州 325035；3.金風科技股份有限公司，北京 100176）

我國近海岸以大規模的砂質沉積體為主，大直徑單樁基礎是海上風電系統常用的基礎形式之一。海上風電單樁基礎在服役期間，長期承受風荷載、浪荷載等引起的循環傾覆力矩，荷載循環作用次數可達到108次[1]。單樁基礎長期承受這種低頻的循環荷載作用，會導致樁周土體剛度發生變化，造成樁基側向累積變形增大，嚴重危及海上風電系統的正常運行。因此，海上風電單樁基礎在長期荷載作用下的累積變形是樁基設計中重點考慮的因素之一。

目前，海上風電單樁基礎累積變形特性的研究主要采用修正p-y曲線法、室內或現場試驗和數值模擬分析方法等。MCCLELLAND B等[2]于1958年基于溫克爾地基梁模型提出了計算水平荷載作用下單樁非線性變形的方法，即p-y曲線法。2006年，LESNY K等[3]發現美國石油學會（American Petroleum Institute，API）推薦的p-y曲線法高估了地基初始剛度，針對該問題對p-y曲線法進行了簡單的修正。HYUNSUNG L等[4]于2018年發現荷載作用后地基剛度變化影響風機系統頻率，在不同的加載頻率下，利用模型試驗建立擬靜力分析的簡化動力p-y曲線。2020年，胡安峰等[5]基于剛度衰減模型對p-y曲線進行修正，綜合考慮了荷載特性、循環次數、砂土剛度循環弱化等因素。2021年，王衛等[6]針對p-y曲線過于保守等問題，提出一種綜合考慮樁徑和地層深度影響的黏土p-y修正模型。盡管修正p-y曲線法在一定程度上考慮了土體條件、樁基特點和荷載環境對預測樁基累積變形的影響，但考慮荷載的循環作用次數與實際樁基受荷次數相差較大，對于預測大直徑樁基受長期荷載作用的循環累積變形相對保守。

為了進一步研究海上風電單樁的循環累積變形問題，不少學者基于室內外模型試驗建立循環累積變形與靜態荷載作用下的變形、循環次數的經驗函數關系，提出不同類型的水平循環受荷樁的累積變形預測模型。VERDURE L等[7]、BIENEN B等[8]、LI W等[9]基于水平循環加載模型試驗研究了長期荷載作用下單樁的動力響應特性，認為樁基的循環累積變形與樁基初始變形呈對數關系。LEBLANC C等[10]、RASMUS等[11]、吳金標等[12]通過開展離心機模型循環加載試驗，認為冪函數模型能夠較好地預測砂土場地中水平受荷樁的循環累積變形。

由于試驗條件、荷載特性等差異，上述學者所建立的經驗預測模型不盡相同，對數形式和冪函數形式的模型預測結果存在差異。為此，眾多學者基于循環試驗，通過建立樁—土相互作用數值模擬，研究長期荷載作用下海上風機單樁基礎的累積變形特性。ACHMUS M等[13]于2008年基于循環三軸試驗結果提出長期荷載作用下的剛度衰減模型，利用該模型預測N次循環荷載作用下單樁累積側向變形，并與HETTLER A等[14]、LITTLE R L等[15]分別提出的對數函數型和冪函數型樁基累積變形預測模型進行對比；2015年，DEPINA I等[16]將三維有限元樁土模型與剛度退化模型耦合，模擬了砂土中單樁基礎在長期循環荷載作用下的變形、轉角和彎矩響應；羅如平等[17]、董愛民等[18]分別于2016年、2017年利用剛度衰減模型，分析長期荷載作用下大直徑單樁的循環累積變形特性，并對現有的樁基累積變形預測模型進行修正；2018年，YANG M[19]基于Martin Achmus提出的剛度退化模型運用三維有限差分模型，研究了單樁基礎的長期循環特性，定量地分析了循環荷載特性、樁基埋置深度、荷載偏心率對單樁長期累積變形的影響；2020年，LUO R等[20]基于剛度衰減模型建立了修正的樁基累積變形預測方法，開展不同樁徑、荷載等對樁基累積變形的影響規律研究。

綜上所述，剛度衰減模型是分析單樁基礎累積側向變形較為簡便的方法。但目前該模型參數的確定尚存在困難。為此，本文利用豎向—扭轉雙向耦合循環剪切試驗（模擬波浪等循環荷載作用）結果獲得剛度衰減模型，并基于Flac3D有限差分平臺開展三維樁—土相互作用數值模擬，分析循環荷載幅值、場地密實度和樁徑等對長期荷載作用下單樁基礎累積變形的影響規律；根據數值計算結果，對常用的冪函數型和對數函數型樁基累積變形預測模型的適用性進行評估，為海上風電樁基設計提供參考。

1 剛度衰減模型

1.1 剛度衰減模型介紹

HUURMAN M[21]基于室內循環三軸試驗指出砂土在長期循環荷載下其模量會發生衰減，并將土體的累積塑性應變增量比改寫為土體的割線模量增量比，具體表達式如式（1）所示，示意圖如圖1所示。

圖1 砂土剛度衰減模型示意圖

式中，EsN為循環荷載作用N次后的割線模量；Es1第一次循環荷載作用后的割線模量；εp，N=1為第一次循環荷載作用后的應變；εp，N為循環荷載作用N次后的應變。

ACHMUS M等[13]基于Huurman提出的土體割線模量衰減公式，建立了剛度衰減模型來表征不同循環次數下的樁基累積變形發展規律，具體公式如下。

式中，N為荷載循環次數；X為定義的循環應力比；b1、b2是根據試驗得到的回歸參數；C是由參數b1、b2和循環應力比X控制的指數。

由于在實際海上風電樁—土相互作用系統中，土體存在非均勻初始應力，即使未施加循環應力，也會存在初始的非零循環應力比，使式（3）得到一個不符合現實的應變增量，對地基模量計算存在影響。考慮非均勻初始應力的影響，定義特征循環應力比Xc，在數值模型計算中用以替換X。具體表達如下。

式中，X(1)為受荷狀態下土體的循環應力比；X(0)為初始應力狀態下土體循環應力比；σ1為土體最大主應力；σ3為土體最小主應力；σ1，sf為土體在靜態力下破壞時的主應力；σ1，cyc為施加循環力時土體的最大主應力；φ為土體內摩擦角。

1.2 土體剛度衰減特性研究

本文利用豎向—扭轉雙向耦合剪切儀[22]開展相對密實度Dr=35%、Dr=50%、Dr=70%的飽和砂土（福建標準砂）循環排水剪切試驗，得到不同相對密實度砂土在長期循環荷載作用下的模量衰減規律。為了模擬實際海床樁周土的受力狀態，在豎向—扭轉耦合循環剪切試驗中，對試樣同時施加軸向荷載及扭矩，并保持二者的相位差為90°；設置固結應力比Kc=2，試樣初始固結圍壓σ1=100 kPa、σ2=σ3=50 kPa。海上風機樁基周圍土體主要受風、波浪荷載的影響。常見波浪和風荷載的典型頻率為0.1 Hz[23]。砂性土在這種低頻的循環荷載作用下易排水。因此，循環試驗中的加載頻率為0.1 Hz。通過簡單的數值計算，選取樁周土的循環應力比X為0.10、0.12、0.15，試驗方案見表1。

表1 試驗方案

利用式（1）計算每一次循環荷載作用后的割線模量，得到土體的割線模量隨循環加載次數的變化規律，如圖2所示。由圖可知，割線模量在循環荷載作用前100個周期內衰減最為顯著，此后逐漸衰減至某一數值，且與荷載幅值沒有明顯的依賴關系。

圖2 不同相對密實度砂土的割線模量衰減規律

根據前文介紹的剛度衰減模型，基于式（4）對不同相對密實度砂土的割線模量進行擬合，得到割線模量對數比與循環加載次數對數之間的關系，如圖3所示。由圖可知，兩者之間呈近似的線性關系，且不同相對密實度砂土的割線模量衰減曲線斜率C不同。此外，由式（5）可知，斜率C由循環荷載影響指數b1、b2和循環應力比X控制。

圖3 不同相對密實度砂土的割線模量衰減比

為了將循環試驗結果應用于樁基的循環累積變形分析，基于試驗結果利用式（6）對參數b1、b2進行反算，得到b1、b2和土體相對密實度Dr之間的關系，如圖4所示。由圖可知，參數b1、b2與相對密實度Dr呈線性關系式，b1、b2與Dr擬合公式如下（R2>0.99）。

圖4 參數b1、b2與砂土相對密實度關系曲線

1.3 剛度衰減模型介紹

應用FLAC3D有限差分平臺內置的FISH語言，根據試驗結果編制相應的剛度衰減程序，并建立樁—土相互作用整體數值分析模型。將循環荷載作用下土體模量的衰減規律嵌入的數值分析計算模型中，實現土體模量隨應力狀態及循環次數的更新。在整體數值分析模型中分3個步驟實現樁基累積變形的計算，示意圖如圖5所示。剛度衰減模型的參數b1、b2根據實際場地的土體密實度，按照式（11）、式（12）計算取值。根據圖5中的步驟，更新土體模量后計算的樁基側向變形即為不同循環次數荷載作用后樁基的循環累積變形。

圖5 樁基累積變形計算步驟圖

2 數值模型驗證

基于上述剛度衰減數值分析模型，模擬ACHMUS M等[24]開展的離心機模型試驗及張紀蒙等[25]開展的水平加載模型試驗；并與當前兩種主要形式的樁基累積變形預測模型計算結果進行對比，具體樁基累積變形預測模型如表2所示。

表2 樁基累積變形預測模型

2.1 數值模型的建立

數值模型中，土體采用摩爾—庫倫本構模型，同時將式（9）土體模量與小主應力之間的關系引入到數值模型中，考慮土體彈性模量沿深度的非線性變化，土體模量非線性增長公式如下。

式中，σat為大氣壓力；σm為土體單元平均應力；κ與λ為土體剛度參數。

樁體采用彈性模型，樁基的建模采用等效剛度法。樁—土之間的相互作用采用接觸單元來模擬，接觸單元的法向和切向剛度的取值方法按照式（14）計算。

式中，K、G為接觸面兩側材料體積模量和剪切模量；ΔZmin為接觸面兩側單元最小尺寸。

按上述各式和模型尺寸計算得到kn=ks=1 012 Pa/m[26]，樁—土接觸單元參數的粘聚力和內摩擦角取相鄰土體0.3~0.5倍[27]。樁基和土體的具體參數見表3、表4。

表3 ACHMUS M等[24]模型試驗參數

表4 張紀蒙等[25]模型試驗參數

樁—土相互作用數值計算模型水平方向上的土體長度取為12倍樁徑，樁端底部的土體模型長度取為10倍樁徑；保證在該模型尺寸下樁基性能不受邊界的影響；為提高計算效率，在樁基附近劃分較密的網格，外側網格相對疏松；以ACHMUS M等[24]的樁基試驗數值模型為例，如圖6所示。

圖6 ACHMUS M等[24]樁基試驗數值計算模型

2.2 樁基累積變形計算對比

圖7、圖8分別為ACHMUS M等[24]、張紀蒙等[25]模型試驗中的樁基累積變形計算結果。本文數值模型計算的樁基累積位移與模型試驗實測值基本吻合，較好地反映了單樁基礎長期受荷載作用的累積特性，可以用來分析長期循環荷載作用下海上風機單樁基礎的循環累積變形。

圖7 ACHMUS M等[24]模型試驗對比

圖8 張紀蒙等[25]模型試驗對比

此外，兩種預測模型計算的樁基累積變形與試驗實測值的趨勢基本一致；冪函數模型的預測結果與實測值更加接近；而對數型模型的預測結果在循環次數較大時，與試驗實測值有一定偏差。

3 樁基累積變形特性分析

3.1 數值模型建立

基于前文建立的數值模型，考慮中密、密實、緊密砂土3類場地中，鋼管樁樁長L=30 m，樁徑D=4 m、5 m、6 m、7 m、8 m的單樁基礎，樁基受水平側向荷載與泥面高1 m，數值模型如圖9，具體的樁、土參數如表5、表6所示。

圖9 單樁基礎三維數值模型

表5 數值模型土參數[28]

表6 數值模型樁參數

3.2 樁基水平極限承載力

為了確定數值模型中對樁基施加的荷載大小，先計算靜力荷載下樁基水平極限承載力，通過引入無量綱荷載幅值系數φb以確定不同工況下的水平荷載大小。荷載幅值系數φb的表達式如式（15）所示。

式中，Hmax為施加的水平靜力側向荷載；Hu為計算的樁基水平極限承載力。

樁基水平極限承載力不是樁基破壞時的承載力，而是被定義為發生某一變形值時的水平荷載。AHMED S S等[29]、BARARI A等[30]提出單樁泥面水平位移為0.1D時所受的水平力為極限承載力。本文采用該極限承載標準，同時計算了不同條件下單樁基礎的水平極限承載力。數值計算結果如圖10所示。

圖10 不同樁徑樁基水平極限承載力

當場地土體的密實程度相同時，樁基水平極限承載力隨樁徑增加而增大，且二者呈現近線性增長的關系；這是由于隨著樁徑的增大，土體的破壞模式逐漸向半剛性樁，剛性樁的破壞模式發展，即土體受力區范圍增大，因而樁基的水平承載力增加。此外，當樁徑相同時，樁基的水平極限承載力隨著場地密實度的增加而增大；例如，當D=5 m時，密實砂土、緊密砂土場地中樁基水平極限承載力分別為中等密實砂土場地中的1.25倍、1.52倍；這是因為土體的強度、剛度與密實度緊密相關，砂土越密實其承載能力越強，因而樁基極限承載力越大。

3.3 荷載幅值的影響

以中等密實場地中D=4 m的工況為例，計算樁基不同荷載大小下樁基的累積變形。設置4組不同的荷載幅值，φb=0.10、0.20、0.25、0.30，最大循環次數N=106，該工況下的樁基水平極限承載Hu=26.8 MN。

以荷載幅值系數φb=0.20、0.25、0.30的3種工況為例，繪制樁身側向變形隨循環次數變化曲線，如圖11所示。由圖可知，當樁基受靜態側向力作用時（N=1），荷載幅值系數越大，樁基的初始側向變形越大；隨著循環次數的增加，樁基側向變形不斷增長，且荷載幅值系數越大樁基位移累積越快；此外，在不同大小的荷載長期作用下，樁基轉動中心始終位于泥面下24.8 m處，約為0.7倍的樁基埋深。

圖11 不同荷載幅值下樁身位移發展

為了對比不同條件下樁基累積變形隨循環荷載作用次數增加的發展規律，評價樁基累積變形預測模型的可靠性，對樁基泥面位移進行無量綱處理。根據數值模型計算結果，得到樁基泥面位移累積率yN/y1與循環荷載次數N的關系。最后，對比分析冪函數型、對數函數型樁基變形預測模型和數值模型計算的樁基累積位移發展規律。不同循環荷載幅值下樁基泥面位移累積率的發展規律如圖12所示。

圖12 不同荷載幅值下樁基泥面位移發展

由圖12可知，樁基累積變形隨荷載作用次數先增加隨后逐漸趨于穩定，在循環加載前期（N<105），冪函數模型的預測結果偏小，與數值計算結果的誤差在18.88%以內，對數函數模型的預測結果與數值計算結果吻合較好，其最大誤差在13.61%以內；隨著循環加載次數不斷增加，當N>105時，冪函數模型的計算結果與數值計算值最大誤差僅為1.42%，其預測結果更加準確；而對數函數模型預測值的最大誤差為4.68%；隨著荷載幅值系數增大，冪函數模型預測結果的偏差逐漸減小，而對數函數模型的預測結果偏差逐漸增大。分析表明，從樁基累積變形的長期發展來看，在分析不同荷載大小作用下的樁基累積變形時，冪函數模型的預測結果更加準確。

根據式（16）、（17）反算得到不同荷載大小下的模型參數α與tb，分析冪函數和對數函數預測模型中參數α與tb受荷載幅值影響的規律。

式中，y1為樁基的初始側向位移；yN為循環荷載作用N次后樁基的側向位移。

圖13繪制了α和tb與φb的關系曲線。由圖可知，參數α和tb均與荷載幅值系數φb呈正比。此外，對數函數模型參數擬合曲線的斜率大于冪函數模型參數擬合曲線的斜率，這說明對數函數預測模型的準確性受荷載幅值的影響更大，與圖12所得結論一致，預測不同荷載幅值下樁基的循環累積變形規律時，冪函數預測模型的計算結果較為準確。

圖13 模型參數與荷載幅值系數的關系

3.4 場地密實度的影響

為研究場地土體密實度對樁基循環累積變形的影響，考慮一樁徑D=5 m的單樁設置于中密、密實和緊密砂土3種場地中，在樁頂施加循環的水平側向荷載，分析樁基長期荷載作用下樁基累積變形規律。樁基的循環累積變形計算屬于疲勞承載狀態，基于DNV（Det Norske Veritas）規范，樁基的循環荷載幅值取為水平極限承載力的25%，即φb=0.25[31-32]。水平荷載幅值大小如表7所示。

表7 不同砂土中樁基承受的荷載幅值

圖14為不同密實度砂土中樁身變形發展。由圖可知，3種場地中的樁基側向變形均隨荷載作用次數增加而增大；場地土體越密實，樁基受靜態荷載（N=1）作用后的初始位移越小，受循環荷載作用后的位移累積越小。此外，3種密實度砂土中樁基轉動中心均位于泥面以下24.8 m位置處，約為0.7倍的樁基埋深，這說明樁基轉動中心受場地密實度的影響較小。

圖14 不同密實場地中樁身位移發展

為分析不同密實場地中樁基累積變形預測模型的適用性，根據式（16）和式（17）反算得到不同密實度砂土中的模型參數α與tb，并對比中密、密實、緊密砂土中，兩種預測模型計算的樁基泥面位移累積率，如圖15所示。由圖可知，3種密實度砂土中樁基的泥面位移累積率在循環加載初期（N<105）增加顯著隨后逐漸趨于穩定，緊密砂土中樁基的位移累積率最小。

圖15 不同密實度下樁身泥面位移發展

兩種模型的預測結果顯示，當N<105時，冪函數模型的預測結果與數值計算結果的最大誤差約為32.42%，對數函數預測結果與數值計算結果的最大誤差約為27.47%。當N>105時，對數函數預測模型結果與數值計算結果的最大誤差為13.86%；冪函數模型的預測結果與數值結果的最大誤差為4.23%，其預測結果更為準確。由上述分析可知，在不同密實程度的砂土場地中，使用冪函數模型預測樁基累積變形的長期發展更具優勢。

為了分析土體密實程度對冪函數模型參數α和對數函數模型參數tb的影響規律，圖16繪制了α和tb與土體密實度Dr的關系曲線。由圖可知，冪函數模型和對數函數模型中的參數α和tb均隨場地土體密實度的增加而降低。這是因為場地土體越密實，樁基位移累積率越小。此外，α與Dr關系曲線的斜率接近為0，說明密函數預測模型受土體密實度的影響較小。

圖16 模型參數與土體相對密實度的關系

3.5 樁徑的影響

考慮樁徑D=4 m、5 m、7 m、6 m、8 m的單樁設置于中密砂土場地中，計算不同樁徑的樁基在長期荷載作用下的累積變形，不同樁徑樁基所施加的荷載大小如表8所示。

表8 不同樁徑樁基的荷載大小

圖17以樁徑D=4 m、6 m、8 m的樁基為例，繪制了長期荷載作用下樁身側向變形的發展。由圖可知，樁徑越大，樁基的側向變形越大；當樁基受靜力荷載作用時，3種樁徑樁基的泥面初始位移分別為5.3 cm、7.7 cm、9.1 cm；當N=106時，樁徑D=4 m、6 m、8 m的樁基側向變形分別增加了83.5%、203.3%、334%。此外，3種樁徑樁基的轉動中心分別位于泥面下23 m、24.5 m、25 m，樁基轉動中心隨樁徑的增大而呈現下移趨勢，但不隨循環次數的增加而變化。

圖17 不同樁徑單樁基礎的樁身側向位移發展

根據式（16）和式（17）反算得到樁徑不同時的模型參數α與tb，根據兩種預測模型計算不同樁徑樁基的位移累積率，將其計算結果與數模擬結果進行對比，如圖18所示。由圖可知，隨荷載作用次數的增加，不同樁徑樁基的泥面位移累積率均先增加后逐漸趨于穩定；且樁基的泥面位移累積率隨樁徑的增加而增大，當N=106時，較小直徑（D=4 m、5 m）樁基的泥面位移累積率逐漸穩定于某一值，而直徑較大（D=6 m、7 m、8 m）的樁基泥面位移累積率仍有進一步增加的趨勢。對比樁基累積預測模型的計算結果發現，在樁徑不同的條件下，冪函數型預測模型的計算結果與數值模型的結果在整個加載周期內的最大誤差在2.92%以內，與數值計算值基本一致；在樁徑較大的情況下（D=7 m、D=8 m），當N<105時，對數函數模型的預測結果與數值計算結果最大誤差為4.01%，存在明顯的偏差，在長期荷載作用后（N>105），對數函數模型的預測結果與數值計算結果的最大誤差為7.73%。分析表明，在計算較大直徑樁基的循環累積變形時，冪函數模型仍然具有較大優勢。

圖18 不同樁徑單樁基礎泥面累積變形

為進一步分析樁徑對預測模型中參數α和tb的影響規律，圖19分別繪制了兩種形式預測模型中的參數α和tb與樁徑D的關系。由圖可知，冪函數模型和對數函數模型中的參數α和tb均與荷載幅值系數呈正比，但兩條關系曲線的斜率不同，對數函數預測模型受樁基直徑的影響偏大。對于樁徑較大的樁基，利用對數函數模型預測樁基累積變形可能偏大。

圖19 模型參數與樁徑的關系

4 結 論

本文基于FLAC3D有限差分平臺，根據循環試驗結果編制剛度衰減程序，建立海上風機三維樁—土相互作用數值分析模型，研究單向水平循環荷載作用下荷載幅值、砂土相對密度和樁徑大小對單樁基礎累積變形特性的影響；根據整體數值分析模型的計算結果評價了冪函數和對數函數形式的樁基變形預測模型的適用性。相關結論如下。

（1）在長期水平循環荷載作用下，樁基的累積變形隨循環荷載幅值的增加而增大；隨場地砂土相對密實度的增加而降低；當無量綱荷載幅值系數相同時，樁基累積變形隨樁徑的增加而增大；總體上，樁基的累積變形隨循環次數的增加先增加后逐漸趨于穩定。

（2）當循環荷載幅值、場地密實度增大時，樁基轉動中心的埋深位置基本不變；當荷載無量綱系數相同時，樁基轉動中心位置隨樁徑的增大呈現下移趨勢。

（3）冪函數和對數函數模型的模型參數隨荷載幅值和樁徑的增大而增加，隨場地土體密實度的增加而減小。

（4）在循環荷載作用前期（N<105），對數函數模型預測樁基變形較為準確，冪函數模型的計算結果偏小；隨著循環荷載加載次數的增加，冪函數預測模型的預測結果更為準確；從樁基累積變形的長期發展來看，冪函數型預測模型更具有優勢。