基于物質點法的白格滑坡堰塞壩形成過程數(shù)值模擬

鐘啟明 吳昊 單熠博 趙鯤鵬

摘要：為準確模擬堰塞壩的形成過程，在物質點方法中引入了滑動速度依賴型的摩擦系數(shù)模型，實現(xiàn)了滑面接觸算法的求解，并與試驗結果對比驗證方法的可靠性。在此基礎上，反演了白格滑坡堰塞壩的形成過程，探討了滑床底摩擦對堰塞壩成壩模式的影響規(guī)律。研究結果表明：白格滑坡堰塞壩形成過程可分為滑坡體加速運動、河谷制動和堆積成壩三階段。白格滑坡速度最大達到45 m/s且最大速度位于滑坡體后緣，而滑坡體前緣受到河床底部摩擦影響，速度不穩(wěn)定，變化幅度大且比滑坡體后緣早20 s制動。不同滑面接觸模式?jīng)Q定了滑坡最終堆積成壩的狀態(tài)：若采用滑動速度依賴型的摩擦系數(shù)模型模擬，則能夠真實反映白格滑坡堰塞壩堆積狀態(tài)；采用滑面的峰值摩擦系數(shù)計算，則導致計算的白格滑坡運動距離減小了32%，不能形成完全堵塞金沙江的堰塞壩；而采用殘余摩擦系數(shù)計算，則導致計算的白格滑坡運動距離增大了12%。研究成果對加深堰塞壩成壩模式的認識具有意義，對于滑坡堰塞壩鏈生災害預測及防災減災應急搶險處置具有一定的參考價值。

關鍵詞：物質點法； 數(shù)值模擬； 接觸模型； 摩擦系數(shù)； 白格滑坡堰塞壩

中圖法分類號： TV122.4

文獻標志碼： A

DOI：10.16232/j.cnki.1001-4179.2024.04.004

0引 言

在地震、降雨和冰川融雪等外部營力作用下河谷兩岸斜坡發(fā)生失穩(wěn)、大變形并運動堆積于谷底，形成自然的擋水體，即滑坡堰塞壩［1-3］?；卵呷麎蔚男纬煽梢l(fā)一連串的次生地質災害［4-6］，如造成上游形成堰塞湖洪澇、下游潰壩洪水泥石流、河道失穩(wěn)、堰塞湖區(qū)二次滑坡并產(chǎn)生涌浪及其級聯(lián)災害效應，這給公路、鐵路工程、管線工程、水電工程等的建設和運營帶來極大威脅，嚴重的，會對所在流域人民生命財產(chǎn)造成災難性影響［7］。

滑坡堰塞壩的形成涉及復雜的多物理過程，特別是滑坡體與滑面的接觸作用，直接影響滑坡體的運動距離，決定了其后所形成的堰塞壩狀態(tài)，如是否完全堵塞河谷、壩體長度和高度等。數(shù)值模擬方法為研究這一復雜過程提供了高效的研究手段。目前數(shù)值方法主要分為基于網(wǎng)格的方法和無網(wǎng)格方法。網(wǎng)格方法，如有限單元法，在模擬滑坡等大位移、大變形、高速運動和碰撞等問題時存在網(wǎng)格畸變問題，故難以再現(xiàn)滑坡堰塞壩形成過程這類復雜的物理過程［8］。相比之下，無網(wǎng)格方法解決這類問題優(yōu)勢更加明顯，如：以離散單元法（DEM）［9］為代表的無網(wǎng)格粒子類方法、光滑粒子流體動力學方法（SPH）［10］、非連續(xù)變形分析方法（DDA）［11-12］等均得到了廣泛的應用。其中，物質點法（MPM）作為一種兼具歐拉和拉格朗日描述的計算方法，綜合了傳統(tǒng)網(wǎng)格類方法和無網(wǎng)格方法的優(yōu)勢，近年來廣泛應用于模擬滑坡失穩(wěn)大變形運動等復雜力學問題上［13-16］。然而，考慮滑動面摩擦系數(shù)隨著滑坡體速度變化的物質點計算方法目前研究尚不足。

