LAMMPS SPH压力与人工黏性

在SPH里，粒子密度、状态方程和相互作用共同决定受力。本篇保持四点结构不动，用静止、相向和背向速度三个实际条件，拆开看Tait压力和人工黏性怎样进入原子力。这里的“原子”是LAMMPS记录名，物理上仍是粗粒子。

1. 先运行完整小例#

下载四个实际条件和可复跑绘图包，进入 lammps-sph，在CPU作业内执行：

bash run.sh

如果同名 raw已存在，脚本拒绝覆盖。重新练习请另解压目录。所有条件使用LAMMPS30Sep2026 CPU/oneAPI2024.2、一个线程；模板按56核整节点分配和计费，不在登录节点做生产模拟。

2. 完整相向速度输入#

cat inputs/approach.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/taitwater
pair_coeff * * 1000 10 1 0.03
variable vx atom "-0.2*(x-0.05)"
variable vy atom "-0.2*(y-0.05)"
variable vz atom "-0.2*(z-0.05)"
velocity all set v_vx v_vy v_vz units box
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

pair_coeff * * 1000 10 1 0.03在此风格中依次表示参考密度ρ₀、声速c₀、人工黏性系数和核半径。不同pair_style的参数位置不同，不能沿用上一页rhosum的解释。

本篇没有调用rhosum重新求密度，而是固定 sph/rho 1100，这样能单独观察状态方程和速度项。它与上一页的核密度算例是分开的条件。

3. Tait压力如何计算#

p = ρ₀ c₀² / 7 × [(ρ/ρ₀)^7 - 1]

本例ρ₀=1000 kg/m³、ρ=1100 kg/m³、c₀=10 m/s。图中压力由这个状态方程与原生密度计算，单位为Pa；不是普通原子MD的维里压力，也不用于推断真实水的状态方程。

静止条件最大粒子力约 8.46225 N，相向条件约 8.48245 N；背向条件与静止条件相同。每个条件的12个力分量均与独立公式逐项对照。

4. 为什么背向运动没有这项人工黏性#

相向条件取 v=-0.2(r-0.05)，背向条件取相反符号。对每对粒子，程序检查相对速度与相对位置的点积；只有接近时的负点积触发当前人工黏性项。因此即使两组速度大小相同，额外力也可以不同。

这个参数用于稳定相应SPH近似，不是本例测出的材料剪切黏度。如果研究真实流体黏度，还要选相应流体模型、边界与流动方式，积累统计并检查空间和时间步收敛。

head -n 14 raw/pressure/state.dump
head -n 14 raw/approach/state.dump
head -n 14 raw/recede/state.dump

5. 与密度核求和有什么区别#

此版本 sph/rhosum三维密度采用quadric核，而 sph/taitwater压力力采用Lucy核导数。不要把二者混成同一个能量函数，再用错误的能量梯度去判断受力。程序应按各自定义分别检查密度、压力和完整力向量；本包 verify.py也按这个方式实现。

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

完整粒子数组在 results/full-fields.npz，CSV包含坐标、核密度、静止压力力和额外相向力。本例四个 run 0没有推进时间，不能用来报告流量、扩散、温度平衡或黏度。准备做动态SPH前，应继续核对质量与分辨率、边界粒子、时间积分和CFL条件；官方SPH指南与Tait参数说明给出了相应入口。