3.1. Evolution of bubbles

Figure 2.

Instantaneous snapshots of bubbly upflow in a vertical channel. The density ratio is $\rho _l/\rho _g=1000$ , 30 and 10 for the top row, middle row and bottom row, respectively. The viscousity ratio is $\mu _l / \mu _g = 100$ for all cases.

We start the results analysis from the evolution of bubbles. Figure 2 presents instantaneous snapshots of bubbles for three cases. The top row shows the results for case R1000 with $\rho _l/\rho _g = 1000$ , which approximates a water–air system. Deformable bubbles aggregate in the core region, an observation qualitatively consistent with several previous studies at lower density ratios (Lu & Tryggvason Reference Lu and Tryggvason2008; Bolotnov et al. Reference Bolotnov, Jansen, Drew, Oberai, Lahey and Podowski2011; Dabiri et al. Reference Dabiri, Lu and Tryggvason2013; Du Cluzeau et al. Reference du Cluzeau, Bois and Toutant2019). Additionally, a large bubble forms in the core region, moving faster than the surrounding smaller bubbles. This large bubble continuously absorbs small bubbles and persists throughout the simulation duration. Similar bubble behaviour is also observed in the middle row for case R30 with $\rho _l/\rho _g = 30$ . However, such a large bubble is not seen in the bottom row, which shows the results for case R10 with the lowest density ratio $\rho _l/\rho _g = 10$ under investigation. Instead, the size variation among bubbles is less significant in case R10 than in cases R30 and R1000, and the number of bubbles is smaller in case R10 than in cases R30 and R1000. The bubble distribution in case R10 is similar to the results reported by Lu & Tryggvason (Reference Lu and Tryggvason2018). A more detailed evolution process of bubbles is provided in the ementary movies.

Figure 3.

Evolution of bubble size for cases with density ratio (a) $\rho _l/\rho _g=10$ and (b) $\rho _l/\rho _g=30$ . The contours represent the void fraction corresponding to the effective radius $\alpha _r$ (scales on the left vertical axis), normalised by the total void fraction of gas $\alpha _t$ . The superimposed blue solid line denotes the total number of bubbles $N_b$ (scales on the right vertical axis).

To visualise the evolution of bubbles with different sizes, we define an effective radius $r_e$ for each bubble based on its volume $V_b$ as $r_e = (0.75V_b/\pi )^{1/3}$ . The volume of each bubble is determined using a connected component labelling algorithm (Herrmann Reference Herrmann2010). Figure 3 illustrates the evolution of bubble sizes for case R10N8 with $\rho _l/\rho _g = 10$ and case R30N8 with $\rho _l/\rho _g = 30$ . The contours represent the void fraction corresponding to the effective radius $\alpha _r$ , normalised by the total void fraction of gas $\alpha _t$ . The evolution of the total number of bubbles, $N_b$ , is superimposed as the solid line to demonstrate the convergence of the bubble number.

The evolution of the total number of bubbles shown by the solid line in figure 3 converges at $t/t^* \approx 50$ . However, a longer time is needed for the bubbles with different sizes to evolve into a steady state. In these two cases, eight bubbles of equal diameter are initialised in the channel, and there is no change in their size during the early stages. After these initial bubbles break up into smaller ones around $t/t^* \approx 12$ , bubbles with different sizes appear. For the low-density ratio case $\rho _l/\rho _g = 10$ (figure 3 a), bubbles first break into smaller ones with effective radii mainly ranging from $r_e/H = 0.1$ to $0.3$ during $t/t^* = 12{-}30$ . Subsequently, smaller bubbles tend to merge into larger ones during $t/t^* = 30{-}100$ , with the largest bubble reaching a size of $r_e/H \approx 0.55$ . The largest bubbles break into smaller ones during $t/t^* = 100{-}125$ , and after $t/t^* = 125$ , the bubble size evolves to the statistically stationary state, with bubbles appearing across all scales from $r_e/H = 0.02$ to $0.5$ . For the case with $\rho _l/\rho _g = 30$ (figure 3 b), the initial eight bubbles also start to break up into smaller ones after $t/t^* \approx 12$ . Subsequently, the size of the largest bubble continuously increases during $t/t^{\ast} = 20 - 70$ . The size of the largest bubble reaches $r_e/H \approx 0.65$ , and it persists for the entire simulation duration. By contrasting figures 3(a) and 3(b), it is seen that the number of smaller bubbles (with sizes ranging from $r_e/H = 0.02$ to $0.5$ ) for $\rho _l/\rho _g = 30$ is smaller than that for $\rho _l/\rho _g = 10$ . A scale separation is evident in figure 3(b), where bubbles with $ r_e/H \approx 0.4$ almost disappear.