本文針對滑坡堰塞壩形成過程，在MPM方法中實現(xiàn)了滑面接觸算法的求解，發(fā)展了滑動速度依賴型的摩擦系數(shù)模型，通過兩個典型的物理模擬試驗驗證了方法的準確性，在此基礎上，模擬了2018年第一次白格堰塞壩形成過程，并探討了不同滑面接觸摩擦系數(shù)條件下的堰塞壩成壩堆積狀態(tài)，進一步豐富了對滑坡堰塞壩形成機制的認識。

1物質點法原理

MPM方法結合了固定的歐拉網(wǎng)格和移動的材料點（拉格朗日粒子）的優(yōu)點，將材料視為連續(xù)體，并將其離散化為材料點［17］。粒子的狀態(tài)，包括質量、位置、動量和變形梯度，通過歐拉網(wǎng)格進行更新，在此網(wǎng)格上計算粒子的運動。本文采用顯式MPM算法，其步驟如下：在每個時間步的開始，將材料點的信息使用常見的形狀函數(shù)傳遞到網(wǎng)格節(jié)點；然后，求解控制方程以獲得節(jié)點加速度，通過這些加速度計算當前時間材料點的加速度、速度和位移；最后，在歐拉網(wǎng)格中更新材料點的位置，以準備計算下一個時間步。

本文采用單點模式模擬干顆粒崩滑體的運動過程，不考慮水相對土體運動過程的影響，土體被視為單一介質，并通過一組材料點進行離散化，每個材料點代表固體和流體相［18］。

1.1平衡方程

由物質點求解顆粒的運動，整個系統(tǒng)滿足質量和動量平衡方程。

1.2本構模型

對于土體材料，本文采用具有莫爾-庫侖屈服準則的彈塑性本構模型。有效應力表達式為

1.3接觸模型

本文MPM算法采用如下思路：首先，對每個物體基于動量方程分別求解其速度，稱之為下一時間步的預測速度；然后，將接觸對象視作整體，基于動量方程求解耦合體的整體速度和受力；最后，通過耦合體受力關系，根據(jù)接觸關系判斷，在此基礎上確定是否對每個接觸對象的預測速度進行校正。

若物體的接觸關系是黏結接觸，則不需要對速度預測值進行修正；若物體的接觸關系是滑移接觸，則需要對速度預測值進行校正。下一時間步的校正速度按式（9）

2模型驗證

2.1顆粒柱崩塌試驗驗證

本節(jié)采用Lube等［20］開展的顆粒柱坍塌物理模擬試驗，驗證基于莫爾庫倫強度準則的MPM模擬得到的顆粒柱坍塌速度和位移。顆粒柱初始堆積的長度和高度均為1 m，顆粒摩擦角為31°，黏聚力為零。試驗開始前，滑坡體右側被閘門限制移動，當打開閘門后，滑坡體開始坍塌、滑動。數(shù)值模擬的設置與物理模擬試驗相同，MPM的網(wǎng)格尺寸為0.2 m。圖1為MPM計算得到的顆粒柱坍塌演化過程。Fern和Soga［21］開展了同樣的數(shù)值模擬研究，并在MPM模型的右上角設置了一個監(jiān)測點記錄顆粒柱崩塌的速度和位移。圖2對比了本文計算得到的顆粒柱坍塌速度和位移結果與Fern和Soga［21］的數(shù)值模擬結果?？梢钥闯觯疚哪M得到的監(jiān)測點速度和位移變化趨勢與前人結果一致。0.5 s時刻，監(jiān)測點滑動達到最大速度2.53 m/s，與前人結果相比誤差在5%以內。但本文模擬得到的監(jiān)測點運動時間略大于前人結果，監(jiān)測點最大位移達到1.5 m，相較于前人結果高了7%?？傮w來說，顆粒柱崩塌試驗結果定量誤差在10%以內，驗證了數(shù)值模擬結果的合理性。

2.2考慮底部摩擦作用的滑坡試驗驗證

