Ind. Eng. Chem. Res. 2005, 44, 4959-4967
4959
Prediction of Backmixing and Mass Transfer in Bubble Columns Using a Multifluid Model D. Wiemann and D. Mewes* Institute of Process Engineering, University of Hannover, Callinstrasse 36, 30167 Hannover, Germany
Three-dimensional, instationary flow fields in bubble columns are numerically calculated using a Euler multifluid model. The local bubble size is calculated depending on the flow field using a transport equation for the mean bubble volume originating from a population balance equation. The resulting set of equations enables the simultaneous calculation of mass transfer and backmixing in the gas and liquid phases. Therefore, a perfect tracer is added to both the gaseous and liquid phases. From the calculations, the time-dependent concentration fields of the tracer are obtained. Assuming the one-dimensional dispersion model, the axial dispersion coefficients are calculated for comparison with experimental results. The obtained dispersion coefficients for the gas and liquid phases are in good agreement with the experimental investigations of several authors. On the basis of the reasonable prediction of the interfacial area and backmixing, the multifluid model is extended to consider mass transfer. The absorption of the gas phase reduces the overall momentum transfer; thus, the flow pattern is changed. Introduction The flow pattern in bubble columns is induced by the rising-bubble swarm. For low superficial gas velocities, the homogeneous flow regime occurs. The bubbles are of almost equal size and rise up straightforward. For higher superficial gas velocities, the enhanced coalescence of smaller bubbles leads to the formation of larger bubbles, which rise up much faster. In this heterogeneous flow regime, the liquid flow field is characterized by large-scale vortices.1 Several research groups performed experimental investigations of the flow structure in bubble columns. The local phase distribution in cylindrical bubble columns is measured by Warsito et al.2,3 using ultrasound tomography. Devanathan et al.4 use a radioactive tracer particle to determine the flow structure in bubble columns. A summary over the experimental investigations and modeling aspects in bubbly flow is given by Dudukovic et al.,5 Dudukovic,6 and Joshi.7 The conventional design approach for bubble-column reactors is based on one-dimensional reactor models. Therefore, several empirical and semiempirical correlations have been proposed to determine the overall interfacial area density, the integral gas holdup, and mixing parameters. The applicability of these correlations is, in principle, limited by the experimental conditions for which they are obtained. Thus, the scaleup from the laboratory scale to the industrial scale is difficult. In particular, these models only give information about integrated quantities, whereas the local quantities, which govern the transport processes of energy, momentum, and mass, remain unresolved. With increasing computer power, computational fluid dynamics (CFD) is an alternative way to gain further insight into the complex two-phase flow. Several calculations of the large-scale flow field in bubble columns have been published. On the local scale, detailed simulations of the * To whom correspondence should be addressed. Tel.: +49 511 762 3828. Fax: +49 511 762 3031. E-mail: dms@ ifv.uni-hannover.de.
bubble motion8 lead to a further understanding of the bubble-bubble interaction and the influence of bubbles on the turbulence in the continuous phase. However, these calculations are restricted to a few bubbles and small dimensions. Therefore, these calculations are mainly used to investigate local phenomena in bubbly flow. Numerical calculations of the two-phase flow in two-dimensional bubble columns have been performed by Sokolichin and Eigenberger,9 Pan et al.,10,11 and Pfleger et al.,12 to name only a few of the many publications. Flow fields in three-dimensional bubble columns have been calculated based on Eulerian descriptions by Torvik and Svendsen,13 Grienberger,14 Jakobsen et al.,15 and Ranade16 under the assumption of an axisymmetric and stationary flow field. For the calculation of the instationary flow fields in threedimensional cylindrical bubble columns, Sokolichin et al.17 and Lain et al.18 used a Euler-Lagrange approach. For the same purpose, Krishna et al.19 and Sanyal et al.20 used a Euler-Euler approach. Lehr et al.21 also used a Euler-Euler model, but they also considered the breakup and coalescence of bubbles depending on the local flow field. The flow field is significantly influenced by the coalescence and breakup of bubbles. For the description of the bubble-size distribution, the population balance equation is widely used. On the basis of the population balance equation, Ishii and Kim22 proposed a transport equation for the interfacial area density in bubbly flow, which predicts the interfacial area density in good agreement with experimental results. The coupling of the population balance equation with CFD enables the calculation of the local interfacial area density. One approach is the discretization of the population balance equation by several classes, which represents a certain bubble diameter. However, the population balance equation has to be discretized by a large number of classes because a wide range of bubble diameters has to be considered. The multiple interactions between the different bubble classes limit this approach to only a small number of classes. Therefore, the need for an
10.1021/ie049163c CCC: $30.25 © 2005 American Chemical Society Published on Web 03/12/2005
4960
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005
Table 1. Transport Equations for the Mean Bubble Volume (Lehr et al.21) I
∂vj 1 µg vj 1 1 DFg +b ug1∇(vj 1) ) -Z1vj 1 + 0.3463r2(vj 1,vj 1) R1 vj + n˘ ap1 ∂t Fg Dt 1 Fg R1
II
∂vj 2 µg vj 2 1 DFg +b ug2∇(vj 2) ) -Z2vj 2 + 0.4250r2(vj 2,vj 2) R2 + 0.9024r2(vj 2,5vj 1) R1 vj 2 + n˘ ap2 ∂t Fg Dt F g R2
III
∂R1Fg R2R1 R12 ug1) ) Z2FgR2 - 0.9024r2(vj 2,5vj 1) Fg - 3.1043r2(vj 1,vj 1) Fg + n˘ ap1µg + ∇(R1Fgb ∂t vj 2 vj 2
IV
∂R2Fg R2R1 R12 + ∇(R2Fgb ug2) ) -Z2FgR2 + 0.9024r2(vj 2,5vj 1) Fg + 3.1043r2(vj 1,vj 1) Fg + n˘ ap2µg ∂t vj 2 vj 2 Z1 ) 0.6082 Z2 ) 0.2545
r2(vi,vj) )
() Ff σ
7/5
() σ Fl
[(
[
19/15vj 15/9 exp 1 vj 24/9
()
0.8565 σ Ff vj
x
1
1 3/5
9/10
]
2/5 1/15
π 6 v 4 π i
1/3
)
+
6 v π j
1/3
( )
]
2
[(
min(u′,ucrit) exp -
efficient coupling of population balances and CFD arises. In this paper, a transport equation for the mean bubble volume originating from a population balance equation approach is coupled with a Euler multifluid model. For the calculation of reactive flows, not only the interfacial area but also the mixing properties of the column has to be reasonably predicted. Therefore, a perfect tracer is injected into the liquid and the risingbubble swarm to determine backmixing. Governing Equations for Bubbly Two-Phase Flow For the modeling of bubbly flow, the Euler-Lagrange and Euler-Euler approaches are common. In addition, the direct numerical calculation can be used to calculate the motion of single bubbles and the movement of the gas-liquid interface. For large-scale flow fields, the Euler-Euler approach is most suitable because it is numerically more efficient than the Euler-Lagrange approach. The Euler approach assumes all phases to be continuously distributed in the computational domain. Therefore, also the dispersed gas phase is described by a quasi-continuous fluid. The interphase transfer of energy, mass, and momentum is considered via interphase transfer terms. The interfacial area density is obtained by solving the transport equations for the mean bubble volume of the large and small bubble fractions. This transport equation originates from a population balance equation for the bubble number density. The breakup and coalescence of bubbles is described by a physical-based model. From the numerical solution of the population balance equation, selfsimilar bubble-size distributions are obtained. For the heterogeneous flow regime, a bimodal bubble-size distribution is calculated; thus, the gas phase can be divided into one bubble class containing the small bubbles and a second class containing the large bubbles. From the numerical solution, the transport equations for the mean bubble volume of the small and large bubbles are obtained, as given in Table 1.21 The number of bubbles, the volume fraction, and the bubble size are obtained from the transport equations. Because both
Rmax1/3 R1/3
)] 2
-1
bubble fractions are coupled by interphase transfer terms for breakup and coalescence, the volume fractions of the small and large bubbles are not fixed. The transport equations are coupled with the balance equations for mass and momentum; thus, the local bubble size is calculated depending on the local flow field. The solution of these transport equations also enables the calculation of the Sauter diameter of the small and large bubble fractions. The momentum equation for the liquid phase, the small bubble phase, and the large bubble phase is written as
∂ u ) + ∇[Ri(Fib uib ui)] ) -Ri∇pi + ∇{Riηi[3u bi + (R F b ∂t i i i g+F Bij i ) l, 1, 2 (1) (3u bi)T]} + RiFib The terms on the left-hand side of eq 1 describe the local and convective change of momentum. These terms are balanced by the forces on the right-hand side, namely, the forces due to the bulk pressure gradient, molecular transport, gravitational force, and interphase transfer of momentum. The interphase momentum transfer is mostly influenced by the drag force. Other forces such as added mass force, lift force, or extra dispersion forces are not considered. In particular, the modeling of the dispersion forces and the lift force is doubtful because a rigorous modeling of these forces has not been arrived at yet, as discussed by Loth,23 Joshi,7 and Sokolichin et al.24 Therefore, the interphase momentum transfer is calculated as
3 Ri b -b ul|(u bi - b u l) F Bil ) CD Fl |u 4 di i
(2)
The drag coefficient CD is calculated following Clift et al.25 as
CD ) max
24 (1 + 0.1Re [Re
0.75
);
{ (
) }] (3)
2 8 min max 0.44, Eo1/2 , 3 3
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 4961
The Reynolds and Eo¨tvo¨s numbers are defined as follows:
ul|di |u bi - b νl
(4)
g(Fl - Fi)di2 σ
(5)
Rei )
Eo )
In eq 5, the surface tension between the liquid and gas phases is σ. The Sauter bubble diameter di is calculated from a transport equation for the mean bubble volume. Because of bubble coalescence and breakup, second-order fluxes arise in the momentum equation for the small and large bubble fractions. These fluxes are considered in addition to the terms given in eq 1. In the case of mass transfer, the additional term
{
-M ˙ ifjb ui, phase i F Bmass ) M ˙ ifjb ui, phase j
{
(6)
{
∂(RiζAiFi) -ζ M ˙ ifl, i ) 1, 2 ui) ) +ζAiM + ∇(RiζAiFib ∂t Aj ˙ jfl, i ) l; j ) 1, 2 (8) has to be solved in addition to eq 7. Again source and sink terms due to bubble coalescence and breakup are considered in addition to eq 8. In the liquid phase, four components are present: two identical components (which are the liquid counterparts of the gaseous components), a tracer component, and a background fluid. The calculations also consider mass transfer between the gaseous and liquid phases. Therefore, an additional transport equation for the transferring component is considered in both phases, accounting for convective and molecular transport. For the calculation of the mass transfer rate, the phase equilibrium at the gas-liquid interface is described following Henry’s law. The phase equilibrium is described with
( )
(10)
with the bulk molar concentration of the liquid phase cl and the bulk molar concentration of the transferred component cA,l. The mass-transfer coefficient is calculated from a Sherwood number depending on the bubble volume. For small Reynolds numbers, the bubbles are spherical-shaped, whereas with increasing Reynolds number, turbulent motions at the surface become more and more important and the bubbles lose their shape. Therefore, for small Reynolds numbers below
Recrit ) 3.73(Flσ3/gη4)0.209
(11)
the Sherwood number is calculated according to Brauer26
with
Sh∞ ) 2 +
0.651(Re × Sc)1.72 1 + (Re × Sc)1.22
(9)
In eq 9, the molar masses of the liquid and gas phases are labeled µl and µg. The Henry constant is H. The mass transfer rate is calculated as
(13)
For larger Reynolds numbers, the correlation
Sh ) 2 + 0.015Re0.89Sc0.7 (7)
considering mass transfer from the gaseous phase to the liquid phase. In addition, source terms due to bubble coalescence and breakup, as given in Table 1, are considered. The gas phase consists of two components that are identical in their properties. One of these components is used as the tracer. Therefore, the mass balance equation for one component
µl FA,pl ) pA,gFl/ H µg
cl β (F - FA,l) cl - cA,l l A,pl
Sh ) Sh∞[(1 + 0.433Re2)-1 + 0.0000423]-0.055 (12)
arises in the momentum equation considering secondorder momentum fluxes due to mass transfer from phase i to phase j. The mass balance equation is written as
∂(RiFi) -M ˙ ifl, i ) 1, 2 + ∇(RiFib ui) ) +M ˙ jfl, i ) l; j ) 1, 2 ∂t
m ˘ i,l )
(14)
is used instead of eq 12. In bubbly two-phase flow, the turbulent velocity fluctuations in the liquid phase are caused by the shearinduced turbulence as well as the presence of bubbles. Several experimental investigations deal with the influences of a dispersed phase on the continuous-phase turbulence. In this work, the influence of bubbles on the liquid turbulence is described following the proposal of Lopez de Bertodano et al.27 Therefore, in addition to the well-known transport equations for the turbulent energy dissipation rate and the turbulent kinetic energy, additional terms are calculated depending on the interphase drag. The resulting kinetic energy and dissipation rate are calculated as
1 1 b1|2 + R2|u b2 | 2 kres ) kl + kb ) kl + R1|u 4 4 F Bd,1,l F Bd,2,l res ) l + b ) l + |u b1 | + |u b2 | Fl Fl
(15)
The second and third terms on the right-hand side of eq 15 represent the bubble-induced dissipation rate and turbulent energy for the small and large bubbles. The gas-phase viscosity is described by the molecular viscosity of the gas. Equations 1-15 are solved with the method of finite volumes. Therefore, the code CFX-5.7 is used. Dispersion Models In bubble columns, a gaseous phase is dispersed in a continuous liquid; thus, the rising-bubble swarm induces a circulating flow field. The resulting flow field is three-dimensional and instationary. Because of the vortical flow field, bubble coalescence, and breakup, a residence time distribution arises in the liquid and gas phases. Thus, for the dimension of bubble columns, backmixing has to be taken into account. Several
4962
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005
experimental investigations of backmixing in bubble columns have been performed so far. For instance, backmixing in the gas phase is investigated by Shetty et al.,28 and liquid backmixing is investigated by Degaleesan et al.,29 to name only two of the various investigations that have been performed. Recently, Yang and Fan30 investigated the axial liquid mixing in pressurized bubble columns. In most of the experimental investigations, the deviation from an ideal plug-flow reactor is described with the axial dispersion model. Backmixing in the liquid and gas phases is caused by several mechanisms. The formation of liquid vortices, the occurrence of stagnation zones, and the turbulent velocity fluctuations cause a residence time distribution in the liquid phase. In addition, the coalescence and breakup of bubbles and the nonuniform distribution of the gas phase enhance the formation of a velocity profile in both phases. Changes in the bubble-size distribution and coalescence and breakup of bubbles also enhance backmixing in the gaseous phase. For the description of backmixing, several models have been developed, as summarized by Levenspiel31 and Fogler.32 One of the most popular ones for backmixing in bubble columns is the axial dispersion model, which describes the deviation from ideal plug flow. The mass flux due to dispersion is described in analogy to the molecular transport, although both phenomena differ in their physical nature. The resulting molar flux density
n˘ ) -(E + D)
∂c ∂x
(16)
is described depending on the concentration gradient. The dispersion coefficent is labeled E, and the molecular diffusion coefficient is D. The limiting cases for the dispersion coefficient are zero for an ideal plug-flow reactor and infinity in the case of an ideal stirred reactor. For bubbly flow, the molar flux due to dispersion is higher than that due to molecular transport; thus, the mass balance for a cross-sectional element of the bubble column can be written as
(
)
∂c ∂2 c ∂2 c ∂c ji ∂c + + E , i ) g, l + ) Ei,rad i,ax ∂t Ri ∂x ∂r2 r ∂r ∂x2 (17) neglecting molecular transport. The first term on the right-hand side of eq 17 describes the radial dispersion flux, whereas the second term considers the axial dispersion. For a given concentration field, the solution of eq 17 enables the calculation of the dispersion coefficients. The axial dispersion is obtained from
∂2c ∂c ji ∂c + ) Ei,ax 2, i ) g, l ∂t Ri ∂x ∂x
(18)
For a batch-operated liquid phase and a pulsewise tracer injection, a homogeneous distribution of the tracer is reached for long times. The liquid dispersion coefficient is calculated from eq 18 considering the following boundary and initial conditions:
∂c )0 ∂x c(x,o) ) co c(x,o) ) 0
x ) 0, L
(19)
for 0 e x e δ
(20)
for x g δ
(21)
Assuming a batch-operated liquid, the solution of eqs 19-21 is approximated by the following summation of infinite terms following Ohki and Inoue33 and Qi et al.: 34
c c∞
∞
)1+2
∑ n)1
{ ( ) [ ( ) ]} cos
nπ nπ 2 Eax,lt x exp H H
(22)
The stationary concentration of the tracer is labeled c∞, and the instantaneous concentration is c. For practical use of eq 22, only the first 10 terms have to be included. The left-hand side of eq 22 can be obtained by experiments or, as in this work, from a numerical calculation of the three-dimensional, instationary flow field. The dispersion coefficient is then calculated by varying the value of the dispersion coefficient on the right-hand side of eq 22; thus, the deviation between the quotient of concentrations and the infinity terms is minimized using the least-squares method. The tracer in the gas phase is injected continuously in the rising-bubble swarm. Therefore, a sinusoidal tracer inlet concentration is used. From the mathematical point of view, the application of this frequency response method requires a linear system response; thus, only the amplitude of the tracer signal is changed, whereas the frequency remains constant. The experimental investigations of Mangartz and Pilhofer35 confirm this assumption for bubble columns. The solution of eq 18 for the case of a sinusoidal input signal leads to the following equation for the axial dispersion coefficient:
Eax,g )
jg3 Rg3ω2H
( )
ln
cin cout
(23)
The angular frequency of the sinusoidal input signal is ω, the dispersion height is H, the integral volume fraction of the gas is Rg, and the superficial gas velocity is jg. The amplitudes of the tracer concentrations are cin at the inlet and cout at the outlet. Results Equations 1-15 are solved with the method of finite volumes. Therefore, the code CFX-5.7 is used. The computational domain is discretized with a blockstructured grid of hexahedral volumes. The edge length is about 1 cm. For the temporal discretization, time steps of 0.02 s are made. The liquid and gas phases enter the column through an inlet at the bottom. The boundary conditions are uniformly over the inflow area. At the top of the column, the gas phase flows out of the column, whereas the liquid phase flows out through an orifice in the circumferential direction. The convective terms are discretized with second-order accuracy. Calculations of the three-dimensional, time-dependent flow field are performed for cylindrical bubble columns considering the homogeneous and heterogeneous flow regimes. In Figure 1, the calculated instantaneous flow field in a bubble column with 3 m height and 40 cm diameter is shown. The superficial gas velocity is 0.1 m/s; thus, the heterogeneous flow regime is present. The liquid streamlines are colored with the specific interfacial area density, the volume fraction of the gas, and the axial liquid velocity.
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 4963
Figure 3. Liquid tracer distribution. Figure 1. Instantaneous flow field in a bubble column.
Figure 4. Calculated tracer concentration in the liquid.
Figure 2. Dispersion of ink in a laboratory-scale bubble column.
The flow field is characterized by several large-scale vortices in the order of the column diameter. Large values of the interfacial area density are calculated, in particular, in the core region of the column, whereas a low interfacial area density is calculated near the walls. The axial liquid velocity reaches from -1.0 to +1.0 m/s. In correspondence with the regions of high gas volume fractions, the liquid is transported upward mainly in the core region and flows downward near the wall of the column. Backmixing in the Liquid Phase In bubbly two-phase flow, backmixing leads to deviations from plug flow. The principle influence of the vertical flow structure on the dispersion of a tracer is visualized by the injection of ink into water. Therefore, a laboratory-scale bubble column with an inner diameter of 15 cm is operated at a low superficial gas velocity. The tracer substance is blue ink. In Figure 2, the timedependent spreading of the tracer is shown for several time steps. First the ink is transported downward, in particular near the column walls. With increasing time, the ink is distributed over the whole cross section of the column and then transported downward. For the numerical calculation of the tracer dispersion, a similar arrangement is used. The liquid tracer is injected pulsewise at the top of the column, and the time-dependent dispersion of the tracer is calculated. For a bubble column of 2 m height and 20 cm diameter, the view from the outside onto the column is depicted in Figure 3 for several time steps. Immediately after the injection, the tracer is transported downward near the column wall. For larger
times, the tracer is distributed over the whole cross section of the column because of the vortical flow structure. Finally, an almost homogeneous distribution of the tracer is reached. These results agree, in principle, with the tracer dispersion as shown in Figure 2. The time-dependent concentration for several points at two axial positions is shown in Figure 4. In addition, the calculated concentrations are given for different radial positions. The tracer concentrations first increase at the upper axial positions 5s after the injection. The tracer concentrations reach maximum values and then decrease to constant values for longer times. In accordance with their longer distance to the injection point, the concentrations at the lower position increase later. The concentrations increase slightly and reach the constant value 68 s after the injection. Depending on the radial position, the fluctuations of the tracer concentration vary. In particular, at the upper position, fluctuations increase from the center position toward the column wall. With increasing distance from the injection point and increasing time, the fluctuations are damped. For comparison with experimental results, the axial dispersion model is applied to calculate the axial dispersion coefficients. Therefore, the tracer response curves, as depicted in Figure 1, are approximated using eq 22. The method of least squares is applied to determine the axial dispersion coefficient. For the dimensionsless representation, the Bodenstein and Froude numbers
Bo ) jgD/Eax,l
(24)
Fr ) jg2/gD
(25)
are used. The calculated Bodenstein and Froude numbers are shown in Figure 5. In addition, the experimental results of Deckwer are given.
4964
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005
Figure 5. Calculated Bodenstein number depending on the Froude number.
Figure 7. Time-dependent concentration of the gas tracer.
Figure 6. Calculated response curves for several superficial liquid velocities.
With increasing superficial gas velocity, the dispersion coefficients increase; thus, enhanced backmixing in the liquid phase occurs. For several applications, the liquid is fed continuously into the reactor rather than operated batchwise. Therefore, a superficial liquid velocity is considered in the calculations. The tracer is then injected pulsewise at the liquid inlet at the bottom of the column. The threedimensional, time-dependent concentration field is calculated, and the response curves are obtained for several positions in the column. In Figure 6, the response curves at the top of the column are shown depending on the time for several superficial liquid velocities. The column height is 1.4 m and the diameter 0.15 m. The time delay after the injection increases with a decrease in the superficial liquid velocity. All mass fraction curves reach a maximum value and then decrease. For low superficial velocities, the maximum is shifted toward longer times and the maximum value is reduced. Thus, for a high superficial liquid velocity, the residence time is reduced and a narrow residence time distribution is obtained. For the gas phase, the tracer is injected continuously into the rising-bubble swarm. Therefore, the concentration of the tracer component at the inlet is varied with time using a sinusoidal profile. The initial signal is then disturbed because of backmixing in the gas phase. The calculated tracer concentration in the gas phase is shown in Figure 7 for several time steps. The tracer concentration varies in the axial and radial directions. In the case of plug flow, a sinusoidal tracer concentration with the initial amplitude would be obtained without a gradient in the radial direction. In particular, in the upper part of the column, the influence of backmixing is obvious. In Figure 8, typical response
Figure 8. Calculated gas tracer response.
Figure 9. Axial dispersion coefficient for the gas phase.
curves for the gaseous phase are shown for points along the column axis. The amplitude of the tracer concentration decreases with an increase in the axial position. However, the frequency remains constant; thus, the application of the frequency response method is reasonable. Using eq 23, the axial dispersion coefficients for the gas phase are obtained. Calculations are performed for a cylindrical bubble column in accordance with the experimental investigations of Mangartz and Pilhofer.35 In Figure 9, the calculated and experimental dispersion coefficients are shown. The dispersion coefficients increase with an increase in the superficial gas velocity. For low superficial gas velocities corresponding to the homogeneous flow regime, the increase is slighter than that in the heterogeneous flow regime. Thus, the enhanced formation of vortices and the formation of a bimodal bubble-size
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 4965
Figure 10. Time-dependent volume fraction of gas for the physical absorption.
distribution significantly influence the mixing properties of the bubble column. The proposed multifluid model reasonably predicts the large-scale mixing behavior of the bubble columns. The multifluid model is then extended to account for mass transfer between the gas and liquid phases. The physical absorption of the gas phase is considered for the system of carbon dioxide and water. The volume fraction of gas considering the physical absorption of carbon dioxide into water is shown in Figure 10 for several time steps. In addition to the cross-sectional areas, a surface with a constant volume fraction of 0.05 is given. The absorption process reduces the volume fraction of gas along the column height; thus, the momentum transfer is reduced also. As a result, the flow pattern changes with a decrease in the volume fraction. In the lower part of the column, liquid vortices arise, whereas in the upper part of the column, the amount of vortices is reduced. The absorption influences the bubble diameter as well. In Figure 11, the mass fraction of carbon dioxide in the liquid, the bubble diameter, and the liquid velocity are shown for the same column. With increasing column height, the mass fraction of carbon dioxide in the liquid increases. Because of absorption, the bubble diameter decreases, in particular in the upper part of the column. In contrast, for nonreactive bubbly flow, an increase has to be expected. The time-averaged flow field is axisymmetric similar to nonreactive flow. The decrease in the gas volume fraction occurs along the column height and in the radial
direction. For a cross section at h ) 1.6 m, the calculated time-averaged profile of the gas volume fraction and experimental results from Hillmer36 are given in Figure 12. The absorption process influences the flow pattern in bubble columns. Therefore, the flow fields are calculated for several mass-transfer rates. The calculated flow fields for bubble columns with 2 m height and 20 cm diameter are shown in Figure 13 with and without considering mass transfer. The calculations consider bubbly flow without absorption, the physical absorption, and the enhanced absorption of the gas phase. The absorption of gas reduces the volume fraction of gas with an increase in the column height. The overall gas holdup for the cases A-C are 0.318, 0.278, and 0.085. In particular, for the case of an enhanced absorption, the volume fraction is significantly reduced. Thus, the interphase momentum transfer is reduced according to the reduced gas volume fraction, and the flow pattern in the bubble column is changed. With increasing masstransfer rates, the amount of vortices is reduced. For the enhanced absorption, only one large circulation loop arises. Because of the absorption, the bubble diameter decreases. The enhanced absorption significantly changes the bubble diameter, as shown in Figure 14. The bubble diameters of the large and small bubble fractions are given for the case without mass transfer and for the enhanced absorption. In particular, in the downward-flowing part of the circulation loop, small bubbles are calculated in correspondence with the low gas volume fraction. Conclusion A multifluid model is applied to calculate threedimensional, instationary flow fields in bubble columns considering the local bubble size. Therefore, the gas phase is decribed by two interpenetrating phases, which represent the small and large bubble fractions. Considering multiple components in the gas and liquid phases, the dispersion of a tracer in both the gas and liquid is calculated. As a result, the time-dependent concentration fields are obtained. For a comparison with experimental investigation, the axial dispersion coefficients are calculated. The multifluid model predicts the axial dispersion coefficients for the gas and liquid phases in good agreement with experimental results. Thus, the model is then extended to include mass transfer from
Figure 11. Flow field in a bubble column considering a physical absorption process.
4966
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005
Nomenclature
Figure 12. Time-averaged volume fraction of gas.
a ) interfacial area density [m-1] CD ) drag coefficient c ) molar concentration DS ) column diameter [m] D ) diffusion coefficient [m2 s-1] d ) bubble diameter [m] E ) dispersion coefficient [m2 s-1] Eo ) Eo¨tvo¨s number ) g∆Fd/σ g ) acceleration due to gravity [m s-2] h ) column height [m] H ) Henry coefficient [kg m-1 s-2] j ) superficial velocity [m s-1] k ) turbulent kinetic energy [m2 s-2] M ˙ ) mass flux [kg s-1] m ˘ ) mass flux density [kg m-2 s-1] N˙ ) molar flux [mol s-1] n˘ ) molar flux density [mol m-2 s-1] p ) pressure [kg m-1 s-2] Re ) Reynolds number ) ud/ν Sh ) Sherwwod number ) βd/D r ) coalescence kernel [m3 s-1] t ) time [s] u, v, w ) Cartesian velocity component [m s-1] vi ) bubble volume [m3] Z ) dimensionless breakage frequency Greek Symbols
Figure 13. Instantaneous flow fields for several mass-transfer rates.
R ) volume fraction β ) mass-transfer coefficient [m s-1] ) turbulent kinetic energy dissipation rate [m2 s-3] η ) dynamic viscosity [kg m-1 s-1] λ ) enhancement factor µ ) molecular weight [kg kmol-1] F ) density [kg m-3] ν ) kinematic viscosity [m2 s-1] ξ ) mass fraction σ ) surface tension [kg s-2] ω ) angular frequency [s-1]
Literature Cited
Figure 14. Influence of mass transfer on the bubble diameter.
the gas phase to the liquid phase. For large masstransfer rates, the reduction of the volume fraction of gas influences the flow pattern significantly. Acknowledgment The authors gratefully acknowledge the financial support of the German Research Foundation (DFG). For the calculations, the high-performance computer of the North-German High-Performance Computing Combine (HLRN) in Hannover and Berlin is used.
(1) Deckwer, W.-D. Bubble Column Reactors; John Wiley and Sons: New York, 1992. (2) Warsito, W.; Ohkawa, M.; Maezawa, A.; Uchida, S. Flow Structure and Phase Distribution in a Slurry Bubble Column. Chem. Eng. Sci. 1997, 21/22, 3941. (3) Warsito, W.; Ohkawa, M.; Maezawa, A.; Uchida, S. CrossSectional Distributions of Gas and Solid Holdups in a Slurry Bubble Column Investigated by Ultrasonic Computed Tomography. Chem. Eng. Sci. 1999, 21, 4711. (4) Devanathan, D.; Moslemian, D.; Dudukovic, M. P. Flow mapping in bubble columns using CARPT. Chem. Eng. Sci. 1990, 8, 2285. (5) Dudukovic, M. P.; Larachi, F.; Mills, P. L. Multiphase ReactorssRevisited. Chem. Eng. Sci. 1999, 13-14, 1975. (6) Dudukovic, M. P. Opaque multiphase reactors: Experimentation, modeling and troubleshooting. Oil Gas Sci. Technol. 2000, 2, 135. (7) Joshi, J. B. Computational flow modeling anddesign of bubble column reactors. Chem. Eng. Sci. 2001, 21-22, 5893. (8) Sankaranarayanan, K.; Kevrekidis, I. G.; Sundaresan, S.; Lu, J.; Tryggvason, G. A comparative study of lattice Boltzman and front-tracking finite-difference methods for bubble simulations. Int. J. Multiphase Flow 2003, 1, 109. (9) Sokolichin, A.; Eigenberger, G. Applicability of the standard k- turbulence model to the dynamic simulation of bubble columns: Part I. Detailed numerical simulations. Chem. Eng. Sci. 1999, 13-14, 2273. (10) Pan, Y.; Dudukovic, M. P.; Chang, M. Dynamic simulation of bubbly flow in bubble columns. Chem. Eng. Sci. 1999, 13-14, 2481.
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 4967 (11) Pan, Y.; Dudukovic, M. P.; Chang, M. Numerical investigation of gas-driven flows in 2-D columns. AIChE J. 2000, 3, 434. (12) Pfleger, D.; Gomes, S.; Gilbert, N.; Wagner, H.-G. Hydrodynamic simulation of laboratory scale bubble columns fundamental studies of the Eulerian-Eulerian modeling approach. Chem. Eng. Sci. 1999, 21, 5091. (13) Torvik, R.; Svendsen, H. F. Modelling of Slurry Reactors, a Fundamental Approach. Chem. Eng. Sci. 1990, 8, 2325. (14) Grienberger, J. Untersuchung und Modellierung von Blasensa¨ulen. Ph.D. Thesis, University of Erlangen-Nu¨rnberg, Erlangen, Germany, 1992. (15) Jakobsen, H. A.; Sannaeas, B. H.; Grevskott, S.; Svendsen, H. F. Modeling of vertical bubble-driven flows. Ind. Eng. Chem. Res. 1997, 10, 4052. (16) Ranade, V. V. Modelling of Turbulent Flow in a Bubble Column Reactor. Chem. Eng. Res. Des. 1997, 14. (17) [17] Sokolichin, A.; Eigenberger, G.; Lapin, A.; Lu¨bbert, A. Dynamic numerical simulation of gas-liquid two-phase flows. Euler/Euler versus Euler/Lagrange. Chem. Eng. Sci. 1997, 4, 611. (18) [18] Lain, S.; Bro¨der, D.; Sommerfeld, M.; Go¨z, M. Modelling hydrodynamics and turbulence in a bubble column using the Euler-Lagrange procedure. Int. J. Multiphase Flow 2002, 8, 1381. (19) Krishna, R.; Urseanu, M. I.; van Baten, J. M.; Ellenberger, J. Influence of Scale on the Hydrodynamics of Bubble Columns Operating in the Churn-Turbulent Regime: Experiments vs Eulerian Simulations. Chem. Eng. Sci. 1999, 21, 4903. (20) Sanyal, J.; Vasquez, S.; Roy, S.; Dudukovic, M. P. Numerical simulation of gas-liquid dynamics in cylindrical bubble column reactors. Chem. Eng. Sci. 1999, 21, 5071. (21) Lehr, F.; Millies, M.; Mewes, D. Bubble-size distributions and flow fields in bubble columns. AIChE J. 2002, 11, 2426. (22) Ishii, M.; Kim, S. Development of One-Group Interfacial Area Transport Equation. Nucl. Sci. Eng. 2004, 3, 257. (23) Loth, E. Numerical approaches for motion of dispersed particles, droplets and bubbles. Prod. Energy Comb. Sci. 2000, 3, 161.
(24) Sokolichin, A.; Eigenberger, G.; Lapin, A. Simulation of Buoyancy Driven Bubbly Flow: Established Simplifications and Open Questions. AIChE J. 2004, 1, 24. (25) Clift, R.; Grace, J. R.; Weber, M. E. Bubbles, Drops and Particles: Academic Press: New York, 1978. (26) Brauer, H. Particle/Fluid Transport Processes. Prog. Chem. Eng. 1981, 19, 81. (27) Lopez de Bertodano, M.; Lahey, R. T.; Jones, O. C. Development of a k- model for bubbly two-phase flow. J. Fluids Eng. 1994, 1, 128. (28) Shetty, S. A.; Kantak, M. V.; Kelkar, B. G. Gas-phase backmixing in bubble column reactors. AIChE J. 1992, 7, 1013. (29) Degaleesan, S.; Roy, S.; Kumar, S. B.; Dudukovic, M. P. Liquid mixing based on convection and turbulent dispersion in bubble columns. Chem. Eng. Sci. 1996, 51, 1967. (30) Yang, G. Q.; Fan, L.-S. Axial Liquid Mixing in HighPressure Bubble Columns. AIChE J. 2003, 8, 1995. (31) Levenspiel, O. Chemical Reaction Engineering, 3rd ed.; John Wiley & Sons: New York, 1999. (32) Fogler, H. S. Elements of Chemical Reaction Engineering, 3rd ed.; Prentice Hall: Englewood Cliffs, NJ, 1999. (33) Ohki, Y.; Inoue, H. Longitudinal mixing of the liquid phase in bubble columns. Chem. Eng. Sci. 1970, 1, 1. (34) Qi, M.; Lorenz, M.; Vogelpohl, A. Mathematical solution of the two-dimensional dispersion model. Chem. Eng. Technol. 2002, 7, 693. (35) Mangartz, K. H.; Pilhofer, T. Untersuchungen zur Gasphasendispersion in Blasensa¨ulenreaktoren. Verfahrenstechnik 1980, 1, 40. (36) Hillmer, G. Experimentelle Untersuchung und fluiddynamische Modellierung von Suspensionsblasensa¨ulen. Ph.D. Thesis, University of Erlangen-Nu¨rnberg, Erlangen, Germany, 1993.
Received for review September 2, 2004 Revised manuscript received November 22, 2004 Accepted November 24, 2004 IE049163C
