颗粒流沿斜面越过墙体：光滑与非光滑离散元方法的交叉对比

《Computational Particle Mechanics》：Granular flow past a wall down an incline: Cross-comparison between smooth and non-smooth discrete element methods

　　

基于离散元方法（DEM）的模拟研究了干燥颗粒材料沿粗糙斜面冲击刚性墙体的稳态流动。该颗粒流-障碍物装置中同时存在准静态、稠密惯性及稀薄流动状态，被用作严格测试以交叉对比采用YADE的标准DEM与采用SICONOS的非光滑DEM。研究人员应用相同的后处理流程（包括粗粒化及时间平均）以评估不同DEM实现方式的后果。在三个不同位置（给料槽出口处的流动入口、未受干扰区域及障碍物附近）比较了两种方法获得的运动学变量和应力场。在宏观尺度上，两种方法总体吻合良好；特别是μ(I)流变学方面恢复出相似趋势，且刚性墙上的冲击力相当。然而，在接触尺度上出现显著差异，导致两种模拟中的能量耗散不同。在底部及壁面附近的边界处观察到局部差异，此处接触力实现方式的差异最为显著，从而导致基底滑移速度和应力分布的不同。

颗粒流沿斜面冲击刚性墙体的交叉对比研究：光滑与非光滑离散元方法的系统性分析

**研究背景与意义**

颗粒流广泛存在于自然灾害（如山体滑坡、雪崩、泥石流）及工业过程（如农业、制药、土木工程中的物料处理与运输）中。尽管成分看似简单，颗粒流却表现出结合固体、液体和气体特征的复杂行为，挑战经典连续介质描述。过去数十年，对颗粒流的理解从经验描述发展为统一物理框架，通过实验和数值模拟识别出不同流动状态，并构建了如μ(I)流变学（描述颗粒材料从类固态到类液态转变的本构模型）等本构模型。数值模拟中，离散元方法（DEM）被广泛采用，但存在两大流派：标准（光滑）DEM基于显式时间积分和阻尼弹性接触律，而非光滑DEM采用隐式时间积分和基于单边约束的接触模型。尽管两者目标一致，但其接触公式和时序方案的差异可能导致结果偏差。当前，缺乏对两种方法在复杂流动-障碍物交互场景中系统性的交叉对比，该场景同时存在准静态、稠密惯性与稀薄流动状态，可作为严格测试案例。本研究旨在通过对比标准DEM（采用YADE软件）与非光滑DEM（采用SICONOS软件）在模拟颗粒流冲击刚性墙时的表现，评估两者可靠性与局限性，探究宏观行为与颗粒尺度量的异同，为颗粒流数值模拟方法选择提供依据。论文发表在《Computational Particle Mechanics》。

**主要关键技术方法**

研究人员采用离散元方法（DEM）进行二维数值模拟：使用YADE（标准DEM）和SICONOS（非光滑DEM）两种软件。模拟几何为倾斜角Θ=28°的粗糙平面，含下游固定刚性墙状障碍物，n=25,000个颗粒直径均一分布（d±0.15d, d=1mm）。底部粗糙度由固定球形颗粒生成。后处理采用粗粒化（coarse-graining）和时间平均方法，提取连续场量（体积分数、速度、应力等）。时间步长基于能量耗散收敛分析确定：YADE取Δt_Y=1×10^-7s，SICONOS取Δt_S=2.5×10^-5s。模拟参数：颗粒密度ρ_p=2500 kg/m³，摩擦系数μ=0.5，恢复系数e=0.5。在标准DEM中，杨氏模量E=10⁸Pa，泊松比ν=0.3，确保刚体极限满足（重叠δ/d<10^-3）。非光滑DEM中采用θ=1/2的θ-method和100次高斯-赛德尔迭代。

**研究结果**

**3.1 接触尺度比较**

*3.1.1 库仑准则（Coulomb’s criterion）：* 通过分析单个时间步（t=2.1s）所有接触的粘-滑条件发现，YADE中大量接触在力处于库仑锥内时仍存在非零滑移速度，而SICONOS呈现近乎理想的粘-滑剖面。引入度量Δ_ss量化偏差，YADE均值(Δ?_ss)_Y≈10^-2，SICONOS均值(Δ?_ss)_S≈4×10^-5，表明标准DEM未严格满足粘-滑条件，导致后续宏观差异。

