2. Three-dimensional bubble tracking at high image densities

Figure 2. Experimental facility and 3-D Lagrangian bubble tracking. (a) Schematic of the experimental set-up, (b) schematic of the test section including the jet array and the bubble generator, (c) the second-order longitudinal structure function of turbulence in the test section, (d) trajectories of a colliding bubble pair coloured by velocity magnitude in the presence of background bubbles tracked at the same time and (e) consecutive bubble images before and after the moment of coalescence along with the 3-D reconstruction of the two bubbles and the merging interface in between.

Bubble collision and coalescence can be affected by both turbulence and buoyancy. To simplify the problem and focus on the turbulence effect, we constructed a unique setup that features intense turbulence, as depicted in figures 2(a) and 2(b), known as V-ONSET (Masuk et al. Reference Masuk, Salibindla and Ni2019b ; Tan et al. Reference Tan, Xu, Qi and Ni2023a ). The jet array injects powerful jets, producing turbulence with a high Taylor Reynolds number, Re $_\lambda = u'\lambda /\nu$ = 493, where $u'$ is the root-mean-squared fluctuation velocity, the Taylor microscale $\lambda$ is $\sqrt {15\nu /\epsilon }u'$ and $\nu$ is the kinematic viscosity of the fluid, enabling a wide scale separation from the Kolmogorov scale $\eta$ of 0.07 mm to the integral scale $L$ of 110 mm. The turbulent kinetic energy is $k=0.29$ m $^2$ s⁻ $^2$ . Figure 2(c) illustrates the second-order longitudinal structure function, $D_{LL}(r)$ , which displays a distinct scaling law of $r^{2/3}$ , aligning with the well-established Kolmogorov theory for fully developed turbulence (Kolmogorov Reference Kolmogorov1949). The turbulence energy dissipation rate estimated from the structure function is 0.03 m $^2$ s⁻ $^3$ and the Taylor microscale is calculated to be 2.9 mm. Bubbles were produced by an ultra-fine bubble generator at the bottom of the octagonal test section. The bubble size distribution is presented in figure 3(a) and the diameter ranges from 0.2 to 2 mm ( $3\eta-30\eta$ ). The rise velocity for bubbles in this range is approximately 0.1 m s⁻¹ (Salibindla et al. Reference Salibindla, Masuk, Tan and Ni2020), and its ratio to $u'$ is around 1, yielding a bubble Reynolds number of the order of 100. The void fraction, $\alpha$ , in our experiments is approximately 0.04 %. The detailed setup parameters and flow control specifications are provided in Appendix A.

Figure 3. (a) Bubble size distribution for all bubbles, and only those in collision and coalescence events. The shaded probability density function (PDF) is fitted with a log-normal distribution function for better illustration. (b) The collision kernel representing the clustering effect as a function of two Stokes numbers based on two different time scales, the $\textit{St}_\eta$ Kolmogorov scale and $\textit{St}_{d_b}$ based on the bubble-sized eddies.

At the moderate Reynolds number, the Strouhal number, Str, is approximately 0.18, leading to a vortex shedding time scale of $\tau _f=d_b/(v_b\textit{Str} )\approx 55$ ms. In comparison, the eddy turnover time at the same length scale can be estimated as $\tau _{d_b}=d_b^{2/3}/\epsilon ^{1/3}\approx 3.2$ ms, which is an order of magnitude smaller than the bubble vortex shedding time scale. This suggests that the wake structures generated by the bubbles are rapidly convected by the surrounding turbulent eddies, and the bubble motion is mostly governed by the turbulence. Furthermore, the rate of potential energy injection to turbulence by bubble buoyancy, given by $\alpha g w_b$ (where $w_b$ is the bubble rise velocity), is estimated to be approximately 0.0003 m $^2$ s⁻ $^3$ as the void fraction $\alpha \approx 0.04\,\%$ is very small. In comparison, the energy dissipation rate is 0.03 m $^2$ s⁻ $^3$ , much larger than the buoyancy effect, suggesting the modulation by bubbles to turbulent background flows is negligible.

Four cameras positioned at the periphery of the test section, as illustrated in figure 2(a), are used to capture the bubble motion. As the bubble concentration increases to acquire sufficient collision and coalescence events, the occurrence of overlapping bubble images increases, leading to challenges in tracking bubbles pairs in close proximity (Katz & Sheng Reference Katz and Sheng2010; Xue, Qu & Wu Reference Xue, Qu and Wu2013; Lebon et al. Reference Lebon, Perret, Coëtmellec, Godard, Gréhan, Lebrun and Brossard2016; Mathai et al. Reference Mathai, Huisman, Sun, Lohse and Bourgoin2018; Masuk et al. Reference Masuk, Salibindla and Ni2019a ; Salibindla et al. Reference Salibindla, Masuk, Tan and Ni2020; Shao et al. Reference Shao, Mallery, Kumar and Hong2020; Song et al. Reference Song, Qian, Zhang, Yin and Wang2022; Wu et al. Reference Wu, Zhang, Wu and Cen2020). To address this challenge, a new method is developed to utilise temporal information (bubble images may not overlap prior to the moment when they do) to predict the bubble location and then cross-correlated the predicted images against the captured ones to identify the most probable location of the bubble (a brief overview of this new method is provided in Appendix B). The tracked bubble locations correspond to the centroids of the bubbles, as those within the specified size range are spherical. Velocity and acceleration are determined by convolving the position data with a Gaussian kernel (Mordant, Crawford & Bodenschatz Reference Mordant, Crawford and Bodenschatz2004; Ni, Huang & Xia Reference Ni, Huang and Xia2012) to eliminate noise.

