Aurore Loisy
email: aurore.loisy@bristol.ac.uk

Aurore Naso

Peter D M Spelt

Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction

Various expressions have been proposed previously for the rise velocity of gas bubbles for homogeneous steady bubbly flows, generally a monotonically decreasing function of the bubble volume fraction. For suspensions of freely moving bubbles, some of these are of the form expected for ordered arrays of bubbles, and vice versa, as they do not reduce to the behaviour expected theoretically in the dilute limit. The microstructure of weakly inhomogeneous bubbly flows not being known generally, the effect of microstructure is an important consideration. We revisit this problem here for bubbly flows at small to moderate Reynolds number values for deformable bubbles, using direct numerical simulation and analysis. For ordered suspensions, the rise velocity is demonstrated not to be monotonically decreasing with volume fraction due to cooperative wake interactions. The fore-and-aft asymmetry of an isolated ellipsoidal bubble is reversed upon increasing the volume fraction, and the bubble aspect ratio approaches unity. Recent work on rising bubble pairs is used to explain most of these results; the present work therefore forms a platform of extending the former to suspensions of many bubbles. We adopt this new strategy also to support the existence of the oblique rise of ordered suspensions, the possibility of which is also demonstrated analytically. Finally, we demonstrate that most of the trends observed in ordered systems also appear in freely evolving suspensions. These similarities are supported by prior experimental measurements, and attributed to the fact that free bubbles keep the same neighbours for extended periods of time.

Introduction

Bubble columns are widely employed in industry because they can offer excellent heat and mass transfer characteristics without requiring any additional mechanical stirring. They are conceptually simple: a gas is sparged at the bottom of a liquid-filled vessel, and the bubbles rise under the effect of buoyancy. High transfer rates can then be attained owing to the increased contact area between the gas and the liquid phases, and to the liquid agitation induced by the bubbles' motion. Since a reliable prediction of bubble residence time and available interfacial area is crucial for an accurate and successful design of industrial columns, the understanding of bubble flow dynamics is essential.

In this context, a fundamental problem is the accurate prediction of the average bubble rise velocity in a statistically homogeneous buoyancy-driven flow of monodisperse bubbles, with the intention of using the results to systems with weak confinement and/or large-scale gradients. The average bubble rise velocity relative to the average velocity of the entire suspension is termed herein the drift velocity and the magnitude of the drift velocity is denoted by U . In the dilute limit (wherein the bubble volume fraction is vanishingly small), if no clustering occurs, bubbles behave as if they were isolated, and the drift velocity approaches the terminal velocity U 0 of a single bubble in unbounded liquid under otherwise the same conditions, for which a number of correlations is available [START_REF] Clift | Bubbles, drops, and particles[END_REF][START_REF] Loth | Quasi-steady shape and drag of deformable bubbles and drops[END_REF]. As volume fraction increases, U generally departs from U 0 . We therefore define G(φ) = U/U 0 , where φ is the gas volume fraction; the dependencies of G on various other dimensionless groups (termed herein as 'flow conditions') are suppressed in the notation for G for brevity, but are not ignored.

A common form of the dependency on volume fraction used in empirical correlations is the Richardson-Zaki expression G(φ) = (1 -φ) n [START_REF] Richardson | Sedimentation and fluidisation: Part I[END_REF][START_REF] Ishii | Drag coefficient and relative velocity in bubbly, droplet or particulate flows[END_REF], where n is an empirical parameter that depends on flow conditions. Various experiments have been carried out previously to determine G(φ) by injecting bubbles of millimetric size in a vertical column filled with clean water. In their respective experiments, [START_REF] Zenit | Measurements of the average properties of a suspension of bubbles rising in a vertical channel[END_REF] found that their measurements could be approximately described by the Richardson-Zaki formula with exponent n = 2.8, [START_REF] Garnier | Measurement of local flow characteristics in buoyancy-driven bubbly flow at high void fraction[END_REF] obtained instead G(φ) = 1 -φ 1/3 , and [START_REF] Colombet | Dynamics and mass transfer of rising bubbles in a homogenous swarm at large gas volume fraction[END_REF] fitted their data by G(φ) = (0.28+0.72 exp(-15φ)) 1/2 . Experimental complexities may have arisen that could explain these differences. For example, it is challenging to purify water of surfactants, to rule out wall effects, and to obtain a truly monodisperse suspension. A main inconvenience is that usually the bubble diameter cannot be kept constant if the volume fraction is changed, and the terminal velocity (and shape) of a single bubble strongly depends on the bubble size [START_REF] Clift | Bubbles, drops, and particles[END_REF]). This renders difficult disentangling the dependency of drift velocity on hydrodynamic interactions -along with the microstructure -from that on bubble size, if the latter is varied simultaneously with volume fraction.

Theoretical predictions of the drift velocity of bubbles are available for asymptotic and ideal systems. The drift (or rather, sedimentation) velocity derived by [START_REF] Batchelor | Sedimentation in a dilute dispersion of spheres[END_REF] for rigid particles in Stokes flow, wherein a uniformly random sedimenting suspension of particles is considered, has been generalized to droplets and bubbles, yielding G(φ) = 1 -4.44φ + O(φ 2 ) for bubbles whose viscosity can be ignored compared to that of the liquid [START_REF] Keh | Slow motion of multiple droplets in arbitrary threedimensional configurations[END_REF]; earlier work cited therein used various approximations for the mobility matrix). This linear correction includes a contribution -φ due to the back flow induced by the bubble motion, to maintain the mixture velocity of the entire suspension. It also includes a contribution -3φ from a volume-exclusion effect: a test bubble drags some liquid along with it which is compensated by a downflow further away; other bubbles cannot access the fluid in the direct vicinity of the test bubble (as this would require bubbles to overlap) and therefore sample a net fluid downflow (see [START_REF] Batchelor | Sedimentation in a dilute dispersion of spheres[END_REF] and for the application to bubbles/droplets see [START_REF] Wacholder | Sedimentation in a dilute emulsion[END_REF], albeit that an approximation for the mobility was used by the latter). The remainder of the linear term in G(φ) is due to near-bubble interactions.

The drift velocity strongly depends already on the microstructure in the dilute limit, though. If the probability of finding a particulate near a test particulate is uniform outside the excluded volume, as is appropriate for colloidal dispersions [START_REF] Glendinning | A pairwise additive description of sedimentation and diffusion in concentrated suspensions of hard spheres[END_REF] and as is assumed in the work cited above, interactions with nearby particulates result in G(φ) ∼ 1 -O(φ), whereas for a microstructure wherein a dominant contribution comes from particulates at a distance that corresponds to the mean separation (V p /φ) 1/3 , with V p the particulate volume, one expects G(φ) = 1 -O(φ 1/3 ) (e.g., [START_REF] Davis | Sedimentation of Noncolloidal Particles at Low Reynolds Numbers[END_REF]). A well-known example of the latter is a microstructure wherein particulates are arranged in a regular, ordered array. For bubbles arranged in an ordered array, G(φ) = 1 -bφ 1/3 + O(φ 2 ), where the coefficient b is known for several cubic arrays, and the O(φ) term cancels if the gas viscosity is ignored compared to the liquid viscosity [START_REF] Sangani | Creeping flow through cubic arrays of spherical bubbles[END_REF][START_REF] Sangani | Sedimentation in ordered emulsions of drops at low Reynolds numbers[END_REF]). The microstructure is not known a priori, therefore these limiting cases of random and ordered arrays may provide the means to determine the possible magnitude of the effects of order of the microstructure.

