2. Mathematical formulation

2.1. Problem statement

We investigate the dynamics of a monodisperse suspension of bubbles rising under the action of buoyancy in a fluid initially at rest. Several physical parameters characterise the problem: the gas volume fraction $\phi$; the number of bubbles $N_b$; the diameter of the bubbles $d_b$ calculated as the diameter of the sphere of equivalent volume; the gravity acceleration $g$; the viscosity of the two fluids $\mu _b,~\mu _l$; their densities $\rho _b,~\rho _l$; and the surface tension $\sigma$. We use the subscripts $b$ for bubbles and $l$ for liquid. The density, the viscosity and the surface tension of each fluid are considered constant during each numerical experiment. Four dimensionless groups can be formed in addition to the number of bubbles and the volume fraction. Two are the density and viscosity ratio, ${\rho _b}/{\rho _l}$ and ${\mu _b}/{\mu _l}$, respectively. We briefly analyse the impact of the density ratio but in almost all simulations we have fixed ${\rho _b}/{\rho _l}=10^{-3}$ and ${\mu _b}/{\mu _l}=10^{-2}$, which are typical values for air bubbles in water.

The other two dimensionless groups can be characterised by the Galileo number

(2.1)\begin{equation} Ga\equiv\frac{\rho_l\vert {\rm \Delta} \rho\vert g d_b^3}{\mu_l^2}, \end{equation}

where ${\rm \Delta} \rho = \rho _b-\rho _l$, or equivalently the Archimedes number $Ar\equiv \sqrt {Ga}$ and the Eötvös (or Bond) number

(2.2)\begin{equation} Eo\equiv \frac{\vert {\rm \Delta} \rho\vert g d_b^2}{\sigma}. \end{equation}

These numbers indicate the relative importance of buoyancy and surface tension.

When bubbles move, the flow is also characterised by a velocity scale, which is typically given by the average bubble velocity $\langle U_b \rangle$, which we have computed averaging in space. It is then possible to define the bubble Reynolds number based on this velocity

(2.3)\begin{equation} Re\equiv\frac{\langle U_b \rangle d_b}{\nu_l}, \end{equation}

where $\nu _l$ is the kinematic viscosity of the liquid. It is also possible to use another group which compares inertial effects with surface tension, the Weber number

(2.4)\begin{equation} We \equiv \frac{\rho_l \langle U_b \rangle ^2 d_b}{\sigma} = \frac{Eo Re^2}{Ga}. \end{equation}

It is important to note that the average bubble velocity may or may not reach a stationary state in our numerical experiments, so that in general the dynamic dimensionless numbers are dependent on time $Re=Re(t)$.

2.2. Governing equations

Both fluids are governed by the Navier–Stokes equations, which we take here in the incompressible limit

(2.5)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}

(2.6)\begin{gather}\frac{\partial {\boldsymbol u}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} ({\boldsymbol u} \otimes {\boldsymbol u}) = \frac{1}{\rho} (-\boldsymbol{\nabla} p + \boldsymbol{\nabla} \boldsymbol{\cdot} (2 \mu \boldsymbol{ D}) + {\boldsymbol f} + {\boldsymbol f}_{\sigma}\delta_s), \end{gather}

here the viscosity $\mu$ and the density $\rho$ varies across the two phases; $\boldsymbol { D}= [\boldsymbol {\nabla } u + (\boldsymbol {\nabla } u)^{\textbf {T}}]/2$ is the symmetric deformation tensor; ${\boldsymbol f}$ represents the volumetric forces, which in the present case are the gravity $\boldsymbol {f}_i=\rho {\boldsymbol g}$; ${\boldsymbol f}_{\sigma }$ is the force exerted by the surface tension; $\delta _s = \delta _s ({\boldsymbol x} - {\boldsymbol x}_s)$ is a Dirac delta function that identifies the presence of the surface. The volumetric surface tension force is expressed as (Tryggvason et al. Reference Tryggvason, Scardovelli and Zaleski2011)

(2.7)\begin{equation} {\boldsymbol f}_{\sigma} = \sigma \kappa {\boldsymbol n} + \nabla_s \sigma. \end{equation}

The first term depends on the surface tension coefficient (a material property), the local curvature $\kappa = \boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol n}$ and the surface normal, while the last term is different from zero only if a non-constant surface tension is present. In the present work, we shall deal with constant surface tension and, therefore, the second term is zero. In practice the surface tension balances the jump in pressure across the interface and jump relations can be derived analogously to shock waves. It is worth remarking that since the surface force acts in the plane of the surface, if we integrate it over the whole closed surface it should give a null contribution. More details about the numerical representation of the surface tension are given in a recent review (Popinet Reference Popinet2018).

This set of equations is solved with the Basilisk library with the numerical methods described in the following section.

3. Numerical method

