- Email: [email protected]
Hydrodynamics in bubble columns
101
Hydrodynamics
in Bubble Columns
Hydrodynamik von Blasenshlen H.-E. GASCHE,
Ch. EDINGER,
H. KGMPEL
and H. HOFMANN
Institut fiir Technische Chemie I, Universitiit Erlungen-Niirnberg, 8520 Erlangen (F.R.G.) Dedicated
to Prof. Dr.-lng.
(Received
March
Dr. h.c.llNPL
E. U. Schb%&r
Egerlandstr.
3,
on the occasion of his #h
birthaky
28, 1989; in final form May 25, 1989)
Abstract A heterogeneous fluid dynamic model has been developed to describe the complex flow structure of two-phase flow in bubble columns. The equation of continuity and the momentum balances are the basis of the model. The counling of the two phases is performed by a force of interaction which is deduced by-a force balance around by introducing the a single rising bubble. Multiphase flow mixing processes are taken into consideration turbulent viscosities of the two phases involved. The model equations were implemented successfully by applying a tridiagonal matrix algorithm.
Kurzfassung Zur Beschreibung der komplexen Hydrodynamik der Mehrphasenstrijmung in Blasensiiulen wurde ein heterogenes Model1 entwickelt, welches sowohl die Gas- als such die Fltissigkeitsstromung in diesem Reaktor beschreibt. Kontinuitiitsund Impulsbilanz sind die Basis flir das Striimungsmodell. Die Kopplung der beiden Phasen erfolgt tiber einen Wechselwirkungsterm, der sich aus der Krlftebilanz an einer Einzelblase ableitet. Das Verstandnis der Durchmischungsvorgiinge in mehrphasigen Striimungen wird durch die hier eingefiihrte Turbulenzmodellierung auf eine physikalisch besser interpretierbare Grundlage gestellt. Die Losung der Modellgleichungen erfolgt iiber einen Tridiagonal-Matrix-Algorithmus.
1. Introduction There is widespread use of gas-liquid reactors in the chemical industry. The potential range of applications is even expanding, as biotechnological processes are growing more and more important. The economic importance emphasizes the need for effective reactor models [ 11. However, in the field of bubble columns there is still a remarkable lack of physically based models which connect basic fluid dynamic principles with absorption and reaction problems, so that the influence of fluid dynamics on the design of the gas-liquid chemical reactor may be investigated in detail. Numerous problems in the prediction of bubble column performance are due to the fact that the construction of this reactor with no internals is very simple and therefore does not provide the operator with tools to control the fluid dynamic behaviour. But, as this is one of the key parameters in the performance of bubble columns, a detailed descrip0255-2701/89/%3.50
Gem.
Ena. Process..26(1989)
tion of the complex hydrodynamics in this reactor is important and should be based on physical principles controlling the turbulent mixing patterns in bubble columns. 1.1. Experimental
set-up
To approximate industrial practice, the reactor used for the experiments was constructed on the scale of a pilot plant. Figure 1 shows a flowsheet of the set-up. All data were acquired by computer. The experimental instrumentation comprised a conductivity microprobe for the fluid dynamic quantities of the gas phase and a hot-film anemometer to measure the liquid velocities. Measurements were performed at five different levels in the axial direction and 12-15 levels in the radial direction to obtain axial and radial profiles of the measured values. Semiconductor pressure transducers were used to measure the integral gas hold-up. 101-109
0 Ekvier
SeauoiaPrinted
in The Netherlands
.* .1 D.
,
-1 -2
Fig. 1. Experimental set-up (signal flowsheet). Column dimensions: height, 4.52 m; diameter, 0.29 m. Height above gas distributor: level 1, 0. IS m; level 2, 0.75 m; level 3, 1.50 m; level 4, 2.50 m; level 5, 3SOm.
-3 -‘
L -!A
-,*
I
-10
.
-8
-6
-4
-*
0
4
6
0
x1
Radial Position
1.2. Experimental
Summary of the experimental results From Fig. 2 it can clearly be seen that the axial bubble rise velocity develops very rapidly from the
p .s 3 _ .J.
.z 1
3
.2
.l
.O -11-1240-8
-6 -4 -2
0
2
6
4
Radial
8
10 12 U
Position
(cm)
Fig. 2. Axial development of the bubble rise velocity V,. (System: water/air; u, = 2.55 cm s-r; uoL = 0.00 cm s-r; 0, level 1; A, level 3; +, level 5.)
1
7t
-IL-12
-lO -8
-6
-4
-2
0 -
2
1*
u
km)
Fig. 4. Radial velocity distribution of the liquid phase. (System: water/air;u,,=lcms-‘;A,~=8cms-‘;0,aco=4cm~-’.)
results
The main objective of these investigations is the fluid dynamic description of bubble columns. Figures 2-4 show some typical results of the very extended experimental program, to indicate the type and precision of the measurements [2].
*
.A
*
2
I.
6
6
0
mdaol positim
12 14 Icm
I
Fig. 3. Radial velocity distribution of the gas phase. (System: water/air; T = 20 “C; u, = 6.7 cm s-r; countercurrent, 0, ucr_= -2cms-r; batch, A, u,,,=Ocms-‘; cocurrent, +, UoL=2cms-‘.)
bottom of the reactor up to level 3, but remains almost unchanged after crossing the distance between levels 3 and 5. Figures 3 and 4 also show the occurrence of recirculation caused by increasing superficial velocities [ 31. 1.3. Conclusions from the experimental
results
The rising velocities of the bubbles show a very strong radial dependence. From Fig. 3, a pronounced velocity decrease (even down to negative values) is immediately apparent close to the wall. The reason for this is due to the recirculating liquid which is flowing down in the peripheral region of the reactor (see Fig. 4). Therefore, the region of the steep velocity gradient, where velocities drop from + 50 to -20 cm s-i, must be a region of high turbulence. Miyauchi et al. were the first to establish a fluid dynamic model including recirculation of the liquid phase [4,5]. A further development of this model has been achieved recently, taking into account the slip velocities of the two phases [6]. However, up to now there exists no model that is able to describe the complex flow patterns of the gas and the liquid phase simultaneously in bubble columns. A well-known proof for this lack of knowledge is the upper part of Fig. 5. Different flow regimes are depicted in a flowsheet combining key design parameters to derive design rules for the chemical engineer. In contrast, the lower part of the Figure represents the result of a theoretical study performed with the hydrodynamic model that will be presented below. This study of the influence of the reactor diameter on the mean velocity profile reveals whether recirculation of the liquid phase can be expected or not for a given system, which should be the object of the upper part of Fig. 5, too. The basis of the model is the equation of continuity and the momentum balances. These equations connect global operation parameters with the local flow structure. The coupling of the two phases is performed by a force of interaction, which is deduced from a force balance around a single rising bubble. Physically, the mixing behaviour of the multiphase system is represented by the occurrence of turbulence.
103
25
(a)
5
10
Reactor
20
50
100
Diameter DR km)
tive use of this approach. Moreover, numerical problems (memory capacity of the computer) set their limits. Furthermore, from the experimental point of view Euler’s approach can be realized better than Lagrange’s, because it is easy to introduce into a reactor a stationary probe which measures continuously. Thus the present work utilizes Euler’s approach [8, 91. Both flowing phases are balanced heterogeneously, that is, each phase separately, but are considered to be quasi-continuous. The coupling of both phases is performed by a force of interaction which will be deduced subsequently. 2.1. Force balance on a rising bubble [8]
(b)
Relative
Radial Position (r/R)
Fig. 5. Flow
regimes in bubble columns: (a) flowsheet; (b) theoretical study (tloG = 7.2 cm SC’; q,,_ = 2.0 cm S- I., A , D,=O.lm; +, D,=0.3m; x, D,=0.6m).
Newton’s approach The resultant force on a bubble follows:
de
PC dt The model equations were implemented successfully in a computer simulation program by applying a finite-volume method. This method has many advantages because, by balancing across finite volumes, mass conservation solutions of the differential equations for each mesh point of the numerical grid as well as economic solutions for the whole of the grid can be achieved [7]. The results of the calculations are local values of the phase velocities which enable local rates for the chemical reaction to be determined. The connection of the hydrodynamic model with the mass and eventual energy balance of a given chemical reaction yields radial and axial profiles of reactant concentrations and thus provides the tool for a detailed simulation of bubble columns. Hence, the influence of hydrodynamics on the conversion and selectivity of complex gas-liquid reactions can be investigated and understood better. 2. Modelling
of the fluid dynamics
As the geometric shape of bubble columns is similar to a vertical cylinder, all the model equations are formulated in terms of cylindrical coordinates. Axial symmetry reduces the three spatial dimensions (r, z, cp) to two coordinates (r, z). The momentum balance in the q-direction is dropped by assuming that u, = 0 and d/dq = 0. The physical system to be treated with these equations is a two-phase flow system: gas distributed in a liquid. In principle, flowing fluids can be looked upon either from the Lagrange or from the Euler point of view. The dispersed two-phase system we deal with renders it possible to identify exactly defined fluid elements (e.g. bubbles), and thus seems to be an ideal case for Lagrange’s approach to the problem. However, as bubbles might contact each other and then lose their identity, the obstacle of coalescence hinders the effec-
=
-a%c&G B
-- dP dxi
(pressure
ULIbG
is formulated
- UL)
gradient
mass)
’
c1L o.5 d(u, - u,)/dz
(>I
(t
VL
(Basset forcz
force)
in flow direction)
(virtual +g!!5 4
(drag
as
dz
90.5
-
for non-stationary
processes)
(external forces, e.g. F, = pGg, gravitational force and F. = p,_g, buoyancy force)
+-Fe
Equation (1) means that the resultant on the bubble is composed of drag force + pressure
gradient
+ Basset force + external
+ virtual
(1)
force FrCS mass
forces
Now we regard the force balance in the main flow direction (xi = z). Our basic assumption is that we can neglect the virtual mass because of stationarity and slip free acceleration. Rearranging (1) gives F,, = -+c&G
-
uL\(uG
-
UL)
+
PLg
-
PGg
-
$
B
= -drag
force + buoyancy
- gravitational
force
force - pressure
gradient
(2)
Equation (2) incorporates the force a bubble is exposed to in a surrounding liquid under the above assumptions. For consistency, this force must disappear if one of the two interacting phases disappears. Equation (2) then becomes FI, k = f -
;F
B
&k(1 - &k)c&G - u&G dP
%pLg
+
&kPGg
+
x
1
- #,)
104 The sign of Fk is ‘minus’ for the gas phase (loses momentum) and ‘plus’ for the liquid (gains momentum). 2.2. Equation of continuity and momentum balances Equation
of continuity:
~(w,u,,~) Momentum
=0
balance
(4)
in z-direction:
(14)
* k312/1, = Ediss
(15)
* PLt=
(16)
cflPL(1eEdiss)1'3
(17) Assumptions (7)-( 17) lead to the final z-direction momentum balance to be implemented in a simulation program: *
+~$(ww,~)
Cpk’I&ctissP
* it = and
Pea
=
Pt
+
hn
+f$(~krpku~,ku,~k)
$[&k~k(u,,k)21
= -g +
ia
(Ek%,)+ ;
Ek + &kPk& +;
z
(Ekrd (5)
Fz,k
Momentum
balance
=--ap&k + i%
&kPkgr
in r-direction:
+
za
(& k t rz )
+;+k%)
(k = phase index; d/dt = 0 [9]). 2.3. Assumptions
&kPGg
(18) 1
ti; = mk,i, - rirz
The flow structure is assumed to be in the steady state and fully developed axially. The mean radial pressure gradient is negligible. The shear stresses are expressed as turbulent Reynolds shear stresses. In detail the assumptions can be described as follows. axially fully developed
flow
* u,=o
(7) (8)
Under the above assumptions nuity may be dropped. * g,=o * dP/dr
-
Solving the momentum balances might produce results which are contradictory to the continuity equation, because the axial pressure gradient used in the calculations is an estimated pressure gradient, as the true value is unknown. An algorithm has to be developed, therefore, to find the right pressure gradient to maintain continuity. As this is one of the crucial points to reach convergent behaviour of the system of equations a short survey of the solution strategy is given here [2, lo]. The fundamental idea is simple:
--kh+F,,k r
* d/dz = 0,
+ &kPLg
the equation
of conti-
(19) where the prime indicates a corrected value, and the asterisk an uncorrected value. In terms of velocities, R R U:, k, i&iri dr = tik, i,R2/2pk -
riU:
0
(20)
Equation (20) needs to be0 connected with the z-direction momentum balance in order to link pressure gradient and divergent mass fluxes. As the derivation of the appropriate expression is long and time consuming the deduction is omitted and only the final equations are presented [lo]:
(10)
+!&_d~!!!_~,,fkckrdr) 0
Jekr dr
(21)
0
1
dr
(9) = 0
These additional assumptions allow the r-direction momentum balance to be dropped too.
a,= ;
k, iEi
I
I
Numerous test runs of the simulation program proved that this correction of the pressure gradient works very satisfactorily. In order to accelerate convergence we used eqn. (21) to couple a subsequent correction of the velocities. u:, k is then evaluated as
i=j
i#j
zii is formulated in terms of Reynolds shear stresses; k represents the turbulent kinetic energy and is defined as: * k = (z&2+ uoz + u,2)/2
R
(13)
To evaluate the eddy viscosity w an approach similar to that of the b model (principle of eddy viscosity) is used:
c Ekr& d
AZ
R r u:* k6kr
2Pk J
0
dr
105 where (dP/dz)’ is the correction gradient according to the algorithm
* phase fractions are described by linear functions within the f&rite volume; * the force of interaction is considered to be constant within the finite volume. These assumptions give
of the pressure of eqn. (21).
2.4. Numerical procedure The discretization of the model equations is performed by a finite-volume method to guarantee mass conservation and economic solutions. We used an implicit tridiagonal matrix algorithm to solve the solution matrix with minimum computer time and memory requirements [lo, 111. The first step to develop the numerical solution procedure of the differential equations is to derive equivalent difference equations: for this the considered and balanced area is divided into a number of finite volumes which are centred around a central point (Fig. 6). The geometry of the grid is orthogonal and equidistant which means that all distances between any neighbouring grid points are equal in each spatial direction. Integration across one finite volume delivers the difference equations we are looking for. Optimum grid selection is conditioned by stability and economy. This point will be referred to in $3 [ lo].
Nj ----% --WPE-L !F - -- -- -- j-1
z
r
-
-_
5 __
__
i-l
i
i+l
‘w
----2
j+l
X (UC,
i -
uL, i)
-
&k. iPL,
u i, k, +
Ek. ErEk
-
&k, W rW k
ig
i + 1 -
k, E
+
&k, iPG,
%,
k, i
ig
II
ri Ar
Ar U I,
k, i
-Uz,k,i-I
k, W
(23)
Ar
Solving eqn. (23) in terms of u,, k, i and rearranging it according to the nomenclature of eqns. (25)+28) yields A P, k, i%,
k, i --Ae,k,iU=,k,i+l+AW,k,iU=,k,i-_++k,i
(24) Coefficients Ap, k, i (grid point P), Aw, k, i (west side), Ae, t, i (east side) and the source term Sk,, are calculated according to (24): A P, k, i AE.
(25)
--(Ek,ErE~~,ff,k,E+~k,WrW~,~.k,W)/Ar
k, i =
@k,
Erd+,
(26)
k, dlAr
S
AW,
k, i =
bk,
wrw&ff,
(27)
k, w)lAr
f,
Fig. 6. The numerical grid.
Discretization of the momentum balances The starting point of the discretization is eqn. (18). To demonstrate the procedure we introduce a section of the numerical grid (Fig. 6). Because of the assumption of an axially fully developed flow we do not need to integrate from south to north. For this special case we consider the radial direction only. The assumptions for integration are: * derivatives are approximated by central differences;
‘AP. Aw,
k, 2 k, 3
A E,
k. 2
A P,
k, 3
AE,
A W,
x (%, i -
AP,
k, i
AS,
ig
+
&k, iPG,
uk, 2
sk.2
uk, 3
S k, 3
uk, i
A P,k. A’
A E. k,
A W,
A P, k, NPM
k, NPM
&k, iPL,
ig
11 r, Ar
(28)
k. i
A W,k.N
-
This system of difference equations defined by (24) can be transformed into matrix shape (eqn. 29). This matrix is a so-called ‘tridiagonal matrix’. A possible and easy way to handle the solution of this matrix is by using the Gauss algorithm. Referring to the shape of this matrix this algorithm is called the ‘tridiagonal matrix algorithm’. Rearranging, we have eqn. (30).
k, 3
k, i
uL, i)
N
=
sk,
i
uk, N
Sk.
N
uk, NPhj
Sk,
NPM
(29)
106 4 P. k. 2
A E. k, 2 A E,k, 3
k, 3
A,.
A P,
k, i
A,.
The source term S&i is evaluated by sk,
i -
Aw,
k.
fUk,i -
(31)
I
The last line of matrix (30) enables us to calculate the axial velocity at the grid position NPM: uk, NPM
=
S’k. NPM
(32)
k, NPM
1-4~.
By introducing this value into the last-but-one line of the matrix system we obtain the algorithm to evaluate uk,,,,. Repeating this procedure we can determine the whole velocity vector according to *A, i =
SL,
i IAP~
k, i -
AE
k, i”k,
i+
I IAP,
k. i
S’k, 2
uk. 3
s;
3
sL,
i
=
4.1
k, i
A P,k, N
Sk, f =
uk, 2
A E, k, N
uk, N
S’k, N
A P,
uk. NPh
%.
NPN I
(30)
NPll I
Another crucial question related to the numerical procedure is whether this numerical solution of the convergent region is stable against disturbances generated by different initial or boundary conditions. Or, in other words, whether this solution is a unique solution or just one of several possible solutions. To elucidate this problem the initial velocity profiles were varied in several runs from constant to parabolic and finally to turbulent velocity profiles to study their influence on the convergent solutions. The results of this study are depicted in Figs. 8 and 9.
(33)
I
1
5
Boundary conditions Each phase velocity has to fulfil the no-slip condition at the wall of the reactor, so the boundary condition at the wall is: reactor wall (r = I?):
%, k, (r
= R)
-0 -
(34)
On the axis of the reactor (r = 0) symmetry is assumed: reactor axis (r = 0):
(du, k /dr) ((I= 0) = 0
(35)
0.F
‘..d
0
,......
1000
3000
LOO0
Number of Iteration
3. Calculated results
_I
.,I.....
2000
Steps
Fig. 7. Cbnvergence of the numerical procedure: effect of iteration steps on the velocity profile. (Da = 0.15 m; u, = 7.3 cm SC’;
3.1. Checks for plausibility Each model has to be checked to see whether the calculated results are reliable and reasonable within the given limits of the restricting assumptions indispensable in the course of model derivation. First of all, numerical artefacts must be .avoided. The presented solutions of the model equations should be stable and convergent. In this context the question of the minimum number of iteration steps arises. For the critical judgement of this problem we follow the iterative development of a selected value of the calculated solution. Of course, this should be that point of the numerical grid which gets the starting information last. In our case, as we start to integrate from the reactor wall, this point is represented by the axial velocity on the axis of the reactor. Figure 7 shows how the solution runs into the convergence zone after about 2000 iteration steps.
a
3
2
3
.L
.5
Relative
5
.7
.t3
.9
1.
Radial Position (r/R)
Fig. 8. Influence of different initial profiles on u,,~(T). (DR = 0.29 m; u, = 7.2cm SK’; q,L= Z.Ocms-‘; the points for u,, = constant, turbulent and parabolic lie on the same curve.)
107 ‘;i
.3E
2
.3
2 3 .2
.l
0.
-.l 1
2
.3
.L
5
6
Relative
Radial
7
8 Position
9
1
(r/R)
Fig. 9. Influence of different initial profiles on u,, L(r). (DR = 0.29 m; u, = 7.2 cm s-l; u,,, = 2.Ocm s-‘; the points for u,, = constant, turbulent and parabolic lie on the same curve.)
As these two Figures clearly prove, the calculated profiles are not affected at all by different initial conditions. This is very satisfying from a practical point of view, because this means, in reverse, that we do not need precise foreknowledge to start the calculations.
.6
e
.5
cl
.L
a
r* “j 3
.3 .2 .l
0. 1 -2 .l
2
3
.L
.5
Relative Fig.
10.
0,
u=,~(T);
experiment 2.0 cm s-l.)
0
[12, 131.
.l
.2
.3
.b Radial
7
.FJ
Position
9
1
(r/R)
A, u=,~(T): comparison of calculation with (D,=O.lSm; u,=~.~c~Is-‘; u,,,_=
.L Relative
.5
.6 Radial
.7
.0
Position
.9
In the next Figures a comparison of calculated and measured velocity profiles is shown. Figure 10 shows a comparison of our calculations with experimental results found in the literature [ 12, 131. Figure 11 represents some actual investigations of this paper. Both Figures clearly demonstrate the fundamental suitability of the model for multiphase flow systems. Both gas and liquid velocity proties represent fairly well the basic features of the measured velocity profiles. The radial locations of maximum, minimum and turning point of the calculated curves are found at the right positions. The good agreement between calculation and experiment is emphasized by taking into account the fact that no mathematical routine for the optimization of model parameters has been applied yet. 3.2. Parameter
studies
Once a mathematical model is implemented on a computer it is a mighty tool for the simulation and prediction of the behaviour of the reactor. To yield an optimum design of a bubble column a simulation program helps to find out the sensitive parameters which control and dominate the fluid dynamic behaviour. Figures 12 and 13 depict the response of the model on variation of the superficial gas velocity. They prove the applicability of the model within a broad range of operation modes as the trends are represented in an appropriate manner. The slight decrease of the central curvature with high superficial gas velocities is an indication of the influence of the bubble diameter on the phase interaction force via the drag coefficient. For this reason a more detailed investigation of the effect of the bubble diameter was carried out (Figs. 14 and 15). The analysis of the momentum balance reveals the drag force to be dependent on the bubble diameter. Figure 14 shows the effect of different bubble sizes on the drag force. The bigger the bubbles, the smaller the contribution to the force balance. Small bubbles suffer higher loss of momentum and thus decreasing bubble diameters are correlated with
1.
(r/R)
0.
3
.2
.3
.I
.5
Relative Fig. 11. A, ~,,~(r); x, ZJ,,L(r): comparison of calculation and experiment (21. (DR = 0.29 m; conditions for u,, o measurements: z&=4.ocms-‘, uoL = 0.0 cm SC’; conditions for U, L measurements: u, = 4.0 cm s-l, u,, = 1.0 cm s-l.)
6 Radial
3
.8
9
1.
Position (r/R)
Fig. 12. u,,o(r) as a function of supertkial gas &=0.29m; %,=O.Ocms-‘; A, u,=2.0cmP’; &,o=4.0cms-‘; 0, ~~=66.0cms-‘.)
velocity. +,
0.
.1
2
.3
.I
.5
Relative
.6
3
Radial
.6
9
Position
1.
Fig. 13. u,, L(r) as a function of superficial (D,=0.29m; u,,,=O.Ocms-‘; A, ~~=2.Ocms-‘; 4.0 cm s-‘; u, = 6.0 cm SC’.)
a
.l
.2
.3
1
.5
Relative
.6
3
Radial
B Position
0.
1
.2
.3
k/RI
9
.L
.5
Relative
gas
velocity. +, u,=
1.
(r/R)
Fig. 14. Effect of bubble diameter on the drag force (FG/sL). (D,=O.lSm; u,=7.2cms-‘; u0,=2.0cm-‘; A, D,= 0.25cm; +,D,=0.5cm; x,D,=l.Ocm; O,D,=2.Ocm.)
.6 Radial
.7
.8
Position
9
1.
(r/R)
Fig. 16. Effect of different bubble size distributions on velocity profiles. (Da = 0.15 m; u, = 7.2 cm SC’; u,,~ = 2.0 cm S-'; u,, G: D, = 0.55 cm; A, constant; +, parabolic; u,, L: D, = 0.55 cm; x , constant; 0 parabolic.)
Another Figure may emphasize the crucial role of the bubble diameter. Up to this point of the calculations the radial distribution of the bubble diameter was assumed to be constant. Measurements of the bubble size distribution and a theoretical study of the effect of different kinds (e.g. parabolic) of radial bubble size distributions on the fluid dynamics are depicted in Fig. 16. The gas-phase velocity profile is not very sensitive to the type of radial size distribution, but the liquid-phase velocity profile is affected significantly. We observe a more pronounced profile by keeping the bubble diameters constant. In the case of a parabolic distribution we find very small bubbles in the peripheral region of the reactor which hinder the backflowing liquid in a more effective way than big bubbles would do, as can be deduced again from Fig. 14. In order to maintain continuity over a cross-section of the reactor the total sum of the liquid velocity profile is shifted to lower values.
4. Conclusions
-2
2
Relative
Radial
Position
(r/RI
Fig.
15. Effect of bubble diameter on the liquid velocity profile, u,=7.2cms-‘; uoL=2.0cm-‘; A, ur,L. (D,=O.l5m; D,=O.25cm; +,DB=0.5cm; x,D,=l.Ocm; 0,&,=2.O~m.)
decreasing gas-phase velocities. The same proceedings can be observed in the central region of the liquid velocity profile. For bubble diameters smaller than 2.5 mm the curves of Fig. 15 become physically meaningless, because the drag force dominates all other source terms of the momentum balance.
The present work demonstrates a physical approach to multiphase flow systems. The heterogeneous model is able to predict the occurrence of recirculation and therefore is a more precise tool than flowcharts. The agreement of experiment and theory is fairly good. Parameter studies of the contributions to the momentum balance make clear the role of the decisive parameters, such as bubble diameters. With regard to the reciprocal effect of fluid dynamics on absorption/chemical reactions this approach seems to be very promising. At present it is being extended to ‘hot systems’, including homogeneous as well as heterogeneous chemical reactions.
Acknowledgement The authors gratefully acknowledge the financial support of the Deutsche Forschungsgemeinschaft for this work.
109
Nomenclature
4 -4 Aw CD
z$n DR g
F 4
m m* m’ P x S t uOk %,
k
%,
k
V,X xi Z
6 &k
coefficient for grid point P, kg m-’ s-’ coefficient for grid point E, kg m-’ s-’ coefficient for grid point W, kg m-’ s-’ drag coefficient eddy viscosity factor bubble diameter, m reactor diameter, m acceleration due to gravity, m se2 reactor height, m turbulent kinetic energy, mz se2 characteristic length scale, m superficial mass flux, kg m--2 s-l estimated superficial mass flux, kg m -2 s ’ corrected superficial mass flux, kg mm2 s-’ pressure, N m-’ radial coordinate, m reactor radius, m source term, kg m-’ SC’ time, s superficial phase velocity, m s- 1 mean axial phase velocity, m S-I mean radial phase velocity, m SC’ mean bubble rise velocity, m s-l general space coordinate, m axial coordinate, m Kronecker’s delta phase fraction of phase k dissipation rate, m2 ss2 viscosity, kg m- 1 s- ’ kinematic viscosity, m2 s- ’ density, kg rn-’ shear stress, N me2 tangential direction (cylindrical rad
Subscripts B E eff G i
bubble east side of finite volume effective gas radial grid position
coordinates),
i k L lam N N NPM P r S I, Z v
axial grid position phase index (gas or liquid) liquid laminar north side of finite volume number of grid points maximum number of grid points grid point (position i, j) radial south side of finite volume turbulent west side of finite volume axial tangential direction
( = N + 1)
References 1 W.-D. Deckwer, Reaktionstechnik in Blasenstiulen, Salle + Sauerliinder, Aarau, 1985. 2 H.-E. Gasche, DissPrtation,UniversitPt Erlangen-Ntirnberg, 1989. 3 H. Hofmann and H.-E. Gasche, Investigations of the fluid dynamics of bubble columns, Proc. German/Japanese Symp. on Bubble Columns, 1988, Schwerte, F.R.G. 4 T. Miyauchi and C. N. Shyu, Kagaku Kogaku, 34 (1970) 985. 5 K. Ueyama and T. Miyauchi, AZChE J., 25 (1979) 258. 6 Th. Menzel and U. Gnken, Deterministische und stochastische Komponenten der Fliissigkeitsstr6mung in Gas/Fltissigkeitsreaktoren, Proc. 2% Tutzing Symp: Znstationdre Proresse und Mejlmethoden in Katalyse, Etektrochemie und Korrosion, 1989. 7 G. Scheuerer, Vorlesungsskript: Numerische LGsung der strdmungsmechanischen Grundgleichungen, Lehrstuhl fur Striimungsmechanik, Univ. Erlangen-Nihnberg. 1986. 8 M. Sommerfeld, Vorlesungsskript: Mehrphasenstrtimungen, Lehrstuhl fiir Stromungsmechanik, Univ. Erlangen-Numb&g, 1987. 9 R. B. Bird, W. E. Stewart and E. N. Lightfoot, Transport Phenomena, Wiley, New York, 1960. 10 H. Kompel, Diplomarbeit, Inst. Tech. Chem. I, Univ. Erlangen-Niirnberg, 1989. 11 C. Edinger, Diplomarbeit, Inst. Tech. Chem. I, Univ. ErlangenNiimberg, 1988. 12 J. Steinemann, Dissertation, Univ. Dortmund, 1985. 13 K. Franz, Dissertation, Univ. Dortmund, 1983.
Hydrodynamics
in Bubble Columns
Hydrodynamik von Blasenshlen H.-E. GASCHE,
Ch. EDINGER,
H. KGMPEL
and H. HOFMANN
Institut fiir Technische Chemie I, Universitiit Erlungen-Niirnberg, 8520 Erlangen (F.R.G.) Dedicated
to Prof. Dr.-lng.
(Received
March
Dr. h.c.llNPL
E. U. Schb%&r
Egerlandstr.
3,
on the occasion of his #h
birthaky
28, 1989; in final form May 25, 1989)
Abstract A heterogeneous fluid dynamic model has been developed to describe the complex flow structure of two-phase flow in bubble columns. The equation of continuity and the momentum balances are the basis of the model. The counling of the two phases is performed by a force of interaction which is deduced by-a force balance around by introducing the a single rising bubble. Multiphase flow mixing processes are taken into consideration turbulent viscosities of the two phases involved. The model equations were implemented successfully by applying a tridiagonal matrix algorithm.
Kurzfassung Zur Beschreibung der komplexen Hydrodynamik der Mehrphasenstrijmung in Blasensiiulen wurde ein heterogenes Model1 entwickelt, welches sowohl die Gas- als such die Fltissigkeitsstromung in diesem Reaktor beschreibt. Kontinuitiitsund Impulsbilanz sind die Basis flir das Striimungsmodell. Die Kopplung der beiden Phasen erfolgt tiber einen Wechselwirkungsterm, der sich aus der Krlftebilanz an einer Einzelblase ableitet. Das Verstandnis der Durchmischungsvorgiinge in mehrphasigen Striimungen wird durch die hier eingefiihrte Turbulenzmodellierung auf eine physikalisch besser interpretierbare Grundlage gestellt. Die Losung der Modellgleichungen erfolgt iiber einen Tridiagonal-Matrix-Algorithmus.
1. Introduction There is widespread use of gas-liquid reactors in the chemical industry. The potential range of applications is even expanding, as biotechnological processes are growing more and more important. The economic importance emphasizes the need for effective reactor models [ 11. However, in the field of bubble columns there is still a remarkable lack of physically based models which connect basic fluid dynamic principles with absorption and reaction problems, so that the influence of fluid dynamics on the design of the gas-liquid chemical reactor may be investigated in detail. Numerous problems in the prediction of bubble column performance are due to the fact that the construction of this reactor with no internals is very simple and therefore does not provide the operator with tools to control the fluid dynamic behaviour. But, as this is one of the key parameters in the performance of bubble columns, a detailed descrip0255-2701/89/%3.50
Gem.
Ena. Process..26(1989)
tion of the complex hydrodynamics in this reactor is important and should be based on physical principles controlling the turbulent mixing patterns in bubble columns. 1.1. Experimental
set-up
To approximate industrial practice, the reactor used for the experiments was constructed on the scale of a pilot plant. Figure 1 shows a flowsheet of the set-up. All data were acquired by computer. The experimental instrumentation comprised a conductivity microprobe for the fluid dynamic quantities of the gas phase and a hot-film anemometer to measure the liquid velocities. Measurements were performed at five different levels in the axial direction and 12-15 levels in the radial direction to obtain axial and radial profiles of the measured values. Semiconductor pressure transducers were used to measure the integral gas hold-up. 101-109
0 Ekvier
SeauoiaPrinted
in The Netherlands
.* .1 D.
,
-1 -2
Fig. 1. Experimental set-up (signal flowsheet). Column dimensions: height, 4.52 m; diameter, 0.29 m. Height above gas distributor: level 1, 0. IS m; level 2, 0.75 m; level 3, 1.50 m; level 4, 2.50 m; level 5, 3SOm.
-3 -‘
L -!A
-,*
I
-10
.
-8
-6
-4
-*
0
4
6
0
x1
Radial Position
1.2. Experimental
Summary of the experimental results From Fig. 2 it can clearly be seen that the axial bubble rise velocity develops very rapidly from the
p .s 3 _ .J.
.z 1
3
.2
.l
.O -11-1240-8
-6 -4 -2
0
2
6
4
Radial
8
10 12 U
Position
(cm)
Fig. 2. Axial development of the bubble rise velocity V,. (System: water/air; u, = 2.55 cm s-r; uoL = 0.00 cm s-r; 0, level 1; A, level 3; +, level 5.)
1
7t
-IL-12
-lO -8
-6
-4
-2
0 -
2
1*
u
km)
Fig. 4. Radial velocity distribution of the liquid phase. (System: water/air;u,,=lcms-‘;A,~=8cms-‘;0,aco=4cm~-’.)
results
The main objective of these investigations is the fluid dynamic description of bubble columns. Figures 2-4 show some typical results of the very extended experimental program, to indicate the type and precision of the measurements [2].
*
.A
*
2
I.
6
6
0
mdaol positim
12 14 Icm
I
Fig. 3. Radial velocity distribution of the gas phase. (System: water/air; T = 20 “C; u, = 6.7 cm s-r; countercurrent, 0, ucr_= -2cms-r; batch, A, u,,,=Ocms-‘; cocurrent, +, UoL=2cms-‘.)
bottom of the reactor up to level 3, but remains almost unchanged after crossing the distance between levels 3 and 5. Figures 3 and 4 also show the occurrence of recirculation caused by increasing superficial velocities [ 31. 1.3. Conclusions from the experimental
results
The rising velocities of the bubbles show a very strong radial dependence. From Fig. 3, a pronounced velocity decrease (even down to negative values) is immediately apparent close to the wall. The reason for this is due to the recirculating liquid which is flowing down in the peripheral region of the reactor (see Fig. 4). Therefore, the region of the steep velocity gradient, where velocities drop from + 50 to -20 cm s-i, must be a region of high turbulence. Miyauchi et al. were the first to establish a fluid dynamic model including recirculation of the liquid phase [4,5]. A further development of this model has been achieved recently, taking into account the slip velocities of the two phases [6]. However, up to now there exists no model that is able to describe the complex flow patterns of the gas and the liquid phase simultaneously in bubble columns. A well-known proof for this lack of knowledge is the upper part of Fig. 5. Different flow regimes are depicted in a flowsheet combining key design parameters to derive design rules for the chemical engineer. In contrast, the lower part of the Figure represents the result of a theoretical study performed with the hydrodynamic model that will be presented below. This study of the influence of the reactor diameter on the mean velocity profile reveals whether recirculation of the liquid phase can be expected or not for a given system, which should be the object of the upper part of Fig. 5, too. The basis of the model is the equation of continuity and the momentum balances. These equations connect global operation parameters with the local flow structure. The coupling of the two phases is performed by a force of interaction, which is deduced from a force balance around a single rising bubble. Physically, the mixing behaviour of the multiphase system is represented by the occurrence of turbulence.
103
25
(a)
5
10
Reactor
20
50
100
Diameter DR km)
tive use of this approach. Moreover, numerical problems (memory capacity of the computer) set their limits. Furthermore, from the experimental point of view Euler’s approach can be realized better than Lagrange’s, because it is easy to introduce into a reactor a stationary probe which measures continuously. Thus the present work utilizes Euler’s approach [8, 91. Both flowing phases are balanced heterogeneously, that is, each phase separately, but are considered to be quasi-continuous. The coupling of both phases is performed by a force of interaction which will be deduced subsequently. 2.1. Force balance on a rising bubble [8]
(b)
Relative
Radial Position (r/R)
Fig. 5. Flow
regimes in bubble columns: (a) flowsheet; (b) theoretical study (tloG = 7.2 cm SC’; q,,_ = 2.0 cm S- I., A , D,=O.lm; +, D,=0.3m; x, D,=0.6m).
Newton’s approach The resultant force on a bubble follows:
de
PC dt The model equations were implemented successfully in a computer simulation program by applying a finite-volume method. This method has many advantages because, by balancing across finite volumes, mass conservation solutions of the differential equations for each mesh point of the numerical grid as well as economic solutions for the whole of the grid can be achieved [7]. The results of the calculations are local values of the phase velocities which enable local rates for the chemical reaction to be determined. The connection of the hydrodynamic model with the mass and eventual energy balance of a given chemical reaction yields radial and axial profiles of reactant concentrations and thus provides the tool for a detailed simulation of bubble columns. Hence, the influence of hydrodynamics on the conversion and selectivity of complex gas-liquid reactions can be investigated and understood better. 2. Modelling
of the fluid dynamics
As the geometric shape of bubble columns is similar to a vertical cylinder, all the model equations are formulated in terms of cylindrical coordinates. Axial symmetry reduces the three spatial dimensions (r, z, cp) to two coordinates (r, z). The momentum balance in the q-direction is dropped by assuming that u, = 0 and d/dq = 0. The physical system to be treated with these equations is a two-phase flow system: gas distributed in a liquid. In principle, flowing fluids can be looked upon either from the Lagrange or from the Euler point of view. The dispersed two-phase system we deal with renders it possible to identify exactly defined fluid elements (e.g. bubbles), and thus seems to be an ideal case for Lagrange’s approach to the problem. However, as bubbles might contact each other and then lose their identity, the obstacle of coalescence hinders the effec-
=
-a%c&G B
-- dP dxi
(pressure
ULIbG
is formulated
- UL)
gradient
mass)
’
c1L o.5 d(u, - u,)/dz
(>I
(t
VL
(Basset forcz
force)
in flow direction)
(virtual +g!!5 4
(drag
as
dz
90.5
-
for non-stationary
processes)
(external forces, e.g. F, = pGg, gravitational force and F. = p,_g, buoyancy force)
+-Fe
Equation (1) means that the resultant on the bubble is composed of drag force + pressure
gradient
+ Basset force + external
+ virtual
(1)
force FrCS mass
forces
Now we regard the force balance in the main flow direction (xi = z). Our basic assumption is that we can neglect the virtual mass because of stationarity and slip free acceleration. Rearranging (1) gives F,, = -+c&G
-
uL\(uG
-
UL)
+
PLg
-
PGg
-
$
B
= -drag
force + buoyancy
- gravitational
force
force - pressure
gradient
(2)
Equation (2) incorporates the force a bubble is exposed to in a surrounding liquid under the above assumptions. For consistency, this force must disappear if one of the two interacting phases disappears. Equation (2) then becomes FI, k = f -
;F
B
&k(1 - &k)c&G - u&G dP
%pLg
+
&kPGg
+
x
1
- #,)
104 The sign of Fk is ‘minus’ for the gas phase (loses momentum) and ‘plus’ for the liquid (gains momentum). 2.2. Equation of continuity and momentum balances Equation
of continuity:
~(w,u,,~) Momentum
=0
balance
(4)
in z-direction:
(14)
* k312/1, = Ediss
(15)
* PLt=
(16)
cflPL(1eEdiss)1'3
(17) Assumptions (7)-( 17) lead to the final z-direction momentum balance to be implemented in a simulation program: *
+~$(ww,~)
Cpk’I&ctissP
* it = and
Pea
=
Pt
+
hn
+f$(~krpku~,ku,~k)
$[&k~k(u,,k)21
= -g +
ia
(Ek%,)+ ;
Ek + &kPk& +;
z
(Ekrd (5)
Fz,k
Momentum
balance
=--ap&k + i%
&kPkgr
in r-direction:
+
za
(& k t rz )
+;+k%)
(k = phase index; d/dt = 0 [9]). 2.3. Assumptions
&kPGg
(18) 1
ti; = mk,i, - rirz
The flow structure is assumed to be in the steady state and fully developed axially. The mean radial pressure gradient is negligible. The shear stresses are expressed as turbulent Reynolds shear stresses. In detail the assumptions can be described as follows. axially fully developed
flow
* u,=o
(7) (8)
Under the above assumptions nuity may be dropped. * g,=o * dP/dr
-
Solving the momentum balances might produce results which are contradictory to the continuity equation, because the axial pressure gradient used in the calculations is an estimated pressure gradient, as the true value is unknown. An algorithm has to be developed, therefore, to find the right pressure gradient to maintain continuity. As this is one of the crucial points to reach convergent behaviour of the system of equations a short survey of the solution strategy is given here [2, lo]. The fundamental idea is simple:
--kh+F,,k r
* d/dz = 0,
+ &kPLg
the equation
of conti-
(19) where the prime indicates a corrected value, and the asterisk an uncorrected value. In terms of velocities, R R U:, k, i&iri dr = tik, i,R2/2pk -
riU:
0
(20)
Equation (20) needs to be0 connected with the z-direction momentum balance in order to link pressure gradient and divergent mass fluxes. As the derivation of the appropriate expression is long and time consuming the deduction is omitted and only the final equations are presented [lo]:
(10)
+!&_d~!!!_~,,fkckrdr) 0
Jekr dr
(21)
0
1
dr
(9) = 0
These additional assumptions allow the r-direction momentum balance to be dropped too.
a,= ;
k, iEi
I
I
Numerous test runs of the simulation program proved that this correction of the pressure gradient works very satisfactorily. In order to accelerate convergence we used eqn. (21) to couple a subsequent correction of the velocities. u:, k is then evaluated as
i=j
i#j
zii is formulated in terms of Reynolds shear stresses; k represents the turbulent kinetic energy and is defined as: * k = (z&2+ uoz + u,2)/2
R
(13)
To evaluate the eddy viscosity w an approach similar to that of the b model (principle of eddy viscosity) is used:
c Ekr& d
AZ
R r u:* k6kr
2Pk J
0
dr
105 where (dP/dz)’ is the correction gradient according to the algorithm
* phase fractions are described by linear functions within the f&rite volume; * the force of interaction is considered to be constant within the finite volume. These assumptions give
of the pressure of eqn. (21).
2.4. Numerical procedure The discretization of the model equations is performed by a finite-volume method to guarantee mass conservation and economic solutions. We used an implicit tridiagonal matrix algorithm to solve the solution matrix with minimum computer time and memory requirements [lo, 111. The first step to develop the numerical solution procedure of the differential equations is to derive equivalent difference equations: for this the considered and balanced area is divided into a number of finite volumes which are centred around a central point (Fig. 6). The geometry of the grid is orthogonal and equidistant which means that all distances between any neighbouring grid points are equal in each spatial direction. Integration across one finite volume delivers the difference equations we are looking for. Optimum grid selection is conditioned by stability and economy. This point will be referred to in $3 [ lo].
Nj ----% --WPE-L !F - -- -- -- j-1
z
r
-
-_
5 __
__
i-l
i
i+l
‘w
----2
j+l
X (UC,
i -
uL, i)
-
&k. iPL,
u i, k, +
Ek. ErEk
-
&k, W rW k
ig
i + 1 -
k, E
+
&k, iPG,
%,
k, i
ig
II
ri Ar
Ar U I,
k, i
-Uz,k,i-I
k, W
(23)
Ar
Solving eqn. (23) in terms of u,, k, i and rearranging it according to the nomenclature of eqns. (25)+28) yields A P, k, i%,
k, i --Ae,k,iU=,k,i+l+AW,k,iU=,k,i-_++k,i
(24) Coefficients Ap, k, i (grid point P), Aw, k, i (west side), Ae, t, i (east side) and the source term Sk,, are calculated according to (24): A P, k, i AE.
(25)
--(Ek,ErE~~,ff,k,E+~k,WrW~,~.k,W)/Ar
k, i =
@k,
Erd+,
(26)
k, dlAr
S
AW,
k, i =
bk,
wrw&ff,
(27)
k, w)lAr
f,
Fig. 6. The numerical grid.
Discretization of the momentum balances The starting point of the discretization is eqn. (18). To demonstrate the procedure we introduce a section of the numerical grid (Fig. 6). Because of the assumption of an axially fully developed flow we do not need to integrate from south to north. For this special case we consider the radial direction only. The assumptions for integration are: * derivatives are approximated by central differences;
‘AP. Aw,
k, 2 k, 3
A E,
k. 2
A P,
k, 3
AE,
A W,
x (%, i -
AP,
k, i
AS,
ig
+
&k, iPG,
uk, 2
sk.2
uk, 3
S k, 3
uk, i
A P,k. A’
A E. k,
A W,
A P, k, NPM
k, NPM
&k, iPL,
ig
11 r, Ar
(28)
k. i
A W,k.N
-
This system of difference equations defined by (24) can be transformed into matrix shape (eqn. 29). This matrix is a so-called ‘tridiagonal matrix’. A possible and easy way to handle the solution of this matrix is by using the Gauss algorithm. Referring to the shape of this matrix this algorithm is called the ‘tridiagonal matrix algorithm’. Rearranging, we have eqn. (30).
k, 3
k, i
uL, i)
N
=
sk,
i
uk, N
Sk.
N
uk, NPhj
Sk,
NPM
(29)
106 4 P. k. 2
A E. k, 2 A E,k, 3
k, 3
A,.
A P,
k, i
A,.
The source term S&i is evaluated by sk,
i -
Aw,
k.
fUk,i -
(31)
I
The last line of matrix (30) enables us to calculate the axial velocity at the grid position NPM: uk, NPM
=
S’k. NPM
(32)
k, NPM
1-4~.
By introducing this value into the last-but-one line of the matrix system we obtain the algorithm to evaluate uk,,,,. Repeating this procedure we can determine the whole velocity vector according to *A, i =
SL,
i IAP~
k, i -
AE
k, i”k,
i+
I IAP,
k. i
S’k, 2
uk. 3
s;
3
sL,
i
=
4.1
k, i
A P,k, N
Sk, f =
uk, 2
A E, k, N
uk, N
S’k, N
A P,
uk. NPh
%.
NPN I
(30)
NPll I
Another crucial question related to the numerical procedure is whether this numerical solution of the convergent region is stable against disturbances generated by different initial or boundary conditions. Or, in other words, whether this solution is a unique solution or just one of several possible solutions. To elucidate this problem the initial velocity profiles were varied in several runs from constant to parabolic and finally to turbulent velocity profiles to study their influence on the convergent solutions. The results of this study are depicted in Figs. 8 and 9.
(33)
I
1
5
Boundary conditions Each phase velocity has to fulfil the no-slip condition at the wall of the reactor, so the boundary condition at the wall is: reactor wall (r = I?):
%, k, (r
= R)
-0 -
(34)
On the axis of the reactor (r = 0) symmetry is assumed: reactor axis (r = 0):
(du, k /dr) ((I= 0) = 0
(35)
0.F
‘..d
0
,......
1000
3000
LOO0
Number of Iteration
3. Calculated results
_I
.,I.....
2000
Steps
Fig. 7. Cbnvergence of the numerical procedure: effect of iteration steps on the velocity profile. (Da = 0.15 m; u, = 7.3 cm SC’;
3.1. Checks for plausibility Each model has to be checked to see whether the calculated results are reliable and reasonable within the given limits of the restricting assumptions indispensable in the course of model derivation. First of all, numerical artefacts must be .avoided. The presented solutions of the model equations should be stable and convergent. In this context the question of the minimum number of iteration steps arises. For the critical judgement of this problem we follow the iterative development of a selected value of the calculated solution. Of course, this should be that point of the numerical grid which gets the starting information last. In our case, as we start to integrate from the reactor wall, this point is represented by the axial velocity on the axis of the reactor. Figure 7 shows how the solution runs into the convergence zone after about 2000 iteration steps.
a
3
2
3
.L
.5
Relative
5
.7
.t3
.9
1.
Radial Position (r/R)
Fig. 8. Influence of different initial profiles on u,,~(T). (DR = 0.29 m; u, = 7.2cm SK’; q,L= Z.Ocms-‘; the points for u,, = constant, turbulent and parabolic lie on the same curve.)
107 ‘;i
.3E
2
.3
2 3 .2
.l
0.
-.l 1
2
.3
.L
5
6
Relative
Radial
7
8 Position
9
1
(r/R)
Fig. 9. Influence of different initial profiles on u,, L(r). (DR = 0.29 m; u, = 7.2 cm s-l; u,,, = 2.Ocm s-‘; the points for u,, = constant, turbulent and parabolic lie on the same curve.)
As these two Figures clearly prove, the calculated profiles are not affected at all by different initial conditions. This is very satisfying from a practical point of view, because this means, in reverse, that we do not need precise foreknowledge to start the calculations.
.6
e
.5
cl
.L
a
r* “j 3
.3 .2 .l
0. 1 -2 .l
2
3
.L
.5
Relative Fig.
10.
0,
u=,~(T);
experiment 2.0 cm s-l.)
0
[12, 131.
.l
.2
.3
.b Radial
7
.FJ
Position
9
1
(r/R)
A, u=,~(T): comparison of calculation with (D,=O.lSm; u,=~.~c~Is-‘; u,,,_=
.L Relative
.5
.6 Radial
.7
.0
Position
.9
In the next Figures a comparison of calculated and measured velocity profiles is shown. Figure 10 shows a comparison of our calculations with experimental results found in the literature [ 12, 131. Figure 11 represents some actual investigations of this paper. Both Figures clearly demonstrate the fundamental suitability of the model for multiphase flow systems. Both gas and liquid velocity proties represent fairly well the basic features of the measured velocity profiles. The radial locations of maximum, minimum and turning point of the calculated curves are found at the right positions. The good agreement between calculation and experiment is emphasized by taking into account the fact that no mathematical routine for the optimization of model parameters has been applied yet. 3.2. Parameter
studies
Once a mathematical model is implemented on a computer it is a mighty tool for the simulation and prediction of the behaviour of the reactor. To yield an optimum design of a bubble column a simulation program helps to find out the sensitive parameters which control and dominate the fluid dynamic behaviour. Figures 12 and 13 depict the response of the model on variation of the superficial gas velocity. They prove the applicability of the model within a broad range of operation modes as the trends are represented in an appropriate manner. The slight decrease of the central curvature with high superficial gas velocities is an indication of the influence of the bubble diameter on the phase interaction force via the drag coefficient. For this reason a more detailed investigation of the effect of the bubble diameter was carried out (Figs. 14 and 15). The analysis of the momentum balance reveals the drag force to be dependent on the bubble diameter. Figure 14 shows the effect of different bubble sizes on the drag force. The bigger the bubbles, the smaller the contribution to the force balance. Small bubbles suffer higher loss of momentum and thus decreasing bubble diameters are correlated with
1.
(r/R)
0.
3
.2
.3
.I
.5
Relative Fig. 11. A, ~,,~(r); x, ZJ,,L(r): comparison of calculation and experiment (21. (DR = 0.29 m; conditions for u,, o measurements: z&=4.ocms-‘, uoL = 0.0 cm SC’; conditions for U, L measurements: u, = 4.0 cm s-l, u,, = 1.0 cm s-l.)
6 Radial
3
.8
9
1.
Position (r/R)
Fig. 12. u,,o(r) as a function of supertkial gas &=0.29m; %,=O.Ocms-‘; A, u,=2.0cmP’; &,o=4.0cms-‘; 0, ~~=66.0cms-‘.)
velocity. +,
0.
.1
2
.3
.I
.5
Relative
.6
3
Radial
.6
9
Position
1.
Fig. 13. u,, L(r) as a function of superficial (D,=0.29m; u,,,=O.Ocms-‘; A, ~~=2.Ocms-‘; 4.0 cm s-‘; u, = 6.0 cm SC’.)
a
.l
.2
.3
1
.5
Relative
.6
3
Radial
B Position
0.
1
.2
.3
k/RI
9
.L
.5
Relative
gas
velocity. +, u,=
1.
(r/R)
Fig. 14. Effect of bubble diameter on the drag force (FG/sL). (D,=O.lSm; u,=7.2cms-‘; u0,=2.0cm-‘; A, D,= 0.25cm; +,D,=0.5cm; x,D,=l.Ocm; O,D,=2.Ocm.)
.6 Radial
.7
.8
Position
9
1.
(r/R)
Fig. 16. Effect of different bubble size distributions on velocity profiles. (Da = 0.15 m; u, = 7.2 cm SC’; u,,~ = 2.0 cm S-'; u,, G: D, = 0.55 cm; A, constant; +, parabolic; u,, L: D, = 0.55 cm; x , constant; 0 parabolic.)
Another Figure may emphasize the crucial role of the bubble diameter. Up to this point of the calculations the radial distribution of the bubble diameter was assumed to be constant. Measurements of the bubble size distribution and a theoretical study of the effect of different kinds (e.g. parabolic) of radial bubble size distributions on the fluid dynamics are depicted in Fig. 16. The gas-phase velocity profile is not very sensitive to the type of radial size distribution, but the liquid-phase velocity profile is affected significantly. We observe a more pronounced profile by keeping the bubble diameters constant. In the case of a parabolic distribution we find very small bubbles in the peripheral region of the reactor which hinder the backflowing liquid in a more effective way than big bubbles would do, as can be deduced again from Fig. 14. In order to maintain continuity over a cross-section of the reactor the total sum of the liquid velocity profile is shifted to lower values.
4. Conclusions
-2
2
Relative
Radial
Position
(r/RI
Fig.
15. Effect of bubble diameter on the liquid velocity profile, u,=7.2cms-‘; uoL=2.0cm-‘; A, ur,L. (D,=O.l5m; D,=O.25cm; +,DB=0.5cm; x,D,=l.Ocm; 0,&,=2.O~m.)
decreasing gas-phase velocities. The same proceedings can be observed in the central region of the liquid velocity profile. For bubble diameters smaller than 2.5 mm the curves of Fig. 15 become physically meaningless, because the drag force dominates all other source terms of the momentum balance.
The present work demonstrates a physical approach to multiphase flow systems. The heterogeneous model is able to predict the occurrence of recirculation and therefore is a more precise tool than flowcharts. The agreement of experiment and theory is fairly good. Parameter studies of the contributions to the momentum balance make clear the role of the decisive parameters, such as bubble diameters. With regard to the reciprocal effect of fluid dynamics on absorption/chemical reactions this approach seems to be very promising. At present it is being extended to ‘hot systems’, including homogeneous as well as heterogeneous chemical reactions.
Acknowledgement The authors gratefully acknowledge the financial support of the Deutsche Forschungsgemeinschaft for this work.
109
Nomenclature
4 -4 Aw CD
z$n DR g
F 4
m m* m’ P x S t uOk %,
k
%,
k
V,X xi Z
6 &k
coefficient for grid point P, kg m-’ s-’ coefficient for grid point E, kg m-’ s-’ coefficient for grid point W, kg m-’ s-’ drag coefficient eddy viscosity factor bubble diameter, m reactor diameter, m acceleration due to gravity, m se2 reactor height, m turbulent kinetic energy, mz se2 characteristic length scale, m superficial mass flux, kg m--2 s-l estimated superficial mass flux, kg m -2 s ’ corrected superficial mass flux, kg mm2 s-’ pressure, N m-’ radial coordinate, m reactor radius, m source term, kg m-’ SC’ time, s superficial phase velocity, m s- 1 mean axial phase velocity, m S-I mean radial phase velocity, m SC’ mean bubble rise velocity, m s-l general space coordinate, m axial coordinate, m Kronecker’s delta phase fraction of phase k dissipation rate, m2 ss2 viscosity, kg m- 1 s- ’ kinematic viscosity, m2 s- ’ density, kg rn-’ shear stress, N me2 tangential direction (cylindrical rad
Subscripts B E eff G i
bubble east side of finite volume effective gas radial grid position
coordinates),
i k L lam N N NPM P r S I, Z v
axial grid position phase index (gas or liquid) liquid laminar north side of finite volume number of grid points maximum number of grid points grid point (position i, j) radial south side of finite volume turbulent west side of finite volume axial tangential direction
( = N + 1)
References 1 W.-D. Deckwer, Reaktionstechnik in Blasenstiulen, Salle + Sauerliinder, Aarau, 1985. 2 H.-E. Gasche, DissPrtation,UniversitPt Erlangen-Ntirnberg, 1989. 3 H. Hofmann and H.-E. Gasche, Investigations of the fluid dynamics of bubble columns, Proc. German/Japanese Symp. on Bubble Columns, 1988, Schwerte, F.R.G. 4 T. Miyauchi and C. N. Shyu, Kagaku Kogaku, 34 (1970) 985. 5 K. Ueyama and T. Miyauchi, AZChE J., 25 (1979) 258. 6 Th. Menzel and U. Gnken, Deterministische und stochastische Komponenten der Fliissigkeitsstr6mung in Gas/Fltissigkeitsreaktoren, Proc. 2% Tutzing Symp: Znstationdre Proresse und Mejlmethoden in Katalyse, Etektrochemie und Korrosion, 1989. 7 G. Scheuerer, Vorlesungsskript: Numerische LGsung der strdmungsmechanischen Grundgleichungen, Lehrstuhl fur Striimungsmechanik, Univ. Erlangen-Nihnberg. 1986. 8 M. Sommerfeld, Vorlesungsskript: Mehrphasenstrtimungen, Lehrstuhl fiir Stromungsmechanik, Univ. Erlangen-Numb&g, 1987. 9 R. B. Bird, W. E. Stewart and E. N. Lightfoot, Transport Phenomena, Wiley, New York, 1960. 10 H. Kompel, Diplomarbeit, Inst. Tech. Chem. I, Univ. Erlangen-Niirnberg, 1989. 11 C. Edinger, Diplomarbeit, Inst. Tech. Chem. I, Univ. ErlangenNiimberg, 1988. 12 J. Steinemann, Dissertation, Univ. Dortmund, 1985. 13 K. Franz, Dissertation, Univ. Dortmund, 1983.