Computational Fluid Dynamics Simulations in Bubble-Column

Ind. Eng. Chem. Res. 2005, 44, 1413-1423

1413

Computational Fluid Dynamics Simulations in Bubble-Column Reactors: Laminar and Transition Regimes Kalekudithi Ekambara and Jyeshtharaj B. Joshi* Institute of Chemical Technology, University of Mumbai, Matunga, Mumbai 400 019, India

In the present work, a computational fluid dynamics (CFD) model was developed to describe both of the extreme regimes (viscous and turbulent), including the transition regime. The second objective was to examine the extent to which CFD models are able to describe quantitatively the variation of εG with VG as a function of flow regimes. This study helps to underline the distinguishing characteristics of both regimes: homogeneous and heterogeneous. An extensive comparison of predicted mean axial liquid velocity profiles and fractional gas hold-up profiles with the experimental data has been presented. The agreement has been shown to be excellent. The CFD model has also been compared with the simplified analytical solutions. In a fully developed heterogeneous regime and under turbulent conditions, the reversal point for the axial velocity was found to be in the range of 0.7-0.75R. With a decrease in the Reynolds number (DVCFL/µL), the reversal point was found to shift toward the center up to 0.5R. Further, the CFD simulation was found to reveal a large number of characteristic features of flow in viscous and transition regimes. 1. Introduction Bubble-column reactors are now being extensively used in a wide variety of chemical and biochemical processes. The present design practice of the bubblecolumn reactors is still closer to an art than science because of the complexity of fluid mechanics. Over the past few years, there have been sustained efforts to improve our understanding of the governing equations of change (equations of continuity and motion) for twophase flows. Both Euler-Euler and Euler-Lagrangian approaches have been extensively used, and a considerable research effort has been focused on improving the computational fluid dynamics (CFD) simulations for reliable predictions. Recent publications have developed transient 3D models for analyzing the detailed flow structures in bubble columns. However, all of these papers are still at the exploratory stage, and substantial enhancement in computational power is needed to make progress in the understanding of detailed 3D hydrodynamics. Until then, 3D models are not suitable for developing relationships between the flow pattern and the design objectives. Further, the 3D transient simulations of real systems (cylindrical column, high VG, high HD/D ratio, and large column diameter) are computationally demanding and are very difficult. However, in the past, practically the entire simulation effort has concentrated on the turbulent regime, and there is a need to investigate the role of molecular transport properties such as viscosity on the flow pattern. In view of this, it was thought desirable to understand the hydrodynamics of bubble columns operating in laminar and transition flow regimes using CFD. The objective of the present work is to develop a CFD model to describe both of the extreme regimes (laminar and turbulent), including the transition regime, and is to examine the extent to which CFD models are able to describe quantitatively the variation of εG * To whom correspondence should be addressed. Tel.: + 9122-414 5616. Fax: +91-22-414 5614. E-mail: [email protected].

with VG as a function of flow regimes. This study also helps to underline the distinguishing characteristics of both regimes. 2. Literature Review The flow pattern in a bubble column was qualitatively described by deNever1 and by Freedman and Davidson,2 who identified the existence of two streams in the column: one heading upward driven by the buoyancy of gas bubbles and the other carrying the liquid down. Crabtree and Bridgwater3 studied theoretically and experimentally the liquid motion induced by a chain of bubbles in a viscous liquid. Experiments were performed in 50-, 76-, and 145-mm-i.d. columns, and the liquid height was 900 mm in all cases. The measurement location was between 400 and 500 mm from the column base. They showed a good comparison and developed a model of the flow supposing the gas to be equivalent to a line force acting vertically upward along the axis at the center of the column. Further, they found the maximum downward velocity to be around r/R ) 0.7 and the reversal point around r/R ) 0.45-0.5. For a very gentle bubbling rate (the superficial gas velocity was below 0.3 mm/s), the velocity distribution in the central plume was found to be of parabolic shape. Crabtree and Bridgwater showed a good agreement. Rietema and Ottengraf4 used a flow visualization technique by tracking tiny air bubbles. They also made observations similar to those of Crabtree and Bridgwater.3 Yoshitome and Shirai5 measured the velocities in the central plume by a spherical float. At relatively high superficial gas velocities (between 3 and 75 mm/s), they found the axial liquid velocity profile to be essentially flat. Because the diameter of the float was in the range of 0.167-0.334D, this experimental technique could not provide point velocities. Ulbrecht and Baykara6 measured both the axial plume velocity and the axial dispersion coefficient in non-Newtonian fluids. They measured the plume velocity by tracking the motion of a neutrally buoyant

10.1021/ie0492606 CCC: $30.25 © 2005 American Chemical Society Published on Web 01/19/2005

1414

Ind. Eng. Chem. Res., Vol. 44, No. 5, 2005

styrofoam sphere and the dispersion coefficient using a transient tracer response technique. They found the dispersion coefficient to vary linearly with the average axial plume velocity. However, these observations were restricted to very low gas holdup and to relatively small diameter columns. Despite the limitations, the study clearly suggests the importance of liquid circulation on axial mixing. Walter and Blanch7 measured axial liquid velocity profiles in bubble column (270-mm i.d. and 2.7-m height) using a hot-film anemometer. They observed that, at low VG, the boundary layer thickness was large and was found to control the liquid velocity profile. In the case of turbulent flow, the boundary layer thickness was found to be small and did not have any significant effect on the liquid velocity profile. They also reported the effect of the liquid velocity profile on the gas holdup and the liquid-phase axial dispersion. Ulbrecht et al.8 extended the work of Ulbrecht and Baykara6 to higher VG and hence high Reynolds number. They also examined the influence of the column geometry and the liquid-phase physical properties. They showed the applicability of the earlier model6 for both laminar and turbulent regimes of flow. In both cases, they observed a close agreement between the measured velocity and the model predictions in the bulk region. However, in the near-wall region, the model resulted in underpredictions. Deen et al.9 used a bubble column of a cross section (W × D) of 0.15 m × 0.15 m and a height of 1 m with uniform sparging. The noncoalescing system of an airaqueous salt solution was studied. The CFD simulation was carried out for a superficial gas velocity of 4.9 mm/s using the commercial software CFX (version 4.3). They compared the performance of the two turbulence models, k- and large-eddy simulation (LES), with and without inclusion of the bubble-induced turbulence. The effect of different interface forces on the simulated results was also investigated. The effect of the virtual mass force on the simulated results was found to be negligible. With only a drag force, the bubble plume did not show any transverse spreading. Out of the two turbulence models, they have found that LESs captured the strong transient movement of the bubble plume as observed in the experiment. The LES model resolved many more details of the flow, and large vortices were observed along the bubble plume. Because of the high turbulent viscosity, a quasi-stationary state was obtained for the k- model and transient details were not resolved but were implicitly contained in the turbulent kinetic energy. Further, they have observed that consideration of bubble-induced turbulence has a marginal effect on the predictions. van Baten and Krishna10 studied the hydrodynamics of bubble columns operating in the homogeneous and heterogeneous flow regimes. Simulations have been carried out for superficial gas velocities (VG) in the range of 6-80 mm/s spanning both regimes. The heterogeneous flow regime was assumed to consist of two bubble classes: “small” and “large” bubbles. The large bubbles were found to concentrate in the central core of the bubble column, whereas the small bubbles were distributed throughout the column. Further, they found that the small bubble holdup is practically constant in the heterogeneous flow regime and its value corresponds to that at the regime transition point.