Basilisk is a library of solvers written using an extension of the C programming language, called Basilisk C, adapted for discretisation schemes on Cartesian grids (see http://basilisk.fr). Space is discretised using a Cartesian (multilevel or tree-based) grid where the variables are located at the centre of each control volume (a square in 2-D, a cube in 3-D) and at the centre of each control surface. The possibility to adapt the grid dynamically is key to efficiently simulate multiphase flows (Popinet Reference Popinet2009). Two primary criteria are used to decide where to refine the mesh. They are based on a wavelet decomposition of the velocity and volume fraction fields, respectively (van Hooft et al. Reference van Hooft, Popinet, van Heerwaarden, van der Linden, de Roode and van de Wiel2018). The velocity criterion is mostly sensitive to the second-derivative of the velocity field and guarantees refinement in developing boundary layers and wakes. The volume fraction criterion is sensitive to the curvature of the interface and guarantees the accurate description of the shape of bubbles. Both criteria are usually combined with a maximum allowed level of refinement. As demonstrated in previous work, using the earlier code Gerris (Cano-Lozano et al. Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016a), this strategy leads to very large savings in computational cost compared with fixed Cartesian grid approaches.

The numerical scheme implemented in Basilisk is very close to that used in Gerris as described in Popinet (Reference Popinet2009). The Navier–Stokes equations are integrated by a projection method (Chorin Reference Chorin1969). Standard second-order numerical schemes for the spatial gradients are used (Popinet Reference Popinet2003, Reference Popinet2009; Lagrée, Staron & Popinet Reference Lagrée, Staron and Popinet2011). In particular, the velocity advection term $\partial _j (u_j u_i )^{n+1/2}$ is estimated by means of the Bell–Colella–Glaz second/third-order unsplit upwind scheme (Popinet Reference Popinet2003). In this way, the problem is reduced to the solution of a 3-D Helmholtz–Poisson problem for each primitive variable and a Poisson problem for the pressure correction terms. Both the Helmholtz–Poisson and Poisson problems are solved using an efficient multilevel solver (Popinet Reference Popinet2003, Reference Popinet2015).

Time is advanced using a second-order fractional-step method with a staggered discretisation in time of the velocity and scalar fields (Popinet Reference Popinet2009): one supposes the velocity field to be known at time $n$ and the scalar fields (pressure, temperature, density) to be known at time $n-1/2$, and one computes velocity at time $n+1$ and scalars at time $n+1/2$.

The interface between the fluids is tracked with a geometric VOF method (Hirt & Nichols Reference Hirt and Nichols1981; Scardovelli & Zaleski Reference Scardovelli and Zaleski1999; Tryggvason et al. Reference Tryggvason, Scardovelli and Zaleski2011). The surface tension term is computed using an accurate well-balanced, height-function method (Popinet Reference Popinet2018). In this formulation, the surface tension in (2.6) is expressed as a gradient, and may thus be included in the pressure term.

Periodic, no-slip and free-slip boundary conditions will be imposed in the different computations considered.

In the present work, we always consider flows with a bubble concentration of a few per cent, $\phi <5\,\%$. It is known that in this case, coalescence and breakup effects are negligible (Jha & Govardhan Reference Jha and Govardhan2015). We have checked that the resolution and the physical set-up are always consistent to avoid spurious effects, as briefly described in Appendix B.

4. Preliminary tests

To assess the accuracy of the numerical code for the simulation of two-phase bubbly flows, we have reproduced several literature test cases (Sangani Reference Sangani1987; Esmaeeli & Tryggvason Reference Esmaeeli and Tryggvason1998, Reference Esmaeeli and Tryggvason1999; Loisy, Naso & Spelt Reference Loisy, Naso and Spelt2017). In particular, we have focused on the configuration of regular arrays of bubbles rising due to buoyancy. Previous numerical studies were carried out using different approaches, namely front tracking (Esmaeeli & Tryggvason Reference Esmaeeli and Tryggvason1998) and level-set with diffuse interface (Loisy et al. Reference Loisy, Naso and Spelt2017). The present comparison thus allows a mutual validation of the different methods. After an initial transient where bubbles accelerate, they eventually reach a quasi-steady-state regime. Depending on bubble size, surface tension and density, they may follow non-rectilinear paths, with periodic or chaotic lateral oscillations (Cano-Lozano et al. Reference Cano-Lozano, Martinez-Bazan, Magnaudet and Tchoufag2016a). A regular array of bubbles is reproduced numerically using a single bubble in a periodic cell. Changing the cell size with respect to the bubble size, we can adjust the volume fraction of the array. Note that since the computational domain is unbounded in all directions, an additional body force $- \langle \rho \rangle g$ must be added to avoid the system accelerating in the vertical downward direction. In this section, we present briefly only the most significant results, while more details are given in Appendix A.

We have first compared our simulations with the theory of Sangani (Reference Sangani1987) for the Stokes flow regime. The configuration consists of a cubic array of spherical bubbles at different volume fractions. The non-dimensional numbers of the simulation are the same as in previous DNS studies, namely

(4.1a–d)\begin{equation} Ar = 0.15, \quad Eo= 0.38, \quad \rho_b / \rho_l = 0.005, \quad \mu_b / \mu_l = 0.01. \end{equation}

Although at very low Reynolds number, this is a severe test case since it is 3-D and the number of points required may increase rapidly when varying the concentration. We have carried out simulations at different resolution, asking for a relative adaptation error less than $5\,\%$. In table 1 we show the steady-state velocity of the bubble array normalised with the velocity of a single isolated bubble, and the quantitative numerical difference from the analytical solution. A satisfactory agreement is obtained between the numerical and the analytical solution ${U}/{U_0} = 1 - 1.1734 \mu ^* \phi ^{1/3} + O(\phi)$, where $\mu ^*=(\mu _l+3/2\mu _b)/(\mu _l+\mu _b)$, $U$ is the drift velocity and $U_0$ is the terminal velocity of a single bubble. More specifically, $U$ is the vertical component of ${\boldsymbol U} = \langle \boldsymbol u \rangle _b - \langle \boldsymbol u \rangle$, where $\langle \rangle$ means an average over the entire cell, while $\langle \rangle _b$ denotes the average over the volume occupied by the bubble only.

Table 1. Direct numerical simulations and analytical results for the test analysed by Sangani (Reference Sangani1987). The volume of the cell is kept always the same, and in the last column we display the resolution used.

We have then considered test cases at finite Reynolds numbers. In figure 1(a), we show the results for the 3-D test case proposed by Esmaeeli & Tryggvason (Reference Esmaeeli and Tryggvason1999). In this case the flow parameters are

(4.2a–d)\begin{equation} Ar = 29.9, \quad Eo= 2, \quad \rho_b / \rho_l = 0.1, \quad \mu_b / \mu_l = 0.1. \end{equation}

Our simulations are compared against both the original DNS of Esmaeeli & Tryggvason (Reference Esmaeeli and Tryggvason1999), and the more recent results of Loisy et al. (Reference Loisy, Naso and Spelt2017). We have analysed the grid convergence. Results are in good agreement, while convergence is achieved with a slightly larger number of points ($d_b/{\rm \Delta} x\approx 40$) than in previous works (Esmaeeli & Tryggvason Reference Esmaeeli and Tryggvason1999; Loisy et al. Reference Loisy, Naso and Spelt2017), where the authors indicate that $30$ points per diameter are sufficient.

Figure 1. (a) Time evolution of the bubble Reynolds number at different grid resolutions for the 3-D configuration proposed by Esmaeeli & Tryggvason (Reference Esmaeeli and Tryggvason1999). (b) Time evolution of the components of the bubble Reynolds number for the case of steady oblique rise, compared against the results by Loisy et al. (Reference Loisy, Naso and Spelt2017). The Reynolds number is defined here as $Re=U d_b/\nu$, recalling that $U$ is the vertical component of the drift velocity ${\boldsymbol U} = \langle \boldsymbol u \rangle _b -\langle \boldsymbol u \rangle$, where $\langle \rangle$ means an average over the entire cell, while $\langle \rangle _b$ denotes the average over the volume occupied by the bubble only.

The last set of moderate $Re$ test cases is the oblique rise of periodic arrays of bubbles performed by Loisy et al. (Reference Loisy, Naso and Spelt2017). For this test the numerical set-up is the same as for previous tests, i.e. a single periodic lattice to mimic a regular array. Loisy et al. (Reference Loisy, Naso and Spelt2017) pointed out that for certain values of the non-dimensional parameters, bubbles can have an oblique trajectory (not aligned with gravity) at certain volume fractions, although a single bubble in the same parameter regime would follow a straight vertical path. Analytical considerations support the possibility of a non-trivial path indicating a possible transition for $Ar\approx 20$. In particular three different oblique regimes have been found: (a) a steady oblique rise; (b) an oscillatory oblique rise, with a bubble oscillating around a straight oblique path; and (c) a chaotic oblique rise. Such a behaviour had been previously noticed numerically (Sankaranarayanan et al. Reference Sankaranarayanan, Shan, Kevrekidis and Sundaresan2002), but using a diffuse interface method and a small density ratio. In the present work, we have simulated the configurations corresponding to the three regimes in Loisy et al. (Reference Loisy, Naso and Spelt2017). The density ratio and the viscosity ratio between the two phases are the same for all the cases, $\rho _b/ \rho _l = 0.005$, $\mu _b / \mu _l = 0.01$. The number of points is varied together with the domain size in order to always get the same bubble resolution $d_b / \varDelta = 40$. Global agreement between present simulations and those by Loisy et al. (Reference Loisy, Naso and Spelt2017) is excellent. In particular, final values are the same within numerical errors. In figure 1(b) we show as an example the comparison for regime a. The details of these simulations together with other results are given in Appendix A.

Before analysing complex flows at high Reynolds numbers, we have also carried out a specific quantitative analysis on the effect of two crucial numerical issues: (i) resolution; (ii) density ratio. It is worth emphasising that there is a strong link between physical properties and numerical parameters and that this cannot be overlooked. While the simulation of a single bubble remains feasible even with a very fine grid thanks to the adaptive mesh, it would not be possible to tackle a problem with many bubbles with the same grid. Moreover, without the adaptive mesh even the single bubble case appears desperate at large Reynolds numbers. In contrast, using a coarse grid may make the computation easy but the results might be largely unreliable. We summarise here our findings, details are given in Appendix B. In order to simulate bubble flows quantitatively and in detail, it turns out to be key: (1) to have a number of points per diameter increasing with the $Ar$ number (we have found that convergence is obtained with approximately $N_{points}\approx Ar/2$); (2) using an adaptive mesh, it is sufficient to have such a fine resolution inside the bubble and in the wake; (3) a large density ratio, namely $\rho _l/\rho _b > 100$, is mandatory to avoid spurious effects which are similar to those found with a too-coarse grid, leading to a too-high rate of coalescence. Our results confirm the previous results by Cano-Lozano et al. (Reference Cano-Lozano, Tchoufag, Magnaudet and Martínez-Bazán2016b).

5. Pseudo-turbulence in two dimensions

In this section, we discuss the results of a 2-D bubble column configuration (Biswas & Tryggvason Reference Biswas and Tryggvason2007). We consider a square domain with the vertical direction $z$ aligned with gravity, acting downward. The tank, of size $50 d_b\times 50 d_b$, is filled with a liquid and $32$ initially spherical bubbles are placed at the bottom, in a region confined between $z = 0$ and $z = 8 d_b$, and are homogeneously distributed in the lateral direction $x$, while avoiding any initial bubble overlap, and with a minimum distance between them of one diameter. This results in a local volume fraction in the region $0 \le z \le 8$ of $\phi \simeq 5\,\%$. The domain is closed at the bottom by a wall (no-slip boundary condition), and an outflow boundary condition is used at the top, while on the lateral sides the domain is periodic. At $t=0$ both the liquid and the bubbles are at rest.

The viscosity and density ratios are constant in all the simulations and their values are $\mu _l / \mu _b = 100$ and $\rho _l / \rho _b = 1000$. Three different simulations have been carried out, and the corresponding parameters are reported in table 2. In particular, the $Ar$ number is within the range $Ar \simeq 100\text {--}300$, which corresponds to typical 3-D experiments (Riboux et al. Reference Riboux, Risso and Legendre2010). For 2-D cases, in all the three cases we have used regularly spaced grids with different resolutions depending on the increasing bubble Reynolds number. In any case, the resolution requirements to get physically sound results have always been fulfilled, as highlighted in table 2.

Table 2. Non-dimensional parameters for the 2-D bubble column. Here $N$ represents the number of points and $d_b/\varDelta$ the grid resolution in terms of points per bubble diameter. The Reynolds number usually defined as $Re_b=\langle U_b \rangle d_b/\nu$ is not defined a priori. Since the present test case is not steady, it is not possible to identify it clearly. We have computed it by averaging over the time range where it is approximately steady ($t\in [13-20])$ to obtain $Re_b\approx 200, 280, 470$ for case a, case b, case c, respectively.

We have focused in this work on the liquid agitation induced by bubbles within the swarm. Yet, since the problem is non-homogeneous in the vertical direction and non-stationary, particular care must be taken in the procedure used to compute the observables and we have, therefore, performed a careful analysis in two dimensions to prepare the 3-D case. As in the experiment of Riboux et al. (Reference Riboux, Risso and Legendre2010), the spectra $S_{ii}(k) = \langle | \hat {u}_i(k) |^2 \rangle$, where $\hat {u}_i(k)$ is the Fourier transform of the velocity fluctuation in the $i$ direction, are evaluated separately for the vertical and the horizontal components. The transform is performed in the $x$ direction, which can be considered as statistically homogeneous, for both components of the velocity. We have in particular verified that $\langle U \rangle _x=0$. We have computed the statistics at each $z$ inside the region $z\in [15-25]$, and at different times when bubbles are inside this interrogation window. It is worth remarking that in the statistical analysis of the velocity fluctuations we have used all the grid data available in the window, therefore belonging to both phases, as done by Dodd & Ferrante (Reference Dodd and Ferrante2016) for droplets. The results have been found to be statistically homogeneous to a good degree of approximation (less than 5 %) over windows of length $5d_b$. For this reason, in the following we show results averaged over $5d_b$, to improve statistics. In particular, we show statistics only computed between $z=15$ and $z=20$, because results obtained in the second window are practically indistinguishable.

As detailed in Appendix C, we have found that the spectra are independent of time, over almost the whole time-window considered. In particular the spectral slope appears rather constant, when the bubbles have entered and not yet left the interrogation window. Furthermore, no appreciable difference is found between the horizontal and the vertical spectrum, showing that both components dynamically distribute the energy in a similar manner. We can then write that the one-dimensional spectrum is $E(k)=S_{ii}$ without compromise. We compare the spectra at different $Ar$ numbers, at the same time $t=15$, as shown in figure 2. Time is always made non-dimensional with the bubble buoyancy time $\sqrt {d_b/g}$. In all cases spectra are compatible with a scaling $E(k)\sim k^{-3}$ in a range around the diameter scale. At small scales, a steeper scaling $E(k)\sim k^{-4}$ is found also in all cases, which can be related to a range where viscous effects are important (Monin & Yaglom Reference Monin and Yaglom1975). However, for case a (as shown in table 2) this dissipative range appears to dominate over almost the whole range of scales smaller than the diameter. In case b, the spectrum displays a $-3$ slope over roughly a decade, while for the highest $Ar$ number the range appears even larger. Moreover, we observe for cases b and c that around the bubble diameter there is a crossover and the spectrum is flatter at larger scales with a slope close to $-5/3$. To check the statistical robustness of our analysis, we have repeated the simulation of the case at $Ar=313$ with periodic conditions in both directions. In this case, the flow is statistically homogeneous in all directions, and after a transient a steady-state is attained. Therefore, both spatial and time averages are taken. The periodic simulation confirms the results obtained in the unsteady case. In particular, a $k^{-5/3}$ scaling is obtained at scales larger than the bubble diameter. The $k^{-3}$ scaling appears to be present at scales smaller than the bubble diameter and then a steeper slope typical of a viscous range is found. The $k^{-5/3}$ suggests the presence of an inverse cascade, typical of 2-D turbulence (Boffetta & Ecke Reference Boffetta and Ecke2012), as confirmed by the negative kinetic-energy flux displayed in Appendix C. In figure 3 we show the vorticity field in the space window that has been used for the evaluation of the spectra at a fixed time $t = 15$. The visualisation allows us to link the statistical spectral properties to the actual dynamics of the flow. The bubbles are a source of vorticity, which then creates the trailing wakes. We observe that at $Ar=100$, the interaction between the wakes exists but is small, notably in the upper part of the window. The plot for $Ar=140$ clearly suggests a stronger interaction between bubble wakes, and the vorticity field is diffused through nonlinear interactions. The case at $Ar=313$ is similar to the $Ar=140$, but the strong interaction between wakes and the presence of dynamics at smaller scales are even more visible, with thin unstable vorticity filaments released behind the bubbles. The nonlinear wake interactions are clearly dominant here and bubbles follow quite intricate paths. Although a $k^{-3}$ scaling has been found in all cases, the present results show that in case a the spectrum is basically related to the coherent structures of the wakes. In contrast to the other two cases, because of the higher Reynolds number, the agitation induced by bubbles starts to play an important role. Notably bubble dynamics lead to an injection of energy and vorticity at the scale of the bubble diameter, and energy is transferred towards different scales. In both cases at $Ar=140$ and $Ar=313$ these interactions are significant enough that an inverse cascade of energy towards large scales could be triggered, as suggested by the $-5/3$ scaling of the spectrum. In figure 4, the p.d.f.s of the velocity fluctuations for the different cases are shown together with those obtained in the steady case. The velocity fluctuation field is computed at each $z$ subtracting the average velocity computed over the corresponding plane ${\boldsymbol u^\prime }={\boldsymbol U}-\langle \boldsymbol U \rangle _z$. We have verified that keeping only the liquid phase does not change appreciably the results. From a physical point of view, p.d.f.s are clearly not Gaussian with exponential tails, and while the horizontal one is symmetric, the vertical one is skewed, showing anisotropy of fluctuations and the particular status of the vertical direction. The p.d.f.s obtained are similar for all the $Ar$ studied, although it has been observed that the dynamics is different. In particular, it has been observed that stronger interactions at higher $Ar$ lead to more intricate paths. While the vertical p.d.f. appears a little less skewed at higher $Ar$, the difference is within statistical error. It is worth noting nonetheless that this p.d.f. is a global one-point statistical observable, and the link between it and instantaneous geometrical differences is not straightforward.

Figure 2. (a) Spectra of the vertical component of the velocity of liquid fluctuations for different $Ar$ evaluated at time $t=15$. The vertical line corresponds to the bubble diameter. The dot–dashed line indicates the $k^{-5/3}$ slope and the dashed line the $k^{-3}$ slope. (b) Energy spectrum of vertical fluctuations against $k$ for the case $Ar=313$ for both the unsteady and the steady configurations. For the steady case, the spectrum is obtained by averaging over time between $t=13$ and $t=23$. Lines are the same as in panel (a).

Figure 3. Vorticity field displayed in the domain between $15$ and $25$ bubble diameters in the vertical direction, at time $t=15$ for the different $Ar$ cases: (a) $Ar = 100$ and $Eo= 0.12$; (b) $Ar = 140$ and $Eo= 0.2$; (c) $Ar = 313$ and $Eo= 0.56$. The colour bar is the same for the three cases and is displayed laterally.

Figure 4. Probability density functions of the velocity fluctuations in the lateral $x$ (a) and vertical $z$ (b) directions for different $Ar$. The steady simulation at $Ar=313$ is plotted for comparison. The unsteady p.d.f.s are shown at $t=14$ for the cases at $Ar=100,140$, and at $t=20$ for $Ar=313$. The time are chosen such that observables are computed well within the swarm. Yet, the results obtained at different times are very similar. As in figure 2, in the steady case, time averages have been taken in the window $t=13\text {--}23$.

From a statistical point of view, the p.d.f.s show unambiguously that results obtained in the unsteady regime are statistically robust, provided the analysis is performed well within the swarm. In our case, this happens starting at approximately $t=13$ for all $Ar$ numbers, for the region $z=[15\text {--}20]$. After that time, results are basically frozen for some characteristic times, that is up to the early decay regime, that is when all bubbles have left the region of observation. Furthermore, we have verified that results are statistically the same if the window $z=[20\text {--}25]$ is used, as for spectra. Of course, smoother profiles are obtained in the steady case because of the time averaging.

These results have been used to build the 3-D simulation described in the following section.

6. Three-dimensional bubble column

The 3-D bubble column is a direct extension of the previous 2-D numerical experiments: the cubic tank, of size $50 d_b \times 50 d_b \times 50 d_b$, is filled with a liquid and $256$ initially spherical bubbles are placed at the bottom within a region whose height is approximately $5 d_b$. The bubbles are homogeneously distributed in the lateral directions $x,y$, while avoiding any initial bubble overlap, and with a minimum separating distance of one diameter. This results in a local volume fraction in the region $0 \le z \le 7$ of $\phi \simeq 1\,\%$. The domain is closed at the bottom by a wall (no-slip boundary condition), and an outflow boundary condition is used at the top, while on the lateral sides the domain is periodic. At $t=0$ both the liquid and the bubbles are at rest. The dimensionless characteristic numbers of our numerical experiment are the following: $Ar =185;\ Eo=0.28;\ {\rho _l/\rho _b=800};\ \mu _l/\mu _b=100$. The configuration is in many respects very close to that investigated experimentally by Riboux et al. (Reference Riboux, Risso and Legendre2010). Following the 2-D analysis, the $Ar$ has been chosen large enough to trigger important nonlinear interactions. The Reynolds number is not well defined because of the unsteadiness, still a typical order of magnitude is $Re_b\sim 500$. This is in line with the results obtained with a single-bubble and imposing the same parameters (Appendix B).

From the numerical point of view, an adaptive mesh has been used with a maximum possible refinement of $N=4096$ cells in each direction, meaning a maximum resolution in terms of the bubble diameter of $d_b / \varDelta = 82$. The grid is refined or coarsened relying on the errors on the volume fraction and on the velocity components, using as absolute thresholds for the refinement the values $e_f = 0.01$ and $e_v = 0.003$, based on the analysis detailed in Appendix B. With such criteria of refinement it is possible to have the desired grid resolution in the regions where bubbles are present and where wakes develop, while in the remaining part of the domain where there is no agitation the grid is left coarser. The total number of computational cells grows in time because of the elongation of the wakes, starting from $N_{tot} \simeq 10^7$, and attaining $N_{tot} \simeq 9 \times 10^8$ at $t= 12$. Note that using a non-adaptive mesh would require $4096^3\approx 69\times 10^9$ grid points, which is beyond present computational capabilities. With respect to experiments we analyse the dynamics of a thin layer of bubbles rather than of a full swarm. As anticipated in the introduction, it turns out to be computationally too heavy to follow more bubbles than that. The present numerical experiment is, therefore, basically the best that can be done in simulating bubble column configurations today.

At a qualitative level, figure 5(a) shows the instantaneous motion of the bubbles at an early stage of the rising. The vorticity generated by the bubbles is included in elongated wakes. At this time, the transient has approximately finished and the thin layer of bubbles has stabilised to a width of approximately $7d_b$. A video is also given as supplementary material available at https://doi.org/10.1017/jfm.2021.288 to help the reader to have a clearer idea of the set-up.

Figure 5. (a) Snapshot of the 3-D simulation at $t=6$ after the release (a). Time is made non-dimensional with the bubble buoyancy time $\sqrt {d_b/g}$. The VOF field is shown with blue isosurfaces and the $\lambda _2$ vorticity field is shown with grey contours. The right-hand wall displays the level of mesh adaptation, while the left-hand wall displays the vertical component of the velocity field. (b) Bubble positions at different times, $t=6$ and $t=9$. Positions on the $y$ direction are collapsed on the plane. The lines correspond to the interrogation window.

In the same figure 5(b) we can see a lateral 2-D projection of the domain showing bubble positions at $t=6$ and at $t=9$, that is the last time used for the statistical analysis. The same procedure as in 2-D has been followed to acquire data and compute observables, with an interrogation window composed by the horizontal planes between $z=22$ and $z=25$, see figure 5(b).

We have acquired the data from each cell in this domain, i.e. considering both phases, and used them to compute the statistics in the time range $t\in [8.6,9.2]$. That approximately corresponds to the range over which bubbles are present in the whole window. More specifically, at $t=9.4$ only few bubbles are still present, and the statistics computed at $t=9.6$ turn out to be already strongly damped, as in earlier studies fluctuations are rapidly (exponentially) attenuated behind the swarm (Risso Reference Risso2018). In this time range, statistics are found to be roughly homogeneous in the interrogation window. We have, therefore, averaged over this space window to get better statistics, as done in the 2-D case. Statistics have been found to be also approximately steady in the time range considered, but with significant fluctuations and we have preferred to avoid time averaging. Such a choice for the statistical analysis region excludes the initial transient regime that concerns only the first few times. We display in figure 6 the p.d.f. of the vertical $z$ and horizontal $x$ velocity fluctuations (the $y$ component does not present appreciable statistical differences with respect to the $x$ component). As in the 2-D case, the velocity fluctuation field is computed at each $z$ subtracting the average velocity computed over the corresponding plane ${\boldsymbol u^\prime }={\boldsymbol U}-\langle \boldsymbol U \rangle _z$. We find the same characteristics reported in experiments (Riboux et al. Reference Riboux, Risso and Legendre2010), against which results are compared. The vertical velocity is strongly skewed, indicating a more important probability of having positive fluctuations, while the horizontal components are symmetric. Furthermore, both components are non-Gaussian, which is related to the complex features of the bubble-induced agitation. The agreement between numerical simulations and experiments is globally good. Yet in the numerical experiment the extreme events tend to be more frequent than in experiments, and exponential tails are found for $\sigma >rsim 3$. This may be related to the fact that the flow is here unsteady, and experiments may be under-sampling extreme events because they plausibly filter the smallest scales, and to the different measurement protocol. Since the p.d.f. of both components have been measured only in the experiments by Riboux et al. (Reference Riboux, Risso and Legendre2010), it is difficult to conclude. Finally, with respect to the 2-D case, figure 4, the p.d.f. is more skewed in the 3-D case. That is consistent with what is observed in experiments in a swarm confined in a thin gap (Bouche et al. Reference Bouche, Roig, Risso and Billet2014), although the wall friction effect is important there and thus the comparison is not conclusive. In figure 7(a), we present the spectrum of the kinetic energy computed in the same window. To compute the spectra, we have interpolated the data on a regular grid. To avoid spurious errors, we have eliminated the highest wave modes, so that spectra are calculated for 512 modes, although the maximum refinement is up to 4096 points. As for 2-D simulations, we have computed the spectrum at different times (not shown here), and we have found very little difference if spectra are computed at those times when bubbles are present in the plane used for the computation. When the bubbles have left the spatial region under investigation for a few characteristic times, agitation then decays rapidly, and an exponential fall-off is recorded. Figure 7(a) shows that bubbles are able to generate significant fluctuations in the length range $\lambda \in [10d_b,0.1 d_b]$, before being dissipated. After the energy range $\lambda \in [10d_b,1 d_b]$, the spectrum follows a power law with a scaling $E(k) \sim k^{-3}$ in an inertial range over the decade $\lambda \in [d_b,0.1 d_b]$. No hint of an inertial range with slope $k^{-5/3}$ is observed, neither at large or small scales.

Figure 6. Probability density functions of the velocity fluctuations in the lateral $x$ direction (a), and in the vertical $z$ (b), at time $t=9$, averaged over the space window $z=22\text {--}25$. Results are compared with the experimental data by Riboux et al. (Reference Riboux, Risso and Legendre2010), for a concentration of $\phi =1.7\,\%$, close to that of the numerical configuration. The points have been extracted directly from Risso (Reference Risso2018), and show the results obtained from measurement of the liquid agitation within the swarm. The error bars indicate the fluctuations recorded in the time range $t\in [8.6\text {--}9.2]$.

Figure 7. (a) Normalised spectrum of the kinetic energy. Present results ($Ar = 185$ and $Eo= 0.28$) are displayed at $t=8.8$ averaged over the space window $z=22\text {--}25$. The dashed line represents the $-3$ slope. Data obtained by Pandey et al. (Reference Pandey, Ramadugu and Perlekar2020) are shown for comparison at two different $Ar$ numbers. (b) Vorticity field for bubbles at $t=9.0$ and $z=25d_b$.

For comparison, the very recent numerical results presented by Pandey et al. (Reference Pandey, Ramadugu and Perlekar2020) are also shown for two $Ar$ numbers. It is worth recalling that these results have been obtained with a lower resolution ($24$ points per diameter instead of $82$ for the present simulation), and using a much lower density ratio of approximately $20$, instead of $800$ for the present simulation as for water/air. Despite these important differences, the results obtained in the present DNS are in quite good agreement with those obtained in Pandey et al. (Reference Pandey, Ramadugu and Perlekar2020) at $Ar=358$. A departure from DNS only appears at small scales around $1/10 d_b$, plausibly because of the lower resolution (256 points used in LES against 4096 used in the DNS here). The other simulation at $Ar=113$ seems instead to decay faster at all scales. We compare the results also against experiments. It can be observed that numerical and actual experiments do not investigate the same scales. Experiments are able to access a much larger domain, and they suggest that the $k^{-3}$ scaling might be valid over a larger range than that displayed in our results. On the other hand, numerical simulations appear to be more adequate to analyse accurately the small scales, where a possible change of slope from the Kolmogorov one is not recorded. Moreover, it is important to note that in experiments spectra are computed just behind the swarm, and not within it as in simulations. This may have an effect at small scales.

In order to qualitatively complement this analysis, we show in figure 7(b) the vorticity field on the same plane used to compute the spectra. This field highlights the position of the bubbles and the generation of vorticity at the scale of the diameter and slightly more. Several bubbles are still present at this time. In some cases it is apparent that different vortices have interacted, producing more complex structures.

While energy spectra contain key information about the flow, they cannot be used to disentangle the different mechanisms leading to the observed scalings, and a scale-by-scale analysis can be particularly useful (Alexakis & Biferale Reference Alexakis and Biferale2018). For that purpose, we apply a coarse-graining approach (Duchon & Robert Reference Duchon and Robert2000; Eyink & Sreenivasan Reference Eyink and Sreenivasan2006b) linked to the filtering approach used in LES (Germano Reference Germano1992), and recently applied to different turbulent configurations (Chen et al. Reference Chen, Ecke, Eyink, Rivera, Wan and Xiao2006b; Xiao et al. Reference Xiao, Wan, Chen and Eyink2009; Faranda et al. Reference Faranda, Lembo, Iyer, Kuzzay, Chibbaro, Daviaud and Dubrulle2018; Dubrulle Reference Dubrulle2019; Valori et al. Reference Valori, Innocenti, Dubrulle and Chibbaro2020). More specifically, we have applied this methodology to the velocity field, obtaining information about the energy flux and the dissipation. The advantage with respect to a spectral approach is that one can gain details also on the locality of the cascade, differentiating regions with positive or negative fluxes. Moreover, this spatial filtering approach is positive definite and local in space and can, therefore, be applied also in non-homogeneous flows.

In this filtering approach, the dynamic velocity field ${\boldsymbol u}$ is spatially (low-pass) filtered over a scale $\ell$ to obtain a filtered value $\bar {{\boldsymbol u}}_\ell ({\boldsymbol x})$ as follows:

(6.1)\begin{equation} \bar{{\boldsymbol u}}_\ell({\boldsymbol x}) = \int \textrm{d}^3 \,r G_\ell({\boldsymbol r}) {\boldsymbol u}({\boldsymbol x}+{\boldsymbol r}),\end{equation}

where $G_\ell$ is a smooth filtering function, spatially localised and such that $G_\ell (\boldsymbol r) = \ell ^{-3}G(\boldsymbol r/\ell )$ where the function $G$ satisfies $\int \textrm {d}\boldsymbol r \, G(\boldsymbol r)=1$, and $\int \textrm {d}\boldsymbol r \, \vert \boldsymbol r \vert ^2 G(\boldsymbol r) = O(1)$. By applying the filtering to the Navier–Stokes equations for the liquid phase we obtain the coarse-grained dynamics

(6.2)\begin{equation} \partial_{t} \bar{{\boldsymbol u}}_\ell+ (\bar{{\boldsymbol u}}_\ell \cdot {\boldsymbol{\nabla}})\bar{{\boldsymbol u}}_\ell ={-}\dfrac{1}{\rho}{\boldsymbol{\nabla}}\bar{p}_{\ell} -{\boldsymbol{\nabla}}\cdot{\boldsymbol{\tau}}_\ell+\nu\nabla^{2}\bar{{\boldsymbol u}}_\ell . \end{equation}

Since we focus on the liquid agitation, we neglect the gravity contribution, which acts as power injection through bubbles. In the same vein, at interfaces surface tension and density effects play a role. However, few bubbles are present in the analysed region, so that the impact should be small and we may retain the single-phase formulation given by (6.2). Furthermore, we are mainly interested at understanding whether a cascade process is active. To address this issue, the key term is given by ${\boldsymbol {\tau }}_\ell$, subscale stress tensor (or momentum flux), which describes the force exerted on scales larger than $\ell$ by fluctuations at scales smaller than $\ell$. It is given by

(6.3)\begin{equation} ({\boldsymbol{\tau}}_\ell)_{i,j} =\overline{(u_i u_j)}_\ell - (\bar{u}_\ell)_i (\bar{u}_\ell)_j . \end{equation}

The corresponding pointwise kinetic energy budget reads

(6.4)\begin{equation} \partial_t \left(\tfrac{1}{2}|\bar{\boldsymbol u}|^2\right) + \partial_j \widetilde{{ q}}_j={-}\varPi_\ell - \nu|{\boldsymbol{\nabla}}\bar{\boldsymbol u}|^2,\end{equation}

where we have dropped the $\ell$ subscript whenever unambiguous for the sake of clarity, and

(6.5a,b)\begin{equation} \widetilde{{ q}}_j=\left[ \left(\frac{1}{2}|\bar{\boldsymbol u}|^2 + \dfrac{1}{\rho}\bar{p} \right)) \bar{u}_j + \tau_{ij}\bar{u}_i - \nu\partial_j\left(\frac{1}{2}|\bar{\boldsymbol u}|^{2}\right)\right];\quad \varPi_\ell({\boldsymbol x}) \equiv{-}(\partial_{j}\bar{u}_{i})\tau_{ij}, \end{equation}

