LAMMPS MACHDYN 包 `pair_style smd/ulsph` 详解：更新拉格朗日 SPH 的密度、速度梯度与材料本构配置

• 科研
• 科学计算
• 高性能计算

smd/ulsph是 LAMMPS MACHDYN（Smooth Mach Dynamics）包中基于连续介质力学本构定律与更新拉格朗日（Updated Lagrangian）平滑粒子流体动力学（SPH）算法计算粒子相互作用的 pair 风格。它把压力（应力张量对角项）与剪切应力（非对角项）解耦计算，适合模拟水柱坍塌、流固耦合、颗粒流经复杂 STL 边界等大变形流体与固体问题。阅读本文后，你将掌握pair_style smd/ulsph的完整语法、六个风格级关键字的选择依据、pair_coeff的关键字化本构参数组织方式，以及如何搭配fix smd/integrate_ulsph完成可运行的 ULSPH 算例。

1. 概述：ULSPH 在 MACHDYN 中的定位

MACHDYN 包（源码位于 src/MACHDYN，README 见 src/MACHDYN/README）为 LAMMPS 提供两类 SPH 求解器：

• 传统 SPH 处理液体问题（如水柱、漏斗流动）；
• 带沙漏控制（hourglass control）的稳定无网格方法处理固体力学问题。

smd/ulsph属于前者，采用更新拉格朗日框架：粒子的参考构型随运动不断更新，通过连续介质本构关系计算相互作用。其核心思路是：

• 用状态方程（EOS）计算压力，即应力张量的对角项；
• 用材料强度模型 / 黏性模型计算剪切应力，即应力张量的非对角项。

这种"静水压力与偏应力分开指定"的设计（README 中明确说明 "can specify hydrostatic interactions independently from material strength models"），使得同一套 ULSPH 代码既能模拟无黏流体，也能模拟弹塑性固体。

2. 语法与三个必填关键字

2.1 命令语法

pair_style smd/ulsph args

args必须按顺序给出三个关键字，每个关键字以*开头（星号是关键字的一部分，不是通配符）：

keyword = *DENSITY_SUMMATION 或 *DENSITY_CONTINUITY 且 *VELOCITY_GRADIENT 或 *NO_VELOCITY_GRADIENT 且 *GRADIENT_CORRECTION 或 *NO_GRADIENT_CORRECTION

命令格式对参数个数有严格要求。在 pair_smd_ulsph.cpp 的settings()中，若参数个数不等于 3 会直接报错 "Illegal number of arguments for pair_style ulsph"，三个位置分别按固定关键字匹配，任何拼写偏差都会抛出 "Illegal settings keyword ..." 错误。

典型调用（文档示例）：

pair_style smd/ulsph *DENSITY_CONTINUITY *VELOCITY_GRADIENT *NO_GRADIENT_CORRECTION

2.2 第一组关键字：密度求解方式

实现层面（pair_smd_ulsph.cpp）：

• PreCompute_DensitySummation()（src/MACHDYN/pair_smd_ulsph.cpp#L127-L202）在*DENSITY_SUMMATION下被调用：每个粒子先计入自身贡献（h = 2.0*radius[i]处的 Poly6 核函数值乘自身质量），再遍历邻居累加wf * rmass[j]，最后以vfrac[i] = rmass[i] / rho[i]更新粒子体积分数。
• 若不使用密度求和，则密度变化量在时间积分中被累积（代码注释明确：otherwise, change in mass density is time-integrated）。

简单流体算例常选*DENSITY_CONTINUITY；*DENSITY_SUMMATION更稳健但每步都需要完整邻居循环重建密度，计算开销更大。

2.3 第二组关键字：是否计算速度梯度

速度梯度是强度模型与黏性模型的硬性前置条件。源码中，*STRENGTH_LINEAR_PLASTIC、*STRENGTH_LINEAR、*VISCOSITY_NEWTON分支都会先检查：

if (!velocity_gradient) error->all(FLERR, "A strength model was requested but *VELOCITY_GRADIENT is not set");

（见 pair_smd_ulsph.cpp#L1097-L1178）。因此只要pair_coeff里出现了强度或黏性模型卡片，pair_style就必须选*VELOCITY_GRADIENT，否则直接报错。

2.4 第三组关键字：核梯度修正

文档特别强调：*GRADIENT_CORRECTION可能导致严重的数值不稳定，一般流体模拟建议使用*NO_GRADIENT_CORRECTION。

