Theory of two-phase cocurrent downflow in networks of passages

Ind. Eng. Chem. Res. 1991, 30,951-969

95 I

Theory of Two-Phase Cocurrent Downflow in Networks of Passages? Tomis R. Mellit and L. E. &riven* Department of Chemical Engineering and Materials Science, University of Minnesota, Minneapolis, Minnesota 55455

Soluble substances and heat are often transferred between liquid and gas on a commercial scale by means of coarse porous media such as packed beds. Two-phase cocurrent downflow in packed beds finds wide application in devices such as trickle-bed reactors. The first step in understanding the mechanisms of transfer is to understand how the two fluid phases are distributed and flow within the pore space, or “void”, of the packing. This paper shows that fluid distribution and flow in packed beds can be understood in terms of well-defined flow processes a t the level of individual passages and junctions. It reports further that, by considering interactions in the network, understanding can be formulated as a theory that explains observed flow regimes and transitions between them, and from which predictions of macroscopic behavior, including hysteresis and pulsing, can be computed.

In trod uc tion Countercurrent flow of gas and liquid in packed beds has been widely studied since early this century (White, 1935), owing to its applications in mass-transfer devices such as absorbers and distillation columns. Sherwood et al. (1938), on the basis of experiments with a glass-walled column filled with rings, described “flooding”,the buildup of water in the upper part of a packed bed because it cannot descend fast enough against the countercurrent flow of air; he also correlated pressure drop needed to drive the gas up the bed with the superficial gas and liquid velocities. Since Sherwood’s paper, much research has been addressed to understanding and predicting pressure drop, liquid holdup, and flooding in packed beds. Ergun (1952) developed a semiempirical correlation of pressure drop for single-phase flow. Today’s state-of-the-art methods of predicting these things are empirical and semiempirical correlations based on Ergun’s (1952) and Sherwood et al.’s (1938) landmark works. The study of two-phase cocurrent downflow parallels that of countercurrent flow. From the landmark paper on two-phase cocurrent downflow in packed beds of Larkins et al. (1961)until the beginning of the 19809, pressure drop and liquid holdup were correlated by using experimental data obtained with bench- or pilot-scale prototypes (Rao et al., 1983a; Kan and Greenfield, 1978). Ergun’s equation was the basis of correlations to predict pressure drop in two-phase cocurrent downflow (Rao et al., 1983b; Saez and Carbonell, 1985). More recently Saez et al. (1986) and Levec et al. (1986) modeled the hydrodynamics of two-phase flow in packed beds by representing the porous medium as an array of parallel conduits. This led to a relationship between the bed-scale pressure drop and the average interstitial velocities which modeled successfully the trickling regime. However this approach did not lead to predictions of the pulsing and bubbling regimes or the boundaries between them. Other efforts have been made to find a way to predict the transition regions between flow regimes, holdup, and pressure drop in trickle bed reactors (Talmor, 1977; Tosun, 1984). Scattering of data (Saez et al., 1986) and contradictory results (Kan and Greenfield, 1978)make the information that has been available inadequate to ‘This paper is dedicated to the memory of Professor Robert L. Pigford, Advisor and Advisor’s Advisor. *Aafim6iat x o v r c ( Gradwaoua1u aXXqXo1r“.

* Research Fellow of CONICET, Repfiblica Argentina. 0888-5885/91/2630-0951$02.50/0

understand the causes of different flow regimes. To this lack of understanding can be attributed the unreliability of existing design methods and the need for “safety factors”, Le., overdesign (Westerterp et al., 1984). Over the past decade, many authors have stressed the need for fundamental research in two-phase fluid dynamics of packed beds (Westerterp et al., 1984; the United States National Research Council, Panel on Natural Resources and Energy, 1986). The present paper is founded on the proposition that the causes of different macroscopic flow regimes are to be found at the level of packing particles or elements, in the passages where the gas and liquid flow interact.

A Modular Approach to a Complex Problem The theory reported here was built of elements that fit together as pieces in a puzzle. The approach was modular and proceeded in three successive stages: the formulation stage, the solution stage, and the interpretation stage; see Figure 1. In the formulation stage the theory is assembled from the basic elements. The macroscale (bed scale) and the microscale (particle scale) were connected by means of a statistical network theory. This approach-used for modeling widely different phenomena such as forest fires, diffusion in disordered media, and gelation processes (Stauffer, 1985)-relates the elements with the collection. Its successful application to capillarity-dominated multiphase displacement processes in porous sedimentary rocks (Heiba et al., 1984, 1986) was greatly encouraging. The “void space” in a packed bed is a collection of enlargements, or pore bodies, connected by passages that are necessarily constricted and so can be thought of as pore throats. The flow route map, or topology, of a packed bed, like that of almost any other porous media, is a network (Figure 2). The junctions are usually in the enlargements and can be represented as the nodes or sites; the flow passages, which contain the constrictions, are the branches or bonds connecting the sites. On a topological map so constituted the salient geometrical and compositional properties-the volumes of the sites and the cross sections of the constrictions-can be recorded as a roadmap maker would do. These properties are called the decoration of the network. Thus a decorated network model records the relevant topology, geometry, and composition of a packed bed. In the individual elements, i.e., the sites and bonds, the fluid phases that are present compete for space. The 0 1991 American Chemical Society

952 Ind. Eng. Chem. Res., Vol. 30, No. 5 , 1991

I

I

MASS CONSERVATION

+

EXPERIMENTS IN NETWORKS OF PASSAGES (MELLI et al. 1989)

I

(NETWORK SITES)

I1 It

I

GEOMETRY & TOPOLOGY OF THE PORE SPACE

(NETWORK BONDS)

I

i

EXPERIMENTS IN SINGLE PASSAGES (DE SANTOS et al. 1989)

I

CONSTITUTIVE LAWS & MASS

I

t

1

BOUNDARY CONDITIONS

STATISTICAL NETWORK THEORY

-------

FORMULATION STAGE

_ _ _ _ _ _ _ _ - - _ _ _ _ .

\ THEORETICAL MODELLING (MELLI 1989)

LINEARIZATION ALLOWABILITY CRITERIA

NEWTON-RAPHSON ITERATIVE SOLUTION

1

CRITERIA

SOLVER J

TIME STEP UPDATING PREDICTION-CORRECTION STRATEGY

*

FINE P---OKOUS MEDIA

v MATHEMATICAL SOLUTION STRATEGY

SOLUTION STAGE _ __--------------------------I

DISCRETE FAST FOURIER TRANSFORM

PULSE FREQUENCY ANALYSIS

competition takes place through the action of several forces, namely, those of gravity, viscosity, inertia, and capillarity (those of wetting are put aside here, the liquid being regarded as wetting the solid surfaces), which generate different local flow regimes, depending on the state of occupancy and flow within an element and its neighbors. The flow regimes in individual passages have been found experimentally and described by de Santos et al. (1989). The same flow regimes were found by Melli (1989) in the passages of small model packed beds, Le., in networks of passages. Different states of sites were also found experimentally and classified in terms of the relative amounts of gas and liquid present (Melli, 1989). The mass conservation equations as written for sites and passages relate the gas and liquid flow rates and the rate of change of liquid saturation (the fraction of the total volume occupied by liquid). The mass conservation equations together with constitutive equations for the flow rates generate a differential algebraic system of equations for which boundary and initial conditions are needed. Commonly in practice the boundary condition at the bottom of a packed bed is, approximately at least, a uni-

GRAPHIC TECHNICS

TRICKLING REGIME HYSTERESIS ANALYSIS

form and constant pressure in the great enlargement to which all of the last passages lead, and at the top, again approximately at least, fixed total gas flow and totalliquid flow leading into the passages there. So the boundary conditions of the network were drawn from the common practice: the same constant pressure in all the sites at the bottom of the network. The initial conditions were drawn from distributions of pressures and states of sites and bonds or from the final states of a previously calculated prediction. The solution stage is composed of two main parts, as follow. The first part was a strategy to follow the evolution of states of occupancy and flows in the elements-the sites and the passages. The theory of quasi-steady-state multiphase flow in fine porous media developed by Heiba et al. (1986) was based on the postulate that the fluid phase that occupies a given segment of a pore space is governed by definite allowability and accessibility rules. A similar postulate was used in the present theory. There are four basic flow regimes in single passages (de Santos et al., 1989): gas-continuous, bridged, flooded, and bubbling. An individual passage performs in one or an-

Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991 953

* NETWORK

RELATES MACROSCALE TO MICROSCALE.

*S E A C C O U N T MICROSCALE :FLOW THROUGH PASSAGES, OR BONDS; HOLD-UP IN ENLARGEMENTS, OR SITES.

FOR HOLD-UP CAPACITIES.

BONDS FOR FLQW CONDUCTANCE, BOTH DRAWN FROM REALISTIC DISTRIBUTIONS.

Figure 2. Void space and flow distribution in a packed bed.

other of these flow regimes when certain rules are satisfied. For example, a single passage performs in the gas-continuous regime when the liquid flow rate is lower than a maximum value given by the geometry of the passage, properties of the fluid phases, and gas flow rate. When these rules are satisfied, a passage in the network is allowed to perform in the gas-continuous regime, though it may not actually do so, as explained shortly. The particular rules for each of the four flow regimes, described theoretically and experimentally by de Santos et al. (1989) and Melli (1989), constitute the allowability criteria of the theory. Allowability criteria are not enough to determine the current flow regimes in the passages of the network. The history of fluid migrations in the network sets the states of occupancy in the sites, and those states sway the flow regimes in adjacent passages. For example, the gas-continuous regime is precluded from a passage if its neighbor sites are occupied by liquid. The rules dictated by the state of occupancy are the accessibility criteria of the theory. A passage in the network performs in a particular flow regime if appropriate accessibility and allowa bility rules are satisfied. The second part of the solution stage was actually solving the differential algebraic system of nonlinear equations. The backward Euler method is used to integrate in time, and Newton’s method is used to solve iteratively the nonlinear algebraic equations at each time step. In each iteration the linearized system, which produces a large, sparse, banded, nonsymmetric matrix, is solved with a frontal solver. The time step used proved critical because states frequently change in the elements and countless iterations would have been needed to predict the much lower frequency oscillatory states at the macroscale. A predictor-corrector strategy was employed to select the largest feasible time step that still captured all the changes in states of individual sites and bonds. The interpretation stage was the analysis of results. Graphic displays (computer graphics) were used to illustrate effects of operational history, e.g., hysteresis in the trickling regime. Differences in passage and site occupancy stood revealed as the germs of differences observed at the macroscale, i.e., in overall network behavior. Oscillatory states showed characteristic oscillatory frequencies found in packed-bed experiments by Talmor (1977) and in sim-

pler networks of passages by Melli et al. (1990). The fast Fourier transform was applied to the time-dependent outputs in order to study the effects of the main variables-overall gas and liquid flow rates-on the pulsing frequency in the pulsing regime.

The Network Model To describe two-phase flow in a packed bed completely would require the velocity and pressure fields in each phase and the meniscus or interface position between the phases everywhere in the void space. Because of the extent and shape of the void space, this is impractical. Moreover, the description of two-phase flow in single constricted passages is still an area of intensive research (de Santos et al., 1989). The crucial aspects of packed-bed behavior appear to be insensitive to local details of the velocity and pressure fields. However, knowing how rates of gas and liquid flow through an individual passage vary with the amount of liquid present and the pressure difference across the passage should suffice for describing two-phase flow in a packed bed adequately. Thus what are needed are constitutive equations that relate flow rates to occupancies and pressure differences. The pressure field can be approximated by the pressures in the sites, viz., the pressure of the gas in gas-dominated sites, of liquid in liquid-dominated sites, and the pressure of both phases in contested sites. This approach leads to a connection between the pore scale and the macroscale, which is detailed enough that macroflows can be related to microscale mechanisms; yet it is simple enough to allow mathematical analysis. Visualization of flow in a network of well-wetted passages revealed that the cocurrent downflow regimes, in wetting situations, at least, are statistical combinations of the pore-scale flow regimes (Melli et al., 1990). The macroscale regimes that were found-trickling, spray, bubbling, and pulsing-are also found in coarse porous media such as trickle-bed reactors. The statistical network theory presented below accounts for steady, transient, and oscillatory regimes in a network of passages and provides understanding of the basic mechanisms that generate these regimes. In well-wetted individual bonds, the basic flow regimes are gas-continuous or trickling, in which the liquid flows as a film on the wall and the gas flows unobstructed through the core of the passage; bridged, in which the

954 Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991 TYPE OF SITES

A FRACTION OF INCOMING LIQUID ACCUMULATES WHEN THE PENDANT DROP IS SMALL

GASDOMINATED SITE CONTESTED SITE LOW LIQUID HOLDUP INTERMEDIATE LIQUID HOLD-UP

LIQUID-DOMINATED SITE HIGH LIQUID H O L D U P

TYPE OF BONDS

a GAS-CONTINUOUS

b BRIDGED

c FLOODED

d BUBBLING ALL LIQUID LEADING INTO THE SITE FLOWS OUT DROP SUE IS CONSTANT WHEN THE PENDANT DROP IS BIG

A COMBINATION OF REGIMES a AND b GIVES THE LOCAL PULSING REGIME

Figure 3. Local flow regimes in the passages and types of sites.

-~~ [ ~ ~ ~) ~ ~ " , %

6 A S FLOW I N BONDS

GAS FLOW

f

VoLGE

d(

Figure 6. Liquid accumulation alternatives in a site.

1

CONSTITUTIVE

LIOUIDFLOW INTO

EQUATION

\

I a - Bubble-generating bond

b - Bubble-passing bond

Figure 7. Types of bonds in bubbling regime. VOLUME

I

Figure 4. Variables and equations in a site and adjoining passages.

make a bioconical unit. Such a unit is characterized by its throat radius r,, body radius rb, length lo, and inclination from the flow direction 8. The capillary, gravity, viscous, and inertia forces are related to the relative magnitudes of rb/rt, io/rt, 8 and play the main role in setting the limits between the four basic flow regimes. In order to mimic the more-or-less random distribution of pore sizes and shapes, the throat radius r, and the passage inclination 8 are assigned at random according to normalized distribution functions a(r)and p(O), where a(r) is the probability that a passage has a pore throat radius between r and r + dr and @(e) is the probability that it has an inclination between 8 and 8 de. As a convenient approximation, rb and 1, can be regarded as perfectly coordinated with r so that rb = rb (r);and lo = lo (r) (cf. Heiba et al., 1984). The adjustable functions of the theory, a(r) and p(S), can be used to describe the particular geometry of a packed bed or porous medium. The site is gas-dominated when the gas occupies almost entirely the site space and the liquid is present as films and sometimes also as pendant drops from the solid surfaces. The site is liquid-dominated when the liquid fills entirely the site space and the gas, if present, is dispersed as bubbles. The site is contested when the gas and liquid holdups are comparable (see Figure 3). The sites account for the holdup capacity. The site volume influences the pulse frequency at the macroscale; bigger sites require longer times to be filled and emptied by an invading phase, whether it is liquid or gas. The ratio between the volume of a site and the volume of the passages connected to it influences the pulsing frequency in the passages; larger site volume/passage volume ratios give smaller passage pulsing frequencies. To mimic the random distribution of site volumes in the network, the site volumes are assigned a t random according to a normalized volume distribution function V), which becomes another adjustable function of the theory. In the network to which the theory was applied the average throat radius was r, =

+

a

-

Liquid is distributed evenly I n both outflow bonds l a n d 2

c

-

Bond 2 Is p r e f e r r e d o u t f l o w o f t h e Ilquid

Bond 2 carrier almost all the liquid flowing out o f the site

Figure 5. Liquid distribution to the bonds leading out of a site.

liquid clogs the constriction in the passage and the gas is accessible to one or both ends of the passage; flooded, in which the liquid f i i completely the passage; and bubbling, in which the liquid flows as a continuous phase and in it the gas is dispersed as bubbles (see Figure 3). Each flow regime has its own relation of gas and liquid fluxes to occupancy and pressure difference, i.e., its own constitutive equation. Each flow regime arises in a particular range of occupancy and of gas and liquid flow rates. The limiis between flow regimes are set by the occupancy, the gas and liquid flow rates, the fluid properties, and the size and shape of the passage. A convenient approximation to the geometry of the passage is a converging cone fused to a diverging cone to

r(

Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991 955

(a) Network Model. Numbers in plain characters identify passages; numbers in italic characters identify sites.

(b) Mesh Model. Numbers in plain characters identify nodes with four variables (passages); numbers in italic characters identify nodes with two variables (sites).

Figure 8. Network model and mesh model of the network.

1 mm and the average ratio of site volume to passage volume was 15. In both cases a truncated Gaussian distribution around the average was used to set the throat radius and the site volume in each individual site and passage. In this way the network of passages and sites can be laid out. The next steps are to identify appropriate constitutive equations for the local flow regimes within the network, allowability criteria that set the limits on those regimes, and accessibility rules about when an allowable local regime can actually arise. These are elements of any network model, including assemblages of spheres (Zimmerman e t al., 1987) and cylinders (Chu and Ng, 1989) employed to simulate trickling cocurrent downflow. Zimmerman et al. (1987) employed a two-dimensional network created by spheres in a four-coordinated diamond array. The allowability criterion was that when the sum of liquid flow rates arriving at a site fell below a critical value, the outflow passages were taken to be only partially wetted. In this case the left outflow was made equal to the left inflow and the right outflow to the right inflow; when the total flow exceeded the critical value half of it was assigned to each outflow. This constitutive relation of site action is called a splitting rule: use of an entire distribution of splitting rules is described below. Chu and Ng (1989) used a three-dimensional, six-coordinated, simple cubic network of cylindrical channels. They considered two gas-continuousflow regimes: annular flow and segregated flow. To each were assigned relations between pressure drop, holdup, and flow rates. In the junctions or sites no mass accumulation was allowed. The allowability criterion was simply that when overall liquid rate was raised, all channels performed in the segregated regime and when overall liquid rate was lowered, they all performed in the annular regime. (From this criterion sprang hysteresis loops of pressure drop versus liquid rate at fixed gas rates in the trickling regime.) The statistical network theory presented here incorporates to a much greater extent the realities exposed by laboratory studies of the microflow behavior in individual passages and sites (de Santos et al., 1989) and in small networks of passages (Melli et al., 1990). All the local flow regimes, gas-discontinuous as well as gas-continuous, are taken into account. Each is modeled by a constitutive equation drawn from theory and experiment with representative single passages; allowability criteria are likewise founded. Simple analysis of occupancy states are allowed to fill and empty with liquid and gas, a key figure of the

l9

L

LON ST ANT GAS REYNOLDS NUMBER

0

A t

01 0

I

I zoo0

I

Raising liquid flow raw Reducing liquid flow raw hRr,,= 3090 Reducing liquid flow rate from R e = 1717 Reducing liquid flow rate from R%= 687 Reducing liquid flow raw from Re,,= 515

I 2000

I

I 3000

I

4000

2 LPL

LIQUID REYNOLDS NUMBER,

Rk= Z r t N P P L

Figure 9. Hysteresis in trickling regime in the 6 by 10 network uniform distribution at the top of the network.

pulsing regime, as confirmed by the laboratory studies. What bears emphasizing is what visualizations by means of magnifying lenses and high-speed video cameras revealed: that the various combinations of local regimes predicted by the theory can indeed occur in packed beds-small, nearly two-dimensional ones at least (Melli, 1989). The statistical network theory is formulated as a set of governing equations and boundary conditions.

Governing Equations Figure 4 shows a site and adjacent bonds in a four-coordinated network of passages. Two of the bonds lead in and two lead out of each site. In its present state, the theory is for isothermal systems without mass transfer between phases and without chemical reactions. Thus, momentum and mass balance equations suffice to describe the behavior fully. Pressure drops are assigned to the passages and not the sites, on the grounds that sites are

956 Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991

"I

,

, mono TIZ~HYMtlONLEJS)

0

0

$1

I

Moo0

0

, 8001

, low0

nw (DIWWJIONLLSS)

I

-0

;I 0

,

,

,

6000

loo00

IWOO

TlME (DIMENSIONLESS)

a)Low liquid load. Evolution from Re,= 1029 to Re L= 1374. Output liquid flow and pressure difference climb steadily. Output gas rate first falls and later grows up to the steady state value.

