多材料介质中无压渗流问题的离散元求解

The solution of unconfined seepage problem in multi-material medium using a discrete element method

摘要

为了更好地模拟多材料介质无压渗流问题，提出了一种基于颗粒离散元的离散颗粒渗流模型（discrete seepage element model，DSEM）。该模型将连续介质离散化为密集堆积的颗粒系统，赋予颗粒流体性质，并通过颗粒间压力差模拟渗流过程。DSEM模型能够自动计算自由面位置和水头分布，无须预设自由面位置。为了验证DSEM模型的有效性，首先，模拟了土坝渗流问题，模拟结果与室内模型试验结果一致。随后，模拟了双材料矩形坝和分层土体渗流问题，并与其他数值方法进行比较，证明了DSEM在多种材料问题中的有效性。最后，将DSEM应用于分区梯形坝的渗流问题，成功计算了不同渗透系数下的自由面位置。研究结果表明，DSEM模型适用于多材料介质的无压渗流问题，为工程渗流问题提供了一种新的有效计算方法。

Abstract

To better simulate the unconfined seepage problems in multi-material media, A discrete seepage element model (DSEM) based on the discrete element method for particles was proposed. This model discretizes the continuous medium into a densely packed system of particles, endows the particles with fluid properties, and simulates the seepage process through the pressure difference between particles. The DSEM model can automatically calculate the position of the free surface and the distribution of the water head without the need to pre-set the position of the free surface. To verify the effectiveness of the DSEM model, the seepage problem of an earth dam was first simulated, and the simulation results were consistent with the results of laboratory experiments. Subsequently, the seepage problems of a double-material rectangular dam and a stratified soil were simulated, and when compared with other numerical methods, the effectiveness of the DSEM in problems involving multi-materials was demonstrated. Finally, the DSEM was applied to the seepage problem of a zoned trapezoidal dam, successfully calculating the positions of the free surface under different permeability coefficients. The research results indicate that the DSEM model is suitable for solving seepage problems in multi-material media and provides a new and effective computational method for solution of seepage problems in engineering.

Graphical abstract

关键词

离散元 / 无压渗流 / 自由面 / 多层材料介质

Key words

discrete element method / unconfined seepage / free surface / multi-material medium

地下水流动对地质构造稳定性、地下水资源管理以及环境保护等方面都具有重要影响，研究地下水渗流问题一直是岩土工程和地质科学领域中的一个重要研究课题。考虑到自然土体的成层特点，以及坝体的功能分区设计原则，多种材料介质中无压渗流问题是经常遇到的情形，也是分析的难点。因此，对多种材料介质中无压渗流行为进行深入研究具有重要的理论和实际意义。

数值分析方法具有成本低、计算准确的特点，近年来在岩土工程和地质领域得到了广泛的应用^[1]。目前常用于渗流问题的数值分析方法主要为有限元法^[2-3]、边界元法^[4]、有限差分法^[5-6]、有限体积法^[7]和数值流形法^[8-10]以及这些方法相互结合的混合方法^[11-13]。有限元法主要包括自适应网格法和固定网格法，自适应网格法^[14]需要在计算过程中不断更新网格，当自由面计算过程变化较大时，可能导致网格变形，而且在非均匀或者复杂结构中，自适应网格会改变介质边界，很难达到稳定解。固定网格法自Neuman提出用不变的网络Galperin法以来，在国内外学者的研究中，出现了剩余流量法、虚单元法、单元渗透矩阵调整法^[15]等方法。

近年来，计算机技术的不断发展为国内外学者提供了很多计算软件。目前计算无压渗流问题的软件主要有Midas GTSNX、ABAQUS、COMSOL等，这些软件大多基于有限元方法，在单元网格中迭代计算自由面。本文借鉴孔隙网络模型^[16]，发展了基于离散元(discrete element method, DEM)的DSEM。该方法在无压渗流问题的求解中具有能够考虑复杂结构和非均质性的优势。该方法适用于瞬态和稳态渗流的模拟，也适用于求解自由面的位置，边界条件设置较其他方法更简单，不需要额外考虑网络的问题，具有精度高、数值求解稳定性好的特点。与其他方法相比，对于不连续的复合材料和结构，模型也不需要进行额外的调整。

