CFD Simulation of Bubble Column Reactor Using Population Balance

Ind. Eng. Chem. Res. 2008, 47, 8505–8516

8505

CORRESPONDENCE CFD Simulation of Bubble Column Reactor Using Population Balance Kalekudithi Ekambara,† Kumar Nandakumar,† and Jyeshtharaj B. Joshi*,‡ Department of Chemical and Materials Engineering, UniVersity of Alberta, Edmonton, AB, Canada T6G 2G6, and Institute of Chemical Technology, UniVersity of Mumbai, Matunga, Mumbai, India 400 019

In this paper, we have presented a comprehensive analysis of the development of flow pattern in a bubble column reactor by the introduction of a population balance equation combined with the three-dimensional two-fluid model (Reynolds stress model). The multiple size group (MUSIG) model has been used to account for the nonuniform bubble size distribution in a gas-liquid mixture. The coalescence and breakage effects of the gas bubbles are modeled according to the coalescence by the random collision driven by turbulence and wake entrainment while for bubble breakage by the impact of turbulent eddies. Local radial distributions of the gas hold-up, Sauter mean bubble diameter, axial liquid velocity, turbulent kinetic energy, turbulent energy dissipation rate, and Reynolds stresses for superficial gas velocity of 20 mm/s are compared against experimental data in a bubble column reactor. The development of flow pattern were examined at six axial locations H/D ) 0.2, 1.4, 2.6, 3.9, 5.0, and 6.2. Good quantitative agreement with the experimental data is obtained with three different models (i.e., k-ε, RSM with constant bubble size, and RSM with population balance). The model prediction shows better agreement with the experimental data with population balance than constant bubble diameter predictions. 1. Introduction Bubble columns are intensively utilized as multiphase contactors and reactors in chemical, petrochemical, biochemical, and metallurgical industries. They are used especially in chemical processes involving reactions such as oxidation, chlorination, alkylation, polymerization, and hydrogenation; in the manufacture of synthetic fuels by gas conversion processes; and in biochemical processes such as fermentation and biological wastewater treatment. Bubble column reactors owe their wide application area to a number of advantages they provide in both design and operation as compared to other reactors. First of all, they have excellent heat and mass transfer characteristics, meaning high heat and mass transfer coefficients. Little maintenance and low operating costs are required because of the lack of moving parts and compactness. The operation of bubble column reactors is affected by global operating parameters such as gas superficial velocity, operating pressure and temperature, and liquid height. Many hydrodynamic variables that influence bubble column performance are: gas holdup distribution, bubble breakup, coalescence and dispersion rate, bubble rise velocity, bubble size distribution, gas-liquid interfacial area distribution, gas-liquid mass/heat transfer coefficients, and the extent of liquid phase backmixing.1-5 To develop design tools for engineering purposes, a large amount of research has been carried out in the area of computational fluid dynamics (CFD) modeling of bubble column reactors. However, most of the CFD studies carried out use single mean bubble size.6-11 To generate a simulation result that resembles available data, the “mean” bubble size is customarily adjusted by a trial-and-error procedure and the value so chosen often is far from reality. Available estimates of the local interfacial area are based on predicted hold-up and assumed * Corresponding author. E-mail: [email protected]. Tel.: 91-22-2414 0865. Fax: 91-22-2414 5614. † University of Alberta. ‡ University of Mumbai.

mean bubble size, which introduces error in the calculation. To overcome this, knowledge of bubble size distribution is essential. Further, this also helps in modeling the drag force, which is one of the key closures. Population balance is a well-established method in computing the size distribution of the bubbles and accounting for the breakage and coalescence effects in bubbly flows. This approach is concerned with maintaining a record for the number of bubbles, whose presence or occurrence may dictate the behavior of the system under consideration. Lately, mounting interest on population balances have resulted in a number of significant development especially toward modeling bubbly flows. The most important issue of the computational fluid dynamics-population balance model (CFD-PBM) coupled model is therefore its ability to predict the hydrodynamics in the heterogeneous regime where bimodal bubble size distributions are observed, especially in high superficial gas velocities. However, most reported results by the CFD-PBM coupled model are either only for the homogeneous regime12-14 or for both the homogeneous and the heterogeneous regimes but do not successfully predict bimodal bubble size distribution in the heterogeneous regime.15-18 The present status on population balance modeling of bubbly flows has been examined by Jakobsen et al.19 Some of the recent attempts using LES9,20,21 and DNS22,23 of bubbly flow are encouraging. Nevertheless, this approach will require certain time to be able to handle the complexities in bubbly turbulent flows. In this work, an attempt has been made to combine population balance with computational fluid dynamics (CFD) for the case of gas-liquid flow in the bubble column reactor. The flow hydrodynamics is simulated considering the bubble coalescence and break-up. The gas hold-up, Sauter mean bubble diameter, axial liquid velocity, turbulent kinetic energy, turbulent energy dissipation rate, and Reynolds stress have been predicted at different height to diameter ratios. Further, the simulations were carried out with k-ε and Reynolds stress models (RSM) with

10.1021/ie071393e CCC: $40.75  2008 American Chemical Society Published on Web 10/01/2008

8506 Ind. Eng. Chem. Res., Vol. 47, No. 21, 2008

the constant bubble size (i.e., 6 mm). The models predictions were compared with experimental data.

N

∑ Ω(V : V )n

BB )

j

i

(6)

j

j)i+1

2. Mathematical Modeling 2.1. Population Balance Model. The type of flow encountered in bubble columns is referred to as polydispersed multiphase flow, i.e., flow in which the dispersed phase covers a broad range of size groups. One of the attributes of polydispersed multiphase flow is that the different sizes of the dispersed phase interact with each other through the mechanisms of breakup and coalescence. To deal with this type of flow, a population balance equation must be formulated. Population balance is a well-established method in computing the size distribution of the dispersed phase and accounting for the breakage and coalescence effects in bubbly flows. A general form of the population balance equation is ∂ni + ∇ · (uGni) ) BB + BC - DB - DC (1) ∂t where uG is the gas velocity, ni represents the number density of size group i, and terms on the right-hand side BB, BC, DB, and DC are, respectively, the “birth” and “death” due to breakup and coalescence of bubbles. The bubble number density ni is related to the gas volume fraction ∈G by ∈Gfi ) niVi

(2)

