颗粒流体系统的格子 Boltzmann 数值方法研究进展

2022-07-04王利民付少童

计算力学学报 2022年3期

王利民，付少童

(1.中国科学院过程工程研究所多相复杂系统国家重点实验室，北京 100190；2.中国科学院大学化学工程学院，北京 100049)

1 引言

颗粒流体系统广泛存在于自然界和工业过程，由于其涉及复杂的多相流动过程和时空多尺度结构，通过实验或理论分析存在较大困难，因此基于数值模拟展开颗粒流体系统的研究已成为不可或缺的重要手段。根据对流体相和颗粒相描述精度的不同，颗粒流体系统的模拟方法通常分为颗粒解析直接数值模拟PR-DNS(particle-resolved direct numerical simulation)、离散颗粒模拟DPS(discrete particle simulation)和双流体模型TFM(two -fluid model)三类[1]。

颗粒流体系统模拟中,TFM的颗粒相处理为连续介质的拟流体，PR-DNS和DPS的颗粒相通常是离散处理，根据牛顿第二定律对颗粒进行 Lagrange 跟踪。目前在颗粒流体系统的PR-DNS中，常忽略或简单处理颗粒之间的作用。然而，大多数体系中颗粒间的相互作用至关重要，需要准确和高效地对其进行描述。颗粒碰撞处理主要有软球模型、硬球模型和直接Monte Carlo方法DSMC(direct simulation Monte Carlo)。软球模型的典型代表是离散单元法DEM(discrete element method)[2]，该方法基于弹性阻尼和滑移等机制计算颗粒间的碰撞力，再根据牛顿第二定律采用时间驱动算法更新颗粒的速度和位置；硬球模型将颗粒碰撞处理成二体碰撞，其不求解颗粒受力，基于动量守恒原理直接获得颗粒碰撞后的速度信息；DSMC[3]同样采用硬球模型更新碰撞后颗粒速度与角速度，但运用概率抽样确定颗粒是否碰撞，因此该方法也可视为硬球模型的一种变化形式。其中，软球模型由于能够获得颗粒详尽的受力信息而得到广泛应用。

颗粒流体系统中的流体相主要基于连续介质假设框架下的N-S(Navier-Stokes)方程或其简化形式，采用常规的有限差分、有限体积和有限元等方法进行求解。由于该类方法流程复杂且包含隐式算法，在并行化和大规模计算上还存在困难[4]。格子Boltzmann方法LBM(lattice Boltzmann method)[5]是20世纪80年代末逐渐发展起来的一种流体系统建模和模拟的新方法。从离散的网格来说，其具有欧拉方法的属性；从离散的粒子观点来看，又具备拉格朗日方法的特点。这种特性使得该方法的研究思路与传统的流体模拟方法完全不同，其属于介于宏观连续介质模型和微观分子动力学模型之间的介观模型，且物理背景清晰，受到国内外学者的广泛关注。随着LBM基础理论不断完善，目前其研究已经从理论向工程应用发展，并在许多复杂物理现象和过程进行了大量的应用[6]，如湍流DNS、大涡模拟和耦合各种湍流模型的模拟等；针对多组分系统和不混溶的多流体系统，分别建立了颜色模型、伪势模型、自由能模型和动力学模型等；针对化学反应系统，建立扩散和燃烧模型。此外，LBM已拓展到模拟非牛顿流体、磁流体、聚合物流动、血液流、光子运动、声场和微流动等复杂体系[6]。

工业流化床反应器的最大流动结构可达到米量级，但这些结构却受发生在毫米量级的颗粒-颗粒间碰撞和颗粒-流体间相互作用的直接影响。目前很难通过某一具体的模拟方法涵盖所有的时间和空间尺度，根据关注重点问题的不同，应用格子Boltzmann框架可在不同时空尺度和精度上开展颗粒流体系统模拟(图1)，如在微观层次，LB -based PR-DNS拟获得的详细信息深入探索颗粒流体界面的流动、传递和反应以及两相相互作用的本构关系；在介观层次，利用 LB -based DPS探索实验室规模的颗粒流体流动，建立更为准确的宏观多相流动模型；在宏观层次，发展格子Boltzmann双流体模型(LB -based TFM)有望实现工业级别的运算。

