王淑云，葉天力，魯曉兵，聶雪媛

(中國科學(xué)院力學(xué)研究所流固耦合系統(tǒng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，北京 100190)

0 引言

降雨是誘發(fā)斜坡破壞的主要因素之一。降雨所致的泥石流和滑坡災(zāi)害，會嚴(yán)重威脅人民生命安全和帶來經(jīng)濟(jì)損失。近幾年來，我國每年都會有較大的泥石流發(fā)生，尤其在滇西北和東北、川西、甘南及白龍江流域等地區(qū)。典型的例子是在2008年9月25日，四川省北川縣城，因持續(xù)降雨引發(fā)多處山體滑坡和泥石流，30多人死亡或失蹤，老縣城一半以上被泥石流掩埋，最深處的污泥達(dá)40m以上。2010年9月9日，北川再遭泥石流沖擊，橋梁被沖毀，雷禹路全線癱瘓。2013年7月16日，北川的桂溪鄉(xiāng)陳家壩因暴雨突發(fā)泥石流，50余萬m3泥石流沖入青林溝，將河床整體抬高近4m，形成堰塞湖。

降雨誘發(fā)泥石流是雨水不斷入滲至地表土層內(nèi)部的結(jié)果。在降雨過程中，地表土層的飽和度增加，孔壓上升，基質(zhì)吸力減少，抗剪強(qiáng)度下降。當(dāng)降雨達(dá)到一定時(shí)間，雨水入滲到土層的深度增加，土體強(qiáng)度大大降低，土體失穩(wěn)，泥石流起動(dòng)。目前，在降雨引發(fā)滑坡和斜坡破壞研究方面，已有不少文獻(xiàn)報(bào)道［1-9］。但因?yàn)槊總€(gè)地區(qū)的降雨強(qiáng)度和歷時(shí)、現(xiàn)場地形和坡度以及土質(zhì)情況的復(fù)雜性，對其在降雨下的斜坡破壞的模擬和分析結(jié)果會有所不同，尤其反應(yīng)在定量方面。因此對特定地區(qū)土的斜坡破壞的分析和預(yù)測必須建立在其所處的外界環(huán)境和土本身的特性基礎(chǔ)上。

本文采用FLAC兩相流模塊，針對北川地區(qū)地形(主要是坡角)和土質(zhì)(現(xiàn)場土)，對北川地區(qū)土在不同降雨強(qiáng)度和降雨歷時(shí)情況下的斜坡破壞進(jìn)行了數(shù)值模擬，并對降雨過程中土體位移、速度矢量等變化規(guī)律進(jìn)行了分析。同時(shí)開展了室內(nèi)降雨模型實(shí)驗(yàn)，并將數(shù)值模擬計(jì)算結(jié)果與室內(nèi)實(shí)驗(yàn)結(jié)果進(jìn)行比較，驗(yàn)證了數(shù)值模擬計(jì)算和分析的可信度和可行性。

1 數(shù)值模擬和模型試驗(yàn)

1.1 數(shù)值模擬介紹

土體模型采用修正的Mohr-Columb模型，二維的飽和-非飽和土中水的滲流控制方程采用Richards方程。有限差分網(wǎng)格如圖1所示，其中X軸為沿土坡方向，Y軸為垂直于坡面方向。模型形狀為梯形體型。坡頂350mm，坡底520mm，坡面土層厚度100mm。在數(shù)值模擬中，土坡坡角分別取27°、30°、33°，本文分析以坡角27°時(shí)的數(shù)值計(jì)算結(jié)果為主。

圖1 有限差分網(wǎng)格Fig.1 Element mesh

北川土的主要土性參數(shù)為:干容重1.58g/cm3，粘聚力20kPa，體積模量800MPa，剪切模量600MPa，滲透系數(shù)1.0×10-5m/s。土體初始飽和度設(shè)為0.5，并假定殘余飽和度為0。降雨參數(shù)的選取基于北川當(dāng)?shù)亟涤昵闆r，分別采用5種雨強(qiáng):1×10-5m/s、2×10-5m/s、3 ×10-5m/s、4 ×10-5m/s和 5 ×10-5m/s。