-d l

I

moo0 nuE (DIUENSIONUSS)

0

BO00

16000

0

I

6000

(0001

ntaE (DIMENSIONIXSS)

la000

6000

0

lO0OO

woo

TlME (DIMENSIONLESS)

b)High liquid load. Evolution from ReL=2400 to Re,=2745. Formation and grow of liquid-rich clusters change gas and liquid distribution and induce oscillations in pressure difference and output flow rates during transient evolution

Figure 10. Time evolution of pressure difference and output liquid and gas flow rates in a network of 6 by 10 sites, as liquid flow rate to the distribution manifold is raised. Gas flow rate to the manifold is constant such that Rec = 115. Units of measurement are given in the text.

represented by junction points and that in reality the pressure varies much less across enlargements than along the passages. The liquid phase is taken to be incompressible. Compressibility of the gas phase is accepted in the distribution manifold and in the sites of the network, but in the passages the gas flow is considered incompressible. Mass Balance in the Sites. The variables at each site-holdup, or liquid accumulation, and pressure-are calculated from the liquid mass balance and the overall mass balance in the site. The mass balance equations for each of the phases in a site are

Equations 1can be written in terms of the site variables, pressure and liquid accumulation (in dimensionless form): d In Ps (VS - VL) -- e ( ~ iG') dt i=l

+

-d -V -~ 5 L i dt

i=l

The units of measurements are flow rate

2?r(PL - PG)r;g Q 3

PL

v

volume time

vL

Here is the volume of liquid accumulated in the site, L' is the volumetric liquid flow rate in the ith passage, z is the site coordination number, i+, the number of passages connecting there, MG = PG,VG is the mass of gas accumulated in the site, _Ft; = pGG,is the mass flow rate of gas in the ith passage, t is time, G is the gas volumetric flow rate in ths ith passage, PG = (WPG)/(RT) is the ideal gas density, T is the a_bsolute temperature, W the gas molecular weight, and VG the volyme o,f gas accumulated in the site. The site volume is Vs = VL VG.

+

pressure

t=

3

2*r,3 PL

(PL -

(PL

PG)gro

- PG)glo

Here ro is the pore constriction radius, + is the liquid density, pG is the gas density, and pL is the liquid viscosity. Constitutive Equations and Mass Balance in the Passages. In steady state, the two-phase cocurrent flow in constricted passages is fully described by two constitutive equations that relate the gas and liquid flows to the volume average saturation and with the potential differ-

Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991 957

b)Intermediale liquid flow rate, ReL= 1372, reached from a lower flow rate. Small clusters of liquid- dominated sites appear