Figure 1. Comparison between the prediction and experimental profiles of the axial liquid velocity available in the literature (refer to Table 4 for experimental details). Experimental data. Crabtree and Bridgwater:3 [, 50 mm i.d.; 2, 76 mm i.d.; ∆, 145 mm i.d.; ], Rietema and Ottengraf;4 0, Walter and Blanch.7 Correlations: (1) Crabtree and Bridgwater;3 (2) Ulbrecht and Baykara;6 (3) Walter and Blanch;7 (4) Ulbrecht et al.8

From the foregoing discussion, it can be seen that a number of models have been proposed in the literature3,6-8 for the investigation of the flow pattern in the laminar/transition regimes. Crabtree and Bridgwater3 and Walter and Blanch7 assumed a parabolic velocity profile in the central upflow region, whereas Ulbrecht6,8 and co-workers assumed a flat velocity profile. The predictive capability of these models against the experimental data is shown in Figure 1. It can be seen that the predictive capability of all the models have limited applicability when compared with the experimental data of liquid velocity profile. Further, it can be seen that the first two models heavily overestimate the liquid velocities in the central region and all of the models underestimate the velocities in the annular downflow region. The design parameters (heat-transfer coefficient, dispersion coefficient, etc.) depend strongly on the flow pattern; the model that gives good flow predictions can be expected to give good predictions of the residence time distribution as well. In view of the importance of the accurate description of the velocity profiles, it was thought desirable to develop a comprehensive CFD model for the laminar and turbulent regimes, including the transition regime. 3. Model Formulation for Flow Regimes For two-phase gas-liquid flow in bubble columns, the equations of continuity and motion (steady state) for a 2D cylindrical coordinate system can be represented in the following generalized form:

1 ∂ ∂ ∂ 1 ∂ (rεFvΦ)K + [εFuΦ)K ) (rµεΦ) + r ∂r ∂z ∂r r ∂r K 2 ∂ (εµΦ)K + SΦ,K (1) ∂z2

(

]

where Φ is the transport variable. K denotes the phase (K ) G or L), and SΦ is the source term for the respective dependent variable. Values of Φ and SΦ for different transport variables are given for laminar and turbulent flows in Tables 1 and 2, respectively. 3.1. Boundary Conditions. (i) Along the axis: axisymmetry in uG, uL, vG, vL, k, , and εG. (ii) Along the

Ind. Eng. Chem. Res., Vol. 44, No. 5, 2005 1415 Table 1. 2D Laminar Model for Gas-Liquid Flows Governing Equationsa governing equations written in a general form

(

conservation of mass axial velocity momentum

Φ 1 u

radial velocity momentum

v

FLr ) -CLεGFL(uG - uL)

[1r ∂r∂ r(v 1 ∂ ) -C ε F [ r(u r ∂r

FVz ) -CVεGFL FVr a

1 ∂ ∂ 1 ∂ ∂Φ (rεFvΦ)K + (εFuΦ)K ) rε r ∂r ∂z r ∂r ∂r

V G L

+ K

∂ ∂Φ εµ ∂r ∂z

(

)

K

+ SΦ,K

SΦ,K ) source terms ∂P + εKFKg ( F/DZεL ( F/VZεL ∂z ∂P -εK ( F/DrεL ( F/VrεL ( F/LrεL ∂z -εK

∂uL ∂r FDz ) -4.9 × 104εG(uG - uL)

] ]

∂ (v - vL) ∂z G ∂ (u - uL) G - uL) + ∂z G

G

)

- vL) +

FDr ) -4.9 × 104εG(vG - vL)

Terms marked with asterisks are force terms. Force terms are positive for the gas phase and negative for the liquid phase.

Table 2. 2D k-E Model for Gas-Liquid Flows Governing Equationsa governing equations written in a general form conservation of mass

1 ∂ ∂ ∂Φ ∂ ∂Φ (rεFvΦ)K + rεΓ + εΓ + SΦ,K r ∂r ∂z ∂r K ∂z ∂z K Φ σ Φ σf SΦ,K ) source terms 1 ∞ 1 to ∞ 1 ∂ µ1,K ∂εK ∂ µ1,K ∂εK r + r ∂r σf ∂r ∂z σf ∂z

(

)

(

)

(

u 1.0 1 to ∞

axial velocity momentum

-εK

) (

∂P 1 ∂ ∂v ∂u ∂ + εKFKg ( F/DzεL ( F/VzεL + rεµt + εµt ∂z r ∂r ∂z ∂z ∂z

[ (

v

1.0 1 to ∞

)

[ ( ) ( )] ( ) (

uK radial velocity momentum

)

∂ µt ∂ε 1 ∂ µt ∂ε r + r ∂r σf ∂r ∂z σf ∂r

+

K

µt ∂ε σf ∂z

K

∂P 1 ∂ ∂v ∂u ∂ rεµt + εµt ( F/DrεL ( F/VrεL ( F/LrεL + -εK ∂r r ∂r ∂r ∂z ∂z

[ (

)

µt,k ) 0.09FK(k2/);

G ) µt,L a

1.0 1.3

PB ) CB(FDrVS + FDzVS)

[1r ∂r∂ r(v 1 ∂ ) -C ∈ F [ r(v r ∂r

] ]

∂ (v - vL) ∂z G ∂ (u - uL) V G L G - vL) + ∂z G ∂vL 2 ∂uL 2 ∂vL ∂uL 1 vL 2 + 2 + + + ∂r r r ∂z ∂z ∂r

FVz ) -CV∈GFL FVr

k 

G

- vL) +

{ [( ) ( ) ( ) ] (

v 1 ∂ µt ∂ε ∂ µt ∂ε 2εµt 2 + vK r + r ∂r σ ∂r ∂z σf ∂r r K f εL(G + PB - FL)  εL [C1(G + PB) - C2FL] k ∂uL FLr ) -CLεGFL(uG - uL) ∂r FDz ) -4.9 × 104εG(uG - uL)

K

+ K

)

∂u 1 ∂ (rv) + r ∂r ∂z

( ) [ ( ) ( )] (

turbulent kinetic energy turbulent dissipation energy

)]

(

(

)]

K

K

)

