1 Introduction
Bubble columns are commonly used in a broad range of technologies, notably in the chemical and biochemical industries. Simple bubble columns do not require active stirring and can therefore operate without interior moving parts. The large surface area between gas and liquid is useful for mass and species transfer, possibly involving chemical reactions (e.g. Deckwer Reference Deckwer1992) such as in air-lift bioreactors, an example of which is in the treatment of wastewater. Bubble columns are also used for this reason in direct contact heat transfer (e.g. Hewitt, Shires & Bott Reference Hewitt, Shires and Bott1994). Besides offering a large surface area, rising bubbles agitate the liquid flow, which results in enhanced mixing that is usually desired, but this also poses a modelling difficulty. Similar mixing arises also in the diffusion through porous media in the presence of flow, which has been well studied previously, but mostly for fixed beds of particulates, often under creeping flow (e.g. Batchelor Reference Batchelor1974; Koch & Brady Reference Koch and Brady1985). Mixing in bubble columns is complicated further by the fact that liquid velocity fluctuations are coupled with the dynamics of (deformable) bubbles, usually beyond creeping flow (e.g. Alméras et al. Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015).
In the present study, we consider transport of a scalar (such as the concentration of a chemical species, or the temperature) through incompressible bubbly flows. Gradients of temperature and concentration may, in general, induce fluid motion and influence the velocity field through changes in density and viscosity, or through the interface rheology. If these effects are small, as assumed herein, temperature and solute concentration can be considered as passive scalars. Although the arbitrary choice was made in this study to use the terminology of the mass transfer problem, the results carry over to thermal applications (upon assuming that effects of viscous heating can be ignored).
Our present main interest is the formulation and closure of conservation equations and constitutive relations governing the dispersion of such a scalar in a bubbly suspension over scales (termed hereinafter the ‘macroscale’) that are much larger than the bubble size (termed hereinafter the ‘microscale’). Under the assumption of macroscale homogeneity and stationarity, scalar dispersion in multiphase systems can be described by a macroscale version of Fick’s (or Fourier’s) law, which relates the macroscale scalar flux to the macroscale scalar gradient through an effective diffusivity tensor (or effective conductivity tensor in thermal applications) (Batchelor Reference Batchelor1974; Koch & Brady Reference Koch and Brady1985, Reference Koch and Brady1987). This effective diffusivity is defined from an Eulerian perspective. Experimentally, scalar dispersion is usually investigated from a Lagrangian point of view. In the Lagrangian framework, the effective diffusivity is defined as the long-time limit of the time rate of change of a fluid tracer’s mean-square displacement, that is, as a measure of spread about the mean position. Koch & Brady (Reference Koch and Brady1987) demonstrated that the Lagrangian effective diffusivity is equivalent to the symmetric part of the Eulerian effective diffusivity, and that the antisymmetric part of the Eulerian effective diffusivity is associated with anisotropic microstructures.
Scalar dispersion in a suspension of particulates (bubbles, drops or rigid particles) results from two processes of very different nature: the diffusion by Brownian motion of the molecules, and the convection by the fluid velocity disturbances induced by the particulate motion. The relative importance of these two processes is measured by the Péclet number
$Pe=Ud_{b}/D$
, where
$U$
is the characteristic velocity of the particulates relative to that of the system (defined in § 2),
$d_{b}$
is the characteristic size of the particulates and
$D$
is the diffusivity of the bulk. In the limit
$Pe=0$
, the effective diffusivity is purely diffusive and depends only on the particulate-to-bulk diffusivity ratio, possible discontinuity of the scalar at the interface, particulate volume fraction and suspension microstructure (i.e. the positions, shapes and orientations of the inclusions). This particular situation is essentially relevant to heat and electricity conduction in composite materials. When
$Pe\gg 1$
, the dominant contribution to the effective diffusivity is due to convective mixing. This last regime is that generally encountered in bubbly flows.
Recently Alméras et al. (Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015) investigated experimentally the dispersion of a low-diffusive dye within a homogeneous swarm of high-Reynolds-number rising bubbles at
$Pe=O(10^{6})$
; herein we define the Reynolds number as
$Re=Ud_{b}/\unicode[STIX]{x1D708}_{c}$
, where
$\unicode[STIX]{x1D708}_{c}$
is the kinematic viscosity of the liquid. They showed that scalar mixing primarily results from pseudo-turbulence, i.e. from the liquid agitation produced by bubble wake interactions, and can be modelled in a manner analogous to dispersion in shear-induced turbulence (Taylor Reference Taylor1921). Apart from the work of Alméras et al. (Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015), the only other experimental investigation of mixing in homogeneous bubbly flows reported in the literature is the preliminary study of Mareuge & Lance (Reference Mareuge and Lance1995), which consists of a single data point. To the best of our knowledge, neither theoretical nor numerical investigations of scalar mixing in homogeneous bubbly flows have been reported thus far. Theoretical work is, however, available for other types of multiphase systems, and we shall review these now.
The determination of such an effective diffusivity, at the macroscale, necessitates consideration of the conditions at the microscale. One class of analytical work is devoted to the study of dilute systems with fixed random microstructure, for instance, as a model of a porous medium. In the absence of convection (
$Pe=0$
), the analytical expression of the effective diffusivity is available in the dilute limit from analysis of the corresponding problem in conduction of heat or electricity through a dispersed medium (e.g. Maxwell Reference Maxwell1873; Jeffrey Reference Jeffrey1973). The problem of scalar dispersion in the presence of a bulk convective motion (
$Pe>0$
) has been analysed by Koch & Brady (Reference Koch and Brady1985) for Stokes flow through a random bed of fixed solid spheres. Using the method of conditional averaging pursued earlier by Hinch (Reference Hinch1977), they carried out an asymptotic analysis in low volume fraction of the effective diffusivity for all values of the Péclet number. Three mechanisms causing dispersion at high Péclet number were identified: mechanical dispersion resulting from the stochastic velocity field in the bulk, which is independent of Brownian diffusion and grows as
$Ud_{b}$
; holdup dispersion in stagnant and recirculating regions, which is proportional to
$U^{2}d_{b}^{2}/D$
; and boundary-layer dispersion, which grows as
$Ud_{b}\ln (Ud_{b}/D)$
near the solid particle surfaces.
Another class of analytical studies assumes a periodic microstructure. For the pure diffusion problem (
$Pe=0$
), analytical solutions have been derived for a composite material consisting of regularly arranged spheres embedded in a homogeneous matrix (Rayleigh Reference Rayleigh1892; Sangani & Acrivos Reference Sangani and Acrivos1983), and the effect of anisotropy has been investigated by considering periodic arrangements of spheroidal inclusions (Kushch Reference Kushch1997; Harfield Reference Harfield1999). In the presence of convection (
$Pe>0$
), the general theory of dispersion developed by Brenner (Reference Brenner1980) and Brenner & Adler (Reference Brenner and Adler1982) provides a consistent framework for determining the effective diffusivity in spatially periodic media. Koch et al. (Reference Koch, Cox, Brenner and Brady1989) carried out explicit calculations for a periodic porous medium consisting of fixed solid particles arranged in a cubic lattice and embedded in a continuous phase under Stokes flow conditions. They showed that, in ordered systems, the mechanical dispersion encountered in random media is absent, and that at high Péclet number, either Taylor dispersion, growing as
$U^{2}d_{b}^{2}/D$
, or enhanced diffusion, which is proportional to
$D$
, is obtained depending on the direction of the mean flow relative to the lattice structure.
In bubbly flows, the spatial arrangement of the inclusions evolves in time, the microstructure of the suspension is unknown a priori and Stokes flow is usually not applicable. For these reasons, previous analyses are, a priori, not applicable to bubbly suspensions. Nevertheless, we showed in earlier work (Loisy, Naso & Spelt Reference Loisy, Naso and Spelt2017) that the dynamics of freely evolving bubbly suspensions at moderate Reynolds number shares some common features with that of ordered arrays of bubbles. It is therefore of fundamental interest to investigate, contrast and compare the mixing properties of ordered and freely evolving bubbly suspensions in light of prior asymptotic analyses for ordered and random arrangements of rigid particles.
In this paper we investigate scalar dispersion, by determining the effective diffusivity in ordered and freely evolving bubbly suspensions, specifically, the contribution of bubble-induced velocity disturbances thereof. The prior work outlined above has established that, in the systems studied therein, the effective diffusivity can be much larger than that in each of the fluids involved, even if the diffusivity in the two media is the same and the scalar is continuous at the surface of particulates. In view of the already significant number of parameters involved, we shall therefore adopt this restriction here. Such a simplified approach will not provide an accurate description of real bubbly flows, but should shed some light on the fundamental mechanisms of mixing in these systems.
The paper is organized as follows. The theoretical framework and problem statement are provided in § 2. Our numerical approach to compute the effective diffusivity is presented, and followed by a description of the regimes and the range of parameter values that are investigated herein, in § 3. The first objective (in § 4) is to elucidate the role played by liquid inertia in ordered suspensions, using direct numerical simulation (DNS) and analysis. The second objective (in § 5) is to investigate the effective diffusivity of freely evolving suspensions for a wide range of Péclet numbers, to compare it with that obtained for ordered suspensions, and to evaluate the effect of introducing additional degrees of freedom in the system. Finally, the main results and perspectives of this work are provided in § 6.
2 Problem statement
The local evolution of the passive scalar
$c$
in each fluid is governed by
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{q}=0,\end{eqnarray}$$
where
$\boldsymbol{q}$
is the flux of scalar given by
$$\begin{eqnarray}\boldsymbol{q}=\boldsymbol{u}c-D\unicode[STIX]{x1D735}c,\end{eqnarray}$$
with
$\boldsymbol{u}$
the fluid velocity and
$D$
the constant scalar diffusivity. We assume that the scalar and its gradient are continuous across the interface, and phase change is not considered in this study. Under these assumptions, no distinction between the phases is needed for the scalar transport, which is described by (2.1) in the entire system. We return to these restrictions in § 2.1 and in the Conclusions section (§ 6); the objective here is to study this key basic reference problem.
In the context of heat transfer, equation (2.1) derives from the energy balance upon neglecting viscous heating; in this case
$c$
would represent the temperature, continuous at the interface, and
$D$
the thermal diffusivity as defined by Fourier’s law, assumed to be equal in both gas and liquid. In the context of mass transfer, equation (2.1) describes the transport of a chemical species present at very low concentration
$c$
so that Fick’s law describes the conservation of mass, neglecting any difference in molecular diffusivity
$D$
and solubility of the species in the two phases. While the assumption of equal molecular diffusivities is never satisfied in real systems, the assumption of a unit dimensionless Henry’s constant is reasonably applicable to, for example, carbon dioxide dispersion in the air–water system.
The fluid motion in the gas and liquid is governed by the incompressible Navier–Stokes equations, which are coupled at the interface by the appropriate jump conditions, namely the continuity of velocity and of tangential traction across the interface, and a jump in normal traction due to surface tension.
2.1 Macroscale description
The problem we are concerned with here is the modelling of scalar transport at a macroscale, that is, at the scale whereat the suspension may be seen as a homogeneous continuum, without distinction between the two phases. In order to obtain such a macroscopic description, we consider an ensemble of realizations of the suspension, these realizations having the same macroscopic conditions (e.g. fluid properties, gas volume fraction) but different microscopic configurations (e.g. bubble individual positions, shapes and velocities), and average over those realizations. In concrete terms, ensemble averaging would be realized by averaging over a large number of experiments run under identical macroscopic conditions. The ensemble-averaged transport equation is obtained from ensemble averaging the local transport equation (2.1). It reads
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\langle c\rangle }{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\langle \boldsymbol{q}\rangle =0,\end{eqnarray}$$
where
$\langle ~\rangle$
denotes the ensemble average operator, and where the ensemble-averaged flux is given by
$$\begin{eqnarray}\langle \boldsymbol{q}\rangle =\langle \boldsymbol{u}\rangle \langle c\rangle -D\unicode[STIX]{x1D735}\langle c\rangle +\langle \boldsymbol{u}^{\prime }c^{\prime }\rangle ,\end{eqnarray}$$
where the velocity fluctuations are defined by
$\boldsymbol{u}^{\prime }=\boldsymbol{u}-\langle \boldsymbol{u}\rangle$
and the scalar fluctuations by
$c^{\prime }=c-\langle c\rangle$
. Under the restrictions set out above, the average flux consists of three contributions: (i)
$\langle \boldsymbol{u}\rangle \langle c\rangle$
is the advection of the average scalar field at the average system velocity; (ii)
$-D\unicode[STIX]{x1D735}\langle c\rangle$
is the diffusion of the average scalar field directly by the average scalar gradient; and (iii)
$\langle \boldsymbol{u}^{\prime }c^{\prime }\rangle$
corresponds to the advection of the scalar fluctuations by the velocity fluctuations induced by bubble motion.
When the suspension is statistically homogeneous and in a statistically stationary state, the linearity in
$c$
of the local flux (2.1b
) results, in the presence of an imposed constant average scalar gradient, in a macroscale constitutive relation of the form (Koch & Brady Reference Koch and Brady1985, Reference Koch and Brady1987)
$$\begin{eqnarray}\langle \boldsymbol{q}\rangle =\langle \boldsymbol{u}\rangle \langle c\rangle -\unicode[STIX]{x1D63F}^{\,eff}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\langle c\rangle ,\end{eqnarray}$$
where
$\unicode[STIX]{x1D63F}^{\,eff}$
is a constant effective diffusivity tensor. Comparison of the effective diffusivity definition (2.4) with the average flux expression (2.3) yields the expression for the effective diffusivity. In order to reflect the contributions to the scalar flux identified above, it is customary to write the effective diffusivity as
$$\begin{eqnarray}\unicode[STIX]{x1D63F}^{\,eff}=D\unicode[STIX]{x1D644}+\unicode[STIX]{x1D63F}^{\,conv},\end{eqnarray}$$
where
$$\begin{eqnarray}\unicode[STIX]{x1D63F}^{\,conv}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\langle c\rangle =-\langle \boldsymbol{u}^{\prime }c^{\prime }\rangle\end{eqnarray}$$
is the convective contribution arising from bubble-induced velocity fluctuations. For this model to be complete, one must find a closure relation for
$\unicode[STIX]{x1D63F}^{\,conv}$
only in terms of macroscopic quantities appearing in the problem statement. We recall here that further contributions to the average flux (2.3) and hence to the effective diffusivity (2.5) arise if the diffusivities in the fluids are not the same, or if the concentration is discontinuous at fluid–fluid interfaces (e.g. Batchelor & O’Brien Reference Batchelor and O’Brien1977; Koch & Brady Reference Koch and Brady1985). We return to the significance of this in § 6 below.
2.2 Effective transport properties
To determine the effective diffusivity for (unbounded) homogeneous bubbly suspensions, we represent such flows by the periodic repetition of a cubic unit cell containing a finite number
$N_{b}$
of freely moving bubbles of equal volume, building on our prior work on the dynamics of bubbles for this model system (Loisy et al.
Reference Loisy, Naso and Spelt2017). In the limit
$N_{b}=1$
, one obtains a simple cubic array of bubbles, which is of interest as a model of perfectly ordered suspensions. The opposite limit of large
$N_{b}$
is of interest as a model of real suspensions, although convergence with the number of bubbles would have to be verified. We shall refer hereinafter to this set-up with one bubble in the cell as an ordered array, and to that with more than one bubble in the unit cell as a free array.
The bubbles rise under the sole effect of buoyancy. Herein, an upward-pointing primary axis
$\boldsymbol{e}_{3}$
of the periodic arrangement is taken to be aligned with gravity (with the exception of the more general analysis presented in § 4.1). From symmetry arguments, and adopting a Cartesian coordinate system,
$$\begin{eqnarray}\unicode[STIX]{x1D63F}^{\,conv}=\left[\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D60B}_{\bot }^{\,conv} & \unicode[STIX]{x1D60B}_{12}^{\,conv} & \unicode[STIX]{x1D60B}_{13}^{\,conv}\\ \unicode[STIX]{x1D60B}_{12}^{\,conv} & \unicode[STIX]{x1D60B}_{\bot }^{\,conv} & \unicode[STIX]{x1D60B}_{13}^{\,conv}\\ \unicode[STIX]{x1D60B}_{31}^{\,conv} & \unicode[STIX]{x1D60B}_{31}^{\,conv} & \unicode[STIX]{x1D60B}_{\Vert }^{\,conv}\\ \end{array}\right],\end{eqnarray}$$
where we have introduced the longitudinal and transverse components of the convective contribution to the effective diffusivity, which are denoted by
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
and
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
, respectively, and defined by
$$\begin{eqnarray}\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}=\unicode[STIX]{x1D60B}_{33}^{\,conv}\quad \text{and}\quad \unicode[STIX]{x1D60B}_{\bot }^{\,conv}=\unicode[STIX]{x1D60B}_{11}^{\,conv}=\unicode[STIX]{x1D60B}_{22}^{\,conv}.\end{eqnarray}$$
Our first goal is to characterize the effects of liquid inertia (through
$Re$
) on the dependence of
$\unicode[STIX]{x1D63F}^{\,conv}$
on
$Pe$
for ordered suspensions (
$N_{b}=1$
), thereby extending prior work on dilute ordered arrays of rigid spheres in Stokes flow conditions (Koch et al.
Reference Koch, Cox, Brenner and Brady1989). Our second goal is to evaluate the effect of introducing additional degrees of freedom in the system (through increasing
$N_{b}$
), and to investigate the dependence of
$\unicode[STIX]{x1D63F}^{\,conv}$
on
$Pe$
in freely evolving suspensions (sufficiently large
$N_{b}$
). As we found the off-diagonal components to be zero in all configurations that we investigated, only results for the longitudinal and the transverse components of
$\unicode[STIX]{x1D63F}^{\,conv}$
will be presented.
In dimensionless groups, we shall use as characteristic length scale the bubble size
$d_{b}$
, which is defined, since bubbles are deformable, as the (equivalent) diameter of a sphere of the same volume. The characteristic velocity
$U$
is taken here as the bubble rise velocity in the frame of the suspension (the so-called drift velocity
$\langle \boldsymbol{U}\rangle =\langle \boldsymbol{u}\rangle _{d}-\langle \boldsymbol{u}\rangle$
, where the first term is the volume average of velocity on the disperse phase only and the second one is the same average in the entire system). As already mentioned, a key dimensionless group appearing in the scalar transport problem is the Péclet number
$Pe=Ud_{b}/D$
, which compares advective and diffusive transport. Our main objective is to elucidate the effect of the value of
$Pe$
on the effective diffusivity using analytical and numerical methods.
The effective diffusivity necessarily also depends on the gas volume fraction
$\unicode[STIX]{x1D719}=(N_{b}\unicode[STIX]{x03C0}d_{b}^{3})/(6h^{3})$
, where
$h$
is the linear size of the unit cell. The analytical and computational methods used here pose some restrictions on the range of
$\unicode[STIX]{x1D719}$
values that can be studied herein; we postpone discussion of that to the pertinent sections below. We also consider the effects of the number of bubbles in the periodic cell,
$N_{b}$
, which affects the order in the suspension:
$N_{b}=1$
corresponds to a cubic array, whereas more bubbles results in a different microstructure (the latter term encompasses all the information about the statistical distribution of the bubble positions, shapes, orientations, etc.). Since scalar transport is coupled to momentum transport, the bubble Reynolds number
$Re=Ud_{b}/\unicode[STIX]{x1D708}_{c}$
may also play a significant role that will be investigated here as well. The ranges of
$\unicode[STIX]{x1D719}$
,
$N_{b}$
and
$Re$
studied here are summarized in table 1.
As the bubbly flows we consider are buoyancy-driven, a difficulty arises from the fact that
$U$
is a priori unknown, and depends in a complex manner on
$N_{b}$
,
$\unicode[STIX]{x1D719}$
, the density and viscosity ratios between both phases, the Archimedes (or Galileo) number
$Ar=\sqrt{\unicode[STIX]{x1D70C}_{c}|\unicode[STIX]{x1D70C}_{d}-\unicode[STIX]{x1D70C}_{c}|gd_{b}^{3}}/\unicode[STIX]{x1D707}_{c}$
, and the Bond (or Eötvös) number
$Bo=|\unicode[STIX]{x1D70C}_{d}-\unicode[STIX]{x1D70C}_{c}|gd_{b}^{2}/\unicode[STIX]{x1D6FE}$
, where the subscripts
$d$
and
$c$
refer to the disperse (gas) and continuous (liquid) phases, respectively,
$g$
is the magnitude of the gravitational acceleration,
$\unicode[STIX]{x1D70C}$
denotes density,
$\unicode[STIX]{x1D707}$
is the dynamic viscosity and
$\unicode[STIX]{x1D6FE}$
is the surface tension. In most bubbly flows of practical relevance, the gas-to-liquid density and viscosity ratios are vanishingly small. Their precise values are not important from a physical point of view as long as they are small enough; in the simulations, the gas-to-liquid density and viscosity ratios were set to
$\unicode[STIX]{x1D70C}_{d}/\unicode[STIX]{x1D70C}_{c}=10^{-3}$
and
$\unicode[STIX]{x1D707}_{d}/\unicode[STIX]{x1D707}_{c}=10^{-2}$
, respectively. The dependence of
$U$
on (
$Ar,Bo,\unicode[STIX]{x1D719},N_{b}$
) has been addressed in Loisy et al. (Reference Loisy, Naso and Spelt2017) and is not further discussed here. In the present study, we shall therefore assume that
$U$
is known.
Simulated flow configurations:
$Bo$
and
$Ar$
define the flow regime,
$N_{b}$
is the number of free bubbles in the unit cell and
$\unicode[STIX]{x1D719}$
is the gas volume fraction. The resulting bubble Reynolds number (
$Re$
) and shape are also provided.
3 Methodology
For convenience of numerical implementation, we reorganize the problem formulation by introducing the decomposition
$$\begin{eqnarray}c=\bar{c}+\tilde{c},\end{eqnarray}$$
where
$\bar{c}$
is the imposed constant linear scalar field
$$\begin{eqnarray}\bar{c}=\unicode[STIX]{x1D735}\langle c\rangle \boldsymbol{\cdot }\boldsymbol{x}.\end{eqnarray}$$
The advantage of this decomposition is that the disturbance field
$\tilde{c}$
is then spatially periodic. The governing equation for this disturbance field is
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\tilde{c}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}\tilde{c})-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(D\unicode[STIX]{x1D735}\tilde{c})=-\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\langle c\rangle ,\end{eqnarray}$$
which is the equation we integrate numerically. The convective contribution to the effective diffusivity is then calculated from
$$\begin{eqnarray}\unicode[STIX]{x1D63F}^{\,conv}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\langle c\rangle =-\langle \boldsymbol{u}^{\prime }\tilde{c}\rangle ,\end{eqnarray}$$
which can be shown to be equivalent to (2.6). In this expression,
$\langle ~\rangle$
is defined as an ensemble average operator, as above. For statistically homogeneous and stationary systems, as considered here, it is inferred from ergodicity that ensemble averaging is identical to volume and time averaging. As a consequence,
$\unicode[STIX]{x1D63F}^{\,conv}$
is computed from (3.4) with the ensemble average being replaced in practice by a volume average combined with a time average over an appropriate time period.
3.1 Numerical method
Thus, the components of
$\unicode[STIX]{x1D63F}^{\,conv}$
are obtained from DNS by imposing a constant linear scalar field
$\bar{c}$
and determining the resulting periodic disturbance scalar field. Two distinct simulations are required to fully determine the five independent components of
$\unicode[STIX]{x1D63F}^{\,conv}$
: in one simulation,
$\unicode[STIX]{x1D735}\bar{c}=\boldsymbol{e}_{3}$
, which yields
$\unicode[STIX]{x1D60B}_{13}^{\,conv}$
and
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
; in the other simulation,
$\unicode[STIX]{x1D735}\bar{c}=\boldsymbol{e}_{1}$
, which yields
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
,
$\unicode[STIX]{x1D60B}_{12}^{\,conv}$
and
$\unicode[STIX]{x1D60B}_{31}^{\,conv}$
. The off-diagonal components of
$\unicode[STIX]{x1D63F}^{\,conv}$
were found to be zero (up to computer accuracy for ordered arrays, and statistical uncertainty for free arrays) for all the sets of parameters we considered, and therefore will not be shown.
The numerical methods employed to solve the two-phase flow have been described in detail in Loisy et al. (Reference Loisy, Naso and Spelt2017). In short, we employ a standard projection method (Chorin Reference Chorin1968) to integrate the incompressible Navier–Stokes equations, a level-set method (e.g. Sussman, Smereka & Osher Reference Sussman, Smereka and Osher1994) to capture the moving gas–liquid interface, and surface tension is accounted for using the continuum surface force model (Brackbill, Kothe & Zemach Reference Brackbill, Kothe and Zemach1992).
Our algorithm proceeds iteratively through the following steps:
-
(i) The position of the interface is first advanced in time according to the modified level-set method of Sabelnikov, Ovsyannikov & Gorokhovski (Reference Sabelnikov, Ovsyannikov and Gorokhovski2014) using a third-order total-variation-diminishing (TVD) Runge–Kutta scheme. The level-set function is then reinitialized using the procedure of Russo & Smereka (Reference Russo and Smereka2000), and a correction is finally applied to enforce volume conservation.
-
(ii) The scalar transport equation (3.3) is advanced by using a mixed Crank–Nicolson/third-order Adams–Bashforth time-stepping scheme.
-
(iii) The time integration of the incompressible Navier–Stokes equations is then carried out using a mixed Crank–Nicolson/third-order Adams–Bashforth scheme and consists of the combination of a predictor step, where a temporary velocity field is estimated by ignoring the effect of pressure, and a corrector step, where the velocity field is corrected by the pressure gradient term computed from the divergence-free condition.
Spatial discretization relies on a mixed finite difference/finite volume approach on a fixed, staggered, Cartesian grid. Second-order centred schemes are generally employed, except for advective terms, which are discretized using fifth-order weighted essentially non-oscillatory (WENO) schemes.
Results of numerical tests are presented in the Appendix.
3.2 Parametric study
Four different flow regimes, as defined by the set (
$Ar,Bo$
), are considered here. These are described in table 1, and have been studied in Loisy et al. (Reference Loisy, Naso and Spelt2017) (the same case code names are used). In case S0, the bubbles are spherical and the Reynolds number is vanishingly small, which approaches Stokes flow conditions. In case S1, the bubbles are (nearly) spherical and
$Re\gtrsim 1$
. In case C, the bubbles are skirted, and
$Re\approx 8$
. In case E1, the bubbles are ellipsoidal, and
$Re\approx 30{-}40$
.
Ordered arrays of bubbles in these four flow regimes have been considered for the smallest volume fraction numerically accessible (value provided in table 1). After a transient regime, all ordered suspensions considered here are in a strictly steady state (for the flow and the scalar) during which the results presented in § 4 were obtained. Simulations of scalar transport in free arrays have been performed for
$2\leqslant N_{b}\leqslant 12$
in case E1 at
$\unicode[STIX]{x1D719}=2.4\,\%$
. In these conditions, coalescence is indeed absent (it does occur at larger
$\unicode[STIX]{x1D719}$
), whereas simulations at lower
$\unicode[STIX]{x1D719}$
for free arrays are excessively expensive for the method and facilities used. In this regime, the system is in an unsteady but statistically stationary state (for the flow and the scalar), during which the statistics presented in § 5 have been measured. For each of these configurations (
$Ar,Bo,\unicode[STIX]{x1D719},N_{b}$
), the drift velocity (and thereby the Reynolds number) is known from Loisy et al. (Reference Loisy, Naso and Spelt2017). This allowed us to impose the Péclet number a priori.
The numerical simulation results for ordered arrays are compared with the results of analysis at small (but possibly finite) Reynolds number and small volume fraction.
4 Ordered suspensions
We examine in this section the dispersion of a passive scalar in ordered suspensions of deformable bubbles. Our main objective here is to elucidate the effects of inertia on dispersion, using theoretical analysis and numerical simulation.
4.1 Asymptotic analysis
We first determine analytically the convective contribution to the effective diffusivity of ordered suspensions of spherical fluid particulates (bubbles or drops). The Reynolds number of the particulates is assumed to be small so that the Navier–Stokes equations can be approximated by the Oseen equations.
4.1.1 General solution
An ordered array of particulates translating at a drift velocity
$\boldsymbol{U}$
is equivalent to an ordered array of fixed particulates immersed in a viscous fluid moving with an average system velocity
$\langle \boldsymbol{u}\rangle =-\boldsymbol{U}$
. The centres of the particulates are located on the nodes of a simple cubic lattice:
$$\begin{eqnarray}\boldsymbol{r}_{\boldsymbol{n}}=h\,(n_{1}\boldsymbol{e}_{1}+n_{2}\boldsymbol{e}_{2}+n_{3}\boldsymbol{e}_{3}),\quad n_{1},n_{2},n_{3}=0,\pm 1,\pm 2,\ldots ,\end{eqnarray}$$
where
$h$
is the lattice spacing and
$\boldsymbol{e}_{i}$
are the unit vectors aligned with the primitive axes of the cubic lattice. In the dilute limit (
$d_{b}/h\ll 1$
), the action of these particulates on the fluid can be represented by point forces
$-\boldsymbol{f}$
. The convective contribution to the effective diffusivity arising from the far field has been derived by Koch et al. (Reference Koch, Cox, Brenner and Brady1989) for an ordered array of rigid spheres in the Stokes flow regime. In what follows we extend their result to the case of spherical fluid particulates at small but finite
$Re$
.
When
$Pe\ll 1$
, the convective contribution to the effective diffusivity arising from the far field can be approximated by (Koch et al.
Reference Koch, Cox, Brenner and Brady1989)
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D63F}^{\,conv}}{D}=\mathop{\sum }_{\boldsymbol{k}\neq \mathbf{0}}\frac{k^{2}\hat{\boldsymbol{u}}^{\prime }(\boldsymbol{k})\hat{\boldsymbol{u}}^{\prime }(-\boldsymbol{k})}{(2\unicode[STIX]{x03C0})^{2}k^{4}D^{2}+(\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{k})^{2}},\end{eqnarray}$$
where the summation is over all vectors
$\boldsymbol{k}$
in the reciprocal lattice,
$$\begin{eqnarray}\boldsymbol{k}=\frac{1}{h}\,(n_{1}\boldsymbol{e}_{1}+n_{2}\boldsymbol{e}_{2}+n_{3}\boldsymbol{e}_{3}),\end{eqnarray}$$
and where
$\hat{\boldsymbol{u}}^{\prime }$
is the three-dimensional Fourier transform of the velocity disturbance
$\boldsymbol{u}^{\prime }=\boldsymbol{u}-\langle \boldsymbol{u}\rangle$
. In Oseen flow past an ordered array of point particulates,
$\hat{\boldsymbol{u}}^{\prime }$
is given by
$$\begin{eqnarray}\hat{\boldsymbol{u}}^{\prime }(\boldsymbol{k})=\frac{\boldsymbol{f}\boldsymbol{\cdot }(\boldsymbol{k}\boldsymbol{k}/k^{2}-\unicode[STIX]{x1D644})}{(2\unicode[STIX]{x03C0}k)^{2}h^{3}\unicode[STIX]{x1D707}_{c}+\text{i}2\unicode[STIX]{x03C0}h^{3}\unicode[STIX]{x1D70C}_{c}\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{k}},\quad \boldsymbol{k}\neq \mathbf{0},\end{eqnarray}$$
where
$\boldsymbol{f}$
is the hydrodynamic force exerted by the ambient fluid on a particulate. In the dilute limit,
$\boldsymbol{f}$
can be approximated by the Oseen drag exerted on a single spherical fluid particulate,
$$\begin{eqnarray}\boldsymbol{f}=F\boldsymbol{f}_{0,Stokes},\end{eqnarray}$$
where
$\boldsymbol{f}_{0,Stokes}$
is the Stokes drag on that particulate (Hadamard Reference Hadamard1911; Rybczynski Reference Rybczynski1911),
$$\begin{eqnarray}\boldsymbol{f}_{0,Stokes}=-2\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}^{\ast }\unicode[STIX]{x1D707}_{c}d_{b}\boldsymbol{U},\quad \text{with}~\unicode[STIX]{x1D707}^{\ast }=\frac{\unicode[STIX]{x1D707}_{c}+3\unicode[STIX]{x1D707}_{d}/2}{\unicode[STIX]{x1D707}_{c}+\unicode[STIX]{x1D707}_{d}},\end{eqnarray}$$
and where
$F$
accounts for the finite-
$Re$
correction to the Stokes drag (Brenner & Cox Reference Brenner and Cox1963),
$$\begin{eqnarray}F=1+{\textstyle \frac{1}{8}}\unicode[STIX]{x1D707}^{\ast }Re.\end{eqnarray}$$
The convective contribution to the effective diffusivity of a dilute ordered array of fluid particulates in Oseen flow conditions is therefore
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D63F}^{\,conv}}{D}=\frac{\unicode[STIX]{x1D707}^{\ast 2}}{(2\unicode[STIX]{x03C0})^{2}}\frac{d_{b}^{2}}{h^{2}}F^{2}\unicode[STIX]{x1D63E},\end{eqnarray}$$
where
$\unicode[STIX]{x1D63E}$
is the dimensionless tensor
$$\begin{eqnarray}\unicode[STIX]{x1D63E}=\mathop{\sum }_{\boldsymbol{k}^{\ast }\neq \mathbf{0}}\frac{\left[\boldsymbol{U}^{\ast }\boldsymbol{\cdot }\left({\displaystyle \frac{\boldsymbol{k}^{\ast }\boldsymbol{k}^{\ast }}{k^{\ast 2}}}-\unicode[STIX]{x1D644}\right)\right]^{2}}{k^{\ast 2}\left[{\displaystyle \frac{(2\unicode[STIX]{x03C0})^{2}k^{\ast 4}}{Pe_{h}^{2}}}+(\boldsymbol{U}^{\ast }\boldsymbol{\cdot }\boldsymbol{k}^{\ast })^{2}\right]\left[1+{\displaystyle \frac{Re_{h}^{2}(\boldsymbol{U}^{\ast }\boldsymbol{\cdot }\boldsymbol{k}^{\ast })^{2}}{(2\unicode[STIX]{x03C0})^{2}k^{\ast 4}}}\right]},\end{eqnarray}$$
with
$\boldsymbol{U}^{\ast }=\boldsymbol{U}/U$
,
$\boldsymbol{k}^{\ast }=\boldsymbol{k}h$
,
$Re_{h}=\unicode[STIX]{x1D70C}_{c}Uh/\unicode[STIX]{x1D707}_{c}$
and
$Pe_{h}=Uh/D$
. The solution given by Koch et al. (Reference Koch, Cox, Brenner and Brady1989) (equation (4.5) therein) for rigid spheres and Stokes flow is recovered in the limit
$Re\rightarrow 0$
and
$\unicode[STIX]{x1D707}_{d}/\unicode[STIX]{x1D707}_{c}\rightarrow \infty$
.
The tensor
$\unicode[STIX]{x1D63E}$
only depends on
$Pe_{h}$
,
$Re_{h}$
and the orientation of
$\boldsymbol{U}$
relative to the reciprocal lattice (the structure of which is, for cubic arrays, identical to that of the direct lattice). As highlighted by Koch et al. (Reference Koch, Cox, Brenner and Brady1989), the asymptotic behaviour of
$\unicode[STIX]{x1D63E}$
, and hence of
$\unicode[STIX]{x1D63F}^{\,conv}$
, depends on whether there exists any
$\boldsymbol{k}$
such that
$\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{k}=0$
, that is, on whether there exists any separation vector
$\boldsymbol{r}_{\boldsymbol{n}}$
in real space that is perpendicular to
$\boldsymbol{U}$
. The asymptotic behaviour of
$\Vert \unicode[STIX]{x1D63F}^{\,conv}\Vert$
, where
$\Vert ~\Vert$
denotes the tensorial Frobenius norm, is provided in table 2. The results show that the dependence of
$\Vert \unicode[STIX]{x1D63F}^{\,conv}\Vert$
on
$Pe$
in the limits
$Pe_{h}\ll 1$
and
$Pe_{h}\gg 1$
is, qualitatively, not affected by (weak) inertial effects.
Asymptotic order of
$\Vert \unicode[STIX]{x1D63F}^{\,conv}\Vert$
depending on
$Pe_{h}$
,
$Re_{h}$
and the orientation of the mean flow relative to the real lattice, based on the solution (4.8), derived for an ordered array of point particulates in Oseen flow conditions (
$F$
is the Oseen drag divided by the Stokes drag).
4.1.2 Application to ordered arrays rising vertically
Let us now come back to our original problem of an ordered array of particulates rising under the effect of buoyancy. The gravitational acceleration is oriented along a primary axis of the array,
$\boldsymbol{g}=-g\boldsymbol{e}_{3}$
, and although this is not the only possible solution (see e.g. Loisy et al.
Reference Loisy, Naso and Spelt2017), we restrict the analysis to the simplest case of bubbles rising vertically. In this case the hydrodynamic force exerted by the fluid on a particulate is parallel to the drift velocity, and, since this force balances the buoyancy force at steady state,
$F$
is related to
$U$
through
$$\begin{eqnarray}F=\frac{U_{0,Stokes}}{U},\end{eqnarray}$$
where
$U_{0,Stokes}$
is the terminal velocity of an isolated spherical fluid particulate in Stokes flow,
$$\begin{eqnarray}U_{0,Stokes}=\frac{1}{12}\frac{|\unicode[STIX]{x1D70C}_{c}-\unicode[STIX]{x1D70C}_{d}|gd_{b}^{2}}{\unicode[STIX]{x1D707}^{\ast }\unicode[STIX]{x1D707}_{c}},\quad \text{with }\unicode[STIX]{x1D707}^{\ast }=\frac{\unicode[STIX]{x1D707}_{c}+3\unicode[STIX]{x1D707}_{d}/2}{\unicode[STIX]{x1D707}_{c}+\unicode[STIX]{x1D707}_{d}}.\end{eqnarray}$$
Note that
$F$
can also be expressed in terms of commonly employed dimensionless groups:
$$\begin{eqnarray}F=\frac{1}{12\unicode[STIX]{x1D707}^{\ast }}\frac{Ar^{2}}{Re}.\end{eqnarray}$$
In the ‘sedimentation’ problem considered here,
$F$
is generally not known (as
$U$
is generally not known): it is a non-trivial function of the flow regime and volume fraction, which reduces to (4.7) when
$\unicode[STIX]{x1D719}\rightarrow 0$
and when the Oseen flow approximation is applicable.
The longitudinal and transverse components of the convective contribution,
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
and
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
, respectively, have been calculated from (4.8) for
$d_{b}/h=10^{-6}$
as a function of
$Pe_{h}$
for various
$Re<1$
. This very low value of
$d_{b}/h$
is required to allow
$Pe_{h}\gg 1$
while satisfying the condition
$Pe=Pe_{h}\,d_{b}/h\ll 1$
under which the analytical solution has been derived. The results, shown in figure 1, indicate that the asymptotic dependences of
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
and
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
on
$Pe$
are independent of
$Re$
. The sole effect of inertia is to modify the proportionality constants (by a substantial amount for the transverse component, though).
Longitudinal (a) and transverse (b) components of
$\unicode[STIX]{x1D63F}^{\,conv}$
as a function of the Péclet number based on the lattice spacing (
$Pe_{h}=Uh/D$
) for ordered arrays of point particulates at various small but finite Reynolds numbers (
$\boldsymbol{U}=U\boldsymbol{e}_{3}$
,
$d_{b}/h=10^{-6}$
and
$F$
is given by (4.9)). Note that, in (a),
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
is compensated by
$Pe^{2}$
.
In the limit of low
$Pe_{h}$
(say,
$Pe_{h}<10^{1}$
), both the transverse and the longitudinal components of
$\unicode[STIX]{x1D63F}^{\,conv}$
exhibit a quadratic dependence on the Péclet number (
$\unicode[STIX]{x1D60B}_{\bot ,\Vert }^{\,conv}\propto D\,Pe^{2}$
). In this regime, diffusion is much faster than convection. As the scalar is advected by velocity disturbances, it rapidly spreads out owing to diffusion, and convective dispersion (measured through
$\unicode[STIX]{x1D63F}^{\,conv}$
) is influenced by both mechanisms. This regime corresponds to the ‘convectively enhanced dispersion’ regime in Koch et al. (Reference Koch, Cox, Brenner and Brady1989).
In the limit of high
$Pe_{h}$
(say,
$Pe_{h}>10^{3}$
), the transverse component of
$\unicode[STIX]{x1D63F}^{\,conv}$
is independent of the Péclet number (
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}\propto D$
) whereas its longitudinal component grows quadratically with the Péclet number (
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}\propto D\,Pe^{2}$
). In this regime, convection dominates, but owing to the spatial periodicity of the flow, convective dispersion is obtained only if molecular diffusion across streamlines is considered (Koch et al.
Reference Koch, Cox, Brenner and Brady1989). This regime is termed ‘Taylor dispersion’ owing to the formal analogy, pointed out by Brenner (Reference Brenner1980), with one-dimensional shear-induced Taylor dispersion in a capillary tube.
We emphasize that the expression (4.8) has been derived from the approximation (4.2), the validity of which is established only for
$Pe\ll 1$
(which is, in practice, of limited use). Using symmetry arguments, Koch et al. (Reference Koch, Cox, Brenner and Brady1989, § 4.2) showed that, in the limit
$Pe\gg 1$
, Taylor dispersion is obtained if the average flow is perpendicular to a set of planes of both translational and reflectional symmetry, such as Stokes flows parallel to the primary axis of an ordered array of spheres. Taylor dispersion is then easily understood by remarking that, owing to the symmetries of the flow, a fluid tracer particle entering the unit cell at one point, say
$\boldsymbol{x}$
, exits the cell at the equivalent point in the next cell, that is,
$\boldsymbol{x}+h\boldsymbol{e}_{3}$
, so that dispersion can only occur if diffusion across streamlines is present (Koch et al.
Reference Koch, Cox, Brenner and Brady1989). In the presence of inertial effects, the reflectional symmetry is lost, hence this argument does not hold. Koch et al. (Reference Koch, Cox, Brenner and Brady1989) also demonstrated that, for Stokes flow, the solution for
$Pe\ll h/d_{b}$
is identical, at lowest order, to that obtained for
$Pe\ll 1$
(§ 4.3 therein, note that their
$Pe$
corresponds to
$Pe_{h}$
in our notations). Such a demonstration for Oseen flow will not be attempted here. Instead, the range
$Pe\geqslant 1$
will be explored using DNS.
4.2 Numerical results
The above analysis provides explicit expressions for
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
and
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
. These are valid for spherical bubbles rising at
$Re<1$
(strictly speaking, at a Reynolds number sufficiently small to assume Oseen flow, in terms of Archimedes and Bond numbers, this regime would be reached for
$Bo<1$
and
$Ar\lesssim 1$
) and in the limits
$\unicode[STIX]{x1D719}\rightarrow 0$
and
$Pe\ll 1$
. We shall now determine using numerical simulations whether these restrictions can be relaxed, and, if so, to what extent.
Longitudinal (a) and transverse (b) components of
$\unicode[STIX]{x1D63F}^{\,conv}$
as a function of the Péclet number based on the lattice spacing (
$Pe_{h}=Uh/D$
) for ordered arrays in various flow regimes at small volume fraction (
$\unicode[STIX]{x1D719}=0.2\,\%$
). The normalizations of
$\unicode[STIX]{x1D60B}_{\Vert ,\bot }^{\,conv}$
are those suggested by the asymptotic analysis (identical to those used in figure 1), and
$F$
is given by (4.9). The lines are drawn to guide the eye. Note that, in (a),
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
is compensated by
$Pe^{2}$
.
We examine the case of suspensions at low (but not vanishing) volume fraction in order to approach the dilute limit assumption, and to focus on the sole effect of inertia. The longitudinal and transverse components of the convective contribution to the effective diffusivity have been computed for
$h/d_{b}=6.4$
, which corresponds to a gas volume fraction of
$\unicode[STIX]{x1D719}=0.2\,\%$
(the smallest volume fraction accessible with the method and facilities used), for each of the four flow regimes listed in table 1, and
$Pe$
has been varied from
$10^{-1}$
to
$10^{3}$
. The results are shown in figure 2 as a function of
$Pe_{h}$
, the Péclet number based on the lattice spacing, which is the parameter governing the transition between the two asymptotic limits (see § 4.1 and figure 1). The different colours, symbols and line styles depict the different flow regimes (the lines are drawn to guide the eye). Qualitatively, figure 2 bears a striking resemblance to figure 1, even for case C (skirted bubbles): analysis and simulations yield similar dependences of
$\unicode[STIX]{x1D60B}_{\Vert ,\bot }^{\,conv}$
on
$Pe$
and qualitatively comparable effects of increasing
$Re$
. At low Péclet number (
$Pe_{h}\lesssim 10^{1}$
), dispersion occurs primarily by molecular diffusion and convective mixing grows quadratically with
$Pe$
in both the longitudinal and transverse directions (
$\unicode[STIX]{x1D60B}_{\Vert ,\bot }^{\,conv}\propto D\,Pe^{2}$
). At high Péclet number (
$Pe_{h}\gtrsim 10^{3}$
), Taylor dispersion is the dominant process, with very efficient mixing in the flow direction (
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}\propto D\,Pe^{2}$
) and negligible mixing in the transverse one (
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}\propto D$
). Inertial effects and bubble deformation only affect the proportionality constants, rather weakly for
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
but substantially for
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
.
Numerical solution
$\unicode[STIX]{x1D63F}^{\,conv}$
divided by the analytical solution
$\unicode[STIX]{x1D63F}^{conv,anal}$
as a function of the Péclet number based on the bubble diameter (
$Pe=Ud_{b}/D$
) for ordered arrays in various flow regimes at small volume fraction (
$\unicode[STIX]{x1D719}=0.2\,\%$
): longitudinal (a) and transverse (b) components. Diffusivity
$\unicode[STIX]{x1D63F}^{conv,anal}$
is given by (4.8).
Scalar fluctuation field
$c^{\prime }$
associated with
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
, shown in a vertical symmetry plane passing through the centre of a bubble, for ordered arrays in various flow regimes at
$Pe=10^{-1}$
(a) and
$Pe=10^{3}$
(b). The imposed scalar field
$\bar{c}$
increases linearly within the cell from bottom to top (
$\unicode[STIX]{x1D719}=0.2\,\%$
, the entire cell is shown, and gravity is pointing downward).
Scalar fluctuation field
$c^{\prime }$
associated with
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
, shown in a vertical symmetry plane passing through the centre of a bubble, for ordered arrays in various flow regimes at
$Pe=10^{-1}$
(a) and
$Pe=10^{3}$
(b). The imposed scalar field
$\bar{c}$
increases linearly within the cell from left to right (
$\unicode[STIX]{x1D719}=0.2\,\%$
, the entire cell is shown, and gravity is pointing downward).
To allow a quantitative comparison between the DNS and the analysis, we present in figure 3 the ratio of
$\unicode[STIX]{x1D60B}_{\Vert ,\bot }^{\,conv}$
to
$\unicode[STIX]{x1D60B}_{\Vert ,\bot }^{conv,anal}$
, where
$\unicode[STIX]{x1D60B}_{\Vert ,\bot }^{conv,anal}$
is given by (4.8) with
$F$
computed directly from its definition (4.9). As the range of validity of the analysis is defined in terms of
$Pe$
(
$Pe\ll 1$
, with
$Pe$
the Péclet number based on the bubble diameter), the data are presented here as a function of
$Pe$
rather than
$Pe_{h}$
. For the longitudinal component, the numerical solution does not deviate by more than 5 % from the theoretical prediction, as can be seen from figure 3(a). The fact that the low-
$Pe$
, Oseen flow analysis yields accurate predictions for
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
at
$Pe=10^{3}$
and
$Re=O(10)$
is not surprising, as the behaviour of
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}/(DF^{2}Pe^{2})$
is rather insensitive to both the flow regime and the Péclet number (as shown in figures 1(a) and 2(a), this quantity does not vary more than 15 % for the cases studied). We conclude that, at small volume fraction,
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
can be predicted within
$\pm 5\,\%$
from (4.8) at any Péclet number up to
$10^{3}$
and any Reynolds number up to 40, even when the bubbles are strongly deformed. For the transverse component, the asymptotic analysis underpredicts the value of
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
at high Péclet number, even for
$Re\lesssim 1$
. As a consequence,
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
cannot be accurately estimated from our analytical solution when the assumptions underlying its derivation are not satisfied. It must be kept in mind though that this component varies much more than the longitudinal one between the regimes of small and large Péclet numbers, and is much more sensitive to the flow regime (
$Re$
, shape), which means that its value is more difficult to predict. In all, it is worth stressing that the asymptotic analysis yields the correct qualitative behaviour and order of magnitude for
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
at least up to
$Pe=10^{3}$
and
$Re\approx 10$
, even for strongly deformed bubbles. Finally, we emphasize that we found
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}/\unicode[STIX]{x1D60B}_{\bot }^{\,conv}\gtrsim 10^{2}$
, so the most important component of the effective diffusivity tensor is the longitudinal one, except in situations where there is no longitudinal component of the gradient of the scalar on the macroscale.
To illustrate the dispersion regimes at low and high Péclet number, we present in figures 4 and 5 visualizations of the scalar fluctuation field
$c^{\prime }$
used to compute
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
and
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
, respectively. In each of these figures, the field of
$c^{\prime }$
is represented for each flow regime in a vertical symmetry plane passing through the centre of a bubble for
$Pe=10^{-1}$
(a) and
$Pe=10^{3}$
(b), and the Reynolds number increases from top to bottom. The field of
$c^{\prime }$
associated with
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
, shown in figure 4, exhibits similar features at low and high
$Pe$
. In contrast, the field of
$c^{\prime }$
associated with
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
, represented in figure 5, is qualitatively different in these two limits. This illustrates qualitatively why the regimes at low and high
$Pe$
are similar for
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
(
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}\propto Pe^{2}$
), whereas the scaling laws identified for
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
are different in both limits (see figure 2). In addition, the Reynolds number and the bubble shape affect the fore-and-aft symmetry and the details of
$c^{\prime }$
, but not its essential features, which results in quantitative but not qualitative effects on
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
and
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
.
5 Freely evolving suspensions
We examine in this section scalar mixing in freely evolving suspensions as represented by the periodic repetition of a unit cell containing several independent bubbles (‘free arrays’). Our objectives here are threefold: (i) to investigate the effective diffusivity of freely evolving suspensions at small and high Péclet numbers; (ii) to compare and contrast these results with those obtained in ordered systems; and (iii) to evaluate the effect of the system size (number of bubbles in a unit cell,
$N_{b}$
).
For that purpose, we considered a single flow regime (ellipsoidal bubbles at
$Re=O(10)$
, corresponding to case E1 in table 1) at intermediate volume fraction (
$\unicode[STIX]{x1D719}=2.4\,\%$
) and explored the effect of varying the number of free bubbles
$N_{b}$
on the dependence of
$\unicode[STIX]{x1D63F}^{\,conv}$
on the Péclet number. Owing to the multiplicity of simulations involved and to their duration (typically several months on 64 cores), only a few different values of
$N_{b}$
belonging to a rather limited range have been considered (namely
$N_{b}=\{2,3,5,8,12\}$
in the simulations for the determination of
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
, and
$N_{b}=\{2,8\}$
in those for
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
). For the same reason, investigations of the effects of volume fraction and flow regime could not be undertaken.
The longitudinal and transverse components of
$\unicode[STIX]{x1D63F}^{\,conv}$
are plotted in figure 6 as a function of the Péclet number for various values of
$N_{b}$
. Note that a very wide range of Péclet numbers is considered. Convergence of
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
with the system size is very fast: the values of
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
are essentially independent of the number of free bubbles for
$2\leqslant N_{b}\leqslant 12$
at all Péclet numbers. This suggests that
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
is independent of the system size
$N_{b}$
, although this would need to be confirmed by considering larger values of this parameter. Our data for
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
suggest that convergence with
$N_{b}$
is slower for this quantity, especially at high Péclet number, although conclusions can hardly be drawn on this point due to the few values of
$N_{b}$
considered.
Longitudinal (a) and transverse (b) components of
$\unicode[STIX]{x1D63F}^{\,conv}$
as a function of the Péclet number for various numbers of free bubbles
$N_{b}$
in the unit cell (
$N_{b}=1$
corresponds to an ordered array). Symbols other than purple stars, DNS (
$Re\approx 30$
,
$\unicode[STIX]{x1D719}=2.4\,\%$
); purple stars, experimental data of Alméras et al. (Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015) (
$Re\approx 700$
,
$\unicode[STIX]{x1D719}\approx 2.4\,\%$
). A spatial resolution of
$d_{b}/\unicode[STIX]{x0394}x=20$
was used for
$N_{b}>1$
; the effect of increasing resolution to
$d_{b}/\unicode[STIX]{x0394}x=30$
is illustrated by the filled red squares for
$N_{b}=8$
and
$Pe\approx 10^{3}$
(
$d_{b}$
is the bubble volume-equivalent diameter and
$\unicode[STIX]{x0394}x$
is the grid spacing).
Instantaneous scalar fluctuation field
$c^{\prime }$
associated with
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
for a free array of eight bubbles, at
$Pe=10^{-1}$
(a) and
$Pe=10^{6}$
(b). The gradient of
$\bar{c}$
is vertical (the entire cell is shown, and gravity is pointing downward).
Instantaneous scalar fluctuation field
$c^{\prime }$
associated with
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
for a free array of eight bubbles, at
$Pe=10^{-1}$
(a) and
$Pe=10^{6}$
(b). The gradient of
$\bar{c}$
is horizontal (the entire cell is shown, and gravity is pointing downward).
We first examine the dependence of
$\unicode[STIX]{x1D63F}^{\,conv}$
on the Péclet number in free arrays of bubbles (
$N_{b}>1$
). At small
$Pe$
,
$\unicode[STIX]{x1D60B}_{\Vert ,\bot }^{\,conv}\propto D\,Pe^{2}$
, whereas at high
$Pe$
,
$\unicode[STIX]{x1D60B}_{\Vert ,\bot }^{\,conv}\propto D\,Pe=Ud_{b}$
. Note that the scaling at high
$Pe$
is expected from a simple dimensional analysis in a convection-dominated regime where diffusion plays no role. This regime corresponds to the ‘mechanical dispersion’ regime in Koch & Brady (Reference Koch and Brady1985). The different dispersion regimes at low and high
$Pe$
can also be identified from the features of the scalar fluctuation field
$c^{\prime }$
. Instantaneous snapshots of
$c^{\prime }$
associated with
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
and
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
are shown in figures 7 and 8, respectively, for an array of eight free bubbles at
$Pe=10^{-1}$
(a) and at
$Pe=10^{6}$
(b). For a given component, the isocontours of
$c^{\prime }$
follow markedly different patterns at low and high
$Pe$
.
We now compare these results with those obtained for ordered arrays (black crosses in figure 6) and discuss the effect of the microstructure. At small
$Pe$
,
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
and
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
grow quadratically with
$Pe$
in both free and ordered arrays. This scaling was also obtained by Koch & Brady (Reference Koch and Brady1985) for low-
$Pe$
dispersion in porous media with random microstructure (albeit in the Stokes flow limit). Since, in the low-
$Pe$
regime, diffusion by the random motion of molecules is much faster than convection by the flow, the microstructure has only a quantitative incidence on
$\unicode[STIX]{x1D63F}^{\,conv}$
, and dispersion is qualitatively identical in ordered and freely evolving suspensions. Note that similar features in the spatial distribution of
$c^{\prime }$
can be identified in ordered and free arrays at low
$Pe$
(see tubular structures on figures 4(a) and 7(a) for
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
, and quadrupolar ones on figures 5(a) and 8(a) for
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
). We emphasize, however, that precise quantitative agreement between the results for one and for many bubbles at low Péclet number in figure 6 is not expected, as the flows and the microstructures in the two systems are different (Bunner & Tryggvason Reference Bunner and Tryggvason2002; Loisy et al.
Reference Loisy, Naso and Spelt2017).
At high
$Pe$
, the Taylor dispersion scaling obtained in ordered arrays is replaced, in both directions, by a scaling similar to the one characterizing mechanical dispersion, as soon as the relative motion between bubbles is allowed. In this regime, the transverse dispersion is indeed governed by mechanical dispersion. Irrespective of the value of
$Pe$
, any Taylor dispersion in the vertical direction is limited by transverse diffusion or dispersion, the latter becoming more significant at large
$Pe$
. This results in a scaling similar to that of mechanical dispersion in the longitudinal direction as well, such that a distinction between these two mechanisms (pure mechanical dispersion, or Taylor dispersion limited by transverse mechanical one) cannot be made. Incidentally, mechanical dispersion is also obtained at high
$Pe$
in random media in Stokes flow conditions (Koch & Brady Reference Koch and Brady1985). Although the microstructure of the present bubbly suspensions has not been evaluated quantitatively, visual inspection and prior results on their dynamics (Loisy et al.
Reference Loisy, Naso and Spelt2017) showed that it is not random, but rather characterized by a certain ‘organization’. Despite the fact that freely evolving suspensions resemble ordered ones with respect to their dynamics, scalar dispersion is extremely sensitive to the presence of disorder, and is fundamentally different in perfectly ordered and weakly disordered suspensions at high Péclet number. It does not, however, seem to be sensitive to the degree of disorder, as suggested by the fact that the same scalings with
$Pe$
are obtained for random porous media and weakly disordered suspensions. We stress that this last statement is purely speculative, and would require a quantitative study of the effect of the microstructure to be confirmed.
We finally attempt a comparison of our results with the experimental data of Alméras et al. (Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015), who measured the effective diffusivity of a homogeneous swarm of high-Reynolds-number rising bubbles at
$Pe\approx 1.75\times 10^{6}$
for gas volume fractions ranging from 1 % to 13 %. It is important to stress that, in these experiments,
$Re\approx 700$
, whereas in the simulations,
$Re\approx 30$
, so the comparison is only indicative. Interpolation (by eye) of their data at
$\unicode[STIX]{x1D719}\approx 2.4\,\%$
(figure 10 in their paper) yields
$\unicode[STIX]{x1D60B}_{\Vert }^{\,eff}/D=1\times 10^{6}$
and
$\unicode[STIX]{x1D60B}_{\bot }^{\,eff}/D=5\times 10^{5}$
. These experimental values are represented by purple stars in figure 6. Note that, at such high Péclet number, the dominant contribution to
$\unicode[STIX]{x1D63F}^{\,eff}$
is due to
$\unicode[STIX]{x1D63F}^{\,conv}$
, so it seems reasonable to assume that these are equivalent. The orders of magnitude of
$\unicode[STIX]{x1D60B}_{\Vert }^{\,eff}/D$
are comparable in the experiment and in the simulation, whereas
$\unicode[STIX]{x1D60B}_{\bot }^{\,eff}/D$
is much higher in the experiment. This difference can be explained from the different properties of the numerical and experimental flows considered: partition coefficient (the dye concentration in the gas is presumably zero in the experiments from Alméras et al. (Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015)), diffusivity ratio and bubble-induced liquid agitation in the horizontal direction. In our simulations of free arrays at moderate
$Re$
, the bubbles were indeed observed to rise along nearly straight vertical lines, and the anisotropy ratio characterizing the liquid velocity variance,
$2\langle u_{3}^{\prime }u_{3}^{\prime }\rangle /\langle u_{1}^{\prime }u_{1}^{\prime }+u_{2}^{\prime }u_{2}^{\prime }\rangle$
, is approximately 8 (for
$N_{b}=8$
), whereas in the experiment at high
$Re$
, the bubble motion is fully three-dimensional, and the anisotropy ratio is approximately 2. Finally, as only one value of the Péclet number was considered in the experiments of Alméras et al. (Reference Alméras, Risso, Roig, Cazin, Plais and Augier2015), no comparison of their data with our results can be offered regarding the dependence of the effective diffusivity on the Péclet number.
6 Conclusions
In this study we investigated scalar dispersion in homogeneous bubbly suspensions as described by an effective diffusivity tensor. The longitudinal and transverse components of the convective contribution to the effective diffusivity, denoted
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
and
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
, respectively, have been computed for bubbly suspensions in various flow regimes. This convective contribution is that associated with bubble-induced agitation, and is the dominant contribution to the effective diffusivity in commonly encountered bubbly flows.
The dispersion theory of Koch et al. (Reference Koch, Cox, Brenner and Brady1989) indicates that convective mixing mechanisms in ordered suspensions in Stokes flow conditions differ at low and high Péclet numbers. According to this theory, when the bulk flow is aligned with a primary axis of a simple cubic lattice of spheres, convectively enhanced dispersion is expected at low Péclet number, whereas Taylor dispersion should dominate at high Péclet number. In the present study, we have extended this theory to account for weak inertial effects, and we have shown that these two dispersion regimes are qualitatively unchanged in the presence of (weak) inertia. This result has been confirmed by DNS for values of the Reynolds number ranging from vanishingly small to moderate. In all investigated cases,
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
was found to be significantly larger than
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
, and theoretical predictions have been shown to yield the correct qualitative behaviour and order of magnitude of both
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
and
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
in a variety of flow regimes (spherical to strongly deformed bubbles with Reynolds numbers up to 10) at small volume fraction.
DNS of scalar transport in freely evolving bubbly suspensions, as represented by free arrays of bubbles, have been carried out for a wide range of Péclet numbers, and the effect of introducing additional degrees of freedom in the system has been evaluated. At low Péclet number, dispersion in free arrays is convectively enhanced, as in ordered ones. At high Péclet number, in freely evolving suspensions wherein at least two bubbles are present in a unit cell, the longitudinal component of the effective diffusivity exhibits a scaling that is similar to that characterizing mechanical dispersion. This suggests that the limiting role of molecular diffusion to Taylor dispersion is taken over by mechanical dispersion, or that mechanical dispersion itself dominates. Besides, the effective diffusivity seems to be weakly sensitive to the number of bubbles present in a unit cell. This last assertion requires more thorough investigations to be confirmed, but is encouraging regarding the possibility of computing the effective diffusivity of homogeneous bubbly flows from DNS of systems of relatively small size. This would allow in particular a thorough investigation of the roles played by the volume fraction and the flow regime, which could not be undertaken as part of the present study.
The results presented in this paper are restricted to bubbles having the same diffusivity as that of the surrounding liquid, and to scalar fields that are continuous across the interface, and therefore cannot be straightforwardly compared to those obtained in real bubbly flows. A jump in the scalar field, which represents the difference in solubilities given by Henry’s law in the context of chemical species transport, as well as a difference in diffusivities, would introduce a diffusive contribution to the effective diffusivity tensor (2.5) in addition to the convective one considered in this study. The present results show the convective contribution at large Péclet numbers and modest volume fraction to be substantially larger than the diffusive contribution from non-equal diffusivities or solubilities (Maxwell Reference Maxwell1873; Jeffrey Reference Jeffrey1973; Koch & Brady Reference Koch and Brady1985). A difference in diffusivities or solubilities would, however, also have some indirect effect on the convective contribution, the magnitude of which should be investigated in the future.
Besides the effective diffusivity, another quantity of practical importance is the rate of interfacial scalar transport in the presence of an average scalar gradient between the disperse phase and the bulk. Heat and mass exchanges across phase boundaries are traditionally expressed as dimensionless transfer coefficients called the Nusselt and Sherwood numbers, respectively. Their functional dependences on suspension properties, in particular the volume fraction, have been the subject of analytical (Acrivos, Hinch & Jeffrey Reference Acrivos, Hinch and Jeffrey1980), numerical (Aboulhasanzadeh & Tryggvason Reference Aboulhasanzadeh and Tryggvason2014) and experimental (Colombet et al. Reference Colombet, Legendre, Cockx, Guiraud, Risso, Daniel and Galinat2011, Reference Colombet, Legendre, Risso, Cockx and Guiraud2015) studies. Formally, the Nusselt and Sherwood numbers are closure coefficients for the conditionally averaged scalar transport equation, where the conditional average is defined as an ensemble average over the subset of realizations wherein a particulate is present at a given position. Less formally, the Nusselt and Sherwood numbers are related to a ‘mesoscale’ description of scalar transfer between the two phases, whereas the effective diffusivity is associated with a ‘macroscale’ description of scalar transport through a two-phase mixture seen as a continuum. They correspond to different closure problems, and one cannot be inferred from the other. Nevertheless, the present work will be primarily important for mass transfer processes in bubbly flows that are liquid-phase controlled. This is because then the mixture concentration distribution is key, whereas if it is gas-phase controlled, the concentration in the liquid will be almost uniform and one is primarily concerned by the circumstances inside each bubble.
Acknowledgements
This work benefited from the financial support of the French Research Agency (grant ANR-12-BS09-0011), and was performed using the HPC resources provided by GENCI-CINES and GENCI-IDRIS (grant x20162b6893), PSMN (École Normale Supérieure de Lyon), P2CHPD (Université Claude Bernard Lyon 1) and PMCS2I (École Centrale de Lyon).
Appendix. Spatial convergence tests
Spatial convergence for an ordered array of bubbles in case E1 at
$Pe=10^{3}$
: relative error in
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
and
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
as a function of the grid spacing
$\unicode[STIX]{x0394}x$
(
$d_{b}$
is the bubble volume-equivalent diameter;
$\unicode[STIX]{x1D60B}_{\unicode[STIX]{x0394}x=0}^{\,conv}$
is extrapolated assuming
$\unicode[STIX]{x1D60B}^{\,conv}=\unicode[STIX]{x1D60B}_{\unicode[STIX]{x0394}x=0}^{\,conv}-k\unicode[STIX]{x0394}x^{n}$
, where the values of the three parameters
$\unicode[STIX]{x1D60B}_{\unicode[STIX]{x0394}x=0}^{\,conv}$
,
$k$
and
$n$
are fitted from numerical data).
We present the results of some spatial convergence tests of the algorithm solving the scalar transport equation. The results of similar tests for the algorithm solving the flow are shown in Loisy (Reference Loisy2016) and Loisy et al. (Reference Loisy, Naso and Spelt2017).
The effect of the grid spacing on
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
and
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
has been assessed for case E1 at
$Pe=10^{3}$
for one value of the volume fraction (
$\unicode[STIX]{x1D719}=2.4\,\%$
), in both ordered and free configurations. For ordered arrays, three different resolutions were tested, namely
$d_{b}/\unicode[STIX]{x0394}x=\{20,40,60\}$
with
$\unicode[STIX]{x0394}x$
the grid spacing. The results are shown in figure 9. The error in the values of
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
and
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
arising from spatial discretization is less than 1 % when a resolution of 40 grid cells per bubble diameter is used. This resolution is the same as that used for the simulation of the corresponding bubbly flow in Loisy et al. (Reference Loisy, Naso and Spelt2017). In practice, we used for each configuration the same resolution as that selected for the simulation of the corresponding ordered bubbly suspensions (see Loisy et al.
Reference Loisy, Naso and Spelt2017), namely 60 grid cells per diameter for case C and 40 grid cells per diameter for the other cases.
For free arrays, due to the computational cost of the simulations, only two different resolutions were tested, namely 20 and 30 grid cells per bubble diameter, for an array of eight bubbles. Simulations at higher resolution were too expensive to be continued over sufficiently long times to allow a quantitative estimate of the uncertainty. However, the values of
$\unicode[STIX]{x1D60B}_{\Vert }^{\,conv}$
and
$\unicode[STIX]{x1D60B}_{\bot }^{\,conv}$
obtained with the finer grid, depicted by filled red squares in figure 6, are nearly indistinguishable from those obtained with the coarser grid. A resolution of 20 grid cells per diameter was therefore concluded to be sufficient for free arrays in view of the present purposes.
For a given case, the same resolution was used for all Péclet numbers. Note that, when the gas diffusivity differs from that of the liquid (a situation not considered here but frequently encountered in practice), finer resolutions may be required, as thin scalar boundary layers around the bubbles would then need to be resolved.
