3.1.1. Bubble density distribution

Bubbles are initially randomly released into the fully developed turbulent flow. Due to turbulent transport and forces acting on the bubbles, they tend to accumulate in specific flow regions. Figure 4(a) shows the time evolution of the Shannon entropy $S$ , which quantifies the spatial non-uniformity of bubbles in the wall-normal direction:

(3.1) \begin{align} S&=\zeta /\zeta _{\textit{max}},\nonumber \\ \zeta &=-\sum _{j=1}^{N_{y}}\frac {N_{b,j}}{N_{b}}\ln \frac {N_{b,j}}{N_{b}}, \end{align}

where $N_y=400$ is the number of wall-parallel layers, $N_{b,j}$ is the number of bubble in the $j{\textrm{th}}$ layer, $N_b$ is the total number and $\zeta _{max}=\ln N_y$ is the normalisation factor. Once $S$ reaches convergence, the bubble distribution is considered to be statistically steady. For the smallest 120 ${\unicode{x03BC}} \textrm{m}$ bubbles, stabilisation occurs at $t\approx 2\times 10^4$ , with a Shannon entropy of $S\approx 0.75$ , indicating strong non-uniformity. For the 125 and 130 ${\unicode{x03BC}} \textrm{m}$ bubbles, stabilisation is achieved more rapidly ( $t\approx 1\times 10^4$ ), with final $S$ of approximately 0.94 and 0.99, respectively. Larger bubbles $(d\geqslant 140\,{\unicode{x03BC}} \textrm{m})$ stabilise even faster, at $t\sim O(10^3)$ , and exhibit a Shannon entropy close to unity, suggesting an almost uniform wall-normal distribution.

Figure 3. Comparison of bubble velocity statistics with the results of Molin et al. (Reference Molin, Marchioli and Soldati2012): (a) mean streamwise bubble velocity profile $\langle u^+\rangle$ ; (b) r.m.s. of bubble velocity fluctuation in the streamwise $ u^{\prime +}_{b,rms}$ , wall-normal $v^{\prime +}_{b,rms}$ and spanwise $w^{\prime +}_{b,rms}$ directions.

Figure 4. (a) Time evolution of Shannon entropy $S$ for bubble diameter range from 120 $\rm {\unicode{x03BC}} \textrm{m}$ to 240 ${\unicode{x03BC}} \textrm{m}$ ; (b) bubble density distribution $c$ , normalised by the bulk density $c_0$ .

The density distribution of bubbles after reaching a steady state is shown in figure 4(b). The mean bubble concentration is defined as $c=\Delta N_b/ \Delta V$ , where $\Delta N_b$ is the number of bubbles within a thin wall-parallel layer of volume $\Delta V$ at a specific $y^+$ , normalised by the bulk-averaged concentration $c_0$ .

Two distinct wall-normal bubble distribution patterns are observed based on bubble diameter. The first pattern corresponds to the smaller bubble $(d\leqslant 130\,{\unicode{x03BC}} \textrm{m})$ , which exhibits strong wall accumulation with a significant number of bubbles confined on the wall. This pattern aligns with the ‘Type I distribution profile’ in the experiment study by Felton & Loth (Reference Felton and Loth2001). The smallest bubbles (120 ${\unicode{x03BC}} \textrm{m}$ ) exhibit a peak concentration of $c/c_{0\ max}\approx 800$ near the lowest height $y^+\approx d^+/2$ , primarily due to the bubble–wall collision. Beyond $y^+\gt 10$ , the concentration stabilises at $c/c_0\approx 0.55$ , indicating that nearly half of these bubbles reside within the near-wall region ( $y^+\lt 10$ ). Strictly speaking, the peak concentration corresponds to a local value of $c=2.4\times 10^{-2}$ , which falls within the regime of four-way coupling (Balachandar & Eaton Reference Balachandar and Eaton2010), where bubble–bubble interactions should be taken into account. However, this extremely high concentration ( $c\gt 1.0\times 10^{-3}$ ) is confined to a very thin layer near the wall ( $0.4\lt y^+\lt 0.8$ ), which is not the primary region of interest in the present study. Moreover, given that the inertia of these bubbles is negligible ( $St\ll 1$ ), the effect of bubble–bubble interactions is expected to be very limited and will not spread to the outer region ( $y^+\gt 1$ ), which is different from the case of inertial particles. Therefore, neglecting the four-way coupling is not expected to significantly affect the main conclusions of this study.

As the bubble diameter increases $(d\lt 140\,{\unicode{x03BC}} \textrm{m})$ , the near-wall peak concentration decreases sharply, while a second peak gradually emerges. For the $130\, {\unicode{x03BC}} \textrm{m}$ bubbles, the second peak appears at $y^+\approx d^+$ and becomes comparable to the near-wall peak. The second distribution pattern is observed for larger bubbles $(d\geqslant 140\,{\unicode{x03BC}} \textrm{m})$ , which do not accumulate in the vicinity of the wall, but instead exhibit a peak concentration at $y^+\approx 2.5\sim 3.0d^+$ . As the bubble diameter increases, the $c/c_{0\ max}$ decreases and eventually stabilises. This pattern shows some similarity with the previous experimental and numerical studies, despite the much smaller bubble diameters considered here (Hibiki et al. Reference Hibiki, Goda, Kim, Ishii and Uhle2004; Lu & Tryggvason Reference Lu and Tryggvason2008).

Similar wall-accumulation behaviour is widely observed in inertial particle-laden turbulent flows due to the turbophoresis phenomenon (Reeks Reference Reeks1983; Soldati & Marchioli Reference Soldati and Marchioli2009), which is particularly pronounced at intermediate Stokes number ( $St=10{-}50$ ), where the particle response time is comparable to the characteristic flow time scale (Picano et al. Reference Picano, Sardina and Casciola2009). However, bubbles generally exhibit much lower Stokes numbers ( $St\sim O (10^{-2} )$ ), making turbophoresis an unlikely primary accumulation mechanism. In upward turbulent channel flow, the positive streamwise slip velocity generates a shear-induced lift force directed towards the wall, promoting wall-ward bubble migration and resulting in the substantial wall-accumulation (Giusti et al. Reference Giusti, Lucci and Soldati2005; Molin et al. Reference Molin, Marchioli and Soldati2012). When bubbles are very close to the wall, the wall-lift force acts in the opposite direction, preventing them from directly colliding with the wall (Antal, Lahey & Flaherty Reference Antal, Lahey and Flaherty1991; Kim & Kim Reference Kim and Kim2022).

Figure 5(a) presents the mean force components acting on bubbles in the wall-normal direction. The dominant forces in the wall-normal direction are the lift, wall-lift and drag force. As for other force components, the Basset force $\langle F_{y,B}\rangle ^+$ is minor, constituting approximately $1/3$ of drag force $\langle F_{y,D}\rangle ^+$ near the wall and decaying to zero as $y^+$ increases; the pressure gradient force $\langle F_{y,P}\rangle ^+$ and added mass force $\langle F_{y,A}\rangle ^+$ are nearly zero across all heights due to the low Stokes number. The wall-lift force $\langle F_{y,W}\rangle ^+$ acts exclusively in the wall-normal direction and is always positive. Since the slip velocity $u_s^+$ remains nearly constant (as will be discussed in detail in § 3.2.1), the magnitude of wall-lift force primarily depends on the wall-lift coefficient $C_W$ , which increases sharply as the wall distance decreases and vanishes beyond $y^+\gt 10$ . As the bubble diameter increases, $\langle F_{y,W}\rangle ^+$ also rises. The lift force $\langle F_{y,L}\rangle ^+$ is more complicated compared with the wall-lift force. Its maximum occurs at the lowest positions, where the spanwise shear rate $\partial u^+/\partial y^+$ is greatest. As $y^+$ increases, it decreases and approaches zero in the logarithmic layer. For the 120 ${\unicode{x03BC}} \textrm{m}$ bubbles, $\langle F_{y,L}\rangle ^+$ decreases gradually across the channel, while for the larger bubbles ( $d\geqslant 140\,{\unicode{x03BC}} \textrm{m}$ ), inflection points are located near the peak concentration ( $y^+\approx 3$ ). Above the inflection points, $\langle F_{y,L}\rangle ^+$ changes slowly, whereas below the inflection point, it increases rapidly. This is due to the bubble’s preference for regions with higher spanwise shear rates, which will be elaborated on in § 3.1.2.

