2.1 General formulation

In this section, we consider that the liquid fluctuations $\boldsymbol{u}$ generated in a bubble swarm are described by the sum of the velocity disturbances $\boldsymbol{v}$ generated by each individual bubble. A bubble located at location $x_{i}$ and rising at velocity $V$ induces a velocity field $V\tilde{\boldsymbol{v}}(\boldsymbol{x}_{\boldsymbol{i}})$ at the origin.

We assume that the swarm is dilute enough to neglect any interactions between the bubbles. Consequently, (i) the flow disturbance $\tilde{v}$ is a given function of only the bubble location, and (ii) the bubble positions $\boldsymbol{x}_{\boldsymbol{i}}$ are random variables that are statistically independent from each other and uniformly distributed in space. For $N$ bubbles in a volume $\unicode[STIX]{x1D717}$ , the total velocity fluctuation

(2.1) $$\begin{eqnarray}\boldsymbol{u}=V\mathop{\sum }_{i=1}^{N}\tilde{\boldsymbol{v}}(\boldsymbol{x}_{\boldsymbol{ i}})\end{eqnarray}$$

is a random variable. In the limit of $\unicode[STIX]{x1D717}$ tending towards infinity, the statistical distribution of $\boldsymbol{u}$ is fully characterized by the knowledge of its cumulants, which can be obtained from the summation over space of the flow disturbance generated by a single bubble (Rice Reference Rice1944),

(2.2) $$\begin{eqnarray}k_{n}={\tilde{N}}V^{n}\int _{x\in \unicode[STIX]{x1D717}}\tilde{\boldsymbol{v}}^{n}(\boldsymbol{x})\,\text{d}\unicode[STIX]{x1D717},\end{eqnarray}$$

where ${\tilde{N}}=N/\unicode[STIX]{x1D717}$ is the bubble number density. This result is valid provided the velocity disturbance $\tilde{\boldsymbol{v}}$ generated by each bubble decays fast enough with the distance to the bubble such that the integrals in (2.2) converge. Physically, that means that the kinetic energy of the steady flow induced by a single moving bubble has to be finite.

The moments $m_{n}$ of the statistical distribution are related to the cumulants by the classic recurrence relation,

(2.3) $$\begin{eqnarray}k_{n}=m_{n}-\mathop{\sum }_{i=1}^{n-1}C_{i-1}^{n-1}k_{i}m_{n-i},\end{eqnarray}$$

where $C_{i}^{\,j}$ denotes the binomial coefficients. This relation is particularly simple for the moments of the first three orders: the mean is equal to the first cumulant ( $m_{1}=k_{1}$ ), the variance to the second cumulant ( $\unicode[STIX]{x1D70E}^{2}=m_{2}-m_{1}^{2}=k_{2}$ ), and the third-order centred moment to the third cumulant ( $\unicode[STIX]{x1D707}_{3}=m_{3}-3m_{1}m_{2}+2m_{1}^{3}=k_{3}$ ). Even if the cumulants (or the moments) are sufficient to fully define the statistical distribution, we are presently interested by obtaining the p.d.f. This is particularly simple in the limiting case where ${\tilde{N}}$ tends towards infinity. Equation (2.2) shows that all $k_{n}$ are proportional to ${\tilde{N}}$ and hence $k_{n}/k_{2}^{n/2}$ tends towards zero for all $n$ greater than 2. Consequently, only the first two cumulants of the limit distribution are non-zero, which means that it is Gaussian. The fact that the statistics tend towards a Gaussian distribution as the number density is increased can be seen as a direct consequence of the central limit theorem.

