科研科学计算高性能计算【免费下载链接】lammpsPublic development project of the LAMMPS MD software package项目地址https://gitcode.com/gh_mirrors/la/lammps点击查看免费下载LAMMPS 的 GRANULAR 包提供了一组经典的颗粒接触力场用于模拟球形颗粒之间的碰撞、摩擦与阻尼。本文围绕pair_style gran/hooke、gran/hooke/history、gran/hertz/history三个 pair 风格系统讲解其命令语法、力模型公式、系数单位与换算关系、库仑摩擦准则、历史剪切存储机制并结合仓库源码src/GRANULAR剖析底层实现最后给出构建要求、restart 限制与官方测试用例帮助读者在离散元DEM仿真中正确配置与调参。概述三种风格与适用场景gran系列 pair 风格根据 Brilliantov 等人1996、Silbert2001和 Zhang 与 Makse2005的文献为两个半径分别为 R_i、R_j 的球形颗粒计算接触摩擦力。当两颗粒中心距离 r 小于接触距离 d R_i R_j 时产生作用力当 r d 时两颗粒之间无作用力。三种风格的核心区别在于法向接触力的形式以及是否计入切向位移历史gran/hooke法向接触力是重叠量的线性函数Hookean 线性弹簧模型不含切向历史剪切项gran/hooke/history在线性法向模型基础上计入切向位移的历史效应历史剪切力gran/hertz/history法向接触力与重叠面积成正比Hertz 非线性模型含粒径相关的前置因子同时计入切向历史剪切力。其中gran/hooke/history是gran/hertz/history的基类见源码 pair_gran_hertz_history.cppPairGranHertzHistory直接继承自PairGranHookeHistory三者共享同一套参数解析与邻域遍历框架。命令语法与参数详解三种风格使用完全相同的命令行语法pair_style style Kn Kt gamma_n gamma_t xmu dampflag keyword参数含义单位 / 取值范围stylegran/hooke、gran/hooke/history或gran/hertz/history—Kn法向排斥弹性常数力/距离Hookean或压强Hertzian见下文单位讨论Kt切向接触弹性常数力/距离Hookean或压强Hertziangamma_n法向碰撞阻尼系数1/时间Hookean或 1/(时间·距离)Hertziangamma_t切向碰撞阻尼系数1/时间Hookean或 1/(时间·距离)Hertzianxmu静态屈服准则库仑摩擦系数无量纲0.0 ~ 1.0e4dampflag0 或 1表示排除或包含切向阻尼力0 / 1keywordlimit_damping无取值限制阻尼避免产生吸引作用limit_damping关键字没有附加数值当两颗粒处于接触中且正在相互远离时阻尼项可能造成等效的吸引力启用该关键字后一旦法向合力出现等效吸引即法向系数 ccel 0就把法向力置零。官方示例文档给出三个可直接使用的配置示例pair_style gran/hooke/history 200000.0 NULL 50.0 NULL 0.5 1 pair_style gran/hooke 200000.0 70000.0 50.0 30.0 0.5 0 pair_style gran/hooke 200000.0 70000.0 50.0 30.0 0.5 0 limit_damping注意第一个示例中Kt和gamma_t使用了NULL当Kt为 NULL 时采用默认值Kt 2/7 * Kn当gamma_t为 NULL 时采用默认值gamma_t 1/2 * gamma_n。实测示例unittest/force-styles/tests/atomic-pair-gran_hooke_history.yaml也验证了pair_style gran/hooke/history 1000.0 NULL 50.0 NULL 0.5 1 limit_damping这种 NULL 组合的合法性。力模型公式法向项与切向项两种 Hookean 风格使用如下公式F_hk (k_n δ n_ij − m_eff γ_n v_n) − (k_t Δs_t m_eff γ_t v_t)Hertzian 风格使用如下公式F_hz √δ √(R_i R_j / (R_i R_j)) F_hk即 Hertz 模型是 Hooke 模型整体乘以一个与重叠量 δ 和两颗粒半径相关的前置因子√δ · √(R_iR_j/(R_iR_j))。两个方程中第一个括号为法向力含接触力与阻尼力两项第二个括号为切向力含剪切力与阻尼力两项剪切力是历史效应代表两颗粒在整个接触持续时间内累积的切向相对位移仅存在于hooke/history与hertz/history中hooke不包含此项切向阻尼力在三个风格中均存在但仅当dampflag 1时才计入。各物理量定义符号含义δ d − r两颗粒重叠距离K_n法向接触弹性常数K_t切向接触弹性常数γ_n法向接触黏弹性阻尼常数γ_t切向接触黏弹性阻尼常数m_eff M_i M_j / (M_i M_j)两颗粒质量 M_i、M_j的有效质量Δs_t两颗粒间切向位移矢量按摩擦屈服准则截断n_ij沿两颗粒中心连线的单位矢量v_n两颗粒相对速度的法向分量v_t两颗粒相对速度的切向分量系数单位的两种解读Hookean 与 Hertzian 的本质差异文档特别强调四个系数K_n、K_t、γ_n、γ_t在 Hookean 与 Hertzian 方程中的量纲解释完全不同这是配置时最容易出错的地方。Hookean 模型法向回推push-back力是重叠距离的线性函数因此 K_n 的单位是力/距离。该回推力与颗粒绝对尺寸无关单分散情形、与两颗粒相对尺寸无关多分散情形其余各项同理γ_n 单位为1/时间K_t 单位为力/距离γ_t 单位为1/时间。在 Hookean 情形下K_n 可视为线弹性弹簧常数。Hertzian 模型法向回推力与两颗粒的重叠面积成正比是重叠距离的非线性函数因此 K_n 的单位是力/面积压强。粒径的绝对效应单分散与相对效应多分散都体现在半径相关的前置因子中将这些前置因子代入方程其余各项后γ_n 单位为1/(时间·距离)K_t 单位为压强γ_t 单位为1/(时间·距离)。从材料参数换算 Hertzian 系数据 Zhang Makse 2005 文献K_n 4G / (3(1−ν))其中 G E / (2(1ν)) 为剪切模量E 为杨氏模量ν 为泊松比 K_t 4G / (2−ν)文档特别更正早期手册版本曾错误地把 K_t 写为 8G/(2−ν)。因此 Hertzian 情形下 K_n、K_t 可以直接对应被模拟材料的真实物性参数而 Hookean 情形下则必须为模型中颗粒的绝对尺寸挑选合适的弹簧常数。由于 Hookean 不反映相对粒径效应Hookean 风格可能不适合多分散体系多分散体系优先选择gran/hertz/history。9Jan09 版本变更警告LAMMPS 2009 年 1 月 9 日之前的版本中Hertzian 方程不含√(r_i r_j / (r_i r_j))项对多分散体系不够准确。该项引入后即使对单分散体系要复现旧版本的计算结果也必须重新标定四个系数——因为对单分散体系而言√(r_i r_j/(r_ir_j))是常数因子等价于把四个系数整体缩放。文档给出的换算示例对常见的颗粒直径 1的单分散体系四个系数现在都应设置为旧版本的 2 倍。xmu 库仑屈服准则与 limit_dampingxmu是切向力上限的静态屈服准则通过库仑准则 Ft xmu·Fn 生效Ft、Fn 分别为公式中的总切向与法向力分量。Hookean 情形下两颗粒间的切向力按照切向弹簧 阻尼器模型增长直到 Ft/Fn xmu之后被钳制在 Ft Fn·xmu直到颗粒失去接触Hertzian 情形下原理类似只是弹簧不再是线性的。xmu通常应取 0.0 到 1.0 之间的分数值LAMMPS 允许最大值高达 1.0e4用于模拟能承受极大切向力的体系例如需要冻结相对切向滑移的场景。limit_damping关键字针对两颗粒接触且相互远离时阻尼项可能产生等效吸引力的问题启用后若法向合力出现等效吸引则将该法向分量清零源码中对应if (limit_damping (ccel 0.0)) ccel 0.0;见 pair_gran_hooke.cpp。官方测试 unittest/commands/test_restart_settings.cpp 也验证了gran/hooke/history ... limit_damping在 restart 场景下的可用性。有效质量与刚性体处理有效质量 m_eff 按两个孤立颗粒的公式M_i M_j / (M_i M_j)计算。若任一颗粒属于某个刚性体其质量替换为该刚性体的总质量这是通过搜索fix rigid命令及其变体实现的。从源码看pair_gran_hooke.cpp实现细节是若检测到fix_rigid则通过fix_rigid-extract(body, ...)与extract(masstotal, ...)提取每个原子所属的刚体序号与刚体总质量缓存到mass_rigid[]数组再通过comm-forward_comm(this)把刚体质量同步到 ghost 原子。另外若颗粒属于冻结组fix freeze见 pair_granular.cpp 对 freeze 的检测m_eff 退化为另一颗粒的质量即meff mji 冻结时或meff mij 冻结时这在 pair_gran_hooke.cpp 中直接体现。pair_coeff 的全局与混合用法与大多数 pair 风格不同gran 系列没有针对每个原子类型对的额外系数——所有设置都是全局的通过pair_style命令完成。但你仍然必须为所有涉及颗粒的原子类型对给出pair_coeff命令例如pair_coeff * *当所有原子都通过上述某个 gran 风格相互作用时使用* *。若 gran 风格作为pair_style hybrid的子风格使用则可以在pair_coeff中指定具体原子类型以决定哪些类型组合通过颗粒势相互作用。源码级实现剖析力计算主循环以 pair_gran_hooke.cpp 的compute()为例每对接触颗粒的计算流程为计算相对位移与距离delx/dely/delz、rsq接触判据rsq radsum*radsum把相对平移速度分解为法向分量vnnr投影到连线方向与切向分量计算相对转动速度wr由角速度 omega 与半径加权合成计算法向力ccel kn*(radsum-r)*rinv - damp其中damp meff*gamman*vnnr*rsqinv为法向阻尼项计算切向阻尼力并按库仑准则截断fn xmu*fabs(ccel*r)fs meff*gammat*vrel最终ft MIN(fn,fs)/vrel——这正是 Ft xmu·Fn 的逐对实现叠加力并计算力矩tor1/tor2/tor3最后按newton_pair规则把反作用力与力矩施加到 j 颗粒。历史剪切量的存储与更新hooke/history与hertz/history通过fix neigh/history机制源码 pair_gran_hooke_history.cpp 中构造占位的FixDummy随后在init_style中替换为FixNeighHistory为每个接触对存储 3 个历史量size_history 3。每个时间步内剪切位移按切向相对速度累积shear[] vtr[]*dt由于颗粒会转动剪切位移需做旋转变换去除沿连线方向的分量rsht shear·del · rsqinv投影后扣除见 pair_gran_hooke_history.cpp剪切力为kt*shear[]与切向阻尼之和当剪切力幅值fs fn xmu*|ccel*r|时对剪切位移与剪切力整体按fn/fs比例重新标定rescaling实现库仑摩擦滑移见 pair_gran_hertz_history.cpp。值得注意的是颗粒一旦脱离接触rsq radsum*radsumtouch标志被清零、剪切历史被归零见 pair_gran_hertz_history.cpp保证历史量只反映本次接触过程。Hertzian 的 polydispersity 因子hertz/history与hooke/history的唯一差别在于法向项乘以polyhertz sqrt((radsum-r)*radi*radj/radsum)见 pair_gran_hertz_history.cpp该因子同时作用于切向力fs1 -polyhertz * (kt*shear[0] ...)这与文档公式中整个 F_hk 乘以 √δ√(R_iR_j/(R_iR_j))完全一致也解释了为何旧版本需要重新标定四个系数。与新一代 pair_style granular 的关系文档明确指出pair_style granular见 doc/src/pair_granular.rst在法向、切向、滚动、扭转等子模型上大大扩展了pair gran/*系列的能力。现代统一入口 pair_granular.cpp 通过pair_coeff按类型对组合选择子模型normal/tangential/damping/rolling/twisting/heat实现同一套体系内不同材料类型使用不同接触模型。例如 examples/granular/in.chute 使用pair_style granular pair_coeff * * hertz 1 0.5 tangential linear_history 1.0 0.5 0.1而旧式gran/hooke等风格则在需要全局统一参数、快速上手的场景下依然简洁有效例如 examples/granular/in.restitution 中评估颗粒碰撞恢复系数的实验即使用pair_style granularhertz/material子模型。从源码结构看pair_gran_hooke、pair_gran_hooke_history、pair_gran_hertz_history属于相对早期的经典实现其计算内核法向弹簧-阻尼、切向历史剪切、库仑截断与 granular 的默认子模型逻辑一脉相承。加速变体三种风格均提供加速变体需要相应构建支持gran/hookegran/hooke/kkKOKKOS、gran/hooke/ompOpenMPgran/hooke/historygran/hooke/history/omp、gran/hooke/history/kkgran/hertz/historygran/hertz/history/kk、gran/hertz/history/omp。mixing / shift / table / tail、restart 与 rRESPA 行为pair_modify的 mix、shift、table、tail 选项对 gran 系列不适用这些 pair 风格会把自身信息写入二进制 restart 文件因此读取 restart 文件的输入脚本无需重复指定pair_style命令在run_style respa中gran 系列只能通过pair关键字使用不支持inner、middle、outer关键字这些风格的single()函数返回 0.0 作为成对相互作用能量——因为这类耗散势不守恒能量它只返回成对作用力的法向分量。single() 返回的 10 个额外量single()还会计算 10 个额外的成对量可通过compute pair/local命令以p1、p2、…、p10访问索引内容p1–p3作用于颗粒 I 的切向力三个分量p4该切向力的幅值p5–p7法向相对速度的三个分量沿两球心连线方向p8–p10切向相对速度的三个分量对应源码见 pair_gran_hooke.cppsvector[0..2]为切向力分量svector[3]为其模svector[4..6]为法向相对速度svector[7..9]为切向相对速度。构建与运行限制所有 gran 风格都属于GRANULAR 包必须使用该包构建 LAMMPS 才能启用见 Build package 页面包源码位于 src/GRANULAR这些风格要求原子存储力矩与角速度omega由 atom_style 定义并要求存储每个颗粒的半径——atom_style sphere同时满足这两点必须使用comm_modify vel yes命令使 ghost 原子也保存速度信息使用read_restart时不能精确复现运行结果统计上相似即可因为力依赖原子速度但单个接触的累积历史量会写入 restart 文件历史量不会写入 data 文件因此通过read_data重新加载体系时由于丢失接触历史力可能有显著差异。测试与示例验证仓库提供了针对三种风格的回归测试用于校验力与应力计算的确定性unittest/force-styles/tests/atomic-pair-gran_hooke.yamlpair_style gran/hooke 1000.0 350.0 50.0 30.0 0.5 1 limit_dampingunittest/force-styles/tests/atomic-pair-gran_hooke_history.yamlunittest/force-styles/tests/atomic-pair-gran_hertz_history.yaml测试前置命令见 atomic-pair-gran_hooke.yaml给出了运行 gran 风格的最小骨架units lj、atom_style sphere、comm_modify vel yes、atom_modify map array以及set group all diameter、set group all mass颗粒质量与半径是必需的。examples/granular 目录下的in.restitution双球对撞测恢复系数与in.chute二维溜槽流动可作为直接运行的实战起点。相关命令与默认值相关命令pair_coeff。默认值none即所有系数都必须显式给定Kt、gamma_t用 NULL 时分别取2/7 Kn、1/2 gamma_n的默认值。参考文献Brilliantov, Spahn, Hertzsch, Poschel, Phys Rev E, 53, p 5382-5392 (1996)Silbert, Ertas, Grest, Halsey, Levine, Plimpton, Phys Rev E, 64, p 051302 (2001)Zhang and Makse, Phys Rev E, 72, p 011301 (2005)赞分享科研科学计算高性能计算【免费下载链接】lammpsPublic development project of the LAMMPS MD software package项目地址https://gitcode.com/gh_mirrors/la/lammps点击查看免费下载相关推荐LAMMPS fix wall/gran/region 命令详解用任意几何 region 构造颗粒摩擦壁LAMMPS fix wall/gran/region 命令详解用任意几何 region 构造颗粒摩擦壁 导读 fix wall/gran/region 是科研科学计算高性能计算LAMMPS pair_style buck 系列Buckingham 势的语法、参数与源码实现解析LAMMPS pair_style buck 系列Buckingham 势的语法、参数与源码实现解析 导读 本文基于 LAMMPS 官方文档 doc/src/科研科学计算高性能计算LAMMPS pair_style colloid 胶体相互作用势语法、参数与源码实现详解LAMMPS pair_style colloid 胶体相互作用势语法、参数与源码实现详解 导读 本文围绕 LAMMPS 的 pair_style collo科研科学计算高性能计算上一篇WinUtil告别Windows繁琐设置一键搞定系统优化与软件安装下一篇immich-go 更优资产选择算法解析文件去重、日期比对与画质择优的上传判定机制创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考 
长河君 · 2026/10/16 1:55:20 · 阅读 1.2k LAMMPS 颗粒接触力场 pair_style gran/hooke 与 gran/hertz/history 全解析：从语法参数到源码实现
科研科学计算高性能计算【免费下载链接】lammpsPublic development project of the LAMMPS MD software package项目地址https://gitcode.com/gh_mirrors/la/lammps点击查看免费下载LAMMPS 的 GRANULAR 包提供了一组经典的颗粒接触力场用于模拟球形颗粒之间的碰撞、摩擦与阻尼。本文围绕pair_style gran/hooke、gran/hooke/history、gran/hertz/history三个 pair 风格系统讲解其命令语法、力模型公式、系数单位与换算关系、库仑摩擦准则、历史剪切存储机制并结合仓库源码src/GRANULAR剖析底层实现最后给出构建要求、restart 限制与官方测试用例帮助读者在离散元DEM仿真中正确配置与调参。概述三种风格与适用场景gran系列 pair 风格根据 Brilliantov 等人1996、Silbert2001和 Zhang 与 Makse2005的文献为两个半径分别为 R_i、R_j 的球形颗粒计算接触摩擦力。当两颗粒中心距离 r 小于接触距离 d R_i R_j 时产生作用力当 r d 时两颗粒之间无作用力。三种风格的核心区别在于法向接触力的形式以及是否计入切向位移历史gran/hooke法向接触力是重叠量的线性函数Hookean 线性弹簧模型不含切向历史剪切项gran/hooke/history在线性法向模型基础上计入切向位移的历史效应历史剪切力gran/hertz/history法向接触力与重叠面积成正比Hertz 非线性模型含粒径相关的前置因子同时计入切向历史剪切力。其中gran/hooke/history是gran/hertz/history的基类见源码 pair_gran_hertz_history.cppPairGranHertzHistory直接继承自PairGranHookeHistory三者共享同一套参数解析与邻域遍历框架。命令语法与参数详解三种风格使用完全相同的命令行语法pair_style style Kn Kt gamma_n gamma_t xmu dampflag keyword参数含义单位 / 取值范围stylegran/hooke、gran/hooke/history或gran/hertz/history—Kn法向排斥弹性常数力/距离Hookean或压强Hertzian见下文单位讨论Kt切向接触弹性常数力/距离Hookean或压强Hertziangamma_n法向碰撞阻尼系数1/时间Hookean或 1/(时间·距离)Hertziangamma_t切向碰撞阻尼系数1/时间Hookean或 1/(时间·距离)Hertzianxmu静态屈服准则库仑摩擦系数无量纲0.0 ~ 1.0e4dampflag0 或 1表示排除或包含切向阻尼力0 / 1keywordlimit_damping无取值限制阻尼避免产生吸引作用limit_damping关键字没有附加数值当两颗粒处于接触中且正在相互远离时阻尼项可能造成等效的吸引力启用该关键字后一旦法向合力出现等效吸引即法向系数 ccel 0就把法向力置零。官方示例文档给出三个可直接使用的配置示例pair_style gran/hooke/history 200000.0 NULL 50.0 NULL 0.5 1 pair_style gran/hooke 200000.0 70000.0 50.0 30.0 0.5 0 pair_style gran/hooke 200000.0 70000.0 50.0 30.0 0.5 0 limit_damping注意第一个示例中Kt和gamma_t使用了NULL当Kt为 NULL 时采用默认值Kt 2/7 * Kn当gamma_t为 NULL 时采用默认值gamma_t 1/2 * gamma_n。实测示例unittest/force-styles/tests/atomic-pair-gran_hooke_history.yaml也验证了pair_style gran/hooke/history 1000.0 NULL 50.0 NULL 0.5 1 limit_damping这种 NULL 组合的合法性。力模型公式法向项与切向项两种 Hookean 风格使用如下公式F_hk (k_n δ n_ij − m_eff γ_n v_n) − (k_t Δs_t m_eff γ_t v_t)Hertzian 风格使用如下公式F_hz √δ √(R_i R_j / (R_i R_j)) F_hk即 Hertz 模型是 Hooke 模型整体乘以一个与重叠量 δ 和两颗粒半径相关的前置因子√δ · √(R_iR_j/(R_iR_j))。两个方程中第一个括号为法向力含接触力与阻尼力两项第二个括号为切向力含剪切力与阻尼力两项剪切力是历史效应代表两颗粒在整个接触持续时间内累积的切向相对位移仅存在于hooke/history与hertz/history中hooke不包含此项切向阻尼力在三个风格中均存在但仅当dampflag 1时才计入。各物理量定义符号含义δ d − r两颗粒重叠距离K_n法向接触弹性常数K_t切向接触弹性常数γ_n法向接触黏弹性阻尼常数γ_t切向接触黏弹性阻尼常数m_eff M_i M_j / (M_i M_j)两颗粒质量 M_i、M_j的有效质量Δs_t两颗粒间切向位移矢量按摩擦屈服准则截断n_ij沿两颗粒中心连线的单位矢量v_n两颗粒相对速度的法向分量v_t两颗粒相对速度的切向分量系数单位的两种解读Hookean 与 Hertzian 的本质差异文档特别强调四个系数K_n、K_t、γ_n、γ_t在 Hookean 与 Hertzian 方程中的量纲解释完全不同这是配置时最容易出错的地方。Hookean 模型法向回推push-back力是重叠距离的线性函数因此 K_n 的单位是力/距离。该回推力与颗粒绝对尺寸无关单分散情形、与两颗粒相对尺寸无关多分散情形其余各项同理γ_n 单位为1/时间K_t 单位为力/距离γ_t 单位为1/时间。在 Hookean 情形下K_n 可视为线弹性弹簧常数。Hertzian 模型法向回推力与两颗粒的重叠面积成正比是重叠距离的非线性函数因此 K_n 的单位是力/面积压强。粒径的绝对效应单分散与相对效应多分散都体现在半径相关的前置因子中将这些前置因子代入方程其余各项后γ_n 单位为1/(时间·距离)K_t 单位为压强γ_t 单位为1/(时间·距离)。从材料参数换算 Hertzian 系数据 Zhang Makse 2005 文献K_n 4G / (3(1−ν))其中 G E / (2(1ν)) 为剪切模量E 为杨氏模量ν 为泊松比 K_t 4G / (2−ν)文档特别更正早期手册版本曾错误地把 K_t 写为 8G/(2−ν)。因此 Hertzian 情形下 K_n、K_t 可以直接对应被模拟材料的真实物性参数而 Hookean 情形下则必须为模型中颗粒的绝对尺寸挑选合适的弹簧常数。由于 Hookean 不反映相对粒径效应Hookean 风格可能不适合多分散体系多分散体系优先选择gran/hertz/history。9Jan09 版本变更警告LAMMPS 2009 年 1 月 9 日之前的版本中Hertzian 方程不含√(r_i r_j / (r_i r_j))项对多分散体系不够准确。该项引入后即使对单分散体系要复现旧版本的计算结果也必须重新标定四个系数——因为对单分散体系而言√(r_i r_j/(r_ir_j))是常数因子等价于把四个系数整体缩放。文档给出的换算示例对常见的颗粒直径 1的单分散体系四个系数现在都应设置为旧版本的 2 倍。xmu 库仑屈服准则与 limit_dampingxmu是切向力上限的静态屈服准则通过库仑准则 Ft xmu·Fn 生效Ft、Fn 分别为公式中的总切向与法向力分量。Hookean 情形下两颗粒间的切向力按照切向弹簧 阻尼器模型增长直到 Ft/Fn xmu之后被钳制在 Ft Fn·xmu直到颗粒失去接触Hertzian 情形下原理类似只是弹簧不再是线性的。xmu通常应取 0.0 到 1.0 之间的分数值LAMMPS 允许最大值高达 1.0e4用于模拟能承受极大切向力的体系例如需要冻结相对切向滑移的场景。limit_damping关键字针对两颗粒接触且相互远离时阻尼项可能产生等效吸引力的问题启用后若法向合力出现等效吸引则将该法向分量清零源码中对应if (limit_damping (ccel 0.0)) ccel 0.0;见 pair_gran_hooke.cpp。官方测试 unittest/commands/test_restart_settings.cpp 也验证了gran/hooke/history ... limit_damping在 restart 场景下的可用性。有效质量与刚性体处理有效质量 m_eff 按两个孤立颗粒的公式M_i M_j / (M_i M_j)计算。若任一颗粒属于某个刚性体其质量替换为该刚性体的总质量这是通过搜索fix rigid命令及其变体实现的。从源码看pair_gran_hooke.cpp实现细节是若检测到fix_rigid则通过fix_rigid-extract(body, ...)与extract(masstotal, ...)提取每个原子所属的刚体序号与刚体总质量缓存到mass_rigid[]数组再通过comm-forward_comm(this)把刚体质量同步到 ghost 原子。另外若颗粒属于冻结组fix freeze见 pair_granular.cpp 对 freeze 的检测m_eff 退化为另一颗粒的质量即meff mji 冻结时或meff mij 冻结时这在 pair_gran_hooke.cpp 中直接体现。pair_coeff 的全局与混合用法与大多数 pair 风格不同gran 系列没有针对每个原子类型对的额外系数——所有设置都是全局的通过pair_style命令完成。但你仍然必须为所有涉及颗粒的原子类型对给出pair_coeff命令例如pair_coeff * *当所有原子都通过上述某个 gran 风格相互作用时使用* *。若 gran 风格作为pair_style hybrid的子风格使用则可以在pair_coeff中指定具体原子类型以决定哪些类型组合通过颗粒势相互作用。源码级实现剖析力计算主循环以 pair_gran_hooke.cpp 的compute()为例每对接触颗粒的计算流程为计算相对位移与距离delx/dely/delz、rsq接触判据rsq radsum*radsum把相对平移速度分解为法向分量vnnr投影到连线方向与切向分量计算相对转动速度wr由角速度 omega 与半径加权合成计算法向力ccel kn*(radsum-r)*rinv - damp其中damp meff*gamman*vnnr*rsqinv为法向阻尼项计算切向阻尼力并按库仑准则截断fn xmu*fabs(ccel*r)fs meff*gammat*vrel最终ft MIN(fn,fs)/vrel——这正是 Ft xmu·Fn 的逐对实现叠加力并计算力矩tor1/tor2/tor3最后按newton_pair规则把反作用力与力矩施加到 j 颗粒。历史剪切量的存储与更新hooke/history与hertz/history通过fix neigh/history机制源码 pair_gran_hooke_history.cpp 中构造占位的FixDummy随后在init_style中替换为FixNeighHistory为每个接触对存储 3 个历史量size_history 3。每个时间步内剪切位移按切向相对速度累积shear[] vtr[]*dt由于颗粒会转动剪切位移需做旋转变换去除沿连线方向的分量rsht shear·del · rsqinv投影后扣除见 pair_gran_hooke_history.cpp剪切力为kt*shear[]与切向阻尼之和当剪切力幅值fs fn xmu*|ccel*r|时对剪切位移与剪切力整体按fn/fs比例重新标定rescaling实现库仑摩擦滑移见 pair_gran_hertz_history.cpp。值得注意的是颗粒一旦脱离接触rsq radsum*radsumtouch标志被清零、剪切历史被归零见 pair_gran_hertz_history.cpp保证历史量只反映本次接触过程。Hertzian 的 polydispersity 因子hertz/history与hooke/history的唯一差别在于法向项乘以polyhertz sqrt((radsum-r)*radi*radj/radsum)见 pair_gran_hertz_history.cpp该因子同时作用于切向力fs1 -polyhertz * (kt*shear[0] ...)这与文档公式中整个 F_hk 乘以 √δ√(R_iR_j/(R_iR_j))完全一致也解释了为何旧版本需要重新标定四个系数。与新一代 pair_style granular 的关系文档明确指出pair_style granular见 doc/src/pair_granular.rst在法向、切向、滚动、扭转等子模型上大大扩展了pair gran/*系列的能力。现代统一入口 pair_granular.cpp 通过pair_coeff按类型对组合选择子模型normal/tangential/damping/rolling/twisting/heat实现同一套体系内不同材料类型使用不同接触模型。例如 examples/granular/in.chute 使用pair_style granular pair_coeff * * hertz 1 0.5 tangential linear_history 1.0 0.5 0.1而旧式gran/hooke等风格则在需要全局统一参数、快速上手的场景下依然简洁有效例如 examples/granular/in.restitution 中评估颗粒碰撞恢复系数的实验即使用pair_style granularhertz/material子模型。从源码结构看pair_gran_hooke、pair_gran_hooke_history、pair_gran_hertz_history属于相对早期的经典实现其计算内核法向弹簧-阻尼、切向历史剪切、库仑截断与 granular 的默认子模型逻辑一脉相承。加速变体三种风格均提供加速变体需要相应构建支持gran/hookegran/hooke/kkKOKKOS、gran/hooke/ompOpenMPgran/hooke/historygran/hooke/history/omp、gran/hooke/history/kkgran/hertz/historygran/hertz/history/kk、gran/hertz/history/omp。mixing / shift / table / tail、restart 与 rRESPA 行为pair_modify的 mix、shift、table、tail 选项对 gran 系列不适用这些 pair 风格会把自身信息写入二进制 restart 文件因此读取 restart 文件的输入脚本无需重复指定pair_style命令在run_style respa中gran 系列只能通过pair关键字使用不支持inner、middle、outer关键字这些风格的single()函数返回 0.0 作为成对相互作用能量——因为这类耗散势不守恒能量它只返回成对作用力的法向分量。single() 返回的 10 个额外量single()还会计算 10 个额外的成对量可通过compute pair/local命令以p1、p2、…、p10访问索引内容p1–p3作用于颗粒 I 的切向力三个分量p4该切向力的幅值p5–p7法向相对速度的三个分量沿两球心连线方向p8–p10切向相对速度的三个分量对应源码见 pair_gran_hooke.cppsvector[0..2]为切向力分量svector[3]为其模svector[4..6]为法向相对速度svector[7..9]为切向相对速度。构建与运行限制所有 gran 风格都属于GRANULAR 包必须使用该包构建 LAMMPS 才能启用见 Build package 页面包源码位于 src/GRANULAR这些风格要求原子存储力矩与角速度omega由 atom_style 定义并要求存储每个颗粒的半径——atom_style sphere同时满足这两点必须使用comm_modify vel yes命令使 ghost 原子也保存速度信息使用read_restart时不能精确复现运行结果统计上相似即可因为力依赖原子速度但单个接触的累积历史量会写入 restart 文件历史量不会写入 data 文件因此通过read_data重新加载体系时由于丢失接触历史力可能有显著差异。测试与示例验证仓库提供了针对三种风格的回归测试用于校验力与应力计算的确定性unittest/force-styles/tests/atomic-pair-gran_hooke.yamlpair_style gran/hooke 1000.0 350.0 50.0 30.0 0.5 1 limit_dampingunittest/force-styles/tests/atomic-pair-gran_hooke_history.yamlunittest/force-styles/tests/atomic-pair-gran_hertz_history.yaml测试前置命令见 atomic-pair-gran_hooke.yaml给出了运行 gran 风格的最小骨架units lj、atom_style sphere、comm_modify vel yes、atom_modify map array以及set group all diameter、set group all mass颗粒质量与半径是必需的。examples/granular 目录下的in.restitution双球对撞测恢复系数与in.chute二维溜槽流动可作为直接运行的实战起点。相关命令与默认值相关命令pair_coeff。默认值none即所有系数都必须显式给定Kt、gamma_t用 NULL 时分别取2/7 Kn、1/2 gamma_n的默认值。参考文献Brilliantov, Spahn, Hertzsch, Poschel, Phys Rev E, 53, p 5382-5392 (1996)Silbert, Ertas, Grest, Halsey, Levine, Plimpton, Phys Rev E, 64, p 051302 (2001)Zhang and Makse, Phys Rev E, 72, p 011301 (2005)赞分享科研科学计算高性能计算【免费下载链接】lammpsPublic development project of the LAMMPS MD software package项目地址https://gitcode.com/gh_mirrors/la/lammps点击查看免费下载相关推荐LAMMPS fix wall/gran/region 命令详解用任意几何 region 构造颗粒摩擦壁LAMMPS fix wall/gran/region 命令详解用任意几何 region 构造颗粒摩擦壁 导读 fix wall/gran/region 是科研科学计算高性能计算LAMMPS pair_style buck 系列Buckingham 势的语法、参数与源码实现解析LAMMPS pair_style buck 系列Buckingham 势的语法、参数与源码实现解析 导读 本文基于 LAMMPS 官方文档 doc/src/科研科学计算高性能计算LAMMPS pair_style colloid 胶体相互作用势语法、参数与源码实现详解LAMMPS pair_style colloid 胶体相互作用势语法、参数与源码实现详解 导读 本文围绕 LAMMPS 的 pair_style collo科研科学计算高性能计算上一篇WinUtil告别Windows繁琐设置一键搞定系统优化与软件安装下一篇immich-go 更优资产选择算法解析文件去重、日期比对与画质择优的上传判定机制创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考