降雨入滲條件下邊坡濕陷變形的數(shù)值模擬

王柳江，劉斯宏，李 卓，白福清

(河海大學(xué)水利水電學(xué)院，南京 210098)

降雨入滲是引起邊坡產(chǎn)生滑動破壞的主要因素之一，且非飽和土邊坡的變形和穩(wěn)定也一直是巖土工程界研究的熱點和難點之一．由于我國分布有廣泛的土質(zhì)邊坡，每年降雨引起的滑坡都對國家造成了重大的經(jīng)濟損失和人員傷亡．因此，深入研究降雨引起斜坡失穩(wěn)的規(guī)律并建立定量的分析模型對于滑坡的預(yù)防有重要的指導(dǎo)意義．

目前，研究降雨入滲對邊坡穩(wěn)定性影響的手段主要有 2類：①采用原位監(jiān)測，通過采集大量的監(jiān)測數(shù)據(jù)建立相應(yīng)的統(tǒng)計模型；②建立系統(tǒng)的數(shù)學(xué)物理計算模型對降雨入滲下的邊坡進(jìn)行定量分析．迄今為止，人們主要采用第 2種方法對邊坡穩(wěn)定性進(jìn)行分析．其中，通過滲流計算和極限強度平衡法相結(jié)合分析降雨入滲下邊坡穩(wěn)定性的方法被廣泛使用[1-4]，然而該法通常根據(jù)非飽和土的強度特性來反映邊坡的失穩(wěn)機制而更能直觀反映邊坡失穩(wěn)的變形規(guī)律卻無法得到．由于降雨入滲下的邊坡變形實質(zhì)上是典型的非飽和土固結(jié)問題，因此，它的解決必須基于非飽和土的滲流和應(yīng)力耦合理論[5-6]．Cho等[7]開發(fā)了水、氣、固三相耦合程序，且對降雨入滲下的均質(zhì)土坡進(jìn)行了分析，但沒有對邊坡變形進(jìn)行過深入研究．徐晗等[8]建立了一個考慮水力滲透系數(shù)特征曲線、土水特征曲線和修正 Mohr-Coulomb破壞準(zhǔn)則的非飽和土流固耦合計算模型，然而所采用的非飽和土本構(gòu)模型并不是真正意義上的彈塑性模型．因此，在非飽和土多孔介質(zhì)力學(xué)理論的基礎(chǔ)上采用非飽和土彈塑性模型對降雨入滲下邊坡的濕陷變形進(jìn)行分析具有重要意義．為了更準(zhǔn)確地反映非飽和土的變形特性，很多學(xué)者在彈塑性理論框架內(nèi)建立了非飽和土本構(gòu)模型，Alonso等[9-11]提出的 BBM 模型目前已被廣泛應(yīng)用.但該模型也存在一些缺陷：如通過 BBM 模型中的Load-Collapse(LC)屈服函數(shù)計算得到的濕化變形隨著圍壓的增大而增大，這與大量的試驗結(jié)果不符；模型采用凈應(yīng)力作為應(yīng)力狀態(tài)變量，無法很好地反映非飽和土的水力-力學(xué)耦合特性，且在數(shù)值計算中運用困難．為此，在 BBM 模型的基礎(chǔ)上對其修正進(jìn)行建模的思路已經(jīng)得到了很多學(xué)者的關(guān)注[12-13]．

筆者在非飽和土多孔介質(zhì)力學(xué)理論的基礎(chǔ)上，采用能夠考慮吸力影響以及非飽和土水力-力學(xué)耦合特性的修正BBM 模型，開發(fā)相應(yīng)的有限元程序，模擬了降雨入滲下邊坡的濕陷變形，研究了降雨持時、降雨強度以及飽和滲透系數(shù)對邊坡滲流和變形的影響．

1 非飽和土彈塑性本構(gòu)模型

