基于SPH 無(wú)網(wǎng)格法的紡錘形樁靴連續(xù)貫入過(guò)程模擬

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

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

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

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

1 基本理論

1.1 SPH 的基本原理

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

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

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

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

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

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

式中：下標(biāo)i，j 和a，b 分別用來(lái)標(biāo)識(shí)土體和孔隙水；υ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)的流動(dòng)法則；假定孔隙水為弱可壓縮性牛頓流體，由孔隙水的狀態(tài)方程求得；由動(dòng)力黏滯系數(shù)和剪應(yīng)變率求得。土、水兩相的粒子與樁靴之間的接觸力采用摩擦滑移算法[14]計(jì)算。

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

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

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

模擬中密、密實(shí)砂土兩種工況，模型參數(shù)見(jiàn)表1。其中cf為孔隙水中聲音傳播速度；υ 為土體的泊松比；E 為楊氏模量；c 為黏聚力；?，ψ 分別為土體內(nèi)摩擦角和剪脹角；μ 為液相及水的動(dòng)力黏滯系數(shù)；ID為土體相對(duì)密實(shí)度；其他參數(shù)物理意義同上。Qiu 等[8]在研究樁靴基礎(chǔ)穿透破壞問(wèn)題中，采用CEL 對(duì)樁靴基礎(chǔ)的連續(xù)貫入過(guò)程進(jìn)行模擬。為便于對(duì)比分析，采用的模型與之保持一致，最大貫入深度為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 中密砂土中的樁靴貫入過(guò)程。初始狀態(tài)下，樁靴底部完全浸沒(méi)在自由液面以下、土體表面以上。貫入深度d=2.0 m 時(shí)，樁靴底部土體向下、向外發(fā)生位移并發(fā)生隆起變形，周圍土體和水的SPH 粒子遵循各自控制方程運(yùn)動(dòng)，位置發(fā)生錯(cuò)動(dòng)。d=3.6 m 時(shí)，樁靴頂部形成深度約為3.0 m 的柱狀孔，未發(fā)生土體回淤現(xiàn)象?？梢钥闯?，計(jì)算過(guò)程中水面始終未發(fā)生變化，與直觀認(rèn)識(shí)一致。采用SPH 可以很好地模擬土、水混合物大變形過(guò)程而不會(huì)發(fā)生網(wǎng)格畸變等問(wèn)題。

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

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

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

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

由于黏土滲透系數(shù)較小，樁靴貫入過(guò)程中周圍土體內(nèi)部會(huì)產(chǎn)生較大的超孔隙水壓力，而后隨著時(shí)間逐漸消散，對(duì)樁靴工作階段的安全性有很大影響。早期Hossain 等[4]、Teh 等[9]采用LDFE 對(duì)樁靴連續(xù)貫入過(guò)程進(jìn)行了模擬，但由于采用的是總應(yīng)力分析方法，無(wú)法有效計(jì)算孔隙水壓力。Purwana 等[16]采用離心機(jī)試驗(yàn)方法對(duì)樁靴在正常固結(jié)飽和黏土中的連續(xù)貫入過(guò)程進(jìn)行了模擬。圖4 中，T1～T4 為總應(yīng)力測(cè)點(diǎn)，P1～P8 為孔隙水壓力測(cè)點(diǎn)。試驗(yàn)在100g加速度環(huán)境下進(jìn)行，根據(jù)相似準(zhǔn)則等效于直徑12.5 m的樁靴連續(xù)貫入土中深度達(dá)19 m。SPH 模型與Purwana 等[16]的離心機(jī)模型試驗(yàn)對(duì)應(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ù)意義同上。為了減少計(jì)算量，僅模擬樁靴的前10 m 貫入過(guò)程。

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

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

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

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

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

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

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

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

3 結(jié) 語(yǔ)

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

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

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

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