LAMMPS颗粒接触与球堆下落

颗粒模型的粒子代表有限尺寸球，不是一个原子。半径决定何时接触，密度和体积决定质量，接触刚度、阻尼与摩擦共同决定碰撞行为。这里先用两个球把法向力算清楚，再看27球下落。

1. 运行流程#

下载完整输入、运行脚本和绘图包，解压进入 lammps-mesoscale，在已分配的CPU作业内执行：

bash run.sh

也可 sbatch job.slurm 提交。CPU分区按56核整节点分配和计费；本例只用一个线程，使用 lammps/30Sep2026-8173m 的普通CPU样式、oneAPI2024.2构建，绘图用 sisl/0.16.4。本例不是GPU加速测试。CPU队列处于维护时请待恢复后提交，不在登录节点做生产模拟。

脚本按输入逐项运行，在 raw/ 保留原生输出，在 results/ 输出全数组和汇总，在 plots/ 生成本次图片。三个教程共用一个包，完整运行一次后只读结果和重画即可。

2. 先看完整的两球输入#

cat inputs/granular-contact.lmp

units si
atom_style sphere
boundary f f f
region box block -0.03 0.03 -0.03 0.03 -0.03 0.03 units box
create_box 1 box
create_atoms 1 single 0 0 0 units box
create_atoms 1 single 0.009 0 0 units box
set group all diameter 0.01
set group all density 2500
group left id 1
group right id 2
velocity left set 0 0 0 units box
velocity right set 0 0 0 units box
comm_modify vel yes
pair_style gran/hooke 1000 0 20 0 0.5 0
pair_coeff * *
neighbor 0.001 bin
thermo 1
thermo_modify norm no
run 0
write_dump all custom state.dump id type x y z vx vy vz fx fy fz radius mass modify sort id format float %.17g

units si 对应m、kg、s、N；直径0.01m、密度2500kg/m³给出单球质量约 0.001308997kg。两中心相距0.009m，重叠0.001m，法向刚度1000N/m，因此静态排斥力是 1N。

comm_modify vel yes 让相互作用需要的速度传给相邻域。run 0 只计算当前状态的力，不是碰撞轨迹。

3. 比较间隙、接近和远离#

cat results/SUMMARY.json
cat raw/granular-approach/state.dump

阻尼项与相对法向速度以及有效质量有关，接近时增强阻力、分离时降低排斥力。这里关闭切向弹性/阻尼以便单独观察法向项；不能据这四个状态推断摩擦系数或恢复系数。

4. 27球如何下落并与地面接触#

cat inputs/granular-pile.lmp

该输入改用 gran/hooke/history、fix nve/sphere、重力和 fix wall/gran，对27个球运行15000步，每步10⁻⁵s，总时长0.15s。球面接触允许弹性重叠；刚度增大时要检查积分步长。

fix integrate all nve/sphere
fix grav all gravity 9.81 vector 0 0 -1
fix floor all wall/gran hooke/history 10000 2857 100 50 0.5 1 zplane 0 NULL
timestep 0.00001
run 15000

末帧最低球心约 9.419mm；球仍会反弹，本例没有把0.15s末态称为完全沉降的平衡球堆。接触地面前的全部球位置和速度与自由落体关系相符。

5. 复查并绘图#

module load sisl/0.16.4
python -B verify.py --reference
python -B plot.py

这些命令只处理已有输出。结构图用球心表示位置，标记大小不表示严格按比例画出的球半径。真正做堆积密度或流动问题，还需要更大体系、边界条件和足够的沉降/采样时间。