Energy shielding by cavitation bubble clouds in burst wave lithotripsy

Kazuki Maeda

, Adam D. Maxwell

, Tim Colonius

, Wayne Kreider

, and

Michael R. Bailey

Citation:

The Journal of the Acoustical Society of America

144

, 2952 (2018); doi: 10.1121/1.5079641

View online:

https://doi.org/10.1121/1.5079641

View Table of Contents:

http://asa.scitation.org/toc/jas/144/5

Published by the

Acoustical Society of America

Articles you may be interested in

Translational dynamics of individual microbubbles with millisecond scale ultrasound pulses

The Journal of the Acoustical Society of America

144

, 2859 (2018); 10.1121/1.5063353

On heating of tissues by shear waves generated by ultrasound

The Journal of the Acoustical Society of America

144

, 2962 (2018); 10.1121/1.5079646

Ultrasonic sector imaging using plane wave synthetic focusing with a convex array transducer

The Journal of the Acoustical Society of America

144

, 2627 (2018); 10.1121/1.5065391

Spatial filters suppress ripple artifacts in the computation of acoustic fields with the angular spectrum method

The Journal of the Acoustical Society of America

144

, 2947 (2018); 10.1121/1.5079637

Cavitation bubble nucleation induced by shock-bubble interaction in a gelatin gel

Physics of Fluids

, 051904 (2018); 10.1063/1.5026713

A description of echolocation clicks recorded in the presence of True's beaked whale (Mesoplodon mirus)

The Journal of the Acoustical Society of America

144

, 2691 (2018); 10.1121/1.5067379

Energy shielding by cavitation bubble clouds in burst wave

lithotripsy

Kazuki

Maeda,

Adam D.

Maxwell,

Tim

Colonius,

Wayne

Kreider,

and Michael R.

Bailey

Division of Engineering and Applied Science, California Institute of Technology, 1200 East California

Boulevard, Pasadena, California 91125, USA

Center for Industrial and Medical Ultrasound, Applied Physics Laboratory, University of Washington,

1013 Northeast 40th Street, Seattle, Washington 98105, USA

(Received 20 August 2018; revised 30 October 2018; accepted 31 October 2018; published online

26 November 2018)

Combined laboratory experiment and numerical simulation are conducted on bubble clouds nucle-

ated on the surface of a model kidney stone to quantify the energy shielding of the stone caused by

cavitation during burst wave lithotripsy (BWL). In the experiment, the bubble clouds are visualized

and bubble-scattered acoustics are measured. In the simulation, a compressible, multi-component

flow solver is used to capture complex interactions among cavitation bubbles, the stone, and the

burst wave. Quantitative agreement is confirmed between results of the experiment and the simula-

tion. In the simulation, a significant shielding of incident wave energy by the bubble clouds is quan-

tified. The magnitude of shielding can reach up to 90% of the energy of the incoming burst wave

that otherwise would be transmitted into the stone, suggesting a potential loss of efficacy of stone

comminution. There is a strong correlation between the magnitude of the energy shielding and the

amplitude of the bubble-scattered acoustics, independent of the initial size and the void fraction of

the bubble cloud within a range addressed in the simulation. This correlation could provide for real-

time monitoring of cavitation activity in BWL.

V

2018 Acoustical Society of America

https://doi.org/10.1121/1.5079641

[CCC]

Pages: 2952–2961

I. INTRODUCTION

Cavitation bubble clouds nucleated in the human body

during the passage of tensile components of acoustic waves

are of critical importance for the safety and efficacy of the

treatment of lithotripsy (

Coleman

et al.

, 1987

;

Ikeda

et al.

2006

;

Lingeman

et al.

, 2009

;

Miller

et al.

, 2012

;

Pishchalnikov

et al.

, 2003

;

Stride and Coussios, 2010

;

Tanguay, 2004

). Shock wave lithotripsy (SWL) uses a shock

waveform of pulses with a peak amplitude of

(10–100)

MPa. The strong tensile tail results in the formation of a

large bubble cluster that can violently collapse to cause tis-

sue injury and enhance stone comminution, the former of

which has been seen as a major limitation of SWL (

Bailey

et al.

, 2006

;

Evan

et al.

, 2002

;

McAteer

et al.

, 2005

). Burst

wave lithotripsy (BWL) is an alternative to SWL and uses

high-intensity, focused ultrasound in the form of

(10)

cycles at an amplitude of

(1–10) MPa and a frequency of

(100) kHz for stone comminution (

Maxwell

et al.

, 2015

;

Thoma, 2014

). Due to the low peak amplitude of the wave,

resulting cloud cavitation in BWL is expected to involve

smaller bubbles and consequently less violent bubble

dynamics. A recent experimental study

in vitro

has shown

that cavitation bubble clouds with a size of

(1) mm are

formed during a passage of the burst wave (

Maeda

et al.

2015

). A subsequent study has shown that the bubble clouds

can scatter a large portion of incident wave and reduce the

wave energy further downstream (

Maeda and Colonius,

2018a

). Therefore, it has been hypothesized that the bubble

clouds in BWL may cause energy shielding of nearby kidney

stones. The energy shielding could lower the efficacy of

stone comminution, and thus it is desirable to characterize

the process as a function of the incident wave and expected

size and number of bubble nuclei.

To investigate this phenomenon, we study bubble clouds

nucleated on the surface of an epoxy stone model during the

passage of a burst wave by a medical transducer in water

through combined experimental measurements

in vitro

and

companion numerical simulations. In the experiment, we

visualize the evolution of bubble clouds using a high-speed

camera and measure the backscattered acoustics from the

bubbles with the transducer array elements. In the simula-

