Fate of bubble clusters rising in a quiescent liquid

Abstract

We use experiments to study the evolution of bubble clusters in a swarm of freely rising, deformable bubbles. A new machine learning-aided algorithm allows us to identify and track bubbles in clusters and measure the cluster lifetimes. The results indicate that contamination in the carrier liquid can enhance the formation of bubble clusters and prolong the cluster lifetimes. The mean bubble rise velocities conditioned on the bubble cluster size are also explored, and we find a positive correlation between the cluster size and the rise speed of the bubbles in the cluster, with clustered bubbles rising up to $20\,\%$ faster than unclustered bubbles.

1. Introduction

Bubbles rising though a liquid is a frequently occurring situation in nature (e.g. bubble plumes rising from the bottom of a lake), daily life (e.g. bubble chains rising in carbonated drinks) and technology (e.g. waste water treatment). For the case of bubble swarms rising freely in a quiescent fluid, the bubbles tend to distribute inhomogeneously, spontaneously forming clusters. When the bubbles are large enough, the rising bubbles generate turbulence in the liquid, a phenomenon referred to as bubble-induced turbulence (BIT). This BIT in turn influences the clusters and understanding this is important both for its own sake, and also due to its impact on many aspects of bubble motion, including their rise velocities, collisions, dispersion and the intensity of the BIT generated (Takagi, Ogasawara & Matsumoto 2008; Tagawa et al. 2013; Lohse 2018; Liao et al. 2019; Ma, Lucas & Bragg 2020a; Ma et al. 2021).

Early numerical investigations (e.g. Smereka 1993) on bubble clustering in the presence of BIT assumed spherical bubbles and considered potential flow, and found that the bubbles form large horizontal clusters as they rise. However, subsequent experiments and more realistic numerical simulations did not clearly observe such clusters, but 2-bubble clusters with a wide range of orientations were identified as the most commonly occurring clusters (Zenit, Koch & Sangani 2001; Bunner & Tryggvason 2003; Esmaeeli & Tryggvason 2005). The prevalence of 2-bubble clusters has motivated the community to explore the dynamics of bubble pairs with varying separations and orientations, usually for the case of clean bubbles (Hallez & Legendre 2011). The two extreme cases are bubbles aligned side by side (Legendre, Magnaudet & Mougin 2003) and in line (Zhang, Ni & Magnaudet 2021). The results indicate that the side-by-side configuration is more stable than the in-line case, and this is because for the in-line configuration a slight transverse movement of the trailing bubble relative to the leading bubble makes it ‘feel’ a shear flow, that can drive the trailing bubble out of the leading bubble's wake. Nevertheless, stable in-line bubble chains are often observed in carbonated drinks. This contradiction was recently explained by Atasi et al. (2023) as being due to the combined effects of bubble deformation and contamination in such liquids that can result in a reversal of the lift force and stable chain. In addition to exploring the stability of nearby bubble pairs, the impact of neighbouring bubbles on their rise velocity has been investigated and compared with the case of isolated bubbles. Hallez & Legendre (2011) showed that the side-by-side configuration maximizes the drag force acting on a pair of bubbles, while the in-line bubble configuration minimizes the drag due to wake entrainment for the trailing bubble.

These studies on bubble pair dynamics have provided much insight, however, there are many open questions concerning the behaviour of bubble swarms where two or more bubbles may be clustered together, whose motion may also be affected by the wakes of other bubbles and bubble clusters in the flow. Indeed, while the rise velocity of bubble pairs in a quiescent liquid is well understood, their behaviour in the context of bubble swarms is debated. For example, the experiments of Stewart (1995) and Brücker (1999) for large deformable bubbles in a swarm found that the mean rise velocity was considerably larger than that for a single bubble. However, this contradicts other experimental (Ishii & Zuber 1979) and numerical (Roghair et al. 2011) studies for large bubbles, which argue that the mean bubble rise velocity decreases monotonically as the gas void fraction is increased.

