Numerical Investigation of Gas Holdup and Phase Mixing in Bubble

May 25, 2011 - Chemical Engineering Science 2017 161, 350-359 ... Johan T. Padding , Niels G. Deen , E.A.J.F. (Frank) Peters , J.A.M. (Hans) Kuipers. ...
0 downloads 0 Views 1MB Size
ARTICLE pubs.acs.org/IECR

Numerical Investigation of Gas Holdup and Phase Mixing in Bubble Column Reactors Wei Bai,†,‡ Niels G. Deen,*,† and J. A. M. Kuipers† †

Department of Chemical Engineering and Chemistry, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands ABSTRACT: A discrete bubble model (DBM) has been used to study the overall gas holdup and the phase mixing in bubble column reactors. The EulerianLagrangian approach has the advantage that it is possible to study the overall gas holdup and the gas phase mixing in a direct way. Furthermore, by introducing tracer particles, also the liquid phase mixing can be studied in a Lagrangian manner. Comparisons suggest that the overall gas holdups obtained from the discrete bubble model agree very well with the correlations of Kumar et al., Heijnen and Van’t Riet, and Ruzicka et al. within the applied range of the superficial gas velocity. The gas phase dispersion coefficients from the simulation results agree pretty well with Wachi and Nojima’s correlation, which is derived on the basis of the recirculation theory. In addition, the turbulent diffusion coefficient calculated from the tracer particles velocities are very close to those from the literature within the applied range of the superficial gas velocity.

’ INTRODUCTION Bubble column reactors are frequently utilized in the chemical, petrochemical, biochemical, and metallurgical industries. A bubble column reactor is basically a liquid filled vessel with a gas distributor at the bottom. The gas is sparged through the distributor in the form of bubbles and comes into contact with the liquid. The resulting buoyancy driven flow creates strong liquid recirculation. Bubble column reactors have the advantage that little maintenance and low operating costs are required due to lack of moving parts and compactness. Moreover, they have excellent heat and mass transfer characteristics. There are some important parameters, such as gas holdup, gasliquid interfacial area, interfacial mass transfer coefficient, dispersion coefficient, and heat transfer coefficient and so on, which are essential to characterize, scale up, and design the bubble column reactors. For instance, the gas holdup gives the volume fraction of the phases present in the reactor. On the other hand, the gas holdup combined with knowledge of the mean bubble diameter allows determination of the specific interfacial area and the related volumetric mass transfer rates between the gas and liquid phase. Moreover, the gas holdup is usually used to identify the flow regime in bubble column reactors. During the past decades, extensive studies on the gas holdup have been carried out in several contexts, such as flow regime identification and factors that may influence the gas holdup in bubble column reactors, such as gas flow rate, liquid flow rate, geometry of the bubble column, operating conditions (such as pressure and temperature), physical properties of both phases, sparger type, and so on. The related work has been reviewed by many authors.17 Dispersion coefficients are of significance in chemical reactors as well. There are two main model reactors and related assumptions regarding prevailing state of mixing, ideally mixed flow, and plug flow in either single phase or multiphase chemical reactors.810 In the ideally mixed flow system, the composition of the reaction mixture is uniform and is typically used to describe either batch r 2011 American Chemical Society

reactors or continuously mixed tank reactors. On the other hand, mixing or diffusion of the reacting mixture does not occur in the flow direction in plug flow reactors. Real chemical reactors, however, may neither exhibit perfect mixed nor plug-flow behavior. For this reason, axial dispersion models have been developed to describe the deviations from these simplified model reactors. In addition, the concept of residence time distribution (RTD) has proven a powerful tool to diagnose and characterize nonideal reactors since Danckwerts.11 In bubble column reactors, investigations on the mixing of both the gas and liquid phase have been reviewed by several authors, such as Heijnen and Van’t Riet,1 Shah et al.,2 Joshi,12,13 and Deckwer and Schumpe.3 Computational fluid dynamics (CFD) has gained considerable interest in recent years in the field of chemical engineering.14 In addition to experimental routes, it provides an alternative way to study and characterize chemical reactors in detail. On the other hand, chemical reactor modeling is also crucial for scaleup and design of chemical reactors. EulerianLagrangian models are one of the available computational fluid dynamics tools in the field of chemical reactor modeling. Delnoij et al.15,16 used the EulerianLagrangian approach and then extended it to to a three-dimensional model to study the time-dependent motion of bubbles in bubble columns operating in the homogeneous regime. Sommerfeld et al.17 developed a bubble coalescence model on the basis of the Lagrangian approach, where bubbles are traced through the turbulent flow field along their trajectories. Darmana et al.18 applied a parallel algorithm to the EulerianLagrangian model embedding four-way coupling. Hu and Celik19 used the EulerianLagrangian approach to study gasliquid bubbly flow Special Issue: Nigam Issue Received: December 22, 2010 Accepted: May 24, 2011 Revised: May 20, 2011 Published: May 25, 2011 1949

dx.doi.org/10.1021/ie102557h | Ind. Eng. Chem. Res. 2012, 51, 1949–1961

Industrial & Engineering Chemistry Research

ARTICLE

Kato and Nishiwaki26 proposed the following correlation for the gas holdup in an airwater system:2 2:51ug 0:78 þ βug 0:8 ð1  eγ Þ

ð1Þ

β ¼ 4:5  3:5  2:548ϕ1:8

ð2Þ

γ ¼ 717ug 1:8 =β

ð3Þ

E¼ where

and 27

24

Figure 1. Gas holdup curve (reproduced from Fair et al. ).

in a flat bubble column by means of large-eddy simulation (LES) with two-way coupling. Darmana et al.20 and Gong et al.21 simulated the mass transfer with the EulerianLagrangian method, in which the mass transfer rate is calculated from the information of individual bubble directly. Gong et al.22 investigated the effect of microbubbles on the enhancement of the mass transfer and the residence time with the EulerianLagrangian model. Bai et al.23 adopted the EulerianLagrangian approach and studied the effect of gas sparger properties on the hydrodynamics and the gas phase mixing via analyzing the residence time distribution in a square bubble column. In this work, the EulerianLagrangian approach will be further discussed on predicting the parameters that have significant influence on the performance of bubble column reactors, such as gas holdup and mixing of both the liquid and gas phases. The simulation results will be compared with corresponding correlations reported in the literature. The structure of the present work is organized as follows. The next section contains the introduction including a literature review of the correlations about gas holdup and dispersion coefficient of both the gas and liquid phases, as well as the characterization of phase mixing in bubble columns. The subsequent sections give a brief description of the discrete bubble model and deal with the numerical implementation in the discrete bubble model. The next section presents the simulation details, the simulation results, and comparison with literature correlations. Finally, the conclusions are presented. Correlations of the Gas Holdup. In bubble column reactors, the gas holdup is an important parameter that can be used to characterize the hydrodynamics. It is defined as the overall volume fraction of the gas phase present in the form of bubbles. The gas holdup mainly depends on the superficial gas velocity. It can be used to estimate the interfacial area and is thus required to estimate the volumetric mass transfer rates between the gas and liquid phase. Fair et al.24 carried out measurements of the gas holdup and heat transfer in two commercial scale bubble columns with an airwater system. The diameters of the columns were 0.54 and 1.26 m, respectively. By comparing with results obtained in those columns with smaller diameter, they found that the gas holdup varies directly with the superficial gas velocity and is less in the column with a large diameter than in the column with a small diameter. However, when the diameter is larger than 0.54 m, the gas holdup is no longer influenced by the column diameter. The experimental data reproduced from Fair et al.,24 a linear fitted curve and the gas holdup curve from Shulman25 for a bubble column with a diameter of 0.12 m, are shown in Figure 1.

