降雨條件下非飽和土邊坡滲流-應(yīng)力耦合分析

呂雨樺, 梁德賢, 王 瑩, 黃 翔

(桂林理工大學(xué) 地球科學(xué)學(xué)院, 廣西 桂林 541006)

自然環(huán)境中邊坡土體大多數(shù)處于非飽和狀態(tài), 非飽和土由于基質(zhì)吸力的存在， 其抗剪強(qiáng)度較高, 使邊坡處于穩(wěn)定狀態(tài)。受降雨入滲影響, 非飽和土的體積含水量增加, 力學(xué)強(qiáng)度顯著降低, 進(jìn)而導(dǎo)致滑坡的發(fā)生， 如千將坪滑坡[1]、 東鄉(xiāng)縣何坊村滑坡、 陳家灣溝滑坡[2]、 水城縣滑坡等, 均揭示了降雨入滲是導(dǎo)致非飽和土邊坡失穩(wěn)的主要因素之一。

針對(duì)這一問題, 以往是先分析降雨入滲引起坡體內(nèi)部滲流場(chǎng)的變化, 再結(jié)合極限平衡法或有限單元法來研究邊坡穩(wěn)定性[3-6], 忽略了土-水之間的相互影響。降雨入滲邊坡巖土體的實(shí)質(zhì)是滲流作用下滲流場(chǎng)與應(yīng)力場(chǎng)耦合作用, 雙場(chǎng)的影響具有相互性、 可逆性。目前, 雙場(chǎng)耦合分析主要分為兩大類: 一是間接耦合法, 先將兩場(chǎng)分開計(jì)算, 通過兩場(chǎng)的交叉迭代達(dá)到耦合目的[7-8]; 二是直接耦合法, 以滲流場(chǎng)與應(yīng)力場(chǎng)為未知值建立數(shù)學(xué)模型, 求得解析解來達(dá)到完全耦合目的[9-12]。對(duì)于土體而言, 直接耦合法通常是基于Biot固結(jié)理論, 該方法無需進(jìn)行兩場(chǎng)的反復(fù)迭代, 只需有限單元法求解即可得到全部結(jié)果, 具有簡(jiǎn)單、 直觀的特點(diǎn)。同時(shí), 邊坡變形過程本質(zhì)上是其內(nèi)部滲流場(chǎng)和應(yīng)力場(chǎng)直接耦合的結(jié)果, 因此本文選用直接耦合法分析邊坡穩(wěn)定性問題更具實(shí)際意義。

為此, 本文以廣西興業(yè)縣洛陽(yáng)滑坡為研究對(duì)象, 基于非飽和土流固耦合理論及有限元極限平衡原理, 通過Geo-Studio有限元數(shù)值分析軟件建立非飽和土邊坡數(shù)值模型, 對(duì)非飽和土邊坡滲流場(chǎng)、 應(yīng)力場(chǎng)、 變形場(chǎng)進(jìn)行直接耦合分析, 并在此基礎(chǔ)上研究非飽和土邊坡穩(wěn)定性演化規(guī)律。

1 耦合原理及基本方程

降雨入滲邊坡是典型的非飽和流固耦合現(xiàn)象, 其耦合原理: 降雨入滲使得滲流場(chǎng)以滲透力的形式改變?cè)械膽?yīng)力場(chǎng), 而應(yīng)力場(chǎng)的變化則通過體積應(yīng)變、 孔隙比的變化來改變土體滲透系數(shù), 從而影響土體的滲流場(chǎng)[13]。Geo-Studio軟件求解過程是以位移增量和孔隙水壓力增量為場(chǎng)變量, 依據(jù)初始條件、 邊界條件將平衡方程與滲流方程同時(shí)求解的過程。

1.1 非飽和土本構(gòu)方程

根據(jù)畢肖普有效應(yīng)力原理, 非飽和土應(yīng)變-應(yīng)力增量形式可表示為[14]

Δσ=DΔε-DmH(ua-uw)+Δua,

(1)

