摘要:对装有水的容器施加垂直振动, 水中的气泡可能出现下沉的现象. 针对球形气泡构建了基于理想气体方程的气泡运动学模型, 该模型考虑了Basset力对气泡运动的影响. 观察到气泡的运动呈振荡形式, 存在一个临界深度, 气泡在临界深度时稳定振荡, 在临界深度以上时上升, 在临界深度以下时下沉. 采用分离气泡运动和引入收敛因子的方法对Basset力进行理论求解, 并通过等步长复合梯形公式对Basset力进行数值求解. 模型的数值模拟结果表明: 附加质量力是气泡下沉的关键因素, 而Basset力对气泡下沉的临界深度和气泡下沉的前期轨迹无明显影响, 但对气泡下沉的后期轨迹有较大影响. 采用去噪、二值化、填充图像等数字图像处理的办法提取气泡的特征尺寸, 可提高实验参数测量的准确性, 使气泡下沉的临界深度与理论值相匹配.
关键词: 垂直振动/
气泡下沉/
稳定振荡/
Basset力

English Abstract

全文HTML

--> --> -->

当水柱在外界激励下做垂直振动时, 水中气泡可能不会上浮, 而是出现下沉的现象. 气泡下沉的现象主要分为两种: 当振动产生的压强远小于大气压强时, 气泡以相对分散的形式下沉; 当振动产生的压强远大于大气压强时, 气泡以聚集成团的形式下沉^[1,2]. 文献[3, 4]给出了气泡运动的微分方程, 认为在非饱和流体中气泡下沉是由于气泡体积变化导致附加质量变化, 从而产生一个累积效应向下的力所引起的, 并发现存在一个临界深度, 使得气泡在临界深度以下下沉, 在临界深度以上上浮. 文献[5-7]从共振的角度研究了振荡流体中气泡悬浮的现象, 并通过实验验证了气泡悬浮的位置. 文献[8-11]认为在饱和流体中气泡下沉是由流体中驻波推动一个比共振尺寸小的气泡作高速直线运动所引起的. 文献[12, 13]研究了振荡流体中的固体颗粒, 发现在流体的驻波场中固体颗粒会上升. 文献[14, 15]对垂直振动下气泡产生的机理进行研究, 指出当振动加速度大于临界值时气泡将分裂成更小的部分, 并最终在形态和数量上趋于稳定. 文献[16-19]讨论了气液传质条件和Basset力对上浮气泡运动学行为的影响, 发现Basset力和传质条件对于易溶气泡有明显影响. 当振动产生的瞬时压强远大于大气压时, 会出现空化的现象, 空化的程度取决于气泡形状、大小、数量等因素^[20-22]. 因此对于非饱和流体中气泡下沉、Sorokin临界深度进行验证、以及Basset力对气泡下沉的影响是研究气泡下沉的热点问题. 本文通过构建非饱和流体中气泡的运动方程, 阐明气泡下沉的物理本质, 并研究Basset力的理论求解及Basset力对振荡非饱和流体中气泡运动的影响, 希望为火箭、船舶箱体中气泡的去除提供理论基础.

2

-->

2.1.模型准备

以静止时的液面所处的位置为0点, 竖直向下为正建立坐标系, 以x表示气泡和液面之间的距离, 在正弦函数的激励下, 气泡在x处的内压强可以表示为^[4]

${P_x} = {P_0} + \rho gx + \frac{{2\sigma }}{R} + \rho xA{\omega ^2}\sin \omega t,$

式中, ${P_0}$为大气压强为, $\rho $为水的密度, g为重力加速度, $\sigma $为液体表面张力系数, R为气泡半径, $\omega $为振动的角频率, A为振动的振幅. 考虑到气泡内气体的温度和压强都不太高, 且所观测的气泡半径R在0.5—1.5 mm, 忽略表面张力所产生的压强, 将其看作理想气体, 并忽略温度的变化, 满足

${P_0}{V_0} = {P_x}{V_x},$

${V_0}$为气泡的初始体积, 则在x处的气泡体积可以表示为

$\begin{split}{V_x} \;&= \frac{{{P_0}{V_{\rm{0}}}}}{{{P_0} + \rho x(g + A{\omega ^2}{\rm{sin}}\omega t)}} \\\;& = \frac{{{V_{\rm{0}}}}}{{{{1 + }}\dfrac{{\rho x}}{{{P_0}}}(g + A{\omega ^2}{\rm{sin}}\omega t)}}. \end{split}$

