徐 晗，饒錫保，汪明元

（長江科學(xué)院水利部巖土力學(xué)與工程重點(diǎn)實(shí)驗(yàn)室，武漢 430010）

1 概 述

膨脹巖邊坡失穩(wěn)是工程中經(jīng)常遇到并難于處理的問題，正在設(shè)計(jì)的南水北調(diào)中線工程的膨脹巖土渠段長達(dá)340 km，如何正確進(jìn)行邊坡的穩(wěn)定性分析是合理設(shè)計(jì)的前提。在實(shí)際工程中，膨脹巖邊坡大多是在降雨時發(fā)生破壞的，因此有必要研究降雨條件下膨脹巖邊坡的失穩(wěn)過程。目前數(shù)值模擬降雨入滲引起邊坡失穩(wěn)較為常見的方法是先用滲流理論求解瞬態(tài)的滲流場，再用技術(shù)已經(jīng)成熟的極限平衡法進(jìn)行土坡暫態(tài)穩(wěn)定性分析［1－4］；但降雨條件下邊坡的失穩(wěn)是一個應(yīng)力場和滲流場不斷變化、相互作用的漸進(jìn)過程，因此，在考慮降雨入滲邊坡穩(wěn)定性分析中，引入土體的本構(gòu)關(guān)系，用非線性有限元流固耦合的分析方法模擬邊坡的穩(wěn)定性，具有重要的意義。

經(jīng)過長時間的研究與探討，作者認(rèn)為對膨脹巖邊坡的穩(wěn)定性數(shù)值分析應(yīng)正確模擬膨脹巖的膨脹性、裂隙性，以及非飽和膨脹巖的吸濕軟化特性。目前多數(shù)學(xué)者在膨脹巖邊坡穩(wěn)定分析技術(shù)中，不考慮其非飽和強(qiáng)度基本特性及膨脹性、裂隙性［5］，使得計(jì)算結(jié)果偏離實(shí)際，因此，有必要對膨脹巖邊坡穩(wěn)定技術(shù)進(jìn)行改進(jìn)，使之能考慮膨脹巖與普通土的不同之處，揭示膨脹巖邊坡破壞的淺層性與漸進(jìn)性特點(diǎn)。為此，采用流固耦合的非線性有限元分析方法，對膨脹巖邊坡降雨失穩(wěn)現(xiàn)場試驗(yàn)進(jìn)行了數(shù)值模擬研究，改進(jìn)了膨脹巖邊坡穩(wěn)定分析技術(shù)，使之能考慮膨脹巖的膨脹性、裂隙性，以及非飽和膨脹巖的吸濕軟化特性。

2 流固耦合的分析方法

采用Galerkin有限元格式，將節(jié)點(diǎn)位移和孔隙水壓力作為節(jié)點(diǎn)自由度進(jìn)行空間離散，則力學(xué)平衡方程及滲流連續(xù)方程可寫成如下矩陣形式［6］：

平衡方程為

式中：［K］為通常的剛度矩陣；｛Δ珋δ｝為位移增量；［L］為節(jié)點(diǎn)孔隙水壓力所對應(yīng)的節(jié)點(diǎn)力；｛Δ珋p｝為孔隙水壓力增量；｛F｝為節(jié)點(diǎn)外荷載；｛I｝為增量迭代過程中上一增量步中的不平衡力。

滲流連續(xù)方程為

直接解耦合方程，引入時間積分，差分公式為

式中：ξ為一系數(shù)，0≤ξ≤1；Δt為時間增量。

為了確保數(shù)值穩(wěn)定，可以選擇ξ＝1（后向差分法），此時方程（3）可寫為

因此，在t＋Δt時刻的滲流連續(xù)方程為

其中｛R｝＝Δt［－｛Q｝t＋Δt＋］T｛珋v｝t＋Δt＋［H∧］T｛p｝t＋Δt］，為流體體積變化修正量；［B］為節(jié)點(diǎn)變形對應(yīng)的流體體積改變；［H］為孔壓變化對應(yīng)的流體體積改變。

3 膨脹巖的特殊性質(zhì)

3．1 吸濕強(qiáng)度軟化特性