式中: Δσ為應(yīng)力增量; Δε為應(yīng)變?cè)隽?D為本構(gòu)矩陣; Δua為孔隙氣壓力矢量增量；H為與基質(zhì)吸力(ua-uw)有關(guān)的非飽和土模量(ua為孔隙氣壓力,uw為孔隙水壓力)。假設(shè)孔隙氣壓力恒等于大氣壓力, 上式可表示為

Δσ=DΔε+uwDmH。

(2)

1.2 有限元本構(gòu)方程

當(dāng)只有一個(gè)點(diǎn)荷載F施加于系統(tǒng)時(shí), 虛功方程為

(3)

式中:ε*為虛應(yīng)變;δ*為虛位移;σ為內(nèi)應(yīng)力。

將式(1)代入式(3)進(jìn)行數(shù)值積分, 則有限方程可寫為

∑BTDBΔδ+∑BTDmHNΔuw=∑F;

(4)

K=BTDB;Ld=BTDmHN,

式中: Δδ為位移矢量增量; Δuw為孔隙水壓力矢量增量;B為應(yīng)變矩陣;K為剛度矩陣;Ld為耦合矩陣;N為形函數(shù)行矢量。化簡(jiǎn)式(4)可寫成

KΔδ+LdΔuw=F。

(5)

1.3 滲流方程

基于達(dá)西定律, 通過一個(gè)單元體積土體的二維滲流方程為[13]

(6)

式中:kx、ky分別為x和y方向的滲透系數(shù);γw為單位水的容重;θw為體積含水量;t為時(shí)間。基于虛功原理, 滲流方程可以表示為孔隙水壓力和體積應(yīng)變的形式, 則得到有限元方程為

(7)

式中:Kf為單元?jiǎng)偠染仃?MN為質(zhì)量矩陣;Lf為滲流耦合矩陣;Q為邊界節(jié)點(diǎn)的滲流。聯(lián)立平衡方程(5)與滲流方程(7)可得有限元分析的耦合方程, 通過相應(yīng)的邊界條件即可對(duì)該耦合方程進(jìn)行求解。

2 工程概況

洛陽(yáng)滑坡位于廣西興業(yè)縣洛陽(yáng)鎮(zhèn)政府后山, 屬構(gòu)造侵蝕—?jiǎng)兾g丘陵地貌, 自然坡度22°～33°。如圖1所示, 滑坡平面形態(tài)為“簸箕形”, 主滑方向295°, 滑體平均坡度約35°, 軸線長(zhǎng)度150 m。前緣高程105 m, 后緣高程155 m, 相對(duì)高差為50 m?；潞缶壈l(fā)育多條弧形裂縫如圖2a所示, 主縫兩側(cè)垂向錯(cuò)距0.3～1.0 m; 坡中部橫向張性裂縫發(fā)育, 表面局部呈陡坎狀, 少量樹木歪斜(圖2b)。

圖1 滑坡平面圖

圖2 滑坡宏觀變形跡象

圖3 滑坡1—1′剖面圖

3 降雨條件下邊坡滲流-應(yīng)力耦合分析

選用Geo-Studio中的SEEP/W、 SIGMA/W、 SLOPE/W三種模塊相互結(jié)合的方式進(jìn)行滲流-應(yīng)力耦合分析。利用SEEP/W模塊從土-水特征曲線入手, 建立基質(zhì)吸力與滲透系數(shù)的關(guān)系, 計(jì)算初始滲流場(chǎng), 以文件形式保存, 作為SIGMA/W模塊初始條件; 選用SIGMA/W模塊耦合應(yīng)力/孔隙水壓力選項(xiàng), 進(jìn)行兩場(chǎng)耦合分析, 得出各時(shí)步的滲流場(chǎng)、 應(yīng)力場(chǎng)及位移場(chǎng)的計(jì)算結(jié)果; 將節(jié)點(diǎn)水頭信息、 應(yīng)力數(shù)據(jù)導(dǎo)入SLOPE/W模塊, 采用有限元法對(duì)邊坡穩(wěn)定性進(jìn)行動(dòng)態(tài)分析。

