2.1. Governing equations

Figure 1. Schematic of the dimensionless two-bubble problem. The fluid domain is denoted by $\Omega$ and the the bubble surfaces are $\partial \Omega _{1,2}$ . We supply a uniform outer flow far from the bubbles. The bubble centre–centre distance is $\sigma$ and the angle the bubbles make to the direction of the outer flow is $\phi$ .

As in Booth et al. (Reference Booth, Griffiths and Howell2023), we consider the motion of two bubbles in a Hele-Shaw cell of height $\hat {h}$ parallel to the $(\hat {x}, \hat {y})$ -plane. Here, $\hat {h}$ is assumed to be much smaller than the horizontal dimensions of the cell and the bubbles, so we can employ lubrication theory. The bubbles are flattened by the cell walls above and below into pancake-like shapes with approximately circular profiles when viewed from above (figure 1), whose radii are denoted by $\hat {R}_1$ and $\hat {R}_2$ , where $\hat {R}_{1,2}\gg \hat {h}$ . We prescribe a uniform unidirectional flow with velocity $\hat {U}\boldsymbol {i}$ in the far field (where $\boldsymbol {i}$ denotes the unit vector in the $\hat {x}$ -direction). The viscosity of the liquid and the liquid–air surface tension are denoted by $\hat {\mu }$ and $\hat {\gamma }$ , respectively.

We non-dimensionalise the system by scaling lengths with $\hat {R}_1$ , velocities with $\hat {U}$ , the fluid pressure $\hat {p}$ with $12\hat {\mu }\hat {U}\hat {R}_1/{\hat {h}}^2$ and the pressure inside the $k$ th bubble, $\hat {p}_k$ , with $2\hat {\gamma } / \hat {h}$ , where $\hat {\gamma }$ is the surface tension. This procedure gives the following dimensionless model, in which dimensionless variables are denoted without hats:

(2.1a) \begin{align} \nabla ^2 p & = 0 \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\text {in} \quad \Omega , \end{align}

(2.1b) \begin{align} \boldsymbol {n}\boldsymbol {\cdot } \boldsymbol {\nabla } p & = -U_{n,k} \qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\!\text{on} \quad\!\partial \Omega _{k}, \end{align}

(2.1c) \begin{align} p_k - \frac {3\mathrm {Ca}}{\epsilon } p & = 1 + \mathrm {Ca}^{2/3}\beta \left (U_{n,k}\right ) \left (U_{n,k}\right )^{2/3} + \frac {\epsilon \pi }{4}\kappa _k \quad\quad\,\,\,\mbox {on} \quad\!\partial \Omega _k, \end{align}

(2.1d) \begin{align} \boldsymbol {\nabla }p&\rightarrow -\boldsymbol {i} \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \text {as} \quad x^2 +y^2 \rightarrow \infty . \end{align}

Here, $\Omega$ is the fluid domain, while $\partial \Omega _k$ , $\kappa _k$ and $U_{n,k}$ are the boundary, in-plane curvature and local normal velocity of the interface of the $k$ th bubble, respectively ( $k=1,2$ ), and $\beta$ is the Bretherton coefficient, whose value depends on whether the meniscus is advancing or retreating (Bretherton Reference Bretherton1961; Wong, Radke & Morris Reference Wong, Radke and Morris1995; Halpern & Jensen Reference Halpern and Jensen2002)

(2.2) \begin{equation} \beta \left (U_{n,k}\right )=\begin{cases} \beta _1 \approx 3.88 & \text { when }U_{n,k} \gt 0,\\ \beta _2 \approx -1.13 & \text { when }U_{n,k} \lt 0. \end{cases} \end{equation}

The boundary condition (2.1c) was proposed by Meiburg (Reference Meiburg1989) and later derived by Burgess & Foster (Reference Burgess and Foster1990). In (2.1b ) we neglect to include the contribution due to leakage into the thin films because this effect is always subdominant. However, this effect could easily be included in the model (see, for example Burgess & Foster Reference Burgess and Foster1990; Peng et al. Reference Peng, Pihler-Puzović, Juel, Heil and Lister2015; Wu et al. Reference Wu, Booth, Griffiths, Howell, Nunes and Stone2024).

The system (2.1) contains two dimensionless parameters: the bubble aspect ratio and the capillary number, defined by

(2.3) \begin{align} \epsilon & = \frac {\hat {h}}{2\hat {R}_1}, & \mathrm {Ca} &= \frac {\hat {\mu } \hat {U}}{\hat {\gamma }}, \end{align}

respectively, both of which are assumed to be small. Specifically, in the distinguished limit $\mathrm {Ca}=O(\epsilon ^3)$ , the viscous pressure balances the pressure drop across the menisci (the second and fourth terms in (2.1c )). In this regime, both bubbles remain circular to leading order, so $U_{n,k}=\boldsymbol {U}_k\boldsymbol {\cdot }\boldsymbol {n}$ and $p$ is therefore fully determined by the problem (2.1a ), (2.1b ) and (2.1d ) (up to an irrelevant constant) once the bubble velocities $\boldsymbol {U}_k$ are specified. As a shortcut to determine the bubble velocities we perform an effective net force balance by integrating (2.1c ) around each bubble (see, for example, Booth et al. Reference Booth, Griffiths and Howell2023), to obtain

(2.4) \begin{equation} \frac {\boldsymbol {U}_k}{|\boldsymbol {U}_k|^{1/3}} =\frac {\delta }{\pi R_k}\oint _{\partial \Omega _k} -p\boldsymbol {n}\, \mathrm {d}s, \end{equation}

where $R_k$ is the dimensionless radius of the $k$ th bubble. The resulting problem contains a single dimensionless group, the Bretherton parameter, defined by

(2.5) \begin{equation} \delta = \frac {1}{\eta }\, \frac {\mathrm {Ca}^{1/3}}{\epsilon } = \frac {2}{\eta }\,\frac {\hat {R}_1}{\hat {h}}\,\left (\frac {\hat {\mu }\hat {U}}{\hat {\gamma }}\right )^{1/3}, \end{equation}