Figure 5. (a) All the force components acting on the bubble in the wall-normal direction; (b) sum of lift force and wall-lift force in wall-normal direction $\langle F_{y,L+W}\rangle ^+$ . In panel (a), $\langle F_{y,D}\rangle ^+$ , drag; $\langle F_{y,L}\rangle ^+$ , lift; $\langle F_{y,W}\rangle ^+$ , wall-lift; $\langle F_{y,P}\rangle ^+$ , pressure-gradient; $\langle F_{y,A}\rangle ^+$ , add-mass; $\langle F_{y,B}\rangle ^+$ , Basset force. Solid lines, 120 ${\unicode{x03BC}} \textrm{m}$ bubble; dashed lines, $140\,{\unicode{x03BC}} \textrm{m}$ bubble; dash-dotted lines, 180 ${\unicode{x03BC}} \textrm{m}$ bubble; and dotted lines, $240\,{\unicode{x03BC}} \textrm{m}$ .

In figure 5(b), we focus on the total lift force, defined as the sum of lift force and wall-lift force $\langle F_{y,L+W}\rangle ^+=\langle F_{y,L}\rangle ^++\langle F_{y,W}\rangle ^+$ , which directly governs the bubble’s wall-normal dynamics and determines their concentration distribution. The drag force $\langle F_{y,D}\rangle ^+$ is nearly the inverse of $\langle F_{y,L+W}\rangle ^+$ . For all bubble sizes, $\langle F_{y,L+W}\rangle ^+$ reaches a minimum around $y^+=7$ , while it shows different patterns for different bubble diameters. For 120 ${\unicode{x03BC}} \textrm{m}$ bubbles, $\langle F_{y,L+W}\rangle ^+$ remains negative throughout the channel, resulting in strong wall accumulation. The negative value of $\langle F_{y,L+W}\rangle ^+$ at the lowest position is associated with the bubble–wall collision. When the collision occurs, the lift force increases suddenly. This phenomenon will be further explained in § 3.1.2.

3.1.2. Effect of lift and wall-lift forces

The wall-normal force on bubbles affects not only their distribution in the wall-normal direction, but also their preferential concentration in the wall-parallel plane and their transportation behaviours in the near-wall region. In this section, we first explore the preferential concentration of bubbles based on the instantaneous snapshots of bubble positions, followed by an analysis of the underlying mechanisms.

Figure 6 illustrates the preferential concentration of 120 and 180 ${\unicode{x03BC}} \textrm{m}$ bubbles at different $y^+$ . The contours show the instantaneous streamwise velocity fluctuation, $u^{\prime +}=u^+-\langle u\rangle ^+$ , where $u^+$ and $\langle u\rangle ^+$ are the instantaneous local velocity and mean velocity, respectively. Circles marking bubbles within a thin layer around this height. Bubbles of different diameters exhibit distinct patterns of preferential concentration. As shown in figure 6(a, c), the smallest 120 ${\unicode{x03BC}} \textrm{m}$ bubbles tend to accumulate in the regions with low streamwise velocities (referred to as low-speed regions) within the viscous sublayer, which corresponds to the low-speed streaky structure, a typical near-wall turbulent coherent structure, and is usually accompanied by the ejection motion. In the viscous sublayer, these low-speed streaks are just forming, which are elongated and relatively straight near the wall (Kim, Kline & Reynolds Reference Kim, Kline and Reynolds1971). In figure 6(a), at $y^+=2$ , the concentration of 120 ${\unicode{x03BC}} \textrm{m}$ bubbles is extremely high, forming long and narrow bubble streaks that align with the fluid streaks, but extend over much longer lengths, approximately $1000{-}2000\delta _\nu$ , nearly half the computational domain length. This is similar to the particle streaks in the inertial particle-laden flows (Sardina et al. Reference Sardina, Picano, Schlatter, Brandt and Casciola2011; Jie et al. Reference Jie, Cui, Xu and Zhao2022). As $y^+$ increases, the preferential concentration of 120 ${\unicode{x03BC}} \textrm{m}$ bubbles gradually diminishes.

Figure 6. Instantaneous streamwise velocity fluctuation $u^{\prime +}$ in the wall-parallel planes, along with corresponding bubble distributions: (a, c, e) 120 ${\unicode{x03BC}} \textrm{m}$ bubbles; (b, d, f) 180 ${\unicode{x03BC}} \textrm{m}$ bubbles. (a–b) $y^+=2$ ; (c–d) $y^+=5$ ; (e–f) $y^+=20$ . Contours, instantaneous streamwise velocity fluctuations of fluid $u^{\prime +}$ , normalised by the friction velocity; circles, bubble positions; circle colour, wall-normal velocity of bubbles $v_b^+$ .

In contrast, the larger 180 ${\unicode{x03BC}} \textrm{m}$ bubbles exhibit unique clustering behaviour near the wall. As shown in figure 6(b), below the peak of concentration $(y^+\lt 3)$ , these larger bubbles tend to cluster in the regions with high streamwise velocities (referred to as high-speed regions), which are usually caused by the sweep events that bring high-speed fluid from the outer region towards the wall (Corino & Brodkey Reference Corino and Brodkey1969). The structure of the high-speed regions is fragmented, preventing the formation of continuous bubble streaks. However, as shown in figure 6(d), above the concentration peak, while still within the viscous sublayer, 180 ${\unicode{x03BC}} \textrm{m}$ bubbles tend to accumulate in the low-speed streaks, similar to the 120 ${\unicode{x03BC}} \textrm{m}$ bubbles. Figures 6(e) and 6(f) show that in the middle of the buffer layer at $y^+=20$ , both kinds of bubbles exhibit only a weak preference for the low-speed regions. Above the buffer layer, the bubble distribution becomes nearly uniform (not shown here).

The wall-normal bubble velocity $v_b^+$ indicates the transport characteristics and is represented by the colour of the circles in figure 6. As shown in figure 6(b), at $y^+=2$ , nearly all the 180 ${\unicode{x03BC}} \textrm{m}$ bubbles are located in the high-speed sweep regions. Bubbles within relatively stronger high-speed regions (stronger sweep events) exhibit negative $v_b^+$ , while those in relatively weaker high-speed regions show positive $v_b^+$ . This wall-vicinity preferential concentration behaviour will be discussed in detail in § 3.1.3. In figure 6(a), although the 120 ${\unicode{x03BC}} \textrm{m}$ bubbles accumulate in the low-speed ejection regions, opposite to the 180 ${\unicode{x03BC}} \textrm{m}$ bubbles, their wall-normal velocity $v_b^+$ still correlates negatively with $u_{@b}^+$ . This indicates that the wall-normal transportation of bubbles does not simply follow the coherent structures as it does for inertial particles.

Figure 7 shows the mean relative velocity of the flow at the bubble location, $\langle \Delta u_{@b,i}\rangle ^+=\langle u_{@b,i}\rangle ^+-\langle u_{i}\rangle ^+$ , which represents the difference between the mean fluid velocity at bubble locations and that over the whole plane. This serves as a quantitative measure of preferential concentration of bubbles in the velocity sample plane ( $u^\prime {-}v^\prime$ ). Focusing first on the wall-normal direction, figure 7(a) shows the mean relative wall-normal velocity at bubble location $\langle \Delta v_{@b}\rangle ^+$ , which equals $\langle v_{@b}\rangle ^+$ . A comparison with figure 5(b) reveals that $\langle \Delta v_{@b}\rangle ^+$ is nearly inversely proportional to the mean total wall-normal lift force $\langle F_{y,L+W}\rangle ^+$ . This relationship arises because the total lift force is balanced by the drag force, leading to $\langle F_{y,L+W}\rangle ^+\propto \langle v_{s}\rangle ^+$ . Once the two-phase system reaches a steady state, the mean wall-normal bubble velocity $\langle v_{b}\rangle ^+=\langle v_{@b}\rangle ^++\langle v_{s}\rangle ^+$ approaches zero due to mass conservation, which means that $\langle v_{@b}\rangle ^+\approx -\langle v_{s}\rangle ^+$ . Accordingly, we obtain $\langle F_{y,L+W}\rangle ^+\propto -\langle \Delta v_{@b}\rangle ^+$ .

