基于VOF方法的可调参数对气泡聚并过程计算精度与成本的影响

周萍，蒋怡，廖义香，李家栋

(1. 中南大学 能源科学与工程学院，湖南 长沙，410083；2. 亥姆霍茨德累斯顿-罗森多夫研究中心，德国 德累斯顿，01328)

气泡聚并现象广泛存在于自然界中，并且在石油、化工、冶金、生物制药以及食品加工等行业中得到广泛应用[1]。气泡聚并过程通过改变气液两相间界面面积，影响相间的传热和传质，从而影响鼓泡塔、流化床等工业生产装置的产能与效率[2-4]。因此，对气泡聚并过程如气泡形状变化、聚并时间等进行详细了解，可以更好地理解气泡流动行为，对工业生产设备的设计和操作优化至关重要。

气泡聚并过程从机理上可分为气泡碰撞、薄膜排干和液膜破裂3 个阶段。在气泡碰撞阶段，2个气泡由于速度不同而互相靠近，如在同轴上升情况下，后面的气泡受前面气泡尾涡的影响而加速，碰撞后气泡之间形成一层液膜；在薄膜排干阶段，2个气泡之间的表面张力和毛细管力造成的虹吸使得液膜中液体流出，液膜变薄，气泡更加接近；在液膜破裂阶段，2个足够接近的气泡间液膜破裂，2个气泡聚并成一个更大的气泡[5-8]。

对气泡聚并过程进行研究时，重点关注的是当2 个气泡之间的液膜厚度为1.0×10-8～1.0×10-4m时液膜的排干速率[9]，这个阶段的时间间隔在毫秒甚至微秒级。这种微尺度气液两相流现象的实验研究对实验检测设备的采样频率与分辨率要求非常高，因此，受实验设备性能的限制，对气泡聚并过程的实验研究较少[8,10-11]。另一方面，随着计算流体力学(CFD)技术的快速发展，数值模拟方法在气泡聚并研究中的应用日益广泛。

气液两相相界面的准确描述是实现气泡聚并过程数值模拟的关键。相界面处理的方法分为两大类：欧拉界面捕捉(Front-Capturing)和拉格朗日界 面 追 踪(front-tracking)[12-13]。 前 者 包 括VOF[1,14-17]、Level-Set[18-20]以及Phase Field[21]等方法，而后者则包括Cell-Type、SPH(smoothedparticle-hydrodynamics)， PIC(particle-in-cell) 等 方法。近年来，还有学者将VOF 方法和Level-Set 方法这两种具有代表性的界面捕捉方法进行耦合，得到了CLSVOF 模型[22]。值得一提的是，格子玻尔兹曼方法(LBM)是一种将宏观流体力学和微观分子动力学二者关联起来的介观理论，当其被应用于气泡聚并行为研究时，在处理相界面附近的混合流体方面表现出明显优势[23]。由于LBM 方法受到计算量的限制，实际工程中的气泡聚并行为研究仍以宏观参数的求解为主。而VOF 方法因其固有的质量守恒特性、处理表面拓扑结构改变问题的合理性等特点，得到了广大学者的青睐[24-25]。

采用VOF 方法对气泡聚并的非稳态问题进行数值模拟时，常常会遇到计算量及精度之间的平衡问题，即如何设置计算流程中可调参数(网格尺寸、时间步长以及动量方程的迭代次数等)，可以实现在精度不受太大影响的前提下减少计算量。此外，在现有的基于VOF 方法的气泡聚并过程研究中，有关最大库朗数、动量方程循环次数、控制方程循环次数等可调参数的协调设置的研究较少[26-29]。

为此，本文以同轴2个气泡的聚并过程为研究对象，采用OpenFOAM 开源软件、VOF 方法以及PIMPLE算法对气泡的聚并过程进行数值模拟。在网格尺寸满足界面厚度精确度与网格独立性要求的前提下，探索计算流程中的可调参数(最大库朗数Comax、相方程循环次数nα和控制循环迭代次数npimple)对计算量与精度间平衡性问题的影响规律，以期为气泡聚并过程数值模拟时可调参数的设置提供指导。

1 数学模型与求解方法