实现层面，PreCompute()（src/MACHDYN/pair_smd_ulsph.cpp#L211-L351）在gradient_correction_flag为真时累加形状矩阵K[i] += jvol * (-g * dx^T)，随后用 SVD 伪逆求逆并做特征值限幅（LimitEigenvalues(K[i], 2.0)），最后对速度梯度 L 与变形梯度 F 施加修正。若关闭修正，K 直接置为单位阵（K[i].setIdentity()）。修正后的应力张量还会在AssembleStressTensor()末尾乘以 K（src/MACHDYN/pair_smd_ulsph.cpp#L815-L817）。

3.pair_coeff：关键字组织的本构参数卡片

与 LAMMPS 通常"纯浮点数 + 整数"的 pair_coeff 定义不同，smd/ulsph的pair_coeff用关键字组织参数集，分别对应粒子属性、材料本构模型与损伤模型，并以*END关键字收尾。文档推荐使用续行符&分行书写。

3.1 基本形式

pair_coeff i j *COMMON rho0 c0 Q1 Cp hg & *END

• i j为定义该 pair 的 LAMMPS 粒子类型对；
• *COMMON之后的 5 个参数是所有 ULSPH 粒子共享的基础属性；
• *END必须出现，用于终止pair_coeff行。

3.2*COMMON五个参数详解

源码 pair_smd_ulsph.cpp#L978-L998 严格校验*COMMON之后必须恰好跟 5 个参数（iNextKwd - ioffset != 5 + 1时报错 "expected 5 arguments following *COMMON"）。参数含义：

源码还会据此自动推导体模量：BULK_MODULUS = c0 * c0 * rho0（pair_smd_ulsph.cpp#L987-L988），因此*EOS_LINEAR不需要单独再给体模量参数。

3.3 状态方程（EOS）卡片

AssembleStressTensor()中的 switch（pair_smd_ulsph.cpp#L676-L697）展示了三种 EOS 的实际使用：Tait 与 Perfect Gas 会同时更新局部声速c0[i]（用于后续 CFL 时间步估计），Linear EOS 直接采用参考声速。若某种类型未指定任何 EOS，源码会报错 "must specify an EOS for particle type"。

3.4 强度模型与黏性模型卡片

约束规则（源码强制检查）：

• 一个粒子类型不能同时指定强度模型与黏性模型（cannot have both a strength and viscosity model，见 pair_smd_ulsph.cpp#L1271-L1272）；
• 使用强度/黏性模型必须开启*VELOCITY_GRADIENT。

强度模型的偏应力积分在 AssembleStressTensor 中完成：先由速度梯度分解出应变率 D 与自旋张量 W，偏应力率经过Jaumann 客观应力率（StressRateDevJaumann = stressRateDev - W*old + old*W）修正后时间积分，保证大转动下的客观性；塑性模型通过LinearPlasticStrength计算塑性应变增量并累加到eff_plastic_strain属性中。材料模型声明见 smd_material_models.h。

3.5 稳定化卡片（人工压力 / 人工应力）

人工压力在力计算循环中对正压力迹（张应力）按核函数比值加权削减（pair_smd_ulsph.cpp#L528-L536）；人工应力则对正特征值进行缩放（pair_smd_ulsph.cpp#L542-L556）。

3.6 交叉相互作用（cross interaction）

i与j可以不同，即允许不同粒子类型之间的 ULSPH 交叉作用。规则如下：

• 必须先给出i-i与j-j的 pair_coeff 行，然后才能定义交叉作用；
• 交叉行必须以*CROSS关键字标识：

  pair_coeff i j *CROSS

• 若i-i或j-j属性尚未定义，源码会报错 "properties of type ... have not yet been specified"（pair_smd_ulsph.cpp#L1289-L1293）；
• 交叉项的人工压力/人工应力取两端类型值的算术平均，否则为 0（pair_smd_ulsph.cpp#L1298-L1310）。

4. 一个完整的 ULSPH 水柱算例（源码示例）

仓库中 examples/PACKAGES/machdyn/fluid_structure_interaction/in.fluid_structure_interaction 提供了"水柱坍塌 + 多固体浮体 + 容器"的完整算例（units si，单位制 GPa/mm/ms），水用 ULSPH、固体用 TLSPH、接触用 Hertz 势：

pair_style hybrid/overlay smd/ulsph *DENSITY_CONTINUITY *VELOCITY_GRADIENT *NO_GRADIENT_CORRECTION & smd/tlsph smd/hertz ${contact_scale} pair_coeff 1 1 smd/ulsph *COMMON ${rho_water} ${c0} ${q1} ${Cp} 0 & *EOS_TAIT 7.0 & *END pair_coeff 2 2 smd/tlsph *COMMON ${rho_obj} ${E} ${nu} ${q1} ${q2} ${Hg} ${Cp} & *STRENGTH_LINEAR_PLASTIC ${sigma_yield} ${hardening_parameter} & *EOS_LINEAR & *END pair_coeff 3 3 none pair_coeff 1 2 smd/hertz ${contact_stiffness} pair_coeff 1 3 smd/hertz ${contact_stiffness} pair_coeff 2 3 smd/hertz ${contact_stiffness} pair_coeff 2 2 smd/hertz ${contact_stiffness}

