2.1. Continuous phase

Assuming that the continuous phase is an incompressible Newtonian fluid, its dynamics is described by the Navier–Stokes equations,

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

(2.2)$$\begin{gather}\frac{\partial\boldsymbol{u}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} (\boldsymbol{u} \boldsymbol{u}) ={-}\frac{1}{\rho}\boldsymbol{\nabla} P + \boldsymbol{\nabla}\boldsymbol{\cdot}\left( \nu \left[ \boldsymbol{\nabla}\boldsymbol u + \left(\boldsymbol{\nabla}\boldsymbol u \right)^{\textrm{T}} \right]\right) + \boldsymbol{g}, \end{gather}$$

where $\boldsymbol {u}$ and $p$ are the velocity and pressure fields of the continuous phase, $\rho$ and $\nu$ its density and kinematic viscosity, and $\boldsymbol {g}$ the gravity acceleration. Nevertheless, these equations are to be complemented with very complex boundary conditions that evolve as the bubbles move in response to the fluid stress applied on their interface. When the number of dispersed objects remains limited, one can rely on fully resolved numerical simulations that accurately track the interfaces (Tenneti & Subramaniam Reference Tenneti and Subramaniam2014; Maxey Reference Maxey2017). To date, the number of objects that can be simulated is around 1000 (Lucci, Ferrante & Elghobashi Reference Lucci, Ferrante and Elghobashi2010; Uhlmann & Chouippe Reference Uhlmann and Chouippe2017; Hasslberger, Klein & Chakraborty Reference Hasslberger, Klein and Chakraborty2018; Morente, Laviéville & Legendre Reference Morente, Laviéville and Legendre2018; du Cluzeau, Bois & Toutant Reference du Cluzeau, Bois and Toutant2019; Lai et al. Reference Lai, Fraga, Chan and Dodd2019; du Cluzeau et al. Reference du Cluzeau, Bois, Toutant and Martinez2020; Pandey, Ramadugu & Perlekar Reference Pandey, Ramadugu and Perlekar2020), which, for a given number density of the dispersed phase, automatically limits the size of the fluid domain that can be simulated. However, the flow being set by the interactions of multiple wakes, obtaining a fully developed flow requires a large number of bubbles and a large-size domain.

In the point-particle approach, a coarse-grain description of the flow is adopted and one gives up the fine resolutions of the flow in the bubble vicinity. One can define the local smoothed dispersed-phase indicator function as

(2.3)\begin{equation} \chi_b{(\boldsymbol{x},t)} = \sum_b \mathcal{V}_b G_{\sigma}{(\boldsymbol{x}-\boldsymbol{x}_b(t))}, \end{equation}

where $G_{\sigma }$ is a homogenization kernel of characteristic size $\sigma$ with $\int G_{\sigma }(\boldsymbol {x})\,{\rm d}{\kern 0.2em}x^{3}=1$, $\mathcal {V}_b$ is the volume of a bubble and $\boldsymbol {x}_b$ is the position of the $b$th bubble. The coarse graining results in a great simplification of the boundary conditions, since the dispersed phase simply appears as a volume source of momentum in the Navier–Stokes equations,

(2.4)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}_f = 0, \end{gather}$$

(2.5)$$\begin{gather}\frac{D \boldsymbol{u}_f}{Dt} ={-} \frac{1}{\rho_f} \boldsymbol{\nabla} p + \nu {\rm \Delta} \boldsymbol{u}_f + \boldsymbol{f}, \end{gather}$$

where $\boldsymbol {u}_f$ is the coarse-grained velocity field and $D_t = \partial _t + \boldsymbol {u}_f\boldsymbol {\cdot }\boldsymbol {\nabla }$ is the material derivative involving advection velocity $\boldsymbol {u}_f$. The momentum source field is expressed as

(2.6)\begin{equation} \boldsymbol{f}{(\boldsymbol{x},t)} = \frac{1}{\rho_f}\sum_b \boldsymbol{F}_{b \rightarrow f} G_{\sigma}{(\boldsymbol{x}-\boldsymbol{x}_b(t))}, \end{equation}

where $\boldsymbol {F}_{b \rightarrow f}$ is the momentum source due to bubble $b$. In (2.5), $p$ is the pressure variation relative to the hydrostatic pressure $P_h$ with $-1/\rho _f \boldsymbol {\nabla } P_h = \langle \boldsymbol {f} \rangle$. Note that we have neglected the variations of the mixture density, as well as additional flux terms due to the unresolved scales in (2.4) and (2.5).

Typically, the kernel $G_{\sigma }$ can be considered Gaussian (Maxey et al. Reference Maxey, Patel, Chang and Wang1997; Zamansky et al. Reference Zamansky, Coletti, Massot and Mani2016),

(2.7)\begin{equation} G_{\sigma}(\boldsymbol{r}) = \frac{1}{(2 {\rm \pi}\sigma^{2})^{3/2}} \exp \left( - \frac{\boldsymbol{r}^{2}}{2\sigma^{2}} \right), \end{equation}

where the kernel width $\sigma$ remains as a free parameter. The question is therefore to relate $\sigma$ to the length scales of the problem, which are the bubble diameter $d_b$ and smallest scales $\eta$ of fluctuations that would exist in the vicinity of each dispersed object when the influence of this particular object is not considered. Note that this scale $\eta$ can result from the presence of the other objects constituting the dispersed phase or an imposed forcing. In addition, $\sigma$ must be larger than the mesh resolution ${\rm \Delta} x$ to ensure a sufficiently accurate numerical description of $G_{\sigma }$.

When the resolution is fine enough (${\rm \Delta} x < \eta$) and the dispersed objects are small (${d_b \ll \eta}$), the kernel width is generally based on the dissipative scale ($\sigma \approx \eta$) and the simulation is considered as a DNS of the carrier phase. It is worth mentioning that in this case, Gualtieri et al. (Reference Gualtieri, Picano, Sardina and Casciola2015) proposed to improve the accuracy further by introducing a time lag $\tau _{\sigma } \approx \sigma ^{2} /\nu$ in the calculation of the force to account for the diffusion of the momentum from scale $d_b$ to scale $\sigma$.

In contrast, when the bubbles are larger than the grid size ($d_b > {\rm \Delta} x$), it is relevant to choose $\sigma \approx d_b$. Of course, the velocity field in the vicinity of the bubbles is thus not well represented. Nevertheless, let us remark that it is also the case with small objects ($d_b \ll \eta$), since it is intrinsic to the point-particles approach and consistent with the coarse graining procedure. In any case, this is valid provided the effect of the dynamics of the flow in the vicinity of the dispersed objects on the dynamics of the larger scales is weak or can be modelled by means of an effective description.

2.2. Dispersed phase

The evolution of the dispersed phase is obtained by computing the trajectory of each individual dispersed body by means of Newton's equation,

