基于離散元法的鐵路沿線非飽和邊坡流動(dòng)與沖擊分析

余朝陽 楊彥海 孫大奇 羅剛 謝偉

1.川藏鐵路有限公司，成都 610041；2.四川大學(xué)建筑與環(huán)境學(xué)院，成都 610011；3.中國國家鐵路集團(tuán)有限公司工程管理中心，北京 100038；4.中國鐵道科學(xué)研究院集團(tuán)有限公司鐵道建筑研究所，北京 100081；5.西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，成都 610036；6.西南交大（上海）智能系統(tǒng)有限公司，成都 610031

云南普洱地區(qū)鐵路修建過程中降雨觸發(fā)的崩塌、泥石流災(zāi)害屢見不鮮。在邊坡表面土體達(dá)到飽和之前，往往會(huì)在非飽和狀態(tài)下發(fā)生小規(guī)模崩塌，而此類崩塌很可能發(fā)生在陡峭的鐵路沿線，因此需要對(duì)這些邊坡進(jìn)行適當(dāng)?shù)木S護(hù)。通常可在邊坡表面使用復(fù)合材料［1］或者底部設(shè)置邊坡防護(hù)措施，防止泥沙倒流［2-3］。已有研究認(rèn)為可以將此類山體崩塌看作固體行為，進(jìn)而探討流動(dòng)和沖擊力產(chǎn)生機(jī)制［4］。非飽和土發(fā)生流動(dòng)沖擊現(xiàn)象時(shí)通常在低約束壓力下伴隨大變形和斷裂現(xiàn)象，因此很難尋找滿足這些條件的本構(gòu)模型［5-6］。然而，采用運(yùn)動(dòng)方程追蹤粒子單位行為的離散元法［7］（Discrete Element Method，DEM）進(jìn)行數(shù)值模擬可以有效解決此類問題。當(dāng)前巖土工程中針對(duì)土體的各項(xiàng)參數(shù)分析以及不同破壞形態(tài)的離散元分析也在持續(xù)發(fā)展［8-10］。

本文以云南普洱一鐵路修建現(xiàn)場典型降雨導(dǎo)致的山體崩塌為例，對(duì)鐵路沿線非飽和土的邊坡流動(dòng)和沖擊行為進(jìn)行模型試驗(yàn)，并根據(jù)試驗(yàn)結(jié)果提出一個(gè)簡便計(jì)算最大沖擊力的公式。研究中再現(xiàn)非飽和土在黏聚力作用下產(chǎn)生的邊坡流動(dòng)特性和沖擊行為，并評(píng)估黏附-再黏附粒子接觸模型的適用性，在明確DEM法對(duì)非飽和土邊坡沖擊試驗(yàn)有效性的基礎(chǔ)上進(jìn)行崩塌土的數(shù)值模擬，并對(duì)邊坡實(shí)物尺度崩塌土的流動(dòng)行為和沖擊力產(chǎn)生機(jī)理進(jìn)行研究。

1 地質(zhì)災(zāi)害概況

普洱市坐落于云南省西南部，云貴高原邊緣，海拔高度不均，地形非常復(fù)雜，立體氣候明顯。區(qū)域大部分屬于亞熱帶濕潤氣候，每年平均氣溫在10～15℃，年均降水量在1 500 mm左右，雨水充足。普洱市的基礎(chǔ)建設(shè)中存在大量修建工程，諸如太達(dá)村隧道等。雨季誘發(fā)的鐵路邊坡失穩(wěn)滑塌屢見不鮮，2016年10月及2018年8月，受強(qiáng)降雨的影響，鐵路沿線均發(fā)生大量的山體崩塌，造成了鐵路修建現(xiàn)場職工工棚沖毀，多人被困，人員傷亡慘重。

鐵路沿線屬低中山地貌，洞身分布下第三系漸新至始新統(tǒng)（E2?3）礫巖、砂巖夾泥巖（分布于進(jìn)口至斜井段約4 000 m），白堊系下統(tǒng)曼崗組中段（K1m2）泥巖夾砂巖（分布于斜井至出口段約1 800 m）等。洞身發(fā)育普洱斷裂西支和蠻帕山斷層。雨季正常涌水量為1.2×104m3/d，最大涌水量為1.44×104m3/d。風(fēng)險(xiǎn)區(qū)域廣泛分布于人工活動(dòng)區(qū)及鐵路沿線，主要為粉質(zhì)黏土、黏性土、砂土等，力學(xué)性能較差，結(jié)構(gòu)松散。該區(qū)域邊坡崩塌堆積土堆積于邊坡前緣，剪出口掩埋，一定程度上增加了邊坡穩(wěn)定性，減緩了邊坡蠕滑變形速度。

