基于數(shù)值模擬的滑帶土蠕變特性研究

李 翔，程 聰

（中國地質(zhì)大學(xué)工程學(xué)院，湖北 武漢430074）

滑坡是一種重要的地質(zhì)災(zāi)害，給人類的生命財(cái)產(chǎn)帶來重大威脅?；瑤腔伦钪匾慕Y(jié)構(gòu)，滑坡動(dòng)態(tài)演化規(guī)律的三個(gè)階段：蠕動(dòng)緩滑階段、等速變形階段和加速滑動(dòng)階段的內(nèi)因直接與滑帶土的蠕變特性有關(guān)，而滑帶土的長期強(qiáng)度對(duì)滑坡的穩(wěn)定性起著至關(guān)重要的控制作用。目前室內(nèi)巖石三軸蠕變?cè)囼?yàn)主要用于滑坡滑帶土的蠕變特性及其長期穩(wěn)定性研究，通過大量的試驗(yàn)可以揭示滑坡的蠕變力學(xué)屬性，確定其參數(shù)，并建立力學(xué)模型［1—5］。但是由于室內(nèi)巖石蠕變?cè)囼?yàn)費(fèi)時(shí)長、耗資大，在實(shí)際操作中比較困難。而數(shù)值模擬技術(shù)的經(jīng)濟(jì)性和快速性，使巖石蠕變?cè)囼?yàn)數(shù)值模擬得到了廣泛的應(yīng)用［6—8］。因此，通過實(shí)驗(yàn)技術(shù)及蠕變理論分析，確定滑帶土的蠕變特性，建立其蠕變模型，進(jìn)而對(duì)滑坡演化趨勢進(jìn)行分析，具有重要的理論意義和應(yīng)用價(jià)值。

1 滑帶土的室內(nèi)三軸蠕變?cè)囼?yàn)

本次三軸蠕變?cè)囼?yàn)的滑帶土取自位于黃土坡臨江Ⅱ號(hào)崩滑堆積體上TP4平硐內(nèi)，滑帶土樣品的物質(zhì)組成為粉質(zhì)黏土夾（含）碎石（角礫），其中黏土礦物以蒙脫石、高嶺石、綠泥石、伊利石等親水礦物為主。為了研究黃土坡滑坡滑帶土長期蠕變特性，選擇了Burger’s蠕變方程，并對(duì)參數(shù)進(jìn)行辨識(shí)，同時(shí)利用FLAC 3D內(nèi)置的cvisc模型進(jìn)行了三軸蠕變?cè)囼?yàn)數(shù)值模擬，研究了滑帶土的蠕變特性及其長期強(qiáng)度。

1.1 室內(nèi)三軸蠕變?cè)囼?yàn)

黃土坡滑坡滑帶土的三軸蠕變?cè)囼?yàn)在中科院武漢巖土所經(jīng)過改裝的三軸流變儀上進(jìn)行。軸向加載采用砝碼加載，并保持原有的圍壓系統(tǒng)、反壓系統(tǒng)和孔隙水壓力系統(tǒng)不變，加載增量一般根據(jù)土樣的常規(guī)力學(xué)指標(biāo)推算出土樣達(dá)到破壞時(shí)的應(yīng)力值，并確定蠕變?cè)囼?yàn)的加載分級(jí)，一般取n＝4～8，則每級(jí)應(yīng)力水平的增量為（qf/n），其中qf為在一定的圍壓σ3下的破壞應(yīng)力值，加載方法采用陳氏加載法；固結(jié)壓力一般根據(jù)滑帶土的埋深來確定，垂向固結(jié)壓力σ3分別設(shè)定為200kPa、300kPa、400kPa。本次三軸固結(jié)不排水蠕變?cè)囼?yàn)的穩(wěn)定標(biāo)準(zhǔn)規(guī)定為24h內(nèi)變形量小于0.01mm，則可進(jìn)入下一級(jí)應(yīng)力水平試驗(yàn)。

