2.1 Bubble cloud dynamics

Here we briefly summarize the physical model and numerical method for simulation of cloud cavitation utilized in the present study. Further details are available in Maeda & Colonius (Reference Maeda and Colonius2017). The method employs a two-way coupled Eulerian–Lagrangian approach. We describe the dynamics of bubbly mixture using volume-averaged equations of motion (Caflisch et al. Reference Caflisch, Miksis, Papanicolaou and Ting1985; Commander & Prosperetti Reference Commander and Prosperetti1989):

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\overline{\unicode[STIX]{x1D70C}\boldsymbol{u}})=0, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}(\overline{\unicode[STIX]{x1D70C}\boldsymbol{u}})}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\overline{\unicode[STIX]{x1D70C}\boldsymbol{u}}\otimes \overline{\boldsymbol{u}}+p{\mathcal{I}}-{\mathcal{T}})=0, & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\overline{E}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }((\overline{E}+p)\overline{\boldsymbol{u}}-{\mathcal{T}}\boldsymbol{\cdot }\overline{\boldsymbol{u}})=0, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}$ is the density, $\boldsymbol{u}=(u,v,w)^{\text{T}}$ is the velocity, $p$ is the pressure and $E$ is the total energy. Here, $\overline{(\cdot )}$ denotes the volume-averaging operator: $\overline{(\cdot )}=(1-\unicode[STIX]{x1D6FD})(\cdot )_{l}+\unicode[STIX]{x1D6FD}(\cdot )_{g}$ , where $\unicode[STIX]{x1D6FD}\in [0,1)$ is the volume fraction of gas (void fraction), and subscripts $l$ and $g$ denote the liquid and gas phase, respectively. Here, ${\mathcal{T}}$ is the effective viscous stress tensor of the mixture, which we approximate as that of the liquid phase: ${\mathcal{T}}\approx {\mathcal{T}}_{l}$ . We invoke two approximations valid at the limit of low void fraction: the density of the mixture is approximated by that of the liquid: $\overline{\unicode[STIX]{x1D70C}}\approx (1-\unicode[STIX]{x1D6FD})\unicode[STIX]{x1D70C}_{l}$ ; and the slip velocity between the two phases is zero: $\overline{\boldsymbol{u}}\approx \boldsymbol{u}_{l}=\boldsymbol{u}_{g}$ .

Equations (2.1)–(2.3) can be rewritten as conservation equations in terms of the mass, momentum and energy of the liquid with source terms, as an inhomogeneous hyperbolic system:

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{l}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{l}\boldsymbol{u}_{l})=\frac{\unicode[STIX]{x1D70C}_{l}}{1-\unicode[STIX]{x1D6FD}}\left[\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}_{l}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FD}\right], & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D70C}_{l}\boldsymbol{u}_{l})}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{l}\boldsymbol{u}_{l}\otimes \boldsymbol{u}_{l}+p{\mathcal{I}}-{\mathcal{T}}_{l})=\frac{\unicode[STIX]{x1D70C}_{l}\boldsymbol{u}}{1-\unicode[STIX]{x1D6FD}}\left[\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}_{l}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FD}\right]-\frac{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(p{\mathcal{I}}-{\mathcal{T}}_{l})}{1-\unicode[STIX]{x1D6FD}}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}E_{l}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left((E_{l}+p)\boldsymbol{u}_{l}-{\mathcal{T}}_{l}\boldsymbol{\cdot }\boldsymbol{u}_{l}\right)=\frac{E_{l}}{1-\unicode[STIX]{x1D6FD}}\left[\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}_{l}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FD}\right]-\frac{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(p\boldsymbol{u}_{l}-{\mathcal{T}}_{l}\boldsymbol{\cdot }\boldsymbol{u}_{l})}{1-\unicode[STIX]{x1D6FD}}. & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

For a thermodynamic closure for the liquid, we employ stiffened gas equation of state:

(2.7) $$\begin{eqnarray}p=(\unicode[STIX]{x1D6FE}-1)\unicode[STIX]{x1D70C}\unicode[STIX]{x1D700}-\unicode[STIX]{x1D6FE}\unicode[STIX]{x03C0}_{\infty },\end{eqnarray}$$

where $\unicode[STIX]{x1D700}=E_{l}-\unicode[STIX]{x1D70C}_{l}\boldsymbol{u}_{l}^{2}/2$ and is the internal energy of liquid, $\unicode[STIX]{x1D6FE}$ is the specific heat ratio and $\unicode[STIX]{x03C0}_{\infty }$ is the stiffness. In the present study we use $\unicode[STIX]{x1D6FE}=7.1$ and $\unicode[STIX]{x03C0}_{\infty }=3.06\times 10^{8}$  Pa for water. Stiffened gas equation of state well models the thermodynamics of liquid under high pressure and has been widely used in simulation of fully compressible gas–liquid flows that involve shock and strong pressure waves (Menikoff & Plohr Reference Menikoff and Plohr1989; Shyue Reference Shyue1998; Coralic & Colonius Reference Coralic and Colonius2014; Tiwari et al. Reference Tiwari, Pantano and Freund2015; Rasthofer et al. Reference Rasthofer, Wermelinger, Hadijdoukas and Koumoutsakos2017).

At the limit of small change in the density of liquid, the equation of state can be linearized as

(2.8) $$\begin{eqnarray}p=p_{0}+c_{0}^{2}(\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{0}),\end{eqnarray}$$

where

