Prediction of far-field acoustic emissions from cavitation clouds during shock wave lithotripsy for development of a clinical device

T. G. Leighton, C. K. Turangan, A. R. Jamaluddin, G. J. Ball, P. R. White


This study presents the key simulation and decision stage of a multi-disciplinary project to develop a hospital device for monitoring the effectiveness of kidney stone fragmentation by shock wave lithotripsy (SWL). The device analyses, in real time, the pressure fields detected by sensors placed on the patient's torso, fields generated by the interaction of the incident shock wave, cavitation, kidney stone and soft tissue. Earlier free-Lagrange simulations of those interactions were restricted (by limited computational resources) to computational domains within a few centimetres of the stone. Later studies estimated the far-field pressures generated when those interactions involved only single bubbles. This study extends the free-Lagrange method to quantify the bubble–bubble interaction as a function of their separation. This, in turn, allowed identification of the validity of using a model of non-interacting bubbles to obtain estimations of the far-field pressures from 1000 bubbles distributed within the focus of the SWL field. Up to this point in the multi-disciplinary project, the design of the clinical device had been led by the simulations. This study records the decision point when the project's direction had to be led by far more costly clinical trials instead of the relatively inexpensive simulations.

1. Introduction

During shock wave lithotripsy (SWL), thousands of shock waves are directed into the patient at a rate of about one per second in order to fragment kidney stones or reduce them to a size whereby they can subsequently be dissolved using drugs [1,2]. With current apparatus, the clinician is ill-equipped to determine in-theatre whether the treatment has been successful, with the result that 30–50% of patients need to return for re-treatment, and an unknown number receive a greater exposure to shock waves than is necessary for stone fragmentation [2]. Overexposure carries the potential for adverse side effects [36]. There is therefore an imperative for the development of a non-invasive device for monitoring the progress of the therapy in real time in the clinic [7]. This study describes the simulation that supported the development of a new passive acoustic sensor, which is placed by a nurse on the patient's skin, and passively monitors the scattering and reverberation within the patient's body of the pressure fields generated by the interaction of the SWL, bubbles and tissue.

Interactions between shock waves and bubbles include cases where the bubbles are collapsed by shock waves, cases where bubbles emit shock waves (either by generating them on collapse, or scattering pre-existing shocks), and cases that include all these phenomena. Such a case is studied here, simulating the responses of bubbles subjected to shock waves generated by SWL [1,2]. The technique is however equally applicable to scenarios involving underwater explosions or industrial erosion, which might be undesirable [8] or required, as with ultrasonic cleaning [9,10]. One feature of these studies is that, while simulations of such events usually focus on the responses of single bubbles, the practical phenomena often involve clouds of interacting bubbles, and the collective response of the cloud might be key to the resulting erosion that is generated [1113].

This study forms the third part in a series of simulation studies [14,15] that complement a series of laboratory and clinical studies [1618], which together form a coherent research programme aimed at the development of a passive acoustic device for monitoring the effectiveness of SWL. Following successful clinical trials where its automated response proved to have 94.7 per cent sensitivity in predicting a successful patient outcome, compared with 36.8 per cent achieved by the clinician performing the treatment [18], this device is currently being used for real-time assessment of treatment in the clinic (a typical lithotripsy treatment lasts about an hour). It is also being tested for its ability to predetermine whether a kidney stone will be susceptible to SWL or whether the patient should be referred instead for some other form of stone treatment [2,19].

The first stage of the work developed a free-Lagrange simulation for single bubbles, which offered advantages over the existing boundary-element method of capturing blast waves, and over the Eulerian methods, in that it could follow the bubble collapse after the moment when the liquid jet formed by bubble involution impacts the downstream bubble wall [14,20,21]. Inclusion of the events after this moment into the simulation is critical because the evolution and detection of the blast wave generated by that impact (and its subsequent interaction with the kidney stone and tissue) is important both to kidney stone fracture and also to the operation of the device, which detects this blast wave in the far field via a receiver placed on the patient's torso.

However, computational limitations meant that the simulation was not capable of extending the computational domain out to the far field in order to simulate the signal such receivers would detect. Therefore, the second stage of the work developed techniques for predicting the pressure field generated in the far field by the interactions of the incident lithotripter shock, the blast wave and signals reflected from the stone [15]. All of the work in the first and second stages simulated single bubbles. This third paper extends those methods to multiple bubbles because, in most practical cases, the erosive effect and the far-field emission will depend on the response of a cloud that can not only generate a summation of single-bubble effects, but can also generate important effects through the interactions between bubbles (such as cooperative bubble effects in the generation of erosion) [2224].

The degree of influence of a bubble on its neighbours is dependent on their separation distance, size, composition clustering, the incident pressure field and the location and composition of any boundaries present (such as a kidney stone might represent) [2533]. Recently, Lauer et al. [34] simulated the collapse of two-dimensional arrays of cavities by a planar shock wave using the so-called ‘conservative interface interaction’ method where the fluid–fluid interface was tracked using a level-set approach to study the collapse mechanisms of arrays of cavities and the effects of separation distance. They compared their predictions with those of Ball et al. [20] and the experimental data of Dear & Field [26] and Bourne & Field [35].

Here, simulations of the response of an array of air bubbles to a lithotripter shock wave using the free-Lagrange numerical code Vucalm are presented. The first part, which is presented in §3a, is a study of the interaction of a lithotripter shock with two stable spherical bubbles. The aims of this study are to observe the reflection and transmission of the shock waves, the collapse of the initially spherical bubbles following their interaction with the lithotripter shock wave, the formation of the high-speed liquid jet and its velocity–time history, the effect of the bubble separation distance and degree of influence of the neighbouring bubble and the bubble wall position–time history. The second part, which is presented in §3b, is a study to predict the far-field acoustic signature (pressure wave) emitted from a cloud of bubbles. The reason for making this prediction is to assist in the interpretation of the far-field acoustic signature from a cloud of cavitation bubbles generated within the patient during clinical SWL. The implications that these findings had on the direction of the project, and the clinical relevance of the resulting decisions, are presented in §4.

2. Numerical schemes

(a) The free-Lagrange method

