ENGINEERING DEVELOPMENT OF SLURRY BUBBLE COLUMN REACTOR

(SBCR) TECHNOLOGY

Quarterly Technical Progress Report No. 16

For the Period 1 January – 31 March 1999

FINAL

Contractor
AIR PRODUCTS AND CHEMICALS, INC.
7201 Hamilton Blvd.
Allentown, PA 18195-1501

Bernard A. Toseland, Ph.D.

Program Manager and Principal Investigator

Robert M. Kornosky
Contracting Officer’s Representative

Prepared for the United States Department of Energy

Under Cooperative Agreement No. DE-FC22-95PC95051
Contract Period 3 April 1995 – 2 April 2000

NOTE: AIR PRODUCTS DOES NOT CONSIDER ANYTHING IN THIS

REPORT TO BE CONFIDENTIAL OR PATENTABLE.
 ENGINEERING DEVELOPMENT OF SLURRY BUBBLE COLUMN REACTOR
(SBCR) TECHNOLOGY

Quarterly Technical Progress Report No. 16

For the Period 1 January – 31 March 1999

Contract Objectives

The major technical objectives of this program are threefold: 1) to develop the design tools and a
fundamental understanding of the fluid dynamics of a slurry bubble column reactor to maximize
reactor productivity, 2) to develop the mathematical reactor design models and gain an
understanding of the hydrodynamic fundamentals under industrially relevant process conditions,
and 3) to develop an understanding of the hydrodynamics and their interaction with the
chemistries occurring in the bubble column reactor. Successful completion of these objectives
will permit more efficient usage of the reactor column and tighter design criteria, increase overall
reactor efficiency, and ensure a design that leads to stable reactor behavior when scaling up to
large diameter reactors.

1
WASHINGTON UNIVERSITY IN ST. LOUIS
The report from Washington University for the period follows.

ENGINEERING DEVELOPMENT OF
SLURRY BUBBLE COLUMN REACTOR (SBCR) TECHNOLOGY

Sixteenth Quarterly Report

for
January 1 - March 31, 1999

(Budget Year 4: October 1, 1998 – September 30, 1999)

Objectives for the Fourth Budget Year

The main goal of this subcontract from the Department of Energy via Air Products is to study the
fluid dynamics of slurry bubble columns and address issues related to scaleup and design. The
objectives set for the fourth budget year (October 1, 1998 – September 30, 1999) are listed
below.

• Extension of the CARPT/CT database to conditions of industrial interest such as high

superficial gas velocity (up to 30-50 cm/s).
• Examination of the improved gas mixing phenomenological model against LaPorte tracer
data.
• Critical evaluation of the developed phenomenological models for liquid and gas mixing
against the newly obtained data.
• Testing of the 4-points optical probe for bubble sizes and bubble rise velocity measurements.
• Further improvement in Computational Fluid Dynamics (CFD) using CFDLIB and FLUENT
through development of improved closure schemes and comparison of two-dimensional (2D)
and three-dimensional (3D) model predictions with 2D and 3D data.

In this report, the research progress and achievements accomplished in the sixteenth quarter
(January 1 - March 31, 1999) are summarized.

2
HIGHLIGHTS

Characterization of Gas Phase Mixing Via Gas-Liquid Recirculation Model (GLRM)

• A one-dimensional (1D) model based on the two-fluid approach and the momentum
balances for the gas and liquid phase has been developed. The resulting methodology
allows the de-coupling of the solution for the liquid/slurry phase momentum balance
from the gas phase velocity, which is subsequently calculated based on a suitable drag
formulation.

• The computation of the gas phase velocity profile is based on rigorously satisfying the
gas phase continuity; and is thereby free from shortcomings resulting from the use of an
arbitrary slip velocity correlation for gas velocity computations. Use of such a
correlation may violate the integral mass balance for the gas phase.

• The mathematical description for the physically based two bubble-class model for
representation of gas phase transport in bubble columns has been developed and is
presented in this report.

• The procedure for the estimation of the parameters required in the simulation of the
above model equations has been established and is outlined in this report.

• The developed model presents an intermediate step between the empirical axial
dispersion models and the full-blown 2D and 3D computational fluid dynamic models for
representation of gas flow and transport in bubble columns.

• The model indicates that the time-averaged centerline liquid velocity in the LaPorte
reactor during methanol synthesis may be on the order of 1 m/s, while the gas centerline
velocity could be as large as 1.5 m/s.

Numerical Simulation of 2D Bubble Column and Quantitative Comparison with PIV

Measurements

• The simulations were performed for the 11.2, 15.2, and 10.16 cm columns and for the
flow conditions of 1 to 5.49 cm/s superficial gas velocity, 110 and 160 cm static liquid
height, 2 and 3 gas injectors, and 7.2 to 15.8 two-phase dispersion height to column
diameter aspect ratio.

• Reasonably good agreement between the computed mean quantities--which include

liquid velocity, turbulence intensities and Reynolds shear stress—and the PIV data were
obtained.

3
 • The characteristics of the large structures, i.e., the wave length and frequency, varied with
superficial gas velocity and column dimension. The numerical values were comparable
to those obtained by PIV measurements.

1. CHARACTERIZATION OF GAS PHASE MIXING VIA GAS-LIQUID

RECIRCULATION MODEL (GLRM)

The interpretation of radioactive tracer runs to characterize mixing of the gas phase in a slurry
bubble column during liquid phase methanol synthesis in the LaPorte AFDU reactor was
accomplished earlier using the Axial Dispersion Model (ADM) (Degaleesan et al., 1996a). That
study clearly showed that the gas tracer data interpretation was more difficult than that for the
liquid tracer due to the finite solubility of the tracer (Ar-41). This necessitated modeling of both
the gas and liquid phases simultaneously. The gas dispersion coefficients fitted to the tracer
responses collected at various elevations did not exhibit a consistent trend, and the values were
widely scattered around the mean. This situation was analogous to that observed when the ADM
was used to interpret liquid phase tracer runs. Moreover, attempts to extract other parameters
from the tracer data, such as volumetric mass transfer coefficients, did not produce consistent
results.

In order to remedy this situation, the development of a two-bubble-class, gas-liquid recirculation

