Numerical studies of bubble and particle dynamics in a three-phase fluidized bed at elevated pressures

Powder Technology 112 Ž2000. 46–56 www.elsevier.comrlocaterpowtec

Numerical studies of bubble and particle dynamics in a three-phase fluidized bed at elevated pressures Jianping Zhang, Yong Li, Liang-Shih Fan ) Department of Chemical Engineering, The Ohio State UniÕersity, 140 West 19th AÕenue, Columbus, OH 43210, USA

Abstract A discrete phase simulation is conducted to study the bubble and particle dynamics in a three-phase fluidized bed at high pressures. The Eulerian volume-averaged method, the Lagrangian dispersed particle method, and the volume of fluid ŽVOF. method are employed to describe, respectively, the motion of liquid, solid particles, and gas bubbles. A bubble-induced force model, a continuum surface force ŽCSF. model, and Newton’s third law are applied to illustrate, respectively, the coupling effect of particle–bubble, gas–liquid, and particle–liquid interactions. A close-distance interaction ŽCDI. model is included in the particle–particle collision analysis, which considers the liquid interstitial effect on colliding particles. Effects of the pressure and solids holdup on the bubble rise characteristics such as the bubble rise velocity, bubble shape and trajectory are examined. Simulations of the bubble rise velocity at various solids holdups and pressures are conducted along with the maximum stable bubble size and the particle–bubble interactions. The simulated results compare favorably with the experimental data and calculation from a mechanistic model. q 2000 Elsevier Science S.A. All rights reserved. Keywords: Discrete phase simulation; Gas–liquid–solid fluidization; Bubble; Particle; Volume of fluid method; Dispersed particle method; High pressure

1. Introduction Gas–liquid–solid fluidization systems have been widely utilized as reactors in chemical and petrochemical industries. Many bubble column reactors are operated at highpressure conditions, e.g., methanol synthesis Žat P s 5.5 MPa and T s 08C., resid hydrotreating Žat P s 5.5–21 MPa and T s 300–4258C., Fischer–Tropsch synthesis Žat P s 1.5–5.0 MPa and T s 2508C. and benzene hydrogenation Žat P s 5.0 MPa and T s 1808C. w10,14,23,25,27x. The design and scale-up of these reactors require knowledge of the hydrodynamics as well as heat and mass transfer of the system. Studies in the literature have indicated significant effects of pressure on the hydrodynamics and transport phenomena in bubble columns w13,32x, slurry bubble columns w29x and three-phase fluidized beds w16,19x. It is found that elevated pressures lead to higher gas holdups in these systems. Experimental studies have found that the increased gas holdup can be related to the smaller bubble size w16,17x, slower bubble rise velocity w20x, reduced bubble coalescence rate and increased bubble ) Corresponding author. Tel.: q1-614-292-4935; fax: q1-614-2923769. E-mail address: [email protected] ŽL.-S. Fan..

breakup rate w21x at higher pressures. The variation in physical properties of the gas and liquid with pressure is considered to be the most fundamental reason for bubble size reduction. At a higher pressure, a larger gas density leads to a smaller initial bubble size from the gas distributor, a lower surface tension and a larger liquid viscosity contribute to the reduced bubble coalescence rate, and a larger gas density and a smaller surface tension contribute to the increased bubble breakup rate. The characteristics of a rising bubble can be described in terms of the rise velocity, shape, and motion of the bubble. These rise characteristics are closely associated with the flow and physical properties of the surrounding medium as well as the interfacial properties. The bubble rise velocity is the single most critical parameter in characterizing the hydrodynamics and transport phenomena in liquids and liquid–solid suspensions w8x. Under certain conditions, the rise velocities of single bubbles in liquid– solid fluidized beds were found to be similar to those in highly viscous liquids. The liquid–solid suspension was, thus, represented by an effective viscosity, which accounts for the effect of particles on the rise of single bubbles w5,6,11,30x. However, when considering the rise velocity of smaller bubbles Ž d B - 12–17 mm., the effective viscosity of the liquid–solid medium deviates from the viscosity

0032-5910r00r$ - see front matter q 2000 Elsevier Science S.A. All rights reserved. PII: S 0 0 3 2 - 5 9 1 0 Ž 9 9 . 0 0 3 0 4 - 6

J. Zhang et al.r Powder Technology 112 (2000) 46–56