a) Low liquid flow rate, ReL= 516. Sites are dominated by gas, Flow regimes in passages are gas-continuous and bridged

~-~

~-

c) ReL= 1372, reached from a higher flow rate. Liquid dominated clusters are bigger than those at the same rates reached raising the liquid flow rate

-

d> ReL= 1372, centered liquid distribution. Liquid cluster concentrates in the center. Gas channeling is noticeahle in the sides of the network

Figure 11. Snapshots of steady states in trickling regime: effect of liquid distribution and history of operation. Gas Reynolds number ReG = 115.

ence between adjacent sites. On the basis of experiment and theory, these equations were developed by Melli (1989) for the four basic flow regimes-gas-continuous, bridged, flooded, and bubbling. In time-dependent states, the flows, pressure difference, and average saturation in the passages are function of time. The liquid flow rate leading out of a passage, L, adjusts, with a delay given by the passage holdup capacity, to changes in the liquid flow rate leading into the passage, Li,. The gas flow rate throughout the passage is taken as adjusting instantaneously to changes in the adjacent sites. The average passage saturation, s, is given by the mass conservation equation:

-1, -ds-- Lin - Lout ro dt

(3)

Here, 1, is the passage’s length and rt its throat radius, as shown in Figure 4. In this approximation, four equations are needed in each passage: the mass conservation equation (3) for saturation, a constitutive relation for the liquid flow rate into and another for the flow rate out of the passage, and a constitutive equation for the gas flow rate. The equations chosen for the four flow regimes in the passages are shown in Table I, summarized in the following paragraphs, and explained in the Appendix. In the gas-continuous regime the liquid phase forms a film and flows by gravity, whereas the gas phase is driven by the pressure difference across the passage and by the liquid drag at the gas-liquid interface. The gas-continuous

regime prevails at low gas and liquid flow rates; therefore, the constitutive equations for gas and liquid flow rates are approximated by those of viscous flow in biconical tubes where a volumetric average saturation (see analysis in the Appendix) and average radius account for the particular shape of the passage (see eqs 4 and 5 in Table I and eqs A-3 and A-4 in the Appendix). In the bridged regime, liquid forms a bridge in the passage constriction; the gas phase occupies one or both ends of the passage. Capillary, viscous, and gravity forces set the position of the gas-liquid meniscus in a way detailed elsewhere (Melli, 1989) and outlined in the Appendix. The constitutive laws are shown in eqs 6 and 7 in Table I. In the flooded regime, which is the limit of bridged and gas-continuous regimes when the liquid saturation in the passage is unity, the liquid occupies the entire passage. The constitutive law (eq 8 in Table I) is that of potential-gradient-driven flow in a passage. The bubbling regime appears at moderate and high liquid flow rates in all the ranges of gas flow rate. The constitutive equations in the bubbling regime are approximated by relating the potential difference to a quadratic function of the overall flow rate, and adopting the drift velocity relation (Ishii, 1977) between the gas and liquid flow rates. In terms of the characteristic units the constitutive equations are (9) and (10) of Table I. The liquid that flows into each passage is a fraction of the liquid flowing out of the site that feeds that passage. The liquid flow out of a site, or pore body, divides over

958 Ind. Eng. Chem. Res., Vol. 30, No. 5 , 1991 h

2

Y

h

m

v

h

n

2

13 v

Y

2 ,a-'1

2

m g

el

ew;

gw 2

I

I

$

el

II

II

a

2

el

-1

n

s

w

n

n

a

c

h

3

I

2

+

-1

h YI

I

vv

b

44

el II

n

h

c

h

2

2

v

Y

.

n

r

s

el

5

3

-1

n

a

II

0 II

u

0

0

-.

Q, .

?-

?-

A

A

4

4

Q,

Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991 959 the passages available, and does so according to local sizes, shapes, and orientations that vary more or less chaotically from place to place (see Figure 5 ) , so that there is a distribution of divisions over the network. To approximate the liquid distribution, the network is decorated with a random distribution of liquid split coefficients, 1,so that the liquid flow leading into each passage is a fraction of the total liquid flow leading out of the site. When two passages lead out of a site, one of them carries a fraction q of the liquid and the other a fraction 1 - q. Therefore, the sum of the split coefficients of the passages leading out of each site is unity, i.e., = 1. A fraction of the liquid flow arriving to a gas-dominated site, p, accumulates as a drop hanging from the upper solid surface (see Figure 6a); the rest of the liquid flows out. The size of the pendant drop grows as the liquid accumulates. When the drop reaches a maximum volume (found by experiment to be about 40% of the site volume), liquid drips from the drop, as Figure 6b shows, and all the liquid leading into the site flows out of it. For the examples when liquid shown here, the function cp is equal to accumulates in the drop and cp is equal to 1when the drop reaches its maximum volume. Equation 14 relates the liquid flow rate leading into each passage with the liquid flow rate arriving to the site and with the functions q and cp.

x$t+

gas arriving at the site is distributed in the bubbling passages by a distribution law:

Here 4 is the liquid velocity in the j t h passage, G = Gin + dVG/dt is the total gas flow rate out of a site, and nbb is the total number of outlet passages performing in the bubbling regime. The minimum liquid flow rate able to entrain bubbles can be approximated as L = y/9 (see Melli (1989)), where y is the ratio of the bubble radius to the throat radius of the passage. Accordingly, the function @ vanishes at liquid flow rates below the limit and is unity above it. Boundary Conditions. At the bottom of the network the pressure in the great enlargement to which all of the last passages lead is fixed. At the top the total gas and liquid flow rates leading into the passages there are fixed. With all quantities measured in the units listed above, the set of equations for the distribution manifold (top of the bed) becomes qn= F'L *m

C@=1

j=l

212

Lin

= +ELfn,siteVsite i=l

d In PG

(14)

