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

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

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

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

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

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

1 耦合原理及基本方程

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

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

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

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

(1)

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

Δσ=DΔε+uwDmH。

(2)

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

當只有一個點荷載F施加于系統(tǒng)時, 虛功方程為

(3)

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

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

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

(4)

K=BTDB;Ld=BTDmHN,

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

KΔδ+LdΔuw=F。

(5)

1.3 滲流方程

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

(6)

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

(7)

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

2 工程概況

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

圖1 滑坡平面圖

圖2 滑坡宏觀變形跡象

圖3 滑坡1—1′剖面圖

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

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

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

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

圖4 計算模型

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

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

圖5 土水特征曲線

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

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

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

圖7 降雨強度與入滲強度隨時間變化

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

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

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

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

圖9 特征點孔隙水壓力與體積含水量隨時間變化

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

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

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

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

圖11 特征點有效應(yīng)力隨時間變化

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

圖12 特征點位移隨時間變化

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

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

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

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

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

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

圖15 54 h危險滑移面位置圖

4 結(jié) 論

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

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

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