本文采用离散颗粒渗流模型，主要计算了包含多层材料介质的渗流自由面问题。通过与室内模型试验结果和其他数值方法结果进行比较，验证了该模型的有效性。通过本研究，期望为地质工程中多种材料介质的渗流问题提供新的数值模拟方法，并为未来相关研究提供有益的参考。

1 离散颗粒渗流模型

1.1 模型假设

DSEM是在DEM的基础上，将原本存属于岩土体颗粒缝隙之间液体的性质赋予离散颗粒，岩土体中的渗流通过颗粒接触对进行传递，具体假设如下：

1) 将离散元模型中的颗粒赋予流体的性质。如图1(a)所示，工程中的连续多孔介质模型被简化为不同大小密集排布的颗粒，流体在颗粒间孔隙的流动被简化到颗粒自身流体性质的传递，颗粒的流体性质（孔隙压力P、渗流量q等）通过颗粒间的接触或者颗粒与边界的接触对进行传递，图1(b)中展示模型当中的两种基本接触对，即颗粒-颗粒接触对和颗粒-边界接触对。

2) 颗粒间的流体传递条件。相邻颗粒之间存在孔隙压力差，并且流体传递的方向总是与孔隙压力降的方向相同。

1.2 模型边界条件

在离散渗流单元模型中，模拟无压渗流问题时，不考虑基质吸力，自由面定义为耦合颗粒孔隙水压力为0的分界线。

图2为自由面渗流问题的示意图。在模型经过渗流迭代计算达到稳定后，提取孔隙水压力为0的面便是自由面，不需要额外调整设定。

在离散渗流单元模型中，需要设置两种边界条件，即恒定水头边界和出渗水头边界。由于不透水面不参与渗流的传递，因此不需要额外设置。恒定水头边界条件为

\phi =\overline{\phi }(1)

出渗水头边界条件为

\phi =y(2)

式中：\phi为计算域内任意一点的总水头；\overline{\phi }为已知水头边界值；y为计算点的垂直坐标。式(2)表明出渗面上孔隙压力为零，总水头等于位置水头。

图2中的边界，对于恒定水头边界，AB边界采用下游水头边界条件{\phi }_{AB}=y_{\mathrm{1}}，CD边界采用上游水头边界条件{\phi }_{CD}=y_{\mathrm{2}}。E为未知的出渗点，不需要人为设定。只需要将DF设为出渗面，出渗水头边界DF上条件为{\phi }_{DF}=y_{DF}。

值得注意的是，不同于之前的自适应网格法和固定网格法，本文中自由面的求解方法不需要额外设置和调整自由面，只要保证初始条件正确，模型在计算完成之后，即可得到设定条件下的自由面分布。

1.3 流体传递方式

在DSEM中，每个颗粒可能会因为来自不同周围颗粒和周围边界的流体流入和流出，设t时刻后，经过一个时步∆t颗粒i的水头为{\phi }_i^{t+∆t}，计算公式如下：

{\phi }_i^{t+∆t}={\phi }_i^t+∆{\phi }_p^t+∆{\phi }_w^t(3)

式中：{\phi }_i^t为t时刻的水头；∆{\phi }_p^t为t时刻周围的颗粒p对颗粒i水头的增量；∆{\phi }_w^t为t时刻周围的边界w对颗粒i水头的增量。

流体在连续介质中的流动遵循达西定律，可知流体渗流量与水力梯度的关系为

q=kS\frac{∆\phi }{L}(4)

式中：q为流体流量；k为渗透率；S为渗流方向截面面积；∆\phi为总水头损失；L为渗流路径长度。在渗流问题中，\phi为总水头，是位置水头和压力水头之和。

\phi =y+\frac{P}{\gamma }(5)

式中：y为参考基准面以上的高度，一般取模型纵坐标的最小值为参考基准面；P为孔隙压力；\gamma为流体的容重。

设流体增量∆V，流入给定体积V_{\mathrm{0}}中，体积V_{\mathrm{0}}内的孔隙压力变化可由式(6)确定：

∆P=k_v\frac{∆V}{V_{\mathrm{0}}}(6)

式中：k_v为流体的体积模量。

