砂砾土潜蚀演化及宏微观力学响应的CFD-DEM研究

CFD-DEM Study on Suffusion Evolution and Macro-Micro Mechanical Responses of Sandy Gravel

摘要

库水周期性涨落极易引发消落带内部碎石土的细颗粒流失，进而导致岸坡深层变形甚至失稳，其内在机制极其复杂。针对这一渗流破坏敏感介质，采用计算流体力学-离散元耦合（CFD-DEM）方法，构建3种不同初始细颗粒含量的间断级配砂砾土模型，开展多水力梯度条件下的渗流潜蚀与三轴剪切全过程模拟试验。基于宏微观指标演变规律的对比，重点量化揭示了孔隙比对土体宏观力学性能劣化的控制机制。结果表明：累计侵蚀率与水力梯度呈正相关，与初始细颗粒含量呈负相关；细颗粒流失导致强力链断裂与重组，孔隙比与渗透系数永久性增大，土体微观结构向亚稳态演化；潜蚀诱发的孔隙结构损伤导致土体宏观强度衰减与剪胀弱化，且峰值内摩擦角与孔隙比呈显著的线性负相关；而各试样的残余内摩擦角受潜蚀影响不显著。研究揭示了间断级配土体潜蚀进程与强度劣化的宏微观关联的特征，可为库岸边坡的灾变防控提供理论支撑。

Abstract

Periodic reservoir water-level fluctuations easily trigger the loss of fine particles within the gravelly soil of the hydro-fluctuation belt, which can induce deep-seated deformation or even instability of bank slopes, yet the underlying mechanism is extremely complex. Focusing on this seepage-sensitive medium, this study employs a coupled computational fluid dynamics-discrete element method (CFD-DEM) approach. Models of gap-graded sandy gravel with three distinct initial fine particle contents were constructed to conduct whole-process simulations of seepage-induced suffusion and triaxial shear under varying hydraulic gradients. Based on the comparison of macro- and micro-scopic index evolution, this study specifically quantifies and reveals the controlling mechanism of the void ratio on the deterioration of the soil’s macroscopic mechanical properties. The results indicate that the cumulative erosion ratio is positively correlated with the hydraulic gradient but negatively correlated with the initial fine particle content. The loss of fine particles leads to the rupture and reorganization of strong force chains, causing a permanent increase in void ratio and permeability, thereby driving the soil microstructure toward a metastable state. Suffusion-induced pore structure damage results in soil strength degradation and dilatancy weakening, and the peak internal friction angle exhibits a significant linear negative correlation with the void ratio. However, the residual internal friction angle of various specimens is not significantly affected by suffusion. This study elucidates the macro-micro correlations between the suffusion process and strength degradation in gap-graded soils, providing theoretical support for the disaster prevention and control of reservoir bank slopes.

Graphical abstract

关键词

间断级配砂砾土 / 潜蚀 / 力学性能劣化 / 细颗粒含量 / 水力梯度 / CFD-DEM

Key words

gap-graded sandy gravel / suffusion / mechanical property deterioration / fine particle content / hydraulic gradient / CFD-DEM

0　 引 言

地下水渗流引发潜蚀对土体物理力学性能影响的研究，前期多集中在土石坝和堤防等工程^[1-2]，随后逐渐向堆积层边坡、滑坡堰塞坝、深厚覆盖层等天然级配不良土体方面发展^[3-4]。三峡库区消落带土体在30 m水位落差的周期性涨落驱动下，物理状态发生变化、力学性能持续降低^[5-6]。其中，巫山段多见由崩塌堆积形成的碎石土质岸坡，秭归段多见颗粒呈棱角状的花岗岩风化砂质岸坡，因其土体颗粒级配的不连续性与内部结构的不稳定性，往往成为岸坡渗流破坏的优势通道与内在薄弱区^[7]。因此，深入探索潜蚀的演化进程，揭示其对土体物理力学性能劣化的作用机制，对于库岸边坡的长期稳定性评价具有重要的理论价值。

潜蚀的发生与演化受几何条件、力学状态与水力条件的共同控制，三者相互关联、彼此影响^[8-9]。几何条件决定了土体的潜在可侵蚀性，基于颗粒级配曲线形态与粒径比等参数，研究者已提出多种判别准则，用于评估土体发生潜蚀的敏感性^[10-12]。对于间断级配土而言，细颗粒作为粗颗粒骨架孔隙的填充物，其含量直接决定了土体的骨架结构，进而影响土体的抗侵蚀性阈值^[13-15]。在水力驱动方面，临界水力梯度是控制潜蚀启动的关键指标，其数值受应力状态^[16]、渗流方向^[17]及加载历史^[18]的显著影响。宏观力学响应层面，大量室内试验已证实，潜蚀会导致土体抗剪强度、内摩擦角及变形模量等宏观参数的显著劣化^[19-21]。例如，Ouyang等^[22]定量描述了侵蚀程度对土体应力-应变关系的弱化效应。然而，受限于室内试验的观测尺度，传统手段难以捕捉潜蚀过程中颗粒骨架的动态重构及力链网络的演化特征。

近年来，计算流体力学-离散元耦合方法（Computational Fluid Dynamics-Discrete Element Method，CFD-DEM）因能有效捕捉颗粒-流体相互作用及微观结构演化，被广泛应用于岩土多相流研究^[23-25]。已有学者利用该方法探讨了粒径比^[26]、颗粒形态^[27]及相对密实度^[28]对潜蚀过程的影响，并将该方法拓展应用于边坡失稳^[29]、大坝渗漏^[30]、隧洞涌水^[31]等工程场景分析。然而，针对“渗流潜蚀-结构演化-力学劣化”全过程的连续性耦合分析尚有待加强。同时，针对间断级配砂砾土这一特定介质，其从微观结构损伤到宏观强度衰减的跨尺度定量关联机制，以及潜蚀造成的孔隙比改变如何定量控制强度衰减的内在规律，仍需进一步解析。