3.1 模型構(gòu)建及邊界條件

根據(jù)滑坡地質(zhì)原型和工程地質(zhì)條件, 本文選取1—1′地質(zhì)剖面建立數(shù)值模型。模型分為4層: ① 第四系殘坡積層, ② 全風(fēng)化花崗巖, ③ 強(qiáng)風(fēng)化花崗巖, ④ 中風(fēng)化花崗巖。有限元網(wǎng)格模式采用四邊形和三角形進(jìn)行劃分, 共計(jì)1 022個(gè)節(jié)點(diǎn)、 966個(gè)單元， 如圖4所示。初始滲流邊界條件: 模型兩側(cè)地下水位以上部分及底部邊界為零流量邊界, 左、 右側(cè)總水頭分別為94、 128 m; 浸潤(rùn)線主要沿基巖分布, 坡腳處趨于水平分布。瞬態(tài)滲流邊界條件: 模型坡面為單位流量邊界。位移邊界條件: 約束模型兩側(cè)的水平位移, 約束模型底部水平和豎向位移。為了更加方便分析降雨過程中邊坡失穩(wěn)的內(nèi)在規(guī)律, 選取特征點(diǎn)A～G進(jìn)行分析研究。

圖4 計(jì)算模型

3.2 參數(shù)選取及模擬方案

模型介質(zhì)為理想彈塑性材料, 服從莫爾-庫(kù)倫強(qiáng)度準(zhǔn)則, 物理力學(xué)性質(zhì)參數(shù)見表1。運(yùn)用SEEP/W模塊自帶的樣本函數(shù)擬合得出各巖土體的土水特征曲線圖5； 已知飽和滲透系數(shù)及土水特征曲線, 可依據(jù)Van Genuchten模型估算滲透系數(shù)曲線圖6。

圖5 土水特征曲線

圖6 滲透系數(shù)曲線

表1 巖土體物理力學(xué)參數(shù)取值

根據(jù)現(xiàn)場(chǎng)降雨?duì)顩r及當(dāng)?shù)亟涤陻?shù)據(jù), 模擬大暴雨(150 mm/d)工況下邊坡滲流場(chǎng)、 應(yīng)力場(chǎng)的變化。圖7a為降雨強(qiáng)度與時(shí)間的關(guān)系, 降雨強(qiáng)度先隨時(shí)間線性增加至最大值, 穩(wěn)定15 h后線性減小為0, 整個(gè)計(jì)算時(shí)長(zhǎng)為72 h, 其中， 降雨時(shí)長(zhǎng)為45 h, 累計(jì)降雨量為187.5 mm。當(dāng)降雨強(qiáng)度為6.25 mm/h, 設(shè)置坡腳水平面、 頂面降雨入滲6.25 mm/h, 坡面降雨入滲強(qiáng)度為6.25cos35°=5.2 mm/h。入滲強(qiáng)度隨降雨強(qiáng)度變化而變化, 詳見圖7b。

圖7 降雨強(qiáng)度與入滲強(qiáng)度隨時(shí)間變化

3.3 滲流場(chǎng)數(shù)值分析

圖8為模擬過程中邊坡孔隙水壓力分布演化?？傮w上, 浸潤(rùn)線以上的非飽和區(qū), 負(fù)孔隙水壓力隨著高程的增加而增大, 結(jié)合土-水特征曲線可知, 負(fù)孔隙水壓力與體積含水量呈反比關(guān)系, 距離浸潤(rùn)線越近負(fù)孔隙水壓力越小, 反之越大。隨著降雨入滲, 非飽和區(qū)域不斷減小, 飽和區(qū)域逐漸增大。在接近坡面位置處等值線密集, 負(fù)孔隙水壓力持續(xù)降低; 坡頂處等值線閉合, 由外到內(nèi)孔隙水壓力等值線呈現(xiàn)高→低→高的現(xiàn)象, 主要是因?yàn)檫吰卤韺油馏w最先對(duì)降雨作出響應(yīng), 而下部土體還未作出響應(yīng)所致。降雨未對(duì)深處浸潤(rùn)線造成影響。

