LAMMPS SPH粒子与密度求和

SPH把连续介质表示为带质量、密度、内能等字段的粒子。一个SPH粒子通常代表一小块流体，不是一个水分子。先弄清单位和密度怎样得到，再接压力、积分和边界，比直接套一份流体大算例更容易理解。

本篇用鸢算科技集群实际运行的四粒子静态例子，演示 sph/rhosum。它不推进物理时间，不用于测量水的密度、扩散或黏度。

1. 下载并运行四个小例子#

下载完整SPH输入、原生结果整理和绘图包，解压进入 lammps-sph，在分配好的CPU作业内执行：

bash run.sh

也可 sbatch job.slurm；CPU分区按56物理核整节点分配和计费，脚本本例只使用一个线程。采用LAMMPS30Sep2026 CPU版、oneAPI2024.2；绘图使用 sisl/0.16.4。本例没有使用GPU。

脚本依次运行密度核求和、固定密度压力、相向速度和背向速度四个条件，再输出完整CSV和两张图。

ls inputs
cat results/SUMMARY.json
head -n 5 results/particle-fields.csv

结构图坐标为cm，输入文件坐标为m，两者已明确换算；每个粒子质量0.01 kg。

2. 完整密度输入#

cat inputs/density.lmp

units si
atom_style sph
boundary f f f
region box block 0 0.1 0 0.1 0 0.1 units box
create_box 1 box
create_atoms 1 single 0.03 0.04 0.05 units box
create_atoms 1 single 0.05 0.04 0.05 units box
create_atoms 1 single 0.04 0.06 0.05 units box
create_atoms 1 single 0.04 0.05 0.07 units box
mass 1 0.01
set group all sph/rho 1100
set group all sph/e 0
set group all sph/cv 1
comm_modify vel yes
neighbor 0.002 bin
neigh_modify delay 0 every 1 check yes
velocity all set 0 0 0 units box
pair_style sph/rhosum 1
pair_coeff * * 0.03
compute rho all sph/rho/atom
compute u all sph/e/atom
fix int all sph
timestep 0.00001
thermo 1
thermo_style custom step atoms
thermo_modify format float %.17g
run 0
write_dump all custom state.dump id type x y z vx vy vz fx fy fz c_rho c_u modify sort id format float %.17g

units si使长度为m、质量为kg、时间为s、力为N。不要从普通金属原子模拟里直接复制Å或原子质量数过来。

atom_style sph提供密度、内能等字段。开始时的 sph/rho 1100是初值；本例实际调用 sph/rhosum 1后，会按结构和质量重新求和，结果不是继续保持1100。

pair_coeff * * 0.03里的0.03是此风格的核截断半径h，单位m。它不是普通LJ势参数，也不是“最近3个邻居”。

3. 核求和包括自己#

密度为 ρi = Σj mj W(rij,h)，其中包括i自身的 mi W(0,h)。当前版本三维实现使用：

W(r,h) = 3465 / (512πh³) × (1-r²/h²)^4，r<h
W(r,h) = 0，r≥h

四点原生密度约 848.678–914.831 kg/m³；独立逐粒子求和得到相同结果。这是稀疏四点与截断核的结果，不能据此得出“水密度只有这个数值”的结论。实际流体计算还要考虑粒子间距、边界、核覆盖和空间分辨率。

4. 读取原生输出#

head -n 14 raw/density/state.dump

c_rho是 compute sph/rho/atom得到的原生密度；c_u为内能字段。该密度例子只求核密度，没有压力相互作用，因此受力为零。

输入中的 run 0只初始化邻居和计算当前字段，不做MD积分。虽然列出了 fix sph和时间步，粒子并没有真的运动，更不能把零步输出当作平衡流体轨迹。

5. 后处理与下一步#

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

这两步重读已经完成的原生结果。密度核求和与压力力使用的核导数在此版本并不相同，接下来阅读SPH压力与人工黏性。更多核求和参数见官方sph/rhosum说明。