Besides Stokes flows, a weakly-viscous theory based on potential-flow interactions has also been developed [START_REF] Spelt | Properties and averaged equations for flows of bubbly liquids[END_REF]. Exceptionally, the probability density function for a configuration of identical spherical bubbles is known in that case. The pair probability shows a peak around the horizontal plane and a deficit for bubbles aligned vertically, both vanishing at large separation, if the drift velocity is not small compared to bubble velocity fluctuations, as may be expected from a Bernoulli effect. The averaged drag coefficient could therefore be determined analytically in the dilute limit whilst accounting for the microstructure. In the present notation, this yielded G(φ) = 1 -17 8 + 9 20 A φ + O(φ 2 ), where A is the ratio of U 2 and the root-mean-square bubble velocity; a Padé approximation for use beyond the dilute limit is presented in [START_REF] Spelt | Properties and averaged equations for flows of bubbly liquids[END_REF], and an extension for non-spherical bubbles at small A is also available [START_REF] Kushch | Finite Weber number motion of bubbles through a nearly inviscid liquid[END_REF]. The dependency on A enters there because it affects the pair-probability density function in that analysis, this being nearly isotropic at low A whilst showing a preference for bubbles rising nearly side by side at large A. Thus, in this description, the microstructure is determined from the significance of the drift velocity compared to that of bubble velocity fluctuations, the latter being the result of local inhomogeneities such as the shear rate (see [START_REF] Spelt | Properties and averaged equations for flows of bubbly liquids[END_REF].

Beyond these theoretical approaches, direct numerical simulations (DNS) of unbounded buoyancy-driven flows, in the sense of resolving the full Navier-Stokes equations coupled with the bubble dynamics and deformation, have been performed in prior work for cubic domains that contain a finite number of freely-moving bubbles, subject to periodic boundary conditions. Conveniently, this setup allows variation between microstructures. On the one hand, for a given volume fraction, using a large number of bubbles in the unit cell is of interest as a model of real suspensions, although convergence with number of particulates would have to be verified, since for Stokes interactions this may be slow [START_REF] Phillips | Hydrodynamic transport properties of hardsphere dispersions. I. Suspensions of freely mobile particles[END_REF], or a dependency on system size may even persist [START_REF] Guazzelli | Fluctuations and Instability in Sedimentation[END_REF]; in the studies cited below, typically O(10) to O(100) bubbles are used in a unit cell, and the effect of system size is found to be small. In the other extreme, the special case of one freely-rising bubble in the unit cell, one recovers a simple cubic array. We shall herein refer to this setup with more than one bubble in the unit cell as a free array, and to that with one bubble in the cell as an ordered array.

For spherical bubbles rising at low (O(1)) Reynolds numbers, the DNS results of [START_REF] Esmaeeli | Direct numerical simulations of bubbly flows. Part 1. Low Reynolds number arrays[END_REF], suggest that G(φ) for free arrays may be similar to that predicted for ordered arrays, but the system studied was concluded to be too small to draw definitive conclusions. [START_REF] Bunner | Effect of bubble deformation on the properties of bubbly flows[END_REF] found that their results at moderate (O(10)) Reynolds numbers could be represented by G(φ) = 1 -φ 1/3 for spherical bubbles in free arrays, and a Richardson-Zaki expression with exponent n = 3 for deformable bubbles, but that no theoretical justification for these scalings could be offered, beyond an observed difference in preferential spatial configurations of bubbles, discussed further below. Also, the expressions for G(φ) inferred in these pioneering studies could only be fitted from just a few different values of the volume fraction comprised between 2 and 12 %, their validity outside this range seems unclear. Further results for free arrays at moderate Reynolds numbers were obtained by [START_REF] Yin | Lattice-Boltzmann simulation of finite Reynolds number buoyancydriven bubbly flows in periodic and wall-bounded domains[END_REF], for volume fractions ranging from 1 to 25 %, using a lattice-Boltzmann method. They used O(100) bubbles in a periodic cell, rather than O(10) in most of the early studies cited above, and imposed a spherical shape, facilitated by a force balance that includes the surface integral of the traction acting on the bubble and the buoyancy force. Their results demonstrate that G(φ) is not well fitted by the Richardson-Zaki formula, and they suggest that this is associated with the anisotropic microstructure of bubbly suspensions in this regime. [START_REF] Gillissen | A lattice Boltzmann study on the drag force in bubble swarms[END_REF] conducted similar simulations using a combination of lattice-Boltzmann and immersed boundary methods, and obtained G(φ) = 1-O(φ 1/3 ) for spherical bubbles rising at small to moderate Reynolds number, in qualitative agreement with earlier studies [START_REF] Esmaeeli | Direct numerical simulations of bubbly flows. Part 1. Low Reynolds number arrays[END_REF][START_REF] Bunner | Effect of bubble deformation on the properties of bubbly flows[END_REF].

From these prior studies, the microstructure, along with the drift velocity, is known to vary significantly with bubble Reynolds number and shape. For spherical bubbles rising at O(100) Reynolds number, strong preference for horizontal alignment is observed in the simulated pair probability [START_REF] Esmaeeli | A direct numerical simulation study of the buoyant rise of bubbles at O(100) Reynolds number[END_REF], in agreement with the trends reviewed above for the idealized potential-flow interactions. The anisotropy in microstructure is larger than that observed in the experiments of [START_REF] Zenit | Measurements of the average properties of a suspension of bubbles rising in a vertical channel[END_REF], possibly because of bubble deformation: indeed the simulations of [START_REF] Esmaeeli | A direct numerical simulation study of the buoyant rise of bubbles at O(100) Reynolds number[END_REF] revealed that oblate ellipsoidal bubbles do not form horizontal rafts but instead are rather uniformly distributed. At O(10) Reynolds number, the dynamics of bubblebubble interactions is dominated by wake effects. A vertical pair of spherical bubbles changes its orientation to horizontal through a drafting-kissing-tumbling mechanism, resulting in preferential side by side alignment [START_REF] Esmaeeli | Direct numerical simulations of bubbly flows. Part 2. Moderate Reynolds number arrays[END_REF]Bunner & Tryggvason 2002a;[START_REF] Yin | Lattice-Boltzmann simulation of finite Reynolds number buoyancydriven bubbly flows in periodic and wall-bounded domains[END_REF], whereas deformable bubbles tend to organize in vertical structures owing to the reversed lift force which attract a bubble in the wake of its preceding neighbour [START_REF] Bunner | Effect of bubble deformation on the properties of bubbly flows[END_REF]. These effects decrease with decreasing Reynolds number, and nearly no preference is observed at O(1) Reynolds number [START_REF] Esmaeeli | Direct numerical simulations of bubbly flows. Part 1. Low Reynolds number arrays[END_REF][START_REF] Cartellier | Bubble-induced agitation and microstructure in uniform bubbly flows at small to moderate particle Reynolds numbers[END_REF].

In the DNS studies cited thus far, the microstructure is allowed to develop naturally. For use of the results in general flows, wherein microstructure can be affected by weak gradients, it is necessary to know the role of and sensitivity to the microstructure. DNS results for rising deformable bubbles in an ordered arrangement have been conducted by Sankaranarayanan et al. (2002), using a lattice-Boltzmann method. An empirical correlation of Richardson-Zaki form was used to represent the results for a vast range of flow regimes, albeit for 5 φ 25% and with the reservation that the correlation does not reduce to the analytical result discussed above for creeping flows of ordered arrays.

Despite their apparent artificiality, ordered arrays of bubbles, as well as their relevance to real bubbly flows, certainly deserve further investigation. Firstly, a number of prior simulations and experiments reviewed above for 1 Re 1000 found a bubble rise velocity scaling as φ 1/3 [START_REF] Esmaeeli | Direct numerical simulations of bubbly flows. Part 1. Low Reynolds number arrays[END_REF][START_REF] Bunner | Effect of bubble deformation on the properties of bubbly flows[END_REF][START_REF] Gillissen | A lattice Boltzmann study on the drag force in bubble swarms[END_REF][START_REF] Garnier | Measurement of local flow characteristics in buoyancy-driven bubbly flow at high void fraction[END_REF]. This scaling is the same as that obtained assuming a periodic arrangement of the bubbles, albeit under the Stokes flow approximation [START_REF] Sangani | Creeping flow through cubic arrays of spherical bubbles[END_REF]. The study of ordered arrays beyond the Stokes-flow limit is therefore of fundamental interest in order to connect theoretical, numerical, and experimental work. Secondly, prior experimental and numerical work on bubbly flows at moderate to high Reynolds number has shown that the magnitude of the bubble velocity fluctuations is substantially smaller than the bubble rise velocity (Bunner & Tryggvason 2002b;[START_REF] Esmaeeli | A direct numerical simulation study of the buoyant rise of bubbles at O(100) Reynolds number[END_REF][START_REF] Zenit | Measurements of the average properties of a suspension of bubbles rising in a vertical channel[END_REF][START_REF] Martinez-Mercado | Measurement of pseudoturbulence intensity in monodispersed bubbly liquids for 10 < Re < 500[END_REF], at least when the gas volume fraction remains below approximately 10 %. This further motivates a study of a representation of bubbly suspensions by ordered arrays (zero bubble velocity fluctuations). Thirdly, only experimental investigations can assess the relevance (or lack thereof) of the ordered model to describe real bubbly suspensions. By re-examining prior experimental data [START_REF] Garnier | Measurement of local flow characteristics in buoyancy-driven bubbly flow at high void fraction[END_REF][START_REF] Martinez-Mercado | Measurement of pseudoturbulence intensity in monodispersed bubbly liquids for 10 < Re < 500[END_REF][START_REF] Riboux | Experimental characterization of the agitation generated by bubbles rising at high Reynolds number[END_REF][START_REF] Colombet | Dynamics and mass transfer of rising bubbles in a homogenous swarm at large gas volume fraction[END_REF], we will show that available measurements support the idea that ordered arrays are indeed relevant to bubbly flows of practical interest.

In this paper, we investigate the ordered and free rise of bubbles at low and moderate Reynolds numbers over a wide range of volume fractions, using DNS and analysis; the problem statement is presented in section 2, and the numerical methods in section 3. The first objective (in section 4) is to determine the connection between the DNS results and theory for dilute ordered systems and, beyond the dilute limit, the connection between the DNS results and prior work on bubble pairs, which have been studied analytically (e.g., [START_REF] Harper | On bubbles rising in line at large Reynolds numbers[END_REF][START_REF] Harper | Bubbles rising in line: why is the first approximation so bad?[END_REF]), experimentally [START_REF] Katz | Wake-induced relative motion of bubbles rising in line[END_REF] and computationally [START_REF] Yuan | On the in-line motion of two spherical bubbles in a viscous fluid[END_REF][START_REF] Legendre | Hydrodynamic interactions between two spherical bubbles rising side by side in a viscous liquid[END_REF][START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF]. For that purpose, the first effects of inertia are determined analytically, and a comprehensive parametric DNS study is presented at low and moderate Reynolds number in section 4.1. The study includes ordered suspensions of strongly deformed (skirted) bubbles that have not been studied thus far. The results are summarized in a practical relation that covers a wide range of flows and reduces to analytical results in the dilute limit. In particular, the DNS and the theory in the dilute limit demonstrate that G(φ) is not monotonic at low φ; Bunner & Tryggvason (2002a) and [START_REF] Roghair | On the drag force of bubbles in bubble swarms at intermediate and high Reynolds numbers[END_REF] remarked that G(φ) increases at small volume fractions, but did not investigate this further. The surprising dynamics of ordered suspensions is then investigated in section 4.2. A steady oblique rise was observed in the DNS in some cases, as a precursor to unsteady (periodic and nonperiodic) behaviour. Oblique motion of two-dimensional square arrays of bubbles has been observed previously (Sankaranarayanan et al. 2002;Sankaranarayanan & Sundaresan 2002;[START_REF] Theodoropoulos | Coarse bifurcation studies of bubble flow lattice Boltzmann simulations[END_REF]), but remains essentially unexplained. We support the existence of oblique rise for ordered suspensions with analysis in the dilute limit, and demonstrate that the inferred lift coefficient is similar to that obtained for the rise of bubble pairs. The second objective of the present work is to revisit (in section 5) arrays of free bubbles in light of these findings for ordered arrays, and to contrast and compare these two systems. In section 5.2, numerical results are presented for freely evolving bubbly suspensions at small and intermediate volume fractions, and similarities with ordered systems are highlighted. These are supported by a comparison with experimental data in section 5.3. Finally, conclusions are presented in section 6.

Mathematical formulation

Problem statement

We consider an infinite, homogeneous, monodisperse suspension of bubbles rising under the effect of buoyancy in otherwise quiescent liquid. The density and viscosity of each fluid, as well as the surface tension, are assumed to be constant. The suspension is represented by the periodic repetition of a cubic unit cell containing a given number of bubbles. The gravity is aligned with a primary axis of the periodic array (due to the large number of parameters already involved in the problem, the influence of the orientation of gravity is not investigated here).

The behaviour of this system depends on nine parameters: the number of bubbles N b in the cell, the gas volume fraction φ, the gravitational acceleration g, the bubble volume or, more conveniently, its characteristic size d b defined as the diameter of the volumeequivalent sphere, and the physical properties of the two fluids, namely their densities (ρ d , ρ c ), their viscosities (µ d , µ c ), and the surface tension (γ). The subscripts d and c refer to the disperse (gaseous) and continuous (liquid) phases, respectively.

In addition to the gas volume fraction and to the number of bubbles, four independent dimensionless groups can be constructed from the remaining parameters. Two of these are the ratios of the gas density and viscosity to those of the surrounding liquid. These are usually very small and of the same order for most gas-liquid systems of practical interest. As a consequence their influence will not be investigated, and unless otherwise mentioned, these parameters will be set to ρ d /ρ c = 10 -3 and µ d /µ c = 10 -2 , which roughly corresponds to air bubbles in water. The last two dimensionless numbers are the Archimedes number

Ar = ρ c |ρ d -ρ c |gd 3 b µ c , (2.1)
or equivalently the Galilei number Ga = Ar 2 , and the Bond number (also known as the Eötvös number),

Bo = |ρ d -ρ c |gd 2 b γ . (2.2)
The Archimedes and Bond numbers can be defined a priori, without the knowledge of the bubble velocity, and are therefore traditionally employed to describe the macroscopic conditions of buoyancy-driven bubbly flow (numerical) experiments. At time zero, the bubbles are released from rest and start rising. The time evolution of the system is monitored through U (t), defined as the average drift velocity of the bubbles and computed at any time from

U = u d -u , (2.3)
where denotes a volume average over the entire unit cell and d denotes a volume average over the disperse phase only. In most situations U is parallel to gravity, so there is no need to distinguish between |U | and the vertical component of U . For simplicity, and unless mentioned otherwise, U is used to denote the (positive) vertical component of U . The drift velocity is used as the characteristic velocity scale to define the dynamic counterparts of the Archimedes and Bond numbers: the Reynolds number

Re = ρ c U d b µ c (2.4)
and the Weber number

We = ρ c U 2 d b γ = BoRe 2 Ar 2 , (2.5)
which compare the effects of inertia, viscosity and surface tension. In a system at equilibrium for vertical rise, the hydrodynamic force acting on a bubble, whose magnitude is denoted f , equals the buoyancy force. It follows that the Reynolds number is related to the Archimedes number through

C D = 4 3 Ar 2 Re 2 with C D = 8f πd 2 b ρ c U 2 , (2.6)
where C D is the drag coefficient.

Assuming that a (possibly quasi-)steady state is reached independently of the initial conditions (which is not necessarily the case, but we will come to that later), the (quasi-)steady average bubble drift velocity can be written as U = U (N b , φ, Ar , Bo). Similarly the (quasi-)steady average bubble shape, as described by a parameter χ (which will be specified later, typically an aspect ratio), reads χ = χ(N b , φ, Ar , Bo). Our first goal is to characterize U (φ, Ar , Bo) and χ(φ, Ar , Bo) when the bubbles have a fixed position relative to their neighbours (N b = 1), and to understand how the imposed flow conditions (Ar ,Bo) affect the dependency of these quantities on the volume fraction. Our second goal is to assess the effect of introducing additional degrees of freedom (N b > 1) in the system, and to compare the behaviour of freely evolving suspensions (sufficiently large N b ) with that of ordered suspensions (N b = 1).

Flow regimes

Since we want to assess the effect of volume fraction under various conditions of Archimedes and Bond numbers, it seems natural to refer to the limiting case of a single bubble released in an unbounded quiescent liquid under the same conditions. At steadystate, this bubble is characterized by its shape (and an associated aspect ratio χ 0 ), and its terminal velocity U 0 , usually expressed in the form of a terminal Reynolds number Re 0 = ρ c U 0 d b /µ c . The subscript 0 will be used hereinafter when an isolated bubble is considered.

A rather general description of the equilibrium state reached by a buoyancy-driven bubble is given in the shape regime diagram of [START_REF] Grace | Shapes and velocities of bubbles rising in infinite liquids[END_REF]. This diagram splits the (Bo, Re 0 ) parameter space in a number of subregions and maps them onto the corresponding shape regimes. It also provides a graphical correlation between the Bond number, the Reynolds number, and the Morton number Mo = (|ρ d -ρ c |gµ 4 c )/(ρ 2 c γ 3 ) = Bo 3 /Ar 4 , which is often used in experimental work in place of the Archimedes number. In a simplistic manner, the terminal Reynolds number increases (non-linearly) with the Archimedes number, while the bubble departs from a spherical shape as the Bond number increases.

We considered nine different cases defined by the pair (Ar ,Bo). A complete description of these cases and of the corresponding flow regimes is provided in table 1. They cover Reynolds numbers ranging from 0 to 60 and several shape regimes: spherical (cases "S"), ellipsoidal (cases "E"), and dimpled ellipsoidal-cap (case "C"). The parameters for case C correspond to a single-bubble experiment of [START_REF] Bhaga | Bubbles in viscous liquids: shapes, wakes and velocities[END_REF], which was later reproduced numerically by [START_REF] Hua | Numerical simulation of 3D bubbles rising in viscous liquids using a front tracking method[END_REF]. The terminal Reynolds number and shape of the equivalent isolated bubble have therefore been determined directly from their data. For the other cases, the single-bubble terminal Reynolds number and aspect ratio have been estimated using the correlations for spherical and ellipsoidal bubbles recommended in the review of [START_REF] Loth | Quasi-steady shape and drag of deformable bubbles and drops[END_REF].

Governing equations

In both phases the fluid motion is governed by the incompressible Navier-Stokes equations

∇ • u n = 0, (2.7) ∂ ρ n u n ∂t + ∇ • ρ n u n u n = ∇ • T n + G n where T n = -p n I + µ n (∇u n + ∇u T n ),
(2.8) where n = {c, d} is used here to denote either phases, u is the velocity field, T the stress tensor, I denotes the identity tensor, p the pressure field, and G is the sum of external forces per unit volume, given by G n = (ρ n -ρ )g. The first term in this last expression, ρ n g, is the weight of a unit volume of fluid. Because the suspension we consider is infinite, i.e. not bounded by walls, an additional body force -ρ g, with ρ the system average density, is required to prevent the entire system from accelerating in the downward vertical direction †.

These equations are coupled through the appropriate boundary conditions at the ). The shapes, Re0, We0 = BoRe 2 0 /Ar 2 , and χ0 of an isolated buoyancy-driven bubble at steady-state are also given. Shapes are predicted by the diagram of [START_REF] Grace | Shapes and velocities of bubbles rising in infinite liquids[END_REF]. The values of Re0 are estimated from the correlation of [START_REF] Mei | A note on the history force on a spherical bubble at finite Reynolds number[END_REF] for spherical bubbles (cases S0 to S5) and from the correlation of [START_REF] Loth | Quasi-steady shape and drag of deformable bubbles and drops[END_REF] for ellipsoidal bubbles (cases E1 and E2); the experimental value measured by [START_REF] Bhaga | Bubbles in viscous liquids: shapes, wakes and velocities[END_REF] is reported for case C. The aspect ratio χ0 is estimated from the correlation of [START_REF] Loth | Quasi-steady shape and drag of deformable bubbles and drops[END_REF] for all cases except case C, for which it is directly measured from visualizations of [START_REF] Hua | Numerical simulation of 3D bubbles rising in viscous liquids using a front tracking method[END_REF].

interface. In the presence of viscous effects a no-slip condition is applied, which, combined with the absence of mass flux across the interface, leads to

[u] = 0 (2.9)
where [X] = X c -X d denotes the jump of a variable across the interface. Neglecting any variation of surface tension along the interface, the shear stress is continuous across it, and the jump of normal stress is balanced by the curvature force per unit area:

[n • T ] = γκn, (2.10)
where n is the unit vector normal to the interface and directed outward from the bubbles, and κ is the interface curvature defined by κ = ∇ • n (e.g., Tryggvason et al. (2011), Appendix A.4). This set of equations is solved numerically within a periodic unit cell using the methods described in the next section.

Numerical methods

A comprehensive description of our numerical approach is provided in appendix A.1. A brief overview of its salient features is given hereinafter.

Our approach relies on the one-fluid formulation of the governing equations. In this formulation, the different fluids are treated as a single phase with discontinuous density and viscosity, and surface tension is incorporated as a singular source term. This results in the standard continuum surface force model of [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]. To circumvent numerical difficulties due to the introduction of discontinuous and singular functions, the interface is given a finite thickness proportional to the grid spacing. Surface tension is therefore treated as a volume force distributed over several mesh points, and material properties vary continuously from one phase to the other.

The incompressible Navier-Stokes equations are integrated in their one-fluid form by a projection method [START_REF] Chorin | Numerical Solution of the Navier-Stokes Equations[END_REF]), and the moving interface separating the two fluids is captured by a level-set method [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Sussman | A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow[END_REF]. The velocity field is then solution of the system of equations:

∂ ρu ∂t + ∇ • ρuu = -∇p + ∇ • µ(∇u + ∇u T ) + (ρ -ρ )g -γκ∇H ε (3.1) ∇ • u = 0, (3.2)
where the interface curvature is calculated from (e.g., [START_REF] Prosperetti | Computational Methods for Multiphase Flow[END_REF], section 3.5)

κ = ∇ • ∇ψ |∇ψ| , (3.3)
the variable density and viscosity are given by

ρ = H ε ρ c + (1 -H ε )ρ d , µ = H ε µ c + (1 -H ε )µ d , (3.4)
and H ε denotes the smoothed Heaviside function

H ε (ψ) =        1 if ψ > ε, 0 if ψ < -ε, 1 2 1 + ψ ε + 1 π sin πψ ε if |ψ| ε, (3.5)
where ε is half the interface thickness (ε = 1.5∆x, with ∆x the grid spacing). ψ denotes the level-set function, positive in the continuous phase and negative in the disperse one. It is the solution of the following advection equation, in which the zeroth-order approximation of the additional source term proposed by [START_REF] Sabelnikov | Modified level set equation and its numerical assessment[END_REF] is embedded:

∂ ψ ∂t + u • ∇ψ = A(u, ψ) ψ, with A(u, ψ) = ∇ i ψ∇ i u j ∇ j ψ (3.6) (A(u, ψ)
is the local zeroth-order approximation of the source term in the region close to the interface). The level-set function is reinitialized as a signed distance function at each time step using the procedure devised by [START_REF] Russo | A remark on computing distance functions[END_REF]. These modifications yield better numerical efficiency and accuracy compared to the original level-set method.

The major drawback of basic level-set methods is their poor ability to conserve the mass (volume) of each phase. When using high-order schemes and sufficient resolution, the volume change between two successive timesteps is often negligible. It is, however, not exactly zero, and may even become substantial when accumulated over very long integration times. For this reason, we enforce volume conservation using the correction proposed by [START_REF] Sussman | A computational study of the spreading of oil underneath a sheet of ice[END_REF] (also used by, e.g., [START_REF] Spelt | Shear flow past two-dimensional droplets pinned or moving on an adhering channel wall at moderate Reynolds numbers: a numerical study[END_REF]): at the end of each timestep, the iso-contours of the level-set function are slightly shifted such that the volume of each phase is conserved exactly. We demonstrate in appendix A.2 that the loss of accuracy induced by this correction is negligible compared to the overall numerical error made in the interfacial region, and has therefore no adverse effect on the flow dynamics.

Our time integration algorithm is based on third-order and second-order TVD Runge-Kutta schemes for the level-set advection and reinitialization equations, respectively, and on a mixed Crank-Nicolson/third-order Adams-Bashforth scheme for the Navier-Stokes equations. The spatial discretization is carried out using a standard finite difference/finite volume discretization on a uniform Cartesian staggered grid: fifth-order WENO schemes are used for advection terms, and second-order centered schemes are used otherwise. Periodic boundary conditions are imposed at the bounds of the computational domain.

A standard benchmark test for the simulation of disperse flows consists in comparing the terminal velocity and shape of an isolated bubble to those obtained experimentally in various regimes. We do not simulate here the rise of a single bubble in an unbounded fluid, but of an array of bubbles. A tempting idea to approach this ideal situation would be to introduce a single bubble in a very large unit cell, so that the influence of periodicity could be neglected. We will see in section 4 that even at very low volume fractions (very large domains, in the limit of what is computationally feasible), the bubbles' rise velocity and shape are still significantly affected by their interactions, making such a comparison to experiments irrelevant. The code has therefore been validated against available solutions for regular (N b = 1) and free (N b > 1) arrays of rising bubbles.

We reproduced the numerical simulations of [START_REF] Esmaeeli | Direct numerical simulations of bubbly flows. Part 2. Moderate Reynolds number arrays[END_REF] who considered arrays of deformable bubbles rising at moderate Reynolds number (case E1). As shown in figure 1, the transient rise of both regular and free arrays of bubbles is reproduced accurately by our code. We also compared our simulations with the theory of [START_REF] Sangani | Sedimentation in ordered emulsions of drops at low Reynolds numbers[END_REF] for cubic arrays of spherical bubbles in the creeping flow limit (case S0). The evolution of the steady drift velocity with volume fraction is shown in figure 2. Excellent agreement between the numerical and analytical solutions has been obtained. The effect of resolution is also shown in the same figures. Additional benchmark and sensitivity tests can be found in appendix A.2. Further comparison against prior work is included where prior work is available in subsequent sections.

Grid convergence tests have been carried out systematically for each of the cases reported in table 1 in an ordered array configuration and for one value of the volume fraction. A resolution of 20 grid cells per bubble diameter was found to be sufficient for all regimes except for case S5, which requires a resolution of 30 grid cells per diameter because of the higher Reynolds numbers associated to this regime, and for case C, for which a resolution of 60 grid cells per diameter was needed for capturing the thin skirts of the bubbles. With these resolutions, the error in the steady rise velocity due to the grid spacing is not larger than 2 % in case C and 1 % in the other cases. For a given case, the same resolution is used for all volume fractions, and for both ordered and free arrays. The choice of the time step is constrained by the condition of numerical stability, and the error due to the time discretization is smaller than that due to the spatial one. q q q q q q q q q q q q q Figure 2. Steady drift velocity of an ordered array of spherical bubbles normalized by that of an isolated bubble in the Stokes flow regime (case S0) as a function of volume fraction, compared with the analytical solution of [START_REF] Sangani | Sedimentation in ordered emulsions of drops at low Reynolds numbers[END_REF]. The effect of resolution is shown for φ 1/3 = 0.4 (d b is the bubble volume-equivalent diameter, ∆x the grid spacing).

Ordered arrays

We examine in this section the dynamics of cubic arrays of deformable bubbles ("ordered arrays") in the presence of liquid inertia. The main objective here is to connect DNS results, theoretical analysis for dilute systems, and prior work on bubble pairs. Specifically, the effect of volume fraction on the rise velocity and shape of deformable bubbles in a simple cubic array is revisited here. Direct numerical simulations have been performed for the nine sets of flow conditions summarized in table 1. For each of these cases, the volume fraction φ = π/6(d b /h) 3 (where h is the linear size of the unit cell of the array) was varied from 0.1 to 30 % by changing the size of the computational domain (i.e., the lattice spacing) while keeping the bubble size constant.

Initially, both fluids were at rest, the bubbles were spherical, and gravity was switched on at time zero. After a transient regime, various forms of bubble motion could be observed: steady vertical rise, steady oblique rise, or unsteady oblique rise. The steady vertical rise is first examined in section 4.1. Other types of motions are then discussed in section 4.2.

Simulations were run until the bubble drift velocity became either constant or statistically stationary. This steady state is independent of the initial oblateness of the bubbles, and is reached when the velocity disturbances induced by bubbles' motion have diffused in all directions throughout the liquid, i.e., in a time of order O(h 2 ρ c /µ c ). As a consequence, from a numerical point of view, the investigation of small volume fractions (large domain sizes) is limited both by the needed number of grid points (∼ h 3 ) and by the computation time (∼ h 2 ).

Steady vertical rise of bubbles

The cubic lattice of bubbles is not only convenient from a computational point of view, it is also attractive from a theoretical standpoint since the solution only needs to be determined in a unit cell. When the bubbles rise steadily along straight paths parallel to an axis of the periodic array (as is the case in most of the cases presented here, since gravity is oriented along a lattice axis), the symmetries of the problem greatly simplify the analysis. We determine in this context an analytical expression accounting for the first effect of inertial interactions in cubic arrays of spherical bubbles (at small Reynolds numbers). Outside this narrow range of validity, the influence of the volume fraction on the steady rise velocity and shape of deformable bubbles will be determined from our numerical simulations.

Spherical bubbles at low to moderate Reynolds number

The correction to the drift velocity due to finite volume fraction in the Stokes-flow regime has been determined by [START_REF] Sangani | Sedimentation in ordered emulsions of drops at low Reynolds numbers[END_REF] for cubic arrays of spherical fluid particles (bubbles or drops). The first term arises from a point-force approximation of the particles and reads, for a simple cubic array,

U U 0,Stokes -1 = -1.1734µ * φ 1/3 + O(φ), (4.1)
where U 0,Stokes is the terminal velocity of a spherical fluid particle translating through an unbounded ambient fluid in Stokes-flow conditions [START_REF] Hadamard | Mouvement permanent lent d'une sphere liquide et visqueuse dans un liquide visqueux[END_REF][START_REF] Rybczynski | Uber die fortschreitende Bewegung einer flussigen Kugel in einem zahen Medium[END_REF]):

U 0,Stokes = 1 12 |ρ c -ρ d |gd 2 b µ * µ c , with µ * = µ c + 3/2µ d µ c + µ d . (4.2)
The case of a rigid sphere (µ d /µ c → ∞) is recovered as µ * → 3/2, whereas the case of a clean bubble (µ d /µ c → 0) corresponds to the limit µ * → 1. Numerical simulations have been carried out for spherical bubbles rising at very small Reynolds numbers (case S0, Re 0 = 1.94 × 10 -3 ), and excellent agreement with the theory has been obtained over a wide range of volume fractions, as shown in figure 2. Note that the relation between U and φ 1/3 seems linear even at high volume fraction: indeed even though the analytical solution of [START_REF] Sangani | Sedimentation in ordered emulsions of drops at low Reynolds numbers[END_REF] includes O(φ) and O(φ 2 ) terms, for clean bubbles the O(φ) correction is zero and the O(φ 2 ) correction is negligible compared to the O(φ 1/3 ) term (whereas these corrections are substantial for solid particles).

For small but non-zero Reynolds numbers, the Stokes equations are still valid near and inside the fluid particles, but should be replaced by the Oseen equations farther away, since inertial effects become comparable to viscous ones at distances from the particle of order O(d b /Re). The first correction to the drag force arising from inertial effects has been determined by [START_REF] Hill | The first effects of fluid inertia on flows in ordered and random arrays of spheres[END_REF] for a cubic array of solid spheres. The extension of their result to bubbles and drops is straightforward and is provided in appendix B. We show there that the correction to the bubble drift velocity due to liquid inertia and hydrodynamic interactions can be approximated at any φ 1 by

U U 0,Stokes -1 ≈ - 1 8 µ * Re -1.1734µ * φ 1/3 + 25 8 µ * Reφ 1/3 Re + 25φ 1/3 . (4.3)
The first term accounts for the effect of liquid inertia on an isolated bubble, the second term results from Stokes interactions, and the last term captures the effect of inertial interactions. The significance of each of these terms as a function of volume fraction can be understood as follows.

At zero volume fraction, the drag exerted on a single bubble normalized by the Stokes drag increases linearly with the Reynolds number [START_REF] Brenner | The resistance to a particle of arbitrary shape in translational motion at small Reynolds numbers[END_REF]. This results in the negative correction to the drift velocity

U 0 U 0,Stokes -1 = - 1 8 µ * Re, (4.4)
where U 0 is the terminal velocity of the isolated bubble. At small volume fraction, inertial interactions result in a positive O(φ 1/3 ) correction which overwhelms the negative O(φ 1/3 ) present analysis present simulation isolated bubble prior simulation + RZ fit q q q q q q q q q q q q q q q

φ 1 3 U U 0,

Stokes-flow correction

U -U 0 U 0,Stokes ≈ -1.1734µ * φ 1/3 + 25 8 µ * φ 1/3 ≈ 2.0µ * φ 1/3 when φ 1/3 Re, (4.5)
so the net result is a drift velocity that increases with φ 1/3 . At large volume fraction (with respect to the Reynolds number), the drift velocity correction due to hydrodynamic interactions reads

U -U 0 U 0,Stokes = -1.1734µ * φ 1/3 + 1 8 µ * Re 1 - Re 25φ 1/3
when Re

φ 1/3 1. (4.6)
The O(Re) contribution from inertia is negligible compared to the Stokes O(φ 1/3 ) correction: the drift velocity therefore overall decreases linearly with φ 1/3 , as for creeping flows.

The drift velocity U can be computed, for any φ 1, by finding the positive root of (4.3) (quadratic in U ). The solution for various Archimedes numbers is shown in figure 3(a). Note that higher Archimedes number corresponds to higher isolated-bubble Reynolds number: it can be shown from (2.1), (2.4), (4.2) and (4.4) that in this regime

Re 0 = 1 12 Ar 2 µ * 1 + Ar 2 96 -1 . (4.7)
The non-monotonicity of the function U (φ)/U 0,Stokes at finite Archimedes number contrasts with the case of Stokes flow, for which this function is strictly decreasing. This behaviour results from the competition between "cooperative" long-range inertial interactions, which increase the drift velocity, and "hindering" viscous interactions which reduce it. At small volume fraction, inertial effects dominate, whereas at large volume fraction the liquid is more confined, inertial forces therefore cannot prevail over viscous ones and a Stokes-flow behaviour is recovered. The Oseen approximation is limited to Re < 1, which for an isolated clean bubble roughly corresponds to Ar < 3.5. For Ar = 3, the maximum of U/U 0,Stokes is obtained for q q q q q q q q q q q q q q Figure 4. Steady drift velocity of an ordered array of spherical bubbles, normalized by the terminal velocity of an isolated bubble in Stokes flow conditions, as a function of volume fraction for a large range of Archimedes numbers. Symbols: DNS; --: analytical Oseen-flow solutions obtained from (4.3) for small Ar ; ----: numerical fits of the form of (4.8) matching DNS and isolated-bubble data for larger Ar . In case S5, for which the bubbles motion is not steady and/or not parallel to gravity, the crosses are time-averaged vertical drift velocities. φ = 6×10 -6 . Direct numerical simulation of such a small volume fraction is prohibitively expensive, so our analysis cannot be confirmed by numerical experiments in its expected range of validity. Nevertheless a comparison between the solution obtained from (4.3) and DNS for Ar = 5.03 is shown in figure 3(b), together with the numerical data obtained by Sankaranarayanan et al. (2002) for the identical flow regime using the lattice Boltzmann method. In their study, they found that the effect of volume fraction could be captured by a Richardson-Zaki type of (empirical) correlation Re = 1.58(1 -φ) 4.72 . Although their data is well-fitted by this relation over the narrow range of volume fractions they investigated (0.05 < φ < 0.12), our DNS results show that this expression cannot be used to extrapolate the effect of volume fraction outside this range. In addition, their correlation gives a drift velocity at φ = 0 that differs from the terminal velocity of an isolated bubble by more than 10 %. In contrast, the functional dependency of the drift velocity on volume fraction given by our analysis is in very good agreement with numerical simulations; the modest difference at small volume fractions arises from the limitation of Oseen theory to Reynolds numbers less than unity: for an isolated bubble rising in still liquid, the Oseen-flow solution yields Re 0 = 1.66 whereas the empirical correlation of [START_REF] Mei | A note on the history force on a spherical bubble at finite Reynolds number[END_REF] gives Re 0 = 1.80.

We now turn to ordered arrays of spherical bubbles rising at moderate Reynolds numbers. The bubble drift velocity has been determined for Archimedes numbers ranging from 0 to 40 (case S0 to S5). The numerical results for U (φ)/U 0,Stokes are shown in figure 4 (symbols) together with those of our analysis for small Archimedes numbers (solid lines). It is remarkable that the evolution of the drift velocity with volume fraction for Archimedes numbers up to approximately 30 is consistent with the Oseen-flow analysis carried out for Archimedes numbers that are, at best, O(1). In particular, for Ar = 27.2 (case S4), the predicted increase of the drift velocity at low volume fraction is confirmed numerically. For Archimedes numbers greater than 30 (case S5), the drift velocity of a cubic array of spherical bubbles does not necessarily remain parallel to gravity; we postpone discussion of this to section 4.2.

Figure 5 shows the vertical component of the liquid velocity in a vertical symmetry plane passing through the center of a bubble at φ = 0.2 %. The first row corresponds to spherical bubbles (Bo = 0.38) with increasing Archimedes numbers from left to right. It reveals that the region of liquid dragged up by each bubble extends quite far downstream. Since the bubble motion is parallel to a primary axis of the array, each bubble benefits from this upwards motion by its 'upstairs' neighbour(s). This effect is stronger at larger Archimedes numbers, corresponding to higher Reynolds numbers and for which the wakes of the bubbles therefore extend further downstream. Cooperative rise is thus due to the strong wake interactions between vertically aligned bubbles. Wake interactions are also visible in the transient evolution of the drift velocity shown in figure 6. This quantity first levels off after an initial transient, then the bubbles experience a significant acceleration at t = O(h/U 0 ), that is as they enter into the wake of their first preceding neighbour. At the smallest volume fraction considered, the time scales separation and the wakes strength are sufficient to distinguish the same phenomenon at t ≈ 2h/U 0 : the bubbles rise is then influenced by the wake of their second preceding neighbour, they accelerate again, and so on until convergence.

To complete this analysis we now evaluate how a simple prediction based on pair interaction compares with our results. We estimate for this the drag coefficient of the trailing bubble of a vertically-aligned pair separated by a distance h (our lattice spacing) and translating with a velocity U identical to that of the array of bubbles at the corresponding volume fraction φ = π/6(d b /h) 3 using the model of [START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF] (equation (6.7) therein). Their expression, which accounts for potential and wake interactions, has been established for Re 20, so we show in figure 7 the results obtained for case S4 and φ 1/3 < 0.55, where this condition is met. At very small volume fraction (φ 1/3 0.13), the drag acting on a bubble of the array is comparable to that exerted on the trailing bubble of pair rising in line, as expected since in dilute conditions, wake interaction between vertically-aligned neighbours dominates. At elevated volume fraction, the dimensionless distance between vertically-aligned bubbles h/d b is smaller, but the drag is no longer governed primarily by such pair interactions, as the departure from the pair-interaction results is seen to be substantial in figure 7.

Deformed bubbles at moderate Reynolds numbers

We now examine the effect of the Bond number on bubble deformation and on hydrodynamic interactions. The effect of volume fraction on U/U 0 is first shown, for different values of the Bond number and comparable Archimedes numbers, in figure 8(a) (Ar ≈ 30) and figure 8(b) (Ar ≈ 15). The data points that are apparently missing at some intermediate volume fractions for case E1 actually correspond to bubbles that do not rise steadily and vertically (discussion of these is postponed to section 4.2 and figure 15(a)), and only small volume fractions are shown for case C because bubbles cannot exist at higher φ (instead, unsteady elongated bodies of gas are obtained). We (a)

q Bo = 5.0 (E2) Bo = 2.0 (E1) Bo = 0.38 (S4) Bo = 0.98 (prior) q q q q q q q q q q q q φ 1 3 U U 0 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.8 0.9 have also included in figure 8(a) the numerical data obtained by Bunner & Tryggvason (2002a), who noticed that normalized drift velocities could be larger than unity at low volume fraction (black crosses). In their study, they estimated Re 0 = 36 by interpolating the data of [START_REF] Ryskin | Numerical solution of free-boundary problems in fluid mechanics. Part 2. Buoyancy-driven motion of a gas bubble through a quiescent liquid[END_REF]. Using the correlation of [START_REF] Loth | Quasi-steady shape and drag of deformable bubbles and drops[END_REF], as we did for our own sets of parameters, we obtained Re 0 = 33. For consistency we kept this value for re-plotting their data. Their results follow approximately the same trend as ours, although their normalized drift velocities are slightly lower. This is probably because the effect of the gas viscosity is assumed to be zero when estimating Re 0 , a hypothesis better approached by our DNS than by that from Bunner & Tryggvason (2002a), who used a gas viscosity twice as large.

1 1.1 1.2 1.3 (b) q Bo = 243 (C) Bo = 0.38 (S3) q q q q q q q q Figure 8.
The shape of U (φ)/U 0 , specifically its non-monotonicity, is similar to that obtained for spherical bubbles, but the faster rise at small volume fraction is more pronounced at higher Bond numbers. The origin of this behaviour becomes clear if one examines the effect of the Bond number on bubbles wakes in figure 5, in which the vertical component of the liquid velocity is represented in a vertical symmetry plane, at small volume fraction (φ = 0.2 %). Each column corresponds to comparable Archimedes numbers (Ar ≈ 15 and Ar ≈ 30 for the second and third columns, respectively) with increasing Bond numbers from top to bottom. As the Bond number increases, the bubbles flatten (as discussed below) and their drag coefficient increases (for a given Re) as a result of the increase of their frontal area (e.g., [START_REF] Loth | Quasi-steady shape and drag of deformable bubbles and drops[END_REF]). This induces greater upward liquid velocities in their wakes and therefore stronger cooperative interactions between in-line objects. The transient evolution of the drift velocity is shown in figure 6. It is similar to that of spherical bubbles, with accelerations at time intervals O(h/U 0 ) but a slightly different initial transient in which the time dependence of acceleration is non-monotonic, a feature related to the bubble deformation from a sphere to an ellipsoid.

These results can be directly compared with those of Sankaranarayanan et al. (2002), who found that the evolution of the drift velocity with the volume fraction follows a Richardson-Zaki relation Re = Re 0 (1 -φ) n [START_REF] Richardson | Sedimentation and fluidisation: Part I[END_REF], where n is given by an empirical closure relation in terms of Re 0 Bo 1/4 /Ar (in the present notation). The trends predicted by this relation strongly disagree with our numerical results, and are (a) therefore not shown. Once again, it appears that this power law dependency on volume fraction may be used to obtain a coarse estimate of the drift velocity at high volume fraction, but does not capture the complex influence of hydrodynamic interactions on the rise of cubic arrays of deformable bubbles for low values of φ.

q C 0 C 1 C m1 q q q q q Bo 10 0 10 1 10 2 -1 -0.5 0 0.5 (b) q q C 0 C 1 C m1 q q q q q Figure 9.
To formulate a semi-empirical law for the function U (φ)/U 0 consistent with the Oseenflow analysis and that would be valid for smaller φ values and very deformed bubbles, we note that the positive root of (4.3) for U can be written in the form

U U 0 = 1 + U 0,Stokes U 0 µ * C 0 -(1.1734 + C 1 )φ 1/3 1 + C m1 φ -1/3 , (4.8)
where U 0,Stokes is given by (4.2). We have introduced in this expression three fitting parameters C 0 , C 1 , and C m1 that we have computed for each case by a least-squares fit of DNS data at finite volume fraction and isolated-bubble data at zero volume fraction.

The fitted values are shown in figure 9. These parameters account for the effect of inertial interactions (they are zero in the Stokes-flow regime, Ar = 0), and are monotonic functions (increasing and decreasing, respectively) of the Archimedes and Bond numbers.

The fitted expression of the normalized drift velocity is shown with dashed lines for each case in figure 3(b), figure 4 and figure 8, which show that the effect of volume fraction is well described by this law for all the flow regimes considered. We now investigate the bubbles' shape. At low Bo (not shown here), the bubbles remain approximately spherical as volume fraction is varied; the aspect ratio χ (the maximum bubble width W divided by the maximum bubble height H, see figure 10) does not deviate from 1 by more than 5 %, due to the low value of the Weber number. At Bo O(1), the bubble shape strongly depends on volume fraction, and is investigated below. For intermediate Reynolds numbers, say 1 < Re < 100, no theoretical expression of χ is available. At low Re and We a theoretical result is available [START_REF] Taylor | On the deformation and drag of a falling viscous drop at low Reynolds number[END_REF] for the shape modes introduced below, but this has been found not to predict accurately results of numerical simulations for an isolated bubble if We is increased to unity [START_REF] Ryskin | Numerical solution of free-boundary problems in fluid mechanics. Part 2. Buoyancy-driven motion of a gas bubble through a quiescent liquid[END_REF]. The analysis by [START_REF] Moore | The rise of a gas bubble in a viscous liquid[END_REF] is for values of Re well over 100. Therefore, as a starting point in the following, the results of numerical simulations of an isolated bubble by [START_REF] Ryskin | Numerical solution of free-boundary problems in fluid mechanics. Part 2. Buoyancy-driven motion of a gas bubble through a quiescent liquid[END_REF] are used.

We first focus on ellipsoidal bubbles (cases E1 and E2) †. In all the simulations reported here, bubbles are virtually axisymmetric, but may exhibit significant fore-andaft asymmetry. In figure 10, we present the aspect ratio as a function of volume fraction for these cases. Also shown are the first two shape coefficients a 2,3 defined by writing the local bubble radius R(θ) as

R(θ) = n a n P n (cosθ),
(4.9)

where P n is the Legendre polynomial of order n and θ is the angle between the position vector at the bubble surface and the bubble velocity (a 0 is the radius of the sphere with the same volume). The coefficients were obtained from a distribution of points (at least 500) on the bubble surface and integrating the orthogonality relation for each Legendre polynomial.

In figure 10(a-d), for both cases E1 and E2, as the volume fraction is increased, χ decreases monotonically to unity, and a 2 goes to zero, from about their respective values for isolated bubbles. The corresponding single bubbles (for which results are shown in the figure with filled symbols, with shape coefficients obtained from [START_REF] Ryskin | Numerical solution of free-boundary problems in fluid mechanics. Part 2. Buoyancy-driven motion of a gas bubble through a quiescent liquid[END_REF]), are of oblate-ellipsoid shape. This shape is expected for isolated bubbles at large Re, through a Bernoulli suction effect in the vicinity of the bubble rim, and is expected also at low Re [START_REF] Taylor | On the deformation and drag of a falling viscous drop at low Reynolds number[END_REF]. It may be anticipated that the demise of this shape at elevated volume fraction is partly due to the dependencies of Re and We on volume fraction. We have verified, however, that the empirical correlation by [START_REF] Loth | Quasi-steady shape and drag of deformable bubbles and drops[END_REF] for the aspect ratio of isolated bubbles, using the values of these dimensionless groups obtained from the simulations, although yielding good agreement at φ = 0, gives a very poor prediction of the results presented in figure 10 (and is therefore not shown). The reduction of suction at the bubble rim is therefore due to the detailed bubble interactions. In their study of the hydrodynamic interactions between two spherical bubbles rising side by side, [START_REF] Legendre | Hydrodynamic interactions between two spherical bubbles rising side by side in a viscous liquid[END_REF] showed that at small to moderate Reynolds number (Re 30, as encountered in our study), the transverse force is repulsive and increases when the separation decreases. Such a reduction or elimination of suction between bubbles suggests the liquid downflow due to a bubble pair occurs around the pair as a whole. In a 3D cubic array, although the room for liquid to flow down with little opposition is reduced further (at four sides along a bubble rim), some remains present along vertical edges of each cell in the array. Therefore, any suction effect normally arising at the rim of a bubble would be reduced in between bubbles lying in the same horizontal plane. This may be somewhat countered by an increase elsewhere along the bubble periphery (if not at a greater distance from the bubble), but a variation in curvature is opposed by surface tension.

In addition to this reduction in aspect ratio, the fore-and-aft asymmetry of the shape of an isolated bubble is altered significantly by the finite volume fraction. For single bubbles, cases E1 and E2 are near the boundary between a low-Re regime of bubbles with a blunt tail and, at the same We but larger Re, a regime of bubbles with a flat nose [START_REF] Ryskin | Numerical solution of free-boundary problems in fluid mechanics. Part 2. Buoyancy-driven motion of a gas bubble through a quiescent liquid[END_REF]; only a slightly flattened nose is observed, mostly in case E2, resulting in a positive value of a 3 in (4.9). The results in figure 10(b,d) show that already at small but finite volume fraction, this asymmetry is reversed. The bubble nose becomes rather pointed and the tail blunt as the volume fraction is increased further. This tendency for oblate ellipsoidal bubbles arranged in regular arrays to have their nose pulled upwards at finite shape coefficients a n q a 2 a 0 a 3 a 0 q q q q q q q q q q φ 1 3

χ = W H 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 (c) H W φ 1 3 0 0.2 0.4 0.6 0.8 -0.4 -0.3 -0.2 -0.1 0 0.1 (d)
shape coefficients a n q a 2 a 0 a 3 a 0 q q q q q q q q q q q q q φ 1 3 χ = W H volume fraction has been observed previously by Sankaranarayanan et al. (2002), and attributed to a wake effect. Indeed [START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF] showed that in the present range of Reynolds numbers, two spherical bubbles rising in line are attracted toward each other for separation distances greater than approximately 1.3 bubble diameter, which would be equivalent to φ 1/3 = 0.62, a value close to the upper bound of the range of volume fractions we consider. To investigate this further, the amplitude of the P 3 mode in (4.9) is included in figure 10. It is seen that in concentrated arrays this becomes as significant as that of the P 2 mode. The presence of successive bubbles in each others' wakes does reduce the variation in velocity magnitude between them, as can be seen in figure 11 (top row). The significance of the stagnation-point flow at the bubble nose is thereby reduced as volume fraction is increased, and the large dynamic pressure at a stagnation point in the liquid is reduced (the pressure field is shown in figure 11, bottom row), along with the magnitude of normal deviatoric stress. Both these result in an increase in the jump in normal stress and hence an increase in interface curvature. Finally we note that the results for aspect ratio and shape coefficients for cases E1 and E2 differ by a factor of approximately two, which roughly corresponds to the ratio of the bubble Weber numbers at all volume fractions. As we have not undertaken to extend our parametric study even further to confirm, it is concluded that the results for (χ -1) and a 2,3 versus Weber number, in the present range of We < 7, are consistent with a linear dependency.

We have also investigated a regime characterized by a very high Bond number (case C, Bo = 243). The evolution of the bubbles' "aspect ratio" and shape with volume fraction is shown in figure 10(e). The corresponding bubbles in isolation are indented ellipsoidal caps [START_REF] Bhaga | Bubbles in viscous liquids: shapes, wakes and velocities[END_REF][START_REF] Hua | Numerical simulation of 3D bubbles rising in viscous liquids using a front tracking method[END_REF]. As the volume fraction increases, the upside-down crown of gas issuing from the bubble rim becomes thinner and longer to form a skirt with an inward curvature (in the direction of the bubble axis of symmetry).

The range of φ values that can be considered in this case is relatively narrow, since for φ 6 % the bubbles become so elongated that they coalesce. The theory of [START_REF] Ray | On skirted drops in an immiscible liquid[END_REF] indicates that the finite length of the skirt is dictated by the thinning of the skirt downstream of its point of formation. According to their model, the skirt thickness is proportional to √ -u s , where u s is the (negative) vertical component of the relative velocity (that is, the liquid velocity in the bubble's frame of reference) at the outer side of the skirt at a given altitude (the inward curvature of the skirt being neglected). We show in figure 12 the vertical relative velocity u rel 3 = u 3 -u 3 d in the vicinity of bubble (with gravity pointing in the -e 3 direction). It can be observed that as the distance from the rim increases, |u s | ( u rel 3 along the outer side of the skirt) decreases and the skirt tapers, until the skirt ends for a critical value of |u s |, in (qualitative) agreement with the model of [START_REF] Ray | On skirted drops in an immiscible liquid[END_REF]. By comparing figure 12(a) and (b) one remarks that at high volume fractions a significant downflow of liquid develops outside of bubbles wakes. This backflow of liquid, which is particularly strong because the bubbles rise velocity is not substantially reduced for large values of φ, increases the value of |u s | at a given altitude. Therefore, at higher volume fraction, the skirt must extend further downstream to reach the critical value of |u s | at which the skirt ends, as observed in our simulations.

The aspect ratio of the skirted bubble decreases towards unity as the volume fraction is increased, even if the skirt is not included in the height. Indeed given the moderate value of the Reynolds number (Re ≈ 10), the same reasoning as above for ellipsoidal bubbles is expected to apply, that is, a decrease of the suction effect as the size of the gap between side neighbours decreases. At low volume fraction, an extrapolation of the results to zero volume fraction appears consistent with the corresponding result for a single bubble.

Steady and unsteady oblique rise of bubbles

In the range of parameters considered thus far, the motion of a single bubble in unbounded liquid is straight, steady, and parallel to gravity. Bubble motion that is oblique (i.e., not aligned with gravity) was observed, however, for cases E1 and S5 at certain volume fractions. Such oblique motions have been reported previously for twodimensional square arrays of bubbles (Sankaranarayanan et al. 2002;Sankaranarayanan & Sundaresan 2002;[START_REF] Theodoropoulos | Coarse bifurcation studies of bubble flow lattice Boltzmann simulations[END_REF]), but their triggering and their stability remain essentially unexplained. We analyse oblique rise further here. In this subsection the bubble drift velocity vector is denoted by U = U i e i and the gravity vector by g = -ge 3 , where e i are the unit primary vectors of the periodic array. First, existence of such oblique solutions is demonstrated at Reynolds numbers that are small but finite. The Oseen analysis in appendix B yields the system of equations (B 3) that involve the bubble velocity and the force exerted by the fluid on the bubble. As this force is prescribed (it balances buoyancy), (B 3) yields the bubble velocity. The main solution is naturally a velocity vector aligned with gravity, as studied in section 4.1. Equation (B 3) does, however, allow for other solutions which satisfy the nonlinear system of equations (B 12). We have found these non-trivial solutions at values of Ar around 20. The most convenient way to obtain these solutions was found to be, for a given inclination of the bubble velocity (with respect to the upward vertical direction), to reduce the problem to a single nonlinear equation for the Reynolds number based on the lattice spacing and the magnitude of the bubble velocity (Re h in the Appendix), and to obtain the volume fraction from the remainder of the system of equations. Two types of non-trivial solutions were studied: either the horizontal bubble velocity component was aligned with one of the lattice unit vectors, or it was diagonal to the lattice. The results are presented in figure 13. It is seen that these exist below a critical value of the volume fraction (which we have found to increase rapidly with the value of Ar ), the inclination angle strongly increasing as the volume fraction is reduced. Oblique rise of single light solid particles has already been observed in numerical simulations [START_REF] Jenny | Instabilities and transition of a sphere falling or ascending freely in a Newtonian fluid[END_REF]) and experiments [START_REF] Veldhuis | An experimental study of the regimes of motion of spheres falling or ascending freely in a Newtonian fluid[END_REF], also arising as a non-trivial further solution [START_REF] Fabre | The steady oblique path of buoyancy-driven disks and spheres[END_REF]), but the present solutions crucially involve bubble interactions, we investigate this further below.

