Simulations of pressure pulse–bubble interaction using boundary element method

Comput. Methods Appl. Mech. Engrg. 195 (2006) 4287–4302 www.elsevier.com/locate/cma

Simulations of pressure pulse–bubble interaction using boundary element method Evert Klaseboer a, Cary Turangan a, Siew Wan Fong a, Tie Gang Liu a, Kin Chew Hung a, Boo Cheong Khoo b,c,* a

Institute of High Performance Computing, 1 Science Park Road, # 01-01 The Capricorn, Singapore Science Park II, Singapore 117528, Singapore b Department of Mechanical Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore 119260, Singapore c Singapore-MIT Alliance, 4 Engineering Drive 3, Singapore 117576, Singapore Received 10 January 2005; received in revised form 31 May 2005; accepted 10 August 2005

Abstract We propose a methodology based on the boundary element method (BEM) to simulate pressure pulse–bubble interaction. The pulse resembles a shock wave and is in the form of a step pulse function incorporated into the Bernoulli equation. Compressibility effects of the water surrounding the bubble are neglected, and the dynamic response of the bubble to the impinging pulse is assumed to be mainly inertia controlled. The interaction induces the formation of a high-speed jet that penetrates the bubble. Results show that bubble shape, collapse time and jet velocity are in good agreement with other numerical models and experiments, and the method is more computationally efficient.  2005 Elsevier B.V. All rights reserved. Keywords: Boundary element method; Shock wave; Bubble collapse; Jet impact; Jet velocity

1. Introduction The study of bubbles [33] can be considered as a special case of multiphase flows. Unlike the symmetric collapse of a gasor vapor-filled bubble, where the bubble surface contracts spherically (as first introduced by Lord Rayleigh [34]), the asymmetric collapse is characterized by the formation of a high-speed jet as the bubble surface evolves. The behaviour of the asymmetric bubble collapse is still a subject of intensive investigation owing to its applications in many areas of practical interest particularly in engineering and science. Kornfeld and Suvorov in 1944 [25] were probably the first to suggest the formation of a high-speed jet in an asymmetric collapsing bubble. However, the phenomenon was only observed experimentally nearly two decades later [30]. This phenomenon is caused by a non-uniform fluid flow around the bubble, either due to the presence of a boundary (a solid or a free surface) in the vicinity of the bubble or due to the bubble interaction with a shock wave. The jetting collapse induced by the shock–bubble interaction (shock-induced collapse) is more violent and destructive than the one caused by the presence of a boundary. This is because the velocity of the accelerating high-speed jet that penetrates the bubble interior can be in the order of 2.6 km/s, (or even higher in certain experiments [8]), and the impact of the tip of the jet on the opposite bubble surface emits an intense blast wave [2,4,7,12,37]. If the bubble is attached to a solid boundary, the jet will directly

*

Corresponding author. Address: Department of Mechanical Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore 119260, Singapore. Tel.: +65 68742889; fax: +65 67791459. E-mail address: [email protected] (B.C. Khoo). 0045-7825/$ - see front matter  2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2005.08.014

4288

E. Klaseboer et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4287–4302

impact onto the boundary surface and may subsequently produce surface indentation. This article deals with shockinduced collapsing bubbles. The phenomena of asymmetric bubble collapse, i.e., high-speed jet formation and blast waves emitted from the jet impact and bubble rebound, are perceived to be the main causes of erosion in hydraulic machinery [26]. Bubble jetting is also found to play an important role in initiation of explosives [2,6,7,20], in underwater explosions [10] and in light emission in acoustic cavitation known as sonoluminiscence [8,9,26,32]. The erosive effects of high-speed jets of collapsing bubbles are also exploited in ultrasonic cleaning and disinfection for medical purposes [19,26]. In extracorporeal shock wave lithotripsy (ESWL), i.e., a modern clinical technique used to fragment stone deposits in patientsÕ organs (such as kidney, bladder or liver parenchyma) through the appliance of shock waves, cavitation bubbles are found to be present near the focal point of the lithotripter. These bubbles originate from the pre-existing micro-gas pockets floating inside the liquid or from gas pockets attached to stone surface or solid micro-particles called ‘‘motes’’ [26]. The gas pockets or cavities interact with lithotripsy shocks that are normally characterized by 9–114 MPa pressure amplitude of 1 ls duration compressive wave followed by 2.8–9.9 MPa pressure amplitude of 6 ls duration tensile wave [11,24]. The interaction causes the cavities to undergo motional activities (formation, expansion and collapse) that are argued to correlate with successful ESWL treatment as well as evidence of collateral damage to kidney and surrounding soft tissues [17]. When the compressive wave of a lithotripsy shock hits the pre-existing micro-bubbles and/or the bubbles produced by previous lithotripsy pulses the bubbles collapse violently. The interaction of this wave with the bubbles induces the formation of high-speed jets whose penetration and subsequent impact onto the opposite bubble surface generates intense blast waves. Together with the lithotripter shock and the shock wave emitted when the bubbles rebound, the high-speed jet and the impact blast wave impose mechanical stresses that have been identified as possible mechanisms of kidney stone fragmentation [11,13,24,45]. The high-speed jet penetration to the surrounding tissue has also been associated with injuries during ESWL treatment; typical injuries include cell lysis, DNA damage in human leukosytes [26], vascular damage, perirenal and intrarenal hematomas [24]. A limited amount of data concerning shock-induced bubble collapse is available in literature. Also, most studies concentrate on a single bubble because it is argued that by understanding the simpler dynamics of a single bubble better, the collapse phenomena of multiple bubbles will be easier to predict. In experiments, the interaction of a planar pressure pulse that may resemble a shock wave with a bubble focuses primarily on the study of a ‘‘disc bubble’’ (2D). Only recently experiments with fully three-dimensional bubbles have appeared in literature [31]. This approach is attractive from an experimental point of view as it offers benefits such as the ease of controlling initial conditions (e.g., cavity size, number, positions and compositions) and the ease of visualization [2,7]. Bourne and Field [7] investigated the collapse of a 6 mm cylindrical air cavity in gelatin/water by a planar shock of high amplitudes (1.88 GPa) in order to study the role of hydrodynamic and adiabatic heating in ignition associated with shock–bubble interaction. For applications in ESWL, Kodama and Takayama [24] used silver-azide pellets as micro-explosives to generate a spherical shock of 10.2 ± 0.5 MPa, which interacts with a 0.8 mm diameter bubble and studied the destructive effect of jet penetration on nearby biological tissue specimens. Numerical simulations offer many advantages in exploring the shock-induced bubble collapse, particularly the ability to capture specific events that signify the collapse phenomena. These events are usually difficult or sometimes impossible to observe experimentally due to limitations in temporal and spatial resolution of experimental diagnostics. Ball et al. [2] implemented a 2D free-Lagrange method (FLM) [1,2,16], to study the collapse of a cylindrical bubble of 3 mm radius, imploded by a step shock of 1.9 GPa amplitude. The simulation used the initial conditions of the experiments by Bourne and Field [7,8]. The results successfully captured many important events, including shock transmission inside the bubble, a ring vortex and the prediction of local heating inside the bubble. Although the 2D bubble collapse is easier to produce from an experimental point of view, and yet reveals several interesting flow physics, it is the axisymmetric or 3D collapse that is most important for practical implications. Jamaluddin [17] modeled the interaction of a spherical bubble with a step shock of various strengths in an axisymmetric geometry using FLM. The simulations are similar to those conducted by Ding and Gracewski [14] who utilized an arbitrary Lagrangian–Eulerian (ALE) method. Besides capturing the collapse events of shock–bubble interaction, Jamaluddin [17] highlighted the limitations of the approach adopted by Ding and Gracewski [14], particularly in capturing the transmitted shock inside the bubble, where the wave adjacent to the bubble interior surface is seen to travel faster than the main shock in the water. The transmitted shock that precedes the main shock front is nonphysical because the water has higher acoustic impedance than that of air. The boundary element method (BEM) has emerged to be the most commonly used numerical method for simulating transient bubble dynamics with asymmetric collapse. Previous studies on BEM include the growth and collapse of a single spherical bubble near a planar rigid boundary [3,42,44], near a curved boundary [36] and a single or multiple bubbles near a free-surface [5,43]. When the tip of the jet impacts onto the opposite bubble wall, two numerical treatments exist: a doublyconnected region [3] or a vortex ring is introduced [44]. This permits the simulation to go beyond the jet impact and a toroidal bubble evolution is possible. Recent applications with BEM also include the indirect boundary element method [38], the

