4.1. Turbulent kinetic energy and dissipation rate

The evolution of the turbulent velocity fields along the duct axis is critical for understanding bubble dynamics in developing flows. This spatial variation in turbulence directly governs the interplay between bubble breakup and coalescence, making its quantification essential for predicting bubble size evolution and optimising industrial multiphase flow systems (Li & Liao Reference Liao and Lucas2024). During the characterisation of our experimental set-up, near the inlet of the duct, we find that the incoming flow from the regenerative pump is nearly homogeneous, consistent with our large eddy simulation of the same pump (Javadi et al. Reference Kantarci, Borak and Ulgen2026). We find that the mean streamwise velocity shows an almost flat radial profile across most of the cross-section, followed by a sharp decline confined to the near-wall region. The r.m.s. turbulent fluctuation profiles exhibit a similar trend: fluctuations remain relatively uniform in the core and increase toward the wall, consistent with classical observations of developing turbulent duct and pipe flows (Tavoularis & Corrsin Reference Tomiyama1981; Pope Reference Prince and Blanch2000).

Figure 6.

Axial variation of $k$ along the duct centreline for three bulk velocities $V(Q)$ and three void fractions $\phi$ . Results are shown for (a) $V$ = 6.1 m s⁻¹, (b) $V$ = 7.4 m s⁻¹, (c) $V$ = 8.4 m s⁻¹.

Figure 7.

Axial variation of $\varepsilon$ along the duct centreline for three bulk velocities $V(Q)$ and three void fractions $\phi$ . Note the use of the log–log scale. The range of the coefficient of determination ( $R^2$ ) is provided. Results are shown for (a) $V$ = 6.1 m s⁻¹ ( $R^2$ : 0.94–0.97), (b) $V$ = 7.4 m s⁻¹ ( $R^2$ : 0.93–0.97), (c) $V$ = 8.4 m s⁻¹ ( $R^2$ : 0.93–0.95).

In all our experimental cases presented here, the $\mathcal{I}$ imparted by the pump rapidly decreases as the flow develops. Figures 6 and 7 present the axial evolution of the estimates of centreline normalised turbulent kinetic energy ( $k/V^2$ ) and normalised turbulent dissipation rate ( $\varepsilon {D} / \bar {V}^{3}$ ) with $\mathcal{L}$ at different bulk velocities ( $V$ ). The centreline k and $\varepsilon$ both decrease by an order of magnitude (over ${\sim}$ 90 % reduction) from the inlet to the farthest measurement station. This steep initial decay reflects the absence of sustaining turbulence production in the core flow immediately after the pump, so the injected turbulent eddies dissipate energy quickly. Beyond $\approx$ 22 hydraulic diameters downstream, the rate of decline becomes much gentler (almost levelling off for the lowest bulk velocity case), indicating the approach toward a quasi-equilibrium turbulence state. At the lowest bulk velocity $V=6.1$ m s⁻¹, for instance, $\varepsilon$ drops very rapidly in the first few hydraulic diameters and then flattens by $\mathcal{L}\approx 30$ (suggesting that the core turbulence has largely equilibrated by that distance). In contrast, at the highest velocity ( $V=8.4$ m s⁻¹), centreline dissipation is still gradually decreasing even at the farthest location ( $\mathcal{L}\approx 40$ ), implying a longer development length before fully developed conditions are reached. These observed trends – an initial exponential-like decay of turbulence followed by a slower power-law-like tail – are consistent with classical decaying turbulence behaviour downstream of intense turbulence sources, and they align with prior findings that high core turbulence decays rapidly in developing pipe flows until wall-driven production begins to dominate (Doherty et al. Reference Durst, Jovanović and Sender2007; Wilson & Smith Reference Wilson and Smith2007). The locally steeper decay of $\varepsilon$ for $\mathcal{L} \approx 12{-}18$ indicates a transition region where the inlet-imprinted turbulence rapidly adjusts to the duct confinement and developing shear, leading to a temporarily dissipation-dominated energy budget. Similar decay behaviour has been reported in other duct and pipe flows (Zikanov et al. Reference Zikanov, Krasnov, Boeck, Ziegler and Thess2019; Ding et al. Reference Du Cluzeau, Bois and Toutant2021).

Increasing $V$ significantly elevates both the centreline k and $\varepsilon$ across the board. Flows at higher Re inject more energy into turbulence, yielding larger initial k and accordingly, higher $\varepsilon$ values at the inlet and downstream. For instance, at $\mathcal{L}\approx 5$ the dissipation in the $V=8.4$ m s^–1 case is noticeably greater than in the $V=6.1$ m s⁻¹ case for the same air void fraction $\phi$ , and this gap persists downstream, although all cases eventually decay to significantly lower levels. Such a trend has been reported in other systems – e.g. Shawkat, Ching & Shoukri (Reference Shi and Magnaudet2008) observed that the turbulence dissipation in horizontal bubbly flows rises with increasing $V$ . Notably, the higher Re cases not only start with larger $\varepsilon$ near the inlet, but also retain a higher dissipation level further downstream (since a more energetic flow takes longer to dissipate). The fact that the $V=8.4$ m s⁻¹ data in figure 7(c) still show a slight downward slope at $\mathcal{L}=40$ , whereas the $V=6.1$ m s⁻¹ data in figure 7(a) have nearly plateaued by that distance, suggests that the development length to reach a fully developed turbulence profile grows with Re. This is consistent with single-phase pipe flow literature – higher inertia flows require longer distances for turbulence redistribution and stabilisation, underscoring the significant influence of Re on turbulence evolution in our experiments (Doherty et al. Reference Durst, Jovanović and Sender2007; Wilson & Smith Reference Wilson and Smith2007).

For the relatively low void fractions ( $\phi$ = 0.5–2 %) in the current study, the influence of $\phi$ on turbulence statistics is quite subtle. Any changes in the measured turbulent kinetic energy $k$ and dissipation rate $\varepsilon$ across this range remain within experimental uncertainty, in part because the PSV system lacks the fine-scale resolution to capture such small variations. Indeed, as shown in figures 6 and 7, the $k$ and $\varepsilon$ profiles at these low void fractions $\phi$ nearly overlap. This suggests that a threshold void fraction is required before bubbles appreciably augment the turbulence – a notion supported by prior studies. Lance & Bataille (Reference Lee1991) identified a regime at low void fraction where bubble–turbulence interactions are negligible, and Serizawa et al. (Reference Shi1975) likewise observed minimal turbulence modification at void fractions of the order of 1 %. Similarly, Shawkat et al. (Reference Shi and Magnaudet2008) reported that at void levels of 1 %, any bubble-induced changes in turbulence are minor and often within the measurement scatter. Under the dilute conditions of the present study, therefore, the net impact of $\phi$ on $k$ and $\varepsilon$ is modest, and subtle bubble-induced enhancements cannot be reliably distinguished from experimental noise.

At the duct inlet, where the flow emerges directly from the pump, the turbulent kinetic energy imparted by the pump is initially distributed relatively uniformly across the pipe cross-section, resulting in elevated turbulent dissipation rates both at the centre and near the wall (refer to Appendices A and B). This high core dissipation arises from the injection of large-scale, isotropic turbulent eddies that permeate the bulk flow, a phenomenon frequently observed in turbulent flows downstream of energetic mechanical sources (Durst, Jovanović & Sender Reference Elghobashi1995; Pope Reference Prince and Blanch2000). In this region, turbulence production is not yet dominated by wall-shear effects, and the dissipation profile lacks the classical near-wall peak that is characteristic of fully developed pipe flows (Eggels et al. Reference Estevadeordal and Goss1994). As the flow proceeds downstream, the imposed turbulence in the pipe core decays rapidly due to the absence of a sustaining production mechanism away from the wall (Pope Reference Prince and Blanch2000). In contrast, near-wall turbulence decays slower and is withheld by the strong mean velocity gradient at the solid boundary, which becomes the primary source of turbulent kinetic energy production and subsequent dissipation (Durst et al. Reference Elghobashi1995). Consequently, the radial dissipation profile evolves to exhibit a pronounced peak near the wall and a significant reduction at the centreline, reflecting the gradual transition from pump-driven, spatially uniform turbulence to a fully developed, wall-dominated turbulent regime (Eggels et al. Reference Estevadeordal and Goss1994). This observed evolution aligns with established theoretical and experimental findings on the development of turbulence in wall-bounded shear flows and underscores the interplay between initial isotropic turbulence and the emergence of wall-driven turbulent structures.

Figure 7 shows that the turbulent dissipation decays approximately as a power law ( $\varepsilon \propto \mathcal{L}^{-m}$ ) along the duct, appearing linear on a log–log plot with a best-fit slope $m$ narrowly confined to $m\simeq 1.85$ – $2.10$ . The spread is small enough that the corresponding error bars largely overlap. This behaviour holds over most axial stations, except near the end where the growing wall boundary layer sustains turbulence and the decay becomes noticeably slower. As we see weak radial variations in our measurements, this scaling can be compared with theory by assuming HIT. In decaying HIT, the $\varepsilon$ scaling follows a power law whose exponent is dependent on the large-scale invariant: where $k\ell ^{p}=$ constant with $p\in \{3,5\}$ (Batchelor & Proudman Reference Batchelor and Proudman1956; Saffman Reference Sajjadi, Raman, Shah and Ibrahim1967). This implies that