本節(jié)采用Mangeney等［22］開展的滑坡侵蝕基底試驗，驗證考慮底部摩擦作用下滑坡演化模擬結果。試驗槽長3 m，高1.2 m，整體呈22°傾角；滑坡體初始堆積高度為0.14 m，長度為0.2 m，滑坡底部有一層厚度為0.004 6 m的可侵蝕層鋪滿整個試驗槽；試驗槽、滑坡體和底部可侵蝕層寬度均為0.1 m。試驗開始前，滑坡體右側被閘門限制移動，當打開閘門后，滑坡體開始坍塌、滑動。

本節(jié)采用物質點法模擬再現(xiàn)了上述試驗，分析不同時刻滑坡的滑動過程以及滑坡體在底部摩擦作用下的運動距離?；麦w采用莫爾庫倫強度準則。材料性質如下：彈性模量2.0×104 kPa，泊松比0.3，顆粒密度滑坡體的高度最初沿水槽壁方向為0.14 m；在重力作用下，滑坡體滑移，水槽壁方向的高度降低，沿水槽方向運動距離增加。在t=0.32 s時，本文數(shù)值模擬得到的滑坡體沿水槽壁方向的高度為0.10 m，低于試驗結果16%，本文數(shù)值模擬得到的滑坡體沿著斜坡滑動距離為0.55 m，高于試驗結果10%；在t=0.64 s時，數(shù)值模擬得到的滑坡體沿水槽壁方向高度降至約0.07 m，低于試驗結果20%，數(shù)值模擬得到的滑坡體沿著斜坡滑動了0.9 m，基本與試驗結果相同。總體來看，數(shù)值模擬計算得到的沿試驗槽高度方向的位置低于試驗結果20%以內，分析主要原因是試驗槽側壁與滑坡體之間存在一定的摩擦作用力，而數(shù)值模擬分析忽略了側壁位置處的接觸關系。但滑坡體整體的運動距離和斷面與Mangeney等［22］試驗結果相符，證明了該方法適用于模擬考慮底部摩擦作用的滑坡演化過程。

3白格滑坡堰塞壩案例與數(shù)值模型

2018年10月10日22：06，西藏自治區(qū)昌都市江達縣波羅鄉(xiāng)白格村和四川省甘孜藏族自治州白玉縣交界處（31°04′56.41″ N，98°42′17.98″ E）的金沙江河道右岸發(fā)生大規(guī)模山體滑坡（后文簡稱為“10.10”白格滑坡），堵塞金沙江干流河道，形成白格堰塞壩［23-24］。

“10.10”白格滑坡發(fā)生時，金沙江江水面高程約為2 880 m，滑坡后緣高程約3 680 m，前緣高程約2 900m，高差達780 m。白格滑坡滑動面傾角在海拔3 400 m以上為25°～38°，平均約為31°，在海拔3 400 m以下較陡，坡度范圍為34°～50°，平均約為39°?；挛镔|主要為巖石和礫石，巖石主要有片麻巖和蛇紋巖兩種類型，兩者具有不同程度的風化［25］。根據(jù)現(xiàn)場調查［26］，金沙江白格滑坡及堰塞壩典型縱斷面如圖5所示，白格滑坡體的主要物理力學參數(shù)列于表1［27］，其中參與計算的基巖參數(shù)按照工程類比法選定［15］。

白格滑坡堰塞壩形成過程的物質點法數(shù)值計算模型如圖6所示。采用三角形網(wǎng)格，基巖網(wǎng)格尺寸為60 m，滑坡體網(wǎng)格為20 m，共生成了8 186個網(wǎng)格和4 202個節(jié)點。白格滑坡體與基巖之間的摩擦作用是本文研究的一個重要方面。Hu等［28］試驗研究發(fā)現(xiàn)，當滑坡體剪切速度超過1 m/s時，滑面的摩擦系數(shù)非常低，甚至能達到0.05。因此在本研究中，參考前人對大光包滑坡的研究成果［15］，如圖6所示定義了兩個滑床接觸表面。接觸表面A位于滑坡體下方，而接觸表面B是原始的山體表面，由于滑坡體滑移，該表面摩擦由峰值摩擦系數(shù)迅速降低為殘余的摩擦系數(shù)。本研究中，為使滑坡體觸發(fā)移動，參考前人方法［29］，取滑坡體內摩擦角的70%作為接觸表面A的殘余摩擦系數(shù)（即0.4）；而接觸表面B的摩擦系數(shù)根據(jù)基巖摩擦系數(shù)設定為初始值0.6。Li等［15］研究表明，滑坡體啟動速度增加至2 m/s后，基巖接觸面摩擦系數(shù)降低，本文參考前人研究取25%的初始值作為接觸表面B的殘余摩擦系數(shù)（即0.15）。如圖6所示，分別在白格滑坡體的前緣、中部和后緣設置了一個速度監(jiān)測點，以分析在滑坡演化過程中的速度變化規(guī)律。數(shù)值計算時間采用顯示積分，時間步長為0.001 s，共計算70 s。

