引言
砂土、矿石、药片、粮食——颗粒材料既像固体又能流动，其力学行为由颗粒间的离散接触主导，连续介质假设在颗粒尺度上失效。Cundall 与 Strack 于 1979 年提出的离散元方法（Discrete Element Method, DEM）干脆放弃连续化：直接追踪每个颗粒的平动与转动，用接触力学模型计算颗粒间作用力，显式推进牛顿第二定律。这一”自下而上”的思路使 DEM 成为岩土工程、矿山机械、化工流化床、制药压片等领域分析颗粒流的核心工具。本文介绍软球接触模型的力学基础、时间步长选取的稳定性约束，并给出一个二维颗粒堆积过程的完整 DEM 实现。
理论基础
软球接触模型
DEM 的主流是软球模型（soft-sphere model）：允许颗粒发生微小重叠，以重叠量 $\delta_n$ 表征接触变形。法向力采用弹簧-阻尼器（spring-dashpot）组合：
$$F_n = k_n \delta_n - \eta_n v_n$$
其中 $k_n$ 为法向接触刚度，$\eta_n$ 为法向阻尼系数，$v_n$ 为两颗粒相对速度的法向分量（接近时为负）。阻尼耗散动能，宏观表现为碰撞恢复系数 $e$；对线性弹簧-阻尼模型，二者有解析关系 $\eta_n = -2\ln e \sqrt{m_{\text{eff}} k_n}\big/\sqrt{\pi^2 + \ln^2 e}$，$m_{\text{eff}} = m_i m_j/(m_i + m_j)$ 为等效质量。切向力模拟摩擦，标准实现是”切向弹簧 + 库仑滑移”：切向力随相对切向位移增量累积，且幅值受库仑条件截断：
$$|\mathbf{F}_t| \le \mu |F_n|$$
$\mu$ 为颗粒间摩擦系数。更精细的法向接触可用 Hertz 理论 $F_n = \frac{4}{3} E^* \sqrt{R^*}, \delta_n^{3/2}$（$E^*$、$R^*$ 为等效模量与等效半径），配合相应的非线性阻尼；线性模型则以参数标定直观、计算廉价取胜。
运动方程与时间步长
每个颗粒满足刚体动力学方程：
$$m_i \frac{d\mathbf{v}_i}{dt} = \sum_c \mathbf{F}^c + m_i \mathbf{g}, \qquad I_i \frac{d\boldsymbol{\omega}_i}{dt} = \sum_c \mathbf{T}^c$$
显式积分的稳定性要求时间步长远小于接触振动的固有周期。估计接触持续时间 $t_c \approx \pi\sqrt{m_{\text{eff}}/k_n}$，工程上通常取 $\Delta t \le t_c/10 \sim t_c/50$；更保守的判据以 Rayleigh 波在颗粒表面的传播时间为上限。这一约束（典型 $\Delta t$ 在 $10^{-6}\sim10^{-5}$ s 量级）正是 DEM 计算昂贵的根源：模拟 1 秒物理时间往往需要 $10^5$ 步以上的接触搜索与力计算。
数值方法与实现
下面给出 25 个玻璃珠在二维容器内下落、碰撞并堆积成床的完整模拟。接触检测采用全对搜索（$O(N^2)$，演示规模足够；大规模问题需换用 Verlet 列表或网格分箱），切向摩擦采用粘性正则化近似——生产级代码会存储每个接触对的切向位移历史以再现真正的弹塑性摩擦：
import numpy as np
# ---------------- 物理参数（玻璃珠） ----------------
N = 25 # 颗粒数
R = 0.01 # 半径 [m]
rho_p = 2500.0 # 密度 [kg/m^3]
kn = 1.0e5 # 法向接触刚度 [N/m]
e_rest = 0.3 # 恢复系数
mu_fric = 0.5 # 摩擦系数
g = 9.81
box_w = 0.20 # 容器宽度 [m]
dt = 1.0e-5
nsteps = 50000 # 模拟 0.5 s
m = rho_p * (4.0 / 3.0) * np.pi * R**3
# 由恢复系数换算法向阻尼系数（线性弹簧-阻尼模型）
ln_e = np.log(e_rest)
eta_n = -2.0 * ln_e * np.sqrt(0.5 * m * kn) / np.sqrt(np.pi**2 + ln_e**2)
# ---------------- 初始条件：容器下部规则排布 + 随机扰动 ----------------
rng = np.random.default_rng(7)
nx = 5
rows = []
for i in range(nx): # 相邻层错开半格，促使颗粒滚入间隙
xs = np.linspace(0.04, 0.15, nx) + (i % 2) * 0.014
rows.append(np.column_stack([xs, np.full(nx, 0.04 + i * 2.4 * R)]))
pos = np.concatenate(rows)[:N]
pos += rng.uniform(-0.3 * R, 0.3 * R, pos.shape) # 打破完美对称
vel = np.zeros((N, 2))
def contact_forces(pos, vel):
"""软球模型：法向弹簧-阻尼 + 简化库仑摩擦（粘性正则化），全向量化"""
d = pos[:, None, :] - pos[None, :, :] # r_i - r_j
dist = np.linalg.norm(d, axis=2)
safe = np.where(dist > 1e-12, dist, 1.0)
nij = d / safe[:, :, None]
overlap = np.clip(2 * R - dist, 0.0, None)
mask = (overlap > 0) & (dist > 1e-12) # 仅在真实接触时作用
vij = vel[:, None, :] - vel[None, :, :]
vn = np.einsum('ijk,ijk->ij', vij, nij)
vt = vij - vn[:, :, None] * nij
vtn = np.linalg.norm(vt, axis=2)
fn = (kn * overlap - eta_n * vn) * mask
tij = vt / (vtn[:, :, None] + 1e-30)
ft = np.minimum(0.1 * kn * vtn, mu_fric * np.abs(fn)) * mask
F = (fn[:, :, None] * nij - ft[:, :, None] * tij).sum(axis=1)
# 颗粒-墙：底面 y=0 与两侧 x=0, x=box_w
ov = np.clip(R - pos[:, 1], 0.0, None)
F[:, 1] += kn * ov - eta_n * vel[:, 1] * (ov > 0)
ov = np.clip(R - pos[:, 0], 0.0, None)
F[:, 0] += kn * ov - eta_n * vel[:, 0] * (ov > 0)
ov = np.clip(R - (box_w - pos[:, 0]), 0.0, None)
F[:, 0] += -(kn * ov) + eta_n * vel[:, 0] * (ov > 0)
return F
# ---------------- 半隐式欧拉时间积分 ----------------
for step in range(nsteps):
F = contact_forces(pos, vel) + m * np.array([0.0, -g])
vel += (F / m) * dt
pos += vel * dt
if step % 10000 == 0:
print(f"t={step*dt:.2f}s KE={0.5*m*np.sum(vel**2):.3e} J "
f"mean_h={pos[:,1].mean()*1000:.1f} mm")
# ---------------- 结果统计 ----------------
ke = 0.5 * m * np.sum(vel**2)
d = np.linalg.norm(pos[:, None, :] - pos[None, :, :], axis=2)
n_contacts = int(((d < 2.02 * R) & (d > 0)).sum() // 2)
print(f"nt = {nsteps*dt:.2f} s 末态：")
print(f" 总动能 = {ke:.3e} J")
print(f" 平均高度 = {pos[:,1].mean()*1000:.1f} mm")
print(f" 接触对总数 = {n_contacts}")
print(f" 平均配位数 = {2*n_contacts/N:.2f}")稳定性校核：本例 $m \approx 1.05\times10^{-2}$ kg，$t_c \approx \pi\sqrt{m/2k_n} \approx 7\times10^{-4}$ s，$\Delta t = 10^{-5}$ s 约为 $t_c/70$，显式积分稳定。
工程案例分析
DEM 的典型工程应用是料仓与溜槽的卸料流动分析：颗粒级配、摩擦系数与仓壁倾角共同决定流型是”整体流”（mass flow）还是”漏斗流”（funnel flow），后者会造成死料区与偏析。另一重要场景是边坡与堆石体的稳定性评估——DEM 能自然再现颗粒体的剪胀、应变局部化与渐进破坏，这是连续介质模型难以刻画的。本文代码实测的堆积过程数据如下（玻璃珠参数见下表）：
| 参数 | 符号 | 数值 | 单位 |
|---|---|---|---|
| 颗粒半径 | $R$ | 10 | mm |
| 颗粒密度 | $\rho_p$ | 2500 | kg/m³ |
| 法向接触刚度 | $k_n$ | $1.0\times10^{5}$ | N/m |
| 恢复系数 | $e$ | 0.3 | — |
| 摩擦系数 | $\mu$ | 0.5 | — |
| 时间步长 | $\Delta t$ | $1.0\times10^{-5}$ | s |
模拟结果：颗粒群从平均高度 88 mm 自由下落，0.3 s 内动能从峰值 $6.6\times10^{-2}$ J 衰减至 $10^{-7}$ J 量级，最终形成平均高度 26.7 mm 的稳定料堆；末态接触对 29 个，平均配位数 2.32——与二维摩擦颗粒堆近等静压状态下配位数 2～3 的经验范围一致。值得强调的是，颗粒体系的宏观强度根植于微观接触网络（力链）的各向异性：强接触构成的力链骨架承担大部分外载，其方向随加载历史不断旋转重组，DEM 输出的逐接触力矢量可直接定量刻画这一细观机制，这是连续介质模型无法提供的信息维度。实践中，DEM 的可靠性取决于参数标定：宏观响应（堆积角、卸料流率）对 $\mu$ 与 $e$ 敏感，常需用实验堆积角反演接触参数，而非直接引用材料手册值。
结论与展望
本文从软球接触力学出发，建立了”弹簧-阻尼法向力 + 库仑摩擦切向力 + 显式积分”的 DEM 基本框架，并给出了可复现的二维堆积模拟。DEM 以接触参数换取对颗粒尺度物理的忠实描述，代价是微小时间步长带来的高昂计算量。当前发展方向包括：CFD-DEM 耦合模拟气固两相流（流化床、气力输送）、非球形颗粒（粘结球团、多面体、超二次曲面）建模、GPU 并行将颗粒规模推向 $10^8$ 量级，以及与实验数据同化结合的自动参数标定。对于任何涉及”颗粒”二字的工程问题，DEM 都是连接微观接触与宏观流动之间最直接的桥梁。
发表回复