2.1. Bubble oscillation equation

The physical model of this study is characterized by a spherical bubble with radius $R$ oscillating in the compressible liquid. The bubble oscillation is coupled with phenomena such as bubble migration and phase transition. The fluid domain is treated as weakly compressible and satisfies the linear wave equation (Zhang et al. Reference Zhang, Li, Cui, Li and Liu2023). With the centre of the bubble as the coordinate origin $o$, the wave equation in the spherical coordinate system $o\unicode{x2013} r \theta \phi$ is expressed as

(2.1)\begin{equation} \frac{1}{{{C}^{2}}}\,\frac{{{\partial }^{2}}\varphi }{\partial {{t}^{2}}}=\frac{1}{{{r}^{2}}}\,\frac{\partial }{\partial r}\left( {{r}^{2}}\,\frac{\partial \varphi }{\partial r} \right)+\frac{1}{{{r}^{2}}\sin \theta }\,\frac{\partial }{\partial \theta }\left( \sin \theta\,\frac{\partial \varphi }{\partial \theta } \right)+\frac{1}{{{r}^{2}}{{\sin }^{2}}\theta }\,\frac{{{\partial }^{2}}\varphi }{\partial {{\phi }^{2}}}, \end{equation}

where $\varphi$ is the velocity potential of liquids, and $C$ is the sound speed.

Assuming that the bubble keeps spherical oscillation, we define here that the bubble migrates along the direction $\theta =0$. Thus $\partial \varphi /\partial \phi =0$, and the third term on the right-hand side of (2.1) vanishes because of the axisymmetry. The velocity potential $\varphi _{f_{s}}$ induced by a source at the origin with strength $f_{s}(t)$ is

(2.2)\begin{equation} \varphi_{f_{s}}({\boldsymbol r},t) ={-}\frac{1}{{|\boldsymbol r|}}\, f_{s}\left(t-\frac{|\boldsymbol r|}{C}\right). \end{equation}

According to the linearity characteristics of the wave equation, the superposition of a set of $\phi _{{f}_{s}}$ can produce the solution of the wave equation $\varphi _{f}$ considering the source movement:

(2.3)\begin{equation} {\varphi_{f}(\boldsymbol{r},t) ={-}\frac{C}{(C-\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{r}_{t}/ |\boldsymbol{r}_{t}|)\,|\boldsymbol{r}_{t}|}\, f\left(t-\frac{|\boldsymbol{r}_{t}|}{C}\right)}, \end{equation}

where $\boldsymbol {r}_{t}$ is the vector pointing from the source at $t-|\boldsymbol {r}_{t}|/C$ to $\boldsymbol {r}$, and $f$ is a function whose second-order derivative exists. The relative velocity vector $\boldsymbol {v}$ represents the velocity difference between the bubble migration velocity ${\boldsymbol v}_{m}$ and the ambient flow velocity ${\boldsymbol u}_{a}$. Once (2.3) is obtained, the velocity potential of the moving dipole $\varphi _{q}$ can be calculated as

(2.4)\begin{equation} \varphi_{q} = \lim_{D \rightarrow0}\frac{1}{2D}\, (\varphi_{f}(\boldsymbol{r}+\boldsymbol{e}D,t)- \varphi_{f}(\boldsymbol{r}-\boldsymbol{e}D,t)), \end{equation}

where the unit vector $\boldsymbol {e}$ indicates the direction along which the bubble migrates, and $D$ denotes the distance that is halfway between the point source and the sink.

Considering that the migration velocity is small relative to sound speed, the location variation of singularities could be ignored during the short time when the influences propagate from the singularities to the bubble surface. Thus (2.3) could be simplified, and the linear superposition of the velocity potential of the point source and dipole can be expressed as

(2.5)\begin{equation} \varphi{(r,\theta,t)} ={-}\frac{1}{r}\,f{\left(t-\frac{r}{C}\right)} - \frac{\mathcal{\cos\theta}}{r^2}\,q{\left(t-\frac{r}{C}\right)} -\frac{1}{C}\,\frac{\cos\theta}{r}\,q'{\left(t-\frac{r}{C}\right)}, \end{equation}

where $q$ is a function whose second-order derivative exists, and $q'$ represents the derivative of $q$ with respect to its argument.

The time derivative of $\varphi$ and the velocity of the bubble surface in the $r$ direction are

(2.6)\begin{equation} \frac{\partial \varphi}{\partial t} ={-}\frac{1}{r}\,f'{\left(t-\frac{r}{C}\right)} - \frac{\mathcal{\cos\theta}}{r^2}\,q'{\left(t-\frac{r}{C}\right)} - \frac{1}{C}\,\frac{\cos\theta}{r}\,q''{\left(t-\frac{r}{C}\right)} \end{equation}

and