E. Klaseboer et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4287–4302

4289

modeling of a bubble oscillating near an elastic material [23] and the use of jetting collapse as a means of transferring a very small amount of liquid through a hole, thus creating a micro-pump [21]. In order to study a bubble interacting with some form of external disturbance, Blake et al. [4] modeled the response of a single bubble near an oscillating table (low frequency paint mixer). Their model incorporated the influence of the oscillating table into the reference pressure by adding a modified gravitational term associated with the sinusoidal table displacement in the Bernoulli equation. In their work, the reference pressure was a function of time (and not of the spatial co-ordinates). In our work presented here, we model the interaction of an external pressure disturbance with a single spherical bubble in water using BEM. Unlike the simulations carried out by Blake et al. [4], the disturbance is in the form of a step pressure pulse, rather than a sinusoidal wave. The main impetus of this work is therefore to exploit the capabilities of BEM in predicting the interaction of a planar pressure pulse with a spherical bubble and to compare the results with the simulations of other models and experimental observations. The main focus of this work is on the jet formation in such an asymmetric shock-induced collapse. This paper is arranged as follows. In Section 2, a detailed mathematical formulation of the BEM method, including its non-dimensionalisation procedure and some numerical details, is presented. Numerical results obtained using this method, parametric studies of pulse strength, bubble size, the influence of pressure pulse width on the bubble jet tip velocity and their comparisons with those of free-Lagrange method and experiments are given in Section 3. A further analysis and discussions of the physics involved is given in Section 4. Finally, the conclusions are presented in Section 5. 2. Boundary element method and pressure pulse–bubble interaction model 2.1. Boundary element equations and pressure pulse model In this section the mathematical framework for a collapsing bubble in water due to a pressure impulse is developed. Consider a stationary spherical bubble in an infinite medium with no other boundaries. A pressure pulse is moving with a constant velocity Us and a width Ws towards the bubble in the direction of the z-axis as depicted in Fig. 1. The pressure pulse is assumed to have the simplest possible form with amplitude Ps along the width Ws and Pref anywhere else. Pref is assumed to be the standard atmospheric pressure in our simulations. At a fixed point with time the pressure perturbation would look like as in Fig. 2. The duration of the pressure pulse ts is related to Ws as W s ¼ U s ts .

ð1Þ

Fig. 1. Schematic diagram of a pressure pulse with width Ws moving with a velocity Us towards the bubble in the downwards z-direction (vertical). The initial bubble radius is Rm.

Fig. 2. Schematic diagram of a pressure pulse with duration ts and strength Ps as a function of time t. At all other times the pressure equals the reference pressure Pref.

4290

E. Klaseboer et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4287–4302

The pressure pulse at a distance far from the bubble is assumed to maintain the same shape and velocity during the whole simulation, but it can deviate from this shape in the vicinity of the bubble. The content of the bubble is assumed to behave adiabatically. This is so since the typical timescale of the bubble oscillation is very short and therefore the rate of heat transfer across the bubble interface is exceedingly small. In that case, the pressure pg inside the bubble obeys:  c V0 pg ¼ p 0 ; ð2Þ V where p0 and V0 are the initial pressure and volume, respectively, and V is the bubble volume. At the beginning of the simulation, p0 is equal to Pref if a bubble in equilibrium is assumed and surface tension effects are neglected. The ratio of specific heats c is chosen to be 1.4. If a velocity potential U is introduced, the velocity vector u can be written as u ¼ rU.

ð3Þ

In the whole fluid domain the Laplace equation is assumed to be valid: r2 U ¼ 0.

ð4Þ

Furthermore, the unsteady Bernoulli equation can be applied anywhere in the fluid, as well as on the bubble surface: pg ¼ p 1  q

DU 1 2 þ qjuj . Dt 2

ð5Þ

D ¼ oto þ u  r is used in (5). The density of the medium in which the bubble resides is given by q, and the material derivative Dt Surface tension and gravity effects have been neglected for simplicity, but they could easily be incorporated in (5). It is reckoned and seemed reasonable to use two different reference pressures p1 in the Bernoulli Eq. (5); one for the region within the pressure pulse and one outside this region. In our case problem, it is obvious that there is no disturbance imposed on the bubble surface before the front edge of the pressure pulse makes impact on the bubble. As a result the reference pressure should be equal to p1 = Pref, the pressure far upstream of the bubble; the bubble will thus stay at rest. Once the bubble is fully immersed under the pressure pulse, it seems logical that the reference pressure should take on a value reflecting the strength of the shock, thus p1 = Ps. After the pressure pulse has passed, the reference pressure should be set equal to p1 = Pref once more. In other words, the pressure p1 will be set equal to Pref, except when the pulse passes by, where p1 is set equal to Ps, that is 8 > < p1 ¼ P ref ; z < z0 þ tU s  W s ; ð6Þ p1 ¼ P s ; z0 þ tU s  W s < z < z0 þ tU s ; > : p1 ¼ P ref ; z > z0 þ tU s .

The proposed pressure profile can be construed as an extension of the one used by Blake et al. [4]; at any point on the bubble surface the reference pressure p1 is a function of time only (Fig. 2). The Laplace equation is an elliptic equation and the potential anywhere in the fluid domain can always be computed provided that the potential is given on the boundary of the bubble. The boundary element formulation, which relates U and oU on the surface of the bubble, S, can be written on as Z Z oGðy; xÞ oUðyÞ dS ¼ dS. ð7Þ cðxÞUðxÞ þ UðyÞ Gðy; xÞ on on S S Here, ono ¼ n  r is the normal derivative on the bubble boundary with n directed out of the fluid. Eq. (7) is the boundary element equation where x is a fixed point and y is the integration variable, both are located on S. The solid angle at location x is represented by c(x) and G is the Green function or kernel function defined in a three-dimensional (or axisymmetrical) domain as Gðy; xÞ ¼

1 . jx  yj

ð8Þ

2.2. Non-dimensionalisation of the equations The non-dimensional form of the equations developed in Section 2.1 is presented here. The initial radius Rm, which is also the maximum bubble radius, is used to scale all lengths. The pressure Pref is used to non-dimensionalise all pressure terms. Time is non-dimensionalised with: rffiffiffiffiffiffiffiffi q ð9Þ tref ¼ Rm P ref

E. Klaseboer et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4287–4302

and the potential, U, with: sffiffiffiffiffiffiffiffi P ref . Uref ¼ Rm q With the help of (2) and (5) then takes on the dimensionless form:  0 c DU0 V0 1 2 0 ¼ p  þ ju0 j . 1 2 Dt0 V0

4291

ð10Þ

ð11Þ

In the above and subsequent equations, variables with a prime (0) refer to dimensionless variables. And (7) becomes: Z Z oGðy0 ; x0 Þ 0 oU0 ðy0 Þ 0 cðx0 ÞU0 ðx0 Þ þ U0 ðy0 Þ Gðy0 ; x0 Þ ð12Þ dS ¼ dS . on on S0 S0 The equation describing the pressure pulse (6) can be written in dimensionless form as 8 0 0 0 0 0 0 > < p1 ¼ 1; z < z0 þ t U s  W s ; 0 0 0 p01 ¼ P s ; z00 þ t0 U s  W s < z0 < z00 þ t0 U 0s ; > : 0 p1 ¼ 1; z0 > z00 þ t0 U 0s ;

ð13Þ

where the following dimensionless parameters are used: P 0s ¼ P s =P ref ; 0 s