*3.1.2 重叠（Overlaps）：* 重叠分布分析显示，YADE中重叠值显著小于SICONOS，且受时间步影响较小。SICONOS中较大重叠源于非光滑DEM接触律直接基于速度与冲量而非重叠，旨在控制接触功率交换，不构成根本问题。

*3.1.3 接触力做功（Contact force works）：* 累计耗散能量D_N和D_T在YADE中更大，尤其是法向分量，表明YADE能量耗散更高，可能与库仑准则未完全执行及法向阻尼系数仅对等径接触准确恢复e值有关。

**3.2 深度平均流动尺度比较**

*3.2.1 体积分数（Volume fraction）：* 两种方法均显示流动中部密度低于高应力区（如储料门附近和障碍物上游），准静态死区（dead zone）底部体积分数最高。差异主要源于重叠值不同：高速区（区B）非光滑DEM重叠较大，导致体积分数略高和自由表面略低；标准DEM中重叠随压力增大，在死区底部增加体积分数。全局差异低于5%。

*3.2.2 速度（Velocity）：* 速度剖面显示特征Bagnold型形状，基底存在有限滑移速度。YADE预测的流动速度整体低于SICONOS，基底差异达25%，归因于库仑准则执行不严格及YADE较高能量耗散，导致速度降低和自由表面升高。

*3.2.3 压力和剪应力（Pressure and shear stress）：* 压力与剪应力场几何模式相似，与体积分数对应。在死区与流动层过渡区（y/d≈8-12）剪应力因强剪切而增强。基底附近薄层内应力下降，可能受粗粒化核半径影响。SICONOS在死区内预测略高压力值，可能与低频内部运动和少量接触相关。

*3.2.4 障碍物受力（Force on the obstacle）：* 法向力剖面上部几乎恒定，下部增大但峰值位置两种方法不同。切向力在上部向上、下部向下，可能反映死区内次级流动结构。峰值位置和大小高度敏感于局部力链结构，由于死区准静态性质，粗粒化和时间平均不足以平滑波动，导致差异。

*3.2.5 μ(I)流变学：* 在μ(I)关系中，两种方法均复现准静态与碰撞态单一曲线，YADE数据点更离散。拟合参数（μ₀, μ_∞, I₀）有差异（SICONOS: 0.35, 0.70, 0.27; YADE: 0.30, 0.66, 0.13），主趋势一致但低I区偏离，可能因仅单一倾角Θ限制I范围。?(I)关系拟合也高度吻合（?₀≈0.78, C≈-0.33）。偏离主趋势的数据点位于喷出流（airborne jet）起始处，表明μ(I)流变学在应力与应变率张量非共线区域（障碍物附近及边界）不适用。

**结论讨论与总结**

本研究对比分析了光滑与非光滑离散元方法（DEM）在模拟颗粒流越过斜面刚性墙时的表现。研究结论表明：尽管数值公式存在根本差异，两种方法在宏观流动行为描述上总体一致，接触模型采用几乎相同e和μ值能给出相似结果，证实该倾斜颗粒流-壁面交互场景中宏观物理主要由几何驱动而非精确接触行为。然而，在接触尺度出现显著差异，导致能量耗散不同，并在边界条件附近（粗糙底部和障碍物壁面）引起运动学变量与应力场的局部偏差，尤其是基底滑移速度的明显差异。死区内接触数量少增强了接触实现方式的敏感性，使得SICONOS模拟中死区应力和障碍物受力略高。两种方法均观察到稳健共性：体积分数、速度、应力场高度相似，且?(I)和μ(I)流变关系吻合良好。但两种模拟也一致指出，在应力张量与应变率张量非共线区域（如死区上方高度非均匀流动区、储料门附近、底部边界附近），μ(I)流变学不适用。这些结果为更广泛的参数研究（不同倾角、材料及多分散性）奠定了基础，有助于量化趋势的稳健性并确定μ(I)流变学的有效范围，最终为稠密颗粒流越过障碍物建立更全面的预测框架。