Bubble cloud dynamics in an ultrasound field

Kazuki Maeda

†

Tim Colonius

Division of Engineering and Applied Science, California Institute of Technology, Pasadena, CA

91125, USA

Abstract

The dynamics of bubble clouds induced by high-intensity focused ultrasound are investigated in a

regime where the cloud size is similar to the ultrasound wavelength. High-speed images show that

the cloud is asymmetrical; the bubbles nearest the source grow to a larger radius than the distal

ones. Similar structures of bubble clouds are observed in numerical simulations that mimic the

laboratory experiment. To elucidate the structure, a parametric study is conducted for plane

ultrasound waves with various amplitudes and diffuse clouds with different initial void fractions.

Based on an analysis of the kinetic energy of liquid induced by bubble oscillations, a new scaling

parameter is introduced to characterize the dynamics. The new parameter generalizes the cloud

interaction parameter originally introduced by

d’Agostino & Brennen (1989)

. The dynamic

interaction parameter controls the energy localization and consequent anisotropy of the cloud.

Moreover, the amplitude of the far-field, bubble-scattered acoustics is likewise correlated with the

proposed parameter. Findings of the present study not only shed light on the physics of cloud

cavitation, but may also be of use to quantification of the effects of cavitation on outcomes of

ultrasound therapies including HIFU-based lithotripsy.

1. Introduction

The dynamics of cavitation bubble clouds excited in an intense ultrasound field are of

critical importance for the safety and efficacy of lithotripsy and high-intensity focused

ultrasound (HIFU). In such therapy, cavitation bubbles can be formed in the human body

during the passage of the tensile part of ultrasound pulses. Bubbles can scatter and absorb

subsequent pulses, and the violent collapse of bubbles can cause cavitation damage

(

Coleman

et al.

1987

;

Pishchalnikov

et al.

2003

;

Matsumoto

et al.

2005

;

McAteer

et al.

2005

;

Ikeda

et al.

2006

;

Bailey

et al.

2006

;

Stride & Coussios 2010

;

Miller

et al.

2012

). Due

to the short time scale and three-dimensional nature of cloud cavitation, precise

measurement of individual bubbles has been challenging. Numerical simulations using

mixture-averaging approaches (

vanWijngaarden 1968

;

Biesheuvel & vanWijngaarden 1984

)

have remained central tools for quantification of the dynamics of bubble clouds.

Early studies of bubble cloud dynamics focused on assessment of cavitation noise and

erosion on materials.

Mørch (1980

1982

) theoretically modeled the inward-propagating

collapse of spherical bubble clusters and quantified the resulting collapse pressure.

Omta

†

Email address for correspondence: kazuki.e.maeda@gmail.com;.

Present address: Department of Mechanical Engineering, University of Washington, Seattle, WA 98195, USA.

HHS Public Access

Author manuscript

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Published in final edited form as:

J Fluid Mech

. 2019 March 10; 862: 1105–1134. doi:10.1017/jfm.2018.968.

Author Manuscript

(1987)

studied acoustic emission from the spherical bubble cloud excited by step change of

the pressure in the surrounding liquid.

d’Agostino & Brennen (1989)

formulated the

linearized dynamics of monodisperse, spherical bubble clouds under weak, long wavelength

pressure excitation and identified that the cloud interaction parameter,

βR

, dictates

the linear dynamics of the cloud, where

is the void fraction,

and

are the initial

radius of the cloud and the bubbles, respectively.

Wang & Brennen (1994

1999

) extended

the study to the nonlinear regime, further characterizing the strong collapse of bubble clouds

accompanied by a shockwave.

Shimada

et al.

(2000)

used a similar approach to assess the

effect of the polydispersity of nuclei on the nonlinear dynamics of spherical bubble clouds.

Later, numerical studies of cavitation have gained interest for medical applications.

Tanguay

& Colonius (2003)

extended the mixture-averaging approach to simulate and characterize

the dynamics of cavitation bubble clouds induced in extra-corporeal shock wave lithotripsy

(ESWL).

Matsumoto & Yoshizawa (2005)

extended the method of

Shimada

et al.

(2000)

quantify amplifications in the pressure due to bubble cloud collapse under excitation by

resonant HIFU waves and discussed applications of the collapse energy to kidney stone

comminution as an alternative method of ESWL.