tion, we combine numerical methods for modeling com-

pressible, multi-component flows to capture complex

interactions among bubbles, the stone, and the burst wave.

Simulated evolution of the bubble cloud and the bubble scat-

tered acoustics quantitatively agree with the results of the

experiment. In the simulation, we vary the initial void frac-

tion and the size of bubble cloud to assess the magnitude of

the shielding as well as its correlation with the bubble-

scattered acoustics. Results indicate that the magnitude of

the shielding reaches up to 90% of the energy of the incident

burst wave that otherwise transmits into the stone. We

Current address: Department of Mechanical Engineering, University of

Washington, Seattle, WA 98195, USA. Electronic mail: kazuki.e.maeda@

gmail.com

Also at: Department of Urology, University of Washington School of

Medicine, 1959 Northeast Pacific Street, Seattle, WA 98195, USA.

2952 J. Acoust. Soc. Am.

144

(5), November 2018

V

2018 Acoustical Society of America

0001-4966/2018/144(5)/2952/10/$30.00

further discover a strong correlation between the magnitude

of the shielding and the amplitude of the backscattered

acoustics, independent of the initial condition of the cloud.

These results can be used to further identify appropriate

parameters and improve the efficacy of stone comminution

in BWL.

II. EXPERIMENTAL SETUP

Figure

1(a)

shows a schematic of the experimental

setup. We generate pulses of ten cycles of a sinusoidal pres-

sure wave with a frequency of 340 kHz and a peak maximum

focal pressure of 7.0 MPa from a multi-element array

focused transducer with an aperture of 180 mm and focal

length of 150 mm [Fig.

1(b)

] toward a cylindrical shape of

epoxy stone model with its axis aligned on the acoustic axis

of the transducer. The transducer is composed of 18 circular

elements made from a ring-shaped piezo-ceramic plate with

an outer diameter of 38.1 mm and inner diameter of

12.7 mm. Each of the elements is designed to generate a

spherical wave front that propagates inward to the center

corresponding to the focal point. Further details of the trans-