近年來，非飽和土力學(xué)的發(fā)展為定量計(jì)算因水分入滲而引起的土體軟化的強(qiáng)度變化提供了理論基礎(chǔ)。Fredlund（1978年）提出一個抗剪強(qiáng)度公式為

式中：τf為非飽和土抗剪強(qiáng)度＇，φ＇分別為有效凝聚力和有效內(nèi)摩擦角；σn為法向應(yīng)力；φb為隨吸力變化的內(nèi)摩擦角。根據(jù)非飽和膨脹巖三軸試驗(yàn)，在400 kPa的基質(zhì)吸力范圍內(nèi)，非飽和原狀黏土巖的抗剪強(qiáng)度隨著基質(zhì)吸力的增大而線性增加，該線段的斜率為基質(zhì)吸力引起的強(qiáng)度增量tanφb，其中＝20°。

因此，可在400 kPa的基質(zhì)吸力范圍內(nèi)取φb＝20°，但非飽和土的強(qiáng)度不可能隨基質(zhì)吸力無限制提高，設(shè)定當(dāng)基質(zhì)吸力大于400 kPa時，抗剪強(qiáng)度按400 kPa選取。

3．2 土水特征曲線

泥灰?guī)r的土水特征曲線采用室內(nèi)試驗(yàn)成果。水分運(yùn)動數(shù)值分析中的有關(guān)非飽和土參數(shù)將利用實(shí)測的土水特征曲線確定，并采用VG模型進(jìn)行計(jì)算。Van Genuchten模型能夠在較大的水頭范圍內(nèi)表征水分特征數(shù)據(jù)，應(yīng)用較為廣泛

式中：Ks為飽和滲透系數(shù)；Se＝θ－θr／（θs－θr）為有效飽和度；α，m和n為模型的形狀參數(shù)；參數(shù)m和n的關(guān)系可表達(dá)為m＝1－1／n。

將泥灰?guī)r原狀樣非飽和試驗(yàn)得到的土水特征曲線進(jìn)行擬合，所得公式的具體形式為

式中：θr＝0．12，θs＝0．266 5；α＝0．061；m＝0．14；n＝1．163；h的單位為 kPa。

將黏土巖原狀樣非飽和試驗(yàn)得到的土水特征曲線進(jìn)行擬合，所得VG模型公式的具體形式為

式中：θr＝0．12，θs＝0．312 016；α＝0．002 179；m＝0．178 1；n＝1．216 7；h的單位為 kPa。

3．3 吸濕密度增加

對任意時刻計(jì)算出來的單元積分點(diǎn)的飽和度，根據(jù)土的三相比例指標(biāo)換算關(guān)系可以求出該情況下的含水量，然后根據(jù)干密度可以求出該時刻非飽和土的濕密度ρ。

式中：w為含水量；e為孔隙比；Gs為相對密度；Sr為飽和度；ρd為干密度。

3．4 體積應(yīng)力、含水量變化相關(guān)的膨脹性

假設(shè)膨脹應(yīng)變是由于應(yīng)力第一不變量的改變所引起的，膨脹應(yīng)變僅與體積應(yīng)力有關(guān)，而與應(yīng)力偏斜張量無關(guān)，采用GDS公司的動三軸儀器，對中膨脹土試樣開展了三軸膨脹實(shí)驗(yàn)研究，得到了以體積應(yīng)力、含水量變化為主要變量的膨脹應(yīng)變公式，公式如下。

式中：ε為膨脹體積應(yīng)變，按百分?jǐn)?shù)表示；Δw為含水量變化百分?jǐn)?shù)；p為體積應(yīng)力；Pa為大氣壓力。

3．5 裂隙性

膨脹巖的裂隙性是影響膨脹巖邊坡穩(wěn)定失穩(wěn)的最直接最主要的原因，且受含水量影響顯著，在膨脹巖邊坡穩(wěn)定性分析中，可采用用殘余強(qiáng)度作為風(fēng)化層的最主要的強(qiáng)度參數(shù)，這為膨脹巖邊坡穩(wěn)定評價提供了一定的強(qiáng)度儲備［7］。但在邊坡一定深度范圍內(nèi)，如超過大氣影響深度，膨脹巖可視為保持原狀土的性質(zhì)，在選擇強(qiáng)度參數(shù)時，可采用峰值強(qiáng)度。

