孟 翔，高路恒，吳 浩，Atangana Njock Pierre Guy

(1.江蘇工程職業(yè)技術(shù)學(xué)院 建筑工程學(xué)院，江蘇 南通 226001；2.上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海200240)

紡錘形樁靴是巖土工程、港口、海岸及海洋工程等領(lǐng)域中常見的深基礎(chǔ)形式，為了承擔(dān)上部結(jié)構(gòu)自重及工作載荷，保證平臺的安全使用，樁靴底部一般設(shè)計成擴大的倒圓錐形或圓盤形，樁靴底部直徑一般為10～20 m，貫入土體中深度可達數(shù)倍直徑以上[1-3]。對樁靴貫入過程進行數(shù)值模擬對于提高設(shè)計水平、保證施工安全具有重大意義。然而，由于樁靴特殊的幾何形態(tài)，其貫入過程常伴隨極大的土體變形及超孔隙水壓力的變化，采用傳統(tǒng)的網(wǎng)格類數(shù)值方法（如FEM）時會遭遇嚴(yán)重的網(wǎng)格畸變而導(dǎo)致計算失敗。因此，早期Hossain[4]，Tho[5]及Zhou[6]等的研究均假定樁靴預(yù)先埋置于一定深度，而后基于擬定的初始應(yīng)力場進行小變形分析。然而，Hossain 等[4]已經(jīng)證明，這種簡化后的小變形分析結(jié)果與連續(xù)貫入大變形分析以及室內(nèi)離心機試驗結(jié)果存在很大誤差。近年來，國內(nèi)外學(xué)者嘗試采用CEL，ALE，RITSS 等大變形有限元分析方法（LDFE）對樁靴的連續(xù)貫入過程進行模擬[4,7-9]。從工程應(yīng)用角度而言，LDFE 的使用極大地促進了人們對樁靴貫入過程中物理現(xiàn)象的理解和問題的解決，但以往LDFE 研究均未能考慮到真實貫入過程中土、水兩相耦合作用，并且頻繁地進行網(wǎng)格重分和物理量映射也會造成計算精度損失[10]。

光滑粒子流體動力學(xué)（SPH），作為一種完全的拉格朗日無網(wǎng)格方法，由于其無網(wǎng)格屬性、粒子屬性及自適應(yīng)性，越來越多地被用于求解各類大變形問題。Maeda 等[11]、Bui 等[12]最早建立土、水耦合SPH 方程并用于模擬射流破土問題。其后，Huang 等[13]在此基礎(chǔ)上將土、水耦合SPH 方程推廣到非飽和土即氣-液-固三相耦合作用模型，并對滲流破壞問題進行研究。Wu 等[14]在SPH 框架下首次建立了土-水-結(jié)構(gòu)耦合SPH 算法并應(yīng)用于貫入問題的模擬。本文采用Wu 等建立的土-水-結(jié)構(gòu)耦合SPH 算法對紡錘形樁靴連續(xù)貫入過程進行數(shù)值模擬，在考慮土體大變形以及土-水-樁靴三者耦合作用條件下研究樁靴貫入阻力及樁靴底部超孔隙水壓力的變化規(guī)律。

1 基本理論

1.1 SPH 的基本原理

SPH 本質(zhì)上是一種插值方法，其主要思想是借助于有限數(shù)量的攜帶場變量和材料屬性的粒子，離散化偏微分方程組（PDEs）所定義的問題域，粒子與粒子之間無固定聯(lián)系。場函數(shù)及其梯度首先通過光滑函數(shù)轉(zhuǎn)化為積分近似表達式，而后通過對支持域內(nèi)相鄰粒子插值轉(zhuǎn)化為粒子近似表達式，從而將PDEs 轉(zhuǎn)化為一系列常微分方程（ODEs）進行求解[10]。場函數(shù)及其梯度的最終粒子近似式[10]為：

式中：f 為三維空間坐標(biāo)向量x 的函數(shù)；i，j 為粒子的標(biāo)識；mj為粒子質(zhì)量；ρj為粒子密度；N 為相鄰粒子總數(shù)；W 為光滑函數(shù)；h 為光滑長度。本文采用Wendland 型[15]光滑函數(shù)，形式如下：

式中：αd為正則化參數(shù)，對于二維問題取值7/（4πh2）；q 是點x 和x'之間的相對距離，q=|x-x'|/h。

1.2 土-水混合物控制方程及SPH 離散化