(4.1) \begin{align} k\ell ^{p}=\text{constant}\quad \Rightarrow \quad \ell \propto k^{-1/p}, \end{align}

where $\ell$ is the integral length scale. For HIT,

(4.2) \begin{align} \varepsilon \sim C_\varepsilon \,\frac {k^{3/2}}{\ell }\ \Rightarrow \ \varepsilon {\propto} k^{3/2+1/p},\qquad \frac {{\textrm d}k}{{\textrm d}t}=-\varepsilon \ \Rightarrow \ \frac {{\textrm d}k}{{\textrm d}t}{\propto} -\,{k^{\,{m}_{t}}}, \end{align}

where $t$ is time and $m_t= {3}/{2}+ {1}/{p}$ . Integration gives

(4.3) \begin{align} k(t) \propto \ t^{-\,\frac {1}{m_t-1}} &= t^{-\,\frac {2p}{p+2}}, \end{align}

(4.4) \begin{align} \ell (t)\ \propto \ k^{-1/p}\ \propto \ t^{\,\frac {2}{p+2}}, & \qquad \varepsilon (t)\ \propto \ \frac {k^{3/2}}{\ell }\ \propto \ t^{-\,\frac {3p+2}{p+2}}. \end{align}

Simplifying for the case $t\rightarrow \mathcal{L}/\mathcal{V}$ gives

(4.5) \begin{equation} \varepsilon (t)\propto t^{-m_t}\propto \mathcal{L}^{-m_t},\qquad m_t=\frac {3p+2}{p+2}. \end{equation}

The dimensionless time ( $t$ ) is defined as $\mathcal{L}/\mathcal{V}$ , where $\mathcal{V} = V/\bar {V}$ and $\bar {V}= {V_1 + V_2 + V_3}/{3}$ . Thus, for a Saffman (Reference Sajjadi, Raman, Shah and Ibrahim1967) $\kappa ^2$ spectrum ( $p=3$ ; permanence of big eddies for isotropic energy spectrum, $E(\kappa )\sim \kappa ^{2}$ as $\kappa \!\to \!0$ ; here $\kappa$ is the wavenumber), one gets $\varepsilon \sim \mathcal{L}^{-2.2}$ , while for a Loitsyanskii (Reference Lucas, Beyer and Szalinski1939) and Batchelor & Proudman (Reference Batchelor and Proudman1956) $\kappa ^4$ spectrum ( $p=5$ ), one gets $\varepsilon \sim \mathcal{L}^{-2.43}$ . These theoretical scalings are broadly consistent and closely match with experimental measurements (see figure 7). A small but systematic discrepancy remains: the experimentally observed decay is slightly slower than the ideal HIT prediction, which is primarily due to wall-bounded effects in the duct that sustain dissipation and retard the streamwise decay, more pronounced near the end of the duct.

4.2. Bubble size distribution and flow regime evolution

The modified Gaussian function, $f(d)$ , is widely used to compare and characterise the bubble size population (Friedlander et al. Reference Garrett and Li2000). Controlling the mean/mode, $\sigma$ , and $d_{\max}$ of the bubble size distribution in the flow is critical. Here $f(d)$ is represented by (4.6) and (4.7):

(4.6) \begin{equation} f(d) = \frac {\partial \varPhi }{\partial \ln d} = \frac {1}{\sqrt {2\pi } \ln \sigma _g} \exp \left [ -\frac {1}{2} \left ( \frac {\ln (d/d_g)}{\ln \sigma _g} \right )^2 \right ]\!, \end{equation}

where

(4.7) \begin{equation} d_g = \left ( \prod _{i=1}^{K} (d_i)^{n_i} \right )^{\frac {1}{K}} \,\,\,\,\, \text{and} \,\,\,\,\, \ln \sigma _g = \sqrt { \frac { \sum _{i=1}^{K} n_i (\ln d_i - \ln d_g)^2 }{N} }, \end{equation}

with $d_i$ the representative diameter for bin $i$ , $n_i$ the number of bubbles in that bin, $K$ the total number of bins, $d_g$ and $\sigma _g$ the geometric mean diameter and geometric standard deviation (Razzaque et al. Reference Risso and Fabre2003b ) and $\varPhi$ the cumulative fraction of bubbles with diameters smaller than $d$ . Quantifying $f(d)$ is technically essential because, as established by the statistical framework of Friedlander et al. (Reference Garrett and Li2000), multiplicative coalescence–breakup processes impose integral constraints on the geometric mean $d_g$ and variance $\sigma _g$ , thereby allowing the net interaction rates in population-balance models to be inferred from the measured distribution. Figure 8 presents the experimentally measured bubble size distributions $f(d)$ , constructed from the detected bubble diameters and shown on a semi-logarithmic scale. The curves represent log-normal fits evaluated using (4.6) and (4.7), with the parameters $d_g$ and $\sigma _g$ obtained directly from the measured data near the centre of the duct. Figure 8 ( $a\rightarrow b \rightarrow c$ , $d\rightarrow e\rightarrow f$ ) shows that increasing velocity ( $V$ ), or Re, shifts the $f(d)$ toward a smaller mean and standard deviation of bubble size distribution. At the highest $V$ (8.4 m s⁻¹), $f(d)$ exhibits a tall, sharp peak at a relatively small diameter, with a greatly diminished tail of large bubbles. In contrast, at the lowest $V$ (6.1 m s⁻¹) the distribution is broader and more right-sided, with a noticeable tail consisting of larger bubbles. Physically, the more energetic turbulence at higher $\varepsilon$ breaks bubbles into smaller sizes and rapidly fragments any large bubbles, yielding a narrower size spectrum (Prince & Blanch Reference Prince and Blanch1990). In fact, classic coalescence–breakup models of Prince & Blanch (Reference Prince and Blanch1990) treat turbulence as the driver of bubble collisions, but sufficient contact time is required for coalescence. Extremely intense turbulence shortens bubble contact times, suppressing coalescence efficiency, so breakup prevails and limits the upper bubble size. This is in agreement with the observed disappearance of the large-diameter tail at the highest velocity.

Figure 8.

Bubble size distribution with modified Gaussian function $f(d)$ near the centre at different $V$ and $\phi$ . Results are shown for (a) $V$ = 6.1 m s⁻¹ and $\phi$ = 0.5 %, (b) $V$ = 7.4 m s⁻¹ and $\phi$ = 0.5 %, (c) $V$ = 8.4 m s⁻¹ and $\phi$ = 0.5 %, (d) $V$ = 6.1 m s⁻¹ and $\phi$ = 2 %, (e) $V$ = 7.4 m s⁻¹ and $\phi$ = 2 %, (f) $V$ = 8.4 m s⁻¹ and $\phi$ = 2 %.

Moving downstream in the duct (from $\mathcal{L} \simeq 3.6$ to $40$ after injection) shown in figure 8(a,b,c,d,e,f) at given $\phi$ and $V$ , the bubble size distribution progressively broadens and develops a less pronounced peak. Near the injection point ( $\mathcal{L} \simeq 3.6$ downstream), $f(d)$ is relatively narrow with a small tail of both very small and a few large bubbles – a signature of the initial breakup of the air flow by the pump’s turbulence. With an increase in $\mathcal{L}$ (shown for $\mathcal{L} = 40$ in figure 8), the distribution becomes broader around a modal diameter and the extreme (specifically very large bubbles) are more prevalent. This broadening of the distribution reflects the combined effects of decaying turbulence and ongoing bubble–bubble interactions. As the turbulence energy dissipates downstream, fewer new microbubbles are generated by breakup, and the small bubbles present begin to coalesce into mid-sized ones. The net result is an evolving equilibrium toward coalescence: the overall bubble count decreases with distance, indicating that many small bubbles are merging into larger ones. Notably, the large bubble tail of the distribution grows significantly by $\mathcal{L}\sim$ 40, suggesting that large bubbles are formed via coalescence in the still turbulent (though weakening) flow. Thus, with downstream distance the bubble population shifts toward larger characteristic diameters, as an initially narrow distribution broadens and the peak amplitude decreases. Similar trends have been noted in spatially developing, buoyancy-driven bubbly flows (e.g. Ruiz-Rus et al. Reference Saffman2022), but the present results arise in a spatially decaying, very-high-Reynolds-number turbulent flow, where coalescence likewise produces distributions biased toward larger bubbles, consistent with observations in fully developed horizontal pipe flows (Razzaque et al. Reference Risso and Fabre2003b ) and in decaying turbulence (Serizawa et al. Reference Shi1975). Our findings here similarly indicate that in the axial direction, the turbulence-decaying, coalescence-dominated regime yields a lower peak and broader $f(d)$ distribution.

