基于有效應力法的不排水條件下大直徑管樁水平循環動力響應數值分析

曹光偉，歐強，丁選明，周鵬

(1. 重慶大學 土木工程學院，重慶，400045；2. 重慶大學 山地城鎮建設與新技術教育部重點實驗室，重慶，400045)

大直徑單樁基礎是海上風電基礎的首選[1-2]。隨著風機發電功率的不斷增大，風電樁基礎逐漸大直徑化，其直徑甚至可達10 m，然而，目前針對大直徑海洋管樁的設計方法尚不成熟。我國東部沿海廣泛分布深厚軟黏土[3]，歐洲現有設計經驗無法直接用于我國風電建設，水平動載下黏土大直徑管樁的基礎變形及承載力特性有待進一步研究。此外，與陸地不同，海床地基始終處于飽和狀態，需考慮孔壓對土體力學性能的影響：1) 波浪直接作用于海床，在海床表面形成動孔壓邊界，造成土體強度弱化[4]；2) 波浪作用在樁身，樁周土受到擾動，引起土體孔壓累積；此外，隨樁徑增加，排水路徑變長，這種累積更為明顯[5]。鑒于此，基于孔壓-有效應力分析的風電基礎變形及承載力特性研究就顯得尤為重要[6-7]。

數值模擬作為重要的研究手段在海上風電基礎領域得到了廣泛應用[2-4,8-11]，但往往受限于土體本構模型的選擇。以ABAQUS 為例，其內置的M-C模型、D-P模型和MCC模型均無法有效地反映土體動力循環特性(如剛度退化和變形累積)，不適用于水平循環荷載下的樁-土非線性問題。此外，該軟件的孔壓-有效應力分析功能僅限于Standard 靜力分析，但Standard 隱式模塊對復雜接觸和大變形問題(如樁-土接觸和樁的貫入)收斂性極差。ABAQUS/Explicit子塊非常適用于強烈非線性和大變形問題的研究。目前，ABAQUS/Explicit僅具有總應力分析功能，無法直接進行基于孔壓-有效應力的不排水分析。

在以往的不排水分析中，通常采用M-C 模型或D-P模型的總應力進行分析，并取泊松比υ=0.5來考慮不排水過程體變為0的條件。但該方法無法反映動孔壓累積及土體有效應力變化情況，不滿足有效應力分析的要求。為此，YI 等[12-13]基于MCC 模型開發了適用于ABAQUS/Explicit 的不排水分析VUMAT 子程序。MCC 模型不能有效地反映土體動力循環特性(剛度退化和變形累積)，因此，上述VUMAT 子程序得到土體孔壓是不真實的，也不適用于研究海洋大直徑管樁基礎的變形累積問題。有研究表明，海上風電基礎在其生命周期內需要承受約107次負載循環[14]。風電結構對基礎剛度和變形有著極其嚴格的要求，因此，在孔壓-有效應力分析條件下，土體動本構模型二次開發對近海風電工程具有非常顯著的現實意義。本文基于有效應力原理和黏土有效應力本構模型[15]，開發適用于不排水下的有效應力分析的VUMAT子程序。通過單元試驗和離心機試驗的驗證，發現所開發程序可有效模擬循環動載下土體孔壓、累積變形及剛度退化特性，適用于復雜海洋環境荷載下風電基礎的受力分析。在此基礎上，對水平循環動載下軟黏土大直徑單樁響應進行數值研究，以進一步揭示水平循環動載下大直徑單樁基礎變形及剛度退化規律，并給出相應的經驗公式。

1 Explicit 模塊下不排水條件有效應力分析的數值方法

1.1 黏土本構模型

目前考慮黏土循環特性的本構模型中，以DAFALIAS[16]為代表的邊界面模型因參數意義明確和土動力學行為表征合理等優點而備受關注[17-18]。本文采用的模型為ZHOU 等[15]提出的飽和黏土邊界面本構模型。該模型可正確預測正常固結與超固結飽和黏土在小應變與大應變下的彈塑性響應及循環荷載下剛度退化和變形累積特征。在本文研究中，彈性剪切行為取決于當前的平均應力p'和當前的空隙比e，彈性應變增量可表示為：

邊界面方程是控制土體屈服和塑性模量的關鍵，在p'-q平面內的邊界面方程表示為

式中：F為屈服面函數；n和r為邊界面形狀參數，當分別取為1.6 和2 時與修正劍橋模型一致；p0為先期固結壓力；M為臨界狀態應力比。

邊界面初始大小由p0控制。p0的增量取決于塑性體積應變，并應與修正劍橋模型保持一致。基于非關聯流動法則和一致性條件，塑性應變增量可表示為：

