2. Simulation of the bubble swarm

Although the equations describing precisely this type of flow are relatively well known, their numerical simulation remains out of reach, due to the large spectrum of temporal and spatial scales involved. The smallest scales are a priori associated with the interfacial dynamics and the development of a very thin boundary layer around the bubbles, while the largest scales are related to the length of the wakes and the evolution of the collective dynamics of the flow that takes place. In order to simulate these flows, we use the approach proposed by Le Roy De Bonneville et al. (Reference Le Roy De Bonneville, Zamansky, Risso, Boulin and Haquet2021). This approach abandons the precise description of the flow around the bubbles as well as the capillary effects while keeping a realistic dynamic of the downstream part of each wake and enables a straightforward analysis of the structure of the liquid velocity field and the dynamics of the bubble swarm. This modelling, based on the Euler–Lagrange approach, allows us to simulate flows with a large number of bubbles and to focus on the interactions between wakes. As we briefly recall below, the main difficulty of this type of calculation comes from the self-interaction of a bubble with its own wake. Le Roy De Bonneville et al. (Reference Le Roy De Bonneville, Zamansky, Risso, Boulin and Haquet2021) proposed a method enabling taking into account this effect and to accurately calculate the trajectory of each bubble. This method allowed us to obtain numerical simulations of the turbulence induced by a swarm of bubbles as illustrated in figure 1. We show later, in § 3, that the flow structure predicted by this approach is in good agreement with experiments.

Figure 1. Snapshots of the vertical component of the liquid velocity field in a vertical plane for (a) $\alpha =2\,\%$ and (b) $\alpha =10\,\%$. The blue points represent the positions of the bubbles.

2.1. Modelling

The use of the Euler–Lagrange approach amounts to considering a filtering of the flow field near the bubbles. In this approach, the action of the dispersed phase on the flow is introduced as a volume source of momentum localized around the bubbles. The liquid velocity is given by the Navier–Stokes equations:

(2.1)\begin{equation} D_t \boldsymbol{u}_f ={-}\dfrac{1}{\rho_f}\boldsymbol{\nabla} p_f+\nu {\rm \Delta} \boldsymbol{u}_f + \dfrac{\boldsymbol{f}}{\rho_f} ; \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}_f = 0,\end{equation}

where $\boldsymbol {u}_f$ represents a filtered velocity field around the bubbles, $\nu$ the kinematic viscosity of the liquid and $\rho _f$ its density. Note that $p_f$ is the pressure variation relative to the hydrostatic pressure $P_h$ with $-1/\rho _f \boldsymbol {\nabla } P_h = \langle \boldsymbol {f} \rangle$, which ensures that the computation is carried out in the frame where the average fluid velocity is zero. These equations are supplemented by tri-periodic boundary conditions and their numerical solutions are obtained by a spectral method as detailed in Le Roy De Bonneville et al. (Reference Le Roy De Bonneville, Zamansky, Risso, Boulin and Haquet2021) and Zamansky (Reference Zamansky2022).

The volume forcing of the liquid phase in (2.1) is given by

(2.2)\begin{equation} \boldsymbol{f}(\boldsymbol{x},t) ={-} \sum_{b=1}^{N_b} \boldsymbol{F}_{f\rightarrow b}(t) G_\sigma(\boldsymbol{x}-\boldsymbol{x}_b(t)),\end{equation}

where $\boldsymbol {F}_{f\rightarrow b}$ is the momentum exchange rate between the fluid and the bubble $b$, and $N_b$ is the number of bubbles. Term $G_\sigma$ is the Gaussian kernel of the projection and $\sigma$ is its characteristic size which is of the order of the diameter of the bubble. The latter is thus much larger than the mesh size: $\sigma \approx d > {\rm \Delta} x$. Indeed, although the details close to the bubbles are filtered, the flow presents a priori scales much smaller than the bubble size. These small-scale fluctuations result from the evolution of the turbulent wakes and their interactions.

The trajectory of the bubbles is obtained by solving Newton's equation for each bubble. It involves the hydrodynamic force which depends on the velocity of the liquid (and its derivatives). Le Roy De Bonneville et al. (Reference Le Roy De Bonneville, Zamansky, Risso, Boulin and Haquet2021) considered that the bubble is subject to the drag force, the added-mass force, the inertia force of the fluid as well as the buoyancy. We have not retained the history force because it has a priori a negligible effect for large Reynolds numbers. On the other hand, when the velocity gradient is large at the scale of the bubble, the lift forces can certainly play a role. Similarly, the anisotropic effects of drag and added mass related to a non-spherical bubble are also important. We aim to reproduce the experiments of Riboux et al. (Reference Riboux, Risso and Legendre2010) for millimetric air bubbles in water. Given the Reynolds number of the bubbles and the Morton number, the bubbles clearly adopt a non-spherical shape (Maxworthy et al. Reference Maxworthy, Gnann, Kürten and Durst1996). However, to simplify the modelling of the problem, we consider that the bubbles are spherical, assuming that in the case of the homogeneous swarm, the anisotropic aspects are not essential. Consistently, we have as well not retained the lift force. Note that the value of the lift coefficient, and even its sign, being very dependent on the shape of the object, it would be delicate to choose its value anyway. Finally, considering that the density of the gas is very low, we obtain for the dynamics of the bubble:

(2.3)\begin{equation} C_M \dfrac{{\rm d} \boldsymbol{v}_b}{{\rm d} t} = \dfrac{3C_D}{4 d} (\boldsymbol{v}_b - \tilde{\boldsymbol{u}}_{f,b} )|\boldsymbol{v}_b - \tilde{\boldsymbol{u}}_{f,b}| +(1+C_M) \dfrac{{\rm D} \tilde{\boldsymbol{u}}_{f,b}}{{\rm D} t} - \boldsymbol{g} + \boldsymbol{F}_{{I},b}.\end{equation}

The drag coefficient is chosen, in agreement with the experiment of an isolated bubble, at $C_D = 0.35$ and the added-mass coefficient at $C_M=0.5$ in coherence with the spherical bubble hypothesis. In this equation, the bubble force is calculated from the corrected liquid velocity $\tilde {\boldsymbol {u}}_{f,b}$, as proposed in Le Roy De Bonneville et al. (Reference Le Roy De Bonneville, Zamansky, Risso, Boulin and Haquet2021), and briefly explained below. Note that interactions between bubbles are accounted for through the liquid disturbances generated by the bubbles, which is well suited for the study of homogeneous bubbly flows, but should probably be improved by adjusting of the drag and added-mass coefficients according to the local bubble distribution, in the spirit of the method proposed by Akiki, Jackson & Balachandar (Reference Akiki, Jackson and Balachandar2017).

Finally, the term $\boldsymbol {F}_{{I},b}$ is a repulsive force between bubbles. It is introduced to prevent the bubbles from overlapping and to ensure that the distance between bubbles remains greater than the characteristic size $\sigma$ of the momentum source. The force depends on the distance between each pair of bubbles $r_{b,b'}= | \boldsymbol {x}_b-\boldsymbol {x}_{b'}|$ and is given by $\boldsymbol {F}_{{I},b} = \sum _{b'{\neq} b } -C_I (({\boldsymbol {x}_b-\boldsymbol {x}_{b'}})/{r_{b,b'}})\exp (- r_{b,b'}^2/2r_I^2)$. In the numerical simulation presented in the paper, we use the same value as prescribed in Le Roy De Bonneville et al. (Reference Le Roy De Bonneville, Zamansky, Risso, Boulin and Haquet2021). It has been verified that the modification of the values of these parameters does not modify significantly the simulations.

The momentum exchanged between the bubble and the fluid via the volume force $\boldsymbol {f}$ in (2.2) is given by the sum of the drag force and the added-mass force, the contributions of the Tchen and Archimedes forces being already taken into account in the pressure term in a way consistent with the zero divergence of the flow (Climent & Magnaudet Reference Climent and Magnaudet1999; Le Roy De Bonneville et al. Reference Le Roy De Bonneville, Zamansky, Risso, Boulin and Haquet2021).

The fluid velocity, corrected for the influence of bubble $b$, is defined by introducing the perturbation due to the bubble:

(2.4)\begin{equation} \tilde{\boldsymbol{u}}_{f,b}(x,t)= \boldsymbol{u}_{f}-\boldsymbol{u}^*_{f,b} .\end{equation}

Because of the nonlinearity of the system, this immediately raises the question of the definition of the perturbation $\boldsymbol {u}^*_{f,b}$. It is indeed not trivial to isolate the influence of a bubble among the fluctuations of the flow which include the effect of all the other bubbles. We propose here to define the perturbed field $\boldsymbol {u}^*_{f,b}$ as the flow generated by an isolated fictitious bubble, in a liquid at rest, which would have followed the same trajectory and exchanged as much momentum with the liquid phase as the actual bubble $b$. Le Roy De Bonneville et al. (Reference Le Roy De Bonneville, Zamansky, Risso, Boulin and Haquet2021) proposed an integral model to calculate $\boldsymbol {u}^*_{f,b}$, and its derivatives, valid for the case of bubbles at large Reynolds number.

In this model, the main assumptions to obtain $\boldsymbol {u}^*_{f,b}$ are that in the vicinity of the bubble (i) given the importance of the Reynolds number, the viscous term is neglected and (ii) the flow is considered quasi-parallel. The details of the derivation can be found in Le Roy De Bonneville et al. (Reference Le Roy De Bonneville, Zamansky, Risso, Boulin and Haquet2021), but after a few steps we obtain the following integral expression for $\boldsymbol {u}^*_{f,b}$:

(2.5)\begin{equation} \boldsymbol{u}^*_{f,b}{\scriptstyle(\boldsymbol{x},t)} = \dfrac{1}{\rho_f} \int_{0}^t \boldsymbol{F}_{f\rightarrow b}{\scriptstyle(s)} G_{\sigma}{\scriptstyle(\boldsymbol{x}-\boldsymbol{x}_b(s) + \boldsymbol{\ell}_{adv}(t,s) )}\, {\rm d} s ,\end{equation}

where $\boldsymbol {\ell }_{adv}(t,s) = \int _s^t \tilde {\boldsymbol {u}}_{f,b} {(\boldsymbol {x}=\boldsymbol {x}_b(s^{\prime }),s^{\prime })}\,{\rm d} s^{\prime }$. The velocity perturbation at a given position and time is obtained by integrating, over all previous instants, the momentum supplied by the bubble at the material point of the liquid at that specific position. The material point corresponding to the injection of momentum at an instant $s$ can be advected by the undisturbed flow, and will be found at a distance $\boldsymbol {\ell }_{adv}(t,s)$ at the instant $t>s$. This $\boldsymbol {\ell }_{adv}$ term is essential to guarantee the Galilean invariance of the model.