1.2 室內(nèi)模型實(shí)驗(yàn)介紹

模型實(shí)驗(yàn)簡圖如圖2所示，土坡坡度為27°，均勻鋪土，土層長 ×寬 ×厚為1000mm×1200mm×100mm，土的顆分結(jié)果參見圖3。

圖2 模型實(shí)驗(yàn)示意圖Fig.2 Sketch of model test

圖3 北川土的顆分曲線Fig.3 Grain size distribution

實(shí)驗(yàn)選取的雨強(qiáng)分別為:1.75×10-5m/s、2.25 ×10-5m/s、3.25 × 10-5m/s、3.5 × 10-5m/s、4 ×10-5m/s、5 ×10-5m/s，記錄降雨歷時(shí)和斜坡破壞時(shí)間。

2 數(shù)值模擬和模型試驗(yàn)結(jié)果分析

2.1 數(shù)值模擬結(jié)果分析

2.1.1 土體水平位移分析

圖4為雨強(qiáng)1×10-5m/s時(shí)土體水平方向位移降雨歷時(shí)的變化情況。從圖4可以看出，小雨強(qiáng)情況下，降雨初期土體的水平位移相對較小，最大水平位移發(fā)生在坡腳處，這與室內(nèi)試驗(yàn)現(xiàn)象是比較一致，在試驗(yàn)中發(fā)現(xiàn)，在降雨初期，坡腳處的土體首先發(fā)生坍塌滑動(dòng)。數(shù)值模擬結(jié)果證實(shí)坡腳處的這種坍塌滑動(dòng)主要是由土體的水平方向的位移引起的。土體水平方向的位移主要發(fā)生在模型左側(cè)，如果取一沿坡面方向的截面，則自左至右土體水平位移呈逐漸減小的趨勢，最右側(cè)土體的水平位移很小甚至為零。如果在坡面左側(cè)部分取一豎直截面進(jìn)行分析，則土體的水平位移自上而下先增加后減小，在該截面中部的某一點(diǎn)土體水平位移取得最大值，且隨著該截面的向右移動(dòng)，最大值點(diǎn)的位置將不斷的向上移動(dòng)，這些最大值點(diǎn)構(gòu)成了一條線，這條線所處的位置即為潛在的滑移帶。

圖4 土體水平方向(沿坡面方向)位移分布圖Fig.4 Distribution of horizontal displacement

2.1.2 土體豎直位移分析

圖5為雨強(qiáng)1×10-5m/s時(shí)土體豎直方向位移隨降雨歷時(shí)的變化情況。從圖5可以看出，小雨強(qiáng)情況下，土體發(fā)生了不均勻沉降。土體的豎向位移隨降雨歷時(shí)的增加不斷增大，且土體上層的豎向位移遠(yuǎn)超過下層土體的豎向位移。土體豎向位移最大值出現(xiàn)在土體表層，且隨降雨歷時(shí)的增加，土體豎向位移的最大值位置逐漸向左側(cè)移動(dòng)，最終土體豎向位移的最大值點(diǎn)大約位于距坡頂三分之二處。結(jié)合前文對土體水平位移的分析，可以得出，土體之所以出現(xiàn)不均勻的沉降主要是由于土體發(fā)生了水平位移引起的。隨降雨歷時(shí)的增加，左側(cè)土體不斷的向左下方滑動(dòng)，土體總體積是一定的，這就引起土體的不均勻沉降，出現(xiàn)圖中的現(xiàn)象。正是由于土體豎向位移的不均勻性，隨豎向位移最值點(diǎn)的位移值不斷增大，該點(diǎn)土體將與上部土體產(chǎn)生錯(cuò)位，這為滑移帶的產(chǎn)生在表層形成了突破口。

