0 引　言

气泡在复杂微通道内的流动现象广泛存在于自然界和工业生产过程，如沸腾现象、页岩气开采、燃料电池的使用等。例如在燃料电池中气泡间的碰撞和合并时常发生，气泡在液体中的运动和相互作用直接影响相间的接触时间和接触面积，并改变气液间的传热传质^[1]，进而影响其使用性能^[2-4]。因此对气泡进行微通道内的流动研究以及工业生产过程中的优化具有十分重要的意义。

学者们对复杂微通道内的气泡流动问题进行了大量的实验研究。Chen等^[5]研究了高黏度流体中气泡在孔板通道内上升时的形状变化，发现气泡穿过孔板通道后变为“月牙状”。Corapcioglu等^[6]提出采用定量方法预测气泡在粗糙多孔介质中的上升速度，其预测结果表明气泡的上升速度经过很短的时间后会达到平衡。Dawson等^[7]研究了在恒定体积流量中孔板微通道内气泡的流动过程，发现随着流速的增加，气泡的后尖端速度发散且气泡尾部变薄向中间收缩回主体。

以上实验结果得到了气泡运动过程中的直观现象和宏观规律，为人们认识气泡在通道内的运动规律提供了理论基础。然而由于实验方法难以捕获微孔内的流动细节，采用数值模拟来研究两相流动问题也逐渐成为人们研究气泡运动的一种常用方法^[8-9]。Sussman等^[10]应用Level-set方法，研究了单个气泡上升过程中的界面形变以及破裂现象，发现随着表面张力和黏度比的增加，气泡越不容易发生变形和破裂现象；而随着雷诺数的增加气泡的变形会越来越剧烈。Sun等^[11]利用格子Boltzmann方法模拟了单个气泡在含三个半圆形障碍物的通道内上升的问题，发现了障碍物的存在可以影响气泡的运动形态，使其运动轨迹发生偏转。Lou等^[12]通过两相流格子Boltzmann模型，模拟研究了气泡在浮力作用下上升并通过复杂微通道的问题。研究发现障碍物表面润湿性越强，气泡变形越剧烈、气泡终端速度越小；且随着障碍物半径的增大，气泡的终端速度呈线性减小。

以上研究揭示了一些气泡在复杂微通道内的运动机理，但主要集中于单气泡的研究，较少涉及多气泡在复杂微通道内的运动行为，且大部分研究只能考虑较小的气液密度比。因此，本文选用Liang等^[13]基于Allen-Cahn相场理论提出的大密度比格子Boltzmann模型，对两个上下排列的气泡在含孔板结构的微通道内的运动形态进行研究；主要调查了 $ Eo $ 数、气泡之间的相对大小、两气泡之间的距离以及气泡和孔板之间的距离等因素对两个气泡变形、合并、分裂的动力学行为以及气泡穿过孔板的质量和通过时间的影响。

1 数值方法

Liang等^[13]提出的大密度比LBM模型采用了两个分布函数描述指标参数和速度/压力的演化过程，其中指标函数满足如下Allen-Cahn方程：

$ \begin{array}{c}\dfrac{\partial \phi }{\partial t}+\nabla \cdot \left(\phi \boldsymbol{u}\right) = \nabla \cdot \left[M\left(\nabla \phi -\mathrm{\lambda }\boldsymbol{n}\right)\right]\end{array} $

(1)

式中： $ \phi $ 为序参数，气相 $ \phi = 0 $ ，液相 $ \phi = 1 $ ， $ 0 ＜ \phi ＜ 1 $ 为气液间的界面；M为迁移率； $\boldsymbol{u}$ 为气相运动速度； $ t $ 表示时间； $\boldsymbol{n} = \nabla \phi /|\nabla \phi |$ 表示界面法向单位方向； $ \lambda = 4\phi \left(1-\phi \right)/W $ ，其中W为界面厚度。速度和压力等流动宏观量满足Navier-Stokes方程。对应Allen-Cahn和Navier-Stokes方程的LBM演化方程 ${f}_{i}\left(\boldsymbol{x},t\right)$ 和 ${g}_{i}\left(\boldsymbol{x},t\right)$ 可分别表示为：

$ \begin{array}{c}{f}_{i}\left(\boldsymbol{x} + \boldsymbol{e}_{i}{\delta }_{t},t + {\delta }_{t}\right)-{f}_{i}\left(\boldsymbol{x},t\right) = -\dfrac{{f}_{i}\left(\boldsymbol{x},t\right)-{f}_{i}^{{\rm{eq}}}\left(\boldsymbol{x},t\right)}{{\tau }_{f}}+{\delta }_{t}{{{F}}}_{i}\left(\boldsymbol{x},t\right)\end{array} $

(2)