图1 基于格子Boltzmann的颗粒流体系统计算框架

本文将重点介绍LBM在颗粒流体系统模拟方面的研究进展，包括 LB -based PR-DNS方法、LB -based DPS方法及其最新应用，以及关于开发 LB -based TFM研究的一些进展。

2 颗粒解析直接数值模拟

颗粒流体系统PR-DNS与单相湍流DNS相类似，其中流体相的控制方程直接采用数值计算求解，无任何湍流模型，流场中颗粒周围计算网格缩小到颗粒尺度以下，能够分辨包含Kolmogorov耗散尺度以内的所有空间尺度，而颗粒的受力则通过对其表面的黏性力与压力进行积分获得，不引入经验模型[7]。根据流体相的描述方式不同，颗粒流体系统的PR-DNS又可以分成三类(图2)，一是基于N-S方程的PR-DNS方法；二是基于格子的PR-DNS方法；三是基于粒子的PR-DNS方法。LB -based PR-DNS方法属于基于格子的PR-DNS方法。根据流固作用处理方式不同，基于LBM的颗粒流体系统PR-DNS分为边界链法和格点法。

图2 颗粒流体系统的颗粒解析直接数值模拟

2.1 边界链法

Ladd[8]率先提出了基于有限体积颗粒的LBM模拟液固悬浮，通过修正的回弹边界保证运动边界的无滑移，通过动量交换计算流体和固体颗粒之间的作用力。Ladd方法中固体颗粒实际是由一个阶梯状的壳组成，固体颗粒边界像许多流体-固体格子链组成的，因此称为边界链法。

Ladd方法[8]的流固作用通过流体粒子和固体颗粒之间的动量交换实现，其优点是计算直观简洁，而且保证了流体的质量和动量局部守恒。但这种阶梯状的壳近似表示真实固体颗粒，往往会引入误差，而且在大多数情况下，这种误差会导致颗粒的动力学半径略大于其真实半径。只有在固体颗粒半径远大于LBM网格时，固体颗粒的动力学半径才和其真实半径相近。此外，Ladd方法中，固体颗粒内部由流体粒子占据，要求固体与流体密度比大于1。针对这一情况，Aidun等[9]提出不使用颗粒内部流体的方法，能够处理固体密度小于流体密度的问题。在Ladd方法的基础上，有学者进一步提出了满足局部质量守恒的方法[10]，并模拟了多颗粒沉降问题。然而，以Ladd为基础的方法均存在着边界形状与真实边界并不相符的问题，为了解决阶梯状近似，不少学者提出了改进方法，大部分都是采用处理曲面边界的方法,如曲面边界插值[11]。而这些处理方法虽然提高了边界处理的精确性，但在计算效率上却付出了巨大的代价，通常处理颗粒数目较少。但基于边界链法已可实现模拟含17845个颗粒的湍管流[12]。

2.2 格点法

Feng等[13]抛弃了固定颗粒边界的想法，把颗粒的边界当作可变形，结合浸入边界法，提供了一种新的处理运动流固边界的方法，这种方法避开了颗粒边界无法精确描述的困难，但对颗粒表面刚性系数的选择具有随意性。Noble等[14]提出了浸入运动边界法来处理流固耦合，通过在LBM演化方程中加入附加碰撞项，实现了对固体颗粒边界相对准确的描述，同时沿用了LBM的计算体系，保留了原算法的优点。Noble等在浸入运动边界法中引入了格子控制体和格子固含率的概念。格子固含率为格子控制体内颗粒体积占格子控制体体积的比重，固含率越高表示格子内包含的颗粒组分越多，相应的流体格子与颗粒的相互作用越强。将格子固含率作为重要参数进一步得出权重系数,实现对颗粒边界的光滑描述。权重系数存在多种构建形式[15-17]，均在一定程度上改进了模拟的精度。

