3.1. Characteristics of the scalar dynamics

Figure 1 shows instantaneous snapshots from case B4 of the full periodic computational domain and the bubble interfaces coloured by the total scalar field $c$ with $\nabla ^v \langle c \rangle$ in figure 1(a) and $\nabla ^h \langle c \rangle$ in figure 1(b). Vertical cross-sectional views of the total scalar field $c$ and the scalar disturbance field $c'$ from figure 1(a) are shown in figure 2 and from figure 1(b) in figure 3.

Figure 1.

Istantaneous snapshot from case B4 with bubble interfaces and periodic domain boundaries coloured by the total scalar field $c$. The volume fraction is $\phi =5.2\,\%, Sc_l = 7$ and $\nabla ^v \langle c \rangle =1$ is imposed in panel (a) and $\nabla ^h \langle c \rangle =1$ in panel (b).

Figure 2.

Vertical cross-sectional view of the total $c$ and disturbance $c'$ scalar fields from case B4 with $\phi =5.2\,\%, \nabla ^v \langle c \rangle =1$ and $Sc_l = 7$. The thin black lines represent the bubble interfaces.

Figure 3.

Vertical cross-sectional view of the total $c$ and disturbance $c'$ scalar fields from case B4 with $\phi =5.2\,\%, \nabla ^h \langle c \rangle =1$ and $Sc_l = 7$.

The generation of the scalar disturbances stems from the last term on the right-hand side of (2.4) that (given a positive $\nabla ^v \langle c \rangle$) predicts large negative disturbances in high upward velocity regions and positive disturbances in downward velocity regions. Consequently, we speculate that the scalar disturbances in the case of $\nabla ^v \langle c \rangle$ (shown in figure 2b) are generated by several interacting mechanisms. The first is the generation by the gas phase that preferentially moves in the vertical direction (generation at scales $k/k_{d_b}\approx 1$). The second is the associated generation by the upward flowing liquid in the bubble wakes that is related to the capture–release mechanism described in Alméras et al. (Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015, Reference Alméras, Risso, Roig, Plais and Augier2018) (generation at scales $k/k_{d_b} \gtrapprox 1$). The third is the generation by the average downward flowing liquid in between the rising bubbles due to continuity (generation at scales $k/k_{d_b} \lessapprox 1$). The last is the production by the bubble-induced turbulence that was found to dominate the dispersion of a low-diffusive dye in a homogeneous bubbly flow similar to the flow considered in the present study (Alméras et al. Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015).

In contrast, the scalar disturbances in the case of $\nabla ^h \langle c \rangle$ (shown in figure 3b) show different characteristics. Here, the scalar disturbances are mainly generated by the bubble-induced turbulence since the average bubble motion in the horizontal directions equals zero. However, since the bubble-induced turbulence is characterised by anisotropy (due to the preferential motion of the bubble in the vertical direction), the influence of the bubble-induced turbulence on the scalar dynamics is not, in general, the same in cases with different mean gradient orientations. Given the different characteristics of the scalar disturbances in the cases of $\nabla ^v \langle c \rangle$ and $\nabla ^h \langle c \rangle$, we thus expect some differences in the scalar statistics that are analysed in § 3.4.

3.2. Bubble clustering

Following the methodology of Tagawa et al. (Reference Tagawa, Roghair, Prakash, van Sint Annaland, Kuipers, Sun and Lohse2013) and Pandey et al. (Reference Pandey, Ramadugu and Perlekar2020), we analyse the bubble clustering in our DNS with volume fractions $\phi =1.7\,\%$ and $\phi =5.2\,\%$. We compute the centre of mass of each bubble in the DNS at a statistically steady state and compute Voronoi tessellations using the Voro++ library (Rycroft Reference Rycroft2009). For each volume fraction, we compute the standard deviation $\varSigma$ of the Voronoi volumes. This $\varSigma$ is compared with the standard deviation $\varSigma _{{rnd}}$ of the Voronoi volumes in 200 configurations of randomly positioned, non-overlapping bubbles in the same cubic domain with length $L$ as used in the DNS. Tagawa et al. (Reference Tagawa, Roghair, Prakash, van Sint Annaland, Kuipers, Sun and Lohse2013) shows that the clustering indicator $\mathcal {C}=\varSigma /\varSigma _{{rnd}}$ quantitatively identifies different clustering morphologies. A value of $\mathcal {C}<1$ indicates a regular lattice arrangement, $\mathcal {C}=1$ a random bubble distribution and $\mathcal {C}>1$ irregular clustering. For our case with $\phi =1.7\,\%$, we obtain $\mathcal {C}=0.95$ and for $\phi =5.2\,\%$ we have $\mathcal {C}=1.2$. These values indicate random or weakly irregular clustering in our DNS (Pandey et al. Reference Pandey, Ramadugu and Perlekar2020). Consequently, we observe no significant effect of clustering on our results.

