气泡流不仅存在于日常生活中, 也拥有广泛的工业应用背景,例如二氧化碳的采集和储存、可燃冰的开采、石油的炼化以及核反应的堆冷却等.以上所有工业过程的准确实施,均以了解认识黏性液体中气泡的物理机制为前提,其中一个重要的问题即是气泡与外部边界的相互作用问题^[1-7].

过去数十年来,学者对于气泡流动力学的研究主要集中在开放流场(即无界流域)中的单泡动力学问题,这是泡状流问题的一个最基础构型.对单泡问题的研究主要在形状、速度、路径和尾流结构等方面,前人已经开展了大量的实验^[8-11]或数值^[12-14]研究,代表性的综述可参见文献^[15,16]. 总的来说,气泡的动力学特性由表面张力、黏性、惯性、气泡表面杂质(如果存在)等因素共同决定.

根据$\pi $定理可知,开放流场中的单泡动力学问题可以用4个无量纲参数描述^[17-18], 分别是艾特华数(Eo= $\rho $$_{\rm l}$gd$^{2}$/$\sigma $)、伽利略数Ga= $\rho$$_{\rm l}g^{1 / 2}d^{3 / 2}$/$\mu $$_{l})$、密度比($\rho $$_{\rm r}$= $\rho $$_{l}$/ $\rho $$_{\rm g})$和黏度比( $\mu $$_{\rm r}$= $\mu $$_{l}$/ $\mu $$_{\rm g})$.艾特华数代表重力与表面张力的比值,而伽利略数表示重力与黏性力的比值. 从定义上, 可以看出, Eo和 Ga分别是韦伯数( We = $\rho $$_{l}U^2_{\rm t} d$/ $\sigma $)和雷诺数( Re = $\rho$$_{\rm l}gU_{\rm t}$/$\mu _{\rm l})$的变种. 不过,前面两者是以$\sqrt {gd} $为特征速度,而后面两者是以气泡终态速度$U_{\rm t}$为特征速度. 因此, 更严谨地讲, We和 Re是给定( Eo, Ga, $\rho $$_{r}$, $\mu $$_{\rm r})$条件下的后验参数或结果参数,而非先验参数或控制参数. 此外, 也可用莫顿数( Mo = Eo$^{3}$/ Ga$^{4})$来描述. 对于给定液体而言, Mo为常数.

近年来, 越来越多的学者们开始关注气泡与表面的耦合作用. 此种情况下,气泡的动能会转移到表面能上, 当满足一定条件时,气泡会出现回弹(弹跳)运动.

Tsao等^[19]在实验中观测到气泡会在黏性耗散动能之前在水平壁面处完成数次回弹运动.Klaseboer等^[20]进行了液滴对水平壁面冲击的实验研究并提出了一个考虑了液膜排液效应的回弹轨迹模型.而后,Legemdre研究小组^[21-22]开展了一系列对气泡和液滴在固壁面附近回弹运动的实验研究,他们提出了回复系数($C_{\rm r}$, the coefficent ofrestitution)的概念,用来量化气泡碰撞前后的速度变化.Zawala的研究小组^[23-25]先后研究了气泡与气液交界面和固液交界面的碰撞回弹运动.他们发现, 在固液交界面的情况下,碰撞过程中的能量耗散总量远小于气液交界面的情况.这是因为气液表面可以产生形变.Fujasová-Zedníková等^[26]观测发现,对于固体边界的情况,气泡回弹运动特性与固壁面的材料特性无关,而表面活性剂(杂质)显著影响气泡回弹运动和变形.邱超等^[2]研究发现, 黏度相同时,气泡碰撞过程中的最大变形程度随雷诺数(Re)增大而增大.Pelletier等^[27]通过测量两种初始条件下自来水中气泡撞击水平壁面的过程发现,气泡长宽比是影响回弹过程的关键参数.

以上研究均是通过实验手段完成的. 然而, 由于条件的局限和限制,时至今日用实验手段探究气泡动力学问题仍极具挑战性.随着建模方法、数值算法以及计算能力的高速发展,数值模拟成为研究这一问题的一种极具前景的研究手段^[28-29]. 目前, 这方面的文献还相对较少.