4計算結果分析

白格滑坡堰塞壩形成過程的速度場演化如圖7所示，圖中虛線表示白格堰塞壩實際斷面。白格堰塞壩形成過程大約經(jīng)歷70 s，可分為3個基本階段：t=0～20 s是滑坡體加速運動階段；t=20～30 s，滑坡體運動至河谷中，由于河谷具有較高的粗糙度以及對岸岸坡的空間限制，滑坡體進入河谷制動階段；t=30～70 s，此階段，滑坡體前緣運動至最遠距離，后緣滑坡體逐漸堆積于前緣滑坡體上，速度逐漸減為零，滑坡體在此階段堆積成壩。t=70 s時，數(shù)值模擬得到的白格堰塞壩最終斷面與實際白格堰塞壩斷面對比，可以發(fā)現(xiàn)數(shù)值模擬得到的滑坡體運動至河谷最遠端的距離，即堰塞壩長度，與實際情況較為一致，但數(shù)值模擬得到的白格堰塞壩高度高于實際值。這主要是由于本文采用的是二維的物質點法模擬，而滑坡體運動實際是三維的，在堆積成壩階段，滑坡體由于碰撞對岸斜坡，會向著河谷上下游展寬，目前二維的數(shù)值模擬無法考慮這一情況。此外，白格滑坡堰塞壩的形成實際上存在滑坡體沖向對岸刮鏟對岸崩坡積體后崩回過程，但這一過程涉及新的物質加入運動，需要在數(shù)值分析的質量守恒方程和動量方程中加以考慮，而本文模型尚未考慮這一過程，這也是計算結果存在差異的原因之一。

白格滑坡體前緣、中間及后緣位置處速度特征如圖8所示，最大速度約45 m/s，出現(xiàn)在滑坡體后緣，這一結果基本與Li等［30］基于離散元模擬得到的白格滑坡最大滑速50 m/s吻合?；麦w前緣受到底部摩擦影響，速度變化幅度較大；而滑坡體后緣及中間速度基本呈現(xiàn)三個階段，即快速增大后維持一段較快的運動速度，最后逐漸減小。受滑坡體前緣牽引與后緣擠壓作用，滑坡體中間位置在滑坡體加速運動階段的速度增長幅度最大?；麦w前緣與中部速度在t=50 s前降為零，而滑坡體后緣速度在此后20 s逐漸降低。Hu等［31］基于計算流體動力學研究也揭示出白格滑坡體運動41 s后抵達對岸山體，與本文模擬得到的滑坡體前緣50 s后停止運動的認識基本一致。

進一步對比分析不同滑面接觸模式下白格滑坡最終堆積成壩狀態(tài)。分別采用了峰值摩擦系數(shù)0.4和殘余摩擦系數(shù)0.15模擬白格滑坡成壩過程，計算結果如圖9所示。圖9（a）是考慮速度狀態(tài)相關的滑面接觸模式下白格滑坡最終堆積成壩狀態(tài)，圖9（b）是滑面接觸為峰值摩擦系數(shù)的白格滑坡最終堆積成壩狀態(tài)，圖9（c）是滑面接觸為殘余摩擦系數(shù)的白格滑坡最終堆積成壩狀態(tài)。

