坡面泥石流啟動(dòng)機(jī)理的細(xì)觀數(shù)值模擬研究

王樂榮， 石 崇， 翟淑花， 張一平， 馬金城

（1.河海大學(xué)巖土力學(xué)與堤壩工程教育部重點(diǎn)實(shí)驗(yàn)室，南京 210024；2.河海大學(xué)巖土工程科學(xué)研究所，南京 210024； 3.北京市地質(zhì)研究所，北京 100120）

泥石流是指在山區(qū)或者其他溝谷深壑、地形險(xiǎn)峻的地區(qū)，因?yàn)楸┯?、暴雪或其他自然?zāi)害引發(fā)的山體滑坡并攜帶有大量泥沙以及石塊的特殊洪流. 作為一種突發(fā)性自然災(zāi)害，泥石流暴發(fā)時(shí)往往會(huì)伴隨著巨大的人員傷亡和財(cái)產(chǎn)損失［1-3］. 對泥石流的啟動(dòng)機(jī)理進(jìn)行研究可為泥石流災(zāi)害的預(yù)測和預(yù)防提供參考.

現(xiàn)有的研究泥石流啟動(dòng)機(jī)理的方法主要有三種：現(xiàn)場勘察研究、室內(nèi)模型試驗(yàn)、數(shù)值模擬. Zhao等［4］根據(jù)九寨溝泥石流溝的地形特征分析了泥石流災(zāi)害的歷史和成因，提出了由穩(wěn)定結(jié)構(gòu)、阻塞結(jié)構(gòu)和沉積結(jié)構(gòu)組成的泥石流綜合防治體系. 陳鵬宇和彭祖武［5］應(yīng)用現(xiàn)場勘察研究的方法分析了欒川縣境內(nèi)泥石流的物源特征及啟動(dòng)模式. 甘建軍和儲(chǔ)小東［6］通過野外調(diào)查及室內(nèi)模型試驗(yàn)分析了江西省德安縣泥石流的發(fā)育條件和形成機(jī)理，發(fā)現(xiàn)影響泥石流形成的主要因素是強(qiáng)降雨、豐富的松散物源和陡峻的地貌，同時(shí)利用極限平衡法計(jì)算了泥石流的穩(wěn)定性，并提出了治理的方法和對策. 胡林等［7］通過資料收集、現(xiàn)場調(diào)查等手段，計(jì)算了泥石流的流速、總量、爬高等動(dòng)力學(xué)參數(shù)并按照物源供給對泥石流進(jìn)行分類. 王麗娟等［8］通過野外調(diào)查分析了滑坡型坡面泥石流和沖蝕型坡面泥石流的形態(tài)特征和顆粒組成差異，并結(jié)合地形陡峭指數(shù)和堆積區(qū)流動(dòng)路徑分析了兩種坡面泥石流的流動(dòng)特征. 周健等［9］通過室內(nèi)模型試驗(yàn)研究了不同降雨強(qiáng)度條件下細(xì)粒含量對泥石流形態(tài)的影響. 胡言鵬和劉恩龍［10］通過室內(nèi)模型試驗(yàn)初步探究了降雨特性對黏性泥石流啟動(dòng)過程的影響. 高波等［11］通過室內(nèi)模型試驗(yàn)探究了不同規(guī)模的泥石流啟動(dòng)時(shí)所需的降雨量. 張登等［12］通過室內(nèi)模型試驗(yàn)研究了泥石流啟動(dòng)時(shí)間與坡度、降雨強(qiáng)度之間的關(guān)系. 張永雙等［13］通過室內(nèi)模型試驗(yàn)分析了暴雨作用下寧陜城區(qū)高陡斜坡帶坡面型泥石流的發(fā)育特征和形成過程，探討了泥石流啟動(dòng)的臨界狀態(tài). 隨著計(jì)算機(jī)技術(shù)的發(fā)展，越來越多的學(xué)者［14-16］開始采用數(shù)值模擬方法對泥石流的啟動(dòng)機(jī)理進(jìn)行研究. 姚志雄［17］通過室內(nèi)模型試驗(yàn)和有限元數(shù)值模擬法分析了細(xì)砂型泥石流的誘發(fā)機(jī)制. 張奮翔等［18］采用FLO-2D軟件模擬了不同降雨頻率以及極端降雨條件下泥石流的暴發(fā)過程. 左正金等［19］采用PFC2D軟件模擬了泥石流在降雨條件下的啟動(dòng)機(jī)理及運(yùn)動(dòng)過程.