which is assumed to be $O(1)$ while $\epsilon$ and $\mathrm {Ca}$ both tend to zero. The numerical constant

(2.6) \begin{equation} \eta =\frac {(\beta _1-\beta _2)\Gamma (4/3)}{3\sqrt {\pi }\Gamma (11/6)}\approx 0.894 \end{equation}

incorporates the Bretherton pressure drops (2.2) across the advancing and retreating menisci (Bretherton Reference Bretherton1961; Wu et al. Reference Wu, Booth, Griffiths, Howell, Nunes and Stone2024).

The Bretherton parameter is a dimensionless parameter that compares the magnitudes of the viscous pressure from the flow around the bubble and of the Bretherton pressure, or the pressure drop across the thin films surrounding the bubble. As $\delta$ increases to infinity, the viscous pressure dominates over the Bretherton pressure. In this limit, we recover the result due to Taylor & Saffman (Reference Taylor and Saffman1959) that the bubble moves at twice the background flow velocity, which was obtained while disregarding the thin film drag. Booth et al. (Reference Booth, Griffiths and Howell2023) showed that an isolated circular bubble travels parallel to the background flow with velocity $\boldsymbol {U}_b=U_b\boldsymbol {i}$ , where $U_b$ is a monotonically increasing function of $\delta$ , satisfying $U_b\rightarrow 0$ as $\delta \rightarrow 0$ and $U_b\rightarrow 2$ as $\delta \rightarrow \infty$ . Importantly for this work, it follows that larger bubbles travel faster than smaller ones when all other parameters are fixed, since $\delta \propto \hat {R}_1$ . This conclusion may also be drawn using dimensional analysis, through which it can be shown that the driving force due to the background flow is proportional to $\hat {R}_1^2$ and the drag force due to the thin films is proportional to $\hat {R}_1$ .

2.2. Complex variable formulation

We now reformulate the problem (2.1) in terms of complex variables. At leading order we have two circular bubbles whose centroids are at positions $(x_1,y_1)$ and $(x_2,y_2)$ in the $(x,y)$ -plane, with a uniform velocity in the far field of unit magnitude. We label the bubbles such that the smaller bubble is located at $(x_1,y_1)$ and the dimensionless radii of the two bubbles are thus $R_1=1$ and $R_2 =R$ , where $R\geqslant 1$ is the radius ratio of the two bubbles. As shown schematically in figure 1, the problem is instantaneously characterised by the length $\sigma$ of the vector joining the smaller bubble centre to the larger bubble centre, the angle $\phi$ that it makes with the $x$ -axis (which is parallel to the background flow direction), and the radius ratio $R$ .

Since the flow is governed by Laplace’s equation, we can formulate this as a problem for the complex potential $w(z)=-p+\mathrm {i}\psi$ , where $\psi$ is the streamfunction, and $z=x+\mathrm {i}y$ . Then $w(z)$ is holomorphic in the region $\Omega$ outside the two bubbles and satisfies the boundary conditions

(2.7a) \begin{align} \quad\operatorname {Im}[w(z)] & = Q_1 + \operatorname {Im}\left (\overline {\mathcal {U}}_1 z\right ) \qquad\text {on} \quad |z-z_1| = R_1=1, \end{align}

(2.7b) \begin{align} \operatorname {Im}[w(z)]& = Q_2 + \operatorname {Im}\left (\overline {\mathcal {U}}_2 z\right ) \qquad \text {on} \quad |z-z_2| = R_2=R\geqslant 1, \end{align}

(2.7c) \begin{align} & w(z)\sim z +o(1) \qquad\quad \text {as} \quad z\rightarrow \infty , \end{align}

where, for $k \in \{1,2\}$ , we denote by $z_k = x_k +\mathrm {i}y_k$ and $\mathcal {U}_k = U_k +\mathrm {i}V_k$ the complex representations of the $k$ th bubble position and velocity, respectively, and the $Q_k$ are a priori unknown constants. The over-bars denote complex conjugation. Note that (2.7a ) and (2.7c ) are the complex representations of the kinematic boundary conditions (see, for example, Crowdy Reference Crowdy2008).

Once we have solved for $w(z)$ , to close the system we evaluate the effective force balance (2.4) on each bubble, which in complex variables becomes

(2.8a) \begin{align} &\quad \frac {1}{\mathrm {i}\pi }\oint _{\partial \Omega _1} w(z)\,\mathrm {d}z=-\mathcal {U}_1+\frac {1}{\pi }\oint _{\partial \Omega _1}p\mathrm {i}\,\mathrm {d}z =-\mathcal {U}_1+\frac { \mathcal {U}_1}{\delta \left |\mathcal {U}_1\right |^{1/3}}, \end{align}

(2.8b) \begin{align} &\frac {1}{\mathrm {i}\pi }\oint _{\partial \Omega _2} w(z)\,\mathrm {d}z=-R^2\mathcal {U}_2+\frac {1}{\pi }\oint _{\partial \Omega _2}p\mathrm {i}\,\mathrm {d}z =-R^2\mathcal {U}_2+\frac {R \mathcal {U}_2}{\delta \left |\mathcal {U}_2\right |^{1/3}}. \end{align}

Here $\partial \Omega _k$ is the boundary of the $k$ th bubble, given by $|z-z_k|=R_k$ .

The problem for the pressure field generated by two bubbles of unequal radii was solved by Sarig et al. (Reference Sarig, Starosvetsky and Gat2016) using a bipolar coordinate transformation, resulting in infinite series solutions for the interaction forces between the bubbles. Instead, using our complex variable formulation facilitates the evaluation of the integrals (2.8) in the force balance in closed form.

5.1. Observed behaviour