For the 120 ${\unicode{x03BC}} \textrm{m}$ bubbles, $\langle F_{y,L+W}\rangle ^+$ is always negative, resulting in $\langle \Delta v_{@b}\rangle ^+\gt 0$ , indicating that these bubbles prefer accumulating on the outward regions ( $v\gt 0$ ). Since the instantaneous slip velocity $v_{s}^+$ is more likely to be negative, the bubbles must reside in the outward-moving region ( $v\gt 0$ ), where the positive fluid velocity $v_{@b}^+$ can offset the slip velocity and maintain its wall-normal height. Conversely, if it is located in the inward-moving region ( $v\lt 0$ ), it will move towards the wall more rapidly.

Figure 7. Mean relative fluid velocity at bubble location in the (a) wall-normal direction $\langle \Delta v_{@b}\rangle ^+$ and (b) streamwise direction $\langle \Delta u_{@b}\rangle ^+$ . Solid lines, cases of two-way coupling; dashed line, case of one-way coupling for the 180 ${\unicode{x03BC}} \textrm{m}$ bubbles.

It has been mentioned already that the ejection and sweep events are two typical turbulent coherent structures, where ejections correspond to low-speed streaks that migrate away from the wall $(u^\prime \lt 0,v^\prime \gt 0)$ and are located in the second quadrant (Q2) of the $u^\prime {-}v^\prime$ velocity fluctuation plane. Sweeps, however, correspond to high-speed fluid moving towards the wall and are located in the fourth quadrant (Q4) (Corino & Brodkey Reference Corino and Brodkey1969). The 120 ${\unicode{x03BC}} \textrm{m}$ bubbles preferentially cluster in the outward regions ( $v\gt 0$ ), which corresponds to the ejection regions with low-speed streaks. Therefore, $\langle \Delta u_{@b}\rangle ^+$ in figure 7(b) is negative. Similarly, for the larger 180 ${\unicode{x03BC}} \textrm{m}$ bubbles, when they locate below $y^+\approx 2$ , the wall-lift force dominates in the wall-normal direction, leading to $\langle F_{y,L+W}\rangle ^+\gt 0$ and $\langle \Delta v_{@b}\rangle ^+\lt 0$ . This indicates a preferential concentration in the inward region ( $v\lt 0$ ), which corresponds to the high-speed sweep regions, where $\langle \Delta u_{@b}\rangle ^+\gt 0$ . Above this height, the total lift force becomes negative, causing the 180 ${\unicode{x03BC}} \textrm{m}$ bubbles to accumulate in the ejection regions, similar to the 120 ${\unicode{x03BC}} \textrm{m}$ bubbles.

It is interesting to see in figure 7 that the wall-normal range where bubbles prefer the high-speed region ( $\langle \Delta u_{@b}\rangle ^+\lt 0$ ) is higher and shorter than that of the outward region ( $\langle \Delta v_{@b}\rangle ^+\gt 0$ ). For example, in the case of 180 ${\unicode{x03BC}} \textrm{m}$ bubbles, $\langle \Delta u_{@b}\rangle ^+$ is negative only within $3.6\lt y^+ \lt 16$ , whereas $\langle \Delta v_{@b}\rangle ^+$ becomes positive at $y^+\gt 2.0$ . This discrepancy can be attributed to the two-way coupling between the bubbles and the fluid.

In the streamwise direction, the feedback force from bubble to fluid (mainly the buoyancy force) is always positive, accelerating the local flow around the bubble itself and other bubbles nearby. The dashed line in figure 7 shows the results for the 180 ${\unicode{x03BC}} \textrm{m}$ bubbles under the same conditions but in one-way coupling (without the feedback force on fluid). By comparing the two cases, it is evident that the feedback force increases $\langle \Delta u_{@b}\rangle ^+$ throughout the domain, with the maximum increment at $y^+\approx 4$ , where bubbles are highly concentrated and strongly clustered in low-speed regions. Here, $\langle \Delta u_{@b}\rangle ^+$ is slightly larger than zero above the buffer layer, although not strictly zero. The two-way coupling also exerts an indirect impact in the wall-normal direction. As seen in figure 7(a), it increases $\langle \Delta v_{@b}\rangle ^+$ within the viscous sublayer ( $y^+\lt 5$ ), while the outer regions remain relatively unaffected. This occurs because the feedback force increases the local streamwise velocity $u^+$ , which in turn amplifies the velocity gradient $\partial u^+/\partial y^+$ . As a result, the lift force increases, whereas the wall-lift force remains nearly unchanged. The increased lift force enhances the bubble’s preference for the outward regions, thereby increasing $\langle \Delta v_{@b}\rangle ^+$ .

Comparing figures 7(b) and 4(a), we observe that the peak concentration of bubbles in figure 4(a) is close to the zero point of $\langle \Delta u_{@b}\rangle ^+$ . This is because $\langle \Delta u_{@b}\rangle ^+=0$ implies that the number of bubbles in the high-speed and low-speed regions is nearly equal, indicating a balance between outward- and inward-moving bubbles. In the one-way coupling case, the zero points of both the $\langle \Delta u_{@b}\rangle ^+$ and $\langle \Delta v_{@b}\rangle ^+$ are almost identical and located at $y^+=3$ , where the balance between lift and wall-lift force determines the peak bubble concentration. Within the two-way coupling case, however, the zero point of $\langle \Delta v_{@b}\rangle ^+$ is lowered, while that of $\langle \Delta u_{@b}\rangle ^+$ is raised. This effect of two-way coupling leads to the peak bubble concentration no longer matching the balance point between lift force and wall-lift force.

It should be noted that, at the lowest height ( $y^+\approx 0.4$ ), $\langle \Delta u_{@b}\rangle ^+$ of the 120 ${\unicode{x03BC}} \textrm{m}$ bubbles becomes positive. This is attributed to bubble–wall collision, which causes the sign of $\langle v_{s}\rangle ^+$ to be reversed, from negative to positive. As a result, after collision, the 120 ${\unicode{x03BC}} \textrm{m}$ bubbles transiently migrate to high-speed sweep regions, maintaining a positive slip velocity. Because the bubbles have negligible inertia, the reversed slip velocity decays rapidly. Consequently, this preference for the high-speed region is confined to an extremely thin layer ( $0.4 \lt y^+ \lt 0.5$ ).

Figure 8 further quantitatively shows the preferential concentration of bubbles at various $y^+$ in the $u^\prime {-}v^\prime$ plane. The quadrant proportion of the whole flow field (or flow at bubble location) is defined as $\varPhi _i=N_{i}/N_{\textit{total}}$ , where $N_{i}$ is the number of grid points (or bubbles) located in the ith quadrant at a given height, and $N_{\textit{total}}=\sum _{j=1}^{4}N_{j}$ represents the total number of grid points (or bubbles) in all the four quadrants. Figure 8(a) illustrates the preferential concentration of 120 ${\unicode{x03BC}} \textrm{m}$ bubbles. It is observed that a higher proportion of bubbles resides in Q2 (ejection events) within the viscous sublayer and the lower buffer layer ( $y^+\lt 20$ ), while the proportions in Q4 (sweep events) are lower than in the whole flow field. The proportions in Q1 and Q3 show small differences compared with the flow field. At $y^+\gt 20$ , the proportions of bubbles in each quadrant become nearly identical to those of the flow, confirming that the bubbles are uniformly distributed above the near-wall region. Due to the bubble–wall collision, the proportions of Q1 and Q4 increase rapidly as $y^+$ decreases.

Figure 8. Fraction of the four quadrants in the $u^\prime {-}v^\prime$ plane: (a) 120 ${\unicode{x03BC}} \textrm{m}$ bubbles; (b) 180 ${\unicode{x03BC}} \textrm{m}$ bubbles. Solid lines, fluid at grid point in bubble-laden flow; dashed lines, fluid at bubble locations in bubble-laden flow; dotted lines, fluid at grid point in unladen flow.

