A theoretical model for compressible bubble dynamics considering phase transition and migration

Abstract

A novel theoretical model for bubble dynamics is established that simultaneously accounts for the liquid compressibility, phase transition, oscillation, migration, ambient flow field, etc. The bubble dynamics equations are presented in a unified and concise mathematical form, with clear physical meanings and extensibility. The bubble oscillation equation can be simplified to the Keller–Miksis equation by neglecting the effects of phase transition and bubble migration. The present theoretical model effectively captures the experimental results for bubbles generated in free fields, near free surfaces, adjacent to rigid walls, and in the vicinity of other bubbles. Based on the present theory, we explore the effect of the bubble content by changing the vapour proportion inside the cavitation bubble for an initial high-pressure bubble. It is found that the energy loss of the bubble shows a consistent increase with increasing Mach number and initial vapour proportion. However, the radiated pressure peak by the bubble at the collapse stage increases with decreasing Mach number and increasing vapour proportion. The energy analyses of the bubble reveal that the presence of vapour inside the bubble not only directly contributes to the energy loss of the bubble through phase transition but also intensifies the bubble collapse, which leads to greater radiation of energy into the surrounding flow field due to the fluid compressibility.

1. Introduction

From the eruption of submarine volcanoes and underwater explosions (Klaseboer et al. 2005; Lyons et al. 2019; Wang et al. 2021) to the snapping of pistol shrimps (Versluis et al. 2000; Lohse, Schmitz & Versluis 2001), from targeted drug delivery and ultrasonic lithotripsy (Lokhandwalla et al. 2001; Ferrara, Pollard & Borden 2007; Maeda & Colonius 2019) to ultrasonic cleaning (Verhaagen & Rivas 2016; Oh et al. 2018; Landel & Wilson 2021), bubble dynamics holds significant importance across various academic areas and practical applications. The behaviour of oscillating bubbles involves a complex interplay of factors such as fluid compressibility, bubble migration, and mass and heat transfer (Fujikawa & Akamatsu 1980; Wang & Blake 2011; Brujan et al. 2022; Preso et al. 2024). Understanding these complex physical mechanisms not only advances fundamental knowledge but also drives innovations in technologies reliant on bubble phenomena.

The Rayleigh–Plesset equation (Rayleigh 1917; Plesset 1949) stands as a classical framework that is widely utilized to predict the oscillation behaviour of spherical cavitation bubbles. While rooted in the assumption of incompressible fluids, it has provided foundational insights into various aspects of bubble dynamics, including nonlinear bubble oscillations (Hicks 1970; Best 1991; Storey & Szeri 2001; Oratis et al. 2024), and linear interactions between multiple bubbles (Best 1991; Harkin, Kaper & Nadim 2001; Bremond et al. 2006). Over the years, researchers have developed numerous compressible models to address the limitations of the Rayleigh–Plesset equation (Herring 1941; Gilmore 1952; Keller & Kolodner 1956; Prosperetti & Lezzi 1986). Examples include the Keller–Miksis equation (Keller & Kolodner 1956; Keller & Miksis 1980), known for its robust theoretical foundation, and Ma, Hsiao & Chahine (2018) incorporated the influence of bubble migration by integrating an incompressible migration term. Geers & Hunter (2002) employed the doubly asymptotic approximation approach to develop equations that capture bubble oscillation and migration within a compressible flow environment for underwater explosion bubbles, and Zhang et al. (2023) derived the oscillation and migration equations in a compressible flow field under various environmental conditions based on the wave equation. However, it is difficult for their models to calculate the bubble dynamics in which the bubble contents are composed of both condensable and non-condensable gases.