Here Lin,siteis the liquid flow leading into the site through the ith passage. Flooded passages are unable to carry all the liquid flow given by eq 14; i.e., the liquid flow rate into a flooded passage is equal to the liquid flow rate out of the passage, and there is an overflow of liquid that has to be distributed among the nonflooded passages. When both flooded and nonflooded passages coexist at a site, the constitutive law (14) for nonflooded passages has to be modified to account for the overflow of liquid arriving from flooded passages. The modification is (12) of Table I. In eq 12 noois the number of outlet passages not flooded and nor is the number of outlet passages flooded. In a fourcoordinated network the maximum value of nooand nofis 2. In the bubbling regime two cases arise. The first, called bubble generation (Figure 7a), appears when the gas phase is accessible to the passage. In this case, if the liquid flow rate is high enough to preclude the gas-continuous regime, gas is dragged as bubbles dispersed in the continuous liquid phase; the bubbles are generated at the passage entrance. The liquid flow into the passage is given by eq 12; the two steady-state constitutive laws (eq 9 and 10) and the liquid mass balance (eq 3) complete the set of four equations needed. It bears emphasizing that, in visualizations of the nearly two-dimensional packed bed, bubble generation was encountered, in different instances, in sites at the bottom, middle, and top of the network. At high liquid flow rates it was the most common case in the sites connected to the distribution manifold. The second case, called bubble passing (Figure 7b), appears when the neighbor sites are dominated by liquid. In this case bubbles cannot be generated at the passage entrance because gas is not accessible as continuous phase. Nevertheless, the bubbles arriving at the site can be carried out of it by the flowing liquid and dragged along through the passage. Because the passage is full of liquid excepting the bubbles, which are taken to be incompressible, the inlet and outlet liquid flows are regarded as the same. The flow rate of gas entrained by the liquid is taken to be proportional to the liquid velocity in the passage; therefore the

(VC

-V k ) T

Po

- -G PG

nm

- L - C(yn- Gin) (16) j= 1

Here, nbtis the total number of passages leading out of the distribution manifold, qnand Lfnare the gas and liquid flow rates in the jth passage, G and L are the gas and liquid flow rates delivered to the network, 9 is the split coefficient for the j t h passage, Vc is the volume of the distribution manifold, Vk is the volume of liquid accumulated in the distribution manifold, and Pc is the pressure in the distribution manifold. The initial condition was the same pressure in all sites, and all sites and passages filled with gas, Le., no liquid present in the network.

Allowability and Occupancy Statistics Central to the statistical theory are the laws that govern the state of each site and flow passage. The occupancy of sites is simply given by the immediate past history of the operation: gas-dominated, liquid-dominated, and contested sites are given by the relative amounts of gas and liquid phases present in each site. Allowability rules in the flow passages are given by (1) the number of phases present in the passage, (2) the size and shape of the passage, (3) fluid properties, and (4) relative liquid and gas flow rates. Occupancy statistics are given by the accessibility of the phases to the neighbor sites, which is determined by the history of fluid flow in the pore space. The limits of the basic flow regimes are given in Table 11; they are described in detail by Melli (1989) and summarized in the Appendix. Outside of its limits a flow regime is either impossible or unstable and evolves into one or another new regime. For example, the gas-continuous regime is allowable at low liquid flow rates. But as the liquid rate in the passage is raised (the gas rate being held constant), the liquid film on the passage wall thickens until it collapses and liquid bridges the narrowest part, the throat. Then the gas-continuous regime is no longer allowable in the passage. A passage is considered to have

960 Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991

]

become allowable to a given flow regime when the following rule is satisfied: -bond allowablj if

previously in regime j and its new state violates no limit of that regime

to flow

regime j

or

i 1 bond previously in another regime and its new state falls within the limits of regime j

This rule summarizes all the limits given in Table I1 and explained in the Appendix; the limits in turn come from observations of flow regime transitions in single constricted tubes (cf. de Santos et al. (1989)). Allowability criteria suffice for single constricted tubes, but not for passages of a network. Some of the passages allowed to a given flow regime may not be performing in this flow regime. To perform in a particular flow regime, the passage should satisfy allowability and accessibility criteria. The accessibility rules depend on both the flow regime and the occupancy of neighboring sites. These rules are passage accessible to] gas-continuous and bridged regimes

if

passage accessible to bubbling regime

if