土-水混合物的控制方程由土、水兩相各自的質(zhì)量守恒方程（4），（5）及動量守恒方程（6），（7）構(gòu)成[14]：

式中：k 為土體的滲透系數(shù)。

式中：下標(biāo)i，j 和a，b 分別用來標(biāo)識土體和孔隙水；υji=υj-υi，即粒子j，i 的速度差；υab=υa-υb；土-水兩相之間的黏滯性拖拽力fia的SPH 粒子近似式為[14]：

方程組（9）～（12）含有未知量：孔隙率n，土體速度νs，水的速度νf，水的密度ρf，土體有效應(yīng)力σ'，孔隙水壓力及黏性力。σ'由土體的本構(gòu)方程求得，本文假定土體為理想彈塑性材料，采用Drucker-Prager 屈服準(zhǔn)則和非相關(guān)聯(lián)的流動法則；假定孔隙水為弱可壓縮性牛頓流體，由孔隙水的狀態(tài)方程求得；由動力黏滯系數(shù)和剪應(yīng)變率求得。土、水兩相的粒子與樁靴之間的接觸力采用摩擦滑移算法[14]計算。

2 樁靴貫入過程模擬及結(jié)果分析

采用上述SPH 算法分別對軸對稱條件下砂土和黏土中的樁靴貫入過程進行模擬，研究土體存在極大變形情況下貫入阻力以及孔隙水壓力的變化規(guī)律。如圖1 所示，樁靴為剛體并以恒定速度0.1 m/s向下貫入土中，與混合物之間采用摩擦滑移接觸。模型左側(cè)為對稱軸邊界，采用鏡像虛粒子模擬。模型右側(cè)和底部分別采用兩排、四排虛粒子進行模擬。

2.1 砂土中貫入模擬及貫入阻力分析

模擬中密、密實砂土兩種工況，模型參數(shù)見表1。其中cf為孔隙水中聲音傳播速度；υ 為土體的泊松比；E 為楊氏模量；c 為黏聚力；?，ψ 分別為土體內(nèi)摩擦角和剪脹角；μ 為液相及水的動力黏滯系數(shù)；ID為土體相對密實度；其他參數(shù)物理意義同上。Qiu 等[8]在研究樁靴基礎(chǔ)穿透破壞問題中，采用CEL 對樁靴基礎(chǔ)的連續(xù)貫入過程進行模擬。為便于對比分析，采用的模型與之保持一致，最大貫入深度為3.6 m。

圖1 樁靴貫入模型Fig.1 Schematic diagram of the simulation arrangement

表1 土體、水的物理參數(shù)以及模型幾何尺寸Tab.1 Dimensions and parameters for soil and water phases

圖2 為工況1 中密砂土中的樁靴貫入過程。初始狀態(tài)下，樁靴底部完全浸沒在自由液面以下、土體表面以上。貫入深度d=2.0 m 時，樁靴底部土體向下、向外發(fā)生位移并發(fā)生隆起變形，周圍土體和水的SPH 粒子遵循各自控制方程運動，位置發(fā)生錯動。d=3.6 m 時，樁靴頂部形成深度約為3.0 m 的柱狀孔，未發(fā)生土體回淤現(xiàn)象?？梢钥闯?，計算過程中水面始終未發(fā)生變化，與直觀認(rèn)識一致。采用SPH 可以很好地模擬土、水混合物大變形過程而不會發(fā)生網(wǎng)格畸變等問題。

圖2 樁靴在砂土（中密）中貫入過程Fig.2 Penetration process of spudcan into medium dense sand

圖3 （a）和（b）分別為樁靴在中密、密實砂土層貫入過程中受到的貫入阻力隨深度變化規(guī)律。從圖3（a）可以看出，SPH 值與CEL 值呈現(xiàn)較為一致的規(guī)律，在正則化深度d/D 達到0.3 時，貫入阻力達到峰值550 kPa，兩者間無明顯差異；隨著貫入深度繼續(xù)增加，CEL 預(yù)測值稍小于SPH 值。這是因為樁靴貫入過程中底部將產(chǎn)生相當(dāng)可觀的超孔隙水壓力并對貫入阻力產(chǎn)生較大影響，忽略孔壓影響可導(dǎo)致最大誤差接近10%；圖3（b）為樁靴在密實砂土中貫入過程受到的阻力。樁靴阻力隨正則化貫入深度逐漸增大至d/D=0.2，而后逐漸減小并穩(wěn)定在550～600 kPa。兩種方法預(yù)測的貫入阻力峰值約為680～700 kPa，均在d/D=0.2 時出現(xiàn)。然而，初始貫入階段CEL 的結(jié)果呈現(xiàn)較為劇烈的振蕩，原因是未能考慮密實砂土發(fā)生剪脹時孔隙率變化的影響。本文采用的算法對孔隙率建立單獨的控制方程并求解，故可以考慮到這一點。

