Micro-mechanism of contact erosion between soil layers under cyclic loading based on CFD-DEM coupling

doi:10.11918/202505067

循环荷载作用下土体接触侵蚀的CFD-DEM耦合模拟

doi: 10.11918/202505067

蔡国庆^1,2 ，刁显锋^1,2 ，杨芮^1,2 ，李舰^1,2

1. 城市地下工程教育部重点实验室（北京交通大学），北京 100044

2. 北京交通大学土木建筑工程学院，北京 100044

基金项目: 国家自然科学基金（52425805，U2569208,52278324，42307264）

详细信息

作者简介

蔡国庆（1983―），男，教授，博士生导师

通讯作者

刁显锋，22121053@bjtu.edu.cn

中图分类号: TU443

文献标识码: A

文章编号: 0367-6234(2025)12-0229-16

Micro-mechanism of contact erosion between soil layers under cyclic loading based on CFD-DEM coupling

CAI Guoqing^1,2 ， DIAO Xianfeng^1,2 ， YANG Rui^1,2 ， LI Jian^1,2

1. Key Laboratory of Urban Underground Engineering of Ministry of Education (Beijing Jiaotong University), Beijing 100044 , China

2. School of Civil Engineering, Beijing Jiaotong University, Beijing 100044 , China

摘要

为研究循环荷载作用下土体层间中常见的渗流侵蚀特性，揭示水-力耦合作用下的接触侵蚀的细观机理，考虑循环荷载幅值的影响，建立了一种基于计算流体动力学(CFD)与离散单元法(DEM)耦合的土体接触侵蚀三维计算模型。首先，分析了不同循环荷载作用周次内的颗粒运动规律与空间分布特征。其次，探讨了由于颗粒侵蚀引起的宏观变形响应特性。同时，选取两个局部变形区域，研究了两种典型的颗粒迁移模式。最后，结合力链分析，讨论了侵蚀过程中的颗粒接触力学演化机制。研究结果表明：在一次循环荷载作用的过程中，加载导致的粗颗粒层挤压与卸载产生的应力松弛均是导致细颗粒发生迁移的主要因素；循环荷载作用会使得土层接触面处产生剧烈的颗粒迁移运动，进而导致试样产生较大的轴向变形，同时渗流场会产生向上的水力梯度以促进细颗粒发生泵送迁移。存在一个幅值阈值，当循环荷载幅值超过该阈值后，土体迅速密实，引起侵蚀变形减弱。

关键词

土体接触侵蚀 / CFD-DEM耦合 / 颗粒迁移 / 宏观变形 / 细观机理

Abstract

To investigate the seepage erosion characteristics commonly found between soil layers under cyclic loading and to reveal the micromechanical mechanisms of contact erosion under water-flow coupling, a three-dimensional computational model for soil contact erosion was developed, based on the coupling of computational fluid dynamics (CFD) and discrete element method (DEM), considering the influence of cyclic loading amplitude. Firstly, the movement patterns and spatial distribution characteristics of the particles were analyzed within different cyclic loading periods. Secondly, the macroscopic deformation characteristics resulting from particle erosion were explored. At the same time, two localized deformation regions were selected to study two typical particle migration modes. Finally, the evolution mechanism of particle contact mechanics during the erosion process was discussed, combined with force chain analysis. The results show that during a single cycle of loading, the compression of coarse particles caused by loading and the stress relaxation induced by unloading are the primary factors responsible for the migration of fine particles. Cyclic loading induces intense particle migration at the soil layer interface, resulting in significant axial deformation of the sample. Simultaneously, the seepage field generates an upward hydraulic gradient that promotes the pump-driven migration of fine particles. A threshold value for the loading amplitude exists, and the soil rapidly compacts when the cyclic loading amplitude exceeds this threshold, leading to a reduction in erosion-induced deformation.

Keywords

contact erosion / CFD-DEM coupling / particle migration / macroscopic deformation / mesoscopic mechanism

1 CFD-DEM耦合原理 1.1 DEM控制方程 1.2 CFD控制方程 1.3 颗粒-流体相互作用力方程 1.4 CFD-DEM耦合流程 2 数值模型建立 2.1 模型数值模型 2.2 模型验证 2.3 模拟工况设置 3 颗粒迁移规律 4 宏观变形 4.1 滞回曲线与累积轴向变形 4.2 颗粒侵蚀水平 4.3 孔隙率演化 4.4 渗流场演化 5 细观变形 5.1 局部堵塞演化 5.2 颗粒泵送 6 接触力链演化 7 讨论 8 结论

交通基础设施由于长期暴露在自然环境中，受自然降雨、地下水位上升，以及人类活动等复杂水-力耦合因素影响，土层间极易发生以细颗粒迁移流失为主的接触侵蚀现象，进而引发构筑物变形甚至失稳破坏^[1-3]。因此，透彻认识水-力耦合作用下接触侵蚀的颗粒迁移过程，全面探究宏-细观力学变形特征，对保障土工构筑物防渗设计与排水结构优化具有重要意义。

目前，围绕接触侵蚀机理的既有研究着重采用室内试验方法。早期的关于接触侵蚀的研究主要聚焦于以水力荷载为主的单一影响因素作用。Jin等^[4]采用大型高压侵蚀仪对冲积层的4种代表性级配进行了一系列水平渗流试验，以此评价某土石坝坝基悬挂式防渗墙下砂砾石冲积层的内部稳定性和内部侵蚀行为。Wang等^[5]应用颗粒流失量测试装置进行了土-石填料的渗流试验，监测了不同降雨条件下填料的渗透性、细颗粒流失量和沉降量的变化过程，分析了土壤结构和沉降特性的演化过程。之后，随着对接触侵蚀研究的逐渐深入，研究人员开始研究水-力耦合作用下的接触侵蚀问题。例如，Nguyen等^[6]采用大型循环三轴试验，研究发现相较于土体细颗粒含量，含水量对翻浆冒泥发生过程中道砟层沉降变形的影响更大。Han等^[7]利用循环加载-润湿物理模型系统对铁路路基翻浆冒泥病害进行了研究，得出孔隙水压力是导致路基土层产生颗粒迁移现象的关键因素。

尽管室内试验成果对土体接触侵蚀病害的防治具有一定的指导意义，但对颗粒尺度上的运动规律等微观信息难以获得。因此，从细观尺度上研究土颗粒渗流迁移及其力学特性问题的数值模拟方法成为一种有效的手段^[8-9]。目前，研究人员常采用基于渗流理论的有限元方法（finite element method，FEM）对土体侵蚀问题展开研究^[10-12]。在采用FEM方法建模时，研究对象的尺寸须远大于其每种土体颗粒的尺寸，这使得FEM在模拟宏观尺度下连续介质的变形具有无法比拟的优势^[13]。但是这一假设使得FEM模拟在细观尺度下的离散颗粒的迁移过程与运动规律方面存在不足^[14]。与有限元方法相比，离散单元法（discrete element method，DEM）不仅从颗粒尺度认识土体发生非连续变形与流动的细观机制具有显著优势，还可以从细观视角对孔隙结构和接触力链演化进行可视化^[15]。基于DEM，已开发了诸多流-固耦合计算模型^[16-18]。其中，计算流体力学-离散单元法（CFD-DEM）耦合方法因计算效率高，成本低等优势，已广泛应用于土体渗流侵蚀的数值模拟研究中^[19]。譬如，Zou等^[20]通过CFD-DEM方法研究了土体级配以及水力梯度对渗流侵蚀的影响，分析渗流过程中土体内部孔隙比及颗粒的分布随时间的变化规律。Yin等^[21]为了在微观尺度上更好地理解精细运移过程，基于CFD-DEM模拟了多孔介质中的流体流动并计算了土壤颗粒迁移。Yang等^[22]采用CFD-DEM耦合方法开展了一系列土-结构界面渗流数值试验，分析了渗流装置中不同土层的渗流机理和特性，讨论了界面临界水力坡降的计算公式。