[

[

[

1

either site connected to the passage is accessible to gas-continuous phase

I[

[

1

either site connected to the passage is accessible to dispersed or gas-continuous phase

I [ =

1

always (regardless of the state of connected sites)

The number of combinations of flow regimes in the outlet passages that satisfy the accessibility and allowability requirements depends on the number of basic flow regimes and the coordination number of the site. In general, if p is the number of possible flow regimes in the passages and z, is the number of outlet passages at each site, then the number of alternative combinations is

number of alternatives =

The nonlinear algebraic system of equations that describes the system consists of mass balance equations, constitutive laws, and boundary and initial conditions and is of the form f(Y,Y') = 0 (18) Here y is the vector of unknowns: P and V , in sites, G, Li,, Lout,and s in the passages; y' is the time derivative of y. This differential algebraic system can be written in the form

M

dY

= g(Y)

(19)

The matrix M, sometimes called the "mass matrix", consists of the coefficients of the time-dependent equations (holdup capacity in the mass balance equations). The structure of g(y) is given by the mass balance equations in sites and passages and the constitutive laws for the currently performing flow regimes. M is a singular matrix and any rational, implicit, one-step or multi-step difference scheme can be used for integrating (19). A good choice is the first-order backward difference formula: (Y"+' - Y") = dY"+') At

(20)

In each time step y" is a known vector (equal to yoor initial condition at t = 0), so that (20) can be rewritten as

(If accessible to continuous gas phase, the passage can perform in the bubble generation regime; if accessible to dispersed gas phase only, the bubble generation regime is precluded and the passage can perform in the bubble passing regime.) passage accessible to flooded regime

Mathematical Model and Solution

(17)

There are 10 combinations of four flow regimes @ = 4) in a four-coordinated network with two outlet passages at each site (in a six-coordinated network with three outlets per site the number of combinations is 20 and in a eight-coordinated network with four outlets per site the number of combinations is 35). The combinations for any network with two outlet passages per site are as follows (passage 1, passage 2): (1)gas-continuous,gas-continuous; (2) gas-continuous, bridged; (3) gas-continuous, bubbling; (4) gas-continuous, flooded; (5) bridged, bridged; (6) bridged, bubbling; (7) bridged, flooded; (8) bubbling, bubbling; (9) bubbling, flooded; (10) flooded, flooded. All of these have been identified in visualizations (Melli et al., 1990).

where b = (M/At)y"is known from the previous time step. Equation 21 is nonlinear in yn+l. It can be written h(y"+') = 0, where h(y"+') g - (M/At)yn+'- b, and is best solved by Newton iteration:

where yn+'J"l is the vector of variables in the mth iteration of the ( n + 1)st time step and J (dh/dy)y=yn-l,m is the Jacobian matrix in the mth iteration of the ( n + 1)st time step. The first estimate of the vector of variables (y"+',O) is the vector y". The solution of the algebraic and differential system of (21) and (22) poses several challenges. The Jacobian matrix J is large, nonsymmetric, banded, and sparse, and its structure depends on time because the constitutive laws change when the flow regimes in individual passages change. For example, in a 10 X 24 network such as that used by Melli et al. (1989), there are 2229 variables, and in any row or column of the Jacobian matrix there are no more than six nonzero entries. The changes in constitutive equations of flows in passages are sharp and alter the structure of the matrix. For example, a gas-continuous passage near the maximum allowed saturation carries both phases simultaneously, whereas immediately above the maximum saturation the gas flow drops to zero. The system of (22) can be solved with band solvers, such as LINPACK, sparse solvers, such as Yale Sparse Matrix Package, or frontal solvers, such as Hood's. A frontal solver is the best choice because only a part of the Jacobian is in the central memory at any time, thus reducing the size of memory needed. To use it, the network is transformed into a mesh (see Figure 8), the sites and bonds of the network mapping onto mesh nodes which serve to order the equations for the solver. Some nodes have two variables (P and V , in the sites) and other nodes

:nd. Eng. Chem. Res., Vol. 30, No. 5, 1991 96 1 CONSTANT GAS REYNOLDS NUMBER

~~

CONSTANT GAS REYNOLDS NUMBER

Uniform liquid disEibution: Raising liquid flow rate 0 Reducing liquid flow rate

Centered liquid dishibution: Raising liquid flow rate 0 Reducing liquid flow rate from Re = 3090 0 Reducing liquid flow rate from Re:= 1717 A Reducing liquid flow rate from Re,= 687 ReducingliquidflowratefromReL=515

- .

Raising lisuid flow rate Reducing liquid flow rate

+

1000

0

2000

3000

LIQUID REYNOLDS NUMBER, R%= 0

I000

LIQUID REYNOLDS NUMBER, R e L s

2 LPL xr,NPPL

Figure 13. Hysteresis in trickling regime: effect of liquid distribution in a network of 6 by 10 sites.

2000

2 LPL r(

4000

NP P L

Figure 12. Hysteresis in trickling regime: scanning loops in a network of 6 by 10 sites.

have four variables each (s,Lb,Lout,and G in the passages). The solution gives the overall pressure drop in the network, the total holdup, and the state of each bond and site as well as the local flow regimes and saturations in sites and passages. The strategy used to find the time step is based upon the estimation of the earliest time at which a change in the Jacobian structure is expected. At the end of a time step the allowability and accessibility rules are applied to all the passages, and the constitutive equations are changed in the passages where the rules are violated. The structure of the function g changes. The predicted changes in the vector of variables in the ( n 1)st time step, Ag = gn+l- y", is calculated without iteration, by means of forward Euler difference approximation

+

(23) Here At is the tentative time step, selected as the arithmetic average of the last three time steps. The function g(yn)consists of the mass balance equations and constitutive equations of the flow regimes in the passages that will perform during the next time step, evaluated with the results of the current time step. The actual time step for the ( n + 1)st iteration is calculated from

where AyTaXis the maximum change allowed the jth variable. Frequency Analysis of the Pulsing Regime On the bed scale in the pulsing regime, slugs or pulses of enhanced liquid holdup alternate with slugs of enhanced

gas holdup. The frequency of alternation tends to be regular. Talmor (1977) and Christensen et al. (19861, among others, studied the characteristic frequencies in industrial- and pilot-scale trickle-bed reactors. Kolb et al. (1990) used acoustic frequency analysis to probe the relation of pulse frequency to gas and liquid flow rates in a network of passages. At any fixed place in a packed bed the pressure oscillates as slugs of entrained liquid holdup and enhanced gas holdup pass by. Such pressure oscillations are also predicted by the statistical network theory. The discrete Fourier transform applied to these results gives the characteristic frequency of each oscillatory state. Discrete Fast Fourier Transform. The Fourier transform of a continuous time-dependent function f(t) is

X ( m ) = Smf(t)ci2lmt -m dt

(25)

For a discrete function sn(t)of N sample points, for example the pressures predicted by the statistical theory at fixed places of a network of passages, the integral (25) is replaced by a summation approximation, where the coefficients are N

C2m-2

=

cs, cos n=l N

=

Cs, sin n=l

[

[

(m-1)(;-

1) 27r

I;

m = 2, ...,N / 2

( m - 1);

- 1)27r

; m = 2,

+1

...,N / 2

The power spectrum S(n)= ( c ~ +~czm-12)/N - ~ ~ of the discrete Fourier transform can be used to compare the relative magnitudes of different frequencies and to determine the characteristic frequencies of the pulses. The minimum record length, t,, and the minimum number of points, N , required in the sample are set by the highest

962 Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991

(a)First liquid pulse in the center, second liquid pulse appears at the top, t = 0 msec

(b)Roth pulses travel down the network, t = 9 msec

(cpirst pulse reaches the network exit, t = 13.5 msec

4 4 @I . .

(d)Rotlom: l a i l of first ~ I I I S C , Top: Third pulse forms, t = 18 msec

(e)Rottom: Second pulse reaches the network exit, t = 27 msec

(f)Bottom: Tail of the second pulse, t = 36 msec

Figure 14. Pulsing regime sequence predicted by the theory in a network of 3 by 10 sites. Gas Reynolds number ReG = 3166; liquid Reynolc number ReL = 3000.

frequency that is desired, Viz., fh, and the minimum desired frequency resolution, F (Stanley et al., 1984): t, = 1/F

De Santos et al. (1989) found, experimentally, that the pulses in passages have high-frequency oscillations, on the order of hundreds of hertz, and these set the highest frequency that needs to be detected. Melli et al. (1990) found that pulses at the bed scale have low-frequency oscillations, on the order of hertz, that set the lowest frequency that needs to be resolved. In order to analyze bed-scale pulses with the statistical theory, long samples with numerous points are needed. For example, if fh = 200 Hz and F = 0.1 Hz, the minimum record length is t, = 10 sec and the minimum number of points is N = 4000. Examples of small networks of passages are illustrated in the following sections. Two-Phase Flow Examples Trickling Regime. When the liquid flows down the network at low rates, it moves as films along the passage walls while gas flows down the center. Occasional sites and passages may be filled with liquid; there are clusters of liquid-dominated passages and sites. As the liquid flow

3

rate is raised, liquid clusters form and grow, and the di tribution of gas and liquid changes. When the liquid flow rate is lowered, clusters shrink through drainage of sites in a different sequence than that in which they were formed. More than one steady state at a given gas rate and given liquid rate is common; in other words, there is hysteresis. Hysteresis in the trickling regime was first reported by Kan and Greenfield (1978) in a 10-cm-diameter column filled with glass beads 2 mm in diameter. Christensen et al. (1986) reported hysteresis in the trickling regime when the liquid flow rate was raised and lowered at constant gas flow rate in a packed bed of square cross section. Figure 9 shows the theoretical prediction of the pressure difference at a constant gas flow rate, corresponding to a gas Reynolds number Rec = 115. The liquid flow rate is raised from zero to a maximum corresponding to a liquid Reynolds number ReL = 3090, and subsequently lowered to zero in a network 6 sites wide and 12 sites long, in which the average pore radius is r,, = 1mm. Here the Reynolds numbers are defined as Re, = 2Q,p,/?rr,AvNPp,, where pa is the density, pa is the viscosity, Q, is the volumetric flow rate of phase CY, and NP the total number of passages in each row. Figure 10a shows the predicted transient evolution of the liquid and gas flow rates at the exit of the network, and the overall pressure difference, when the liquid flow

Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991 963 rate in the distribution manifold is raised so that ReL goes from 1029 to ReL = 1374. The gas flow delivered to the distribution manifold is constant (ReG = 115). The overall pressure difference and the liquid outflow rate increase steadily, as does the liquid saturation in the passages. The output gas flow rate first falls and then passes through a minimum as the liquid saturation in the passages grows and the area for gas flow shrinks. A fraction of the gas flow arriving at the distribution manifold, unable to flow down, accumulates, and the gas pressure rises. The buildup of pressure in the distribution manifold forces more gas to flow through the network; the gas outflow mounts until it reaches the steady-state value. In the steady state the sites are gas dominated, and the flow regime in the passages i s gas-continuous. Both phases share the void space. Figure 10b shows more of the transient evolution as the liquid load increases (Rec = 115 still). When the liquid rate to the distribution manifold is raised from ReL = 2400 to ReL = 2745, the liquid films in the passages become unstable. The flow regimes in some passages evolve from the gas-continuous regime to the bridged regime and back again in a cyclic sequence-the local pulsing regime; other passages evolve to the flooded regime. The sites evolve from gas dominated to contested to liquid dominated. Liquid-dominated clusters form and grow, and the predicted liquid outflow rate falls momentarily as liquid clusters form. As the liquid clusters appear, the number of flow passages available to the gas phase shrinks and the local gas flow through these passages is heavier. Heavier local gas flows produce drainage in neighboring liquid clusters. As the liquid clusters change in size and position within the network, fluctuations in local liquid and gas flow rates are noticeable until the next steady state is reached. When the liquid flow to the distribution manifold is reduced, the gas displaces the liquid from some sites at the margins of liquid-dominated clusters. When the gas phase invades the sites, adjacent liquid-filled passages become accessible to gas, liquid bridges in them become unstable, and the gas phase invades those passages. Liquid clusters move down and fluctuations in both the gas and liquid flow rates at the bottom of the network are noticeable as the liquid clusters shrink, descend, and exit. Finally the steady state is reached at a lower pressure difference and the fluctuations in gas and liquid flow rates cease. Some passages, though allowable to gas-continuous regime, continue to perform in bridged or flooded regimes because the requisite accessibility criteria are not satisfied. The phenomenon of hysteresis is illustrated by the snapshots of predicted steady states shown in Figure 11. The effect of the highest liquid load attained in a scan of liquid rate is shown in Figure 12. The higher the liquid load, the more and the larger are the liquid clusters that form. These clusters shrink as the liquid flow rate is subsequently reduced, but the pressure difference then is higher than when the liquid flow rate is raised to the same value. The effect of liquid distribution at the top of the network is shown in Figure 13. When the liquid is fed to the two sites at the center of the first row, the gas flows through the sides of the network and the liquid flows through the center in parallel paths; in other words channeling, a phenomenon much noted in the literature, is also noticeable here. The overall pressure drops predicted when the liquid feed is concentrated in the center of the network are all lower than the pressure drop predicted with uniform liquid distribution. Pulsing Regime. The network employed in experiments by Melli et al. (1990) had 24 rows of sites with 10 sites in each row. In this network the pulse frequency

/I n

"ii, ,

I i

L;,=

2250

0.25

0

0 FREQUENCY (HERTZ)

11

B

s4 vi

3

I

OTlME (DIMENSIONLESS)

1000

1

0.751

l

0.25

q1

0

L

R e c = 1875

0 FREQUENCY (HERTZ)

It

2s -

300

I

L,- J

OTlME (DIMENSIONLESS)

300

IOW

I

0.i5

,

~,

_,

ROG,

I

1312

025

0

0 FREQUENCY (HERTZ) '1 I

300

e,- 1

OTIME (DIMENSIONLESSJ

LOW

ReG = 900 035

O U C h 0 FREQUENCY (HERTZJ

,

Joo

Figure 15. Evolution of pulse frequency with gas flow rate in a network of 3 by 10 sites. Constant liquid flow rate, ReL = 1583.

ranged between 40 and 70 Hz and the length of pulses ranged between four and eight rows, depending on the gas and liquid rates. Lower frequencies and pulse sizes were found near the transition between the pulsing and gascontinuous regimes-at low gas and liquid rates. The highest frequencies and pulse sizes were found in the transition between pulsing and bubbling regimes-at high gas and liquid loads. However, instances of individual passages changing regime with much higher frequency were common. That is, at the microscale there were time scales much shorter than at the macroscale. The total number of variables needed to model twophase flow through this network is 2229. The computational effort required to predict completely from the theory the pulsing regime in this network of passages can be estimated from the experimental results together with the minimum record length required for frequency analysis given by (27). Visualization of flows in networks of passages, accomplished with the aid of a high-speed video camera, revealed that a characteristic pulse travels along a test section nine rows long in 18 ms (Melli et al. 1990). Acoustic sensing of pressure fluctuations at the bottom of the network revealed that the range of characteristic pulse frequencies spans 30 Hz (Kolb et al., 1990) in the same network. If the minimum desired frequency resolution, F, is of the order of 1 Hz, the minimum record length sampled, from (27), is t, = 1 s. Thus, in 1 s 55 pulses exit the test section. Within slugs of enhanced gas holdup, the passages perform in local trickling flow and the sites are dominated by gas. Within the slugs of enhanced liquid holdup, the

964 Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991

ja)Center: First pulse travels ' down the network, t = 0 msec

(1, tiot t om : First p ti I be reactie\ the nctwork exit, t =7.2 msc'c

(b)lop: Second pulse forms t = 2.4 msec

(c)Top: Second pulse grows t = 4.8 msec

'ei('enter: Second pulse grows and travels down, t = 9.6 niser

Figure 16. Overtaking mechanism predicted by statistical network theory in a network of 3 by 10 sites. Gas Reynolds number Rec = 2000; liquid Reynolds number ReL = 1424.

passages perform in local bubbling and the sites are dominated by liquid. In the front and back of each slug are contested sites. Ahead of slugs of enhanced liquid holdup, passages evolve from trickling to pulsing to bubbling, and conversely behind such slugs. To predict accurately saturation and flow rates in the passages and pressures and accumulation in the sites with the statistical network theory, the liquid and gas invasion processes should be so discretized in time that the theoretical predictions are independent of the time step selected. The time step is calculated by (24) and is set by the maximum change allowed for the jth variable, Aymax. When Ay" for the saturation in the passage is fixed at 5 % or lower and Ay- for the liquid accumulation in the sites is fixed at 10% or lower, the theoretical predictions turn out to be time-step independent. With these maximum allowable changes in the variables, the minimum number of time steps for bonds performing in gas-continuous and bubbling regimes is 10 but it is 40 for those performing in the local pulsing regime. Ten time steps are also required to fill and empty the sites. Thus the overall cyclic evolution in a site and related bonds requires a minimum of 140 time steps, a characteristic of the microscale. Therefore, for a network of 24 rows the total number of time steps needed to resolve a single macroscale pulse is 3360-provided the changes in sites and passages located in the same row occur simultaneously. These 3360 time

steps should be repeated at least 55 times to complete 1 s of operation (the minimum record length required for frequency analysis). Now each time step requires at least 0.5 CPU s in a Cray I1 to compute the solution. Therefore with the present form of the theory, 25 or more CPU hours would be required to evaluate the pulsing regime at a single set of gas and liquid flow rates and a single decoration of site and bond properties in a network the size of the experimental one of Melli et al. (1990),which is itself a tiny model of an industrial packed bed. The available computational resources, even today's advanced supercomputers, fall short. However, important conclusions can be drawn by studying the pulsing regime in networks of passages smaller than 24 X 12 sites. Pulsing Regime in a Long Narrow Network. In this example the evolution of the pulse frequency predicted by the theory is studied in a network of 10 rows with 3 sites in each row. The passage throat radius rt = 1 mm, the passage inclination is 0 = 30°, the ratio of body radius to throat radius is rb/rt = 1.2, and the ratio of passage length to throat radius is Zo/rt = 4; these are the same in all passages. The average site volume is 15 times the average passage volume. Figure 14 shows the predicted pulsing regime at intermediate gas and liquid flow rates, corresponding to ReG = 2000 and ReL = 1424. In the pulsing regime, waves or slugs of enhanced liquid holdup and enhanced gas holdup chase each other down the network. The size of the waves

Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991 965 Constant Cos Reynolds Number, Reo

Re~z2374

5

3000

I: 4

tu

A

0

FREQUENCY(HERTZ)

300

OTIME (DIMENSIONLESS)

p

~

IWtO

Figure 19. Converging-diverging passage and its characteristic dimensions when it performs in bridged regime at maximum gas pressure difference that a liquid bridge can support.

0.75

e

2 2

025

n 1)

FREQUENCY (HERTZ) 3(nl

Figure 17. Evolution of pulse frequency with liquid flow rate in a network of 3 by 10 sites at constant gas flow rate, Rec = 3000.

Figure 18. Converging-diverging passage and its characteristic dimensions a t maximum gas pressure difference that a liquid drop can support.

is virtually constant and spans two or three rows of sites in this example. Further insight into the pulsing regime comes from analysis of the evolution of pulse frequency a t constant gas rate and varying liquid flow rate, and the evolution of pulse frequency at constant liquid flow rate and varying gas rate. The pressure difference oscillates between a maximum value that is attained just before the gas-rich slug reaches the bottom of the network and a

minimum value attained as the liquid-rich slug starts to form. The gas rate peaks as the gas breaks through the network and passes through a very low minimum when the gas travels as dispersed bubbles as the liquid slug reaches the bottom of the network. The liquid flow rate changes sharply when the lowest row of bonds changes from gascontinuous to bubbling and vice versa. At a constant liquid flow rate corresponding to ReL = 1583, the pulsation frequency climbs as the gas flow rate is raised. When pulses of different sizes flow down the network, multiple pulse frequencies are detected, as recorded in Figure 15. Pulses of different size travel down the network at different velocities. Melli et al. (1990) found in experiments with a network 12 sites wide by 24 sites long that the larger a liquid pulse the faster it travels as it descends, and so it can overtake and merge with the pulse ahead to produce a yet bigger and faster pulse; the pulse frequency is not constant when this is happening in the network. Figure 16 shows this overtaking mechanism, as predicted by the theory in a network 3 sites wide by 10 sites long. A t a constant gas rate corresponding to ReG = 3000, multiple sizes of pulses and multiple frequencies appear at ReLbetween 1424 and 1741. Pulse size mounts as liquid flow rate is raised to the value at which the pulse is large enough to span the network length; beyond this value the pulse frequency decreases as the liquid flow rate is raised; see, for example, Figure 17 at ReL = 2374 and 3166. Examples of the theoretical predictions summarized in this section are recorded in a computer-generated videotape; comparable experimental observations of a model network of passages are recorded in another videotape (cf. Melli (1989)).

Conclusions The hydrodynamics of two-phase flow in a network of passages, or bonds, and enlarged junctions, or sites, is the statistical combination of four local flow regimes in the individual passages. The flow regimes in the passages are outcomes of the local competition of the two phases for

966 Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991

the void space. The constitutive equations of the local flow regimes, the conservation equations in sites and bonds, and the allowability and accessibility rules governing changes in local regimes together provide a framework of theory with which to predict the time evolution of the macroscale flow regimes most often important in trickle-bed reactors: the trickling regime, the pulsing regime, and their transitions. In the trickling regime, the way the liquid is introduced at the top of the network can influence strongly the distribution of gas and liquid in the void space, e.g., the occurrence of channeling. Moreover, the history of the gas and liquid flows sways that distribution. So it is that more than one steady state exists at a given gas and given liquid rate; in other words, there is hysteresis in the trickling regime. When there are multiple steady states, they depend on the history of the operation. Hysteresis is the outcome of different local flow regimes performing in the individual passages. In the pulsing regime, the local flow regimes predicted by the theory in the liquid-rich and gas-rich slugs are the same as those found experimentally in networks of passages by Melli et al. (1990) and Kolb et al. (1990). The theory predicts accurately the qualitative behavior of two-phase flow in trickle beds operating in steady, transient, and oscillatory states. It makes plain the roles that gas and liquid flow rates, the necessarily constricted passages, the fluid properties, and the history of the operation have in the macroscale flow regimes. Extensions and improvements are needed before the theory can be used for the quantitative prediction of flow maps, pressure drop, and holdup in trickle beds. One obvious possibility is to formulate constitutive equations that average over local pulsing in passages, in order to reduce the number of time steps required to follow macroscale pulsing. The sensitivity of results to changes in the constitutive equations has to be studied: for one thing, the constitutive equations were derived from experiments using air and water at atmospheric pressure. The theory must be extended to predict the behavior of larger three-dimensional networks and to include the effect of less-than-complete wetting of passage walls by the liquid. Nevertheless, the conclusion from the predictions of the version of the theory reported in this paper is inescapable: the macroscale flow regimes are rooted in the microscale flow mechanisms, that is, in the mechanisms by which gas and liquid compete for individual passages and sites-the "void space". Acknowledgment The research reported here and in the companion papers by Melli et al. and Kolb et al. was partially supported by grants from the Chevron Research Corporation and the Minnesota Supercomputer Institute. We are indebted to J. M. de Santos for help in bringing this paper to publication. Appendix. Derivation of Equations and Limits of Flow Regimes in Individual Passages Two-phase flow in constricted passages is governed by the Navier-Stokes system of equations. This system is made more complicated than usual by the convergingdiverging shape of the passage and the free surfaces present. The solution of the governing equations simply cannot be expressed in terms of tabulated functions, i.e., put in "analytic form", unless approximations are made.

Axial symmetry is adequate for the purpose at hand. The constitutive equations developed must be simple enough relations between the fluxes and potentials to be readily useful in building a statistical theory to relate macroscopic behavior to what goes on in individual passages. Gas-Continuous Regime. For steady flow of incompressible fluids in steady state, the Navier-Stokes system reduces to P ~ V ~ * -V psg V ~ - V*T, = 0

(A-1)

v.v, = 0

(-4-2)

Here vs is the velocity in phase p, g is the gravitational body force, and T is the total stress tensor. The four boundary conditions are symmetry at the centerline, no slip at the wall, and stress and velocity continuity at the gas-liquid interface. Solving the system (A-1) and (A-2) for a constant radius passage, or a converging-diverging passage with average radius, and integrating the velocity distribution in the gas from the centerline to the gas-liquid interface and that in the liquid from the gas-liquid interface to the wall yields the constitutive equations for gas and liquid flow. The constitutive equations in terms of the units of flow rate and pressure listed in the Governing Equations section are

Lout= k L ( s 2 ( u- COS e) + 2(1 - s)[s + (1 - s) In (1- s)] cos BJ(rAv/ro)4(A-3) G = kL{2(1 - s ) cos e[(l - s) In (1 - s ) + s] + [ ( l sP/m + 2 0 - s ) I U l ( r ~ v / r o )(A-4) ~ Here s is the liquid saturation, m is the gas-to-liquid viscosity ratio, and kL = 1/16 is the passage conductivity. The key to using this approximation is the choice of length and radius, rc, that predicts accurately the pressure drop and liquid saturation in the trickling regime for any given constricted tube. To make the length of the cylinder equal to the length of the constricted tube, I , is a natural choice. The best agreement of pressure drop with theoretical predictions of de Santos et al. (1989) was reached when r, = rAvand s = st (rt/rAVY,where st is the saturation in the throat. This relation represents the liquid inventory; s is an average saturation in the passage. These expressions were used in (4) and ( 5 ) in Table I. In the limiting case of s = 1and G = 0, (A-3) reduces to the equation that describes the flooded regime, namely, (8) in Table I. Bubbling Regime. No simple constitutive equations of the bubbling regime in small constricted passages are known to the authors. However, the bubbling regime in large cylindrical tubes, typical of nuclear technology, has been studied extensively. Vasilev and Maiorov (1979) idealized the two-phase gas and liquid mixture as a single fluid with average properties. An equation often used relates the potential difference to the sum of a linear and a quadratic function of the two-phase velocity: A P / l = apij + pp02 (A-5) Another equation, the drift velocity relation (Ishii, 1977), relates gas and liquid flow rates: UG = vGmS3 + ij (3-6) Here p pG + SAP is the density of the mixture; ij e (L+ G ) / A is the volumetric average velocity of the mixture and A rrAV2; UG 1 G / A (1 - s) is the volumetric average velocity of the gas; uGm is_the drift velocity of the bubbles relative to the mixture; P'= P + pgl is the potential; and (Y and 8 are the viscous and inertial resistance coefficients.

Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991 967 In terms of the characteristic units defined earlier, (A-5) is (see (10) in Table I) .

-AF

Stokes' law, derived for a bubble rising in an infinite pool of still liquid, gives the functional form of the bubble velocity relative to the mixture, i.e., the drift velocity: UGm

(A-8)

= Apgrb2/9h

To account roughly for the effect of the passage size on bubble drift velocity, (A-8) can be modified by multiplying the right side by y = (rb/rAv)2, the ratio between the bubble and passage radii. The presence of other bubbles, especially when a bubble moves in the wake of another bubble, influences the relative velocity of a bubble. From experiments of Ishii and Chawla (1979),the effect of higher bubble concentration on drift velocity is proportional to the cube of the liquid saturation, s3. The expression relating gas flow rate and liquid flow rate is obtained by substituting the drift velocity from (A-8) and the definitions of the average velocity of gas and of the mixture in (A-6). In terms of the units in the main text (see (9) in Table I) 7 G = --s2(1 9

- S) + L-

1-s S

(A-9)

Bridged -me. In the bridged regime the liquid clogs the passage's throat but does not fill the entire space: gas is also present within the passage. The bridged regime is the link between trickling and flooded regimes, and so it is appropriate to choose the functional form of its constitutive equation similar to the form of those of the trickling and flooded regimes: (A-10) L = kj(S)(AFL + P L d

average saturation, s, in elongated constricted tubes is about 0.7. The third limit appears when no liquid flows through the passage, Le., L = 0, but a liquid drop clogs the throat. This saturation, smin, corresponds to the largest excess of upstream gas pressure over downstream gas pressure that a liquid bridge can support, viz., APG = AP- In this m e the upper and lower menisci kiss (see Figure 18). The liquid saturation smh is the ratio between the volume of liquid in the passage, VL, and the volume of the passage, Vp. The volume of liquid in the passage is calculated as the volume of the diverging cone R3 minus the volumes of gas above and below the throat, viz., R1and R2. The volume of the passage is the volume of two fused cones of junction radius rt, base radius rb, and height 1

Vp = ( ~ / 3 ) r 2 ( 1+ 5

+ f2)u

(A-13)

where 5 = rb/rt is the ratio of body radius to throat radius and u = l / r , is the ratio of passage length to throat radius. The volume of gas in region R1is that of the hemisphere of radius rc Vl = ( 2 / 3 ) ~ r ? (A-14) The volume of gas in region R2 is that of the spherical cap of radius rL and height rL (1- sin 8): (A-15) V2 = (.rr/3)rL3(3- sin 6) The volume of the gas in the region R3 is the volume of the cone frustum of height 112 - r, - rL sin 6, base radius rb, and top radius rL cos 8:

v, = (A-16) The lower meniscus radius rL, the throat radius r,, and angle 8 can all be written in terms of 5 = rb/rt and u = Z/rt. From (A-12) to (A-15) and simple geometry, the minimum saturation in terms of 5 and u is

At each of the gas-liquid interfaces, the gas and liquid pressures differ by the capillary pressure; accordingly, the relation between APL and AF', is AFL = b p ~ 2[%u(S) + %L(S)] (A-11) Spherical caps give the simplest approximation to the mean curvatures of the upper and lower menisci of the bridge, namely, Su(s) = l / r u and 7tL(s)= l / r L (as in Figure 19, which shows only the limiting case in which the menisci kiss, however); both are functions of the liquid saturation. From (A-10) and (A-11) the form of the constitutive equation of the bridged regime is (see (7) in Table 1) (A-12) QL = k L [ f l ( S ) ( h P + 1) + fib)] Here f z = 2uf(s)[%&) + PL(s)]and f i b )= f ( s ) depend on f ( s ) and the positions of the lower and upper menisci and are approximated by quadratic equations because three limiting cases of the bridged regime are known. The first limit of the bridged regime appears at full saturation, i.e., s = 1. This is the flooded regime. Equation A-12 should match the liquid rate predicted by (A-5). The second limit appears at a saturation, sm, at which the liquid film in trickling regime collapses and clogs the throat. At this saturation the liquid rate in the trickling and bridged regimes should match. From experiments of de Santos et al. (1989) and theoretical predictions about the trickling regime by de Santos (1989) the maximum

With the three limits, the functions fib) and fib) are readily evaluated. Flow Regime Boundaries: Allowability Criteria. Experimental results were used to locate the boundaries between the four basic flow regimes. Experiments with elongated and squat tubes made by de Santos et al. (1989) showed that the gas-continuous regime at low gas flow rates is stable at throat saturations up to about 0.9 (the average volumetric saturation in passages with rblrt = 1.2, l / r , = 13, and throat saturation t = 0.9 is roughly, ,s = 0.7). The transition from the gas-continuous regime to the bubbling and bridged regimes occurs around this value of throat saturation, s = sm,. At saturation s L s, and liquid flow rate such that bubbles can be carried by the liquid flow (see eq A-9), the gas-continuous regime evolves into bubbling regime; otherwise the passage performs in the bridged regime. The limit between the flooded regime and the bubbling regime is also given by the liquid flow able to drag bubbles. In

968 Ind. Eng. Chem. Res., Vol. 30, No. 5,1991

the bridged and flooded regimes, there is no gas flow at all, i.e., G = 0. The bridged regime gives way to gas-continuous regime when the pressure difference exerted by the gas phase is enough to break the liquid bridge. The maximum pressure difference that the two menisci of a liquid bridge in the throat of a passage can support is hp = 2[%u(s) SL(s)]= 2/rU + 2/rL (A-18)

+

At any position x j in the passage the local liquid saturation, ti, is related to the position of the wall, rj, and the position of the interface, 3, through

ti = 1 -

(rj/rj)2

(A-19)

Therefore, the upper meniscus radius is related to the saturation t at the throat, and to the throat radius, rt, by (cf. Figure 19) ru = r,(l - t ) 1 / 2

(A-20)

When the cone angle 6 is low, the saturation in any position j is related to the saturation in the throat through the ratio of cross-sectional areas: tj =

t(rt/rj),

Here u = 1f rt is the ratio of the length of the passage to the throat radius, and 7 rL/rt is the ratio between the lower meniscus radius and the throat radius. The lower meniscus radii, rL, is the radius of the circle that osculates the upper meniscus and is tangent to the film below the throat at a distance xL from the throat and yL from the axis of the tube: xL = r,(l - t ) 1 / 2+ rL(l - sin 6)

(A-23)

yL = rL cos 6 = rL

(A-24)

Simple geometry gives the relation between the position of the wall, rj, the throat radius, r,, the cone angle, 6, and the distance from the throat, x,: rj = r, + x j tan 6 (A-25) From eqs A-19-A-21 and A-23-A-25, the lower meniscus radius satisfies rL2cos2 6 = (r, + xL tan - tr,? (A-26) Finally, eq A-27, which relates q with the cone angle and the liquid saturation in the throat, is obtained when eq A-23 is replaced in eq A-26 and the result is divided by the square of the throat radius: q2[cos26 - (1 - sin @,tan2 61 1[2(1 - sin 6 ) tan 6 ( l + (1 - t ) 1 / 2tan e)] [ ( l - t ) ( l + tan2 6) + 2 0 - t V 2 tan 61 = 0 (A-27) When t is known-through its relation with the volumetric average saturation given below--? is the positive root of eq A-27. This relation is calculated by a procedure similar to that used to estimate the minimum saturation smin.The saturation, s, is (see Figure 19) vp

-

v,- vi - v,- v, VP

(A-28)

(A-29)

Here the integral is readily evaluated, by substituting for t j ( x ) , the local saturation at the position j , with eq A-21, which relates ti to the saturation at the throat and the throat radius. The volume of the pore, Vp, and the volume of gas in the regions R1, R2, and R3, are those derived in the calculation of smin.The radius of the upper meniscus, ru, is related to the throat saturation through (A-20). When all equations are substituted in (A-28),the volumetric average saturation in the passage in terms of 17 = rLf r,, 5 = rbf r,, u = 1f r,, 6, and t is found to be

[ ( l - t)'/, + ~ (- lsin 6)][l

+~

( -1sin

-

+ 2(1 - t)'/, + q(1 - sin 6 )

(1 - t ) 3 / 2- v3(l - sin 6)(2 + sin 6 ) / 3 ] (A-30)

(A-21)

The maximum pressure difference in terms of the unit dimensions defined in the text (see Table 11) is

s=

The volume of gas in the region R,, is 1/2 1 V, = sVp rr?(x) t j ( x ) dx

Concluding Remarks Equations A-22, A-21, and A-30 are the criteria used to calculate the largest excess of upstream gas pressure over downstream pressure that a liquid bridge can support, P,,. The shape of the passage is described in terms of its throat radius, r,, alone; 5 = rbf r, and u = l / r , are considered to be the same in all the passages. Given the throat radius, for example for rt = 1 mm, the maximum pressure difference hp-, is a function only of the passage saturation, s (see eq A-22). The saturation s is readily evaluated for different values of the throat saturation, t , by solving first (A-27) for 7 and then computing s from (A-30). In this way hpJrt=l m m is tabulated versus s. For tubes with rt # 1 mm, the maximum pressure difference is immediately evaluated as U I r t = S I r , = l mm/rt2 (A-31) Literature Cited Chu, C. F.; Ng, K. M. Model for Pressure drop hysteresis in trickle-beds. AIChE J . 1989,35, 1365-1369. Christensen, G.; McGovern, S. J.; Sundaresan, S. Cocurrent Downflow of Air and Water in a Two-Dimensional Packed Column. AIChE J . 1986,32, 1677-1682. de Santos, J. M.; Melli, T. R.; Scriven, L. E. Cocurrent Downflow in Packed Beds: Basic Flow Mechanisms in Individual Passages. Presented at the AIChE Spring National Meeting, Houston, TX, Apr 2-6, 1989. Ergun, S. Fluid Flow through Packed Beds. Chem. Eng. Prog. 1952, 48,93-99. Heiba, A. A.; Davis, H. T.; Scriven, L. E. Statistical Network Theory of Three-phase Relative Permeabilities. SPE Paper 12690. 4th SPEIDOE Forum Symposium, Tulsa, OK, Apr 15-18,1984; Society of Petroleum Engineers: Dallas, TX, 1984. Heiba, A. A.; Jerauld, G. R.; Davis, H. T.; Scriven, L. E. Mechanism-Based Simulation of Oil Recovery Processes. SPE Paper 15593.61st Annual Technical Conference SPE, New Orleans, LA, Oct 5-41986; Society of Petroleum Engineers: Dallas, TX, 1986. Ishii, M. Drift Flux Model and Derivation of Kinematic Constitutive Laws. In Proceedings of the 1976 Meeting on Two-Phase Flows and Heat Transfer; Mayinger, F., Ed.; NATO Adv. Study Inst.; Hemisphere Publishing: Washington, DC, 1977; Vol. 1, pp 187-206. Ishii, M.; Chawla, T. C. 'Local Drag Laws in Dispersed Two-Phase Flow". Report ANL-79-105; Division of Reactor Safety Research, Nuclear Regulatory Commission: Washington, DC, 1979. Ishii, M.; Zuber, N. Drag Coefficient and Relative Velocity in Bubbly, Droplet or Particulate Flows. AIChE J. 1979, 25, 843-849. Kan, K. N.; Greenfield, P. F. Pressure Drop and Holdup in Two-

Ind. Eng. Chem. Res. 1991, 30, 969-978 Phase Cocurrent Trickle Flows through Beds of Small Packings. Ind. Eng. Chem. Process Des. Deu. 1978,18,740-746. Kolb, W. B.; Melli, T. R.; de Santos, J. M.; Scriven, L. E. Cocurrent Downflow in Packed Beds. Flow Regimes and Their Acoustic Signatures. Ind. Eng. Chem. Res. 1990,29,2380-2389. Larkins, R. P.; White, R. R.; Jeffrey, D. Two-Phase Cocurrent Flow in Packed Beds. AIChE J. 1961,7,231-239. Levec, J.; Saez, A. E.; Carbonell, R. G. Holdup and Pressure Drop in Trickle Bed Reactors. Inst. Chem. Eng. Symp. Ser. 1986,87, 185-194. Melli, T. R. Two-Phase Cocurrent Downflow in Packed-Beds; Macroscale from Microscale. Ph.D. Thesis, University of Minnesota, Minneapolis, MN, 1989. Melli, T. R.; de Santos, J. M.; Kolb, W. B.; Scriven, L. E. Cocurrent Downflow in Networks of Passages. Microscale Roots of Macroscale Flow Regimes. Ind. Eng. Chem. Res. 1990,29,2367-2379. Rao, V. G.; Ananth, M. S.; Varma, Y. B. G. Hydrodynamics of Two-Phase Cocurrent Downflow through Packed Beds. Part I. Macroscopic Model. AIChE J. 1983a,29,467-472. Rao, V. G.; Ananth, M. S.; Varma, Y. B. G. Hydrodynamics of Two-Phase Cocurrent Downflow through Packed Beds. Part 11: Experiments and Correlations. AIChE J. 1983b,29, 473-477. Saez, A. E.; Carbonell, R. G. Hydrodynamic Parameters for GasLiquid Cocurrent Flow in Packed Beds. AZChE J. 1985, 31, 52-62. Saez, A. E.; Carbonell, R. G.; Levec, J. The Hydrodynamic of Trickling Flow in Packed Beds. Part I Conduit Models. AIChE J. 1986,32, 353-368.

969

Sherwood, T. K.; Shipley, G. H.; Holloway, G. A. L. Flooding Velocities in Packed Columns. Ind. Eng. Chem. 1938,30,765-779. Stanley, W. D.; Dougherty, G. R.; Dougherty, R. Digital Signal Processing; Reston: 1984,Reston, VA. Stauffer, D. Introduction to Percolation Theory; Taylor and Francis: London, England, 1985. Talmor, E. Two-Phase Downflow through Catalyst Beds. Part I; Flow Maps. AIChE J. 1977,23,868-873. Tosun, G. A Study of Cocurrent Downflow of Nonfoaming GasLiquid Systems in a Packed Bed. 1. Flow Regimes: Search for a Generalized Flow Map. Jnd. Eng. Chem. Process Des Dev. 1984, 23,29-34. Vasiliev, L. L.; Maiorov, V. A. An Analytical Study of Resistance, Heat Transfer and Stability in Evaporative Cooling of a Porous Heat-Producing Element. Int. J. Heat MQSSTransfer 1979,22, 301-306. Westerterp, W. P. M.; Van Swaaij, W. P. M.; Beenackers, A. A. C. M. Chemical Reactor Design and Operation; Wiley: New York, NY, 1984. White, A. M. Pressure Drop and Loading Velocities in Packed Towers. Trans. AIChE 1935,391-397. Zimmerman, S.P.; Chu, C. F.; Ng, K. M. Axial and Radial Dispersion in Trickle-Bed Reactors with Trickling Gas-Liquid Downflow. Chem. Eng. Commun. 1987,50,213-240. Received for review December 27, 1989 Revised manuscript receiued December 4, 1990 Accepted December 14, 1990

Design of Monolith Catalysts for Power Plant NO, Emission Control Jean W. Beeckman and L. Louis Hegedus* W.R. Grace & Go.-Conn., 7379 Route 32, Columbia, Maryland 21044

Research Diuision,

A mathematical model was developed to describe the reduction of NO with NH, in simulated power plant stack gases, over the internal surfaces of monolith-shaped catalysts of the vanadia-titania type. The model indicated that up to 50% improvement in the volumetric activity of these catalysts would be possible if they could be prepared with the computed optimum pore structure. Based on the guidance of the computer calculations, a new type of catalyst was developed. It consists of vanadia, titania, and silica and shows the predicted 50% improvement in activity in laboratory experiments.

Introduction The selective catalytic reduction (SCR) of nitrogen oxides in power plant stack gases is in large-scale commercial practice in Japan and West Germany. Tightening nitrogen oxide emission standards are expected to stimulate the application of SCR processes in the United States as well (Boer et al., 1990). The overall chemical reaction for NO reduction in the SCR process can be expressed as follows: 4N0 + 4NH3 + O2 4N2 + 6H20 The above reaction is equimolar with respect to NO and NH,, and control of NO removal in commercial units can easily be accomplished by injecting the required amount of NHD The SCR process employs either ceramic monoliths or catalyst-coated metal plates as catalytic elements. Typically, 1-1.5 m3 of catalyst is required for each megawatt of electric power generating capacity. A typical coal-fired power plant has a capacity of about 700 electric MW, and the United States has a total fossil-fuel-firedelectric power generating capacity of about 530 000 MW. On the basis of recent German experience, the SCR process increases the cost of electricity by about 5%,while providing about 80% NO, removal (Gouker et al., 1989). Approximately 50% of the "levelized" cost (i.e., the sum of operating charges and prorated capital charges) of SCR is catalyst-related (Boer et al., 1990). Therefore, a powerful

-

economic incentive exists to improve the performance of present-day SCR catalysts. Briefly, the process employs vanadia or tungsten-vanadia as active components deposited on porous anatase-type titania which is extruded into the shape of square-channeled monoliths (Ando, 1983). Figure 1shows one such element. The catalyst can be placed either before or after the fly ash precipitators. Typical operating temperatures are in the range of 300-400 "C, and gas hourly space velocities (GHSV) are in the 2000-7000 h-'(STP) range. A typical catalyst block is 70-100 cm long with a face area of 15 cm by 15 cm. The channels are 0.3-1.0 cm wide, and the porous catalyst walls are between 0.05 and 0.15 cm thick. Several outstanding reviews of this technology are available (Ando, 1983, 1985; Bosch and Janssen, 1987). The current SCR process is the result of about 15 years of development, primarily in Japan. The objective of the present authors was to investigate what further potential may exist for catalyst improvements and to identify avenues that would lead to the indicated improvements, if any. For this, we have employed the methods of chemical reaction engineering. Specifically, we have set out to determine the intrinsic kinetics of the key SCR reactions over typical catalytic surfaces, to determine the pore structure of typical SCR catalysts, to determine the effective diffusivity of the key reactant and product species in the porous walls of SCR

0888-5885/91/2630-0969$02.50/00 1991 American Chemical Society