Figure 8(b) shows four different preferential concentration patterns of the 180 ${\unicode{x03BC}} \textrm{m}$ bubbles. In figure 9, we select four typical heights representing these preferential patterns ( $y^+=1.5,2.5,5$ and $20$ ) and plot the joint probability density function (p.d.f.) of the relative fluid velocity at bubble location ( $u^{\prime +}_{@b},v^{\prime +}_{@b}$ ) (contours) and fluid velocity fluctuation of all grid points ( $u^{\prime +},v^{\prime +}$ ) for comparison (isoline), where $\boldsymbol{u}_{@b}^{\prime +}$ is defined as $\boldsymbol{u}_{@b}^{\prime +}=\boldsymbol{u}_{@b}^+-\langle \boldsymbol{u}\rangle ^+$ . As $y^+$ increases, four distinct preferential concentration patterns are observed. (i) In the vicinity of the wall ( $y^+\lt 2$ ), as seen in figure 8(b), bubbles tend to cluster in Q4 regions. As shown in figure 9(a), at $y^+=1.5$ , bubbles concentrate in the strong sweep event, with $\langle \Delta u_{@b}\rangle ^+\approx 3u_{\textit{rms}}^{\prime +}$ . As the bubbles approach the wall, the proportion of Q4 increases, reaching approximately $90\,\%$ at the lowest point ( $y^+\approx 1$ ), where nearly all bubbles are in the sweep events. (ii) As the height increases to $2\lt y^+\lt 3$ , as seen in figure 9(b), bubbles tend to concentrate in Q3 regions. As discussed before, this is attributed to the two-way coupling effect, a pattern that does not appear in the one-way coupling case. (iii) Above the concentrate peak ( $3 \lt y^+ \lt 15$ ), as shown in figure 8(b), the proportion of bubbles in Q2 regions increases rapidly, while the proportion of Q4 decreases. In figure 9(c), at $y^+=5$ , the peak of joint p.d.f. is at ( $u_{@b}^{\prime +}\approx -0.9u_{\textit{rms}}^{\prime +}$ , $v_{@b}^{\prime +}\approx 0.3v_{rms}^{\prime +}$ ), indicating that most bubbles are concentrated in strong ejection events. Although a small proportion of bubbles are still in the high-speed Q1 and Q4 regions with a combined proportion of less than $30\,\%$ , this makes the $\langle \Delta u_{@b}\rangle ^+$ increase to approximately $-0.2u_{\textit{rms}}^{\prime +}$ , away from the peak point. (iv) Above $y^+=15$ , as shown in figure 9(d), the bubbles are nearly uniformly distributed, with their quadrant proportion distribution closely matching that of the flow field. Similarly, for the smallest 120 ${\unicode{x03BC}} \textrm{m}$ bubbles, only the patterns (iii) and (iv) are observed (not shown here).

Figure 9. Joint p.d.f. of fluid velocity fluctuations $(u^{\prime +},v^{\prime +})$ and those at bubble location $(u^{\prime +}_{@b},v^{\prime +}_{@b})$ in the case of 180 ${\unicode{x03BC}} \textrm{m}$ bubbles: (a) $y^+=1.5$ ; (b) $y^+=2.5$ ; (c) $y^+=5$ and (d) $y^+=20$ . Isolines, fluid velocity; contours, fluid velocity at bubble location; dashed lines, mean relative fluid velocity at bubble location $(\langle \Delta u_{@b}\rangle ^+, \langle \Delta v_{@b}\rangle ^+)$ .

3.1.3. Preferential concentration in the wall vicinity

We have already discussed the influence of lift and wall-lift forces on the preferential concentration of bubbles at different heights across the channel. Furthermore, even at the same $y^+$ , bubbles experience varying magnitudes of lift and wall-lift forces depending on the local flow region, which further affects their preferential concentration behaviour. Here, we take 180 ${\unicode{x03BC}} \textrm{m}$ bubbles as an example to illustrate this mechanism. In figure 10(a), the joint p.d.f. of wall-normal lift force $F_{y,L}^+$ and streamwise velocity fluctuation at bubble locations $u_{@b}^{\prime +}$ shows a strong negative correlation in the vicinity of the wall ( $y^+=1.5$ ). This arises because, within the viscous sublayer, the spanwise shear rate $\partial u^+/\partial y^+$ is well proportional to the local streamwise velocity $u^+$ , while the slip velocity $u^+_s$ remains nearly constant. Consequently, the lift force follows $F_{y,L}^+\approx C_Lu^+_s\partial u^+/\partial y^+ \propto u_{@b}^{\prime +}$ , indicating that bubbles located in higher streamwise velocity regions experience a stronger wall-normal lift force that pushes them towards the wall. Meanwhile, the wall-lift force $F_{y,W}^+$ remains nearly unchanged (not shown here), leading to a reduction in the total lift force $F_{y,L+W}^+$ and, in turn, a lower wall-normal slip velocity $v^+_s$ . As shown in figure 10(b), $v^+_s$ exhibits a strong negative correlation with $u_{@b}^{\prime +}$ , while the wall-normal velocity fluctuation at bubble location $v^{\prime +}_{@b}$ also shows a weaker negative correlation (figure 9 a). These two effects collectively result in a strong negative correlation between the overall wall-normal bubble velocity $v^{\prime +}_b=v^{\prime +}_{@b}+v^+_s$ and $u_{@b}^{\prime +}$ (figure 10 c). This implies that even when all bubbles are located in high-speed regions, those in stronger high-speed regions will experience a larger lift force and are more likely to migrate towards the wall ( $v^+_b\lt 0$ ), while bubbles in relatively weaker high-speed regions tend to move away from the wall $(v^+_b\gt 0)$ .

Figure 10. Joint p.d.f. of bubble statistics in wall-normal direction and streamwise relative fluid velocity at bubble location $u_{@b}^{\prime +}$ in the case of 180 ${\unicode{x03BC}} \textrm{m}$ bubbles: (a–c) $y^+=1.5$ ; (d–f) $y^+=5$ and (g–i) $y^+=20$ . (a, d, g) Lift force $F_{y,L}$ ; (b, e, h) bubble slip velocity $v^+_s$ and (c, f, i) bubble velocity $v^{\prime +}_b$ . Dashed lines represent the mean value of each item.

In $y^+\lt 2$ , the mean wall-lift force is greater than the mean lift force (figure 5 b), and the overall trend of bubbles is away from the wall. However, only those in higher speed regions experience sufficient lift to balance the wall-lift force. Consequently, as bubbles close to the wall, they are more likely to accumulate in higher speed regions (figure 7 b) and are more likely to move towards the wall by following stronger sweep events (figure 8 b).

This wall-vicinity preferential mechanism is most pronounced in the viscous sublayer and diminishes with increasing distance from the wall. There are two main reasons for this. First, as the wall distance increases, the flow becomes more dominated by larger-scale turbulent structures, causing the correlation between $\partial u^+/\partial y^+$ and $u^+$ to be rapidly weakened. This, in turn, reduces the correlation between $F_{y,L}^+$ and $u_{@b}^{\prime +}$ . As shown in figure 10(d), at $y^+=5$ , there is still a noticeable correlation between $F_{y,L}^+$ and $u_{@b}^{\prime +}$ , though weaker compared with the correlation at $y^+=1.5$ (figure 10 a). At $y^+=20$ , this correlation nearly disappears (figure 10 g). Second, $F_{y,L}^+$ modulates $v^+_s$ , which in turn affects the overall wall-normal velocity $v^+_b$ , while the magnitude of $v^+_s$ remains at the order of $O(0.01)$ throughout the viscous sublayer and buffer layer (figure 10 b, e, h). In the vicinity of the wall, the wall-normal fluid velocity $v^+$ is small and its magnitude is comparable to the slip velocity, resulting in a strong correlation between $v^+_b$ and $u_{@b}^{\prime +}$ . However, as the wall-normal velocity increases rapidly with wall distance ( $v_{rms}\propto y^2$ ) (Kim, Moin & Moser Reference Kim, Moin and Moser1987), by $y^+=5$ , the magnitude of $v^+_s$ becomes much smaller than that of $v^+_{@b}$ , weakening the correlation between $v^+_b$ and $u_{@b}^{\prime +}$ . For bubbles in the buffer layer, the effect of $v^+_s$ on $v^+_b$ becomes negligible (figure 10 i).

3.1.4. Preferential concentration in flow topologies

We further analysed the preferential concentration of bubbles in the flow topologies, using the local topology classification scheme by Chong, Perry & Cantwell (Reference Chong, Perry and Cantwell1990). The three eigenvalues $\lambda _k(k=1,2,3)$ of the velocity gradient tensor $\partial u_i/\partial x_j$ satisfy the characteristic equation $\lambda ^3+P\lambda ^2+Q\lambda +R=0$ , where $P$ , $Q$ and $R$ are the tensor invariants:

(3.2) \begin{align} &P=-\left (\lambda _1+\lambda _2+\lambda _3\right ) \nonumber ,\\ &Q=\lambda _1\lambda _2+\lambda _1\lambda _3+\lambda _2\lambda _3, \nonumber \\ &R=-\lambda _1\lambda _2\lambda _3. \end{align}

