LMFD应用｜LBM求解多相混合物模型，实现含滑移速度的混溶两相流模拟
摘要：本文提出了一种基于压力的格子Boltzmann方法(PLBM)直接求解混合物模型(MM)中的方程，即PLBM-MME，该模型将各相视为相互渗透的连续体，与二次相体积分数方程耦合后即可实现各相速度不同的混溶两相流模拟。与有限体积方法(FVM)相比，该方法在模拟扩散界面时能够更严格地保证质量守恒性，且能够发挥LBM在处理均匀网格时高并行的优势，具备向更复杂封闭理论拓展的潜力。
研究背景
研究背景
多相流从流型的角度可以被划分为分层流和分散流及其混合形式，其中分散流具有更加充分的相间接触和传递作用，以及更加显著的非均匀性和多尺度特征。颗粒悬浮液流动作为一种重要的分散两相流形式，其模拟方法包括欧拉-欧拉(EE)方法和欧拉-拉格朗日(EL)方法，其中EE能够模拟更稠密的悬浮液流动。具体而言，EE方法包括混合物模型(MM)和双流体模型(TFM)。当相间规律未知或其可靠性受到质疑时，TFM可能不可行，MM作为一种更简单的模型可以像TFM一样执行，同时求解的变量数量比TFM少。
混相模型及传统求解
等温不可压两相流的混相模型守恒方程[1]包括混相质量和动量守恒方程及二次相体积分数守恒方程：
其中ρm = αpρp + αgρg表示混合物密度; ρp和ρg分别为颗粒相p和连续相g的常数密度; αp和αg表示两相的体积分数; cp = αpρp/ ρm为颗粒相的质量分数; um = (αpρpup+ αgρgug) / ρm表示混合物的质量平均速度. uslip =up − ug表示二次相p相对于主相g的滑移速度; 压力p被视为各相和混合物的统一压力; DMp为二次相体积分数方程的扩散常数，常被视为0; τGm表示粘性应力张量。
当不考虑滑移速度时，MM被简化为均匀多相流模型，在各相速度不同的情况下，需额外预先插值求解滑移速度矩阵[2]，或滑移速度的经验关联代数表达式来封闭MM。Fluent软件内默认的滑移速度代数表达式是基于Manninen等人[1]的局部瞬时平衡假设，其中uslip与二次相颗粒的加速度a方向一致，大小上基于Schiller和Naumann[3]的阻力关系和颗粒雷诺数雷诺数Rep = dpρg |uslip| / μg，其中滑移速度需要迭代求解。
此外，由式(1)和式(3)可得，混合物的质量平均速度需满足以下条件：
该条件也可以使用体积平均速度表示，即为两相连续性方程：
传统FVM求解MM时无法使用基于密度的求解器，必须使用基于压力的求解器，因而必须求解压力泊松方程(PPE)。PPE的离散求解依赖连续性方程，其在不可压单相流中常用的零散度速度假设，在MM中已不严格成立，因此需要额外考虑引入混相密度的瞬态项或改换考虑两相连续性方程。
基于混合物模型的压力基格子Boltzmann法
含外力项的格子BGK模型，采用标准D3Q19的格子速度和权重，其演化方程为：
其中松弛时间τ满足：
含Dtρm的PLBM-MME完整恢复了混相模型，在[4]工作的基础上，Fi和压力计算将各相表观密度重新定义为混合物密度：
基于混合物速度零散度假设的混相模型不含Dtρm，其Fi和压力分别为：
所提出的PLBM-MME在耦合二次相体积分数方程时，采用显式的FVM进行离散。MM中体积分数的空间离散方案通常包括两种类型的界面类型，即扩散和锐利/扩散界面。锐利/扩散界面适用于捕捉介于锐利选项和扩散选项之间的轻度锐利界面。具体而言，采用一阶迎风方案模拟扩散界面，采用基于通量校正（FCT）的QUICK格式模拟锐利/扩散界面。
该方法已集成在由中国科学院过程工程研究所EMMS团队开发的多相流模拟软件LMFD（Lattice-based Multi-Fluids Dynamics）中，并基于GPU做高并性的实现。
收敛性测试
通过制造解方法（MMS）验证所提出含Dtρm的PLBM-MME二阶收敛性。
引入L2范数来评估数值和理论混合速度之间的误差。
不同晶格尺度和混合密度分布的误差如图1所示。混合物速度在所有情况下都表现出二阶收敛性，除了在较高粘度和较大混合物密度梯度时略有精度损失。
图1制造解瞬态算例的收敛阶数
开闸式异重流
该流动为在移除闸门后两种不同密度的静态流体在封闭通道中混合过程，该过程涉及前缘运动和界面不稳定性等独特物理现象。使用800×1×100的网格计算这一过程，初始时刻左右两侧密度分别为0.9992 和 1.0008 kg/m3.运动粘度, ν, νg, νp,均为1.0 × 10-4 m2/s.重力为 1.0 × 10-4 m/s2。
图2 开闸式异重流算例的初始条件
图3 在t=10Ts后，不同方法模拟开闸式异重流的比较
(a) Fringer等人[5] FVM；(b) PLBM-MME不含Dtρm；(c) PLBM-MME含Dtρm
与其他两种方法相比，没有Dtρm的PLBM-MME表现出略大的混合物密度耗散。使用这两种PLBM捕获的Kelvin–Helmholtz现象与FVM获得的非常相似，表明所提出的LB方案，当与二次相的体积分数方程结合时，即使在运动粘度足够小的情况下，也可以准确地模拟低密度比条件下的两相混合过程。
Rayleigh-Taylor不稳定性
Rayleigh-Taylor不稳定性（RTI）不仅是研究不混溶流体分层流动行为的基准算例，也是EE模型界面评估的重要测试之一，该算例使用64×1×320的网格实现拟二维模拟。为了便于比较，LBM的格子单位和传统FVM的物理单位之间直接应用1:1的转换比，即一个晶格单元和时间步长分别对应于1.0m和1.0s，所有其他物理量都根据这种关系相应地导出。两相密度和重力加速度分别设置为3.0, 1.0, 3.742218×10-7，其余条件如下表所示：
表一 RTI各算例不同参数设置
基于传统FVM的结果来自商业软件Ansys Fluent 2024。Fluent®求解器的设置配置为双精度，LBM和传统FVM程序输出的物理量的有效数字设置为10位。传统FVM的体积分数参数被指定为隐式的，时间步长为1000s，每个时间步长最多500次迭代，扩散界面采用一阶迎风方案，锐利/扩散界面采用可压缩方案。通过测试更小的时间步长和更多的迭代步数，验证了FVM的当前时间推进设置，以最大限度地减少隐式公式引起的数值误差。图4显示了使用三种不同方法的四种情况下二次相体积分数的演化过程。
图4 三种方法模拟四种RTI算例中二次相体积分数的演化过程对比
(a) FVM; (b) PLBM-MME不含Dtρm; (c) PLBM-MME含Dtρm
算例（1）和（2）之间的比较表明，独立降低主相的运动粘度会导致下降液滴两侧形成更薄、更细长的尖峰。此外，传统FVM和不含Dtρm的PLBM-MME在算例(1), (2)和(3)下的演化过程在界面下降速度方面表现出相对相似性。这说明了所提出的LB方案恢复的混合物方程与FVM离散求解的混合物方程之间的一致性。然而，含Dtρm的PLBM-MME演化过程中，密度较大的流体的界面下降速度明显慢于其他两种方法。上述差异理论上可以用式(5)来解释。使用一阶迎风格式对体积分数进行空间离散化会产生数值扩散，导致DMp变得显著，从而导致混合物速度的非无散度条件。然而，严格来说，混合物模型中三个守恒方程的正确耦合不应直接将混合物速度的散度视为零。
此外基于算例（4）中的锐利/扩散界面锐度特征，图5进一步提取了三种方法的二次相体积分数为0.5的界面轮廓。可以看出，三种方法的重液界面下降速度基本一致。这仍然需要从式（5）的角度来解释。由于这三种方法都有效地抑制了数值耗散，因此可以忽略方程（5）中的DMp。此时，式(1)等价于混合物速度的无散度条件，这是三种方法模拟结果整体一致的重要原因。此外，由于使用了更精确的三阶精度格式，FCT方案比可压缩方案更好地抑制了界面耗散，而可压缩方案则基于二阶精度方法。
图5 算例（4）中不同方法得出的RTI界面的演变过程
为了进一步验证上述关于式(5)的解释，在算例(1)中整场每10000个时间步分别计算混合速度散度和质量守恒误差∂tρm+∇·(ρmum)的Frobenius范数。范数表示为：
其中，∂t和∇算子都使用二阶中心差分方案进行近似。FVM的∂tρm是通过改变时间步长来获得的，以确保FVM和LBM之间误差计算的一致性。具体来说，在采样时间周围的6秒窗口内，FVM计算的时间步长被修改为1秒，而不是1000秒。误差随时间的演化过程如图6所示。
图6 算例(1)中RTI整场误差范数随时间演化过程
可以观察到，含Dtρm的PLBM-MME的结果与不含Dtρm的FVM和PLBM-MEM的结果存在显著差异。在整个演化过程中，含Dtρm的PLBM-MME的混合速度零散度误差明显高于其他两个。然而，含Dtρm的PLBM-MME质量守恒误差明显较低。上述差异在重流体从上壁分离阶段尤其明显，即100000到200000步。这是由于含Dtρm的PLBM-MME完全恢复了混合物密度的守恒方程，这是传统FVM和不含Dtρm的PLBM-MEM所缺乏的特征。综上所述，含Dtρm的PLBM-MME在混合物模型中表现出更严格的质量守恒性质。
此外，比较了在RT算例下采用双精度算法的两个PLBM-MME的计算性能(MLUPS)。很明显，时间导数的计算对计算速度的影响很小，而数值扩散的抑制和滑移速度的迭代求解显著影响了计算效率。在PLBM中，在GPU框架上实现的FCT算法需要至少两个额外的同步措施来纠正一阶方案。总体而言，在均匀网格上，LBM的并行显式推进方法在运动粘度变化、非线性滑移速度的迭代求解和伪扩散抑制的条件下保持了稳定的并行性能。
表二计算机性能对比
创新点总结
提出了一种新的基于压力的格子Boltzmann方法来求解混合模型方程。
PLBM-MME支持将任意的滑移速度关联式纳入混合物模型。
PLBM-MME在不连续的密度间断中实现了更严格的守恒特性。
PLBM-MME有效地捕获了Kelvin–Helmholtz和Rayleigh-Taylor不稳定性。
参考文献：
[1] M. Manninen, V. Taivassalo, S. Kallio, On the Mixture Model For Multiphase Flow, VTT Technical Research Centre of Finland, 1996.
[2] H. Zhang, L. Wang, Dynamic energy-minimization multi-scale heterogeneous continuum evolution model for gas-solid fluidization, Chem. Eng. J. (2024) 480.
[3] Z. Naumann, L. Schiller, A drag coefficient correlation, Z. Ver, Dtsch., Ing 77 (1935) e323.
[4] S. Fu, Z. Hao, L. Wang, A pressure-based lattice Boltzmann method for the volume-averaged Navier-Stokes equations, J. Comput. Phys. (2024) 516.
[5] O.B. Fringer, M. Gerritsen, R.L. Street, An unstructured-grid, finite-volume, nonhydrostatic, parallel coastal ocean simulator, Ocean Model. 14 (2006) 139–173.
本文参考自：Z. K. Hao and L. M. Wang, A pressure-based lattice Boltzmann method for mixture model equations of multiphase flow, J. Comput. Phys. 539, 114247, (2025).
