0.   引言

向后侵蚀管涌(backward erosion piping)是每年汛期常发生于土石坝地基中的一种渗透破坏形式，它通常形成于透水砂层与不透水黏土层的接触面处。当下游的黏土层因外力扰动形成集中的出水口时，其附近的砂土在集中渗流作用下会被侵蚀带走，并在出水口外部形成砂沸现象，最终堆积为一个锥形沙丘。当上下游水头差足够大时，冲刷通道将沿着强/弱透水层的界面处从下游出水口向上游不断贯通，从而导致坝体结构的破坏。近30 a来，学术界针对向后侵蚀管涌现象进行了大量的室内试验研究，研究内容主要涉及了管涌的启动机制^[1-3]、演化机制^[4-9]以及颗粒运移机制^[10]等方面。

目前，关于向后侵蚀管涌的机理性研究已经较为成熟，但多数结论都是基于小尺度的模型试验获得的，对于实际工程问题而言，仍然有待发展出一个准确、高效的数值分析方法。目前，相关的数值方法大致可被概括为以下3种：①基于传统的渗流方程，采用合理地增大通道内渗透系数的方式分析管涌对于流场的影响^[11-15]。该方法的优点在于简单易用，通常可以结合泥沙运输法则^[13-14]来反映管道内颗粒的侵蚀作用。为了使计算更加精确，Robbins等^[14]、Liang等^[15]使用自适应网格技术加密管涌通道内的计算网格，使管涌区域的数值精确性得到了极大的提升。但是这一方法的缺点在于仍然采用达西定律描述管道内的流体运动，与真实管道内的流体运动存在着较大差异。②离散元和计算流体力学的耦合求解方法^[16-18]。离散元法在模拟颗粒在管道内的侵蚀、搬运、沉积等问题上有着极大的优势，加深了研究者对于管涌微观行为的认识。但是，受困于离散元法在计算效率上所需的巨大开销，使得该方法难以应用于实际工程。③将被侵蚀的颗粒、土骨架、孔隙水分为不同的三相，并分别建立相应的颗粒侵蚀质量守恒方程^[19-23]。这种方法的优势在于能够考虑管涌随时间的演化过程，但是缺点是需要定义土颗粒侵蚀的本构法则，且相关的计算参数难以标定。

近年来，物质点法(material point method，MPM)作为一种能够求解大变形问题的数值方法而受到了岩土力学界的广泛关注。其中，也不乏一些学者将MPM用于土体渗流破坏的计算。Lei等^[24]通过在非饱和MPM中考虑细颗粒的迁移方程，计算了边坡降雨入渗条件下因细粒流失而产生的失稳问题。Liang等^[25]使用两套物质点分别描述流体和土体，模拟了不同防渗类型堤坝的漫顶和渗透破坏。Ceccato等^[26]模拟了库水位升降后堤坝的大变形失稳破坏问题。但是，以上研究中并未涉及后退型管涌的计算，且存在着流体压力不稳定、压力边界难以施加的困难^[27]。

基于笔者所提出的耦合物质点-特征有限元方法^[27](coupled material point-characteristic finite element method，下文简称“耦合MPM-FEM”)，提出了一种适用于计算向后侵蚀管涌问题的新算法。新算法使用局部水力梯度临界值作为管涌触发的判别条件，对水力梯度大于临界值的物质点进行删除，能够同时描述待侵蚀的土体和管涌通道内的流体。通过数值算例初步地证明了算法在向后侵蚀管涌问题中的适用性。

1.   耦合物质点-特征有限元方法

1.1   控制方程

假设流体和土体颗粒具有不可压缩性，则饱和孔隙介质的连续性方程可记为^[27]

式中：ϕ为孔隙率；u_i为速度；上标s和f分别为固体和流体。当ϕ = 1时，式(1)退化为不可压流体的质量守恒方程。

固体和流体的动量守恒方程分别为