W ¼ W s =Rm ; rffiffiffiffiffiffiffiffi q U 0s ¼ U s . P ref

ð14Þ ð15Þ ð16Þ

Eqs. (14)–(16) give the dimensionless parameters for the pressure perturbation. In (16), U 0s is independent of the maximum bubble radius Rm and only depends on Us, Pref and the reference density of the fluid q. Here Pref = 1 · 105 Pa and q = 1000 kg/m3 (under atmospheric conditions). The parameter z00 is chosen as z00 ¼ 1:0, so that t 0 = 0 coincides with the time the pulse first hits the bubble. 2.3. Pressure pulse (shock) velocity Us If the pressure pulse represents a shock wave, the amplitude of the shock wave Ps and its speed Us are not independent. For a given shock pressure, the particle velocity behind the shock and hence the shock velocity can be computed. In water, it is assumed that the Tait equation of state applies. The density behind the shock, or post-shock density, is therefore:  1=n Ps  q ¼ þ1 qRd ; ð17Þ B where Ps is the shock pressure, qRd = 999.96 kg/m3 is the water reference density at zero pressure and B and n are constants with values of 3.31 · 108 and 7, respectively. As described by Ball et al. [2], the post-shock water particle velocity, u*, can be calculated from s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi   q q u ð18Þ u ¼ uu þ ðP s  P u Þ; q qu where the subscript u refers to the pre-shock condition. In our simulations, for an initially quiescent state of water, uu = 0.0 m/s and Pu = Pref = 1 · 105 Pa. It is known that for most solids and liquids the shock velocity Us is near-linearly proportional to the velocity of the material (relative to the undisturbed medium) behind the shock u*, provided that there is no phase change or existence of porosity [16]. The near-linear relation of experimental shock data is given in the form U s ¼ ak þ Ak u ;

ð19Þ

where ak is the local isentropic sound velocity and Ak is the shock density ratio parameter. For water ak = 1480 m/s and the linear fit of Us vs u* for u* < 3.5 km/s gives Ak = 1.815 (based on shock Hugoniot data [29]). The dimensionless pulse speed U 0s as used in the numerical simulations can now be determined using (16) and (19). For a weak pressure disturbance under atmospheric conditions, the wave speed in water approaches the local sound velocity, which is Us = ak  1480 m/s [26], and the dimensionless value of U 0s therefore assumes the value U 0s ¼ 148, regard-

4292

E. Klaseboer et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4287–4302

less of the size of the bubble. Thus the only parameters that can be chosen freely for this problem are the amplitude of the pressure perturbation P 0s and its width W 0s . U 0s is either constant for a pressure perturbation or given by (16) and (19) for a shock wave. 2.4. Numerical details of the BEM model The surface of the bubble is represented with N nodes. The number of nodes used on the bubble surface is N = 50 for all cases simulated in this work. For each node, (12) gives an equation relating this particular node to all the other nodes. An axisymmetrical formulation as described in Wang et al. [39,40] was used. The integrations are performed with a linear representation of the potential and the normal velocity. A system of equations with size N · N is then obtained with influence 0 matrices G and H and vectors U 0 and oU representing the potential and its normal derivative at each node: on G

oU0 ¼ H  U0 . on

ð20Þ

In the time evolution of the bubble dynamics as enforced by (11), the new potential at the next time step can be determined using the potential and velocities at the current time step. The right hand side of the system of Eq. (20) is now known and 0 (20) can be solved for the unknown normal velocity oU at each node using Gauss elimination. The velocity vector u 0 can be on determined with this normal velocity and the potential distribution along the surface of the bubble. The position vector of a node on the bubble surface r 0 can be updated for each time step using the velocity vector u 0 i.e., Dr0 ¼ u0 . Dt0

ð21Þ

The initial condition at the bubble surface which is initially at rest is U0 ¼ 0

at t0 ¼ 0.

ð22Þ

0

The time step dt is not a constant during the simulation, but is calculated using:  0 V 0 0 dt ¼ dt0 . V 00

ð23Þ

Eq. (23) assures that a smaller time step is employed when the bubble becomes smaller. In all our simulations, the value of dt00 ¼ 1  105 has been chosen. This simple time stepping scheme differs from the one used for example in [39], but the results are similar for both schemes. A smoothing scheme that prevents numerical surface instabilities has been applied at every 10 time steps [22]. Variations of this code have been rigorously and extensively tested through simulations of a bubble near a free surface [39,40] the merging of two bubbles [35] and more recently a bubble near a fluid–fluid interface [22] or elastic material [23]. 3. Results of the BEM calculations, comparison with other methods and experiments 3.1. Preliminaries In this section, the BEM simulations of the response of a gas bubble to a step pressure pulse that resembles a shock wave are presented. The problem under consideration involves a single air bubble of various initial spherical size immersed in water under standard atmospheric conditions and imploded by a planar pressure pulse. Nine cases have been simulated: bubbles with initial bubble radii of 0.1, 1.0 and 10.0 mm are collapsed by a step pressure pulse of 0.528 GPa, 1.011 GPa and 2.06 GPa strength. These cases are similar to those examined by Jamaluddin [17] using a free-Lagrange method (FLM) and those of Ding and Gracewski [14] using an arbitrary Lagrangian–Eulerian method (ALE). As both methods are used as benchmarks in our simulations, they warrant a more detailed description, which is presented below. One case example is given special attention in Section 3.2; the case with initial bubble radius of 1.0 mm and pressure pulse strength of 0.528 GPa with a very large width. In addition, other cases of bubbles interacting with various pulse strengths are conducted and the results concerning the jet velocities are presented in Section 3.3. The free-Lagrange method (FLM) inherits the conventional Lagrangian principle whereby the computational meshes convect in a strict Lagrangian fashion as the flow evolves, but allows changing of ‘‘neighbours’’ during the advance of the simulation. This method prevents mesh distortion that hinders conventional Lagrangian methods from simulating flow simulations with high deformations, such as the phenomenon of bubble collapse under investigation. A FLM code (Vucalm) developed by Ball [1] solves the 2D, unsteady, compressible Euler equations on an unstructured Lagrangian frame using Godunov-type solvers [2,15,16,37]. In the FLM simulation of bubble dynamics, the computational domain is initially filled with computational ‘‘particles’’ and the domain is discretised by constructing a Voronoi mesh, such as

E. Klaseboer et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4287–4302

4293

Fig. 3. (a) Typical Voronoi mesh used in the FLM simulation at t 0 = 0. Due to flow symmetry, only half of the domain is represented. The solid line represents the approximate location of the bubble surface. For the shock–bubble interaction simulations, approximately 19 · 103 cells are initially used (t 0 = 0) with about 980 cells located inside the bubble [17]. (b) Typical node distribution used in the BEM simulation at t 0 = 0; 50 nodes are equally distributed along the bubble surface of the bubble. Both methods use an axisymmetric configuration (only r 0 > 0 is given here), the initial bubble radius is R 0 = 1.