By applying this method, we were able to track bubbles with high image densities with up to 87 % of the area overlapped with neighbours (Tan et al. Reference Tan, Zhong and Ni2023b ). Figure 2(d) shows trajectories, coloured by velocity magnitude, of a sample pair of bubbles that collide with each other, in the presence of a high concentration of background bubbles. The time series of the coalescence moment is shown in figure 2(e). The reconstructed position of bubbles can achieve a high level of precision at around 10 $\unicode{x03BC} \textrm {m}$ , reaching sub-pixel accuracy through the utilisation of bilinear interpolation (Press et al. Reference Press, Teukolsky, Vetterling and Flannery2007) in close proximity to the peak of the image cross-correlation. This technique allows us to determine the bubble pairs when they are as close as 20 $\unicode{x03BC} \textrm {m}$ , even allowing for 3-D reconstruction of the merging interface, as shown in figure 2(e). So the separation that can be resolved is close to four orders of magnitude from 20 $\unicode{x03BC} \textrm {m}$ to 5 cm.

From the experimental images of bubble motion, we identified around 3.85 $\times\, 10^6$ individual bubbles and their corresponding trajectories. Any two bubbles appearing simultaneously within the view volume were considered a potential bubble pair, resulting in approximately 7 $\times\, 10^8$ pairs in our datasets, which is sufficient for calculating the bubble structure function and approach velocity of collision, as detailed in the following section.

3. Size distribution and size ratio

The bubble size distribution, computed from $3.85\times 10^6$ tracked bubbles, is presented as a black solid line in figure 3(a), spanning from $3\eta$ to $30\eta$ with a peak located at $7\eta$ . Based on the trajectories of these bubbles, collision events were identified when the film thickness between two bubbles, $|r - d_b|$ , fell below $0.2d_b$ (Gong et al. Reference Gong, Gao, Han and Luo2020). Finally, only those events when the film thickness was below $10\, \unicode{x03BC} \textrm {m}$ were checked manually against experimental images to identify bubble pairs that ultimately coalesce.

The representative size of bubble pairs involved in bubble collisions and coalescence is characterised by the equivalent bubble diameter, defined as $d_{{eq}} = 2 d_{b1} d_{b2} / (d_{b1} + d_{b2})$ . As we will show later, the collision and coalescence are dominated by events with two similar-sized bubbles. For simplicity, $d_b$ is used in place of $d_{{eq}}$ wherever $d_{{eq}}$ is intended. The size distribution of $d_b$ of bubble pairs undergoing collision events is shown in figure 3(a) as blue dots. The blue curve with the shaded area represents the log-normal fit to the PDF. It is evident that the PDF for collision events peaks at a larger size compared with the overall bubble distribution. This shift may be attributed to preferential concentration, where bubbles of certain sizes tend to cluster due to inertial effects, leading to locally elevated concentrations in turbulence – especially in the intermediate Stokes number range (Rensen, Luther & Lohse Reference Rensen, Luther and Lohse2005; Calzavarini et al. Reference Calzavarini, Van den, Thomas, Toschi and Lohse2008). To quantify this clustering behaviour, we introduce a size-dependent clustering function $ g_{\textit{CL}}(d_b)$ , which describes how the likelihood of collisions between bubbles of diameter $ d_b$ deviates from what would be expected under a random, uniform distribution. Specifically, we relate the measured size distribution of collided bubble pairs, $ f_d^c(d_b)$ , to the overall bubble size distribution $ f_d(d_b)$ via

(3.1) \begin{equation} f_d^c(d_b) = \frac {f_d^2(d_b) \, g_{\textit{CL}}(d_b)}{\int _0^\infty f_d^2(d_b) \, g_{\textit{CL}}(d_b) \, \mathrm{d}d_b}. \end{equation}

