Chemical Engineering Science 57 (2002) 3157 – 3183

www.elsevier.com/locate/ces

Review Paper

Coherent "ow structures in bubble column reactors J. B. Joshi ∗ , V. S. Vitankar, A. A. Kulkarni, M. T. Dhotre, K. Ekambara Institute of Chemical Technology, University of Mumbai, Matunga, Mumbai 400 019, India Received 26 November 2001; received in revised form 7 January 2002; accepted 14 January 2002

Abstract This paper reviews the "ow patterns in bubble columns with a focus on transient "ow structures. The subject of mean "ow pattern has not been covered in view of the earlier publication (Joshi, J.B. Computational "ow modelling and design of bubble column reactors. Chemical Engineering Science 56 (2001) 5893–5933), which also deals with the subject of the formulation of governing equations and the description of interface force terms for gas–liquid dispersions. The published literature in the last ten years has been analyzed on a coherent basis and the present status has been brought out for (i) Euler–Euler versus Euler–Lagrange approach, (ii) closure formulation for eddy di;usivity in gas–liquid dispersions, (iii) numerical issues such as grid size, discretization scheme and the solution algorithms and (iv) two-dimensional versus three-dimensional formulation. The following geometries of bubble columns have been investigated in the past: (i) two-dimensional column with o;-center sparging, (ii) two-dimensional column with central sparging, (iii) two-dimensional column with uniform sparging and (iv) cylindrical columns with uniform sparging. All these cases have been covered in the present paper. A critical account of the predictive capability of all the models for transient "ow structures has been presented. This paper also includes the application of multiresolution analysis of velocity–time series for the identi>cation of coherent "ow structures. ? 2002 Elsevier Science Ltd. All rights reserved. Keywords: Bubble column; Computational "uid dynamics; Coherent structures; Multiresolution analysis

