A compressible multi-scale model to simulate cavitating flows

Abstract

We propose a compressible multi-scale model that (i) captures the dynamics of both large vapour cavities (resolved vapour) and micro-bubbles (unresolved vapour), and (ii) accounts for medium compressibility. The vapour mass, momentum and energy in the compressible homogeneous mixture equations are explicitly decomposed into constituent resolved and unresolved components that are independently treated. The homogeneous mixture of liquid and resolved vapour is tracked as a continuum in an Eulerian sense. The unresolved vapour terms are expressed in terms of subgrid bubble velocities and radii that are tracked in a Lagrangian sense using a novel ‘$kR$-$RP$ equation’ (k, constant multiple; R, bubble size; RP, Rayleigh-Plesset). The $kR$-$RP$ equation is formally derived in terms of the pressure at a finite distance ($kR$) from the bubble while accounting for the effects of neighbouring bubbles; $p(kR)$ may therefore be either a near-field or far-field pressure. The equation exactly recovers the classical Rayleigh–Plesset and Keller–Miksis equations in the limits that $k$ and $c$ (speed of sound) become very large. Also, the results are independent of $k$ for a single bubble for all $k$, and for multiple bubbles when $kR < d$ (where $d$ denotes separation distance). Numerical results show this robustness of the model to the choice of $k$, which can be different for each bubble. The multi-scale model is validated for the collapse of a single resolved/unresolved bubble. Its ability to capture inter-bubble interactions is demonstrated for multiple bubbles exposed to an acoustic pulse. The model is then applied to a problem where resolved and unresolved bubbles co-exist. Finally, it is validated using a cluster of $1200$ bubbles exposed to a strong acoustic pulse. The results show the impact of the bubble cluster on the transmitted and reflected waves and the shielding effect where bubbles at the edge of the cluster shield the interior bubbles by dampening the incident acoustic wave.

1. Introduction

Cavitation refers to the phase change from liquid to vapour due to the liquid pressure dropping below the vapour pressure. It is commonly observed in turbomachinery, the wake of hydrofoils and propeller blades where it is often associated with noise and loss of efficiency. Violent bubble collapse during cavitation can even result in structural damage.

Cavitating flows possess vapour pockets of varying size, large variations in sound speed, turbulence and even shock waves. Accounting for all these features makes the numerical simulation of cavitating flows a very challenging task. As figure 1 shows, large vapour cavities can be captured on the computational grid while small bubbles are likely to be smaller than the grid elements. The mixture speed of sound reduces significantly in the presence of vapour/gas (Karplus 1957), and hence the presence of large cavities results in a highly compressible medium. Micro-bubble clusters can undergo violent collapse when subjected to a large external pressure, even generating shock waves (Reisman, Wang & Brennen 1998; Wang & Brennen 1999). Hence, models that can (i) accurately capture both resolved (large-scale vapour cavities) and unresolved phases (micro-bubbles), and (ii) account for the compressibility of the medium are necessary to capture the essential dynamics of a general cavitating flow.

Figure 1. A sketch showing a large vapour cavity and unresolved micro-bubbles found in a general cavitating flow. The right plot shows the magnified view of the cell where the unresolved bubbles lie.

The homogeneous mixture model (HMM) is one of the commonly used models for simulating cavitating flows (Bensow & Bark 2010; Gnanaskandan & Mahesh 2015; Asnaghi, Feymark & Bensow 2017; Schenke & van Terwisga 2017; Budich, Schmidt & Adams 2018). The HMM represents the mixture of liquid and vapour as a single continuum, and both are assumed to be in thermodynamic equilibrium. Past studies suggest the importance of medium compressibility. Bensow & Bark (2010) studied cavitation over a hydrofoil using incompressible HMM and observed a discrepancy in the cavity length that they attributed to the incompressible approximation. Gnanaskandan & Mahesh (2016) used compressible HMM to study the re-entrant jet mechanism during the sheet to cloud transition over a wedge and obtained good agreement with experiments. They also highlighted the role of medium compressibility in regulating the vorticity in the sheet cavity region. Budich et al. (2018) and Bhatt & Mahesh (2020) studied the same configuration in the periodic regime and observed a bubbly shock wave produced due to the cloud collapse downstream of the wedge. They found the bubbly shock to be locally supersonic, in agreement with experiments. Brandao, Bhatt & Mahesh (2020) studied the flow over a cylinder and concluded that the condensation shock is responsible for cavity collapse. Although HMM has successfully captured different physical mechanisms associated with large cavities, it is less accurate for vapour regions of the order of the computational cell size. For example, Bhatt & Mahesh (2020) have observed noticeable differences in mean vapour void fraction between HMM and experiment (Ganesh, Makiharju & Ceccio 2016) in the incipient regime for the flow over a wedge. Asnaghi, Feymark & Bensow (2018) observed a discrepancy in the cavity shedding location on a hydrofoil and concluded that finer computational meshes were needed to capture the shedding accurately. The HMM can indeed accurately predict the behaviour of vapour cavities (irrespective of their size) provided they are well resolved on the grid. This implies that regions with tiny bubbles need a very fine mesh that can lead to computationally expensive simulations.

Two popular approaches to track such tiny bubbles are Euler–Euler (EE) and Euler–Lagrange (EL). In the EE approach (Zhang & Prosperetti 1994, 1997; Pan, Dudukovic & Chang 1999; Park, Drew & Lahey 1999; Ando 2010), both phases are tracked in an Eulerian sense, and the governing equations are developed using the ensemble averaging method. The bubbles are grouped into bins in each computational cell based on their size distributions, and their statistics are computed. While it is cheaper than EL approach for monodisperse bubbles, it can become computationally expensive for a large number of polydisperse bubbles because the number of bins increases. Also, most of these models assume the carrier liquid to be incompressible, constraining them from being applied to highly compressible cavitating flows. In the EL approach (Seo, Lele & Tryggvason 2010; Fuster & Colonius 2011; Ma, Chahine & Hsiao 2015; Ghahramani, Arabnejad & Bensow 2019; Pakseresht & Apte 2019), the Navier–Stokes equations govern the ‘carrier’ liquid continuum, and each bubble is tracked in a Lagrangian sense using Newton's law of motion and the Rayleigh–Plesset (RP) equation. The RP equation (and most of its variants) assumes the bubble to be spherical and have spatially uniform properties. Most EL models developed for cavitation simulations assume the liquid to be incompressible (Ma et al. 2015; Ghahramani et al. 2019; Pakseresht & Apte 2019). For single-bubble collapse, Ghahramani et al. (2019) show the incompressible EL model to be accurate and allow for a higher time step and coarser grid compared with a fully resolved simulation. Fuster & Colonius (2011) developed an EL model where the compressible Navier–Stokes equations govern the carrier liquid. Maeda & Colonius (2018) used it to study burst wave lithotripsy. However, the liquid pressure is governed by an isentropic equation rather than a full equation of state. While such EL and EE formulations can accurately estimate the dynamics of a bubbly flow, they fare poorly when the large cavities either have an arbitrary shape or undergo asymmetric collapse.

Hybrid models aim to capture both large cavities and micro-bubbles and the transition between them. Hsiao, Ma & Chahine (2017) developed a hybrid model where the mixture of liquid and resolved vapour is modelled as an incompressible homogeneous mixture. This was coupled to a Lagrangian solver for the micro-bubbles. This hybrid model was applied to model sheet cavity dynamics and cloud shedding for a two-dimensional hydrofoil. Ghahramani, Strom & Bensow (2021) adopted a similar approach in developing a hybrid model. Their approach also accounts for effects such as bubble–bubble collision and bubble breakup. They found the hybrid model to perform better than the Eulerian approach in predicting cavitation inception in the wake of a bluff body. They acknowledge the need to account for the compressibility of the mixture medium while developing hybrid models for events such as cavitation erosion, where intense bubble-cloud collapse generates shock waves that damage the nearby surface material. A model that can capture both resolved and unresolved vapour accurately while accounting for medium compressibility does not exist, to the best of our knowledge. Developing such a multi-scale model is the first goal of this paper.

Another significant focus of this paper is the modelling of the bubble dynamics. The original RP equation (Rayleigh 1917; Plesset 1949) describes single-bubble behaviour in a quiescent incompressible liquid. Numerous modifications to this equation have been proposed to account for liquid compressibility (Gilmore 1952; Keller & Miksis 1980) and inter-bubble interactions (Seo et al. 2010; Fuster & Colonius 2011; Maiga, Coutier-Delgosha & Buisine 2018). In the original RP equation, the external pressure a bubble experiences is the pressure at infinity ($p_\infty$). Hsiao, Chahine & Liu (2000) modified $p_\infty$ to be the average liquid pressure close to the bubble surface. This modified version, the surface averaged pressure model, revealed scaling effects for cavitation inception, which the original RP equation could not show. They note that their proposed modification may be inaccurate for bubbles of the order of the cell size or larger. Seo et al. (2010) developed the locally volume-averaged Rayleigh–Plesset (LVARP) model where the volume average of local mixture pressure models local flow effects on the bubble. This volume averaging is performed over a region extending from the bubble surface ($r = R$) to $d/2$ (where $d$ is the inter-bubble separation distance). Its derivation requires a closed-form expression for pressure at a finite distance. Hence, accounting for medium compressibility becomes a difficult task. In addition, the dependence of the volume averaging process on $d$ makes it computationally expensive for non-uniform distribution of bubbles. Fuster & Colonius (2011) developed a compressible RP equation that accounts for inter-bubble interactions using potential flow theory. For a bubble, the impact of $N$ neighbouring bubbles is modelled explicitly in terms of their velocity potential. Obtaining the velocity potential of these $N$ bubbles, in turn, requires solving a system of $N$ linear equations resulting in $O(N^2)$ complexity. Ghahramani et al. (2019) developed an incompressible RP variant by integrating the incompressible spherical momentum equation from the bubble surface ($r = R$) to a finite distance $r = 2R$. The accuracy of pressure at $r = 2R$ being the external pressure has not been discussed for different bubble sizes and bubble separation distances. To summarize, the models mentioned either incorrectly assume the local pressure to be $p_{\infty }$ in the RP equation, do not account for medium compressibility or have a high computational cost associated with the explicit modelling of the inter-bubble interactions. In this paper, we derive a novel RP variant we term ‘$kR$-$RP$ equation’, that accounts for medium compressibility, local flow effect, inter-bubble interactions, and is derived in terms of the external pressure at an arbitrary finite distance from the bubble. The $kR$-$RP$ equation is shown to yield the classical incompressible RP and Keller–Miksis equations in the limits that the distance from the bubble centre and the speed of sound become infinitely large.