Canot等^[30]在势流理论中加入润滑理论,模拟了二维气泡的回弹运动. 李帅等^[1]基于势流理论, 将Canot等^[30]的二维气泡模型推广到了轴对称的情况,研究了气泡与水平壁面碰撞的平衡构型,并指出韦伯数是影响气泡最终稳态构型的重要影响因素.Omori等^[31]采用界面追踪法(fronttracking)分析了气泡与壁面碰撞过程中的能量耗散过程. Qin等^[32]采用任意拉格朗日--欧拉方法, 研究了大莫顿数(0.1$\le$Mo $\le$16)条件下(即液体黏性很高时)气泡与水平固壁面相互作用的情况.根据气泡进入稳态后的形状的不同, 耦合类型被分成了三类.在他们的研究中, 流场是以有限元的模式计算的, 忽略了泡内的流体流动.气泡表面的捕捉依赖于动网格,液膜排液是通过自适应网格来确保壁面和气泡表面之间至少有3个单元.这样就不能模拟出液膜撕裂的过程,液膜厚度低于气泡半径的1/100时计算便无法继续.Zawala等^[33]成功应用流体体积法(volume of fluid,VOF)研究了顶部壁面做正弦振动时的气泡碰撞动力学问题,侧重分析了液膜的稳定性. 他们发现,液膜稳定性不仅依赖于壁面亲疏水特性,也与多相流系统中的能量传递有关.Albadawi等^[34]评估了VOF方法对气泡碰撞壁面问题的适用性.他们指出了精细地捕捉相间表面的液膜的重要性,特别是对于具有强疏水性能表面的材料而言. 然而,他们采用的是局部加密的静止网格,所以对大接触角时的三相接触线形成的问题对网格数量和计算量要求特别高.最近,Denner^[35]采用VOF法研究了蠕动流条件下的气泡与壁面的碰撞运动,揭示了壁面与气泡碰撞对于毛细管数 Ca的自相似特性.

时至今日, 人们对气泡在水平表面处的回弹动力学的认识还很初步. 例如,已有的实验一般是从有量纲参数出发, 研究水中气泡等少数的情况,而数值研究方兴未艾^[36-39]. 虽然,李帅等^[1]研究了浮力和表面张力对气泡回弹特性的影响,但目前为止, 尚未见到对气泡惯性影响的报道. 此外, 现实情况中,气泡在与顶部壁面碰撞或接触之前并不一定进入了终态,即壁面位置并非处于理想高度,所以研究气泡接触速度的影响也极具现实意义. 有鉴于此,本工作拟采用计算流体力学(computational fluid dynamis, CFD)方法,从无量纲的视角出发, 研究伽利略数 Ga和接触速度$U_{\rm a}$对气泡在壁面回弹的影响, 力求得到更具普遍性的结论.

针对气液两相流问题, 假设流动为不可压的层流流动, 流体等热且不相融.不考虑由于流体剪切而使黏性发生改变, 即关注的流体是牛顿的.为描述这样的两相流问题, 控制方程是密度可变的Navier-Stokes方程

(1)$$\delta\cdot {{ u}} = 0 $$

(2)$$\rho (\frac{\partial }{\partial t}{{ u}} + {{ u}} \cdot \nabla {{ u}})\mbox{ = } - \nabla p + \nabla (2\mu {D} ) + \rho {{ g}}\mbox{ + }\sigma \kappa \delta _s {{ n}} $$

其中, ${ u}$是速度矢量, $\rho $是密度, $t$是时间,$p$是压力, $\mu $是动力黏度, ${{ D}} = {\left( {\nabla {{ u}} + \nabla {{ u}}^{\rm T}} \right)} /2$是变形张量,${ g}$是重力加速度. 表面张力通过动量方程中右边最后一项加入. $\delta $$_{\rm s}$是狄拉克分布函数,表征表面张力集中在两相的表面上. $\sigma $, $\kappa $,${ n}$分别是相间表面张力系数、曲率以及表面法向单位向量.

气液两相流中的物性参数是根据体积加权平均获得的

(3)$$\rho = c\rho _{\rm l} + (1 - c)\rho _{\rm g} $$

(4)$$\mu =c\mu _{\rm l} + (1 - c)\mu _{\rm g} $$

其中, 下标l和g分别表征液体和气体.通过引入体积分数$c$和求解体积分数输运方程来记录气液表面

