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

导读

fix wall/gran/region是 LAMMPS GRANULAR 包提供的命令，它把region命令定义的任意几何区域表面当作有摩擦的刚性壁，与有限尺寸的球形颗粒发生接触力学相互作用。与只能描述平面/圆柱的fix wall/gran不同，该命令支持 block、sphere、cylinder、cone 等基本体以及 union/intersect 组合体，还能跟随 region 的move/rotate/变量参数动态运动，因此常被用来搭建漏斗、搅拌器、料仓等复杂颗粒容器。读完本文，你将掌握该命令的完整语法、力模型参数（hooke/hooke/history/hertz/history/granular）、接触几何与有效半径的计算方式、contacts逐原子输出格式，以及运动 region 与 restart 协同使用的关键注意事项。

颗粒在漏斗形 region 壁面内流动的模拟截图 颗粒在搅拌器式 intersect region 壁面内混合的模拟截图

一、命令语法与参数

fix ID group-ID wall/gran/region fstyle fstyle_params wallstyle regionID keyword values ...

• ID、group-ID：与所有 fix 命令一致；
• wall/gran/region：本 fix 的 style 名；
• fstyle：颗粒与壁之间的力相互作用风格；
• fstyle_params：与该力风格对应的参数；
• wallstyle：此处固定为region（其它 wallstyle 选项见 fix wall/gran）；
• region-ID：充当壁面的区域 ID；
• keyword：contacts或temperature（见后文）。

1.1 经典力风格参数

对于hooke、hooke/history、hertz/history，fstyle_params依次为 6 个参数：