当ρ = 998 kg/m³, x = 0.17 m, A = 0.0015 m, P₀ = 101.325 kPa, ω = 100π rad/s, g = 9.8 m/s²时, $\dfrac{{\rho x}}{{{P_0}}}(g + A{\omega ^2}) \approx {\rm{0}}.{\rm{26}}$, 对(3)式进行Taylor展开, 保留一次项, 得到

${V_x} \approx {V_{\rm{0}}}\left[ {{\rm{1}} - \frac{{\rho x}}{{{P_0}}}(g + A{\omega ^2}{\rm{sin}}\omega t)} \right]. $

气泡体积的舍去误差约为3%. 考虑到惯性力、浮力、表面张力和黏性阻力都会对气泡的形状产生影响, 这里采用$ {\rm{Tadaki}} $^[23]提出的无量纲数$Ta$来描述气泡的形态, 其表达式为

$Ta = {\rho _{\rm{l}}}vd{\left[ {\frac{{({\rho _{\rm{l}}} - {\rho _{\rm{g}}})g}}{{{\sigma ^3}\rho _{\rm{l}}^{\rm{2}}}}} \right]^{0.23}} \approx 0.704,$

其中, d为气泡的当量直径, v为气泡运动的相对速度, ${\rho _{\rm{l}}}$表示液相密度, ${\rho _{\rm{g}}}$表示气相密度. 当$Ta$小于1时, 根据刘柳^[24]的实验结果, 气泡的横纵比不低于0.9, 因此可以近似认为气泡呈球形. 将气泡假设为球形, 得到位于x处的气泡半径:

${R_x} = \sqrt[3]{\frac{3}{\rm{4\pi }}{V_x}} = {R_0}{\left[ {\frac{1}{{{\rm{1}} \!+\! \dfrac{{\rho x}}{{{P_0}}}(g \!+\! A{\omega ^2}{\rm{sin}}\omega t)}}} \right]^{\textstyle\frac{{\rm{1}}}{{\rm{3}}}}}. $

对(6)式进行Taylor展开, 保留一次项, 得到

${R_x} \approx {R_{\rm{0}}}\left[ {{\rm{1}} - \frac{{\rho x}}{{3{P_0}}}(g + A{\omega ^2}{\rm{sin}}\omega t)} \right], $

式中, ${R_{\rm{0}}}$为气泡初始半径. 当ρ = 998 kg/m³, g = 9.8 m/s², ω = 100π rad/s, P₀ = 101.325 kPa, x = 0.17 m, A = 0.0015 m时, $\dfrac{{\rho x}}{{{\rm{3}}{P_0}}}(g + A{\omega ^2}) \approx {\rm{0}}.{\rm{09}}$, 气泡半径的舍去误差约为1%.
2 -->

2.2.气泡的运动分析

对于水中变速运动的气泡, 我们主要考虑黏性阻力、浮力、水的推力、附加质量力和Basset力, 忽略气泡自身重力.
1) 黏性阻力
当气泡与周围的水存在相对运动时, 由于气泡表面与水之间存在摩擦, 气泡表面对液相有阻力作用, 反过来气泡也会受到周围水带来的阻力, 气泡黏性阻力公式可以由下式给出^[4]:

${F_{\rm{D}}} = {\rm{4}}\psi \rho {R^{\rm{2}}}{v^{\rm{2}}}{\rm{sgn}} v.$

对于气泡阻力系数的确定, 需要对雷诺数Re进行估计:

$Re = {{2\rho Rv}}/{\mu }, $

其中μ是液体动力黏度系数. 在本实验中, 雷诺数的取值范围为1000—5000, 根据标准阻力系数曲线, 在该范围内阻力系数变化不大, $\psi \approx {\rm{0}}.{\rm{2}}$.
2) 浮力
由于水中压强随深度的增大而增大, 气泡下表面的压强大于上表面的压强, 因此存在一个向上的浮力, 浮力等于气泡排开水的重力:

${F_{\rm{v}}} = \rho g{V_x}.$

3) 水流推力
瓶子中的水在激振器的激励下做正弦函数的上下运动, 产生一个周期性的推力, 其表达式为

${F_{\rm{p}}} = \rho {V_x}\frac{{{{\rm{d}}}{u_{\rm{l}}}}}{{{{\rm{d}}}t}} = \rho {V_x}A{\omega ^{\rm{2}}}\sin \omega t, $