本研究在典型坡面泥石流室內(nèi)模型試驗(yàn)的基礎(chǔ)上，利用三維顆粒流數(shù)值模擬軟件PFC3D建立了坡面泥石流的數(shù)值模型，并提出等效滲透壓力的簡化方法，然后對不同坡角、不同坡表飽和區(qū)厚度及不同土體抗剪強(qiáng)度條件下坡面泥石流的啟動(dòng)過程和啟動(dòng)機(jī)理進(jìn)行了分析，以期為泥石流地質(zhì)災(zāi)害的預(yù)防提供理論參考.

1 室內(nèi)模型試驗(yàn)

1.1 試驗(yàn)材料

構(gòu)成泥石流物源體的顆粒級配會(huì)直接影響土體的物理特性和水理特性. 為保證試驗(yàn)結(jié)果具有真實(shí)性和代表性，從北京西山某坡面泥石流處采取原狀土試樣，采用振篩機(jī)對試樣進(jìn)行顆粒篩分與級配實(shí)驗(yàn). 篩分后的顆粒級配曲線如圖1 所示.通過計(jì)算得到試樣的不均勻系數(shù)Cu=114.226，曲率系數(shù)Cc=0.541. 根據(jù)土的分類標(biāo)準(zhǔn)可知，該試樣為級配不連續(xù)的土.

圖1 顆粒級配曲線Fig.1 Particle grading curve

1.2 室內(nèi)模型試驗(yàn)?zāi)M坡面泥石流的啟動(dòng)過程

采用泥石流模擬槽體裝置進(jìn)行室內(nèi)模型試驗(yàn)，槽體角度可調(diào)節(jié)，試驗(yàn)斜坡體長2 m，寬1.5 m，用人工降雨機(jī)模擬降雨，試驗(yàn)中降雨以高強(qiáng)度連續(xù)降雨為主. 進(jìn)行室內(nèi)模型試驗(yàn)時(shí)，分別設(shè)置三種坡角的斜坡體，如圖2 所示，左側(cè)斜坡體坡角為30°，中間斜坡體坡角為40°，右側(cè)斜坡體坡角為50°.

由圖2 可知，在室內(nèi)模型試驗(yàn)?zāi)M的整個(gè)過程中，右側(cè)斜坡體先出現(xiàn)坡腳臨空面土體的局部坍塌，而后在持續(xù)降雨作用下，右側(cè)斜坡體的上部土體率先出現(xiàn)裂縫，裂縫的出現(xiàn)加速了土體與雨水的耦合作用，隨著降雨入滲量的增加，裂縫逐漸侵蝕、擴(kuò)展，土體的含水率逐漸加大、孔隙水壓逐漸升高、抗剪強(qiáng)度逐漸降低，當(dāng)土體的下滑力大于抗滑力時(shí)，右側(cè)斜坡體發(fā)生較大規(guī)模的滑塌. 在室內(nèi)模型試驗(yàn)?zāi)M的整個(gè)過程中，中間斜坡體不僅出現(xiàn)了坡腳臨空面土體的局部坍塌現(xiàn)象，還出現(xiàn)了細(xì)顆粒運(yùn)移、坡體表面粗化以及坡體坡度變緩等現(xiàn)象. 左側(cè)斜坡體的整體穩(wěn)定性較好，在降雨雨力濺蝕下，左側(cè)斜坡體出現(xiàn)了坡體表面粗化以及坡體坡度變緩等現(xiàn)象，但未出現(xiàn)坡腳臨空面土體的局部坍塌現(xiàn)象.對比三種斜坡體的坡面泥石流啟動(dòng)過程發(fā)現(xiàn)，坡角與坡面泥石流形成的難易程度呈正相關(guān)，且坡角越大，坡面泥石流造成的危害也越大.