We now return to numerical results. After an initial transient during which they accelerate from rest under the effect of buoyancy, the bubbles may be deflected from their original vertical trajectories. At this point, the horizontal components of the bubble velocity grow in magnitude while the rise velocity drops off. After that, velocity fluctuations set in. Three types of dynamic behaviours have been identified depending on the evolution of these fluctuations: (i) the fluctuations may rapidly dampen out, and the bubbles finally rise steadily on a straight (but skewed to gravity) path; (ii) they may take the form of oscillations, so that the bubbles motion is oscillatory around a straight oblique path; (iii) they may become aperiodic, so the bubble rise is chaotic and, on average, not aligned with gravity. These regimes are exemplified in figure 14, and will be respectively referred to as steady oblique rise (a,b), oscillatory oblique rise (c), and chaotic oblique rise (d).

The bubble drift velocities are in all cases either steady or statistically stationary, so mean drift velocities can be defined by averaging over a sufficient time period. The horizontal and vertical components of the (statistically-)steady drift velocity are plotted against volume fraction in figure 15 for cases E1 and S5. Filled symbols are used for steady bubbles motion, vertical or oblique. For unsteady bubble motion, the time-averaged drift velocity is shown with open symbols, and the standard deviation is represented using vertical bars. In case E1 (figure 15(a)), steady vertical rise is obtained at low and high volume fractions, whereas at intermediate volume fractions the bubbles rise steadily along an oblique path with an inclination angle of about 3 • . This figure shows that three solutions exist in this regime: a symmetric vertical solution (U 1 = U 2 = 0), and two asymmetric oblique solutions consisting of horizontal velocity components of equal q Re 1 × 10 Re 2 × 10 Re 3 q q q q q q q q q q q q q φ 1 3 Re 1 Re 2 Re 3