In incompressible flow, $P=0$ , allowing the eigenvalues to be mapped onto the $Q{-}R$ plane. The discriminant $D=(27/4)R^2+Q^3$ determines the nature of the eigenvalues: for $D\gt 0$ , there is one real eigenvalue and two complex conjugates, corresponding to focal topologies (high-enstrophy and low-dissipation). For $D\lt 0$ , all eigenvalues are real, corresponding to nodal topologies (low-enstrophy and high-dissipation) (Hasslberger et al. Reference Hasslberger, Cifani, Chakraborty and Klein2020). As shown in figure 11(a), the curves $D=0$ and $R=0$ (solid lines) divide the $Q{-}R$ plane into four topological quadrants (I–IV), proceeding counter-clockwise from the upper right: (I) unstable focus/ compressing; (II) stable focus/ stretching; (III) stable node/saddle/saddle and (IV) unstable node/saddle/saddle.

Figure 11. Joint p.d.f. of $Q$ and $R$ in the case of 180 ${\unicode{x03BC}} \textrm{m}$ bubbles: (a–c) $y^+=1.5$ ; (d–f) $y^+=5$ and (g–i) $y^+=20$ . (a, d, g) Fluid at grid point; (b, c, e, f, h, i) fluid at bubble location. Contours, all bubbles; solid lines, outward bubbles ( $v_b\gt 0.1v_{b,rms}$ ); dashed lines, inward bubbles ( $v_b\lt -0.1v_{b,rms}$ ).

Figure 11(a, d, g) shows the joint p.d.f. of $(Q^+,R^+)$ for the flow at grid points of the 180 ${\unicode{x03BC}} \textrm{m}$ bubble case. In the viscous sublayer, the flow is more evenly distributed across the four quadrants (figure 11 a,d), As the height increases, in the buffer layer (figure 11 g), the proportion of topological quadrants II and IV increases, and the distribution elongates along the curve $D=0$ , forming a well-known teardrop shape (Blackburn, Mansour & Cantwell Reference Blackburn, Mansour and Cantwell1996). Figure 11(b,c,e,f,h,i) shows joint p.d.f.s of $(Q^+,R^+)$ at bubble locations $(Q^+_{@b},R^+_{@b})$ , with the contours representing the whole bubbles. We also examined the outward and inward bubbles separately, excluding those with near-zero wall-normal velocity, using $0.1v_{b,rms}^{\prime +}$ as the threshold. The solid and dashed isolines represent the outward ( $v_b\gt -0.1v_{b,rms}$ ) and inward ( $v_b\lt -0.1v_{b,rms}$ ) bubbles, respectively.

In the vicinity of the wall ( $y^+=1.5$ ), the contours in figure 11(b, c) show that nearly all bubbles are located in the convergence region ( $D^+\lt 0$ ) and tend to accumulate in the topological quadrant IV, stretching along the $D^+=0$ curve. This behaviour is consistent with the results of Bijlard et al. (Reference Bijlard, Oliemans, Portela and Ooms2010) and Hasslberger et al. (Reference Hasslberger, Cifani, Chakraborty and Klein2020), but with a much stronger preference. This is because, at this $y^+$ , the 180 ${\unicode{x03BC}} \textrm{m}$ bubbles prefer to cluster in strong sweep events, which belong to the topological quadrant IV (unstable node/saddle/saddle) (Chacin & Cantwell Reference Chacin and Cantwell2000). The inward bubbles are more likely to be observed in stronger sweep events and, thus, they are more concentrated in the more active part of topological quadrant IV with larger $Q^+$ and $R^+$ (dashed lines in figure 11 c). However, outward bubbles are still within the sweep event, while it is relatively weaker and thus corresponds to the more stable part of topological quadrant IV (solid lines in figure 11 b).

At a larger distance from the wall ( $y^+=5$ ), as shown in the contours of figure 11(e, f), bubbles predominantly cluster in the topological quadrants III and IV. At this height, bubbles prefer the low-speed streaks generated by the ejection events, which correspond to nodal topologies ( $D^+\lt 0$ ), but do not possess a clear classification in the $Q{-}R$ plane (Picciotto, Marchioli & Soldati Reference Picciotto, Marchioli and Soldati2005). The outward bubbles, which are more numerous, are primarily located in stronger ejection events and tend to concentrate in the more stable topological quadrant III (figure 11 e). Meanwhile, inward bubbles are mainly found in regions where $u^\prime \lt 0$ and $v^\prime \lt 0$ , with a small proportion in sweep events, and these bubbles tend to be concentrated in topological quadrant IV. Similar results have also been observed in particle-laden flows by Bijlard et al. (Reference Bijlard, Oliemans, Portela and Ooms2010).

As the height further increases, as shown in the contours of figure 11(h, i), in the middle of the buffer layer at $y^+=20$ , the distribution of flow at bubble locations in the $Q{-}R$ plane becomes almost identical to that at grid points, similar to the results of Rouson & Eaton (Reference Rouson and Eaton2001). The outward bubbles show a weak preference for ejection events, leading to an increased proportion of the topological quadrant IV, while the distribution of inward bubbles appears slightly stretched, with both the focal and nodal topologies increasing slightly in proportion.

Figure 12 quantitatively illustrates the preferential distribution of bubbles in the four topological quadrants of the $Q{-}R$ plane, comparing the proportions of the flow at grid points (solid lines) with those at bubble locations (dashed lines). As shown in figure 12(a), the 120 ${\unicode{x03BC}} \textrm{m}$ bubbles tend to accumulate in nodal topologies ( $D^+\lt 0$ ) throughout both the viscous sublayer and buffer layer, corresponding to their preference for ejection events, as depicted in figure 8. In figure 12(b), the 180 ${\unicode{x03BC}} \textrm{m}$ bubbles in the vicinity of the wall ( $y^+\lt 3$ ) primarily concentrate in sweep events, corresponding to the topological quadrant IV (unstable node/saddle/saddle). While at $3\lt y^+\lt 30$ , they also show a preference for ejection events, similar to the 120 ${\unicode{x03BC}} \textrm{m}$ bubbles, exhibiting a preference for nodal topologies. In summary, within the viscous sublayer and buffer layer ( $y^+\lt 30$ ), bubbles of all sizes show a strong preference for the nodal topologies, which are located in topological quadrants III and IV. This indicates that due to the force in the wall-normal direction, bubbles prefer regions of low-enstrophy and high-dissipation. This is consistent with the result of light particles, which tend to avoid high-enstrophy regions and, instead, concentrate in nodal topologies (Rouson & Eaton Reference Rouson and Eaton2001; Picciotto et al. Reference Picciotto, Marchioli and Soldati2005).

Figure 12. Fraction of the four quadrants in the $Q{-}R$ plane: (a) 120 ${\unicode{x03BC}} \textrm{m}$ bubbles; (b) 180 ${\unicode{x03BC}} \textrm{m}$ bubbles. The same notation is used in figure 8.

Furthermore, we compare the quadrant proportions in the $u^{\prime }{-}v^{\prime }$ plane (figure 8) and $Q{-}R$ plane (figure 12) between the flow field in bubble-laden cases (solid lines) and unladen case (dotted lines) to analyse the modulation of bubbles on the turbulent structures. Despite the strong preferential concentration behaviours exhibited by the bubbles, the impact of bubbles on the flow is found to be marginal, especially for large bubbles. In the viscous sublayer ( $y^+\lt 5$ ), the introduction of the 120 ${\unicode{x03BC}} \textrm{m}$ bubbles slightly reduces the proportion of Q2 and Q4 events (figure 8 a), indicating the sweep/ejection events are slightly suppressed. This leads to a slight increase in the proportion of nodal topologies (topological quadrants III and IV), as shown in figure 12(a), suggesting a modest increase in dissipation in the viscous sublayer. Above the viscous sublayer, there is virtually no difference between the bubble-laden and unladen flows in both figures 8(a) and 12(a). Similarly, in figures 8(b) and 12(b), the 180 ${\unicode{x03BC}} \textrm{m}$ bubbles show almost no effect on the flow; the proportions of the 180 ${\unicode{x03BC}} \textrm{m}$ bubbles are almost identical to those in the unladen flow throughout the entire domain. This suggests that while bubbles exhibit strong preferential concentration in coherent structures such as the ejection and sweep events, their overall impact on these coherent structures is quite small. This conclusion is somewhat counterintuitive. In the next section, we will explore the underlying mechanisms of how bubbles interact with the flow field.

3.2.1. Velocity statistics