Higher $\phi$ leads to a broader and flatter bubble size distribution as shown in figure 8( $a\rightarrow d$ , $b\rightarrow e$ , $c\rightarrow f$ ), with an enhanced probability of larger bubbles. In these higher void fraction cases, the distribution’s large-diameter tail becomes more pronounced, meaning a greater presence of big bubbles compared with lower void conditions. The physics driving this trend is the increased frequency of bubble–bubble collisions at higher $\phi$ . Even though the flow is turbulent, a higher collision rate allows many small bubbles to merge into larger ones and decaying turbulence provides stability. Although rapidly decaying, developing turbulent multiphase flows are scarcely examined, our observations align with results for fully developed, coalescence-prone turbulence: increasing the void fraction $\phi$ raises the mean bubble size and broadens the distribution (Lehr et al. Reference Li, Liu, Xu and Li2002). At high $\phi$ , collision probabilities, and hence, collision/coalescence frequencies, increase (Prince & Blanch Reference Prince and Blanch1990; Razzaque et al. Reference Risso and Fabre2003b ), directly shaping the size–distribution function $f(d)$ . By contrast, in strongly breakup-dominated regimes, bubble size can remain nearly independent of $\phi$ up to a threshold. In our measurements, coalescence plays a measurable role: as $\phi$ increases, enhanced bubble–bubble interactions shift $f(d)$ toward larger diameters and yield a broader, right-sided distribution, consistent with Lance & Bataille (Reference Lee1991).

Figure 9.

Scaling of compensated probability density function (PDF $* (d/d_{\textit{H}})^{({3}/{2})}$ ). The bubble size ( $d$ ) is normalised with the Hinze scale ( $d_{\textit{H}}$ ). Sub-Hinze power-law scaling: $d^{-3/2}$ ; Super-Hinze power-law scaling: $d^{-10/3}$ . (a) PDF near centre ( $\mathcal{L}$ = 3.6) at $\phi$ = 0.5 %, (b) PDF near centre ( $\mathcal{L}$ = 40) at $\phi$ = 0.5 %, (c) PDF near centre ( $\mathcal{L}$ = 3.6) at $\phi$ = 2 %, (d) PDF near centre ( $\mathcal{L}$ = 40) at $\phi$ = 2 %.

While mean diameters and distribution widths characterise bulk bubble evolution, the probability density function (PDF) of sizes captures the complete statistical structure of the population. The PDF provides physical characterisation of the distribution: the presence and slope of power-law tails and their axial evolution – identify the dominant mechanism (inertial/capillary breakup versus coalescence) and delineate regime transitions (Deane & Stokes Reference Diaz, Gonzalez and Vera2002; Chan et al. Reference Chen, Jhuang, Wu, Yang, Wang and Chen2018). Figure 9 compares the bubble size PDF measured near the duct centreline at two axial locations ( $\mathcal{L}$ = 3.6 and 40) for two void fractions ( $\phi = 0.5\,\%$ and $2\,\%$ ) and varying bulk velocities (or $Q$ ). Near the inlet ( $\mathcal{L} = 3.6$ ), both void fraction cases (figures 9 a and 9 c) exhibit a similar dual-slope PDF behaviour on log–log plots. Bubble diameters in sub-Hinze follow an approximate $d^{-3/2}$ power law, indicating a no active breakup regime (mostly coalescence) (Crialesi-Esposito, Chibbaro & Brandt Reference Dabiri, Lu and Tryggvason2023). In this region ( $\mathcal{L}$ = 3.6 for all $\phi$ and $V$ ) virtually all of the bubble population lies under the Hinze diameter ( $d/d_{\textit{H}} \lt 1$ ), so turbulent eddies are too weak to overcome surface tension and break the bubbles. Consistently, breakup is relatively small at $\mathcal{L} = 3.6$ but bubble growth is governed by both coalescence and breakup. The observed $-3/2$ scaling is consistent with prior findings by Deane & Stokes (Reference Diaz, Gonzalez and Vera2002) in flows lacking active fragmentation, where this exponent is classically associated with sub-Hinze bubble generation via capillary mechanisms observed in breaking waves or turbulent-induced pinch-off events (Deane & Stokes Reference Diaz, Gonzalez and Vera2002). For example, oceanic wave-breaking experiments showed that sub-Hinze bubbles obey a $d^{-3/2}$ size spectrum, similar to our inlet measurements (Deane & Stokes Reference Diaz, Gonzalez and Vera2002; Farsoiya et al. Reference Garcia-Ochoa and Gomez2023; Roccon, Zonta & Soldati Reference Ruiz-Rus, Ern, Roig and Martínez-Bazán2023; Jain & Elnahhas Reference Javadi, Kumar and Aidun2025). This finding also accords with the notion that as long as bubbles remain smaller than the critical Hinze size for breakup, successive collisions will merge them into larger bubbles and broaden the distribution without any competing fragmenting mechanism (Guo et al. Reference Guo, Zhou, Li and Chen2016). Our data indeed provide direct evidence of such a coalescence-driven cascade: an initially narrow bubble size distribution at the duct inlet rapidly evolves into a broader $d^{-3/2}$ distribution before any significant breakup occurs. Notably, this near inlet $d^{-3/2}$ PDF holds for all tested flow conditions – all three liquid velocities and both void fractions collapse to the same slope at $\mathcal{L} = 3.6$ – underscoring that the early developing flow is a robust coalescence-dominated regime largely independent of $\phi$ or $\mathcal{I}$ (Prince & Blanch Reference Prince and Blanch1990; Farsoiya et al. Reference Garcia-Ochoa and Gomez2023). Although most of the bubbles are sub-Hinze scale, there are considerable amounts of bubbles with super-Hinze scale ( $d \gt d_{\textit{H}}$ ) at $\mathcal{L} = 3.6$ due to still high $\varepsilon$ . These super-Hinze scale bubbles are susceptible to turbulent breakup, causing the PDF’s upper tail to steepen toward a ${\sim} d^{-10/3}$ slope. Figure 9 shows that once bubbles exceed the local $d_{\textit{H}}$ , they indeed begin to fragment and the distribution tail assumes the characteristic ${\sim} d^{-10/3}$ form associated with inertial breakup. This $d^{-10/3}$ scaling for fragmenting bubbles is consistent with classical fragmentation cascade models (Garrett & Li Reference George and Hussein2000) and has been observed in experiments and simulations of turbulent dispersions (Soligo, Roccon & Soldati Reference Tan, Zhong, Qi, Xu and Ni2019). Near the tail of figure 9(a,c), the distributions for $d/d_{\textit{H}}\gt 1$ at $\mathcal{L}=3.6$ exhibit slopes steeper than $-10/3$ , ranging from $-3.8$ to $-4.5$ , i.e. close to but slightly steeper than the inertial-breakup scaling. Owing to the low void fraction, the occurrence of large bubbles is sparse, which can make the tail appear artificially steep, consistent with the trend reported by Farsoiya et al. (Reference Garcia-Ochoa and Gomez2023).

Far downstream ( $\mathcal{L} = 40$ ), the bubble size distribution has shifted and developed a pronounced single-slope scaling. All bubbles remain sub-Hinze and continue to follow the coalescence-dominated $d^{-3/2}$ trend. Unlike an equilibrium regime where breakup balances coalescence effects, here coalescence still dominates the majority of the duct region, while limited breakup occurs only near the inlet ( $\mathcal{L} \lt 5{-}8$ , depending on $V$ ) where $\varepsilon$ is very high. In fact, beyond the $\mathcal{L}\gt 5{-}8$ downstream of the inlet, turbulent energy has decayed so substantially that super-Hinze scale bubble breakup becomes non-existent – the flow effectively transitions into a nearly purely coalescence-driven regime. The disappearance of the super-Hinze tail with increasing $\mathcal{L}$ (from $\mathcal{L}=3.6$ ) indicates that while some bubbles briefly exceed the Hinze scale, they are rapidly broken up, preventing a persistent super-Hinze population. This behaviour is fully consistent with the concept of the Hinze scale acting as a cutoff: a critical diameter $d_{\textit{H}}$ exists at which turbulent dynamic pressure balances capillary pressure ( $\rho \mathcal{V}^2 \sim \gamma /d_{\textit{H}}$ ), corresponding to a Weber number of order unity ( $We \sim \mathcal{O}(1)$ ) for breakup onset (Hinze Reference Hosokawa and Tomiyama1955). This equilibrium size concept has long been utilised in population-balance models to predict stable bubble size distributions in steady-state bubbly flows (Lehr et al. Reference Li, Liu, Xu and Li2002). In our decaying flow, however, the continuous reduction in turbulence intensity with increasing $\mathcal{L}$ ensures that bubbles rarely exceed the Hinze scale. Bubbles that do exceed the Hinze scale are not actively fragmented. Consequently, beyond a few hydraulic diameters from the inlet, bubble coalescence mostly proceeds unopposed by breakup.