The governing equations of the free-Lagrange numerical code Vucalm are derived by applying the conservation laws for continuous media in a Lagrangian reference frame. The resulting unsteady and compressible system of equations in a swirl-free axisymmetric form is Embedded Image2.1The conserved variables vector Embedded Image, flux vector Embedded Image and source terms vector Embedded Image are given by Embedded Image2.2where A(t) is the computational cell area, S(t) is the computational cell boundary length, and both are functions of time t. Variables rc and rm are the radial distances of the cell centre and the boundary middle point from the axis of symmetry, respectively. The unit vector of the cell boundary is Embedded Image, ρ is the fluid density and Embedded Image is the fluid velocity vector (i.e. Embedded Image where ux and ur are the velocity components in the axis of symmetry (x) and radial (r) directions, respectively). Variable E is the energy per unit mass (internal and kinetic energy), Embedded Image is the unit tensor, p is the absolute pressure, Embedded Image is the unit radial vector and Pf is the pressure acting on a control volume associated with the axisymmetric formulation (details of the derivation of the governing equations are presented in Ball [36] and Turangan [37]). The superscript represents the transpose of the vector. The numerical scheme uses a Voronoi mesh, where each computational cell represents a control volume bounded by cell boundaries (figure 1). To calculate the numerical fluxes across the cell boundaries and update the flow properties (density, velocities, temperature and pressure), approximate Riemann solvers (Godunov-type) are used [20,36,38].

Figure 1.

Voronoi mesh used in Vucalm simulations. Each cell represents a control volume bounded by cell boundaries. The cell encloses a cell-centred particle. The target particle has neighbour particles that share common cell boundaries. All particles are assigned with a fluid type, thermodynamic properties, coordinate and flow condition. The material interface is determined by the cell boundaries shared between dissimilar fluids, e.g. ‘fluid 1’ and ‘fluid 2’.

As convective flux across the control volume boundaries (e.g. velocity) is absent in Lagrangian schemes [36], physical features such as material interfaces are sharply resolved. This quality makes Lagrangian schemes suitable for multi-phase flow simulations. In free-Lagrange schemes such as Vucalm, the severe mesh tangling associated with conventional Lagrangian schemes that use fixed mesh connectivity is prevented by allowing the computational mesh to connect freely. In order to suppress mesh-induced high wavenumber instabilities associated with free-Lagrangian schemes that produce non-smooth dissimilar fluid interfaces, an artificial interface smoothing algorithm has been implemented in Vucalm [39]. With this algorithm, both the amplitude and growth rate of small-scale perturbations of the mesh-induced wavenumber instabilities are significantly reduced.

(b) The Kirchhoff acoustic emission scheme

The available computational resources do not allow calculation of the far-field pressures that would be detected by sensors because they will not enable the computational domain to extend to 10 cm from the lithotripter focus. Pressures at such distances for propagation through water can readily be calculated using the Gilmore model [40], although this study highlights that, when doing so, it is important to assess whether the assumptions inherent in the Gilmore model are germane to the scenario being simulated (see §4). Even for propagation in water, if the bubble collapses asymmetrically as in this study, the source of the pressure wave is incorrectly modelled by the Gilmore model as being the spherically compressed gas, rather than the blast wave [15,41]. Furthermore, propagation through tissue generates higher absorption than occurs in water. There are standard methods for derating predictions made for propagation in water to give estimates of what form the pulse would take if it propagated through tissue, and these have been applied to this problem [15,41]. However, the appropriateness of using the standard method to correct for absorption by the intervening tissue has not been addressed for cases such as the shocks presented here, and this will be also discussed in §4.

Apart from the free-Lagrange studies cited above and the work by Hawker & Ventikos [42] that used a high-resolution front-tracking method to study a shock–bubble interaction, earlier simulations had also included loss of sphericity and jets, but not compressibility or bubble fragmentation. Consequently, they operated on a time scale that was restricted to the early periods of the cavitation event, before these effects became important [4345]. The advantage of the free-Lagrange technique developed for this project [14,15,21] is that it can capture all these aspects, and follow the collapse through to the formation of the blast wave (which dominates the far-field emission from such a cavitation) and beyond. It is however computationally intensive, and with a limited computational domain, the current resources allow only approximate estimates to be made of the far-field emissions, with clear understanding of the ways by which such methods fall short of a full calculation [15]. A control surface can be placed around the bubble, and the pressures and densities on that surface are used to calculate the emissions to far field assuming linear propagation outside of that surface. However, account must be taken of absorption in the tissue and the effect of neglecting any nonlinear propagation in this way. Providentially, tissue absorption is the simpler feature to include if standard techniques are followed and are valid (a question discussed further in §4), and its inclusion reduces the importance of nonlinear propagation because high absorption rapidly attenuates the amplitude of the wave [15,41].

Following the arguments of Jamaluddin et al. [15], in this study, therefore, no account will be taken of nonlinear propagation outside the control surface. This compromises the validity of the estimated far-field pressures but records the principle for such a time when computational resources allow the extension of the control surface such that only linear propagation occurs outside of the control surface. Such a compromise would be more difficult to justify were the propagation to occurr in a material with low acoustic absorption and a high degree of nonlinearity, such as water. Justification for this necessary limitation is discussed by Jamaluddin et al. [15] for propagation in the more highly absorbing human body tissue that occurs when the device is used in the clinic. While limited resources necessitate this approach, it was sufficient to allow development of the device that was then validated as useful in the clinic through empirical clinical trials [16,17].

There are two methods for estimating the far-field emissions following linear propagation outside the control surface, and previously two numerical schemes based on the Kirchhoff and Ffowcs William–Hawkings formulations were developed and tested for this SWL problem [15,41]. The Kirchhoff method assumes that nonlinear waves are confined within the control surface, and the input parameters for the calculation of far-field acoustic emissions are compatible with linear wave propagation. The Ffowcs William–Hawkings method does not require the flow to obey the linear wave equation at the control surface. However, it is more difficult to implement because it requires a volume integral calculation of the input parameters. Tests of the two options for predicting the far-field emissions for single bubbles in the lithotripsy case studied here showed that the computationally simpler approach (the Kirchhoff method) gave results close enough for design of the eventual hospital device [15,41], and so only that method is used in the multi-bubble calculations presented in this study.