鉴于此，本文聚焦岸坡典型的间断级配砂砾土，基于CFD-DEM耦合框架建立了考虑不同初始细颗粒含量（FC）的离散元模型。通过在多水力梯度下开展渗流潜蚀与三轴剪切模拟，系统揭示了侵蚀率、力链网络及渗透系数等关键参数的演变规律，阐明了细颗粒含量及水力梯度对潜蚀发展过程的耦合作用机制。此外，通过对比分析潜蚀前后土体的孔隙比、峰值偏应力、体积应变及内摩擦角等宏微观参量的变化，定量揭示了细颗粒流失对土体物理力学性能的劣化机制，并构建了孔隙比与峰值内摩擦角的定量关系模型。

1　 计算流体力学-离散元耦合方法

本文采用的CFD-DEM耦合框架基于3个开源C++库构建^[32]，包括用于在DEM域中模拟颗粒运动的LIGGGHTS、用于在CFD域中求解流体压力与速度场的OpenFOAM以及用于实现CFD与DEM双向数据交互与耦合计算的CFDEMcoupling接口。

1.1　颗粒运动控制方程

在离散元框架中，土颗粒被视为独立的刚性单元，其平动与转动控制方程分别遵循牛顿第二定律和欧拉第二定律。为准确描述颗粒间的接触力学响应，本研究采用非线性Hertz-Mindlin模型模拟颗粒间的相互作用力。考虑到传统球形颗粒难以再现天然砂砾土不规则形状产生的咬合效应，本研究在接触模型中引入了基于弹塑性弹簧-阻尼模型^[33]的滚动阻力模型。该模型通过在接触点处施加额外的滚动力矩，表征非规则几何形状对颗粒转动的约束作用，从而更真实地反映砂砾土颗粒间的接触阻抗行为。颗粒i的运动控制方程为：

m_i\frac{\mathrm{d}{\mathit{v}}_i}{\mathrm{d}t}={\mathit{f}}_{\mathrm{p}\mathrm{f},i}+{\displaystyle \sum _{j=\mathrm{1}}^{k_c}}({\mathit{f}}_{c,ij}+{\mathit{f}}_{d,ij})+m_i\mathit{g}（1）

I_i\frac{\mathrm{d}{\mathit{\omega }}_i}{\mathrm{d}t}={\displaystyle \sum _{j=\mathrm{1}}^{k_c}}({\mathit{M}}_{t,ij}+{\mathit{M}}_{r,ij})（2）

其中：m_i和I_i分别为颗粒i的质量和惯性矩；{\mathit{v}}_i和{\mathit{\omega }}_i分别为颗粒i的平移速度和角速度；k_{\mathrm{c}}为与颗粒i发生接触的颗粒总数；{\mathit{f}}_{\mathrm{p}\mathrm{f},i}、{\mathit{f}}_{\mathrm{c},ij}和{\mathit{f}}_{\mathrm{d},ij}分别为颗粒-流体相互作用力、颗粒间弹力和黏性阻尼力； g 为重力加速度；{\mathit{M}}_{\mathrm{t},ij}和{\mathit{M}}_{\mathrm{r},ij}分别为机械扭矩和黏性阻尼扭矩。

1.2　流体相控制方程

本文采用未解析的CFD-DEM耦合方法。由于流体网格尺寸大于颗粒粒径，流体相运动通过求解连续性方程与基于局部平均的Navier-Stokes方程描述。此外，模型中引入局部颗粒体积分数与流体-颗粒相互作用力，以实现流-固两相间的动量耦合，具体方程为：

\frac{\partial {\alpha }_{\mathrm{f}}}{\partial t}+\nabla \cdot ({\alpha }_{\mathrm{f}}{\mathit{U}}_{\mathrm{f}})=\mathrm{0}（3）

\frac{\partial ({\alpha }_{\mathrm{f}}{\rho }_{{}_{\mathrm{f}}}{\mathit{U}}_{{}_{\mathrm{f}}})}{\partial t}+\nabla \cdot ({\alpha }_{{}_{\mathrm{f}}}{\rho }_{{}_{\mathrm{f}}}{\mathit{U}}_{{}_{\mathrm{f}}}{\mathit{U}}_{{}_{\mathrm{f}}})=-{\alpha }_{{}_{\mathrm{f}}}\nabla \overline{p}-{\mathit{f}}_{\mathrm{p}\mathrm{f}}+\nabla \cdot ({\alpha }_{\mathrm{f}}{\mu }_{{}_{\mathrm{f}}}\nabla {\mathit{U}}_{{}_{\mathrm{f}}})+{\rho }_{{}_{\mathrm{f}}}{\alpha }_{{}_{\mathrm{f}}}\mathit{g}（4）

其中：{\alpha }_{\mathrm{f}}为流体的体积分数；{\mathit{U}}_{\mathrm{f}}为流体的矢量速度（{\mathit{U}}_{\mathrm{f}}{\mathit{U}}_{\mathrm{f}}为并矢）；\overline{p}为局部流体平均压力；{\mathit{f}}_{\mathrm{p}\mathrm{f}}为流体对颗粒的作用力；{\mu }_{\mathrm{f}}和{\rho }_{\mathrm{f}}分别为流体的动力黏度系数和流体密度。

1.3　流体-颗粒耦合方程

流体作用于颗粒的合力主要由拖曳力、压力梯度力、黏性力及浮力组成。鉴于本研究涉及的渗流过程雷诺数较低，虚拟质量力与巴塞特力等非定常项相对于拖曳力而言量级极小，故予以忽略^[34]。针对间断级配土体内孔隙结构非均匀性的特点，本文采用文献[35]的拖曳力模型。该模型通过孔隙率函数对拖曳力进行动态修正，能准确反映细颗粒在复杂孔隙空间内运移时所受的阻力作用，方程为：