(2.8a,b)\begin{equation} \frac{\textrm{d}{\kern 0.2em}\boldsymbol{x}_b(t)}{\textrm{d} t} = \boldsymbol{v}_b;\quad m_b \frac{\textrm{d} \boldsymbol{v}_b(t)}{\textrm{d} t} = m_b \boldsymbol{a}_b = \boldsymbol{F}_{{fluid},b} + \boldsymbol{F}_{{ext},b}+ \boldsymbol{F}_{{I},b}, \end{equation}

where $\boldsymbol {v}_b$ and $\boldsymbol {a}_b$ are the body velocity and acceleration, and $m_b$ is the mass of the body. The force applied on the body is decomposed into the hydrodynamic force exerted by the carrier fluid ($\boldsymbol {F}_{{fluid},b}$), the external volume force which reduces here to gravity ($\boldsymbol {F}_{{ext},b}$) and the force resulting from the direct interactions between bodies ($\boldsymbol {F}_{{I},b}$). The main problem of the Euler–Lagrange approach is to compute the hydrodynamic force, since the fluid stresses at the interface are not known because of the coarse graining approach. Formally, $\boldsymbol {F}_{{fluid},b}$ could be expressed as a functional depending on the coarse-grained fluid velocity field and the trajectory of all the bodies and their derivatives,

(2.9)\begin{equation} \boldsymbol{F}_{{fluid}, b} = \varPhi\left[ \boldsymbol{\boldsymbol{u}}_f(\boldsymbol{x},t) ; \boldsymbol{x}_1(t) \cdots \boldsymbol{x}_{N_b}(t) \right]. \end{equation}

Here, $\varPhi$ is a functional that relates the hydrodynamic force to the velocity field, which could in principle be obtained from the Navier–Stokes equation, but remains generally unknown. For instance, the force on a fixed object immersed in a flow that is uniform and stationary at infinity involves drag and lift components, which are described by theoretical or empirical relations. In the most general case, the orientation and rotation of the particles should also be accounted for in (2.9), but is disregarded here as it does not appear essential for a first attempt to simulate bubble-induced agitation.

When $\eta$ is large compared to the dispersed bodies ($\eta \gg d_b$), there exists an intermediate scale $\ell$ in between $\eta$ and $d_b$, at which the background flow around the objects can be considered as uniform. Based on this scale separation, the expression of the fluid force should only depend on the local value of the fluid velocity field and its derivatives described at scale $\ell$. Instead of introducing explicitly the scale $\ell$, one chooses to express the force from the value of the fluid velocity field at the bubble position, $\boldsymbol {u}_f{(\boldsymbol {x}=\boldsymbol {x}_b(t),t)}$, since the liquid velocity $\boldsymbol {u}_f$ obtained by (2.5) is defined everywhere, even in the region that is actually occupied by gas. However, because of the presence of the bubble, the liquid velocity close to the bubble differs significantly from its value at a distance of order $\ell$ from the bubble (Saffman Reference Saffman1973; Gatignol Reference Gatignol1983; Maxey & Riley Reference Maxey and Riley1983). Therefore the liquid velocity and its derivatives at the bubble position need to be corrected to reliably compute the hydrodynamic force. We note that $\tilde {\boldsymbol {u}}_{f,b}$ is the fluid velocity field corrected from the perturbation due to the presence of bubble $b$. We will detail in the next section how such a correction is obtained. In principle, one could express the hydrodynamic force as a functional of the corrected velocity field at the bubble position:

(2.10)\begin{equation} \boldsymbol{F}_{{fluid},b} = \phi \left[\boldsymbol{v}_b(t)-\tilde{\boldsymbol{u}}_{f,b} {(\boldsymbol{x}=\boldsymbol{x}_b(t),t)},\ \boldsymbol{\nabla} \tilde{\boldsymbol{u}}_{f,b}|_{\boldsymbol{x}_b},\ \partial_t \tilde{\boldsymbol{u}}_{f,b}|_{\boldsymbol{x}_b}, d_t\boldsymbol{v}_b(t) \right]. \end{equation}

Note that the hydrodynamic force must only involve velocity differences to preserve Galilean invariance. In addition to the scale separation, if the bubble Reynolds number is small enough, $Re_b = |\boldsymbol {v}_b{(t)}-\tilde {\boldsymbol {u}}_{f,b}{(\boldsymbol {x}=\boldsymbol {x}_b(t),t)}|d_b / \nu \ll 1$, the hydrodynamic force can be described as a sum of independent contributions, as expressed by the classical decomposition into steady drag, unsteady drag, added mass, fluid inertia, buoyancy and lift (Gatignol Reference Gatignol1983; Maxey & Riley Reference Maxey and Riley1983; Magnaudet & Eames Reference Magnaudet and Eames2000). In the present study, we are considering bubbles subject to a significant buoyancy force inducing an important slip relative to the carrier phase. Therefore, the Reynolds number of the bubbles is large (typically several hundred) and the abovementioned force decomposition cannot be rigorously applied. Nevertheless, as it is commonly accepted and validated by numerical simulations (Merle et al. Reference Merle, Legendre and Magnaudet2005), we will assume that this decomposition is relevant in the present case. The hydrodynamic force thus writes

(2.11)\begin{align} \mathcal{V}_b \rho_b \frac{\textrm{d} \boldsymbol{v}_b}{\textrm{d} t} &={-}\frac{1}{2} C_D \rho_f \mathcal{S}_b (\boldsymbol{v}_b - \tilde{\boldsymbol{u}}_{f,b} )|\boldsymbol{v}_b - \tilde{\boldsymbol{u}}_{f,b}|- \mathcal{V}_b \rho_f C_M \left(\frac{\textrm{d} \boldsymbol{v}_b}{\textrm{d} t} - \frac{\textrm{D} \tilde{\boldsymbol{u}}_{f,b}}{\textrm{D} t} \right) \nonumber\\ &\quad + \mathcal{V}_b \rho_f \frac{\textrm{D} \tilde{\boldsymbol{u}}_{f,b}}{\textrm{D} t} + \mathcal{V}_b (\rho_b-\rho_f) \boldsymbol{g} + \boldsymbol{F}_{\text{I},b}, \end{align}

where $\rho _b$ is the gas density, and $\mathcal {V}_b$ and $\mathcal {S}_b$ are the volume and cross-sectional area of a bubble, respectively. Here $C_D$ and $C_M$ are respectively the drag coefficient and the added-mass coefficient. The drag, the added-mass and the Tchen forces are expressed as a function of the local undisturbed fluid velocity at the position of the bubble, $\tilde {\boldsymbol {u}}_{f,b} = \tilde {\boldsymbol {u}}_{f,b}{(\boldsymbol {x} = \boldsymbol {x}_b(t),t)}$. This point is discussed in detail in the next section.

