1 数值计算方法

湿颗粒碰撞过程的直接数值模拟通过流体体积模型（volume of fluid，VOF）模型追踪气液界面，连续表面张力模型（continuum surface force，CSF）模型计算界面的表面张力，两者共同作用，控制液膜的形变过程。颗粒运动则采用六自由度模型（six degrees of freedom，6 DOF）结合重叠网格技术进行描述，为了模拟颗粒的反弹运动并提取分析数据，本文开发并应用了用户自定义函数。

1.1 VOF模型

采用VOF模型来模拟气液两相流体的流动问题，通过求解体积分数的对流方程实现相界面的瞬态追踪。在VOF模型中，通过跟踪不同相的体积分数来描述不同相的位置和界面形状[式(1)]。

∂α∂t+∇⋅(αu)=0(1)

式中，α为液体的体积分数，%，α=1处的网格充满液体，α=0处的网格充满气体，0<α<1的网格区域含有自由面；t为时间，s；u为流体的速度矢量，m/s；∇为梯度算子。

假设流体为不可压缩流体，且气液两相的密度保持不变，涉及质量与动量守恒。

连续性方程为式(2)。

∇⋅u=0(2)

动量方程为式(3)。

∂(ρu)∂t+u⋅∇(ρu)=-∇p+∇⋅[μ⋅(∇u+∇uT)]+g+FSV(3)

式中，ρ为流体的局部平均密度，kg/m3；p为流体的压力，Pa；μ为动力黏度，Pa·s；g为流体自身重力，N；FSV为表面张力的等价体积力，N。

局部平均密度ρ和动力黏度μ表达式为式(4)和式(5)。

ρ=αρl+(1-α)ρg(4)

μ=αμl+(1-α)μg(5)

式中，下角标l为液体相；下角标g为气体相。

1.2 连续表面张力模型

连续表面张力模型能够将液体界面上的力转换为周围的体积内连续分布的力，被广泛应用于液滴流动模拟中。界面处表面张力可以表示为式(6)。

Fsf=σκ∇φ(6)

式中，Fsf为表面张力，N；σ为表面张力系数；∇φ为界面切向的梯度算子；κ为表面曲率。

界面处的单位法向量计算如式(7)、式(8)。

n=-∇φ|∇φ|(7)

n=nwcosθ+twcosθ(8)

式中，θ为气-液界面在壁面处的接触角，(°)；nw和tw分别为与壁面法向和切向的单位向量。得到界面的曲率为式(9)。

κ=∇⋅∇φ|∇φ|(9)

1.3 颗粒运动模型

控制颗粒运动位移的方程为式(10)。

dupdt=1mpFp+Fv+Fsf+Fg(10)

式中，Fp为压差阻力，N；Fv为黏性阻力，N；Fg为重力，N；mp为颗粒质量，kg；up为颗粒速度，m/s。

为进一步分析湿颗粒碰撞过程中的能量损失，可以将其分为以下四部分：压差阻力引起的能量损失比Ep、黏性阻力引起的能量损失比Ev、由湿颗粒间非完全弹性碰撞引起的能量损失比Ec、由表面张力引起的动能耗损比Esf。并将颗粒1刚接触液膜时刻的初始动能定义为E0。计算见式(11)~式(13)。

Ep=(∫h0htFpdh+∫hthmFpdh)/E0Fp=-∫Spiz⋅ndS(11)

Ev=(∫h0htFvdh+∫hthmFvdh)/E0Fv=μ∫S∇uz⋅ndS(12)

Ec=(12mut2-12mus2)/E0(13)

式中，h0为颗粒1接触液膜时质心高度，μm；ht为颗粒1与颗粒2碰撞时质心高度，μm；hm为颗粒1碰撞后反弹到最高点的质心高度，μm；ut为颗粒1碰撞前的速度，m/s；us为颗粒1碰撞后的速度，m/s。

1.4 重叠网格

本文使用的重叠网格方法广泛应用于模拟移动网格的边界运动问题。重叠网格技术的工作原理包括挖孔、重叠最小化和供体搜索。首先，生成互不约束的背景网格和重叠网格，满足N-S求解器的三阶离散化。在求解的迭代过程中，流场数据通过内嵌边界进行传输。其次，确定执行各自职责的单元，包括挖孔单元、重叠边界单元（接收单元）和计算单元（供体单元）。挖孔单元是指通过挖掘操作丢弃的一些烦琐单元。重叠边界单元主要围绕挖孔单元，并被视为每个重叠网格中的边界。重叠边界单元与计算单元之间的关系通过插值确定。具体实现流程如图1所示。

图1   重叠网格计算求解流程