基于物质点-特征有限元耦合方法的向后侵蚀管涌模拟

王兆南; 王刚; 金伟

doi:10.11779/CJGE20230207

基于物质点-特征有限元耦合方法的向后侵蚀管涌模拟 English Version

王兆南^{1, 2,},
王刚^{1, 2},
金伟^{1, 3}

1.
山地城镇建设与新技术教育部重点实验室, 重庆 400045
2.
重庆大学土木工程学院, 重庆 400045
3.
中国电建集团成都勘测设计院有限公司, 四川成都 610072

基金项目:

国家自然科学基金面上项目 52079126

详细信息

作者简介:
王兆南(1994—)，男，博士研究生，主要从事岩土力学数值方法和本构模型等方面的研究工作。E-mail: znwang@cqu.edu.com

中图分类号: TU411
计量
- 文章访问数: 326
- HTML全文浏览量: 72
- PDF下载量: 83
出版历程
- 收稿日期: 2023-03-08
- 网络出版日期: 2024-06-04
- 刊出日期: 2024-05-31

Simulation of backward erosion piping based on coupled material point-characteristic finite element method

WANG Zhaonan^{1, 2,},
WANG Gang^{1, 2},
JIN Wei^{1, 3}

1.
Key Laboratory of New Technology for Construction of Cities in Mountain Area, Chongqing University, Chongqing 400045, China
2.
School of Civil Engineering, Chongqing University, Chongqing 400045, China
3.
Chengdu Engineering Corporation Limited, Power China, Chengdu 610072, China

摘要

摘要: 向后侵蚀管涌是汛期堤坝中一种常见的渗透破坏形式，多发生在下游存在无保护出水口的二元结构堤基中。由于土体渗流在出水口处存在较高的水力梯度，使得固、液界面附近的土体极易被侵蚀带走，一旦侵蚀土体的上层存在不透水的黏土层时，便会形成不断向上游迎水侧发展的管涌通道，最终导致堤坝出现失稳破坏。基于向后侵蚀管涌的发展过程，采用局部水力梯度作为管涌发展的判别准则，结合饱和孔隙介质的耦合物质点-特征有限元法，发展了一个能够模拟向后侵蚀管涌的新方法。新方法将求解域内的物质点分为3种粒子类型，对满足管涌触发条件的粒子进行删除，以此表示被侵蚀带走的土颗粒。由于算法中的流体部分采用了广义的Navier-Stokes方程进行描述，因此新方法能够同时计算孔隙水渗流和管道流体的自由流动。通过对小尺度模型试验的计算，验证了新方法在向后侵蚀管涌问题中的适用性。
- 向后侵蚀管涌 /
- 耦合物质点-特征有限元法 /
- 临界水力梯度
Abstract: The backward erosion piping is a common form of seepage failure in embankments during flood seasons, and it mostly occurs in the dual-structure foundations with unprotected outlet downstream. Due to the high hydraulic gradient of the soil seepage at the water outlet, the soil near the solid-liquid interface is easy to be eroded away. Once there is an impervious clay layer on the upper layer of the eroded soil, a piping channel will be formed to continuously develop to the upstream side, and eventually lead to the instability and failure of the embankment. Based on the coupled material point - characteristic finite element method, a novel modeling approach for the backward erosion piping is developed by employing the local hydraulic gradient as the triggering criterion of piping. The novel approach divides the particles within the solution domain into three types, and deletes the particles that meet the triggering conditions of piping to represent the granular taken away by erosion. Since the fluid phase is described by the generalized Navier-Stokes equation, the proposed approach can simultaneously calculate the seepage of pore water and the free flow of piping channel. Finally, the small-scale erosion experiments are provided to perform the applicability of the proposed approach.
- backward erosion piping /
- coupled material point-characteristic finite element method /
- critical hydraulic gradient

HTML全文

0. 引言