of the corresponding Newtonian liquid w4,8,30x. The deviations that mark the reduction in the bubble rise velocity reflect significant close-range interactions of the bubble with the liquid–solid media or with individual particles. A more elaborate analysis is required to account for the heterogeneous characteristics of liquid–solid suspensions with regard to the rising bubble. Jean and Fan w15x developed a mechanistic model based on a force balance on a rising bubble involving the net gravity, liquid drag, and particley bubble collision forces. The model can predict the bubble rise velocity in liquid– solid media for small particles Ž d p - 500 mm., low-to-intermediate solids holdups Ž ´s - 0.45., and large spherical bubbles Ž d B ) 15 mm.. Luo et al. w20x, in studying the bubble rise velocity in liquid–solid suspensions at elevated pressure and temperature, extended Jean and Fan’s model to cover a smaller bubble size range. However, their results have shown that the model overestimates the bubble rise velocity at low temperatures, while better agreement with the experimental data is obtained at higher temperatures. Due to the complex nature of the hydrodynamics of bubble flow in liquid–solid media, it is difficult to obtain a mechanistic model that can be used to calculate the bubble rise velocity in various physical properties and system parameters. For wide-range predictions of the bubble rising characteristics in liquid–solid suspensions and of the dynamic behavior of gas–liquid–solid fluidization systems including the interactions of individual bubbles and particles, numerical simulations based on the computational fluid dynamics are required. In this study, the bubble and particle dynamics in a gas–liquid–solid fluidized bed at elevated pressures are numerically studied based on a discrete phase simulation. The Eulerian volume-averaged method, the Lagrangian dispersed particle method and the volume of fluid ŽVOF. method are employed to describe, respectively, the motion of liquid, solid particles, and gas bubbles. A bubble-induced force model, a continuum surface force ŽCSF. model, and Newton’s third law are applied to illustrate, respectively, the coupling effect of particle–bubble, gas–liquid, and particle–liquid interactions. A close-distance interaction ŽCDI. model is included in the particle–particle collision analysis, which considers the liquid interstitial effects between colliding particles. The flow characteristics studied include single-bubble rise velocity at various solids holdups, the stable maximum bubble size, the bubble breakup, and the bubble–particle interactions. Comparisons of computational results with experimental data are also made. 2. Computational models 2.1. Volume-aÕeraged equations and VOF method In this study, the volume-averaged Navier–Stokes equations are used to solve the liquid phase flow Žoutside the

47

gas bubble. and the gas phase flow Žinside the gas bubble. in the presence or absence of dispersed particles. The liquid phase turbulence is not included in this study. However, as a simulation program for dynamic flow, the current method is able to simulate the temporal fluctuations in three-phase fluidization systems since very small time step, ; 10y6 s, is used in simulations. The VOF method is used to simulate the motion and the topological change of gas bubbles w12x. The volume-averaged mass and momentum equations are given as the continuity equation: E´ Et

q = P Ž ´n . s 0

Ž 1.

and the momentum equation:

r

E Ž ´n . Et

q r= P Ž ´nn . s y= p q ´= P t q ´r g q f b Ž 2 .

where n is the velocity vector; ´ is the holdup of liquid or gas phase for the liquid–solid mixture ´s q ´ l s 1, for the gas phase with solid particles ´s q ´g s 1; r is the fluid density; p is the fluid pressure; t is the viscous stress tensor; g is the acceleration of gravity; and f b is the total volumetric body force, which includes the forces acting on the liquid–gas from the particles and bubbles Žgas–liquid interface., f b s f pf q f bf . Based on Newton’s third law of motion, the forces acting on a particle from the liquid phase, Ffp , yield a reaction force on the liquid. Therefore, the momentum transfer from particles to the liquid–gas phase is taken into account in Eq. 2 by adding the volumetric particle–fluid interaction force, f pf , given below to the body force term: f pf s y

Ý Ffpk DVi kj

,

x pk g V i , j

Ž 3.

where subscripts i, j define the location of a computational cell; V and DV are the domain and volume of this cell, respectively; x pk is the location vector of particle k; Ffp is the fluid–particle interaction force acting on any individual particle, which includes the drag, added mass, and Basset force: Ffp s F D q FAM q FBA

Ž 4.

It is noted that the pressure gradient is not included in Eq. 4 as, for consistency w33x, this force has been accounted for in the liquid phase momentum equation given in Eq. 2. The liquid properties on a particle are obtained by the area-weighted averaging over the properties at the four grid points of the computational cell containing the particle. For solving the volume-averaged equations for the liquid phase, the liquid holdup, ´ l , is obtained as the

J. Zhang et al.r Powder Technology 112 (2000) 46–56

48

volume fraction of the liquid in the computational cell. For evaluating particle drag coefficient in the liquid–solid medium, the liquid holdup based on the particle-centered averaging over a prescribed area is, however, used in the computational scheme. The bubble–fluid interaction force, f bf , is obtained by using a CSF model w2x: f bf s "d Ž x ,t . sk Ž x ,t . =a Ž x ,t . ,

Ž 5.

where the plus sign is used for the liquid phase and the minus sign is used for the gas phase; d Ž x,t . is the Delta function, which equals 1 at the gas–liquid interface and equals 0 elsewhere; s is the surface tension of the interface; k Ž x,t . is the curvature of the free surface; a Ž x,t . is the volume fraction of the fluid a Ž x,t . s 1 in the liquid or liquid–solid mixture, 0 - a Ž x,t . - 1 at the free surface, and a Ž x,t . s 0 in the gas phase. The advection equation for a Ž x,t . is given as: Ea Et

q Ž n P = . a s 0.

Ž 6.

At the gas–liquid interface, the property discontinuity is approximated by a smooth variation of mixture properties from one phase to the other within a finite interface thickness. The mixture properties are then calculated using a volume fraction weighted averaged method as:

r s r 1 P a Ž x ,t . q rg P Ž 1 y a Ž x ,t . . m s m 1 P a Ž x ,t . q mg P Ž 1 y a Ž x ,t . .

The drag force acting on a suspended particle is proportional to the relative velocity between the phases as follows: FD s

1 2

CD s

24 Re p

mp

d np dt