The influences of $\phi$ and $V$ on the PDF scaling trends are also evident in figure 9. Near the inlet ( $\mathcal{L} \leq 8.2$ ), cases with a higher void fraction ( $\phi = 2\,\%$ vs $0.5\,\%$ ) exhibit a broader distribution due to accelerated coalescence-driven bubble growth, along with a more pronounced tail of large bubbles following the $d^{-10/3}$ scaling (comparing figures 9 c and 9 a). Although at $\mathcal{L} = 40$ , the $2\,\%$ case exhibits predominantly a single-slope regime, it still displays a longer large bubble tail compared with the $0.5\,\%$ case – again attributable to the enhanced coalescence at higher $\phi$ (comparing figures 9 d and 9 b). This trend aligns with theoretical expectations: a higher concentration of bubbles increases the frequency of collisions and consequently the coalescence rate, allowing bubbles to reach the critical Hinze scale more rapidly (Guo et al. Reference Guo, Zhou, Li and Chen2016). In practical terms, increasing the void fraction enhances the extent of the distribution and leads to an earlier formation of large, unstable bubbles susceptible to breakup within a few hydraulic diameters from the inlet – manifesting as an extended $d^{-10/3}$ fragmentation tail. Downstream of this region, even as the turbulence continues to decay, the higher $\phi$ case maintains a larger population of bubbles and a longer $d^{-3/2}$ tail relative to the lower $\phi$ case.

Bulk liquid velocity exerts a subtler yet important influence on the PDF scaling. Remarkably, for a given void fraction $\phi$ , the normalised distributions, plotted as PDF vs $d/d_{\textit{H}}$ , collapse onto one another, showing strong overlap across different axial locations $\mathcal{L}$ and velocities $V$ . This collapse suggests a self-similar evolution of the bubble size distribution when scaled by the local Hinze diameter. The behaviour can be understood by examining the role of $\varepsilon$ : higher liquid velocities inject more turbulent energy, leading to greater dissipation rates ( $\varepsilon$ ) and consequently smaller Hinze scales ( $d_{\textit{H}}$ ). In these high- $V$ cases, strong turbulence suppresses the coalescence-driven growth (due to reduced coalescence time scales) of large bubbles by fragmenting them early, effectively ‘nipping in the bud’ any transition toward a broader distribution (Prince & Blanch Reference Prince and Blanch1990). As a result, the bubble population in high- $V$ flows remains narrowly distributed and biased toward smaller diameters, with minimal spread around the mean. Conversely, lower velocity flows (e.g. $V = 6.1$ m s⁻¹) exhibit weaker turbulence (smaller $\varepsilon$ ) and, thus, larger Hinze diameters (Yang et al. Reference Yang, Sun, Yang, Xing and Gui2025). Under such conditions, bubbles can grow to larger sizes before turbulent stresses are sufficient to induce breakup, resulting in a broader size distribution and a more prominent large bubble population. Although the low- $V$ flows still experience decaying turbulence and correspondingly smaller $d_{\textit{H}}$ farther downstream, the diminished fragmentation enables broader coalescence-driven growth. These results confirm that the bubble size distribution in decaying turbulent flows scales robustly with the local Hinze diameter, and the PDFs exhibit a self-similar trend across a range of flow conditions.

Figure 10.

Cumulative bubble size distribution ( $\varPhi$ ) with bubble size ( $d/d_{32}$ ) at $V=$ 6.1 m s⁻¹. The slopes $\theta (\mathcal{L}_1,\mathcal{L}_2,\mathcal{L}_3)$ are reported at three axial locations, where $\mathcal{L}_1 = 3.6$ , $\mathcal{L}_2 = 21.8$ and $\mathcal{L}_3 = 40$ . Here $\sigma _g$ is calculated using (4.7). Results are shown for (a) $\phi$ = 0.5 %, $\theta \in$ (3.63,2.11,1.59), and $\sigma _g \in$ (1.48,1.76,1.91); (b) $\phi$ = 1 %, $\theta \in$ (2.61,1.76,1.55), and $\sigma _g \in$ (1.54,1.84,1.96); (c) $\phi$ = 2 %, $\theta \in$ (2.52,1.45,1.39), and $\sigma _g \in$ (1.57,1.89,2.03).

The cumulative bubble size distribution $\varPhi$ , plotted as a function of the normalised diameter $d/d_{32}$ , provides a compact description of the evolving bubble population in the decaying turbulent flow (Razzaque et al. Reference Risso and Fabre2003b ). Figure 10 shows $\varPhi$ at $V=6.1\,\text{m}\,{\textrm s}^{-1}$ for $\phi =0.5\,\%,\,1\,\%,\,2\,\%$ , comparing measurements at $\mathcal{L}=3.6,\,21.8$ and $40$ . For each case, we report a fitted slope $(\theta )$ , together with the $\sigma _g$ , which directly characterises the width of the corresponding log-normal distribution. At $\phi =0.5\,\%$ , the distribution near the inlet ( $\mathcal{L}=3.6$ ) is relatively steep ( $\theta =3.63$ ), whereas downstream the curves flatten ( $\theta =2.11$ at $\mathcal{L}=21.8$ and $\theta =1.59$ at $\mathcal{L}=40$ ). Consistently, $\sigma _g$ increases from 1.48 to 1.91 with $\mathcal{L}$ , indicating progressive broadening of the distribution as the flow convects downstream. These trends are consistent with earlier pipe flow studies (Hesketh et al. Reference Hibiki and Ishii1987; Razzaque et al. Reference Risso and Fabre2003b ). Hesketh et al. (Reference Hibiki and Ishii1987) showed that, in breakup-dominated conditions, the distribution approaches a self-preserving log-normal form with nearly invariant width ( $\sigma _g \approx 1.2$ – $1.4$ ). In contrast, Razzaque et al. (Reference Risso and Fabre2003b ) demonstrated that, under coalescence-dominated conditions, the distribution broadens downstream, with $\sigma _g$ typically increasing to ${\approx} 1.7$ – $2.0$ . Thus, the distribution width is not universal, but reflects the local competition between breakup and coalescence. For $\phi = 0.5\,\%$ at $\mathcal{L}=3.6$ , $\sigma _g \approx 1.48$ , reflecting a slight widening beyond the canonical breakup range. This is consistent with the fact that the strongest turbulent fragmentation occurs within the pump, whereas measurements are taken in the duct where turbulence has already decayed. The inlet region therefore retains significant breakup, but is no longer fully breakup dominated. As turbulence decays further downstream, breakup weakens, coalescence becomes increasingly effective and the distributions migrate into the coalescence-dominated range ( $\sigma _g \gtrsim 1.7$ ), in agreement with observations that turbulence decay promotes bubble growth through coalescence (Chan et al. Reference Chen, Jhuang, Wu, Yang, Wang and Chen2018; Chieco & Durian Reference Coulaloglou and Tavlarides2023). Increasing $\phi$ produces flatter cumulative curves (smaller $\theta$ ) and larger $\sigma _g$ . The stronger broadening at higher $\phi$ reflects enhanced collision frequency and coalescence probability. With increasing $\phi$ ( $0.5\,\% \rightarrow 1\,\% \rightarrow 2\,\%$ ), $\sigma _g$ approaches $1.7$ – $2.0$ faster, characteristic of strongly coalescence-dominated conditions, whereas lower void fractions require longer axial development to reach similar states.

Figure 11.

Axial variation of $d_{99.8}, d_{32}$ and $d_{\textit{H}}$ in the developing turbulent duct flow. The reported fit parameters $\beta _{\textit{exp}}$ (power law coefficient) and $R^2$ (coefficient of determination) are obtained from data for $\mathcal{L}\gt 3.6$ , excluding the first data point, and are reported in table 1. Results are shown for (a) $V$ = 6.1 m s⁻¹, (b) $V$ = 7.4 m s⁻¹, (c) $V$ = 8.4 m s⁻¹, (d) the legend.

The axial evolution of the bubble size distribution, combining $d_{32}$ with high quantiles such as $d_{99.8}$ (99.8th percentile bubble diameter) and referencing the Hinze scale $d_{\textit{H}}$ , captures both bulk and tail dynamics as turbulence decays. Identifying shifts between breakup-limited and coalescence-dominated regimes is critical, since the former suppresses large bubbles while the latter promotes tail broadening and enhanced interfacial renewal. Regime-resolved distribution tracking therefore provides stringent calibration targets for population-balance models and a mechanistic basis for reliable scale-up.

Figure 11 presents the axial variation of bubble size statistics near the centre of the duct ( $6 \times 6\,\text{mm}^2$ ) on a log–log scale, namely $d_{32}$ , $d_{99.8}$ , and the Hinze scale $d_{\textit{H}}$ , for each void fraction $\phi$ and bulk velocity $V$ . The Hinze scale is evaluated locally at each axial position using (1.1). As $\mathcal{L}$ increases, the Hinze scale $d_{\textit{H}}$ grows linearly on a log–log scale (following the power law). Downstream, the dissipation $\varepsilon$ decreases markedly (by approximately $90\,\%$ over the test section; see figure 7). Through (1.1), this gives $d_{\textit{H}} \propto \varepsilon ^{-2/5}$ , so a rapid fall in $\varepsilon$ produces only a gradual increase in $d_{\textit{H}}$ . If $\varepsilon \sim \mathcal{L}^{-m}$ with $m \approx 1.85$ to $2.2$ (see figure 7), then $d_{\textit{H}} \sim \mathcal{L}^{2m/5}$ , that is, $d_{\textit{H}} \sim \mathcal{L}^{0.74}$ to $\mathcal{L}^{0.88}$ . Thus, the growth of $d_{\textit{H}}$ with $\mathcal{L}$ becomes sublinear on the log–log scale.

