魚雷錨上拔承載力的物質(zhì)點法數(shù)值分析

程 鈺，邱長林

(天津大學 建筑工程學院，天津 300072)

魚雷錨是一種用于系泊深海設(shè)施的基礎(chǔ)，形似魚雷，其外表為錐形的圓柱形鋼管，內(nèi)用混凝土和廢金屬填充。安裝時將魚雷錨在一定水深處釋放，其依靠自身重力下落獲得的動能上拔海床，安裝快速簡便，是傳統(tǒng)海上錨固件經(jīng)濟性的替代方案。2001年Medeiros C J[1]概述了魚雷錨在巴西近?？财账篂?Campos Basin)進行魚雷錨現(xiàn)場原位試驗，水深范圍200～1 000 m，并分析魚雷錨的貫入深度與其上拔承載力，為魚雷錨的工程施工規(guī)范和理論分析發(fā)展奠定了基礎(chǔ)。國內(nèi)目前一般海洋工程項目基礎(chǔ)的安全承載力設(shè)計主要參考API規(guī)范[2]。但由于理論估算受參數(shù)影響大，結(jié)果存在不穩(wěn)定性，近年來越來越多的學者采用模型試驗和數(shù)值分析方法研究魚雷錨上拔承載力。

Wenkai Wang等[3]進行了240組實驗室試驗，將11種不同形狀的魚雷錨從不同類型、不同埋深的粘性土體中垂直拉出，提出魚雷錨不排水單調(diào)抗拔承載力估算公式。Hossain M S[4]等，通過200 g的離心機模型試驗探究黏土和淤泥中魚雷錨的單向抗拔承載力，試驗結(jié)果表明，承載力與安裝后固結(jié)時間、錨埋深、土體不排水的抗剪強度呈正相關(guān)。相比于現(xiàn)場試驗或模型試驗，數(shù)值分析方法可更細致直觀研究魚雷錨的承載能力。目前魚雷錨承載力數(shù)值分析方法主要為有限元法，包括常規(guī)小變形有限元方法和可以考慮大變形的幾何非線性有限元方法，如CEL、ALE等。Chen[5]通過ABAQUS有限元方法對動力錨在15°～80°不同拉拔傾角下承載力，并分析其影響因素。瑜璐等[6]應(yīng)用ABAQUS軟件分析特定埋深下魚雷錨的抗拔承載力，探討了錨型、土體類型、拉拔荷載傾角、拉拔荷載水平分量與錨翼夾角等多種因素對拉拔承載力的影響并得出規(guī)律性結(jié)論，最終提出預(yù)測魚雷錨承載力的歸一化V-H包絡(luò)公式。Raie[7]通過CEL方法分析研究魚雷錨安裝成功后魚雷錨周圍地基土體二次固結(jié)情況，結(jié)果直觀反映了超孔隙水壓力的消散速率以及土體強度恢復(fù)情況。關(guān)于接觸問題，CEL方法無法直接模擬正常固結(jié)土中錨-土之間的摩擦作用，大多數(shù)學者采用True提出的牛頓第二定律來模擬錨-土之間的相互作用[8]，但該方法無法真正反映土體的率效應(yīng)和應(yīng)變軟化效應(yīng)對摩擦力的影響。

魚雷錨安裝后的承載力對海上工程設(shè)計施工至關(guān)重要，目前研究表明，魚雷錨的承載力還缺乏有效的數(shù)值模擬方法，這導致魚雷錨的研究還很不成熟。物質(zhì)點易于追蹤物質(zhì)運動及邊界情況，規(guī)避了網(wǎng)格大變形問題，在解決魚雷錨上拔過程土體大變形問題上具有十分明顯的優(yōu)點。針對上述問題，在已有研究基礎(chǔ)上，通過三維顯式物質(zhì)點法開源程序[9]模擬魚雷錨在土體中上拔過程，采用接觸界面法向量與API相結(jié)合的接觸算法模擬錨-土之間摩擦力，并且通過對現(xiàn)場原位試驗結(jié)果進行可靠性分析，與現(xiàn)場試驗結(jié)果相互驗證補充，為魚雷錨上拔承載力研究提供一種新的數(shù)值算法，此外，還對上拔過程周圍土體變形規(guī)律和上拔承載力影響因素進行分析，為工程實際提供指導性建議。