In experiments,

Reisman

et al.

(1998)

used high-speed imaging to observe cloud cavitation

collapse on a finite-span hydrofoil, and analyzed acoustic signals from the cloud collapse,

and associated the results with the inward propagating shockwave predicted in the

aforementioned studies.

Arora

et al.

(2007)

observed collapse in bubble clouds of controlled

nuclei concentration, and observed an inward-propagating collapse with high nuclei

concentrations.

et al.

(2013)

studied the spatial distribution and the translational motions

of bubble clouds in standing acoustic waves in a time scale longer than the collapse.

The aforementioned theoretical studies focus on bubble clouds in an otherwise

incompressible liquid so that the wavelength of the pressure excitation is much larger than

the size of the cloud. In practical conditions of ultrasound therapies, however, the scale

separation invoked above does not hold. In fact,

Maeda

et al.

(2015)

observed bubble clouds

with a size of

(1) mm

in vitro

during the passage of a strong ultrasound wave with a

wavelength as short as the cloud size.

As a representative example of such bubble clouds, in figure 1 we show the evolution of an

isolated bubble cloud nucleated in a pulse of focused ultrasound consisting of 10

wavelengths with a carrier frequency of

= 335 kHz, thus a wavelength of

= 4.4 mm,

generated by a medical transducer in water. The experiment, extended from a setup

documented in

Maeda

et al.

(2015)

, is designed to characterize the dynamics of bubble

clouds in a recently proposed HIFU-based lithotripsy,

burst wave lithotripsy

(

Maxwell

et al.

2015

). Details of the experimental setup are described in appendix A. The cloud is growing

up through the 10th frame, as the ultrasound wave is propagating through the focal region,

and then decays in the subsequent frames. The bubble cloud occupies an approximately

sphereical volume with

≈

2.5 mm. The size of the bubble cloud is thus at the same scale

as the wavelength of the incident wave. Notably, the cloud possesses an anisotropic structure

in that the proximal bubbles grow to a larger radius than the distal bubbles. These dynamics

are significantly different from bubble clouds in the long wavelength regime.

Maeda and Colonius

Page 2

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Author Manuscript

Advanced interface capturing methods are capable of simulating detailed dynamics of each

bubble in a cloud in a compressible liquid at fine spatial scales, and have been applied to

bubble cloud collapse in a free field by

Rossinelli

et al.

(2013)

;

Rasthofer

et al.

(2017)

and

near a wall by

Tiwari

et al.

(2015)

. Yet, such methods are still computationally intensive and

applications are limited to the dynamics within a short time scale, typically that of a single

cycle of bubble collapse. For more complex problems, modeling assumptions have to be

made to reduce the computational cost. In this spirit, we have recently developed a method

that enables simulation of cloud cavitation in an intense ultrasound field with a fine spatial

resolution without a constraint of the scale separation (

Maeda & Colonius 2018

). The

method solves mixture-averaged equations using an Eulerian-Lagrangian approach (

Kameda

& Matsumoto 1996

;

Fuster & Colonius 2011

;

et al.

2015

2018

). The bubbly-mixture is

discretized on an Eulerian grid, while the individual bubbles are tracked on Lagrangian

points and their radial dynamics are solved using the Keller-Miksis equation. By doing so,

the dynamics of each bubble and the bubble-scattered pressure wave are accurately captured

with a reasonable computational expense.

Our simulation successfully reproduces the aforementioned inward-propagating spherical

cloud collapse in the long wavelength regime. Figure 2 shows the evolution of bubbles with

an initial radius of 10

m that are randomly distributed within a sphere of a radius 2.5 mm

and excited by a single cycle of a planer sinusoidal wave with a frequency of 50 kHz, thus a

wavelength of 30 mm, and an amplitude of 1 MPa. The wavelength is much longer than the

cloud size and this effectively models the scale separation. Figure 3 shows the evolution of

the radii of representative bubbles that are initially located at three distinct radial coordinates

of the cloud: cloud center; mid-point between the center and the cloud periphery; periphery.

The peripheral bubble grows to a larger maximum radius and collapse faster than the other

bubbles, while the inner bubbles are subsequently collapsed during the arrival of the inward

propagating bubbly shockwave. The result qualitatively reproduces the numerical simulation