Kumar et al. correlated 382 gas holdup data points (ɛ < 0.35 and ug < 0.15 m/s) with the following equation: 

E ¼ 0:728u  0:485u 2 þ 0:0975u

3

ð4Þ

where u* is defined as a dimensionless gas velocity given by !0:25 Fl 2  u ¼ ug ð5Þ σΔFg According to the circulation theory and an assumption for the single bubble rise velocity, Heijnen and Van’t Riet1 proposed the following expression for the homogeneous flow regime: E ¼ 4ug

ð6Þ

which is consistent with the experimental observations by Krishna et al.28 On the basis of the concept of a characteristic turbulent kinematic viscosity in bubble columns, Kawase and Moo-Young29 developed a theoretical model for the gas holdup with Newtonian and non-Newtonian fluids. For Newtonian fluids, the equation is given as !1=3 ug 2 ð7Þ E ¼ 1:07 ϕg On the basis of hydrodynamic coupling between both phases, Ruzicka et al.30 developed a simple physical model for the homogeneous to heterogeneous regime transition in a bubble column. For the homogeneous regime, the formula for the gas holdup is given by pffiffiffi ug þ uT  B ð8Þ E¼ 2ð1 þ aÞuT where B ¼ ug 2 þ uT 2  2ð2a þ 1Þug uT

ð9Þ

For a bubble column with dimensions (H = 1 m, ϕ = 0.15 m), the two parameters, uT and a, are 0.22 m/s and 0.35, respectively, according to experimental data. Correlations for the Gas Phase Dispersion Coefficient. For the axial dispersion coefficient of the gas phase in a bubble column, Diboun and Sch€ugerl31 suggested: ur ϕ Bo ¼ ¼ 0:2 ð10Þ Dg where Bo is Bodenstein number and ur is the relative velocity between the gas and the liquid phase. Men’shchikov and Aerov32 studied the gas phase axial mixing in a bubble column (ϕ = 0.3 m) with the superficial gas velocity 1950

dx.doi.org/10.1021/ie102557h |Ind. Eng. Chem. Res. 2012, 51, 1949–1961

Industrial & Engineering Chemistry Research

ARTICLE

varying in the range of 0.0080.1 m/s and proposed the following correlation: Dg ¼ 1:47ug 0:72

ð11Þ

Towell and Ackermann33 correlated their own experimental data and those from the literature with a correlation as follows: Dg ¼ 19:7ϕ2 ug

ð12Þ

The superficial gas velocity applied for the correlation ranges from 0.009 to 0.13 m/s, whereas the diameter of the bubble column ranges from 0.09 to 1 m. Joshi13 reviewed six experimental investigations for the gas phase dispersion in bubble columns. On this basis, he proposed the following correlation, which covers all the experimental data presented in the six experimental studies: ug 2 ð13Þ E Heijnen and Van’t Riet1 suggested that the dispersion coefficient is constant in the homogeneous flow regime according to experimental data in the literature. They gave the following dispersion coefficient correlation in the heterogeneous flow regime for the airwater system: Dg ¼ 110ϕ2

Dg ¼ 78ðug ϕÞ1:5

ð15Þ

They also measured the gas phase dispersion coefficients in two bubble columns (ϕ = 0.2 and 0.5 m) using pulse experiments at superficial gas velocities in the range of 0.0290.456 m/s. Comparison between the theoretical correlation and the measurements revealed that the derived correlation is able to reflect the dependencies of the axial dispersion of the gas phase on the column diameter and the superficial gas velocity. Correlations for the Liquid Phase Dispersion Coefficient. Ohki and Inoue35 measured the axial liquid phase dispersion coefficient in batch type bubble columns of three different diameters (ϕ = 0.04, 0.08, and 0.16 m) and derived the correlations for homogeneous bubbly flows and heterogeneous flows respectively. The following correlation for the homogeneous bubbly flows was established on the basis of velocity distribution model: Dl ¼ 0:30ϕ2 ug 1:2 þ 170δ

Dl ¼ 1:23ϕ1:5 ug 0:5

ð16Þ

where δ is hole diameter of the perforated plate. Note that the units in the equation are non-SI units. That is, the unit of dispersion coefficient is cm2/s, superficial gas velocity has units of cm/s, and the hole diameter is in cm. Furthermore, the range of superficial gas velocity in his experiments ranges 0.02 and 0.25 m/s. Towell and Ackermann33 developed a correlation for dispersion coefficient as a function of the column diameter and

ð17Þ

The diameters of the bubble columns used in this study are 0.04 and 1.07 m, while the superficial gas velocity ranged from 0.016 to 0.13 m/s and 0.009 to 0.034 m/s, respectively.13 By measuring the dispersion coefficient in two tall bubble columns with (ϕ = 0.15 and 0.2 m), Deckwer et al.37 correlated his own experimental data with those in the literature and proposed an expression that is similar to the correlation proposed by Towell and Ackermann:33 Dl ¼ 0:678ϕ1:4 u0:3 g

ð18Þ

The superficial gas velocity used in the two bubble columns ranged from 0.004 to 0.13 m/s and 0.01 to 0.14 m/s, respectively. Hikita and Kikukawa38 studied the liquid dispersion coefficient in two bubble columns and presented a correlation for liquid dispersion coefficient: Dl ¼ ð0:15 þ 0:69ug 0:77 Þϕ1:25

ð19Þ

Baird and Rice39 adopted dimensional analysis to derive the dispersion coefficient based on isotropic turbulence model: Dl ¼ 0:35ϕ4=3 ðug gÞ1=3

ð14Þ

However, they also mentioned that the transition from the homogeneous to heterogeneous flow regime is difficult to distinguish due to geometric parameters, such as sparger location and uniformity of gas distribution. On the basis of the recirculation theory for a bubble column in the turbulent flow regime, Wachi and Nojima34 derived the following equation for the axial dispersion coefficient of the gas phase: Dg ¼ 20ϕ3=2 ug

the gas velocity:36

ð20Þ

Even though the turbulence in a bubble column is not necessarily isotropic, the correlation using an isotropic turbulence model as the basis gives a dimensionally consistent and simple equation. Zehner40 proposed a correlation based on a vortex cell model: Dl ¼ 0:368ϕ4=3 ðug gÞ1=3

ð21Þ

where the constant 0.368 is an adjustable parameter. Heijnen and Van’t Riet1 derived the liquid dispersion coefficient from either the liquid circulation velocity or the Peclet number according to Joshi.12 They found out that the two correlations are in good agreement with each other except for the constant C. The correlations have the following general form: Dl ¼ Cϕ4=3 ðug gÞ1=3

ð22Þ

where C equals either 0.33 or 0.36. By relating the axial dispersion coefficient to the longitudinal displacement of a fluid element during a certain period, Kawase and Moo-Young41 derived the following correlation valid for Newtonian liquids: Dl ¼ 0:343ϕ4=3 ðug gÞ1=3

ð23Þ

Note that most of the above correlations that are either obtained from experiments or derived theoretically in the literature show proportional dependencies of the liquid phase axial dispersion coefficient on the diameter of the bubble column and the superficial gas velocity. The only difference is that the powers of these two parameters and the constant vary slightly from each other. Phase Mixing. There are several factors that may induce longitudinal mixing in a bubble column, such as turbulent eddies, a nonuniform velocity distribution over the cross section of the reactor and molecular diffusion. In most of the practical cases, the former two factors dominate the molecular diffusion.8 Axial Dispersion Model. Axial dispersion models are often used to characterize the extent of phase mixing of chemical 1951