\left\{\begin{array}{l}{\mathit{f}}_{\mathrm{d}}=\frac{\mathrm{1}}{\mathrm{8}}{C}_{\mathrm{d}}{\rho }_{\mathrm{f}}\mathrm{\pi }{d}_{\mathrm{p}}^{\mathrm{2}}({\mathit{U}}_{\mathrm{f}}-{\mathit{U}}_{\mathrm{p}})\left|{\mathit{U}}_{\mathrm{f}}-{\mathit{U}}_{\mathrm{p}}\right|n^{\mathrm{1}-\chi }\\ C_{\mathrm{d}}={(\mathrm{0.63}+\frac{\mathrm{4.8}}{\sqrt[]{\mathrm{R}{\mathrm{e}}_{\mathrm{p}}}})}^{\mathrm{2}}\\ \mathrm{R}{\mathrm{e}}_{\mathrm{p}}=\frac{n{\rho }_{\mathrm{f}}{d}_{\mathrm{p}}\left|{\mathit{U}}_{\mathrm{f}}-{\mathit{U}}_{\mathrm{p}}\right|}{{\mu }_{\mathrm{f}}}\\ \chi =\mathrm{3.7}-\mathrm{0.65}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{(\mathrm{1.5}-\mathrm{l}\mathrm{g}\mathrm{R}{\mathrm{e}}_{\mathrm{p}}{)}^{\mathrm{2}}}{\mathrm{2}}\right]\end{array}\right.（5）

其中：{\mathit{f}}_{\mathrm{d}}为流体对颗粒作用的拖曳力；C_{\mathrm{d}}为模型中单颗粒拖曳力系数；d_{\mathrm{p}}为颗粒的直径；{\mathit{U}}_{\mathrm{f}}和{\mathit{U}}_{\mathrm{p}}分别为流体和颗粒的速度；n为局部孔隙率；{\mu }_{\mathrm{f}}为流体的动力黏度系数；\mathrm{R}{\mathrm{e}}_{\mathrm{p}}为颗粒雷诺数；\chi为考虑局部孔隙率影响的经验修正指数。

2　 数值模拟方案

2.1　模型构建

间断级配砂砾土广泛存在于坝基深厚覆盖层和岸坡堆积体中。天然原状砂砾土的粒径跨度极广，若在CFD-DEM耦合计算中直接采用全级配模拟，计算单元数量将呈指数级增长，将导致计算效率显著降低且流-固耦合难以收敛。鉴于此，本文依据Wang等^[36]室内试验的原状级配曲线特征，提取骨架与填充颗粒的主导粒径，并参照间断级配土微观模拟中的粒径简化方法^[26,37-38]，构建了由粗（12 mm）、细（2 mm）两组特征粒径组成的离散元试样。该模型设定的粗颗粒与细颗粒粒径比为6∶1，严格满足潜蚀发生的几何准则，这一设定有效复现了原状内部不稳定土^[36]易发生内部侵蚀的宏观物理本质。同时，该简化排除了次要过渡粒径的干扰，有助于在数值模拟中更纯粹、更精准地揭示“骨架-填充”结构在渗流作用下的演化机制。不同初始细颗粒含量（FC）土体试样的颗粒级配分布如图1所示。

根据Kezdi^[39]和Kenney、Lau^[1]提出的判定准则，所构建试样均属于内部不稳定土，具备发生潜蚀的潜在几何条件。为确保数值模型与物理试验在材料参数上的一致性，采用与实测数据相同的颗粒密度（2 720 kg/m³）。试验过程中，须保持粗颗粒数量恒定，设置细颗粒质量分数分别为10%、20%和30%，该范围覆盖了从欠填充结构到过渡型结构的典型状态^[40]，并经压缩制样后得到初始孔隙比。参照潜蚀从孕育、发展到破坏所对应的不同水力环境^[12,19]，试验中分别施加0.8、1.6和3.2三种水力梯度，以模拟库区水位在低速、中速及高速变化条件下的差异化水力强度。试样经渗流侵蚀后，分别经历100、200、400 kPa三种围压等向固结，再进行轴向加载剪切。详细的数值模拟方案如表1所示。

DEM试样尺寸为90 mm×90 mm×135 mm, CFD计算域尺寸（90 mm×90 mm×225 mm）在渗流方向上显著大于DEM域，有利于形成均匀稳定的入流和出流边界条件，降低边界效应对流场的影响，同时确保土体处于饱和状态，如图2所示。在局部平均的CFD-DEM耦合方法中，为了保证流体求解的精度，Tao等^[37]和Hu等^[38]建议CFD网格宜设置为颗粒平均粒径的1~4倍。此外，时间步长的大小也直接影响模拟结果的计算效益与准确性。Li等^[41]建议DEM时间步长∆t_{\mathrm{D}\mathrm{E}\mathrm{M}}应满足临界瑞利时间步（式（6）），而CFD时间步长∆t_{\mathrm{C}\mathrm{F}\mathrm{D}}通常需满足库朗数（即流体质点在一个时间步长内运移距离与网格特征尺寸比值的无量纲数，用于评估流体求解的数值稳定性）小于1（式（7））。

（6）

（7）

其中：R_{\mathrm{m}\mathrm{i}\mathrm{n}}为最小颗粒的粒径；{\rho }_{\mathrm{p}}、v_{\mathrm{p}}、E分别为颗粒的密度、泊松比和弹性模量；Co为库朗数；{\Phi }_{\mathrm{m}}表示通过网格m的流体通量；V_{\mathrm{m}}为该网格的体积。