其中${u_{\rm{l}}}$表示水流的速度.
4) 附加质量力
对于水中变速运动的气泡, 推动气泡运动的力除了要推动气泡本身, 还要将周围的水一起推动, 这意味着推动气泡运动的力F就会满足, $ F > m{\ddot x}$, 可以将其写成$ F=({m}_{0}+m){\ddot x}={F}_{\rm{m}}+m{\ddot x}, \;{F}_{\rm{m}}$定义为附加质量力, ${m_{\rm{0}}}$定义为附加质量, 且它们可以表示为

$ {F}_{\rm{m}}={C}_{\rm{m}}\rho {V}_{x}\frac{{\rm{d}}{\dot x}}{{\rm{d}}t},$

$ {m_0} = {C_{\rm{m}}}\rho {V_x}, $

式中${C_{\rm{m}}}$为附加质量系数, 它跟气泡的形状有关^[25]:

${C_{\rm{m}}} = 0.62\chi - 0.12,~~{{ 1 < }}\chi {{ < 2}}.{\rm{5}}, $

其中$\chi $为气泡的横纵比. 考虑将气泡假设为球形, 故$\chi {{ = 1}}, {\;}{C_{\rm{m}}} = 0.5$. 此外, 当附加质量发生变化时, 还会产生一个额外的附加力$\dot m \dot x $^[4]. 将附加质量力和额外的附加力合在一起, 记成:

${F}_{\rm{k}}={m}_{0}{\ddot x}+ {{\dot m}_{0}} {\dot x}.$

附加质量力的意义在于突然加速分散相周围的连续相的表面质量所需要的力, 如果分散相的密度与连续相的密度之比远大于1, 附加质量力的影响可以忽略不计^[26]. 模型中分散相为空气, 其密度远小于连续相为水的密度, 因此附加质量力不能忽略. 可以看到附加质量力的大小为气泡所受浮力的一半.
5) Basset力
Basset力是指当气泡与周围的水存在相对加速度时, 气泡的运动与气泡表面附着层的运动存在滞后, 由此产生一个黏性附加力, 与气泡的加速过程有关, 其表达式采用Basset的研究结果:

$ {F}_{\rm{B}}={C}_{\rm{B}}{R}^{2}\sqrt{{\rm{\pi}}{\rho }_{\rm{l}}\mu }{\displaystyle \int_{t_0}^{t}\frac{{\rm{d}}{\dot x}/{\rm{d}}\tau }{\sqrt{t-\tau }}}{\rm{d}}\tau , $

其中, ${t_0}$是气泡开始加速时刻; ${C_{\rm{B}}}$为Basset力系数, Basset通过理论计算得到${C_{\rm{B}}}{{ = 6}}$.
对于水中的气泡而言, 在加速运动的初期Basset力较为重要, 当运动达到稳定后可以忽略此力. 本实验中气泡以振荡的形式运动, 会出现很多个加速运动的初期, 所以对Basset力进行考虑.
给出气泡的运动方程:

$ m{\ddot x}=-{F}_{\rm{D}}-{F}_{\rm{v}}-{F}_{\rm{p}}-{F}_{\rm{k}}-{F}_{\rm{B}}.$

此外, 对于任意时刻气泡的内压强${P_x}$应满足:

${P_x} > {P_{{\rm{d}}}}, $

其中P_d为水的饱和蒸汽压, 否则会出现空化的现象, (17)式就不再成立.

2

-->

3.1.Basset力

对于Basset力的数值求解, 注意到Basset力是以奇异积分的形式出现, 黄社华和程良骏^[17]证明了Basset力的收敛性, 并通过等步长复合梯形公式将Basset力离散成以下形式:

$\begin{split}{F_{\rm{B}}} =\;& \frac{{\rm{1}}}{{\rm{2}}}Bh\left[ {\frac{{g(0)}}{{\sqrt t }} + 2\sum\limits_{i = 1}^{n - 2} {\frac{{g(ih)}}{{\sqrt {t - ih} }}} + \frac{{g(t - h)}}{{\sqrt h }}} \right] \\&+ 2Bg(t)\sqrt h ,\\[-10pt] \end{split}$

其中, $h \!=\! t/n$是步长, $g(\tau )\!=\!{\rm{d}}{\dot x}/{\rm{d}}\tau$, $B \!=\! 6{R^2}\sqrt {{\rho _l}\mu {\text{π}}}$.
对于Basset力的理论求解, 根据振动力学中直接分离运动的方法, 将气泡的运动x分解成^[4]:

$x(t) = X(t) + \eta (t),$

