基于強度折減理論的海底斜坡穩(wěn)定性分析

楊林青，王忠濤，李家鋼，劉樂軍，胡光海

(1.大連理工大學 建設(shè)工程學部，遼寧 大連116024；2.中海石油研究中心工程設(shè)計部，北京100027；3.國家海洋局第一海洋研究所，山東青島266100)

21世紀是全面開發(fā)利用和保護海洋的世紀。中國具有豐富的油氣儲量，海洋油氣資源已成為我國能源開發(fā)的一個重要方向。海底油氣輸送管道是海上油氣田開發(fā)中油氣輸送的主要方式[1]。然而隨著沿海地區(qū)經(jīng)濟迅速發(fā)展，海洋災(zāi)害造成的經(jīng)濟損失日益增多。海底滑坡是一種常見的海洋災(zāi)害，應(yīng)當對海底滑坡的機理研究和災(zāi)害預(yù)防引起足夠的重視[2]。我國正逐步加大對海洋科學的投入，1996年與法國合作開展東海陸坡的穩(wěn)定性研究項目[3]，近年將南海作為“十二五”期間深海研究的重點海域。然而由于客觀條件的限制，對海底滑坡的認識定性描述多，定量研究少且不夠準確。今后應(yīng)著重針對滑坡土體性質(zhì)、失穩(wěn)區(qū)影響范圍預(yù)測以及斜坡失穩(wěn)機理等方面進行研究。

強度折減有限元法[4]是一種數(shù)值分析方法，以描述滑坡體體內(nèi)應(yīng)力、應(yīng)變特征的本構(gòu)模型為基礎(chǔ)分析邊坡的變形和穩(wěn)定問題。眾多學者針對這種方法做出各類探索[4-7]，研究表明，強度折減有限元法可有效分析邊坡的穩(wěn)定性[6]。但是這些研究多是針對陸地邊坡，很少涉及海底斜坡不穩(wěn)定性。本文運用ABAQUS二次開發(fā)穩(wěn)定性分析模塊，應(yīng)用強度折減技術(shù)探討海底斜坡坡度、土體不排水強度、土體密度、摩擦角等參數(shù)對海底斜坡穩(wěn)定性的影響。本模塊并不試圖追蹤海底滑坡的完整過程，而是將強度折減有限元法應(yīng)用到海底邊坡穩(wěn)定評價，為將來的海底滑坡模擬儲備先期技術(shù)。

1 強度折減理論

強度折減法基本過程是不斷對巖土體抗剪強度參數(shù)粘聚力 c和內(nèi)摩擦角φ進行折減:c*=c/F，tanφ*=tanφ/F，直到邊坡處于失穩(wěn)極限破壞狀態(tài)，此時對應(yīng)的折減系數(shù)F為邊坡穩(wěn)定的強度儲備安全系數(shù)[7]。安全系數(shù)搜索列表中前一個折減系數(shù)一般作為邊坡穩(wěn)定安全系數(shù)。

失穩(wěn)判據(jù)的選擇在強度折減法中存在爭議，現(xiàn)有的邊坡失穩(wěn)判據(jù)大致可以分為四類:①折減后的土體強度參數(shù)使得有限元計算在規(guī)定的迭代次數(shù)內(nèi)不能收斂[8-9]；②塑性區(qū)從坡腳貫通至坡頂[10]；③超過某一幅值的廣義剪應(yīng)變或等效塑性應(yīng)變形成貫通帶[11]。④以邊坡坡面特征點位移陡增作為失穩(wěn)判據(jù)[12]。本文強度折減計算中，采用數(shù)值迭代不收斂作為邊坡失穩(wěn)判據(jù)。

2 基于強度折減理論ABAQUS二次開發(fā)實現(xiàn)自動搜索安全系數(shù)

2.1 ABAQUS

本文所研究的深?；滤苄詤^(qū)開展的過程中會產(chǎn)生非常強烈的非線性并伴隨土體弱化，因此采用該有限元軟件是合理的。同時ABAQUS攜帶有功能強大的用戶二次開發(fā)接口，并內(nèi)置了腳本語言Python[13]，ABAQUS可以通過編寫Python腳本直接操縱ABAQUS內(nèi)核，實現(xiàn)ABAQUS自動前處理過程?；诖藘?yōu)勢發(fā)展具備安全系數(shù)自動搜索功能的邊坡穩(wěn)定性分析模型，計算海底二維邊坡的安全系數(shù)，并輸出邊坡失穩(wěn)破壞模式，為輔助決策系統(tǒng)提供相關(guān)技術(shù)支持。