In this section, we consider situations involving two nearly circular bubbles in which the larger bubble is initially far behind the smaller one and offset in the $y$ -direction by a distance less than the sum of the two bubble’s radii, such that $x_1-x_2\gg 1$ and $0\lt |y_2-y_1|\lt 1+R$ . As explained in § 2.1, the larger bubble at the rear is expected to travel faster than the bubble at the front (Booth et al. Reference Booth, Griffiths and Howell2023). Thus, the bubbles would collide if they only moved parallel to the background flow. However, for a range of starting positions, we find that the nonlinear hydrodynamic interaction between the bubbles allows them to avoid collision by rotating around one another. Lubrication forces prevent the collision of the nearly circular bubbles in a manner that is analogous to how they cause rigid spheres or cylinders to rotate around and pass each other without contact in shear flow (see, e.g. Arp & Mason Reference Arp and Mason1977b ). As the larger bubble approaches from behind, it continues along a relatively straight trajectory. It overtakes the smaller bubble, which manoeuvres out of the way to let the larger bubble pass.

Figure 4. Two-bubble rollover with $\delta =1.17$ and $R=2.05$ at different dimensionless times $t = \hat {t}\hat {U}/\hat {R}_1$ . (top) Experimental images are compared with (bottom) simulations of the dimensionless dynamical system (3.17) with the same initial conditions at $t=0$ . The background flow is from left to right. Experimental images have been rescaled by the rear bubble radius, $\hat {R}_1 =$ 2.6 mm, for comparison with the theory.

In figure 4, we show experimental images demonstrating this rollover effect for a system of two approximately circular bubbles with $\delta =1.17$ and $R=2.05$ (see movie 1 provided in the Supplementary Material). The larger bubble catches up to the smaller one, which evades contact by rolling over the larger one. In the lower plots, we demonstrate good qualitative agreement with solutions of the dynamical system (3.17) for the same parameter values and initial conditions. Movies 2–7 in the Supplementary Material show additional instances of the two-bubble rollover phenomenon, serving as evidence that it is reproducible for various combinations of bubble sizes and initial conditions.

Figure 5. The instantaneous bubble velocity components $(U_k,V_k)$ (top and bottom, respectively) versus dimensionless time $t$ for (a) $\delta =1.17$ and $R=2.05$ , (b) $\delta =0.90$ and $R=2.32$ . Solid lines show theoretical predictions and points show experimental data. The bubble of unit dimensionless radius ( $k=1$ ) is shown in blue (circles), and the bubble of radius $R$ ( $k=2$ ) is shown in red (triangles). In each plot, the time at which $x_1=x_2$ is shown with a vertical line. Error bars are comparable to the size of the markers and are thus omitted.

Figure 6. The positions of the bubble centres $(x,y)$ (top and bottom, respectively) versus dimensionless time $t$ for (a) $\delta =1.17$ and $R=2.05$ , (b) $\delta =0.90$ and $R=2.32$ . Solid lines show theoretical predictions, and points show experimental data. The bubble of unit radius ( $k=1$ ) is shown in blue (circles), and the bubble of radius $R$ ( $k=2$ ) is shown in red (triangles). In each plot, the time at which $x_1=x_2$ is shown with a vertical line. Error bars are comparable to the size of the markers and are thus omitted.

Figure 7. Trajectories for the two-bubble dynamical system (3.17) in the reference frame of the smaller bubble, with (a) $\delta =1.17$ and $R=2.05$ , (b) $\delta =0.90$ and $R=2.32$ . The blue vectors show the predicted trajectories of the centre of the larger bubble relative to the smaller one, and the red points show the experimentally measured bubble positions. Error bars are comparable to the size of the markers and are thus omitted. Any trajectories entering the solid grey region $|z_2-z_1| \leqslant (1+R)$ are such that the two bubbles will collide. The solid black region $|z_2-z_1| \leqslant 1$ represents the smaller bubble.

In figure 5, the instantaneous bubble velocity components $(U_k,V_k)$ are plotted for the same experiment as shown in figure 4 and for another example in which $\delta =0.90$ and $R=2.32$ (see movie 2 provided in the Supplementary Material). The model predicts that the smaller bubble decelerates in the $x$ -direction as the large bubble approaches from behind, while also translating in the $y$ -direction such that $|y_2-y_1|$ is increasing. The time at which the pair of bubbles is aligned perpendicularly to the background flow (i.e. when $x_1=x_2$ ) coincides with when the axial velocity of the smaller bubble, $U_1$ , reaches a minimum and when $V_1=V_2=0$ . We observe reasonable agreement between theory and experiment. However, in experiments, the velocity components $U_1$ and $V_1$ of the smaller bubble are generally biased to reduce the distance between the bubble centres.

In figure 6, we compare the experimental and theoretical results for the $(x,y)$ -positions of the bubble centres. In both cases, we observe that the motion of the larger bubble is largely unaffected by the interaction while the smaller bubble is displaced in the $y$ -direction such that the bubbles avoid contact as the larger one passes. The final distance in the $y$ -direction between the bubbles in experiments is significantly smaller than what is theoretically predicted. The bubbles also become slightly closer in the $x$ -direction in experiments as compared with theory. The small discrepancies between the theoretical and experimental velocities shown in figure 5 accumulate over time and lead to noticeable differences between the theoretical and experimental bubble trajectories.

Finally, in figure 7 we plot the trajectories of the centre of the larger bubble relative to that of the smaller one (i.e. $z_2-z_1$ ) calculated using (3.17). Any trajectory entering the solid grey region $|z_2-z_1|\leqslant 1+R$ corresponds to a collision between the bubbles. Points extracted from the experiments are superimposed on the theoretically determined bubble trajectories. We observe that in experiments, the larger bubble initially follows a trajectory, then departs from that trajectory when the two bubbles are close. This departure is likely to be due to interactions between the bubbles that are not included in the model. Finally, as the bubbles separate, the larger bubble once again closely follows a trajectory, albeit a different trajectory from the one on which the bubble started.