3.3. Velocity spectra

We define the velocity spectrum $E_u$ of the suspension at wavenumber $k=\lvert {\boldsymbol k}\rvert$ as

(3.1)\begin{equation} E_u(k) =\frac{1}{2} \sum_{k-{\rm \Delta} k/2 < \lvert{\boldsymbol k}\rvert < k+{\rm \Delta} k/2} \dfrac{\langle \hat{\boldsymbol{u}}'({\boldsymbol k}) \boldsymbol{\cdot} \hat{\boldsymbol{u}}'^*{(\boldsymbol k)} \rangle}{{\rm \Delta} k}, \end{equation}

where ‘$\hat { }$’ represents the Fourier mode, ‘$^*$’ is the complex conjugate and $\langle {\cdot } \rangle$ is the ensemble average. Similarly, we can define the velocity spectrum of the liquid phase considering the liquid velocity ${\boldsymbol {\tilde {u}_l}=\boldsymbol {u}\lvert {(\,f=1)}}$; however, we regularise the liquid velocity directly with the volume fraction field to avoid discontinuities that can lead to possible Gibbs phenomena in the Fourier transform defining the regularised liquid velocity as $\boldsymbol {u_l}=f\boldsymbol {u}$

(3.2)\begin{equation} {E}_{u,l}(k) =\frac{1}{2} \sum_{k-{\rm \Delta} k/2 < \lvert{\boldsymbol k}\rvert < k+{\rm \Delta} k/2}\dfrac{\langle \widehat{\boldsymbol{u_l}}'({\boldsymbol k}) \boldsymbol{\cdot} \widehat{\boldsymbol{u_l}}'^*{(\boldsymbol k)} \rangle}{{\rm \Delta} k}.\end{equation}

Figure 4(a) shows the suspension and liquid velocity spectra for our cases, and previous bubbly flow DNS by Pandey et al. (Reference Pandey, Ramadugu and Perlekar2020) at similar governing parameters (case R7 in that study). All cases show a peak at scales close to the bubble diameter ($k/k_{d_b}\approx 1$, with $k_{d_b}=2{\rm \pi} /d_b$) whereas at $k/k_{d_b} > 1$, the velocity spectra of the bubbly flows scale approximately as $k^{-3}$, in agreement with previous studies. The single-phase velocity spectra show an even steeper slope since, here, the velocity fluctuations are only produced at the forcing length scale. This is opposed to bubbly flows where velocity fluctuations are produced by both large scales, and at $k/k_{d_b} > 1$ where the fluctuations are continuously produced and directly dissipated in the bubble wakes (Lance & Bataille Reference Lance and Bataille1991). The single-phase velocity spectrum does not show the Kolmogorov $-$5/3 scaling in a large interval since the Taylor–Reynolds number is only 39. The $-3$ scaling in the bubbly flows is more pronounced in the liquid than in the suspension, so we focus on the scalar dynamics in the liquid phase.

Figure 4.

Normalised velocity spectra $E_u(k)$ and scalar spectra $E_{c,l}(k)$ against the normalised wavenumber $k/k_{d_b}$. (a) Velocity spectra of single-phase, bubble suspensions and liquid phases at two different volume fractions and comparisons with Pandey et al. (Reference Pandey, Ramadugu and Perlekar2020). (b) Value of $E_c(k)$ for the single phase isotropic turbulence case. (c) Bubbly flow cases with $\nabla ^v \langle c \rangle$. (d) Bubbly flow cases with $\nabla ^h \langle c \rangle$. The gas volume fraction is $\phi =1.7\,\%$ in all cases A1–A5. Case B4 with $\phi =5.2\,\%$ is included to illustrate the effects of $\phi$ on $E_{c,l}(k)$.