the one shown in Fig. 3a for the initial mesh at t 0 = 0. For the simulations concerning the shock–bubble interaction, an initial number of approximately 1.9 · 104 particles is used. For full details of this method, the reader is referred to [1,17,18]. Another method, the arbitrary Lagrangian–Eulerian method (or ALE), was used by Ding and Gracewski [14]. It involves the use of a finite volume method to solve the Euler equations in a two-dimensional geometry. An artificial viscosity, which introduces a dissipative term in the shock front, was adopted to capture the shock. The fluid surrounding the bubble was assumed to be inviscid but compressible and the process is treated as adiabatic except within the shock front. An adaptive mesh generation technique was developed for the numerical computations whereby a weight function depending on the pressure gradient was used. 3.2. A typical case: a 0.528 GPa pressure pulse acting on a bubble of 1.0 mm radius The case shown in this section corresponds to a strong shock wave interacting with a bubble. A bubble with an initial radius of 1.0 mm at t = 0.0 is selected. This example has been specially chosen in order to show the capabilities of the developed numerical scheme with a very strong shock wave (less severe weaker shock waves and pressure pulses have also been simulated; results are not shown here). The reference pressure is Pref = 105 Pa. The calculations were carried out using both BEM and FLM methods with a shock profile as shown in Figs. 1 and 2. The shock expanse is assumed to be very wide for both the BEM and FLM simulations and with a width W 0s set to a value of 1000. The shock propagates through the water from top to bottom and the elapsed time is measured from the first contact between the bubble and the shock front (i.e., t 0 = 0). For a shock strength of Ps = 0.528 GPa, the particle velocity behind the shock reaches a value of u* = 259 m/s (18), and therefore P 0s ¼ 5280 (14) and Us = 1950 m/s (19) or U 0s ¼ 195:0 (16). In Fig. 3b the mesh for the BEM method is shown. Only 50 nodes are equally distributed along the bubble surface. When comparing Fig. 3a and b, it is clear that the much higher computational speed of BEM is a result of significantly much less number of nodes being used. When the pressure pulse strikes the upstream bubble wall, the impact subsequently imparts a large momentum on the gas-water interface, creates a high-pressure region in the water around the impact location and forces the bubble boundary to form an accelerating jet in the direction of the pulse. Owing to flow symmetry due to the axisymmetrical nature of the problem, Jamaluddin [17] presented only half domains. In our simulations complete pictures of the domain are obtained by transposing the left hand side from the right hand image (mirror image). Fig. 4 depicts the BEM results on the bubble evolution as a result of its interaction with the pressure pulse. At all times prior to t 0 = 0, the velocity anywhere on the bubble surface is identically zero. The results are plotted with a time interval of 0.22 ls. In Fig. 5 the corresponding results are shown for the FLM method (with a time interval of 0.2 ls); these bubble shapes were taken directly from Jamaluddin [17]. On the FLM method, the pressure pulse front represents the Mach contour front in Fig. 5, whereas for the BEM method, the horizontal line in Fig. 4 represents the pulse front, which is always maintained as a straight line. Physically, the transmitted shock inside the bubble propagates slower than the pressure pulse front as shown correctly by the FLM

4294

E. Klaseboer et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4287–4302

Fig. 4. BEM simulation of the interaction of a very wide pressure pulse of 0.528 GPa with a bubble of radius 1.0 mm. The dimensionless parameters are: P 0s ¼ 5280, W 0s ¼ 1000 and U 0s ¼ 195. The time interval between each plot is 0.22 ls. The horizontal line represents the shock front. Also indicated are velocity vector plots. Each picture is surrounded by a frame of 4.0 mm · 4.0 mm, and the sequence goes from left to right and top to bottom. The top bubble surface accelerates first and it accelerates to form a jet, which attains a maximum velocity of around 2 km/s upon jet impact.

Fig. 5. FLM simulation of shock–bubble interaction [17], using the same parameters as in Fig. 4; with a very wide pressure pulse with strength 0.528 GPa acting on a bubble with radius 1.0 mm. The dimensionless parameters are: P 0s ¼ 5280, W 0s ¼ 1000 and U 0s ¼ 195. The time interval between each plot is 0.2 ls. Besides the bubble surface, the other line indicates the Mach contour front (inside and outside the bubble). Each picture is surrounded by a frame of 4.0 mm · 4.0 mm, and the sequence goes from left to right and top to bottom.

simulations due to the lower acoustic impedance of air compared to that of water. Due to its simplicity, these features are not reproduced with the BEM method. At time t = 0 (first image of Fig. 4), the pulse has just reached the upper bubble interface. In the following few images, the bubble surface is flattening at the top as the fluid behind the shock responds to the pressure pulse. As the pulse propagates downwards, the upstream bubble wall deforms rapidly to produce an accelerating jet for both methods (Figs. 4 and 5). In the sixth image of Fig. 4, the leading edge of the pressure pulse has completely passed the bubble. The whole bubble is now submerged in a region of very high pressure; the shock width is very large for this particular case. Due to this pressure difference and liquid inertia, the bubble contracts further while the jet continues to approach the opposite bubble wall, as shown in the last images of Figs. 4 and 5. Although the simulations by BEM omitted the wave propagation inside the bubble, very reasonably similar topographical features of the bubble surface with those of FLM is obtained, providing support

E. Klaseboer et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4287–4302

4295

to the use of BEM in simulating transient bubble dynamics such as pressure pulse–bubble interaction problems. The jet impacts at t = 1.6 ls for the FLM method and at t = 1.74 ls for the BEM method. Apart from bubble shapes, velocity vector plots are also shown in Fig. 4. Similar plots are also presented by Jamaluddin [17], and are in good agreement with Fig. 4 (not reproduced here). Some non-physical (small) velocity vectors are observed in the region preceding the shock front. These are due to the incompressible nature of the BEM method. These velocities are small with respect to the jet tip velocity. Remarkably, in spite of its simplicity, the BEM method produces fairly similar results as obtained in the FLM simulation. Fig. 6 shows a typical mesh as used by the FLM method. A high concentration of nodes can be observed in the region where the shock advances. This is in contrast to the BEM method, where still the original 50 nodes are being used for the whole simulation. In order to make a more quantitative comparison between the FLM and BEM methods, the time history of the jet tip velocity ujet (the velocity at the upper part of the bubble) for this particular case is compared to that of [17]. Jamaluddin [17] took the average of the absolute velocity of a few particles at the tip of the jet; in our simulations, the velocity of a single node on the tip of the jet (on the axis of symmetry) is recorded as a function of time. The results are presented in Fig. 7 for comparison. There are limited differences between the two methods and the general trend is similar. The jet tip velocity reaches a maximum value of about ujet,max = 2000 m/s for the BEM method and about 2200 m/s

Fig. 6. Voronoi mesh of the FLM simulation, showing a highly complex mesh structure near the top and inside the upper part of the bubble following its response to the shock. The top bubble surface evolves to form an accelerating high-speed jet.

2500

ujet (m/s)

2000 1500

BEM FLM

1000

ALE

500 0 0

0.5

1 t (µs)

1.5

2

Fig. 7. Plots of the jet tip velocity ujet as a function of time for the BEM, FLM and ALE methods. The shock wave first hits the bubble at t = 0.0 ls for all simulations. Jet impact occurs at t = 1.74 ls in the BEM simulation and around t = 1.6 ls in both the FLM and ALE simulations. The maximum jet tip velocity reaches a value of ujet,max = 2000 m/s and around 2200 m/s for the BEM and FLM/ALE simulations, respectively.

4296

E. Klaseboer et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4287–4302 0.8 0.7

zc (mm)

0.6 0.5 0.4 0.3

BEM

0.2

FLM

0.1 0 0

0.5

1

1.5

2

t (µs)

Fig. 8. The position of the center of the bubble, zc, as a function of time plotted for both BEM and FLM method for the simulations of Figs. 4 and 5.