1 物質(zhì)點法方法基本理論

1.1 物質(zhì)點離散

圖1 物質(zhì)點法示意圖[10]Fig.1 Schematic diagram of the material point method[10]

在物質(zhì)點法(MPM)中，連續(xù)體和其運動空間分別用物質(zhì)點和網(wǎng)格兩套不同的系統(tǒng)表示，如圖1所示。MPM先將連續(xù)體離散為一組質(zhì)點，物體材料區(qū)域內(nèi)的物質(zhì)信息均由質(zhì)點攜帶，并由質(zhì)點代表其所在區(qū)域，在每個計算時間步內(nèi)，質(zhì)點和計算背景網(wǎng)格不發(fā)生相對運動，質(zhì)點會根據(jù)具體受力情況在背景網(wǎng)格中運動，背景網(wǎng)格只用于求解動量方程和計算空間導數(shù)，計算網(wǎng)格結(jié)點與質(zhì)點間通過單值映射函數(shù)進行運動量信息傳遞。

MPM將連續(xù)體離散后，其密度為

(1)

式中：mp是質(zhì)點p的質(zhì)量；np為離散質(zhì)點總數(shù)；δ(x)是Dirac Delta函數(shù)；xip是質(zhì)點p的坐標。

1.2 動量方程及其求解

在不考慮熱量交換情況下，連續(xù)體更新拉格朗日的動量方程和邊界條件分別為

(2)

(3)

為求解動量方程，取虛位移δuj∈R0和權(quán)函數(shù)R0={δuj|δuj∈C0，δuj|Γu=0}，結(jié)合給定面力邊界條件(3)式，得出動量方程(2)式等效積分弱形式為

(4)

將(1)式代入(4)式中，得到連續(xù)體離散后動量方程等效積分弱形式

(5)

由于應(yīng)用物質(zhì)點法求解動量方程時，質(zhì)點和背景網(wǎng)格結(jié)點之間信息傳遞是通過在背景網(wǎng)格結(jié)點上建立的有限元函數(shù)NI(Xi)來實現(xiàn)的。因此，質(zhì)點p的位移uip和虛位移δuip為

uip=NIpuiI

(6)

δuip=NIpδuiI

(7)

式中：uiI為結(jié)點位移；下標I、J、K表示網(wǎng)格結(jié)點。將(6)、(7)式代入(5)式中，得到背景網(wǎng)格結(jié)點的運動方程

(8)

(9)

(10)

網(wǎng)格結(jié)點I在i方向的動量piI為

(11)

式中：網(wǎng)格質(zhì)量矩陣mIJ為

(12)

應(yīng)用集中質(zhì)量陣(ng為結(jié)點總數(shù))

(13)

可以將結(jié)點運動方程(8)式簡化為

(14)

之后，利用中心差分法對(14)式進行求解，得到當前時間步背景網(wǎng)格結(jié)點的位移增量，再借助形函數(shù)插值映射到質(zhì)點，得到質(zhì)點相關(guān)的其他物理量。運動方程求解完成，計算步結(jié)束，拋棄變形后的網(wǎng)格。

1.3 接觸算法

解決無滑移黏著接觸問題，標準MPM不需采取額外的處理，當需要考慮對象間的相對滑動和分離時，需要將接觸算法引入到MPM中。采用接觸界面法向量計算方法[10]與API方法結(jié)合模擬錨-土間摩擦作用力。

通常來說，當兩個物體r和s間滿足接觸判定條件式(15)時

(15)

考慮到兩物體法向單位向量不完全共線時，會導致動量不守恒和界面穿透情況，故兩物體外法向單位向量選取至關(guān)重要。當物體r比物體s更硬時，應(yīng)選取物體r來計算接觸面的公法線方向，即取

(16)

(17)

(18)

(19)

(20)

1.4 材料模型