Figure 4.

Bubble size spectra for cases with different density ratios. The two black solid lines indicate the $-3/2$ and $-10/3$ power law proposed by Garrett et al. (Reference Garrett, Li and Farmer2000).

To further quantify the number of bubbles with different sizes, figure 4 shows the bubble size spectra for different cases after the bubbles evolve to a statistically stationary state. We note here that although the bubble size is sensitive to the initial condition at an early stage, the bubble statistics and turbulent statistics in the stationary steady-state are independent of the initial conditions (see Appendix B). Therefore, in the following content, the results for cases R1000, R100, R30 and R10 are analysed, while the results for cases R30N8 and R10N8 are not presented. It is seen from figure 4 that the bubbles at small scale show similar behaviour, while the number of large-scale bubbles is more sensitive to the density ratio. As the density ratio increases, the size of the largest bubble increases, whereas the number of bubbles with intermediate sizes between $0.07H$ and $0.2H$ decreases. This decrease in the number of intermediate-sized bubbles occurs because these bubbles in the channel core region are absorbed by the largest bubble, as shown in figure 3(b) and the supplementary movies.

The bubble size spectra are also compared with the power-law scaling proposed by Garrett et al. (Reference Garrett, Li and Farmer2000) for bubbles in the upper ocean in figure 4. A similar scaling law was also observed in a droplet-laden homogeneous isotropic turbulence (Cannon et al. Reference Cannon, Soligo and Rosti2024) and in a bubbly channel flow without gravity (Soligo et al. Reference Soligo, Roccon and Soldati2019). The result for $\rho _l / \rho _g = 10$ is consistent with the empirical law, but it deviates as the density ratio increases due to the combined effect of solid wall and gravity. In contrast to an upper ocean or a homogeneous condition, under which bubbles are free to cascade without the constraints of solid boundaries, bubbles are confined between two solid walls in the present study. The combined effect of gravity drives the bubbles to aggregate towards the channel core (Ervin & Tryggvason Reference Ervin and Tryggvason1997), ultimately resulting in the coalescence of smaller bubbles into larger ones. As the density ratio increases, the difference in the gravitational acceleration between gas and liquid increases, and as a result, the deviation of the bubble size spectra from the $-10/3$ law becomes more significant. More discussions about the effect of gravity on the bubble size are given in § 3.3 through a slip velocity between two phases (see figure 10 and corresponding discussions).

3.2. Turbulent statistics

Figure 5.

Profiles of mean momentum $\langle \rho u \rangle$ in (a) wall units and (b) outer-layer units for single-phase turbulence and bubbly turbulence with different density ratios. The single-phase flow result of Tsukahara et al. (Reference Tsukahara, Seki, Kawamura and Tochio2005) is superimposed for validation.

Figure 5 compares the profiles of mean momentum $\left \langle \rho u \right \rangle$ for different density ratios. The angular brackets $\langle \ \rangle$ denote time and plane averaging as

(3.1) \begin{align} \langle f \rangle (y) = \frac {1}{L_x L_z \Delta T} \int _0^{\Delta T} \int _{0}^{L_z} \int _{0}^{L_x} f(x, y, z, t) \mathrm {d}x \mathrm {d}z \mathrm {d}t , \end{align}

where $f$ represents an arbitrary variable and $\Delta T$ is the time duration of the data used for performing time averaging. Although the definition of averaging remains mathematically consistent with single-phase turbulent channel flow, what is analysed here is the mean momentum $\langle \rho u \rangle$ . With the density involved in the averaging process, the mean momentum can be also regarded as a density-weighted mean velocity. This is different from single-phase turbulence. In two-phase flows, the momentum is a conservative quantity while the velocity is not. Therefore, comparing the mean momentum for different cases is more rational.