In reality, the finite size of the bubbles limits the value of ${\tilde{N}}$ because bubbles cannot penetrate into each other. Gaussian distributions are therefore not expected to be observed. To be approximated by a Gaussian distribution, the velocity at a given point must result from the addition of a large enough number of random bubble disturbances. A good parameter to quantify the convergence towards a Gaussian distribution is hence the fraction of the total volume that is disturbed by the bubbles, which is proportional to both the number density ${\tilde{N}}$ and the volume of the region perturbed by a single bubble. At low Reynolds number, the slow decay of the velocity disturbance allows long-range interactions. That is probably the reason why a Gaussian distribution has been predicted for particle volume fractions of approximately 35 %. Except that the maximum of the p.d.f.s becomes slightly more rounded, the convergence towards a Gaussian distribution is not observed in figure 1 for a swarm of large-Reynolds-number bubbles at volume fractions up to 8 %.

We therefore need to determine the shape of the p.d.f.s for finite values of ${\tilde{N}}$ . Unfortunately, the p.d.f. cannot be easily obtained analytically even when the expressions of all the cumulants are known. In particular, it is worth recalling that a truncated sequence of cumulants does not represent a valid probability law and therefore cannot be used as a relevant approximation (except for either the uniform or the Gaussian distribution, for which only the first two cumulants are non-zero). In the following, we will determine the p.d.f.s numerically by adding disturbances that are randomly distributed over space.

2.2 Contribution of the average wakes

Let us consider an isolated axisymmetric body rising at constant velocity. In the laminar regime, the flow profile in the wake is Gaussian in the radial direction $r$ and decreases as the reciprocal of the distance $z$ to the body: $\tilde{\boldsymbol{v}_{\boldsymbol{z}}}(\boldsymbol{r},\boldsymbol{z})\propto (L/z)\exp (-r^{2}/w^{2})$ , where $L$ and $w$ are two length scales that depend on the body size, drag coefficient and Reynolds number (Batchelor Reference Batchelor1967). For a turbulent wake, the decrease is even slower (Tennekes & Lumley Reference Tennekes and Lumley1972): $\tilde{v_{z}}(r,z)\propto (L/z)^{2/3}\exp (-r^{2}/w^{2})$ . As remarked by Parthasarathy & Faeth (Reference Parthasarathy and Faeth1990), if we sum the wakes of randomly distributed bodies, integral (2.2) diverges for $n=2$ , leading to an infinite velocity variance in both the laminar and turbulent cases. However, it has been observed (Risso & Ellingsen Reference Risso and Ellingsen2002) that the wake of a high-Reynolds-number bubble belonging to a swarm of bubbles vanishes at a distance of a few bubble diameters provided the volume fraction is at least a few tenths of one per cent. The wake of an isolated bubble is therefore not relevant for bubbly flows. Indeed, at large Reynolds number, the wake behind a body immersed in a swarm of similar bodies turns out to decrease exponentially with the distance (Risso et al. Reference Risso, Roig, Amoura, Riboux and Billet2008): $\tilde{v_{z}}(z)\propto \exp (-z/L)$ .

The mechanism responsible for this strong wake attenuation is not an enhanced diffusion induced by the turbulent fluctuations that develop within the swarm. First, it has been observed that the wakes do not spread significantly before they vanish (Amoura Reference Amoura2008). Second, the action of an external turbulence leads to a power-law decay (Wu & Faeth Reference Wu and Faeth1994, Reference Wu and Faeth1995; Bagchi & Balachandar Reference Bagchi and Balachandar2004; Legendre, Merle & Magnaudet Reference Legendre, Merle and Magnaudet2006; Risso et al. Reference Risso, Roig, Amoura, Riboux and Billet2008), which is not able to cause the strong exponential attenuation observed in bubbly flows. The wake cancellation is caused by interactions between the wakes and can probably be explained by the mechanisms proposed by Hunt & Eames (Reference Hunt and Eames2002).

For a homogeneous swarm of rising bubbles, $L$ has been found to scale as the ratio between the bubble diameter and its drag coefficient, and to be independent of the gas volume fraction in the range ${\sim}$ 0.5–10 % (Risso et al. Reference Risso, Roig, Amoura, Riboux and Billet2008). It is therefore reasonable to model the velocity disturbance generated by each individual bubble located at the point of coordinates $(x,y,z)$ by

