非飽和流固耦合物質(zhì)點方法（CMPM的原理及應(yīng)用

王斌 周傲 陸盟 王佳俊

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

關(guān)鍵詞：耦合物質(zhì)點方法;控制方程;離散過程;數(shù)值實現(xiàn);降雨滑坡

中圖分類號：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

我國是滑坡災(zāi)害多發(fā)的國家，近年來，滑坡造成的年均死亡人數(shù)已連續(xù)多年超過1000人，滑坡災(zāi)害 不僅給當?shù)鼐用竦纳敭a(chǎn)造成極大損失，有的還 嚴重影響鐵路、公路、水運及水電站等基礎(chǔ)設(shè)施的安全運營[1-4].降雨一直是誘發(fā)滑坡災(zāi)害的一個重要因 素，全國 290個縣市地質(zhì)災(zāi)害調(diào)查結(jié)果表明，滑坡在地質(zhì)災(zāi)害中所占比例高達51%，而降雨誘發(fā)的滑坡比例竟達到滑坡總數(shù)的90%[5].因此，準確描述與預(yù) 測滑坡的運動過程及其致災(zāi)范圍，對于滑坡災(zāi)害的防治具有重要的意義.

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

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

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

1耦合物質(zhì)點法

基于連續(xù)介質(zhì)力學(xué)的質(zhì)量以及動量守恒原理，本節(jié)針對基于速度場（v-w formulation）的物質(zhì)點流

固耦合系統(tǒng)控制方程進行了詳細的推導(dǎo).

1.1控制方程

基于連續(xù)介質(zhì)力學(xué)的質(zhì)量以及動量守恒原理，任意連續(xù)體的系統(tǒng)控制方程可以表示為

1.2 控制方程的弱形式

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

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

1.3控制方程的離散化

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

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

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

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

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

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

1）將儲存在物質(zhì)點上的信息映射到背景網(wǎng)格 上，初始化節(jié)點上所有變量.

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

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

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

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

6）基于物質(zhì)點上新的速度，返回映射更新節(jié)點 上的速度，基于速度梯度，計算物質(zhì)點上的應(yīng)力和水壓力.

7）重置背景網(wǎng)格，開始下一個計算周期.

2 數(shù)值算例

2.1一維固結(jié)試驗

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

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

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

2.2 非飽和土的入滲試驗

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

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

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

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

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

2.3 降雨滑坡算例

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

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

2.3.1第一類降雨（持續(xù)降雨）直至產(chǎn)生深層破壞

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

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

3結(jié)論與展望

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

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

1）針對流固耦合問題，程序的穩(wěn)定性依然是物質(zhì)點方法發(fā)展的一個瓶頸，引入B-bar 形函數(shù)，使用分步法（fractional step algorithm），抑或?qū)ξ镔|(zhì)點上孔 壓進行光滑化等數(shù)值手段，需要進一步探究;

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

參考文獻

[1]黃潤秋.20世紀以來中國的大型滑坡及其發(fā)生機制[J].巖石力學(xué)與工程學(xué)報，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]王思敬.工程地質(zhì)學(xué)的任務(wù)與未來[J].工程地質(zhì)學(xué)報，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]蔣承菘.中國地質(zhì)災(zāi)害的現(xiàn)狀與防治工作[J].中國地質(zhì)，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]殷躍平.中國地質(zhì)災(zāi)害減災(zāi)回顧與展望：從國際減災(zāi)十年到國際減災(zāi)戰(zhàn)略[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. 楊峰，林斐然，魯千國.深圳一工業(yè)園山體滑坡 59人失聯(lián)[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]李芬，鄭威，郭銳.土石混雜高邊坡穩(wěn)定性分析數(shù)值方法研究[J].湖南大學(xué)學(xué)報（自然科學(xué)版），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]彭文哲，趙明華，肖堯，等.抗滑樁加固邊坡的穩(wěn)定性分析及最優(yōu)樁位的確定[J].湖南大學(xué)學(xué)報（自然科學(xué)版），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.