Contents 1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2. Mathematical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3. Flow driven by a bubble plume . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1. O;-center single-point sparging (two-dimensional column) . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1. Becker et al. (1994) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2. Sokolichin et al. (1997) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.3. Delnoij et al. (1997b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.4. Mudde and Simonin (1999) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.5. Sokolichin and Eigenberger (1999) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2. Single-point central sparging (two-dimensional column) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1. Borchers et al. (1999) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2. P"eger et al. (1999) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3. Becker et al. (1999) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. Uniform sparging in two-dimensional columns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1. Experimental observations of Chen, Jamialahmadi, and Li (1989) . . . . . . . . . . . . . . . . . . . . . 4.2. Lapin and Lubbert (1994) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3. Sokolichin and Eigenberger (1994) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4. Sokolichin et al. (1997) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5. Deen et al. (2001) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

∗

Corresponding author. Tel.: +91-22-414-5616; fax: +91-22-414-5614. E-mail address: [email protected] (J. B. Joshi).

0009-2509/02/$ - see front matter ? 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 0 9 - 2 5 0 9 ( 0 2 ) 0 0 1 9 2 - 6

3158 3162 3166 3166 3166 3168 3168 3169 3170 3171 3172 3173 3174 3174 3174 3174 3175 3175 3175

3158

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

5. Cylindrical columns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1. Experimental observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.1. Devanathan, Moslemian, and Dudukovic (1990) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.2. Drahos, Zahradnik, Fialova, and Bradka (1992) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.3. Tzeng, Chen, and Fan (1993) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.4. Chen, Reese, and Fan (1994) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.5. Reese and Fan (1994) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.6. Lin, Reese, Hong, and Fan (1996) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2. Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.1. Lapin and Lubbert (1994) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.2. Millies and Mewes (1994) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.3. Grevskott et al. (1996) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.4. Becker et al. (1999) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.5. P"eger and Becker (2001) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6. Application of multiresolution analysis of velocity–time data for the characterization of . . . . . . . . turbulent "ow structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1. Preamble . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7. Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1. Introduction Bubble columns are widely used in industry because of their simple construction and operation. Important applications include oxidation, hydrogenation, halogenation, hydrohalogenation, ammonolysis, hydroformylation, Fischer– Tropsch reaction, ozonolysis, carbonylation, carboxylation, alkylation, fermentation, waste water treatment, hydrometallurgical operations, steel ladle stirring, column "otation, etc. In bubble columns, the gas phase exists as a dispersed bubble phase in a continuous liquid phase. The gas phase moves in one of the two characteristic regimes depending upon the nature of dispersion. The two regimes are: homogeneous and heterogeneous. The homogeneous regime occurs at relatively low super>cial gas velocities (less than about 50 –80 mm=s). Almost uniformly sized bubbles characterize this regime and the concentration of bubbles is uniform, particularly in the transverse direction. The heterogeneous regime occurs at relatively high super>cial gas velocities. This regime is characterized by the presence of a radial hold-up pro>le against a "at pro>le in the homogeneous regime. The gas hold-up pro>le results in a pro>le of static pressure, which is lower in the central region as compared to the near-wall region. As a consequence, intense liquid circulation is developed which is upward in the central region and downwards near the column wall as shown in Fig. 1A. The liquid circulation velocities are of the order of 0.1–3:0 m=s depending upon the super>cial gas velocity and the column diameter. However, the depicted "ow pattern is the average "ow pattern over a long time. The instantaneous "ow pattern consists of a spectrum of "ow structures. The largest ones appear as circulation cells (Joshi, 1980; Zehner, 1982). These are of di;erent sizes and shapes and usually do not remain stationary in time and space. These cells form one or two arrays in the vertical direction (as seen in the

3175 3175 3175 3176 3176 3176 3177 3177 3177 3177 3177 3177 3177 3178 3178 3178 3180 3181

cross-sectional view), the latter being more common and appear as staggered or unstaggered arrays. The maximum size of circulation cell is of the order of column diameter and can be considered as the largest eddy in the column. The smaller eddies are also present covering a very wide range of wave number typically represented in an energy spectrum (Kolmogorov, 1941). The intensity of turbulence in bubble columns is much higher (15 – 40%, Kulkarni, Joshi, Ravi Kumar, & Kulkarni, 2001b) as compared to that in the single-phase pipe "ow (3–10%, Davies, 1972). Though the bubble columns are simple in construction consisting of only two major elements (the cylindrical vessel and the gas distributor), the design procedures are still not reliable. The most important hurdle has been the lack of understanding of the underlying physics of "ow and hence the quantitative relationships between the reactor parameters (hardware and operating) and the performance objectives. Therefore, the available relationships are empirical and cannot be considered reliable over a wide range of applications encountered in practice. The combination of such empirical relationships and the accumulated operating experience has been the present state of the art rather than the most desired state of the science. For this transformation, vigorous attempts have been made during the past 25 years by focussing on the physics of "uid "ow by using two tools, namely experimental "uid dynamics (EFD) and the computational "uid dynamics (CFD). There have been two principal objectives: (i) the prediction of "ow pattern in terms of three-dimensional pro>les of holdup, axial and radial velocities of gas and liquid phases, eddy di;usivity, mean and turbulent kinetic energy (K and k), mean and turbulent energy dissipation rates, Reynolds stresses, etc. In addition, attempts have also been made to predict eddy size distribution and eddy velocity distribution; (ii) to develop relationships between the "ow pattern and the design objectives such as mixing time, residence time distribution

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

Fig. 1. (A) Typical gas hold-up and static pressure pro>les in a bubble column. (B) Schematic representation of the circulation "ow pattern in bubble column.

(RTD), bubble size distribution, gas–liquid e;ective interfacial area and mass transfer coeLcient, two-phase pressure drop, wall heat transfer coeLcient and wall mass transfer coeLcient. CFD de>nitely has great potential (Technology Vision 2020, 1996). However, the present status of CFD (which includes all commercial and the known homemade codes) is primitive, though very promising. In the >rst stage of the development of the CFD codes, the entire attention was given for the development of the mean pro>les, and the next stage of transient analysis is being followed by the progress in computational power. The mean "ows in bubble columns have been extensively investigated in the published literature and reviewed by Stewart and Wendro; (1984), Jakobsen, Sannaes, Grevskott, and Svendsen (1997) and Joshi (2001). However, these review papers have mainly considered the average "ow pattern practically in the form of a single circulation cell, which is schematically shown in Fig. 1B. The up"ow occurs in the central region, down"ow near the column wall and the "ow reversals take place at the bottom and top of the column. Therefore, in a major portion of the column, the "ow is considered to be essentially axial and the radial "ow is seen only at the ends of the column. Such a one-cell "ow pattern is known to be unstable for shallow (HD =D ¡ 0:5) and tall columns (HD =D ¿ 1:0). In the former case, multiple cells occur in the transverse direction whereas in the latter case, multiple cells occur in the axial direction. In practice, higher HD =D ratios are common. For these cases, Joshi and Sharma (1979) and Joshi (1980) proposed multiple circulation cell models and

3159

Fig. 2. Schematic representation of multiple circulation cells in bubble columns after (A) Joshi and Sharma (1979), (B) Joshi (1980) and (C) Zehner (1982).

these are schematically shown in Fig. 2. Fig. 2A shows non-interacting cells. Fig. 2B shows interacting cells with considerable intercirculation. Further, Joshi (1982, 1992) has pointed that the cell structure as shown in Fig. 2A and B is instantaneous "ow structure and does not remain stationary in space and time. The second model satisfactorily predicted the liquid velocity pro>les and the fractional gas holdup in the column diameter range of 0.14 –5:5 m and the super>cial gas velocity in the range of 19 –920 mm=s. The second model also developed a rational correlation for axial mixing in bubble columns. It also developed a uni>ed correlation for variety of multiphase reactors (Joshi, 1980). Zehner (1982) also found it necessary to assume the secondary "ow structures for favorable comparison with the experimental liquid velocity pro>les. The arrangement of the cells is shown in Fig. 2C. Lee and Korpela (1983) have analyzed the problem of natural convection in a two-dimensional slot. In this case, the driving force for circulation was the density di;erence arising out of the temperature pro>le. This problem is similar to bubble columns where the density di;erence is generated by the gas hold-up pro>le. Lee and Korpela (1983) have solved the equations of motion together with the equation of energy. They have shown that the multiple cells occur when HD =D ratio and the density di;erence exceed a certain limit. Interestingly, some of the "ow patterns obtained by Lee and Korpela (1983) are identical to those proposed by Joshi (1980). A large number of the published evidence (Wu & Gidaspow, 2000) also support the "ow pattern suggested by Joshi (1980). These authors have pointed out that the multiple circulation cells of the type proposed by Joshi and Sharma (1979) are

3160

Table 1 Summary of previous work: experimental details, formulation of the interface force terms and remarks

Author(s)

Column and sparger design

Measurement techniques

Forces

Limitations

Drag

Lift

Virtual

Lapin and Lubbert (1994)

Rectangular, Cylindrical 4:0 m × 1:0 m; 1:5 m × 1:0 m

Uniform point

—

(a)

NC

NC

1, 2, 3, 4, 6c, 7, 8, 9, 10, 11, 12, 15, 16, 17, 18, 19, 20, 21, 22, 23

Sokolichin and Eigenberger (1994)

Rectangular 0:75 m × 0:15 m

Uniformly aerated

—

Eq. (2.7)

NC

NC

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

Becker et al. (1994)

Rectangular 2:0 m × 0:5 m × 0:08 m

Firt sparger, tube sparger

LDA Double-channel >ber optical probe

Eq. (2.7)

NC

NC

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

Sokolichin et al. (1997)

Rectangular 2:0 m × 0:5 m × 0:08 m

Single point, sparger and uniformly aerated

—

(a)

NC

NC

2, 3, 4, 6a, 7, 8, 9, 10, 11, 12, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23

Delnoij et al. (1997b)

Rectangular 1:5 m × 0:5 m × 0:05 m

Single point, located 15 cm from left wall

LDA and microturbine

(b)

(c)

(d)

1, 2, 4, 6a, 7, 8, 9, 13, 14, 16, 17, 18, 19, 20, 22, 23

Delnoij et al. (1997a)

Rectangular 2:0 m × 0:25 m × 0:05 m; 2:0 m × 0:25 m × 0:02 m

Single point,

Video imaging technique, streak photography

(b)

(c)

(d)

1, 2, 4, 6a, 7, 8, 9, 13, 14, 16, 17, 18, 19, 20, 22, 23

Delnoij (1999)

Rectangular 1:5 m × 0:5 m × 0:05 m

Single point, located 15 cm from left wall

LDA and microturbine

(b)

(c)

(d)

1, 2, 4, 6a, 7, 8, 9, 13, 14, 15, 16, 17, 18, 19, 20, 22, 23

Sokolichin and Eigenberger (1999)

Rectangular 2:0 m × 0:5 m × 0:08 m

Single point, sparger

—

(a)

NC

NC

2, 3, 4, 5a, 5b, 5c, 6a, 7, 8, 9, 10, 11, 12, 14, 15, 16, 17, 18,19, 20, 21, 22, 23

Borchers et al. (1999)

Rectangular 2:0 m × 0:5 m × 0:08 m

Frit sparger

LDA

(a)

NC

NC

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

Mudde and Simonin (1999)

Rectangular 2:0 m × 0:5 m × 0:08 m

Single point, (circular sparger)

—

(e)

NC

(f)

1, 2, 3, 5b, 5c, 6b, 7, 8, 9, 12, 15, 16, 17, 18, 19, 20, 21, 22, 23

P"eger et al. (1999)

Rectangular 0:45 m × 0:2 m × 0:05 m

Single point, (rectangular)

—

(g) CD = 0:66

NC

NC

1, 2, 5b, 5c, 6b, 7, 8, 9, 11, 12, 15, 16, 17, 18, 19, 20, 21, 22, 23

P"eger and Becker (2001)

Cylindrical 2:6 m × 0:288 m

Uniform sparger, ring sparger

LDA

(g) CD = 0:44

NC

NC

1, 2, 5b, 5c, 6b, 7, 8, 9, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23

Deen et al. (2001)

Square 1:0 m × 0:15 m × 0:15 m

Uniform sparger

PIV

(g) CD = 1:0

(h)

(i)

1, 2, 5b, 5c, 6b, 7, 8, 9, 13, 15, 16, 17, 18, 19, 20, 21, 22, 23

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

Column type (H × W × D)

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

(a) Considered constant slip velocity (b) FD = 12 CD L R2b |uG − uL |(uG − uL ) CD = 24=Re(1 + 0:15Re0:687 ) Re ¡ 1000 =0:44 Re ¿ 1000 (c) FL = −CL L Vb (uG − uL ) × P P = ∇ ×
u  DI (d) FV = − + I · ∇u Dt I = CV L Vb (uG − uL )  CV = 0:5[1 + 2:78(1 − jL )] CD (e) FD = − 43 jG L ur |ur | db  |ur | = ur · ur + ur · ur G 24 CD = if Re ¡ 1000 (1 + 0:15Re0:687 )j−1:7 L Re j L L | ur | d Re = L   @ur (f ) FV = −CV jG L + uG · ∇ur + ∇ · jG L CV uG uL G @t 1 (g) FD = − 43 CD jG L |uG − uL |(uG − uL ) db (h) FL = −CL jG L (uG − uL ) × ∇ × uL CL = 0:5
 DuG DuL − (i) FV = −jG L CV Dt Dt CV = 0:5 NC : not considered: Assumptions made in the transient "ow simulations: (1) All bubbles have been assumed to be spherical. (2) All bubbles have been assumed to be of same size. (3) Gas phase momentum equations are not solved and a constant slip velocity has been assumed. (4) Liquid phase is assumed to be laminar. (5) Liquid phase is assumed to be turbulent: (a) Correlations of the velocity "uctuations are neglected in the equation of motion. (b) Correlations of the hold-up=velocity "uctuations are neglected in the equation of motion. (c) Correlations of the hold-up "uctuation of the two phases are neglected in the equation of motion. (6) (a) Drag coeLcient has been taken only for creeping "ow. (b) Drag coeLcient has been taken only for turbulent "ow. (c) Drag force was not formulated. (7) Drag coeLcient has been estimated for the single isolated bubble and presence of bubble swarm has been neglected. (8) Radial variation of the slip velocity has not been included. (9) Axial variation of slip velocity has not been included. (10) A value of slip velocity has been considered without proper formulation of drag force. (11) Virtual mass force has been neglected. (12) Lift force has been neglected. (13) The movement of the bubbles due to dispersion mechanism is not considered. (14) The formulation of buoyancy force needs re-examination or the formulation of buoyancy force was not involved. (15) Bubble–bubble interaction is not considered. (16) A coalescence phenomenon is not considered. (17) Gas–liquid mass transfer is not considered. (18) Overall energy balance is not established. (19) Bubble break-up phenomenon is not considered. (20) The di;erences between the various gas–liquid systems (because of physical properties and coalescing nature) are included. (21) A wide range of column diameter (D), height to diameter ratio (HD =D), and super>cial gas velocity (VG ) is not considered. (22) The e;ect of sparger design is not considered. (23) A complete correspondence between a real system and the predicted "ow pattern is not established.

3161

3162

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

common in the physics of waves (Elmore & Heald, 1969, p. 291). For example, starting with the Navier–Stokes equation, Tolstoy (1993) showed that standing waves in a narrow tube set up a system of steady vortices or closed loops, spaced a quarter wavelength apart. This type of behavior was observed experimentally more than a century ago in the well-known Kundt’s tube. Wu and Gidaspow (2000) have themselves computed such multiple circulation cells for the "ow in tall circulating "uidized beds. In view of the above observations, there has been a long-felt need to further our understanding of circulation cells and other coherent structures together with their dynamics. Such investigations received momentum when suLcient computation power was available. Therefore, after the mid-90s more than 15 papers have been published to investigate coherent "ow structures in bubble columns using CFD (Lapin & Lubbert, 1994; Sokolichin & Eigenberger, 1994; Becker, Sokolichin, & Eigenberger, 1994; Grevskott, Sannaes, Dudukovic, Hjarbo, & Svendsen, 1996; Sokolichin, Eigenberger, Lapin, & Lubbert, 1997; Delnoij, Lammers, Kuipers, & van Swaaij, 1997b,c; Delnoij, Kuipers, & van Swaaij, 1997a; Delnoij, 1999; Sokolichin & Eigenberger, 1999; Borchers, Busch, Sokolichin, & Eigenberger, 1999; Mudde & Simonin, 1999; P"eger, Gomes, Gilbert, & Wagner, 1999; Becker, De Bie, & Sweeney, 1999; Wu & Gidaspow, 2000; P"eger & Becker, 2001; Deen, Solberg, & Hjertager, 2001). The present paper focusses on the critical analysis of these published papers and makes a coherent presentation of the present status on the coherent structures and their dynamics. In addition to the CFD, the paper brie"y describes the application of multiresolution analysis of velocity–time series (measured at one or more locations in the column) for the identi>cation of coherent "ow structures. 2. Mathematical model The most important characteristic of a multiphase "ow is the existence of an interface separating the phases and the associated discontinuities of properties across the phase interface. The various transfer mechanisms between phases and between a two-phase mixture and a surrounding wall strongly depend upon the two-phase "ow regime. The two-phase "ows are generally classi>ed into two "ow regimes—separated "ows and dispersed "ows. In this review, we will focus on only dispersed "ows because of their most frequent occurrence in the chemical process industry. Transport phenomena in dispersed phase "ows depend upon the collective dynamics of bubbles interacting with each other and with the surrounding continuous phase. The descriptions of these interactions are still incomplete and are currently under progress. The present status pertaining to the formulation of governing equations and the interphase force terms have been given by Jakobsen et al. (1997); Jakobsen (2001) and Joshi (2001) have elaborated on the

interface force terms. Table 1 gives some formulations of the drag, lift, and virtual mass forces. The transient "ow patterns in bubble columns have been investigated by Lapin and Lubbert (1994), Sokolichin and Eigenberger (1994), Becker et al. (1994), Sokolichin et al. (1997), Delnoij et al. (1997b,c), Delnoij (1999), Sokolichin and Eigenberger (1999), Borchers et al. (1999), Mudde and Simonin (1999), P"eger et al. (1999), Becker et al. (1999), P"eger and Becker (2001) and Deen et al. (2001). These investigators have made several simplifying assumptions and these have been summarized in Table 1. It must be pointed out that the problem of transient "ow simulation is numerically an order of magnitude more diLcult than the steady-state "ow simulations. Therefore, for the transient "ow simulations, simplifying assumptions have been made in the past. These formulations have been described below in a chronological manner. Before embarking into the following subject of the transient "ow formulations, the readers may wish to see the derivation of governing equations for the mean "ow (Elghobashi & Abou-Arab, 1983; Joshi, 2001), description of interphase force terms (Jakobsen, 2001, Joshi, 2001), description of interphase energy transfer as well as overall energy balance (Thakre & Joshi, 1999). Lapin and Lubbert (1994) have proposed Euler–Euler and Euler–Lagrange models to predict the "ow structure in bubble columns. The proposed macro-model does not distinguish between the two phases in the gas–liquid medium and assumes a quasi-single-phase "ow system. Thus, the equations of continuity and momentum, respectively, were written for the gas–liquid phase as follows: @m + ∇ · (m um ) = 0; @t

(2.1)

@(m um ) (2.2) + ∇ · (m um um ) = ∇ ·  − ∇p + m g: @t The parameters, e;ective viscosity and density of the dispersion, required by the macro-models were calculated by a separate model, gas phase model. The authors have considered the gas phase in two ways, continuous Euler and discrete Lagrangian. In the continuous model, the gas phase was characterized by the bubble distribution density function w(r; m) in the phase space de>ned by the location vector and the bubble mass. The size has been assumed to remain constant in the domain. The governing equation in the continuous model has been written as
 dm @w @ + ∇ · (ub ; w) + w = 0: (2.3) @t @m dt The density in every control volume element was determined as a function of temperature, Tp and density function w: m = f(Tp ; w): The experimentally determined e;ective viscosity (in two-phase capillary "ows) has been used in the macro-model. In the discrete model for the gas phase, the in"uence of the

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

bubbles or the clusters on the density was evaluated using probability density function    Vk (r)fk (r) : (2.4) (r) = "uid 1 − k

Vk is the volume of bubbles in the cluster, fk (r) the probability density function characterizing the number of bubbles of class k per unit volume. The bell-shaped probability density function was assumed. The phenomena of the bubble breakup and coalescence that a;ect the bubble distribution density function was neglected in both the Eulerian and Lagrangian gas phase models. Due to the simpli>ed assumption of quasi-single-phase "ow, the problem of interface force formulation was eliminated. The values of slip velocity were assumed at several levels: 50, 100, 200 mm=s. The turbulence modeling was not considered and the experimentally determined e;ective viscosity was used in the macro-model. Such an assumption is a drastic simpli>cation and the estimated values of t are an order of magnitude di;erent than the real values observed in columns operated in a heterogeneous regime. Further, the e;ective viscosity has been assumed to be constant in the entire column. In addition, the important aspect of energy balance was also neglected. Sokolichin and Eigenberger (1994) have presented a considerably simpli>ed model. An extreme assumption of the laminar "ow was included in the model. The liquid phase continuity and momentum balance equations were written in Eulerian frame. The liquid phase continuity equation was simpli>ed by neglecting the variation of the volume fraction of the dispersed phase in time and space. The continuity and momentum equations thus get transformed into: ∇ · uL = 0;

(2.5)

@(jL L uL ) + ∇ · (jL L uL uL ) @t =jL ∇ ·  − jL ∇p + jL L g + FD :

(2.6)

They have considered the important interfacial drag force but have neglected the virtual mass and lift forces. Following Schwarz and Turner (1988), they have also assumed slip velocity to be constant (200 mm=s) and not related to the actual nature of gas–liquid systems. The drag force was de>ned by the following equation: FD = 5 × 104 jG (uG − uL ):

(2.7)

The dispersion mechanism in the gas phase was not considered and hence the dispersion term does not appear in the gas phase continuity equation: @(jG G ) + ∇ · (jG G uG ) = 0: @t

(2.8)

The quasi-steady form of the gas phase momentum balance was used. The inertia force, the lift force, the added mass force were neglected in the formulation of the dispersed

3163

phase and the simpli>ed gas momentum balance reduced to the following slip relationship: uG = uL + uslip :

(2.9)

No attempt has been made to incorporate the bubble breakup and coalescence phenomena in the gas phase momentum description. The model does not give any insight of energy transfer. Sokolichin et al. (1997) further simpli>ed the model of Sokolichin and Eigenberger (1994). They have assumed that, at low gas velocities, the dispersed phase elements do not change the overall density and hence neglected the momentum transfer due to the bubbles. The gas–liquid mixture density was calculated in the following simpli>ed manner:  = j L L :

(2.10)

Due to this simpli>cation, the dispersion continuity equation reduced to liquid phase continuity equation, which was given by @(jL L ) + ∇ · (jL L uL ) = 0: @t

(2.11)

It further avoided the formulation of momentum exchange terms, namely drag, virtual mass and lift force. Further simpli>cation was done in the momentum balance of gas–liquid dispersion by replacing dispersion velocity with the liquid phase velocity. The resulting equation of motion was given by @(L jL uL ) + ∇ · (L jL uL uL ) @t =∇ ·  − ∇p + jL L g:

(2.12)

It may be noted that due to this simpli>cation, the error was introduced in the right-hand side terms in Eq. (2.12). The pressure and stress tensor terms were made to in"uence liquid phase alone. The pressure was not assumed to be shared between the phases. The interfacial forces were also not formulated. The gas phase momentum balance was not solved instead the simpli>ed slip relation (Eq. (2.9)) was used to calculate the gas phase velocity. Assuming that the gas phase is incompressible, the slip velocity was prescribed externally as 200 mm=s. In addition to the continuity equation (Eq. (2.11)), a simpli>ed continuity equation was written for the gas phase:
 @jG @jG @ + ∇ · (jG uG ) = DG : (2.13) @t @xi @xi The dispersion coeLcient, DG , was estimated from the experimental data. The gas phase equation was solved by both the UPWIND and an implicit higher order total variation diminishing (TVD) schemes. The Eulerian simulations were compared with Lagrangian approach. In the Lagrangian approach, the gas phase was represented by a >nite number of dispersed gas phase particles (GPPs) and the gas holdup was calculated from the size and the position of the GPPs in the solution domain.

3164

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

The gas phase holdup obtained either by Eulerian or Lagrangian approach was used to update the local dispersion density in the liquid phase equations. The authors observed that the UPWIND scheme su;ered from serious numerical di;usion whereas implicit higher order TVD and LAGRANGE with large number of GPPs had similar long-time behavior. The authors have shown excellent agreement between the Euler–Euler approach (implicit higher order TVD) and the Lagrangian approach. However, the results of both the approaches did not compare favorably with the experimental data of Becker et al. (1994) for which the simulation was attempted. Delnoij et al. (1997b) followed the Euler–Lagrange approach to predict the vortical structures in the twodimensional bubble columns as observed by Becker et al. (1994). They also made an extreme assumption of liquid "ow pattern to be laminar. The remarkable feature of the work of Delnoij et al. (1997b) is their rigorous attempt to establish force balance for the bubble phase. The liquid phase continuity and momentum balance equations were written in Eulerian framework and are given by Eqs. (2.11) and (2.6), respectively. All the three forces, drag, virtual mass and lift force were considered. For the gas phase, Lagrangian approach was considered. The force balance around a single bubble was given by mb

duG = FD + F V + F L + F G + F P ; dt

(2.14)

where mb is the mass of a bubble. The >rst three terms on the right-hand side represent the three forces acting on the bubbles, namely drag, virtual mass, and the lift force. All the bubbles were considered to be spherical and hard. Following this assumption, the standard drag curve was used to specify the drag coeLcient: 24 (1 + 0:15 Re0:687 ); Re ¡ 1000; CD = Re (2.15) 0:44; Re ¿ 1000: The above equations are for rigid interfaces such as solid particles. Though, the general drag curve for bubbles could have been used, the focus of the work was on the completeness of the force balance. Delnoij et al. (1997a–c) have formulated positive lift force. However, experimental observations suggest radial inward movement of the bubbles. The virtual mass and the lift force coeLcients were assumed to be 0.5. One of the key features of the model is the establishment of complete force balance and the detailed bubble dynamics incorporated in it. The bubble–bubble interaction was considered in minute details. The processing of sequence of collisions was based on the method developed by Hoomans, Kuipers, Briels, and van Swaaij (1996). Sokolichin and Eigenberger (1999) have made a simplifying assumption that the two-phase "ow in bubble column can be considered as quasi-homogeneous single-phase "ow with averaged density given by m = jG G + jL L :

The continuity and momentum equations were written for the gas–liquid mixture as given by Eqs. (2.1) and (2.2), respectively. It may be noted that Sokolichin and Eigenberger (1994) considered the two "uid models and solved the complete momentum balance for the liquid phase whereas simpli>ed slip relation was solved for the gas phase. This simpli>cation in the gas phase momentum balance was continued by Sokolichin and Eigenberger (1999). The slip velocity was externally prescribed as 200 mm=s. To obtain the phase velocities and holdup, additional closures of gas phase continuity equation (2.8), the following equation representing the relationship between the velocities of both phases and the gas–liquid mixture were used: m um = jG G uG + jL L uL :

(2.16)

Due to the consideration of single-phase dispersion, the interfacial forces were not required. The second-order central di;erence "ux approximation was used to descretize the diffusion terms in the liquid phase equation and the convective terms were treated with near-second-order accurate TVD approach. The set of governing equations were solved for both the cases, laminar and turbulent. In the laminar approach, the e;ective viscosity of the dispersion was assumed to be equal to the liquid viscosity whereas in turbulent case, standard k–$ model was used for the liquid phase. The notable and remarkable feature of this work was three-dimensional simulations and its comparison with two-dimensional simulations (details will be given in Section 3). Mudde and Simonin (1999) have used transient two- and three-dimensional k–$ turbulence model to simulate the experimental data of Becker et al. (1994) with Euler–Euler approach. The governing equation for the two "uid models was derived by ensemble averaging of the local instant conservation equation of the single-phase "ows. The gas and liquid phase continuity equations were given by Eqs. (2.8) and (2.11), respectively. The momentum equation for the phase K was given by @uK + (jK K uK ) · ∇uK @t = − jK ∇p + jK K g − ∇ · (jK uK

uK

K ) + FK :

(jK K )

(2.17)

The contribution of molecular viscosity was neglected. The interfacial momentum exchange was considered to take place by drag force, virtual mass force and the correlation between instantaneous distribution of the particles and the undisturbed "uid pressure "uctuation. The lift force was not considered in the force formulation. Due to the averaging procedure, the drag force was not only considered to be proportional to the slip velocity but also involved a contribution of drifting velocity due to correlation between the distribution of particles and the turbulent "uid motion. The drifting velocity was represented for the dispersion of the particles due to transport by the "uid turbulence. However, the concept of drifting velocity approach su;ers from the

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

following limitations: (i) The drifting velocity arising out of dispersion mechanism cannot be considered as a part of relative velocity between the gas and liquid and authors may like to reconsider their Equation (9). (ii) The mechanism of drifting has not been considered in the continuity equation. (iii) In bubble columns, the hold-up gradients occur mainly in the radial direction. (iv) The authors have given the following equation for drag coeLcient: CD =

24 (1 + 0:15 Re0:687 )jL−1:7 Re

for Re ¡ 1000:

It may be pointed out that the >rst part of the right-hand side is applicable for solid particles and not for bubbles. The inclusion second part (namely jL−1:7 ) is a novel feature. However, the physical signi>cance of (−1:7) and its constancy over the wide Re range needs to be explained. The total force (FK ) was de>ned as   @ur FG = −FL = jG L FD ur + jG L CV + uG · ∇ur @t + ∇ · (jG L CV uG

uL

G ) − L uG

uL

G · ∇jG : (2.18) The k–$ model equations have been modi>ed to include the e;ect of the interface on the turbulence of the carrier phase. The complete discussion of k–$ model has been given by Thai Van, Minier, Simonin, Freydier, and Oliver (1994) and Bel F’Dhila and Simonin (1992). The turbulent predictions of the dispersed bubble phase were achieved by an extension of Tchen’s theory (Simonin, 1990). Three important time scales were considered for this purpose: the characteristic time of the energetic turbulent eddies (3=2C k=$), characteristic time of "uid turbulence viewed by the bubbles, and the characteristic time of the particle entrainment by the continuous "uid motion. P"eger et al. (1999) have used commercial CFD software to simulate the two-phase gas–liquid "ow in bubble columns. Unlike the earlier attempts of Lapin and Lubbert (1994), Sokolichin and Eigenberger (1994), Sokolichin et al. (1997) and Sokolichin and Eigenberger (1999), they have solved Reynolds averaged momentum balance equation for both the phases given as @ (jK K uK ) + ∇ · (jK K uK uK ) @t = − jK ∇p + jK ∇ ·  + jK K g + FD :

(2.19)

The coupling between phases has been modeled by the interfacial drag force. The drag coeLcient was arbitrarily selected as 0.66 for the uniform bubble diameter of 2 mm. The virtual mass force and the lift force were neglected. The continuity equations of both the phases were solved. The dispersion mechanism in the gas phase was considered in

the continuity equation: 
t @(jK K ) + ∇ · (jK K uK ) = jK ∇ · ∇jK : @t %&

3165

(2.20)

They have studied the e;ect of turbulent dispersion on the prediction of "ow patterns (details given in Section 3). Further, they have relaxed the laminar assumption and have considered turbulence in the two-phase "ows. The standard k–$ turbulence model was used for this purpose. However, they did not consider the contribution of the bubble-generated turbulence in the formulation. The gas phase viscosity has been de>ned following Grienberger and Hofmann (1992): G t; G = t:L : (2.21) L Following Sokolichin et al. (1997), an implicit higher order TVD scheme was used for the discretization of the convective "uxes. P"eger and Becker (2001) have used the hydrodynamic model of P"eger et al. (1999) (Eq. (2.19)). The drag coeLcient was now selected as 0.44 for uniform bubble diameter of 4 mm to de>ne constant slip velocity of 200 mm=s. The dispersion mechanism in the gas phase was not considered in the present case. Accordingly, the continuity equation given by Eqs. (2.8) and (2.11) were solved for the gas and the liquid phase, respectively. Turbulence was considered only for the continuous liquid-phase- and gas-phase-modeled laminar. However, its in"uence on the turbulence in the continuous phase was considered by a bubble-induced turbulence model. The production term in the turbulence model was modi>ed to include the bubble-induced turbulence as GL = L; e; (∇uL + ∇uLT ) : ∇uL + Ck=$ FD · |uG − uL |:

(2.22)

The constant Ck=$ correspond either to the k or the $ equation and were combined with standard k–$ constant. For the turbulent kinetic energy transport equation, Ck was assigned the value of C$1 (1.44) whereas for dissipation rate transport equation, C$ was taken to be same as of C$2 = 1:92. This approach is di;erent from Joshi (2001), where the bubble-induced turbulence in calculated (modeled) based on energy transfer coeLcient (CB ). Deen et al. (2001) have simulated the gas–liquid "ow in bubble columns using large eddy simulation (LES) and the standard k–$ model. The continuity equations for gas and liquid phases were given by Eqs. (2.8) and (2.11), respectively. The momentum balance equation were also written for both the phases as given by Eq. (2.19). The ensemble averaging technique was used in the case of k–$ turbulence model and >ltering technique was adopted for LES model. The instantaneous and "uctuating components refer to grid and sub-grid scale velocities after >lteration. The e;ective viscosity of the liquid phase was considered to be composed of three contributions; the molecular viscosity, the turbulent viscosity, and the extra term due to bubble induced turbulence. The model of Sato and

3166

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

Sekoguchi (1975) was used for the bubble-induced turbulence: BI; L = L C; BI jG dB |uG − uL |;

(2.23)

where a model constant C; BI was selected arbitrarily as 0.6. For the case of LES, the sub-grid scale viscosity model proposed by Smagorinsky (1963) was used to calculate the turbulent viscosity T; L = L (CS *)2 |S|;

(2.24)

where a model constant CS was taken to be 0.1. For both the models, the e;ective viscosity of the gas phase was calculated by using Eq. (2.21). The interaction between the phases was modeled by the interfacial forces, drag, lift and virtual mass. The lift and virtual mass force coeLcients were arbitrarily set to 0.5. Further, the lift force has been considered to be positive (i.e. bubbles move radially outwards). Following Ishii and Zuber (1979), the drag force coeLcient has been calculated as unity for distorted bubbles of 4 mm. 3. Flow driven by a bubble plume A single-point sparger gives a bubble plume and rest of the column is practically gas free. The "ow pattern under steady-state conditions have been analyzed by Harlow and Amsden (1975), Szekely, Wang, and Kiser (1976), Liles and Reed (1978), Deb Roy and Majumdar (1981), Carver (1982,1984), Sahai and Guthrie (1982a–c), Grevet, Szekely, and El-Kaddah (1982), Salcudean, Low, Hurda, and Guthrie (1983), Elghobashi and Abou-Arab (1983), Salcudean, Lai, and Guthrie (1985), Neti and Mohamed (1989), Lai and Salcudean (1987), Johansen and Boysan (1988), Lopez de Bertodano, Lee, Lahey, and Drew (1990), Lapin and Lubbert (1994), Sokolichin and Eigenberger (1994, 1999), Becker et al. (1994), Sokolichin et al. (1997), Delnoij et al. (1997a,b), Borchers et al. (1999), Mudde and Simonin (1999), P"eger et al. (1999), P"eger and Becker (2001) and Deen et al. (2001). The investigations up to 1990 were restricted to the simulation of steady-state "ows. Therefore, the following sections review the developments after 1990 and which were concerned with the transient "ows. 3.1. O1-center single-point sparging (two-dimensional column) 3.1.1. Becker et al. (1994) The pioneering contributions to the transient "ow simulations have been made by Lapin and Lubbert (1994) and Becker et al. (1994). The latter authors have documented the detailed observations from 0:5 m (width)×2:0 m (height)× 0:08 m (depth) two-dimensional column (Fig. 3). LDA was used for the measurement of the "ow pattern whereas photographic technique was used for the documentation of hold-up pattern. They have employed two gas "ow rates of 8 l=min (VG = 3:3 mm=s) and 1:6 l=min (VG = 0:66 mm=s).

Fig. 3. Schematic representation of a bubble column used by Becker et al. (1994).

For the gas distribution, frit sparger (8 mm diameter with 0:3 m pores) located at 0:1 m left to the column center was used. The authors’ observations have been reproduced below in somewhat greater detail because these observations have played a major role in the subject of transient "ow simulations. (i) Due to the uneven gas distribution, a gross circulation "ow over the whole height of the column develops, which pushes the bubble swarm on one side. At the high gas "ow rate of 8 l=min, it was observed that, as the bubbles approach the bed surface, their upward "ow is in"uenced by a secondary vortex near the left corner. Part of the bubbles sharply rise, change their path and move downward towards the column center. The liquid "ows upward in the bubble swarm and move downward along the sidewalls. Further, small vortices are formed in three corners of the two-dimensional column. (ii) At a lower gas "ow rate of 1:6 l=min, the "ow pattern is markedly di;erent from that obtained at the high "ow rate. Several liquid circulation cells are formed which continuously change their size and location. As a result, the bubble swarm partially follows these vortices in a meandering shape (Fig. 4). It may be

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

3167

Fig. 4. Photographs of oscillating bubble swarm in locally aerated "at bubble column (Becker et al., 1994; Sokolichin & Eigenberger, 1999).

Fig. 5. Transient measurement of the vertical component of liquid velocity at position A (x = 35 mm; y = 900 mm) and B (x = 450 mm; y = 1050 mm) (Sokolichin & Eigenberger, 1999). For coordinate system see Fig. 3.

emphasized that, the direction of the lower part of the bubble swarm is stable and directed against the near sidewall. However, the upper part changes its appearance and location corresponding to the transient liquid circulation "ows. Because of the transient "ow characteristics, the LDA measuring time of 30 –60 s is not suLcient to get time-averaged steady-state mean velocities. Long-time measurements at the point A (Fig. 5) show a periodical change of the vertical and the horizontal velocity component with a period time of 41 s. Further, a long-time average of liquid velocity at point A was found to be −38 mm=s. The above experimental observations of Becker et al. (1994) have been used for the transient CFD simulations by Becker et al. (1994), Sokolichin et al. (1997), Delnoij et al. (1997a–c), Sokolichin and Eigenberger (1999), Borchers et al. (1999), and Mudde and Simonin (1999). These investi-

gations have made substantial contributions in advancing the subject of transient "ow simulations. These developments have been presented below to bring out the present status in a coherent manner. Becker et al. (1994) used Euler–Euler approach with 18 × 25 grid size. Other details are given in the Section 2, Tables 1 and 2. When the liquid "ow was assumed laminar, the three small vortices (observed at high "ow rate of 8l=min) could not be obtained with the CFD simulation. For better production of vortices, the liquid viscosity had to be increased by a factor of 100 in the whole simulation domain. However, the turbulent liquid viscosity as calculated by the k–$ model resulted in a viscosity enhancement up to a factor of 20,000, which dampened the small vortices completely. At the low "ow rate of 1:6 l=min, the "ow simulation based on the laminar model and a coarse spatial grid (18 × 25) gave only stationary results which qualitatively agreed with the experimental results. Re>ning the spatial grid (36 × 50)

3168

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

Table 2 Summary of numerical simulation of the "ow pattern: previous work

Author(s)

Lapin and Lubbert (1994)

Approach

Solution procedure Model

Algorithm

Scheme

Grid (H × W × D)

Population balance 2D Laminar 2D Laminar 2D with 100 factor Laminar 2D

SIMPLER

—

200 × 60

SIMPLER

—

75 × 25

SIMPLER

—

SIMPLER

Laminar-2D Laminar 2D

— —

UPWIND TVD — —

25 × 18 50 × 36 150 × 50

Sokolichin and Eigenberger (1994)

Euler–Euler Euler–Lagrange Euler–Euler

Becker et al. (1994)

Euler–Euler

Sokolichin et al. (1997) Delnoij et al. (1997b) Delnoij et al. (1997a)

Euler–Euler Euler–Lagrangian Euler–Lagrange Euler–Lagrange

Delnoij et al. (1997c)

Euler–Lagrange

Laminar 3D

—

—

Delnoij (1999)

Euler–Lagrange

Laminar 3D

—

—

Sokolichin and Eigenberger (1999)

Euler–Euler

Laminar—2D, 3D Std. k–$ model

—

Higher order TVD

Borchers et al. (1999)

Euler–Euler

Std. k–$ model

SIMPLER

Higher order TVD

Mudde and Simonin (1999)

Euler–Euler

—

—

P"eger et al. (1999)

Euler–Euler

Std. k–$ model 2D, 3D Std. k–$ model 2D, 3D

—

Higher order TVD

P"eger and Becker (2001)

Euler–Euler

SIMPLEC

Deen et al. (2001)

Euler–Euler

Std. k–$ model 3D Std. k–$ model and LES 3D

Higher order TVD QUICK

could give transient solution with periodical movement of the bubble swarm. The simulated period time at point A (Fig. 5) was found to be 40 s which excellently agreed with the experimental value of 41 s. 3.1.2. Sokolichin et al. (1997) Sokolichin et al. (1997) used Euler–Euler and Euler– Lagrange approaches for the simulation of "ow pattern in the two-dimensional bubble column of Fig. 3. Other details of the model (Section 2) and the numerical simulation are given in Tables 1 and 2. The simulation results (50 × 150 grids) obtained with three di;erent numerical algorithms were compared with each other: (i) the Eulerian approach with upwind discretization, (ii) the Eulerian approach with TVD discretization and (iii) the Lagrangian approach for the gas equation. As a >rst case for comparison, the dispersion in bubble path was neglected (DG = 0 in Eq. (2.13)). The comparisons of the evolution of the vertical liquid velocity component in time at point A (Fig. 5) is shown in Fig. 6A. It can be seen that there is an excellent agreement between the implicit higher order TVD and the Lagrange solutions. However, the upwind solution shows a markedly

–

100 × 50 20 × 20; 40 × 20; 96 × 20 154 × 20; 228 × 20; 55 × 25 20 × 20 × 20; 40 × 20 × 20 96 × 20 × 20 20 × 20 × 20; 40 × 20 × 20 96 × 20 × 20; 154 × 20 × 20 50 × 36; 150 × 50 150 × 50 × 8; 75 × 25 × 4 225 × 75 × 12 50 × 36; 150 × 50 150 × 50 × 8; 75 × 25 × 4 225 × 75 × 12 52 × 27; 52 × 38 52 × 27 × 10; 52 × 38 × 18 90 × 44 × 15; 45 × 22 × 5 60 × 33 × 7; 90 × 44 × 7 90 × 44 × 1 80 × 10 × 17; 60 × 8 × 14 60 × 5 × 9; 30 × 5 × 9 —

di;erent behavior, which was attributed to the tremendous in"uence of the numerical di;usion in the Eulerian solution obtained with the upwind discretization. It may be emphasized however that, though the TVD and Lagrange solutions agreed with each other (Fig. 6A), they did not compare favorably with the experimental pro>le. Further, the spread of the bubble plume (Fig. 7, lower row) calculated with the implicit higher order TVD or Lagrange method was found to be much smaller than those observed experimentally. These di;erences were attributed to the assumption of zero gas phase dispersion due to random motion. Therefore, the authors included the term DG in Eq. (2.13). It can be seen from Fig. 6B that the implicit higher order TVD and the Lagrange solutions again show very similar behavior; however, do not agree with the experimental observations. The prediction of bubble plume by Lagrange approach gave a qualitative comparison with the experimentally observed shape and size of the bubble plume. 3.1.3. Delnoij et al. (1997b) Delnoij et al. (1997b) used Euler–Lagrange model which describes the time-dependent two-dimensional motion of

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

Fig. 6. Vertical liquid velocity at position A calculated with di;erent methods (Sokolichin et al., 1997): (A) Without dispersion of gas phase and (B) with dispersion of gas phase.

small, spherical gas bubbles in a bubble column. The importance of the added mass force, the lift force, the drag force and the hydrodynamic interaction force acting on the bubbles was investigated theoretically. It was argued that the added mass force is of particular importance near the gas distributor and should be accounted for to realistically predict the behavior of bubbles in the gas distributor region. The lift force acting on the bubbles disperses the bubbles over the cross section of the column. It was shown that a model that neglects the lift force acting on a bubble could not predict the experimentally determined wall peaking of the void fraction. It was also argued that the particular approximation (Eq. (2.15)) for the drag coeLcient used in the computations does not signi>cantly a;ect the "ow pattern ultimately predicted by the model. Finally, the model was used to study the e;ects of the force due to hydrodynamic interaction between bubbles on the overall "ow pattern. It was shown that the hydrodynamic interaction force exerted by the bubbles on each other does not signi>cantly a;ect the macroscopic "ow pattern observed in a bubble column operated in the homogeneous "ow regime at relatively low void fractions (¡ 5%). Delnoij et al. (1997b) applied their model to predict the experimental behavior reported by Becker et al. (1994) at high (8 l=min) as well as low (1:6 l=min) "ow rates. In the former case, the model has been shown to compare reasonably well the behavior of bubble plume. The only exception was that Delnoij et al. (1997a–c) predicted unsteady move-

3169

Fig. 7. Comparison of the simulation results at >ve di;erent times of (A) laminar Euler–Lagrange model by Delnoij et al. (1997a–c) with (B) laminar Euler–Euler model by Sokolichin and Eigenberger (1999). (Grid size for both the cases: 50 × 100 and Ut = 10 s).

ment of secondary vortex, which was found to move up and down in the upper part of left column wall. For the case of low "ow rate also, the behavior of bubble plume was qualitatively predicted. The form of the undulation was found to di;er from that observed in the experiment. In particular, the stable lower part of the bubble swarm, which is always directed against the near sidewall in the experiment was not correctly reproduced in the simulations. Quantitatively, large di;erences were observed. At point A of Fig. 5, the predicted average axial liquid velocity was found to be −77 mm=s as against the observed velocity of −38 mm=s. Further, the period of oscillation was predicted to be 30 s which is much lower than the experimental value of 41 s. 3.1.4. Mudde and Simonin (1999) Mudde and Simonin (1999) have used transient twoand three-dimensional k–$ turbulence model, to simulate the experimental data of Becker et al. (1994) with Euler– Euler approach. The oscillatory nature of "ow was not observed with two-dimensional k–$ model. The steady state was obtained in about 20 s with one liquid circulation cell by two-dimensional simulations. The meandering movement of the bubble plume was dampened completely due to high turbulent viscosity calculated by k–$ model. In the case of three-dimensional simulations the "ow becomes transient, however, the simulation heavily underpredicted the

3170

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

oscillation period (4.44 against 41 s), the amplitude (63 against 200 mm=s) and also the mean velocity (334 against −38 mm=s). The grid size variation could not improve the results. The meandering movement of the bubble plume was not observed in the simulated results. The plume remains attached to the left wall. The oscillations were caused by the "uctuations in the width of the bubble plume at height of about 400 mm. While moving upwards, the thickness of the bubble plume was observed to increase and took the shape of a ‘sawtooth’ accompanied by a small liquid circulation cell with an ‘eye’ at the plume boundary. In three-dimensional simulations, the di;usion of the bubble plume was observed to be much lesser than in two-dimensional case. The authors attribute this fact to the lower turbulent viscosity obtained. It may be noted that Delnoij et al. (1997a) obtained the spread in the plume by considering the lift force in their model formulation. With the inclusion of virtual mass force, the oscillation period increased to 34 s, the amplitude of the vertical liquid velocity at the monitor point reduced to 200 mm=s and obtained slightly negative mean liquid velocity. The oscillatory motion of the plume is induced by the vortex in the upper left corner. The vortex moves downwards along the wall and pushes the plume to the middle. The vortex gets destroyed as it reaches the lower part of the column by the large liquid circulation at the lower right part of the vessel. Delnoij et al. (1997b) have demonstrated a di;erent sensitivity of the virtual mass force. They observed the in"uence of the virtual mass force in the near sparger region. It also contributed to the stability of the solution. Further, with the low Reynolds k–$ model, signi>cant changes were not observed as compared to the k–$ model simulation. In the forgoing discussion, the experimental results of Becker et al. (1994) have been simulated by Becker et al. (1994), Sokolichin et al. (1997), Delnoij et al. (1997b) and Mudde and Simonin (1999). The numerical solution was found to depend strongly on grid size, time step, discretization method, liquid viscosity, turbulence model and the description of bubble phase random motion. The solution was also found to depend upon whether Euler–Euler or Euler– Lagrange approach was employed. In none of the cases, real satisfactory predictions were made to simulate all the observations summarized in Figs. 4 and 5. For these observed discrepancies, the underlying reasons have been elegantly analyzed recently by Sokolichin and Eigenberger (1999). For this purpose, they have used laminar, two- and three-dimensional k–$ model. Further, they have investigated the e;ects of grid size, time step and the discretization method. 3.1.5. Sokolichin and Eigenberger (1999) 3.1.5.1. Simulation results for laminar model. Simulation results for the test case described above were >rst presented by Becker et al. (1994). Their results were based upon a mathematical model, which is very similar to the laminar model of Euler–Euler type, described in Section 2.

The convective terms in all equations were discretized with a >rst-order upwind scheme. Numerical results for two di;erent spatial grids were presented. The "ow simulation based on a coarse spatial grid (18 × 25 points) gave a stationary result. The calculated liquid velocity >eld agreed qualitatively, however, with the long-term-averaged results of LDA measurements. Re>ning the spatial grid to 36 × 50 points lead to a transient solution. The simulation showed the observed periodic movement of the bubble swarm. The calculated period time, the long-time velocity "uctuation at point A (see Fig. 5) and even the long-time averaged velocity >eld were in very good agreement with experiments. One question, however, still remained open: how good did the presented simulation results correspond to the underlying mathematical model, and what will happen if the spatial resolution would be further re>ned? Delnoij et al. (1997b) tried to simulate the two-phase "ow in a "at bubble column in the frame of the Euler–Lagrange model. They also neglected turbulence e;ects and obtained on a 50×100 space grid the quasi-periodic dynamic solution with an undulating bubble swarm (Fig. 7, upper row). The form of the undulation, however, was found to di;er from that observed in the experiment. In particular, the stable lower part of the bubble swarm, which is always directed against the near sidewall in the experiment (Fig. 4), was not correctly reproduced in the simulation. Also, the period time was undersimulated in the calculation by a factor of 1.9. In Delnoij et al. (1997b), the motion of gas bubbles was calculated from a force balance for individual bubble, accounting for ‘all relevant forces acting on them’. In the lower row of Fig. 7, simulation results obtained with laminar model with constant slip velocity are presented. The resolution in space and in time is the same as in Delnoij et al. (1997a) and the TVD-scheme (which is second-order fully implicit) is used for the convective terms of the gas continuity equation in order to reduce the numerical di;usion. Despite the strongly simpli>ed version of gas momentum balance of Sokolichin and Eigenberger (1999) (described in Section 2) the results obtained with both models are in good agreement (Fig. 7). Both the solutions can be seen to be stationary and periodic. The period time and the amplitude of the "uctuation of the vertical liquid velocity at point A from Fig. 5 were also found to be in a very good quantitative agreement. This comparison shows that, the main reason for the di;erences in the form of the undulation obtained by Becker et al. (1994) and Delnoij et al. (1997c) is the >ner spatial grid used in the latter publication, and not the di;erent treatment of the gas momentum balance. Sokolichin and Eigenberger (1999) also investigated the e;ect of grid size for the laminar model by taking >ve different grids with 25 × 75, 50 × 150, 100 × 300, 200 × 600 and 400 × 1200 grid points. The grid re>nement did not lead to a converging solution. Instead, the number of circulation cells resolved by the simulation was found to grow continuously with the spatial resolution (Figure 3 of Sokolichin and Eigenberger). This is a typical result of a laminar calculation

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

of a turbulent "ow. This observation is truly a signi>cant result. The following conclusions can be drawn from the numerical experiments with the laminar model. The laminar model neglects the e;ects of turbulence. The >ner the spatial grid the more are the vortices resolved. Therefore, the result of the simulation depends strongly on the number of grid points used and cannot be considered reliable. A very good agreement between the simulations and experiments as reported by Becker et al. (1994) was pointed out (by Sokolichin and Eigenberger) as a pure coincidence. It is due to the fact that, the numerical di;usion of the upwind discretization has a similar in"uence as the turbulent eddy viscosity in the turbulence models as will be shown in the following. If the space grid is re>ned or an implicit higher-order TVD scheme is used, the numerical di;usion decreases, and the solution of the laminar model becomes more and more unstable. 3.1.5.2. Simulation results for the turbulent model. All the results presented so far were obtained with two-dimensional models. The prime reason for using two-dimensional models is the assumption that, in a "at bubble column with a small depth, the "ow structure has an essentially two-dimensional character. The results of LDA measurements con>rm, that the long-time-averaged velocities do not change remarkably in depth. The same is, however, not true for the high-frequency component of the velocity "uctuations and consequently for the turbulent kinetic energy. The results of three-dimensional simulations with the turbulent Euler–Euler model show that the front and the back walls indeed dampen the intensity of turbulence inside the bubble column, so that the turbulent eddy viscosity becomes about one order of magnitude smaller than that in the two-dimensional simulation. The overestimation of the effective viscosity in two-dimensional simulation (in the middle of the column values up to 5 kg=ms are reached; this value is 5000 times higher than the laminar viscosity of water) is the main reason for the non-existence of the dynamic solution in the two-dimensional case. The overprediction of the width of the bubble swarm by Sokolichin and Eigenberger (1999) and also by Sommerfeld, Deckwer, and Kohnen (1997) is also due to the higher values of the turbulent kinetic energy, because in both the models, the intensity of the bubble dispersion is related to the intensity of turbulence in the liquid phase. The results of the time-dependent three-dimensional simulation showed good agreement with experiments (Fig. 4). The long-time vertical liquid velocity "uctuations at positions A and B (see Fig. 5), calculated with the three-dimensional turbulent model are shown in Fig. 8. The high-frequency part was found to be >ltered out by the turbulence model. Both curves from Fig. 8 can be described as quasi-periodic. Their amplitude and frequency can be seen to be in good agreement with the low-frequency component of the measured velocity (shown in Fig. 5).

3171

Fig. 8. Comparison between the experimental data of Becker et al. (1994) (thick line) and predicted vertical liquid velocity calculated with three-dimensional k–$ model of Sokolichin and Eigenberger (1999) (thin line).

Also the comparison between the long-time-averaged velocities obtained from the LDA measurements and from three-dimensional simulations, especially in the bubble-free region, is quite convincing (Figure 10 in Sokolichin and Eigenberger, 1999). In order to estimate the error of the space discretization the simulations were carried out with a space resolution of 75 × 25 × 4, 150 × 50 × 8 and 225 × 75 × 12 grid points. Because of the chaotic character of the solution, it was not possible to make a direct quantitative comparison between instantaneous "ow pattern results obtained on di;erent space grids. However, the authors have presented an indirect comparison through the calculation of the long-time average velocity patterns. In Figure 10 of Sokolichin and Eigenberger (1999), the vertical liquid velocity pro>les, averaged over 4 and 20 min of real time, are shown for three grid resolutions. One can see that after 4 min of averaging, the results for all the three grids di;er remarkably, whereas for 20 min of averaging, the pro>les obtained on the 150 × 50 × 8 and 225 × 75 × 12 space grids show a very good agreement. The dependence of the solution from the time step was also studied at three levels: 0.05, 0.1 and 0:2 s. The long-time averaged velocity and gas hold-up pro>les for all three-time steps were found to be practically undistinguishable. Sokolichin and Eigenberger (1999) have also investigated the real contribution of high-order discretization and have recommended a strategy of combination of high- and low-order schemes. 3.2. Single-point central sparging (two-dimensional column) In order to check the validity of the conclusions of previous section of the o;-center sparging to the case of central sparging, EFD and CFD results have been reported by Borchers et al. (1999), P"eger et al. (1999) and Becker

3172

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

Fig. 9. E;ect of aspect ratio on the number of circulation cells and shape of bubble plume (Borchers et al., 1999): (A) H=W = 1; (B) H=W = 1:5; (C) H=D = 2 and (D) H=W = 3.

et al. (1999). The numerical results of Delnoij et al. (1997a) have not been included because laminar "ow has been considered in all the cases. 3.2.1. Borchers et al. (1999) The experimental studies have been performed in a two-dimensional column of 0:5 m width, 2:0 m height and 0:08 m depth (0:5 m × 2:0 m × 0:08 m), the same column as that used by Becker et al. (1994) (Fig. 3). However, Becker et al. (1994) used o;-center sparging, whereas, Borchers et al. (1999) have used central sparging (a frit sparger of 8 mm diameter with 0:3 m holes). Further, Borchers et al. (1999) have used the height to width (H=W ) ratios of 1, 2 and 3. In the case of H=W =1, the "ow was found to have an essentially stationary character in which only one circulation cell prevailed over the full width of the bubble column (Fig. 9A). As the H=W ratio was increased slightly over 1.5, an unsteady structure was found to develop with two staggered rows (Fig. 9B) of vortices moving downwards in a periodic way. The number of vortices increased with the liquid height. Three vortices were observed for H=W = 2 (Fig. 9C) and four vortices in the case of H=W = 3 (Fig. 9D). The measured time series at point A (Fig. 9A) and point C (Fig. 9C) are shown in Fig. 10. The time series at point A (shown in Fig. 10A) illustrates that the "ow in the column with H=W = 1 has statistically stationary behavior with a constant mean velocity. High-frequency turbulent "uctuations around the mean velocity are of the order of the mean velocity itself. The time series at point C (Fig. 10B) corresponds to the quasi-periodic "ow in the bubble column with H=W = 2. In addition to the high-frequency "uctuations, low-frequency oscillations of the liquid velocity can be seen, which correspond to the periodic movement

Fig. 10. Time series of vertical liquid velocity at points A and C of Fig. 9 (Borchers et al., 1999).

of the bubble swarm. The radial pro>les of axial velocity (long-time average) were measured at 4, 8 and 13 axial locations for H=W = 1, 2 and 3, respectively. The CFD model developed by Sokolichin and Eigenberger (1999) was found to be useful for the prediction of average mean axial velocity and periodic "ows. For instance, a comparison of measurements and simulation at points B (x=10 mm; z=500 mm) and C (x=250 mm; z=500 mm) of Fig. 9C was carried out for a gas "ow rate of 1 l=min and aspect ratio of 2. The time series of the vertical liquid velocity at these points over a measurement time of 200 s are shown in the top row of Fig. 11 (A and B). A Fast Fourier Transfer was applied to the data in order to divide the low-frequency motion caused by the bubble swarms periodic oscillation from the high-frequency "uctuations caused by turbulence. In the column center (point C), the period of the vertical velocity oscillation (17 s, Fig. 11B) can be seen to be half as long as that at point B (34 s, Fig. 11A). The reason is that the liquid vortices of the right and left sides both pass the central position, while closer to the column side only one of the vortices passes by. The phase trajectories (a plot of vertical velocity component versus horizontal velocity component) at points B and C have been shown in Figs. 11C and D, respectively. The FFT->ltered data was used in order to visualize the low-frequency, quasi-periodic development of the liquid velocity. At point B, the velocity oscillation of both the components is nearly sinusoidal. Further, horizontal and vertical ◦ velocity is shifted in phase by about 90 . Therefore, the corresponding phase trajectory has the shape of an ellipse. In the column center (point C), the vortices pass on both the sides so that the phase trajectory resembles at points B and C result from the same vortex movement, their mean cycle times are identical (34 s). For comparison with experimental data, the simulated time-dependent velocity components

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

3173

Fig. 11. Analysis of velocity time data at points B and C of Fig. 9(C). (A, B) Time series of vertical liquid velocity. (C, D) Phase-plane plot of measured vertical liquid and horizontal liquid velocities. (E, F) Corresponding simulation results. (Borchers et al., 1999).

were low-pass >ltered as the measurement data. The resulting phase trajectories at points B and C are shown in Fig. 11E and F for a time of 200 s. It can be seen that the shape and velocity values are in reasonable agreement with the experiment. Though Figs. 11C–F do not show quantitative comparison of the periodic component, the authors have shown an excellent comparison for the average axial velocity at all the 4, 8 and 13 locations for H=W = 1, 2 and 3, respectively. From the foregoing description of the investigations of Borchers et al. (1999), the following points emerged:

sented for the high-frequency "uctuating velocity. Perhaps, the possibility existed for the comparison of k, $ and t pro>les or the complete energy balance. It may be emphasized that the motivation of transient "ow simulation has been the better representation of "ow pattern so that the actual design (for instance, from mixing point of view) is reliable. This objective is possible only when the success of prediction of t pro>les is achieved.

(i) for the transient "ow simulation the recommendations of Sokolichin and Eigenberger (1999) hold. For the appropriate estimation of viscosity, three-dimensional: k–$ model needs to be used with proper grid size particularly in the depth direction. (ii) Borchers et al. (1999) as well as Sokolichin and Eigenberger (1999) have shown an excellent agreement between the predicted and the experimental values of axial velocity at practically all the locations. However, for the low-pass >ltered periodic component the comparison is only reasonable. (iii) with the well-known limitations of the present day computation facility, no comparison could be pre-

3.2.2. P7eger et al. (1999) P"eger et al. (1999) used 0:2 m (width) ×0:45 m (height) ×0:05 m (depth) two-dimensional column with central sparging. The system was air–water and the gas "ow rate was varied in the range of 0.33–1:5 l=min. The measurements of "ow pattern (velocity–time series, radial pro>les of axial velocity and bubble plume behavior) were made by LDA and PIV. The CFD simulation was carried out using commercial software CFX (version 4.2). The authors have made observations similar to Sokolichin and Eigenberger regarding the need of three-dimensional: k–$ simulation with adequate number of grids. It is very interesting and satisfying to note that two altogether di;erent investigating

3174

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

groups have arrived at the same conclusions. This work also su;ers from the limitation described under (iii) of Section 3.2.1. Further, P"eger et al. (1999) have made an additional observation that; the inclusion of gas phase dispersion in the CFD formulation is undesirable. This is perhaps due to the need of proper quanti>cation of gas phase dispersion. 3.2.3. Becker et al. (1999) Becker et al. (1999) used 0:2 m (width)×0:45 m (height) × 0:04 m (depth) two-dimensional column with central sparging of 1 mm diameter hole. The system was air–water with a volumetric gas "ow rate of 0:8 l=min(VG =1:7 mm=s). Becker et al. (1999) have not reported CFD simulation and the experimental observations have been similar to those of Borchers et al. (1999). 4. Uniform sparging in two-dimensional columns 4.1. Experimental observations of Chen, Jamialahmadi, and Li (1989) Chen et al. (1989) studied the "ow pattern in two-dimensional columns. They used two columns (0:76 m width × 1:5 m height × 0:05 m gap and 0:175 m × 1:8 m × 0:015 m gap). The >rst column had a single pipe gas distributor made up of 45 mm diameter tube with upward facing evenly spaced 99 holes of 1 mm diameter. In the second column, wire gauze acted as a gas distributor. The "ow visualization was carried out using photography. When the ratio of height to width HD =W was ¡ 1, a single circulation cell (gulf stream "ow pattern) was observed similar to Freedman and Davidson (1969) and Whalley and Davidson (1974). When the HD =W was ¿ 1, multiple cells were found to appear with a staggered orientation. 4.2. Lapin and Lubbert (1994) Lapin and Lubbert (1994) have used Euler–Euler approach for the simulation of two-dimensional column with 1 m width and 1:5 m height (1 m × 1:5 m). Symmetry was assumed and hence half of the column (0:5 m × 1:5 m) was simulated with a grid structure of 60 × 200. For the liquid phase "ow, laminar conditions were assumed. The gas phase momentum balance was not solved and a constant value of slip velocity (50 mm=s) was assumed. The gas holdup at the bottom of the column was taken to be 3%. The gas phase dispersion was not considered. Other details are given in Section 2, Tables 1 and 2. From the simulation results, authors concluded that the two-"uid Eulerian simulations are very sensitive to false numerical di;usion. This observation is due to the assumption of laminar conditions and the non-inclusion of random motion of bubbles and the resulting dispersion. Since the Euler–Euler approach was found to su;er from the problem of numerical di;usion, Lapin and Lubbert

alternatively tried Euler–Lagrange approach. The uniform gas distribution was selected in such a way that the local gas holdup along the bottom was kept constant at 4%. The slip velocity was assumed to be 50 mm=s and to be constant throughout the column. When simulation was started at t =0 (s), the liquid was stagnant and free of bubbles. Initially, the bubbles rise in the liquid in the shape of a plane front, but after about 2 s, the front becomes unstable. After 5 s, the front becomes slightly wavy. However, a considerable instability develops near the wall leading to a local circulation in "ow due to noticeable pressure inhomogeneities. This enhances the instability and leads to >ngering in the bubble front (Figure 7(a) of Lapin & Lubbert, 1994). The density inhomogeneities shown in the >gure result in buoyancy e;ects and, hence, convective "ows, which shortly lead to a chaotic appearance of the "ow pattern and the bubble position in the "ow. The situation after 20 s when 16,952 bubble clusters are in the system can be characterized by randomly appearing local circulation cells, and the bubble positions show an irregular structure leading to considerable density di;erences. Interestingly, this highly irregular motion becomes more orderly with time. Furthermore, the inhomogeneities in the gas holdup, as manifested by bubble location plot, become insigni>cantly smaller. At the same time, the density also becomes more homogeneous. About 70 s after the start of numerical experiment (26,124 bubble clusters), velocity pattern approaches fairly regular and three circulation cells appear >lling the entire column. Then 26,124 bubble clusters are more or less homogeneously distributed over the whole wafer. The "ow structures seen are not stable. The centers of the circulation cells move and the size of the cell changes. Eventually, one of the cells becomes so small that it cannot even be recognized as a separate cell. Then, the outer appearance of the "ow pro>le is one containing only two major circulation cells, which depicts the "ow situation 110 s after startup. The appearance of the 28,067 bubble position plot and the local density plot does not change signi>cantly. The "ow remains instationary, sometimes showing two, sometimes three circulation cells in the velocity pattern. Interestingly, the two cells circulate in the same sense. This is similar to the idea behind Joshi and Sharma’s (1979) assumption. Such "ow models have often been criticized in literature with the argument that this would lead to complications at the interface. As shown in Fig. 7(d) of Lapin and Lubbert (1994), the local "ow is easily providing a sliding interface, however, the numbers of circulation cells as well as their positions considerably change with time. Only few operational conditions were found in which the pattern is fairly stationary. As most of the experimental information currently available in the literature is a long-time averaged information, it is of interest to simulate "ow patterns averaged over the period which correspond to integration time constants of the local anemometer measurements. The averaged velocity

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

pattern depicts much less structures and the secondary pattern becomes smoothed out with time leading to the simple "ow structure measured in nearly all the experiments performed in bubble columns (e.g. Yang, Rustemeyer, Buchholz, & Onken, 1986). 4.3. Sokolichin and Eigenberger (1994) Sokolichin and Eigenberger (1994) have used Euler– Euler approach for the two-dimensional simulation of the bubble column of 0:15 m width and 0:75 m height (0:15 m × 0:75 m). The column with uniform aeration was simulated with a grid structure of 25 × 75. The laminar conditions and a constant value of slip velocity (200 mm=s) were assumed. The gas phase dispersion was not considered. The authors have observed that, above a minimum value of the super>cial gas velocity (about 20 mm=s), an unsteady "ow structure develops where a sequence of vortices is created near the gas distributor. The vortices were found to have a mean diameter corresponding to the column diameter. Contrary to the approximation of Millies and Mewes (1994), the vortices are not stationary but move upwards at alternating sides of the column walls. It was observed that the mean number of vortices present at any time corresponds to the aspect ratio of the bubble column. If the calculated local velocities and the gas holdups were averaged over longer time period, a regular "ow structure with one overall circulation cell was obtained. Further, the authors concluded that only dynamic simulation leads to a converged solution and no steady-state solution could be obtained above a certain gas velocity if a suLciently >ne space grid was used. 4.4. Sokolichin et al. (1997) Sokolichin et al. (1997) have used Euler–Euler and Euler– Lagrange approaches for the two-dimensional simulations of the "ow pattern in the bubble column of Fig. 3, however, provided with uniform sparging. The simulation results over the entire bottom obtained by the Lagrangian approach with di;usion have been presented. It was observed that, above a minimum value of the super>cial gas velocity (about 20 mm=s), an unsteady "ow structure was found to develop. When the calculated local velocities were averaged over a longer time period, a regular "ow structure with one overall circulation cell resulted. The long-time averaged velocity patterns obtained with three di;erent numerical algorithms (discussed in Section 3) were also compared with each other at three di;erent heights. A very good quantitative agreement between the implicit higher order TVD and the LAGRANGE solutions was seen and the UPWIND scheme was found to lead to quantitatively di;erent results. The authors also have observed that the in"uence of the numerical di;usion in the case of a uniformly aerated bubble column is not so high as

3175

that in the case of a locally aerated bubble column. This may be because of the smoother distribution of the gas phase. 4.5. Deen et al. (2001) Deen et al. (2001) have used bubble column of a cross-section (W × D) of 0:15 m × 0:15 m and a height of 1m with uniform sparging. The non-coalescing system of air–aqueous salt solution was studied. The CFD simulation was carried out for super>cial gas velocity of 4:9 mm=s using the commercial software, CFX (version 4.3). They have compared the performance of the two turbulence models, k–$ and LES, with and without inclusion of the bubble-induced turbulence. The e;ect of di;erent interface forces on the simulated results was also investigated. The e;ect of virtual mass force on the simulated results was found to be negligible. With only drag force, the bubble plume did not show any transverse spreading. It may be recalled that Delnoij et al. (1997b) have reported a strong in"uence of virtual mass force on the simulations (Section 3.1.3). Out of the two turbulence models, they have found that LES simulations captured the strong transient movement of the bubble plume as observed in the experiment. The LES model resolved much more details of the "ow and large vortices were observed along the bubble plume. Due to the high turbulent viscosity, a quasi-stationary state was obtained for the k–$ model and transient details were not resolved but implicitly contained in the turbulent kinetic energy. Further, they have observed that the consideration of bubble-induced turbulence has a marginal e;ect on the predictions. It may be noted that P"eger and Becker (2001) have observed a strong in"uence of bubble-induced turbulence on simulation results. This may be because of di;erence in formulation. 5. Cylindrical columns 5.1. Experimental observations 5.1.1. Devanathan, Moslemian, and Dudukovic (1990) Devanathan et al. (1990) measured the "ow patterns by a computer-automated radioactive particle tracking (CARPT). The technique consists of monitoring the motion of a single neutrally buoyant radioactive particle by an online computer to map the "ow >eld. Mean and turbulent "ow velocities were measured in a 0:212 m i.d and 0:584 m long column. A porous stainless steel plate (average pore size 40 m) was used as a sparger. Visually, the column appeared to be operating in the heterogenous regime with a wide bubble size distribution. The "ow pattern was found to consist of two circulation cells. The primary "ow pattern consisted of a single circulation cell. A small secondary cell was observed near the distributor. In addition, several secondary axi-symmetric "ows were detected. The one-dimensional velocity pro>le was found to agree with that reported by Hills (1974). The

3176

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

Reynolds stress pro>le was found to be similar to that of single-phase pipe "ow, however, the magnitude of the stress was high. 5.1.2. Drahos, Zahradnik, Fialova, and Bradka (1992) Drahos et al. (1992) measured the liquid "ow patterns in 0.14 and 0:292 m i.d columns. Flow structures were identi>ed by analysis of wall-pressure "uctuations and by "ow visualizations. They observed a single-cell circulation pattern when HD =D ¡ 2. At higher aspect ratios, multicell "ow pattern was observed. A vertical scale of a single circulation loop equaled approximately to two-column diameters in the entire range of turbulent bubbling conditions. They also observed the existence of back"ow between the adjacent circulation cells. 5.1.3. Tzeng, Chen, and Fan (1993) Tzeng et al. (1993) used a two-dimensional column (0:483 m×1:6 m×0:0127 m) to investigate the macroscopic structures and the mechanisms of the liquid circulation. The gas distributor was composed of multiple injectors, which were individually regulated to generate desired gas "ow rate, bubble injection frequencies and bubble sizes. Colored bed particles and neutrally buoyant particles as solid and liquid tracers, respectively, were used for "ow visualization through video photography. The bubble streams injected near both the sidewalls were observed to migrate towards the bed vertical axis and vortices appeared along the sidewalls when the super>cial gas velocity exceeded 4 –6 mm=s. Based on bubble dynamics and the local liquid "ow patterns, they found four distinct "ow regions when a gross "ow circulation occurred in the system. As shown in Fig. 12, the four regions were central plume region, fast bubble "ow region, vertical "ow region and descending "ow region. The hold-up distribution of each phase was highly non-uniform. Interactions between the fast bubble "ow region and the vertical "ow region were observed to govern the instantaneous "ow >eld. 5.1.4. Chen, Reese, and Fan (1994) Chen et al. (1994) investigated the macroscopic "ow structures in a cylindrical bubble column (0:102 m ×2:2 m). A tube ori>ce type of distributor was used to provide uniformity of gas and liquid "ows. Particle image velocimetry (PIV) was used for the quanti>cation of "ow structures. Flow visualization was also conducted with the help of laser sheeting technique. Under low gas "ow velocities, dispersed bubble regime prevailed. Bubbles were found to rise rectilinearly and the liquid was found to fall between the bubble streams. Further increase in the gas velocity resulted in the transition of "ow regime into vortical spiral "ow regime. In this regime, clusters of bubbles or coalesced bubbles form the central bubble stream moving in a spiral manner with the liquid phase not only moving in a vortical pattern but also spirally downward in the region between central bub-

Fig. 12. Four "ow regions observed by Tzeng et al. (1993) in two-dimensional bubble column.

Fig. 13. Flow pattern observed by Chen et al. (1994) in a cylindrical bubble column.

ble stream and the column wall. Four "ow regions (similar to those in a two-dimensional column)—descending region, vortical spiral "ow region, fast bubble "ow region and central plume region were identi>ed (Fig. 13). The four "ow regions were instantaneous phenomena, which deviated from the "ow characteristics obtained when time averaging procedures were used to quantify the "ow properties. The bubble coalescence became dominant and formed large bubbles

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

as the "ow regime moved into the turbulent "ow regime at high gas velocities. Local chaotic motion of the liquid phase (caused by bubble wake and drift e;ects due to bubble motions) was found to destroy the vortical and spiral "ow structures and resulted into turbulent "ow structures. 5.1.5. Reese and Fan (1994) Reese and Fan (1994) investigated the hydrodynamic behavior in the entrance region of a cylindrical bubble column (0:102 m × 2:2 m). Two di;erent gas distributors (porous plate and perforated plate) were studied to elucidate their e;ects on the overall hydrodynamic behavior of the column and the entrance region. The instantaneous and the average "ow information were obtained using PIV system. In the dispersed bubble regime, the liquid phase was found to be more turbulent than the "ow occurring in the bulk region. The averaged liquid velocity and the "ow in the entrance and bulk regions consisted of a single circulation cell in the axial direction, with the liquid ascending in the center of the column and descending along the column wall. The strength of this re-circulation decreased as the "ow developed throughout the entrance region and into the bulk region. In the coalesced bubble region (including vortical– spiral and turbulent "ow), the "ow development was found to be quicker and the entrance e;ects were found to diminish. The overall and the entrance regions were found to be independent of the distributor design. 5.1.6. Lin, Reese, Hong, and Fan (1996) Lin et al. (1996) investigated the e;ect of the scale of a two-dimensional bubble column on the macroscopic hydrodynamic characteristics. When the column width was greater than 200 mm, the transition from the dispersed bubble "ow regime to the 4 and then to 3 region "ow in the coalesced bubble regimes was directly related to measurable coherent "ow structures. Thus, the dispersed bubble regime and the coalesced bubble regime were demarcated based on the behavior of the vortex size while the variation of the wavelength could di;erentiate between the 4 region and the 3 region "ow condition which are the characteristic features of the coalesced bubble regime.

3177

quality of the structures in two-dimensional and the cylindrical columns have been described in Section 4.2. 5.2.2. Millies and Mewes (1994) Millies and Mewes (1994) have presented a detailed derivation of circulating cell model. They showed that a stability analysis of the simpli>ed equations of motion exhibit instabilities that can be attributed to the recycle of liquid and thus mark the height of a circulation cell and was found to be equal to one column diameter. 5.2.3. Grevskott et al. (1996) Grevskott et al. (1996) observed the existence of multiple cell structure even under steady-state simulation conditions. These authors included improved formulation for the radial inward movement of bubbles on the basis of the investigation of Jakobsen, Svendsen, and Hjarbo (1993). In the three-phase "ow experiments, the existence of two primary solid circulation cells was established. Analysis of the Reynolds stress pattern indicated secondary cell structures extending from the bottom of the upper primary circulation cell up to the top of the dispersion layer. The cell size in this region was of the order of column diameter and evenly distributed along the vertical direction. The authors have opined that, these ‘cells’ might have been the result of recurring transient cells in these speci>c positions. They have expressed a need for future studies. 5.2.4. Becker et al. (1999) Becker et al. (1999) carried out LDA measurements in 0:295 m i.d, 2:5 m height bubble column. The time-dependent "ow behavior in a cylindrical column was found to be chaotic and not predictable (as compared to near periodic "ow in two-dimensional columns as shown in Fig. 14). Some regular structures were detected in the velocity–time series in axial and tangential directions that exhibited some kind of periodicity. However, no low-frequency "ow structure was observed that can be considered as characteristic for the whole column. The velocity–time series typically contain a broad range of frequency contributions. In view of this, the velocity–time series measured in a three-dimensional column (Fig. 15) does

5.2. Numerical simulations 5.2.1. Lapin and Lubbert (1994) In a cylindrical geometry, Lapin and Lubbert (1994) have shown the situation to be di;erent from the two-dimensional bubble column. The average velocity pattern shows a single circulation cell only, covering the entire column, while the instantaneous "ow >eld depicts some sub-motions, which can be interpreted as a secondary circulatory motion. The structure of the momentaneous "ow pattern changes continuously. The individual cells change their positions with time, thus their in"uence is averaged out in the long-time average. Further details pertaining to the di;erences in the

Fig. 14. Time series of the lateral and vertical velocities in an exemplary measuring point of the "at bubble column (Becker et al., 1999).

3178

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

Fig. 15. Typical instantaneous axial velocity versus time data in a bubble column obtained using LDA.

not appear to be similar to that in a two-dimensional column (Fig. 5), where clear periodicity in velocity is observed. From the foregoing discussion of the numerical results and the experimental observations, it can be seen that the "ow pattern in the three-dimensional column is much complex than the two-dimensional column. Further, there is a suLcient evidence of the existence of the multiple cells of di;erent sizes. Further, these cells are not stationary. In addition, turbulent eddies are also present with a wide size and velocity distribution. Therefore, a velocity time plot gives broad frequency distribution. In view of such a complexity, it is obvious that, substantial additional progress is needed for the CFD simulations for transient "ow structures in cylindrical columns. As an alternative the technique of multiresolution analysis, which has shown a potential for the resolution of coherent structures, will be described below. 5.2.5. P7eger and Becker (2001) P"eger and Becker (2001) used 0:288 m diameter and 2:6 m height of a bubble column with uniform sparging. The system was air–water and gas "ow rate was varied in the range of 0.15 –20 mm=s. The measurement (velocity–time series, radial pro>les of axial velocity) were made by LDA. The CFD simulation was carried out using commercial software CFX (version 4.3). A remarkable feature of the model is the inclusion of the bubble-induced turbulence. This modi>cation had a stronger in"uence and showed a positive impact on the correct simulation of both the radial pro>les of the time-averaged axial velocities but it deteriorates the prediction of the local and overall gas holdup. On comparing with the LDA measurements, the di;erence in the two time series clearly showed the inability of the simulations to include even high-frequency liquid phase variation in instantaneous velocity. The comparison of the time series through the probability distribution function of the liquid velocity showed good agreement in terms

of the range of velocity on X -axis while the exactness in the histogram showed the necessity of taking into account the high-frequency velocity "uctuations. A better agreement with experimental gas hold-up data was observed without the bubble-induced turbulence model. The disagreements in hold-up pro>le (steep experimental and "at simulated) suggest the need of incorporation of lift force that pushes the bubbles towards the center and the dispersion mechanism in the model. In all the geometries covered so far, i.e. (i) twodimensional columns with o;-center or central sparging, (ii) two-dimensional columns with uniform sparging and (iii) cylindrical bubble columns, a good amount of liquid circulation has been observed together with non-homogeneities in the hold-up structures. However, it must be pointed out that the transient simulations have been successful only for the simple case of two-dimensional columns with single-point sparging. As regards to the turbulent structures in three-dimensional columns, though some aspects have been covered by Kulkarni et al. (2001b), in order to reveal the intricacies in the turbulent "ow >eld it is required to continue the similar analysis.

6. Application of multiresolution analysis of velocity–time data for the characterization of turbulent &ow structures 6.1. Preamble Sokolichin and Eigenberger (1999) and Borchers et al. (1999) have simulated two-dimensional column with single-point sparging. The "ow developed in such a geometry has a characteristic feature that the coherent structures have bimodal frequency distribution. The authors have modeled the high-frequency part using three-dimensional k–$ formulation and predicted the low-frequency structures.

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

This methodology can be considered to be akin to LES. In fact, most of the coherent structures are basically turbulent eddies and for the resolution of all the eddies, it is well known that the DNS is useful. However, in case of a cylindrical bubble column, the frequency distribution of coherent structures is multimodal. Since the "ow is very complex, the applicability of DNS is diLcult considering the present status of the available computational power. Under these conditions, the principal problem is the identi>cation of the frequency range for which the transient knowledge is desired or the range for which turbulence modeling may be adopted. In this context the applicability of LES seems to be promising. However, for the use of LES, we need to know at least approximately, the eddy size distribution so that suitable >ltering can be done. In this context, the approach of multiresolution analysis is expected to be useful and some results have been described by Kulkarni, Joshi, Ravi Kumar, and Kulkarni (2001a). They have used the technique of wavelet transform and typical distribution of

Fig. 16. Schematic representation of local turbulent structures resolved by multiresolution analysis.

3179

coherent structures at a measurement location in the column is shown in Fig. 16. The multiresolution analysis helps in classi>cation of turbulent eddies on the basis of their frequency. The behavior of eddies in each class is described in detail by Kulkarni et al. (2001b). For some of these scales, phase-plane plots are also shown in Fig. 17. It can be seen that the low-frequency structures have distinctly di;erent behavior than the high-frequency structures. As discussed in Kulkarni et al. (2001b) the velocity–time data was >rst resolved into 13 scales using wavelet transforms. Further, from the knowledge of the scalewise power and their variation with respect to each other, scales corresponding to speci>c features in the data were identi>ed. Since the local information in time and frequency domain at each scale can be understood by comparing each variation in energy corresponding to data point (t) with respect to the mean energy for scale (j). This ratio represents the intermittency in the local turbulence and is known as the local intermittency measure and is denoted as LIMj (t). LIMj (t) was calculated for each scale and plotted with time for each data set. The distribution of the data points in each scale was analyzed using phase-plane plots. The phase-plane plots were obtained by plotting the LIMj (t) versus LIMj (t + dt) by Taken’s time-delay embedding method for an estimated delay-time dt, for t de>ned in real time. The method of choice for estimating a reasonable value for dt is to locate the time at which the >rst minimum in the delayed mutual information occurs. The >rst minimum signi>es the delay time where LIMj (t + dt) adds the maximum amount of information to LIMj on increasing the embedding dimension. It was observed that, for the >ner scales, a dense distribution of data points [LIMj (i)] resulted in a complex phase plot indicating the non-deterministic nature of the data. On the other hand, for the coarser scales (viz: j = 9–11) the plots were well developed and indicated the presence of large-scale structures in the continuous phase. The regularity in the plots was distinctly di;erent for the coarser scales than >ner ones. Since the large-size coherent structures are observed at low

Fig. 17. Phase-plane plots of scalewise resolved local liquid velocity.

3180

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

frequency, the data from these scales can be used for the further analysis of the large-scale structures in turbulence. The nature of LIMj (i) was hence used in choosing the scale where large-scale structures exist in the monitored data. For orthogonal WT of a signal and their subsequent scalewise inversion xˆj (i), the periodic features of the wavelet basis function in the scalewise reconstructed signals could be present in the xˆj (i) at coarser scales. The procedure of minimizing the in"uence of wavelet basis function is explained in Kulkarni et al. (2001b), and it facilitated a way to compare the magnitudes of the LIMj (i) at all positions i on an equivalent basis. Since the structures from every individual scale have their characteristic size, their sizes were also estimated for di;erent experimental conditions for di;erent data sets. It was seen that the sizes of coherent structures varied from scale to scale. Also, the variation in the sizes obtained from the di;erent velocity components for same scale at the same position and operating conditions were negligible, indicating the universality in the size distribution in di;erent directions. In summary, the technique of multiresolution analysis can provide a suitable >lter for LES and thus enables us to simulate the low-frequency structures, which usually appear in the form of circulation cells. 7. Conclusions (i) The experimental measurement of Becker et al. (1994) has triggered a systematic development in the CFD and has formed a benchmark for comparisons in a two-dimensional geometry and o;-center sparging. (ii) Eigenberger and co-workers have made sustained and pioneering e;orts to conclude the following: (a) two-dimensional solution procedures have a limited scope and three-dimensional solution procedures are essential, (b) laminar "ow assumption is very restrictive and a suitable turbulence model needs to be used, (c) care has to be exercised for proper grid structure and (d) Euler–Euler and Euler– Lagrange approaches give practically the same results when due attention is given to all the related issues. (iii) Conclusion (ii) is for two-dimensional columns with a single-point sparging. In this case, the "ow pattern consists of (a) low-frequency circulation cells and (b) high-frequency eddy motion with a clear demarcation in the range of frequencies. Excellent agreement has been shown between the predicted and experimental values of periodicity. Additional work is needed for the prediction of further details of circulation cells and the eddy di;usivity (or=and k and $) pro>les in the column. (iv) The transient "ow pattern in cylindrical columns consists of "ow structures with a broad range of

frequencies. For the resolution of all these transient structures, substantial additional work is needed in the development of CFD codes. (v) Low-frequency circulation cells have been observed experimentally by all the investigators in the two- as well as three-dimensional cylindrical column. The size, number and locations change in the column as pointed out by Joshi (1980, 1982). Further, in the cylindrical column, the frequency of such structures has a broad band. (vi) The principal outstanding problem is the clear identi>cation of the frequency range for which the transient knowledge is necessary or the range for which the turbulence modeling may be adopted. To begin with, the above problem may be investigated for the simulation of blending (or homogenization) and residence time distribution (RTD) for which experimental data are available in the published literature over a wide range of D; HD ; VG , sparger design and various gas–liquid systems. For such a identi>ed objective, the low-frequency coherent structures (circulation cells) need to be simulated for their transient behavior and the range of frequency necessary for the simulation of mixing and RTD needs to be established. The exercise should be continued for various design objectives. For this purpose, the technique of large eddy simulation appears to be useful. (vii) The technique of multiresolution analysis has a potential for the resolution of coherent "ow structures of all the frequencies. The technique needs further exploitation for its application in the LES model.

Notation CB CD Ck Ck=$ CL CS CV C$ C; BI C , C$1 , C$2 db dt D DG fk (r)

energy transfer coeLcient drag coeLcient in presence of other bubbles constant in bubble-induced turbulence model, dimensionless constant correspond either to the k- or $-equation. lift force coeLcient sub-grid scale model constant, dimensionless virtual mass coeLcient constant in bubble-induced turbulence model, dimensionless model constant bubble-induced turbulence, dimensionless constant in k–$ model, dimensionless bubble diameter, m delay time, s column diameter, m axial dispersion coeLcient, m2 =s probability density function characterizing the number of bubble of class k per unit volume

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

FD FG FK FL FP FV g G H HD i I j k LIM(i) m mb p r R Rb Re S t Tp u ub ud uG uG

uL

G uK uK

uK

uK

K uL um ur ur
ur

uslip Vb VG VK VL VS w W x z

drag force on a particle, N=m3 parameter de>ned in Eq. (2.18), N=m3 interfacial momentum transfer, N=m3 lift force on bubble, N=m3 force on bubble due to pressure gradient, N=m3 virtual mass part of the interaction force, N=m3 acceleration due to gravity, m=s2 turbulence production term, J=m3 s height of the column, m height of the dispersion in the bubble column, m time index kelvin impulse deformable body, Ns=m3 wavelet scale turbulent kinetic energy, m2 =s2 wave number local intermittency measure hold-up pro>le index mass of the bubble pressure, N=m2 radial coordinate, m column radius, m bubble radius, m Reynolds number characteristic >lter rate of strain, 1=s time, s ◦ mean axial temperature, C velocity vector, m=s bubble rise velocity, m=s drifting velocity, m=s velocity of the gas phase, m=s Reynolds stress of gas phase, N=m2 velocity of the respective Kth phase, m=s "uctuating velocity of Kth phase, m=s Reynolds stress of the Kth phase, N=m2 velocity of the liquid phase, m=s velocity of the gas–liquid mixture, m=s relative velocity, [(uG − uL ) − ud ], m=s instantaneous relative velocity, m=s "uctuating part of relative velocity, m=s slip velocity, m=s volume of the bubble, m3 super>cial gas velocity, m=s volume of bubbles in the cluster, m3 super>cial liquid velocity, m=s slip velocity, m=s probability density function width, m time series data axial coordinate

Greek letters * $ jG

>lter width, m turbulent dissipation rate, m3 =s3 gas phase holdup

jK jL  BI; L k L; e; t t; k T; L  t  "uid G L m % %f %&  P

3181

instantaneous holdup of Kth phase liquid phase holdup viscosity, Pa s bubble-induced turbulent viscosity, Pa s molecular viscosity of phase k; Pa s e;ective viscosity of liquid phase, Pa s turbulent viscosity, Pa s turbulent viscosity of phase k; Pa s liquid turbulent viscosity, Pa s molecular kinematic viscosity, m2 =s turbulent kinematic viscosity, m2 =s density, kg=m3 "uid density, kg=m3 density of gas phase, kg=m3 density of liquid phase, kg=m3 mixture density, kg=m3 surface tension, N=m gas phase dispersion number gas phase dispersion number stress tensor, N=m2 vorticity in liquid phase, s−1

Subscripts b e; G j K L

bubble e;ective gas phase wavelet scale phase liquid phase

Superscript ∧

denoised data

References Becker, S., De Bie, H., & Sweeney, J. (1999). Dynamic "ow behavior in bubble columns. Chemical Engineering Science, 54, 4929–4935. Becker, S., Sokolichin, A., & Eigenberger, G. (1994). Gas–liquid "ow in bubble columns and loop reactors: Part II. Comparison of detailed experiments and "ow simulations Chemical Engineering Science, 49, 5747–5762. Bel F’Dhila, R., & Simonin, O. (1992). Eulerian prediction of a turbulent bubbly "ow down stream a sudden pipe expansion. Proceedings of the sixth workshop on two phase 7ow predictions, Erlangel, Germany. Borchers, O., Busch, C., Sokolichin, A., & Eigenberger, G. (1999). Applicability of the standard k–$ turbulence model to the dynamic simulation of bubble columns. Part II: Comparison of detailed experiments and "ow simulation Chemical Engineering Science, 54, 5927–5935. Carver, M. B. (1982). A method of limiting intermediate values of volume fraction in iterative two-"uid computations. Journal of Mechanical Engineering Science, 24, 221–224. Carver, M. B. (1984). Numerical computation of phase separation in two "uid "ow. Journal of Fluids Engineering, 106, 147–153. Chen, J. J. J., Jamialahmadi, M., & Li, S. M. (1989). E;ect of liquid depth on circulation in bubble columns. Chemical Engineering Research and Design, 67, 203–207.

3182

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183

Chen, R. C., Reese, J., & Fan, L. S. (1994). Flow structure in a three– dimensional bubble column and three phase "uidized bed. American Institution of Chemical Engineers Journal, 40, 1093–1104. Davies, J. T. (1972). Turbulence phenomena. New York: Academic Press. Deb Roy, T., & Majumdar, A. K. (1981). Predicting "uid "ow in gas– stirred systems. Journal of Metals, 42– 47. Deen, N. G., Solberg, T., & Hjertager, B. H. (2001). Large eddy simulation of gas–liquid "ow in a square cross–sectioned bubble column. Chemical Engineering Science, 56, 6341–6349. Delnoij, E. (1999). Fluid dynamics of gas–liquid bubble columns. Ph.D. thesis, Twente University, Enschede, The Netherlands. Delnoij, E., Kuipers, J. A. M., & van Swaaij, W. P. M. (1997a). Computational "uid dynamics applied to gas–liquid contactors. Chemical Engineering Science, 52, 3623–3638. Delnoij, E., Lammers, F. A., Kuipers, J. A. M., & van Swaaij, W. P. M. (1997b). Dynamic simulation of dispersed gas–liquid two-phase "ow using a discrete bubble model. Chemical Engineering Science, 52, 1429–1458. Delnoij, E., Lammers, F. A., Kuipers, J. A. M., & van Swaaij, W. P. M. (1997c). Dynamic simulation of gas–liquid two-phase "ow: E;ect of column aspect ratio on the "ow structure Chemical Engineering Science, 52, 3759. Devanathan, N., Moslemian, D., & Dudukovic, M. P. (1990). Flow mapping in bubble columns using CARPT. Chemical Engineering Science, 45, 2285. Drahos, J., Zahradnik, J., Fialova, M., & Bradka, F. (1992). Identi>cation and modelling of liquid "ow structures in bubble column reactors. Chemical Engineering Science, 47, 3313–3320. Elghobashi, S., & Abou-Arab, T. W. (1983). A two-equation turbulence model for two-phase "ows. Physics of Fluids, 26(4), 931–938. Elmore, W. C., & Heald, M. A. (1969). Physics of waves. New York: McGraw–Hill. Freedman, W., & Davidson, J. F. (1969). Holdup and liquid circulation in bubble columns. Transactions of Institution of Chemical Engineers, 47, T251. Grevet, J. H., Szekely, J., & El-Kaddah, N. (1982). An experimental and theoretical study of gas bubble driven circulation systems. International Journal of Heat and Mass Transfer, 25(4), 487–497. Grevskott, S., Sannaes, B. H., Dudukovic, M. P., Hjarbo, K. W., & Svendsen, H. F. (1996). Liquid circulation, bubble size distributions, and solids movement in two- and three-phase bubble columns. Chemical Engineering Science, 51, 1703–1713. Grienberger, J., & Hofmann, H. (1992). Investigations and modelling of bubble columns. Chemical Engineering Science, 47, 2215–2220. Harlow, F. H., & Amsden, A. A. (1975). A numerical calculation of multiphase "uid "ow. Journal of Computational Physics, 17, 19–52. Hills, J. H. (1974). Radial non-uniformity of velocity and voidage in a bubble column. Transactions of Institution of Chemical Engineers, 52, 1–9. Hoomans, B. P. B., Kuipers, J. A. M., Briels, W. J., & van Swaaij, W. P. M. (1996). Discrete particle simulation of bubble and slug formation in a two dimensional gas–"uidized bed: A hard sphere approach Chemical Engineering Science, 51, 99. Ishii, M., & Zuber, N. (1979). Drag coeLcient and relative velocity in bubbly, droplet or particulate "ows. American Institution of Chemical Engineers Journal, 25, 843–855. Jakobsen, H. A. (2001). Phase distribution phenomena in two-phase bubble column reactors. Chemical Engineering Science, 56, 1049– 1056. Jakobsen, H. A., Sannaes, B. H., Grevskott, S., & Svendsen, H. F. (1997). Modelling of vertical bubble-driven "ows. Industrial and Engineering Chemistry Research, 36, 4052–4074. Jakobsen, H. A., Svendsen, H. F., & Hjarbo, K. W. (1993). On the prediction of local "ow structures in internal loop and bubble column reactors using a two-"uid model. Proceedings of the European symposium on computer aided process engineering, Vol. 2. (pp. 531– 536).

Johansen, S. T., & Boysan, F. (1988). Fluid dynamics in bubble stirred ladles: Part-II Mathematical modelling Metallurgical Transactions B, 19B, 755–764. Joshi, J. B. (1980). Axial mixing in multiphase Contactors—A uni>ed correlation. Transactions of Institution of Chemical Engineers, 58, 155–165. Joshi, J. B. (1982). Gas phase dispersion in bubble columns. Chemical Engineering Journal, 24, 213–216. Joshi, J. B. (1992). Comments on "ow mapping in bubble columns using CARPT. Chemical Engineering Science, 47, 508–509. Joshi, J. B. (2001). Computational "ow modelling and design of bubble column reactors. Chemical Engineering Science, 56, 5893–5933. Joshi, J. B., & Sharma, M. M. (1979). A circulation cell model for bubble columns. Transactions of Institution of Chemical Engineers, 57, 244–251. Kolmogorov, A. N. (1941). Dissipation of energy in locally isotropic turbulence. Doklady Akadamii Nauk SSSR, 32, 19–21. Kulkarni, A. A., Joshi, J. B., Ravi Kumar, V., & Kulkarni, B. D. (2001a). Application of multiresolution analysis for simultaneous measurements of gas and liquid velocities and fractional gas hold–up in bubble column using LDA. Chemical Engineering Science, 56, 5037–5048. Kulkarni, A. A., Joshi, J. B., Ravi Kumar, V., & Kulkarni, B. D. (2001b). Wavelet transform of velocity–time data for the analysis of turbulent structures in a bubble column. Chemical Engineering Science, 56, 5305–5315. Lai, K. Y. M., & Salcudean, M. (1987). Computer analysis of multi-dimensional, turbulent, buoyancy-induced, two-phase "ows in gas-agitated-liquid reactors. Computers and Fluids, 15(3), 281–295. Lapin, A., & Lubbert, A. (1994). Numerical simulations of the dynamics of two-phase gas–liquid "ows in bubble columns. Chemical Engineering Science, 49(21), 3661–3674. Lee, Y., & Korpela, S. A. (1983). Multicellular natural convection in a vertical slot. Journal of Fluid Mechanics, 126, 91–121. Liles, D. R., & Reed, W. H. (1978). A semi–implicit method for two-phase "uid dynamics. Journal of Computational Physics, 26, 390–407. Lin, T. J., Reese, J., Hong, T., & Fan, L. S. (1996). Quantitative analysis and computation of two-dimensional bubble columns. American Institute of Chemical Engineers Journal, 42, 301–318. Lopez de Bertodano, M., Lee, S. J., Lahey, R. T., & Drew, D. A. (1990). The prediction of two phase distribution phenomena using a Reynolds stress model. Journal of Fluid Engineers, 112, 107–113. Millies, M., & Mewes, D. (1994). Calculation of circulating "ows in bubble columns. Chemical Engineering Science, 50(13), 2093–2106. Mudde, R. F., & Simonin, O. (1999). Two and three dimensional simulation of a bubble plume using a two "uid model. Chemical Engineering Science, 54, 5061. Neti, S., & Mohamed, O. E. E. (1989). Prediction of two-phase turbulent jet and vertical pipe reactor. American Institution of Chemical Engineers Symposium Series, 86, 280–281. P"eger, D., & Becker, S. (2001). Modelling and simulation of the dynamic "ow behavior in a bubble column. Chemical Engineering Science, 56, 1737–1747. P"eger, D., Gomes, S., Gilbert, N., & Wagner, H. G. (1999). Hydrodynamic simulation of laboratory scale bubble columns: Fundamental studies of Eulerian–Eulerian modelling approach Chemical Engineering Science, 54, 5091–5099. Reese, J., & Fan, L. S. (1994). Transient "ow structure in the entrance region of a bubble column using particle image velocimetry. Chemical Engineering Science, 49(24B), 5623–5636. Sahai, Y., & Guthrie, R. I. L. (1982a). E;ective viscosity models for gas stirred ladles. Metallurgical Transactions B, 13B, 125–127. Sahai, Y., & Guthrie, R. I. L. (1982b). Hydrodynamics of gas stirred melts: Part I. Gas=liquid coupling Metallurgical Transactions B, 13B, 193–202. Sahai, Y., & Guthrie, R. I. L. (1982c). Hydrodynamics of gas stirred melts: Part II. Axisymmetric "ows Metallurgical Transactions B, 13B, 203–211.

J. B. Joshi et al. / Chemical Engineering Science 57 (2002) 3157–3183 Salcudean, M., Lai, K. Y. M., & Guthrie, R. I. L. (1985). Multidimensional heat, mass and "ow phenomena in gas stirred reactors. Canadian Journal of Chemical Engineering, 63, 51–61. Salcudean, M., Low, C. H., Hurda, A., & Guthrie, R. I. L. (1983). Computation of three-dimensional "ow and heat transfer in gas-agitated reactor. Chemical Engineering Communication, 21, 89–103. Sato, Y., & Sekoguchi, K. (1975). Liquid velocity distribution in two-phase bubble "ow. International Journal of Multiphase Flow, 2, 79. Schwarz, M. P., & Turner, W. J. (1988). Applicability of the standard k–$ turbulence model to gas stirred baths. Applied Mathematical Modelling, 12, 273–279. Simonin, O. (1990). Eulerian formulation for particle dispersion in turbulent two-phase "ows. In M. Sommerfeld, & P. Wennerberg (Eds.), Proceedings of the Dfth workshop on two-phase 7ow predictions. (pp. 156 –166). Erlangen, Germany. Smagorinsky, J. (1963). General circulation experiments with the primitive equations. Monthly Weather Review, 91, 99–165. Sokolichin, A., & Eigenberger, G. (1994). Gas–liquid "ow in bubble columns and loop reactors: Part—I. Detailed modelling and numerical simulation Chemical Engineering Science, 49, 5735–5746. Sokolichin, A., & Eigenberger, G. (1999). Applicability of the standard turbulence model to the dynamic simulation of bubble columns: Part— I. Detailed numerical simulations Chemical Engineering Science, 54, 2273–2284. Sokolichin, A., Eigenberger, G., Lapin, A., & Lubbert, A. (1997). Dynamic numerical simulation of gas–liquid two-phase "ows Euler=Euler versus Euler=Lagrange. Chemical Engineering Science, 52(4), 611–626. Sommerfeld, M., Deckwer S., & Kohnen, G. (1997). Time-dependent calculation of bubble columns based on Reynolds stress Navier Stokes equations with k–$ turbulence model. Proceedings of the second Japanese–German symposium on multi-phase 7ow. (pp. 323–334). Paper no. 3–1– 4.

3183

Stewart, H. B., & Wendro;, B. (1984). Two-phase "ow: Models and methods Journal of Computational Physics, 56, 363–409. Szekely, J., Wang, H. J., & Kiser, K. M. (1976). Flow pattern velocity and turbulence energy measurements and predictions in a water model of an argon-stirred ladle. Metallurgical Transactions B, 7B, 287–295. Technology Vision 2020, The US Chemical Industry. (1996). American Chemical Society, American Institute of Chemical Engineers, Chemical Manufacturers Association, Council for Chemical Research, Synthetic Organic Chemical Manufacturers Association, Washington DC. Thai Van, D., Minier, J., Simonin, O., Freydier, P., & Oliver, J. (1994). Multidimensional two-"uid model computation of turbulent dispersed two-phase "ows. ASME, Numerical Methods in Multiphase Flows, FED-185, 277–291. Thakre, S. S., & Joshi, J. B. (1999). CFD simulation of bubble column reactors: Importance of drag force formulation Chemical Engineering Science, 54, 5055–5060. Tolstoy, I. (1993). Wave propagation. New York: McGraw-Hill. Tzeng, J. W., Chen, R. C., & Fan, L. S. (1993). Visualization of "ow characteristics in a 2–D bubble column and three phase "uidized bed. American Institution of Chemical Engineers Journal, 39, 733–744. Whalley, P. B., & Davidson, J. F. (1974). Liquid circulation in bubble columns. Proceedings of symposium on two phase 7ow systems. Institution of Chemical Engineers Symposium Series, Vol. 38. (pp. 1–29). Paper J5. Institute of Chemical Engineers, London. Wu, Y., & Gidaspow, D. (2000). Hydrodynamics simulation of methanol synthesis in gas–liquid slurry bubble column reactors. Chemical Engineering Science, 55, 573–587. Yang, Z., Rustemeyer, U., Buchholz, R., & Onken, U. (1986). Pro>le of liquid "ow in bubble columns. Chemical Engineering Communication, 49, 51. Zehner, P. (1982). Momentum, mass and heat transfer in bubble columns. International Chemical Engineers, 26, 22–35.