for the FLM method, just before jet impact. This particular case was also compared to the numerical results of Ding and Gracewski [14], who implemented an arbitrary Lagrangian–Eulerian (ALE) method. The results in terms of bubble shape are very similar to those obtained with the BEM and FLM methods. The results concerning the jet tip velocity with the ALE method, which depict some ÔoscillationsÕ in the initial period, are also indicated in Fig. 7. Although there are observed differences between the three methods, the general trend is similar and clear. The motion of the bubbleÕs centroid, zc, as a function of time is given in Fig. 8 for both the BEM and the FLM simulations. It can be seen that the general trend for both methods is quite similar. Near the jet impact, the centroid has moved 0.63 mm for the FLM simulation, while it has moved 0.70 mm for the BEM simulation. The centroid moves slightly faster in the FLM simulation as compared to the BEM simulation. The reason for this might be that the shock wave reflects from the bubbleÕs surface in the FLM model, thus generating some momentum transfer to the bubble, while it does not happen to be so in the BEM model. 3.3. Comparison of jet tip velocity for several pressure pulse strengths and bubble sizes After having studied the effects of a shock wave with a strength Ps = 0.528 GPa as presented in Section 3.2, the results of a bubble interacting with pressure pulses of different strengths are presented next. Ding and Gracewski [14] conducted a numerical study of the behaviour of a single bubble with initial bubble radii of Rm = 0.1, 1.0 and 10.0 mm and applied step pressure pulses with strengths Ps = 0.528 GPa, 1.011 GPa and 2.06 GPa for all three bubble sizes using the ALE method. Jamaluddin [17] simulated the same nine cases, using the FLM method. All nine cases have been compared to the numerical results obtained with the present BEM method. In order to have a quantitative comparison, the jet tip velocity is plotted as a function of time. Since the dimensionless BEM model is independent of the bubble radius, the results will only depend on the dimensionless shock strength P 0s . It appears that the numerical results of both Ding and Gracewski [14] and Jamaluddin [17] also scale reasonably according to the same condition: if their results are plotted dimensionless, the velocity–time history plots will be identical for the different bubble sizes (as far as can be judged from their respective figures). Since the results are in essence fairly indistinguishable for the different bubble sizes, it is therefore reckoned that only the results for a bubble radius of Rm = 1.0 mm will be shown. Fig. 9 shows the dimensionless jet tip velocity–time history plot for BEM, FLM and ALE. The comparisons show a reasonably good agreement between the BEM and FLM/ALE results. All curves start at t 0 = 0 when the pulse first hits the bubble. The simulations are stopped when the jet impacts on the opposite bubble surface. Jet impact occurs at a smaller value of t 0 for larger Ps (in all three methods). The jet tip velocity u0jet also reaches higher values for an increasing shock strength Ps. This further strengthens the contention that the BEM method is capable of capturing the essence of the associated bubble dynamics due to shock impact comparable to the other methods, even in the case of a very strong shock wave. The BEM seems to predict a (slightly) lower jet velocity in the early stages of the bubble evolution; we will defer discussion on this difference to Section 4. 3.4. Influence of pressure pulse width on jet tip velocity and collapse time In all the simulations that were presented in Sections 3.2 and 3.3, the bubble had already collapsed before the tail of the pressure pulse passed over the bubble. As a consequence, the width of the shock had essentially no influence on the flow dynamics. This is no longer true for smaller (shock) widths. The aim of this section is therefore to investigate the effects of pressure pulse widths on the bubble dynamics, particularly on the jet tip velocity. This could have implications for shock lithotripsy where the pressure waves are delivered in pulses to ÔdestroyÕ the kidney stones. A pressure pulse of 0.528 GPa magnitude is used to implode a bubble with an initial radius of 1.0 mm (similar to the case of Section 3.2). The dimension-

E. Klaseboer et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4287–4302

4297

u'jet

BEM 0.528 GPa 400

BEM 1.011 GPa BEM 2.06 GPa

350

FLM 0.538 GPa

300

FLM 1.011 GPa FLM 2.06 GPa

250

ALE 0.528 GPa

200

ALE 1.011 GPa ALE 2.06 GPa

150 100 50 0 0

0.005

0.01

0.015

0.02

t' Fig. 9. Comparisons of dimensionless jet tip velocity u0jet as a function of dimensionless time t 0 for the BEM simulations. Results are given for three values of the strength parameter: Ps = 0.528 GPa, Ps = 1.011 GPa and Ps = 2.06 GPa. Results for the FLM and ALE are also indicated. The dimensionless results do not seem to depend on the bubble radius (0.1 mm, 1.0 mm and 10.0 mm). The jet tip velocities reach maximum dimensional values of about ujet,max = 2000 m/s for Ps = 0.528 GPa, 2500 m/s for Ps = 1.011 GPa and 3000 m/s for Ps = 2.06 GPa just before jet impact.

less strength is P 0s ¼ 5280. The dimensionless shock velocity is U 0s ¼ 195 (i.e., Us = 1950 m/s) as obtained from (16), (18) and (19), and is assumed to be independent of W 0s . The dimensionless shock widths are W 0s ¼ 0:2, 0.3, 0.5, 0.9 and 2.0, corresponding to a dimensional shock width, Ws, ranging from 0.2 to 2.0 mm (15) or a shock duration ts ranging roughly from 0.1 to 1 ls (1). The results showing the jet tip velocity as a function of time are given in Fig. 10. It is observed that the dimensionless jet tip velocity reaches a maximum value of about u0jet;max ¼ 10 for the smallest shock width of W 0s ¼ 0:2. It attains a value of u0jet;max ¼ 33 for the shock width of W 0s ¼ 0:3. For the width of 0.5 it reaches a value of u0jet;max ¼ 60, jet impact occurs at t 0 = 0.033 and is indicated with an arrow in Fig. 10. The jet tip velocity increases further for W 0s ¼ 0:9 and reaches a maximum value of about u0jet;max ¼ 100. The jet tip velocity seems to remain fairly constant during the last phase of the collapse (beyond t 0 = 0.014). Jet impact occurs at t 0 = 0.022. In the last case with the highest shock width of W 0s ¼ 2:0, the jet tip velocity reaches its highest value of u0jet;max ¼ 200. This means the dimensional jet tip velocity reaches a value of ujet,max = 2000 m/s. Jet impact occurs at t 0 = 0.0174 or t = 1.74 ls. The results of this last case correspond rather closely to the case presented in Section 3.2. Any further increase in shock width will give rise to essentially the same results; this is so because the bubble has already ÔcollapsedÕ before the back of the pulse has passed the bubble. The minimum width at which this occurs is henceforth referred to as critical width or W 0s;cr (to be further discussed in Section 3.5). It may be noted that all the curves in Fig. 10 overlap initially. From the above, it is apparent the shock width has a very profound influence on the collapse time as well as on the jet tip velocity.

0.2

200

0.3 160

0.5 0.9

120

u'jet

2.0

80

40 0 0

0.01

0.02

0.03

0.04

t' Fig. 10. Plots of jet tip velocity ðu0jet Þ–time (t 0 ) histories for five pressure pulse widths W 0s with a pulse strength P 0s ¼ 5280 and velocity U 0s ¼ 195. The pulse strengths are chosen to be W 0s ¼ 0:2, 0.3, 0.5, 0.9 and 2.0. Pulse widths larger than 2.0 give similar results to those of 2.0. The arrows indicate the moment of jet impact for larger values of W 0s .

4298

E. Klaseboer et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4287–4302

3.5. Bubble collapse time for large pressure pulse width The collapse time of a bubble is defined here as the time measured from the moment the shock hits the bubble to the moment the jet impacts on the opposite bubble interface. If the width Ws of the pressure pulse is relatively large, a simple approximation for the collapse time of the bubble can be found. If the bubble collapses immersed in a liquid with pressure Ps, the theoretical collapse time is proportional to rffiffiffiffiffi q . ð24Þ tcoll;th ¼ Rm Ps This equation is closely related to the Rayleigh collapse time [41] for a spherical vapor bubble collapsing from its maximum size Rm to a radius of 0.0 (though the moment of jet impact is chosen as a reference time here which may cause some differences). For example, for the case study of Section 3.2, a bubble with Rm = 1.0 mm and a shock strength Ps = 0.528 GPa collapses in about 1.74 ls. The above equation gives tcoll,th = 1.38 ls. The (slight) difference can easily be explained in this way. In (24) it is assumed that the high pressure Ps is applied at all points on the bubble at the same instance, while in reality it still takes some time for the pulse to pass the bubble and thereby yielding a larger collapse time than that predicted by (24). Also the Rayleigh collapse time is strictly only valid for spherical vapor bubbles, but it is applied here to non-spherical gas bubbles. The so-called critical shock width Ws,cr for which the subsequent results will not change any more, if Ws is still increased further, can be easily calculated. This occurs if the back of the pressure pulse has yet to pass the bubble location, when the bubble has already collapsed. The pressure pulse moves with a velocity Us. From this Ws,cr is defined as (24): rffiffiffiffiffi q W s;cr  U s tcoll;th ¼ U s Rm ð25Þ Ps or, in dimensionless form, with (14)–(16): rffiffiffiffiffiffiffiffi P ref U 0s ffi. W 0s;cr  U 0s ¼ pffiffiffiffi Ps P 0s

