李海亮，黃潤秋，吳禮舟，方 正

（1.地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610059；2.四川省交通運(yùn)輸廳交通勘察設(shè)計(jì)研究院，四川 成都 610059）

滑坡是一種在世界范圍內(nèi)普遍存在的地質(zhì)災(zāi)害現(xiàn)象，降雨是其失穩(wěn)的主要誘發(fā)因素［1］。降雨入滲坡體是一種典型的非飽和流固耦合現(xiàn)象［2］。要科學(xué)合理的分析評價邊坡降雨入滲過程和穩(wěn)定性，非飽和土流-固耦合是首先需要解決的重大問題。

在地質(zhì)、環(huán)境因素等的作用下，土體大多呈現(xiàn)層狀結(jié)構(gòu)，形成層狀土的原因有很多，但是都可將層狀土分成2大類［3］:1)滲透系數(shù)較小的細(xì)質(zhì)土覆蓋滲透系數(shù)較大的粗質(zhì)土；2)滲透系數(shù)較大的粗質(zhì)土覆蓋滲透系數(shù)較小的細(xì)質(zhì)土。眾多學(xué)者對降雨入滲和非飽和土流固耦合的問題進(jìn)行了研究，Thomas(1997年)提出一考慮非飽和土變形的熱、氣、水運(yùn)動耦合理論［4］；Wong(1998年)進(jìn)行了非飽和土的固結(jié)耦合數(shù)值計(jì)算研究［5］；朱偉等［6］研究設(shè)計(jì)室內(nèi)降雨入滲土柱試驗(yàn)，在初步揭示降雨入滲過程和規(guī)律的基礎(chǔ)上，探討準(zhǔn)確反映降雨入滲量的有限元計(jì)算方法；詹良通等［7］建立了降雨入滲條件下無限長斜坡內(nèi)水分運(yùn)移模型，在合理簡化的基礎(chǔ)上得到了解析解；徐晗等［8］建立一個考慮水力滲透系數(shù)特征曲線、土-水特征曲線以及修正的Mohr-Coulomb破壞準(zhǔn)則的非飽和土流固耦合有限元計(jì)算模型。同時，田東方等［9］、張曉詠等［10］、徐炎兵等［11］、吳禮舟等［12～13］在流固耦合數(shù)值分析和邊坡穩(wěn)定性問題方面也進(jìn)行了研究。然而，大多數(shù)研究者假設(shè)邊坡為均質(zhì)的，對非均質(zhì)邊坡研究較少；或者采用了流固準(zhǔn)耦合理論，并未真正實(shí)現(xiàn)流固的真耦合，這與實(shí)際情況是截然不同的，這可能導(dǎo)致研究的結(jié)果與實(shí)際情況產(chǎn)生偏差。在降雨條件下非均質(zhì)邊坡的入滲特性及變形等性狀不同于均質(zhì)邊坡。本文擬在前人研究的基礎(chǔ)上，通過數(shù)值分析方法，采用流固真耦合理論，以兩種典型的雙層層狀土坡為例，分析非均質(zhì)邊坡在降雨條件下的入滲過程及穩(wěn)定性，揭示了降雨誘發(fā)非均質(zhì)土坡的變形破壞機(jī)制，并對一些影響因素進(jìn)行分析，以期為研究多層土的降雨入滲以及更加準(zhǔn)確地評價邊坡穩(wěn)定性提供借鑒。

1 流固耦合理論及強(qiáng)度折減法

1.1 非飽和土流-固耦合理論

土體滲流場和應(yīng)力場之間是相互作用、相互影響的，其實(shí)質(zhì)為應(yīng)力場使土體發(fā)生變形，從而改變了介質(zhì)滲透系數(shù)；滲流場的變化又會引起孔隙水壓力分布的改變，使得滲流力和有效應(yīng)力發(fā)生改變。本文采用的大型有限元軟件ABAQUS可以對土中水的滲流和土體變形進(jìn)行直接耦合分析。

根據(jù)虛功原理［14］，固相介質(zhì)的應(yīng)力平衡為:

式中:δε——虛變形速率，δε=sym?σv/?()x；

σ'——有效應(yīng)力；

δv——虛應(yīng)變場；