We do not detail here the discretization of (2.5). The details of its numerical implementation can be found in Le Roy De Bonneville et al. (Reference Le Roy De Bonneville, Zamansky, Risso, Boulin and Haquet2021). We just mention that we have developed an efficient algorithm to compute the history integral (see also Zamansky Reference Zamansky2022), so the extra cost of computing the correction term is negligible.

It is interesting to note that, although the coupling between the two phases conserves momentum, it does not conserve energy (Xu & Subramaniam Reference Xu and Subramaniam2007; Subramaniam et al. Reference Subramaniam, Mehrabadi, Horwitz and Mani2014; Le Roy De Bonneville et al. Reference Le Roy De Bonneville, Zamansky, Risso, Boulin and Haquet2021). Indeed, the power $P_b$ of the hydrodynamic force working at the bubble velocity is greater than the power $P_f$ of the diffuse force working at the fluid velocity:

(2.6)\begin{equation} P_b=\sum_b \boldsymbol{F}_{f\rightarrow b }\boldsymbol{\cdot} \boldsymbol{v}_b > P_f=\int {{\rm d}\kern 0.06em x}^3 \boldsymbol{f}\boldsymbol{\cdot}\boldsymbol{u}_f . \end{equation}

From a physical perspective, it is acceptable that energy is dissipated during the coupling. We consider a coarse description of the continuous phase, in which the strong velocity gradients in the close vicinity of the bubbles are not described. The dissipation of kinetic energy into heat that occurs in the region surrounding the bubble at scales smaller than $\sigma$ cannot be calculated from the resolved velocity $\boldsymbol {u}_f$. In the case of an isolated bubble, Le Roy De Bonneville et al. (Reference Le Roy De Bonneville, Zamansky, Risso, Boulin and Haquet2021) estimated analytically that

(2.7)\begin{equation} \dfrac{P_f}{P_b} = \dfrac{C_D }{64} \left( \dfrac{d_b}{\sigma}\right)^2. \end{equation}

Consequently, an important part of the mechanical energy is dissipated around the bubbles, in the boundary layer and the near wake.

3. Comparison with experiments

Figure 2 shows the PDFs of the horizontal and vertical components of the liquid velocity obtained by simulations for different $\alpha$ and by the experiments of Riboux et al. (Reference Riboux, Risso and Legendre2010). It is found that for both components the PDFs present an exponential decay and that the PDFs of the vertical component are clearly asymmetric. We further observe that the normalized PDFs are nearly invariant with $\alpha$ as the central part of the experimental PDFs. This behaviour is the signature of the turbulence induced by the interactions between wakes. For large fluctuations, the experimental PDFs show a second region characterized by a less steep exponential decay. This behaviour has been attributed to the large localized fluctuations very close to the bubbles (Risso Reference Risso2016). As this region is not described in our modelling, we indeed find that the second exponential part of the PDFs is not reproduced by the simulations.

Figure 2. The PDF liquid velocity for the horizontal velocity component (a) and the vertical velocity component (b) from simulations at various $\alpha$ and comparison with the experiments of Riboux et al. (Reference Riboux, Risso and Legendre2010) for $\alpha =0.54\,\%$, 1.1 %, 1.7 %, 3.0 %, 4.4 % and $8.2\,\%$.

For the same reason, the mean velocity of the bubbles decreases only very slightly with $\alpha$ according to the simulations, whereas experimentally, it is observed to decrease as $\langle v \rangle /v_0 \approx 0.6 \alpha ^{-0.1}$ (Riboux et al. Reference Riboux, Risso and Legendre2010). Also, the kinetic energies of the variances of the fluctuations of both the liquid and the bubbles are underestimated compared with the experiments.

Figure 3(a) compares the longitudinal spectra (that is to say $E_x(k_x)=\int 1/2 \phi _{xx}(\boldsymbol {k}') \delta (\boldsymbol {k}' \boldsymbol {\cdot } e_x -k_x)d^3\boldsymbol {k}'$ and $E_z(k_z)=\int 1/2 \phi _{zz}(\boldsymbol {k}') \delta (\boldsymbol {k}'\boldsymbol {\cdot } e_z-k_z)d^3\boldsymbol {k}'$ with $\phi _{ij}(\boldsymbol {k}) = \sum _{\boldsymbol {k}'}\langle \hat {u}_{i}(\boldsymbol {k}')\hat {u}_{i}^* (\boldsymbol {k}') \rangle \delta (\boldsymbol {k}-(\boldsymbol {k}'))$ and $\hat {\boldsymbol {u}}$ the coefficients of the Fourier series of the velocity field $\boldsymbol {u}_f(\boldsymbol {x})= \sum _{\boldsymbol {k}} \,{\rm e}^{{\rm i}\boldsymbol {k}\boldsymbol {\cdot } \boldsymbol {x}}\hat {\boldsymbol {u}}(\boldsymbol {k})$) of the vertical and horizontal velocities $E_z(k_z)$ and $E_x(k_x)$ obtained experimentally and by simulations. Note that, with this simulation approach, the spectra of the velocity field are very easy to obtain because the $\boldsymbol {u}_f$ field is smooth, whereas with a DNS-type approach, as well as in experiments, the velocity of the liquid is not defined everywhere, which poses a number of problems for spectral analysis. Here the approximations are made prior to the simulation, at the modelling phase, and there is no particular precaution to take for the calculation of the spectra. We can see that the spectra of the vertical and horizontal components are in fairly good agreement with the experiment. In particular, the simulations seem to reproduce a $k^{-3}$ evolution of the spectra as experimentally observed on small scales (large wavenumbers) and a $k^{-1}$ decay at large scales.