In (2.11) the Basset term has been ignored because it is expected to have a negligible effect for large-Reynolds-number objects. In non-homogenous flow conditions, when velocity gradients at large scales are present, accounting for the lift force is certainly important, as well as considering anisotropic drag and added mass associated with non-spherical bubbles that can rotate. In the present work, we focus on the fluctuations generated by a homogeneous swarm of rising bubbles for which such effects are of secondary importance. The lift force is thus not taken into account in (2.11) either.

The momentum transferred from the bubble to the fluid via the force $\boldsymbol {f}$ in (2.5) is equal to the sum of the drag and added-mass forces (first two terms on the right-hand side of (2.11)). The contributions of Tchen and Archimedes forces (third and fourth terms on the right-hand side of (2.11)) are accounted for in the pressure to be consistent with divergence free of the fluid (Climent & Magnaudet Reference Climent and Magnaudet1999). Owing to that, the momentum source writes

(2.12)\begin{align} \boldsymbol{F}_{b \rightarrow f} &={+} \frac{1}{2} C_D \rho_f \mathcal{S}_b (\boldsymbol{v}_b - \tilde{\boldsymbol{u}}_{f,b})|\boldsymbol{v}_b - \tilde{\boldsymbol{u}}_{f,b}| + \mathcal{V}_b \rho_f C_M \left(\frac{\textrm{d} \boldsymbol{v}_b}{\textrm{d} t} - \frac{\textrm{D} \tilde{\boldsymbol{u}}_{f,b}}{\textrm{D} t} \right) \nonumber\\ &= \mathcal{V}_b \rho_f \left( \frac{\textrm{D} \tilde{\boldsymbol{u}}_{f,b}}{\textrm{D} t} - \boldsymbol{g} \right) + \mathcal{V}_b \rho_b \left( \boldsymbol{g} - \frac{\textrm{d} \boldsymbol{v}_b}{\textrm{d} t} \right) + \boldsymbol{F}_{{I},b}. \end{align}

It is worth noting that, in the case of a bubbly flow ($\rho _b \ll \rho _f$) without direct interactions between the bubbles ($\boldsymbol {F}_{{I},b}=0$), the coupling term reduces to $\boldsymbol {F}_{b \rightarrow f} = \mathcal {V}_b \rho _{f} (({\textrm {D}}/{\textrm {D} t})\tilde {\boldsymbol {u}}_{f,b} - \boldsymbol {g} )$.

2.3. Model to correct the self-induced force

As seen from (2.12), when a bubble has a non-zero velocity or acceleration relative to the surrounding liquid, it disturbs the liquid velocity field. This alters in return the hydrodynamic force applied to the bubble. However, the expression of the hydrodynamic force used in the bubble motion equation (2.11) already accounts for the interaction of the bubble with its own perturbation. Therefore, the fluid velocity that intervenes in the calculation of the hydrodynamic force needs to be corrected to avoid the introduction of a spurious self-induced force exerted on the bubble.

The fluid velocity corrected from the influence of bubble $b$ is

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

Because of the nonlinearity of the problem, this immediately calls into question the definition of $\boldsymbol {u}^{*}_{f,b}$, as it appears non-trivial to isolate the effect of one bubble from the background fluctuations that include the effect of all other bubbles. Here, we propose to define the perturbation field $\boldsymbol {u}^{*}_{f,b}$ as the flow generated by an isolated bubble in a quiescent flow, which would have followed the same trajectory and exchanged as much momentum with the liquid phase as bubble $b$. Thus, $\boldsymbol {u}^{*}_{f,b}$ is formally obtained by solving the following Navier–Stokes equation:

(2.14)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}^{*}_{f,b}=0, \end{gather}$$

(2.15)$$\begin{gather}\partial_t\boldsymbol{u}^{*}_{f,b} +\boldsymbol{u}^{*}_{f,b}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}^{*}_{f,b} ={-}\boldsymbol{\nabla} p^{*} + \nu {\rm \Delta} \boldsymbol{u}^{*}_{f,b} + \frac{1}{\rho_f}\boldsymbol{F}_{b\rightarrow f}{(t)}G_{\sigma}{(\boldsymbol{x}-\boldsymbol{x}_b(t))}, \end{gather}$$

where the last term of the right-hand side is a forcing calculated by accounting for the bubble position and momentum exchange between bubble $b$ and the liquid phase.

This definition only describes the disturbance that is directly generated by bubble $b$, but ignores indirect perturbations resulting from the fact that bubble $b$ may also have affected other bubbles and from nonlinear interactions within the liquid. As a consequence, other definitions of $\boldsymbol {u}^{*}_{f,b}$ can be proposed. For example, instead of isolating the trajectory of the bubble $b$, one could define $\tilde {\boldsymbol {u}}_{f,b}$ as the flow generated by the motion of all bubbles except bubble $b$. This would give

