無黏性土坡滑動過程與接觸網(wǎng)絡(luò)演化分析

蔣煉　劉恩龍　田健秋　姜曉瓊

摘要：基于離散元法建立邊坡模型，將顆粒間的黏結(jié)強度降低為零來模擬邊坡的滑動過程。通過記錄不同時間步下邊坡的破壞形態(tài)，從復(fù)雜網(wǎng)絡(luò)的視角分別對邊坡坡頂、坡中和坡腳的平均度、平均最短路徑和網(wǎng)絡(luò)密度進行了分析。結(jié)果表明：靠近邊坡表面的顆粒沿著斜坡向下運動，且越靠近坡面的顆粒位移變化越大;失穩(wěn)后的斜坡通過初始邊坡的中點，形成新的穩(wěn)定邊坡的坡面傾角略小于初始邊坡雙軸壓縮試驗得到的峰值應(yīng)力內(nèi)摩擦角;在初始滑動階段坡頂、坡中和坡腳顆粒的平均度均減小，接觸網(wǎng)絡(luò)結(jié)構(gòu)的穩(wěn)定性減弱，初始階段邊坡的抗剪強度降低，很容易發(fā)生失穩(wěn)變形;隨著時間步的增加，坡頂顆粒的平均最短路徑減小，坡腳顆粒的平均最短路徑增大，坡頂顆粒比坡腳顆粒更容易進行信息交流;隨著時間步的增加，坡頂?shù)慕佑|網(wǎng)絡(luò)變得越來越稠密，坡腳的接觸網(wǎng)絡(luò)變得越來越稀疏，這是滑坡發(fā)生后坡腳處土體變得松散的原因。

關(guān)鍵詞：接觸網(wǎng)絡(luò);離散元法：邊坡失穩(wěn);細觀尺度;復(fù)雜網(wǎng)絡(luò)

中圖分類號：TU449

文獻標(biāo)志碼：A

doi：10.3969/j.issn.1000-1379.2019.04.032

士坡失穩(wěn)是一種常見的地質(zhì)災(zāi)害，其失穩(wěn)破壞是一個發(fā)展變化的過程，可分為局部失穩(wěn)、整體失穩(wěn)滑動、劇烈滑動破壞和滑動后穩(wěn)定4個階段。數(shù)值模擬方法被廣泛應(yīng)用于評估邊坡工程的穩(wěn)定性，研究邊坡的破壞形態(tài)和坡體滑動規(guī)律，常用的方法有極限平衡法[1]、有限單元法和離散元法。基于連續(xù)介質(zhì)力學(xué)的有限單元法很難描述土坡滑動過程，而基于離散單元法的研究僅給出坡體滑動的宏觀規(guī)律[2]，比如運動速度、滑動距離等，沒有提供坡體滑動過程中滑塊之間的接觸關(guān)系以及滑塊的接觸網(wǎng)絡(luò)演化規(guī)律。傳統(tǒng)的極限平衡法，如Bishop法[1]、Sarma法[3]等被廣泛應(yīng)用于分析邊坡的穩(wěn)定性，但有一定局限性。Clough等[4]最早將有限元法用于邊坡的穩(wěn)定性分析，鄭穎人等[5]將有限元強度折減法應(yīng)用于分析邊坡穩(wěn)定性。離散元法是一種基于非連續(xù)介質(zhì)理論的數(shù)值方法，是Cundall[6]最早提出的，這種方法通過定義細觀參數(shù)，不僅可以處理宏觀問題，而且儲存大量土體顆粒信息，因此被廣泛應(yīng)用于巖土問題分析。

盡管有學(xué)者對土坡的滑動過程和規(guī)律進行了大量研究[6-17]，但是已有研究的重點在坡體滑動的宏觀現(xiàn)象描述上?，F(xiàn)階段基于復(fù)雜網(wǎng)絡(luò)的分析方法僅用于分析密實試樣壓縮過程，而從復(fù)雜網(wǎng)絡(luò)視角對邊坡滑動過程中接觸網(wǎng)絡(luò)的演化分析鮮有研究。土坡作為一種復(fù)雜地質(zhì)體，其不同部位的力學(xué)性質(zhì)、位移規(guī)律等不同?；诖耍\用離散元法建立土質(zhì)邊坡，從復(fù)雜網(wǎng)絡(luò)的視角對邊坡滑動過程中坡頂、坡中和坡腳的平均度、平均路徑長度、網(wǎng)絡(luò)密度進行研究，以期為坡體滑動的宏觀規(guī)律進行有益補充。