Wang & Brennen (1994

1999

) as well as the experimental observation of

Arora

et al.

(2007)

Motivated by the high-speed imaging and the preliminary simulation, in this paper we aim to

use numerical simulations to provide a first insight into the bubble clouds excited by HIFU

where the cloud radius is commensurate with the wavelength.

The paper is organized as follows. § 2 provides a summary of the modeling and numerical

methods. In § 3 we introduce metrics to quantify the dynamics of bubble clouds, including

the cloud interaction parameter introduced by

d’Agostino & Brennen (1989)

, and moments

of the volume and the kinetic energy. In § 4 we simulate the dynamics of bubble clouds

excited by a focused ultrasound wave with various polydispersities and populations of nuclei

in a setup that mimics the experimental condition. We quantitatively compare results with

the experimental high-speed images shown in figure 1 and evaluate the anisotropic structure.

To further elucidate the dynamics in more generalized conditions, in § 5 we conduct a

parametric study of bubble clouds excited by a plane ultrasound wave, varying the nuclei

populations and the amplitudes of the wave. In § 5.2 we quantitatively analyze the

anisotropic structure, and in § 5.3 we propose a new scaling parameter to characterize the

dynamics of the clouds by generalizing the cloud interaction parameter of

d’Agostino &

Maeda and Colonius

Page 3

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Author Manuscript

Brennen (1989)

. In § 5.4 we collapse the moments of bubble-induced kinetic energy in terms

of the proposed parameter and identify the mechanisms by which energy is localized in the

proximal side of the cloud. In § 5.5 the amplitude and directionality of the scattered acoustic

field are evaluated and collapsed by the proposed parameter. The energy localization and the

scattered acoustics are directly correlated. In § 6 we discuss implications of the numerical

results to the effects of cloud cavitation on outcomes of HIFU-based lithotripsy. In § 7 we

state conclusions.

2. Model formulation

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 (2017)

. 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.

1985

;

Commander & Prosperetti 1989

∂

+ ∇ ⋅ (

) = 0,

(2.1)

∂(

)

∂

+ ∇ ⋅ (

⊗

ℐ −

풯

) = 0,

(2.2)

∂

+ ∇ ⋅ ((

)

−

풯

⋅

) = 0,

(2.3)

where

is the density,

= (

)

is the velocity,

is the pressure and

is the total energy,

respectively.

( ⋅ )

denotes the volume averaging operator:

( ⋅ )

= (1 −

)( ⋅ )

( ⋅ )

, where

∈

[0,1) is the volume fraction of gas (void fraction), and subscripts

and

denote the liquid

and gas phase, respectively.

풯

is the effective viscous stress tensor of the mixture, that we

approximate as that of the liquid phase:

풯

≈

풯

. We invoke two approximations valid at the

limit of low void fraction: the density of the mixture is approximated by that of the liquid:

≈ (1 −

)

; the slip velocity between the two phases is zero:

≈

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:

∂

+ ∇ ⋅

1 −

∂

⋅ ∇

(2.4)

Maeda and Colonius

Page 4

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Author Manuscript

∂

+ ∇ ⋅

⊗

ℐ −

풯

1 −

∂

⋅ ∇

−

∇ ⋅

ℐ −

풯

1 −

(2.5)

∂

+ ∇ ⋅

−

풯

⋅

1 −

∂

⋅ ∇

−

∇ ⋅

−

풯

⋅

1 −

(2.6)

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

= (

− 1)

ρε

−

γπ

∞

(2.7)

where

−

and is the internal energy of liquid,

is the specific heat ratio, and

∞

is the stiffness, respectively. In the present study we use

= 7.1 and

∞

= 3.06×10

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 1989

;

Shyue 1998

;

Coralic &

Colonius 2014

;

Tiwari

et al.

2015

;

Rasthofer

et al.

2017

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

−

(2.8)

where

∞

(2.9)

is the speed of sound in liquid and the subscript 0 denotes reference states. With

= 1000

−1

, we recover an ambient speed of sound in water,

= 1475 ms

−3

To model the gas phase, we employ a Lagrangian point-bubble approach, in that the gas

phase is treated as spherical, radially oscillating cavities consisted of a non-condensible gas

and liquid vapor. The center of

th bubble (

∈ ℤ:

∈ [1,

]

), with a radius of

and a radial