(2.9) $$\begin{eqnarray}c=\sqrt{\unicode[STIX]{x1D6FE}(p+\unicode[STIX]{x03C0}_{\infty })/\unicode[STIX]{x1D70C}}\end{eqnarray}$$

is the speed of sound in liquid and the subscript 0 denotes reference states. With $\unicode[STIX]{x1D70C}_{0}=1000~\text{m}^{3}~\text{kg}^{-1}$ , we recover an ambient speed of sound in water, $c_{0}=1475$  m s $^{-1}$ .

To model the gas phase, we employ a Lagrangian point-bubble approach, in that the gas phase is treated as spherical, radially oscillating cavities consisting of a non-condensible gas and liquid vapour. The centre of the $n$ th bubble ( $n\in \mathbb{Z}:n\in [1,N]$ ), with a radius of $R_{n}$ and a radial velocity of ${\dot{R}}_{n}$ , is initially defined at the coordinate $\boldsymbol{x}_{n}$ and tracked as Lagrangian points during simulations. To define the continuous field of the void fraction in the mixture at coordinate $\boldsymbol{x}$ , we smear the volume of bubble using a regularization kernel $\unicode[STIX]{x1D6FF}$ :

(2.10) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}(\boldsymbol{x})=\mathop{\sum }_{n=1}^{N}V_{n}(R_{n})\unicode[STIX]{x1D6FF}(d_{n},h),\end{eqnarray}$$

where $V_{n}$ is the volume of bubble $n$ , $V_{n}=4\unicode[STIX]{x03C0}/3R_{n}^{3}$ , and $d_{n}$ is the distance of the coordinate $\boldsymbol{x}$ from the centre of the bubble, $d_{n}=|\boldsymbol{x}-\boldsymbol{x}_{n}|$ . We discretize equations (2.4)–(2.6) on an axisymmetric grid. The Lagrangian bubbles are distributed in three-dimensional space. The kernel $\unicode[STIX]{x1D6FF}$ maps the volume of bubbles onto the void fraction field defined on axisymmetric coordinates with a kernel width of $h$ . In the present study, we use $h=1.1\unicode[STIX]{x0394}x$ , where $\unicode[STIX]{x0394}x$ is the characteristic width of the computational cell that encloses the bubble. The averaged inter-bubble distance is larger than the size of bubbles as well as $h$ . The bubbles, while three-dimensional, are forced by the axisymmetric pressure field, which is appropriate to obtain spatially averaged quantities of the bubbles addressed in the present study (Maeda & Colonius Reference Maeda and Colonius2017). Spatial integration is realized using a fifth-order finite-volume weighted essentially non-oscillatory scheme (Coralic & Colonius Reference Coralic and Colonius2014). A fourth-/fifth-order Runge–Kutta–Cash–Karp algorithm (Cash & Karp Reference Cash and Karp1990) is employed for time integration of solutions.

To model the dynamics of volumetric oscillations of bubbles, we employ the Keller–Miksis equation (Keller & Miksis Reference Keller and Miksis1980):

(2.11) $$\begin{eqnarray}\displaystyle & \displaystyle \left(R_{n}\left(1-\frac{{\dot{R}}_{n}}{c}\right)\right)\ddot{R}_{n}+\frac{3}{2}{\dot{R}}^{2}\left(1-\frac{{\dot{R}}_{n}}{3c}\right)=\frac{p_{n}-p_{\infty }}{\unicode[STIX]{x1D70C}}\left(1+\frac{{\dot{R}}_{n}}{c}\right)+\frac{R_{n}{\dot{p}}_{n}}{\unicode[STIX]{x1D70C}c}, & \displaystyle\end{eqnarray}$$

(2.12) $$\begin{eqnarray}\displaystyle & \displaystyle p_{n}=p_{Bn}-\frac{4\unicode[STIX]{x1D707}_{l}{\dot{R}}_{n}}{R_{n}}-\frac{2\unicode[STIX]{x1D70E}_{s}}{R_{n}}, & \displaystyle\end{eqnarray}$$

where $p_{n}$ is the pressure at the bubble wall, $p_{Bn}$ is the pressure inside the bubble, $\unicode[STIX]{x1D70E}_{s}$ is the surface tension and $p_{\infty }$ is the component of the pressure that forces the radial oscillations of the bubble. Pressure $p_{\infty }$ is obtained by a sub-grid-scale model as a function of $p_{n}$ and the pressure of the grid cells surrounding bubble $n$ . To close the equation, we consider conservation equations inside the bubble. To do so, we utilize a reduced-order model introduced by Preston, Colonius & Brennen (Reference Preston, Colonius and Brennen2007) to account for the effect of heat transfer, and mass transfer due to evaporation and condensation of vapour across the bubble–liquid interface. A brief remark on the closure is given in appendix B. The model does not account for fusion, fission or non-spherical deformation of bubbles. In diffuse bubble clouds considered in the present study, the averaged inter-bubble distance is so long that fusion of bubbles may be a rare event, while fission and non-spherical deformation may occur at the moment of violent collapse of bubbles (Brennen Reference Brennen2002). We know of no way to account for such effects save directly simulating all the bubbles, which would be prohibitively expensive for the present cloud.

2.2 Acoustic source

In simulations we excite volumetric oscillations of bubbles using plane and focused pressure waves. In order to generate the waves in the computational domain, we utilize a source-term approach introduced by Maeda & Colonius (Reference Maeda and Colonius2017). The method can generate uni-directional acoustic waves from an arbitrarily chosen source surface, by forcing the mass, momentum and energy equations (2.4)–(2.6) in a thin volume enclosing the surface.

3.1 Cloud interaction parameter