2 考慮非飽和土黏附強(qiáng)度的簡易黏附模型

2.1 粒子間黏附-再黏附接觸模型

考慮到非飽和土的黏聚力作用，土體崩塌一般存在三種可能性：①土體在崩塌時(shí)破裂。②破裂的土體在與邊坡支護(hù)結(jié)構(gòu)沖擊后堆積在一起。③邊坡底部與邊坡支護(hù)之間存在平滑帶時(shí)，土體堆積在此區(qū)域。研究發(fā)現(xiàn)有必要建立邊坡崩塌后可再次產(chǎn)生黏附效果的模型。因此，研究考慮了黏附-再黏附模型：當(dāng)顆粒之間產(chǎn)生一定的拉伸力時(shí)會(huì)分離，但當(dāng)顆粒再次接觸時(shí)，又會(huì)產(chǎn)生黏聚力［11］。由于本文主要目的不是在微觀層面上對(duì)非飽和土進(jìn)行分析，而是從模型試驗(yàn)和數(shù)值模擬的角度分析非飽和土的行為，因此，主要針對(duì)土體顆粒之間相互接觸時(shí)發(fā)生的宏觀最大拉力作為黏附現(xiàn)象進(jìn)行建模。破壞分析參數(shù)見表1。

表1 破壞分析參數(shù)

DEM允許顆粒間黏附-再黏附的接觸模型（簡化的黏附模型）如圖1所示，圖中與黏聚力有關(guān)的接觸模型用紅色表示，其他基本接觸模型用黑色表示?；窘佑|模型是一個(gè)在法線方向上彈簧和沖床并聯(lián)的模型，并串聯(lián)了一個(gè)分壓器，當(dāng)顆粒之間的距離超過一定值時(shí)，分壓器就會(huì)阻止接觸力發(fā)揮作用。在圖1（a）所示的顆粒半徑相對(duì)應(yīng)的小范圍內(nèi)，使用了接觸模型左側(cè)所示的帶有拉伸彈簧的無量綱DEM附著力影響半徑rc。兩個(gè)解析粒子在半徑尺度的小范圍內(nèi)接近時(shí)，未發(fā)現(xiàn)二者發(fā)生接觸。在實(shí)際的粒子尺度上，接觸瞬間會(huì)產(chǎn)生一個(gè)極大的拉力，此時(shí)被認(rèn)為進(jìn)入接觸狀態(tài)，為了簡化黏附現(xiàn)象，拉伸方向的彈簧系數(shù)被暫定為與壓縮方向的彈簧系數(shù)相同。在黏附過程中，切向的接觸力隨著法向拉伸力的增加而增加。

圖1 顆粒間的黏附-再黏附接觸模型顆粒間接觸狀態(tài)

2.2 基于簡易黏附模型的DEM崩塌分析

為研究附著力影響半徑rc對(duì)臨界高度的影響，設(shè)置了1 m（寬）×1 m（高）的地基模型，從左右元件初始狀態(tài)開始，通過將左側(cè)的邊界面下降0.1 m并釋放應(yīng)力來進(jìn)行崩塌分析。崩塌分析時(shí)的參數(shù)參見表1。粒度分布是最大最小粒徑比為2的均勻分布，在本文中，不僅通過改變黏聚力影響半徑來進(jìn)行模擬，還通過改變顆粒半徑進(jìn)行模擬分析。垂直邊坡崩塌模擬結(jié)果見圖2?？芍街绊懓霃絩c越大，臨界高度就越大。此外，從活動(dòng)塌陷面的角度估計(jì)內(nèi)摩擦角約為25°，該值等于粒子之間摩擦系數(shù)。

圖2 垂直邊坡崩塌的結(jié)果

