非飽和流固耦合物質點方法（CMPM的原理及應用

王斌 周傲 陸盟 王佳俊

摘 要：基于連續介質力學的質量及動量守恒原理，考慮了固-液兩相材料的水力耦合相互作用，推導并提出了基于速度場（v-w formulation）的耦合物質點方法（CMPM），可考慮水力耦 合作用下，飽和/非飽和結構的大變形力學行為.同時，詳細描述了耦合物質點方法的公式推導、矩陣離散過程以及數值實現步驟.隨后，通過對比一維太沙基飽和土固結理論解，以及Liakopou-los非飽和砂土的入滲試驗，初步驗證了耦合物質點方法在水力耦合問題上的準確性.最后，結合 某邊坡在降雨作用下的失穩破壞全過程，分別模擬了持續降雨以及短暫降雨作用下，邊坡的深層破壞及淺層剝蝕現象，進一步驗證了方法在巖土流固全耦合大變形問題中的適用性.

關鍵詞：耦合物質點方法;控制方程;離散過程;數值實現;降雨滑坡

中圖分類號：TU443 文獻標志碼：A

Formulations and Applications of Coupled Material Point Method for Unsaturated Soils

WANG Bin1?，ZHOU Ao1，2，LU Meng3，WANG Jiajun1，4

（1.State Key Laboratory of Geomechanics and Geotechnical Engineering，Chinese Academy of Sciences，Wuhan 430071，China;

2.School of Civil and Environmental Engineering，Hubei University of Technology，Wuhan 430068，China;

3.College of Civil Engineering，Tongji University，Shanghai 200092，China;

4.Faculty of Engineering，China University of Geosciences，Wuhan 430074，China）

Abstract：Based on the mass and momentum conservations of the continuum，this paper presents a novel coupled material point method based on the v-w formulation，which fully considers the large deformation mechanics of saturate and unsaturated structures under the interactions between the solids and the fluids.For simplicity，the gas phase is neglected in the formulation，while an additional item，i.e.degree of saturation，is incorporated in the gov-erning equations to simply study the saturated/unsaturated soils.The detailed derivation process of the governing equations，discretization process in the matrix forms and the computational cycles of CMPM is introduced.Via two benchmark examples，i.e.one-dimensional Terzaghi consolidation solution，and the Liakopoulos test，the validity ofthe CMPM is proven.In the end，a slope failure analysis due to the rainfall infiltration is presented，where both deep progressive and superficial slope failures are shown，further demonstrating that CMPM is a promising tool in simulat-ing hydro-mechanical problems.

Key words：coupled material point method（CMPM）;governing equations;discretization process;numerical simulation;rainfall-induced slope failure

我國是滑坡災害多發的國家，近年來，滑坡造成的年均死亡人數已連續多年超過1000人，滑坡災害 不僅給當地居民的生命財產造成極大損失，有的還 嚴重影響鐵路、公路、水運及水電站等基礎設施的安全運營[1-4].降雨一直是誘發滑坡災害的一個重要因 素，全國 290個縣市地質災害調查結果表明，滑坡在地質災害中所占比例高達51%，而降雨誘發的滑坡比例竟達到滑坡總數的90%[5].因此，準確描述與預 測滑坡的運動過程及其致災范圍，對于滑坡災害的防治具有重要的意義.

目前，對于滑坡運動過程的描述一般采用數值試驗進行仿真模擬[6-7].傳統的有限元方法因為在處理大變形時，容易造成網格的奇異，雅可比矩陣出現 負值，因此不能有效模擬坡體在產生初始破壞后的漸進大變形過程.相比之下，物質點方法（MPM）則 通過使用拉格朗日以及歐拉兩套網格，避免了有限 元中由于大變形而出現的網格畸形問題[8]，同時，其 計算精度及計算效率都比較高，因此，近年來在國內 外受到了諸多研究者的推崇[9-12].