Re 1 Re 2 Re 3 q q q q q q q q q q q Figure 15. Quasi-steady drift velocity components (U = Uiei) in the form of Reynolds numbers (Rei = ρcUid b /µc) as a function of volume fraction for an ordered array of bubbles in the cases where bubbles exhibit non-vertical motions. Gravity is pointing in the -e3 direction. (a) Case E1 (Bo = 2.0, Ar = 29.9, note that the horizontal components have been multiplied by 10 for clarity). (b) Case S5 (Bo = 0.38, Ar = 40.7). Filled symbols are used when the bubbles motion is steady (oblique or vertical). Open symbols and vertical bars are used when the bubbles motion is unsteady: symbols indicate the mean drift velocity, and bars show its root mean square. The dashed line is a numerical fit of the form of (4.8) matching DNS and isolated-bubble data for the vertical rise only.

magnitude (|U

1 | = |U 2 | = 0),
as predicted from the Oseen-flow analysis (figure 13). In case S5 (figure 15(b)), steady vertical rise, steady oblique rise, oscillatory oblique rise, and chaotic oblique rise are obtained in that order as volume fraction is decreased. Inclination of the velocity with respect to the upward vertical direction is between 6 and 13 • (maximum for φ 1/3 = 0.4). As volume fraction approaches zero, the steady vertical rise of the isolated bubble must be recovered, although the occurrence of this transition cannot be evidenced by numerical simulations owing to their prohibitive cost.