This paper proposes a novel multi-scale model that aims to simulate a wide range of complex problems such as: ($a$) sheet-to-cloud cavitating flows, ($b$) bubbly flows where dense bubble clusters are exposed to strong acoustic pulses and ($c$) interaction between micro-bubbles and large-scale vapour cavities. The multi-scale model in its current form does not account for the transition between resolved and unresolved vapour scales as well as phenomena such as bubble break-up, coalescence and deformation. The key idea in developing this model is to split the vapour mass, momentum and energy in the compressible homogeneous mixture equations into constituent resolved and unresolved vapour components. These components are treated differently in that the homogeneous mixture of liquid and the resolved vapour is tracked as a continuum entity in an Eulerian sense, while the unresolved vapour is tracked using the $kR$-$RP$ equation developed in this work. The $kR$-$RP$ equation is formally derived in terms of external pressure at a finite distance from the bubble to capture the local flow effects and inter-bubble interactions. Standard EL models do not perform such explicit decomposition; also they typically assume the cell pressure to be $p_\infty$ in the RP equation. The derivation of the multi-scale model and the $kR$-$RP$ equation and their key properties are discussed in §§ 2 and 3. Sections 4 and 5 summarize the multi-scale model and discuss the numerical method, respectively. In § 6, the model is validated for some benchmark problems. The model is then applied to a complex problem where a planar acoustic wave interacts with a cloud of $O(10^3)$ bubbles. A brief summary concludes the paper in § 7.

2. Derivation of the multi-scale model

The compressible Navier–Stokes equations for the homogeneous mixture of liquid and vapour are

(2.1)\begin{equation} \left.\begin{gathered} \frac{\partial{\rho}}{\partial{t}} + \frac{\partial{(\rho u_j)}}{\partial{x_j}} = 0, \\ \frac{\partial{(\rho u_i)}}{\partial{t}} + \frac{\partial{(\rho u_i u_j)}}{\partial{x_j}} =-\frac{\partial p}{\partial{x_i}} + \frac{\partial \sigma_{ij}}{\partial{x_j}}, \\ \frac{\partial{(\rho e_s)}}{\partial t} + \frac{\partial{(\rho e_s u_j)}}{\partial{x_j}} = \frac{\partial{Q_j}}{\partial{x_j}} -p\frac{\partial{u_j}}{\partial{x_j}} + \sigma_{ij}\frac{\partial{u_i}}{\partial{x_j}}, \end{gathered}\right\} \end{equation}

where $\rho$, $\rho u_i$, $\rho e_s$ and $p$ are the mixture density, velocity, internal energy and pressure, respectively; $\rho = \rho Y_l + \rho Y_v$, where $Y_l$ and $Y_v$ are the liquid and vapour mass fraction, respectively. The mixture pressure depends on the liquid and vapour mass ($p \equiv p(\rho Y_l, \rho Y_v)$). Hence, a governing equation is needed for the vapour mass to close the system of equations. The HMM uses the following transport equation to track the vapour mass:

(2.2)\begin{equation} \frac{\partial{(\rho Y_v)}}{\partial{t}} + \frac{\partial{(\rho Y_v u_j)}}{\partial{x_j}} = S_e - S_c, \end{equation}

where $S_e$ and $S_c$ are the empirical source terms responsible for the change in phase; $S_e$ causes the liquid to evaporate into vapour when its pressure drops below the saturation vapour pressure, and $S_c$ causes the vapour to condense into liquid when its pressure goes above the saturation vapour pressure. The accuracy of the transport equation depends on the resolution of the vapour region. Hence, HMM performs well for large vapour cavities and under-performs when only micro-bubbles are present. On the other hand, EL models use an RP equation (or its variant) to track the vapour mass. The RP equation (and most of its variants) assumes the bubble to be spherical and possess spatially uniform properties. Hence, the EL models accurately capture the behaviour of micro-bubbles but not massive cavities that are often non-spherical and have spatially varying properties. The implication is that the resolved and the unresolved vapour regions must be tracked differently to capture their dynamics precisely. Hence, the vapour mass is split into constituent resolved and unresolved components

(2.3)\begin{equation} \rho Y_v = \rho Y_{v_{res}} + \rho Y_{v_{un}}, \end{equation}

where $\rho Y_{v_{res}}$ and $\rho Y_{v_{un}}$ represent the resolved vapour mass and the unresolved vapour mass, respectively. Figure 1 schematically illustrates the vapour regions represented by $\rho Y_{v_{res}}$ and $\rho Y_{v_{un}}$. Such splitting enables the independent treatment of the resolved and unresolved vapour. Now $\rho Y_{v_{res}}$ is governed by the transport equation and $\rho Y_{v_{un}}$ is governed by an RP variant. By doing so, both HMM and EL models are brought together to develop the multi-scale model. The mixture density, momentum and internal energy can be written as

(2.4)\begin{equation} \left.\begin{gathered} \rho = \rho Y_l + \rho Y_{v_{res}} + \rho Y_{v_{un}} = \rho_l \alpha_l + \rho_{v_{res}} \alpha_{res} + \rho_{v_{un}} \alpha_{un}, \\ \rho u_i = \rho Y_l u_{l_i} + \rho Y_{v_{res}} u_{res_i} + \rho Y_{v_{un}} u_{un_i}, \\ \rho e_s = \rho Y_l e_{s_l} + \rho Y_{v_{res}} e_{s_{res}} + \rho Y_{v_{un}} e_{s_{un}}, \end{gathered}\right\} \end{equation}

where $\alpha$ is the volume fraction and $Y$ is the mass fraction. The subscripts $l$, $res$ and $un$ refer to the liquid, resolved vapour and unresolved vapour, respectively; $\rho$, $u$ and $e_s$ refer to the density, velocity and specific internal energy of the corresponding phases, respectively.

Note that only the liquid and resolved vapour are assumed to be in thermodynamic equilibrium, i.e. there is no temperature difference or slip velocity between these phases. This implies $u_l = u_{v_{res}} = u_{lr}$ and $T_l = T_{v_{res}} = T_{lr}$ where $T$ denotes the temperature. However, the unresolved vapour is not constrained to be in thermodynamic equilibrium with the liquid or resolved vapour. With these assumptions, the mixture momentum and energy in (2.4) can be rewritten as follows:

(2.5)\begin{equation} \left.\begin{gathered} \rho u_i = \rho Y_{lr}u_{lr_i} + \rho Y_{v_{un}} u_{un_i},\\ \rho e_s = \rho e_{s_{lr}} + \rho Y_{v_{un}} e_{s_{un}}, \end{gathered}\right\} \end{equation}

where $Y_{lr} = Y_l +Y_{v_{res}}$, $e_{s_{lr}} = Y_l e_{s_l} + Y_{v_{res}}e_{s_{res}}$. Substituting (2.4) and (2.5) in (2.1), we obtain the governing equations for the homogeneous mixture of liquid and resolved vapour

(2.6)\begin{equation} \left.\begin{gathered} \frac{\partial{(\rho Y_{l})}}{\partial{t}} + \frac{\partial{(\rho Y_{l} u_{lr_j})}}{\partial{x_j}} =-\frac{\partial{(\rho Y_{v_{un}})}}{\partial{t}} - u_{un_j}\frac{\partial{(\rho Y_{v_{un}})}}{\partial{x_j}} - \rho Y_{v_{un}} \frac{\partial{u_{un_j}}}{\partial{x_j}} - S_e + S_c, \\ \frac{\partial{(\rho Y_{lr}u_{lr_i})}}{\partial{t}} + \frac{\partial{(\rho Y_{lr} u_{lr_i} u_{lr_j})}}{\partial{x_j}} =-\frac{\partial p}{\partial{x_i}} + \frac{\partial \sigma_{ij}}{\partial{x_j}} - \frac{\partial{(\rho Y_{v_{un}} u_{un_i})}}{\partial{t}} - u_{un_j}\frac{\partial{(\rho Y_{v_{un}}u_{un_i})}}{\partial{x_j}} \\ \quad - \rho Y_{v_{un}}u_{un_i} \frac{\partial{u_{un_j}}}{\partial{x_j}}, \\ \frac{\partial{(\rho e_{s_{lr}})}}{\partial t} + \frac{\partial{(\rho e_{s_{lr}}u_{lr_j})}}{\partial{x_j}} = \frac{\partial{Q_j}}{\partial{x_j}} -p\frac{\partial{u_j}}{\partial{x_j}} + \sigma_{ij}\frac{\partial{u_i}}{\partial{x_j}} - \frac{\partial{(\rho Y_{v_{un}} e_{s_{un}})}}{\partial{t}}\\ \quad - u_{un_j}\frac{\partial{(\rho Y_{v_{un}}e_{s_{un}})}}{\partial{x_j}} - \rho Y_{v_{un}}e_{s_{un}} \frac{\partial{u_{un_j}}}{\partial{x_j}},\\ \frac{\partial{(\rho Y_{v_{res}})}}{\partial{t}} + \frac{\partial{(\rho Y_{v_{res}} u_{lr_j})}}{\partial{x_j}} = S_e - S_c, \end{gathered}\right\} \end{equation}

where $\sigma _{ij}$ and $Q_j$ are the viscous stress and heat flux of the mixture. The mixture pressure ($p$) is defined as

(2.7)\begin{equation} p = (1-\alpha_{un})p_{lr} + \sum_{k=1}^{N} \alpha_k p_k, \end{equation}

where $p_{lr}$ is the pressure of homogeneous mixture of the liquid and resolved vapour, $N$ is the number of unresolved bubbles and $p_{k}$ and $\alpha _k$ are the pressure and volume fraction of the $k\textrm {th}$ unresolved bubble, respectively; $p_{lr}$ is expressed in terms of the liquid and resolved vapour's equation of state, as shown below:

(2.8)\begin{equation} \left.\begin{gathered} p_{lr} = \rho_l \alpha_l' K_l T_{lr} \frac{p_{lr}}{p_{lr} + P_c} + \rho_{v_{res}} \alpha_{res}'R_vT_{lr}\\ \text{where} \quad \alpha_l' = \frac{\alpha_l}{1-\alpha_{un}} \quad \text{and} \quad \alpha_{res}' = \frac{\alpha_{res}}{1-\alpha_{un}}. \end{gathered}\right\} \end{equation}

Here, $R_v = 461.6\ \textrm {J}\ \textrm {kg}^{-1}\ \textrm {K}^{-1}$, $K_l = 2684.075\ \textrm {J}\ \textrm {kg}^{-1}\ \textrm {K}^{-1}$ and $P_c = 786.333 \times 10^6$ Pa are the constants associated with these equations of state. For an unresolved vapour bubble, its pressure ($p_k$) is the saturated vapour pressure that exclusively depends on the temperature. However, the pressure of an unresolved gas bubble is a function of both its temperature and density. A common approach to compute it is to assume the bubble behaviour to be isentropic. However, thermal damping can cause significant energy loss during the nonlinear oscillation of the bubbles. Prosperetti, Crum & Commander (1988) developed an equation for gas pressure that accounts for such heat losses in terms of the heat flux across the bubble surface. Preston, Colonius & Brennen (2007) developed a reduced-order model to compute the heat flux efficiently. The final equation that is used to compute the gas pressure in this paper is

