简介：一份关于LIGGGHTS与CFDEM耦合的专题技术报告PDF，源自2013年DEM6国际会议，作者来自奥地利林茨大学颗粒流体建模系与DCS Computing公司。报告从工业颗粒过程建模痛点切入，指出颗粒物料的处理能耗与浪费问题，并系统梳理了解析型与未解析型两套CFD-DEM建模方法，以及粗粒化、MP-PIC等加速模拟策略。在工程实现部分，报告详细介绍了分布式与共享内存并行化、基准测试和代码设计优化，助力大规模颗粒系统的高性能仿真。应用场景涵盖钢铁制造、散装物料处理、环境工程、流化床、矿物加工与农业等领域，对从事离散元仿真、多相流模拟和工业过程优化的科研人员与工程师均有较强参考价值。资源共1个PDF文件，压缩包大小为5.15MB，已有402人浏览学习；内容浓缩了会议演讲的全部幻灯片，包含大量原理图示、数据对比和参考文献线索，可作为快速掌握LIGGGHTS+CFDEM耦合技术架构与实际应用的导览资料。
1. 从单相 DEM 到流固耦合：LIGGGHTS 与 CFDEM 到底在耦合什么
做颗粒输运、流化床或者料仓卸料仿真的人，迟早会撞上这么一个问题：离散元（DEM）能算出每个颗粒的碰撞和运动，但算不了气体或液体给颗粒的曳力；CFD 能算出流场的涡和压降，但把颗粒当成连续相处理又丢失了颗粒尺度的真实行为。LIGGGHTS 与 CFDEM coupling 的完整建模方案，就是在这个交叉点上被逼出来的——它把 LAMMPS 技术栈里的颗粒求解器（LIGGGHTS）和开源 CFD 求解器（OpenFOAM）通过 CFDEMcoupling 框架拼成一个双向数据回路：CFD 算流动，把速度和压力场映射到颗粒上；DEM 算碰撞和位移，把颗粒位置和速度映射回流场网格，更新孔隙率和动量交换源项。这个方案解决的典型问题是：料仓卸料时气流从料口倒灌、密相气力输送的堵塞预测、以及喷动床的颗粒循环频率。适合的人群是有 LAMMPS 或 OpenFOAM 基础、准备做颗粒-流体两相仿真的工程师和研究生。先说结论：这套耦合的关键不在两个求解器各自多强，而在它们之间的映射模型和耦合时间步的配合——这也是本文后面反复出现的主线。
2. 搞清耦合骨架：LAMMPS 技术栈里 LIGGGHTS 与 CFDEM 的分工边界
2.1 LIGGGHTS 为什么是 LAMMPS 的颗粒分支而非替代品
LIGGGHTS 是 LAMMPS 的一个开源扩展版本，全称 LAMMPS Improved for General Granular and Granular Heat Transfer Simulations。它保留了 LAMMPS 的基本架构——原子风格、邻居列表、时间积分器这些核心机制都在，但把原子类型替换成了适合颗粒力学的球体颗粒，并加入了 Hertzian 接触模型、JKR 粘附模型、颗粒传热模型等一套 DEM 专用工具。所以如果你已经写过 LAMMPS 的 in 脚本，看 LIGGGHTS 脚本会非常亲切：atom_style granular、pair_style granular、fix nve/sphere这些命令都是 LAMMPS 语法的直接延伸。
安装 LIGGGHTS 时常见的一个坑是把它当作 LAMMPS 的补丁来打——直接下载 LAMMPS 源码再加 LIGGGHTS 补丁包。这种"从 LAMMPS 做加法"的思路在三五年前还能走通，但现在的 LIGGGHTS 已经是一个独立维护的代码库，src目录下自带完整的颗粒求解器实现，不需要也不建议与上游 LAMMPS 混编。正确做法是把 LIGGGHTS 当作一个独立的可执行文件（通常命令行就是liggghts）来编译和调用，与 OpenFOAM 的耦合也全部通过这个可执行文件完成。
2.2 CFDEMcoupling 的求解器-颗粒引擎对接协议
CFDEMcoupling 不是一个独立软件，而是一套 OpenFOAM 的求解器框架和库集合。它在 OpenFOAM 的src目录下新增了一组耦合库，并预置了cfdemSolverPiso、cfdemSolverPimple等求解器，这些求解器本质上是在 OpenFOAM 标准 PISO/PIMPLE 算法里嵌入了对 LIGGGHTS 的调用逻辑。
耦合对接的核心是"分区迭代"而非"全隐式联立"。在每个耦合时间步内，CFD 求解器先推进一个时间步，算出流场；然后把流场信息（速度、压力、粘度）通过 MPI 或文件传递给 LIGGGHTS 进程；LIGGGHTS 根据这些信息计算颗粒受到的流体力，推进颗粒位置和速度；最后把新的颗粒数据映射回 CFD 网格，更新孔隙率和动量交换源项，进入下一个时间步。这个协议天然支持双向耦合，但代价是时间步长必须取 CFD 和 DEM 中较小者的约束——实际工程中通常会让 DEM 子循环多步，而 CFD 每步同步一次。
相互通信的数据通过共享内存和文件两种模式传递。文件模式适合调试，因为每一步的颗粒位置（particles文件）和耦合输入（cfdemCoupling文件）都落在磁盘上，出问题可以直接查看；共享内存模式适合生产计算，省去磁盘 IO 开销，但排查问题时要靠日志输出来定位。
2.3 耦合数据流：体积分数、曳力与三向映射
2.3.1 从颗粒到网格的映射
每一轮 DEM 计算结束后，CFDEM 要把颗粒位置"翻译"到 CFD 网格上。这个翻译包含两个量：孔隙率（void fraction）和动量交换源项。孔隙率计算是即使标题里只写了 coupling 也必须理解的核心概念——某个网格单元的体积中，被颗粒占据的体积占比是多少。CFDEM 提供了多种孔隙率模型，最常用的是中心点法（颗粒完全计入其中心所在的网格）和分割法（按颗粒与网格的重叠体积比例分摊）。
# 以颗粒中心点法为例的伪码，说明孔隙率映射逻辑 for cell in mesh_cells: void_fraction[cell] = 1.0 for particle in particles: cell_id = locate_particle_center(particle.position, mesh) # 颗粒体积除以所在网格体积，得到该颗粒贡献的固相体积分数 void_fraction[cell_id] -= particle.volume / mesh.cell_volume[cell_id]这段伪码展示了最朴素的映射策略：先定位颗粒中心所在网格，再用颗粒体积直接扣减该网格的孔隙率。locate_particle_center在 CFDEM 内部是通过网格索引树实现的，OpenFOAM 的mesh.findCell()负责这一步。分割法更精确，但每个颗粒需要遍历其覆盖的所有网格单元，计算几何重叠体积，计算量大约高出 3 到 5 倍。对于密相输送这类颗粒体积分数超过 30% 的场景，建议直接用分割法，因为中心点法在颗粒直径接近网格尺寸时会产生孔隙率剧烈跳变，导致压力场振荡。表面平滑处理也是工程中常用的手段，CFDEM 里的smoothing子字典可以指定对孔隙率场做几轮拉普拉斯平滑——注意，这只有在网格比颗粒大时才有效，网格比颗粒还细时平滑容易抹掉真实物理。
2.3.2 从网格到颗粒的映射
反向映射的主要内容：每个颗粒所在位置的流场速度、压力梯度和湍流粘度。常用的映射方式是从颗粒中心所在网格直接取单元中心值，不做插值。这样做在网格比颗粒细的时候会带来明显的数值噪声——颗粒跨越网格边界时，感受到的流场速度会发生阶跃。改进方案是距离加权插值（inverse distance weighting），取颗粒周围若干个网格的流场值按距离倒数加权平均。
曳力模型是第三个映射的承载者。CFDEM 里常见的曳力模型包括 WenYu、DiFelice、Gidaspow 和 KochHill，它们都是基于颗粒雷诺数和局部孔隙率的经验关联式。选型原则很简单：颗粒体积分数低于 10% 用 WenYu；高于 10% 用 Gidaspow；不知道选什么先用 DiFelice，它对孔隙率的依赖形式最平滑，不容易发散。无论选哪种，曳力最终会以体积力的形式加回 CFD 的动量方程源项，完成闭环。
3. 把 LAMMPS 技术底座的耦合环境装起来：LIGGGHTS 与 OpenFOAM 安装
3.1 版本匹配是安装的首要矛盾
LIGGGHTS 和 CFDEMcoupling 的版本与 OpenFOAM 版本之间存在严格的对应关系。社区里最稳定的组合是 LIGGGHTS 3.8.0 配合 CFDEMcoupling 3.8.0 和 OpenFOAM 5.x——这套组合被大量论文和算例验证过，编译报错最少。OpenFOAM 7 以上配合 CFDEMcoupling 4.x 也可以跑通，但需要手动处理一些头文件路径变化。安装前先在终端确认编译器版本。
# 检查编译工具链和 MPI 环境 g++ --version mpirun --version echo $WM_PROJECT_VERSION # 若已安装 OpenFOAM，会输出版本号3.2 编译安装的最小步骤
先装 OpenFOAM，再装 LIGGGHTS，最后编译 CFDEMcoupling——这个顺序不能乱。LIGGGHTS 的编译非常轻量，解压源码后直接 make 即可：
# LIGGGHTS 编译（3.8.0 版本为例） cd LIGGGHTS-PUBLIC-3.8.0/src make -j4 auto # 编译完成后生成 liggghts 可执行文件 ./liggghts -hmake auto是 LIGGGHTS 提供的自动检测选项，会探测 MPI 和编译环境。如果之后要让 LIGGGHTS 和 CFDEM 通过共享内存通信，还需要在make auto前额外导出CFDEM_COUPLING环境变量并启用-DCFDEM_COUPLING编译宏。这一步漏掉的话，后面 CFDEM 和 LIGGGHTS 之间无法建立连接，最常见的报错是cannot open shared object file。
CFDEMcoupling 的编译是通过其自带的Allwmake脚本完成的：
# 编译 CFDEMcoupling cd CFDEMcoupling-PUBLIC-3.8.0 ./Allwmake -j4Allwmake会依次编译耦合库和预置求解器。编译过程常见的失败点是找不到 LIGGGHTS 的头文件路径，解决办法是在环境变量CFDEM_LIGGGHTS_SRC_PATH中显式指定 LIGGGHTS 的源码目录。
3.3 安装后的连通性验证
安装完成后，做一次最小的连通性验证比直接上手算例要省时间得多。以下命令检查四个关键文件是否就位：
# 验证安装完整性 which cfdemSolverPiso # CFD 求解器 which liggghts # DEM 求解器 ls $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc # 耦合模板目录 ls $CFDEM_SRC_DIR/lagrangian/cfdemParticle/src # 耦合库源码如果cfdemSolverPiso命令找不到，多半是$CFDEM_SRC_DIR环境变量没有 source。登录 shell 的配置文件（~/.bashrc）里需要写入 OpenFOAM 的环境变量和 CFDEM 的环境变量，顺序是先 OpenFOAM 后 CFDEM。
4. 跑通耦合的最小算例与参数调优
4.1 最小算例的目录结构
一个完整的 CFDEM-DEM 耦合算例有三个组成部分：CFD 侧（OpenFOAM 标准算例结构）、DEM 侧（LIGGGHTS 输入脚本）、耦合配置文件（couplingProperties）。目录结构如下：
case/ ├── constant/ │ ├── couplingProperties/ # 耦合参数主文件 │ │ ├── couplingProperties │ │ ├── dragModel/ │ │ └── voidFractionModel/ │ ├── dynamicMeshDict # 动网格配置（静态算例可忽略） │ ├── transportProperties # 流体物性 │ └── turbulenceProperties # 湍流模型 ├── system/ │ ├── controlDict # 时间步长、库加载 │ ├── fvSolution # 求解器设置 │ └── fvSchemes # 离散格式 ├── constant/LIGGGHTS/ # DEM 脚本和颗粒模板 │ ├── in.liggghts │ └── mesh.stl └── 0/ # 初始场这个结构里couplingProperties是 CFDEM 区别于普通 OpenFOAM 算例的关键文件。没有任何一个标准 OpenFOAM 算例会有这个目录——它是 CFDEM 编译时通过system/controlDict里的libs字段加载耦合库后自动识别的配置入口。
4.2 CFD 侧求解器配置与耦合开关
controlDict中需要显式加载耦合库：
// system/controlDict 关键配置 application cfdemSolverPiso; libs ("libcfdemParticle.so"); // 加载耦合库 deltaT 0.0001; // CFD 时间步长（秒） writeInterval 0.01;deltaT的选择直接决定耦合稳定性。经验公式是 CFD 时间步长不超过颗粒穿过一个网格所需时间的 20%：
deltaT < 0.2 * hexCellSize / particleVelocity如果颗粒速度是 1 m/s，网格尺寸 1 mm，那么deltaT应小于 0.0002 秒。CFDEM 本身不会校验这个条件，超过之后的表现是颗粒体积分数场出现"幽灵颗粒"——颗粒移动过快，从一个网格直接穿越到另一个网格，中间没有中间步骤，导致孔隙率更新滞后。
fvSolution里需要关注的是 PISO 的迭代次数设置：
// system/fvSolution 关键配置 PISO { nCorrectors 2; // PISO 修正次数，颗粒耦合建议至少 2 nNonOrthogonalCorrectors 1; pRefCell 0; pRefValue 0; }nCorrectors决定每个时间步内压力修正方程被求解的次数。纯 CFD 算例设 1 或 2 都行，但耦合算例中颗粒源项会把动量方程变得更 stiff，nCorrectors设为 2 才能让压力场收敛得足够好，否则在颗粒密集区会出现负压力振荡。
4.3 DEM 侧 LIGGGHTS 脚本与 couplingProperties
LIGGGHTS 侧的最小脚本如下：
# in.liggghts —— DEM 侧最小耦合脚本 atom_style granular boundary m m m newton off communicate single vel yes region box block -0.005 0.005 -0.005 0.005 -0.01 0.01 units mm create_box 1 box create_atoms 1 single 0 0 -0.008 units mm pair_style gran hertz tangential history cohesion sjkr pair_coeff * * 100000 0.3 0.4 0.5 1.0e-6 0.005 1.0 fix g1 all gravity 9.81 vector 0 0 -1 fix cfd all couple cfdCoupling # 开启与 CFDEM 的耦合 fix integr all nve/sphere run 1脚本中fix cfd all couple cfdCoupling是耦合的开关命令，它让 LIGGGHTS 在每一步计算前等待 CFDEM 传入的曳力和压力梯度力，计算完颗粒运动后再把颗粒坐标返回给 CFDEM。pair_style gran hertz tangential history后面的参数是等效杨氏模量、恢复系数、摩擦系数、粘附系数和截止距离——这些值取决于你的颗粒材料，钢球和玻璃珠的杨氏模量差了三个数量级，不要直接照搬。
耦合的核心配置在constant/couplingProperties/couplingProperties文件中：
// constant/couplingProperties/couplingProperties couplingProperties { dragModel DiFelice; // 曳力模型 voidFractionModel default; // 孔隙率模型：中心点法 smoothing { smoothingStep 2; // 孔隙率平滑次数，0 为不平滑 } maxNumberOfParticles 1000000; // DEM 模拟的最大颗粒数上限 numberOfParticles 1000; // 实际颗粒数 }maxNumberOfParticles和numberOfParticles的差值会让 CFDEM 预分配颗粒数组的内存。保守做法是让上限设为实际颗粒数的 1.5 倍，太大占用内存，太小会在颗粒增多时直接报错退出。
4.4 3 个必调耦合参数
第一个是deltaT的比值。CFDEM 支持 CFD 步与 DEM 步不同步推进，通过DEM子字典控制：
DEM { deltaT 0.00001; // DEM 时间步长 numberOfSubSteps 10; // 每个 CFD 步内 DEM 推进 10 步 }numberOfSubSteps的计算依据是颗粒接触时间。Hertz 接触的接触时间约为颗粒半径除以瑞利波速的量级，DEM 时间步必须小于接触时间的 1/10。如果你的颗粒弹性模量很高，接触时间极短，就需要增大子步数——这是耦合算例里最常见的 DEM 侧发散来源。
第二个是曳力模型的参数修正。DiFelice模型需要指定颗粒雷诺数的参考尺度：
dragModel { DiFelice { C1 0.63; C2 4.8; } }C1和C2是 DiFelice 曳力关联式中的经验常数，默认值适配球形颗粒，如果颗粒是非球形或粗糙表面，需要根据实验沉降速度反推校准。
第三个是孔隙率模型的平滑系数。前面提过smoothingStep，这个参数对压力场稳定性影响巨大。颗粒直径与网格尺寸比在 1:2 到 1:3 之间时，建议从smoothingStep 1起步，逐步增加，直到压力场不再出现锯齿状振荡。注意：平滑意味着孔隙率场的分辨率下降，平滑次数过多时颗粒在 CFD 网格里看起来像一团模糊的云，失去 DEM 的意义。
5. 收在工程技巧上：后处理与耦合稳定性判断
5.1 颗粒位置与流场的时空对齐检查
耦合计算结束后的第一步不是画云图，而是核对颗粒位置和流场在时间轴上是否真的对齐了。CFDEM 默认把颗粒信息写在postProcessing/liggghts目录下，而 CFD 场量写在时间步目录中。检查的方法是取同一时间步的颗粒坐标，确认它们落在 CFD 计算域的有效网格内——如果颗粒坐标超过了几何边界且没有触发碰撞，多半是 DEM 侧时间步太大导致颗粒穿透了壁面。另一个常见问题是文件输出的时间精度：CFDEM 的颗粒输出时间步长由controlDict的writeInterval决定，而颗粒位置写出的瞬间并不总是与 CFD 场完全同步，后处理时要做插值对齐。
5.2 速度松弛与湍流耦合的稳定性联动
工程中做气固两相流时，CFD 的湍流模型对耦合稳定性影响非常大。用kEpsilon湍流模型做密相颗粒流，颗粒源项会把湍动能方程顶到发散。常见做法是改用力学粗糙度更低的层流假设——把turbulenceProperties里的simulationType设为laminar，先验证耦合逻辑，再加湍流模型。这个方法在流化床和料仓卸料算例中都适用，不属于物理妥协，而是分步验证的手段。另外一个联动因素是速度松弛（velocity relaxation）：在couplingProperties里可以给颗粒速度加一个松弛因子，让颗粒速度更新时保留一部分旧值。松弛因子在 0.3 到 0.7 之间，值越小越稳定但收敛越慢。启动阶段建议从 0.5 起步，颗粒流场建立后再逐步调到 0.8 以上。
5.3 判断耦合发散的两类直接信号
耦合发散与纯 CFD 发散的信号不尽相同。第一类是压力场出现棋盘式振荡，即相邻网格单元的压力值正负交替，且振荡区域恰好与颗粒密集区重合。这时优先检查孔隙率映射是否过于粗糙——把voidFractionModel从default切换到divided（分割法），同时增加smoothingStep。第二类是 DEM 侧的温度或动能爆炸式增长——颗粒速度在几个时间步内跳增几个数量级，对应的日志里会出现NaN或inf。这时检查点集中在numberOfSubSteps是否太小，以及pair_style里的刚度系数是否与deltaT匹配。一个便捷的自检方法是把耦合关掉，单独用 LIGGGHTS 跑相同的颗粒碰撞算例，如果单跑 DEM 自己能稳定，问题就出在 CFDEM 映射的曳力上；如果单跑 DEM 也不稳定，问题回到 DEM 时间步长本身。这两类信号配合两个侧面的排查方向，基本上能覆盖耦合算例里 80% 的发散问题。
本文还有配套的精品资源，点击获取