In figure 5(a), the results are shown in wall units $y^+ = \rho _b y u_\tau / \mu _b$ , with the mean momentum non-dimensionalised using $\rho _l u_\tau$ . The result for single-phase flow is compared with that of Tsukahara et al. (Reference Tsukahara, Seki, Kawamura and Tochio2005) for validation. The single-phase mean momentum agrees with the result in the literature. The mean momentum for single-phase flow aligns with the linear law in the viscous sub-layer and the logarithmic law in the outer layer. In bubbly flows, although the linear law still holds in the viscous sub-layer, the outer layer no longer conforms to the logarithmic law. As the density ratio increases from 10 to 30, the magnitude of $\left \langle \rho u \right \rangle / \rho _l u_\tau$ continues to decrease in the core region. However, further increasing the density ratio from 30 to 1000 does not lead to significant changes in $\left \langle \rho u \right \rangle / \rho _l u_\tau$ . The effect of density ratio on the mean momentum can be equivalently shown using an outer-layer scaling (figure 5 b), with the bulk mean velocity $u_b=M_b/\rho _b$ and half-width of channel $H$ being the characteristic velocity and length scales, respectively. The gradient of mean momentum is steeper in bubbly flows in the near-wall region. The velocity gradient further increases as the density ratio increases from $\rho _l / \rho _g = 10$ to 30. The profiles of $\left \langle \rho u \right \rangle / \rho _l u_b$ for $\rho _l / \rho _g = 30$ , 100 and 1000 collapse, indicating a saturated effect of density ratio on the mean momentum for higher density ratios. We also examined the effect of viscosity ratio on the turbulent statistics. The viscosity ratio effect is less significant than the density ratio effect (see Appendix C). Therefore, the following analyses focus on the effects of the density ratio.

To further explore the momentum transfer in bubbly turbulence, we integrated the mean momentum equation from the bottom wall to an arbitrary wall-normal coordinate $y$ to yield the following balance equation of the mean shear stress:

(3.2) \begin{align} \Pi y = \underbrace {\left \langle \mu \left (\frac {\partial u}{\partial y } + \frac {\partial v}{\partial x}\right )\right \rangle }_{{S}} \underbrace {-\left \langle \rho u'v' \right \rangle }_{{C}} \underbrace {-\int _{-H}^{y}{\left \langle \rho -{{\rho }_{b}} \right \rangle g{\textrm{d}}y}}_{ {G}}+\underbrace {\int _{-H}^{y}{\left \langle {\sigma \kappa \frac {\partial H(\phi )}{\partial {x}}} \right \rangle {\textrm{d}}y}}_{ {T}} . \end{align}

On the right-hand side, S, C, G and T represent the shear stresses corresponding to the viscosity, convection, gravity and surface tension, respectively. The convection shear stress C is also known as Reynolds shear stress in single-phase turbulence. On the left-hand side, the total shear stress $\Pi y$ shows a straight line with a slope of $\Pi$ . Figure 6 compares the profiles of the mean stresses between single-phase turbulence and bubbly turbulence with $\rho _l / \rho _g = 1000$ . In the single-phase flow (figure 6 a), the total shear stress is dominated by the viscosity shear stress $S$ and Reynolds stress $C$ in the near-wall and core regions, respectively. In the bubbly flow, as shown in figure 6(b), the viscous stress still dominates the near-wall region. In the core region, the gravity term $G$ exhibits an opposite sign to the total shear stress $\Pi y$ . Compared with the single-phase turbulence, the magnitude of the convection shear stress $C$ is larger in bubbly turbulence, with its peak shifting towards the core region, to balance the negative effect of $G$ . The surface tension term $T$ provides shear stress when bubbles are stretched in the shear direction (Lu & Tryggvason Reference Lu and Tryggvason2019; Du Cluzeau et al. Reference du Cluzeau, Bois and Toutant2019; Cannon et al. Reference Cannon, Izbassarov, Tammisola, Brandt and Rosti2021; Lee et al. Reference Lee, Chang, Kim and Choi2024). In the present cases, it is small compared with other terms.

Figure 6.

Balance of mean shear stress given by (3.2) for (a) single-phase turbulence and (b) bubbly turbulence with $\rho _l/\rho _g = 1000$ . All terms are non-dimensionalised using $\rho _bu^2_b$ .