The Kirchhoff scheme requires time histories of flow variables, i.e. the pressure p and its normal and time derivatives (∂p/∂n and ∂p/∂t, where n is the normal vector on the Kirchhoff control surface under consideration) are recorded on a control surface that surrounds the acoustic source (in this case, the bubble collapse; see Jamaluddin et al. [15,41] for details). These variables are calculated directly from Vucalm simulations and become the input parameters for the Kirchhoff scheme to approximate the far-field acoustic emissions. The scheme assumes that the linear wave equation is valid outside the control surface, which requires that the control surface must be large enough to contain the region of all nonlinear behaviour. In order not to delay the development of the device until sufficient resources become available to extend the computational domain this far, a method of providing suitable estimates was developed and tested whereby the error associated with neglecting nonlinearity is reduced sufficiently by taking into consideration the absorption in the tissue [15,41].

3. Simulation results

(a) Lithotripter shock wave interaction with multiple bubbles

 Figure 2a shows the first example scenario chosen from the several cases simulated in this study. It comprises two spherical air bubbles immersed in water in a free field, which interact with a planar lithotripter shock wave. The water is represented by the Tait equation of state (EOS) and the air bubbles by the ideal gas EOS. At the present setting, the free-Lagrange code requires a level of precision in its input values for absolute pressure p, temperature T and density ρ that does not reflect the accuracy to which these parameters could measurably be said to be constant throughout a tissue sample. Consequently, the water and air bubbles are stated to be initially at the following (high-precision) standard atmospheric conditions: pw=pb=101 325 Pa, Tw=Tb=288.15 K, ρw=1000 kg m−3 and ρb=1.2246 kg m−3. The subscripts ‘w’ and ‘b’ refer to the water and air bubble, respectively. The specific heat ratio for the air bubble, γ, is 1.4. The initial air bubble radius, R0, is 0.04 mm, which is based on the interpretation of passive acoustic in vivo patient data using the Gilmore model [46]. The lithotripter shock profile shown in figure 2b follows the ideal lithotripter shock profile given by Church [40].

Figure 2.

(a) Two spherical air bubbles separated by a distance L are collapsed by a lithotripter shock wave. Figure is not to scale. (b) Time history of the lithotripter shock used in the simulations, which is based on the idealized lithotripter shock profile given by Church [41], i.e. Graphic where P+=90 MPa is the compressive wave pressure amplitude of the lithotripter shock wave, αd=9.1×105 s−1 is a decay constant, ω=2πf is the radial frequency, f=83.3 kHz is frequency and t is time.

All elapsed times in the simulations are measured from the moment the lithotripter shock impacts the first bubble (bubble 1). The lithotripter shock is introduced by imposing a time-dependent pressure boundary condition on the left computational boundary. The top and right computational boundaries are non-reflecting at all times. The separation distance between the initial centres of the bubbles is denoted by L. Simulations for four different separation distances, i.e. L={0.085,0.09,0.10,0.20} mm, were studied and are discussed.

The results for the lithotripter shock wave–bubble array interaction problem in the free field are shown in figures 3 and 4. The bubbles are separated by a distance of 0.09 mm. At t=0.111 μs, the incident shock has traversed the two bubbles (figure 3a). The interaction between the lithotripter shock and expansion waves originating at the bubble surface results in weakening and curvature of the shock. The dynamic behaviour of the collapse of bubble 1 is nearly identical to the problem for a single air bubble in the free field [14]. Bubble 1 is collapsed by the shock wave, and a strong air shock propagates in bubble 1, whereas a weak pressure wave is transmitted in the air of bubble 2. During the time that bubble 1 is collapsing, bubble 2 has been shielded from the initial incident shock wave and has experienced only a slight lateral compression (see figure 3c onwards). The liquid jet begins to form in bubble 1 in figure 3b, developing in amplitude by figure 3c. The jet travels across the bubble and reaches the downstream wall at approximately t=0.183 μs (figure 3d). At this time, bubble 2 shows no sign of liquid jet formation.

Figure 3.

(af) Pressure contours from an array of two air bubbles impacted by a lithotripter shock. Contour intervals in the water, bubble 1 and bubble 2 are denoted with ΔPw, ΔPb1 and ΔPb2, respectively. The left wall of bubble 1 involutes to form a high-speed jet (whose direction is indicated by an arrow on frames (c) and (d)) that pierces the bubble and impacts the opposite wall emitting a blast wave (the leading edge of which is at the tips of the arrows in frames (d) and (e)). The dashed lines indicate the initial locations of the surface of the bubbles and the black and grey dots on the axis of symmetry indicate their centres. Bubble separation distance, L, is 0.09 mm.

Figure 4.

(af) Pressure contours from an array of two air bubbles impacted by the lithotripter shock described in the text, showing the collapse of the bubbles, the formation of a high-speed jet of bubble 2, blast emission when the jet impacts the opposite wall of the bubble and interaction of the collapse of bubble 1 and bubble 2. Horizontal arrows in frames (c,d) show the direction of the movement of bubble 2's jet, and the tip of each arrow is vertically aligned with the tip of the jet. Bubble separation distance, L, is 0.09 mm.

At t=0.191 μs (figure 3e), the air cavity of bubble 1 is drawn into the vortex core, while the blast wave continues to propagate outwards radially from the bubble. The blast wave arising from the liquid–liquid jet impact of bubble 1 hits the upstream wall of bubble 2, causing it to collapse and produce a jet (figure 4a,b). On impact, the strength of the blast wave is calculated to be approximately 0.5 GPa. This is much greater than the amplitude of the incident lithotripter pulse, although its amplitude will decay as it propagates away from the bubble. In a cavitation cloud, such a monotonic decay is offset by the generation of new blast waves if this one causes the collapse of neighbouring cavities, a phenomenon related to the use of cavities to sensitize explosives [47,48]. A strong air shock is also transmitted into bubble 2. The strength of the blast wave decreases as it propagates into the surrounding water. It is clear that the collapse of bubble 2 is greatly amplified by the blast wave originating from the collapse of bubble 1.