4 膨脹巖現(xiàn)場試驗(yàn)數(shù)值模擬

4．1 數(shù)值分析的基本思路

（1）根據(jù)現(xiàn)場試驗(yàn)報(bào)告描述建立自重作用下的初始應(yīng)力場，設(shè)定地下水位，建立初始含水量場。

（2）進(jìn)行降雨數(shù)值模擬，期間考慮如下過程（整個過程為連續(xù)計(jì)算）：入滲引起的邊坡含水量變化場；隨著邊坡含水量的變化引起的土體重度增加（由濕重度變?yōu)轱柡椭囟龋?；土體非飽和強(qiáng)度指標(biāo)隨含水量增加而軟化降低；孔壓與應(yīng)力場的耦合作用。

（3）對于膨脹力的考慮，根據(jù)三軸膨脹試驗(yàn)提出的本構(gòu)模型，分階段或者僅在計(jì)算的最后階段考慮膨脹力的影響（當(dāng)計(jì)算達(dá)到某一階段后，根據(jù)當(dāng)前計(jì)算點(diǎn)的應(yīng)力與含水量場的變化，假定在該階段膨脹過程中應(yīng)力場不變，采用三軸膨脹試驗(yàn)提出的本構(gòu)模型計(jì)算膨脹引起的體積應(yīng)變，轉(zhuǎn)化為溫度荷載施加到模型上）。

（4）邊坡穩(wěn)定計(jì)算，對任一降雨時刻，提取該時刻的單元的強(qiáng)度指標(biāo)、重度、孔壓、含水量變化引起的膨脹應(yīng)變，采用強(qiáng)度折減法來計(jì)算邊坡穩(wěn)定。

4．2 現(xiàn)場試驗(yàn)簡介

長江科學(xué)院巖土力學(xué)與工程重點(diǎn)實(shí)驗(yàn)室在新鄉(xiāng)膨脹巖進(jìn)行了大型的膨脹巖降雨失穩(wěn)現(xiàn)場試驗(yàn)研究，該邊坡從坡底到坡頂17 m高，一級馬道以下高9 m。降雨試驗(yàn)主要在一級馬道以下進(jìn)行，坡比為1∶1．5，降雨范圍具體如圖1斜方向箭頭所示。一級馬道以下大致可分為5層，位于頂部為壤土層，厚約0．2 m；下為厚約6．8 m的泥灰?guī)r、其間有0．8 m厚黏土巖夾層，最下層為黏土巖。除表層壤土外，各地層自由膨脹率分別為21%～29%，63%，40%～54%，48%～70%。

圖1 邊坡有限元網(wǎng)格及降雨示意圖Fig．1 Slopemesh and rainfall schematic diagram

該試驗(yàn)從2008年9月24日開始，在3 d的正式監(jiān)測期間，共進(jìn)行了5次降雨過程，現(xiàn)場試驗(yàn)降雨強(qiáng)度為16 mm／h，共歷時16 h，在數(shù)值計(jì)算中假定降雨為連續(xù)不間斷過程。

當(dāng)現(xiàn)場試驗(yàn)累計(jì)降雨約6 h左右時，一級馬道以下開始發(fā)生大面積滑坡，發(fā)生整體的塑性變形破壞，導(dǎo)致表層水平位移急劇增加至12 mm，在滑坡后的第2 d繼續(xù)降雨1 d后，水平位移持續(xù)以較大的幅度增加。

4．3 初始狀態(tài)考慮

在計(jì)算中將地下水位線設(shè)在渠坡底上1 m左右。由于該邊坡開挖后表層裸露在大氣作用下半年后才開始降雨，因此，根據(jù)泥灰?guī)r、黏土巖入滲深度設(shè)定表層2 m范圍內(nèi)左右為裂隙發(fā)育區(qū)，強(qiáng)度參數(shù)按照有裂隙情況下的殘余強(qiáng)度選取。表層2 m以下為未風(fēng)化層，強(qiáng)度參數(shù)按照峰值強(qiáng)度選取。

4．4 計(jì)算參數(shù)取值