3.4. Scalar spectra

We define the suspension and liquid scalar spectra similar to (3.1) and (3.2)

(3.3a,b)\begin{align} E_c(k) = \sum_{k-{\rm \Delta} k/2 < \lvert{\boldsymbol k}\rvert < k+{\rm \Delta} k/2} \dfrac{\langle \hat{{c}}'({\boldsymbol k}) \hat{{c}}'^*{(\boldsymbol k)} \rangle}{{\rm \Delta} k}, \quad E_{c,l}(k) = \sum_{k-{\rm \Delta} k/2 < \lvert{\boldsymbol k}\rvert < k+{\rm \Delta} k/2} \dfrac{\langle \widehat{{c_l}}'({\boldsymbol k}) \widehat{{c_l}}'^*{(\boldsymbol k)} \rangle}{{\rm \Delta} k}. \end{align}

The single-phase scalar spectra are shown in figure 4(b). Below the forcing scale, the classical $-5/3$ scaling emerges, especially for the highest $Sc$. For lower $Sc$, the diffusion significantly influences the dynamics below the forcing scale, like in the velocity spectra, decreasing exponentially due to viscous dissipation (Pope Reference Pope2001).

The normalised liquid scalar spectra for the bubbly flow simulations are shown in figure 4(c) for the cases with $\nabla ^v \langle c \rangle$ and in figure 4(d) for the cases with $\nabla ^h \langle c \rangle$. For brevity, we only show the spectra for all cases at $\phi =1.7\,\%$ since the general trends are the same at $\phi =5.2\,\%$. Case B4 is, however, included to illustrate the minor effects of a higher $\phi =5.2\,\%$ on $E_c/(c^2_{{std}}d_b)$. We do not see the different scaling laws at different volume fractions as in the experiments by Dung et al. (Reference Dung, Waasdorp, Sun, Lohse and Huisman2023). However, in that experiments, the liquid turbulence results from the interaction between the bubble-induced turbulence and the incident turbulence generated by an active grid. This experimental system corresponds to a finite value of the bubblance parameter (ratio of the kinetic energy produced by the bubbles to the turbulent kinetic energy in the absence of bubbles, Rensen, Luther & Lohse Reference Rensen, Luther and Lohse2005). This is in contrast with the present DNS results with only bubble-induced turbulence and an infinite bubblance parameter. Indeed, Prakash et al. (Reference Prakash, Mercado, van Wijngaarden, Mancilla, Tagawa, Lohse and Sun2016) showed experimentally how the $-3$ scaling of the velocity spectra emerged as the bubblance parameter was modified. This difference between the present DNS and the experiment by Dung et al. (Reference Dung, Waasdorp, Sun, Lohse and Huisman2023) is a possible reason why we do not observe different scaling laws at different volume fractions. In addition, the latter experimental study considers a thermal mixing layer with a nonlinear mean gradient, whereas the present study imposes a constant mean scalar gradient. The influence of a nonlinear gradient on the scalar spectrum is not clear but may also explain the different observations.

In figures 4(c) and 4(d), we observe a maximum of the scalar spectra at the lowest wavenumber, a signature that the scalar is directly forced by the mean gradients. Consistent with the results of Dung et al. (Reference Dung, Waasdorp, Sun, Lohse and Huisman2023), the liquid scalar spectra show two different scalings. The larger scales show a standard $-5/3$ scaling only in a narrow range of wavenumbers. However, after a transition scale, the spectrum scaling changes towards the $-$3 scaling or lower for all the Schmidt numbers. When the gradient is in the horizontal direction (figure 4d), we observe a faster decay of the scalar spectrum below the bubble diameter wavenumber. The transition wavenumber increases with the Schmidt number, depending on the molecular diffusivity. It coincides approximately with the bubble wavenumber for the smaller $Sc$ since, at $Sc$ of order one, we expect that the dynamics of the passive scalar transport is similar to the momentum transport. Case A5 is included to assess the effects of increasing $Sc_g$ on the scalar dynamics while keeping $Sc_l=0.7$ as in case A1; this case corresponds to the same scalar diffusivity in the gas and liquid phases. We observe that also, in this case, the normalised spectrum shows a $-$3 scaling after the bubble diameter wavenumber. Figures 4(c) and 4(d) show only minor differences between cases A1 and A5, indicating that $Sc_g$ does not significantly influence the liquid scalar fluctuations. This is, however, not true for the statistics of the suspension, as discussed later.