圖5 土體豎向方向(垂直坡面方向)位移分布圖Fig.5 Distribution of vertical displacement

2.1.3 土體速度矢量分析

圖6為雨強(qiáng)1×10-5m/s情況下土體速度矢量隨降雨歷時(shí)的變化情況。從圖6可以看出，在小雨強(qiáng)情況下，隨降雨歷時(shí)的增加，土體的速度矢量不斷增大，且在降雨初期，土體速度增加的速率比較小，但隨著降雨的持續(xù)，增加的速率也越來越大。從圖7可以看到，土體速度矢量在水平方向分量的最大值位于坡腳處，這與前文分析中得到的土體的水平最大位移發(fā)生與坡腳處是一致的。土體速度矢量在豎直方向分量的最大值大約在坡體表層距坡頂大約三分之二處，這與前文的分析結(jié)果也是一致的。從整個(gè)坡體來看，速度矢量在豎直方向上的分量明顯小于其在水平方向上的分量，也就是說土體的水平運(yùn)動(dòng)處于主導(dǎo)地位。從土體產(chǎn)生較大運(yùn)動(dòng)的面積上來看，隨降雨歷時(shí)的增加，面積逐漸減小，且產(chǎn)生較大運(yùn)動(dòng)的區(qū)域不斷向左側(cè)偏移。

圖6 土體速度矢量與降雨歷時(shí)關(guān)系Fig.6 Distribution of velocity vector

圖7所示為土體速度矢量最大值隨降雨歷時(shí)變化曲線?？梢钥吹?，土體速度矢量的最大值也隨著降雨歷時(shí)的增加而增大，降雨歷時(shí)越長，最大速度矢量值增長的越快。

圖7 土體速度矢量最大值隨降雨歷時(shí)變化曲線Fig.7 Maximum velocity vector versus rainfall duration

2.2 斜坡破壞及與室內(nèi)試驗(yàn)結(jié)果的分析和比較

圖8、9為四種雨強(qiáng)下土體水平向最大位移以及豎向最大位移隨時(shí)間變化的關(guān)系曲線?？梢钥闯?相同雨強(qiáng)條件下，整個(gè)土體位移場的最大位移值隨時(shí)間不斷增大，在降雨初期土體位移最大值比較小，當(dāng)降雨達(dá)到一定時(shí)間土體位移的最大值將有一突變。隨雨強(qiáng)增大這種突變所需要的時(shí)間越短。這一現(xiàn)象也在坡腳30°和33°情況中得到了證實(shí)。

圖8 不同雨強(qiáng)下土體水平向最大位移變化曲線Fig.8 Maximum horizontal displacement under different rainfall density conditions

圖9 不同雨強(qiáng)下土體豎向最大位移變化曲線Fig.9 Maximum vertical displacement under different rainfall density conditions

對于坡角的影響，選取雨強(qiáng)為4×10-5m/s情況，針對坡腳27°、30°和33°進(jìn)行模擬計(jì)算。圖10、11 為不同坡角下土體水平向最大位移以及豎向最大位移隨時(shí)間的變化曲線。從圖中可以看出，相同雨強(qiáng)下，隨坡角的增大，土體最大水平位移和最大豎向位移發(fā)生突變所需時(shí)間大大縮短。

圖10 不同坡角下土體水平最大位移曲線Fig.10 Maximum horizontal displacement under different slope angle conditions

本文模擬主要采用坡面特征點(diǎn)的位移突變作為泥石流是否起動(dòng)的判據(jù)，因?yàn)樨Q向最大位移發(fā)生在距坡頂三分之二處，該處是土體滑移帶與坡面的交匯點(diǎn)，所以我們?nèi)∪鐖D1所示的A點(diǎn)作為特征點(diǎn)。