在模型中，当流体从颗粒i流向颗粒j时，假想一个渗流通道，如图3(a)所示，渗流通道是图中红色虚线框矩形，渗流通道的长度L_{ij}取颗粒i与颗粒j圆心之间的距离，渗流通道的截面可以粗略地看作两个颗粒直径的平均值S_{ij}=R_i+R_j，其中R_i和R_j分别为颗粒i和j的半径。通道的渗透系数k_{ij}定义为两个颗粒渗透系数的最小值，即k_{ij}=\mathrm{m}\mathrm{i}\mathrm{n}(k_i,k_j)，其中k_i和k_j分别为颗粒i和颗粒j的渗透系数，这种处理可以方便地考虑多种介质材料问题。

需要注意的是，渗流通道的面积是人为定义的，与真实渗流通道的面积有一定的差异，因此需要引入修正系数\alpha对其修正。\alpha对瞬态渗流至关重要，但本文探讨的是自由面渗流，属于稳态渗流，\alpha的修正可以忽略。由于修正主要影响渗流到达稳定的时长，而不会影响自由面的最终位置，因此本文统一取\alpha =\mathrm{1}。

同样地，当流体从边界w流向颗粒i时，如图3(b)所示，渗流通道的长度L_{iw}取颗粒i到边界的垂直距离，渗流通道的截面可以取颗粒i的直径S_{iw}=\mathrm{2}{R}_i。通道的渗透系数k_{iw}定义为两个颗粒渗透系数的最小值，即k_{iw}=\mathrm{m}\mathrm{i}\mathrm{n}(k_i,k_w)，其中k_w为边界w的渗透系数。

因此，颗粒i与周围n个颗粒之间的流体交换量V_{j\to i}为

V_{j\to i}={\displaystyle \sum _{j=\mathrm{1}}^n}\frac{\alpha k_{ij}\left({\phi }_j-{\phi }_i\right)S_{ij}∆t}{L_{ij}}(7)

流体交换量所引起颗粒i孔隙压力的变化量∆P_i为

∆P_i=k_v\frac{V_{j\to i}}{V_i}=\frac{\alpha k_v∆t}{V_i}{\displaystyle \sum _{j=\mathrm{1}}^n}\frac{k_{ij}{S}_{ij}}{L_{ij}}\left({\phi }_j-{\phi }_i\right)(8)

式中：{\phi }_i和{\phi }_j分别为颗粒i和颗粒j的总水头；\alpha为修正系数，用于修正渗流通道截面面积。

由于颗粒的位置不发生变化，所以颗粒所处位置总水头的变化量只取决于压力水头，即

{\phi }_i^{t+∆t}-{\phi }_i^t=\frac{\mathrm{1}}{\gamma }\left(P_i^{t+∆t}-P_i^t\right)(9)

式中：{\phi }_i^{t+∆t}、P_i^{t+∆t}分别表示在t时刻经过时步∆t后的颗粒i的总水头值和孔隙压力值。颗粒i因周围颗粒流体流入导致水头的变化量为

∆{\phi }_p^t=\frac{\alpha k_v∆t}{\gamma V_i}{\displaystyle \sum _{j=\mathrm{1}}^n}\frac{k_{ij}{S}_{ij}}{L_{ij}}\left({\phi }_j^t-{\phi }_i^t\right)(10)

类似地，可以得到颗粒i因相邻的m个边界流入导致水头的变化量为

∆{\phi }_w^t=\frac{\alpha k_v∆t}{\gamma V_i}{\displaystyle \sum _{w=\mathrm{1}}^m}\frac{k_{iw}{S}_{iw}}{L_{iw}}\left({\phi }_w^t-{\phi }_i^t\right)(11)

综上，颗粒在模型当中更新水头的计算公式为

\begin{array}{l}{\phi }_i^{t+∆t}={\phi }_i^t+\frac{\alpha k_v∆t}{\gamma V_i}\left[{\displaystyle \sum _{j=\mathrm{1}}^n}\frac{k_{ij}{S}_{ij}}{L_{ij}}\left({\phi }_j^t-{\phi }_i^t\right)\right.\\ \left.+{\displaystyle \sum _{w=\mathrm{1}}^m}\frac{k_{iw}{S}_{iw}}{L_{iw}}\left({\phi }_w^t-{\phi }_i^t\right)\right]\end{array}(12)

1.4 计算流程

算例的计算流程如下：

1) 建立离散渗流颗粒模型及赋予材料参数。首先采用动态平衡的方式生成颗粒堆积模型，然后赋予颗粒材料参数，包括颗粒密度\rho、渗透系数k、体积模量k_v等。另外，根据渗流边界条件，设置水头边界和出渗水头边界。