While the $x$ -positions of the bubble centres are well captured by the theory, there is a significant disagreement between the predicted and observed $y$ -positions of the smaller bubble during and after the rollover (see figure 6). In the experiments, the smaller bubble appears to be entrained behind the larger bubble such that the distances between their centres in both the $x$ - and $y$ -directions are smaller than the corresponding theoretical trajectories. This process breaks the fore–aft symmetry that is predicted by (3.17), and indeed which would be expected in Stokes flow for circular bubbles. However, as noted in § 4, there are perturbations to the bubble shape due to the background flow, which also happen to be fore–aft asymmetric due to the differences between the advancing and retreating menisci (Wu et al. Reference Wu, Booth, Griffiths, Howell, Nunes and Stone2024). Deformations due to interactions between bubbles are known to cause asymmetric trajectories for unconfined bubbles rising due to buoyancy. Experiments and numerical simulations performed at low Reynolds numbers have shown that a smaller bubble tends to align itself behind a larger bubble and even accelerate towards it so that the two bubbles collide, all while both bubbles undergo significant deformations (Manga & Stone Reference Manga and Stone1993, Reference Manga and Stone1995). It is possible that small inertial effects also play a role: experiments and numerical simulations have shown that a deformable bubble rising due to buoyancy behind another one tends to get drawn into the wake of the latter (Crabtree & Bridgwater Reference Crabtree and Bridgwater1971; Katz & Meneveau Reference Katz and Meneveau1996; Bunner & Tryggvason Reference Bunner and Tryggvason2003; Huisman, Ern & Roig Reference Huisman, Ern and Roig2012). In § 6, we investigate the deviations in the bubble shape of two bubbles in a line parallel to the background flow.

5.2. Do the bubbles collide?

5.2.1. Conditions for a bubble collision

In § 5.1 we found that the bubbles can avoid colliding by rolling over one another. By analysing the dynamical system (3.17), we can predict when or if the bubble rollover effect will occur. We note that the following analysis of bubble–bubble collisions is conducted within the context of the Hele-Shaw model. The Hele-Shaw model will break down when the bubbles are close to contact, in which case a three-dimensional analysis would be needed to achieve a complete description of the dynamics.

At each instant in time, (3.16) determines $\mathcal {U}_1$ and $\mathcal {U}_2$ as functions of $\sigma$ and $\phi$ . We can then update $\sigma$ and $\phi$ using $\mathcal {U}_2-\mathcal {U}_1= (\dot {\sigma }+\mathrm {i}\sigma \dot {\phi } )\mathrm {e}^{\mathrm {i}\phi }$ (where the dot represents differentiation with respect to $t$ ). In figure 8, we plot the phase space showing the resulting trajectories of the larger bubble relative to the smaller one, i.e. $z_2-z_1=\sigma \mathrm {e}^{\mathrm {i}\phi }$ . In this figure, we take $R=2$ for illustration. The solid grey region, $1\lt |z_2-z_1|\leqslant (1+R)$ , corresponds to the region of intersection between the bubbles. The rollover effect occurs on any trajectory that starts from $x_1-x_2\gg 1$ with $0\lt |y_2-y_1|\lt 1+R$ and that does not enter the solid grey region, and the likelihood of observing the effect depends strongly on the value of $\delta$ . In figure 8(a), we show a case where $\delta$ is large, and all suitable initial conditions satisfying the inequalities stated above will give rise to the rollover effect. In this case, the bubbles repel each other so strongly that collision between the bubbles is impossible. On the other hand, in figure 8(b), we show that a smaller value of $\delta$ leads to much weaker interaction between the bubbles, such that the trajectories remain almost parallel to the flow. In this case, the rollover effect can occur only for a very narrow band of initial conditions, and we are much more likely to observe the bubbles colliding with each other. To understand the underlying physical mechanisms, we recall that $\delta$ is a dimensionless parameter that compares the relative magnitudes of the viscous pressure and of the Bretherton pressure, and that in this system $\delta$ is defined using the radius of the smaller bubble, whose motion is essential to successful rollover. As $\delta$ increases, the magnitude of the viscous pressure dominates that of the Bretherton pressure, so the motion of the smaller bubble is less hindered by the Bretherton drag, and collision is less likely.

Figure 8. Trajectories for the two-bubble dynamical system (3.17) in the reference frame of the smaller bubble, with $R=2$ and (a) $\delta =5$ , (b) $\delta =1/2$ . Any trajectories entering the solid grey region $|z_2-z_1| \leqslant (1+R)$ are such that the two bubbles will collide. Stationary points are shown in red. The solid black region $|z_2-z_1| \leqslant 1$ represents the smaller bubble.

In this section, we consider conditions under which the bubbles will collide. First, we observe that there are stationary points (saddle points, located at $\phi = 0$ and $\phi =\pi$ , shown in red) in figure 8(a) but not in figure 8(b). The existence of such stationary points outside of the solid grey region as in figure 8(a) implies that two aligned bubbles (i.e. with $y_1=y_2$ ) will never collide. The stable manifolds of the two saddle points coincide with the horizontal axis, $y_1-y_2 = 0$ , so a point on the horizontal axis also lies on the stable manifold of one of the stationary points. Therefore trajectories beginning on the horizontal axis converge to a stationary point without entering the solid grey region. Furthermore, we find that, in figure 8(a), the trajectories on the surface $|z_2-z_1|=1+R$ with $x_2\gt x_1$ (the larger bubble in front) are directed inwards (into the solid grey region) and for $x_2\lt x_1$ are directed outwards. In this case, bubbles may only collide if they are initially close to each other, and the larger bubble is ahead of the smaller one when the collision occurs. The reverse is true in figure 8(b), in which the surface $|z_2-z_1|=1+R$ is entirely outside of the separatrix connecting the two stationary points.

Motivated by these observations, we examine the following two conditions on the flow:

1. (i) The stationary points of the dynamical system (3.17) in the reference frame of the smaller bubble are in the region $|z_2-z_1|\geqslant 1+R$ .

2. (ii) In a neighbourhood of $x_1=x_2$ , the trajectories point into the region $|z_2-z_1|\leqslant 1+R$ for $x_2\gt x_1$ and out of the region $|z_2-z_1|\leqslant 1+R$ for $x_2\lt x_1$ .