Recent research by Zhong et al. (2020) and Han, Chen & Guo (2023) has revealed that in addition to fluid compressibility and viscosity, the condensation and evaporation processes of vapour inside laser-induced and spark bubbles also significantly affect the dynamic characteristics of bubbles, particularly concerning the issue of energy loss after the second cycle of bubble oscillation. Furthermore, previous compressible bubble models using the adiabatic gas equation of state have struggled to accurately reproduce the energy loss during the multi-cycle oscillation of bubbles, regardless of how the initial conditions are configured (Zeng et al. 2018; Cerbus et al. 2022; Fan et al. 2024). This suggests indirectly that these compressible bubble models lack certain crucial physical mechanisms. In fact, the phase transition model of bubbles has been studied extensively in previous works, including the state equation of gases (Abbondanza, Gallo & Casciola 2023; Gallo et al. 2023), the rate of phase transition (Fuster, Hauke & Dopazo 2010; Yasui 2018), and the temperature boundary layer near the bubble surface (Fujikawa & Akamatsu 1980; Hauke, Fuster & Dopazo 2007; Tian et al. 2022). Moreover, previous studies predominantly analysed the bubble migration under the assumption of incompressible fluids (Hicks 1970; Best 1991; Seo, Lele & Tryggvason 2010), overlooking the impact of fluid compressibility. In this study, we will specifically examine the roles of phase transition and bubble migration in formulating a comprehensive bubble oscillation equation within the compressible fluid domain, and a new migration equation that accounts for the effects of fluid compressibility and condensation/evaporation will be deduced.

Furthermore, quantitative analyses of bubble content remain challenging due to the complex mixture of water vapour and non-condensable gases within bubbles. To address this challenge, we manipulate the composition of bubble content by changing the initial vapour proportion inside the bubble at a constant initial internal bubble pressure, allowing for a systematic analysis of its influence on bubble dynamics. This study seeks to provide comprehensive insights into bubble dynamics, fostering advancements in both fundamental understanding and practical applications.

The structure of this paper is as follows. We first derive the theoretical model in detail in § 2, including the bubble oscillation equation, state equation of mixed gases, bubble migration equation, multiple bubble equation, and bubble equation of boundary effect. In § 3, the theoretical model is fully validated by several bubble experiments in the free field, near boundaries, and near multiple bubbles. In § 4, parametric studies on the effects of initial vapour proportion inside a bubble are conducted for an initially high-pressure bubble. Finally, this study is summarized and conclusions are made in § 5.

2. Theory

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. 2023). 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 1980), 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 (1980), 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 2002), 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 1997; Kyuichi 2021) 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 1953; Akhatov et al. 2001) 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 1998), 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. (2023), $\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 2014; Li et al. 2020), 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 1980; Yasui 1999). 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 (2016), 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 1969; Yasui 1995), 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 1980; Tian et al. 2022; Dai et al. 2024), 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 (1996) 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 1996, 1997) 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 1980) 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. (2020) and Yasui (2001), 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. 2014; Nagalingam et al. 2023; Chen et al. 2024), 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. 2020) 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. 2023) 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. Validation of the present theoretical model

In this section, we conduct experiments on bubbles with different sources and environmental conditions, capturing the bubble oscillation and migration processes. The experimental results are compared with the theoretical values to validate the present theoretical model.

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. 2024). 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. 2022), 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. (2022). 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. (2023). 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. 2023): ${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 2013; Zhang et al. 2023). 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.

4. Discussion on the energy loss of bubbles