Here, $ f_d^2(d_b)$ is the expected distribution of randomly selected bubble pairs, and $ g_{\textit{CL}}(d_b)$ accounts for the enhanced probability of collisions due to clustering. From (3.1), $ g_{\textit{CL}}(d_b)$ can be directly computed using the measured distributions $ f_d(d_b)$ (black solid line) and $ f_d^c(d_b)$ (blue) shown in figure 3(a). The resulting values are plotted as blue dots in figure 3(b). In addition, the blue solid line represents a smoothed curve of $g_{\textit{CL}}$ by using the fitted PDF of collision events. Since the Stokes number is widely used in the analysis of the clustering effect (Loth Reference Loth2000), the $x$ -axis is presented in terms of the Stokes number, St $=d_b^2/36\nu \tau _{fb}$ , where $\tau _{fb}$ represents the characteristic time scale of the flows around bubbles. There are two possible time scales for $\tau _f$ , the Kolmogorov time scale, $\tau _\eta =(\nu /\epsilon )^{1/2}$ , which is used in the study by Mathai et al. (Reference Mathai, Calzavarini, Brons, Sun and Lohse2016) for micro-bubbles, and the bubble-sized eddy turnover time, $\tau _{d_b}=d_b^{2/3}/\epsilon ^{1/3}$ . The Stokes numbers based on both flow time scales, $\textit{St}_\eta$ and $\textit{St}_{d_b}$ , are also shown as the top and bottom $x$ -axes in figure 3(b), respectively.

Interestingly, instead of $\textit{St}_\eta$ , the curve of $g_{\textit{CL}}(\textit{St}_{d_b})$ exhibits a peak near $\textit{St}_{d_b}=1$ , which suggests that $\textit{St}_{d_b}$ is more appropriate to characterise the preferential concentration effect of bubbles in turbulence. The dimensionless function $g_{\textit{CL}}(\textit{St}_{d_b})$ expressed in terms of $\textit{St}_{d_b}$ should be applicable to other bubbly flows, provided that bubble–turbulence modulation remains negligible, as discussed in § 2.

By using $g_{\textit{CL}}$ , the bubble collision frequency in turbulence can be estimated as follows:

(3.2) \begin{equation} \mathcal{F}^T_c(d_b)=n^2f_d^c(d_b)\mathrm{d}d_b\pi d_b^2\delta v_c^\parallel (d_b) ,\end{equation}

where $n$ is the number density of all bubbles in the flows. The parameter $f_d^c(d_b)\mathrm{d}d_b$ represents the probability of colliding bubble pairs with size $d_b\pm \mathrm{d}d_b$ , and $f_d^c(d_b)$ is obtained through (3.1) in terms of $g_{\textit{CL}}$ and $f_d(d_b)$ .

The PDF for coalescence is also shown as yellow symbols in figure 3(a), which is considerably narrower than that for collision events, indicating that the efficiency of coalescence is probably also size dependent. It will be further examined in § 6.

Figure 4. The PDF of the bubble size ratio $\xi$ for bubble pairs with sizes $d_{b1}$ and $d_{b2}$ , respectively, in collision and coalescence events, as well as distant bubble pairs.

In addition to the size distribution, since collision and coalescence may occur between bubble pairs of different sizes, the size ratio, defined as $\xi =d_{b1}/d_{b2}$ , was also computed. The resulting distributions are shown as blue and yellow histograms in figure 4. It is evident that both collision and coalescence are dominated by similar-sized pairs. To ensure that this is not a result of the size distribution alone, a PDF of bubble pairs that are at 80 $\eta$ away is shown as green dots, which shows a flatter distribution, indicating more uneven sized pairs. At this far distance, the dynamics of the bubbles becomes less correlated and the distribution becomes almost random. If we assume the bubble locations are completely random, the probability of finding a pair with a certain size ratio based solely on the bubble size distribution, is represented by the green solid line in figure 4, which is close, although not perfect, to the experimental measurement, suggesting that these pairs are close to random at large enough separations.

The bias of the PDF for collisions from that of distant bubble pairs is likely due to the preferential concentration effect in turbulence, where bubbles of similar sizes tend to cluster within eddies of similar size, increasing their likelihood of collision. The distribution for coalescence events (yellow histogram) closely follows that of collisions but is slightly more constrained, suggesting that coalescence is even more probable when the size ratio is closer to unity.

To be more quantitative, 98 % of collision events occur within $\xi \lt 2$ . Consequently, in the subsequent analysis, the discussion focuses on the collision and coalescence of bubble pairs with common size.

4. Bubble structure function and bubble approach velocity

For $7\times 10^8$ bubble pairs, we project their velocity differences in the direction of their separation, $\delta _r v^\parallel (r)=v_\parallel (\boldsymbol {x}+\boldsymbol {r}) - v_\parallel (\boldsymbol {x})$ . The collision velocity is defined as $\delta _r v^\parallel (r)$ at $r=1.2d_b$ , which is the separation when the film between two bubbles starts to be drained (Gong et al. Reference Gong, Gao, Han and Luo2020). Nevertheless, to fully understand the entire bubble dynamics, in this section, we extend the relative velocity to a wide range of separations, which is similar to the structure function studied in single-phase turbulence.