值得注意的是，以往的 CFD-DEM耦合方法研究大多聚焦于粗、细颗粒混合的间隙级配土的侵蚀，对于粗-细颗粒界面分明的土工结构（如土石坝防渗滤层、路基基床排水结构层）在侵蚀过程中接触面处发生的接触侵蚀关注甚少^[23-24]，尤其针对循环荷载作用下土层间接触面的颗粒迁移现象与力学特征的研究更是鲜见报道。另外，采用CFD-DEM方法模拟分析间隙级配土体颗粒侵蚀迁移问题的研究中，流场设置大部分聚焦于固定水力梯度^[25-26]，而对于循环荷载变化引起的初始静水状态流场变化的特征尚不清楚。因此，借助于CFD-DEM方法全面研究土层间接触侵蚀过程中的颗粒迁移现象、粗细颗粒接触面力学特性、流场变化规律以及宏观变形响应的工作亟须开展。

综上所述，本文采用CFD-DEM耦合方法，考虑不同循环荷载幅值作用，对土层间接触侵蚀过程进行模拟，研究颗粒迁移运动规律，分析循环荷载幅值对宏观变形响应特征、细观颗粒迁移过程以及粗细颗粒接触特性的影响，模拟结果可为防渗滤层设计、路基翻浆冒泥病害的防治提供参考。

1 CFD-DEM耦合原理

1.1 DEM控制方程

颗粒运动采用离散单元法模拟计算，其运动方程包括平动方程与转动方程，分别服从牛顿第二定律和欧拉第二定律。颗粒运动的控制方程分别为：

m_{i} \frac{d u_{p, i}}{d t} = F_{i}^{g} + F_{i}^{f} + \sum_{j = 1}^{n_{i}^{c}} F_{i j}^{c}

(1)

I_{i} \frac{d ω_{i}}{d t} = \sum_{j = 1}^{n_{i}^{c}} (M_{t i j} + M_{r i j})

(2)

式中：m_i、I_i分别为颗粒i的质量与转动惯量，u_p，_i、ω_i分别为颗粒i的速度与角速度，

F_{i}^{g} 、 F_{i j}^{c} ， F_{i}^{f}

分别为作用在颗粒上的重力、颗粒间作用力以及流体作用力，

M_{t i j} 、 M_{r i j}

分别为颗粒间产生相对滑动时产生的切向力矩以及滚动摩擦产生的力矩。

颗粒间接触力采用Hertz-Mindlin接触模型计算，该模型基于Hertz理论^[27]与Mindlin理论^[28]，模型参数的物理含义明确，可以反映颗粒接触的非线性变化特征^[29-30]。此外，颗粒由于形状不规则而产生的力学效应是不容忽视的，但是直接构建不规则形状颗粒往往会增加计算负担^[31]。近年来已有研究人员通过在颗粒间增加滚动力矩来模拟不规则颗粒形状效应^[32-35]，因此本文在Hertz-Mindlin接触模型的基础上引入 EPSD（elastic-plastic spring dashpot）模型^[35]来考虑颗粒因形状不规则产生的力学效应。在EPSD模型中，由颗粒形成产生的滚动阻力矩可分为弹簧力矩和黏性阻力力矩，可表示为

M_{r} = M_{r}^{k} + M_{r}^{d}

(3)

式中：

M_{r}^{k}

为弹簧转矩，

M_{r}^{d}

为黏性阻尼转矩。

M_{r}^{k}

在通过一个方向上存在一个最大值

M_{r}^{max}

，在未到达

M_{r}^{max}

之前，呈现出线性变化。当到达或超过最大值时，其值不再发生变化。具体计算式如下：

(4)

式中：k_r为滚动刚度，Δθ_r为接触的两颗粒之间的角速度增量，μ_r为滚动摩擦系数，R^*为相接触两颗粒的等效半径，r_i、r_j分别为相互接触粒子i、j的半径，F_n为颗粒法向接触力。

黏性阻尼转矩

M_{r}^{d}

与相接触两颗粒滚动的瞬时相对角速度

θ_{r e l}^{d}

及阻尼系数

C_{r}^{d}

有关。当

M_{r}^{k}

被充分调动时，

M_{r}^{d}

为0，即在

M_{r}^{d}

到达

M_{r}^{m a x}

之前

M_{r}^{d}

才会起作用，原因在于

M_{r}^{d}

的作用为防止颗粒产生滚动振荡以保证运动稳定的同时又耗散了能量。t+Δt时刻的阻尼转矩

M_{r ， t + Δ t}^{d}

可表示为

(5)

1.2 CFD控制方程

流体相由Navier-Stokes（N-S）方程来描述，控制方程满足质量守恒定律与动量守恒定律，具体为：

\frac{\partial ε ρ_{w}}{\partial t} + \nabla (ε ρ_{w} u_{w}) = 0

(6)

\begin{matrix} \frac{\partial (ε ρ_{w} u_{w})}{\partial t} + \nabla \cdot (ε ρ_{w} u_{w} u_{w}) = \\ - ε \nabla p + \nabla \cdot ε τ + ε ρ_{w} g - F_{p} \end{matrix}

(7)

式中：ρ_w为流体密度，u_w为流体速度，ε为孔隙率，p为压力张量，τ为流体黏性应力张量，

τ = μ (\nabla u_{w} + {(\nabla u_{w})}^{T})

；μ为流体的动力黏度系数，g为重力加速度，F_p为颗粒与流体的相互作用。

1.3 颗粒-流体相互作用力方程

在未解析的CFD-DEM耦合方法中，颗粒-流体间的相互作用力可表示为

F_{p} = F_{d} + F_{g} + F_{v} + F_{f}

(8)

式中：压力梯度力是由流体中的压强差产生的，

F_{g}

- V_{p i} \nabla p

，其中V_p_i为颗粒i的体积；黏滞力是由流体的黏性产生的，

F_{v}

V_{p i} \nabla \cdot τ

；浮力是影响土粒在流体中浮沉的重要因素之一，

F_{f}

\frac{1}{6} π ρ_{w} d_{p}^{3} g

，其中d_p为颗粒的直径。

拖曳力是由流体运动和粒子运动之间的速度差引起的作用力，是粒子在流体中运动时所受到的阻力。本文采用Di Felice 拖曳力计算模型^[36]，该模型不仅具有连续函数形式的拖曳力计算公式，同时还兼具计算精度高的优势^[37-39]。其表达式为

F_{d} = f (ε) \frac{C_{d} π ρ_{w} d_{p}^{2} |u_{w} - u_{p}| (u_{w} - u_{p})}{8}

(9)

式中：u_w、u_p分别为液体和颗粒的运动速度，C_d为拖曳力系数，其与雷诺数Re_p相关，具体为：

(10)

{R e}_{p} = \frac{2 ε ρ_{w} r_{p} |u_{w} - u_{p}|}{μ_{w}}

(11)

式（9）中，f（ε）为与孔隙率相关的函数，对于多颗粒体系，其可表示为

f (ε) = ε^{- χ}

(12)

上式表征了颗粒之间相互影响的孔隙率修正系数，指数χ可计算为

χ = 3.7 - 0.65 e x p (- \frac{{(1.5 - {l o g}_{10} R e_{p})}^{2}}{2})

(13)

1.4 CFD-DEM耦合流程

计算颗粒的运动方程采用颗粒离散元软件LIGGGHTS，求解流场N-S方程使用计算流体力学软件OpenFOAM，最后采用开源软件CFDEM实现DEM与CFD的并行计算，该程序可以实现多核CPU并行计算，以提高求解效率。CFD-DEM耦合的计算流程见图1。

图1CFD-DEM耦合流程

Fig.1CFD-DEM coupling process

2 数值模型建立

2.1 模型数值模型

离散元试样包括下层细颗粒与上层粗颗粒，且上、下层颗粒尺寸分布差异较大，因此在循环荷载、水力渗流条件等作用下细颗粒极易发生向上迁移。鉴于计算机的计算量限制，模型所研究的间隙级配试样为实际分布土样的简化^[40]，根据TB 10001―2016《铁路路基设计规范》，路基填料C2组良好级配细砂，考虑1~3 mm粒径。粗颗粒土选择基床表层Ⅱ型级配，颗粒级配曲线见图2（a）。制样尺寸为90 mm×90 mm×180 mm，其中下层细颗粒土高度为60 mm，制样过程采用分层欠压法，使土样均匀^[41]，最终模型见图2（b）。制样完成后根据实际工况饱和固结，并采用伺服加压的方式施加30 kPa的轴向压力以模拟上部设施质量^[42]，四周墙体保持固定，使试样处于均匀应力状态。