velocity of

, is initially defined at the coordinate

and tracked as Lagrangian points

during simulations. To define the continuous field of the void fraction in the mixture at

coordinate

, we smear the volume of bubble using a regularization kernel

(

) =

∑

= 1

(2.10)

Maeda and Colonius

Page 5

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Author Manuscript

where

is the volume of bubble

= 4

, and

is the distance of the coordinate

from the center of the bubble,

= |

−

|. We discretize equations (2.4)–(2.6) on an axi-

symmetric grid. The Lagrangian bubbles are distributed in three-dimensional space. The

kernel

maps the volume of bubbles onto the void fraction field defined on axisymmetric

coordinates with a kernel width of

. In the present study, we use

= 1.1

Δx

, where

Δ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

. 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 2017

). Spatial integration is realized using a fifthorder finite volume Weighted

Essentially Non-Oscillatory (WENO) scheme (

Coralic & Colonius 2014

). A 4th/5th order

Runge-Kutta-Cash-Karp (RKCK) algorithm (

Cash & Karp 1990

) is employed for time

integration of solutions.

To model the dynamics of volumetric oscillations of bubbles, we employ the KellerMiksis

equation (

Keller & Miksis 1980

1 −

−

∞

1 +

ρc

(2.11)

−

(2.12)

where

is the pressure at the bubble wall,

is the pressure inside the bubble,

is the

surface tension, and

∞

is the component of the pressure that forces the radial oscillations of

the bubble.

∞

is obtained by a sub-grid scale model as a function of

and the pressure of

the grid cells surrounding bubble

. To close the equation, we consider conservation

equations inside the bubble. To do so, we utilize a reduced-order model introduced by

(

Preston

et al.

2007

) to account for the effect of heat transfer, and mass transfer due to

evaporation and condensation of vapor 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 2002

). 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 (2017)

. The method can generate uni-directional

Maeda and Colonius

Page 6

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Author Manuscript

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. Theory and scaling for the dynamics of bubble clouds

3.1 Cloud interaction parameter

d’Agostino & Brennen (1989)

(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)

termed as the

cloud interaction parameter

. DB found that, when

≪

1, the effect of inter-

bubble interaction is weak and each bubble in the cloud behaves like a single, isolated

bubble. When

≫

1, inter-bubble interactions cause bubbles to oscillate coherently at a

lower frequency than an isolated single bubble.

Wang & Brennen (1999)

simulated the

dynamics of a spherical bubble cloud with various values of

in nonlinear regime.

The cloud interaction parameter can be interpreted in different ways. Substituting

into equation (3.1),

can be rewritten as

(3.2)

where

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

spherical bubbles can be expressed using a multipole

expansion (

Takahira

et al.

1994

;

Ilinskii

et al.

2007

) as

= 2

πρ

∑

(3.3)

where

i,j

is the distance between the centers of bubble

and bubble

. The first term in the

bracket 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

Maeda and Colonius

Page 7

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Author Manuscript

the velocity,

and

. The characteristic inter-bubble distance can be scaled as

r ~ R

. Then

can be scaled as

1 +

(3.4)

In the limit of small amplitude oscillations we have

≈

, and therefore we obtain

1 +

(3.5)

We see that the interaction parameter dictates the kinetic energy induced by bubbles. With

= 0 the kinetic energy is that of

isolated bubbles, while with

1 there is an

additional contribution from the inter-bubble interactions. Based on equation (33), an

extended R-P equation for the dynamics of the bubbles can be derived (

Chahine 1983

;

Takahira

et al.

1994

;

Doinikov 2004

;

Bremond

et al.

2006

;

Ilinskii

et al.

2007

;

Zeravcic

et al.

2011

). In fact, the scaling of kinetic energy in terms of

was mentioned by

Ilinskii

al.

(2007)

, but was not associated with the parameter derived by DB.

In what follows we specify the size distribution of nuclei to be log-normal distribution given

by ln(

b,ref

) ~

). Therefore we employ the following expression for

≈

re f

(3.6)

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 center of the cloud (hereafter denoted as the moment of volume and moment of

kinetic energy, respectively). The

-th moments of the bubble volume and the bubble-

induced kinetic energy of liquid are thus respectively defined as:

Vcn

∑

bubble

∑

bubble

and

Kcn

∑