坡面泥石流啟動(dòng)過程的剖面堆積示意圖如圖3所示. 在降雨初期，因?yàn)橥馏w初始含水率低，所以雨水全部下滲進(jìn)入土體，從而導(dǎo)致土體中孔隙水壓力升高，土體抗剪強(qiáng)度降低. 在降雨中期，坡體上緣出現(xiàn)微小的裂縫，裂縫的出現(xiàn)加速了土水的相互作用，使得表層土體逐漸飽和，產(chǎn)生坡面徑流，坡體表面的細(xì)顆粒被表面徑流沖刷至坡腳. 隨著降雨的持續(xù)，即在降雨后期，坡體產(chǎn)生的裂縫逐漸增多、擴(kuò)展，最終貫通形成滑動(dòng)面，當(dāng)土體破壞面上的剪應(yīng)力高于抗剪強(qiáng)度時(shí)，坡面泥石流啟動(dòng).

圖3 坡面泥石流啟動(dòng)過程的剖面堆積示意圖Fig.3 Schematic diagram of profile accumulation of slope debris flow startup process

2 數(shù)值模擬試驗(yàn)

2.1 數(shù)值模型的建立

基于北京西山某坡面泥石流的室內(nèi)模型試驗(yàn)，利用三維顆粒流數(shù)值模擬軟件PFC3D建立坡面泥石流試驗(yàn)槽數(shù)值模型. 先建立滑槽，滑槽坡角α為30°，滑槽水平長2 m，寬1 m，如圖4（a）所示；然后在滑槽內(nèi)生成0.4 m高的土體顆粒，如圖4（b）所示，土體顆粒的最小粒徑為5 mm，坡面最大粒徑為15 mm.

圖4 滑槽和坡面泥石流試驗(yàn)槽的數(shù)值模型Fig.4 Numerical model of chute and slope debris flow test chute

2.2 細(xì)觀參數(shù)的選取

在顆粒流方法中，通過球（ball）來模擬土體顆粒，并通過接觸（contact）來模擬土體顆粒之間的相互作用，通過賦予接觸不同的細(xì)觀參數(shù)來模擬土體顆粒之間的不同力學(xué)性質(zhì).為確定數(shù)值模型中采用的細(xì)觀參數(shù)的合理性，需先建立標(biāo)準(zhǔn)試樣的數(shù)值模型，并設(shè)定數(shù)值模型的初始細(xì)觀參數(shù)，在此條件下進(jìn)行標(biāo)準(zhǔn)試樣的直剪試驗(yàn)，從而可獲得初始的土體抗剪強(qiáng)度，然后通過不斷調(diào)整細(xì)觀參數(shù)，使數(shù)值模擬中標(biāo)定出的土體抗剪強(qiáng)度與室內(nèi)模型試驗(yàn)確定的土體抗剪強(qiáng)度相對應(yīng)，最終即可獲得數(shù)值模型所采用的細(xì)觀參數(shù)，如表1所示.土體顆粒間的接觸模型均選用接觸黏結(jié)模型.

表1 坡面泥石流試驗(yàn)槽數(shù)值模型中采用的細(xì)觀參數(shù)Tab.1 Mesoscopic parameters used in numerical model of slope debris flow test chute

3 降雨影響的簡化分析

在降水作用下雨水入滲至坡體表面的堆積體內(nèi)，降雨過程涉及飽和-非飽和滲流過程.在顆粒流數(shù)值模型中，因降雨引起的飽和-非飽和滲流涉及水-力相互耦合作用，通常需要進(jìn)行復(fù)雜的模擬及計(jì)算，為了簡化分析，本研究僅考慮因降雨引起的坡表飽和區(qū)的滲流作用.