BBM模型最早由西班牙學(xué)者Alonso提出，該模型能通過吸力變化來反映非飽和土干濕循環(huán)下的脹縮特性．但其僅是基于特定試驗結(jié)果的本構(gòu)模型，尚有一些地方不夠完善，例如：LC曲線在描述高圍壓下土體的濕化變形時與試驗結(jié)果不符；模型中參考應(yīng)力 pc的物理概念不明確，且當(dāng) pc等于飽和狀態(tài)下土體的前期固結(jié)應(yīng)力 p0*時，LC曲線在 p-s平面上為一條豎直線，此時無法反映吸力變化引起的土體塑性體應(yīng)變；模型采用凈應(yīng)力作為應(yīng)力狀態(tài)變量，在描述非飽和土的飽和-非飽和狀態(tài)過渡時存在困難．為了使其更好地適用于非飽和土的流固耦合分析，這里從應(yīng)力狀態(tài)變量、壓縮線斜率和LC屈服面3方面對其進(jìn)行了修正．

1.1 修正BBM模型簡介

修正后的BBM模型描述如下．

屈服函數(shù)采用了修正劍橋模型的形式，即

式中：p′為平均主應(yīng)力；q為廣義剪應(yīng)力；sp′為隨吸力增大的附加黏聚力；0p′為非飽和土的前期最大固結(jié)應(yīng)力，其值隨吸力的變化而變化，在 p-s平面上繪制出來即為LC屈服線．

與初始BBM模型不同，本文中的LC曲線方程是在高應(yīng)力狀態(tài)下推導(dǎo)得到的，當(dāng)圍壓增大到一定值時，所有吸力下的壓縮曲線交于一點，此時土體吸力降低不引起變形，則LC曲線函數(shù)為

式中：pn為不同吸力下土體壓縮曲線交點所對應(yīng)的應(yīng)力；p0*為飽和狀態(tài)下土體的前期最大固結(jié)應(yīng)力；λ( 0 )為飽和狀態(tài)下土體的壓縮系數(shù)；κ為土體回彈系數(shù)；s為吸力；pat為大氣壓；κs為圍壓等于 1,MPa時土體干化或濕化時由吸力變化引起的體積壓縮系數(shù)；λ( s )為非飽和狀態(tài)下吸力等于 s時對應(yīng)的土體壓縮曲線斜率．根據(jù)Sun等[14]試驗結(jié)果可知在高圍壓下，非飽和土體的壓縮曲線斜率隨吸力的增大而增大，其關(guān)系式為

式中λs為試驗參數(shù)，由吸力趨向于無窮大時的土體壓縮曲線斜率與吸力等于 0時的土體壓縮曲線斜率相減即可得到．圖1為根據(jù)式(2)和式(3)在p-s平面上繪制的 LC曲線，其中參數(shù)λ(0)=0.2，κ=0.03，pn= 1 MPa ，κs= 0 .002保持不變，變化參數(shù) p0*和λs得到．

為更好地描述非飽和土模型中水力-力學(xué)的相互作用關(guān)系，使模型更適用于有限元流固耦合分析，該模型采用了平均骨架應(yīng)力σi′j作為應(yīng)力狀態(tài)變量，即

式中：ijδ為Kronecker符號；rS為飽和度；au和wu分別為孔隙氣壓力和孔隙水壓力，通常巖土工程中認(rèn)為土體與大氣相通，則孔隙氣壓為大氣壓，ua=0．

圖1 不同參數(shù)下的LC曲線Fig.1 LC curves at different values of parameters p0* and λs

1.2 應(yīng)力-應(yīng)變關(guān)系式推導(dǎo)

有限元計算過程中需采用應(yīng)力-應(yīng)變的增量關(guān)系，因此進(jìn)行有限元計算前需將屈服函數(shù)先進(jìn)行變換．首先，根據(jù)相關(guān)聯(lián)流動法則得到塑性應(yīng)變增量為

式中Λ為比例系數(shù)，可根據(jù)屈服函數(shù)的連續(xù)性條件獲得．將式(1)寫成

其次，根據(jù)式(2)可得

由于塑性體應(yīng)變作為硬化參數(shù)，則屈服面上塑性體應(yīng)變增量相等，可采用飽和狀態(tài)下土體的應(yīng)力狀態(tài)及其應(yīng)力增量計算得到，即

然后，將式(5)代入式(8)，得

非飽和土的應(yīng)力-應(yīng)變增量關(guān)系可表示為

式中：De為彈性剛度矩陣；We為與吸力對應(yīng)的彈性剛度向量；dεe為彈性應(yīng)變增量，可表示為

將式(9)代入式(7)后再代入式(6)，同樣，將式(11)代入式(10)后再代入式(6)，即可得到比例系數(shù)Λ，即

