同轴式微通道内气泡生成特性数值模拟及试验验证

王 月，徐 博,2，王艳丽

(1.北华大学机械工程学院，吉林吉林 132021； 2.燕山大学机械工程学院，秦皇岛河北 066000；3. 吉林省工业技师学院，吉林吉林 132021)

引言

微细气泡的等效直径通常小于100 μm相对常规气泡体积甚小，具有许多优良特性，如气泡比表面积较大、在两相间传质效率高、上浮速率慢、表面自带活性自由基等[1-2]。微细气泡因具有如上所述特性而逐渐被广泛应用于臭氧投加处理污水、医学检测、美容保健、去果蔬农残等诸多领域，且在各领域中微细气泡都表现出很好的潜在应用及科研价值[3-5]。

近年来，国内外有许多学者进行了大量的研究，并得到大量的研究成果。谭立新[6]研究了水利工程中水-气二相流的一些基本理论，并提出了水-气二相流的数值模拟方法与实测数值验证的方式。赵铎[7]通过在大型多相流试验环道上进行气-水二相流试验及不同流量的数值模拟得出气-液二相的分布及流速。柏超[8]通过理论分析方法解释微细泡(液滴)生成过程的一些现象与流型转换机理，并得出预测微细气泡(液滴)大小的数量级关系式。袁希钢等[9]针对T形微通道，运用数值模拟与试验的方式研究了黏度和表面张力对气-液二相流型分布的影响。杨丽等[10]针对流动聚焦模型，通过数值模拟研究了微管道内液滴的形成机理及影响因素。冯俊杰[11]研究了2个共轴线上升气泡之间的连接和聚并行为，以及气泡在湍流场中的上升和破碎行为。

国外研究主要以剪切式研究为主，较少涉及同轴式的研制。CHEN W B等[12]运用界面描述法研究微细气泡(液滴)形成时的表面形状，仿真研究了在浸没(同轴)式喷嘴处的两相气泡的形成和相变。整个气泡由蒸气和液相组成，假定内部气相被薄的冷凝层包围，当两相气泡增长时蒸气部分冷凝。界面描述法用于描述气泡形成的动力学，气泡界面处的气液传热效率与不同流量下的气液压差有关。QIAN D等[14]在针对T形管的研究中得出，气体段塞长度随着表面气体速度的提高而增加，并且表面液体速度降低。LI X等[15]通过数值模拟与试验验证研究T形管道中微液滴形成机制。WONG V L等[16]运用数值模拟方法研究T形管道微流体黏度对气泡体积的影响，得出选择聚合物浓度和分子量来控制剪切变形和弹性，可以调整液滴尺寸和形成机制。

本研究以同轴式微通道中的微细气泡生成过程为研究对象，基于COMSOL开展微通道中的气、液二相流相互作用过程数值模拟，获得气泡生成过程中不同时刻的轮廓图像，分析不同气体压强对气泡生成频率、生成体积等生成特性的影响趋势，并基于微细气泡生成特性测试平台，验证微细气泡生成过程数值模拟结果的准确性。

1 微细气泡生成过程仿真建模

1.1 COMSOL简介

COMSOL对模型定义是非常灵活的，可以用常数、通用函数、表达式以及插值函数等来表达材料属性和边界条件等。COMSOL具有多物理场应用模式，可以解决气液二相间多物理场耦合问题，用户可以自主选择需要的物理场进行数值计算。针对流体仿真，该软件为用户提供了许多优化完善的控制方程，供用户计算使用。

1.2 控制方程

由于微细气泡流的雷诺数很小，因此将连续相设置为层流流动，在COMSOL Multiphysics中可用二相流水平集方法研究两种流体之间的变化。数值模拟中涉及到的控制方程包括3个，即连续性方程、Navier-Stokes方程、水平集方程。

=▽·[-pI+μ(▽u+(▽u)T)+σkδn]

(1)

▽·u=0

(2)

(3)

式中，ρ—— 流体密度

u—— 流体速度向量

t—— 时间

p—— 总压力

I—— 单位矩阵

μ—— 流体动力黏度

σ—— 表面张力系数

k—— 流体界面的曲率

δ—— 二相流界面函数

n—— 指向液滴的界面单位法线向量

σkδn—— 作用在二相流界面上的表面张力

φ—— 水平集函数

γ—— 重新初始化参数

ε—— 界面厚度控制参数