where fi represents the volume fraction of bubbles of group i and Vi is the corresponding volume of a bubble of group i. 2.1.1. Bubble Break-up Model. Bubble break up is analyzed in terms of bubble interactions with turbulent eddies. The turbulent eddies increase the surface energy of the bubbles through deformation. Break up occurs if the increase in surface energy reaches a critical value. The break-up of bubbles in turbulent dispersions employs the model developed by Luo and Svendsen.24 Binary break-up of the bubbles is assumed and the model is based on the theories of isotropic turbulence. For binary breakage, a dimensionless variable describing the sizes of daughter drops or bubbles (the breakage volume fraction) can be defined as fBV )

Vi d3i d3i ) 3) 3 V d d + d3 i

(3)

j

where di and dj are diameters (corresponding to Vi and Vj) of the daughter bubbles in the binary breakage of a parent bubble with diameter d (corresponding to volume V). The value interval of the breakage volume fraction is between 0 and 1. The breakup rate of bubbles of volume Vj into volume sizes of Vi ()VfBV) can be obtained as

( )

Ω(Vj : Vi) ε )C 2 (1 - ∈G)nj dj

1/3

∫

1

ζmin

(

)

12cfσ (1 + ζ)2 exp dζ 11⁄3 ζ βFLε2/3dj5/3ζ11/3 (4)

where ζ ) λ/dj is the size ratio between an eddy and a particle in the inertial subrange and consequently ζmin ) λmin/dj; C and β are determined, respectively, from fundamental consideration of drops or bubbles breakage in turbulent dispersion systems to be 0.923 and 2.0 in Luo and Svendsen;24 and cf is the increase coefficient of surface area given by cf ) [fBV2/3 + (1 - fBV)2/3 - 1]

(5)

The birth rate of group i bubbles due to break-up of larger bubbles is

The death rate of group i bubbles due to break-up into smaller bubbles is N

DB ) Ωini with Ωi )

∑Ω

ki

(7)

i)1

2.1.2. Bubble Coalescence Model. Bubble coalescence is modeled by considering bubble collision because of turbulence buoyancy and laminar shear. The Prince and Blanch25 coalescence model was used in this work. Collisions between bubbles take place by three mechanisms, which arise from random motion of bubbles due to turbulence (θijT), differences in rising velocity induced by the buoyancy force (θBij ), and laminar shear (θLS ij ). The calculations showed that laminar shear collisions are negligible because of the low superficial gas velocities considered in this investigation. The coalescence rate of bubbles in the bubble size groups i and j is expressed as

( )

χ ) (θTij + θBij )exp -

tij τij

(8)

where τij is the contact time for two bubbles given by (dij/2)2/3/ε1/3 and tij is the time required for two bubbles to coalesce having diameter di and dj estimated to be {(dij/2)3Fl/ 16σ}1/2 ln(h0/hj). The equivalent diameter dij is calculated as suggested by Chesters and Hoffman:26 dij ) (2/di + 2/dj)-1. The parameters h0 and hj represent the film thickness when collision begins and critical film thickness at which rupture occurs, respectively. The values of the above parameters depend mainly on the physical properties of the liquid phase and have been experimentally computed for the air-water system. According to Prince and Blanch,25 for air-water systems, experiments have determined h0 and hj to be 1 × 10-4 m (Kirkpatrick and Lockett27) and 1 × 10-8 m (Kim and Lee28), respectively. It has been proven, that even 10-fold variation of the ratio h0/hj causes less than 10% difference in dS value (Pohorecki et al.29). The turbulent collision rate θijT for two bubbles of diameter di and dj is given by π θTij ) [di + dj]2(uti2 + utj2)1⁄2 4

(9)

where the turbulent velocity ut in the inertial subrange of isotropic turbulence30 is ut ) 1.4ε1⁄3d1⁄3

(10)

The buoyancy contribution to collision frequency is modeled as θBij ) FCBSij|uij - uri|

(11)

where FCB is a calibration factor and uri )

2.14σ + 0.505gdi Fcdi

(12)

The birth rate of group i due to coalescence of group k and group l bubbles is N

BC )

N

∑∑

1 χ nn 2 k)1 l)1 i,kl i j

(13)

The death rate of group i due to coalescence with other bubbles is

Ind. Eng. Chem. Res., Vol. 47, No. 21, 2008 8507

In eq 16, Si is a source term due to coalescence and break-up. The right term of eq 17 describes all the forces acting on the phase k in each control volume: pressure gradient, gravity, the viscous stress term and the Fkm represents the sum of the interfacial forces that include the drag force FD, lift force FL, virtual mass force FVM, wall lubrication force FWL and turbulent dispersion force FTD. In the present hydrodynamic model, all forces (the drag force, lift force, and turbulence dispersion forces) except the virtual mass force have been used. According to Hunt et al.,31 the contribution of the virtual mass force becomes negligible for column diameters greater than 0.15 m. Thakre and Joshi;32 Deen et al.,9 and Sokolichin and Eigenberger33 have also pointed out the negligible effect of the virtual mass force. A brief description of each interfacial force component is presented below. The origin of the drag force is due to the resistance experienced by a body moving in the liquid. Viscous stress creates skin drag and pressure distribution around the moving body creates form drag. The later mechanism is due to inertia and becomes significant as the particle Reynolds number becomes larger. Re )

dS|ur| νL

(18)

In eq 18, ur ) uG - uL is the slip velocity, dS is the disperse phase Sauter mean bubble diameter,and νL is the liquid kinetic viscosity. The interphase momentum transfer between gas and liquid due to drag force is given by 1 3 FD ) CD∈GFL |uG - uL|(uG - uL) 4 dS

(19)

Figure 1. Grid and sparger configurations used in numerical simulations [black regions in (B) denote 25 source points corresponding to sparger representation].

where CD is the drag coefficient taking into account the character of the flow around the bubble. Hydrodynamic interactions of the bubble with other particles also influence the drag coefficient. This phenomenon can be taken into account.34 The lift force arises from the interaction between bubble and shear stress in liquid. The lift force in terms of the slip velocity and the curl of the liquid phase velocity can be described as

N

FL ) CL∈GFL(uG - uL) × ∇ × uL

DC )

∑χ nn

ij i j

(14)

j)1