目前，耦合的物質點方法主要可分為兩種形 式[13]：一種是基于v-p格式，也就是以固體顆粒的速度和孔壓作為未知量，這種格式的結構相對簡單，但是不能考慮流體-固體間的相對速度;另一種則是基于v-w 格式，也就是以固體的速度以及流體的速度作為未知量，這種格式考慮更為全面，但是由于流體速度和固體速度同階積分，使得程序的穩定性比較 欠缺.Zhang等[14]最先將物質點方法應用到孔隙介質分析中，但是由于針對流體、固體采用了同階的積分形式，在考慮整體程序穩定性的前提下，這種方法 僅能分析小變形問題，卻犧牲了物質點在大變形分析中的優勢和特點.Jassim等人[15]則基于混合物理論提出了v-w 形式下的流固耦合方法，目前它也是 物質點流固耦合方法中使用最為普遍的，它通過依次求解基于固體以及混合物的動量守恒方程，先后 獲得固體顆粒及流體顆粒的速度，并且通過高斯積分來提高應力計算精度，但是由于物質點在背景網 格中的流動性，很難保證網格的質量和網格中的物質點的質量相等，因此質量難以保證守恒.Bandara和Soga[16]則基于混合物理論，推導了基于兩套顆粒（一套表征固體，另一套則表征流體）的物質點流固耦合控制方程，并且通過引入B-bar型函數在物質點 上直接求解孔壓以防止出現應力自鎖的現象，但是 由于兩套顆粒的使用，程序的穩定性以及計算效率都大幅降低.Yerro等[17]則在隨后的研究中又重新 使用了一套顆粒，通過流體、固體之間的比例來確 定每個顆粒的組分組成，并且進一步推導了在非飽和情況下液相-氣相-固相的全耦合控制方程，但由于氣相的引入，使得這種方法在較大（極端）變形的情況下會引起程序的不穩定，從而降低了方法的適用性.

本文基于連續介質力學的質量以及動量守恒原理，同樣使用了一套顆粒，在考慮極端變形的情況下，通過引入飽和度簡化考慮了非飽和土的力學行為，推導并提出了流固耦合物質點方法（CMPM），可求解流固耦合下巖土結構的極端變形問題.在構建 這樣一個流固耦合控制方程的時候，雖然使用飽和度簡化考慮了氣相，但在實用性上增加了程序的穩 定性.另一方面，降雨滑坡問題則是一個典型的非飽和問題.因此，本文首先比較了經典的一維非飽和土 滲流試驗（Liakopoulos test），以及Terzaghi一維固結試驗（飽和度設置成1的時候，程序即可從非飽和的情況退化為飽和的情況），以證明 CMPM在飽和/非 飽和土力學行為模擬中的適用性，最后，進一步結合 降雨滑坡算例，驗證了方法在解決降雨誘發滑坡模擬問題中的可行性.

1耦合物質點法

基于連續介質力學的質量以及動量守恒原理，本節針對基于速度場（v-w formulation）的物質點流

固耦合系統控制方程進行了詳細的推導.

1.1控制方程

基于連續介質力學的質量以及動量守恒原理，任意連續體的系統控制方程可以表示為

1.2 控制方程的弱形式

至此，基于速度場的流固耦合控制方程已經全 部得出，即流體的動量守恒方程（公式（12）），以及混 合物的動量守恒方程（公式（15））.在相應的公式中，兩邊同時乘以加權函數wh，并在相應積分區域進行積分即可得到動量平衡方程的弱形式.

采用分部積分，將應力項展開，并結合散度定理，水的動量平衡方程為

1.3控制方程的離散化

不同于有限元方法，應力應變選擇在高斯點上 進行積分，物質點法中通常直接在物質點上進行積分.為描述清晰，物質點方法離散示意圖如圖1所示，其中上標“0”表示連續體初始狀態，“1”表示變形 后狀態.物質點用于表征結構本身，儲存和結構本身相關的一切信息以及狀態變量;背景（歐拉）網格僅用于求解控制方程，網格節點信息不進行存儲.通過背景網格上的形函數，物質點上的信息與節點上的信息可以直觀地聯系起來，構造于積分區域上的控 制方程也就可以轉化為物質點上信息的求和.

以物質點質量mp 求解為例，物質點方法將原始構形中的連續體分解成Np個材料點集合，同時認為，每個物質點在空間上占據有相應的體積區域，但不一定需要明確具體的形狀，在對應區域內，存在質量，而在區域外，質量即為零.因此，任一物質點的質量即可表達為連續體的空間密度在其相應區域上的積分，即

值得注意的是，飽和狀態下的流固耦合控制方程，其原理與非飽和情況下的推導計算過程大體相同，只需將飽和度設置為1，控制方程即可直觀地退化成飽和狀態下的動量守恒方程.這里不做詳細介紹，可留給讀者做簡要推導.

1.4 耦合物質點方法（CMPM）的計算步驟

為了求解控制方程，不僅在空間上需要對公式（23）（24）進行離散，同時也需要對其在時間維度上 進行離散.因為將質量矩陣離散為對角形式，其在存 儲時就可存儲為向量的形式.于是，在整個控制方程 求解時，避免了大型的矩陣運算，可以通過向量運算直接求解，而在時間上的迭代則選擇了中心點插值的顯示積分形式.