式(1)中的流体密度ρ和动力黏度μ可表示为:

ρ=ρl+(ρg-ρl)φ

(4)

μ=μl+(μg-μl)φ

(5)

式中，ρl—— 液体密度

ρg—— 气体密度

μl—— 液体动力黏度

μg—— 气体动力黏度

针对微通道内二相流体的流动和相互作用，利用COMSOL的水平集模块求解方程。在利用COMSOL进行数值模拟过程中，要求用户针对润湿壁边界条件设定滑移长度，根据摩擦力定义，在此暂不考虑液体在润湿壁面的相对滑动，遂将滑移长度初值取0。

1.3 几何模型建立

在实际的试验过程中，气体入口是直径60 μm的毛细金属管，液体通道是边长为10 mm的方形玻璃管。为减轻仿真计算量和仿真周期，先建立二维仿真模型，再在仿真中分别运用轴对称旋转和阵列得到三维立体仿真计算结果，其模型的每部分尺寸如图1所示。

图1 几何模型尺寸

1.4 边界条件设置

在同轴式微通道中，忽略能量交换，设定初始流场中A区域为空气、B区域为液体。由于微细通道中流体雷诺数很小，Re<<1，因此流体流态为层流，针对同轴式微通道中气泡生成的仿真模型，设置边界条件如下：

1) 入口边界条件

气液二相流体在入口边界的速度均匀，无梯度分布，气体和液体均以各自恒定的流速进入微通道。

2) 出口边界条件

气液二相混合流体从固定出口流出时的压力恒定，为1个大气压，pout=0 Pa或pout=1 Pa。

3) 对称边界条件

对于整体同轴微通道仿真模型而言，其中进行的流体流动数值模拟计算量较大，因此为减轻计算量，取整体模型的1/4部分作为仿真模型。另外，为方便控制方程识别，需要将仿真模型设置为轴对称。

4) 润湿壁边界条件

将同轴微通道壁面添加润湿性属性，壁面润湿程度可通过设置接触角θ值确定。润湿壁的滑移长度可设置为1 μm，如图2所示。

图2 同轴聚焦结构边界条件设置

1.5 初始条件设置

将毛细管道设置为气体，外部管道设置为液体，图2的初始界面为气液分界面。根据实际试验中获得的液体流速和气体流速设置仿真中的初始条件，液流流进管道的速度u1=0.049 m/s，气流流入毛细管道的速度u2=2.97 m/s。设置接触角为90°，重新初始化参数设置为0.5，界面厚度控制参数为ls.ep_default(表示细化网格后最小网格尺寸的1/2)，初始条件参数如表1所示。

表1 同轴聚焦模型初始条件

表1中取用的值为室温(20 ℃)下的数值，接触角选择为气泡最初形成时气液二相分界面的切线与同道壁面的夹角。此外，仿真过程中的气体流量由气泡生成体积除以气泡的生成时间求得。

1.6 网格单元划分

由于连续相与分散相交界处网格单元尺寸较小，不能以较大尺寸的网格单元进行划分，而全局网格若均设置成很小的网格单元，会增加仿真的计算量，致使仿真效率降低。综合考虑以上条件，在此全局上采用物理场控制网格、单元格大小选择及细化。局部分散相与连续相的交界处选择自动细化网格，减少不必要的计算。在仿真之前，预先划分一次网格，共59695个三角形网格单元，30311个网格顶点。在仿真过程中，系统自动细化网格，最终有116676个三角形网格单元和58944个网格顶点。图3表示仿真前后模型的网格局部示意图。

图3 仿真前后网格对比图

1.7 求解器设置

在同轴聚焦模型的求解中主要运用瞬态求解器，方法选择向后迭代求解器，一致初始化部分选用向后欧拉法。在误差估计结束后，完成数值模拟部分相关设置。

2 微细气泡生成特性试验

2.1 试验设计及试验平台搭建

为探究同轴式微通道中微细气泡的生成效果和生成特性，搭建了一套微细气泡制备试验系统。试验中气泡生成的载体为同轴机械式微流控芯片，普通氮气瓶为芯片提供气体，并利用YQD-4双级氮气减压阀对气体压强进行调节；Harvard 704500微量注射泵作为芯片的液体供给单元，按需调整芯片液体入口流量。将Olympus CX21显微镜和Phantom v12.1高速摄像机组合使用，同时具备高速拍摄和显微观测的功能，用于观察和采集气泡生成过程中的图像，试验系统的原理示意图如图4所示。