bubble

πρR

∑

bubble

πρR

(3.7)

We will treat the first moment (

= 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 1st moments

satisfy

V cn

Kcn

= −0.375. Therefore, moments smaller than this value indicate a large

bias in the volume or kinetic energy toward the proximal side of the cloud.

Maeda and Colonius

Page 8

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Author Manuscript

4. Cloud cavitation in a focused ultrasound wave

4.1. Setup

In order to investigate the dynamics of bubble clouds, we conduct numerical simulations that

mimic the laboratory setup. Figure 4 shows the schematic of numerical setup. The size of the

simulation domain is 500×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 center located at the origin of

−

axi-

symmetric coordinates. The grid size is uniform near the region of bubble cloud with a

characteristic grid size of 100

m. Symmetry boundary condition is used on the axis of

symmetry. The grid is smoothly stretched toward the other domain boundaries, where

characteristic boundary conditions are used to reduce spurious reflections of waves. The

transducer used in the experiment is modeled 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.

The modeled 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 the 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 the root-mean-square pressure at the focal point of 1.8 MPa during

∈

[0,50]

s. Details of the calibration are described in appendix A.

The parameters of bubble clouds used in the simulations are summarized on 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

within a range of

∈

[0.625,5]. To assess the effect of polydispersity, for each value of

we simulate

monodisperse and polydisperse clouds. For the polydisperse case, the initial radii of bubbles

follow a log-normal distribution given as ln(

b,ref

) ~

) (

Ando

et al.

2011

), where

b,ref

is the most probable bubble size, chosen as

b,ref

= 10

m. In the monodisperse and

polydisperse cases, we use

= 0 and 0.7, respectively.

= 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

,ref

= 5

m with various values of

∈

[0.625,5]. The results did not show a

significant difference from cases with

,ref

= 10

m, thus they are omitted in this paper.

Maeda and Colonius

Page 9

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Author Manuscript

4.2. Comparisons with the high-speed image

Figure 6 compares the high-speed image (the 7th image of fig 1) and the image of bubbles

obtained in run-F7 at

= 13.6. A similar anisotropic structure is evident in the simulated

cloud; the proximal bubbles are larger than the distal bubbles.

Figure 7 (a) and (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

πR

(4.1)

where

is the area occupied by the bubbles on the two-dimensional images. In all the

cases, the projected area steadily grows and reaches its maximum value within the range of

0.1–0.15 at around

* = 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

. 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.

For quantification of the anisotropic structure, we compute evolutions of the moments of

volume and the kinetic energy in each cloud during the course of simulations. Figure 8

shows the result. In all the clouds, the moment of volume oscillates around −0.25 during the

passage of the wave until around

* = 17 then grows back to zero. This suggests that the size

of the proximal bubbles are larger than the distal bubbles for all

* and the structural

anisotropy is the most significant around at

* = 17. After the initial transient, the moment of

kinetic energy oscillates between −0.25 and −0.5 for all

*. 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 are 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. Parametric simulations using plane ultrasound waves

5.1. Setup

In the setup 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

Maeda and Colonius

Page 10

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Author Manuscript

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.

The set of parameters addressed in the simulations is summarized in table 2. The radius of

clouds and variations of

follow 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

ref

= 10

m and

= 0.7, but we expect only small differences

with monodisperse clouds in the present cases.

The mesh size follows the previous section. We excite 10 cycles of a plane, sinusoidal

pressure wave from a source plane located at

= −20 mm to the positive

direction, that

gives the pressure at the origin, without bubble cloud, of

1 +

10 −

sin

πt

(5.1)

where

is the Heaviside step function. The frequency of the wave is

= 300 Hz, 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

) we simulate

= 5 clouds with distinct, random spatial

distributions of nuclei. In the highest amplitude case with

= 10, the root-mean-square

pressure without bubble during the period of excitation is 0.72 MPa. This pressure is in 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

bubble clouds with the same

set of (

A,B

) as those of

ensemble averaged cloud

. We denote the ensemble average of

arbitrary quantity

ens

∑

= 1

(5.2)

where

is obtained from

-th realization of the bubble cloud. In the present simulations,

= 5 is sufficient to obtain ensemble averaged quantities.

5.2. Anisotropic structure

Here we analyze the volumetric evolution and the anisotropic structure of the clouds. We