Table 1.

Power-law coefficient ( $\beta _{\textit{exp}}$ ) and coefficient of determination ( $R^2$ ) for the fits for the axial variation of $d_{99.8}, d_{32}$ and $d_{\textit{H}}$ shown in figure 11.

In comparison, $d_{99.8}$ also grows linearly (on log–log; general power law) with a smaller slope along the axial direction and remains below $d_{\textit{H}}$ at all axial locations except near the inlet ( $\mathcal{L}\approx 3.6$ ). Near the inlet, the size distribution straddles the Hinze scale, so bubbles exist on both sides of $d_{\textit{H}}$ ; those with $d\gt d_{\textit{H}}$ are susceptible to turbulence-induced breakup. Downstream, as turbulence decays, $d_{\textit{H}}$ increases faster than $d_{99.8}$ , so $d_{99.8}$ falls below and stays below $d_{\textit{H}}$ . After $\mathcal{L}\gt 8.2$ , roughly no bubbles exceed $d_{\textit{H}}$ , which indicates a coalescence–dominated regime with active breakup effectively suppressed. The $d_{32}$ also shows a power-law growth (linear on log–log) with $\mathcal{L}$ (except near the inlet), but it remains well below both $d_{99.8}$ and $d_{\textit{H}}$ . The monotonic rise of $d_{32}$ is consistent with coalescence-driven growth.

For $\mathcal{L}\ge 8.2$ , both $d_{99.8}$ and $d_{32}$ follow power-law trends with $\mathcal{L}$ in figure 11, and the corresponding fit parameters are listed in table 1. The fits are strong, with $R^2$ values close to $0.95$ . The fitted exponents $\beta _{\textit{exp}}$ lie within $0.45$ to $0.51$ , indicating that the mean and the upper-tail bubble sizes increase at comparable rates in the coalescence-dominated regime. The theoretical scaling for the growth of $d_{32}$ and $d_{99.8}$ follows from the effective coalescence rate $\varGamma$ , defined as the product of the collision frequency $h$ and the coalescence efficiency $\lambda$ (Coulaloglou & Tavlarides Reference Nambiar, Kumar and Gandhi1977; Prince & Blanch Reference Prince and Blanch1990):

(4.8) \begin{align} \varGamma = h\,\lambda . \end{align}

Classical inertial-range scaling (Kolmogorov Reference Kulkarni1991), as presented in standard treatments (Batchelor Reference Batchelor1953; Monin & Yaglom Reference Monin and Yaglom2013), underpins collision models for bubbly dispersions. In these measurements, bubble interactions are driven by inertial-subrange eddies, consistent with the observed length-scale ordering $\eta \ll d \sim \lambda _T \ll \ell$ , for which $d/\eta \gg 1$ and $d/\lambda _T = \mathcal{O}(1)$ . In the inertial subrange of turbulence, the velocity increment $u_{\textit{rel}}$ across a bubble scale $d$ follows Kolmogorov similarity, resulting in (Batchelor Reference Batchelor1953; Monin & Yaglom Reference Monin and Yaglom2013)

(4.9) \begin{align} u_{\textit{rel}}(d)\sim (\varepsilon d)^{1/3}, \end{align}

where $\varepsilon$ is the mean turbulent dissipation rate. The relative-velocity scaling is based on interactions between bubbles and inertial-range eddies of comparable size. Accordingly, the collision velocity is governed by the turbulent velocity difference evaluated at the inter-bubble separation scale (Coulaloglou & Tavlarides Reference Nambiar, Kumar and Gandhi1977; Prince & Blanch Reference Prince and Blanch1990; Luo & Svendsen Reference Matiazzo, Decker, Bastos, Silva and Meier1996; Kamp et al. Reference Khodaparast, Borhani, Tagliabue and Thome2001). Combining the inertial–range relative velocity with the geometric cross-section ${\sim} d^{2}$ yields the turbulent collision kernel

(4.10) \begin{align} h \propto u_{\textit{rel}}(d) \times d^2 \propto \varepsilon ^{1/3}\, d^{7/3}. \end{align}

As shown in Appendix C, the assumption of spherical bubbles is justified because the bubble interfacial relaxation time is comparable to or shorter than the integral-scale eddy turnover time, allowing the interface to re-equilibrate between successive collisions (Kumar et al. Reference Laakkonen, Alopaeus and Aittamaa2026). Further simplification results in

(4.11) \begin{equation} \varGamma = 8C\,\lambda \,(\varepsilon )^{1/3}\,d^{\frac {7}{3}} . \end{equation}

The kernel has units of volume per unit time and the coalescence efficiency satisfies $\lambda \in [0,1]$ . The evolution of the number density $n$ can then be expressed using Smoluchowski’s binary coagulation equation (accounting for double counting and assuming a monodisperse distribution) (Kreer & Penrose Reference Kumar, Jain, Upadhyay, Bharath and Waghmare1994):

(4.12) \begin{equation} \frac {\text{d}n}{\text{d}t} = -\frac 12\,\varGamma \,n^2. \end{equation}

The rationale for this monodisperse simplification is that, in the present high-turbulence regime, both collisions and successful coalescence are dominated by similar-sized bubbles, with the coalescence efficiency highest when bubble sizes are comparable, leading to a kernel that is sharply peaked near unity size ratios (Liao & Lucas Reference Loitsyanskii2010; Tan et al. Reference Tavoularis and Corrsin2025). For duct flow, substituting $t=\mathcal{L}/\mathcal{V}$ and $n=6\phi /(\pi d^{3})$ into (4.12), and simplifying for decaying turbulence with $\varepsilon =\varepsilon (\mathcal{L})$ (refer Appendix C), we obtain

(4.13) \begin{equation} \frac {\text{d}d}{\text{d}\mathcal{L}} \equiv K_o [\varepsilon (\mathcal{L})]^{1/3}d^{\frac {1}{3}}, \end{equation}

where $K_o$ depends on $\lambda$ . Under the present flow conditions, $\lambda$ is evaluated using the commonly employed film drainage model and is found to be a constant (see Appendix C). Using a power-law decay referenced to a finite start location $\mathcal{L}\gt 0$ (see (4.5)):

(4.14) \begin{equation} d(\mathcal{L}) \propto \mathcal{L}^{\frac {3-m}{2}} \, \propto \, \mathcal{L}^{\beta _{\textit{th}}}, \quad \text{where} \quad \beta _{th} = \frac {3-m}{2}. \end{equation}

The theoretical power-law coefficient $\beta _{\textit{th}}$ in (4.14) is listed in table 2. For $m$ in the range $1.85$ – $2.10$ , $\beta _{\textit{th}}$ varies from $0.45$ to $0.57$ , which closely matches the experimentally fitted exponents for $d_{32}$ and $d_{99.8}$ , from $0.44$ to $0.53$ , in table 1. Thus, a steeper dissipation decay (larger $m$ ) corresponds to a larger $\beta _{\textit{th}}$ . At fixed $m$ , higher sustained $\varepsilon$ increases the collision frequency and accelerates bubble size growth. For applications that require tight control of the size distribution (for example, targeting $d_{32}$ or $d_{99.8}$ ), one may adjust $m$ by passive means such as modifying wall conditions or by active means such as acoustic or ultrasonic forcing. To make $d$ grow in parallel with the Hinze scale, one could equate the scalings $d \sim \mathcal{L}^{(3-m)/2}$ and $d_{\textit{H}} \sim \mathcal{L}^{2m/5}$ , which gives $m=5/3$ . At $m=5/3$ , both $d_{32}$ and $d_{99.8}$ remain self-similar and grow in parallel with $d_{\textit{H}}$ as $\mathcal{L}$ increases.

Table 2.

Bubble diameter scaling in highly decaying turbulent flow. Here $ \beta _{th}= {3-m}/{2}$ .

Near the inlet, the pure-coalescence scaling does not apply because breakup remains active. Figure 11 shows that, owing to breakup and an initially unstable distribution, $d_{99.8}$ is capped and grows at a different rate than $d_{32}$ ; consequently, the ratio $\mathcal{D}=d_{99.8}/d_{32}$ evolves through the entry region before saturating downstream (once the two growth curves become roughly parallel). Consistent with this behaviour, Razzaque et al. (Reference Risso and Fabre2003b ) reported that monodisperse injections ( $\mathcal{D}\approx 1$ ) in fully developed, coalescence-dominated turbulence relax to a constant $\mathcal{D}\approx 2.2$ . Monitoring $\mathcal{D}$ is operationally important in mixing, heat-transfer and reactor applications because it governs interfacial area and distribution shape (Lehr et al. Reference Li, Liu, Xu and Li2002; Risso Reference Roccon, Zonta and Soldati2018).