從圖12中可以看出，在坡角27°時(shí)，雨強(qiáng)越小，斜坡破壞所需時(shí)間越長;當(dāng)雨強(qiáng)從較大開始變小時(shí)，斜坡破壞所需時(shí)間急劇減小，隨著雨強(qiáng)的進(jìn)一步減小，斜坡破壞所需時(shí)間減小的速度逐漸變小，最終斜坡破壞時(shí)間趨于一個(gè)很大值。相同雨強(qiáng)下，坡角越大，斜坡破壞所需時(shí)間越短。這些規(guī)律跟室內(nèi)的人工降雨試驗(yàn)是一致的。從圖13可以看出，在坡角為27°時(shí)，A點(diǎn)水平速度隨降雨歷時(shí)的變化趨勢與圖12相同，降雨初期速度比較小，當(dāng)降雨達(dá)到一定時(shí)間速度也發(fā)生突變.從曲線還可以看出，速度突變時(shí)間超前于相應(yīng)的位移突變時(shí)間。

圖11 不同坡角下土體豎向最大位移曲線Fig.11 Maximum vertical displacement under different slope angle conditions

圖12 坡角27°時(shí)A點(diǎn)水平位移曲線Fig.12 Horizontal displacement of A when slope angle is 27°

圖13 坡角27°時(shí)A點(diǎn)水平速度曲線Fig.13 orizontal velocity of A when slope angle is 27°

以特征點(diǎn)位移是否發(fā)生突變作為泥石流起動(dòng)與否的判斷標(biāo)準(zhǔn)，從相同坡角不同雨強(qiáng)情況下的特征點(diǎn)的位移曲線，可以得到不同雨強(qiáng)下所對應(yīng)的斜坡破壞所需的降雨歷時(shí)，并將其與試驗(yàn)所得結(jié)果進(jìn)行對比，結(jié)果如圖14所示?？梢钥吹剑瑪?shù)值計(jì)算得到的斜坡破壞-降雨歷時(shí)曲線與試驗(yàn)得到的結(jié)果規(guī)律相同。即在坡角一定時(shí)，當(dāng)雨強(qiáng)較小，斜坡破壞所需時(shí)間較長，隨雨強(qiáng)的不斷增大，斜坡破壞所需時(shí)間越來越短，且所需時(shí)間減小的速率越來越大。

圖14 斜坡破壞時(shí)間隨雨強(qiáng)變化關(guān)系(室內(nèi)試驗(yàn)和數(shù)值模擬結(jié)果對比)Fig.14 Slope failure time versus rainfall duration

3 結(jié)論

初步結(jié)論如下:

(1)土體水平方向最大位移發(fā)生在坡腳處，且土體水平方向位移隨降雨歷時(shí)的增加不斷增大。在降雨初期土體位移較小，當(dāng)降雨達(dá)到一定時(shí)間，土體的位移將有一突變。無論是坡腳處的坍塌，還是斜坡破壞過程中滑移帶的形成，土體水平方向的位移都起到了至關(guān)重要的作用。

(2)土體豎向位移具有不均勻性，在大約距坡頂三分之二處取得最大值，這種不均勻性主要是由于土體發(fā)生水平方向的運(yùn)動(dòng)引起的，也正是因?yàn)檫@種不均勻性，為斜坡破壞在表層形成了突破口。

(3)可以采用特征點(diǎn)位移是否發(fā)生突變作為斜坡破壞與否的判斷標(biāo)準(zhǔn)，其數(shù)值模擬結(jié)果與室內(nèi)模型試驗(yàn)結(jié)果吻合較好。

［1］Iverson R M.The physics of debris-flows［J］.Reviews of Geophysics，1997，35(3):245-296.

［2］Klubertanz G.，Laloui，L.，Vulliet，L.，Identification of mechanisms for landslide type initiation of debris flows［J］.Workshop on Mechanics and Velocity of Large Landslides，Engineering Geology，2009，109(1-2):114-123.

［3］Tang，C.，Zhu，J.，Li，W.L.，Rainfall-triggered debris flows following the Wenchuan earthquake［J］.Bulletin of Engineering Geology and the Environment，2009，68:187-194.