ð26Þ

If W 0s > W 0s;cr , it is expected that the dynamics of the bubble conform to the characteristics of W 0s;cr . The variable W 0s;cr is independent of the bubble radius Rm. In reality, (26) seems to overestimate W 0s;cr slightly. For example the pressure pulse with strength P 0s ¼ 5280 and velocity U 0s ¼ 195 as discussed in Section 3.4 and shown in Fig. 10, W 0s;cr would be 2.68 according to (26). However, the results for W 0s > 2:0 are shown to be fairly independent of W 0s already. This is so, since the effect of the translational movement of the bubble is not taken into consideration in (26), and (26) will therefore give an upper limit for W 0s;cr . With (1) the critical duration of the pressure pulse ts,cr can also be easily obtained. 3.6. Comparisons with available experimental results Thus far, only comparisons with other numerical models are presented. It would be most interesting to see if experimental data yield comparable results. Unfortunately, experimental results such as those conducted by Bourne and Field [8] focused only on the shock-induced collapse of two dimensional bubbles, where an amount of air was trapped in air filled cavity in gelatin. Nevertheless, comparisons with their results can still be relevant in terms of trends observed and the magnitude of collapse times and jet velocities. They used two glass plates to clamp the gelatin and created a 2D bubble of 3 mm thickness. The experiments were carried out to observe the details of the bubble collapse and jet formation. According to them: ‘‘the shock wave easily overcame the weak polymeric bonding in the gel, which flows in a manner indistinguishable from water’’. The density of the gel was slightly lower than that of water: 950 kg/m3. Four cases will be considered: two bubbles with a radius of 6 mm each are imploded with a shock wave with strength of 0.26 GPa and 1.88 GPa, and two bubbles with radius of 3 mm each are imploded with 1.88 GPa and 3.49 GPa shock wave. The first case, as considered by Bourne and Field [8] and labeled as Case 1 here, has an initial radius Rm = 6.0 mm and a shock strength of Ps = 0.26 GPa. The experimental collapse time is about 100 ls with a maximum jet tip velocity ujet,max = 100 m/s. In Table 1 the results are tabulated with the other cases. Case 2 is a bubble with a similar radius (Rm = 6.0 mm) but with an increased shock strength of Ps = 1.88 GPa. This resulted in a considerable faster collapse time of about 7 ls and a jet with a much higher velocity, i.e., ujet,max = 1700 m/s. Another bubble with radius Rm = 3.0 mm was considered in Cases 3 and 4. In Case 3 the shock strength is Ps = 1.88 GPa as in Case 2, and in Case 4 the shock pressure is increased to Ps = 3.49 GPa. The collapse time in Cases 3 and 4 is about 3 ls and 2 ls, respectively. The jet tip velocity reaches extremely high values of about ujet,max = 5000 m/s in Case 3. The jet tip velocity in Case 4 could not be retrieved from the experiments due to the limitation in the cameraÕs temporal resolution. The shock width Ws was not recorded in the experiments. As a first estimate the dimensionless shock width W 0s is taken to be 1000, which is much larger than W 0s;cr , as defined in (26). This case is referred to as Case 1a and the parameters used

E. Klaseboer et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4287–4302

4299

Table 1 Experimental results for four cases according to Bourne and Field [8] Case

Rm (mm)

Ps (GPa)

tcoll (ls)

ujet,max (m/s)

1 2 3 4

6.0 6.0 3.0 3.0

0.26 1.88 1.88 3.49

100 7 3 2

100 1700 5000 Unknown

The cases are numbered 1–4. The collapse times and jet velocities were measured (except for Case 4, where the jet tip velocity is unknown).

Table 2 Numerical parameters and results concerning collapse times and jet velocities for the cases in Table 1 Case

P 0s

U 0s

W 0s

tcoll (ls)

ujet,max (m/s)

u00jet

1a 1b 1c 2a 2b 3 4

2600 2600 2600 18,800 18,800 18,800 34,900

175 175 175 269 269 269 332

1000 ð> W 0s;cr Þ 1 0.2 1000 0.5 1000 1000

14 19 80 5.8 7.9 2.9 2.2

1600 800 250 3100 1700 3100 4200

3.1 – – 2.3 – 2.3 2.2

For Case 1, three different W 0s were applied and for Case 2 two values of W 0s are employed.

(P 0s , U 0s and W 0s ) are listed in Table 2. The results concerning the collapse time and maximum jet velocity ujet,max are also listed in Table 2. The results seem to depend considerably on the value of W 0s . For example, for Case 1a the collapse time is tcoll = 14 ls for W 0s ¼ 1000; it is tcoll = 19 ls for W 0s ¼ 1 (Case 1b) and tcoll = 80 ls for W 0s ¼ 0:2 (Case 1c). The last value corresponds best with the experimental value of tcoll = 100 ls. The experimental jet tip velocity was reported to be about ujet,max = 100 m/s. The numerical maximum jet tip velocity just before jet impact was 1600 m/s, 800 m/s and 250 m/s for W 0s ¼ 1000, W 0s ¼ 1 and W 0s ¼ 0:2, respectively. The value for W 0s ¼ 0:2 is then again the best approximation for the observed numerical data, but is about twice as large compared to the observed experimental result. For Case 2, two values of W 0s have been investigated: W 0s ¼ 1000 (Case 2a) and W 0s ¼ 0:5 (Case 2b). The experimental collapse time is tcoll = 7 ls. It is tcoll = 5.8 ls and tcoll = 7.9 ls for the numerical results pertaining to W 0s ¼ 1000 and W 0s ¼ 0:5, respectively. The jet tip velocity is ujet,max = 1700 m/s for the experiments. The maximum jet tip velocity for W 0s ¼ 0:5 is also about 1700 m/s, while the jet tip velocity for W 0s ¼ 1000 is nearly twice this value at ujet,max = 3100 m/s. In this case the value W 0s ¼ 0:5 gives very reasonable results for both the collapse time and the jet tip velocity. What is interesting to note from the numerics for Cases 1 and 2 is that at the best W 0s arising from the closest matching value of tcoll to the experiments, the same can be inferred for ujet,max. The experimental collapse time for Case 3, tcoll  3 ls, corresponds well with the numerical value of tcoll = 2.9 ls. The experimental jet tip velocity ujet = 5000 m/s is higher than the one predicted by the simulation (i.e., 3100 m/s). No jet tip velocity is available for the experiment of Case 4, but the numerical simulation gives about 4200 m/s. The experimental collapse time for Case 4, tcoll  2 ls, corresponds very well with the numerical value of tcoll = 2.2 ls. The observed trend is that the maximum jet tip velocity for Case 4 is larger than that for Case 3 is true both the experiments and numerics. It must be reiterated however, that the results of Bourne and Field [8] were obtained with essentially two-dimensional bubbles whereas the simulations presented here dealt with spherical bubbles. Also the shock width of their experiments is not given and the properties of gelatin and water are different. In view of the above results, it is probable that the shock width increases with increasing shock strength. This can possibly explain the differences and trends observed between the experimental and numerical results. Although the results differ, the order of magnitude of the values of tcoll and ujet,max between experiments and numerics are comparable. 4. Discussion The main advantage of the boundary element method (BEM) is that only the boundaries of the domain, in this case the bubble surface, need to be meshed. This is the reason for the enormous speed up of the BEM simulations compared to those of FLM or ALE. This means that results can be obtained relatively easy in a few second to a few minutes as compared to several hours or even days for the other methods under investigation. There is a speed up of at least several orders of magnitude. The jet tip velocity has been observed to be very high in the above numerical and experimental results. Values in the order of 1000 m/s or even higher are attained for higher values of Ps (Ps  Pref) and values of Ws > Ws,cr. For a jet, associated with the non-equilibrium bubble collapse in the presence of a nearby solid surface, the maximum jet tip velocity is