∂u 1 ∂ (rv) + r ∂r ∂z

K

FDr ) -4.9 × 104εG(vG - vL)

)} 2

Terms marked with asterisks are force terms. Force terms are positive for the gas phase and negative for the liquid phase.

wall: k ) 0,  ) 2ν(dxk/dr)2, uL ) uG ) vL ) vG ) 0. (iii) At the inlet: the values of vG and vL are set to zero for steady-state flows, and uL is specified for continuous flows. (iv) At the top surface: the gradients of uG, vL, vG, k, and  are set to zero, and uL ) 0 for steady-state and continuous flows, where the gradient of uL is set to zero. For initiation of the numerical solution, εG and uG are specified at the inlet. At all of the other locations, uG, uL, vG, vL, and εG are taken to be zero at t ) 0. 3.2. Interphase Force Term. The interfacial forces arise because of momentum transfer across the interface. If the slip velocity is constant, the force is called the drag force. If the relative motion is unsteady, a virtual mass force prevails in addition to the drag force. When the liquid-phase flow pattern is nonuniform in the radial direction, the rising bubble experiences a radial (or lateral) force. The formulation of these forces has been discussed in detailed by Joshi.11 In the present

work, all three forces (drag, lift, and virtual mass forces) have been incorporated. The detailed formulations of these forces are given in Tables 1 and 2. 3.3. Quantitative Representation of a Bubble Column. Correspondence between Real Systems and Predicted Flow Patterns. The flow pattern in bubble columns mainly depends on the superficial gas velocity, column diameter, height-to-diameter ratio, sparger design, and nature of the gas-liquid system. The task seems to be formidable. However, it is possible to represent the combined effect of VG, D, HD, etc., through two parameters: the radial hold-up profile and the average bubble-rise velocity (VS). Any change in VG, D, HD, and physical properties affects the εG profile and VS. Such a relationship is usually expressed through the drift flux model of Zuber and Findlay12 and Ranade and Joshi13 and is given by the following equation:

VG/εjG ) C0(VG + VL) + C1

(2)

1416

Ind. Eng. Chem. Res., Vol. 44, No. 5, 2005

where

C0 )

〈εG(εGuG + VL)〉 〈εG〉〈εGuG + VL〉

+

〈εGεLuG〉 〈εG〉〈εGuG + VL〉

(3)

flow region) a net energy supply from the gas phase to the liquid phase. The rate of energy release from the gas phase to the liquid phase is given by the following equation:

E ) π/4D2(FL - FG)gHDεjL(VG - εjGVS)

and

C1 )

〈εGεLVS〉 〈εG〉

(4)

The parameters C0 and C1 are the drift flux constants. C0 represents the hold-up profile and C1 the bubble-rise velocity. The most fortunate characteristic feature of the bubble column is that the values of C0 and C1 are practically independent of the column diameter (of course, when D > 150 mm and the sparger region is exceeded). Therefore, for a given gas-liquid system, a few measurements of εjG with respect to VG and VL (over the range of interest) in a small-diameter column (∼150 mm) enable the estimation of C0 and C1. 3.4. Energy Balance. Consider a single bubble (of volume vB) rising in a pool of liquid. The bubble at any location has pressure energy and potential energy. As the bubble rises, the pressure energy decreases and the potential energy increases. The rate of change of energy is given by the following equation:

vB(FL - FG)g

dh ) v(FL - FG)gVS ) fDVS dt

(5)

where VS is the slip velocity and fD is the drag force. When the bubble rise is in the creeping flow regime, the energy released during the rise is dissipated in a viscous manner. In other words, the gas-phase energy gets directly converted into internal energy. When the bubble rise is in the turbulent flow regime (Re > 500), the energy released during the rise is first of all converted to turbulence in the liquid phase and finally into internal energy at the Kolmogorov scale. Therefore, the fraction of gas-phase energy used for generating liquid-phase turbulence depends on the bubble Reynolds number. There is one more aspect of the energy balance. In bubble columns, the bubble-rise velocity with respect to the wall can be much different from the slip velocity because of liquid circulation. When the actual velocity is higher than the slip velocity (as in the central upward region), the rate of energy release is different from that given by eq 5.

eR ) vB(FL - FG)g(VS + uL)

(6)

where uL is the liquid velocity at the location of bubble rise. The rate of energy converted into turbulence and spent in the viscous dissipation depends only on the slip velocity:

eD ) vB(FL - FG)gVS

(7)

Therefore, when uL is positive, the rate of energy release is higher than the dissipation. The extra energy is used for maintaining the liquid circulation. In contrast, when the bubble moves downward in the nearwall region, the liquid-phase energy is supplied to bubbles. The overall effect of upflow and downflow is (because of the gas hold-up profile, more bubbles are present in the upflow region as compared to the down-

(8)

The eq 8 assumes that the slip velocity is the same for all bubbles. The rate of energy converted into turbulence and spent in viscous energy dissipation is given by

E ) π/4D2(FL - FG)gHDεjGεjLVS

(9)

It may be, however, noted that the turbulence generated by bubbles has a scale corresponding to that of bubbles, whereas the turbulence generated in liquid circulation (using the energy given by eq 8) has a scale corresponding to that of column dimensions. Therefore, this gives rise to a very important question regarding the quantitative role of the bubble-generated turbulence to the transport phenomena. Sato and Sekoguchi14 and Theofanous and Sullivan15 have addressed this problem. It was assumed that a fraction of CB is used in the transport and this fraction forms a source term for the liquid-phase k (turbulent kinetic energy per unit mass) equation:

ES ) CBπ/4D2(FL - FG)gHDεjLεjGVS

(10)

where CB takes values between 0 and 1. Johansen and Boysan16 included the radial variation of the slip velocity and used CB in the range of 0.1-0.2. Later, Svendsen et al.17 used CB as a tuning parameter to match their experimental results. Kataoka et al.18 reported a detailed analysis of the extra terms due to the presence of the dispersed phase in the source of turbulent kinetic energy. Their analysis suggests that the extra generation of turbulence due to large bubbles is almost compensated for by the extra dissipation due to the small-scale interfacial structures. The studies by Hillmer et al.19 and Ranade20,21 also suggest a zero value for CB. Hjertager and Morud22 used a very low value of 0.02. The total energy released from the gas phase to the liquid phase is given by eq 8. Furthermore, the quantitative role of bubble-generated turbulence is given by eq 10. Therefore, the sum of eqs 8 and 10 gives the total energy received by the liquid phase and is expressed as

π E + ES ) D2(FL - FG)gHDεjL[VG + (CB - 1)εjGVS] 4 (11) It must be emphasized at this stage that, whatever may be the value of CB, the energy balance must be satisfied. When the k- model is used for the prediction of the flow pattern, we get radial and axial variation of  (energy dissipation rate per unit mass) as one of the answers. From this  field, the total energy dissipation rate can be calculated by suitable volume integration. Further details pertaining to the energy balance have been given by Joshi11 and Dhotre and Joshi.23 In all of the simulations reported in this work, the predicted energy dissipation rate equals the energy input rate given by eq 11. 4. Method of Solution The set of steady-state governing equations given in Tables 1 and 2 were solved numerically, which consisted

Ind. Eng. Chem. Res., Vol. 44, No. 5, 2005 1417 Table 3. Flow Regimes in Bubble Columns D (m)

VG (m/s)

εG

FL (kg/m3)

0.22 0.22 0.22 0.22 0.22 0.22 0.22 0.22

0.00005 0.0001 0.0005 0.001 0.005 0.01 0.02 0.04

0.0026 0.0038 0.0092 0.0134 0.0322 0.0469 0.0689 0.1015

998 998 998 998 998 998 998 998

0.22 0.22 0.22 0.22 0.22 0.22

0.0005 0.0005 0.0005 0.0005 0.0005 0.0005

0.00475 0.00509 0.00555 0.00693 0.00750 0.00924

998 998 998 998 998 998

Effect of the Viscosity 0.164 0.00016 0.097 9.7× 10-5 0.050 5× 10-5 0.009 9× 10-6 0.005 5× 10-6 0.001 1 × 10-6

0.05 0.1 0.15 0.22 0.4

0.0005 0.0005 0.0005 0.0005 0.0005

0.0092 0.0089 0.0091 0.0090 0.0092

998 998 998 998 998

νt (m2/s)

% νt

%ν

Re

1.06 × 10-6 1.51× 10-6 2.65× 10-6 6.65× 10-6 1.24× 10-5 9.30× 10-5 6.15× 10-4 2.31× 10-3

51.41 60.18 72.56 86.91 92.52 98.93 99.84 99.96

48.59 39.82 27.44 13.09 07.47 01.06 00.16 00.04

6977 9048 16546 21458 39237 50884 65989 85577

0.0398 0.0425 0.0462 0.0573 0.0616 0.0754

9.83× 10-8 4.4× 10-7 5.9× 10-7 1.39× 10-6 3.13× 10-6 2.65× 10-6

0.06 0.454 1.166 13.36 38.45 72.56

99.94 99.55 98.83 86.64 61.55 27.44

53 96 202 1396 2706 16546

Effect of the Column Diameter 0.001 1 × 10-6 0.0359 0.001 1 × 10-6 0.0508 0.001 1 × 10-6 0.0622 0.001 1 × 10-6 0.0753 0.001 1 × 10-6 0.1016

1.01× 10-6 1.25× 10-6 2.12× 10-6 2.65× 10-6 2.21× 10-5

50.19 55.45 67.91 72.56 95.66

49.81 44.55 32.09 27.44 03.34

1792 5070 9315 16546 40565

µL (kg/ms)

ν (m2/s)

VC (m/s)

Effect of the Superficial Gas Velocity 0.001 1 × 10-6 0.0318 0.001 1 × 10-6 0.0412 0.001 1 × 10-6 0.0754 0.001 1 × 10-6 0.0977 0.001 1 × 10-6 0.1787 -6 0.001 1 × 10 0.2318 0.001 1 × 10-6 0.3006 0.001 1 × 10-6 0.3898

of the following steps: (i) generation of a suitable grid system; (ii) conversion of governing equations into algebraic equations; (iii) selection of discretization schemes; (iv) formulation of a discretized equation at every grid location; (v) formulation of a pressure equation; (vi) development of a suitable iteration scheme for obtaining a final solution. A finite-control-volume technique of Patankar24 was employed for the solution of these equations. A staggered grid arrangement proposed by Patankar and Spalding,25 consisted of 50 × 100 grid points, with 50 grid points in the radial direction and 100 grid points in the axial direction. The same grid size was used for both flow conditions (laminar and turbulent regimes). The velocity components were calculated for the points that lie on the faces of the control volume, while all of the other variables were calculated at the center of the control volume. A power law scheme was used for discretization of the governing equations. A SIMPLE algorithm was used to solve the pressure velocity coupling term for the turbulent regime. The set of algebraic equations obtained after discretization was solved by TDMA. Relaxation parameters and internal iterations for the variables were tuned to optimize the balance between the convergence criteria (1.0 × 10-3) and the computation time. For this purpose, an in-house CFD code was developed. The detailed procedure was given by Ekambara and Joshi.26 To define the three viscous, transition, and turbulent flow regimes, the effect of the viscous term in the definition of the Reynolds number (DVCFL/µL) has been considered. Simulations were carried out for a wide range of superficial gas velocities, liquid viscosities, and column diameters. The simulation results are given in Table 3. In the viscous flow regime, the contribution of νt to momentum transfer is small, whereas in the turbulent flow regime, the contribution of ν is small. To investigate the transition Reynolds number, the value of volume average νt was estimated in all cases. The flow condition at which the ratio νt/ν was less than 0.1 was termed the viscous flow regime, and that at which it was greater than 10 was termed the turbulent flow regime.