根据颗粒粒径，本文粗、细颗粒均属于石英砂颗粒^[43]。因此颗粒的物理参数根据已有的研究进行选取：1）天然石英砂的密度采用常用值2 650 kg/m³^[44-46]。2）由于泊松比对土体宏观变形的影响较小^[47]，因此本文采用常用值0.3^[44-45，48]。3）Chand等^[49]的研究发现，在保证颗粒的重叠率不超过2%的前提下，弹性模量可采用常用值7 GPa（7.0×109 Pa）^{[45，48，50]}。此外，颗粒的接触参数通过试验研究的标定结果，具体为：1）Zhou 等^[51]的试验结果建议在DEM建模过程中石英砂的恢复系数应取0.90；2）根据Mitchell等^[52]的试验结果，石英砂的摩擦系数取0.50；3）Cheng等^[47]通过试验标定结果，建议石英砂的滚动摩擦系数取0.26，滚动刚度系数取0.19。试样的物理参数见表1。

图2计算模型

Fig.2Calculation model

表1颗粒材料的物理参数与接触参数^[23-24]

Tab.1Physical parameters and calculation parameters of granular materials^[23-24]

为保证整个侵蚀过程流体域始终覆盖固体域，同时为保证流体流动的均匀性，在荷载施加方向上的流体域尺寸大于固体域。CFD区域的顶面与底面分别设置为出口与入口，四周边界设置为不透水壁面。其中，进、出口边界设置为均匀压力条件，均为0。壁面处的压力梯度设置为0。速度边界条件是将侧面边界设置为无滑移边界条件，在进、出口边界处设置速度梯度为0。为防止在出口处发生流体的回流而影响计算稳定性，将回流速度设置为0。最终整个流体域为静水状态。流体相采用密度为1 000 kg/m³，动力黏滞系数μ=0.001 Pa·s的流体水。

在DEM中，由于采用了显式积分的方案，时间步长在小于式（14）计算的临界Rayleigh条件时能够保证进行稳定求解计算^[45]。本文采用的DEM时间步长为1×10^-6 s。CFD的时间步长应满足式（15）定义的Courant数小于1的条件^[1]，确保压力-速度耦合流体相方程得到收敛解。此外，Tao等^[53]建议，CFD的时间步长是DEM时间步长的50~100倍时，既可以满足耦合计算的求解精度，又可以保证具有较高的计算效率。因此，CFD的时间步长设置为1×10^-4 s。

Δ t_{Rayleigh} = \frac{π r \sqrt{\frac{ρ}{G}}}{0.1631 ν + 0.8766}

(14)

C_{Courant} = \frac{u_{w} Δ t_{C F D}}{Δ l}

(15)

式中：ρ为颗粒的密度，G为颗粒的剪切模量，ν为颗粒的泊松比，Δl为流体网格的长度。

2.2 模型验证

为了验证本文建立的CFD-DEM耦合模型的准确性，选择Zhang等^[54]的试验研究作为验证案例，其研究的重点为循环荷载作用下接触侵蚀过程中的颗粒迁移，与本文具有一定的相似性。采用双向循环加载，循环荷载示意图见图3（a）。循环荷载的幅值为3 kPa，循环荷载频率分别为5、15、25 Hz，平均应力为12 kPa。共进行3组模拟，每一组模拟情况的循环荷载加载20 s。

图3（b）为浊度增长率与频率的关系。采用Zhang 等^[54]的试验结果进行对比分析验证。观察发现，本文所建立的CFD-DEM耦合模型的计算结果与Zhang等^[54]的试验结果相差不大。同时本模型的模拟结果与试验拟合曲线的相关系数R2=0.829 1，拟合程度比较好，一定程度上说明了CFD-DEM耦合模型的准确性。但是需要说明的是，随着频率的增大，数值结果与试验结果及其拟合曲线的差距之间增大，尤其是当频率为25 Hz时。因此，基于上述验证结果，本文选择5 Hz作为后续研究的荷载加载频率。

图3CFD-DEM耦合计算结果与Zhang等^[54]的试验结果对比

Fig.3Comparison of the CFD-DEM coupling simulation results and the experimental results of Zhang et al. ^[54]

2.3 模拟工况设置

为研究循环荷载作用下的接触侵蚀的细观发生机理，采用如图4所示的更符合实际交通循环荷载形式的正弦循环荷载加载方案^[55-58]，循环荷载输入公式为

F = A + p s i n (10 π t + 1.5 π)

(16)

式中：F为循环荷载，A为静止荷载与循环荷载幅值之和，p为施加的荷载幅值，t为模拟时间。循环荷载频率为5 Hz，动应力幅值分别为40、60、80、100 kPa。

图4循环荷载加载方案

Fig.4Schematic of the cyclic load of the CFD-DEM model

3 颗粒迁移规律

通过CFD-DEM方法对循环荷载作用下土体接触侵蚀变形过程进行探究。本文以f=5 Hz，p=40 kPa工况为例进行分析，循环加载次数为0、25、50、75、100次后的试样整体变形图见图5。如图5（a）所示，上部设施静载作用下，粗颗粒层与细颗粒层的变形发生在两者相接触处，粗颗粒层挤压下部细颗粒层产生少量相对位移，导致细颗粒层表面发生不均匀变形，部分表层细颗粒进入底层粗颗粒孔隙通道，但未发生粗颗粒土层的整体堵塞。施加25次循环荷载后，如图5（b）所示，粗、细颗粒层发生明显的相对运动，粗颗粒层的底层部分颗粒被细颗粒覆盖，细颗粒层表层颗粒侵入粗颗粒层的孔隙通道。如图5（c）~图5（e）所示，随着循环荷载的持续作用，粗、细颗粒层之间的相对运动加剧，粗颗粒层向下挤压细颗粒层，使得其底部颗粒完全侵入细颗粒层，细颗粒向上发生明显的迁移现象，试样整体变形明显。

图5试样总体变形

Fig.5Overall deformation of the sample

为分析在循环荷载作用下颗粒迁移的运动规律，绘制了不同循环次数周期内细颗粒速度图，见图6。在第1次循环周期内，开始加载后，如图6（a）、（b）所示，原本稳定的土体由于粗颗粒的向下挤压，使得细颗粒层边缘处部分较为松散的颗粒开始产生运动，发生迁移。且随着加载的持续，发生迁移的颗粒数量增多。当开始卸载时，如图6（c）所示，观察发现细颗粒层中产生速度的颗粒数量急速增大，发生迁移的细颗粒主要集中分布于粗颗粒层底部周围及孔隙处，并在孔隙下方沿细颗粒层向下延伸形成颗粒集中迁移带。随着卸载的持续，如图6（d）所示，细颗粒的运动更加剧烈，发生迁移的颗粒数量进一步增多，原本多个颗粒集中迁移带连通形成更大颗粒集中迁移区域。

如图6（e）、（f）所示，在第25次循环荷载的加载过程中，细颗粒层只有少量颗粒具有较小速度，没有较为明显的颗粒迁移现象发生。如图6（g）、（h）所示，卸载后，在粗颗粒侧面的孔隙通道处颗粒迁移速度较大，细颗粒集中迁移区域缩小。分析原因为：随着循环荷载的施加，粗颗粒向下挤压使得细颗粒层逐渐密实，抑制了颗粒迁移的发生。