1.2 不排水有效應力分析

雖然ABAQUS/Explicit求解得到的均是土體總應力，但總應力{σ}為有效應力{σ′ }與土體孔壓{u}之和。因此，有效應力{σ′} 與土體孔壓{u} 可分別通過定義土骨架的彈塑性剛度矩陣[D′ep]和水體積模量矩陣[Dw]進行求解，進而求得土體總應力{ }σ，見式(9)～(11)。將上述過程編譯到用戶自定義的材料子程序VUMAT中，便可實現不排水條件下的孔壓-有效應力分析[12]。

式中：{dε}為應變增量；[De]為彈性剛度矩陣；g為塑性勢函數；Kw為水的體積模量，本文取2.18 GPa；n為土體孔隙率；E和0分別表示為各元素為1和0的3階矩陣。

1.3 數值實現

VUMAT子程序的主要任務是基于上一步應力狀態及當前給定的應變增量，直接對材料點的應力進行更新并反饋給主程序，無需求解材料的雅可比矩陣[19]。基于此，在VUMAT子程序內部，對所有與土體骨架相關參數均按有效應力狀態進行更新求得有效應力，再對孔壓進行單獨求解，并根據疊加原理求解總應力。對于土體單元的孔壓與有效應力輸出可通過輸出相關狀態變量來實現。

數值積分算法是求解土體有效應力的關鍵。目前，數值積分算法主要有2大類：顯式積分算法和隱式積分算法。考慮到對高度非線性問題的收斂性及計算效率，本文采用帶有誤差控制的改進歐拉應力積分算法[20]來求解積分點處的應力。該算法引入子步概念，根據主程序輸入參數的初值，將問題劃分成多個子步進行逐步求解。當應力預測值誤差滿足給定的控制誤差時(為盡可能控制求解誤差，本文取控制誤差為10-6)，對子步中的材料點應力和狀態變量進行更新。當所有子步的總步長等于1時，計算土體孔壓，更新材料積分點處的總應力和狀態變量，并反饋至主程序，具體流程見圖1，其中：steptime 為ABAQUS 軟件中分析時間；ΔT為子步步長；T為累積步長；Tn為第n次修正后子步步長；ξ為步長修正系數。本程序采用Fortran語言進行雙精度編程，共有9 個材料常數和19 個狀態變量，可實現顯式模塊下飽和黏土的完全排水及不排水的動力分析。

圖1 VUMAT計算流程Fig. 1 Computing flow of VUMAT

2 數值方法驗證

2.1 單元試驗

2.1.1 Itsukaichi海洋黏土不排水動三軸

為了評估本程序對飽和黏土靜力與動力特性的描述性能，對HYODO等[21]的Itsukaichi高塑性海洋黏土動三軸試驗結果進行數值重現。土樣直徑為50 mm，高度為100 mm，在平均圍壓200 kPa下固結，其土體參數見表1。三軸不排水單調加載試驗為應變控制，加載應變速率為0.001 min-1；單向循環試驗為應力控制，加載頻率為0.02 Hz。圖2 所示為土體單調與循環響應的試驗結果與模擬結果對比，其中：pc為初始固結壓力；qs為初始偏應力；qcyc和qd均為循環偏應力。由圖2 可以看出：偏應力預測曲線與試驗曲線基本吻合，VUMAT子程序可較好地模擬Itsukaichi 黏土在單調荷載下不排水過程中的有效應力。對于循環荷載，偏應力計算曲線可以很好地反映黏土在加卸載過程中的位移累積特性，與試驗曲線吻合較好。

圖2 數值模擬結果與試驗結果對比Fig. 2 Comparisons between simulation and test

2.1.2 飽和高嶺黏土不排水動三軸

為了進一步驗證程序的有效性，分別對LI等[22]得出的高嶺黏土單向與雙向循環三軸試驗結果進行了模擬。黏土材料基本參數見表1。試樣直徑為98 mm，高為110 mm，循環加載頻率為0.1 Hz。單向加載試樣平均圍壓為450 kPa，循環動應力幅值為116 kPa；雙向加載試樣平均圍壓為350 kPa，循環動應力幅值為130 kPa。圖3 所示為孔壓-循環次數關系曲線。由圖3 可以看出：無論是雙向還是單向循環加載，孔壓均隨循環次數的增加逐漸累積并趨于穩定，數值曲線與試驗曲線差異較小，總體基本一致。

圖3 孔壓-循環次數曲線Fig. 3 Pore pressure changes versus number of cycles

2.2 單樁離心機試驗