5. Results and Discussion 5.1. Simulations of Flow Patterns. The predicted flow pattern should satisfy (i) the values of drift flux constants C0 and C1 for the given system and (ii) the energy balance. Mass and momentum balances for the phases get satisfied during the numerical integration of the equations of continuity and motion. The following stepwise procedure was used: (i) Parameter C1: The value of C1 corresponds to VS. Therefore, the value of C1 needs to be adjusted through the interface drag force:

FD ) εG(FL - FG)g ) -CDεG(uG - uL)

(12)

CD ) (FL - FG)g/VS

(13)

where

The value of VS was selected in such a way that the value of C1 obtained from the CFD flow pattern and eq 4 agrees with the value of C1 for the system under consideration. The values of C1 and VS are given in Table 6. The values of VS were found to be higher than the typical terminal rise velocities and found to increase with an increase in the superficial gas velocity. In the heterogeneous regime of bubble columns, upflow occurs in the central region and downflow in the near-wall region. Therefore, reduced pressure gradients prevail in the addition to static pressure gradients. The reduced pressure gradient exerts a drag force on the bubble that is upward in the central region. The usual drag force is downward. Therefore, for a given slip velocity (VS), the net drag force gets reduced in the field of the reduced pressure gradient. In other words, the slip velocity gets enhanced in the central upflow region. In the annular downflow region, the form drag is also downward, and hence the slip velocity gets reduced. Because the number of bubbles in the central region is far more than that in the annular region, the net effect of reduced pressures is in the enhancement of the slip velocity value. Further, all of the reasons that increase the liquid circulation also correspondingly increase the slip velocity. Thus, an increase in VG and/or C0 increases VS as

152/608

Ulbrecht and Baykara6

viscosity (mPa s)/ density (kg/m3)

100-540/-

102, 127, 16.8-25.9 1860-10200/ 152/1000 1001-1274

0-40

0.5-1.3

0.03-0.15 1008-1386/ 1260

range of VG (mm/s)

micro-hot film anemometer

photocolorimeter

injecting a dye tracer

measurement technique

()

{

2

1-

2

(

)

]

}

}

2

]

()

()

()

for turbulent flow r 2 r U ) 1.373 - 0.98 ln + 0.323 R R

[( )

for viscous flow r 2 r U ) 2.35 - 1 - 2.55 ln R R

{

for turbulent flow (-2γ2 + 1)(N + 4) r2 U)1- 2+ (-γNr2 + rN+2) γ [-(N + 4)γ2 + 4]2

[

2

N + 4 N+2 (r - r2) U) 1-r + N

2

(Rt ) ln r 1 + [( ) ] ln t (R) (R) r R

r2 R + ln -1 r R2

for slow flow

U)

U)

analytical solution of their mathematical model

(i) 1D; (ii) a flat velocity profile was assumed; (iii) the liquid is not moving radially or circumferentially; (vi) the flow was upward at the center and downward at the wall; (v) footnotes a-f

(i) 1D; (ii) the eddy or bulk viscosity is assumed to be a constant; (iii) wall effects are neglected; (iv) the flow was upward at the center and downward at the wall; (v) a parabolic velocity profile was assumed; (vi) the liquid is not moving radially or circumferentially; (vii) footnotes a-f

(i) 1D; (ii) the liquid is Newtonian and incompressible and the system is considered to be in the steady state; (iii) the liquid is not moving radially or circumferentially; (iv) the pressure gradient other than the hydrostatic gradient must exist to counteract the forcesexerted on the liquid; (v) the flow was upward at the center and downward at the wall; (vi) there is no net flow across a horizontal plane through the column; (vii) parabolic velocity profile; (viii) footnotes a-f (i) 1D; (ii) flat velocity profile; (iii) the liquid is not moving radially or circumferentially; (vi) the flow was upward at the center and downward at the wall; (v) footnotes a-f

assumptions and limitations of the model

a Gas-phase mass balance was not established. b Drag force was not considered. c Lift force was not considered. d Virtual mass force was not considered. e The predictions do not hold over a wide range of column diameters, superficial gas velocities, height-to-diameter ratios, and types of gas-liquid systems. f Energy balance was not established.

Ulbrecht et

al.8

Walter and Blanch7 270/2700

50, 76, 145/900

Crabtree and Bridgwater3

author

column diameter/ height (mm)

Table 4. Experiemental Details and Correlations of Previous Work

1418 Ind. Eng. Chem. Res., Vol. 44, No. 5, 2005

Ind. Eng. Chem. Res., Vol. 44, No. 5, 2005 1419 Table 5. Comparison of Radial Gas Holdup and Mean Axial Liquid Velocity Profiles: Experimental Data of Previous Work

researchers

column diameter/height (m)

Hills29

0.138/1.37

Menzel et al.30

0.6/3.45

Yu and Kim31

0.254/2.50

sparger

superficial meas. gas-liquid gas velocity location system (mm/s) HD/D

measurement techniques employed holdup liquid velocity

sieve plate: air-water hole diameter ) 0.4 mm, 61 nos. air-water

19-169

4.34

electroconductivity needle probe

Pavlov tube

12-96

4.86

perforated pipes: air-water 6 mm, 8 nos., hole diameter ) 1 mm, 78 nos.

10-140

0.475

electroconductivity hot film anemometry microprobe two-channel optical two Pt-Rh electrodes, fiber probe two U-shaped optical fibers, and a tracer injection

Table 6. Comparison between CFD Predictions and the Experimental Observationsa material balance authors Hills29 Menzel et al.30 Yu and Kim31

energy balance

VG (m/s)

εG

CB

CL

CV

C0

C1

LHS

RHS

0.0206 (0.019) 0.038 (0.038) 0.012 (0.012) 0.026 (0.024) 0.052 (0.048) 0.011 (0.01) 0.020 (0.02) 0.040 (0.04)

0.065 (0.067) 0.106 (0.107) 0.029 (0.029) 0.049 (0.048) 0.105 (0.103) 0.025 (0.025) 0.35 (0.034) 0.049 (0.048)

0.43 0.44 0.36 0.38 0.41 0.28 0.30 0.32

0.35 0.42 0.32 0.37 0.46 0.30 0.34 0.43

0.34 0.37 0.42 0.35 0.37 0.38 0.36 0.41

3.84 (2.0) 3.34 (2.0) 3.13 (2.5) 2.80 (2.5) 2.30 (2.5) 3.17 (2.1) 2.4 (2.1) 2.5 (2.1)