It is possible to obtain insight in this behaviour by using prior results for bubble pairs. At steady-state, the integral of fluid stresses over the bubble surface, denoted by f , is balanced by the buoyancy force f buoy : f = -f buoy , with f buoy = f buoy e 3 = 1 6 πd 3 b (ρ c -ρ d )ge 3 . The total surface force f acting on the bubble can be decomposed into a drag force f drag and a lift force f lift , defined by

f drag = (f • U ) U |U | 2 , f lift = f -f drag , (4.10)
and corresponding to longitudinal and transverse components of f with respect to the direction of motion, respectively (these definitions can be used for unsteady but statistically stationary systems by replacing U and f by their time averages). The persistence of a (possibly average) oblique motion implies the existence of a net (average) lift force exerted on the bubble. The magnitude of this lift force is classically presented in the form of a dimensionless lift coefficient C L defined by

C L = |f lift | 0.125πd 2 b ρ c |U | 2 . (4.11)
The (average) lift coefficient is plotted as a function of volume fraction in figure 16 for case E1 (open grey triangles) and case S5 (open black squares).

Since the tilt angle remains small (not larger than 13 • in our simulations), each bubble is in the wake of its predecessor, and the oblique path is expected to originate from Filled symbols: for the trailing bubble of a pair of spherical bubbles within the same conditions (orientation, Reynolds number, separation distance) as two vertically-aligned bubbles of the ordered array, with (large symbols) and without (small symbols) accounting for the interaction with the wake of the leading bubble, from [START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF].

the vorticity produced by the preceding bubble. We therefore investigate whether the lift force induced by the preceding bubble can be estimated from prior work on bubble pairs separated by a fixed distance equal to the present lattice spacing, both rising at a constant velocity U and where the angle between U and the vertical line joining their centers is the inclination angle measured from our simulations. The model proposed by [START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF] for bubble pairs (equation (6.12) therein) is used here for this purpose. A spherical bubble shape appears a reasonable approximation given that in our simulations the aspect ratio does not exceed 1.3 for case E1, and 1.1 for case S5.

The results are shown in figure 16 with larger filled symbols. To assess the influence of the leading bubble's wake, the lift coefficient obtained for bubble pairs without the contribution from the wake is shown on the same figure with smaller filled symbols. The trend obtained by considering the pair interaction is in excellent adequacy with our numerical results for periodic arrays when the wake of the top bubble is accounted for, thereby demonstrating that oblique rise is essentially a wake-induced effect. The lift coefficient in ordered suspensions is found to be larger than that due to the interaction with the wake of a single bubble, and the difference is more pronounced at higher volume fractions, since in the periodic configuration the bubble has an infinite number of top neighbours that may contribute to the lift force. This reasoning can even be made more precise by considering the expression of the lift force acting on a single spherical bubble moving in a (e.g. wake-induced) rotational flow [START_REF] Auton | The lift force on a spherical body in a rotational flow[END_REF][START_REF] Legendre | The lift force on a spherical bubble in a viscous linear shear flow[END_REF][START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF]:

f lift ∝ d 3 b ρ c Ω × U . (4.12)
In this expression, Ω = |Ω|e h is the liquid vorticity "seen" by the bubble and produced by the motion of all the other bubbles. No clear definition of this quantity is available if vorticity is not uniform at the bubble scale, as is the case in the present study, but it is reasonable to assume that it can be qualitatively estimated by examining vorticity profiles in the bubble vicinity. In turbulent flows, [START_REF] Merle | Forces on a high-Reynolds-number spherical bubble in a turbulent flow[END_REF] and Naso & Prosperetti (2010) approximated the velocity and vorticity "seen" by bubbles and solid particles respectively, by the average of these quantities over shells of different sizes. In order to show that our results are qualitatively consistent with (4.12), it is thus convenient to introduce an orthonormal direct basis (e , e ⊥ , e h ), defined by the unit vectors

e = U |U | , e ⊥ = f lift |f lift | , e h = e × e ⊥ . (4.13)
The bubble's steady motion is contained in the vertical plane defined by (e , e ⊥ ). We now examine the sign and magnitude of the liquid vorticity component ω h = ω • e h ahead of the bubble, where it should overall give a positive contribution to |Ω| for the above model to be correct. We show in figure 17 the liquid vorticity field projected onto e h in the vertical plane normal to e h for two examples of steady oblique rise. The bubble in the right panel experiences a stronger lift force (indicated by the thick orange arrow) than the bubble in the left panel. This is consistent with (4.12) and the fact that, upstream of the bubble (its drift velocity being shown with the black arrow), ω h is positive, and its magnitude is larger than that in the left panel (although only two examples are shown here the same result holds for all our simulations). In addition it is seen from figure 17 that ω h is transported from the surface of preceding bubbles, hence confirming the key role played by the wakes and the associated lift force for the stability of oblique motion. We will now examine the time dependence of bubbles motion in the unsteady regimes. As illustrated in figure 18, which depicts the single-sided amplitude spectra of the discrete Fourier transforms of the drift velocity components time signals in the oscillatory (left) and chaotic (right) rise regimes, a spectral analysis of the unsteady velocity signals reveals clear peaks at a frequency equal to f h = (U 1 + U 2 )/(2h) (where the bar denotes a timeaverage). Normalizing frequencies by f 3 ≡ U 3 /h (not shown here) does not lead to a collapse of the curves. This suggests that the force fluctuations experienced by a bubble are also driven by the interaction with the wakes of the preceding bubbles that are not on the same vertical axis. As a consequence, the dynamic behaviour of a bubble in an ordered array, although greatly influenced by the direct interaction with its top neighbour, is also dictated by longer-range nonlinear interactions with other bubbles located in above horizontal planes.

In the light of these results, we are now in a position to propose the following scenario for explaining the transitions between the various regimes of motion reported in figure 15.

frequency f h amplitude U 3 0 1 2 3 4 5 0 1 2 3 4 5 (a) q q q U 1 U 2 U 3
× 10 -2 q qqqqqqqqqqq q q q qq q q q q q qq q q qqqqq qqq qqqq q q qq qqq qqqqqqqqqqqqqq q qqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqq q qqqq q q q qq qq qqqqqqqqqqqqq q qqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q qqqqqqqqqqqqq qq q q q q qq qq qq qqqqqqqqqq qqq q qq qqq qqqqqqqqqqqqqq q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq frequency f h × 10 -2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 18. Single-sided amplitude spectra of discrete Fourier transforms of the drift velocity components time signals. Unsteady oblique rise of ordered arrays of bubbles: (a) oscillatory motion (case S5 at φ = 3.8 %), (b) chaotic motion (case S5 at φ = 1.6 %). The amplitude is normalized by the average vertical drift velocity U 3, the frequency is normalized by

f h = (U 1 + U 2)/(2h).
First, non-vertical motion can only occur when the flow conditions allow sufficiently high Reynolds numbers to be attained (here, cases E1 and S5). Vorticity then becomes significant in the vicinity of each bubble due to the wake of its predecessor; an infinitesimal asymmetry can then result in a lift force that is sufficient to result in oblique motion [START_REF] Koch | Hydrodynamic diffusion in dilute sedimenting suspensions at moderate Reynolds numbers[END_REF]. If each bubble is only influenced by the wake of its immediate predecessor, this motion is steady. When the wakes extend horizontally over distances large compared to the lattice spacing, each bubble interacts with the wakes of a great number of neighbours, including some that are not located on the same vertical axis, and the motion becomes chaotic. Then, for a given flow regime, the volume fraction is in the first place related to the distance between the bubbles, but also affects the Reynolds number in a nonmonotonic manner. At low volume fraction, when the Reynolds number increases with volume fraction, steady vertical rise, steady oblique rise, and unsteady oblique rise occur in that order at increasing volume fraction. At higher volume fraction, the situation becomes more complex because the Reynolds number decreases with φ. It appears that the dominant effect of increasing volume fraction is then not to bring the bubbles closer to each other, but to reduce their velocity, so that steady oblique rise is first recovered, and is replaced by steady vertical rise at the highest volume fractions.

Free arrays

We examine in this section the behaviour of freely evolving bubbly flows as represented by the repetition of a unit cell containing several independent bubbles ("free arrays"). This problem, studied previously by several groups (e.g., Bunner & Tryggvason (2002a); [START_REF] Esmaeeli | A direct numerical simulation study of the buoyant rise of bubbles at O(100) Reynolds number[END_REF]; [START_REF] Yin | Lattice-Boltzmann simulation of finite Reynolds number buoyancydriven bubbly flows in periodic and wall-bounded domains[END_REF]) but only at moderate and high volume fractions, is revisited here following the insights gained in the previous section for ordered arrays. Our main objective in this section is to investigate the dynamics of free arrays at small and intermediate volume fractions, and to compare the observed trends with those obtained for ordered suspensions.

Simulations of free arrays of bubbles have been undertaken for cases E1 and C (table 1). The two fluids were initially at rest, N b identical spherical bubbles were introduced in a cubic periodic unit cell of size h (which results in a gas volume fraction φ = (πN b d 3 b )/(6h 3 )), and gravity was switched on at time zero. It was found that the transient evolution of the system can follow either of two routes: one with successive pair coalescence events until an ordered array configuration is recovered, and the other one in which the number of bubbles remains constant throughout the simulation. The former was observed for case C, which corresponds to highly deformable bubbles of dimpled or skirted ellipsoidal-cap shapes, even at relatively low volume fractions (the lowest volume fraction considered in that case was φ = 0.8 %). In case E1 (weakly ellipsoidal bubbles), by contrast, coalescence was never observed, provided that the bubbles interfaces are initially sufficiently separated from each other and that the volume fraction remains below approximately 5 % (we shall elaborate on this last point in section 5.2). Examination of the suspension evolution revealed that bubbles never come into close contact in that case, as previously observed by [START_REF] Esmaeeli | Direct numerical simulations of bubbly flows. Part 2. Moderate Reynolds number arrays[END_REF] in a similar flow regime. After a transient regime, the flow was found to become independent of the initial bubbles positions and a well-defined statistically steady state was reached.