計(jì)算參數(shù)均根據(jù)泥灰?guī)r與黏土巖原狀樣室內(nèi)試驗(yàn)統(tǒng)計(jì)選取。如表1與表2所示。

表2 壤土與砂巖的材料參數(shù)Tabale 2 Material parameters of loam and sand rock

5 計(jì)算結(jié)果分析

5．1 初始狀態(tài)

由于該邊坡高度僅17 m，因此可將飽和度初始狀態(tài)設(shè)定為沿地下水位的自重分布，如圖2所示。位移場采用自重作用下的位移分布，水平位移與豎向位移如圖3所示?？芍患夞R道坡面的水平位移約在2 mm左右，豎向位移在1～3 cm之間，處于穩(wěn)定狀態(tài)。

圖2 降雨試驗(yàn)開始前飽和度Fig．2 Saturation degree before rainfall

表1 材料參數(shù)取值（由室內(nèi)試驗(yàn)統(tǒng)計(jì)得到）Table 1 Material parameters（acquired by indoor test）

圖3 降雨試驗(yàn)開始前水平位移、豎向位移Fig．3 Horizontal displacement and vertical displacement before rainfall

5．2 變形分析

為了監(jiān)視邊坡在降雨過程中發(fā)生的變化，在坡面設(shè)置了3個監(jiān)控點(diǎn)，分別在坡腳A、坡面B、坡頂處C。具體示意圖如圖1所示。圖4為降雨后邊坡飽和度分布圖，可知在邊坡2 m范圍內(nèi)接近于飽和狀態(tài)。

圖4 降雨后飽和度Fig．4 Saturation degree after rainfall

圖5 （a）為降雨引起的等效塑性應(yīng)變，圖5（b）為考慮膨脹力后的等效塑性應(yīng)變，可知膨脹力會加劇塑性區(qū)的擴(kuò)散作用。圖6（a）至圖6（d）分別為降雨后未考慮膨脹力與考慮膨脹力位移分布圖，基本規(guī)律為坡頂水平位移最小，坡底水平位移最大，坡面中間漸變?？芍纯紤]膨脹力坡底水平位移為9．3 cm，在考慮膨脹力后，位移場基本上分布在邊坡坡面附近，坡底水平位移達(dá)38 cm，表明邊坡已發(fā)生較大的滑動變形。

圖5 等效塑性應(yīng)變Fig．5 Equivalent plastic strain

圖7 （a）為降雨過程中監(jiān)控點(diǎn)飽和度過程曲線，可知8 h左右（30 000 s）邊坡坡面會形成一層飽和狀態(tài)，此時降雨入滲率會逐步下降并維持穩(wěn)定。圖7（b）為降雨過程中監(jiān)控點(diǎn)（未考慮膨脹）的合位移過程曲線，可知在降雨過程中坡腳A點(diǎn)位移變化迅速，從剛開始的1 cm左右迅速變化到12 cm以上；位于坡面中心B點(diǎn)位移變化比坡腳小，但也達(dá)到了8 cm；位于坡頂?shù)腃點(diǎn)位移變化很小，說明邊坡的破壞過程是由坡腳逐步向坡頂延伸的漸進(jìn)過程，具有明顯的漸進(jìn)破壞特點(diǎn)。

根據(jù)變形分析可知：

（1）降雨對膨脹巖邊坡產(chǎn)生較大的影響，降雨過程中坡腳范圍會先飽和并軟化，而且坡腳是應(yīng)力集中的部位，因此在該處首先出現(xiàn)較大的塑性應(yīng)變。

（2）隨著降雨的進(jìn)行，飽和度逐漸向邊坡內(nèi)部延伸，重度進(jìn)一步增加，下滑力逐漸增大，同時，其抗剪強(qiáng)度又因含水率增加而急劇下降，位移也因此非線性地增大；塑性應(yīng)變逐步沿著坡腳向上部擴(kuò)展開來，這是降雨引起邊坡漸進(jìn)性破壞的主要原因。

（3）當(dāng)考慮膨脹特性后，等效塑性應(yīng)變范圍由坡頂?shù)狡碌仔纬韶灤┑膮^(qū)域，表明膨脹力會加劇邊坡下滑趨勢，是滑坡的重要影響因素。