In single-phase turbulence, the well-known FIK identity (Fukagata, Iwamoto & Kasagi Reference Fukagata, Iwamoto and Kasagi2002) links the wall friction to the turbulent motion. To investigate the contributions of the different components of the mean shear stress in (3.2) to the wall friction in bubbly turbulence, we derive the following generalised form of the FIK identity in a bubbly channel flow:

(3.3) \begin{align} \int _{-H}^{0}{yC^{{eff}}{\textrm{d}}y}+\int _{-H}^{0}{y\left \langle \mu \left (\frac {\partial u}{\partial y } + \frac {\partial v}{\partial x}\right )\right \rangle {\textrm{d}}y}=\int _{-H}^{0}{\Pi y^2{\textrm{d}}y}=\frac {\tau _wH^2}{3}, \end{align}

where $C^{{eff}}=G+C+T$ denotes a defined effective Reynolds stress to maintain the form consistent with the original FIK identity, and $\tau _w = \Pi H$ denotes the wall shear stress. The skin friction coefficient $C_f$ is then decomposed as follows:

(3.4) \begin{align} C_f=\frac {2\tau _w}{\rho u^2_b}=\underbrace {\frac {6}{\rho u^2_bH^2}\int _{-H}^{0}{yC^{{eff}}{\textrm{d}}y}}_{C^T_f}+\underbrace {\frac {6}{\rho u^2_bH^2}\int _{-H}^{0}{y\left \langle \mu \left (\frac {\partial u}{\partial y } + \frac {\partial v}{\partial x}\right )\right \rangle {\textrm{d}}y}}_{C^L_f}, \end{align}

where $C^T_f$ and $C^L_f$ represent the turbulent and laminar contributions, respectively, to the total wall friction coefficient $C_f$ . For single-phase turbulent channel flows, the laminar part $C^L_f=6\mu /\rho u_bH=6/Re_b$ , with $Re_b=\rho _bu_bH/\mu _b$ , can be obtained analytically. In bubbly flows, analytical integration is not feasible since $\mu$ varies with local flow conditions. Consequently, we performed numerical integration of this term using simulation data. Figure 7 illustrates the laminar and turbulent contributions to the total wall friction coefficient. For comparison, the decomposed wall friction coefficients $C^T_f$ and $C^L_f$ were normalised using the analytical laminar value, $6/Re_b$ . The presence of bubbles does not alter the laminar part of $C^L_f$ , whereas the turbulent part of $C^T_f$ is higher in bubbly flows than in single-phase flows. The friction Reynolds number $Re_\tau$ for all cases is superimposed in figure 7, and it can be observed that $Re_\tau$ follows the behaviour of $C^T_f$ as the density ratio increases. Specifically, both $Re_\tau$ and $C^T_f$ were larger in bubbly turbulence than in single-phase turbulence, whereas as the density ratio increased from 10 to 30, both increased further. The increase in $Re_\tau$ and $C^T_f$ slowed down as the density ratio further increased from 30 to 1000, which is consistent with the observations in figure 5 for the saturation of the mean momentum. The above observation suggests that the inclusion of bubbles induces additional drag compared with a single-phase flow.

Figure 7.

Contributions of laminar ( $C^L_f$ , triangles) and turbulent ( $C^T_f$ , circles) parts to the total wall friction coefficient ( $C_f$ ) in cases with different density ratios (scales on the left $y$ -axis). Both $C^L_f$ and $C^T_f$ are normalised using the laminar contribution $6/Re_b$ for single-phase flow. The diamond shows the friction Reynolds number in different cases (scales on the right $y$ -axis). Solid and hollow symbols represent single-phase and bubble flows, respectively.

Figure 8.

Profiles of (a) Reynolds stress $C$ , (b) effective Reynolds stress $C^{{eff}}\!,$ (c) bubble internal force $B$ and (d) averaged bubble void fraction $\langle \alpha \rangle$ for single-phase turbulence and bubbly turbulence with different density ratios.

The results shown in figure 7 suggest that the effect of density ratio on the wall friction is primarily determined by the effective Reynolds stress $C^{{eff}}$ . This reminds us to decompose the convection term $C$ into two parts in the bubbly turbulence to examine the effect of density ratio separately. One part is $B=-(G+T)$ that balances the internal forces exerted by the bubbles, including the effects of gravity ( $G$ ) and surface tension ( $T$ ). This component has no net contribution to wall friction. The other part is the effective Reynolds stress $C^{{eff}}$ that influences the wall friction.

