1. Introduction
The collapse of a gas or vapour bubble near a solid boundary forms a liquid jet that impinges on the boundary (Plesset & Chapman Reference Plesset and Chapman1971). The impinging jet can damage, or clean, the boundary surface due to the pressure and wall shear stress induced at the boundary (Dijkink & Ohl Reference Dijkink and Ohl2008; Koukouvinis et al. Reference Koukouvinis, Strotos, Zeng, Gonzalez-Avila, Theodorakakos, Gavaises and Ohl2018). In addition to the jet impingement, shock waves emitted during collapse events may contribute to the effects on the boundary, however, it is not clear in which regimes the jet or shock waves are more dominant (van Wijngaarden Reference van Wijngaarden2016). Studying the collapse of single bubbles is a valuable tool for understanding the mechanics involved as the effects of the jet are the result of individual bubble collapses (Benjamin & Ellis Reference Benjamin and Ellis1966). An understanding of individual bubble collapses can also inform investigation of cavitation erosion caused by many individual collapse events building up over time (Fernandez Rivas et al. Reference Fernandez Rivas, Betjes, Verhaagen, Bouwhuis, Bor, Lohse and Gardeniers2013). Some common areas in which these bubbles occur are cavitation damage (Sreedhar, Albert & Pandit Reference Sreedhar, Albert and Pandit2017), hydraulic systems (Luo, Ji & Tsujimoto Reference Luo, Ji and Tsujimoto2016) and ultrasonic cleaning (Verhaagen & Fernández Rivas Reference Verhaagen and Fernández Rivas2016).
Developing a better understanding of the effects of cavitation collapse can increase the lifespan of components such as ship propellers (van Terwisga et al. Reference van Terwisga, van Wijngaarden, Bosschers and Kuiper2007) and tidal turbines (Kumar & Saini Reference Kumar and Saini2010) or determine the effectiveness of ultrasonic cleaning for complex geometries (Verhaagen, Zanderink & Fernández Rivas Reference Verhaagen, Zanderink and Fernández Rivas2016; Reuter et al. Reference Reuter, Lauterborn, Mettin and Lauterborn2017). There are also numerous applications in biomedical fields such as reducing tissue damage during surgery (Palanker, Vankov & Miller Reference Palanker, Vankov and Miller2002), investigating mechanisms of cell death in cases where cavitation could be used for drug delivery (Dijkink et al. Reference Dijkink, Le Gac, Nijhuis, van den Berg, Vermes, Poot and Ohl2008), using cavitation to facilitate needle-free injections (Oyarte Gálvez et al. Reference Oyarte Gálvez, Fraters, Offerhaus, Michel, Hunter and Fernández Rivas2020) and studying the contribution of cavitation in traumatic brain injuries (Canchi et al. Reference Canchi, Kelly, Hong, King, Subhash and Sarntinoranont2017) where cavitation is induced by high accelerations (Pan et al. Reference Pan, Kiyama, Tagawa, Daily, Thomson, Hurd and Truscott2017). More novel applications of cavitation include producing high-speed liquid jets at small scales (Karri et al. Reference Karri, Avila, Loke, O'Shea, Klaseboer, Khoo and Ohl2012), using cavitation as an ice-breaking mechanism (Cui et al. Reference Cui, Zhang, Wang and Khoo2018) and understanding biological mechanisms such as those employed by the snapping shrimp (Versluis et al. Reference Versluis, Schmitz, Von der Heydt and Lohse2000; Shimu et al. Reference Shimu, Yang, Junqi and Changchun2019).
The majority of prior experimental research has focused on the jet dynamics near simple flat boundaries (Benjamin & Ellis Reference Benjamin and Ellis1966; Plesset & Chapman Reference Plesset and Chapman1971; Kucera & Blake Reference Kucera and Blake1990; Dijkink & Ohl Reference Dijkink and Ohl2008; Supponen et al. Reference Supponen, Obreschkow, Tinguely, Kobel, Dorsaz and Farhat2016) or cases with limited complexity such as axisymmetric boundaries (Tomita et al. Reference Tomita, Robinson, Tong and Blake2002) or parallel boundaries (Han et al. Reference Han, Zhu, Guo, Liu and Ni2018; Gonzalez-Avila et al. Reference Gonzalez-Avila, van Blokland, Zeng and Ohl2020). Although understanding how bubbles behave near simple geometries is important, in reality geometries have many more complex features such as corners, indentations, slots and surface imperfections. For this reason further investigation of complex geometries is important.
Some more complex geometries have been studied, such as parallel boundaries closed at one end (Brujan, Takahira & Ogasawara Reference Brujan, Takahira and Ogasawara2019), semi-infinite boundaries (Kucera & Blake Reference Kucera and Blake1990), near combinations of a free surface and an inclined flat boundary (Zhang et al. Reference Zhang, Zhang, Wang and Cui2017), inside a set of concave corners (Kucera & Blake Reference Kucera and Blake1990; Brujan et al. Reference Brujan, Noda, Ishigami, Ogasawara and Takahira2018; Tagawa & Peters Reference Tagawa and Peters2018; Wang et al. Reference Wang, Mahmud, Cui, Smith and Walmsley2020) and inside rectangular and triangular channels (Molefe & Peters Reference Molefe and Peters2019). A recent experimental study has investigated how a ridge-patterned structure can affect the bubble surface evolution and jet dynamics when a bubble is generated on the ridge and between the ridges (Kim & Kim Reference Kim and Kim2020). There have also been investigations into related phenomena such as the behaviour of ultrasonically driven bubble clouds near larger bubbles trapped on pit geometries (Stricker et al. Reference Stricker, Dollet, Fernández Rivas and Lohse2013) and bubble collapse dynamics in microfluidic systems of channels with various shapes (Zwaan et al. Reference Zwaan, Le Gac, Tsuji and Ohl2007).
Studies on more complex geometries typically focus on interior geometries and exploit analytical tools such as the method of images (Kucera & Blake Reference Kucera and Blake1990; Tagawa & Peters Reference Tagawa and Peters2018; Molefe & Peters Reference Molefe and Peters2019) or fully resolve the flow using complicated simulations (Wang et al. Reference Wang, Mahmud, Cui, Smith and Walmsley2020). There are also some studies that only investigate the geometry experimentally (Zhang et al. Reference Zhang, Zhang, Wang and Cui2017). Harris (Reference Harris1996) presents a simplified numerical model to determine the Kelvin impulse for the treatment of general complex geometries, however, the results are not validated. In this research we present and validate a simple numerical model that can be used to treat complex geometries and present results for a complex exterior geometry.
A liquid jet impinging on a boundary can be characterised by its strength and direction. In this paper we investigate the effect of a slot in a flat plate (shown in figure 1a) on the jet direction. Slot geometries are common and some other geometries in flat surfaces could be modelled as a series of slots. Examples of geometries that could be approximated by slots are scratches in flat surfaces, trenches in semiconductor manufacturing and various 3D printed objects. These slots could impact how well the surfaces are cleaned or cause concentrations of cavitation damage.
Figure 1. (a) Schematic view of a slot in the surface of a flat boundary constructed from layered acrylic. (b) The parameters defining a slot, bubble position and jet direction.
The paper is structured as follows. The problem is defined in § 2 and some qualitative predictions are made. The experimental procedure is outlined in § 3.1 and the numerical method is defined in § 3.2, including some key mathematical derivations. The numerical method is then employed in § 4.1 to make some quantitative predictions of the jet direction. Experimental results are presented in § 4.2, and subsequently compared to the numerical predictions in § 4.3.
2. Problem definition
We define a slot as a rectangular channel in a surface, as shown in figure 1(a). The slot has width 
$W$ and height 
$H$. A bubble is positioned horizontally at a distance 
$X$ from the slot centre and vertically at a distance 
$Y$ from the boundary surface. The bubble radius can be neglected when considering only the jet direction, which is discussed further in § 3.2.3. The jet direction is measured anticlockwise from the downwards direction such that a positive angle is towards the right side of the slot. These definitions are shown in figure 1(b). The jet angle is a function of the other four parameters, as defined in (2.1). From dimensional analysis, this function can be reduced to a function of the three non-dimensional variables 
$x = 2X/W$, 
$y = Y/W$ and 
$h = H/W$. The horizontal position of the bubble is normalised with respect to half of the width of the slot so that the bubble is directly above the edge of the slot at 
$x = 1$. The vertical position of the bubble is normalised with respect to the width as this provides the most versatility. A slot with an infinite height can still be regarded as a slot, whereas a slot with infinite width is no longer a slot. Normalising with only the width allows the non-dimensional vertical position to retain relevance for all geometries that could be considered as slots. The jet angle, defined in both forms, is thus
\begin{equation} \theta = f(X, Y, H, W) = g(x, y, h). \end{equation} The physical mechanisms affecting the jet direction can be qualitatively understood. When the bubble collapses it draws in fluid from the surroundings. When the bubble collapses in infinite fluid, with no nearby boundaries, the fluid is drawn in completely symmetrically, leading to a spherically symmetrical collapse. If a solid boundary is present, the fluid cannot be drawn directly through the boundary so the boundary is ‘impeding’ the flow (Blake Reference Blake and Antonia1983). This impedance means that the fluid on the boundary side of the bubble is slower than on the open fluid side. For this reason the bottom of the bubble collapses slower than the top, creating an impulse towards the boundary known as the Kelvin impulse (Blake Reference Blake and Antonia1983). Using this idea of relative impedance, some predictions can be made about how the bubble should behave near the slot. The slot contains fluid, so it should be easier for the bubble to draw fluid from the slot than from a solid boundary, effectively having a lower relative impedance. Thus, the fluid on the slot side of the bubble will move more quickly than on the solid boundary side. This means that the jet should be directed away from the slot. An example of a bubble collapse near a slot is shown in figure 2; note how the jet is directed away from the slot. Based on symmetry and limiting behaviour, two more predictions can be made. First, based on symmetry, when the bubble is above the centre of the slot, the jet should be directed straight down (
$\theta = 0$). Second, as the bubble moves infinitely far from the slot (
$x \to \infty$) the boundary becomes a simple flat boundary and the jet must also be directed straight down. Based on these predictions, it is expected that there will be some maximum 
$\theta$ between 
$x=0$ and 
$x \to \infty$ with a negative minimum of an equal magnitude on the opposite side of the slot.
Figure 2. Snapshots from a high-speed recording of a bubble collapse near a slot with width 
$W = 2.2$ mm and 
$H = 2.7$ mm with the bubble positioned at a vertical distance 
$Y = 2.29$ mm and horizontal distance 
$X = -2.03$ mm. The jet angle is measured to be 
$\theta = -0.099$ radians (
$-5.7^{\circ }$). The slot outline is indicated by the black dashed line. This movie may be found in the supplementary material available at https://doi.org/10.1017/jfm.2020.552.
3. Methods
3.1. Experimental method
Experiments were conducted using laser-induced cavitation (Lauterborn Reference Lauterborn1972; Noack & Vogel Reference Noack and Vogel1999). Slot geometries for a range of slot widths and heights were created by layering laser cut acrylic, as shown in figure 1(a). The experiment configuration is shown in figure 3. The slot geometry was connected to an arm supported by a translation stage. The translation stage was able to translate the geometry in three dimensions to an accuracy of 
$5\ \mathrm {\mu }\textrm {m}$. The slot geometry was submerged approximately 
$25$ mm deep in a 
$180\ \textrm {mm}\times 180\ \textrm {mm} \times 100\ \textrm {mm}$ acrylic tank of water. The plates were 50 mm wide and the slots extended 50 mm along the plate. Bubbles ranged in diameter from approximately 1 mm to 3 mm. Around the jet angle peaks the bubbles were at least 10 bubble radii away from the plate edges. The relatively large distances to the tank walls, plate edges, and plate holder meant that their effects on the jet angle were negligible. Bubbles were at most 16 bubble radii away from the slot centre horizontally, and ranged between 1.25 and 8.6 bubble radii from the slot surface vertically. Prior to conducting experiments, the water was degassed by subjecting it to a near vacuum in a vacuum chamber for approximately 30 min. A Q-switched Nd:YAG laser (‘Bernoulli PIV’ from Litron Lasers) was used to generate a 6 ns pulse at a wavelength of 532 nm. The laser output energy was set to 140 mJ per pulse, and then modulated down by adjusting the attenuator setting (typically between 30 % and 35 %) to use the minimal level that could produce a bubble. The laser pulse was passed through a beam expander as it left the laser and reflected downwards by a mirror. It then passed through a Nikon Plan Fluor 10X microscope objective (
$\text {Numerical Aperture NA} = 0.30$) to focus the laser to the position where the bubble was created. The high numerical aperture, in combination with the expanded laser pulse, ensured that the laser could use minimal power and would not create any secondary bubbles along the laser path (Sinibaldi et al. Reference Sinibaldi, Occhicone, Alves Pereira, Caprini, Marino, Michelotti and Casciola2019). A Photron FASTCAM SA-X2 high-speed camera was used to record the bubble collapse. The camera recorded at 100 000 frames per second and used a 105 mm Nikon Micro-Nikkor lens. A 550 nm longpass filter was used to protect the camera from the laser. The camera triggered the laser to ensure accurate synchronising of the recording and the bubble collapse. For each bubble collapse, 100 frames were recorded, spanning 1 ms of time. A 100 W LED panel was used to backlight the bubble collapse such that the bubble appears dark on a light background in the recorded images. Figure 2 shows snapshots from one such recording.
Figure 3. Diagram of the experiment configuration. In addition to the equipment shown, a high-speed camera observes the bubble collapse, which is backlit by a 100 W LED panel.
Recordings were post-processed in Python to measure key characteristics such as the bubble size variation over time and the direction of the resulting jet. Initially the background was removed from each frame of the recording, and a threshold filter was applied to isolate the bubble. The number of pixels provided a measure of bubble size, and the geometric centre of the pixels provided the position of the bubble centroid. To determine the jet direction, a vector was taken between the centroid positions of the first and second maxima of the bubble size (at 
$99\ \mathrm {\mu }\textrm {s}$ and 
$219\ \mathrm {\mu }\textrm {s}$ in figure 2) which is generally a good approximation of the nominal jet direction at a position and is straightforward to measure (Tagawa & Peters Reference Tagawa and Peters2018). From this analysis, and position measurements from the translation stage, graphs were produced characterising how the jet direction varied with bubble position and geometry characteristics.
3.2. Numerical method
The collapse of a gas or vapour bubble occurs on a time scale where viscous effects are not significant and the water can be considered incompressible for the majority of the collapse. This means that the collapse can be modelled by a velocity potential 
$\phi$. Various methods have been used to solve for this velocity potential. Simple potential flow models using mirror sinks have been used to predict the jet direction in limited geometries (Kucera & Blake Reference Kucera and Blake1990; Tagawa & Peters Reference Tagawa and Peters2018; Molefe & Peters Reference Molefe and Peters2019). More complex boundary integral methods and boundary element methods have been applied, typically to simulate how the surface of the bubble moves during collapse (Kucera & Blake Reference Kucera and Blake1990; Harris Reference Harris1996; Li et al. Reference Li, Han, Zhang and Wang2016; Brujan et al. Reference Brujan, Takahira and Ogasawara2019). In this research a simplified boundary element method is used to predict the jet direction without resolving how the surface of the bubble moves over time. This is similar to the mirror sink model employed by Tagawa & Peters (Reference Tagawa and Peters2018) and Molefe & Peters (Reference Molefe and Peters2019) to compute jet direction but with an infinite distribution of ‘mirror’ sinks along the geometry boundary as in boundary element methods. This method combines the simplicity of the mirror sink model for finding the jet direction with the geometric versatility of boundary element methods. Using a simple model keeps the computational cost low, allowing for analysis of many different geometries across large parameter spaces. The method employed in this research is similar in form to that used by Harris (Reference Harris1996) but without solving for the bubble variation with time. Three key derivations are included here for completeness: the numerical method used (§ 3.2.1), the derivation of the effect of a panel on itself (§ 3.2.2) and the derivation of average bubble surface velocity for considering the effect of bubble size (§ 3.2.3).
3.2.1. Boundary element method
The boundary is modelled as an infinite distribution of sinks on the boundary surface with varying strength. The infinite distribution of sinks is divided into panels such that the strength of the sinks within a panel can be approximated by a constant strength across the whole panel. The strength of each panel is determined by asserting no flow through the boundary surface at the centre of every panel and solving the resulting system of linear equations. The derivation for these equations is shown here.
A sink at position 
$\boldsymbol {x}_s$ induces a velocity 
$\boldsymbol {\nabla } \phi = \boldsymbol {u}$ at any position 
$\boldsymbol {x} \neq \boldsymbol {x}_s$
\begin{equation} \boldsymbol{u} = \frac{m (\boldsymbol{x} - \boldsymbol{x}_s)}{4 {\rm \pi}| \boldsymbol{x} - \boldsymbol{x}_s | ^ 3}, \end{equation}
where 
$m$ is the volume flow rate of the sink, typically referred to as the sink strength.
The flow through a surface with normal vector 
$\boldsymbol {n}$ is the scalar product of the velocity at that point, 
$\boldsymbol {u}$, with the normal vector
\begin{equation} \frac{\partial \phi}{\partial n} = \frac{\partial \phi}{\partial x}n_x + \frac{\partial \phi}{\partial y}n_y + \frac{\partial \phi}{\partial z}n_z = \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}. \end{equation} The velocity through the wall is a combination of the normal velocity from the bubble, 
$\partial \phi _b/\partial n$, and the normal velocity from the wall, 
$\partial \phi _w / \partial n$, which must sum to zero for the no-through-flow condition to be met
\begin{equation} \frac{\partial \phi}{\partial n} = \frac{\partial \phi_b}{\partial n} + \frac{\partial \phi_w}{\partial n} = 0. \end{equation}Thus, the normal component of the velocity from the wall can be expressed as the negative of the normal component of the velocity potential from the bubble.
\begin{equation} \frac{\partial \phi_w}{\partial n} = - \frac{\partial \phi_b}{\partial n}. \end{equation} To determine the normal velocity from the bubble at a position 
$\boldsymbol {x}$, the results from (3.1) and (3.2) can be combined to give
\begin{equation} \frac{\partial \phi_b}{\partial n} = \boldsymbol{u}_b \boldsymbol{\cdot} \boldsymbol{n} = \frac{m_b (\boldsymbol{x} - \boldsymbol{x}_b)}{4 {\rm \pi}| \boldsymbol{x} - \boldsymbol{x}_b | ^ 3} \boldsymbol{\cdot} \boldsymbol{n}, \end{equation}
where 
$\boldsymbol {x}_b$ is the position of the bubble.
Similarly, the velocity contribution from the wall can be expressed with the integral of the boundary sink strength density 
$\sigma$ over the wall surface 
$W$.
\begin{equation} \frac{\partial \phi_w}{\partial n} = \iint_W \sigma (\boldsymbol{x}_w) \frac{\boldsymbol{x} - \boldsymbol{x}_w}{4 {\rm \pi}|\boldsymbol{x} - \boldsymbol{x}_w|^3} \boldsymbol{\cdot} \boldsymbol{n} \,\mathrm{d}W, \end{equation}
where 
$\boldsymbol {x}_w$ is a position on the wall surface.
The wall surface is expressed as a series of 
$N$ panels, each with a constant sink strength density, 
$\sigma _j$, centroid position, 
$\boldsymbol {x}_j$, and an area, 
$A_j$. Thus, (3.6) is approximated by
\begin{equation} \frac{\partial \phi_w}{\partial n} = \sum_{j=1}^N \sigma_j \frac{ A_j(\boldsymbol{x} - \boldsymbol{x}_j) \boldsymbol{\cdot} \boldsymbol{n}}{4 {\rm \pi}|\boldsymbol{x} - \boldsymbol{x}_j|^3}. \end{equation}Substituting (3.5) and (3.7) into (3.4) yields
\begin{equation} \sum_{j=1}^N \sigma_j \frac{ A_j(\boldsymbol{x} - \boldsymbol{x}_j) \boldsymbol{\cdot} \boldsymbol{n}}{4 {\rm \pi}|\boldsymbol{x} - \boldsymbol{x}_j|^3} = - m_b\frac{ (\boldsymbol{x} - \boldsymbol{x}_b)\boldsymbol{\cdot} \boldsymbol{n}}{4 {\rm \pi}| \boldsymbol{x} - \boldsymbol{x}_b | ^ 3}. \end{equation} This can be rewritten in terms of factors relating to the relative positions of each point. In this research, 
$R$ denotes these factors such that
\begin{gather} R_{j} = \frac{ A_j(\boldsymbol{x} - \boldsymbol{x}_j) \boldsymbol{\cdot} \boldsymbol{n}}{4 {\rm \pi}|\boldsymbol{x} - \boldsymbol{x}_j|^3}, \end{gather}
\begin{gather}R_{b} = \frac{(\boldsymbol{x} - \boldsymbol{x}_b) \boldsymbol{\cdot} \boldsymbol{n}}{4 {\rm \pi}|\boldsymbol{x} - \boldsymbol{x}_b|^3}, \end{gather}and
\begin{equation} \sum_{j=1}^N R_j \sigma_j = - m_b R_b. \end{equation} In order to calculate the sink strength densities, 
$\sigma _j$, (3.11) must be solved at a number of points on the boundary equal to the number of panels. For a reasonable distribution of points, and convenience, these points are selected to be the panel centroids. Defining 
$R_{12}$ as the effect of the second panel at the centroid of the first panel, and so on, we write
\begin{gather} \left[\begin{array}{cccc} R_{11} & R_{12} & \ldots & R_{1N} \\ R_{21} & R_{22} & \ldots & R_{2N} \\ \vdots & \vdots & \ddots & \vdots \\ R_{N1} & R_{N2} & \ldots & R_{NN} \end{array} \right] \boldsymbol{\sigma} = -m_b \left[\begin{array}{c} R_{b1} \\ R_{b2} \\ \vdots \\ R_{bN} \end{array} \right], \end{gather}
\begin{gather}\boldsymbol{\mathsf{R}} \boldsymbol{\sigma} = -m_b \boldsymbol{R}_b, \end{gather}
where each term of the 
$\boldsymbol{\mathsf{R}} $ matrix is
\begin{equation} R_{ij} = \frac{ A_j(\boldsymbol{x}_i - \boldsymbol{x}_j) \boldsymbol{\cdot} \boldsymbol{n}_i}{4 {\rm \pi}|\boldsymbol{x}_i - \boldsymbol{x}_j|^3}, \end{equation}
which is undefined for any value 
$i = j$. This is resolved in § 3.2.2.
As with any system of linear equations, the system in (3.13) could be solved in many ways. However, it is noted that the 
$\boldsymbol{\mathsf{R}} $ matrix is fixed for any given geometry and is entirely independent of bubble position. Thus, if multiple bubble positions need to be evaluated for a single boundary geometry, a single inversion of the 
$\boldsymbol{\mathsf{R}} $ matrix can be used to solve the system for all bubble positions using
\begin{equation} \boldsymbol{\sigma} = -m_b \boldsymbol{\mathsf{R}} ^{-1}\boldsymbol{R}_b. \end{equation}Although other methods would be faster for solving the system for a single position, this method is far more efficient when multiple bubble positions need to be solved.
This method generally performs well, as will be demonstrated in later sections, but is vulnerable to ill-conditioned systems in some cases. However, such systems can be identified using the condition number of the 
$\boldsymbol{\mathsf{R}} $ matrix and avoided.
3.2.2. Panel integral
The normal velocity induced by a panel at its own centroid is required in order to solve (3.12). For the majority of the panel this poses no issue; the velocity induced by a point that is not at the centroid is entirely tangential to the panel and so does not contribute to the normal velocity. However, the velocity induced by the centroid itself is undefined. Following the example of Brebbia & Dominguez (Reference Brebbia and Dominguez2001), the singularity can be resolved.
The panel is deformed such that the centroid is expanded into a hemisphere with radius 
$\epsilon$, as shown in figure 4. The velocity normal to the panel at the centroid becomes the integral of the normal velocity induced by the sink at the centroid position over the whole hemisphere.
Figure 4. A cross-section of a panel (dashed line) deformed such that the centroid is a hemisphere (solid line), with radius 
$\epsilon$, centred on the panel centroid.
For any point on the hemisphere, the velocity induced by the centroid sink is always normal to the surface. The normal velocity at any such point is thus
\begin{equation} \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n} = \frac{\sigma}{4 {\rm \pi}\epsilon^2}. \end{equation}
For a circle on the hemisphere aligned parallel to the panel, with all points at an angle 
$\varphi$ from the horizontal, the sum of velocities is
\begin{equation} u_{circle} = \frac{\sigma}{4 {\rm \pi}\epsilon^2} 2 {\rm \pi}\epsilon \cos(\varphi) = \frac{\sigma}{2 \epsilon} \cos(\varphi). \end{equation}This is integrated over the whole hemisphere to give
\begin{equation} u_{total} = \int_0^{{\rm \pi} / 2} \frac{\sigma}{2 \epsilon} \cos(\varphi) \epsilon \,\mathrm{d}\varphi = \frac{\sigma}{2} \int_0^{{\rm \pi} / 2} \cos(\varphi) \,\mathrm{d}\varphi = \frac{\sigma}{2}, \end{equation}
which is independent of the hemisphere radius, 
$\epsilon$. As 
$\epsilon \to 0$, the velocity induced by the panel at its own centroid is a constant 
$\sigma / 2$ and so the relative position factor is 
$R = 0.5$.
3.2.3. Average surface velocity
Previous research (Tagawa & Peters Reference Tagawa and Peters2018; Molefe & Peters Reference Molefe and Peters2019) has shown that the jet direction is accurately predicted by the velocity that the mirror sinks induce at the position of the bubble centroid.
The jet direction is aligned with the bubble translation velocity which is defined as the velocity of the centroid. We define the centroid as the average position of the bubble surface, so the velocity of the centroid is the average velocity of the bubble surface. In this section, it is shown that the induced velocity at the centroid of a spherical bubble is equal to the average velocity of the bubble surface, independent of bubble size. A more rigorous version of this derivation may be found in appendix A.
The velocity induced by any individual sink on a point is given by
\begin{equation} \boldsymbol{u} = -\frac{m\boldsymbol{r}}{4 {\rm \pi}|\boldsymbol{r}|^3}, \end{equation}
where 
$\boldsymbol {r}$ is the vector from the sink to the point.
By integrating 
$\boldsymbol {u}$ over a sphere it can be shown that the total velocity on the surface of the sphere is
\begin{equation} \boldsymbol{u}_{total} = -m \frac{\mathcal{R}^2\boldsymbol{d}}{|\boldsymbol{d}|^3}, \end{equation}
where 
$\mathcal {R}$ is the radius of the sphere and 
$\boldsymbol {d}$ is the vector from the sink to the bubble centroid.
The average velocity on the surface of the sphere is thus the total velocity divided by the surface area of the sphere,
\begin{equation} \boldsymbol{u}_{avg} = \frac{\boldsymbol{u}_{total}}{4 {\rm \pi}\mathcal{R}^2} = -\frac{m \mathcal{R}^2 \boldsymbol{d}}{4 {\rm \pi}\mathcal{R}^2 |\boldsymbol{d}|^3} = -\frac{m \boldsymbol{d}}{4 {\rm \pi}|\boldsymbol{d}|^3}, \end{equation}which is equal to the velocity induced by the sink at the bubble centroid.
Due to the properties of potential flow, the velocity induced by every sink is simply a linear summation of the velocity induced by each individual sink. Thus, the average surface velocity of a sphere induced by a combination of sinks is equal to the velocity at the centre of the sphere. It is noted that the bubble sink itself does not contribute to the average surface velocity as all components cancel out.
4. Results and discussion
4.1. Numerical results
Using the boundary element method described in § 3.2, the jet direction can be modelled for given geometric parameters. Up to 20 000 panels were used to model the boundary, with a distribution such that there were more panels near the slot than towards the edges of the plate. The Python implementation of the boundary element method developed for this research runs on a standard desktop computer, where the primary limitation is the memory required to store the matrices.
Figure 5 shows a contour plot of the jet angle for a slot with a square cross-section (
$h = 1$). As predicted qualitatively, the jet is directed towards the boundary but angled away from the slot. Three main regimes are revealed. Close to the boundary but far from the slot (low 
$y$, high 
$x$) the jet angle is dominated by the flat boundary and is not significantly affected by the slot. This leads to a very low jet angle. Near the centre of the slot the geometric symmetry dominates, again leading to low jet angles. The third regime is between these two, where the bubble is far enough from the centre that symmetry does not dominate, and positioned such that the slot has a significant effect relative to the flat boundary.
Figure 5. (a) A contour plot of 
$\theta$ as a function of 
$x$ and 
$y$. The slot, with 
$h = 1$, is represented by the black line, with points plotted at the panel centroids. Very close to the boundary the jet angle is more strongly affected by individual panel sinks rather than the boundary as a whole leading to the more rapid variations visible near the boundary. (b) A plot of 
$\theta$ against 
$x$ for 
$y = 1$, corresponding to the black dashed line on the contour plot in (a).
In order to more easily quantify this shape, and for simpler comparison to experimental results, slices of the 
$\theta$ contour are taken with constant 
$y$ and a range of 
$x$. One such slice is shown in figure 5(b) and is the form in which data will be presented hereafter. As in the contour plot, this plot shows that the jet tends to be directed straight at the boundary when the bubble is far from the slot. As the bubble approaches the slot the jet is angled away from the slot centre down to a negative peak jet angle. After the peak the jet angle tends back towards zero, crossing zero at the centre point. The jet angle then increases to an equal and opposite peak on the other side of the slot due to the symmetry of the geometry.
Using the normalised geometric parameters, two sets of numerical predictions have been plotted in figures 6 and 7 to show how the jet angle curve varies with 
$h$ and 
$y$.
Figure 6. (a) Value of 
$\theta$ against 
$x$ for a range of slot geometries characterised by 
$h$ and a fixed value of 
$y = 1$. The slot boundaries are indicated by the dashed vertical lines (
$x = -1$, 
$x = 1$). (b) Normalised 
$\theta$, 
$\hat {\theta }$, against normalised 
$x$, 
$\hat {x}$, for the same values of 
$h$. The curves collapse well onto one curve when normalised. The 
$h = 1$ curves here are the same as the 
$y = 1$ curves in figure 7.
Figure 7. (a) Value of 
$\theta$ against 
$x$ for a range of slot geometries characterised by 
$y$ and a fixed value of 
$h = 1$. The slot boundaries are indicated by the dashed vertical lines (
$x = -1$, 
$x = 1$). (b) Normalised 
$\theta$, 
$\hat {\theta }$, against normalised 
$x$, 
$\hat {x}$, for the same values of 
$y$. The curves collapse well onto one curve when normalised. The 
$y = 1$ curves here are the same as the 
$h = 1$ curves in figure 6.
Figure 6(a) shows that the peak jet angle increases as the normalised slot height, 
$h$, increases. This is because the slot contains more fluid and so its relative impedance is decreased. The position of the peak moves closer to the slot as the height is increased. These curves have a similar shape, which can be characterised by two parameters defining the peak. For a given 
$y$ and 
$h$, the maximum value of 
$\theta$ is 
$\theta ^\star$ and occurs at 
$x = x^\star$. Thus, 
$\theta$ can be normalised with 
$\theta ^\star$, and 
$x$ can be normalised with 
$x^\star$.
\begin{equation} \hat{\theta} = \frac{\theta}{\theta^\star},\quad \hat{x} = \frac{x}{x^\star}. \end{equation}
When these curves are normalised, as shown in figure 6(b), they collapse down to being very close to the same curve, with the exception of the data with the lowest height (
$h = 0.5$).
Figure 7(a) shows that the peak jet angle, 
$\theta ^\star$, increases as the normalised vertical distance away, 
$y$, decreases. This is because the bubble is closer to the boundary and so is more strongly affected. The position of the peak, 
$x^\star$, moves further from the slot as the vertical distance away is increased. When these curves are normalised, as shown in figure 7(b), they also collapse down to being very close to the same curve, with the exception of the data closest to the boundary (
$y = 0.5$).
It is noted that the 
$h = 1$ curves in figure 6 are the same as the 
$y = 1$ curves in figure 7 showing that the curves from both figures collapse to the same curve.
As mentioned, the curves collapse except for bubbles very close to the boundary, or with slots that have a low height. For these cases, the effect of the slot can be treated as two opposing steps with a relatively large separation. A single step would cause a single jet angle peak. When the effects of two peaks with opposite signs are combined, the gradient in the middle depends on the separation of the peaks. If the peaks are very close, then the peaks blend together resulting in the collapsed shapes seen in figures 6(b) and 7(b). If, however, the peaks are far apart then a pronounced kink in the curve appears at the middle of the slot, as shown by the 
$h = 0.5$ curve in figure 6 and the 
$y = 0.5$ curve in figure 7. When bubbles are close to the boundary the effects of each side of the slot are relatively more separated, as is also the case for low height slots. Most slots have 
$h \geq 1$ and 
$y \geq 1$ where the peaks merge giving the collapsed curve shape; only these curves will be considered in further analysis.
As the curves in figures 6 and 7 collapse well onto the same curve, 
$\theta ^\star$ and 
$x^\star$ can be used to characterise the variation of the jet angle with the parameters 
$h$ and 
$y$. Two contour plots in figure 8 show 
$\theta ^\star$ and 
$x^\star$ plotted as functions of 
$h$ and 
$y$.
Figure 8. Contour plots of peak jet angle, 
$\theta ^\star$, and peak jet angle position, 
$x^\star$, as functions of non-dimensional vertical distance, 
$y$, and non-dimensional slot height, 
$h$.
Figure 8(a) shows that the maximum jet angle, 
$\theta ^\star$, depends most strongly on the dimensionless vertical bubble position, 
$y$. As 
$y$ increases, with the bubble moving far from the slot, the maximum jet angle tends towards zero. As 
$y$ tends towards zero, the maximum jet angle increases very rapidly, but is bounded by the solution for a convex right-angle corner as would be the case for a bubble at 
$y = 0$ as 
$x \to 1^-$.
Figure 8(a) shows a weaker dependence of maximum jet angle on the dimensionless slot height, 
$h$, compared to the dependence on 
$y$. As 
$h$ increases, the maximum jet angle tends towards a limit, as can be observed in the contour, showing only a dependence on 
$y$ for large values of 
$h$. The figure also shows that the maximum jet angle tends to zero as 
$h$ tends to zero which is expected as the slot becomes a flat plate.
Figure 8(b) shows the dependence of the position of the maximum jet angle, 
$x^\star$, on both 
$y$ and 
$h$. 
$x^\star$ increases approximately linearly with 
$y$, as might be predicted from figure 5(a). For low values of 
$h$ the gradient of the linear relationship changes, but as 
$h$ increases the gradient tends towards 
$1$ (this can be seen in figure 13b where the gradient of the 
$x^\star -y$ line is approximately 
$1$). It is noted that, for very low jet angles, the position of the peak becomes more sensitive to numerical errors and imperfect boundary conditions, such as the edges of the plate. Very low jet angles are found at low 
$h$ and high 
$y$ values, in the top left corners of the contour plots in figure 8.
Both 
$\theta ^\star$ and 
$x^\star$ exhibit limiting behaviour as 
$h$ increases. This suggests that at some point increasing the height of the slot will have a negligible effect on the jet angle. Even for lower values of 
$h$, the variation with 
$h$ is typically much less significant than with 
$y$.
4.2. Experimental results
Experiments were conducted by performing a horizontal sweep over a selected slot at fixed vertical distances 
$Y$. The horizontal positions tested were selected to focus most of the data around both 
$\theta$ peaks and to observe the behaviour at a large horizontal distance on at least one side of the slot. Each position was tested multiple times; for most experiments there were five repeats.
In order to understand the stochastic variation of 
$\theta$ that occurs at each position, two series of experiments were conducted with 50 repeats at each position. From these data it was observed that the standard deviation of jet angle is reasonably consistent for all positions and conforms well to the normal distribution. We did not observe a dependence of jet angle variation on the magnitude of the bubble displacement(see supplementary material). This standard deviation was therefore applied to the remainder of the data to provide statistical error bars that show a 99 % confidence interval of the mean at each position based on the number of repeats at those positions.
A second-order polynomial curve fit was applied around each of the two peaks of the 
$\theta -x$ curve. As the geometry is symmetrical the two peaks should be equal and opposite; the peak on the negative 
$x$ side should have a negative 
$\theta$ value of the same magnitude as on the positive 
$x$ side. Thus, if the polynomial curve fits on each side are slightly offset from being symmetrical, all of the data can be shifted by the offset in both 
$\theta$ and 
$x$ to achieve symmetry. Offsets in 
$\theta$ are most likely caused by misalignments in the experiment. Other possible sources of asymmetry would cause offsets of much lower magnitudes than even minor misalignments. In addition to correcting for misalignments, low quality data can be identified by very large curve fit offsets. The curve fit from one side, mirrored in both axes, should fit the data from the opposite side. If this is not the case then it can be concluded that the data are of a low quality and are thus neglected in further analysis. Of the data collected in this research, horizontal sweeps at three 
$y$ values for the W2H6 geometry from table 1 were neglected on this basis and thus not presented here.
Table 1. Measurements of geometries used. Labels refer to the nominal width and height.
Data were gathered for both peaks so that this symmetry analysis could be conducted. However, only one side of the slot was tested to a greater horizontal distance as the behaviour in this region is already well understood. An example sweep with the curve fit plotted is shown in figure 9. In this example the peaks are slightly offset from being symmetrical, so the data would be shifted before conducting further analysis.
Figure 9. An example curve fit on data using slot ‘W2H3a’ from table 1 (
$W = 2.2$ mm, 
$H = 2.7$ mm). The bubble was positioned at a vertical distance 
$Y = 2.81$ mm. The points are the experimental data and the lines are the two curve fits. The slot boundaries are indicated by the dashed vertical lines (
$x = -1$, 
$x = 1$).
The geometries tested are shown in table 1 and the jet angle data are summarised in figure 10. The experimental results follow the same qualitative trends as the numerical predictions. There is a negative jet angle peak on the left side of the slot, and a positive jet angle peak on the right. The magnitude of jet angle peak increases as 
$y$ decreases and 
$h$ increases. The position of the jet angle peak increases as 
$y$ increases, but it is more difficult to discern how the position varies with 
$h$.
Figure 10. Value of 
$\theta$ plotted against 
$x$ for all experimental data from table 1. The slot boundaries are indicated by the dashed vertical lines (
$x = -1$, 
$x = 1$). Data are coloured by 
$y$ value; peak 
$\theta$ decreases as 
$y$ increases. (a) W1H3: blue 
$y = 1.58$, orange 
$y = 2.37$, green 
$y = 3.16$. (b) W2H3a: blue 
$y = 0.80$, orange 
$y = 1.04$, green 
$y = 1.28$, red 
$y = 1.51$, purple 
$y = 1.75$. (c) W2H3b: blue 
$y = 1.21$, orange 
$y = 1.67$. (d) W2H6: blue 
$y = 0.69$, orange 
$y = 0.90$. (e) W2H9: blue 
$y = 0.78$, orange 
$y = 1.24$. (f) W2H12: blue 
$y = 1.20$. (g) W4H12: blue 
$y = 0.58$, orange 
$y = 0.82$.
4.3. Comparison
We will now proceed to directly compare our experimental data to the numerical results.
Figure 11 shows a direct comparison between four experimental data curves and boundary element method predictions for the same geometric parameters. These plots show a good agreement between predicted curves and experimental results. The most significant difference is for the W4H12 geometry where the experimental data have a steeper gradient on the outer sides of the peaks, but the magnitude and position of the peak are well predicted. It is noted that the W4H12 data in general have a steeper peak curve than other experimental data when compared to numerical results, as can be observed in the supplementary material where all experimental data are compared to numerical results. There is a tendency for the numerical model to under-predict the magnitude of the jet angle peak, although the peak position is generally predicted well. On average, across all 17 horizontal sweeps presented here, the numerical model under-predicted 
$\theta ^\star$ by 
$13$ % and 
$x^\star$ by 
$2.6$ %.
Figure 11. Experimental data for a series of slots from table 1. These data are compared to numerical predictions for the same configurations. Experimental data are plotted as points with error bars, numerical model predictions are plotted as solid lines. The slot boundaries are indicated by the dashed vertical lines (
$x = -1$, 
$x = 1$). Geometries are ordered by 
$h$ value. (a) W2H3a, 
$h = 1.23$, 
$y = 1.28$. (b) W1H3, 
$h = 2.23$, 
$y = 2.37$. (c) W4H12, 
$h = 2.73$, 
$y = 0.58$. (d) W2H12, 
$h = 5.23$, 
$y = 1.20$. Comparison plots for all geometries may be found in the supplementary material.
All experimental data collapse onto a single curve when the experimental data are normalised with the peak values 
$x^\star$ and 
$\theta ^\star$, which were determined using the previously described curve fitting method. Figure 12 shows all of the data normalised and compared to the normalised prediction curve that matches all predictions with 
$h \geq 1$ and 
$y \geq 1$. The experimental data collapse very well, validating the collapse observed from the numerical predictions. Although the numerical prediction curve has a slightly higher normalised jet angle on the outer sides of the peaks, the collapsed experimental data curve matches the numerical curve remarkably well. The variations in the normalised curve observed numerically at very low 
$h$ and at very low 
$y$ would likely not be visible in the experimental data due to the magnitude of the error compared to the magnitude of the variations.
Figure 12. Experimental results from all geometries, each with a range of 
$y$, showing normalised jet angle, 
$\hat {\theta }$, against normalised horizontal position, 
$\hat {x}$. This is compared to the collapsed numerical curve that matches all numerical curves for 
$y \geq 1$ and 
$h \geq 1$. The experimental results collapse down onto one curve and match the numerical prediction.
A comparison between numerically predicted 
$\theta ^\star$ and 
$x^\star$ trends and experimental results is shown in figure 13. Here, 
$\theta ^\star$ and 
$x^\star$ have been calculated using the curve fitting method described in § 4.2. The error bars in this figure are based on the error distribution from large amounts of data, synthesised using the numerical model, with similar properties to experimental data. It is also noted from the synthesised data that the curve fit has a tendency to over-estimate the peak position. The results are compared for the same 
$h$ values as the experimental data and the same range of 
$y$ values as figure 8. These results generally show a good agreement between the numerical and experimental results. The most significant difference is that the numerical prediction underestimates the 
$\theta ^\star$ values for W1H3. This is to be expected as a similar discrepancy is observed in figure 11.
Figure 13. Experimental data points from two slots compared with boundary element method predictions for the same configurations. (a) Peak jet angle, 
$\theta ^\star$, as a function of 
$y$. (b) Position of peak jet angle, 
$x^\star$, as a function of 
$y$. Slot W1H3 has 
$h = 2.23$. Slot W2H3a has 
$h = 1.23$. Experimental data are plotted as dots, numerical model predictions are solid lines.
Overall, the model tends to under-predict the jet angle, but performs well on the curve shape and trends. The position of the maximum jet angle, 
$x^\star$, is especially well predicted.
5. Conclusion
We have investigated the collapse of bubbles near a slot geometry using experiments and a simple numerical model. Our main observation is the variation of the jet angle with the horizontal position of the bubble. At the centre of the slot the jet is directed straight downwards due to symmetry and far from the slot the surface acts as a flat boundary so the jet is also directed downwards. Between these two limits there is a peak jet angle deflection, angled away from the slot. The peak jet angle and position of the peak jet angle both tend to a limiting value as the slot height increases. This shows that, for slots of sufficient height, the jet angle depends only on the bubble position. As the vertical position of the bubble increases, the peak jet angle decreases and the horizontal position of the peak increases. For 
$h \gg 1$, the position of the peak jet angle is directly proportional to the vertical distance: 
$x^\star \propto y$. When the jet angle and horizontal position are normalised by their respective peak values, we find that all jet angle curves collapse.
In this research we do not present any quantitative analysis of the strength of the jet, however, we can make some qualitative predictions. For a simple flat boundary, jet strength depends on the distance between the bubble and the boundary. For a slot there is the distance to the slot surface and the distance to the slot floor. Far from the slot the jet strength would depend on the distance from the slot surface. Above an infinitely wide slot the jet strength would depend on the distance from the slot floor. We may conclude that near a finite slot the jet strength would be between these two limiting cases but quantifying the exact behaviour will require further experiments.
The numerical model has a tendency to under-predict the jet angle compared to experimental data gathered in this research, although it is often within reasonable error. The numerical model predictions very closely follow the shape and overall trends of the data, and provide a good prediction of the position of the peak jet angle, 
$x^\star$. In addition, the collapsed curve predicted by the numerical model matches the collapsed curve found from the experimental data. These comparisons serve to validate the numerical model presented in this research and provide a good basis from which to continue the study of complex geometries using this model. The velocity profile of the fluid in and around the slot can also be predicted by this model, although further investigation would be required to validate the predicted profiles.
In the context of cleaning with bubbles, particularly with ultrasonic cavitation, this research suggests that slots in surfaces to be cleaned would likely experience fewer jet impacts because bubbles would be drawn away from the slots rather than into them, and thus be cleaned less rigorously than the rest of the surface. Where cavitation damage is a problem this property could also be used to protect sensitive components by recessing them within slots, although this technique would likely have other implications depending on the flow conditions and requirements of the components.
Although we do not quantify the cleaning or damage effects here, these areas have been extensively studied previously using simple geometries. Our research provides a valuable step towards applying this knowledge to the more complex geometries that are found in many real-world applications.
Acknowledgements
Many thanks to L. Molefe for invaluable discussions and comments. We acknowledge financial support from the EPSRC under Grant No. EP/P012981/1. D.F.R. acknowledges the funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 851630). Data and code supporting this study are openly available from the University of Southampton repository at https://doi.org/10.5258/SOTON/D1454.
Declaration of interests
The authors report no conflict of interest.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2020.552.
Appendix A. Full average surface velocity derivation
The translation velocity of the bubble is defined as the average velocity of the bubble surface and can be used to calculate the jet direction as the jet direction follows the translation of the bubble. However, many evaluations of the velocity field would be required to determine the average velocity at the bubble surface numerically, so a different method is preferred. From the method employed by Tagawa & Peters (Reference Tagawa and Peters2018) and Molefe & Peters (Reference Molefe and Peters2019), and from numerical testing, it is found that the velocity induced by mirror sinks at the bubble centre accurately describes the jet direction, independent of bubble radius. It is of interest to more fully understand this relationship and how it relates to the average surface velocity.
The potential flow model is a linear summation of multiple sinks, including the bubble and the wall or mirror sinks. Thus, the average velocity on the bubble surface for the whole system can be described as the summation of the average velocity on the bubble surface for each individual sink. It is assumed that the bubble surface is a perfect sphere, centred on the bubble position, with a radius 
$\mathcal {R}$. The average velocity on the bubble surface induced by its own sink is zero as each point has an equal and opposite point that sum to zero so the bubble sink can be neglected from further analysis.
The velocity induced by a sink on an arbitrary point is directed at the sink and has a calculable magnitude. In the present analysis each sink is evaluated in isolation, so it can be assumed that the sink is positioned at the origin and the bubble is positioned at a distance 
$k$ along an axis from the sink to the bubble centre, as shown in figure 14(a). This reduces the problem to an axisymmetric configuration where the off-axis components of the velocity will cancel out. Thus, only the velocity in the axial direction need be evaluated. For this analysis the axis is denoted by 
$y$ and the axes normal to it are denoted by 
$x$ and 
$z$ in the usual manner. It is assumed that the bubble radius, 
$\mathcal {R}$, is less than the distance from the sink, 
$|\boldsymbol {r}|$, such that the bubble is not in contact with the sink. The velocity along the 
$y$ axis is
\begin{equation} u_y = -\frac{m}{4 {\rm \pi}} \frac{y}{|\boldsymbol{r}|^3}, \end{equation}
where 
$m$ is the sink strength, 
$y$ is the 
$y$ position of a point to be evaluated and 
$\boldsymbol {r}$ is the position vector of the point to be evaluated.
Figure 14. (a) The arrangement of the bubble, sink and the direction of the 
$y$ axis. (b) The definition of 
$\varphi$ and 
$\mathcal {R}$. (c) A section of the sphere surface used for integration.
For this derivation 
$\varphi$ is defined as the angle anticlockwise from the horizontal towards the 
$y$ axis, as shown in figure 14(b). At 
$\varphi =-{\rm \pi} /2$ the point is at the bottom of the sphere, and at 
$\varphi ={\rm \pi} /2$ it is at the top of the sphere. Thus, for a point on the surface of the bubble, the 
$y$ position is defined as
\begin{equation} y = k + \mathcal{R} \sin(\varphi). \end{equation} The magnitude of the 
$\boldsymbol {r}$ vector is
\begin{align} |\boldsymbol{r}| & = \sqrt{y^2 + \mathcal{R}^2 \cos^2(\varphi)} \nonumber\\ & = \sqrt{k^2 + 2k\mathcal{R}\sin(\varphi) + \mathcal{R}^2\sin^2(\varphi) + \mathcal{R}^2\cos^2(\varphi)} \nonumber\\ & = \sqrt{k^2 + \mathcal{R}^2 + 2k\mathcal{R}\sin(\varphi)}. \end{align}
Note that 
$\mathcal {R}^2\cos ^2(\varphi )$ includes both the 
$x$ and 
$z$ components of the 
$r$ vector.
Thus, (A 1) can be expanded as
\begin{equation} u_y = - \frac{m}{4 {\rm \pi}} \frac{k + \mathcal{R}\sin(\varphi)}{(k^2 + \mathcal{R}^2 + 2k\mathcal{R}\sin(\varphi))^{3/2}}. \end{equation} The sphere surface can be broken into a series of circles, aligned parallel to the 
$xz$ plane. The sum of the velocities around any such circle is described by
\begin{equation} u_{yc} = 2 {\rm \pi}\mathcal{R}_c u_y = 2 {\rm \pi}\mathcal{R} \cos(\varphi) u_y, \end{equation}
where 
$\mathcal {R}_c = \mathcal {R}\cos (\varphi )$ is the radius of the circle.
By considering a section of the sphere surface, as shown in figure 14(c), with a radius 
$\mathcal {R}$ and height 
$\mathcal {R}\,\textrm {d}\varphi$, (A 5) can be integrated over the whole sphere to find the total velocity 
$u_{ytotal}$.
\begin{align} u_{ytotal} & = \int_{-{\rm \pi} / 2}^{{\rm \pi} / 2} u_{yc} \mathcal{R} \,\mathrm{d}\varphi \end{align}
\begin{align} & = - \frac{m}{2} \mathcal{R}^2 \int_{-{\rm \pi} / 2}^{{\rm \pi} / 2} \frac{\cos(\varphi) (k + \mathcal{R}\sin(\varphi))}{(k^2 + \mathcal{R}^2 + 2k\mathcal{R}\sin(\varphi))^{3/2}} \,\mathrm{d}\varphi \end{align}
\begin{align} & = \left.- \frac{m}{2} \mathcal{R}^2 \frac{\mathcal{R} + k\sin(\varphi)}{k^2 \sqrt{k^2 + \mathcal{R}^2 + 2k\mathcal{R}\sin(\varphi)}} \right|_{-{\rm \pi} / 2}^{{\rm \pi} / 2} \end{align}
\begin{align} & = - \frac{m}{2} \mathcal{R} ^2 \frac{2}{k^2} = - m \frac{\mathcal{R}^2}{k^2}. \end{align} The average velocity is thus the total velocity 
$u_{ytotal}$ divided by the total area 
$4 {\rm \pi}\mathcal {R}^2$. Substituting in the expression from (A 9) yields
\begin{equation} u_{yavg} = \frac{u_{ytotal}}{4 {\rm \pi}\mathcal{R}^2} = - \frac{m \mathcal{R}^2}{4 {\rm \pi}\mathcal{R}^2 k^2} = - \frac{m}{4 {\rm \pi}k^2}, \end{equation}
which is notably independent of the bubble radius 
$\mathcal {R}$. Equation (A 10) is the same as the 
$y$ velocity evaluated at the centre of the bubble, found by substituting 
$y = k$ and 
$|\boldsymbol {r}| = k$ into (A 1).
Thus, the velocity induced at the centre of the bubble is equal to the average surface velocity. This can be extended to the whole system where the average surface velocity is equal to the velocity induced at the centre of the bubble by all sinks.
\begin{equation} \boldsymbol{u} = \sum_{k=1}^N \sigma_k \frac{ A_k(\boldsymbol{x}_b - \boldsymbol{x}_k)}{4 {\rm \pi}|\boldsymbol{x}_b - \boldsymbol{x}_k|^3}. \end{equation}
Only the direction of 
$\boldsymbol {u}$ is considered, as the magnitude is effectively arbitrary in this research.