近年来,盾构施工技术已经在地铁隧道中得到了广泛的应用,如何正确预测城市地下工程施工对周围地层变形尤其是地表位移显得异常重要^[1-2]。纵观目前的地下施工技术,理论研究明显不足,实际施工大多数依赖于工程经验,因此,许多学者对盾构隧道施工引起的地层变形进行了研究。Peck^[3]于1969年提出了一种经验公式来预测地表沉降变形,许多工程实例已经证明了它的使用效果,但该公式中参数的选取往往取决于实际工程经验,未考虑土层条件,也不能预测深部土层位移情况；吴昌胜等^[4]研究了不同盾构直径对隧道地层损失率和地表沉降的影响；张冬梅等^[5]在研究隧道地表沉降时考虑了衬砌局部渗流作用的影响。以上研究均没有考虑土体的小应变特性,而大量工程实测数据表明,实际工程诸如地铁、基坑变形等中岩土体应变都很小,主要集中在0.01%～0.3%^[6-7],在此情况下,土体的剪切刚度是随着应变的增加而减小的,土体在小应变阶段的刚度要远大于较大应变阶段的刚度。众多的岩土试验^[8-9]也表明,土体的剪切模量与剪应变之间存在一种非线性关系,图1是相应的归一化刚度退化曲线。

图 1 刚度衰减曲线和应变范围

Figure 1. Curves of stiffness decreasing with strain and ranges of strain

下载: 全尺寸图片幻灯片

Lee等^[10-11]的研究表明,利用理想弹塑性本构模型（莫尔-库仑模型）来预测实际盾构施工变形,其结果偏于保守。因此,土体小应变特性会对预测盾构施工地层尤其是地表变形产生重要的影响。

鉴于此,本文从小应变塑性硬化本构（PH-SS模型）^[12]和莫尔-库仑本构（MC模型）出发,对盾构施工地层变形规律展开了研究。首先基于三轴模拟试验证实了PH-SS模型反映土体小应变特征的有效性；依托上海某一软土隧道工程,比较PH-SS模型和MC模型在模拟盾构施工中地表及深部土体变形规律的异同；最后结合现场监测数据,说明PH-SS模型在盾构施工地表变形预测中的先进性和合理性。

1. PH-SS模型的介绍

1.1 有限差分程序中的小应变本构模型

针对有限差分程序中的小应变特性研究,尹骥^[13]在室内试验的基础上,利用FLAC^3D中内嵌的fish语言将试验结果应用到有限差分程序的二次开发中,但这往往难度较大,同时自定义的本构模型往往受到试验结果较大的影响。本文利用FLAC^3D内置的Plastic Hardening-Small Strain model（简称PH-SS模型）^[12]来反映土体的小应变特性。