t——單位面積的表面力；

f——單位體積的體積力(不含流體質(zhì)量)；

s、n——固相材料的飽和度、孔隙率；

ρf——流體密度；

g——重力加速度。

用拉格朗日公式將虛功方程離散化得到固相材料的有限元網(wǎng)格，此時流體還需要滿足連續(xù)方程:

式中:vf——滲流速度；

n——S面外法線方向。

ABAQUS能將滲流場和應(yīng)力場直接耦合，同時求解應(yīng)力平衡與滲流連續(xù)性方程，并采用Newton迭代法求解方程。

1.2 有限元強(qiáng)度折減法的原理

1975年Zienkiewicz等首次提出了抗剪強(qiáng)度折減的概念。基于強(qiáng)度儲備的概念可將穩(wěn)定安全系數(shù)Fs定義為當(dāng)土體材料的抗剪強(qiáng)度參數(shù)c和φ分別用其臨界抗剪強(qiáng)度參數(shù)cf和φf所代替后，結(jié)構(gòu)將處于臨界破壞狀態(tài)，其中:

強(qiáng)度折減法的基本思想是在外荷載不變的情況下，將土體強(qiáng)度參數(shù)按式(3)進(jìn)行折減，得到新的土體強(qiáng)度參數(shù)cf、φf進(jìn)行分析，通過不斷增加折減系數(shù)，直至土坡達(dá)到極限狀態(tài)，此時的折減系數(shù)就是所要求解的邊坡穩(wěn)定性系數(shù)。

1.3 有限元強(qiáng)度折減法的失穩(wěn)判據(jù)

彈塑性有限元強(qiáng)度折減法用于邊坡失穩(wěn)的判據(jù)大致可以分為三類:(1)有限元數(shù)值計(jì)算迭代的不收斂性；(2)特征部位(如坡頂、邊坡內(nèi)位移最大處等)位移(包括水平位移、豎向位移或總位移)發(fā)生突變；(3)廣義塑性應(yīng)變或等效塑性應(yīng)變從坡腳到坡頂貫通。目前對每種失穩(wěn)判據(jù)的有效性存在著不同的看法和解釋。鑒于此，本文將根據(jù)特征部位位移發(fā)生突變并結(jié)合塑性區(qū)貫通作為失穩(wěn)判據(jù)。

2 實(shí)驗(yàn)概況