1 無黏性土坡的顆粒流數(shù)值模擬

1.1 邊坡模型建立

數(shù)值模擬采用接觸黏結(jié)模型，顆粒的最小半徑為3 mm、最大半徑為6 mm，法向剛度為50 MPa、切向剛度為50 MPa.顆粒間初始接觸強度為50 MPa、摩擦系數(shù)為0.3、初始孔隙率為0.16。相似于離心機模型，本次模擬采用放大的重力場法[18]，在200 9重力場下以0為坐標(biāo)原點，向右為x軸正方向，向上為y軸正方向，建立平面直角坐標(biāo)系，生成一個長1.8 m、高0.6 m的矩形模型。根據(jù)開挖的邊坡形狀，削去多余顆粒，生成邊坡的斜率為0.8.坡頂為自由面，底部和側(cè)面的墻無摩擦。為了觀察顆粒的位移情況，對邊坡劃分網(wǎng)格并染色，網(wǎng)格邊長約為0.05 m，一個網(wǎng)格含有約30個顆粒，生成邊坡模型見圖1。

邊坡失穩(wěn)的原因可以分為內(nèi)摩擦角降低、黏聚力降低、孔隙水壓力變化。為了簡化模型，采用顆粒間接觸強度為0來模擬無黏性土坡失穩(wěn)[18]，并對坡體的滑動過程和接觸網(wǎng)絡(luò)的演化進行分析。

1.2 參數(shù)標(biāo)定

離散元法基于顆粒的細觀力學(xué)參數(shù)進行計算，而細觀力學(xué)參數(shù)不能直接與宏觀力學(xué)聯(lián)系，因此需要進行相應(yīng)的數(shù)值試驗將細觀力學(xué)參數(shù)和宏觀力學(xué)聯(lián)系起來，進行雙軸壓縮試驗。顆粒間接觸強度為0的偏應(yīng)力一軸向應(yīng)變和體積應(yīng)變一軸向應(yīng)變關(guān)系曲線見圖2。雙軸壓縮試驗顆粒的尺寸與邊坡顆粒的相同，雙軸壓縮試驗的試樣長度為0.4 m、高為0.8 m，初始孔隙率為0.16，生成約4 000個顆粒，在固結(jié)壓力為50、100、200、400 kPa下的壓縮應(yīng)變速率為1%/min，當(dāng)試樣的軸向應(yīng)變達到200-/0時停止壓縮。雙軸試驗結(jié)果表明，試樣在4種固結(jié)壓力下均表現(xiàn)為應(yīng)變軟化。峰值應(yīng)力和殘余應(yīng)力包絡(luò)曲線見圖3（φ為內(nèi)摩擦角）。試樣的峰值應(yīng)力內(nèi)摩擦角為23.30°殘余應(yīng)力內(nèi)摩擦角為15.8°。

2 邊坡失穩(wěn)分析

2.1 變形特征

在顆粒間黏結(jié)強度降低為0的情況下，邊坡滑動面的破壞形態(tài)見圖4。圖4（a）坡頂和坡腳發(fā)生了變形;圖4（b）坡頂和坡腳發(fā)生了較大變形，坡腳的顆粒向前移動;圖4（c） -圖4（f）邊坡的變形越來越大，邊坡上部的顆粒向下滑動并帶動下部的顆粒發(fā)生變形，部分網(wǎng)格被完全損壞，其中圖4（f）失穩(wěn)的邊坡達到新的穩(wěn)定，并可根據(jù)網(wǎng)格發(fā)生的變形，粗略地得到滑動面的位置。

2.2 位移分析

坡頂線1和坡中線2的位置見圖1，對線上網(wǎng)格在邊坡滑動過程中的位移路徑以O(shè)為原點進行統(tǒng)計，見圖5。由圖5可知，失穩(wěn)后網(wǎng)絡(luò)節(jié)點的位置都低于初始位置。線1，x<1.35 m的網(wǎng)格節(jié)點沿著斜坡向下運動，x>1.53 m的網(wǎng)格節(jié)點向下輕微運動，越靠近坡面的網(wǎng)格節(jié)點位移越大：線2上的網(wǎng)格節(jié)點運動規(guī)律與線1類似。

2.3 失穩(wěn)形狀分析