In this section, we examine the influence of phase transition on the energy loss of bubbles. First, we present the distribution of feature parameters of bubbles across different Mach numbers $Ma =\sqrt {P_{g0}/\rho }/C$ and initial vapour proportions $M_{v}/M$, as depicted in figure 6. The Mach number serves to quantify the impact of fluid compressibility, and the vapour proportion represents the influence of phase transition. In this section, all physical quantities, with the exception of the temperature, are rendered dimensionless by using the maximum radius of the bubble $R_{max}$, the density of the liquid $\rho$, and the hydrostatic pressure at the bubble's initial location $P_{\infty }$. The dimensionless physical quantities are indicated by the superscript * in the latter descriptions. The studied characteristic parameters are the radius ratio during the first two cycles $R_{max2}/R_{max1}$, and the scaled internal bubble pressure $P^*=P^*_{max} r^*$ (where $P^*_{max}$ and $r^*$ are the peak pressure inside the bubble and the minimum bubble radius at the first collapse stage, respectively). Here, $R_{max2}/R_{max1}$ is used to measure the energy loss of the bubble in the first two cycles. Also, $P^*$ is utilized to characterize the energy at the moment of minimum bubble volume, and it equals the pressure peak induced by the bubble at the unit distance. The initial conditions of bubbles are ${{R}_{0}^*}=0.19$, ${\dot {R}_{0}^*}=0$, ${{P}_{g0}^*}=50$. The value of $\alpha _{m}$ is 0.041. The initial vapour proportion inside the bubble is altered under the condition of a constant initial internal pressure. In the calculations, the effect of bubble migration is removed in order to analyse the effect of vapour proportion and Mach number more accurately. Figure 6(a) illustrates the distribution of $R_{max2}/R_{max1}$ for varying $Ma$ and proportions of vapour $M_{v}/M$. The energy loss of bubbles shows a steady increase as the Mach number and vapour proportion increase. Note that when the vapour proportion approaches 1, the energy loss of bubbles is substantially greater compared to other cases. By examining the relationship between bubble energy and radius $E_1/E_2=(R_{max1}/R_{max2})^3$ (where $E_1$ and $E_2$ denote the bubble energies in the first and second cycles, respectively), it is observed that the energy loss in vapour bubbles is typically over 80 %, which is at least twice the energy loss observed in the bubbles formed by non-condensable gases at the same Mach number. Figure 6(b) shows the distribution of $P^*$ for varying $Ma$ and proportions of vapour $M_{v}/M$. Here, $P^*$ decreases as the Mach number increases because the fluid compressibility damps the pressure induced by the bubble oscillation. However, the increase in the vapour proportion results in a progressive increase in $P^*$, indicating that the vapour bubble can generate higher pressure peaks in the flow field compared to those composed purely of non-condensable gases.

Figure 6. Distribution of (a) radius ratio during the first two cycles of the bubble, and (b) scaled internal bubble pressure at the moment of minimum bubble volume for varying Mach numbers and vapour proportions (${{R}_{0}^*}=0.19$, ${\dot {R}_{0}^*}=0$, ${{P}_{g0}^*}=50$).

To elucidate the underlying mechanism behind the variation of $P^*$ with $M_{v}/M$ in figure 6, we conduct an energy analysis for a special case where the bubble contents consist purely of vapour ($M_{v}/M=1$). The formulas for calculating the internal energy $E_{i}$, potential energy $E_{p}$, kinetic energy $E_{k}$ and radiated acoustic energy $E_{a}$ of a bubble (Wang 2016; Li et al. 2020) are given below:

(4.1)$$\begin{gather} {E_{i}} = 3\,\frac{{{n_{t}}}}{{{N_{A}}}}\,{R_{g}}T, \end{gather}$$

(4.2)$$\begin{gather}{E_{p}} = {P_\infty }V, \end{gather}$$

(4.3) \begin{gather} {E_{k}} ={-}\frac{1}{2}\,\rho \int_S {{\varphi _{f}}\,\frac{{\partial {\varphi _{f}}}}{{\partial n}}\,{\rm d}S} = \frac{1}{2}\,\rho 4{\rm \pi} {R^2}\,\frac{{f( {t - R/C} )}}{R}\left( {\dot R - \frac{{\dot m}}{\rho }} \right)\nonumber\\ = 2\rho {\rm \pi}{R^2}\left( {\dot R - \frac{{\dot m}}{\rho }} \right)\left[ {R\left( {\dot R - \frac{{\dot m}}{\rho }} \right) - \frac{1}{C}\left( {\ddot R - \frac{{\ddot m}}{\rho }} \right){R^2} - \frac{{2R}}{C}{{\left( {\dot R - \frac{{\dot m}}{\rho }} \right)}^2}} \right], \end{gather}

(4.4)\begin{gather} {E_{a}} = \frac{\rho }{{4{\rm \pi} C}}\left[ { \int_0^t {{{\ddot V}^2}( t )\,{\rm d}t+\dot V( 0 )\,\ddot V( 0 ) - \dot V( t )\,\ddot V( t )} } \right], \end{gather}

where $\dot V$ and $\ddot V$ are the first- and second-order derivatives of bubble volume with respect to time, respectively, $\ddot m$ is the second-order derivative of $m$ with respect to time, and $f(t-R/C)$ in (4.3) is solved by conducting the perturbation method on (2.9).