The axial evolution of $\mathcal{D}$ quantifies tail heaviness relative to the interfacial-area–weighted mean, providing a robust, outlier-resistant proxy for the largest statistically reliable bubbles (preferable to raw $d_{\text{max}}$ ) and for population broadening as turbulence decays (Razzaque et al. Reference Risso2003a ; Risso & Fabre Reference Rong, Li, Zhang and Sun1998). Figure 12 presents the streamwise evolution of the normalised bubble size ratio $\mathcal{D}^*$ (normalised by $\mathcal{D}_o$ , which is $\mathcal{D}\, \text{at}\, \mathcal{L}=0$ ) for three $\phi$ at three $V$ . The initial $\mathcal{D}_o$ ( ${\approx} 1.91\pm 0.02$ ) is nearly invariant across $\phi$ and $V$ because intense inlet turbulence rapidly breaks larger bubbles and promotes the swift coalescence of smaller ones, yielding a narrow distribution. The figure shows $\mathcal{D}^*$ is virtually the same at the duct centre and near the wall (differences $\lt$ 5 %), indicating that the bubble size distribution width is almost uniform across the cross-section. This is consistent with the expectation that in a fully developed (dynamic equilibrium) bubbly flow, the size distribution becomes independent of position (Kleinstreuer & Agarwal Reference Kolmogorov2001). Across all the variables, $\mathcal{D}^*$ stays within a range of approximately 1–1.15. In the developing region (near the inlet), $\mathcal{D}^*$ rises with axial distance $\mathcal{L}$ , reflecting the broadening of the bubble size distribution due to coalescence. Higher void fractions (1 % and 2 %) reach this plateau earlier (at smaller t), whereas at 0.5 % the ratio requires a longer time to level off. This trend indicates a coalescence-dominated flow development – as bubbles travel downstream, frequent coalescence creates larger bubbles and a wider size spread until a balance is achieved. The eventual ‘saturation’ value of $\mathcal{D}^*$ falls around 1.15 (where $\mathcal{D} \approx$ 2.2) for all velocity and void fraction combinations (with the caveat that the 0.5 % case may not have fully reached its plateau in the available length). Once this value is attained, the distribution width no longer grows with $\mathcal{L}$ , implying that breakup and coalescence rates have equilibrated.

Figure 12.

Variation of $\mathcal{D}^*$ (normalised) in a developing turbulent duct flow. Results are shown for (a) $\phi$ = 0.5 %, (b) $\phi$ = 1 %, (c) $\phi$ = 2 %.

The observed plateau ratio of $\mathcal{D}^*\approx$ 1.15 ( $\mathcal{D}\approx$ 2.2) is consistent with the results of Razzaque et al. (Reference Risso and Fabre2003b ) for fully developed turbulent flow fields dominated by coalescence. The value $\mathcal{D}_o \approx 1.9$ is characteristic of a breakup-controlled regime: Evans, Jameson & Atkinson (Reference Farsoiya, Liu, Daiss, Fox and Deike1992) reported $\mathcal{D} \approx 1.8$ in liquid-jet bubble columns where turbulent breakup primarily set the bubble sizes. The initial slightly larger value reflects the short delay between bubble formation (due to breakup) in the pump and the first viewing station, during which limited decay and early coalescence shift the upper tail upward. Thus, the inlet condition sampled here corresponds to a near-breakup-controlled distribution, from which the flow subsequently evolves into a coalescence-dominated regime downstream. Notably, Razzaque et al. (Reference Risso and Fabre2003b ) also found that in the coalescence-dominated regime this ratio becomes nearly independent of axial position, bulk velocity and void fraction, mirroring our observation that once the distribution equilibrates, further axial evolution of $\mathcal{D}^{*}$ is weak. Prince & Blanch (Reference Prince and Blanch1990)’s classic bubble-column experiments, together with their population-balance model, showed that beyond a certain column height (the ‘dynamic equilibrium region’) the bubble size distribution becomes independent of axial position and reflects a self-preserving distribution, and Lehr et al. (Reference Li, Liu, Xu and Li2002) consistently demonstrated in simulations that once coalescence and breakup reach balance, the distribution shape no longer evolves with height. In summary, the plateau in $\mathcal{D}^{*}$ marks the attainment of a self-preserving bubble size distribution: the overall shape no longer evolves with $\mathcal{L}$ , even though individual bubbles continue to merge and break. The common asymptotic value of $\mathcal{D}\approx 2.2$ reported in the literature therefore serves as a practical reference for identifying equilibrium conditions and for evaluating breakup–coalescence models.

The exponential broadening of the size ratio arises from several coalescence mechanisms acting together. Random turbulent fluctuations in the liquid velocity field induce stochastic bubble–bubble encounters and bring bubbles from different streamlines into contact, increasing the collision frequency beyond laminar drift (Prince & Blanch Reference Prince and Blanch1990). In our high-Re flow, the random bubble–bubble encounters are overwhelmingly driven by the strong background turbulence; estimates of the bubblance parameter ( $b \approx 10^{-5}$ – $10^{-3}$ , ratio of bubble-induced turbulence to background turbulence) and the associated time scales confirm that the bubble-induced agitation is negligible by comparison. Turbulence-induced coalescence has therefore been modelled as a dominant mechanism in similar systems (Kamp et al. Reference Khodaparast, Borhani, Tagliabue and Thome2001; Bröder & Sommerfeld Reference Bunner and Tryggvason2004). A second contribution arises from velocity shear and differential buoyancy. Faster-rising bubbles overtake slower ones, and velocity gradients near the wall promote collisions, consistent with Prince and Blanch’s identification of turbulence, buoyancy-induced velocity differences, and laminar shear as primary collision sources (Prince & Blanch Reference Prince and Blanch1990). All three are relevant here: even in co-current flow, larger bubbles have higher rise velocities relative to the liquid and can catch up to smaller ones (a buoyancy-driven overtaking collision), while shear in the liquid especially near walls creates lateral motion that brings bubbles into contact. Third, radial migration redistributes bubbles in a way that further enhances coalescence. Small nearly spherical bubbles tend to move toward the wall, whereas larger deformable bubbles migrate toward the centreline (Lucas, Beyer & Szalinski Reference Luo and Svendsen2023). As bubbles grow, they accumulate in the core where collision frequency increases, while smaller bubbles remain near the periphery with fewer encounters. Wake entrainment behind large bubbles further increases collision likelihood (Ruiz-Rus et al. Reference Saffman2022). The resulting radial structure reflects turbulent dispersion, shear-induced collisions, lateral lift-induced migration, wake entrainment and differential buoyancy effects (negligible here) (Tomiyama Reference Valente and Vassilicos1998; Tomiyama et al. Reference van den, Thomas, Luther, Lathrop and Lohse2002). In flows containing small bubbles and very high carrier-flow turbulence, bubble–bubble collisions are dominated by the carrier turbulence by more than an order of magnitude, whereas bubble-induced turbulence, shear, buoyancy-driven motion, lift-induced migration and wake capture primarily act to redistribute bubbles rather than generate additional collisions (Lance & Bataille Reference Lee1991; Prakash et al. Reference Prakash, Martinez Mercado, van Wijngaarden, Mancilla, Tagawa, Lohse and Sun2016; Alméras et al. Reference Alméras, Mathai, Lohse and Sun2017; van den Berg et al. Reference Bouche, Roig, Risso and Billet2020). Together, these mechanisms produce an initial burst of coalescence when many small bubbles are present, followed by slower coalescence as the bubble number decreases and the distribution approaches equilibrium.

Figure 13.

Transient variation of ( $\mathcal{D}$ ) $^*$ (curve fitted) in a developing turbulent duct flow. Results are shown for (a) $\phi$ = 0.5 %, (b) $\phi$ = 1 %, (c) $\phi$ = 2 %.

As discussed, the $\mathcal{D}^*$ increases monotonically with $\mathcal{L}$ . Empirically, the evolution follows an exponential relaxation towards an asymptote as seen in figure 13 and fitted with

(4.15) \begin{align} \mathcal{D}^* &= 1.15 - 0.15 {\textrm e}^{- \frac {t}{\tau }}; \quad \text{where} \quad \mathcal{D}^* = \frac { \mathcal{D}}{ \mathcal{D}_o} , \end{align}

(4.16) \begin{align} \Delta \mathcal{D}^* (\,\%) &= \frac { \mathcal{D}_{\text{max}}^* - \mathcal{D}^*} {\mathcal{D}_{\text{max}}^*} \times 100 \approx 13\,{\textrm e}^{- \frac {t}{\tau }} , \end{align}

where $\tau$ is a characteristic relaxation time. Near the duct inlet at $\mathcal{L}=0$ , $\mathcal{D}^*=1$ by definition; downstream it rises toward an asymptotic value of 1.15 (a 15 % increase). This exponential trend indicates a first-order process: the rate of change of the size ratio is highest near the inlet (where the distribution is narrow) and decays progressively as the distribution broadens and approaches a new equilibrium state. Such behaviour is consistent with classical coagulation kinetics, wherein the population’s evolution can be characterised by a single dominant time scale (Ruiz-Rus et al. Reference Saffman2022). The form of the above relation is analogous to Smoluchowski’s theory for Brownian coalescence (with $\tau$ playing the role of a ‘half-life’ of the initial population) (Ruiz-Rus et al. Reference Saffman2022).