0.222 (0.25) 0.24 (0.25) 0.39 (0.40) 0.36 (0.40) 0.42 (0.40) 0.39 (0.40) 0.554 (0.73) 0.739 (0.73)

0.105 0.127 0.057 0.130 0.223 0.045 0.122 0.187

0.115 0.205 0.07 0.157 0.263 0.054 0.100 0.194

a The numbers in parentheses indicate the experimental values. LHS is the volume integral of the energy dissipation rate, and RHS is eq 11.

Figure 3. Comparison between the predictions and experimental data: (A) Walter and Blanch;7 (B) Rietema and Ottengraf.4 Correlations: (1) Crabtree and Bridgwater;3 (2) Walter and Blanch;7 (3) Ulbrecht et al.;8 (4) CFD.

Figure 2. Comparison between the predictions and experimental data of Crabtree and Bridgwater3 for different column diameters (D): (A) 5.0 cm; (B) 7.6 cm; (C) 14.5 cm. Correlations: (1) Crabtree and Bridgwater;3 (2) Walter and Blanch;7 (3) Ulbrecht et al.;8 (4) CFD.

compared with the characteristic VB∞. Similar observations have been made by Harper27 and Ishii and Zuber.28 (ii) Parameter C0: The value of C0 represents the steepness of the hold-up profile, which in turn depends on the radial movement of the bubbles. The bubbles move radially because of the radial lift force, turbulent dispersion, and radial slip. Therefore, the values of εjG and C0 depend on the combined effect of the radial force

1420

Ind. Eng. Chem. Res., Vol. 44, No. 5, 2005

Figure 4. Comparison between the predictions and experimental data of (A) Hills,29 (B) Menzel et al.,30 and (C) Yu and Kim.31

and turbulent dispersion. Therefore, an iterative procedure was implemented by the repeated use of the following two steps: (i) the value of CL was adjusted in such a way that the simulated value of εjG was equal to the experimental value, and (ii) in keeping CL from step i constant, the value of σf was varied in such a way that the simulated value of C0 was equal to the desired value of C0. Finally, the desired values of εjG and C0 were obtained at a particular combination of σf and CL. (iii) Average gas holdup (εjG): The average value of εjG of the CFD simulation must agree with the εjG of the real gas-liquid system. (iv) Energy balance: The value of CB was adjusted in such a way that the energy balance was satisfied. It may be noted that the values of CD, σφ, and CB were selected in such a way that the values of εjG, C0, and C1 correspond to the experimental data given in Table 6. The energy balance was also satisfied. Thus, a complete correspondence was established between the experimental and simulated systems. 5.2. Comparison between CFD Predictions and Experimental Observations. The CFD simulations were carried out for the experimental conditions of Crabtree and Bridgwater,3 Rietema and Ottengraf,4 Walter and Blanch,7 Hills,29 Menzel et al.,30 and Yu and Kim.31 Tables 4 and 5 give the details pertaining to the column diameter, column height, range of superficial gas velocity, range of viscosity, and measurement

techniques employed for the measurement of the mean axial liquid velocity and the fractional gas holdup. Figure 2 shows the comparison between the predictions of Crabtree and Bridgwater,3 Walter and Blanch,7 Ulbrecht et al.8 and the present CFD model with the experimental data of Crabtree and Bridgwater3 for different column diameters (50 < Re < 600). It can be seen that the models of Walter and Blanch7 and Ulbrecht et al.8 result in overpredictions. The model of Crabtree and Bridgwater3 gives a close agreement when r/R is greater than 0.2. The present CFD predictions can be seen to be in good agreement with the experimental data for all of the column diameters and over the entire range of radial positions. Further, it was observed that, in the CFD results, the reversal point changed from r/R ) 0.45 to 0.52 when the column diameter was increased from 50 to 145 mm (which means Re increases from 50 to 550). Figure 3 shows the comparison of the model predictions with the experimental data of Walter and Blanch7 and Rietema and Ottengraf.4 It can be observed that the predictions of Crabtree and Bridgwater3 and Ulbrecht et al.8 differ substantially from the experimental observations. The model of Walter and Blanch7 gives a close agreement with experiments. The CFD predictions can be seen to be in good agreement with the experimental data for both systems. Further, it can be observed from Figure 3 that the reversal point gets

Ind. Eng. Chem. Res., Vol. 44, No. 5, 2005 1421

Figure 6. Effect of the liquid viscosity on the average fractional gas holdup: (1) VG ) 0.03 m/s; (2) VG ) 0.10 m/s; (3) VG ) 0.15 m/s.

Figure 5. Effect of the liquid viscosity on the axial liquid velocity and gas holdup at a constant superficial gas velocity of 0.3 mm/s: (1) 1.08 Pa s; (2) 0.622 Pa s; (3) 0.350 Pa s; (4) 0.164 Pa s; (5) 0.097 Pa s.

changed from r/R ) 0.5 to 0.55 with a change in the column diameter from 220 to 270 mm. To confirm the validity of the present CFD model, it was thought desirable to establish a comparison at relatively high superficial gas velocity and/or low viscosity (Re > 30 000). For this purpose, the model predictions were compared with the experimental data of Hills,29 Menzel et al.,30 and Yu and Kim.31 The results are shown in Figure 4 and Table 6. The agreement between the predicted and experimental profiles can be seen to be good. The point of reversal in this case was found to be at r/R ) 0.7. The agreement over the wide range of D and VG and the types of gas-liquid systems ensures the applicability of the present model for the estimation of the flow pattern. In the foregoing simulations (Figures 2-4), the column diameter (50 mm < D < 600 mm), superficial gas velocity (0.03 mm/s < VG < 169 mm/s), and liquid viscosity (1 mPa s < µ < 10 200 mPa s) have been varied over a wide range. The Reynolds number at which the ratio (νjt/ν) is less than 0.9 was found to be 100 ( 10, whereas the ratio was found to be greater than 10 when Re was greater than 1000 ( 100. Thus, using the same criterion, the transition Reynolds number range can be seen to be 100 < Re < 1000. After the validity of the CFD model was confirmed, the effect of the liquid viscosity (Re e 100) was investigated on the axial liquid velocity and the gas hold-up profiles. The simulation results are shown in Figure 5. It can be observed that the reversal point of the axial velocity shifts toward the column center with an increase in the liquid viscosity even up to 0.5R. This means that, under the viscous flow regime (Re < 100), the upflow occurs only in the central 25% of the cross section and the downflow occurs in the remaining 75% of the cross section. In regards to the location of the

Figure 7. Effect of the superficial gas velocity on holdup and the axial liquid velocity: (1) VG ) 0.001 m/s; (2) VG ) 0.005 m/s; (3) VG ) 0.009 m/s; (4) VG ) 0.013 m/s; (5) VG ) 0.017 m/s; (6) VG ) 0.021 m/s; (7) VG ) 0.025 m/s; (8) VG ) 0.03 m/s; (9) VG ) 0.04 m/s; (10) VG ) 0.05 m/s.