本文的颗粒群体接触刚度和密度远高于流体，颗粒时间步须明显小于流体时间步，设定DEM时间步长为5.0×10^-7 s，CFD时间步长为5.0×10^-5 s，流固耦合数据交换间隔设定为100个时间步，其他参数如表2所示。模拟试验过程主要包括5个连续步骤：生成颗粒、试样制备、渗流潜蚀、等向固结和三轴剪切，如图3所示。所有数值模拟均在搭配Intel Xeon E5-2696 V4处理器（22核心44线程）的工作站上完成，针对单一试样，模拟持续渗流30 s、剪切过程20 s，计算总耗时约为140 h。

2.2　模型验证

为进一步评估所建模型在模拟潜蚀过程中的可靠性，同时兼顾岸坡实际埋深环境的应力特征，本研究选取Wang等^[36]针对砾石土开展的大型室内渗透试验作为基准对照组。数值模拟中的颗粒密度（2 720 kg/m³）、围压（100 kPa）及水力条件均与室内试验保持一致，以最大限度提高模型模拟的准确度。考虑到室内试验与数值模拟在试样尺寸和渗流时间尺度上的客观差异，对数据进行了归一化处理。分别确定了颗粒累计侵蚀率（定义为累计流失细颗粒质量与初始细颗粒质量之比）及体积应变，并将数值模拟结果与试验数据进行比较，结果如图4所示。分析表明，数值模拟结果与试验值在关键的颗粒启动阶段吻合较好，但在较高的水力梯度区间存在一定偏差，推测是颗粒形状的接触效应以及水力梯度的加载方式差异所致。尽管如此，数值模拟结果所呈现的细颗粒侵蚀率和体积应变演化趋势与试验结果基本一致，表明所建试样模型和模拟方法在研究潜蚀行为方面具有良好的可靠性。

3　 潜蚀过程演化规律分析

3.1　累计侵蚀率

不同细颗粒含量试样，在不同水力梯度下潜蚀程度随渗流时间的变化规律如图5所示。随着渗流过程的推进，土体剩余颗粒数量呈现“先快速下降，随后下降速率逐渐减缓甚至趋于稳定”的趋势，见图5(a)。在相同水力梯度下，初始细颗粒含量较高的试样流失的颗粒总数更多。例如，在水力梯度为3.2时，FC30组试样累计流失细颗粒最多（28 814个），约为FC10组试样（12 137个）在同等条件下的2.4倍。分析图5(b)可知，各试样呈现的规律是“渗流初期侵蚀速率较快，随后逐渐减缓”，不同之处在于细颗粒含量越低且水力梯度越高，潜蚀趋稳的历时越短，而细颗粒含量较高或水力梯度较小时，计算时长内潜蚀并未完全稳定。在恒定水力梯度下，细颗粒含量越高，初始侵蚀速率越低，最终累计侵蚀率也越低，且影响效应相对显著。这一现象主要归因于细颗粒的填充效应：细颗粒含量越高，土体的初始孔隙比越小（见表1），细颗粒填充粗颗粒间孔隙的程度越高，从而形成更为密实的骨架结构，使土体渗透性降低，进而抑制细颗粒的迁移与流失。在相同细颗粒含量条件下，水力梯度成为控制侵蚀过程的关键因素，水力梯度越大，渗流前期侵蚀速率越快，最终累计侵蚀率也越高。其主要原因是，渗流力随水力梯度增大而增强，更大的流体驱动力使细颗粒更易从骨架中剥离，从而加剧了侵蚀发展。

不同工况下的最终累计侵蚀率如图6所示。总体而言，最终累计侵蚀率与水力梯度呈正相关，而与细颗粒含量呈负相关，但均不满足线性关系，说明潜蚀发展的内在机制较为复杂。

3.2　接触力链网络

潜蚀结束后试样的接触力链分布如图7所示。一方面，细颗粒含量对内部力链网络结构的分布产生显著影响。以水力梯度0.8为例，随着细颗粒含量从10%增至30%（图7(a)(b)(c)），试样中的力链总数明显增加，同时力链类型的比例却发生了根本性转变：在低细颗粒含量试样（FC10i0.8）中，承载结构主要由粗颗粒间接触形成的强力链主导，随着细颗粒含量的增加，强力链比例逐渐降低，弱力链（包括细-细和细-粗颗粒接触）则相应增多。其内在本质是：当细颗粒含量较低时，粗颗粒之间直接接触形成稳定的骨架结构，外部荷载主要通过强力链传递；而在细颗粒含量增加后，细颗粒逐渐填充粗颗粒间的孔隙，不仅导致弱力链数量增加，同时也阻碍了粗颗粒间的直接接触，致使强力链数量减少。另一方面，在相同细颗粒含量条件下，水力梯度的增大亦会引起力链网络的重构。以30%细颗粒含量试样为例，当水力梯度从0.8升高至3.2时（图7(c)(d)(e)），尽管力链总数有所减少，但强力链却显著增加。该现象说明，在较高水力梯度条件下，增强的渗流力优先破坏接触强度较低的弱力链，引发细颗粒迁移与重排，导致粗颗粒间直接接触比例上升，从而促使强力链数量增加、弱力链相应减少。

4　 土体物理力学性能演变规律

4.1　孔隙比的演化规律

孔隙比是连接微观结构损伤与宏观力学响应的核心状态变量。试样在潜蚀与三轴剪切过程中不同阶段的孔隙比变化规律如图8所示。结果显示，潜蚀结束后所有试样的孔隙比均显著上升，且初始细颗粒含量越高，孔隙比的绝对增量越大，FC30i3.2试样的孔隙比增量（∆e=0.117）约为FC10i3.2试样（∆e=0.062）的2倍。这是因为细颗粒流失会引起孔隙空间扩大，尽管高细颗粒含量试样的累计侵蚀率较低，但其在强渗流作用下流失的细颗粒数量更多（见图5(a)），因此表现出更大的孔隙比变化。此外，孔隙体积的增大，不仅直接改变了土体的物理状态，更为其后续的渗透性突变与力学强度劣化提供了条件。随后的等向固结过程中，孔隙比自潜蚀后的高点明显回落，但所有试样均未能恢复至潜蚀前的初始状态，说明由潜蚀引起的结构损伤具有不可逆性。其根本原因在于，细颗粒流失引发的孔隙比增大无法通过机械压缩完全恢复，导致土体骨架存在永久性缺陷。