单位说明：当fstyle为hooke系列时，Kn、Kt采用力/距离单位；hertz/history因 Hertz 接触公式与压入深度相关，Kn、Kt采用压力单位，gamma_n、gamma_t对应为 1/时间-距离单位。具体讨论见 pair_style gran/*。

其中Kt、gamma_t都可以写NULL：Kt缺省为2/7 * Kn，gamma_t缺省为1/2 * gamma_n。另有可选关键字limit_damping可追加在参数之后。

1.2 granular 风格

当fstyle为granular时，fstyle_params采用与 pair_style granular 的pair_coeff完全相同的语法，可自由组合法向（hooke/hertz/jkr/dmt 等）、切向（linear_nohistory/mindlin 等）、滚动、扭转、阻尼、热传导等子模型。

1.3 可选关键字

• contacts：无值，为每个颗粒生成接触信息（逐原子数组输出，见第五节）；
• temperature：值为温度（数值或v_name形式 equal-style 变量），指定壁温，需配合 granular 热传导子模型使用（见第四节）。

二、官方示例逐条解读

fix wall all wall/gran/region hooke/history 1000.0 200.0 200.0 100.0 0.5 1 region myCone fix 3 all wall/gran/region granular hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 damping velocity region myBox fix 4 all wall/gran/region granular jkr 1e5 1500.0 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall region myCone fix 5 all wall/gran/region granular dmt 1e5 0.2 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall damping tsuji region myCone fix wall all wall/gran/region hooke/history 1000.0 200.0 200.0 100.0 0.5 1 region myCone contacts

• 示例 1：经典hooke/history壁，Kn=1000、Kt=200、gamma_n=200、gamma_t=100、xmu=0.5、dampflag=1；
• 示例 2：granular风格，法向hooke（Kn=1000，恢复系数阻尼参数 50.0），切向linear_nohistory（1.0、0.4），阻尼子模型velocity；
• 示例 3/4：granular风格下分别使用 jkr/dmt 内聚法向模型，叠加 mindlin 切向、sds 滚动、marshall 扭转，示例 4 还显式指定damping tsuji；
• 示例 5：在示例 1 基础上追加contacts关键字输出逐原子接触信息。

可以看出，同一命令既可写"一行 6 参数"的经典风格，也可写出与pair_coeff等价的完整子模型组合——颗粒与壁的力模型可以不同于颗粒间的 pair 模型。

三、工作原理：接触几何、法向与有效半径

3.1 距离与力的方向

颗粒与区域边界的距离定义为到区域表面上最近点的距离；壁对颗粒的力沿该点与颗粒中心的连线方向，也就是该点的表面法向。若区域表面由多个"面"构成，则每个面在足够近时都可对颗粒施加力。例如 region_style block 内部、靠近角落的颗粒可能同时感受到 block 的 1 个、2 个甚至 3 个面的壁力——这正是源码中"单颗粒可同时与多面接触"（ncontact、tmax）设计的物理来源。

3.2 力方程与有效半径

壁/颗粒间作用力的方程与 pair_style gran/*、pair_style granular 文档页上对应的颗粒间方程相同，区别在于有效半径的取法：

• 一般曲面：delta = radius - r为颗粒与壁的重叠量，m_eff为颗粒质量，有效半径用Ri*Rj/(Ri+Rj)计算，其中Rj用壁在接触点的曲率半径替代；
• 凹面（如圆柱内壁）：壁曲率半径为负，会减小有效半径；
• 平面：有效半径就是颗粒自身的半径。

在源码 src/GRANULAR/fix_wall_gran_region.cpp 中可以看到，每个接触会把region->contact[ic]里预先算好的delx/dely/delz、壁曲率半径radj、中心距r直接写入GranularModel对象，随后由model->check_contact()与model->calculate_forces()完成接触判定与受力计算——也就是说，壁接触与颗粒对接触共用同一套granular_model内核。

3.3 运动 region 的壁速度

若 region 是动态的（move/rotate关键字或变量驱动的形状变化），壁速度由region->velocity_contact(vwall, x[i], ic)计算并赋给model->vj（源码 post_force），从而切向摩擦与阻尼力会正确反映壁面本身的运动——这是搅拌器、旋转圆筒等场景的关键。

四、区域定义与使用要点

4.1 内部/外部与 open 选项

region 体积可以是边界面内部或外部，这决定了表面法向（即壁力作用方向）。例外情况是：若 region 命令指定了open选项（一个或多个），颗粒会同时与区域的内部表面和外部表面相互作用。漏斗示例in.granregion.funnel中，锥形料斗就同时用了open 1 open 2打开上下底面：

region conereg cone z ${xc} ${yc} ${radconelo} ${radconehi} ${zconelo} ${zconehi} side in units box & open 1 open 2 #Bottom and top are open

4.2 基本体与组合体

region 可以是 block、sphere、cylinder、cone 等基本体，也可以通过union/intersect组合构造复杂容器。搅拌器示例in.granregion.mixer就是用三个子区域intersect构造出带旋转叶片的料筒：

region cylreg cylinder z 10 10 10 0 20 side in units box variable theta equal (step/400000)*2*PI region b1 block 2 18 9 11 0 4 side out rotate v_theta 10 10 0 0 0 1 units box region b2 block 9 11 2 18 0 3.99999 side out rotate v_theta 10 10 0 0 0 1 units box region mixer intersect 3 cylreg b1 b2 side in fix mixwall all wall/gran/region hertz/history ${kn} ${kt} ${gamma_n} ${gamma_t} ${coeffFric} 1 region mixer

注意b2的高度故意写成3.99999而非4，正是为了规避下文"共面重叠"问题。

4.3 动态 region

region 可通过region命令的move、rotate关键字或变量输入改变形状。若被本 fix 使用，区域表面会按相应方式随时间运动。in.granregion.box示例演示了在运行中途反复region container delete后以不同角速度/转轴rotate v_theta ...重建容器，从而让颗粒容器"翻转"。

4.4 周期性边界注意

区域在 LAMMPS 中不会跨周期性边界包裹（详见 region 文档）。当用 region 作为约束颗粒运动的壁时，必须自行通过region与boundary命令把区域相对周期/非周期边界的位置设置正确。

4.5 尖锐角/边与共面重叠的严格限制

文档明确给出两条必须遵守的警告（源码层面同样会体现）：

1. 凸角不精确：基本体自带尖角/棱边（如 block、cylinder）对内部和外部 region 的壁力计算都是精确的；但union/intersect产生的交线/交角中，凹角（如立方体内角）计算正确，凸角（外角，顶点伸入颗粒所在的体积）LAMMPS 不直接定位，法向在该点不连续，会导致颗粒"看不见壁"→移动 ε 后突然受大力，最坏情况会把颗粒弹出模拟盒。因此不要用带凸角/凸边的 union/intersect region 配本命令（单个子区域自带的凸角/凸边仍可接受）。
2. 禁止共面重叠：union/intersect的子区域不得共享属于整体外边界的重合面（内边界可以），因为 LAMMPS 在计算壁/颗粒相互作用时会丢弃属于多个子区域的点以避免重复计数，共面可能导致该面"隐形"。解决办法：把两个面错开 ε 位置（如 mixer 示例的3.99999）。

五、contacts逐原子输出

指定contacts后，本 fix 为每个本地颗粒（每处理器上的 nlocal）生成至少 8 列的逐原子数组，无接触时全部为 0。前 8 列定义如下：

如果某个 granular 子模型还会计算额外接触量（例如热子模型计算换热量），这些量会追加到数组末尾，顺序为：法向子模型 → 阻尼 → 切向 → 滚动 → 扭转 → 热。各子模型额外量的含义见 pair_style granular。

从源码看，接触点坐标正是x[i] - dx（fix_wall_gran_region.cpp），即颗粒中心减去接触位移矢量；model->svector携带各子模型额外量。测试 unittest/granular/tests/dem11-regionwall-3d-si.yaml 使用单颗粒撞击圆柱壁（fix wreg all wall/gran/region granular hertz ${kn} ${en} tangential linear_nohistory 0.0 0.0 damping coeff_restitution region rcyl）并与wall_restitution解析模型对比位置/速度，验证了壁接触力学的正确性。

六、temperature关键字与热传导

temperature值可以是数值或 equal-style 变量（写成v_name）。变量形式每个时间步都会被求值，其值作为壁温。该选项必须与以下三件套配合：

1. pair_style granular 中定义热传导子模型（heat）；
2. fix property/atom 存储颗粒温度和热流；
3. fix heat/flow 积分热流。

源码中heat_flag置位后，壁温写入model->Tj = Twall，颗粒温度写入model->Ti，计算得到的热流model->dq累加到heatflow[i]（fix_wall_gran_region.cpp）。

七、restart、fix_modify 与运行限制

• restart 支持：与fix wall/gran类似，本 fix 会把与壁相互作用的原子的剪切摩擦历史写入二进制 restart 文件，保证使用带 history 效应的颗粒势时续跑正确；此外还会在 restart 文件中保存运动 region 的信息。
• region 定义不在 restart 中：region 本身不会被写入 restart（见 read_restart），因此续跑时必须重新定义 region；若它是运动 region，其运动属性要与写 restart 的模拟一致。特别是想改变 region 运动属性（如速度）时，要保证初始 restart 时间步上 region 的位置/朝向与写文件时一致。做不到时，可以在 restart 脚本中用不同的 fix ID新定义一个fix wall/gran/region来忽略旧信息；若只想保留剪切历史而丢弃 region 运动信息，则用相同 fix ID但换一个不同的 region ID。
• 源码中FixWallGranRegion::init()（fix_wall_gran_region.cpp）会在两次 run 之间检查 region 属性是否变化、restart 信息是否不一致，必要时发出警告并region->reset_vel()重置运动。
• fix_modify：无任何选项对本 fix 相关。
• start/stop：本 fix 无参数可与run命令的 start/stop 关键字连用；能量最小化（minimize）期间本 fix 不被调用。

八、dump image 可视化与相关命令

• 本 fix不支持dump image 的fix关键字；要可视化区域壁面，请用dump image的region关键字直接绘制 region 几何。
• 相关命令：fix_move、fix wall/gran、fix wall/region、pair_style granular、region。
• 默认值：无。

九、源码实现纵深：多接触剪切历史管理

与平面壁（单接触）不同，region 壁允许一个颗粒同时接触多个面，因此剪切历史必须按"接触"维度存储。实现类FixWallGranRegion（继承自FixWallGran）在 src/GRANULAR/fix_wall_gran_region.h 中声明了如下数据结构：

• tmax：单个颗粒最多可接触的壁面数（来自region->tmax）；
• ncontact[i]：颗粒 i 当前的剪切接触数；
• walls[i][*]：每个接触对应哪个壁面；
• history_many[i][*][*]：每个颗粒、每个接触的剪切历史向量；
• c2r[]：接触到 region 接触列表的映射。

核心流程在post_force（fix_wall_gran_region.cpp）：

1. 先由region->match()判断颗粒是否在区域内，再由region->surface()枚举颗粒周围（膨胀pulloff_distance后的截断距离内）的接触面并填充region->contact[]；
2. use_history时，单接触直接维护、多接触调用update_contacts()：遍历旧接触，若不在新接触列表中则删除并回收剪切历史；遍历新接触，若不在压缩后的旧列表中则追加零历史（update_contacts）；
3. 对每个接触依次填充model->dx/radi/radj/r/meff/vi/omegai/history/Ti，调用model->calculate_forces()得到力与力矩，累加进f[i]、torque[i]，热模型下累加heatflow[i]；
4. contacts选项置位时把接触标志、三分力、接触点、半径及子模型额外量写入逐原子数组。

内存估算（memory_usage，fix_wall_gran_region.cpp）约等于nmax*(4 + 4*tmax + 8*tmax*size_history)字节量级，tmax与size_history越大，多接触历史开销越高；交换（pack_exchange/unpack_exchange）与 restart（pack_restart/unpack_restart）均按实际ncontact[i]动态打包，避免为每个颗粒都预留tmax空间。此外，接触类型通过GranularModel::contact_type = WALLREGION（granular_model.h 中PAIR/WALL/WALLREGION/SURFACE枚举）告知力内核，使其按壁接触的几何约定（曲率半径替代 Rj、radj=0时退化为平面壁）计算。

十、完整实战：三个 granregion 输入脚本

示例目录 examples/granregion 包含三个可直接运行的脚本：

1.in.granregion.box（最简入门）：单位 lj，atom_style sphere，边界 fff，用 block 容器壁：

region container block -6 6 -6 6 -6 6 units box fix container all wall/gran/region hooke/history 4000.0 NULL 100.0 NULL 0.5 1 region container

配合fix pour向盒内注球，随后多次region container delete后以不同转轴rotate v_theta重建，演示容器旋转/翻转时壁面的动态响应。注意 Kt、gamma_t 均用NULL，即采用Kt=2/7*Kn、gamma_t=1/2*gamma_n的默认值。

2.in.granregion.funnel（漏斗流）：单位 si，region conereg cone ... open 1 open 2与cylreg（先开顶、后开底）union 成料斗hopreg，fix hopper3 all wall/gran/region hertz/history ${kn} ${kt} ${gamma_n} ${gamma_t} ${coeffFric} 1 region hopreg。脚本第一阶段注球静置，第二阶段unfix后重建开口向下的cylreg（open 1 open 2）并重新定义 union，实现"拔掉塞子、颗粒流出"。脚本头部还给出了标准的 Hertz 参数标定公式：

kn = 4*G/(3*(1-nu))，kt = 4*G/(2-nu) gamma_n = sqrt( a*2*kn/min_mass/(1+0.25*a) )，其中 a = (-2*ln(恢复系数)/PI)^2 gamma_t = 0.5*gamma_n dt = 0.05 * PI / sqrt( 2*kn/min_mass - gamma_n/4 )

3.in.granregion.mixer（搅拌器）：如上文 4.2 节所示，用intersect 3 cylreg b1 b2构造带旋转叶片的料筒，叶片转速由变量theta = (step/400000)*2*PI驱动的rotate决定；该脚本也演示了如何规避共面重叠问题（b2 高度 3.99999）。

三者共同点：atom_style sphere、comm_modify vel yes、neighbor skin multi、平衡负载的fix balance，以及fix gravity+fix nve/sphere驱动颗粒运动。

结语

fix wall/gran/region把 LAMMPS 的几何 region 系统与颗粒接触力学无缝对接：底层复用与pair_style granular相同的GranularModel内核，上层通过多接触剪切历史管理、动态 region 速度传递与 restart 兼容，支撑漏斗、搅拌器、料仓等真实工业场景的离散元建模。使用时请务必遵守两条红线——避免 union/intersect 的凸角区域、避免子区域共面重叠，并妥善处理周期性边界与 restart 中的 region 重定义。若仅需平面/圆柱等简单壁面，可退回更轻量的 fix wall/gran。

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