其中, $X(t)$表示气泡的平均运动, 也称为气泡的慢运动; $\eta (t)$表示气泡的脉动, 也称为气泡的快运动. 将(20)式代入(16)式得到:

$ \begin{split}{F}_{\rm{B}}=\;&6{R}_{\rm{0}}^{{2}}{\left[{1}-\frac{\rho }{3{P}_{0}}(X(t)+\eta (t))(g+A{\omega }^{2}{\rm{sin}}\omega t)\right]}^{\rm{2}}\\&\times\sqrt{{\rm{\pi}}{\rho }_{\rm{l}}\mu }{\displaystyle \int _{t_0}^{t}\frac{{\ddot X}(\tau )+{\ddot \eta}(\tau )}{\sqrt{t-\tau }}}{\rm{d}}\tau .\\[-20pt] \end{split}$

考虑到${\dot \eta}\gg {\dot X}, {\ddot \eta}\gg {\ddot X}, X\gg \eta, {\left[\dfrac{\rho x}{{P}_{0}}(g \!+\! A{\omega }^{2})\right]}^{2} \!\ll\! 1$, 舍去(21)式中的较小量进一步得到:

$ {F}_{\rm{B}}=6{R}_{\rm{0}}^{\rm{2}}\sqrt{{\rm{\pi}}{\rho }_{\rm{l}}\mu }{\displaystyle \int _{{t}_{0}}^{t}\frac{{\ddot \eta}(\tau )}{\sqrt{t-\tau }}}{\rm{d}}\tau . $

求解(22)式的关键在于求解式中的奇异积分项, 将奇异积分项单独提出, 并假设气泡的脉动方程具有以下形式的解:

$\eta = B\sin (\omega t + \varphi ), $

其中$\varphi $为相位, 则

$\begin{split}{\displaystyle \int_{{t}_{0}}^{t}\frac{{\ddot \eta}(\tau )}{\sqrt{t-\tau }}}{\rm{d}}\tau =\;&-B{\omega }^{2}{\displaystyle \int_{{t}_{0}}^{t}\frac{{\rm{sin}}(\omega t+\varphi )}{\sqrt{t-\tau }}}{\rm{d}}\tau \\ {=}\;&-B{\omega }^{2}\sqrt{\frac{2{\rm{\pi}}}{\omega }}\!\left[\!{\rm{cos}}(\omega t\!+\!\varphi )S\!\left(\!\sqrt{\frac{2\omega }{{\rm{\pi}}}t}\!\right)\right. \\&\left.-{\rm{sin}}(\omega t+\varphi )C\left(\sqrt{\frac{2\omega }{{\rm{\pi}}}t}\right)\right],\\[-20pt]\end{split} $

其中$S(\alpha )$和$C(\alpha )$为菲涅耳积分, 且菲涅耳积分无法用初等函数的形式表示:

$S(\alpha ) = \int_0^\alpha {\sin {q^2}{{\rm{d}}}q} , $

$C(\alpha ) = \int_0^\alpha {\cos {q^2}{{\rm{d}}}q} . $

将收敛因子${{\rm{e}}^{ - \beta {q^2}}}$引入(25)和(26)式^[27], 得到:

$G(\beta ) = \int_0^\alpha {{{\rm{e}}^{ - \beta {y^2}}}\sin {y^2}{{\rm{d}}}y} , $

$F(\beta ) = \int_0^\alpha {{{\rm{e}}^{ - \beta {y^2}}}\cos {y^2}{{\rm{d}}}y} . $

则可以通过计算二重积分来间接计算菲涅耳积分:

$\begin{split} {(F(\beta ))^2} - {(G(\beta ))^2} =\;& \int_0^\alpha {\int_0^\alpha {{{\rm{e}}^{ - \beta ({y^2} + {z^2})}}\cos ({y^2} + {z^2})} } {{\rm{d}}}y{{\rm{d}}}z \\ =\;&\left( {\int_0^{\tfrac{{\rm{\pi }}}{2}} {\int_0^\alpha {{{\rm{e}}^{ - \beta {r^2}}}r\cos {r^2}} } {{\rm{d}}}r{{\rm{d}}}\theta + \int_0^{\tfrac{{\rm{\pi }}}{2}} {\int_0^{\sqrt 2 \alpha } {{{\rm{e}}^{ - \beta {r^2}}}r\cos {r^2}} } {{\rm{d}}}r{{\rm{d}}}\theta } \right){{\Big/2}} \\ =\;& \frac{{{\rm{\pi }}{{\rm{e}}^{ - \beta {\alpha ^2}}}\beta }}{{{\rm{8}}({\beta ^2} + 1)}}\left[ {{{\rm{e}}^{\beta {\alpha ^2}}} - \cos {\alpha ^2} + \frac{{\sin {\alpha ^2}}}{\beta } + {{\rm{e}}^{ - \beta {\alpha ^2}}}\Big({{\rm{e}}^{2\beta {\alpha ^2}}} - \cos 2{\alpha ^2} + \frac{{\sin 2{\alpha ^2}}}{\beta }\Big)} \right] = {k_1}(\beta ), \end{split}$