尽管9个试样在经历渗流潜蚀和等向固结作用后，细颗粒含量和孔隙比各不相同，但在三轴剪切阶段，各试样的孔隙比均呈现一致的显著增大趋势，说明固结作用导致砂砾土密实度有所提高，在剪切阶段发生了体积膨胀。

4.2　渗透系数的变化特征

渗流作用下，随着细颗粒的流失与局部堵塞的形成，土体孔隙结构不断被重塑，导致其渗透性能持续非线性增大。各试样在不同水力梯度下渗透系数随渗流时间的变化规律如图9所示。以FC10组试样为例，3.2水力梯度下，渗透系数在初始阶段（0~6 s）快速上升，随后逐渐趋于稳定；而在0.8和1.6的水力梯度下，渗透系数的增长较为缓慢。这一变化趋势与图5(b)中累计侵蚀率的演变密切相关：当累计侵蚀率上升至一定值后趋于平稳，说明随着细颗粒的运移与流失，土体内部孔隙逐渐扩大，并在颗粒侵蚀基本稳定后，渗流通道结构亦趋于稳定，从而使渗透系数变化趋缓。总体来看，在相同细颗粒含量条件下，渗透系数随水力梯度的增大而显著增大；而在相同水力梯度条件下，渗透系数则随细颗粒含量的增加而降低。例如，FC10组试样在3.2水力梯度下的渗透系数（0.93 cm/s）最大，是其在0.8水力梯度下的渗透系数（0.10 cm/s）的9倍多，也是FC30试样在3.2水力梯度下（0.47 cm/s）的近2倍。为了揭示其物理机制，进一步分析了渗透系数与孔隙比的关联，发现：渗透系数的终值与孔隙比呈正相关性。这一规律符合经典的Kozeny-Carman渗流理论：潜蚀导致的细颗粒流失直接增大了土体的孔隙率，降低了流体的沿程阻力，从而在宏观上表现为渗透系数的显著增大。

4.3　偏应力-轴向应变关系

潜蚀后土样内部结构发生显著变化，其应力-应变关系的演化依赖于潜蚀后的孔隙比状态。以100 kPa围压为例，试样在潜蚀前后的偏应力-轴向应变关系如图10所示。各试样均呈现典型的应变软化特征：在加载初期，偏应力q随轴向应变{\epsilon }_{\mathrm{a}}的增加呈近似线性增长；当达到峰值强度后，随着剪切位移的累计，土体结构发生剪胀破坏，偏应力逐渐回落，最终趋于稳定。

对于初始细颗粒含量相同的试样，随着细颗粒流失率增大，试样的峰值偏应力显著衰减，应变软化趋势逐渐向应变硬化或弱软化过渡。这表明细颗粒流失不仅降低了土体的抗剪切承载力，还促使其破坏特性由脆性向延性转化。对于恒定水力梯度的试样，峰值偏应力随初始细颗粒含量的增加而增大。此外，峰值偏应力与孔隙比呈现显著的负相关关系：孔隙比越大，意味着土体骨架越松散，初始弹性模量衰减越剧烈，峰值强度越低。上述宏观力学响应的差异本质上归因于土体颗粒体系的结构状态及颗粒咬合程度：在间断级配土（FC10和FC20）中，细颗粒的剥蚀直接导致土体孔隙比增大，削弱了粗颗粒间的约束效应。剪切过程中，高孔隙比状态下的粗颗粒骨架具有更高的自由度，无需克服强烈的剪胀势能与咬合阻力，因此宏观上表现为峰值强度的衰减与软化特征的弱化。反之，当粗颗粒骨架间填充的细颗粒含量较高时（即潜蚀后孔隙比相对较小），颗粒可活动空间受限，且形成更多有效接触（图11），从而构建出更为致密的内部结构，这不仅增强了颗粒间的摩擦与咬合作用，还有助于抑制剪切带的早期发展，从而维持骨架的整体稳定性。

值得注意的是，尽管潜蚀前后试样的细颗粒含量和孔隙比均发生较大变化，但其在大剪切变形阶段的残余强度却基本相同，内在原因在于：各试样中粗颗粒的数量恒定，而砂砾土的残余强度主要由粗颗粒之间的咬合作用所控制。

4.4　体积应变-轴向应变关系

体积应变可反映砂砾土剪切过程中的颗粒重排特征，以100 kPa围压为例，潜蚀前后各试样体积应变随轴向应变的演化规律见图12。结果表明，所有试样的体积应变-轴向应变曲线均表现出先略有压缩、后持续膨胀的特征；随着轴向应变持续增大，体积应变逐渐趋于稳定，试样达到临界破坏状态。对比发现，细颗粒含量越少，潜蚀后孔隙比越大，试样膨胀响应越弱，最终体积应变值越小。其内在机制可从两个方面分析：一方面，在初始结构中，部分细颗粒参与力链传递。渗流作用引发可移动细颗粒流失，导致细颗粒参与形成的接触点数量减少。随着剪切变形的持续发展，颗粒发生重排，低细颗粒含量试样的承载体系发生了转变，即由粗-细颗粒共同承载体系逐渐转向粗颗粒主导的力链结构。因此，此时试样的剪胀变形主要源于粗颗粒骨架间的相互错动。较高细颗粒含量的试样，潜蚀后仍有较多细颗粒保留于土体内部（孔隙比相对较低），试样在剪切时需要更大的体积膨胀来实现颗粒重排，因此表现出更显著的膨胀变形。另一方面，潜蚀后孔隙比的增大为剪切过程中的颗粒滚动提供了更充裕的空间，从而抑制整体体积膨胀的发展。这一机制与松散砂土以剪缩为主的内在原因一致。