在黏彈性體假設(shè)的前提下，利用Boltzmann迭加原理對(duì)蠕變?cè)囼?yàn)原始數(shù)據(jù)進(jìn)行處理，可以得到各級(jí)剪應(yīng)力作用下的剪應(yīng)變（ε）隨時(shí)間（t）變化的蠕變曲線，即分別加載蠕變曲線，詳見圖1至圖3。

圖1 圍壓σ3為200kPa時(shí)各級(jí)剪應(yīng)力下滑帶土的蠕變曲線Fig.1 Results of creep experiments under differentσ levels with 200kPaσ3

由圖1至圖3可以看出：當(dāng)剪應(yīng)力較小時(shí)，剪應(yīng)變速率逐漸減小，最后趨向于零，試樣主要以瞬時(shí)變形為主，而沒有達(dá)到破壞，即發(fā)生衰減蠕變；隨著剪切力的增加，試樣開始出現(xiàn)等速蠕變階段，最后發(fā)生非衰減蠕變，試樣發(fā)生破壞。

此外，從各蠕變曲線可見，在任何一個(gè)剪切荷載作用下，試樣都存在著瞬時(shí)變形，因此可使用一個(gè)彈性元件來模擬該現(xiàn)象。當(dāng)剪應(yīng)力較低（小于長期強(qiáng)度）時(shí)，隨著時(shí)間的增加，剪應(yīng)變速率逐漸降低，近似呈負(fù)指數(shù)的形式趨于某一漸近線，這種性質(zhì)可以采用Kelvin體來模擬。將彈性元件和Kelvin體串聯(lián)起來，可以反映滑帶土蠕變變形的瞬時(shí)變形和衰減變形階段。當(dāng)剪應(yīng)力大于長期強(qiáng)度時(shí)，試樣將出現(xiàn)等速蠕變，并最終進(jìn)入加速蠕變階段而破壞，因此在上述模型的基礎(chǔ)上再串聯(lián)一個(gè)黏性元件就可以模擬滑帶土的等速蠕變階段。這樣，模擬瞬時(shí)變形的彈性元件和模擬等速蠕變的黏性元件串聯(lián)，則構(gòu)成了Maxwell體。

圖2 圍壓σ3為300kPa時(shí)各級(jí)剪應(yīng)力下滑帶土的蠕變曲線Fig.2 Results of creep experiments under differentσ levels with 300kPaσ3

圖3 σ3為400kPa時(shí)各級(jí)剪應(yīng)力下滑帶土的蠕變曲線Fig.3 Results of creep experiments under differentσ levels with 400kPaσ3

1.2 Burger’s模型及參數(shù)辨識(shí)

Burger’s模型是由Kelvin體和 Maxwell體串聯(lián)而成，如圖4所示。該模型的蠕變方程為

式中：ε（t）為三軸蠕變?cè)囼?yàn)過程中某級(jí)圍壓σ3及剪應(yīng)力σ作用下的軸向應(yīng)變；E1、E2、β1、β2分別為Maxwell和Kelvin單元的彈性模量（kPa）及黏滯系數(shù)（kPa·min）。

在剪切蠕變?cè)囼?yàn)中，對(duì)于每一個(gè)固定剪切荷載值τ，在固定的圍壓下可將式（1）改寫為如下形式：

圖4 Burger’s模型Fig.4 Burger’s model

式（2）可以進(jìn)一步簡化為

令t＝0，式（3）右端等于A，即可以根據(jù)瞬時(shí)變形來確定A；當(dāng)時(shí)間足夠大時(shí)，式（3）右端最后一項(xiàng)趨于常數(shù)C，則式（3）可看作是直線方程，直線斜率即為B。當(dāng)土體中的應(yīng)力小于長期強(qiáng)度時(shí)將產(chǎn)生衰減蠕變，此時(shí)蠕變過程可能減速進(jìn)行，蠕變速率逐漸減小，最后趨向于零，變形值也趨向于穩(wěn)定值。由于衰減蠕變過程中沒有等速蠕變階段，因此式（3）可簡化為