其中水的参数为：rho_water = 1000、c0 = 10.0、q1 = 1.0、Cp = 1.0、EOS_TAIT指数 7.0、沙漏幅值 0。注意 ULSPH 水粒子只配了 EOS 与人工黏性（Q1），没有强度模型——这正是"压力与剪切应力解耦"的典型用法。

funnel_flow示例（examples/PACKAGES/machdyn/funnel_flow/in.funnel_flow）则展示了 ULSPH 与smd/tri_surface三角面片墙耦合、让流体流经 STL 复杂几何体的场景，其pair_style同样为smd/ulsph *DENSITY_CONTINUITY *VELOCITY_GRADIENT *NO_GRADIENT_CORRECTION，流体配*EOS_TAIT 2。

这两个示例的配套run脚本（*.sh）与输出日志可在同目录下查看，便于对照验证数值结果。

5. 配套时间积分：fix smd/integrate_ulsph

ULSPH 粒子必须用专门的积分器推进，即 fix smd/integrate_ulsph：

fix ID group-ID smd/integrate_ulsph keyword

可选的 keyword/value 对：

官方建议的adjust_radius典型取值为 factor=1.02、min_nn=15、max_nn=20（文档原文），水柱算例中实际使用adjust_radius 1.01 10 15，funnel_flow中使用adjust_radius 1.01 10 15。该 fix 要求atom_style同时含能量与体积属性（源码 fix_smd_integrate_ulsph.cpp#L50-L51 检查esph_flag与vfrac_flag），对应atom_style smd。

文档还提示：自 29Aug2024 版本起，该 fix 与会重映射速度的变形控制（如fix deform的remap v选项）不兼容。

6. 输出与诊断：ULSPH 配套 computes

smd/ulsph与其集成 fix 通过extract()（pair_smd_ulsph.cpp#L1463-L1485）暴露内部量，并有配套 compute 命令供 dump/thermo 输出：

• compute smd/ulsph/stress：应力张量；
• compute smd/ulsph/strain 与 compute smd/ulsph/strain_rate：应变与应变率；
• compute smd/ulsph/effm：有效剪切模量；
• compute smd/ulsph/num_neighs：每个粒子的邻居数，用于监测光滑核半径是否合适。

水柱算例中即用compute nn water smd/ulsph/num/neighs监控邻居数，配合fix smd/adjust_dt（动态时间步）与thermo_style输出f_dtfix dt观察 CFL 步长。ULSPH 的 CFL 时间步由声速主导：dtCFL = min(2*radius / p_wave_speed, dtCFL)，其中纵波速p_wave_speed = sqrt((K_eff + 4/3·G_eff) / rho)，黏性模型还会追加dtCFL = min(4·radius²·rho / mu, dtCFL)（pair_smd_ulsph.cpp#L800-L810）。

7. 限制与注意事项

依据文档与源码，smd/ulsph存在以下明确限制：

• 不做自动混合（No mixing is performed automatically），所有类型对必须显式给出pair_coeff；
• MACHDYN 全包不支持 restart 与 minimize，该 pair style 同样如此；fix smd/integrate_ulsph也没有输出量；
• rRESPA（多时间步）不适用于该 pair style；
• 该 fix/pair 属于 MACHDYN 包，仅在编译时启用该包后可用，构建方式见 Build_package；
• 力计算通过前向通信（comm_forward = 18，18 个 double：体积分数、声速、6 个对称应力分量、9 个变形梯度分量、等效塑性应变）把局部量同步到 ghost 原子（pair_smd_ulsph.cpp#L1388-L1457），并行下务必保证comm_modify vel yes与atom_modify map array等设置（见示例文件头部）；
• 核函数、Kernel_Cubic_Spline、Spiky 核等底层实现在 smd_kernels.h，力计算中人工黏性采用 Monaghan-Balsara 形式，速度梯度/变形梯度/应变增量累积全部基于 Eigen 张量运算。

从源码结构看，smd/ulsph与smd/tlsph（总拉格朗日 SPH）共享smd_material_models与smd_kernels，二者可在一个hybrid/overlay中混合使用——这正是在水柱算例中"流体 ULSPH + 固体 TLSPH + Hertz 接触"三类相互作用的组合方式，也是 MACHDYN 处理流固耦合问题的标准套路。

• 科研
• 科学计算
• 高性能计算