2.2 模塊組織結(jié)構(gòu)

模塊調(diào)用有限元軟件ABAQUS(Version 6.5)進行數(shù)值模擬，整體框架由Fortran主程序和Python用戶程序兩部分組成，體系架構(gòu)如圖1所示。

2.3 Fortran主程序

本模塊采用兩分法自動搜索安全系數(shù)，搜索流程如圖2所示，針對每一折減系數(shù)，主控制程序調(diào)用其他3個子模塊進行計算，以判斷是否在規(guī)定的迭代次數(shù)內(nèi)滿足收斂條件和精度要求。

2.4 Python用戶程序

開發(fā)的Python用戶程序為Slope.py，實際是ABAQUS/CAE建模流程的程序化實現(xiàn)，用戶只需要修改幾何外形、土體本構(gòu)和網(wǎng)絡(luò)剖分密度等少數(shù)參數(shù)，便可自動創(chuàng)建有限元模型。將修改好的用戶程序文件放置到主程序SSFEA3.exe所在的文件夾里，雙擊SSFEA3.exe運行，運動期間不需要任何干預(yù)便可自動完成計算。

圖1 模塊的體系架構(gòu)

圖2 兩分法自動搜索安全系數(shù)流程

3 算例驗證

取有代表性的土坡[9]為算例1。均質(zhì)土坡如圖3所示(以 β=40°為例)，坡高 H=20 m，土的重度 γ=20 kN/m3，彈性模量 E=100 MPa，泊松比 ν=0.3，粘聚力c=42 kPa，內(nèi)摩擦角 φ=17°。求坡角 β=30°，35°，40°，45°，50°時邊坡的穩(wěn)定安全系數(shù) 。

圖3 算例1邊坡幾何外形(m)

3.1 有限元模型建立

按照平面應(yīng)變建立模型，邊坡體左右邊界為水平約束，底面為水平及垂直約束，其他為自由約束。土的本構(gòu)模型均選用基于Mohr-Coulomb屈服準則和非關(guān)聯(lián)流動法則(剪脹角ψ=0)的理想彈塑性模型。采用二階縮減積分單元剖分網(wǎng)格，整個計算區(qū)域劃分為四邊形單元。一次性施加重力荷載，即荷載分析步設(shè)置為1。折減系數(shù) F的精度要求達到0.001，最大迭代次數(shù)為500次。

3.2 安全系數(shù)計算結(jié)果

表1為各屈服準則采用非關(guān)聯(lián)流動法則(膨脹角 ψ=0)時的安全系數(shù)結(jié)果[9]，表1中，DP1是外接圓DP準則，DP2是摩爾-庫侖等面積圓DP準則，DP3是摩爾-庫侖匹配DP準則(平面應(yīng)變條件下)，S是Spencer法。

表1 用不同方法求得的穩(wěn)定安全系數(shù)

從表1可以看出，采用平面應(yīng)變條件下自動搜索程序計算的安全系數(shù)與傳統(tǒng)Spencer法、摩爾-庫侖匹配DP準則(DP3)求得的安全系數(shù)非常接近，誤差在1%左右，而采用其他方法的計算結(jié)果誤差較大，摩爾-庫侖等面積圓DP準則(DP2)的計算結(jié)果比Spencer法計算的結(jié)果大約6%。外接圓DP準則條件(DP1)下的安全系數(shù)比傳統(tǒng)極限平衡法大約25%，結(jié)果證明，在各種坡度下，本方法的計算結(jié)果都是符合實際的，是一種有效的計算方法。

4 海底斜坡算例及計算結(jié)果

基于工程實例，建立算例2。均質(zhì)土坡坡長 L=100 m，坡角為5°，材料參數(shù)為飽和密度 ρ=1 521 kg/m3，不排水強度 c=8.5 kPa，彈性模量 E=3 MPa，泊松比 ν=0.49，內(nèi)摩擦角忽略不計，外荷載為重力荷載。邊坡幾何尺寸及有限元模型如圖4(a)所示。邊界約束條件、本構(gòu)模型、單元類型同算例1，圖4(b)為有限元網(wǎng)格圖。折減系數(shù)F的精度達到0.001，表2給出了安全系數(shù)自動搜索過程，最終確定斜坡安全系數(shù)為0.995。

