Prediction of initial motion of a gas bubble in liquids

Prediction of initial motion of a gas bubble in liquids M. Nilmani Department

of Metallurgy,

University

of Strathclyde,

Glasgow,

Scotland,

UK

T. T. Maxwell Department

of Mechanical

Engineering,

Auburn

University,

Auburn,

Alabama,

USA

D. G. C. Robertson Department of Metallurgy, London SW7, UK

imperial

College of Science

and Technology,

D. B. Spalding Department of Mechanical London SW7, UK (Received

January

Engineering,

Imperial

1980; revised September

College of Science

and Technology,

1980)

Numerical predictions are presented for the motion and distortion of a single gas bubble rising through the liquid. The computations were made with an implicit finite-difference procedure which solves the transient equations of motion throughout the bubble and the liquid, such that the free surface between the gas bubble and the liquid is not a boundary of the computational domain. The predictions compare well with the experimental results of others. Computations are presented for bubble sizes from 0.02 to 0.05 m radius and for bubbles of different gas densities rising in liquids of different densities. Surface tension effects are neglected. Introduction A large number of processes in chemical, mechanical and metallurgical engineering are based on mass transfer between a liquid and a gas across their common interface. This common interface is often provided by gas bubbles rising through a liquid. Experimentally3ps it has been found that under conditions of high gas flow rate large gas bubbles are generated at submerged nozzles. After severance from the nozzle these disintegrate within a rise path of about three bubble diameters to produce a column of smaller bubbles in the liquid. Disintegration of the bubble results from the asymmetric pressure distribution around the bubble as it accelerates upwards after detaching from the nozzle. The present paper predicts the distortion of large spherical bubbles of different gas densities during their initial motion in liquids of different densities.

Past work The change in the shape of an accelerating and freely rising bubble has been studied by Walters and Davidson.1*2 The theory proposed for the bubble rise and distortion was

24

Appl. Math. Modelling, 1981, Vol 5, February

based on the assumption of irrotational motion in the liquid around the bubble and constant pressure within the bubble. The theoretical approach was developed first for a two-dimensional1 system and then extended to a threedimensional2 system. Walters and Davidson’, 2 carried out experiments to measure the vertical distortion of rising bubbles with time and data are available for both the two-dimensional and three-dimensional systems. The equations derived predict the initial rate of bubble motion and that a tongue of liquid will rise into the base of the accelerating bubble. However, it fails to predict the observed behaviour in the later stages when the tongue approaches the top surface.

Present contribution The present paper reports the results of a theoretical analysis for a gas bubble in a liquid rising freely under gravity from rest. Both axisymmetric and plane flows are considered. It compares predictions with the available experimental data and studies the effect of gas and liquid densities on the distortion behaviour of the bubble, and so provides an insight into the physical processes involved. The main assumption is that the surface tension effects have been 0307-904X/81/010024-05/$02.00 0 1981 IPC Business Press

Gas bubble motion: ignored in the calculations since in large bubbles the effect of surface tension can be neglected compared to the dominant forces of buoyancy and form drag. In the later stages of bubble deformation, when part of the bubble interface acquires a small radius of curvature, surface tension may become important. A spherical bubble centrally placed in a cylindrical container 0.3 m in height and 0.3 m in diameter (Figure I) was considered. The water depth in the container was 0.26 m and the rest of the container was occupied by air at a pressure of one atmosphere. The bubble was allowed to rise freely from rest in the stagnant liquid under the influence of gravity. The effect of liquid and gas densities were studied.

Mathematical

formulation

and solution

to be solved are:

az

+ ; -=R a@)

r

(1)

at-

CL

-

(2) where B, is the body force in the axial direction represents any other sources of axial momentum.

where S, represents any additional momentum.

OLit-

Gas

-

________________--_-----

’

3.26b

0 PortKle

____

string 2

I

Llauld

1

F‘article ‘s ,trmg

