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 (Reference Prosperetti, Crum and Commander1988) 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 (Reference Preston, Colonius and Brennen2007) 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. (Reference Saito, Takami, Nakamori and Ikohagi2007). 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 Reference Apte, Mahesh and Lundgren2008; Raju et al. Reference Raju, Singh, Hsiao and Chahine2011). 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 (Reference Horne and Mahesh2013) 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. Reference Mettin, Akhatov, Parlitz, Ohl and Lauterborn1997; Doinikov Reference Doinikov2001) 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. Reference Bremond, Arora, Ohl and Lohse2006). 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 (Reference Gnanaskandan and Mahesh2015).

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 (Reference Fuster and Colonius2011).

5. Numerical method

The left-hand side of (2.6) is solved using a predictor–corrector approach developed by Gnanaskandan & Mahesh (Reference Gnanaskandan and Mahesh2015). 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 (Reference Gnanaskandan and Mahesh2015). 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.

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 Reference Gnanaskandan and Mahesh2016; Bhatt & Mahesh Reference Bhatt and Mahesh2020) and cylinders (Brandao et al. Reference Brandao, Bhatt and Mahesh2020). 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.