相比較于有限元分析中更廣泛使用的v-p 格 式，耦合物質點方法采用了基于速度場的v-w 格式，其有兩個優點：1）積分過程中，時間步長的選擇更為寬松[17]，可以提高計算效率，節省計算時間;2）避免了大型矩陣運算，求解過程更為直接.綜合起來，耦 合物質點方法算法步驟總結如下.

1）將儲存在物質點上的信息映射到背景網格 上，初始化節點上所有變量.

2）求解流體（水）相的加速度，見公式（23）.

3）求解固相的加速度，見公式（24）.

4）更新節點上流、固兩項的速度.

5）更新物質點上流、固兩相的速度，然后根據計算得出的速度更新物質點的位置.

6）基于物質點上新的速度，返回映射更新節點 上的速度，基于速度梯度，計算物質點上的應力和水壓力.

7）重置背景網格，開始下一個計算周期.

2 數值算例

2.1一維固結試驗

圖2（a）給出了一維固結試驗示意圖，砂柱高1m、寬0.1m，處于完全飽和狀態.在砂柱上表面處（由于變形小而在邊界節點處）施加10kPa的壓力.假定 兩側和底部不排水，僅允許水從頂面排出.數值計算中，將砂柱離散成10個等大小的四邊形單元，尺寸為0.1m×0.1m.使用各向同性的線彈性材料對砂土 顆粒力學行為進行模擬，材料參數選擇見表1，時間 步長為1.0×10-7? s，共計算1.0×107 步.

圖2（b）給出了耦合物質點法（CMPM）的計算結果，并與Terzaghi理論解進行了對比.對于各種時間 因子 Tv 繪制圖的等時線，如圖2所示，其中

由于應用了低階線性單元，孔隙水壓力在每一個單元內是恒定的，而每個單元使用了上下兩排共 4個積分點，因此孔隙水壓力呈階梯狀分布特征，但總 體而言，孔隙水壓力均勻分布在太沙基理論解的兩 邊，因此，物質點方法在解決飽和土體的流固耦合問題上是適用的.

2.2 非飽和土的入滲試驗

為了驗證考慮非飽和狀態下的耦合物質點法（CMPM）的適宜性，本節采用并分析了Liakopoulos[19]一維非飽和砂土的入滲試驗.試驗裝置如圖3所示，高1m，寬0.1m，試樣采用了Del Monte 砂，材料特性參數參見表2.

試驗具體分為兩個部分，試驗準備階段和試驗 本身.試驗準備階段，在柱樣容器內均勻布滿砂子，然后從砂柱頂部不斷注水，直至整個試樣處于完整 飽和的狀態;試驗開始后，停止注水，兩個側壁不可透水，而水能夠從底部自由排出，這樣，砂柱上半部分逐漸呈現非飽和狀態，沿砂柱高度均勻布置吸力計，用于測量砂柱內部吸力分布.為計算簡便，假定氣體壓力等于大氣壓力.

對于土水特征曲線的選擇，Lewis等[20]認為當飽和度Sw >0.91時，其可采取以下格式：

數值計算中，同樣將砂柱離散成10個大小一樣、尺寸為0.1m×0.1m的四節點四邊形單元，時間 步長為5.0×10-6 s，共計時長120min.試驗過程中，記錄了砂樣在不同時刻的飽和度、垂直位移、毛細壓力隨高度的變化，并與Liakopoulos 試驗進行了對比，結果如圖4所示.

整體而言，計算結果與試驗結果吻合較好，說明了該方法在非飽和土中的適用性.同時，如圖4所示，隨著水從砂柱樣中逐漸滲出，砂柱頂部逐漸形成非飽和狀態，形成毛細水吸力，飽和度變化逐漸減 慢，相應地，水流出速度也逐漸降低.

2.3 降雨滑坡算例

最后，為探究耦合物質點法（CMPM）在滑坡模擬中的適宜性，本節探討了邊坡在兩種不同降雨條件下的失穩破壞過程.

邊坡的幾何尺寸如圖5所示，高度為10m，坡度為45°，土體特性參數如表3所示.數值模型的背景網格尺寸為0.5m×0.5m，每個單元網格基于高斯點位置插入4個物質點，邊坡結構本身表征使用了3720個物質點.初始孔壓設置為-50kPa.降雨邊界條件通過孔壓邊界條件施加，即所有的降雨量都滲入邊坡內部，不形成表面徑流.邊坡內部土體假設為理想彈塑性材料，破壞準則服從摩爾-庫倫強度準則.

2.3.1第一類降雨（持續降雨）直至產生深層破壞