Cook等[18]耦合浸入运动边界法和DEM模拟了颗粒相互作用占主导的颗粒流。冯云田等[19]在此基础上，引入大涡模拟LES(large eddy simulation)，成功地对高雷诺数颗粒流体系统进行了模拟。王利民等[20]发展了一种基于粒子-格子耦合方法的PR-DNS方法，该方法中颗粒模型采用结合硬球模型和软球模型各自优点的时间驱动硬球模型TDHS(time -driven hard-sphere)[21]进行求解，该算法将颗粒/流体密度比从2.5提高到了1500，成功复现了散式流态化和聚式流态化的典型现象[20]。在此基础上，周国峰等[22]成功模拟了二维和三维DKT(Drafting-Kissing-Tumbling)过程，验证了该算法的有效性，此外，还实现了Lee-Edwards边界条件与基于格点 LB -based PR-DNS的耦合[23]，克服了边界链接法无法满足伽利略不变性的问题，并成功预测了低雷诺数下颗粒悬浮体系的表观粘度。熊勤钢等[7]借鉴基于粒子-格子耦合方法的颗粒流体系统PR-DNS的思想，使用简化DEM代替硬球算法，利用中国科学院过程工程研究所研制的千万亿次超级计算系统 Mole -8.5，使用672个GPU进行了大规模并行计算，实现了全球最大规模颗粒流体系统PR-DNS，其中二维模拟包含1166400个颗粒，三维模拟包含129024个颗粒(图3)。

图3 颗粒悬浮系统大规模直接模拟粒子和流场分布[7]

周国峰等[24]分析了非均匀气固系统中颗粒受到的曳力情况，对气固流态化系统进行了三维大规模PR-DNS，成功复现了流动的非均匀结构，揭示了没有考虑非均匀结构的Wen &Yu曳力公式预测偏高，且受力方向存在偏差。

刘晓雯等[25,26]在中国科学院过程工程研究所研制的千万亿次超级计算系统 Mole -8.5E 上，使用192个Tesla K80 GPU运行17天，对4.83×109个流体网格和115200个固体颗粒的周期悬浮进行了直接数值模拟(图4)，探索了介尺度结构对颗粒统计属性、气固相间作用力以及传递特性的影响。

图4 大规模气固悬浮系统直接模拟瞬态结果[26]

刘晓雯等[25,26]发现颗粒脉动速度和局部平均无量纲曳力具备很强的尺度依赖性，提出了一种双峰分布形式以修正颗粒脉动速度分布，阐明了颗粒脉动速度偏离Maxwell分布的内在机理，并从微观层次上揭示了介尺度曳力的尺度依赖性，提出了以Fr数反映结构效应的曳力关联式[26]。

此外，耦合直接力浸入边界法的LBM应用于三维颗粒悬浮PR-DNS，如1024个热颗粒在流化床中的运动[27]；Rubinstein等[28]模拟了2000个颗粒在流化床中的运动；Eshginejadfard等[29]对305个颗粒槽道流进行了模拟。LB -based PR-DNS也应用于分析多颗粒沉降[30,31]或颗粒自旋及振动对流动的影响[32-34]、圆柱与柔性体耦合[35]、管道泄漏时土壤的流化现象[36]、泥石流[37]以及耦合传热的颗粒流问题[38]。通过对LBM中的弛豫时间进行修正，LB -based PR-DNS扩展到研究非牛顿流体颗粒流体系统中的颗粒流动规律[39-44]。通过重构颗粒的形式和颗粒间作用力，基于浸入边界法的LBM可以非常容易地处理诸如多边形颗粒[45]、可变形颗粒[46]和粘性颗粒[47]的迁移等，甚至可以分析岩土工程中的水力压裂问题[48]。以上均显示了基于 LB -based PR-DNS的计算策略在模拟颗粒流体系统中的巨大优势和潜力。