Figure 13(a) shows the normalised mean streamwise fluid velocity $\langle u\rangle ^+$ , mean streamwise bubble velocity $\langle u_b\rangle ^+$ and mean streamwise fluid velocity at bubble location $\langle u_{@b}\rangle ^+$ , compared with those of the unladen flow. All terms are normalised by the friction velocity $u_\tau$ . For the smallest bubbles (120 ${\unicode{x03BC}} \textrm{m}$ ), $\langle u\rangle ^+$ is reduced relative to the unladen flow in the buffer layer and logarithmic layer, indicating a lower bulk velocity and, hence, suggesting an increase in drag. In contrast, for the larger bubbles of 180 and 240 ${\unicode{x03BC}} \textrm{m}$ , the $\langle u\rangle ^+$ profiles remain nearly indistinguishable from that of the unladen flow, suggesting that the turbulence modulation is minimal, which aligns with the result of Molin et al. (Reference Molin, Marchioli and Soldati2012).

Figure 13. (a) Mean streamwise velocity profiles. Solid lines, fluid velocity $\langle u\rangle ^+$ ; dashed lines, bubble velocity $\langle u_b\rangle ^+$ ; dotted lines, fluid velocity at bubble location $\langle u_{@b}\rangle ^+$ . (b) Mean streamwise slip velocity. Circle, slip velocity $\langle u_s\rangle ^+$ ; dashed lines, slip velocity estimated from (3.4).

The mean velocity at bubble location $\langle u_{@b}\rangle ^+$ reveals the tendency of the bubble to concentrate in either high- or low-speed regions, which has been discussed in detail in § 3.1. In the upward channel flow, the bubble experiences a positive buoyancy force. Therefore, the bubble mean velocity $\langle u_b\rangle ^+$ always exceeds $\langle u\rangle ^+$ , with larger bubbles having higher $\langle u_b\rangle ^+$ .

The mean streamwise bubble slip velocity is defined as $\langle u_s\rangle ^+=\langle u_b\rangle ^+-\langle u_{@b}\rangle ^+$ . In the streamwise direction, the drag force and buoyancy force in (2.6) are dominant, leading to

(3.3) \begin{align} \langle u_s\rangle ^+\approx \sqrt {\dfrac {4\left (\rho _f-\rho _b\right )g^+d^+}{3C_D\rho _f}}. \end{align}

Considering the contribution of spanwise and wall-normal slip velocity is minor, the mean bubble Reynolds number can be approximated as $ Re_b\approx \langle u_s\rangle ^+d^+$ . Using the drag force model given in (2.7) (Schiller & Naumann Reference Schiller and Naumann1935), the mean streamwise slip velocity can be solved from

(3.4) \begin{align} \langle u_s\rangle ^+ \left (1+0.15\left (\langle u_s\rangle ^+d^+\right )^{0.687}\right )\approx \left (St/Fr\right ), \end{align}

where $St/Fr$ is a prior parameter as listed in table 1 and scales with the bubble diameter as $St/Fr\propto d^{+\,2}$ . As seen in figure 13(b), the mean streamwise slip velocity $\langle u_s\rangle ^+$ remains approximately constant across the channel, i.e. $\langle u_s\rangle ^+ \approx [1.16,1.25,1.34,1.52,2.31,3.60 ]$ for all six cases. These correspond to the bubble Reynolds numbers $ Re_b$ ranging from 0.83 to 5.15. The numerical solutions of mean slip velocity closely align with the analytical prediction from (3.4). For small bubbles $(Re_b\lt 1)$ , the relative motion lies in the Stokesian region, so the finite Reynolds number correction in (2.7) remains close to unity, yielding the scaling $\langle u_s\rangle ^+\propto d^{+\,2}$ . As $Re_b$ increases, this nonlinear correction becomes significant and the scaling is changed to $\langle u_s\rangle ^+\propto d^{+\,1.7}$ .

Figure 14 illustrates the root mean square (r.m.s.) of fluid velocity fluctuations $u_{i,rms}^{\prime +}$ , bubble velocity fluctuations $u_{bi,rms}^{\prime +}$ and fluid velocity fluctuations at bubble location $u_{@bi,rms}^{\prime +}$ , as well as the Reynolds shear stress. For the larger bubbles (180 and 240 ${\unicode{x03BC}} \textrm{m}$ ), both $u_{i,rms}^{\prime +}$ and $-\langle u^\prime v^\prime \rangle ^+$ almost entirely overlap with those of unladen flow, which aligns with the results in figure 13. This indicates that, in these cases, the bubble has minimal influence on the turbulent structure, as noted by Molin et al. (Reference Molin, Marchioli and Soldati2012). In contrast, the smaller bubbles (120 ${\unicode{x03BC}} \textrm{m}$ ) lead to a reduction in both the turbulence intensity and Reynolds shear stress, suggesting that these bubbles suppress the overall flow field, decreasing both the bulk velocity and the turbulence intensity.

Figure 14. R.m.s. of fluid velocity fluctuations $u_{i,rms}^{\prime +}$ , bubble velocity fluctuations $u_{bi,rms}^{\prime +}$ and fluid velocity fluctuation at bubble location $u_{@bi,rms}^{\prime +}$ : (a) streamwise; (b) wall-normal; (c) spanwise components and (d) Reynolds shear stress.

Owing to their small Stokes number, bubbles in all cases exhibit small slip velocity fluctuations. Consequently, the r.m.s. of bubble velocity fluctuations ( $u_{bi,rms}^{\prime +}$ ) closely match with the r.m.s. of fluid velocity fluctuations at bubble locations ( $u_{@bi,rms}^{\prime +}$ ), as seen in figure 14, indicating that $u_{bi,rms}^{\prime +}$ is primarily governed by the local flow regions where the bubbles preferentially accumulate.

Below the concentration peak, the 180 and 240 ${\unicode{x03BC}} \textrm{m}$ bubbles show strong preference for high-speed regions, resulting in reduced values of $u_{@b,rms}^{\prime +}$ and $u_{b,rms}^{\prime +}$ , as shown in figure 14(a). This behaviour aligns with the joint p.d.f. of $(u_{@b}^{\prime +},v_{@b}^{\prime +})$ in figure 9(a, b), where the $u_{@b}^{\prime +}$ exhibits a narrower distribution than $u^\prime$ due to its preference of high-speed region. Similarly, above the concentration peak but still within the viscous sublayer $(y^+\lt 5)$ , the preference shifts towards the low-speed region. However, $u_{@b,rms}^{\prime +}$ are still lower, which is also consistent with figure 9(c). These suggest that in the near-wall region, due to the strong preferential concentration, regardless of whether high- or low-speed region, $u_{@b,rms}^{\prime +}$ always remains lower than $u_{\textit{rms}}^{\prime +}$ .

In the spanwise direction, as seen in figure 14(c), the extremely high values of $w_{b,rms}^{\prime +}$ and $w_{@b,rms}^{\prime +}$ are observed for the 180 and 240 ${\unicode{x03BC}} \textrm{m}$ bubbles in the vicinity of the wall ( $y^+\lt 2$ ). Almost the same level of enhancement of $w_{b,rms}^{\prime +}$ in the vicinity of the wall is observed in the one-way coupling simulation (not shown here), indicating that it is not caused by the bubble–fluid feedback force. This phenomenon occurs as bubbles are entrained by the sweep events, which transport them towards the wall. These sweep events, located on the side of quasi-streamwise vortices, also carry strong spanwise momentum. In other words, as bubbles are driven towards the wall by the quasi-streamwise vortices, both the wall-normal and spanwise velocity fluctuations of the bubbles, $v_{b,rms}^{\prime +}$ and $w_{b,rms}^{\prime +}$ , are significantly amplified.

As $y^+$ increase into the buffer layer, both $u_{@b,rms}^{\prime +}$ and $u_{b,rms}^{\prime +}$ become slightly larger than those of the fluid, while $v_{b,rms}^{\prime +}$ and $w_{b,rms}^{\prime +}$ closely match with those of the fluid. This behaviour is consistent with the characteristics of low-inertia particles ( $St\lt 1$ ) (Marchioli, Picciotto & Soldati Reference Marchioli, Picciotto and Soldati2007). It may be attributed to the interactions between bubbles and turbulent structures, as bubbles are driven by the lift force in the wall-normal direction to move perpendicularly to the mean flow, between regions of high- and low-speed regions (Portela, Cota & Oliemans Reference Portela, Cota and Oliemans2002). Beyond the buffer layer, all terms of bubble velocity fluctuations converge with those of the fluid phase, as the preferential concentration diminishes and the effect of the lift force becomes negligible.