The second-order bubble structure function is defined as $D^B_{LL}(r) \equiv \langle \delta _r v_\parallel ^2 \rangle$ . The measured bubble structure functions are presented in figure 5(a). Bubbles are divided into four groups based on their size. Each curve in figure 5(a) represents the statistics of bubble pairs in the same bubble group.

Two distinct scalings can be observed. When two bubbles are sufficiently far apart, $r\gt 2d_b$ , the structure function $D^B_{LL}$ follows $r^{2/3}$ , similar to the inertial-range scaling observed in single-phase turbulence. This similarity arises because the motion between bubbles with a large separation is predominantly influenced by large-scale turbulent transport. As the separation between bubbles decreases below $2d_b$ , $D^B_{LL}$ , transitions to a scaling of 2. Unlike the Kolmogorov scale $\eta$ , which marks the transition between the dissipative and inertial ranges, the bubble structure function transitions between two scalings at $2d_b$ . This transition scaling with bubble size suggests that the $r^2$ scaling is not due to the turbulence’s dissipative range but rather to the linear dependence of the bubble relative velocity on separation. A plausible explanation for this linear dependence is the entrainment of bubble pairs by a common vortex, causing their relative acceleration to scale as $\delta a \sim r$ . The value of $\delta a$ can be expressed in terms of $\delta v$ as $\delta a=\delta v \text{d}\delta v/\text{d}r$ , which finally leads to $\delta v^2\sim r^2$ .

Figure 5. The second-order longitudinal bubble velocity structure function between bubbles with similar size, normalised (a) with the turbulent fluctuation velocity $u'$ and (b) with bubble-scaled kinematic energy $(d_b/\lambda )^{2/3}(\epsilon d_b)^{2/3}$ , versus the bubble separation normalised (a) with the Kolmogorov scale and the integral length scale and (b) with the bubble diameter.

More importantly, a noticeable observation is the upward shift of the curves as the bubble size increases, despite the fact that the structure functions are normalised by turbulence scales, which indicates that the bubble dynamics does not simply mirror that of the surrounding eddies, as incorrectly assumed in many models. The key missing piece is the biased sampling, i.e. bubbles preferentially get entrained in eddies of similar sizes (Wang et al. Reference Wang, de, Xander and Toschi2024). Subsequently, by further non-dimensionalising $D^B_{LL}$ with the eddy energy that scales as $(d_b/\lambda )^{2/3}$ , we obtain the normalised plot presented in figure 5(b). The curves exhibit a reasonable collapse in both limits, i.e. $r^2$ and $r^{2/3}$ , providing support to the biased sampling argument. By utilising the collapse observed in figure 5(b), we can readily formulate the bubble structure function as follows:

(4.1) \begin{equation} D^B_{LL}(r, d_b)= \begin{cases} C_0(d_b/\lambda )^{2/3}(\epsilon d_b)^{2/3}(r/d_b)^2 & \text{for $r\leqslant 2d_b$},\\ C_1(d_b/\lambda )^{2/3}(\epsilon d_b)^{2/3}(r/d_b)^{2/3} & \text{for $2d_b\lt r\ll L$}, \end{cases} \end{equation}

where $C_0$ and $C_1$ are constants. From the data plotted in figure 5, $C_0$ and $C_1$ are estimated to be 1.73 and 4.36, respectively, both of which are the order of unity, suggesting that the additional term of biased sampling is reasonable.

Equation (4.1) readily provides the estimation of the approach velocity which is defined as the relative bubble velocity before the lubrication pressure slows them down, which begins at around $r=1.2d_b$ (Gong et al. Reference Gong, Gao, Han and Luo2020). The approach velocity, $\delta v_c^\parallel$ , can then be determined by setting $r=1.2d_b$ in the bubble structure function

(4.2) \begin{equation} \big\langle \delta v_c^\parallel (d_b) \big\rangle \equiv \Big [\big\langle \delta _r v^2_\parallel (r=1.2d_b) \big\rangle \Big ]^{1/2} = 1.2C_0^{1/2}(\epsilon / \lambda )^{1/3}d_b^{2/3} .\end{equation}

Equation (4.2) clearly shows that the approach velocity, instead of following $d_b^{1/3}$ as proposed in the classical framework, scales as $d_b^{2/3}$ . The additional 1/3 scaling highlights the missing biased sampling effect that was overlooked before.