By t=0.239 μs (figure 4c), bubble 1 has expanded to a volume greater than the collapsing bubble 2. The pressure gradient near the upstream of bubble 2 continues to increase with time. A high-speed liquid jet is formed and impacts on the downstream wall of bubble 2 at 0.239 μs. As a result of the lateral compression experienced by bubble 2 earlier in the collapse process, the liquid jet that is formed is narrower than that of bubble 1. On impact, the jet produces an intense blast wave in the surrounding water (figure 4d). This blast wave will interact with the expanded bubble 1, and cause the latter to undergo a secondary collapse. At t=0.248 μs, the jet has penetrated through the bubble, isolating a lobe of trapped air and highly compressed gas that resembles a tear-drop (figure 4f). This is in agreement with experimental observation by Bourne & Field [49], who went on to examine the influence of particles on jet direction [50].

By taking the physical insight provided by these simulations with high-speed photography of the collapse of arrays of bubbles [26], it can be summarized that if a third bubble were to be positioned downstream of bubble 2, a chain reaction would occur, and the third bubble would collapse in a similar manner to the collapse and rebound of the second bubble. These situations for bubble collapse and jet formation are likely to take place in many practical circumstances where shock waves are generated in bubble clouds (including lithotripsy, ultrasonic cleaning, etc. [13,51]) because pressure waves from the collapse and rebound of some bubbles will pass over neighbouring bubbles. Therefore, it is necessary to carry out a study to investigate how the collapses of neighbouring bubbles are affected by their mutual interactions. It is shown here, that for L=0.09 mm, bubble 2 is shielded from the incident lithotripter shock wave. However, the blast wave originating from the collapse of bubble 1 interacts with bubble 2. This is analogous to an incident shock wave passing over bubble 2, causing a jet in the direction of the shock.

These trends conform with the high-speed Schlieren images taken by Dear & Field [26] (figure 5) where air discs cut into a gelatin slab are acted upon by the incident shock wave (S, of about 0.1 GPa amplitude) that approaches from the left, generating refraction waves (visible in frame 1) and causing bubble involution and the formation of a jet (J, labelled in frame 2). The middle cavity is not collapsed by the incident shock S, from which it is partially shielded by the first bubble (by frame 3, it has experienced only a slight lateral compression). Instead, the middle bubble is collapsed by the blast wave S′, which is emitted when the jet in the first bubble has crossed it and impacted on the downstream wall of bubble 1. The third bubble is partially shielded from S by the second bubble, and collapses only by the blast emitted by the second bubble when its jet impacts its downstream wall. Dear & Field [26] also examined two-dimensional arrays of cavities and showed how the presence of neighbouring cavities affected the collapse of cavities placed off the axis of rotational symmetry of the array, amending the direction of the jet away from the centreline of the array. Such geometries are not amenable to simulation with the current axisymmetric code.

Figure 5.

From Dear & Field [26], which uses Schlieren high-speed photography (with an interframe time of 4.25 μs) to show the collapse of three air discs (each of diameter 3 mm) cut into a 200×200×3 mm3 slab of gelatin (sandwiched between two glass blocks, the photograph being taken through the upper glass block). See text for details.

In addition to the separation distance L=0.09 mm, three other cases were studied, i.e. L=0.085, 0.10, 0.20 mm [41]. Figure 6a,b shows the bubble wall position–time histories for the upstream and downstream walls of both bubbles, for these four different initial separation distances. In figure 6a,b, B1 and B2 denote bubble 1 and bubble 2, respectively, and UW and DW denote the bubble upstream and downstream walls, respectively. Note that in figure 6a, the only bubble wall position–time history for bubble 1 that is plotted is the data from the case where L=0.085 mm. This allows the moment of jet impact (where UW and DW meet, labelled X in figure 6b, which is a magnified section of figure 6a) to be placed in perspective, before the other B1 data (for the other values of L) are shown in figure 6c. The high-speed jet impact in bubble 1 occurs earlier as the separation distance L increases. For L=0.085 mm, the jet impacts occurs at t=0.31 μs at a distance x=0.54 mm, where x is the bubble wall position from the left boundary of the computational domain, and for L=0.20 mm, the jet impact occurs at t=0.29 μs at a distance x=0.53 mm.

Figure 6.

(a) The time histories of bubble wall positions for bubble 1 and bubble 2. (b) Magnified section of frame (a). (c) The time history of the position of the downstream walls (DWs) and upstream walls (UWs) of bubble 1 (measured from the points where these walls cross the axis of rotational symmetry) for various initial separation distances L. (d) Correlation between the normalized parameter AL and separation distance L/R0 for bubbles of initial radius 40 μm.

As shown in figure 6b, for L=0.085 mm, the UW of bubble 2 remains nearly stationary until the intense blast wave from the high-speed jet impact of bubble 1 causes a sudden increase in velocity of the UW of bubble 2. This is because bubble 2 is nearly completely shielded by bubble 1 from the incident shock wave. As the separation distance, L, increases, the effects of the incident shock on the collapse of bubble 2 become more profound as bubble 2 is only partially shielded by bubble 1. However, for L=0.09 and L=0.10 mm, the blast wave from jet impact in bubble 1 still causes a significant acceleration of the UW of bubble 2, as marked by a distinct increase in the curve of the UW wall position–time histories (figure 6b), and the jetting collapse of bubble 2 is said to be induced to a large extent by the blast wave from bubble 1.