The 120 ${\unicode{x03BC}} \textrm{m}$ bubble in the viscous sublayer prefers accumulating in the low-speed streak region, resulting in lower $u_{b,rms}^{\prime +}$ and slightly lower $w_{b,rms}^{\prime +}$ compared with the fluid. As the distance from the wall increases, within the buffer layer ( $5\lt y^+\lt 30$ ), the $u_{b,rms}^{\prime +}$ gradually switches to be higher than $u_{\textit{rms}}^{\prime +}$ . This is first because the tendency of bubble accumulation decreases with $y^+$ , causing $u_{b,rms}^{\prime +}$ to converge towards $u_{\textit{rms}}^{\prime +}$ . Additionally, due to the significant streamwise slip velocity, bubbles exhibit a crossing-trajectory effect, where they rapidly pass through turbulent structures, leading to greater bubble streamwise velocity fluctuations. The crossing-trajectory effect is most pronounced in the buffer layer, as this region has the most active turbulence and moderate fluid velocity. Beyond the buffer layer, all terms of bubble velocity in figure 14 converge with those of the fluid phase since there is nearly no preferential concentration, and the cross-trajectory effect gradually weakens.

3.2.2. Energy transfer

In an upward vertical channel, bubbles experience a buoyancy force in the streamwise direction, which performs positive work on the flow and consequently accelerates the flow. This is the most important feedback effect of bubbles on the flow field. We first explore the effect of bubbles on flow from the perspective of the mean flow. The transport equation of the mean kinetic energy (MKE) $K=\langle u \rangle ^2/2$ in the channel can be expressed as (Pope Reference Pope2000)

(3.5) \begin{align} \frac {{\bar DK}}{{\bar Dt}} =\,&\, \underbrace {\nu \frac {\partial }{{\partial y}}\frac {{\partial K}}{{\partial y}}}_{{T_{\nu ,K}}}+\underbrace {\frac {\partial }{{\partial y}}\left ( { - \langle {u'v'} \rangle \langle u \rangle } \right )}_{{T_{\textit{t,K}}}} + \underbrace {\left ( { - \nu \frac {{\partial \langle u \rangle }}{{\partial y}}\frac {{\partial \langle u \rangle }}{{\partial y}}} \right )}_{{\varepsilon _K}}\nonumber \\ - & \underbrace {\left ( { - \langle {u'v'} \rangle \frac {{\partial \langle u \rangle }}{{\partial y}}} \right )}_{{P_k}}+\underbrace {\left ( { - \frac {1}{\rho }\frac {{{\textrm{d}}P}}{{{\textrm{d}}x}}\langle u \rangle } \right )}_{{P_{P,K}}} + \underbrace {\langle u \rangle \langle f_x \rangle }_{{P_{B,K}}}, \end{align}

where $\bar D/\bar Dt=\partial /\partial t+\langle u_i \rangle \partial /\partial x_i$ is the mean material derivatives, $P_{P,K}$ is the mean production term due to the mean pressure gradient and $P_{B,K}$ is the production term of bubble feedback force acting on flow. The total energy production is $P_{K}=P_{P,K}+P_{B,K}$ . As mentioned before, the mean pressure gradient in the unladen flow is increased to ${\textrm{d}}P/{\textrm{d}}x={\textrm{d}}P_0/{\textrm{d}}x+\alpha _v(\rho _f-\rho _b)g$ , meaning that the total driving force is the same for both bubble-laden flow and unladen flow. Here, $T_{\textit{t,K}}$ and $T_{\nu ,K}$ are the diffusion terms due to the turbulence and viscosity, respectively, which transport $K$ towards the near-wall region; $\varepsilon _K$ is the viscous dissipation term and $P_k$ is the turbulent production term, which represents the transfer of mean kinetic energy $K$ to turbulent kinetic energy $k$ .

Figure 15(a) shows the budget terms in (3.5) for flows with 120, 125 and 180 ${\unicode{x03BC}} \textrm{m}$ bubbles, along with the unladen flow. Figure 15(b) provides a magnified view of the mean production terms $P_{P,K}$ , $P_{B,K}$ and $P_{K}$ . For the 120 ${\unicode{x03BC}} \textrm{m}$ bubbles, all terms are decreased compared with the unladen flow. The viscous-related terms $T_{\nu ,K}$ and $\varepsilon _K$ in the vicinity of the wall gradually approach those of the unladen flow as $y^+$ decreases. We focus primarily on the production terms. In the vicinity of the wall ( $y^+\lt 1$ ), the bubble concentration $c$ is extremely high (figure 4 a), while the local velocity $\langle u \rangle$ is very low. Even though the mean streamwise bubble feedback force $\langle f_x \rangle$ is substantial, their contribution $P_{B,K}$ remains limited (figure 15 b). However, since nearly half the bubbles are concentrated within the very near-wall region ( $y^+\lt 10$ ), the bubble volume fraction decreases in the logarithmic layer and outer layer, where $\langle u \rangle$ is high. This results in lower work done on the flow, thus reducing the total production term $P_{K}$ . This reduction in $P_{K}$ further leads to the reduction in both the mean velocity profile (figure 13) and the velocity fluctuation profiles (figure 14) for the 120 ${\unicode{x03BC}} \textrm{m}$ case. The decreased mean velocity further reduces the work done by the mean pressure gradient $P_{P,K}$ .

Figure 15. Budget of the mean kinetic energy $K$ transportation equation: (a) all budget terms; (b) production terms only. Solid lines, unladen flow; dashed lines, 120 ${\unicode{x03BC}} \textrm{m}$ bubbles; dash-dotted lines, $125\,{\unicode{x03BC}} \textrm{m}$ bubbles; dotted lines, 180 ${\unicode{x03BC}} \textrm{m}$ bubbles.

For the $125\,{\unicode{x03BC}} \textrm{m}$ bubbles, although the diameter is only slightly increased, the peak concentration near the wall significantly decreases, resulting in a reduction of $P_{B,K}$ in the vicinity of the wall. However, the increase in $c$ in the outer region leads to an increase in $P_{B,K}$ , resulting in the total production term $P_{K}$ being closer to that of the unladen flow. For the larger 180 ${\unicode{x03BC}} \textrm{m}$ bubbles, although a concentration peak is observed in figure 4(a), it is relatively low and the concentration $c$ at $y^+\gt 10$ is almost equal to the bulk concentration $c_0$ . Consequently, as shown in figure 15(b), $P_{K}$ shows only a very slight increase near the peak of $c$ at $y^+\approx 4$ , with a slight decrease for $y^+\lt 2$ . Beyond the near-wall region, it nearly matches the unladen flow. The terms of the $K$ transportation equation in figure 15(a) are also closely aligned with those of the unladen flow. In general, the smallest bubbles, which are predominantly concentrated in the vicinity of the wall, contribute less work to the mean flow, resulting in a lower bulk velocity and reduced turbulence statistics. In contrast, the larger bubbles, despite accumulating near the wall as well, exhibit a small peak of concentration. However, the overall concentration distribution is relatively uniform across the channel. As a result, the work done by these bubbles on the mean flow is nearly equivalent to an additional mean pressure gradient. This finding is consistent with the previous studies (Molin et al. Reference Molin, Marchioli and Soldati2012; Asiagbe et al. Reference Asiagbe, Fairweather, Njobuenwu and Colombo2019).