(5)$$\frac{\partial c}{\partial t} + \nabla \cdot (c{{ u}}) =0 $$

对于CFD而言, 由于气泡回弹过程涉及快速的局部拓扑变化,这对网格分辨率以及计算量要求很高,如使用均匀的静态网格无疑耗费巨大的计算量.本文采用动态自适应网格加密技术(adaptive mesh refiment, AMR).采用结构的笛卡尔网格, 空间域的离散采用分级的八叉树方法^[40].二维问题则采用四叉树方法. 如图1所示, 初始时,计算域都是由正方形构成. 每个原始正方形可以再被细分(加密)成4个子块.根据设定的加密标准, 每个原始正方形可以不断细分. 针对气泡流问题,考虑到体积分数和尾流对流场的决定性影响,同时设置了如下两种加密 标准

(6)$$\frac{\left| {\nabla \times {{ u}}} \right|\Delta }{\max(\left| {{ u}} \right|)} < \varepsilon $$

(7)$$\left| {\nabla c} \right|\Delta < \varepsilon\quad~ $$

其中, $\varepsilon $是标准阈值, 本文统一设定为0.01.以上两种标准的设定可以确保在气泡表面和大涡量处有很高的网格分辨率.自适应的网格加密方法为捕捉流场细节提供了更多的灵活性,平衡了计算量与精度的需求.

时间域的离散上采用的是经典的时间分裂投影算法^[41],得到如下的离散方程组

(8)$$\rho _{n + 1 /2} (\frac{{{ u}}_\ast - {{ u}}_n }{\Delta t} + {{ u}}_{n + 1 /2} \cdot \nabla {{ u}}_{n + 1 /2} - {{ g}})\mbox{ = } \nabla [ {2\mu _{n + 1 /2} ({{ D}}_n + {{ D}}_\ast )} ]\mbox{ + }(\sigma \kappa \delta _s {{ n}})_{n + 1 /2} $$

(9)$$\frac{c_{n + 1 /2} - c_{n - 1 /2} }{\partial t} + \nabla \cdot (c_n {{ u}}_n ) = 0 $$

(10)$${{ u}}_{n{ + }1} = {{ u}}_{\ast } - \frac{\Delta t}{\rho _{n + 1 /2} }\nabla p_{n + 1 / 2}$$

(11)$$\nabla \cdot {{ u}}_{n + 1 /2} \mbox{ = }0 $$

要获得$n$+1时间步的校正后的速度场, 需要先求解泊松方程

(12)$$\nabla \cdot \left( {\frac{\Delta t}{\rho _{n + 1 /2} }\nabla p_{n + 1/2} } \right) = \nabla \cdot {{ u}}_{\ast } $$

然后, 动量方程(8)可写成如下的Helmholtz形式

(13)$$\frac{\rho _{n + 1 /2} {{ u}}_\ast }{\Delta t} - \nabla \cdot \left( {\mu _{n + 1 /2} {{ u}}_\ast } \right) = \nabla \cdot \left( {\mu _{n + 1 /2} {{ u}}_\ast } \right) + (\sigma \kappa \delta _s {{ n}})_{n + 1 /2} + \rho _{n + 1 /2} \left( {\frac{{{ u}}_n }{\Delta t} - {{ u}}_{n + 1 /2} \cdot \nabla {{ u}}_{n + 1 /2} + {{ g}}} \right) $$

上式中,黏性项的Crank-Nicholson离散是无条件稳定的并在时间域上具有标准的二阶精度.速度对流项的离散采用Bell-Colella-Glaz不分裂迎风格式^[42],当库朗数$f_{\rm c}<1$时满足稳定性要求.为精细捕捉气泡变形即气泡周围流场变化, 往往需要大量的时间步.本文取$f_{\rm c}$= 0.3, 时间步由下式控制

(14)$$f_{\rm c} =\frac{\Delta t}{\min \left(\dfrac{\varDelta }{\left| {{{ u}}_i} \right|},\dfrac{\rho \Delta ^2}{\mu },\sqrt {\dfrac{\rho \Delta^3}{\sigma }} \right)} $$

其中, $\varDelta $是局部网格尺寸, ${ u}_{i}$是单元中心或面心在$i$方向的速度.