PH-SS模型以硬化土模型（PH模型）为基础,充分考虑了BENZ^[8]提出的卸载—再加载应变相关刚度关系,将土体的小应变特性考虑到有限差分程序中,最终用小应变参数G₀（初始剪切模量）和γ₇₀（剪应变阈值）来反应土体的小应变特性。除此之外,PH-SS还包含了11个硬化土模型参数：有效黏聚力c′、膨胀角ψ、有效内摩擦角 $φ^{'}$ 、三轴固结排水剪切试验的参考割线模量 $E_{50}^{ref}$ 、固结试验的参考切线模量 $E_{oed}^{ref}$ 、三轴固结排水卸载再加载试验的参考加卸载模量 $E_{ur}^{ref}$ 、与模量应力水平相关的幂指数m、参考应力p^ref、破坏比R_f、泊松比 $ν$ 、正常固结时的静止侧压力系数K₀。

PH模型以增量弹性定律来描述弹性行为；采用的是剪切屈服流动准则,确定剪切硬化开始和发展的剪切屈服函数可以表示为

$f^{s} = \frac{E_{ur}}{E_{i}} \cdot \frac{q_{a} q}{(q_{a} - q)} - q - \frac{E_{ur} γ^{p}}{2} = 0$ 。

(1)

式中 E_i=2E₅₀/(2-R_f)是初始变形模量； $q_{a}$ = $q_{f}$ /R_f,破坏比R_f是一个比1还小的值（通常取0.9）；而最终偏应力q_f的定义和莫尔-库仑准则是一致的,q是偏应力,PH模型还引入一种截止因子f_cut（Cut-off）以防止偏应力过大； $γ^{p}$ 是剪切硬化参数（内部变量之一）, $Δ γ^{p} = - (Δ ε_{1}^{p} + Δ ε_{2}^{p} + Δ ε_{3}^{p})$ ,由此可见, $γ^{p}$ 与3个塑性应变有关。

PH模型中的体积塑性应变和剪切塑性应变之间的流动规律如下：

$Δ ε_{v}^{p} = \sin ψ_{m} Δ γ^{p}$ 。

(2)

式中, $\sin ψ_{m}$ 为动力膨胀角,应该小于或等于定义的最终膨胀角ψ。

PH模型的剪切势函数为

$g^{s} = (- 1 + sin ψ_{m}) + (1 + \sin ψ_{m}) σ_{3} / 2$ 。

(3)

式中, $σ_{1}$ 为最大主应力, $σ_{3}$ 为最小主应力。

为了避免土体达到临界孔隙状态e_max时发生过度膨胀,对膨胀角进行了稍微修改。主要是采用Schanz等^[14]提出的一种设定分割线方法：

$\begin{array}{l} \sin ψ_{m}^{cut} = 0 & (e \geq) e_{max} \end{array}$ 。

(4)

这时为避免膨胀角的突然变化,对膨胀角进行平滑处理：

$\sin ψ_{m}^{smooth} = 100 (1 - e / e_{\max}) \sin ψ_{m} (e \geq 0.99 e_{\max})$ 。

(5)

对有截止因子f_cut（Cut-off）,光滑处理（Smoothing）和没有截止因子f_cut（No Cut-off）3种情况下的土体膨胀规律进行了如图2所示的对比,可以看出,光滑处理（Smoothing）情况下的数值计算最为稳定（渐变过程）,与实际情况也最为贴近。

图 2 标准三轴压缩试验的体积应变曲线

Figure 2. Curves volumetric strain of standard triaxial compression tests

下载: 全尺寸图片幻灯片

PH模型的体积屈服函数定义为

$f^{v} = g^{v} = {\tilde{q}}^{2} / a^{2} + p^{2} - p_{c}^{2} = 0$ 。

(6)

式中 $α$ 为一个常量,是基于数值固结仪测试推导得到的； $\tilde{q}$ 为一个剪应力量度；p为硬化参数,而硬化参数初值 $p_{c}^{}$ （代表是先前固结压力）可以由初始应力状态和超固结比OCR来确定。此处的体积硬化采用相关联流动法则。

PH-SS模型中的小应变特性是由小应变参数 $G_{0}$ 和 $γ_{70}$ 来反映的,其相应的关系式为

$G / G_{0} = {(1 + 0.385 γ / γ_{70})}^{^{- 2}}$ ,

(7)

式中 $γ_{70}$ 为剪切模量衰减到初始剪切模量的70%所对应的剪应变；初始剪切模量G₀可以由参考压力p^ref下的模量 $G_{0}^{ref}$ 在极小的应变水平下得

$G_{0} = G_{0}^{ref} {(\frac{c^{'} \cos ϕ^{'} + {σ^{'}}_{3} \sin ϕ^{″}}{c^{'} \cos ϕ^{'} + p^{ref} \sin ϕ^{'}})}^{m}$ 。

(8)

考虑到土体的实际情况,有限差分程序在PH-SS模型中充分考虑了卸载剪切模量 $G_{u r}$ 和卸载弹性模量 $E_{u r}$ 之间的关系,并设置了一个下限值：

$G_{u r} = E_{u r} / [3 (1 - 2 v)]$ 。

(9)

1.2 基于三轴模拟试验的PH-SS模型分析

利用FLAC^3D构建三轴压缩数值计算模型,模型尺寸为1 m×1 m×1 m,围压为100 Pa。为说明PH-SS模型在反映软土加载和卸载过程的优越性,该三轴模拟试验利用不同的本构模型（MC模型,PH模型和PH-SS模型）进行计算,且每种情况下都采用相同的循环荷载过程。

图3为三轴试样偏应力随轴向应变的变化曲线。可以看出,MC模型的卸载模量和加载模量是相同的；PH模型能够反映加载和卸载模量的差异性,但是不区分加载路径和卸载路径；PH-SS模型不仅能够区分加载和卸载模量的差异,而且能较好地区分加载路径和卸载路径,这种情况最符合相关试验^[15]的结果,同时也和谢东武等^[16]的模拟结果保持一致。PH-SS模型中出现回滞环进一步说明,随着剪应变的增大,土体的剪切模量迅速衰减,在加载—卸载的过程中产生塑性积累。

图 3 三轴试验中加卸载的应力应变曲线

Figure 3. Stress-strain curves of loading and unloading in triaxial tests

下载: 全尺寸图片幻灯片

图4为三轴试样的剪切模量随剪应变的变化曲线。由此可知：MC模型和PH模型中的模量始终保持恒定；而PH-SS模型能够较好地反应土体模量随剪应变增大逐渐衰减的特性,并达到一个下限值,这与式（9）是保持一致的,说明了PH-SS模型能够很好地表现土体的小应变特性。

图 4 剪切模量随剪应变的关系曲线

Figure 4. Curve of shear modulus with shear strain

下载: 全尺寸图片幻灯片

2. 基于不同模型的盾构隧道地层变形分析

2.1 计算模型的建立

根据上海某隧道^[17]的实际工程情况,建立二维数值模型。其中,隧道埋深为11.9 m,隧道直径为D为6.2 m。为避免边界的影响,X两侧方向取6倍的直径尺寸,下底边界取3倍的直径尺寸,即模型尺寸为76 m×40 m（宽×高）。隧道周围取1倍的隧道半径进行加密处理,单元尺寸为0.5 m,其他单元尺寸为1 m,数值模型共计有11942个节点、5880个单元。

本文分别利用PH-SS模型与MC模型对土体的力学行为进行模拟。参照文献[17～20],确定上海地区典型土层的小应变模型参数见表1。为保证数值计算结果更具有对比性,MC模型参数选取与PH-SS模型中相同的数值（其中,弹性模量取对应土层的加卸载模量）。衬砌材料利用shell单元来模拟,衬砌结构为线弹性材料, $γ$ =24.5 kN·m^-3, $ν$ =0.2；E= 24.44×10³ MPa,衬砌厚度d=0.35m。

表 1 PH-SS模型的土体参数^[17-20]

Table 1. Parameters of soil of PH-SS model

介质	$γ$ /(kN·m^-3)	$ν$	φ/(°)	c/kPa	/MPa	$E_{ur}^{ref}$ /MPa	/MPa	$p_{}^{ref}$ /kPa	m	$γ_{70}$ /10^-4
1素填土	18.0	0.35	18.0	10.0	4.0	23.2	92.8	100	0.9	2
2粉质黏土	18.7	0.35	15.0	18.0	6.0	34.8	139.2	100	0.9	2
3淤泥质粉质黏土	18.0	0.35	24.0	10.0	5.0	29.0	116.0	100	0.9	2
4淤泥质黏土	17.0	0.35	19.5	10.0	4.0	23.2	92.8	100	0.9	2
5-1粉质黏土	18.3	0.35	18.5	12.0	5.0	29.0	116.0	100	0.9	2
5-2粉质黏土	18.3	0.35	21.0	15.0	6.5	37.7	150.8	100	0.9	2
6粉质黏土	19.9	0.35	17.0	30.0	6.4	37.1	148.4	100	0.9	2
7-粉质黏土夹粉砂	20.1	0.30	26.3	6.0	10.2	59.2	236.6	100	0.9	2

下载: 导出CSV

| 显示表格

2.2 计算方案

采用应力释放法模拟盾构隧道的开挖和支护过程^[21]。利用不同的应力释放系数（ $λ$ 为0.1,0.2,0.3,0.4）反映隧道周围的地应力因开挖所引起的扰动程度。本文分别采用PH-SS和MC两种本构模型,总计八组计算工况,如表2。

表 2 计算工况设计

Table 2. Design of numerical cases

工况名称	是否考虑小应变特性	本构模型	应力释放系数λ
MC-0.1	否	MC	0.1
MC-0.2			0.2
MC-0.3			0.3
MC-0.4			0.4
PH-SS-0.1	是	PH-SS	0.1
PH-SS-0.2			0.2
PH-SS-0.3			0.3
PH-SS-0.4			0.4

下载: 导出CSV

| 显示表格

为研究深部土体的水平位移情况,取5个剖面（编号为^#1,^#2,^#3,^#4和^#5,相对位置分别如图5所示）为研究对象。对于地层竖向位移情况,本文重点关注地表、地表以下5 m处、地表以下10 m处、以及隧道以下5 m处的地层隆沉规律。

图 5 深层位移分析剖面位置示意图

Figure 5. Profile locations for displacement analysis of deep soils layers

下载: 全尺寸图片幻灯片

2.3 地表位移分析

图6给出了利用PH-SS模型计算得到的地表位移曲线,图7是利用MC模型计算得到的地表位移曲线。

图 6 利用PH-SS模型计算得到的地表位移曲线

Figure 6. Curves of surface displacement based on PH-SS model

下载: 全尺寸图片幻灯片

图 7 利用MC模型计算得到的地表位移曲线

Figure 7. Curves of surface displacement based on MC model

下载: 全尺寸图片幻灯片

由图6,7可知,两种地表沉降规律大致相似,但是在较大应力释放系数情况下（如 $λ$ =0.4）,PH-SS模型下的最大地表位移明显比MC模型下的地表位移大很多；在同一种本构模型计算工况下,应力释放系数越大,地表变形（竖向沉降、水平位移）越大,这是因为应力释放系数较大时,地层的应力扰动较大,造成地表变形较大。由图8可以得到,当应力释放系数小于0.3时,采用MC模型和PH-SS模型计算得到的地层损失率和最大地表沉降差别相对较小；但是当应力释放系数大于0.3时,采用PH-SS模型比MC模型计算的结果大得多,此时考虑模型之间的差异是必要的,这是因为盾构开挖应力释放系数较小（地层扰动较小）时,地层变形的小应变行为表现得不明显,当应力释放系数较大（地层扰动较大）时,小应变特征则十分明显,此时因剪应变而引起的刚度变化则十分显著。

图 8 地表最大沉降值、地层损失率与应力释放系数之间的关系

Figure 8. Relationship among maximum surface settlement, volume loss rate and stress release coefficient

下载: 全尺寸图片幻灯片

2.4 深部土层位移分析

以应力释放系数为0.3为例进行两种本构模型下的深部土层位移分析。图9分别为地表以下5 m、地表以下10 m的土层沉降的对比曲线,可以得到,深部土层位移场中的考虑小应变特性的沉降槽宽度比不考虑小应变特性的情况要小一些。图10是剖面^#1,^#2,^#4,^#5剖面的水平位移曲线,鉴于^#3剖面处于中央的位置,水平位移的差异性不明显,所以不单独列出^#3剖面的曲线。对图10的结果进行分析,沿着隧道对称的剖面（^#1和^#4,^#2和^#5）的水平位移变化规律沿隧道轴线对称,当埋深较浅时,由于地层沉降槽的形成产生了趋向隧道轴线方向的水平位移,而当埋深达到开挖面附近时,由于盾构管片的挤压作用,水平位移变成正值（土体沿着坐标轴正向移动）,随着分析剖面与隧道轴向距离的增加,岩土体受开挖扰动的影响逐渐减小。由此可见,盾构推进引起隧道周围土体向开挖面临空移动,最终水平位移大致关于中轴线对称。同时,考虑到隧道开挖造成周围地层在不同方向都产生了卸载作用,因而考虑土体加载模量和卸载模量的差异性时,土体模量变化较大,PH-SS模型相较MC模型在隧道轴线埋深处的地层水平位移更大。

图 9 地表以下地层剖面的沉降曲线

Figure 9. Settlement curves of soil layers below surface

下载: 全尺寸图片幻灯片

图 10 深部土体剖面的水平位移曲线

Figure 10. Curves of horizontal displacement of deep soil profile

下载: 全尺寸图片幻灯片

2.5 隧道以下地层隆起分析

盾构开挖会引起隧道以下地层发生卸载,PH-SS模型能区分土体加载和卸载模量的差异性。图11给出了不同应力释放系数情况下的隧道以下地层的位移规律。可以看出,两种本构模型下的地层位移规律差异较大。图12给出了对应地层处的最大隆起量（隆起为正、沉降为负）变化规律。可以看出,PH-SS模型计算隆起量比MC模型计算值较小,并且当应力释放系数增大时,其差异性逐步增大,说明应力释放系数越大（地层扰动越大）,应力相关模量的力学响应越明显,这也与前述保持一致,说明了PH-SS模型在地层隆起研究中的适宜性。

图 11 不同应力释放系数下的隧道以下5 m地层处的隆沉规律

Figure 11. Heave or settlement curves of layer 5 m below tunnel under different stress release coefficients

下载: 全尺寸图片幻灯片

图 12 地层隆起量随应力释放系数的变化曲线

Figure 12. Curves of soil layer uplift with stress release coefficient

下载: 全尺寸图片幻灯片

3. 与实测数据对比

根据上述研究,地层变形受到应力释放系数的影响较大。而实际工程情况受很多因素的影响,本文从简单问题出发,以地层损失为衡量标准^[22],通过三维数值计算结果反演二维模型的应力释放系数。三维模拟按照现场工况进行,每步推进1个管片环宽即1.0 m,利用强度不同的等代层来模拟注浆层及超挖层部分,最终确定二维数值计算的应力释放系数为0.38。因此,将 $λ$ =0.38情况下基于PH-SS模型和MC模型的数值计算结果与文献[17]中的监测数据进行对比,如图13所示。

图 13 数值计算结果与实测数据对比

Figure 13. Comparison between numerical results and measured data

下载: 全尺寸图片幻灯片

由此可知,利用PH-SS模型计算的地表沉降值更大,其地表沉降曲线相较MC模型更加符合实际监测结果,因此验证了预测盾构隧道地层变形中考虑小应变特性的必要性,同时也说明了PH-SS模型在预测地表沉降变形中的适用性。

4. 结论

软土地区盾构施工中的大部分区域处于小应变的状态,本文基于FLAC^3D内置的MC模型和PH-SS模型对盾构施工地层变形规律展开研究。首先基于三轴模拟试验证明了PH-SS在表现土体小应变特性方面的合理性和优越性；依托上海某软土隧道工程,分别利用MC模型和PH-SS模型进行二维盾构施工数值计算研究,着重研究了地层变形尤其是地表沉降随应力释放系数的演化规律,以及利用两种本构模型计算所得结果之间的差异性；最后将数值结果与现场数据进行了对比。

（1）相较MC模型,PH-SS模型不仅可以区分软土加载和卸载模量的差异性,而且可以较好地反映土体剪切模量随剪应变增大而逐渐衰减的特性,在表现软土行为时更为合理。

（2）当应力释放系数小于0.3时,预测地表沉降时考虑小应变特性与不考虑时的情况区别不大；但是当应力释放系数大于0.3时,考虑土体小应变特性是十分必要的。

（3）应力释放系数越大,土体应力相关模量的力学响应越明显,地层变形越大。

（4）数值计算结果与实测数据进行对比的结果表明,考虑小应变特性的PH-SS模型地表沉降更大,与实测数据更为贴近,说明了PH-SS模型在预测盾构隧道地表变形规律更为合理。

图 1 管涌出口处不同的粒子类型

Figure 1. Different particle types at downstream outlet

下载: 全尺寸图片幻灯片

图 2 向后侵蚀管涌计算方法的流程图

Figure 2. Algorithm procedure of backward erosion piping

下载: 全尺寸图片幻灯片

图 3 数值模型的几何示意图

Figure 3. Geometric conditions of numerical model

下载: 全尺寸图片幻灯片

图 4 计算结束时的管涌通道和总水头分布

Figure 4. Piping channel and total head contours at the end of simulation

下载: 全尺寸图片幻灯片

图 5 不同网格尺寸和粒子数量下压力水头的分布

Figure 5. Pressure water heads under different mesh sizes and particle numbers

下载: 全尺寸图片幻灯片

图 6 含防渗墙砂土地基的几何示意图

Figure 6. Geometric conditions of sandy soil foundation with cutoff wall

下载: 全尺寸图片幻灯片

图 7 不同上游水头下管涌通道的长度

Figure 7. Lengths of piping channel under different water heads without cut-off wall

下载: 全尺寸图片幻灯片

图 8 含防渗墙地基中管涌通道和总水头分布

Figure 8. Contours of piping channels and total water head of foundation with cut-off wall

下载: 全尺寸图片幻灯片

图 9 含防渗墙地基顶部边界上的压力水头

Figure 9. Pressure water heads on top boundary of foundation with cut-off wall

下载: 全尺寸图片幻灯片

参考文献(30)

[1]	VAN BEEK V M, BEZUIJEN A, SELLMEIJER J B, et al. Initiation of backward erosion piping in uniform sands[J]. Géotechnique, 2014, 64(12): 927-941. doi: 10.1680/geot.13.P.210
[2]	ROBBINS B A, VAN BEEK V M, LÓPEZ-SOTO J F, et al. A novel laboratory test for backward erosion piping[J]. International Journal of Physical Modelling in Geotechnics, 2018, 18(5): 266-279. doi: 10.1680/jphmg.17.00016
[3]	VANDENBOER K, DOLPHEN L, BEZUIJEN A. Backward erosion piping through vertically layered soils[J]. European Journal of Environmental and Civil Engineering, 2019, 23(11): 1404-1412. doi: 10.1080/19648189.2017.1373708
[4]	ROBBINS B A, MONTALVO-BARTOLOMEI A M, GRIFFITHS D V. Analyses of backward erosion progression rates from small-scale flume experiments[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2020, 146(9): 04020093. doi: 10.1061/(ASCE)GT.1943-5606.0002338
[5]	VANDENBOER K, VAN BEEK V M, BEZUIJEN A. 3D character of backward erosion piping[J]. Géotechnique, 2018, 68(1): 86-90. doi: 10.1680/jgeot.16.P.091
[6]	POL J C, KANNING W, BEEK V M, et al. Temporal evolution of backward erosion piping in small-scale experiments[J]. Acta Geotechnica, 2022, 17(10): 4555-4576. doi: 10.1007/s11440-022-01545-1
[7]	AKRAMI S, BEZUIJEN A, VAN BEEK V, et al. Analysis of development and depth of backward erosion pipes in the presence of a coarse sand barrier[J]. Acta Geotechnica, 2021, 16(2): 381-397. doi: 10.1007/s11440-020-01053-0
[8]	VANDENBOER K, VAN BEEK V M, BEZUIJEN A. Analysis of the pipe depth development in small-scale backward erosion piping experiments[J]. Acta Geotechnica, 2019, 14(2): 477-486. doi: 10.1007/s11440-018-0667-0
[9]	VAN BEEK V M, VAN ESSEN H M, VANDENBOER K, et al. Developments in modelling of backward erosion piping[J]. Géotechnique, 2015, 65(9): 740-754. doi: 10.1680/geot.14.P.119
[10]	BRIAUD J L, GOVINDASAMY A V, SHAFII I. Erosion charts for selected geomaterials[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2017, 143(10): 04017072. doi: 10.1061/(ASCE)GT.1943-5606.0001771
[11]	周晓杰, 丁留谦, 姚秋玲, 等. 悬挂式防渗墙控制堤基渗透变形发展模型试验[J]. 水力发电学报, 2007, 26(2): 54-59. doi: 10.3969/j.issn.1003-1243.2007.02.011 ZHOU Xiaojie, DING Liuqian, YAO Qiuling, et al. Laboratory model test for evolution of seepage deformation controlled by means of suspended cut-off wall in foundation of dike[J]. Journal of Hydroelectric Engineering, 2007, 26(2): 54-59. (in Chinese) doi: 10.3969/j.issn.1003-1243.2007.02.011
[12]	ZHOU X J, JIE Y X, LI G X. Numerical simulation of the developing course of piping[J]. Computers and Geotechnics, 2012, 44: 104-108. doi: 10.1016/j.compgeo.2012.03.010
[13]	FASCETTI A, OSKAY C. Dual random lattice modeling of backward erosion piping[J]. Computers and Geotechnics, 2019, 105: 265-276. doi: 10.1016/j.compgeo.2018.08.018
[14]	ROBBINS B A, GRIFFITHS D V. A two-dimensional, adaptive finite element approach for simulation of backward erosion piping[J]. Computers and Geotechnics, 2021, 129: 103820. doi: 10.1016/j.compgeo.2020.103820
[15]	LIANG Y, YEH T C J, WANG J J, et al. An auto-adaptive moving mesh method for the numerical simulation of piping erosion[J]. Computers and Geotechnics, 2017, 82: 237-248. doi: 10.1016/j.compgeo.2016.10.011
[16]	WANG Y A, NI X D. Hydro-mechanical analysis of piping erosion based on similarity criterion at micro-level by PFC^3D[J]. European Journal of Environmental and Civil Engineering, 2013, 17(S1): 187-204.
[17]	TRAN D K, PRIME N, FROIIO F, et al. Numerical modelling of backward front propagation in piping erosion by DEM-LBM coupling[J]. European Journal of Environmental and Civil Engineering, 2017, 21(7/8): 960-987.
[18]	FROIIO F, CALLARI C, ROTUNNO A F. A numerical experiment of backward erosion piping: kinematics and micromechanics[J]. Meccanica, 2019, 54(14): 2099-2117 doi: 10.1007/s11012-019-01071-7
[19]	RAHIMI M, SHAFIEEZADEH A. Coupled backward erosion piping and slope instability performance model for levees[J]. Transportation Geotechnics, 2020, 24: 100394. doi: 10.1016/j.trgeo.2020.100394
[20]	ROTUNNO A F, CALLARI C, FROIIO F. A finite element method for localized erosion in porous media with applications to backward piping in levees[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2019, 43(1): 293-316. doi: 10.1002/nag.2864
[21]	FUJISAWA K, MURAKAMI A, NISHIMURA S I. Numerical analysis of the erosion and the transport of fine particles within soils leading to the piping phenomenon[J]. Soils and Foundations, 2010, 50(4): 471-482. doi: 10.3208/sandf.50.471
[22]	ZHANG X S, WONG H, LEO C J, et al. A thermodynamics-based model on the internal erosion of earth structures[J]. Geotechnical and Geological Engineering, 2013, 31(2): 479-492. doi: 10.1007/s10706-012-9600-8
[23]	WEWER M, AGUILAR-LÓPEZ J P, KOK M, et al. A transient backward erosion piping model based on laminar flow transport equations[J]. Computers and Geotechnics, 2021, 132: 103992. doi: 10.1016/j.compgeo.2020.103992
[24]	LEI X Q, HE S M, CHEN X Q, et al. A generalized interpolation material point method for modelling coupled seepage-erosion-deformation process within unsaturated soils[J]. Advances in Water Resources, 2020, 141: 103578. doi: 10.1016/j.advwatres.2020.103578
[25]	LIANG D F, ZHAO X Y, SOGA K. Simulation of overtopping and seepage induced dike failure using two-point MPM[J]. Soils and Foundations, 2020, 60(4): 978-988. doi: 10.1016/j.sandf.2020.06.004
[26]	CECCATO F, YERRO A, GIRARDI V, et al. Two-phase dynamic MPM formulation for unsaturated soil[J]. Computers and Geotechnics, 2021, 129: 103876. doi: 10.1016/j.compgeo.2020.103876
[27]	王兆南, 王刚. 饱和孔隙介质的耦合物质点-特征有限元方法[J]. 岩土工程学报, 2023, 45(5): 1094-1102. doi: 10.11779/CJGE20220332 WANG Zhaonan, WANG Gang. Coupled material point method and characteristic finite element method for saturated porous media[J]. Chinese Journal of Geotechnical Engineering, 2023, 45(5): 1094-1102. (in Chinese) doi: 10.11779/CJGE20220332
[28]	ZIENKIEWICZ O C, TAYLOR R L, NITHIARASU P. The Finite Element Method for Fluid Dynamics (Seventh Edition)[M]. Oxford: Butterworth-Heinemann, 2014.
[29]	ABAQUS G. Abaqus 6.11[M]. Providence: Dassault Systemes Simulia Corporation, 2011.
[30]	ROBBINS B A, VAN BEEK V M, POL J C, et al. Errors in finite element analysis of backward erosion piping[J]. Geomechanics for Energy and the Environment, 2022, 31: 100331. doi: 10.1016/j.gete.2022.100331