為了深入討論和研究非均質(zhì)土坡在降雨條件下的入滲過程、規(guī)律和穩(wěn)定性變化，本文建立了兩種典型的非均質(zhì)土質(zhì)邊坡模型，即粗-細(xì)型(上層為#2砂土，下層為#1砂土)和細(xì)-粗型(上層為#2砂土，下層為#3砂土)。為了進(jìn)行對比分析，兩種模型中上層土均為1m厚的#2砂土，下層分別為#1和#3砂土(厚度為1.5m)。土的物理力學(xué)參數(shù)如表1所示，邊坡數(shù)值分析模型見圖1，邊坡坡比為1:1.5，AB、BC和CD段為降雨邊界；DE段為排水邊界，即當(dāng)孔壓為正時，假定流速與孔壓成正比，孔壓為負(fù)時，假定流速為零；AG和DE段固定水平方向的位移，EG段固定水平和豎直方向位移。L1、L2、L3和 L4為孔隙水壓力監(jiān)測剖面線。降雨強(qiáng)度考慮了20mm/h和10mm/h。

表1 土體的物理力學(xué)參數(shù)Table 1 Soil physical and mechanical properties

圖1 邊坡數(shù)值分析模型Fig.1 Numerical model of slope

滲透系數(shù)是非飽和土的重要參數(shù)之一，但是其變化范圍很大且不易測量，而從土-水特征曲線能獲取滲透系數(shù)，土-水特征曲線是非飽和土中體積含水量(飽和度)和吸力之間的關(guān)系曲線，非飽和土的行為特征與土-水特征曲線有密切關(guān)系。目前很多經(jīng)驗(yàn)公式可以描述土-水特征曲線，本文采用應(yīng)用很廣泛的van Genuchten模型。模型表達(dá)式為:

式中:θ——土的體積含水量；

ψ——壓力水頭(kPa)；

θs、θr——土的飽和、殘余體積含水量；

α、n——經(jīng)驗(yàn)擬合參數(shù)，m=1 －1/n。

非飽和土滲透系數(shù)表達(dá)式為:

式中:k——非飽和土滲透系數(shù)；

ks——飽和土滲透系數(shù)；

Se——有效飽和度，Se=1/[1 +(αψ)n]m。

本次實(shí)驗(yàn)砂土的土－水特征曲線如圖2所示，滲透系數(shù)如圖3所示。

圖2 土體的土－水特征曲線Fig.2 Soil-water characteristic curve of soils

圖3 土體的滲透系數(shù)曲線Fig.3 Hydraulic conductivity curve of soils

為了得到較為實(shí)際的邊坡初始條件，分析中假定整個坡體孔隙水壓力為－20kPa，在無流量邊界的條件下，運(yùn)用ABAQUS的Geostatic分析步進(jìn)行地應(yīng)力平衡，平衡后邊坡最大位移達(dá)到了10－8m，可以認(rèn)為邊坡達(dá)到了初始平衡狀態(tài)。獲取的邊坡初始狀態(tài)包括初始孔隙水壓力、初始應(yīng)力及飽和度等分布情況，可作為邊坡降雨入滲暫態(tài)分析的初始條件。計(jì)算得出的滲流場見圖4。

3 模擬結(jié)果與分析

圖4 初始時刻粗－細(xì)型邊坡滲流場孔壓等值線圖(單位:kPa)Fig.4 Coarse-fine type slope pore water pressure isoclines of initial seepage(kPa)

圖5為兩種類型邊坡在降雨強(qiáng)度10mm/h、歷時12h、耦合情況下的孔隙水壓力等值線圖，相應(yīng)的非耦合情況見圖6。等值線圖清晰地顯示了整個土坡孔隙水壓力的空間變化。耦合和非耦合情況下滲流場差異較大，主要原因是，考慮耦合作用時，應(yīng)力場會導(dǎo)致土體滲透系數(shù)變小，所以降雨時耦合較非耦合入滲將會慢很多。因此，考慮耦合作用更符合實(shí)際情況?；诖?，下面的分析將考慮非飽和滲流場和應(yīng)力場的耦合作用。

圖5 降雨強(qiáng)度10mm/h、歷時12h兩種類型邊坡耦合孔隙水壓力等值線圖(單位:kPa)Fig.5 Two type of slope isoclines of coupled pore water pressure with rainfall 12 hour and its intensity of 10 mm/h(kPa)

圖7 降雨強(qiáng)度10mm/h時，不同時刻L2剖面線上兩種類型邊坡孔隙水壓力－深度分布曲線Fig.7 Two type of slope pore water pressure-depth curve with rainfall intensity of 10 mm/h in profile L2

圖7顯示了降雨強(qiáng)度10mm/h時不同時刻L2剖面線上孔隙水壓力pw隨深度z的分布曲線。由等值線圖可以看出，表層孔壓變化很大，降雨約9h內(nèi)，由孔壓－深度分布曲線看出，下層土體孔壓基本未發(fā)生變化，這說明雨水還未入滲到下層土體，隨后，雨水入滲到下層土體，孔壓不斷增大，但是這種增大的趨勢對于兩種類型的非均質(zhì)土坡并不一樣。由此，降雨入滲過程可以分解為幾個階段。

對于細(xì)－粗型土坡，入滲過程可分解為2個過程:(1)雨水未進(jìn)入到土層接觸面，所有雨水均進(jìn)入到土體中，并且只在上層土中滲流；(2)雨水入滲到下層土，由于下層土滲透系數(shù)大于上層土，因此在土層接觸帶處不會產(chǎn)生積水情況。這時下層土體孔壓增大，沿深度呈近似線性分布，上層土體孔壓增加速率明顯減小。

對于粗－細(xì)型土坡，入滲過程將變得復(fù)雜一些，可以分解為3個過程:(1)這個過程與細(xì)－粗型土坡過程類似；(2)雨水進(jìn)入到下層土，由于土體滲透系數(shù)變小，在土層接觸面處將發(fā)生積水；(3)隨著降雨進(jìn)行，土層接觸面發(fā)生積水后，積水面將不斷向上遷移，最終將導(dǎo)致坡體表面積水。

3.1 降雨強(qiáng)度和時間對滑坡的影響

圖8顯示了降雨強(qiáng)度20mm/h時不同時刻L2剖面線上孔隙水壓力pw隨深度z的分布曲線。與圖7對比分析發(fā)現(xiàn)，隨著降雨強(qiáng)度的增大，雨水到達(dá)土層接觸面的時間約為6h，提前了3h。上下層土體在接觸面的孔壓差值大大增加，將分別達(dá)到18kPa和20kPa，這將容易造成更大的塑性變形。

圖8 降雨強(qiáng)度20mm/h時，不同時刻L2剖面線上孔隙水壓力－深度分布曲線Fig.8 Two type of slope pore water pressure-depth curve with rainfall intensity of 20 mm/h in profile L2

3.2 穩(wěn)定性分析

運(yùn)用強(qiáng)度折減法對不同降雨強(qiáng)度和歷時的邊坡穩(wěn)定性進(jìn)行了計(jì)算，以坡體位移突變并結(jié)合塑性區(qū)貫通作為失穩(wěn)判據(jù)。圖9為粗-細(xì)型邊坡降雨強(qiáng)度10mm/h的折減系數(shù)和水平位移的關(guān)系曲線，可以明顯的看出不同降雨歷時條件下，位移均存在突變點(diǎn)，位移突變的時間即為邊坡破壞時刻。降雨強(qiáng)度10mm/h、歷時12h的邊坡，當(dāng)折減系數(shù)為1.63時，等效塑性應(yīng)變等值云圖見圖10，此時塑性區(qū)已經(jīng)貫通到坡頂，已達(dá)到失穩(wěn)的臨界狀態(tài)。綜合分析位移突變曲線和塑性區(qū)貫通云圖，確定邊坡穩(wěn)定性系數(shù)為1.63。由圖10看出在兩種土層接觸面處產(chǎn)生了較大的塑性應(yīng)變，隨著變形不斷變大，沿著接觸面形成了滑動帶。

圖9 折減系數(shù)-水平位移分布曲線Fig.9 Reduction factor-horizontal displacement curve

圖10 等效塑性應(yīng)變等值云圖Fig.10 Nephogram of equivalent plastic strain

邊坡穩(wěn)定性系數(shù)隨降雨強(qiáng)度和降雨歷時的變化曲線見圖11。兩種典型非均質(zhì)土坡穩(wěn)定性系數(shù)均隨降雨強(qiáng)度和歷時增大而降低；降雨強(qiáng)度越大，穩(wěn)定性系數(shù)降低越快；隨后的降雨過程中，穩(wěn)定性降低速率受降雨強(qiáng)度影響減小，而且曲線近于平行。

圖11 兩種類型邊坡穩(wěn)定性系數(shù)隨時間變化曲線Fig.11 Two type of slope stability factortime curve

由圖11可見，無論降雨強(qiáng)度多大，曲線均存在一個拐點(diǎn)。降雨強(qiáng)度為10mm/h時，曲線拐點(diǎn)橫坐標(biāo)為約9h；降雨強(qiáng)度為20mm/h時，曲線拐點(diǎn)橫坐標(biāo)為約6h。由前述可知，雨強(qiáng)為10mm/h，經(jīng)過9h，雨水剛好進(jìn)入到下層土體；同樣，雨強(qiáng)為20mm/h，雨水入滲到下層土體的時間是6h。也就是說，雨水入滲到達(dá)下層土體時，邊坡穩(wěn)定性變化速率將減小。同時，粗-細(xì)型邊坡相對細(xì)-粗型邊坡受降雨影響較大，穩(wěn)定性降低較快；隨降雨時間增長，細(xì)-粗型邊坡穩(wěn)定性雖然在不斷降低，但是變化較小。

圖10顯示滑帶基本沿著兩土層的接觸面部分，為進(jìn)一步分析降雨誘發(fā)非均質(zhì)土坡變形失穩(wěn)機(jī)制，選取粗-細(xì)型邊坡，首先對抗剪強(qiáng)度參數(shù)c和φ折減1.63，然后進(jìn)行雨強(qiáng)10mm/h、歷時12h的降雨，降雨結(jié)束時，邊坡穩(wěn)定性系數(shù)為1。對圖1的4條剖面線孔隙水壓力進(jìn)行分析(圖12)。可見，降雨12h之后，孔隙水壓力變化很大，同時4條曲線均存在一個突變，突變位置分別為:L1、L2在 z=1.0m，L3在 z=1.5m，L4在z=2.0m。不難發(fā)現(xiàn)，這4個剖面中，孔隙水壓力突變位置均為土層接觸面位置，也是土坡滑帶所在位置。對于非均質(zhì)土坡，土層接觸面的位置容易產(chǎn)生孔隙水壓力的突變，造成有效應(yīng)力的迅速降低，使得土體產(chǎn)生較大的塑性變形，不斷累進(jìn)性變形之后便形成了滑帶。

圖12 不同剖面孔隙水壓力-深度分布曲線Fig.12 Pore water pressure-depth curve of different profiles

3.3 滲透系數(shù)對滑坡的影響

以圖1的粗-細(xì)型土質(zhì)邊坡模型為例，假設(shè)其他條件不變，將#1土體滲透系數(shù)取為2E-006m/s，以分析滲透性強(qiáng)弱對孔隙水壓力和邊坡穩(wěn)定性的影響，不同滲透系數(shù)時的穩(wěn)定性計(jì)算結(jié)果見圖13。

圖13 穩(wěn)定性系數(shù)隨時間變化曲線Fig.13 Slope stability factor-time curve

可見，在雨水未入滲到土層接觸面時，底層土體滲透系數(shù)對邊坡穩(wěn)定性影響很小。當(dāng)雨水到達(dá)接觸面之后，穩(wěn)定性系數(shù)降低速率隨著#1土體的滲透系數(shù)變小而增大，這是因?yàn)橄聦油馏w滲透系數(shù)越小，接觸帶部位越容易積水，進(jìn)而導(dǎo)致了穩(wěn)定性系數(shù)降低更大。

4 結(jié)論

(1)運(yùn)用流固真耦合理論分析降雨條件下的非均質(zhì)土坡的滲流變形更符合實(shí)際。細(xì)-粗型土坡的降雨入滲過程分為2個階段，粗-細(xì)型土質(zhì)邊坡的降雨入滲過程可以分為3個階段。

(2)當(dāng)雨水未入滲到底層土體時，穩(wěn)定性下降很快，降低速率基本不變；當(dāng)雨水入滲到下層土體時，邊坡穩(wěn)定性系數(shù)降低速率會變?。环€(wěn)定性系數(shù)和時間分布曲線存在一個拐點(diǎn)，拐點(diǎn)的橫坐標(biāo)即為雨水到達(dá)土坡土層接觸面時間。

(3)降雨誘發(fā)非均質(zhì)土坡變形破壞機(jī)制具有特殊性，在土層接觸帶處，產(chǎn)生了較大的孔隙水壓力突變，土體產(chǎn)生較大塑性變形，最終在接觸面處形成滑帶。

(4)雨水未到達(dá)土層接觸面時，下層土體的滲透系數(shù)對邊坡穩(wěn)定性影響很?。坏竭_(dá)之后，穩(wěn)定系數(shù)降低速率隨著下層土體的滲透系數(shù)變小而增大。

［1］王世梅.非飽和滑坡土體力學(xué)特性試驗(yàn)及其數(shù)值模擬［D］.武漢:武漢大學(xué)，2007.［WANG S M.Test and numerical simulation of the mechanical characteristics of the unsaturated Landslide［D］.Wuhan:WuHan University，2007］

［2］陳正漢，盧再有，朱元青.非飽和土的理論和實(shí)踐［J］.力學(xué)與實(shí)踐，2001，23(5):8-15.［CHEN Z H，LU Z Y， ZHU Y Q. Theory and practiceon unsaturated soils［J］.Mechanics and Engineering，2001，23(5):8 －15.(in Chinese)］

［3］韓同春，黃福明.雙層結(jié)構(gòu)土質(zhì)邊坡降雨入滲過程及穩(wěn)定性分析［J］.浙江大學(xué)學(xué)報:工學(xué)版，2012，46(1):39 －45.［HAN T C，HUANG F M.Rainfall infiltration process and stability analysis of two-layerd slope［J］.Journal of Zhejiang University:Engineering Science，2012，46(1):39 －45.(in Chinese)］

［4］Thomas H R，He Y.Coupled heat-moisture transfer theory for deformable unsaturated soil and its algorithmic implementation［J］.International Journal for Numerical Methods in Engineering，1997，40(18):3421－3441.

［5］Wong T T，F(xiàn)redlund D G，Krahn J.Numerical study of coupled consolidation in unsaturated soils［J］.Canadian Geotechnical Journal，1998，35(6):926 －937.

［6］朱偉，陳學(xué)東，鐘小春.降雨入滲規(guī)律的實(shí)測與分析［J］.巖土力學(xué)，2006，27(11):1873 － 1879.［ZHU W，CHEN X D，ZHONG X C.Observation and analysis of rainfall infiltration［J］.Rock and Soil Mechanics，2006，27(11):1873 － 1879. (in Chinese)］

［7］詹良通，賈官偉，陳云敏.考慮土體非飽和特性的無限長斜坡降雨入滲解析解［J］.巖土工程學(xué)報，2010，32(8):1214 －1220.［ZHAN L T，JIA G W，CHEN Y Mc.Analytical solution for rainfall infiltration into infinite long slopes considering properties of unsaturated soil［J］.Chinese Journal of Geotechnical Engineering，2010，32(8):1214 －1220.(in Chinese)］

［8］徐晗，朱以文，蔡元奇，等.降雨入滲條件下非飽和土邊坡穩(wěn)定分析［J］.巖土力學(xué)，2005，26(12):1957－1962.［XU H，ZHU Y W，CAI Y Q，et al.Stability analysis of unsaturated soil slopes under rainfall infiltration［J］.Rock and Soil Mechanics，2006，27(11):1873－1879.(in Chinese)］

［9］田東方，劉德富，王世梅，等.土質(zhì)邊坡非飽和滲流場與應(yīng)力場耦合數(shù)值分析［J］.巖土力學(xué)，2009，30(3):810－814.［TIAN D F，LIU D F，WANG S M，et al.Coupling numerical analysis of unsaturated seepage and stress fields for soil slope［J］.Rock and Soil Mechanics，2009，30(3):810 － 814.(in Chinese)］

［10］張曉詠，戴自航.應(yīng)用ABAQUS程序進(jìn)行滲流作用下邊坡穩(wěn)定分析［J］.巖石力學(xué)與工程學(xué)報，2010，29(1):2927－2934.［ZHANG X Y，DAI Z H.Analysis of slope stability under seepage by using ABAQUS program［J］.Chinese Journal of Rock Mechanics and Engineering，2010，29(1):2927 －2934.(in Chinese)］

［11］徐炎兵，韋昌富，李幻，等.非飽和土滲流與變形耦合問題的有限元分析［J］.巖土力學(xué)，2009，30(5):1490 － 1496.［XU Y B，WEI C F，LI H，et al.Finite element analysis of coupling seepage deformation in unsaturated soils［J］.Rock and Soil Mechanics，2009，30(5):1490 － 1496. (in Chinese)］

［12］Wu L.Z.，＆ Zhang L.M.Analytical solution to 1D coupled water infiltration and deformation in unsaturated soils.International Journal for Numerical and Analytical Methods in Geomechanics，2009，33(6):773－790.

［13］吳禮舟，黃潤秋.非飽和土滲流－變形耦合的數(shù)值分析［J］.土木建筑與環(huán)境工程，2011，33(3):63－67.［WU L Z，HUANG R Q.Numerical analysis of seepage and deformation in unsaturated soil［J］.Journal of Civil，Architectural ＆ Environmental Engineering，2011，33(3):63 －67.(in Chinese)］

［14］朱以文，蔡元奇，徐晗.ABAQUS與巖土工程分析［M］.香港:中國圖書出版社，2005.［ZHU Y W，CAI Y Q， XU H. ABAQUS and geotechnical engineering analysis［M］.Hongkong:China Tushu Publishing House，2005.(in Chinese)］