4.5　力学性能参数

通过绘制不同围压下的莫尔包络线，可得到潜蚀前后试样内摩擦角的变化规律（见表3）。分析表明，潜蚀作用显著削弱了土体的峰值内摩擦角，且该宏观强度指标的劣化与潜蚀引起的孔隙比增大存在显著的定量关联。为从理论上量化孔隙结构变化对宏观强度的控制机制，本研究引入基于孔隙率变化的强度演化模型（公式8）进行验证：

\phi =\frac{\mathrm{1}-n}{\mathrm{1}-n_{\mathrm{0}}}{\phi }_{\mathrm{0}}（8）

其中：{\phi }_{\mathrm{0}}与n_{\mathrm{0}}为初始试样的内摩擦角和孔隙率；\phi与n分别为潜蚀后试样的内摩擦角和孔隙率。

模拟得到的峰值内摩擦角与潜蚀后孔隙比之间呈现明显的负相关关系（图13）：随着孔隙比的增大，峰值内摩擦角单调下降，线性拟合结果表明孔隙比是控制土体抗剪强度劣化的主导状态变量。同时，模拟结果与理论预测值在演化趋势和数值量级上均高度吻合。这一结果深刻揭示了潜蚀导致的跨尺度力学劣化机理：水力驱动下的细颗粒流失直接导致土体孔隙比升高，孔隙结构的扩容效应削弱了粗颗粒骨架间的咬合与互锁作用，减少了颗粒间的有效接触点数，从而促使宏观内摩擦角随孔隙比增大而显著降低。该规律与梁越等^[19]针对间断级配散粒土的研究结论相符。

值得注意的是，各试样的残余内摩擦角基本恒定，且显著低于对应的峰值内摩擦角。这是由于在残余剪切阶段，土体已进入临界状态，大变形使细颗粒对力链的贡献显著减弱，强度主要由粗颗粒间的摩擦特性控制，因而受初始孔隙比和细颗粒含量的影响较小。Xu等^[42]与Chang等^[43]的试验结果同样表明，潜蚀虽明显改变土体的孔隙状态和峰值强度，但对大变形下的残余强度影响甚微。

5　 结 论

本文聚焦库岸堆积体中的间断级配砂砾土，采用CFD-DEM耦合方法，从宏观和微观角度探究了不同初始细颗粒含量和水力梯度对间断级配砂砾土潜蚀演化进程的影响机制，并引入潜蚀后孔隙比作为核心状态变量，量化揭示了微观结构损伤对宏观力学性能劣化的控制机制。主要结论如下：

1）间断级配砂砾土的潜蚀演化受初始细粒含量与水力梯度的协同控制。细颗粒对粗颗粒骨架孔隙的填充效应显著增强了土体的抗侵蚀能力，导致累计侵蚀率与初始细颗粒含量呈显著的非线性负相关；反之，水力梯度增大显著加剧了细颗粒运移，使累计侵蚀率与水力梯度呈非线性正相关。

2）潜蚀对土体微观结构造成了不可逆损伤。细颗粒流失导致土体强力链断裂与重构、孔隙比增大，诱发土体微观结构向亚稳态演化。同时，潜蚀后孔隙比的增大直接导致实际渗流截面扩大，使渗透系数呈非线性增长，该规律与Kozeny-Carman渗流理论相符。

3）潜蚀诱发的微观结构损伤直接导致宏观力学性能的劣化，其中孔隙比被证实是控制土体抗剪强度劣化与变形特性的主导状态变量。随着潜蚀后孔隙比的增大，土体剪胀效应逐渐弱化，峰值内摩擦角单调线性下降，且演化规律与基于孔隙率的理论预测值高度吻合。其物理本质在于：潜蚀导致的孔隙扩容为剪切过程中的颗粒重排提供了充足的调整空间，削弱了粗颗粒骨架间的互锁咬合作用，导致宏观强度衰减，并显著减弱土体的剪胀行为。

4）潜蚀前后各试样的残余内摩擦角未发生显著变化，表明潜蚀主要削弱了由密实颗粒骨架提供的结构性强度，而在临界状态下，土体的残余强度主要依赖于粗颗粒间的摩擦作用，因此受细颗粒流失的影响相对较小。

参考文献

[1] KENNEY T C， LAU D. Internal stability of granular filters［J］. Canadian Geotechnical Journal， 1985， 22（2）：215-225. DOI： 10.1139/t85-029 .

[2] FOSTER M， FELL R， SPANNAGLE M. A method for assessing the relative likelihood of failure of embankment dams by piping［J］. Canadian Geotechnical Journal， 2000， 37（5）：1025-1061. DOI： 10.1139/t00-029 .

[3] WANG F W， OKEKE A C， KOGURE T， et al. Assessing the internal structure of landslide dams subject to possible piping erosion by means of microtremor chain array and self-potential surveys［J］. Engineering Geology， 2018， 234：11-26. DOI： 10.1016/j.enggeo.2017.12.023 .

[4] LEI X Q， YANG Z J， HE S M， et al. Numerical investigation of rainfall-induced fines migration and its influences on slope stability［J］. Acta Geotechnica， 2017， 12（6）：1431-1446. DOI： 10.1007/s11440-017-0600-y .

[5] 邓华锋， 肖瑶， 方景成， 等. 干湿循环作用下岸坡消落带土体抗剪强度劣化规律及其对岸坡稳定性影响研究［J］. 岩土力学， 2017， 38（9）： 2629-2638. DOI：10.16285/j.rsm.2017.09.021 .

