桶形基礎二次加載的大變形有限元模擬

陳霄，鄭敬賓，王棟

(中國海洋大學 山東省海洋環境地質工程重點實驗室，山東 青島 266100)

近年來，負壓桶形基礎(以下簡稱為桶基)作為單樁基礎和導管架樁基的替代形式，已用于支撐多個海上風電場的風機和升壓站，如英國的Hornsea 3項目[1]。與風機的一次性吊裝不同，桶基支撐的升壓站平臺分兩個階段安裝：第1階段，桶基與平臺貫入土中預定深度，施工船返回港口，桶基靜置數周至數月；第2階段，各種大型電氣模塊運輸至場地，并吊放到平臺上，桶基實質上承受二次加載。第1階段的桶基貫入在黏土中引起超靜孔壓，之后孔壓部分或完全消散[2-3]，導致第2階段開始時桶基周圍黏土的不排水強度高于第1階段結束時。

類似于一般的淺基礎承載力問題，常規有限元方法被廣泛用于分析不排水條件下桶基的豎向承載力及復合承載力[4-6]，但這些研究均假定桶基預埋就位，不考慮桶基安裝過程對承載力的影響。這是因為桶基長距離貫入導致側壁端部周圍土體單元的嚴重扭曲，常規有限元方法無法捕捉這一大變形過程。

研究者借助離心機試驗[7-11]或大變形有限元模擬[12-13]再現了桶基的貫入，探究了土層分布、安裝方式、基礎尺寸等一系列影響桶基貫入阻力的因素，提出了確定第1階段貫入阻力的方法。少量的數值模擬還進一步給出了桶基安裝后的孔壓消散過程[3]，但這些結果僅評估了地基土的固結程度，不能預測升壓站平臺第2階段的抗壓承載力。Jeanjean[14]在墨西哥灣海域進行了現場試驗，并總結了歷史上該區域的施工記錄，得到桶基安裝就位時的貫入阻力以及固結一段時間后的抗拔力，發現桶基在安裝后存在顯著的固結效應，最終抗拔力最高可達初始貫入阻力的1.6倍，但其試驗針對的是抗拔力，而不是抗壓承載力。

鑒于與桶基二次加載有關的研究較為有限，筆者發展了有效應力形式的大變形有限元方法，模擬黏性土中桶基的“貫入-固結-二次加載”完整過程，探索桶基初始貫入和部分排水固結后的超靜孔壓分布，進而考察固結效應對二次加載的影響，合理評估二次加載的豎向承載力。

1 研究方法

1.1 大變形有限元方法

桶形基礎的安裝是典型的巖土工程大變形問題，因此使用一種名為RITSS的大變形有限元方法(Remeshing and interpolation technique by small strain)。這種方法最初由Hu等[15]提出。Wang等[16-17]將其與有限元商業軟件Abaqus結合并擴展至有效應力分析。RITSS方法中，整個計算過程被劃分為一系列小變形計算步，通過網格重劃分及場變量映射，避免了單元的嚴重扭曲，具體細節可參見文獻[17-18]。

采用基于Abaqus的RITSS方法，將桶基的“貫入-固結-二次加載”過程簡化為二維軸對稱問題進行分析，建立如圖1所示的數值模型。桶基的長度、直徑和厚度分別為L、D和t，支撐電氣平臺的桶基長徑比L/D一般為1左右。為避免邊界效應，土體的徑向與豎向邊界分別取5D與8D。約束底邊徑向邊界的水平與豎向位移，豎向邊界的水平位移。桶基周圍附近土體的網格尺寸為1.5t。為提高計算效率，并避免貫入初期土體表面單元的劇烈變形，將桶基預埋一定深度。后文將證明該設置對后續貫入阻力的影響可以忽略。采用二階軸對稱縮減積分單元劃分土體，以保證映射精度[17]。

圖1 有限元模型范圍與網格