如图6（i）、（j）所示，第50次加载时，在细颗粒层的表层有部分及与粗颗粒作用处有少量颗粒发生迁移。如图6（k）、（l）所示，在卸载过程中细颗粒集中迁移区域较第25次循环时更大，粗、细颗粒接触处细颗粒运动更为剧烈，并再次形成较为明显的细颗粒迁移区域。在此次循环荷载作用下，加-卸载过程中细颗粒迁移再次加剧原因为：在循环荷载作用下，使得原本静态的流场发生扰动，在渗流路径上产生水压差，进而形成较大的水力梯度驱动颗粒发生迁移。在第75次循环荷载作用时，如图6（m）~图6（p）所示，已有大量细颗粒沿着粗颗粒层孔隙通道发生迁移并逐渐填充孔隙，粗颗粒层底部颗粒因向下挤压而完全进入细颗粒层中。加载过程中迁移的表层细颗粒数量增多，卸载过程中粗、细颗粒接触面处的相互作用程度进一步加剧。第100次循环荷载作用下，如图6（q）、（r）所示，加载过程中，细颗粒层表面产生运动的颗粒数量进一步增多。如图6（s）、（t）所示，在卸载过程中，深层颗粒发生迁移的颗粒数量减少，细颗粒的迁移运动有所减弱，分析原因为：粗颗粒上部孔隙通道内填充的细颗粒数量较多，形成了颗粒的堆积堵塞，抑制了颗粒迁移。

图6颗粒迁移过程（f=5 Hz，p=40 kPa，截面x=20 mm）

Fig.6Particle migration process (f=5 Hz， p=40 kPa) with fine particle velocity (x=20 mm section)

4 宏观变形

4.1 滞回曲线与累积轴向变形

图7为不同循环荷载幅值作用下的试样在不同加载阶段的滞回曲线以及累积轴向变形图。随着循环次数的增大，滞回曲线的逐渐闭合围成滞回圈，滞回圈的面积越来越小，表明随着循环次数的发展，模型逐渐被压实，加载所消耗的能量逐渐降低。循环荷载不断作用下，滞回圈随着轴向应变的累积而逐渐向右移动，同时沿逆时针方向发展，斜率逐渐增大，尤其在1~25次循环荷载周期内，斜率的变化最为明显。随着循环荷载作用次数的增大，斜率的变化不再显著，说明在加载过程中，试样在模拟初期较大的轴向变形，试样的刚度逐渐增大，发生应变硬化，并在模拟后期逐渐趋于稳定。如图7（a）~（e）所示，随着循环荷载幅值的增大，滞回圈的面积增大。试样在不同荷载幅值作用下的累积轴向变形如图7（f）所示，试样的累积轴向变形随荷载幅值的增大而增大。需要注意的是，p=100 kPa时，试样的轴向变形小于p=80 kPa时的轴向变形，原因在于荷载幅值较大时，荷载作用传递给土体的能量较大，使得土体快速密实发生硬化，从而抑制了土体变形。

图7滞回曲线演化与累积轴向变形

Fig.7Hysteresis curve evolution and cumulative axial strain under different cyclic loading amplitude

4.2 颗粒侵蚀水平

不同循环荷载幅值作用下，侵蚀前、后细颗粒各层的颗粒质量对比见图8。相较于初始状态，各层细颗粒在模拟结束时均发生了明显的迁移现象。其中细颗粒表层的质量差异较大，中层次之，说明接触侵蚀在粗、细颗粒交界面处极易发生。

图8侵蚀前、后质量对比

Fig.8Mass comparison before and after erosion

参照Bendahmane等^[59]的研究，定义每秒每平方米面积上的细粒流失量为细粒侵蚀率。如图9（a）所示为不同荷载幅值作用下细颗粒侵蚀率随时间的变化。细颗粒侵蚀率总体呈现出先增大后减小的变化趋势，说明该现象与过往的实验现象^[59]接近。细颗粒侵蚀率在t=1 s时达到最大，说明在循环荷载加载的初始阶段，细颗粒极易发生流失，该阶段粗颗粒的向下挤压为导致细颗粒迁移发生的主要因素。随后细颗粒侵蚀率下降并保持在较低水平，此时主要发生因水力梯度作用引发的颗粒迁移。细颗粒侵蚀率随荷载幅值的增大呈现出递增的趋势。为描述粗、细颗粒相互作用下细颗粒侵蚀情况，将细颗粒侵蚀占比定义为最终侵入粗颗粒层的细颗粒质量占渗流开始前细颗粒总质量的比值。各工况下颗粒迁移比如图9（b）所示，粗颗粒层中的细颗粒质量占比随荷载幅值增大，当荷载幅值过大时，饱和砂土被压密实，使得细颗粒的迁移现象被抑制，且荷载幅值作用下的变化规律相对明显，说明荷载幅值对颗粒迁移的影响较大。

图9细颗粒侵蚀率与侵蚀占比

Fig.9Fine particle erosion rate and erosion ratio

为评估粗颗粒层的污染程度，采用Ebrahimi等^[60]定义考虑污染材料种类的粗颗粒污染指数，其表达式为

F_{1} = F_{C} (G_{S O} / G_{S})

(17)

式中：G_S0为参考相对质量密度，取2.6；G_S为污染材料真实相对密度，F_C为Selig等^[61]定义的粗颗粒污染物含量，其表达式为

F_{C} = w_{4} + w_{200}

(18)

式中：w₄为污染粗颗粒中直径小于4.750 mm的颗粒质量分数，w₂₀₀为污染粗颗粒中粒径小于0.075 mm的颗粒质量分数。数值试验所有颗粒均大于0.075 mm且小于4.750 mm，故取F_C=1。不同工况下粗颗粒层污染指数计算结果见表2。模拟结果已发生严重的污染，尤其在荷载幅值较大时，对粗颗粒层的污染更为严重。

表2不同工况下粗颗粒层污染指数

Tab.2Pollution index of coarse stratum granulosum under different working conditions

4.3 孔隙率演化

现有试验条件无法实时测量孔隙率变化，但是CFD-DEM耦合方法则可解决上述问题。图10为不同工况下，模型细粒层孔隙率变化时程曲线。从图10中可以看出，孔隙率均呈现出先减小后增大趋势。其原因在于初始阶段，粗颗粒嵌入作用剧烈，细颗粒层受到挤压，细颗粒层受挤压而逐渐密实，因此孔隙率减小。随后由于下部细颗粒层逐渐密实并形成稳定性较强的土骨架，粗颗粒的嵌入作用减弱，此时由于流场开始产生水力梯度，颗粒迁移过程开始发生以泵送迁移为主的细颗粒流失，导致孔隙率开始逐渐增大。另外，在高应力幅值作用下，p=100 kPa时，土体孔隙率下降幅度最大，印证了土体密实进而抑制土体轴向变形及颗粒迁移。

图10细颗粒层孔隙率演化

Fig.10Variation curve of porosity in fine-grained soil layers

4.4 渗流场演化

定义细颗粒层水力梯度，其计算公式为

i = \frac{Δ p}{ρ g Δ h}

(19)

式中：Δp为细颗粒层顶面与底面处孔隙水压力差，ρ为流体密度，Δh为所选截面高差。图11为细颗粒各层水力梯度随时间变化曲线。曲线整体呈现出增大的变化趋势。流场初始为静水状态，由于粗颗粒嵌入发生较大的轴向变形，此时粗颗粒层向下挤压，土体孔隙变小导致水压力开始累积。循环荷载的持续施加引起水压力的增大，从而产生向上的水压力梯度，水流在压力梯度的作用下产生向上的流速。随着时间增大，水力梯度增大。4种荷载幅值作用下对应的最大水力梯度分别为3.39、3.78、4.50、4.05。即当p＜100 kPa时，呈现出随荷载幅值增大而整体增大的变化趋势，促进颗粒泵送漂浮发生。根据已有研究^[62-64]，该水力梯度满足诱发土体颗粒发生迁移的条件。在p=100 kPa时，压力梯度减小的原因在于：在高应力水平下，土体更加密实从而抑制了孔隙变形，进一步使得水压力的累积减弱，导致细颗粒层上、下截面压差减小。

图11细颗粒层水力梯度演化

Fig.11Hydraulic gradient evolution in the fine particle layer

5 细观变形

5.1 局部堵塞演化