2.2. Flow Equations. The numerical simulations presented are based on the two-fluid model Eulerian-Eulerian approach. The Eulerian modeling framework is based on ensembleaveraged mass and momentum transport equations for each phase. Regarding the liquid phase (∈L) as continuum and the gaseous phase (bubbles) as disperse phase (∈G), these equations without mass transfer can be written as Continuity equation of the liquid phase ∂ (F ∈ ) + ∇ (FL∈LuL) ) 0 ∂t L L

(15)

Continuity equation of the gas phase ∂ (F ∈ f ) + ∇ (FG∈GuGfi) ) Si ∂t G G i

(16)

Momentum equation ∂ (F ∈ u ) + ∇ (Fk∈kukuk) ) -∈k ∇ F + Fk∈kgi + ∇ [(µk + ∂1 k k k µk) ∇ (∈kuk)] + Fkm(k, m ) l, g)(17)

(20)

The sign of this force depends on the orientation of slip velocity with respect to the gravity vector. The CL has been modeled using Tomiyama et al.35 lift coefficient relationship. The lift coefficient can be expressed as

CL )

{

min[0.288tanh(0.121Re);f(Eod)]

Eod < 4

f(Eod) ) 0.00105Eod 4 e Eod e 10 0.0159Eod2 - 0.0204Eod + 0.474 Eod > 10 -0.29 (21) 3

where Re is the local Reynolds number of the gas phase (bubbles) and Eod is the Eotvos number. The origin of the wall lubrication force is due to the fact that liquid flow rate between bubble and the wall is lower than between the bubble and outer flow. This results in a hydrodynamic pressure difference driving bubble away from the wall. This force density is approximated as FWL ) -∈GFL

[

]

(ur - (urnw)nw) dS max Cw1 + Cw2 , 0 dS yw

(22)

where nw is the outward unit vector perpendicular to the wall


Figure 2. Comparison between the simulated and experimental profiles of gas hold-up at different axial positions in a 150 mm i.d. bubble column with multipoint sparger at VG ) 20 mm/s (A) H/D ) 0.2; (B) H/D ) 1.4; (C) H/D ) 2.6; (D) H/D ) 3.9; (E) H/D ) 5.0; and (F) H/D ) 6.2. ], Experimental; (1) RSM-PBM model; (2) RSM model with constant bubble size; and (3) k-ε model with constant bubble size.

and yw is the distance from the wall to the bubble. The wall lubrication constants Cw1 and Cw2 as suggested by Antal et al.36 are -0.01 and 0.05, respectively. The turbulent dispersion force, derived by Lopez de Bertodano,37 is based on the analogy with molecular movement. It approximates a turbulent diffusion of the bubbles by the liquid eddies. It is formulated as FTD ) -CTDFLk ∇ ∈L

(23)

where k is the turbulent kinetic energy per unit of mass. The turbulent dispersion coefficient of CTD ) 0.5 was found to give the good results which is in the recommended range of 0.1-1.0. By definition, the interfacial area aij for bubbly flows can be determined through the relationship aij )

6∈G dS

(24)

where dS is the Sauter mean bubble diameter. The local Sauter mean bubble diameter based on the calculated values of the scalar fraction fi and discrete bubble sizes di can be deduced from dS )

1

∑

fi id i

(25)

From the drag and nondrag forces above, it is evident that the interfacial area aij as well as the Sauter mean bubble diameter

in eq 24 are essential parameters that link the interaction between the liquid and gas (bubbly) phases. In most two-phase flow studies, the common approach of prescribing constant bubble sizes through the Sauter mean bubble diameter is still prevalent. Such an approach does not allow dynamic representation of the changes in the interfacial structure. 2.3. Turbulence Equations. For turbulence modeling, k-ε model, Reynolds stress model (RSM) have been used. In the RSM model, individual Reynolds stresses u′iu′j are computed via a differential transport equation. Thus, the RSM model solves six Reynolds stress transport equations. Along with these, an equation for dissipation rate is also solved. The exact form of Reynolds stress transport equations is derived by taking moments of exact momentum equation. This is a process wherein the exact momentum equations are multiplied by a fluctuating property, the production then being Reynolds averaged. The exact transport equations for the transport of Reynolds stresses Fu′iu′j are given by

(( )

∂ ∂ ∂ (F u′ u′ ) + (F u u′ u′ ) ) Pij + φij + µ+ ∂t L i j ∂xk L k i j ∂xk

)

2 k2 ∂u′iu′j 2 CSF - δijFLε(26) 3 ε ∂xk 3 where φij is the pressure-strain correlation, and P, the exact production term, is given by P ) -F(u′iu′j(∇u)T + (∇u)u′iu′j)

(27)


As the turbulence dissipation appears in the individual stress equations, an equation for ε is computed with the model transport equation:

(( ) ) ( )

µt ∂ε ∂ ∂ ∂ (F ∈ ε) + (FL∈Luiε) ) ∈ µ+ + ∂t L L ∂xi ∂xj L σε ∂xj ∂ui ε ε2 - Cε2FL∈L (28) ∈LCε1FL u′iu′k ∂xk k k When the k-ε model is used, the turbulent viscosity is formulated as follows µt ) FLCµ

k2 ε

µt 2 F βg ∇ T G ) µt ∇ u(∇u + ∇ uT) - ∇ u(3µt ∇ u + Flk) 3 Fσp l (32) For the continuous liquid phase, a k-ε model is applied with its standard constants: Cε1 ) 1.44, Cε2 ) 1.92, Cµ ) 0.09, σk ) 1.0, σε ) 1.3. No turbulence model is applied on the dispersed gas phase but the influence of the dispersed phase on the turbulence of the continuous phase is taken into account with Sato’s additional term.38 3. Numerical Solution

(29)

The turbulent kinetic energy k and its energy dissipation rate ε are calculated from their conservation equations

(( ) )

µt ∂k ∂ ∂ ∂ (FL∈Lk) + (FL∈LuLk) ) ∈L µ + + ∈L(G ∂t ∂xi ∂xi σk ∂xi ∈Lε)(30)

(( ) )

µt ∂ε ∂ ∂ ∂ (FL∈Lε) + (FL∈LuLε) ) ∈L µ + + ∂t ∂xi ∂xi σε ∂xi ε ∈L (Cε1G - Cε2∈Lε)(31) k where G is the turbulence production due to viscous and buoyancy forces, which is modeled using

