Cavitation-induced damage of soft materials by focused ultrasound bursts

13 Pages·2016·1.45 MB·English

by  Pooya Movahed

Cavitation-induced damage of soft materials by focused
ultrasound bursts: A fracture-based bubble dynamics model
PooyaMovaheda)
DepartmentofMechanicalScienceandEngineering,UniversityofIllinoisatUrbana–Champaign,
1206WestGreenStreet,Urbana,Illinois61801,USA
WayneKreiderandAdamD.Maxwellb)
CenterforIndustrialandMedicalUltrasound,AppliedPhysicsLaboratory,UniversityofWashington,
1013Northeast40thStreet,Seattle,Washington98105,USA
ShelbyB.HutchensandJonathanB.Freundc)
DepartmentofMechanicalScienceandEngineering,UniversityofIllinoisatUrbana–Champaign,
1206WestGreenStreet,Urbana,Illinois61801,USA
(Received 20 May 2016; revised 27 July 2016; accepted 5 August 2016; published online 30
August2016)
A generalized Rayleigh–Plesset-type bubble dynamics model with a damage mechanism is devel-
opedforcavitationanddamageofsoftmaterialsbyfocusedultrasoundbursts.Thisstudyislinked
to recent experimental observations in tissue-mimicking polyacrylamide and agar gel phantoms
subjectedtoburstsofakindbeingconsideredspecificallyforlithotripsy.Theseshowbubbleactiva-
tionatmultiplesitesduringtheinitialpulses.Morecavitiesappearcontinuouslythroughthecourse
oftheobservations,similartowhatisdeducedinpigkidneytissuesinshock-wavelithotripsy.Two
different materialmodelsareusedtorepresentthedistinctproperties ofthetwogelmaterials. The
polyacrylamide gel is represented with a neo-Hookean elastic model and damaged based upon a
maximum-straincriterion;theagargelisrepresentedwithastrain-hardeningFungmodelanddam-
aged according to the strain-energy-based Griffith’s fracture criterion. Estimates based upon inde-
pendently determined elasticity and viscosity of the two gel materials suggest that bubble
confinementshouldbesufficienttopreventdamageinthegels,andpresumablyinjuryinsometis-
sues.Damageaccumulationisthereforeproposedtooccurviaamaterialfatigue,whichisshownto
beconsistentwithobserveddelaysinwidespreadcavitationactivity.
VC 2016AcousticalSocietyofAmerica.[http://dx.doi.org/10.1121/1.4961364]
[MFH] Pages:1374–1386
I. INTRODUCTION observations. BWL is an attractive alternative to SWL
because it yields more uniform stone fragmentation in labo-
A. Bubblesintherapeuticultrasound
ratory models, with fragment sizes related to the fundamen-
Eitherasanobjectiveorasanundesirableconsequence, tal ultrasound frequency.8 Its lower pressure amplitudes and
application of therapeutic ultrasound often entails cavitation straightforward integration with real-time ultrasound imag-
and subsequent bubble dynamics within soft tissue. Such ing are also hoped to reduce kidney injury9 even while
bubble action is thought to cause injury1–3 in nominally potentially faster clinical treatments are delivered. Still,
noninvasive lithotripsy for comminution of kidney stones, experimental evidence in both animals3,10 and tissue-
both with commonly used shock waves in current clinical mimicking phantom gels9 indicate that cavitation-induced
applications [shock-wave lithotripsy (SWL)] or with pro- damagecanstillpresentachallengeforBWL.
posed high-intensity focused ultrasound bursts [burst-wave While cavitation is an undesirable collateral effect in
lithotripsy (BWL)]. Studies of pig kidneys exposed to SWL BWL, which provides our specific motivation, it is an essen-
pulses—compressive shocks followed by rarefactions—sug- tialfeatureinhistotripsy,inwhichtheviolentactionofcavita-
gest that the damage can be extensive and seems to occur tion bubbles provides the primary mechanism for removing
duringthetreatment,4,5inconcertwiththeacousticsignature tissue,oftentumors.11–14Inthiscase,bettercontrolofcavita-
of cavitation.6,7 We consider SWL waveforms, though the tion might, for example, allow treatments at higher pulse-
present analysis is motivated primarily by BWL, for which repetition frequency (PRF), i.e., faster clinical treatment, yet
we also report some specific recent experimental also suppress tissue injury outside of the target. We investi-
gate possible mechanisms for injury to explain observations
andhelpguideprotocolsandprocedurestoavoidinjury.
a)Electronicmail:[email protected]
b)Also at: Department of Urology, University of Washington School of
Medicine,1959NEPacificStreet,Seattle,WA98195. B. Bubbledynamics
c)Also at: Department of Aerospace Engineering, University of Illinois at
The basic theory for bubble mechanics in liquids is
Urbana–Champaign, 1206 West Green Street, Urbana, Illinois 61801,
USA. foundedontheRayleigh–Plessetequation,whichdescribesthe
1374 J.Acoust.Soc.Am.140(2),August2016 0001-4966/2016/140(2)/1374/13/$30.00 VC 2016AcousticalSocietyofAmerica
FIG.1.Bubblesinpolyacrylamidegel.Theimagedregionis5.41cm(cid:3)2.48cm.TheBWLpulsehas20cycleswithbasefrequency335kHz,peakamplitudeof
7.51MPaandrepeatsat200Hz.Imagesareshownattimes10cyclesintothelabeledpulsenumber.
dynamics of a bubble with finite surface tension in a viscous sufficiently parameterized, but experimental data are cur-
incompressibleNewtonianfluidandsubjectedtofar-fieldpres- rently insufficient to do this for the present gels. Thus, the
sure variation.15–18 In the long-wavelength limit, this pressure Kelvin–Voigt description is suited to our objectives since it
variationcancorrespondtoapassingacousticwave.Thebasic providesthesimplestdescriptionthatwillfacilitateexamina-
Rayleigh–Plesset formulation has been augmented to include tionofbubbledynamicsinourtissue-mimickingphantoms.
first-order liquid compressibility,19–21 which provides for
energy loss via acoustic radiation. This is important in some C. Gel-basedtissuephantoms
regimessowelikewiseincludecompressibilityeffects.
OneofthemainquestionsregardingtissueinjuryinSWL
Describing the corresponding bubble dynamics in soft
concernswhyittakessomany((cid:2)1000)shockwavestoinitiate
tissues,even in this simplified single-bubble case, introduces
significant tissue injury, as observed for example, by Bailey
additionalchallenges.Softtissueishighlydissipativeandcan
etal.6Thistimealsocorrespondstotheonsetoftheacousticsig-
have sufficiently strong elastic resistance to counter bubble
natures of cavitation in pig kidneys in vivo.6,41 Recent experi-
motionawayfromitsstress-freereferencesize.22–26Torepre-
ments,investigatingcavitationinpolyacrylamideandagargels
sentthesemechanisms,modelshavebeenfurthergeneralized
showasimilardelay(Figs.1and2).Theprincipalobservations
to viscoelastic material, often using the Kelvin–Voigt model
thatmotivatethepresentanalysisarethatbubblesbecomevisu-
foritssimplicity.Theinfinitesurroundingmediumwasmod-
allyactiveatmultiplethoughfewsitesduringtheinitialpulses
eled as a linear elastic material by Yang and Church27 and
andthatnewcavitiesbecomeapparent,eachsuddenly,overthe
generalized to include finite-deformation elasticity with a
courseofhundredsofbursts,similartowhatisdeducedinpig
neo-Hookean model by Gaudron et al.28 These models have
kidney tissues during a SWL treatment.6 Furthermore, once
been valuable for quantifying the relative importance of dif-
activethebubblesareallofroughlythesamesize;theyseemto
ferent mechanisms for confining bubbles within tissue,
share a similar maximum radius upon delivery of subsequent
despitetheirsignificantsimplifications.
acousticwaves.Toassesspotentialmechanismsleadingtothese
The two main additions to the Kelvin–Voigt model in
observations, we develop a Rayleigh–Plesset-like bubble
thepresentstudyarestrain-hardening, whichisknowntobe
dynamicsmodelwithmaterialpropertiesbasedonindependent
important in many tissues as well as in some gel models of
quantification of the material properties of agar and polyacryl-
tissues, and a fracture-like disruption based on reported
amide gels, including models for viscous dissipation, strain
resultsforothersoftmaterials.29,30Onlywiththeseadditions
hardening,materialfatigue,andfracture.
do the models we develop reproduce the principal observa-
tions from recent BWL experiments in gel-based tissue
D. Organizationofthispaper
phantoms.31 Because of its simplicity, the Kelvin–Voigt
model has been used to model the viscoelastic behavior of The overarching hypothesis guiding the investigation is
hydrogels in conjunction with several experiments.32–36 As that during the initial pulses, some combination of viscous
with any model, it simplifies the description and is under- and elastic resistance confines the bubble growth. A bubble-
stoodtobeinsufficienttopreciselydescribeallthemechani- dynamics model that incorporates material properties known
cal behavior of hydrogels. There are alternative descriptions to be important in similar gels is developed in Sec. II. Two
available, which potentially provide a better quantitative constitutivemodelsaremotivatedandimplemented,withand
model for some types of deformation, such as linear and without strain hardening as appropriate for each gel consid-
nonlinear Maxwell,22,23 upper-convected Maxwell,37 Zener, ered. Two damage mechanisms are proposed in Sec. IV: one
standard linear solid,38 Phan–Thien–Tanner,39 Jeffreys, isbasedonthewidelyusedGriffith’sfracturecriterion,appro-
Oldroyd-B, and Giesekus40 models. These are expected to priateforagargels;theotherisbasedonobservedmaximum
provide better representations of the viscoelastic behavior if allowed stretch ratio, appropriate for polyacrylamide gels.
J.Acoust.Soc.Am.140(2),August2016 Movahedetal. 1375
FIG. 2.Bubbles in agar gel. The
imaged region is 4.24cm(cid:3)2.48cm.
The BWL pulse has 20 cycles with
base frequency 335kHz, peak ampli-
tude of 7.51MPa and repeats at
200Hz.Imagesareshownat10cycles
intothepulsenumberaslabeled.
(cid:5)(cid:3) (cid:4) (cid:6)
When damage conditions are met, they predict rapidly grow- 2 @U @U @U
r ¼p þI B (cid:4) B B
ing damage volume, which is assumed to disrupt the gel’s e;ij ffiIffiffiffi @I 1@I ij @I ik kj
3 1 2 2
elasticandviscouscharacteristics.ThedrivingBWLpressure
pffiffiffiffi@U
pulseisbaseduponthewaveformmeasurementsdescribedin þ2 I3@I dij; (2)
Sec. V. It is shown in Sec. VI that with independently mea- 3
sured material properties, in conjunction with a standard where I , I , and I are the principal invariants of the left
1 2 3
materialfatiguemodelduetocyclicloading,thebasicexperi- Cauchy–Green tensor B¼FFT for deformation gradient F.
mental observations are reproduced. The bubble response to For an incompressible material, detF¼1, so I ¼1 and Eq.
3
BWL and SWL pulses is compared in Sec. VII. Conclusions (2)becomes
andimplicationsaresummarizedinSec.VIII.
(cid:5)(cid:3) (cid:4) (cid:6)
@U @U @U
r ¼2 þI B (cid:4) B B : (3)
e;ij @I 1@I ij @I ik kj
1 2 2
II. BUBBLEDYNAMICSWITHSTRAINHARDENING
We consider two models for U, both of which are com-
Thoughthe bubbles might appear to be denselypositioned monly used to represent soft tissue-like materials and are
in space in Figs. 1 and 2, this is in part an artifact of the two- motivated and parameterized for the specific gels we con-
dimensional projection of bubbles distributed within the trans- siderinSec.III:theneo-Hookeanmodelwith
parent three-dimensional material. The cloud is approximately
g
axisymmetric,sothe bubblesseen inthe images are distributed U ¼ ðI1(cid:4)3Þ; (4)
2
alongthelineofsightoveradepthcomparabletotheobserved
verticalextent.Evenattheirgreatestexpansionanddensity,void wheregistheshearmodulus,andtheFungmodel,43,44
fractions evaluated in an area of 1cm(cid:3)1cm corresponding to
g
the center of bubble activityinthe two-dimensional images are U ¼ expaðI1(cid:4)3Þ; (5)
2a
onlyabout6%.Wethereforeconsiderasinglesphericalbubble
in an infinite material, which before damage is homogeneous
where g and a are material properties. The Fung model Eq.
and isotropic. This will be an especially good approximation
(5) simplifies to the neo-Hookean model Eq. (4) in the limit
whenthebubblesareintheirnascentstagesofdevelopment.
of small (I (cid:4) 3), so g in the Fung model can be considered
1
The material stress is described by a nonlinear
the low-strain shear modulus, as described in more detail in
Kelvin–Voigtmodel,forwhichtheelasticr andviscousr
e v Sec. III and the Appendix. The strain energy function U is
stressessumas
onlyafunctionofthefirstprincipalinvariantI forthemod-
1
r¼r þr : (1) elsweconsider,soEq.(3)simplifiesto
v e
(cid:5)(cid:3) (cid:4) (cid:6)
The material is hyperelastic under the assumptions that (i) @U
r ¼2 B : (6)
theelasticresponseisreversible,(ii)thestressdependsonly e;ij @I ij
1
on the current strain, (iii) the material is isotropic, and (iv)
the material is effectively incompressible. In this case, a Both models have a simple form for the Cauchy stress
strainenergydensityUdeterminestheelasticstresstensor42 whenitisexpressedforanincompressiblematerial:
1376 J.Acoust.Soc.Am.140(2),August2016 Movahedetal.
r ¼gdevðBÞg(cid:4)pI; (7) wherep~isrelatedtothehydrostaticstresspby28
e
r þr þr
wherepisthehydrostaticpressure,and p¼(cid:4) e;rr e;hh e;//
3
" #
g¼(cid:7)1 neo–HookeanEq: ð4Þ (8) ¼p~(cid:4)g (cid:3)ro(cid:4)4þ2(cid:3)r(cid:4)2 g: (17)
expaðI1(cid:4)3Þ FungEq: ð5Þ: 3 r r
o
Theviscousstressrv ispurelydeviatoricandistakento The viscous contribution to the stress tensor based on
belinearlyrelatedtothestrainrate,so
Eq.(9)is
rv ¼2lS¼lðruþruTÞ; (9) @u R2R_
r ¼2l ¼(cid:4)4l and
v;rr @r r3
wherelistheviscosityanduisthevelocity.
u R2R_
Assuming spherical symmetry and neglecting the r ¼r ¼2l ¼2l : (18)
v;hh v;// r r3
bubble-interior low-density gas dynamics, the governing
equationscanbesimplifiedintheusualway18to Asappropriateforanearlyincompressiblemedium,dissipa-
(cid:3) (cid:4) tionassociatedwithdilatationisneglected.46
@u @u @ 2r (cid:4)r (cid:4)r
q @t þu@r ¼@rðrrrÞþ rr rhh //; (10) Integrating Eq. (10) from the bubble radius r¼R to
(Ref.18)r!1yields
wherFeoqr tihsetheeladsetnicsictyo,natrnidbuutiiosnthtoertahdeiasltrveesslotceintyso.r, symme- RR€þ3R_2 ¼pB(cid:4)p1ðtÞ(cid:4)4lR_ (cid:4) 2S (cid:4)E; (19)
2 q qR qR q
try dictates that F is diagonal, with components in spherical
coordinates45
which now includes the time-dependent far-field pressure
p (t). This is approximately uniform in the present applica-
@r r 1
Frr ¼@r and F// ¼Fhh ¼r ; (11) tionssinceacousticwavelengths((cid:6)4.4mm)arelongrelative
o o to the bubble size ((cid:6)0.1mm). In Eq. (19), S is the surface
tension, and the internal bubble pressure p is related to the
where r¼r(r , t) is the current position corresponding to B
o
bubblevolumeassuminganisentropicprocess
stress-free reference configuration r , R(R , t) is the bubble
o o
radiusgivenreference-configurationbubbleradiusR ,andR_
o (cid:3) (cid:4)3c
R
isitstimerateofchange.SincedetF¼1foranincompress- p ¼p þp o ; (20)
iblematerial, B v B;o R
@r (cid:3)r (cid:4)2 wherepvisthevaporpressure,c¼1.4istheratioofthespe-
¼ o ; (12) cific heats, and p is the equilibrium bubble pressure for
@r r B,o
o R¼R :
o
whichintegratesto 2S
p ¼p þ (cid:4)p : (21)
B;o 1 R v
rðr ;tÞ¼½r3þkðtÞ(cid:5)1=3; (13) o
o o
The final term in Eq. (19) includes all the effects of
where the condition kðtÞ¼R3(cid:4)R3o corresponds to having elasticity,
r¼Ratthebubblesurface.Thus,theradialvelocityis
ð1 (cid:5)2r (cid:4)r (cid:4)r (cid:6)
@r R2R_ E¼ (cid:4) e;rr e;hh e;// dr
u¼ ¼ ; (14) RðtÞ r
@t r2
" #
ð1 r4 r
¼ (cid:4)2g o (cid:4) gdr; (22)
andthefirstprincipalinvariantoftheleftCauchy–Greenten- r5 r2
RðtÞ o
sorBis
(cid:3)r (cid:4)4 (cid:3)r(cid:4)2 where ro ¼ðr3(cid:4)R3þR3oÞ1=3 from Eq. (13). By defining
I1 ¼ ro þ2 r : (15) b¼r/roandrecognizingthat
o
!
1 r3
From Eq. (7), the corresponding contribution to the elastic db¼ (cid:4) dr; (23)
r r4
stresstensoris o o
r ¼g(cid:3)ro(cid:4)4g(cid:4)p~; r ¼r ¼g(cid:3)r(cid:4)2g(cid:4)p~; Eq.(22)withk¼R=roðRÞ¼R=Ro becomes
e;rr r e;hh e;// r
o ð1
(16) E¼ (cid:4)2g½b(cid:4)5þb(cid:4)2(cid:5)gdb; (24)
k
J.Acoust.Soc.Am.140(2),August2016 Movahedetal. 1377
pressures in Eq. (18),it is neglected. A small amount of gas
is expected to be transferred into the bubble from the sur-
roundingmediumduetorectifieddiffusion,50,51whichisre-
absorbed in the time (5ms) between subsequent BWL
bursts.52 Thus, for our purposes, we neglect mass diffusion,
assumenochangeinR ,andfocusonviscouseffects,elastic
o
effects, and damage all of which seem to be of primary
importance for the present regimes and materials. This
assumptionisfurthersupportedbytheobservationstowhich
wecompare.Residualbubbles,suchaswouldbeexpectedif
significant gas diffused into the bubbles, are not seen
between the applicationof BWL burst, so anygas exchange
isexpectedtobesmall.
The sound speed of hydrogels are generally very close
to water; thus, we take c¼1480m/s. The surface tension of
hydrogels are likewise similar to water, so we take
FIG. 3.(Color online) Confinement effect of elasticity vs stretch ratio for
theFung—(red)andneo-Hookean---(green)models. S¼0.073N/m (Ref. 53). Equation (26) is solved with an
adaptivefifth-orderCash–KarpRunge–Kuttascheme,54with
initialconditions
which is a preferred form for its simplicity. For a neo- Rðt¼0Þ¼R ¼0:1lm and R_ðt¼0Þ¼0: (28)
o
Hookean material, Eq. (24) can be integrated directly to
obtain
g III. MATERIALPROPERTIESFORGELPHANTOMS
E¼ ½5(cid:4)4k(cid:4)1(cid:4)k(cid:4)4(cid:5); (25)
2 Hydrogels,suchasthetissuephantomsinFigs.1and2,
consist of a sparse polymer network filled with a solvent of
fortheFungmodelanadaptivequadratureisusedtoapprox-
mostlywater.Theirelasticbehaviorhasbeenstudied,andis
imate the integral numerically up to the double-precision.
understood to depend upon many factors, including concen-
The obvious and significant effect of the strain hardening in
tration and how the polymer chains are interlinked.
the Fung model can be seen by comparing E for different
Microindentation indicates Young’s moduli are 175kPa for
stretchratioskinFig.3.
the stiffer 2% agar and 30kPa for the 5% polyacrylamide
First-order compressibility effects (based on a Mach
number expansion in R_=c) are added following the gels.31,55Forlarger deformations, torsiontests show distinct
Keller–Miksisapproach20,27,47 strain-hardening for agar gel,56 which provides sufficient
data to estimate a¼1 in Eq. (5), as described in the
(cid:3)1(cid:4)R_(cid:4)RR€þ3(cid:3)1(cid:4) R_(cid:4)R_2 Appendix. In contrast, polyacrylamide gels display linear
c 2 3c stress-strainbehaviorbeforerupture,andarethereforerepre-
1(cid:3) R_(cid:4) R d sentedwiththeneo-Hookeanmodel.57,58
¼ 1þ Qþ Q; (26) The mass densities of both gels are close to that of
q c qcdt water. Given the densities of water (q¼1000kg/m3), agar
(q¼1640kg/m3), and acrylamide (q¼1130kg/m3) (Ref.
wherecisthesoundspeed,and
59),yieldsq¼1012.8kg/m3for2%agarandq¼1006.5kg/
2S 4lR_ m3for5%polyacrylamidegel.
Q¼pB(cid:4)p1(cid:4) R (cid:4) R (cid:4)E: (27) There are no reported measurements of viscosity at our
specific conditions, which are challenging because they
We confirm that first-order compressibility effects obtained include both high strains and high strain rates. However,
in this approach match those that would be obtained more thereissufficientdataavailabletoenableustoselectreason-
rigorously following the matched asymptotic analysis by able ranges of values to assess the potential importance of
YangandChurch27forbothneo-HookeanandFungmodels. viscous dissipation in Sec. VI. Assuming the same
In theabsence ofstrain hardening, which corresponds to the Kelvin–Voigtdescriptionweemploy,Cathelineetal.33used
neo-Hookeanmodel,thecurrentformulationEq.(26)recov- transient elastography measurements of shear-wave speed at
ers the model developed by Gaudron et al.28 for polyacryl- 400Hztodeduceviscosityl¼0.22Pa sforagar.Theselow
amidegel. frequency measurements are complemented by measure-
Although initial void sizes cannot be measured in the ments of shock attenuation in ballistic gelatin, which is a
currentgels,theyareexpectedtohaveadistributionbetween similar hydrogel tissue-mimicking phantom made of gelatin
0.01 and 10lm (Refs. 48 and 49). In Sec. VI, we will typi- powderandwater.Liuetal.34reportl¼1Pa sforthestrain
cally assume an initial radius of R ¼0.1lm, though the rates up to 5.2(cid:3)104s–1, also assuming a Kelvin–Voigt
o
conclusions are insensitive to this choice. Since vapor pres- model. These rates match those anticipated in cavitation,
sure will be small (p <0.1MPa) relative to the other although at smaller strains. Typically, higher viscosities are
v
1378 J.Acoust.Soc.Am.140(2),August2016 Movahedetal.
measured for higher strains,60,61 up to the failure point. isusedfortheagargel,64thoughagainourprincipalconclu-
Thus, we take these values as lower bounds for the present sionsarenotstronglysensitivetothisparticularvalue.
conditions. Other measurements suggest significantly higher Bythismodel,cracksinitiatewhenabubblegrowsfrom
viscosities. Dynamic indentation measurements suggest that its initial radius R [Fig. 4(a)] to the onset of fracture at R
o f
the viscosity for agar gel varies between l¼5 and 900Pas [Fig.4(b)].However,subsequentcrackgeometriescannotbe
for frequencies in the range 20–200Hz (Ref. 35). Similarly, describedjustbyanenergybalance.Sincenolongcracksare
large deformation torsion experiments suggest that the vis- obviousinthevisualization,weassumethatfracturedamage
cosity for polyacrylamide gel varies between l¼1 and is localized near the bubble, in some number of small dam-
300Pasforthefrequencyrange0.01–30Hz(Ref.36). agefeatures.Thelocaldetailsoffracturearelikelytodisrupt
InSec.VI,weconsiderthebehaviorforthiswholerangeof the spherical symmetry, though this is not prominent in the
lfrom0.001to10.0Pasforbothpolyacrylamideandagargels. visualization either, so consistent with our overall descrip-
Fortunatelyourconclusionsdonothingeonthespecificvalues. tion, we homogenize this damaged region into a D-thick
Indeed,assuminganyofthevaluesreporteddoesnotfundamen- spherical layerof damaged material beyond the bubble [Fig.
tally alter our principal findings. It should also be recognized 4(c)]. This assumption is expected to apply for the majority
that there is no expectation that a single viscosity parameter ofbubblesthatseemtoremainsphericalasinFigs.1and2.
couldfullydescribethedissipativebehaviorofthegels.Thekey Forthisgeometry,thestoredstrainenergyis
fact we rely upon is that the gel is expected to be significantly
dissipative, with an effective viscosity that far exceeds that of ðR
W ¼ E4pr2dr; (29)
water(l¼0.001Pas)inalldeformationregimes.
Ro
IV. DAMAGEMODEL andtheenergyreleaseasthecrackgrowsis
A. Criteriafortheonsetofdamage (cid:8)
@W(cid:8)
The two criteria utilized to specify the onset of damage (cid:4) @A (cid:8)(cid:8) ¼2R0FðkÞ; (30)
V
are introduced in the following subsections. The first crite-
rionisbasedonGriffith’stheorem,whichisanenergybased where A¼pR2 is the surface area of the bubble.29,30 With
thresholdmechanism, andappropriate, aswedescribeinthe k¼R/R representingthestretchratioatthebubblesurface,
o
following, for the agar gel. The other is an ultimate strain
!(cid:8)
limitthatisbaseduponthegelmicrostructure,whichisused @ ðb (cid:8)
to describe the polyacrylamide gel, which is likewise based FðkÞ¼k4@b b(cid:4)3 Ek2ddkd (cid:8)(cid:8) ; (31)
1 (cid:8)
onreportedresultsdiscussedsubsequently. b¼k
1. Agarfailure:Griffith’stheorem with k and b appearing as dummy variables. The energy
d
released increases with k, and fracture occurs when it
Elastic materials store energy when deformed, and this
reaches fracture toughness G. Thus, the resulting critical
energy provides the work that can reversibly return them to
stretchratiok ¼R/R is30
their initial configuration when applied stress is released. f f o
However, there is a finite amount of energy that a material G
(cid:9) (cid:10)
F k ¼ ; (32)
canstore,beyondwhichitfailsviafracture.Thislimitcorre- f 2R
o
spondstotheenergyrequiredtogeneratenewsurfaces(e.g.,
cracks), and the corresponding threshold energy between withEinEq.(31)embodyingtheelasticbehaviorofthespe-
these regimes constitutes Griffith’s criterion: a crack will cificmaterial.Foraneo-HookeanmaterialEq.(4),
form and propagate when the elastic strain energy released
by crack growth exceeds the surface energy of the new FðkfÞ¼gð1þk2f (cid:4)2k(cid:4)f 1Þ: (33)
crack.62,63 Fracture toughness G is the corresponding mate-
rialpropertythatquantifiesamaterial’sresistancetofracture. Thecorrespondingk (orR)fortheFungmodelisevaluated
f f
An accepted quasi-static fracture toughness of G¼6.5J/m2 numerically for any fracture toughness G, elastic parameter
FIG. 4.(Color online) Schematic of
bubble growth with a fracture-like
damagemodel.R istheinitialbubble
o
radius, R is the radius at which frac-
f
ture occurs initially, R is the bubble
radiusattimet,andDisthethickness
ofthedamagedmaterial.
J.Acoust.Soc.Am.140(2),August2016 Movahedetal. 1379
TABLEI.Thefracturestretchratiok basedonEq.(32)fortheFungmodel which can be solved for the time dependent damage thick-
f
Eq.(5)witha¼1,G¼6.5J/m2fordifferentshearmoduligandinitialcav- ness with respect to the current bubble radius. Thus, new
ityradiusR.
o damagethicknessis
g(Pa) 8
0; ifRðtÞ=R (cid:8)k ;
Ro(lm) 102 103 104 105 106 D0 ¼>>< R3(cid:4)R3!1=3 o f (35)
k o (cid:4)R; otherwise;
1.0 2.65 2.41 2.16 1.85 1.45 >>: f k3f (cid:4)1
0.1 2.86 2.65 2.41 2.16 1.85
0.01 3.05 2.86 2.65 2.41 2.16
where RðtÞ=R is the stretch ratio at the bubble surface. The
o
overall damaged volume is tracked in time and increased
g,andinitialcavityradiusRo.Someresultsforthisaretabu- whenever D0 >D, where D is the maximum damage thick-
latedinTableIforRo¼0.01to1lm.Overthislargerange, nessuptothecurrenttime.Theextentofdamageisthuspre-
Griffith’s theorem predicts only a modest variation of kf, servedfromoscillationtooscillation.
from kf¼1.45 to 3.05, in agreement with experiments for In the dynamic model, the lack of elastic resistance of
agar gel.30,56 The primary initial condition we consider, damaged material is accounted for by excluding the dam-
Ro¼0.1lm, yields kf¼2.2, so the onset of fracture is agedregionintheintegrationboundsofEq.(22):
R ¼0.22lm in most cases. Applying Griffith’s theorem
f " #
with fracture toughness as measured under quasi-static con- ð1 r4 r
E¼(cid:4) 2g o (cid:4) gdr: (36)
ditions, where viscous effects are negligible, constitutes a r5 r2
RðtÞþDðtÞ o
necessary assumption if we are to use available estimates of
fracture toughness in our dynamic model. Extension of this
Inaddition,theviscosityofthedamagedmaterialisreduced
to viscoelastic solids under dynamic loading would presum-
totheviscosityofthesolvent,inthiscasewater.Theeffects
ably allow a more accurate measure of the onset of fracture
ofthisdamageonbubbleevolutionarereportedinSec.VI.
andimproveourdamagemodel.
C. Fatiguedamagewithstresssoftening
2. Polyacrylamidefailure:Ultimatestrain
Cyclic large deformation of gels such as we consider is
All gels must fracture (or somehow fundamentally fail)
expected to lead to stress softening due to fatigue, which
whenstrainexceedsthemaximumallowedstretchkfoftheir weakens the material.68–70 We shall see in Sec. VI that in
constituent polymer chains. Measurements for polyacryl-
manycasestheinitialloadingisinsufficienttocausedamage
amide gel suggest that it is this mechanism that sets the
based on thresholds described in Sec. IVA. However, for
stretch ratio to be k ¼2.1 for polyacrylamide gel.65,66
f these same loadings, the material is expected to degrade
Applying Griffith’s theorem to polyacrylamide gel Eq. (33),
such that multiple pulses are required for material failure as
with G reported between 2 and 3.1J/m2 (Ref. 67) for
wouldcorrespondtotheobservedbehavior.
R ¼0.001to1lm,suggestsk between9.9and122.5.Thus,
o f Severalstudiesofhydrogelssupportsuchamechanism.
polyacrylamidegelisexpectedtofailbasedonitsparticular Martins dos Santos et al.70 quantified fatigue due to strain-
microstructurerather thanthe Griffith’stheoremenergy bal-
controlled loading in agar gels and observed that the maxi-
ance since its constituent molecules are disrupted before the
mum stress resistance to compression decreases by about
Griffiththresholdisreached.
50%formaximumstrainof22.8%,and75%oftheultimate
strain after 1000 cycles. Stress softening has also been
B. Damagezonethickness
reported for synthetic biopolymers used in medical devi-
Since we are particularly interested in how subsequent ces.71–73 Ultra-high molecular weight polyethylene shows
bubbledynamicsareaffectedbydamage,wealsoneedtorep- distinct stress softening due to cyclic loading.74 Similarly,
resent the effect of failure beyond its onset. For this, we fatiguehasbeenshownintestsofsofttissues.Schulmanand
assume that damage is localized to the neighborhood of the Bader75reporteda50%decrease inultimate tensilestrength
expanding void, and constrained by the same spherical sym- under cyclic loading for human extensor digitorum longs
metry of the bubble dynamics model. For both gel materials, tendons. Several studies also reported fatigue phenomenon
whenthebubbleradiusR(cid:7)R,weassumenewlyformeddam- incardiovasculartissues.76–79
f
age of thickness D0, measuredfrom the current bubble radius. Complicated phenomenological fatigue models have
The boundary between the damaged and undamaged zones is been proposed for fatigue damage of soft tissues with stress
thusatr ¼RðtÞþD0ðtÞ,andcorrespondstok¼k whendam- softening.80However, specific quantitativemodels for stress
f
ageisforming.UsingEq.(13),whichcanbeinvertedtodeter- softening under cyclic loading for agar and polyacrylamide
mine the referenceradiusr (r)for materialcurrentlyat radius are not currently available due to the lack of experimental
o
r,thiskðRþD0Þ¼k conditioncorrespondsto data.Still,thefactthattheywillfatigueiscertain,andfortu-
f
natelyasimplestandardmodelissufficientforourpurposes.
RþD0 RþD0 Fatigueundercyclicloadingisassumedtodecreasebothvis-
¼ ¼k ; (34)
roðRþD0Þ (cid:11)ðRþD0Þ3(cid:4)R3þR3o(cid:12)1=3 f cboyuasafancdtoerlaast,i7c5,8s0trewshseicshwiinthcrreeasspeescwttioththeeacuhndlaomadaignegdcsytcaltee
1380 J.Acoust.Soc.Am.140(2),August2016 Movahedetal.
where
pþ
max(cid:4)1
p(cid:4)
cðtÞ¼1þ max ð1þcosð2pftÞÞ; (40)
2
is introduced to account for the asymmetric peak amplitude.
The shown signal has n¼10 cycles, principal frequency of
f¼335 kHz, peak negative pressure of p(cid:4) ¼ 1MPa and
max
FIG. 5.(Color online) Driving pressure pulse: (a) SWL pulse with maxi- wasusedtosetthering-upandring-downtimeforthepulse
mum of 35MPa and minimum of (cid:4)10MPa; (b) BWL pulse with n¼10 envelope s ¼5ls and s ¼12ls, respectively, in Eq. (39).
1 2
cycles,f¼335kHzandp(cid:4) ¼1MPa;modeledbasedonEq.(39)—(red),
max The treatments of the gels in Figs. 1 and 2 and other animal
measuredexperimentally—(blue).
studies3 used the same BWL source but operated at a higher
from an initial a¼0 up to a¼1, at which point the gel intensity.Sinceitwasnotpossibletoprovideanaccuratemea-
affords no resistance, though in our cases the other damage sure of the signal at these higher pressures, extrapolation was
mechanisms become activated before a¼1 would be usedtoestimatep(cid:4) ¼7.51MPaandpþ ¼10.09MPa.Full
max max
achieved, switching from the slow fatigue to fast fracture details and justification of this procedure are reported else-
mechanism.Toincorporatefatiguemodeledinthiswayintothe where.Thenumberofcyclesinthecorrespondingexperiments
bubbledynamicsmodel,Eq.(27)ismodifiedtoincludeaas didnotshowsignificantsensitivitytonrangingfromn¼10to
20; our present conclusions are insensitive to our typical
2S (cid:3)4lR_ (cid:4)
Q¼p (cid:4)p (cid:4) (cid:4)ð1(cid:4)aÞ þE : (37) choiceofn¼10.
B 1 R R
VI. BWLMODELPREDICTIONRESULTS
Since fatigue is a probabilistic phenomenon, depending for
example, on the local polymer network at each nucleation A. Polyacrylamidegel
site, the specific rate at which a increases is expected to be
We first consider the polyacrylamide gel, which is sim-
different at each cavitation site leading to the typical scatter
plerinthesensethatitdoesnotshowstrainhardeningbefore
observed in fatigue experimental data.69 This is proposed
fracture, the fracture itself initiates based on a simple strain
herefortheseeminglyrandomappearanceofnewbubblesas
criterion, and the elasticity is relatively weak. Given the
variabilityin thereducingconfinement due to fatigueallows
uncertainty of the viscosity discussed in Sec. III, we also
fracturedamagetooccuratdifferentsitesatdifferenttimes.
consider a wide range of viscosities between 0.001 and
10.0Pas. The maximum radius reached for this range is
V. DRIVINGPRESSUREPULSE
shown in Fig. 6(a). We see that the damage criterion
Themodelweconsiderisgeneral,andcouldbeusedfor k>kf¼2.1, is not satisfied for l>lc¼1.8Pas, so even a
any strong pressure signal. We consider two waveforms: a viscositythatislowcomparedsomemeasurementsforpoly-
standard SWL pulse and a BWL pulse. The SWL profile is acrylamide suppresses the bubblemotionsufficientlytopre-
showninFig.5(a)andrepresentedby clude damage. Below this critical value, however, the
damageissignificantandbubblesgrowtonearly1000times
p ðtÞ¼p þ2pþ exp(cid:4)ctcos(cid:3)xtþp(cid:4); (38) their initial size. Elastic confinement plays little role in this
1 amb max 3 case,whichcanbeanticipatedsinceeventhemaximumelas-
tic confinement at k, as calculated from Eq. (25), is
wherep ¼0.1MPaistheambientpressure,pþ ¼35MPa f
amb max E¼15.2kPa,whichisonlyabout1/500timesthepeaknega-
isthepeakpositivepressure,c¼9:1(cid:3)105s(cid:4)1 andx¼p=4
tiveBWLpressure.
(cid:3)106s(cid:4)1, which were chosen to match the waveform in a
Illustrative R(t) histories are shown in Fig. 6(b). For
DornierHM3lithotripter.81,82TheBWLpulseisbasedupon
l¼0.01Pas,viscosityhaslittleeffectandthebubbleradius
focal pressure measurements of a BWL transducer using a
fiber-optichydrophoneinadegassedwaterbath.64Themea-
suredBWLpressurepulseisshowninFig.5(b)andiswell-
fittedby
p ðtÞ¼p þp(cid:4) cðtÞcosð2pftÞ
1 amb max
2 tf
(cid:3)4(cid:3)1(cid:4)exp(cid:5)(cid:4) t(cid:6)(cid:4)1þtanh2
s 2
1
0 2n 31 3
(cid:4)t
(cid:4)@1(cid:4)exp4f 5A1þtanhðtf (cid:4)nÞ7; FIG.6.(Coloronline)Predictedbubblebehaviorinpolyacrylamidegel:(a)
s 2 5 the maximum bubble radius achieved during one BWL pulse for a wide
2
rangeofpossiblegelviscosities(seeSec.III);and(b)thebubbleradiushis-
toryforl¼0.01Pas(cid:9)(cid:9)(cid:9)(blue),l¼1.0Pas– –(red),andl¼3.0Pas—,
(39)
therectanglescorrespondtotheonsetofdamage.
J.Acoust.Soc.Am.140(2),August2016 Movahedetal. 1381
pixel-size) defects in the gels. By the tenth pulse, both radii
reachalong-timebehavior.We note,however, thatthefinal
R (cid:6) 554R after 25 pulses is still smaller than R
max o max
(cid:6) 1026R reported for l¼0.01Pas in Fig. 6, which indi-
o
cates thatthehighviscosityofthesurroundingmediumout-
sidethedamagedzonestillprovidessomeconfinement.
B. Agargel
Wecananticipatethatadditionalfactorswilldictatethe
FIG.7.(Coloronline)Fatiguebehaviorofpolyacrylamidegel:(a)Thecriti-
behavior of the agar gel. At k ¼2.22, as indicated by
calstresssofteninga fordifferentgelviscosities:bubbledynamicsmodel(cid:10) f
c
(red)and––(blue)Eq.(41);and(b)Themaximumbubbleradius—(red), GriffithcriterionindicatesinSec.IVA1,themaximumelas-
anddamagezonethickness––(blue)forl¼3.0Pasanda¼0.4afterstress tic confinementfrom (24)is E¼3.22MPa, which iscompa-
softeningwithexposures. rable to the negative pressure peaks of BWL pulse p(cid:4)
max
¼ 7.51MPa. Thus we expect elasticity to play a role in the
R(t)increasesbyafactorofabout1000almostimmediately.
bubble response. The bubble dynamics are again considered
For l¼1.0Pas, the bubble growth is initially suppressed. forarangeofviscosities.Forl¼0.01Pa selasticityiscapa-
However, during the second cycle, the damage model acti-
ble of confining the bubble growth sufficiently to preclude
vates, which reduces subsequent resistance and allows bub-
damageduringthefirsttwocycles(Fig.8).However,subse-
blestogrowstilllargerandcausestillmoredamageoneach quently R(t) does exceed R leading to run-away damage.
f
subsequent cycle until it is essentially unconfined. For Interestingly, even for the much larger l¼1.0Pas, also
l>lc¼1.8Pas, no damage occurs, and the bubble stays showninFig.8,thedamagemodelalsoactivatesyetbubble
stronglyconfined,asshownforl¼3.0Pas.Thisisthemost growth remains modest, with R /R ¼12.1. The stiffness
max o
likelyscenariobasedontheviscosityestimatesforthepoly- and particularly the strain hardening of the agar are thus
acrylamidegelsdiscussedinSec.III. protective. Finally, for l¼3.0Pas the radius R(t) stays
Still,evenwhentheconfinementissuchthatthebubble belowthedamageonsetthreshold.ThecorrespondingR
max
remains small, fatigue-induced material degradation is for the full range of viscosity is also plotted in Fig. 8(b),
expected upon repeated BWL pulses, as discussed in Sec. and it is clear that the requisite viscosity l ¼1.5 Pas to
c
IVC.Sinceelasticcontributiontothetotalstressissmallfor suppress run-away damage and bubble growth is smaller
polyacrylamidegel,thecriticalaccanbeestimatedtobe than for polyacrylamide gel. The corresponding amount of
strain softening necessary for fatigue to initiate fracture-
l
ac ¼1(cid:4) lc; (41) likefailureisshowninFig.8(c).Theincreasedelasticcon-
finement for agar gel causes deviation from viscosity-only
modelEq.(41).
whichisplottedinFig.7(a)along withresultsfromthecor-
responding numerical solution, which includes elastic
VII. SWLANDCOMPARISONWITHBWL
effects. When this critical a is reached, presumably within
around1000pulsesbasedonthegelfatigueexperimentsdis- In this section, we investigate the model response to
cussed in Sec. IVC, significant growth occurs as shown for theSWLpulseandcompareitwithBWLresults.Themaxi-
l¼3.0Pasanda¼0.4inFig.7(b).Toconstructthisplot,a mumbubbleradiusachievedduringonepulseforarangeof
greaterthana ischosensuchthatk¼1.02k isreacheddur- possible effective viscosities between l¼0.1 and 10.0Pas
c f
ing the first subsequent exposure, which leads to rapidly isshowninFig.9.Asmallermaximumbubbleradiuscorre-
increasing damage and unconstrained bubble growth after sponds to a smaller damage thickness within our damage
onlyaboutsevenmoreexposures.Thisisconsistentwiththe model, and might indicate less potential tissue damage
rapid appearance of bubbles from unseen (presumably sub- in clinical lithotripsy treatments as well. For SWL, the
FIG.8.(Coloronline)Predictedbubblebehaviorinagargel:(a)thebubbleradiushistoryforl¼0.01Pas(cid:9)(cid:9)(cid:9)(blue),l¼1.0Pas––(red),andl¼3.0Pas—
,therectanglescorrespondtotheonsetofdamage;(b)themaximumbubbleradiusachievedduringoneBWLpulseforawiderangeofpossiblegelviscosities
(seeSec.III);and(c)thecriticalstresssofteninga vsviscosity;bubbledynamicsmodel(cid:10)(red)and––(blue)Eq.(41).
c
1382 J.Acoust.Soc.Am.140(2),August2016 Movahedetal.
(i) polyacrylamide(linear):neo-Hookeanmodel,and
(ii) agar(strainhardening):Fungmodel;
• criteria for onset of damage based on reported results for
eachmaterial:
(i) polyacrylamide:ultimatestrain,and
(ii) agar:Griffith’sfracturetheorem;
• incorporationofdamageinthebubbledynamics;and
• stress softening based upon a simple standard approach
forfatigueundercyclicloading
FIG.9.(Coloronline)Themaximumbubbleradiusforpolyacrylamideand
The resulting model reproduces the key features of the
agargelsforarangeofpossibleeffectiveviscosities.
experimental observations. It suggests that a combination of
criticalno-damageviscosities arel ¼4.67and3.67Pasfor viscousandelasticresistanceiscapableofconfiningthebubble
c
polyacrylamide and agar gels, respectively. The larger l growth. For polyacrylamide, the main mechanical resistance
c
compared to BWL suggests that a for a given viscosity comesfrom viscosity,whileelasticresistance issignificantfor
c
agar due to its strain hardening. Above a critical viscosity,
basedonEq.(41)islowerfortheSWLpulse.Consequently,
withintherangeofthosereportedforthesematerials,thebub-
lessfatigueisneededforSWLcomparedtoBWL.Whilethe
blegrowthissuppressedsuchthatthedamagemodeldoesnot
current results suggest that fewer pulses are needed to initi-
activate,andconsequentlynoimmediatedamageorlargebub-
ate damage due to fatigue in SWL, damage in a treatment
blegrowthoccurs.Fortheserangeofviscosities,theamountof
will depend on many factors, including specific amplitudes
fatigue driven stress softening required to activate the damage
andnumber ofpulses neededfor stonecomminution.A bet-
modelispredicted,andverifiedbythebubbledynamicsmodel.
terunderstandingofBWL’sabilitytobreakkidneystonesis
The overall mechanism, consistent with the observations
thereforerequiredtomakeafaircomparison.
thatmotivatedthismodel,issummarizedinFig.10.Initially,
some nuclei are large enough, or defects make the material
VIII. SUMMARYANDADDITIONALDISCUSSION weakenough,suchthatsomesmallnumberofbubblesappear
with the first burst. Subsequently, bubbles appear after cyclic
A high-speed camera showed particular bubble activity
fatigue increases a to a . Once this occurs, a fracture mecha-
insidepolyacrylamideandagargelssubjectedtoBWLpulses: c
nism becomes activated. Unlike the fatigue, which is a slow
• bubbleactivationatmultiplesitesduringtheinitialpulses, accumulation of modest damage, both of the fracture criteria
with we introduced in Sec. IVA lead to an almost immediate and
• further continuous appearance of suddenly large bubbles complete weakening of the material. This process completes
atnewlocationsthroughthecourseoftheobservations to its ultimate bubble size after only a few cycles, but in this
same period the bubbles become nearly unconfined. The sto-
ARayleigh–Plesset-typebubbledynamicsmodelwitha
chastic character of the fatigue mechanisms explains their
damagemechanismisdevelopedtoevaluatepotentialmech-
continuous appearance as different sites reach the critical
anisms of viscous confinement, elastic confinement,
pointatdifferenttimes;therapidityofthefracturemechanism
fracture-like damage and fatigue in both gel phantoms. The
explains the apparent sudden (seemingly immediate) appear-
mainfeaturesofthemodelinclude:
anceofnewlylargebubblesduringatreatment.
• constitutive models based on independent experimental Despitetheapparentsuccessofthismodelatexplaining
dataforthetwomaterials: these basic observations, it should be reiterated that this
FIG. 10.(Color online) Schematic of
bubble activation process by focused
ultrasound bursts in the presence of
fracture-likeandfatiguedamagemech-
anisms.The mean and maximumbub-
ble radii in the experimental images
such as Figs. 1 and 2 (in Ref. 31) are
about60and120lm,respectively.
J.Acoust.Soc.Am.140(2),August2016 Movahedetal. 1383