$\begin{split} {\rm{2}}F(\beta )G(\beta ) =\;& \int_0^\alpha {\int_0^\alpha {{{\rm{e}}^{ - \beta ({y^2} + {z^2})}}\sin ({y^2} + {z^2})} } {{\rm{d}}}y{{\rm{d}}}z \\ =\;&\left( {\int_0^{\tfrac{{\rm{\pi }}}{2}} {\int_0^\alpha {{{\rm{e}}^{ - \beta {r^2}}}r\sin {r^2}} } {{\rm{d}}}r{{\rm{d}}}\theta + \int_0^{\tfrac{{\rm{\pi }}}{2}} {\int_0^{\sqrt 2 \alpha } {{{\rm{e}}^{ - \beta {r^2}}}r\sin {r^2}} } {{\rm{d}}}r{{\rm{d}}}\theta } \right)\Big/2 \\=\;&\frac{{{\rm{\pi }}{{\rm{e}}^{ - \beta {\alpha ^2}}}\beta }}{{{\rm{8}}({\beta ^2} + 1)}}\left[ {\frac{{{{\rm{e}}^{\beta {\alpha ^2}}}}}{\beta } - \sin {\alpha ^2} - \frac{{\cos {\alpha ^2}}}{\beta } + {{\rm{e}}^{ - \beta {\alpha ^2}}}\Big(\frac{{{{\rm{e}}^{2\beta {\alpha ^2}}}}}{\beta } - \sin 2{\alpha ^2} - \frac{{\cos 2{\alpha ^2}}}{\beta }\Big)} \right] ={k_2}(\beta ),\end{split} $

$\begin{split}{\left[ {{{(F(\beta ))}^2} + {{(G(\beta ))}^2}} \right]^2} =\;& {\left[ {{{(F(\beta ))}^2} - {{(G(\beta ))}^2}} \right]^2} + 4{(F(\beta )G(\beta ))^2}.\end{split}$

当$ \beta = 0$时,

$\begin{split} G(0) =\;& S(\alpha ) = \sqrt {\frac{{\sqrt {{{({k_1}(0))}^2} + {{({k_2}(0))}^2}} - {k_1}(0)}}{2}} = \sqrt {\frac{{\rm{\pi }}}{{16}}} \sqrt {\sqrt {6 - 2\cos {\alpha ^2} - 4\cos 2{\alpha ^2}} - \sin {\alpha ^2} - \sin 2{\alpha ^2}}, \end{split} $

$\begin{split} F(0) =\;& C(\alpha ) = \sqrt {\frac{{\sqrt {{{({k_1}(0))}^2} + {{({k_2}(0))}^2}} + {k_1}(0)}}{2}} = \sqrt {\frac{{\rm{\pi }}}{{16}}} \sqrt {\sqrt {6 - 2\cos {\alpha ^2} - 4\cos 2{\alpha ^2}} + \sin {\alpha ^2} + \sin 2{\alpha ^2}} .\end{split} $

将(24)式、(32)式、(33)式代入(22)式中即得到垂直振动下气泡所受Basset力的近似解析表达式.
当气泡运动到某一深度时, 气泡会出现稳定振荡的现象, 这一深度被称为临界深度, 当气泡在临界深度以上时, 气泡上升, 在临界深度以下时, 气泡下沉, Sorokin等^[4]给出了不考虑Basset力时气泡的慢运动方程:

$\begin{split}& (m+{m}_{01}){\ddot X}+\frac{1}{T}\!{\displaystyle \int_{0}^{T}\!{F}_{D}}{\rm{d}}t \\={}& \gamma \rho {V}_{0}w\frac{A{\omega }^{2}X}{2{H}_{0}}-(\rho {V}_{0}-m)g, \end{split}$