將式(12)代入式(11)，然后聯(lián)合式(10)，移項后得到土體的彈塑性應(yīng)力-應(yīng)變增量表達(dá)式為

2 計算原理

式中：Δσi′j為平均骨架應(yīng)力張量的增量；xi為坐標(biāo)；Δuw為孔隙水壓力增量；Δbi為體力增量．

非飽和孔隙水連續(xù)方程為

2.1 非飽和土固結(jié)理論數(shù)學(xué)模型

在忽略溫度的影響以及假設(shè)土體與大氣相通的條件下，根據(jù)非飽和土多孔介質(zhì)理論，非飽和土固結(jié)方程由土骨架應(yīng)力平衡微分方程和非飽和土孔隙水連續(xù)流動方程組成．

土骨架的應(yīng)力平衡微分方程為

式中：n為孔隙率；Kw為水的體積壓縮模量；h為水頭；C ( s)為容水度，C =?θw?s，θw為體積含水量；εii為土骨架單元體應(yīng)變；kisj為飽和滲透系數(shù)張量；kr( s)為相對滲透系數(shù)；γw為水的容重．

2.2 土-水特征曲線模型

土體土-水特征曲線指基質(zhì)吸力和飽和度之間的關(guān)系，代表土體孔隙系統(tǒng)的持水特性．下文對文獻(xiàn)[15]中的邊坡進(jìn)行數(shù)值模擬，根據(jù)土-水特征實測值曲線擬合情況，其結(jié)果與文獻(xiàn)[7]介紹的土-水特征曲線模型一致，則本文采用的飽和度與吸力之間的關(guān)系為

式中：Sre為殘余水飽和度；Sm為最大飽和度；as、bs和 cs為參數(shù)．下文計算中土體水力滲透系數(shù)與吸力之間也滿足文獻(xiàn)[7]中的水力滲透特性曲線函數(shù)，即

式中：rk為相對滲透系數(shù)；aw、bw和 cw為擬合參數(shù)．圖 2(a)和(b)為文獻(xiàn)[15]中土-水特征及水力滲透系數(shù)的實測值采用式(18)和式(19)的擬合情況．

圖2 土-水特征曲線和水力滲透特性曲線Fig.2 Soil-water characteristic curve and hydraulic conductivity characteristic curve

3 邊坡降雨入滲變形的有限元模擬

為研究降雨入滲下非飽和土邊坡的變形特性，根據(jù)上文介紹的計算原理和應(yīng)力-應(yīng)變本構(gòu)關(guān)系，開發(fā)了有限元計算程序，對一降雨邊坡進(jìn)行了數(shù)值模擬[15]．

3.1 有限元計算模型及條件參數(shù)

圖3為原型邊坡斷面的有限元網(wǎng)格，圖中黑點代表監(jiān)測點，在計算結(jié)果分析中應(yīng)用．其中 a～e在坡面以下 5～10,m，而 f和 g分布較深，在坡面以下50～100,m．其邊界條件如下：底邊界為位移固定和不透水邊界；左右兩側(cè)為水平向位移約束及零流量邊界；邊坡表面為降雨入滲邊界，入滲量的確定將在下節(jié)介紹．初始條件包括初始應(yīng)力和初始孔壓，其中初始應(yīng)力根據(jù)邊坡自重確定，而初始孔壓的確定采用如下方法：設(shè)定右側(cè)水位為350,m 作為常水頭邊界，模擬未降雨條件下邊坡滲流至穩(wěn)定，將計算得到的孔壓和飽和度作為初始值，圖4為初始孔壓分布情況．根據(jù)文獻(xiàn)[15]論述，邊坡土體為礫類土，容重γ= 1 6.7 kN/m3，孔隙比 e = 0 .62，比重 Gs= 2 .51，飽和滲透系數(shù) ks= 4× 1 0-5m/s ，彈性模量 E = 1 0 MPa ，泊松比ν=0.3，黏聚力c′= 4 .0 kPa，內(nèi)摩擦角φ′=36.9°，其中土體的土-水特征曲線及水力滲透特性曲線如圖2(a)和(b)所示．由于文獻(xiàn)[15]并未針對修正BBM模型專門進(jìn)行參數(shù)測定，因此本文通過已知參數(shù)變換以及工程類比的方法近似獲取修正BBM模型參數(shù)．首先，根據(jù)彈性模量E與回彈系數(shù)κ之間的關(guān)系式為