体积分数输运方程的求解采用的是分段线性几何VOF方法,可分为以下两步：

(1) 相间表面重构;

(2) 几何通量计算和相间表面对流.

在步骤(1)中, 计算单元里的相间表面采用下式

(15)$${{ m}} \cdot {{ x}} = \alpha $$

其中, ${ m}$是相间表面法向量, ${ x}$是位置向量.给定$c$和${ m}$,通过确定计算单元里在平面之下流体体积就可以得到 $\alpha $的值.法向向量${ m}$可以通过近似相邻单元里的体积分数来获得.

完成相间表面重构后,笛卡尔网格上的方向分裂集合通量也可轻松得到^[43].而后得到的二阶精度的对流项可确保形成锋利的相间表面.

对于表面张力的计算是基于平衡力表面张力离散格式^[43],该格式纳入了高度函数(height-function)曲率计算. 经此处理,曲率半径与局部网格尺寸相当,可以得到具有二阶精度和数值平衡的相间表面构型. 这种表面张力的处理,可有效避免VOF模拟中常常遇到的寄生流(parasitic currents)问题.

本工作的建模和模拟都是在基于有限体积法的开源代码Gerris^[40,43]中实现的, 并在Dell Precision Tower 7910工作站上完成并行计算的.

Sharaf等^{[9, 17]}和Cano-Lozana等^[44-45]两个 研究小组对单泡动力学的系统性研究表明, 对于Eo和 Ga的值均较小的气泡而言, 其形状是轴对称的,上升轨迹也是垂直向上的. 本工作将分析限制在轴对称条件下, 即讨论Eo和Ga的值均较小的孤立浮升气泡与顶部水平固壁面碰撞、回弹的动力学过程.

如图2所示, 本工作的计算在一个轴对称域中进行. 除轴边界外,其余边界均采用无滑移边界条件. 为避免侧边面的影响,计算域宽度即圆柱形容器半径取4^[3]. 为壁面底部边界的影响,至少在底部壁面以上2.5$d$位置静止释放, 气泡距离顶部壁面距离为$L$.为保证气泡在与顶部壁面碰撞之前可到达末态, 计算域高度设置为48$d$.初始时刻, 流场静止, 气泡以球形形状在轴线上静止释放.下文中的所有计算结果均以无量纲的形式呈现, 分别以$d$, $\sqrt {d /g} $, $\sqrt {gd} $为长度、时间和速度特征量.

本章将通过模拟自由上升的浮力气泡的浮升过程,来验证本文所采用数值方法的可靠性和准确性. 所谓自由气泡,即气泡不受边界影响或处于开放流域中的情况. 本章中,密度比和黏度比分别固定为1000和100.

气液交界面和气泡尾流都是影响气泡动力学的关键因素, 均需精确捕捉.为此, 本节设计了一个网格收敛性测试,采用了5组不同的网格来进行收敛性分析. 条件参数 Eo = 116, Ga = 23.07. 根据表1的结果可知, 采用网格B时计算收敛. 也即是说,在气泡表面附近, 每单位长度(气泡直径$d)$有128个计算单元;而在大涡量处, 每单位长度有32个计算单元.而在远离气泡表面或者涡量很小的区域,网格密度为每单位长度16个计算单元.

首先, 本文模拟了开放的静止流场中单气泡浮升运动的过程,并将得到的终态雷诺数和气泡构型与前人的实验结果相比较.

表2给出了选取了13组工况的条件参数, 全部来自Bhaga等^[8]的实验.图3比较了本文计算得到的雷诺数Re和前人实验测量的结果$Re_{\exp}$可以看到,虽然本文模拟了很多种工况, 覆盖了大范围的条件参数,但是每组模拟结果与实验结果均吻合良好, 相对误差在3%以内.

图4比较了气泡进入终态后的稳定构型,其中实线曲线表示模拟得到的气泡轮廓. 可以看到,本文的模拟结果与实验观测的气泡形状吻合良好.

注意到, 表2中, 工况3$\sim $5, 7, 9$\sim $11其艾特华数 Eo相同, 而伽利略数 Ga不同. 本节将分析Ga的变化对自由气泡的浮升动力学动力学特性.