Figure 3. (a) Longitudinal velocity spectra of the vertical component in the vertical direction ($E_z(k_z)$) (continuous lines) and of the horizontal component in the horizontal direction (dashed lines) from the simulations at various $\alpha$ ($E_x(k_x)$). Comparison with the experiments of Riboux et al. (Reference Riboux, Risso and Legendre2010) for $\alpha =2.5\,\%$ and $d=2.5$ mm in blue, and with the power law $k^{-1}$ (grey dot-dashed line), $k^{-3}$ (grey dashed line) and $k^{-5/3}$ (grey dotted line). (b) Frequency spectra of the liquid velocity at the bubble position from the numerical simulation and comparison with the experimental spectra of the liquid velocity of the flow past a random array of fixed spheres (Amoura et al. Reference Amoura, Besnaci, Risso and Roig2017) in blue and with the $\omega ^{-1}$ and $\omega ^{-3}$ power laws.

However, we note, on the one hand, that the simulations underestimate the kinetic energy at small scales. We attribute this to the lack of near-bubble resolution, which leads to an underestimation of the power injected at the bubble scale (as discussed in the previous section). On the other hand, we also notice that the larger scales of the horizontal component are also underestimated, due to the absence of bubble trajectory oscillations, which are expected to enhance the redistribution of the fluctuating energy between the vertical and horizontal components.

It has also been reported that experimentally the spectra are invariant with the volume fraction and the bubble diameter. When the spectra of the numerical simulations are normalized by the injected power and by the viscosity, we observe the same invariance of the spectra of the numerical simulation with $\alpha$.

In figure 3(b), we present the frequency spectra of the vertical liquid velocity measured at the position of the bubbles (the true velocity, not the corrected one). These spectra are compared with the frequency spectrum of the velocity of the flow passing through an array of spheres held at a fixed position obtained experimentally by Amoura et al. (Reference Amoura, Besnaci, Risso and Roig2017) for a Reynolds number, based on the sphere diameter, of 600. Although the two flows are different, in both cases these spectra can be considered as characterizing the fluctuations of the liquid in the frame of reference moving with the bubbles. It can be seen that the simulations and the experiment show again a remarkable agreement. For frequencies lower than $d/v_0$ the spectrum shows an $\omega ^{-1}$ behaviour, while at high frequencies, the cut-off is much stronger with a slope close to $\omega ^{-3}$. It is interesting to note that this behaviour seems invariant with $\alpha$ in the simulations, while experiments have reported that it is invariant with the Reynolds number, provided that $Re>200$ (Amoura et al. Reference Amoura, Besnaci, Risso and Roig2017).

To summarize, while the bubble rising velocity and the total energy production are underestimated because of the filtering of the energy transferred from the bubble to the fluid, the interactions between wakes are well reproduced. Since this essential physical mechanism is correctly accounted for, the normalized spectra of the velocity are representative of real flows.

4. Characteristic scales