In § 5.2.2 and § 5.2.3, for each condition $k\in \{1,2\}$ , we will find a critical minimum value of $\delta = \delta _k(R)$ . Then, for $\delta \lt \delta _1$ , we argue that there is always a range of initial conditions with $x_1-x_2\gg 1$ and $|y_2-y_1|\lt 1+R$ such that the bubbles collide (including the case $y_2=y_1$ where the bubbles are aligned). On the other hand, for $\delta \gt \delta _2$ , it is impossible for bubbles that start far apart in the $x$ -direction to collide, regardless of their initial transverse separation.

Note that there exists a third critical value of $\delta =\delta _c(R)$ satisfying $\delta _1\leqslant \delta _c\leqslant \delta _2$ , at which the separatrix connecting the two stationary points is tangent to $|z_2-z_1|=1+R$ . This critical value provides a sharp bound on $\delta$ above which collision between two bubbles that are initially well separated in the $x$ -direction (the direction of the background flow) is impossible. However, $\delta _c$ is delicate to compute numerically as it depends on the global properties of the flow whereas, as we will show, the critical values $\delta _1$ and $\delta _2$ can be determined from purely local information about the normal velocity $U_n$ at the collision boundary $|z_1-z_2|=1+R$ .

Figure 9. Position of the stationary point, $\sigma _s$ , as a function of the Bretherton parameter, $\delta$ , for radius ratios $R = 1.5$ (red), 2 (blue), 2.5 (purple). The dashed black curve shows where $\sigma _s=1+R$ .

5.2.2. Condition i: stationary points

If this condition is satisfied, then two aligned bubbles will never collide. By analysing (3.16), we can find the stationary points by solving for $\mathcal {U}_1=\mathcal {U}_2\equiv \mathcal {U}$ and for $(\sigma ,\phi )\equiv (\sigma _s,\phi _s)$ at a fixed $\delta$ . Since each $f_k$ in (3.14) is real, by symmetry we find that $\mathcal {U}=U\in \mathbb {R}$ and $\phi _s =0$ or $\pi$ . We focus on the case $\phi _s=0$ since by symmetry the stationary points are at $(\pm \sigma _s,0)$ . Thus, for given $\delta$ and $R$ we find $U$ and $\sigma _s$ by solving the nonlinear algebraic equations

(5.1a) \begin{align} \quad \bigl (f_1(\sigma _s,R)-f_2(\sigma _s,R)\bigr )(U-1)&=-U+\frac {U^{2/3}}{\delta }, \end{align}

(5.1b) \begin{align} \bigl (f_1(\sigma _s,R)-f_3(\sigma _s,R)\bigr )(U-1)&=-R^2U+\frac {R U^{2/3}}{\delta } ,\end{align}

numerically, using Newton’s method.

The position of the stationary point, $\sigma _s$ , is plotted as a function of the Bretherton parameter, $\delta$ , in figure 9. The black dashed curve shows where $\sigma _s=1+R$ . For each fixed value of $R$ we observe that, for suitably small $\delta$ , there are no stationary points in the region $|z_2-z_1|\geqslant 1+R$ . As $\delta$ is increased, there exists a first value $\delta = \delta _1(R)$ at which a stationary point appears at $\sigma _s = 1+R$ . Then, for $\delta \gt \delta _1$ , $\sigma _s$ is a monotonically increasing function of $\delta$ .

We can find $\delta _1(R)$ by substituting $\sigma _s=1+R$ in (5.1) and solving for $U$ and $\delta _1$ ; the details of this calculation may be found in Appendix A. We plot $\delta _1$ as a function of the bubble radius ratio, $R$ , in figure 10. We observe that $\delta _1$ is a monotonically decreasing function of $R$ , which means that for larger values of $R$ the stationary points are present for smaller values of $\delta$ . We also observe that, as $R\rightarrow 1^+$ , $\delta _1(R)$ tends to a finite value that is approximately $2.37$ .

In figure 11(a), we plot the phase space showing the resulting trajectories of the larger bubble relative to the smaller bubble with $R=2$ for $\delta =\delta _1(2)$ . We observe that the stationary points of the system occur on the real axis at $\sigma _s=1+R$ (shown by red points); however, there are still trajectories that enter the solid grey region $|z_2-z_1|\leqslant 1+R$ . Hence the bubbles can still collide.

Figure 10. Minimum values $\delta _1(R)$ (dashed) and $\delta _2(R)$ (solid) of the Bretherton parameter, $\delta$ , satisfying conditions i (see § 5.2.2) and ii (see § 5.2.3), respectively.

Figure 11. Trajectories for the two-bubble dynamical system (3.17) in the frame of the smaller bubble, with $R=2$ and (a) $\delta =\delta _1(2)$ , at which the stationary points (shown as red points) lie on the surface $|z_2-z_1| =1+R$ , (b) $\delta =\delta _2(2)$ , above which the separatrix encloses the region $|z_2-z_1| \lt 1+R$ (solid grey fill). The solid black region $|z_2-z_1| \leqslant 1$ represents the smaller bubble.

5.3. Do the bubbles collide in finite time?

In figure 8(b), we observe trajectories that enter the solid grey region $|z_2-z_1|\leqslant 1+R$ , which suggests that the bubbles collide. To show that a collision occurs in finite time, we calculate the relative normal velocity $U_n$ of the two bubbles in in the limit when they are touching as $\sigma \rightarrow 1+R$ (see Appendix A for the behaviour of the functions $f_k$ given by (3.14) in this limit). If $U_n\lt 0$ , the bubbles collide in finite time if they start sufficiently close. We plot $U_n$ as a function of $\phi$ in figure 12(a) for $R=2$ and various values of $\delta$ . Figure 12(b) shows a schematic of the two bubbles touching with the definitions of $\boldsymbol {n}$ and $\phi$ .

We find three possible regimes:

1. (i) If $\delta \geqslant \delta _c$ (see § 5.2), then when a trajectory starts inside the separatrix with a non-zero offset in the $y$ -direction, it will result in a collision in finite time (see figure 8 a).

2. (ii) If $\delta _1\lt \delta \lt \delta _c$ , we are in an intermediate regime where $U_n \gt 0$ for parts of both $\phi \in (0,\pi /2)$ and $\phi \in (\pi /2,\pi )$ and the separatrix does not completely enclose the region $|z_2-z_1|\leqslant 1+R$ . In this regime, the stationary points of (3.17) are in the region $|z_2-z_1|\gt 1+R$ . Hence, there exist trajectories with the larger bubble beginning far behind the smaller one ( $x_1- x_2\gg 1$ ) that result in collision in finite time.

3. (iii) If $\delta \leqslant \delta _1$ , we have $U_n\lt 0$ for $\phi \in (\pi /2,\pi )$ . Thus, in configurations where the larger bubble is behind the smaller one ( $x_1 \gt x_2$ ), they collide in finite time provided that the initial value of $|y_2-y_1|$ is not too large. Example trajectories of this kind are observed in figure 8(b).

In figure 12(a) we observe that the values of $\delta _1$ and $\delta _2$ can be determined by the local information about the normal velocity $U_n$ . The first critical value, $\delta _1$ is the value of $\delta$ at which $U_n =0$ at $\phi =0$ and $\pi$ , and the second critical value, $\delta _2$ , is the value of $\delta$ at which $\partial U_n/\partial \phi = 0$ at $\phi =\pi /2$ .

It should be noted that we would expect our model to break down in the moments preceding the collision because the squeezing and drainage of liquid out from between the bubbles significantly influences the bubble dynamics (see, for example, Crabtree & Bridgwater Reference Crabtree and Bridgwater1971; Chauhan & Kumar Reference Chauhan and Kumar2020; Ohashi, Toramaru & Namiki Reference Ohashi, Toramaru and Namiki2022). Furthermore, when the distance between the bubble interfaces is of the order of the gap height, we expect additional three-dimensional effects to become important.

Figure 12. (a) The relative normal velocity, $U_n$ , of the two bubbles as a function of the polar angle, $\phi$ , for a fixed $R=2$ and $\delta$ shown by the colour bar. The dotted and dashed curves show $U_n$ as a function of $\phi$ at $\delta =\delta _1(2)$ , and $\delta =\delta _2(2)$ , respectively (see § 5.2). (b) Schematic of two bubbles touching showing the definitions of $\boldsymbol {n}$ and $\phi$ .

6.1. Asymptotic expansions

In this section, we calculate the first-order corrections in $\epsilon$ , the bubble aspect ratio (2.3), to the shapes of a pair of bubbles, each of which undergoes deformations induced by the presence of the other. For simplicity, we consider two bubbles aligned in the direction of the flow with centres at positions $(0,0)$ and $(\sigma ,0)$ , respectively, in the $(x,y)$ -plane (see figure 13). At leading order, the bubbles are circles of radii $ R_1=1$ , and $R_2=R$ , with velocities $\mathcal {U}_1=U_1$ and $\mathcal {U}_2=U_2$ given by (3.16), respectively. We assume that the deformations occur faster than the time scale, $1/|U_1-U_2|$ for the relative motion of the two bubbles, which means we can treat the deformations as quasi-steady, with $\sigma$ assumed to be a known constant.

Figure 13. Schematic of the two-bubble deformation problem. The background flow is from left to right.

To find the corrections to the bubble shapes, we return to the dynamic boundary condition (2.1c ) and expand the curvatures and bubble pressures in powers of $\epsilon$ as

(6.1a) \begin{align} &\qquad \kappa _k \sim \frac {1}{R_k}+\epsilon \kappa _{k1}+\cdots , \end{align}

(6.1b) \begin{align} &p_k \sim 1+\frac {\pi \epsilon }{4R_k}+\epsilon ^2 p_{k2}+\cdots , \end{align}

for $k\in \{1,2\}$ . Note that, for completeness, one should also expand the complex potential, $w(z)$ , and the bubble velocities, $U_1$ and $U_2$ , as asymptotic series in powers of $\epsilon$ . However, to find the first-order shape correction we only need the leading-order solutions (3.9) and (3.16), and so for ease of notation we do not include an additional subscript 0 for these variables. We note that our analysis does not at present determine the first corrections to the bubble velocities due to the deformations.

6.2. Deformation of the rear bubble

For the first bubble, the dynamic boundary condition (2.1c ) at $O(\epsilon ^2)$ reads

(6.2) \begin{align} \kappa _{11} = \frac {4p_{12}}{\pi } +\frac {12\delta ^3\eta ^3}{\pi }\operatorname {Re}\left [z+W\left (\frac {1-a z}{z-a}\right )\right ] - \frac {4\delta ^2\eta ^2U_1^{2/3}}{\pi }\beta (\boldsymbol {i}\boldsymbol {\cdot }\boldsymbol {n})|\boldsymbol {i}\boldsymbol {\cdot }\boldsymbol {n}|^{2/3}, \end{align}

on $|z|=1$ . We define polar coordinates centred at $(0,0)$ , so the bubble surface is given by $r=1+\epsilon g_1(\theta )$ , where $\theta$ is the polar angle. The dynamic boundary condition (6.2) in polar coordinates is given by

(6.3) \begin{align} -g_1^{\prime\prime}\! - g_1 = \frac {4p_{12}}{\pi } +\frac {12\delta ^3\eta ^3}{\pi }\operatorname {Re}\left [\mathrm {e}^{\mathrm {i}\theta }+W\!\left (\frac {1-a \mathrm {e}^{\mathrm {i}\theta }}{\mathrm {e}^{\mathrm {i}\theta }-a}\right )\!\right ] - \frac {4\delta ^2\eta ^2U_1^{2/3}}{\pi }\beta (\cos \theta )|\cos \theta |^{2/3}. \end{align}