3 离散颗粒模拟

离散颗粒模拟(或颗粒轨道模型)将流体相处理为连续介质，固相处理为离散颗粒。在欧拉坐标系考察流体运动，在拉格朗日坐标系以及颗粒层次分析固相运动[49]。由流固耦合的实现方式不同，离散颗粒模拟一般可分三类，(1)单向耦合(one-way coupling)，只考虑流体对颗粒的作用，不考虑颗粒对流体的影响；(2)双向耦合(two-way coupling)，考虑流体和颗粒之间的相互作用；(3)四向耦合(four-way coupling)，除了考虑流体和颗粒之间的相互作用外，还考虑颗粒间碰撞。现在通常指的离散颗粒模拟，都是四向耦合的确定性颗粒轨道模型，在考虑颗粒间相互作用的情况下，采用曳力公式对颗粒与流体的相互作用进行描述。同时，基于牛顿第三定律可实现流场与颗粒的相互关联，即在各个控制体中，颗粒受到流体的作用力等于颗粒对流体施加的作用力。

在传统的离散颗粒模拟中，流体相的控制方程是基于体积平均Navier-Stokes (volume-averaged Navier-Stokes,VANS)方程或者Navier-Stokes方程的简化形式，通常采用SIMPLE(semi-implicit method for pressure-linked equations)算法等求解，这些处理方式都不易实现并行化,难以应用于大规模的计算中[4]。基于LBM的算法具有良好的并行性，王利民等[1,50]将LBM应用于离散颗粒模拟中流体相的流动(图5)，提出了基于格子 Boltzmann 的离散颗粒模拟( LB -based DPS)，其使用一种既考虑孔隙率又考虑气固相对滑移速度的修正LBM求解流体流动，取代了基于VANS方程的求解传统，固体运动采用TDHS求解，即固体颗粒的碰撞为硬球碰撞，在固定的时间步长内检索颗粒间是否发生碰撞以及更新颗粒运动。由于在颗粒浓度较高或颗粒温度较低的场景下应用硬球模型还有一定的应用局限性，王利民等[51]进一步发展了一种介尺度LBM-DEM耦合模型。熊勤钢等[52]也提出了一种快速的LBM-DEM耦合模型，但对力源项的添加方式有所不同。

图5 D2Q9计算流体格点和固体颗粒[1]

(1)

式中v为格点处的流体速度，vs为颗粒在当前格点的平均速度，-i表示与i方向相反。修正后的格子Boltzmann方程为[1]

(2)

式中εs和εg分别为格子控制体内部固体颗粒和气体的体积分率，且εs+εg=1。

LB -based DPS成功模拟了二维单孔射流鼓泡床[1,50,52]和三维鼓泡床(图6)的流态化过程，得到了与经验关联式和实验数据相吻合的模拟结果。李斌等[53,54]将 LB -based DPS进一步应用于二维喷动床和鼓泡床的流化过程，获得了与实验数据相吻合的结果。以上均验证了 LB -based DPS能够得到流固系统更多的细节信息和内在特性，此外，LB -based DPS在计算精度、耗费时间和效率之间能达到很好的平衡,可以直接扩展到与其他固相求解器耦合计算。

图6 三维流化床离散颗粒模拟流动结构的演化

4 双流体模型