如图5所示, Ga改变, 最后得到的终态速度$U_{\rm t}$也不同.总的来说, Ga越大, $U_{\rm t}$也越大.这是由于黏性力的影响减弱, 而惯性力增大引起的. 与此同时,气泡在进入终态之前, 瞬时上升速度的波动在大 Ga时更加明显.当 Ga = 23.07时, 气泡从释放后逐渐加速到最大速度,再呈现微弱的减速到达其终态速度. 而当 Ga = 62.39时,气泡在大约$t$= 1.2时刻到达最大速度后, 会呈现剧烈的波动.气泡速度呈现下降、上升、下降趋势, 在$t$=4.8时刻左右再次得到极大值后, 气泡速度逐渐平稳下来. 相反, Ga$<$10时, 气泡速度变化过程平缓, 从初始释放后, 逐渐加速到终态速度.

对于不同的 Ga的气泡而言, 气泡浮力固定, 黏性阻力随Ga增大而减小, 所以气泡在释放后获得的加速度会随着Ga增大而增大, 这样 Ga更大时可以获得更大的终态速度.

以上现象可以根据气泡的形状衍化过程解释. 从图3可知, Ga增大,气泡的形状变形更为明显. 为表征气泡的形状特性, 本文采用长宽比$A_{\rm r}$(气泡径向宽度$d_{\rm w}$与轴向高度$d_{h}$之比)来定量地分析气泡的形状.图6给出了气泡进入终态后的长宽比与伽利略数的关系. 可以看到,随着 Ga的增大, 气泡长宽比变小, 气泡扁平程度增大. 此外,根据前面的定义可知, $Re = {U_{\rm t} Ga} /{\sqrt {\mbox{g}d} }$,而从图5可知 Ga增大使得$U_{\rm t}$也增大, $\sqrt {\mbox{g}d}$为定值1, 这样可以解释 Ga增大使得 Re增大的原因.

根据得到的气泡终态速度, 可通过下式$C_{\rm d} = {4\left( {\rho_{\rm l} - \rho _{\rm g} } \right)\left| g \right|d} /{\left({3\rho _{\rm l} U_{\rm t}^2 }\right)}$得到气泡进入终态后的阻力系数^[46]. 结果如图6所示,伽利略数越大, 气泡所受阻力系数更小.

前一章研究了开放流域中自由上升的气泡的动力学特性,本章将讨论顶部边界为固体壁面的情况.

首先, 对比Kosior等^[23]的实验结果,进一步验证代码的准确性和适用性. 条件参数为： Ga = 178.33, Eo = 0.2984, $\rho $$_{\rm r}$= 769, $\mu$$_{\rm r}$= 42.74, 顶部壁面平衡接触角 $\theta $$_{\rm c}$= 0.

图7和图8分别给出了气泡瞬时速度和位置的时历曲线. 注意$t$=0时刻代表气泡开始受到壁面影响的时刻(根据速度判定). 可以看到,模拟结果尤其是速度与实验数据非常接近,而气泡质心位置的轨迹也与实验结果基本吻合.误差主要体现在气泡的回弹间隔上,即模拟中回弹比实验中发生地稍早一些. 换言之,模拟中墙壁对气泡的阻尼更大一些,所以气泡与壁面相撞后的回弹距离较实验值稍小一些,但该差异随着回弹次数的增加而逐渐减小.

Albadawi等^[34]的研究表明,气泡与固壁面之间处的网格分辨率对于准确模拟气泡与壁面的碰撞至关重要.因此,本文也对近壁面区域(距壁面的2$d$而距轴线也是2$d$的矩形区域)进行均匀的网格加密,以期缩小与实验数据的误差. 结果表明,即使近壁面区域加密到每单位长度256个计算单元时,与原有结果也无明显变化. 因此,作者认为本文模拟与实验数据之间的误差是由于数学模型、实验误差等不可避免的一些因素造成,是可接受的.

图9给出了回弹过程中的长宽比 Ar的时历曲线. 同图7和图8一样,不论是数值上还是趋势上, 模拟结果与实验结果吻合度很高. 整体而言,依然可以清晰辨认出至少5次明显的回弹过程. 在气泡靠近壁面的过程中,气泡在轴向上被不断挤压(参见图10),泡内气体向径向扩散使得径向宽度增加而轴向高度减小,导致长宽比急剧增大. 在气泡逐渐脱离壁面的过程中,气泡形状呈相反的变化趋势, 致使$A_{\rm r}$迅速减小.