圖3 樁靴貫入阻力隨正則化貫入深度變化Fig.3 Penetration resistance versus normalized penetration depth

2.2 黏土中貫入模擬及孔隙水壓力分析

由于黏土滲透系數(shù)較小，樁靴貫入過程中周圍土體內(nèi)部會產(chǎn)生較大的超孔隙水壓力，而后隨著時間逐漸消散，對樁靴工作階段的安全性有很大影響。早期Hossain 等[4]、Teh 等[9]采用LDFE 對樁靴連續(xù)貫入過程進行了模擬，但由于采用的是總應(yīng)力分析方法，無法有效計算孔隙水壓力。Purwana 等[16]采用離心機試驗方法對樁靴在正常固結(jié)飽和黏土中的連續(xù)貫入過程進行了模擬。圖4 中，T1～T4 為總應(yīng)力測點，P1～P8 為孔隙水壓力測點。試驗在100g加速度環(huán)境下進行，根據(jù)相似準(zhǔn)則等效于直徑12.5 m的樁靴連續(xù)貫入土中深度達19 m。SPH 模型與Purwana 等[16]的離心機模型試驗對應(yīng)的原型尺寸保持一致且材料參數(shù)完全一致，水的相關(guān)參數(shù)為μ=0.001 Pa·s,000 kg/m3, cf=142.0 m/s，黏土相關(guān)參數(shù) 為600 kg/m3, kE=520 kPa/m, E=6.5 MPa,n=0.63, k=2.0×10?7m/s, υ=0.33, c=1 kPa, ?=23°, ψ=0°，幾何尺寸為HW=5.0 m, HSW=30.0 m, WSW=36.0 m,D=12.5 m, Lb=2.0 m, Lm=1.0 m, Lt=1.5 m，其中，kE為楊氏模量沿深度方向變化梯度，其他參數(shù)意義同上。為了減少計算量，僅模擬樁靴的前10 m 貫入過程。

圖4 樁靴貫入離心機模型（單位：mm）Fig.4 Centrifuge model set-up for spudcan penetration (unit:mm)

圖5 和6 分別為樁靴貫入不同深度時土體內(nèi)產(chǎn)生的孔隙水壓力及超孔隙水壓力云圖。可以看出，d/D=0 時，土體內(nèi)部孔隙水壓力呈現(xiàn)出靜水壓力分布規(guī)律，超孔隙水壓力為零。d/D=0.4 時，樁靴周圍產(chǎn)生最大值約為60 kPa 的超孔隙水壓力，并且沿深度方向逐漸衰減至18 m 處，沿水平方向逐漸衰減至13 m 處。地表產(chǎn)生約1.0 m 的隆起變形，并隨著貫入過程繼續(xù)形成柱狀的孔洞直至達到某一臨界深度之后發(fā)生回淤現(xiàn)象。d/D=0.8 時，回淤的土體完全覆蓋樁靴底部，這與Purwana 等[16]室內(nèi)離心機試驗觀測到的現(xiàn)象一致。樁靴周圍產(chǎn)生超孔隙水壓力最大約為130 kPa，隨深度逐漸衰減至26.0 m 處。

圖5 樁靴貫入過程中孔隙水壓力分布云圖（單位：kPa）Fig.5 Contours of pore water pressure during spudcan penetration (unit: kPa)

圖6 樁靴貫入過程中超孔隙水壓力分布云圖（單位：kPa）Fig.6 Contours of excess pore water pressure during spudcan penetration (unit: kPa)