where $\widetilde {{\boldsymbol q}}_\ell$ is the transport term and $\varPi _\ell$ is the sub-grid scale (SGS) energy flux. This term is key since it represents the space-local transfer of energy among large and small scales across the scale $\ell$. The term $\varPi _\ell$ identifies the presence of a local direct (positive) or inverse (negative) energy cascade according to its sign. The last term in (6.4) represents the coarse-grained dissipation $\epsilon _\ell =\nu |{\boldsymbol {\nabla }}\bar {\boldsymbol u}|^2$. If a spatial average is done for different values of the filter width, one can find the average transfer of energy at each scale. Since our configuration is non-homogeneous in the vertical direction, the transport term $q_j$ in (6.4) is not zero. Yet, this term is related to spatial redistribution of energy and not directly linked to the cascade process, contrarily to $\varPi _\ell$ and $\epsilon _\ell$. For this reason, we have not analysed those terms.

In this work, we have applied a Gaussian filter defined as

(6.6)\begin{equation} G ({\boldsymbol r}) = \sqrt{\frac{6}{\rm \pi}} \exp({-}6 {\boldsymbol r}^2), \end{equation}

as typically used in LES (Pope Reference Pope2000). Since the flow is homogeneous in the horizontal direction, the filtering can be efficiently performed in spectral Fourier space, multiplying the quantity to be filtered by the Fourier transform of the filter