臨界高度可以通過極限分析法［11-12］或者極限平衡法［13］計(jì)算出的土體黏聚力c得到。垂直邊坡崩塌分析結(jié)果，見圖3?？芍焊街τ绊懓霃絩c和臨界高度幾乎呈線性關(guān)系，并且通過將粒徑減小到一定水平以下，臨界高度會(huì)收斂。由此可說明，DEM顆粒水平的附著力影響半徑（微黏聚力）與土體水平的宏觀黏聚力之間存在線性關(guān)系。

圖3 附著力影響半徑和黏聚力的關(guān)系

在引入簡易黏附模型的情況下，與不考慮黏附模型的普通DEM分析結(jié)果相比，塌陷后的沉積角更陡，沉積的粒子群中存在許多大間隙，研究再現(xiàn)了非飽和土的沉積行為。

2.3 不同坡角的邊坡崩塌分析

對(duì)坡度60°、坡高0.1 m的邊坡進(jìn)行崩塌分析。最大顆粒半徑為5 mm，附著力影響半徑為0.002，邊坡的崩塌分析結(jié)果見圖4?？芍哼吰赂叨刃∮诘扔?.5 m的邊坡沒有崩塌，邊坡高度0.6 m時(shí)邊坡變形顯著，邊坡高度0.7 m時(shí)邊坡出現(xiàn)崩塌。若以有無滑移作為崩塌的指標(biāo)，可以認(rèn)為邊坡的臨界高度為0.5 m。

圖4 不同邊坡高度時(shí)60°邊坡的崩塌解析結(jié)果

不同剪脹角ψ下，邊坡角度和邊坡穩(wěn)定系數(shù)的關(guān)系見圖5。可知：在剪脹角為0°的情況下，當(dāng)邊坡坡度從90°下降到60°時(shí)，穩(wěn)定系數(shù)增加了約1.4倍。由于穩(wěn)定系數(shù)與坡度高度成正比，在黏聚力和單位體積重量保持不變的條件下，穩(wěn)定系數(shù)增加1.4倍意味著臨界坡度高度增加1.4倍。

圖5 邊坡角度和邊坡穩(wěn)定系數(shù)的關(guān)系

本文的結(jié)果顯示，坡度60°邊坡在高度0.4 m時(shí)沒有崩塌，而是在高度0.6 m時(shí)崩塌，這也證實(shí)了采用簡化黏附模型的DEM可以大致再現(xiàn)邊坡土體的黏聚力。

3 山體崩塌模型和非飽和試驗(yàn)分析

根據(jù)對(duì)現(xiàn)場的勘察顯示，區(qū)域內(nèi)平均崩塌深度為1.5 m，坡度多為25°～45°，尤其集中在30°～45°內(nèi)更易觸發(fā)。

3.1 模擬試驗(yàn)的概要

勘察發(fā)現(xiàn)，崩塌堆積體堆積于邊坡前緣，對(duì)破壞起到了抑制作用，一定程度上增加了邊坡穩(wěn)定性，減緩了滑坡蠕滑變形速度。本文進(jìn)行了非飽和土邊坡的沖擊試驗(yàn)，其中鋁板安裝在滑動(dòng)面的底部。模型采用黏土，統(tǒng)一初始條件，含水率為10%，密度為1 560 kg/m3。試驗(yàn)中，測量了崩塌土初始厚度、初始長度和流動(dòng)長度變化時(shí)的壁面沖擊力。使用測力傳感器以100 Hz的頻率對(duì)壁面沖擊力進(jìn)行采樣。

3.2 分析參數(shù)的設(shè)置

參數(shù)參見表1，無量綱粒子附著力影響半徑為0.003。

3.2.1 彈簧系數(shù)

如果整個(gè)顆粒體的泊松比為1/3，彈性波P波與彈性波S波的速度比為2。對(duì)多質(zhì)量點(diǎn)彈簧耦合系統(tǒng)中一維波傳播速度的研究結(jié)果顯示，法線和切線方向的彈簧系數(shù)與彈性P波、S波的平方成正比。根據(jù)這些結(jié)果，確定切線方向的彈簧系數(shù)是法線方向的彈簧系數(shù)的1/4。