堆积堵塞在粗、细颗粒接触面的不同位置随机发生，研究试样内部局部区域的堵塞行为有助于进一步揭示接触侵蚀的机理。如图12所示为选定的试样中局部堵塞区域、粗细颗粒接触面堆积堵塞演化过程及其接触力链的演化过程。如图12（b）所示，在第5次循环荷载作用时，粗颗粒的孔隙通道内仅有少量细颗粒，颗粒之间仍存在较大孔隙。通过分析如图12（c）的接触力链的演化过程发现，此时孔隙通道内的因颗粒间接触形成的接触力链网络稀疏，在加载过程中，部分接触力链极易发生断裂破坏，形成颗粒迁移缺口。虽然在卸载过程中该处力链会重塑，但是由于此时力链的强度较低，在下一个循环荷载作用过程中仍会发生断裂破坏，为颗粒的迁移、堆积、堵塞提供条件。

如图12（d）所示，在第100次循环荷载作用时，细颗粒已填充粗颗粒孔隙通道，并在一次循环荷载作用下并未发现有明显的颗粒迁移现象，说明此时该孔隙通道内已发生了颗粒的堆积堵塞。进一步分析接触力链发现，如图12（e）所示，此时孔隙通道内由细颗粒组成的接触力链已形成密布的网络，接触力链之间缝隙较小，且在加载过程中，细颗粒参与粗颗粒之间力的传递，使得接触力强度增大，形成强度与密度更高的接触力网络。虽然在卸载过程中有少量强接触力链的削弱，但是在孔隙通道内因颗粒的堆积已形成稳定的传力结构，进而阻碍了颗粒的迁移。

为对比分析不同工况下粗颗粒层侵蚀过程中的堵塞情况，绘制了粗颗粒层的流速变化曲线，见图13。平均流速随时间呈现出先增大后减小的变化规律。初始阶段，流速不断增大，原因在于水压力的累积效应，导致产生水力梯度，流体开始流动。随着水力梯度的增大，在粗颗粒孔隙通道内逐渐形成细颗粒优势渗流通道，此时孔隙通道未发生堵塞。随后粗颗粒孔隙内的细颗粒数量增多，孔隙内发生细颗粒的堆积，并逐渐发生堵塞，导致粗颗粒层表观流速降低。在实际工程中，防渗滤层渗透性损伤、铁路路基排水不畅等问题均是由孔隙通道堵塞引起。

图12局部堵塞演化

Fig.12Evolution of the local packing profile

图13粗颗粒层平均流速变化

Fig.13Average flow velocity in coarse-particle layers

5.2 颗粒泵送

如图14（a）所示，循环荷载作用下，粗、细颗粒的接触面处极易发生细颗粒的迁移流失。图14（b）所示为两次循环荷载周期内的泵送过程，在两次循环荷载作用下，该孔隙通道内均发生了明显的颗粒迁移现象。颗粒迁移过程的接触力链演化如图14（c）所示，在泵送过程中，颗粒在孔隙通道内迁移过程中产生较弱的接触力链甚至不产生力链，说明细颗粒在孔隙通道内接触较少或者不发生接触，即颗粒在迁移过程中为悬浮状态，此时颗粒迁移主要受水力梯度影响。颗粒随加-卸载过程而呈现出“抽-吸”的“泵送”迁移规律，而该孔隙通道为泵送通道。XOY平面内的接触力链演化过程如图14（d）所示，观察发现：在泵送通道的底部，接触力链形成了一个稳定的“缺口”，该“缺口”空隙较大且允许细颗粒通过，为颗粒泵送提供了有利条件。该“缺口”的稳定性较高，在第1次循环荷载作用下，该通道未发生明显的变形。但是在第2次循环荷载加-卸载过程中，该“缺口”发生变形破坏，并在卸载的过程中发生了重塑，形成了更大的 “缺口”，进一步导致更多颗粒发生泵送流失。

图14泵送过程演化

Fig.14Evolution of the pumping

不同工况下，泵送发生的时间及泵送导致颗粒迁移的累积质量见图15。荷载幅值对试样因泵送而产生的累积效应较为明显，随着荷载幅值的增大，泵送的颗粒质量增大，发生泵送的时间缩短。虽然增大荷载幅值在一定程度上抑制了泵送，但也会导致土体颗粒流失，影响构筑物的安全服役。

图15颗粒泵送质量与发生时间

Fig.15Fine particle pumping under different load amplitudes

6 接触力链演化

图16为不同循环次数下f=5 Hz，p=40 kPa工况为例的接触力网络图。在施加循环荷载之前，如图16（a）所示，粗颗粒层接触力链强度较大，说明上部静载大部分粗颗粒层承担，并在粗、细颗粒交界面处分散传递给细颗粒层，同时细颗粒层的接触力网络分布较为均匀，此时试样处于相对稳定的状态。如图16（b）所示，在施加25次循环荷载之后，粗颗粒层接触力链的强度增大，细颗粒层接触力网络发生了明显的重新分布，其中在粗颗粒下方产生强接触力链集中分布，说明在此处产生了应力集中，发生了较大的塑性变形。同时在粗、细颗粒交界面的孔隙处接触力链的强度削弱，说明在该处更易形成颗粒堆积，且该处的细颗粒更易发生侵蚀迁移。随着循环次数的增多，如图16（c）、（d）所示，原本粗颗粒下方的强接触力链以及细颗粒层更深层的接触力链被削弱。当循环荷载作用第100次时，对比图16（d）、（e），接触力链的分布变化不大，说明此时接触力链的分布逐渐趋于稳定，土体仍能形成相对稳定的受力结构。

图16试样接触力链演化

Fig.16Contact force chain network after different cyclic loads

土体作为散体材料，将力的传递方式通常划分为强力接触力网和弱力接触力网。参考已有研究^[65]，定义特征力F^*=1.2〈F_n〉作为强、弱接触力的界限值，其中〈F_n〉为细颗粒层所有法向力的平均值，大于F^*则认为为强接触力链，反之则为弱接触力链。不同工况下，细颗粒层接触力链数量图及其占比曲线图见图17。

图17试样强接触力占比变化

Fig.17Change in the proportion of strong contact force

从图17可以发现，模拟开始阶段，粗颗粒的侵入压实作用，使得强接触力链占比增大，土体骨架发生变形，宏观表现为土体轴向位移较大。之后粗颗粒的不断嵌入，粗颗粒下部细颗粒层逐渐密实，产生新的强接触力，但是嵌入作用使得粗颗粒周围的细颗粒发生向上迁移，导致部分强力链发生断裂。随着水力泵送逐渐剧烈，部分表层强接触力链发生断裂，宏观表现为泵送现象。随着荷载幅值的增大，初始状态强接触力占比减小，说明荷载幅值的增大导致土体发生更大的塑性变形。

7 讨论

本文针对循环荷载作用下土体接触侵蚀这一关键岩土工程问题，采用CFD-DEM耦合数值模拟方法，系统揭示了循环荷载条件下土体侵蚀的宏、细观演化机制。在工程实践中，此类侵蚀现象广泛存在于路基、河岸边坡以及地下结构等重大基础设施中，其诱发因素主要包括地震、火车与汽车行驶等复杂应力环境。研究成果不仅为路基不均匀沉降、翻浆冒泥等典型侵蚀灾害的防控提供了理论依据，同时也为工程安全预警体系的完善奠定了科学基础。

需要指出的是，本文在数值建模过程中存在若干需要改进的方面：首先，为确保计算收敛性对循环荷载进行了简化处理；其次，尽管采用EPSD模型表征颗粒运动层面的力学行为，但颗粒几何形状对侵蚀和流体力学响应的影响未被完全还原；最后，受计算资源限制，采用的网格分辨率尚不足以精确捕捉孔隙尺度下的渗流场演化特征。基于上述局限性，建议后续研究可着重从以下方向展开：①融合计算几何学方法或数字扫描技术，建立真实颗粒形态的三维重构模型；②综合考虑水力梯度、不同孔隙率以及围压耦合等多物理场作用；③开发自适应网格算法与并行计算技术，以提升大规模颗粒-流体耦合系统的计算精度与效率。

8 结论

针对土层间的土接触侵蚀问题，运用CFD-DEM耦合模型以及计算方法，分析了不同循环荷载幅值作用下颗粒运动规律与宏观变形响应特征，探讨了颗粒细观迁移行为与接触力网络变化规律，得到以下结论：

