简介LIGGGHTS与CFDEM耦合是基于LAMMPS技术扩展的颗粒-流体耦合模拟方案主要面向颗粒技术、散体物料与多相流方向的工程师和研究人员解决宏观颗粒过程建模、仿真效率及工业放大中的实际问题。这份PDF是DEM6会议的技术报告讲义压缩包内共1个PDF文件整体大小约5.15MB内容紧凑适合作为技术入门或快速查阅。讲义从工业视角引入颗粒处理的普遍性与能耗浪费问题系统梳理了解析型与未解析型CFD-DEM建模策略、粗粒化与MP-PIC方法并围绕并行化、代码设计与基准测试讨论可扩展性和可维护性还覆盖钢铁冶炼、散装物料处理、流化床、环境工程、矿物加工及农业等典型应用。读者可借此完整了解LIGGGHTSCFDEM的框架组成、建模选择与常见应用场景目前已有402人学习浏览是一份面向耦合仿真的精炼参考资料。1. 用 LAMMPS 技术栈做的 CFD-DEM 耦合到底解决什么问题一条反直觉的结论先放在前面LIGGGHTS 并不是给 LAMMPS 装个颗粒插件而是直接 fork 了 LAMMPS 的分子动力学内核把原子间的短程势替换成颗粒接触力学再让它和 OpenFOAM 通过 CFDEM 耦合层做内存级数据交换。CFDEM 耦合建模要做的事情是在同一个模拟里让流体场和颗粒场互相演化颗粒受曳力、压力梯度力被流体推动反过来通过孔隙率和动量源项改写流体的 N-S 方程。这种双向反馈是流化床、气力输送、料仓卸料这些场景的核心也是纯 CFD 的欧拉-欧拉双流体模型难以精确描述的场合。这套方案对两类人最有用。一类是做散料输送和设备优化的工程师需要预测颗粒的轨迹、堆积角和磨损分布另一类是研究颗粒两相流的科研人员需要把自己的接触模型或曳力模型快速写进框架里验证。LAMMPS 的技术底子保证了颗粒求解器的高性能和可扩展性OpenFOAM 的开源生态则让流体侧的离散格式、湍流模型和网格工具直接可用。接下来的内容按架构→安装→算例→调参的顺序展开每一步都给出可以直接执行的命令和参数依据。2. 从 LAMMPS 到 LIGGGHTSCFDEM 耦合层到底交换了什么LIGGGHTS 的命名已经说明了它的血统LAMMPS Improved for General Granular and Granular Heat Transfer Simulations。它在代码层面保留了 LAMMPS 的 neighbor 列表、域分解和 MPI 通信框架但把 atom_style atomic 换成了 granular 专用风格。这里有个关键差异必须点破分子动力学里的原子碰撞是瞬时的二体事件而颗粒堆是持久的多体接触——一个颗粒可以同时被六个邻居压着接触力需要历史信息切向位移、法向重叠量所以 LIGGGHTS 重写了力计算的内核而不是复用 LAMMPS 的 pair_style。2.1 颗粒侧atom_style granular 与接触力模型的对应关系在 in.liggghts 文件里建模的第一步就是声明颗粒的属性维度。与 LAMMPS 安装完默认的原子风格不同granular 风格在创建原子时就要带上半径、密度和取向信息atom_style granular atom_modify map array create_box 1 box boundary create_atoms 1 random 1000 12345 box units box set group all density 2500 # 颗粒密度流化床石英砂常用 2500 kg/m3 set group all diameter 0.0005 # 颗粒直径 0.5 mmset group all diameter是 LIGGGHTS 扩展出来的命令LAMMPS 原生没有。它会在内部改写每个颗粒的 radius、mass 和初始转动惯量。接触模型通过 pair_style 选择常用的三种写法和适用场景如下pair_style 写法力学模型适合场景pair_style gran hertz/historyHertz-Mindlin 法向 切向历史接触干颗粒、无粘性堆积、流化床的默认选择pair_style gran jkrJohnson-Kendall-Roberts 粘附力湿颗粒、粉尘团聚、微米级颗粒pair_style gran hooke/history线性弹簧-阻尼调试期快速验证、接触刚度未知的预研代码块里pair_coeff * *后面的参数需要手动指定至少要给剪切模量、泊松比、恢复系数和摩擦系数pair_style gran hertz/history pair_coeff * * 1e7 0.45 0.5 0.5 # 剪切模量 1e7 Pa泊松比 0.45恢复系数 0.5滑动摩擦 0.5注意pair_coeff的四组数值不写关键字LIGGGHTS 按位置解析参数。剪切模量不要直接拿真实石英砂的 30 GPa 去算那会让瑞利波时间步落到纳秒级数天内跑不完。工程上常用软化接触技巧把剪切模量降低两到三个数量级只要保证颗粒重叠量不超过粒径的 1%结果仍具有工程精度。2.2 流体侧带孔隙率的 N-S 方程与动量交换源项CFDEM 耦合不是简单地把两个求解器轮流跑它的建模基础是局部平均化的 N-S 方程。OpenFOAM 这边求解的是带颗粒体积分数的控制方程∂(α_f ρ_f u_f)/∂t ∇·(α_f ρ_f u_f u_f) -α_f ∇p ∇·(α_f τ) α_f ρ_f g - F_fp其中 α_f 是流体体积分数F_fp 是颗粒对流体的反作用力源项。这个方程写出来耦合的物理意图就清楚了颗粒占掉的体积要从孔隙率里扣掉颗粒受到的流体合力要以相反方向还回流体动量方程。alpha场的计算发生在每个耦合时间步的开始CFDEM 把颗粒的包围球体积分摊到重叠的网格单元内累计得到alpha再传递给 OpenFOAM 的孔隙率场。曳力是这个方程里最需要斟酌的项。CFDEM 针对不同颗粒雷诺数区间提供了多种曳力模型建模时必须根据流态选择Gidaspow 模型组合了 WenYu 稀相公式与 Ergun 密相公式在流化床的鼓泡段表现最好是 coupling 算例的默认选择WenYu 模型只适用颗粒雷诺数小于 1000 的稀相流颗粒体积分数超过 0.5 时会严重低估曳力DiFelice 模型用一个体积分数修正指数覆盖从稀相到密相的过渡适合宽粒径分布的气力输送。选择依据不是哪个更先进,而是你的颗粒浓度范围在哪个区间。如果算例里某些区域的 α_f 低到 0.4 以下Gidaspow 的 Ergun 项是唯一能给出合理压降的。2.3 双向耦合的时间序列与守恒关系CFDEM 采用先流体后颗粒的顺序。一个耦合循环里OpenFOAM 先推进一个流体时间步得到收敛的速度场和压力场然后 CFDEM 通过插值把流场信息映射到每个颗粒位置计算曳力和压力梯度力接着 LIGGGHTS 把这些力加到颗粒上推进若干个子时间步最后把颗粒对流体产生的反作用力回传给 OpenFOAM 的动量源项。整个循环内数据都存在内存里不经过磁盘文件这也是 CFDEM 比自写的文本交换式耦合快得多的原因。力守恒的检查口径在两侧是不同的。OpenFOAM 侧看的是动量方程残差是否降到 1e-7 以下LIGGGHTS 侧看的是run结束时thermo输出的动能统计。CFDEM 在日志里会每步打印一条力平衡信息包含颗粒总受力和流体受到的反馈力两者的差值应该在浮点误差范围内。实际操作时不要只盯收敛残差要同时观察这条差值——它才是耦合正确性的直接证据。3. 把 LIGGGHTS、CFDEM 和 OpenFOAM 在本地装通的标准做法LIGGGHTS 的安装套路和 LAMMPS 安装大体一致迷你版说辞是选对 MPI、选对编译器、再在源码里开 cfdem 耦合开关。但 CFDEM 不是 LAMMPS 的一个内置包它是一套独立编译的 OpenFOAM 库组运行时才把 LIGGGHTS 的共享库链接进来。这决定了安装顺序必须严格先编译 LIGGGHTS再编译 CFDEMcoupling最后用自带算例验证。3.1 版本匹配的三角关系CFDEM 的代码对 OpenFOAM 内部类名的变动很敏感不同分支对应不同 OpenFOAM 版本。公开仓库的分支命名规则一般是CFDEM-版本号-OpenFOAM版本号比如针对 OpenFOAM 8 或 ESI 版的分支。在开始编译之前先确认三件事编成一张表组件推荐取值说明编译器gcc 8.3 及以上老版本对 C17 支持不完整wmake 会报模板错误MPIOpenMPI 4.xMPICH 也可以但和 OpenFOAM 的 compat 层偶有冲突OpenFOAM 版本Foundation 版或 ESI 版按你拿到的 CFDEM 分支决定不可混用# 进入 OpenFOAM 环境以 ESI 版为例 source /opt/openfoam8/etc/bashrc # 设置 CFDEM 相关环境变量 export CFDEM_VERSION6.0 export CFDEM_DIR$HOME/CFDEM/CFDEMcoupling-PUBLIC export LIGGGHTS_DIR$HOME/CFDEM/LIGGGHTS-PUBLICCFDEM_VERSION这个环境变量会被编译脚本用来生成库文件的版本后缀。如果你用的是某个较早的 release 分支版本号要写成分支名里的数字。判断标准很简单进入 CFDEMcoupling 源码目录看src/Makefile里引用$(CFDEM_VERSION)的地方对照写即可。3.2 编译 LIGGGHTS 并确认 cfdem 命令组存在LIGGGHTS-PUBLIC 源码根目录下就是标准的 LAMMPS 风格 src 目录编译目标名为auto它会自动探测 MPI 并生成带共享库支持的二进制cd $LIGGGHTS_DIR/src make -j4 auto # 安装到 LIGGGHTS 的默认 bin 目录 sudo make installauto目标的关键在于它默认开启了-fPIC编译选项CFDEM 的库在运行时需要以动态方式加载 LIGGGHTS 的符号表。如果用make serial这类目标编译后续链接会报 undefined reference 错误。编译结束后检查生成的二进制里是否包含耦合命令$LIGGGHTS_DIR/src/lmp_auto -h | grep cfdCoupling输出里能看到fix cfdCoupling这条命令才算通过。这一步很多人会漏掉装上的是一个普通 LAMMPS后面 CFDEM 一启动就提示无法识别命令。LIGGGHTS 从 LAMMPS 继承了大量分子动力学命令所以-h输出的绝大部分内容看起来和 LAMMPS 一模一样只有多出来的fix cfdCoupling、fix cfdDataExchange这些命令才是耦合真正依赖的扩展。3.3 编译 CFDEMcoupling 与安装验证CFDEMcoupling 的源码用 OpenFOAM 的 wmake 工具链编译不写 Makefile.am而是沿用 OpenFOAM 的Allwmake脚本规范cd $CFDEM_DIR/src ./Allwmake -j4编译过程会生成两类产物。一类是 OpenFOAM 的库文件放进$FOAM_USER_LIBBIN另一类是可执行求解器比如cfdemSolverPiso、cfdemSolverPimple放进$FOAM_USER_APPBIN。编译日志里如果出现wmake error: no such file ... foamVersion.H八成是前三步的source /opt/openfoam8/etc/bashrc没执行干净重新打开一个 shell 按顺序执行。# 验证库是否被正确加载 cfdemSolverPiso -help # 跑一个自带的最小算例做端到端验证 cd $CFDEM_DIR/tutorials ls自带的 tutorials 目录会列出多个算例挑一个名为mpicrudedBed或类似的最小流化床用例跑几分钟。正常情况 stdout 里会出现耦合日志、alpha场的平均值和 LiGggHTs 的每步颗粒信息。出现任何段的FOAM FATAL ERROR先看是出现在流体求解阶段还是颗粒推进阶段——前者多半是 OpenFOAM 版本不匹配后者多半是 LIGGGHTS 编译时没开共享库。4. 写一个最小可复现的流化床耦合算例装通的下一步不应该是立刻改大算例而是把一个最小流化床从零跑起来。这里给出一个 0.01m × 0.01m × 0.02m 的二维等效盒式算例颗粒直径 0.5 mm数量 1000 个气体从底部进入。这个规模在 4 核机器上十分钟内能跑完一个流化周期适合做耦合配置的调试平台。4.1 网格粗细与颗粒直径的比例约束CFDEM 对网格有一个硬性经验约束网格单元的特征尺寸要大于颗粒直径的 3 倍最好在 3 到 5 倍之间。原因是alpha场的计算需要对每个网格单元做颗粒体积的几何叠加颗粒太大、单元太小会让某个单元同时被好几个颗粒的边界跨过导致体积分数出现锯齿状分布进而让曳力振荡。# constant/polyMesh/blockMeshDict 的关键块 vertices ( (0 0 0) (0.01 0 0) (0.01 0.01 0) (0 0.01 0) (0 0 0.02) (0.01 0 0.02) (0.01 0.01 0.02) (0 0.01 0.02) ); blocks ( hex (0 1 2 3 4 5 6 7) (5 5 20) simpleGrading (1 1 1) );这里 5 × 5 × 20 的网格把每个方向单元尺寸拉到 2 mm是颗粒直径的 4 倍处于推荐区间。网格轴向 (z 向) 分层 20 层负责捕捉气泡和床层膨胀水平方向 5 × 5 的粗网格说明这个算例只追求定性正确和流程跑通不是定量研究定量研究需要把这三个方向的网格同时加密。4.2 双份求解器配置couplingProperties 与 in.liggghtsCFDEM 的核心配置在constant/couplingProperties文件里它用 OpenFOAM 的字典格式定义了耦合选项。最关键的三个条目是曳力模型、插值方案和颗粒数据交换方式particleCloud { couplingProperties { solverModel HI; dragModel Gidaspow; turbulenceModel off; regularization off; velocityInterpolationType cellCenter; } }solverModel HI表示采用 Hirt 和 Ishii 的隐式耦合方法把曳力对颗粒速度的依赖隐式处理到流体动量方程里这对密相颗粒流是必要的。velocityInterpolationType cellCenter使用单元中心处的流场值直接赋给颗粒计算最快但会在网格较粗的空间产生阶梯效应更平滑的选项是cellCenter对应的梯度插值调试初期不必追求。颗粒侧由constant/DEM目录下的in.liggghts文件描述它的运行方式是由 CFDEM 在每次耦合循环内调用 LIGGGHTS 推进颗粒子步atom_style granular boundary m m m newton off neighbor 0.003 bin neigh_modify delay 0 every 1 check yes region box block 0 0.01 0 0.01 0 0.02 units box create_box 1 box create_atoms 1 random 1000 12345 box units box set group all density 2500 set group all diameter 0.0005 pair_style gran hertz/history pair_coeff * * 1e7 0.45 0.5 0.5 timestep 1e-6 fix gravity all gravity 9.81 vector 0 0 -1 # 先自由下落形成堆积床 run 5000 # 开启与 CFDEM 的双向耦合 fix cfd all cfdCoupling run 50000颗粒下落 5000 步即 5ms 物理时间在 0.01m 高的盒子里足够重力和接触力取得平衡形成稳定的初始床层。随后fix cfd all cfdCoupling把 LIGGGHTS 挂到 CFDEM 的调度循环上这时开始跑的 50000 步由 CFDEM 在每个耦合步里以子步方式驱动。注意timestep容不得粗心1e-6 秒对于 0.5mm 颗粒配合 1e7 Pa 的剪切模量是安全的如果把剪切模量改成 1e8 Pa瑞利波时间步会缩短一个量级1e-6 就不再收敛表现为颗粒互相穿透或能量爆炸。4.3 启动命令与运行时日志的判断顺序blockMesh -case . # 准备耦合子目录CFDEM 会在运行时创建颗粒输出 mkdir -p postProcessing/dem # 启动耦合求解器 cfdemSolverPiso -case . log.cfd 21 cfdemSolverPiso是 CFDEM 提供的压力速度耦合求解器名字里的 PISO 表明它用 PISO 算法推进流体。启动后肉眼判断运行状态要按三层顺序看第一层看log.cfd最后几行有没有FOAM FATAL ERROR第二层看进程是否还活着CFDEM 在 LIGGGHTS 无输出时会显得像死锁实际是颗粒子步在跑第三层看固定的时间间隔有没有生成postProcessing下的场输出。用tail -f log.cfd观察每步的alpha平均值如果数值一开始接近 0.061000 个颗粒体积除以盒体总容积,随后随着床层膨胀逐渐下降说明耦合的空间插值链路是通的。5. 把耦合调稳的三个参数与三个常见坑耦合参数之间有强烈的联动关系单独调一个很难奏效。核心是先确定 DEM 子步与 CFD 时间步的比例再确定每个耦合间隔内 CFD 走几步最后用松弛因子把两个求解器的相互反馈压住。参数推荐取值调参信号DEM 时间步1e-6 s0.5mm 颗粒颗粒穿透边长或总动能发散CFD 时间步1e-4 s耦合步间隔 10库朗数超过 1 就减半曳力松弛因子0.6 ~ 0.8压力场逐周期振荡就往下调网格/粒径比3 ~ 5alpha 场呈棋盘分布就减小颗粒半径或加密网格第一个多发的坑是颗粒穿透量现象是颗粒在重力作用下漏出计算域底部。原因通常不是接触模型参数问题而是boundary m m m三个方向都用了周期性边界盒子底部没有物理墙面。流化床的底面必须是固定壁面在 in.liggghts 里要显式添加几何墙fix wallbase all wall/gran model hertz tangential history primitive type 1 zplane 0这行fix在 z0 处定义了一个无限平面墙颗粒与它的接触力由与颗粒对颗粒接触完全相同的 Hertz-Mindlin 模型计算。忘了加这行颗粒就会像穿过空气一样漏走。第二个坑是压力场出现周期性振荡但速度残差却收敛得不错。这通常发生在把 CFD 时间步从 1e-5 直接加到 5e-4 之后CFD 步长过大一个耦合步内曳力源项的变化幅度超过了流体压力修正能处理的量级PISO 的外循环数已经压不住。CFDEM 的couplingProperties里对源项做松弛处理couplingProperties { underRelaxation 0.7; momentumTransferModel kTGF; }underRelaxation取 0.7 意味着每次耦合循环只把新的动量源项和上一步的值做七三开切断振荡的正反馈。如果振荡仍然明显先把 CFD 时间步降低两倍再调松弛顺序不可颠倒。第三个坑不是发散而是假收敛计算结果看起来很正常但床层膨胀高度明显低于经验值。一个常见诱因是velocityInterpolationType用了默认的cellCenter把包含颗粒的单元中心速度直接赋给了颗粒。当网格偏粗时这个插值会低估颗粒位置处的真实流速等效于曳力打了折扣。把插值方式换成cellCenter对应的梯度插值选项或在阻力模型上改成DiFelice做交叉验证两者的差值能反映出插值误差的边界。验证耦合是否正确比看收敛残差更有效的是做力平衡核对在log.cfd里打印每个耦合步的全场曳力对颗粒的作用力与相同入口条件下纯压降乘以床层截面积得到的总压差力比较。一个稳态流化床的压降应当约等于床层颗粒的总浮重扣除浮力误差超过 15% 时优先检查alpha场而非曳力模型。最终把dragModel依序从 Gidaspow 换成 WenYu 再换成 DiFelice观察同一边界条件下的压降和床层空隙率曲线差异就能定量看到模型假设对耦合结果的影响边界。本文还有配套的精品资源点击获取 CFDEM耦合模拟深度解析：从LIGGGHTS与OpenFOAM原理到流化床算例实践
简介LIGGGHTS与CFDEM耦合是基于LAMMPS技术扩展的颗粒-流体耦合模拟方案主要面向颗粒技术、散体物料与多相流方向的工程师和研究人员解决宏观颗粒过程建模、仿真效率及工业放大中的实际问题。这份PDF是DEM6会议的技术报告讲义压缩包内共1个PDF文件整体大小约5.15MB内容紧凑适合作为技术入门或快速查阅。讲义从工业视角引入颗粒处理的普遍性与能耗浪费问题系统梳理了解析型与未解析型CFD-DEM建模策略、粗粒化与MP-PIC方法并围绕并行化、代码设计与基准测试讨论可扩展性和可维护性还覆盖钢铁冶炼、散装物料处理、流化床、环境工程、矿物加工及农业等典型应用。读者可借此完整了解LIGGGHTSCFDEM的框架组成、建模选择与常见应用场景目前已有402人学习浏览是一份面向耦合仿真的精炼参考资料。1. 用 LAMMPS 技术栈做的 CFD-DEM 耦合到底解决什么问题一条反直觉的结论先放在前面LIGGGHTS 并不是给 LAMMPS 装个颗粒插件而是直接 fork 了 LAMMPS 的分子动力学内核把原子间的短程势替换成颗粒接触力学再让它和 OpenFOAM 通过 CFDEM 耦合层做内存级数据交换。CFDEM 耦合建模要做的事情是在同一个模拟里让流体场和颗粒场互相演化颗粒受曳力、压力梯度力被流体推动反过来通过孔隙率和动量源项改写流体的 N-S 方程。这种双向反馈是流化床、气力输送、料仓卸料这些场景的核心也是纯 CFD 的欧拉-欧拉双流体模型难以精确描述的场合。这套方案对两类人最有用。一类是做散料输送和设备优化的工程师需要预测颗粒的轨迹、堆积角和磨损分布另一类是研究颗粒两相流的科研人员需要把自己的接触模型或曳力模型快速写进框架里验证。LAMMPS 的技术底子保证了颗粒求解器的高性能和可扩展性OpenFOAM 的开源生态则让流体侧的离散格式、湍流模型和网格工具直接可用。接下来的内容按架构→安装→算例→调参的顺序展开每一步都给出可以直接执行的命令和参数依据。2. 从 LAMMPS 到 LIGGGHTSCFDEM 耦合层到底交换了什么LIGGGHTS 的命名已经说明了它的血统LAMMPS Improved for General Granular and Granular Heat Transfer Simulations。它在代码层面保留了 LAMMPS 的 neighbor 列表、域分解和 MPI 通信框架但把 atom_style atomic 换成了 granular 专用风格。这里有个关键差异必须点破分子动力学里的原子碰撞是瞬时的二体事件而颗粒堆是持久的多体接触——一个颗粒可以同时被六个邻居压着接触力需要历史信息切向位移、法向重叠量所以 LIGGGHTS 重写了力计算的内核而不是复用 LAMMPS 的 pair_style。2.1 颗粒侧atom_style granular 与接触力模型的对应关系在 in.liggghts 文件里建模的第一步就是声明颗粒的属性维度。与 LAMMPS 安装完默认的原子风格不同granular 风格在创建原子时就要带上半径、密度和取向信息atom_style granular atom_modify map array create_box 1 box boundary create_atoms 1 random 1000 12345 box units box set group all density 2500 # 颗粒密度流化床石英砂常用 2500 kg/m3 set group all diameter 0.0005 # 颗粒直径 0.5 mmset group all diameter是 LIGGGHTS 扩展出来的命令LAMMPS 原生没有。它会在内部改写每个颗粒的 radius、mass 和初始转动惯量。接触模型通过 pair_style 选择常用的三种写法和适用场景如下pair_style 写法力学模型适合场景pair_style gran hertz/historyHertz-Mindlin 法向 切向历史接触干颗粒、无粘性堆积、流化床的默认选择pair_style gran jkrJohnson-Kendall-Roberts 粘附力湿颗粒、粉尘团聚、微米级颗粒pair_style gran hooke/history线性弹簧-阻尼调试期快速验证、接触刚度未知的预研代码块里pair_coeff * *后面的参数需要手动指定至少要给剪切模量、泊松比、恢复系数和摩擦系数pair_style gran hertz/history pair_coeff * * 1e7 0.45 0.5 0.5 # 剪切模量 1e7 Pa泊松比 0.45恢复系数 0.5滑动摩擦 0.5注意pair_coeff的四组数值不写关键字LIGGGHTS 按位置解析参数。剪切模量不要直接拿真实石英砂的 30 GPa 去算那会让瑞利波时间步落到纳秒级数天内跑不完。工程上常用软化接触技巧把剪切模量降低两到三个数量级只要保证颗粒重叠量不超过粒径的 1%结果仍具有工程精度。2.2 流体侧带孔隙率的 N-S 方程与动量交换源项CFDEM 耦合不是简单地把两个求解器轮流跑它的建模基础是局部平均化的 N-S 方程。OpenFOAM 这边求解的是带颗粒体积分数的控制方程∂(α_f ρ_f u_f)/∂t ∇·(α_f ρ_f u_f u_f) -α_f ∇p ∇·(α_f τ) α_f ρ_f g - F_fp其中 α_f 是流体体积分数F_fp 是颗粒对流体的反作用力源项。这个方程写出来耦合的物理意图就清楚了颗粒占掉的体积要从孔隙率里扣掉颗粒受到的流体合力要以相反方向还回流体动量方程。alpha场的计算发生在每个耦合时间步的开始CFDEM 把颗粒的包围球体积分摊到重叠的网格单元内累计得到alpha再传递给 OpenFOAM 的孔隙率场。曳力是这个方程里最需要斟酌的项。CFDEM 针对不同颗粒雷诺数区间提供了多种曳力模型建模时必须根据流态选择Gidaspow 模型组合了 WenYu 稀相公式与 Ergun 密相公式在流化床的鼓泡段表现最好是 coupling 算例的默认选择WenYu 模型只适用颗粒雷诺数小于 1000 的稀相流颗粒体积分数超过 0.5 时会严重低估曳力DiFelice 模型用一个体积分数修正指数覆盖从稀相到密相的过渡适合宽粒径分布的气力输送。选择依据不是哪个更先进,而是你的颗粒浓度范围在哪个区间。如果算例里某些区域的 α_f 低到 0.4 以下Gidaspow 的 Ergun 项是唯一能给出合理压降的。2.3 双向耦合的时间序列与守恒关系CFDEM 采用先流体后颗粒的顺序。一个耦合循环里OpenFOAM 先推进一个流体时间步得到收敛的速度场和压力场然后 CFDEM 通过插值把流场信息映射到每个颗粒位置计算曳力和压力梯度力接着 LIGGGHTS 把这些力加到颗粒上推进若干个子时间步最后把颗粒对流体产生的反作用力回传给 OpenFOAM 的动量源项。整个循环内数据都存在内存里不经过磁盘文件这也是 CFDEM 比自写的文本交换式耦合快得多的原因。力守恒的检查口径在两侧是不同的。OpenFOAM 侧看的是动量方程残差是否降到 1e-7 以下LIGGGHTS 侧看的是run结束时thermo输出的动能统计。CFDEM 在日志里会每步打印一条力平衡信息包含颗粒总受力和流体受到的反馈力两者的差值应该在浮点误差范围内。实际操作时不要只盯收敛残差要同时观察这条差值——它才是耦合正确性的直接证据。3. 把 LIGGGHTS、CFDEM 和 OpenFOAM 在本地装通的标准做法LIGGGHTS 的安装套路和 LAMMPS 安装大体一致迷你版说辞是选对 MPI、选对编译器、再在源码里开 cfdem 耦合开关。但 CFDEM 不是 LAMMPS 的一个内置包它是一套独立编译的 OpenFOAM 库组运行时才把 LIGGGHTS 的共享库链接进来。这决定了安装顺序必须严格先编译 LIGGGHTS再编译 CFDEMcoupling最后用自带算例验证。3.1 版本匹配的三角关系CFDEM 的代码对 OpenFOAM 内部类名的变动很敏感不同分支对应不同 OpenFOAM 版本。公开仓库的分支命名规则一般是CFDEM-版本号-OpenFOAM版本号比如针对 OpenFOAM 8 或 ESI 版的分支。在开始编译之前先确认三件事编成一张表组件推荐取值说明编译器gcc 8.3 及以上老版本对 C17 支持不完整wmake 会报模板错误MPIOpenMPI 4.xMPICH 也可以但和 OpenFOAM 的 compat 层偶有冲突OpenFOAM 版本Foundation 版或 ESI 版按你拿到的 CFDEM 分支决定不可混用# 进入 OpenFOAM 环境以 ESI 版为例 source /opt/openfoam8/etc/bashrc # 设置 CFDEM 相关环境变量 export CFDEM_VERSION6.0 export CFDEM_DIR$HOME/CFDEM/CFDEMcoupling-PUBLIC export LIGGGHTS_DIR$HOME/CFDEM/LIGGGHTS-PUBLICCFDEM_VERSION这个环境变量会被编译脚本用来生成库文件的版本后缀。如果你用的是某个较早的 release 分支版本号要写成分支名里的数字。判断标准很简单进入 CFDEMcoupling 源码目录看src/Makefile里引用$(CFDEM_VERSION)的地方对照写即可。3.2 编译 LIGGGHTS 并确认 cfdem 命令组存在LIGGGHTS-PUBLIC 源码根目录下就是标准的 LAMMPS 风格 src 目录编译目标名为auto它会自动探测 MPI 并生成带共享库支持的二进制cd $LIGGGHTS_DIR/src make -j4 auto # 安装到 LIGGGHTS 的默认 bin 目录 sudo make installauto目标的关键在于它默认开启了-fPIC编译选项CFDEM 的库在运行时需要以动态方式加载 LIGGGHTS 的符号表。如果用make serial这类目标编译后续链接会报 undefined reference 错误。编译结束后检查生成的二进制里是否包含耦合命令$LIGGGHTS_DIR/src/lmp_auto -h | grep cfdCoupling输出里能看到fix cfdCoupling这条命令才算通过。这一步很多人会漏掉装上的是一个普通 LAMMPS后面 CFDEM 一启动就提示无法识别命令。LIGGGHTS 从 LAMMPS 继承了大量分子动力学命令所以-h输出的绝大部分内容看起来和 LAMMPS 一模一样只有多出来的fix cfdCoupling、fix cfdDataExchange这些命令才是耦合真正依赖的扩展。3.3 编译 CFDEMcoupling 与安装验证CFDEMcoupling 的源码用 OpenFOAM 的 wmake 工具链编译不写 Makefile.am而是沿用 OpenFOAM 的Allwmake脚本规范cd $CFDEM_DIR/src ./Allwmake -j4编译过程会生成两类产物。一类是 OpenFOAM 的库文件放进$FOAM_USER_LIBBIN另一类是可执行求解器比如cfdemSolverPiso、cfdemSolverPimple放进$FOAM_USER_APPBIN。编译日志里如果出现wmake error: no such file ... foamVersion.H八成是前三步的source /opt/openfoam8/etc/bashrc没执行干净重新打开一个 shell 按顺序执行。# 验证库是否被正确加载 cfdemSolverPiso -help # 跑一个自带的最小算例做端到端验证 cd $CFDEM_DIR/tutorials ls自带的 tutorials 目录会列出多个算例挑一个名为mpicrudedBed或类似的最小流化床用例跑几分钟。正常情况 stdout 里会出现耦合日志、alpha场的平均值和 LiGggHTs 的每步颗粒信息。出现任何段的FOAM FATAL ERROR先看是出现在流体求解阶段还是颗粒推进阶段——前者多半是 OpenFOAM 版本不匹配后者多半是 LIGGGHTS 编译时没开共享库。4. 写一个最小可复现的流化床耦合算例装通的下一步不应该是立刻改大算例而是把一个最小流化床从零跑起来。这里给出一个 0.01m × 0.01m × 0.02m 的二维等效盒式算例颗粒直径 0.5 mm数量 1000 个气体从底部进入。这个规模在 4 核机器上十分钟内能跑完一个流化周期适合做耦合配置的调试平台。4.1 网格粗细与颗粒直径的比例约束CFDEM 对网格有一个硬性经验约束网格单元的特征尺寸要大于颗粒直径的 3 倍最好在 3 到 5 倍之间。原因是alpha场的计算需要对每个网格单元做颗粒体积的几何叠加颗粒太大、单元太小会让某个单元同时被好几个颗粒的边界跨过导致体积分数出现锯齿状分布进而让曳力振荡。# constant/polyMesh/blockMeshDict 的关键块 vertices ( (0 0 0) (0.01 0 0) (0.01 0.01 0) (0 0.01 0) (0 0 0.02) (0.01 0 0.02) (0.01 0.01 0.02) (0 0.01 0.02) ); blocks ( hex (0 1 2 3 4 5 6 7) (5 5 20) simpleGrading (1 1 1) );这里 5 × 5 × 20 的网格把每个方向单元尺寸拉到 2 mm是颗粒直径的 4 倍处于推荐区间。网格轴向 (z 向) 分层 20 层负责捕捉气泡和床层膨胀水平方向 5 × 5 的粗网格说明这个算例只追求定性正确和流程跑通不是定量研究定量研究需要把这三个方向的网格同时加密。4.2 双份求解器配置couplingProperties 与 in.liggghtsCFDEM 的核心配置在constant/couplingProperties文件里它用 OpenFOAM 的字典格式定义了耦合选项。最关键的三个条目是曳力模型、插值方案和颗粒数据交换方式particleCloud { couplingProperties { solverModel HI; dragModel Gidaspow; turbulenceModel off; regularization off; velocityInterpolationType cellCenter; } }solverModel HI表示采用 Hirt 和 Ishii 的隐式耦合方法把曳力对颗粒速度的依赖隐式处理到流体动量方程里这对密相颗粒流是必要的。velocityInterpolationType cellCenter使用单元中心处的流场值直接赋给颗粒计算最快但会在网格较粗的空间产生阶梯效应更平滑的选项是cellCenter对应的梯度插值调试初期不必追求。颗粒侧由constant/DEM目录下的in.liggghts文件描述它的运行方式是由 CFDEM 在每次耦合循环内调用 LIGGGHTS 推进颗粒子步atom_style granular boundary m m m newton off neighbor 0.003 bin neigh_modify delay 0 every 1 check yes region box block 0 0.01 0 0.01 0 0.02 units box create_box 1 box create_atoms 1 random 1000 12345 box units box set group all density 2500 set group all diameter 0.0005 pair_style gran hertz/history pair_coeff * * 1e7 0.45 0.5 0.5 timestep 1e-6 fix gravity all gravity 9.81 vector 0 0 -1 # 先自由下落形成堆积床 run 5000 # 开启与 CFDEM 的双向耦合 fix cfd all cfdCoupling run 50000颗粒下落 5000 步即 5ms 物理时间在 0.01m 高的盒子里足够重力和接触力取得平衡形成稳定的初始床层。随后fix cfd all cfdCoupling把 LIGGGHTS 挂到 CFDEM 的调度循环上这时开始跑的 50000 步由 CFDEM 在每个耦合步里以子步方式驱动。注意timestep容不得粗心1e-6 秒对于 0.5mm 颗粒配合 1e7 Pa 的剪切模量是安全的如果把剪切模量改成 1e8 Pa瑞利波时间步会缩短一个量级1e-6 就不再收敛表现为颗粒互相穿透或能量爆炸。4.3 启动命令与运行时日志的判断顺序blockMesh -case . # 准备耦合子目录CFDEM 会在运行时创建颗粒输出 mkdir -p postProcessing/dem # 启动耦合求解器 cfdemSolverPiso -case . log.cfd 21 cfdemSolverPiso是 CFDEM 提供的压力速度耦合求解器名字里的 PISO 表明它用 PISO 算法推进流体。启动后肉眼判断运行状态要按三层顺序看第一层看log.cfd最后几行有没有FOAM FATAL ERROR第二层看进程是否还活着CFDEM 在 LIGGGHTS 无输出时会显得像死锁实际是颗粒子步在跑第三层看固定的时间间隔有没有生成postProcessing下的场输出。用tail -f log.cfd观察每步的alpha平均值如果数值一开始接近 0.061000 个颗粒体积除以盒体总容积,随后随着床层膨胀逐渐下降说明耦合的空间插值链路是通的。5. 把耦合调稳的三个参数与三个常见坑耦合参数之间有强烈的联动关系单独调一个很难奏效。核心是先确定 DEM 子步与 CFD 时间步的比例再确定每个耦合间隔内 CFD 走几步最后用松弛因子把两个求解器的相互反馈压住。参数推荐取值调参信号DEM 时间步1e-6 s0.5mm 颗粒颗粒穿透边长或总动能发散CFD 时间步1e-4 s耦合步间隔 10库朗数超过 1 就减半曳力松弛因子0.6 ~ 0.8压力场逐周期振荡就往下调网格/粒径比3 ~ 5alpha 场呈棋盘分布就减小颗粒半径或加密网格第一个多发的坑是颗粒穿透量现象是颗粒在重力作用下漏出计算域底部。原因通常不是接触模型参数问题而是boundary m m m三个方向都用了周期性边界盒子底部没有物理墙面。流化床的底面必须是固定壁面在 in.liggghts 里要显式添加几何墙fix wallbase all wall/gran model hertz tangential history primitive type 1 zplane 0这行fix在 z0 处定义了一个无限平面墙颗粒与它的接触力由与颗粒对颗粒接触完全相同的 Hertz-Mindlin 模型计算。忘了加这行颗粒就会像穿过空气一样漏走。第二个坑是压力场出现周期性振荡但速度残差却收敛得不错。这通常发生在把 CFD 时间步从 1e-5 直接加到 5e-4 之后CFD 步长过大一个耦合步内曳力源项的变化幅度超过了流体压力修正能处理的量级PISO 的外循环数已经压不住。CFDEM 的couplingProperties里对源项做松弛处理couplingProperties { underRelaxation 0.7; momentumTransferModel kTGF; }underRelaxation取 0.7 意味着每次耦合循环只把新的动量源项和上一步的值做七三开切断振荡的正反馈。如果振荡仍然明显先把 CFD 时间步降低两倍再调松弛顺序不可颠倒。第三个坑不是发散而是假收敛计算结果看起来很正常但床层膨胀高度明显低于经验值。一个常见诱因是velocityInterpolationType用了默认的cellCenter把包含颗粒的单元中心速度直接赋给了颗粒。当网格偏粗时这个插值会低估颗粒位置处的真实流速等效于曳力打了折扣。把插值方式换成cellCenter对应的梯度插值选项或在阻力模型上改成DiFelice做交叉验证两者的差值能反映出插值误差的边界。验证耦合是否正确比看收敛残差更有效的是做力平衡核对在log.cfd里打印每个耦合步的全场曳力对颗粒的作用力与相同入口条件下纯压降乘以床层截面积得到的总压差力比较。一个稳态流化床的压降应当约等于床层颗粒的总浮重扣除浮力误差超过 15% 时优先检查alpha场而非曳力模型。最终把dragModel依序从 Gidaspow 换成 WenYu 再换成 DiFelice观察同一边界条件下的压降和床层空隙率曲线差异就能定量看到模型假设对耦合结果的影响边界。本文还有配套的精品资源点击获取