$ \begin{array}{c}{g}_{i}\left(\boldsymbol{x}+\boldsymbol{e}_{i}{\delta }_{t},t + {\delta }_{t}\right) - {g}_{i}\left(\boldsymbol{x},t\right) = -\dfrac{{g}_{i}\left(\boldsymbol{x},t\right)-{g}_{i}^{{\rm{eq}}}\left(\boldsymbol{x},t\right)}{{\tau }_{g}}+{\delta }_{t}{\mathit{G}}_{i}\left(\boldsymbol{x},t\right)\end{array} $

(3)

${f}_{i}\left(\boldsymbol{x},t\right)$ 为粒子的密度分布函数， $i = \mathrm{0,1},2,\cdots ,Q$ ； $ Q $ 为离散速度方向个数， $\boldsymbol{x}$ 和 $ t $ 分别表示粒子的空间位置和演化时间； $ {\tau }_{f} $ 表示无量纲时间， $ {\tau }_{f} $ 与迁移率M之间的关系为 $ M = {c}_{s}^{2}\left({\tau }_{f}-0.5\right){\delta }_{t} $ ； ${g}_{i}\left(\boldsymbol{x},t\right)$ 为粒子流场的分布函数； ${g}_{i}^{{\rm{eq}}}\left(\boldsymbol{x},t\right)$ 为相应的平衡态分布函数； ${\mathit{F}}_{i}\left(\boldsymbol{x},t\right)$ 为力项； ${\mathit{G}}_{i}\left(\boldsymbol{x},t\right)$ 为外力项； $ {\tau }_{g} $ 为与运动黏度 $\upsilon$ 相关的无量纲松弛时间，其关系式为： $\upsilon = c_s^2 \left({\tau }_{g}-0.5\right){\delta }_{t}$ 。

在两相流系统中，由于气液相界面的存在，运动黏度 $\upsilon$ 不再是定值，其由序参数和气液两相的运动黏度的线性方程得到（ ${\upsilon }_{l}$ 为液相运动黏度， ${\upsilon }_{g}$ 为气相运动黏度），计算式： $\upsilon = \phi ({\upsilon }_{l}-{\upsilon }_{g})+{\upsilon }_{g}$ 。

平衡态分布函数 ${f}_{i}^{{\rm{eq}}}\left(\boldsymbol{x},t\right)$ 和 ${g}_{i}^{{\rm{eq}}}\left(\boldsymbol{x},t\right)$ 为：

$ \begin{array}{c}{f}_{i}^{{\rm{eq}}} = {\omega }_{i}\phi \left(1+\dfrac{\boldsymbol{c}_{i}\cdot \boldsymbol{u}}{{c}_{s}^{2}}\right)\end{array} $

(4)