圖4 邊坡幾何外形尺寸及網(wǎng)格剖分

表2 安全系數(shù)范圍搜索過程

從圖5可以清楚看到隨著分析步的增加，荷載逐步施加，土體中塑性應(yīng)變從坡頂向坡角逐步發(fā)展，當達到臨界破壞狀態(tài)時，塑性區(qū)已發(fā)展到基坑頂面，形成經(jīng)典的圓弧狀貫通區(qū)域，此滑裂面的形狀與極限平衡法中的圓弧是一致的。

本文建立了大量計算模型，探討各海底斜坡土體參數(shù)對海床穩(wěn)定性的影響。上述算例中其他土體材料參數(shù)不變，分別計算斜坡坡度為 2.0°、2.5°、3.0°、3.2°、3.3°、3.4°、3.5°、4.0°、5.0°時的安全系數(shù)，計算結(jié)果見表3。從表3中可以看出，隨著坡角增加安全系數(shù)減小，土坡首先局部發(fā)生破壞，隨著土體應(yīng)力的轉(zhuǎn)移，土坡逐漸破壞，直至最后達到失穩(wěn)臨界狀態(tài)。土體坡角為2.00°時該區(qū)段海底斜坡是穩(wěn)定的，如果坡角增加5.00°，斜坡處于臨界狀態(tài)。

4.1 強度參數(shù)的影響

取用算例2，泊松比 ν=0.3，分別取摩擦角 8°，膨脹角4°和摩擦角2°，膨脹角1°這兩種情況進行計算。圖6顯示了斜坡臨界坡度與土體摩擦角的關(guān)系，定義安全系數(shù)約為1時的斜坡坡角為臨界坡度。從圖6中可以看出，摩擦角對安全系數(shù)影響較大，同樣的工況條件下摩擦角越大，對應(yīng)的安全系數(shù)越大。同時安全系數(shù)隨著坡度增大而變小，即斜坡坡度越大，其穩(wěn)定性越差，這是符合常規(guī)規(guī)律的。

圖5 塑性應(yīng)變分布隨分析步的變化

表3 土體不排水強度為C=8.5 kPa時的坡角與安全系數(shù)關(guān)系

圖7顯示了在兩種摩擦角取值情況下，斜坡臨界坡度與土體不排水強度的關(guān)系，可以看出土體不排水強度越大，斜坡達到臨界破壞狀態(tài)的臨界坡度越大，不排水強度是斜坡穩(wěn)定的有利因素。

圖6 臨界坡度與摩擦角的關(guān)系

圖7 臨界坡度與不排水強度的關(guān)系

圖8顯示了在兩種摩擦角取值情況下，斜坡臨界坡度與土體密度的關(guān)系，可以看出，土體密度大小對斜坡體穩(wěn)定性影響不大。土體密度對邊坡穩(wěn)定性的影響是一分為二的。土體密度越大即土體容重越大，一方面沿滑面切向的下滑力增加，這會誘發(fā)滑坡，另一方面，垂直于滑面方向的正應(yīng)力也隨之增加，即增加了滑面上的抗滑力，有利于邊坡穩(wěn)定。因而不能簡單的評價土體密度大小對邊坡穩(wěn)定性的利弊，而應(yīng)根據(jù)具體的工程地質(zhì)條件和邊坡形狀(坡度、坡高等)進行綜合評判。

圖8 臨界坡度與土體密度的關(guān)系圖

4.2 海底斜坡的穩(wěn)定性定量估計方法

對于海底土體，摩擦角幾乎可以忽略不計，土體不排水強度以及斜坡坡角就成了海底斜坡穩(wěn)定性的決定因素，經(jīng)過大量模擬計算，得到土體在不同不排水強度下坡角與安全系數(shù)的關(guān)系圖9。本文為判斷海底斜坡的穩(wěn)定性提供了一種定量估計方法，只要獲取海底某段土體的不排水強度、斜坡坡度，便可以根據(jù)圖9初步估計該段海底坡段的穩(wěn)定性情況。

圖9 不排水強度和臨界坡度組合

通過圖9的分析，可以得出如下規(guī)律:

(1)當不排水強度固定時，坡角越小，邊坡穩(wěn)定安全系數(shù)越大，符合邊坡穩(wěn)定的常規(guī)規(guī)律；如不排水強度為25 kPa時，安全系數(shù)取1.0時，坡角為14.7°，要使安全系數(shù)達到2.0以上，坡角需為7.3°以下。