2) 进行颗粒与颗粒接触对检索、颗粒与边界接触对检索，生成虚拟渗流通道。

3) 分别判断每个接触对（两个耦合颗粒或者耦合颗粒与边界）之间是否存在孔隙压力差以及水头差。

4) 在颗粒之间和颗粒与边界之间分别进行流体流量计算、孔隙压力变化量的计算以及总水头的变化量的计算。

5) 根据步骤4)计算的结果更新耦合颗粒位置处的孔隙压力以及总水头。

6) 等待计算达到稳定（相邻两步的水头差小于阈值），绘制自由面以及进行其他后续处理。

2 模型验证

多材料介质渗流是一个复杂的非线性过程，受多种因素的影响。在进行多材料介质渗流模拟之前，需要先验证所提离散颗粒渗流模拟方法的有效性。相较于采用解析解或其他数值解作为参考解，更倾向于采用试验作为参考解，用以验证所提方法是否能应用于工程问题。为此，选取两个土石坝试验作为参考解，对所提方法进行验证。下面依次介绍模型试验结果和数值模拟结果，并分析所提方法的有效性。

2.1 试验模型

Al-Janabi等^[17]设计了两个土坝试验模型，图4^[17]为土坝渗流试验的模型，放置在长2.4 m、宽0.5 m、高0.5 m的透明玻璃箱内。图5^[17]为两个模型的截面图，两个模型的排水条件不同，一个是棱体排水，另一个是褥垫排水。试验中上游水位保持0.3 m，坝体模型保持目标水位9 h，以达到稳定渗流状态。为了观察潜水线测量水位高度，在大坝下游斜坡上钻了4个孔，得到了4个不同点位的自由面位置。

2.2 模拟结果

根据试验建立数值模拟计算模型，模型a、b的颗粒数分别为41 276和38 251，上游设置恒定水头边界，水头为0.3 m；排水边界设置恒定水头边界，水头为0 m。经过计算得到结果如图6所示。

图6中，黑色点为试验测得的4个自由面位置。表1给出了模型a试验测得自由面位置和数值计算结果以及相对误差。

结果表明，数值模拟计算结果与试验结果基本吻合，证明所提渗流算法是有效的。由于试验所提供的有效数据有限，本文后续算例中，将模型与其他数值模拟结果进一步对比分析。

3 算例

上一节验证了所提渗流模拟方法针对单种材料渗流问题的有效性，接下来将进一步模拟多层材料介质的渗流问题。首先，3.1节模拟非均质矩形坝的渗流问题，该问题是多材料介质渗流问题的经典案例，被很多文献用来验证所提方法对复杂渗流问题的鲁棒性；其次，3.2节模拟了分层材料的渗流问题，并验证本文绘制的等水头线的正确性；最后，3.3节模拟分区梯形坝无压渗流问题，验证所提方法对于工程问题的适用性。

3.1 非均质矩形坝无压渗流问题

对于不连续材料的工程渗流问题，为了测试模型处理自由面剧烈变化的能力，按照计算模型进行测试，如图7所示。

