LAMMPS `fix wall/gran/region` 命令详解：用任意几何 region 构造颗粒摩擦壁

科研科学计算高性能计算【免费下载链接】lammpsPublic development project of the LAMMPS MD software package项目地址https://gitcode.com/gh_mirrors/la/lammps点击查看免费下载导读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/granregion-ID充当壁面的区域 IDkeywordcontacts或temperature见后文。1.1 经典力风格参数对于hooke、hooke/history、hertz/historyfstyle_params依次为 6 个参数参数含义单位Kn法向排斥弹性常数力/距离 或 压力见下Kt切向接触弹性常数力/距离 或 压力gamma_n法向碰撞阻尼系数1/时间 或 1/时间-距离gamma_t切向碰撞阻尼系数1/时间 或 1/时间-距离xmu静态屈服准则无量纲0.0 ~ 1.0e4dampflag0 或 1是否包含切向阻尼力—单位说明当fstyle为hooke系列时Kn、Kt采用力/距离单位hertz/history因 Hertz 接触公式与压入深度相关Kn、Kt采用压力单位gamma_n、gamma_t对应为 1/时间-距离单位。具体讨论见 pair_style gran/*。其中Kt、gamma_t都可以写NULLKt缺省为2/7 * Kngamma_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壁Kn1000、Kt200、gamma_n200、gamma_t100、xmu0.5、dampflag1示例 2granular风格法向hookeKn1000恢复系数阻尼参数 50.0切向linear_nohistory1.0、0.4阻尼子模型velocity示例 3/4granular风格下分别使用 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/(RiRj)计算其中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 open4.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 动态 regionregion 可通过region命令的move、rotate关键字或变量输入改变形状。若被本 fix 使用区域表面会按相应方式随时间运动。in.granregion.box示例演示了在运行中途反复region container delete后以不同角速度/转轴rotate v_theta ...重建容器从而让颗粒容器翻转。4.4 周期性边界注意区域在 LAMMPS 中不会跨周期性边界包裹详见 region 文档。当用 region 作为约束颗粒运动的壁时必须自行通过region与boundary命令把区域相对周期/非周期边界的位置设置正确。4.5 尖锐角/边与共面重叠的严格限制文档明确给出两条必须遵守的警告源码层面同样会体现凸角不精确基本体自带尖角/棱边如 block、cylinder对内部和外部 region 的壁力计算都是精确的但union/intersect产生的交线/交角中凹角如立方体内角计算正确凸角外角顶点伸入颗粒所在的体积LAMMPS 不直接定位法向在该点不连续会导致颗粒看不见壁→移动 ε 后突然受大力最坏情况会把颗粒弹出模拟盒。因此不要用带凸角/凸边的 union/intersect region 配本命令单个子区域自带的凸角/凸边仍可接受。禁止共面重叠union/intersect的子区域不得共享属于整体外边界的重合面内边界可以因为 LAMMPS 在计算壁/颗粒相互作用时会丢弃属于多个子区域的点以避免重复计数共面可能导致该面隐形。解决办法把两个面错开 ε 位置如 mixer 示例的3.99999。五、contacts逐原子输出指定contacts后本 fix 为每个本地颗粒每处理器上的 nlocal生成至少 8 列的逐原子数组无接触时全部为 0。前 8 列定义如下列号值单位11.0颗粒与壁接触否则 0.0—2壁施加的力 f_x力单位3壁施加的力 f_y力单位4壁施加的力 f_z力单位5壁上接触点的 x 坐标距离单位6壁上接触点的 y 坐标距离单位7壁上接触点的 z 坐标距离单位8原子半径 r距离单位如果某个 granular 子模型还会计算额外接触量例如热子模型计算换热量这些量会追加到数组末尾顺序为法向子模型 → 阻尼 → 切向 → 滚动 → 扭转 → 热。各子模型额外量的含义见 pair_style granular。从源码看接触点坐标正是x[i] - dxfix_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。变量形式每个时间步都会被求值其值作为壁温。该选项必须与以下三件套配合pair_style granular 中定义热传导子模型heatfix property/atom 存储颗粒温度和热流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-tmaxncontact[i]颗粒 i 当前的剪切接触数walls[i][*]每个接触对应哪个壁面history_many[i][*][*]每个颗粒、每个接触的剪切历史向量c2r[]接触到 region 接触列表的映射。核心流程在post_forcefix_wall_gran_region.cpp先由region-match()判断颗粒是否在区域内再由region-surface()枚举颗粒周围膨胀pulloff_distance后的截断距离内的接触面并填充region-contact[]use_history时单接触直接维护、多接触调用update_contacts()遍历旧接触若不在新接触列表中则删除并回收剪切历史遍历新接触若不在压缩后的旧列表中则追加零历史update_contacts对每个接触依次填充model-dx/radi/radj/r/meff/vi/omegai/history/Ti调用model-calculate_forces()得到力与力矩累加进f[i]、torque[i]热模型下累加heatflow[i]contacts选项置位时把接触标志、三分力、接触点、半径及子模型额外量写入逐原子数组。内存估算memory_usagefix_wall_gran_region.cpp约等于nmax*(4 4*tmax 8*tmax*size_history)字节量级tmax与size_history越大多接触历史开销越高交换pack_exchange/unpack_exchange与 restartpack_restart/unpack_restart均按实际ncontact[i]动态打包避免为每个颗粒都预留tmax空间。此外接触类型通过GranularModel::contact_type WALLREGIONgranular_model.h 中PAIR/WALL/WALLREGION/SURFACE枚举告知力内核使其按壁接触的几何约定曲率半径替代 Rj、radj0时退化为平面壁计算。十、完整实战三个 granregion 输入脚本示例目录 examples/granregion 包含三个可直接运行的脚本1.in.granregion.box最简入门单位 ljatom_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即采用Kt2/7*Kn、gamma_t1/2*gamma_n的默认值。2.in.granregion.funnel漏斗流单位 siregion conereg cone ... open 1 open 2与cylreg先开顶、后开底union 成料斗hopregfix hopper3 all wall/gran/region hertz/history ${kn} ${kt} ${gamma_n} ${gamma_t} ${coeffFric} 1 region hopreg。脚本第一阶段注球静置第二阶段unfix后重建开口向下的cylregopen 1 open 2并重新定义 union实现拔掉塞子、颗粒流出。脚本头部还给出了标准的 Hertz 参数标定公式kn 4*G/(3*(1-nu))kt 4*G/(2-nu) gamma_n sqrt( a*2*kn/min_mass/(10.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 gravityfix nve/sphere驱动颗粒运动。结语fix wall/gran/region把 LAMMPS 的几何 region 系统与颗粒接触力学无缝对接底层复用与pair_style granular相同的GranularModel内核上层通过多接触剪切历史管理、动态 region 速度传递与 restart 兼容支撑漏斗、搅拌器、料仓等真实工业场景的离散元建模。使用时请务必遵守两条红线——避免 union/intersect 的凸角区域、避免子区域共面重叠并妥善处理周期性边界与 restart 中的 region 重定义。若仅需平面/圆柱等简单壁面可退回更轻量的 fix wall/gran。赞分享科研科学计算高性能计算【免费下载链接】lammpsPublic development project of the LAMMPS MD software package项目地址https://gitcode.com/gh_mirrors/la/lammps点击查看免费下载相关推荐LAMMPS fix wall/ees 与 fix wall/region/ees 完全指南椭圆体粒子与 EES 势墙壁的相互作用LAMMPS fix wall/ees 与 fix wall/region/ees 完全指南椭圆体粒子与 EES 势墙壁的相互作用 导读 fix wall/e科研科学计算高性能计算LAMMPS 壁面Walls建模完全指南从粗糙颗粒壁到理想光滑墙与摩擦墙LAMMPS 壁面Walls建模完全指南从粗糙颗粒壁到理想光滑墙与摩擦墙 本文是 LAMMPS 分子动力学模拟中壁面建模的实战指南围绕官方 Howto科研科学计算高性能计算LAMMPS compute temp/region 命令详解几何区域温度计算与分区恒温器实战指南LAMMPS compute temp/region 命令详解几何区域温度计算与分区恒温器实战指南 导读 compute temp/region 是 LAMM科研科学计算高性能计算上一篇DeepSeek-LLM推理框架对比vLLM vs HuggingFace Transformers下一篇如何快速掌握OrcaSlicer从基础设置到高级功能的完整指南创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考