Figure 8 compares the profiles of $C$ , $C^{{eff}}$ , $B$ and bubble void fraction $\langle \alpha \rangle$ for various density ratios. As shown in figure 8(a), the magnitude of $C$ increases as the density ratio increases from 10 to 30, while the difference of $C$ between R30 an R1000 is less significant. Figure 8(b) shows that the peak of the effective Reynolds stress $C^{{eff}}$ is located in the near-wall region for all cases. As the density ratio increases, the magnitude of this peak increases, with its location shifting towards the wall. This observation indicates that the bubbles influence the wall friction by altering the near-wall structures. Figure 8(c) shows that $B$ is trivial in a single-phase flow. The significant increase in the magnitude of $C$ with the inclusion of bubbles is primarily caused by the bubble internal force component, $B$ . The peak of $B$ is located at approximately $y/H = -0.5$ , resulting in the shift of the convective shear stress, $C$ , towards the channel centre. The bubble internal forces are also closely related to the mean gas void fraction. Figure 8(d) presents the vertical profiles of the averaged bubble void fraction, $\langle \alpha \rangle$ , for different cases, with the horizontal dash-dotted line representing bulk void fraction $\alpha _t$ . The intersection points between this horizontal line and the profiles of $\langle \alpha \rangle$ is close to the peak locations of $B$ in figure 8(c) for all cases. The greater number of deformable bubbles with sizes $r_e/H = 0.1\sim 0.4$ (figure 4) in the case with $\rho _l/\rho _g = 10$ results in a larger proportion of gas being distributed in the channel core region ( $y/H \geqslant - 0.75$ ). Consequently, the peak position of $B$ for $\rho _l/\rho _g = 10$ is closer to the core region than that for other cases. Another notable observation is that $\langle \alpha \rangle$ is larger near $y/H = -0.5$ in cases with higher density ratios compared with the lowest density ratio case. This is because the dominant large bubble is substantial in size and extends to the $y/H = -0.5$ region, resulting in a higher void fraction at this location. These observations suggest that while the internal forces do not directly contribute to wall friction, they significantly influence the flow dynamics in the core region.

3.3. Interaction between turbulence and bubbles

Figure 9.

Instantaneous snapshots of the bubbles and vortex structures coloured by $u$ . The vortex structures are shown using the isosurface of Q = 1 and Q = 10 for single-phase and bubbly flow cases, respectively. The upper and lower rows show the vortex structures conditioned by $u\lt u_b$ and $u\gt u_b$ , respectively. Each column shows the results for one case: (a, b) single-phase flow, (c, d) $\rho _l / \rho _g = 10$ , (e, f) $\rho _l / \rho _g = 30$ , (g, h) $\rho _l / \rho _g = 100$ and (i, j) $\rho _l / \rho _g = 1000$ .

From the above analysis, the influence of the density ratio on the flow dynamics differs between the near-wall and core regions. In this section, we further analyse the influences of turbulence and bubbles on each other. Figure 9 compares instantaneous snapshots of the vortex structures and bubbles in different cases. Because the streamwise velocity $u$ is in general larger in the core region than in the near-wall region, its value is used as a condition to distinguish the vortex structures in different regions. Specifically, in the upper and lower rows of figure 9, vortex structures with $u\lt u_b$ and $u\gt u_b$ are shown. Therefore, vortex structures in the near-wall and core regions are observed in the upper and lower rows, respectively.

The upper row in figure 9 shows that small bubbles dominate the near-wall region. The vortex structures are broken into smaller structures by these small bubbles (figures 9 c, 9 e, 9 g and 9 i). This effect enhances the wall-normal fluctuation in the near-wall region, leading to an increase in the wall friction. The distribution of small bubbles in the near-wall region remains similar for cases with $\rho _l/\rho _g \geqslant 30$ , such that the wall friction saturates for higher density ratios.

In the core region, large bubbles tend to generate vortex structures in their wakes (figures 9 d, 9 f, 9 h and 9 j), and the sizes of these large bubbles increases with an increasing density ratio. Furthermore, more vortex structures are observed in the wake of large bubbles in the cases with higher density ratios. Several small bubbles are trapped in the wake regions of the large bubbles. These small bubbles are continuously absorbed by larger bubbles, increasing the size of the large bubble (figures 9 f, 9 h and 9 j).