The spherically averaged spectra of the velocity (i.e. $E(k)=\int 1/2 \phi _{ii}(\boldsymbol {k}') \delta (|\boldsymbol {k}'|-k)d^3\boldsymbol {k}'$) are shown in figure 4. Contrary to the longitudinal spectra presented in figure 3(a), the three-dimensional spectra show a more complicated evolution with $k$ as well as a rather clear dependence on $\alpha$ at large scales. Several regions can be distinguished. The local maximum, located around $k \eta \approx 2$ in the figure, coincides, as we will see, with the scale of the bubbles which gives the cut-off scale of the energy injection. We see that for larger wavenumbers, a region in $k^{-3}$ clearly develops as $\alpha$ increases. The local minimum, located at large scales, corresponds to the wake scale. Between these two scales, the energy spectrum depends on $\alpha$ and corresponds to the scales directly influenced by the wakes. The fact that the spherically averaged spectra $E$ show such a qualitative difference at large scales with the spectra averaged over planes obviously indicates that the flow has a strong anisotropy at these scales. We will come back to the characterization of the anisotropy below.

Figure 4. Spherically averaged spectra of the velocity field $E$ from the numerical simulation at various $\alpha$. The $k^{-3}$ power law is shown by a grey dashed line. (The grey dotted line corresponding to a $k^{2/3}$ power law is just here as a guide for the eyes).

The flow being presumably dominated by the wakes of the bubbles and their interactions, an essential scale of the flow is the characteristic length of the wakes. To determine the latter, we consider the mean field conditioned on the position of a bubble (equivalent to a spatial phase average):

(4.1)\begin{equation} \langle \boldsymbol{u}_f \rangle_b(\boldsymbol{x}) = \dfrac{1}{T}\int \,{\rm d} t \dfrac{1}{N_b}\sum_{b=1}^{N_b} \boldsymbol{u}_f(\boldsymbol{x}-\boldsymbol{x}_b(t),t). \end{equation}

This mean field is illustrated in figure 5 for the case $\alpha =5\,\%$. This figure also shows the evolution of the mean vertical velocity along the vertical axis passing through the bubble for the different volume fractions, as well as for an isolated bubble. The first observation is that the wakes are much shorter in the case of the bubble swarm than in the case of an isolated bubble. By plotting the logarithmic derivative of the wakes, also shown in figure 5, we see that the wakes present a self-similar evolution for all $\alpha$ and that the mean velocity decreases exponentially with $z$. This remarkable feature is in agreement with the experimental results presented by Risso et al. (Reference Risso, Roig, Amoura, Riboux and Billet2008). This exponential decay of the wakes is likely due to the cancelling of the vorticity between neighbouring wakes, as proposed by Hunt & Eames (Reference Hunt and Eames2002).

Figure 5. (a) Cross-section of the vertical velocity conditionally averaged to the bubble position (equation (4.1)): $\langle u_{f,z} \rangle _b$ for $\alpha =5\,\%$. (b) Evolution of $\langle u_{f,z} \rangle _b$ along the vertical axis passing through the origin for the various $\alpha$ and for an isolated bubble (in grey dashed line). Comparison with the exponential decay with a rate $L_w(\alpha )$ (in grey dotted line). For ease of reading, we take advantage of the periodicity of the domain, and move the increasing part of $\langle u_{f,z} \rangle _b(z)$ upstream of the bubble. (c) Inverse of the logarithmic derivative, $({1}/{\langle u_z \rangle _b })({{\rm d} \langle u_z \rangle _b }/{{\rm d} z})$, of the former quantity, normalized by the characteristic wake length $L_w$. For this panel the vertical position is shifted by $z_0$ corresponding to the vertical position of the maximum of $\langle u_{f,z} \rangle _b$.

We choose the relaxation length of the exponential as the characteristic scale of the wakes $L_w$. The evolution of the ratio $L_w/d$ as a function of $\alpha$ is presented in figure 6. We can see that the length of the wakes shows evolution in $L_w \sim d \alpha ^{-1/3}$. One can interpret this evolution as a simple geometrical relation, considering that the wakes tend to screen each other. It should be noted, however, that this is quite a notable difference compared with the experiments of Risso et al. (Reference Risso, Roig, Amoura, Riboux and Billet2008) where the characteristic length of the wakes was observed to be independent of $\alpha$. This certainly reflects that there is an additional dependence of $L_w$ on $C_D$, since in the experiments the average speed of the bubbles decreases with $\alpha$.

Figure 6. Evolution of the characteristic scales $L_w$, $\lambda$, $\eta$ and $L_{int}$ in the simulations of the bubble swarm with $\alpha$. Comparison with the power law $\alpha ^{-1/3}$ in dashed lines, $\alpha ^{-1/6}$ in dotted lines and $\alpha ^{-2/3}$ in dot-dashed lines.

From this characteristic wake length, we define an inverse time scale $f=v_0/L_w$. This frequency $f$ can be considered as imposing a shear-rate scale to small scales $k\gtrsim 1/d$. This assumption allows us to estimate the average dissipation rate in the simulations as

(4.2)\begin{equation} \langle \varepsilon \rangle = \nu f^2.\end{equation}

Equivalently, we can interpret $L_w$ as a Taylor length scale based on the velocity $v_0$, $\lambda = \sqrt {\nu v_0^2/\langle \varepsilon \rangle }$. This is confirmed in figure 6 which shows that the evolution of $\lambda /d$ varies as $L_w/d$ in $\alpha ^{-1/3}$.

The volume-averaged power injected in the system corresponds to $\langle P_{tot} \rangle =n_b\langle P_b\rangle$ with $n_b$ the average number of bubbles per unit of volume and $P_b=\boldsymbol {F}_{f\rightarrow b} \boldsymbol {\cdot } \boldsymbol {v}_b$ the power supplied by bubble $b$. In the steady regime, this quantity is approximately given by $\langle P_{tot} \rangle = \alpha g v_0$ and as we have discussed above, it is larger than the power effectively received by the fluid in our simulations: $\langle P_{tot} \rangle >\langle P_f \rangle = \langle \varepsilon \rangle$. We will thus interpret $\langle P_f \rangle$ as the mechanical energy effectively injected in the wakes. Combining the previous relations, we find $\langle P_{tot} \rangle / \langle \varepsilon \rangle \sim \alpha Re_0 C_D (L_w/d)^2$ with $C_D = 4 gd/3v_0^2$. Therefore at $C_D$ and $Re_0$ constant, the proportion of energy injected in the wakes decreases as $\alpha ^{-1/3}$.

From the estimate (4.2) of the mean dissipation rate, we compute the dissipative scale $\eta = \nu ^{3/4} \langle \varepsilon \rangle ^{-1/4}= \sqrt {\nu /f} \sim d Re_0^{-1/2}\alpha ^{-1/6}$. It is this scale that is used to normalize the spectra presented in figures 3 and 4.

For $\alpha \geq 5\,\%$ we observed (not shown here, for brevity) that the kinetic energy of the liquid is invariant with $\alpha$ and is commensurate with $v_0^2$. Consequently we estimate the integral scale $L_{int} = \langle K\rangle ^{3/2}/ \langle \varepsilon \rangle$ to vary as $L_{int} \sim d Re_0 \alpha ^{-2/3}$. This behaviour is observed for $\alpha \geq 5\,\%$ in figure 6. Note that, in the experiments of Riboux et al. (Reference Riboux, Risso and Legendre2010), the liquid kinetic energy was observed to vary roughly as $\alpha v_0^2$. The discrepancy of the numerical simulations with the experiments is once again attributed to the absence of the fluctuations localized in the vicinity of the bubbles, which scale with $\alpha$.

5. Spectral analysis of bubble-induced turbulence

In order to identify the different regions of the spectra and to explain the observed scaling laws, we are interested in the spectral decomposition of the energy balance:

(5.1)\begin{equation} \frac{{\rm d}}{{\rm d} t} E(k) = T(k) -D(k) + P(k).\end{equation}

The terms of the right-hand side correspond, respectively, to the inter-scale energy transfer from a scale $k$ ($T$), the kinetic energy dissipation at a scale $k$ ($D$) and the rate of energy injected by the bubbles ($P$). The expressions of these different terms are obtained from the Navier–Stokes equation (2.1). The transfer term $T$ is the contribution of the nonlinear terms: $T(k) = \int \sum _{k'} [ -{\rm i} k'_j ( \widehat {u_iu_j}\hat {u}_i^*)] \delta (\boldsymbol {k'}-\boldsymbol {k})\delta (|\boldsymbol {k}|-k) d^3\boldsymbol {k} + {\rm c.c.}$, where $+ {\rm c.c.}$ denotes the complex conjugate terms, $D(k)= 2 \nu {k}^2 E(k)$, $P(k)$ is the integral over the wavenumbers $|\boldsymbol {k}|=k$ of the real part of $\hat {f}_i \hat {u}_i^*$ and the ‘hat’ denotes the Fourier transform. At steady state, the left-hand term of (5.1) is zero, so $T= D-P$. We present in figure 7 the terms $P(k)$ and $D(k)$ for the various $\alpha$. We can see that the production term presents a cut-off for $k>1/\sigma$ (we recall that $\sigma /d = 0.28$), and that on large scales it grows as $k^2$ for the largest $\alpha$, while it is roughly constant for small $\alpha$. Concerning the dissipation term, we notice that it also presents a peak around $k\sim 1/\sigma$. At large scales, the production dominates compared with the dissipation, which implies that $P(k) \approx -T(k)$. On the other hand, the dissipation dominates on small scales which means that $D(k)\approx T(k)$. The absence of scale separation between production and dissipation peaks means that this flow does not present an inertial zone. These budgets also show that there is no range of scales in which there is an equilibrium between $P$ and $D$. This contradicts the hypothesis made by Lance & Bataille (Reference Lance and Bataille1991) to explain the presence of a $k^{-3}$ zone in the velocity spectra. Furthermore, we notice that the $k^{-1}$ region of $D$, which corresponds to the $k^{-3}$ region of the velocity spectra, is observed in the crossover between the production-dominated scales and the dissipation-dominated scales.

Figure 7. Evolution of $P$ in continuous line and of $-D$ in dashed lines for the various $\alpha$. Comparison with $\alpha g v_0 \sqrt {2/{\rm \pi} } (k\sigma )^2\exp (-k^2\sigma ^2)$ in dotted line. (a) The power density and wavenumber are normalized by $\alpha g v_0d$ and $\sigma$, respectively. (b) The power density and wavenumber are normalized by the dissipative scales $\varepsilon \eta$ and $\eta$, respectively.

To interpret the behaviour of the production term $P$, we study the spectrum of the force applied to the flow, $E_f(k)$, corresponding to the spherical integration of $\hat {f}_i\hat {f}_i^*$. From the expression of the coupling force between the phases (2.2), we can obtain the following analytical expression for $E_f$:

(5.2)\begin{equation} E_f(k) / \alpha \sigma g^2 = \frac{1}{12{\rm \pi}}(d/\sigma)^3 (k\sigma)^2\, {\rm e}^{- (k\sigma)^2} . \end{equation}

This expression is obtained by assuming that (i) the positions of the bubbles are independent from each other and that (ii) the fluctuations of the rate of momentum exchanged between the bubble and the liquid are small: $\langle \boldsymbol {F}_{f\rightarrow b}^2 \rangle \approx \langle \boldsymbol {F}_{f\rightarrow b} \rangle ^2 = (\rho g {\rm \pi}d^3/6 )^2$. The spectra of the force are presented in figure 8. We note that, at all volume fractions, the agreement with the proposed expression is relatively good. We distinguish two regions: a region that grows as $k^2$ which reflects the equipartition of the fluctuations of the forces at large scales ($k<1/\sigma$) and an exponential decrease imposed by the Gaussian projection kernel $G_\sigma$ for $k>1/\sigma$. Note that the oscillations observed at the end of the spectra are due to the sharp cut-off of the kernel $G_\sigma (\boldsymbol {x}-\boldsymbol {x_b})$ for $|\boldsymbol {x}-\boldsymbol {x_b}|> 3 \sigma$. We notice that for high volume fractions, the positions of the bubbles are no longer really independent, because they cannot overlap, which explains that the prefactor of the $k^2$ increase at small $k$ is reduced compared with (5.2).

Figure 8. Spectra $E_f$ of the force at the various $\alpha$ and comparisons with (5.2) in dotted line and with the power law $k^2$ in dashed line.

It is more difficult to propose an analytical estimate of the spectrum $P$ of the work of $\boldsymbol {F}_{f\rightarrow b}$. However, the spectrum $E$ of $u$ is much less steep than the spectrum $E_f$ of $\boldsymbol {F}_{f\rightarrow b}$. At large scales, $E$ is rather flat compared with the $k^2$ evolution of $E_f$ and, at small scales, $E$ shows a power-law decay compared with the exponential cut-off $E_f$. It is thus reasonable to expect $P$ to behave similarly to $E_f$. Considering that $P$ is dominated by buoyancy, we obtain the following estimate:

(5.3)\begin{equation} P/ \alpha g v_0 \sigma \sim (k\sigma)^2 \,{\rm e}^{- (k\sigma)^2}. \end{equation}

Indeed, $P$ presents an exponential damping for $k>1/\sigma$, that overlaps for all $\alpha$ when normalized by $\alpha g v_0 \sigma$ as can be seen in figure 7(a). For $k<1/\sigma$ and $\alpha >2\,\%$, $P$ increases roughly as $k^{2}$ in agreement with the previous relation. The underestimation of $P$ at large scales is attributed to the the large fluctuations of the injected power $\langle P_f^2 \rangle \gg \langle P_f \rangle ^2$. For small $\alpha$, $P$ no longer follows the proposed relation on large scales due to the presence of large structures in the flow leading to a significant correlation of the liquid velocity between distant bubbles.

For completeness, we present in figure 7(b) the production and dissipation terms normalized by the dissipative scales $\langle \varepsilon \rangle$ and $\eta$ for the various $\alpha$. Consistently with the velocity spectra shown previously in figure 4, with this normalization, we observed that the values of $D$ of all $\alpha$ overlap at high wavenumbers (typically $k>1/2\eta$). From the estimates of the characteristic scales of the simulations presented in the previous section, we have $d/\eta = Re_0^{1/2}\alpha ^{1/6}$, indicating that the gap between the production-dominated scales and the dissipation-dominated scales increases, slowly, with $\alpha$. It seems that the $k^{-1}$ subrange of $D$ is observed in this gap of scales, provided that $\alpha$ is large enough.

In conclusion, we consider that for $k < 1/d$ the flow structure is driven by the interactions between wakes, while in the range $1/d< k < 1/\eta$ the strong damping of the wakes imposes a shear scale.

This assumption of constant shear rate $f$ across scales allows us to explain the presence of a power law in $k^{-3}$ for the flow, because of a matching argument similar to that proposed by Kolmogorov in 1941. It is assumed that at scales that are small compared with the length of the wakes ($k \gg 1 /L_w$) the structure of the flow depends only on the diameter of the bubbles $d$, the viscosity and the shear rate $f$ :

(5.4)\begin{equation} E=E(k; d, \nu, f ). \end{equation}

It is assumed that at scales larger than $\eta = \sqrt {\nu /f}$, one can neglect the effect of viscosity. Therefore, in this limit we can write

(5.5)\begin{equation} E=E_I(k; d,f ) = d^3 f^2 \varPhi_I(kd), \end{equation}

where $\varPhi _I$ is a dimensionless function. Conversely, at scales much smaller than the bubble size, we will suppose that the diameter no longer plays a role, and we will make the hypothesis that

(5.6)\begin{equation} E=E_S(k; \nu, f) = \nu^{3/2} f^{1/2} \varPhi_S(k\eta), \end{equation}

where $\varPhi _S$ is another dimensionless function. Finally, if we assume that for a range of intermediate scales ($1/d \ll k \ll 1/\eta$) the two previous relations remain valid, we have $E_S(k)=E_I(k)$. Since $kd$ and $k\eta$ can vary independently, the previous equality can hold only if the following expressions are constant:

(5.7)\begin{equation} (kd)^3 \varPhi_I(kd) = (k\eta)^3 \varPhi_S(k\eta) = c. \end{equation}

This gives us for the velocity spectra:

(5.8)\begin{equation} E(k)=c f^2 k^{{-}3},\end{equation}

in a range of scales where the shear rate can be considered constant.

The temporal spectra of the velocity seen by the bubbles (presented in figure 3b) are influenced by the fact that the bubbles cut the wakes of other bubbles. Thus the high frequencies of the temporal spectrum are dominated by the Doppler shift due to the high-speed crossing, of the order of $v_0$, of the dissipative structures of the flow. So using (5.8) and taking the argument of Tennekes (Reference Tennekes1975), with $\omega \sim v_0 k$, we can estimate the high-frequency behaviour of the frequency spectra:

(5.9)\begin{equation} E(\omega) = E(k) \dfrac{k}{\omega} = f^2v_0^2 \omega^{{-}3}. \end{equation}

The temporal spectra from both the experiments and the simulations present an $\omega ^{-3}$ zone at high frequency.

We have seen that at large scales ($k\ll 1/d$), where the flow is dominated by wake interactions, there is a balance between production and inertia. At these scales we notice that the one-dimensional spectra of the velocity present a $k^{-1}$ dependence. This means that each decade contains an equal amount of energy. This behaviour can be explained by the intermittence of the wake passages, giving rise to an alternation between periods of activity and calm (Mandelbrot et al. Reference Mandelbrot, Berger, Kahane and Peyriere1999). At these scales, the characteristic velocity no longer depends on a specific length scale and corresponds to the typical velocity of the bubbles $v_0$. These periods of activity (the wakes) are characterized by their self-similar character (Marinari et al. Reference Marinari, Parisi, Ruelle and Windey1983; Bak, Tang & Wiesenfeld Reference Bak, Tang and Wiesenfeld1987) and present a variable intensity and duration, whereas the calm periods follow a Poissonian distribution reflecting the quasi-uniform distribution of the bubbles. The absence of a characteristic length leads directly to the absence of a characteristic time for the fluctuations. Hence the frequency spectrum of the velocity also shows a decay close to $\omega ^{-1}$ at low frequency. As pointed out by Mandelbrot (Reference Mandelbrot1967), these behaviours in $\omega ^{-1}$ and $k^{-1}$ must also be connected to the non-Gaussianity of the velocity distributions as well as to a long-range correlation of the velocity.

6. Characterization of the anisotropy

To extend the discussion on the large scales, we need to take into account the anisotropy of the flow. For the characterization of the anisotropy, it is necessary to distinguish the fact that the energy can be carried mainly by one component of the velocity vector (anisotropy between components) from the fact that the fluctuations in certain directions can carry more energy (directional anisotropy) (Sagaut & Cambon Reference Sagaut and Cambon2008). To characterize the latter, we consider a spherical coordinate system of the wave vector space, as schematized in figure 9. The angle $\theta$ characterizes the orientation of the wave vector with respect to the vertical direction: $\sin \theta = k_z / |k|$ ($\sin \theta =0$ corresponds to fluctuations in the horizontal direction and $\sin \theta =\pm 1$ to fluctuations in the vertical direction). We then consider the directional spectra $E(k,\theta )$ which allows decomposing the energy of the fluctuations according to the wavelength and the orientation with respect to the vertical direction. More specifically, $E(k, \theta )$ is defined by integration on all ‘longitudes’ for a fixed ‘latitude’ and modulus of the wave vector: $E(k, \theta ) = \int 1/2 \phi _{ii}(\boldsymbol {k}') \delta (|\boldsymbol {k}'|-k)\delta (k'_z/|\boldsymbol {k}'|-\sin \theta )d^3\boldsymbol {k}' / 2{\rm \pi} \cos \theta = \int _0^{2{\rm \pi} }1/2 \phi _{ii}(k,\theta,\phi )k^2\sin \theta d \phi / 2{\rm \pi} \cos \theta$. The normalization factor $2{\rm \pi} \cos \theta$ is introduced to correct the geometrical effect due to the fact that a band near the poles covers a less important surface than a band near the equator.

Figure 9. (a) Scheme of spherical coordinates in Fourier space. Directional spectra of (b) the velocity vector, (c) the vertical velocity component and (d) a horizontal velocity component for $\alpha =10\,\%$. Note that the spherical spectra correspond to an integral along a dashed grey line, while the longitudinal spectra in the vertical direction are obtained by integration along a dashed black line.

We show in figure 9 the directional spectra of the velocity for $\alpha =10\,\%$. It can be seen that at large scales ($k<1/d$) the energy is concentrated in the horizontal direction. This concentration is characteristic of vertically aligned tubular structures (Cambon & Jacquin Reference Cambon and Jacquin1989), which can be seen on the flow visualization in figure 1. At smaller scales ($k \gg 1/d$), one can see that the directional spectra become invariant with $\theta$ indicating that the directional anisotropy tends to vanish.

To characterize the anisotropy between components, we present as well in figure 9 the directional spectrum for both the vertical and the horizontal components of the velocity. Note that the axisymmetry of the flow imposes $E(k,\theta ) = E_z(k,\theta )+2 E_x(k,\theta )$. At large scales, we notice that the spectrum of the vertical velocity is very similar to that of the total kinetic energy, which indicates that at these scales the vertical component carries almost all the kinetic energy. This can be explained simply by the fact that the forcing due to the bubbles is essentially vertical. We can also note that the horizontal component of the velocity presents a very weak directional anisotropy. Finally, at small scales, we notice that the flow tends to become much more isotropic and presents a decrease of the difference both between the components and between the directions. This indicates that at scales where the energy injection is zero, there is a redistribution between the components, which ensures a return to isotropy as $k$ increases.

It should be noted that when considering the spectra with angular dependence for both the horizontal and vertical components, no $k^{-3}$ zone is distinguished. Thus, as noticed by Bellet et al. (Reference Bellet, Godeferd, Scott and Cambon2006) for the decaying turbulence under strong rotation, the $k^{-3}$ region for the spherically averaged spectra results from an average between the different directions and the different components.

If the structure of the flow becomes locally isotropic at scales much smaller than $d$, one should expect the appearance of a $k^{-5/3}$ inertial range for the velocity spectra, if the local Reynolds number is large enough, since at these scales there is no more energy injection. Such inertial range is not present in the numerical simulations reported here, but seems visible in the spectra obtained from the experiments of Riboux et al. (Reference Riboux, Risso and Legendre2010) (see figure 3a). In the case where the velocity spectra first present evolution as $E(k)\sim f^2k^{-3}$ followed by an inertial range $E(k)\sim \langle \varepsilon \rangle ^{2/3} k^{-5/3}$, the characteristic length of the crossover between these two regimes would be given by $\ell _I = \sqrt {\langle \varepsilon \rangle /f^3}$. It is interesting to note that this length scale corresponds to the classical estimates of the scale from which a turbulent flow subject to mean shear can be considered as locally isotropic (Champagne, Harris & Corrsin Reference Champagne, Harris and Corrsin1970; Pope Reference Pope2000). Taking the usual Kolmogorov scale $\eta =\nu ^{3/4}\langle \varepsilon \rangle ^{3/4}$, the extension of the inertial regime is given by $\ell _I/\eta =(\langle \varepsilon \rangle /\nu f^2)^{3/4}$. In the simulations, as mentioned previously, no inertial range is present ($\ell _I=\eta$) and we have indeed $\langle \varepsilon \rangle = \nu f^2$. It is likely that a larger power injection in the simulations (e.g. increasing $C_D$) would allow obtaining a separation between the scales of return to isotropy and the dissipative scales and thus to obtain a $-5/3$ range in agreement with the experiments. Nevertheless, even in the absence of the $k^{-5/3}$ regime, this indicates that the rate of shear imposed by the bubble wakes controls the relaxation to small-scale isotropy.