The simulations were carried out as 3D transient flow pattern in a cylindrical bubble column using the commercial software ANSYS-CFX-10.0. The diameter of the column is 150 mm and the clear liquid height in the cylindrical column was 900 mm. Water was considered as the continuous phase, and air was considered as the dispersed phase. In the present study, bubbles ranging from 1m to 30 mm diameter are equally divided into 15 classes. The MUltiple Size Group (MUSIG) model has been used in CFX-10.0 to account for the nonuniform bubble size distribution in a gas-liquid mixture. The MUSIG model has been used extensively for different systems.15,17,39-42 The discrete bubble sizes prescribed in the dispersed phase were tracked by solving an additional set of 15 transport equations, which these equations were progressively coupled with the flow equations during the simulations. Instead of considering 16

Figure 3. Comparison between the simulated and experimental profiles of Sauter mean bubble diameter at different axial positions in a 150 mm i.d. bubble column with multipoint sparger at VG ) 20 mm/s (A) H/D ) 0.2; (B) H/D ) 1.4; (C) H/D ) 2.6; (D) H/D ) 3.9; (E) H/D ) 5.0; and (F) H/D ) 6.2. ], Experimental; (1) RSM-PBM model; (2) RSM model with constant bubble size; and (3) k-ε model with constant bubble size.


different complete phases, it was assumed that each bubble class travels at the same mean algebraic velocity to reduce the computational time and resource. This therefore results in 15 continuity equations for the gas phase coupled with a single continuity equation for the liquid phase. Sensitivity study on the number of size groups was performed through the consid-

Figure 4. Bubble classes holdup-based probability distribution. ], Experimental; (1) H/D ) 0.2; (2) H/D ) 1.4; (3) H/D ) 2.6; (4) H/D ) 3.9; (5) H/D ) 5.0; and (6) H/D ) 6.2.

eration of equally dividing the bubble diameters into 10, 15, and 20 size groups. The analysis revealed that no appreciable difference was found for the predicted maximum Sauter bubble mean diameter between the 15 and 20 bubble size groups. For the subdivision into 10 size groups, the maximum Sauter bubble diameter was under predicted by a maximum difference of 2%. In view of computational resources and times, it was therefore concluded that the subdivision of the bubbles sizes into 15 size groups were sufficient and hereafter in the following the computational results are all based on the discretization of 15 bubble size groups. In addition, simulations were carried out with constant bubble size of 6 mm using k-ε and RSM models. A solution to the two sets of governing equations for the balances of mass and momentum of each phase was sought. The conservation equations were discretized using the control volume technique. The velocity-pressure linkage was handled through the SIMPLE procedure. The hybrid-upwind discretization scheme was used for the convective terms. The computational grid is based on the unstructured set of blocks each containing structured grid. The structured grid within each block is generated using general curvilinear coordinates ensuring accurate representation of the flow boundaries. In order to select an adequate grid resolution, the effect of changing grid size was investigated. Several simulations were carried out using

Figure 5. Comparison between the simulated and experimental profiles of axial liquid velocity at different axial positions in a 150 mm i.d. bubble column with multipoint sparger at VG ) 20 mm/s: (A) H/D ) 0.2; (B) H/D ) 1.4; (C) H/D ) 2.6; (D) H/D ) 3.9; (E) H/D ) 5.0; and (F) H/D ) 6.2. ], Experimental; (1) RSM-PBM model; (2) RSM model with constant bubble size; and (3) k-ε model with constant bubble size.


Figure 6. Comparison between the simulated and experimental profiles of turbulent kinetic energy at different axial positions in a 150 mm i.d. bubble column with multipoint sparger at VG ) 20 mm/s (A) H/D ) 0.2; (B) H/D ) 1.4; (C) H/D ) 2.6; (D) H/D ) 3.9; (E) H/D ) 5.0; and (F) H/D ) 6.2: ], Experimental; (1) RSM-PBM model; (2) RSM model with constant bubble size; and (3) k-ε model with constant bubble size.

three different grid sizes; fine (454 330), medium (340 880), and coarse (133 570). It can be observed that result changes with insignificant difference if we increase the grid points by 20%. Further, it can be observed that there is practically no change in the gas volume fraction and axial liquid velocity profiles when the grid size increased beyond 340 880. In view of the observed effect of grid size, the simulations have been carried out by using 454 330 grid points. The time step 0.005 s was used. The default convergence criterion of CFX is that scaled residual of all equations fall below 1 × 10-4. Typical grid distribution and sparger representation is shown in A and B in Figure 1. Initially, there is no gas and no flow in the column (εL ) 1, Vk ) 0). At the column inlet (z ) 0), the sparger was incorporated by creating 25 source points at the specified position to mimic the exact sparger, at each source point, the mass flow rate based on the superficial gas velocity (0.02m/s) and nonuniform bubble size in the range from 1 mm to 30 mm diameter are equally divided into 15 classes and each class bubble volume fraction of 6.66% was specified. The turbulent kinetic energy (k ) (3/2)I2V2) and dissipation rates (ε ) (k3/2/ 0.3dh)) are low but nonzero. At the column outlet, the

atmospheric pressure was specified as boundary condition. No slip boundary conditions were used at all the impermeable walls. The simulations are carried out at superficial gas velocity of 20 mm/s and compared after getting steady state results. The fluids data are taken at room temperature (25 °C) and are treated as isothermal and saturated. Therefore, heat and mass transfer is neglected in the simulations. 4. Results and Discussion The CFD simulations were carried out using three different models [k-ε, RSM model with constant bubble size of 6 mm and RSM-PBM model with bubble size distribution (1-30 mm)] for the experimental conditions of Kulkarni et al.43,44 Kulkarni43 has carried out the experiments in a cylindrical plexiglass bubble column of 150 mm i.d. with a multipoint sparger having 0.25% free area (d0 ) 1.96 mm). Simultaneous measurements of axialradial and axial-tangential velocity measurements were made [laser Doppler anemometer (LDA)] at five axial levels (H/D ) 1.4, 2.6, 3.2, 3.9, and 5.0 away from the sparger). In addition, measurements of the axial, radial and tangential components


Figure 7. Comparison between the simulated and experimental profiles of turbulent energy dissipation rate at different axial positions in a 150 mm i.d. bubble column with multipoint sparger at VG ) 20 mm/s (A) H/D ) 0.2; (B) H/D ) 1.4; (C) H/D ) 2.6; (D) H/D ) 3.9; (E) H/D ) 5.0; and (F) H/D ) 6.2. ], Experimental; (1) RSM-PBM model; (2) RSM model with constant bubble size; and (3) k-ε model with constant bubble size.