(2.9a,b)\begin{equation} \frac{1}{R^{3\gamma}}\frac{{\rm d}(\,p_kR^{3\gamma})}{{\rm d}t} = \frac{3(\gamma-1)}{R^2} k_w \beta (T_k - T_0), \quad T_k = p_k/{\rho R_g}, \end{equation}

where $R$ and $\rho$ are the bubble's radius and density, respectively; $T_k$ is the bubble temperature (assumed to be uniform across the bubble), and $T_0$ is the liquid temperature at the bubble wall; $\gamma$, $k_w$, $R_g$ and $\beta$ are the constant parameters defined as the ratio of specific heats, thermal conductivity, gas constant and heat transfer coefficient of the gas, respectively. Their values are $\gamma = 1.4$, $R_g = 287.0\ \textrm {J}\ \textrm {kg}^{-1}\ \textrm {K}^{-1}$, $\beta = 5$ and $\kappa = 0.02479\ \textrm {W}\ \textrm {m}^{-1}\ \textrm {K}^{-1}$.

The heat flux is defined as

(2.10)\begin{equation} Q_j = (\alpha_l k_l + \alpha_{res} k_{res}) \frac{\partial T_{lr}}{\partial x_j} + \alpha_{un} k_{un} \frac{\partial T_{un}}{\partial x_j}, \end{equation}

where $k_l$, $k_{res}$ and $k_{un}$ are the thermal conductivities of the liquid, resolved vapour and unresolved vapour, respectively; $S_e$ and $S_c$ are the evaporation and condensation source terms for vapour, obtained from Saito et al. (2007). They are shown below:

(2.11)\begin{equation} \left.\begin{gathered} S_e = C_e \alpha_{res}^2 (1-\alpha_{res})^2\frac{\rho_l}{\rho_v}\frac{\max(\,p_v-p_{lr},0)}{\sqrt{2{\rm \pi} R_vT_{lr}}} \\ S_c = C_c \alpha_{res}^2 (1-\alpha_{res})^2\frac{\max(\,p_{lr} - p_v,0)}{\sqrt{2{\rm \pi} R_vT_{lr}}} \\ p_v = p_a \exp\left(\left(1-\frac{T_m}{T_{lr}}\right)(a+(b-cT_{lr})(T_{lr}-d)^2)\right) \end{gathered}\right\}, \end{equation}

where $C_e$ and $C_c$ are the empirical constants with units $\textrm {m}^{-1}$; $p_a = 22.130$ MPa, $T_m = 647.31$ K, $a = 7.21$, $b = 1.152 \times 10^{-5}$, $c = -4.787 \times 10^{-9}$ and $d = 483.16$. The resolved component of the mixture internal energy is defined as follows:

(2.12)\begin{equation} \left.\begin{gathered} \rho e_{s_{lr}} = \rho_l \alpha_l'e_l + \rho_{v_{res}}\alpha_{res}'e_{v_{res}}, \\ e_l = C_{v_l}T_{lr} + \frac{P_c}{\rho_l}, \quad e_{v_{res}} = C_{v_v}T_{lr}, \\ \rho e_{s_{lr}} = \frac{1}{1-\alpha_{un}}\left(\rho Y_{v_{res}}C_{v_{v}}T_{lr} + \rho Y_lC_{v_{l}}T_{lr} + \frac{\rho Y_lP_cK_lT_{lr}}{p_{lr} + P_c}\right), \end{gathered}\right\} \end{equation}

where $T_{lr}$ is the temperature of the homogeneous mixture of the liquid and resolved vapour; $C_{v_l}$ and $C_{v_v}$ are the specific heats at constant volume for the liquid and vapour, respectively. Note that the latent heat effects during the phase transfer have not been considered here. The pressure term on the right-hand side of the energy equation ($p({\partial {u_j}}/{\partial {x_j}})$) can be expressed as follows:

(2.13)\begin{equation} p\frac{\partial{u_j}}{\partial{x_j}} = (1-\alpha_{un})p_{lr}\frac{\partial{u_{lr_j}}}{\partial{x_j}} + \alpha_{un}p_{un}\frac{\partial{u_{un_j}}}{\partial{x_j}}. \end{equation}

Note that, if the resolved bubble contained gas instead of vapour, it can still be tracked using the transport equation. However, since gas does not undergo a phase change, the source terms will be inactive, i.e. $S_e = S_c = 0$.

2.1. Source terms

The unresolved vapour terms on the right-hand side of (2.6) act as source terms. The unresolved vapour mass ($\rho Y_{v_{un}}$), momentum ($\rho Y_{v_{un}} u_{un_i}$) and internal energy ($\rho e_{s_{un}}$) can be expressed in terms of the unresolved bubble properties as shown below

(2.14a–c)\begin{equation} \rho Y_{v_{un}} = \sum_{k=1}^{N}\rho_{v_k}\alpha_k, \quad \rho Y_{v_{un}}u_{un_i} = \sum_{k=1}^{N}\rho_{v_k}\alpha_ku_k, \quad \rho e_{s_{un}} = \sum_{k=1}^{N}\rho_{v_k} \alpha_k C_{v_v} T_k, \end{equation}

where $\alpha _k$, $\rho _{v_k}$, $u_k$ and $T_k$ are the volume fraction, density, translational velocity and temperature of the $k\textrm {th}$ unresolved bubble; $u_k$ may be obtained by interpolating the local Eulerian velocity field to the bubble's centre of mass (examples in this paper) or from a bubble momentum equation. Some of the terms have a divergence rate $({\partial {u_{un_j}}}/{\partial {x_j}})$ that is a measure of the expansion/collapse rate of the unresolved bubbles. Computing this divergence by differentiating the velocity field will not be accurate for unresolved bubbles. A better approach is to obtain ${\partial {u_{un_j}}}/{\partial {x_j}}$ in terms of the bubble quantities. The velocity divergence for a single bubble can be written as

(2.15)\begin{equation} \frac{\partial{u_{un_j}}}{\partial{x_j}} =-\frac{1}{\rho_v}\frac{{\rm D}\rho_v}{{\rm D}t} = \frac{1}{V}\frac{{\rm D}V}{{\rm D}t} = \frac{3\dot{R}}{R}, \end{equation}

where $\rho _v$, $V$ and $R$ are the density, volume and radius of the bubble, respectively. Summing over $N$ unresolved bubbles, the divergence term on the right-hand side of the liquid transport equation becomes

(2.16)\begin{equation} \rho Y_{v_{un}}\frac{\partial{u_{un_j}}}{\partial{x_j}} = \sum_{k=1}^{N}\frac{3\rho_{v_k}\alpha_k\dot{R_k}}{R_k}. \end{equation}

Similarly, the divergence terms on the right-hand side of the momentum and energy equation become

(2.17)\begin{equation} \left.\begin{gathered} \rho Y_{v_{un}}u_{un_i}\frac{\partial{u_{un_j}}}{\partial{x_j}} = \sum_{k=1}^{N}\frac{3\rho_{v_k}u_k\alpha_k\dot{R_k}}{R_k}, \quad \rho Y_{v_{un}}e_{s_{un}}\frac{\partial{u_{un_j}}}{\partial{x_j}} = \sum_{k=1}^{N}\frac{3\rho_{v_k}C_{v_v}T_k\alpha_k\dot{R_k}}{R_k}, \\ \alpha_{un}p_{un}\frac{\partial{u_{un_j}}}{\partial{x_j}} = \sum_{k=1}^{N}\frac{3p_k\alpha_k\dot{R_k}}{R_k}. \end{gathered}\right\} \end{equation}

Here, $p_k$ and $T_k$ are obtained from (2.9a,b). The computation of $\alpha _k$ is shown in the following section. The bubble size ($R_k$) and velocity ($\dot {R_k}$) are obtained from the novel $kR$-$RP$ equation derived in § 3.

2.2. Computing volume fraction ($\alpha _k$)

While computing $\alpha _{k}$ as the ratio of the volume of the bubble to the cell volume seems intuitive, it often gives rise to unstable solutions due to the discontinuous nature of $\alpha _{k}$ (Apte, Mahesh & Lundgren 2008; Raju et al. 2011). Using a kernel to smoothen the bubble volume distribution has been a popular approach in developing stable EL models. The commonly used Gaussian kernel has a standard distribution ($\sigma$) which depends exclusively on the cell size. Horne & Mahesh (2013) have shown that such kernels can result in unphysical volume fractions especially when bubble size is bigger and grows to become much larger than a computational cell. They proposed a new kernel function, which depends on both the bubble size and the computational cell size, and demonstrated its superior performance for bubbles bigger than the cell. The kernel function ($\,f$) is given by

(2.18)\begin{equation} f = \sum_{k=1}^{N} \frac{{\rm e}^{-{{r_k}^2}/{2\sigma^2}}}{2{\rm \pi}^{3/2}m^3}, \end{equation}

where $N$ is the total number of bubbles, $r_k$ is the distance between the cell centre and the $k{\textrm {th}}$ bubble and $\sigma$ (standard distribution) is defined as $\sigma = m({4{\rm \pi} }/{3})^{1/3}R_k$, where $m$ is a constant parameter and $R_k$ is the size of the $k{\textrm {th}}$ bubble. The volume fraction is then computed as follows:

(2.19)\begin{equation} \alpha_{un} = \frac{\displaystyle\int f \,{\rm d}V}{V_{cv}}, \end{equation}

where $V_{cv}$ is the volume of the cell.

3. Derivation of $kR$-$RP$ equation

The RP equation (and its several variants) are derived in terms of the ambient pressure at large distances ($p_\infty$); $p_\infty$ may be unambiguously defined for a single bubble or a cluster of bubbles with large inter-separation distances in a quiescent medium. However, cavitating flows often possess dense bubble clusters or bubbles in the vicinity of large vapour cavities. The local pressure and not $p_\infty$, determines the dynamics of the bubbles in such scenarios. Most simulations improperly use the standard RP equation with $p_\infty$ defined as the carrier fluid pressure interpolated to the bubble location. Here, we formally derive a novel $kR$-$RP$ equation in terms of the pressure at finite distance $kR$; $p(kR)$ may therefore be a near- or far-field pressure. The proposed equation is very attractive for simulations since the resolved mixture yields this pressure; the $kR$-$RP$ equation ensures that the effect on bubble radius is independent of $kR$.