1）导致细颗粒迁移的原因有两种，即粗颗粒的嵌入作用以及水力泵送作用，初期颗粒迁移以嵌入作用为主。之后随水力梯度的增大，主要发生水力泵送为主的颗粒迁移；颗粒迁移过程中粗、细颗粒接触面相互作用剧烈，细颗粒沿粗颗粒层孔隙通道发生颗粒迁移，造成表层细颗粒流失严重。

2）循环荷载作用导致原本稳定的土体接触力网络发生破坏，宏观表现为试样发生较大轴向变形；粗颗粒的向下挤压使得细颗粒层稳定土骨架形成的同时又抑制了嵌入作用；另外，水力作用使得细颗粒表面及粗颗粒孔隙处的接触力链强度降低，甚至发生断裂，宏观表现为颗粒漂浮现象的产生。

3）循环荷载作用下，土体孔隙结构的变化引发初始静水流场发生扰动而产生水压力，形成向上的水力梯度是导致流体运动并携带细颗粒发生泵送的主要原因。细颗粒在粗颗粒通过优势渗流孔隙通道迁移运动，随着细颗粒的不断增多并发生堆积，形成粗颗粒层孔隙通道的局部堵塞。

4）增大荷载幅值对接触侵蚀过程中微观颗粒迁移与接触力学特性、流场变化以及模型宏观变形响应的影响程度显著。但是土体存在一个循环荷载幅值的阈值，当荷载幅值超过该阈值后，土体因高应力作用而快速密实，进而抑制了土体变形流失与渗流破坏。

图1CFD-DEM耦合流程

Fig.1CFD-DEM coupling process

下载: 全尺寸图片

图2计算模型

Fig.2Calculation model

下载: 全尺寸图片

图3CFD-DEM耦合计算结果与Zhang等^[54]的试验结果对比

Fig.3Comparison of the CFD-DEM coupling simulation results and the experimental results of Zhang et al. ^[54]

下载: 全尺寸图片

图4循环荷载加载方案

Fig.4Schematic of the cyclic load of the CFD-DEM model

下载: 全尺寸图片

图5试样总体变形

Fig.5Overall deformation of the sample

下载: 全尺寸图片

图6颗粒迁移过程（f=5 Hz，p=40 kPa，截面x=20 mm）

Fig.6Particle migration process (f=5 Hz， p=40 kPa) with fine particle velocity (x=20 mm section)

下载: 全尺寸图片

图7滞回曲线演化与累积轴向变形

Fig.7Hysteresis curve evolution and cumulative axial strain under different cyclic loading amplitude

下载: 全尺寸图片

图8侵蚀前、后质量对比

Fig.8Mass comparison before and after erosion

下载: 全尺寸图片

图9细颗粒侵蚀率与侵蚀占比

Fig.9Fine particle erosion rate and erosion ratio

下载: 全尺寸图片

图10细颗粒层孔隙率演化

Fig.10Variation curve of porosity in fine-grained soil layers

下载: 全尺寸图片

图11细颗粒层水力梯度演化

Fig.11Hydraulic gradient evolution in the fine particle layer

下载: 全尺寸图片

图12局部堵塞演化

Fig.12Evolution of the local packing profile

下载: 全尺寸图片

图13粗颗粒层平均流速变化

Fig.13Average flow velocity in coarse-particle layers

下载: 全尺寸图片

图14泵送过程演化

Fig.14Evolution of the pumping

下载: 全尺寸图片

图15颗粒泵送质量与发生时间

Fig.15Fine particle pumping under different load amplitudes

下载: 全尺寸图片

图16试样接触力链演化

Fig.16Contact force chain network after different cyclic loads

下载: 全尺寸图片

图17试样强接触力占比变化

Fig.17Change in the proportion of strong contact force

下载: 全尺寸图片

表1颗粒材料的物理参数与接触参数^[23-24]

Tab.1Physical parameters and calculation parameters of granular materials^[23-24]

下载: 全尺寸图片

表2不同工况下粗颗粒层污染指数

Tab.2Pollution index of coarse stratum granulosum under different working conditions

下载: 全尺寸图片

图1CFD-DEM耦合流程

Fig.1CFD-DEM coupling process

图2计算模型

Fig.2Calculation model

图3CFD-DEM耦合计算结果与Zhang等^[54]的试验结果对比

Fig.3Comparison of the CFD-DEM coupling simulation results and the experimental results of Zhang et al. ^[54]

图4循环荷载加载方案

Fig.4Schematic of the cyclic load of the CFD-DEM model

图5试样总体变形

Fig.5Overall deformation of the sample

图6颗粒迁移过程（f=5 Hz，p=40 kPa，截面x=20 mm）

Fig.6Particle migration process (f=5 Hz， p=40 kPa) with fine particle velocity (x=20 mm section)

图7滞回曲线演化与累积轴向变形

Fig.7Hysteresis curve evolution and cumulative axial strain under different cyclic loading amplitude

图8侵蚀前、后质量对比

Fig.8Mass comparison before and after erosion

图9细颗粒侵蚀率与侵蚀占比

Fig.9Fine particle erosion rate and erosion ratio

图10细颗粒层孔隙率演化

Fig.10Variation curve of porosity in fine-grained soil layers

图11细颗粒层水力梯度演化

Fig.11Hydraulic gradient evolution in the fine particle layer

图12局部堵塞演化

Fig.12Evolution of the local packing profile

图13粗颗粒层平均流速变化

Fig.13Average flow velocity in coarse-particle layers

图14泵送过程演化

Fig.14Evolution of the pumping

图15颗粒泵送质量与发生时间

Fig.15Fine particle pumping under different load amplitudes

图16试样接触力链演化

Fig.16Contact force chain network after different cyclic loads

图17试样强接触力占比变化

Fig.17Change in the proportion of strong contact force

表1颗粒材料的物理参数与接触参数^[23-24]

Tab.1Physical parameters and calculation parameters of granular materials^[23-24]

表2不同工况下粗颗粒层污染指数

Tab.2Pollution index of coarse stratum granulosum under different working conditions

图1CFD-DEM耦合流程

Fig.1CFD-DEM coupling process

图2计算模型

Fig.2Calculation model

图3CFD-DEM耦合计算结果与Zhang等^[54]的试验结果对比

Fig.3Comparison of the CFD-DEM coupling simulation results and the experimental results of Zhang et al. ^[54]

图4循环荷载加载方案

Fig.4Schematic of the cyclic load of the CFD-DEM model

图5试样总体变形

Fig.5Overall deformation of the sample

图6颗粒迁移过程（f=5 Hz，p=40 kPa，截面x=20 mm）

Fig.6Particle migration process (f=5 Hz， p=40 kPa) with fine particle velocity (x=20 mm section)

图7滞回曲线演化与累积轴向变形

Fig.7Hysteresis curve evolution and cumulative axial strain under different cyclic loading amplitude

图8侵蚀前、后质量对比

Fig.8Mass comparison before and after erosion

图9细颗粒侵蚀率与侵蚀占比

Fig.9Fine particle erosion rate and erosion ratio

图10细颗粒层孔隙率演化

Fig.10Variation curve of porosity in fine-grained soil layers

图11细颗粒层水力梯度演化

Fig.11Hydraulic gradient evolution in the fine particle layer

图12局部堵塞演化

Fig.12Evolution of the local packing profile

图13粗颗粒层平均流速变化

Fig.13Average flow velocity in coarse-particle layers

图14泵送过程演化

Fig.14Evolution of the pumping

图15颗粒泵送质量与发生时间

Fig.15Fine particle pumping under different load amplitudes

图16试样接触力链演化

Fig.16Contact force chain network after different cyclic loads

图17试样强接触力占比变化

Fig.17Change in the proportion of strong contact force

表1颗粒材料的物理参数与接触参数^[23-24]

Tab.1Physical parameters and calculation parameters of granular materials^[23-24]

表2不同工况下粗颗粒层污染指数

Tab.2Pollution index of coarse stratum granulosum under different working conditions

CAI Guoqing, DIAO Xianfeng, HE Xuzhen,et al. Micromechanical analysis of contact erosion under cyclic loads using the coupled CFD-DEM method[J]. Géotechnique,2025,75(5):649. DOI:10.1680/jgeot.23.00458