(2.4) $$\begin{eqnarray}\left.\begin{array}{@{}l@{}}\tilde{v}_{x}(x,y,z)=0,\\ \tilde{v}_{y}(x,y,z)=0,\\ \tilde{v}_{z}(x,y,z)=0\quad \text{for }z<0,\\ \tilde{v}_{z}(x,y,z)=\exp (-z/L)\exp (-(x^{2}+y^{2})/w^{2})\quad \text{for }z\geqslant 0,\end{array}\right\}\end{eqnarray}$$

where the characteristic length $L$ and width $w$ of the wake are independent of $\unicode[STIX]{x1D6FC}$ . It is important to stress that expression (2.4) represents the wake of a bubble within the bubble swarm and accounts for the attenuation due to interactions with the other wakes.

Substituting this expression in (2.3) yields the expression of the cumulants of the fluctuations generated by summing randomly distributed wakes,

(2.5) $$\begin{eqnarray}k_{n}={\tilde{N}}\frac{\unicode[STIX]{x03C0}Lw^{2}}{n^{2}}V^{n}.\end{eqnarray}$$

The corresponding p.d.f.s can be estimated analytically in the limits of ${\tilde{N}}$ tending towards either zero or infinity. In the limit of vanishing bubble number density, the probability that two wakes intersect is very low. Consequently, we can assume that the liquid velocity at a given point is due only to the wake of a single bubble. Under this assumption, the analytical expression of the probability density $f_{z}$ of the normalized vertical velocity fluctuations, $\tilde{u_{z}}=u_{z}/V$ , is found to be

(2.6) $$\begin{eqnarray}f_{z}(\tilde{u_{z}})=\unicode[STIX]{x1D705}\,\frac{1}{\tilde{u_{z}}+M}\log \left(\frac{1}{\tilde{u_{z}}+M}\right)\quad \text{for }\unicode[STIX]{x1D716}-M\leqslant \tilde{u_{z}}\leqslant 1-M,\end{eqnarray}$$

where $\unicode[STIX]{x1D705}=6\unicode[STIX]{x1D6FC}\tilde{w}^{2}\tilde{L}$ , $\unicode[STIX]{x1D716}=\exp (-\sqrt{2/\unicode[STIX]{x1D705}})$ , $M=\unicode[STIX]{x1D705}(1-\unicode[STIX]{x1D716}-\unicode[STIX]{x1D716}\sqrt{2/\unicode[STIX]{x1D705}})$ , and $\tilde{L}=L/d$ and $\tilde{w}=w/d$ are the normalized length and width of the wake (see appendix A for a detailed derivation).

In the limit of large bubble number density, the p.d.f. tends towards a Gaussian function,

(2.7) $$\begin{eqnarray}f_{z}(\tilde{u} _{z})=\frac{1}{\tilde{\unicode[STIX]{x1D70E}_{w}}\sqrt{2\unicode[STIX]{x03C0}}}\exp \left(-\frac{\tilde{u} _{z}^{2}}{2\tilde{\unicode[STIX]{x1D70E}_{w}}^{2}}\right),\end{eqnarray}$$

with a variance $\unicode[STIX]{x1D70E}_{w}^{2}$ obtained by taking $n=2$ and $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x03C0}d^{3}{\tilde{N}}/6$ in (2.5),

(2.8) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D70E}}_{w}^{2}=\left(\frac{\unicode[STIX]{x1D70E}_{w}}{V}\right)^{2}=\frac{3}{2}\unicode[STIX]{x1D6FC}\,\tilde{w}^{2}\tilde{L}.\end{eqnarray}$$