將上式中常數(shù)項(xiàng)移至等式左邊，在等式兩邊取對(duì)數(shù)，該式在對(duì)數(shù)坐標(biāo)系中也為一直線方程，A已確定，D為直線的斜率，再將D值代入式（4）即可解出C。由此可確定一組初值，然后在 Matlab中利用lsqcurvefit非線性回歸工具可回歸出A、B、C、D的穩(wěn)定解，這樣Burger’s模型所有的參數(shù)均可獲得。按上述方法，得到了各固結(jié)壓力和不同剪切力作用下的Burger’s模型的參數(shù)，詳見表1。

表1 Burger’s模型的參數(shù)Table 1 Parameters of Burger’s model

根據(jù)蠕變?cè)囼?yàn)數(shù)據(jù)，本文僅對(duì)圍壓為200kPa時(shí)各級(jí)剪應(yīng)力下的蠕變曲線進(jìn)行分析，Burger＇s模型的擬合曲線見圖5和圖6。由圖5和圖6可以看出：Burger’s模型能較好地?cái)M合出試樣在低剪應(yīng)力情況下的蠕變特性，即試樣主要以瞬時(shí)變形為主，發(fā)生衰減蠕變，隨著剪應(yīng)力的增加，開始出現(xiàn)等速蠕變階段。

圖5 圍壓為200kPa、剪應(yīng)力為30kPa時(shí)Burger’s模型的擬合曲線Fig.5 Fitting curve of Burger’s model with 200kPaσ3 and 30kPaσ

圖6 圍壓為200kPa、剪應(yīng)力為180kPa時(shí)Burger’s模型的擬合曲線Fig.6 Fitting curve of Burger’s model with 200kPaσ3 and 180kPaσ

2 滑帶土室內(nèi)三軸蠕變?cè)囼?yàn)的數(shù)值模擬

本文利用Burger＇s模型擬合參數(shù)并借助FLAC 3D對(duì)室內(nèi)三軸蠕變?cè)囼?yàn)進(jìn)行了數(shù)值模擬研究。模型的運(yùn)行條件與實(shí)際的室內(nèi)試驗(yàn)條件相同，三軸試樣尺寸為：d＝3.91cm，h＝8cm。模型體為圓柱體，網(wǎng)格采取均勻劃分，劃分的網(wǎng)格為10×10×10，單元底面為豎向位移約束，周圍采用應(yīng)力邊界條件；模型頂部荷載為分別加載，荷載為實(shí)際室內(nèi)試驗(yàn)的多級(jí)加載值。蠕變模型選用FLAC 3D內(nèi)置的cvisc模型，即Burger’s（伯格斯）蠕變模型和摩爾－庫侖模型合成的伯格斯蠕變粘塑性模型，且采用蠕變大變形模式，根據(jù)室內(nèi)三軸蠕變?cè)囼?yàn)的結(jié)果以及室內(nèi)試驗(yàn)時(shí)各級(jí)荷載下軸向位移量的大小，對(duì)滑帶土試樣的蠕變參數(shù)進(jìn)行必要的調(diào)整后取值，并通過對(duì)室內(nèi)三軸蠕變?cè)囼?yàn)進(jìn)行數(shù)值模擬，可以得出不同圍壓、不同軸向應(yīng)力時(shí)滑帶土分別加載的蠕變模擬曲線，見圖7至圖9。

圖7 圍壓σ3為200kPa時(shí)滑帶土分別加載的蠕變模擬曲線Fig.7 Simulation curves of the creep of the slip soil under different loads with 200kPaσ3

圖8 圍壓σ3為300kPa時(shí)滑帶土分別加載的蠕變模擬曲線Fig.8 Simulation curves of the creep of the slip soil under different loads with 300kPaσ3

圖9 圍壓σ3為400kPa時(shí)滑帶土分別加載的蠕變模擬曲線Fig.9 Simulation curves of the creep of the slip soil under different loads with 400kPaσ3

