TY - JOUR AU - Cheong Khoo,, Boo AB - Abstract Professor John Blake spent a considerable part of his scientific career on studying bubble dynamics and acoustic cavitation. As Blake was a mathematician, we will be focusing on the theoretical and numerical studies (and much less on experimental results). Rather than repeating what is essentially already known, we will try to present the results from a different perspective as much as possible. This review will also be of interest for readers who wish to know more about the boundary element method in general, which is a method often used by Blake and his colleagues to simulate bubbles. We will, however, not limit the discussion to bubble dynamics but try to give a broad discussion on recent advances and improvements to this method, especially for potential problems (Laplace) and wave equations (Helmholtz). Based on examples from Blake’s work, we will guide the reader and show some of the mysteries of bubble dynamics, such as why jets form in collapsing bubbles near rigid surfaces. Where appropriate, we will illustrate the concepts with examples drawn from numerical simulations and experiments. 1. Introduction John Robert Blake was one of the pioneering scientists to explore numerically the rich phenomena associated with oscillating bubbles (in which the bubble is interpreted as a volume filled with either gas or vapour in a medium of much higher density, such as water). One of the most interesting phenomena is the occurrence of high-speed re-entrant jets, which can form when a bubble collapses near an object (such as a solid plate, a free surface, another bubble etc.). These jets are believed to be, at least partly, responsible for the damage occurring due to cavitation on propeller blades. Blake already started this research subject when he was still in Australia (Gibson & Blake, 1982; Blake & Gibson, 1987). Blake’s most cited article is his review with Gibson from 1987 on cavitation bubbles near boundaries (Blake & Gibson, 1987), with more than 650 citations according to Google Scholar (July 2018). As mentioned in that review (Blake & Gibson, 1987), no theoretical solution exists for this problem and we must employ numerical techniques, among which the boundary element method (BEM) turns out to be the most suitable candidate. Readers interested in the subject of bubble dynamics are also referred to the book of Brennen (1995). The work that Blake performed on ciliary propulsion (in Stokes flow) is not reviewed here. Other groups that were/are working on the numerical simulation of bubble dynamics with boundary elements were Harris (1992), see e.g. Zhang et al., (1993) and Ma et al., (2018), Sato et al. (1994) and Wang et al. (1996), to mention a few. It is perhaps time, after many years of research in this area, to take a step back and see what can be learned from all this research, both in terms of the physics concerning bubble dynamics but also in terms of its main numerical tool, the BEM. This is the goal of the current work. Usually, the BEM is considered a difficult method to implement; Becker (1992) even refers to it as ‘A mathematical monster that leaps out of every page’. We will try to convince the reader that this does not need to be the case, by showing some recent enhancements to this method. In particular, we will challenge the widespread conventional belief that singularities are unavoidable. The outline of the article is as follows: since Blake was a mathematician, in Section 2, the BEM is briefly introduced first (but the equations are kept to a minimum in Section 2.1). Then an application of this method, bubble dynamics, with an emphasis on the physics of jets occurring in collapsing bubbles is given in Section 2.2. The physics of jets is discussed in Section 3, together with some examples (both experimental and numerical in Section 3.1) and the concept of the Kelvin impulse (Section 3.2). In Section 3.3, we show that it is possible to predict the formation of a jet using the potential on the surface of the bubble and the Bernoulli equation. We also investigated when the bubble starts to deviate from its spherical shape. Blake’s contribution to bubble dynamics is discussed in chronological order in Section 3.4. Related subjects such as image methods (Blake’s 1971 second most cited article in Stokes flow deals with this) follow in Section 4, with an example of the sound generated by a parabolic piezo electric transducer near an interface. We will also discuss some new developments in the BEM for acoustic simulations. A discussion on the future of bubble dynamics simulations is given in Section 5, followed by conclusions in Section 6. 2. The BEM 2.1 General theory The boundary element (or integral) method is a numerical method used for elliptic partial differential equations. It uses the fact that for elliptic equations (such as the Laplace equation, the Helmholtz equation and the equations describing Stokes flow or even electromagnetic scattering, Klaseboer et al., 2017), the variables on the boundaries determines the solution in the whole domain. The BEM exploits this fact and only requires nodes and elements on the boundaries of the problem. With meshing only the boundary, the dimension of the simulation is reduced by one, and thus BEM is efficient in terms of computational speed and memory requirements. It is noted that there are limitations to BEM flow simulations. Typically, flow irrotationality, incompressibility and inviscidity are assumed. Also heat transfer between the bubble and the fluid is ignored because of the small time scale of the bubble dynamics in comparison with the heat transfer rate. However, BEM is capable of capturing most of the bubble dynamics. A historical overview of the development of the BEM can be found in Cheng & Cheng (2005). Although BEMs were developed from the early 1960s onwards (Jaswon, 1963; Symm, 1963), the 1977 review of Plesset & Prosperetti (1977) on bubble dynamics and cavitation does not mention the use of this method for oscillating bubbles yet, even though other numerical methods trying to simulate re-entrant jets started to emerge mainly based on expansions in terms of spherical harmonics (see e.g. their Fig. 11 or Hsieh, 1972). By the time of Blake & Gibson (1987), BEM had ‘become the principal numerical technique used because of the ease with which it can follow the contortions of the bubble shape’. Let us first focus on potential theory where a potential |$\phi =\phi (\boldsymbol{x})$|⁠, with |$\boldsymbol{x}$| the position vector, satisfies the Laplace equation as $$\begin{align} \nabla^2 \phi = 0. \end{align}$$ (1) Viscosity, compressibility and heat effects have been ignored here. The two following general identities are valid for any smooth function |$\phi $| and |$\psi $|⁠: $$\begin{eqnarray} \phi \nabla^2 \psi = - \nabla \phi \cdot \nabla \psi + \nabla \cdot (\phi \nabla \psi)\ \end{eqnarray}$$ (2) $$\begin{eqnarray} \psi \nabla^2 \phi = - \nabla \psi \cdot \nabla \phi + \nabla \cdot (\psi \nabla \phi). \end{eqnarray}$$ (3) If we subtract (3) from (2), then the first terms on the right-hand side cancel each other out. If we integrate the resulting expression over the volume |$V$| in which they are valid and use the divergence theorem |$\int _V \nabla \cdot \boldsymbol{A} dV= \int _S \boldsymbol{A} \cdot \boldsymbol{n} dS$|⁠, with |$\boldsymbol{n}$| the normal vector pointing out of the volume |$V$| at the surface |$S$|⁠, then we obtain $$\begin{align} \int_V \phi \nabla^2 \psi \textrm{d}V -\int_V \psi \nabla^2 \phi \textrm{d}V = \int_S \phi \nabla \psi \cdot \boldsymbol{n} \textrm{d}S - \int_S \psi \nabla \phi \cdot \boldsymbol{n} \textrm{d}S.\end{align}$$ (4) The next step is to define |$\phi $| as the function we wish to evaluate and to choose |$\psi =G=1/ |\boldsymbol{x}-\boldsymbol{x_0}|$| as the so-called 3D free space Green’s function. Here |$\boldsymbol{x}$| is situated on the surface and |$\boldsymbol{x_0}$| can be situated in the fluid domain or on the surface. If the point |$\boldsymbol{x_0}$| is not situated on the surface, then the first term on the left-hand side of (4) can be evaluated as $$\begin{align} \int_V \phi \nabla^2 G \textrm{d}V = \phi (\boldsymbol{x_0})\int_V \nabla^2 G \textrm{d}V = -4 \pi \phi (\boldsymbol{x_0}). \end{align}$$ (5) The last identity can be proven by taking a small sphere surrounding |$\boldsymbol{x_0}$| and letting the sphere radius approach zero. If |$\boldsymbol{x_0}$| is situated on the surface, then the term |$4 \pi $| must be replaced by the solid angle |$c(\boldsymbol{x_0})$| at the location |$\boldsymbol{x_0}$|⁠. With (1), the second term in (4) is zero. Also, in boundary element terminology, it is customary to replace the normal derivatives such as |$\nabla \phi \cdot \boldsymbol{n}$| by |$\partial \phi / \partial n$|⁠. Then we finally get $$\begin{align} c(\boldsymbol{x_0})\phi (\boldsymbol{x_0}) + \int_S \phi \frac{\partial G} {\partial n} \textrm{d}S = \int_S G \frac{\partial \phi} {\partial n} \textrm{d}S.\end{align}$$ (6) This is now referred to as the boundary integral method. See any book on BEMs (e.g. Becker, 1992) for a more detailed derivation. Usually, |$\phi $| (or |$\partial \phi / \partial n$|⁠) is given and |$\partial \phi / \partial n$| (or |$\phi $|⁠) needs to be calculated on the surface. This can be done numerically by dividing the surface into elements with |$N$| nodes and solving an |$N \times N$| matrix system. Note that even |$\partial G/ \partial n$| is singular as |$1/|\boldsymbol{x}-\boldsymbol{x_0}|$|⁠, due to the fact that |$(\boldsymbol{x}-\boldsymbol{x_0})$| is perpendicular to |$\boldsymbol{n}$| when |$\boldsymbol{x}$| approaches |$\boldsymbol{x_0}$| (contrary to the popular belief that this is a Cauchy principal value integral). Once |$\phi $| and |$\partial \phi / \partial n$| are known on the surface, the value of |$\phi $| anywhere in the domain can be calculated with (6), by putting |$\boldsymbol{x_0}$| at the desired location and taking |$c=4\pi $|⁠. Thus far we have followed the standard textbook approach. In their review, Blake & Gibson (1987) state that ‘Considerable care needs to be taken with evaluating first-kind Fredholm equations, which have singular kernels..... that require special quadrature formulae’. We will now deviate from the classical approach and show how the singularities in (6) can be eliminated entirely analytically and that the above quote essentially no longer applies. We will do so by noting that (6) is not only valid for |$\phi $| but for any function that satisfies the Laplace equation. If we choose a simple analytical function, say |$\chi $|⁠, which satisfies the Laplace equation, then $$\begin{align} c(\boldsymbol{x_0})\chi (\boldsymbol{x_0}) + \int_S \chi \frac{\partial G} {\partial n} \textrm{d}S = \int_S G \frac{\partial \chi} {\partial n} \textrm{d}S.\end{align}$$ (7) And on top of that, we require that the following limits are satisfied: $$\begin{eqnarray} \lim_{\boldsymbol{x}\to \boldsymbol{x_0}} \chi(\boldsymbol{x})=\phi (\boldsymbol{x_0})\, \, \end{eqnarray}$$ (8) $$\begin{eqnarray} \lim_{\boldsymbol{x}\to \boldsymbol{x_0}} \frac{\partial\chi(\boldsymbol{x})}{\partial n}= \frac{\partial\phi(\boldsymbol{x_0})}{\partial n}. \end{eqnarray}$$ (9) If (7) is subtracted from (6), then a fully desingularized BEM (Klaseboer et al., 2012) will emerge: $$\begin{align} \int_S (\phi -\chi) \frac{\partial G} {\partial n} \textrm{d}S = \int_S G \left(\frac{\partial \phi} {\partial n}-\frac{\partial \chi} {\partial n} \right) \textrm{d}S.\end{align}$$ (10) Note, that the term with the solid angle |$c$| has disappeared. There is thus no need to calculate this quantity anymore. As an example, we can use a combination of a constant and a linear function $$\begin{align} \chi(\boldsymbol{x}) = \phi(\boldsymbol{x_0})+\frac{\partial \phi (\boldsymbol{x_0})}{\partial n} \boldsymbol{n}(\boldsymbol{x_0}) \cdot(\boldsymbol{x}-\boldsymbol{x_0})\ \ \, \, \end{align}$$ (11a) $$\begin{align} \frac{\partial \chi(\boldsymbol{x})}{\partial n} = \boldsymbol{n(\boldsymbol{x})} \cdot \nabla \chi (\boldsymbol{x}) = \frac{\partial \phi (\boldsymbol{x_0})} {\partial n} \boldsymbol{n(\boldsymbol{x_0})}\cdot \boldsymbol{n(\boldsymbol{x})}, \end{align}$$ (11b) but there are many other possible choices. For the particular choice of (11), we will end up with $$\begin{equation} \int_S \Big[\phi(\boldsymbol{x}) -\phi(\boldsymbol{x_0})-\frac{\partial \phi (\boldsymbol{x_0})} {\partial n}(\boldsymbol{x}-\boldsymbol{x_0})\Big] \frac{\partial G} {\partial n} \textrm{d}S = \int_S G \Big[\frac{\partial \phi(\boldsymbol{x})} {\partial n}-\frac{\partial \phi(\boldsymbol{x_0})} {\partial n}\boldsymbol{n(\boldsymbol{x})}\cdot\boldsymbol{n(\boldsymbol{x_0})} \Big] \textrm{d}S.\end{equation}$$ (12) Equation (12) is now the totally desingularized boundary integral method. It can easily be seen that all terms within the ‘[.. ]’ on the left- and right-hand side approach zero when |$\boldsymbol{x}$| approaches |$\boldsymbol{x_0}$| and effectively cancel out the singularities in the Green’s function and its normal derivative. Further proof that (10) is indeed non-singular can be found in Klaseboer et al. (2012). Early attempts to desingularize the BEM equations can also be found in Liu (1999)and Liu & Rudolphi (2000). It should be noted that in the final matrix implementation of (10), each node has its own corresponding (11), since for each node |$\phi (\boldsymbol{x_0})$| and |$\partial \phi (\boldsymbol{x_0})/\partial n$| are different. After discretizing the surface into nodes and elements, standard Gauss quadrature integration methods can be employed (there is no need anymore for special quadrature formulae developed for singular integrals). The implementation of quadratic shape functions and quadratic triangular elements is straightforward (essentially the same as in the finite element method). If the problem is discretized with N nodes, finally an |$N \times N$| matrix system can be constructed, which can be solved for the unknowns (usually either |$\phi $| or |$\partial \phi / \partial n$| is unknown). In the implementation for external problems (such as in bubble dynamics), the use of (11) will result in a contribution of the integrals over a surface at infinity (actually |$S$| should be |$S+S_\infty $|⁠). In the conventional BEM as in (6), this contribution will always be zero. However, for the particular implementation of (11), a term with |$4\pi \phi (\boldsymbol{x_0})$| will appear on the left-hand side of (12). Actually, ‘half desingularized’ forms of (11) in the shape of |$\chi (\boldsymbol{x}) = \phi (\boldsymbol{x_0})$| have appeared earlier in the Laplace BEM literature, (see e.g. Harris (1992), his Eq. 7 or Becker (1992)), it is then often referred to as ‘constant value subtraction’. It can also be found in the Ph.D. thesis by Taib (1985) under Blake. The formulation for an axial symmetric problem is slightly more elaborate, but it can be desingularized in a very similar manner as the 3D case; interested readers are referred to Sun et al. (2014). The desingularization techniques as described above are relatively new and have not been used in the majority of the BEM works concerning bubble dynamics. The classical approach is for the integral with the Green’s function to isolate the singular element and to treat this element with a polar coordinate transform, usually combined with the above mentioned half-desingularization for the other singular integral with the normal derivative of the Green’s function. The integrals on all other elements are calculated with the standard Gaussian quadrature method. Finally, it should be mentioned that the fast multipole method (Liu, 2009) can greatly increase the speed and efficiency of BEMs, not only in potential problems but also for Stokes flows and other areas where the BEM can be applied. 2.2 Specific application to oscillating bubbles In the previous section, the general theory of the BEM was discussed. Since there are a number of issues that are very specific for the simulation of bubble dynamics, which are probably not that relevant for the reader who is more interested in the general BEM framework, we will briefly discuss some of those here. The usual potential flow conditions are assumed, i.e. incompressible, inviscid and irrotational flow. Although the BEM is linear, the boundary conditions are non-linear, since they are based on Bernoulli’s equation. Since the goal of the numerical simulation is to follow the bubble (by following its surface), a Lagrangian description is being used. The velocity potential |$\phi $| (⁠|$\boldsymbol{u }=\nabla \phi $|⁠) is updated with the unsteady Bernoulli equation as $$\begin{align} \rho \frac{D\phi}{Dt}=p_0-p_g+\frac{1}{2} \rho |\boldsymbol{u}|^2, \end{align}$$ (13) where D/Dt denotes the material derivative and |$p_0$| is a reference pressure (usually 1 bar). Thus, the potential on the surface is known for all nodes and |$\partial \phi / \partial n$| is searched for using (12). If the response of a bubble in a sound field is required, then this pressure can be incorporated in |$p_0$| as a function of time (provided that the wavelength of the sound wave is much larger than the size of the bubble). The pressure inside the liquid, just outside the bubble, is assumed to be equal to the pressure in the bubble |$p_g$| (for simplicity surface tension and vapour pressure effects are neglected here). The pressure in the bubble is usually assumed to behave adiabatically as |$p_g=p_{g0} (V_0/V)^\gamma $|⁠, with |$p_{g0}$| the initial pressure inside the bubble (usually a value around 100 bar, essentially the driving force for explosion bubbles), and |$\gamma $| the polytropic constant. Finally, all nodes of the mesh themselves are updated with $$\begin{align} \frac{D\boldsymbol{x_0}}{Dt}=\boldsymbol{u}(\boldsymbol{x_0}). \end{align}$$ (14) The velocity |$\boldsymbol{u}$| can be calculated from the normal component |$\partial \phi / \partial n$| and the tangential component (which is obtained from the potential distribution along the surface). After mesh discretization, the actual integration is very similar to the ones usually performed in finite element methods (shape functions, transformation to a unit triangle with a Jacobian etc.). Since usually a few thousand time steps are needed for a full simulation, some kind of smoothing and node redistribution is required along the way. Special care needs to be taken after the jet has impacted on the opposite bubble wall, since the computational domain has now become multiple connected and the bubble assumes a toroidal (donut) shape. For more details on the numerical implementation, including shape functions, linear elements, bubble mesh generation, Jacobians, time step choice, initial conditions etc., see the work of Wang (1998). 3. Jets in bubbles Since Lord Rayleigh (1917) first studied the collapse of a spherical bubble, many more studies on bubble dynamics have followed. The Rayleigh–Plesset equation (Ohl et al., 2015) describes the dynamical response of a spherical bubble. There are quite a number of varieties of this equation, where surface tension, viscosity, external pressure excitation etc., are being added. From a physics point of view, the collapse of a bubble becomes even more interesting when it is no longer spherical symmetric. This can occur when the bubble is situated near a surface or another bubble. The title of this section ‘Jets in bubbles’ is identical to an article by Pearson et al. (2004a), which focused on the origins of the re-entrant jets and their numerical simulation. With this in mind, this section will revisit this issue. 3.1 Examples of re-entrant jets Perhaps the first clear image of a re-entrant jet was obtained by Crum (1979) (reproduced in greater resolution in Blake et al., 1999 and also in Leighton, 2004), in which a bubble was producing reproducible jets towards a table oscillating at 60 Hz with amplitude around 2–3 mm in a tank with water height around 10 cm in a low-pressure environment of 0.04–0.05 bar (near the vapour pressure of water) in which the liquid consisted of a mixture of glycerol and water (L.A. Crum: personal communication, 24 February 2005). In Fig. 1, an oscillating bubble jetting towards a plate situated on top of the figure can be seen. For a general 3D collapsing bubble, the re-entrant jets are usually very difficult to observe, only from the fluid motion after jet impact it can be deduced that such a jet must have been there. In order to see the jetting and subsequent vortex behaviour more clearly, we have generated the bubble in between two transparent narrow plates (which are placed in front and at the back of the figure). In this way, we can clearly see the ‘rolling’ of the vortices on the upper plate. The combined effect of jet impact and shear from the vortices is most likely responsible for the cleaning effect of ultrasonic baths. The bubble was generated with a 60 V low-voltage spark bubble (the electrodes, which are touching, can be seen in Frame 1). This way of generating bubbles is much more convenient than the use of non-touching electrodes, where very high voltages are needed. Fig. 1. Open in new tabDownload slide The experimental expansion and collapse of an oscillating spark generated bubble (at the location of the touching electrodes in the first frame). In order to see the developing vortex structures better, the bubble was placed in between two acrylic plates (plate distance, 4 mm). The maximum bubble radius is about 4 mm. The experimental details are similar to the ones in Azam et al. (2013). The corresponding video is available online as Video 1. Fig. 1. Open in new tabDownload slide The experimental expansion and collapse of an oscillating spark generated bubble (at the location of the touching electrodes in the first frame). In order to see the developing vortex structures better, the bubble was placed in between two acrylic plates (plate distance, 4 mm). The maximum bubble radius is about 4 mm. The experimental details are similar to the ones in Azam et al. (2013). The corresponding video is available online as Video 1. When discussing an earlier paper by Hsieh (1972), Plesset & Prosperetti (1977) mention that ‘as soon as the property of spherical symmetry is lost, the analysis of the various aspects of bubble dynamics becomes exceedingly complex both from the theoretical and experimental point of view’, while Hsieh (1972) himself says that ‘the free boundary value problem becomes extremely difficult’. A quick comparison of the theoretical framework of Hsieh (1972) and the theory developed in Section 2 reveals why the BEM has become the method of choice when dealing with simulations of non-spherical bubbles. A typical BEM simulation of an oscillating bubble near a very large solid plate in axial symmetric configuration with an image method (see also Section 4.2) to account for the solid infinite surface is carried out next (Fig. 2). This bubble has an initial dimensionless radius of |$0.14851R_m$|⁠, with the maximum bubble radius |$R_m=2.6$| mm and an initial pressure of |$p_{g0} = 100$| bar. The bubble is placed at 2.0 |$R_m$| (the so-called stand-off distance) from the solid plate, which is situated at |$z=0$|⁠. Figure 2 shows the simulation results and depicts the bubble’s expansion phase from 0 to 276 |$\mu $|s. The high pressure in the bubble causes the bubble to expand. It reaches its maximum size at 276 |$\mu $|s. By this time, it has over-expanded due to inertial effects, and the reference pressure outside the bubble is now higher than the gas pressure inside the bubble. This pressure difference causes the bubble to collapse. The presence of the solid plate influences the bubble dynamics. The bubble collapses from 276 to 554 |$\mu $|s. In the collapse phase, the bubble moves slightly towards the solid wall. A jet is formed on the top surface of the bubble at 543 |$\mu $|s. This jet develops and hits the bottom surface of the bubble at 554 |$\mu $|s. Fig. 2. Open in new tabDownload slide The (a) expansion and (b) collapse of an oscillating bubble with its centre initially located at 2.0 |$R_m$| from a solid plate (located at |$z=0$|⁠, |$R_m=2.6$| mm). The time for each bubble shape is indicated besides each bubble shape. (a) shows the bubble expanding from its initial size to its maximum size at 276 |$\mu $|s. (b) shows the collapse of the bubble from 276 |$\mu $|s onwards. The bubble moves slightly towards the wall as a re-entrant jet is formed penetrating the bubble from the top surface of the bubble (543–554 |$\mu $|s) towards the plate. Fig. 2. Open in new tabDownload slide The (a) expansion and (b) collapse of an oscillating bubble with its centre initially located at 2.0 |$R_m$| from a solid plate (located at |$z=0$|⁠, |$R_m=2.6$| mm). The time for each bubble shape is indicated besides each bubble shape. (a) shows the bubble expanding from its initial size to its maximum size at 276 |$\mu $|s. (b) shows the collapse of the bubble from 276 |$\mu $|s onwards. The bubble moves slightly towards the wall as a re-entrant jet is formed penetrating the bubble from the top surface of the bubble (543–554 |$\mu $|s) towards the plate. The maximum jet velocity obtained is 113 m/s. From a dimensionless analysis, one can deduce that this velocity is independent of the actual maximum size of the bubble. Since the time is roughly proportional to |$R_m \sqrt{\rho / p_0}$| (with |$\rho $| the density of the liquid and |$p_0$| the reference pressure) and all lengths scale with |$R_m$|⁠, thus the velocity scales with |$\sqrt{p_0 / \rho }$| and is thus independent of |$R_m$|⁠. Results similar to Fig. 2 have appeared in many articles, not only by Blake, but by many others as well (Zhang et al., 1993; Wang, 1998). The first numerical jet simulation was done by Plesset & Chapman (1971) in 1971, using potential theory, the exact numerical implementation is not explained in the article, but should have been a finite difference approach according to Zhang et al. (1993). Figure 3 shows the evolution of a spark-generated bubble (at the centre of Frame 1, very near to the an elastic surface). In Frame 2, the bubble expands and collapses again in Frame 3. Instead of generating a jet towards the surface (as will happen in the case of a rigid wall), the jet is directed away from this elastic surface in Frame 4. The behaviour of an oscillating bubble near an elastic material can be very complex. Apart from the elasticity of the medium on the other side of the water interface, the density of the medium influences the jet direction as well. In another experiment (Gonzalez Avila et al., 2009), a liquid, hydrofluoroether (HFE), which is 1.5 times denser than water (but similar in viscosity), is used. As seen in Fig. 4, the laser-generated bubble, which is in water (liquid on top), expands and collapses with a jet towards the denser HFE (liquid at the bottom). After which (at Frame 7), a liquid jet of HFE is seen ‘jetting’ into the water medium on top. Fig. 3. Open in new tabDownload slide The experimental expansion and collapse of an oscillating bubble near an elastic material. A jet directed away from the boundary is observed. The maximum bubble radius is about 2 mm. Fig. 3. Open in new tabDownload slide The experimental expansion and collapse of an oscillating bubble near an elastic material. A jet directed away from the boundary is observed. The maximum bubble radius is about 2 mm. Fig. 4. Open in new tabDownload slide High speed photography of a bubble oscillating near a water-HFE interface. The bubble is generated by a 15 mJ Nd:YAG laser. The photo sequence was taken with a DSLR camera (D200 Nikon) with an inter-frame exposure time of 5 ms. The pixel resolution is 6 |$\mu m$|/pixel. Reproduced from https://www.aps.org/units/dfd/pressroom/gallery/2009/ohl.cfm with the permission of AIP Publishing. Fig. 4. Open in new tabDownload slide High speed photography of a bubble oscillating near a water-HFE interface. The bubble is generated by a 15 mJ Nd:YAG laser. The photo sequence was taken with a DSLR camera (D200 Nikon) with an inter-frame exposure time of 5 ms. The pixel resolution is 6 |$\mu m$|/pixel. Reproduced from https://www.aps.org/units/dfd/pressroom/gallery/2009/ohl.cfm with the permission of AIP Publishing. 3.2 The Kelvin impulse, application to bubble dynamics If we try to calculate the force on an oscillating bubble, by integrating the pressure just outside the bubble, then a surprising result is obtained. Since the pressure just outside the bubble surface must be equal to that of the internal contents of the bubble |$p_0$| (neglecting surface tension effects) and this pressure is assumed to be constant throughout the bubble due to the large density difference of the bubble’s contents and the surrounding liquid (usually an adiabatic relationship is assumed), the following expression can be obtained for the force on the bubble. |$\boldsymbol{F}=\int p \boldsymbol{n} dS = p_0\int \boldsymbol{n} dS = \boldsymbol{0}$|⁠. Thus, the force on the bubble is zero Best & Blake (1994). We therefore must come up with another concept to describe the dynamics of the bubble in some respect. Blake used the concept of Kelvin impulse to get around this dilemma. The title of this section is ‘The Kelvin impulse, application to cavitation bubble dynamics’, as one of Blake’s sole author articles from 1988 (Blake, 1988). A conference paper with a very similar title appeared already as early as 1983 (again as a sole author contribution, Blake, 1983). As noted in Blake & Gibson (1987), ‘the concept of Kelvin impulse... is likely to provide valuable indicators as to the physical properties required of boundaries in order to reduce or eliminate cavitation damage’. The Kelvin impulse is a recurring subject in many of Blake’s papers (Gibson & Blake, 1982; Blake, 1988; Blake et al., 1986, 2015; Wilson et al., 2008). The Kelvin impulse |$\boldsymbol{I}$| (a vector quantity) was apparently first introduced by Benjamin & Ellis (1966). $$\begin{align} \boldsymbol{I}=\rho\int_S \phi \boldsymbol{n} \textrm{d}S \end{align}$$ (15) In (15), |$\boldsymbol{n}$| is the normal vector pointing out of the fluid domain (thus into the bubble) and the integral is performed over the whole bubble surface |$S$|⁠. It can easily be seen that for a spherically oscillating bubble far from any surface and without gravity, this quantity is zero, since then the potential on the bubble surface, |$\phi $|⁠, does not depend on the spatial coordinates and |$\int \boldsymbol{n} \textrm{d}S = \boldsymbol{0}$| for any closed surface. The main idea behind the Kelvin impulse vector is that its direction will indicate the direction of the jet, which would be useful in predicting damage in e.g. ship propellers, underwater explosions and in countless biomedical applications. The potential in the region with the jet will be relatively high, resulting in the Kelvin impulse being directed in the direction of the jet. In 2008, the Kelvin impulse concept was applied by Blake and his co-workers to predict the behaviour of bubble clouds (Wilson et al., 2008). In 2015, the concept was revisited (Blake et al., 2015) and was one of the last papers Blake worked on. 3.3 Re-entrant jets: an alternative explanation 3.3.1 Introduction Even though we now know from experiments and numerical simulation that re-entrant jets can form in collapsing bubbles near rigid walls, a more intuitive explanation is still lacking. It is well known that a high-pressure zone occurs behind the jet region. However, this will raise the question why a high-pressure region occurs there in the first place. In the spirit of the mathematical background of Blake, let us try to investigate this problem from a potential point of view. We know that the Laplace equation (1) is elliptic in nature. Therefore, whatever happens in the flow field must be determined by the potential on the surface of the bubble, including the formation of a jet. In the next section a seemingly simple but wrong theory will be presented. Then the effect of the Bernoulli equation (13) will be investigated, which acts as a boundary condition for the Laplace equation, since it determines what the potential will be on the bubble surface |$S$| as time progresses. 3.3.2 A seemingly simple, but wrong explanation For simplicity, let us focus on an oscillating bubble near a flat solid infinite interface. Imagine a collapsing bubble near such an interface (Fig. 5). In the absence of the interface, the bubble would have collapsed with velocity vectors indicated in red colour. Due to the presence of the interface, the bubble ‘sees’ it own reflection (we will discuss more on image methods in Section 4), and this will induce additional (blue) velocity vectors in the original collapsing bubble. Upon inspection, it is seen that the part of the bubble closest to the wall experiences a velocity reduction (by adding the velocity vectors), while the part of the bubble that is situated farthest from the wall experiences a higher velocity than would be the case in the absence of the wall. It thus appears that we have found a physical acceptable explanation for the jetting. Fig. 5. Open in new tabDownload slide Schematic of a collapsing bubble near a flat solid rigid infinite interface (left) and the principle of the image method (right). In the potential flow framework, the rigid wall can be represented by an image bubble. Apart from the velocities induced by itself (red arrows), the collapsing image bubble creates additional velocities (indicated with blue arrows) directed towards the image bubble. Fig. 5. Open in new tabDownload slide Schematic of a collapsing bubble near a flat solid rigid infinite interface (left) and the principle of the image method (right). In the potential flow framework, the rigid wall can be represented by an image bubble. Apart from the velocities induced by itself (red arrows), the collapsing image bubble creates additional velocities (indicated with blue arrows) directed towards the image bubble. However, this simple explanation cannot be the full story. This can easily be seen if we not only consider the collapse phase but also the expansion phase of a bubble. If we start with a small spherical bubble, then during expansion the part farthest from the wall will also experience a greater outward velocity according to the same reasoning as above. Since the process described above is totally symmetric with respect to expansion and collapse (in essence, all vectors in Fig. 5 will have the opposite sign for an expanding bubble), the bubble would end up in its original spherical shape after a full cycle of expansion and collapse. In order to understand what is going on, we will have to invoke one more piece of physics, the Bernoulli equation. 3.3.3 The role of the Bernoulli equation If we take the Bernoulli equation (13) on the surface of the bubble and integrate with respect to time, then we get $$\begin{align} \phi(t) = \phi(t=0) +\int \Bigg[\frac{p_0 - p_g}{\rho} + \frac{1}{2}|\boldsymbol{u}|^2 \Bigg]\textrm{d}t. \end{align}$$ (16) The first part of the integral |$(p_0 - p_g)/\rho $| is the same for any position on the surface of the bubble and thus cannot cause any anisotropies such as a jet (although the internal pressure of the bubble, |$p_g$|⁠, does change with time). Note that we have neglected surface tension effects. Thus, the remaining term with the squared velocity must be responsible for any jetting behaviour (assuming that |$\phi (t=0)$| is the same for the whole bubble surface). Based on the fact that the absolute value of the velocity near the wall is lower than that at the farthest point, we can see that a specific value of the potential, say |$\phi =0$|⁠, will be reached earliest for the point situated at the farthest side. In practical terms, this means that the velocity will have reached zero here and that this part of the bubble surface will start to collapse before the rest of the bubble does. This will ultimately result in the re-entrant jet. Most likely, Blake had this in mind when he made his remarks on high and low impedance, as given in Blake et al. (1986): ‘Collapse occurs from the low-impedance side yielding the instability cum liquid jet on the far side of the bubble’. In order to check the above hypothesis, some BEM calculations were performed for a bubble with initial internal gas pressure |$p_{g0}=500p_0$|⁠, |$\gamma =1.25$| and its initial centre placed at a distance of |$1.5R_m$| from an infinite flat wall (the standoff distance). The normal velocities and potentials in Fig. 6 for a node nearest to the wall, a node placed at the opposite pole (’far’ away from the wall) and a node in between those two nodes, dubbed here ‘Equator’. Note that for general potential problems, always a constant potential can be added without changing the velocity field. However, for the BEM framework, the potential at infinity must necessarily be zero (this can easily be observed from (6), by taking |$\boldsymbol x_0$| to infinity, the Green’s function will automatically guarantee a zero potential at infinity). Fig. 6. Open in new tabDownload slide (a) Normal velocity as a function of time for three nodes, two on each pole (the initially farthest and nearest node to the wall) and one on the equator. The parameters are |$p_{g0}=500p_0$|⁠, |$\gamma =1.25$| and standoff distance 1.5 |$R_m$|⁠. Time is normalized as |$t^{\prime}=t/(R_m \sqrt{\rho /p_0})$| and the normal velocity as |$(\partial \phi /\partial n)/\sqrt{p_0/\rho }$|⁠. (b) Potential as a function of time. Time is normalized as |$t/(R_m \sqrt{\rho /p_0})$| and the potential is normalized as |$\phi /(R_m \sqrt{p_0/\rho })$|⁠. The jet starts to form just before |$t^{\prime}=2.0$|⁠. Fig. 6. Open in new tabDownload slide (a) Normal velocity as a function of time for three nodes, two on each pole (the initially farthest and nearest node to the wall) and one on the equator. The parameters are |$p_{g0}=500p_0$|⁠, |$\gamma =1.25$| and standoff distance 1.5 |$R_m$|⁠. Time is normalized as |$t^{\prime}=t/(R_m \sqrt{\rho /p_0})$| and the normal velocity as |$(\partial \phi /\partial n)/\sqrt{p_0/\rho }$|⁠. (b) Potential as a function of time. Time is normalized as |$t/(R_m \sqrt{\rho /p_0})$| and the potential is normalized as |$\phi /(R_m \sqrt{p_0/\rho })$|⁠. The jet starts to form just before |$t^{\prime}=2.0$|⁠. Very soon after the creation of the bubble at |$t^{\prime}=0.000$|⁠, the minimum dimensionless normal velocity of |$\partial \phi /\partial n = -10.3$| is reached for almost all nodes on the surface simultaneously at |$t^{\prime}=3.9\times 10^{-3}$|⁠. Very shortly after at |$t^{\prime}=1.09\times 10^{-2}$|⁠, the dimensionless potential reaches its minimum value of |$\phi =-1.405$|⁠. The potential |$\phi $| for the location farthest away from the wall changes sign from negative to positive at dimensionless time |$t^{\prime}=1.020$|⁠, while the equator changes sign at |$t^{\prime}=1.041$| and the nearest point at |$t^{\prime}=1.124$| (all times are non-dimensionalized with |$R_m \sqrt{\rho /p_0}$| again and named |$t^{\prime}$|⁠). For the normal velocity |$\partial \phi /\partial n$| (non-dimensionalized with |$\sqrt{p_0/\rho }$|⁠), the farthest node changes sign at |$t^{\prime}=0.999$|⁠, the equator at |$t^{\prime}=1.030$| and the nearest node at |$t^{\prime}=1.378$|⁠. This indicates that the farthest node starts to collapse first eventually resulting in a jet. Jet impact occurs at |$t^{\prime}=2.091$|⁠. The plots for both normal velocity and potential start to deviate at larger times. Very similar plots can be obtained for bubbles at other distances from the wall or other values of the initial internal bubble pressure. A plot very similar to Fig. 7 showing the potential appeared in Zhang et al. (1993). Fig. 7. Open in new tabDownload slide Phase plot, showing the potential versus the normal velocity. The parameters are |$p_{g0}=500p_0$|⁠, |$\gamma =1.25$| and standoff distance 1.5 |$R_m$|⁠. The potential is normalized as |$\phi /(R_m \sqrt{p_0/\rho })$| and the normal velocity as |$(\partial \phi /\partial n)/\sqrt{p_0/\rho }$|⁠. The solution of the spherical Rayleigh Plesset solution (18) is also indicated. The arrow shows the evolution with time. Towards the end of the simulation, a jet is formed and the phase plot for the farthest node (now located inside the jet) reaches a constant |$\partial \phi / \partial n$| value while |$\phi $| still increases. Fig. 7. Open in new tabDownload slide Phase plot, showing the potential versus the normal velocity. The parameters are |$p_{g0}=500p_0$|⁠, |$\gamma =1.25$| and standoff distance 1.5 |$R_m$|⁠. The potential is normalized as |$\phi /(R_m \sqrt{p_0/\rho })$| and the normal velocity as |$(\partial \phi /\partial n)/\sqrt{p_0/\rho }$|⁠. The solution of the spherical Rayleigh Plesset solution (18) is also indicated. The arrow shows the evolution with time. Towards the end of the simulation, a jet is formed and the phase plot for the farthest node (now located inside the jet) reaches a constant |$\partial \phi / \partial n$| value while |$\phi $| still increases. 3.3.4 Comparison with a spherically oscillating bubble We will now investigate to which extent the solution deviates from that of a spherically oscillating bubble with an adiabatically behaving interior. In order to do so, we can use the equation of motion for a spherical symmetrically oscillating bubble (the Rayleigh–Plesset equation), where surface tension and viscosity of the liquid have been ignored: $$\begin{align} \rho \Bigg[\frac{3}{2}\Bigg(\frac{dR}{dt}\Bigg)^2 + R \frac{d^2R}{dt^2} \Bigg]=-p_0 + p_{g0}\Bigg(\frac{R_0}{R}\Bigg)^{3\gamma}. \end{align}$$ (17) This equation cannot be solved analytically if |$R$| is required as a function of time |$t$|⁠. However, the analytical solution of |$dR/dt$| as function of |$R$| can be obtained as (with |$R_m$| the maximum bubble radius and |$R_0$| the initial bubble radius): $$\begin{align} \rho\frac{3}{2}\left(\frac{dR}{dt}\right)^2=p_0 \Bigg[-1+\frac{R_m^3}{R^3} \Bigg]-\frac{p_{g0} R_0^{3\gamma}}{(\gamma-1)R_m^{3\gamma}}\Bigg[\frac{R_m^{3\gamma}}{R^{3\gamma}}-\frac{R_m^3}{R^3}\Bigg], \end{align}$$ (18) as can be verified by calculating |$d^2R/dt^2$| from (18) and back substitution of |$dR/dt$| and |$d^2R/dt^2$| into (17). We have assumed here that at |$R=R_0$|⁠, |$dR/dt=0$| and at |$R=R_m$| again |$dR/dt=0$|⁠. Since the potential for such a spherical bubble is |$\phi =-\frac{R^2}{r}\frac{dR}{dt}$| and thus at |$r=R$|⁠, we get |$\phi =-R\frac{dR}{dt}$| and |$\frac{\partial \phi }{\partial n} = -\frac{dR}{dt}$|⁠. In Fig. 7, the 8-shaped figure shows the solution of the spherical bubble, (18), together with the solution of the jetting bubble of Fig. 6, plotted in terms of |$\phi $| and |$\partial \phi /\partial n$| as a phase plot. The point (⁠|$\phi =0$|⁠, |$\partial \phi /\partial n=0$|⁠) corresponds to both the initial bubble size |$R=R_0$| as well as the maximum bubble size |$R=R_m$|⁠. Again, as in Fig. 6, the solution at both poles and the equator are shown. The phase plot shows that the numerical solution starts off in the lower left quadrant still closely following the Rayleigh–Plesset solution. We can see that the solution of the farthest node starts to deviate strongly from the Rayleigh–Plesset solution in the upper right quadrant and |$\partial \phi /\partial n$| becomes almost constant, while |$\phi $| increases in amplitude. This is why the jet starts to develop. Another way to express the above findings is with a quote from Plesset & Prosperetti (1977): ‘the asymmetry of the flow that the boundary itself introduces is sufficient to give rise to highly distorted bubble shapes’. As we have seen from the examples in Section 3.1, oscillating bubbles with impinging re-entrant jets are an ideal way to generate vortices in a initially quiescent liquid. Finally, for completeness, we should mention that there is still a completely different way of creating jets in bubbles, by exposing these bubbles to a shock wave (Ding & Gracewski, 1996) or an acoustic traveling wave of sufficient amplitude (Calvisi et al., 2007). Very high-speed jets of several km/s can be observed (even when there is no nearby surface), much higher than the ones obtained under ‘normal circumstances’ (for cavitation bubbles or underwater explosions the jet speeds are usually around 100 m/s, see also Section 3.1). The mechanism behind jet formation is believed to be the fact that the side facing the shock wave receives a very high pressure and starts to collapse before the back of the bubble ‘feels’ the shock wave, resulting in an uneven collapse and finally jet formation in the travel direction of the shock wave. Such a shockwave–bubble interaction has been successfully simulated using a BEM framework (Klaseboer et al., 2007). The jet speed is function of the shock strength (Ohl et al., 2015). 3.4 Blake’s contribution to bubble dynamics When bubbles oscillate near surfaces (even multiple), jets can develop as we have seen in the previous sections. These jets can be directed towards or away from the surface, depending on the surface properties (solid surface, free surface, elastic surface, flat or rounded surface, other nearby bubbles etc.). Blake, together with his coworkers extensively studied such bubbles near all kinds of surfaces. In this section, we will give a brief overview of their work. Blake & Gibson (1981) studied the behaviour of a bubble near a free surface. A method based on potential flow but not on the BEM was used. Instead of the now familiar boundary element framework, a singularity distribution was used inside the bubble and on top of the free surface. A relatively simple model to predict the behaviour of bubbles near deformable surfaces appeared in 1982 (Gibson & Blake, 1982). In Blake et al. (1986), oscillating bubbles near rigid flat interfaces were investigated but now with a ‘real’ BEM as we know it today. Particularly interesting are the path-lines that fluid particles situated on the bubble surface follow (their Fig. 3). It turns out that most of these fluid particles eventually end up in the jet region. Also shown are the high pressure regions in the fluid that essentially drive the jet. A follow up (Part 2) was published in 1987 (Blake et al., 1987), now dealing with bubbles near free surfaces. From a physics point of view, it is interesting that a high-pressure region develops in between the bubble and the free surface, causing two jets in the opposite direction; one inside the bubble directed downwards and one on the free surface directed upwards. The effect of buoyancy was also studied (mainly of interest for very large bubbles such as in underwater explosions). In 1989, a purely experimental paper on bubbles near composite surfaces appeared (Shima et al., 1989). Depending on the surface properties, the jets of nearby oscillating bubbles can be directed towards or away from the surface. Two cavitation bubbles near a rigid boundary (in an axial symmetric configuration) were studied in 1993 (Blake et al., 1993) and compared to experimental results produced with a laser. A bit off topic work in the current context, but nevertheless worth mentioning, is the effect of viscous boundary layers in gas bubbles bursting near free surfaces which was investigated by Boulton-Stone & Blake (1993). Returning to the topic of oscillating bubbles, the concept of Kelvin impulse (see Section 3.2) was studied in 1994 (Best & Blake, 1994). A ‘manual’ of how to utilize BEMs for bubble simulations appeared in 1995 (Blake et al., 1995). In it appeared a detailed discussion on the numerical techniques used, as well as a 3D example of jet formation. More 3D examples of collapsing bubbles due to a nearby flat vertical surface and gravity are shown in Blake et al. (1997). As a last example in this article, an oscillating bubble is compared against the Rayleigh–Plesset equation. Due to numerical instabilities, the numerical solution starts to deviate considerably from the theoretical value after around 10 oscillations and develops a ring jet (while it should have stayed spherical). A cluster of 13 collapsing bubbles (essentially a simple model for a bubble cloud) due to an acoustic pulse was simulated in 1999 (Blake et al., 1999). The bubbles situated at the outer parts collapse first before the more central bubbles start their collapse; they are essentially shielded for a while by the outer bubbles. The famous Crum (1979) bubble-jet was also simulated in the same article. In the same year, the role of splashing was discussed in another article (Tong et al., 1999). In 2001, an article followed on bubbles interacting with a free surface (Robinson et al., 2001) followed by another one in 2004 (Pearson et al., 2004b). This kind of bubble often occurs in underwater explosions, in which the bubbles can have radii up to 10 m. Besides a jet inside the bubble, which is now directed away from the free surface, a jet can also be observed on the free surface away from the bubble. Some simulations on a system of two bubbles near a free surface were also studied. Bubbles near vortices arising in hydraulic and naval applications were studied in 2012 (Wilson et al., 2014). The Kelvin impulse concept (see also Section 3.2) was revisited once more in 2015 (Blake et al., 2015). 4. The method of images In potential theory, a powerful tool is the use of a Green’s function that has an image. Suppose that we add an image of the Green’s function then the potential on the ‘mirror surface’ will have a zero normal derivative (this will now be a convenient way to introduce an infinite flat rigid surface, without actually having to mesh an infinite domain). Similarly, the subtraction of an image Green’s function leads to a zero-potential surface, such as an infinitely flat free surface. Blake has considerably contributed to this area by finding the image for a Stokeslet in low Reynolds number flow. 4.1 Blake’s note on the image system for a Stokeslet Blake not only studied bubble dynamics, in which usually potential flow is used (i.e. the Reynolds number of the flow is very high and only inertial effects are considered), but also did quite some work on Stokes flows, in which viscous effects are dominant and inertial effects can be neglected. An example that should be mentioned here (although it is not directly related to bubble dynamics) is the image method for Stokes flow (Blake, 1971) (Stokeslet near a flat solid wall), which appeared in 1971 with the title ‘A note on the image system for a Stokeslet in a no-slip boundary’. Although it is a ‘note’, nevertheless it has become Blake’s second most cited article with around 630 citations according to Google Scholar (July 2018). Unlike for the potential problem, the placement of an image Stokeslet is not sufficient to guarantee that the velocity on the flat surface is zero (not only the normal velocity but also the tangential velocity must be zero). Blake therefore introduced higher-order singularities to ensure that the velocity on the wall is indeed zero (see also Blake & Chwang, 1974). Without going into further details, we will mention here that the Stokes BEM (which gives a relation between the velocity and traction vectors in Stokes flow and uses tensor notation) can also be desingularized as was shown in Klaseboer et al. (2012), essentially according to the same principle as shown in (10). Readers interested to further study the BEM in Stokes flow are referred to Pozrikidis (2002). 4.2 Image methods for the Laplace equation The BEM works best for problems that have a relatively small surface but a large volume surrounding it. In other words, a small bubble surrounded by an ‘infinite ocean’ is a perfect example. Where other numerical methods need to mesh the whole fluid domain and take care of boundary conditions at infinity, the BEM only needs a mesh on the surface. When this bubble is situated near a flat infinite surface, potential theory allows us to retain its simplicity by using the method of images. This eliminates the need to use a mesh on the flat surface and keeps the numerical procedure elegant. For example, for a flat rigid surface, we can use (Blake et al., 2015) $$\begin{align} G=\frac{1}{|\boldsymbol{x}-\boldsymbol{x_0}|} + \frac{1}{|\boldsymbol{x}-\boldsymbol{x_0^I}|}, \end{align}$$ (19) in which |$\boldsymbol{x_0^I}$| represents the image of node |$\boldsymbol{x_0}$| with respect to the flat surface. The first term on the right of (19) is the ‘normal’ free space 3D Green’s function; the second term never becomes singular since |$\boldsymbol{x}$| never becomes equal to |$\boldsymbol{x_0^I}$|⁠. Blake and co-workers used this technique extensively in their bubble dynamics simulations, see e.g. Blake et al. (1995). This method cannot always be used, e.g., if an oscillating bubble is situated very near to a free surface, then there is no other choice than to model this surface as well since the free surface will no longer remain flat. Rather than repeating the examples for oscillating bubbles that can all be found in Blake’s articles, we will show some examples taken from acoustics in the next section. Even though Blake did not work with the acoustic (small perturbation) wave Helmholtz equations as far as we know, he did study ultrasonic cavitation in 2013 (Curtiss et al., 2013) for example. These acoustic simulations explain e.g. what happens with acoustic transducers (and their interactions with bubbles). This will be investigated in the next section. Due care must be taken when implementing two parallel walls with images; this will change the nature of the mathematical problem. For an oscillating bubble (with volume change), the main contribution to the potential will be a source/sink term. The pressure will then no longer decay as 1/r, but increase with ln(r), due to the presence of an infinite set of mirror images. 4.3 Image methods for the Helmholtz equation, examples from acoustics near boundaries A subject closely related to the topic of bubbles is the topic of insonification, i.e. creating cavitation with ultrasound. Blake (1999) writes that ‘Intensive insonification of liquids can generate cavitation bubbles.’ Although cavitation occurs usually at elevated sound levels, some insight can still be obtained by investigating the linearized sound equations as manifested in the Helmholtz equation. As is well known, the linear sound wave approximation can be obtained from the Navier–Stokes equations by neglecting viscosity and higher-order terms with |$\rho \partial \boldsymbol{u} / \partial t = - \nabla p $| and the continuity equation |$\partial \rho /\partial t + \nabla \cdot (\rho \boldsymbol{u}) = 0 $|⁠. Using the speed of sound |$c$| defined as |$c^2 = \partial p / \partial \rho $| and if we assume all variables have an |$e^{i\omega t}$| time dependence, we finally get the wave equation for the velocity potential |$\phi $| in the frequency domain as $$\begin{align} \nabla^2 \phi(\boldsymbol{x}) + k^{2} \phi(\boldsymbol{x}) = 0. \end{align}$$ (20) The parameter |$k$| is the wave number, |$k=\omega /c$|⁠, with |$\omega $| the angular frequency and |$c$| the speed of sound. Equation (20) is also an elliptic equation and the framework developed in Section 2.1 can be straightforwardly applied once more. One small difference is that in the Helmholtz equation, |$\phi $| is now a complex variable. It can also easily be seen that if |$\phi $| is a solution then |$\phi e^{i \alpha }$| is also a solution of Eq.20, with |$\alpha $| a constant. We can use this to reconstruct the sound field as a function of time. The Green’s function for (20) is |$G^k = e^{i k r}/r$|⁠, with |$r=|\boldsymbol{x}-\boldsymbol{x_0}|$|⁠, instead of |$G=1/r$| for the potential problem. What remains to be done is to find an equivalent of (11) for the Helmholtz problem. $$\begin{align} \chi(\boldsymbol{x}) = \phi(\boldsymbol{x_0}) \cos{y} +\frac{1}{k}\frac{\partial \phi (\boldsymbol{x_0})} {\partial n} \sin{y}\qquad\qquad\quad\ \ \end{align}$$ (21a) $$\begin{align}\ \frac{\partial \chi(\boldsymbol{x})}{\partial n} = \big[-k \phi(\boldsymbol{x_0}) \sin{y} + \frac{\partial \phi (\boldsymbol{x_0})}{\partial n} \cos{y} \big] \boldsymbol{n}(\boldsymbol{x_0})\cdot \boldsymbol{n}(\boldsymbol{x}) \end{align}$$ (21b) $$\begin{align} y = k\boldsymbol{n}(\boldsymbol{x_0}) \cdot (\boldsymbol{x} -\boldsymbol{x_0})\qquad\qquad\qquad\qquad\ \ \, \end{align}$$ (21c) One possible choice is shown in (21), which essentially is the solution for two standing waves, one with the node and the other with the anti-node at |$\boldsymbol{x_0}$|⁠, both standing waves are normal to the surface at |$\boldsymbol{x_0}$| (and satisfy the Helmholtz equation). In the limit |$k\rightarrow 0$|⁠, the solution of (21) approaches that of (11), as can be shown by a simple Taylor expansion. Other choices are possible, see Sun et al. (2015). The advantage of the use of (21) is that, when used for external problems, again a term |$4 \pi \phi (\boldsymbol{x_0})$| appears on the left-hand side of the boundary element equations, due to the contributions from the integral at infinity, as was the case for the Laplace equation. In analogy to the ‘bubbles near boundaries’ articles of Blake et al. (1986, 1987), we will discuss some results of similar Helmholtz problems. An acoustic transducer is often used to generate ultrasound (and acoustic cavitation). In Fig. 8, such a bowl-shaped transducer with a parabolic internal surface (vibrating up and down) is simulated using a fully desingularized 3D BEM. As in any other BEM method, only the surface of the transducer is meshed. The field values are obtained using post-processing. The instantaneous as well as the absolute values for the potential and pressure are shown. The pressure is proportional to the potential for acoustic problems (up to an imaginary constant). The profiles of the absolute pressure are relatively smooth. Let us now investigate the influence of a nearby free surface. This can be implemented with a modified Green’s function as $$\begin{align} G^k=\frac{e^{ik|\boldsymbol{x}-\boldsymbol{x_0}|}}{|\boldsymbol{x}-\boldsymbol{x_0}|} - \frac{e^{ik|\boldsymbol{x}-\boldsymbol{x_0^I}|}}{|\boldsymbol{x}-\boldsymbol{x_0^I}|}. \end{align}$$ (22) Fig. 8. Open in new tabDownload slide a) The wave field surrounding a parabolic or bowl shaped 3D acoustic transducer facing upwards. The transducer is oscillating up and down in this particular implementation, i.e. |$ \partial \phi / \partial n = \boldsymbol{n} \cdot \boldsymbol{e_z}$| a) The instantaneous real potential (a movie file, Video 2, can be found in the supplementary materials) b) The absolute pressure, note the focal point and the relatively smooth pressure distribution. The simulations were done with a fully desingularized 3D BEM with quadratic elements. Fig. 8. Open in new tabDownload slide a) The wave field surrounding a parabolic or bowl shaped 3D acoustic transducer facing upwards. The transducer is oscillating up and down in this particular implementation, i.e. |$ \partial \phi / \partial n = \boldsymbol{n} \cdot \boldsymbol{e_z}$| a) The instantaneous real potential (a movie file, Video 2, can be found in the supplementary materials) b) The absolute pressure, note the focal point and the relatively smooth pressure distribution. The simulations were done with a fully desingularized 3D BEM with quadratic elements. Note that a minus sign appears in the above equation to account for the fact that the potential on the free surface is zero. Figure 9 shows the results for this configuration. The free surface is placed at the top of the figure. As can be seen, even though the free surface is relatively far from the transducer, the interference patterns caused by this free surface modify the sound field quite considerably, resulting in interference patterns, which are spaced at half a wavelength from each other. Fig. 9. Open in new tabDownload slide As in Fig. 8, but now near a free surface (situated right at the top of the figure) (a) the instantaneous potential profile (which is proportional to the pressure) and (b) the absolute pressure profile. Note the influence of the reflections that cause interference patterns, which are absent in Fig. 8. The simulations were done using an image method to account for the free surface, clearly showing the influence of a nearby surface. A movie file showing the instantaneous potential as a function of time can be found in the supplementary materials (Video 3). Fig. 9. Open in new tabDownload slide As in Fig. 8, but now near a free surface (situated right at the top of the figure) (a) the instantaneous potential profile (which is proportional to the pressure) and (b) the absolute pressure profile. Note the influence of the reflections that cause interference patterns, which are absent in Fig. 8. The simulations were done using an image method to account for the free surface, clearly showing the influence of a nearby surface. A movie file showing the instantaneous potential as a function of time can be found in the supplementary materials (Video 3). The relevance of these simulations with respect to the subject of bubble dynamics becomes obvious in Fig. 10, which shows that bubbles can be ‘trapped’ inside the interference patterns created by the interaction of the transducer and the free surface. It is highly likely that similar phenomena might occur for certain biomedical related applications of focused ultrasound. Fig. 10. Open in new tabDownload slide Experiment with an ultrasound transducer near a free surface with ‘rings of bubbles’ trapped in the interference patterns. On the left the transducer, on the right a zoom-in of the trapped bubbles. The distance between two rings is about half a wavelength (3–4 mm); the diameter of the rings is about 50 mm. Waveform Generator: Agilent 33500B Series (250kHz, 500 mVpp, 0 V offset and 0 deg phase). Amplifier: AG series, T&C Power Conversion Inc. The setup is illuminated with a 250 lumen light source. An acrylic water tank with dimensions |$245\times 245\times 288$| mm was used with water height 160 mm; the distance of the transducer to the water surface was 115 mm. A movie file is available as supplementary material (Video 4). Fig. 10. Open in new tabDownload slide Experiment with an ultrasound transducer near a free surface with ‘rings of bubbles’ trapped in the interference patterns. On the left the transducer, on the right a zoom-in of the trapped bubbles. The distance between two rings is about half a wavelength (3–4 mm); the diameter of the rings is about 50 mm. Waveform Generator: Agilent 33500B Series (250kHz, 500 mVpp, 0 V offset and 0 deg phase). Amplifier: AG series, T&C Power Conversion Inc. The setup is illuminated with a 250 lumen light source. An acrylic water tank with dimensions |$245\times 245\times 288$| mm was used with water height 160 mm; the distance of the transducer to the water surface was 115 mm. A movie file is available as supplementary material (Video 4). Even though the simulations shown in this section are strictly speaking only valid for small pressure perturbations, they can still give us some insight into the complex phenomena that occur. 5. The future of bubble dynamics simulations In their 1987 review, Blake & Gibson (1987) state that ‘Theoretical calculations yield far greater detail on the dynamics of cavitation bubble collapse than can be recorded through experiment’. Even though this was written in 1987, and experimental high speed cameras have improved in frame rate, resolution and ease of use, the same statement is essentially still valid today. Numerical simulations can reveal if and where a jet develops and at which speed (experimentally, jets are very hard to see). Furthermore, pressure contours and other flow details can easily be obtained from a numerical method and lead to greater insight of the underlying physics. Although single bubble collapse is now relatively well understood, there still remain some discoveries to be made. For example, the peculiar stretching of blood cells situated near oscillating bubbles. This phenomenon was investigated numerically with BEM (Tandiono et al., 2013) and it was found that the elasticity of the blood cell plays a crucial part in this strange behaviour. The study of multiple bubbles (e.g. bubble clouds, Ma et al., 2018; Wilson et al., 2008; Klaseboer & Bruin, 1992) is still a subject that can benefit from numerical simulations, e.g. in the research of cavitation on ship propellers (Wilson et al., 2014; PF Pelz & Gross, 2017) and kidney stone lithotripsy treatment (Klaseboer et al., 2007; Sankin et al., 2005; Pishalnikov et al., 2003). At the time of this writing, there are many emerging medical treatments involving bubbles, with perhaps lithotripsy (Klaseboer et al., 2007) the most widely used. Histotripsy, the use of high-intensity ultrasound for localized homogenization of tissues, is used for clinical trials or animal studies in treating prostate cancer and in the removal of renal masses or calculi (Roberts et al., 2006; Roberts, 2014). The bubbles generated by the high-intensity ultrasound at the region of focus oscillate and collapse violently. As a result, the targeted tissues are fractioned or homogenized. The same technique is now being explored for non-invasive surgeries such as thrombolysis (Maxwell et al., 2009), trans-cranial tissue removal (Kim et al., 2014), liver ablation (Vlaisavljevich et al., 2013, 2016)and tissue engineering and regeneration (Khokhlova et al., 2015). Blake himself contributed to an article on ultrasonic cavitation near a tissue layer towards the end of his career (Curtiss et al., 2013). Another medical area where bubbles are involved is the use of ultrasound contrast agents (UCAs) for imaging, diagnosis and treatment (Goldberg et al., 1994). UCAs (Stride & Saffari, 2003) are micron-sized gas or air bubbles with a coating (such as lipid or polymer) that will oscillate, move or collapse in an ultrasound field. These bubbles are approved to be used to image blood perfusion in organs and to measure the blood flow rate in the heart (or other organs). There are emerging techniques to use UCA to obtain super-resolution of blood flow in vasculature (Christensen-Jeffries et al., 2015). A resolution of under 20 |$\mu $|m can be achieved (Christensen-Jeffries et al., 2015). The main idea is to use the bubbles as point scatterers in ultrasound sub-diffraction imaging. A similar technique in optical sub-diffraction imaging uses localization of the bleached fluorophores (Munck et al., 2012). While high resolution is achieved using standard clinical ultrasound equipment with sub-diffraction imaging techniques, there are other attempts to modify the hardware to get higher resolutions. One method is to manipulate the ultrasound beams by various beamforming methods (e.g. minimum variance beamformerand delay-and-sum beamformer) to achieve over 20 fold resolution gains (Diamantis et al., 2017). High-resolution ultrasound imaging may be a good clinical alternative for diagnosis because it is relatively safe, cheap and has small footprint (in comparison with magnetic resonance imaging or computed tomography scans). Many studies attempt to employ targeted UCA for diagnostics or treatment (Wang et al., 2015; Hsiao et al., 2010; Manmi, 2015; Tsigklifis & Pelekasis, 2013). Ligands or protein coatings can be applied on these bubbles. They can then be used for drug and gene delivery (Ferrara et al., 2007), to detect inflammation (Klibanov et al., 2006) or to detect or lyse cancer cells (Frauscher et al., 2001; Yong et al., 2014). To manipulate the UCA in these applications, the bubble dynamics in tissue or elastic materials have to be understood. Viscoelastic models to represent biological tissues have been developed (Yang & Church, 2005; Gaudron et al., 2015; Doinikov & Marmottant, 2018). These models are based on various forms of Voigt viscoelastic models. Recently, Hamaguchi & Ando (2015) studied the bubble oscillation in a gelatin gel experimentally with high-speed photography. They focus on the linear oscillation of a single laser-generated gas bubble. However, in most clinical applications, bubbles occur in clouds and they are not oscillating linearly. Therefore, more experimental and numerical works have to be done to fully understand the complex bubble behaviours in biological systems. From a physics point of view, the forces acting on bubbles in an acoustic field (acoustic radiation forces or Bjerknes forces) are of great importance and remain an active research area (Spratt et al., 2015). The BEM framework as described in Section 4.3 would be very suited to simulate such forces. Detailed understanding of the forces allows better manipulations of the bubbles. This may lead to breakthroughs in many medical applications involving sound waves and bubbles. 6. Conclusions Even though the behaviour of a single bubble is now understood reasonably well, clusters of bubbles are still poorly understood. Also, there is still a lot of work to be done concerning the emergence of a variety of exciting biomedical applications involving bubbles and sound waves. In many aspects, our understanding of biomedical treatments involving bubbles is still in its infancy. This is probably where the pioneering work of Blake and his co-workers will be of great impact to further enhance our knowledge. Most likely the BEM will continue to play an important role in the understanding of such systems. References Azam , F. I. , Karri , B. , Ohl , S. W. , Klaseboer , E. & Khoo , B. C. ( 2013 ) Dynamics of an oscillating bubble in a narrow gap . Phys. Rev. E , 88 , 043006 . Google Scholar Crossref Search ADS WorldCat Becker , A.A. The Boundary Element Method in Engineering: A Complete Course . McGraw-Hill Companies , 1992 . Benjamin , T. B. & Ellis , A. T. ( 1966 ) The collapse of cavitation bubbles and the pressures thereby produced against solid boundaries . Philos. Trans. Roy. Soc. A , 260 , 221 – 240 . Google Scholar Crossref Search ADS WorldCat Best , J. P. & Blake , J. R. ( 1994 ) An estimate of the Kelvin impulse of a transient cavity . J. Fluid Mech. , 261 , 75 – 93 . Google Scholar Crossref Search ADS WorldCat Blake , J. R. ( 1971 ) A note on the image system for a stokeslet in a no-slip boundary . Math. Proc. Camb. Philos. Soc , 70 , 303 – 310 . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Blake , J. R. ( 1983 ) The Kelvin impulse: applications to bubble dynamics . Eight Australasian Fluid Mechanics Conference , University of Newcastle, N.S.W. , New South Wales, Australia. 28 November/2 December 1983 , pp. 10B.1 – 10B.4 . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Blake , J. R. ( 1988 ) The Kelvin impulse: application to cavitation bubble dynamics . J. Aust. Math. Soc. , 30 , 127 – 146 . Google Scholar Crossref Search ADS WorldCat Blake , J. R. ( 1999 ) Preface to acoustic cavitation and sonoluminescence. A theme compiled and edited by J. R. Blake . Philos. Trans. Roy. Soc. A , 357 , 201 – 201 . Google Scholar Crossref Search ADS WorldCat Blake , J. R. , Boulton-Stone , J. M. & Tong , R. P. ( 1995 ) Boundary integral methods for rising, bursting and collapsing bubbles . BE Applications in Fluid Mechanics , vol. 4 , pp. 31 – 71 . Computational Mechanics Publications Southampton , Southampton, England . Blake , J. R. & Chwang , A. T. ( 1974 ) Fundamental singularities of viscous flow . J. Engrg. Math. , 8 , 23 – 29 . Google Scholar Crossref Search ADS WorldCat Blake , J. R. & Gibson , D. C. ( 1981 ) Growth and collapse of a vapour cavity near a free surface . J. Fluid Mech. , 111 , 123 – 140 . Google Scholar Crossref Search ADS WorldCat Blake , J. R. & Gibson , D. C. ( 1987 ) Cavitation bubbles near boundaries . Annu. Rev. Fluid Mech. , 19 , 99 – 123 . Google Scholar Crossref Search ADS WorldCat Blake , J. R. , Hooton , M. C. , Robinson , P. B. & Tong , R. P. ( 1997 ) Collapsing cavities, toroidal bubbles and jet impact . Philos. Trans. Roy. Soc. A , 355 , 537 – 550 . Google Scholar Crossref Search ADS WorldCat Blake , J. R. , Keen , G. S. , Tong , R. P. & Wilson , M. ( 1999 ) Acoustic cavitation: the fluid dynamics of non-spherical bubbles . Philos. Trans. Roy. Soc. A , 357 , 251 – 267 . Google Scholar Crossref Search ADS WorldCat Blake , J. R. , Leppinen , D. M. & Wang , Q. X. ( 2015 ) Cavitation and bubble dynamics: the Kelvin impulse and its applications . Interface Focus , 5 , 20150017. WorldCat Blake , J. R. , Robinson , P. B. , Shima , A. & Tomita , Y. ( 1993 ) Interaction of two cavitation bubbles with a rigid boundary . J. Fluid Mech. , 255 , 707 – 721 . Google Scholar Crossref Search ADS WorldCat Blake , J. R. , Taib , B. B. & Doherty , G. ( 1986 ) Transient cavities near boundaries. Part 1. Rigid boundary . J. Fluid Mech. , 170 , 479 – 497 . Google Scholar Crossref Search ADS WorldCat Blake , J. R. , Taib , B. B. & Doherty , G. ( 1987 ) Transient cavities near boundaries. Part 2. Free surface . J. Fluid Mech. , 181 , 197 – 212 . Google Scholar Crossref Search ADS WorldCat Boulton-Stone , J. M. & Blake , J. R. ( 1993 ) Gas bubbles bursting at a free surface . J. Fluid Mech. , 254 , 437 – 466 . Google Scholar Crossref Search ADS WorldCat Brennen , C. E. ( 1995 ) Cavitation and Bubble Dynamics . Oxford University Press , Oxford, England. Google Preview WorldCat COPAC Calvisi , M. L. , Lindau , O. , Blake , J. R. & Szeri , A. J. ( 2007 ) Shape stability and violent collapse of microbubbles in acoustic traveling waves . Phys. Fluids , 19 , 047101 . Google Scholar Crossref Search ADS WorldCat Cheng , A. H. D. & Cheng , D. T. ( 2005 ) Heritage and early history of the boundary element method . Eng. Anal. Bound. Elem. , 29 , 268 – 302 . Google Scholar Crossref Search ADS WorldCat Christensen-Jeffries , K. , Browning , R. J. , Tang , M. X. , Dunsby , C. & Eckersley , R. J. ( 2015 ) In vivo acoustic super-resolution and super-resolved velocity mapping using microbubbles . IEEE Trans. Med. Imag. , 34 , 433 – 440 . Google Scholar Crossref Search ADS WorldCat Crum , L. A. ( 1979 ) Surface oscillations and jet development in pulsating bubbles . J. Phys. Colloques , 40 , C8–285 – C8–288 . Google Scholar Crossref Search ADS WorldCat Curtiss , G. A. , Leppinen , D. M. , Wang , Q. X. & Blake , J. R. ( 2013 ) Ultrasonic cavitation near a tissue layer . J. Fluid Mech. , 730 , 245 – 272 . Google Scholar Crossref Search ADS WorldCat Diamantis , K. , Greenaway , A. , Anderson , T. , Jensen , J. A. & Sboros , V. ( 2017 ) Experimental performance assessment of the sub-band minimum variance beamformer for ultrasound imaging . Ultrasonics , 79 , 87 – 95 . Google Scholar Crossref Search ADS PubMed WorldCat Ding , Z. & Gracewski , S. M. ( 1996 ) The behaviour of a gas cavity impacted by a weak or strong shock wave . J. Fluid Mech. , 309 , 183 – 209 . Google Scholar Crossref Search ADS WorldCat Doinikov , A. A. & Marmottant , P. ( 2018 ) Natural oscillations of a gas bubble in a liquid-filled cavity located in a viscoelastic medium . J. Sound Vibration , 420 , 61 – 72 . Google Scholar Crossref Search ADS WorldCat Ferrara , K. , Pollard , R. & Borden , M. ( 2007 ) Ultrasound microbubble contrast agents: fundamentals and application to gene and drug delivery . Annu. Rev. Biomed. Eng. , 9 , pp. 415 – 447 Google Scholar Crossref Search ADS PubMed WorldCat Frauscher , F. , Klauser , A. , EJ Halpern , W. H. & Bartsch , G. ( 2001 ) Detection of prostate cancer with a microbubble ultrasound contrast agent . Lancet , 357 , 1849 – 1850 . Google Scholar Crossref Search ADS PubMed WorldCat Gaudron , R. , Warnez , M. T. & Johnsen , E. ( 2015 ) Bubble dynamics in a viscoelastic medium with nonlinear elasticity . J. Fluid Mech. , 766 , 54 – 75 . Google Scholar Crossref Search ADS WorldCat Gibson , D. C. & Blake , J. R. ( 1982 ) The growth and collapse of bubbles near deformable surfaces . Appl. Sci. Res. , 38 , 215 – 224 . Google Scholar Crossref Search ADS WorldCat Goldberg , B. B. , Liu , J. B. & Forsberg , F. ( 1994 ) Ultrasound contrast agents: a review . Ultrasound Med. Biol. , 20 , 319 – 333 . Google Scholar Crossref Search ADS PubMed WorldCat Gonzalez Avila , S. R. , Lim , S. Y. , Betari , J. , Ohl , S. W. , Klaseboer , E. & Khoo , B. C. ( 2009 ) Interfacial jets: an oscillating bubble collapses near a fluid-fluid interface . Gallery of Fluid Motions, American Physical Society Meeting , New York, USA: Physical Review Fluids. Nov. 22–24, Minneapolis, USA. Google Preview WorldCat COPAC Hamaguchi , F. & Ando , K. ( 2015 ) Linear oscillation of gas bubbles in a viscoelastic material under ultrasound irradiation . Phys. Fluids , 27 , 113103 . Google Scholar Crossref Search ADS WorldCat Harris , P. J. ( 1992 ) A numerical model for determining the motion of a bubble close to a fixed rigid structure in a fluid . Internat. J. Numer. Methods Engrg. , 33 , 1813 – 1822 . Google Scholar Crossref Search ADS WorldCat Hsiao , C. T. , Lu , X. & Chahine , G. ( 2010 ) Three-dimensional modeling of the dynamics of therapeutic ultrasound contrast agents . Ultrasound Med. Biol. , 36 , 2065 – 2079 . Google Scholar Crossref Search ADS PubMed WorldCat Hsieh , D. Y. ( 1972 ) On the dynamics of nonspherical bubbles . J. Basic Eng. , 94 , 655 – 665 . Google Scholar Crossref Search ADS WorldCat Jaswon , M. A. ( 1963 ) Integral equation methods in potential theory. I . Proc. Royal Soc. A , 275 , 23 – 32 . Google Scholar Crossref Search ADS WorldCat Khokhlova , V. A. , Fowlkes , J. B. , Roberts , W. W. , Schade , G. R. , Xu , Z. , Khokhlova , T. D. , Hall , T. L. , Maxwell , A. D. , Wang , Y. N. & Cain , C. A. ( 2015 ) Histotripsy methods in mechanical disintegration of tissue: towards clinical applications . Int. J. Hyperthermia , 31 , 145 – 162 . Google Scholar Crossref Search ADS PubMed WorldCat Kim , Y. , Hall , T. L. , Xu , Z. & Cain , C. A. ( 2014 ) Transcranial histotripsy therapy: a feasibility study . IEEE Trans. Ultrason., Ferroelectr., Freq. Control , 61 , 582 – 593 . Google Scholar Crossref Search ADS WorldCat Klaseboer , E. & de Bruin , G. J. ( 1992 ) Behaviour of a cloud of bubbles filled with vapour and a small amount of noncondensable gas (a theoretical and numerical study) . La Houille Blanche , 545 – 550 . WorldCat Klaseboer , E. , Fong , S. W. , Turangan , C. K. , Khoo , B. C. , Szeri , A. J. , Calvisi , M. L. , Sankin , G. N. & Zhong , P. ( 2007 ) Interaction of lithotripter shockwaves with single inertial cavitation bubbles . J. Fluid Mech. , 593 , 33 – 56 . Google Scholar Crossref Search ADS PubMed WorldCat Klaseboer , E. , Sun , Q. & Chan , D. Y. C. ( 2012 ) Non-singular boundary integral methods for fluid mechanics applications . J. Fluid Mech. , 696 , 468 – 478 . Google Scholar Crossref Search ADS WorldCat Klaseboer , E. , Sun , Q. & Chan , D. Y. C. ( 2017 ) Nonsingular field-only surface integral equations for electromagnetic scattering . IEEE Trans. Antennas Propag. , 65 , 972 – 977 . Google Scholar Crossref Search ADS WorldCat Klibanov , A. L. , Rychak , J. J. , Yang , W. C. , Alikhani , S. , Li , B. , Acton , S. , Lindner , J. R. , Ley , K. & Kaul , S. ( 2006 ) Targeted ultrasound contrast agent for molecular imaging of inflammation in high-shear flow . Contrast Media Mol. Imaging , 1 , 259 – 266 . Google Scholar Crossref Search ADS PubMed WorldCat Leighton , T.G. From seas to surgeries, from babbling brooks to baby scans: The acoustics of gas bubbles in liquids. . 2004 , International Journal of Modern Physics B , 18 ( 25 ), 3267 – 3314 . Liu , Y. J. ( 1999 ) On the simple-solution method and non-singular nature of the bie/bem—a review and some new results . Eng. Anal. Bound. Elem. , 24 , 286 – 292 . WorldCat Liu , Y. J. ( 2009 ) Fast Multipole Boundary Element Method: Theory and Applications in Engineering . Cambridge University Press , Cambridge, England. Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Liu , Y. J. & Rudolphi , T. J. ( 2000 ) New identities for fundamental solutions and their applications to non-singular boundary element formulations . Comput. Mech. , 24 , 789 – 795 . WorldCat Lord Rayleigh , O. M. F. R. S. ( 1917 ) VIII. On the pressure developed in a liquid during the collapse of a spherical cavity . Philos. Mag. , 34 , 94 – 98 . Google Scholar Crossref Search ADS WorldCat Ma , J. S. , Hsiao , C. T. & Chahine , G. L. ( 2018 ) Numerical study of acoustically driven bubble cloud dynamics near a rigid wall . Ultrason. Sonochem. , 40 , 944 – 954 . Google Scholar Crossref Search ADS PubMed WorldCat Manmi , K. M. A. ( 2015 ) Three dimensional acoustic microbubble dynamics near rigid boundary . Ph.D. Thesis , University of Birmingham . Google Preview WorldCat COPAC Maxwell , A. D. , Cain , C. A. , Duryea , A. P. , Yuan , L. , Gurm , H. S. & Xu , Z. ( 2009 ) Noninvasive thrombolysis using pulsed ultrasound cavitation therapy—histotripsy . Ultrasound Med. Biol. , 35 , 1982 – 1994 . Google Scholar Crossref Search ADS PubMed WorldCat Munck , S. , Miskiewicz , K. , Sannerud , R. , Menchon , S. A. , Jose , L. , Heintzmann , R. , Verstreken , P. & Annaert , W. ( 2012 ) Sub-diffraction imaging on standard microscopes through photobleaching microscopy with non-linear processing . J. Cell Sci. , 125 , 2257 – 2266 . Google Scholar Crossref Search ADS PubMed WorldCat Ohl , S. W. , Klaseboer , E. & Khoo , B. C. ( 2015 ) Bubbles with shock waves and ultrasound: a review . Interface Focus , 5 , 20150019 . Google Scholar Crossref Search ADS PubMed WorldCat Pearson , A. , Blake , J. R. & Otto , S. R. ( 2004a ) Jets in bubbles . J. Engrg. Math. , 48 , 391 – 412 . Google Scholar Crossref Search ADS WorldCat Pearson , A. , Cox , E. , Blake , J. R. & Otto , S. R. ( 2004b ) Bubble interactions near a free surface . Eng. Anal. Bound. Elem. , 28 , 295 – 313 . Google Scholar Crossref Search ADS WorldCat PF Pelz , T. K. & Gross , T. F. ( 2017 ) The transition from sheet to cloud cavitation . J. Fluid Mech. , 817 , 439 – 454 . Google Scholar Crossref Search ADS WorldCat Pishalnikov , Y. A. , Sapozhnikov , O. A. , Bailey , M. R. , Williams Jr., J. C. , RO Cleveland , T. C. & Crum , L. A. ( 2003 ) Cavitation bubble cluster activity in the breakage of kidney stones by lithotripter shockwaves . J. Endourol. , 17 , 435 – 446 . Google Scholar Crossref Search ADS PubMed WorldCat Plesset , M. S. & Chapman , R. B. ( 1971 ) Collapse of an initially spherical vapour cavity in the neighbourhood of a solid boundary . J. Fluid Mech. , 47 , 283 – 290 . Google Scholar Crossref Search ADS WorldCat Plesset , M. S. & Prosperetti , A. ( 1977 ) Bubble dynamics and cavitation . Annu.Rev. Fluid Mech. , 9 , 145 – 185 . Google Scholar Crossref Search ADS WorldCat Pozrikidis , C. ( 2002 ) A Practical Guide to Boundary Element Methods with the Software Library BEMLIB . CRC Press . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Roberts , W. W. ( 2014 ) Development and translation of histotripsy: current status and future directions . Curr. Opin. Urol. , 24 , 104 . Google Scholar Crossref Search ADS PubMed WorldCat Roberts , W. W. , Hall , T. L. , Ives , K. , Wolf , J. S. , Fowlkes , J. B. & Cain , C. A. ( 2006 ) Pulsed cavitational ultrasound: a noninvasive technology for controlled tissue ablation (histotripsy) in the rabbit kidney . J. Urol. , 175 , 734 – 738 . Google Scholar Crossref Search ADS PubMed WorldCat Robinson , P. B. , Blake , J. R. , Kodama , T. , Shima , A. & Tomita , Y. ( 2001 ) Interaction of cavitation bubbles with a free surface . J. Appl. Phys. , 89 , 8225 – 8237 . Google Scholar Crossref Search ADS WorldCat Sankin , G. N. , Simmons , W. N. , Zhu , S. L. & Zhong , P. ( 2005 ) Shock wave interaction with laser-generated single bubbles . Phys. Rev. Lett. , 95 , 034501. WorldCat Sato , K. , Tomita , Y. & Shima , A. ( 1994 ) Numerical analysis of a gas bubble near a rigid boundary in an oscillatory pressure field . J. Acoust. Soc. Am. , 95 , 2416 – 2424 . Google Scholar Crossref Search ADS WorldCat Shima , A. , Tomita , Y. , Gibson , D. C. & Blake , J. R. ( 1989 ) The growth and collapse of cavitation bubbles near composite surfaces . J. Fluid Mech. , 203 , 199 – 214 . Google Scholar Crossref Search ADS WorldCat Spratt , K. S. , Lee , K. M. , Wilson , P. S. & Hamilton , M. F. ( 2015 ) Radiation damping of an arbitrary shaped bubble . J. Acoust. Soc. Am. , 137 , 2254 – 2254 . Google Scholar Crossref Search ADS WorldCat Stride , E. & Saffari , N. ( 2003 ) Microbubble ultrasound contrast agents: a review . Proc. Inst. Mech. Eng. H , 217 , 429 – 447 . Google Scholar Crossref Search ADS PubMed WorldCat Sun , Q. , Klaseboer , E. , Khoo , B. C. & Chan , D. Y. C. ( 2014 ) A robust and non-singular formulation of the boundary integral method for the potential problem . Eng. Anal. Bound. Elem. , 43 , 117 – 123 . Google Scholar Crossref Search ADS WorldCat Sun , Q. , Klaseboer , E. , Khoo , B. C. & Chan , D. Y. C. ( 2015 ) Boundary regularized integral equation formulation of the Helmholtz equation in acoustics . Royal Society open science , 2 (1), 140520. WorldCat Symm , G. T. ( 1963 ) Integral equation methods in potential theory. II . Proc. Royal Soc. A , 275 , 33 – 46 . Google Scholar Crossref Search ADS WorldCat Taib , B. B. ( 1985 ) Boundary integral method applied to cavitation bubble dynamics . Ph.D. Thesis , University of Wollongong, Australia . Google Preview WorldCat COPAC Tandiono , T. , Klaseboer , E. , Ohl , S. W. , Ow , D. S. W. , Choo , A. B. H. , Li , F. & Ohl , C. D. ( 2013 ) Resonant stretching of cells and other elastic objects from transient cavitation . Soft Matter , 9 , 8687 – 8696 . Google Scholar Crossref Search ADS WorldCat Tong , R. P. , Schiffers , W. P. , Shaw , S. J. , Blake , J. R. & Emmony , D. C. ( 1999 ) The role of splashing in the collapse of a laser-generated cavity near a rigid boundary . J. Fluid Mech. , 380 , 339 – 361 . Google Scholar Crossref Search ADS WorldCat Tsigklifis , K. & Pelekasis , N. A. ( 2013 ) Simulations of insonated contrast agents: saturation and transient break-up . Phys. Fluids , 25 , 032109 . Google Scholar Crossref Search ADS WorldCat Vlaisavljevich , E. , Greve , J. , Cheng , X. , Ives , K. , Shi , J. , Jin , L. , Arvidson , A. , Hall , T. , Welling , T. H. , Owens , G. , et al. ( 2016 ) Non-invasive ultrasound liver ablation using histotripsy: chronic study in an in vivo rodent model . Ultrasound Med. Biol. , 42 , 1890 – 1902 . Google Scholar Crossref Search ADS PubMed WorldCat Vlaisavljevich , E. , Kim , Y. , Allen , S. , Owens , G. , Pelletier , S. , Cain , C. , Ives , K. , Ives , K. , & Xu , Z. ( 2013 ) Image guided non-invasive ultrasound liver ablation using histotripsy: feasibility study in an in vivo porcine model. . Ultrasound Med. Bio. , 39 , 1398 – 1409 . Google Scholar Crossref Search ADS WorldCat Wang , Q. X. ( 1998 ) The evolution of a gas bubble near an inclined wall . Theor. Comput. Fluid Dyn. , 12 , 29 – 51 . Google Scholar Crossref Search ADS WorldCat Wang , Q. X. , Manmi , K. & Calvisi , M. L. ( 2015 ) Numerical modeling of the 3D dynamics of ultrasound contrast agent microbubbles using the boundary integral method . Phys. Fluids , 27 , 022104 . Google Scholar Crossref Search ADS WorldCat Wang , Q. X. , Yeo , K. S. , Khoo , B. C. & Lam , K. Y. ( 1996 ) Nonlinear interaction between gas bubble and free surface . Comput. & Fluids , 25 , 607 – 628 . Google Scholar Crossref Search ADS WorldCat Wilson , M. , Blake , J. R. & Haese , P. M. ( 2008 ) A potential multipole theory for the hydrodynamics of bubble clouds . IMA J. Appl. Math. , 73 , 556 – 577 . Google Scholar Crossref Search ADS WorldCat Wilson , M. , Leppinen , D. M. , Haese , P. M. & Blake , J. R. ( 2014 ) Modelling of solid spheres and bubbles near a tip-vortex and a vortex ring . IMA J. Appl. Math. , 79 , 124 – 162 . Google Scholar Crossref Search ADS WorldCat Yang , X. & Church , C. C. ( 2005 ) A model for the dynamics of gas bubbles in soft tissue . J. Acoust. Soc. Am. , 118 , 3595 – 3606 . Google Scholar Crossref Search ADS PubMed WorldCat Yong , C. L. L. , Ow , D. S. W. , Tandiono , T. , Heng , L. L. M. , Chan , K. K.-K. , Ohl , C.-D. , Klaseboer , E. , Ohl , S.-W. , and Choo , A. B.-H. ( 2014 ) Microbubble-mediated sonoporation for highly efficient transfection of recalcitrant human B-cell lines . Biotechnol. J. , 9 , 1081 – 1087 . Google Scholar Crossref Search ADS PubMed WorldCat Zhang , S. , Duncan , J. H. & Chahine , G. L. ( 1993 ) The final stage of the collapse of a cavitation bubble near a rigid wall . J. Fluid Mech. , 257 , 147 – 181 . Google Scholar Crossref Search ADS WorldCat © The Author(s) 2019. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) TI - Blake, bubbles and boundary element methods JO - IMA Journal of Applied Mathematics DO - 10.1093/imamat/hxz032 DA - 2020-04-26 UR - https://www.deepdyve.com/lp/oxford-university-press/blake-bubbles-and-boundary-element-methods-dFySmFScly SP - 1 VL - Advance Article IS - DP - DeepDyve ER -