The $kR$-$RP$ equation is derived using the spherical momentum equation (assuming the bubble to be spherical) along with the linear wave equation (to account for medium compressibility). The momentum and wave equations are

(3.1)\begin{equation} \left.\begin{gathered} \frac{\partial}{\partial r}\Bigg(\frac{\partial{\phi}}{\partial{t}} + \frac{1}{2}\left(\frac{\partial{\phi}}{\partial{r}}\right)^2\Bigg) =-\frac{1}{\rho}\frac{\partial{p}}{\partial{r}},\\ \frac{\partial^2{\phi}}{\partial{t^2}} = c^2 \Delta \phi ,\end{gathered}\right\} \end{equation}

where $\phi$ and $p$ are the velocity potential and pressure respectively. Instead of integrating the momentum equation from the bubble surface ($r = R$) to infinity ($r = \infty$), we integrate it from $r = R$ to $r = kR$ (where $k$ is a constant parameter) to account for the effect of local flow on the bubble dynamics. This yields

(3.2)\begin{equation} \phi_t(kR) - \phi_t(R) + \frac{1}{2}(\phi_r^2(kR) - \phi_r^2(R)) + \frac{p(kR) - p(R)}{\rho} = 0. \end{equation}

The general solution of the linear wave equation is

(3.3)\begin{equation} \left.\begin{gathered} \phi = \phi^{ou} + \phi^{in} \\ \text{where} \quad \phi^{ou} = \frac{1}{r}f(t - r/c), \quad \phi^{in} = \frac{1}{r}g(t + r/c). \end{gathered}\right\} \end{equation}

Here, $\phi ^{ou}$ is the potential of the outgoing wave generated by the bubble, and $\phi ^{in}$ is the net potential of the incoming waves generated by the neighbouring bubbles. Substituting (3.3) in (3.2) yields

(3.4)\begin{equation} \left.\begin{gathered} \phi^{ou}_t(kR) - \phi^{ou}_t(R) + \frac{1}{2}({\phi^{ou}_r}^2(kR) - {\phi^{ou}_r}^2(R)) = \frac{p(R) - p(kR)}{\rho} + I^{**}\\ \text{where}\quad I^{**} =- (\phi^{in}_t(kR) - \phi^{in}_t(R)) - ({\phi^{in}_r}(kR){\phi^{ou}_r(kR)} - {\phi^{in}_r(R)}{\phi^{ou}_r}(R))\\ - \frac{1}{2}({\phi^{in}_r}^2(kR)- {\phi^{in}_r}^2(R)). \end{gathered}\right\} \end{equation}

From (3.3), $\phi ^{in}_t$ and $\phi ^{ou}_t$ can be expressed in terms of $\phi ^{in}_r$ and $\phi ^{ou}_r$, respectively.

(3.5a,b)\begin{equation} \phi^{in}_t = \frac{c}{r}\phi^{in} + c \phi^{in}_r, \quad \phi^{ou}_t =-\frac{c}{r}\phi^{ou} - c \phi^{ou}_r. \end{equation}

The kinematic boundary condition at the bubble surface is $\phi _r(R) = \dot {R}$, i.e. $\phi _r^{in}(R) + \phi _r^{ou}(R) = \dot {R}$. Also, it is reasonable to assume that this velocity is mainly due to that bubble itself. Hence, $\phi ^{in}_r(R) \approx 0$ and $\phi ^{ou}_r(R) \approx \phi _r(R) = \dot {R}$. Substituting the kinematic boundary condition and (3.5a,b) in (3.4), and then taking its temporal derivative yields

(3.6)\begin{align} &R\ddot{R}\left(1-\frac{\dot{R}}{c}\right) + \dot{R}^2\left(1-\frac{\dot{R}}{2c}\right) -\frac{1}{k}\frac{{\rm d}}{{\rm d}t}\phi^{ou}(kR) - \left(\dot{R} + R\frac{{\rm d}}{{\rm d}t}\right)\phi^{ou}_r(kR)\nonumber\\ &\quad + \frac{1}{2}\left(\frac{\dot{R}}{c} + \frac{R}{c}\frac{{\rm d}}{{\rm d}t}\right){\phi^{ou}_r}^2(kR) + \frac{{\rm d}}{{\rm d}t}\phi^{ou}(R) = \frac{1}{\rho}\left(\frac{\dot{R}}{c} + \frac{R}{c}\frac{{\rm d}}{{\rm d}t}\right) \left(\,p(R) - p(kR)\right) + \frac{{\rm d}I}{{\rm d}t}^{**}. \end{align}

Also, $p(R) = p_b - 2\sigma /R - 4\mu \dot {R}/R$, where $p_b$ is the bubble pressure, $\sigma$ is the surface tension of the bubble and $\mu$ is the viscosity of the liquid. Substituting $p(R)$ in (3.6) yields (please refer to Appendix A for the details)

(3.7)\begin{align} \left.\begin{gathered} R\ddot{R} \left(1-\frac{\dot{R}}{c} \right) + \frac{3}{2}{\dot{R}}^2\left(1-\frac{\dot{R}}{3c}\right) = \frac{1}{\rho} \left(1+\frac{\dot{R}}{c}+\frac{R}{c}\frac{\partial}{\partial t}\right)\left(\,p_b - \frac{2\sigma}{R} - 4\mu\frac{\dot{R}}{R} - p(kR)\right)\\ - \frac{1}{2}\left(1+\frac{\dot{R}}{c}+\frac{R}{c}\frac{\partial}{\partial t}\right){\phi_r^{ou}}^2(kR) - \phi_t^{ou}(kR) - \frac{R(k\dot{R}-c)}{c^2}{\phi^{ou}_{tt}}(kR) + I \\ \text{where}\quad I = \frac{{\rm d}I}{{\rm d}t}^{**} =- \frac{1}{2}\left(1+\frac{\dot{R}}{c}+ \frac{R}{c}\frac{\partial}{\partial t}\right)\left({\phi_r^{in}}^2(kR) + 2\phi_r^{in}(kR)\phi_r^{ou}(kR)\right) \\ -\phi_t^{in}(kR) - \frac{R(k\dot{R}+c)}{c^2}\phi^{in}_{tt}(kR) + \phi^{in}_t(R). \end{gathered}\right\} \end{align}

Equation (3.7) is referred to as the $kR$-$RP$ equation. Two important terms in this equation are $p(kR)$ and $I$; $p(kR)$, the pressure of the resolved field at a finite distance, is defined as follows:

(3.8)\begin{equation} p(kR) = \frac{\displaystyle\sum\nolimits_{i=1}^{N_{cell}} p_i(r=kR)}{N_{cell}}, \end{equation}

where $N_{cell}$ is the number of cells which are at a distance $r = kR$ from the bubble centre and $p_i(r=kR)$ is the pressure of the resolved phase ($p_{lr}$) in the $i{\textrm {th}}$ cell. Hence, $p(kR)$ can be either near-field or far-field pressure depending on the value of $k$; $I$ represents the variations in the velocity potential and kinetic energy of the neighbouring bubbles. While $p(kR)$ captures the local flow effects, both $p(kR)$ and $I$ together contribute to the inter-bubble interactions.

3.1. Closure to the $kR$-$RP$ equation

The unknown terms in (3.7) are $\phi ^{ou}(kR)$, $\phi ^{in}(kR)$ and their derivatives. Since the goal is to capture the local flow effects, regions in the bubble vicinity are considered and a small $k$ represents such regions. The flow can be assumed to be incompressible in the vicinity of the bubble as the time taken for sound to travel is very small relative to the bubble timescale (acoustic compactness). Hence, the potential of the outgoing wave can be approximated as

(3.9)\begin{equation} \phi^{ou}(r) \approx-\frac{\dot{R}R^2}{r}. \end{equation}

Therefore

(3.10a–c)\begin{align} \phi^{ou}_r(kR) \approx \frac{\dot{R}}{k^2}; \quad \phi^{ou}_t(kR) \approx-\frac{2R\dot{R}^2 + \ddot{R}R^2}{kR}; \quad \phi^{ou}_{tt}(kR) \approx-\frac{2\dot{R}^3 + 6R\dot{R}\ddot{R} + R^2\dddot{R}}{kR}. \end{align}

The net potential of the incoming waves is trickier to compute because of its non-uniform nature. Consider the sketch shown in figure 2 where the bubble $i$ has only one neighbour (bubble $j$) in its vicinity. Let $\phi _j(P)$ denote the potential of the $j\textrm {th}$ bubble at point $P$ on the $kR$ surface of the $i\textrm {th}$ bubble. It can be approximated under the incompressible assumption as

(3.11)\begin{equation} \phi_j(P) =-\frac{\dot{R_j}R_j^2}{\sqrt{d_{ij}^2 + k^2R^2 -2kRd_{ij}\cos\theta}}, \end{equation}

where $\theta$ is the orientation of $P$ with respect to the coordinate system attached to the $i^{th}$ bubble. Note that $\phi _j$ varies with $\theta$, i.e. $\phi _j \equiv \phi _j(kR,\theta )$. Hence, it needs to be integrated over the $kR$ surface to obtain an average value (please refer to Appendix B for details of the integration). The final expression is shown below

(3.12)\begin{equation} \phi^{in}(kR) = \frac{R_j^2\dot{R_j}}{2kRd}(d_{ij} + kR - |d_{ij} - kR|). \end{equation}

For $N$ neighbouring bubbles in the vicinity, the net potential of the incoming wave is

(3.13)\begin{equation} \phi^{in}(kR) =-\sum_{\substack{j=1}}^{N} \frac{R_j^2\dot{R_j}}{2kRd_{ij}}(d_{ij} + kR - |d_{ij} - kR|). \end{equation}

The other terms such as $\phi ^{in}_t(kR)$, $\phi ^{ou}_r(kR)\phi ^{in}_r(kR)$, ${\phi ^{in}_r}^2(kR)$ and $\phi ^{in}_{tt}(kR)$ are obtained in a similar way. Their final expressions are given below:

(3.14)\begin{align} \left.\begin{gathered} \phi^{in}_t(kR) =-\sum_{\substack{j=1}}^{N}\frac{2R_j{\dot{R}_j^2} + R_j^2\ddot{R_j}}{2kRd_{ij}}\left(d_{ij} + kR - |d_{ij} - kR|\right), \\ \phi^{in}_{tt}(kR) =-\sum_{\substack{j=1}}^{N}\frac{2\dot{R}_j^3 + 6R_j\dot{R}_j\ddot{R}_j + R_j^2\dddot{R_j}}{2kRd_{ij}}(d_{ij} + kR - |d_{ij} - kR|), \\ \phi^{ou}_r(kR)\phi^{in}_r(kR) = \sum_{\substack{j=1}}^{N}\frac{R_j^2\dot{R_j}\dot{R}}{2d_{ij}Rk^3}(1 + \frac{d_{ij}-kR}{|d_{ij}-kR|} - \frac{1}{kR}(d_{ij} + kR - |d_{ij} - kR|)), \\ {\phi^{in}_r}^2(kR) = \sum_{\substack{j=1}}^{N} \sum_{\substack{k=1 \\ k \neq j}}^{N} \frac{R_j^2\dot{R_j}R_k^2\dot{R_k}}{4d_{ij}d_{ik}k^2R^2_i}\left(1 + \frac{d_{ij}-kR}{|d_{ij}-kR|} - \frac{1}{kR}\left(d_{ij} + kR - |d_{ij} - kR|\right)\right) \\ \quad \left(1 + \frac{d_{ik}-kR}{|d_{ik}-kR|} - \frac{1}{kR}\left(d_{ik} + kR - |d_{ik} - kR|\right)\right) + \sum_{\substack{j=1}}^{N} \frac{(R_j^2\dot{R_j})^2}{(d_{ij}+kR)^2(d_{ij}-kR)^2},\\ \phi^{in}_t(R) =-\sum_{\substack{j=1}}^{N}\frac{2R_j{\dot{R}_j^2} + R_j^2\ddot{R_j}}{d_{ij}} \end{gathered}\right\}. \end{align}

Substituting (3.14) in (3.7) gives a closed-form expression for $I$. It can be further simplified based on the sign of ($d_{ij} - kR$) as shown below

(3.15)\begin{align} \left.\begin{gathered} I = \sum_{\substack{j=1}}^{N} I^*_j \\ I^*_j =- \frac{1}{2}\left(1+\frac{\dot{R}}{c}+\frac{R}{c}\frac{\partial}{\partial t}\right) \left(- \frac{R_j^2\dot{R_j}\dot{R}}{k^4R^2} + \frac{(R_j^2\dot{R_j})^2}{(d_{ij}+kR)^2(d_{ij} -kR)^2} + \sum_{\substack{k=1 \\ k \neq j}}^{N} \frac{R_j^2\dot{R_j}R_k^2\dot{R_k}}{k^4R^4}\right) \\ \quad + (2R_j{\dot{R}_j^2} + R_j^2\ddot{R_j})\left(\frac{1}{kR} - \frac{1}{d_{ij}}\right) + \frac{2\dot{R}_j^3 + 6R_j\dot{R}_j\ddot{R}_j + R_j^2\dddot{R_j}}{kc} \quad \text{if } d_{ij} < kR\\ = \frac{(R_j^2\dot{R_j})^2}{(d_{ij}+kR)^2(d_{ij}-kR)^2} + \frac{R(2\dot{R}_j^3 + 6R_j\dot{R}_j\ddot{R}_j + R_j^2\dddot{R_j})}{kcd_{ij}} \qquad\quad \text{if } d_{ij} > kR \end{gathered}\right\}, \end{align}

where $I_j^*$ represents the impact of the $j{\textrm {th}}$ neighbouring bubble. The $kR$-$RP$ equation, upon substituting (3.10a–c) in (3.7) and neglecting $c^{-2}$ terms, becomes

(3.16)\begin{align} &R\ddot{R}\left(1-\frac{1}{k}-\frac{\dot{R}}{c}\left(1+\frac{6}{k}-\frac{1}{k^4}\right)\right) + \frac{3}{2}{\dot{R}^2}\left(1-\frac{4}{3k}-\frac{1}{3k^4}-\frac{\dot{R}}{3c}\left(1+\frac{4}{k}-\frac{1}{k^4}\right)\right)\nonumber\\ & \quad = \frac{1}{\rho} \left(1+\frac{\dot{R}}{c}+\frac{R}{c}\frac{\partial}{\partial t}\right)(\,p_b - 2\sigma/R - 4\mu\dot{R}/R - p(kR_i)) + \frac{ R^2\dddot{R}}{kc} + I, \end{align}

where $I$ is given by (3.15).

Figure 2. A sketch of two bubbles that gives a physical sense of the non-uniform effect of bubble $j$ on the $kR$ surface of bubble $i$.

3.2. Significance of the interaction term $I$

If the inter-bubble separation distance is much larger than the bubble radius (i.e. $d_{ij} > R$, $d_{ij} > R$) and $|{d_{ij}}/{kR}| \gg 1$, then

(3.17)\begin{equation} \frac{(R_j^2\dot{R_j})^2}{(d_{ij}+kR)^2(d_{ij}-kR)^2} \sim 0. \end{equation}

Also, $c^{-1}$ terms can be neglected for gentle bubble oscillations. Thus, $I$ would become

(3.18)\begin{equation} \left.\begin{gathered} I = \sum_{\substack{j=1}}^{N}\left(2R_j{\dot{R}_j^2} + R_j^2\ddot{R_j}\right)\left(\frac{1}{kR} - \frac{1}{d_{ij}}\right) \\ \quad - \sum_{\substack{j=1}}^{N}\frac{R_j^2\dot{R_j}\dot{R}}{k^4R^2} + \sum_{\substack{j=1}}^{N} \sum_{\substack{k=1 \\ k \neq j}}^{N} \frac{R_j^2\dot{R_j}R_k^2\dot{R_k}}{k^4R^4} \quad \text{if } d_{ij} < kR \\ = 0 \quad \text{if } d_{ij} > kR. \end{gathered}\right\} \end{equation}

The important implication is that, for large inter-bubble separation distances, $p(kR)$ alone can capture the inter-bubble interactions when $kR < d_{ij}$. Else, both $p(kR)$ and $I$ contribute to the inter-bubble interactions.

3.3. Special cases of the $kR$-$RP$ equation

Well-known RP variants such as the incompressible RP equation, Keller–Miksis equation and the RP equation with inter-bubble interaction can be obtained from the $kR$-$RP$ equation with appropriate assumptions as shown below.

Case $1$: $k \to \infty$

For large $k$

(3.19)\begin{equation} I \sim-\sum_{\substack{j=1}}^{N}\frac{2R_j{\dot{R}_j^2} + R_j^2\ddot{R_j}}{d_{ij}}. \end{equation}

Hence, (3.16) becomes

(3.20)\begin{align} R\ddot{R}\left(1-\frac{\dot{R}}{c}\right) + \frac{3}{2}{\dot{R}^2}\left(1-\frac{\dot{R}}{3c}\right) &= \frac{1}{\rho}\left(1+\frac{\dot{R}}{c}+\frac{R}{c}\frac{\partial}{\partial t}\right)\left(\,p_b - \frac{2\sigma}{R} - 4\mu\frac{\dot{R}}{R} - p_\infty\right)\nonumber\\ & \quad -\sum_{\substack{j=1}}^{N}\frac{2R_j{\dot{R}_j^2} + R_j^2\ddot{R_j}}{d_{ij}}. \end{align}

This equation has been widely used to study the primary and secondary Bjerknes forces (Mettin et al. 1997; Doinikov 2001) during the interaction between an acoustic pulse and a pair of bubbles. If $N = 1$, $I = 0$. Equation (3.20) then reduces to

(3.21)\begin{align} R\ddot{R}\left(1-\frac{\dot{R}}{c}\right) + \frac{3}{2}{\dot{R}}^2\left(1-\frac{\dot{R}}{3c}\right) = \frac{1}{\rho} \left(1+\frac{\dot{R}}{c}+\frac{R}{c}\frac{\partial}{\partial t}\right)\left(\,p_b - \frac{2\sigma}{R} - 4\mu\frac{\dot{R}}{R} - p_\infty\right). \end{align}

This is the well-known Keller–Miksis equation for a single bubble that accounts for medium compressibility.

Case $2$: $c \to \infty$

For large $c$, all the $c^{-1}$ terms will become $0$. Hence (3.16) becomes

(3.22)\begin{align} R\ddot{R}\left(1-\frac{1}{k}\right) + \frac{3}{2}{\dot{R}^2}\left(1-\frac{4}{3k}-\frac{1}{3k^4}\right) = \frac{1}{\rho}\left(\,p_b - 2\sigma/R - 4\mu\dot{R}/R_i-p(kR)\right) + I, \end{align}

where

(3.23)\begin{align} I &= \sum_{\substack{j=1}}^{N} I^*_j \end{align}

(3.24)\begin{align} I^*_j &=- \frac{1}{2}\left(- \frac{R_j^2\dot{R_j}\dot{R}}{k^4R^2} + \frac{(R_j^2\dot{R_j})^2}{(d_{ij}+kR)^2(d_{ij} -kR)^2} + \sum_{\substack{k=1 \\ k \neq j}}^{N} \frac{R_j^2\dot{R_j}R_k^2\dot{R_k}}{k^4R^4}\right) \nonumber\\ &\quad + (2R_j{\dot{R}_j^2} + R_j^2\ddot{R_j})\left(\frac{1}{kR} - \frac{1}{d_{ij}}\right) \quad \text{if } d_{ij} < kR \nonumber\\ & = \frac{(R_j^2\dot{R_j})^2}{(d_{ij}+kR)^2(d_{ij}-kR)^2} \qquad\qquad\enspace\text{if } d_{ij} > kR. \end{align}

This is the incompressible version of the $kR$-$RP$ equation.

Case $3$: $c \to \infty$ and $k \to \infty$

For large $k$ and $c$, (3.16) becomes

(3.25)\begin{equation} R\ddot{R} + \frac{3}{2}{\dot{R}^2} = \frac{1}{\rho}\left(\,p_b - 2\sigma/R - 4\mu\dot{R}/R - p_\infty \right) -\sum_{\substack{j=1}}^{N}\frac{2R_j{\dot{R}_j^2} + R_j^2\ddot{R_j}}{d_{ij}}. \end{equation}

This is the incompressible RP equation with inter-bubble interaction terms (Bremond et al. 2006). For $N = 1$, (3.25) would reduce to the classic incompressible RP equation.

4. Summary of the model

To summarize the multi-scale model, the governing equations for the homogeneous mixture of liquid and resolved vapour are