We further analyse the effects of bubbles on turbulence by examining the transportation equation of the turbulent kinetic energy (TKE), defined as $k=\langle {u_i^{\prime }u_i'} \rangle /2$ (Lee & Lee Reference Lee and Lee2015; Wang, Xu & Zhao Reference Wang, Xu and Zhao2021):

(3.6) \begin{align} \frac {{\bar Dk}}{{\bar Dt}} = \underbrace { - \langle {u'v'} \rangle \frac {{\partial \langle u \rangle }}{{\partial y}}}_{{P_k}} + \frac {\partial }{{\partial y}}\left \{ {\underbrace { - \frac {{\langle {p'v'} \rangle }}{\rho }}_{{T_{p,k}}} + \underbrace {\nu \frac {{\partial k}}{{\partial y}}}_{{T_{\nu ,k}}}\underbrace { - \langle {kv'} \rangle }_{{T_{t,k}}}} \right \}\underbrace { - \nu \left \langle {\frac {{\partial {u^{\prime}_i}}}{{\partial {x_j}}}\frac {{\partial {u^{\prime}_i}}}{{\partial {x_j}}}}\right \rangle }_{{\varepsilon _k}} + \underbrace {\langle {u^{\prime}_i{f^{\prime}_i}} \rangle }_{{B_{k}}}, \end{align}

where $P_k$ is the turbulent production term that transfers the mean kinetic energy to turbulence (3.5); $T_{p,k}$ , $T_{\nu ,k}$ and $T_{t,k}$ are the diffusion terms due to pressure, viscosity and turbulence, respectively; $\varepsilon _k$ is the dissipation term and $B_{k}$ represents the bubble feedback on the TKE. As shown in figure 16(a), for both the 120 and 125 ${\unicode{x03BC}} \textrm{m}$ bubble cases, all budget terms decrease by nearly the same scale, except the bubble-induced term $B_{k}$ . This is the same as the trends observed in the budget of $K$ in figure 15(a), indicating that bubbles have a suppression effect on turbulence. Near the peak of concentration at $y^+\approx 0.4$ , the bubble-induced term $B_{k}$ reaches the maximum for the 120 ${\unicode{x03BC}} \textrm{m}$ bubbles. This leads to an increase in the viscous-related terms $\varepsilon _k$ and $T_{\nu ,k}$ below the peak and a decrease above the peak. As the bubble diameter increases, all the TKE budget terms increase rapidly. For the large 180 ${\unicode{x03BC}} \textrm{m}$ bubbles, the curves almost collapse with those of the unladen flow.

Figure 16(b) provides a magnified view of the bubble-induced term $B_{k}$ . Due to the small inertia of bubbles, the feedback force on the flow is almost entirely due to the streamwise buoyancy force, while the feedback forces in the wall-normal and spanwise directions are negligible. As a result, $B_{k}=\langle {{u^{\prime}_i}{f^{\prime}_i}}\rangle \approx \langle {{u'}{f^{\prime}_x}}\rangle$ , meaning that the sign of $B_{k}$ depends on the preference of bubbles. For the 120 ${\unicode{x03BC}} \textrm{m}$ bubbles, wall collisions cause them to accumulate in the high-speed regions at the lowest bubble height ( $y^+\approx d^+/2 \approx 0.36$ ), where the concentration also reaches its maximum. This results in a positive shark peak at $y^+\approx 0.3$ . This peak appears to be slightly lower than the lowest bubble height due to the projection method (§ 2.1.3), which distributes the feedback force to the nearby grids. As $y^+$ increases ( $y^+\gt 1$ ), the preference shifts to the low-speed regions, causing $B_{k}$ to become negative. Similarly, for larger bubbles, $B_{k}$ is positive when bubbles cluster in high-speed regions and negative when they concentrate in low-speed regions. Beyond $y^+\gt 20$ , where bubbles are more uniformly distributed, $B_{k}$ approaches zero.

However, compared with other terms in figure 16(a), the $B_{k}$ term remains relatively small. For larger bubbles ( $d\gt 130\,{\unicode{x03BC}} \textrm{m}$ ), $B_{k}$ is approximately two orders of magnitude smaller than the other terms, making its contribution nearly negligible. This is primarily because the regions where bubble preferential concentration is most pronounced are confined to the near-wall viscous sublayer and the bottom of the buffer layer ( $y^+\lt 10$ ), where the velocity is very low. As a result, the total bubble-induced production term $B_{\textit{tot}}=P_{B,K}+B_{k}=\langle {{{u}_i}{{f}_i}}\rangle$ remains very small. The main region of the bubble-induced term is too low to reach the centre of the buffer layer, where the turbulence is most active. It can thus be concluded that, despite the clear preference for bubbles in sweep and ejection events, the overall impact of bubbles on turbulent flow is insignificant, resulting in only slight modulation of the turbulent structures across the channel.

Figure 16. Budget of the TKE transportation equation. (a) All the budget terms, where the same notation of lines is used in figure 15. Solid lines, unladen flow; dashed lines, 120 ${\unicode{x03BC}} \textrm{m}$ bubbles; dash-dotted lines, $125\,{\unicode{x03BC}} \textrm{m}$ bubbles; dotted lines, 180 ${\unicode{x03BC}} \textrm{m}$ bubbles. (b) Bubble-induced term $B_{k}$ .

In addition, the turbulence modulation induced by bubbles can be further elucidated by examining the energy spectra of fluid velocity, as presented in detail in Appendix C. It is interesting to observe an enhancement of the energy spectra at high wavenumbers in the near-wall region where the preferential concentration is most pronounced, while further investigation is still required to elucidate its physical origin.

3.2.3. Skin friction

The mean skin friction coefficient is defined as $C_{\!f} = \tau _w/ ({0.5\rho u_m^2} ) = \mathop {{ Re}}\nolimits _\tau ^2/ (0.5{\mathop { {Re}}\nolimits } _m^2 )$ , where $\tau _w$ is the wall shear stress. As mentioned previously, all the bubble-laden flows and unladen flow share the same friction Reynolds number $\mathop {{ Re}}\nolimits _\tau$ . Therefore, any decrease/increase in the bulk velocity $u_m$ leads to a corresponding increase/decrease in $C_f$ . The drag increase rate is defined as $DI= (C_f-C_{f0} )/C_{f0}$ , where $C_{f0}$ is the skin friction coefficient of the unladen flow. As shown in table 2, drag increase is observed for all the cases. The smallest bubbles (120 ${\unicode{x03BC}} \textrm{m}$ ) exhibit the greatest increase in skin friction, by approximately 16.7 %. As the bubble diameter increases, the $DI$ rapidly decreases. For the large bubbles ( $d\gt 130\,{\unicode{x03BC}} \textrm{m}$ ), the $DI$ is only approximately 1 %, indicating that the skin friction of these cases is almost the same as that of unladen flow.

Table 2. Drag increase rate.

To further analyse the drag increase, we conduct a decomposition of mean skin friction, following the RD identity proposed by Renard & Deck (Reference Renard and Deck2016) and extending this decomposition to the two-phase system by incorporating an additional bubble-induced term. This allows us to better quantify the contributions of bubbles to the whole drag. The RD identity has a physical interpretation based on the budget of mean kinetic energy $K$  (3.5) and can be expressed as

(3.7) \begin{align} {C_f} = \frac {2}{{u_m^3}}\left ( {\underbrace {{{\int _0^h {\nu \left ( {\frac {{\partial \langle u \rangle }}{{\partial y}}} \right )} }^2}{\textrm{d}}y}_{{C_{f,\nu }}} + \underbrace {\int _0^h {\overline { - u'\nu '} \frac {{\partial \langle u \rangle }}{{\partial y}}} {\textrm{d}}y}_{{C_{f,t}}} + \underbrace {\int _0^h {\left ( {{u_m} - \langle u \rangle } \right )\langle {{F_x}} \rangle } {\textrm{d}}y}_{{C_{f,b}}}} \right ), \end{align}

where $C_{f,\nu }$ is the viscous dissipation term, $C_{f,t}$ is the turbulence production term and $C_{f,b}$ is the bubble-induced term. As shown in figure 17, for the large bubbles, $C_{f,b}$ is almost zero, and $C_{f,\nu }$ and $C_{f,t}$ are almost identical to those of the unladen flow. In the case of 120 ${\unicode{x03BC}} \textrm{m}$ bubbles, the increase in $C_{\!f}$ is almost entirely due to the bubble-induced term $C_{f,b}$ , while $C_{f,\nu }$ only shows a slight increase and $C_{f,t}$ remains almost unchanged. This suggests that, although all velocity statistics decrease in the 120 ${\unicode{x03BC}} \textrm{m}$ case, the rate of decrease aligns with the reduction in bulk velocity. The strong wall accumulation of 120 ${\unicode{x03BC}} \textrm{m}$ bubbles at $y^+\lt 1$ generates a significant positive feedback force in the vicinity of the wall. However, this effect is confined to the extremely near-wall region ( $y^+\lt 1$ ) and only leads to an increase in the wall shear, while flow structures in the main region of the flow ( $y^+\gt 1$ ) are only slightly affected. This further confirms that smaller bubbles tend to accumulate near the wall, contributing less work to the flow, which results in lower bulk velocity and reduced turbulent intensity. However, larger bubbles maintain a more uniform concentration across the entire channel and their effect on the flow is almost equivalent to adding an additional mean pressure gradient.

Figure 17. Skin friction coefficient and its decomposition using the RD identity (3.7). Dashed lines, unladen flow.