[6] DENG H F， XIAO Y， FANG J C， et al. Shear strength degradation and slope stability of soils at hydro-fluctuation belt of river bank slope during drying-wetting cycle［J］. Rock and Soil Mechanics， 2017， 38（9）： 2629-2638. DOI：10.16285/j.rsm.2017.09.021（Ch ）.

[7] 陈勇， 夏雨， 闵泽鑫， 等. 往返渗流与干湿循环诱发土壤结构损伤的协同效应分析［J］. 农业工程学报， 2025， 41（5）： 57-67. DOI：10.11975/j.issn.1002-6819.202408066 .

[8] CHEN Y， XIA Y， MIN Z X， et al. Analysis on soil structural damage induced by synergistic effects of upward-downward seepage and drying-wetting cycle［J］. Transactions of the Chinese Society of Agricultural Engineering， 2025， 41（5）： 57-67. DOI：10.11975/j.issn.1002-6819.202408066（Ch ）.

[9] MA C Y， PENG M， ZHANG L M， et al. Erosion， deposition and breach evolution of landslide dams composed of various dam material types based on flume tests［J］. Engineering Geology， 2024， 337：107598. DOI： 10.1016/j.enggeo.2024.107598 .

[10] LI S J， RUSSELL A R， MUIR WOOD D. Internal erosion of a gap-graded soil and influences on the critical state［J］. Acta Geotechnica， 2024， 19（8）：5363-5381. DOI： 10.1007/s11440-024-02249-4 .

[11] LIU L L， CHEN R， ZHOU C， et al. Microscopic mechanism of the combined effects of confining pressure and fines content on suffusion in gap-graded underfilled soils［J］. Journal of Hydrology， 2023， 626：130370. DOI： 10.1016/j.jhydrol.2023.130370 .

[12] FONSECA J， SIM W W， SHIRE T， et al. Microstructural analysis of sands with varying degrees of internal stability［J］. Géotechnique， 2015， 65（7）： 620-623. DOI：10.1680/geot.14.D.006 .

[13] LEE H J， KIM I H， CHUNG C K. Evaluation of the internal stability of well-graded silty sand through the long-term seepage test［J］. International Journal of Geo-Engineering， 2021， 12（1）：21. DOI： 10.1186/s40703-021-00151-6 .

[14] CHANG D S， ZHANG L M. Critical hydraulic gradients of internal erosion under complex stress states［J］. Journal of Geotechnical and Geoenvironmental Engineering， 2013， 139（9）：1454-1467. DOI： 10.1061/（ASCE）GT.1943-5606.0000871 .

[15] KE L， TAKAHASHI A. Strength reduction of cohesionless soil due to internal erosion induced by one-dimensional upward seepage flow［J］. Soils and Foundations， 2012， 52（4）：698-711. DOI： 10.1016/j.sandf.2012.07.010 .

[16] PRASOMSRI J， SHIRE T， TAKAHASHI A. Effect of fines content on onset of internal instability and suffusion of sand mixtures［J］. Géotechnique Letters， 2021， 11（3）：209-214. DOI： 10.1680/jgele.20.00089 .

[17] NAQI A L， WATANABE K. Effects of internal erosion on mechanical properties of different soil compositions and its impact on their cyclic response［J］. Soils and Foundations， 2025， 65（1）： 101540. DOI：10.1016/j.sandf.2024.101540 .

[18] MOFFAT R， FANNIN R J. A hydromechanical relation governing internal stability of cohesionless soil［J］. Canadian Geotechnical Journal， 2011， 48（3）：413-424. DOI： 10.1139/T10-070 .

[19] PACHIDEH V， MOHAMMAD HOSSEINI S M MIR. A new physical model for studying flow direction and other influencing parameters on the internal erosion of soils［J］. Geotechnical Testing Journal， 2019， 42（6）：1431-1456. DOI： 10.1520/GTJ20170301 .

[20] ROCHIM A， MAROT D， SIBILLE L， et al. Effects of hydraulic loading history on suffusion susceptibility of cohesionless soils［J］. Journal of Geotechnical and Geoenvironmental Engineering， 2017， 143（7）：04017025. DOI： 10.1061/（ASCE）GT.1943-5606.0001673 .

[21] 梁越， 龚胜勇， 杨咏梅， 等. 间断级配散粒土侵蚀过程及强度演变机理研究［J］. 岩土工程学报， 2024， 46（3）： 632-639. DOI：10.11779/CJGE20221400 .

[22] LIANG Y， GONG S Y， YANG Y M， et al. Erosion process and strength evolution mechanism of gap-graded cohesionless soil［J］. Chinese Journal of Geotechnical Engineering， 2024， 46（3）： 632-639. DOI：10.11779/CJGE20221400（Ch ）.

[23] CHEN L， HE J J， YANG S T， et al. Experimental study on the evolution of the drained mechanical properties of soil subjected to internal erosion［J］. Natural Hazards， 2020， 103（1）： 1565-1589. DOI：10.1007/s11069-020-04049-x .

[24] 吴梦喜， 叶发明， 张琦. 细颗粒流失对砂砾石土本构关系的影响研究［J］. 岩土力学， 2017， 38（6）：1550-1556. DOI： 10.16285/j.rsm.2017.06.002 .

[25] WU M X， YE F M， ZHANG Q. Effect of fine grain loss on the stress-strain relationship of sand and gravel soils［J］. Rock and Soil Mechanics， 2017， 38（6）：1550-1556. DOI： 10.16285/j.rsm.2017.06.002（Ch ）.

[26] OUYANG M， TAKAHASHI A. Influence of initial fines content on fabric of soils subjected to internal erosion［J］. Canadian Geotechnical Journal， 2016， 53（2）： 299-313. DOI：10.1139/cgj-2014-0344 .

[27] RANJBARI P， EBRAHIMI M， EIN-MOZAFFARI F， et al. A critical review of the coupled CFD–DEM method for the simulation of two-phase liquid-solid systems［J］. Powder Technology， 2025， 454：120677. DOI： 10.1016/j.powtec.2025.120677 .