(4.1)\begin{align} \left.\begin{gathered} \frac{\partial{(\rho Y_{l})}}{\partial{t}} + \frac{\partial{(\rho Y_{l} u_{lr_j})}}{\partial{x_j}} =-\frac{\partial{(\rho Y_{v_{un}})}}{\partial{t}} - u_{un_j}\frac{\partial{(\rho Y_{v_{un}})}}{\partial{x_j}} - \rho Y_{v_{un}} \frac{\partial{u_{un_j}}}{\partial{x_j}} -S_e + S_c, \\ \frac{\partial{(\rho Y_{lr}u_{lr_i})}}{\partial{t}} + \frac{\partial{(\rho Y_{lr} u_{lr_i} u_{lr_j})}}{\partial{x_j}} =-\frac{\partial p}{\partial{x_i}} + \frac{\partial \sigma_{ij}}{\partial{x_j}} - \frac{\partial{(\rho Y_{v_{un}} u_{un_i})}}{\partial{t}} - u_{un_j}\frac{\partial{(\rho Y_{v_{un}}u_{un_i})}}{\partial{x_j}}, \\ \frac{\partial{(\rho e_{s_{lr}})}}{\partial t} + \frac{\partial{(\rho e_{s_{lr}}u_{lr_j})}}{\partial{x_j}} = \frac{\partial{Q_j}}{\partial{x_j}} -p\frac{\partial{u_j}}{\partial{x_j}} + \sigma_{ij}\frac{\partial{u_i}}{\partial{x_j}} - \frac{\partial{(\rho Y_{v_{un}} e_{s_{un}})}}{\partial{t}} \\ \quad -u_j\frac{\partial{(\rho Y_{v_{un}}e_{s_{un}})}}{\partial{x_j}} - \rho Y_{v_{un}}e_{s_{un}} \frac{\partial{u_j}}{\partial{x_j}}, \\ \frac{\partial{(\rho Y_{v_{res}})}}{\partial{t}} + \frac{\partial{(\rho Y_{v_{res}} u_{lr_j})}}{\partial{x_j}} = S_e - S_c, \end{gathered}\right\} \end{align}

where the terms on the right-hand side with divergence rate ($\partial {u_{un_j}}/\partial {x_j}$) are modelled in terms of unresolved bubble properties

(4.2)\begin{equation} \left.\begin{gathered} \rho Y_{v_{un}}\frac{\partial{u_{un_j}}}{\partial{x_j}} = \sum_{k=1}^{N}\frac{3\rho_{v_k}\alpha_k\dot{R_k}}{R_k}, \quad \rho Y_{v_{un}}u_{un_i}\frac{\partial{u_{un_j}}}{\partial{x_j}} = \sum_{k=1}^{N}\frac{3\rho_{v_k}u_k\alpha_k\dot{R_k}}{R_k}, \\ \rho Y_{v_{un}}e_{s_{un}}\frac{\partial{u_{un_j}}}{\partial{x_j}} = \sum_{k=1}^{N}\frac{3\rho_{v_k}C_{v_v}T_k\alpha_k\dot{R_k}}{R_k}, \quad \alpha_{un}p_{un}\frac{\partial{u_{un_j}}}{\partial{x_j}} = \sum_{k=1}^{N}\frac{3p_k\alpha_k\dot{R_k}}{R_k}. \end{gathered}\right\} \end{equation}

The dynamics of the unresolved bubbles is governed by the following $kR$-$RP$ equation:

(4.3)\begin{align} &R\ddot{R}\left(1-\frac{1}{k}-\frac{\dot{R}}{c}\left(1+\frac{6}{k}-\frac{1}{k^4}\right)\right) + \frac{3}{2}{\dot{R}^2}\left(1-\frac{4}{3k}-\frac{1}{3k^4}-\frac{\dot{R}}{3c}\left(1+\frac{4}{k}-\frac{1}{k^4}\right)\right) \nonumber\\ & \quad = \frac{1}{\rho} \left(1+\frac{\dot{R}}{c}+\frac{R}{c}\frac{\partial}{\partial t}\right)(\,p_b - 2\sigma/R - 4\mu\dot{R}/R - p(kR_i)) + \frac{ R^2\dddot{R}}{kc} + I \end{align}

where

(4.4)\begin{equation} I = \sum_{\substack{j=1}}^{N} I^*_j \end{equation}

and

(4.5)\begin{align} I^*_j &=- \frac{1}{2}\left(1+\frac{\dot{R}}{c}+\frac{R}{c}\frac{\partial}{\partial t}\right) \left(- \frac{R_j^2\dot{R_j}\dot{R}}{k^4R^2} + \frac{(R_j^2\dot{R_j})^2}{(d_{ij}+kR)^2(d_{ij} -kR)^2} + \sum_{\substack{k=1 \\ k \neq j}}^{N} \frac{R_j^2\dot{R_j}R_k^2\dot{R_k}}{k^4R^4}\right) \nonumber\\ & \qquad + (2R_j{\dot{R}_j^2} + R_j^2\ddot{R_j})\left(\frac{1}{kR} - \frac{1}{d_{ij}}\right) + \frac{2\dot{R}_j^3 + 6R_j\dot{R}_j\ddot{R}_j + R_j^2\dddot{R_j}}{kc} \quad \text{if } d_{ij} < kR \nonumber\\ & \quad = \frac{(R_j^2\dot{R_j})^2}{(d_{ij}+kR)^2(d_{ij}-kR)^2} + \frac{R(2\dot{R}_j^3 + 6R_j\dot{R}_j\ddot{R}_j + R_j^2\dddot{R_j})}{kcd_{ij}} \quad \text{if } d_{ij} > kR. \end{align}

For a resolved gas bubble, $S_e = S_c = 0$; $p_b = p_{v_{sat}}$ for an unresolved vapour bubble. For a $k{\textrm {th}}$ unresolved gas bubble, $p_k$ is obtained from the following equations (accounting for thermal damping):

(4.6a,b)\begin{equation} \frac{1}{R^{3\gamma}}\frac{{\rm d}\left(\,p_kR^{3\gamma}\right)}{{\rm d}t} = \frac{3(\gamma-1)}{R^2} k_w \beta (T_k - T_0), \quad T_k = p_k/{\rho R_g}. \end{equation}

4.1. Special cases

4.1.1. Case $\alpha _{un} = 0$

If $\alpha _{un} = 0$, all the unresolved source terms on the right-hand side of (2.6) become 0, and the equation reduces to

(4.7)\begin{equation} \left.\begin{gathered} \frac{\partial{(\rho Y_l)}}{\partial{t}} + \frac{\partial{(\rho Y_l u_j)}}{\partial{x_j}} =-S_e + S_c, \\ \frac{\partial{(\rho u_i)}}{\partial{t}} + \frac{\partial{(\rho u_i u_j)}}{\partial{x_j}} =-\frac{\partial p}{\partial{x_i}} + \frac{\partial \sigma_{ij}}{\partial{x_j}}, \\ \frac{\partial{(\rho e_s)}}{\partial t} + \frac{\partial{(\rho e_s u_j)}}{\partial{x_j}} = \frac{\partial{Q_j}}{\partial{x_j}} -p\frac{\partial{u_j}}{\partial{x_j}} + \sigma_{ij}\frac{\partial{u_i}}{\partial{x_j}},\\ \frac{\partial{(\rho Y_v)}}{\partial{t}} + \frac{\partial{(\rho Y_v u_j)}}{\partial{x_j}} = S_e - S_c, \end{gathered}\right\} \end{equation}

with the mixture pressure $p$ as

(4.8)\begin{equation} p = \rho Y_l K_l T \frac{p}{p + P_c} + \rho Y_v R_v T. \end{equation}

This is the HMM derived by Gnanaskandan & Mahesh (2015).

4.1.2. Case $\alpha _{res} = 0$

If $\alpha _{res} = 0$ and $\rho _{v_{un}}$ is neglected, then $\rho Y_l = \rho _l(1-\alpha _{un})$ and $\rho Y_{un} = \rho _{v_{un}} \alpha _{un} \sim 0$. In addition, if the liquid is assumed to have an isentropic behaviour, (2.6) would then become

(4.9)\begin{equation} \left.\begin{gathered} \frac{\partial{\rho_l}}{\partial{t}} + \frac{\partial{(\rho_l u_j)}}{\partial{x_j}} = \frac{\rho_l}{1-\alpha_{un}}\left(\frac{\partial{\alpha_{un}}}{\partial{t}} + u_j\frac{\partial{\alpha_{un}}}{\partial{x_j}}\right), \\ \frac{\partial{\rho_l u_i}}{\partial{t}} + \frac{\partial{(\rho_l u_i u_j)}}{\partial{x_j}} = \frac{\rho_l u_i}{1-\alpha_{un}}\left(\frac{\partial{\alpha_{un}}}{\partial{t}} + u_j\frac{\partial{\alpha_{un}}}{\partial{x_j}}\right) - \frac{1}{1-\alpha_{un}}\left(\frac{\partial p}{\partial{x_i}} - \frac{\partial \sigma_{ij}}{\partial{x_j}}\right), \\ p = p_0 + c^2(\rho_l - \rho_{l_{0}}). \end{gathered}\right\} \end{equation}

These are the Eulerian equations for liquid in the EL model developed by Fuster & Colonius (2011).

5. Numerical method

The left-hand side of (2.6) is solved using a predictor–corrector approach developed by Gnanaskandan & Mahesh (2015). The predictor step uses a non-dissipative scheme to compute the convective fluxes at the faces of the cell, and the time advancement is computed using a second-order explicit Adams–Bashforth scheme. The corrector step uses a characteristic-based filtering method to add dissipation locally and capture discontinuities in the flow. This method has been validated for a variety of problems and for details on implementation, refer to Gnanaskandan & Mahesh (2015). We initialize the calculations with background vapour/gas $\alpha = 10^{-6}- 10^{-9}$ in liquid regions and liquid fraction of $10^{-6}/{\text zero}$ in vapour/gas bubbles. The temporal terms on the right-hand side of (2.6) are computed using a second-order backward difference scheme. For the Lagrangian model, we use the fourth-order Runge–Kutta method to compute the temporal derivative term of the $kR$-$RP$ equation. Using such higher-order time-stepping scheme ensures that events such as violent bubble collapse are accurately captured.

6. Results

The HMM, a special case of the current multi-scale model, has been used to study problems such as sheet to cloud cavitation transition over cylinders (Brandao et al. 2020) and wedges (Gnanaskandan & Mahesh 2016; Bhatt & Mahesh 2020) where phenomena such as re-entrant jet formation and bubbly shocks were observed. Hence, in this paper, we choose cases that emphasize the other features of the multi-scale model such as accurately estimating the behaviour of single unresolved bubbles, inter-bubble interactions between the micro-bubbles and the mutual interaction between resolved cavities and unresolved bubbles.

Section 6.1 demonstrates the ability of the model to capture the collapse dynamics of both resolved and unresolved vapour bubbles. The model is validated for gas bubbles in § 6.2. Inter-bubble interactions between gas bubbles are considered in § 6.3. Section 6.4 demonstrates the ability to capture the interaction between resolved and unresolved gas bubbles. Finally, we discuss a complex problem in § 6.5 where $1200$ gas bubbles are exposed to a strong acoustic pulse.

6.1. Resolved and unresolved vapour bubbles

The multi-scale model is first validated for its ability to capture the collapse of both, resolved and unresolved vapour bubbles subjected to a high external pressure. The relevant parameters are as follows:

(6.1) \begin{equation} \left.\begin{gathered} \textrm{resolved bubble}\quad R_{0} = 1\ \textrm{mm}; \quad R_{0}/\Delta x = 50; \quad P_{0} = 0.02\ \textrm{atm} \\ \textrm{unresolved bubble}\quad R_{0} = 0.1\ \textrm{mm}; \quad R_{0}/\Delta x = 0.5; \quad P_{0} = 0.02 \ \textrm{atm}, \end{gathered}\right\} \end{equation}

where $R_0$, $\Delta x$ and $P_0$ are the initial bubble size, cell size and initial bubble pressure, respectively; $k$ is the constant parameter that allows $p(kR)$ to capture the effect of local disturbances on the unresolved bubble – $k=5$ here. The external pressure ($P_{ext}$) for both cases is $1$ atm. For the unresolved bubble case, $\alpha _{res} = 0$ and $\alpha _{un} > 0$ in (2.6). Similarly, $\alpha _{res} > 0$ and $\alpha _{un} = 0$ for the resolved bubble case. Figures 3(a) and 3(b) compare the unresolved and resolved vapour bubble size with the RP equation, respectively. We observe good agreement between the multi-scale model and the reference solution for both cases. Also, the collapse time for both cases agree with the reference solution. The $kR$-$RP$ equation is expected to be independent of $k$ for a single bubble. We choose $k = 5$ and $40$ for this purpose, where $k = 5$ represents regions close to the bubble and $k = 40$ represents regions far away. The comparison is shown in figure 3(c), and we observe a good agreement.

Figure 3. Bubble size comparison between the multi-scale model and the RP ODE for the collapse of a (a) resolved vapour bubble and an (b) unresolved vapour bubble. (c) The independence of the multi-scale model from $k$ is shown for the collapse of the unresolved vapour bubble.

6.2. Unresolved gas bubble

We simulate the collapse of an unresolved gas bubble subjected to a high external pressure. A key difference between vapour and gas bubbles is that the gas pressure varies inversely with the bubble volume. Hence the gas bubble rebounds following its initial collapse. For this problem, $R_0 = 0.1$ mm, $R_0/\Delta x = 0.5$, $P_0 = 0.5$ atm, $P_{ext} = 1$ atm and $k=5$. Figure 4(a) shows the bubble radius for three oscillation cycles and the good agreement between the multi-scale model and the reference RP equation solution.

Figure 4. (a) Comparison of the bubble size between the multi-scale model and the RP ODE for an unresolved gas bubble. (b) The instantaneous pressure at three different locations from the bubble centre ($r = 5R, r = 10R$ and $r = 40R$) is compared with the analytical solution ($\square$). (c) The bubble size is compared for $k = 5$, $10$ and $40$ to demonstrate the insensitivity of the multi-scale model to $k$.

Following its initial collapse, the gas bubble rebounds and continues to oscillate. During these oscillations, it generates pressure waves that propagate away from its surface. The simulated values of pressure are compared with the analytical solution at three different locations: $r$ = $5R$ (close to the bubble), $r = 10R$ (intermediate distance from the bubble) and $r = 40R$ (far away from the bubble). Figure 4(b) shows the good agreement between the analytical and multi-scale model solutions at all three locations. Since the pressure wave decays with distance from the bubble, its amplitude is the highest at $r = 5R$ and the lowest at $r = 40R$. This example illustrates the ability of the multi-scale model to capture the two-way interaction between the bubble and the surrounding liquid. We also compare the solutions for different $k$ ($k = 5$, $10$ and $40$) to test the robustness of the model, as done for the unresolved vapour bubble. Figure 4(c) shows the comparison, and we indeed observe a good agreement.

6.3. Bubble–bubble interaction

6.3.1. A pair of bubbles

When a gas bubble is exposed to an acoustic wave, its initial size and the wavelength of the acoustic wave determine the intensity of its oscillations. However, this intensity can be damped or amplified when another bubble oscillates in the vicinity. The ability of the multi-scale model to capture such bubble–bubble interactions is demonstrated in this section. For this purpose, we compare the solutions of the $kR$-$RP$ equation and the analytical equation (3.20). Equation (3.20) has been widely used to understand the nature of primary and secondary Bjerknes forces that arise during the interaction between bubbles and an acoustic wave. Hence, it is chosen here for comparison. For a pair of bubbles, (3.20) becomes

(6.2)\begin{align} R_1\ddot{R_1}\left(1-\frac{\dot{R_1}}{c}\right) + \frac{3}{2}{\dot{R}_1^2}\left(1-\frac{\dot{R_1}}{3c}\right) &= \left(1+\frac{\dot{R_1}}{c}+\frac{R_1}{c}\frac{\partial}{\partial t}\right)\left(\frac{p_{b_1}-p_\infty}{\rho}\right)\nonumber\\ & \quad - \frac{2{\dot{R}_2^2} R_2 + R_2^2 \ddot{R_2}}{d}, \end{align}

where the subscripts $1$ and $2$ refer to the first and second bubbles, respectively; $d$ denotes the distance between the bubbles. The last term in this equation (${1}/{d}(2{\dot {R}_2^2} R_2 + R_2^2 \ddot {R_2})$) accounts for the interaction between bubbles. We consider two cases for validation. The parameters of the problem are

(6.3)\begin{equation} \left.\begin{gathered} \textrm{Acoustic wave:} \ P = P_0 + \Delta P \sin{(2 {\rm \pi}f t)}; \quad \Delta P = 1.32\ \textrm{bar};\quad f = 75\ \textrm{kHz};\quad P_0 = 1\ \textrm{bar}.\\ \textrm{Initial bubble size:} \ R_{1_0} = R_{2_0} = 50\ \mathrm{\mu} {\rm m}.\\ \textrm{Distance between the bubbles ($d$):} \ 4\ \textrm{mm} \ - \textrm{ first case;} \quad 0.4\ \textrm{mm} \ - \textrm{ second case.} \end{gathered}\right\} \end{equation}

The bubbles are closer to each other in the second case; hence the nature of the interaction is expected to be different; $k = 4$ for both cases, and since the bubbles are identical, results are shown for one of the bubbles. Figure 5(a) shows the first case (bubbles far apart), and figure 5(b) shows the second case (bubbles close to each other). Good agreement is obtained between the multi-scale model and the reference solution for both cases. While the first cycle of oscillation is similar for both cases, the subsequent oscillations get damped when the bubbles are closer, indicating a strong interaction between them. It is important to test the sensitivity of the solution to $k$ as both $p(kR)$ and $I$ are strongly dependent on $k$. Figure 5(c) shows the bubble size for the second case for $k = 2, 4$ and $10$ and we observe a good agreement among them. This strengthens the robustness feature of the $kR$-$RP$ equation as the insensitivity of the bubble dynamics for a cluster of bubbles is non-trivial.

Figure 5. The bubble–bubble interaction feature of the multi-scale model is demonstrated for a pair of unresolved gas bubbles. The bubble radius ($R$) is compared between the multi-scale model and the analytical equation (3.20) when they are (a) far apart ($d = 4$ mm) and (b) close to each other ($d = 0.4$ mm). (c) The independence of the bubble dynamics from $k$ is shown by comparing the bubble size for $k = 2$, $4$ and $10$.

We also perform the simulations with and without $I$ for each $k$ to understand its significance; $kR < d$ for $k = 2$ and $4$, and $kR > d$ for $k = 10$. Hence, $I$ could be negligible for the former but significantly important for the latter, as per the discussion in § 3.2. Figures 6(a), 6(b) and 6(c) show the plots for $k = 2, 4$ and $10$, respectively. The solutions indeed agree for $k = 2$ and $4$; however, they differ for $k = 10$ as expected. In addition, we also analyse the ratio of $I$ to $p(kR)$ for each $k$ as shown in figure 6(d). Firstly, $I$ is at least an order of magnitude larger for $k = 10$ than that for $k = 2$. Secondly, at some time instances (between $T = 40$ and $60\ \mathrm {\mu } \textrm {s}$), $I$ is larger than $p(kR)$ and hence cannot be neglected when $kR > d$. Therefore, $p(kR)$ alone can capture these interactions when $kR < d$. This implies that for $kR < d$, a bubble's dynamics is well represented by the pressure at a finite distance from the bubble when this pressure is in turn affected by the incoming waves from neighbouring bubbles.

Figure 6. The simulations are performed for (a) $k = 2$, (b) $k = 4$ and (c) $k = 10$ with and without $I$ for a pair of bubbles with the separation distance, $d = 0.4\ \textrm {mm}$. (d) The ratio of $I$ to $p(kR)$ is compared for $k = 2$ and $10$.

6.3.2. Multiple bubbles

We consider sixteen bubbles with initial size ($R_0 = 50\ \mathrm {\mu } \textrm {m}$) and distance between adjacent bubbles ($d$) being $0.6$ mm (figure 7a). The dimensions of the domain are $L_x = 100$ mm, $L_y = 100$ mm and $L_z = 500$ mm. The acoustic wave is the same as that used in § 6.3. Since the bubbles lie across different planes along the direction of the wave propagation, the bubbles oscillate out of phase. As a result, the interaction dynamics is more complicated than that of the pair of bubbles discussed in § 6.3.1. We investigate the sensitivity of the bubble dynamics to $k$ by considering three values of $k$: $3,6$ and $12$. Figure 7(b), which shows a two-dimensional sketch of the set-up for four bubbles, gives a physical sense of the chosen $k$. We focus on one bubble (black) for which the spherical surfaces represented by the different $k$ are shown. Note that $k=3$ and $6$ surfaces do not include neighbouring bubbles (i.e. $kR < d$) whereas $k=12$ intersects some neighbouring bubbles (i.e. $kR \leqslant d$). Figure 7(c) shows the average size of the bubbles, and we observe a good agreement for all $k$. This underlines the robustness feature of the multi-scale model.

Figure 7. (a) Sketch of the set-up for the multi-bubble case to study the insensitivity of bubble dynamics to $k$. (b) A two-dimensional sketch of the bubble cluster giving a physical sense of the surfaces represented by $k$. (c) Comparison of the average size of the bubble cluster ($\bar {R}$) among $k = 3$, $6$ and $12$.

6.4. Unresolved and resolved bubble pair

Here, we consider a situation where both unresolved and resolved cavities co-exist. A pair of resolved and unresolved gas bubbles are initially subjected to high external pressure, and their subsequent interaction is studied. Figure 8 shows the schematic of the problem, defined as follows:

(6.4) \begin{equation} \left.\begin{gathered} \textrm{resolved bubble}\quad R_{1_0} = 1\ \textrm{mm}; \quad R_{1_0}/\Delta x = 25; \quad P_{1_{0}} = 0.5\ \textrm{atm} \\ \textrm{unresolved bubble}\quad R_{2_0} = 0.1\ \textrm{mm}; \quad R_{2_0}/\Delta x = 1.25; \quad P_{2_{0}} = 0.5\ \textrm{atm} \\ \textrm{domain size}\quad -10\ \textrm{mm} \leqslant x,y \leqslant 10\ \textrm{mm}; \quad -10\ \textrm{mm} \leqslant z \leqslant 15\ \textrm{mm} \end{gathered}\right\}. \end{equation}

Figure 8. Problem set-up for the interaction between a resolved and an unresolved gas bubble.

The unresolved bubble is much smaller than the resolved bubble. Hence, it will have a negligible impact on the resolved bubble, whereas the resolved bubble can significantly affect its smaller counterpart. Figure 9(a) shows the radius of the unresolved bubble with and without the resolved bubble in its vicinity. Since the initial pressure of the resolved bubble is much less than the ambient pressure, shock waves are generated that travel towards its centre and expansion waves propagate outward. These expansion waves cause the unresolved bubble to experience lower pressure than the ambient pressure. As a result, the intensity of its oscillations decreases during the initial phase, as observed in figure 9(a) (before $t = 0.1$ ms). As the resolved bubble starts to collapse, it generates compression waves that propagate outward. This can be observed in figure 9(b), which is a snapshot taken during the collapse of the resolved bubble. As these compression waves hit the unresolved bubble, its oscillation intensity gets dampened further as observed in figure 9(a) (post $t = 0.15$ ms). Eventually, the unresolved bubble comes to rest, whereas the resolved bubble oscillates unperturbed. Capturing such impact of the resolved bubble on the unresolved counterpart would not have been possible if the Eulerian equations governed only the liquid. This result demonstrates the importance of splitting the net vapour quantities into the constituent resolved and unresolved components, which allows for tracking their behaviour independently and accurately.

Figure 9. (a) The radius ($R$) of the unresolved bubble ($b_{un}$) is shown with and without the resolved bubble ($b_{res}$) in its vicinity. (b) Instantaneous snapshot showing the pressure wave generated by the resolved bubble during its collapse. The bubbles are represented by the iso-contour values of volume fraction ($\alpha$), coloured with pressure ($p$): $\alpha = 0.65$ for the resolved bubble and $\alpha = 0.015$ for the unresolved bubble.

6.5. Bubble-cloud interaction with an acoustic wave

We apply the multi-scale model to a problem with a large number of gas bubbles: $1200$ bubbles exposed to an incoming pressure pulse $P = P_0 + \Delta P \sin (2 {\rm \pi}f t)$ where $P_0 = 1$ atm, $\Delta P = 10$ atm and $f = 300$ kHz (Maeda & Colonius 2018). The initial size of all the bubbles is $R_0 = 10\ \mathrm {\mu } \textrm {m}$ and bubble to cell size ratio ($R_0/\Delta x)$ is $1/10$. The bubbles are distributed randomly in a box with the following dimensions: $-5$ mm $\leqslant x,y,z \leqslant 5$ mm.

Figure 10(a–f) shows the interaction between the pressure pulse and the bubbles. The pressure pulse approaches the bubble cloud from the left (figure 10a). The small size of the bubbles implies a small initial volume fraction ($\alpha _0 \sim 0.001$) as seen in figure 10(b). As the pulse impinges on the cloud (figure 10c), the exterior bubbles (close to the left edge of the cloud) undergo acoustic cavitation first. A large value of $\Delta p$ causes a strong collapse and rebound of these bubbles, with the peak volume fraction reaching ten times the initial volume fraction ($\alpha \sim 0.015$) as observed in figure 10(d). The pulse then passes through the cloud, causing the bubbles in the interior to undergo acoustic cavitation (figure 10f). Meanwhile, the acoustic pulse splits into transmitted and reflected components as it interacts with the bubble cloud (figure 10e). As these pressure waves propagate away from the cloud, the oscillations decay with time, and the bubbles will eventually come to rest.

Figure 10. The interaction between the acoustic pulse and a cluster of $1200$ bubbles is shown via the instantaneous pressure ($p$) and volume fraction ($\alpha$) contour plots at three different instances. (a,b) The acoustic pulse travelling towards the bubble cloud, (c,d) occurrence of acoustic cavitation as the pulse impinges on the bubble cloud and (e,f) reflection and transmission of the pulse by the bubble cloud.

Figures 10(d) and 10(f) show that the expansion of the exterior bubbles ($\alpha \sim 0.015$) appear to be more intense than those in the interior ($\alpha \sim 0.007$). To understand this difference, we compare the size ($R$) of a bubble lying at the edge of the cloud ($z = -5L/10$) with the size of two bubbles lying in the interior of the cloud ($z = -3L/10$ and $z = -L/10$) in figure 11(a). The amplitude of the oscillations is relatively smaller for the interior bubbles implying that they are exposed to a weaker pressure pulse as compared with the exterior bubble. This is the classic shielding phenomenon (Wang & Brennen 1999) wherein the bubbles lying at the edge of the cloud dampen the pressure waves impinging the cloud and thus shield the interior bubbles.

Figure 11. (a) The shielding effect is shown for bubbles at three different locations. The outermost bubble is located at $z = -5L/10$, and the other two bubbles are in the interior at $z = -3L/10$ and $z = -L/10$ ($L$ – size of the bubble cloud). (b) The volume fraction ($\alpha$) of the bubble cloud estimated by the incompressible $kR$-$RP$ equation ($- - -$) and the $kR$-$RP$ equation (——–) are compared with the reference solution ($\square$).

The initial strong expansion of the exterior bubble is followed by an intense collapse where its size decreases almost by a factor of $10$ (figure 11a). At such a large collapse velocity, the liquid compressibility effects on the bubble can be significant. To assess this impact, we perform the simulation with both the $kR$-$RP$ equation and its incompressible version (3.22). We compare the volume fraction of the cloud for these two cases with the reference solution (Maeda & Colonius 2018) in figure 11(b). First, we observe good agreement between the multi-scale model and the reference solution. In addition, reasonable agreement is observed among the three curves for the first two oscillation cycles. For the subsequent oscillations, the incompressible version estimates higher peaks for $\alpha$ than the $kR$-$RP$ equation. This is because the latter accounts for medium compressibility that dampens violent oscillations of the bubbles. As a result, the incompressible version estimates the bubbles to oscillate longer, causing a slower decay in $\alpha$ during the later stages. This difference in the solutions suggests the need to account for medium compressibility to simulate such flows. Note that Maeda & Colonius (2018) used an RP equation derived by Fuster & Colonius (2011) where they expressed $p_\infty$ in terms of the local liquid pressure and velocity potential of the $N$ neighbouring bubbles to capture local flow effects. The velocity potential of the bubbles was modelled separately using Bernoulli's equation. Obtaining the velocity potential of these $N$ bubbles required solving a system of $N$ additional equations which can become computationally expensive. In contrast, the derivation of $kR$-$RP$ equation in terms of $p(kR)$ and the near-field approximation of the velocity potential allows for the local flow effect and the inter-bubble interactions to be taken into account without being computationally expensive.

7. Summary

Cavitating flows possess a wide range of length scales of vapour, ranging from massive vapour cavities to micro-bubbles. Both vapour cavities and micro-bubbles can display highly compressible behaviour. Hence, we propose a compressible multi-scale model that captures the dynamics of vapour cavities that can be resolved on the computational grid alongside subgrid micro-bubbles. A key idea behind the model is to split the net vapour mass, momentum and energy in the compressible homogeneous mixture equations into constituent resolved and unresolved components (e.g. $\rho Y_v = \rho Y_{v_{res}} + \rho Y_{v_{un}}$). Use of the mixture equation allows the large-scale vapour dynamics to be captured, in contrast to approaches that treat the liquid as a carrier fluid for subgrid bubbles. The explicit splitting of variables enables independent treatment of the resolved and unresolved vapour components. The resolved variables are tracked in an Eulerian sense using a transport equation. The subgrid variables are obtained from a Lagrangian representation using a novel RP variant developed in this paper, termed the ‘$kR$-$RP$ equation’. The impact of the unresolved micro-bubbles on the resolved phase is explicitly modelled in terms of the bubble radii, velocities and volume fraction.

The basic idea behind the $kR$-$RP$ equation is that far-field pressure is easily defined for a single bubble but is not as obvious for multiple bubbles. The $kR-RP$ model is therefore formally derived in terms of the pressure at a finite distance, $kR$ from the bubble; $p(kR)$ may therefore be either a near-field or far-field pressure. Such choice is attractive in our opinion for numerical simulations. We show that the model exactly recovers the classical RP equations in the limits that $k$ and sound speed $c$ become very large. Also we show that the results are independent of $k$ for a single bubble for all $k$, and for multiple bubbles when $kR < d$. Numerical results are discussed to demonstrate the robustness of the model to the choice of $k$ which can be different for each bubble if necessary. The proposed multi-scale model reduces to HMM in the absence of unresolved bubbles and to a purely EL model in the absence of resolved vapour cavities. The HMM has been validated against experiment for sheet to cloud cavitation transition over wedges (Gnanaskandan & Mahesh 2016; Bhatt & Mahesh 2020) and cylinders (Brandao et al. 2020). Hence, this paper emphasizes the ability of the multi-scale model to predict other problems such as the behaviour of single resolved and unresolved bubbles, inter-bubble interactions and the mutual interaction between resolved cavities and unresolved bubbles.

The model is able to capture bubble collapse and oscillation for both unresolved and resolved vapour and gas bubbles. The behaviour is independent of $k$ for single bubbles. Inter-bubble interactions are considered first for a pair of unresolved bubbles exposed to an acoustic pulse, and then for a cloud of 16 bubbles. The interaction term, $I$ is shown to be negligible when $kR < d$ ($d$ being the inter-bubble separation distance), and comparable to $p(kR)$ when $kR > d$. This suggests that $p(kR)$ alone can capture the inter-bubble interactions for an appropriate $k$ ($kR < d$). We also demonstrate the independence of bubble dynamics from $k$ when $kR \leqslant d$ for bubble clusters. The model's ability to capture both resolved and unresolved vapour dynamics simultaneously is demonstrated by simulating a problem where a resolved bubble interacts with an unresolved bubble. The collapse of the resolved bubble generates compression waves that dampen the oscillations of the unresolved bubble. Finally, the model is applied to a cluster of $1200$ bubbles exposed to an acoustic pulse. The rapid expansion and violent collapse of the bubbles, as well as the bubble–bubble interactions, are captured accurately by the multi-scale model. Furthermore, the shielding effect is observed where the incoming acoustic wave weakens as it passes through the exterior bubbles.