The case where L=0.20 mm, however, looks different as it appears that the jetting collapse of bubble 1 has an almost negligible effect on the collapse of bubble 2, as shown by figure 6a. For this case, assuming that the incident lithotripter shock moves with the speed of sound in water (approx. 1480 m s−1; indeed direct measurement from figure 6a gives a shock speed of 1481.5 m s−1), the lithotripter shock will reach the UW of bubble 2 at about t=0.14 μs (i.e. obtained from L/Us, where L=0.20 mm and Us=1480 m s−1), and it completely passes by bubble 2 at about t=0.19 μs (i.e. from (L+2R0)/Us, where R0=0.04 mm). After the interaction of bubble 2 with the incident lithotripter shock, the UW of bubble 2 will start to move slightly. Figure 6a demonstrates that the UW of bubble 2 only starts to accelerate at about t=0.25 μs and the jet impact in bubble 2 occurs at about t=0.43 μs. The question to answer is whether the jetting collapse of bubble 2 is induced by the blast wave from the high-speed jet impact in bubble 1, by the incident lithotripter shock wave or by a combination of the two. Given that the distance of the initial location of the UW of bubble 2 (x=0.66 mm) from the location of the jet impact of bubble 1 (x=0.53 mm) is 0.13 mm (i.e. 0.66–0.53 mm), the blast wave from the jet impact in bubble 1 will reach the initial location of the UW of bubble 2 within about 0.09 μs (i.e. 0.13 mm/1480 m s−1). In other words, the blast wave from the jetting collapse of bubble 1 reaches the initial location of the UW of bubble 2 at about t=0.38 μs. By this time, the UW of bubble 2 has already moved away from its original location by about 0.02 mm (figure 6a). The blast wave from bubble 1 will probably catch up the UW of bubble 2 at about t=0.41 μs. However, from figure 6a, it appears that for this case of L=0.20 mm, there is no evidence of a sudden increase in the UW curve of bubble 2 associated with interaction of the blast wave from bubble 1 with the jetting of bubble 2, which is different from the behaviours seen in the cases where L=0.085, 0.09 and 0.10 mm. For this case of L/R0=5 (for the 40 μm bubble radius used here), the bubbles are considered to be far apart.

For the particular conditions under consideration here, there is a critical value of L above which bubble 2 behaves like an isolated bubble in terms of the inducement of jet collapse in that the influence of the collapse of bubble 1 is of secondary importance. This is quantified by plotting a dimensionless parameter AL against the normalized separation distance L/R0 (where R0=40 μm is considered here). The parameter AL is given by AL=(L/Us)/(tj2tj1), where L is the bubble initial separation distance, Us is the incident shock wave velocity, and tj2 and tj1 are the time to jet impact for bubble 2 and bubble 1, respectively. From figure 6d, AL approaches the asymptotic value of unity for large values of separation distance L/R0. If L/R0 is greater than 5 (recognizing that these simulations are axisymmetric and so do not cover two-dimensional or three-dimensional arrays of bubbles) as a working hypothesis, it will be assumed that the two bubbles are sufficiently far apart that they do not affect the dynamics of one another for the kind of lithotripter-induced collapses of initially stationary bubbles being discussed here. This information was used for the calculations of far-field pressure wave signatures from a cloud of cavitation bubbles (in the following subsection), in which the shortest distance between two adjacent bubbles was set to L/R0=5.

(b) Far-field acoustic emissions from cavitation cloud

In this section, the Kirchhoff solution for a single bubble in the free field [15] is extended to a multi-bubble problem. A Gaussian normal distribution is used randomly to distribute the bubbles in water, with a high bubble density concentrated around the focal point of the converging incident lithotripter shock wave (figure 7a). The focal size is commonly used to describe the spatial pressure distribution of the acoustic field of a lithotripter.

Figure 7.

(a) The cloud of 1000 bubbles that is generated in a cigar-shaped lithotripter shock wave focal point (contained within a 60×10×10 mm3 box). Random bubble distribution in panel (b) the xz plane and (c) the xy plane.

Indeed, in vitro passive acoustic and luminescence [52] and high-speed photographic observations [53] show that a bubble cloud is generated inside in the focal point of a lithotripter shock wave in a cigar-shaped volume. From the photographic evidence of a bubble cloud produced using an HM3 lithotripter as observed by Sokolov et al. [53], Krimmel et al. [54] estimated the number density of the bubble cloud in water to be of less than 40 bubbles per cm3.

In this study, the focal region of the lithotripter shock wave is approximated to be of a cigar-shaped volume (60 mm long and 10 mm in diameter) [52,55]. The ‘random’ distribution of bubbles is therefore appropriately shaped by spatially distributing the Kirchhoff control surface for an isolated air bubble using a normal distribution about the focal point of the shock wave (figure 7b,c). Hence, each control surface represents a single bubble in the free field. In the study here, 1000 bubbles are randomly generated, and the observer is positioned 500 mm from the bubble cloud centre, where ϕ=45 °.

A number of assumptions have been made for the prediction of the far-field bubble cloud pressure signature. They are as follows. (i) The greatest density of bubbles is near the focal point of the lithotripter shock wave (in keeping with the assumption that the cloud was generated by the previous SWL pulse) and that the distribution is approximated as a Gaussian distribution. (ii) There is no shielding effect between bubbles, as discussed in the previous section. (iii) The collapse of the bubbles is caused by the incident lithotripter shock wave and not the blast wave emanated from the liquid jet impact of neighbouring bubbles. (iv) The pressure peak positive amplitude seen by each bubble is identical at P+=90 MPa, i.e. the strength of the shock wave remains constant as it traverses through the cloud of bubbles. (v) The cloud of bubbles is initiated from cavitation nuclei left over from a preceding lithotripter shock wave (these are fired into the patient at a rate of 1 or 2 per second). The bubbles then undergo a series of expansion and collapse cycles before reaching a stable equilibrium bubble radius of 40 μm. The far-field acoustic wave of the bubble cloud predicted here is for the interaction of these stable bubbles with the subsequent incident lithotripter shock wave. These assumptions are not inherent in the process, and the scheme could readily be used to include a range of initial bubble conditions (initial bubble sizes and wall speeds, collapse by bubble-generated pressure fields, inter-bubble effects) because these can all be included in the separate computational fluid dynamics (CFD) simulations of the type shown in figure 4. However, each separate bubble response would require separate computation that is beyond the resources available. Hence, in this cloud simulation, all the bubbles are assumed to be at rest with an initial radius of 40 μm. (vi) The pressure signature of the bubble cloud at the observer point is given by a linear summation of the pressure signatures of the individual bubbles. (vii) The collapse time of the bubbles is staggered to simulate the finite time taken for the shock to sweep through the cloud.