式中p為圍壓，通常在100～150,kPa之間．將已知參數(shù)代入式(20)即可得到κ，且κ比λ小，一般κ= ( 0 .1 ～ 0 .2)λ，則λ可近似得到．

在p-q平面上臨界狀態(tài)線的斜率M可采用內(nèi)摩擦角表示，即

根據(jù)工程類比的方法，其余參數(shù)參考文獻(xiàn)[16]，則得到的修正 BBM 模型計算參數(shù)見表 1．本文共模擬了 30,h的降雨過程，平均降雨量為 1.2,m/d，在本文中將該情況設(shè)為基本工況．為研究降雨強度和土體飽和滲透系數(shù)對邊坡滲流和變形的影響，本文根據(jù)基本工況，通過改變降雨強度和飽和滲透系數(shù)，分別計算了飽和滲透系數(shù)保持 4.0×10-5,m/s不變，降雨強度為0.6,m/d、1.8,m/d、2.4,m/d的情況；以及降雨強度保持 1.2,m/d不變，飽和滲透系數(shù)為 1.0× 10-5,m/s、5.0×10-6,m/s、1.0×10-6,m/s的情況．

圖3 有限元計算網(wǎng)格及其監(jiān)測點

Fig.3 FEM mesh and selected nodes for result analysis

圖4 初始孔壓分布等值線(單位：kPa)Fig.4 Contour of initial pore-water pressure in slope(unit：kPa)

表1 修正BBM模型參數(shù)Tab.1 Parameters of improved BBM constitutive model

3.2 降雨入滲邊界的處理

降雨入滲水流在土體非飽和區(qū)是變化的，屬于水分通過巖土體包氣帶運移的兩相流問題，且雨水入滲量通常受降雨強度、降雨持時、巖土初始含水量以及入滲面幾何特征等因素的影響，因此降雨入滲是一個復(fù)雜的邊界非線性滲流問題．在滲流計算中，入滲面通常簡化為邊界條件，當(dāng)降雨強度大于入滲能力產(chǎn)生地表徑流時，為水頭邊界；當(dāng)降雨強度小于飽和滲透系數(shù)時，為流量邊界．本文采用了以積水點為判斷標(biāo)準(zhǔn)的降雨入滲計算方法[17]，其地表入滲能力通過地表水頭和入滲率進(jìn)行判斷．當(dāng)?shù)乇砉?jié)點非飽和時，以地表節(jié)點為臨界水頭(h＝0)時的入滲量為入滲能力來判斷是否由流量邊界轉(zhuǎn)化為水頭邊界；飽和時，以地表節(jié)點的入滲率為入滲能力來判斷是否由水頭邊界轉(zhuǎn)化為流量邊界．

4 計算結(jié)果及其分析

圖5 不同時刻水平位移分布等值線(單位：cm)Fig.5 Contours of horizontal displacement in slope at different time(unit：cm)

圖6 不同時刻豎向位移分布等值線(單位：cm)Fig.6 Contours of vertical displacement in slope at different time(unit：cm)

圖 5和圖 6為不同時刻邊坡水平位移和沉降等值線．由圖可見，水平位移主要發(fā)生在地下水位線附近及其以上的非飽和區(qū)，其中坡腳位置土體的水平位移向右，而地下水位線以上的向左，最大值出現(xiàn)在邊坡坡面拐角處，且邊坡順坡向的水平位移隨著降雨持續(xù)而增大；同時，降雨導(dǎo)致了邊坡內(nèi)出現(xiàn)了沉降和隆起，沉降主要發(fā)生在地下水位以上非飽和區(qū)，而地下水位以下飽和區(qū)的豎向位移朝上．其中，沉降主要是由非飽和區(qū)吸力減小導(dǎo)致土體強度降低軟化引起，而隆起則是由飽和區(qū)土體孔壓增大導(dǎo)致有效應(yīng)力減小而引起的土體回彈．由圖可知，最大濕陷變形發(fā)生在坡頂及邊坡坡面拐角處，且沉降值及隆起量均隨降雨持續(xù)而增大．因此，根據(jù)邊坡的變形分布可知坡面拐角和坡頂是決定該邊坡穩(wěn)定性的主要部位，由降雨引起的邊坡滑動面極有可能貫穿這2個部位．