(6.7)\begin{equation} \hat{G}_\ell ({\boldsymbol k}) = \exp({-}k^2 \ell^2/24), \end{equation}

and then transforming back into physical space. In figure 8(a), we show the fluxes computed from the coarse-grained quantities defined in (6.5a,b). It is worth emphasising that fluxes are averaged in space but not in time. The physical features which unfold are the following.

1. (i) The scale-by-scale fluxes show variability in time, pointing out the statistical unsteadiness of the transfer processes.

2. (ii) Both inverse (negative flux) and direct (positive flux) cascades are found at different times. Both cascades involve more than a decade of scales. The direct cascade is more significant at scales smaller than the diameter, and the inverse cascade at scales larger than it.

3. (iii) Energy is transferred from scales around $\ell \approx d_b/2$ in both directions.

4. (iv) At around the same length scale the dissipation becomes significant.

5. (v) At smaller scales dissipation and direct flux become comparable.

Figure 8. (a) Mean energy flux at different filter lengths, the energy flux (solid lines) and the dissipative flux (dashed lines) are displayed at different times. Fluxes are computed at $z=25d_b$, the same as for the computations of spectra. The length scale displayed on the $x$ axis is normalised with the diameter $d_b$, so that $\ell =1$ corresponds to the initial bubble diameter. (b) Local energy flux at $t=9$ and $z=25d_b$. In the upper half-panel the filter length is $l=0.5 d_b$, in the lower half-panel $l=0.1 d_b$. The colour scale is the same in both half-panels. Each slice is made by $450$ points taken in the horizontal direction, and $180$ in the vertical direction.