The far-field pressure signature shown in figure 8 is calculated by taking a linear summation of the pressure signatures of individual bubbles. The high-pressure region of the far-field signal (between 335 and 355 μs) corresponds to the region of high bubble density. Pressure waveforms originating from a single bubble can also be clearly seen in figure 8 around t=320 μs and t=370 μs. Peaks at the extremities have similar amplitudes because of the uniform initial bubble radius, which would not be a feature in reality. Figure 8a shows the prediction for propagation in water (which, as stated earlier, is compromised by the assumption of no nonlinear propagation outside of that emission's control surface). Such compromise is less serious for propagation through tissue, and figure 8b shows the result when the same data have been derated for tissue (necessarily assumed to be homogeneous) in the manner described by Jamaluddin et al. [15]. The use of the same bubble collapse data in figure 8a,b is for comparative purposes, and not meant to imply that bubble nucleation in tissue will be as likely as that which occurs in water (an issue that is discussed further in §4).

Figure 8.

The predicted acoustic (pressure) signature emanated from a bubble cloud consisting of 1000 bubbles located as given by figure 7b,c, assuming propagation in (a) water and (b) human tissue, with an absorption of 0.3 dB cm−1 MHz−1. The observer position was 50 cm from the lithotripter focal point. The power spectral density for each panel is shown as an inset, and derating clearly causes high-frequency attenuation (as expected) in the spectrum in the plot in frame (b) compared with that in frame (a). Although bubble fragmentation is included in the free-Lagrange bubble wall simulations, these blast wave emissions occur prior to significant fragmentation [14], and so while fragmentation can, in principle, affect cavitation noise [56], it does not do so here because the initial number and size of bubbles are stated within the initial conditions.

4. Discussion

The results of §3b provided, in figure 8b, a prediction of the far-field acoustic emissions from a cloud of 1000 cavitation bubbles produced by SWL, where the propagation between the lithotripter focus and the sensor was through homogeneous tissue with attenuation of 0.3 dB cm−1 MHz−1, attenuation which can be placed in the context of the neglect of nonlinear propagation outside the control surface [15]. The inter-bubble spacings that were assumed allowed bubble–bubble interactions to be neglected, based on the findings from §3a. With greater computational resources, this method would directly allow inclusion of such interactions as were simulated in §3a from bubbles that were closer than those in §3b. For specific data for cavitation in living tissue, the bubble population of figure 7b,c was used, which was based on data taken in water with a dual-pulse lithotripter system. It is generally accepted that, for the same sound field, cavitation in tissue is likely to be much less than that in water [57,58], although there is very little specific experimental data on such direct comparisons for living human tissue [59].

Most simulations undertaken to predict or interpret the far-field pressure emissions from cavitation during SWL understandably use the Gilmore model because of the vastly reduced computational load, which also means that it can be run for the hundreds of microseconds required to encapsulate the complete response of a bubble to a lithotripter pulse (i.e. the time history of the multiple expansions and collapses shown in fig. 25 of Turangan et al. [14]). The free-Lagrange approach is too computationally intensive to continue the bubble time series beyond the first few microseconds after the lithotripter shock wave first impacts the bubble. This is sufficient to show the initial collapse (figure 4), and predict the emissions from this (figure 8b), but unlike the Gilmore model, it cannot (with current computational resources) go on to simulate the subsequent expansion of the compressed bubbles to many times their original size, followed by a second violent collapse some hundred microseconds or more later. This pair of violent collapses (first experimentally identified from the timing of the associated cavitation luminescence and correlated to the acoustic emissions from cavitation by Coleman et al. [60]) dominates the far-field acoustic emission from cavitation during SWL, the first collapse of the cloud generating a first ‘burst’ of amplitude m1, and the second collapse generating a second burst of amplitude m2 some time tc seconds later (exact definitions of m1, m2 and tc are given by Leighton et al. [17]). These are clearly shown, recorded from a patient using the sensor developed by this research programme, in figure 9a.

Figure 9.

(a) A typical signal recorded from a patient during the monitoring of lithotripsy treatment in hospital. The label (‘a’) indicates the electromagnetic triggering signal generated when the lithotripter is fired, producing pickup on the device which is used to trigger data acquisition. The signal after the label (‘b’) shows the passive acoustic emission with two noticeable bursts of emissions (‘1’) and (‘2’) (of amplitudes m1 and m2, respectively, separated by time tc). The peak positive and peak negative pressures from this lithotripter setting, measured in degassed water, were 31 and 6.4 MPa, respectively, and the sensor was placed on the patient's skin 10 cm from the lithotripter focus. Data obtained during use of the device discussed in this study to monitor a patient (data provided by Guys & St Thomas' Foundation Trust, London, courtesy of F. Fedele). (b) Pressure–time history of the Gilmore model to simulate the nominal insonification conditions, for air bubbles in water, presented in frame (a). The assumed bubble size distribution (unknown in vivo) is from 20 to 60 μm radius with equal numbers of bubbles in each micrometre-sized size bin between these radius limits. The spatial distribution of the 1000 bubbles in the simulation is as given in figure 7b,c. The assumed pressure profile of the lithotripter shock wave resembles the shape shown in figure 2b, but has peak positive and peak negative pressures of 31 and 6.4 MPa, respectively. The far-field pressures have been derated following the standard manner for an observation point 10 cm from the source (the appropriate form of derating is discussed in §4).

Comparing the predictions of the Gilmore model with in vitro experimental data, Coleman et al. [52,60] were able to relate the interval between these peaks (tc) with key characteristics of the event, such as the initial bubble size and the magnitude of the lithotripter pulse, as coded into the Gilmore model parameters. In the years since its discovery, several laboratories around the world have found ways to exploit this interpretation of the two peak structure by the Gilmore model to characterize responses as a result of lithotripsy [6166]. Indeed, application of the technique to in vivo data [46] identified the initial bubble size to be used in the CFD simulations in this and earlier papers in the series [14,15,21].