d’Agostino & Brennen (Reference d’Agostino and Brennen1989) (hereafter DB) studied the linear response of monodisperse, spherical bubble clouds subjected to harmonic, long-wavelength pressure excitation. DB deduced that the response of the bubble cloud with a low void fraction is characterized by a non-dimensional parameter,

(3.1) $$\begin{eqnarray}B_{0}=\frac{\unicode[STIX]{x1D6FD}_{0}R_{c}^{2}}{R_{b0}^{2}},\end{eqnarray}$$

termed as the cloud interaction parameter. DB found that, when $B_{0}\ll 1$ , the effect of inter-bubble interaction is weak and each bubble in the cloud behaves like a single, isolated bubble. When $B_{0}\gg 1$ , inter-bubble interactions cause bubbles to oscillate coherently at a lower frequency than an isolated single bubble. Wang & Brennen (Reference Wang and Brennen1999) simulated the dynamics of a spherical bubble cloud with various values of $B_{0}$ in the nonlinear regime.

The cloud interaction parameter can be interpreted in different ways. Substituting $\unicode[STIX]{x1D6FD}=N_{b}R_{b0}^{3}/R_{c}^{3}$ into (3.1), $B_{0}$ can be rewritten as

(3.2) $$\begin{eqnarray}B_{0}=\frac{N_{b}R_{b0}^{3}}{R_{c}^{3}}\frac{R_{c}^{2}}{R_{b0}^{2}}=\frac{N_{b}R_{b0}}{R_{c}},\end{eqnarray}$$

where $N_{b}$ is the number of bubbles in the cloud. We notice that this scaling parameter can be independently derived from the Lagrangian mechanics of spherical bubbles under mutual interactions. The global kinetic energy of potential flow of an incompressible liquid induced by volumetric oscillations of $N_{b}$ spherical bubbles can be expressed using a multipole expansion (Takahira, Akamatsu & Fujikawa Reference Takahira, Akamatsu and Fujikawa1994; Ilinskii, Hamilton & Zabolotskaya Reference Ilinskii, Hamilton and Zabolotskaya2007) as

(3.3) $$\begin{eqnarray}K=2\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}_{l}\left[\mathop{\sum }_{i}^{N_{b}}R_{i}^{3}{\dot{R}}_{i}^{2}+\mathop{\sum }_{i}^{N_{b}}\mathop{\sum }_{j}^{N_{b}}\frac{R_{i}^{2}R_{j}^{2}{\dot{R}}_{i}{\dot{R}}_{j}}{r_{i,j}}+O\left(\frac{R^{7}{\dot{R}}^{2}}{r^{4}}\right)\right],\end{eqnarray}$$

where $r_{i,j}$ is the distance between the centres of bubble $i$ and bubble $j$ . The first term in the brackets represents the kinetic energy induced by direct contributions from each bubble and the second term represents the energy induced by the inter-bubble interactions. When bubbles have an approximately uniform size distribution and experience simultaneous change in pressure, we can assume that each bubble takes the same characteristic radius and velocity, $R$ and ${\dot{R}}$ . The characteristic inter-bubble distance can be scaled as $r\sim R_{c}$ . Then $K$ can be scaled as

(3.4) $$\begin{eqnarray}K\sim \unicode[STIX]{x1D70C}_{l}N_{b}R^{3}{\dot{R}}^{2}\left(1+\frac{N_{b}R}{R_{c}}\right).\end{eqnarray}$$

In the limit of small-amplitude oscillations we have $R\approx R_{b0}$ , and therefore we obtain

(3.5) $$\begin{eqnarray}K\sim \unicode[STIX]{x1D70C}_{l}N_{b}R_{b0}^{3}{\dot{R}}^{2}\left(1+B_{0}\right).\end{eqnarray}$$

We see that the interaction parameter dictates the kinetic energy induced by bubbles. With $B_{0}=0$ the kinetic energy is that of $N_{b}$ isolated bubbles, while with $B_{0}>1$ there is an additional contribution from the inter-bubble interactions. Based on equation (3.3), an extended Rayleigh–Plesset equation for the dynamics of the bubbles can be derived (Chahine Reference Chahine1983; Takahira et al. Reference Takahira, Akamatsu and Fujikawa1994; Doinikov Reference Doinikov2004; Bremond et al. Reference Bremond, Arora, Ohl and Lohse2006; Ilinskii et al. Reference Ilinskii, Hamilton and Zabolotskaya2007; Zeravcic, Lohse & Van Saarloos Reference Zeravcic, Lohse and Van Saarloos2011). In fact, the scaling of kinetic energy in terms of $N_{b}R_{b}/R_{c}$ was mentioned by Ilinskii et al. (Reference Ilinskii, Hamilton and Zabolotskaya2007), but was not associated with the parameter derived by DB.

In what follows we specify the size distribution of nuclei to be a log-normal distribution given by $\text{ln}(R_{b}/R_{b,ref})\sim N(0,\unicode[STIX]{x1D70E}^{2})$ . Therefore we employ the following expression for $B_{0}$ :

(3.6) $$\begin{eqnarray}B_{0}\approx \frac{N_{b}R_{b0,ref}}{R_{c}}.\end{eqnarray}$$

3.2 Moments

In order to quantify the anisotropic structure and associated bubble dynamics, we use the moments of either bubble volume or kinetic energy of the liquid, both measured with respect to the initial centre of the cloud (hereafter denoted as the moment of volume and moment of kinetic energy, respectively). The $n$ th moments of the bubble volume and the bubble-induced kinetic energy of liquid are thus respectively defined as