Figure 13 shows that the axial evolution of the normalised diameter ratio, $\mathcal{D}^*$ , as a function of residence time ( $t$ ), demonstrates a universal trend. Symbols represent experimental data measured at different bulk velocities ( $V$ = 6.1–8.4 m s⁻¹) and radial positions (centreline and near-wall) for initial void fractions of $\phi = 0.5\,\%$ , $1\,\%$ and $2\,\%$ in figure 13. The curve in each panel is the best-fit exponential (4.15), illustrating the rapid rise from unity at $t=0$ toward an asymptote of about $1.15$ . A higher void fraction yields a steeper rise (smaller $\tau$ ), consistent with more frequent bubble coalescence at larger $\phi$ . The $\mathcal{D}^{*}$ plotted versus time $t$ , for all $V$ measurements, near-wall and centreline, collapse onto a single normalised curve at fixed initial $\phi$ (see figure 13), implying a degree of universality in the coalescence dynamics. The evolution of $\mathcal{D}^{*}$ is therefore controlled primarily by $\phi$ and $t$ , and is effectively independent of radial position. The influence of different bulk velocities is absorbed by the $t$ scaling – faster flows simply convect bubbles further in a given time, but the coalescence progression per unit time is unchanged. Such behaviour suggests that the coalescence process is governed by the cumulative interaction time between bubbles, and that our normalisation captures the essential time scale of those interactions. In our data, the master curve is well described by the above exponential fit, reinforcing the notion that a first-order kinetics governs the axial broadening of the bubble size distribution.

Figure 14.

Transient variation of $\mathcal{D}^*$ and its deviation from equilibrium with residence time.

Each master curve is, however, parameterised by the initial void fraction $\phi$ . The relaxation length (or e-folding length) $\tau$ in the exponential fit depends strongly on $\phi$ , whereas other parameters (e.g. bulk velocity, radial position) have negligible influence on $\tau$ . Quantitatively, increasing $\phi$ from 0.5 % to 2 % reduces $\tau$ from 131 to 76.7 (as inferred from the curves in figure 14), indicating that bubble populations with more initially crowded conditions coalesce and relax to their new size distribution much sooner. This indicates that $\tau$ is essentially an increasing function of the initial bubble spacing (or inverse function of $\phi$ ). Importantly, $\tau$ in our flow appears to be independent of the mean velocity $V$ , consistent with the notion that, in the reference frame of the flow, the intrinsic coalescence time scale is set by bubble collision kinetics (dictated by $\phi$ and turbulence levels) rather than convective transport.

4.3. Void fraction ( $\phi$ )

In order to explore the radial variation of the local void fraction, the observation window was divided into coloured regions (see figure 2 b), corresponding to near-wall (red) and central (blue) zones. For each region, several hundred 2D images were analysed to extract the bubble size distribution, repeating until the results became statistically stationary. The local void fraction in each region was then determined as the volume-weighted mean bubble volume per image, normalised by the total region volume. The local void fraction is denoted by $\phi _{wall}$ for near-wall and $\phi _{centre}$ for centreline. While these averages are assumed uniform within each coloured zone, minor radial variations may exist (Bunner & Tryggvason Reference Burns, Frank, Hamill and Shi2002, Reference Burns, Frank, Hamill and Shi2003; Lu & Tryggvason Reference Luo2006; Balachandar & Eaton Reference Balachandar and Eaton2010; Du Cluzeau, Bois & Toutant Reference Eggels, Unger, Weiss, Westerweel, Adrian, Friedrich and Nieuwstadt2019).

Figure 15.

Variation of $\phi _{\textit{wall}}$ and $\phi _{\textit{centre}}$ with axial position ( $\mathcal{L}$ ). Results are shown for (a) $\phi$ = 0.5 %, (b) $\phi$ = 1 %, (c) $\phi$ = 2 %.

Figure 15 shows the change in local void fraction $\Delta \phi = {\phi _{\textit{centre/wall}} - \phi _0}/{\phi _0}$ at the centre and wall in the axial direction of the duct across three $\phi$ ( $0.5\,\%, 1\,\%, 2\,\%$ ). At the inlet ( $\mathcal{L} = 0$ ), the void fraction is nearly uniform, with both centreline and near-wall values equal to the cross-sectional average (input value). As the flow develops downstream, a distinct segregation emerges: the near-wall void fraction decreases steadily, while the centreline void fraction increases commensurately, indicating progressive bubble accumulation in the core. For the lowest $\phi$ of $0.5\,\%$ , this redistribution is gradual but substantial; by the end of the measured domain, the centreline void fraction rises to ${\sim} 1.6 \times$ its inlet value, while the near-wall void fraction drops to ${\sim} 0.4\times$ its inlet value, resulting in a growing centre-wall void fraction difference ( $\Delta \phi$ ) that increases nearly linearly with time and shows no plateau. At higher $\phi$ ( $= 2\,\%$ ), a similar trend is observed with larger absolute changes: the centreline void fraction rises from 2 % to nearly 3 %, while the near-wall region is depleted to about 1 %. However, when normalised by the initial average, the relative enrichment at the core is slightly less pronounced at higher $\phi$ , indicating that bubble–bubble interactions and enhanced mixing at high void fractions can temper the relative segregation (Bunner & Tryggvason Reference Burns, Frank, Hamill and Shi2002, Reference Burns, Frank, Hamill and Shi2003; Burns et al. Reference Chan, Klaseboer and Manica2004b ; Dabiri, Lu & Tryggvason Reference Deane and Stokes2013).

The influence of mean flow velocity is also evident. At lower bulk velocities ( $V \approx 6.1 \ {\textrm m\,s}^{-1}$ ), the transition toward a core-peaked profile occurs more rapidly and distinctly, with bubbles segregating inward sooner and to a greater extent. At higher velocities ( $V \approx 8.4 \ {\textrm m\,s}^{-1}$ ), void fraction profiles remain closer to uniform for longer axial distances and the centre-wall divergence grows more slowly. This arises because stronger turbulence at higher $V$ is sustained farther downstream, keeping bubbles smaller and better mixed, and thus delaying the onset of segregation (Lu & Tryggvason Reference Luo2006; Burns et al. Reference Chan, Klaseboer and Manica2004b ; Balachandar & Eaton Reference Balachandar and Eaton2010). At lower $V$ , turbulence decays more quickly, enabling earlier bubble coalescence and growth, which amplifies the effect of lift and buoyancy on radial migration (Bunner & Tryggvason Reference Burns, Frank, Hamill and Shi2003; Lu & Tryggvason Reference Luo2006; Balachandar & Eaton Reference Balachandar and Eaton2010).

These observations represent a significant departure from classical duct/pipe flow studies, which consistently report wall-peaked void fraction profiles in turbulent flows (Delnoij et al. Reference Delnoij, Kuipers, van Swaaij and Akker1997; Hibiki & Ishii Reference Hinze1999). In these flows, small bubbles – owing to positive lift coefficients – accumulate near the wall, producing pronounced wall peaks. For our flow, however, even with small bubble sizes ( $d_{32} \sim$ 200–600 ${\unicode{x03BC}}$ m), a strong and persistent core peak emerges. This divergence from established patterns indicates that bubble migration mechanisms in high Reynolds number, vertical bubbly flows are fundamentally distinct. A key novelty of the present study is the demonstration that core-peaking arises for bubbles an order of magnitude smaller than the classical lift-reversal diameter (typically 5–6 mm in air–water systems) (Tomiyama et al. Reference van den, Thomas, Luther, Lathrop and Lohse2002; Hidman et al. Reference Hinze2022). Traditional lift models predict that only sufficiently large, deformable bubbles undergo a reversal in the lateral lift force, driving them from the wall toward the core and generating centre-peaked profiles. Our experiments, however, show that even micron-scale bubbles experience strong inward migration under diminished turbulence, challenging established predictions and aligning with recent DNS results that reveal early lift reversal in deformable bubbles at lower Reynolds numbers (Tomiyama Reference Valente and Vassilicos1998; Du Cluzeau et al. Reference Eggels, Unger, Weiss, Westerweel, Adrian, Friedrich and Nieuwstadt2019; Hidman et al. Reference Hinze2022).

Studies by Adoua, Legendre & Magnaudet (Reference Adoua, Legendre and Magnaudet2009), Shi (Reference Soligo, Roccon and Soldati2024) show that the direction of bubble migration is governed by the balance between the shear-induced lift force ( $F_L$ ) and the pressure-gradient (or ‘potential’) force ( $F_P$ ). The pressure-gradient force arises from the curvature of the mean carrier-flow velocity profile and always acts toward the centreline in duct and pipe flows. This is consistent with the transverse force decomposition in Shi & Magnaudet (Reference Sreenivasan2020) and Shi (Reference Soligo, Roccon and Soldati2024), where the potential component never reverses sign. The lift force, in contrast, may change sign depending on the lift coefficient $C_L$ . In a developing turbulent duct, we use the turbulence–Weber-based lift coefficient proposed by Salibindla et al. (Reference Sciacchitano and Wieneke2020), which is suitable for millimetric bubbles with moderate deformability. The corresponding lift force is

(4.17) \begin{align} F_L(r)=\rho _\ell \,C_L({We})\,|\Delta w(r)|\,|\gamma (r)|\,V_b, \end{align}