Figure 10.

Profiles of (a) mean velocity in the gas phase $\langle u \rangle _g$ , (b) mean velocity in the liquid phase $\langle u \rangle _l$ and (c) slip velocity $\langle u \rangle _s = \langle u \rangle _g - \langle u \rangle _l$ for bubbly turbulence with different density ratios.

The vortex structures in the wake region of large bubbles are attributed to the slip velocity between the liquid and gas phases. Figure 10 shows a comparison of the mean velocities of the gas phase $\langle u \rangle _g$ , liquid phase $\langle u \rangle _l$ and the slip velocity $\langle u \rangle _s = \langle u \rangle _g - \langle u \rangle _l$ for various density ratios. As shown in figure 10(a), owing to the stronger buoyancy effect, the mean velocity of the gas phase is higher for the case with a larger density ratio. Meanwhile, the change in liquid velocity is insignificant (figure 10 b), resulting in a higher slip velocity between the two phases in cases with higher density ratios (figure 10 c). The relative motion between bubbles and the surrounding liquid generates vortex structures in the wake of the large bubble, and more vortex structures are observed in the wake of bubbles with an increasing density ratio (figures 9 c, 9 e, 9 g and 9 i) because of the larger slip velocity (figure 10 c). However, as observed in figures 5 and 8, the differences in mean momentum $\langle \rho u \rangle$ and Reynolds shear stress $-\langle \rho u' v'\rangle$ are insignificant for density ratios larger than 30. This further confirms that the wall friction is mainly influenced by the near-wall structures.

Figure 11.

Profiles of (a) $u_{\mathrm {rms}}$ , (b) $v_{\mathrm {rms}}$ and the density-weighted correlation coefficient $R_{uv}$ between streamwise and wall-normal velocity fluctuations for single-phase turbulence and bubbly turbulence with different density ratios. .

This analysis of the slip velocity explains the observation in figure 9 that the vortex structures in the core region are increasingly intensified for a higher density ratio. To quantitatively demonstrate the contribution of wake vortices, we examined the density weighted root-mean-square velocities $u_{\mathrm {rms}} = (\langle \rho u^\prime u^\prime \rangle / \rho _b )^{1/2}$ and $ v_{\mathrm {rms}} = ( \langle \rho v^\prime v^\prime \rangle / \rho _b )^{1/2}$ , and $R_{uv} = -\langle \rho u^\prime v^\prime \rangle / (\rho _b u_{\mathrm {rms}} v_{\mathrm {rms}} )$ , the correlation coefficient between $u'$ and $v'$ . We weighted the correlation coefficient by density because $-\langle \rho u^\prime v^\prime \rangle$ appears as the Reynolds shear stress in the balance equation (3.2) of the mean momentum $\langle \rho u \rangle$ for two-fluid flow rather than $-\langle u^\prime v^\prime \rangle$ in single-phase flow. The results of $u_{\mathrm {rms}}$ , $v_{\mathrm {rms}}$ and $R_{uv}$ for various density ratios are shown in figure 11. From these figures, dual effects of wake vortices on the Reynolds shear stress are observed. On the one hand, the occurrence of wake vortices enhances the velocity fluctuations in both the streamwise and vertical directions. Specifically, the magnitudes of both $u_{\mathrm {rms}}$ and $v_{\mathrm {rms}}$ increase monotonically as the density ratio increases, with larger values occurring in the channel core for higher density ratio cases. This observation suggests that at higher density ratios, the stronger vortex structures (as shown in figure 9) in the wakes of large bubbles are responsible for intensifying velocity fluctuations in the core region. On the other hand, the magnitude of $R_{uv}$ decreases monotonically as the density ratio increases. This indicates that the wake vortices also impose a decorrelation effect on the streamwise and wall-normal velocity fluctuations. The Reynolds stress can be also written as $-\langle \rho u^\prime v^\prime \rangle = \rho _b R_{uv} u_{\mathrm {rms}} v_{\mathrm {rms}}$ . The enhancement effects of wake vortices on $u_{\mathrm {rms}}$ and $v_{\mathrm {rms}}$ are balanced by their suppression effect on $R_{uv}$ . Consequently, the Reynolds stress remains almost unchanged for cases with different density ratios in the core region.