在降雨過程中，雨水下滲進(jìn)入土體，當(dāng)坡面的雨水入滲速度大于土體的雨水入滲速度時(shí)，會(huì)率先在坡表形成飽和區(qū).顆粒流數(shù)值模擬方法中，考慮降雨對坡體的影響時(shí)，會(huì)通過施加附加力來考慮土體顆粒與流體的作用.在飽和土體內(nèi)，由于水頭壓力差的存在，流體對土體產(chǎn)生動(dòng)水壓力和拖曳力，這兩部分力可簡化為等效滲透壓力作用在土體顆粒上［20］，每個(gè)土體顆粒受到流體施加的等效滲透壓力F為：

式中：γw為水的重度，N·m-3；i為水力梯度，其值為土體中兩點(diǎn)水頭之差（H1-H2）與其滲透距離L的比值；V為單個(gè)土體顆粒的體積，m3.

在坡表一定厚度的飽和區(qū)內(nèi)可能形成如圖5（a）和5（b）所示的兩種基本的滲透壓力形式，即矩形滲透壓力和三角形滲透壓力.滲透的邊界情況不同產(chǎn)生的滲透壓力形式也會(huì)不同，同時(shí)還會(huì)存在兩種基本滲透壓力的組合形式，即梯形滲透壓力.對于多數(shù)的邊坡而言，水有向坡腳運(yùn)動(dòng)的趨勢，均會(huì)產(chǎn)生如圖5（c）所示的梯形滲透壓力.

圖5 坡表飽和區(qū)的滲透壓力形式Fig.5 The forms of seepage pressure in the saturated area of the slope surface

土體顆粒在飽和狀態(tài)下還會(huì)受到浮力作用，其大小通過公式（2）計(jì)算得出：

式中：ρw為流體密度，kg·m-3；g為重力加速度，m·s-2；V為土體顆粒的體積，m3.

由于降雨入滲會(huì)引起土體的重度增加，因此降雨影響下的土體重度按照公式（3）計(jì)算：

式中：ρd為土體骨架密度，kg·m-3；ρw為流體密度，kg·m-3；θ為體積含水率，%.

除此之外，土體在飽和狀態(tài)下，水對土體具有軟化作用，會(huì)使其物理力學(xué)參數(shù)存在一定的折減，根據(jù)工程經(jīng)驗(yàn)近似折減15%.

4 坡面泥石流啟動(dòng)過程分析

根據(jù)以上理論，基于如圖4所示的坡面泥石流試驗(yàn)槽數(shù)值模型，采用PFC3D軟件對坡面泥石流的啟動(dòng)過程進(jìn)行數(shù)值模擬.采用施加流體作用力方法考慮降雨對土體的作用，并采用公式（1）～（3）計(jì)算土體的等效滲透壓力、浮力及土體重度，同時(shí)對各參數(shù)進(jìn)行15%的折減來考慮降雨作用的影響，飽和區(qū)滲透壓力采用如圖5（c）所示的梯形滲透壓力，坡表飽和區(qū)厚度為0.05～0.15 m，坡角為30°，摩擦系數(shù)為0.3，數(shù)值模型采用的細(xì)觀參數(shù)如表1所示，模擬結(jié)果如圖6所示.

圖6 坡面泥石流的啟動(dòng)過程Fig.6 Startup process of slope debris flow

從圖6（a）可以看出，受到降雨影響的土體在流體的滲透壓力作用下首先會(huì)在坡體上方形成少量裂隙，此時(shí)坡面形態(tài)完整，如圖6（b）所示.隨著降雨的持續(xù)影響，坡體產(chǎn)生的裂隙逐漸增多，且裂隙主要集中在坡頂和坡中，如圖6（c）所示，裂隙的存在加快了水土作用，使得坡頂?shù)耐馏w向坡腳流動(dòng)，并在坡底堆積，坡頂變薄，如圖6（d）所示.降雨后期，這些裂隙逐步貫通，如圖6（e）所示，在坡體表層形成滑動(dòng)面，處于滑動(dòng)面以上的土體在流體作用下向下滑動(dòng)，最終形成坡面泥石流，如圖6（f）所示.