dx.doi.org/10.1021/ie102557h |Ind. Eng. Chem. Res. 2012, 51, 1949–1961

Industrial & Engineering Chemistry Research

ARTICLE

reactors. It assumes a diffusion-like process superimposed on plug flow. The unsteady state one-dimensional convection diffusion equation is given by9

where N represents the number of fluid elements. In turbulent flows, a turbulent diffusion coefficient can be expressed as42

DC D2 C DðUCÞ ¼D 2 Dt Dz Dz

DT ¼ ν 0 2 T L

ð24Þ

where C represents the species concentration, D is the axial dispersion coefficient, and U is a characteristic velocity that depends on the phase. For instance for the gas phase, the characteristic velocity is defined as U = ug/ɛ. The axial dispersion coefficient D is directly related to the process of mixing in reactors. For instance, a large D means rapid mixing while a small D implies slow mixing, whereas D = 0 corresponds to plug flow conditions. When the boundary conditions are known, i.e., openopen boundary conditions, the above equation can be solved analytically. From a pulse tracer experiment, the dispersion coefficient can then be obtained from the tracer curve at the outlet of the reactor, that is, from the mean and variance of the residence time distribution of the tracer. Residence Time Distribution. Since Danckwerts11 introduced the concept of residence time distribution, it has been used as an important tool in the analysis of chemical reactor performance. The residence time distribution E(t) normally has the property that the area under the curve is unity: Z ¥ EðtÞ dt ¼ 1 ð25Þ 0

In addition, the mean residence time tm and the variance σ2 of the residence time distribution can be calculated as follows: Z ¥ tEðtÞ dt ð26Þ tm ¼ 0

Z

¥

σ ¼ 2

ðt  tm Þ2 EðtÞdt

ð27Þ

0

Consequently, the axial dispersion coefficient can be obtained with the mean residence time tm and the variance σ2. For instance, the Peclet number Pe can be calculated in a closedclosed system:9 σθ 2 ¼

σ2 2 2 ¼  2 ð1  ePe Þ 2 Pe Pe tm



UH Pe

ð29Þ

where H is the height of the system. Lagrangian Description of the Mixing. In a Lagrangian view, the diffusion coefficient can be evaluated from the change in mean-square displacement of fluid elements z2 with time t:42 D¼

1 dz2 2 dt



0

where RL ðτÞ ¼

vðtÞ vðt þ τÞ ν02

ð34Þ

’ DISCRETE BUBBLE MODEL The discrete bubble model (DBM) used in this study was originally developed by Delnoij et al.15,16 In this model, the liquid phase hydrodynamics is represented by the volume-averaged continuity and NavierStokes equations while the motion of the each individual bubble is tracked in a Lagrangian way. Bubble Motion. The motion of each individual bubble is computed with Newton’s second law. The liquid phase contributions are taken into account by the interphase momentum transfer experienced by each individual bubble. For an individual bubble, the equations can be written as Fg V

dv ¼ dt

∑F

ð35Þ

dr ¼v ð36Þ dt The net force acting on each individual bubble is calculated by considering all contributing forces. It is assumed that the net force is composed of separate, uncoupled contributions due to gravity, pressure, drag, lift, virtual mass, and wall force, respectively:

∑F ¼ FG þ FP þ FD þ FL þ FVM þ FW

ð37Þ

The gravity force acting on a bubble in a liquid is given by FG ¼ Fg V g

ð38Þ

The bubble will experience a force due to buoyancy, inertial forces, and viscous strain. The pressure gradient reflects all of these forces as can be inferred from the NavierStokes equations. The far field pressure force incorporating contributions of all the above forces is given by FP ¼  V rP

ð39Þ

The drag force exerted on a bubble rising in a liquid is expressed as 1 FD ¼  CD Fl πd2 jv  ujðv  uÞ 8

ð30Þ

where the mean-square displacement z2 is calculated as " #2 1 N 1 N 2 z ¼ ðzi ðtÞ  zi ð0ÞÞ  ðzi ðtÞ  zi ð0ÞÞ ð31Þ N i¼1 N i¼1



where ν0 is the root-mean-square value of the fluctuating part of the instantaneous fluid velocity and the Lagrangian time scale TL is defined as Z ¥ RL ðτÞ dτ ð33Þ TL ¼

ð28Þ

The dispersion coefficient D can hence be determined from the Peclet number:

ð32Þ

ð40Þ

where the drag coefficient CD is determined according to Rusche43 which takes swarm effects into account: CD ¼ CD¥ ½expð3:64Rg Þ þ Rg 0:864  1952

ð41Þ

dx.doi.org/10.1021/ie102557h |Ind. Eng. Chem. Res. 2012, 51, 1949–1961

Industrial & Engineering Chemistry Research

ARTICLE

where the drag coefficient of a single bubble CD¥ is taken from Tomiyama et al.:44     16 Eo 0:687 48 8 ð1 þ 0:15Re CD¥ ¼ max min Þ; ; ð42Þ Re Re 3 Eo þ 4 where Re is the bubble Reynolds number, Re = (Fl|v  u|d)/μl. A bubble rising in a nonuniform liquid flow field experiences a lift force due to vorticity or shear in the flow field. The shear induced lift force acting on a bubble is usually written as45 FL ¼  CL Fl V ðv  uÞ  r  u

ð43Þ

where the lift coefficient is calculated according to Tomiyama et al.:46 8 > < min½0:288 tanhð0:121ReÞ;f ðEoH Þ EoH < 4 4 e EoH e 10 CL ¼ f ðEoH Þ > : 0:29 EoH > 10 ð44Þ where f ðEoH Þ ¼ 0:00105EoH  0:0159EoH  0:0204EoH þ 0:474 3

2

ð45Þ The modified E€otv€os number, EoH, is defined by using the maximum horizontal dimension of a bubble as a characteristic length as follows: ðFl  Fg ÞgdH 2

ð46Þ σ The maximum horizontal diameter of the bubble is obtained from the bubble aspect ratio E according to Wellek et al.:47 EoH ¼



dV 1 ¼ 1 þ 0:163Eo0:757 dH

ð47Þ

where dV is the maximum vertical diameter of the bubble and Eo is the E€otv€os number, Eo = (Fl  Fg)gd2)/σ. The relation between the above two diameters and the diameter of the bubble d in the discrete bubble modeling is as follows: d ¼ ðdV dH 2 Þ1=3

ð48Þ

Accelerating bubbles experience a resistance, which is described as the virtual mass force:45   Dv Du  FVM ¼  CVM Fl V ð49Þ Dt Dt