model, which is based on experimental evidence, has now been completed to overcome some of
the inherent inadequacies of the ADM. The details of the initial model formulation and
parameter estimation were reported as part of the semiannual report under UCR grant DE FG22
95 PC95212 (Dudukovic' et al., 1996). However, the model development reported there was
only preliminary, and no systematic study was undertaken to examine model predictions as a
function of the input parameters. In view of the possibilities of using the model as a future
scaleup tool, the model has been revisited to incorporate the best available physics to estimate
parameters in order to render the model as predictive as possible. We outline here the model
structure and some of the developments in estimating model parameters. We also present the
detailed model equations, along with the physical basis for this mechanistic model.

1.1 Two-Bubble-Class, Gas-Liquid Recirculation Model

The compartmentalization of the developed mechanistic (phenomenological) model is shown in
Figure 1.1 for completeness. The liquid is a batch system, with up-flow (L1) in the core and
down-flow (L2) near the walls. The extension to situations with net co-current flow of the
liquid/slurry is straightforward. The gas phase also has a similar recirculation pattern; however,
only small bubbles (SB2), which do not possess sufficient momentum, get recirculated along
with the liquid. The large bubbles (LB), with smaller bubbles (SB1) trapped in the wakes of
these larger bubbles, rise up through the central core due to large buoyancy forces and drag the
liquid along with them. The top (disengagement) zone and the bottom (distributor) zone are
modeled as CSTRs with uniform (small) bubbles only. The height of these zones is taken to be
the same as the column diameter Dc. Changing the height of these zones between 0.5 and 1.0 Dc
does not have any noticeable effect on the liquid backmixing [Degaleesan et al., 1996b], as long

4
as the height of the gas-liquid mixture is much greater than the column diameter (L/Dc of at least
6). The effect of these zones on the gas phase model still needs to be verified.

For the developed part of the flow, which occupies most of the column, five transient
convection-diffusion equations with mass transfer between regions as source terms can be
written. Additionally, reaction rates appear as source terms in the liquid phase equations if one is
simulating a column under reaction conditions. The other source terms come from the large-
small bubble interactions in the core of the column.

The developed phenomenological model, schematically shown in Figure 1.1, is based on the
currently available physical evidence. The hydrodynamic phenomena observed experimentally
through the unique CARPT and CT techniques are represented schematically in Figure 1.2 and
form the basis of the mechanistic model described earlier. In a bubble column in the time
average sense, a single liquid recirculation loop is established with the liquid rising in the center
and flowing downwards by the walls (Degaleesan, 1997). Superimposed on this recirculation is
eddy diffusion caused by the bubble wakes. Gas travels in the column via small and large
bubbles (Krishna and Ellenberger, 1996), which exchange their contents in the up-flow region.
However, only small bubbles are dragged down by the downward flowing liquid. A soluble-gas-
tracer experiences mass transfer to the liquid phase, and the mass transfer coefficients for the
small and large bubbles are different and have to be evaluated separately.

Figure 1.2 shows the representative liquid and gas velocity profiles. r' and r'' are the radial
locations at which the time-averaged liquid and gas velocities change sign, respectively, i.e., go
to zero. As the gas and liquid velocities reach zero at different radial locations, there is a region
between r' and r'' where small bubbles are moving up but the liquid is moving down.

1.2 Two-Fluid Model for Gas and Liquid Phase Momentum Exchange

The basis for the derivation of the 1D model for the liquid- and gas-phase-velocity profiles are
the two-fluid model equations presented below. These are the result of the ensemble-averaged
approach of Drew (1983). The subscript 'c' denotes the continuous liquid/slurry phase, whereas
the subscript 'd' denotes the dispersed gas phase. Both phases are considered incompressible.

Equations of Continuity

∂ε c
Liquid/Slurry + ∇ .(ε c u c ) = 0 (1.1)
∂t
∂ε d
Gas + ∇ .(ε d u d ) = 0 (1.2)
∂t

5
Momentum Equations

Liquid/Slurry:
 ∂u 
ρ c ε c  c + u c .∇ u c  = ρ c ε c g − ε c ∇p − (M d + M vm ) + ∇.(ε c σ c ) (1.3)
 ∂t 
Gas:
 ∂u 
ρ d ε d  d + u d .∇ u d  = ρ d ε d g − ε d ∇p + (M d + M vm ) + ∇.(ε d σ d ) (1.4)
 ∂t 
In the momentum balance equations, Md is the drag force term, while Mvm is the virtual mass
term defined below (Drew, 1983).

6ε cε d
Md = Fd (1.5)
πd p3

with
1
Fd = ρ cπd p2 C D u c − u d (u c − u d ) (1.6)
8
 24
C D = max  (
1 + 0.15 Re 0.687 , f )
8 Eo 
3 Eo + 4 
(1.7)
 Re
2
 1 + 17.67ε c9 7 
f =  (1.8)
 18.67ε c
32

where
Eo ≡ Eotvos Number ≡ gρ c d p2 γ

and Re ≡ Bubble Re ynolds Number

≡ ρ c d p uc − ud µ c

1  Du c Du d 
M vm = ε c ε d C vm  −  (1.9)
2  Dt Dt 

C vm = 1 + 3.32ε d + O( ε d2 ) (1.10)

In the well-developed flow region of the column, the flow in the time-averaged sense is known
to be axisymmetric, with only the axial velocities being non-zero. Hence, the time-averaged

6
liquid flow pattern is represented by a single radial profile of the axial velocity. These
assumptions are fair in view of the holdup profile database available at CREL via Computed
Tomography (CT); and the liquid velocity profile database via Computer Automated Radioactive
Particle Tracking (CARPT).

Under these assumptions, the equations of continuity for both the phases (Equations 1.1 and 1.2)
are identically satisfied, and one cannot use the traditional approach of solving the Poisson
equation for the pressure correction via these continuity equations (as is done in 2D and 3D CFD
codes). In addition, the left-hand side of the momentum equations for both the gas and liquid
phase becomes zero, as does the virtual mass term. Finally, due to no flow in the radial and
azimuthal directions, the pressure becomes independent of the radial and azimuthal coordinates,
and the pressure gradient term in the momentum equations reduces to dP/dz.

After retaining the non-zero gradients and velocity components in the momentum equations for
the two phases, one obtains the following simplified equations:

dP 1 d
Liquid/Slurry 0 = − ρ cε c g − ε c − Md + ( rε cσ c ,rz ) (1.11)
dz r dr
dP 1 d
Gas 0 = −ρdε d g − ε d + Md + ( rε d σ d ,rz ) (1.12)
dz r dr

Adding Equations 1.11 and 1.12, one obtains

 negligible
678  dP 1 d
0 = − ρ d ε d + ρ cε c  g −
 + ( rτ rz ) (1.13)
  dz r dr
 

where

τ rz = ε cσ c ,rz + ε d σ d ,rz (1.14)

is the effective shear stress of the two-phase mixture. However, since the gas phase viscosity
and density are two to three orders of magnitude lower than the liquid/slurry values, the effective
shear stress is primarily due to the shear stress of the liquid /slurry phase. Thus, τrz is the
liquid/slurry phase shear stress. A model is needed for this stress, and the simplest closure in
terms of turbulent kinematic viscosity is employed in Equation 1.16.

1 d dP
( rτ rz ) = + ρ c (1 − ε d )g (1.15)
r dr dz

7
 2
du du  du 
τ rz ( r ) = − ρ c (ν c ,m + ν c ,t ) c = − ρ cν c ,m c − ρ c l 2  c  (1.16)
dr dr  dr 

The turbulent eddy viscosity, νc,t, is closed by a modified mixing length, l(ξ), as given by Kumar
(1994).

du c
ν c ,t = l 2 (1.17)
dr
a(1 − ξ )
l (ξ ) = + d (1 − ξ )
e
(1.18)
(ξ + b ) c

The parameters a, b, c, d and e were obtained by Kumar et al. (1994) after extensive data on
liquid recirculation velocities from CARPT were considered, as well as results from experiments
of other researchers who have made measurements of the liquid recirculation velocity by other
experimental means.

At this point, the computation of the radial profile of the axial liquid velocity requires only the
gas holdup profile as the input. The circumferentially averaged radial gas holdup profile is
frequently represented as follows:

 m+2 
ε d (ξ ) = ε d   1 − cξ
m
( ) (1.19)
 m + 2 − 2c 

With the gas holdup profile as input, the liquid velocity profile can be readily computed by a
procedure described elsewhere (Kumar, 1994; Kumar et al., 1994).

Once the liquid velocity profile and dP/dz are determined from the converged 1D liquid
circulation model (Kumar et al., 1994), one returns to the gas phase momentum equation and
iterates on the bubble diameter in the drag formulation. In this way, the converged gas phase
velocity profile that satisfies the gas phase continuity is obtained.

dP
Md = εd + ρdε d g (1.20)
dz
6ε cε d
Md = Fd (1.21)
πd p3

1
Fd = − ρ cπd p2 C D (u c − u d )
2
(1.22)
8

8
Substitution of Equation (1.21) into Equation (1.20) with the help of Equation (1.22) results in
the equation from which the gas velocity profile, ud, is determined. This calculation requires
iteration, since the bubble size, dp, must be guessed and the drag coefficient must be evaluated by
Equation 1.7. The guessed bubble size, dp, is used along with the available gas holdup profile,
εd(ξ), from Equation 1.19 and the liquid recirculation velocity calculated from the model,
consisting of Equations 1.15-1.19, to obtain the updated ud. Iteration on the bubble diameter, dp,
continues until the integral gas phase continuity is satisfied.

It should be mentioned that the use of a single bubble size during iteration for the gas velocity
profile is to demonstrate the applicability of this novel technique. To be consistent with the
overall two-bubble class approach, the existence of two different bubble classes needs to be
reflected in the computation of the gas velocity profile as well. A variant to this model is being
developed in which the gas phase momentum will be split into separate contributions from the
small and large bubbles, along with a momentum interchange term for the two bubble classes.
This work is still under development and will be reported upon completion.

1.3 Gas Phase Mixing Model

1.3.1 Assumptions and Model Development
1) The process is semi-batch, but the development can be extended readily for continuous
liquid/slurry flow.
2) The interaction between small bubbles moving up (SB1) and small bubbles moving down
(SB2) is due to radial turbulent diffusion alone.
3) The small bubble coalescence and breakup in the SB2 region where bubbles are dragged
down by the liquid is negligible, as no large bubbles are observed near the wall. To
communicate with the large bubbles (LB), the bubbles from SB2 have to be transported
radially via turbulent diffusion to the region SB1, where they can interact with large
bubbles.
4) The axial and radial turbulent diffusivities for the small bubbles and the liquid are
assumed to be the same. (If the instantaneous gas velocities are known, one can do away
with this assumption, and calculate the gas turbulent diffusivities in the same manner as
for the liquid.)
5) The interaction between large bubbles in the upflow region, LB, and small bubbles in the
upflow region, SB1 is due to bubble breakup and coalescence. The mass transfer
coefficients (k) for the two bubble classes are calculated based on Higbie’s relation, by
calculating the contact times based on average bubble diameters and average gas
velocities through each bubble phase.
6) Reactions, if any, occur only in the liquid phase. Gas side mass transfer resistance is
negligible.

With these assumptions, the 1D mass balances written for all the regions, excluding the CSTRs
at the top and bottom, yield five transient, turbulent diffusion-convection partial differential
equations (PDEs) with interphase transport, which will be described in the next section. Here we
outline the details of the mass balances for one of the compartments of the column into which the
fully developed region of the reactor is sectioned, viz., the portion of the column where the small
9
bubbles move upwards. The source terms are either due to gas-liquid mass transfer or are from
exchange due to bubble coalescence and re-dispersion. The exchange between the liquid going
up and the liquid going down is expressed in terms of the radial turbulent diffusivity, which is
also assumed to characterize the exchange between small bubbles rising and small bubbles
travelling downwards.

Figure 1.3 presents the various inputs and outputs to an element of the fully developed region,
SB1, of height dx. The various terms are mathematically represented below (Refer to Figures 1.2
and 1.3).

 ∂C g 1 
(Convective + turb. diffusive input)x (
: u g 1π (r' ' ) ε g 1C g 1
2
) − D xx1 π (r' ' ) ε g 1 
2

 ∂x  x
x

 ∂Cg1 
(
(Convective + turb. diffusive output)x + dx : u g1π (r' ' ) ε g1Cg1
2
) x + dx
− Dxx1 π( r '') 2
ε g1  
 ∂x  x + dx

: k su a sulu ( HCg 1 − Cl 1 )( π( r ' ') dx )

2
Mass transfer to liquid moving upwards

Mass transfer to liquid moving downwards : k su a suld ( HCg 1 − Cl 2 )( π( r ' ') dx )

2

Net radial turbulent diffusional exchange between SB1 and SB2 (output):

[ D ε ] (C
rr g
r = r ''
g1 )
− Cg 2 dx

Net source term for SB1 from bubble breakage and coalescence (input):

 π *3  π 3
 d S  [ − DSCL ] ( π(r ' ') 2 dx ) C g1 +  d L*  BSBL ( π(r ' ') 2 dx ) C g 3
