科研科学计算高性能计算【免费下载链接】lammpsPublic development project of the LAMMPS MD software package项目地址https://gitcode.com/gh_mirrors/la/lammps点击查看免费下载本文围绕 LAMMPS 中由三角形3D或线段2D构成的颗粒表面granular surfaces展开系统讲解其两种存储形式全局global与局部local、输入文件格式molecule、STL、data 文件、表面属性与原子风格要求、拓扑与连通性规则、表面运动控制以及底层基于连通性网络的粒子-表面受力计算模型。读完本文你将能够为颗粒sphere模拟构造任意复杂形状的边界几何如螺旋输送器、搅拌桨、容器壁理解何时该用fix surface/global或fix surface/local并掌握如何在 2D/3D 模型中正确建立可移动、可连通的颗粒表面。引言什么是颗粒表面正如 Howto granular 文档页所解释的颗粒系统由具有显式直径的有限尺寸球形或非球形粒子组成区别于点粒子因此它们拥有角速度可以被施加力矩而发生转动。颗粒表面则是一组三角形3D 系统或线段2D 系统的集合它们作为边界与颗粒相互作用特别适合用来定义复杂的壁面或边界几何。粒子与表面之间的相互作用可以使用与粒子-粒子相互作用类似的选项进行指定。下图展示了两个典型应用实例均使用 dump image 命令渲染螺旋输送器screw feeder一个圆柱形容器从上方接收颗粒流颗粒被旋转的螺旋桨向前输送见 examples/gransurf/in.screwfeeder.local糸川小行星Itokawa asteroid形状容器依据真实小行星外形构建的表面作为容器被多分散颗粒填充。在 examples/gransurf 目录中包含了上述及其他边界几何的多个示例脚本。颗粒表面功能由GRANSURF包实现核心源码位于 src/GRANSURF主要文件包括fix_surface_global.cpp、fix_surface_local.cpp、pair_surf_granular.cpp与surf_extra.cpp。全局表面global与局部表面local理解 LAMMPS 颗粒表面的关键是知道它支持两种形式且同一模拟中不能同时使用两者。全局表面global每个处理器都存储全部三角形/线段的副本。适用于三角形/线段数量适中modest number的场景。它们可以是很大的三角形/线段任意一个都可以在某个或多个维度上跨越模拟盒尺寸的很大一部分。对应的启用命令是 fix surface/globalfix ID group-ID surface/global input args input args ... model args model args ... keyword value ...全局表面通过两种方式读入molecule 文件由 molecule 命令预先定义。文件应使用头关键字header keywords定义三角形或线段并包含Triangles或Lines节sectionSTLstereolithography文件文本或二进制格式均可命令自动探测格式定义一组三角形只能用于 3D 模拟。局部表面local三角形/线段的集合以与粒子相同的方式分布在处理器之间。每个处理器被分配模拟盒的一个子域拥有其中心点落在该子域内的三角形/线段与粒子类似每个处理器还可以拥有有限尺寸与该子域重叠的三角形/线段的 ghost 副本。此时系统中的三角形/线段总数可以非常大。对应的启用命令是 fix surface/localfix ID group-ID surface/local input args input args ... keyword value ...局部表面有三种定义方式molecule 文件与全局表面相同三角形/线段列表按中心点分布到各处理器STL 文件同上LAMMPS data 文件通过 read_data 命令读入。若文件中包含Triangles或Lines节则三角形/线段与 data 文件中的粒子一起被读入并分布前提是使用了合适的 atom_style详见下文。从 fix_surface_local.cpp 的源码可以看到input关键字支持molmolecule 模板与stl两种来源多个input可以重复出现其三角形/线段的并集构成局部表面集合data 文件方式则由atom_style tri/line在读取 data 时直接创建三角形/线段粒子。两种形式的选型准则维度全局表面global局部表面local存储方式每个处理器存全部副本按子域分布在处理器间含 ghost 副本适用规模适中数量可非常大尺寸要求可很大甚至跨越模拟盒应小不超过几个粒子直径速度影响—若存在哪怕一个大的三角形/线段邻表截断与通信截断会被相应放大拖慢模拟启用命令fix surface/globalfix surface/local对于局部表面如果确实需要大三角形/线段可以将其拆分为多个相连的小三角形/线段而不改变整体表面的拓扑此时 neighbor 的multi选项可能有用。关于复用 3D 打印 STL 文件的重要提醒互联网上有大量为 3D 打印生成的 STL 文件例如上图水泥搅拌车玩具模型中用于打印搅拌筒的 STL。原则上它们可用于 LAMMPS但有几点需要注意退化三角形degenerate triangles对 3D 打印无碍但对 LAMMPS 是问题。LAMMPS 会尝试忽略这些三角形或报错终止。实际上从 fix_surface_global.rst 可知自 2Sep2026 版本起读入 STL 时会跳过具有重复或共线角点零面积的退化三角形并打印被跳过的数量警告重复三角形duplicate triangles会使 LAMMPS 中力翻倍这通常不是你想要的细长三角形long and narrow triangles会迫使 LAMMPS 使用很大的邻表与很长的通信截断两者都会被检测并打印警告或错误。因此多数情况下建议重新网格化remeshSTL 文件用更规则、更少的三角形重建表面。一种做法是使用 MeshLab/PyMeshLabimport pymeshlab ms pymeshlab.MeshSet() ms.load_new_mesh(input.stl) ms.meshing_isotropic_explicit_remeshing(iterations5, targetlenpymeshlab.PercentageValue(1.0) # target edge length as % of bbox diagonal ) ms.save_current_mesh(output.stl)表面属性type 与 molecule ID对全局与局部表面每个三角形/线段都被赋予一个type和一个molecule ID这是在从 molecule、STL 或 data 文件读入时完成的。type用于定义每个三角形/线段的颗粒相互作用风格。同一个由相连三角形/线段构成的物体内可以混用不同的 type因而不同的风格。颗粒-表面接触模型通过 fix surface/global 的model关键字指定语法为model ptype stype fstyle fstyle-args其中ptype/stype支持通配符*、*n、n*、m*nfstyle可选hooke、hooke/history、hertz/history或通用granular。力的计算公式与对应 pair_style granular 页面相同只是取其中一个粒子趋于无穷大半径与质量即平坦表面的极限delta radius - r为粒子与三角形/线段的几何重叠m_eff为粒子质量接触有效半径RiRj/(RiRj)取为粒子半径。Kn、Kt、gamma_n、gamma_t、xmu、dampflag与可选limit_damping参数含义与单位均与 pair_style gran/* 一致Kt可用 NULL 取默认2/7*Kngamma_t可用 NULL 取默认1/2*gamma_nmolecule ID用于限制哪些三角形/线段被视为相连connected因此意图是对每一组相互连接的三角形/线段赋予唯一的值仿佛每个物体是一个分子。详见下文表面连通性一节。由于 STL 文件不包含 type 与 molecule IDfix surface/global 与 fix surface/local 命令通过input stl stype smol stlfile指定赋予读入三角形的 type 与 molecule ID。全局表面还支持用fix_modify type/region命令按几何区域赋予 type用fix_modify mol/region命令按区域赋予 molecule ID。对于局部表面由于它们实际上就是粒子可以用 set 命令修改任一三角形/线段的 type 或 molecule ID。对于局部表面molecule ID 还可用于定义 group例如group case molecule 1从而为不同表面物体分配不同运动其他操作 group 或 molecule 的 LAMMPS 命令也可用于统计或输出单个表面物体的信息。在 fix_surface_global.cpp 源码中可以看到表面 type 与 molecule ID 的数量上限MAXSURFTYPE 1024与MAXSURFMOL 1024命令还提供了smaxtype与smaxmol关键字来显式扩大上限这在计划事后用fix_modify type/region、fix_modify mol/region赋更大的值时有用。颗粒表面所需的原子风格atom styles局部表面的三种定义方式中三角形/线段在 LAMMPS 内部都以三角形风格或线段风格粒子存储因此必须使用包含下列子风格的混合原子风格atom_style tri3D 模拟atom_style line2D 模拟同时atom_style hybrid命令还必须为与表面相互作用的颗粒定义 atom_style sphere 子风格。示例见 examples/gransurf/in.screwfeeder.localatom_style hybrid tri sphere 1 atom_modify map array对于 molecule/STL 文件输入fix surface/local 命令读取文件后为每个表面创建一个新的三角形或线段粒子对于 data 文件输入三角形/线段粒子在读取 data 文件时创建。从 fix_surface_local.cpp 中#include atom_vec_line.h、#include atom_vec_tri.h可以看出该 fix 依赖三角形/线段原子向量实现存储。全局表面的颗粒模拟中不应使用定义了三角形/线段风格的混合原子风格。通常只需 atom_style sphere 定义模拟中的粒子属性三角形/线段由 fix surface/global 命令存储而非作为三角形/线段粒子存在对应源码 fix_surface_global.cpp 的构造函数中检查了radius与omega属性。表面拓扑规则对于全局与局部表面颗粒粒子与每个三角形或线段的两侧都发生相互作用——这意味着搅拌桨叶片之类的表面可以是无限薄的。三角形与线段被视为相连并构成连续表面当且仅当它们 (1) 具有相同的 molecule ID 且 (2) 共享公共边或公共角点三角形或共享端点线段。每个三角形边或角点可以被多个相邻三角形共享若一条边的两个端点同时是两个三角形的角点则该边被这两个三角形共享若一个点同时是两个三角形的角点则该角点被两个三角形共享线段端点同理。不要求每条边/角点/端点必须与另一个三角形/线段相连。没有连接的边或点称为自由未连接边/点粒子与自由边或角点的相互作用方式与粒子重叠于三角形/线段内部所产生的力保持一致。LAMMPS不检查两个三角形或线段是否相交——只要对表面集合的几何有意义允许相交。但若缺失连通性导致力不准确相交可能引发问题。有效相交示例2D 搅拌器考虑一个 2D 颗粒搅拌模拟用全局线段同时定义箱形容器与中心的搅拌器。4 个搅拌桨叶呈大 X 形用fix_modify move命令见下驱动旋转示意图见 doc/src/img/gransurf_mixer.png两个桨叶可定义为在中心相交的 2 条线段左或用共享中心公共端点的 4 条线段中或用 4 条不相接触、各自端点靠近中心与中心距离小于颗粒半径的线段右。无论哪种方式当粒子非常靠近搅拌器中心时它会如预期与两条邻近线段相互作用。无效相交示例2D T 形物体考虑由仅 2 条线段定义的 T 形物体竖直线段蓝终止于水平线段红的中点。缺少正确连通性时当粒子灰水平移动跨越 T 顶部时无法屏蔽censor几何上被隐藏的竖直线段产生的力。同样的 T 形若用 3 条线段红、绿、蓝定义三者共享 T 顶部中心的公共端点连通性将正确屏蔽竖直线段蓝的力。示意见 doc/src/img/gransurf_T_geom.png。表面连通性Surface Connectivity当多个三角形/线段用于定义平坦、平缓弯曲或带有尖锐棱角/拐角的连续表面时LAMMPS 会检测同一 molecule 中两个及以上线段2D共享同一端点或两个及以上三角形3D共享同一边或同一角点。这一连通性在内部存储用于在粒子同时与 2 个及以上相连三角形/线段重叠时计算精确的力。以 2D 搅拌器改用局部线段定义的模型为例每个搅拌桨叶的平坦表面以及容器箱壁面由多个小线段定义。这些线段必须相连这样当粒子同时接触两条相邻线段时其受力与接触一条长线段中部时相同。保证连通性被正确检测的方法molecule 或 STL 文件输入全局与局部表面均可文件为每个三角形/线段列出 3 个角点三角形或 2 个端点线段。LAMMPS只在两个三角形/线段共享的点**完全相同EXACT**时才建立连接角点/端点连接要求两个对象中的同一点一致边连接要求两个三角形中的两个点均一致data 文件输入局部表面同样列出角点/端点但此时 LAMMPS允许共享点的不精确INEXACT匹配来建立连接。INEXACT 意味着两个点可以相差 EPSILON其中 EPSILON 定义为系统中最小三角形/线段的尺寸的极小比例1.0e-4。允许不精确匹配的原因是 data 文件可以由多种方式生成包括 LAMMPS 运行期间通过 write_data 命令自行写出。内部而言三角形/线段粒子并不直接存储角点而是存储三角形/线段的中心点、取向以及从中心到每个角点的位移矢量。因此为两个不同的三角形/线段写出新的角点值时可能因有限精度算术的舍入而略有差异。注意若两个三角形/线段的角点之间距离很小但大于EPSILON它们不会被视为相连若一个三角形/线段的角点非常靠近甚至位于另一个三角形表面或另一个线段中部也不会相连。这类颗粒表面通常是有问题的应尽量避免但 LAMMPS 不检查这些条件连通性仅在相同 molecule ID的两个三角形/线段之间定义——这样两个 molecule 的表面可以独立运动见下节若三角形/线段具有自由边或自由角/端点未与其他三角形/线段连接当球形粒子中心的最邻近接触点位于自由边或就是自由角/端点上时颗粒粒子仍会与该三角形/线段相互作用。表面运动Surface Motion默认情况下无论全局还是局部表面三角形/线段在模拟中静止不动。三角形/线段对颗粒粒子施加力与力矩但粒子施加于三角形/线段的反作用力/反作用力矩不会使其移动。不过三角形/线段可以按预定方式运动例如前述 2D 搅拌桨的旋转。使用以下两个命令全局表面fix_modify move该命令专属于fix surface/global局部表面fix movefix_modify move可以按多种方式平移、旋转等移动指定子集的三角形/线段。哪些三角形移动基于每个三角形的molecule ID指定molecule ID 在fix surface/global定义表面时指定也可由fix_modify mol/region命令指定。语法示例来自 fix_surface_global.rstfix_modify 1 move 2 rotate 0 0 0 0 0 1 25 fix_modify 1 move 1,3*5,8* rotate 0 0 0 0 0 1 25 fix_modify 1 type/region 3 myBlockfix move对局部表面按 group-ID 参数指定哪些三角形移动局部表面的 group 可用 group 命令结合molecule关键字定义。示例examples/gransurf/in.tri.gran.data中三角形粒子被归入tri组后fix 1 tri move rotate 0 0 0 0 0 1 25在 examples/gransurf/in.screwfeeder.local 中螺旋输送器案例即用 group 加fix move驱动旋转fix 1 all surface/local input stl 2 1 surf.container.stl input stl 3 2 surf.helix.stl flat 10.0 group case molecule 1 group corkscrew molecule 2 fix 5 corkscrew move rotate 0 0 0 0 1 0 1000 units box注意对于由两个及以上相连三角形/线段构成的物体为其中一部分而非全部相连三角形/线段分配运动是错误的因为这会破坏连接。LAMMPS 对全局表面会检查这一点对局部表面只检查任意fix move实例的 fix group 是否包含一组相连三角形/线段中的全部或全部不包含。力的计算Calculation of Forces连接的分类与点的分类生成表面连通性后LAMMPS 依据法向量之间的角度差将每个连接分类为平坦flat或非平坦non-flat再根据法向量将两侧分为凹concave或凸convex。对于完全平坦法向量平行的表面凹/凸的指定是任意的。随后每个线段/三角形的点或边被分类为内部internal、外部external或未连接unconnected线段端点仅具有平坦连接 → 内部至少一个非平坦凹或凸连接 → 外部无连接 → 未连接。三角形的边同理三角形角点继承于汇合于该角点的两条边的分类——若任一边未连接则角点未连接否则若任一边外部则角点外部若两边均为内部则角点内部。对应源码见 fix_surface_global.cppenum{INTERNAL 0, EXTERNAL, UNCONNECTED}FLATTHRESH 1.0-cos(MY_PI/180.0)即 1 度与flat关键字默认的最大平坦角一致。接触集合与有效重叠/法向为计算粒子受力LAMMPS 找出所有与粒子几何接触的线段/三角形集合 S_all。对每个线段/三角形i计算几何重叠 δ_i 与从该线段/三角形接触点指向粒子的法向量 n̂_{r,i}。接触模型决定粒子-表面力依赖的变量相对速度、角速度、滑动历史等。其中最关键的是有效重叠距离 δ_f与有效接触点等价地法向力方向 n̂_f粒子仅接触单个线段/三角形i时接触点与重叠由几何唯一确定δ_f δ_in̂_f n̂_{r,i}粒子同时接触两个及以上未连接的线段/三角形时分别施加两个力重叠与方向由各自的几何值决定粒子同时接触两个及以上相连的线段/三角形时δ_f 与 n̂_f 的计算较为复杂是本节剩余部分的主题。一致接触侧面的确定首先LAMMPS 需要为 S_all 中每个线段/三角形确定一组一致的接触侧面。单个线段/三角形时接触哪一侧是明确的但当粒子接触两个及以上相连的线段/三角形时这取决于连通性网络。下图展示 2D 系统粒子灰圆接触菱形物体中两条绿色线——仅凭局部判断会认为粒子接触两条线的内侧/外侧蓝色法向其一指向物体内部而通过凸连接的上下文才能识别物理正确的表面法向红色。LAMMPS 从主接触出发遍历所有接触表面的连接评估一致的侧面集合设定 n̂_s 的符号若仍有未检查表面则取重叠最大的未检查表面并重复该过程示意见 doc/src/img/gransurf_consistent_side.png。复合集合聚类与力合成接下来LAMMPS 将所有接触且相连的线段/三角形聚类为若干互不相同的复合集合composite sets每个集合由相互平坦相连、作为一个物理对象作用的线段/三角形表面构成记为 S_n。例如粒子接触两个平坦相连表面与第三个凹相连表面时会分成大小分别为 2 与 1 的两个复合表面。每个复合集合计算一个有效接触点对粒子施加单一力。聚类过程从主接触出发检查其所有连接平坦连接加入当前集合 S_n 并进入待遍历列表非平坦连接跳过但重叠更小的凸连接首先被标记为隐藏hidden。随后逐个遍历待遍历列表直至为空——这等价于依次搜索主接触的 1 阶平坦邻域、2 阶平坦邻域等。然后取重叠最大的下一个未处理表面作为新的主接触重复以创建新的复合集合 S_m直至 S_all 中所有接触表面都归入某个复合集合可退化为单一表面。遍历中任何隐藏标记都会传递给随后遍历的平坦连接表面使隐藏状态沿凸转角级联传播。对每个平坦表面复合集合 S_nLAMMPS 计算单一力先清零集合中任何隐藏表面的单独重叠δ_i 0有效重叠取最大值δ_f max(δ_i | i ∈ S_n)力的方向为集合内所有表面的加权平均n̂_f A · Σ W_i δ_i n̂_{f,i}其中 n̂_{f,i} 是该表面计算出的力方向W_i 为逐表面权重A 为归一化因子。在计算各表面 n̂_{f,i} 与 W_i 之前LAMMPS 计算一般性的外部/内部接触权重W_ext ≡ (δ_max,ext / δ_max)² 与 W_int ≡ 1 − W_ext其中 δ_max,ext 是复合集合 S_n 中外部接触表面的最大重叠。该加权用于强调粒子沿凸转角移动时凸边界上表面的贡献。2D 中的方向规则在 2D 中表面i的权重 W_i 依据其状态取 W_ext 或 W_int。默认 n̂_{f,i} 就是从该表面局部接触点指向粒子的方向 n̂_{r,i}例如接触点在线段内部时n̂_{r,i} 等价于表面法向 n̂_{s,i}。有两个例外(a) 接触在凹连接点时n̂_{f,i} n̂_{s,i}(b) 接触在凸连接点且 n̂_{r,i} 有指向相邻线段j方向的分量时n̂_{f,i} n̂_{s,j}。这些规则限制合力指向相连线段j内部的程度以保证粒子转过弯角时 n̂_{f,i} 平滑变化。细节注意非平坦连接表面在凸转角后被隐藏δ_i 0故凸场景中的限制只对既平坦又凸的连接相关若两条平坦线段完全平行n̂_{f,i} n̂_{s,i} n̂_{s,j}凹/凸指定无效果若一个点被多于两条线段共享LAMMPS 找出其法向最接近 n̂_{s,i} 的连接线段以判定该连接为凹或凸。下图左面板展示粒子在不同位置红、绿、蓝、紫接触由两条线构成的凹弯珊瑚棕色的场景左线段标记为i右线段标记为j力方向 n̂_{f,i} 由箭头指示整个接触过程中 n̂_{f,i} n̂_{s,i}。右面板展示接触点位于凸角n̂_{f,i} n̂_{r,i}红、绿除非 n̂_{r,i} 有指向相邻线段j的分量此时 n̂_{f,i} n̂_{s,j}蓝、紫。这些规则本质上保证了原子跨越线段运动时力的合理连续性见 doc/src/img/gransurf_nonflat_turn.png。3D 中的方向规则3D 中首先考虑无自由/未连接边接触的场景所有边至少有一个连接。每个复合集合 S_n 计算单一力接触点在三角形内部n̂_{f,i} n̂_{s,i}接触点在边上行为类似 2D 点——凹连接 ⇒ n̂_{f,i} n̂_{s,i}凸连接 ⇒ n̂_{f,i} n̂_{r,i}除非 n̂_{r,i} 指向相邻三角形j此时 n̂_{f,i} n̂_{s,j}。这是通过比较 n̂_{s,i}、n̂_{s,j}、n̂_{r,i} 三者之间的三个点积确定的。角点接触时角点首先计算粒子分别接触两条边a、b时 n̂_{f,i} 的方向记为 n̂_{f,a}、n̂_{f,b}。设这两条边的归一化线向量为 l̂_a 与 l̂_b取向指向所考虑的角点。若 n̂_{f,a} · l̂_b 0即该边的力有指向另一边的分量则 n̂_{f,a} 被替换为 n̂_{s,i}边b同理以避免力指向三角形内部。最终的 n̂_{f,i} 通过对这两个边贡献做加权平均得到使结果在粒子从接触一条边、经过角点、到接触另一条边的过程中在两个极限之间插值。两个边权重的计算较复杂简要说明如下以便理解代码上下文每个权重主要是若干点积的函数包括 C_a n̂_{r,i} · l̂_a、D_a n̂^p_{r,i} · l̂_a 以及边b的对应量其中 n̂^p_{r,i} 是 n̂_{r,i} 在三角形平面内的投影。若接触点恰好位于角点正上方n̂_{r,i} 平行于 n̂_{s,i}则直接令 n̂_{f,i} n̂_{s,i} 并跳过此过程。基础权重为 W_a (1 − C_a)·D_b、W_b (1 − C_b)·D_a。这一构造保证例如边a的权重在 n̂_{r,i} 与 l̂_a 对齐此时该边难以算出合理力或与 l̂_b 垂直使其贡献占主导时趋于零。评估角点接触时若两条边之一如a具有凸连接且正在隐藏另一三角形则该边需要额外权重以保证 n̂_{f,i} 在粒子跨过该凸转角移动到表面j时连续过渡。它以点积 E_a n̂^a_{r,i} · (δ_{s,i} δ_{s,j})/2 计算其中 n̂^a_{r,i} 是 n̂_{r,i} 减去沿 l̂_a 的分量后的向量。W_b 随后乘以 (1 − E_a)当粒子逼近切换到接触共享边a的表面j的阈值时该项趋于零。未连接自由边与角点的力粒子也会受到三角形未连接自由边或角点的力。例如若用三角形构造一个圆板未连接边构成板的周界。目前 LAMMPS 不计算也不存储沿此类边界的未连接边之间的连通性信息因此粒子沿其运动时力不保证以物理真实的方式变化例如边界有凸弯时粒子绕行时无法屏蔽远处边的力。在此约束下下述计算旨在优先保证力方向的连续性并确保力始终如预期地指向背离表面边界的方向。不过一般建议若表面存在裸露边界应将其设计得相对简单平滑。先假设复合集合中只有一个被接触的三角形边的力按与凸连接边相同的方式计算未连接角点至少有一条未连接边的角点也平均来自其两条边的贡献若两边均连接或均未连接此平均相同若一条连接、一条未连接则连接边如a的贡献会被修改。为在未连接与连接极限之间插值使用一般权重 W_c max(0.0, 1.0 − dr_uc / δ_max)其中 dr_uc 是复合集合内离三角形未连接角点/边的最大面内距离。若粒子位于集合有效表面正上方该权重为零粒子移离该有效表面时权重趋近 1。所有未在未连接特征处接触的表面的 W_i 随后乘以 W_c此外这些表面还将 W_i 乘以因子 W_ip min(1.0, n̂_{s,i} · n̂_{r,i} / (1 − W_c))该项在粒子移入表面平面内时趋于零实际上只应在粒子处于表面外侧且远离表面时发生避免内部连接无法判定接触表面哪一侧的复杂情况。最后当接触一个具有一条未连接边与一条连接边的未连接角点时连接边的贡献要乘以这两个因子 W_c 与 W_ip。几何设计建议Recommendations for Geometries设计颗粒表面几何时有几件事需要牢记以避免意外的力贡献并提高效率。凸转角与薄壁特征虽然凸角被识别并用于屏蔽壁面物理隐藏部分的力但如果粒子未接触凸转角的全部表面力就无法被正确屏蔽。例如考虑由 3 条线段定义的 2D U 形壁见 doc/src/img/gransurf_U_geom.png若 U 的宽度大于典型粒子-壁面重叠左预计无问题但若 U 的宽度相对典型重叠非常薄中、右粒子可能同时接触 U 的两条竖直腿。若粒子同时接触 U 的底边中绿色则可识别凸转角并屏蔽最右竖直腿蓝的力但若粒子位于 U 上部右且未接触底边绿则无法识别凸转角、无法屏蔽远处竖直壁蓝的力。因此间隙小于粒子-表面预期重叠距离的极薄壁特征可能导致意外的附加力。未连接边界如上节所述2D 中与线段未连接端点的接触力的方向符合预期且粒子绕过端点时无间断。但 3D 中与未连接边的接触只能产生大致背离边的合理方向力且接触沿一系列未连接边移动时力的确切方向可能摆动凸转角也可能无法被恰当屏蔽。因此建议避免沿未连接边界设计复杂几何例如相对于边长度的面内/面外快速振荡如褶皱或锯齿。邻表构建与性能为构建粒子与线段/三角形之间的邻表LAMMPS 为每个线段/三角形分配一个与包围该物体的圆/球半径对应的半径。因此须注意三角形/线段半径与粒子半径差异巨大的系统会减慢邻表构建可能受益于为局部表面使用 neighbor 的multi风格。此外高纵横比aspect ratio的三角形应尽量避免因为会导致很大的邻表包含许多实际并不接近接触的粒子误报。对应实现上全局表面邻表构建支持neighbor nsq与neighbor bin两种算法源码 fix_surface_global.cpp 中enum{NSQ, BIN}bin默认将原子与表面空间分箱后只检查邻近箱nsq则对所有原子与所有表面循环。经验法则是表面数量越多bin越快但若表面数量很大fix surface/local 可能更具性能优势。示例脚本Example Scriptsexamples/gransurf 目录包含同时使用global与local表面的示例输入脚本覆盖 2D 与 3D 模型。每个脚本使用 dump image 命令生成一系列快照图像同时可视化粒子与颗粒表面并可合成动画观看模拟过程。该目录同时包含surf.tri、surf.line、surf.tribox等表面定义文件、surf.container.stl、surf.helix.stl、surf.tri.stl等 STL 文件以及data.tri、data.line、data.tribox等 data 文件。全局表面 3D 示例molecule 输入examples/gransurf/in.tri.gran.global 演示向三角形构成的箱体倒入球体并让箱体旋转units lj atom_style sphere boundary f f f dimension 3 pair_style granular pair_coeff * * hooke 4000.0 100.0 tangential linear_history 1000.0 1.0 0.5 damping mass_velocity neighbor 0.3 bin molecule tris surf.tri fix 20 all surface/global input mol tris model * * granular hooke 4000.0 100.0 tangential linear_history 1000.0 1.0 0.5 damping mass_velocity fix 2 all nve/sphere fix 3 all gravity 1.0 vector 0 0 -1 region slab block -2 2 -2 2 7 9 units box fix ins all pour 100 1 4767548 vol 0.4 10 diam one 1.0 region slab ignore run 5000 fix_modify 20 move * rotate 0 0 0 0 0 1 25 run 5000 fix_modify 20 move * none run 5000其中molecule tris surf.tri定义表面文件 examples/gransurf/surf.tri包含 4 个三角形每个三角形行格式为index type molID x1 y1 z1 x2 y2 z2 x3 y3 z3fix_modify 20 move * rotate ...让所有表面绕 z 轴以角速度 25 旋转move * none则停止运动。全局表面 3D 示例STL 输入examples/gransurf/in.tri.gran.global.stl 与上例等价区别在于表面来自 STL 文件fix 20 all surface/global input stl 2 1 surf.tri.stl model * * granular hooke 4000.0 100.0 tangential linear_history 1000.0 1.0 0.5 damping mass_velocityinput stl 2 1 surf.tri.stl中2是赋予所有三角形的表面 type1是赋予的 molecule ID。局部表面 3D 示例data 文件输入examples/gransurf/in.tri.gran.data 演示 data 文件输入方式——三角形直接作为atom_style tri粒子从 examples/gransurf/data.tri 读入粒子-粒子与粒子-表面相互作用通过混合 pair 风格分开设置atom_style tri atom_modify map array boundary f f f dimension 3 comm_modify vel yes cutoff 12.5 read_data data.tri pair_style hybrid granular surf/granular pair_coeff 1 1 granular hooke 4000.0 100.0 tangential linear_history 1000.0 1.0 0.5 damping mass_velocity pair_coeff 1 2 surf/granular hooke 4000.0 100.0 tangential linear_history 1000.0 1.0 0.5 damping mass_velocity pair_coeff 2 2 none group sphere type 1 group tri type 2 fix fsl all surface/local fix 2 sphere nve/sphere fix 3 sphere gravity 1.0 vector 0 0 -1 fix 1 tri move rotate 0 0 0 0 0 1 25这里 data 文件 examples/gransurf/data.tri 同时包含 4 个 sphere 粒子Atoms节与 4 个三角形Triangles节三角形行格式为index x1 y1 z1 x2 y2 z2 x3 y3 z3。fix surface/local无input关键字即表示表面来自已读入的 data 文件fix 1 tri move rotate ...用 fix move 驱动tri组即所有三角形粒子旋转。局部表面 3D 示例molecule 输入examples/gransurf/in.tri.gran.surf 与上例结构相同但表面来自 molecule 文件atom_style tri ... molecule tris surf.tri fix 20 all surface/local input mol tris此时fix surface/local读取 examples/gransurf/surf.tri 中的 4 个三角形每个生成一个三角形粒子。局部表面螺旋输送器示例examples/gransurf/in.screwfeeder.local 综合演示了局部表面、STL 双输入、flat关键字与分子级运动控制units lj atom_style hybrid tri sphere 1 atom_modify map array boundary f f f dimension 3 comm_modify vel yes cutoff 1.0 pair_style hybrid granular surf/granular pair_coeff 1 * surf/granular hertz 1.0 1.0 tangential linear_history 1.0 1.0 0.1 pair_coeff 1 1 granular hertz 1.0 1.0 tangential linear_history 1.0 1.0 0.1 pair_coeff 2*3 2*3 none group sphere type 1 fix 1 all surface/local input stl 2 1 surf.container.stl input stl 3 2 surf.helix.stl flat 10.0 fix 2 sphere nve/sphere fix 3 sphere gravity 1e-5 vector 0 0 -1 fix 4 sphere pour 400 1 7892374 region insert_region diam range 0.5 1.0 group case molecule 1 group corkscrew molecule 2 fix 5 corkscrew move rotate 0 0 0 0 1 0 1000 units box关键点容器molecule ID 1与螺旋molecule ID 2作为两个独立 STL 输入flat 10.0将螺旋表面连接的最大平坦角阈值设为 10 度group case molecule 1/group corkscrew molecule 2按 molecule ID 分组后仅螺旋组被fix move驱动旋转——这是对同一物体相连表面须整体运动规则的正确应用。2D 线段表面示例examples/gransurf/in.line.gran.global 演示 2D 全局线段表面箱体直接旋转无需独立搅拌器dimension 2 atom_style sphere boundary f f p molecule lines surf.line fix 20 all surface/global input mol lines model * * granular hooke 4000.0 100.0 tangential linear_history 1000.0 1.0 0.5 damping mass_velocity fix_modify 20 move * rotate 0 0 0 0 0 1 25其中 examples/gransurf/surf.line 定义了 9 条线段格式index type molID x1 y1 x2 y2构成箱形边界整个边界绕 z 轴旋转模拟转鼓混料。该目录还有对应的局部 data 输入版本in.line.gran.data与局部 molecule 输入版本in.line.gran.surf。总结从选型到落地的决策路径选型三角形/线段数量少、允许大型表面 → 全局fix surface/global配atom_style sphere数量大、表面应细小均匀 → 局部fix surface/local配atom_style hybrid ... tri/line sphere输入优先考虑可携带 type/molecule ID 的 molecule 文件或 data 文件STL 文件需用stype smol显式赋予类型与分子 ID并留意退化、重复与细长三角形问题必要时用 PyMeshLab 重网格化连通性确保共享角点/端点精确一致molecule/STL或允许 EPSILON 不精确匹配data 文件同一物体共用 molecule ID赋予运动时须覆盖整个相连物体运动全局用fix_modify move按 molecule ID 或mol/region指定局部用fix move按 group 指定可用group ... molecule定义调试与可视化用 dump image 的fix关键字渲染表面2D 渲染为圆柱状线段3D 渲染为三角形/线框颜色与透明度可经dump_modify fcolor/ftrans调节用neighbor nsq/bin权衡邻表构建性能。赞分享科研科学计算高性能计算【免费下载链接】lammpsPublic development project of the LAMMPS MD software package项目地址https://gitcode.com/gh_mirrors/la/lammps点击查看免费下载相关推荐QuantLib波动率建模完全指南从Black模型到局部波动率曲面QuantLib波动率建模完全指南从Black模型到局部波动率曲面 在金融衍生品定价和风险管理领域 QuantLib波动率建模 是构建精确金融模型的核心技术金融科技科学计算全局流局部注意力(GFLA)模型使用指南全局流局部注意力 GFLA 模型使用指南 本指南旨在提供一个清晰的路径帮助开发者和研究人员快速了解并上手 全局流局部注意力 https://github.coNetworkX 聚类算法完全指南从三角形计数到广义度的网络局部结构分析NetworkX 聚类算法完全指南从三角形计数到广义度的网络局部结构分析 在复杂网络分析中物以类聚是一个核心观察社交网络中朋友的朋友往往也是朋友蛋白图计算数据分析科学计算上一篇M3u8Downloader_H批量下载指南3个步骤完成多视频M3U8下载下一篇Sketch MeaXure 智能标注插件指南3 分钟装好把设计交付交给快捷键创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考