wall in that direction. Finally, the wall force coefficient CW is given by ( expð  0:933Eo þ 0:179Þ 1 e Eo e 5 CW ¼ ð51Þ 0:007Eo þ 0:04 5 < Eo e 33 Liquid Phase Dynamics. The liquid phase hydrodynamics are described by a set of volume-averaged equations, which consists of the continuity and the NavierStokes equations. The presence of the bubbles is reflected by the local volume fraction of the liquid phase Rl and the interphase momentum transfer rate Φ (i.e., the sum of the interfacial forces):

D ðRl Fl Þ þ r 3 ðRl Fl uÞ ¼ 0 Dt

ð52Þ

D ðRl Fl uÞ þ r 3 ðRl Fl uuÞ ¼  Rl rp  r 3 ðRl τl Þ þ Rl Fl g þ Φ Dt ð53Þ The liquid phase is assumed to be Newtonian, thus the stress tensor τl can be expressed as   2 T ð54Þ τl ¼  μeff ;l ðruÞ þ ðruÞ  Iðr 3 uÞ 3 where μeff,l is the effective viscosity. In the present model, the effective viscosity is composed of two contributions, the molecular viscosity and the turbulent viscosity: μeff ;l ¼ μL;l þ μT;l

ð55Þ

where μT,l is the turbulent viscosity (or eddy viscosity), which is determined from turbulence modeling of the liquid phase. In the present work, large eddy simulation (LES) is used to simulate the turbulence induced by the rising bubbles. Therefore, the above continuity and NavierStokes equations are resolved for the field representing the “large” eddies, while the effect of the subgrid part of the velocity representing the “small scales” is included through an eddyviscosity subgrid-scale model. The eddyviscosity model proposed by Vreman49 is used to calculate the eddy viscosity: sffiffiffiffiffiffiffiffiffiffi Bβ 2 ð56Þ μT;l ¼ 2:5Fl CS Rij Rij where Rij = ∂uj/∂xi, βij = Δm2RmiRmj and Bβ = β11β22  β122 þ β11β33  β132 þ β22β33  β232. Δi is the filter width in the i direction. Collision model. In this work, a hard sphere collision model developed by Hoomans et al.50 is adopted to describe the bubblebubble interaction. It mainly consists of two parts. The first part is processing the sequence of collisions and the second part is dealing with the collision dynamics. A detailed description about the model is referred to Darmana51 and Hoomans.52

where the D/Dt operators denote the substantiative derivatives pertaining to the respective phases. In the present work, bubbles are assumed to have a spherical shape and a virtual mass coefficient of 0.5 is used. Bubbles in the vicinity of a solid wall experience a force referred to as the wall force:48 " # 1 1 1 FW ¼  CW d 2  Fl jðv  uÞ 3 nz j2 nW ð50Þ 2 y ðL  yÞ2

’ NUMERICAL IMPLEMENTATION OF PHASE MIXING STUDY IN DBM

where nz and nW, respectively, are the normal unit vectors in the vertical and wall normal direction, L is the dimension of the system in the wall normal direction, and y is the distance to the

Residence Time Distribution of the Gas Phase. In EulerianLagrangian modeling of bubbly flows, each bubble is tracked individually according to Newton’s second law. Hence, the 1953

dx.doi.org/10.1021/ie102557h |Ind. Eng. Chem. Res. 2012, 51, 1949–1961

Industrial & Engineering Chemistry Research

ARTICLE

Figure 2. Sketch of the square bubble column.

residence time of each bubble in the system can be determined directly from the simulation and from this information, the axial dispersion coefficient can be determined according to the axial dispersion model. Tracer Particles of the Liquid Phase. Massless tracer particles have been used to visualize unsteady flows in computational fluid dynamics (CFD).5355 In the present work, tracer particles are used to study the mixing of the liquid phase. For the liquid phase, the NavierStokes equations are solved on an Eulerian grid and the liquid velocities are obtained on the faces of each computational cell. The velocity of a tracer particle v at a certain moment t can then be interpolated according to its position r(t) in the grid. Accordingly, after a time Δt, the new position of the tracer particle can be calculated as Z rðt þ ΔtÞ ¼ rðtÞ þ

t þ Δt

vðrðtÞ;tÞ dt

ð57Þ

t

Equation 57 can be evaluated using a numerical integration scheme. Lane54 has reported that the first-order integration scheme could lead to erroneous results and suggested a higherorder integration to be used. Hence, the RungeKuttaFehlberg method is used in the present simulations. Details of this method are given in the Appendix. Yeung and Pope56 have shown that the accuracy of the tracer particle trajectories depends on the accuracy of the interpolation scheme used to calculate the tracer particle velocities rather than on the time-stepping method. In the present work, two interpolation methods, i.e., trilinear interpolation and tricubic interpolation,57 are compared with each other. More information on these interpolation methods is given in the Appendix.

’ SIMULATION DETAILS AND RESULTS The bubble column studied here is shown in Figure 2. The cross-sectional area of the column is 0.15 m  0.15 m (W  D). The column is initially filled with water to a height of 0.6 m (H). Air is used as the dispersed phase and introduced into the column through a perforated plate at the bottom of the column. The material properties of both phases are taken according to their room temperature values. The bubble column is operating under atmospheric pressure. A perforated plate with 576 holes (24  24) is used as gas distributor. All the holes have a diameter of 1 mm and are positioned in the central region of the plate with a square pitch of 6.25 mm. A grid with (20  20  80) cells is adopted in the simulations. Five superficial gas velocities (ug = 0.005, 0.01, 0.015, 0.02, and 0.025 m/s) are adopted. The diameter of the bubbles in the simulations is kept as a constant, 0.004 m. Due to the fact that the perforated plate adopted has many holes and bubbles injected from each hole have very low velocities, the coalescence and breakup of bubbles can be neglected. The time step for the liquid phase equations is set to 1.0  103 s and for the bubbles it is set to 5.0  105 s. The total simulation time is 150 s. There are two sets of tracer particles being used in the simulations. The two sets of tracer particles are released uniformly over the cross-section of the bubble column. The set of tracer particles released from the bottom of the bubble column is denoted as “Tracer0” and the other set of tracer particles released from the top of the bubble column is denoted as “Tracer1”. Moreover, both sets of tracer particles are released simultaneously at the simulation time t = 30 s at which the flow in the bubble column is fully developed. The grid dependency was checked by comparing simulation results for different grids (15  15  45, 20  20  60, 30  30  90, and 40  40  120 cells) to experimental data.58 In this case a square column with a diameter of 0.15 m and a height of 0.45 m was used with a central inlet and a superficial gas velocity of ug = 0.005 m/s. More details on this case can be found in the work of Deen et al.58 The differences between the different grids are small, as can be seen in Figure 3. These may further be reduced by running the simulation for longer periods. For these simulations a total time of 200 s was used, where the first 20 s was discarded because of start up effects. On the basis of these results a computational grid was selected with 20 cells along the width and depth of the column for the simulations presented in this work. Gas Holdup. The instantaneous gas holdup is recorded at every 0.04 s during the simulation. Examples of gas holdup histories at a superficial gas velocity ug = 0.005 m/s and ug = 0.025 m/s are plotted in Figure 4. The plots show that the gas holdup increases monotonically at the beginning of the simulation. After a short time, the gas holdup reaches a peak and then starts to fluctuate during the simulation. Bubbles are injected into the bubble column continuously and rise in the vertical direction, which is reflected in the increasing initial part of the curves. Once the bubbles reach the top of the bubble column, most of them escape from the column directly. This leads to the small decrease of the gas holdup after the first peak. As soon as a liquid circulation pattern in the bubble column has developed, some bubbles may be trapped in the circulating liquid and remain in the column for a longer time. Hence, fluctuations in the gas holdup curves result. 1954

dx.doi.org/10.1021/ie102557h |Ind. Eng. Chem. Res. 2012, 51, 1949–1961

Industrial & Engineering Chemistry Research

ARTICLE

Figure 3. Time averaged velocity profiles of the vertical liquid velocity (left) and vertical liquid velocity fluctuations (right) in the vertical center plane of a square bubble column, at a height of 0.25 m employing different grids.

Figure 4. Instantaneous gas holdup [] at ug = 0.005 m/s (left) and ug = 0.025 m/s (right).

Table 1. Overall Gas Holdup ug [m/s] 0.005

0.01

0.015

0.02

0.025

ɛ []

0.0234

0.0471

0.0671

0.0916

0.1119

σ []

0.0004

0.0017

0.0040

0.0056

0.0060

To obtain the overall gas holdup in the bubble column for each superficial gas velocity, the instantaneous gas holdups are averaged after 30 s. The resulting overall gas holdup for each superficial gas velocity is listed in Table 1. From this table it can be seen that the standard deviation of the averaged gas holdup σ is quite small compared to the overall gas holdup. Furthermore, within the adopted range of the superficial gas velocity, the overall gas holdup is related to the superficial gas velocity with a multiplication factor of 4 approximately. In Figure 5, the overall gas holdups obtained from the simulations are compared with the literature correlations introduced earlier. It is seen that the simulation results agree with the correlations of Kumar et al.,27 Heijnen and Van’t Riet,1 and Ruzicka et al.30 quite well. These results indicate that the overall gas holdup increases nearly linearly with the superficial gas velocity with a factor of 4 in the adopted range of the superficial gas velocity. It can also be seen that the gas holdups obtained from the correlation of Shulman25 are higher than those obtained from others, which may be due to the diameter of the bubble column as reported by Fair et al.24 In addition, the correlations of Fair et al.24 and Kato and Nishiwaki26 are close to each other when the superficial gas velocity ug < 0.025 m/s. As the superficial gas velocity increases, the gas holdup obtained according to Kato and Nishiwaki26 increases faster than that obtained from Fair et al.24

Figure 5. Comparison of the computed gas holdup with literature correlations.

Moreover, Kawase and Moo-Young’s29 correlation possesses a more flat trend. Gas Phase Dispersion. The residence times of all individual bubbles in the bubble column are recorded during the simulation, and from this data set, the residence time distribution is determined. In Figure 6, the residence time distribution for each superficial gas velocity is shown. The plot shows a single peak in the residence time distribution, which means that most of the bubbles in the bubble column exit the column at the corresponding time. Moreover, the long tail of the curve indicates that some bubbles remain in the column for longer residence times due to reasons mentioned before. As the superficial gas velocity increases, it is found that the peak of the residence time distribution shifts to smaller residence time. That can be explained by the fact that at higher gas loading a stronger liquid circulation is created, leading to a 1955

dx.doi.org/10.1021/ie102557h |Ind. Eng. Chem. Res. 2012, 51, 1949–1961

Industrial & Engineering Chemistry Research

ARTICLE

Figure 6. Residence time distributions of the gas phase.

Figure 7. Comparison of the simulated gas dispersion coefficient with literature correlations.

Table 2. Mean and Variance of the Simulated Residence Time Distributions ug [m/s] 0.005

0.01

0.015

0.02

0.025

tm [s]

2.88

2.94

2.90

2.99

3.01

σ2 [s2] σθ2 []

0.81 0.10

1.35 0.16

2.06 0.25

2.77 0.31

3.38 0.37

Pe []

19.50

11.74

7.00

5.25

4.05

Dg [m2/s]

0.007

0.011

0.019

0.025

0.033

faster transport of the bubbles. In addition, the distribution becomes wider at higher superficial gas velocity. Relevant moments of the residence time distributions, such as the mean residence time tm and the variance σ2 are determined from the curves. The computed results are shown in Table 2. A priori the mean residence time can be estimated from the ratio of the column height and the terminal rise velocity, e.g., tm ≈ H/u¥ = 0.6/0.2 = 3 s. It is seen that the simulated mean residence time is very close to this value. By using the solution for a closed closed system, the Peclet number and the gas phase dispersion coefficients can be determined. From Table 2, it can be seen that the variance σ2 and the dimensionless variance σθ2 of the residence time distribution increase with increasing superficial gas velocity, which is reflected by the widespread of the residence time distribution. Furthermore, the dispersion coefficient increases with increasing superficial gas velocity. This implies that the extent of gas phase backmixing (i.e., bubbles getting trapped in liquid vortices and/or liquid down flow) becomes more pronounced at higher superficial gas velocity. The simulated dispersion coefficients of the gas phase will subsequently be compared with a number of correlations from the literature. Note that for the earlier introduced correlations, the correlation of Joshi13 involves the relation between the superficial gas velocity and the overall gas holdup in the bubble column. Due to the fact that the overall gas holdups from the simulations agree with the correlations in the literature quite well, we used the overall gas holdup data obtained from the simulations to calculate the gas phase dispersion coefficient from the correlation proposed by Joshi.13 A comparison is shown in Figure 7. The comparison suggests that the simulation results agree very well with the correlation of Wachi and Nojima,34 which is derived on the basis of the recirculation theory. The gas phase dispersion

coefficients from the simulations, however, are far below the correlation of Men’shchikov and Aerov.32 Men’shchikov and Aerov32 proposed the correlation according to the experimental data from a bubble column with the diameter of 0.3 m and did not consider the effect of the column diameter in the correlation. The dependence of the gas phase dispersion on the column diameter, however, is used in most of the correlations. Moreover, the effect of the column diameter on the gas phase dispersion coefficient has been proven by Mangartz and Pilhofer.59 In addition, the simulation results are close to the correlations of Diboun and Sch€ugerl,31 Towell and Ackermann,33 Joshi,13 and Heijnen and Van’t Riet1 at low superficial gas velocities (i.e., ug e 0.01 m/s), while the difference becomes larger when the superficial gas velocity increases. Liquid Phase Dispersion. The two methods for interpolating the Eulerian velocity to the Lagrangian velocity, tricubic interpolation, and trilinear interpolation are compared with each other. For instance, the mean-square displacement of the tracer particles and the autocorrelation of vertical component of the tracer particle velocities at superficial gas velocity ug = 0.005 m/s are shown in Figure 8. It can be seen that there are some minor differences between the two methods with respect to both the mean-square displacement and the autocorrelation, which suggests that the interpolation scheme only slightly influences the resulting Lagrangian time series. In the following discussions, only results obtained with the tricubic interpolation scheme are considered. The turbulent diffusion coefficient of the liquid phase in the bubble column is determined according to eq 32. The Lagrangian time scale TL is determined as the time for the autocorrelation to decrease to 1/e of its initial value.60 The results obtained from both sets of tracer particles are listed in Table 3. The turbulent diffusion coefficient determined from the tracer particles “Tracer0” is denoted as Dl0 and that determined from the tracer particles “Tracer1” is represented as Dl1. Moreover, for each set of tracer particles, the turbulent diffusion coefficient is determined four times and the average of these four samples is presented in Table 3 along with the standard deviation which are denoted as σ0 and σ1. The table shows that the turbulent diffusion coefficients obtained from the two set of tracer particles are quite close to each other. However, the standard deviation of the turbulent diffusion coefficient is relatively large, particularly at low superficial gas velocity (i.e., ug = 0.005 m/s). Moreover, it can be found 1956

dx.doi.org/10.1021/ie102557h |Ind. Eng. Chem. Res. 2012, 51, 1949–1961

Industrial & Engineering Chemistry Research

ARTICLE

Figure 8. Comparison of interpolation methods at ug = 0.005 m/s; mean-square displacement (left) and autocorrelation (right).

Table 3. Turbulent Diffusion Coefficient of the Liquid Phase ug [m/s] 0.005 2

0.01

0.015

0.02

0.025

Dl0 [m /s]

0.006

0.010

0.019

0.016

0.021

σ0 [m2/s]

0.002

0.004

0.002

0.002

0.003

Dl1 [m2/s]

0.008

0.011

0.020

0.018

0.023

σ1 [m2/s]

0.003

0.004

0.003

0.003

0.001

Figure 9. Comparison of the simulated liquid phase dispersion coefficient with literature correlations.

that the turbulent diffusion coefficient of the liquid phase nearly increases with the superficial gas velocity. The turbulent diffusion coefficient of the liquid phase obtained from the simulations are compared with the correlations in the literature and shown in Figure 9. According to the dispersion coefficient correlations of the liquid phase introduced earlier, the correlations of Deckwer et al.,37 Baird and Rice,39 Heijnen and Van’t Riet,1 Zehner,40 and Kawase and Moo-Young41 are very close to each other. All of them exhibit similar dependence on the column diameter and the superficial gas velocity, and the constants of the correlations are also quite similar. Hence, only two of them are used for the comparison in the plot (i.e., Heijnen and Van’t Riet’s1 correlation with the constant equal to 0.33 and Zehner’s40 correlation). Moreover, Towell and Ackermann’s33 correlation possesses a similar trend with the above correlations because they use a very similar form of the expression with only slightly different powers of the column diameter and the superficial gas velocity. Hikita and Kikukawa’s38 correlation reveals a weaker dependence on the superficial gas velocity than the others. Moreover, the dispersion coefficients calculated from Hikita and Kikukawa’s38 correlation

are higher than those obtained from the others at low superficial gas velocity (i.e., ug < 0.02 m/s). Ohki and Inoue’s35 correlation increases with the superficial gas velocity faster than the others and is larger than the others when ug > 0.02 m/s. It is clearly found that the turbulent diffusion coefficients obtained from the previous section are very close to the dispersion coefficients calculated from the correlations in the literature within the applied range of the superficial gas velocity.

’ CONCLUSIONS The present work utilized the discrete bubble model (DBM) to investigate the overall gas holdup and the mixing of the both phases in bubble columns. One of the advantages of the EulerianLagrangian approach is that the dispersed phase is treated as individual elements, i.e., bubbles that are tracked individually. The total number and the residence times of bubbles in the bubble column reactors can be determined conveniently during simulations and thus, the overall gas holdup and the residence time distribution of the gas phase can be obtained accordingly. Hence, the dispersion coefficient describing the extent of the gas phase mixing can be calculated. Furthermore, by introducing liquid phase tracer particles, it is possible to study the liquid phase mixing as well. In the first part of this work, the overall gas holdups from the discrete bubble model are obtained within the range from ug = 0.005 m/s to ug = 0.025 m/s. The results are then compared with correlations in the literature. The comparison reveals that the simulation results agree well with the correlations of Kumar et al.,27 Heijnen and Van’t Riet,1 and Ruzicka et al.30 The results suggest that the overall gas holdup increases nearly linearly with the superficial gas velocity with a scaling constant of 4 in the adopted range of the superficial gas velocity. Second, the residence time distributions (RTD) of the gas phase are determined in the residence times of bubbles from the discrete bubble model. A sharp peak in the residence time distribution is obtained, which means that most of the bubbles in the bubble column exit the column at the corresponding time. Moreover, the long tail of the curve suggests that some bubbles remain in the column for longer residence times. As the superficial gas velocity increases, one can find that the peak of the residence time distribution shifts to smaller residence time, which is attributed to the stronger liquid circulation. In addition, the distribution becomes wider when the superficial gas velocity is higher, which is reflected by the increase of the variance of the residence time distribution. The comparison with the gas dispersion correlations in the literature shows that the simulation results agree very well with 1957

dx.doi.org/10.1021/ie102557h |Ind. Eng. Chem. Res. 2012, 51, 1949–1961

Industrial & Engineering Chemistry Research

ARTICLE

Wachi and Nojima’s34 correlation which is derived on basis of the recirculation theory. Finally, the liquid phase mixing is studied through the introduction of liquid phase tracer particles. Two methods for the interpolation scheme of the Lagrangian velocity from the Eulerian velocity of the liquid phase are compared with each other. Only small differences have been found between the different methods. The turbulent diffusion coefficient is calculated from the Lagrangian velocities and compared with correlations in literature. The comparison shows that the obtained turbulent diffusion coefficients are very close to those from the literature within the applied range of the superficial gas velocity.

’ APPENDIX RungeKuttaFehlberg Method. The RungeKutta

Fehlberg method (denoted RKF45) is a numerical technique used to solve ordinary differential equations of the form dr ¼ vðrðtÞ;tÞ dt

rðt þ ΔtÞ ¼ rðtÞ þ

In addition, a better value for the solution is determined using a RungeKutta method of order 5: r0kþ1 ¼ rk þ

ð58Þ

Formal integration of eq 58 yields Z

Figure 10. Element for interpolation in three dimensions.

t þ Δt

vðrðtÞ;tÞ dt

ð59Þ

t

The RungeKuttaFehlberg method has a procedure to determine if the proper step size h is being used. At each step, two different approximations from both O(h4) and O(h5) methods for the solution are made and compared. If the two results are in close agreement, the approximation is accepted. If the two results do not agree to a specified accuracy, the step size is reduced. If the results agree to more significant digits than required, the step size is increased. A brief description is given below. When integrating eq 59 with the step size h using the RungeKuttaFehlberg method, let rk be the current position of a certain particle. Each step requires the use of the following six values: k1 ¼ hvðrk ;tk Þ   1 1 k2 ¼ hv rk þ k1 ;tk þ h 4 4   3 9 3 k3 ¼ hv rk þ k1 þ k2 ;tk þ h 32 32 8   1932 7200 7296 12 k1  k2 þ k3 ;tk þ h k4 ¼ hv rk þ 2197 2197 2197 13   439 3680 845 k1  8k2 þ k3  k4 ;tk þ h k5 ¼ hv rk þ 216 513 4104   8 3544 1859 11 1 k3 þ k4  k5 ;tk þ h k6 ¼ hv rk  k1 þ 2k2  27 2565 4104 10 2

The two approximations are then compared with each other. If the error of approximation falls within the appropriate error bound, the approximation will be accepted. Otherwise, an optimal step size will be computed based on the two approximations. Interpolation Methods. Trilinear interpolation is a fast and simple interpolation scheme for a three dimensional Cartesian grid. It approximates the value of a point within the rectangular grid cell linearly, using data on the corners of the cell. If point p is in the cell shown in Figure 10 and has the fractional offsets (dx, dy, dz) from the grid point p1, then f ðdx;dy;dzÞ ¼ f ðp1 Þð1  dxÞð1  dyÞð1  dzÞ þ f ðp2 Þdxð1  dyÞð1  dzÞ þ f ðp3 Þð1  dxÞdyð1  dzÞ þ f ðp4 Þdxdyð1  dzÞ þ f ðp5 Þð1  dxÞð1  dyÞdz þ f ðp6 Þdxð1  dyÞdz þ f ðp7 Þð1  dxÞdydz þ f ðp8 Þdxdydz

f ðdx;dy;dzÞ ¼

Then an approximation to the solution is made using a RungeKutta method of order 4: 25 1408 2197 1 k1 þ k3 þ k4  k5 216 2565 4101 5

ð61Þ

ð63Þ

In addition, Lekien and Marsden57 introduced a local tricubic interpolation scheme in three dimensions that is both C1 and isotropic. The algorithm is based on a specific 64  64 matrix that gives the relationship between the derivatives at the corners of the elements and the coefficients of the tricubic interpolant for this element. A brief description of this method is as follows. As shown in the Figure 10, the function f at the position p in the regular grid can be represented by

ð60Þ

rkþ1 ¼ rk þ

16 6656 28561 9 2 k1 þ k3 þ k4  k5 þ k6 135 12825 56430 50 55 ð62Þ

N

N

N

aijk dxi dyj dzk ∑ ∑ ∑ i¼0 j¼0 k¼0

ð64Þ

The coefficients aijk are determined in such a way that the function f is C1, that is, f is continuous and its 3 first derivatives are also continuous. 1958

dx.doi.org/10.1021/ie102557h |Ind. Eng. Chem. Res. 2012, 51, 1949–1961

Industrial & Engineering Chemistry Research

ARTICLE

The 64 coefficients aijk are stacked in a vector a by defining R1þiþ4jþ16k ¼ aijk " i; j; k ∈ 0; 1; 2; 3

ð65Þ

Similarly, the constraints on f and its derivatives are stacked in a vector b by defining 8 > if 1 e i e 8 f ðp Þ > > i > > D > > f ðpi8 Þ if 9 e i e 16 > > Dx > > > D > > > f ðpi16 Þ if 17 e i e 24 > > Dy > > > > D > > f ðpi24 Þ if 25 e i e 32 > > > Dz < 2 D ð66Þ bi ¼ if 33 e i e 40 f ðpi32 Þ > > DxDy > > > 2 > > > D f ðpi40 Þ if 41 e i e 48 > > > DxDz > > > D2 > > > f ðp Þ if 49 e i e 56 > > DyDz i48 > > > > D3 > > > f ðpi56 Þ if 57 e i e 64 : DxDyDz The derivatives of f can be computed and evaluated for the 8 points pi. This gives a linear system in the 64 unknown coefficients Ri: BR ¼ b

ð67Þ

The matrix B is 64  64 and has only integer entries. The determinant of B is 1. As a result, B is invertible and the coefficients Ri can be computed using the linear relationship: 1

R¼B b

ð68Þ

1

The matrix B is the core of the tricubic interpolator and can be found in Lekien and Marsden.57

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]. Present Addresses ‡