6  6 

where d S* and d *L are the Sauter mean diameter of the small and large bubble classes,
respectively, and may vary along the height of the reactor.

Therefore, Accumulation = ( Inpout ) − (Output ) + Generation gives,

10
  ∂ 2 C g1 [
∂C g 1 Drr ε g r = r '']( ) k su a suld
( )

 D xx 1 − u − C − C −. HC − C 
∂C g 1 
=
∂x 2
g 1
∂x 2
( )
π r '' ε g 1
g 1 g 2
ε g1
g 1 l 2

 (1.23)
∂t  k su a sulu  π * 3  C g 1   π * 3  C g 3 
− ε (HC g1 − C l1 ) − DSCL  d S   + BSBL  d L   
 g1 6  ε g 1  6  ε g 1  

Similarly, for the other four regions we obtain the following four transient convection-diffusion
equations.

For the small bubbles moving downwards (SB2)


Dxx
∂2Cg2
+ ug2
∂Cg2
+
1 Drr ε g ' '
r =r
(
[ ] 
Cg1 − Cg2 )
∂Cg2 
=
2
∂ x 2
∂x 2
(
π R − (r'' ) 2
ε g2 ) 
 (1.24)
∂t  k a  
−  sd sdld (( HCg2 − Cl 2 )) 
  ε g2  

For the large bubbles moving upwards (in plug flow)(SB3)

 ∂C g 3  k Lu a Lu  
− ug 3
∂C g 3  ∂x
−
ε
 g 3 
(
 ( HC g 3 − Cl 1 ) 

)
=  (1.25)
∂t
 +  π d * 3 B C − π d * 3 D C  1 
  6 S LCS g 1 6 L LBS g 3  ε g 3 
 

For the liquid moving upwards (L1)


∂2C D ε [ ] 

( )
 ∂C