,,-$1’ [ Gas

,-

T 2a

1

015M

7

and S,

sources of radial

where rp is the position vector of a particle and Up is the particle velocity.

Solution

f

Figure

Axial momentum

dr;,_- UP dt

Volume continuity F

of fluid volume,

Particle movement The free surface marker particles move in accordance with the equation:

equations

The equations

et al.

procedure

An Eulerian finite-difference procedure was used to compute the hydrodynamics in which the free-surface is defined by the motion of a set of ‘Langrangian particles’. The conventional mass continuity equation is rearranged to form a volumetric or bulk-continuity equation. The use of this bulk-continuity relation allows the hydrodynamic variables to be computed over the entire flow domain. Thus, the free-surface boundary conditions are imposed implicitly and the formulation of the problem is greatly simplified.

Governing

where R is the rate of reduction - (llp)(WW

M. Nilmani

System considered

for numerical

predictions

procedure

The solution procedure is based on an implicit finite difference scheme SIMPLE7>a (Semi Implicit Method for Pressure-Linked Equations). SIMPLE is an iterative procedure in which the velocity fields are calculated from an estimated pressure field, and pressure corrections are then computed from the continuity relations. These ‘pressure corrections’ are used to adjust both the pressure and the velocity fields. The new fields are used as ‘guesses’ for the next iterations; after several iterations, the hydrodynamic fields are obtained. In the present case of a transient flow, this iterative procedure is repeated at each time step and the guessed values are taken as the results of the previous step. To compute free-surface flows, two modifications to the SIMPLE procedure are necessary. The first modification is the introduction of GALA9 (Gas And Liquid Analyser). The central concept of GALA is the replacement of the conventional mass continuity equation with a volume continuity equation. The use of GALA permits the computation of the hydrodynamic variables over the entire flow domain (throughout the liquid and gas regions), so that the pressure and velocity are continuous across the free surface. Thus, because the free surface is not a boundary of the computation, the complex free-surface boundary conditions are implicitly imposed. The second modification is the addition of a set of massless particles to ‘track’ the free surface. These particles are moved, at the end of each step in time,

Appl.

Math.

Modelling,

1981,

Vol

5, February

25

Gas bubble motion:

/VI. Nilmani

et al.

with the local fluid velocities and the particle locations are then used to determine the local fluid densities. Introducing GALA enhances the stability and convergence and so no problems are encountered in this respect. To complete the problem specification both initial conditions and boundary conditions are needed. The necessary initial conditions include the initial values of the flow variables and the initial positions of the free surface marker particles. The necessary boundary conditions are the values of the flow variables along the boundary of the calculation domain as functions of time. To simulate the gas compressibility:

which is the rate of reduction of fluid volume within a control volume, must be represented in the R terms of equation (1). The substantial derivative of the density describes the density of the moving fluid elements, but is being integrated over a fixed control volume. This practice is not entirely consistent; however, the important contribution of this term is the net expansion of the fluid in the vicinity of the control volume. If the location of this expansion is displaced a small amount from its true location, the effect will be secondary. R is computed from the reversible adiabatic relationship, pVk = constant.

Results and discussion Figure 2 shows the experimental

points of Walters and Davidson L * plotted as the change in dimensionless vertical diameter of the bubble D/2a vs. dimensionless time N1’2 = t(g/o)“* for a three-dimensional system. Computed values obtained from grid 22 x 17 and grid 37 x 22 are represented by curves A and B respectively. Details of the finite grid distributions are given in Tables I and 2. The computed values depend on the grid size but since the program is expensive to run it was not economical to use a still finer grid, where an even closer fit can be expected. The computer time to run the calculations on the CDC 6400 at Imperial College were 903.6 s for the 22 x 17 grid to solve 0.3040 s

0

1

I

I

I

02

04

06

08

:

I

I

10

12

N’k

Figure 2 Vertical predictions; f---j, (01, experimental’

26

Appl.

distortion Walters

of bubble. PL = 1000 kg/m”. I--), et a/.’ A, grid 22 X 17; B, grid 37 X 22;

Math. Modelling,

1981,

Vol 5, February

Table 1

Details

of the finite

grid

distribution

(grid

size 37 X 22)

Coordinates 2, (ml 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37

Tab/e 2

r, (m)

0. 6.000E-03 1.800E-02 2.600E-02 3.400E-02 4.000E-02 4.600E-02 5.200E-02 5.800E-02 6.400E-02 7.000E-02 7.600E-02 8.200E-02 8.800E-02 9.400E-02 1 .OOOE-01 l.O60E-01 1 .I 20E-01 l.l80E-01 1.240E-01 1.300E-01 1.360E-01 1.420E-01 3.500E-01 1.580E-01 1.660E-01 1.740E-01 1.820E-01 1.920E-01 2.020E-01 2.100E-01 2.200E-01 2.340E-01 2.500E-01 2.700E-01 2.900E-01 3.000E-01

Details

of finite

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

grid distribution

(grid

0. 2.000E-03 6.000E-03 1 .OOOE-02 1.400E-02 1.800E-02 2.200E-02 2.800E-02 3.400E-02 4.000E-02 4.600E-02 5.200E-02 5.800E-02 6.400E-02 7.000E-02 8.000E-02 9.000E-02 1 .OOOE-01 1 .140E-01 1.300E-01 1.500E-01 1.600E-01

size 22 X 17)

Coordinates I, (ml

2, (ml 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

0. 7.500E-03 2.250E-02 3,75OE-02 5.250E-02 6.750E-02 8.250E-02 9.750E-02 1 .125E-01 1.275E-01 1.425E-01 1.575E-01 1.725E-01 1.875E-01 2.025E-01 2.175E-01 2.325E-01

18 19 20 21 22

2.475E-01 2.625E-01 2.775E-01 2.925E-01 3.000E-01

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

0. 5.000E-03 1.500E-02 2.500E-02 3.500E-02 4.500E-02 5.500 E-02 6.500E-02 7.500E-02 8.500E-02 9.500E-02 l.O50e-01 1 .150E-01 1.250E-01 1.350E-01 1.450E-01 1.500E-01

of real time, and 1033 s for the 37 x 22 grid to solve 0.093 of real time. Figure 3 shows the effect of the size of bubble on the bubble distortion. Computed values of all three air bubbles

s

Gas bubble motion:

M. Nilmani

et al.

2

.8 F

0

0.4

08

1.2

16

I

2

Trne

N ‘12

~0 000

5

Time=

00125

Time

= 0024s

Figure 3 Effect of bubble volume on vertical distortion. (0). (a), 179.6 ml; (*I, 268 ml; (01,523.O ml. Grid used 37 X 22

I

0 51

’

TlmezO048s

Figure 4 Distortion of freely rising gas bubble (7 cm dia.) in liquids. Effect of gas-liquid density ratio on vertical distortion bubble after 0.096 s is shown

of

in water fall on the same curve, showing that the vertical distortion of the bubble is independent of the size of the bubble. Walters and Davidson’s experimental data agree with the prediction. The effect of gas-liquid densities on bubble distortion is shown in Figure 4. Both gas and liquid densities were varied in the computation. Calculation shows the gas-liquid density ratio to be an important quantity, governing distortion behaviour of the bubble. The gas-liquid density ratio has been plotted against the dimensionless vertical diameter of the bubble at N”’ = 1.6069. Below density ratio values of - 1000 plot shows that with decrease of density ratio the rate of bubble distortion decreases. However, at density ratio above - 1000 program fails to predict the expected behaviour. D&a obtained for an air bubble in water has been given in the plots of bubble shapes and the velocity vectors at different time intervals (Figure 5). The figure confirms that the tongue of liquid slowly rising into the bubble touches the top interface of the bubble after a displacement of two bubble diameters. This has been found experimentally5 during jetting of air into water where the large rising bubbles shatter within a height of two to three bubble diameters, although of course the two cases are not strictly comparable. Figure 6 shows the displacement of the bubble centroid against the dimensionless time IV”‘. Predictions fit well with the experimental results for a density of 1000 kg/m3. Bubbles rise initially with twice the acceleration due to gravity. Figures 7 and 8 show the distortion and acceleration behaviour of a two-dimensional bubble 0.025 m in diameter in the form of a cylinder with its axis horizontal. Calcula-

Time

Figure (179.6

Time

: 0072s

I

0

= 00845

Tfme

5 Bubble shape and position at different ml) in water

02

04

L

06

,

08

1

z 0096s

times. Air bubble

I

10

12

14

N”z Figure 6 Displacement of bubble centroid as function of dimensionless time. (-_), predictions; (o), experimental’; C---j, Walters et al.‘. Grid used 37 X 22

Appl.

Math.

Modelling,

1981,

Vol 5, February

27

Gas bubble motion:

M. Nilmani

et al.

Conclusions The present numerical technique referred to as SPLASH (Simple Program For Liquid Air System Hydrodynamics) developed for two-phase flow systems has been found very useful in predicting the behaviour of bubbles inside liquids.

Acknowledgements Concentration, Heat and Momentum for the use of the SPLASH computer

A I

06

1

0

02

I

04

I

I

06

05 N "2

07

References I 08

Walters, J. K. and Davison, J. F. J. FluidMech. 1962, 12,408 Walters, J. K. and Davidson, J. F. J. Fluid Mech. 1963, 17, 322 Leibson,l.etaL A.I.Ch.E.J. 1956,2, 296 Wraith, A. E. ‘Advances in extractive metallurgy and refining’ (M. J. Jones, Ed.), IMM, London, 1972, 303-316 Nilmani, M. PhD Thesis, University of London, 1978 Maxwell, T. T. PhD Thesis, University of London, 1977 Caretto, L. S. et al. In Proc. third Int. Conf. Numer. Meth. Vol 2, Springer-Verlag, Heidelberg, 1973 (Ed. J. Ehlers et al.) Patankar, S. V. and Spalding, D. B. Int. J. Heat Mass Transfer 1972,15,1787 Spalding, D. B. In ‘Heat transfer and turbulent buoyant conPublishvection’ (Ed. D. B. Spalding and N. Afgan), Hemisphere ing Corporation, Washington, D.C., 1977

Figure 7 Vertical distortion of two-dimensional bubbles. PL = 1000 kg/m3. (---_), predictions; (01, experimental; t---j, Walters eta/.’ Grid used 37 X 22

Notation

N

‘12

Figure 8 Displacement of centroid of two-dimensional bubbles as function of dimensionless time. pi = 1000 kg/m3. (Key as in Figure 7)

tions have been done for a plane flow situation on a grid 37 x 22. The computed results are displaced by a certain value from the experimental data. Waters and Davidson have noted a zero error of 0.098 on the abscissa of Figures 7 and 8 for the experimental data. If the plot from computed values is displaced on the time axis a good fit with the experimental points is obtained.

28

Appl.

Math.

Modelling,

1981,

Vol 5, February

Limited, are thanked code.

; g k

P r s V W

Z P IJ PL Pg

Radius of spherical bubble Vertical diameter of bubble Acceleration due to gravity Ratio of specific heats of gas Pressure Radial coordinate Displacement of bubble centroid Radial velocity Velocity in z-direction Axial coordinate Density Viscosity Liquid density Gas density