Note that the approach velocity refers to the relative velocity between two bubbles at a particular separation scale. In most classical models – and in our work here – it is taken to be the relative velocity at $r = 1.2d_b$ , i.e. when the bubbles are nearly in contact. However, an alternative choice is to define it at larger separations, such as at the mean inter-bubble distance determined by the local bubble concentration. In the idealised case where particles move ballistically – with constant velocity magnitude and direction – the relative velocity would be the same regardless of the separation $r$ at which it is measured, and the two definitions would yield identical results. However, this assumption does not hold in turbulent flows, where the relative velocity between bubbles varies with separation $r$ , as shown in figure 5. In turbulence, the approach velocity should therefore be interpreted as a weighted average of the relative velocity over a range of separations – from the mean separation (approximately $12d_b$ ) down to contact. Since this range lies within the inertial subrange, the averaging does not affect the scaling with respect to bubble diameter $d_b$ , although it introduces a prefactor. This prefactor depends on a weighting function that reflects the probability that two bubbles remain on a collision course over a given separation $r$ . This probability is near one at small separations $r \approx d_b$ , and rapidly decays at larger values of $r$ , since the likelihood of continued approach over long distances becomes negligible. Consequently, evaluating the approach velocity at $r = 1.2d_b$ provides a good approximation, as the weighting function gives higher importance to bubble pairs that are already close to contact. This makes the choice of selecting relative velocity at $r=1.2d_b$ both physically reasonable and consistent with prior modelling approaches. Future work is needed to determine the exact correction prefactor introduced by the weighting function, in order to improve quantitative accuracy.

5. Contact time

Once the bubble pairs get sufficiently close and the film begins to drain, the contact time can be calculated based on two different methods, as depicted in figure 6. Figures 6(a), 6(b) and 6(c) show the trajectories, distance $r$ and longitudinal relative acceleration $\delta _r a_\parallel$ of a bubble pair with bubble diameters of 0.44 and 0.46 mm, respectively. One definition of contact time, $t^d_{ct}$ , is based on the distance, which measures the duration between bubbles when the gap between the two interfaces falls within 0.2 $d_b$ . Another definition of the contact time, denoted as $t^a_{ct}$ , focuses on measuring the duration over which the longitudinal relative acceleration $\delta _r a_\parallel$ becomes positive, indicating a repulsive force between the two bubbles from the lubrication pressure, until $r$ reaches the minimum.

Figure 6. Bubble contact time. (a) Schematic of trajectories of a bubble pair through collision, (b–c) the time series of the separation between these two bubbles and their relative acceleration along with the extracted contact times $\tau _d$ and $\tau _a$ , (d) the contact time normalised by the Kolmogorov time scale as a function of the bubble diameter and (e) the relative longitudinal acceleration as a function of the bubble diameter.

The contact time results based on both the distance ( $t^d_{ct}$ ) and the longitudinal relative acceleration ( $t^a_{ct}$ ) are presented in figure 6(d). The contact time determined by either criterion shows no clear trend with the bubble size. These results are clearly different from the prediction either based on the dimensional analysis, i.e. $d_b^{2/3}/\epsilon ^{1/3}$ , or another relation, $(\rho d_b^3/\sigma )^{1/2}$ , where $\rho$ and $\sigma$ are the fluid density and bubble surface tension, derived by Chesters (Reference Chesters1991) through the balance between the surface energy and kinetic energy. This difference suggests that larger bubble sizes do not necessarily result in a longer contact time or increased chances for coalescence in turbulence. In the analysis of coalescence efficiency in the following section, $t^a_{ct}$ is adopted as the representative contact time, as it is more directly related to the drainage process.

To explain this difference, an alternative way to estimate the contact time is through $\langle \delta _r v_\parallel \rangle / \langle \delta _r a_\parallel \rangle$ , which is the time scale that it takes for the approach velocity to decelerate to zero. With both bubbles tracked simultaneously, the relative acceleration between two bubbles at a given separation can be computed. The separation $r=1.1 d_b$ , or film thickness $0.1 d_b$ , is chosen, which is the midpoint between the start of the film drainage ( $r=1.2d_b$ ) to the point of contact $r=d_b$ . The dependence of $\langle \delta _r a_\parallel (r=1.2d_b)\rangle$ on bubble size is illustrated in figure 6(e). It is evident that $\langle \delta _r a_\parallel \rangle$ scales with $d_b^{2/3}$ , suggesting that larger bubble pairs experience stronger repulsive forces. The large repulsive force indicates a high pressure within the film which likely results in local deformation for larger bubbles. With the access to both $\langle \delta _r v_\parallel \rangle$ and $\langle \delta _r a_\parallel \rangle$ experimentally, it is then straightforward to derive the relationship

(5.1) \begin{equation} \tau _{ct} \approx \big\langle \delta _r v_\parallel \big\rangle \big/ \big\langle \delta _r a_\parallel \big\rangle \propto r^{2/3}/r^{2/3} \sim r^0 ,\end{equation}

which explains the finding in figure 6(d) that the contact time remains nearly constant regardless of the variation in bubble size.

So far, we have measured the critical parameters used in traditional models to predict collision and coalescence. The scaling of the collision velocity and contact time with the bubble size was hypothesised to be $d_b^{1/3}$ and $d_b^{2/3}$ , respectively, in traditional models (Coulaloglou Reference Coulaloglou1975; Lee et al. Reference Lee, Erickson and Glasgow1987; Prince & Blanch Reference Prince and Blanch1990a ; Luo Reference Luo1995). However, we have shown that these hypotheses are not correct. The corrected scalings should be $d_b^{2/3}$ and $d_b^0$ . Such deviation stems from the additional dimensionless number that accounts for the size ratio $d_b/\lambda$ between bubble and turbulence, which measures the effect of biased sampling that was missed in the models.