初始邊坡失穩(wěn)后，顆粒發(fā)生滑動，經(jīng)過不斷調(diào)整，最后達到新的穩(wěn)定，失穩(wěn)后形成新的穩(wěn)定邊坡的形狀見圖6。失穩(wěn)后形成新的穩(wěn)定邊坡并沒有經(jīng)過初始邊坡的坡腳，而是通過初始邊坡的中點。失穩(wěn)后形成新的穩(wěn)定邊坡的坡角為20.0°，略小于圖3中雙軸壓縮試驗的峰值應(yīng)力內(nèi)摩擦角23.3°，大于殘余應(yīng)力內(nèi)摩擦角15.8°。

3 坡體滑動過程中細觀參數(shù)演化

研究邊坡內(nèi)的顆粒時，可以將單個顆?？闯梢粋€點，顆粒間的接觸看成一條邊，需要關(guān)注的是顆粒之間是否有接觸，即兩點之間是否有線連接，因此引入圖論的概念對邊坡進行分析。圖論中的點稱為頂點，連線稱為邊[19]。以y方向位置將邊坡劃分為坡頂（0.50 m

點度[20]是網(wǎng)絡(luò)最基本的特征，指的是與點相關(guān)聯(lián)的邊數(shù)。網(wǎng)絡(luò)中存在一小部分點度值較高的點被稱為高點度點。高點度點越多，意味著網(wǎng)絡(luò)中高點度點與其他點的聯(lián)系越多，結(jié)構(gòu)越穩(wěn)定。平均度（k），也稱為配位數(shù)，代表網(wǎng)絡(luò)中點度的平均值，它和顆粒數(shù)目的關(guān)系為

坡頂、坡中和坡腳的平均度隨時間步的演化見圖7。由圖7可知，坡頂、坡中和坡腳在初始滑動段（時間步為0-1萬步）平均度均減小，接觸網(wǎng)絡(luò)的穩(wěn)定性減弱，初始階段邊坡的抗剪強度降低，邊坡易發(fā)生失穩(wěn)變形。隨著時間步的增加，由于邊坡發(fā)生滑動，顆粒重新排列，因此平均度出現(xiàn)一定波動;當(dāng)邊坡達到新的穩(wěn)定狀態(tài)后，平均度不再隨時間步發(fā)生變化。在發(fā)生滑動前，坡頂?shù)钠骄葹?.63，坡中的平均度為4.03，坡腳的平均度為4.17：邊坡達到新的穩(wěn)定狀態(tài)后，坡頂?shù)钠骄冉档蜑?.59，坡中的平均度降低為3.82，坡腳的平均度降低為3.95。即邊坡發(fā)生滑動時，坡頂?shù)钠骄茸兓钚?，坡中次之，坡腳的平均度變化最大。出現(xiàn)這一現(xiàn)象的原因可能是，邊坡滑動過程中，坡頂?shù)念w粒向下滑動，未發(fā)生滑動的顆粒受擾動程度較小，因此坡頂?shù)钠骄茸兓钚。侯w粒間的黏結(jié)強度為0，顆粒的平衡靠顆粒間的接觸力，前面的顆?；氯ズ螅竺娴念w粒失去平衡，馬上填補前面顆粒的空缺，顆粒在坡腳不斷堆積，坡腳發(fā)生剪切破壞，顆粒間失去部分接觸，大量顆粒重新排列，坡腳平均度變化最大。

3.2 平均最短路徑

最短路徑是指連接兩個節(jié)點i和j之間邊數(shù)最短的路徑[21]，網(wǎng)絡(luò)的平均最短路徑（L）為任意兩個節(jié)點之間距離的平均值，即

接觸網(wǎng)絡(luò)中平均最短路徑可以反映接觸網(wǎng)絡(luò)之間信息的傳遞情況，平均最短路徑越近接觸網(wǎng)絡(luò)之間信息傳遞越容易，平均最短路徑越遠接觸網(wǎng)絡(luò)之間信息傳遞越困難，因此研究坡頂、坡中和坡腳平均最短路徑變化，見圖8。從圖8可以看出，隨著時間步的增加，坡頂平均最短路徑減小，坡腳平均最短路徑不斷增大，坡中平均最短路徑略有增大但整體變化不大，這一現(xiàn)象與顆粒數(shù)目有關(guān)。顆粒數(shù)目越多，任意兩個顆粒之間的距離越遠，接觸網(wǎng)絡(luò)的信息傳遞越困難。坡頂顆粒數(shù)目減少，接觸網(wǎng)絡(luò)的信息傳遞變得容易;坡腳顆粒數(shù)目增多，接觸網(wǎng)絡(luò)信息傳遞變得困難。