We determine $p_{12}$ by enforcing conservation of bubble area, i.e.

(6.4) \begin{align} \int _0^{2\pi }g_{1}\,\mathrm {d}\theta =-\int _0^{2\pi }\kappa _{11}\,\mathrm {d}\theta = 0. \end{align}

We solve (6.3) by expanding $g_1$ as the Fourier cosine series

(6.5) \begin{align} g_1(\theta ) = \frac {c_0}{2}+\sum _{n=1}^{\infty } c_n \cos n \theta . \end{align}

By the area conservation condition (6.4), we find that $c_0=0$ . We further fix the centroid of the bubble at the origin, which corresponds to $c_1=0$ . The remaining coefficients ( $n\geqslant 2)$ are determined by

(6.6) \begin{align} c_n = \frac {1}{(n^2-1)}\int _0^{2\pi }\frac {12\delta ^3\eta ^3}{\pi ^2}\operatorname {Re}\left [W\left (\frac {1-a \mathrm {e}^{\mathrm {i}\theta }}{\mathrm {e}^{\mathrm {i}\theta }-a}\right )\right ] \cos n \theta \,\mathrm {d}\theta - \frac {4\delta ^2\eta ^2U_1^{2/3}b_n}{\pi (n^2-1)}, \end{align}

where the $b_n$ are the Fourier coefficients of $\beta (\cos \theta )|\cos \theta |^{2/3}$ and are given by

(6.7) \begin{align} b_n &= \frac {\Gamma \left (\frac {5}{3}\right )\Gamma \left (\frac {n}{2}-\frac {1}{3}\right )}{4\pi 2^{2/3}\Gamma \left (\frac {n}{2}+\frac {4}{3}\right )}\left [\left ((\sqrt {3}+1)\beta _1+(\sqrt {3}-1)\beta _2\right )(-1)^{\left\lfloor \frac {n-1}{2}\right\rfloor }\right. \nonumber \\ &\quad -\left. \left ((\sqrt {3}-1)\beta _1+(\sqrt {3}+1)\beta _2\right )(-1)^{\left\lfloor \frac {n}{2}\right\rfloor }\right ]. \end{align}

Equation (6.5) then determines the first-order shape correction of the rear bubble $\partial \Omega _1$ .

6.3. Deformation of the front bubble

We proceed similarly with the second bubble. By defining polar coordinates centred at $(\sigma ,0)$ , we find that the bubble surface is given by $r=R+\epsilon g_2(\theta )$ , where $g_2(\theta )$ is given by the Fourier series

(6.8) \begin{align} g_2(\theta ) = \frac {d_0}{2}+\sum _{n=1}^{\infty } d_n \cos n \theta . \end{align}

By area conservation, we find that $d_0=0$ . We further fix the centroid of the bubble at $(\sigma ,0)$ , which corresponds to $d_1=0$ . The remaining coefficients ( $n\geqslant 2)$ are determined by

(6.9) \begin{equation} d_n = \frac {R^2}{(n^2-1)}\int _0^{2\pi }\frac {12\delta ^3\eta ^3}{\pi ^2}\operatorname {Re}\left [W\left (\frac {1-a (\sigma +R\mathrm {e}^{\mathrm {i}\theta })}{ (\sigma +R\mathrm {e}^{\mathrm {i}\theta })-a}\right )\right ] \cos n \theta \,\mathrm {d}\theta - \frac {4R^2\delta ^2\eta ^2U_1^{2/3}b_n}{\pi (n^2-1)}. \end{equation}

This then determines the first-order shape correction of the front bubble $\partial \Omega _2$ .

Figure 14. Experimental bubble shapes (black solid), asymptotic solution (6.5) and (6.8) (red dashed) dashed for $R=1$ , $\delta =2.86$ and $\sigma =$ (a) 2.68, (b) 2.56, (c) 2.43. The corresponding different dimensionless times $t = \hat {t}\hat {U}/\hat {R}_1$ are shown above for the experiments. The background flow is from left to right. Experimental images have been rescaled by the rear bubble radius, $\hat {R}_1 =$ 5.4 mm, for comparison with the theory. The bubble shapes from experiment and asymptotics are aligned so that the centroids of the bubble pairs coincide.

6.4. Results

6.4.1. Identical bubbles ( $R=1$ )

In figure 14, we show example solutions for the bubble shapes at different separations, $\sigma$ , calculated using (6.5) and (6.8), with $R=1$ , $\delta =2.86$ and $\epsilon =0.027$ , alongside experimental measurements under the same conditions. We observe good agreement between theory and experiments. The bubble in front flattens in the direction of motion (left to right), and the bubble behind elongates. In the theoretical plots, we use $\sigma$ as a proxy for time, because we assume that the deformations are quasi-steady.

Figure 15. The in-plane bubble aspect ratios, $A_{k}$ , versus separation, $\sigma$ , for the rear bubble ( $k=1$ , dashed curve and open markers) and the front bubble ( $k=2$ , solid curve and filled markers), with $\delta =2.86$ and $\epsilon =0.027$ . The points show experimental measurements and the curves are the asymptotic predictions (6.10). The different marker shapes (triangle, circle, diamond) represent distinct pairs of bubbles that were tracked and measured as the rear bubble caught up and collided with the front bubble. The error between experiment and theory is approximately $6$ %– $10\,\%$ .

To quantify our results, we define the in-plane bubble aspect ratios as

(6.10) \begin{equation} A_k = \frac {2R_k+\epsilon (g_k(0)+g_k(\pi ))}{2R_k+2\epsilon g_k(\pi /2)}\sim 1+\frac {\epsilon }{2R_k}(g_k(0)+g_k(\pi )-2g_k(\pi /2)), \end{equation}

