十亿颗粒算不动？粗粒化模型让流化床DEM仿真走进工业级
粗粒化模型（Coarse Grain Model, CGM）提供了一条可行的路径。本文介绍 Sakai 等人提出的 CGM，重点阐述该模型如何从非黏性颗粒扩展到黏性细颗粒（考虑范德华力），并通过二维流化床数值实验验证其有效性。
一、为什么需要粗粒化？
离散元法（DEM）是一种拉格朗日方法，对每个颗粒单独求解牛顿第二定律，天然适合研究颗粒尺度的流动特征。自 Tsuji 等人将 DEM 与 CFD 结合以来，该方法已被广泛应用于喷动床、鼓泡流化床、循环流化床、气力输送、堆积、湿法球磨等气固流动问题中。
然而，DEM 存在一个关键瓶颈：可计算的颗粒数量受限于过高的计算成本。即便使用最新的 PC 或 PC 集群，已有研究中可处理的颗粒数量最多也只有几十万。因此，现有 DEM 难以直接应用于工业规模系统。
粗粒化模型的核心思想是：用一个粗粒粒子代表一组原始粒子。通过严格的尺度推导，使粗粒粒子的运动学与动力学行为在统计意义上等价于原始粒子群。这样，粒子总数可减少为原来的 1/l³（l 为粗粒化比例），而计算时间步长却无需改变。此前，CGM 已成功应用于气力输送和鼓泡流化床中的非黏性颗粒模拟。本文则将其进一步拓展至黏性细颗粒。
二、粗粒化模型的力学框架
2.1 平动与转动的等价性
如图1所示，一个边长为 l 的粗粒化立方体中包含 l³ 个原始颗粒。粗粒粒子的平动速度取为原始粒子群的平均速度；转动则假设原始粒子围绕各自质心以相同角速度旋转。
粗粒粒子的角速度满足以下关系，以确保转动动能与原始粒子群一致：
2.2 接触力与曳力的缩放
当两个粗粒粒子发生碰撞时，模型假设内部所有原始粒子同时发生二元碰撞。基于动能等价原则，接触力（法向与切向）和曳力均按 l³ 缩放。法向接触力表示为：
切向接触力则根据是否滑动，分别采用弹性恢复或库仑摩擦形式。曳力表达式为：
其中下标 CGM 表示粗粒模型，O 表示原始颗粒。这一处理确保了粗粒粒子在碰撞和流体作用下的响应与原始粒子群一致。
2.3 范德华力的关键缩放：l²
本文的核心创新在于将 CGM 扩展到黏性颗粒。作者从势能等价出发推导范德华力的缩放关系。如图2所示，当两个粗粒粒子靠近时，假设内部原始粒子均匀分布且间距同步变化，粗粒粒子表面间距 hCGM 与原始粒子表面间距 hO 满足：
通过对势能积分，得到粗粒粒子间范德华力的表达式：
= l³ · HA · ds,CGM / (6 · hCGM²)
注意，范德华力按 l² 缩放，而非接触力和曳力的 l³。这一差异源于范德华力是表面力，其衰减特性与体积力不同。为避免数值奇点，粒子表面间距设置了截断值（0.4 nm）。
2.4 完整的运动方程
综合曳力、压力梯度力、接触力、重力和范德华力，粗粒粒子的运动方程可写为：
该模型在二维和三维模拟中均可使用。即使在二维模拟中，颗粒形状仍按球形（三维物体）处理，因此方程保持一致。
2.5 时间步长不变
由于质量和刚度均按 l³ 缩放，粗粒模型的时间步长与原始系统完全相同：
这意味着粗粒化带来的加速完全来自粒子数量的减少，而非时间步长的放大。
三、流化床验证方案
为验证模型，作者在二维流化床中进行了 DEM-CFD 耦合模拟。如图3所示，计算域为 30 mm × 360 mm 的矩形区域，颗粒初始随机堆积。底部均匀进气，表观气速从 0.07 m/s 开始阶梯式递减，如图4所示。
三组工况分别对应不同的黏性水平：
每组工况包含 4 个算例：原始系统（l=1，90,000 颗粒）、CGM l=2（22,500 颗粒）、CGM l=3（10,000 颗粒），以及一个对照算例——直接使用 600 μm 大颗粒但不采用粗粒化缩放（10,000 颗粒，即 Case X-4）。
四、结果与讨论
4.1 非黏性颗粒（Case 1）
图5展示了 Case 1 在 5 s、10 s、15 s 和 20 s 时刻的典型快照。随着表观气速降低，床层高度逐渐下降，气泡在上升过程中不断长大，且气速越高气泡越大。Case 1-1（原始系统）、Case 1-2（l=2）和 Case 1-3（l=3）的流态化形态高度一致。
然而，对照算例 Case 1-4（无粗粒化的大颗粒）表现完全不同：如图6所示，颗粒几乎停留在初始堆积位置，几乎没有流化。这证明，简单地放大粒径而不做力缩放，无法复现细颗粒行为。
定量对比显示，床层高度随气速的变化曲线（图7）和压降随气速的变化曲线（图8）几乎重合。压降在最小流化速度附近开始下降，估计最小流化速度约为 0.020 m/s。
4.2 弱黏性颗粒（Case 2）
图9展示了 Case 2 的典型快照。引入弱范德华力后，床层整体仍表现出良好的流态化特征，少量颗粒黏附在壁面上，气泡行为与 Case 1 类似。CGM 各算例的流态化形态与原始系统定性一致。
定量上，图10的床高-气速曲线和图11的压降-气速曲线再次验证了 CGM 的准确性，Case 2-1 至 2-3 的结果几乎等价，最小流化速度保持一致。同样，无粗粒化的大颗粒算例 Case 2-4 几乎不发生运动。
4.3 强黏性颗粒（Case 3）
当 Hamaker 常数增大到 1.0×10⁻¹⁹ J 时，颗粒间黏性显著增强。如图12所示，流化床内出现明显的沟流（channeling）现象，大量颗粒黏附在壁面。由于颗粒流动具有随机性，沟流的具体通道在原始系统与 CGM 之间未必同步，但宏观流动形态（沟流存在、壁面黏附程度、床层膨胀行为）高度一致。
定量上，图13的床高-气速曲线和图14的压降-气速曲线再次验证了 CGM 的准确性，Case 3-1 至 3-3 的结果几乎等价，最小流化速度也得到了精确复现。这充分说明，本文提出的范德华力 l² 缩放律是正确的。
五、结论与工程意义
本文将粗粒化模型从非黏性颗粒推广到黏性细颗粒，核心在于推导了范德华力的 l² 缩放关系。通过二维流化床数值实验，得到以下关键结论：
对于工程实践而言，粗粒化模型意味着：在保持物理可信度的前提下，Geldart A 类甚至 C 类细颗粒的大规模流化床仿真，终于有望从实验室规模走向工业规模。粗粒化比例的选择受限于 CFD 网格尺寸（计算颗粒直径应远小于网格尺寸），但在该约束内，CGM 为十亿级颗粒系统的数值模拟打开了大门。
