- Email: [email protected]
CFD Simulation of Mixing and Dispersion in Bubble Columns
0263–8762/03/$23.50+0.00 # Institution of Chemical Engineers Trans IChemE, Vol 81, Part A, September 2003
www.ingentaselect.com=titles=02638762.htm
CFD SIMULATION OF MIXING AND DISPERSION IN BUBBLE COLUMNS K. EKAMBARA and J. B. JOSHI Institute of Chemical Technology, University of Mumbai, Mumbai, India
L
iquid phase mixing time was measured in 0.2 and 0.4 m i.d. columns over a wide range of super cial gas velocity (0.07–0.295 m s 1) and height-to-diameter ratio (1–10). A CFD model was developed for the prediction of ow pattern in terms of mean velocity and eddy diffusivity pro les. The predictions agree favourably with all the experimental data published in the literature. Complete energy balance was established in all cases. The validated CFD model was extended for the simulation of the macro-mixing process by incorporating the effects of both the bulk motion and the eddy diffusion. Excellent agreement was observed between the CFD predicted and experimental values of the mixing time over the entire range of D, VG, and HD=D covered. The model was further extended for the prediction of residence time distribution and hence the axial dispersion coef cient (DL). The predicted values of DL were found to agree very well with all the experimental data published in the literature. Keywords: bubble column; gas–liquid; mixing; residence time distribution; computational uid dynamics; axial dispersion coef cient.
the fact that it contains a single parameter. Many correlations have been presented for the axial dispersion coef cient such as those by Reith et al. (1968), Ohki and Inoue (1970), Towell and Ackermann (1972), Kato et al. (1972), Hikita and Kikukawa (1974), Deckwer et al. (1974), Shah et al. (1978), Joshi (1980), Field and Davidson (1980), Riquarts (1981), Krishna et al. (2000) and Moustiri et al. (2001). The results show that the axial dispersion coef cient depends very strongly on column diameter (D) and to a lesser extent on the super cial gas velocity (VG). Some publications also mention a slight in uence of kinematic viscosity n of the liquid phase. The following is a typical form of the correlation:
INTRODUCTION Bubble column reactors nd widespread application in industrial processes such as absorption, catalytic slurry reactions, bio-reactions, coal liquefaction and waste water treatment. Bubble columns offer many advantages as gas– liquid contactors due to their simple construction, low operating cost and the ease with which the liquid residence time can be varied. However, the scale-up of bubble columns is still poorly understood because of the complexity of ow patterns, and their unknown behaviour under different sets of design parameters such as diameter, heightto-diameter ratio, gas sparger type, etc. The knowledge of mixing and axial mixing in the liquid phase is essential for the modelling, design and optimization of bubble column reactors. Accordingly, there have been several experimental studies and the results have been presented in the form of empirical correlations. The important process parameter for the batch mixing process is the mixing time. Experimental determination of the mixing time involves introduction of a tracer at some location in the column (usually a pulse input) and measuring the tracer concentration as a function of time with the help of a sensor. Mixing time is typically considered as the time required for the tracer concentration to reach the desired degree of homogeneity. In a continuous operation of bubble columns, the liquid phase axial mixing and the axial dispersion model is commonly used. The liquid is assumed to be in plug ow with axial dispersion superimposed on it. This model lumps the different mechanisms of liquid mixing into a single axial dispersion coef cient. The popularity of the model is due to
DL ˆ nDb VGc nd
(1)
Table 1 gives the summary of experimental measurements and correlation of axial dispersion coef cient. The value of the power law index b varies between 1 and 1.5, that on the super cial gas velocity ranges between 0.25 and 0.5, while viscosity exerts such a slight in uence that it is either neglected altogether or quoted with a very small exponent. The above-mentioned correlations are semi-empirical and it is most desirable to develop predictive procedures on the basis of understanding of the hydrodynamics. An initial attempt in this direction has shown that the dispersion coef cient is proportional to the average liquid circulation velocity (VC; Joshi and Sharma, 1979; Field and Davidson, 1980), the value of VC was estimated using energy balance. The mixing process occurs due to the transport at three levels: molecular, eddy and bulk (convection). Usually, the bulk motion (or bulk diffusion) is superimposed on either molecular or eddy diffusion or both. The convective 987
DL ˆ 0:31DVL (0) DL ˆ VL HD =Pe(1 ¡ EG )
DL ˆ 0:068g0:375 D1:5 VG0:375 n¡0:125
0.52–5.5 4.25, 4.50 150, 200
Aqueous NaCl solution
Pulse Aqueous NaCl solution 4.0 174, 380, 630
@ 1@ @ (EL c) ‡ (rvEL c) ‡ (uEL c) r @r @t @z´ ³ ³ ´ 1@ @ @ @ ˆ rDm (EL c) ‡ Dm (EL c) r @r @r @z @z
(2)
c ˆ c· ‡ c0
v ˆ v· ‡ v0 u ˆ u· ‡ u0 EL ˆ E· L ‡ E0L
(3)
Substituting the above terms in equation (2) and incorporating the Reynolds averaging procedure (further details about solution have been given Joshi, 2001), we get: @ 1@ @ (·E c· ) ‡ (r·EL c· v· ‡ E· L c0v0 ) ‡ (·EL c· u· ‡ E· L c0 u0) r @r @t L @z ³ ´ ³ ´ 1@ @ @ @ ˆ rDm (·EL c· ) ‡ Dm (·EL c· ) (4) r @r @r @z @z @ 1@ @ (·E c· ) ‡ (r·EL c· v· ) ‡ (·EL c· u· ) r @r @t L @z ³ ´ ³ ´ 1@ @ @ @ ˆ rDm (·EL c· ) ‡ Dm (·EL c· ) r @r @r @z @z 1@ @ ¡ (r·E c0 v0) ¡ (·EL c0 u0) (5) r @r L @z The correlation of uctuating concentration and velocity can be modelled using Boussinesq’s approach: µ ¶ @·c c0 v0 ˆ ¡Dtr @r µ ¶ (6) @·c c0 u0 ˆ ¡Dtz @z substitution of equation (6) in equation (5) and subsequent rearrangement gives:
Moustiri et al. (2001)
Krishna et al. (2000)
5.0–35.0
Pulse
— — — — — Riquarts (1981)
Field and Davidson (1980)
4.0–5.5 18.9 3200
Radioactive
UV recorder
—
DL ˆ 0:31DVC
DL ˆ 0:678D1:4 VG0:3 DL ˆ 0:33DVC
DL ˆ 0:66D1:25 VG0:38 Transient state measurement
Steady state Electrolyte, dye and heat —
Aqueous KCl solution
7.23 — 200 —
4.3–33.8 1.5, 2.4 100, 190
Deckwer et al. (1974) Joshi (1980)
Hikita and Kikukawa (1974)
5.0–12.0 —
DL ˆ nD1:4 VG0:3 DL ˆ nDb VGc nd Steady state, transient pulse
Pulse Aqueous KCl solution 2.0, 4.05
Methylene blue dye 1.7–26.8 1.5, 1.9, 2.84 406, 1067
6.6, 12.2, 21.4 Kato et al. (1972)
Towell and Ackermann (1972)
4.0–18.0
Steady state Transient state measurement 7.5–45.0 2.0–25.0 3.53, 3.8 2.0, 3.0 Reith et al. (1968) Ohki and Inoue (1970)
14, 29 40, 80, 160
Aqueous NaCl solution Potassium chloride
Technique of measurement Tracer Column height (m)
Range of VG (cm s¡1)
transport can be characterized by the mean velocity components in the liquid phase. The dispersive transport can be characterized by eddy diffusivity in the liquid phase. In view of the controlling role of the convection and eddy diffusion, the solution of convection-diffusion equation for the liquid phase can give insight into the mechanism of axial mixing and hence a rational expression for the effective axial dispersion coef cient. The transient mass balance for a tracer substance in a two-dimensional axisymmetric cylindrical co-ordinate system is given by the following equation:
The instantaneous concentration (c) may be written as the sum of the time-averaged concentration c· and a concentration uctuation c0. Similarly, the velocity components and hold-up can also be written:
Column diameter (mm) Author
Table 1. Axial dispersion coef cient: experimental details and correlation of previous work.
DL ˆ 0:30D2 VG1:2 ‡ 170d DL ˆ 0:30D2 VG1:2 ‡ 170d
EKAMBARA and JOSHI
Correlation
988
@ 1@ @ (·E c· ) ‡ (·EL c· v· ) ‡ (·EL c· u· ) r @r @t L @z µ ¶ 1@ @ ˆ r(Dm ‡ Dtr ) (·EL c· ) r @r @r µ ¶ @ @ ‡ (Dm ‡ Dtz ) (·EL c· ) @z @z
(7)
LDA measurements in bubble column have shown that (Deshpande et al., 2000) the rms velocities in the radial and axial direction are practically equal with deviation at some location within 0.56%. In view of this the sensitivity of Trans IChemE, Vol 81, Part A, September 2003
CFD SIMULATION OF MIXING AND DISPERSION eddy diffusivity in axial (Dtz) and radial (Dtr) dimensions was analysed and it will be shown later that the variation in Dtr and Dtz (within 0.63%) has a nominal effect. In view of this it will be assumed that Dtr ˆ Dtz ˆ Dt. With these simpli cations, equation (7) takes the following form: @ 1@ @ (·EL c· ) ‡ (r·vEL c· ) ‡ (u·EL c· ) r @r @t @z¶ µ µ ¶ 1@ @ @ @ ˆ rDeff (·EL c· ) ‡ Deff (·EL c· ) r @r @r @z @z
MODEL FORMULATION Equation of Continuity and Motion and the k±e Model For the case of bubbly ow, the equations of continuity and motion for a two-dimensional cylindrical coordinate system can be represented in the following generalized form: 1@ @ (rErvF)K ‡ (EruF)K r @r @z ³ ´ ³ ´ 1@ @F @ @F ˆ rEG ‡ EG ‡S r @r @r K @z @z K F,K
(8)
where Deff denotes the effective (molecular and eddy) diffusion coef cient. The above equation assumes that (i) the solute does not undergo any chemical reaction and (ii) all of the physical properties of the system remain constant. Thakre et al. (1999) have used computational uid dynamics (CFD) to obtain the ow pattern (radial and axial pro les of eddy diffusivity and axial, radial mean velocities) and used the ow pattern for the prediction of mixing time by solving the above equation. They compared the predictions with the experimental data from a 0.2 m i.d. column with height to diameter ratio of 8 and found excellent agreement. In the present work, the CFD model has been extended for the prediction of axial dispersion coef cient over a wide range of D, VG and HD=D. A comprehensive comparison has also been presented with the experimental results published in the literature.
989
(9)
where F is the transport variable (for instance, F is one for equation of continuity or it is a velocity component for the respective equation of motion, and F ˆ k and e for the conservation equation for the turbulent kinetic energy and turbulent energy dissipation rate respectively). K denotes the phase (K ˆ G or L) and SF is the source term for the respective dependent variable. The values of F and SF for different transport variables have been given in Table 2. From Table 2, it can be seen that most of the terms are derived from gas and liquid velocities. In addition, in twophase ows, momentum and energy transfer occurs across the interface. Therefore, the drag force (FDZ and FDR), virtual mass force (Fvz and Fvr), and lift force (FLr) appear in the axial and radial components of the momentum balance (Table 2). Furthermore, the interphase transfer of
Table 2. Two-dimensional k–e model for gas–liquid ows governing equations. The governing equations written in a general form: ³ ´ ³ ´ 1 @ @ 1 @ @F @ @F EG (rErvF)K ‡ (EruF)K ˆ rEG ‡ ‡S r @r @z r @r @r K @z @z K F,K Conservation of:
F
sF
sf
Mass
1
1
1 to 1
Axial velocity momentum
u
1.0
1 to 1
Radial velocity momentum
v
1.0
1 to 1
Turbulent kinetic energy
k
1.0
—
Turbulent dissipation energy
e
1.3
SF,K ˆ source terms
³ ´ ³ ´ 1 @ mt,K @EK @ mt,K @EK r ‡ r @r sf @r @z sf @z µ ³ ´ ³ ´¶ @P 1 @ @v @ @u ¤ ¤ EL § FVS EL ‡ Emt ¡EK ‡ EK rK g § FDZ rEmt ‡ @z r @r @z @z @z K µ ³ ´ ³ ´¶ ³ ´ ³ ´ mt @E 1 @ mt @E @ mt @E 1 @ @u ‡ uK r ‡ ‡ (rv) ‡ sf @z K r @r r @r sf @r @z sf @z K @z K µ ³ ´ ³ ´¶ @P 1 @ @v @ @u ¤ ¤ ¤ EL § FVS EL § FLr EL ‡ Emt ¡EK § FDr rEmt ‡ @r r @r @z @z @r K µ ³ ´ ³ ´¶ ³ ´ µ ¶ ± ² mt @E v 1 @ mt @E @ mt @E 1 @ @u ¡ 2Emt 2 ‡vK r ‡ ‡ (rv) ‡ r K r @r sf @r @z sf @z K sf @r K r @r @z K
—
EL (G ‡ PB ¡ rL e) e EL [Ce1 (G ‡ PB ) ¡ Ce2 rL e] k @u ˆ ¡CL EG rL (uG ¡ uL ) L @r
mt,K ˆ 0:09rk (k 2 =e); PB ˆ CB (FDr VSr ‡ FDz VS ); FLr » ¼ @ @ Fvz ˆ ¡CV EG rL (u ¡ uL ) ‡ (uG ¡ uL ) ; FDz ˆ ¡4:9 £ 104 EG (uG ¡ uL ) @r G @z µ ¶ 1 @ @ Fvr ˆ ¡CV EG rL r(vG ¡ vL ) ‡ (vG ¡ vL ) ; FDr ˆ ¡4:9 £ 104 EG (vG ¡ vL ) r @r @z " Á» ¼ ¼ 2 ! ³» ¼ » ¼ ´2 # 2 n o2 » @vL v @uL @vL @uL G ˆ mt,L 2 ‡ L ‡ ‡ ‡ @r r @z @z @r Force terms are positive for the gas phase and negative for the liquid phase.
Trans IChemE, Vol 81, Part A, September 2003
990
EKAMBARA and JOSHI
energy (PB) appears in the equations for turbulent kinetic energy and turbulent energy dissipation rates.
eD ˆ vB (rL ¡ r G)gVS
Boundary and Initial Conditions Boundary and initial conditions for equation (9) (i) Along the axis: axisymmetry in uG, uL, EG, vG, vL, k and e. (ii) Along the wall: the velocities satisfy the no-slip boundary conditions (the wall function method based on the log law of the wall is invoked to calculate the wall shear stress and the values of k and e close to the wall). (iii) At the inlet: gradients of vG, vL, k and e are set to zero. (iv) At the top surface of the computational domain, the gradient of the dependent variable is set to zero. For initiating the numerical solution, EG and uG are speci ed at the inlet. At all the other locations uG, uL, vG, vL and EG are taken to be zero at t ˆ 0. For k and e, the initial guess values were found to be important. Boundary and initial conditions for equation (8) The solution to the diffusion-convection equation (8) must satisfy the following conditions: c· (r, z, 0) ˆ 0 c· (r, 0, 0) ˆ Ci @·c(0, z, t) @·c(R, z, t) ˆ ˆ0 @r @r
for all z and t
Energy Balance Consider a single bubble (of volume vB) rising in a pool of liquid. The bubble at any location has pressure energy and potential energy. As the bubble rises, the pressure energy decreases and the potential energy increases. The rate of change of energy is given by the following equation: dh ˆ v(rL ¡ rG )gVS ˆ fD VS (10) dt where VS is the slip velocity and fD is the drag force. When the bubble rise is in a creeping ow regime, the energy released during rise is dissipated in a viscous manner. In other words, the gas phase energy gets directly converted into internal energy. When the bubble rise is in the turbulent ow regime (Re1 > 500), the energy released during the rise is rst of all converted to turbulence in the liquid phase and nally into internal energy at the Kolmogorov scale. Therefore, the fraction of gas phase energy used for generating liquid phase turbulence depends upon the bubble Reynolds number. There is one more aspect of the energy balance. In bubble columns, the bubble rise velocity with respect to the wall can be very different from the slip velocity because of liquid circulation. When the actual velocity is higher than the slip velocity (as in the central upward region), the rate of energy release is different from that given by equation (10). vB (rL ¡ rG )g
eR ˆ vB (r L ¡ rG)g(VS ‡ uL )
where uL is the liquid velocity at the location of bubble rise. The rate of energy converted into turbulence and spent in the viscous dissipation depends only on the slip velocity:
(11)
(12)
Therefore, when uL is positive, the rate of energy release is higher than the dissipation. The extra energy is used for maintaining the liquid circulation. In contrast, when the bubble moves downward in the near-wall region, the liquid phase energy is supplied to bubbles. The overall effect of up ow and down ow (because of the gas hold-up pro le, more bubbles are present in the up ow region as compared to the down ow region) is a net energy supply from gas phase to liquid phase. The rate of energy release from the gas phase to the liquid phase is given by the following equation: p E ˆ D2 (rL ¡ rG )gHD E·L (VG ¡ E·G VS ) (13) 4 The above equation assumes that the slip velocity is same for all the bubbles. The rate of energy converted into turbulence and spent in viscous energy dissipation is given by: p ED ˆ D2 (rL ¡ rG )gHD E·G E· L VS (14) 4 It may be, however, noted that the turbulence generated by bubbles has a scale corresponding to bubbles, whereas the turbulence generated in liquid circulation [using the energy given by equation (13)] has a scale corresponding to column dimensions. Therefore, this gives rise to a very important question regarding the quantitative role of the bubble generated turbulence to the transport phenomena. Sato and Sekoguchi (1975) and Theofanous and Sullivan (1982) have addressed to this problem. It was assumed that a fraction of CB is used in the transport and this fraction forms a source terms for the liquid phase k (turbulent kinetic energy per unit mass) equation: p ES ˆ CB D2 (rL ¡ r G)gHD E· L E· G VS (15) 4 where CB takes values between 0 and 1. Johansen and Boysan (1988) have included the radial variation of slip velocity and CB was in the range of 0.1–0.2. Later Svendsen et al. (1992) used CB as a tuning parameter to match their experimental results. Kataoka et al. (1992) reported a detailed analysis of the extra terms due to the presence of the dispersed phase in the source of turbulent kinetic energy. Their analysis suggests that the extra generation of turbulence due to large bubbles is almost compensated by the extra dissipation due to the small-scale interfacial structures. The studies by Hillmer et al. (1994) and Ranade (1992, 1997) also suggest zero value for CB. Hjertager and Morud (1993) have used a very low value of 0.02. The total energy released from the gas phase to the liquid phase is given by equation (13). Furthermore, the quantitative role of bubble-generated turbulence is given by equation (15). Therefore, the sum of equations (13) and (15) gives the total energy received by the liquid phase and is expressed as: p E ‡ ES ˆ D2 (rL ¡ rG )gHD E·L [VG ‡ (CB ¡ 1)·EG VS ] 4 (16) Trans IChemE, Vol 81, Part A, September 2003
CFD SIMULATION OF MIXING AND DISPERSION It must be emphasized at this stage that, whatever may be the value of CB, the energy balance must be satis ed. When the k–e model is used for the prediction of ow pattern, we get radial and axial variation of e (energy dissipation rate per unit mass) as one of the answers. From this e eld, the total energy dissipation rate can be calculated by suitable volume integration. The total energy dissipation rate must equal the energy-input rate given by equation (16). METHOD OF SOLUTION The mathematical model was solved in two steps. In the rst step, equations of continuity and motion were solved (together with k–e equations) for getting complete ow pattern in terms of distribution of gas and liquid velocities, eddy diffusivity and gas holdup. This ow pattern was then used in step two, for the estimation of mixing time and residence time distribution. The set of equations given in Table 2 was solved numerically, which consisted of the following steps: (i) the generation of suitable grid points; (ii) the conversion of governing equations into algebraic equations; (iii) the selection of discretization schemes; (iv) the formulation of discretzed equations at every grid location; (v) the formulation of pressure equation; and (vi) the development of suitable iteration scheme for use in obtaining the nal solution. The nite control volume technique of Patankar (1981) was employed for the solution of equations given in Table 2. A staggered grid arrangement proposed by Patankar and Spalding (1972) was employed. It consisted of 50 £ 100 grid points with 50 grid points in the radial direction and 100 grid points in axial direction. The velocity components were calculated for the points that lie on the faces of the control volume, while all other variables were calculated at the centre of the control volume. Several discretization schemes (upwind, hybrid, exponential and the power law) were tested. Out of these, the power law scheme was found to be the best and was used for the discretization of the governing equations. The most important numerical concern is the manner in which the pressure eld is computed and its effect on the continuity of each eld. For multidimensional analyses of single phase and homogenous two phase ows, many algorithms like SIMPLE, SIMPLER (Patankar, 1981), SIMPLEST (Spalding, 1980) and SIMPLEC (Van Doormal and Raithby, 1984) have been proposed for solving the pressure velocity coupling. All the above-mentioned algorithms incorporate the pressure equation, thus using the staggered grid concept, but differ principally in the exact form of the pressure equation, the manner in which the equations are formulated over the control volume, the degree of implicitness used in time advancement, and the solution sequence (Craver, 1984). Although other numerical schemes have been found to give faster convergence than SIMPLE, one iteration of these schemes involves more computational effort than SIMPLE (Patankar, 1981). Moreover, as compared with other schemes, the SIMPLE algorithm has been extensively used and has served well in a number of CFD simulation studies. In an iterative type of solution, it is a general practice to adopt some kind of convergence criteria. In fact, it is this which will decide not only the accuracy but also the Trans IChemE, Vol 81, Part A, September 2003
991
Table 3. Effect of internal iterations on the rate of convergence. Pressure
vr
vz
ur
uz
k
e
Iterations
3 4 5 6 10
2 2 2 4 6
2 2 2 4 6
2 3 1 4 6
2 3 1 4 6
10 4 5 8 20
10 6 5 10 20
624 624 624 624 624
computational time required and stability of the given problem. However, if convergence criteria are given without prior knowledge of the numerical problem, it often ends up in using excessive computational time. Hence, in the present simulation the role of the convergence criteria on the ow variable was critically analysed. Effect of Internal Iterations For an iterative process, it is a common practice that two or three iterations are used for solving the set of algebraic equations for vr, vz, ur, uz, k and e. The number of iterations used for pressure correction equation is typically four to ve times than that for the above-mentioned variables when the line-by-line method TDMA method is used. However, the above-mentioned practice does not give better convergence. The coef cients of the velocity variables are in uenced by the eddy viscosity. Therefore, the solution of the k and e equations should be iterated more times than the pressure correction equation. However, as seen from Table 3, the variation of number of internal iterations for k and e as well as other variable gives the same rate of convergence of the solution (i.e. 624 iterations), indicating its independence of the interal iterations. Effect of the Under-relaxation Parameter It has been shown by Peric et al. (1988) that a nearly optimum convergence is found for the SIMPLE algorithm when the under-relaxation parameters for pressure and velocity variables are selected by the following relation ap ˆ 1 7 au. In addition, in the present study it has been observed that the under-relaxation parameters for k and e are also imporatant to obtain a converged solution. It was found that, by decreasing the under-relaxation parameter for pressure correction equation, ap from 0.75 to 0.1, the total number of iterations required for convergence increased by 18%. A similar observation was found with the axial liquid velocity relaxation equation. However, no improvement in the rate of convergence was found by increasing the under-relaxation parameters of the radial and axial gas phase momentum equations and turbulent energy dissipation rate equation. In Table 4, the effect of under-relaxation Table 4. Effect of under relaxation parameters on the rate of convergence. ap 0.75 0.10 0.75 0.75 0.75
avr
avz
aur
auz
ak
ae
Iterations
0.75 0.75 0.75 0.75 0.75
0.75 0.75 0.75 0.75 0.75
0.75 0.75 0.75 0.5 0.75
0.75 0.75 0.5 0.75 0.75
0.75 0.75 0.75 0.75 0.75
0.75 0.75 0.75 0.75 0.95
624 818 630 620 567
0.002 0.002 0.002 0.002 0.002 0.002 0.002 0.002 0.012 0.2 0.2
684 684 684 684 684 684 684 671 634 616 622
0.6=3.45
Five-point conductivity probe Height measurement 3.50 20–80 Air–water
0.475 10–140
Hole diameter, 0.2 mm Perforated pipes: 6 mm, 8 nos, hole diameter 1 mm, 78 nos. Hole diameter 0.2 mm
Air–deionized water Air–water
20–80
3.49
Electroconductivity microprobe Electroconductivity microprobe Two channel optical bre probe 4.86 12–96 Air–water
Optical probe — 53–1452 Air–water
Liquid velocity
Pavlov tube Electroconductivity needle probe 4.34 19–169 Air–water
Holdup Sparger
Seive plate: Hole diameter 0.4 mm, 61 nos Seive plate hole diameter, 1 mm —
Measurement location HD=D
Fly-wheel anemometer Hot lm anemometry Hot lm anemometry Two Pt–Rh electrodes, two U shaped optical bres and a tracer injection
0.29=4.50
0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.005 0.015 0.015 0.15
Grienberger and Hofmann (1992)
0.005 0.005 0.005 0.005 0.005 0.005 0.05 0.005 0.005 0.005 0.005
0.254=2.50
¡0.02 ¡0.04 ¡0.5 0.5 ¡0.02 ¡0.02 ¡0.02 ¡0.02 ¡0.02 ¡0.02 ¡0.02
Yu and Kim (1991)
0.33 0.33 0.33 0.53 0.33 0.33 0.33 0.33 0.33 0.33 0.33
0.29=4.50
Iterations
Yao et al. (1991)
e
Menzel et al. (1990)
k
0.45=5.60
uz
Nottenkamper et al. (1983)
vz
0.138=1.37
0.02 0.02 0.02 0.02 0.02 0.05 0.02 0.02 0.02 0.02 0.02
ur
Hills (1974)
vr
Researchers
Table 5. Effect of initial guesses of the ow variables on the rate of convergence for a given set of under-relaxation parameters.
Super cial gas velocity (mm s¡1)
The effect of the initial guess values on the rate of convergence is tabulated in Table 5. It is interesting to note that the initial guess values of vr, vz, ur, uz and p have no in uence on that rate of convergence. The rate of convergence of the iterative process primarily depends on the guess values of k and e. However, it is dif cult to establish any mathematical relation between the guess values of k and e that yield an optimum rate of convergence. In the present work, we have investigated the performance of SIMPLE, SIMPLER and SIMPLEC, and SIMPLE was found to be the best. The set of algebraic equations obtained after discretization was solved by TDMA. The relaxation parameters, initial guess values and internal iterations for the variables were tuned to optimize the balance between the convergence criteria (1.0 £ 10¡3) and the number of iterations required. While the solution was found to be independent of the initial guesses and internal iterations (for vr, vz, ur, uz and p), the relaxation parameters had a dominant effect on the convergence of solution. The effect was prominent for lower height to diameter ratios (HD=D) and higher super cial gas velocities (VG), enforcing the use of lower relaxation parameters. Further details about method of solution have been given Joshi (2001). After getting the ow pattern, the model was extended to the understanding of the mixing process. It has been assumed that the tracer concentration has no in uence on the ow characteristics. Therefore, ow results computed under fully developed and steady-state conditions were used to solve the transient scalar mixing problem. It was assumed that the eddy diffusivities for momentum and mass are equal. The transport equation for concentration in a cylindrical coordinate system is given in equation (8). A separate code was developed to simulate transient mixing in
Table 6. Comparison of radial gas hold-up and mean axial liquid velocity pro les: experimental data of previous work.
Effect of Guess Values of the Variables on the Rate of Convergence
Gas–liquid system
parameters of different variable has been given in terms of number of iterations. It can be seen from the above table that an optimum convergence can be achieved for the values of ak ˆ 0.5 and ae ˆ 0.95 with other parameters taking a constant value of 0.75. It may be emphasized that the above optimum combination is for a particular set of VG, D and HD=D. For another set the optimum of a is also different and the selection gets facilitated by the experience in CFD.
Measurement techniques employed
EKAMBARA and JOSHI
Column diameter= height (m)
992
Trans IChemE, Vol 81, Part A, September 2003
CFD SIMULATION OF MIXING AND DISPERSION bubble columns. This code uses the ow patterns results generated in the rst step. In the present case of scalar mixing, there is no source term.
993
to extend the model for the prediction of mixing time and the axial dispersion coef cient. Comparison Between the Predicted and Experimental Data for Liquid Phase Mixing Time
RESULT AND DISCUSSION Comparison Between the Predicted and Experimental Data for Flow Pattern As a rst step, it is important to establish the validity of the model for ow pattern. Therefore, comparison has been made with experimental data of Hills (1974), Nottenkamper et al. (1983), Menzel et al. (1990), Yao et al. (1991), Yu and Kim (1991) and Grienberger and Hofmann (1992). Table 6 gives the details pertaining to column diameter, column height, sparger design, gas–liquid system, range of VG and VL, the measurement location and the measurement techniques employed for the measurement of mean axial liquid velocity and fractional holdup. It may be noted that practically all the major experimental investigations published in the literature have been covered in the table. The results are given in Table 7 and Figure 1A–F. The agreement between the predicted and the experimental pro les can be seen to be excellent. The agreement over the wide range of D, VG and the gas–liquid systems ensures the applicability of the model (for the estimation of ow pattern). In view of the abovementioned favourable comparison, it was thought desirable
The knowledge of ow pattern was extended for the prediction of liquid phase mixing time (ymix). A systematic study was undertaken to investigate the effect of various parameters such as the locations of tracer addition and detection, the quantity of tracer, the effect of time step, column diameter, super cial gas velocity (VG) and heightto-diameter ratio on the mixing time. In order to understand the role of detector position, the detectors were placed at the grid locations of 3 £ 100, 3 £ 50, 3 £ 3, 25 £ 60, 50 £ 50 and 50 £ 3. For all these cases, the tracer was introduced at the 50 £ 90 location. The results of concentration vs time pro les are shown in Figure 2. It can be seen that all the pro les originate from zero and reach a constant value of 0.162, the nal concentration. The concentration pro le at 3 £ 100 detector (close to the tracer input location) position shows a peak overshoot. For all the other detectors, the response was overdamped. However, the time required for homogenization can be seen to be independent of detector position. In second step, the effect of the location of tracer introduction was investigated. For this purpose, tracer was introduced at 50 £ 90, 3 £ 100, 3 £ 50, 3 £ 3, 25 £ 50 and 50 £ 3
Table 7. Comparison between CFD predictions and the experimental observations. Material balance Authors
¡1
Energy balance
VG (m s )
E·G
CB
C0
C1
LHS
RHS
0.0206 (0.019) 0.038 (0.038) 0.067 (0.064) 0.106 (0.095) 0.171 (0.169)
0.065 (0.067) 0.106 (0.107) 0.161 (0.160) 0.197 (0.194) 0.218 (0.212)
0.43 0.44 0.37 0.28 0.34
3.84 (2.0) 3.34 (2.0) 2.58 (2.0) 2.10 (2.0) 1.91 (2.0)
0.222 (0.25) 0.24 (0.25) 0.25 (0.25) 0.32 (0.25) 0.465 (0.25)
0.105 0.127 0.210 0.355 0.857
0.115 0.205 0.221 0.366 0.831
Nottenkamper et al. (1983)
0.054 (0.053) 0.106 (0.105) 0.331 (0.324) 0.601 (0.823) 1.209 (1.452)
0.121 (0.120) 0.175 (0.168) 0.312 (0.271) 0.405 (0.409) 0.494 (0.497)
0.36 0.38 0.41 0.21 0.23
2.64 (2.0) 2.5 (2.0) 1.64 (2.0) 1.63 (2.0) 1.61 (2.0)
0.375 (0.40) 0.35 (0.40) 0.54 (0.40) 0.39 (0.40) 0.40 (0.40)
0.190 0.502 1.815 2.38 6.59
0.198 0.513 1.836 2.29 5.80
Menzel et al. (1990)
0.012 (0.012) 0.026 (0.024) 0.052 (0.048) 0.096 (0.096) 0.021 (0.02) 0.041 (0.04) 0.062 (0.06) 0.081 (0.08)
0.029 (0.029) 0.049 (0.048) 0.105 (0.103) 0.137 (0.128) 0.057 (0.057) 0.108 (0.106) 0.156 (0.153) 0.186 (0.180)
0.36 0.33 0.32 0.44 0.37 0.44 0.43 0.46
3.13 (2.5) 5.80 (2.5) 5.30 (2.5) 3.02 (2.5) 2.96 (2.0) 2.71 (2.0) 2.28 (2.0) 2.69 (2.0)
0.39 (0.40) 0.36 (0.40) 0.42 (0.40) 0.35 (0.40) 0.285 (0.25) 0.236 (0.25) 0.251 (0.25) 0.223 (0.25)
0.057 0.130 0.223 0.593 0.075 0.208 0.299 0.417
0.07 0.157 0.263 0.588 0.078 0.233 0.312 0.447
Yu and Kim (1991)
0.011 (0.01) 0.02 (0.02) 0.040 (0.04) 0.061 (0.06) 0.101 (0.1) 0.142 (0.14)
0.025 (0.025) 0.035 (0.034) 0.049 (0.048) 0.062 (0.060) 0.090 (0.084) 0.136 (0.130)
0.28 0.30 0.21 0.20 0.27 0.22
3.17 (2.1) 2.4 (2.1) 2.5 (2.1) 2.23 (2.1) 1.87 (2.1) 1.77 (2.1)
0.39 (0.40) 0.554 (0.73) 0.739 (0.73) 0.853 (0.73) 0.76 (0.73) 0.80 (0.73)
0.045 0.122 0.187 0.322 0.526 0.661
0.054 0.100 0.194 0.288 0.549 0.668
Grienberger and Hofmann (1992)
0.021 (0.02) 0.081 (0.08)
0.057 (0.054) 0.156 (0.150)
0.37 0.46
2.96 (2.0) 2.70 (2.0)
0.285 (0.25) 0.22 (0.25)
0.075 0.417
0.078 0.447
Khare (1990)
0.072 (0.07) 0.173 (0.17) 0.292 (0.295)
0.136 (0.135) 0.223 (0.223) 0.304 (0.282)
0.30 0.32 0.39
2.17 (3.01) 1.99 (3.01) 1.76 (3.01)
0.38 (0.30) 0.43 (0.30) 0.51 (0.30)
0.362 0.821 1.292
0.329 0.815 1.488
Hills (1974)
Yao et al. (1991)
The bracketed numbers indicate the experimental values; LHS is volume integral of energy dissipation rate and RHS is equation (16).
Trans IChemE, Vol 81, Part A, September 2003
994
EKAMBARA and JOSHI
locations. For all these cases, the detector position was at 50 £ 100. The results are given in Table 8. It can be observed that the mixing time is practically independent of the point of tracer introduction. The same result was found to hold for many other detector locations. The effect of tracer quantity on mixing time was also investigated and was found to have no effect on mixing time. These numerical results (for the effects of tracer quantity and locations of tracer introduction and detection) were also con rmed experimentally. The effect of height to diameter ratio (HD=D) on mixing time was studied in 0.2 and 0.4 m i.d. column for a gas–liquid system with variation in VG from 0.07 to 0.295 m s¡1, respectively. For 0.2 m i.d. column HD=D ratio was varied from 3 to 10 and from 1 to 5 for 0.4 m i.d. column. The results have been shown in Figure 3A–C for column diameter 0.2 m and super cial gas velocity from 0.07 to 0.295 m s¡1. The comparison between the predicted and experimental data (Khare, 1990) of mixing time for different column diameters, super cial gas velocities and height to diameter ratios are given in Table 9. It can be observed from Table 9 that, for a particular VG, the mixing time increases with an increase in HD=D ratio for both the columns. From Table 9, it can also be noted that, at the same HD=D ratio the mixing time increases with an increase in column diameter (D). Further, it can be seen that the mixing time decreases with an increase in VG. Comparison Between the Predicted and Experimental Data for Axial Dispersion Coef cient Residence time distribution curves were obtained from the CFD simulation. For this purpose, concentrations were estimated at the topmost axial location and at all the radial locations. The average cross-sectional concentration at the exit gives the mixing cup concentration and the values were estimated with respect to time. The resulting residence time distribution curves are shown in Figures 4 and 5 for 0.2 and 0.4 m i.d. column with variation in VG from 0.07– 0.295 m s¡1, respectively. For 0.2 m i.d. column HD=D ratio was varied from 3 to 10 and from 1 to 5 for 0.4 m i.d. column. The range of super cial liquid velocity has been given in Table 12. These curves were used for the estimation of axial dispersion coef cient following the stepwise procedure suggested by Levenspiel (1999). The mean solute concentration is given by: … 1 R C ˆ 2 2p·cr dr (17) pR 0 When the reactor is modelled as a closed vessel, the rst and second moments of the residence time distribution curve are related to the mean residence time (t), and the axial dispersion coef cient, DL, respectively (Levenspiel, 1999) as: P tC(t)Dt t P (18) C(t)Dt s2
P2 t C(t)Dt P ¡ t2 C(t)Dt
(19)
Figure 1. Comparison between the simulated and experimental pro les of fractional gas holdup and axial liquid velocity (refer to Table 3 for experimental details). (A) Hills (1974); (B) Nottenkamper et al. (1983); (C) Menzel et al. (1990); (D) Yao et al. (1991); (E) Yu and Kim (1991); (F) Grienberger and Hofmann (1992).
Trans IChemE, Vol 81, Part A, September 2003
CFD SIMULATION OF MIXING AND DISPERSION
Figure 1. Continued.
Trans IChemE, Vol 81, Part A, September 2003
995
996
EKAMBARA and JOSHI
s2y
³ ´µ ³ ´¶ s2 DL DL V L HD ˆ 2 ˆ2 ¡2 1 ¡ exp ¡ t VL HD VL HD DL (20) VL HD 2 sy 2
DL
Figure 2. Concentration pro les at various detector positions for VG 0.17 m s 1, D 0.2 m and HD=D 10.
(21)
where t is the mean residence time, DL is the axial liquid dispersion coef cient, HD is the height of the column and s represents the variance. The effect of column diameter (D) and super cial gas velocity (VG) on axial dispersion coef cient was studied. The comparison between the CFD predicted results and experimental data of different investigators are given in Table 10. It can be observed that the axial dispersion coef cient increases with an increase in VG with a proportionality in the range of VG0:25±0:42 . Further, it can also be
Table 8. Effect of tracer introduction positions on liquid phase mixing time. Liquid phase mixing time for various tracer introduction positions in the column (s) VG (m s¡1) 0.07 0.17 0.295
E· G
50 £ 90
3 £ 100
3 £ 50
3£3
25 £ 50
50 £ 3
0.135 0.223 0.282
71.89 67.81 56.63
71.89 67.81 56.63
71.89 67.81 56.63
71.89 67.81 56.63
71.89 67.81 56.63
71.89 67.81 56.63
0.07; (B) VG
0.17; and (C) VG 0.295m s 1; 1—3
Figure 3. Effect of HD=D ratio on mixing time for column diameter 0.20 m. (A) VG 3—25 60; 4—50 50; 5—3 3; 6—50 3.
100; 2—3
50;
Trans IChemE, Vol 81, Part A, September 2003
CFD SIMULATION OF MIXING AND DISPERSION
Figure 3. Continued.
Trans IChemE, Vol 81, Part A, September 2003
997
998
EKAMBARA and JOSHI Table 9. Mixing time (s). HD=D
Column diameter ˆ 0.2 m VG (m s¡1) 0.070 0.170 0.295
3
8
10
E· G
Expt ymix
CFD ymix
Expt ymix
CFD ymix
Expt ymix
CFD ymix
Expt ymix
CFD ymix
0.135 0.223 0.282
33.5 29.3 25.2
29.8 30.0 24.2
38.4 35.2 31.2
39.0 36.2 32.0
51.6 45.6 39.6
52.0 45.1 40.3
72.5 67.2 —
71.9 67.8 56.6
Column diameter ˆ 0.4 m 0.070 0.170 0.295
5
0.135 0.223 0.282
2 56.5 51.6 47.0
3 57.1 52.0 45.8
66.0 61.2 55.5
Figure 4. Effect of HD=D ratio on RTD for column diameter 0.20 m. (A) VG 0.07; (B) VG 0.17; and (C) VG 0.295 m s 1.
4 65.2 59.7 53.8
78.0 71.5 64.0
5 75.9 69.4 65.2
88.0 79.5 —
86.3 75.8 69.4
Figure 5. Effect of HD=D ratio on RTD for column diameter 0.40 m. (A) VG 0.07; (B) VG 0.17; and (C) VG 0.295 m s 1.
Trans IChemE, Vol 81, Part A, September 2003
CFD SIMULATION OF MIXING AND DISPERSION
999
Table 10. Comparison between CFD predictions and the experimental observations. Material balance Authors
VG (m s¡1)
Energy balance
E· G
C0
C1
LHS
RHS
Experimental DL (m2 s¡1)
Predicted DL (m2 s¡1)
Ohki and Inoue (1970)
0.02 (0.02) 0.048 (0.05) 0.098 (0.10) 0.148 (0.15) 0.21 (0.20)
0.063 (0.062) 0.123 (0.122) 0.205 (0.218) 0.301 (0.307) 0.380 (0.382)
4.223 3.699 2.568 2.042 1.923
0.226 0.238 0.262 0.257 0.213
0.103 0.249 0.485 0.629 0.928
0.101 0.247 0.497 0.638 0.957
0.0023 0.0056 0.0097 0.0123 0.0147
0.0031 0.0042 0.0106 0.01315 0.01856
Towell and Ackermann (1972)
0.021 (0.020) 0.023 (0.022) 0.09 (0.088)
0.053 (0.05) 0.06 (0.06) 0.178 (0.18)
1.958 1.754 1.548
0.347 0.365 0.391
0.086 0.091 0.351
0.091 0.106 0.364
0.0412 0.0670 0.0930
0.0489 0.0668 0.0977
Hikita and Kikukawa (1974)
0.069 (0.07) 0.098 (0.10) 0.158 (0.16) 0.22 (0.23) 0.33 (0.32) 0.048 (0.05) 0.074 (0.075) 0.11 (0.10) 0.119 (0.12)
0.15 (0.152) 0.176 (0.175) 0.20 (0.207) 0.264 (0.267) 0.32 (0.317) 0.141 (0.143) 0.166 (0.165) 0.179 (0.180) 0.21 (0.212)
2.462 2.224 2.064 1.983 1.864 2.546 2.212 2.043 1.998
0.257 0.274 0.367 0.401 0.423 0.236 0.267 0.358 0.361
0.377 0.562 0.876 1.168 1.569 0.239 0.387 0.474 0.534
0.367 0.549 0.902 1.089 1.468 0.251 0.394 0.451 0.548
0.0302 0.0352 0.0411 0.0468 0.0542 0.0312 0.0358 0.0387 0.0423
0.0286 0.0331 0.0413 0.0473 0.0518 0.0321 0.0360 0.0365 0.0415
Krishna et al. (2000), D ˆ 0.174 m
0.033 (0.034) 0.092 (0.09) 0.158 (0.16) 0.22 (0.23) 0.27 (0.27) 0.31 (0.30)
0.090 (0.087) 0.168 (0.165) 0.21 (0.218) 0.256 (0.250) 0.265 (0.262) 0.268 (0.270)
3.398 2.948 2.643 2.421 2.182 1.986
0.286 0.306 0.321 0.372 0.423 0.486
0.164 0.477 0.885 1.260 1.452 1.566
0.182 0.513 0.898 1.326 1.483 1.624
0.0252 0.0291 0.0356 0.0401 0.0443 0.0498
0.0245 0.0278 0.0291 0.0377 0.0393 0.0411
van Baten and Krishna (2001), D ˆ 0.63 m Moustiri et al. (2001)
0.089 (0.09) 0.21 (0.23) 0.32 (0.30) 0.032 (0.0319) 0.047 (0.0472) 0.057 (0.0567)
0.075 (0.089) 0.254 (0.265) 0.288 (0.295) 0.076 (0.076) 0.125 (0.126) 0.155 (0.154)
2.846 2.392 1.992 2.467 2.146 2.013
0.323 0.382 0.476 0.347 0.371 0.333
0.477 1.260 1.566 0.137 0.150 0.166
0.513 1.326 1.624 0.148 0.161 0.173
0.1746 0.2483 0.2743 0.0032 0.0062 0.0086
0.1283 0.1674 0.2043 0.0051 0.0059 0.0091
Deckwer et al. (1974)
The bracketed numbers indicate the experimental values; LHS is volume integral of energy dissipation rate and RHS is equation (16).
observed that, at a particular VG, the value of DL increases with an increase in column diameter with a proportionality in the range of D1.2–1.6. The validity of the CFD model for axial dispersion coef cient was established. Therefore, comparison has been made with experimental data of Ohki and Inoue (1970), Towell and Ackermann (1972), Hikita and Kikukawa (1974), Deckwer et al. (1974), Krishna et al. (2000), Moustiri et al. (2001) and Khare (1990). Table 1 gives the detail pertaining to column diameter, column height, gas– liquid system, range of VG, tracer and the measurement techniques employed for the measurement of axial dispersion coef cient. It may be noted that practically all the major experimental investigations published in the literature have been covered in the table. The results are given in Table 10 and Figure 6. The agreement between the predicted and the experimental data can be seen to be excellent. The agreement over the wide range of D, VG and the gas–liquid systems ensures the applicability of the model (for the estimation of axial dispersion coef cient). Now we would like to bring out one observation regarding the relative importance of convection and eddy diffusivity (nt) on axial dispersion coef cient (DL). Table 11 shows the values of DL and average nt in the column. It can be seen that, in all the cases of VG and D, the value of DL is at least 10 times higher than that of nt. This observation indicates the overriding contribution mean circulatory ow on axial mixing and mixing as compared with that of nt. The sensitivity analysis of the effects of mean ow and nt Trans IChemE, Vol 81, Part A, September 2003
(results not shown) con rms this observation. In fact, the axial dispersion model assumes the mean ow to be unidirectional ow with a at pro le with axial dispersion superimposed. The ideal situation is one directional at pro le with eddy diffusivity superimposed. In comparison with such a ideal situation, the real dispersed plug ow model accounts for any nonidealities from the at velocity pro le such as bypass, recycle, channeling, dead zones, etc.
Figure 6. Parity plot for axial dispersion coef cient. ( ) Ohki and Inoue (1970); ( ) Towell and Ackermann (1972); ( ) Hikita and Kikukawa (1974); ( ) Deckwer et al. (1974); ( ) Khare (1990); ( ) Krishna et al. (2000); ( ) Moustiri et al. (2001); and ( ) van Baten and Krishna (2001).
1000
EKAMBARA and JOSHI Table 11. Comparison between CFD predictions and the experimental observation (Khare, 1990).
VG (m s¡1)
Material balance
Energy balance
Column diameter (m)
E· G
C0
C1
LHS
RHS
nt (m2 s¡1)
Dispersion coef cient (DL)
0.20 0.20 0.20 0.40 0.40 0.40
0.136 (0.135) 0.223 (0.223) 0.304 (0.282) 0.133 (0.135) 0.227 (0.223) 0.311 (0.282)
2.17 1.99 1.76 2.155 2.135 1.814
0.38 0.43 0.51 0.367 0.399 0.431
0.362 0.821 1.292 0.399 0.877 1.374
0.329 0.815 1.488 0.329 0.856 1.487
0.0018 0.0026 0.0032 0.0048 0.0067 0.0074
0.0287 0.0386 0.0469 0.0660 0.0956 0.1196
0.072 (0.07) 0.173 (0.17) 0.292 (0.295) 0.073 (0.07) 0.173 (0.17) 0.302 (0.295)
The bracketed numbers indicate the experimental values; LHS is volume integral of energy dissipation rate and RHS is equation (16).
Table 12. Axial dispersion coef cient (m2 s 1) VG (m s¡1)
VL (m s¡1)
HD=D
E· G
Column diameter ˆ 0.2 m 0.070 0.170 0.295
0.0066 0.015 0.035
0.135 0.223 0.282
Column diameter ˆ 0.4 m 0.070 0.170 0.295
0.0066 0.015 0.035
0.135 0.223 0.282
However, such nonidealities are expected to be within limits for the applicability of dispersed plug ow model. In fact, the average recirculation ow (recycle) in bubble columns is 10–1000 times higher than the mean net ow across the column. Such a major deviation is an alerting signal for checking the applicability of dispersed plug ow model for bubble columns. Let us review some of the characteristic features of the dispersed plug ow model visa-vis the performance of a bubble column. (i) For a dispersed plug ow model, the relationship between the mixing time, column height and DL is given by the following equation (Towell and Ackermann, 1972): ymix ˆ
aHD2 DL
(22)
3
5
8
10
0.0034 0.0139 0.0209
0.0125 0.0225 0.0340
0.0198 0.0316 0.0423
0.0287 0.0386 0.0469
2
3
4
5
0.0166 0.0345 0.0366
0.0260 0.0485 0.0657
0.0409 0.0690 0.0926
0.0660 0.0956 0.1196
shown in Figures 4 and 5. The CFD-predicted results of axial dispersion coef cient are given in Table 12. From Table 12, it can be observed that the axial dispersion coef cient increases with an increase HD=D ratio for both the columns. The proportionality relationship was found to be DL / HD0:65±1:72 . Since the effect of HD on DL has not been considered in the past, all the published correlations do not consider HD as a variable parameter. In order to understand the performance of the published correlations, the values of DL were calculated using different correlations. Figure 7 shows wide differences in the correlations. It is believed that the source of these differences is nonconsideration of the effect of HD on DL. Therefore, for all these published investigations, the CFD approach (of this paper) was employed and the predicted
where the value of constant a depends upon the extent of homogeneity for which the mixing time is measured. It may be noted that the mixing time is proportional to the square of the dispersion height (HD) assuming that DL is independent of HD. However, the experimental observations of mixing time (Table 9) conform to ymix HD0:40±0:66 . This experimental is quantitatively supported by the CFD (the predicted ymix is given in Table 9). (ii) As mentioned earlier, the axial dispersion coef cient is considered to be independent of dispersion height. In fact, this relationship has not been systematically investigated in the published literature. Therefore, it was thought desirable to probe into such a relationship using computational uid dynamics. The effect of height to diameter ratio (1–10) on axial dispersion coef cient was studied for column diameters 0.2 and 0.4 m and super cial gas velocity of 0.07– 0.295 m s¡1. The residence time distribution curves are
Figure 7. Dispersion coef cient from different correlations for column diameter 0.2 m. 1, Towell and Ackermann (1972); 2, Hikita and Kikukawa (1974); 3, Shah et al. (1978); 4, Joshi (1980); 5, Field and Davidson (1980); 6, Khare (1990); and 7, Krishna et al. (2000).
Trans IChemE, Vol 81, Part A, September 2003
CFD SIMULATION OF MIXING AND DISPERSION values were shown in Figure 6. The excellent agreement with all the published information shows the usefulness of CFD approach over the semi-empirical correlations.
FDz FL FVr
CONCLUSION A two-dimensional CFD code has been developed for two-phase gas–liquid dispersion in bubble columns. A detailed comparison has been presented between the CFD simulation and the experimental data reported by Hills (1974), Nottenkamper et al. (1983), Menzel et al. (1990), Yao et al. (1991), Yu and Kim (1991) and Grienberger and Hofmann (1992). Excellent agreement was observed between the predicted and the measured pro les of axial liquid velocity for all investigations. A complete energy balance was also established for all the cases. The ow model has been extended to the process of homogenization of an inert tracer. A systematic attempt has been made to investigate the effect of tracer addition and detector position and the column diameter, super cial gas velocity and height to diameter ratio on mixing time. A good agreement was observed between the predicted and the experimental values of mixing time for two column diameters (0.2 and 0.4 m i.d.) and three super cial gas velocities (0.07, 0.17 and 0.295 m s¡1). Further, the model has been extended for the estimation of axial dispersion coef cient (DL). Excellent agreement was found between the predicted and the experimental values of the axial dispersion coef cient measured over a wide range of D, VG and HD=D covered by various investigators in the past four decades. The parametric sensitivity of dispersion height has been investigated and in this regards, the usefulness of the CFD approach has been brought out as compared with the conventional semi-empirical correlations. NOMENCLATURE a b, c, d c c0 c· c0 u0 c0 v0 C Ci CB CD CL CV C0, C1 Ce1 Ce2 D Deff DL Dm Dt Dtr Dtz eR E ES fD FDr
constant in equation (22) empirical constant in equation (1) instantaneous concentration of the tracer uctuating concentration of the tracer time-averaged concentration of the tracer term de ned in equation (6) term de ned in equation (6) mean concentration of the tracer, de ned in equation (17) initial concentration of the tracer interface energy transfer factor drag force coef cient lift force coef cient virtual mass force coef cient drift ux constants model parameter in turbulent dissipation energy equation (ˆ1.44) model parameter in turbulent dissipation energy equation (ˆ1.92) diameter of the column, m effective diffusion coef cient, m2 s¡1 axial dispersion coef cient, m2 s¡1 molecular diffusivity, m2 s¡1 eddy diffusivity, m2 s¡1 radial eddy diffusivity, m2 s¡1 axial eddy diffusivity, m2 s¡1 rate of energy release for a single bubble, J s¡1 rate of energy released by all the bubble in the column, J s¡1 fraction of dissipation energy, J s¡1 drag force on a single particle in uidized bed, N m¡2 frictional force in the radial direction per unit volume of dispersion, N m¡3
Trans IChemE, Vol 81, Part A, September 2003
1001
frictional force in the axial direction per unit volume of dispersion, N m¡3 lift N m¡3 £ force per ¡unit volume ¢ of dispersion, ¤ ˆCL EL EG rL uG ¡ uL @uL =@r radial virtual mass force per unit volume of dispersion, N m¡3 µ » ¼ ¶ 1@ @ ˆCV EL EG rL r(vG ¡ vL ) ‡ (vG ¡ vL ) r @r @z
FVz
axial virtual mass force per unit volume of dispersion, N m¡3 µ » ¼ ¶ 1@ @ ˆCV EL EG rL r(uG ¡ uL ) ‡ (uG ¡ uL ) r @r @z
g
gravitational constant, m s¡2 µ Á» ¼ 2 n o » ¼ @vL v 2 @uL ˆ mt,L 2 ‡ L ‡ @r r @z
G H HD k n P PB Pe r R Re SF,K t u uG uL u· u0 v vB vG vL v· v0 VC VG VL VL(0) VSr VS z
2
! ³» ¼ » ¼ ´2 ¶ @vL @uL ‡ ‡ @z @r
height of the bubble column, m height of the gas dispersion, m turbulent kinetic energy per unit mass, m2 s¡2 empirical constant in equation (1) pressure, N m¡2 CB[FDrVSr ‡ FDzVS] Peclet number, (DL=VLHD) radial distance, m column radius, m Reynolds number source term of phase K time, s instantaneous axial velocity, m s¡1 axial component of gas velocity, m s¡1 axial component of liquid velocity, m s¡1 axial time-averaged velocity, m s¡1 axial uctuating velocity, m s¡1 instantaneous radial velocity, m s¡1 volume of single bubble, m3 radial component of gas velocity, m s¡1 radial component of liquid velocity, m s¡1 radial time-averaged velocity, m s¡1 radial uctuating velocity, m s¡1 circulation velocity, m s¡1 super cial gas velocity, m s¡1 super cial liquid velocity, m s¡1 center line liquid velocity, m s¡1 radial slip velocity between gas and liquid, m s¡1 axial slip velocity between gas and liquid, m s¡1 axial distance along the column, m
Greek symbols mK ‡ mt,K=sF GK a under relaxation parameter E fractional phase holdup EG fractional gas holdup EL fractional liquid holdup E· L time-averaged fractional liquid holdup E0L uctuating fractional liquid holdup e turbulent energy dissipation rate per unit mass, m2 s¡3 mK molecular viscosity of phase K, Pa s mt,K turbulent viscosity of phase K, Pa s n molecular kinematic viscosity of liquid, m2 s¡1 nt turbulent kinematic viscosity of liquid, m2 s¡1 ymix mixing time, s r density, kg m¡3 rG density of the gas, kg m¡3 rL density of the liquid, kg m¡3 s2 variance, de ned in equation (19), s2 sF turbulent Prandtl number for momentum transfer sf turbulent Prandtl number for bubble motion s2y dimensionless variance, de ned in equation (20) t mean residence time, de ned in equation (18), s Subscripts G K L
gas phase phase, K ˆ G, gas phase; L, liquid phase liquid phase
1002
EKAMBARA and JOSHI REFERENCES
Craver, M.B., 1984, Numerical computation of phase separation in two phase ow, J Fluids Eng, 106: 147–153. Deckwer, W.D., Burckhart, R. and Zoll, G., 1974, Mixing and mass transfer in tall bubble columns, Chem Eng Sci, 29: 2177–2188. Deshpande, N.S., Prasad, Ch.V., Kulkarni, A.A. and Joshi, J.B., 2000, Hydrodynamic characterization of dispersed two-phase ows by laser Doppler velocimeter, Trans IChemE, Part A, Chem Eng Res Des 78: 903–910. Field, R.W. and Davidson, J.F., 1980, Axial dispersion in bubble columns, Trans IChemE, 58: 228–236. Grienberger, J. and Hofmann, H., 1992, Investigation and modeling of bubble columns, Chem Eng Sci, 47: 2215–2220. Hillmer, G., Weismantel, L. and Hofmann, H., 1994, Investigations and modelling of slurry bubble columns, Chem Eng Sci, 49: 837–843. Hills, J.H., 1974, Radial nonuniformity of velocity and voidage in the bubble column, Trans IChemE, 2: 1–9. Hikita, H. and Kikukawa, H., 1974, Liquid phase mixing in bubble columns: effect of liquid properties, Chem Eng J, 8: 191–197. Hjertager, B.H. and Morud, K., 1993, Computational uid dynamics simulation of bioreactors, CFD Simul, 47–61. Johansen, S.T. and Boysan, F., 1988, Fluid dynamics in bubble stirred ladles: Part II. Mathematical modelling, Metall Trans B, 19B: 755– 764. Joshi, J.B., 1980, Axial mixing in multiphase contactors—a uni ed correlation, Trans IChemE, 58: 155–165. Joshi, J.B., 2001, Computational ow modelling and design of bubble column reactors, Chem Eng Sci, 56: 5893–5933. Joshi, J.B. and Sharma, M.M., 1979, A circulation cell model for bubble columns, Trans IChemE, 57: 244–251. Kataoka, I., Besnard, D.C. and Serizawa, A., 1992, Basic equation of turbulence and modelling of interfacial transfer terms in gas–liquid two-phase ow, Chem Eng Commun, 118: 221–236. Kato, Y., Nishiwaki, A., Fukuda, T. and Tanaka, S., 1972, The behaviour of suspended solid particles and liquid in bubble columns, J Chem Eng Japan, 5: 112–118. Khare, A.S., 1990, Development of radio active tracer techniques and its application to the design of multiphase reactors, Ph.D. thesis, University of Mumbai, India. Krishna, R., Urseanu, M.I., van Baten, J.M. and Ellenberger, J., 2000, Liquid phase dispersion in bubble columns operating in the churn-turbulent ow regime, Chem Eng J, 78: 43–51. Levenspiel, O., 1999, Chemical Reaction Engineering, 3rd edition (Wiley, New York). Menzel, T., Weide, T., Staudacher, O. and Onken, U., 1990, Reynolds stress model for bubble column reactor, Ind Eng Chem Res, 29: 988– 994. Moustiri, S., Hebrard, G., Thakre, S.S. and Roustan, M., 2001, A uni ed correlation for predicting liquid axial dispersion coef cient in bubble columns, Chem Eng Sci, 56: 1041–1047. Nottenkamper, R., Steiff, A. and Weinspach, P.M., 1983, Experimental investigation of hydrodynamics of bubble columns, Ger Chem Eng, 6: 147–155. Ohki, Y. and Inoue, H., 1970, Longitudinal mixing of the liquid phase in bubble columns, Chem Eng Sci, 25: 1–16. Patankar, S.V., 1981, Numerical Heat Transfer and Fluid Flow (McGrawHill, New York, USA). Patankar, S.V. and Spalding, D.B., 1972, A calculation procedure for heat, mass and momentum transfer in three dimensional parabolic ows, Int J Heat and Mass Transfer, 15: 126–129.
Peric, M., Kessler, R. and Scheuerer, G., 1988, Comparison of nite volume numerical methods with staggered and collocated grids, Comput Fluids, 16: 389–403. Ranade, V.V., 1992, Flow in bubble columns: some numerical experiments, Chem Eng Sci, 47(8): 1857–1869. Ranade, V.V., 1997, Modelling of turbulent ow in a bubble column reactor, Trans IChemE, Part A, Chem Eng Res Des, 75: 14–23. Reith, T., Renken, S. and Israel, B.A., 1968, Gas holdup and axial mixing in the uid phase of bubble columns, Chem Eng Sci, 23: 619. Riquarts, H.P., 1981, A physical model for axial mixing of the liquid phase for heterogeneous ow regime in bubble columns, Ger Chem Eng, 4: 18–23. Sato, Y. and Sekoguchi, K., 1975, Liquid velocity distribution in two-phase bubble ow, Int J Multiphase Flow, 2: 79–95. Shah, Y.T., Stiegel, G.J. and Sharma, M.M., 1978, Back-mixing in gas– liquid reactors, AIChEJ, 24: 369–400. Spalding, D.B., 1980, Mathematical Modelling of Fluid Mechanics, Heat Transfer and Chemical Reaction Processes, Imperial College Report HTS=80=1, Canbay, UK. Svendsen, H.F., Jakobsen, H.A. and Torvik, R., 1992, Local ow structures in internal loop and bubble column reactors, Chem Eng Sci, 47: 3297–3304. Thakre, S.S., Phanikumar, D.V., Khare, A.S. and Joshi, J.B., 1999, CFD modeling of ow, macro-mixing and axial dispersion in bubble column, Can J Chem Eng, 77: 826–837. Theofanous, T.G. and Sullivan, J., 1982, Turbulence in two-phase dispersed ows, J Fluid Mech, 116: 343–362. Towell, G.D. and Ackermann, G.H., 1972, Axial mixing of liquids and gas in large bubble reactor, in Proceeding of 2nd International Symposium on Chemical Reaction Engineering, Amsterdam, The Netherlands, B3.1–B3.13. van Baten, J.M. and Krishna, R., 2001, Eulerian simulations for determination of the axial dispersion of liquid and gas phases in bubble columns operating in the churn-turbulent regime, Chem Eng Sci, 56: 503–512. Van Doormal, J.P. and Reithby, G.D., 1984, Enhancement of the SIMPLE method for predicting incompressible uid ows, Numer Heat Transfer, 7: 147–163. Yao, B.P., Zheng, C., Gasche, H.E. and Hoffmann, H., 1991, Bubble behaviour and ow structure of bubble columns, Chem Eng Proc, 29: 65–75. Yu, Y.H. and Kim, S.D., 1991, Bubble properties and local liquid velocity in the radial direction of cocurrent gas–liquid ow, Chem Eng Sci, 46: 313–320.
ACKNOWLEDGEMENTS One of us (Ekambara) acknowledges the fellowship support given by the University Grants Commission (UGC), Government of India.
ADDRESS Correspondence concerning this paper should be addressed to Professor J. B. Joshi, Institute of Chemical Technology, University of Mumbai, Matunga, Mumbai-400 019, India. E-mail: [email protected] The manuscript was received 5 August 2002 and accepted for publication after revision 23 May 2003.
Trans IChemE, Vol 81, Part A, September 2003
www.ingentaselect.com=titles=02638762.htm
CFD SIMULATION OF MIXING AND DISPERSION IN BUBBLE COLUMNS K. EKAMBARA and J. B. JOSHI Institute of Chemical Technology, University of Mumbai, Mumbai, India
L
iquid phase mixing time was measured in 0.2 and 0.4 m i.d. columns over a wide range of super cial gas velocity (0.07–0.295 m s 1) and height-to-diameter ratio (1–10). A CFD model was developed for the prediction of ow pattern in terms of mean velocity and eddy diffusivity pro les. The predictions agree favourably with all the experimental data published in the literature. Complete energy balance was established in all cases. The validated CFD model was extended for the simulation of the macro-mixing process by incorporating the effects of both the bulk motion and the eddy diffusion. Excellent agreement was observed between the CFD predicted and experimental values of the mixing time over the entire range of D, VG, and HD=D covered. The model was further extended for the prediction of residence time distribution and hence the axial dispersion coef cient (DL). The predicted values of DL were found to agree very well with all the experimental data published in the literature. Keywords: bubble column; gas–liquid; mixing; residence time distribution; computational uid dynamics; axial dispersion coef cient.
the fact that it contains a single parameter. Many correlations have been presented for the axial dispersion coef cient such as those by Reith et al. (1968), Ohki and Inoue (1970), Towell and Ackermann (1972), Kato et al. (1972), Hikita and Kikukawa (1974), Deckwer et al. (1974), Shah et al. (1978), Joshi (1980), Field and Davidson (1980), Riquarts (1981), Krishna et al. (2000) and Moustiri et al. (2001). The results show that the axial dispersion coef cient depends very strongly on column diameter (D) and to a lesser extent on the super cial gas velocity (VG). Some publications also mention a slight in uence of kinematic viscosity n of the liquid phase. The following is a typical form of the correlation:
INTRODUCTION Bubble column reactors nd widespread application in industrial processes such as absorption, catalytic slurry reactions, bio-reactions, coal liquefaction and waste water treatment. Bubble columns offer many advantages as gas– liquid contactors due to their simple construction, low operating cost and the ease with which the liquid residence time can be varied. However, the scale-up of bubble columns is still poorly understood because of the complexity of ow patterns, and their unknown behaviour under different sets of design parameters such as diameter, heightto-diameter ratio, gas sparger type, etc. The knowledge of mixing and axial mixing in the liquid phase is essential for the modelling, design and optimization of bubble column reactors. Accordingly, there have been several experimental studies and the results have been presented in the form of empirical correlations. The important process parameter for the batch mixing process is the mixing time. Experimental determination of the mixing time involves introduction of a tracer at some location in the column (usually a pulse input) and measuring the tracer concentration as a function of time with the help of a sensor. Mixing time is typically considered as the time required for the tracer concentration to reach the desired degree of homogeneity. In a continuous operation of bubble columns, the liquid phase axial mixing and the axial dispersion model is commonly used. The liquid is assumed to be in plug ow with axial dispersion superimposed on it. This model lumps the different mechanisms of liquid mixing into a single axial dispersion coef cient. The popularity of the model is due to
DL ˆ nDb VGc nd
(1)
Table 1 gives the summary of experimental measurements and correlation of axial dispersion coef cient. The value of the power law index b varies between 1 and 1.5, that on the super cial gas velocity ranges between 0.25 and 0.5, while viscosity exerts such a slight in uence that it is either neglected altogether or quoted with a very small exponent. The above-mentioned correlations are semi-empirical and it is most desirable to develop predictive procedures on the basis of understanding of the hydrodynamics. An initial attempt in this direction has shown that the dispersion coef cient is proportional to the average liquid circulation velocity (VC; Joshi and Sharma, 1979; Field and Davidson, 1980), the value of VC was estimated using energy balance. The mixing process occurs due to the transport at three levels: molecular, eddy and bulk (convection). Usually, the bulk motion (or bulk diffusion) is superimposed on either molecular or eddy diffusion or both. The convective 987
DL ˆ 0:31DVL (0) DL ˆ VL HD =Pe(1 ¡ EG )
DL ˆ 0:068g0:375 D1:5 VG0:375 n¡0:125
0.52–5.5 4.25, 4.50 150, 200
Aqueous NaCl solution
Pulse Aqueous NaCl solution 4.0 174, 380, 630
@ 1@ @ (EL c) ‡ (rvEL c) ‡ (uEL c) r @r @t @z´ ³ ³ ´ 1@ @ @ @ ˆ rDm (EL c) ‡ Dm (EL c) r @r @r @z @z
(2)
c ˆ c· ‡ c0
v ˆ v· ‡ v0 u ˆ u· ‡ u0 EL ˆ E· L ‡ E0L
(3)
Substituting the above terms in equation (2) and incorporating the Reynolds averaging procedure (further details about solution have been given Joshi, 2001), we get: @ 1@ @ (·E c· ) ‡ (r·EL c· v· ‡ E· L c0v0 ) ‡ (·EL c· u· ‡ E· L c0 u0) r @r @t L @z ³ ´ ³ ´ 1@ @ @ @ ˆ rDm (·EL c· ) ‡ Dm (·EL c· ) (4) r @r @r @z @z @ 1@ @ (·E c· ) ‡ (r·EL c· v· ) ‡ (·EL c· u· ) r @r @t L @z ³ ´ ³ ´ 1@ @ @ @ ˆ rDm (·EL c· ) ‡ Dm (·EL c· ) r @r @r @z @z 1@ @ ¡ (r·E c0 v0) ¡ (·EL c0 u0) (5) r @r L @z The correlation of uctuating concentration and velocity can be modelled using Boussinesq’s approach: µ ¶ @·c c0 v0 ˆ ¡Dtr @r µ ¶ (6) @·c c0 u0 ˆ ¡Dtz @z substitution of equation (6) in equation (5) and subsequent rearrangement gives:
Moustiri et al. (2001)
Krishna et al. (2000)
5.0–35.0
Pulse
— — — — — Riquarts (1981)
Field and Davidson (1980)
4.0–5.5 18.9 3200
Radioactive
UV recorder
—
DL ˆ 0:31DVC
DL ˆ 0:678D1:4 VG0:3 DL ˆ 0:33DVC
DL ˆ 0:66D1:25 VG0:38 Transient state measurement
Steady state Electrolyte, dye and heat —
Aqueous KCl solution
7.23 — 200 —
4.3–33.8 1.5, 2.4 100, 190
Deckwer et al. (1974) Joshi (1980)
Hikita and Kikukawa (1974)
5.0–12.0 —
DL ˆ nD1:4 VG0:3 DL ˆ nDb VGc nd Steady state, transient pulse
Pulse Aqueous KCl solution 2.0, 4.05
Methylene blue dye 1.7–26.8 1.5, 1.9, 2.84 406, 1067
6.6, 12.2, 21.4 Kato et al. (1972)
Towell and Ackermann (1972)
4.0–18.0
Steady state Transient state measurement 7.5–45.0 2.0–25.0 3.53, 3.8 2.0, 3.0 Reith et al. (1968) Ohki and Inoue (1970)
14, 29 40, 80, 160
Aqueous NaCl solution Potassium chloride
Technique of measurement Tracer Column height (m)
Range of VG (cm s¡1)
transport can be characterized by the mean velocity components in the liquid phase. The dispersive transport can be characterized by eddy diffusivity in the liquid phase. In view of the controlling role of the convection and eddy diffusion, the solution of convection-diffusion equation for the liquid phase can give insight into the mechanism of axial mixing and hence a rational expression for the effective axial dispersion coef cient. The transient mass balance for a tracer substance in a two-dimensional axisymmetric cylindrical co-ordinate system is given by the following equation:
The instantaneous concentration (c) may be written as the sum of the time-averaged concentration c· and a concentration uctuation c0. Similarly, the velocity components and hold-up can also be written:
Column diameter (mm) Author
Table 1. Axial dispersion coef cient: experimental details and correlation of previous work.
DL ˆ 0:30D2 VG1:2 ‡ 170d DL ˆ 0:30D2 VG1:2 ‡ 170d
EKAMBARA and JOSHI
Correlation
988
@ 1@ @ (·E c· ) ‡ (·EL c· v· ) ‡ (·EL c· u· ) r @r @t L @z µ ¶ 1@ @ ˆ r(Dm ‡ Dtr ) (·EL c· ) r @r @r µ ¶ @ @ ‡ (Dm ‡ Dtz ) (·EL c· ) @z @z
(7)
LDA measurements in bubble column have shown that (Deshpande et al., 2000) the rms velocities in the radial and axial direction are practically equal with deviation at some location within 0.56%. In view of this the sensitivity of Trans IChemE, Vol 81, Part A, September 2003
CFD SIMULATION OF MIXING AND DISPERSION eddy diffusivity in axial (Dtz) and radial (Dtr) dimensions was analysed and it will be shown later that the variation in Dtr and Dtz (within 0.63%) has a nominal effect. In view of this it will be assumed that Dtr ˆ Dtz ˆ Dt. With these simpli cations, equation (7) takes the following form: @ 1@ @ (·EL c· ) ‡ (r·vEL c· ) ‡ (u·EL c· ) r @r @t @z¶ µ µ ¶ 1@ @ @ @ ˆ rDeff (·EL c· ) ‡ Deff (·EL c· ) r @r @r @z @z
MODEL FORMULATION Equation of Continuity and Motion and the k±e Model For the case of bubbly ow, the equations of continuity and motion for a two-dimensional cylindrical coordinate system can be represented in the following generalized form: 1@ @ (rErvF)K ‡ (EruF)K r @r @z ³ ´ ³ ´ 1@ @F @ @F ˆ rEG ‡ EG ‡S r @r @r K @z @z K F,K
(8)
where Deff denotes the effective (molecular and eddy) diffusion coef cient. The above equation assumes that (i) the solute does not undergo any chemical reaction and (ii) all of the physical properties of the system remain constant. Thakre et al. (1999) have used computational uid dynamics (CFD) to obtain the ow pattern (radial and axial pro les of eddy diffusivity and axial, radial mean velocities) and used the ow pattern for the prediction of mixing time by solving the above equation. They compared the predictions with the experimental data from a 0.2 m i.d. column with height to diameter ratio of 8 and found excellent agreement. In the present work, the CFD model has been extended for the prediction of axial dispersion coef cient over a wide range of D, VG and HD=D. A comprehensive comparison has also been presented with the experimental results published in the literature.
989
(9)
where F is the transport variable (for instance, F is one for equation of continuity or it is a velocity component for the respective equation of motion, and F ˆ k and e for the conservation equation for the turbulent kinetic energy and turbulent energy dissipation rate respectively). K denotes the phase (K ˆ G or L) and SF is the source term for the respective dependent variable. The values of F and SF for different transport variables have been given in Table 2. From Table 2, it can be seen that most of the terms are derived from gas and liquid velocities. In addition, in twophase ows, momentum and energy transfer occurs across the interface. Therefore, the drag force (FDZ and FDR), virtual mass force (Fvz and Fvr), and lift force (FLr) appear in the axial and radial components of the momentum balance (Table 2). Furthermore, the interphase transfer of
Table 2. Two-dimensional k–e model for gas–liquid ows governing equations. The governing equations written in a general form: ³ ´ ³ ´ 1 @ @ 1 @ @F @ @F EG (rErvF)K ‡ (EruF)K ˆ rEG ‡ ‡S r @r @z r @r @r K @z @z K F,K Conservation of:
F
sF
sf
Mass
1
1
1 to 1
Axial velocity momentum
u
1.0
1 to 1
Radial velocity momentum
v
1.0
1 to 1
Turbulent kinetic energy
k
1.0
—
Turbulent dissipation energy
e
1.3
SF,K ˆ source terms
³ ´ ³ ´ 1 @ mt,K @EK @ mt,K @EK r ‡ r @r sf @r @z sf @z µ ³ ´ ³ ´¶ @P 1 @ @v @ @u ¤ ¤ EL § FVS EL ‡ Emt ¡EK ‡ EK rK g § FDZ rEmt ‡ @z r @r @z @z @z K µ ³ ´ ³ ´¶ ³ ´ ³ ´ mt @E 1 @ mt @E @ mt @E 1 @ @u ‡ uK r ‡ ‡ (rv) ‡ sf @z K r @r r @r sf @r @z sf @z K @z K µ ³ ´ ³ ´¶ @P 1 @ @v @ @u ¤ ¤ ¤ EL § FVS EL § FLr EL ‡ Emt ¡EK § FDr rEmt ‡ @r r @r @z @z @r K µ ³ ´ ³ ´¶ ³ ´ µ ¶ ± ² mt @E v 1 @ mt @E @ mt @E 1 @ @u ¡ 2Emt 2 ‡vK r ‡ ‡ (rv) ‡ r K r @r sf @r @z sf @z K sf @r K r @r @z K
—
EL (G ‡ PB ¡ rL e) e EL [Ce1 (G ‡ PB ) ¡ Ce2 rL e] k @u ˆ ¡CL EG rL (uG ¡ uL ) L @r
mt,K ˆ 0:09rk (k 2 =e); PB ˆ CB (FDr VSr ‡ FDz VS ); FLr » ¼ @ @ Fvz ˆ ¡CV EG rL (u ¡ uL ) ‡ (uG ¡ uL ) ; FDz ˆ ¡4:9 £ 104 EG (uG ¡ uL ) @r G @z µ ¶ 1 @ @ Fvr ˆ ¡CV EG rL r(vG ¡ vL ) ‡ (vG ¡ vL ) ; FDr ˆ ¡4:9 £ 104 EG (vG ¡ vL ) r @r @z " Á» ¼ ¼ 2 ! ³» ¼ » ¼ ´2 # 2 n o2 » @vL v @uL @vL @uL G ˆ mt,L 2 ‡ L ‡ ‡ ‡ @r r @z @z @r Force terms are positive for the gas phase and negative for the liquid phase.
Trans IChemE, Vol 81, Part A, September 2003
990
EKAMBARA and JOSHI
energy (PB) appears in the equations for turbulent kinetic energy and turbulent energy dissipation rates.
eD ˆ vB (rL ¡ r G)gVS
Boundary and Initial Conditions Boundary and initial conditions for equation (9) (i) Along the axis: axisymmetry in uG, uL, EG, vG, vL, k and e. (ii) Along the wall: the velocities satisfy the no-slip boundary conditions (the wall function method based on the log law of the wall is invoked to calculate the wall shear stress and the values of k and e close to the wall). (iii) At the inlet: gradients of vG, vL, k and e are set to zero. (iv) At the top surface of the computational domain, the gradient of the dependent variable is set to zero. For initiating the numerical solution, EG and uG are speci ed at the inlet. At all the other locations uG, uL, vG, vL and EG are taken to be zero at t ˆ 0. For k and e, the initial guess values were found to be important. Boundary and initial conditions for equation (8) The solution to the diffusion-convection equation (8) must satisfy the following conditions: c· (r, z, 0) ˆ 0 c· (r, 0, 0) ˆ Ci @·c(0, z, t) @·c(R, z, t) ˆ ˆ0 @r @r
for all z and t
Energy Balance Consider a single bubble (of volume vB) rising in a pool of liquid. The bubble at any location has pressure energy and potential energy. As the bubble rises, the pressure energy decreases and the potential energy increases. The rate of change of energy is given by the following equation: dh ˆ v(rL ¡ rG )gVS ˆ fD VS (10) dt where VS is the slip velocity and fD is the drag force. When the bubble rise is in a creeping ow regime, the energy released during rise is dissipated in a viscous manner. In other words, the gas phase energy gets directly converted into internal energy. When the bubble rise is in the turbulent ow regime (Re1 > 500), the energy released during the rise is rst of all converted to turbulence in the liquid phase and nally into internal energy at the Kolmogorov scale. Therefore, the fraction of gas phase energy used for generating liquid phase turbulence depends upon the bubble Reynolds number. There is one more aspect of the energy balance. In bubble columns, the bubble rise velocity with respect to the wall can be very different from the slip velocity because of liquid circulation. When the actual velocity is higher than the slip velocity (as in the central upward region), the rate of energy release is different from that given by equation (10). vB (rL ¡ rG )g
eR ˆ vB (r L ¡ rG)g(VS ‡ uL )
where uL is the liquid velocity at the location of bubble rise. The rate of energy converted into turbulence and spent in the viscous dissipation depends only on the slip velocity:
(11)
(12)
Therefore, when uL is positive, the rate of energy release is higher than the dissipation. The extra energy is used for maintaining the liquid circulation. In contrast, when the bubble moves downward in the near-wall region, the liquid phase energy is supplied to bubbles. The overall effect of up ow and down ow (because of the gas hold-up pro le, more bubbles are present in the up ow region as compared to the down ow region) is a net energy supply from gas phase to liquid phase. The rate of energy release from the gas phase to the liquid phase is given by the following equation: p E ˆ D2 (rL ¡ rG )gHD E·L (VG ¡ E·G VS ) (13) 4 The above equation assumes that the slip velocity is same for all the bubbles. The rate of energy converted into turbulence and spent in viscous energy dissipation is given by: p ED ˆ D2 (rL ¡ rG )gHD E·G E· L VS (14) 4 It may be, however, noted that the turbulence generated by bubbles has a scale corresponding to bubbles, whereas the turbulence generated in liquid circulation [using the energy given by equation (13)] has a scale corresponding to column dimensions. Therefore, this gives rise to a very important question regarding the quantitative role of the bubble generated turbulence to the transport phenomena. Sato and Sekoguchi (1975) and Theofanous and Sullivan (1982) have addressed to this problem. It was assumed that a fraction of CB is used in the transport and this fraction forms a source terms for the liquid phase k (turbulent kinetic energy per unit mass) equation: p ES ˆ CB D2 (rL ¡ r G)gHD E· L E· G VS (15) 4 where CB takes values between 0 and 1. Johansen and Boysan (1988) have included the radial variation of slip velocity and CB was in the range of 0.1–0.2. Later Svendsen et al. (1992) used CB as a tuning parameter to match their experimental results. Kataoka et al. (1992) reported a detailed analysis of the extra terms due to the presence of the dispersed phase in the source of turbulent kinetic energy. Their analysis suggests that the extra generation of turbulence due to large bubbles is almost compensated by the extra dissipation due to the small-scale interfacial structures. The studies by Hillmer et al. (1994) and Ranade (1992, 1997) also suggest zero value for CB. Hjertager and Morud (1993) have used a very low value of 0.02. The total energy released from the gas phase to the liquid phase is given by equation (13). Furthermore, the quantitative role of bubble-generated turbulence is given by equation (15). Therefore, the sum of equations (13) and (15) gives the total energy received by the liquid phase and is expressed as: p E ‡ ES ˆ D2 (rL ¡ rG )gHD E·L [VG ‡ (CB ¡ 1)·EG VS ] 4 (16) Trans IChemE, Vol 81, Part A, September 2003
CFD SIMULATION OF MIXING AND DISPERSION It must be emphasized at this stage that, whatever may be the value of CB, the energy balance must be satis ed. When the k–e model is used for the prediction of ow pattern, we get radial and axial variation of e (energy dissipation rate per unit mass) as one of the answers. From this e eld, the total energy dissipation rate can be calculated by suitable volume integration. The total energy dissipation rate must equal the energy-input rate given by equation (16). METHOD OF SOLUTION The mathematical model was solved in two steps. In the rst step, equations of continuity and motion were solved (together with k–e equations) for getting complete ow pattern in terms of distribution of gas and liquid velocities, eddy diffusivity and gas holdup. This ow pattern was then used in step two, for the estimation of mixing time and residence time distribution. The set of equations given in Table 2 was solved numerically, which consisted of the following steps: (i) the generation of suitable grid points; (ii) the conversion of governing equations into algebraic equations; (iii) the selection of discretization schemes; (iv) the formulation of discretzed equations at every grid location; (v) the formulation of pressure equation; and (vi) the development of suitable iteration scheme for use in obtaining the nal solution. The nite control volume technique of Patankar (1981) was employed for the solution of equations given in Table 2. A staggered grid arrangement proposed by Patankar and Spalding (1972) was employed. It consisted of 50 £ 100 grid points with 50 grid points in the radial direction and 100 grid points in axial direction. The velocity components were calculated for the points that lie on the faces of the control volume, while all other variables were calculated at the centre of the control volume. Several discretization schemes (upwind, hybrid, exponential and the power law) were tested. Out of these, the power law scheme was found to be the best and was used for the discretization of the governing equations. The most important numerical concern is the manner in which the pressure eld is computed and its effect on the continuity of each eld. For multidimensional analyses of single phase and homogenous two phase ows, many algorithms like SIMPLE, SIMPLER (Patankar, 1981), SIMPLEST (Spalding, 1980) and SIMPLEC (Van Doormal and Raithby, 1984) have been proposed for solving the pressure velocity coupling. All the above-mentioned algorithms incorporate the pressure equation, thus using the staggered grid concept, but differ principally in the exact form of the pressure equation, the manner in which the equations are formulated over the control volume, the degree of implicitness used in time advancement, and the solution sequence (Craver, 1984). Although other numerical schemes have been found to give faster convergence than SIMPLE, one iteration of these schemes involves more computational effort than SIMPLE (Patankar, 1981). Moreover, as compared with other schemes, the SIMPLE algorithm has been extensively used and has served well in a number of CFD simulation studies. In an iterative type of solution, it is a general practice to adopt some kind of convergence criteria. In fact, it is this which will decide not only the accuracy but also the Trans IChemE, Vol 81, Part A, September 2003
991
Table 3. Effect of internal iterations on the rate of convergence. Pressure
vr
vz
ur
uz
k
e
Iterations
3 4 5 6 10
2 2 2 4 6
2 2 2 4 6
2 3 1 4 6
2 3 1 4 6
10 4 5 8 20
10 6 5 10 20
624 624 624 624 624
computational time required and stability of the given problem. However, if convergence criteria are given without prior knowledge of the numerical problem, it often ends up in using excessive computational time. Hence, in the present simulation the role of the convergence criteria on the ow variable was critically analysed. Effect of Internal Iterations For an iterative process, it is a common practice that two or three iterations are used for solving the set of algebraic equations for vr, vz, ur, uz, k and e. The number of iterations used for pressure correction equation is typically four to ve times than that for the above-mentioned variables when the line-by-line method TDMA method is used. However, the above-mentioned practice does not give better convergence. The coef cients of the velocity variables are in uenced by the eddy viscosity. Therefore, the solution of the k and e equations should be iterated more times than the pressure correction equation. However, as seen from Table 3, the variation of number of internal iterations for k and e as well as other variable gives the same rate of convergence of the solution (i.e. 624 iterations), indicating its independence of the interal iterations. Effect of the Under-relaxation Parameter It has been shown by Peric et al. (1988) that a nearly optimum convergence is found for the SIMPLE algorithm when the under-relaxation parameters for pressure and velocity variables are selected by the following relation ap ˆ 1 7 au. In addition, in the present study it has been observed that the under-relaxation parameters for k and e are also imporatant to obtain a converged solution. It was found that, by decreasing the under-relaxation parameter for pressure correction equation, ap from 0.75 to 0.1, the total number of iterations required for convergence increased by 18%. A similar observation was found with the axial liquid velocity relaxation equation. However, no improvement in the rate of convergence was found by increasing the under-relaxation parameters of the radial and axial gas phase momentum equations and turbulent energy dissipation rate equation. In Table 4, the effect of under-relaxation Table 4. Effect of under relaxation parameters on the rate of convergence. ap 0.75 0.10 0.75 0.75 0.75
avr
avz
aur
auz
ak
ae
Iterations
0.75 0.75 0.75 0.75 0.75
0.75 0.75 0.75 0.75 0.75
0.75 0.75 0.75 0.5 0.75
0.75 0.75 0.5 0.75 0.75
0.75 0.75 0.75 0.75 0.75
0.75 0.75 0.75 0.75 0.95
624 818 630 620 567
0.002 0.002 0.002 0.002 0.002 0.002 0.002 0.002 0.012 0.2 0.2
684 684 684 684 684 684 684 671 634 616 622
0.6=3.45
Five-point conductivity probe Height measurement 3.50 20–80 Air–water
0.475 10–140
Hole diameter, 0.2 mm Perforated pipes: 6 mm, 8 nos, hole diameter 1 mm, 78 nos. Hole diameter 0.2 mm
Air–deionized water Air–water
20–80
3.49
Electroconductivity microprobe Electroconductivity microprobe Two channel optical bre probe 4.86 12–96 Air–water
Optical probe — 53–1452 Air–water
Liquid velocity
Pavlov tube Electroconductivity needle probe 4.34 19–169 Air–water
Holdup Sparger
Seive plate: Hole diameter 0.4 mm, 61 nos Seive plate hole diameter, 1 mm —
Measurement location HD=D
Fly-wheel anemometer Hot lm anemometry Hot lm anemometry Two Pt–Rh electrodes, two U shaped optical bres and a tracer injection
0.29=4.50
0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.005 0.015 0.015 0.15
Grienberger and Hofmann (1992)
0.005 0.005 0.005 0.005 0.005 0.005 0.05 0.005 0.005 0.005 0.005
0.254=2.50
¡0.02 ¡0.04 ¡0.5 0.5 ¡0.02 ¡0.02 ¡0.02 ¡0.02 ¡0.02 ¡0.02 ¡0.02
Yu and Kim (1991)
0.33 0.33 0.33 0.53 0.33 0.33 0.33 0.33 0.33 0.33 0.33
0.29=4.50
Iterations
Yao et al. (1991)
e
Menzel et al. (1990)
k
0.45=5.60
uz
Nottenkamper et al. (1983)
vz
0.138=1.37
0.02 0.02 0.02 0.02 0.02 0.05 0.02 0.02 0.02 0.02 0.02
ur
Hills (1974)
vr
Researchers
Table 5. Effect of initial guesses of the ow variables on the rate of convergence for a given set of under-relaxation parameters.
Super cial gas velocity (mm s¡1)
The effect of the initial guess values on the rate of convergence is tabulated in Table 5. It is interesting to note that the initial guess values of vr, vz, ur, uz and p have no in uence on that rate of convergence. The rate of convergence of the iterative process primarily depends on the guess values of k and e. However, it is dif cult to establish any mathematical relation between the guess values of k and e that yield an optimum rate of convergence. In the present work, we have investigated the performance of SIMPLE, SIMPLER and SIMPLEC, and SIMPLE was found to be the best. The set of algebraic equations obtained after discretization was solved by TDMA. The relaxation parameters, initial guess values and internal iterations for the variables were tuned to optimize the balance between the convergence criteria (1.0 £ 10¡3) and the number of iterations required. While the solution was found to be independent of the initial guesses and internal iterations (for vr, vz, ur, uz and p), the relaxation parameters had a dominant effect on the convergence of solution. The effect was prominent for lower height to diameter ratios (HD=D) and higher super cial gas velocities (VG), enforcing the use of lower relaxation parameters. Further details about method of solution have been given Joshi (2001). After getting the ow pattern, the model was extended to the understanding of the mixing process. It has been assumed that the tracer concentration has no in uence on the ow characteristics. Therefore, ow results computed under fully developed and steady-state conditions were used to solve the transient scalar mixing problem. It was assumed that the eddy diffusivities for momentum and mass are equal. The transport equation for concentration in a cylindrical coordinate system is given in equation (8). A separate code was developed to simulate transient mixing in
Table 6. Comparison of radial gas hold-up and mean axial liquid velocity pro les: experimental data of previous work.
Effect of Guess Values of the Variables on the Rate of Convergence
Gas–liquid system
parameters of different variable has been given in terms of number of iterations. It can be seen from the above table that an optimum convergence can be achieved for the values of ak ˆ 0.5 and ae ˆ 0.95 with other parameters taking a constant value of 0.75. It may be emphasized that the above optimum combination is for a particular set of VG, D and HD=D. For another set the optimum of a is also different and the selection gets facilitated by the experience in CFD.
Measurement techniques employed
EKAMBARA and JOSHI
Column diameter= height (m)
992
Trans IChemE, Vol 81, Part A, September 2003
CFD SIMULATION OF MIXING AND DISPERSION bubble columns. This code uses the ow patterns results generated in the rst step. In the present case of scalar mixing, there is no source term.
993
to extend the model for the prediction of mixing time and the axial dispersion coef cient. Comparison Between the Predicted and Experimental Data for Liquid Phase Mixing Time
RESULT AND DISCUSSION Comparison Between the Predicted and Experimental Data for Flow Pattern As a rst step, it is important to establish the validity of the model for ow pattern. Therefore, comparison has been made with experimental data of Hills (1974), Nottenkamper et al. (1983), Menzel et al. (1990), Yao et al. (1991), Yu and Kim (1991) and Grienberger and Hofmann (1992). Table 6 gives the details pertaining to column diameter, column height, sparger design, gas–liquid system, range of VG and VL, the measurement location and the measurement techniques employed for the measurement of mean axial liquid velocity and fractional holdup. It may be noted that practically all the major experimental investigations published in the literature have been covered in the table. The results are given in Table 7 and Figure 1A–F. The agreement between the predicted and the experimental pro les can be seen to be excellent. The agreement over the wide range of D, VG and the gas–liquid systems ensures the applicability of the model (for the estimation of ow pattern). In view of the abovementioned favourable comparison, it was thought desirable
The knowledge of ow pattern was extended for the prediction of liquid phase mixing time (ymix). A systematic study was undertaken to investigate the effect of various parameters such as the locations of tracer addition and detection, the quantity of tracer, the effect of time step, column diameter, super cial gas velocity (VG) and heightto-diameter ratio on the mixing time. In order to understand the role of detector position, the detectors were placed at the grid locations of 3 £ 100, 3 £ 50, 3 £ 3, 25 £ 60, 50 £ 50 and 50 £ 3. For all these cases, the tracer was introduced at the 50 £ 90 location. The results of concentration vs time pro les are shown in Figure 2. It can be seen that all the pro les originate from zero and reach a constant value of 0.162, the nal concentration. The concentration pro le at 3 £ 100 detector (close to the tracer input location) position shows a peak overshoot. For all the other detectors, the response was overdamped. However, the time required for homogenization can be seen to be independent of detector position. In second step, the effect of the location of tracer introduction was investigated. For this purpose, tracer was introduced at 50 £ 90, 3 £ 100, 3 £ 50, 3 £ 3, 25 £ 50 and 50 £ 3
Table 7. Comparison between CFD predictions and the experimental observations. Material balance Authors
¡1
Energy balance
VG (m s )
E·G
CB
C0
C1
LHS
RHS
0.0206 (0.019) 0.038 (0.038) 0.067 (0.064) 0.106 (0.095) 0.171 (0.169)
0.065 (0.067) 0.106 (0.107) 0.161 (0.160) 0.197 (0.194) 0.218 (0.212)
0.43 0.44 0.37 0.28 0.34
3.84 (2.0) 3.34 (2.0) 2.58 (2.0) 2.10 (2.0) 1.91 (2.0)
0.222 (0.25) 0.24 (0.25) 0.25 (0.25) 0.32 (0.25) 0.465 (0.25)
0.105 0.127 0.210 0.355 0.857
0.115 0.205 0.221 0.366 0.831
Nottenkamper et al. (1983)
0.054 (0.053) 0.106 (0.105) 0.331 (0.324) 0.601 (0.823) 1.209 (1.452)
0.121 (0.120) 0.175 (0.168) 0.312 (0.271) 0.405 (0.409) 0.494 (0.497)
0.36 0.38 0.41 0.21 0.23
2.64 (2.0) 2.5 (2.0) 1.64 (2.0) 1.63 (2.0) 1.61 (2.0)
0.375 (0.40) 0.35 (0.40) 0.54 (0.40) 0.39 (0.40) 0.40 (0.40)
0.190 0.502 1.815 2.38 6.59
0.198 0.513 1.836 2.29 5.80
Menzel et al. (1990)
0.012 (0.012) 0.026 (0.024) 0.052 (0.048) 0.096 (0.096) 0.021 (0.02) 0.041 (0.04) 0.062 (0.06) 0.081 (0.08)
0.029 (0.029) 0.049 (0.048) 0.105 (0.103) 0.137 (0.128) 0.057 (0.057) 0.108 (0.106) 0.156 (0.153) 0.186 (0.180)
0.36 0.33 0.32 0.44 0.37 0.44 0.43 0.46
3.13 (2.5) 5.80 (2.5) 5.30 (2.5) 3.02 (2.5) 2.96 (2.0) 2.71 (2.0) 2.28 (2.0) 2.69 (2.0)
0.39 (0.40) 0.36 (0.40) 0.42 (0.40) 0.35 (0.40) 0.285 (0.25) 0.236 (0.25) 0.251 (0.25) 0.223 (0.25)
0.057 0.130 0.223 0.593 0.075 0.208 0.299 0.417
0.07 0.157 0.263 0.588 0.078 0.233 0.312 0.447
Yu and Kim (1991)
0.011 (0.01) 0.02 (0.02) 0.040 (0.04) 0.061 (0.06) 0.101 (0.1) 0.142 (0.14)
0.025 (0.025) 0.035 (0.034) 0.049 (0.048) 0.062 (0.060) 0.090 (0.084) 0.136 (0.130)
0.28 0.30 0.21 0.20 0.27 0.22
3.17 (2.1) 2.4 (2.1) 2.5 (2.1) 2.23 (2.1) 1.87 (2.1) 1.77 (2.1)
0.39 (0.40) 0.554 (0.73) 0.739 (0.73) 0.853 (0.73) 0.76 (0.73) 0.80 (0.73)
0.045 0.122 0.187 0.322 0.526 0.661
0.054 0.100 0.194 0.288 0.549 0.668
Grienberger and Hofmann (1992)
0.021 (0.02) 0.081 (0.08)
0.057 (0.054) 0.156 (0.150)
0.37 0.46
2.96 (2.0) 2.70 (2.0)
0.285 (0.25) 0.22 (0.25)
0.075 0.417
0.078 0.447
Khare (1990)
0.072 (0.07) 0.173 (0.17) 0.292 (0.295)
0.136 (0.135) 0.223 (0.223) 0.304 (0.282)
0.30 0.32 0.39
2.17 (3.01) 1.99 (3.01) 1.76 (3.01)
0.38 (0.30) 0.43 (0.30) 0.51 (0.30)
0.362 0.821 1.292
0.329 0.815 1.488
Hills (1974)
Yao et al. (1991)
The bracketed numbers indicate the experimental values; LHS is volume integral of energy dissipation rate and RHS is equation (16).
Trans IChemE, Vol 81, Part A, September 2003
994
EKAMBARA and JOSHI
locations. For all these cases, the detector position was at 50 £ 100. The results are given in Table 8. It can be observed that the mixing time is practically independent of the point of tracer introduction. The same result was found to hold for many other detector locations. The effect of tracer quantity on mixing time was also investigated and was found to have no effect on mixing time. These numerical results (for the effects of tracer quantity and locations of tracer introduction and detection) were also con rmed experimentally. The effect of height to diameter ratio (HD=D) on mixing time was studied in 0.2 and 0.4 m i.d. column for a gas–liquid system with variation in VG from 0.07 to 0.295 m s¡1, respectively. For 0.2 m i.d. column HD=D ratio was varied from 3 to 10 and from 1 to 5 for 0.4 m i.d. column. The results have been shown in Figure 3A–C for column diameter 0.2 m and super cial gas velocity from 0.07 to 0.295 m s¡1. The comparison between the predicted and experimental data (Khare, 1990) of mixing time for different column diameters, super cial gas velocities and height to diameter ratios are given in Table 9. It can be observed from Table 9 that, for a particular VG, the mixing time increases with an increase in HD=D ratio for both the columns. From Table 9, it can also be noted that, at the same HD=D ratio the mixing time increases with an increase in column diameter (D). Further, it can be seen that the mixing time decreases with an increase in VG. Comparison Between the Predicted and Experimental Data for Axial Dispersion Coef cient Residence time distribution curves were obtained from the CFD simulation. For this purpose, concentrations were estimated at the topmost axial location and at all the radial locations. The average cross-sectional concentration at the exit gives the mixing cup concentration and the values were estimated with respect to time. The resulting residence time distribution curves are shown in Figures 4 and 5 for 0.2 and 0.4 m i.d. column with variation in VG from 0.07– 0.295 m s¡1, respectively. For 0.2 m i.d. column HD=D ratio was varied from 3 to 10 and from 1 to 5 for 0.4 m i.d. column. The range of super cial liquid velocity has been given in Table 12. These curves were used for the estimation of axial dispersion coef cient following the stepwise procedure suggested by Levenspiel (1999). The mean solute concentration is given by: … 1 R C ˆ 2 2p·cr dr (17) pR 0 When the reactor is modelled as a closed vessel, the rst and second moments of the residence time distribution curve are related to the mean residence time (t), and the axial dispersion coef cient, DL, respectively (Levenspiel, 1999) as: P tC(t)Dt t P (18) C(t)Dt s2
P2 t C(t)Dt P ¡ t2 C(t)Dt
(19)
Figure 1. Comparison between the simulated and experimental pro les of fractional gas holdup and axial liquid velocity (refer to Table 3 for experimental details). (A) Hills (1974); (B) Nottenkamper et al. (1983); (C) Menzel et al. (1990); (D) Yao et al. (1991); (E) Yu and Kim (1991); (F) Grienberger and Hofmann (1992).
Trans IChemE, Vol 81, Part A, September 2003
CFD SIMULATION OF MIXING AND DISPERSION
Figure 1. Continued.
Trans IChemE, Vol 81, Part A, September 2003
995
996
EKAMBARA and JOSHI
s2y
³ ´µ ³ ´¶ s2 DL DL V L HD ˆ 2 ˆ2 ¡2 1 ¡ exp ¡ t VL HD VL HD DL (20) VL HD 2 sy 2
DL
Figure 2. Concentration pro les at various detector positions for VG 0.17 m s 1, D 0.2 m and HD=D 10.
(21)
where t is the mean residence time, DL is the axial liquid dispersion coef cient, HD is the height of the column and s represents the variance. The effect of column diameter (D) and super cial gas velocity (VG) on axial dispersion coef cient was studied. The comparison between the CFD predicted results and experimental data of different investigators are given in Table 10. It can be observed that the axial dispersion coef cient increases with an increase in VG with a proportionality in the range of VG0:25±0:42 . Further, it can also be
Table 8. Effect of tracer introduction positions on liquid phase mixing time. Liquid phase mixing time for various tracer introduction positions in the column (s) VG (m s¡1) 0.07 0.17 0.295
E· G
50 £ 90
3 £ 100
3 £ 50
3£3
25 £ 50
50 £ 3
0.135 0.223 0.282
71.89 67.81 56.63
71.89 67.81 56.63
71.89 67.81 56.63
71.89 67.81 56.63
71.89 67.81 56.63
71.89 67.81 56.63
0.07; (B) VG
0.17; and (C) VG 0.295m s 1; 1—3
Figure 3. Effect of HD=D ratio on mixing time for column diameter 0.20 m. (A) VG 3—25 60; 4—50 50; 5—3 3; 6—50 3.
100; 2—3
50;
Trans IChemE, Vol 81, Part A, September 2003
CFD SIMULATION OF MIXING AND DISPERSION
Figure 3. Continued.
Trans IChemE, Vol 81, Part A, September 2003
997
998
EKAMBARA and JOSHI Table 9. Mixing time (s). HD=D
Column diameter ˆ 0.2 m VG (m s¡1) 0.070 0.170 0.295
3
8
10
E· G
Expt ymix
CFD ymix
Expt ymix
CFD ymix
Expt ymix
CFD ymix
Expt ymix
CFD ymix
0.135 0.223 0.282
33.5 29.3 25.2
29.8 30.0 24.2
38.4 35.2 31.2
39.0 36.2 32.0
51.6 45.6 39.6
52.0 45.1 40.3
72.5 67.2 —
71.9 67.8 56.6
Column diameter ˆ 0.4 m 0.070 0.170 0.295
5
0.135 0.223 0.282
2 56.5 51.6 47.0
3 57.1 52.0 45.8
66.0 61.2 55.5
Figure 4. Effect of HD=D ratio on RTD for column diameter 0.20 m. (A) VG 0.07; (B) VG 0.17; and (C) VG 0.295 m s 1.
4 65.2 59.7 53.8
78.0 71.5 64.0
5 75.9 69.4 65.2
88.0 79.5 —
86.3 75.8 69.4
Figure 5. Effect of HD=D ratio on RTD for column diameter 0.40 m. (A) VG 0.07; (B) VG 0.17; and (C) VG 0.295 m s 1.
Trans IChemE, Vol 81, Part A, September 2003
CFD SIMULATION OF MIXING AND DISPERSION
999
Table 10. Comparison between CFD predictions and the experimental observations. Material balance Authors
VG (m s¡1)
Energy balance
E· G
C0
C1
LHS
RHS
Experimental DL (m2 s¡1)
Predicted DL (m2 s¡1)
Ohki and Inoue (1970)
0.02 (0.02) 0.048 (0.05) 0.098 (0.10) 0.148 (0.15) 0.21 (0.20)
0.063 (0.062) 0.123 (0.122) 0.205 (0.218) 0.301 (0.307) 0.380 (0.382)
4.223 3.699 2.568 2.042 1.923
0.226 0.238 0.262 0.257 0.213
0.103 0.249 0.485 0.629 0.928
0.101 0.247 0.497 0.638 0.957
0.0023 0.0056 0.0097 0.0123 0.0147
0.0031 0.0042 0.0106 0.01315 0.01856
Towell and Ackermann (1972)
0.021 (0.020) 0.023 (0.022) 0.09 (0.088)
0.053 (0.05) 0.06 (0.06) 0.178 (0.18)
1.958 1.754 1.548
0.347 0.365 0.391
0.086 0.091 0.351
0.091 0.106 0.364
0.0412 0.0670 0.0930
0.0489 0.0668 0.0977
Hikita and Kikukawa (1974)
0.069 (0.07) 0.098 (0.10) 0.158 (0.16) 0.22 (0.23) 0.33 (0.32) 0.048 (0.05) 0.074 (0.075) 0.11 (0.10) 0.119 (0.12)
0.15 (0.152) 0.176 (0.175) 0.20 (0.207) 0.264 (0.267) 0.32 (0.317) 0.141 (0.143) 0.166 (0.165) 0.179 (0.180) 0.21 (0.212)
2.462 2.224 2.064 1.983 1.864 2.546 2.212 2.043 1.998
0.257 0.274 0.367 0.401 0.423 0.236 0.267 0.358 0.361
0.377 0.562 0.876 1.168 1.569 0.239 0.387 0.474 0.534
0.367 0.549 0.902 1.089 1.468 0.251 0.394 0.451 0.548
0.0302 0.0352 0.0411 0.0468 0.0542 0.0312 0.0358 0.0387 0.0423
0.0286 0.0331 0.0413 0.0473 0.0518 0.0321 0.0360 0.0365 0.0415
Krishna et al. (2000), D ˆ 0.174 m
0.033 (0.034) 0.092 (0.09) 0.158 (0.16) 0.22 (0.23) 0.27 (0.27) 0.31 (0.30)
0.090 (0.087) 0.168 (0.165) 0.21 (0.218) 0.256 (0.250) 0.265 (0.262) 0.268 (0.270)
3.398 2.948 2.643 2.421 2.182 1.986
0.286 0.306 0.321 0.372 0.423 0.486
0.164 0.477 0.885 1.260 1.452 1.566
0.182 0.513 0.898 1.326 1.483 1.624
0.0252 0.0291 0.0356 0.0401 0.0443 0.0498
0.0245 0.0278 0.0291 0.0377 0.0393 0.0411
van Baten and Krishna (2001), D ˆ 0.63 m Moustiri et al. (2001)
0.089 (0.09) 0.21 (0.23) 0.32 (0.30) 0.032 (0.0319) 0.047 (0.0472) 0.057 (0.0567)
0.075 (0.089) 0.254 (0.265) 0.288 (0.295) 0.076 (0.076) 0.125 (0.126) 0.155 (0.154)
2.846 2.392 1.992 2.467 2.146 2.013
0.323 0.382 0.476 0.347 0.371 0.333
0.477 1.260 1.566 0.137 0.150 0.166
0.513 1.326 1.624 0.148 0.161 0.173
0.1746 0.2483 0.2743 0.0032 0.0062 0.0086
0.1283 0.1674 0.2043 0.0051 0.0059 0.0091
Deckwer et al. (1974)
The bracketed numbers indicate the experimental values; LHS is volume integral of energy dissipation rate and RHS is equation (16).
observed that, at a particular VG, the value of DL increases with an increase in column diameter with a proportionality in the range of D1.2–1.6. The validity of the CFD model for axial dispersion coef cient was established. Therefore, comparison has been made with experimental data of Ohki and Inoue (1970), Towell and Ackermann (1972), Hikita and Kikukawa (1974), Deckwer et al. (1974), Krishna et al. (2000), Moustiri et al. (2001) and Khare (1990). Table 1 gives the detail pertaining to column diameter, column height, gas– liquid system, range of VG, tracer and the measurement techniques employed for the measurement of axial dispersion coef cient. It may be noted that practically all the major experimental investigations published in the literature have been covered in the table. The results are given in Table 10 and Figure 6. The agreement between the predicted and the experimental data can be seen to be excellent. The agreement over the wide range of D, VG and the gas–liquid systems ensures the applicability of the model (for the estimation of axial dispersion coef cient). Now we would like to bring out one observation regarding the relative importance of convection and eddy diffusivity (nt) on axial dispersion coef cient (DL). Table 11 shows the values of DL and average nt in the column. It can be seen that, in all the cases of VG and D, the value of DL is at least 10 times higher than that of nt. This observation indicates the overriding contribution mean circulatory ow on axial mixing and mixing as compared with that of nt. The sensitivity analysis of the effects of mean ow and nt Trans IChemE, Vol 81, Part A, September 2003
(results not shown) con rms this observation. In fact, the axial dispersion model assumes the mean ow to be unidirectional ow with a at pro le with axial dispersion superimposed. The ideal situation is one directional at pro le with eddy diffusivity superimposed. In comparison with such a ideal situation, the real dispersed plug ow model accounts for any nonidealities from the at velocity pro le such as bypass, recycle, channeling, dead zones, etc.
Figure 6. Parity plot for axial dispersion coef cient. ( ) Ohki and Inoue (1970); ( ) Towell and Ackermann (1972); ( ) Hikita and Kikukawa (1974); ( ) Deckwer et al. (1974); ( ) Khare (1990); ( ) Krishna et al. (2000); ( ) Moustiri et al. (2001); and ( ) van Baten and Krishna (2001).
1000
EKAMBARA and JOSHI Table 11. Comparison between CFD predictions and the experimental observation (Khare, 1990).
VG (m s¡1)
Material balance
Energy balance
Column diameter (m)
E· G
C0
C1
LHS
RHS
nt (m2 s¡1)
Dispersion coef cient (DL)
0.20 0.20 0.20 0.40 0.40 0.40
0.136 (0.135) 0.223 (0.223) 0.304 (0.282) 0.133 (0.135) 0.227 (0.223) 0.311 (0.282)
2.17 1.99 1.76 2.155 2.135 1.814
0.38 0.43 0.51 0.367 0.399 0.431
0.362 0.821 1.292 0.399 0.877 1.374
0.329 0.815 1.488 0.329 0.856 1.487
0.0018 0.0026 0.0032 0.0048 0.0067 0.0074
0.0287 0.0386 0.0469 0.0660 0.0956 0.1196
0.072 (0.07) 0.173 (0.17) 0.292 (0.295) 0.073 (0.07) 0.173 (0.17) 0.302 (0.295)
The bracketed numbers indicate the experimental values; LHS is volume integral of energy dissipation rate and RHS is equation (16).
Table 12. Axial dispersion coef cient (m2 s 1) VG (m s¡1)
VL (m s¡1)
HD=D
E· G
Column diameter ˆ 0.2 m 0.070 0.170 0.295
0.0066 0.015 0.035
0.135 0.223 0.282
Column diameter ˆ 0.4 m 0.070 0.170 0.295
0.0066 0.015 0.035
0.135 0.223 0.282
However, such nonidealities are expected to be within limits for the applicability of dispersed plug ow model. In fact, the average recirculation ow (recycle) in bubble columns is 10–1000 times higher than the mean net ow across the column. Such a major deviation is an alerting signal for checking the applicability of dispersed plug ow model for bubble columns. Let us review some of the characteristic features of the dispersed plug ow model visa-vis the performance of a bubble column. (i) For a dispersed plug ow model, the relationship between the mixing time, column height and DL is given by the following equation (Towell and Ackermann, 1972): ymix ˆ
aHD2 DL
(22)
3
5
8
10
0.0034 0.0139 0.0209
0.0125 0.0225 0.0340
0.0198 0.0316 0.0423
0.0287 0.0386 0.0469
2
3
4
5
0.0166 0.0345 0.0366
0.0260 0.0485 0.0657
0.0409 0.0690 0.0926
0.0660 0.0956 0.1196
shown in Figures 4 and 5. The CFD-predicted results of axial dispersion coef cient are given in Table 12. From Table 12, it can be observed that the axial dispersion coef cient increases with an increase HD=D ratio for both the columns. The proportionality relationship was found to be DL / HD0:65±1:72 . Since the effect of HD on DL has not been considered in the past, all the published correlations do not consider HD as a variable parameter. In order to understand the performance of the published correlations, the values of DL were calculated using different correlations. Figure 7 shows wide differences in the correlations. It is believed that the source of these differences is nonconsideration of the effect of HD on DL. Therefore, for all these published investigations, the CFD approach (of this paper) was employed and the predicted
where the value of constant a depends upon the extent of homogeneity for which the mixing time is measured. It may be noted that the mixing time is proportional to the square of the dispersion height (HD) assuming that DL is independent of HD. However, the experimental observations of mixing time (Table 9) conform to ymix HD0:40±0:66 . This experimental is quantitatively supported by the CFD (the predicted ymix is given in Table 9). (ii) As mentioned earlier, the axial dispersion coef cient is considered to be independent of dispersion height. In fact, this relationship has not been systematically investigated in the published literature. Therefore, it was thought desirable to probe into such a relationship using computational uid dynamics. The effect of height to diameter ratio (1–10) on axial dispersion coef cient was studied for column diameters 0.2 and 0.4 m and super cial gas velocity of 0.07– 0.295 m s¡1. The residence time distribution curves are
Figure 7. Dispersion coef cient from different correlations for column diameter 0.2 m. 1, Towell and Ackermann (1972); 2, Hikita and Kikukawa (1974); 3, Shah et al. (1978); 4, Joshi (1980); 5, Field and Davidson (1980); 6, Khare (1990); and 7, Krishna et al. (2000).
Trans IChemE, Vol 81, Part A, September 2003
CFD SIMULATION OF MIXING AND DISPERSION values were shown in Figure 6. The excellent agreement with all the published information shows the usefulness of CFD approach over the semi-empirical correlations.
FDz FL FVr
CONCLUSION A two-dimensional CFD code has been developed for two-phase gas–liquid dispersion in bubble columns. A detailed comparison has been presented between the CFD simulation and the experimental data reported by Hills (1974), Nottenkamper et al. (1983), Menzel et al. (1990), Yao et al. (1991), Yu and Kim (1991) and Grienberger and Hofmann (1992). Excellent agreement was observed between the predicted and the measured pro les of axial liquid velocity for all investigations. A complete energy balance was also established for all the cases. The ow model has been extended to the process of homogenization of an inert tracer. A systematic attempt has been made to investigate the effect of tracer addition and detector position and the column diameter, super cial gas velocity and height to diameter ratio on mixing time. A good agreement was observed between the predicted and the experimental values of mixing time for two column diameters (0.2 and 0.4 m i.d.) and three super cial gas velocities (0.07, 0.17 and 0.295 m s¡1). Further, the model has been extended for the estimation of axial dispersion coef cient (DL). Excellent agreement was found between the predicted and the experimental values of the axial dispersion coef cient measured over a wide range of D, VG and HD=D covered by various investigators in the past four decades. The parametric sensitivity of dispersion height has been investigated and in this regards, the usefulness of the CFD approach has been brought out as compared with the conventional semi-empirical correlations. NOMENCLATURE a b, c, d c c0 c· c0 u0 c0 v0 C Ci CB CD CL CV C0, C1 Ce1 Ce2 D Deff DL Dm Dt Dtr Dtz eR E ES fD FDr
constant in equation (22) empirical constant in equation (1) instantaneous concentration of the tracer uctuating concentration of the tracer time-averaged concentration of the tracer term de ned in equation (6) term de ned in equation (6) mean concentration of the tracer, de ned in equation (17) initial concentration of the tracer interface energy transfer factor drag force coef cient lift force coef cient virtual mass force coef cient drift ux constants model parameter in turbulent dissipation energy equation (ˆ1.44) model parameter in turbulent dissipation energy equation (ˆ1.92) diameter of the column, m effective diffusion coef cient, m2 s¡1 axial dispersion coef cient, m2 s¡1 molecular diffusivity, m2 s¡1 eddy diffusivity, m2 s¡1 radial eddy diffusivity, m2 s¡1 axial eddy diffusivity, m2 s¡1 rate of energy release for a single bubble, J s¡1 rate of energy released by all the bubble in the column, J s¡1 fraction of dissipation energy, J s¡1 drag force on a single particle in uidized bed, N m¡2 frictional force in the radial direction per unit volume of dispersion, N m¡3
Trans IChemE, Vol 81, Part A, September 2003
1001
frictional force in the axial direction per unit volume of dispersion, N m¡3 lift N m¡3 £ force per ¡unit volume ¢ of dispersion, ¤ ˆCL EL EG rL uG ¡ uL @uL =@r radial virtual mass force per unit volume of dispersion, N m¡3 µ » ¼ ¶ 1@ @ ˆCV EL EG rL r(vG ¡ vL ) ‡ (vG ¡ vL ) r @r @z
FVz
axial virtual mass force per unit volume of dispersion, N m¡3 µ » ¼ ¶ 1@ @ ˆCV EL EG rL r(uG ¡ uL ) ‡ (uG ¡ uL ) r @r @z
g
gravitational constant, m s¡2 µ Á» ¼ 2 n o » ¼ @vL v 2 @uL ˆ mt,L 2 ‡ L ‡ @r r @z
G H HD k n P PB Pe r R Re SF,K t u uG uL u· u0 v vB vG vL v· v0 VC VG VL VL(0) VSr VS z
2
! ³» ¼ » ¼ ´2 ¶ @vL @uL ‡ ‡ @z @r
height of the bubble column, m height of the gas dispersion, m turbulent kinetic energy per unit mass, m2 s¡2 empirical constant in equation (1) pressure, N m¡2 CB[FDrVSr ‡ FDzVS] Peclet number, (DL=VLHD) radial distance, m column radius, m Reynolds number source term of phase K time, s instantaneous axial velocity, m s¡1 axial component of gas velocity, m s¡1 axial component of liquid velocity, m s¡1 axial time-averaged velocity, m s¡1 axial uctuating velocity, m s¡1 instantaneous radial velocity, m s¡1 volume of single bubble, m3 radial component of gas velocity, m s¡1 radial component of liquid velocity, m s¡1 radial time-averaged velocity, m s¡1 radial uctuating velocity, m s¡1 circulation velocity, m s¡1 super cial gas velocity, m s¡1 super cial liquid velocity, m s¡1 center line liquid velocity, m s¡1 radial slip velocity between gas and liquid, m s¡1 axial slip velocity between gas and liquid, m s¡1 axial distance along the column, m
Greek symbols mK ‡ mt,K=sF GK a under relaxation parameter E fractional phase holdup EG fractional gas holdup EL fractional liquid holdup E· L time-averaged fractional liquid holdup E0L uctuating fractional liquid holdup e turbulent energy dissipation rate per unit mass, m2 s¡3 mK molecular viscosity of phase K, Pa s mt,K turbulent viscosity of phase K, Pa s n molecular kinematic viscosity of liquid, m2 s¡1 nt turbulent kinematic viscosity of liquid, m2 s¡1 ymix mixing time, s r density, kg m¡3 rG density of the gas, kg m¡3 rL density of the liquid, kg m¡3 s2 variance, de ned in equation (19), s2 sF turbulent Prandtl number for momentum transfer sf turbulent Prandtl number for bubble motion s2y dimensionless variance, de ned in equation (20) t mean residence time, de ned in equation (18), s Subscripts G K L
gas phase phase, K ˆ G, gas phase; L, liquid phase liquid phase
1002
EKAMBARA and JOSHI REFERENCES
Craver, M.B., 1984, Numerical computation of phase separation in two phase ow, J Fluids Eng, 106: 147–153. Deckwer, W.D., Burckhart, R. and Zoll, G., 1974, Mixing and mass transfer in tall bubble columns, Chem Eng Sci, 29: 2177–2188. Deshpande, N.S., Prasad, Ch.V., Kulkarni, A.A. and Joshi, J.B., 2000, Hydrodynamic characterization of dispersed two-phase ows by laser Doppler velocimeter, Trans IChemE, Part A, Chem Eng Res Des 78: 903–910. Field, R.W. and Davidson, J.F., 1980, Axial dispersion in bubble columns, Trans IChemE, 58: 228–236. Grienberger, J. and Hofmann, H., 1992, Investigation and modeling of bubble columns, Chem Eng Sci, 47: 2215–2220. Hillmer, G., Weismantel, L. and Hofmann, H., 1994, Investigations and modelling of slurry bubble columns, Chem Eng Sci, 49: 837–843. Hills, J.H., 1974, Radial nonuniformity of velocity and voidage in the bubble column, Trans IChemE, 2: 1–9. Hikita, H. and Kikukawa, H., 1974, Liquid phase mixing in bubble columns: effect of liquid properties, Chem Eng J, 8: 191–197. Hjertager, B.H. and Morud, K., 1993, Computational uid dynamics simulation of bioreactors, CFD Simul, 47–61. Johansen, S.T. and Boysan, F., 1988, Fluid dynamics in bubble stirred ladles: Part II. Mathematical modelling, Metall Trans B, 19B: 755– 764. Joshi, J.B., 1980, Axial mixing in multiphase contactors—a uni ed correlation, Trans IChemE, 58: 155–165. Joshi, J.B., 2001, Computational ow modelling and design of bubble column reactors, Chem Eng Sci, 56: 5893–5933. Joshi, J.B. and Sharma, M.M., 1979, A circulation cell model for bubble columns, Trans IChemE, 57: 244–251. Kataoka, I., Besnard, D.C. and Serizawa, A., 1992, Basic equation of turbulence and modelling of interfacial transfer terms in gas–liquid two-phase ow, Chem Eng Commun, 118: 221–236. Kato, Y., Nishiwaki, A., Fukuda, T. and Tanaka, S., 1972, The behaviour of suspended solid particles and liquid in bubble columns, J Chem Eng Japan, 5: 112–118. Khare, A.S., 1990, Development of radio active tracer techniques and its application to the design of multiphase reactors, Ph.D. thesis, University of Mumbai, India. Krishna, R., Urseanu, M.I., van Baten, J.M. and Ellenberger, J., 2000, Liquid phase dispersion in bubble columns operating in the churn-turbulent ow regime, Chem Eng J, 78: 43–51. Levenspiel, O., 1999, Chemical Reaction Engineering, 3rd edition (Wiley, New York). Menzel, T., Weide, T., Staudacher, O. and Onken, U., 1990, Reynolds stress model for bubble column reactor, Ind Eng Chem Res, 29: 988– 994. Moustiri, S., Hebrard, G., Thakre, S.S. and Roustan, M., 2001, A uni ed correlation for predicting liquid axial dispersion coef cient in bubble columns, Chem Eng Sci, 56: 1041–1047. Nottenkamper, R., Steiff, A. and Weinspach, P.M., 1983, Experimental investigation of hydrodynamics of bubble columns, Ger Chem Eng, 6: 147–155. Ohki, Y. and Inoue, H., 1970, Longitudinal mixing of the liquid phase in bubble columns, Chem Eng Sci, 25: 1–16. Patankar, S.V., 1981, Numerical Heat Transfer and Fluid Flow (McGrawHill, New York, USA). Patankar, S.V. and Spalding, D.B., 1972, A calculation procedure for heat, mass and momentum transfer in three dimensional parabolic ows, Int J Heat and Mass Transfer, 15: 126–129.
Peric, M., Kessler, R. and Scheuerer, G., 1988, Comparison of nite volume numerical methods with staggered and collocated grids, Comput Fluids, 16: 389–403. Ranade, V.V., 1992, Flow in bubble columns: some numerical experiments, Chem Eng Sci, 47(8): 1857–1869. Ranade, V.V., 1997, Modelling of turbulent ow in a bubble column reactor, Trans IChemE, Part A, Chem Eng Res Des, 75: 14–23. Reith, T., Renken, S. and Israel, B.A., 1968, Gas holdup and axial mixing in the uid phase of bubble columns, Chem Eng Sci, 23: 619. Riquarts, H.P., 1981, A physical model for axial mixing of the liquid phase for heterogeneous ow regime in bubble columns, Ger Chem Eng, 4: 18–23. Sato, Y. and Sekoguchi, K., 1975, Liquid velocity distribution in two-phase bubble ow, Int J Multiphase Flow, 2: 79–95. Shah, Y.T., Stiegel, G.J. and Sharma, M.M., 1978, Back-mixing in gas– liquid reactors, AIChEJ, 24: 369–400. Spalding, D.B., 1980, Mathematical Modelling of Fluid Mechanics, Heat Transfer and Chemical Reaction Processes, Imperial College Report HTS=80=1, Canbay, UK. Svendsen, H.F., Jakobsen, H.A. and Torvik, R., 1992, Local ow structures in internal loop and bubble column reactors, Chem Eng Sci, 47: 3297–3304. Thakre, S.S., Phanikumar, D.V., Khare, A.S. and Joshi, J.B., 1999, CFD modeling of ow, macro-mixing and axial dispersion in bubble column, Can J Chem Eng, 77: 826–837. Theofanous, T.G. and Sullivan, J., 1982, Turbulence in two-phase dispersed ows, J Fluid Mech, 116: 343–362. Towell, G.D. and Ackermann, G.H., 1972, Axial mixing of liquids and gas in large bubble reactor, in Proceeding of 2nd International Symposium on Chemical Reaction Engineering, Amsterdam, The Netherlands, B3.1–B3.13. van Baten, J.M. and Krishna, R., 2001, Eulerian simulations for determination of the axial dispersion of liquid and gas phases in bubble columns operating in the churn-turbulent regime, Chem Eng Sci, 56: 503–512. Van Doormal, J.P. and Reithby, G.D., 1984, Enhancement of the SIMPLE method for predicting incompressible uid ows, Numer Heat Transfer, 7: 147–163. Yao, B.P., Zheng, C., Gasche, H.E. and Hoffmann, H., 1991, Bubble behaviour and ow structure of bubble columns, Chem Eng Proc, 29: 65–75. Yu, Y.H. and Kim, S.D., 1991, Bubble properties and local liquid velocity in the radial direction of cocurrent gas–liquid ow, Chem Eng Sci, 46: 313–320.
ACKNOWLEDGEMENTS One of us (Ekambara) acknowledges the fellowship support given by the University Grants Commission (UGC), Government of India.
ADDRESS Correspondence concerning this paper should be addressed to Professor J. B. Joshi, Institute of Chemical Technology, University of Mumbai, Matunga, Mumbai-400 019, India. E-mail: [email protected] The manuscript was received 5 August 2002 and accepted for publication after revision 23 May 2003.
Trans IChemE, Vol 81, Part A, September 2003