begin by looking at the highest amplitude case,

= 10, in detail.

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

after initial transient until

* = 10. After

* = 10 the range of moment takes on a

wider spread in values. In order to assess variability associated with the random position of

Maeda and Colonius

Page 11

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Author Manuscript

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

* = 5.7, obtained from one of the realizations

from runs A6v[1–4]. As expected, the anisotropic structure is similar to the clouds excited

by HIFU.

We now consider the effect of varying the excitation amplitude,

. In order to simplify the

discussion, we concentrate on the dynamics during the excitation phase where 2

< t

and time-average (denoted by

( ⋅ )

) the corresponding moments. Figure 11 shows the time

averaged moment of volume plotted against the normalized amplitude of the incident wave.

Regardless of

, the moment of volume is small and nearly constant with

up to around

= 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

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

and

experiences a sub-harmonic bifurcation at

≈

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.

(2007)

. At

≈

2.85, the radius returns to a

quasi-periodic behavior, then at

≈

4 it re-transits to a chaotic regime with an amplitude

growing with

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

≈ 5

with

then discontinuously

grows to

≈ 12

≈

2.85. Then it grows with much faster rate with

, toward

≈ 50

= 10. The nonlinear growth of the volume with

A >

1 corresponds to

cavitation.

Figure 13 (a) shows the evolution of the moment of kinetic energy of clouds A6v. After the

initial transient until

* = 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 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

Maeda and Colonius

Page 12

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Author Manuscript

number of cycles of the pressure excitation may not largely affect the structure of the clouds.

When the number of cycle is as small as 1, however, the bubble dynamics do not reach the

stationary state, as shown 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 center of

volume, and becomes visible as the anisotropic structure.

The evolutions of the moments of kinetic energy in figure 8 (c) and (d) and those in figure