were also made in the disengagement zone (H/D ) 6.2) and also in the sparger zone (H/D ) 0.2). In the radial direction, the measurements were made every 2 mm from the center and at the axial locations. Further, the measurements were made after a steady flow pattern was achieved (more precisely, 300 s after the beginning of sparging). The radial profiles of the fractional gas hold-up, Sauter bubble mean diameter, liquid velocity, turbulent kinetic energy, turbulent energy dissipation rate,and Reynolds stresses are compared with the predictions using the two-fluid and population balance model. The comparison was made at six different axial locations of H/D ) 0.2, 1.4, 2.6, 3.9, 5.0, and 6.2, in which the H/D ) 0.2 is close to the sparger region where the flow develops, the H/D ) 2.4 is nearly developing region, H/D ) 3.9 and 5.0 locations are developed regions, and H/D ) 6.2 is the disengagement region of the bubbly flow pattern. The following sections describe the simulation results and the experimental observation. 4.1. Fractional Gas Hold-up. In a bubble column reactor, the hold-up profile has different nature from sparger and the liquid phase circulation velocity also changes. In fact, the distance over which the hold-up profile (and liquid velocity profile) undergoes a change is called a sparger region. In this region as we go away from the sparger, the more number of bubbles are brought toward the center and less toward the wall,

which causes an increase in the steepness in the holdup profile. In the sparger region (H/D ) 0.2), the bubbling occurs continuously and the profile does not follow any trend. At higher axial levels, where the hold-up profile is developed, it takes a parabolic shape. Thus, with an increase in the hold-up in the center of column, the centerline mean axial liquid velocity also increases. This can be observed in Figure 2, which shows the comparison of the predicted fraction gas hold-up profiles with three different models (k-ε, RSM models with constant bubble size and RSM-PBM model with bubble size distribution) and the experimental data obtained from laser Doppler anemometer43 at different aspect ratios for superficial gas velocity of 20 mm/ s. Figure 2 shows that in the near sparger zone, the profile is relatively flat indicating small hold-up gradient from center to wall, whereas in the bulk region of the column, a parabolic profile can be observed, with maximum hold-up in the center and minimum at the wall. Also, the centerline hold-up was found to increase with an increase in the distance from the sparger. Further, the bubbling occurs continuously and there is a transient variation of bubbling area. In this region, all the turbulence (k-ε, RSM & RSM-PBM) models were not able to capture the effect of random variation of bubble formation and the flow pattern as reported experimentally. All the model predictions show several obscured peaks in the hold-up profile as shown in


Figure 8. Comparison between the simulated and experimental profiles of Reynolds stresses at different axial positions in a 150 mm i.d. bubble column with multipoint sparger at VG ) 20 mm/s (A) H/D ) 0.2; (B) H/D ) 1.4; (C) H/D ) 2.6; (D) H/D ) 3.9; (E) H/D ) 5.0; and (F) H/D ) 6.2: ], Experimental; (1) RSM-PBM model; and (2) RSM model with constant bubble size.

Figure 2A. This could be because of numerical instability and the use of finer grid could not change this situation. The predictive capability improves at higher axial locations (H/D > 2), where the flow gets developed, and stabilizes the oscillations. The prediction of hydrodynamics near the sparger requires further investigation. 4.2. Sauter Mean Bubble Diameter. The comparison of predicted and the experimental data of the local Sauter mean bubble diameter distribution is shown at different axial locations in Figure 3. It can be observed from Figure 3A that the model overestimates the bubble size near the sparger. It is well know that the gas distribution can play an important role on the evolution of the bubble size distribution in the gas-liquid reactors. Recently, Polli et al.45 is observed that, in the sparger region, increasing the distance from the gas distributor at constant superficial gas velocity, the bubble size distribution shifts from larger bubble size to smaller bubble sizes. This trend is the same that was observed far from the sparger by Colella et al.12 Further, Otake et al.,46 Colella et al.,12 and Polli et al.45 have found that the coalescence to be the dominant phenomenon only in the region immediately above the sparger. Good agreement was achieved against the measured bubble size for all the above the sparger locations with RSM-PBM model. It

can be observed that the local mean bubble size in the core region is larger than in the wall region. In the core region, gas holdup is higher, whereas the dissipation rate is lower, which results in larger bubbles because of the higher coalescence rate and lower breakup rate. Further, it can be seen that the bubble size distribution is almost uniform in the column cross section except near the wall region. The mean bubble size tends to reduce close to the column wall. This can be attributed to the strong velocity gradient, which causes further break-up of the bubbles into smaller sizes. Figure 4 shows the effect of height to diameter ratio on bubble class holdup based probability function in bubble column reactor using RSM-PBM model. It can be observed from the figure that the bubble class distribution gets stabilized after the H/D ) 3.0. It may be emphasized that the results reported in Figures 4 and 5 are possible by incorporating the population balance model. The conventional k-ε as well as the RSM models are incapable of predicting either the interfacial area or the bubble size distribution. 4.3. Axial Liquid Velocity. Figure 5 shows the comparison of predicted and experimental data of axial liquid velocity profiles at different aspect ratios using three different models. It can be observed from the figure that the liquid had a central


