CFD-DEM模拟简介

CFD-DEM  

CFD-DEM模拟是一种基于欧拉-拉格朗日参考系的方法，是离散模拟的典型代表(Ge et al., 2017)。该方法的特征是：

1 将流体视为连续介质，对计算空间进行网格划分，然后求解纳维斯托克斯方程得到每个网格上流体的速度、密度、和压力等；

2 对颗粒相则接计算其受力，然后根据牛顿定律计算其加速度、速度和位置的变化。

流体相的求解一般采用有限体积等传统计算流体力学方法，相对成熟，并不是制约该类方法的难点。颗粒相的求解，虽然在算法上较为简单，但巨大的颗粒数目导致计算量巨大，往往限制了该类方法的应用。如何准确计算流体体积分数，如何计算曳力，如何在非规则网格中计算颗粒所在网格等两相耦合计算也是该方法的难点。

本文简要介绍CFD-DEM模型方程、算法流程、研究热点以及典型应用。

1 研究模型

目前在CFD-DEM学术论文中最常见的控制方程都基于A类模型，即颗粒受力直接包含了压力梯度力，而B类模型只考虑了颗粒排开流体体积所带来的阿基米德浮力。这导致B类模型无法正确的模拟旋风分离器等旋转流场(Zhouet al., 2010)。在A类模型中，Ishii形式的动量方程应用较多，其具体形式如下：

• 方程1

• 方程2

其中为流体体积分数，为其密度，u_f 为其速度，为应力，f_drag表示流体对颗粒的曳力。应力的计算方式如下：

• 方程3
  

其中和分别表示动力粘度和第二粘性系数。曳力的计算方式如下：

• 方程4

该公式的含义为，先插值计算颗粒质心处的流体速度，然后计算颗粒受到的曳力，最后根据牛顿第三定律插值计算CFD网格内颗粒对流体的反作用力。

颗粒的运动包括平动和转动，其平动控制方程如下：

• 方程5
  

公式右侧各项分别为重力、压力梯度力、曳力和碰撞力。

颗粒转动的控制方程如下：

• 方程6
  

其中，I_p为转动惯量，其数值为 m_pd_p²/10. W为转动速度，T为扭矩。

注：碰撞力和扭矩的具体计算方法，可以参考DEM专业书籍或文献。

2 算法流程

CFD-DEM的计算流程主要包括：初始化、流体相计算和颗粒相计算。

初始化

初始化操作主要包括读入计算参数，申请内存空间，初始化流场和颗粒，构建颗粒邻居列表，构建颗粒流体网格映射关系。在初始化颗粒位置时，要保证颗粒-颗粒以及颗粒-壁面无重叠。

流体相计算

流体相的计算采用压力耦合方程组的半隐式方法（SIMPLE）即根据假定的速度场与压力场对动量方程离散求解，然后求解压力修正方程，对压力和速度进行修正，迭代直到流场收敛。该计算流程和单相流的计算一致，流体网格中的流体体积分数和平均颗粒速度由颗粒位置信息显式插值计算得到。在进行流体计算时，颗粒的速度、位置等保持不变。

颗粒相计算

当流场收敛后，进行颗粒相的计算，主要包括以下8个子步骤：

1 颗粒-颗粒、颗粒-壁面碰撞力计算；

2 计算流体压力梯度力和曳力；

3 更新颗粒速度和位置；

4 多进程并行计算时，传递进程边界颗粒信息；

5 更新颗粒邻居列表；

6 更新颗粒-流体网格映射关系；

7 更新颗粒-流体网格插值系数；

8 更新流体网格中流体体积分数。

一般情况下，由于流体的计算时间步长大于颗粒计算时间步长，以上1-8步需要重复多次，但在每一步的计算中，流体的速度、密度和压力分布保持不变。当颗粒和流体的时间同步后，停止颗粒的计算转为流体的计算。重复以上步骤，直到达到预先设定的模拟时间。

3  研究热点

CFD-DEM研究热点集中在提高其计算速度和精度。