為驗證本程序對側向荷載下樁-土不排水行為預測的性能，選取文獻[23]中的水平加載單樁離心機試驗作為研究對象。試驗縮尺比為1∶48，原型樁外徑為0.91 m，壁厚為50.8 mm，單樁預埋長度為20.2 m，位于泥線以上的偏心距為4.3 m。試驗所使用土為Alwhite 高嶺土，其材料特性見表2。原型鋼樁密度為7 800 kg/m3，彈性模量為206 GPa，屈服強度為414 MPa。由于模型試驗的對稱性，取樁土模型的一半進行建模。土體采用增強沙漏控制的C3D8R實體單元建模，單樁采用C3D8I單元。為消除邊界的影響，土體半徑取20 倍樁徑。土體底部采用固定邊界，采用外側弧面約束側向位移；此外，限制樁和土對稱平面的Y向位移。對淺層和靠近樁身(5D內，D為單樁直徑)的土體網格進行加密處理，保證樁-土分析的計算精度。采用庫侖摩擦法則和考慮接觸后分離的硬接觸的面-面接觸屬性表征樁-土界面行為。根據文獻[24]推薦公式計算樁-土界面的摩擦因數，本文取0.3。試驗循環加載為單向循環位移控制，荷載頻率為1 Hz。圖4所示為單調與循環荷載下數值與試驗結果對比。從圖4(a)可以看出：數值曲線可以很好地模擬單樁荷載-位移曲線的加工硬化特性；此外，在深度為6D處的歸一化循環P-y曲線的模擬曲線中，隨著循環次數的增長，滯回曲線逐漸下移，即存在剛度與承載力退化現象，這與試驗結果是一致的。由于文獻[23]并未給出高齡黏土的邊界面形狀參數n與r，因此，本文采用了與修正劍橋一致的形狀參數(即n=1.6和r=2)，這可能也是引起計算與試驗誤差的主要原因。雖然模擬荷載-位移曲線略低估了單樁承載性能且單圈滯回圈面積比試驗結果略大，但數值與試驗結果的整體趨勢基本一致。

表2 Alwhite高嶺土參數[23]Table 2 Parameters for Alwhite kaolin[23]

圖4 離心機與數值結果對比Fig. 4 Comparisons between centrifuge and simulation results

3 計算結果及分析

在涉及樁-土相互作用的風機動力分析中，常采用子結構法，因此，水平循環荷載下大直徑管樁在泥線處的剛度和位移累積規律研究是十分重要的。本節針對不同循環幅值和樁徑的大直徑管樁，討論樁的位移累積、剛度退化及孔壓累積規律。樁身錨固長度Lp為30 m，壁厚取樁徑的1.1%。樁土模型的建模方法與2.2節的相同，土體與管樁的基本材料參數也與文獻[23]中的一致。樁頂荷載采用力控制的單向循環水平荷載，加載頻率為0.5 Hz，最大循環次數為100 次。樁土相對剛度系數Kr的計算及類型采用Poulos 準則[25](即剛性和柔性樁的Kr上限和下限分別為0.208 0和0.002 5)進行判別。樁研究參數見表3。循環幅值比定義ξ為加載幅值Fcyc比單樁極限荷載Fu，其中極限荷載Fu取泥線0.2D側向位移所對應荷載[2]。

表3 樁研究參數Table 3 Study parameters for piles

3.1 荷載-位移曲線

圖5 所示為循環荷載下樁泥面處荷載-位移曲線(考慮到原荷載-位移過于密集，對其進行了過濾處理)。由圖5 可以看出：在整個加載過程中，樁泥線處的側向位移出現不可恢復的塑性變形，有明顯的棘輪效應。循環幅值比對單樁的荷載-位移響應影響非常顯著，即循環荷載比越大，由加載引起的土體擾動越大，樁泥線處的位移累積越大。類似地，對于同一循環荷載比下不同直徑的樁，樁徑越大，累積位移越大。

圖5 泥線處的荷載-位移曲線Fig. 5 Load-displacement curves at mudline

3.2 累積位移

圖6所示為循環荷載下樁泥面處位移隨循環次數的演化過程。由圖6可以看出：樁泥線處的側向位移ym隨循環次數N的增加而增加；樁的累積位移與循環次數基本呈線性趨勢，且在100個循環內無衰減跡象，這與文獻[26]中軟黏土的離心機試驗結果較為類似；另外，隨循環幅值的增加，樁泥線處的位移累積速率顯著增加；單樁直徑越大，累積速率也越大，這與3.1節中的結果一致。使用首次循環卸載后的殘余變形(荷載卸至0時的位移)yres,1對殘余側向位移進行歸一化處理，結果見圖7。從圖7 可知：歸一化后的累積位移受幅值比控制，但不受直徑影響。為預測水平荷載下大直徑樁的累積變形情況，采用對數坐標對數據進行擬合。式(12)～(14)為循環次數與累積位移的經驗公式，預測公式所得結果與數值結果吻合較好。