VOF 方法基于欧拉法描述输运方程，假设多种流体(或相)互不相溶，引入相体积分数来实现对每一个计算单元相界面的追踪。

1.1 数学模型

不可压缩牛顿流体的连续性方程和动量方程分别如式(1)、(2)所示。

式中：∇为哈密顿算子；U为速度矢量，m/s；ρ为密度，kg/m3；P为压力，Pa；μ为黏度，Pa·s；t为时间，s；g为重力加速度，m/s2；F为表面张力，N/m3。

式中：σ为表面张力系数，N/m；α为相体积分数；下标l表示液体；k为界面的曲率。

式中：n为相界面单位法向量。

式中：δ为极小的非零项，以确保分母不为零，可由计算得到；N为网格数量；为所有网格的体积之和，下标i表示第i个网格。

在VOF 方法中，相分数αg和αl分别表示气相与液相在网格单元中的体积分数。建立并求解相分数输运方程获得气液两相体积分数，可以实现对气液相界面的捕捉。当系统中只存在气液两相时，仅考虑其中一相的相体积分数方程即可，另一相体积分数通过关系式αg+αl=1 求得。以液相的相体积分数输运方程为例，其计算式为

当网格单元内同时存在气液两相时，混合流体的密度和黏度计算式为：

1.2 迭代算法及可调参数