采用罰函數法定義“桶-土”界面的摩擦力，在有效應力分析中，根據庫倫摩擦定律，摩擦力取決于法向有效應力及界面的庫倫摩擦系數。安裝造成桶基周圍土體的應變軟化，設計時一般假定周圍土體完全擾動。Koh等[3]采用大變形方法模擬桶基貫入時取摩擦系數為0.3tanφ′/St，其中φ′為有效內摩擦角，St為靈敏度，得到了與離心機試驗結果相近的貫入阻力。筆者的數值分析也采用相同的摩擦系數。

以往離心機試驗[7-11]比較了吸力式安裝(即通過氣泵從桶蓋與土體之間抽水以施加吸力)與純壓入式安裝(即采用位移控制使桶基貫入)兩種桶基貫入過程，發現安裝方式對黏土中桶基貫入阻力的影響可以忽略不計。因此，通過位移控制模式實現桶基的初始貫入和二次加載。在固結階段，桶基上的豎向力保持為定值，豎向力大小等于當前深度的初始貫入阻力。在整個“貫入-固結-二次加載”過程，僅允許土體表面邊界排水，其余各邊均為不排水邊界。

1.2 土體本構模型

采用經典的修正劍橋模型描述黏性土的有效應力與應變之間的關系，該模型已被成功用于分析靜力觸探儀、樁靴基礎的大變形貫入[18-19]。考慮兩種廣泛用于離心機試驗中的黏性土，即UWA高嶺土和Malaysian高嶺土，相關計算參數總結于表1。側壓力系數K0=1-sinφ′=0.61，浮重度γ′取6.76 kN/m3。

表1 土體參數[9,20-21]

根據臨界狀態土力學理論，軸對稱條件下的不排水強度su為[22]

(1)

基于修正劍橋模型，K0狀態下，正常固結土的孔隙比e與深度z的關系可以表示為

(2)

式中：C為等向固結線與K0固結線的豎向距離,m。根據臨界狀態理論與表1中的參數，UWA高嶺土與Malaysian高嶺土的C分別為0.048與0.057，而滲透系數k與孔隙比e的關系可表示為[18]

(3)

2 數值模擬驗證

目前，還沒有公開發表的桶形基礎二次加載模型試驗或現場監測數據。通過對比桶基首次安裝時的貫入阻力，驗證數值模擬的可靠性。Westgate等[9]完成了6個桶基在UWA高嶺土中貫入的離心機試驗。使用RITSS方法模擬其中名為Type1的試驗。桶基的原型尺寸為D=11.3 m，L=6 m，t=0.05 m。試驗前和試驗后分別進行兩次T-bar試驗，得到土體的平均抗剪強度su=10+2.8zkPa。利用兩次循環插拔T-bar試驗獲得土體的靈敏度St為2.5和3。大變形有限元分析中使用St=3，確定桶-土摩擦系數。試驗中測得γ′=5.9 kN/m3，根據式(1)推導得到修正劍橋模型中的su=1.44zkPa。為便于比較，將貫入阻力F歸一化為F/suA，A為桶形基礎豎向投影面積，A=πD2/4。圖2比較了離心模型試驗[9]、有效應力形式的RITSS法、總應力形式的耦合歐拉-拉格朗日有限元法(Coupled Eulerian-Lagrangian，CEL)結果[13]、DNV規范的預測[23]。圖中的縱坐標為歸一化貫入深度d/D，即貫入深度與桶基直徑之比。

DNV規范與CEL的結果均基于總應力分析，假定安裝時地基土完全不排水。其中DNV規范建議周圍土體的強度折減系數α=1/St=0.33，表示桶的最大側壁摩阻力為αsu。文獻[13]中CEL分析使用的系數為α=0.2，故圖2也給出了α=0.2時的DNV方法預測結果。由圖2可見，有效應力形式的RITSS方法和總應力形式的CEL方法給出的結果都與試驗數據接近，且均在DNV公式預測的α=0.2和α=0.33兩條曲線范圍內。采用RITSS方法時，預埋使得計算曲線從0.2D開始，因此，無法比較0～0.2D范圍內的結果。然而，當桶基的貫入深度超過0.25D后，貫入阻力曲線與試驗非常接近，表明一旦土體破壞模式被完全激發，預埋對貫入阻力不會產生明顯影響，這與已有結論一致[18-19]。圖2的比較也表明大變形有限元分析中的摩擦系數取值合理，可以在很大程度上考慮安裝擾動對貫入阻力的影響。