彈簧系數(shù)設(shè)定為2×107N/m（法線方向的彈簧系數(shù)），再現(xiàn)了在低封閉壓力下沙子的應(yīng)力傳播速度約為100 m/s。應(yīng)力傳播速度被定義為從自重下落沖擊試驗(yàn)中自重接觸砂子緩沖材料的時(shí)間（安裝在自重上的沖擊加速度計(jì)的上升時(shí)間）到在緩沖材料底部測得的傳遞沖擊力達(dá)到1 kN以上的時(shí)間，也就是應(yīng)力波傳播時(shí)間以及緩沖材料的初始層厚度除以該時(shí)間。

3.2.2 摩擦系數(shù)

采用tanφ=0.466（φ為內(nèi)摩擦角，φ=25°），顆粒與坡面之間的摩擦系數(shù)參照摩擦試驗(yàn)得到。在0.5 m（長）×0.2 m（高）的框架中構(gòu)建了一個(gè)模型土體，當(dāng)土體以大約0.01 m/s的速度在模型的邊坡上滑行時(shí)，根據(jù)稱重傳感器的數(shù)值測量動(dòng)摩擦系數(shù)。在摩擦試驗(yàn)中，土與坡面的摩擦系數(shù)為0.63，在DEM分析中取值一致。作用在底部元件上的切向接觸力的時(shí)間波形與DEM分析結(jié)果一致。

3.2.3 附著力影響半徑

試驗(yàn)時(shí)沒有測量模型土的黏聚力，在崩塌土的初始長度為1.0 m、滑動(dòng)長度為3.0 m的條件下，進(jìn)行了不同初始厚度的邊坡試驗(yàn)，即使崩塌土的初始厚度大于0.4 m也會(huì)出現(xiàn)崩塌，且與支護(hù)撞擊時(shí)的土層厚度大致相同。據(jù)這個(gè)試驗(yàn)和圖4的結(jié)果，本文采用了對(duì)應(yīng)于臨界高0.4 m的模型附著力影響半徑0.003 m。

3.2.4 顆粒大小

最大顆粒直徑與最小顆粒直徑之比設(shè)為2.0，并采用了均勻的顆粒尺寸分布。顆粒尺寸越小，緩沖材料沖擊時(shí)作用在重物上的最大沖擊力越?。?4］。由此可知當(dāng)粒徑小于一定值時(shí)，最大沖擊力收斂到一個(gè)恒定值。邊坡初始厚度0.2 m，邊坡初始長度0.75 m，流動(dòng)長度2.0 m，不同附著力影響半徑時(shí)顆粒大小和沖擊力關(guān)系見圖6。

圖6 不同附著力影響半徑時(shí)顆粒大小和沖擊力關(guān)系

由圖6（a）可知：當(dāng)最大粒子半徑為0.005 m或更小時(shí)，最大沖擊力趨于收斂，最大粒子半徑為0.005 m的沖擊力波形特征與0.002 5 m的沖擊力波形特征幾乎相同。由圖6（b）可知：附著力影響半徑rc＝0的沖擊力波形特征與rc＝0.003 m的趨勢基本相同。

顆粒大小與沖擊壁面崩塌土形狀的關(guān)系見圖7，圖中t為數(shù)值模擬解析時(shí)間。可知，當(dāng)最大粒子半徑小于0.005 m時(shí)，沖擊時(shí)塌陷土趨于收斂。這一結(jié)果也與流體不影響沉積物流動(dòng)試驗(yàn)的分析一致。基于上述分析，本文采用的最大粒子半徑為0.005 m。

圖7 顆粒大小與沖擊壁面的崩塌土的形狀的關(guān)系（解析結(jié)果）

3.3 模型試驗(yàn)和再現(xiàn)分析結(jié)果的比較

對(duì)模型試驗(yàn)結(jié)果和不同外力條件下再現(xiàn)分析結(jié)果進(jìn)行比較，分析參數(shù)見表1。在試驗(yàn)中，負(fù)載是在100 Hz下測量的，因此進(jìn)行了寬度為0.01 s的移動(dòng)平均處理，使分析結(jié)果與其對(duì)應(yīng)。

3.3.1 崩塌土的流動(dòng)距離（沖擊速度）

崩塌土沖擊速度和最大沖擊力的關(guān)系見圖8?？芍鹤畲蟮臎_擊力與沖擊前的土體向支護(hù)的滑動(dòng)速度有關(guān)。最大沖擊力隨著被撞土體的沖擊速度增加呈增加趨勢。

圖8 崩塌土沖擊速度和最大沖擊力的關(guān)系