在瞬态气泡聚并数值模拟时，常用的控制方程迭代算法包括PISO(pressure-implicit splitoperator，压力隐式分裂算子)算法[30-31]和PIMPLE(结合PISO 与SIMPLE 算法[32-33]，SIMPLE 算 法为压力耦合方程组的半隐式方法算法[34-35]。PISO 算法采用速度预测—压力校正的计算策略，即在求解动量方程(速度预测)之后，基于压力泊松方程对速度场进行迭代校正(压力校正)。由于PISO 算法只进行一次动量方程求解，为避免对流项线性化处理时(面)通量信息滞后带来的影响，通常采用较小的时间步长来减少时间步进后速度场的变化幅度。PIMPLE算法采用类似的计算策略，但在每个时间步长内对整套控制方程进行多次迭代计算以更新通量信息，因此，该算法允许使用更大的时间步长。

图1 所示为基于OpenFOAM 的VOF 求解器interFoam 中PIMPLE 算法的计算流程。其中，对计算精度以及计算量有显著影响的主要可调参数包括最大库朗数Comax、控制方程循环次数npimple以及相方程循环次数nα。

图1 PIMPLE算法流程Fig.1 Flow chart of PIMPLE algorithm

Comax建立了时间步长与网格内流体速度及网格尺寸的约束关系[36]：

式中：Co为库朗数；Δt为时间步长，s；ΔV为网格单元体积，m3；uf为网格面心速度，m/s；下标f为网格单元面编号；Sf为网格单元法向面积矢量，m2。在模拟气泡聚并时，气液界面处的网格分辨率往往要求较高(即网格尺寸较小)，通过Comax限制Δt，有助于提高计算精度，但过小的Δt会增大计算量。

对于PIMPLE循环，首先根据前一时间步获得的速度更新面通量，之后对相方程进行求解：将Δt划分为nα个子时间步长Δtsub，即Δtsub=Δt/nα；在Δtsub至nαΔtsub时间范围内对相体积分数进行求解(以Δtsub为计算时间步长)，并更新气、液相体积分数；对Δtsub至nαΔtsub内网格单元的面质量通量ρUf·Sf进行积分，用于后续动量方程计算。可以看到，在不影响其他方程计算时间步长的情况下，相方程在Δt范围内嵌套了时间步长更小的瞬态计算；nα为该瞬态计算的迭代次数，nα越大，相方程的数值稳定性及精度越高，但计算成本也会相应增加。

PIMPLE 算法的动量预测采用PISO 算法的策略，通过一步预测、多次校正实现对速度场以及压力场的计算。npimple控制了PIMPLE 算法中整套控制方程的循环求解次数，当npimple=1 时，PIMPLE 算法等同于PISO 算法。若控制方程的循环次数达到npimple时，则PIMPLE 算法结束，计算进入下一时间步长。与nα类似，npimple的取值也需要从平衡计算精度与计算量的角度考虑。

2 物理模型与定解条件

本研究选取同轴气泡聚并过程开展数值实验，定量评价上述可调参数，包括Comax、nα及npimple对数值模拟的影响，同时得到较优的可调参数组合。本研究模拟的聚并过程为圆柱体容器中轴线上2个等径气泡的聚并过程，由于球形气泡及圆柱体容器均具有轴对称性，因此，可以假设气液两相的流动具有轴对称性，简化模型为二维计算域[37]。图2所示为同轴等径气泡聚并的物理模型。为了减小壁面对模拟结果的影响(壁面效应)，物理模型的宽度设置为dc=10d0，同时高度取h=20d0，以保证足够的气泡上升空间。在初始阶段，将直径d0=0.56 mm 的2 个静止气泡置于计算区域底部，两气泡中心的距离为1.18d0，下部气泡中心距离底面d0。

图2 计算区域示意图Fig. 2 Schematic diagram of computational domain

主要边界条件设置见表1，壁面为无滑移条件，顶部为压力出口。重力加速度g为9.81 m/s2，与y方向相反。

表1 主要边界条件Table 1 Main boundary conditions

表2所示为数值模拟工况的主要物性参数，系统处于大气压(101 325 Pa)条件中。在计算的初始阶段，设置时间步长为2×10-5s，同时采用自适应时间步长，根据网格尺寸、流体速度以及库朗数的相关设置，在计算过程中自动调整时间步长；此外，设置残差为1×10-8来保证控制方程的收敛性。

表2 气液两相主要物性参数Table 2 Main physical properties of gas and liquid

3 网格独立性分析与自适应网格

3.1 网格独立性分析

在网格独立性分析的计算工况中均采用均匀布置的结构化网格(uniform mesh, UM)。针对基准工况，设置4组不同尺寸的网格来验证网格的无关性，网格的尺寸及数量如表3 所示，4 组方案的分辨率均达到10-5m量级。本节模拟工况中可调参数取值分别为Comax=2、nα=1及npimple=2。

表3 均匀网格算例的网格尺寸及数量Table 3 Mesh size and number of cases with uniform grid

将体积分数云图中2个气泡边界连接成为一个封闭整体的这一时刻定义为气泡聚并时刻。图3所示为均匀网格算例的气泡聚并时刻与网格数量的关系。当网格数量从156 800 个(UM3)增加到278 258 个(UM4)时，2 个工况的气泡形状基本相同，且均在0.011 s 发生聚并，气泡聚并时刻趋于稳定。

图3 均匀网格算例的气泡聚并时刻与网格数量的关系Fig. 3 Relationship between bubble coalescence time of cases with uniform grid number

图4所示为不同网格数量下聚并时刻前上、下2 个气泡的上升速度随时间的变化曲线。从图4 可以看出：当网格尺寸从3.0×10-5m(UM1)减小到2.0×10-5m(UM3)，上、下2 个气泡上升速度的变化趋势以及大小都趋于一致；使用网格方案UM2时，上方气泡在0.006 s 后的上升速度与网格方案UM3、UM4 存在明显的区别。因此，综合考虑网格精度与计算成本，选择网格方案UM3 进行后续的数值模拟研究。

图4 上、下2个气泡聚并时刻前的上升速度Fig. 4 Rising velocity of the upper and lower bubbles before coalescence

3.2 自适应网格

液膜排干阶段是气泡聚并过程中的关键阶段，其发生的尺度在微米量级，对气泡周围网格的细化尺度提出了更高的要求，若采用均匀网格，则势必导致计算量巨大。为了更好地捕捉这一过程且尽可能减少计算量，本节在均匀布置的最优网格方案(网格方案UM3)基础上，使用自适应网格(adaptive mesh refinement, AMR)探究最佳网格细化次数，确定计算精度与计算资源间的平衡点。

自适应网格的基本思想是根据需要动态改变计算域的网格结构，在物理量变化剧烈的计算区域，采用空间尺度较小的精细网格予以计算；在物理量变化缓慢的区域，采用空间尺度较大的粗网格来计算，从而提高计算效率。针对气泡聚并过程，在气液界面处使用精细网格，气泡内部及气泡周围的液相区域使用较粗的网格。设计4组自适应网格方案，通过多次细化最小网格尺寸可达微米级，同时网格数量达到近200万个，具体参数见表4。

表4 自适应网格算例配置参数Table 4 Configuration parameters of cases with AMR

图5 所示为以液体体积分数αl≤0.5 为判别条件截取出的两气泡整体的中心位置处速度随时间的变化关系。从图5 可以看到，在计算初期(0～3 ms)，中心点速度为0 m/s，这是因为初期两气泡间存在一定距离，中心位置处于两气泡间的液相中，速度为0 m/s；到计算末期(17～20 ms)，中心位置速度再次降为0 m/s，这是因为气泡在上升过程形变为“帽状”，到后期形变幅度大，中心点处于两气泡下方的区域中，速度仍为0 m/s。当网格细化次数达到3 次以上(AMR3 和AMR4)，两气泡中心位置上升速度在各个时刻都十分接近，而细化次数为1 时中心位置处速度存在较大的差异，工况AMR2的气泡中心上升速度随时间的变化趋势大体上与AMR3、AMR4的相同。

图5 自适应网格算例两气泡中心位置速度Fig. 5 Velocity in the middle of two bubbles of cases with AMR

图6所示为4种网格细化条件下，气泡形状随时间变化的示意图。从图6可见，5 ms时4个工况的气泡形态基本相同，细化次数多的工况气液边界相对来说更加锐利。在25 ms时，AMR1中气泡已经发生聚并，因而与其他工况下的气泡形态区别较大。AMR3 和AMR4 的气泡形状无明显区别，但与AMR2相比，前两者的气液边界更加平滑。

图6 自适应网格细化次数对气泡形态的影响Fig. 6 Effect of refinement times of the AMR on bubble shape

综合考虑计算精度和计算成本，选择3次网格细化操作，并在此基础上开展可调参数对计算精度与计算成本影响的研究。

4 不同可调参数下聚并过程的数值模拟

在气泡聚并过程中，气泡间液膜厚度的变化是气液相间相互作用的体现，直接反映了两气泡相对速度的变化。另外，液膜厚度是气泡聚并模型中的关键参数[38]。因此，选择液膜厚度作为评价可调参数影响的指标，其值取截出的2个气泡整体的竖直中心线上液面间距离。计算效率采用数值模拟的CPU 时间进行衡量，即CPU 完成0.03 s(物理时间)的工况数值模拟计算及内核系统调用所耗的时间总和。需要注意的是，本文4.1 节与4.2节的计算是在2个不同的计算平台上完成的，计算平台的性能参数分别如下：1) CPU 主频3.4 GHz，内存储器容量为128 GB；2) CPU主频2.5 GHz，内存储器容量192 GB。针对Comax、nα及npimple这3个可调参数的数值实验方案如表5所示。