Figure 5 shows the compensated scalar spectra to more clearly illustrate the scalar spectrum scalings and transitions. Figure 5(a) shows that the classical $-5/3$ scaling in the single-phase cases emerges in a narrow range of wavenumbers (where the range increase with $Sc_l$) after the forcing length scale. This limited range is due to an insufficient Taylor–Péclet number. Defining the Taylor length scale for the scalar field as $\lambda _c = \sqrt {6D_l/\varepsilon _c}c_{l,{std}}$ we obtain Taylor–Péclet numbers ${Pe}_\lambda = u_{l,{std}} \lambda _c / D_l$ in the range of 33–143 and Taylor length scales $\lambda _c$ in the range of 0.15–0.6$d_b$. These values indicate that molecular diffusion significantly affects the scalar dynamics at scales comparable to or slightly smaller than the bubble size $d_b$. It is thus reasonable that the inertial range with the $-5/3$ scaling is not very pronounced. The $-5/3$ scaling is, however, well established for single-phase DNS with wider inertial ranges (Corrsin Reference Corrsin1951; Gotoh & Watanabe Reference Gotoh and Watanabe2012). Since the forcing length scale in the single phase cases is close to the bubble size, it is reasonable that the $-5/3$-scaling is not pronounced in the bubbly flow cases either, as shown in the compensated spectra of figure 5(b).

Figure 5.

Compensated and normalised scalar spectra $E_{c,l}(k)$ against the normalised wavenumber $k/k_{d_b}$. (a) Single-phase isotropic turbulence cases compensated by $k^{5/3}$. (b) Bubbly flow cases with $\nabla ^h \langle c \rangle$ compensated by $k^{5/3}$. (c) Bubbly flow cases with $\nabla ^v \langle c \rangle$ compensated by $k^{3}$.

Figure 5(c) shows the compensated scalar spectra for the bubble flow cases with $\nabla ^v \langle c \rangle$. Here, it is clear that the spectra change scaling towards approximately $-3$ at scales $k/k_{d_b} \geq 1$ and that the transition occurs at higher $k$ for increasing $Sc_l$.

3.5. Scalar spectral budget

The transition to the $-$3 scaling is clearly a footprint of the bubble-induced agitation mechanisms that influence the scalar dynamics also at scales below the bubble diameter $k/k_{d_b}\gg 1$. To investigate the origin of this scaling, we compute the spectral budget of the passive scalar by manipulating (2.4) in Fourier space (Monin & Yaglom Reference Monin and Yaglom1975)

(3.4)\begin{equation} \dfrac{\partial E_c(k,t)}{\partial t} = T_c(k,t) - \mathcal{E}_c(k,t) + P_c(k,t),\end{equation}

where $T_c(k,t) = \sum _{(k-{\rm \Delta} k/2 < \lvert {\boldsymbol k}\rvert < k+{\rm \Delta} k/2)} -\langle \widehat {h_c}\hat {c}'^* + \widehat {h_c}^*\hat {c}' \rangle / {\rm \Delta} k$ is the local transfer term, $\mathcal {E}_c(k,t) = \sum _{(k-{\rm \Delta} k/2 < \lvert {\boldsymbol k}\rvert < k+{\rm \Delta} k/2)} \langle \widehat {d_c}\hat {c}'^* + \widehat {d_c}^*\hat {c}' \rangle / {\rm \Delta} k$ is the dissipation term and $P_c(k,t) = \sum _{(k-{\rm \Delta} k/2 < \lvert {\boldsymbol k}\rvert < k+{\rm \Delta} k/2)} \boldsymbol {\nabla } \langle c \rangle \boldsymbol {\cdot } \langle (\hat {\boldsymbol {u}}'\hat {c}'^* + \hat {\boldsymbol {u}}'^*\hat {c}')\rangle / {\rm \Delta} k$ is the production at wavenumber $k$.