6. Coalescence rate in turbulence

Figure 7. The coalescence efficiency and the PDF of collision velocity. (a) The PDF of the coalescence and collision velocity calculated from all collision and coalescence events, along with the ratio in between which is defined as the coalescence efficiency (red symbols). Our model is plotted with a red solid line, to be compared with the one by Coulaloglou (Reference Coulaloglou1975) plotted with a blue solid line. (b) The PDF of $\delta v_c^\parallel$ for different bubble sizes normalised by the velocity scale $\langle \delta v_c^\parallel \rangle$ predicted by (4.2). They all collapse and can be modelled based on the log-normal distribution of the turbulence energy dissipation rate $f_\epsilon$ (black solid line).

The PDFs of the collision velocity normalised by the model prediction $\langle \delta v_c^\parallel \rangle$ in (4.2) for both collision and coalescence events are computed and shown as coloured bars in figure 7(a). Although the PDFs appear similar, the peak of the coalescence PDF is shifted towards a higher relative velocity than that of the collision PDF. The bubble coalescence efficiency, $\mathcal{E}_{cc}^T$ , is calculated by dividing the number of coalescence events ( $N_{cc}$ , the subscript $cc$ represents ‘collision then coalescence’) by the number of collision events ( $N_{c}$ ), i.e. $\mathcal{E}_{cc}^T=N_{cc}/N_c$ . The value of $\mathcal{E}_{cc}^T$ is plotted as red dots in figure 7(a) as a function of the collision velocity $\delta v_c^\parallel$ ; $\mathcal{E}_{cc}^T(\delta v_c^\parallel )$ exhibits a peak at an intermediate collision velocity, a Goldilocks zone for coalescence, approaching zero for both small and large collision velocities.

This Goldilocks zone for coalescence efficiency contrasts sharply with the well-adopted traditional model (Coulaloglou Reference Coulaloglou1975; Chesters Reference Chesters1991; Doubliez Reference Doubliez1991; Duineveld Reference Duineveld1996; Lehr & Mewes Reference Lehr and Mewes2001; Lehr et al. Reference Lehr, Millies and Mewes2002), which is shown as the blue line in figure 7. The traditional model is expressed as $\mathcal{E}_{cc}^T(\delta v_c^\parallel )=\text{exp}(-t_{\textit{dr}}/t_{ct})$ (Ross Reference Ross1971), where $t_{ct}\approx 1.7$ ms is taken from our experimental result (figure 6 d) and the drainage time model proposed by Chesters (Reference Chesters1975, Reference Chesters1991) and Luo (Reference Luo1995) is adopted:

(6.1) \begin{equation} t_{\textit{dr}}=C_{\textit{dr}}\rho _c \delta v_c^\parallel d_b^2/\sigma ,\end{equation}

where $\rho _c$ represents the density of the carrier fluid and $\sigma$ the surface tension of bubbles. The prefactor $C_{\textit{dr}}$ is introduced in the model to account for detailed interfacial properties, such as bubble surface mobility, the onset time of the drainage process and surfactant effects. In our study, we tried our best to use reverse osmosis water but the surfactant concentration in the entire tunnel may not be perfectly zero so we determined $C_{\textit{dr}}$ by fitting the coalescence efficiency model $\mathcal{E}^T_{cc}(\delta v_c^\parallel )$ to our experimental measurements, yielding a value of approximately 7.5. The model prediction is shown as a blue solid line. Although the prediction clearly departs from the measured $\mathcal{E}_{cc}^T(\delta v_c^\parallel )$ , the classical model seems to converge with the data at large collision velocities, i.e. $\delta v_c^\parallel \gt 5\langle \delta v_c^\parallel \rangle$ . At low collision velocities, however, the model predicts nearly 100 % coalescence efficiency yet the data suggest that the coalescence efficiency should continue to decline down to zero.

This mismatch resides in the assumption of the model that the collision velocity is persistent over $t_{\textit{dr}}$ for the liquid film to be drained. In turbulence, however, the relative velocity between bubbles is affected by the surrounding eddies carrying their own characteristic scales. If the collision velocity is smaller than the eddy velocity, the collision takes too long, longer than the decorrelation time of the flow. In this case, the drainage process is interrupted by the chaotic motion of surrounding eddies in random directions, causing bubble pairs to be driven apart.