For intermediate gas volume fractions, the p.d.f. has been obtained numerically by the following method. A cubic domain $\unicode[STIX]{x1D717}_{T}$ of size $L=40d$ is considered. The coordinates $(x_{i},y_{i},z_{i})$ of $N=6\unicode[STIX]{x1D6FC}L^{3}/\unicode[STIX]{x03C0}d^{3}$ bubbles are chosen randomly and independently of each other such that their distribution be uniform over $\unicode[STIX]{x1D717}_{T}$ . The velocity induced by each bubble, given by (2.4), is computed at the middle of the domain. (Note that, because the wakes only develop behind the bubbles, only the bubbles located in the upper half of the domain contribute to the fluctuation.) A velocity sample is obtained by summing the contributions of the $N$ bubbles. The operation is repeated a large number, $N_{s}$ , of times in order to build a statistically representative ensemble of samples. The average velocity is then subtracted from each sample to get the fluctuating velocity. Several values of $N_{s}$ have been tested. For all cases considered in this work (figure 2 and following), increasing $N_{s}$ beyond $10^{5}$ no longer changes the global shape of the p.d.f.s but only reduces the amplitude of the high-frequency scattering around this global shape, which is a residue of an imperfect statistical convergence. Finally, we have chosen $N_{s}=5\times 10^{5}$ , which was found to be a good compromise between computation time and the smoothness of the numerical curves.

Figure 2.

The p.d.f.s of liquid velocity fluctuations induced by average bubble wakes described by (2.4). (a,b) Results at various gas volume fractions for a bubble of 2.5 mm (case C in table 1). The two light (red) curves in (a) correspond to the analytical solution for low bubble number density (2.6) for $\unicode[STIX]{x1D6FC}=0.34\,\%$ and $\unicode[STIX]{x1D6FC}=0.91\,\%$ . The two light (red) curves in (b) correspond to the Gaussian solution valid at large bubble number density (2.7) and (2.8) for $\unicode[STIX]{x1D6FC}=40\,\%$ and $\unicode[STIX]{x1D6FC}=80\,\%$ . (c) Results for various wake parameters at a single gas volume fraction. (d) Effect of bubble path inclination.

Typical p.d.f.s of statistical samples generated by this method are shown in figure 2. In these examples, the model parameters, $\tilde{L}$ and $\tilde{w}$ , have been chosen to correspond to a bubble of $d=2.5$  mm rising in water (case C in table 1). The discussion of these values will be presented in § 4 when comparing the model with experiments. Figure 2(a,b) shows the p.d.f.s for various gas volume fractions from 0.34 % to 80 % for a fixed pair of parameters $\tilde{L}=3$ and $\tilde{w}=0.6$ . Considering values of $\unicode[STIX]{x1D6FC}$ lower than a few tenths of one per cent would not be relevant with the assumption that the wake dimensions are independent of $\unicode[STIX]{x1D6FC}$ . Considering values of $\unicode[STIX]{x1D6FC}$ much greater than 10 % is probably not physically relevant because strong hydrodynamic interactions between bubbles should take place. We have however plotted p.d.f.s up to unreasonably large values of $\unicode[STIX]{x1D6FC}$ to examine the convergence towards a Gaussian function.

As expected, the p.d.f.s are well described by the low-gas-volume-fraction approximation as long as the probability of occurrence of fluctuations of order one or greater remains negligible. In the present case, the agreement of the computed p.d.f.s with (2.6) is excellent at $\unicode[STIX]{x1D6FC}=0.34$  % and good until $\tilde{u} _{z}=0.5$ at $\unicode[STIX]{x1D6FC}=0.91$  %. Conversely, the convergence towards the other asymptotic solution (2.7) and (2.8) for increasing $\unicode[STIX]{x1D6FC}$ is slow. The Gaussian function only starts to be a relevant approximation for a gas volume fraction as large as 80 %, and has therefore no chance of being observed in an experiment. At intermediate gas volume fractions ( $1\,\%\leqslant \unicode[STIX]{x1D6FC}\leqslant 40\,\%$ ), the p.d.f.s are strongly asymmetric. When $\unicode[STIX]{x1D6FC}$ is increased, negative fluctuations, which are initially totally absent, developed very slowly as the p.d.f. maximum gradually becomes more and more rounded. Conversely, positive fluctuations show a long exponential tail, which is still present at $\unicode[STIX]{x1D6FC}=40$  %.