表5 可调参数数值实验Table 5 Numerical experiment for adjustable parameters

4.1 可调参数对液膜厚度及CPU时间的影响

图7 所示为最大库朗数Comax对气泡间液膜厚度以及CPU时间的影响。从图7(a)可以看到，在计算刚开始时，液膜厚度为0.18d0，随着计算过程进行，液膜厚度缓慢减小；随后液膜厚度的减小速度逐渐增大，表明两气泡互相接近的速度也相应变大；随着液膜厚度进一步减小，其变化速度及气泡接近速度逐渐放缓。对比4 个工况发现，Comax＞0.01 的3 组工况间液膜厚度随时间的变化曲线区别较大，而当Comax从0.05 减小到0.01 时，2组工况的液膜厚度十分接近。另一方面，Comax与液膜变薄速度成反比，即Comax越小，相同时间内液膜厚度越小(变薄速度大)。这主要是由于Comax较大时，时间步长也相应较大，在数值扩散的影响下流体在网格单元上的移动出现滞后，从模拟结果上表现为气泡聚并过程变慢，液膜变薄速度减小；当通过减小Comax来提高计算精度时，能有效改善流体移动的滞后现象。从图7(b)可以发现，CPU 时间随最大库朗数增大而减小。综上所述，Comax=0.05可以较好地兼顾计算精度及计算成本。