(3.7a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D707}_{Vcn}=\frac{\displaystyle \mathop{\sum }_{bubble}\frac{4\unicode[STIX]{x03C0}}{3}R_{b}^{3}\left(\frac{x_{b}}{R_{c}}\right)^{n}}{\displaystyle \mathop{\sum }_{bubble}\frac{4\unicode[STIX]{x03C0}}{3}R_{b}^{3}}\quad \text{and}\quad \unicode[STIX]{x1D707}_{Kcn}=\frac{\displaystyle \mathop{\sum }_{bubble}2\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}R_{b}^{3}{\dot{R}}_{b}^{2}\left(\frac{x_{b}}{R_{c}}\right)^{n}}{\displaystyle \mathop{\sum }_{bubble}2\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}R_{b}^{3}{\dot{R}}_{b}^{2}}.\end{eqnarray}$$

We will treat the first moment ( $n=1$ ), unless otherwise noted. The moments are normalized to vary within the range of $[-1,1]$ . In an extreme case, when monodisperse bubbles are distributed in a left hemisphere and oscillate with the same radial velocity, the first moments satisfy $\unicode[STIX]{x1D707}_{Vcn}=\unicode[STIX]{x1D707}_{Kcn}=-0.375$ . Therefore, moments smaller than this value indicate a large bias in the volume or kinetic energy towards the proximal side of the cloud.

4.1 Set-up

In order to investigate the dynamics of bubble clouds, we conduct numerical simulations that mimic the laboratory set-up. Figure 4 shows a schematic of the numerical set-up. The size of the simulation domain is $500~\text{mm}\times 250$  mm, which has been verified to be sufficiently large to effectively mimic free space. For the initial condition, we randomly distribute bubble nuclei in a spherical region with radius 2.5 mm with its centre located at the origin of $x$ – $r$ axisymmetric coordinates. The grid size is uniform near the region of the bubble cloud with a characteristic grid size of 100  $\unicode[STIX]{x03BC}\text{m}$ . Symmetry boundary condition is used on the axis of symmetry. The grid is smoothly stretched towards the other domain boundaries, where characteristic boundary conditions are used to reduce spurious reflections of waves. The transducer used in the experiment is modelled by an acoustic source uniformly distributed on a portion of spherical surface with an aperture of 30 mm and a radius of 50 mm concentric with the bubble cloud. The axis of the spherical surface is aligned with the axis of symmetry of the coordinates.

Figure 4.

Schematic of the numerical set-up. Every fifth point is plotted in the computational mesh.

Figure 5.

Focal pressure evolution of the modelled transducer and that of the experimental measurement. Parameter $t^{\ast }$ denotes non-dimensional time, where $t^{\ast }=tf$ .

The modelled acoustic source is calibrated by comparing the focal pressure evolution with an experimental measurement with a low input voltage, producing the good agreement shown in figure 5. In the simulations of bubble clouds, the peak maximum and negative amplitudes are adjusted to 6.0 and $-4.5$  MPa, respectively. In the simulations, we excite pressure waves at the acoustic source surface with a frequency of 335 kHz, such that peak focal maximum and negative pressures without presence of bubbles are 6.0 and $-4.5$  MPa, respectively. This gives a root-mean-square pressure at the focal point of 1.8 MPa during $t\in [0,50]~\unicode[STIX]{x03BC}\text{s}$ . Details of the calibration are described in appendix A.

Table 1.

List of parameters used in the eight runs for simulations of bubble cloud dynamics in the focused ultrasound wave.

The parameters of bubble clouds used in the simulations are summarized in table 1. It is challenging to measure the population and the initial size distribution of nuclei in the experiment. Therefore, we empirically assess the effects of the nuclei population on the resulting bubble cloud dynamics by varying the value of $B_{0}$ within a range of $B_{0}\in [0.625,5]$ . To assess the effect of polydispersity, for each value of $B_{0}$ we simulate monodisperse and polydisperse clouds. For the polydisperse case, the initial radii of bubbles follow a log-normal distribution given as $\text{ln}(R_{b0}/R_{b,ref})\sim N(0,\unicode[STIX]{x1D70E}^{2})$ (Ando, Colonius & Brennen Reference Ando, Colonius and Brennen2011), where $R_{b,ref}$ is the most probable bubble size, chosen as $R_{b,ref}=10~\unicode[STIX]{x03BC}\text{m}$ . In the monodisperse and polydisperse cases, we use $\unicode[STIX]{x1D70E}=0$ and $0.7$ , respectively. Here, $\unicode[STIX]{x1D70E}=0.7$ models highly polydisperse bubble clouds. This is in order to obtain an upper bound of the variability in the resulting bubble dynamics due to polydispersity. We neglect fission/break-up of bubbles during the simulations. In order to assess the variability of the bubble cloud dynamics due to the initial reference radius of bubbles, we also simulated monodisperse and polydisperse clouds with $R_{b0,ref}=5~\unicode[STIX]{x03BC}\text{m}$ with various values of $B_{0}$ : $B_{0}\in [0.625,5]$ . The results did not show a significant difference from cases with $R_{b0,ref}=10~\unicode[STIX]{x03BC}\text{m}$ , and thus they are omitted in this paper.

4.2 Comparisons with the high-speed image

Figure 6 compares the high-speed image (the seventh image of figure 1) and the image of bubbles obtained in run F7 at $t^{\ast }=13.6$ . A similar anisotropic structure is evident in the simulated cloud; the proximal bubbles are larger than the distal bubbles.