目前 LB -based TFM研究文献较少，传统LBM的颜色模型、伪势模型、自由能模型和动力学模型只适用于不混溶的二元流体模拟，不能应用于流体相互渗透的多相流双流体层次模拟。王铁峰等[55]率先开展利用LBM构建双流体模型的研究，其利用无量纲方法变换格子Boltzmann方程并添加压力修正项来实现，但其恢复的宏观方程中，密度为两相的合密度而非分相密度；Sankaranarayanan等[56]从连续Boltzmann方程出发，考虑颗粒温度方程，提出了一种隐式格子Boltzmann模型，此模型也没有真正解决问题。实现 LB -based TFM的关键在于构造能恢复为VANS方程的格子Boltzmann模型，郭照立等[57]构造了一个可用于求解VANS方程的格子Boltzmann模型，该模型的平衡分布函数和作用力项都包含孔隙率，作用力的计算需要给定组成多孔介质的颗粒粒径，适合模拟均匀孔隙率的多孔介质内部流动，对于具有孔隙率梯度的复杂多相流问题，显式添加的离散力项通常导致模型不易收敛。Sungkon等[58]将VANS方程中带有体积分数的项均移到方程右端作为力源项，变换得到带有源项的N-S方程，并构造能恢复带有源项的N-S方程的LB方程，但是所得LB方程只有一阶精度。宋飞飞等[59]类似地将VANS方程转为带有源项的N-S方程，利用带有源项的作用力LBM模型发展了具有二阶精度的隐式格子Boltzmann方程，处理较复杂，不易并行。王利民等[1]考虑了体积分数和相间滑移速度对流体相流动的影响，基于浸入边界的思想，提出了修正的格子Boltzmann方程，用于描述气相流动，尝试拓展到双流体模拟。张金凤等[4]提出基于VANS方程的格子Boltzmann模型，通过Chapman-Enskog分析证明提出模型确实能够恢复到VANS方程。Blais等[60]指出体系存在孔隙率梯度的情况下，直接添加力源项修正压力梯度的LBM形式是不稳定的，并通过重构LBM中的碰撞算子以还原成VANS方程，其在大孔隙率梯度的情况下具备更好的稳定性。Höcker等[61]则是通过修正迁移步骤以提高稳定性，结合对流扩散的LBM实现了Rayleigh-Taylor不稳定性的模拟，但如何解决TFM下颗粒的紧密堆积还需进一步研究。虽然存在多种思路构建VANS方程的格子Boltzmann模型，但在使用上均存在一定局限性，将其应用到两相流的研究仍需进一步探索。

5 格子多相流体力学模拟软件

LB -based PR-DNS和DPS基于的浸入运动边界法以固含率将流体、颗粒和边界统一描述，保留了LB求解的并行性和迁移的简单性；同时格子Boltzmann方程基于笛卡尔网格求解能够处理复杂边界，避免了传统CFD绘制复杂构体网格的困难，在大规模并行条件下仿真规模已具备由实验室尺度向工业尺度的可行性。基于此，中国科学院过程工程研究所EMMS团队集成相关算法开发了以颗粒流仿真和多GPU并行为特色的高性能格子多相流体力学LMFD(lattice-based multi-fluids dynamics)模拟软件，如图7所示。LMFD基于C/C++开发，目前稳定版本支持二维 LB -based PR-DNS和DPS，以及三维简单构型的 LB -based PR-DNS计算。

图7 LMFD软件界面

颗粒流体仿真的主流商业软件有Barracuda，Fluent和开源的软件如MFiX和OpenFOAM。以上产品均基于传统的有限体积法开发，在并行性上受到很多制约。而LB基的主流求解器，如PowerFlow，XFlow和Palabos均以求解单相流为主，而以多相流仿真为主的LBM求解器还比较少。Seil等[62]耦合了Palabos与LIGGGHTS两个开源软件，实现了 LB -based PR-DNS求解器LBDEMcoupling。Götz等[63]在具备自适应网格加密的单相流程序WaLBerla基础上实现了 LB -based PR-DNS，用于研究河床中的沙丘形成。但以上程序的应用场景较为有限，且均无法实现 LB -based DPS。LMFD软件通过构建通用的LB基高性能气固两相流离散模拟平台，再逐步形成完备的多相流求解模块，将有望弥补这一空缺，在化工、冶金、能源、水利、环保和航空航天等相关领域具有广泛的应用前景。