圖7為a～g點的孔壓增量隨時間變化情況．由圖7可知：降雨初期，a～e的孔壓增量基本一致，隨著降雨持續(xù)，孔壓增速逐漸減小，且與位置有關(guān)，其中a、b和c的減小量依次增大，而d和e基本相同；邊坡淺層b和 d的孔壓增速明顯大于深層土體 f和g.由此可知，上述點的孔壓變化規(guī)律取決于土體所處的位置．對于邊坡表層土體，由于初始吸力隨高度增大而增大，在降雨條件下，吸力越低的表層土體能越快達(dá)到飽和狀態(tài)，且當(dāng)表層土體趨向于飽和時，由于雨水入滲率逐漸減小，則表面土體的孔壓增大速度也相應(yīng)減?。粚τ谏顚油馏w單元，由于降雨入滲是雨水由表及里流動的過程，因此，深層土體的孔壓增大與消散過程較表層土體有明顯的滯后性．圖8為監(jiān)測點塑性體應(yīng)變的變化規(guī)律，而降雨30,h后的塑性剪應(yīng)變等值線分布如圖9所示．由圖可見，表層土體單元的塑性體應(yīng)變隨高度的增大而減小，隨著降雨的持續(xù)而增大，塑性剪應(yīng)變最大值位于坡腳地下水溢出部位，說明塑性變形不僅與位置有關(guān)，還與降雨持續(xù)時間有關(guān)，且由降雨引起的邊坡塑性破壞區(qū)從坡腳往坡頂方向發(fā)展．圖10為 b～e在降雨過程中的應(yīng)力路徑方向．由圖10可見，其應(yīng)力路徑方向朝左上方，這主要與降雨的入滲有關(guān)．可以解釋如下：降雨入滲導(dǎo)致表層土體的孔壓增大而使平均應(yīng)力p′減小，同時，土體模量降低以及土體自重增大又使廣義剪應(yīng)力q′增大．

圖7 監(jiān)測點孔壓增量與時間關(guān)系Fig.7 Pore-water pressure increments of selected nodes versus time

圖8 監(jiān)測點塑性體應(yīng)變與時間關(guān)系Fig.8 Plastic volumetric strains of selected nodes versus time

圖9 降雨結(jié)束后塑性剪應(yīng)變分布等值線(單位：%)Fig.9 Contour of plastic shear strain at the end of rainfall(unit：%)

圖10 監(jiān)測點的應(yīng)力路徑Fig.10 Stress path of selected nodes

圖11～圖13為不同降雨強度下監(jiān)測點孔壓、沉降和水平位移隨降雨持時的變化．由于本文計算采用的降雨強度均小于土體的飽和滲透系數(shù)(4.0×10-5,m/s)，因此雨水基本能夠完全入滲．由圖可知，該邊坡內(nèi)孔壓、沉降和水平位移隨降雨強度的增大而增大，且當(dāng)降雨強度為 2.4,m/d時，降雨持續(xù) 30,h后，深層土體 g點出現(xiàn)了正孔壓，且該位置土體沉降出現(xiàn)了回彈，這與降雨引起的地下水位上升且有效應(yīng)力減小有關(guān)．

圖11 不同降雨強度下監(jiān)測點孔壓與時間關(guān)系Fig.11 Pore-water pressure versus time with different rainfall intensity

圖12 不同降雨強度下監(jiān)測點沉降與時間關(guān)系Fig.12 Settlement versus time with different rainfall intensity

圖13 不同降雨強度下監(jiān)測點水平位移與時間關(guān)系Fig.13 Horizontal displacement versus time with different rainfall intensity

圖14 不同飽和滲透系數(shù)時監(jiān)測點孔壓與時間關(guān)系Fig.14 Pore-water pressure versus time with different saturated permeability