Figure 6.

Images of the bubble cloud obtained in the experiment and simulation (run F7) at $t^{\ast }=13.6$ . The red dashed circle has a radius of 2.5 mm.

Figure 7.

Comparisons of the evolutions of the two-dimensional void fraction of bubbles, defined in (4.1), during the experiment and the simulations with (a) monodisperse (runs F1–F4) and (b) polydisperse (runs F5–F8) clouds.

Figures 7(a) and 7(b) compare the evolutions of the two-dimensional void fraction of bubbles in the experiment and the simulation of initially monodisperse and polydisperse clouds, respectively. The two-dimensional void fraction is obtained as

(4.1) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}_{2D}=\frac{A_{B}}{\unicode[STIX]{x03C0}R_{c}^{2}},\end{eqnarray}$$

where $A_{B}$ is the area occupied by the bubbles in the two-dimensional images. In all cases, the projected area steadily grows and reaches its maximum value within the range of 0.1–0.15 at around $t^{\ast }=10{-}15$ then decays. Overall, trends in the evolution of the void fraction are similar between the monodisperse and polydisperse clouds, though the polydisperse clouds present slightly higher peak values than the monodisperse clouds with the same values of $B_{0}$ . The magnitudes of the slope of the void fraction during the growth and the decay are larger in the experiment than the simulation. The discrepancies could be due to experimental uncertainties, including the size distribution of nuclei in the simulation, non-sphericity of the cloud in the experiment and the finite resolution and/or the noise of the high-speed images. Nevertheless, the results confirm that the simulated bubble clouds quantitatively reproduce the experimental observation with reasonable accuracy.

Figure 8.

Evolution of (a,b) the moment of volume, $\unicode[STIX]{x1D707}_{Vc,1}$ , and (c,d) the kinetic energy, $\unicode[STIX]{x1D707}_{Kc,1}$ , of the simulated clouds. Results of (a,c) monodisperse (runs F1–F4) and (b,d) polydisperse (runs F5–F8) clouds are shown.

For quantification of the anisotropic structure, we compute evolutions of the moments of volume and kinetic energy in each cloud during the course of simulations. Figure 8 shows the result. In all clouds, the moment of volume oscillates around $-0.25$ during the passage of the wave until around $t^{\ast }=17$ then grows back to zero. This suggests that the sizes of the proximal bubbles are larger than those of the distal bubbles for all $t^{\ast }$ and the structural anisotropy is the most significant at around $t^{\ast }=17$ . After the initial transient, the moment of kinetic energy oscillates between $-0.25$ and $-0.5$ for all $t^{\ast }$ . This indicates that the proximal bubbles experience a larger amplitude of pressure excitation and oscillate more actively than the proximal bubbles.

The results above indicate that the bubble dynamics is relatively insensitive to both the population and initial polydispersity of the clouds. Therefore, the anisotropic structure is expected to be observed over a wide range of the nuclei distribution and population.

5.1 Set-up

In the set-up considered in the previous section, bubbles are forced by the pressure wave with a complex waveform generated by a specific transducer. This hinders further generalization of the obtained results, including the anisotropic structure and the bubble-induced kinetic energy, to the bubble cloud dynamics excited in other geometries of pressure fields. For generalization, analysis using a wider range of parameters, but with a simpler geometry of acoustic source, is desirable. To this end, as an idealized problem, we conduct parametric simulations of bubble cloud dynamics excited by plane ultrasound waves of various amplitudes.

Table 2.

List of parameters used in the parametric study. The numbering of 1 to 4 after the symbol v denotes values of nuclei densities, corresponding to distinct values of $B_{0}$ : $B_{0}=[0.625,1.25,2.5,5.0]$ , respectively. For each set of $(A,B_{0})$ , five bubble clouds with distinct initial coordinates of bubbles are simulated.

* A $n$ v[1–4] denotes A $n$ v1, A $n$ v2, A $n$ v3 and A $n$ v4, where $n=1,2,\ldots ,6$ .

The set of parameters addressed in the simulations is summarized in table 2. The radius of clouds and variations of $B_{0}$ follow those of the previous section. It is realistic to assume that the radial distribution of bubbles is polydisperse rather than monodisperse in practical conditions. Thus we assume that the distribution of the initial radius of nuclei follows a log-normal distribution with $R_{b0,ref}=10$ $\unicode[STIX]{x03BC}\text{m}$ and $\unicode[STIX]{x1D70E}=0.7$ , but we expect only small differences from monodisperse clouds in the present cases.

The mesh size follows that of the previous section. We excite 10 cycles of a plane, sinusoidal pressure wave from a source plane located at $x=-20$  mm to the positive $x$ direction, which gives the pressure at the origin, without bubble cloud, of

(5.1) $$\begin{eqnarray}p_{a}=p_{0}[1+H(10-t^{\ast })A\text{sin}(2\unicode[STIX]{x03C0}t^{\ast })],\end{eqnarray}$$

where $H$ is the Heaviside step function. The frequency of the wave is $f=300$  Hz, and thus the wavelength is 4.9 mm and approximately equal to the diameter of the bubble clouds. In order to assess the variability of the bubble cloud dynamics due to spatial distribution of bubbles, with each set of $(A,B_{0})$ we simulate $N_{s}=5$ clouds with distinct, random spatial distributions of nuclei. In the highest amplitude case with $A=10$ , the root-mean-square pressure without bubble during the period of excitation is 0.72 MPa. This pressure is of the same order of magnitude as that of the focused wave considered in the previous section. In what follows, we denote quantities obtained by averaging $N_{s}$ bubble clouds with the same set of $(A,B_{0})$ as those of ensemble-averaged cloud. We denote the ensemble average of arbitrary quantity $f$ as