圖6 泥線處的位移響應Fig. 6 Displacement response at mudline

圖7 歸一化的累積側向位移與循環次數曲線Fig. 7 Normalized cumulative lateral displacement versus number of cycles

式中：yres,1和yres,N分別為第1 次和第N次循環的樁泥線處的水平累積位移；Fcyc和Fu分別為循環荷載和靜力極限荷載；N為循環次數。

3.3 卸載剛度

對于風機的設計，基礎剛度是一個非常重要的參數。定義循環荷載-位移曲線的卸載段最高點與最低點(0 荷載點)之間連線的斜率為卸載剛度。考慮到剛度的退化主要是土體塑性變形引起的，因此，相比于循環次數，樁基礎的殘余塑性位移更適合用于表征單樁剛度的退化規律。圖8所示為歸一化后的卸載剛度kun與yres的演化關系曲線。由圖8可以看出：在循環初期，樁的卸載剛度隨累積變形急劇衰減，然而，隨累積變形的持續增長，剛度衰減逐漸趨于平緩，這點與文獻[23]的報道一致。此外，不同幅值下的歸一化曲線基本一致，歸一化的卸載剛度退化受循環幅值比的影響很小；對于不同直徑的大直徑單樁，直徑增加會使卸載剛度衰減更迅速，存在“直徑效應”。造成這一現象的原因是：樁徑越大，加載引起的土體擾動越大，剛度衰減越迅速。式(15)反映了卸載剛度與累積殘余位移之間的經驗關系，可有效地反映出單樁直徑對于剛度退化的影響，與數值結果整體吻合度較高。

圖8 歸一化的累積位移與卸載剛度關系Fig. 8 Normalized relationship of unloading stiffness with cumulative displacement

式中：kun,1和kun,N分別為第1次和第N次循環的樁泥線處的卸載剛度；Dref為參考直徑，本文取4 m。

3.4 土體孔壓

圖9 給出了直徑為4 m 和8 m 樁在100 個荷載循環后的土體超孔壓(狀態變量SDV5)分布。由圖9可以看到：單樁底部出現明顯的超孔壓累積，且樁徑變大后孔壓區域顯著增加，尤其是樁底。造成這一現象的原因是：在水平荷載下，剛性樁(或半剛性樁)樁體發生剛體轉動，樁底出現“踢腳”位移，對樁底土造成擾動使土體超孔壓增大。因此，在工程實際中，應考慮樁底超孔壓對土體的弱化。圖10 和圖11 所示分別為深度3 m 和樁底處(見圖9 紅色方框)的超孔壓與循環次數關系曲線。由圖10 和圖11 可以看出：加載初期孔壓增長迅速，但后期逐漸趨于穩定；另外，荷載增加，樁周土擾動程度加劇，超孔壓幅值及其累積顯著增加，且初期增長速率增加，但土體超孔壓對樁徑不敏感。雖然土體超孔壓累積量不受樁徑控制，但樁徑增加后孔壓影響區域顯著增加(見圖9)，這意味著樁徑增加會加劇整個土體性能下降，進而影響單樁的受力性能。

圖9 100個循環后的土體超孔壓分布Fig. 9 Distributions of excess pore pressure after 100 cycles

圖10 不同加載幅值下樁側超孔壓與循環次數關系曲線Fig. 10 Excess pore pressure changes with number of cycles considering different cyclic amplitudes

圖11 樁底超孔壓與循環次數曲線Fig. 11 Excess pore pressure changes with number of cycles at pile base

4 結論

1) 不排水單元試驗及水平單樁離心機試驗結果表明所開發VUMAT子程序可有效模擬循環荷載下土體在不排水過程中的孔壓和位移累積特性及剛度退化行為且具有較高精度。

2) 在水平循環荷載下，樁側向位移累積隨加載次數的增加而增加，有明顯的棘輪效應，并與循環次數大致呈線性關系。歸一化的位移累積速率會隨著循環幅值的增加而顯著增加，但與直徑無關。

3) 隨累積變形增加，單樁卸載剛度衰減速率先大后小，大致呈冪函數關系；歸一化的卸載剛度和累積變形演化關系與循環幅值無關，但與直徑有關，即卸載剛度的衰減存在“直徑效應”，樁徑越大，衰減速率越大。

4) 對于大直徑海洋單樁，單樁底部會出現顯著的孔壓累積，工程實際應考慮樁底超孔壓對土體弱化的影響。此外，超孔壓累積量受循環幅值控制；超孔壓對樁徑不敏感，但樁徑增加后其影響區域顯著增加。