upward movement and downward near the wall. In the bulk region, all the models perform well and show good agreement with the experimental data. The center-line liquid velocity was seen to vary with the axial level. This is because, in a bubble column reactor, the hold-up profile has different nature from sparger to the disengagement zone and the liquid phase circulation velocity also changes. This development of the holdup profile (Figure 2) in coalescing liquid viz. tap water mainly occurs because of the axially varying reduced pressure gradient, which increases the bubble free area toward the disengagement zone. As a result, as we go away from the sparger, more of bubbles are brought toward the center and less toward the wall, which causes an increase in the steepness of the hold-up profile. Further, with an increase in the steepness, because of noticeable density gradients, liquid circulation is induced, which is upward in the region where the gas hold-up is high. In the sparger region (H/D ) 0.2), the bubbling occurs continuously and the profile does not follow any trend (Figure 5A). At higher axial levels, where the hold-up profile is developed, it takes a parabolic shape. Thus, with an increase in the hold-up in the center of column, the center-line mean axial liquid velocity also increases. 4.4. Turbulent Kinetic Energy. Figure 6 shows the comparison of different model predictions and experimental data of turbulent kinetic energy. It can be observed that the kinetic energy profiles typically exhibit a maximum around the crossover point for the time-averaged liquid axial velocity, due to large gradients and large fluctuations in the liquid velocity. It can be seen from the figure that the RSM-PBM model shows good agreement with the experimental data at all the locations. The k-ε and RSM (with constant bubble size) models over predict in the near sparger and bulk regimes. Further, the turbulent kinetic energy was found to decrease toward the wall. Also, the turbulent kinetic energy was found to increase away from the sparger. 4.5. Turbulent Energy Dissipation Rate. Figure 7 shows a comparison between the experimentally measured turbulent energy dissipation rate and the corresponding CFD predictions using different models. It can be seen that the turbulent energy dissipation rate increases very nominally from center to r/R ) 0.8 and then increases significantly close to the wall. The predictions were seen to deviate significantly in the near sparger region. The reason behind the higher values of turbulent energy dissipation rate from the measurements is discussed.43 Away from the sparger, the profiles obtained from the measurements were found to show a hump between r/R ) 0.4-0.6, following a decrease until r/R ) 0.7 and a sharp increase toward the wall. The predictions did not show any hump as seen in the experiments, however the increase in energy dissipation rate toward the wall was very well captured. The predicted values of energy dissipation rate showed consistent increase in values across the column cross-section. Further, it can be seen that the RSM-PBM model shows close agreement with the experimental data and the RSM, k-ε models with constant bubble size show under predictions. The predicted values of energy dissipation rate showed consistent increase in ε values across the column cross-section. 4.6. Reynolds Stresses. The numerical (RSM-PBM & RSM with constant bubble size) predictions of Reynolds shear stress, u′V′ were compared with the experimental data in Figure 8. It can be observed that the stress value continuously increases from center and then decreased rapidly well below zero just before the wall. In the literature,11,47-49 it has been shown that the point, at which stress profile attains a negative corresponds to the location where the highest downward axial liquid velocity

prevails. The RSM-PBM model prediction shows close agreement with the experimental data and the RSM model with constant bubble size gives over predictions. 5. General Remarks Recently, we have published50 on the same subject entitled CFD simulation of bubble columns incorporating population balance modeling. The paper uses in house code51-53 and the k-ε model accounts for size specific gas velocities for the bubbles, incorporates dispersion mechanism in the continuity equation and does not include dispersion force in the momentum equation. Furthermore, the paper uses the original break-up model of Luo and Svendson24 and modified coalescence model of Price and Blanch.25 Majority of these provisions could be made because the code was inhouse. The present paper uses the commercial code (CFX 10.0), and hence size-specific bubble velocities could not be given. In addition, the paper needs to include dispersion in the momentum equation and not in the continuity equation. The paper also uses the break-up and coalescence models in their original forms.24,25 In both the papers, the same experimental data are used and they show good agreement between the experimental data and the CFD predictions. For such an agreement, we feel that the following points are responsible: (i) in the present paper, which is not using the size-specific bubble velocities, the effect is compensated by the high coalescence rate; (ii) though an identical break-up model is used in both the cases, the bubble size boundary conditions are different. In the present paper, which uses commercial code, uniform size distribution can not be given, similar to that given by Bhole et al.50 The foregoing discussion brings out an urgent need for development of comprehensive break-up and coalescence models with the understanding of respective physics. Substantial additional work is needed for the description of mechanism of bubble formation in a situation like a bubble column. Most of the bubble formation models are available for a single orifice and a stagnant pool of liquid. The future work also needs to address the formulation of drag and lift forces for bubbles in a real situation of bubble column. 6. Conclusions A two-fluid model (Reynolds stress model) coupled with population balance approach is presented in this paper to handle gas-liquid flows in bubble column reactor. To demonstrate the population balance approach, the Multiple size group (MUSIG) model has been used in CFX-10.0 to account for the nonuniform bubble size distribution in a gas-liquid mixture. Population balance combined with coalescence and break-up models were taken into consideration. A detailed comparison has been presented between the CFD predictions and the local experimental data reported by Kulkarni et al.43,44 The development of flow pattern were examined at six different height to diameter ratios of H/D ) 0.2, 1.4, 2.6, 3.9, 5.0, and 6.2, in which the H/D ) 0.2 is close to the sparger region where the flow develops, the H/D ) 2.6 is nearly developing region, H/D ) 3.9 and 5.0 locations are developed regions, and H/D ) 6.2 is the disengagement region of the bubbly flow patterns. The results were compared the local fractional gas hold-up, Sauter mean bubble diameter, axial liquid velocity, turbulent kinetic energy, turbulent energy dissipation rate and Reynolds stress for superficial gas velocity of 20 mm/s. Overall, the CFD prediction gives good agreement in bulk region, some deviation was seen in the near sparger and disengagement zone. However,


in the gas phase, given that the assumption was invoked where each bubble class traveled at the same mean algebraic velocity to reduce the computational time and resources, a weakness of the model was evidenced at the near sparger and disengagement zone. Research is currently ongoing to consider addition momentum equations or develop an algebraic slip model to account for bubble separation, to yield a more realistic prediction. Further, the k-ε and RSM model prediction were compared with the experimental data. An attempt has been made to incorporate the real sparger in the CFD model. Appendix Nomenclature aij ) interfacial area, 1/m BB ) the “birth” due to breakup of bubbles, defined in eq 6, 1/ m3 s BC ) the “birth” due to coalescence of bubbles, defined in eq 13 1/ m3 s cf ) the increase coefficient of surface area defined in eq 5 Cµ, Cε1, Cε2 ) k-ε model constants, dimensionless CD ) drag coefficient, dimensionless CL ) lift coefficient, dimensionless CTD ) the turbulent dispersion coefficient, dimensionless Cw1, Cw2 ) the wall lubrication constants, dimensionless D ) column diameter, m d ) bubble diameter, m d, di, dj ) diameters (corresponding to Vi and Vj) of the daughter bubbles, m dh ) maximum bubble horizontal dimension of bubble ) (dS + 0.163Eod0.757)1/3, m dij ) The equivalent diameter ) (2/di + 2/dj)-1, m dS ) Sauter bubble mean diameter, m DB ) the “death” due to breakup of bubbles, defined in eq 7, 1/ m3 s DC ) the “death” due to coalescence of bubbles, defined in eq 14, 1/ m3 s Eod ) Eotvos number, ) (FL - FG)dh2/σ fi ) volume fraction of bubbles of group i, dimensionless fBV ) breakage volume fraction, dimensionless FCB ) calibration factor, dimensionless FD ) drag force, defined in Eq. (19), N Fkm ) the sum of the interfacial forces, N FL ) lift force, defined in eq 20, N FTD ) turbulent dispersion force, defined in eq 23, N FVM ) virtual mass force, N FWL ) wall lubrication force, defined in Eq. (22), N g ) gravitational acceleration, m/s2 h0, hj ) parameters represent the film thickness when collision H ) clear liquid height, m I ) turbulence intensity k ) turbulent kinetic energy, m2/s2 ni ) the number density of size group i, m3 nj ) the number density of size group j, m3 nw ) the outward unit vector perpendicular to the wall p ) pressure, Pa P ) the exact production term, defined in eq 27 r ) radial coordinate, m Re ) Reynolds number, defined in eq 18, dimensionless Si ) source term due to coalescence and break-up, kg/m3 s t ) time, s tij ) is the time required for two bubbles to coalesce having diameter di and dj, s uG ) gas velocity, m/s ur ) slip velocity, m/s