Dung et al. (Reference Dung, Waasdorp, Sun, Lohse and Huisman2023) speculated that the passive scalar transport mechanisms are not due to direct scalar production and simultaneous dissipation at scales smaller than the bubble diameter, as for the momentum transport. Instead, scalar production, generated at larger scales by the mean gradient, decays faster so that at small scales just a balance between liquid scalar transfer $T_{c,l}$ and spectral dissipation $\mathcal {E}_{c,l}(k,t)$ occurs. By dimensional arguments, they show that, if the liquid scalar transfer is only a function of scalar dissipation and wavenumber, then the scaling is a $-1$ power law for the wavenumber $T_{c,l}=T_{c,l}(\varepsilon _{c,l},k)\propto \varepsilon _{c,l}k^{-1}$. Since the liquid spectral dissipation is defined as $\mathcal {E}_{c,l}(k,t)=2D_lk^2E_{c,l}$, it is trivial to show that if a balance exists with the scalar transfer term, then the scalar spectra should scale as a $-$3 power law in $k, E_{c,l} \propto k^{-3}$. The different terms of (3.4) are challenging to measure in experiments, so, as suggested in Dung et al. (Reference Dung, Waasdorp, Sun, Lohse and Huisman2023), DNS are needed to estimate them. Here, we assess their hypotheses, extracting the relevant statistics from our numerical simulations.

In figure 6(a,b,c), we show the normalised production $P_{c,l}(k)/(\epsilon _{c,l} d_b)$ of (3.4) for the isotropic single-phase (a) and bubbly flow cases with $\nabla ^v \langle c \rangle$ (b) and $\nabla ^h \langle c \rangle$ (c). The panels show that a significant portion of the scalar fluctuations is produced at approximately $k/k_{d_b} \leq 1$ for both the single-phase and bubbly flow cases. At $k/k_{d_b} > 1, P_{c,l}(k)$ scales as approximately $-7/3$ in a narrow range of $k$ for the single-phase cases (as predicted by Lumley (Reference Lumley1964) with dimensional analysis and found in the DNS of Gotoh & Watanabe Reference Gotoh and Watanabe2012). Bubble suspensions, instead, show that $P_{c,l}(k)$ decays with a faster power law (approximately $k^{-3}$). Therefore, at a statistically steady state, (3.4) implies that the scalar spectra $E_{c,l}(k)$ at large wavenumbers are governed by the balance of the net transfer and the diffusive dissipation. Interestingly, we observe a self-similarity in the normalised production spectra for all the cases at different liquid diffusivity. This self-similarity can probably be explained by the fact that the production at large scales is due to the mean scalar gradient, while the liquid diffusivity only influences the smaller scales where the production is small.

Figure 6.

Normalised scalar production term $P_{c,l}(k)/(\epsilon _{c,l} d_b)$ for the single-phase cases (a) and bubbly flow cases A1–A5 with $\phi =1.7\,\%$ for $\nabla ^v \langle c \rangle$ (b) and $\nabla ^h \langle c \rangle$. (c) Normalised scalar transfer term $|T_{c,l}(k)|/(\epsilon _{c,l} d_b)$ for the corresponding single-phase cases (d) and bubbly flow cases with $\nabla ^v \langle c \rangle$ (e) and $\nabla ^h \langle c \rangle$ (f).

The bottom panels of figure 6 show the absolute value of the scalar transfer spectra. When the production becomes significantly smaller than the transfer, the latter scales close to $-1$ in a limited range of wavenumbers. In the same range, the scalar spectra show approximately a $-3$ scaling. These trends are highlighted in the zoomed-in plots of figure 7 that show the scalar spectrum, absolute transfer and the production spectrum for each case A1–A4 with $\nabla ^v \langle c \rangle$. These plots are consistent with the proposed hypothesis in Dung et al. (Reference Dung, Waasdorp, Sun, Lohse and Huisman2023) that, at large wavenumbers ($k/k_{d_b} \gtrapprox 3$ in figure 7), the scalar spectrum is governed by the balance of the net transfer and the diffusive dissipation (as the production is negligible) and that the $-1$ scaling of the transfer term induces the $-3$ scaling of the scalar spectra. The trends shown in figure 7 are observed for all our bubble suspension cases.

Figure 7.

Zoomed-in plots of the scalar, transfer and production spectra for our cases A1–A4 with $\nabla ^v \langle c \rangle$. In the same range where $|T_{c,l}(k)| \propto k^{-1}$, the scalar spectra show a transition to approximately $E_{c,l} \propto k^{-3}$.

3.6. Effective scalar diffusivity

