DEM-CFD耦合

没错。DEM（离散元素法：Discrete Element Method）追踪各个粒子的运动，CFD求解流体场。通过双向耦合两者，再现粒子-流体间的相互作用。

决定性的区别在于粒子间接触的处理。DPM中粒子间碰撞采用简化模型（概率性碰撞）或忽略，而DEM用弹性弹簧-阻尼器模型严格计算粒子间的接触力。因此，对于粉体和颗粒这样粒子间力重要的系统更合适。

追踪各粒子 $i$ 的平移运动和旋转运动。

$\mathbf{F}_{c,ij}$ 是与粒子 $j$ 的接触力，Hertz-Mindlin模型是代表。法向接触力表示为：

其中 $E^*$ 是等效杨氏模量，$R^*$ 是等效半径，$\delta_n$ 是重叠量，$\gamma_n$ 是阻尼系数。

求解局部平均化的Navier-Stokes方程。考虑粒子存在导致的孔隙率 $\varepsilon_f$ 。

$\mathbf{S}_p$ 是粒子对流体的反力（动量源项），即CFD网格单元内所有粒子所受流体力之和除以体积。

最重要的是阻力，采用依赖局部孔隙率的阻力模型。

压力梯度力、虚拟质量力、Saffman升力、Magnus力等也可以考虑，但对于密度比 $\rho_p / \rho_f \gg 1$（粉体-空气系等），阻力是支配性的，其他力在多数情况下可以省略。

我来说明基本耦合方案。

1. 在CFD时间步 $\Delta t_{CFD}$ 内求解流体场

DEM时间步由于接触力计算需要非常小，Rayleigh时间的20～30%是目安。

典型情况下，CFD的1步对应DEM的100～1000个子步骤。这是DEM-CFD计算成本的主要原因。

从CFD单元内的粒子体积计算。对于跨越单元边界的粒子，有以下处理方法。

CFD单元大小必须是粒子径的3～5倍以上，这是unresolved DEM-CFD的前提条件。单元小于粒子时孔隙率定义会崩溃。

unresolved是粒子径小于CFD网格，用阻力模型表现相互作用。resolved是粒子径大于网格，直接解析粒子表面的流场。通过浸入边界法或Overset Mesh实现，但粒子数限制在数百个左右。

以流动床反应器为例来说明。

1. 粒子物性确定：粒径分布、密度、杨氏模量、反弹系数、摩擦系数

Hertz-Mindlin模型的参数设置很重要。

使用实际材料杨氏模量（玻璃：70 GPa）会导致DEM时间步极小，计算变得不现实。降至$10^7$～$10^8$ Pa，粒子流动行为在很多研究中基本不变。但反弹系数和接触时间会改变，所以需要标定。

最常见的是安息角（angle of repose）的标定。实验中从圆筒排出粒子形成的堆积角度，用DEM摩擦系数调整使其相符。

仅从安息角无法唯一确定参数，需要结合排出流量、圆筒转动试验等多个实验来验证。

大致分为原生实现和外部耦合两种。

原生实现（Fluent DEM、STAR-CCM+ DEM）配置简单，数据交换开销少。但DEM功能成熟度不如EDEM或Rocky DEM。

EDEM粒子形状自由度高（多面体、纤维等），Rocky DEM大规模GPU计算强。粒子数超过百万时，Rocky DEM的GPU实现速度压倒优势。

LIGGGHTS + OpenFOAM组合在学术界应用最广。定制自由度高，易于实现新的阻力模型和耦合方案。但没有GUI，全部用脚本配置。

我介绍几个研究方向。

为了减少粒子数，粗粒化法受关注。用1个"代表粒子"表示多个实粒子。代表粒子径为实粒子的$k$倍时，粒子数减为$k^3$分之一。

需要适当的缩放律。Sakai等人（2009）的粗粒化模型中，对阻力进行$k^3$倍的修正。有Exact Scaling法、Statistical Scaling法等多个流派。

实际粉体不是球形。非球形粒子处理是重要研究课题。

EDEM标准采用Multi-sphere法，用重叠球的集合近似复杂形状。Rocky DEM支持多面体法，可实现更精确的形状表示。

用DNS或LBM（格子Boltzmann法）计算的详细流体力数据作为教师数据，用神经网络构建阻力模型的研究在进行。Beetstra等人、Tang等人的模型是代表。

特别是在高孔隙率过渡区和非球形粒子阻力上，精度大幅高于传统相关式。但DNS数据生成本身计算成本很高是课题。

逐个看。

症状：粒子穿过壁面或跑出CFD域。

对应：

对应：

原因和对应：

症状：粒子僵化不动。

对应：