FOSTER M, FELL R. Assessing embankment dam filters that do not satisfy design criteria[J]. Journal of Geotechnical and Geoenvironmental Engineering,2001,127(5):398. DOI:10.1061/(asce)1090-0241(2001)127:5(398)

LIU Kaiwen, SU Qian, YUE Fei,et al. Effects of suffosion-induced contact variation on dynamic responses of saturated roadbed considering hydro-mechanical coupling under high-speed train loading[J]. Computers and Geotechnics,2019,113:103095. DOI:10.1016/j.compgeo.2019.103095

JIN Wei, DENG Zezhi, WANG Gang,et al. Internal erosion experiments on sandy gravel alluvium in an embankment dam foundation emphasizing horizontal seepage and high surcharge pressure[J]. Water,2022,14(20):3285. DOI:10.3390/w14203285

WANG Yueyue, MAO Xuesong, WU Qian,et al. Study on the influence of seepage conditions with different rainfall intensities on the structural evolution of soil-rock mixture filler[J]. Geofluids,2022,2022(1):7694663. DOI:10.1155/2022/7694663

NGUYEN T T, INDRARATNA B. Influence of varying water content on permanent deformation of mud-fouled ballast[J]. Transportation Geotechnics,2023,38:100919. DOI:10.1016/j.trgeo.2022.100919

HAN Bowen, CAI Guoqing, SU Yanlin,et al. Investigating the particle migration and dynamic characteristics in ballasted track subgrade under intermittent load-wetting coupling[J]. Transportation Geotechnics,2025,51:101501. DOI:10.1016/j.trgeo.2025.101501

CAI Guoqing, SU Yanlin, ZHOU Annan,et al. An elastic-viscoplastic model for time-dependent behavior of unsaturated soils[J]. Computers and Geotechnics,2023,159:105415. DOI:10.1016/j.compgeo.2023.105415

SU Yanlin, CAI Guoqing, YIN Fengjie,et al. Elastic-viscoplastic model for coarse-grained soil considering particle breakage[J]. Computers and Geotechnics,2024,165:105962. DOI:10.1016/j.compgeo.2023.105962

MEKONNEN A. Train induced pore water pressure generation model:numerical comparison[J]. Journal of Vibroengineering,2019,21(4):952. DOI:10.21595/jve.2018.20253

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

ZHAO H Y, INDRARATNA B, NGO T. Numerical simulation of the effect of moving loads on saturated subgrade soil[J]. Computers and Geotechnics,2021,131:103930. DOI:10.1016/j.compgeo.2020.103930

RAVICHANDRAN N, MACHMER B, KRISHNAPILLAI H,et al. Micro-scale modeling of saturated sandy soil behavior subjected to cyclic loading[J]. Soil Dynamics and Earthquake Engineering,2010,30(11):1212. DOI:10.1016/j.soildyn.2010.05.002

WANG Yin, CHENG Kuang, YANG Yefeng,et al. Microscopic mechanical analysis of sand production using a new arbitrary resolved-unresolved CFD-DEM model[J]. International Journal of Multiphase Flow,2022,149:103979. DOI:10.1016/j.ijmultiphaseflow.2022.103979

WANG Xiukai, HUANG Bo, TANG Yao,et al. Microscopic mechanism and analytical modeling of seepage-induced erosion in bimodal soils[J]. Computers and Geotechnics,2022,141:104527. DOI:10.1016/j.compgeo.2021.104527

LOMINÉ F, SCHOLTS L, SIBILLE L,et al. Modeling of fluid-solid interaction in granular media with coupled lattice Boltzmann/discrete element methods:application to piping erosion[J]. International Journal for Numerical and Analytical Methods in Geomechanics,2013,37(6):577. DOI:10.1002/nag.1109

MA Gang, ZHOU Wei, CHANG Xiaolin,et al. A hybrid approach for modeling of breakable granular materials using combined finite-discrete element method[J]. Granular Matter,2016,18(1):7. DOI:10.1007/s10035-016-0615-3

CHEN Tianhao, HU Zheng, YANG Zhongxuan,et al. A resolved CFD-DEM investigation into the onset of suffusion:effect of confining pressure and stress anisotropy[J]. International Journal for Numerical and Analytical Methods in Geomechanics,2023,47(16):3018. DOI:10.1002/nag.3611

蔡国庆, 刁显锋, 杨芮, 等. 基于CFD-DEM的流-固耦合数值建模方法研究进展[J]. 哈尔滨工业大学学报,2024,56(1):17. CAI Guoqing, DIAO Xianfeng, YANG Rui,et al. Research progress of fluid-solid coupling model based on CFD-DEM coupling[J]. Journal of Harbin Institute of Technology,2024,56(1):17. DOI:10.11918/202309001

ZOU Yuhua, CHEN Qun, CHEN Xiaoqing,et al. Discrete numerical modeling of particle transport in granular filters[J]. Computers and Geotechnics,2013,47:48. DOI:10.1016/j.compgeo.2012.06.002

YIN Yanzhou, CUI Yifei, TANG Yao,et al. Solid-fluid sequentially coupled simulation of internal erosion of soils due to seepage[J]. Granular Matter,2021,23(2):20. DOI:10.1007/s10035-020-01076-0

YANG Xuemin, XU Zengguang, CHAI Junrui,et al. Numerical investigation of the seepage mechanism and characteristics of soil-structure interface by CFD-DEM coupling method[J]. Computers and Geotechnics,2023,159:105430. DOI:10.1016/j.compgeo.2023.105430

CAI Guoqing, DIAO Xianfeng. A coupled VOF-DEM investigation into contact erosion under the effects of cyclic loading and seepage path[J]. Canadian Geotechnical Journal,2025,62:1. DOI:10.1139/cgj-2024-0303

DIAO Xianfeng, CAI Guoqing, YANG Rui,et al. Contact erosion of soil layers with different water table levels under cyclic loading using VOF-DEM coupled method[J]. Computers and Geotechnics,2025,180:107090. DOI:10.1016/j.compgeo.2025.107090

GUO Yuan, WANG Chen, LU Mengmeng,et al. Numerical investigations on evolution mechanism of contact erosion in layered soils considering particle-scale parameters[J]. Transportation Geotechnics,2023,41:101023. DOI:10.1016/j.trgeo.2023.101023

GUO Yuan, XIAO Jinqin, WANG Chen. A coarse graining CFD-DEM simulation on contact erosion of layered cohesionless soils[J]. Canadian Geotechnical Journal,2024,61(3):534. DOI:10.1139/cgj-2022-0408

HERTZ H R. Uber die beruhrung fester elastischer korper und uber die harte[J]. Verhandlung des Vereins zur Beforderung des GewerbefleiBes, Berlin,1882:449. DOI:10.1515/crll.1882.92.156

MINDLIN R D, DERESIEWICZ H. Elastic spheres in contact under varying oblique forces[J]. Journal of Applied Mechanics,1953,20(3):327. DOI:10.1115/1.4010702

CIANTIA M O, ARROYO M, BUTLANSKA J,et al. DEM modelling of cone penetration tests in a double-porosity crushable granular material[J]. Computers and Geotechnics,2016,73:109. DOI:10.1016/j.compgeo.2015.12.001

SALAZAR A, SÁEZ E, PARDO G. Modeling the direct shear test of a coarse sand using the 3D Discrete Element Method with a rolling friction model[J]. Computers and Geotechnics,2015,67:83. DOI:10.1016/j.compgeo.2015.02.017

JIANG M J, YU H S, HARRIS D. A novel discrete model for granular material incorporating rolling resistance[J]. Computers and Geotechnics,2005,32(5):340. DOI:10.1016/j.compgeo.2005.05.001

JIANG Mingjing, SHEN Zhifu, WANG Jianfeng. A novel three-dimensional contact model for granulates incorporating rolling and twisting resistances[J]. Computers and Geotechnics,2015,65:147. DOI:10.1016/j.compgeo.2014.12.011

KONDIC L. Dynamics of spherical particles on a surface:collision-induced sliding and other effects[J]. Physical Review E, Statistical Physics, Plasmas, Fluids,and Related Interdisciplinary Topics,1999,60(1):751. DOI:10.1103/physreve.60.751