圖9a為特征點(diǎn)A～G孔隙水壓力隨時(shí)間變化曲線。 0～15 h降雨強(qiáng)度持續(xù)增長(zhǎng), 坡體表層土體孔隙水壓力最先對(duì)降雨作出響應(yīng), 特征點(diǎn)A、B、D、F孔隙水壓力顯著增長(zhǎng), 其中點(diǎn)F增幅最大, 且接近飽和狀態(tài); 15～30 h降雨強(qiáng)度持穩(wěn), 點(diǎn)F孔隙水壓力回落, 點(diǎn)A、B、D孔隙水壓力線粘合, 增幅放緩, 表明此階段邊坡表層土體中的水分有向坡腳開始遷移的趨勢(shì); 30～72 h降雨強(qiáng)度線性減小至降雨停止后, 整體上, 孔隙水壓力出現(xiàn)回落減小態(tài)勢(shì), 減幅較小, 72 h(雨停后27 h)孔隙水壓力依然沒有完全消散的滯后現(xiàn)象, 如圖8d。坡體內(nèi)部特征點(diǎn)C、E、G孔隙水壓力基本不變, 與邊坡表層土體反應(yīng)迅速、 反應(yīng)劇烈形成鮮明對(duì)比。圖9b為特征點(diǎn)A～G體積含水量隨時(shí)間變化曲線, 與孔隙水壓力變化具有相似的規(guī)律。邊坡表層土體體積含水量與降雨強(qiáng)度呈正相關(guān)關(guān)系, 坡體內(nèi)部基本不受降雨影響。淺層土體的體積含水量持續(xù)增長(zhǎng)至30 h達(dá)到最大值, 30～72 h發(fā)生回彈現(xiàn)象。根據(jù)非飽和土力學(xué)理論可知, 非飽和土體基質(zhì)吸力對(duì)邊坡穩(wěn)定性產(chǎn)生積極的影響, 而隨著降雨的持續(xù)進(jìn)行, 非飽和區(qū)體積含水率增加, 導(dǎo)致負(fù)孔隙水壓力減小, 即基質(zhì)吸力降低, 基質(zhì)吸力對(duì)土體抗剪強(qiáng)度的貢獻(xiàn)隨之減少, 影響邊坡的穩(wěn)定性。自然環(huán)境中, 非飽和土邊坡基質(zhì)吸力在旱季升高, 雨季降低, 土體強(qiáng)度劣化, 致使滑坡災(zāi)害的發(fā)生。

圖8 非飽和土邊坡孔隙水壓力分布演化

圖9 特征點(diǎn)孔隙水壓力與體積含水量隨時(shí)間變化

3.4 應(yīng)力場(chǎng)數(shù)值分析

圖10為降雨前、 雨停后3 h最大剪應(yīng)力等值線。剪應(yīng)力分布整體上呈從坡面向坡內(nèi)逐漸增大的趨勢(shì)。降雨后, 邊坡表面土體最大剪應(yīng)力增大, 且坡腳處剪應(yīng)力集中區(qū)范圍增大, 特征點(diǎn)A的剪應(yīng)力由140 kPa增至160 kPa; 坡面中上部出現(xiàn)新剪應(yīng)力集中帶, 表層土體的剪應(yīng)力明顯提高至140 kPa。應(yīng)力集中區(qū)通常是滑坡發(fā)生剪切破壞最有利的部位, 可見坡中上部及坡腳處易發(fā)生剪切破壞。

圖10 非飽和土邊坡最大剪應(yīng)力等值線演化