While the ability to estimate the far-field signature from the collapse of a cloud of bubbles (which undergo lithotripter-induced jetting to form blast waves) represents a technical advance, the project had reached the limit of the extent to which simulation-driven insights could progress to a technical device. The results of Coleman et al. [60] had produced the expectation that the eventual clinical sensor would record two significant acoustic bursts subsequent to each lithotripter shock wave, and indeed figure 9a shows that it did. With current computational resources, the free-Lagrange method could only simulate the first of the two bursts, and only for a single bubble size. Therefore, neither of the two key parameters in the diagnosis (m2/m1 and tc) can be estimated using the free-Lagrange method, although in future, this should be possible.

The free-Lagrange results enabled the clinical sensor to be designed so that it has suitable sensitivity, bandwidth and noise floor to measure the expected far-field pressures that had been predicted (e.g. through simulations such as figure 8). However, the slowness of the simulations indicated that the expected point of handover to the next stage in the sequence of simulation and experiment towards a clinical device had been reached. The question was then whether to proceed to empirical data from the sensor that the free-Lagrange results had allowed to be designed, or use the Gilmore model to interpret the clinical data in a more directed manner. As indicated in the preceding paragraph, use of the Gilmore model is by far the most popular method for simulating lithotripter-induced bubble collapses, and it could be used to estimate both m2/m1 and tc. However, the question as to whether such estimations would be useful in guiding the sensor design needed to be addressed because the Gilmore model assumes spherical symmetry at all times, and therefore the source of the pressure wave (at least for the first burst) is incorrect when the Gilmore model is used. Rather than attribute the source to a blast wave, as the free-Lagrange results of this study showed should be done, the Gilmore model launches a pressure pulse from the bubble wall because of the boundary condition across that wall, which generates high pressures in the liquid as the bubble collapses and compresses the gas (with spherical symmetry) within the bubble. The timing of the collapses (giving tc) is likely to be accurate, but the value of m2/m1 could only be valid by coincidence because the mechanism for generation of m1 at least is incorrect in the Gilmore model.

Furthermore, the acoustic radiation pattern that the Gilmore model assumes is spherically symmetric, and only rarely is this adapted to the environment by the introduction of reflecting boundaries, or the imposition of a finite lithotripter shock travel time across a spatially distributed cloud of bubble nuclei (as was done for figure 9b). In figure 9b, the Gilmore model is used to simulate, as near as possible, the conditions of figure 9a. This Gilmore simulation includes a range of initial bubble sizes from 20 to 60 μm radius with equal numbers of bubbles in each micrometre-sized size bin between these radius limits. It derates the signals emitted by the bubbles, but does not derate the form of the incident SWL pulse (so that peak positive and peak negative pressures are 31 and 6.4 MPa, respectively). The latter will lead to an overestimate of the pulse amplitude at the lithotripter focus in the patient, which is very difficult to quantify for the tissue path in question [67]. However, it is preferable to applying the standard derating procedure along the propagation path of the incident lithotripter pulse because during that path, intense focusing occurs and the propagation is highly nonlinear.

As expected, although the prediction of tc is close to that measured, the ratio m2/m1 gives very poor agreement, being in fact greater than 1 (in common with most Gilmore simulations in the literature for lithotripter-induced cavity collapses; Church [40] and Leighton et al. [17]). This is an important point. In introducing figure 9b, it was noted that the only source of asymmetry in the radiated field came from the finite time for the passage of the lithotripter shock wave across a spatially distributed bubble cloud (and even then it lacks bubble–bubble interactions that can affect emissions in clouds; [13,68]). In the free-Lagrange simulations, the ability to introduce a nearby boundary with real material properties was demonstrated [14], but even this cannot mimic the acoustical effects of the stone because the far-field pressures will also include emissions generated by the finite size of the stone (internal reflections—including resonances—and surface waves [17,69]). The passive sensor was designed to detect these in m1 so that they could be included in the diagnosis. Its high-pass filter has a low-frequency cut-off set at 297 kHz [17] to allow it to pick up stone scattering in m1, and to allow m1 to be linked to respiration. None of these features are in the Gilmore simulations. Fedele [69] pointed out the contradiction that if one simulates the far-field emission of only a single bubble in an infinite body of liquid using the Gilmore method, the amplitude of the second burst will be greater than that of the first, which contrasts with the data of the sensor in clinic. Currently, we have data on the emissions of about half a million far-field pressure signatures from patients with kidney stones, and although for some passive emission traces, the amplitude of the second burst (m2) is sometimes close to that of the first burst (m1) (an allowable result if the shock was not effective at contributing to stone fracture, see later), we have only rarely seen m2 exceed m1 (and then only when the signal-to-noise ratio is low). As the vast majority of Gilmore simulations lack stones and directionality and more than one bubble, and the majority of in vivo data includes them, it is important that if the two are compared directly, there is a clear statement of the experimental conditions (the bandwidth of the sensor and its position and orientation with respect to the lithotripter shock, the inclusion or not of a stone and/or patient, etc.).

The final problem with simulations (free-Lagrange, Gilmore and others) of in vivo sound fields is one of the specific account that needs to be taken of absorption by tissue. So far in this study, the standard derating technique has been applied, but even during linear propagation, this cannot provide more than a conservative estimate of the field amplitude. The value 0.3 dB cm−1 MHz−1 was chosen by the Food and Drug Administration to apply to diagnostic ultrasound equipment as part of the procedure to authorize them for sale. It does not specifically represent most of the soft tissue in the propagation path used here, and assumes a linear frequency dependence within the bandwidth of a typical diagnostic device: blast waves close to the source will contain energies outside of that band, and the nonlinear propagation that occurs during SWL and blast wave propagation are not taken into account. For the purposes of this study, standard derating techniques to allow simulations to account for tissue absorption have reached the limit of their usefulness.

Having therefore assessed the limitations of continuing with a simulation-led route through the project, the decision was made to use these simulations to inform the design (sensitivity, signal-to-noise requirements, bandwidth, etc.) of the sensor, and to outline the expected key parameters (limits on m2/m1 and tc), but not to drive the diagnostic interpretations of the waveform (i.e. the key threshold values of m2/m1 and tc). In such a complex system as the nonlinear response of living human tissue to such extreme pressure fields, putative empirical criteria were drawn to make use of the following three simple deductions: (i) poor targeting would cause a low value of m1 (because the stone would not scatter and ring from the incident lithotripter shock if said shock missed it), (ii) weak cavitation would cause a low value of m2, and (iii) weak cavitation would cause a small value for tc (points (ii) and (iii) follow from the comparison by Coleman et al. [52,60] of experimental data with the predictions of the Gilmore method). Further details are given by Leighton et al. [17].