(5.2) $$\begin{eqnarray}f_{ens}=\frac{1}{N_{s}}\mathop{\sum }_{i=1}^{N_{s}}f_{i},\end{eqnarray}$$

where $f_{i}$ is obtained from the $i$ th realization of the bubble cloud. In the present simulations, $N_{s}=5$ is sufficient to obtain ensemble-averaged quantities.

5.2 Anisotropic structure

Here we analyse the volumetric evolution and the anisotropic structure of the clouds. We begin by looking at the highest amplitude case, $A=10$ , in detail.

Figure 9.

Evolution of the moment of volume of (a) all clouds and (b) ensemble-averaged clouds with $A=10$ and various values of $B_{0}$ from runs A6v[1–4].

Figure 9(a) shows evolutions of the moment of volume of bubble clouds from run A6v during the course of simulation. The moment of volume oscillates between $-0.3$ and 0 for all values of $B_{0}$ after initial transient until $t^{\ast }=10$ . After $t^{\ast }=10$ , the range of moment takes on a wider spread in values. In order to assess variability associated with the random position of bubbles, figure 9(b) shows the same quantities of the ensemble-averaged clouds. The similarity of the moments in the two plots indicates small incoherence among the the dynamics of bubble clouds of distinct realizations. The clouds share the same anisotropic structure regardless of the initial population and spatial distribution of nuclei.

Figure 10 shows images of bubble clouds at $t^{\ast }=5.7$ , obtained from one of the realizations from runs A6v[1–4]. As expected, the anisotropic structure is similar to that of the clouds excited by HIFU.

Figure 10.

Images of the bubble clouds with $A=10$ and various values of $B_{0}$ at $t^{\ast }=5.7$ from runs (a) A6v1, (b) A6v2, (c) A6v3 and (d) A6v4.

Figure 11.

Time-averaged moment of volume of all the ensemble-averaged clouds against the normalized amplitude of the incident wave.

We now consider the effect of varying the excitation amplitude, $A$ . In order to simplify the discussion, we concentrate on the dynamics during the excitation phase where $2<t^{\ast }<9$ , and time-average (denoted by $\overline{(\cdot )}$ ) the corresponding moments. Figure 11 shows the time-averaged moment of volume plotted against the normalized amplitude of the incident wave. Regardless of $B_{0}$ , the moment of volume is small and nearly constant with $A$ up to around $A=1$ . For $A>1$ the moment decreases, indicating larger anisotropy. Thus anisotropy is observed with high-amplitude excitation.

To understand this dependency of the structure on the pressure amplitude, we use the Keller–Miksis equation to examine the nonlinear response of a single, isolated spherical bubble with an initial diameter of 10  $\unicode[STIX]{x03BC}\text{m}$ under periodic far-field pressure excitation with a frequency of 300 kHz. Figure 12(a) shows a bifurcation diagram of the radius of the bubble sampled at every period of forcing pressure with a slowly increasing forcing amplitude within the range addressed in the parametric study. The computed radius monotonically grows with $A$ and experiences a sub-harmonic bifurcation at $A\approx 1.65$ , then transits to a chaotic regime with a growing amplitude of radius. The bifurcation diagram in this range of the excitation amplitude was also reported by Preston et al. (Reference Preston, Colonius and Brennen2007). At $A\approx 2.85$ , the radius returns to a quasi-periodic behaviour, then at $A\approx 4$ it re-transits to a chaotic regime with an amplitude growing with $A$ .

Figure 12(b) shows the time-averaged volume of the same bubble during the period of forcing. The growth of the averaged volume follows a similar trend to that of the radius, but with a larger slope. The volume smoothly grows to $\overline{V}_{b}/V_{b0}\approx 5$ with $A$ then discontinuously grows to $\overline{V}_{b}/V_{b0}\approx 12$ at $A\approx 2.85$ . Then it grows with much faster rate with $A$ , towards $\overline{V}_{b}/V_{b0}\approx 50$ at $A=10$ . The nonlinear growth of the volume with $A>1$ corresponds to cavitation.

Figure 12.

(a) Bifurcation diagram of the bubble radius and (b) averaged volume of a single, spherical bubble under periodic pressure excitation with varying amplitude.

Figure 13.

Evolution of the moment of kinetic energy of (a) all clouds and (b) ensemble-averaged clouds with $A=10$ and various values of $B_{0}$ from runs A6v[1–4].

Figure 13(a) shows the evolution of the moment of kinetic energy of clouds A6v. After the initial transient until $t^{\ast }=10$ , the moments of kinetic energy oscillate between $-0.5$ and 0 around an approximately constant level, after which oscillations with larger amplitudes occur. Figure 13(b) shows the same quantities of the ensemble-averaged clouds. The result is similar to that of figure 13(a), further confirming that the trend of the moment results from the coherent dynamics of the cloud. The plots indicate that the oscillations of proximal bubbles are more energetic than the distal bubbles during the course of excitation, regardless of the initial nuclei population. Since the moment of volume and the moment of kinetic energy reach quasi-stationary states during the 10 cycles of pressure excitation, increasing the number of cycles of the pressure excitation may not largely affect the structure of the clouds. When the number of cycles is as small as 1, however, the bubble dynamics does not reach the stationary state, as shown in $t^{\ast }\in [0,1]$ in figure 13, and the structure may not be observed.