滲流場(chǎng)引起孔隙水壓力變化, 破壞了其原始應(yīng)力平衡狀態(tài), 引起有效應(yīng)力變化, 如圖11所示。0～30 h降雨強(qiáng)度增大至定值, 坡面點(diǎn)B、D、F的有效應(yīng)力大幅度降低, 其中坡頂處點(diǎn)F降低幅度最大; 30～72 h降雨強(qiáng)度減小至雨停一段時(shí)間內(nèi), 坡面處有效應(yīng)力發(fā)生少量回彈。基于Fredlund非飽和土抗剪強(qiáng)度理論, 有效應(yīng)力減少對(duì)非飽和土抗剪強(qiáng)度產(chǎn)生消極影響。

圖11 特征點(diǎn)有效應(yīng)力隨時(shí)間變化

由于有效應(yīng)力降低, 邊坡表層應(yīng)力場(chǎng)改變, 產(chǎn)生變形響應(yīng), 主要體現(xiàn)為各土單元產(chǎn)生一定的位移, 如圖12所示。特征點(diǎn)水平、 豎向位移在降雨作用下均呈增加的態(tài)勢(shì), 淺層土體位移變化幅度較大。特征點(diǎn)F豎向位移變幅最大, 約是其水平位移的3倍, 主要因?yàn)槠马斖馏w吸水, 在豎直方向發(fā)生膨脹, 產(chǎn)生豎向變形; 特征點(diǎn)B的水平位移增長(zhǎng)速率最快, 豎向位移僅次于F點(diǎn), 這是由于坡體前緣較陡, 缺少臨空面方向的約束, 易產(chǎn)生水平位移, 導(dǎo)致小規(guī)模的土體滑塌; 特征點(diǎn)A僅發(fā)生水平位移, 降雨導(dǎo)致邊坡土體的自重增大, 繼而滑動(dòng)力增大, 該滑動(dòng)力傳遞于邊坡下部對(duì)坡腳處產(chǎn)生擠壓作用, 向邊坡臨空面產(chǎn)生水平變形, 可見坡腳處在降雨作用下易出現(xiàn)滑動(dòng)失穩(wěn)現(xiàn)象。

圖12 特征點(diǎn)位移隨時(shí)間變化

圖13為72 h位移等值線, 位移分布由坡內(nèi)向坡外逐漸增大。降雨大量入滲主要發(fā)生在淺層土體, 體積含水量增大, 淺層土體基質(zhì)吸力、 有效應(yīng)力降低, 應(yīng)力場(chǎng)的變化引起淺層土體相對(duì)深層土體更大的變形。位移量持續(xù)變化使得坡體處于持續(xù)運(yùn)動(dòng)的狀態(tài), 當(dāng)位移達(dá)到一定值, 坡體發(fā)生嚴(yán)重的變形破壞。前緣土體缺少水平方向的約束力最先破壞, 之后臨空面上部土體因失去下部土體的支撐而導(dǎo)致滑塌失穩(wěn), 出現(xiàn)更大規(guī)模的滑坡災(zāi)害。結(jié)合現(xiàn)場(chǎng)調(diào)查, 滑坡后緣處有多條拉張裂縫, 坡面局部呈陡坎狀, 判斷此次滑坡為牽引式滑坡。

圖13 72 h X-Y位移等值線

3.5 邊坡穩(wěn)定性演化分析

有限元極限平衡法是將邊坡有限元應(yīng)力分析結(jié)果與極限平衡法相結(jié)合進(jìn)行邊坡穩(wěn)定性分析的方法。該方法避免了傳統(tǒng)極限平衡法中的許多假設(shè), 既能考慮邊坡巖土體變形對(duì)穩(wěn)定性的影響, 又能給出明確的穩(wěn)定系數(shù)及破壞面。本文選取該方法對(duì)邊坡進(jìn)行穩(wěn)定性演化分析, 其原理是以有限元計(jì)算的應(yīng)力分布結(jié)果為基礎(chǔ), 通過應(yīng)力張量變換, 求出指定滑動(dòng)面上的應(yīng)力分布, 結(jié)合極限平衡法計(jì)算公式, 求解出該滑動(dòng)面的穩(wěn)定系數(shù)[15]。