(2)對應(yīng)固定的安全系數(shù)，土體不排水強度越大，所對應(yīng)的臨界坡角就越大。如取安全系數(shù)為1.0，不排水強度為15 kPa時，斜坡坡角大于8.8°才能穩(wěn)定；當不排水強度為35 kPa時，臨界坡角達到20.5°。

(3)對應(yīng)固定的安全系數(shù)，土體不排水強度與穩(wěn)定坡角基本呈線性對應(yīng)關(guān)系。同時安全系數(shù)越大，穩(wěn)定坡角與不排水強度的比值(斜率)越小。如不排水強度從25 kPa增大到45 kPa:安全系數(shù)為1時，關(guān)系曲線斜率約為0.49，安全系數(shù)為1.5時，關(guān)系曲線斜率約為0.36；安全系數(shù)為2時，關(guān)系曲線斜率約為0.285。

5 結(jié) 論

通過以上分析，可以得出如下的結(jié)論:

(1)運用基于ABAQUS的二次開發(fā)模塊分析海底斜坡穩(wěn)定性，實現(xiàn)了強度折減法的自動搜索安全系數(shù)功能，極大地免去了手工操作的麻煩。

(2)坡體的c，φ值對邊坡穩(wěn)定性影響顯著，不論對陸地邊坡還是海底斜坡都是敏感性因素。土體強度越大、摩擦角越大，邊坡越穩(wěn)定。對海底斜坡穩(wěn)定性研究中，往往視摩擦角為零，土體不排水強度起著決定性作用，影響最大。土體密度對海底斜坡的穩(wěn)定性影響不大。

(3)對于同一安全系數(shù)，土體不排水強度與對應(yīng)坡角基本呈線性對應(yīng)關(guān)系。同時安全系數(shù)越大，臨界狀態(tài)對應(yīng)坡角與不排水強度的比值(斜率)越小。

(4)可以根據(jù)圖9初步估計該段海底坡段的穩(wěn)定性情況，其結(jié)果與工程實踐是相符合的。

[1]周延東，劉日柱.我國海底管道的發(fā)展狀況與前景[J].中國海上油氣(工程)，1998，10(4):1-5.

[2]吳時國，陳珊珊，王志君，等.大陸邊緣深水區(qū)海底滑坡及其不穩(wěn)定性風險評估[J].現(xiàn)代地質(zhì)，2008，22(3):430-437.

[3]YIN P，BER NE B，VA GNER P，et al.Mud volcanoes at the shelf margin of the East China Sea[J].Marine Geology，2003，194(3-4):135-149.

[4]鄭穎人，陳祖煜，王恭先，等.邊坡與滑坡工程治理[M].北京:人民交通出版社，2009.

[5]秦 帆，王正中.基于有限元強度折減理論的邊坡穩(wěn)定分析方法探討與改進[J].水利與建筑工程學報，2012，10(1):43-47.

[6]萬少石，年廷凱，蔣景彩，等.邊坡穩(wěn)定強度折減有限元分析中的若干問題討論[J].巖土力學，2010，31(7):2283-2288.

[7]鄭穎人，趙尚毅，王恭先，等.邊(滑)坡工程設(shè)計中安全系數(shù)的討論[J].巖石力學與工程學報，2006，25(9):1937-1940.

[8]Zienkiewicz O C，Humpheson C，Lewis R W.Associated and non-associated viscoplasticity and plasticity in soil mechanics[J].Geotechnique，1975，25(4):671-689.

[9]趙尚毅，鄭穎人，張玉芳.極限分析有限元法講座-Ⅱ有限元強度折減法中邊坡失穩(wěn)的判據(jù)探討[J].巖土力學，2005，26(2):332-336.

[10]Matsui T，San K-C.Finite element slope stability analysis by shear strength reduction technique[J].Soils and Foundations，1992 ，32(1):59-70.

[11]鄭 宏，李春光，李焯芬，等.求解安全系數(shù)的有限元法[J].巖土工程學報，2002，24(5):626-628.

[12]遲世春，關(guān)立軍.基于強度折減的拉格朗日差分方法分析土坡穩(wěn)定[J].巖土工程學報，2004，26(1):42-46.

[13]鐘同圣，衛(wèi) 豐，王 鷙，等.Python語言和ABAQUS前處理二次開發(fā)[J].鄭州大學學報，2006，38(1):60-64.