由圖可知，若采用考慮速度狀態(tài)相關的滑面接觸模式計算，白格滑坡運動至水平距離2 129.4 m后停止，最大運動距離為841.5 m，這與實際情況基本一致；若采用滑面接觸為峰值摩擦系數(shù)計算，白格滑坡運動至水平距離1 858.2 m后停止，最大運動距離為570.3 m，與圖9（a）計算結果相比，減小了32%；若采用滑面接觸為殘余摩擦系數(shù)計算，白格滑坡運動至水平距離2 231.5 m后停止，最大運動距離為943.6 m，與圖9（a）計算結果對比，增大了12%。故精確考慮速度狀態(tài)相關的滑面接觸模式，能更加精確地計算滑坡體的運動狀態(tài)以及模擬其堵江成壩的模式。

5結 論

本文在MPM方法中實現(xiàn)了滑面接觸算法的求解，提出了與滑動速度相關的摩擦系數(shù)模型，在此基礎上反演了白格滑坡堰塞壩形成過程，得到如下結論：

（1） 通過顆粒柱崩塌試驗和考慮底部摩擦作用的滑坡試驗，驗證了所提出的數(shù)值模型的正確性，MPM方法對滑坡大變形產(chǎn)生的位移和速度的計算均與試驗結果一致。

（2） 白格滑坡堰塞壩形成過程可分為三個階段，即滑坡體加速運動階段、河谷制動階段和堆積成壩階段?；麦w速度最大達到45 m/s且位于滑坡體后緣，而滑坡體前緣受到底部摩擦影響，速度不穩(wěn)定，變化幅度大比滑坡體后緣早20 s制動。

（3） 不同滑面接觸模式?jīng)Q定了滑坡最終堆積成壩狀態(tài)。若采用考慮速度狀態(tài)相關的滑面接觸模式模擬，則能夠真實反映白格滑坡堰塞壩堆積狀態(tài)；采用滑面接觸為峰值摩擦系數(shù)計算，則計算的白格滑坡運動距離減小了32%，不能形成完全堵塞金沙江的堰塞壩；采用滑面接觸為殘余摩擦系數(shù)計算，則導致計算的白格滑坡運動距離增大了12%。

本文研究還存在一些局限性，未來值得進一步深入研究。如本研究采用的是二維物質點法，而溝梁相間真實地形條件下考慮滑坡體基底侵蝕與碰撞刮產(chǎn)增容機制的三維物質點法值得進一步研究。

參考文獻：

［1］石振明，張公鼎，彭銘，等.非均質結構堰塞壩潰決機理模型試驗研究［J］.工程科學與技術，2023，55（1）：129-140.

［2］蔡耀軍，楊興國，張利民，等.堰塞湖風險評估快速檢測與應急搶險技術和裝備研發(fā)研究構想與成果展望［J］.工程科學與技術，2020，52（2）：10-18.

［3］吳昊，年廷凱，單治鋼.滑坡堵江成壩的形成演進機制及危險性預測方法研究進展［J］.巖石力學與工程學報，2023，42（增1）：3192-3205.

［4］ZHONG Q，WANG L，CHEN S，et al.Breaches of embankment and landslide dams-State of the art review［J］.Earth-Science Reviews，2021，216：103597.

［5］崔鵬，郭劍.溝谷災害鏈演化模式與風險防控對策［J］.工程科學與技術，2021，53（3）：5-18.

［6］鐘啟明，錢亞俊，單熠博.崩滑堰塞湖的形成—孕災—致災機理與模擬方法［J］.人民長江，2021，52（2）：90-99.

［7］吳昊.滑坡堵江成壩過程模擬及危險性預測方法研究［D］.大連：大連理工大學，2021.

［8］杜文杰，盛謙，楊興洪，等.基于兩相雙質點 MPM 的滑坡堵江災害鏈生全過程分析［J］.工程科學與技術，2022，54（3）：36-45.

［9］LI D，NIAN T，WU H，et al.A predictive model for the geometry of landslide dams in V-shaped valleys［J］.Bulletin of Engineering Geology and the Environment，2020，79：4595-4608.

［10］WU H，NIAN T，SHAN Z，et al.Rapid prediction models for 3D geometry of landslide dam considering the damming process［J］.Journal of Mountain Science，2023，20：928-942.