图中非均质矩形坝体是在不透水基础上由两种材料组成的坝，尺寸为\mathrm{5}\mathrm{m}\times \mathrm{10}\mathrm{m}，横向分为两个区域，两个区域的渗透系数之比为1∶10，左侧区域的渗透系数为k_{\mathrm{1}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}\mathrm{m}/\mathrm{s}，右侧区域的渗透系数为k_{\mathrm{2}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{m}/\mathrm{s}，上游总水头为10 m，下游总水头为2 m。AB为上游恒定水头边界，{\phi }_{AB}=\mathrm{10}\mathrm{m}；CG为下游恒定水头边界，{\phi }_{CG}=\mathrm{2}\mathrm{m}；GD为出渗水头边界，{\phi }_{GD}=y_{GD}。

从图8中可以看出，采用DSEM计算得到的自由面与Zheng等^[18]和Yang等^[19]得到的结果比较一致，说明DSEM可以处理非均质材料中的无压渗流问题。一般方法中，如果不对材料界面进行特殊处理，许多传统方法难以实现。Wang等^[20]在界面处加入了额外的场节点，来提高数值解的稳定性。

图9为模型的流网图，发现与文献[18]计算得到的流网基本一致，可以证实该方法在求解这类问题上的准确性。

3.2 多层材料介质问题

为了验证本文中的数值方法求解多层材料介质流网的准确性，采用了Tracy等^[21]设计提出的计算模型，计算模型如图10所示。

图中左边、底边和右上边蓝色实线为不透水边界，不同层的材料用虚线隔开，其他边界为恒定水头边界和出渗水头边界，颗粒数为27 503，图中边界AE上，{\phi }_{AE}=\mathrm{22.86}\mathrm{m}，边界CD上，{\phi }_{CD}=\mathrm{26.67}\mathrm{m}。从下到上各层材料的渗透系数分别为k_{\mathrm{1}}=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{m}/\mathrm{s}，k_{\mathrm{2}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{m}/\mathrm{s}，k_{\mathrm{3}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{m}/\mathrm{s}。

图11为计算得到的水头分布图。对比图12、图13和图14，可以看出等水头线与不排水边界基本相互垂直，并且与文献[21]和文献[22]计算结果基本一致，证明了该方法的可行性。

3.3 分区梯形坝无压渗流问题

为了验证模型求解无压渗流问题的正确性，建立计算模型，如图15所示。

图15中坝体分为3个区域，区域1的渗透系数为k_{\mathrm{1}}，区域2的渗透系数为k_{\mathrm{2}}，考虑分区的影响，将两个区域的渗透系数分两种情况设计，情况一：k_{\mathrm{1}}=k_{\mathrm{2}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{m}/\mathrm{s}；情况二：k_{\mathrm{1}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{m}/\mathrm{s}，k_{\mathrm{2}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}\mathrm{m}/\mathrm{s}。上游AB为恒定水头边界，{\phi }_{AB}=\mathrm{10}\mathrm{m}；下游DE为出渗水头边界，{\phi }_{DE}=y_{DE}；CD为排水边界，{\phi }_{CD}=\mathrm{0}。模型颗粒数为28 254。

图16和图17分别为坝体情况一和情况二时的水头分布情况和自由面位置，图16(a)着色部分为离散渗流单元模型的水头分布计算结果，图16(b)黑色流网为文献[22]计算得到的结果，从图中可以看出，自由面的位置基本一致。对比分析图16和图17，可以看到k_{\mathrm{1}}=k_{\mathrm{2}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{m}/\mathrm{s}时，整个自由面是连续光滑的，而当k_{\mathrm{1}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{m}/\mathrm{s}，k_{\mathrm{2}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}\mathrm{m}/\mathrm{s}时，自由面在右侧界面处发生了明显的转折。这是由于该界面两侧的渗透系数发生了突变，该转折类似于3.1节非均质矩形坝算例。这个算例充分说明了心墙（区域2）对渗流的阻碍作用，体现了心墙的防渗效果。该算例进一步验证了所提渗流算法对解决复杂多材料介质渗流问题的有效性。

4 结论

本文在颗粒离散元方法的基础上借鉴孔隙网络模型，拓展了二维颗粒渗流离散元模型(DSEM)，将渗流的流体性质赋予模型中的紧密堆积的离散颗粒，简化了计算模型。该方法简化了边界条件的设置，只需要设置恒定水头边界和出渗水头边界即可完成对边界条件的设定，自由面不需要进行特殊设定，计算稳定后可以很容易得出自由面的位置，通过多个算例验证了所提DSEM的有效性。首先是以两个土石坝渗流模型试验为参考，验证了DSEM求解均质材料问题的正确性。接着将DSEM应用于多个多材料介质渗流模拟中，包括非均质矩形坝无压渗流问题、多层材料介质的渗流问题，以及分区梯形坝无压渗流问题。模拟的结果与其他数值方法高度一致，验证了DSEM对于多材料介质渗流问题的准确性和可靠性。模拟算例还表明，该方法在多材料介质渗流问题中具有简单、统一、适应性强等优点。在模型的设置和计算过程中，无须预设自由面，不同材料边界处也无须进行特殊处理，具有较强的通用性。这些特点对于处理各种复杂的材料组合和边界条件问题非常有益，可避免人为设定对结果的影响。因此，所提的DSEM有望为复杂渗流相关的工程问题提供有效计算方法。

本文主要研究了二维多材料介质渗流问题，而现实问题往往是三维的，因此，后续有必要研发基于离散元的三维渗流模型。另一方面，本文主要考虑纯无压渗流问题，无法进行渗透破坏模拟。研发鲁棒的流固耦合模块，并进一步进行渗透破坏模拟，是接下来的重点工作。

参考文献

[1] 熊祎滢, 刘春, 刘恒, 等. 水力耦合围岩数值模拟现状分析[J]. 四川建筑, 2023, 43(4): 188-189, 197.

[2] 何朝葵, 速宝玉, 盛金昌, 等. 基于弱有限元法的渗流分析[J]. 岩土工程学报, 2023, 45(7): 1526-1532.

[3] 宁皓, 张向阳, 牟宗琪. 有限元分析在尾矿库渗流及抗滑稳定工程实践中的应用[J]. 岩土工程技术, 2023, 37(3): 357-360.

[4] 方思冬, 程林松, 饶翔, 等. 一种改进格林元方法及在渗流问题中的应用[J]. 计算力学学报, 2021, 38(6): 787-795.

[5] 王睿, 周宏伟, 卓壮, 等. 非饱和土空间分数阶渗流模型的有限差分方法研究[J]. 岩土工程学报, 2020, 42(9): 1759-1764.

[6] 朱帅润, 何博, 吴礼舟, 等. 非饱和渗流模拟中非均匀空间网格的改进方法[J]. 水文地质工程地质, 2023, 50(1): 32-40.

[7] 宿晓辉, 张明亮, 贺英芝. 基于有限体积法的三维稳态渗流计算方法及其应用[J]. 人民长江, 2022, 53(12): 198-203.

[8] 陈远强, 杨永涛, 郑宏, 等. 饱和-非饱和渗流的数值流形法研究与应用[J]. 岩土工程学报, 2019, 41(2): 338-347.

[9] 陈远强, 郑宏, 屈新. 基于数值流形法的降雨入渗与坡面径流耦合算法研究[J]. 应用数学和力学, 2023, 44(12): 1499-1511.

[10] 周凌峰, 王媛, 冯迪. 求解非均质渗流场的改进数值流形方法[J]. 岩土工程学报, 2021, 43(7): 1288-1296, 1377-1379.

[11] 李卓徽, 周宗青, 高成路, 等. 基于近场动力学与有限体积法耦合的裂隙岩体渗流模拟[J]. 同济大学学报(自然科学版), 2022, 50(9): 1251-1263.

[12] 张道生. 工程岩体应力—渗流耦合的PD-FEM-FVM模拟分析方法及应用[D]. 济南: 山东大学, 2022.

[13] 邸元, 吴大卫, WU Y S. 油藏渗流-应力耦合分析的FEM-FVM混合方法的改进[J]. 岩石力学与工程学报, 2020, 39(S1): 2645-2654.

[14] 王敏. 非均质油藏高效数值模拟算法的研究[D]. 合肥: 中国科学技术大学, 2019.

[15] 秦甜甜, 丁国辉, 程勤波. 自由渗流面及渗流量推求的数值计算方法探讨[J]. 地下水, 2018, 40(4): 12-14.

[16] BLUNT M J. Flow in porous media-pore-network models and multiphase flow[J]. Current Opinion in Colloid & Interface Science, 2001, 6(3): 197-207.

[17] AL-JANABI A M S, GHAZALI A H, GHAZAW Y M, et al. Experimental and numerical analysis for Earth-Fill dam seepage[J]. Sustainability, 2020, 12(6): 14.

[18] ZHENG H, LIU F, LI C G. Primal mixed solution to unconfined seepage flow in porous media with numerical manifold method[J]. Applied Mathematical Modelling, 2015, 39(2): 794-808.

[19] YANG Y T, SUN G H, ZHENG H. Modeling unconfined seepage flow in soil-rock mixtures using the numerical manifold method[J]. Engineering Analysis With Boundary Elements, 2019, 108: 60-70.

[20] WANG B, LI J, JIANG Q, et al. An improved FE-Meshfree method for solving steady seepage problems[J]. Computers and Geotechnics, 2020, 119: 103223.

[21] TRACY F T, RADHAKRISHNAN N. Automatic generation of seepage flow nets by finite element method[J]. Journal of Computing in Civil Engineering, 1989, 3(3): 268-284.

[22] 李伟, 郑宏. 基于数值流形法的渗流问题边界处理新方法[J]. 岩土工程学报, 2017, 39(10): 1867-1873.

基金资助

天津市津南区“揭榜挂帅”科技计划项目(2023JB03)