In the following we analyse the statistically stationary rise of free, non-coalescing, deformable bubbles at moderate Reynolds number (case E1). About fifty simulations of free arrays were run in total, corresponding to different initial conditions, numbers of bubbles, and volume fractions. For each of these, the transient evolution of the system was monitored through the time signals of the bubble drift velocity U (defined, as for ordered arrays, by equation (2.3), and therefore equal to the average drift velocity of the N b bubbles) and of the interface surface area A (which is a measure of the average deformation of the N b bubbles). In all simulations, the instantaneous horizontal components of U were found to be negligibly small, so we shall hereinafter simply use U to denote the vertical component of the drift velocity vector. Each simulation was continued until U and A became statistically stationary. Their time averages, denoted by overbars in what follows, were then computed by averaging over a sufficient time interval.

Convergence with the number of bubbles

The influence of the number of free bubbles in the unit cell is evaluated by varying N b from 2 to 27 while keeping the volume fraction (and all other parameters) constant. The evolution of the bubbles drift velocity with the number of bubbles is shown in figure 19 for φ = 2.4 % (filled circles). The main effect of additional degrees of freedom is to slow down the bubbles: the drift velocity drops by 15 % when the relative motion between two bubbles is allowed, and is reduced further (up to ≈ 30 %) if the number of bubbles in the unit cell is increased. For N b O(10), the drift velocity becomes nearly independent of the number of bubbles. The rate of convergence and the maximal relative decrease in drift velocity with the number of bubbles appear to be essentially independent of the volume fraction, at least in the limited range considered here, as shown in figure 19 (open triangles and squares). The drift velocities obtained by Bunner & Tryggvason (2002a) for a similar flow regime have been reported in the same figure (crosses). It is worth mentioning that although the maximum number of bubbles shown in figure 19 is Bunner & Tryggvason (2002a) have performed simulations for 1 N b 216 (see Figure 8a in their paper), and also concluded that the effect of the system size on the drift velocity becomes negligible for N b O(10) in this flow regime. Overall the agreement between the two data sets is excellent, including in the peculiar case N b = 4.

N b = 27,
It can indeed be noticed in figure 19 that convergence is not monotonic, and that the rise is abnormally slow for N b = 4. Visual inspection of the spatial distribution of the bubbles reveals a significant preference for horizontal alignment in that case. This bias is particularly pronounced for values of φ that are not very small and is therefore 0.8 0.9 1 q φ = 1.6 % (present) φ = 2.4 % φ = 3.8 % φ = 6.0 % (prior) q q q q q q q q q q q Figure 19. Influence of the number of bubbles on the average bubbles drift velocity for various volume fractions. The drift velocity is normalized by that obtained in the ordered configuration (N b = 1). Symbols other than crosses: present DNS for case E1 (Ar = 29.9, Bo = 2.0). Crosses: prior DNS of Bunner & Tryggvason (2002a) for a comparable flow regime (Ar = 29.7, Bo = 0.98). illustrated in figure 20 for φ = 3.8 %: while the bubbles are rather uniformly distributed within the periodic cell for N b = 8 (right), as is the case for other typical values of N b , they all remain in the same horizontal plane when N b = 4 (left), so that in the latter case the suspension actually consists of successive horizontal layers of bubbles. As shown by [START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF], side-by-side alignment maximizes the drag force acting on a pair of bubbles, resulting in lower drift velocities than with other types of spatial distributions. This particular behaviour for N b = 4 demonstrates that an ordered microstructure is not always unstable. Such arrangements in horizontal planes are indeed possible if the number of free bubbles possesses an integer square root, due to periodicity and system symmetries. For N b = 9, the bubbles also tend to arrange within a single horizontal plane, but this arrangement rapidly breaks up and is only observed intermittently. This results in the small but noticeable anomalous reduction of the drift velocity visible for N b = 9 in figure 19. For N b = 16, no horizontal layer of bubbles is formed during the simulation, and no anomaly is detectable in figure 19. We conclude that the artificial effects of symmetry and periodicity observed when N b is an exact square rapidly vanish as the number N b of independent bubbles is increased. ordered, Bo = 2.0 free, Bo = 2.0 free, Bo = 1.0 (prior) free, Bo = 4.9 (prior) q q q q q q q q q q q q q q q Figure 21. Influence of volume fraction on the bubbles drift velocity in ordered and freely evolving suspensions, with comparable Ar (≈ 30) and various Bo. The drift velocity is normalized by the terminal velocity of an isolated bubble for the same Archimedes and Bond numbers (estimated from [START_REF] Loth | Quasi-steady shape and drag of deformable bubbles and drops[END_REF]).

• and ----: present DNS data, and their numerical fit, for N b = 1 (ordered arrays) and Bo = 2.0; •: present DNS for free arrays of 8 bubbles with Bo = 2.0; : prior DNS of Bunner & Tryggvason (2002a) for free arrays of 27 bubbles with Bo = 0.98; : prior DNS of [START_REF] Bunner | Effect of bubble deformation on the properties of bubbly flows[END_REF] for free arrays of 27 bubbles with Bo = 4.9.

The effect of volume fraction on bubble drift velocity and deformation

The volume fraction has been varied from 0.2 to 3.8 % by reducing the size of the unit cell while keeping all the other parameters constant, for a number of freely moving bubbles set to N b = 8. For volume fractions greater than 5 %, numerical coalescence (that is, coalescence due to the spacing between bubble interfaces being less than the grid spacing) occurs during the transient evolution of the flow, therefore no data could be obtained for high volume fractions.

Data for fairly high volume fractions are however available from prior studies of Bunner & Tryggvason (2002a) and [START_REF] Bunner | Effect of bubble deformation on the properties of bubbly flows[END_REF] who performed simulations of free arrays of bubbles using a front-tracking method that does not allow coalescence. In their simulations, the Archimedes numbers are similar to ours (Ar ≈ 30) but the Bond numbers (and hence the bubble shapes) are different: the bubbles are nearly spherical (Bo = 0.98) in Bunner & Tryggvason (2002a), they are oblate ellipsoids (Bo = 5.0) in [START_REF] Bunner | Effect of bubble deformation on the properties of bubbly flows[END_REF], our present simulations are for Bo = 2.0. In order to present their results together with ours, we first estimate the terminal velocity of the corresponding isolated bubbles. We proceed for that purpose in the same manner as we did for our own simulations, that is by using the correlation of [START_REF] Loth | Quasi-steady shape and drag of deformable bubbles and drops[END_REF] for single ellipsoidal bubbles, which leads to Re 0 = 33 for Bo = 0.98, and Re 0 = 26 for Bo = 5.0.

The influence of volume fraction on the drift velocity normalized by the terminal velocity of the same bubble in unbounded liquid is shown in figure 21. Squares and triangles correspond to prior simulations of nearly spherical and oblate ellipsoidal bubbles, respectively. Filled circles correspond to our present simulations of weakly ellipsoidal bubbles. For comparison, the results we obtained for the corresponding ordered arrays are shown in the same figure with open circles. Remarkably, the evolution of the drift velocity with φ seems to be different in dilute and fairly dense suspensions of free bubbles, as is the case when bubbles are perfectly ordered. We have checked that neither a linear evolution with φ nor a law of the form U /U 0 = k(1 -φ) n (with k and n free parameters) is compatible with the data presented in figure 21.

At moderate to fairly high volume fractions (say, 0.015 φ 0.25, that is, 0.25 ordered free sphericity = A sphere A bubble q q q q q q q q q q q q q Figure 22. Influence of volume fraction on the bubbles sphericity in ordered (open circles) and free (filled circles) arrays for case E1. The sphericity is defined as the ratio between the total surface areas of a set of N b volume-equivalent spheres and that of the bubbles.

φ 1 3 0 0.
φ 1/3 0.63), the drift velocity of free bubbles decreases approximately linearly in φ 1/3 . This scaling, which is also independent of the bubble oblateness (in the limit of the range of shapes considered here), agrees with that obtained for ordered suspensions in the same conditions.

The drift velocity dependence on volume fraction is radically different at vanishing φ: although we cannot approach the dilute limit in the simulations, it is clear that a simple extrapolation to φ = 0 from results at larger φ is not feasible. In the absence of inertial effects and in this dilute limit, a linear reduction of the drift velocity with φ would be expected, according to the analytical solution from [START_REF] Keh | Slow motion of multiple droplets in arbitrary threedimensional configurations[END_REF], derived for random bubbly suspensions in the Stokes flow regime. At finite Reynolds numbers, however, inertial effects are expected to dominate far from the bubbles. The results for free arrays at low φ in figure 21 suggest strongly that bubbles in dilute suspensions rise faster than their isolated counterpart, as in ordered arrays, due to cooperative wake interactions. Such interactions would be much weaker than in ordered suspensions due to the less likely occurrence of vertical alignments (since a spherical or slightly oblate bubble lying in the wake of one of its neighbour experiences a transverse lift force directed away from the wake so that two weakly ellipsoidal bubbles cannot remain in line), but they might still play a role in the suspension dynamics. The uncertainty of the terminal velocity of isolated bubbles prevents us from drawing definitive conclusions on this point, as that would require simulations at even lower volume fractions, beyond the reach of the computational capabilities at our disposal.

The effect of volume fraction on bubble deformation is now evaluated, both qualitatively from visualizations of the flow, and quantitatively through the measurement of their interfacial surface area, larger surface areas being associated with a stronger departure from the spherical shape. We show in figure 22 the volume fraction dependence of the bubbles sphericity, defined as the ratio between the total surface areas of a set of N b volume-equivalent spheres and that of the bubbles. The trends obtained for ordered and free arrays are qualitatively similar, ellipsoidal bubbles becoming more spherical as volume fraction increases. A plausible explanation of the larger oblateness (smaller sphericity) of freely moving bubbles is the weaker role of wake-induced nose elongation, due to the less likely occurrence of vertically-aligned pairs, as explained above.

An explanation of the observed similarities between freely evolving suspensions and ordered arrays at small to intermediate volume fractions may be sought in the fact that in the former, the bubbles spatial distribution is nonrandom and possesses a certain degree of order. The presence of order in suspensions is classically evaluated using the structure factor or pair distribution function. Due to the level-set method, tracking individual bubbles is not automatically done and would require a significant further coding effort. For this reason the microstructure has not been evaluated quantitatively in our simulations. Nevertheless flow visualizations have been used for a qualitative evaluation of the bubble spatial distribution.

Visualizations of bubble motion reveal that free bubbles rise at comparable velocities with very weak horizontal displacements and never get close to each other for the entire range of volume fractions we considered. This may be seen in supplementary movies 1 and 2, which show top views of the unsteady bubble motion at the maximum (φ = 3.8 %) and minimum (φ = 0.24 %) considered volume fractions, respectively. Their spatial distribution within the cell is fairly uniform, and their relative positions remain more or less constant as they rise, in agreement with prior observations by [START_REF] Esmaeeli | Direct numerical simulations of bubbly flows. Part 2. Moderate Reynolds number arrays[END_REF], who also found, for similar flow conditions, that bubbles dispersion in the horizontal direction is almost absent.

These observations are consistent with prior quantitative evaluations of the microstructure of dilute and moderately concentrated suspensions of (nearly) spherical bubbles at Re = O(10). On the experimental side, [START_REF] Cartellier | Bubble-induced agitation and microstructure in uniform bubbly flows at small to moderate particle Reynolds numbers[END_REF] evidenced that in the range 10 -4 < φ < 10 -2 , a test bubble experiences a deficit of neighbours in its immediate vicinity and an excess of neighbours at the border of the deficit zone, or in other words, that a certain degree of order is present in the suspension. The magnitude and extent of the deficit zone decrease with increasing φ, but a clear nonrandom microstructure has been shown by [START_REF] Cartellier | Induced agitation in homogeneous bubbly flows at moderate particle Reynolds number[END_REF] to persist at least up to φ = 0.08. On the numerical side, Bunner & Tryggvason (2002a) and [START_REF] Yin | Lattice-Boltzmann simulation of finite Reynolds number buoyancydriven bubbly flows in periodic and wall-bounded domains[END_REF] identified analogous deficits of bubbles at short distances and excesses of bubbles farther away from a test bubble for 0.02 φ 0.12 and φ = 0.05 in their respective simulations.

Although bubbles are free to sample the entire liquid, they stay with the same neighbours for long times, which explains why suspensions of free bubbles share some properties with perfectly ordered ones, at least up to moderately high volume fraction and at moderate Reynolds number.

A related point, beyond the scope of the present study, is the properties of the bubbleinduced liquid agitation. Information about this can be found in [START_REF] Loisy | Direct numerical simulation of bubbly flows: coupling with scalar transport and turbulence[END_REF].

Comparison with experiments

Finally, we investigate the relation of our results and prior experimental data. A direct comparison between direct numerical simulations and experiments is often impossible because the typical flow conditions differ strongly between these two approaches (moderate vs. high Reynolds numbers, nearly spherical vs. wobbling bubbles, absence vs. presence of surfactants, monodispersity vs. polydispersity, constant vs. varying bubble diameter at varying φ, etc.). To the best of our knowledge, the only experiments carried out under conditions comparable to those in the present work are those of [START_REF] Martinez-Mercado | Measurement of pseudoturbulence intensity in monodispersed bubbly liquids for 10 < Re < 500[END_REF], who measured the average velocity of nearly monodisperse air bubbles rising in a mixture of water and glycerin (50 % mass fraction), for volume fractions ranging from 0.4 to 6.5 %. Importantly, they found the bubble equivalent diameter to be almost independent of the gas volume fraction, so that comparison with our numerical data is relevant. According to the physical properties of the fluids and bubble equivalent diameter (d b = 1.20 ± 0.05 mm) reported in their paper, their experimental conditions correspond to Ar = 26.3 ± 1.6 and Bo = 0.25 ± 0.02. In this regime, bubbles are nearly spherical (as confirmed by the photographs in their paper), so the terminal Reynolds number of the q present simulations (ellipsoidal bubbles) prior experiments (spherical bubbles) q q q q q q equivalent isolated bubble can be estimated from the correlation of [START_REF] Mei | A note on the history force on a spherical bubble at finite Reynolds number[END_REF] (as we did for our simulations of spherical bubbles), which yields Re 0 = 29.7 ± 3.1. These experimental conditions are therefore comparable to our case E1 (Ar = 29.9, Bo = 2.0, Re 0 = 31) in terms of Archimedes and Reynolds numbers, but not in terms of Bond numbers.

The evolution of the bubble drift velocity with volume fraction is shown in figure 23, wherein experimental measurements are represented by crosses and the present numerical results by circles. Numerical and experimental trends are very similar, both exhibiting two different scaling laws at moderate and low volume fractions. In particular, experimental data are compatible, like the numerical ones, with a linear dependence of the rise velocity on φ 1/3 in the case of moderately concentrated suspensions. This behaviour suggests that ordered arrays are able to capture some properties of real bubbly suspensions. Besides, we note that experimental velocities are systematically lower than that predicted from our simulations. Although perfect agreement is not expected due to the differences in the flow conditions, as discussed below, we speculate this may also be partly due to our work being on perfectly homogeneous suspensions, whereas experiments may be affected by the presence of walls and of weak gradients.