Figure 12.

Profiles of Reynolds stresses $-\langle \rho u'v'\rangle$ conditioned by (a) Q1, (b) Q2, (c) Q3 and (d) Q4 events.

To further investigate the mechanism underlying the saturation of wall friction for $\rho _l / \rho _g \geqslant 30$ , and the decorrelation effect of wake vortices on $u'$ and $v'$ in the core region, we performed a quadrant analysis of the Reynolds stress. Figure 12 depicts the Reynolds stress conditioned by the signs of $u'$ and $v'$ . Four quadrants, Q1–Q4, denote the conditions $(u'\gt 0, v'\gt 0)$ , $(u'\lt 0, v'\gt 0)$ , $(u'\lt 0, v'\lt 0)$ and $(u'\gt 0, v'\lt 0)$ , respectively. It is observed that Reynolds stresses corresponding to all four quadrants are larger in bubbly turbulence than in the single-phase turbulence. This indicates that the wall-normal motion and corresponding momentum convection are more active in a bubbly channel. In the near-wall region ( $y/H \leqslant -0.9$ ), as the density ratio increases from 10 to 30, the increase in Q2 and Q4 events is more significant, leading to the overall increase in $C^{{eff}}$ and leading to a higher wall friction. Additionally, in the near-wall region, the contributions of Q1–Q4 to the total Reynolds stress all saturate for $\rho _l/\rho _g \geqslant 30$ , resulting in the saturation of effective Reynolds stress and wall friction. In the core region ( $y/H \geqslant -0.5$ ), an increase in the contribution of Q1 and Q3 events is comparable to that of Q2 and Q4 events. This results in a negligible net contribution to the total Reynolds stress, as shown in figure 8(a).

It is worth noting that the suppression of turbulence in bubbly upflow mentioned by Takagi & Matsumoto (Reference Takagi and Matsumoto2011) is not observed in the present study. The experimental results of So et al. (Reference So, Morikita, Takagi and Matsumoto2002) indicate that adding a certain amount of surfactant to the liquid can prevent the coalescence of bubbles into larger ones, resulting in bubbles with a diameter of approximately 1 mm occupying the channel. Due to shear-induced lateral forces, bubbles accumulate near the wall, forming a sheet-like structure that pulls the surrounding liquid upward through buoyancy. The buoyant bubble sheet typically spans a few hundred wall units in a turbulent boundary layer, creating two distinct regions with different velocity fluctuation behaviours. In the narrow region between the wall and the bubble sheet, fluctuations are primarily generated by high-shear-rate turbulence. In the outer region, fluctuations are caused by the motion of bubbles rising in a plug-like flow, with significantly suppressed velocity fluctuations compared with single-phase turbulence.

However, a bubble sheet is not observed in the present study, even though the flow parameters for case R1000 with $\rho _l / \rho _g = 1000$ are close to those of the experimental set-up. This discrepancy is presumably because there are no surfactants in the present simulations, as such bubbles coalesce to form larger bubbles. The wakes of larger bubbles interact with the wall layer, preventing smaller bubbles (with $r_e/H = 0.025$ , corresponding to 1-mm bubbles in the experiment) from forming a bubble sheet. Consequently, large buoyancy effects near the wall are not observed. Furthermore, as noted by Takagi & Matsumoto (Reference Takagi and Matsumoto2011), the tangential stress at the gas–liquid interface resulting from the Marangoni effect becomes non-negligible when surfactants are added. Bubble clustering is observed only when a specific amount of surfactants is added; otherwise, the bubbles are uniformly dispersed across the channel. This could also be another important factor that no bubble sheet is formed in the present simulations.

The above discussion on the difference between the present study and the experiments of So et al. (Reference So, Morikita, Takagi and Matsumoto2002) is also supported by several previous numerical studies. By preventing bubble coalescence, Lu & Tryggvason (Reference Lu and Tryggvason2013) observed the occurrence of bubble sheets near the wall and the suppression of velocity fluctuations. However, the suppression disappears when bubbles are allowed to coalesce, as shown in a later study of Lu & Tryggvason (Reference Lu and Tryggvason2019). They also showed that the wakes of larger bubbles in the channel core can disturb the wall layer, even at a relatively low density ratio of 10. Thus, it is reasonable that the bubble sheet is absent in our study.