for $k\in \{1,2\}$ . In figure 15, we plot $A_{1,2}$ versus bubble separation, $\sigma$ , for a fixed value of  $\delta$ . We observe that, as the bubbles become close, the disparity between their aspect ratios increases: the bubble in front becomes more flattened, while the rear bubble develops a more pronounced elongation. There is good agreement between the predicted and experimentally measured aspect ratio $A_2$ of the front bubble, however, there is a constant offset of approximately $0.06$ , which induces an approximate error between the theory and experiments. The model generally over-predicts the degree of flattening of the front bubble. For the aspect ratio $A_1$ of the rear bubble, there is a discrepancy between theory and experiments. In the experiments, $A_1$ is approximately constant, however, our model predicts this to be a monotonically decreasing function, and thus under-predicts the elongation of the rear bubble. In the experiments, the two bubbles become very close and, in this limit, we expect the theory may break down due to the three-dimensional effects in the fluid flow between the two bubbles. In addition, our dynamic boundary condition (2.1c ) is strictly valid only when the normal velocities at corresponding points on the front and rear menisci are equal and opposite; when the bubbles deform significantly this is no longer true and we should incorporate the full Burgess & Foster (Reference Burgess and Foster1990) boundary conditions on the bubble surface. Furthermore, bubbles in Hele-Shaw cells that are approaching or separating experience additional stresses due to their relative motion (Bremond, Thiam & Bibette Reference Bremond, Thiam and Bibette2008; Lai, Bremond & Stone Reference Lai, Bremond and Stone2009; Chan, Klaseboer & Manica Reference Chan, Klaseboer and Manica2010) that have not been included in our analysis.

In § 3 we found that, if $R=1$ , the bubbles travel at the same velocity at leading order in $\epsilon$ . However, in experiments, we observe that the bubbles approach each other while deforming, due to $O(\epsilon )$ corrections to the velocities which we currently do not calculate. Ultimately, the bubbles collide and coalesce when $\sigma \lt 1+R$ , a range that is inaccessible with our current analytical methods.

Wu et al. (Reference Wu, Booth, Griffiths, Howell, Nunes and Stone2024) found that, if an isolated bubble is flattened in the direction of motion, then the leading-order solution over-predicts the bubble velocity, and vice versa if the bubble is elongated. The same line of reasoning here would suggest that the velocity of the bubble at the front is over-predicted by (3.16), while the velocity of the bubble behind is under-predicted by (3.16). Thus, the bubble behind would travel faster than the bubble in front, resulting in the collision of bubbles of equal size. Similar behaviour has been observed experimentally and computationally for a pair of unconfined bubbles rising at low Reynolds numbers due to buoyancy (Manga & Stone Reference Manga and Stone1993, Reference Manga and Stone1995). As a result of the interaction between the bubbles, the leading bubble flattens in the direction of motion while the bubble behind elongates, and the distance between them decreases until they collide. Our results establish that there is an analogous mechanism for bubble collision in Hele-Shaw cells.

Figure 16. Experimental bubble shapes (black solid), asymptotic solution (6.5) and (6.8) (red dashed) for (a–c) $R=1.23$ , $\delta =2.55$ and $\sigma =$ (a) 2.39, (b) 2.34, (c) 2.28, (d–f) $R=1.65$ , $\delta =1.94$ and $\sigma =$ (d) 3.45, (e) 3.23, (f) 2.94. The corresponding different dimensionless times $t = \hat {t}\hat {U}/\hat {R}_1$ are shown above for the experiments. The background flow is from left to right. Experimental images have been rescaled by the rear bubble radii, $\hat {R}_1 =$ (a–c) 2.9 mm and (d–f) 4.8 mm, for comparison with the theory. The bubble shapes from experiment and asymptotics are aligned so that the centroids of the bubble pairs coincide.

6.4.2. Bubbles of different radii ( $R \neq 1$ )

In the absence of shape deformation, larger bubbles are expected to travel faster than smaller ones (Booth et al. Reference Booth, Griffiths and Howell2023). For this case, conditions were derived in § 5.2 under which a larger bubble can catch and collide with a smaller bubble in finite time. In § 6.4.1, we presented suggestions of a further mechanism arising from shape deformation by which two bubbles of equal size can collide. Here, we show that shape deformations and the resulting effects on the surrounding flow can be strong enough to enable a smaller bubble to catch a larger bubble.

Figure 17. The bubble aspect ratios, $A_{k}$ , versus separation, $\sigma$ , for the rear bubble ( $k=1$ , dashed curve and open markers) and the front bubble ( $k=2$ , solid curve and filled markers), with (a) $R=1.23$ , $\delta =2.55$ and $\epsilon =0.03$ (b), $R=1.65$ , $\delta =1.94$ and $\epsilon =0.05$ . The points show experimental measurements, and the curves are the asymptotic predictions (6.10). The error between experiment and theory is approximately (a) $5$ %– $7\,\%$ and (b) $10$ %– $13\,\%$ .

We show example solutions for the bubble shapes given by (6.5) and (6.8) alongside experimental images with, $R=1.23$ and $\delta =2.55$ in figures 16(a)–16(c) and $R=1.65$ and $\delta =1.94$ in figures 16(d)–16(f). Similarly to the examples of the bubbles with the same leading-order radius (see figure 14), the leading bubble flattens in the direction of motion, whereas the rear bubble elongates. To quantify this observation, we plot the bubble aspect ratios $A_{1,2}$ versus separation, $\sigma$ , in figure 17. We observe good agreement between theory and experiments. In particular, we correctly predict that $A_1\gt A_2$ . Again, there is a discrepancy between the experimentally measured aspect ratios and theoretical predictions, which we attribute to the same reasons as discussed in § 6.4.1. Nevertheless, these results hint that, although the smaller rear bubble is expected to lag behind the larger front bubble when they are both circular, deformations may allow for a region of parameter space in which a smaller bubble can catch up to a larger one. Several collisions of this type have been observed experimentally, and the progression of shape deformation for a few examples is shown in figure 16. To establish this result theoretically, one would need to find the perturbation to the bubble speeds, for example by performing a complex variable analysis similar to that done by Wu et al. (Reference Wu, Booth, Griffiths, Howell, Nunes and Stone2024). We leave such analysis for future work.