GONIVA C, KLOSS C, DEEN N G,et al. Influence of rolling friction on single spout fluidized bed simulation[J]. Particuology,2012,10(5):582. DOI:10.1016/j.partic.2012.05.002

AI Jun, CHEN Jianfei, ROTTER J M,et al. Assessment of rolling resistance models in discrete element simulations[J]. Powder Technology,2011,206(3):269. DOI:10.1016/j.powtec.2010.09.030

DI FELICE R. The voidage function for fluid-particle interaction systems[J]. International Journal of Multiphase Flow,1994,20(1):153. DOI:10.1016/0301-9322(94)90011-6

WANG Tianyu, WANG Shuai, SHEN Yansong. Particle-scale study of gas-solid flows in a bubbling fluidised bed:effect of drag force and collision models[J]. Powder Technology,2021,384:353. DOI:10.1016/j.powtec.2021.02.034

AGRAWAL V, SHINDE Y, SHAH M T,et al. Effect of drag models on CFD-DEM predictions of bubbling fluidized beds with Geldart D particles[J]. Advanced Powder Technology,2018,29(11):2658. DOI:10.1016/j.apt.2018.07.014

MARCHELLI F, HOU Qinfu, BOSIO B,et al. Comparison of different drag models in CFD-DEM simulations of spouted beds[J]. Powder Technology,2020,360:1253. DOI:10.1016/j.powtec.2019.10.058

ZHANG Peiyun, MU Linlong, HUANG Maosong. A coupled CFD-DEM investigation into hydro-mechanical behaviour of gap-graded soil experiencing seepage erosion considering cyclic hydraulic loading[J]. Journal of Hydrology,2023,624:129908. DOI:10.1016/j.jhydrol.2023.129908

JIANG M J, KONRAD J M, LEROUEIL S. An efficient technique for generating homogeneous specimens for DEM studies[J]. Computers and Geotechnics,2003,30(7):579. DOI:10.1016/S0266-352X(03)00064-8

ISRAR J, INDRARATNA B. Mechanical response and pore pressure generation in granular filters subjected to uniaxial cyclic loading[J]. Canadian Geotechnical Journal,2018,55(12):1756. DOI:10.1139/cgj-2017-0393

FANNIN R J, MOFFAT R. Observations on internal stability of cohesionless soils[J]. Géotechnique,2006,56(7):497. DOI:10.1680/geot.2006.56.7.497

CHEN Fan, XIONG Hao, WANG Xiang,et al. Transmission effect of eroded particles in suffusion using the CFD-DEM coupling method[J]. Acta Geotechnica,2023,18(1):335. DOI:10.1007/s11440-022-01568-8

MU Linlong, ZHANG Peiyun, SHI Zhenhao,et al. Coupled CFD-DEM investigation of erosion accompanied by clogging mechanism under different hydraulic gradients[J]. Computers and Geotechnics,2023,153:105058. DOI:10.1016/j.compgeo.2022.105058

WANG Tuo, WANG Pei, YIN Zhenyu,et al. DEM-DFM modeling of suffusion in calcareous sands considering the effect of double-porosity[J]. Computers and Geotechnics,2022,151:104965. DOI:10.1016/j.compgeo.2022.104965

CHENG Kuang, WANG Yin, YANG Qing,et al. Determination of microscopic parameters of quartz sand through tri-axial test using the discrete element method[J]. Computers and Geotechnics,2017,92:22. DOI:10.1016/j.compgeo.2017.07.017

LIU Yajing, WANG Lizhong, YIN Zhenyu,et al. Micromechanical investigation of suffusion and its consequence in gap-graded granular soil at the opening of underground structures using the coupled CFD-DEM method[J]. Computers and Geotechnics,2023,158:105395. DOI:10.1016/j.compgeo.2023.105395

CHAND R, ALI KHASKHELI M, QADIR A,et al. Discrete particle simulation of radial segregation in horizontally rotating drum:effects of drum-length and non-rotating end-plates[J]. Physica A: Statistical Mechanics and Its Applications,2012,391(20):4590. DOI:10.1016/j.physa.2012.05.019

LIU Yajing, WANG Lizhong, HONG Yi,et al. A coupled CFD-DEM investigation of suffusion of gap graded soil:coupling effect of confining pressure and fines content[J]. International Journal for Numerical and Analytical Methods in Geomechanics,2020,44(18):2473. DOI:10.1002/nag.3151

ZHOU Y C, XU B H, YU A B,et al. Numerical investigation of the angle of repose of monosized spheres[J]. Physical Review E, Statistical, Nonlinear,and Soft Matter Physics,2001,64(2 Pt 1):021301. DOI:10.1103/physreve.64.021301

MITCHELL J K, SOGA K, O’SULLIVANC. Fundamentals of soil behavior[M].3rd ed. Hoboken: John Wiley & Sons Inc.,2005. DOI:10.1002/9781119832348

TAO Hui, TAO Junliang. Quantitative analysis of piping erosion micro-mechanisms with coupled CFD and DEM method[J]. Acta Geotechnica,2017,12(3):573. DOI:10.1007/s11440-016-0516-y

ZHANG Sheng, GAO Feng, HE Xuzhen,et al. Experimental study of particle migration under cyclic loading:effects of load frequency and load magnitude[J]. Acta Geotechnica,2021,16(2):367. DOI:10.1007/s11440-020-01137-x

DING Yu, JIA Yu, WANG Xuan,et al. The influence of geotextile on the characteristics of railway subgrade mud pumping under cyclic loading[J]. Transportation Geotechnics,2022,37:100831. DOI:10.1016/j.trgeo.2022.100831

DUONG T V, CUI Yujun, TANG A M,et al. Effects of water and fines contents on the resilient modulus of the interlayer soil of railway substructure[J]. Acta Geotechnica,2016,11(1):51. DOI:10.1007/s11440-014-0341-0

GU Chuan, YE Xingchi, WANG Jun,et al. Resilient behavior of coarse granular materials in three-dimensional stress state[J]. Canadian Geotechnical Journal,2020,57(9):1280. DOI:10.1139/cgj-2019-0353

HAN Bowen, CAI Guoqing, LI Jian,et al. Investigating particles migration caused by mud pumping in ballasted track subgrade under cyclic loading-wetting coupling[J]. Transportation Geotechnics,2022,37:100830. DOI:10.1016/j.trgeo.2022.100830

BENDAHMANE F, MAROT D, ALEXIS A. Experimental parametric study of suffusion and backward erosion[J]. Journal of Geotechnical and Geoenvironmental Engineering,2008,134(1):57. DOI:10.1061/(asce)1090-0241(2008)134:1(57)

EBRAHIMI A, TINJUM J M, EDIL T B. Deformational behavior of fouled railway ballast[J]. Canadian Geotechnical Journal,2015,52(3):344. DOI:10.1139/cgj-2013-0271

SELIG E T, WATERS J M. Track geotechnology and substructure management[M]. London: Thomas Telford Publishing,1994. DOI:10.1680/tgasm.20139

GAO Feng, HE Xuzhen, ZHANG Sheng. Pumping effect of rainfall-induced excess pore pressure on particle migration[J]. Transportation Geotechnics,2021,31:100669. DOI:10.1016/j.trgeo.2021.100669

INDRARATNA B, ISRAR J, RUJIKIATKAMJORN C. Geometrical method for evaluating the internal instability of granular filters based on constriction size distribution[J]. Journal of Geotechnical and Geoenvironmental Engineering,2015,141(10):04015045. DOI:10.1061/(asce)gt.1943-5606.0001343

INDRARATNA B, SINGH M, NGUYEN T T,et al. Internal instability and fluidisation of subgrade soil under cyclic loading[J]. Indian Geotechnical Journal,2022,52(5):1226. DOI:10.1007/s40098-022-00616-0

THORNTON C, ZHANG L. On the evolution of stress and microstructure during general 3D deviatoric straining of granular media[J]. Géotechnique,2010,60(5):333. DOI:10.1680/geot.2010.60.5.333

Publication Statement

Journal Subscription

1 CFD-DEM耦合原理

2 数值模型建立

3 颗粒迁移规律

4 宏观变形

5 细观变形

6 接触力链演化

7 讨论

8 结论