ut ) turbulent velocity, defined Eq. (10), m/s u′, V′ ) fluctuation velocity, m/s Vi ) volume of a bubble of group i, m3 Vj ) break-up rate of bubbles of volume, m3 yw ) the distance from the wall to the bubble, m z ) axial coordinate, m Greek letters β ) measured constant in eq 4, dimensionless F ) density, kg/m3 FC ) density of the continuous phase, kg/m3 σ ) surface tension ε ) turbulent energy dissipation rate, m2/s3 νL ) liquid kinetic viscosity, m2/s µ ) viscosity, Pa s µt ) turbulent viscosity, defined in eq 29, Pa s ∈G ) gas volume fraction, dimensionless ∈L ) liquid volume fraction, dimensionless φij ) the pressure-strain correlation θijB ) buoyancy collision rate, defined in eq 11, 1/ m3 s θijLS ) laminar shear collision rate, 1/ m3 s θijT ) turbulence collision rate, defined in eq 9, 1/m3 s λ ) size of an eddy, m ζ ) the size ratio between an eddy and a particle in the inertial subrange () λ/dj) τij ) the contact time for two bubbles, s χ ) parameter defined in eq 8, 1/m3 s Ω ) break-up rate, 1/m3 s Subscripts G ) gas L ) liquid min ) minimum w ) wall β ) measured constant in eq 4, dimensionless