The results of single bubble dynamics in figure 12 and the moment of energy in figure 13 may explain the mechanism of the anisotropic structure. The bubbles nearest the source are exposed to an incoming pressure wave, while the distal bubbles experience smaller amplitudes of pressure fluctuations due to the scattering of the wave by the proximal bubbles. This results in larger amplitude of oscillations of bubbles locally in the region near the proximal surface of the cloud, seen as the bias in the moment of kinetic energy. With a pressure amplitude larger than $A>1$ , the proximal bubbles can grow to much larger radius than the distal bubbles due to local cavitation, which results in the bias in the centre of volume, and becomes visible as the anisotropic structure.

The evolutions of the moments of kinetic energy in figure 8(c,d) and those in figure 13(a,b) are similar during the period of ultrasound excitation; they oscillate around 0.4. This indicates that the clouds excited by the focused and plane waves present similar dynamical features including structural anisotropy, and confirms that the parametric study using the plane wave can represent the dynamics of the bubble cloud excited by the focused wave.

5.3 Dynamic cloud interaction parameter

In order to further quantify the bubble cloud dynamics, we seek to generalize the definition of the cloud interaction parameter introduced by DB. The critical difference in the bubble cloud dynamics in the present study and that considered by DB lies in the wavelength and the amplitude of the pressure excitation. As discussed in § 3, the original interaction parameter can be interpreted as a scaling parameter of the global kinetic energy of liquid induced by small-amplitude oscillations of a bubble cluster under weak pressure excitation with long wavelength. Meanwhile, in the bubble clouds considered in the present study, the wavelength is as small as the size of a cloud. Due to the strong amplitude of the pressure, bubbles experience cavitation growth and their radii can deviate from their initial values. The radius of bubbles can also vary in space.

Figure 14.

Evolution of the mean radius of bubbles in ensemble-averaged clouds with $B_{0}=1.25$ and various values of $A$ from runs A[3–6]v2.

Figure 14 shows the evolutions of the spatial mean of the radius of bubbles in the cloud with $B_{0}=1.25$ , normalized by its initial value, with various values of the excitation amplitude: $A=[10^{-0.5},1.0,10^{0.5},10]$ . For $A>1$ , the mean radius grows rapidly on arrival of the wave, then oscillates around an approximately constant value larger than 1 until $t^{\ast }=10$ , while for $A<1$ the mean radius oscillates around 1. After $t^{\ast }>10$ the radius decays to the initial value in all cases. This indicates that the spatial mean of the bubble radius oscillates around the quasi-stationary equilibrium whose value is unique to the pressure amplitude during the course of excitation.

Motivated by this result, we extend the definition of the cloud interaction parameter as

(5.3) $$\begin{eqnarray}B=\frac{N_{b}\langle \overline{R}_{b}\rangle }{R_{c,L}},\end{eqnarray}$$

where $\overline{R}_{b}$ is the time-averaged radius of the bubble during the pressure excitation. Hereafter we denote this parameter as dynamic cloud interaction parameter. A detailed discussion motivating the specific form of $B$ is given in appendix C.

Figure 15.

Dynamic interaction parameter of all the ensemble-averaged clouds plotted against the amplitude of pressure excitation.

In figure 15, we plot the average value of $B$ (over the ensemble) obtained for all the runs in table 1. For all the values of $B_{0}$ , $B$ monotonically grows and deviates from $B_{0}$ for $A>1$ . This is due to deviation of $\langle \overline{R}_{b}\rangle$ from $R_{b,ref}$ and thus can be associated with the cavitation growth of the mean bubble radius in the cloud shown in figure 15 with $A>1$ .

5.4 Scaling of the moment of kinetic energy

The dynamic interaction parameter is proposed as an appropriate scaling parameter for bubble cloud dynamics excited in the short-wavelength regime. In this section, to examine the extent to which $B$ controls the dynamics, we correlate the moment of kinetic energy against both $B$ and the original cloud interaction parameter, $B_{0}$ , and compare the results.

Figure 16.

Scatter plots of the time-averaged moments of kinetic energy against $B_{0}$ and $B$ . (a,b) All realizations; (c,d) the ensemble-averaged values. Circle and diamond denote the first and third moments, respectively.

Figures 16(a) and 16(b) show scatter plots of the time-averaged, first and third moments of kinetic energy against $B_{0}$ and $B$ , respectively. The first and third moments show negative correlations against both $B_{0}$ and $B$ , while the data points are vertically more spread against $B_{0}$ than against $B$ . The dynamic parameter does a somewhat better job of collapsing the dynamics than the original one. However, once we ensemble-average the data in figure 16(c,d), we see that much of the variation is associated with the randomized positions of the bubbles, and, in general, the dynamic interaction parameter collapses the moment of kinetic energy of the clouds. This confirms that the moments of kinetic energy can been seen as monotonic, decreasing functions of $B$ . Overall, the results indicate that $B$ is a more appropriate parameter to scale the moments than $B_{0}$ . The similarity of figures 16(b) and 16(d) indicates small variability of the spatial bias in the energy due to initial spatial distribution of bubbles.

The moments in figure 16(d) approach zero for small $B$ , which confirms that in the limit of $B=0$ inter-bubble interactions are negligible and the resulting spatial bias in the mean kinetic energy is statistically zero, since bubbles experience the same amplitude of pressure excitation at any location in the cloud. The plots also indicate that as $B$ increases, the slope of the curve monotonically decreases and thus the moment saturates. This indicates that the distribution of energetic bubbles in the cloud becomes more localized in the proximal side of the cloud with increasing pressure amplitude, while the magnitude of energy localization eventually becomes invariant to the amplitude.