物質(zhì)點法中，本構(gòu)方程在質(zhì)點上計算。由于鋼的強度和模量遠大于土體，魚雷錨的本構(gòu)模型采用線彈性模型，其應(yīng)力應(yīng)變關(guān)系為

(21)

(22)

土體本構(gòu)模型采用理想彈塑性模型，其破壞形式包括剪切破壞和拉伸破壞。土體剪切失效模型采用Drucker-Prager模型，其屈服函數(shù)為

(23)

式中：J2為偏應(yīng)力張量的第二不變量；kΦ、qΦ與材料的黏聚力c和摩擦角φ的關(guān)系由下式確定

(24)

(25)

式中：加號表示D-P屈服面在π平面上內(nèi)接Mohr-Coulomb屈服面，減號則表示外接，計算采用內(nèi)接圓來進行計算[11]。土體拉伸失效屈服函數(shù)為

ft=σm-σt

(26)

式中：σt為材料抗拉強度。剪切失效和拉伸失效相應(yīng)的流動法則均采用關(guān)聯(lián)流動法則。

為了保證應(yīng)力率不受剛體轉(zhuǎn)動的影響，每個物質(zhì)點的應(yīng)力率都采用焦曼應(yīng)力率進行應(yīng)力控制。剪切失效和拉伸失效相應(yīng)的流動法則均采用關(guān)聯(lián)流動法則。

圖2 魚雷錨上拔示意圖Fig.2 Schematic diagram of torpedo anchor pulling up

2 魚雷錨上拔過程物質(zhì)點法分析

以巴西近?？财账篂臭~雷錨上拔現(xiàn)場原位試驗[1]為研究對象，通過物質(zhì)點法建立三維模型對其進行模擬，從而研究無翼魚雷錨在飽和黏土中上拔承載力。

2.1 分析模型

巴西近海坎普斯灣現(xiàn)場試驗[1]將在正常固結(jié)黏土中的泥線角θ0= 0°的無翼魚雷錨從20 m深處拉出，但并未介紹錨鏈拉出角θa，拉出角θa定義為負載方向和水平方向之間的角度，考慮到現(xiàn)場試驗實際情況，與豎向拉出情況類似，故取尾翼拉出角θa=90°(θ0，θa如圖2所示)。

圖3 魚雷錨幾何尺寸Fig.3 Torpedo anchor geometry

巴西近?？财账篂超F(xiàn)場試驗[1]的T-40魚雷錨，直徑D=0.76 m，長L=12 m，錨干重為240 kN，無尾翼，魚雷錨形態(tài)及其參數(shù)如圖3所示，下端部為圓錐形狀，在模型中其離散間距為0.3 m的物質(zhì)點。考慮到魚雷錨上拔影響范圍，土體模型尺寸取10 m×10 m×60 m，離散間距為0.5 m的物質(zhì)點。整個模型總共55 390個物質(zhì)點。模型計算區(qū)域為10 m×10 m×60 m，網(wǎng)格單元類型為八節(jié)點六面體，網(wǎng)格間距為0.5 m。魚雷錨和土體物質(zhì)點離散模型如圖4所示。

圖4 魚雷錨-土物質(zhì)點離散模型Fig.4 Discrete model of torpedo anchor and soils

魚雷錨材料為鋼，其本構(gòu)參數(shù)見表1。結(jié)合已有的研究結(jié)果[12]，土體為軟黏土，其內(nèi)摩擦角為0°，不排水抗剪強度Su取(5+2z)kPa，其中z為地基深度，彈性模量E0取500Su。土體本構(gòu)參數(shù)見表2。

表1 魚雷錨本構(gòu)參數(shù)Tab.1 Torpedo anchor constitutive parameters

表2 土體本構(gòu)參數(shù)Tab.2 Soil constitutive parameter

模型初始狀態(tài)為魚雷錨埋入土體中，魚雷錨的初始埋深H為20 m[12]。邊界條件施加在計算網(wǎng)格上，4個側(cè)面和底面均為固定。