13 (a) and (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 the structural anisotropy, and confirms that the parametric

study using the plane wave can represent the dynamics of 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 those 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 a small amplitude oscillations of 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 shows the evolutions of the spatial mean of the radius of bubbles in the cloud with

= 1.25, normalized by its initial value, with various values of the excitation amplitude:

= [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

* = 10, while with

A <

1 the mean radius oscillates around 1. After

10 the radius decays to the initial value in

all cases. This indicates that the spatial mean of the bubble radius oscillates around their

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

Maeda and Colonius

Page 13

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Author Manuscript

(5.3)

where

is the time averaged radius of bubble during the pressure excitation. Hereafter we

denote this parameter as

dynamic cloud interaction parameter

. A detailed discussion

motivating the specific form of

is given in appendix C.

In figure 15, we plot the average value of

(over the ensemble) obtained for all the runs in

Table 1. For all the values of

monotonically grows and deviates from

for

A >

This is due to deviation of

from

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 >

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

controls the dynamics, we correlate the moment of kinetic energy

against both

and the original cloud interaction parameter,

, and compare the results.

Figure 16(a) and (b) show scatter plots of the time-averaged, first and third moments of

kinetic energy against

and

, respectively. The 1st and 3rd moments show negative

correlations against both

and

, while the data points are vertically more spread against

than against

. 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

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

. Overall, the results indicate that

is more

appropriate parameter to scale the moments than

. The similarity of figure 16 (b) and (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

, which confirms that in the limit of

= 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

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 the 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, the energy

localization occurs to the cloud, and this happens regardless of the amplitude of pressure

Maeda and Colonius

Page 14

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Author Manuscript

excitation. Meanwhile, the anisotropic structure becomes visible only when the energetic,

proximal bubbles cavitate and reach a large radius 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

, the moment is scattered between −0.3 and 0.1 for both all realizations and 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, bubblescattered

acoustics that result from the bubble cloud dynamics.

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

), while figure 18 (b) shows the most

dilute cloud with the lowest amplitude of excitation (lowest value of

). 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 than the dilute one for all

*. The small, rapid

fluctuations at large

* are due to the bubble oscillations after the passage of the incident

wave. The small amplitude of these fluctuations indicate the absence of a strong, coherent

cloud collapse.

Figure 19 (a) shows an contour plot of the bubble-scattered component of the pressure field

* = 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 the dilute one.

Figure 20 is the analog to figure 16 but with the root-mean-square pressure plotted versus

the original and dynamic cloud interaction parameters. Shown by the different colors are the

3 scattering angles considered in figure 18. The scattered pressure shows positive

correlations with both

and

. The data points are widely spread against

, but collapse

better with

. As was the case with the kinetic energy moment, ensemble averaging of the

clouds remove 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

dynamics interaction parameter. Due to the large mismatch in the acoustic impedance across

Maeda and Colonius

Page 15

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Author Manuscript

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

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, scatters 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 1978

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. Figure 21 (a)

and (b) show data for all clouds considered at the 3 observer angles, with and without

ensemble averaging, resepectively. 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

6. Implications for cavitation in lithotripsy

As the central application and motivation of the present study, it is worth discussing

implications of the present results of numerical experiments to HIFU-based lithotripsy.

In ESWL, bubbles in a cloud experience a nearly identical amplitude of pressure excitation

during the passage of the tensile component of the wave since the tensile tail typically has a

much larger width than the cloud size. The dynamics of bubble cloud consist of spherically

symmetric structures, similar to what was shown in figure 2. The inward-propagating

shockwave causes violent cloud collapse, that results in erosion of surrounding materials.

This injurious effect has been seen as a major disadvantage of ESWL (

McAteer

et al.

2005

;

Bailey

et al.

2006

The use of HIFU for lithotripsy has been proposed as an alternative to ESWL due to their

potentials for safer and more efficient stone comminution (

Ikeda

et al.

2006

;

Yoshizawa

al.

2009

;

Maxwell

et al.

2015

). Cavitation bubble clouds in the HIFU-based lithoripsy, by

contrast to those in ESWL, have a cloud size commensurate with the ultrasound wavelength.

The resulting energy localization of the cloud and scattering of the incoming waves

identified in the previous sections indicate that the bubble clouds with a size at an order of

the incident pressure wavelength can result in strong scattering of the incident wave, with

strong implications for HIFU-based lithotripsy.

Figure 22 compares contours of the maximum pressure on the cross plane over the course of

the simulations without bubble cloud, and with bubble clouds with

= 0.625 and 5.0,

respectively. Without the bubble cloud, the region of high-pressure (

5 MPa) is localized to

the focal region of

∈

[−10,5] mm. With the bubble clouds, the region of high pressure does

not penetrate into the cloud except near the proximal surface. This can be interpreted that the

Maeda and Colonius

Page 16

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Author Manuscript

energetic proximal bubbles scatter the incoming wave to prevent the wave from penetrating

into the cloud and suppress excitation and oscillation of the distal bubbles. There exists

energy shielding

. The small values of the maximum pressure within the cloud shown in

figure 22 (b) and (c) also indicate that strong cloud collapse does not occur during or after

passage of the wave. This agrees with the absence of strong acoustic signals from bubble

clouds after the passage of an incident plane wave confirmed in the numerical experiment,

shown in figure 15. The results are reminiscent of a bubble screen(

Carstensen & Foldy 1947

;

Commander & Prosperetti 1989

) that provides a similar shielding effect. Though the high-

pressure region penetrates into the bubble cloud with

= 0.625 (figure 22 (b)) more than

into the cloud with

= 5.0 (figure 22 (c)), the distribution of the maximum pressure outside

the clouds are similar in the two cases. This similarity indicates that the energy shielding can

be caused by bubble clouds with wide ranges of the nuclei population and the initial void

fraction.

In practical conditions of HIFU-based lithotripsy, cavitation bubble clouds can be nucleated

on the surface of a kidney stone. It can be conjectured that such bubble clouds may have

both positive and negative effects on outcomes of the therapy; they can be less injurious due

to the absence of violent cloud collapse, but they could reduce lower the efficacy of stone

comminution by scattering the incident radiation. Meanwhile, it is apparent that a presence

of kidney stones may complicate the resulting bubble cloud dynamics. For instance, non-

spherical bubble collapse may occur on the surface of a stone to cause erosion (

Tomita &

Shima 1986

;

Johnsen & Colonius 2009

), an effect not considered in the present study. For

future research, simulations of bubble cloud dynamics in the presence of a stone are

desirable.

7. Conclusion

We investigated the dynamics of cavitation bubble clouds excited by strong ultrasound

waves in a regime where the cloud size is similar to the ultrasound wavelength. In a first set

of simulations, we excite bubble clouds by a focused ultrasound wave to mimic the

laboratory setup of HIFU-based lithtoripsy. An anisotropic cloud structure was observed in

both experiments and simulations. The proximal bubbles grow to larger radius than the distal

bubbles. In a second series of simulations, we elucidated the underlying mechanisms leading

to the anisotropy of the observed structure and dynamics. In these simulations, we varied the

amplitude of (plane-wave) excitation, the number density of bubbles, and we considered an

ensemble of five runs for each case with different locations and populations of bubbles.

Based on the kinetic energy of liquid induced by oscillations of a bubble cluster, we

proposed a new scaling parameter, namely a dynamic cloud interaction parameter, that

scales the observed anisotropy and dynamics. The parameter is generalized from the cloud

interaction parameter introduced by

d’Agostino & Brennen (1989)

for linearized bubble

cloud dynamics in the long wavelength regime. We likewise showed that the scattered

acoustic field collapses with the same dynamic interaction parameter, and thus can serve as a

surrogate measure for the extent of energy localization in the cloud. This correlation may be

of use to diagnose

in situ

, via acoustic monitoring, the state of cavitation during ultrasound

therapy as well as in various fluid flows that involve cloud cavitation phenomena.

Maeda and Colonius

Page 17

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Author Manuscript

Acknowledgements

The authors thank Adam Maxwell, Wayne Kreider and Michael Bailey for their support in the companion

experiments. K.M acknowledges the Funai Foundation for Information Technology, for the Overseas Scholarship.

This work was supported by the National Institutes of Health under grant P01-DK043881 and ONR Grant

N00014-17-1-2676. The computations presented here utilized the Extreme Science and Engineering Discovery

Environment, which is supported by the National Science Foundation grant number CTS120005.

Appendix A.: Experimental setup

In this section we describe the experimental setup used to obtain images of bubble cloud

shown in figure 1. Figure 23 (a) shows the schematic of the setup. The setup is designed to

capture the evolution of a single, isolated cavitation bubble clouds excited in a focused,

traveling ultrasound wave. The temperature and pressure are ambient. The water is degassed

by a vacuum pump to realize the oxygen level of 75%. A medical transducer composed of

six piezo-ceramic array elements (figure 23 (b)) is immersed in water. The transducer has an

aperture of 110 × 104 mm, and a focal length of 120 mm. An imaging probe is attached at

the center of the transducer. We excite burst waves at the transducer with a pulse-repetition-

frequency (PRF) of 200 Hz. A high-speed camera captures a rectangular region with a

dimension of 15.3×12.5 mm around the focal point of the transducer. The camera captures

14 consecutive frames with a frame rate of 6

s and an exposure time of 50 ns, with a

resolution of 1200 × 980 pixels. A focused passive cavitation detector (PCD) with

Polyvinylidene difluoride (PVDF) membrane, with ROC of 150 mm and an aperture of

50mm, is positioned confocal to the transducer. We use acoustic signals captured by the

imaging probe and the PCD to map the location of cavitation site to confirm that the bubble

cloud captured by the camera is isolated and no other cloud is present outside the window of

the camera. Note that all the high-speed images presented in this paper are vertically

reflected for consistency with the simulations.

Due to the stochastic nature of cavitation inception, it is challenging to nucleate a perfectly

spherical, isolated, single bubble cloud in the experiment. The image of the bubble cloud

shown in figure 1 is not perfectly circular, which represents the typical uncertainty in the

shape of bubble clouds in the experiment. Nevertheless, even though their shapes are not

perfectly circular/spherical, bubble clouds commonly possess the right-left anisotropy, in

that proximal bubbles grow to larger radius than the distal bubbles.

With the input of

cycles of a sinusoidal voltage, the output of the transducer is modeled

by the following formula:

trans

cos(2

π f t

)

1 −

−

1 −

−

(A1)

where

and

are the ring-up and ring-down time, respectively. In the simulations of

focused waves, we excite this expression of the pressure at the source plane, with

= 4.0

and

= 8.0

s. The validity of acoustic source model is confirmed by comparing the focal

pressure with experimental measurement as shown in figure 5. In the experiment, the focal

Maeda and Colonius

Page 18

J Fluid Mech

. Author manuscript; available in PMC 2019 September 26.

Author Manuscript