The effective scalar diffusivity describes the macroscale scalar flux in the suspension seen as a continuum. Following the detailed derivation of Loisy (Reference Loisy2016), the contributions to the effective diffusivity are obtained by ensemble averaging the scalar transport equation

(3.5)\begin{equation} \dfrac{\partial \langle c \rangle }{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \langle \boldsymbol{q} \rangle = 0 , \end{equation}

where the average flux can be formulated as

(3.6)\begin{align} \langle \boldsymbol{q} \rangle = \langle\boldsymbol{u}\rangle \langle c\rangle - D_{{mol,s}} \boldsymbol{\nabla}\langle c\rangle + (1-\phi)\langle \boldsymbol{u}'c' \rangle_l + \phi\langle \boldsymbol{u}'c' \rangle_g - (1-\phi)D_l \langle \boldsymbol{\nabla} c' \rangle_l -\phi D_g \langle \boldsymbol{\nabla} c' \rangle_g . \end{align}

Here, we have used that the phase-ensemble average of an arbitrary variable $A_s$ of the suspension is $\langle A_s \rangle = \langle f A_l + (1-f)A_g \rangle = (1-\phi )\langle A_l \rangle + \phi \langle A_g \rangle$. The effective diffusivity tensor is defined as

(3.7)$$\begin{gather} {\boldsymbol{\mathsf{D}}}_{{eff}} = D_{{mol,s}} {\boldsymbol{\mathsf{I}}} + {\boldsymbol{\mathsf{D}}}_{{conv}} + {\boldsymbol{\mathsf{D}}}_{{diff}}, \end{gather}$$

(3.8)$$\begin{gather}{\boldsymbol{\mathsf{D}}}_{{conv}} \boldsymbol{\cdot} \boldsymbol{\nabla} \langle c \rangle =-(1-\phi)\langle \boldsymbol{u}'c' \rangle_l - \phi\langle \boldsymbol{u}'c' \rangle_g, \end{gather}$$

(3.9)$$\begin{gather}{\boldsymbol{\mathsf{D}}}_{{diff}}\boldsymbol{\cdot} \boldsymbol{\nabla} \langle c \rangle = (1-\phi)D_l \langle \boldsymbol{\nabla} c' \rangle_l +\phi D_g \langle \boldsymbol{\nabla} c' \rangle_g, \end{gather}$$

so that (3.5) becomes

(3.10)\begin{equation} \dfrac{\partial \langle c \rangle }{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\langle\boldsymbol{u}\rangle \langle c\rangle) - \boldsymbol{\nabla} \boldsymbol{\cdot} ({\boldsymbol{\mathsf{D}}}_{{eff}} \boldsymbol{\cdot} \boldsymbol{\nabla} \langle c\rangle) =0 . \end{equation}

For all cases in this study ${\boldsymbol{\mathsf{D}}}_{{diff}} < 0.02{\boldsymbol{\mathsf{D}}}_{{conv}}$. The value of ${\boldsymbol{\mathsf{D}}}_{{conv}}$ is larger also in Loisy et al. (Reference Loisy, Naso and Spelt2018) for lower bubble Reynolds numbers ($Re\approx 30$) and with the same diffusivity in both phases. We define ${\boldsymbol {D}}^v_{{conv}}$ as the non-zero diagonal component of ${\boldsymbol{\mathsf{D}}}_{{conv}}$ for an imposed scalar gradient $\nabla ^v \langle c \rangle$ and ${\boldsymbol {D}}^h_{{conv}}$ for $\nabla ^h \langle c \rangle$. All off-diagonal components of ${\boldsymbol{\mathsf{D}}}_{{conv}}$ are zero (Loisy Reference Loisy2016).

We define the Sherwood number as ${Sh} = {\boldsymbol {D}}_{{conv}} / D_{{mol,s}}$ that represents the relative importance of the convective contribution to the molecular contribution in the effective diffusivity of the scalar (in heat transfer, the same definition applies to the Nusselt number). In single-phase isotropic turbulence ${Sh} \propto {Pe}_{{std}} = u_{{std}}l_u/D_{{mol,s}}$ according to theory and DNS (Gotoh & Watanabe Reference Gotoh and Watanabe2012) at high Péclet numbers, where diffusion is negligible ($l_u$ is the velocity integral length scale). The same scaling is found at high ${Pe_d=V_d d_b / D_{{mol,s}}}$ (for cases with both $\nabla ^v \langle c \rangle$ and $\nabla ^h \langle c \rangle$) for bubbly flows without fully developed bubble-induced turbulence and with a constant scalar diffusivity (Loisy et al. Reference Loisy, Naso and Spelt2018). Figure 8(a) shows that, indeed, ${Sh} \propto {Pe}_{{std}}$ for our single-phase and bubbly flow simulations with $\nabla ^h \langle c \rangle$. Here, we have used $l_u = d_b$ as the length scale for ${Pe}_{{std}}$ in the bubbly flow cases and the forcing length scale in the single-phase simulations. The higher ${Sh}$ for cases with $\nabla ^v \langle c \rangle$ (compared with cases with $\nabla ^h \langle c \rangle$) can be explained by the preferential vertical motion of the bubbles due to buoyancy. The ratio of the liquid velocity components in the vertical and horizontal directions is $u^v_{l,{std}}/u^h_{{l,std}}\approx 1.4$ (where we obtain for $\phi =1.7\,\%: u^v_{l,{std}}=0.29, u^h_{l,{std}}=0.21$ and for $\phi =5.2\,\%: u^v_{l,{std}}=0.42, u^h_{l,{std}}=0.31$) for all our bubbly flow cases and where the ratios of the vertical to horizontal components are in excellent agreement with experiments (Riboux et al. Reference Riboux, Risso and Legendre2010). Contrarily, in single-phase isotropic turbulence we have $u_{l,{std}}=u^v_{l,{std}}=u^h_{l,{std}}=0.3$. These values suggest that the scalar flux in the single-phase turbulence should be similar to or higher than that in the bubbly flow with $\nabla ^h \langle c \rangle$ (where the scalar flux is proportional to $u^h_{l,{std}}$). The scalar flux in the single-phase turbulence should, however, be less than that in the bubbly flow with $\nabla ^v \langle c \rangle$ (where the total scalar flux is governed by $u^v_{l,{std}}$ and the additional mechanisms acting in the vertical direction discussed in § 3.1).

Figure 8.

Convective contribution to the effective diffusivity for all our simulation cases. The darker blue squares in (a) and (b) indicate the cases A5 and B5 with $\nabla ^v \langle c \rangle$ and $Sc_g=7$. (a) Value of ${Sh}$ against ${Pe}_{{std}}$ based on the total standard deviation of velocity fluctuations in the liquid; (b) ${Sh}$ against ${Pe}_{\phi }$ based on a priori known parameters of the suspension. (c) Non-dimensional ${\boldsymbol {D}}^*_{{conv}} = {\boldsymbol {D}}_{{conv}} / (\sqrt {g d_b}d_b)$ against $\phi$ for $\nabla ^v \langle c \rangle$ (top panel) and $\nabla ^h \langle c \rangle$ (bottom panel). In both panels, the colour of the points represents the same ${Sc}$ as in figure 6. In the top panel, the proposed scalings for ${\boldsymbol {D}}^v_{{conv}}$ are the dash-dotted line (lowest $D_{{mol,s}}$, case A4) and the dashed line (highest $D_{{mol,s}}$, case B1). The dashed line in the bottom panel is the proposed scaling for ${\boldsymbol {D}}^h_{{conv}}$.

In the bubbly flow cases with $\nabla ^v \langle c \rangle$ we observe approximately ${Sh} \propto {Pe}^{1.15}_{{std}}$, indicating that molecular diffusion is not negligible. Using the scaling $u_{{std}} \propto V_0\phi ^{0.4}$ found in Risso & Ellingsen (Reference Risso and Ellingsen2002) and Riboux et al. (Reference Riboux, Risso and Legendre2010) and noting that the rise velocity of a single bubble $V_0 / (\sqrt {g d_b})=O(1)$ we can define $Pe_{\phi } = \sqrt {g d_b}\phi ^{0.4}d_b/D_{{mol,s}}$ based on a priori known parameters. Figure 8(b) shows the ${Sh}$ against ${Pe}_\phi$ for our bubbly flow cases where we again observe ${Sh}^v \propto {Pe}^{1.15}_{\phi }$ and ${Sh}^h \propto {Pe}^{1}_{\phi }$. Solving for ${\boldsymbol {D}}_{{conv}}$ in ${Sh}$ of the latter scalings gives ${\boldsymbol {D}}^v_{{conv}} \propto (\sqrt {g d_b}d_b)^{1.15}\phi ^{0.46}/D_{{mol,s}}^{0.15}$ and ${\boldsymbol {D}}^h_{{conv}} \propto \sqrt {g d_b}\phi ^{0.4}d_b$. These are similar dependencies on $\phi$ as found experimentally in Alméras et al. (Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015) for a low-diffusive dye (${\boldsymbol {D}}_{{conv}}\propto \phi ^{0.4}$ at lower $\phi$) and by Gvozdić et al. (Reference Gvozdić, Alméras, Mathai, Zhu, van Gils, Verzicco, Huisman, Sun and Lohse2018) for heat transport (${\boldsymbol {D}}_{{conv}} \propto \phi ^{0.45}$ up to $\phi =5\,\%$).

The value of ${\boldsymbol {D}}_{{conv}}$ for each simulation case are shown in figure 8(c). In cases A1–A4 and B1–B4 (coloured dots), $Sc_l$ is increased from $0.7$ to $7$ while $Sc_g=0.7$ is maintained. For these cases there is a monotonic increase of both ${\boldsymbol {D}}^v_{{conv}}$ and ${\boldsymbol {D}}^h_{{conv}}$ with $Sc_l$. However, in cases A5 and B5 (grey crosses) we increase $Sc_g=7$ while specifying $Sc_l=0.7$. These parameters result in the highest ${\boldsymbol {D}}^v_{{conv}}$ for the respective cases A1–5 and B1–5 but a ${\boldsymbol {D}}^h_{{conv}}$ closer to lower $Sc_l$.

The influence of $Sc_g$ on ${\boldsymbol{\mathsf{D}}}_{{conv}}$ can be explained by considering the bubbles as a source of scalar disturbance $c'$. When a bubble moves in the direction of the imposed scalar gradient, the scalar disturbance in the bubble increases proportional to the last term on the right-hand side of (2.4). The time it takes for the scalar in the bubble to reach an equilibrium with the liquid is proportional to the characteristic equalisation time $t_e = R^2/D_{g}$. A higher $Sc_g$ thus implies a longer $t_e$ and therefore a higher average value of $c'$ in the gas phase contribution $\phi \langle \boldsymbol {u}'c' \rangle _g$ to ${\boldsymbol{\mathsf{D}}}_{{conv}}$ in (3.8). In our cases, the liquid phase contribution to (3.8) is almost independent of $Sc_g$. The influence of $Sc_g$ on ${\boldsymbol{\mathsf{D}}}_{{conv}}$ is evident in the cases with $\nabla ^v \langle c \rangle$ since the bubbles have an average velocity in the vertical direction due to buoyancy. The high average gas phase velocity in the gradient direction induces large average gas phase scalar fluctuations. The high gas phase velocity and scalar fluctuations cause a significant gas phase contribution $\phi \langle \boldsymbol {u}'c' \rangle _g$ to ${\boldsymbol {D}}^{v}_{{conv}}$ although the parameter $\phi$ is small.

The proposed scalings for ${\boldsymbol {D}}_{{conv}}$ are shown in figure 8(c) with a proportionality constant of $0.94$ for ${\boldsymbol {D}}^v_{{conv}}$ and $0.25$ for ${\boldsymbol {D}}^h_{{conv}}$. In the top panel, the dash-dotted line represents the scaling for ${\boldsymbol {D}}^v_{{conv}}$ using the lowest $D_{{mol,s}}$ (case A4) and the dashed line is the same scaling using the highest $D_{{mol,s}}$ (case B1). All other cases with $Sc_g=0.7$ fall between these limits. Because of the previously discussed effects of $Sc_g$ on ${\boldsymbol {D}}^v_{{conv}}$, the results for $Sc_g=7$ (grey crosses) are not in such good agreement with the proposed scalings. More data are, however, needed to study the effects of $Sc_g$ on ${\boldsymbol {D}}^v_{{conv}}$. The bottom panel of figure 8(c) shows the scaling for ${\boldsymbol {D}}^h_{{conv}}$ and all cases show that the effects of molecular diffusion (${Sc}$) are not as significant.