在气泡向下回弹的末端, 气泡长宽比时历曲线会经历几次迅速的波动,这比图7中速度轨迹的波动更为剧烈. 这主要是因为, 在某些时刻,气泡内部轴向速度的方向并不完全一致造成的,而气泡速度$U$的是由气相的轴向速度的积分平均而定义的. 相似的是,不论是图7或是图9中, 这种波动随着回弹次数的增多而逐渐消失.

图10给出了气泡在第一次与墙壁碰撞并回弹的发展过程,其中粗直线表示顶部壁面. 从气泡的瞬时形状可见,数值模拟很好的重构了实验中气泡的变形过程.

为分析气泡碰撞和回弹过程中的涡结构变化,图11就图10中4个代表时刻给出了气泡周围流场的速度矢量图,其中粗直线表示顶部固壁面. 从中可以看到,在气泡靠近并碰撞壁面的阶段(即从$t$= 0.1到$t$= 0.4),气泡径向边缘的涡环$V_{1}$逐渐由气泡的中部位置向顶部位置转移,同时$V_{1}$强度逐渐减弱、涡环半径逐渐增大,这使得诱导液体产生的向上的射流强度减弱, 气泡的尾流效应减弱,流场中的循环流动逐渐减弱, 气泡随之渐渐减速. 当$t$= 0.5时,受壁面挤压的作用, 气泡边缘处的$V_{1}$径向上向内侧偏移,即此时$V_{1}$完全处于气泡内. 这种气泡内旋涡裹挟着气体的现象,作者曾在之前对于裙带气泡在极度受到径向壁面约束的情况下观察到^[3].不过, 这里旋涡$V_{1}$的强度已经很弱. 当$t $= 0.6时,气泡顶部即将离开壁面, 呈桃形. 此时,$V_{1}$转移到这个``桃子''的顶端最终消失($t $= 0.7时已不能捕捉) .有趣的是, 与此同时, 在气泡的径向表面周围,一个新的旋涡$V_{2}$逐渐形成.它的产生是由于壁面附近液体径向上向轴线方向冲击气泡并与泡内旋涡相互作用形成的.需要注意的是, $V_{2}$的方向与$V_{1}$相反, 随着气泡远离壁面,$V_{2}$的强度增大, 诱导出向下的液体射流, 气泡开始做向下的加速运动.

此外, 可以注意到, 图11中靠近壁面的旋涡的强度均较小,这即是之前对于壁面附近网格加密而结果并无改变的原因.

本节将研究不同 Ga时气泡在顶部壁面处的回弹动力学特性,分析气泡惯性在回弹动力学中的影响. Eo, $\rho $$_{\rm r}$, $\mu $$_{\rm r}$和 $\theta $$_{\rm c}$均与上节参考算例中的参数保持一致, Ga取70, 110,150和178.33(参考算例).

图12给出了不同伽利略数时的回弹高度$H$和回弹周期$T$的结果(参数定义示意图见图8).结果显示, Ga对这几个回弹参数的影响是单调递增的,即$H$和$T$均随着 Ga的增大而增大, 增长规律是近似线性的.

可以从气泡能量转换的角度解释$H$随 Ga增大的原因.以第一次碰撞为例, Ga增大,气泡碰撞后变形程度更大(从图13可知), 附加质量更大. 而由图14可知, Ga更大时, 气泡的质心位置极大值更大,即碰撞时重力势能更小(以底部边界为参考线). 需要强调的是,虽然变形更大使得气泡在碰撞时能量的黏性耗散增大,但这里可以不用考虑能量黏性耗散的影响. 这是因为,Zawala等^[24]对能量的定量分析表明,气泡能量的黏性耗散几乎在碰壁时存在且量级远低于重力势能和动能.再考虑到, Ga更大时, 气泡获得的接触速度本就更大,动能也就更大. 自然地,当气泡向下运动到动能为0(将动能转化为势能)所需的距离必然更长.

Ga增大时, 液体黏性减弱, 气泡运动所受到的黏性阻力降低.这样气泡碰撞后的加速度数值上会更小, 先加速再减速的过程会更慢.因此, 气泡回弹周期$T$也更大.