由圖7至圖9可見，隨著圍壓的增大，滑帶土蠕變強(qiáng)度具有明顯的提高。較低軸向應(yīng)力水平下，試件存在明顯的瞬時(shí)變形，而后產(chǎn)生等速蠕變；較高軸向應(yīng)力水平下，試件沒有明顯的瞬時(shí)變形，而是一開始即產(chǎn)生等速蠕變，之后巖體的應(yīng)變速率減小，隨著時(shí)間的增加，其變形趨于穩(wěn)定或出現(xiàn)流動(dòng)，直至試件發(fā)生加速蠕變破壞。

3 滑帶土的長期強(qiáng)度

土的非衰減蠕變的發(fā)展可引起具有加速特征的流動(dòng)，并最后導(dǎo)致土體破壞，因此土的長期破壞強(qiáng)度可能小于短期荷載作用下的強(qiáng)度值。在實(shí)際工程中，一些線性土坡其安全系數(shù)達(dá)到或甚至超過設(shè)計(jì)要求，然而在施工完成若干年后，其強(qiáng)度卻降低并產(chǎn)生破壞，因此土的長期強(qiáng)度越來越引起人們的重視。在實(shí)驗(yàn)室一般采用以下方法來確定土的長期強(qiáng)度：一種方法是將梯級(jí)加載的試驗(yàn)結(jié)果在對(duì)數(shù)坐標(biāo)下繪制應(yīng)力－應(yīng)變曲線，找到對(duì)數(shù)坐標(biāo)中曲線上有明顯拐點(diǎn)的應(yīng)力作為長期強(qiáng)度極限值τ∞，這是因?yàn)樵诜撬p蠕變階段變形值急劇增長的結(jié)果；另一種方法是根據(jù)流動(dòng)速度來判斷，長期強(qiáng)度極限值是區(qū)別衰減蠕變和非衰減蠕變的界限值，因此可以通過繪制應(yīng)力和該應(yīng)力下的應(yīng)變速率之間的關(guān)系曲線來判斷，曲線和應(yīng)力軸的交點(diǎn)就是長期強(qiáng)度極限值τ∞。

本文采用等時(shí)曲線來確定滑帶土的長期強(qiáng)度，即為某時(shí)刻的應(yīng)力－應(yīng)變曲線。首先從數(shù)值模擬試驗(yàn)中導(dǎo)出滑帶土分別加載蠕變曲線數(shù)據(jù)，如圖10右半部分τ1～τ6等荷載作用下的分別加載蠕變曲線；接著利用Excel，將對(duì)應(yīng)于t4時(shí)間的各級(jí)應(yīng)力作用下的時(shí)間－應(yīng)變點(diǎn)A、B、C、D、E投影到圖10左半部分，并將橫坐標(biāo)轉(zhuǎn)換成對(duì)應(yīng)的剪應(yīng)力，然后在新的坐標(biāo)軸下將投影點(diǎn)連接起來得出其等時(shí)曲線。當(dāng)試驗(yàn)時(shí)間足夠長時(shí)，可以通過t∞直接獲取滑帶土的長期強(qiáng)度；當(dāng)試驗(yàn)時(shí)間較短時(shí)，也可以通過不同時(shí)刻的應(yīng)力－應(yīng)變曲線的拐點(diǎn)連線來確定長期強(qiáng)度。圖11至圖13為不同圍壓、不同時(shí)刻滑帶土的軸向應(yīng)力－軸向應(yīng)變關(guān)系曲線，即等時(shí)曲線。

由圖11至圖13可以看出：不同軸向應(yīng)力下的等時(shí)曲線由近似線性段和非線性段組成，并且有轉(zhuǎn)折點(diǎn)，將各個(gè)轉(zhuǎn)折點(diǎn)連接起來，并利用指數(shù)函數(shù)進(jìn)行擬合，在軸向應(yīng)變達(dá)到0.15時(shí)可以認(rèn)為該試樣在不同圍壓下的長期強(qiáng)度。采用這種方法可以得出：圍壓為200kPa下滑帶土的長期強(qiáng)度為130kPa，圍壓為300kPa下滑帶土的長期強(qiáng)度為180kPa，圍壓為400kPa下滑帶土的長期強(qiáng)度為200kPa。通過繪制莫爾應(yīng)力圓可以得到滑帶土的長期強(qiáng)度參數(shù)：內(nèi)聚力c∞為21.20kPa，內(nèi)摩擦角φ∞為18.56°。