其中T是周期, ${m_{01}}{{ = 0}}{{.5}}\rho {V_{\rm{0}}}, \gamma = \dfrac{{\rho {H_0}g}}{{{P_0}}}, w = \dfrac{{A{\omega ^2}}}{g}$. 进一步考虑Basset力, 只需在(34)式等号右端加入$- \dfrac{1}{T}\displaystyle \int_0^T {{F_B}{\rm{ d}}t}$这一项, 考虑到附加质量的平均累积效应比Basset力的平均累积效应高出两个数量级, 即:

$\frac{1}{T}\int_0^T {{F_B}{\rm{ d}}t} \ll \frac{1}{T}\int_0^T {{F_m}{\rm{ d}}t},$

因此, 可以忽略Basset力对慢方程的影响, 即仍然采用Sorokin等^[4]提出的临界深度:

${X_0} = \dfrac{{2{H_0}}}{{\gamma {w^2}}}\dfrac{{2\left( {1 + \sqrt {1 + \kappa \dfrac{{{A^2}}}{{R_0^2}}} } \right) + \kappa \dfrac{{{A^2}}}{{R_0^2}}}}{{2\left( {1 + \sqrt {1 + \kappa \dfrac{{{A^2}}}{{R_0^2}}} } \right) + \dfrac{\kappa }{3}\dfrac{{{A^2}}}{{R_0^2}}}},$

其中$\kappa = \dfrac{{{\rm{256}}{\psi ^2}}}{{{{\text{π}}^{\rm{4}}}C_m^4}}$.
2 -->

3.2.模型的数值模拟

为了进一步分析气泡长时间运动后的规律, 以及比较不同初始位置对数值求解的影响, 将求解(17)式转化为求解一阶常微分方程组初值问题, 采用向前差分格式将微分方程组转化为差分方程组, 并进行迭代求解. 选取参数ρ = 998 kg/m³, x = 0.17 m, A = 0.002 m, P₀ = 101325 Pa, ω = 100π rad/s, g = 9.8 m/s², μ = 1.01 × 10^–3 Pa·s, R₀ = 0.001 m, 并给定初始速度$ {\dot x}(0)=0$, 得到图1所示的结果.
图 1 气泡运动方程的数值模拟　(a) 气泡初始位置为0.1 m; (b) 气泡初始位置为0.05 m; (c) 气泡初始位置为0.076 m
Figure1. Numerical simulation of the bubble motion equation$ : $ (a) The initial position of the bubble is 0.1 m; (b) the initial position of the bubble is 0.05 m; (c) the initial position of the bubble is 0.076 m.

对于图1(a), 气泡的初始位置为0.1 m, 气泡从该位置开始运动一直下沉到底部; 对于图1(b), 气泡的初始位置为0.05 m, 气泡从该位置开始运动一直上升至液面; 对于图1(c), 气泡的初始位置为0.076 m, 气泡在该位置附近振荡了一段时间后开始下沉到底部. 数值结果表明, 存在一个临界深度X₀, 使得气泡在临界深度以上时上升, 在临界深度以下时下沉, 在临界深度处稳定振荡. 此外, 从图1(b)可以观察到气泡下沉是一个周期性累积的结果.
将Basset力的数值解(19)式和解析解(22)式进行对比分析, 选取参数R₀ = 0.001 m, ρ_l = 998 kg/m³, μ = 1.01 × 10^–3 Pa·s, ω = 100π rad/s, φ = 0, A = 0.002 m, g = 9.8 m/s², P₀ = 101325 Pa, B = 1/3A, x(0) = 0.1, x(0) = 0进行数值模拟, 结果如图2所示.
图 2 Basset力解析解和数值解的比较
Figure2. Comparison between the analytical and numerical solutions of Basset forces.

本文选取$\varphi {{ = 0}}$, 所以和数值解存在一定的相位差, 但在周期和振幅上两者基本保持一致.
2 -->

3.3.Basset力对气泡运动的影响

为了研究Basset力对气泡运动的影响, 首先采用以下参数对气泡的运动方程进行数值模拟, 分析Basset力对气泡运动的影响, 其中μ =1.01 × 10^–3 Pa·s, R₀ = 0.001 m, g = 9.8 m/s², P₀ = 101325 Pa, ρ_l = 998 kg/m³, ω = 100π rad/s, A = 0.002 m, x(0) = 0.1, x(0) = 0, 结果如图3所示.
图 3 Basset力对气泡运动方程的影响
Figure3. Effect of Basset force on the equation of motion of bubbles.