［11］ZHANG Y，WANG J，XU Q，et al.DDA validation of the mobility of earthquake-induced landslides［J］.Engineering Geology，2015，194：38-51.

［12］NIAN T，ZHANG Y，WU H，et al.Runout simulation of seismic landslides using discontinuous deformation analysis （DDA） with state-dependent shear strength model［J］.Canadian Geotechnical Journal，2020，57：1183-1196.

［13］ZHAO K L，QIU L C，LIU Y.Two-layer two-phase material point method simulation of granular landslides and generated tsunami waves［J］.Physics of Fluids，2022，34：123312.

［14］ZHAO K L，QIU L C，LIU Y.Numerical study of water wave generation by granular-liquid mixture collapse using two-phase material point method［J］.Applied Ocean Research，2023，137：103608.

［15］LI X，TANG X，ZHAO S，et al.MPM evaluation of the dynamic runout process of the giant Daguangbao landslide［J］.Landslides，2020，18：1509-1518.

［16］ZHU Y，ISHIKAWA T，ZHANG Y，et al.A FEM-MPM hybrid coupled framework based on local shear strength method for simulating rainfall/runoff-induced landslide runout［J］.Landslides，2022，19：2021-2032.

［17］SOGA K，ALONSO E，YERRO A，et al.Trends in large-deformation analysis of landslide mass movements with particular emphasis on the material point method［J］.Géotechnique，2016，66：248-273.

［18］TRONCONE A，PUGLIESE L，CONTE E.Analysis of an excavation-induced landslide in stiff clay using the material point method［J］.Engineering Geology，2022，296：106479.

［19］ABE K，SOGAOGA K，BANDARA S.Material point method for coupled hydromechanical problems［J］.Journal of Geotechnical Geoenvironmental Engineering，2014，140：04013033.

［20］LUBE G，HUPPERT H E，SPARKS R S J，et al.Collapses of two-dimensional granular columns［J］.Physical Review E，2005，72：041301.

［21］FERN E J，SOGA K.The role of constitutive models in MPM simulations of granular column collapses［J］.Acta Geotechnica，2016，11：659-678.

［22］MANGENEY A，ROCHE O，HUNGR O，et al.Erosion and mobility in granular collapse over sloping beds［J］.Journal of Geophysical Research Earth Surface，2010，115（F3）.

［23］FAN X，YANG F，SIVA S S，et al.Prediction of a multi-hazard chain by an integrated numerical simulation approach：the Baige landslide，Jinsha River，China［J］.Landslides，2019，17：147-164.

［24］ZHONG Q，CHEN S，WANG L，et al.Back analysis of breaching process ofBaige landslide dam［J］.Landslides，2020，17：1681-1692.

［25］蔡耀軍，欒約生，楊啟貴，等.金沙江白格堰塞體結構形態(tài)與潰決特征研究［J］.人民長江，2019，50（3）：15-22.

［26］陳祖煜，張強，侯精明，等.金沙江 “10·10” 白格堰塞湖潰壩洪水反演分析［J］.人民長江，2019，50（5）：1-4.

［27］陳祖煜，陳生水，王琳，等.金沙江上游 “11.03” 白格堰塞湖潰決洪水反演分析［J］.中國科學：技術科學，2020，50：763-774.

［28］HU W，HUANG R，MCSAVENEY M，et al.Superheated steam，hot CO2 and dynamic recrystallization from frictional heat jointly lubricated a giant landslide：Field and experimental evidence［J］.Earth and Planetary Science Letters，2019，510：85-93.

［29］WANG J，WANG S，SU A，et al.Simulating landslide-induced tsunamis in the Yangtze River at the Three Gorges in China［J］.ActaGeotechnica，2021，16：2487-2503.

［30］LI D，NIAN T，TIONG RLK，et al.River blockage and impulse wave evolution of the Baige landslide in October 2018：Insights from coupled DEM-CFD analyses［J］.Engineering Geology，2023，321：107169.

［31］HU Y，YU Z，ZHOU J.Numerical simulation of landslide-generated waves during the 11 October 2018Baige landslide at the Jinsha River［J］.Landslides，2020，17：2317-2328.

（編輯：鄭 毅）