(2.16)\begin{equation} \partial_t\tilde{\boldsymbol{u}}_{f,b} +\tilde{\boldsymbol{u}}_{f,b}\boldsymbol{\cdot}\boldsymbol{\nabla} \tilde{\boldsymbol{u}}_{f,b} ={-}\boldsymbol{\nabla} \tilde{p} + \nu {\rm \Delta} \tilde{\boldsymbol{u}}_{f,b} + \frac{1}{\rho_f} \sum_{b' \neq b}\boldsymbol{F}_{b'\rightarrow f}{(t)} G_{\sigma}{(\boldsymbol{x}-\boldsymbol{x}_{b'}(t))} . \end{equation}

Insofar as the bubble Reynolds number $Re_b$ and the volume fraction are not both vanishingly small, those two definitions of $\tilde {\boldsymbol {u}}_{f,b}$ are not equivalent. In this case, there is no exact way to define $\boldsymbol {u}^{*}_{f,b}$. Anyway, computing $\boldsymbol {u}^{*}_{f,b}$ directly from (2.14)–(2.15) (or (2.16)) implies solving one additional Navier–Stokes equation per bubble, which is not effectively feasible and we need to rely on approximate models, such as the one proposed in this work.

Previous works have mainly focused on the case of small $Re_b$. Neglecting inertia terms in (2.15), the velocity disturbance can be expressed by means of the Stokeslet Green function $S$ (Saffman Reference Saffman1973) as

(2.17)\begin{equation} \boldsymbol{u}^{*}_{f,b} = \int_{\mathcal{V}} S{(\boldsymbol{x}-\boldsymbol{x}')} \boldsymbol{F}_{b\rightarrow f}{(t)} G_{\sigma}{(\boldsymbol{x}'-\boldsymbol{x}_b(t))}\,\textrm{d}^{3}x' \approx S{(\boldsymbol{x}-\boldsymbol{x}_b(t))} \boldsymbol{F}_{b\rightarrow f}{\scriptstyle(t)}, \end{equation}

where the components of $S$ are given by

(2.18)\begin{equation} S_{ij}{(\boldsymbol{r})} = \frac{1}{8 {\rm \pi}\rho_f \nu} \left(\frac{\delta_{ij}}{|\boldsymbol{r}|} + \frac{r_i r_j}{|\boldsymbol{r}^{3}|}\right). \end{equation}

Note that the last equality in (2.17) holds for vanishingly small bubbles as its finite size is neglected and $G_{\sigma }$ is considered as the Dirac distribution. This leads to the formulation of the well-known force coupling method (see the review by Maxey Reference Maxey2017, for a detailed exposition). Subsequently, Gualtieri et al. (Reference Gualtieri, Picano, Sardina and Casciola2015) accounted for transient effects by introducing a diffusion time lag. Recently, Balachandar et al. (Reference Balachandar, Liu and Lakhote2019) extended further the Stokestlet formalism by proposing a quasisteady Galilean-invariant correction by accounting for finite Reynolds number. Alternatively, the corrected velocity field $\tilde {\boldsymbol {u}}_{f,b}(x,t)$ can be obtained from local spatial averaging of $\boldsymbol {u}_{f}$ around bubble $b$ (Capecelatro & Desjardins Reference Capecelatro and Desjardins2013; Ireland & Desjardins Reference Ireland and Desjardins2017), which is consistent when there exists a scale separation between $\boldsymbol {u}_{f}$ and $\boldsymbol {u}^{*}_{f,b}$, as discussed previously. Horwitz & Mani (Reference Horwitz and Mani2016) proposed a numerical method based on an extrapolation of the velocity field by introducing a parameter $C$ which depends on the body size $d_b$ relative to the grid spacing ${\rm \Delta} x$ and on the body positions $\boldsymbol {x}_i$: $\tilde {\boldsymbol {u}}_{f,b} = \boldsymbol {u}_f + C ( d/{\rm \Delta} x, \boldsymbol {x}_i ) {\rm \Delta} x^{2} \nabla ^{2} \boldsymbol {u}_f$. All the abovementioned approaches considered that the dispersed objects are only subject to the viscous drag force and assumed a small Reynolds number or, in the case of Balachandar et al. (Reference Balachandar, Liu and Lakhote2019), deal with finite Reynolds number by correcting the vanishing-Reynolds-number result from an adjustment of parameters obtained from DNS. In particular, none of them consider a correction of the added-mass force, which is, however, crucial during a transient regime when the bubbles accelerate from rest and still plays a role in a steady bubble swarm where the velocities of the bubbles experience fluctuations.

In this work, we propose a model to compute the velocity and the acceleration of the flow disturbance generated by a moving object, which can be used to obtain the right value of all components of the hydrodynamic force and is valid for the case of large-Reynolds-number bubbles whose size is not small compared to the grid spacing. The starting point is to express (2.15) in a reference frame (denoted by a prime), where the corrected velocity at the bubble position is null at every instant,

(2.19)\begin{equation} \left( \begin{array}{c} t \\ \boldsymbol{x} \\ \boldsymbol{u}^{*}_{f,b} \end{array} \right) \rightarrow \left( \begin{array}{l} t^{{\dagger}} = t \\ \boldsymbol{x}^{{\dagger}}=\boldsymbol{x}+\boldsymbol{\ell}_{adv}(t^{{\dagger}}) \\ \boldsymbol{u}^{{\dagger}*}_{f,b} = \boldsymbol{u}^{*}_{f,b} + \tilde{\boldsymbol{u}}_{f,b} {(\boldsymbol{x}=\boldsymbol{x}_b(t^{{\dagger}}),t^{{\dagger}})} \end{array}\right), \end{equation}

where

(2.20)\begin{equation} \boldsymbol{\ell}_{adv}(t) = \int_0^{t} \tilde{\boldsymbol{u}}_{f,b} {(\boldsymbol{x}=\boldsymbol{x}_b(t^{{\dagger}}),t^{{\dagger}})} \,\textrm{d} t^{{\dagger}} \end{equation}

represents the distance travelled by a materiel point that is advected by the corrected velocity $\tilde {\boldsymbol {u}}_{f,b}$. Similarly to Balachandar et al. (Reference Balachandar, Liu and Lakhote2019), this change of reference frame is convenient to obtain a correction that is Galilean invariant. In the new reference frame, (2.15) reads

(2.21)\begin{align} \partial_t\boldsymbol{u}^{{\dagger}*}_{f,b} +\boldsymbol{u}^{{\dagger}*}_{f,b}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}^{{\dagger}*}_{f,b} &={-}\boldsymbol{\nabla} p^{{\dagger}*} + \nu {\rm \Delta} \boldsymbol{u}^{{\dagger}*}_{f,b} + \frac{1}{\rho_f}\boldsymbol{F}_{b\rightarrow f}{(t)}G_{\sigma}{(\boldsymbol{x}^{{\dagger}}-\boldsymbol{x}^{{\dagger}}_b(t))} \nonumber\\ &\quad + d_t\tilde{\boldsymbol{u}}_{f,b} {(\boldsymbol{x}^{{\dagger}}=\boldsymbol{x}^{{\dagger}}_b(t),t)}. \end{align}

For the computation of the hydrodynamic force on a bubble, we are only interested in knowing the field $\boldsymbol {u}^{{\dagger} *}_{f,b}$ in the vicinity of this bubble. In this region, the Reynolds number is assumed to be large ($|\boldsymbol {u}^{{\dagger} *}_{f,b} | \sigma /\nu \gg 1$), so that the term $\nu {\rm \Delta} \boldsymbol {u}^{{\dagger} *}_{f,b}$ can be neglected in (2.21). It is further assumed that in the direct vicinity of the bubble, the flow is quasiparallel, since it essentially has the structure of an axisymmetric jet, which justifies also neglecting the advection and the pressure gradient terms in (2.21) and leads to

(2.22)\begin{equation} \partial_t\boldsymbol{u}^{{\dagger}*}_{f,b} = \frac{1}{\rho_f}\boldsymbol{F}_{b\rightarrow f}{(t)} G_{\sigma}{(\boldsymbol{x}^{{\dagger}}-\boldsymbol{x}^{{\dagger}}_b(t))} + d_t \tilde{\boldsymbol{u}}_{f,b} {(\boldsymbol{x}^{{\dagger}}=\boldsymbol{x}^{{\dagger}}_b(t),t)}. \end{equation}

After integration in time for a constant position $\boldsymbol {x}^{{\dagger} }$, we get

(2.23)\begin{align} \boldsymbol{u}^{{\dagger}*}_{f,b}{(\boldsymbol{x}^{{\dagger}},t)} - \boldsymbol{u}^{{\dagger}*}_{f,b}{(\boldsymbol{x}^{{\dagger}},t_0)} &= \frac{1}{\rho_f} \int_{t_0}^{t} \boldsymbol{F}_{b\rightarrow f}{(s)} G_{\sigma}{(\boldsymbol{x}^{{\dagger}}-\boldsymbol{x}^{{\dagger}}_b(s))}\,\textrm{d} s \nonumber\\ &\quad + \tilde{\boldsymbol{u}}_{f,b}{(\boldsymbol{x}^{{\dagger}}=\boldsymbol{x}^{{\dagger}}_b(t),t)} - \tilde{\boldsymbol{u}}_{f,b}{(\boldsymbol{x}^{{\dagger}}=\boldsymbol{x}^{{\dagger}}_b(t_0),t_0)}. \end{align}

The disturbance velocity field is then expressed back in the original frame of reference by applying the inverse transformation to (2.19):

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

where $\boldsymbol {\ell }_{adv}(t,s) = \int _s^{t} \tilde {\boldsymbol {u}}_{f,b} {(\boldsymbol {x}=\boldsymbol {x}_b(s^{{\dagger} }),s^{{\dagger} })}\,\textrm {d} s^{{\dagger} }$. Figure 2 is a space–time diagram that illustrates the integral formulation (2.24) and its Galilean invariance. The velocity disturbance at a given position and time is obtained from the summation over all the previous times of the momentum supplied by the bubble to the material point of the liquid that is at this specific position.

Figure 2.

Space–time diagram of the evolution of the bubble position (solid black line) and that of the material point at which momentum has been disposed at $t=t_1$ (black dashed line). (a) In the case of a steady carrier flow $\tilde {\boldsymbol {u}}_{f,b} {(\boldsymbol {x}=\boldsymbol {x}_b(t),t)}= 0$; (b) in the case of a uniform advection by the carrier flow $\tilde {\boldsymbol {u}}_{f,b} {(\boldsymbol {x}=\boldsymbol {x}_b(t),t)}=-v_0 \boldsymbol {e}_z$, with $\boldsymbol {e}_z$ the vertical direction. Galilean invariance imposes length preservation as shown in both cases by the distance $\ell (t_1,t_2)$. The shaded area represents the region of the flow influenced by the momentum disposed at $t=t_1$.

Note that the field $\boldsymbol {u}^{*}_{f,b}$ is not divergence free. Although the effect of pressure can be taken into account in a similar way as was done by Gualtieri et al. (Reference Gualtieri, Picano, Sardina and Casciola2015), the formulation proposed in (2.24) has the advantage of being simpler and allows us to analytically compute the integrals in the case of an isolated bubble, as shown in the next section. Moreover, since the determination of the correct hydrodynamic force only requires the values of the velocity correction and its derivatives at the bubble position, the divergence of $\boldsymbol {u}^{*}_{f,b}$ does not appear as a serious issue. Finally, we consider that the fluid does not reside for a long time in the region of size $\sigma$ around a bubble and thus propose to simply evaluate $\boldsymbol {\ell }_{adv}(t,s)$ as

(2.25)\begin{equation} \boldsymbol{\ell}_{adv}(t,s) \approx \tilde{\boldsymbol{u}}_{f,b} {(\boldsymbol{x}=\boldsymbol{x}_b(t),t)} \times (t-s). \end{equation}

In § 3, we propose an efficient algorithm to compute the integral in (2.24) from the particle path history. Knowing $\boldsymbol {u}^{*}_{f,b}$, it is then possible to evaluate the hydrodynamic forces involved in (2.11). The fluid velocity seen by the bubble is simply obtained by the difference between $\boldsymbol {u}_f$ and $\boldsymbol {u}^{*}_{f,b}$ as expressed by (2.13), while the corresponding fluid material derivative is given by

(2.26)\begin{equation} \left.\frac{\textrm{D} \left( \boldsymbol{u}_f - \boldsymbol{u}^{*}_{f,b} \right) }{\textrm{D} t} \right|_{\substack{\boldsymbol{x}_b(t)}} = \frac{\partial \boldsymbol{u}_f}{ \partial t} + \boldsymbol{u}_f \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}_f - \frac{\partial \boldsymbol{u}^{*}_{f,b}}{ \partial t} - \boldsymbol{u}^{*}_{f,b} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}^{*}_{f,b}.\end{equation}