Overall, the results of the parametric simulation further elucidate the underlying mechanism of the anisotropic structure. When the inter-bubble interaction becomes dominant, energy localization occurs for the cloud, and this happens regardless of the amplitude of pressure excitation. Meanwhile, the anisotropic structure becomes visible only when the energetic, proximal bubbles cavitate and reach a large radius that is a nonlinear function of the amplitude of pressure excitation. It is notable that, in fact, the moment of volume is not collapsed by the dynamic cloud interaction parameter. Figure 17 shows scatter plots of the time-averaged moment of volume against the dynamic cloud interaction parameter. For the entire range of $B$ , the moment is scattered between $-0.3$ and 0.1 for all realizations and the ensemble-averaged values.

Figure 17.

Scatter plots of the time-averaged moments of volume against $B$ . (a) All realizations; (b) the ensemble-averaged values.

5.5 Scaling of the far-field, bubble-scattered acoustics

Given the successful scaling of the moments of kinetic energy in terms of the dynamic interaction parameter, we are motivated to explore scaling of the far-field, bubble-scattered acoustics that result from the bubble cloud dynamics.

Figure 18.

Evolution of the scattered pressure field at a distance $r=8R_{c}$ from ensemble-averaged clouds from runs (a) A6v4 and (b) A1v1.

Figure 18 shows the evolution of the far-field sound at different angles to the direction of incident radiation. The pressure has been normalized by the amplitude of the incident wave. These are plotted for two cases: figure 18(a) shows the densest cloud excited by the highest amplitude wave (thus obtaining the largest value of $B$ ), while figure 18(b) shows the most dilute cloud with the lowest amplitude of excitation (lowest value of $B$ ). The scattered pressure shows sinusoidal oscillations at a retarded time associated with the incident wave scattered to the sampling location. The amplitude of the scattered pressure from the dense cloud is an order of magnitude larger than that from the dilute one for all $t^{\ast }$ . The small, rapid fluctuations at large $t^{\ast }$ are due to the bubble oscillations after the passage of the incident wave. The small amplitude of these fluctuations indicates the absence of a strong, coherent cloud collapse.

Figure 19.

(a) Contours of the scattered pressure field at $t^{\ast }=9.6$ , normalized by the amplitude of the incident wave, obtained from run A6v4. The length unit is $R_{c}=2.5$  mm. (b) Polar plots of the root-mean-square pressure sampled on spherical surface with $r/R_{c}=8$ . Results from runs A6v4 and A1v1 are compared to the scattering of a single spherical air bubble of the same size as the cloud.

Figure 19(a) shows a contour plot of the bubble-scattered component of the pressure field at $t^{\ast }=9.6$ from the dense cloud. The scattered component is obtained by subtracting the contribution of the incident pressure wave from the total pressure field. The scattered wave propagates radially outward from the bubble cloud. Figure 19(b) shows a polar plot of the scattered waves from both clouds averaged over the period of direct scattering. The linear scattering from a single spherical air cavity with the same radius as the clouds is also shown for reference. With both clouds, scattering is dominant over angles in the forward direction. The amplitude of scattering is larger at all angles from the dense cloud than from the dilute one.

Figure 20.

Scatter plots of the root-mean-square pressure against $B_{0}$ and $B$ . (a,b) All realizations; (c,d) the ensemble-averaged values.

Figure 20 is the analogue to figure 16 but with the root-mean-square pressure plotted versus the original and dynamic cloud interaction parameters. Shown by the different colours are the three scattering angles considered in figure 18. The scattered pressure shows positive correlations with both $B_{0}$ and $B$ . The data points are widely spread against $B_{0}$ , but collapse better with $B$ . As was the case with the kinetic energy moment, ensemble averaging of the clouds removes additional scatter associated with the randomized bubble positions and distribution. Overall, the results confirm that the proposed interaction parameter scales the amplitude of the bubble-scattered acoustics better than the original parameter.

The polar plots shown in figure 19 may help explain the saturation of both the moments of kinetic energy and the amplitude of bubble-scattered acoustics with a large value of the dynamic interaction parameter. Due to the large mismatch in the acoustic impedance across the air–water interface, a cavity can scatter the most portion of the incident wave energy. The dense cloud gives a similar magnitude of scattering as a single large bubble. A subsequent increase in either the excitation amplitude or cloud volume fraction (thus increasing $B$ ) yields no further effect; the scattered acoustics saturate at a level similar to a single bubble of the cloud dimension. The smaller directionality of the scattered acoustics by the cloud than by the air cavity is associated with the spatially random distributions of bubbles. In multiple scattering theory, scatterers with a random, disordered distribution may act as a rough surface and result in randomized angles of scattering of the incoming wave, compared to a smooth surface like that of the air cavity (Ishimaru Reference Ishimaru1978).

Overall, it has been shown that the dynamic interaction parameter scales both the amplitude of the scattered acoustic field as well as the moment of kinetic energy. Furthermore, this indicates direct correlations between the far-field acoustics and the moment. Figures 21(a) and 21(b) show data for all clouds considered at the three observer angles, with and without ensemble averaging, respectively. For applications, the result indicates that the measurement of the far-field, bubble-scattered pressure waves can be used as a surrogate for the magnitude of the energy localization in the bubble cloud as well as a means to estimate the value of $B$ .

Figure 21.

Scatter plots of the moments of kinetic energy against $p_{rms}$ for (a) all clouds and (b) ensemble-averaged clouds.