圖 14～圖 16為不同飽和滲透系數(shù)時監(jiān)測點孔壓、沉降和水平位移隨降雨持時的變化．設(shè)降雨強度保持 1.2,m/d，改變土體飽和滲透系數(shù)．結(jié)果表明飽和滲透系數(shù)對邊坡內(nèi)孔壓和變形的影響較大．由圖14可見，對于坡面淺層土體，飽和滲透系數(shù)越小，吸力降低速度越快；而深層土體則是在飽和滲透系數(shù)接近降雨強度時吸力降低速度較快，而當(dāng)飽和滲透系數(shù)遠(yuǎn)大于或遠(yuǎn)小于降雨強度時，吸力減小速度較慢，該現(xiàn)象與降雨強度和土體滲透系數(shù)之間的關(guān)系有關(guān)．通常降雨初期邊坡表層土體的入滲量等于降雨量，所以根據(jù)滲流連續(xù)性條件以及達(dá)西定律可知：當(dāng)表層土體具有相同入滲量時，滲透系數(shù)越小，孔壓增量越大；對于深層土體，當(dāng)滲透系數(shù)遠(yuǎn)小于降雨強度時，由于入滲速度緩慢，導(dǎo)致孔壓增速明顯減慢，而當(dāng)飽和滲透系數(shù)遠(yuǎn)大于降雨強度時，由于土體透水性極強，孔壓消散速度同樣較快，則降雨引起的孔壓增量較?。蓤D15和圖16可見，當(dāng)土體的飽和滲透系數(shù)為1.0×10-5～5.0×10-6,m/s時，邊坡的變形遠(yuǎn)大于飽和滲透系數(shù)為 4.0×10-5和 1.0×10-6,m/s的變形，說明當(dāng)降雨強度大于且接近飽和滲透系數(shù)時，邊坡發(fā)生滑動破壞的可能性較大．通過不同飽和滲透系數(shù)下的耦合計算表明：當(dāng)邊坡的滲透系數(shù)遠(yuǎn)大于或遠(yuǎn)小于降雨強度時，邊坡失穩(wěn)的可能性較小，這與邊坡加固工程中設(shè)置坡內(nèi)排水系統(tǒng)或表面隔水措施的原理相似．

圖15 不同飽和滲透系數(shù)時監(jiān)測點沉降與時間關(guān)系Fig.15 Settlement versus time with different saturated permeability

圖16 不同飽和滲透系數(shù)時監(jiān)測點水平位移與時間關(guān)系Fig.16 Horizontal displacement versus time with different saturated permeability

5 結(jié) 論

(1) 降雨入滲時，最大水平位移位于坡面拐角處，而最大濕陷變形位于坡面拐角和坡頂位置，說明坡面拐角和坡頂是決定邊坡穩(wěn)定性的主要部位．則在降雨過程中，這2個部位極有可能率先發(fā)生塑性破壞而誘發(fā)局部滑動，因此，對該部位進(jìn)行防/排水或加固處理對邊坡的穩(wěn)定性具有重要意義．

(2) 降雨期間邊坡內(nèi)孔壓的變化與土體位置有關(guān)，邊坡底部孔壓增量較頂部小，深層較表層?。瑫r，在初始應(yīng)力狀態(tài)接近的情況下，土體的塑性變形與初始吸力有關(guān)，初始吸力越大，塑性變形越小，且降雨引起的塑性破壞區(qū)從坡腳往坡頂發(fā)展．

(3) 在飽和滲透系數(shù)大于降雨強度的條件下，邊坡內(nèi)孔壓和變形的發(fā)展與降雨持時和降雨強度成正比．

(4) 降雨條件下，邊坡內(nèi)滲流和變形的變化受土體飽和滲透系數(shù)的影響顯著，當(dāng)土體的飽和滲透系數(shù)小于且接近降雨強度時，邊坡表層土體位移較大；而當(dāng)飽和滲透系數(shù)遠(yuǎn)大于或小于降雨強度時，邊坡表層土體的位移較?。?/p>

(5) 通過對該邊坡實例的分析，說明了本文建立的計算模型以及開發(fā)的程序能夠客觀地反映非飽和土邊坡在降雨入滲下的濕陷變形特性，為降雨入滲下邊坡穩(wěn)定性的評價及其滑坡預(yù)測提供了新的手段．

［1］吳宏偉，陳守義，龐宇威. 雨水入滲對非飽和土坡穩(wěn)定性影響的參數(shù)研究[J]. 巖土力學(xué)，1999，20(1)：1-14.