It is worth mentioning that no cross-term between $\boldsymbol {u}_f$ and $\boldsymbol {u}^{*}_{f,b}$ is involved in (2.26) despite the nonlinear nature of the material derivative. The reason is that such terms must cancel to preserve the Galilean invariance of the added-mass force. The temporal and spatial derivatives of $\boldsymbol {u}^{*}_{f,b}$ can be easily computed by applying their definition to (2.24):

(2.27)\begin{align} \frac{\partial }{\partial t} \boldsymbol{u}^{*}_{f,b}{(\boldsymbol{x},t)} &= \frac{\boldsymbol{u}^{*}_{f,b} {\left(\boldsymbol{x}, t \right)} - \boldsymbol{u}^{*}_{f,b} {\left( \boldsymbol{x}, t-{\rm \Delta} t \right)} }{{\rm \Delta} t} \nonumber\\ &= \frac{1}{{\rm \Delta} t} \frac{1}{\rho_f} \int_{t-{\rm \Delta} t}^{t} \boldsymbol{F}_{b\rightarrow f}{\scriptstyle (s)} G_{\sigma}{\left(\boldsymbol{x} - \boldsymbol{x}_b(s) + \ell_{{adv}}(t-s)\right) }\,\textrm{d} s \nonumber\\ &= \frac{1}{\rho_f} \boldsymbol{F}_{b\rightarrow f}{(t)} G_{\sigma} {\left(\boldsymbol{x} - \boldsymbol{x}_b(t) \right)} \end{align}

and

(2.28)\begin{equation} \boldsymbol{\nabla} \boldsymbol{u}^{*}_{f,b}{(\boldsymbol{x},t)} = \frac{1}{\rho_f}\int_{t_0}^{t} \boldsymbol{F}_{b\rightarrow f}{(s)} \boldsymbol{\nabla} G_{\sigma} {\left(\boldsymbol{x}-\boldsymbol{x}_b(s)+\ell_{adv}(t-s)\right)} \,\textrm{d} s + \boldsymbol{\nabla}\boldsymbol{u}^{*}_{f,b} {(\boldsymbol{x},t_0)}, \end{equation}

where $\boldsymbol {\nabla } G_{\sigma }$ is the gradient of the homogenization kernel, which is analytically known for a Gaussian kernel.

3.1. A single bubble rising in an unbounded liquid at rest

We now consider the case of a single bubble rising in a liquid otherwise at rest. This simple case has the great advantage that the undisturbed fluid velocity field is trivially zero ($\tilde {\boldsymbol {u}}_b=0$) since, in the absence of the bubble, the liquid would stay at rest. First, values of the bubble drag coefficient $C_D$ and of the kernel width $\sigma$ leading to realistic simulated wakes are estimated from comparisons with experimental results. Then, we present a first validation of the method to compute the hydrodynamic force by verifying that the undisturbed liquid velocity and acceleration provided by the model (2.24) remain zero, i.e. $\boldsymbol {u}^{*}_{f,b}=\boldsymbol {u}_{f}$.