圖7 為樁靴貫入過程中頂部（P1，P2）及底部（P3，P4）4 個測點的孔隙水壓力，縱坐標(biāo)為正則化貫入深度d/D，橫坐標(biāo)為正則化孔隙水壓力，其中γw為水的重度，D 為樁靴的直徑。為便于對比分析，圖7 中給出了Purwana 等[16]的離心機試驗測值及Yi 等[17]的不排水模擬結(jié)果。從圖7（a）和（b）可以看出，d/D<0.4 時，由于土體尚未發(fā)生回淤，樁靴頂部無超孔隙水壓力產(chǎn)生，SPH 與Yi 等的結(jié)果無明顯差異且均與靜水壓力曲線重合。離心機測值稍大于理論值，可能是由于測量誤差。d/D>0.4 時，土體逐漸回淤覆蓋樁靴并受到其拖拽作用向下位移，產(chǎn)生較小的負超孔隙水壓力，故圖7（a）中SPH 結(jié)果及Yi 等的結(jié)果均稍小于靜水壓力。相比較之下，樁靴底部在貫入過程中孔隙水壓力遠大于靜水壓力且大致隨深度d/D 線性增加，如圖7（c）和（d）所示。d/D=0.8 時，離心機測值及SPH 計算得到的均為2.25，而Yi 等的結(jié)果約為2.7。原因是Yi 等采用了完全不排水的假定，故計算結(jié)果偏大。

圖7 貫入過程中樁靴頂部及底部正則化孔壓Fig.7 Normalized pore water pressure at top and bottom of spudcan during penetration

圖8 為樁靴貫入過程中頂部（T1，T2）及底部（T3，T4）的總應(yīng)力，縱坐標(biāo)為正則化貫入深度d/D，橫坐標(biāo)為正則化豎向總應(yīng)力σzz/(γD)，其中γ 為混合物的重度?？倯?yīng)力σzz根據(jù)回淤土體的厚度及混合物的重度計算。由于土體回淤時刻無顯著差異，均在d/D=0.13～0.15 時發(fā)生，故總應(yīng)力數(shù)值相差不大。樁靴頂部σzz/(γD)隨貫入深度可近似看作線性變化，d/D=0.8 時約為0.6；樁靴底部σzz/(γD)在貫入過程初期即d/D<0.3 時呈非線性變化，而后隨d/D 線性增大，d/D=0.8 時約為1.4～1.5 。樁靴頂部T1，T2 處，SPH、離心機測值及Yi 等的結(jié)果無顯著差異，而樁靴底部T3，T4 處差異較大。d/D=0.8 時，SPH 與Purwana 等結(jié)果相差近0.25，原因可能是本文土體采用了Drucker-Prager 屈服準(zhǔn)則，不能反映出土體的體積屈服現(xiàn)象。然而需要說明的是，本文的目的是驗證土-水-結(jié)構(gòu)耦合SPH 算法對大變形貫入問題的適用性，而非提出或改進土體的本構(gòu)模型，為了簡化起見采用了經(jīng)典的理想彈塑性土體本構(gòu)模型。相比較之下，Yi 等[17]的計算結(jié)果與離心機試驗測值誤差幾乎達到30%，說明本文所采用的算法相較于CEL 對于大變形問題更加適用。

圖8 貫入過程中樁靴頂部及底部正則化總應(yīng)力Fig.8 Normalized total stress at top and bottom of spudcan during penetration

3 結(jié) 語

為了克服傳統(tǒng)網(wǎng)格類數(shù)值方法遭遇的網(wǎng)格畸變問題，基于土-水-結(jié)構(gòu)耦合SPH 無網(wǎng)格方法對紡錘形樁靴在飽和砂土及黏土中連續(xù)貫入過程進行模擬，研究了貫入阻力、孔隙水壓力、超孔隙水壓力及總應(yīng)力在土體發(fā)生極大變形情況下的變化規(guī)律。主要結(jié)論如下：

（1）該算法能有效捕捉大變形情況下的土-水混合物自由面特征及貫入阻力，不會發(fā)生網(wǎng)格畸變等問題，并且能夠?qū)⑼馏w、水及結(jié)構(gòu)物三者之間的耦合作用影響計入其中。

（2）樁靴頂部孔隙水壓力在貫入過程中近似靜水壓力分布，樁靴底部正則化孔隙水壓力近似線性分布，d/D=0.8 時，約為2.25；樁靴頂部土體中豎向總應(yīng)力在貫入過程中線性增加，d/D=0.8 時，底部正則化豎向總應(yīng)力σzz/(γD)約為1.4～1.5。相同條件下，不考慮土-水耦合作用的CEL 近似解誤差接近30%。

（3）SPH 無網(wǎng)格法可以作為求解巖土工程大變形問題的有力工具。