外部荷載方面，魚雷錨上拔模擬是通過施加豎直向上恒定錨鏈拉力將魚雷錨拔出，阻力僅考慮側(cè)摩阻力。接觸選取魚雷錨作為主體、主平面，取魚雷錨的表面法線為公法線，摩擦系數(shù)μ取0，取API公式計算由土體粘聚力產(chǎn)生的側(cè)摩阻力作為切向強度。

API規(guī)范[2]中規(guī)定對于黏土中嵌入式樁，沿樁長向側(cè)摩阻力Fs為

Fs=f(z)As=αsuAs

(27)

式中：不排水抗剪強度Su取錨-土接觸點埋深處對應(yīng)值，界面摩擦比α取靈敏度的倒數(shù)1/St，仍采用已有的研究結(jié)果St取3，即α取1/3[12]。

算法采用顯示積分格式，計算時間步長0.167E-04，計算時間為2.0 s，共60 010步。

圖5 魚雷錨豎向上拔阻力-歸一化位移關(guān)系曲線Fig.5 Vertical pull-up resistance-normalized displacement relationship curve of torpedo anchor

2.2 魚雷錨上拔承載力驗證

在巴西近?？财账篂尺M行的無翼魚雷錨現(xiàn)場試驗中，平均貫入深度為20 m。如圖3所示魚雷錨的設(shè)計承載力為1 400 kN[1]，現(xiàn)場原位試驗實測極限承載力在900～1 100 kN范圍內(nèi)變化，平均極限承載力為1 000 kN。根據(jù)現(xiàn)場原位試驗數(shù)據(jù)，對埋深H=20 m的無翼魚雷錨進行物質(zhì)點法上拔數(shù)值模擬，得到魚雷錨豎向上拔阻力與歸一化位移關(guān)系曲線如圖5所示。通過計算發(fā)現(xiàn)，錨歸一化位移(錨尖位移與錨尖初始埋深比值)U/H為0.08～0.16范圍內(nèi)豎向上拔阻力達到最大值，故取位移達0.08H～0.16H所對應(yīng)的承載力平均值為豎向抗拔承載力，即最終豎向上拔承載力為1 187 kN，處于現(xiàn)場試驗測得的抗拔承載力范圍內(nèi)，與現(xiàn)場平均極限承載力相差18.7%，物質(zhì)點法數(shù)模結(jié)果與現(xiàn)場測試數(shù)據(jù)的合理一致性證實了該物質(zhì)點法數(shù)模在評估魚雷錨在黏土中的抗拔承載力的可靠性。

2.3 上拔過程地基變形特性

為便于觀察研究結(jié)果，其中計算輸出區(qū)域為錨兩側(cè)各6D，錨在土中運動時深度向取26D，錨拔出土體時深度向取32D。

物質(zhì)點法可以通過土體質(zhì)點位移變化，得出土體顆粒行動方向，進而得出某一時刻土體具體變形情況。取錨拔出4個典型時刻土體豎向位移進行分析，如圖6所示(圖中位移正值表示位移向上)。由圖6可看出，錨上拔過程地基土體豎向位移以向上為主，錨尾部附近地基土體變形大于錨側(cè)附近地基土體，地基土體產(chǎn)生向上的豎向位移的范圍為以初始錨尖埋深為長軸、12D為短軸的半橢圓形范圍；地基土體產(chǎn)生向下的豎向位移的范圍集中在初始錨尖埋深下方6D，且距錨軸線位置2D～4D范圍，范圍幾乎不變。由于錨附近地基土體受到擠壓作用，從徑向向外擴張、流動，距錨軸線越遠，受影響程度越小，故遠錨側(cè)土體豎向位移也就越小。此外，錨上拔過程中帶動地基土體運動，運動軌跡產(chǎn)生的空腔在錨完全拔出后未被土體完全填充。由圖6-c～圖6-d可看出，錨雖然在不斷上拔，但從錨體部分離開土體開始，地基土體擾動區(qū)域幾乎沒有發(fā)生改變，且豎向位移數(shù)值基本一致。由圖6可看出，錨上拔完成變形土體區(qū)域為以初始錨尖埋深下方6D為高，12D為底邊的等腰三角形范圍。

