Main

The VOF model resolves the transient motion of the gas and liquid phases using the Navier-Stokes equations, and accounts for the topology changes of the gas-liquid interface induced by the relative motion between the dispersed gas bubble and the surrounding liquid. The finite-difference VOF model uses a donor-acceptor algorithm2 to obtain and maintain an accurate and sharp representation of the gas-liquid interface. The VOF method defines a fractional volume or ‘colour’ function c (x, t), which is a function of position vector xand time, t. The colour function indicates the fraction of the computational cell that is filled with liquid, and varies between 0, if the cell is completely occupied by gas, and 1, if it consists only of liquid.The location of the bubble interface is tracked in time by solving a balance equation for this function, (∂c (x, t)/∂t)+∇˙(u c (x, t) = 0, where u is the velocity vector.

The liquid and gas velocities are assumed to equilibrate over a very small distance, and essentially uk = u at the bubble interface, where the subscript k corresponds to either liquid or gas. The mass and momentum conservation equations are considered to be homogeneous. The continuum surface force model3 is used to model the force arising from surface tension acting on the gas-liquid interface. In this model, the surface tension is modelled as a body force Fsf, which is non-zero only at the bubble interface and is given by the gradient of the colour function, Fsf = σκ(x)∇ c (x, t), where σ is the surface tension and κ(x) is the local mean curvature of the bubble interface κ(x, t) = −∇(n/|n|), where n is the vector normal to the bubble interface.

All our simulations were carried out in a rectangular column using a uniform two-dimensional cartesian-coordinate grid. The front of the two-dimensional rectangular grid is formed by the x -z plane. The front and rear faces of the column are modelled as symmetry planes. At the two walls, the no-slip boundary condition is imposed. The column is modelled as an open system, so the pressure in the gas space above the initial liquid column is equal to the ambient pressure (101 kPa). Hybrid differencing was used for the convective terms in the equations, and upwind differencing was used for time integration. The time steps used in the simulations were usually 0.0003 s or smaller. To counteract excessive smearing of the liquid-gas interface by numerical diffusion, a surface-sharpening routine was invoked. This routine identifies gas and liquid on the ‘wrong’ side of the interface, and moves it back to the correct side, while conserving the volume of the respective phases.

To avoid dissolution of the bubble by surface sharpening, we ensured that the area of each bubble encompassed a few hundred cells. For simulations of bubbles smaller than 12 mm, we found that the rich dynamic features could be captured only by allowing at least 800 grid cells per bubble cross-section, resulting in grid sizes as small as 0.125 mm for the smallest bubble size we could simulate with adequate accuracy. To simulate the rise of spherical cap bubbles, typically found with sizes above 17 mm, a grid size of 1 mm was found to be small enough; in this case, there were more than 300 grid cells per bubble cross-section.

All simulations were carried out using the parallel version of CFX 4.1c (AEA, Harwell, UK) running on a Silicon Graphics Power Challenge machine with six R8000 processors. This package is a finite volume solver that uses body-fitted grids, which were non-staggered, and all variables were evaluated at the cell centres. We used an improved version of the Rhie-Chow algorithm4 to calculate the velocity at the cell faces, and the SIMPLEC algorithm5 to obtain the pressure-velocity coupling. Simulating the rise of a bubble 7 mm in diameter for 0.75 s in a column 0.025 m wide and 0.09 m high, involving 144,000 grid cells, required about two weeks of computer processing time.

Typical bubble trajectories are shown in Fig. 1. Bubbles of 4 and 5 mm in diameter meander up the column with large lateral movements, conforming with the observations of Leonardo da Vinci. The 7-mm bubble oscillates from side to side when moving up the column. Bubbles 8 and 9 mm in diameter behave rather like jelly-fish. The 12-mm bubble moved like a bird flapping its wings; we could find no mention of this type of bubble motion in the literature, but we were able to observe it in our experiments.

Figure 1: Two-dimensional volume-of-fluid simulations of the rise trajectories of bubbles in the size range 4-20 mm.
figure 1

The bubble sizes are: a, 4 mm; b, 5 mm; c, 7 mm; d, 9 mm; e, 12 mm; f, 20 mm. In all simulations, the gas phase was air (density, 1.29 kg m-3; viscosity, 1.7×10-5 Pa s) and the liquid was water (density, 998 kg m-3; viscosity, 10-3 Pa s; surface tension, 0.072 N m-1). The Morton number6 for all simulations was 2.63×10-11. The Eötvös numbers6 for the various bubble sizes were: a, 2.2; b, 3.4; c, 6.7; d, 11.0; e, 19.6; f, 54.5. Animations of the simulations can be seen on our website (http://ct-cr4.chem.uva.nl/single_bubble/).

Full size image

The two-dimensional VOF simulations give velocity values for the rising bubbles that are roughly half those measured experi-mentally in cylindrical columns. A three-dimensional simulation would be required to determine rise velocity quantitatively and to reproduce precisely the experimental observations of Da Vinci. Such a simulation would require a grid with a few million grid cells and the use of a massively parallel computer.