Department of Chemical and Biological Engineering, Iowa State University.

’ ACKNOWLEDGMENT This work was carried out at the University of Twente. ’ NOMENCLATURE C model coefficient []; concentration of a species [kg m3] Bo Bodenstein number [] d bubble diameter [m] D column depth m]; dispersion coefficient [m2 s1] f function [] Eo E€otv€os number, Eo = ((Fl  Fg)gd2)/σ [] E(t) residence time distribution [s 1] F force vector [N] g gravity acceleration [m s2] H column height [m]

I n p Pe r R Re S t tm u U v V W

unit tensor [] unit normal vector [] pressure [N m2] Peclet number [] position vector [m] bubble radius [m] bubble Reynolds number, Re = (Fl|v  u|d)/μl [] characteristic filtered strain rate [s1] time [s] mean residence time [s] liquid velocity vector [m s1] mean bubble velocity [m s1] bubble velocity vector [m s1] volume [m3] column width [ m]

Greek Letters

R Δ δ ɛ μ ϕ Φ F σ σ2 σ θ2 τ

local volume fraction [] subgrid length scale [m] hole diameter of perforated plate [m] gas holdup [] viscosity [kg m1 s1] column diameter [m] volume averaged momentum transfer due to interphase forces [kg m2 s2] density [kg m3] surface tension [N m1]; standard deviation variance [s2] dimensionless variance [] stress tensor [N m2]