As a result, the coalescence efficiency has to account for an additional probability, $\mathcal{P}$ , of the bubble collision velocity $\delta v_{c}^\parallel$ being large enough to keep bubble pairs together for coalescence to occur before eddies with velocity scale of $|u_l|$ separate them. If $|u_l| \ll \delta v_{c}^\parallel$ , $\mathcal{P}$ should approach one since the eddy velocity is too slow to interrupt coalescence and the process proceeds as if the surrounding flow is frozen. At the opposite limit when $|u_l| \gg \delta v_{c}^\parallel$ , $\mathcal{P}$ approaches zero, suggesting that the coalescence cannot happen without sufficient time to drain the film, as is the case for the left tail of curve for $\mathcal{E}^T_{cc}$ in figure 7(a). Between these two limits, there should exist a critical velocity, $\kappa \delta v_c^\parallel$ . Then $\mathcal{P}$ can be expressed as

(6.2) \begin{equation} \mathcal{P} \big(|u_l| \lt \kappa \delta v_{c}^\parallel \big)=\int _{-\kappa \delta v_c^\parallel }^{\kappa \delta v_c^\parallel } f_u(u_l)\text{d} u_l ,\end{equation}

where $f_u(u_l)$ is the PDF of the eddy velocity at the bubble scale which can be approximated by a Gaussian distribution with zero mean and a deviation $C_2 (\epsilon d_b)^{2/3}$ since $d_b$ is in the inertial range (Kolmogorov Reference Kolmogorov1949). The final model of $\mathcal{E}_{cc}^T$ by including $\mathcal{P}$ and the Ross (Reference Ross1971) model is

(6.3) \begin{equation} \mathcal{E}_{cc}^T \big(\delta v_c^\parallel \big) = \mathcal{P} \big(|u_l| \lt \kappa \delta v_{c}^\parallel \big) \text{exp} \big({-}C_{\textit{dr}}\rho _c \delta v_c^\parallel d_b^2/\sigma t_{ct} \big) .\end{equation}

The model prediction using the mean bubble diameter from experiments $d_b=12\eta$ is shown as a red curve in figure 7(a), which matches well with the experimental results using the fitting parameters $\kappa = 0.14$ constrained mostly by the left tail and $C_{\textit{dr}}=7.5$ constrained by the right tail.

Since the coalescence efficiency shows non-trivial dependence on the collision velocity, the PDF of the collision velocity in turbulence becomes important for modeling of coalescence rate, which is also calculated and plotted with symbols in figure 7(b) for three different groups of bubble sizes $d_b$ = 8.46 $\eta$ , 12.32 $\eta$ and 16.72 $\eta$ . The collision velocities are normalised by their respective mean based on (4.2), $v^*=\delta v_c^\parallel /\langle \delta v_c^\parallel \rangle$ . Since bubbles are transported through eddies across multiple scales, the PDF of the collision velocity can be derived based on the log-normal distribution of the normalised energy dissipation rate $\epsilon ^*=\epsilon _r/\langle \epsilon \rangle$ given by (Kolmogorov Reference Kolmogorov1962; Meneveau & Sreenivasan Reference Meneveau and Sreenivasan1987)

(6.4) \begin{equation} f_\epsilon (\epsilon ^*) = \frac {1}{\epsilon ^* \sqrt {2\pi \sigma ^2_{\textit{ln}\epsilon }}} \text{exp}\left [{-\frac {(\text{ln}\epsilon ^*-\mu _{\textit{ln}\epsilon })^2}{2\sigma ^2_{\textit{ln}\epsilon }}}\right]\! ,\end{equation}

where $\sigma ^2_{\textit{ln}\epsilon }=A+\mu \text{ln}(L/r)$ and $\mu _{\textit{ln}\epsilon }=-\sigma ^2_{\textit{ln}\epsilon }/2$ such that the mean of $\epsilon ^*$ given by this distribution is 1. The parameter $\mu =0.25$ is the intermittency exponent and $A$ is a parameter determined by the specific flow conditions. After normalisation, $v^*$ solely scales with $\epsilon ^{*1/3}$ , which can be modelled as $v^*=\zeta \epsilon ^{*1/3}$ , where $\zeta$ is a coefficient that is left out in the non-dimensionalisation. The PDF for $v^*$ can be readily derived as $f_v^c(v^*)=3v^{*2}f_\epsilon ((v^*/\zeta )^{3})/\zeta ^3$ , where $f_\epsilon$ is evaluated at $r=d_b$ . We determined both $\zeta$ and $A$ by fitting $f_v^c(v^*)$ to the PDF of $v^*$ from experiments, which are 2.8 and 1, respectively. The model of the PDF of $v^*$ based on (6.4) is plotted as a solid black line in figure 7(b), which shows a good agreement with experimental data. We also plotted the dimensional $\delta v_c^\parallel$ in the inset of figure 7(b) to confirm the agreement between the model and the measured PDF at different sizes.