Ž 1 q 0.15 Rep0.687 . ,

C D s 0.44,

Re p G 1000

Ž 12 .

The added mass force accounts for the resistance of the fluid mass that is moving at the same acceleration as the particle. For a spherical particle, the volume of the added mass is equal to one-half of the particle volume, so that: 1 2

r Vp

d dt

Ž n y np .

Ftotal s FD q FAM q FG q FB q FBA q Ý Fi i

Ž 9.

Ž 13 .

The Basset force accounts for the effect of past acceleration. The original formulation of the Basset force is derived based on the creeping flow condition. For a particle moving in a liquid with a finite Reynolds number, the modified Basset force is given as w22x: FBA s 3pm d p

t

H0 K Ž t y t .

d Ž n y np . dt

dt ,

Ž 14 .

K Ž t y t . in Eq. 14 is given as:

° p Ž tyt .n ¢ r

K Ž t y t . s~

1r4

2 p

q

Ž 8.

The total force acting on a particle is composed of all applicable forces, including drag Ž FD ., added mass Ž FAM ., gravity Ž FG ., buoyancy Ž FB ., Basset history force Ž FBA ., and other forces ŽÝ Fi .:

Ž 11 .

CXD s C D ´y4 .7

1 s Ftotal .

Re p - 1000

where Re p s r l < n < p d prm l . In the liquid–solid or gas–solid suspensions, the drag force depends strongly on the local phase holdup in the vicinity of the particle under consideration. The effective drag coefficient can be obtained by the product of the drag coefficient for an isolated particle and a correction factor as given by w31x:

FAM s

The motion of a particle in a flow field can be described in the Lagrangian coordinate system with its origin set at the center of the moving particle. The particle movement in a non-uniform flow field includes acceleration and rotation. Since the particle size is smaller than the grid size in the present method, the particle rotational motion is neglected in the dynamic particle motion equation. The change of particle rotation due to collision is, however, considered, as will be discussed in Section 2.3. The translational motion of a particle in the liquid or gas phase is governed by Newton’s second law of motion:

Ž 10 .

where A is the exposed frontal area of the particle to the direction of the incoming flow, C D is the drag coefficient, which is a function of the particle Reynolds number, Re p . For rigid spherical particles, the drag coefficient C D can be estimated by the following equations w26x:

Ž 7.

2.2. Dispersed particle method

C D r A < n y n p < Ž n y np .

2

p

Ž U q np y n . rp n f H3 Ž Re .

¶• ß

1r2 y2

3

Ž tyt .

2

,

Ž 15 . f H Ž Re . s 0.75 q 0.105Re Re s Ud prn where n is the kinematic viscosity of the fluid; U is the mean stream velocity.

J. Zhang et al.r Powder Technology 112 (2000) 46–56

The sum of the gravity and buoyancy forces has the form: FG q F B s Ž rp y r . Vp g

Ž 16 .

If particles move into the interface area, which is identified as 0 - a Ž x,t . - 1, a bubble-induced force is added to the particle: F bp s Vp d Ž x ,t . sk Ž x ,t . =a Ž x ,t .

Ž 17 .

A hard sphere approach is used for the particle–particle collision analysis. In this approach, it is assumed that collisions between spherical particles are binary and quasi-instantaneous, and further, that there is a sequence of collisions during each time step. The equations, which are similar to the equation of molecular dynamic simulation w1x, are used to locate the minimum flight time of particles before any collision. 2.3.1. Liquid shear effect While the shear force can be neglected in gas–solid flow systems, this is not the case in liquid–solid systems. The liquid shear effect between particles becomes important when two particles move close to each other in liquid–solid systems, especially when the distance between two particles is less than 0.1d p w34,35x. Thus, the CDI model is used to locate the particle contact velocity just before collision, which considers the strong damping effect due to the liquid film before particle contact. The particle normal contact velocity can be described by w34,35x:

ž

1q

s

1 rl

3 r l rp3

q

2 rp

16 rp h

3

9 m f f Ž np y n l . 2 rp2 rp

ž

y 1y

/

d np

rl

32 rp h 4

rp

/

9m

g

np

q

up

f s exp

rp

t

H0 K Ž t y t .

d Ž np y n l . dt

dt

2 rp2 rp np

žž / ž / 1.7

f s 1 q 0.15Re p0.687

r

0 .19

¢m U

se

N N N a a q m b Ub s m a Ua

X

Ž 20 .

q m b UbN

X

where U N is the normal velocity of the particle Ža or b. at the contact point before or after Žwith superscript X . the collision. From the contact theory of Mindlin w24x, there are three kinds of frictional contact during the collision: Ž1. sliding contact; Ž2. non-sliding or sticking contact; and Ž3. torsion of elastic particles in contact. By neglecting the effect of particle torsion during collision, the simplified Mindlin’s contact theory is applied to obtain the tangential components after the collision. If the incident angle, defined as the ratio of the particle–particle relative velocity in the tangential direction to velocity in the normal direction, is less than the critical angle Ž a cr s tany1 Ž2 f k ., where f k is the friction coefficient., the sticking collision occurs: X

UaT s UbT

X

Ž 21 .

Otherwise, the sliding collision occurs, in which w9x: X

rp

ž / h

R e p0 .47

/

,

Ž 19 .