4300

E. Klaseboer et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4287–4302

generally up to an order of magnitude smaller i.e., 100 m/s [22]. A possible explanation for this large value of the jet tip velocity for the former lies in the fact that the bubble partly collapses in a fluid with a higher reference pressure, which is Ps instead of Pref. To highlight this further, a new dimensionless jet tip velocity u00jet can be introduced as follows: rffiffiffiffiffi q 00 . ð27Þ ujet ¼ ujet Ps The parameter u00jet is independent of the radius of the bubble, Rm. The value of u00jet;max associated with ujet,max for several typical examples is: with Ps = 5.28 MPa and ujet,max = 200 m/s, u00jet;max ¼ 2:8; with Ps = 528 MPa and ujet,max = 2000 m/s (Fig. 7), u00jet;max ¼ 2:8; with Ps = 1010 MPa and ujet,max = 2400 m/s (Fig. 9), u00jet;max ¼ 2:4; and finally the example with Ps = 2060 MPa and ujet,max = 3200 m/s (Fig. 9) gives u00jet;max ¼ 2:2. Thus a value of u00jet;max varying from 2.2 to 2.8 is obtained. This means that it is possible for pressure pulses with width Ws > Ws,cr to predict the collapse time with (24) and give a reasonable value for the maximum jet tip velocity as follows: sffiffiffiffiffi Ps ; ð28Þ ujet;max  c q where the constant c assumes a value of about 2–3. In most of our numerical simulations, ujet,max occurs just before jet impact. The above approach leads us to investigate if it is possible to re-plot the results of Fig. 9 with the new dimensionless variables for velocity (27) and time given as sffiffiffiffiffi t Ps 00 t ¼ . ð29Þ Rm q In Fig. 11, the results of the simulations for shock strengths of 0.528 GPa, 1.011 GPa and 2.06 GPa are re-plotted with the new dimensionless variables. The results are independent of the bubble diameter (0.1 mm, 1.0 mm and 10.0 mm are used in the FLM and ALE simulations) as discussed earlier. Indeed, the various curves pertaining to the BEM calculations depict a reasonable overlap under the new variables. Even for the FLM and ALE calculations, there is a broad concurrence albeit wider scattering of the data. The new dimensionless collapse time is about t00coll ¼ 1:2–1:4. The new dimensionless jet tip velocity (at jet impact) is about u00jet;max ¼ 2:2–3:0. This means that, for strong shock strengths, it is possible not only to give an estimation of the collapse time with (24), but also to estimate the jet tip velocity upon impact with (28). An outstanding issue arising from the above results is how or why this approach using the BEM method can work so reasonably well. Strictly, the theory behind the boundary element method is only valid for incompressible or very weakly compressible flows. Thus, possibly the dynamics of the bubble prior to the jet impact is mainly governed by the effects of water inertia rather than compressibility effects. This would then explain the results of the BEM using the potential theory depicting such striking resemblance with those of other methods such as FLM and ALE where the compressibility effects are fully accounted for. It is noted that all the curves of Fig. 11 seem to pass through the point as indicated by t00 = 0.7 and

3 2.5

u''jet

2

BEM 0.528 GPa BEM 1.011 GPa BEM 2.06 GPa FLM 0.528 GPa FLM 1.011 GPa FLM 2.06 GPa ALE 0.528 GPa ALE 1.011 GPa ALE 2.06 GPa

1.5 1 0.5 0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

t'' Fig. 11. Jet tip velocity u00jet as a function of time t00 from the BEM, FLM and ALE simulations similar to Fig. 9. Plotted are the simulations for shock strengths of 0.528 GPa, 1.011 GPa and 2.06 GPa. The results are independent of the bubble diameter (0.1 mm, 1.0 mm and 10.0 mm are used in the FLM and ALE simulations). The new dimensionless collapse time is t00 = 1.2–1.4. The new maximum dimensionless jet tip velocity (at jet impact) u00jet;max varies from 2.2 to 3.0.

E. Klaseboer et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4287–4302

4301

u00jet ¼ 1:8 (besides t00 = 0.0 and u00jet ¼ 0:0). A possible explanation for this observation is that in this region, the BEM, FLM and ALE methods are perhaps reflecting the underlying physics as one driven primarily by inertial effects and the weak secondary importance of compressibility. For t00 < 0.4 some discrepancies have been observed between the BEM on the one hand and the FLM/ALE results on the other. Both FLM and ALE give higher values for u00jet in this initial stage of the bubble evolution. The difference can be attributed to the effect of shock scattering which is responsible for a higher velocity in this region. The results obtained by the FLM/ALE methods are for compressible flow of a planar shock wave impacting on the gas bubble. Besides a shock wave being transmitted inside the bubble, the incident shock wave generates a shock wave refraction on the bubble surface which is reflected back into the liquid [27]. Such complex shock reflection on the bubble surface definitely affects the bubble earlier motion. These effects are not taken into consideration by the BEM method due to the assumption of flow incompressibility. Therefore, it is not too surprising to note that there are differences between the results of BEM and FLM/ALE in this early stage of the bubble evolution. The transmitted shock within the bubble is relatively weak [28] and likely does not affect greatly the development/evolution of the bubble constituent with the BEM model where such an effect is not accounted for. As a result, the flow can be well modelled as incompressible after the initial stage (say t00 > 0.4) and before the bubble collapse. In the last stage of the bubble evolution (t00 > 0.8), the jet velocity can be very high so that the jet tip velocity reaches Mach numbers of 2 or even 3. Flow compressibility effects become increasingly important and have to be taken into account. This is why there is more scattering of the results observed in the late stage of bubble jetting in Fig. 11. The above BEM model was applied for an initially quiescent bubble with bubble wall velocity equal to zero. A slight modification can be made to the model to include bubbles that are expanding or contracting when the pulse hits the bubble by adjusting the initial potential (results not shown here). In addition, if the bubbles are extremely small, surface tension effects must be taken into account. The current model has not incorporated surface tension, but such effects can be easily included. 5. Conclusions The dynamics of a pressure pulse that may resemble a shock wave and its interaction with a spherical bubble were investigated numerically in this study. Several important assumptions were made to simplify the mathematical model for studying the general problems of shock–bubble interaction. One of the key assumptions is the neglect of compressibility of the water surrounding the bubble. The problem hence became amendable to solutions by the boundary element method, which is computationally much more efficient. In other words, it was assumed that the dynamic response of the bubble to an impinging shock wave was mainly inertia controlled. Another assumption in the formulation was related to the effects of internal shock waves in the bubble. These effects were neglected in the current approach, since they were thought to be of secondary importance. The boundary element formulation inherently accounted for these simplifications and when compared with other methods such as the free-Lagrange method (FLM) or the arbitrary Lagrangian–Eulerian formulation (ALE), it showed a phenomenal reduction in terms of computational time and effort. The efficacy and accuracy of the proposed approach were compared with cases reported by Jamaluddin [17] using a FLM method and Ding and Gracewski [14] who used an ALE method. Bubble shapes, collapse times and jet velocities all showed very reasonable agreement for a wide range of parameters such as bubble radii, shock strengths and shock widths. The BEM simulations were also compared with experimental observations from Bourne and Field [8]. A good resemblance with the experiments was obtained, even though they used a two-dimensional bubble rather than the present three-dimensional one. One specific case was studied in more detail relating to a very wide pressure pulse with strength of 0.528 GPa and bubble radius of 1.0 mm. In addition, the influence of pulse width was also observed to have a significant impact on the jet tip velocity–time histories. Some general relationships for the bubble collapse time and jet tip velocity were presented and tested against the numerical and experimental data. The findings of this research could be of interest to practitioners in ESWL applications and shed some light on certain applications in sonoluminescence. The current approach might also be interesting for the simulation of a fully three-dimensional problem, for comparison with 3D FLM (or ALE), yet to be developed. Acknowledgement The authors would like to thank A.R. Jamaluddin for providing the free-Lagrange simulations of shock–bubble interactions for comparisons. References [1] G.J. Ball, A free-Lagrange method for unsteady compressible flow: simulation of a confined cylinder blast wave, Shock Waves 5 (1996) 311–325. [2] G.J. Ball, B.P. Howell, T.G. Leighton, M.J. Schofield, Shock-induced collapse of a cylindrical air cavity in water: a free-Lagrange simulation, Shock Waves 10 (2000) 265–276.