5 坡面泥石流啟動(dòng)過程的影響因素分析

5.1 坡角對坡面泥石流啟動(dòng)的影響

為了分析坡角對坡面泥石流啟動(dòng)的影響，利用PFC3D軟件分別設(shè)置坡角為20°、30°、40°、50°的坡面泥石流試驗(yàn)槽數(shù)值模型進(jìn)行模擬，數(shù)值模型采用的細(xì)觀參數(shù)如表1所示，飽和區(qū)滲透壓力采用梯形滲透壓力，坡表飽和區(qū)厚度為0.05～0.15 m，摩擦系數(shù)為0.3，數(shù)值模擬結(jié)果如圖7所示.

圖7 不同坡角對坡面泥石流啟動(dòng)的影響Fig.7 Effect of different slope angles on the startup of slope debris flow

從圖7可以看出：坡角為20°時(shí)，土體的滑動(dòng)體積為0.245 m3，平均滑動(dòng)距離為0.113 m，此時(shí)坡體表層土體有向下滑移的趨勢，坡頂土體位移很小，坡腳有少量土體顆粒溢出，坡面形態(tài)完整；坡角為30°時(shí)，土體的滑動(dòng)體積為0.257 m3，平均滑動(dòng)距離為1.053 m，坡頂土體顆粒向前流動(dòng)，坡腳出現(xiàn)較多的土體顆粒堆積，坡面變化明顯；坡角為40°時(shí)，土體的滑動(dòng)體積為0.267 m3，平均滑動(dòng)距離為2.243 m，坡頂土體顆粒被全部帶走，坡面形態(tài)幾乎被破壞；坡角為50°時(shí)，土體的滑動(dòng)體積為0.297 m3，平均滑動(dòng)距離為3.084 m，坡頂土體顆粒被全部帶走，坡頂土體顆?；坪蟪蕡A弧形.以上結(jié)果表明，坡角越大，越容易發(fā)生坡面泥石流，且土體的滑動(dòng)距離越遠(yuǎn)、滑動(dòng)體積越大.

5.2 坡表飽和區(qū)厚度對坡面泥石流啟動(dòng)的影響

為了分析坡表飽和區(qū)厚度對坡面泥石流啟動(dòng)的影響，利用PFC3D軟件分別設(shè)置坡表飽和區(qū)厚度為0.02～0.06 m、0.05～0.15 m、0.10～0.30 m的坡面泥石流試驗(yàn)槽數(shù)值模型進(jìn)行模擬，數(shù)值模型采用的細(xì)觀參數(shù)如表1所示，坡角為30°，摩擦系數(shù)為0.3，飽和區(qū)滲透壓力采用梯形滲透壓力，當(dāng)設(shè)定的坡表飽和區(qū)厚度超出堆積體厚度時(shí)，將滲透壓力施加到滑面上，數(shù)值模擬結(jié)果如圖8所示.

從圖8可以看出：坡表飽和區(qū)厚度為0.02～0.06 m時(shí)，土體的平均滑動(dòng)距離為0.759 m，滑動(dòng)體積為0.268 m3；坡表飽和區(qū)厚度為0.05～0.15 m時(shí)，土體的平均滑動(dòng)距離為1.002 m，滑動(dòng)體積為0.279 m3；坡表飽和區(qū)厚度為0.10～0.30 m時(shí)，土體的平均滑動(dòng)距離為1.752 m，滑動(dòng)體積為0.283 m3.以上結(jié)果表明，坡表飽和區(qū)越厚，土體的滑動(dòng)體積越大、滑動(dòng)距離越遠(yuǎn)，即坡面泥石流的規(guī)模隨坡表飽和區(qū)厚度的增大而增大.