圖10 滑帶土的等時(shí)曲線圖解法示意圖Fig.10 Sketch map of isochronal curve graphical method

圖11 圍壓σ3為200kPa時(shí)滑帶土的等時(shí)曲線Fig.11 Isochronal curves of the slip soil with 200kPaσ3

圖12 圍壓σ3為300kPa時(shí)滑帶土的等時(shí)曲線Fig.12 Isochronal curves of the slip soil with 300kPaσ3

圖13 圍壓σ3為400kPa時(shí)滑帶土的等時(shí)曲線Fig.13 Isochronal curves of the slip soil with 400kPaσ3

4 結(jié) 論

（1）本文以室內(nèi)三軸蠕變?cè)囼?yàn)結(jié)果為依據(jù)，通過對(duì)Burger＇s蠕變方程中各參數(shù)進(jìn)行辨識(shí)，從而為描述黃土坡滑坡滑帶土的力學(xué)特性，特別是其長期蠕變特征打下了一定的基礎(chǔ)。

（2）數(shù)值模擬試驗(yàn)結(jié)果與室內(nèi)三軸蠕變?cè)囼?yàn)結(jié)果具有很好的相似性。數(shù)值模擬試驗(yàn)主要依托三軸蠕變?cè)囼?yàn)結(jié)果，在滿足其誤差的基礎(chǔ)上，利用數(shù)值模擬，試驗(yàn)結(jié)果可以獲取滑帶土分別加載下的等時(shí)曲線，從而得到滑帶土的長期強(qiáng)度參數(shù)：內(nèi)聚力c∞為21.20kPa，內(nèi)摩擦角φ∞為18.56°。

［1］汪斌，朱杰兵，唐輝明，等.黃土坡滑坡滑帶土的蠕變特性研究［J］.長江科學(xué)院院報(bào)，2008，25（1）：49－52.

［2］嚴(yán)紹軍，項(xiàng)偉，唐輝明.大巖淌滑坡滑帶土蠕變性質(zhì)研究［J］.巖土力學(xué)，2008，29（1）：58－63.

［3］鄒良超，王世梅.古樹包滑坡滑帶土蠕變經(jīng)驗(yàn)?zāi)Ｐ停跩］.工程地質(zhì)學(xué)報(bào)，2011，19（1）：59－64.

［4］曹運(yùn)江，黃潤秋.某水電站高邊坡煤系軟弱結(jié)構(gòu)面流變特性試驗(yàn)研究［J］.巖石力學(xué)與工程學(xué)報(bào)，2008，27（增2）：3732－3740.

［5］王琛，唐明.三峽泄灘滑坡滑動(dòng)帶土的Singh－Mitchell蠕變方程［J］.四川大學(xué)學(xué)報(bào)（工程科學(xué)版），2003，35（5）：93－95.

［6］王旭東，付小敏.壓縮蠕變力學(xué)試驗(yàn)的數(shù)值模擬研究［J］.工程地質(zhì)學(xué)報(bào)，2009，17（4）：533－537.

［7］鄭俊.某滑坡滑帶土流變模型構(gòu)建及FLAC3D二次開發(fā)［D］.宜昌：三峽大學(xué)，2009.

［8］陳育民，劉漢龍，鄧肯·張.本構(gòu)模型在FLAC3D中的開發(fā)與實(shí)現(xiàn)［J］.巖土力學(xué)，2007，28（10）：2123－2126.

［9］Itasca Consulting Group.FastLagrangianAnalysisofContinua in3Dimensions［M］.MN.USA：Itasca Consulting Group Inc.，Minneapolis，2002.

［10］周德培等.流變力學(xué)原理及其在巖土工程中的應(yīng)用［M］.成都：西南交通大學(xué)出版社，1995.