Since $\tilde {\boldsymbol {u}}_b=0$, (2.11) simplifies and the vertical bubble velocity evolves as

(3.5)\begin{equation} \frac{\textrm{d} v_{b,z}}{\textrm{d} t} = \frac{3}{4} \frac{C_D}{C_M} \frac{v_{b,z}^{2} }{d_b} - \frac{g}{C_M}, \end{equation}

where $\boldsymbol {g} = -g \boldsymbol {e}_z$. In general, the added-mass coefficient of a bubble depends on its shape while the drag coefficient depends both the shape and the Reynolds number (Moore Reference Moore1965; Ishii & Zuber Reference Ishii and Zuber1979; Maxworthy et al. Reference Maxworthy, Gnann, Kürten and Durst1996). In this study, we choose the value of the drag coefficient to be consistent with the average rising speed of the experiments of Ellingsen & Risso (Reference Ellingsen and Risso2001) and Riboux et al. (Reference Riboux, Risso and Legendre2010). In their experiments, an air bubble with a diameter $d_b=2.5$ mm rises in water at an average speed of $v_0=0.305$ m s$^{-1}$, which corresponds to a bubble Reynolds number $Re_b = v_0 d_b /\nu \approx 760$ and an average drag coefficient

(3.6)\begin{equation} C_D = \frac{4}{3} g \frac{d_b}{v_{0}^{2}} \approx 0.35. \end{equation}

When many bubbles are present, their interactions make their velocities fluctuate and cause their drag coefficients to increase. There have been several attempts to determine the drag coefficient in a bubble swarm. In particular, Akiki, Jackson & Balachandar (Reference Akiki, Jackson and Balachandar2017) recently proposed a numerical procedure to account for the dependence of $C_D$ on the flow configuration and particle distribution. In the present large-scale simulations, we focus on the fluctuations resulting from wake interactions, which are expected to weakly depend on the details of the flow at the bubble scale. For this reason, we will not consider the variations of the drag coefficient associated with the presence of the other bubbles and we will adopt $C_D=0.35$ for our subsequent simulations of a rising bubbles swarm. Similarly, we will fix the added-mass coefficient $C_M$ to the value 1/2 corresponding to a sphere, thereby neglecting the variations of shapes and orientations of the bubbles, which would require the introduction of an anisotropic added-mass tensor.

In the experiment of Ellingsen & Risso (Reference Ellingsen and Risso2001), the fluid velocity in the wake of an isolated bubble has been measured. We can thus adjust the value of the kernel width $\sigma$ for a given mesh resolution so as to reproduce the experimental bubble wake. According to (2.12), the momentum forcing transferred to the liquid by a single rising bubble can be expressed as

(3.7)\begin{equation} \boldsymbol{F}_{b \rightarrow f} ={-}\mathcal{V}_b \rho_{f} \boldsymbol{g} = F_0 \boldsymbol{e_z}. \end{equation}

Also, the bubble trajectory is known analytically from (3.5). Therefore, we can integrate the Navier–Stokes equations (2.4)–(2.5) to compute the fluid velocity induced by the rising of a single bubble, without explicitly considering the coupling with the dispersed phase. For such a simulation with a prescribed bubble motion, two numerical parameters need to be set: the grid spacing ${\rm \Delta} x$ and the width $\sigma$ of the Gaussian kernel. (Note that we only consider here a uniform mesh with the same spacing in the three directions.)

Several simulations have been performed to analyse the influence of these parameters upon the induced-velocity field. In any case, the computational domain is a parallelepiped of dimensions $L_x=L_y$ and $L_z$, with periodic boundary conditions in all three directions. The domain height $L_z$ is set such that a bubble reaches its terminal velocity and travels about 125 times its diameter before meeting its own wake because of the domain periodicity. Since the wake vanishes very slowly, it would require a longer box for a bubble not to interact at all with its own wake and a truly stationary state to be observed, which indicates that interactions between bubbles are expected to occur even at very small volume fractions. We consider 5 different mesh resolutions ranging from ${\rm \Delta} x / d_b = 0.55$ to 0.034. The details of the computational grids are given in table 1. For each of these grids, we have tested 7 values of $\sigma$ ranging from $0.071 d_b$ to $0.57 d_b$, as summarized in the parameter map of figure 5. For all simulations, the bubble diameter, the kinematic viscosity and the gravity are kept constant as respectively $d_b=2.5$ mm, $\nu = 10^{-6}$ m$^{2}$ s and $g=9.81$ m s$^{-2}$. The resulting wake is illustrated in figure 6. To determine the appropriate value of $\sigma /d_b$ for a given value of ${\rm \Delta} x/d_b$, we compare the vertical fluid velocity along a vertical axis passing through the centre of the bubble with the measurement of Ellingsen & Risso (Reference Ellingsen and Risso2001). This spatial evolution of the fluid velocity is shown in figure 7 for ${\rm \Delta} x/d_b=0.068$ (mesh G3). It is clearly seen that the value of $\sigma /d_b$ strongly influences the results of the numerical simulations. Owing to the self-similarity property of the wakes, we expect that for distances far enough away from the bubble, the wake depends only on the bubble Reynolds number and is no longer influenced by the details of the momentum source. In the near-bubble region, a smaller kernel width will result in a better prediction by the numerical model. It is indeed observed that for large $\sigma$, the fluid velocity is underestimated compared with the experiment, because, in that case, the momentum source is spread over too large a region, while with small values of $\sigma$, the velocity peak in the vicinity of the bubble tends to be well reproduced. Nevertheless, we observe the apparition of instabilities for $\sigma /d_b<0.2$, which are visible in figure 7 and probably caused by too high a momentum source density, which leads to a fluid velocity at the bubble position that is larger than the bubble velocity. It is worth mentioning that the features of this instability, when it develops, do not depend on the mesh-size spacing or on the time step. Similar trends regarding the influence of $\sigma$ on the wake, not reported here, are also observed for the other mesh resolutions. In particular, we observe that when $\sigma$ is smaller than the grid spacing, the entire momentum source is contained in a single computational cell (as in the ‘particle-in-cell’ approach) and there is thus no more effect of further reducing $\sigma$. To summarize, we observe the best agreement between the numerical simulations and the experiments for $\sigma /d_b \approx 0.25$ with grids G2, G3 and G4, while for grid G1 the wake is better reproduced with $\sigma /d_b \approx 0.21$. Figure 7(b) compares these settings. No difference is observed between the results obtained with meshes G2 and G3, which indicates that grid convergence can be reached.

Figure 5.

Map of the parameters ${\rm \Delta} x/d_b$ and $\sigma /d_b$ considered for the simulations of the single rising bubble. The circles denote cases in which the wake is well reproduced and the model for the disturbance gives accurate estimations, whereas triangles indicate cases where the wake is poorly described due to a lack of resolution, and crosses correspond to unstable simulations.

Figure 6.