ducer are available elsewhere (

Maeda and Colonius, 2017

The pulse repetition frequency (PRF) is 100 Hz. The focal

point of the transducer is located at the center of the stone.

The height and diameter of the stone are 10 and 6.25 mm,

respectively. The water is degassed to approximately 65%

saturation to model the environment of the collecting

space of kidney

in vivo

(

Chaigneau and Le, 1968

;

Hwang

et al.

, 1998

). This condition contrasts with that used in the

initial demonstration of BWL

in vitro

in that water was

highly degassed to suppress cavitation (

Maxwell

et al.

2015

). Preliminary experiments identified that a thin layer of

bubble cloud is formed on the proximal base of the stone,

where the pressure is locally amplified due to reinforcement

of the incoming burst wave by the wave reflected by the

stone. A high-speed camera is triggered to capture an image

of bubble cloud at a specific time during the passage of each

pulse. We concurrently sample the backscattered acoustics

from the bubble cloud and the stone by using the transducer

array elements between the matching network and the

transducer.

III. SIMULATION SETUP

A. Numerical methods

In simulations, we formulate the dynamics of the multi-

component mixture using the compressible, multi-

component, Navier-Stokes equation. We model the stone as

an isotropic, elastic solid with a bulk modulus of

7.14 GPa and a density of

1200 kgm

. The values

of the modulus and the density are chosen as those of the

epoxy resin used for the experimental stone model. The

shear modulus is set to zero. The longitudinal sound speed is

ffiffiffiffiffiffiffiffiffiffiffiffiffiffi

2440 ms

. We neglect the viscosity of the

stone; the acoustic attenuation in the epoxy resin is smaller

than 1 dB/cm and negligible for results considered in the pre-

sent study. This modeling is equivalent to assume that the

stone is an inviscid, compressible fluid with the same values

of the modulus and the density. The coupled dynamics of the

stone and the surrounding water are modeled using an inter-

face capturing method (

Coralic and Colonius, 2014

;

Perigaud and Saurel, 2005

). The interface of the stone is

modeled as a sharp, but continuous region of mixture of

water and the stone. The acoustic impedance mismatch at

the interface is accurately modeled for pressure waves with a

wavelength larger than the grid resolution. The bubbles are

assumed to be spherical and fully immersed in water with a

certain distance from the stone (bubble are not attached on

the stone). For modeling the dynamics of the bubble cloud

excited in an ultrasound field, we use an Eulerian-

Lagrangian method. The method is derived and validated in

detail in

Maeda and Colonius (2018b)

, and applied to para-

metric simulation of the dynamics of spherical bubble clouds

excited by burst waves (

Maeda and Colonius, 2018a

). We

provide a brief summary of the method here. In the method, we

describe the dynamics of bubbly mixture using volume-

averaged equations of motion (

Biesheuvel and vanWijngaarden,

1984

;

Commander and Prosperetti, 1989

;

Fuster and Colonius,

2011

;

Kameda and Matsumoto, 1996

)

þr

ðÞ

;

(1)

ðÞ

þr

ðÞ

;

(2)

FIG. 1. (Color online) (a) Schematic of the experimental setup. (b) Multi-element array medical transducer (same as that modeled in Sec.

J. Acoust. Soc. Am.

144

(5), November 2018

Maeda

etal.

2953

þr

;

(3)

where

is the density,

¼ð

;

v

;

is the velocity,

is the

pressure and

is the total energy.

ðÞ

denotes the volume

averaging operator that acts on arbitrary field variables (

ðÞ ¼ð

ÞðÞ

ðÞ

, where

[0,1) is the volume frac-

tion of gas (void fraction), and subscripts

and

denote the

liquid and gas phase, respectively.

is the effective viscous

stress tensor of the mixture, which we approximate as that of

the liquid phase:

. We invoke two approximations

valid at the limit of low void fraction: the density of the mix-

ture is approximated by that of the liquid:

; the

slip velocity between the two phases is zero:

Equations

(1)

–

(3)

can then be rewritten as conservation

equations in terms of the mass, momentum, and energy of

the liquid with source terms as an inhomogeneous hyper-

bolic system

þr

ðÞ

;

(4)

where

;

(5)

;

(6)

and

ðÞ

(7)

For a thermodynamic closure for the liquid, we employ a

stiffened gas equation of state

¼ð

;

(8)

where

is the internal energy of liquid,

is the specific heat

ratio, and

is the stiffness. In the present study we use

7.1 and

3.06

108 Pa for water.

To model the gas phase, we employ a Lagrangian point-

bubble approach in that the gas is treated as spherical, radi-

ally oscillating cavities consisting of a non-condensible gas

and liquid vapor. The center of

th bubble (

Z

2½

;

), with a radius of

and a radial velocity of

,isini-

tially 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 vol-

ume of bubble using a regularization kernel

Þ¼

;

(9)

where

is the volume of bubble

, and

the distance of the coordinate

from the center of the

bubble,

¼j

. For numerical representation, we dis-

cretize Eq.

(4)

on an axisymmetric grid and spatially inte-

grate using a fifth-order finite volume WENO scheme

(

Coralic and Colonius, 2014

). A fourth/fifth-order Runge-

Kutta-Cash-Karp (RKCK) algorithm (

Cash and Karp, 1990

)

is employed for time integration of solutions.

To model the dynamics of volumetric oscillations of

bubbles, we employ the Keller–Miksis equation (

Keller and

Miksis, 1980

)

€

;

(10)

;

(11)

where

is the pressure at the bubble wall,

is the pres-

sure inside the bubble,

is the surface tension, and

is the

component of the pressure that forces the radial oscillations

of the bubble. We use a reduced-order model to account for

heat and mass transfer across the bubble–liquid interface

(

Preston

et al.

, 2007

). In principle, various similar reduced

order models for heat and mass transfer can be instead used

(

Bergamasco and Fuster, 2017

;

Okita

et al.

, 2013

), while the

choice of the model would not largely affect the results

obtained in this study.

The grid size is uniform with a radial and axial width of

100

m in the stone-wave interaction region. We randomly

distribute spherical bubble nuclei with a radius of 10

min

the cylindrical region with the thickness

and the void frac-

tion

, which faces the proximal base of the stone at the ini-

tial condition, then track the radial evolution of the bubbles

at the sub-grid scale and resolve the bubble-scattered acous-

tics on the grid. The parameters are chosen as discussed

below.

B. Modeling parameters

We model the transducer arrays used in the experiment

as a portion of a spherical surface with an aperture of 60 mm

and a radius of 50 mm with its center located at the origin. In

order to generate focused ultrasound waves from the source

surface, we utilize a source-term approach (

Maeda and

Colonius, 2017

). The method can generate unidirectional

acoustic waves from a surface of arbitrary geometry, in this

case the spherical source surface, by forcing Eq.

(4)

on the

surface in a simulation domain.

For accurate simulation of the bubble dynamics, we cal-

ibrate the acoustic source model using experimental mea-

surements. With the input of

cycles of a sinusoidal

voltage, the output of the transducer is modeled by the fol-

lowing formula:

trans

cos 2

ðÞ

;

(12)

where

and

are the ring-up and ring-down times, respec-

tively. In the simulations of focused waves, we excite this

2954 J. Acoust. Soc. Am.

144

(5), November 2018

Maeda

etal.

expression of the pressure at the source plane, with

1.5

and

2.5

In Fig.

we compare the evolution of the focal pressure

during the passage of the burst wave obtained from the

experiment and a simulation without the stone. The head of

the wave arrives at the focal point at

0. The simulated

waveform and the peak pressure agree well with those of the

measurement, except for the period of small oscillations due

to ring down after

s. In both the experiment and sim-

ulation, the peak negative pressure is

5.5 MPa. The dis-

crepancy between the amplitudes of the peak positive

(7.0 MPa) and the peak negative pressures is due to nonlinear

distortion of the burst wave. In simulation, we capture the

evolution of the bubbles and sample the backscattered acous-

tics at coordinates on the acoustic source surface with polar

angles corresponding to the center of the transducer array

elements used in the experiment.

In the experiment, we do not precisely know the popula-

tion, size distribution, and location of bubbles

apriori

Therefore, we conduct parametric simulations by varying the

thickness and the initial void fraction of the cylindrical bubble

cloud within ranges of

[0.25,1.0] mm and

[1.0,

8.0]

, respectively. For reference, we also simulate a

case without bubble cloud. The bubbles are initially monodis-

perse with radii of

m. For each combination of the

thickness and the void fraction of the bubbly layer, we simu-

late a single realization of the monodisperse bubbles. The

ranges of the parameters were chosen based on the experimen-

tal observation and the results of the previous studies of the

dynamics of a spherical bubble cloud (

Maeda, 2018

;

Maeda

and Colonius, 2018a

). In the studies, we used the same numeri-

cal method to simulate the dynamics of an isolated, spherical

bubble cloud that are nucleated during the passage of a burst

wave with an amplitude of

(10) MPa and a frequency of 335

kHz in water with various values of initial nuclei population

and polydispersity. Through the study, we identified that bub-

ble clouds with an initial mean diameter of 10

mcanmodel

the dynamics of bubble clouds observed in companion experi-

ments well, regardless of the initial polydispersity. For this rea-

son, we have chosen monodisperse bubble clouds with an

initial mean radius of 10

m in the present study.

IV. RESULTS AND DISCUSSION

A. Bubble dynamics

Figure

shows pressure contours at different instants in

time during the simulation with

0.75 mm and

1.0

. A burst wave is generated and focused toward the

FIG. 2. (Color online) Evolution of the focal pressure obtained in the experi-

ment and the simulation without stone.

FIG. 3. (Color online) Snapshots of the pressure contour on the axial-plane during the simulation with

1.0 mm and

1.0

at (a)

14.7, (b)

4.5, (c) 23.7, and (d) 42.9

s. The stone and the void fraction of bubbles are indicated with the black line and by the white shading, respectively. In each image,

the inset shows a zoom-in-view of the stone and the bubble cloud.

J. Acoust. Soc. Am.

144

(5), November 2018

Maeda

etal.

2955

stone [Fig.

3(a)

]. A bubble cloud is excited during the pas-

sage of the wave [Fig.

3(b)

]. A portion of the wave is scat-

tered back from the bubbles and the stone, and the rest of the

wave transmits through or around the stone and is scattered

forward [Fig.

3(c)

]. After the passage of the incident wave,

weak reverberation of the pressure wave is observed in the

stone [Fig.

3(d)

]. Strong collapse of the bubble cloud is not

observed. Figure

shows representative images of a bubble

cloud obtained from the high-speed imaging and the same

simulation at

s. A layer of dispersed bubbles is pre-

sent on the proximal base of the stone in both images.

In Fig.

we compare the evolution of the area occupied

by the bubbles projected on a two-dimensional (2D) plane,

, obtained in the experimental visualization and the sim-

ulation with

0.75 mm and

1.0

. In the experi-

ment, we converted each image into a binary image by using

a post-processing algorithm (Otsu method;

Otsu, 1979

), and

then we count pixels that are occupied by bubbles to com-

pute the area. In the simulation, we output side-view images

from the view angle of the high-speed camera [Fig.

4(b)

We define a 2D Cartesian grid with the same spatial resolu-

tion as the high-speed image, and then map the bubbles onto

the grid (if a bubble overlaps with a certain cell, we regard

the cell as being occupied by the bubble). We then binarize

the image in the same way as we did to the experimental

image. By doing so, the uncertainties are matched between

the experimental image and the snapshots obtained from the

simulation. Each data point of the measurement is obtained

by averaging 440 independent realizations of the bubble

clouds. The uncertainty in the experimental results can be

explained by the stochastic nature of cavitation inception;

the location and the number of nuclei on the surface of the

stone can fluctuate in excitation by distinct pulses.

Nevertheless, we observed statistically consistent behaviors

in the bubble clouds over multiple pulses. The result shows

that the cloud experiences a transient growth of bubbles with

rapid oscillations during the passage of the wave, up to about

s, followed by a smooth decay. The measurement

cannot capture the oscillations, even though they may be

present in the experiments, as the sampling frequency is

lower than that of the oscillation, while the overall growth to

s and subsequent decay after the passage of the wave

agree well with the result of the simulation.

Figure

shows the projected area of bubbles at

obtained in the simulation as a function of the initial thickness

of the bubbly layer. At all values of the initial void fraction,

the projected area increases with the thickness. For most

cases, the area also increases monotonically with void frac-

tion. At all

, initially dilute bubbly layers tend to grow by a

larger factor than dense layers. This can be explained by

mutual suppressions of the volumetric growth among bubbles

in a dense bubble cloud due to inter-bubble interactions

(

Maeda and Colonius, 2018a

;

Wang and Brennen, 1999

). The

data points are spread over the range of

[0.5,2] mm.

To gain further insight into the bubble dynamics, we take

a detailed look at the instantaneous state of bubbles at

s, as a representative time of the maximal bubble cloud

in the experiment and the simulation. Figure

shows the prob-

ability distribution function (PDF) of the projected area of bub-

bles obtained from the experiment at

s. The PDF

approximately follows a log-normal distribution. The range of

the projected area shown in Fig.

covers the wide band of the

experimental PDF shown in Fig.

, within which 92% of the

distribution lies. This suggests that the parameters used in the

simulation represent the experiment well.

B. Scattered acoustics

In order to quantify the amplitude of the backscattered

acoustics, we compute a daylight imaging functional using

FIG. 4. Images of representative bubble clouds obtained at

s from

(a) experiment and (b) simulation. The height and diameter of the stone are

10 and 6.25 mm, respectively.

FIG. 5. (Color online) Evolution of the projected area of bubble cloud.

Results of the experiment and simulation with

0.75 mm and

1.0

are compared. The shaded region corresponds to the stan-

dard deviation of the experimental data points.

FIG. 6. (Color online) The projected area of bubbles at

s as a func-

tion of the cloud thickness obtained in the simulation.

2956 J. Acoust. Soc. Am.

144

(5), November 2018

Maeda

etal.

cross correlations from both the experiment and simulations.

Given

;

, where

1,...,

sensor

, as a set of signals sam-

pled at coordinates

over the time interval [0,

], the cross

correlation between signals

;

and

;

is defined as

;

ðÞ

;

ðÞ

;

ðÞ

(13)

Using the cross correlation, the following imaging functional

at position

can be defined as

Ið

Þ¼

sensors

;

sym

;

Þþ

;

(14)

where

sym

;

Þ¼

;

Þþ

;

(15)

sym

is the symmetric component of the cross correlation.

Further details of the algorithms used to obtain the functional

can be found in

Garnier and Papanicolaou (2009)

Figures

8(a)

and

8(b)

show normalized contours of the

tenth power of imaging functional obtained in a representative

case from the experiment and the simulation with

0.75 mm

and

1.0

, respectively. Positions with a large value

of the imaging functional correspond to estimated locations of

the acoustic source. Note that we use the power to visually

enhance the contrast in the contours (gamma correction;

Poynton, 2012

), while the estimated location of the source is

unchanged. In both plots, a region of a large value of the imag-

ing functional is localized within the stone.

Figure

shows correlations between the area occupied

by the bubbles and the scattering factor at

s. The

scattering factor is defined as

max

ref

;

(16)

where

ref

is the reference value of imaging functional without

bubbles. Therefore,

quantifies the enhancement/decay of

the scattered acoustics due to the presence of bubbles. Weak

positive correlations are observed between the projected area

and the scattering factor in both the experiment and the sim-

ulation. The data points of the simulation show that the scat-

tering is enhanced with the increase in the initial void

fraction and the thickness of the bubble cloud. Most of the

data points obtained from the simulation lie in the interior of

the 92% confidence ellipse obtained from the experimental

data points, indicating a good agreement of the scattered

acoustics measured in the experiment and that obtained in

the simulation.

Notice that in both the experiment and simulation, some

data points are in the range of

1. This indicates that the

amplitude of the backscattering is smaller than that without

bubble clouds in these cases. This mitigation of the ampli-

tude can be explained by the widening of the scattering angle

by the bubble clouds. Without bubbles, the incoming wave

tends to get reflected perpendicular to the flat base of the

stone, and the scattered wave is most intense opposite the

direction of the incident wave. On the other hand, the bubbly

layer results in a more omnidirectional scattering. With

increasing the thickness and the density of the bubbly layer,

thus increaseing

, the intensity of the bubble-scattered

acoustics becomes enhanced and this realizes

The good agreement between the experiment and the

simulation shown in this section and Sec.

IV A

suggest that

the simulation correctly captures the physics of cavitation

near the stone model, and thus can be used to predict quanti-

ties inaccessible in the experiment.

C. Energy shielding

In our previous computational study of the dynamics of

spherical bubble clouds excited by burst waves, it was shown

FIG. 7. Probability distribution function (PDF) of the projected area of bub-

bles obtained in the experiment at

s. Solid line indicates a log-

normal fit to the data.

FIG. 8. (Color online) Normalized contour of the tenth power of imaging

functional near the focal region obtained from (a) a representative case in

the experiment and (b) the simulation with

0.75 mm and

1.0

. Stones are indicated by dotted lines. The length unit is mm.

J. Acoust. Soc. Am.

144

(5), November 2018

Maeda

etal.

2957

that the amplitude of the far-field acoustic waves scattered

by the cloud and the kinetic energy of liquid induced by the

bubble dynamics are strongly correlated (

Maeda and

Colonius, 2018a

). This motivates us in this study to assess

correlation between the kinetic energy of the stone excited

by the burst wave as a metric of the energy shielding and the

imaging functional as a metric of the amplitude of the far-

field, scattered acoustics. Direct measurement of the energy

in the stone is challenging in the experiment, and thus we

compute the quantity by processing the data obtained from

the simulation.

Figure

shows the evolution of the kinetic energy in

the stone,

stone

, induced by the burst wave during three dis-

tinct cases from the simulations conducted in this study:

without bubbles; with the thinnest, initially most dilute cloud

(

0.25 mm and

1.0

); with the thickest, ini-

tially densest cloud (

1.0 mm and

8.0

), where

stone

is defined as the kinetic energy of the volume element

integrated over the stone

stone

v

(17)

In all the cases, the energy steadily grows after arrival of the

wave until around at

s, then oscillates around an

average value during the passage of the wave until around

s. After the passage of the wave, the energy steadily

decays and reaches zero around

s. The energy is

highest for all

when there are no bubbles. The energy level

is slightly decreased with the thin and dilute cloud, while it

is drastically reduced with the thick and dense cloud. This

result indicates that, as might be expected, the energy shield-

ing by bubbles is enhanced by increasing the thickness and/

or the void fraction of the cloud.

To quantify the correlation between the shielding and

the scattered acoustics, in Fig.

we plot the shielding factor

as a function of the scattering factor

. The shielding factor

is defined as

ref

;

(18)

where

is the total work done by the acoustic energy to the

stone during each simulation:

stone

ref

is the refer-

ence value of

obtained in the case without bubbles. Note

that

0 and

1 indicate no shielding and perfect shield-

ing (no energy transmission into the stone), respectively.

Notably, the plot indicates a strong positive correlation

between the shielding factor and the scattering factor over

the global range of the data points, independent of the initial

condition of bubbles. The shielding factor grows with the

increase of the scattering factor and reaches

0.9 at

4. This indicates that the energetic state of the bubble cloud

becomes invariant to the initial void fraction and thickness

of the layer as they increase. The result indicates that up to

90% of the energy of the incident burst wave can be

absorbed/scattered by bubbles that would otherwise be trans-

mitted into the stone.

During the wave propagation, in the experiment, the

elastic potential energy and the kinetic energy are excited in

FIG. 9. (Color online) Correlation between the projected area of bubbles at

s and the scattering factor obtained in both the experiment and the

simulation. The black line draws 92% confidence ellipse for the experimen-

tal data points.

FIG. 10. (Color online) Evolution of the kinetic energy in the stone induced

by the burst wave.

FIG. 11. (Color online) Correlation between the shielding factor and the

scattering factor. Dashed line indicates

0.9.

2958 J. Acoust. Soc. Am.

144

(5), November 2018

Maeda

etal.

the stone. The potential energy is composed of the longitudi-

nal and transverse components. For the expression of the

potential energy, see, for instance,

Landau and Lifshitz

(1986)

. In the present simulation, the transverse component

of potential energy is zero since shear stress is not modeled.

Nevertheless, the time-averaged kinetic energy in the stone

excited by the wave equals that of the total elastic potential

energy due to energy conservation, regardless of the magni-

tude of each component. Therefore, we assume that the

kinetic energy is an appropriate metric of the magnitude of

the wave energy transmitted into the stone. The present anal-

ysis may motivate future work to include modeling of shear

stress and the transverse component of the elastic energy as

well as their contributions to stone communition.

D. Implications for BWL

The results of the present study indicate that a thin bub-

bly layer localized on the surface of a stone can reduce trans-

mission of the large portion of incoming wave energy into

the stone. This energy shielding could potentially result in

loss of efficacy of stone comminution during BWL. This

agrees with the hypothesis made in the previous study on the

dynamics of isolated, spherical bubble clouds (

Maeda and

Colonius, 2018a

). Though the scattering and/or the shielding

effects of bubble clouds excited by shock- and high-intensity

focused ultrasound (HIFU) waves have been reported in

both numerical and experimental studies (

Maxwell

et al.

2011

;

Pishchalnikov

et al.

, 2006

;

Tanguay, 2004

;

Yura

et al.

, 2018

), to our knowledge, quantification of the effects

based on direct comparisons of experiment and correspond-

ing computational fluid dynamics simulation has not been

previously made.

The monotonic correlation between the shielding factor

and the scattering factor (Fig.

) suggests that the scattered

acoustics can be directly used to estimate the magnitude of

the shielding regardless of the bubble dynamics, at least

within the range of parameters of bubble clouds addressed in

the present study. The energy shielding, thus the activities of

cavitation bubbles, could be therefore monitored by ultra-

sound sensing outside the human body during the treatment

of BWL. Such monitoring may enable real-time feedback

control of the therapy to optimize the efficacy of stone com-

minution, for instance, by modulating the wave parameters,

including the waveform and the PRF to suppress cavitation

for mitigation of the shielding (

Maeda, 2018

E. Discussion on future strategies to treat cavitation in

BWL

In SWL, the long tensile tail can initiate intense cloud

cavitation in the focal area. Shortly after the passage of the

wave, violent cloud collapse may cause tissue injury (

Bailey

et al.

, 2005

;

Evan

et al.

, 2002

;

Matlaga

et al.

, 2008

Collapse in the proximity to the stone may also enhance

stone comminution (

Pishchalnikov

et al.

, 2003

). During

growth, non-condensable gas diffuses into the bubbles,

which causes persistence of nuclei long after the passage of

the wave and proliferation of bubble clouds that may cause

undesirable energy shielding during the subsequent pulses

(

Pishchalnikov

et al.

, 2011

;

Tanguay, 2004

). The BWL

waveform, on the other hand, lacks the trailing tensile com-

ponent, and this results in milder growth of cavitation bub-

bles compared to those in SWL (

Maeda

et al.

, 2015

). Bubble

clouds are diffuse and do not tend to undergo violent col-

lapse (

Maeda and Colonius, 2018a

). Strong collapse can lead

to fragmentation of bubbles and thus proliferation of bub-

bles; the absence of it in BWL might enable higher PRF and

thus faster treatment than SWL. On the other hand, results of

the present study suggest that the more diffuse bubble clouds

in BWL still lead to a strong shielding of stones from subse-

quent pulses and it may be desirable to avoid them. Overall,

in BWL, suppressing cloud cavitation using feedback may

be an optimal path to achieve both improved safety (

May

et al.

, 2017

;

Movahed

et al.

, 2017

) and high efficacy of stone

comminution.

In contrast to BWL, there have been proposed prototype

methods of HIFU- and histotripsy-based lithotripsy that

actively makes use of cloud cavitation for stone comminu-

tion (

Duryea

et al.

, 2011

;

Ikeda

et al.

, 2006

;

Matsumoto and

Yoshizawa, 2005

;

Yoshizawa

et al.

, 2009

). Particularly, the

HIFU-based method uses HIFU waves with a waveform that

alternately combines wave packets of high and low frequen-

cies, called cavitation-control (C-C) waveform, whose effec-

tiveness was demonstrated in experiments

in vitro

. The high-

and low-frequency components are designed to excite cloud

cavitation on the surface of a stone and trigger violent col-

lapse of the clouds to erode the stone, respectively. Such a

multi-frequency strategy could also be applied to BWL for

enhancement of stone fragmentation.

We need to recall that, however, cavitation may need

precise control such that efficacy and safety are not sacri-

ficed. To date, understanding of cavitation erosion in litho-

tripsy has remained anecdotal; cavitation bubbles in the

vicinity of or attached on a material’s surface may cause ero-

sion and fracture of the material due to non-spherical col-

lapse (

Crum, 1988

;

Johnsen and Colonius, 2009

;

Philipp and

Lauterborn, 1998

;

Zhong

et al.

, 1993

), while the underlying

physical mechanisms, including fluid–structure interactions

and fracture mechanics, have not been fully elucidated. We

believe that further understanding of the complex physical

mechanisms will be crucial in designing better treatment.

V. CONCLUSION

In this paper, the energy shielding of kidney stones by

bubble clouds nucleated on the surface of the stone during

the passage of a burst wave was quantified through a com-

bined experimental and numerical approach. Simulated evo-

lution of the bubble cloud and the bubble-scattered acoustics

using a compressible, multicomponent flow solver showed

quantitative agreement with the results of high-speed pho-

tography and acoustic measurements. Results of the simula-

tion revealed that the magnitude of the energy shielding by a

thin layer of bubble cloud can reach up to 90% of the total

energy of the burst wave that would otherwise be transmitted

into the stone, indicating a large potential loss of efficacy of

the treatment of BWL due to cloud cavitation. The results

extend the previous studies of bubble cloud dynamics in

J. Acoust. Soc. Am.

144

(5), November 2018

Maeda

etal.

2959

ultrasound fields (

Maeda and Colonius, 2018a

;

Maeda

et al.

2015

). Furthermore, we discovered a strong correlation

between the magnitude of the shielding and the amplitude of

the backscattered acoustics. This correlation could be used,

for example, to monitor the magnitude of the shielding in a

kidney by measuring the scattered acoustics outside the

human body in real time during the treatment of BWL, and

adjust parameters to minimize the shielding during a proce-

dure. Future work includes assessment of the effect of the

waveform on the energy shielding as well as extension of the

combined experimental and numerical setups to model envi-

ronments

in vivo

ACKNOWLEDGMENTS

K.M. would like to acknowledge the Funai Foundation

for Information Technology for the Overseas Scholarship.

This work was supported by the National Institutes of Health

under Grant No. P01-DK043881 and K01-DK104854. The

simulations presented here utilized the Extreme Science and

Engineering Discovery Environment, which is supported by

the National Science Foundation Grant No. CTS120005.

Bailey, M., McAteer, J., Pishchalnikov, Y., Hamilton, M., and Colonius, T.

(

2006

). “Progress in lithotripsy research,”

Acoust. Today

(2), 18–29.

Bailey, M. R., Pishchalnikov, Y. A., Sapozhnikov, O. A., Cleveland, R. O.,

McAteer, J. A., Miller, N. A., Pishchalnikova, I. V., Connors, B. A.,

Crum, L. A., and Evan, A. P. (

2005

). “Cavitation detection during shock-

wave lithotripsy,”

Ultrasound Med. Biol.

(9), 1245–1256.

Bergamasco, L., and Fuster, D. (

2017

). “Oscillation regimes of gas/vapor

bubbles,”

Int. J. Heat Mass Transfer

112

, 72–80.

Biesheuvel, A., and vanWijngaarden, L. (

1984

). “Two-phase flow equations

for a dilute dispersion of gas bubbles in liquid,”

J. Fluid Mech.

148

301–318.

Cash, J., and Karp, A. (

1990

). “A variable order Runge-Kutta method for

initial value problems with rapidly varying right-hand sides,”

ACM Trans.

Math. Software

(3), 201–222.

Chaigneau, M., and Le, G. M. (

1968

). “On the composition of gas dissolved

in human urine,” C. R. Seances Acad. Sci., Ser. D

267

(22), 1893–1895.

Coleman, A., Saunders, J., Crum, L., and Dyson, M. (

1987

). “Acoustic cavi-

tation generated by an extracorporeal shockwave lithotripter,”

Ultrasound

Med. Biol.

(2), 69–76.

Commander, K., and Prosperetti, A. (

1989

). “Linear pressure waves in bub-

bly liquids: Comparison between theory and experiments,”

J. Acoust. Soc.

Am.

(2), 732–746.

Coralic, V., and Colonius, T. (

2014

). “Finite-volume WENO scheme for vis-

cous compressible multicomponent flows,”

J. Comput. Phys.

274

, 95–121.

Crum, L. A. (

1988

). “Cavitation microjets as a contributory mechanism for

renal calculi disintegration in ESWL,”

J. Urol.

140

(6), 1587–1590.

Duryea, A. P., Hall, T. L., Maxwell, A. D., Xu, Z., Cain, C. A., and Roberts,

W. W. (

2011

). “Histotripsy erosion of model urinary calculi,”

J. Endourol.

(2), 341–344.

Evan, A. P., Willis, L. R., McAteer, J. A., Bailey, M. R., Connors, B. A.,

Shao, Y., Lingeman, J. E., Williams, J. C., Fineberg, N. S., and Crum, L.

A. (

2002

). “Kidney damage and renal functional changes are minimized

by waveform control that suppresses cavitation in shock wave lithotripsy,”

J. Urol.

168

(4), 1556–1562.

Fuster, D., and Colonius, T. (

2011

). “Modeling bubble clusters in compress-

ible liquids,”

J. Fluid Mech.

688

, 352–389.

Garnier, J., and Papanicolaou, G. (

2009

). “Passive sensor imaging using

cross correlations of noisy signals in a scattering medium,”

SIAM J. Imag.

Sci.

(2), 396–437.

Hwang, E. Y., Fowlkes, J. B., and Carson, P. L. (

1998

). “Variables control-

ling contrast generation in a urinary bladder model,”

J. Acoust. Soc. Am.

103

(6), 3706–3716.

Ikeda, T., Yoshizawa, S., Masataka, T., Allen, J., Takagi, S., Ohta, N.,

Kitamura, T., and Matsumoto, Y. (

2006

). “Cloud cavitation control for

lithotripsy using high intensity focused ultrasound,”

Ultrasound Med.

Biol.

(9), 1383–1397.

Johnsen, E., and Colonius, T. (

2009

). “Numerical simulations of non-

spherical bubble collapse,”

J. Fluid Mech.

629

, 231–262.

Kameda, M., and Matsumoto, Y. (

1996

). “Shock waves in a liquid contain-

ing small gas bubbles,”

Phys. Fluids

(2), 322–335.

Keller, J., and Miksis, M. (

1980

). “Bubble oscillations of large amplitude,”

J. Acoust. Soc. Am.

(2), 628–633.

Landau, L. D., and Lifshitz, E. M. (

1986

Theory of Elasticity, vol. 7,

Course of Theoretical Physics, 3rd edition (Elsevier, New York).

Lingeman, J. E., McAteer, J. A., Gnessin, E., and Evan, A. P. (

2009

“Shock wave lithotripsy: Advances in technology and technique,”

Nat.

Rev. Urol.

(12), 660–670.

Maeda, K. (

2018

). “Simulation, experiments, and modeling of cloud cavita-

tion with application to burst wave lithotripsy,” Ph.D. thesis, California

Institute of Technology.

Maeda, K., and Colonius, T. (

2017

). “A source term approach for generation

of one-way acoustic waves in the Euler and Navier-Stokes equations,”

Wave Motion

, 36–49.

Maeda, K., and Colonius, T. (

2018a

). “Bubble cloud dynamics in an ultra-

sound field,” arXiv preprint arXiv:1805.00129.

Maeda, K., and Colonius, T. (

2018b

). “Eulerian–Lagrangian method for

simulation of cloud cavitation,”

J. Comput. Phys.

371

, 994–1017.

Maeda, K., Kreider, W., Maxwell, A., Cunitz, B., Colonius, T., and Bailey,

M. (

2015

). “Modeling and experimental analysis of acoustic cavitation

bubbles for burst wave lithotripsy,”

J. Phys. Conf. Ser.

656

(1), 012027.

Matlaga, B. R., McAteer, J. A., Connors, B. A., Handa, R. K., Evan, A. P.,

Williams, J. C., Lingeman, J. E., and Willis, L. R. (

2008

). “Potential for

cavitation-mediated tissue damage in shockwave lithotripsy,”

J. Endourol.

(1), 121–126.

Matsumoto, Y., and Yoshizawa, S. (

2005

). “Behaviour of a bubble cluster in

an ultrasound field,”

Int. J. Numer. Methods Fluids

(6-7), 591–601.

Maxwell, A., Cunitz, B., Kreider, W., Sapozhnikov, O., Hsi, R., Harper, J.,

Bailey, M., and Sorensen, M. (

2015

). “Fragmentation of urinary calculi

in vitro

by burst wave lithotripsy,”

J. Urol.

193

(1), 338–344.

Maxwell, A. D., Wang, T.-Y., Cain, C. A., Fowlkes, J. B., Sapozhnikov, O.

A., Bailey, M. R., and Xu, Z. (

2011

). “Cavitation clouds created by shock

scattering from bubbles during histotripsy,”

J. Acoust. Soc. Am.

130

(4),

1888–1898.

May, P. C., Kreider, W., Maxwell, A. D., Wang, Y.-N., Cunitz, B. W.,

Blomgren, P. M., Johnson, C. D., Park, J. S. H., Bailey, M. R., Lee, D.,

Harper, J. D., and Sorensen, M. D. (

2017

). “Detection and evaluation of

renal injury in burst wave lithotripsy using ultrasound and magnetic reso-

nance imaging,”

J. Endourol.

(8), 786–792.

McAteer, J., Bailey, M., Williams, J., Jr., Cleveland, R., and Evan, A.

(

2005

). “Strategies for improved shock wave lithotripsy,” Minerva Urol.

Nefrol.

(4), 271–287.

Miller, D. L., Smith, N. B., Bailey, M. R., Czarnota, G. J., Hynynen, K.,

Makin, I. R. S., and Bioeffects Committee of the American Institute of

Ultrasound in Medicine. (

2012

). “Overview of therapeutic ultrasound

applications and safety considerations,”

J. Ultrasound Med.

(4),

623–634.

Movahed, P., Kreider, W., Maxwell, A., Dunmire, B., and Freund, J. (

2017

“Ultrasound-induced bubble clusters in tissue-mimicking agarphantoms,”

Ultrasound Med. Biol.

(10), 2318–2328.

Okita, K., Sugiyama, K., Takagi, S., and Matsumto, Y. (

2013

“Microbubble behavior in an ultrasound field for high intensity focused

ultrasound therapy enhancement,”

J. Acoust. Soc. Am.

134

(2),

1576–1585.

Otsu, N. (

1979

). “A threshold selection method from gray-level histograms,”

IEEE Trans. Syst. Man Cybern.

(1), 62–66.

Perigaud, G., and Saurel, R. (

2005

). “A compressible flow model with capil-

lary effects,”

J. Comput. Phys.

209

(1), 139–178.

Philipp, A., and Lauterborn, W. (

1998

). “Cavitation erosion by single laser-

produced bubbles,”

J. Fluid Mech.

361

, 75–116.

Pishchalnikov, Y. A., McAteer, J. A., Williams, J. C., Jr., Pishchalnikova, I.

V., and Vonderhaar, R. J. (

2006

). “Why stones break better at slow shock-

wave rates than at fast rates:

In vitro

study with a research electrohydraulic

lithotripter,”

J. Endourol.

(8), 537–541.

Pishchalnikov, Y., Sapozhnikov, O., Bailey, M., Williams, J., Jr., Cleveland,

R., Colonius, T., Crum, L., Evan, A., and McAteer, J. (

2003

). “Cavitation

bubble cluster activity in the breakage of kidney stones by lithotripter

shockwaves,”

J. Endourol.

(7), 435–446.

Pishchalnikov, Y. A., Williams, J. C., Jr., and McAteer, J. A. (

2011

“Bubble proliferation in the cavitation field of a shock wave lithotripter,”

J. Acoust. Soc. Am.

130

(2), EL87–EL93.

2960 J. Acoust. Soc. Am.

144

(5), November 2018

Maeda

etal.