圖2 數值模擬、離心機試驗、規范預測結果的比較

3 結果分析與討論

實際工程中鋼制桶基的壁厚一般小于0.05 m。試算發現，當桶基側壁厚度較小時，有效應力形式的大變形分析不易收斂，后續分析中統一取壁厚為0.19 m。Xiao等[13]比較了側壁較厚(t=0.2～0.6 m)與較薄(t=0.05 m)的桶基端部承載力系數，發現其范圍在10～12之間，表明一定范圍內壁厚的增加沒有改變桶基貫入時端部以下土體的破壞模式，而側壁附近的土體破壞模式不依賴壁厚。因此，增加壁厚不會顯著影響大變形結果。

表2總結了變參數分析中開展的算例。為探究固結系數cv的影響，算例A～算例D分別模擬了桶基在UWA高嶺土和Malaysian高嶺土中的“貫入-固結-二次加載”過程。這兩種土的本構模型參數接近，但滲透系數相差很大，相同應力水平下Malaysian高嶺土的cv約為UWA高嶺土的10倍。算例A～算例D還將用于考察首次安裝深度zp對二次加載的影響。算例E和算例F將考察桶基直徑D對二次加載的影響。表2中的固結時間tc在7～360 d之間變化，其下限是實際工程中最短的施工間隔，而上限人為取得較大，以便展示充分固結造成的二次加載效應。

表2 變參數分析考慮的工況

修正劍橋模型并不能考慮土體結構性的影響，分析中摩擦系數始終保持為0.3tanφ′/St，這也意味著固結前后的土體靈敏度未改變，即不考慮擾動的土體結構在固結過程中部分恢復，所以，計算結果偏于安全。以下分析中統一取典型值St=2.5。

3.1 固結時間對二次加載的影響

算例A、C(D=9.5 m，首次貫入深度zp=0.9D)，固結7、30、60、120、270 d后，二次加載貫入阻力隨深度的變化如圖3所示。為了直觀顯示固結時間的影響，圖中包含了初始連續貫入阻力曲線作為參考。參考Wang等[19]對樁靴基礎貫入阻力的處理方法，利用峰值貫入阻力比Ncp/Ncr表征貫入阻力的改變，如圖3(a)所示。Ncp為二次加載時歸一化的峰值貫入阻力，Ncr為Ncp出現深度處歸一化的初始貫入阻力。

圖3 固結前和固結后的歸一化貫入阻力(zp=0.9D)

從圖3可以看出：

1)桶基的二次加載阻力在很短位移內即可達到峰值，Ncp出現時的深度僅略大于初始貫入結束時的深度zp。

2)UWA高嶺土中即使固結時間僅有7 d，峰值貫入阻力也會增加約20%，見圖3(a)。如果固結時間為7～120 d，二次加載的Ncp隨固結時間增大，說明固結導致的桶基豎向承載力顯著提高。但120 d后的Ncp基本不再增加。

3)二次加載的貫入阻力曲線達到峰值后迅速降低，并向初始連續貫入阻力曲線靠近。這表明，盡管施工間隔期間桶基周圍土體固結，強度增加，但強度增加土體的范圍有限。這一特征與樁靴基礎的二次貫入類似[19]。

4)Malaysian高嶺土的結果表現出與UWA高嶺土類似的趨勢，如圖3(b)所示，Ncp/Ncr最大約為1.35。然而，固結60 d后Ncp就不再有顯著增長。產生這種差異的原因是兩種土的固結系數cv不同。Malaysian高嶺土的固結系數遠高于UWA高嶺土，初始貫入造成的累積孔壓在固結階段消散得更快，桶基周圍土體更快達到強度增長極限。