对两条运动轨迹进行比较, 在第3 s时两者偏移了6%, 在第4 s时两者偏移了20%. 结果表明, Basset力对气泡下沉前期的运动轨迹无明显影响, 但对后期的运动轨迹会产生较大的偏移. 为了进一步分析气泡轨迹偏移的原因, 将气泡运动前期和后期的附加质量力、Basset力、黏性阻力、浮力和水流推力进行对比, 选取与图3相同的参数进行数值模拟, 结果如图4所示.
图 4 气泡受力分析　(a)气泡下沉前期和后期各个分力的对比; (b) 气泡分力对气泡位移的贡献
Figure4. Force analysis of bubble: (a) Comparison of component forces in the early and late stages of bubble sinking; (b) contribution of bubble component force to bubble displacement.

图4(a)中${F_{{\rm{joint}}}}$表示气泡受到的合力, 图4(b)中${x_{{\rm{joint}}}}$, ${x_{\rm{B}}}$, ${x_{\rm{k}}}$, ${x_{\rm{D}}}$, ${x_{\rm{v}}}$, ${x_{\rm{p}}}$分别表示气泡受到的合力、Basset力、附加质量力、黏性阻力、浮力、水流推力对气泡位移的贡献, 其中合力对气泡位移的贡献就是气泡下沉的位移. 图4(a)可以观察到气泡下沉前期和后期各个分力的变化形式, 图4(b)可以观察到各个分力对位移的贡献, 其中附加质量力对气泡位移贡献的方向向下且贡献量最大, 因此附加质量力是气泡下沉的关键. 而Basset力对气泡位移贡献的大小约为附加质量力的1/20, 方向向下, 使得气泡轨迹相应的向下偏移一段距离, 在数值仿真误差允许范围里, 可以选择忽略此力. 此外, 可以观察到黏性阻力、浮力、水流推力对气泡的下沉起到阻碍作用, 且它们对气泡位移贡献的大小从大到小依次是: 浮力、黏性阻力和水流推力.

2

-->

4.1.实验装置

考虑到液体不同会对实验产生影响, 采用工业蒸馏水作为实验液体; 为了避免容器的振动迟滞对实验产生影响, 采用玻璃瓶作为实验容器, 其直径为30 mm, 高200 mm; 控制瓶子中空气柱的长度为30 mm.
采用型号为GF-100B型功率放大器, 型号为332DA Agilent型的信号示波器以及型号为JZ-5 6$ \perp $2的振动台, 选用恒压模式, 并对容器施加正弦波激励. 实验的频率范围为50—100 Hz, 振幅范围为0.5—5.0 mm; 振幅的测量通过型号为ISV2010的测振仪实现, 对于振动台和瓶子的连接, 自制了连接器, 并通过磁铁吸引, 保证瓶子和振动台的稳固连接, 通过水平仪对瓶子的垂直固定进行校准. 通过型号为DSC-RX100 M4的高速摄像机进行拍摄, 实验像素为1080 × 1920, 帧数范围为60—1000帧, 运用Tracker软件进行视频分析. 实验仪器如图5.
图 5 实验仪器　(a)振动台和瓶子; (b)功率放大器; (c)信号示波器
Figure5. Experimental instrument: (a) Shaker and bottle; (b) power amplifier; (c) signal oscilloscope.

2 -->

4.2.基于图像处理的气泡体积测量

为了更好地提取气泡特征尺寸, 首先要得到清晰的气泡图像, 采用逆光拍摄的方法^[28], 将柔光纸作为背景, 用强光逆向照射柔光纸, 从正面进行拍摄, 并通过数字图像处理来提取气泡特征尺寸.
图像处理的过程如下:
1)预处理. 将背景图像图6(a)和带有气泡的图像图6(b)转化为灰度图像, 对两个灰度图像采用剪影法处理^[26], 即对像素矩阵进行减法运算得到图6(c), 并通过$ {\rm{matlab}} $交互式地截取带有气泡的局部图像得到图6(d).
图 6 气泡图像处理　(a)背景图像; (b)带有气泡的图像; (c)差影法处理; (d)截取带有气泡的局部图像; (e)去噪; (f)$3 \times 3$滤波; (g)局部拉伸; (h)二值化处理; (i)填充图像
Figure6. Bubble image processing: (a) Background image; (b) image with bubbles; (c) subtraction processing; (d) local image with bubbles; (e) denoising; (f) $3 \times 3$ filtering; (g) local stretching; (h) image binarization; (i) image filling.