l1 − u l1 − 1 rr l r = r C − C − Rx 
'
D 
∂C xx ∂x 2 l1 ∂x π(r ') 2 ε l1 l2 l1 
l1 =  1 l1 (1.26)
∂t  
  r ''  2 1    
+     k a ( HC − C ) + k a ( HC − C ) 
  r ' ε   Lu Lu g3 l1 su sulu g1 l1  
l1  
11
For the liquid moving downwards (L2)

 ∂2Cl 2 ∂Cl 2 1 [
Drr εl r =r ' ] 
Dxx + u +
∂x π( R2 − (r')2 )
( C − C ) − Rx 
∂Cl 2  2 ∂x εl 2
2 l2 l1 l2 l2

=  (1.27)
∂t   
ksd
( 
+  ε ( R2 − (r')2 )  (r'') aSulu ( HCg1 − Cl 2 ) + ( R − (r'') ) aSdld ( HCg2 − Cl 2 ) 
2 2 2
)
  l 2  

It must be noted that these equations are valid only in the regions excluding the two CSTRs. The
CSTR at the distributor has a volume Va and that at the top, a volume equal to Vb. Each CSTR
can be split into two volumes, one occupied by the gas and other by the liquid. Thus, applying
the mass balance to these volumes, we obtain the following:

Distributor zone
For the liquid phase:

=
(
dCla ε l 2 R 2 − ( r ' )2 ul 2 )
Cl 2 ,0 −
ε l 1 r ' ul 1
2

Cla + CSTR CSTR (HCga − Cla ) − Rxla

k a
(1.28)
dt εl R 2
DC ε l R DC 2
εl

Ql1 Cla Ql2 Cl2,0 (Qg1 + Qg3) Cga Qg2 Cg2,0

ε lVa ε gVa

Qg Cg,in
For the gas phase (with no gas phase reaction):

(
 ε g 2 R 2 − ( r '' ) 2 u g 2 ) −
(ε g 1u g 1 + ε g 3 u g 3 ) r '' 
2

 C g 2 ,0 C ga 
dC ga ε R 2
D ε D R 2

= 
g C g C
 U  (1.29)
dt
+ G
C g ,in −
k CSTR a CSTR
(HC ga − C la ) 
 ε g DC εg 

12
Disengagement zone Qg Cgb

ε lVb ε gVb

Ql1 Cll,L Ql2 Clb (Qg1 Cg1,L + Qg3 Cg3,L) Qg2 Cgb

For the liquid phase:

=
2
dClb ε l 1 r ' u l 1
C −
(
ε l 2 R 2 − ( r ' )2 ul 2 )
Clb + CSTR CSTR (HC gb − Clb ) − Rxlb
k a
(1.30)
l 1,L
dt ε l R DC2
εl R 2
DC εl

For the gas phase (with no gas phase reaction):

 (ε g 1u g 1C g 1,L + ε g 3 u g 3 C g 3 ,L ) r '' ( )
ε g 2 R 2 − ( r '' ) 2 u g 2 
2

 − C gb 
dC gb  ε g DC R 2
εg R 2
DC 
= (1.31)
dt  
− G C gb − CSTR CSTR (HC gb − C lb )
U k a

 ε g DC εg 

The equations above are ordinary differential equations (ODEs) and require only the initial
conditions. For the fully developed sections, initial as well as boundary conditions are required.
These are discussed next.

Initial Conditions

A step input of gas tracer at the inlet (bottom) of the column is assumed, with no tracer originally
present in the column.

t = 0; Cg,in = H(t) (1.32)

Cla = Clb = Cga = Cgb = Cl1 = Cl2 = Cg1 = Cg2 = Cg3 = 0 (1.33)

13
Alternatively, an experimental impulse input can be simulated with a Gaussian function with a
tail. Details on this alternative can be found in Degaleesan (1997).

Boundary conditions for the fully developed region

Danckwerts boundary conditions are used at inlet and exit, guaranteeing preservation of mass
balance for each phase.

Upflow section of the liquid

At the bottom of the fully developed flow zone, i.e., at the boundary with the CSTR
representing the distributor zone:

∂C l 1
x=0; ul1Cla = ul1Cl1 | x = 0 − Dxx1 | (1.34)
∂x x = 0
At the top of the fully developed flow zone, i.e., at the boundary with the CSTR representing
the disengagement zone:
∂C l 1
x=L; | =0 (1.35)
∂x x = L

Downflow section of the liquid

At the top of the fully developed flow zone, i.e., at the boundary with the CSTR representing
the disengagement zone:
∂Cl 2
x=L; ul 2 Clb = ul 2 Cl 2 |x = L + Dxx2 |x = L (1.36)
∂x
At the bottom of the fully developed flow zone, i.e., at the boundary with the CSTR
representing the distributor zone:
∂C l 2
x=0; | =0 (1.37)
∂x x = 0

Upflow of small bubbles

At the bottom of the fully developed flow zone, i.e., at the boundary with the CSTR
representing the distributor zone:
∂C g 1
x=0; u g1C ga = u g1C g1 | x = 0 − Dxx1 | (1.38)
∂x x = 0

14
 At the top of the fully developed flow zone, i.e., at the boundary with the CSTR representing
the disengagement zone:
∂C g 1
x=L; | =0 (1.39)
∂x x = L

Downflow of small bubbles

At the top of the fully developed flow zone, i.e., at the boundary with the CSTR representing
the disengagement zone:
∂Cg 2
x=L; ug 2 Cgb = ug 2 Cg 2 |x = L + Dxx2 |x = L (1.40)
∂x
At the bottom of the fully developed flow zone, i.e., at the boundary with the CSTR
representing the distributor zone:
∂C g 2
x=0; | =0 (1.41)
∂x x = 0

Upflow of large bubbles

At the bottom of the fully developed flow zone, i.e., at the boundary with the CSTR
representing the distributor zone:

x=0; C g 3 | x = 0 = C ga (1.42)

1.3.2 Parameter Estimation

The estimation of the various model parameters is outlined in this section. Some of these
parameters are estimated based on literature correlations; others are calculated by averaging the
velocity and holdup profiles obtained by CARPT/CT.

I. Parameters Estimated from Correlations

Holdup of small bubbles

(
ε small = 0.5 exp − 193ρ −g 0.61η0l .5σ 0.11 ) (Wilkinson et al., 1992) (1.43)

Sauter mean diameter of small bubbles

dS* = 3g −0.44σ0.34 η0l .22ρl−0.45ρ−g0.11U G−0.02 (Wilkinson et al., 1994) (1.44)

15
Superficial velocity of small bubbles

−0.273 0.03
σ σ ρl 
3
 ρl 
Udf = 2.25ε small     (Wilkinson et al., 1992) (1.45)
ηl  gηl4  ρ 
 g

Sauter mean diameter of large bubbles

( )
25
d L* = U G − U df (h * + h0 ) 4 5 g − 1 5 (de Swart, 1996) (1.46)

Bubble size stabilization height above the gas distributor

h* = 0.73 (U G − U df ) (de Swart, 1996) (1.47)

The gas holdup is assumed to vary radially as (Kumar, 1994)

 m + 2
ε g (ζ ) = ε g   (1 − c ζ m ) (1.48)
 m 

where ζ is the dimensionless radius. With the holdup profile parameter ( ε g , m and c) being
either determined experimentally or estimated by the procedure outlined by Degaleesan (1997),
the gas and liquid velocity profiles are computed by the methodology outlined in the previous
monthly report.

II. Computed Parameters

The fraction of the inlet flow rate of the gas, which travels up the column as large bubbles, is
then computed as

U G − U df
α= (1.49)
UG

16
With the knowledge of the radial profile of gas holdup and gas and liquid velocity profiles, the
following average quantities can be computed:

Average holdup of the liquid ascending

r'
2
ε l1 =
r' 2 ∫(1 − ε
0
g )rdr (1.50)

Average holdup of the liquid descending

R
2
ε l2 = 2
(
R − ( r' ) 2 )∫(1 − ε
r'
g )rdr (1.51)

Average holdup of the gas ascending

r ''
2
ε g1 + ε g 3 =
r ' '2 ∫ε
0
g rdr (1.52)

This equation cannot give us the two holdups separately and, therefore we need another
independent equation. In the absence of independent measurements, the holdup of either large or
of small bubbles can be calculated using a literature correlation. For this work, equation 1.43 is
used to calculate the holdup of small bubbles.

Average holdup of the small bubbles descending

R
2
( R − (r '' )2 ) r∫'' g
εg2 = 2 ε rdr (1.53)

Average holdup of the small bubbles ascending with liquid descending

''
r
2
(r ' ' )2 r∫ g
ε '
g1 = ε rdr (1.54)
'

Average velocity of the liquid ascending

r'

2 ∫ (1 − ε g ) ul r dr
ul1 = 0
2 (1.55)
r ' εl1

17
Average velocity of the liquid descending

2 ∫ (1 − ε g ) ul r dr
ul 2 = r'
(1.56)
( R 2 − (r ' ) 2 ) ε l 2
Average velocity of the small bubbles ascending

(1 − α ) R 2U G + 2 ∫ u g ε g rdr
r ''
ug1 = '' 2
(1.57)
r ε g1

Average velocity of the large bubbles ascending

α R 2U G
ug3 = 2 (1.58)
r '' ε g 3
Average velocity of the small bubbles descending
R

2 ∫ ug ε g rdr
r ''
ug 2 = (1.59)
(R 2
− (r '' )2 ) ε g 2

Interfacial area for mass transfer from small bubbles moving up to liquid moving up

a sulu =
(
6 ε g 1 − ε g' 1 ) (1.60)
*
d S

Interfacial area for mass transfer from small bubbles moving up to liquid moving down

a su ld =
( )
6 ε 'g 1
(1.61)
d S*

Interfacial area for mass transfer from small bubbles moving down to liquid moving down

a sdld =
( )
6 ε g2
(1.62)
d S*

Interfacial area for mass transfer from large bubbles to liquid moving up

a Lu =
( )
6 ε g3
(1.63)
d L*

18
Interfacial area for mass transfer in CSTRs

a C ST R =
( )
6 εg
(1.64)
d S*

Mass transfer coefficient from small bubbles to liquid moving up

2 DL ,m ug1
k su = (1.65)
π d S*

Mass transfer coefficient from small bubbles to liquid moving down

2 DL ,m ug 2
k sd = (1.66)
π d S*

Mass transfer coefficient from large bubbles to liquid moving up

2 DL ,m ug 3
k Lu = (1.67)
π d L*

Mass transfer coefficient in CSTRs

2 D L ,m U G
k CSTR = (1.68)
π ε g d S*

1.4 Remarks

The protocol outlined above has been established to evaluate all parameters required as input in
the solution of the two-bubble class model equations, for simulation of tracer and/or reactor
performance (reaction is considered only in the liquid phase). As can be seen from these
equations, the number of required hydrodynamic inputs increases dramatically the physics of the
flow in such systems is explored in detail. This emphasizes the fact that the coupling between
mass and momentum transfer in such systems is very complicated and is precisely the reason for
the inadequacy of the Axial Dispersion Model (ADM) for predictive purposes. It is therefore
necessary to try to incorporate as much physics as is possible into the mathematical models to
represent the phenomena in slurry bubble columns.

A critical evaluation of some of the parameters outlined above for the case of a typical operating
condition at the LaPorte AFDU is in progress. Some alternative forms and procedures for
19
improved parameter estimation based on the latest reports in the literature are being explored.
This includes an analysis of the consistency of the evaluated parameters and the form of the
convection diffusion PDEs.

1.5 Results and Discussion

In this report, simulation results for the estimated liquid and gas velocity profiles are presented
for the case of the LaPorte reactor operated under high-pressure conditions. The operating
conditions for the various runs are listed in Table 1.1. The column diameter for this case was
0.46 m. At the operating conditions for the runs, the effective slurry viscosity was 0.99 cP, and
the effective slurry density was 995 kg/m3.

Table 1.1. Operating Conditions during Tracer Runs

Run Op. Pressure Op. Temperature Ugsup,avg εd
0
(MPa) ( C) (cm/s)
14.6 5.2 250 22.86 0.39
14.7 5.2 250 12.66 0.33
14.8 3.6 250 32.81 0.38

For these experiments, Differential Pressure (DP) and Nuclear Density Gauge (NDG)
measurements were conducted at each operating condition. The DP measurements provide the
average gas holdup ( ε d ), whereas from NDG measurements the chordal average holdup along
the diameter is obtained. The constant c in the gas holdup profile of Equation (1.19) is obtained
(assuming a value of 2 for the exponent m as reported in the previous study by Degaleesan et al.,
1996). This yields the parameters in the gas holdup profile as listed in Table 1.2.

Table 1.2. Parameters in the Gas Holdup Profile

Run εd M C

14.6 0.39 2 0.844

14.7 0.33 2 0.891
14.8 0.38 2 0.943

Figures 1.4 and 1.5 show the liquid and gas velocity profiles, respectively, computed from the
one-dimensional momentum balance equations for the liquid and gas phases as illustrated above.
The figures show that the 1D model has been able to effectively capture the recirculation of both
the gas and liquid phases. The model predicts a centerline slip velocity of about 50 cm/s, which
appears to be high. Such a high slip velocity could be due to several causes:
20
 • Only two independent measurements are available to calculate the three parameters in the
gas holdup profile (Equation 1.19). Thus, the profile cannot be determined uniquely, and
one of the parameters has to be fixed to evaluate the other two. For the present case the
parameter chosen was the exponent m, whose value was fixed at 2 as is customarily
reported in the literature for columns operated under churn-turbulent flow conditions at
atmospheric pressures. The assumption of m being 2 might not be reasonable based on
our experience with a laboratory-scale, high-pressure bubble column, in which values of
m in the range of 2.9-3.5 have been normally observed. Since the AFDU is operated at
pressures exceeding the highest operating pressure of our laboratory unit, the use of m =
2 becomes questionable. The exponent m determines the gradient of the gas holdup
profile, which is the main driving force for liquid recirculation, and thus, influences the
calculation of the pressure drop as well. Therefore, the calculation of the slip velocity
may be indirectly influenced by an inaccurate choice of the gas holdup profile. However,
the purpose of these results is to demonstrate what the model can do, provided a correct
gas holdup profile is available.

• A single bubble size has been used to represent the gas phase in the entire cross section.
However, this is not the case, as evidenced by experimental observations. As mentioned
earlier, a modified model is being derived to account for the presence of two bubble
classes.

Additionally, further investigations are under way to evaluate the possible importance of other
terms in the interaction forces between the phases. However, at present, it is believed that drag is
the main contribution to the interphase interaction.

The liquid and gas velocity profiles, along with other input parameters, provide the information
necessary for the simulation of tracer responses. The complete mathematical representation of
the two-bubble class gas-liquid recirculation model, as well as the procedure for parameter
estimation, will be provided in subsequent monthly reports.

Nomenclature
a interfacial area, cm-1
B number of bubbles formed per unit volume per unit time, # cm-3 s-1
C concentration, moles/cm3
c parameter in the holdup profile to allow non-zero holdup at the wall
D number of bubbles disappearing per unit volume per unit time, # cm-3 s-1
DC column diameter, cm
DL,m molecular diffusivity, cm2/s
Drr radial turbulent diffusivity, cm2/s
Dxx1 axial turbulent diffusivity of small bubbles and liquid ascending, cm2/s
Dxx2 axial turbulent diffusivity of small bubbles and liquid descending, cm2/s
21
d* mean bubble diameter, cm
g acceleration due to gravity, cm2/s
H Henry’s constant
h* height above the gas distributor, m
h0 parameter determining the initial bubble size, h0 = 0.03 m
k mass transfer coefficient, cm/s
L dispersion height between the two CSTRs, cm
m power law exponent in the radial gas holdup profile
P operating pressure, dyne/cm2
Q gas flow rate, cm3/s
R column radius, cm
Rx reaction rate (in the liquid phase), moles cm-3 s-1
r radial position in the column, cm
r' radius at which the liquid velocity profile inverts
r'' radius at which the gas velocity profile inverts
T operating temperature, K
t time, s
tc contact time of the liquid eddies with the bubbles, s
UG gas superficial velocity, cm/s
UL liquid superficial velocity, cm/s
Udf gas superficial velocity through small bubbles, cm/s
u velocity, cm/s
u radially averaged mean velocity, cm/s
ul local liquid velocity in the liquid upflow region, cm/s
Va volume in CSTR A, cm3
Vb volume in CSTR B, cm3
VbL rise velocity of large bubbles in quiescent liquid, cm/s
VbS rise velocity of small bubbles in quiescent liquid, cm/s
x axial position in the column, cm

Greek Symbols
α fraction of the inlet gas flow rate going through large bubbles
ε phase holdup
ε radially averaged phase holdup
ε g1
'
mean holdup of small bubbles ascending with liquid descending
γ surface tension of the liquid, Pa-m
ηl slurry viscosity, cP
ρ density, kg/m3
σ stress tensor
υg kinematic viscosity of the gas, cm2/s
υm kinematic viscosity of the liquid, cm2/s
υt turbulent viscosity of the liquid, cm2/s
τ effective shear stress

22
ζ dimensionless radius
ζ' dimensionless radius at which the liquid velocity profile inverts
ζ'' dimensionless radius at which the gas velocity profile inverts

Subscripts
b bubble
CSTR CSTRs A and B
g gas
g1 small bubbles ascending
g2 small bubbles descending
g3 large bubbles ascending
L large bubbles
LBS large bubbles breaking into small bubbles
LCS large bubble from coalescence of small bubbles
Lu large bubbles ascending with liquid descending
l liquid
la liquid phase in the distributor zone, CSTR A
lb liquid phase in the disengagement zone, CSTR B
l1 liquid ascending
l2 liquid descending
S small bubbles
LBL large bubbles due to breakage of/into large bubbles
LBS large bubbles due to breakage into small bubbles
LCL large bubbles due to coalescence of/into large bubbles
LCS large bubbles due to coalescence of small bubbles
SBL small bubbles from breakage of large bubbles
SBS small bubbles due to breakage of/into small bubbles
SCS small bubbles due to coalescence of/into small bubbles
SCL small bubbles due to coalescence into large bubbles
sd small bubbles descending
sdld small bubbles descending with liquid descending as well
su small bubbles ascending
suld small bubbles ascending with liquid descending
sulu small bubbles ascending with liquid ascending as well

References

Degaleesan, S., “Fluid dynamic measurements and modeling of liquid mixing in bubble columns,”
D.Sc. Thesis, Washington University (1997).
Degaleesan, S., M. P. Dudukovic’, B. L. Bhatt, and B. A. Toseland, “Slurry bubble column
hydrodynamics: tracer studies of the LaPorte AFDU reactor during methanol synthesis,”
Fourth Quarterly Report for Contract DOE-FC 2295 PC 95051 (1996a).

23
Degaleesan, S., S. Roy, S. B. Kumar, and M. P. Dudukovic’, “Liquid backmixing based on
convection and turbulent dispersion in bubble columns,” Chem. Engng Sci., 51, 1967 (1996b).
Drew, D. A., “Mathematical modeling of two-phase flow,” Ann. Rev. Fluid Mech., 15, 261 (1983).
Dudukovic', M. P., L. S. Fan, Y. B. Yang, S. D. Degaleesan, and P. Gupta, “Novel techniques for
slurry bubble column hydrodynamics,” Second Semiannual Report for Contract DOE-FG
2295 PC 95212 (1996).
Krishna, R. and J. Ellenberger, “Gas holdup in bubble column reactors operating in the churn-
turbulent flow regime,” AIChE J., 42(9), 2627 (1996).
Kumar, S. B., N. Devanathan, D. Moslemian, and M. P. Dudukovic', “Effect of scale on liquid
recirculation in bubble columns,” Chem. Engng Sci., 49(24B), 5637 (1994).
Kumar, S. B., “Computed tomographic measurements of void fraction and modeling of the flow in
bubble columns,” Ph.D. Thesis, Florida Atlantic University (1994).
de Swart, J. W. A., “Scaleup of a Fischer-Tropsch slurry reactor,” Ph. D. Thesis, University of
Groningen, The Netherlands (1996).
Wilkinson, P. M., H. Haringa, and L. L. van Dierendonck, “Design parameters estimation for scaleup
of high pressure bubble columns,” Chem. Eng. Sci., 49, 1417 (1994).
Wilkinson, P. M., A. P. Spek, and L. L. van Dierendonck, “Mass transfer and bubble size in a column
under pressure,” AIChE J., 38, 544 (1992).

Cgb Clb
x=L

LB SB1 SB2 L2 L1
Cg3 Cg1 Cg2 Cl2 Cl1

x=0

Qg1
Qg3 Cga Ql2 Cla Ql1
Qg2

Qg Cg, in
Figure 1.1 Schematic Representation of the Compartmentalization of the Reactor

24
 1-εL(r)

Gas
Velocity Dzz
Profile

Liquid Gas
Velocity Holdup
Profile Profile

uL (r)
r'
Drr
r''

0 R Liquid
Velocity
Profile

-R 0 R
Figure 1.2 Schematic Representation of the Experimentally Observed Phenomena in Bubble Columns

(Convective + turbulent diffusive output)x+dx

Death due to coalescence (output)
Mass transfer to liquid
moving upwards (output)
Turbulent diffusion to SB2 (output)
dx
Mass transfer to liquid
Birth due to breakage (input) moving downwards (output)

(Convective + turbulent diffusive input)x

Figure 1.3 Details of 1D Mass Balance for the Small Bubbles Ascending

25
 Radial Liquid Velocity Profiles
120
Run_14.6
90
Liquid Velocity (cm/s)

Run_14.7
60 Run_14.8
30

-30

-60

-90
0.0 0.2 0.4 0.6 0.8 1.0
Dimensionaless Radius (r/R)

Figure 1.4 Computed 1D Liquid Velocity Profiles for the Three Operating Conditions

Radial Gas Velocity Profiles

180
150 Run_14.6
120 Run_14.7
Gas Velocity (cm/s)

90 Run_14.8
60
30
0
-30
-60
-90
0.0 0.2 0.4 0.6 0.8 1.0
Dimensionaless Radius (r/R)

Figure 1.5 Computed 1D Gas Velocity Profiles for the Three Operating Conditions

26
2. NUMERICAL SIMULATION OF 2D BUBBLE COLUMN AND QUANTITATIVE
COMPARISON WITH PIV MEASUREMENTS

Recently Lin et al. (1996) and Mudde et al. (1997) used PIV in their experimental studies of two-
dimensional bubble columns. They provided the detailed measurements of liquid velocity and
turbulence intensities for columns of different sizes and under different operating conditions.
They also studied the characteristics of the macroscopic flow structures, i.e., the central
meandering plume and the companion vortical regions, by measuring their frequency,
wavelength and moving speed. This information yields a better understanding of the fluid
dynamics in a 2D bubble column and provides a database for further numerical investigations.

Most of the previous numerical studies reported in the literature compare the predictions with the
experimental measurements in a qualitative manner, while only a few quantitative comparisons
are made. Therefore, although qualitative comparisons seem satisfactory in general, limited
conclusions regarding the validation and reliability of numerical predictions can be drawn from
these studies. This is, perhaps, partially due to the difficulties in obtaining reliable
measurements in multiphase systems. It also reflects the fact that most numerical studies used
the Eulerian/Lagrangian approach, which is limited to low-speed/dilute cases, while most
experiments are conducted under more realistic conditions of higher gas velocity and holdup.
With experimental techniques being developed and improved, it seems that the numerical study
is somewhat lagging. While most experiments are limited to laboratory scale, industrial needs
for reliable numerical simulations of large-scale columns are real. A computer code for such
simulations should be able to deal with the situations involving large gas velocity, high gas
holdup and churn turbulent flows. A key step towards this goal is to validate the numerical
results, in a quantitative way, by experimental measurements. As mentioned above, the studies
of Lin et al. (1996) and Mudde et al. (1997) provide reliable and extensive data on bubbly flow
between two narrowly separated plates. Their experiments can be approximately simulated by
solving a set of two-dimensional equations. By doing so, one should be able to utilize their
measurement results, both for time-averaged and for the transient properties, to validate the
numerical prediction. In addition to the mean values, the comparison of the transient properties
is of particular importance since it is related to the back mixing and turbulence in the column.
The testing of various physical models needed for closure can also be accomplished relatively
readily in two-dimensional simulations, in comparison with the three-dimensional simulations.

In this study we present an Eulerian/Eulerian dynamic simulation of a two-dimensional bubble

column. The ensemble-averaged equations are used to solve the velocity and volume fraction
field for both phases. A model of bubble-induced turbulent viscosity is incorporated into the
momentum equations for the liquid phase. The numerical predictions of the macroscopic
structures and mean properties are compared with the experimental data provided by Mudde et
al. (1997) and Lin et al. (1996). Our objective is to verify the applicability and accuracy of the
Eulerian/Eulerian method for the dynamic simulation of bubble-driven two-phase flow in two-
dimensional bubble columns. In this report some of the results are discussed.

27
2.1 Simulation Conditions

The simulations were performed for the 11.2-, 15.2- and 10.16-cm columns. The flow
conditions are listed in Table 2.1. These conditions match the experiments by Lin et al. (1996)
and Mudde et al. (1997). The runs for the 11- and 15-cm columns compare mean properties,
while the runs for the 10-cm column provide the characteristics of large structures and their
variations with superficial gas velocity.

Table 2.1 Column Size and Operation Conditions of the Cases Simulated
Column Width Superficial Gas Static Liquid Number of Gas Aspect Ratio
W (cm) Velocity Usup Height H (cm) Injectors
(cm/sec)
11.2 1.0 110 2 9.8
15.2 1.0 110 3 7.2
10.16 1.22 160 2 15.8
2.44 160 2 15.8
3.66 160 2 15.8
5.49 160 2 15.8

2.2 Governing Equations

In accordance with Drew (1983), we provide the following governing equations for the motion
of gas-liquid flow in bubble columns. In these equations, the subscripts ‘c’ and ‘d’ indicate the
properties for the continuous phase (liquid) and the dispersed phase (gas), respectively. The
expression for the drag coefficient for the inter-phase momentum exchange is adopted from
Katsumi et al. (1997) and Drew (1983). A model for the bubble-induced stress in the liquid
phase, as proposed by Sato et al. (1981), is applied in the simulations.

Equations of Continuity

∂ε c
Liquid + ∇ .(ε c u c ) = 0 (2.1)
∂t
∂ε d
Gas + ∇ .(ε d u d ) = 0 (2.2)
∂t

28
Momentum Equations

 ∂u 
Liquid ρ c ε c  c + u c .∇ u c  = ρ c ε c g − ε c ∇p − (M d + M vm ) + ∇.(ε c σ c ) + ∇ ⋅ (ε cσ cb ) (2.3)
 ∂t 

 ∂u 
Gas ρ d ε d  d + u d .∇ u d  = ρ d ε d g − ε d ∇p + (M d + M vm ) + ∇.(ε d σ d ) (2.4)
 ∂t 

Inter-Phase Momentum Exchange

6ε cε d
Md = Fd (2.5)
πd p3

1
Fd = ρ cπd p2 C D u c − u d (u c − u d ) (2.6)
8
 24
C D = max  (
1 + 0.15 Re 0.687 , f )
8 Eo 
3 Eo + 4 
(2.7)
 Re
2
 1 + 17.67ε c9 7 
f =  (2.8)
 18.67ε c
32

where
Eo ≡ Eotvos Number ≡ gρ c d p2 γ

and Re ≡ Bubble Re ynolds Number ≡ ε c d p u c − u d ν c

1  Du c Du d 
M vm = ε c ε d C vm  −  (2.9)
2  Dt Dt 

C vm = 1 + 3.32ε d + O( ε d2 ) (2.10)

Bubble-Induced Stress

σ cb = ρ cν bt (∇uc + ∇ucT ); ν bt = k b ε d d p uc − ud ; k b = 0.6 (2.11)

29
2.3 Mean Properties
Figure 2.1 compares the vertical and horizontal mean velocity profiles in the 11-cm column at
Usup=1 cm/s. These profiles relate to the middle section of the column, where the mean flow is
usually assumed to be one-dimensional. The numerical prediction of the mean horizontal
velocity, U , is essentially zero, as expected. However, the experimental time-averaged
horizontal velocity is non-zero and exhibits an inward flow. As pointed out by Mudde et al.
(1997), this is attributed to a systematic error due to the difficulty in tracking the particles in the
fast-moving bubble stream and other biases. The comparison of the computed mean vertical
velocity profile with data is good except for the near-wall region. One reason may be that the
flow close to the walls is actually three-dimensional due to the effect of finite thickness of the
columns, which is on the order of the distance to walls in the near wall region. The other reason
for this discrepancy is that the boundary layer is too thin to be resolved in our current
simulations. Notice that the curve of V starts and ends at the nearest points next to the walls on
which the value of V is set to zero; this is due to the no-slip boundary condition imposed on the
liquid phase during the simulation. The fact that V at starting and ending point reaches the
maximal negative value clearly indicates that the boundary layer was not resolved. By
examining the measured data points near the wall, one should realize that the experiment did not
resolve the wall layer either. Unfortunately, the gas holdup profile was not reported in the
experiments of Mudde et al. (1997); therefore, one cannot test to what extent the data satisfy the
continuity equation. The numerical results do satisfy the mass balance.

For the 15-cm column, the comparisons of the computed and experimentally determined mean
vertical liquid velocities are presented for the middle, lower and upper sections in Figure 2.2.
Again the numerical values match the data in magnitude except in the wall region. The
experimental profiles for the middle and upper sections are similar, while the one for the lower
section is somewhat different. Mudde et al. (1997) argued that this is due to the fact that the
flow is not fully developed in the lower section. On the other hand, the computed profiles for the
three sections are almost identical. Thus, the numerical values for the upper and middle sections
are higher in magnitude than the experimental values. The discrepancy between the predictions
and data in the near wall region may be attributed to the effect of finite thickness of the third
dimension. In experiments, the "two-dimensional" flow was realized between two narrowly
separated plates. The distance between the plates was 1.27 cm. In the near-wall region the flow
is really rather three-dimensional. Thus the slowing down of fluid by the wall in the third
dimension, i.e., the parallel plates, becomes significant. However, this effect is not accounted for
in the two-dimensional simulations. Since the thickness of the boundary layer is inversely
proportional to the velocity, we would expect this type of discrepancy to be reduced as velocity
increases. In fact the profiles for the lower section in Figure 2.2, where the velocity is high in
magnitude, exhibit this trend.

The numerical predictions of turbulence intensities, u 'u ' and v 'v ' , and Reynolds shear stress,
u 'v ' , for the 15-cm column operated at Usup= 1 cm/s are compared with data in Figure 2.3. The
u 'u ' reaches a peak at the central portion of the column since u attains its highest magnitude in
30
the center due to the nature of the meandering motion of the central plume. The peaks of v 'v ' at
the near-wall vortical region are consistent with the fact that the flow dynamically changes from
upward to downward in such areas. Although the general trends of the numerical values match
the data, there are significant differences in values. It is clear that the simulation under predicts
the turbulence-related properties. In reality, for the cases of low gas velocity, bubbles retain their
identity, while they rise up through the column. The drag force is thus always concentrated in
the liquid surrounding the bubble. However, in the Eulerian/Eulerian approach, bubbles are not
identified as single identities. Although the drag force is calculated based on the single bubble
formulation, it is not confined to the region adjacent to the bubbles, as it should be; rather it is
more "dispersed" in the liquid. (This is obviously the inherent drawback of using the two-fluid
model to describe two-phase flows at low gas holdup.)

2.4 Characteristics of Large Structures

The dynamic behavior of the large structures is characterized by the wavelength and frequency
of the meandering central plume. It is expected that these quantities vary with column size and
gas superficial velocity. Lin et al. (1996) have conducted an extensive and detailed experimental
investigation of this topic, in 2D columns. They found that, in the same column, the frequency
increases with gas velocity, while the wavelength decreases. At the same gas superficial
velocity, the frequency decreases as the size of the column increases. The wavelength is
proportional to the column size. When the wavelength and frequency are multiplied to give the
vortex descending velocity, they found that it is basically a function of the superficial gas
velocity only. They also found that the size of these vortices is independent of gas velocity when
Usup≤ 1 cm/s and varies with column size only.

To study the variation of large-scale structure, i.e., the frequency and wave length, with gas
velocity we conducted a group of simulations for the 10.16 cm column. The flow conditions, as
listed in Table 2.1, match the experiments done by Lin et al. (1996). Figure 2.4 shows the flow
patterns when the column is operated at four different superficial gas velocities. Consistent with
the experimental observation, wavelength decreases as the gas velocity increases. The overall
gas holdup increases as Usup increases. As the gas velocity increases, one can observe that the
turbulence is intensified and the flow structure becomes less clearly defined. Figure 2.5 shows
the time sequences of the liquid velocity component u at a central point for each case collected at
the sampling frequency of 10 Hz. The frequency clearly increases with the superficial gas
velocity. The flow becomes more chaotic as gas velocity increases, while the primary frequency
still remains distinguishable.

For a quantitative comparison we performed each simulation for 200 seconds after it reached the
quasi steady state and recorded all the quantities. By means of flow visualization and Fourier
analysis, the averaged frequency and wave length are calculated.

Figure 2.6 shows the comparisons of the frequency and wavelength of the large structure with
the experimentally measured values of Lin et al. (1996). The experimental values are understood
31
to have been obtained visually, as implicitly indicated in their paper. We denote the frequency
and wavelength of the larger structures as f0 and λ0, respectively. The computed values of λ0 and
f0 are measured by the visualizations and animations of the numerically generated liquid velocity
field. The other way to estimate the characteristics of the larger structures in a column is to
calculate the mean frequency and wave length, f and λ , as follows,

λ=
∑ λ P (λ )
i i u i
f =
∑ f S (f )
i i u i

∑ P (λ )
i u i ∑ S (f )
i u i

in which Pu (λ ) and Su ( f ) are the spatial and temporal Fourier spectra, respectively, for the
horizontal component of the liquid velocity along the central line of the column. λ and f are
then averaged over time and space, respectively. It can be seen that the numerical values match
fairly well with the measured data, particularly for the two cases of intermediate Usup. We have
noticed that at low and high gas velocity, the structures are not uniform over space and time, as
can be seen in Figures 2.4 and 2.5. Such non-uniformity may have induced the deviations in the
calculated and measured values.

The effect of the distributor, the sensitivity analysis and the effect of mesh size will be discussed
in the next monthly report.

References:
Drew, D. A., “Mathematical modeling of two-phase flow,” Ann. Rev. Fluid Mech., 15, 261
(1983).

Katsumi, T., A. Furumoto, L.-S. Fang and J. Zhang, "Suspension viscosity and bubble rise
velocity in liquid-solid fluidized beds," Chem. Eng. Sci., 52, 3053 (1997).

Lin, T.-J., J. Reese, T. Hong and L.-S. Fan, “Quantitative analysis and computation of two-
dimensional bubble columns," AIChE. J., 42, 301 (1996).

Mudde, R. F., D. J. Lee, J. Reese and L.-S. Fan, “Role of coherent structures on Reynolds
stresses in a 2D bubble column," AIChE. J., 43, 913 (1997).

Sato, Y., M. Sadatomi and K. Sekoguchi, "Momentum and heat transfer in two-phase bubble
flow I," Int. J. Multiphase Flow, 7, 167 (1981).

32
Figure 2.1 Time-Averaged Liquid Velocity Profiles for the 11-cm Column at Usup= 1 cm/s.
 and --------- are the Numerical Predictions of the Vertical (Axial) Velocity Component
V and the Horizontal Velocity Component U , Respectively. • and Υ represent the

Experimental Measurements of V and U , Respectively. The Experimental Data are from

Mudde et al. (1997).

33
 Figure 2.2 Time-Averaged Liquid Velocity Profiles, V , for the 15-cm Column at Usup = 1
cm/s. − and • represent the numerical prediction and experimental data (Mudde et al.
(1997), respectively.
34
 Figure 2.3 The Time-Averaged Profiles of Turbulence Intensities and Reynolds Shear Stress for the
Middle Section of the 15-cm Column at Usup= 1 cm/s. 
35
and • represent the numerical values and data,
respectively..
 Figure 2.4 The Instantaneous Contours of Gas Holdup for the 10 cm Column Operated at
Different Superficial Gas Velocity. From the left to the right, Usup =1.22, 2.44, 3.66 and 5.49 cm/s.

36
Figure 2.5 Time series of the velocity components at a central point of the 10-cm column at
different Usup. From the top to the bottom, Usup= 1.22, 2.44, 3.66 and 5.49 cm/s.

37
 Figure 2.6 The Variation of Frequency and Wave length of the Meandering Structure in a 10-
cm Column with Superficial Gas Velocity. The experimental data is taken from Lin et al.(1996).

38
OHIO STATE UNIVERSITY
The report from Ohio State for the period follows.

INTRINSIC FLOW BEHAVIOR IN A SLURRY BUBBLE COLUMN UNDER HIGH

PRESSURE AND HIGH TEMPERATURE CONDITIONS

Quarter Report

(Reporting Period: January 1 to March 31, 1999)

Highlights

• The maximum stable bubble size in high-pressure slurry bubble columns was studied
experimentally and analytically. The maximum bubble size decreases with an increase in
pressure, and increasing solids concentration leads to larger bubbles. The particle effect is
more significant at ambient pressure.

• A mechanistic model based on the concept of internal gas circulation inside a bubble was
proposed to simulate the maximum stable bubble size in a high-pressure slurry bubble
column. The mechanism and criterion of bubble breakup were illustrated by considering the
internal gas circulation inside the bubble.

• Based on the model, an analytical expression was obtained to predict the maximum stable
bubble size under high-pressure conditions. Comparison between the experimental data and
model predictions showed that the proposed model can reasonably predict the maximum
bubble size at elevated pressures.

Work Conducted

Experimental Study of Maximum Bubble Size

The bubble size and bubble size distribution in high-pressure slurry bubble columns can be
measured by using an optic fiber probe, which was described in previous monthly reports
(August - November, 1997). It should be noted that the probe can only measure the bubble
chord length, and not the true bubble diameter. If a uniform bubble shape (either spherical or
ellipsoidal) is assumed, the bubble chord length distribution can be converted to the bubble size
distribution (Liu et al., 1996). However, bubbles in slurry bubble columns have irregular shapes.
The data from the probe can thus only represent the chord length distribution.

Special attention should be placed on the maximum chord length for the dominant effect of large
bubbles on gas holdup. In addition, the measured maximum chord length can be interpreted as
39
the height of the largest bubble (hmax) in the system if the number of bubbles detected by the
probe is large enough. In general, the maximum bubble height differs from the maximum bubble
size. Nevertheless, large bubbles in 2D slurry bubble columns of different solids concentrations
have approximately the same maximum dimensions in the horizontal and vertical directions,
based on the photos shown by de Swart et al. (1996) and on observations from the 2D
experiments. The maximum bubble size can thus be approximated by the maximum bubble
height in 2D slurry bubble columns. It is further assumed that this approximation also holds for
3D columns; thus, the maximum bubble chord length is approximated as the maximum bubble
diameter in the following analysis.

The maximum bubble size under different conditions is shown in Figure 1. Under otherwise
constant conditions, the maximum bubble size decreases with an increase in pressure, especially
at pressures lower than 1.5 MPa. Increasing solids concentration leads to significantly larger
bubbles at ambient pressure over the entire gas velocity range. On the other hand, at the pressure
of 5.6 MPa, the solids concentration has a significant effect on the maximum bubble size only in
the gas velocity range of 8–23 cm/s. At gas velocities above that range, the maximum bubble
size is virtually independent of solids concentration.

Mechanism of Bubble Breakup

Internal Gas Circulation

An analytical criterion for the bubble breakup can be derived by considering a single large
bubble rising in a stagnant liquid at a velocity of ub, without any disturbances on the gas-liquid
interface. The bubble is subjected to breakup when its size exceeds the maximum stable bubble
size due to the circulation and centrifugal force. Large bubbles normally assume a spherical cap
shape; in this work, the spherical-cap bubble is approximated by an ellipsoidal bubble with the
same volume and the same aspect ratio (height to width). The bubble is described in a
cylindrical coordinate system that moves with the bubble and has the center of the ellipsoidal
bubble as its origin, as shown in Figure 2(a). The bubble surface is formed by rotating an ellipse
around the vertical axis, z:
rc2 z 2
+ = 1(1)
a 2 c2
rc2 z 2
+ =1 (1)
a 2 c2
where rc is the distance from the z axis. The aspect ratio of the ellipse (α = c/a) is the same as
that of the real spherical cap bubble. The internal circulation can be described by Hill’s vortex
(Hill, 1894) because of the high Reynolds number of the gas in the bubble. The flow field is
symmetrical about the z axis and has no azimuth component. The circulation velocity in the rc
and z directions are, respectively,
3u 3ub
u = b2 zrc = zrc (2a)
2c 2α 2 a 2
and
3ub  2 z2 
w = 2 (a − 2rc ) − 2  .
2
(2b)
2a  α 
40
In the presence of contaminants, a small bubble would behave like a rigid particle and thus, the
circulation of the gas inside the bubble could be suppressed. For large bubbles, as encountered
in the model, the high shear force generated by the high relative velocity between the bubble and
the liquid or slurry could readily sweep away the contaminants on the surface; thus, the
contaminants induce negligible effects on the internal circulation of the gas (Levich, 1962).

Evaluation of Centrifugal Force

To model the bubble breakup, it is necessary to evaluate the x-component of the centrifugal
force, Fx, on the entire bubble surface, as shown in Figure 2(c). Based on Eqs. (2a) and (2b),
there is a surface within the bubble, S, on which the z component of the circulation velocity, w, is
zero, and the circulation velocity has only the component in the rc direction, u, as shown in
Figures 2(a) and 2(b). This surface can be described as
2 rc2 z2
+ = 1. (3)
a2 α2a 2
Surface S passes through the vortex center and intercepts all vortex streamlines, as shown in
Figure 2(a). Because the pressure field inside the bubble is symmetrical about the x-y plane, Fx
is simply the rate of change of momentum across the surface S in the x direction, based on a
momentum balance for the gas phase. In addition, the streamlines inside the bubble are
symmetrical about the x-y plane and x-z plane, and thus, Fx is four times the rate of the x-
component of the gas momentum flowing across an octant of surface S shown in Figure 2(c),
( )
Fx = 4 ∫∫ ρ g u ⋅ δS u x = 4ρ g ∫∫ u 2 cos φ rc dφdz (4)

( )
S S

where ρ g u ⋅ δS is the mass flow rate across a surface element of S. Substituting Eq. (2a) into
Eq. (4) and integrating the resulted equation yields Fx:
π/2 c 2
 3u zr 
Fx = 4ρ g ∫ cos φ dφ ∫  b2 c2  rc dz
0 0 
2α a 
. (4a)
c
 3ub zrc 
2
9ρ g ub2 c 3 2
= 4ρ g ∫   rc dz = 4 4 ∫ rc z dz
0
2c 2  αa 0
From Eq. (3), rc can be expressed as
1 2  z   2

rc = a −   . (5)
2  α 
Substituting Eq.(5) into Eq. (4a) and integrating the resulting equation yields Fx:
3/ 2
9ρ g ub2 c
 2  z 2 
Fx = ∫ z a −  α  
2
dz
2 2α 4 a 4 0   . (6)
9π 0.312
= ρ g ub2 a 2 = ρ g ub2 a 2
64 2α α

Criterion of Bubble Breakup

The surface tension force is the product of the surface tension and the circumference of the
ellipse,
41
 Fσ = σL = σ ∫
ellipse
( δrc ) 2 + ( δz ) 2 = 4σa E( 1 − α 2 ) (7)

where E( 1 − α 2 ) is the complete second kind of elliptic integral that decreases with a decrease
in α. Also, the volume equivalent bubble diameter, de, is related to a and α by
π 3 4π 3
de = αa (8)
6 3
or
d
a= e . (8a)
3 8α
Note that the centrifugal force is affected significantly by the aspect ratio of the bubble, as well
as by the bubble size and the bubble rise velocity. The bubble is not stable if Fx is larger than the
surface tension force, Fσ, i.e.,
2
0.312  d   d 
ρ g u b2  e  ≥ 4σ  e  E( 1 − α 2 ) (9)
α 3
 8α 
3
 8α 
8α 4 / 3 E( 1 − α 2 ) σ
or u b2 d e ≥ . (9a)
0.312 ρg
When the centrifugal force is larger than the surface tension force, the bubble should be stretched
in the x direction. During the stretching, the aspect ratio, α, becomes smaller while de and ub can
be assumed to remain constant. As a result, the centrifugal force increases, the surface tension
force decreases, and the bubble stretching becomes an irreversible process. When the bubble is
stretched to an extent at which the necking phenomenon occurs, the bubble splits into two
smaller bubbles. A sequence of events in the photographs shown in Figure 3 confirms the
proposed mechanism of bubble breakup. The bubble images in the figure were obtained at a
pressure of 3.5 MPa. It can be seen that bubble stretching during rising takes place from t = 0 to
68 ms, and the necking starts at t = 68 ms. At t = 85 ms, the bubble splits into two smaller
bubbles.

Equation (9a) indicates that the breakup of bubbles is dictated by various factors, including
bubble size, bubble rise velocity, aspect ratio of the bubble, gas density, and gas-liquid surface
tension. In addition to the direct contribution to the centrifugal force as illustrated by Eq. (9a),
increasing the bubble size leads to a larger bubble rise velocity and a smaller aspect ratio of the
bubble, which both favor bubble breakup. As the pressure increases, the gas density increases
and the surface tension decreases; thus, bubbles are less stable.

Evaluation of Maximum Stable Bubble Size

The bubble rise velocity in high-pressure slurries can be calculated by a correlation developed by
Luo et al. (1997). The simplified form of that correlation for a large bubble is
2.8σ gd e ρ sl − ρ g
ub = + . (10)
ρ sl d e 2 ρ sl

42
Equation (10) has a form similar to the modified Mendelson equation (Mendelson, 1967) of
Maneri (1995), except that the slurry density here is replaced with the liquid density in Maneri’s
equation. By substituting Eq. (10) into Eq. (9a), de can be solved from the resulting equation:

de ≥ 
(
2σ  8α 4 / 3 E 1 − α 2) ρ ls
−
2.8 
≅
( )
16α 4 / 3σ E 1 − α 2
. (11)
g  0.312ρ g ρ sl − ρ g ρ sl − ρ g  0.312 gρ g
Therefore, the maximum stable bubble size is
σ
Dmax ≈ 7.16α 2 / 3 E( 1 − α 2 )1 / 2 . (12)
gρ g
The aspect ratio of bubble, α, is related to the wake angle, θw, by
1 − cos θ w
α= (13)
2 sin θ w
and the wake angle for bubbles in liquids depends on the bubble Reynolds number (Clift et al.,
1978),
(
θ w = 50 + 190 exp − 0.62 Re 0.4 . ) (14)
The wake angle of large bubbles in the liquids can be approximated as 50º, since the Reynolds
number is normally high. Hence, the aspect ratio is about 0.21 for the large bubbles in liquids,
and E( 1 − α 2 )1 / 2 is close to unity. For large bubbles rising in liquid-solid suspensions, the
aspect ratio is approximated as 0.3 (Fan and Tsuchiya, 1990) and E( 1 − α 2 )1 / 2 is 1.018. The
simplified forms of Eq. (12) are
σ
Dmax ≈ 2.53 (for α = 0.21) (15a)
gρ g
or
σ
Dmax ≈ 3.27 (for α = 0.3). (15b)
gρ g
The internal circulation model provides an estimation of the upper limit of the maximum stable
bubble size, since the external stresses in the liquid phase are ignored in the model. In actual
bubble columns or slurry bubble columns, the observed maximum bubble size should be smaller
than the predictions by Eqs. (15a) and (15b).

Prediction Results

Figure 4 compares the experimental data of maximum stable bubble sizes obtained in the present
experimental system (Nitrogen-Paratherm NF heat transfer fluid-100 µm glass beads) with the
predictions by different models. The experimental data of maximum stable sizes in Figure 4 are
approximated by the maximum bubble sizes measured over the entire gas velocity range under
the specified pressure, temperature, and solids concentration conditions. The comparison
between the predictions and the experimental data indicates that the internal circulation model
can reasonably predict the maximum stable bubble size at high pressures in the present system,
while the observed effect of pressure on bubble size cannot be explained by the Rayleigh-Taylor
instability or the Kelvin-Helmholtz instability. For an air/water system, the prediction by Eq.
43
(15a) for the maximum stable bubble size is 5.7 cm at 1.5 MPa, while Wilkinson (1991)
observed a maximum bubble size of about 3 cm. However, Eqs. (15a) and (15b) over-predict the
maximum stable size at ambient pressure, as shown in Figure 4. For an air/water system, the
predicted maximum stable bubble size is 22 cm, and the reported value is 13.2 cm. The over-
prediction by the internal circulation model at ambient pressure suggests that the bubbles are
disintegrated by disturbances in the liquid or slurry through other types of mechanisms, e.g.,
Rayleigh-Taylor instability (Grace et al., 1978) and Kelvin-Helmholtz instability (Kischa and
Kocamustafaogullari, 1989), before the bubble size reaches the value predicted by the internal
circulation model. Figure 4 shows that the maximum stable bubble size can be estimated by the
internal circulation model at pressures higher than 5 atm; at lower pressures, models based on
Rayleigh-Taylor instability or Kelvin-Helmholtz instability should be used. To improve the
prediction of the internal circulation model at low pressures, the physical properties of the liquid
or the slurry and the external stresses would need to be taken into account.

Based on the internal circulation model, the inertia of the gas and gas-liquid surface tension are
the dominant factors dictating the maximum stable bubble size at elevated pressures. The model
illustrates an insignificant effect of solids concentration on the maximum bubble size at high
pressures, since the stresses in the slurry are ignored in the model. However, at ambient
pressure, solids concentration has an indirect effect on the maximum bubble size through
variations of the shape or the aspect ratio. The maximum stable bubble size would be larger in
the slurry than in the liquid, as exhibited by Eqs. (15a) and (15b). Note that the aspect ratio in
the slurry (0.3) is larger than that in the liquid (0.21) for a given pressure at high bubble
Reynolds numbers.

Notations
a Width of ellipse
c Height of ellipse
de Volume equivalent bubble diameter
Dmax Maximum stable bubble size
E Complete second kind elliptic integral
Fx x-component of the centrifugal force generated by gas circulation
Fσ Surface tension force
g Gravitational acceleration
hmax Height of the largest bubble
L Circumference of the ellipse
P System pressure
rc Radius in a cylindrical coordinate system
Re Bubble Reynolds number, ρldeub/µl
S Surface within bubble on which the z-component of the circulation velocity is zero
t Time
u rc-component of the circulation velocity of gas inside a bubble
ub Bubble rise velocity
Ug Gas velocity
ux x-component of the circulation velocity of gas inside a bubble
w z-component of the circulation velocity of gas inside a bubble
z z-coordinate in a cylindrical coordinate system
44
Greek letters
α Aspect ratio of bubbles
φ Azimuthal angle in spherical coordinates or cylindrical coordinates
ρ Density
θw Wake angle
σ Gas-liquid surface tension
φs Solids concentration
µ Viscosity

Subscripts
g Gas phase
l Liquid phase
sl Slurry

References

Clift, R., J.R. Grace, and M.E. Weber, Bubbles, Drops, and Particles, Academic Press, New
York (1978).
de Swart, J.W.A., R.F. van Vliet, and R. Krishna, “Size, Structure and Dynamics of “Large”
Bubbles in a Two-Dimensional Slurry Bubble Column,” Chem. Eng. Sci., 51, 4619 (1996).
Fan, L.-S. and K. Tsuchiya, Bubble Wake Dynamics in Liquids and Liquid-Solid Suspensions,
Butterworth-Heinemann, Stoneham, MA (1990).
Grace, J.R., T. Wairegi, and J. Brophy, “Break-Up of Drops and Bubbles in Stagnant Media,”
Can. J. Chem. Eng., 56, 3 (1978).
Hill, M.J.M., “On a Spherical Vortex, ” Phil. Trans. Roy. Soc. London, 185, 213 (1894).
Kitscha, J. and G. Kocamustafaogullari, “Breakup Criteria for Fluid Particles,” Int. J. Multiphase
Flow, 15, 573 (1989).
Levich, V.G., Physiochemical Hydrodynamics, Prentice Hall, Englewood Cliffs, NJ (1962).
Liu W., N.N. Clark, and A.I. Karamavruc, “General Method for the Transformation of Chord-
Length Data to a Local Bubble-Size Distribution,” AIChE J., 42, 2713 (1996).
Luo, X., J. Zhang, K. Tsuchiya, and L.-S. Fan, “On the Rise Velocity of Bubbles in Liquid-Solid
Suspensions at Elevated Pressure and Temperature,” Chem. Eng. Sci., 52, 3693 (1997).
Maneri, C. C., “New Look at Wave Analogy for Prediction of Bubble Terminal Velocities,”
AIChE J., 41, 481 (1995).
Mendelson, H. D., “The Prediction of Bubble Terminal Velocities from Wave Theory,” AIChE
J., 13, 250 (1967).
Wilkinson, P.M., Physical Aspects and Scale-Up of High Pressure Bubble Columns, Ph.D.
thesis, University of Groningen, The Netherlands (1991).

45
 8.0

φs = 0

P = 0.1 MPa
6.0 φs = 0.081

P = 0.1 MPa
(cm)
(cm)

φs = 0.191
P = 5.6 MPa
4.0
d b, max
max
D

P = 0.1 MPa
P = 1.5 MPa
2.0
P = 5.6 MPa

0.0
0 10 20 30 40 50
Ug (cm/s)

Figure 1. Effects of Pressure and Solids Concentration on the Maximum Bubble Size in a Slurry Bubble Column (T = 78ºC)

46
 z
c y
ub Fσ
z
z ux
c u
Bubble surface
c o a Fx
v φ
w u
z S
-a u a rc -a
o rc rc rc
o Bubble surface
S a y
Vortex
S
o a x
streamline -c φφ

-c Bubble surface ux
S
u rc
-a
Bubble surface
(a) (b) (c)

Figure 2. Schematic of the Internal Circulation Model for Bubble Breakup: (a) Internal and External Flow Fields; (b)
Circulation Velocity on Surface S; (c) Force Balance and 3D View of Surface S and the Flow Pattern on S

47
 (a) t = 0 ms (b) t = 17 ms (c) t = 34 ms

(d) t = 51 ms (e) t = 68 ms (f) t = 85 ms

Figure 3. A Sequence of Bubble Images Showing the Process of Bubble Breakup

(P = 3.5 MPa)

48
 0.20
Solids fraction: 0.191, T=78ºC
Solids fraction: 0.081, T=78ºC

0.15 Rayleigh-Taylor instability

Kelvin-Helmhotz instability (α = 0.21)
(α = 0.30)
Dmax (m)

Internal circulation model

0.10 Internal circulation model

0.05

60
0.00
0 20 40 60
Pressure (atm)
Figure 4. Comparison of the Maximum Stable Bubble Size Between the
Experimental Data and the Predictions by Different Models

49