圖8 不同坡表飽和區(qū)厚度對坡面泥石流啟動(dòng)的影響Fig.8 Effect of different thicknesses of slope surface saturation area on the startup of slope debris flow

5.3 土體的抗剪強(qiáng)度對坡面泥石流啟動(dòng)的影響

土體的抗剪切強(qiáng)度用摩擦系數(shù)來表征，為了分析土體的抗剪強(qiáng)度對坡面泥石流啟動(dòng)的影響，利用PFC3D軟件分別設(shè)置土體摩擦系數(shù)為0.1、0.3、0.5的坡面泥石流試驗(yàn)槽數(shù)值模型進(jìn)行模擬，數(shù)值模型采用的細(xì)觀參數(shù)如表1所示，坡角為30°，飽和區(qū)滲透壓力采用梯形滲透壓力，坡表飽和區(qū)厚度為0.05～0.15 m，數(shù)值模擬結(jié)果如圖9所示.

圖9 不同土體抗剪強(qiáng)度對坡面泥石流啟動(dòng)的影響Fig.9 Effect of different shear strengths of soil mass on the startup of slope debris flow

從圖9可以看出：土體摩擦系數(shù)為0.1時(shí)，土體的平均滑動(dòng)距離為1.350 m，滑動(dòng)體積為0.279 m3；土體摩擦系數(shù)為0.3時(shí)，土體的平均滑動(dòng)距離為0.990 m，滑動(dòng)體積為0.257 m3；土體摩擦系數(shù)為0.5時(shí)，土體的平均滑動(dòng)距離為0.020 m，滑動(dòng)體積為0.001 m3.以上結(jié)果表明，土體的抗剪強(qiáng)度越大（即土體的摩擦系數(shù)越大），發(fā)生坡面泥石流的可能性越小，這是因?yàn)樵诮涤曜饔孟拢Σ料禂?shù)較大的土體產(chǎn)生的位移較小，坡體更穩(wěn)定；土體的抗剪強(qiáng)度越?。赐馏w的摩擦系數(shù)越?。?，土體的滑動(dòng)距離越遠(yuǎn)、滑動(dòng)體積越大，發(fā)生的坡面泥石流規(guī)模也越大.

6 結(jié)論

結(jié)合室內(nèi)模型試驗(yàn)，采用三維顆粒流數(shù)值模擬軟件PFC3D建立了坡面泥石流試驗(yàn)槽數(shù)值模型，并通過顆粒流數(shù)值模擬方法分析了不同條件下坡面泥石流的啟動(dòng)過程和啟動(dòng)機(jī)理，得出結(jié)論如下：

1）基于滲流的流體作用過程，在數(shù)值模擬中，對飽和區(qū)的土體采用等效滲透壓力的簡化方法可以近似地模擬降雨對坡面泥石流的作用，并可以提高計(jì)算效率.

2）坡體在滲透力作用下發(fā)生坡面泥石流時(shí)，坡體上首先會(huì)形成一系列微小的裂隙，這些裂隙逐步侵蝕、擴(kuò)展、貫通，在坡表形成滑動(dòng)面，處于滑動(dòng)面以上的土體顆粒在流體作用下向下滑動(dòng)，形成坡面泥石流.數(shù)值模擬結(jié)果與室內(nèi)模型試驗(yàn)結(jié)果大致吻合.

3）坡角越大，則坡體的穩(wěn)定性越差，在相同降雨條件下更容易發(fā)生坡面泥石流.坡面泥石流規(guī)模隨坡表飽和區(qū)厚度的增大而增大，坡表飽和區(qū)越薄，則雨水越容易形成地表徑流，造成坡表沖刷.土體的抗剪切強(qiáng)度用土體的摩擦系數(shù)來表征，摩擦系數(shù)越小的土體越容易滑動(dòng)，則坡面泥石流越容易發(fā)生，且發(fā)生的坡面泥石流規(guī)模越大.