The parameter values were set by comparing the output of the device with the judgements made by consultant urologists, as described by Leighton et al. [17]. The method was as follows. Immediately following each lithotripter shock wave, the device automatically calculated the values of m2/m1 and tc for the received acoustic emission. Leighton et al. [17] had found that, within the two-dimensional space mapped out by m2/m1 and tc, there exists a restricted set of values (0.4<m2/m1<0.8 and tc>100 μs) such that, if the emission conforms to these, the shock can be designated as ‘effective’ in terms of leading to eventual stone fragmentation. If, during a treatment, greater than or equal to 50 per cent of the shocks had been deemed to be ‘effective’, then the device judged the SWL treatment as a whole to be effective.

Figure 10a compares the automatically computed judgement of the device delivered immediately at the end of therapy, with the ‘gold standard’ assessment made by a consultant urologist at the follow-up appointment some time afterwards. Figure 10b compares the assessment of the clinician who delivered the therapy (using all the equipment available in the approximately £1 M diagnostic suite that comes with the lithotripter) with the same gold standard. The data are for 79 patients, combining the results of the two clinical trials, and show that the device outperforms the clinician, agreeing with the gold standard judgement in 78 out of the 79 patients (the one patient it misjudged had a high body mass index, which caused a low signal-to-noise ratio: the device has been redesigned to warn of such problems). In comparison, using the currently available technology, the in-theatre clinician (the radiographer) provided a treatment score that correctly predicted the outcome of only 61 of the 79 therapies. In particular, the passive acoustic device correctly predicted 18 of the 19 treatments that were successful (i.e. 94.7% sensitivity), while the current technology enabled the in-theatre radiographer to predict only seven of the 19 successful treatments (i.e. 36.8% sensitivity).

Figure 10.

Comparison against the ‘gold standard’ score (shown on the abscissa of both plots) given by the consultant urologist at the follow-up appointment three weeks after treatment (where all data except that of the passive acoustic device is made available to the consultant), of the assessment made immediately after the conclusion of the treatment in the clinic by (a) the device (automatically made), and (b) the in-theatre clinician who conducted the treatment. Filled circles indicate that the ‘gold standard’ assessment indicated a successful treatment, and open circles indicate an ‘unsuccessful’ one. Both clinicians gave an integer or half-integer of between 0–5, where a score of ≥3 indicates a ‘success’ was judged. This quantization means that some points in (b) represent several points overlying, the number being indicated by the integer to the right of the data point. The device gave a continuous score, where a ‘successful’ treatment was allocated to all where ≥50% of the shocks had been deemed to be ‘effective’. Dashed lines demarcate the boundaries between ‘successful’ and ‘unsuccessful’ for each assessment. Data from 79 patients obtained by combining the results of two clinical trials, the detailed methodology of which has been reported by Leighton et al. [17].

Subsequently, Guys & St Thomas' National Health Service (NHS) Foundation Trust, London, purchased a device (the first commercially available clinical SWL treatment monitor: Precision Acoustic Ltd., Dorchester, UK) and is conducting a clinical trial to assess to what extent the percentage of ‘effective’ shocks judged by the device after the first 500 shocks reflects that given at the end of treatment, on which the overall success/failure is assessed [70]. If the preliminary good correlation (figure 11) is maintained as further patients are tested, the device will not just be able to assess in real time the effectiveness of SWL treatments, but enable diagnostic tests to be undertaken to determine whether it is worth a given patient undergoing the full SWL. By determining within the first 100 or so shocks (i.e. before the inception of most adverse affects) whether the stone is of a type that will fragment during lithotripsy, or whether the patient needs instead to be sent for ureteroscopic stone removal, the passive sensor will reduce unnecessary and ineffective treatments (and repeat treatments), and enable the patient to be directed to an appropriate alternative therapy. In terms of the management of a healthcare service, such abilities are highly valued because they ‘condense the patient pathway’. The patient pathway describes the route taken by the patient through the healthcare services. Condensing this pathway (e.g. by reducing inaccurate diagnoses, ineffective treatments, waiting times and the number of times that the patient needs to visit the hospital to see different people) was a key aspiration in the 2004 NHS improvement plan [71].

Figure 11.

Correlation between the percentage of shocks judged to be ‘effective’ in the cumulative treatment score after a test dose of the first 500 shocks, and that at the end of the treatment (roughly 3000 shocks). Preliminary results from the clinical trial reported by Fedele et al. [70].

5. Conclusions

The free-Lagrange code was extended to simulate the lithotripter-induced collapses of pairs of bubbles in axisymmetric distributions, and quantify the dependence of the interaction between them on their initial separation. It was then extended to predict the far-field pressure signature from a cloud of bubbles at the focus of a lithotripter. However, concerns remained about the appropriateness of using standard derating for anything more precise than indicative calculations of the effect of tissue absorption. Another limitation arose because current computational resources are not sufficient to extend the spatial extent of the computational domain to the far field, or extend the time of computation to incorporate more than just the initial collapse. The option of solving this issue using Gilmore simulations was examined and rejected because the Gilmore method wrongly attributes the source of the blast wave pressure signature to the dynamic boundary condition at the wall of a spherical bubble. Consequently, this was the point at which the project had to be led by clinical results from the sensor, the overall specification of which (signal-to-noise ratio, sensitivity and bandwidth) could be informed by the free-Lagrange method. Assessment of the stage at which a project moves from relatively inexpensive simulation, to the greatly more expensive clinical trial phase, is key in engineering a successful product of this type.


A.R.J. was supported through a PhD studentship provided by the EPSRC (grant no. GR/N19243; Principal Investigator: T.G.L). The authors are grateful to Dr F. Fedele, Dr A. Coleman and their co-workers for sharing the data presented in figures 10 and 11. Data can be found at

  • Received September 11, 2012.
  • Accepted November 14, 2012.


View Abstract