(2.7)\begin{equation} u_{r} = {{\left. \frac{\partial \varphi }{\partial r} \right|}_{r=R}} = \frac{f}{R^2}+ \frac{f'}{CR} +\frac{2\cos\theta}{R^3}\,q + \frac{2\cos\theta }{CR^2 }\,q', \end{equation}

respectively, where $f'$ represents the derivative of $f$ with respect to its argument, and $q''$ denotes the second derivative of $q$ with respect to its argument. The terms of magnitude $1/C^2$ are ignored in (2.7) and subsequent derivations.

According to the continuity condition at the bubble surface considering the phase transition (Fujikawa & Akamatsu Reference Fujikawa and Akamatsu1980), the normal velocity of fluids at the bubble surface is $\dot R - \dot m/\rho$, where $\dot {R}$ denotes the time derivative of $R$, $\dot m$ is the net evaporation rate of mass per unit area of the bubble surface, and $\rho$ is the liquid density. Then the averaged kinetic boundary condition for the bubble oscillation can be expressed by

(2.8)\begin{equation} \int_{S_{b}} u_{r}\,\mathrm{d}S = 4{\rm \pi} R^2 \left(\dot R- \frac{\dot m}{\rho}\right), \end{equation}

where $S_{b}$ denotes the area of the bubble surface.

Substituting (2.7) into (2.8), $q$ and $q'$ vanish, so

(2.9)\begin{equation} \frac{f}{R^2}+ \frac{f'}{CR} = \dot R- \frac{\dot m}{\rho}. \end{equation}

Let $F(t) = f(t - R/C)$. Then we have

(2.10)\begin{equation} \frac{{\rm d}F}{{\rm d}t}=\left( {1 - \frac{{\dot R}}{C}} \right)f'{|_{r = R}}. \end{equation}

Taking the derivative of (2.9) with respect to time and combining it with (2.10) yields

(2.11)\begin{equation} \left( {\frac{{C - \dot R}}{R} + \frac{{\rm d}}{{{{\rm d}t}}}} \right)\left( {\frac{R}{{C - \dot R}}\,\frac{{{\rm d}F}}{{{\rm d}t}}} \right) +\frac{\text{d}}{\text{d}t}\left( \frac{{{R}^{2}}\dot{m}}{\rho } \right) = 2R{\dot R^2} + {R^2}\ddot R, \end{equation}

where ${{{\rm d}F}}/{{{\rm d}t}}$ depends on different physical problems and environmental conditions. Once ${{{\rm d}F}}/{{{\rm d}t}}$ and $\dot m$ are determined, (2.11) can be solved to obtain the bubble oscillation dynamics. As the effect of phase transition is neglected, the above equation is simplified to

(2.12)\begin{equation} \left( {\frac{{C - \dot R}}{R} + \frac{{\rm d}}{{{{\rm d}t}}}} \right)\left( {\frac{R}{{C - \dot R}}\,\frac{{{\rm d}F}}{{{\rm d}t}}} \right) = 2R{\dot R^2} + {R^2}\ddot R . \end{equation}

To obtain ${{{\rm d}F}}/{{{\rm d}t}}$, we apply the Bernoulli equation to the bubble surface in the moving coordinate system:

(2.13)\begin{equation} {{\left. \frac{\partial \varphi }{\partial t} \right|}_{r=R}} - \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{u}+ \frac{1}{2}\,{| \boldsymbol{u} |^2} + H = 0, \end{equation}

where $\boldsymbol {v} \boldsymbol {\cdot } \boldsymbol {u}$ can be expressed as the inner product of the vectors $( v\cos \theta, -v\sin \theta )$ and $( u_{r}, u_\theta )$, with $v=|\boldsymbol v|$. Here, $H=\int _{{P_{{a}}}}^{{P_{{b}}}} {\rho ^{-1}} \,\mathrm {d} p$ is the enthalpy difference at the bubble surface, and its leading order term can be expressed as $(P_{{b}}-P_{{a}})/\rho$, where $P_{a}$ represents the ambient pressure at the bubble centre ($P_{{a}}$ includes the hydrostatic pressure at infinity $P_{\infty }$, the acoustic pressure, and the pressures induced by boundaries and other bubbles), and $P_{b}$ denotes the liquid pressure at the bubble surface.

According to (2.5), (2.7) and (2.9), the normal and tangential velocities of the bubble surface can be expressed as

(2.14)\begin{equation} {{u}_{r}}=\dot{R}-\frac{\dot{m}}{\rho }+\frac{2\cos \theta }{{{R}^{3}}}\,q+\frac{2\cos \theta }{C{{R}^{2}}}\,q' \end{equation}

and

(2.15)\begin{equation} u_{\theta} = \frac{1}{R} {{\left. \frac{\partial \varphi }{\partial \theta} \right|}_{r=R}} = \frac{\mathcal{\sin\theta}}{R^3}\,q +\frac{\sin\theta}{C R^2}\,q', \end{equation}

respectively.

Integrating the Bernoulli equation over the bubble surface, all the terms containing $\theta$ could be eliminated. Consequently, both $q$ and $q'$ are eliminated, and an equation containing only the unknown quantity $f'$ can be obtained such that

(2.16) \begin{align} &\int_{S} {\left( {{{\left. {\frac{{\partial \varphi }}{{\partial t}}} \right|}_{r = R}} - {\boldsymbol v} \boldsymbol{\cdot} {\boldsymbol u} + \frac{1}{2}\,{{| {\boldsymbol u} |}^2} + H} \right) {\rm d}S} \nonumber\\ &\quad = \int_0^{\rm \pi} {\left.{\left[\!\begin{aligned} & -\!\frac{1}{r}\,f' - \frac{{\cos \theta }}{{{r^2}}}\,q' - \frac{1}{C}\,\frac{{\cos \theta }}{r}\,q'' - ( {v\cos \theta , - v\sin \theta } ) \\ & \quad \boldsymbol{\cdot} ( {{u_{r}},{u_{\theta} }} )+\frac{1}{2}\left( {u_{r}^2 + u_\theta ^2} \right) + H \end{aligned}\!\right]}\right|_{r=R} 2{\rm \pi} {R^2}\sin \theta \, {\rm d}\theta } \nonumber\\ &\quad = 4{\rm \pi} {R^2}\left( { - \frac{{f'}}{R} + \frac{1}{2}{{\left( {\dot R - \frac{{\dot m}}{\rho }} \right)}^2} + \frac{1}{4}\,{v^2} + H} \right) = 0. \end{align}

Combining (2.10) and (2.16), the expression for ${{{\rm d}F}}/{{{\rm d}t}}$ in (2.11) can be obtained:

(2.17)\begin{equation} \frac{{\rm d}F}{{\rm d}t}= R \left( {1 - \frac{{\dot R}}{C}} \right) \left(\frac{1}{2}{\left({{\dot R}-\frac{\dot m}{\rho}} \right)^2}+\frac{1}{4}\,{{v}^{2}}+H \right). \end{equation}

Substituting the above expression into (2.11), the bubble oscillation equation considering bubble migration and phase transition can be provided as

(2.18)\begin{align} \left( {\frac{{C - \dot R}}{R} + \frac{{\rm d}}{{{\rm d}t}}} \right) \left[{\frac{{{R^2}}}{C}\left( {\frac{1}{2}{\left({{\dot R}-\frac{\dot m}{\rho}} \right)^2} + \frac{1}{4}\,{{v}^2} + H} \right)} \right] +\frac{\text{d}}{\text{d}t}\left( \frac{{{R}^{2}}\dot{m}}{\rho } \right) = 2R\dot R^2 + {R^2}\ddot R . \end{align}

Equation (2.18) is the bubble oscillation equation with a unified mathematical form, which takes into account the multiple physical factors. The first term on the left-hand side represents the coupling force of the bubble oscillation, migration and ambient flow field; the second term is the source term due to the phase transition. The right-hand side represents the volume acceleration of the bubble. The enthalpy difference $H$ exhibits good extensibility, determined by the specific physical problems. Once $H$, $\dot m$ and $v$ are obtained, (2.18) can be solved. As $\dot m = 0$, the above equation is simplified to

(2.19)\begin{equation} \left( {\frac{{C - \dot R}}{R} + \frac{{\rm d}}{{{\rm d}t}}} \right) \left[{\frac{{{R^2}}}{C}\left( {\frac{1}{2}\,{{{\dot R}}^2} + \frac{1}{4}\,{{v}^2} + H} \right)} \right] = 2R\dot R^2 + {R^2}\ddot R . \end{equation}

The above equation simplifies to the Keller–Miksis equation when the bubble migration velocity is removed, and transforms to the Gilmore equation by making simple substitutions in the expansion. It can be simplified to the Rayleigh–Plesset equation if fluid compressibility is further neglected.

2.2. Phase transition modelling

The key to solving the enthalpy difference $H$ at the bubble surface is to obtain the liquid pressure at the bubble surface, which is closely related to the phase transition. According to Fujikawa & Akamatsu (Reference Fujikawa and Akamatsu1980), the pressure balance on the surface of the bubble can be expressed as

(2.20)\begin{equation} {{P}_{b}}={{P}_{g}}-\frac{2\sigma }{R}-\frac{4\mu }{R}\left( \dot{R}-\frac{{\dot{m}}}{\rho } \right)-{{\dot{m}}^{2}}\left( \frac{1}{\rho }-\frac{1}{{{\rho }_{g}}} \right), \end{equation}

in which ${{P}_{g}}$ is the inner gas pressure, $\sigma$ is the surface tension coefficient, $\mu$ is the viscosity coefficient, and ${\rho }_{g}$ is the average gas density.

The bubble contents consist mainly of non-condensable gases and vapour (Brenner, Hilgenfeldt & Lohse Reference Brenner, Hilgenfeldt and Lohse2002), which are considered to be uniformly distributed inside the bubble. Considering that the gases inside the bubble are violently compressed at the end of the collapse, the van der Waals equation (Yasui Reference Yasui1997; Kyuichi Reference Kyuichi2021) is employed to model the uniform inner pressure of the bubble:

(2.21)\begin{equation} \left( {{P}_{g}}+\frac{a}{{{\upsilon_{m} }^{2}}} \right)( \upsilon_{m} -b )={{R}_{g}}T, \end{equation}

where $\upsilon _{m}= N_{A}V/n_{t}$ is the molar volume ($N_{A}$ is the Avogadro number, $V$ is the volume of the bubble, $n_{t}$ denotes the number of molecules inside the bubble), $T$ is the temperature at the bubble centre, $R_{g}$ is the gas constant, and $a$ and $b$ are van der Waals constants with the expressions

(2.22) \begin{equation} \left\{ {\begin{array}{@{}l@{}} {a = {{\left({\sqrt {{a_{a}}}\,\dfrac{{{n_{a}}}}{{{n_{t}}}} + \sqrt {{a_{v}}}\,\dfrac{{{n_{v}}}}{{{n_{t}}}}} \right)}^2}},\\ {b = {{\left( {\sqrt {{b_{a}}}\, \dfrac{{{n_{a}}}}{{{n_{t}}}} + \sqrt {{b_{v}}}\, \dfrac{{{n_{v}}}}{{{n_{t}}}}} \right)}^2}}. \end{array}} \right. \end{equation}

Here, ${n}_{a}$ and ${n}_{v}$ are the numbers of air and vapour molecules, respectively, with ${n_{t}}={n}_{a}+{n}_{v}$, ${a}_{a}$ and ${a}_{v}$ are the van der Waals forces of air and vapour molecules, respectively, and ${b}_{a}$ and ${b}_{v}$ are the volumes occupied by air molecules and vapour molecules, respectively. As the van der Waals constants ($a$ and $b$) equal zero, (2.21) simplifies to the ideal gas equation. The change rate in the number of vapour molecules can be calculated as

(2.23)\begin{equation} {{{\dot{n}}}_{v}}=\frac{4{\rm \pi} {{R}^{2}}\dot{m}{{N}_{A}}}{{{M}_{{mv}}}}, \end{equation}

where $M_{{mv}}$ is the molar mass of vapour.

A modified Hertz–Knudsen–Langmuir relationship (Schrage Reference Schrage1953; Akhatov et al. Reference Akhatov, Lindau, Topolnikov, Mettin, Vakhitova and Lauterborn2001) could be used to compute the net evaporation rate of mass:

(2.24)\begin{equation} \dot{m}=\frac{{{\alpha }_{m}}}{\sqrt{2{\rm \pi} {{R}_{v}}}}\left( \frac{P_{s}}{\sqrt{{{T}_{l}}(R)}}-\frac{\varGamma {{P}_{v}}}{\sqrt{{{T}_{B}}}} \right), \end{equation}

where ${\alpha }_{m}$ is the adaptation factor, and its value is characterized by the evaporation and condensation. (The value of ${\alpha }_{m}$ is taken to be approximately 0.04 according to previous works on acoustic bubbles (Yasui Reference Yasui1998), but it is not a fixed value for bubbles generated by different methods; it is determined depending on the specific properties of the bubble contents.) Also, $R_{v}$ is the gas constant for vapour, $T_{l}(R)$ denotes the liquid temperature at the bubble surface, $T_{B}$ is the gas temperature at the bubble surface, $P_{s}$ is the saturated vapour pressure at the temperature $T_{l}(R)$, and $P_{v}$ is the actual saturated vapour pressure. Following Han et al. (Reference Han, Chen and Guo2023), $\varGamma$ is a correction factor taken as 1.0. Additionally, for bubbles where phase transition effects are minimal but there is a flow of gases into and out of the bubble, such as in the case of air-gun bubbles (de Graaf, Penesis & Brandner Reference de Graaf, Penesis and Brandner2014; Li et al. Reference Li, der Meer, Zhang, Prosperetti and Lohse2020), the value of $\dot m$ in (2.24) can be computed in alternative forms according to different physical problems.

To describe the temperature change at the gas–liquid interface, two thermodynamic boundary layers are introduced according to previous works (Fujikawa & Akamatsu Reference Fujikawa and Akamatsu1980; Yasui Reference Yasui1999). Inside the bubble, the gas temperature varies from the temperature $T$ inside the bubble to $T_{B}$ at the bubble surface. Outside the bubble, the liquid temperature $T_{l}(r)$ changes from the temperature $T_{l}(R)$ at the bubble surface to $T_{\infty }$ at infinity. Here, following the linear model proposed by Yasui, Tuziuti & Kanematsu (Reference Yasui, Tuziuti and Kanematsu2016), the gas temperature distribution near the inner bubble surface could be described as

(2.25)\begin{equation} {{\left. \frac{\partial T}{\partial r} \right|}_{r=R}}=\frac{{{T}_{B}}-T}{n\lambda }, \end{equation}

where $\lambda$ is the mean molecular free range of the gas, and $n$ is a constant that determines the thickness of the thermodynamic boundary layer.

Assuming that there is a temperature jump at the gas–liquid interface from $T_{B}$ to $T_{l}(R)$ (Kogan Reference Kogan1969; Yasui Reference Yasui1995), the gas temperature $T_{B}$ at the bubble surface is calculated by

(2.26)\begin{equation} {{T}_{B}}={{T}_{l}(R)}-\frac{\kappa }{2kn'}\sqrt{\frac{{\rm \pi} \bar{m}}{2k{{T}_{B}}}}\,\frac{2-a'{{\alpha }_{e}}}{{{\alpha }_{e}}}{{\left. \frac{\partial T}{\partial r} \right|}_{r=R}}, \end{equation}

where $k$ denotes the Boltzmann constant, $n'$ denotes the number density of molecules inside the bubble, ${\alpha }_{e}$ is the thermodynamic coefficient of adaptation, $a'=0.827$, $\kappa$ is the thermal conductivity coefficient of water vapour, and $\bar {m}$ is the average mass of the molecules inside the bubble, $\bar {m}={( {{n}_{v}}{{M}_{v}}+{{n}_{a}}{{M}_{a}} )}/{( {{n}_{t}}{{N}_{A}} )}$ (with ${M}_{a}$ and ${M}_{v}$ the masses of air and vapour inside the bubble, respectively).

The temperature distributions within the thermodynamic boundary layer outside the bubble surface are examined extensively (Fujikawa & Akamatsu Reference Fujikawa and Akamatsu1980; Tian et al. Reference Tian, Zhang, Yin, Lv and Zhu2022; Dai et al. Reference Dai, Zhu, Wang, Chu and Wang2024), as they need to satisfy the thermodynamic boundary conditions both at the bubble surface and at infinity. Here, the exponential distribution model proposed by Yasui (Reference Yasui1996) is employed to describe the temperature gradient outside the bubble at the collapse stage:

(2.27) \begin{equation} {T_{l}}( r ) =\left\{ \begin{array}{@{}l@{\quad}l@{}} {T_\infty } + [{{T_{l}}( R ) - {T_\infty }} ]\exp\left(-\dfrac{{r - R}}{{{T_\infty } - {T_{l}}( R )}}{{\left.{\dfrac{{\partial {T_{l}}}}{{\partial r}}}\right|}_{r = R}}\right) & \text{when }[ {{T_{l}}( R ) - {T_\infty }} ]{{\left. {\dfrac{{\partial {T_{l}}}}{{\partial r}}} \right|}_{r = R}} <0,\\ {T_\infty } + A\exp(- B( {r - Y} )^2) & \text{when }[ {{T_{l}}( R ) - {T_\infty }} ]{{\left. {\dfrac{{\partial {T_{l}}}}{{\partial r}}} \right|}_{r = R}} > 0, \end{array} \right. \end{equation}

where $A$, $B$ and $Y$ are computed as

(2.28a–d) \begin{align} A = [ {{T_{l}}( R ) - {T_\infty }} ]\, {{\rm e}^{Be_1^2}}, B = \frac{{{{\left. {\dfrac{{\partial {T_{l}}}}{{\partial r}}} \right|}_{r = R}}}}{{2{e_1}[ {{T_{l}}( R ) - {T_\infty }} ]}}, Y = R + {e_1}, {e_1} = e'\left| {\frac{{{T_{l}}( R ) - {T_\infty }}}{{{{\left. {\dfrac{{\partial {T_{l}}}}{{\partial r}}} \right|}_{r = R}}}}} \right|, \end{align}

with $e'=0.001$.

Then, according to the energy conservation within the thermodynamic boundary layer outside the bubble, the variation of $T_l(R)$ with time could be updated as

(2.29) \begin{align} \frac{4}{3}\,{\rm \pi} \rho {c_{p}}[ {{{( {R \,{+}\, {\delta _{e}}} )}^3} \,{-}\, {R^3}} ]\, \frac{\partial {T_{l}}(R)}{\partial t} \,{=}\, 4{\rm \pi} {R^2}\left(\! { - {\kappa_{l}}{{\left. {\frac{{\partial {T_{l}}}}{{\partial r}}} \!\right|}_{r = R}}} \!\right) \,{-}\, 4{\rm \pi} {( {R + {\delta _{e}}} )^2}\left(\! { - {\kappa_{l}}{{\left. {\frac{{\partial {T_{l}}}}{{\partial r}}} \!\right|}_{r = R + {\delta _{e}}}}}\!\right) , \end{align}

where $c_{p}$ denotes the specific heat of liquids at constant pressure, and $\kappa _{l}$ is the thermal conductivity of liquids. The thickness $\delta _{e}$ of the thermodynamic boundary layer outside the bubble (Yasui Reference Yasui1996, Reference Yasui1997) is estimated by

(2.30) \begin{equation} {\delta _{e}} =\left\{ \begin{array}{ll} {{[ {{T_{l}}( R ) - {T_\infty }} ]}\left/{{{\left. {\dfrac{{\partial {T_{l}}}}{{\partial r}}} \right|}_{r = R}}}\right.} & \text{when }[ {{T_{l}}( R ) - {T_\infty }} ]{{\left. {\dfrac{{\partial {T_{l}}}}{{\partial r}}} \right|}_{r = R}} <0,\\ {e_1} + 1/\sqrt B & \text{when }[ {{T_{l}}( R ) - {T_\infty }} ]{{\left. {\dfrac{{\partial {T_{l}}}}{{\partial r}}} \right|}_{r = R}} > 0, \end{array} \right. \end{equation}

where ${{{({\partial {T_{l}}}/{\partial r})}} |_{r = R + {\delta _{e}}}}$ is calculated using (2.27) and (2.28a–d); the temperature gradient of liquid at the bubble surface ${{{({\partial {T_{l}}}/{\partial r})}} |_{r = R}}$ is determined by the continuity condition (Fujikawa & Akamatsu Reference Fujikawa and Akamatsu1980) of the heat flux:

(2.31)\begin{equation} {\kappa_{l}}{\left. {\frac{{\partial {T_{l}}}}{{\partial r}}} \right|_{r = R}} = {\kappa}{\left. {\frac{{\partial {T}}}{{\partial r}}} \right|_{r = R}} + {\dot m}L, \end{equation}

where $L$ is latent heat.

The temperature change inside the bubble with respective to time can be updated according to the change of internal energy

(2.32)\begin{equation} M {c}_{v} {\dot{T}} = \dot E, \end{equation}

where $M$ is the gas mass inside the bubble, and ${c}_{v}$ is the average specific heat capacity of the gas inside the bubble. Following Zhong et al. (Reference Zhong, Eshraghi, Vlachos, Dabiri and Ardekani2020) and Yasui (Reference Yasui2001), the change rate of energy is computed as

(2.33)\begin{align} \dot{E}&={-}S_{b}\dot{R}{{P}_{g}}+\frac{S_{b}[ {{{\dot{m}}}_{{e}}}\,{{e}_{v}}( {{T}_{l}} )-{{{\dot{m}}}_{{c}}}\,{{e}_{v}}( {{T}_{B}} ) ]{{N}_{A}}}{{{M}_{{mv}}}} +S_{b}\kappa {{\left. \frac{\partial T}{\partial r} \right|}_{r=R}} \nonumber\\ & \quad +S_{b}{{\sigma }_{r}}( T_{B}^{4}-{{T}^{4}} ), \end{align}

in which ${e}_{v}$ is the energy carried by a vapour molecule, and ${\sigma }_{r}$ is the Stefan–Boltzmann constant. The first term on the right-hand side is the work done by the bubble on the surrounding fluids; the second term represents the energy carried by the evaporation and condensation of the fluids; the third term is the energy produced by heat conduction; and the fourth term is the energy produced by heat radiation. In some studies (de Graaf et al. Reference de Graaf, Penesis and Brandner2014; Nagalingam et al. Reference Nagalingam, Raghunathan, Korede, Poelma, Smith, Hartkamp, Padding and Eral2023; Chen et al. Reference Chen, Chen, Geng, Huang and Cao2024), the two thermodynamic boundary layers at the gas–liquid interface are often ignored for simplicity. The gas temperature at the inner surface of the bubble is replaced with the temperature at the bubble centre, while the liquid temperature at the outer surface of the bubble is substituted with the ambient temperature. Then the change rate of energy can be expressed simply as

(2.34)\begin{equation} \dot{E}={-}S_{b}\dot{R}{{P}_{g}}+{S_{b}} ({{c}_{p}}-{{c}_{v}}){{T}}\dot{m}-S_{b}\kappa_{s} ( {{T}}-{{T}_{\infty}} ), \end{equation}

where $\kappa _{s}$ is a heat transfer coefficient.

Here, some values of the parameters at room temperature (Zhong et al. Reference Zhong, Eshraghi, Vlachos, Dabiri and Ardekani2020) involved in the above model are provided as shown in table 1.

Table 1. Values of partial thermodynamic parameters at room temperature (293 K).

2.3. Bubble migration equation

In this subsection, the bubble migration equation is derived to solve the migration velocity in (2.18). As the bubble migrates along the axis $\theta =0$, the kinetic boundary condition for the bubble migration can be expressed as

(2.35)\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}\int_V r\cos\theta \,\mathrm{d}V = \frac{4}{3}\,{\rm \pi} {R^3}{ v}. \end{equation}

According to Reynolds’ transport theorem, the above equation can be expanded as

(2.36)\begin{equation} \int_V {\frac{{\partial ( {r\cos \theta } )}}{{\partial t}}\,{\rm d}V} + \int_S {u_{r}R\cos \theta \,{\rm d}S} = \frac{4}{3}\,{\rm \pi} {R^3}{ v}. \end{equation}

The first term on the left-hand side in (2.36) disappears because the integrand does not change over time. Combining the above equation with (2.7) yields

(2.37)\begin{equation} q+\frac{R}{C}\,q' =\frac{1}{2}\,{ {v} }{{R}^{3}}. \end{equation}

Let $Q(t) = q(t - R/C)$. Then we have

(2.38)\begin{equation} \frac{{\rm d}Q}{{\rm d}t}=\left( {1 - \frac{{\dot R}}{C}} \right)q'{|_{r = R}}. \end{equation}

Combining the above two equations and differentiating (2.37) with respect to time gives

(2.39)\begin{equation} \left( \frac{C-\dot{R}}{R}+\frac{\rm d}{{\rm d}t} \right)\left( \frac{R}{C-\dot{R}}\,\frac{{\rm d}Q}{{\rm d}t} \right) =\frac{1}{2}\,{\dot{ v}} {{R}^{3}}+\frac{3}{2}\,{{R}^{2}}\dot{R}{ {v} }. \end{equation}

Equation (2.39) is the bubble migration equation in a unified mathematical form. The left-hand side of (2.39) represents the migration force exerted on the bubble by the flow field, while the right-hand side is the change rate of the bubble's momentum with respect to time. Similar to the bubble oscillation equation (2.11), once we have determined ${{\rm d}Q}/{{\rm d}t}$, we can obtain the migration equation for the bubble. To determine ${{\rm d}Q}/{{\rm d}t}$, we first set out the momentum equation for the bubble:

(2.40)\begin{equation} \frac{{\rm d}( M{\boldsymbol v}_{m} )}{{\rm d}t}= M{\boldsymbol g} -\int_{S}{{{P}_{b}}{\boldsymbol n}\,\text{d}S}-\frac{1}{2}\,{\rm \pi} {{R}^{2}}\rho{{C}_{d}}\,\mathbb{C}( \boldsymbol{v}), \end{equation}

where the three terms on the right-hand side denote the gravity, the inertial force, and the drag force of the bubble, respectively, with $C_{d}$ the drag coefficient, and $\mathbb {C}({\boldsymbol x})={\boldsymbol x}\,|\boldsymbol x|$.

Multiplying the Bernoulli equation (2.13) by $\boldsymbol n$ and integrating it on the bubble surface with $H$ retaining the zero-order term $(P_{b}-P_{a})/\rho$, the terms containing $\theta$ vanish. Consequently, the inertial force can be obtained as

(2.41) \begin{align} &\int_S {{P_{b}}{\boldsymbol n}\,{\rm d}S} = \int_S {\left[ {{P_{{a}}} - \rho {{\left. {\left( {\frac{{\partial \varphi }}{{\partial t}} + ( {v\cos \theta , - v\sin \theta } ) \boldsymbol{\cdot} ( {{u_{r}},{u_{{\theta}} }} ) + \frac{1}{2}\,( {u_{r}^2 + u_{\theta} ^2} )} \right)} \right|}_{r = R}}} \right]{\boldsymbol n}\,{\rm d}S} \nonumber\\ &\quad = \int_V {\boldsymbol{\nabla} {P_{{a}}}\,{\rm d}V} - \rho \int_0^{\rm \pi} {\left( {\left. {\frac{{\partial \varphi }}{{\partial t}}} \right|_{r = R}} + ( {v\cos \theta , - v\sin \theta } ) \boldsymbol{\cdot} ( {{u_{r}},{u_{\theta} }} ) + \frac{1}{2}\,( {u_{r}^2 + u_{\theta} ^2} ) \right) }\nonumber\\ &\qquad \times {\boldsymbol e} \boldsymbol{\cdot} 2{\rm \pi} {R^2}\cos \theta \sin \theta \, {\rm d}\theta = \frac{4}{3}\,{\rm \pi} \rho \left( {q' + \frac{R}{C}\,q''} \right){\boldsymbol e} + \frac{4}{3}\,{\rm \pi} R^3\,\boldsymbol{\nabla} {P_{{a}}}, \end{align}

where $\boldsymbol {\nabla } {P_{a}}=\rho {\boldsymbol g}$ in the free field with gravity.

Substituting (2.41) into (2.40) and organizing gives

(2.42)\begin{equation} \left(q' + \frac{R}{C}\,q''\right) {\boldsymbol e}={-} \frac{{{R^3}\,\boldsymbol{\nabla} {P_{{a}}}}}{\rho } - \frac{3}{8}\,{R^2}{C_{d}}\,\mathbb{C}( {\boldsymbol v} ) + \frac{{{\rho _{g}}}}{\rho }\,{R^3}{\boldsymbol g} - \frac{{{\rho _{g}}}}{\rho }\,{R^3}{\dot {\boldsymbol v}}_{m} - \frac{{3{R^2}\dot m}}{\rho }\,{\boldsymbol v}_{m}. \end{equation}

Multiply both sides of (2.37) by $\boldsymbol e$ and differentiate it with respect to time. Then, associating it with (2.42), we can eliminate $q''$ to have

(2.43)\begin{align} \frac{{\dot R}}{{C - \dot R}}\,q'{\boldsymbol e} &= \frac{C}{{C - \dot R}}\left( {\frac{1}{2}\,\dot {\boldsymbol v}{R^3} + \frac{3}{2}\,{\boldsymbol v}{R^2}\dot R} \right) \nonumber\\ &\quad - \left( { - \frac{{{R^3}\,\boldsymbol{\nabla} {P_{{a}}}}}{\rho } - \frac{3}{8}\,{R^2}{C_{d}}\,\mathbb{C}( \boldsymbol v ) + \frac{{{\rho _g}}}{\rho }\,{R^3}{\boldsymbol g} - \frac{{{\rho _g}}}{\rho }\,{R^3}{\dot {\boldsymbol v}}_{m} - \frac{{3{R^2}\dot m}}{\rho }\,{\boldsymbol v}_{m}} \right). \end{align}

By associating (2.38) with the above equation, the expression for ${{\rm d}Q}/{{\rm d}t}$ is obtained:

(2.44)\begin{align} \frac{\text{d}Q}{\text{d}t}\,{\boldsymbol e}&=\frac{{C - \dot R}}{{\dot R}}\left( {\frac{1}{2}\,\dot {\boldsymbol v}{R^3} + \frac{3}{2}\,{\boldsymbol v}{R^2}\dot R} \right) \nonumber\\ &\quad + \frac{{{{( {C - \dot R} )}^2}}}{{C\dot R}}\left( { \frac{{{R^3}\,\boldsymbol{\nabla} {P_{{a}}}}}{\rho } + \frac{3}{8}\,{R^2}{C_{d}}\,\mathbb{C}( \boldsymbol v ) - \frac{{{\rho _{g}}}}{\rho }\,{R^3}{\boldsymbol g} + \frac{{{\rho _{g}}}}{\rho }\,{R^3}{\dot {\boldsymbol v}}_{m} + \frac{{3{R^2}\dot m}}{\rho }\,{\boldsymbol v}_{m}} \right). \end{align}

Similarly, multiplying both sides of (2.39) by $\boldsymbol e$ and associating it with (2.44), we can arrive at the bubble migration equation:

(2.45) \begin{align} & \left[ 1-\frac{R\ddot{R}}{( C-\dot{R} )\dot{R}}+\frac{R}{C-\dot{R}}\,\frac{\text{d}}{\text{d}t} \right]\left( \frac{1}{2}\,{{R}^{3}}\dot{\boldsymbol v}+\frac{3}{2}\,{{R}^{2}}\dot{R}{\boldsymbol v} \right) \nonumber\\ &\quad =\left[\! 1-\frac{R\ddot{R}}{( C-\dot{R})\dot{R}}+\frac{R}{C}\,\frac{\text{d}}{\text{d}t} \!\right]\left[\!\frac{{{\rho }_{g}}}{\rho }\,R^3({\boldsymbol g} - {\dot{\boldsymbol v}}_{m} ) -3{{R}^{2}}\,\frac{{\dot{m}}}{\rho }\,{\boldsymbol v}_{m} -\frac{{{R}^{3}}\,\boldsymbol{\nabla} {{P}_{a}}}{\rho }-\frac{3}{8}\,{{C}_{{d}}}{{R}^{2}}\,\mathbb{C}( \boldsymbol v )\!\right]. \end{align}

In the free field, ${\boldsymbol u}_{a}=0$, thus ${\boldsymbol v}_{m}={\boldsymbol v}$. When we consider the non-spherical bubble oscillation in many cases, the added mass coefficient of the bubble $C_{a}$ needs to be introduced in (2.45). In fact, in the above equations, the added mass coefficient of the bubble is implicitly fixed to 0.5 due to the assumption of spherical bubbles. By analogy with the derivation of the added mass force of the bubble in an incompressible flow, the expression in parentheses on the left-hand side of (2.45) can be rewritten as $(C_{a}R^3\dot {{\boldsymbol v}} + 3C_{a}R^2\dot {R} {\boldsymbol v})$ if $C_{a}$ is not equal to 0.5.

When the high-order terms related to fluid compressibility are neglected, the above equation simplifies to

(2.46)\begin{equation} C_{a}{{R}}\dot{\boldsymbol v}+3 C_{a} \dot{R}{\boldsymbol v} +\frac{{{R}}\,\boldsymbol{\nabla} {{P}_{a}}}{\rho }+\frac{3}{8}\,{{C}_{{d}}}\,\mathbb{C}( \boldsymbol v ) +\frac{3{\dot{m}}}{\rho }\,{\boldsymbol v}_{m} -\frac{{{\rho }_{g}}}{\rho }\,R({\boldsymbol g} - {\dot{\boldsymbol v}}_{m} ) = {\boldsymbol 0}. \end{equation}

Equation (2.46) could be simplified to the form in our previous works (Zhang et al. Reference Zhang, Li, Cui, Li and Liu2023) if the last two terms on the left-hand side, representing the phase transition and inertia effect of the internal gas, are ignored.

2.4. Multiple-bubble interaction and boundary effects

In this subsection, we incorporate the effects of multiple bubbles and boundaries into the present theoretical model. The principle involves modifying the ambient pressure $P_{a}$ and velocity ${\boldsymbol u}_{a}$ of the background flow field of the bubble when the effect of multiple bubbles is considered. Also, accounting for boundary effects is achieved by introducing image bubbles, thereby transforming it into a multiple-bubble problem. First, we provide the pressure and velocity in the flow field induced by a single bubble. Differentiating (2.5) with respect to $\boldsymbol {r}$ and $t$ with the velocity potential of bubble migration ignored, and combining it with (2.17), we can establish the correlation between the physical information at $|\boldsymbol r|$ and the bubble surface:

(2.47) \begin{equation} \boldsymbol{u}{(\boldsymbol{r},t)} ={-}\frac{\boldsymbol{o} - \boldsymbol{r}}{|\boldsymbol{o} - \boldsymbol{r}|^3}\left. \left[\!R^2 \left(\!\dot{R}-\frac{\dot m}{\rho} \!\right)-\frac{R}{C}\,(|\boldsymbol{o} - \boldsymbol{r}|-R)\left(\!H+\frac12\left(\!\dot{R}-\frac{\dot m}{\rho}\!\right)^2+\frac14\,{v}^2\!\right) \!\right] \!\right|_{({R},t_{c})} \end{equation}

and

(2.48)\begin{equation} \frac{\partial \varphi{(\boldsymbol{r},t)}}{\partial t} ={-}\left. \frac{R}{|\boldsymbol{o} - \boldsymbol{r}|} \left( H + \frac12 \left(\dot{R}-\frac{\dot m}{\rho}\right)^2+ \frac14\,{v}^2\right) \right|_{({R},t_{c})}, \end{equation}

where $t_{c} = t - (|{\boldsymbol {r}}| - R)/C$, denoting the initiation moment of a disturbance induced by the bubble surface that later arrives at $\boldsymbol {r}$ at $t$. The flow pressure induced by the bubble can be solved by substituting the above two equations into the Bernoulli equation:

(2.49)\begin{equation} p ={-} \rho\,\frac{{\partial \varphi ( {{\boldsymbol r},t} )}}{{\partial t}} - \frac{1}{2}\,\rho\,|{\boldsymbol u}( {{\boldsymbol r},t} )|^2 + {{{P_\infty }}}. \end{equation}

Assuming that there are $U$ bubbles in the flow field, the velocity and pressure of the background field for bubble $N$ can be expressed as

(2.50)\begin{equation} \boldsymbol{u}_{{a}}{( {{{\boldsymbol{o}}_{N}},t} )}= \sum_{\genfrac{}{}{0pt}{2}{G=1,U}{G{\ne} N}} \boldsymbol{u}_{G}{(\boldsymbol{o}_{N},t)} \end{equation}

and

(2.51)\begin{equation} P_{{a}}{(\boldsymbol{o}_{N},t)} ={-}\rho \sum_{\genfrac{}{}{0pt}{2}{G=1,U}{G{\ne} N}} \frac{\partial {\varphi}_{G}{(\boldsymbol{o}_{N},t)}}{\partial t} -\frac12\,\rho\, |\boldsymbol{u}_{{a}}{(\boldsymbol{o}_{N},t)}|^2 + P_{\infty}, \end{equation}

respectively. The dynamics of multiple bubbles can be addressed by incorporating the above two equations into the oscillation equation and migration equation.

Further, assuming that a infinite flat boundary exists near the bubble, defined by ${\boldsymbol r} \boldsymbol {\cdot } {\boldsymbol {e}}_{{b}}+s=0$ (where ${\boldsymbol {e}}_{{b}}$ is the outward unit normal vector of the boundary plane, and $s$ is a constant), the position of the image bubble $N_{i}$ of bubble $N$ about the boundary could satisfy ${{\boldsymbol {o}}_{N_i}} = {{\boldsymbol {o}}_{N}} - 2 ( {{{\boldsymbol {o}}_{N}} \boldsymbol {\cdot } {{\boldsymbol {e}}_{{b}}} + s} ){{\boldsymbol {e}}_{{b}}}$. The size and oscillation velocity of bubbles $N$ and $N_{i}$ always remain exactly the same, while the position and migration velocity of the two bubbles are always symmetric about the boundary plane. A reflection coefficient $\xi$ is used to determine the property of the boundary. Specifically, $\xi = 1.0$ for a rigid boundary, and $\xi = -1$ for an ideal free surface. Therefore, when the bubble $N$ is affected by other bubbles and the boundary, the velocity and pressure of the flow field at $\boldsymbol {o}_{N}$ can be expressed as

(2.52)\begin{equation} \boldsymbol{u}_{{a}}{( {{{\boldsymbol{o}}_{N}},t} )}= \sum_{\genfrac{}{}{0pt}{2}{G=1,U}{G{\ne} N}} \boldsymbol{u}_{G}{(\boldsymbol{o}_{N},t)} + \xi \sum_{G=1,U} \boldsymbol{u}_{{ G_i}}{(\boldsymbol{o}_{N},t)}, \end{equation}

and

(2.53) \begin{align} P_{{a}}{(\boldsymbol{o}_{N},t)} ={-}\rho \sum_{\genfrac{}{}{0pt}{2}{G=1,U}{G{\ne} N}} \frac{\partial {\varphi}_{G}{(\boldsymbol{o}_{N},t)}}{\partial t} -\rho \xi \sum_{G=1,U} \frac{\partial {\varphi}_{{ G_i}}{(\boldsymbol{o}_{N},t)}}{\partial t} -\frac12 \rho |\boldsymbol{u}_{{a}}{(\boldsymbol{o}_{N},t)}|^2 + P_{\infty}, \end{align}

respectively.

The dynamics of bubble $N$ near the boundary and other bubbles can be addressed by incorporating (2.52) and (2.53) into the oscillation equation and migration equation of bubble $N$.

3.1. Bubble dynamics in the free field

First, we validate the present theoretical model through two cavitation bubble experiments in a free field. In the first experiment, the bubble is generated by laser focusing with maximum bubble radius 1.01 mm, and the experimental set-up can be referred to in the previous work (Li et al. Reference Li, Zhao, Zhang and Han2024). Figure 1(a) shows high-speed photography images in the first two oscillation cycles of the bubble, where the bubble remains nearly spherical during the first cycle, and undergoes slight deformation in the second cycle. Figure 1(b) presents a comparison between the computed bubble radius and the experimental data. The theoretical calculations start from the moment the bubble reaches its maximum radius, at which point the oscillation velocity of the bubble is zero. Since the initial conditions of bubbles in the experiments are difficult to determine, we discuss the effects of initial parameters here. The temperature of the fluid domain $T_{\infty }$ is fixed at 293 K in all the cases unless stated otherwise. The initial vapour proportion $M_{v}/M$ and the internal pressure $P_{g0}$ of the bubble are two important parameters that significantly affect the maximum radius of the bubble during the second cycle. Figure 1(c) shows the effect of the initial vapour proportion on the bubble radius at a fixed $P_{g0}$, and figure 1(d) shows the effect of the initial internal pressure at a fixed initial $M_{v}/M$. A higher initial vapour proportion and a lower initial internal pressure significantly reduce the maximum radius of the bubble during the second cycle. Considering the small content of non-condensable gases inside the laser bubble (Liang et al. Reference Liang, Linz, Freidank, Paltauf and Vogel2022), the initial vapour proportion is set at 1.0 in this case. The initial internal pressure of the bubble is set at 1.0 kPa to match the experimental bubble radius in the second cycle, and $\alpha _{m}$ in the Hertz–Knudsen–Langmuir relationship is chosen as 0.064. The number of air and vapour molecules at the initial moment is estimated by the ideal gas equation. The present theoretical model effectively captures the experimental bubble radius in the first two cycles. In addition, figure 1(b) also shows the calculation results without the phase transition and without the fluid compressibility, respectively. The results indicate that both phase transition and fluid compressibility are important factors affecting the energy loss of bubbles, with phase transition having a more significant impact on laser-induced bubbles. Figure 1(e) illustrates the temporal variation of the temperature of the bubble centre and the liquid temperature at the bubble surface. During the majority of the bubble cycle, the temperatures inside the bubble and at the bubble surface remain in close proximity. At the final stages of the bubble collapse, the temperature at the bubble centre rises more rapidly over time compared to that at the bubble surface. The discrepancy between the two temperatures manifests primarily during the intense phase of bubble collapse, where there is a stark rise from the ambient water temperature as the bubble collapses, followed by a precipitous drop towards equilibrium.

Figure 1. Laser bubble experiment in the free field and its comparison with theoretical results. (a) High-speed photography images of the bubble oscillation over time, with frame width 3.45 mm. (b) Comparison of the bubble radius between theory and experiment ($R_0 = 1.01$ mm, $P_{g0} = 1.0$ kPa, $M_{v}/M=1.0$, $\alpha _{m} = 0.064$). (c) Effect of initial proportion of vapour inside the bubble. (d) Effect of initial internal bubble pressure. (e) Temporal variations of the temperature at the bubble centre $T$ and the liquid temperature at the bubble surface $T_{l}$.

In the second experiment, the bubble is generated by underwater electrical discharge, reaching maximum bubble radius 18.1 mm. For details of the experiment method, refer to the work of Han et al. (Reference Han, Zhang, Tan and Li2022). Figure 2(a) shows the temporal evolution of the bubble shape. The bubble is accompanied by a large amount of flocculent impurities at the end of the second cycle, making it difficult to clearly observe the bubble profile. However, in general, the profile of the bubble at the moment of maximum volume is clear enough to accurately obtain the maximum radius of the bubble during the first two cycles, allowing for calculations using the present theoretical model. Figures 2(b) and 2(c) compare the bubble radius and the flow-field pressure induced by the bubble oscillation, respectively. The flow-field pressure is measured by a PCB free-field sensor placed 8.5 cm away from the bubble centre. Similar to figure 1, the initial vapour proportion of the bubble is set at 1.0. The initial internal pressure of the bubble is 12 kPa, and $\alpha _{m}=0.043$. The smaller peak of the flow-field pressure in the experiment may be attributed to the limited sampling frequency of the sensor or slight changes in the sensor's position under the influence of bubble oscillation. Overall, the theoretical calculations well reproduce the bubble radius and the oscillation pressure in the experiment. Neglecting the effect of phase transition and fluid compressibility significantly overestimates the radius of the bubble during the second cycle. Figure 2(d) shows the temporal variation of the gas temperature at the centre of the bubble and the liquid temperature at the surface of the bubble. Compared to the laser-induced bubble, the internal gas of the spark-generated bubble reaches a lower temperature at the moment of minimum volume, which to some extent indicates that the collapse intensity of the bubble is weaker, resulting in less energy loss of the bubble at the end of the first cycle. Note that the maximum bubble displacement of the bubble centre in the first two cycles of the cases in this section is small enough to ignore the bubble migration. Thus the migration features of the bubbles are not discussed here.

Figure 2. Spark-generated bubble experiment in the free field and its comparison with theoretical results. (a) High-speed photography images of the bubble oscillation over time, frame width 4.43 cm. (b) Comparison of the bubble radius between theory and experiment ($R_0=18.1$ mm, $P_{g0} = 12$ kPa, $M_{v}/M=1.0$, $\alpha _{m} = 0.043$). (c) Comparison of the flow-field pressure induced by bubble oscillation between theory and experiment. (d) Temporal variations of the temperature at the bubble centre $T$ and the liquid temperature at the bubble surface $T_{l}$.

3.2. Bubble dynamics under different boundary conditions

In this subsection, we validate the effects of boundaries and multiple bubbles in the bubble equation through two bubble experiments. The first case is an underwater explosion bubble generated by 1.05 g TNT explosives near the free surface, as shown in figure 3. The experiment is conducted in a cubic water tank, which can be referred to in the work of Zhang et al. (Reference Zhang, Li, Cui, Li and Liu2023). The underwater explosion bubble is initially $30$ cm from the free surface, with maximum radius $15.8$ cm. The initial condition of the bubble is calculated according to the shock wave theory in the previous works (Zhang et al. Reference Zhang, Li, Cui, Li and Liu2023): ${R_{0}}=0.026$ m, ${\dot {R}_{0}}=109\ {\rm m}\ {\rm s}^{-1}$, ${P_{g0}}=2.74$ MPa. The value of $\alpha _{m}$ is 0.041. In this case, we set the initial vapour proportion in the bubble content as $1\,\%$ in theory. This can be explained by the presence of a large amount of non-condensable contents inside the underwater explosion bubble. Meanwhile, we provide the calculated results without the effect of phase transition, showing that the phase transition plays a relatively minimal role for underwater explosion bubbles. The effects of the boundary and liquid compressibility play an important role in bubble dynamics. Removing the boundary effect in theory leads to a significantly larger bubble oscillation period and the wrong migration direction because the free surface could accelerate the bubble oscillation and induce the bubble to migrate downwards. As the boundary effect is neglected, the bubble migration is controlled by the buoyancy of the bubble. Removing the liquid compressibility results in a great deviation from the experimental values for both the maximum bubble radius and the energy loss of the bubble.

Figure 3. Experiment of an underwater explosion bubble near the free surface and its comparison with theoretical results. (a) High-speed photography images of the bubble oscillation over time, frame width 0.510 m. (b) Comparison of the bubble radius between the theoretical and experimental results (${R_{0}}=0.026$ m, ${\dot {R}_{0}}=109\ {\rm m}\ {\rm s}^{-1}$, ${P_{g0}}=2.74$ MPa, $M_{v}/M=0.01$, $\alpha _{m}=0.041$). (c) Comparisons of the bubble migration between the theoretical and experimental results.

The second case is a spark-generated bubble below the wall. Figure 4 provides the temporal progression of the bubble shape, accompanied by a comparison between the theoretical and experimental results of the bubble radius and vertical displacements. The bubble is 44 mm from the wall at inception, and the maximum bubble radius is 16.6 mm. The bubble does not have a tendency to migrate in most of the first cycle, but migrates upwards obviously at the end of collapse and the second cycle due to the effect of the wall. To compute the dynamics of a spark-generated bubble from its inception, the initial conditions of the bubble are obtained by integrating backwards from the moment of maximum bubble volume (Wang Reference Wang2013; Zhang et al. Reference Zhang, Li, Cui, Li and Liu2023). The detailed procedure is as follows. The computation begins at the moment of maximum bubble volume, when the bubble radius is known and its oscillation velocity is zero. Next, the internal pressure at the moment of maximum bubble volume depends on the experimental bubble radius in the second cycle. The calculation then proceeds in reverse along the time axis from the moment of maximum bubble volume until the computed time approaches zero. According to the method, the initial conditions of the bubble are ${R_{0}}=3.69$ mm, ${\dot {R}_{0}}=60\ {\rm m}\ {\rm s}^{-1}$, ${P_{g0}}=5.0$ MPa. The initial vapour proportion in theory and the value of $\alpha _{m}$ are the same as in the case of the spark-generated bubble in figure 2. The impact of various physical factors is also analysed in this case. Removing the boundary effect causes the bubble period to decrease due to the presence of the wall. The amplitude of the bubble migration is significantly weaker when the wall is neglected. The energy loss of the bubble is much weaker when the phase transition is not considered compared to the computational results without the fluid compressibility. As indicated by the previous underwater explosion experiments and this case, the relative impact of phase transition and fluid compressibility on the energy loss of bubbles is closely related to the composition of gases.

Figure 4. Spark-generated bubble experiment near a rigid wall, and its comparison with theoretical results. (a) High-speed photography images of the bubble oscillation over time, frame width 40 mm. (b) Comparison of the bubble radius between the theoretical and the experimental results (${R_{0}}=3.69$ mm, ${\dot {R}_{0}}=60\ {\rm m}\ {\rm s}^{-1}$, ${P_{g0}}=5.0$ MPa, $M_{v}/M=1.0$, $\alpha _{m}=0.043$). (c) Comparison of the bubble migration between the theoretical and experimental results.

3.3. Multiple bubble dynamics

Finally, we carry out an experiment with multiple spark-generated bubbles and compare it with the theoretical results, as shown in figure 5. Four bubbles are generated simultaneously at the four vertices of a square plane with side length 60 mm, as shown in figure 5(a). The maximum radius of the upper left and lower right bubbles is 10.0 mm, and that of the remaining two bubbles is 10.3 mm. We denote the upper left bubble as bubble 1, and the lower left bubble as bubble 2. The initial oscillation conditions for the two bubbles are $R_{01}=1.04$ mm, $R_{02}=1.06$ mm, $\dot {R}_{01}=\dot {R}_{02}=10\ {\rm m}\ {\rm s}^{-1}$ and $P_{g01}=P_{g02}=40$ MPa. Bubble 1 shares the same initial conditions as the lower right bubble, while the two other bubbles also have identical initial conditions. This is considering that the corresponding two bubbles remain symmetrical for the majority of the bubble oscillation cycle. The proportion of vapour in the bubble contents is 0.99, and $\alpha _{m}=0.041$. During the bubble oscillation, the four bubbles migrate towards the centre of the square plane due to the mutual attraction among bubbles. Figures 5(b) and 5(c) compare the radius and vertical displacement of the two bubbles, respectively. Overall, our theoretical model well reproduces the bubble radius and displacement in the experiment. It is observed that the experimental bubble radius at the end of the second cycle is larger than the computed values. This discrepancy may be attributed to the measurement errors caused by the frame rate of the high-speed camera and the perturbation of frothy impurities on the bubble.

Figure 5. Comparison of the interaction of four spark-generated bubbles with the theoretical results. (a) High-speed photography images of the bubble oscillation over time, frame width 162 mm. (b) Comparison of the bubble radius between the theoretical and the experimental results ($R_{01}=1.04$ mm, $R_{02}=1.06$ mm, $\dot {R}_{01}=\dot {R}_{02}=10\ {\rm m}\ {\rm s}^{-1}$, $P_{g01}=P_{g02}=40$ MPa, $M_{v}/M=0.99$, $\alpha _{m}=0.041$). (c) Comparison of the vertical bubble migration between the theoretical and experimental results.