6-a H=18 m6-b H=12 m6-c H=4 m6-d H=0 m圖6 地基土體豎向位移等值線云圖(單位:cm)Fig.6 Contour cloud diagram of vertical displacement of foundation soil

上拔過程地基土體累計等效塑性應(yīng)變分布如圖7所示。為便于觀察錨周土體塑性區(qū)分布，土體范圍取錨水平向左右6D，深度向取海床表面以下40D范圍。黏土室內(nèi)三軸試驗累計軸向應(yīng)變達到15%即認為其產(chǎn)生破壞，有研究指出沖擊荷載作用下軟黏土室內(nèi)三軸試驗累計軸向應(yīng)變達到2%時土體產(chǎn)生屈服[13]，故認為累計等效塑性應(yīng)變大于10%為土體塑性變形顯著區(qū)域。由圖7可看出，隨著錨不斷上拔，累計等效塑性應(yīng)變區(qū)域由錨身開始逐漸向周圍擴散，且變形范圍不斷增加，部分錨體拔出土體后變形區(qū)域基本不變。此外，隨著魚雷錨不斷上拔，除了錨上拔路徑周圍土體產(chǎn)生變形外，海床表面兩側(cè)土體同樣產(chǎn)生較大范圍變形，魚雷錨回收利用時要考慮土體變形范圍，以防錨軌跡空洞被兩側(cè)土體逐漸填補出現(xiàn)土體塌陷，從而對其周圍設(shè)施產(chǎn)生影響。

對于變形區(qū)域，由圖7可看出，錨側(cè)附近地基土體向兩側(cè)擴散，徑向距離錨側(cè)越近變形越大，直到部分錨體離開土體，錨上拔路徑空腔區(qū)域出現(xiàn)縮小情況(圖7-c)，主要原因是空腔由于深度過高，土體在重力作用下有向腔內(nèi)垮塌，導致腔體側(cè)壁中下部累計等效塑性應(yīng)變增大。由圖7-c、圖7-d中可以看出，魚雷錨拔出土體后土體塑性區(qū)無明顯變化。圖7結(jié)果表明，錨上拔顯著塑性變形區(qū)為以錨尖初始埋深下5～6倍錨徑位置的地基埋深為頂點、以地基表面12倍錨徑為底邊的等腰三角形。

7-a H=18 m7-b H=12 m7-c H=4 m7-d H=0 m圖7 累計等效塑性應(yīng)變云圖Fig.7 Cumulative equivalent plastic strain cloud diagram

3 豎向上拔承載力影響因素分析

魚雷錨貫入地基后的承載力與錨的質(zhì)量和幾何形狀、土體強度分布和埋深有關(guān)，又考慮到土體強度分布和埋深對魚雷錨安裝后承載力起到關(guān)鍵性的作用，通過物質(zhì)點法數(shù)值模擬，得出埋深和土體強度對豎向上拔承載力的影響。

3.1 埋深對豎向上拔承載力的影響

保持魚雷錨及土體參數(shù)與現(xiàn)場原位試驗數(shù)據(jù)[1]相同，現(xiàn)場原位試驗中平均貫入深度為20 m，在典型的深水黏土沉積物中，錨埋深可達到錨長的3倍，豎向單調(diào)上拔承載力通常小于錨干重的5倍，此外有數(shù)值研究[4]表明，錨的埋深為錨長度的1.5～2.6倍，錨豎向上拔承載力在錨干重的2.4～4.1倍，故埋深分別選取20 m、25 m、30 m、35 m，分別對應(yīng)錨長的1.7倍、2.1倍、2.5倍、2.9倍。為深入研究埋深對豎向上拔承載力的影響，模擬魚雷錨在正常固結(jié)土(Su=5+2zkPa)中上拔，得到上拔過程中不同埋深下魚雷錨歸一化上拔阻力隨歸一化位移(錨埋深與參考埋深H0=20 m之比)變化曲線如圖8所示。由圖8可看出，隨埋深增加，魚雷錨上拔阻力穩(wěn)定所需歸一化位移略有增加，說明埋深增加時，土體強度越大，土體對魚雷錨所產(chǎn)生的阻力相應(yīng)較大。