Several fundamental questions remain mostly unexplored: What is the probability of forming clusters involving $N_b$ number of bubbles? What is the lifetime of these clusters? How does the rise velocity of bubbles in a cluster depend on $N_b$? How do the answers to these questions depend on contaminants in the liquid? In this paper we explore these questions experimentally by tracking thousands of deformable bubbles in a vertical column, using a recently developed machine-learning algorithm to detect and follow the evolution of bubble clusters, and explore how the bubble rise velocities depend on $N_b$. We also consider the effect of surfactants to provide a more complete picture for real systems where contaminants may cause behaviour that differs substantially from that of an idealized clean system.

2. Experimental method

2.1. Experimental set-up

The experimental apparatus is identical to that in Ma et al. (2022), and we therefore refer the reader to that paper for additional details; here, we summarize. The experiments were conducted in a rectangular bubble column (depth 50 mm and width 112.5 mm), with a water fill height of 1000 mm (figure 1). Air bubbles are injected through 11 spargers which are homogeneously distributed at the bottom of the column.

Figure 1. Sketch of the bubble column used in the experiments (note that, in the actual experiment, the number of bubbles in the column is $O(10^3)$). The sketch is not to scale. The insets show the cylindrical cluster search region: view from (a) the top, (b) the side and (c) the front. Note that (a,b) are only schematic representations, showing possible bubble arrangements, while (c) shows a cut from a recorded snapshot. The green outline refers to bubbles that are considered inside the sharp centre region.