3.3 網(wǎng)絡(luò)密度

實際接觸網(wǎng)絡(luò)的規(guī)模隨時間發(fā)生變化，因此引入網(wǎng)絡(luò)密度來研究接觸網(wǎng)絡(luò)疏密性變化，見圖9。網(wǎng)絡(luò)密度（p）指接觸網(wǎng)絡(luò)中實際存在的邊數(shù)與最大可能邊數(shù)之比[22]，即

圖9（a）表明，隨著時間步的增加，坡頂顆粒的接觸網(wǎng)絡(luò)變得越來越稠密，坡腳顆粒的接觸網(wǎng)絡(luò)變得越來越稀疏，坡中顆粒接觸網(wǎng)絡(luò)的疏密性變化較小，且坡頂顆粒的接觸網(wǎng)絡(luò)比坡腳顆粒的接觸網(wǎng)絡(luò)相對稠密，主要原因是堆積在坡腳處的滑坡體多為松軟土體。

由式（3）可知，網(wǎng)絡(luò)密度與顆粒數(shù)目和接觸網(wǎng)絡(luò)的邊數(shù)有關(guān)，因此對坡頂、坡中和坡腳的顆粒數(shù)目與接觸網(wǎng)絡(luò)的邊數(shù)（M）進行統(tǒng)計，見圖9（b） -圖9（d）。

由圖9可以看出，實際網(wǎng)絡(luò)演化過程中，坡頂處顆粒數(shù)目和接觸網(wǎng)絡(luò)的邊數(shù)呈線性變化，斜率為1.01，表明接觸網(wǎng)絡(luò)隨著節(jié)點數(shù)的增加變得越來越稠密：坡中顆粒數(shù)目和接觸網(wǎng)絡(luò)邊數(shù)的變化較為復(fù)雜，一方面坡頂顆粒向下滑動，引起坡中顆粒數(shù)目增加，另一方面其自身向坡腳滑動引起顆粒數(shù)目減少，因此接觸網(wǎng)絡(luò)的變化比較復(fù)雜：坡腳顆粒數(shù)目和接觸網(wǎng)絡(luò)的邊數(shù)呈線性變化，斜率為0.83，表明接觸網(wǎng)絡(luò)隨著節(jié)點數(shù)的增加變得越來越稀疏。

4 結(jié)論

（1）邊坡發(fā)生破壞時，坡頂和坡腳的顆粒都發(fā)生滑動。靠近邊坡表面的顆粒沿著斜坡向下運動，且越靠近坡面的顆粒位移變化越大。失穩(wěn)后形成新的穩(wěn)定邊坡并沒有經(jīng)過初始邊坡的坡腳，而是通過初始邊坡的中點，失穩(wěn)后形成新的穩(wěn)定邊坡的坡角略小于初始邊坡雙軸壓縮試驗得到的峰值應(yīng)力內(nèi)摩擦角，但大于殘余應(yīng)力內(nèi)摩擦角。

（2）在初始滑動階段坡頂、坡中和坡腳顆粒的平均度均減小，接觸網(wǎng)絡(luò)結(jié)構(gòu)的穩(wěn)定性減弱，初始階段邊坡的抗剪強度降低，邊坡很容易發(fā)生失穩(wěn)變形。隨著時間步的增加，顆粒間發(fā)生滑動，平均度出現(xiàn)一定波動，當(dāng)達到新的穩(wěn)定狀態(tài)后，平均度不再隨時間步發(fā)生變化。邊坡發(fā)生滑動時，坡頂?shù)钠骄茸兓钚?，坡中次之，坡腳的平均度變化最大。

（3）隨著時間步的增加，坡頂顆粒的平均最短路徑減小，坡中顆粒的平均最短路徑變化較小，坡腳顆粒的平均最短路徑增大，即坡頂顆粒比坡腳顆粒更容易進行信息交流。

（4）坡頂和坡中接觸網(wǎng)絡(luò)密度隨著時間步的增加而增大，坡腳接觸網(wǎng)絡(luò)密度隨時間步的增加而減小，因此隨著時間步的增加，坡頂和坡中接觸網(wǎng)絡(luò)變得越來越稠密，坡腳接觸網(wǎng)絡(luò)變得越來越稀疏，這是滑坡發(fā)生后坡腳處土體變得松散的原因。

參考文獻：