Ng Wangwai Charles，Chen Shouyi，Pang Yuewai. Parametric study of effects of rain infiltration on unsaturated slopes[J].Rock and Soil Mechanics，1999，20(1)：1-14(in Chinese).

［2］Chen H，Chen R H，Yu F C，et al. The inspection of triggering mechanism for a hazardous mudflow in an urbanized territory[J].Environmental Geology，2004，45(7)：899-906.

［3］Ng C W W，Pang Y W. Influence of stress state on soilwater characteristics and slope stability[J].Journal of Geotechnical and Geoenvironmental Engineering，2000，126(2)：157-166.

［4］吳長富，朱向榮，尹小濤，等. 強降雨條件下土質(zhì)邊坡瞬態(tài)穩(wěn)定性分析[J]. 巖土力學(xué)，2008，29(2)：386-391.

Wu Changfu，Zhu Xiangrong，Yin Xiaotao，et al.Analysis of soil slope’s transient stability under intensive rainfall[J].Rock and Soil Mechanics，2008，29(2)：386-391(in Chinese).

［5］沈珠江. 非飽和土簡化固結(jié)理論及其應(yīng)用[J]. 水利水運工程學(xué)報，2003(4)：1-6.

Shen Zhujiang. Simplified consolidation theory for unsaturated soils and its application [J].Hydro Science and Engineering，2003(4)：1-6(in Chinese).

［6］陳鐵林，沈珠江，楊代泉. 濕陷性黃土渠道浸水變形的數(shù)值模擬[J]. 水利學(xué)報，2005，36(3)：309-313.

Chen Tielin，Shen Zhujiang，Yang Daiquan. Numerical simulation on wetting deformation of collapsible loess ditch[J].Journal of Hydraulic Engineering，2005，36(3)：309-313(in Chinese).

［7］Cho S E，Lee S R. Instability of unsaturated soil slopes due to infiltration[J].Computers and Geotechnics，2001，28(3)：185-208.

［8］徐 晗，朱以文，蔡元奇，等. 降雨入滲條件下非飽和土邊坡穩(wěn)定分析[J]. 巖土力學(xué)，2005，26(12)：1957-1962.

Xu Han，Zhu Yiwen，Cai Yuanqi，et al. Stability analysis of unsaturated soil slopes under rainfall infiltra tion[J].Rock and Soil Mechanics，2005，26(12)：1957-1962(in Chinese).

［9］Alonso E E，Gens A，Josa A. A constitutive model for partially saturated soils[J].Geotechnique，1990，40(3)：405-430.

［10］Bolzon G，Schrefler B A，Zienkiewicz O C. Elastoplastic soil constitutive laws generalized to partially saturated states[J].Geothecnique，1996，46(2)：279-289.

［11］Sheng Daichao，F(xiàn)redlund D G，Gens A. A new modeling approach for unsaturated soils using independent stress variables[J].Can Geotech J，2008，45：511-534.

［12］Wheeler S J，Sivakumar V. An elasto-plastic critical state framework for unsaturated soils[J].Geotechnique，1995，45(1)：35-53.

［13］繆林昌. 非飽和土的本構(gòu)模型研究[J]. 巖土力學(xué)，2007，28(5)：855-860.

Miao Linchang. Research of constitutive model of unsaturated soil[J].Rock and Soil Mechanics，2007，28(5)：855-860(in Chinese).

［14］Sun De’an，Sheng Daichao，Xu Yongfu. Collapse behaviour of unsaturated compacted soil with different initial densities[J].Can Geotech J，2007，44：673-686.［15］Chen R H，Chen H P，Chen K S，et al. Simulation of a slope failure induced by rainfall infiltration[J].Environmental Geology，2009，58(5)：943-952.

［16］Sun D A，Sheng D，Li X，et al. Elastoplastic modelling of hydraulic and stress-strain behaviour of unsaturated soils under undrained conditions[J].Computers and Geotechnics，2008，35(6)：845-852.

［17］朱 偉，程南軍，陳學(xué)東，等. 淺談非飽和滲流的幾個基本問題[J]. 巖土工程學(xué)報，2006，28(2)：235-240.

Zhu Wei，Cheng Nanjun，Chen Xuedong，et al. Some fundamental problems of unsaturated seepage[J].Chinese Journal of Geotechnical Engineering，2006，28(2)：235-240(in Chinese).