《Advanced Powder Technology》：Fully resolved CFD–DEM simulation with three-phase contact line dynamics on irregularly shaped particles
编辑推荐：
Kimiaki Washino|Ei L. Chan|Takuya Tsuji大阪大学机械工程系，日本大阪府吹田市565-0871摘要含有不规则固体颗粒及两种不相溶液体的三相流在许多自然和工业过程中十分常见，因此人们对可靠的数值模拟模型的需求日益增加。然而，现有的大多数模型都存在
Kimiaki Washino|Ei L. Chan|Takuya Tsuji
大阪大学机械工程系，日本大阪府吹田市565-0871
摘要
含有不规则固体颗粒及两种不相溶液体的三相流在许多自然和工业过程中十分常见，因此人们对可靠的数值模拟模型的需求日益增加。然而，现有的大多数模型都存在简化问题，比如假设颗粒为球形，或未考虑三相接触线动态。在本研究中，作者开发了一种全新的全解析CFD–DEM耦合模型，将多球体模型引入到之前提出的先进三相流模型中，用以描述任意形状的颗粒。通过多项验证研究，该模型在雷诺数介于0.57至450之间的流体-固体两相系统中，以及立方体状颗粒间的三相系统中，均展现了对流体动力学相互作用和毛细效应的准确描述能力。随后，该模型被用于模拟二相和三相系统中的岩石颗粒，证明了其能够有效捕捉高度不规则且不对称的颗粒形状以及固体润湿性对流动动力学的影响。
引言
在三相流中，固体颗粒与两种不相溶液体共存、相互流动并发生作用，这类现象广泛存在于各种自然和工业过程中。例如河流中的泥沙输送、化学反应器以及浆料混合/分离过程等。在这些过程中，颗粒往往具有非球形结构，而这会显著影响流动动力学。比如，不规则形状的颗粒容易发生互锁现象[1]、[2]，进而产生极其复杂的流体动力学相互作用[3]、[4]，其表现与球形颗粒有本质不同。在三相流中，各相交汇处的接触线动态也起着至关重要的作用。从物理层面来看，接触线行为是由两种（不相溶的）液体与固体基底之间的分子相互作用决定的。从宏观角度而言，这种分子相互作用表现为固体表面的“润湿性”，它能够对颗粒施加毛细力，促使液体自发渗入颗粒间的孔隙中（即毛细作用）。为了更好地控制与优化相关工艺，了解非球形特征以及接触线动态对三相流的影响至关重要。不过，这类复杂现象往往难以通过实验观测，因此人们对可靠数值模拟模型的需求愈发强烈。
用于模拟颗粒流的最常用方法之一是离散元法（DEM）[5]，该方法通过计算作用在单个颗粒上的力与扭矩，并依据牛顿第二定律追踪其运动轨迹。DEM能够提供颗粒层面的详细信息，且与连续介质模型不同，无需复杂的本构关系。传统上，由于接触检测算法简单高效[6]、[7]，颗粒通常被建模为理想球体。如今，DEM已有多种形状模型，包括多球体（团块）模型[8]、[9]、[10]、[11]，多面体模型[12]、[13]、[14]，超二次曲线模型[15]、[16]，以及有符号距离模型[17]、[18]、[19]。每种模型都有各自的优缺点，相关的综合评述可见于[20]、[21]。
根据[20]的论述，最常用的形状模型是多球体模型，该模型通过将多个子球体以固定相对位置聚集在一起，来近似表示非球形颗粒的形状。这种方法优势明显：可以利用成熟的球体接触检测算法进行接触判断，而且通过调整子球体的数量、大小和位置，即可轻松表示任何形状的颗粒。但其主要缺点是计算成本会随着子球体数量的增加而迅速上升。另一种受到关注的模型是超二次曲线模型，该模型通过超二次曲线函数隐式描述颗粒形状[22]。只需调整5个参数即可描述椭球体、长方体、圆盘和棒状等多种形状，由于理论上可以获取颗粒表面法向量，因此可能比其他形状模型具有更高的计算效率。但该模型的局限性在于只能描述旋转对称形状，无法用于模拟不规则形状的颗粒。
在模拟三相流时，除了考虑颗粒间的相互作用外，还需考虑颗粒与周围流体之间的相互作用。为此，DEM通常与计算流体动力学（CFD）相结合。CFD–DEM耦合模型大致可分为两类：未解析模型和解析模型。未解析模型采用局部平均的流体控制方程[23]、[24]、[25]、[26]、[27]、[28]、[29]，适用于模拟大规模系统；而解析模型[30]、[31]、[32]、[33]、[34]、[35]、[36]、[37]则能精确求解单个颗粒周围的流场，有助于深入理解颗粒与流体之间的流体动力学及毛细相互作用。本研究采用的是解析模型。
解析型CFD–DEM耦合模型通常采用“单一流体”建模思路，在整个计算域内求解同一组控制方程，不同流体则通过界面识别方法（如体积流体法[38]、水平集法[39]、受限插值剖面法[40]以及相场法[41]）或界面追踪方法（包括前沿追踪法[42]）来区分其不同属性。在三相交界处需要特殊处理，以便准确描述接触线动态：对于界面识别方法，可运用扩展的指示函数[33]、[37]；而对于界面追踪方法，则可使用滑移速度[43]。现有的三相流模型主要分为以下两类：
- •
无接触线动态的非球形颗粒模型[44]、[45]、[46]
- •
具有接触线动态的球形颗粒模型[47]、[48]、[49]、[50]
据我们所知，目前唯一能够描述具有接触线动态的非球形颗粒的模型，是作者在之前的研究中开发的模型[51]。不过该模型采用超二次曲线模型来描述颗粒形状，因此仅适用于规则且具有旋转对称性的颗粒。
在本研究中，作者对之前的模型进行了进一步改进，引入多球体模型，用以模拟具有接触线动态的不规则形状颗粒的三相流。研究还通过多项验证研究，证明了该模型在描述流体动力学力和毛细力方面的准确性。此外，还通过一些示范性虚拟实验，分析了颗粒形状和固体润湿性对流动动力学的影响。
章节节选
问题设定
本研究所分析的初始三相系统如图1a所示，其中不规则形状的固体颗粒分散在两种流体（流体0和流体1）中。这两种流体被假设为牛顿流体、不可压缩、不相溶且温度恒定，流体间的尖锐边界处存在表面张力。颗粒与流体可以在彼此的尖锐边界处相互作用并同时移动。颗粒表面的润湿性由接触情况决定
流体模型描述
用于描述流体运动的模型大多基于作者之前的研究成果[51]，本节将对这些模型进行简要介绍。
固相模型描述
在多球体模型中，非球形颗粒的形状被近似为多个子球体的集合[8]。原则上，通过调整子球体的数量、大小和位置，即可表示任意形状的颗粒。如有需要，还可以让子球体相互重叠，以提高颗粒表面的平滑度。子球体之间的相对位置是固定的。在每次模拟开始时，会通过蒙特卡洛方法生成大量随机点
实现与求解流程
前文所述的模型是在开源软件包中实现的：CFD计算采用OpenFOAM[60]，DEM计算则使用LIGGGHTS[61]。这两个软件包都完全支持消息传递接口（MPI），便于实现高效的分布式内存并行计算。本研究中的所有模拟都是在九州大学Genkai超级计算机的A子系统中运行的，每个计算节点都配备了双Intel Xeon Platinum 8490H处理器
验证
本节介绍了四项验证模拟，涵盖内容如下：
- •
不同雷诺数条件下的流体动力学相互作用（第6.1节：简单剪切流中的椭球体——蠕动流动；第6.2节：单个球形颗粒的沉降；第6.3节：三球体团块的沉降）
- •
毛细相互作用（第6.4节）
此外还进行了网格敏感性分析，确认本节中的所有模拟结果均不受网格质量的影响。
虚拟实验
本节通过图20所示的由50个相互重叠的子球体构成的岩石颗粒，开展了两项示范性虚拟实验。该颗粒形状极为不规则且不对称，多球体模板可在LIGGGHTS-3.8公版的一个教学案例中找到。两项实验所使用的模板及颗粒参数分别列于表4和表5中。由于该颗粒的下半部分比上半部分更厚（见图20的侧视图），因此
结论
本研究开发了一种全新的全解析CFD–DEM耦合模型，可用于模拟含有不规则形状颗粒以及三相接触线动态的三相流。该模型弥补了现有数值模拟模型的不足，因为现有模型往往存在简化问题，比如假设颗粒为球形或忽略接触线动态。新模型将多球体模型整合到作者之前的框架[51]中，从而能够表征任意形状的颗粒
CRediT作者贡献说明
Ei L. Chan：撰写——审稿与编辑、可视化、验证、研究分析、数据整理。Kimiaki Washino：撰写——初稿撰写、可视化、验证、软件开发、资源准备、项目管理、方法设计、研究分析、资金筹集、数据整理、概念构思。Takuya Tsuji：撰写——审稿与编辑、监督指导。
利益冲突声明
作者声明，他们不存在任何可能影响本文研究结果的已知财务利益或个人关系。
致谢
作者感谢Hosokawa Powder Technology Foundation为本研究提供的财政支持。此外，本研究还得到了九州大学信息技术研究所提供的计算资源帮助。
