LIGGGHTS与CFDEM耦合：DEM-CFD双向数据映射原理解析

简介一份关于LIGGGHTS与CFDEM耦合的专题技术报告PDF源自2013年DEM6国际会议作者来自奥地利林茨大学颗粒流体建模系与DCS Computing公司。报告从工业颗粒过程建模痛点切入指出颗粒物料的处理能耗与浪费问题并系统梳理了解析型与未解析型两套CFD-DEM建模方法以及粗粒化、MP-PIC等加速模拟策略。在工程实现部分报告详细介绍了分布式与共享内存并行化、基准测试和代码设计优化助力大规模颗粒系统的高性能仿真。应用场景涵盖钢铁制造、散装物料处理、环境工程、流化床、矿物加工与农业等领域对从事离散元仿真、多相流模拟和工业过程优化的科研人员与工程师均有较强参考价值。资源共1个PDF文件压缩包大小为5.15MB已有402人浏览学习内容浓缩了会议演讲的全部幻灯片包含大量原理图示、数据对比和参考文献线索可作为快速掌握LIGGGHTSCFDEM耦合技术架构与实际应用的导览资料。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 都行但耦合算例中颗粒源项会把动量方程变得更 stiffnCorrectors设为 2 才能让压力场收敛得足够好否则在颗粒密集区会出现负压力振荡。4.3 DEM 侧 LIGGGHTS 脚本与 couplingPropertiesLIGGGHTS 侧的最小脚本如下# 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% 的发散问题。本文还有配套的精品资源点击获取