選取SLOPE/W模塊SIGMA/W應(yīng)力分析類型, SIGMA/W模塊計(jì)算出的應(yīng)力值、 孔隙水壓力值為模擬基礎(chǔ), 對(duì)滑坡穩(wěn)定系數(shù)變化及危險(xiǎn)滑移面的位置進(jìn)行探討。邊坡穩(wěn)定系數(shù)隨時(shí)間變化的關(guān)系曲線如圖14所示, 降雨期間(0～45 h), 穩(wěn)定系數(shù)隨降雨強(qiáng)度增大而緩慢降低; 降雨強(qiáng)度維持最大值至減小為零的階段(15～45 h), 穩(wěn)定系數(shù)呈迅速衰減的態(tài)勢(shì); 降雨停止后(45～72 h), 穩(wěn)定系數(shù)繼續(xù)下降一定幅度再回升, 源于降雨入滲具有滯后性, 隨著時(shí)間的推移, 孔隙水壓力逐漸消散, 基質(zhì)吸力上升增強(qiáng)土體的抗剪強(qiáng)度, 致使穩(wěn)定系數(shù)回彈。邊坡穩(wěn)定系數(shù)最小值為0.997出現(xiàn)于54 h, 此時(shí)危險(xiǎn)滑移面的位置如圖15。結(jié)合圖3所示實(shí)際滑移面的位置, 與通過數(shù)值模擬計(jì)算所得危險(xiǎn)滑移面位置非常接近, 且與現(xiàn)場(chǎng)調(diào)查中滑坡后緣下錯(cuò)位置、 前緣剪出口位置基本一致, 符合邊坡在大暴雨工況下發(fā)生失穩(wěn)的實(shí)際情況。因此, 本文運(yùn)用Geo-Studio有限元數(shù)值分析軟件對(duì)滑坡進(jìn)行數(shù)值模擬， 并基于此分析其穩(wěn)定性是可靠的。

圖14 邊坡穩(wěn)定系數(shù)與時(shí)間的關(guān)系

圖15 54 h危險(xiǎn)滑移面位置圖

4 結(jié) 論

(1)降雨入滲過程中邊坡的滲流場(chǎng)、 應(yīng)力場(chǎng)、 位移場(chǎng)變化規(guī)律: 孔隙水壓力與體積含水量受降雨入滲的影響在淺表層土體發(fā)生劇烈變化, 與降雨強(qiáng)度呈正相關(guān)關(guān)系, 且降雨停止后的一段時(shí)間內(nèi)仍能保持, 表現(xiàn)出一定的滯后性; 邊坡表面土體最大剪應(yīng)力增大, 坡腳處剪應(yīng)力集中區(qū)范圍擴(kuò)大; 水平、 豎向位移在降雨作用下均呈逐漸增加的態(tài)勢(shì), 淺層土體位移變化幅度較大, 坡頂處豎向位移最大, 坡腳處水平位移為主。

(2)實(shí)際中, 大多數(shù)的邊坡處于非飽和狀態(tài), 降雨入滲過程本質(zhì)是邊坡從上到下由非飽和狀態(tài)轉(zhuǎn)變?yōu)榫植匡柡蜖顟B(tài)的漸變過程, 對(duì)于非飽和區(qū), 體積含水量的增加造成基質(zhì)吸力、 有效應(yīng)力下降, 弱化土體抗剪強(qiáng)度, 對(duì)邊坡的穩(wěn)定性產(chǎn)生不利影響。

(3)采用有限元極限平衡法研究降雨作用下邊坡的穩(wěn)定性, 穩(wěn)定系數(shù)隨著降雨持時(shí)的增大而減小, 最小穩(wěn)定系數(shù)出現(xiàn)在雨停后一段時(shí)間, 計(jì)算所得滑移面位置和實(shí)際情況接近, 說明本文采用Geo-Studio有限元數(shù)值分析軟件對(duì)滑坡進(jìn)行數(shù)值模擬研究結(jié)論可靠, 可為非飽和土邊坡的變形與穩(wěn)定分析提供一定的參考與借鑒。