式中：${\bar \rho ^{\text{s}}}$，${\bar \rho ^{\text{f}}}$分别为孔隙介质内固体和流体的平均密度，${\bar \rho ^{\text{s}}} = (1 - \phi ){\rho ^{\text{s}}}$，${\bar \rho ^{\text{f}}} = \phi {\rho ^{\text{f}}}$；b_i为外力加速度，${\sigma '_{ij}}$为土骨架之间的有效应力；${\tau _{ij}}$为流体黏性力，${\tau _{ij}} =$ ${\mu ^{\rm{f}}}(\partial u_i^{\rm{f}}/\partial {x_j} + \partial u_j^{\rm{f}}/\partial {x_i})$；μ^f为动力黏度；p^f为流体压力。

由于固体和流体分别采用了拉格朗日视角和欧拉视角进行描述，因此固体加速度${{{D^{\rm{s}}}u_i^{{s}}} \mathord{\left/ {\vphantom {{{D^{\rm{s}}}u_i^{{s}}} {Dt}}} \right. } {Dt}}$是采用物质导数形式表示的，而流体的加速度中则考虑了对流项$u_j^{\text{f}}{{\partial u_i^{\text{f}}} \mathord{\left/ {\vphantom {{\partial u_i^{\text{f}}} {\partial {x_j}}}} \right. } {\partial {x_j}}}$的影响。D_i表示固体和流体之间动量交换的拖曳力项，采用达西定律表示：

式中：k为渗透系数；g为重力加速度。

当ϕ = 1时，系统由饱和孔隙介质变为了单相的纯流体。此时，式(3)中的D_i需要被省略，与式(1)相结合便构成了不可压黏性流体的Navier-Stokes(N-S)方程组。

1.2   特征分裂方法

欧拉视角下，有限元在处理流体动量的方程时首要面对的问题便是如何处理非线性对流项，若直接使用标准伽辽金法推导式(3)的弱形式，会导致计算产生数值震荡。使用特征线法^[28]对式(3)进行自伴随处理，处理后的方程消除了加速度中的对流项，经局部泰勒展开后可以获得相应的时间离散格式：

式中：Δt为时间增量。

耦合MPM-FEM方法中，使用了分裂算法^[28](Fractional step/Chorin’s projection)对流体速度和压力进行交替计算，式(5)可以被分裂成两部分：

式中：$u_i^{{\text{f}} * }$为流体的中间辅助速度。

同理，式(2)也可以被分裂成两部分：

将式(7)，(9)代入式(1)中即可获得流体的压力泊松方程：

式(6)~(10)均可以采用标准的伽辽金法推导相应的弱形式。

在推导空间数值格式时，式(8)，(9)使用MPM离散，在物质点上进行数值积分；式(6)，(7)，(10)则采用FEM进行离散，并在高斯点上进行数值积分。离散后的方程遵循着经典的分裂算法求解过程，即“速度预测-压力求解-速度修正”：①计算式(6)，(8)得到流体和固体的辅助速度$u_i^{{\rm{f}} * }$和$u_i^{s * }$；②求解式(10)得到第n+1时步的压力p^f；③基于式(7)，(9)对流体和固体的速度进行修正，得到n+1时步速度${}^{n + 1}u_i^{\rm{f}}$和${}^{n + 1}u_i^{\rm{s}}$。对于孔隙介质的水力耦合大变形问题而言，在每个时步内完成上述求解过程后还需对物质点的变形梯度、应力、空间位置、孔隙率等变量进行更新，之后进入下一时步的计算。整个过程中物质点使用欧拉有限元网格进行控制方程的求解。相较于多相的弱可压MPM，耦合MPM-FEM具备了压力计算稳定、边界易于施加的特点，详细的算法描述可参考文献[27]。

对于本文所关注的后退型管涌问题而言，主要侧重于管涌发生后土体内部的渗流场计算，因而固体物质点是被完全固定的，并忽略了其应力的更新。其原因有主要有：①土颗粒全程在管涌通道内发生侵蚀与运移，土体所受的应力状态对其影响很小；②若考虑土体应力的计算，则在选取本构模型及其相应模型参数上存在困难。

2.   向后侵蚀管涌的计算方法

管涌的触发需要定义一个临界侵蚀条件，该临界条件可以由流速^[12]、水力梯度^[13-14]或其它参数来度量，用以判别管涌向上游的发展。这里采用判断粒子处局部水力梯度临界值的方式作为管涌触发的判别准则：

式中：i_crit为局部水力梯度的临界值；i^p为物质点处的水力梯度，可以通过所属单元节点的总水头插值得到。下标I为背景网格的节点，${N_I}({\boldsymbol x^{\rm{p}}})$为节点对于粒子的形函数，h_I为节点上的总水头，

式中：$p_I^{\rm{f}}/({\rho ^{\rm{f}}}g)$，${z_I}$分别为节点的压力水头和位置水头。

将渗流场内的物质点分为待删除粒子P_delete、过渡粒子P_transition和一般粒子P_general 3种类型。其中，P_delete表示因流体侵蚀作用而在渗流方向上保持临空，且满足管涌触发水力条件的粒子；P_transition表示在未来具有被侵蚀倾向的粒子；P_general则属于具有稳定性质的一般粒子。图 1给出了计算开始前管涌出口附近3种粒子的示意图。在后退型管涌发生的初始阶段，渗流出口处的高水力梯度会引起过高的流速，使得出口处的粒子被水流搬运走。因此，当出口处的P_delete粒子满足管涌的启动准则时，将被从求解域内删除掉。每一个P_delete粒子均含有一个支持域，在其支持域范围内的粒子则被视为P_transition。其余的物质点被定义为一般粒子P_general。为了方便搜索P_delete和P_transition，可将P_delete所在的格子以及周围格子定义为标记单元，从而达到减少搜索计算量的目的。

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

Figure  1.  Different particle types at downstream outlet

下载: 全尺寸图片 幻灯片

图 2给出了向后侵蚀管涌计算方法的流程图。

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

Figure  2.  Algorithm procedure of backward erosion piping

下载: 全尺寸图片 幻灯片

详细的计算步骤可分为以下5步：①在计算开始之前，将管涌出口处的粒子标定为待删除粒子⁽⁰⁾P_delete，并且给定粒子的支持域半径R_P、局部水力梯度临界值i_crit、以及压力e_p和速度e_u的误差容许值。②采用耦合MPM-FEM计算求解域内的渗流场，直至速度和压力均收敛至稳定。这一过程中物质点被完全固定，不更新其应力、速度、空间位置等变量。③将节点上的总水头插值至所有P_delete粒子上，对P_delete粒子所在、及其临近单元进行标记，作为接下来的查找范围。④根据式(11)判断P_delete粒子是否满足管涌的触发条件。一旦成立，则在标记单元中寻找P_delete粒子支持域内的一般粒子P_general，将它们转换为过渡粒子P_transition。⑤删除当前时步的⁽ⁿ⁾P_delete粒子，并将P_transition粒子转换为新的⁽ⁿ⁺¹⁾P_delete粒子，最后进入第下一时步的循环。其中，P_delete粒子的支持域半径R_P和平均粒子间距相一致。

上述的求解过程类似于ABAQUS^[29]中的“生死单元”方法。当满足管涌触发条件的粒子被删除后，其所在位置处的孔隙率即变为ϕ = 1，而后的计算则相当于在求解纯流体的N-S方程，因此可根据孔隙率的分布来判别管涌发生的位置。

3.   算法的网格敏感性

为验证所提出方法在向后侵蚀管涌模拟中的网格敏感性，选取了Robbins等^[30]所建立的数值模型作为相应的验证算例。如图 3所示，二维平面内存在一个长0.48 m，高0.1 m的矩形砂土试样。试样的左侧边界通过施加水头作为上游迎水侧，顶部边界存在一个直径为0.02 m的出水口，出水口中心距离左侧边界的长度为0.35 m，剩余边界均为不透水层。参考Robbins等^[30]所选取的计算参数：渗透系数k = 3.0×10^-4 m/s、孔隙率ϕ = 0.4、流体密度ρ^f = 1000 kg/m³、动力黏度μ^f =10^-3 Pa·s、局部水力梯度的临界值i_crit = 0.43。上游入水口边界施加0.052 m的压力水头，出水口处的压力水头为零。

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

Figure  3.  Geometric conditions of numerical model

下载: 全尺寸图片 幻灯片

图 4给出了向后侵蚀管涌停止发展后的计算结果，其中图 4(a)，(b)分别是管涌通道(孔隙率)和总水头的分布情况。从图 4可以看出，管涌通道尖端处的总水头产生了局部的集中现象。图 5给出了不同网格尺寸和物质点数量下顶部边界压力水头的分布情况。图 5中，Δx为网格的间距，PT/cell为单个网格内初始物质点的数量。计算稳定时管涌尖端大约位于左侧边界0.16 m处，尖端左侧的压力水头从左至右快速地下降，尖端右侧管道内压力水头的下降逐渐放缓。此外，随着背景网格和初始物质点数量的加密，计算结果总体上表现出与试验结果不断接近的趋势。

图  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

下载: 全尺寸图片 幻灯片

4.   含防渗墙地基的向后侵蚀管涌模拟

如图 6所示，周晓杰等^[11]，Zhou等^[12]构建了长2.3 m、宽0.6 m、高0.8 m的砂土试样，用以进行向后侵蚀管涌的模型试验。模型的右侧和底部为不透水边界，在距离左侧边界1.45 m处的上表面设置了一个直径为0.04 m的圆形出水口，水头恒定为零；地基的左侧为全断面进水边界，水头逐级地被加高，从而促使下游出水口处的流速不断增大，以诱发管涌的形成。该试验分为无防渗墙和有防渗墙两类，其中防渗墙被放置于距上游1.1 m处，采用一个高度为0.06 m的有机玻璃板表示。四边形背景网格的间距为0.02 m，每个网格在初始时放置4个物质点。为了计算简便，采用将砂土渗透系数减小1000倍的方式表示防渗墙，且代表防渗墙的粒子不能够被删掉。砂土的渗透系数k = 8.0×10^-3 cm/s，孔隙率ϕ = 0.35，流体密度ρ^f = 1000 kg/m³、动力黏度μ^f = 10^-3 Pa·s。

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

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

下载: 全尺寸图片 幻灯片

无防渗墙保护的条件下，管涌通道随着上游水头的增加而不断地推进至迎水侧，图 7给出了无防渗墙时不同水头下管涌通道长度的计算结果(局部临界水力梯度i_crit = 0.21)。图 7中，在上游水头施加至0.165 m之前，本文计算的管涌通道长度均未发生改变；当水头增加至0.235 m后，管涌通道快速地贯通至上游进水边界，此时的破坏水头为0.255 m。Zhou等^[12]的计算值与本文计算的管涌发展趋势基本相同，而试验结果则表现出更加平缓的递进式发展过程。

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

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

下载: 全尺寸图片 幻灯片

图 8给出了含防渗墙地基中管涌通道和总水头的分布情况(局部临界水力梯度i_crit = 0.31)。图中，当管涌通道在沿顶面不透水层向上游发展一段距离后，管道尖端在抵达防渗墙时会转向防渗墙底部发展，而相应的水头也会在防渗墙底部出现集中的现象，这与试验的观测结果^[11-12]相吻合。图 9给出了顶部边界上的压力水头分布。可以看出，无论是计算值还是试验值均会在管道尖端处因水头集中而形成一个快速下降的过程。计算值与Zhou等^[12]的计算结果基本相吻合，但是两者与试验值在管道尖端左侧却存在着较大的差别，试验中的压力水头要显著高于计算值。

图  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

下载: 全尺寸图片 幻灯片

5.   结语

基于饱和孔隙介质的耦合MPM-FEM，发展了一个适用于模拟向后侵蚀管涌问题的数值算法。该算法将渗流场中的物质点分别定义为待删除粒子、过渡粒子、一般粒子3种类型，通过比较待删除粒子处水力梯度与其临界值的大小，对满足管涌判别法则的粒子进行删除，从而达到模拟管涌发展过程的目的。算法的网格敏感性通过一个小尺度的模型试验进行了验证，其结果表明算法的精确度随着与网格尺度和物质点数量的细化而不断提高。最后，通过一个中尺度模型试验证明了算法模拟含防渗墙地基管涌破坏的能力。