2)去噪、滤波. 先通过小波分解的方法对处理后的图像做平滑处理, 去除图像噪声得到图6(e), 再使用3 × 3的滤波器对图像进行中值滤波得到图6(f), 处理后气泡边缘一圈的像素灰度值变为0, 但不影响后续处理.
3)增强图像. 对处理后的图像进行线性变换, 以实现将区间内的像素拉伸到[0, 1]区间, 初始阈值选为[0.05, 0.3], 调整步长为0.005, 得到图6(g). 对增强后的图像选取阈值为0.3, 将像素灰度值在阈值以下的变为0, 阈值以上的变为1, 调整步长为0.03, 得到图6(h). 并对图像进行填充得到图6(i).
4)通过对像素灰度值为0的像素点计数得到气泡的像素面积S, 以等效半径来近似计算气泡的半径^[29,30], 其公式为

$R = \sqrt {S/{\rm{\pi }}}. $

2 -->

4.3.实验测量的不确定性

以等效半径计算气泡的半径, 考虑到气泡形状、空间分辨率等因素会影响气泡尺寸的提取, 需要对气泡尺寸测量的不确定度进行说明. 由于气泡的横纵比$\chi $在0.85—1之间, 对应的椭球离心率在0—0.5之间, 采用(37)式计算气泡半径的不确定度为 ± 2.6%^[31]. 受拍摄像素的限制, 气泡边界的像素点提取存在误差, 最大误差为 ± 1像素, 气泡半径的不确定度为 ± 6.5%, 在长度的标定过程中, 刚尺与相机无法完全平行, 标尺读数和实际长度存在误差, 实验中使用的刚尺为100 mm, 钢尺偏移不超过1 mm, 气泡半径的不确定度为 ± 0.3%.
2 -->

4.4.理论模型的检验

当水柱受到外界激励垂直振动时, 水面振动带动空气入水产生气泡, 气泡在水中以振荡的形式运动; 当气泡运动到临界深度X₀时, 其将在该深度附近稳定地上下振荡, 这种稳定的现象在不受到外界干扰时能持续很长时间. 保持水柱的振动频率不变, 当调高水柱振动的振幅, 稳定上下振荡的气泡将向更深处运动; 若降低振幅, 稳定上下振荡的气泡将向上运动.
从附加质量累积效应的角度分析气泡下沉的原因, 振动开始时, 液面的上下晃动给气泡带来一个向下的初速度, 气泡速度第一次变为0时, 其已经运动到了一定深度, 若气泡位于稳定振荡位置X₀以上, 气泡所受的压强变化较小, 气泡的体积变化也较小, 在一个周期内, 附加质量引起的累积效应不足以引起气泡下沉, 气泡开始上升; 若气泡位于稳定振荡位置X₀以下, 其所受的压强变化较大, 气泡的体积变化和附加质量的累积效应增大, 此时气泡开始下沉.
为减少水面扰动和气泡间的相互影响, 选取瓶子中下部的气泡作为观测对象, 测量气泡稳定振荡时的位置, 与(36)式给出的临界深度X₀进行对比, 得到图7. 结果表明, 在受到其他气泡影响较小的情况下, 气泡存在一个稳定振荡的位置, 实验中气泡稳定振荡的位置与(36)式给出的临界深度X₀相符, 相对误差在$ 5\% $以内.
图 7 理论值与实际值的对比
Figure7. Comparison between theoretical value and actual value.

通过对装有水的容器施加垂直振动, 观察气泡下沉的现象, 针对球形气泡构建了基于理想气体方程的气泡运动模型, 与前人的工作相比, 该模型引入了Basset力. 基于该模型, 本文研究了Basset力的理论求解, 与数值求解进行对比, 模拟出气泡所受Basset力的形式, 并详细分析了Basset力存在下气泡下沉的动力学特征. 模型的数值模拟结果表明垂直振动容器的水中气泡存在三种运动状态: 位于临界深度以下时气泡下沉, 位于临界深度以上时气泡上升, 位于临界深度附近时气泡稳定振荡. 阐明了气泡下沉的物理本质: 附加质量力降低气泡下沉的临界深度, 有利于气泡下沉; Basset力对气泡下沉的临界深度和气泡下沉的前期轨迹无明显影响, 但对气泡下沉的后期轨迹有较大影响. 对于气泡体积的测量, 采用数字图像处理的办法, 提高了气泡特征尺寸提取的精确性, 为工程应用提供指导意义.