Our scale-by-scale analysis points to the following physical picture of the pseudo-turbulent agitation induced by bubbles. Energy is injected by buoyancy (the only force at play here) and transferred by the bubbles into the liquid via the interface at scales comparable to the bubble diameter. The energy input $W_b$ must be proportional to the work made by buoyancy: $W_b \sim \phi g U_b$. Dissipation becomes significant at $\ell \sim d_b$, showing that fluctuations are mostly dissipated inside the small structures generated by bubble wake-interactions. At smaller scales than the diameter, there is a range where $W_{diss} \approx W_b$, which means $\nu (\delta u_\ell )^2\ell ^2 \sim \phi g U_b$, where we have considered the two-point quantities $\delta u_\ell =u(x+\ell )-u(x)$. This gives the scaling behaviour $\delta u_\ell ^2 \sim \ell ^2$, which means in spectral space $E(k) \sim k^{-3}$. We obtain here the scaling with an argument similar to that used by Lance & Bataille (Reference Lance and Bataille1991) and Prakash et al. (Reference Prakash, Mercado, van Wijngaarden, Mancilla, Tagawa, Lohse and Sun2016), yet in the physical space rather than in the spectral one. The fact that $\varPi _\ell$ changes sign shows that both a direct and an inverse transfer of energy through nonlinear terms are active, and dominate at different times. On average the transfer is more towards small scales, but the inverse process is not negligible. The copresence of direct and inverse cascades intermittently is a feature also of fluid turbulence (Chorin Reference Chorin2013). While the direct transfer is linked to dissipation, the inverse one is related to the formation of the wakes, which are found to develop up to some characteristic lengths.

To further understand the mechanisms indicated, we show in figure 8(b) a slice of the energy flux at two different scales: half the diameter and a smaller scale. The pictures show that the energy flux and dissipation are concentrated in the wakes generated by the bubbles. Furthermore, these structures, initially at the scale of the diameter, may become a little larger, indicating the generation of larger eddies, and are eventually dissipated at small scales, where the imprint of the bubbles is still detectable.