图13对比了不同 Ga时的气泡长宽比时历曲线. 根据定义可知,液体黏性和 Ga成反比. 因此, 在液体黏性强(小 Ga)时,气泡变化所受阻力大, 气泡在与壁面碰撞后形变程度也随之较小,经历小幅度的回弹后气泡即可稳定地贴附于壁面. 这即是 Ga = 70时,只观察到3次明显的回弹且大致在$t$= 4时刻后{$A_{\rm r}$}值就趋于稳定. 相反地, Ga更大时,气泡由于获得更大惯性在碰撞之前获得了更多的动能,受阻尼消耗所需的时间也就更长.

引言中提到, 回复系数($C_{\rm r}$= $U_{\rm b}$/$U_{\rm a})$是描述气泡对壁面回弹运动特性的一个重要的概括性参数. 这里,$U_{\rm a}$是气泡接触壁面前的速度, $U_{\rm b}$是气泡离开壁面后的速度(示意图可见图7). 从图15可知, Ga对回复系数$C_{\rm r}$的影响是正相关的.

现实情况中, 气泡在与顶部壁面碰撞或接触之前并不一定进入了终态.例如, 气泡所在流域的高度可能是受限的,这样气泡与壁面碰撞前并未获得终态上升速度$U_{\rm t}$.本节通过改变初始时刻气泡的释放位置(气泡与壁面的初始间距$L)$,获得不同接触速度$U_{\rm a}$, 以分析对结果的影响.

图16展示了接触速度的变化对回弹高度和周期的影响. 接触速度更大时,即气泡距离顶部水平壁面的初始距离$L$更大,气泡从释放之后的加速过程更长更充分,在接触碰撞壁面之前获得的动能也就更大. 因此, 可以看到,图16中回弹高度和周期均随$U_{\rm a}$增大而增大. 需要注意,与图14中的变化趋势不同的是, 回弹高度随$U_{\rm a}$变化的规律并非线性的. 当$U_{\rm a}$超过2.34之后,两次回弹高度$H_{1}$和$H_{2}$均开始大幅度增大.

图17用长宽比表征了气泡在整个碰撞{--}回弹过程中的形状变化.虽然{$U_{\rm a}$}增大后气泡的变形程度显著增大, 然而,但长宽比时历曲线变化趋势大致相同, 即使{$U_{\rm a}$}从2.6降低到1.47, 仍可观测到4$\sim $5次明显的回弹过程.对比3.2节的情况可知, 相比于伽利略数的影响,气泡的回弹运动特性受接触速度的影响相对较小. 需要注意的是,当气泡稳定附在壁面上(长宽比保持定值)时, 图14和图16均未有明显变化,这是因为气泡接触壁面后的稳定构型主要是由于接触角决定的.

从图18可见, 接触速度的改变并未对回复系数造成太多改变, 这说明$C_{\rm r}$并不依赖于$U_{\rm a}$,这与Zenit等^[22]的实验观测结果一致. 进一步地讲, 这表明,$U_{\rm a}$并不是决定气泡回弹动力学特性的控制参数.

本文基于动态自适应网格技术结合VOF方法, 模拟了直线上升轨迹(the rectilinear regime^[45]的浮力气泡撞击顶部壁面而回弹的运动过程,研究了气泡对壁面的回弹动力学特性.

研究表明, 伽利略数 Ga对浮力气泡与壁面的碰撞过程具有重要影响. Ga的增大, 气泡受到的黏性阻力减小, 惯性增大,可获得更大的接触速度(以及动能).回弹参数(回弹高度和回弹周期)、回复系数均与 Ga正相关, 所以Ga增大可支持气泡与壁面碰撞后发生更持久的回弹,气泡形变程度更为剧烈.

接触速度$U_{\rm a}$的改变在一定程度上可产生与$Ga$类似的影响,即增大$U_{\rm a}$可促进气泡碰撞产生更为剧烈的形变、更大的回弹高度及周期. 然而,回复系数并不受$U_{\rm a}$的影响.

下一步工作中, 将采用完全三维数值模拟,即研究具有气泡不稳定性的情况,以求进一步认识气泡对于顶部壁面的回弹动力学特性.