We discuss now how the p.d.f. depends on the wake width and length. Figure 2(c) shows some results obtained for various values of $\tilde{w}$ and $\tilde{L}$ . We observe that increasing $\tilde{w}$ at constant $\tilde{L}$ has an effect similar to increasing $\unicode[STIX]{x1D6FC}$ . Increasing $\tilde{L}$ at constant $\tilde{w}$ (not shown here) is found to have the same effect. Moreover, it turns out that a similar p.d.f. is obtained for various values of $\tilde{L}$ and $\tilde{w}$ provided the wake volume, $\tilde{L}\tilde{w}^{2}$ , is kept constant. In fact, the p.d.f. depends on the fraction of the total volume that is occupied by the wake, which is proportional to $\unicode[STIX]{x1D6FC}\tilde{L}\tilde{w}^{2}$ .

Equation (2.4) only accounts for the longitudinal velocity. It is reasonable to neglect the transverse velocity because a wake is a quasi-parallel flow. However, because of wake-induced oscillations (Ern et al. Reference Ern, Risso, Fabre and Magnaudet2012), large-Reynolds-number bubbles do not rise on rectilinear paths. In order to assess a possible effect of bubble path inclinations on the p.d.f.s, additional computations have been carried out by considering that the velocity disturbance generated by each bubble is due to a wake that is inclined relative to the vertical direction. The inclination angles have been randomly chosen and uniformly distributed in the interval $[0,\unicode[STIX]{x1D703}_{max}]$ . The maximum inclination angle measured by Riboux et al. (Reference Riboux, Risso and Legendre2010) was $28^{\circ }$ . Figure 2(d) shows an example of p.d.f.s of the vertical and horizontal velocity fluctuations obtained for $\unicode[STIX]{x1D703}_{max}=30^{\circ }$ . The vertical p.d.f. is almost indistinguishable from the reference case without inclination. The p.d.f. of the horizontal fluctuations is a narrow peak centred around zero and involves fluctuations that are negligible compared to those observed in experiments. We can therefore conclude that wake inclinations have no significant effect on the p.d.f.s and we will no longer consider them.

Let us sum up about the p.d.f.s of the fluctuations resulting from the summation of randomly distributed wakes. For a given bubble rise velocity and a gas volume fraction, the model essentially depends on a single parameter that characterizes the wake volume: $\tilde{L}\tilde{w}^{2}$ . In the relevant range of gas volume fractions in which the model is expected to be relevant ( $0.3\,\%\leqslant \unicode[STIX]{x1D6FC}\leqslant 20\,\%$ ), an exponential tail is observed for upward fluctuations, while neither significant horizontal fluctuations nor vertical downward fluctuations are generated.

2.3 Contribution of the potential flow

In the flow regime considered here, the shape of a rising bubble can be approximated as an oblate ellipsoid of revolution. In the close region located above and next to the bubble, i.e. out of the wake, the velocity field is well described by the potential flow around a single rising bubble (Ellingsen & Risso Reference Ellingsen and Risso2001; Risso & Ellingsen Reference Risso and Ellingsen2002). Under the potential flow approximation, the flow $\tilde{\boldsymbol{v}}(\boldsymbol{x})$ around an oblate ellipsoid of revolution in translation in the direction of its axis is known analytically (Lamb Reference Lamb1932, § 108) and depends only on the aspect ratio $\unicode[STIX]{x1D709}$ between the major and the minor axes. In this section, we consider the liquid velocity fluctuations generated by the summation of such potential disturbances.

The statistical cumulants are still given by (2.3). In the simplest case of a spherical bubble ( $\unicode[STIX]{x1D709}=1$ ), the variances of the three components of the velocity fluctuations have been determined by Biesheuvel & Wijngaarden (Reference Biesheuvel and Van Wijngaarden1984) (see also Lance & Bataille Reference Lance and Bataille1991) as