圖4展示了固結過程中桶基附近3個位置的土體超靜孔壓消散過程。3個位置分別為桶壁端部中心、端部上方0.02zp處與桶壁緊鄰的內側與外側。為避免固結初期孔壓重分布對消散曲線的影響，利用Sully等[24]提出的方法確定每一個位置的名義初始孔壓Δumax，得到歸一化孔壓Δu/Δumax隨固結時間的變化曲線。圖4表明，3個位置處的超靜孔壓消散規律幾乎一致。對于UWA高嶺土，孔壓在最初的120 d迅速消散，120 d后的消散速度逐漸減緩，這與峰值貫入阻力的增長規律一致(圖3)。同樣，Malaysian高嶺土中超靜孔壓在前60 d急劇下降，且消散速度比UWA高嶺土更快。

圖4 孔壓消散曲線(zp=0.9D)

圖5(a)為桶基在首次安裝時的土體位移增量場，可見端部的土體主要流向桶內，這與Xiao等[13]的結論一致。固結后貫入至0.91D時(圖5(b))，端部置換的土體也主要向桶內移動，且此時桶內土塞的高度略大于桶外。通過比較發現，固結并未改變端部土體的流動模式。

圖5 算例A固結前后土體流動模式的比較

3.2 峰值貫入阻力比的變化規律

圖6展示了峰值貫入阻力比隨歸一化固結時間T的變化。T按式(4)定義。

圖6 峰值貫入阻力比與歸一化時間的關系

(4)

式中：cv,z為深度zp處的固結系數。T實際上包含了多個因素的耦合影響：tc代表了固結時間；D代表了桶基幾何形狀；cv,z代表了土體的固結性質；初始貫入過程的影響則通過zp體現。

由圖6可以看出，zp/D在0.7～1.5范圍內時，兩種高嶺土參數計算得到的峰值貫入阻力比隨歸一化時間的變化趨勢一致。曲線整體上可以分成兩段：當T小于0.015時，Ncp/Ncr隨T大致呈線性增長關系；當T大于等于0.015時，Ncp/Ncr趨近于一個穩定值。為了建立Ncp/Ncr與T的關系，在圖6中添加了趨勢線。地基土充分固結完成時，算例A的Ncp/Ncr達到1.36，算例F的Ncp/Ncr約為1.26，其他算例結果介于二者之間。也就是說，完全固結時二次加載得到的峰值貫入阻力最少增長約26%。

實際條件下，峰值貫入阻力的增長可能超出圖6的有限元結果，這是因為在固結過程中孔壓消散以及土體結構性的恢復都可以造成后續承載力的上升，且二者的貢獻程度相當[14]。受限于本構模型，本文的有限元分析還不能反映土體結構性的恢復，所得峰值貫入阻力增長程度偏于保守。

4 結論

采用有效應力形式的大變形有限元方法模擬海上風電升壓站平臺的安裝過程，再現了桶基“貫入-固結-二次加載”的全過程。數值模擬中考慮了兩種力學性質接近，但固結系數相差10倍的高嶺土，考察了固結時間、首次安裝深度、桶基直徑等關鍵變量對二次加載時貫入阻力的影響。主要結論如下：

1)當固結時間相同時，高固結系數土中二次加載的貫入阻力增長更迅速，但固結系數對峰值貫入阻力比極限值的影響不大，完全固結時的峰值貫入阻力比在1.26～1.36之間。

2)式(4)定義的歸一化固結時間較好地反映了峰值貫入阻力比的變化趨勢。采用歸一化固結時間表征固結程度時，峰值貫入阻力比隨歸一化固結時間約為線性增加，直至達到極限值。

3)實際施工中，應充分利用二次加載前的施工間隔時間，發揮地基土固結后的強度，從而提高二次加載時的極限承載力，達到優化基礎尺寸，降低工程造價的目的。