Using the Runge–Kutta method, Eq. 18 can be solved to determine the particle normal contact velocity just before the collision.

Ž 22 .

where U T is the tangential particle velocity Ža or b. at the contact point. The conservation of momentum is given as: X

where h is the distance from the center of the approaching particle to the midpoint between the two particles; rp is the radius of particle; and f w28x and F w34,35x are the correction functions and can be expressed as: 0.44

UbN y UaN

m aUaT q m b UbT s m aUaT q m b UbT

Ž 18 .

Re p

~

X

X

9 r l rp4 Ž np y n l . < np y n l <

q

X

UaN y UbN

Ž UaT y UbT . y Ž UaT y UbT . s 2 f k Ž UaN y UbN .

dh

up

2.3.2. Particle collision analysis When two particles are in contact, collision analysis can be conducted to obtain the velocities of the particles after collision. It is assumed that tangential traction and the resulting displacements have no effect on normal collision. For the collision between particles a and b, the normal components after collision can be obtained by solving the equations for the restitution coefficient and the conservation of momentum:

°

2.3. Particle–particle collision dynamics

49

X

Ž 23 .

The tangential velocities after the collision can be obtained by solving Eq. 21 or Eq. 22 together with Eq. 23. As mentioned in Section 2.2, the collision induces a change in particle rotation. The angular velocities after the collision are determined by: X

½

Ia Ž vaX y va . s m a Ž UaT y UaT . P ra X

I b Ž v Xb y v b . s m b Ž UbT y UbT . P r b

Ž 24 .

where v is the angular velocity of the particle Ža or b., and I is the moment of inertia defined by I s 2r5m p rp2 . The tangential velocities of the particle center are given as:

½

X

X

X

X

UacT s UaT y vaX P ra UbcT s UbT y v Xb P r b

Ž 25 .

J. Zhang et al.r Powder Technology 112 (2000) 46–56

50

The velocity of the particle before each collision is updated within the time duration using Eq. 27. The detailed procedure of calculating the particle movement is shown in Fig. 1.

3. Results and discussion 3.1. Simulation conditions

Fig. 1. The flowchart of the particle phase simulation.

2.4. Particle moÕement Within a time step of advance, D t, a particle moves to a new position according to: x p s x p ,0 q np P D t Ž 26 . if no collision is encountered. The particle velocity is updated using a simple explicit integration formula: d np np s np ,0 q Dt Ž 27 . dt where the acceleration of the particle is obtained from Eq. 8. The equations, which are similar to those in the molecular dynamic simulation w1x, are used to determine the minimum flight time of particles before any collision. If the collisions occur within time step D t, this time step is split into several flight time steps, D t c i , and: D t s Ý D tc i

The computational models presented in Section 2 are incorporated into a two-dimensional code of the discrete phase simulation of gas–liquid–solid fluidization, CVD-2, developed earlier by the present authors w18,36,37x. Details of the simulation scheme and numerical methods are given in Refs. w36,37x. The flowchart of the main program of the discrete phase simulation of gas–liquid–solid fluidization systems is shown in Fig. 2. In this study, the numerical simulations are conducted for the bubble and particle dynamics in three-phase fluidization systems at elevated pressures. Paratherm NF heat transfer fluid, nitrogen, and glass beads are used as liquid, gas, and solid phases, respectively. The hydrodynamics of a bubble rising in a liquid–solid fluidized bed is numerically studied under an elevated pressure of 17.3 MPa and a temperature of 758C. To study the pressure effects, the simulations for a bubble rising in similar conditions but atmospheric pressure are also conducted. The physical properties of the gas and liquid phases under the simulation conditions are given in Table 1. The solid particles are spherical with a density of 2500 kgrm3 and a diameter of 0.088 cm. The solids holdups studied are 0.384 and 0.545. To simulate a bubble rising in a liquid–solid fluidized bed, initially, a spherical bubble is positioned in the computational domain with its center located 1.5 cm above the

Ž 28 .

i

The flight time of a particle is determined as the minimum flight time of this particle to its neighboring particles andror walls. The flight time for two particles, a and b, is obtained by:

(

2

2 yrab P nab y Ž rab P nab . y nab2 rab y Ž Ra q R b .

tab s

2

nab2

Ž 29 . where R a and R b are radii of particles a and b, respectively, and rab s ra y r b and nab s na y n b . For the particle-wall collision, the flight time is given as: Ž < x wall < q R a . y < r x ,a < ta ,wall s Ž 30 . zx ,a

Fig. 2. The flowchart of the main program for discrete phase simulation of gas–liquid–solid fluidization systems.

J. Zhang et al.r Powder Technology 112 (2000) 46–56 Table 1 Physical properties of the gas and liquid used in simulations Pressure ŽMPa. Gas density Žkgrm3 . Liquid viscosity ŽPa s. Liquid density Žkgrm3 . Surface tension coefficient ŽNrm.

0.1 0.94 0.00379 843 0.0253

Table 2 Comparison of the simulation results with the experimental data for the rise velocity of a bubble 7.5 mm in diameter

17.3 154 0.00417 872 0.0192

Pressure ŽMPa. Solids holdup Ž – . Simulated velocity Žcmrs. Velocity from experiments Žcmrs. Relative error Ž%.

bottom; the particles are randomly positioned in a 3 = 24 cm2 area. Particles then settle under gravity in the liquid medium with an evenly distributed inlet velocity. The spherical bubble is treated as a stationary obstacle at this stage of the simulation. At a certain liquid inlet velocity, the desired solids holdup is reached in a 3 = 8 cm2 area. The simulation is then restarted in this 3 = 8 cm2 domain with the liquid flow field, bubble tracking, and particle movement calculated at each time step. The time step of simulation for the liquid and solid phases is 5 = 10y6 s. The inlet liquid velocity is 1.0 and 0.45 cmrs for the solids holdups of 0.384 and 0.545, respectively. The simulations are conducted in a Cray T90 supercomputer at the Ohio Supercomputer Center. The memory required to run the code is six MegaWords and about 20 h of CPU time are required for a typical case presented in this paper. 3.2. Bubble rise Õelocity In this study, the simulation is performed in a two-dimensional column. In order to compare with the experimental data for the bubble rising in liquid–solid media at high pressures obtained in a three-dimensional column, a conversion for the bubble rise velocity of two-dimensional case to that of three-dimensional case is required. It is noted that for large bubbles, the theoretical Davis y Taylor w7x equation gives the two-dimensional bubble rise velocity UB 2 s 1r2Ž gR .1r2 and the three-dimensional bubble rise velocity UB 3 s 2r3Ž gR .1r2 . Therefore, the conversion factor of the bubble rise velocity from two-dimensional to three-dimensional is UB 3 s 4r3UB 2 . In the current study, this conversion factor is employed as an approximation to calculate based on the results of the simulated two-dimensional bubble rise velocity, the three-dimensional rise velocities of bubbles of small diameter Ž d B - 2 cm.. Additionally, the following equation is used to correct the wall effect w3x on the bubble rise velocity: UB UB`

dB D

0.1 0.384 20.6 20.8 1.08

17.3 0.384 15.4 15.7 1.85

17.3 0.545 11.9 12.3 3.2

in the liquid–solid medium under various solids holdups and system pressures. The reduction of the bubble rise velocity with an increase in pressure is a reason for the significant increase in the gas holdup of a three-phase fluidized bed at elevated pressures. When the pressure is increased from 0.1 to 15.6 MPa, a 100% increase of gas holdups was reported at all gas velocities w19x. 3.3. Bubble shape and trajectory Similar to the rise velocity of a bubble, the shape of a bubble is affected by the physical properties and system parameters, including the surface tension and viscosity of the liquid, densities of liquid and particles, solids holdup, bubble size, and system pressure and temperature. The aspect ratio of a bubble, hrb, defined as the ratio of the minor axis over the major axis of the bubble, is the parameter used to characterize the bubble shape. The simulation results of the bubble aspect ratio changing with time during the bubble rising in the liquid–solid fluidized bed, with a solids holdup of 0.384 at pressures of 0.1 and 17.3 MPa, and a solids holdup of 0.545 at pressure of 17.3 MPa, are shown in Fig. 3. As shown in the figure, the bubble aspect ratio is similar at different pressures at the same solids holdup, i.e., 0.384, although the bubble rise velocity is reduced from 20.6 to 14.0 cmrs with the increasing pressure as shown in Table 2. In the high solids holdup, 0.545, the bubble aspect ratio is larger than that in the low solids holdup. When rising in pure liquids, the bubble aspect ratio increases with the liquid viscosity. The presence of solid particles yields a similar effect as in-

2 3r2

ž ž //

s 1y

51

.

Ž 31 .

Eq. 31 is applicable for bubbles and drops for Eo - 40, Re ) 200, and d B rD F 0.6. The comparison of the simulation and the experimental results for the bubble rise velocity is given in Table 2. As shown in the table, the simulation results agree well with the experimental data for the rise velocity of a gas bubble

Fig. 3. The simulated results of the bubble aspect ratio changing with time at various solids holdups and pressures.

52

J. Zhang et al.r Powder Technology 112 (2000) 46–56

creasing liquid viscosity with regard to the bubble rise characteristics. The simulation results indicate that the pressure effect is insignificant to the bubble aspect ratio. On the other hand, the increase of the solids holdup has an appreciable effect on the bubble aspect ratio. From experimental studies of bubble rise characteristics in liquid–solid fluidized systems under various pressures and temperatures w8,20,30x, it is noted that the particle effect is small at low solids holdup Ž ´s - 0.4. and is significant at high solids holdup. The change of physical properties and system parameters can also affect the bubble rising trajectory in addition to affecting the rise velocity and shape of the bubble. The bubble rising trajectory, in turn, can also influence the bubble rise velocity. The bubble motion or its rise path is closely related to the wake shedding phenomena w8x. In the bubble Reynolds number range 20 - Re B - 200, the bubble wake is closed; it becomes open with asymmetric shedding at 200 - Re B - 5000, which leads to a zigzag path or bubble rocking. The simulated trajectories of a bubble rising in the liquid–solid fluidized bed with a solids holdup of 0.384 at pressures of 0.1 and 17.3 MPa, and a solids holdup of 0.545 and a pressure of 17.3 MPa are shown in Fig. 4. The time step between two bubbles in Fig. 4 is 0.05 s. The corresponding Reynolds numbers for bubbles in Fig. 4a, b and c are 343.7, 241.5 and 198.5, respectively. The bubbles in Fig. 4a and b fall in the

asymmetric wake shedding regime while the bubble in Fig. 4c is within the transition regime. The tortuous trajectories are evident in Fig. 4a and b. For the same solids holdup, the bubble trajectory is more tortuous at high pressure than at low pressure. Note that the liquid viscosity increases and the surface tension decreases with the increasing pressure. The difference of the bubble trajectories at various pressures in the same solids holdup may result from the competing effects of the liquid viscosity and surface tension, which affect the bubble shape and, hence, wake shedding phenomena. For bubbles rising at the same pressure but with different solids holdups, it is seen that the trajectory of the rising bubble is more stable at high solids holdup than that at low solids holdup as shown in Fig. 4b,c. This is consistent with the viscosity effect on the bubble rising trajectory since increasing the solids holdup would give a similar effect as increasing liquid viscosity. However, at a high solids holdup, more interactions between particles and the bubbles are encountered. Therefore, more bubbles with irregular shapes are observed at a high solids holdup as shown in Fig. 4c. 3.4. Bubble–particle interactions As the bubble rises in the liquid–solid fluidized bed, the interaction between the bubble and particles takes place. In this study, the bubble–particle interaction is accounted for

Fig. 4. Bubble rising trajectory at different pressures and solids holdups. Ža. p s 0.1 MPa, ´s s 0.384, d B s 7.5 mm; Žb. p s 17.3 MPa, ´s s 0.384, d B s 7.5 mm; Žc. p s 17.3 MPa, ´s s 0.545, d B s 7.5 mm.

J. Zhang et al.r Powder Technology 112 (2000) 46–56

53

Fig. 5. The simulated results of bubble–particle interactions Ž p s 17.3 MPa, ´s s 0.384, d B s 7.5 mm.. Ža. t s t 0 ; Žb. t s t 0 q 0.15 s; Žc. t s t 0 q 0.30 s; Žd. t s t 0 q 0.45 s.

by adding a surface-tension-induced force to the particle motion equation. This force is also added to the source term of the liquid momentum equation for the liquid elements in the interfacial area to account for the particle effect on the interface. The particle movement is determined based on the resulting total force acting on the particle. From the simulation results, it is seen that most of particles contacting the bubble do not penetrate the bubble. Instead, they pass around the bubble surface. When the particles penetrate the bubble, they fall quickly to the bubble base because of the low viscosity and density of the gas phase. Fig. 5 illustrates the particle–bubble interac-

tions in four snapshots of the simulated results. As shown in the figure, most of the particles on the bubble surface do not penetrate the bubble; it can be seen that only one or two particles penetrate. 3.5. Bubble breakage and stability In addition to the bubble rise velocity, the bubble size is another important parameter determining the overall gas holdup in a three-phase fluidized system. Luo et al. w21x measured bubble sizes under various pressures, and found that the maximum stable bubble size decreases with in-

Fig. 6. The simulated sequence of bubble shape change and the bubble breakage in a liquid–solid medium Ž p s 17.3 MPa, ´s s 0.384, d B s 10.0 mm.. Ža. t s t 0 ; Žb. t s t 0 q 0.05 s; Žc. t s t 0 q 0.10 s; Žd. t s t 0 q 0.15 s; Že. t s t 0 q 0.20 s.

J. Zhang et al.r Powder Technology 112 (2000) 46–56

54

Fig. 7. The simulated sequence of bubble shape change and the bubble breakage in a liquid medium Ž p s 17.3 MPa, ´s s 0, d B s 10.0 mm.. Ža. t s t 0 ; Žb. t s t 0 q 0.05 s; Žc. t s t 0 q 0.10 s; Žd. t s t 0 q 0.15 s; Že. t s t 0 q 0.20 s.

creasing pressure. They found that the centrifugal force induced by internal circulation of gas inside a bubble can disintegrate the bubble at high pressures. A mechanistic model was developed by them to account for the maximum stable bubble size, d B ma x, as given below:

(

d B ma x f C

Fig. 8 shows the velocity field of the gas and liquid phases before the bubble breakage in a moving coordinate with the rise of the bubble. The bubble rises in a liquid–

s g rg

i

Ž C s 2.53 for a s 0.21, C s 3.27 for a s 0.3 . Ž 32 . where a is the aspect ratio of bubble. This model reveals that the gas inertia and gas–liquid surface tension dictate the maximum stable bubble size at high pressures. From Eq. 32, the maximum stable diameter of a bubble is predicted as 9.0 mm at 17.3 MPa. The rise trajectory of a bubble with the equivalent diameter of 7.5 mm is shown in Fig. 4, which indicates the bubble rises without significant breakage. The breakup of the bubble is simulated when the bubble diameter is larger than the maximum stable bubble for a given flow condition. Fig. 6 shows a series of bubble shape changes for a bubble of 10.0 mm in diameter rising in a liquid–solid fluidized bed with a solids holdup of 0.384 at 17.3 MPa. As shown in the figure, the rising bubble changes its shape drastically and eventually breaks into three parts. The largest part of the breaking bubble has an equivalent diameter of about 9.0 mm and rises without further breakage. Similar results are obtained for the bubble rising in pure liquid as shown in Fig. 7, which shows the relatively small effect of solids concentration on the maximum bubble size and on the gas holdup as observed by Luo et al. w21x. The results of the simulated bubble breakage agree well with the experimental data and the predictions from the mechanistic model.

Fig. 8. Velocity vector field of gas and liquid phases before the bubble breakup Ž ps17.3 MPa, ´s s 0.384, d B s10.0 mm..

J. Zhang et al.r Powder Technology 112 (2000) 46–56

solid fluidized bed with a solids holdup of 0.384 and a pressure of 17.3 MPa. As shown in the figure, a significant internal circulation is observed in the bubble. Luo et al. w21x assume that the centrifugal force induced by this internal circulation is the main driving force for the bubble breakup at elevated pressures. The simulation results of this study agree with the calculations from the mechanistic model developed based on this assumption. 4. Concluding remarks The simulations indicate that the bubble aspect ratio increases with increasing solids holdup. For the same solids holdup, the bubble trajectory is more tortuous at high pressure than at low pressure. The trajectory of a rising bubble is more stable at high solids holdup than that at low solids holdup. For bubble–particle interactions, it is found that most of the particles in front of a bubble pass around the bubble surface while only very few particles penetrate the rising bubble. The results of the simulated bubble rise velocity, bubble breakage, and the maximum stable bubble size at high pressures agree well with the experimental data and calculations from a mechanistic model. 5. Nomenclature A CD d D DB ma x e Eo F f f fk g h I p r ra rb r x,a R Ra Rb Re t U UB UB 2 UB 3

Area Drag coefficient Diameter Column diameter Maximum stable bubble diameter Restitution coefficient Eotvos ¨ ¨ number Force Volumetric body force Correction function Friction coefficient Gravity Separation distance, minor axis of bubble Moment of inertia Scalar pressure Radius of particle Position vector for particle a Position vector for particle b Position vector from the wall to particle a Volume equivalent bubble radius Radius of particle a Radius of particle b Reynolds number Time Mean stream velocity, particle velocity Bubble rise velocity Two-dimensional bubble rise velocity Three-dimensional bubble rise velocity

UB` u V n x x wall

55

Bubble rise velocity in an infinite container velocity Volume Velocity vector Coordinate vector Position vector for the wall

Greek letters a Volume fraction of fluid a cr Critical angle ´ Holdup F Correction function k Free surface curvature m Dynamic viscosity n Kinematic viscosity r Density s Surface tension t Viscous stress tensor, time v Angular velocity Subscripts AM Added mass a Particle index ac Center of particle a b Particle index bc Center of particle b B Basset, bubble bf Bubble–fluid interaction bp Bubble–particle interaction D Drag fp Fluid–particle interaction GrB Gravityrbuoyancy g Gas phase i, j Cell indices l Liquid phase M Magnus P Pressure p Particle pf Particle–fluid interaction Superscripts k Particle index N Normal direction T Tangential direction Acknowledgements This work is sponsored, in part, by the NSF Grant CTS-9528380 and the Ohio Supercomputer Center. References w1x M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. w2x J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comp. Phys. 100 Ž1992. 335.

56

J. Zhang et al.r Powder Technology 112 (2000) 46–56

w3x R. Clift, J.R. Grace, M.E. Weber, Bubbles, Drops, and Particles, Academic Press, New York, 1978. w4x R.C. Darton, The physical behavior of three-phase fluidized beds, in: J.F. Favidson, R. Clift, D. Harrison ŽEds.., Fluidization, Academic Press, London, 1985, pp. 495–528, Chap. 15. w5x R.C. Darton, D. Harrison, The rise of single gas bubbles in liquid fluidized beds, Trans. Inst. Chem. Eng. 52 Ž1974. 301. w6x S. Dasgupta, R. Jackson, S. Sundaresan, Turbulent gas particle flow in vertical risers, AIChE J. 40 Ž1974. 215. w7x R.M. Davis, G.I. Taylor, The mechanics of large bubbles rising through extended liquids and through liquids in tubes, Proc. R. Soc. London, Ser. A 200 Ž1950. 375. w8x L.-S. Fan, K. Tsuchiya, Bubble Wake Dynamics in Liquids and Liquid–Solid Suspensions, Butterworth, Boston, 1990. w9x L.-S. Fan, C. Zhu, Principles of Gas–Solid Flows, Cambridge Univ. Press, New York, 1998. w10x J.M. Fox, Fisher–Tropsch reactor selection, Catal. Lett. 7 Ž1990. 281. w11x H.K. Henriksen, K. Ostergaard, Characteristics of large two-dimensional air bubbles in liquids and in three-phase fluidized beds, Chem. Eng. J. 7 Ž1974. 141. w12x C.W. Hirt, B.D. Nichols, Volume of fluid ŽVOF. method for the dynamics of free boundaries, J. Comp. Phys. 39 Ž1981. 201. w13x K. Idogawa, K. Ikeda, T. Fukuda, S. Morooka, Behavior of bubbles of the air–water system in a column under high pressure, Int. Chem. Eng. 26 Ž1986. 468. w14x B. Jager, R. Espinoza, Advances in low temperature Fischer–Tropsch synthesis, Catal. Today 23 Ž1995. 17. w15x R.-H. Jean, L.-S. Fan, Rise velocity and gas–liquid mass transfer of a single large bubble in liquids and liquid–solid fluidized beds, Chem. Eng. Sci. 45 Ž1990. 1057. w16x P. Jiang, D. Arters, L.-S. Fan, Pressure effects on the hydrodynamic behavior of gas–liquid–solid fluidized beds, Ind. Eng. Chem. Res. 31 Ž1992. 2322. w17x P. Jiang, X. Luo, T-.J. Lin, L.-S. Fan, High-pressure and high-temperature three-phase fluidization bed expansion phenomena, Powder Technol. 90 Ž1997. 103. w18x Y. Li, J. Zhang, L.-S. Fan, Numerical simulation of gas–liquid–solid fluidization systems using a combined CFD–DPM–VOF method: bubble wake behavior, Chem. Eng. Sci. 54 Ž1999. 5101. w19x X. Luo, P. Jiang, L.-S. Fan, High pressure three-phase fluidization: hydrodynamics and heat transfer, AIChE J. 43 Ž1997. 2432. w20x X. Luo, J. Zhang, K. Tsuchiya, L.S. Fan, On the rise velocity of bubbles in liquid–solid suspensions at elevated pressure and temperature, Chem. Eng. Sci. 52 Ž1997. 3693. w21x X. Luo, D.J. Lee, R. Lau, G. Yang, L.S. Fan, Maximum stable bubble size and gas holdup in high-pressure slurry bubble columns, AIChE J. 45 Ž1999. 665. w22x R. Mei, R.J. Adrian, Flow past a sphere with an oscillation in the

w23x

w24x w25x

w26x w27x w28x w29x w30x

w31x w32x

w33x

w34x

w35x

w36x

w37x

freestream and unsteady drag at finite Reynolds number, J. Fluid Mech. 237 Ž1992. 323. P.L. Mills, J.R. Turner, P.A. Ramachandran, M.P. Dudukovic, The Fischer–Tropsch synthesis in slurry bubble column reactors: analysis of reactor performance using the axial dispersion model, in: K.D.P. Nigam, A. Schumpe ŽEds.., Three-Phase Sparged Reactors, Gordon & Breach, 1996, p. 339. R.D. Mindlin, Compliance of elastic bodies in contact, Trans. ASME, J. Appl. Mech. 16 Ž1949. 259. X.D. Peng, B.A. Toseland, P.J.A. Tijm, Kinetic understanding of chemical synergy in LPDMEe process, 1998, Presented at the 1998 AIChE Annual Meeting, Nov. 15–20, Miami Beach, FL. P.N. Rowe, G.A. Henwood, Drag forces in a hydraulic model of a fluidized bed, Part I., Trans. Inst. Chem. Eng. 39 Ž1961. 43. S.C. Saxena, Bubble column reactors and Fischer–Tropsch synthesis, Catal. Rev. Sci. Eng. 37 Ž1995. 227. L. Schiller, A. Naumann, Uber die grundlegenden Berechnungen bei der Schwerkraftaufbereitung, Ver. Dtsch. Ing. 77 Ž1933. 318. B. Tarmy, M. Chang, C. Coulaloglou, P. Ponzi, Hydrodynamic characteristics of three-phase reactors, Chem. Eng. 407 Ž1984. 18. K. Tsuchiya, A. Furumoto, L.-S. Fan, J. Zhang, Suspension viscosity and bubble velocity in liquid–solid fluidized beds, Chem. Eng. Sci. 52 Ž1997. 3053. C.Y. Wen, Y.H. Yu, A generalized method for predicting the minimum fluidization velocity, AIChE J. 12 Ž1966. 610. P.M. Wilkinson, Physical aspects and scale-up of high-pressure bubble columns, PhD Thesis, Univ. of Groningen, Groningen, the Netherlands, 1991. B.H. Xu, A.B. Yu, Authors’ reply to the comments of B.P.B. Hoomans, J.A.M. Kuipers, W.J. Briels, and W.P.M van Swaaij, Chem. Eng. Sci. 53 Ž1998. 2646. J. Zhang, L.-S. Fan, C. Zhu, R. Pfeffer, D. Qi, Dynamic behavior of collinear collision of elastic spheres in viscous fluids, Proceedings of PTF Topical Conference at AIChE Annual Meeting, Nov. 15–20, Miami Beach, FL, Advanced Technologies for Particle Processing, Particle Technology Forum, AIChE Vol. 2 1998, pp. 44–49. J. Zhang, L.-S. Fan, C. Zhu, R. Pfeffer, D. Qi, Dynamic behavior of collinear collision of elastic spheres in viscous fluids, Powder Technol. 106 Ž1999. 98. J. Zhang, Y. Li, L.-S. Fan, Numerical simulation of gas–liquid–solid fluidization systems using a combined CFD–DPM–VOF method: single bubble rise behavior, Proceedings of PTF Topical Conference at AIChE Annual Meeting, Nov. 15–20, Miami Beach, FL, Advanced Technologies for Particle Processing, Particle Technology Forum, AIChE Vol. 2 1998, pp. 509–514. J. Zhang, Y. Li, L.-S. Fan, Numerical simulation of gas–liquid–solid fluidization systems using a combined CFD–DPM–VOF method: single bubble rise behavior, Powder Technol. Ž2000. in press.