4302

E. Klaseboer et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4287–4302

[3] J.P. Best, The formation of toroidal bubbles upon the collapse of transient cavities, J. Fluid Mech. 251 (1993) 79–103. [4] J.R. Blake, G.S. Keen, R.P. Tong, M. Wilson, Acoustic cavitation: the fluid dynamics of non-spherical bubbles, Philos. Trans. R. Soc. London A 357 (1999) 251–267. [5] J.R. Blake, B.B. Taib, G. Doherty, Transient cavities near boundaries: Part 2. Free surface, J. Fluid Mech. 191 (1987) 197–212. [6] N.K. Bourne, On the collapse of cavities, Shock Waves 11 (2000) 447–455. [7] N.K. Bourne, J.E. Field, Shock-induced collapse of single cavities in liquids, J. Fluid Mech. 244 (1992) 225–240. [8] N.K. Bourne, J.E. Field, Shock-induced collapse and luminescence by cavities, Phil. Trans. R. Soc. London A 357 (1999) 295–311. [9] M. Brenner, S. Hilgenfeldt, D. Lohse, Single bubble sonoluminescence, Rev. Mod. Phys. 74 (2002) 425–484. [10] R.H. Cole, Underwater Explosions, Princeton University Press, 1948. [11] A.J. Coleman, J.E. Saunders, L.A. Crum, M. Dyson, Acoustic cavitation generated by an extracorporeal shock-wave lithotripter, Ultrasound Med. Biol. 13 (1987) 69–76. [12] J.P. Dear, J.E. Field, High-speed photography of surface geometry effects in liquid/solid impact, J. Appl. Phys. 63 (1988) 1015–1021. [13] M. Delius, Lithotripsy, Ultrasound Med. Biol. 26 (Suppl. 1) (2000) S55–S58. [14] Z. Ding, S.M. Gracewski, The behaviour of a gas cavity impacted by a weak or strong shock wave, J. Fluid Mech. 309 (1996) 183–209. [15] B.P. Howell, G.J. Ball, Damping of mesh-induced errors in free-Lagrange simulations of Richtmyer–Meshkov instability, Shock Waves 10 (2000) 253–264. [16] B.P. Howell, G.J. Ball, A free-Lagrange augmented Godunov method for the simulation of elastic–plastic solids, J. Comput. Phys. 175 (2002) 128– 167. [17] A.R. Jamaluddin, Free-Lagrange simulations of shock–bubble interaction in ESWL, Ph.D. Thesis, School of Engineering Sciences, University of Southampton, UK, 2004. [18] A.R. Jamaluddin, C.K. Turangan, G.J. Ball, T.G. Leighton, Free-Lagrange simulations of the jetting collapse of air bubbles with application to extracorporeal shock wave lithotripsy. Proc. R. Soc. London A, in preparation. [19] K.M. Kalumuck, G.L. Chahine, G.S. Frederick, P.D. Aley, Development of high erosity cavitating and acoustically enhanced water jets for well scale removal, 10th American Waterjet Conference, 14–17 August 1999, Houston, Texas, Paper 61. [20] V.K. Kedrinskii, The role of cavitation effects in the mechanisms of destruction and explosive processes, Shock Waves 7 (1997) 63–76. [21] B.C. Khoo, E. Klaseboer, K.C. Hung, A collapsing bubble-induced micro-pump using the jetting effect, Sensors Actuators A (2005) 152–161. [22] E. Klaseboer, B.C. Khoo, Boundary integral equations as applied to an oscillating bubble near a fluid–fluid interface, Comput. Mech. 33 (2004) 129– 138. [23] E. Klaseboer, B.C. Khoo, An oscillating bubble near an elastic material, J. Appl. Phys. 96 (2004) 5808–5818. [24] T. Kodama, K. Takayama, Dynamic behaviour of bubbles during extracorporeal shock-wave lithotripsy, Ultrasound Med. Biol. 24 (1998) 723–738. [25] M. Kornfeld, L. Suvorov, On the destructive action of cavitation, J. Appl. Phys. 15 (1944) 495–506. [26] T.J. Leighton, The Acoustic Bubble, Academic Press, Cambridge, 1994. [27] T.G. Liu, B.C. Khoo, K.S. Yeo, The simulation of compressible multi-medium flow Part II: Applications to 2D underwater shock refraction, Comput. Fluids 30 (3) (2001) 315–337. [28] T.G. Liu, B.C. Khoo, K.S. Yeo, Ghost fluid method for strong shock impacting on material interface, J. Comput. Phys. 190 (2003) 651–681. [29] S.P. Marsh, LASL Shock Hugoniot Data, University of California Press, USA, 1980. [30] C.F. Naude´, A.T. Ellis, On the mechanism of cavitation damage by nonhemispherical cavities collapsing in contact with a solid boundary, Trans. ASME D, J. Basic Engrg. 83 (1961) 648–656. [31] C.D. Ohl, R. Ikink, Shock-wave-induced jetting of micron-sized bubbles, Phys. Rev. Lett. 90 (2003) 214502. [32] C.D. Ohl, T. Kurz, R. Geisler, O. Lindau, W. Lauterborn, Bubble dynamics, shock waves and sonoluminescence, Proc. R. Soc. London A 357 (1999) 269–294. [33] A. Prosperetti, Bubbles, Phys. Fluids 16 (2004) 1852–1865. [34] Lord Rayleigh, On the pressure developed in a liquid during the collapse of a spherical cavity, Philos. Mag. 34 (1917) 94–98. [35] S. Rungsiyaphornrat, E. Klaseboer, B.C. Khoo, K.S. Yeo, The merging of two gaseous bubbles with an application to underwater explosions, Comput. Fluids 32 (2003) 1049–1074. [36] Y. Tomita, P.B. Robinson, R.P. Tong, J.R. Blake, Growth and collapse of cavitation bubbles near a curved rigid boundary, J. Fluid Mech. 466 (2002) 259–283. [37] C.K. Turangan, Free-Lagrange simulations of cavitation bubble collapse, Ph.D. Thesis, School of Engineering Sciences, University of Southampton, UK, 2003. [38] C. Wang, B.C. Khoo, An indirect boundary element method for three-dimensional explosion problem, J. Comput. Phys. 194 (2004) 451–480. [39] Q.X. Wang, K.S. Yeo, B.C. Khoo, K.Y. Lam, Nonlinear interaction between gas bubble and free surface, Comput. Fluids 25 (1996) 607–628. [40] Q.X. Wang, K.S. Yeo, B.C. Khoo, K.Y. Lam, Strong interaction between a buoyancy bubble and a free surface, Theor. Comput. Fluid Dyn. 8 (1996) 73–88. [41] Z. Zapryanov, S. Tabakova, Dynamics of Bubbles, Drops and Rigid Particles, Kluwer Academic Publishers, 1999, p. 360. [42] S. Zhang, J.H. Duncan, G.L. Chahine, The final stage of the collapse of cavitation near rigid wall, J. Fluid Mech. 257 (1993) 147–181. [43] Y.L. Zhang, K.S. Yeo, B.C. Khoo, W.K. Chong, Three-dimensional computation of bubbles near a free-surface, J. Comput. Phys. 146 (1998) 105– 123. [44] Y.L. Zhang, K.S. Yeo, B.C. Khoo, C. Wang, 3D jet impact and toroidal bubbles, J. Comput. Phys. 166 (2001) 336–360. [45] S.L. Zhu, F.H. Cocks, G.M. Preminger, P. Zhong, The role of stress waves and cavitation in stone comminution in shock wave lithotripsy, Ultrasound Med. Biol. 28 (2002) 661–671.