3.1. General features

The results to be discussed below are obtained by solving (2.1) with the open source flow solver Basilisk developed by Popinet (Popinet Reference Popinet2009, Reference Popinet2015). Basilisk (see basilisk.fr) is the successor of Gerris (http://gfs.sourceforge.net) which has been widely employed over the last fifteen years in detailed explorations of interfacial flows. It makes use of Cartesian grids with a collocated discretization of the velocity and scalar fields. The temporal discretization is based on a second-order fractional step method. In particular, the Godunov-type unsplit upwind scheme developed by Bell, Colella & Glaz (Reference Bell, Colella and Glaz1989) is used to discretize the advection term, and a fully implicit scheme is used to compute the viscous term. A second-order projection method is employed to ensure that the computed velocity field is divergence free at the end of each time step. Interfaces are tracked and geometrically reconstructed by a VOF approach in which an accurate well-balanced height function method is used to calculate the interface curvature (Popinet Reference Popinet2009, Reference Popinet2018). An AMR technique makes it possible to locally refine the grid close to interfaces and high-vorticity regions, based on a wavelet decomposition of the gas volume fraction and velocity fields, respectively (Van Hooft et al. Reference Van Hooft, Popinet, van Heerwaarden, van der Linden, de Roode and van de Wiel2018). This strategy greatly enhances the computational efficiency while guaranteeing a high numerical accuracy in flow regions where subtle physical phenomena are likely to take place. In the present study, the spatial resolution is refined down to $\varDelta _{min} = R/68$ close to the bubble interface. Hence, at the highest Reynolds number considered here ($Re\approx 120$), approximately $6$ grid points lie within the boundary layer whose thickness is estimated to $\delta _b \sim Re^{-1/2}$. Figure 2 shows how a typical grid is refined in the vicinity of the two bubbles. In addition to the near-interface zones where $\varDelta =\varDelta _{min}=R/68$, refined regions include the boundary layer and wake of each bubble, where the local cell size is $\varDelta =R/17$. The grid coarsens drastically beyond the region displayed in the figure, the largest cells in the far field corresponding to $\varDelta \approx 7.5R$. Hence the ratio between the largest and smallest cells in the whole domain is $2^9=512$. Previous works have established the capability of Basilisk to accurately simulate the dynamics of isolated rising bubbles. Moreover, since Basilisk succeeded Gerris and essentially makes use of the same algorithms, the numerous validations performed with Gerris over the years also hold for Basilisk. For instance, Popinet (Reference Popinet2017) reproduced successfully with Basilisk the results obtained with Gerris by Cano-Lozano et al. (Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016a) concerning the transition from rectilinear to zigzagging or spiralling paths of air bubbles rising in various liquids. In § 5 we shall compare present predictions for the equilibrium distance of two nearly spherical bubbles ($Bo\ll 1$) rising in line with the findings of Yuan & Prosperetti (Reference Yuan and Prosperetti1994) and Hallez & Legendre (Reference Hallez and Legendre2011). This comparison will establish the relevance and accuracy of the present approach in the axisymmetric case. Nevertheless, this configuration will be proved to be unstable in § 6. Therefore, it is necessary to determine to which extent the fate of the three-dimensional system depends on numerical details driving the onset of this instability. This is the purpose of Appendix A in which several specific tests are reported.

Figure 2. Detail of a typical grid for a bubble pair with $Ga=30$ and $Bo=0.3$ at $t=8$. $(a)$ Grid refinement in the wake and boundary layer regions; four refinement levels are shown, from dark blue ($\varDelta =R/8.5$) to yellow ($\varDelta =\varDelta _{min}=R/68$). $(b)$ Grid detail within the leading bubble and in the neighbourhood of its surface. The gas–liquid interface is marked with a black (respectively red) line in $(a)$ (respectively $b$).

Most of the computations were run on a personal computer with 24 Intel® Xeon® E5-2630 v3 processors. The Intel-MPI library was used to exchange information between the processors. A typical run extending over $50$ time units took approximately $50$ days (e.g. the case $Ga=30,Bo=0.3$ discussed in § 6.3). Note that, since capillary effects impose a specific time step constraint, low-$Bo$ cases require longer computational times. For instance, to reach a given physical time, the case $Ga=30, Bo=0.05$ was $1.5$ more time consuming than the previous case. Due to the AMR procedure, the grid evolves over time. For $Ga=30, Bo=0.3$, its size stabilizes at approximately $2.5$ million cells during the second half of the run. When additional levels of grid refinement are introduced because the bubbles get very close to each other (see below), this size increases significantly. In such cases, e.g. $Ga=20,Bo=0.5$ discussed in § 6.4, the complete grid involves up to $4.2$ million cells.

3.2. Treatment of thin films: numerical vs physical coalescence

In the present problem it is likely that, under certain conditions, the two bubbles get very close to each other. When this happens in a real flow, coalescence may or may not take place, depending on the mobility of the interfaces involved and on the strength and duration of the forces that drive the two bubbles toward each other (Chesters Reference Chesters1991; Chan, Klaseboer & Manica Reference Chan, Klaseboer and Manica2011). Dealing numerically with such situations is particularly challenging, owing to the very small scales involved. Some numerical approaches, especially the front tracking technique (Unverdi & Tryggvason Reference Unverdi and Tryggvason1992; Tryggvason et al. Reference Tryggvason, Bunner, Esmaaeli, Juric, Al-Rawahi, Tauber, Han, Nas and Jan2001), totally prevent coalescence. The same may be achieved in VOF approaches by identifying the two bubbles with separate markers, each of them representing the local volume fraction of the corresponding body. However, this numerical option is not fully appropriate here. Indeed, for the reasons mentioned below, bubbles rising in line in a pure liquid offer one of the physical situations with the highest coalescence probability. If coalescence takes place under real conditions, it is of course desirable to track numerically the post-coalescence dynamics, i.e. the shape and path evolution of the resulting bubble, which the above option would not allow.

That bubbles rising in line in a pure liquid are prone to coalesce in a number of cases is due to the combination of two factors. First, as will be discussed in the next section, the wake of the leading bubble provides a permanent attractive force to the trailing bubble. Under a number of flow conditions, this force is strong enough to make the two bubbles come virtually in contact in the head-on configuration. Second, the mobility of the gas–liquid interfaces when the bubbles are free of any contamination makes the drainage of the interstitial film several orders of magnitude faster than that of liquid–liquid or contaminated gas–liquid interfaces (Vakarelski et al. Reference Vakarelski, Marica, Li, Basheva, Chan and Thoroddsen2018). The coalescence process may take several distinct forms, depending on the fluid characteristics and bubble size. After the two bubbles collide, the interstitial film may rupture quickly, or the bubbles may stay almost in contact during a long time before eventually coalescing. If viscous effects are small enough, bubbles larger than a critical size bounce after the film has been drained partially and coalesce only after one or several bounces. Energetic considerations and detailed experimental data may be employed to determine which of these scenarios takes place. This issue, together with the specificities of the coalescence of ‘clean’ bubbles are discussed in Appendix B.

In situations potentially leading to coalescence, one would ideally like to track numerically all steps of the drainage, until non-hydrodynamic effects such as the London–van der Waals force come into play and rupture the film. This typically occurs when the minimum gap between the two interfaces is of the order of $10$ nm. For millimetre-size bubbles, this would require approximately ten additional grid levels beyond the one corresponding to $\varDelta _{min}$. Since the largest scales to be resolved in the present context are of the order of $10$ cm, i.e. seven orders of magnitude larger, the corresponding computational cost would be prohibitive. At least, it is possible to track the first stages of the drainage based on a suitable grid refinement technique. Then, referring to the available knowledge summarized in Appendix B, one can reasonably predict which coalescence scenario is taking place in the real system, although the grid resolution does not allow all its details to be captured. The main shortcoming of this approach is that bubbles usually merge too early in the computations, i.e. they would merge at a higher vertical position in a real flow. Nevertheless, as we discuss in Appendix B, the theoretical and experimental knowledge available on the coalescence of nearly spherical clean bubbles allows the corresponding temporal shift to be estimated in a number of cases.

To track the first steps of the drainage, we developed a specific topology-based AMR scheme to refine the grid within the gap (Zhang et al. Reference Zhang, Chen and Ni2019). The corresponding algorithm checks whether or not any cell crossed by the gas–liquid interface has at least one neighbouring cell filled with only liquid or gas. If not, the cell crossed by the interface is automatically refined. This adaptive strategy is illustrated in figure 3. At $t=t_1$ (panel $a$), the two interfaces are separated by a ‘pure’ liquid cell, but this is no longer the case in $(b)$ where the cells standing in the yellow region all contain a non-zero gas fraction. Therefore, these cells are refined by a factor of 2 in each direction, as shown in $(c)$. Ultimately, for the reasons discussed above, we let the two interfaces merge numerically if the number $N$ of successive refinements required to satisfy the above ‘pure liquid neighbouring cell’ criterion exceeds a prescribed value. In practice, this is achieved by merging the two markers which identify each bubble in all previous steps. The three-dimensional results to be discussed later were all obtained with $N=2$, so that the minimum cell size in the gap was $\delta _{min}=\varDelta _{min}/2^N=R/272$. For a $0.25$ mm radius bubble, this implies $\delta _{min}\approx 1\ \mathrm {\mu }$m, which is still two orders of magnitude larger than the typical thickness at which the interstitial film actually ruptures.

Figure 3. Topology-based AMR strategy employed to refine the grid from time $t=t_1$ in $(a)$ to $t=t_1+{\rm \Delta} t$ in $(b,c)$, as the liquid film separating the two bubbles gets squeezed.

6.1. Overview of the results

As mentioned earlier, it has been known since Harper (Reference Harper1970) that the in-line configuration of two clean spherical bubbles is unstable with respect to side disturbances when the Reynolds number is large. Experimental investigations carried out under surfactant-free conditions (silicone oils) in the range $10\lesssim Re\lesssim 150$ confirm this prediction (Kusuno & Sanada Reference Kusuno and Sanada2015; Kusuno et al. Reference Kusuno, Yamamoto and Sanada2019). We performed fully three-dimensional time-dependent simulations to assess this stability issue and explore its consequences. Similar to the axisymmetric case, the gas/liquid density and viscosity ratios were set to $10^{-3}$ and $10^{-2}$, respectively. In agreement with the above experimental and theoretical findings, we observed that the bubble pair never maintains a straight vertical path except when the Bond number is large enough for coalescence to eventually happen. We actually identified three drastically different evolutions, depending on the value of the Galilei and Bond numbers. Beyond a critical Bond number, $Bo_c(Ga)$ which increases approximately from $0.2$ for $Ga=10$ to $0.5$ for $Ga=30$, the two bubbles collide and eventually coalesce, this ‘coalescence’ having to be interpreted in the light of the discussion of § 3.2 (see below). Predictions of the three-dimensional and axisymmetric simulations superimpose for $Bo\gtrsim 0.5$ when $Ga\leq 20$ (and for $Bo\gtrsim 1$ when $Ga=30$), indicating that the in-line configuration is stable with respect to azimuthal disturbances for sufficiently deformed bubbles. In contrast, for $Bo< Bo_c(Ga)$, the TB escapes from the wake of the LB at some point, and the two go on rising with their line of centres more or less inclined with respect to the vertical and their centroids widely separated. In such cases, we found that the interplay of the two bubbles after the three-dimensional effects set in may follow two markedly different scenarios. One is clearly a DKT-type mechanism. In this case, both bubbles deviate from their initial trajectory and eventually rise almost side by side along two straight vertical lines distinct from the initial path. In contrast, in the other scenario, which we refer to as asymmetric side escape (hereinafter abbreviated as ASE), the lateral drift of the TB leaves the path of the LB almost unaffected. Hence, this bubble goes on rising virtually along its initial path, while after the system has reorganized itself, the TB follows a markedly distinct vertical path. The structural differences between the configurations corresponding to the DKT and ASE scenarios are well visible in the recent observations of Kusuno et al. (Reference Kusuno, Yamamoto and Sanada2019); see especially their figures 4 and 5.

Figure 6 displays the phase diagrams and typical paths corresponding to the above three evolutions, still for $\bar {S}_0=8$. The influence of initial conditions, i.e. angular inclination and separation, on the borders of the different subdomains will be discussed in §§ 7.3 and 7.4, respectively. The influence of numerical parameters, among which the grid resolution, on the coalescence threshold, i.e. on $Bo_c(Ga)$, is estimated in Appendix A. These technical aspects are found to barely change $Bo_c(Ga)$ by a few per cent.

Figure 6. Set of configurations encountered in the three-dimensional simulations: $(a)$ $(Ga, Bo)$ phase diagram; $(b)$ $(Re, \chi )$ phase diagram (for each $(Ga,Bo)$ pair, $Re$ and $\chi$ are the steady-state values determined with the corresponding isolated bubble); $(c)$ typical trajectories illustrating each of the three configurations. Bullets: DKT scenario; solid squares: ASE scenario; solid triangles: collision followed by coalescence. In $(c)$, the three images from left to right correspond to $(Ga, Bo) = (10, 0.1)$, $(30, 0.3)$ and $(20, 0.5)$, respectively. The thin solid lines in $(a)$ are the iso-$Mo$ lines corresponding to different liquids, with water at the very bottom, then silicone oils T0-T11 of increasing viscosity from bottom to top (see e.g. Zenit & Magnaudet (Reference Zenit and Magnaudet2008) for the corresponding physical properties). The dashed line in $(a,b)$ separates the subdomain of weakly deformed bubbles in which a finite equilibrium separation is reached in the axisymmetric case, from that in which the bubbles eventually coalesce.

As the thin black lines indicate, coalescence only takes place in the present configuration in liquids with a sufficiently high Morton number, from $Mo\approx 8\times 10^{-7}$ for $Ga=10$ to $Mo\approx 1.5\times 10^{-7}$ for $Ga=30$. This means in particular that bubbles rising in pure water ($Mo\approx 2.6\times 10^{-11}$) never coalesce in the $Ga$-range considered here. Interestingly, figure 6$(a),(b)$ also reveals that, under some conditions, such as $(Ga,Bo)=(20,0.2)$ or $(30,0.2)$, bubbles that would coalesce if the system were constrained to remain axisymmetric actually escape from coalescence through the ASE scenario. From figure 6$(b)$, it may also be concluded that for $Re\approx 15$, even an $8\,\%$ departure from sphericity is sufficient to lead to coalescence. Actually, the discussion of Appendix B indicates that, for such modest Reynolds numbers, viscous effects tremendously delay the drainage of the interstitial film. This is why experiments performed in this regime (Sanada, Watanabe & Fukano Reference Sanada, Watanabe and Fukano2006; Watanabe & Sanada Reference Watanabe and Sanada2006) reveal that, after the two bubbles collide, they merely stay ‘glued’ to each other and rise as a single ‘dumbbell’ bubble without coalescing during the time window of the observations. However, since all forces involved in the physical system, including the London–van der Waals force, are attractive, there is no doubt that such bubbles eventually coalesce. Unfortunately, the spatial resolution of the present simulations is clearly insufficient to reproduce this slow coalescence process. That is, the computations correctly predict the final state of the physical system but this state is reached too early. This underestimate of the coalescence time still exists at larger Reynolds number (see § 6.4) but, for reasons discussed in Appendix B, reduces as $Re$ increases. In figure 6$(b)$, the maximum deformation below which the two bubbles do not coalesce is seen to increase significantly with $Re$. For instance, bubbles with $\chi \approx 1.5$ follow an ASE scenario for $Re\approx 70$, but coalesce for $Re\approx 35$ (in these estimates, $\chi$ and $Re$ are evaluated from the steady-state properties of the isolated bubble corresponding to the same $(Ga,Bo)$ set).

Figure 7 depicts the bubble shapes and relative positions during the lateral escape of the TB or the coalescence process. Considering the DKT (green) and ASE (red) regimes at a fixed $Ga$, the figure indicates that the larger $Bo$ the shorter the separation during the escape stage. This may be interpreted as a stabilizing effect of the deformation, since the bubble pair is able to maintain a vertical path during a longer time (i.e. down to a shorter separation) when the Bond number increases. Then, coalescence takes place when $Bo$ exceeds the critical threshold $Bo_c(Ga)$ which is seen to increase significantly with $Ga$, as already noticed in figure 6$(a)$. In most cases, coalescence is reached through a head-on approach, i.e. without any prior escape of the TB, the axisymmetric film in the gap being gradually squeezed. However, for $Ga=30$ and $Bo=0.5$, an intermediate configuration corresponding to an oblique approach is noticed. These different regimes are examined in more detail in the rest of this section.

Figure 7. Snapshots of the bubble shapes and relative positions observed during the lateral escape of the TB or the pre-coalescence process. Red, green and blue pairs correspond to the ASE, DKT and coalescence scenarios, respectively. In the red and green regions, the snapshots are taken by the time the horizontal distance between the two centroids is approximately equal to one bubble initial radius; in the blue region, each snapshot is the last in a run before numerical coalescence occurs. Since successive snapshots are separated by a finite time interval, the remaining time until coalescence may differ among the various blue pairs.

6.2. Drafting–kissing–tumbling

The DKT scenario followed by sedimenting particle pairs has been widely described for spheres (Joseph et al. Reference Joseph, Fortes, Lundgren and Singh1986; Fortes et al. Reference Fortes, Joseph and Lundgren1987; Feng, Hu & Joseph Reference Feng, Hu and Joseph1994), thick disks (Brosse & Ern Reference Brosse and Ern2014) and, under certain conditions, prolate spheroids (Ardekani et al. Reference Ardekani, Costa, Breugem and Brandt2016). In short, owing to the sheltering effect induced by the wake of the leading body, the trailing body first catches up with it (drafting) until the two collide (kissing). Then the resulting prolate compound body becomes unstable to transverse disturbances, which makes it rotate (tumbling) in such a way that the line joining the centroids of the two individual bodies tends to become horizontal, letting them eventually fall/rise separately in a side-by-side configuration. A similar behaviour of bubbles rising in line has been reported experimentally by Kusuno et al. (Reference Kusuno, Yamamoto and Sanada2019) in the range $10\lesssim Re\lesssim 25$, $0.3\lesssim We\lesssim 1.1$. Here, we identified it for bubble pairs corresponding to $Ga\lesssim 12$ and $Bo\lesssim 0.2$. With $(Ga, Bo)=(10, 0.02)$, i.e. $Mo=8\times 10^{-10}$, the final state corresponds to $Re\approx 16$ and $\chi \approx 1.02$, which indicates that the bubbles shape remains close to a sphere. The corresponding evolution is displayed in figure 8.

Figure 8. Path of a bubble pair following a DKT scenario ($Ga=10, Bo=0.02$). $(a)$ Side and top views, with $r$ the radial distance to the initial path and numbers referring to the dimensionless time at the corresponding position; $(b)$ three-dimensional streamlines past the bubbles at different instants of time, in the reference frame of the LB. A zoom of the flow field in the gap is provided for $t=30$ and $t=38$, to highlight the inception of the non-axisymmetric fluid motion.

As revealed in panel $(a)$, the tumbling process starts at $t\gtrsim 37$ when the two bubbles reach the same rise velocity (see figure 9$a$). After this process is completed (say at $t=54$ in figure 8$a$), both paths have experienced large lateral deviations. Their lateral separation is approximately $2.3$, the LB having been slightly more deviated. They stand virtually side by side and keep on rising in this configuration, although the rise velocity of the (formerly) TB is barely larger than that of the LB, inducing a tiny difference in the final altitude reached by the two bubbles. During the side-by-side rise, the lateral separation $\bar {S}_r$ evolves significantly, until an equilibrium value $\bar {S}_{re}\simeq 2.55$ is reached. This is consistent with the findings of Legendre et al. (Reference Legendre, Magnaudet and Mougin2003), which indicate that, in this configuration, the transverse force acting on a pair of spherical bubbles almost vanishes for $\bar {S}_r=2.5$ when $Re=25$ (their figure 13). For $Y>175$, the centre of inertia of the bubble tandem is seen to drift slightly towards the left, a consequence of the tiny inclination of its line of centres. The sign of this drift is consistent with the conclusions of Hallez & Legendre (Reference Hallez and Legendre2011) who found (their figure 9$b$) that in the range $10\leq Re\leq 25$, the centre of inertia drifts laterally toward the position of the higher bubble, provided $40^\circ \lesssim \theta <90^\circ$. The horizontal trace of the TB and LB paths in figure 8$(a)$ indicates that the entire motion remains planar and takes place within a vertical plane. Hence, by breaking the initial axial symmetry of the flow, the DKT mechanism changes the initial one-dimensional path into a two-dimensional path, but the latter keeps track of the former since the preferential direction of the bubble motion remains unchanged.

Figure 9. Evolution of several characteristics of a bubble pair with $(Ga, Bo)=(10, 0.02)$ undergoing a DKT interaction. $(a)$ Vertical velocity component of the LB (solid red line) and TB (dashed blue line); $(b)$ same for the horizontal component; $(c)$ vertical (red line, left axis) and horizontal (blue line, right axis) separations; $(d)$ inclination of the line of centres with respect to the vertical. Open squares and circles in (a,c) refer to the axisymmetric prediction.

Some three-dimensional streamlines past the two bubbles are displayed in figure 8$(b)$. No standing eddy is observed at the back of the bubbles. This is in line with the generating mechanism of vorticity on a curved shear-free surface discussed in §  4, which maintains the flow past a shear-free sphere unseparated whatever the Reynolds number (Blanco & Magnaudet Reference Blanco and Magnaudet1995; Magnaudet & Mougin Reference Magnaudet and Mougin2007). At $t=38$, the separation has almost reached its minimum ($\bar {S}\approx 2.3$) and the flow is still almost axisymmetric, except within the gap where a small left–right asymmetry is discernible in the enlarged view. In the present case, the origin of the axial symmetry breaking stands in tiny numerical asymmetries, in the first place those resulting from the pressure field returned by the multilevel Poisson solver (see Appendix A). Tumbling then starts. At $t=44$ and $t=48$, the flow field exhibits strong left–right asymmetries, which result in transverse forces that act to move the two bubbles in opposite directions. Later, say for $t\geq 54$, the flow field gradually rearranges towards a left–right symmetric configuration corresponding to a side-by-side motion of the two bubbles, with an equilibrium horizontal separation $\bar {S}_{re}\approx 2.55$.

The evolution of the rise speed of the two bubbles is plotted in figure 9$(a)$. The three-dimensional prediction is seen to coincide with its axisymmetric counterpart up to $t\approx 45$, i.e. throughout the time period during which the bubble pair moves in straight line and even during the early stage of the tumbling process. That the three-dimensional and axisymmetric predictions virtually superimpose up to $t=37$ proves the reliability of the former. Note that the red and blue curves superimpose up to $t=4$, which corresponds to the very early stage during which the two bubbles rise independently. During the next stage, say $4< t<20$, the TB goes on accelerating while the rise speed of the LB (hereinafter denoted as $V_{LB}$) stays almost constant. This corresponds to the early stage of the interaction during which the TB is sucked in the wake of the LB while the latter remains unaffected. Then, from $t=20$ to $t=37$, the two bubbles get close enough for the TB to modify the wake of the LB, the rise speed of which increases sharply until the two bubbles rise with the same speed.

Tumbling starts at $t\approx 38$, without any real ‘kiss’ since figure 9$(c)$ indicates that $\bar {S}$ never fell below $2.3$ at previous times. Up to $t=46$, no change is noticed in the rise speed of the two bubbles, nor in their vertical separation, although the tumbling process is going on, as the sharp rise of their lateral velocities (figure 9$b$) and that of the inclination of their line of centres (figure 9$d$) confirm. The situation significantly changes within the next short stage $46\leq t\leq 50$, during which $V_{LB}$ and $V_{TB}$ (the rise speed of the TB) drop sharply and the line of centres rotates by more than $30^\circ$. Not surprisingly, this rotation forces the lateral separation to increase dramatically, from $\bar {S}_{r}\approx 1$ at $t=46$ to $\bar {S}_{r}\approx 2.5$ at $t=50$. Conversely, the vertical separation is reduced to $\bar {S}\approx 2$ at $t=50$. The tumbling motion is also responsible for the slowing down of the two bubbles, as it makes the flow around them fully three-dimensional, which enhances the dissipation in the liquid, hence the drag on each bubble.

Tumbling goes on more slowly until $t\approx 85$, when the side-by-side configuration ($\theta =90^\circ$) and the equilibrium horizontal separation ($\bar {S}_{re}\approx 2.55$) are reached. Ultimately, $\theta$ slightly exceeds $90^\circ$, in line with the tiny difference between the final positions of the two bubbles noticed in figure 8$(a)$.

6.3. Asymmetric side escape

Computations carried out with somewhat larger values of the Galilei number ($15\leq Ga\leq 30$) but still fairly low Bond numbers reveal a drastically different evolution of the bubble pair. In what follows, we select the case $(Ga, Bo)=(30, 0.3)$, i.e. $Mo=3.3\times 10^{-8}$, as typical of this regime to discuss the corresponding dynamics. With these parameters, an isolated bubble follows a rectilinear path (Cano-Lozano et al. Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016a,Reference Cano-Lozano, Tchoufag, Magnaudet and Martinez-Bazanb) and its characteristic rise Reynolds number and aspect ratio at steady state are $Re=66$ and $\chi =1.62$, respectively. The corresponding paths, together with the streamlines past the two bubbles, are shown in figure 10. No standing eddy is observed at the back of the bubbles, in line with the flow structure past an isolated bubble with the same parameters (Blanco & Magnaudet Reference Blanco and Magnaudet1995). In panel $(b)$, only streamlines emanating from the right half-plane ahead of the LB are shown for $10\leq t\leq 16$ in order to help identify the onset of the non-axisymmetric motion. While all streamlines get around the TB within the right half-plane at $t=10$, one of them deviates to the left half-plane at $t=11$, indicating that axial symmetry has just broken. At this stage, the separation between the two bubbles is still large ($\bar {S}\approx 6.8$), in contrast with the situation observed in the DKT scenario. According to panel $(a)$, the TB also departs from its original vertical path at $t\approx 11$, but the LB is left virtually unaffected by this departure until $t\approx 15$. At this point, the TB still stands in the wake of the LB but its equatorial plane has tilted clockwise, which distorts the flow in the gap and makes it significantly asymmetric at the back of the LB. This asymmetry induces a slight anticlockwise tilt of the LB path until $t\approx 20$, whereas the TB goes on drifting laterally and escapes completely from the wake of the LB. Then the TB rotates anticlockwise in such a way that its equatorial plane becomes again horizontal, and its drift stops at $t\approx 23$, when the lateral separation between the two bubbles is approximately $\bar {S}_r\approx 4$. Due to this temporary anticlockwise rotation, the TB slightly drifts back until $t\approx 28$. The motion of the bubble pair has remained planar up to this point. However, at $t\approx 28$, they both start to drift slightly out of their previous plane of rise. Then, each bubble goes on rising along a straight but slightly inclined path. Both paths being tilted in the same direction but the angle being larger in the case of the LB, the lateral separation weakly but consistently increases over time. As the TB spent part of its potential energy to move out from the wake of the LB during the time period $11\lesssim t\lesssim 23$, it never catches up, and the LB remains significantly ahead of the TB in the final configuration.

Figure 10. Path of a bubble pair following an ASE scenario ($Ga=30, Bo=0.3$). $(a)$ Side and top views of the path with numbers referring to the dimensionless time at the corresponding position of the TB; $(b)$ three-dimensional streamlines past the bubbles at different instants of time, in the reference frame of the LB.

Figure 11 describes the evolution of the same four characteristics as in figure 9 for the above bubble pair. According to the axisymmetric prediction reported in panel $(c)$, the vertical separation decreases monotonically and the two bubbles enter the initial stage of coalescence at $t\approx 23$. However, the three-dimensional evolution reveals a totally different evolution beyond $t\approx 11$, although the rise speed of both bubbles and their vertical separation do not exhibit discernible differences with the axisymmetric results up to $t\approx 16$. At this time, $V_{TB}$ starts to drop sharply and equals $V_{LB}$ at $t\approx 18$. In the meantime, the horizontal velocity of the TB has grown tremendously. Its maximum value, approximately $40\,\%$ of $V_{TB}$, is reached at $t\approx 19$. The horizontal velocity of the LB reaches its maximum at the same time. However, the ratio of the two maxima is less than $4\,\%$, which confirms that the LB is only marginally disturbed by the lateral drift of the TB. This situation dramatically differs from that depicted in figure 9$(b)$, where the two $V_r$ maxima have almost the same magnitude. Moreover, in the present ASE scenario, the maximum of $V_r$ for the TB is typically three times larger than the maxima encountered during the DKT-type interaction. As figure 11$(c)$ indicates, the horizontal separation has already grown up to $\bar {S}_r\approx 2$ by the time $V_r$ reaches its maximum (i.e. the TB has left the wake of the LB), and doubles during the next four time units at the end of which the horizontal velocity of the TB vanishes ($t\approx 23$). From $t\approx 18$ to $t\approx 23$, $V_{TB}$ has dropped below $V_{LB}$, which makes the vertical separation again increase, from $\bar {S}\approx 4.5$ at $t=17.5$ to $\bar {S}\approx 5.3$ at $t=23$. The $25\,\%$ drop of $V_{TB}$ from $t\approx 16$ to $t\approx 19.5$ emphasizes the fact that a substantial fraction of the kinetic energy of the liquid displaced by the TB has been spent in the meantime to move it laterally and balance the rate of work of the corresponding sideways force.

Figure 11. Evolution of several characteristics of a bubble pair with $(Ga, Bo)=(30, 0.3)$ undergoing an ASE interaction. $(a)$ Vertical velocity component of the LB (solid red line) and TB (dashed blue line); $(b)$ same for the horizontal component; $(c)$ vertical (red line, left axis) and horizontal (blue line, right axis) separations; $(d)$ inclination of the line of centres with respect to the vertical. Open squares and circles in (a,c) refer to the axisymmetric prediction.

In the next stage ($23\leq t\leq 28$), the horizontal velocity of the TB changes sign and reaches a minimum of approximately $-0.06V_{TB}$. This negative $V_r$ results in the reversed lateral drift already discussed in connection with figure 10. The horizontal velocity still describes some damped oscillations before the vertical and horizontal separations stabilize and reach values close to $5.6$ and $4$, respectively. Meanwhile, the inclination angle of the bubble pair stabilizes at $\theta \approx 36^\circ$. As already pointed out, this configuration is not entirely steady as the slight non-zero slopes of the curves in the right part of figure 11$(c)$ confirm. In other words, the vertical and transverse components of the interaction force driving the relative position of the two bubbles have not completely vanished yet. However, the remaining values of these components are very small, so that it takes an extremely long time for the system to reach a true steady state. This slow final evolution may be connected to the findings of Hallez & Legendre (Reference Hallez and Legendre2011) for spherical bubbles. Their figure 9$(c)$ indicates that, provided $Re>25$, the two bubbles repel each other whatever their separation distance when the inclination of their line of centres is less than a critical angle $\theta _{c}\approx 53^\circ$. However, for $30^\circ <\theta <\theta _c$, this repulsive separation-dependent effect is very weak for $\bar {S}\gtrsim 5.0$, which corresponds to the present situation. That the equilibrium angle is close to $37^\circ$ for $Ga=30,Bo=0.3$ instead of the above value for spherical bubbles is likely an effect of the significant oblateness of the bubbles considered here.

6.4. Head-on collision and coalescence

As figures 6 and 7 revealed, increasing the Bond number beyond the $Ga$-dependent threshold $Bo_c(Ga)$ leads to the numerical coalescence of the two bubbles. Most of the time, this coalescence is initiated by a head-on approach, but in some cases the bubbles may also approach each other in an asymmetric manner; e.g. the case $(Ga, Bo)=(30, 0.5)$ in figure 7, which corresponds to near-threshold conditions. Actually, the approach configuration depends on the time elapsed since the bubble pair was released. As discussed in § 4, the attractive effect of the LB wake increases with the bubble deformation, owing to the direct relation between the interface curvature and the strength of the surface vorticity. This is why, for a given $Ga$, the time at which the two bubbles collide decreases as $Bo$ increases. This leaves less time for non-axisymmetric disturbances to grow before the collision if $Bo$ is significantly larger than $Bo_c(Ga)$, favouring the head-on configuration. In the case of strongly deformed bubbles, the lift reversal mechanisms discussed in Appendix C also act to increase the stability of the in-line configuration. Here, we select conditions $(Ga, Bo)=(20, 0.5)$, i.e. $Mo=7.8\times 10^{-7}$, as an archetype of the head-on coalescence scenario observed in the moderate-Reynolds-number range ($Re\approx 35$). As expected from the previous discussion, figure 12 shows that the two bubbles rise almost in a straight line before they numerically coalesce (a tiny lateral deviation of the LB actually takes place in the late stage of the approach, see figure 13$b$). As far as the bubbles evolve independently, their deformation is quite large, characterized by aspect ratios $\chi \approx 1.55$ and $\chi \approx 1.48$ for the LB and TB, respectively. In both cases, the aspect ratio is smaller than the critical value $\chi _c=1.65$ beyond which a standing eddy develops at the back of an isolated bubble (Blanco & Magnaudet Reference Blanco and Magnaudet1995). This is why no such structure is present in figure 12$(a)$. Beyond $t\approx 15$, the shape of the two bubbles changes significantly. On the one hand, the front part of the LB flattens, due to the proximity of the TB which makes its rise speed increase (see figure 13$a$). On the other hand, the front part of the TB becomes more rounded, owing to the suction induced by the wake of the LB. Finally, the TB catches up with the LB, leaving only a thin liquid film in the gap and making the two bubbles behave essentially as a bluff compound body. After coalescence takes place (see below), figure 13$(a)$ indicates that the resulting bubble first undergoes a series of large-amplitude oscillations until it relaxes to an oblate shape with an aspect ratio close to $2.0$ and a significant fore–aft asymmetry. Due to volume conservation, the radius of this bubble is $2^{1/3}\approx 1.26$ times that of the initial bubbles, so that its characteristic parameters are $Ga=2^{1/2}\times 20\approx 28.3$ and $Bo=2^{2/3}\times 0.5\approx 0.8$, respectively. Under such conditions, the computations of Cano-Lozano et al. (Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016a) predict that an isolated bubble still rises vertically, although it is close to the transition to a non-straight path. Present observations are in line with these earlier conclusions, since figure 13$(a)$ indicates that the small horizontal velocity component present at the time of coalescence subsequently decreases over time.

Figure 12. Pre-coalescence dynamics of a bubble pair with $Ga=20, Bo=0.5$; the streamlines in the cross-sectional plane are defined in the reference frame of the LB.

Figure 13. Evolution of some characteristics of a bubble pair with $Ga=20, Bo=0.5$ undergoing coalescence. $(a)$ Vertical (solid lines, left axis) and horizontal (dash-dotted lines, right axis) components of the bubble velocity, with the red, blue and black lines corresponding to the LB, TB and final bubble, respectively; $(b)$ vertical (red line, left axis) and horizontal (blue line, right axis) components of the separation. Open squares and circles: axisymmetric prediction. In $(a)$, the shape of the resulting bubble is shown at several successive time instants during the transient $24.5\leq t\leq 30$ following coalescence.

Returning to the near-coalescence stage, figure 14 shows how the interface topology evolves in the stages that just precede and follow numerical coalescence. In its late stage, the film is almost flat. No dimple has formed at its periphery, unlike what is customarily observed with coalescing drops (Hartland Reference Hartland1968, Reference Hartland1969; Jones & Wilson Reference Jones and Wilson1978). This is because, for low-to-moderate $Bo$ and fully mobile interfaces, a dimple forms only when the film has thinned down by several orders of magnitude (Yiantsios & Davis Reference Yiantsios and Davis1990). Indeed, with the same grid resolution, we observed a clear dimple for Bond numbers $\gtrsim 1$. Numerical coalescence takes place at $t\approx 24.44$. Considering that the film starts to form when the gap is $0.5R$-thick on the symmetry axis, one can estimate that the computation tracks the drainage process during approximately $1.15$ time units. The arguments discussed in Appendix B may then be used to estimate the time by which the film would actually rupture under real conditions. Figure 13$(a)$ indicates that the dimensionless approach velocity $\bar {V}_a=V_{TB}-V_{LB}$ of the two bubbles is approximately $0.5$ during the early stage of the drainage. From this, the approach Weber number, $We_a=\bar {V}_a^2Bo$, and the approach capillary number, $Ca_a=\bar {V}_aBo/Ga=\mu V_a/\gamma$, are found to be $0.125$ and $0.0125$, respectively. Then, for nearly spherical bubbles, the estimate (B1) predicts a dimensionless inertial drainage time $\bar {T}_{di}\approx 0.27$, while (B2) (once properly transposed to the case of two identical bubbles) predicts a viscous drainage time $\bar {T}_{dv}\approx 2.5$. Hence, the drainage is controlled by viscous effects and the limited resolution would shorten it by $\approx 1.35$ time units, would the initial bubble oblateness be small. Given that the tandem rises with the average velocity $\bar {V}_m=\frac {1}{2}(V_{TB}+V_{LB})\approx 2.4$, this implies that the vertical position at which coalescence actually happens would be underestimated by $3.25$ bubble radii in this limit. The non-negligible bubble oblateness certainly increases the actual drainage time. For instance, the actual value of the parameter $k_i$ in (B1) is $2.05$ for $\chi =1.5$ (Duineveld Reference Duineveld1994), so that the correct estimate for the inertial drainage time is $\bar {T}_{di}\approx 0.5$. The quantitative influence of the bubble oblateness on the viscous drainage time is unknown but presumably increases also significantly $\bar {T}_{dv}$.

Figure 14. Numerical coalescence process. $(a)$ Cross-sectional bubble shapes; $(b)$ zoom on the dashed rectangle in the left image of $(a)$, showing the successive grid refinements and the grid distribution in the near-contact region just before coalescence ($t=24.4$). In the first three images of $(b)$, the thin black line indicates the interface and the colour scale spans six grid levels, from $\varDelta =R/8.5$ (dark blue) to $\varDelta =R/272$ (brown).

As the last two snapshots in figure 14$(a)$ show, the neck radius of the resulting bubble grows very fast after coalescence happens. More specifically, the growth law (not shown) is found to be $(t-t_{coal})^{0.36}$ over the very first stage following the coalescence time instant $t_{coal}$. This is somewhat slower than the $(t-t_{coal})^{1/2}$ self-similar behaviour observed in the experiments of Paulsen et al. (Reference Paulsen, Carmigniani, Kannan, Burton and Nagel2014) and confirmed theoretically and numerically by Munro et al. (Reference Munro, Anthony, Basaran and Lister2015) and Anthony et al. (Reference Anthony, Kamat, Thete, Munro, Lister, Harris and Basaran2017), respectively. However, Paulsen et al. (Reference Paulsen, Carmigniani, Kannan, Burton and Nagel2014) noticed that the growth rate reduces gradually once the neck radius is larger than $0.3R$. Here, due to the limitations inherent to minimum cell size $\delta _{min}=R/272$ and to the numerical procedure used to let the initial two interfaces merge, we barely observe the neck before it reaches this radius, which is certainly the reason for the above slower growth rate.

7.1. Influence of bubble deformation in the DKT and ASE regimes

Figure 15 shows how several characteristics of the path of bubble pairs with $Ga=10$ vary with the Bond number in the range $0.02\leq Bo\leq 0.2$, i.e. $8\times 10^{-10}\leq Mo\leq 8\times 10^{-7}$. The DKT scenario is observed for the lower three values of $Bo$, while the pair with $Bo=0.2$ eventually experiences head-on coalescence. For lower $Bo$, the interaction always leads to the side-by-side configuration and the lateral separation stabilizes at a value close to $2.5$. While the Bond number does not have any noticeable influence on the final path in this regime, it affects the critical time by which the tumbling process sets in. Owing to the connection between the bubble oblateness and the amount of vorticity produced at its surface, the more oblate the LB is the stronger the attractive wake effect is. Therefore, at a given time, the larger $Bo$ the shorter the separation between the two bubbles. From figure 15$(b{-}d)$, it may be inferred that the axial symmetry of the flow breaks down when $\bar {S}=\bar {S}_c\approx 2.6$, a critical value reached in a shorter time as $Bo$ increases.

Figure 15. Variations with the Bond number of some path characteristics for $Ga=10$. $(a)$ Side view of the path (the solid and dashed lines correspond to the LB and TB, respectively); $(b)$ vertical separation; $(c)$ horizontal separation; $(d)$ inclination of the line of centres. The square symbol denotes the position/time at which coalescence takes place for $Bo=0.2$.

Figure 16 shows the same path characteristics for $Ga=30$ and Bond numbers ranging from $0.02$ to $0.5$, i.e. $9.9\times 10^{-12}\leq Mo\leq 1.5\times 10^{-7}$. Here, all cases with $Bo\leq 0.45$ correspond to the ASE scenario, whereas the two bubbles coalesce in an asymmetric way for $Bo=0.5$ (see figure 7). In all non-coalescing cases, it is found that the larger $Bo$ the smaller the long-term lateral deviation of the TB (figure 16$a$). Hence, the long-term lateral separation decreases with $Bo$ (figure 16$b$), since the LB hardly deviates from its initial vertical path except for $Bo=0.45$. The vertical separation $\bar {S}$ is also seen to decrease with the Bond number. This is merely a geometrical consequence of the decrease of $\bar {S}_r$. Indeed, the shorter the time spent by the TB in its lateral motion, the shorter the slowing down of its vertical motion, hence the smaller the increase of the vertical separation during the lateral escape stage.

Figure 16. Variations with $Bo$ of some path characteristics for $Ga=30$. $(a)$ Side view of the path (the solid and dashed lines correspond to the LB and TB, respectively); $(b)$ vertical separation; $(c)$ horizontal separation; $(d)$ inclination of the line of centres. In $(a)$, the square symbol indicates the position/time at which coalescence takes place for $Bo=0.5$.

The long-term inclination of the bubble pair exhibits more complex variations. First, it increases with the Bond number up to $Bo=0.2$. Then it slightly decreases for $Bo=0.3$ until it experiences a sharp increase for $Bo=0.45$. In the whole range $Bo\leq 0.3$, the long-term inclination angle stands in the range $30^\circ <\theta <40^\circ$, in agreement with the experimental findings of Kusuno & Sanada (Reference Kusuno and Sanada2015) obtained under similar conditions ($Re<150$). The situation corresponding to $Bo=0.45$ is specific. Indeed, the two bubbles are close to coalescing asymmetrically at some point ($\bar {S}\approx 2$ at $t=22.5$). This is why the lateral motion of the TB significantly disturbs the subsequent motion of the LB which is seen to experience a series of damped oscillations before rising along a rectilinear, slightly inclined path. Moreover, since the vertical separation of the two bubbles is small just before the TB starts to drift laterally ($\bar {S}\approx 3.8$ at $t=18$), the inclination of tandem after this drift has been completed is significantly larger ($\theta \approx 60^\circ$) than those found with smaller $Bo$. We shall come back to the long-term evolution of this bubble pair later.

Beyond these geometric features, the main information provided by figure 16 is that the smaller $Bo$ is the earlier the side escape of the TB starts. This is in stark contrast with the observations made on figure 15 for the DKT regime, where bubble deformation is found to promote the tumbling process. Moreover, considering the critical time at which $\bar {S}_r$ and $\theta$ depart from zero for the various $Bo$ reveals that the corresponding vertical separation does not keep a constant value (figure 16$b$–$d$). Rather it decreases from $\bar {S}_c\approx 8$ for $Bo=0.02$ to $\bar {S}_c\approx 5.8$ for $Bo=0.3$. These values are significantly larger than the equilibrium separation predicted in the axisymmetric configuration which, according to figure 5, ranges from $\bar {S}_e\approx 6.1$ for $Bo=0.02$ to $\bar {S}_e\approx 4.3$ for $Bo=0.3$. Moreover, in all cases, this critical separation is much larger than that at which tumbling is found to start in the DKT regime. This finding indicates that, unlike the breakdown of the axial symmetry in the DKT scenario, the ASE mechanism is driven by a long-range interaction. This is in line with the conclusion of Yin & Koch (Reference Yin and Koch2008) who simulated the motion of buoyancy-driven suspensions of spherical non-coalescing bubbles at $Re\approx 10$.

To get additional insight into the long-term behaviour of the bubble pair, figure 17 shows the evolution of the separation and inclination angle over longer times for the two sets $(Ga, Bo)=(30, 0.3)$ and $(Ga, Bo)=(30, 0.45)$ , i.e. $Mo=3.3\times 10^{-8}$ and $Mo=1.1\times 10^{-7}$, respectively. In the former case, the two components of the separation are seen to slightly increase until the end of the computation but the inclination angle changes by less than $1^\circ$ over the last forty time units, reaching the ‘asymptotic’ value $\theta \approx 36^\circ$. Conversely, in the latter case, the two components of the separation exhibit significant variations throughout the computation. The vertical separation gradually tends to zero while $\bar {S}_r$ goes on increasing even for $t\approx 100$. At this final time, $\theta$ is close to $83^\circ$ and is still gently increasing, so that there is little doubt that the system eventually reaches a perfect side-by-side configuration. To understand why the two bubble pairs behave so differently on the long term, it must first be noticed that, at the end of the lateral drift of the TB ($t\approx 40$), the distance $\mathscr {S}=(\bar {S}^2+\bar {S}_r^2)^{1/2}$ between the two centroids is approximately $6.8$ for $Bo=0.3$ and $4.5$ for $Bo=0.45$, the Reynolds numbers being approximately $66$ and $60$, respectively. For spherical bubbles, the results of Hallez & Legendre (Reference Hallez and Legendre2011) indicate that, for $Re=100$, the torque acting on the tandem is negligibly small for $\theta \gtrsim 25^\circ$ when $\mathscr {S}\gtrsim 6$ but has significant values whatever $\theta$ for $\mathscr {S}=4.5$. Although bubble deformation is large in the two sets considered in figure 17, this is a strong indication that, after the TB has drifted laterally, only the pair with $Bo=0.45$ still experiences a noticeable torque. Consequently only this pair may reach the side-by-side configuration in a reasonable time. According to figure 16, $\mathscr {S}$ decreases with $Bo$. Hence, it may be concluded that in the ASE scenario, the distance between the two centroids after the TB has completed its drift is usually too large for the torque to remain significant. Therefore the inclination of the tandem remains almost unchanged in the subsequent stages. Only bubbles whose deformation is close to the coalescence threshold (here $Bo\approx 0.5$) escape this rule, since the vertical separation of the corresponding pairs is very small when the side escape of the TB takes place.

Figure 17. Evolution of some characteristics of a bubble pair with $(Ga, Bo)=(30,0.45)$. $(a)$ Vertical (red line) and horizontal (blue line) separation distances; $(b)$ inclination of the line of centres.

Interestingly, for $Bo=0.45$, $\bar {S}_r$ goes on increasing at the end of the computation, suggesting that the interaction force is still non-zero and repels the two bubbles. At the same $Re$ and $\mathscr {S}$, the computational results of Legendre et al. (Reference Legendre, Magnaudet and Mougin2003) indicate that the interaction force between two spherical bubbles rising side by side is attractive. The difference between the two situations may be understood by noting once again that wake effects are much stronger for $Bo=0.45$ than for $Bo=0$. Since these effects are repulsive in the side-by-side configuration, the critical Reynolds number at which the interaction force switches from repulsive to attractive is significantly larger in the present case, which explains the observed behaviour.

7.2. Role of the TB shape in near-critical conditions for coalescence

The global geometric indicators reported in figure 16 reveal how the two bubbles get closer to coalescence as the Bond number increases. Nevertheless it is of interest to examine also how the bubble shapes vary with $Bo$ close to the corresponding threshold, especially with respect to the potential influence of the bubble distortion on the magnitude and even the direction of the sideways force acting on the TB through the two mechanisms reviewed in Appendix C. Figure 18 displays the evolution of the shapes and relative positions of the two bubbles for near-critical conditions at $Ga=30$ and $15$, respectively. It is striking that, for both values, the TB shape exhibits almost no left/right asymmetry until the late stage of the interaction. This is a strong indication that the $A$-mechanism discussed in Appendix C cannot reduce significantly the sideways force. Similarly, the flow characteristics are such that the $S$-mechanism cannot take place. For instance, with $Ga=30$ and $Bo=0.4$, the aspect ratio of the TB at $t=14$ is approximately $1.6$ and its Reynolds number is close to $80$. For these parameters, the path of an isolated bubble rising in a fluid at rest is stable and the numerical results of Adoua, Legendre & Magnaudet (Reference Adoua, Legendre and Magnaudet2009) (their figure 5) indicate that finite-Reynolds-number effects decrease the lift force by only $15\,\%$ with respect to the inviscid prediction. So, it can be concluded that, within the parameter range considered in the present study, the lateral migration of the TB is merely driven by the standard shear-induced mechanism. Hence, close to the threshold, what makes the difference between non-coalescing and coalescing situations is essentially the enhancement of the attractive wake effect as the TB becomes increasingly oblate. For a sufficiently large $Bo$, this attractive effect becomes so strong that it delays the occurrence of the lateral instability of the TB to such an extent that, although this bubble subsequently migrates ‘normally’, its migration lasts for a too short time to prevent it from hitting the LB.

Figure 18. Evolution of the bubble shapes and relative positions for near-critical conditions: $(a)$ $Ga=30$; $(b)$ $Ga=15$.

7.3. Influence of an initial angular deviation

Despite sophisticated bubble release systems, tiny initial lateral deviations can hardly be avoided in laboratory experiments, yielding small-but-non-zero angular deviations of the tandem with respect to the perfect in-line configuration. This calls for the investigation of the influence of an initial non-zero inclination on the subsequent dynamics. For this purpose, we systematically examined the impact of an initial inclination $0^\circ <\theta _0\leq 2^\circ$ throughout the $(Ga, Re)$-range considered in §§ 5 and 6. Note that in all configurations considered in this subsection, the initial horizontal separation between the bubble centres is several times larger than the finest grid cells located on both sides of the bubbles surface ($4$ times larger for the smallest inclination). This initial configuration results in a ‘macroscopic’ asymmetry of the discretized solution, in which the tiny numerical asymmetries of the pressure field returned by the Poisson solver (see § 6.2 and Appendix A) play no role.

For $Ga=10$, we observed that under such conditions, the DKT scenario no longer takes place for $Bo\leq 0.1$. Instead, the system follows a clear ASE evolution. For instance, selecting $(Ga, Bo)=(10, 0.1)$, the TB starts drifting laterally when $\bar {S}\approx 4$ (respectively $\bar {S}\approx 5$) if $\theta _0$ is set to $1^\circ$ (respectively $2^\circ$), leaving the path of the LB almost unaffected in both cases. Instead of ending up in the side-by-side configuration as it does when $\theta _0=0^\circ$, the tandem then reaches an approximate final inclination of $50^\circ$ (respectively $40^\circ$). That $\theta _0$ has such a dramatic influence on the existence of the DKT regime and the late geometry of the tandem is at variance with observations reported for rigid bodies falling in tandem. Indeed, experiments performed with short cylinders (Brosse & Ern Reference Brosse and Ern2014) indicate that the DKT regime is still observed when a small initial lateral offset is imposed to the trailing body. The key difference with a bubble pair is the relative magnitude of attractive effects which is much stronger for rigid bodies, owing to the different vorticity generation mode at the body surface. Moreover, for a given body shape, shear-induced lift effects are significantly weaker for rigid bodies with Reynolds numbers of $O(10)$ or larger, owing to the presence of large separated regions (Kurose & Komori Reference Kurose and Komori1999). These two factors make the DKT configuration much more sensitive to small lateral deviations in the case of nearly spherical bubbles.

Figure 19 provides the evolution of the tandem geometry in the case $(Ga, Bo)=(30, 0.3)$ for initial inclinations $\theta _0$ ranging from $0^\circ$ to $2^\circ$. The system is found to follow the ASE scenario in all cases. However, its final geometry strongly depends on $\theta _0$. In particular, as figure 19($d$) shows, the final inclination of the line of centres decreases from $36^\circ$ for $\theta _0=0^\circ$ to $22^\circ$ for $\theta _0=2^\circ$. This is because when $\theta _0$ is non-zero, the TB starts drifting laterally soon after it is released from rest, which is not the case for $\theta _0=0^\circ$. Indeed, with $\theta _0\neq 0^\circ$, the initial flow configuration is no longer axisymmetric and the TB faces an asymmetric wake as soon as vorticity generated at the LB surface has been advected downstream over an $O(\bar {S})$-distance. Consequently, the TB experiences a non-zero sideways force much earlier than in the perfect in-line configuration, where this force occurs only after the system has become unstable. A non-zero $\theta _0$ shortens the initial time period after which $\bar {S}_r$ starts to grow, as is made clear in figure 19$(c)$. A direct consequence of this shortening is the fact that the vertical separation at which the lateral drift starts is larger for $\theta _0\neq 0^\circ$ ($\bar {S}\approx 8$ for $\theta _0=2^\circ$ instead of $\bar {S}\approx 6.8$ for $\theta _0=0^\circ$). For this reason, weaker velocity gradients across the LB wake exist for $\theta _0\neq 0^\circ$ at the position of the TB when it starts drifting, resulting in a weaker sideways force. This makes the maximum of $V_r$ smaller (figure 19$b$), yielding a shorter final lateral position, hence a smaller final inclination of the tandem.

Figure 19. Evolution of several characteristics of a bubble pair with $(Ga, Bo)=(30, 0.3)$ undergoing initial angular deviations $\theta _0$ from $0^\circ$ to $2^\circ$. $(a)$ Vertical velocity of the LB (solid line) and TB (dashed line); $(b)$ horizontal velocity of the TB; $(c)$ vertical (solid line, left axis) and horizontal (dashed line, right axis) separations; $(d)$ inclination of the line of centres. In $(b)$, the bullets on the $\bar {S}_r$-curve indicate the location where the maxima of $V_r$ are reached.

Figure 20 displays the influence of the Bond number on the evolution of bubble pairs with $Ga=30$ when $\theta _0$ is set to $2^\circ$. This figure is the counterpart of figure 16 discussed above for $\theta _0=0^\circ$. All pairs with $Bo<0.7$, i.e. $Mo<4.2\times 10^{-7}$, are seen to follow a clear ASE scenario. For the aforementioned reason, the final inclinations reported in figure 20 are significantly smaller than those observed in figure 16. In addition, they exhibit a marked and consistent decrease with the Bond number, a trend which is absent from figure 16. For pairs with $Bo\geq 0.8$, $\theta$ is found to decrease over time until the two bubbles perfectly align vertically, which eventually forces them to coalesce. The near-critical case $Bo=0.7$ is quite specific, and in many instances similar to the situation encountered for $Bo=0.45$ in figure 16. In this case, the two bubbles are very close to coalescing at $Y\approx 58$ in figure 21$(a)$. The corresponding gap is so thin that the flow at the back of the LB is strongly disturbed, forcing the latter to deviate abruptly from its vertical path. This allows the lateral separation to increase beyond the critical value $\bar {S}_{rc}\gtrsim 2$ within a short time lapse (not visible in figure 21$(b)$ which is limited to shorter times), allowing the tandem to avoid coalescence. Similar to the evolution displayed in figure 17$(b)$, but through a sharper transition, the inclination of the tandem grows until the side-by-side configuration is reached. Given the small gap at which the transition takes place and the quite symmetric final lateral positions of the two bubbles, this specific evolution is closer to the DKT scenario than to a standard ASE evolution.

Figure 20. Variations with $Bo$ of some path characteristics for $Ga=30$ in the presence of an initial angular deviation $\theta _0=2^\circ$. $(a)$ Front view of the path (the solid and dashed lines correspond to the LB and TB, respectively); $(b)$ horizontal separation; $(c)$ inclination of the line of centres.

Figure 21. Phase diagram showing whether a bubble pair released with an angular deviation $\theta _0=2^\circ$ evolves in the ASE regime (squares) or eventually coalesces (triangles). $(a)$ $(Ga, Bo)$ map; $(b)$ $(Re, \chi )$ representation based on the steady-state values corresponding to the rise of the corresponding isolated bubble.

The critical Bond number beyond which the two bubbles coalesce stands in the range $0.45-0.5$ in figure 16 ($\theta _0=0^\circ$), but is slightly larger than $0.7$ for $\theta _0=2^\circ$. In line with the discussion on figure 19, the reason for this marked increase stems directly from the longer vertical distance over which the shear-induced lift force acts on the TB when $\theta _0\neq 0^\circ$, and the slightly shorter lateral distance this bubble has to drift to avoid hitting the LB. These two factors imply that, compared with the reference case, a weaker positive lift force is sufficient to avoid coalescence when $\theta _0$ is non-zero. The two lift reversal mechanisms discussed in Appendix C being directly linked to the ability of the TB to deform, increasing the Bond number makes the lift force decrease and eventually change sign, other things being equal. Hence, in the presence of an initial angular deviation, coalescence can be avoided over a broader range of distortion of the TB, i.e. up to a larger Bond number. The fate of all bubble pairs released with an initial inclination $\theta _0=2^\circ$ over the range $10\leq Ga\leq 30, 0< Bo\leq 1.0$ is summarized in figure 21. This figure is the counterpart of figure 6 obtained with $\theta _0=0^\circ$. In line with the above discussion, the ASE regime observed when $\theta _0=2^\circ$ exists over a significantly broader range of $Bo$, hence for bubbles with a larger oblateness. For instance, bubbles with aspect ratios $\chi \gtrsim 1.3$ (respectively $1.7$) are found to coalesce for $Ga=15$ (respectively $30$) in the absence of any initial deviation. With $\theta _0=2^\circ$, the corresponding critical aspect ratios raise beyond $1.5$ (respectively $2.0$) and the critical Reynolds numbers are close to $50$ and $120$, respectively. Figure 22 shows how the bubble shapes and relative positions evolve close to the transition to coalescence for $Ga=30$. Unlike those reported in figure 18, the TB shapes now exhibit a marked left/right asymmetry (although less pronounced, this trend still exists for lower $Ga$). The difference with the strict in-line configuration is that here the TB is fully immersed in an asymmetric flow as soon as the wake of the LB has developed within the gap. Consequently, its shape has to adapt to this asymmetric environment throughout the interaction process. While the TB is seen to migrate toward the right for $Bo=0.7$, a tiny migration toward the left may be identified for $Bo=0.8$, leading unavoidably to coalescence. The discussion in Appendix C allows the origin of this reverse migration to be readily identified. The aspect ratio of the TB being only $1.85$ for $Bo=0.8$, the $S$-mechanism is not present, since it only takes place beyond a threshold $\chi _{cS}\approx 2.2$ ( we observed a reversed migration due to this mechanism by increasing the Bond number beyond unity). In contrast, the egg-like shapes of the TB point to the $A$-mechanism. The capillary number based on the rise velocity of this bubble at $t=14$ only increases from $0.051$ for $Bo=0.7$ to $0.058$ for $Bo=0.8$. However, this modest increase turns out to be sufficient for the deformation-induced force directed to the left to exceed the shear-induced lift directed to the right, preventing the lateral escape of the TB.

Figure 22. Evolution of the bubble shapes and relative positions at $Ga=30$ under near-critical $Bo$-conditions in the presence of an initial angular deviation $\theta _0=2^\circ$.

7.4. Influence of initial separation

All simulations considered up to now are based on an initial separation distance $\bar {S}_0=8$. It is relevant to examine how this choice influences the upcoming evolution of the bubble pair. Since the in-line configuration is stable in the head-on coalescence regime, $\bar {S}_0$ is not expected to have an influence in this regime. This is why we consider the influence of $\bar {S}_0$ only in the DKT and ASE regimes.

Figure 23 displays the evolution of some characteristics of the bubble pair obtained by increasing $\bar {S}_0$ from $8$ to $12$ in the case $Ga=10, Bo=0.02$. With $\bar {S}_0=8$, this set of parameters yields the DKT configuration described in figures 8 and 9. Introducing an appropriate $\bar {S}_0$-dependent time shift $t_0$ and setting $t^*=t-t_0(\bar {S}_0)$ allows the $t^*$-evolutions of all quantities for the three initial separations to collapse on a single curve beyond $t^*\approx 10$. Hence, a DKT scenario yielding the same final configuration is observed whatever the initial separation. This is no real surprise since, for this set of parameters, the in-line configuration becomes unstable only when the two bubbles get very close, which happens only for $t^*\approx 40$.

Figure 23. Evolution of some characteristics of a bubble pair with $(Ga, Bo)=(10, 0.02)$ for different initial separations. $(a)$ Rise speed of the LB (solid lines) and TB (dashed lines); $(b)$ vertical (solid lines, left axis) and horizontal (dashed lines, right axis) separations. Red, blue and green lines refer to $\bar {S}_0=8, 10$ and $12$, respectively. A time shift $t_0=11$ (respectively $t_0=22$) is applied for $\bar {S}_0=10$ (respectively $S_0=12$) and evolutions are plotted vs the modified time $t^*=t-t_0$.

Figure 24 shows how the evolution of the same characteristics varies with $\bar {S}_0$ for $Ga=30, Bo=0.3$. In this case, the interaction process yields an ASE scenario whatever $\bar {S}_0$. However, the final configuration now depends on the initial separation. In particular, the shorter $\bar {S}_0$ the larger the final inclination of the tandem, with $\theta \approx 47^\circ$ and $\theta \approx 22^\circ$ for $\bar {S}_0=6$ and $\bar {S}_0=12$, respectively. Not unlikely, the larger $\bar {S}_0$ is the longer it takes for the path of the TB to start bending. Examining the evolutions of the vertical separation and those of the rise speed of each bubble reveals that none of these three quantities exhibits a constant value by the time the TB starts drifting laterally. In contrast, it turns out that the difference $V_{TB}-V_{LB}$ is close to $0.35$ whatever $\bar {S}_0$ when this lateral motion sets in. Therefore we conclude that the in-line configuration becomes unstable when the velocity difference between the two bubbles exceeds the above threshold, irrespective of the separation distance at which this threshold is reached.

Figure 24. Influence of the initial separation $\bar {S}_0$ on several characteristics of a bubble pair with $(Ga, Bo)=(30, 0.3)$. $(a)$ Vertical velocity component of the LB (solid line) and TB (dashed line); $(b)$ same for the horizontal component; $(c)$ vertical (solid line) and horizontal (dashed line) separations; $(d)$ inclination of the line of centres. Red, blue, green and orange lines correspond to $\bar {S}_0=6, 8, 10$ and $12$, respectively.

Closely related to the influence of the initial separation is that of the sequential release of the two bubbles in actual laboratory experiments. Still with $Ga=30, Bo=0.3$, we ran a computation in which the TB was allowed to start rising only after the LB reached a dimensionless height $\bar {S}(t_0)=8$ above the point of release. Figure 25$(a)$ displays the corresponding evolution of the vertical separation. Due to the time lapse the TB needs to reach its final rise speed, $\bar {S}$ reaches a maximum close to $10.2$ before the attractive interaction sets in and the separation starts to decrease. This is why we compare the subsequent evolution with that obtained for $\bar {S}_0=10$ when the two bubbles released simultaneously. For this purpose, we introduce an appropriate time shift $t_0$ and define a modified time $t^*=t+t_0$ to make the two evolutions coincide in the early stage of the interaction. With this time shift, the two further evolutions of $\bar {S}$ almost superimpose, as do those of the lateral velocity of the TB (figure 25$b$) and the inclination angle (figure 25$c$). Consequently, it can be concluded that the time elapsed in between the release of two successive bubbles merely increases their ‘initial’ separation, defined as the distance that separates them before the wake of the LB starts to influence the rise of the TB.

Figure 25. Influence of the sequential release of the two bubbles. $(a)$ Vertical separation; $(b)$ lateral velocity; $(c)$ inclination angle. Solid line: bubbles are released simultaneously with $\bar {S}_0=8$; blue dash-dotted line: same with $\bar {S}_0=10$ and a time shift $t_0=2.6$; red dashed line: the TB is released from rest after the LB has travelled a distance $\bar {S}(t_0)=8$.