CFD计算速度的提高主要依赖于高效的线性方程组求解器和合适的流体方程求解算法。

	目前常用的线性方程求解器为代数多重网格和共轭梯度等迭代方法，以减少求解离散方程的时间。
	流体方程求解算法主要包括SIMPLE和PISO。PISO算法一般采用较大的松弛因子或不进行松弛，因此其压力-速度耦合迭代修正次数小于SIMPLE算法中的迭代次数，该部分的研究一般可直接采用单相流领域的进展，由于CFD计算在CFD-DEM中所占比重较小，多相流研究人员很少直接涉及该领域。

DEM计算速度的提高主要依赖于高效的碰撞搜索算法。

目前应用最广泛的方法是邻居列表法：将颗粒映射到搜索网格中，计算颗粒与相邻网格中其他颗粒的距离并建立邻居列表。在碰撞检测时只需检测其邻居列表中的颗粒。邻居颗粒的截断距离越大，则邻居颗粒数目越多，可以减少邻居列表更新的频率；反之则需要频繁更新。对于多粒径系统，采用均匀分布的单一搜索网格时，由于搜索网格需大于系统中最大的颗粒，这导致网格内颗粒数目巨大，从而降低了搜索效率。因此，多采用树形或多层嵌套搜索网格。除了改进算法外，还可以采用GPU加速(Xu et al., 2011)或者粗粒化方法(Lu et al., 2016)直接提高DEM的计算速度。

CFD-DEM模拟的计算结果受CFD网格质量、离散格式、颗粒碰撞参数、时间步长、曳力模型和耦合方法等多种因素的影响。目前研究热点集中在曳力模型、耦合方法以及颗粒形状。

	CFD-DEM模拟可以采用DNS曳力、EMMS曳力或基于试验数据的曳力。相比于双流体模拟，CFD-DEM模拟可以获得每个网格内颗粒的分布信息，因此该信息也可被用来对曳力进行修正。
	耦合方法包括颗粒信息向CFD网格映射和流场信息向颗粒映射。前者用来计算CFD网格内流体体积分数和平均颗粒速度，后者用来计算颗粒质心处的流速和压力梯度。对于笛卡尔网格，可以直接计算颗粒所在的CFD网格，并采用线性插值的方法计算。对于其他类型的CFD网格，需要根据CFD网格的拓扑结构和颗粒运动速度判定颗粒所在的网格，然后在进行插值计算。此外，还可以在原有CFD网格的基础上重建一套笛卡尔网格，并在这套网格上进行两相的耦合。
	大多数CFD-DEM模拟都将颗粒假定为球形，但实际中常常会有圆柱形、胶囊形、方形等规则形状以及各种非规则形状。非球形颗粒直接的碰撞、与流体的作用计算都是研究的热点(Zhong et al., 2016)。

4  典型应用

CFD-DEM常常用来模拟气固或液固多相流。由于直接跟踪颗粒的运动，该方法可以方便的用来研究：

颗粒在反应器中的停留时间分布(Lu et al., 2017)；

多粒径系统的混合和离析(图1，Lu et al., 2018)；

颗粒之间的静电力(Korevaar et al., 2014)、粘性力(Zhou et al., 2018)。

由于直接计算颗粒-颗粒，颗粒-壁面的碰撞力，该方法还可以研究颗粒和反应器的磨损(图2，Christopher et al., 2015)。

采用并行计算以及粗粒化模型提高计算速度后，该方法还可以用来研究煤气化等(图3，Hu et al., 2018)。

图1. 不同粒径的FCC 颗粒在流化床中的分布，在表观气速为0.05m/s的鼓泡床中模拟0.2s后发现粒径较大的FCC颗粒分布在底部，较小的颗粒分布在床层上部。图片摘自Elsevier (Lu et al., 2018)

图2. 由颗粒- 壁面碰撞引起的弯管壁面磨损。图片摘自Elsevier (Christopher et al., 2015)

图3. 流化床煤气化反应器中煤炭颗粒（红色）和沙子（蓝色）的分布。图片摘自Elsevier(Hu et al., 2018)

本文是面向初步接触CFD-DEM的读者所写的简介类文章，并不是完整的模型综述，因此未包含许多该领域的重要文献，建议对该模型有兴趣的读者自行查阅相关文献，以便对该方法有更为完整的认识，后续会向大家介绍详细的CFD-DEM模拟。

作者： 顾龙夏  

来源：多相流在线