圖6（a）～（f）和圖7（a）～（f）分別展示了不同時刻下以孔隙水壓力和塑性剪切應變不變量為指標的邊坡破壞過程.由圖6可看到，降雨入滲將導致邊坡表面形成濕潤區，孔隙水壓力為0，濕潤區以下的土體未受到影響，孔隙水壓力為初始值-50kPa.隨著入滲時間的增加，濕潤區深度逐漸增加，5 s和10s的入滲深度分別為0.83m 和 1.36 m，最終在邊坡表層土體重度的作用下，約在 12 s時邊坡發生了較大位移，約在15 s時發生顯著變形的失穩破壞，此時入滲深度為1.86 m.隨著進一步入滲的發生，邊坡肩部后緣的土體也發生了坍塌破壞，如圖6（f）所示. 由圖7可看到，降雨入滲過程中邊坡產生了兩條滑移面，一條位于邊坡表層，一條位于邊坡內部約3m的位置.邊坡發生破壞時，土體沿著深層滑面形成了整體滑動趨勢，而后表層滑移面發生加速滑動，塑性應變不變量顯著增加，如圖7（c）（d）所示;直至25 s時，滑體在新位置上達到新的平衡，共產生滑動距離約為12.8 m，如圖7（e）～（f）所示.2.3.2第二類降雨（短暫降雨）產生淺層剝蝕破壞降雨時間設置為8 s，之后任其自由下滲 . 圖8 （a）～（f）和圖 9（a）～（f）分別展示了不同時刻下以孔隙水壓力和塑性剪切應變不變量為指標的邊坡破壞 過程.由圖8和圖9可以看到，在短暫降雨的入滲條 件下，邊坡在降雨結束時首先發生了表層剝蝕破壞，剝蝕土體的滑動距離約為5.8 m;隨著入滲時間的增加，濕潤區逐漸增大，約在12 s時邊坡肩部發生了較大位移，約在15s時土體沿著深層滑移面發生滑動，邊坡發生顯著變形的失穩破壞.25 s時滑坡體的滑動距離約為8.6 m.

綜上所述，本算例中降雨誘發邊坡失穩破壞主 要是由于邊坡表層基質吸力的降低（負孔隙水壓的消散）而引起的.兩類不同的降雨條件很鮮明地揭示了兩種不同的破壞模式，持續降雨條件下，邊坡易形成深層滑面，引起整體的滑動，淺層破壞則更多是由于深層破壞引起的二次破壞，這種破壞模式下的滑動距離主要由初始的深層破壞控制;短暫降雨條件下破壞模式則主要表現為表層的剝蝕破壞，而當初 始的剝落量變得足夠大時，就會進而引起漸進式的深層破壞，這種破壞模式主要依賴于邊坡淺層的破 壞程度.總而言之，CMPM 能夠較好地模擬降雨誘發的邊坡失穩破壞過程.

3結論與展望

基于連續介質力學的質量以及動量守恒原理，考慮固-液兩相非飽和材料的水力耦合作用，本文推 導了基于速度場（v-w formulation）的考慮固-液兩相耦合的系統控制方程.結合 Terzaghi一維固結試驗、Liakopoulos 非飽和土的入滲試驗，以及邊坡在不同 降雨條件下的破壞過程再現，驗證了該計算模型在模擬非飽和土的水力耦合相互作用問題上的可靠 性，對準確評價降雨滑坡的破壞過程具有重要的理論與實際指導意義.

目前的物質點流固耦合模型（CMPM）雖然通過簡化引入飽和度近似模擬了降雨入滲邊坡的過程，可大幅度提高程序的穩定性，但是對于耦合物質點方法的發展及應用依然有很多工作需要開展：

1）針對流固耦合問題，程序的穩定性依然是物質點方法發展的一個瓶頸，引入B-bar 形函數，使用分步法（fractional step algorithm），抑或對物質點上孔 壓進行光滑化等數值手段，需要進一步探究;

2）目前的工作僅討論了孔壓邊界下降雨誘發邊 坡破壞的過程，加入流量邊界，探討不同降雨強度以及降 雨模式下邊 坡 破 壞 過程 也是下一步 需 要研究的.

參考文獻

[1]黃潤秋.20世紀以來中國的大型滑坡及其發生機制[J].巖石力學與工程學報，2007，26（3）：433-454.

HUANG R Q.Large-scale landslides and their sliding mecha-nisms in China since the 20th century [J].Chinese Journal of Rock Mechanics and Engineering，2007，26（3）：433-454.（In Chinese）

[2]王思敬.工程地質學的任務與未來[J].工程地質學報，1999，7（3）：195-199.