[28] GUO Y， YU X B. Comparison of the implementation of three common types of coupled CFD-DEM model for simulating soil surface erosion［J］. International Journal of Multiphase Flow， 2017， 91：89-100. DOI： 10.1016/j.ijmultiphaseflow.2017.01.006 .

[29] EL-EMAM M A， ZHOU L， SHI W D， et al. Theories and applications of CFD-DEM coupling approach for granular flow： A review［J］. Archives of Computational Methods in Engineering， 2021， 28（7）： 4979-5020. DOI：10.1007/s11831-021-09568-9 .

[30] ZHAO Y， HU Z， ZHENG Y W， et al. Influences of particle size ratio and fines content on the suffusion characteristics of gap-graded soils［J］. Granular Matter， 2024， 26（2）： 52. DOI：10.1007/s10035-024-01425-3 .

[31] XIONG H， ZHANG Z M， BAO X H， et al. Micro-mechanical analysis of particle shape effect on suffusion of gap-graded soils［J］. Computers and Geotechnics， 2024， 165：105925. DOI： 10.1016/j.compgeo.2023.105925 .

[32] 张沛云， 木林隆， 黄茂松. 基于计算流体力学-离散元耦合的相对密实度对断级配土体渗蚀效应的影响研究［J］. 岩土力学， 2024， 45（1）：267-283. DOI： 10.16285/j.rsm.2023.0053 .

[33] ZHANG P Y， MU L L， HUANG M S. Influence of soil relative density on suffusion of gap-graded soil based on coupled computational fluid dynamics-discrete element method［J］. Rock and Soil Mechanics， 2024， 45（1）：267-283. DOI： 10.16285/j.rsm.2023.0053（Ch ）.

[34] BAHMANI M， BRIAUD J L. CFD-DEM simulation of some riprap slopes［C］//Geo-Congress 2023. Los Angeles： American Society of Civil Engineers， 2023： 457-468. DOI：10.1061/9780784484654.046 .

[35] QIU W， LI Y L， ZHANG Y， et al. Numerical investigation on the evolution process of cascade dam-break flood in the downstream earth-rock dam reservoir area based on coupled CFD-DEM［J］. Journal of Hydrology， 2024， 635： 131162. DOI：10.1016/j.jhydrol.2024.131162 .

[36] SUN J Q， LI X A， LI J， et al. Numerical investigation of characteristics and mechanism of tunnel erosion of loess with coupled CFD and DEM method［J］. Catena， 2023， 222： 106729. DOI：10.1016/j.catena.2022.106729 .

[37] HAGER A， KLOSS C， GONIVA C. Combining open source and easy access in the field of DEM and coupled CFD-DEM： LIGGGHTS^®， CFDEM^®coupling and CFDEM^®workbench［J］. Computer Aided Chemical Engineering， 2018， 43：1699-1704. DOI： 10.1016/B978-0-444-64235-6.50296-5 .

[38] AI J， CHEN J F， ROTTER J M， et al. Assessment of rolling resistance models in discrete element simulations［J］. Powder Technology， 2011， 206（3）：269-282. DOI： 10.1016/j.powtec.2010.09.030 .

[39] ZHAO J D， SHAN T. Coupled CFD-DEM simulation of fluid-particle interaction in geomechanics［J］. Powder Technology， 2013， 239：248-258. DOI： 10.1016/j.powtec.2013.02.003 .

[40] DI FELICE R. The voidage function for fluid-particle interaction systems［J］. International Journal of Multiphase Flow， 1994， 20（1）：153-159. DOI： 10.1016/0301-9322（94）90011-6 .

[41] WANG G， DENG Z Z， YANG J， et al. A large-scale high-pressure erosion apparatus for studying internal erosion in gravelly soils under horizontal seepage flow［J］. Geotechnical Testing Journal， 2022， 45（6）： 1037-1053. DOI：10.1520/gtj20220045 .

[42] TAO H， TAO J L. Quantitative analysis of piping erosion micro-mechanisms with coupled CFD and DEM method［J］. Acta Geotechnica， 2017， 12（3）： 573-592. DOI：10.1007/s11440-016-0516-y .

[43] HU Z， ZHANG Y D， YANG Z X. Suffusion-induced deformation and microstructural change of granular soils： A coupled CFD–DEM study［J］. Acta Geotechnica， 2019， 14（3）： 795-814. DOI：10.1007/s11440-019-00789-8 .

[44] KEZDI A. Increase of protective capacity of flood control dikes： Report 1［R］. Budapest： Technical University of Budapest， 1969.

[45] SHIRE T， O’SULLIVAN C， HANLEY K J， et al. Fabric and effective stress distribution in internally unstable soils［J］. Journal of Geotechnical and Geoenvironmental Engineering， 2014， 140（12）：04014072. DOI： 10.1061/（ASCE）GT.1943-5606.0001184 .

[46] LI L M， LI B K， LIU Z Q. Modeling of spout-fluidized beds and investigation of drag closures using OpenFOAM［J］. Powder Technology， 2017， 305： 364-376. DOI：10.1016/j.powtec.2016.10.005 .

[47] XU Z D， ZHANG L M， ZARCH M K， et al. Experimental study of internal erosion in granular soil subject to cyclic hydraulic gradient reversal［J］. Journal of Geotechnical and Geoenvironmental Engineering， 2022， 148（5）： 04022014. DOI：10.1061/（ASCE）GT.1943-5606.0002785 .

[48] CHANG D S， ZHANG L M. A stress-controlled erosion apparatus for studying internal erosion in soils［J］. Geotechnical Testing Journal， 2011， 34（6）： 579-589. DOI： 10.1520/GTJ103889 .

基金资助

国家自然科学基金面上项目(42272317)

国家自然科学基金区域发展联合基金重点项目(U21A2031)