图4 微细气泡生成特性测试平台示意图

由于数值模拟结果需要与试验形成对比关系，因此试验中需要收集到规律性较明显的几组参数，防止由于数值模拟算法的问题出现不收敛或无规律现象。试验采取固定液体流量，改变气体压强的试验条件，探究气体压强对微细气泡脱离体积和生成频率的影响规律。

2.2 试验过程及方法

为探究不同液体流量和气体压强下的气泡生成特性，利用所搭建的微细气泡制备试验系统开展气泡生成特性研究试验，所设定气压与液体流量记录在表2中。

试验过程中，液体流量由Harvard 704500微量注射泵提供，该注射泵可精确调节液体流量。气体压强可通过连接氮气瓶的减压阀和精密气压调节器精准调节。设置高速摄像机拍摄帧率为10000 fps，记录每个微细气泡的生成演变时间，试验中的相关参数及数值，如表3所示。

表2 试验变量参数及其数值

表3 试验中的相关参数及其数值

2.3 试验数据的选择

为了可以和数值模拟结果进行清晰地对比，需要选择规律性较为明确的试验数据作为仿真的输入控制参数。表4记录了液体流量为150 mL/h时，不同气体压强pg下气泡的脱离体积Vg和生成时间t。可以明显看出，气体压强变化对气泡的生成时间及体积影响较大，可以作为数值模拟输入量仔细研究气泡形成过程的变化。

表4 不同气体压强下气泡生成时间及体积

3 仿真与试验结果对比分析

3.1 微细气泡形成过程

细致观察气泡的形成过程是研究微细气泡生成的良好方法，而使用数值模拟方法可以更便捷的观察微细气泡的形成过程。图5a为相同的控制量下通过数值模拟方法计算出的气泡形成过程，图5b为试验中相同控制量下气泡的形成过程，其中为了与试验图像形成对比选取了相同气泡生长阶段的图形比较。

图5 不同时刻下微细气泡生成过程

将这些图像进行对比分析，可以发现整个气泡生成过程同样可分为3个阶段：膨胀阶段、轴向拉伸阶段和颈缩脱离阶段。经过对比容易发现，相同的气泡生成阶段，数值模拟中的形成时间较短，并且颈缩脱离阶段的颈缩量相对较长。

3.2 气泡生成频率

由于数值模拟方法需要在实际工程中运用，因此必须找到数值模拟方法与实际气泡生成频率之间的区别。研究中选取了液体流量为140 mL/h，通过改变气体压强探究不同工况下气泡生成频率的变化规律。表5记录了不同气体压强下气泡的生成频率，观察表中气泡生成频率的试验值fgt和仿真值fgs可以发现，随着气体压强的增大，气泡的生成频率变快。

表5 不同气体压强下气泡的生成频率

为了更清晰的观察实际试验与数值模拟计算的气泡生成频率之间的区别，根据表5绘制实际气泡生成频率与数值模拟计算的气泡生成频率随气体压强增加时的变化趋势，如图6所示。可以明显看出，经过数值模拟计算出的气泡生成频率相对实际试验气泡的生成频率较大，因此实际工程运用数值模拟方法分析气泡生成频率时应考虑相应的计算差值。

图6 不同气体压强下的气泡生成频率

3.3 气泡的体积

将液体流量设置为恒定，通过调节气体压强获得固定流量不同气压下微细气泡的脱离体积，并将不同工况下气泡的体积记录在表6中。

表6 恒定流量不同气压下气泡的体积

表6中的数据记录了不同气体压强下气泡脱离体积的试验值Vgt和仿真值Vgs，为了鲜明对比二者的变化趋势，将表中数据进一步处理，得到不同气体压强下气泡体积的变化曲线，如图7所示。

图7 不同气体压强下气泡生成体积

从图7中可以明显看出，恒定液体流量不同气体压强下气泡体积的试验值和仿真值均随气体压强的增大而减小，而且仿真结果普遍略小于实际试验值。

4 结论

通过本研究可以得出：采用COMSOL Multiphysics软件数值模拟机械式同轴型微管道内气泡的形成，与实际试验规律明显一致，但数值上具有一定差值。数值模拟在计算中的优势明确，可以显著的观察气泡形成中的具体变化，经过研究分析得到数值模拟方法在微细气泡生成的研究分析上有着较大的运用前景。