maximum downward velocity, it was found to be in the range of 0.96-0.995R when Re > 10 000. With a decrease in Re, the location shifts toward the center and was found to be at r/R of 0.7 when Re < 100. This can be attributed to the increasing importance of viscosity in the momentum transfer and the concomitant increase in the thickness of the boundary layer. The influence of the liquid viscosity (100 < Re < 75 000) on the average gas holdup (εjG) is illustrated in Figure 6. It can be seen that the value of εjG increases with an increase in µ up to 3 mPa s. A further increase in µ up to 11 mPa s decreases εjG. However, when µ is increased beyond 11 mPa s, εjG was found to remain practically constant. These predictions are in agreement with a large number of experimental observations. The average bubble size increases with an increase in the viscosity because of the increasing coalescing nature of the liquid phase. Further, for the same bubble size, the

1422

Ind. Eng. Chem. Res., Vol. 44, No. 5, 2005

rise velocity decreases with an increase in the liquid viscosity. The overall effect of the above two factors gives maxima in the εjG-µ plots. The effect of the superficial gas velocity (12 000 < Re < 55 000) on the gas holdup and the axial liquid velocity was studied, and simulation results are shown in Figure 7. It was also found that the average liquid velocity increases with an increase in the superficial gas velocity and was found to be practically independent of the liquid height. However, the liquid velocities were found to decrease with an increase in the liquid viscosity. Thus, the CFD simulations reveal a large number of characteristic features of flow in the viscous and transition regions. 6. Conclusions (1) We have developed a CFD model to describe the hydrodynamics in the bubble column operating in both laminar and turbulent flow regimes, including the transition regime. (2) An extensive comparison of the predicted and experimental axial liquid velocities is presented. The agreement is shown to be excellent. The results also compare with other analytical models available in the literature. Further, it was observed that the reversal point of the axial liquid velocity was found to shift toward the column center with a decrease in the Reynolds number. Simultaneously, the locations of the maximum downward velocity also shifted toward the column center. This can be attributed to the increasing importance of the molecular viscosity in the momentum transfer. (3) The values of the average gas holdup were found to increase with an initial increase in the viscosity. However, a further increase in the viscosity resulted in a decrease in the average gas holdup. This is due to the fact that the average bubble size increases and the rise velocity decreases with an increase in the liquid viscosity. (4) The liquid velocity increases with an increase in the superficial gas velocity and was found to be practically independent of the liquid height. However, the liquid velocities were found to decrease with an increase in the liquid viscosity. Thus, the CFD simulations were found to reveal a large number of characteristic features of flow in the laminar and transition regions. Notation CFD ) computational fluid dynamics CB ) interface energy transfer factor CD ) drag force coefficient CL ) lift force coefficient CV ) virtual mass force coefficient C0, C1 ) drift flux constants C1 ) model parameter in the turbulent dissipation energy equation ()1.44) C2 ) model parameter in the turbulent dissipation energy equation ()1.92) D ) diameter of the column, m DL ) axial dispersion coefficient, m2/s eD ) frictional energy dissipation rate for a single bubble, W eR ) rate of energy release for a single bubble, J/s E ) parameter defined in eq 8, J/s ED ) frictional energy dissipation rate for all of the bubbles, W ES ) parameter defined in eq 10, J/s

fD ) drag force on a single particle in a fluidized bed, N/m2 FDr ) frictional force in the radial direction per unit volume of dispersion, N/m3 FDz ) frictional force in the axial direction per unit volume of dispersion, N/m3 FL ) lift force per unit volume of dispersion, N/m3 [)CLεLεGFL(uG - uL) (∂uL/∂r)] FVr ) radial virtual mass force per unit volume of dispersion, N/m3 [)CVεLεGFL{(1/r) (∂/∂r) r(vG - vL) + (∂/∂z) (vG - vL)}] FVz ) axial virtual mass force per unit volume of dispersion, N/m3 [)CVεLεGFL{(1/r) (∂/∂r) r(uG - uL) + (∂/∂z) (uG - uL)}] g ) gravitational constant, m/s2 G ) µt,L{2[(∂vL/∂r)2 + (vL/r)2 + (∂uL/∂z)2] + [(∂vL/∂z) + (∂uL/ ∂r)]2} H ) height of the column, m HD ) height of the gas dispersion, m k ) turbulent kinetic energy per unit mass, m2/s2 LES ) large-eddy simulation N ) constant P ) pressure, N/m2 PB ) CB[FDrVSr + FDzVS] r ) radial distance, m R ) column radius, m Re ) Reynolds number ) DVCµL/FL SΦ,K ) source term of phase K t ) time, s t ) locus of the liquid flow reversal u ) instantaneous axial velocity, m/s uc ) centerline liquid velocity, m/s uG ) axial component of the gas velocity, m/s uL ) axial component of the liquid velocity, m/s U ) dimensionless axial velocity (uL/uc) UL ) average liquid velocity, m/s v ) instantaneous radial velocity, m/s vB ) volume of a single bubble, m3 vG ) radial component of the gas velocity, m/s vL ) radial component of the liquid velocity, m/s VC ) circulation velocity, m/s VG ) superficial gas velocity, m/s VS ) slip velocity between gas and liquid, m/s VSr ) radial slip velocity between gas and liquid, m/s VB∞ ) terminal rise velocity, m/s W ) width of the column, m z ) axial distance along the column, m Greek Symbols ε ) fractional phase holdup εG ) fractional gas holdup εL ) fractional liquid holdup εjL ) time-averaged fractional liquid holdup  ) turbulent energy dissipation rate per unit mass, m2/s3 µK ) molecular viscosity of phase K, Pa s µt,K ) turbulent viscosity of phase K, Pa s ν ) molecular kinematic viscosity of the liquid, m2/s νt ) turbulent kinematic viscosity of the liquid, m2/s F ) density, kg/m3 FG ) density of the gas, kg/m3 FL ) density of the liquid, kg/m3 σφ ) turbulent Prandtl number for momentum transfer σf ) turbulent Prandtl number for bubble motion Subscripts G ) gas phase K ) phase, where K ) G, gas phase; K ) L, liquid phase L ) liquid phase

