摘要:为了研究单液膜气泡单层液膜包裹气体的特殊结构形式和内外壁面均受气液表面张力的特殊力学形式，在拉格朗日框架下采用无网格移动粒子半隐式法并基于表面自由能表面张力模型，建立了单液膜双气液界面表面张力模型，从而实现了单液膜气泡振荡变形过程中的复杂界面计算和捕捉。在此基础上对2个单液膜气泡的聚并和连接过程进行了模拟分析，获得了典型的流动现象和液膜变形特征与规律，发现减小表面张力系数或增大黏性系数均会减弱气泡变形过程中表面张力项的变形主导作用。为此，提出了凹点切线法用于计算连接型气泡的液膜夹角，明晰了连接型气泡的形状。计算结果可为工业消泡技术提供一定的理论依据。

单液膜气泡是由一层很薄的液膜在表面张力作用下包裹一定量气体形成的，其内外壁面均存在气液表面张力作用，力学作用和几何结构均比较特殊。气泡聚并行为是指两个气泡相互接触后合二为一的过程。气泡连接是指两个气泡相互接触且含有一层液膜的稳定存在形式。

理论分析方面，Ramakrishnan等提出了气泡形成的两阶段模型；杨林涛等根据表面吸附和两相流体力学理论，对液体中气泡体积进行理论分析后得出了液体中气泡体积的计算方程，并计算出水中气泡的大小。实验研究方面，张建生等用高速摄影技术对水中气泡运动规律进行了详细研究；Chan等使用原子力显微镜、光干涉条纹图像及测量聚并时间3种方式对液滴和气泡的排液过程及依赖于时间的相互作用关系进行了研究；顾汉洋等利用双平行探针技术和摄像方法对水平直管内气液两相流中单气泡的形态特征进行了研究。数值模拟方面，常用方法有Front-tracking、VOF、Level-set和LBM等。Li等采用VOF方法研究了表面张力系数和黏性系数对液体中上升气泡运动情况的影响；张淑君等采用VOF方法中的界面重构方法模拟研究了液体中不同放置方式气泡之间的相互作用；Sussman等将VOF与Level-set方法相结合，模拟了无黏气泡的上升和球形气泡在近壁处的破裂；Amaya采用LBM方法研究了单气泡或者少量气泡的变形行为。

上述研究都是基于单气液界面的表面张力模型对液体中气泡的形态进行分析，关于空气中单液膜气泡，即液膜双气液界面在表面张力作用下的变形行为涉及较少，这是因为单液膜气泡的聚并和连接等行为过程伴随着剧烈和复杂的相界面变化，对数值方法在相界面捕捉和计算稳定性等方面要求较高。移动粒子半隐式法(MPS)在处理计算区域存在巨大变形的流动问题时具有较大优势，迄今已在诸多工程和科学研究中得到了广泛应用，例如Sun等编程实现了二维情况下水中气泡自由上升运动的数值模拟。

本文基于MPS方法，对二维空气中单液膜气泡在表面张力作用下的振荡变形、聚并及连接过程进行数值模拟和研究。

1数值方法

1.1移动粒子半隐式法

移动粒子半隐式法的基本思想是在求解区域内部及其边界上布置一系列离散的粒子来代表所要研究的宏观流体，采用时间层推进法获得整个流场的流动信息。在拉格朗日坐标系下表征流场的流动信息，不可压缩流体的控制方程为

MPS方法用核函数w形成的粒子作用模型来离散控制方程，本文所使用的核函数为

为2个粒子间的距离；r_{e}为影响域，取值为3.1 l_0，l_0$为粒子间的初始距离。控制方程中的梯度算子和拉普拉斯算子可以基于核函数进行离散。

1.2表面张力模型

采用表面自由能模型来实现对单液膜气泡在表面张力作用下变形过程的数值模拟，其基本思想是通过构建粒子间的相互作用力来模拟粒子所受到的表面张力。定义2个MPS粒子之间的势能为P(r)，则表面张力为

式中:n为2个粒子之间的单位方向矢量。r小于临界值时，粒子间为斥力；r大于临界值时，粒子间为引力。粒子间势能公式为

式中:r_{\min}为斥力、引力的分界位置，即表面自由能的极值点，取值为1.5 l_0；C为修正系数，可根据流体的物性求得。

由式(4)和式(5)得表面张力计算公式为

1.3单液膜双气液界面表面张力模型

针对单液膜气泡的特点，假设气泡外部空气压力恒定为0，模拟时不使用流体粒子模拟外部大气，忽略其对液膜流动的影响；假设气泡内部空气压力均匀恒定，采用不可压缩流体粒子来填充代替，其物性参数为空气。两部分气体中间的液膜采用黏性流体粒子来离散计算，参数设置为具有一定表面张力系数的液体，具体实施方案如图1所示。

（图1 单液膜双气液界面表面张力模型示意图）

假设粒子总数为N，i为粒子编号，i=1\sim m时为液膜粒子，i=(m+1)\sim N时为内部空气粒子，则相对应的密度ρ(i)和运动黏性系数ν(i)的表达式分别为

对于液膜自由面上的表面粒子，当r<3.1 l_0时，其表面张力计算公式为

对于黏性项即拉普拉斯项，其计算公式为

1.4数值方法验证

表面张力总是使液体表面具有收缩至表面积最小的趋势。本文选取典型验证算例之一——二维方形气泡在表面张力作用下收缩振荡、在黏性耗散作用下最终稳定在圆形状态来进行验证，计算参数见表1，初始粒子布置，粒子间距为0.3 mm，气泡边长为18mm，其中内部黑色部分为气体，不同深浅粒子代表液膜(由外到内)，总厚度为0.9mm。计算时间步长取为\Delta t=10 μs。t=0~0.02s时，方形气泡4个角逐渐收缩，t=0.02s时近似收缩为菱形，之后再由菱形逐步变为圆形达到稳定状态。

气泡平均压力随时间的变化，内层液膜的平均压力和内部气体粒子的平均压力在不同时刻其大小和变化频率几乎一致，约为13Pa。中间液膜的平均压力大于内层液膜和气体粒子的平均压力，约为203Pa。外层液膜在程序中被判断为表面粒子，其平均压力为0(大气表压)。上述气泡变形过程及压力计算结果均符合空气中单液膜气泡的物理特征，验证了本文单液膜气泡计算模型的有效性。

表1 方形气泡的粒子物性参数