Illustration of the liquid vertical velocity field $u_{f,b}$ in a $x$–$z$ plane passing through the bubble centre for the simulation with grid G3 (${\rm \Delta} x/d_b = 0.068$) and $\sigma /d_b=0.24$. (a) View of the full numerical domain, (b) zoom-in close to the bubble. The bubble is represented by the red circle.

Figure 7.

Vertical velocity of the fluid in the wake of the bubble normalized by the terminal velocity of the bubble: (a) for different values of $\sigma /d_b$ and mesh G3. Red line, experimental result of Ellingsen & Risso (Reference Ellingsen and Risso2001); black symbols, numerical results for various values of $\sigma /d_b$; (b) for grids G2 and G3, and $\sigma /d_b= 0.24$ and 0.28.

Table 1.

Characteristics of the different computational grids for the test of a single rising bubble. In all simulations the bubble diameter, the kinematic viscosity and the gravity are kept constant as respectively $d_b=2.5$ mm, $\nu = 10^{-6}$ m$^{2}$ s$^{-1}$ and $g=9.81$ m s$^{-2}$.

We now focus on the validation of our model for the disturbed fluid velocity. The simulations just mentioned above, for which the bubble trajectory and momentum forcing were imposed from (3.5) and (3.7), give the reference liquid velocity disturbance, since the undisturbed velocity was identically zero. This reference is used to determine the constants of our model in (3.2), (3.3) and (3.4) as follows. Once the velocity of the bubble is stationary, (3.2)–(3.4) can be evaluated analytically, which leads to relate $c_1$, $c_2$ and $c_3$ to $c_0$ by

(3.8)$$\begin{gather} c_1 = u^{*}_{f,b,z} \frac{ 4{\rm \pi} \sigma^{2} c_0^{2} v_0 }{F_0/\rho_f}, \end{gather}$$

(3.9)$$\begin{gather}c_2 ={-} \partial_z u^{*}_{f,b,z} \frac{(2{\rm \pi} \sigma^{2} c_0^{2})^{3/2} v_0 } {F_0/\rho_f}, \end{gather}$$

(3.10)$$\begin{gather}c_3 = \partial_t u^{*}_{f,b,z} \frac{(2{\rm \pi} \sigma^{2} c_0^{2} )^{3/2}} {F_0/\rho_f}, \end{gather}$$

where the values of the liquid velocity disturbance and its derivatives at the bubble centre, $u^{*}_{fb,z}$, $\partial _z u^{*}_{f,b,z}$ and $\partial _t u^{*}_{f,b,z}$, are taken from the uncoupled simulation, in order to have $\boldsymbol {u}_{f} = \boldsymbol {u}^{*}_{f,b}$.

Then, the performance of our model is evaluated by comparing its prediction to the reference simulation. In figure 8(a), we examine the fluid vertical velocity at the bubble position for grid G2 and $\sigma /d_b=0.25$. In this case the value of $c_0$ has been adjusted to 1.62 by least-square fitting between the reference simulation and (3.2), which leads to the other parameters being $c_1=2.88$, $c_2=1.17$ and $c_3=1.17$ by using (3.8)–(3.10). We present the temporal evolution of the actual fluid velocity at the bubble position, as well as the velocity disturbance computed from the bubble path history according to (3.2). It is seen that both curves overlap very well, which indicates that the undisturbed velocity is indeed zero. Moreover, the comparison with the reference simulation is also very good, which provides a first validation of the model used to calculate the disturbance. In figure 8(b), we show a similar comparison for the vertical component of the total fluid acceleration. It is seen as well that the fluid acceleration at the bubble position is nearly equal to the computed acceleration disturbance and to the reference fluid acceleration at the bubble position. The undisturbed fluid acceleration is thus zero and we can conclude that the perturbed acceleration is accurately computed. It should be noted that, for $\sigma /d_b \leq 0.21$, the velocity and acceleration disturbances show some deviation compared with the reference values.

Figure 8.

Temporal evolution of the fluid vertical velocity at the bubble position normalized by the bubble terminal velocity (a) and the fluid vertical acceleration at the bubble position normalized by gravity (b). Red triangles, results of the reference simulation with prescribed momentum source; blue lines, results from the full simulation including the correction of the hydrodynamic force; black dashed lines, velocity or acceleration disturbance given by the model implemented in the full simulation; green lines, results of the full simulation without accounting for the correction of the hydrodynamic force ($\boldsymbol {u}^{*}_{f,b} = 0$). These tests have been carried out with grid G2 (${\rm \Delta} x / d_b = 0.14$) and $\sigma /d_b=0.28$. The inset of panel (b) is a zoom-out that shows uncorrected results.

Since the unperturbed fluid velocity and acceleration are accurately predicted, the hydrodynamic force on the bubble can also be accurately computed. This is confirmed by figure 9, which compares the reference analytical evolution of the bubble velocity obtained from (3.5) with the evolution obtained from the simulation of the full set of coupled equations. The analytical and simulated bubble velocities are very close to each other, which validates the proposed method to cancel the spurious self-induced force. To emphasize the importance of this correction in the case of a single bubble rising at large Reynolds number, figures 8 and 9 also show the results of the coupled simulation obtained without accounting for the correction of the liquid velocity at the bubble location, i.e. by setting $\boldsymbol {u}^{*}_{f,b}=0$. Without correction, the computed values are indeed largely off from the reference values and can even further lead to the divergence of the numerical simulation.

Figure 9.

Vertical velocity of the bubble normalized by its terminal velocity. Red triangles, reference case ($\tilde {\boldsymbol {u}}_{f,b}=0$); green line, case without correction ($\tilde {\boldsymbol {u}}_{f,b}=\boldsymbol {u}_f$); blue line, results accounting for the correction ($\tilde {\boldsymbol {u}}_{f,b} = \boldsymbol {u}_f - \boldsymbol {u}^{*}_{f,b}$). The inset is a zoom-out that shows the uncorrected results.

Now, it is useful to estimate when the correction of the spurious self-induced force is required. When $| \tilde {\boldsymbol {u}}_{f,b} - \boldsymbol {v}_b| \gg |\boldsymbol {u}^{*}_{f,b}|$, the velocity disturbance has a negligible effect on the drag force. Additionally, when $| D_t \tilde {\boldsymbol {u}}_{f,b} - d_t \boldsymbol {v}_b| \gg | D_t \boldsymbol {u}^{*}_{f,b}|$, the effect on the forces resulting from the liquid inertia is negligible. These inequalities can easily be tested a posteriori, from the knowledge of the bubble path and with the help of the proposed model for $\boldsymbol {u}^{*}_{f,b}$. They can be also be estimated a priori by assuming that the momentum exchange in (2.12) is only due to the drag force (which is exact in the case for a single rising bubble at stationary state and approximate for a bubble in a bubble swarm) and by considering that the integrand in (2.24) is non-negligible only over time of the order of $\sigma / | \tilde {\boldsymbol {u}}_{f,b} - \boldsymbol {v}_b|$ (see figure 2). This leads to the following scaling for the ratio of the velocity disturbance to the relative velocity between the bubble and the liquid:

(3.11)\begin{equation} \frac{|\boldsymbol{u}^{*}_{f,b}| }{ | \tilde{\boldsymbol{u}}_{f,b} - \boldsymbol{v}_b|} \sim C_D \frac {d_b^{2}}{\sigma^{2}}. \end{equation}

From (3.11), we conclude that it is essential to correct the fluid velocity for large-Reynolds-number bubbles, since, in this case, $C_D = O(1)$ and $\sigma /d_b = O(1)$. Note that for objects moving at small Reynolds number, an estimate obtained from the averaging of the Stokeslet (2.17) gives $|\boldsymbol {u}^{*}_{f,b}| / | \tilde {\boldsymbol {u}}_{f,b} - \boldsymbol {v}_b| = O(d_b/{\rm \Delta} x)$ (Saffman Reference Saffman1973; Boivin et al. Reference Boivin, Simonin and Squires1998), which means that the correction can be disregarded for sufficiently small objects and coarse description of the velocity field.

It is worth noting that even if coupling between the two phases through the hydrodynamic force acting on the bubbles and the transfer of momentum to the liquid by diffused sources ensure the conservation of momentum for the mixture, it does not conserve its energy (Xu & Subramaniam Reference Xu and Subramaniam2007; Subramaniam et al. Reference Subramaniam, Mehrabadi, Horwitz and Mani2014). Indeed, the power $P_b$ of the hydrodynamic force working at the bubble velocity is larger than the power $P_f$ of the diffused force working at the fluid velocity:

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

However, it is physically acceptable that a certain amount of energy is dissipated within the interphase coupling, since we consider a coarse-grained description of the continuous phase, in which the strong velocity gradients in the vicinity of the bubbles are not resolved. In other words, the dissipation of kinetic energy into heat that occurs in the region surrounding the bubble at scales smaller than $\sigma$ cannot be computed from $\boldsymbol {u}_f$.

Nevertheless, the energy dissipated in the coupling can be estimated from the model for the perturbed velocity $\boldsymbol {u}^{*}_{f,b}$. Let us do this for the case of the isolated bubble. From (2.6), (2.24) and (3.7), and considering $\boldsymbol {u}_f=\boldsymbol {u}^{*}_{f,b}$, we obtain

(3.13)\begin{equation} P_f = \frac{F_0^{2}}{\rho_f} \int G_{\sigma}( \boldsymbol{x}-\boldsymbol{x}_b(t)) \int_{-\infty}^{t} G_{\sigma}( \boldsymbol{x}-\boldsymbol{x}_b(s))\,\textrm{d} s {\textrm{d}{\kern 0.2em}x}^{3}. \end{equation}

By introducing the non-dimensional variables $\boldsymbol {\chi } = (\boldsymbol {x}-\boldsymbol {x}_b(t)) / \sigma$ and $\tau = (t-s) v_0/\sigma$, we get

(3.14)\begin{equation} P_f = \frac{F_0^{2}}{\rho_f \sigma^{2} v_0} \underbrace{\int \int_{0}^{\infty} G_{1}( \boldsymbol{\chi} ) G_{1}( \boldsymbol{\chi} + \tau\boldsymbol{e}_z )\,\textrm{d}\tau\,\textrm{d}\chi^{3} }_{k}, \end{equation}

with $\boldsymbol {e}_z$ the vertical direction. Owing to $G$ being a Gaussian kernel, the summation can be evaluated analytically and gives $k=1/8{\rm \pi}$. It should be remarked that if we account for the numerical parameters in the computation of the velocity perturbation (that is to say, taking (3.2) instead of (2.24)), we get $k=c_1/c_0^{2} 8{\rm \pi}$. Since the energy given by the bubble is $P_b = F_0 v_0$, we have

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

For $C_D=0.35$ and $\sigma /d_b = 0.28$, this estimate leads to ${P_f}/{P_b} \approx 0.07$, which means that most of the mechanical energy that is exchanged between the bubble and the liquid is dissipated in the vicinity of the bubble.

In Appendix A, we present some tests showing that the proposed model respects Galilean invariance.

3.2. Repulsive force between bubbles

Since the sources of momentum that model the bubbles are spanned over a few mesh-grid steps, it is possible that two bubbles overlap when simulating multiple bubble trajectories. To prevent such a non-physical behaviour, we introduce an interaction force between bubbles that prevents them from overlapping. In the bubble momentum balance equation (2.11), the term $\boldsymbol {F}_{{I},b}$ is expressed as

(3.16)\begin{equation} \boldsymbol{F}_{{I},b} = \sum_{b'\neq b} \boldsymbol{F}_{r_{b', b}}, \end{equation}

with

(3.17)\begin{equation} \boldsymbol{F}_{r_{b', b}} ={-} C_I \frac{\boldsymbol{x}_{b'}-\boldsymbol{x}_b }{r_{b', b}} \exp\left( - r_{b', b}^{2} / 2 r_{I}^{2} \right), \end{equation}

where $r_{b', b}=| \boldsymbol {x}_{b'}-\boldsymbol {x}_b |$ is the distance between the centres of bubbles $b$ and $b'$. The range of the force is limited by setting $\boldsymbol {F}_{r_{b' \rightarrow b}} = 0$ for $r_{b', b} > 3 r_{I}$. The parameter $r_{I}$ is set as $d_b/2$ whereas $C_I$, which has a unit of force, is adjusted so that two rising bubbles reach an equilibrium distance $r_{b', b}$ equal to $d_b$. For grid G3 with ${\rm \Delta} x/d_b=0.068$ and $\sigma /d_b=0.28$, this corresponds to $C_I/(\rho _f {v_0}^{2} d^{2})=0.34$.

The effect of this force is illustrated in figure 10, which shows the temporal evolution of the positions of two vertically aligned rising bubbles. When the repulsive force is implemented, the distance between the bubble centres reaches a stationary value corresponding to $r_{b', b}= d_b$ and the terminal velocity of the pair is approximately 30 % larger than that of a single bubble. The actual interaction between two rising bubbles is more complicated (van Wijngaarden Reference van Wijngaarden1993; Harper Reference Harper1997; Hallez & Legendre Reference Hallez and Legendre2011). However, our motivation to include this repulsive force is just to ensure numerical stability, since the bubble velocity can diverge when two bubbles are allowed to overlap.

Figure 10.

(a) Temporal evolution of the distance between the two bubbles rising in line. Continuous black line, with repulsive force and dashed line with triangles, without repulsive force. (b) Temporal evolution of the vertical velocity of the two bubbles. The first bubble corresponds to the red curve, the second bubble corresponds to the black curve. Comparison of the cases with or without repulsive force shown respectively by the continuous lines with symbols and dashed lines.