where $\Delta w(r)$ is the local slip velocity, $\gamma (r)=|\partial u_z/\partial r|$ is the local shear rate and $V_b$ is the bubble volume. The sign of $F_L$ is therefore controlled by the sign of $C_L$ , which we evaluate using the correlation of Salibindla et al. (Reference Sciacchitano and Wieneke2020):

(4.18) \begin{align} C_L= \begin{cases} 2.671\,{We}^{3/5}, & {We}\lt 0.1,\\[2pt] 1.25-1.608\,{We}^{3/5}, & 0.1\le {We}\lt 0.9,\\[2pt] -0.23, & {We} \ge 0.9. \end{cases} \end{align}

Here ${We}$ is evaluated from (1.1). In all experimental cases and at all axial locations, a substantial fraction of bubbles falls in the regime $\text{We} \gt 0.6$ , for which the lift coefficient becomes non-positive ( $C_L \lesssim 0$ ), implying reversal of the lateral lift force.

Classically, core peaking is not expected for such small bubbles. For instance, Lu & Tryggvason (Reference Luo2006) found that nearly spherical bubbles remain wall peaked downstream in highly turbulent flows, and that a lift reversal – leading to core accumulation – requires ${Eo} \gg 1$ (i.e. significantly deformable bubbles of the order of 5–6 mm). Here, ${Eo} = {(\rho _\ell -\rho _g)\,g\,d^{2}}/{\gamma }$ is the Eötvös number, where $g$ is gravity and $\rho _g$ is air density, representing the ratio of buoyancy to surface-tension forces. However, their study does not address cases where ${Eo} \ll 1$ , particularly in highly decaying flows, where turbulent dispersion has weakened to the extent that even a modest positive $C_L$ drives bubble accumulation toward the core. As in this case, turbulence decays by over 90 % along the duct (see the monotonic drop in $k/V^2$ – see figure 6), and turbulent dispersion is dramatically weakened, allowing lift and wall forces to dominate and drive bubbles toward the core (Lu & Tryggvason Reference Luo2006; Burns et al. Reference Calado and Balaras2004a ).

The observed transition from uniform to core-peaked void fraction profiles arises from three key mechanisms: (i) as $\mathcal{I}$ diminishes (by ${\sim}$ 90 %), turbulent dispersion is weakened, losing its ability to homogenise the bubble distribution and allowing other forces to dominate; (ii) even micron-scale bubbles experience sufficient inward lift under low turbulence, overturning the conventional view that a critical bubble size ( ${\sim}$ 5–6 mm) is necessary for lift reversal and core peaking, as recently confirmed by DNS (Du Cluzeau et al. Reference Eggels, Unger, Weiss, Westerweel, Adrian, Friedrich and Nieuwstadt2019; Hidman et al. Reference Hinze2022); and (iii) as lift reversal takes hold, bubbles migrate toward the duct centreline, where upward liquid velocity is highest and drag is minimised, energetically favouring central accumulation (Lu & Tryggvason Reference Luo2006). This drag-reduction mechanism is further reinforced by wall lubrication, which prevents bubble accumulation near the wall and maintains a depleted boundary region.

These findings not only depart from established wall-peaked paradigms seen in prior studies, but also suggest that the interplay of lift, wall-lubrication and turbulent-dispersion forces govern radial void fraction profiles even under extreme conditions. The present results thus fill an important gap in the literature, demonstrating persistent and robust core peaking even for micron-scale bubbles, a regime not previously reported, and provide new insight into dispersed phase redistribution mechanisms in turbulent multiphase flows.

In particular, the monotonic decrease of $\phi$ near the wall and its growth near the centre with increasing $\mathcal{L}$ raise important questions about the mechanisms driving migration from the wall toward the centre. This behaviour warrants further investigation by examining how radial distributions evolve under varying flow conditions and system parameters.

Figure 16.

Transient radial variation of void fraction ( $\phi (r)$ ) in a decaying turbulent regime. Parameters: sampling instants $t_i$ are defined by $t_i=\mathcal{L}_i/V_1$ ( $i=0,\ldots ,5$ ). Bulk velocities: $V_1=6.1$ , $V_2=7.4$ , $V_3=8.4$ m s^–1. Axial locations: $\mathcal{L}_0=0$ , $\mathcal{L}_1=3.64$ , $\mathcal{L}_2=12.73$ , $\mathcal{L}_3=21.82$ , $\mathcal{L}_4=26.36$ , $\mathcal{L}_5=30.91$ . Results are shown for (a) $\phi (r)/\phi _0$ with $\mathcal{L}$ at $\phi$ = 0.5 % and $V=6.1$ m s⁻¹, (b) $\phi (r)/\phi _0$ with $\phi$ at $V=6.1$ m s⁻¹, (c) $\phi (r)/\phi _0$ with $V$ at $\phi$ = 0.5 % & $\mathcal{L}$ .

As bubbles migrate radially from the wall region (red) toward the duct centre (blue) (see figure 2 b), the available cross-sectional area decreases. This geometric constraint leads to an amplification of local $\phi$ near the centre. Consequently, the radial variation of the local void fraction is inherently nonlinear. To quantitatively characterise this evolution, the measured radial void fraction profiles at each axial location ( $\mathcal{L}$ ), void fraction ( $\phi$ ) and bulk velocity ( $V$ ) were fitted to a Gaussian function (shown in figure 16):

(4.19) \begin{equation} \frac {\phi (r)}{\phi _o} = \varPhi _i \exp \left (-\frac {(r/R - R_o)^2}{2\sigma ^2}\right ) .\end{equation}

Here $\phi _0$ is the cross-sectional mean, $\varPhi _i$ is the peak amplitude (centreline void fraction), $R_o$ is the mean position (fixed at zero by symmetry) and $\sigma$ is the standard deviation describing the distribution width. The fitted parameters of figure 16 for all test cases are given in table 3. Figure 16 shows that near the duct inlet ( $\mathcal{L}=0$ ), the void fraction profile is nearly flat ( $\varPhi _i=1$ , $\sigma \rightarrow \infty$ ), confirming uniform bubble dispersion. Downstream, as $\mathcal{L}$ increases, $\varPhi _i$ rises and $\sigma$ decreases monotonically – indicating enhanced centreline peaking and profile narrowing due to turbulence decay, increased migration time and flow development. These trends are most pronounced at lower $\phi$ and $V$ , as shown in figure 16(a).

Table 3.

Fitted parameters (refer to figure 16) of the normalised radial void fraction profiles for three different cases.

It is interesting to note the effects of $\phi$ and $V$ (figure 16 b,c). For moderate $\phi$ ( $\leq 1\,\%$ ), $\varPhi _i$ and $\sigma$ remain nearly unchanged with increasing $\phi$ , indicating self-similar profile shapes. However, at $\phi =2\,\%$ , the distribution broadens ( $\sigma$ increases) and the centreline amplitude decreases ( $\varPhi _i$ drops). This broadening arises from increased bubble–bubble interactions and, crucially, from stronger turbulent dispersion, which depends directly on the radial gradient $\Delta \phi (r)$ . The dispersion flux is proportional to $-\boldsymbol{\nabla }\phi _g$ ; so larger $\Delta \phi (r)$ (i.e. greater radial contrast) enhances outward mixing and inhibits further accumulation at the centre (Hosokawa & Tomiyama Reference Jassal, Song and Schmidt2009). This self-regulating mechanism, where increased peaking amplifies dispersion, prevents unlimited core enrichment and is rarely captured in classical drift-flux models (Burns et al. Reference Calado and Balaras2004a ; Ooms, Kewen & Mudde Reference Pope2007; Hosokawa & Tomiyama Reference Jassal, Song and Schmidt2009). A few recent DNS studies further confirm that bubble-driven dispersion follows Kolmogorov-like scaling and is modulated by radial gradients in the void fraction (Balachandar & Eaton Reference Balachandar and Eaton2010; Elghobashi Reference Evans, Jameson and Atkinson2019). Similarly, increasing $V$ at fixed $\phi$ leads to lower $\varPhi _i$ and higher $\sigma$ , consistent with higher dissipation, stronger turbulent mixing and reduced migration time. Thus, broader, less-peaked distributions at high $\phi$ or $V$ reflect the dominance of dispersion and bubble interactions over migration (Tomiyama et al. Reference van den, Thomas, Luther, Lathrop and Lohse2002; Du Cluzeau et al. Reference Eggels, Unger, Weiss, Westerweel, Adrian, Friedrich and Nieuwstadt2019; Hidman et al. Reference Hinze2022).

The Gaussian parameterisation offers more than a compact description: it provides a direct link between observable profile shape and the balance between migration and dispersion mechanisms. The fitted $\sigma$ quantifies the effectiveness of turbulent dispersion for given flow conditions, while $\varPhi _i$ encodes the extent of core accumulation. These parameters are directly relevant for the calibration and validation of turbulence-dispersed two-phase flow models, including gradient-based closures (Burns et al. Reference Calado and Balaras2004a ), and enable systematic comparison across modelling approaches. By relating $\sigma$ to $\mathcal{I}$ and $\Delta \phi (r)$ , this approach enables predictive scaling and cross-regime analysis, bridging detailed experiments, numerical simulations and practical engineering models.

Table 4.

List of symbols and their descriptions used in the paper.