(2.9) $$\begin{eqnarray}\left.\begin{array}{@{}l@{}}\tilde{\unicode[STIX]{x1D70E}}_{pot\,x}^{2}={\textstyle \frac{3}{20}}\,\unicode[STIX]{x1D6FC},\\ \tilde{\unicode[STIX]{x1D70E}}_{pot\,y}^{2}={\textstyle \frac{3}{20}}\,\unicode[STIX]{x1D6FC},\\ \tilde{\unicode[STIX]{x1D70E}}_{pot\,z}^{2}={\textstyle \frac{1}{5}}\,\unicode[STIX]{x1D6FC}.\end{array}\right\}\end{eqnarray}$$

Because of the very complex expression of $\tilde{\boldsymbol{v}}(\boldsymbol{x})$ , obtaining the analytical expression of the statistical moments for any $\unicode[STIX]{x1D709}$ is of no use. It is however worth mentioning that the variance of the fluctuations is in any case a linear function of $\unicode[STIX]{x1D6FC}$ .

Statistically representative ensembles of samples of fluctuations induced by potential disturbances are obtained numerically by a method similar to that described previously for the wake disturbances. Figure 3 presents p.d.f.s that have been obtained for $\unicode[STIX]{x1D709}=2$ (case C, table 1) and various gas volume fractions. Figure 3(b) shows the p.d.f.s of the horizontal velocity fluctuations. As expected, they are symmetric. At low $\unicode[STIX]{x1D6FC}$ we observe a peak-shaped curve that is characteristic of fluctuations that each result from a single disturbance. Increasing $\unicode[STIX]{x1D6FC}$ , the maximum becomes rounded and exponential tails develop. Because the potential flow disturbance is closely localized near the bubble, the convergence towards a Gaussian function is very slow and only achieved for an unreasonable gas volume fraction of 100 %. Figure 3(a) shows the p.d.f.s of the vertical velocity fluctuations. Even if the streamlines of the potential flow are symmetric relative to the bubble mid-plane, the vertical velocity is not. The convergence towards a Gaussian function is even slower than for the horizontal fluctuations and the asymmetry needs a very high bubble number density to fade away. For intermediate gas volume fractions, exponential tails are still present.

Figure 3.

The p.d.f.s of liquid velocity fluctuations induced by potential flow disturbances at various gas volume fractions for a bubble aspect ratio $\unicode[STIX]{x1D709}=2$ (case C in table 1). The light (red) curve shows the Gaussian solution valid at large bubble number density for the case $\unicode[STIX]{x1D6FC}=100\,\%$ .

The presence of exponential tails in the p.d.f.s of velocity fluctuations resulting from the summation of randomly distributed disturbances seems to be a general feature, independent of the nature of the disturbances. It was previously found in the Stokes regime (Drazer et al. Reference Drazer, Koplik, Khusid and Acrivos2002; Abbas et al. Reference Abbas, Climent, Simonin and Maxey2006). It is found here for both wake and potential disturbances. It is observed for intermediate number densities between the very dilute regime, in which each fluctuation results from a single disturbance, and the very concentrated regime, in which each fluctuation results from a very large number of individual disturbances and where the p.d.f. is Gaussian. The corresponding range of gas volume fractions depends on the spatial extent of the disturbance around each bubble, which is larger for a Stokes disturbance than for a wake-induced disturbance, which is itself larger than for a potential disturbance. It also depends on the asymmetry of the disturbance, which delays the convergence towards the asymptotic solution at large concentration. Both wake and potential disturbances do not converge to a Gaussian distribution for volume fractions that can be achieved when the finite size of the bubble is taken into account. At large Reynolds number, the summation of individual disturbances thus generates exponential tails in the range of gas volume fractions where they are observed in experiments, from a few tenths of one per cent to a few tens of per cent.