Literature Cited (1) Krishna, R.; Ellenberger, J.; Sie, S. T. Reactor development for conversion of natural gas to liquid fuels: a scale-up strategy relying on hydrodynamic analogies. Chem. Eng. Sci. 1996, 51, 2041. (2) Dudukovic, M. P.; Degaleesan, S.; Gupta, P.; Kumar, S. B. Fluid dynamics in churn-turbulent bubble columns: measurements and modeling. In ASME Fluid Engineering DiVision Summer Meeting, Vancouver, BC, June 22-26, 1997; American Society of Mechanical Engineers: New York, 1997. (3) Joshi, J. B. 2001. Computational flow modeling and design of bubble column reactors. Chem. Eng. Sci. 2001, 56, 5893. (4) Ekambara, K.; Joshi, J. B. CFD simulation of mixing and dispersion in bubble column reactors. Chem. Eng. Res. Des. 2003, 81, 987. (5) Lorenz, O.; Schumpe, A.; Ekambara, K.; Joshi, J. B. Liquid phase axial mixing in bubble columns operated at high pressure. Chem. Eng. Sci. 2005, 60, 3573. (6) Ranade, V. V. Modelling of turbulent flow in a bubble column reactor. Chem. Eng. Res. Des. 1997, 75, 14. (7) Jakobsen, H. A.; Sannaes, B. H.; Grevskott, S.; Svendsen, H. F. Modeling of vertical bubble driven flows. Ind. Eng. Chem. Res. 1997, 36, 4052. (8) Sokolichin, A.; Eigenberger, G. Applicability of the standard turbulence model to the dynamic simulation of bubble columns: Part-I. Detailed numerical simulations. Chem. Eng. Sci. 1999, 54, 2273. (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. (10) 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, 689. (11) Ekambara, K.; Dhotre, M. T.; Joshi, J. B. CFD simulations of bubble column reactors: 1D, 2D and 3D approach. Chem. Eng. Sci. 2005, 60, 6733. (12) Colella, D.; Vinci, D.; Bagatin, R.; Masi, M.; Bakr, E. A. A study on coalescence and breakage mechanisms in three different bubble columns. Chem. Eng. Sci. 1999, 54, 4767.

8516 Ind. Eng. Chem. Res., Vol. 47, No. 21, 2008 (13) Sha, Z.; Laari, A.; Turunen, I., Implementation of population balance into multiphase-model in CFD simulation for bubble column. In Proceedings of the 16th International Congress of Chemical and Process Engineering, Prague, Czech Republic, Aug 22-26; Czhek Society of Chemical Engineers: Praha, Czech Republic, 2004. (14) Sha, Z.; Laari, A.; Turunen, I. Multi-phase-multi-size group model for the inclusion of population balances into the CFD simulation of gasliquid bubbly flows. Chem. Eng. Technol. 2006, 29, 550. (15) Olmos, E.; Gentric, C.; Vial, C.; Wild, G.; Midoux, N. Numerical simulation of multiphase flow in bubble column reactors: Influence of bubble coalescence and breakup. Chem. Eng. Sci. 2001, 56, 6359. (16) Lehr, F.; Millies, M.; Mewes, D. Bubble-size distributions and flow fields in bubble columns. AIChE J. 2002, 48 (11), 2426–2443. (17) Chen, P.; Sanyal, J.; Dudukovic, M. P. Numerical simulation of bubble columns flows: effect of different breakup and coalescence closures. Chem. Eng. Sci. 2005, 60, 1085. (18) Wang, T. F.; Wang, J. F.; Jin, Y. A CFD-PBM coupled model for gas-liquid flows. AIChE J. 2006, 52 (1), 125. (19) Jakobsen, H. A.; Lindborg, H.; Dorao, C. A. Modeling of bubble column reactors: progress and limitations. Ind. Eng. Chem. Res. 2005, 44, 5107. (20) Dhotre, M. T.; Niceno, B.; Smith, B. L. Large eddy simulation of a bubble column using dynamic sub-grid scale model. Chem. Eng. J. 2008, 136, 337. (21) Hu, G.; Celik, I. Eulerian-Lagrangian based large-eddy simulation of partially aerated flat bubble column. Chem. Eng. Sci. 2008, 63, 253. (22) Tryggvason, G.; Esmaeeli, A.; Lu, J.; Biswas, S. Direct numerical simulations of gas/liquid multiphase flows. Fluid Dyn. Res. 2006, 38, 660. (23) Lu, J.; Biswas, S.; Tryggvason, G. A DNS study of laminar bubbly flows in a vertical channel. Int. J. Multiphase Flow. 2006, 32, 643. (24) Luo, H.; Svendsen, H. Theoretical model for drop and bubble breakup in turbulent dispersions. AIChE J. 1996, 42, 1225. (25) Prince, M. J.; Blanch, H. W. Bubble coalescence and break-up in air sparged bubble columns. AIChE J. 1990, 36, 1485. (26) Chesters, A. K.; Hoffman, G. Bubble coalescence in pure liquids. Appl. Sci. Res. 1982, 38, 353. (27) Kirkpatrick, R. D.; Lockett, M. J. , The influence of approach velocity on bubble coalescence. Chem. Eng. Sci. 1974, 29, 2363. (28) Kim, W. K.; Lee, K. L. Coalescence behavior of two bubbles in stagnant liquids. J. Chem. Eng. Jpn. 1987, 20, 449. (29) Pohorecki, R.; Moniuk, W.; Bielski, P.; Zdrojkowski, A. Modelling of the coalescence/redispersion processes in bubble columns. Chem. Eng. Sci. 2001, 56 (21-22), 6157. (30) Rotta, J. C. Turbulente Stromungen; B. G. Teubner: Stuttgart, Germany, 1974. (31) Hunt, J. C. R. In Proceedings of the International Seminar on Transient Phenomena in Multiphase Flow, Dubrovnikm, Yugoslavia, May 25-29, 1987; International Centre for Heat and Mass Transfer: Ankara, Turkey, 1987. (32) Thakre, S. S.; Joshi, J. B. CFD simulation of bubble column reactors: Importance of drag force formulation. Chem. Eng. Sci. 1999, 54, 5055. (33) Sokolichin, A.; Eigenberger, G. Simulation of buoyancy driven bubbly flow: Established simplifications and open questions. AIChE J. 2004, 50, 24.

(34) Ishii, M.; Zuber, N. Drag coefficient and relative velocity in bubbly, droplet or particulate flows. AIChE J. 1979, 25, 843. (35) Tomiyama, A.; Tamai, H.; Zun, I.; Hosokawa, S. Transverse migration of single bubbles in simple shear flows. Chem. Eng. Sci. 2002, 57, 1849. (36) Antal, S.; Lahey, R.; Flaherty, J. Analysis of phase distribution in fully-developed laminar bubbly two phase flow. Int. J. Multiphase Flow 1991, 7, 635. (37) Lopez de Bertodano, M. A. Turbulent bubbly two-phase flow in a triangular duct. Ph. D. dissertation, Rensselaer Polytechnic Institute, Troy, NY, 1992. (38) Sato, Y.; Sekoguchi, K. Liquid velocity distribution in two-phase bubbly flow. Int. J. Multiphase Flow 1975, 2, 79. (39) Lo, S. Application of Population Balance to CFD Modelling of Bubbly Flow Via the MUSIG model; Technical Report AEAT-1096; AEA Technology: Oxford, U.K., 1996. (40) Yeoh, G. H.; Tu, J. Y. Population balance modelling for bubbly flows with heat and mass transfer. Chem. Eng. Sci. 2004, 59, 3125. (41) Yeoh, G. H.; Tu, J. Y. A unified model considering force balances for departing vapour bubbles and population balance in sub-cooled boiling flow. Nucl. Eng. Des. 2005, 59, 3125. (42) Cheung, S. C. P.; Yeoh, G. H.; Tu, J. Y. On the modeling of population balance in isothermal vertical bubbly flows-Average bubble number density approach. Chem. Eng. Process. 2007, 46 (8), 742. (43) Kulkarni, A. A.; Ekambara, K.; Joshi, J. B. On the development of flow pattern in a bubble column reactor: Experiments and CFD. Chem. Eng. Sci. 2007, 62, 1049. (44) Kulkarni, A. A.; Joshi, J. B.; Ramkrishna, D. Determination of bubble size distributions in bubble columns using LDA. AIChE J. 2004, 50, 3068. (45) Polli, M.; Di Stanislao, M.; Bagatin, R.; Abu Bakr, E.; Masi, M. Bubble size distribution in the sparger region of bubble Columns. Chem. Eng. Sci. 2002, 57, 197. (46) Otake, T.; Tone, S.; Nakao, K.; Mitsuhashi, Y. Coalescence and break-up of bubbles in liquids. Chem. Eng. Sci. 1977, 32, 377. (47) Menzel, T.; Weide, T.; Staudacher, O.; Onken, U. Reynolds stress model for bubble column reactor. Ind. Eng. Chem. Res. 1990, 29, 988. (48) Devanathan, N.; Moslemian, D.; Dudukovic, M. P. Flow mapping in bubble columns using CARPT. Chem. Eng. Sci. 1990, 45, 2285. (49) Pan, Y.; Dudukovic, M. P.; Chang, M. Dynamic simulation of bubble flow in bubble columns. Chem. Eng. Sci. 1999, 54, 2481. (50) Bhole, M. R.; Joshi, J. B.; Ramkrishna, D. CFD simulation of bubble columns incorporating population balance modeling. Chem. Eng. Sci. 2008, 63, 2267. (51) Ranade, V. V.; Mishra, V. P.; Saraph, V. S.; Deshpande, G. B.; Joshi, J. B. Comparison of axial flow impellers using LDA. Ind. Eng. Chem. Res. 1992, 31, 2370. (52) Patwardhan, A. W.; Joshi, J. B. Design of gas inducing reactors. Ind. Eng. Chem. Res. 1999, 37, 49. (53) Thakre, S. S.; Joshi, J. B. CFD simulation of flow in bubble column reactors: importance of drag force formulation. Chem. Eng. Sci. 1999, 54, 5055.

IE071393E

CFD Simulation of Bubble Column Reactor Using Population Balance

Recommend Documents