$ \begin{array}{l}{g}_{i}^{{\rm{eq}}}\left\{\begin{array}{l}\dfrac{p}{{c}_{s}^{2}}\left({\omega }_{i}-1\right)+\rho {s}_{i}\left(\boldsymbol{u}\right),\qquad i = 0\\ \dfrac{p}{{c}_{s}^{2}}{\omega }_{i}+\rho {s}_{i}\left(\boldsymbol{u}\right),\;\qquad\qquad i = 1～8\end{array}\right.\end{array} $

(5)

$ \begin{array}{c}{s}_{i}\left(\boldsymbol{u}\right) = {\omega }_{i}\left[\dfrac{\boldsymbol{c}_{i}\cdot \boldsymbol{u}}{{c}_{s}^{2}}+\dfrac{{\left(\boldsymbol{c}_{i}\cdot \boldsymbol{u}\right)}^{2}}{2{c}_{s}^{2}}-\dfrac{\boldsymbol{u}\cdot \boldsymbol{u}}{2{c}_{s}^{2}}\right]\end{array} $

(6)

其中： $ {c}_{s} $ 为格子声速， $ {c}_{s} = c/\sqrt{3} $ ； $\boldsymbol{c}_{i}$ 为离散速度； $ {\omega }_{i} $ 为权重系数； $\boldsymbol{c}_{i}$ 和 $ {\omega }_{i} $ 由不同的格子模型决定。本文采用D2Q9模型^[14]进行数值研究，其速度模型如图1所示，离散速度 $\boldsymbol{c}_{i}$ 为：

$ \begin{array}{l}\boldsymbol{c}_{i} = \left\{ \begin{array}{l}c\left(\mathrm{0,0}\right),\qquad\qquad\qquad\qquad\qquad\qquad\qquad i = 0\\ c\left\{\mathrm{cos}\left[\left(i-1\right)\dfrac{{\text{π}} }{2}\right],\mathrm{sin}\left[\left(i-1\right)\dfrac{{\text{π}} }{2}\right]\right\},\;\qquad\qquad i = 1～4\\ \sqrt{2}c\left\{\mathrm{cos}\left[\left(2i-1\right)\dfrac{{\text{π}} }{4}\right],\mathrm{sin}\left[\left(2i-1\right)\dfrac{{\text{π}} }{4}\right]\right\},\qquad i = 5～8\end{array}\right.\end{array} $

(7)

式中， $ c = {\delta }_{x}/{\delta }_{t} $ 为格子速度， $ {\delta }_{x} $ 和 $ {\delta }_{t} $ 分别表示格子步长和时间步长。本文中取 $ {\delta }_{x} = {\delta }_{t} = 1 $ 。D2Q9模型^[14]对应的权重系数分别为：当 $ i = 0 $ 时， $ {\omega }_{i} = 4/9 $ ；当 $i = 1～4$ 时， $ {\omega }_{i} = 1/9 $ ；当 $i = 5～8$ 时， $ {\omega }_{i} = 1/36 $ 。

力项 $ {\mathit{F}}_{i} $ 和外力项 $ {\mathit{G}}_{i} $ 的表达式为：

$ {\mathit{F}}_{i}\left(\boldsymbol{x},t\right) = \left(1-\dfrac{1}{2{\tau }_{f}}\right)\frac{{\omega }_{i}\boldsymbol{c}_{i}\cdot \left[{\partial }_{t}\left(\phi \boldsymbol{u}\right)+{c}_{s}^{2}\lambda \boldsymbol{n}\right]}{{c}_{s}^{2}} $

(8)

$ {\mathit{G}}_{i}\left(\boldsymbol{x},t\right) = \left(1-\frac{1}{2{\tau }_{g}}\right)\times {\omega }_{i}\left[ \frac{{\boldsymbol{c}_{i}\cdot (\boldsymbol{F}}_{s}+\boldsymbol{G})}{{c}_{s}^{2}}+\left({\rho }_{l}-{\rho }_{g}\right)\boldsymbol{u}\nabla \phi :\dfrac{\boldsymbol{c}_{i}\boldsymbol{c}_{i}}{{c}_{s}^{2}}\right] $

(9)

式（8）中的 ${\partial }_{t}\left(\phi \boldsymbol{u}\right)$ 可根据显式欧拉法得到：

${\partial }{t}\left(\phi \boldsymbol{u}\right) = \frac{\phi \left(x,t\right)\boldsymbol{u}\left(x,t\right)-\phi \left(x,t-{\delta }_{t}\right)\boldsymbol{u}\left(x,t-{\delta }_{t}\right)}{{\delta }_{t}}$

(10)

式（9）中： $ \;{\rho }_{l} $ 和 $ \;{\rho }_{g} $ 表示液相和气相密度； $\boldsymbol{G}$ 表示体积力项即浮力项；表面张力 $ {\boldsymbol{F}}_{s} = {\mu }_{\phi }\nabla \phi $ ，其中化学势 $ \;{\mu }_{\phi } = 4\beta \phi \left(\phi -1\right)\left(\phi -0.5\right)-\kappa {\nabla }^{2}\phi $ 。 $ \kappa $ 为气泡表面张力强度系数； $\; \beta $ 是与界面厚度 $ W $ 和表面张力 $ \sigma $ 相关的物理参数，关系式为： $\kappa = 1.5\sigma W，\beta = 12\sigma /W$ 。此外，通过对分布函数的统计，可计算得到宏观变量序参数、密度和速度：

$ \begin{split}\phi =& \sum _{i}{f}_{i} \\ \rho =& \phi \left({\rho }_{l}-{\rho }_{g}\right)+{\rho }_{g} \\ \rho \boldsymbol{u} =& \sum _{i = 0}^{8}\boldsymbol{c}_{i}{g}_{i}+0.5{\delta }_{t}\left({\boldsymbol{F}}_{s}+\boldsymbol{G}\right)\end{split} $

(11)

2 物理问题

本文所研究的物理问题如图2所示。在计算区域为 $ {L}_{x}\times {L}_{y} = 200\times 450 $ 的微通道内有两块大小相同的水平孔板，其中：孔板下表面距离计算区域下端 $ h = 160 $ ，两孔板间的距离 $ w = 20 $ ，孔板厚度 $ \delta = 12 $ 。在孔板下方放置两个直径分别为 $ {d}_{1} $ 和 $ {d}_{2} $ 且圆心在同一条垂直线上的气泡，处于上方的气泡圆心位置为（ ${x}_{1c}$ ， ${y}_{1c} $ ），其圆心距离孔板的距离为 $ {l}_{1} $ ；下方气泡的圆心位置为（ $ {x}_{2c} $ ， ${y}_{2c} $ ），上下两气泡圆心间距离为 $ {l}_{2} $ 。初始时，气泡密度设置为 $\; {\rho }_{g} = 1.0 $ ，运动黏度设置为 $ {\nu }_{g} = 1.0 $ ，其周围的液滴密度设置为 $\; {\rho }_{l} = 1\;000.0 $ ，运动黏度设置为 $ {\nu }_{l} = 0.1 $ ，气液界面厚度 $ W = 4.0 $ ^[15]。两气泡在浮力的作用下（浮力 $ G = ({\rho }_{g}-{\rho }_{l})g $ ）穿过孔板并流出管道， $ g $ 取 $ -3.91\times {10}^{-6} $ 。通道上下为周期边界，左右壁面以及孔板表面均为无滑移边界。文中未考虑壁面润湿性对气泡运动行为的影响，关于此问题，我们在后续的工作中会进一步研究。本文中所涉及的特征数主要为Eötvös数（ $Eo$ ）和雷诺数（Re），其中 $ Eo $ 数的表达式为 $ Eo = g({\rho }_{l}-{\rho }_{g}){D}^{2}/\sigma $ ，其表示浮力与表面张力的比值；重力作用下的雷诺数的表达式为 $ R{e}_{Gr} = ({\rho }_{l}-{\rho }_{g})\sqrt{g{D}^{3}}/{\mu }_{l} $ ，代表浮力和黏性力的比值。为了定量地分析各个因素对气泡流过多孔介质后的剩余质量的影响，定义气泡剩余质量比 $ {M}_{s} = 1-M/{M}_{t0} $ ，其中 $ M $ 是气泡到达通道顶部后在障碍物下方残留的气泡质量， $ {M}_{t0} $ 是模拟初始的气泡总质量，气泡质量取格点密度小于 $ ({\rho }_{l}+{\rho }_{g})/2 $ 处的气泡质量进行统计。气泡平均速度 $ {v}_{b} $ 统计表达式为：

$ {v}_{b} = \frac{\mathrm{\Sigma }{V}_{z}\phi }{\mathrm{\Sigma }\phi } $

(12)

另外，取 $ \sqrt{{d}_{1}^{2}+{d}_{2}^{2}} $ 为特征长度，定义公式为：

$ {t}^{\mathrm{*}} = t\sqrt{\frac{g}{\sqrt{{d}_{1}^{2}+{d}_{2}^{2}}}} $

(13)

$ {v}^{\mathrm{*}} = \frac{{v}_{b}}{\sqrt{{g}\sqrt{{d}_{1}^{2}+{d}_{2}^{2}}}} $

(14)

对计算时间和速度进行无量纲化，其中： $ t $ 为格子时间， $ {d}_{1}、{d}_{2} $ 为对应气泡的直径。下文中除验证外均使用上述无量纲量。

3 模型验证 3.1 拉普拉斯验证

首先通过Laplace定律验证模型的正确性。Laplace定律指出，当系统达到平衡态时，气泡内部压力与液相压力的差和表面张力成正比，和气泡半径 $ r $ 成反比。对于气泡可表示为：

$ \Delta p = {p}_{{\rm{inside}}}-{p}_{{\rm{outside}}} = \frac{\sigma }{r}$

(15)

其中 $ \sigma $ 为表面张力。

在本文验证中，初始气泡半径选20、25、30、35、40这五种情况，并考虑不同表面张力σ = 0.2、0.4、0.6、0.8、1.0，模拟结果如图3显示。气泡半径的倒数与表面张力成正比，模拟结果与Laplace定律吻合，证明了本文所使用模型的准确性。

3.2 气泡在光滑通道内上升行为

气泡在光滑壁面通道内受浮力上升的问题与本文所要研究的问题相似, 且已有学者对该问题进行了大量研究, 因此本小节将采用该算例验证模型的正确性。模拟计算区域为 $L_x\times L_y = 80\times 300$ ，液相密度 $\; {\rho }_{l} = 1\;000.0 $ ，气泡密度 $ \;{\rho }_{g} = 1.0 $ ，气泡半径 $ r = 20 $ ，界面厚度 $ W = 3.0 $ ，初始时刻气泡圆心设在（ $L_x/\mathrm{2}$ ，40）处。边界条件设置：左右为无滑移固壁，上下为周期边界。

考虑 $Eo = 5、10、20、40$ 四种工况， $ Re $ 数固定为17.7。表1列出了不同工况下，本文所使用的LBM方法和其他数值方法得到的气泡稳定形状对比，可以看出随着 $ Eo $ 数的增大，气泡从球形变为上凹的球帽形。本文得到的结果与Takada等^[16]使用双流体LBM模型得到的模拟结果以及Ren等^[17]使用相场LBM法模拟小密度比气泡得到的模拟结果高度相似，且在大密度比工况下得到了数值稳定的结果。之后统计气泡终端速度，如表2所示，并与Alizadeh 等^[18]和VOF 方法^[19]进行对比。本文结果与两者结果一致，故本模型适用于研究气泡动态行为。

4 数值结果与分析

本文中主要研究Eötvös数（Eo）、气泡之间的相对大小和气泡的初始相对位置对气泡在孔板微通道内流动过程的影响。

4.1 $ \mathit{E}\mathit{o} $ 数的影响

$ Eo $ 数是一个表征流体所受浮力和表面张力之比的无量纲参数，本小节主要研究 $ Eo $ 数对气泡在孔板通道内的动力学行为的影响，主要考虑了四种不同 $ Eo $ 数的情况，分别为 $ Eo = $ 5、10、20、40。在数值模拟中，通过调节表面张力 $ \sigma $ 的大小而得到不同的 $ Eo $ 数，其他参数设置如下：上下两气泡大小相等（ ${d}_{1} = {d}_{2} = 40$ ），初始位置为（ ${x}_{1c} $ ， ${y}_{1c} $ ）=（ $ {L}_{x}/\mathrm{2}$ ，100）和（ ${x}_{2c} $ ， ${y}_{2c} $ ）=（ ${L}_{x}/\mathrm{2} $ ，40）。

图4给出了四种 $ Eo $ 数工况下气泡穿过通道的流动过程瞬时图。图4（a）所示为 $ Eo = 5 $ 的工况，如图所示，在浮力作用下两气泡逐渐靠近孔板，当上方气泡接近孔板时速度下降，随后气泡缓慢挤入孔板通道并发生形变。下方气泡在 $ {t}^{*} = 6.71 $ 时刻与上方气泡接触，此时两气泡的接触位置处于中间；随着两气泡越来越近以及上方气泡的形变越来越严重，除上下气泡的中间以外，上方气泡与下方气泡在两侧也开始合并，导致上下方气泡之间包裹两个极小液泡（ $ {t}^{*} = 7.10 $ ）；随后在表面张力作用下液泡从两侧被排出，此时通道内的两个气泡在孔板附近完全合并形成一个完整的大气泡（ $ {t}^{*} = 9.98 $ ），大气泡通过孔板通道后在孔板正下方残留两个小气泡，而其余部分通过孔板后在孔板上方近似恢复为椭球形（ $ {t}^{*} = 18.13 $ ）。

图4（b）所示为 $ Eo = 10 $ 的工况，此时气泡表面张力相对于浮力减小，因此两个气泡在上升的过程中都有较大的形变。上方气泡下表面向上凹陷，在两端形成尾翼，下方气泡水平拉长且上表面向上凸起与上方气泡下表面近似贴合； $ {t}^{*} = 7.09 $ 时刻下方气泡的顶点与上方气泡中间接触，接着上方气泡两端向下弯曲，下方气泡上表面两侧与其接触，气泡间形成两个液泡，因为表面张力相对浮力减小，所以此时形成的液泡比 $ Eo = 5 $ 的工况下生成的大。由于孔口的阻力小于两侧，下方气泡质量主要从中间的接触点并入上气泡， $ {t}^{*} = 9.19 $ 时刻上气泡与两端断裂，其两端并入下气泡两端，液泡也被气泡从断口排出。随后在 $ {t}^{*} = 11.04 \;$ 时刻两个气泡合成了一个主气泡（能够排出的气泡），在部分质量穿过孔口后，主气泡运动速度增大与孔板正下方部分质量再次断裂，这些质量最终残留在孔板下方。主气泡在穿过孔口后由于拉伸再次破裂出一个小气泡，该小气泡尾随主气泡上升（ $ {t}^{*} = 20.76 $ ）。对比图4（a）和图4（b）可以发现，随着 $ Eo $ 数增加，气泡在穿越孔板前后的形变更明显。

图4（c）为 $ Eo = 20 $ 得到的气泡运动过程。从图中也发现随着 $ Eo $ 数的持续增大导致气泡变形和断裂更严重。例如在气泡的上升过程中，上方气泡尾端形成了更长的尾翼，下方气泡也更扁长。而更多的接触面积和更小的表面张力也让上下两个气泡间形成更多的接触点（ ${t}^{\mathrm{*}} = 7.36$ ），两个气泡中间存在着更多更大的液泡，尽管此时气泡在过孔板通道前依旧可以将液泡挤出，但也导致主气泡和其两端的连接容易断裂（ ${t}^{\mathrm{*}} = 11.30$ ），残余在孔板下的气泡质量增加。主气泡在通过孔板后也分裂出一个更大的尾气泡（ $ {t}^{\mathrm{*}} = 20.76 $ ）。而当 $ Eo $ 数增加到40时（图4（d）），上下气泡发生了多次接触、合并与破裂的过程。在主气泡通过孔板通道后，气泡内部还夹带着一个小液泡与其上升了一段距离才被排出（ $ {t}^{*} = 21.55 $ ）。从图中可以看出四种工况下气泡首次接触的时间分别为 ${t}^{*} = 6.71、7.09、7.36、7.61$ ，接触时上气泡都已进入孔道，主气泡到达计算区域上方的时间如图5所示。因此随着 $ Eo $ 数的增大，气泡接触时间越晚，气泡间接触面积增大，接触点随之增多；接触点的增多并不能促进两气泡的合并，反而导致中间生成更多液泡，加剧气泡分裂，气泡通过孔板通道需要的时间也越长。气泡到达顶部的所需时间随着 $ Eo $ 数的增大而增多，并且当 $ Eo $ 数大于10后呈线性上升的趋势。

上述四种工况下， $ Eo $ 数对两个气泡通过孔板通道的剩余质量比 $ {M}_{s} $ 有影响，为了分析这种影响，本文对不同 $ Eo $ 数下的气泡剩余质量进行统计。由图5对 $ {M}_{s} $ 曲线的拟合发现，气泡剩余质量的变化趋势随 $ Eo $ 数的增加呈负指数级地减小；减小是因为随着 $ Eo $ 数的增大，气泡表面张力相对于浮力减小，气泡难以维持初始形状，其在上升和变形的过程中更容易包裹液泡，并更容易发生破裂，导致残留在孔板下方。而随着 $ Eo $ 数的持续增大，气泡的最终剩余质量趋于稳定；这是因为尽管气泡被拉伸，但气泡初始位置与孔板间的距离导致气泡还未充分形变便到达了孔板位置，有一半多的质量处在孔口正下方且最终通过孔板。

4.2 气泡相对大小的影响

本小节研究了气泡相对大小对气泡通过孔板微通道的影响。数值模拟中参数设置为 $ \sigma = 0.624 $ ；为了方便比较，选择孔口间距 $ w $ 作为参照尺寸，上下两气气泡大小 $ {d}_{1}:{d}_{2} $ 分别设为 $ 2.4w:1.5w $ 、 $ 2.0w:2.0w $ 、 $ 1.5w:2.4w $ ；三种气泡方案的质量和近似；两气泡初始位置分别设为（ ${x}_{1c}$ ， ${y}_{1c} $ ） =（ ${L}_{x}/\mathrm{2} $ ，100）和（ ${x}_{2c} $ ， ${y}_{2c} $ ）=（ $ {L}_{x}/\mathrm{2}$ ，40）。

图6（a）展示了 $ {d}_{1}:{d}_{2} = 2.4w:1.5w $ ，即上方气泡大于下方气泡时气泡穿过孔板的运动过程。如图所示，在 $ {t}^{*} = 9.99 $ 时下方气泡接触到上方气泡，此时上方气泡部分已经穿过孔板通道，下方气泡与上方气泡尾翼接触并将一个直径约为 $ 1.0w $ 宽的液泡包裹在内。液泡随气泡的上升与其一起穿过孔道。当液泡处于孔板通道内时气泡被夹断分成上下两个气泡（ ${t}^{*} = 11.56$ ），之后上气泡虽然已经穿过孔板通道，但由于形变较严重，维持了较低的运动速度。当液泡受到下气泡对其的推动力后从孔口被排出，下气泡加速向上运动与上气泡合并（ $ {t}^{*} = 17.35 $ ），之后尾部再次破裂产生一个小气泡，并在表面张力的作用下，逐渐恢复为球帽状。

图6（b）给出了上下两气泡直径相等时（ $ {d}_{1}:{d}_{2} = 2.0w:2.0w $ ）气泡的上升过程。如图所示，上气泡靠近孔板后拉伸为长条状，而下气泡也被拉长，上气泡下方和下气泡上方贴合，在中间和两端产生接触点（ $ {t}^{*} = 7.36 $ ），气泡内部形成两个液泡。上气泡与其两端断开，接着液泡从断口排出（ $ {t}^{*} = 11.04 $ ），上气泡端部汇入下气泡的端部。但此时下方气泡质量主要在中部并且已并入上方气泡，其与两端的连接也迅速断开。主气泡加速通过孔板通道并继续上升（ $ {t}^{*} = 15.51 $ ），之后在上升过程中气泡尾部又断裂出一个小气泡。

上方气泡直径小于下方气泡（ $ {d}_{1}:{d}_{2} = 1.5w:2.4w $ ）的结果如图6（c）所示。在 $ {t}^{*} = 6.05 $ 时上方气泡靠近孔板后变扁平状中间凸起。下方气泡从中间和两侧同时接触到上方气泡，形成两个小液泡。与气泡直径相等的情况不同的是，下气泡主要从两端汇入上气泡，再从两端向孔板通道中间收缩（ $ {t}^{*} = 11.04 $ ），过程中两个小液泡被挤压合并为一个直径为 $ w $ 的液泡。液泡在孔板通道下方主要受到气泡对其的表面力，方向从向上转为向中间，液泡几乎处于静止状态，未跟随气泡上升，随后液泡被主气泡从下表面排出。主气泡穿过孔板通道继续上升。

为了研究气泡相对大小对气泡通过孔板通道后剩余气泡质量的影响，增加了几组气泡直径比不同的工况以便比较。图7给出了三组相同变化比例但不同初始质量和的情况，统计了气泡穿过孔板后的剩余质量以及气泡运动到计算区域上方的时刻。从图7中间组可以看出，当上气泡质量大于下气泡质量时，上气泡位于孔板正下方部分质量接触到下气泡后在表面张力作用下向中间收缩，最终有较多的剩余质量可以达到顶部。而当下方气泡大于上方气泡时，两气泡发生合并的时间和位置提前，尽管下方气泡是从两侧汇入上气泡，但两个气泡形变程度较小，合并的气泡质量向孔口方向移动，故最终穿过孔板通道的质量比前者情况少。从中间组可以看出上下气泡直径对比度越大，剩余质量越多；而当上下两个气泡大小相等，穿过孔板到达通道上半部分的剩余质量最少，为0.82。在相反的直径比例下，质量大的气泡处在上方的情况最终的剩余质量始终比处在下方的情况多。另外，随着下方气泡相对上方气泡直径增大，气泡通过孔板到达顶部所需时间变短。结合图7上下两组可以看出气泡尺寸相对孔径越小，其通过得越完整，即剩余质量越多；而较小的质量也导致其受到的浮力相对较小，到达顶部时间变慢，在不同直径比下的剩余质量和到达时间也与中间组的变化趋势相似。

4.3 两气泡之间的距离与气泡与孔板的距离的影响

本小节研究两气泡之间的距离以及气泡与孔板的距离对双气泡运动行为的影响。首先研究不同气泡之间距离的情况，为了方便比较，此算例也以孔板间距 $ w $ 作为参照尺寸， $ {l}_{1} $ 和 $ {l}_{2} $ 都表示为 $ w $ 的倍数。在数值模拟中设置了三组工况，首先使上气泡的位置固定不变，此时上方气泡圆心与孔板间下表面距离 $ {l}_{1} = 3.0w $ 为固定值，接着改变下气泡的位置，使上下气泡圆心之间的距离 $ {l}_{2} $ 分别为2.5 $ w $ 、3.0 $ w $ 和3.5 $ w $ 三种情况。此外， $ Eo = 10 $ ， $ {{d}_{1} = d}_{2} = 20 $ 。

当两个气泡圆心距离 $ {l}_{2} = 2.5w $ 时，所得到的结果如图8（a）所示。上气泡下表面和下气泡上表面的最近点距离为 $ 0.5w $ ， $ {t}^{\mathrm{*}} = 3.68 $ 时下气泡向上浮动从中间接触并入上气泡，并在上气泡到达孔板之前完成合并（ ${t}^{*} = 6.04$ ），气泡间未有液泡形成。合并完成后的主气泡横向拉抻，气泡中间部分质量在浮力作用下穿过孔口；在表面张力作用下气泡两侧向中间收缩穿越孔口，最快地穿过了孔板通道（ ${t}^{*} = 12.35$ ），穿过后恢复为球帽形并保持匀速上升运动。

当两个气泡圆心距离 $ {l}_{2} = 3.0w $ 时，上气泡下表面和下气泡上表面的最近距离为 $ 1.0w $ ，此时下气泡并入上方气泡的同时上方气泡准备穿过孔板通道；当 $ {t}^{*} = 8.41 $ 时下气泡先后从中间和两端接触到了上方气泡，中间夹着液泡。液泡把上气泡与其两端挤断，接着液泡从断口排出（ $ {t}^{*} = 10.77 $ ），上方气泡两端与下方气泡的两侧合并；但此时下方气泡多数质量集中在中间并已并入上方气泡。下方气泡与两侧的连接断开，继续上升。

图8（c）给出了两个气泡圆心距离 $ {l}_{2} = 3.5w $ 时的结果。从图中可以看出， $ {t}^{*} = 9.20 $ 时刻在两个气泡合并前，上气泡中间部分已经与两侧脱离后先行穿过孔板通道，穿过孔板通道时速度上升，但其所占整体质量的比重小，平均上升速度提升不明显；两侧残余的质量与下气泡合并， $ {t}^{*} = 12.88 $ 时上气泡离开孔板通道的同时下方气泡进入孔板通道，气泡之间存在液泡，上下两气泡运动速度均下降；当下方气泡头部穿过孔板通道并把液泡从孔板通道挤出后，其速度上升并加速和上气泡合并成一个气泡。

经模拟计算发现，当 $ {l}_{1} = 3.0w $ 、 $ {l}_{2} = 3.3w $ 时，几乎在上方气泡与两侧脱离的同时下方气泡与其两侧接触，而上气泡主体已开始加速上升。相比 $ {l}_{2} = 3.25w $ 的情况，两气泡主体接触的位置从孔口下方约 $ 0.5w $ 处跃迁到了距离孔口上方 $ 1.5w $ 的位置；气泡主体第一次接触的时间也从 $ {t}^{*}\approx 6.0 $ 推迟到了 $ {t}^{*}\approx 14.2 $ 。在 $ {l}_{1} $ 从 $ 2.5w $ 增加至 $ 3.5w $ 的过程中，临界的 $ {l}_{2} $ 从 $ 3.1w $ 变化为 $ 3.45w $ 。

图9展示了两个气泡从不同位置开始上升时统计得到的速度变化曲线。从图9可以看出，在气泡未到达障碍物通道前，气泡的上升速度变化趋势基本相同；但当上气泡与孔板的距离 $ {l}_{1} $ 不变时，在 $ {l}_{2} = 2.5w $ 气泡间距较小的情况下，初始运动的速度提升最大。这是因为两气泡在通过孔板通道前完成了合并，合并后的气泡受到的浮力增大；在完成合并后通过孔板通道时速度变化剧烈，作为整体通过孔道时减速更明显，通过孔道后的加速也越剧烈；之后因表面张力作用与孔板下方质量发生断裂时, 速度略微下降, 最后再度趋于平稳。而当两气泡间距增大（如 $ {l}_{2} = 3.0w $ ），两气泡在孔板下方发生接触但未完全合并，这也使其速度变化与间距较小时相同，但变化幅度更平和。如果两气泡间距进一步增大，则上气泡部分质量先进入孔道，另外部分质量并入下气泡，穿越孔道的过程分成两次，加速过程随之变为两次；其中第二次在 $ {t}^{*} = 17.35 $ 时，出现转折是在下气泡通过孔道后接触到上气泡并汇入的时刻。

接下来研究上方气泡与孔板间的距离对双气泡运动行为的影响，这里固定两气泡间距，通过改变上气泡的位置来改变气泡与孔板的间距 $ {l}_{1} $ 。图9的实线和虚线分别表示 $ {l}_{1} $ 为 $ 3.0w $ 和 $ 3.5w $ 的情况。从图中虚线和实线的对比可以发现，气泡离孔板越远（ $ {l}_{1} $ 越大），气泡的运动相对越迟滞；初始时表现为速度上升得更高，之后在穿越孔道前后时的速度变化规律相同；气泡穿过孔板通道后又恢复到匀速运动。

为了分析气泡初始位置和初始间距对气泡通过孔板通道后剩余气泡质量的影响，本文在上述三种工况的基础上，通过改变 $ {l}_{1} $ 的大小增加了几组工况。在不同 $ {l}_{1} $ 和 $ {l}_{2} $ 的情况下，气泡穿过孔板后到达计算区域上方的剩余质量进行了统计和比较，如图10所示。在上气泡位置不变的情况下，下气泡所处位置越高（气泡之间间距 $ {l}_{2} $ 越短），两气泡在越早的时刻合并，可以穿过孔板到达顶部的质量也越多。在气泡间距相同的情况下，上气泡与孔板的距离 $ {l}_{1} $ 越大，通过的剩余质量反而越少。对于剩余质量，气泡间的距离比气泡与孔板间的距离更重要。结合图10可以发现终端速度和剩余质量也存在着联系：剩余质量越多，气泡受到的浮力越大，其终端速度也越大。

5 结　论

本文基于大密度比两相流LBM模型，对双气泡在含孔板微通道内的运动过程进行了研究。分别考虑了 $ Eo $ 数、气泡之间的相对大小和相对位置等因素对气泡在孔板通道内流动的影响，得到如下结论：

（1）随着 $ Eo $ 数的增大，气泡在通过孔板通道时所包裹的液泡越多，气泡容易发生变形、合并和破裂；另外当气泡在相同的初始位置上升时，随着 $ Eo $ 数的增大，气泡剩余质量 $ {M}_{s} $ 呈负指数级地减小。

（2）大气泡在上、小气泡在下时，运动过程中气泡受到的阻力更大，表现为上方气泡通过孔板的过程中形变更明显，两个气泡在孔口附近合并得更充分，通道内的残留质量较少；但到达通道顶部需要的时间更长。大气泡在下、小气泡在上时，气泡到达孔板附近的时间更快，此种情况两气泡的形变最小，气泡到达顶部所需要的时间最短。而当上下两气泡相等时，穿过孔板到达通道上半部分的剩余质量相对于其他两种情况最少。

（3）当两个气泡初始距离相近时，两气泡之间的合并通常发生在上方气泡到达孔板通道之前，此时通过的剩余质量最多；随着两个气泡之间的距离拉大，两气泡之间的合并更容易发生在上方气泡经过孔板通道时，通过的剩余质量次之；当两气泡之间的距离增加到一定程度，上方气泡穿过孔板通道并发生断裂之后，下方气泡才追上并与之合并，到达计算区域上方的剩余质量最小。另外，上下气泡初始位置同时升高也可以提高通过孔道的剩余质量。