Indices

b cell D eff g G i j l L P S T VM W

bubble computational cell drag effective gas phase gravity i direction j direction liquid lift; molecular viscosity pressure subgrid turbulent virtual mass wall

’ REFERENCES (1) Heijnen, J. J.; Van’t Riet, K. Mass transfer, mixing and heat transfer phenomena in low viscosity bubble column reactors. Chem. Eng. J. 1984, 28, B21–B42. (2) Shah, Y. T.; Kelkar, B. G.; Godbole, S. P.; Deckwer, W. D. Design parameters estimations for bubble column reactors. AIChE J. 1982, 28, 353–379. (3) Deckwer, W. D.; Schumpe, A. Improved tools for bubble column reactor design and scale-up. Chem. Eng. Sci. 1993, 48, 889–911. (4) Kantarci, N.; Borak, F.; Ulgen, K. O. Bubble column reactors. Process Biochem. 2005, 40, 2263–2283. (5) Ribeiro, C. P., Jr.; Lage, P. L. C. Gas-liquid direct-contact evaporation: A review. Chem. Eng. Technol. 2005, 28, 1081–1107. 1959

dx.doi.org/10.1021/ie102557h |Ind. Eng. Chem. Res. 2012, 51, 1949–1961

Industrial & Engineering Chemistry Research (6) Shaikh, A.; Al-Dahhan, M. H. A review on flow regime transition in bubble columns. International Journal of Chemical Reactor Engineering 2007, 5 (7) Joshi, J. B.; Pandit, A. B.; Kataria, K. L.; Kulkarni, R. P.; Sawarkar, A. N.; Tandon, D.; Ram, Y.; Kumar, M. M. Petroleum residue upgradation via visbreaking: A review. Ind. Eng. Chem. Res. 2008, 47, 8960–8988. (8) Westerterp, K. R.; van Swaaij, W. P. M.; Beenackers, A. A. C. M. Chemical Reactor Design and Operation, 2nd ed.; John Wiley & Sons Ltd.: New York, 1987. (9) Levenspiel, O. Chemical Reaction Engineering, 3rd ed.; John Wiley & Sons: New York, 1999. (10) Scott Fogler, H. Elements of Chemical Reaction Engineering, 4th ed.; Prentice Hall PTR: Englewood Cliffs, NJ, 2005. (11) Danckwerts, P. V. Continuous flow systems: Distribution of residence times. Chem. Eng. Sci. 1953, 2, 1–13. (12) Joshi, J. B. Axial mixing in multiphase contactors - A unified correlation. Trans. Inst. Chem. Eng. 1980, 58, 155–165. (13) Joshi, J. B. Gas phase dispersion in bubble columns. Chem. Eng. J. 1982, 24, 213–216. (14) Jakobsen, H. A. Chemical Reactor Modeling: Multiphase Reactive Flows, 1st ed.; Springer: Berlin, 2008. (15) Delnoij, E.; Lammers, F. A.; Kuipers, J. A. M.; van Swaaij, W. P. M. Dynamic simulation of dispersed gas-liquid two-phase flow using a discrete bubble model. Chem. Eng. Sci. 1997, 52, 1429–1458. (16) Delnoij, E.; Kuipers, J. A. M.; van Swaaij, W. P. M. A threedimensional CFD model for gas-liquid bubble columns. Chem. Eng. Sci. 1999, 54, 2217–2226. (17) Sommerfeld, M.; Bourloutski, E.; Br€oder, D. Euler/Lagrange calculations of bubbly flows with consideration of bubble coalescence. Can. J. Chem. Eng. 2003, 81, 508–518. (18) Darmana, D.; Deen, N. G.; Kuipers, J. A. M. Parallelization of an Euler-Lagrange model using mixed domain decomposition and a mirror domain technique: Application to dispersed gas-liquid two-phase flow. J. Comput. Phys. 2006, 220 (19) Hu, G.; Celik, I. Eulerian-Lagrangian based large-eddy simulation of a partially aerated flat bubble column. Chem. Eng. Sci. 2008, 63, 253–271. (20) Darmana, D.; Deen, N. G.; Kuipers, J. A. M. Detailed modeling of hydrodynamics, mass transfer and chemical reactions in a bubble column using a discrete bubble model. Chem. Eng. Sci. 2005, 60 (21) Gong, X.; Takagi, S.; Huang, H.; Matsumoto, Y. A numerical study of mass transfer of ozone dissolution in bubble plumes with an Euler-Lagrange method. Chem. Eng. Sci. 2007, 62 (22) Gong, X.; Takagi, S.; Matsumoto, Y. The effect of bubbleinduced liquid flow on mass transfer in bubble plumes. Int. J. Multiphase Flow 2009, 35, 155–162. (23) Bai, W.; Deen, N. G.; Kuipers, J. A. M. Numerical analysis of the effect of gas sparging on bubble column hydrodynamics. Ind. Eng. Chem. Res. 2011, 50, 4320–4328. (24) Fair, J. R.; Lambright, A. J.; Andersen, J. W. Heat transfer and gas holdup in a sparged contactor. Ind. Eng. Chem. Process Des. Dev. 1962, 1, 33–36. (25) Shulman, H. L.; Molstad, M. C. Gas-Bubble Columns for GasLiquid Contacting. Ind. Eng. Chem. 1950, 42, 1058–1070. (26) Kato, Y.; Nishiwaki, A. Longitudinal dispersion coefficient of a liquid in a bubble column. Int. Chem. Eng. 1972, 12, 182–187. (27) Kumar, A.; Degaleesan, T. E.; Laddha, G. S.; Hoelscher, H. E. Bubble swarm characteristics in bubble columns. Can. J. Chem. Eng. 1976, 54, 503–508. (28) Krishna, R.; Wilkinson, P.; Dierendonck, L. V. A model for gas holdup in bubble columns incorporating the influence of gas density on flow regime transitions. Chem. Eng. Sci. 1991, 46, 2491–2496. (29) Kawase, Y.; Moo-Young, M. Theoretical prediction of gas holdup in bubble columns with newtonian and non-newtonian fluids. Ind. Eng. Chem. Res. 1987, 26, 933–937. (30) Ruzicka, M. C.; Zahradnik, J.; Drahos, J.; Thomas, N. H. Homogeneous-heterogeneous regime transition in bubble columns. Chem. Eng. Sci. 2001, 56, 4609–4626.