Figure 8. The mean coalescence velocity $\langle |\delta v_{cc}^\parallel | \rangle$ and the size distribution of bubble pairs undergoing coalescence. (a) The mean $\langle |\delta v_{cc}^\parallel | \rangle$ calculated from the coalescence efficiency $\mathcal{E}_{cc}^T$ and the PDF of $\delta v_c^\parallel$ . Predictions by our model and the RCC model (Ross Reference Ross1971; Chesters Reference Chesters1975; Coulaloglou Reference Coulaloglou1975) are shown for comparison. The experiment measured $\langle |\delta v_{cc}^\parallel | \rangle$ is also plotted to validate our model. (b) The size distribution of bubble pairs involved in coalescence events predicted by our model and the RCC model.

With the derived PDF of the collision velocity and the coalescence efficiency, the PDF of the coalescence velocity can be formulated as $f_v^{cc}(\delta v_{cc}^\parallel )=\mathcal{E}_{cc}^T(\delta v_c^\parallel )f_v^c(\delta v_c^\parallel )/ \int _0^\infty \mathcal{E}_{cc}^T(\delta v_c^\parallel )f_v^c(\delta v_c^\parallel )\text{d}\delta v_c^\parallel$ , where $f_v^c(\delta v_c^\parallel )$ is the PDF of collision velocity, and $f_v^{cc}(\delta v_{cc}^\parallel )$ the PDF of collision velocity for coalescing bubble pairs, $\delta v_{cc}^\parallel$ . Then the average coalescence velocity can be calculated as $\langle \delta v_{cc}^\parallel \rangle = \int _0^\infty vf_v^{cc}(v)\text{d}v$ . In figure 8(a), the model-predicted $\langle \delta v_{cc}^\parallel \rangle$ is shown as a red solid line. It is important to note that the calculation of $\mathcal{E}_{cc}^T(\delta v_c^\parallel )$ from (6.3) uses a constant contact time for all bubble sizes, as demonstrated in figure 6. To validate our model prediction, we divided the coalescing bubble pairs into groups with four different sizes, and calculated their mean collision velocity and corresponding 95 % confidence intervals, plotted as scattered squares with error bars.

As shown in figure 8(a), our model prediction aligns reasonably with the experimental measurement. Note that symbols represent the statistics of only the coalesced bubbles and the line came from the integration of several quantities, including the new ‘Goldilocks effect’ of coalescence efficiency, the nearly constant contact time and the PDF of the collision velocity derived from the log-normal distribution of the turbulence energy dissipation rate. For comparison, the prediction by the traditional hypothesised model based on the previous studies (Ross Reference Ross1971; Chesters Reference Chesters1975; Coulaloglou Reference Coulaloglou1975), specifically the RCC model, where $\delta v_c^\parallel = [C_2(\epsilon d_b)^{2/3} ]^{1/2}$ and $t_{ct} = d_b^{2/3} /(C_2^{1/2}\epsilon ^{1/3})$ , is also calculated and plotted as a blue solid line. The resulting coalescence velocity is much lower than both our model prediction and experimental measurement because of the lack of bubble–eddy interactions in several key quantities.

Our model can be further validated by calculating the bubble size distribution for bubble pairs involved in coalescence events. The number of coalescence events, $N_{cc}(d_b)$ , for the given number of collision events, $N_{c}(d_b)$ , can be estimated as

(6.5) \begin{equation} N_{cc}(d_b) = \int _0^\infty N_c(d_b)\mathcal{E}^T_{cc} \big(\delta v_c^\parallel \big) f_v^c \big(\delta v_c^\parallel \big) \text{d}\delta v_c^\parallel .\end{equation}

Since $N_{c}(d_b)$ is not a function of $\delta v_c^\parallel$ , it can be factored out of the integral; $N_{c}(d_b)$ can be presented by the bubble size distribution of collision events $f_d^c(d_b)$ as $N_c(d_b)=Nf_d^c(d_b)\text{d}d_b$ , where $N$ is the total number of collision events and $f_d^c(d_b)$ is shown as the blue solid line in figure 3(a). Defining the integral as $\vartheta (d_b)=\int _0^\infty \mathcal{E}^T_{cc}(\delta v_c^\parallel ) f_v^c(\delta v_c^\parallel ) \text{d}\delta v_c^\parallel$ , the size distribution of coalescence events can be expressed as $f_d^{cc}(d_b) =f_d^c(d_b)\vartheta (d_b) / \int _0^\infty f_d^c(d_b)\vartheta (d_b) \text{d}d_b$ . Since we do not have a model for the size distribution of collision events, we fitted $f_d^c(d_b)$ by a log-normal distribution in figure 3(a) with a mean and deviation of 2.6 and 0.12, and focus on the effect of chosen model of $\mathcal{E}^T_{cc}$ , which can uses either the RCC model (blue line) or the revised Goldilocks zone model (red line in figure 8 b).

The RCC model clearly deviates from the experimental measurement. Our model shows good agreement with the experimental measurement shown with the red symbols, further demonstrating the validity of our model in predicting coalescence efficiency and collision velocity.