3.3.2 有無簡易黏附模型對(duì)崩塌土的流動(dòng)和沖擊影響

在崩塌土初始厚度0.2 m，初始長度1.0 m，邊坡長度3.0 m的條件下，在有簡易黏附模型（rc=0.003 m，非飽和模型）和沒有（rc=0，干燥模型）的情況下進(jìn)行了沖擊分析。沖擊波形見圖9?？芍号c非飽和模型相比，干燥模型的沖擊力較小，與試驗(yàn)結(jié)果大致一致，干燥模型的最大沖擊力比非飽和模型小30%左右。

圖9 沖擊波形（簡易黏附模型的影響）

沖擊時(shí)非飽和模型保持了初始層的厚度，而干燥模型則很薄，并向破壞方向延伸。非飽和模型抑制了下坡速度，但由于在支護(hù)撞擊時(shí)土層較厚，干燥模型的最大撞擊力可能比非飽和模型的大。如果黏聚力相對(duì)于崩塌土的初始層厚度小，則干燥模型和非飽和模型沖擊之前的崩塌土層厚度的差相對(duì)較小。由于干燥模型的坡面流動(dòng)速度更大，所以認(rèn)為非飽和模型的沖擊力不一定比干燥模型大。

4 實(shí)物尺寸崩塌土的數(shù)值模擬分析

本文進(jìn)行了基于崩塌案例的全尺寸土體的數(shù)值模擬試驗(yàn)，試驗(yàn)所得沖擊力波形見圖10?？芍?，兩種情況下都產(chǎn)生了大于沉積物壓力的沖擊載荷，最大沖擊力隨著流動(dòng)長度的增加而增加。

圖10 沖擊力波形（實(shí)物尺度的數(shù)值模擬）

依據(jù)試驗(yàn)提出了以下計(jì)算公式

式中：Pmax為最大沖擊力；α為系數(shù)；T為崩塌土的初期厚度；V為崩塌土的沖擊速度；g為重力加速度；L為崩塌土的流動(dòng)距離；θ為斜面的傾斜角度；μ為滑動(dòng)摩擦系數(shù)。

L=10、15 m時(shí)崩塌土的速度分布見圖11。可知，流動(dòng)距離越長，碎土在坡面方向的分布越長，層厚越小。

圖11 崩塌土的速度分布（實(shí)物尺度的數(shù)值模擬）

在L=15.0 m時(shí)，最大沖擊力發(fā)生時(shí)崩塌土的速度和接觸力分布見圖12。可知，超過3 kN的強(qiáng)接觸力分布在距離墻面約1.5 m的地方，而相應(yīng)點(diǎn)的速度分布是一個(gè)速度梯度為5～10 m/s的區(qū)域。這表明最大的沖擊力是由距墻體1.5 m以上區(qū)域的速度梯度產(chǎn)生的接觸力引起的，認(rèn)為沖擊力是由下面的塌陷土與前面沉積的1.5 m塌陷土沖擊力傳遞而來。因此，在總結(jié)墻體最大沖擊力原因時(shí)，最好能考慮到先前沉積的土體的緩沖行為。

圖12 最大沖擊力發(fā)生時(shí)崩塌土的速度與接觸力分布

5 結(jié)論

1）為分析云南普洱某鐵路沿線非飽和土邊坡在持續(xù)降雨下黏聚力產(chǎn)生的影響，基于引入顆粒間的黏附-再黏附模型，可以通過DEM近似地再現(xiàn)黏聚力作用下的非飽和土邊坡的坍陷行為。

2）通過將非飽和土的黏聚力以及土體與坡面的摩擦系數(shù)與各項(xiàng)分析參數(shù)相結(jié)合，可以在坡面流動(dòng)和沖擊試驗(yàn)中獲得最大沖擊力，并提出一個(gè)簡易最大沖擊力計(jì)算公式。

3）當(dāng)支護(hù)上產(chǎn)生最大沖擊力時(shí)，之前堆積在支護(hù)上的滑體可起到緩沖作用，抵御后續(xù)具有較大動(dòng)能土體的沖擊，而現(xiàn)場勘察發(fā)現(xiàn)邊坡前緣的堆積土一定程度上增加了邊坡穩(wěn)定性，減緩了邊坡蠕滑變形速度。