Literature Cited (1) deNever, N. Bubble Driven Fluid Circulations. AIChE J. 1960, 14, 222-226.

Ind. Eng. Chem. Res., Vol. 44, No. 5, 2005 1423 (2) Freedman, W.; Davidson, J. F. Hold-up and Liquid Circulation in Bubble Columns. Trans. Inst. Chem. Eng. 1969, 47, T251T262. (3) Crabtree, J. R.; Bridgwater, J. Chain Bubbling in Viscous Liquids. Chem. Eng. Sci. 1969, 24, 1755-1768. (4) Rietema, K.; Ottengraf, S. P. P. Laminar Liquid Circulation and Bubble Street Formation in a Gas-Liquid System. Trans. Inst. Chem. Eng. 1970, 48, T54-T62. (5) Yoshitoma, H.; Shirai, T. The Intensity of Bulk Flow in a Bubble Bed. J. Chem. Eng. Jpn. 1970, 3, 29-33. (6) Ulbrecht, J. J.; Baykara, Z. S. Significance of the Central Plume Velocity for the Correlation of Liquid Phase Mixing in Bubble Columns. Chem. Eng. Commun. 1981, 10, 165-185. (7) Walter, J. F.; Blanch, H. W. Liquid Circulation Patterns and their Effect on Gas Hold-up and Axial Mixing in Bubble Columns. Chem. Eng. Commun. 1983, 19, 243-262. (8) Ulbrecht, J. J.; Kawase, Y.; Auyeung, K. F. More on Mixing of Viscous Liquids in Bubble Columns. Chem. Eng. Commun. 1985, 35, 175-191. (9) Deen, N. G.; Solberg, T.; Hjertager, B. H. Large Eddy Simulation of Gas-Liquid Flow in a Square Cross-Sectioned Bubble Column. Chem. Eng. Sci. 2001, 56, 6341-6349. (10) van Baten, J. M.; Krishna, R. CFD Simulations of a Bubble Column Operating in the Homogeneous and Heterogeneous Flow Regimes. Chem. Eng. Technol. 2002, 25, 1081-1086. (11) Joshi, J. B. Computational Flow Modelling and Design of Bubble Column Reactors. Chem. Eng. Sci. 2001, 56, 5893-5933. (12) Zuber, N.; Findlay, J. A. Average Volumetric Concentration in Two-phase Flow Systems. J. Heat Transfer 1969, 87, 453. (13) Ranade, V. V.; Joshi, J. B. Transport Phenomena in Multiphase Systems. Proceedings of the Symposium on Transport Phenomena in Multiphase Systems; BHU: Varanasi, 1988; pp 113-196. (14) Sato, Y.; Sekoguchi, K. Liquid Velocity Distribution in Twophase Bubble Flow. Int. J. Multiphase Flow 1975, 2, 79-95. (15) Theofanous, T. G.; Sullivan, J. Turbulence in Two-phase Dispersed Flows. J. Fluid Mech. 1982, 116, 343-362. (16) Johansen, S. T.; Boysan, F. Fluid Dynamics in Bubble Stirred Ladles: Part-II. Mathematical Modelling. Metall. Trans. B 1988, 19B, 755-764. (17) Svendsen, H. F.; Jakobsen, H. A.; Torvik, R. Local Flow Structures in Internal Loop and Bubble Column Reactors. Chem. Eng. Sci. 1992, 47, 3297-3304.

(18) Kataoka, I.; Besnard, D. C.; Serizawa, A. Basic Equation of Turbulence and Modelling of Interfacial Transfer Terms in GasLiquid Two-phase Flow. Chem. Eng. Commun. 1992, 118, 221236. (19) Hillmer, G.; Weismantel, L.; Hofmann, H. Investigations and Modelling of Slurry Bubble Columns. Chem. Eng. Sci. 1994, 49, 837-843. (20) Ranade, V. V. Flow in Bubble Columns: Some Numerical Experiments. Chem. Eng. Sci.1992, 47 (8), 1857-1869. (21) Ranade, V. V. Modelling of Turbulent Flow in a Bubble Column Reactor. Chem. Eng. Res. Des. 1997, 75, 14-23. (22) Hjertager, B. H.; Morud, K. Computational Fluid Dynamics Simulation of Bioreactors. CFD Simul.1993, 47-61. (23) Dhotre, M. T.; Joshi, J. B. Two-dimensional CFD Model for Prediction of Pressure Drop, Heat Transfer Coefficient in Bubble Column Reactors. Chem. Eng. Res. Des. 2004, 82, 689707. (24) Patankar, S. V. Numerical heat transfer and fluid flow; McGraw-Hill: New York, 1981. (25) Patankar, S. V.; Spalding, D. B. A Calculation Procedure for Heat, Mass and Momentum Transfer in Three-dimensional Parabolic Flows. Int. J. Heat Mass Transfer 1972, 15, 126-129. (26) Ekambara, K.; Joshi, J. B. CFD Simulation of Mixing and Dispersion in Bubble Column Reactors. Chem. Eng. Res. Des. 2003, 81, 987-1002. (27) Harper, J. F. On Bubbles Rising Inline at Large Reynolds Numbers. J. Fluid Mech. 1970, 41, 751-758. (28) Ishii, M.; Zuber, N. Drag Coefficient and Relative Velocity in Bubbly, Droplet or Particulate Flows. AIChE J. 1979, 5, 843. (29) Hills, J. H. Radial Nonuniformity of Velocity and Voidage in the Bubble Column. Trans. Inst. Chem. Eng. 1974, 2, 1-9. (30) Menzel, T.; Weide, T.; Staudacher, O.; Onken, U. Reynolds Stress Model for Bubble Column Reactor. Ind. Eng. Chem. Res. 1990, 29, 988-994. (31) Yu, Y. H.; Kim, S. D. Bubble Properties and Local Liquid Velocity in the Radial Direction of Co-current Gas-Liquid Flow. Chem. Eng. Sci. 1991, 46, 313-320.

Received for review August 14, 2004 Revised manuscript received December 19, 2004 Accepted December 23, 2004 IE0492606