魚雷錨歸一化豎向上拔承載力隨埋深比(初始錨尖埋深與錨長之比)變化曲線如圖9所示。由圖9可看出，在實驗范圍內(nèi)，魚雷錨歸一化豎向上拔承載力與埋深比近似成線性關(guān)系，線性比約為8.34，二者線性化公式為

(28)

在實驗范圍內(nèi)，魚雷錨埋深每增加5 m，相應(yīng)豎向上拔承載力增加錨重的3.1～4.0倍，魚雷錨豎向上拔承載力約為錨重的5.5～15.5倍。

圖8 魚雷錨上拔阻力-位移關(guān)系曲線Fig.8 Torpedo anchor normalized uplift resistance-normalized displacement relationship curve圖9 魚雷錨豎向上拔承載力隨埋深比變化曲線Fig. 9 Curve of normalized vertical pull-out bearing capacity of torpedo anchor with depth-to-depth ratio

圖10 魚雷錨歸一化上拔阻力-歸一化位移關(guān)系曲線Fig.10 Torpedo anchor normalized uplift resistance-normalized displacement relationship curve圖11 魚雷錨歸一化豎向上拔承載力隨k值變化曲線Fig.11 Curve of normalized vertical pull-out bearing capacity of torpedo anchor with k value

3.2 土體強度對豎向上拔承載力的影響

為了研究不同性質(zhì)黏土對豎向上拔承載力的影響，保持上述魚雷錨和土體參數(shù)不變，埋深H=20 m不變情況下，不排水抗剪強度Su取(5+kz)kPa，其中k分別取0、1、2、3，得出魚雷錨在不同不排水抗剪強度Su土體中歸一化上拔阻力與歸一化位移關(guān)系曲線如圖10所示。

由圖10可見，隨k值增加，土體強度增加，土體對魚雷錨所產(chǎn)生的阻力相應(yīng)較高，且隨k值增加，魚雷錨上拔阻力增加階段越來越長，即所需位移越來越大。

如圖11為不同k值下魚雷錨歸一化豎向上拔承載力隨k值變化曲線。由圖11可見，歸一化豎向上拔承載力與k值呈正相關(guān)，二者關(guān)系可擬合為

Q=(1.305k+0.697 9)WA

(29)

4 結(jié)論

應(yīng)用物質(zhì)點法對室內(nèi)模型試驗和巴西近海坎普斯灣[1]原位試驗中埋深20 m的無翼魚雷錨上拔過程進行物質(zhì)點法分析，得到魚雷錨豎直上拔承載力，并進一步研究上拔過程地基土體變形規(guī)律、塑性變形區(qū)范圍與豎向上拔承載力影響因素，得出以下結(jié)論：

(1)物質(zhì)點法數(shù)模結(jié)果與現(xiàn)場試驗結(jié)果之間的比較表明，應(yīng)用物質(zhì)點法研究魚雷錨上拔過程規(guī)律是可靠的、有效的。

(2)錨上拔過程地基土體豎向位移以向上為主，錨尾部附近地基土體變形大于錨側(cè)附近地基土體，豎向位移向上地基土體范圍為以初始錨尖埋深為長軸，12D為短軸的半橢圓形范圍；豎向位移向下地基土體集中在初始錨尖埋深下方6D，12D為底邊的等腰三角形范圍，范圍幾乎不變。

(3)錨上拔顯著塑性變形區(qū)為以錨尖初始埋深下5～6倍錨徑位置的地基埋深為頂點、以地基表面12倍錨徑為底邊的等腰三角形。因此，在魚雷錨上拔回收過程中，需要考慮該范圍內(nèi)其對其他結(jié)構(gòu)物的影響。

(4)魚雷錨豎向上拔承載力影響因素分析結(jié)果表明，試驗范圍內(nèi)，埋深、土體強度與豎向上拔承載力呈正相關(guān)。在實驗范圍內(nèi)，不同埋深下，魚雷錨歸一化豎向上拔承載力與深長比近似呈正相關(guān)；不同土體強度下，魚雷錨歸一化豎向上拔承載力與k值呈正相關(guān)。