图7 最大库朗数对气泡间液膜厚度及CPU时间的影响Fig. 7 Influence of the maximum Courant number on liquid film thickness and CPU time

图8所示为不同相方程循环次数nα下气泡间液膜厚度变化及所耗CPU 时间。观察图8(a)可以发现，当nα为1时，液膜厚度变化曲线明显位于其他4个工况下的变化曲线的上方；当nα≥3时可以获得稳定的液膜厚度变化曲线。根据1.2 节的分析可知，增大nα相当于使用更小的时间步长求解相方程，有助于提高计算精度。nα越大，液膜变薄速度越大。结合图8(b)可知，nα=3可以在保证计算精度的同时，节约计算成本。

图8 相方程循环次数对气泡间液膜厚度及CPU时间的影响Fig. 8 Influence of number of cycles of the phase equation on liquid film thickness and CPU time

图9 所示为不同控制方程循环次数npimple对气泡间液膜厚度变化及CPU 时间的影响。从图9(a)可以看到，计算精度随着npimple增大而提高，当npimple≥8 时，可以获得稳定数值解。从计算耗时来看(图9(b))，CPU 时间与npimple成正比，因为增大npimple意味着PIMPLE 算法的外部循环次数增多，因此，消耗的CPU 时间增加。从平衡计算量与计算精度的角度出发，取npimple=8。

图9 控制方程循环次数对气泡间液膜厚度及CPU时间的影响Fig. 9 Influence of number of cycles of the governing equations on liquid film thickness and CPU time

4.2 较优的可调参数组合

从4.1 节可知，3 种可调参数分别取Comax=0.05、nα=3 以及npimple=8 时，能够在各自数值实验组内有效平衡计算精度与成本。本节选取计算精度最高时的3 个可调参数的最值(各参数取值范围参考4.1节)进行组合，即(Comax，nα，npimple)=(0.01，5，10)。基于同轴气泡聚并数值模拟，在求解精度方面将其与(Comax，nα，npimple)=(0.05，3，8)进行比较，数值实验的相关设置见表6，其他设置与4.1节中的设置相同。

表6 可调参数组合数值实验Table 6 Numerical experiment for optimal adjustable parameters

图10、图11 所示分别为Toptimal和Tmaxmin这2 个工况下液膜厚度随时间的变化曲线及气泡形状的比较。以工况Tmaxmin为基准，液膜厚度比在1附近浮动，液膜厚度的平均相对误差为3.21%。观察气泡形状，发现2 组工况几乎没有区别。综上所述，认为(Comax，nα，npimple)=(0.05，3，8)是针对同轴气泡VOF数值模拟的较优的可调参数组合。

图10 Toptimal和Tmaxmin这2个工况中液膜厚度随时间的变化规律Fig. 10 Variation law of liquid film thickness with time in cases of Toptimal and Tmaxmin

图11 Toptimal和Tmaxmin这2个工况的气泡形态比较Fig. 11 Comparison of bubble shape between two cases of Toptimal and Tmaxmin

5 结论

1)Comax与液膜变薄速度成反比，即Comax越小，数值模拟时Δt相应减小，流体在网格单元上以更大速度进行输运，导致液膜变薄速度增大，反之亦然；此外，Comax越小，计算精度越高，相应的计算成本也会提高；基于数值实验获得其最优取值为Comax=0.05。

2) 液膜减薄速度随nα增大而增大，这主要是由于增加nα相当于使用更小的时间步长求解相方程，由此导致流体输运速度增大，加速液膜厚度减小；计算精度随nα增大而提高，计算成本亦然；基于数值实验获得其最优取值为nα=3。

3)npimple与液膜减薄速度成正比；计算精度随着npimple增大而提高，计算成本也随着控制方程的迭代次数增大而增加；基于数值实验获得其最优取值为npimple=8。

4) 通过与计算精度最高的参数组合比较，获得基于VOF 方法模拟气泡聚并时的较优参数组合为(Comax，nα，npimple)=(0.05，3，8)。