圖6 降雨后位移Fig．6 Displacements after rainfall

5．3 穩(wěn)定性分析

圖8 （a）、（b）分別為初始狀態(tài)、降雨后的滑弧示意圖，可知初始狀態(tài)安全系數(shù)為1．45，處于穩(wěn)定狀態(tài)，而降雨后處于滑動狀態(tài)。

由穩(wěn)定性分析可知：

（1）只有考慮膨脹巖的非飽和強(qiáng)度基本特性，才有可能保證邊坡存在淺層裂隙而未降雨情況下仍處于穩(wěn)定狀態(tài)，否則邊坡在初始狀態(tài)下就會失去穩(wěn)定，不可能呈現(xiàn)為深層滑弧。

圖7 監(jiān)控點(diǎn)飽和度和合位移的過程曲線Fig．7 Saturation degree and displacement process curves atmonitoring points（without considering expansion force）

圖8 滑弧示意圖Fig．8 The sketch of the sliding curves

（2）考慮了膨脹巖的膨脹特性后，邊坡的安全系數(shù)會進(jìn)一步降低，加劇邊坡穩(wěn)定破壞。

（3）在降雨后，由于淺層存在裂隙，入滲較快，強(qiáng)度指標(biāo)急劇下降，同時淺層重度也在不斷增加，致使下滑力逐步增大的同時抗滑力逐步減小，從而引起安全系數(shù)急劇下降，穩(wěn)定由深層滑弧逐漸變化到淺層滑弧。如果沒有淺層裂隙的存在，雨水難于入滲，也不會因含水率增加而引起抗剪強(qiáng)度下降。因此大氣影響帶膨脹巖的風(fēng)化、干濕循環(huán)影響的淺層裂隙，是膨脹巖邊坡淺層滑動的根源。

6 結(jié) 語

（1）膨脹巖邊坡在雨水入滲的情況下，重度增加，下滑力增大，同時，其抗剪強(qiáng)度又因含水率增加而急劇下降，這是降雨引起邊坡漸進(jìn)性破壞的主要原因，而膨脹特性進(jìn)一步加劇了這種影響。

（2）大氣影響帶膨脹巖的風(fēng)化和干濕循環(huán)影響導(dǎo)致的淺層裂隙，是膨脹巖邊坡淺層滑動的根源。

（3）本文的方法全面考慮了膨脹巖的基本特性與特殊性質(zhì)，可以完整地描述膨脹巖邊坡降雨引起滑坡發(fā)生的漸進(jìn)破壞過程。通過對膨脹巖現(xiàn)場試驗(yàn)的數(shù)值模擬，完整再現(xiàn)了膨脹巖邊坡降雨失穩(wěn)的整個過程，且材料參數(shù)均為室內(nèi)試驗(yàn)得到的成果，故分析成果具有相當(dāng)?shù)膶?shí)用價值。

［1］ 林魯生，蔣 剛．考慮降雨入滲影響的邊坡穩(wěn)定分析方法探討［J］．武漢大學(xué)學(xué)報(bào)（工學(xué)版），2001，34（1）：42－44．

［2］ 陳善雄，陳守義．考慮降雨的非飽和土邊坡穩(wěn)定性分析方法［J］．巖土力學(xué)，2001，22（4）：447－450．

［3］ 魏保義．用非飽和滲流理論研究降雨人滲的特性及對邊坡穩(wěn)定的影響［D］．南京：河海大學(xué)，2001．

［4］ 劉艷華．膨脹土邊坡的穩(wěn)定性研究［D］．武漢：長江科學(xué)院，1999．

［5］ 韋立德，楊春和，徐衛(wèi)亞，等．考慮飽和一非飽和滲流場和應(yīng)力場耦合的三維強(qiáng)度折減有限元程序研制［J］．水文地質(zhì)工程地質(zhì)，2006，33（3）：16－20．

［6］ Hibbitt，Karlson＆Sorensen，Inc．ABAQUSManualsversion 6．3［M］．USA：Hibbitt，Karlson＆Sorensen，Inc．，2003．

［7］ 陳善雄，陳守義．膨脹土判別與分類方法探討［J］．巖土力學(xué)，2005，26（12）：1895－1900．