The main difference between our simulations and the above-mentioned experiments is the value of the Bond number which directly alters the bubble shape (ellipsoidal in the former, nearly spherical in the latter). It follows that a quantitative comparison of bubble shapes between our numerical results and experimental measurements is not possible. It is however worth mentioning that, in their experiment at high Reynolds number (Re ≈ 400), [START_REF] Zenit | Measurements of the average properties of a suspension of bubbles rising in a vertical channel[END_REF] found the aspect ratio of ellipsoidal bubbles to decrease with increasing volume fraction (from χ = 1.5 at φ ≈ 0 to χ = 1.2 at φ ≈ 0.05). This trend is qualitatively similar to our numerical results presented in figure 10 for an ordered array of ellipsoidal bubbles with comparable aspect ratios, and more generally, is qualitatively similar to what we observe in our simulations of both ordered and free arrays of ellipsoidal bubbles.

To further support the idea that ordered arrays may be relevant to bubbly flows of practical interest, experimental data obtained by [START_REF] Garnier | Measurement of local flow characteristics in buoyancy-driven bubbly flow at high void fraction[END_REF], Riboux et al. 2002), d b0 = 3.4 mm q q q q q q q q q q q q q q q q q Figure 24. Prior experimental data for air bubbles rising in water: U (average drift velocity in meters per second) replotted versus φ 1/3 (with φ the gas volume fraction). In the legend, d b0 is the experimentally-determined bubble diameter at zero volume fraction.

(2010), and [START_REF] Colombet | Dynamics and mass transfer of rising bubbles in a homogenous swarm at large gas volume fraction[END_REF] for the air-water system at high Reynolds number (Re 100) are presented in figure 24 in the form U vs. φ 1/3 . Figure 24 bears a striking resemblance with figure 23: at moderate to fairly high volume fractions (say 0.2 φ 1/3 0.5, that is, 0.008 φ 0.13) the bubble velocity decreases linearly with φ 1/3 whereas in the dilute limit a different scaling law seems to hold. We have checked that neither a linear evolution with φ nor a law of the form U ∝ (1 -φ) n (with n a constant, possibly different for each data set) is compatible with the data presented in figure 24. It is also worth mentioning that some of these data sets are consistent with the notion that, in dilute bubbly suspensions, the bubble drift velocity may be higher than that of a single bubble (although again definitive conclusions can hardly be drawn on this point owing to experimental uncertainty on the values of the isolated bubble velocity).

Overall our numerical results as well as prior experimental data suggest that some properties of bubbly flows are sensitive to the presence of order, and that modelling a bubbly suspension by a cubic lattice of bubbles to investigate such properties is not irrelevant, except maybe at very high volume fractions, and as long as clusters are not formed.

Conclusions

The effect of volume fraction φ on the rise velocity and deformation of bubbles was first investigated when these are arranged in a cubic array. A non-monotonic behaviour of the rise velocity U at increasing volume fraction was obtained in the whole range of parameters considered from the DNS; this was supported by an analysis in the limit of weakly inertial suspensions of spherical bubbles. For low values of φ, "cooperative" wake interactions dominate and lead to an increase of U at increasing volume fraction, whereas the opposite behaviour occurs in the limit of large φ because of the predominance of "hindering" viscous interactions. These findings were supported further by comparison with the drag on a bubble behind another bubble when no other bubbles are present. A semiempirical law for the volume fraction dependence of the rise velocity, consistent with our numerical results even in the case of highly deformed bubbles, was also proposed.

The investigation of bubble shapes showed that ellipsoidal and skirted bubbles tend to become spherical at increasing volume fraction, and that the fore-and-aft asymmetry of isolated ellipsoidal bubbles is reversed for non-vanishing values of φ. An oblique motion of the bubbles was observed for certain parameter values, and supported by the abovementioned analysis. In this regime, the lift coefficient can be approximated by that of bubble pairs that are aligned vertically. The behaviour in this regime can be steady, oscillatory or chaotic, the latter arising if the horizontal extension of the bubbles wakes is large enough to allow interaction of bubbles with the wakes of neighbours which are not vertically aligned with them. A scenario explaining the transitions between these three regimes has been proposed.

The free rise of weakly deformed bubbles at moderate Reynolds number was then investigated for small and intermediate volume fractions. Simulations of free arrays of bubbles revealed that these share some common properties with ordered ones. Most notably, the drift velocity of free bubbles decreases linearly in φ 1/3 at moderately high φ whereas a different scaling law holds in the limit of low φ, as in ordered suspensions. This change of behaviour is compatible with available experimental data, and is believed to be responsible for the confusion in the literature regarding the form of empirical correlations in the context of corresponding asymptotic expressions. In addition, deformable bubbles were observed to become spherical as the volume fraction is increased, as in ordered arrays. We attribute the similarities between ordered and freely evolving suspensions to the fact that free bubbles were observed to keep the same neighbours for a long time, in agreement with prior work indicating that a certain degree of order is present in bubbly flows at comparable Reynolds number and volume fractions.

The present work is restricted to bubbles rising at a Reynolds number that is not more than 40. Beyond this, the dynamics of bubbly suspensions will be enrichened by the possibility of a single bubble already exhibiting path instability (e.g., [START_REF] Ern | Wake-Induced Oscillatory Paths of Bodies Freely Rising or Falling in Fluids[END_REF]). Such a study would require substantially larger computational resources, as comparatively thin boundary layers must be resolved. In addition, the present conclusions apply to perfectly homogeneous systems, perfectly monodisperse suspensions, and perfectly clean bubbles. Weak shear, polydispersity, and interface contamination may all have significant effects in real bubbly flows. The work presented herein forms the basis of an investigation of passive scalar transport in bubbly flows and of turbulent bubbly flows. which allows a stronger coupling between velocity and pressure than with co-colocated grids.

The governing equations are integrated in a coupled manner using a time-staggered discretization: the velocity components are computed at integer time steps while the pressure and the level-set function are computed at half-integer time steps. In what follows, ∆t is the time step, and the superscripts n and * are used to denote the current time iteration and some intermediate iteration, respectively. After initial conditions have been defined for the level-set, velocity, and pressure fields, the time integration algorithm proceeds iteratively through the following steps.

At the beginning of timestep t n , ψ n-1/2 , u n , u n-1 , and u n-2 are known.

Step 1: Advection of the level-set function. ψ is advanced from ψ n-1/2 to ψ n+1/2 according to equation (3.6) using the three-stage third-order TVD Runge-Kutta scheme [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF]:

ψ * = ψ n-1/2 + ∆t -L (u n , ψ n-1/2 ) + A (u n , ψ n-1/2 ) ψ n-1/2 , ( A 1)
ψ * * = 3 4 ψ n-1/2 + 1 4 ψ * + 1 4 ∆t -L (u n , ψ * ) + A (u n , ψ * ) ψ * , (A 2) ψ n+1/2 = 1 3 ψ n-1/2 + 2 3 ψ * * + 2 3 ∆t -L (u n , ψ * * ) + A (u n , ψ * * ) ψ * * , ( A 3)
where L (u, ψ) and A (u, ψ) are finite difference approximations of the advection term u • ∇ψ and of the source term A(u, ψ), respectively. In L (u, ψ), u is interpolated at the cell center with a second-order scheme and ∇ψ is computed using a fifth-order WENO scheme [START_REF] Jiang | Efficient Implementation of Weighted ENO Schemes[END_REF], as recommended by [START_REF] Salih | Some numerical studies of interface advection properties of level set method[END_REF].

In A (u, ψ), ∇ i ψ is calculated through a fourth-order centered scheme and ∇ i u j through a second-order centered scheme.

Step 2: Reinitialization of the level-set function. As already mentioned in section 3, the iterative approach introduced by [START_REF] Sussman | A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow[END_REF] where sgn is the sign function. An interesting feature of equation (A 4) is that the reinitialization of the level-set function starts near the interface and propagates outward: when this equation is solved up to pseudo-time T , d(x, τ = T ) is the signed distance function for all the points within distance T from the interface. Since it is important for ψ to be a signed distance function only inside the interfacial region of thickness 2ε, the reinitialization is not carried out to steady-state but only up to a given pseudo-time which must be at least equal to ε. Our algorithm is based on the second-order TVD Runge-Kutta scheme [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF] for the time integration of equation (A 4), which is carried out until τ = M ∆τ , where ∆τ is the artificial timestep, and M is a fixed number of iterations. (i) Initially This system is solved by an over-relaxed Gauss-Seidel method.