［4］CUI Peng，HU Kaiheng，ZHANG Jianqiang et al..Prediction of debris flow area by combing hydrological and inundation simulation methods［J］.Journal of Mountain Science，2011，8(1):1-9.

［5］文海家，張永興，柳源.滑坡預(yù)報(bào)國內(nèi)外研究動(dòng)態(tài)及發(fā)展趨勢［J］.中國地質(zhì)災(zāi)害與防治學(xué)報(bào)，2004，15(1):1-4.WEN Haijia，ZHANG Yongxing，LIU Yuan.A home and abroad investigative-tendency of landslide forecast［J］.The Chinese Journal of Geological Hazard and Control，2004，15(1):1-4.

［6］李海亮，黃潤秋，吳禮舟，等.非均值土坡降雨入滲的耦合過程及穩(wěn)定性分析［J］.水文地質(zhì)與工程地質(zhì)，2013，40(4):70-76.LI Hailiang， HUANG Runqiu， WU Lizhou， et al.Rainfall infiltration coupling process and stability analysis ofa heterogeneous soilslope［J］.Hydrogeology＆Engineering Geology，2013，40(4):70-76.

［7］Guzzetti，F(xiàn).，Peruccacci，S.，Rossi，M..The rainfall intensity-duration control of shallow landslides and debris flows:an update［J］.Lanslides 2008，5(1):3-17.

［8］Berti，M..Simoni，A..Experimental evidences and numerical modelling of debris flow initiated by channel runoff［J］.Landslides，2005，2(3):171-182.

［9］Tohari，A.，Nishigaki，M..Komatsu，M..Laboratory rainfall-induced slope failure with moisture content measurement［J］.Journal of Geotechnical and Geoenvironmental Engineering， 2007， 133(5):575-587.

［10］Major J J，Iverson R M.Debris-flow deposition:effects of pore-fluid pressure and friction concentrated atflow margins［J］.Geological Society of America Bulletin，1999，111(10):1424-1434.

［11］戚國慶，黃潤秋.泥石流成因機(jī)理的非飽和土力學(xué)理論研究［J］.中國地質(zhì)災(zāi)害與防治學(xué)報(bào)，2003，14(3):4-9.QI Guoqing，HUANG Runqiu.Study on genetic and mechanical analysis of debris flow based on unsaturated soils mechanics［J］.The Chinese Journal of Geological Hazard and Control，2003，14(3):4-9.

［12］陳小念，孫海洋.降雨條件下飽和-非飽和土坡的滲流分析［J］.科學(xué)技術(shù)與工程，2007，7(24):115-118.CHEN Xiaonian，SUN Haiyang.Saturated-unsaturated seepage analysis of slope under rainfall［J］.Science Technology and Engineering，2007，7(24):115-118.

［13］陳守義.考慮入滲和蒸發(fā)影響的土坡穩(wěn)定性分析方法［J］.巖土力學(xué)，1997，18(2):8-12.CHEN Shouyi.A method of stability analysis taken effects of infiltration and evaporation into consideration for soil slopes［J］.Rock and Soil Mechanics，1997，18(2):8-12.

［14］朱文彬，等.降雨條件下土體滑坡的有限元數(shù)值分析［J］.巖石力學(xué)與工程學(xué)報(bào)，2002，21(4):509-512.ZHU Wenbin，LIU Baochen.Forming and development process of soil landslide during rainfall［J］.Chinese Journal of Rock Mechanics and Engineering，2002，21(4):509-512.

［15］Sun，Y.，Sakajo S and Nishimura ，M..Appliaction research on a numerical model of two-phase flow in deformation porous medium［C］.Int.Conf.Computer Methods and Advances in Geomech，Wuhan，1997，Ⅱ:1171-1176.

［16］Sun，Y.，Nishigaki，M.and Kohno，I..A study on stability analysis of shallow layer slope due to raining permeation［C］.Proceedings of the 1stInt.Conf.on unsaturated soils.Rotterdam Balkema， 1995:1135-1141.