WANG S J.Tasks and future of engineering geology[J].Journal of Engineering Geology，1999，7（3）：195-199.（In Chinese）

[3]蔣承菘.中國地質災害的現狀與防治工作[J].中國地質，2000，27（4）：3-5.

JANG C S.Present situation and prevention of geological disasters in China [J].Chinese Geology，2000，27（4）：3-5.（In Chinese）

[4]殷躍平.中國地質災害減災回顧與展望：從國際減災十年到國際減災戰略[J].國土資源科技管理，2001，18（3）：26-29.

YIN Y P.A review and vision of geological hazards in China[J].Management Geological Science and Technology，2001，18（3）： 26-29.（In Chinese）

1. 楊峰，林斐然，魯千國.深圳一工業園山體滑坡 59人失聯[EB/OL].（2015-12-21）[2021-06-29].http：//epaper.bjnews.com.cn/ html/2015-12/21/content_614615.htm？div=-1.YANG F，LIN F R，LU Q G.Fifty-nine people are missing after a landslide at an industrial park in Shenzhen[EB/OL].（2015-12-21）[2021-06-29].

http：//epaper.bjnews.com.cn/html/2015-12/ 21/content_614615.htm？div=-1.（In Chinese）

[6]李芬，鄭威，郭銳.土石混雜高邊坡穩定性分析數值方法研究[J].湖南大學學報（自然科學版），2018，45（S1）：74-79.

LI F，ZHENG W，GUO R.Study on numerical method of stability analysis for earth-rock mixed high slope [J].Journal of Hunan University（Natural Sciences），2018，45（S1）：74-79.（In Chi-nese）

[7]彭文哲，趙明華，肖堯，等.抗滑樁加固邊坡的穩定性分析及最優樁位的確定[J].湖南大學學報（自然科學版），2020，47（5）： 23-30.

PENG W Z，ZHAO M H，XIAO Y，et al.Stability analysis of anti-slide pile reinforced slope and determination of optimal pile posi-tion[J].Journal of Hunan University（Natural Sciences），2020，47（5）：23-30.（In Chinese）

[8]SULSKY D，CHEN Z，SCHREYER H L.A particle method forhistory-dependent materials [J].Computer Methods in AppliedMechanics and Engineering，1994，118（1/2）：179-196.

[9]WI?CKOWSKI Z，YOUN S K，YEON J H.A particle-in-cell solu-tion to the silo discharging problem[J].International Journal for Numerical Methods in Engineering，1999，45（9）：1203-1225.

[10]HU W Q，CHEN Z.Model-based simulation of the synergistic ef-fects of blast and fragmentation on a concrete wall using the MPM[J].International Journal of Impact Engineering，2006，32（12）： 2066-2096.

[11]NAIRN J A.Material point method calculations with explicit

cracks[J].CMES：Computer Modeling in Engineering & Sciences，2003，4（6）：649-664.

[12]WI?CKOWSKI Z.The material point method in large strain engi-neering problems [J].Computer Methods in Applied Mechanics and Engineering，2004，193（39/40/41）：4417-4438.

[13]VAN ESCH J，STOLLE D，JASSIM I.Finite element method for

coupled dynamic flow-deformation simulation [C]//2nd Interna-tional Symposium on Computational Geomechanics.2011.

[14]ZHANG H W，WANG K P，ZHANG Z.Material point method for

numerical simulation of failure phenomena in multiphase porous media [M]//Computational? Mechanics. Berlin，?? Heidelberg： Springer，2009：36-47.

[15]JASSIM I，STOLLE D，VERMEER P.Two-phase dynamic analy-sis by material point method[J].International Journal for Numeri-cal and Analytical Methods in Geomechanics，2013，37（15）： 2502-2522.

[16]BANDARA S，SOGA K.Coupling of soil deformation and pore

fluid flow using material point method[J].Computers and Geotech-nics，2015，63：199-214.

[17]YERRO A，ALONSO E E，PINYOL N M.The material point

method for unsaturated soils [J].Géotechnique，2015，65（3）： 201-217.

[18]VAN GENUCHTEN M T.A closed-form equation for predicting

the hydraulic conductivity of unsaturated soils[J].Soil Science So-ciety of America Journal，1980，44（5）：892-898.

[19]LIAKOPOULOS A C.Transient flow through unsaturated porous

media[D].Berkeley：University of California，Berkeley，1964.

[20]LEWIS R W，SCHREFLER B A.The finite element method in the

static and dynamic deformation and consolidation in porous media[M].Hoboken，NJ，USA：John Wiley and Sons，1998.