We use tap water in the present work as the base liquid and consider two different bubble sizes by using spargers with different inner diameters. Note that the tap water will already be slightly contaminated prior to adding the surfactants or salts, and the bubbles can behave differently in this tap water compared with that in a pure water system (Craig, Ninham & Pashley 1993; Takagi & Matsumoto 2011). For the bubble sizes considered in the present study (with $d_b>2$ mm), however, the slight contamination in the tap water does not have a significant effect on the bubble shape/motion, as confirmed by our previous studies (Hessenkemper et al. 2021; Ma et al. 2023). For each bubble size, we manipulate the gas flow rate and ensure that all cases are not in the heterogeneous regime of dispersed bubbly flows. In total, we have six mono-dispersed cases (see supplementary movies 1–6 available at https://doi.org/10.1017/jfm.2023.807) labelled as SmTapLess, SmTap, SmPen+, LaTapLess, LaTap and LaPen+ in table 1, including some basic characteristic dimensionless numbers for the bubbles. Here, ‘Sm/La’ stand for smaller/larger bubbles, ‘Pen+’ stands for corresponding cases added by 1000 ppm 1-Pentanol and ‘Less’ stands for lower gas void fraction than **Tap/**Pen+ cases for smaller/larger bubbles, respectively. It should be noted that the three cases with larger bubble sizes have higher gas void fractions than the three cases with smaller bubbles. This is because in our set-up it is not possible to have the same flow rate for two different spargers while also maintaining a homogeneous gas distribution for mono-dispersed bubbles. Furthermore, the bubble size is slightly reduced when adding 1-Pentanol for both types of sparger. This is due to the influence of the surfactants that reduce the surface tension and hence affect the bubble formation at the rigid orifice.

Table 1. Selected characteristics of the six bubble swarm cases. Here, $\alpha$ is the averaged gas void fraction, $d_b$ the equivalent bubble diameter, $\chi$ the aspect ratio, $L_b$ the inter-bubble distance, $Ga\equiv \sqrt {|{\rm \pi} _\rho -1|gd_{b}^{3}}/\nu$ the Galileo number and $Eo\equiv \Delta \rho gd_b^2/\sigma$ the Eötvös number. The bubble Reynolds number $Re_b$ and drag coefficient $C_D$ are based on $d_b$ and the bubble to fluid relative velocity.

To identify and track bubble clusters, we use planar shadow images obtained by recording the flow with a high-speed camera and illuminating the set-up with a LED. The measurement resolution in time and space are 250 f.p.s. and 59.9 $\mathrm {\mu }$m Px$^{-1}$, respectively, with a field of view (FOV) of $90\,\mathrm{mm}\times 76\ \mathrm {mm}$. For each case, we record 1000 sequences with each having 70–75 frames – approximately the time in which a bubble passes through the complete image height.

2.2. Bubble identification and tracking

In our study one camera is used and as such the bubbles can be identified and tracked only in two dimensions. However, as will be shown in § 2.3, we are nevertheless able to identify bubble clusters in quasi-three dimensions. Independent of the number of cameras used, the task to track bubbles in image sequences can be done in a detect-to-track or in a track-to-detect fashion. While the former links previously detected bubbles in each frame to form suitable tracks, the latter uses extrapolations of already established tracks to detect bubbles in follow-up images. We use the former detect-to-track strategy, which allows us to incorporate detections among multiple frames to establish tracks with, however, relying more strongly on an accurate detector that finds bubbles in individual frames.

Even for low gas volume fractions, detecting bubbles in individual images is a challenging task since bubbles can overlap in the images. Fully overlapping bubbles cannot be detected, but partially overlapping bubbles can be dealt with and deep-learning-based strategies for this have recently shown very promising results (e.g. Kim & Park 2021). In our previous work (Hessenkemper et al. 2022), we developed such an approach that used a trained convolutional neural network (CNN) to segment overlapping bubbles. Furthermore, the contour of each detected bubble is reconstructed using 64 radial vectors pointing from the segmentation centre to the boundary (figure 2a), and the radial vectors of partly occluded bubbles are corrected using an additional multi-layer perceptron (MLP).

Figure 2. Initialization of the tracking algorithm. (a) Cut out snapshots of three different time steps, with three detected bubbles and their reconstructed radial vectors highlighted in colour. The dashed lines indicate possible connections of the right most bubble at $t_1$ to the successive time steps $t_2$ and $t_3$. (b) Illustrates a constructed graph with all possible connections of these three bubbles.

The subsequent tracking of multiple detected bubbles in close proximity poses further challenges, as the tracker not only has to be robust against inaccuracies of the detector, i.e. missing or false detections, but also has to be able to track bubbles that are fully occluded even for multiple times steps, while at the same time having numerous possible associations in the near vicinity. To solve these issues, a graph-based tracking formalism is used. Specifically, we follow the framework of Brasó & Leal-Taixé (2020), utilizing multiple MLPs to predict valid connections of detections on graph structured data. The four main aspects of this tracking framework are described as follows. Details on the network architectures and the training dataset are provided in Appendix A.

Graph construction: to track the bubbles, each sequence is modelled as a graph, with detections (bubbles) being the nodes of the graph and possible connections in time being the edges of the graph, i.e. a pair of detections forward or backward in time that are possibly from the same bubble (figure 2b). The task is then to classify the edges into active and non-active edges, which at the end form a set of valid tracks that fulfil the so-called ‘flow conservation constraints’ – each node having an active edge to at most one node forward in time and at most one node backward in time.

Feature encoding: for both the nodes and the edges of the graph, features are encoded with two separate MLPs (figure 3a). The node embeddings represent the appearance features of the detections, which are usually encoded with a CNN (Brasó & Leal-Taixé 2020). However, monochromatic bubble images show few distinct features, with the size and the shape of the bubble image being the most relevant ones. Thus, we have chosen the 64 radial vectors from the bubble detector as input for the node feature encoder, providing not only a more accurate description of the features bubble size/shape, but also a better two-dimensional (2-D) bubble contour in the case of overlapping bubbles due to the additional correction of the radial locations. The edge embeddings represent tracking-related features. For each detected pair of frames, the time increment, relative coordinate and size are fed into the edge feature encoding MLP to generate edge embeddings.

Figure 3. Steps of the tracking algorithm: (a) feature encoding with node encoder MLP (NE) and edge encoder MLP (EE) together with edge update MLP (EU); (b) time-aware node update MLP (NU) with past node update MLP (PNU) and future node update MLP (FNU); (c) predicted active edges.

Message passing network: the core of the tracking algorithm is the message passing network (MPN) whose main purpose is to update node and edge embeddings with respect to their surrounding nodes and edges in the graph, and this is done iteratively using message passing steps. First, the edge embeddings are updated by combining their embeddings with the embeddings of the adjacent pair of nodes and feeding them into an edge-update MLP (figure 3a). Then, a time-aware node update step is applied by aggregating edge embeddings of adjacent edges, which already contain information of connected nodes due to the previous edge update step. The time awareness is introduced by first aggregating and updating separately incoming edges, i.e. connections backward in time, and outgoing edges, i.e. connections forward in time, with individual MLPs and then concatenating the outcome to finally update the node embeddings with a node update MLP (figure 3b). For each iteration, information of nodes one step further in time is passed through the network to the node/edge to be updated. Thus, the number of iterations defines the time increment of the information of other nodes that are supplied to the current node.

Edge classification and post-processing: after updating all node and edge embeddings with the MPN, the edges are classified with a classifier MLP (figure 3c). The predictions are post-processed and remaining violations of the flow conservation constraints are treated with an exact rounding scheme (Brasó & Leal-Taixé 2020). Lastly, missing links in the trajectories are interpolated using bilinear interpolation and each trajectory is smoothed with a uniform filter. In supplementary movie 8 we provide qualitative results of the tracking algorithm applied on a sequence of case LaTapLess, while in Appendix B quantitative results are provided with the use of common multi-object tracking metrics. The results demonstrate that the algorithm can track the majority of the detected bubbles, with a slight performance decrease with increasing bubble size and increasing number of bubbles.

2.3. Identification and tracking of bubble clusters

The detection of bubble clusters at each time step follows a distance criterion between neighbouring bubbles whose centres in the 2-D image domain are below a predefined threshold $2d_b$ from each individual case (figure 1c). This value is mainly based on the work of Legendre et al. (2003) that a considerable drag enhancement is observed for a bubble pair rising side by side within this distance. Tests for different thresholds ($2d_b\pm 0.5d_b$) were conducted and the trends of the results in § 3 were found to be insensitive to the choice of this parameter. Furthermore, since we attempt to detect the bubble cluster in a quasi-3-D manner, we keep the in-focus region in the depth direction to also be $2d_b$ (figure 1a,b). To estimate this depth distance to the centre plane we use the grey value gradient of the detected bubbles and consider only sharp bubbles in the shallow depth of field (DoF) region (see Appendix C for more detail). In summary, we utilize a cylindrical search volume, ${\rm \pi} (2d_b)^2\times 2d_b$, for the cluster identification, radial in the 2-D image domain and linear in the depth direction. For all the cases, the mean inter-bubble distance $L_b$ based on the global void fraction (table 1) is much larger than the search radius $2d_b$, indicating that the bubble clusters to be discussed are dynamically significant.

The cluster tracking strategy is inspired by the work of Liu et al. (2020) for characterizing the temporal evolution of inertial particle clusters in turbulence. Considering two clusters identified in two consecutive time steps ($\Delta t=1/250$ s), we take both to be successive realizations of the same cluster if the number of bubbles they share is above a given threshold. The shared bubbles across clusters in successive time steps are termed connections. We consider forward-in-time and backward-in-time connections, and apply thresholds on the fraction of connected bubbles over the total number of bubbles in each cluster. We illustrate in figure 4(a) an example: cluster A (identified in time step 1) shares all its bubbles with cluster C (identified in time step 2), while C shares 2/3 of its bubbles with A. Therefore, the fractions of forward and backward connections between A and C are 1 and 2/3, respectively. On the other hand, B shares 1/3 of its bubbles with C, and C shares 1/3 of its bubbles with B. Thus, the forward and backward connections between B and C are 1/3 and 1/3, respectively. Following Liu et al. (2020), two clusters in consecutive time steps are identified as the same cluster when the fractions of their backward and forward connections are both $\geq 1/2$. In the example of figure 4(a), A and C are recognized as belonging to the same cluster. The cluster lifetime is defined as the time elapsed between birth (the first instance a cluster is identified) and death (the last time it is recognized). Here, we explicitly include the lower threshold of $1/2$, since many 2-bubble clusters appear and require an additional criterion for tracking. In figure 4(b) we give an example where cluster A at time step 2 splits into B and C. To decide whether B or C should be regarded as the continuation of cluster A for the purposes of tracking, we consider whether cluster B or C persists longer into the future. In this example, while cluster C survives until time n, B does not. Therefore, we regard C, D and A as belonging to the same cluster, while cluster B is considered to be a newborn cluster at time step 2. This approach eliminates ambiguities since it ensures that a cluster at any instant can only be associated with at most one cluster either in the past or future.

Figure 4. Method of tracking clusters: (a) illustrates an example of two clusters merging into one cluster; (b) illustrates an example of one cluster separating into two clusters.

3. Results

3.1. Probability of being clustered

We first consider the percentage of bubbles in the flow that are clustered based on the data points listed in table 2. The results in figure 5(a) show that this percentage increases in the order of SmTapLess, SmTap to SmPen+ for the smaller bubbles. While the increase from SmTapLess to SmTap is quite understandable due to the increase of the gas void fraction, the result from SmPen+ (whose $\alpha$ is sightly less than that of SmTap) shows that the surfactant promotes the formation of clusters. We obtained similar results for the three cases with larger bubbles (figure 5b).

Table 2. Number of bubbles tracked and number of clusters in the FOV used for computing the statistics.

Figure 5. Percentage of bubbles in cluster for the cases with smaller (a) and larger (b) bubbles. Probability of number of bubbles within a cluster for the different cases using linear (c) and log (d) plot, respectively.

In figure 5(c,d) we consider the probability of finding a given number of bubbles within a cluster for all cases, plotting the results using both linear and logarithmic vertical axes in order to highlight regions of both high and low probability, respectively. The results show that the probability decreases with increasing $N_b$ and, consistent with previous studies, $N_b=2$ is the most common cluster size for all 6 cases (Zenit et al. 2001; Bunner & Tryggvason 2003). However, the results also show that $N_b=3,4$ clusters occur with non-negligible probability, and there are even rare events with $N_b=8$ clusters. The results also show that adding contaminants decreases the probability of forming $N_b=2$ clusters, and increases the probability of forming larger clusters, although the dependence is not too strong. For a fixed contaminant level, increasing $\alpha$ has the same effect.

3.2. Lifetime

We now turn to consider the mean lifetime of the clusters, $\langle t_{life}\rangle$, as a function of $N_b$ (only the results for $N_b\leq 5$ are shown as the statistics for $N_b>5$ are not converged). Figure 6(a,b) shows $\langle t_{life}\rangle$ normalized by a characteristic time scale of BIT (Ma et al. 2017, 2020b), $t_{BIT}\equiv d_b/U_b$, where $U_b$ is the mean vertical slip velocity between the bubble and liquid at the column centre. The values of $\langle t_{life}\rangle /t_{BIT}$ are order unity, suggesting that $t_{BIT}$ is indeed a dynamically relevant time scale for the cluster lifetimes. The results also reveal a systematic dependence on $N_b$ and the liquid contamination. First, $\langle t_{life}\rangle$ decreases monotonically with increasing $N_b$, such that larger bubble clusters are not only rarer (see § 3.1), but also more unstable. While this may not seem surprising, it is in fact the opposite to what has been observed for inertial particles where the cluster size and its lifetime are positively correlated (Liu et al. 2020). The difference could be simply due to the fact that the most common values of $N_b$ for our clusters are much smaller than those for the inertial particles in Liu et al. (2020), and as a result relatively small changes in the bubble configurations can result in the formation or destruction of a given cluster. The other significant difference is that our bubbles hydrodynamically interact, unlike the numerically simulated inertial particles in Liu et al. (2020), where a one-way coupling assumption is used. Second, increasing $\alpha$ not only leads to the formation of larger clusters, but also slightly longer mean lifetimes for the clusters, although the lifetimes for $N_b=2$ are the least sensitive to $\alpha$. Third, the mean lifetimes of the bubble clusters notably increase with increasing contamination levels. In a recent paper we showed that increased contamination leads to a reduction of $Re_b$ and an increase in BIT (Ma et al. 2023). The reduction in $Re_b$ causes the bubble trajectories to be less chaotic, and this may explain why the cluster lifetimes increase with increasing contamination.

Figure 6. Mean lifetime of 2-, 3-, 4- and 5-bubble clusters (a,b) and PDF of the bubble cluster lifetime (c,d) for smaller bubble cases (a,c) and larger bubble cases (b,d).

Figure 6(c,d) shows the probability density functions (PDFs) for the cluster lifetimes, which have been computed using clusters of all sizes. The general dependence on the flow variables is similar to that observed for the mean cluster lifetime, with the PDF tails becoming increasingly heavy in the order TapLess, Tap and Pen+ for both the small and large bubbles. The majority of the bubble clusters survive for $t_{life}/t_{BIT}=O(1)$, however, there are extreme cases where clusters survive for up to $t_{life}/t_{BIT}\approx 15$ for the Pen+ cases. Such extreme cases reflect the strongly nonlinear and non-equilibrium nature of the flow, with large fluctuations about the mean-field behaviour being an essential feature of turbulent flows. The central regions of the PDFs are approximately described by stretched exponential functions, however, the range of $t_{life}$ over which this applies is quite small.

3.3. Mean $N_b$-bubble cluster rise velocity

We finally consider the role that clustering plays in the mean bubble rise velocity. Figure 7 shows the mean rise velocity of bubbles in clusters $U_b$ conditioned on $N_b$, and the results for unclustered bubbles $N_b=1$ are also shown for reference. Consistent with our previous results based on averaging over all bubbles (Ma et al. 2023), the results show that, for almost all $N_b$, increasing the liquid contamination leads to a reduction in $U_b$, due to the modification of the bubble boundary conditions. For the larger bubbles we also observe a clear increase in $U_b$ with increasing $N_b$, with an increase of up to $20\,\%$ when going from unclustered bubbles ($N_b=1$) to bubble pairs ($N_b=2$), while the increase is more moderate when $N_b$ is increased beyond 2. The enhancement of $U_b$ when going from $N_b=1$ to $N_b=2$ is also observed in the SmPen+ case, with only slight enhancements when $N_b$ is increased beyond 2. However, for the SmTapLess and SmTap cases, $U_b$ varies only weakly with $N_b$, even in going from $N_b=1$ to $N_b=2$.

Figure 7. Mean bubble rise velocity as a function of the number of bubbles $N_b$ in the cluster: (a) smaller bubbles and (b) larger bubbles. Here, $N_b=1$ denotes unclustered bubbles.

What is the physical explanation for why the clustered bubbles rise faster than unclustered bubbles? We begin by considering the case of bubble pairs $N_b=2$ and plot in figure 8 the mean inclination angle $\theta$ of the bubble pair centreline with respect to the vertical direction (see sketch in the figure). (It should be noted that, since our measurements are only quasi-three-dimensional, $\theta =0$ does not necessarily mean that the bubbles are in line because they may nevertheless be separated in the depth direction by up to a distance $2d_b$.) For all cases an almost uniform distribution of $\theta$ is observed, i.e. there is no preferential alignment for a bubble pair. This is consistent with visual inspection of the experimental images (see supplementary movies) which show that the bubble pair orientations are not persistent, but instead the bubbles continually trade places in a ‘leapfrog’ fashion. This observation was also found in many 3-D experiments (Stewart 1995; Riboux, Risso & Legendre 2010) and direct numerical simulation for bubble swarms (Bunner & Tryggvason 2003; Esmaeeli & Tryggvason 2005). It is, however, strikingly different from the behaviour observed for isolated bubble pairs where a stable configuration is observed for two clean spherical bubbles where they rise side by side (Hallez & Legendre 2011), while for contaminated systems their stable configuration is in line due to the lift reversal experienced by the trailing bubble (Atasi et al. 2023). From the standpoint of our experiment, several possible reasons might explain this discrepancy. First, compared with the case of spherical bubbles, a pair of deformable bubbles that are hydrodynamically interacting are likely to generate an asymmetric flow which could destabilize their configuration. Second, in our experiments the bubbles have oscillating and/or chaotic rising paths, for which the probability that two bubbles will rise in a stable arrangement is very low. By contrast, in Hallez & Legendre (2011) the bubbles are fixed at various positions, and in Atasi et al. (2023) $Re_b$ is small enough such that the bubbles have straight rising paths. Third, the bubble pairs in our experiments are not isolated. Given the present bubble sizes/void fractions, they can experience fluctuations and turbulence due to the wakes of other bubbles in the flow, and this will readily suppress any preferential orientation that might have occurred were the bubble pairs isolated.

Figure 8. Orientation of bubble pair for the different cases (e.g. in-line bubble pair for $\theta =0^\circ$ and side by side for $\theta =90^\circ$).

Although the bubble pair orientation is almost random, the impact of their interaction on $U_b$ will, however, depend on $\theta$, especially in the present bubble regime (deformable bubbles with $Re_b\sim O$(100–1000)). For example, in the side-by-side configuration the two bubbles are outside of each other's wakes and the modification to the drag force on each bubble is minimal (Kong et al. 2019). On the other hand, for the in-line configuration the trailing bubble is sheltered by the leading bubble and the reduced pressure behind the leading bubble causes the trailing bubble to be sucked towards it, increasing the rise velocity of the trailing bubble while the leading bubble is almost unaffected (Zhang et al. 2021). Consequently, two mechanisms are in competition: one is vortex interaction (for relatively large $\theta$) and the second is BIT, i.e. wake entrainment (for smaller $\theta$). These effects mean that only the rise velocity of trailing bubbles will be significantly affected by the clustering, and hence, when averaged over all orientations, the increased vertical velocity of the trailing bubbles leads to an overall increase in $U_b$. This explains the increased mean rise velocity for $N_b=2$ compared with the $N_b=1$ results in figure 7. The increase is, however, minimal for the cases SmTapLess and SmTap. This is most likely due to the bubble wakes being weaker for these cases, a result of which is that the bubble interaction and the associated effect on $U_b$ is also weaker.

The results in figure 7 show that $U_b$ further increases as $N_b$ is increased beyond 2. This can be understood in terms of the enhanced opportunity for bubbles to be sheltered by other bubbles as $N_b$ increases. However, the increase when going from $N_b=2$ to $N_b=3$ is not as strong as when going from $N_b=1$ to $N_b=2$ because the greatest effect of sheltering will occur when all the bubbles in the cluster are in line, and the probability of the bubbles in an $N_b=3$ cluster all being in line will be less than that for the bubbles in an $N_b=2$ cluster being in line. It is interesting to note that for experiments on heavy particles settling in a quiescent fluid, similar behaviour was also found, with clustered particles falling faster (Huisman et al. 2016). In that case, the enhanced settling velocity was also attributed to a sheltering effect, i.e. reduced drag on particles that are falling in the wake of other particles. However, in that context, the particle clusters were found to exhibit strong alignment with the vertical direction, unlike our bubble clusters whose orientations are almost random (at least for the $N_b=2$ case).

3.4. The PDF of $N_b$-bubble cluster rise velocity fluctuations

More detailed information about the bubble rise velocities is provided in figure 9, where PDFs of the fluctuating rise velocity $u_b=\tilde {u}_b-U_b$ (where $\tilde {u}_b$ is the instantaneous vertical bubble velocity) of bubbles in an $N_b$-bubble cluster are shown. Only results up to $N_b=3$ are shown since the results for $N_b>3$ are very noisy.

Figure 9. The PDFs of the fluctuating rise velocity of bubbles in an $N_b$-bubble cluster, normalized by their standard deviations $\sigma _{u_b}$, with $N_b = 1, 2, 3$ for all the six cases: (a) SmTapLess, (b) SmTap, (c) SmPen+, (d) LaTapLess, (e) LaTap and (f) LaPen+.

The PDFs for all cases are asymmetric, and in particular are negatively skewed, especially for the three smaller bubble cases. This is consistent with the observations in Martínez et al. (2010) for a bubble swarm rising in still water, and also with a recent study for a single bubble rising in a flow with background turbulence (Ruth et al. 2021). The latter study is relevant to our experiments since any bubble rising in a swarm experiences a turbulent flow due to the turbulence generated by the surrounding bubbles. Other previous studies, however, show contrasting results for the PDF shapes. Riboux et al. (2010) used optical probes to measure the bubble velocity and found its vertical fluctuation was positively skewed. They argued the vertical bubble velocity fluctuations are controlled by the vertical liquid velocity fluctuations which have a positive skewness. In our previous study (Ma et al. 2022) we also found that, for dispersed bubbly flows, the vertical liquid velocity fluctuations were positively skewed, but here we find the opposite skewness for the bubble velocities. Loisy & Naso (2017) simulated a single bubble interacting with homogeneous isotropic turbulence and found that the PDF of the bubble vertical velocity fluctuations was approximately Gaussian. They argued that a skewness of either sign in the PDF of the bubble vertical velocity fluctuations is possible, depending in a complex manner on the background turbulence and presumably on other factors such as the bubble size compared with the turbulence length scales. For example, a trailing bubble located in the wake of a leading bubble would predominantly experience a positive upward liquid velocity due to the wake entertainment effect. On the other hand, the counterflow generated between bubbles leads to a downwards liquid motion that balances the flow being entrained into the bubble wakes. If other bubbles are in these downward regions it could lead to an enhancement of their negative fluctuating velocity. The skewness of the bubble rise velocity fluctuations therefore depends in a complex way on how the bubbles interact with the flow field generated by other bubbles in their vicinity.

4. Conclusions

We conducted experiments on the temporal evolution of bubble clusters with the aid of a new bubble tracking method for crowded swarms. Our results show that 2-bubble clusters are the most common, however, 3- and 4-bubble clusters also often occur. The clusters persist on average for a time of order $d_b/U_b$, although rare clusters persisting for an order of magnitude longer are also observed. Furthermore, surfactants are observed to enhance the cluster sizes and their lifetimes. A positive correlation between cluster size and bubble rise speed is observed, with clustered bubbles rising up to $20\,\%$ faster than unclustered bubbles. We investigated the mean inclination angle $\theta$ of a bubble pair and found that, for all cases, $\theta$ has an almost uniform distribution. However, the impact of the inclination on the bubble rise velocities is different for different $\theta$. For instance, the in-line configuration is expected to have a higher impact on the bubble rise velocity compared with the side-by-side configuration.

We also investigated the PDFs of the fluctuating rise velocities of bubbles in clusters of varying sizes. The PDFs for all cases are asymmetric, and in particular are negatively skewed, especially for the three smaller bubble cases.

Finally, while our cluster tracking method is only quasi-three-dimensional, a fully 3-D method for dense, deformable bubbles can be developed by combining our bubble identification method with the recent tracking algorithm of Tan, Zhong & Ni (2023) that currently applies to spherical bubbles.