ARTICLE

(31) Diboun, M.; Sch€ugerl, K. Eine blasens€aule mit Gleichstrom von Wasser und Luft II—Mischungvorg€ange in der Fl€ussigkeitsphase. Chem. Eng. Sci. 1967, 22, 161–169. (32) Men'shchikov, V. A.; Aerov, J. Longitudinal Mixing of Gas Phase in Bubble-Plate Reactors. Theor. Found. Chem. Eng. 1967, 1, 739–741. (33) Towell, G. D.; Ackermann, G. H. Axial mixing of liquid and gas in large bubble reactors. 2nd International Symposium on Chemical Reaction Engineering, Amsterdam , 1972, p B 3-1 (34) Wachi, S.; Nojima, Y. Gas-phase dispersion in bubble columns. Chem. Eng. Sci. 1990, 45, 901–905. (35) Ohki, Y.; Inoue, H. Longitudinal mixing of the liquid phase in bubble columns. Chem. Eng. Sci. 1970, 25, 1–16. (36) Majumder, S. K. Analysis of dispersion coefficient of bubble motion and velocity characteristic factor in down and upflow bubble column reactor. Chem. Eng. Sci. 2008, 63, 3160–3170. (37) Deckwer, W. D.; Burckhart, R.; Zoll, G. Mixing and mass transfer in tall bubble columns. Chem. Eng. Sci. 1974, 29, 2177–2188. (38) Hikita, H.; Kikukawa, H. Liquid-phase mixing in bubble columns: Effect of liquid properties. Chem. Eng. J. 1974, 8, 191–197. (39) Baird, M. H. I.; Rice, R. G. Axial dispersion in large unbaffled columns. Chem. Eng. J. 1975, 9, 171–174. (40) Zehner, P. Momentum, mass and heat transfer in bubble columns. part 2. axial blending and heat transfer. Int. Chem. Eng. 1986, 26, 29–35. (41) Kawase, Y.; Moo-Young, M. Liquid phase mixing in bubble columns with newtonian and non-newtonian fluids. Chem. Eng. Sci. 1986, 41, 1969–1977. (42) Brodkey, R. S. Mixing: Theory and Practicem; Academic Press Inc.: New York, 1972; Vol. 1. (43) Rusche, H. Computational fluid dynamics of dispersed twophase flows at high phase fractions. Ph.D. thesis, 2002 (44) Tomiyama, A.; Kataoka, I.; Zun, I.; Sakaguchi, T. Drag coefficients of single bubbles under normal and micro gravity conditions. JSME Int. J., Ser. B 1998, 41, 472–479. (45) Auton, T. R. The lift force on a spherical body in a rotational flow. J. Fluid Mech. 1987, 197, 241–257. (46) Tomiyama, A.; Tamai, H.; Zun, I.; Hosokawa, S. Transverse migration of single bubbles in simple shear flows. Chem. Eng. Sci. 2002, 57, 1849–1858. (47) Wellek, R. M.; Agrawal, A. K.; Skelland, A. H. P. Shape of liquid drops moving in liquid media. AIChE J. 1966, 12, 854–862. (48) Tomiyama, A.; Matsuoka, T.; Fukuda, T.; Sakaguchi, T. A simple numerical method for solving an incompressible two-fluid model in a general curvilinear coordinate system. Adv. Multiphase Flow 1995, 241–252. (49) Vreman, A. W. An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications. Phys. Fluids 2004, 16, 3670–3681. (50) Hoomans, B. P. B.; Kuipers, J. A. M.; Briels, W. J.; van Swaaij, W. P. M. Discrete particle simulation of bubble and slug formation in a two-dimensional gas-fluidised bed: A hard-sphere approach. Chem. Eng. Sci. 1996, 51, 99–118. (51) Darmana, D. On the multiscale modelling of hydrodynamics, mass transfer and chemical reactions in bubble columns. Ph.D. thesis, Enschede, 2006 (52) Hoomans, B. P. B. Granular dynamics of gas-solid two-phase flows. Ph.D. Thesis, Twente University, 2000. (53) Hin, A. J. S.; Post, F. H. Visualization of turbulent flow with particles, VIS '93: Proceedings of the 4th conference on Visualization '93; IEEE Computer Society: Washington, DC 1993; pp 4653. (54) Lane, D. A. Scientific visualization of large-scale unsteady fluid flows; IEEE Computer Society: Washington, DC, 1997; pp 125145. (55) Bauer, D.; Peikert, R.; Sato, M.; Sick, M. VIS ’02: Proceedings of the conference on Visualization ’02; IEEE Computer Society: Washington, DC, 2002; pp 525528. 1960

dx.doi.org/10.1021/ie102557h |Ind. Eng. Chem. Res. 2012, 51, 1949–1961

Industrial & Engineering Chemistry Research

ARTICLE

(56) Yeung, P. K.; Pope, S. B. An algorithm for tracking fluid particles in numerical simulations of homogeneous turbulence. J. Comput. Phys. 1988, 79, 373–416. (57) Lekien, F.; Marsden, J. Tricubic interpolation in three dimensions. Int. J. Numer. Methods Eng. 2005, 63, 455–471. (58) Deen, N. G.; Solberg, T.; Hjertager, B. H. Large eddy simulation of the gas-liquid flow in a square cross-sectioned bubble column. Chem. Eng. Sci. 2001, 56, 6341–6349. (59) Mangartz, K. H.; Pilhofer, T. Interpretation of mass transfer measurements in bubble columns considering dispersion of both phases. Chem. Eng. Sci. 1981, 36, 1069–1077. (60) Squires, K. D.; Eaton, J. K. Lagrangian and eulerian statistics obtained from direct numerical simulations of homogeneous turbulence. Phys. Fluids A 1991, 3, 130–143.

1961

dx.doi.org/10.1021/ie102557h |Ind. Eng. Chem. Res. 2012, 51, 1949–1961