d 0 = ψ n+1/2 . ( A
Step 6: Corrector step for the velocity field. The intermediate velocity u * is corrected by the pressure gradient term to obtain u n+1 :

u n+1 = u * - ∆t ρ n+1/2 D pn+1/2 . ( A 23)
The algorithm then proceeds to timestep t n+1 .

A.2. Validation and sensitivity tests

The results of some benchmark tests performed with our code were already presented in section 3. We focus here on the fact that, in our simulations, volume conservation is enforced at each timestep by equation (A 12). This trick allows us to run simulations for a virtually infinite amount of time. However, as it expands the interface position uniformly (whereas volume losses may be localized), the numerical solution accuracy might be deteriorated. The excellent agreement we obtained for ellipsoidal bubbles with the results of [START_REF] Esmaeeli | Direct numerical simulations of bubbly flows. Part 2. Moderate Reynolds number arrays[END_REF] (see figure 1), who do not use such a correction, is therefore reassuring in this respect. Since it would be nonetheless desirable to evaluate the performance of our code without enforcing volume conservation, we disabled this fix and repeated the simulations of an ordered array of skirted bubbles (case C in table 1), a shape much more challenging to capture.

Simulations were carried out for φ = 1.6 %, ρ d /ρ c = 0.1 and µ d /µ c = 0.1. The time evolution of the bubbles drift velocity is shown in figure 25 for resolutions from 10 to 60 grid cells per bubble diameter, with and without volume conservation enforcement. The bubbles' shapes obtained at t mid = 10 d b /g and t end = 60 d b /g are reported in table 2 together with measurements of the volume variation (when the volume fix is not used) and of the volume correction (when the volume fix is used). When volume conservation is not imposed, the bubbles shrink inexorably, thus preventing the system from reaching a steady-state. As expected, the rate of mass gain decreases as the grid refines: at t end , the bubbles have almost completely disappeared with the coarsest grid (10 cells per bubble diameter), while the bubbles volume has reduced by 21 % with the finest grid (60 cells per diameter). But even with (reasonably) high resolutions, volume conservation remains problematic for long-time simulations. The volume correction we use to fix this issue is satisfactory, since it conserves the bubbles volume without affecting their dynamics: at short times (before volume loss becomes large), simulations with and without volume fix yield the same results. The error made when modifying the location of the interface by an amount ∆ψ remains much smaller than the O(∆x) global error in this region due to the finite thickness of the interface (max |∆ψ| ≈ 0.05∆x 2.8 , see table 2).

For completeness, we mention that in the simulations of ellipsoidal bubbles presented in figure 1a, the level-set correction is also found to be negligible (max |∆ψ| ≈ 0.1∆x 3.2 ), and that the volume loss in the absence of correction is very small (between t = 0 and t = 30 d b /g, 2.4 % and 0.1 % of the bubble volume are lost for resolutions of 20 and 40 grid cells per bubble diameter, respectively). Finally, we add that the magnitudes of the level-set correction in the simulations of ordered and free arrays are comparable.

Appendix B. The first effect of inertia in ordered arrays of bubbles and drops

In this section we derive the first effect of inertia on the steady drift velocity of an ordered suspension 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 linearized. Since all the particulates move with the same velocity, this configuration is equivalent to that of a cubic array of fixed particulates immersed in a viscous fluid moving with an average mixture velocity u = -U , and the problem becomes that of determining the hydrodynamic force, denoted f , exerted by the ambient fluid on a representative particulate of the array. It is customary to non-dimensionalize f with the magnitude of the Stokes-flow drag exerted on a single particulate in unbounded fluid to define a normalized force F :

F = f f 0,Stokes , ( B 1)
where f = |f | and f 0,Stokes is the drag force exerted on an isolated spherical fluid particulate in Stokes flow [START_REF] Hadamard | Mouvement permanent lent d'une sphere liquide et visqueuse dans un liquide visqueux[END_REF][START_REF] Rybczynski | Uber die fortschreitende Bewegung einer flussigen Kugel in einem zahen Medium[END_REF]:

f 0,Stokes = -2πµ * µ c d b U with µ * = µ c + 3/2µ d µ c + µ d . ( B 2)
Deviations of F from unity are induced both by hydrodynamic interactions and by the external fluid inertia. [START_REF] Hill | The first effects of fluid inertia on flows in ordered and random arrays of spheres[END_REF] obtained, under the assumption of φ 1, the first correction to F due to a small but non-zero Reynolds number for a cubic array of rigid spheres by matching the far-field fundamental periodic solution of the Oseen equations to the nearfield solution of the Stokes equations past an isolated rigid sphere. Their derivation can be extended to cubic arrays of bubbles and drops by replacing the inner solution for a rigid sphere by that for a fluid sphere, which yields the following linear system from which f is determined:

f - 2.8373 3 µ * d b h f -2πµ * d b h f • S = f 0,Stokes , ( B 3)
where f 0,Stokes is given by equation (B 2) and h is the lattice spacing. We have here introduced S, the dimensionless symmetric tensor given by are the vectors in the reciprocal lattice defined by the primitive vectors b i . Therefore, at finite Reynolds numbers, the force exerted on a particulate and the drift velocity have, in general, different directions. Note that in the limit of Re h → ∞ (that is, φ → 0), one recovers the result of [START_REF] Brenner | The resistance to a particle of arbitrary shape in translational motion at small Reynolds numbers[END_REF] for the first inertial contribution to the normalized drag on a single fluid particulate translating in an unbounded fluid

S = k * =0 Re 2 h (U * • k * ) 2 (I -k * k * /k * 2 ) (2π) 4 k * 6 1 + Re 2 h (2π) 2 k * 4 (U * • k * ) 2 , ( B 4
1 F = 1 - 1 8 µ * Re for φ → 0. (B 6)
We first consider the specific situation where the drift velocity is aligned with a primary axis of the cubic array. In that case the off-diagonal components of S are zero and the total force f is parallel to the drift velocity U . In the limit Re h → 0, the solution of equation (B 3) reads 1 F = 1 -1.1734µ * φ 1/3 -0.0050µ We now look for the existence of non-vertical solutions for the drift velocity in the case where the system is buoyancy-driven. The problem is then reversed: the hydrodynamic force acting on the bubble is prescribed (it opposes buoyancy) and one wants to determine the drift velocity of the bubbles. We further assume that gravity is aligned with an axis of the lattice and writes g = -ge 3 . The buoyancy force is then

f buoy = 1 6 πd 3 b (ρ c -ρ d )ge 3 = 1 6 π µ 2 c ρ c
Ar 2 e 3 . (B 11)

Replacing f = -f buoy in equation (B 3) yields the following nonlinear system of dimensionless equations

U * 1 = 1 96π 3 Ar 2 Re h s 13 (B 12a) U * 3 = Ar 2 12Re h 1 µ * h d b - 2.8373 3 -2π Re 2 h (2π) 4 s 33 (B 12b) U * 1 + U * 2 + U * 3 = 1 (B 12c)
where s ij = (2π) 4 S ij /Re 2 h and from which U = U i e i can be determined.

Figure 1 .

 1 Figure 1. Time evolution of the drift velocity of arrays of deformable bubbles rising at moderate Reynolds numbers (case E1). Solid line: prior DNS of Esmaeeli & Tryggvason (1999). Non-solid lines: present DNS. (a) Ordered array (N b = 1, φ = 13 %, ρ d /ρc = µ d /µc = 0.1), different resolutions (d b is the bubble volume-equivalent diameter, ∆x the grid spacing). (b) Free array (N b = 8, φ = 6.5 %, ρ d /ρc = µ d /µc = 0.05), with three realizations of the flow.

Figure 3 .

 3 Figure 3. Steady drift velocity of an ordered array of spherical bubbles, normalized by the terminal velocity of an isolated bubble in Stokes flow conditions, as a function of volume fraction for small Archimedes numbers. (a) The effect of small but finite Ar , from Oseen-flow analysis. (b) Comparison of analytical and numerical solutions for Ar = 5.03 (case S1). Analytical solutions are obtained from (4.3). In (b), --: analysis; •: DNS; •: isolated bubble, estimated from Mei et al. (1994); ----: numerical fit of the form of (4.8) matching DNS and isolated-bubble data; and • • • • • •: numerical data and fit by a Richardson-Zaki relation from Sankaranarayanan et al. (2002).

q

  Ar = 0.15 (S0) Ar = 3.00 Ar = 5.03 (S1) Ar = 10.0 (S2) Ar = 15.3 (S3) Ar = 27.2 (S4) Ar = 40.7 (S5)

  Figure 5. Vertical component of the liquid velocity normalized by the bubble drift velocity in a vertical symmetry plane passing through the center of a bubble in an ordered array configuration at φ = 0.2 %. Gravity is pointing downward (g = -ge3). Increasing Archimedes numbers from left to right, and increasing Bond numbers from top to bottom.

Figure 6 .Figure 7 .

 67 Figure6. Time signal of the normalized drift velocity of an ordered array of bubbles at small volume fraction (φ = 0.1 %). The dashed vertical lines denote t = nh/U0, the characteristic times at which the bubbles enter into the wake of their n th preceding neighbour. For clarity only short times are shown (the steady-state is reached for t ≈ 50h/U0).

  Figure 8. Normalized steady drift velocity of an ordered array of deformable bubbles as a function of volume fraction for various Bond numbers. (a) Ar ≈ 30: spherical and ellipsoidal bubbles. (b) Ar ≈ 15: spherical and dimpled ellipsoidal-cap bubbles. Open symbols: present DNS. Black crosses: prior DNS of Bunner & Tryggvason (2002a) (Ar = 29.7, Bo = 0.98, ρ d /ρc = µ d /µc = 0.02). Dashed lines: numerical fits of the form of (4.8) matching DNS and isolated-bubble data.

  Figure 9. Fitted coefficients for U (φ)/U0 as given by (4.8): (a) as a function of the Archimedes number for Bo = 0.38, (b) as a function of the Bond number for Ar ≈ 15 (open symbols) and for Ar ≈ 30 (filled symbols). The coefficients are obtained by least-squares fits of DNS and isolated-bubble data.

Figure 10 .

 10 Figure 10. Steady aspect ratio (a,c,e) and shape coefficients (b,d) of bubbles in an ordered array configuration as a function of volume fraction for: (a,b) case E1 (Bo = 2.0, Ar = 29.9); (c,d) case E2 (Bo = 5.0, Ar = 30.0); (e) case C (Bo = 243, Ar = 15.2). Bubbles shapes are shown in grey for the highest and the lowest simulated volume fractions. Open symbols: DNS.Filled symbol: isolated bubble (aspect ratios estimated from[START_REF] Loth | Quasi-steady shape and drag of deformable bubbles and drops[END_REF] for cases E1 and E2 and from[START_REF] Hua | Numerical simulation of 3D bubbles rising in viscous liquids using a front tracking method[END_REF] for case C; shape coefficients estimated from[START_REF] Ryskin | Numerical solution of free-boundary problems in fluid mechanics. Part 2. Buoyancy-driven motion of a gas bubble through a quiescent liquid[END_REF]).

Figure 11 .

 11 Figure 11. Influence of volume fraction on the deformation of an ordered array of oblate ellipsoidal bubbles (case E2): relative velocity (u rel = uu d ) and pressure fields in a vertical symmetry plane passing through the center of a bubble. Increasing volume fractions from left to right: (a,c) φ = 0.1 %, U/U0 = 1.34, and (b,d) φ = 13 %, U/U0 = 0.73. Top row: magnitude of the liquid relative velocity in the entire unit cell. Bottom row: total pressure (including the mixture-average hydrostatic component) near and inside a bubble. Gravity is pointing downward (g = -ge3). The black lines show the interface location.

Figure 12 .

 12 Figure 12. Influence of volume fraction on the deformation of an ordered array of dimpled ellipsoidal-cap bubbles (case C): vertical component of the relative velocity (u rel = uu d ) in a vertical symmetry plane passing through the center of a bubble. Increasing volume fractions from left to right: (a) φ = 0.2 %, U/U0 = 1.22, and (b) φ = 2.4 %, U/U0 = 1.17. Same conventions as in figure 11. Only the region near and inside a bubble is shown.

Figure 13 .

 13 Figure 13. Subset of analytical solutions obtained from solving (B 3) for Ar = 20 and ρ d /ρc = µ d /µc = 0: bubble drift velocity horizontal (a) and vertical (b) components, given as Reynolds numbers Rei = Uiρcd b /µc, as a function of volume fraction. Gravity is pointing in the -e3 direction. • • • • • •: vertical rise (U1 = U2 = 0); ----: oblique rise with U2 = 0; --: oblique rise with |U2| = |U1|.

Figure 14 .

 14 Figure 14. Time signals of the bubble drift velocity components (given as Reynolds numbers Rei = ρcUid b /µc), with U = Uiei and gravity pointing in the -e3 direction, in the three regimes of motion: steady oblique rise (a,b), oscillatory oblique rise (c), and chaotic oblique rise (d). These regimes are obtained for (a) case E1 at φ = 0.8 %, (b) case S5 at φ = 13 %, (c) case S5 at φ = 3.8 %, (d) case S5 at φ = 0.5 %.

Figure 16 .

 16 Figure16. Average lift coefficient as a function of volume fraction for case E1 (grey triangles) and case S5 (black squares). Open symbols: for a bubble of an ordered array, present DNS. Filled symbols: for the trailing bubble of a pair of spherical bubbles within the same conditions (orientation, Reynolds number, separation distance) as two vertically-aligned bubbles of the ordered array, with (large symbols) and without (small symbols) accounting for the interaction with the wake of the leading bubble, from[START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF].

Figure 17 .

 17 Figure 17. Liquid vorticity horizontal component ω h = ω • e h in the vertical plane passing through the center of a bubble and containing its drift velocity, for ordered arrays of bubbles rising steadily in an oblique direction: (a) case E1, φ = 0.8 %; (b) case S5, φ = 6.5 %. Positive (negative) values upstream of the bubble would give a positive (negative) contribution to the lift force as modeled by (4.12). The thin black arrow shows the direction of bubble drift velocity. The thick arrows show the magnitudes and directions of the drag (green arrows) and lift (orange arrows) forces scaled by the buoyancy force.

Figure 20 .

 20 Figure 20. Influence of the number of free bubbles on their spatial distribution within a periodic unit cell: typical instantaneous snapshots for N b = 4 (left) and N b = 8 (right), for case E1 at φ = 3.8 %.

Figure 23 .

 23 Figure 23. Influence of volume fraction on the normalized drift velocity of freely rising bubbles: comparison between simulations and experiments. •: present DNS (N b = 8, Ar = 29.9, Bo = 2.0, Re0 = 31); ×: experiments of Martinez-Mercado et al. (2007) (Ar = 26.3 ± 1.6, Bo = 0.25 ± 0.02, Re0 = 29.7 ± 3.1). The terminal velocities U0 in unbounded liquid have been estimated from the correlation of[START_REF] Loth | Quasi-steady shape and drag of deformable bubbles and drops[END_REF] for numerical simulations (ellipsoidal bubbles) and from the correlation of[START_REF] Mei | A note on the history force on a spherical bubble at finite Reynolds number[END_REF] for experiments (spherical bubbles).

q

  Colombet et al. (2015), d b0 = 2.1 mm Riboux et al. (2010), d b0 = 1.6 mm Riboux et al. (2010), d b0 = 2.1 mm Riboux et al. (2010), d b0 = 2.5 mm Garnier et al. (

  is used to reinitialize ψ. This consists in solving for an artificial time τ∂ d ∂τ = sgn(ψ)(1 -|∇d|), with d(x, τ = 0) = ψ(x) (A 4)

Figure 25 .

 25 Figure 25. Effects of resolution and of volume conservation enforcement on the rise of an ordered array of strongly deformed bubbles (case C at φ = 1.6 % with ρ d /ρc = 0.1 and µ d /µc = 0.1).

  ψ mid = 0 ψ end = 0 max |∆ψ|/∆x 2.8 ψ mid = 0 ψ end = 0 (V end d Sensitivity analysis of the effects of resolution and of volume conservation enforcement on the shape of strongly deformed bubbles (case C at φ = 1.6 % with ρ d /ρc = 0.1 and µ d /µc = 0.1, ordered array). Bubble shape (2D cut-off in a symmetry plane) at t mid = 10 d b /g and t end = 60 d b /g, maximum magnitude of the level-set correction, and volume relative variation between t 0 and t end (∆x is the grid spacing, ∆ψ is the level-set correction, V d is the disperse phase volume).

  ) with U * = U /U , Re h = ρ c U h/µ c , and k * = kh where k = n 1 b 1 + n 2 b 2 + n 3 b 3 n 1 , n 2 , n 3 = 0, ±1, ±2, . . . (B 5)

Table 1 .

 1 Simulated regimes: Bo, Ar , and Mo = Bo 3 /Ar 4 are input parameters (with ρ d /ρc = 10 -3 and µ d /µc = 10 -2

	case Bo Ar Mo	shape	Re0	We0	χ0
	S0 0.38 0.15 1.00×10 2 spherical	1.94×10 -3 6.11×10 -5 1.000
	S1 0.38 5.03 8.60×10 -5 spherical	1.80	4.88×10 -2 1.007
	S2 0.38 10.0 5.49×10 -6 spherical	5.94	0.134	1.015
	S3 0.38 15.3 1.00×10 -6 spherical	12.1	0.236	1.024
	S4 0.38 27.2 1.00×10 -7 spherical	31.4	0.507	1.064
	S5 0.38 40.7 2.00×10 -8 spherical	62.5	0.897	1.124
	E1 2.0 29.9 1.00×10 -5 ellipsoidal	31	2.1	1.32
	E2 5.0 30.0 1.54×10 -4 ellipsoidal	26	3.8	1.62
	C	243 15.2 2.66×10 2 dimpled ellipsoidal-cap 7.77	63.2	1.89

  Re 2 φ -1/3 + O(Re 4 a reasonable estimate of S for any Re h and F can be approximated at any At steady state, the total force exerted by the fluid on the particulate is balanced with the buoyancy force so the solution of the sedimentation problem relates to F through

	φ	1 by
		1 F	≈ 1 -	1 8	µ * Re -1.1734µ * φ 1/3 +	25 8	µ * Reφ 1/3 Re + 25φ 1/3 .	(B 9)
	the identity
					U U 0,Stokes	=	1 F	.	(B 10)
								h )	for φ 1/3	O(Re). (B 7)
					S ≈	Re h 1.53Re h 16π + (2π) 4	(B 8)

*

For intermediate values of Re h , the longitudinal component of S, denoted S = U * •S•U * , is needed and can be computed numerically. In practice, the simple expression provides

† This body force is equivalent to the average hydrostatic pressure gradient that would be generated by the base of a flow container to balance the total gravitational force on the mixture.

† The denomination "ellipsoidal bubble" refers to a bubble that is approximately spheroidal with weak fore-and-aft asymmetry, but does not mean that the bubble shape is strictly ellipsoidal.

The authors acknowledge stimulating discussions with F. Risso. 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 A. Discretization schemes and additional validation tests

Details on the numerical methods used in our study, together with the results of some benchmark and sensitivity tests of our code, are provided hereinafter.

A.1. Algorithm and discretization schemes

Spatial discretization relies on a finite difference/finite volume approach on a fixed, staggered, Cartesian grid. Scalar variables (level-set, pressure) are located at cell centers, which discrete coordinates are denoted with subscripts (i, j, k), and the three components of velocity are stored on cell-face centers (i + 1/2, j, k), (i, j + 1/2, k) and (i, j, k + 1/2), (ii) Then for m = 0 to m = M :

In practice we use the standard value ∆τ = 0.5∆x, and set M = 5. In the above algorithm, R(d) represents the discretization of the spatial term sgn(ψ)(1-|∇d|) devised by [START_REF] Russo | A remark on computing distance functions[END_REF], which reads in three dimensions

Σ ∆x is the set of points located within one grid point from the zero-level of d 0 , where D i,j,k is computed by

) is an upwind discretization of |∇d| -1 computed with a finite-difference second-order ENO scheme [START_REF] Harten | Uniformly high order accurate essentially non-oscillatory schemes[END_REF].

Step 3: Correction of the level-set function. To enforce volume conservation, the iso-contours of ψ n+1/2 are shifted. ψ n+1/2 is then replaced by

where V d is the volume of the disperse phase calculated from (Ω is the computational domain):

and S i is the surface area of the interfaces between the two phases, obtained from:

where

is the smoothed version of the delta function, defined as the derivative of H ε with respect to ψ. The basis for this correction is discussed further in section A.2.

Step 4: Predictor step for the velocity field. A provisional mid-step velocity u * is computed from u n by omitting the pressure gradient term in the momentum conservation equation and by using a mixed Crank-Nicolson/Adams-Bashforth time-stepping scheme:

where C , V and F are spatial discretizations of the advection, viscosity and surface tension terms respectively. The advection term is extrapolated at t n+1/2 using a thirdorder Adams-Bashforth scheme:

where C is the discretization of u • ∇u based on a finite-difference fifth-order WENO scheme for ∇u with a second-order interpolation of u when needed. The contribution V is the discretized version of ∇ • µ(∇u + ∇u T ), its component in the p-direction expands in

where D are discrete spatial derivatives calculated using second-order central-difference and interpolation schemes. The temporal discretization of the p-component of the viscous contribution employs a semi-implicit Crank-Nicolson scheme for the four terms involving the derivatives of u p , and an explicit third-order Adams-Bashforth scheme for the two terms involving the derivatives of u q =p . This writes

and

(A 20) The surface tension term is computed at t n+1/2 directly from ψ n+1/2 , i.e.,

where F is the space discretization of γκ∇H ε : κ is obtained from a second-order centered finite-volume discretization of equation (3.3), and ∇H ε is computed using a secondorder central differencing scheme. Note that the formulation of the singularity as ∇H ε is preferred over the usual form δ ε (ψ)∇ψ because it effectively reduces the amplitude of so-called spurious currents [START_REF] Meland | Reduction of parasitic currents in level-set calculations with a consistent discretization of the surface-tension force for the CSF model[END_REF], some parasitic currents arising from an inconsistent discretization of the surface tension force and the pressure gradient. The resulting linear system is solved iteratively for u * using a hybrid Jacobi/Gauss-Seidel algorithm.

Step 5: Density-weighted Poisson equation for pressure. The pseudo-pressure pn+1/2 is obtained from

where D is the second-order central-difference discretization of the spatial derivatives.