The time evolutions for the mass and radius of the bubble are presented in figure 7(a), along with the radius of the bubble in the absence of phase transition or fluid compressibility. The initial conditions of the bubble are ${{R}_{0}^*}=0.19$, ${\dot {R}_{0}^*}=0$, ${{P}_{0}^*}=50$, $M_{v}/M=1$, $\alpha _{m}=0.041$. The condensation of vapour dominates the phase transition, leading to a consistent decreasing trend in the bubble mass. By comparing the results without considering the phase transition and without considering the compressibility of the fluid, it can be observed that the effect of the phase transition on the energy loss of the bubble is more pronounced. According to the relationship between bubble energy and radius, the energy loss of the bubble caused by the fluid compressibility ($E_2/E_1=0.87$) is less than $1/5$ that due to phase transition ($E_2/E_1=0.21$). Figure 7(b) shows the time histories of the internal energy of the bubble, the kinetic energy and potential energy of fluids induced by bubble oscillation, and the acoustic wave energy that propagates away. The internal energy of the bubble decreases with time in most of the cycles, except for a small increase near the moment of minimum bubble volume due to the increase in temperature inside the bubble. The changes in the potential and kinetic energies of fluids are closely related to the volume and oscillation velocity of the bubble, but their amplitudes decrease significantly in the second bubble cycle compared to the first bubble cycle. The sum of the internal energy, the potential and the kinetic energy characterizes the total energy of the bubble system, as shown by the purple line in figure 7(b). The total energy of the bubble system in the second cycle is reduced by approximately 4.7 near the moment of minimum bubble volume, while the acoustic wave energy radiated into the flow field during this period is approximately 2.4 (approximately half of the bubble energy loss). This contrasts with the trend suggested by the radius curves depicted in figure 7(a). This feature implies that the impact of vapour condensation on bubble energy loss is manifested not only through the decrease in bubble mass but also through the intensification of the bubble collapse. The intensification results in a greater propagation of energy into the flow field in the form of acoustic radiation, which explains the observed increase of the scaled internal pressure with higher vapour proportions.

Figure 7. Time evolutions of mass and energy of a vapour bubble (${{R}_{0}^*}=0.19$, ${\dot {R}_{0}^*}=0$, ${{P}_{0}^*}=50$, $M_{v}/M=1$, $\alpha _{m}=0.041$). (a) Mass and radius of the bubble. (b) Various energies of the system.

5. Conclusion

In this study, the theoretical model for bubble dynamics that accounts for the oscillation, migration, phase transition, fluid compressibility, boundary effect, multiple bubbles, viscosity and surface tension is derived. The oscillation and migration equations are characterized by a unified mathematical form, and the terms in the equations have clear physical meanings. The oscillation equation exhibits good extensibility and could be simplified to the classical Keller–Miksis equation after ignoring the effects of phase transition and bubble migration.

The present theoretical model is validated through comparisons of the theoretical results with experimental values of a laser bubble and a spark-generated bubble in the free field. The effects of the initial vapour proportion and internal pressure on the bubble dynamics are discussed as important parameters for comparing theoretical results with experimental values. Subsequently, the theoretical model is further validated by an underwater explosion bubble experiment and several spark-generated bubble experiments under different boundary conditions. For underwater explosion bubbles, the non-condensable gases are the principal constituents of the bubble contents, and the phase transition does not significantly affect the bubble dynamics; for laser bubbles and spark-generated bubbles, on the other hand, the phase transition is an important cause of bubble energy loss.

Based on the present theoretical model, the effects of the Mach number and the initial vapour proportion inside the bubble on the energy loss of the bubble are investigated for an initially high-pressure bubble. The energy loss of the bubble increases with increasing Mach number and initial vapour proportion. Specifically, the energy loss of a vapour bubble is more than twice that of a bubble composed purely of non-condensable gases. Also, the radiated pressure peak by the bubble increases with the increasing vapour proportion. The vapour not only causes a loss of bubble contents through condensation, but also leads to a more intense collapse. This intensified collapse, in turn, releases more acoustic energy into the surrounding fluid through pressure waves.