[1]BISHOP A W.The Use of the Slip Circle in the Stability A—nalysis of Slopes[J]. Geotechnique， 1955，5（1）：7-17.

[2]WANG C，TANNANT D D，LILLY P A.Numerical Analysis ofthe Stability of Heavily Jointed Rock Slopes Using PFC2n[J].Intemational Joumal of Rock Mechanics&Mining Sciences，2003，40（3）：415-424.

[3] SARMA S K.Stability Analysis of Embankments and Slopes[J]. Joumal of Ceotechnical Engineering Division， 1979，105（12）：1511-1524.

[4] CLOUCH R W， WOODWARD R J.Analysis of EmbankmentStresses and Deformations[J].Soil Mechanics&FoundationDivision， 1967， 93（4）： 529-549.

[5] 鄭穎人，趙尚毅，張魯渝，用有限元強度折減法進行邊坡穩(wěn)定分析[J].中國工程科學(xué)，2002，4（10）：57-61.

[6]CUNDALL P A.A Computer Model for Simulating Progres-sive， Large-Scale Movements in Block Rock Systems[ C]//Nancy： Symposium of Intemational Society of Rock Mechan-ICS.1971： 8-11.

[7]CHEN H，LIU S H.Slope Failure Characteristics and Stabi-lization Methods[J].Canadian Ceotechnical Journal， 2007，44（4）：377-391.

[8]KIM J S，KIM J Y，LEE S R.Analysis of Soil Nailed EarthSlope by Discrete Element Method[ J]. Computers&Ceotechnics， 1997， 20（1）：1-14.

[9]JIANC M， JIANC T，CROSTA G B，et al.Modeling Failureof Jointed Rock Slope with Two Main Joint Sets Using aNovel DEM Bond Contact Model[J].Engineering Ceology，2015， 193：79-96.

[10]SCHOLTES L，DONZEF V.A DEM Analysis of Step-PathFailure in Jointed Rock Slopes[ J]. Comptes RendusMecanique， 2015（2）：155-165.

[11] 周健，王家全，曾遠，等，顆粒流強度折減法和重力增加法的邊坡安全系數(shù)研究[J].巖土力學(xué)，2009， 30（6）：1549-1554.

[12]周健，王家全，曾遠，等，土坡穩(wěn)定分析的顆粒流模擬[J].巖土力學(xué)，2009，30（1）：86-90.

[13] 曾锃，張澤輝，楊宏麗，等，基于邊坡漸進破壞特征對傳統(tǒng)極限平衡法幾點假設(shè)的合理分析[J].巖土力學(xué)，2012，33（增刊1）：146-150.

[14]ROTHENBURC L，KRUYT N P.Critical State and Evolu-tion of Coordination Number in Simulated CranularMaterials[J].Intemational Journal of Solids&Structures，2004， 41（ 21）： 5763-5774.

[15] PETERS J F，MUTHUSWAMY M， WIBOWO J，et al.Characterization of Force Chains in Cranular Material[J].Physical Review E，2005， 72（4）： 138-148.

[16]TORDESILLAS A，0'SULLIVAN P，WALKER D M， et al.Evolution of Functional Connectivity in Contact and ForceChain Networks： Feature Vectors，K- Cores and MinimalCycles[J].Comptes Rendus Mecanique， 2010， 338（ 10- 11）：556-569.

[17] SILBERT L E，CREST G S，LANDRY J W. Statistics ofthe Contact Network in Frictional and Frictionless CranularPackings[J].Physical Review E Statistical Nonlinear&Soft Matter Physics， 2002， 66（1）：114-129.

[18]JLANC M， MURAKAMI A.Distinct Element Method Ana-lyses of Idealized Bonded-Cranulate Cut Slope[J].Cranular Matter， 2012（ 14）： 393-410.

[19]李明哲，金俊，石瑞銀，圖論及其算法[M].北京：機械工業(yè)出版社，2010：4-5.

[20]DOROCOVTSEV S N， MENDES J F F.The Shortest Pathto Complex Networks[J].Physics， 2004， 71： 47-53.

[21]王小帆，李翔，陳關(guān)榮，網(wǎng)絡(luò)科學(xué)導(dǎo)論[M].北京：高等教育出版社，2012：87-88.

[22]何大韌，劉宗華，汪秉宏，復(fù)雜系統(tǒng)與復(fù)雜網(wǎng)絡(luò)[M].北京：高等教育出版社，2009：126-127.

【責(zé)任編輯呂艷梅】