A numerical study of the interactions between viscous flow, transport and kinetics in fixed bed reactors

Computers and Chemical Engineering 26 (2002) 333– 357 www.elsevier.com/locate/compchemeng

A numerical study of the interactions between viscous flow, transport and kinetics in fixed bed reactors Hugo A. Jakobsen *, Ha˚vard Lindborg, Vemund Handeland Department of Chemical Engineering, The Norwegian Uni6ersity of Science and Technology, NTNU, Sem Sælands 6ei 4, NO-7491 Trondheim, Norway Received 29 January 2001; received in revised form 13 August 2001; accepted 13 August 2001

Abstract A numerical method is tailored for the solution of a reduced set of model equations developed for the description of reactive flows in chemical reactors. The numerical method synthesis intends to utilize experience both in chemical engineering solving dispersion models containing complex reaction kinetics, phase equilibria and non-ideal thermodynamics, and in fluid dynamics applying classical numerical algorithms constructed for pure flow calculations. The two types of model equations involved in reactive flow simulations reflect very different physics. A modular method is therefore suggested enabling a split between the flow- and chemistry model parts. In this way more optimal solution methods can be adopted for each operator, in contrast to more traditional computational fluid dynamics (CFD) methods where all equations (and operators) are normally solved by the numerical solvers originally intended for pure flow calculations. In addition, emphasis has been put on modularity enabling the same model framework to be used both for the more traditional 1D and 2D dispersion models with simple 1D flow calculations, as well as for the more elaborate 2D and 3D reactive flow calculations. In the flow part of the algorithm a fractional-step algorithm is constructed for the simulation of dynamic reactive flows in chemical reactors. The numerical scheme is based on the compressible transport equations in the low-Mach-number limit. The method can handle real gas (and liquid) mixtures with variable density as well as constant density fluids. The velocity field is advanced using a projection scheme, which consists of a partial convection– diffusion update followed by a pressure correction step with an intermediate an-elastic filter. A variable density corrector step is implemented for variable density systems in order to couple the evolution of the density and the velocity fields. In the chemistry part of the model all scalar fields are updated using Strang-type operator-split integration steps that combine several explicit convection and semi-implicit diffusion transport operators with a suitable solver tailored for non-linear sink/source terms. Diffusion terms are discretized with a second order central difference scheme in space. Convection terms are discretized with a second order TVD scheme in space. The temperature advection is discretized using a second order upwind scheme. The chemical reaction part of the model is discretized by an implicit Euler approximation, and the resulting set of algebraic equations is solved using a Broyden subroutine. The other source terms are solved using an explicit Euler approximation. The performance and behavior of the operator-split scheme are assessed based on simulations of two industrial chemical processes (i.e. the synthesis gas and methanol production processes) performed in multi-tube fixed bed reactors. These processes are important parts of the Statoil methanol plant at Tjeldbergodden in mid Norway. Both 1D and 2D pseudo-homogeneous dispersion models (i.e. heat and mass balances) with prescribed velocity, mixture density and total pressure profiles are used describing the reactor performance. The predicted profiles for both the methanol and synthesis gas processes were in accordance with results reported in the literature. An oscillatory radial void fraction distribution was then implemented in the 2D model. It was found that the non-uniform void fraction distribution had no significant effect on the temperature and mole fraction profiles. A 1D CFD simulation was performed to evaluate the effect of the variations in velocity, pressure and mixture density on the reactor performance. The changes in these variables significantly effect the composition and conversion in the reactor. A 2D CFD model was then developed for further studies of the multi-tube fixed bed reactor for the synthesis gas process, analyzing the influence of non-uniform void fraction distributions on the flow and chemical conversion. The non-uniform void fraction -distribution induces a significant reduction in axial pressure drop and a much higher fluid velocity close to the wall. However, * Corresponding author. E-mail address: [email protected] (H.A. Jakobsen). 0098-1354/02/$ - see front matter © 2002 Elsevier Science Ltd. All rights reserved. PII: S0098-1354(01)00758-X

334

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

the influence on the temperature and mole fraction profiles was hardly noticeable. Mass- and enthalpy budgets were implemented for the heat and mass quantities in the program. These budgets show that the enthalpy and both the mixture and component mass balances are fulfilled in the simulations. The numerical CFD algorithm developed in this paper has been found to be both computationally stable and accurate. Computational efficiency analysis shows that the Poisson solver requires about 75% of the total computational time. The computational time spend on the chemistry part of the model is about 15% of the total time. About 10% of the cost used on the chemistry part is spend on the chemistry solver. © 2002 Elsevier Science Ltd. All rights reserved. Keywords: Viscous flow; Computational fluid dynamics; Fixed bed reactors

Nomenclature Ar Az Cp Cr CS Cz Da Dea Der Dr Dz dp dti f fD g H h J Ma Peh =usdpzCp/uea Pem =usdp/De Pr = vCp/uea Q R r rc Re =z0L60/v0 Rep =z 7 dp/v Sc = vg/zgDe T Ta t Ui 7 6z 6r yc z

numerical approximation of the radial advection operator (K/s) numerical approximation of the axial advection operator (K/s) heat capacity for gas mixture (J/kg K) numerical approximation of the radial convection operator (kg/m3 s) speed of sound (m/s) numerical approximation of the axial convection operator (kg/m3 s) Damko¨ hler number (– ) effective diffusivity, axial direction (m2/s) effective diffusivity, radial direction (m2/s) numerical approximation of the radial conduction operator (K/s), or numerical approximation of the radial dispersion operator (kg/m3 s) numerical approximation of the axial conduction operator (K/s), or numerical approximation of the axial dispersion operator (kg/m3 s) particle diameter (m) inner diameter of tube (m) fricton factor of porous media (–) Darcy fricton factor (– ) body force/gravity force (N/m3) numerical approximation of the heat of reaction operator (K/s) enthalpy variable (J/m3) mass diffusion flux (kg/m2 s) Mach number (– ) Peclet number for energy (– ) the Peclet number for mass (– ) Prandtl number (– ) numerical approximation of the external heating/cooling operator (K/s) tube radius (m), or numerical approximation of the reaction operator (kg/m3 s) position in the reactor, radial axis direction (m) rate of formation for component c (mole/kg cat) Reynold number ( – ) particle Reynold number (– ) Schmidt number ( – ) temperature (K) temperature at the inside of the tube (K) time (s) overall heat transfer coefficient (J/m2 s K) fluid velocity (m/s) axial velocity component (m/s), or superficial velocity (m3/m2 s). radial velocity component (m/s) mole fraction of component c (–) position in reactor, axial axis direction (m)

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

Greek letters h i DHr, j m k uea uer …c z zc zcat zg

335

isothermal compressibility (per Pa) expansion coefficient (per K) heat of reaction for reaction j (J/mole) void fraction (– ) composition induced expansion coefficient (–) axial effective conductivity (J/m s K) radial effective conductivity (J/m s K) mass fraction of component c (–) density (kg/m3) density of component c (kg/m3) density of catalyst (kg/m3) mixture density of gas (kg/m3)

Mathematical operators D/Dt substantial derivative 9 the ‘del’ or ‘nabla’ operator

1. Introduction The classical approach creating numerical solution algorithms for fluid flow models was based both on firm mathematical understanding of the model equations involved and on detailed analysis of the physical problems to be simulated. A survey of the existing algorithms (e.g. Peyret & Taylor, 1983; Ferziger & Peric´ , 1996) shows that the optimal methods have been developed after introducing a few classical model simplifications, creating systems of equations that are much easier to solve than the complete set of equations. An essential finding is that no algorithm is optimal for all flow situations although several fairly general algorithms have been constructed (e.g. Ferziger & Peric´ , 1996). Panton (1996) stated that in fluid mechanics, the term incompressible flow is applied to any flow situation where changes in the density of a fluid particle are negligible. In gas dynamics (e.g. Thompson, 1972), the flow of compressible fluids is considered incompressible for very low Mach numbers (Ma “0). This is based on an underlying assumption that heat (thermal energy) will only be generated in the flow field by viscous dissipation. Therefore, when the Mach number becomes small, the dynamic processes do not change the temperature or pressure enough to cause any appreciable change in the thermodynamic state, and all thermodynamic properties are considered constants. For incompressible flows of one component fluids under isothermal conditions, having fluid transport properties (i.e. density and viscosity) that can be considered constant, the unknowns are the velocity and pressure only. These variables can then be determined from the condition that 9 · 7 =0 (i.e. divergence free flow) and the momentum equation in an appropriate

form. In these systems, the most direct way to solve for the pressure is to combine the two equations to obtain a Poisson equation for the pressure (e.g. Peyret & Taylor, 1983; Patankar, 1980; Kim & Moin, 1985; Versteeg & Malalasekera, 1995). In this context it is also important to recognize that the energy principles are such that the loss of mechanical energy by dissipation represents an irreversible rate of conversion of energy into heat. Another term that represents reversible rates of conversion of energy from one energy form to another is the product of pressure and the fractional rate of change of volume of a material element. This latter term represents the rate of release of energy per unit volume of the fluid by expansion. In the limit of low Mach number flows, the energy dissipation term is very small and can be neglected (e.g. Majda & Sethian, 1985). For incompressible flows, the second energy conversion term vanishes. Therefore, for low velocity incompressible flows the energy equation reduces to a scalar transport equation for the temperature and only the advection and heat conduction are important. In such flow situations the numerical flow algorithms that were originally constructed for the simpler isothermal flows, sketched above, have been adopted introducing a numerical split between the flow and energy calculations. For low temperature variations, it has been assumed that the fluid transport properties (e.g. viscosity, conductivity, diffusion, heat capacity) are still approximately constant and the energy equation has been solved alongside the flow part of the model (e.g. Anderson, Tannehill & Pletcher, 1984, p. 503). To compute compressible flows, on the other hand, it is necessary to solve not only the continuity and momentum equations but also a transport equation for the energy or heat (e.g. in terms of temperature) and an

336

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

equation of state (EOS). In compressible flows, viscous dissipation may be a significant heat source and conversion of internal energy to kinetic energy (and visa versa) is also important. All terms in the equations must then be retained. For compressible flows it is therefore natural to use the continuity equation to compute the density and to derive the temperature from the energy equation. This leaves the role of determining the pressure to the EOS. The roles of the various equations are thus somewhat different from the ones they play in incompressible flows. The nature of the pressure is also different. In incompressible flows of fluids having constant fluid properties the pressure variable occur only in the pressure gradient term in the momentum equation, and the absolute value of pressure is of minor importance. For compressible flows, however, the absolute value of the thermodynamic pressure occur in the energy equation and is thus of critical importance. Due to the importance of compressible flow in gas dynamics many methods of solving the equations of compressible flow have been developed (e.g. Godunov, 1959; MacCormack, 1969, 1976; Beam & Warming, 1978; Sod, 1978). However, most of these methods are specifically designed for compressible flows and become very inefficient when applied to constant density flows. As the incompressible flow limit is approached for inert fluids under isothermal conditions, the density is constant along any streamline (that is, Dz/Dt = − z(9 · 7)= 0), and the speed of sound approaches infinity (Steger & Kutler, 1977, Anderson et al., 1984, p. 503). According to Ferziger and Peric´ (1996), the Navier –Stokes equations become extremely stiff in the limit of weak compressibility since the time derivative in the continuity equation drops out. As a result, to solve the equations, very small time steps or implicit methods should be used. In addition, the compressible equations support sound waves, which have a definite speed, associated with them. As some information propagates at the flow velocity, the larger of the two velocities determines the allowable time step in an explicit method. This step size may be much smaller than the one a method designed for incompressible constant density flows might allow. Under such conditions, one typically makes an incompressible assumption to remove the sound speed time scale, but recall that in doing so, the solution of a Poisson equation is usually required at each time step to obtain the pressure (e.g. Steger & Kutler, 1977). To handle low compressible flows some of these methods originally developed for incompressible flows (pressure correction methods) have been extended to variable density flows of arbitrary Mach numbers (see Ferziger and Peric´ (1996) and references therein). Alternatively, once a fully implicit scheme is used, the stiffness problem that motivated the incompressible approximation no longer exists. The use

of implicit methods requiring iterative procedures are, however, not always an optimal choice as the convergence rate is often low for very large systems of equations. In many reactive flow situations convergence is even never reached, as proper estimates for initiating the iteration procedure cannot be provided. To understand the numerical method synthesis that will be presented in Section 3, it may be useful to summarize some studies in numerical analysis that have induced important advancements in meteorological dynamics. Marchuk (1974) was the first to suggest that time-split methods (additive-splitting schemes) could be used to integrate numerical models of the atmosphere. Skamarock and Klemp (1992) discussed the more recent time-split schemes used in weather forecast models. As for most analyzes of chemical reactors, sound waves are not important in atmospheric motions of meteorological interest. However, as mentioned earlier, their high propagation speed severely restricts the maximum allowable time step in an explicit numerical integration. To eliminate this restriction several attempts have been made to remove the acoustic modes from the governing equations. It has been found that these modes are removed when the pressure in the momentum equation is decoupled from density and temperature fluctuations arising through the EOS (e.g. Rehm & Baum, 1978). Another group of methods were developed intending to speed up the integration of the fully compressible set of equations. Two equation sets are thus widely used in high-resolution forecast models: an-elastic equations, from which the sound waves have been filtered, and elastic equations, which retain the sound waves. The principal advantage of the an-elastic equations is that a reasonably large time step is allowed, the time step being limited by the fastest propagating gravity waves. However, use of the system requires the solution of a multidimensional elliptic pressure equation with every time step. Ogura and Phillips (1962), for example, introduced an an-elastic approximation by replacing the continuity equation with 9 · (z07)= 0. Rehm and Baum (1978) derived another set of approximate equations which is applicable to non-adiabatic, non-dissipative, buoyant flows of a perfect gas. The approximate equations are characterized by a spatially uniform mean pressure appearing in both the energy equation and the EOS with the spatially non-uniform portion of the pressure only appearing in the momentum equation. Therefore, the pressure remains almost constant in space while significant density and temperature variations are allowed. This formulation thus accounts for the compressibility of the medium but ignores the effect of elastic waves. Thus, it naturally avoids potential sonic CFL like limitations. This low-Mach-number formulation may also be useful for many reactive flows (e.g. Majda & Sethian, 1985; Najm, Wychoff & Knio, 1998;

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

Knio, Habib & Wyckoff, 1999). Using the fully compressible elastic equations where the time step is limited by the sound waves, the explicit integrations are much more costly than integrations of the an-elastic equations. Several methods have been developed to speed up the integration of the full elastic set. Tapp and White (1976) employed an elastic model wherein the terms responsible for the acoustic modes are handled implicitly. Another approach to integrate the elastic equations involves isolating the terms responsible for the acoustic modes and integrating them with a smaller explicit time step (e.g. Klemp & Wilhelmson, 1978). The computational effort required integrating the existing an-elastic and elastic models is roughly equivalent. The numerical solution of chemically reacting flow models are also associated with pronounced difficulties due to the inherently large ranges of spatial and temporal scales involved. The numerical algorithms that can tackle fluid flow problems are often not suitable solving stiff non-linear chemistry schemes, which are coupled with the flow and transport processes. In this context the utility of operator-splitting techniques derives from the advantage of the sequential application of individual operators in the model equations. Thus, the integration procedure over each split time step can be optimized for the individual operators independently, and computational efficiency can be enhanced. For this reason, operator-splitting techniques are widely used in the geophysical sciences (e.g. meteorological- and oceanographic studies) not only for fluid dynamics, but also for the physical and chemical parts of these models. For example, in air pollution modeling studies (e.g. Berge & Jakobsen, 1998; Jonson, Bartnicki, Olendrzynski, Jakobsen & Berge, 1998; Olendrzynski, Jonson, Bartnicki, Jakobsen & Berge, 2000), the approach has been used to decouple reaction from diffusion and convection terms, and to decouple operators in different spatial dimensions. For complex chemistry schemes, stiffness limitations are typically overcome by adopting a stiff-integration scheme (e.g. BDF schemes, or Implicit and semi-implicit Runge– Kutta schemes). The algorithms have been found to be both computationally efficient, stable, and accurate. In this work a numerical method is tailored for the solution of a reduced set of model equations developed for the description of reactive flows in chemical reactors. The performance and behavior of the operatorsplit scheme are assessed based on simulations of two industrial chemical processes (i.e. the synthesis gas and methanol production processes) performed in multitube fixed bed reactors. These processes are important parts of the Statoil methanol plant at Tjeldbergodden in mid Norway. This paper presents the first step in our multiphase reactor-modeling programme applying CFD to reactor technology. The numerical algorithm presented here

337

solving dynamic 2D CFD models for viscous reactive flows within porous fixed bed reactors, will be further developed alongside the continuation and expansion of this programme. The main purpose is to investigate the more complex chemical systems like bubble columns, stirred tanks and fluidized bed reactors for the Norwegian industry. The turbulent multiphase flow patterns observed in these systems require the third space dimension to be introduced into the model, as well as operators considering multiple phases. Further effort will thus continue to analyze the possibilities to enable time dependent solutions for these systems within reasonable time limits.

2. Model formulation

2.1. General reacti6e flow modeling aspects Considering any flow system the fluid or mixture density is governed by the continuity equation and a thermodynamic EOS. For a general multicomponent fluid, a thermodynamic state equation can be formulated assuming density to be a function of P, T and …c (e.g. Bird, Stewart & Lightfoot, 1960; Rosner, 1986, Najm et al., 1998) 1 Dz Dp DT D…c =h −i + % kc z Dt Dt Dt Dt c

(1)

where the isothermal compressibility (e.g. Sandler, 1999) is h(p,T,…c )=

1 (z z (p T,…c

(2)

the bulk expansion coefficient (e.g. Panton, 1996; Sandler, 1999) is i(p,T,…c )= −

1 (z z (T P,…c

(3)

and the composition induced expansion coefficient (e.g. Bird et al., 1960, Rosner, 1986) is kc (p,T,…c )=

1 (z z (…c P,T,…s " c

(4)

For reactive flows the extended criterion for incompressible flow is that the Mach number is low, and that no significant density changes occur due to changes in temperature and chemical composition. The right hand side of Eq. (1) will be small only if the pressure, temperature and chemical composition changes are small enough. In turn, the magnitudes of these variables are governed by the dynamic and chemical processes occurring in the flow field. The energy and momentum equations will play a major role in fixing the pressure and temperature. The species balance equations will determine the chemical composition in

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

338

the flow field. The advantage of writing the EOS in the form of Eq. (1) is that the flow field effects are isolated in Dp/Dt, DT/Dt and D…c /Dt while the thermodynamic character of the fluid is isolated in h, i and kc. The starting point for chemical reactor modeling is the list of the more fundamental equations, found in many well known textbooks (e.g. Bird et al., 1960, Kuo, 1986): Continuity: 1 Dz = −9 · 7 z Dt

(5)

Momentum: z

N D7 = −9p − 9 · ¬¯ + % zcgc Dt c=1

(6)

Net viscous stress: 9 · ¬¯ = − 9 · {v[97 +(97)T]}

!

+9



"

2 v− vB (9 · 7) 3

Energy: N

% zcCp,…c

c=1

(7)

 

DT T (zc = −9 · q− Dt zc (T N

Dp − (¬:9§) p,…c Dt

N

+ % (Jc · gc ) + % hc 9 · Jc

 

c=1 q

c=1

Rr,c + % ( −DHrr,c ) r = 1 M…c

(8)

Species mass: z

D…c = 9 · (zDc 9…c ) +Rc Dt

(9)

with zc = z…c. The boundary conditions at a certain reference location are 6i =(60,0,0), T =T0, zc =zc,0 and z = z0, and on the walls 6i = 0 and niqi =0 (or a known function) or T= TW (or a known function). Not all the terms in these equations have the same importance in determining the flow solution in chemical reactors. Physical considerations justify that a few modeling modifications can be made. First, the only body force considered in most reactor models, gi (per unit mass), is gravitation which is the same for all chemical species, g. The model equations for momentum and energy can then be simplified. In the momentum equaN tion N c = 1 zcgc = c = 1 zcg = zg. In the energy equation: N N c = 1(Jc · gc )=c = 1 Jc · g = 0. Second, in most multicomponent flows, the energy or heat flux contributions from the inter-diffusion processes are in general believed to be small and omitted in most applications, N c = 1 hc 9 · Jc : 0 (e.g. Slattery, 1990, p. 816; Kuo, 1986, p. 198; Bird et al., 1960, p. 566).

To determine whether temperature, pressure and chemical composition changes are small, one normally non-dimensionalizes the pressure, temperature and composition with scales that are governed by the dynamics of the reactive flow. An important question in reactive flow analysis is whether or not the density changes due to non-ideal mixing and chemical composition caused by reaction have significant effects on the flow. Considering the physical processes occurring in this problem, one expects that heat will be generated in the flow field by the heat of reaction, boundary temperatures at the walls, and heating/cooling devices rather than by viscous dissipation. Pressure changes may be due to chemical, thermo-dynamical and/or momentum effects. That is, both chemistry and flow scales may be important. For an ideal gas system many of the non-dimensional variables can be formed in a straightforward manner using the boundary values. In this way, one can define x*= xi /L, t*=t60/L, 6*= 6i /6o, z* = z/z0, h*=h/h0, i i i* = i/i0, k*= k /k , C * = CP /CP0, v*= v/v0, k*= c c c,0 P D /D Rc /Rc,0. The k/k0, Fr =6 20/ g L, D*= c c c,0 and R*= c temperature and pressure variables are more difficult to determine. Panton (1996) derived a suitable pressure variable based on the following arguments. In incompressible flow, pressure will play the role of a force in the momentum equation. Since pressure occurs as a gradient in this equation, a reference level may be substituted without any effect. That is, 9(p −p0)=9p. The first step in finding a pressure scale is to substitute the non-dimensional variable definitions into the momentum equation. Then, one argues that both the pressure and inertia terms would be needed in a typical incompressible flow problem. pS is used temporarily to symbolize the proper scale, that is, let p*= p−p0/pS. Substituting into the momentum equation, it can be seen that the pressure term will be of the same order as the inertia terms if the non-dimensional pressure is defined as p* =p− p0/z06 20. Thus, when pressure changes in the flow are dominated by momentum effects, this is the proper non-dimensional pressure variable. However, if chemistry mechanisms (e.g. volume and/or temperature changes with reaction) rather than momentum effects determine the pressure changes in the reactor, one must redo the non-dimensional scales for the pressure and reanalyze the results. A temperature variable is obtained using the difference between the initial and the adiabatic temperature as a characteristic temperature scale. The adiabatic temperature is an estimate of the maximum temperature (Tad = T0 + (−DHR )…c,0/CP0) that can be obtained, due to the heat generated by a representative reaction, in a reactor operating at adiabatic conditions (e.g. Froment & Bischoff, 1990, p. 409). Using this quantity the problem is supposed to have a characteristic scale: Tad − T0, which is a driving force for conduc-

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

tion of heat from the reactive zones into the bulk fluid field. Thus, one may define: T* = T− Tref/Tscale =T− T0/Tad − T0. When suitable problem dependent non-dimensional variables are substituted into the equations, the following non-dimensional groups occur: Re = z0L60/v0, Pr = v0CP0/k0, k0 = CP0/CV0, Fr =6 20/ g L, Ma =6o/CS, A = h0z0CP0T0, B=i0T0, Da =tflow/tchem =(L/60)/(z0/ Rc,0) (i.e. a Damko¨ hler number (Rosner, 1986, pp. 426–427)), and Q =(L/60)/z0CP,0(Tad −T0)qr= 1(Rr,c / M…c )(− DHrr,c ) (i.e. a dimensionless heat of reaction term (Rosner, 1986, p. 421)). For any given flow problem these parameters have specific fixed values. If they are large or small, they magnify or diminish the effect of the terms in which they appear as coefficients. In non-dimensional variables the mathematical problem is formulated as follows: Thermodynamics:





n



1 Dz* Dp* T −T0 DT* = k0Ma 2 h* −Bi* ad z* Dt* Dt* T0 Dt* +%k*k c c,0 c

D…c Dt*

(10)

The boundary conditions at the reference location are 6*i = (1,0,0), T*=1, z*c = 1 and z*=1, and on the walls 6*i = 0 and niq*i = 0. Inspection of these equations reveals the background for a few of the common reactor model assumptions. First, these systems are low velocity flows so that the terms with a pre factor containing the Mach number can be neglected. The density is not dependent on the pressure changes due to the flow, and the viscous dissipation and pressure terms in the energy equation can be neglected. Second, the density and thermodynamic coefficients are not constants. These quantities may be functions of both temperature and mixture composition. The governing equations for low Mach number flow derived based on the dimensional analysis can then be expressed as: (z + 9 · (z7)= 0 (t

(16)

((z6) + 9 · (z77)= − 9p − 9 · ¬¯ + zg (t

(17)

N

% zcCp…c

Continuity:

c=1

N (T + % zcCp…c7 · 9T (t c = 1 q

1 Dz* = −9* · 7* z* Dt*

(11)

D7* g = − 9*p* − 9* · ¬¯ * +z* Dt* Fr g

z*

(12)

Viscous stress: 9* · ¬¯ *= −

 !

1 9* · {v*[9*7* +(9*7*)T]} Re

+9*



2 v* − v*B (9* · 7*) 3

"

(13)

Energy: z*C*P

=9 · (u9T)+ %

 

r=1

Rr,c (− DHrr,c ) M…c

((z…c ) + 9 · (z7…c )= 9 · (zDm,c 9…c )+ Rc (t

Momentum:

339

(18) (19)

In accordance with the definition given by Panton (1996), the reactive flows in chemical reactors are characterized as low velocity compressible flows if the changes in temperature and composition significantly influence the flow. For dilute liquid systems the temperature dependence is probably most important, whereas both the temperature and composition dependencies may be important for gas systems. In this context many isothermal reactor systems may be considered incompressible, considering both gas and liquid systems.

2.2. Fixed bed dispersion models

DT* 1 1 = 9* · (k*9*T*) Dt* Re Pr +





The 1D axial dispersion model equations are given by (Froment & Bischoff, 1990): Mass balance:

k0 Ma 2 T0 v*b*viscous A Re Tad −T0

+ i*B

k0Ma A

2





T0 Dp* +T* +Q Tad −T0 Dt* (14)

(zc ((zc6z ) ( (… + = (zDea c )−rczcat(1−m) (t (z (z (z Enthalpy balance:

where v*b*viscous =(v/v0)(bviscous/(6 20/L)). Species mass:

zgCp

1 1 D…c = 9* · (z*D*9*… z* c c ) +DaR* c Dt* Re Sc

= (15)

(20)

(T (T +zgCp6z (t (z

 

( (T 4Ui u + % (− DHr, j )rjzcat(1−m)− (T−Ta) (z ea (z di j (21)

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

340

Momentum balance:

Initial conditions (t= 0):

dp z6 = −f dp dz

2 g z

(22)

Several parameterizations for the friction factor, f, are given in the literature. A relation (Eq. (23)) valid for spheres over a relatively broad range of particle Reynolds numbers is used in this work (Froment & Bischoff, 1990): f = 6.8

(1− m)1.2 − 0.2 Re m3

(23)

zc = zc0 T= T0

"

for all r and z

(29)

Boundary conditions (t\ 0) are given as: zc = zc0 T= T0

"

at z= 0

(zc (T = =0 (z (z

05 r5R

at z= L

05r5R

where Rep =dpz 6z /vL. Initial conditions (t = 0):

(zc =0 (r

at r= 0 and r= R

zc = zc0 Â Ã T =T0 Ì Ã pt = pt0 Å

(T =0 (r

at r= 0

for all z

(24)

for z=0

(zc dT = =0 (z dz

(zc ((zc6z ) + (t (z



for z =L

at r=R

all z

(34)

(26)

(zg ((zg6z ) =0 + (t (z







(27) Enthalpy balance: (T (T + zgCp6z (t (z





( (T 1 ( (T = u + u r (z ea (z r (r er (r + % (−DHr, j )rjzcat(1 − m) j

(33)

A 1D CFD model is formulated (e.g. Delhaye & Achard, 1976, 1977; Lahey & Drew, 1989): Continuity:

( (… 1 ( (… zDea c + rzDer c −rczcat(1 − m) (r (z r (z (r



(32)

(25)

Notice that the equation for pressure (Eq. (22)) will be identical for all the dispersion models used here as the cross sectional averaged velocity (i.e. superficial velocity) is used in the relation on the right hand side (i.e. giving rise to a linear pressure drop profile). For CFD models on the other hand, the full momentum equation is solved and the porous friction term will occur as a source term on the right hand side (as will be discussed later). The 2D dispersion model equations yield (Froment & Bischoff, 1990): Mass balance:

zgCp

all z

2.3. Fixed bed CFD models

zc = zc0 Â Ã T =T0 Ì Ã p= p0 Å

=

(31)

all z

(T Ui =− (T−Ta) (r uerzgCp

Boundary conditions (t \0) are given as:

(30)

(35)

Momentum balances: ( ( (zg6z )+ (zg6z6z ) (t (z =−

 

(p z 6z 6z ( (6 1 f −f +2 v z + zggz − zg D 6z 6z (z 2 (z dp (z di (36)

On macroscopic scales, the catalyst bed behaves as a porous media. Within a porous body the flow of a fluid is resisted by viscous and geometric (tortuosity) effects. A porous media friction term (− f(z 6 6z/dp)) is therefore added to the right hand side of the momentum equation. The volumetric wall shear force is parameterized as − 1/2zg( fD/di ) 6z 6z, where fD is the Darcy friction factor. Component mass balances:





(zi ((zi 6z ) ( (…i + = zD − rizcat(1−m) (t (z (z g ea (z

(37)

Enthalpy balance:



zgCp

(28)

=

(T (T +zgCp6z (t (z





( (T 4Ui u + % (− DHr, j )rjzcat(1−m)− (T−Ta) (z ea (z di j (38)

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

Initial conditions (t =0): initially there is no flow in the reactor (6z =6z0 = 0), the tube is filled with stagnant gas having a prescribed composition at a given temperature, as for the 1D dispersion model. Boundary conditions (t \0): uniform plug flow is assumed at the inlet. A prescribed pressure is specified at the reactor outlet. For the scalar variables Dirichlet boundary conditions are used at the inlet, whereas Neumann conditions are used at the other boundaries, as for the 1D dispersion model. A consistent 2D axisymmetrical CFD model is formulated in cylindrical coordinates (e.g. Delhaye & Achard, 1976, 1977; Lahey & Drew, 1989): Continuity: (zg ((zg6z ) 1 ((rzg6r) + + =0 (t (z r (r

+

   

(40)

and ( 1 ( ( (zg6z )+ (rzg6r6z ) + (zg6z6z ) r (r (z (t

+

 

 

(p z 6 6z 1 ( (6 ( (6 −f + zggz + rv z +2 v z (z dp (r (z r (r (z

 

1 ( (6 rv r (z r (r

(41)

Table 1 Tube data

Inner tube diameter (m) Outer tube diameter (m) Tube length (m) Temperature outside the wall (K) Heat coefficient for the metal (J/m s K)





(42) Enthalpy balance: zgCp

(T (T (T +zgCp6z + zgCp6r (t (z (r







( (T 1 ( (T u + uerr (z ea (z r (r (r



(43)

Initial conditions (t=0): initially there is no flow in the reactor (6r = 6r0 = 6z = 6z0 = 0), the tube is filled with stagnant gas having a prescribed composition at a given temperature, as for the 2D dispersion model simulations. Boundary conditions (t\ 0): the standard conditions for axisymmetrical flow in a 2D tube are applied (e.g. see Versteeg & Malalasekera, 1995). That is, no flow through the reactor wall. The normal velocity component is set to zero at the symmetry boundary. Uniform plug flow is assumed at the inlet. A prescribed pressure is specified at the reactor outlet. For the scalar variables Dirichlet boundary conditions are used at the inlet, whereas Neumann conditions are used at the other boundaries, as for the 2D dispersion model simulations.

2.4. Chemistry models, input data and parameter 6alues Syngas process

Methanol process

0.102 0.132 7 1100

0.016 0.026 0.15 523.2

52

52

Table 2 Catalyst data

Catalyst density (kg/m3) Particle diameter (m) Void fraction in the bulk region



j

( (6z v (z (r

=−



( (…i 1 ( (…i zD + rzgDer − rizcat(1−m) (z g ea (z r (r (r

=

+ % (− DHr, j )rjzcat(1−m)

(p z 6 6r 2 ( (6 ( (6 6 −f + rv r + v r −2v 2r (r r (r (z dp (r (z r

 

(zi ((zi 6z ) 1 ((rzi 6r) + + (t (z r (r

=

( 1 ( ( (zg6r)+ (rzg6r6r) + (zg6z6r) (t r (r (z =−

A porous media friction term is added to the right hand side of both components of the momentum equation. In the 2D case: Rep = dpz 6 /vL, where 6 is given by 6 = 6 2z + 6 2r . Component mass balances:

(39)

Momentum balances:

341

Syngas process

Methanol process

1500 0.0173 0.528

1775 0.0005 0.5

The kinetics, physical data, correlations and operating conditions needed for the synthesis gas conversion process are taken from Froment and Bischoff (1990), De Groote and Froment (1995), Xu and Froment (1989a,b), while the data for methanol conversion are taken from Vanden Bussche and Froment (1996). The operating conditions used in the simulations are given in the Tables 1–3. The data for the reactor tubes are given in Table 1. The data for the catalyst are given in Table 2. The data for the gas are given in Table 3. The reaction enthalpies are estimated by linear interpolation between the measured data obtained at T= 298 and 948 K (Xu & Froment, 1989a). The component viscosity is estimated using an expression found in Lydersen (1979). The heat capacities and conductivities for the individual components are calculated using ex-

342

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

Table 3 Gas data Syngas process

Mass fraction of CH4

0.1911

Mass fraction of CH3OH Mass fraction of CO

– 0.0001

Mass fraction of CO2

0.0200

Mass fraction of H2

0.0029

Mass fraction of H2O Mass fraction of N2

0.7218 0.0641

Superficial inlet velocity (m/s) Inlet temperature (K) Inlet pressure (bar)

1.89 793 29

The basis parameterizations used for the radial void fraction variations are taken from Govindarao and Froment (1986). However, to optimize on the model performance the void relations have been approximated by cubic spline interpolations (Faires & Burden, 1993).

Methanol process 0.1490 (6.00 mol%) 0 (0.00 mol%) 0.1734 (4.00 mol%) 0.2044 (3.00 mol%) 0.2564 (82.00 mol%) 0 (0.00 mol%) 0.2168 (5.00 mol%) 0.0067184 493.2 50

pressions found in Reid et al. (1988). The gas mixture values are estimated by a mole fraction weighted sum of the individual component contributions. The overall heat transfer coefficient, Ui, is calculated from the parameterizations given by Xu and Froment (1989b), Froment and Bischoff (1990). Mixing in the reactor due to the turbulence and the presence of packing is accounted for by introducing an ‘effective’ dispersive transport mechanism (Froment & Bischoff, 1990). The proportionality coefficients are the ‘effective’ diffusivities and conductivities. The effective axial diffusivity for the one-dimensional models was estimate based on experimental data given by Froment and Bischoff (1990) (p. 447). It was assumed that the Peclet numbers will be about 2 for these systems. The same value was used for all components. The effective axial conductivity for the one-dimensional model was estimated assuming that the ratio of the effective kinematic viscosity to the effective thermal diffusivity is of order unity, uea =zgCpDea (Bird et al., 1960). In the 2D simulations we adopt the correlations provided by Skaare (1993) for the effective axial and radial diffusivities, and for the effective axial conductivity: 1 1− 1 − m 1 = + Pema RepSc 2

(44)

1 1− 1 −m 1 + = Pemr RepSc 8.8[2 −(1 −2dp/di )2]

(45)

1 u 0 /u = er g + 0.7 Peha RepPr

(46)

The radial conductivity is estimated by (Xu and Froment): uer =u 0er + 0.14ugRepPr

(47)

3. Numerical scheme The simplified model equations (i.e. Eqs. (16)–(19)) given above reflect physical properties that simplify the construction of a suitable numerical solution algorithm. A split between the terms considering acoustic waves and the pressure gradient terms in the momentum equations will be introduced (e.g. Erlebacher, Hussaini, Speziale & Zang, 1992). In the flow calculations only the pressure gradient will be involved, whereas the absolute value of the thermodynamic pressure is needed in the heat (thermodynamics) and chemistry calculations. This split is beneficial as the momentum equations can be solved using larger time steps than allowed by the acoustic waves. This split between compressible flow and compressible fluid process calculations is an advantage as the flow part of the model can be solved by an incompressible an-elastic flow algorithm, whereas the compressibility due to an exposed external pressure or spatial static pressure variations can be imposed in a separate step. A similar split is made between the flowand chemistry parts as the algorithms that can tackle fluid flow problems may not be optimal solving stiff non-linear chemistry schemes which are coupled with flow and transport processes. The governing equations are written in conservative form and solved by a predictor–corrector method, a finite volume semi-implicit scheme. It is well known that the finite volume technique conserves the local heat, mass and momentum better than a simpler finitedifference approximation (e.g. Braza, Chassaing & Ha Minh, 1986; Leonard, 1994). A staggered grid for the velocity and pressure is used, as developed by Harlow and Welch (1965) in the marker and cell, MAC, method. The field variables are discretized using uniform cell sizes along the coordinate directions. The velocity components are specified at cell edges, while the scalar variables (e.g. pressure, densities, mass fractions) are specified at the scalar variable cell centers. This numerical method relies on the projection scheme originally developed by Chorin (1968) for the incompressible Navier–Stokes equations describing systems with constant fluid properties. Recently, several finite difference variants of this scheme have been proposed for variable density and reacting flows (e.g. Bell & Marcus, 1992, Najm et al., 1998, Almgren, Bell, Colella, Howell & Welcome, 1998, Knio et al., 1999). The present method is an alternative finite volume formulation which reduces to the finite volume ana-

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

logue to the simplified MAC method, SMAC (Amsden & Harlow, 1970a,b), in the case of constant density flows. The pressure correction step involves the inversion of a pressure Poisson equation, which is performed using a line-wise Gauss–Seidel iteration (TDMA) and a multigrid solver (Settari & Aziz, 1973; Hutchinson & Raithby, 1986; Kelkar & Patankar, 1989; Hutchinson, Galpin & Raithby, 1988; Sathyamurthy & Patankar, 1994). The fractional steps defining the algorithm are sketched in the following: 1. The Navier–Stokes equations are first split in an an-elastic step followed by an elastic step. The temporally an-elastic step can be formulated z n6*− z n6 n = − C n +D n −9p n +z ng Dt

(48)

z n6**−z n6* = − 9ƒ* Dt

(49) n

9 · (z n7**) −9 · (z n7*) = −92ƒ* Dt

n

(51)

to obtain the temporally an-elastic pressure correction equation − 9 · (z n7*) = −92ƒ* Dt

(52)

At this point both the pressure and velocity are corrected (i.e. calculating 6** and p*). 2. Integrate the transport equations for other scalar variables like T, zc and turbulent quantities by use of suitable discretizations and splitting algorithms (these are further described in a separate paragraph below). 3. A density update is included due to changes in temperature and fluid mixture composition.The density will be updated by the use of a suitable EOS. If an appropriate EOS is available, this procedure can be used for both gas and liquid phases. In the case of an ideal gas: z*=

p*M… n + 1 RT n + 1

= − 9(z*7**)

(54)

Erlebacher et al. (1992) used similar splitting steps performing large-eddy simulations of compressible turbulent flows. 5. A velocity update to account for changes in mixture density. By use of the product derivation rule the transient term in the momentum equation is rewritten as ((z6) (6 (z =z +6 (t (t (t

(55)

which implies that z n + 16***− z n6 n (6 (z = z n n + 6 n n (t (t Dt

(56)

(6 n 6**− 6 n : (t Dt

(57)

and (z n z n + 1 − z n : (t Dt

(58)

(50)

where the an-elastic filter has been used at the intermediate time step (z* = − 9 · (z n7**) = 0 (t

− z* (z* (z n = + = 0−9 · (z*7**) Dt (t (t

where

where ƒ*= p*−p , and C and D represent the numerical approximations of the convection and diffusion terms, respectively. n

z

343

n+1

(53)

4. The density at time step n +1 is obtained from the continuity equation. The elastic step is given by

6. The component densities are updated due to changes in temperature, pressure and mixture composition.At this point the mixture composition is given by the component mass densities, … nc + 1 =z*/ c c z*, c calculated by use of the component transport equations. With a given composition, the component mass densities are updated using the mass fraction definition: z nc + 1 = … nc + 1z n + 1. In this relation the mixture density variable is obtained from the mixture continuity equation, ensuring that the sum of component densities is equal to the mixture density. 7. The pressure is updated to fulfill the EOS, as the density obtained solving the continuity equation in one of the previous steps may not be consistent with the EOS.If an appropriate EOS is available, the following procedure can be used for both gas and liquid phases. In the case of an ideal gas p n+1 =

z n + 1RT n + 1 M… n + 1

(59)

8. Finally, the velocity is corrected due to the pressure update in the previous step (z n + 16 n + 1)− (z n + 16***) = − 9ƒ n + 1 Dt

(60)

where ƒ n + 1 = p n + 1 − p*. At this point only the velocity is corrected. A correction step similar to this one has been suggested by Knio et al. (1999).

344

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

For constant density fluids, only step 1 (and step 2 for reactive flows) is performed which corresponds to a time step in the SMAC algorithm. The scheme has been sketched with first order accuracy in time. In the time integration, an explicit first order Euler scheme has been used to advance the velocity and pressure fields. This can easily be improved by using a second-order Adams– Bashforth scheme. The numerical technique used solving the chemistry part of the model, step 2 in the algorithm sketched above, may need some further attention. The fractional step concept applied consists of successively applying predefined operators determining parts of the transport equation by intermediate time integrations. In many applications, the operators are used to represent different spatial directions, i.e. a complicated 3D problem can be solved by splitting it up into simpler 1D problems. However, the choice of operators is somewhat arbitrary, and in the present model the operators are selected on the basis of the physical and chemical processes according to the individual transport equations. The convective and diffusive terms are further split into their components in the various coordinate directions. The time-truncation error in the splitting scheme is of first order. However, if the order of applying the operators is reversed every time step, the splitting method itself can be shown to be second order accurate (Strang, 1968). Thus, if the individual operators applied are of second order (or higher order) in time, the total time integration procedure will be second order accurate. Applying either the first order or the second order version did not give any significant differences in the results presented in this paper. On each sub-step, solving the component transport equations, the diffusive terms are discretized with a second order central difference scheme in space, and a first order semi-implicit Euler scheme in time. The convective terms are solved using a second order TVD scheme in space, and a first order explicit Euler scheme in time. The TVD scheme applied was constructed by combining the central difference scheme and the classical upwind scheme by adopting the ‘smoothness monitor’ of van Leer (1974) and the ‘Superbee’ limiter of Roe (1986), see also Sweby (1984), Le Veque (1990), Toro (1999). The temperature advection is discretized using a second order upwind scheme in space, and a first order explicit Euler scheme in time. The chemical reaction part of the model is discretized by a first order implicit Euler approximation, and the resulting set of algebraic equations are solved using a Broyden subroutine (Press, Teukolsky, Vetterling & Flannery, 1992). The other source terms are solved using a first order explicit Euler approximation. Solving the component mass balances in the 2D CFD model, this method thus decouples the chemistry (i.e. kinetics) and the transport (i.e. convection and diffu-

sion) terms. The operators are successively integrated in the following order: z*− z nc c = Cz (z nc ) Dt

(61)

z**−z *c c = Cr(z*) c Dt

(62)

z***− z** c c = Dz (z**) c Dt

(63)

z****− z*** c c =Dr(z***) c Dt

(64)

z*****− z**** z nc + 1 − z**** c c c = = R(z****) c Dt Dt

(65)

where Cz(z nc ), Cr(z *), Dz(z **), Dr(z ***) and R(z ****) c c c c represent the numerical approximations of the axial convection, radial convection, axial dispersion, radial dispersion and reaction terms, respectively. Finally, as the method consists of successively applying the operators by intermediate time integrations, z nc + 1 = z *****. c The sum of these equations (or intermediate steps) reveals that: z nc + 1 − z nc = Cz (z nc )+ Cr(z*)+ Dz (z**)+ Dr(z***) c c c Dt +R(z****) c

(66)

The superscript stars indicate that the intermediate solution of the first operator is used in the second operator, and so on. Solving the heat balance in the 2D CFD model this method decouples the chemistry (i.e. heat of reaction), heating/cooling and the transport (i.e. advection and diffusion) terms. The operators are successively integrated in the following order: T* − T n = Az (T n) Dt

(67)

T** −T* = Ar(T*) Dt

(68)

T***− T** = Dz (T**) Dt

(69)

T****−T*** =Dr(T***) Dt

(70)

T*****− T**** =Q(T****) Dt

(71)

T******− T***** T n + 1 − T***** = = H(T*****) Dt Dt (72) n

where Az(T ), Ar(T*), Dz(T**), Dr(T***), Q(T****) and H(T*****) represent the numerical approximations of the axial advection, radial advection, axial dispersion, radial dispersion, heating/cooling and heat

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

of reaction terms, respectively. The sum of these equations (or intermediate steps) yields: T n+1 − T n =Az(T n)+Ar(T*) +Dz (T**) + Dr(T***) Dt +Q(T****) +H(T****)

(73)

4. Numerical analyses The performance and behavior of the algorithm, presented in Section 3, are assessed based on simulations of two industrial chemical processes (i.e. the synthesis gas and methanol production processes) performed in multi-tube fixed bed reactors. The fractional step algorithm is implemented in FORTRAN.

4.1. 1D dispersion model simulations A 1D pseudo-homogeneous dispersion model (Eqs. (20) –(26)) is used simulating the performance of the synthesis gas and methanol production processes. The results predicted by the fractional step method are evaluated and compared with simulations performed by a corresponding model implementation in Matlab. In the Matlab program, the integration in time is performed by use of the function ode15s, while discretization in space is performed by the dss020 and dss042 functions (Schiesser, 1991).

345

In this analysis prescribed velocity and mixture density values were used to enable proper solutions with the Matlab programme using the predefined functions. The results for the synthesis gas process are given in Fig. 1. In the synthesis gas process CO is produced from CH4 and H2O. The reversed water gas shift reaction consumes CO in the z-interval from 0 to 1.5 m. The mole fraction of CO2 is increasing from the reactor entrance and reaches a maximum at a distance of 1.5 m from the inlet. At z= 1.5 m the reaction is reversed, and a maximum occurs in the CO2 profile. Then there is a slight decrease to the point where equilibrium is achieved. The gradient of the CO-profile is not very steep at the reactor entrance. However, it gets steeper in the interval of 0.25–1 m. Further away from the entrance, the profile flattens until equilibrium is achieved about 6–7 m into the reactor. The agreement comparing the results from the two simulations is excellent. However, there is a slight difference in the temperature curves in the interval from 0 to 0.5 m. In the fractional step program, the time increment is fixed at 2× 10 − 3 s. The implicit nature and the automatic step-length control of the ode15s integrator, indicate that the Matlab solution may be more accurate than the one obtained by use of the fractional step algorithm. However, combined with a full CFD model an automatic step-length control will

Fig. 1. 1D mole fraction- and temperature profiles of the synthesis gas process are compared for the Matlab program and the fractional step program. Matlab results are plotted with a solid line, fractional step results with a dash-dotted line.

346

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

Fig. 2. The predicted 1D mole-fraction and temperature profiles for the adiabatic methanol process. Matlab results are plotted with a solid line, fractional step results with a dash-dotted line.

not be appropriate due to the increased computational costs involved. The results from the simulation of the adiabatic methanol process are given in Fig. 2. A comparison between these profiles and the simulated results reported by Vanden Bussche and Froment (1996) shows that the corresponding profiles are hardly distinguishable (not shown). This was expected, since all important model parameterizations used here are the same as the ones used by Vanden Bussche and Froment (1996). In this simulation the time increment in the fractional step program is 4× 10 − 4 s. The plots correspond very well to the results obtained by Vanden Bussche and Froment. All the trends are the same. Methanol is produced at the expense of CO, CO2 and H2. However, in the interval of 0– 0.005 m the amount of CO is still increasing, due to the water gas shift reaction. Then the reaction turns, and CO is consumed until the total equilibrium is achieved. The change in the gradients in all of the mole fraction plots at z= 0.005 m is due to the reversion of the water gas shift reaction. However, the match between the fractional step and the Matlab program results is not as good as for the synthesis gas process, although a smaller time increment was applied in the latter. The temperature in the fractional step simulation is approximately 2 K above the temperature in the Matlab simulation. This indicates that the methanol process is a stiffer problem than the synthesis gas process.

As a secondary effect of this temperature difference, there are also differences in some of the component profiles. Especially, the production of CH3OH as well as the consumption of CO are sensitive to variations in temperature. All the profiles have steep gradients at the reactor entrance. A notable change in these gradients occurs at z =0.005 m. Equilibrium is reached at z= 0.035 m. CO is produced in the interval from the reactor entrance to 0.005 m. Then, there is a considerable decrease in the amount in the interval of 0.005– 0.035 m. When cooling is introduced to the reactor, as presented in Fig. 3, a larger reactor is required to achieve equilibrium. A non-adiabatic reactor of length 0.4 m is simulated. The match between the fractional step and Matlab simulations is improved when cooling is introduced. In the interval of 0–0.1 m, however, the temperature from the fractional step simulation is still a few degrees (K) above the results of the Matlab simulation.

4.1.1. 2D dispersion model simulations 2D dispersion models, i.e. Eqs. (27)–(34), are solved by use of the fractional step algorithm. For the synthesis gas process, the results given in Fig. 4 show that both the temperature and mole fraction profiles are steep at the inlet but become rather flat after about 2 m within the reactor, in accordance with the 1D model predictions. In addition, a radial temperature gradient

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

with a maximum at the wall is also observed at the reactor entrance. Further away from the reactor entrance, the radial profile is flat. The mole fraction profiles also contain marked radial gradients within the first 1.5 m of the reactor. For the methanol process (Fig. 5), the temperature gradient at the reactor entrance is very steep. The temperature increases to about 550 K in the center of the tube, but near the walls the maximum temperature is approximately 530 K. The temperature decreases further away from the reactor inlet, and the radial profiles are almost flat at z =0.3 m downstream from the inlet. The effect of oscillatory void fraction distributions in packed beds were then investigated. Experimental measurements of Benenati and Brosilow (1962), Ridgway and Tarbuck (1968) among others indicate the presence of oscillatory radial variations of the void fraction in packed beds. These variations seem to be fairy welldefined and are due to the confining effect to the wall of the bed. To evaluate the importance of the radial variations in the void fraction profiles embedded in the 2D dispersion models, simulations without radial variations in the bed structure have been performed for both processes (not shown). No significant differences were observed comparing the corresponding results for both processes.

347

4.1.2. 2D dispersion model performance optimization In a first approach solving the model both the reaction terms in the component mass balances and the heat sources/sinks in the enthalpy balance were solved simultaneously using an implicit Euler discretization in time. The resulting set of algebraic equations were solved by a Broyden subroutine. Further analysis reveals that an additional operator split is beneficial. In the final approach solving the model, only the kinetics in the component mass balances were solved using the implicit approach described above. The heat sources/sinks in the enthalpy balance were solved by use of an explicit Euler discretization in time. Results for comparison of these two solution approaches show that the model predictions are hardly distinguishable (not shown). However, the time consumed by the final program version was only 5% of the previous one. To compare the computational costs of the different subroutines in the improved program, the time consumed by each operator was measured. The result is given in Fig. 6. Optimization of the integration method applied on the chemical reactions is recommended. The choice of an optimal solver is usually problem dependent. For reactor modeling, the accuracy required is considered to be in the range of three to four significant digits. This

Fig. 3. Comparison of Matlab and fractional step results predicting the 1D mole fraction- and temperature profiles of the non-adiabatic methanol process. Matlab results are plotted with a solid line, fractional step results with a dash-dotted line.

348

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

Fig. 4. 2D mole fraction- and temperature profiles of the synthesis gas process simulated by the fractional astep method. Oscillating variations in void fraction is included in the model. Note, as a 2D axisymmetric tube has been simulated, the model is solved for positive r-values only.

allows the application of low order methods with high efficiency, in preference to the more accurate but time consuming methods. Explicit methods do not need an iteration procedure to find the solution like implicit methods do. Hence, the most efficient methods for non-stiff problems are often the explicit methods. But, the application of these methods for stiff problems implies constrictions in the step-length, due to stability requirements. Implicit methods have a larger absolute stability region than explicit methods, which implies less stability restrictions. Hence a much larger step-length may be used. However, many implicit methods require that analytical Jacobians are available. In the processes investigated here as well as in many other practical examples, the kinetics are very complicated and we need approximate Jacobians. This implies that the matrix must be evaluated at every iteration and every time step, which can make the computations costly. The computation costs can be reduced, if it is possible to use the same Jacobian matrix in several iterations and over a few integration steps. For large chemistry models, the number of zeros in the Jacobian may be considerable. The sparsity can be exploited to reduce the cost of the linear algebra calculations.

The Runge –Kutta, Rosenbrock and linear multi-step methods may be used solving the present model equations. However, Hundsdorfer and Verwer (1995) showed that the direct use of multi-step methods like BDF for reactive systems give rise to amplified splitting errors and a drop in accuracy. The identification and control of splitting errors has thus been a common subject of investigation (e.g. Dabdub & Seinfeld, 1995; Saylor & Ford, 1995; Hundsdorfer & Verwer, 1995; Verwer, Blom & Hundsdorfer, 1996a; Jacobson, 1996; Jacobson, Tabazadeh & Turco, 1996; Knoth & Wolke, 1995, 1998; Sandu et al., 1997a; Sandu, Verweer, Blom, Spee, Carmichael & Potra, 1997b; Lanser & Verwer, 1999; Verwer & Simpson, 1995; Verwer, Blom, Van Loon & Spee, 1996b; Verwer, Spee, Blom & Hundsdorfer, 1999). Emphasis has been placed on the stability of different operator-splitting schemes and the role of stiffness and stability. In addition, in the context of differential algebraic equations (DAE) -systems (e.g. for systems involving both ordinary differential equations for kinetics and algebraic equilibrium equations), the lack of iteration makes the Rosenbrock methods only suitable for explicit integration schemes where the two types of equation systems are solved separately. The more stable

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

implicit integration schemes require a fully implicit integration method. Hence, for solving the chemistry equations in general, an implicit Runge– Kutta method combined with the Broyden’s method may be applied. A reasonable alternative is to use an implicit discretization scheme combined with the Broyden method, as applied in the present work.

5. Results and discussion The low tube to particle diameter ratio in the synthesis gas process indicate that oscillating radial void variations could occur over the entire cross section of the reactor. If such oscillations occur in the bed structure, they may have a significant influence on the reactor performance (Froment & Bischoff, 1990). However, as discussed above, the non-uniform void fraction distribution have no significant effect on the 2D dispersion model results. It can be argued that the standard dispersion models are not capable of resolving such effects, since the radial variations in the flow are averaged out. In this context the 2D dispersion models are considered incon-

349

sistent, as the basis for such a model is a cross sectional averaged model formulation (e.g. see the more fundamental model derivation given by Delhaye & Achard, 1976, 1977; Lahey & Drew, 1989). In addition, all the dispersion model simulations discussed above were performed adopting prescribed values for the superficial velocity and the mixture density throughout the reactor, and a linear pressure profile was used according to Eq. (22). The superficial velocity and the mixture density values were set to the inlet values. Therefore, in this section, the effect of physical variations in velocity, pressure and mixture density on the reactor performance will be investigated using a 1D CFD model. The influence of radial variations in the bed structure on the chemical systems will then be further analyzed with a 2D CFD model. For the methanol process, on the other hand, the tube to particle diameter ratio is much larger, and the non-uniform oscillating void fraction distribution are only expected to occur at a very thin zone close to the wall. An very fine grid resolution is needed in the radial direction to resolve these effects. No CFD simulations have therefore been run for the methanol process, as the computational costs will be prohibitive.

Fig. 5. 2D mole fraction- and temperature profiles of the non-adiabatic methanol process simulated using the fractional step method. Oscillating variations in void fraction is included in the model. Note, as a 2D axisymmetric tube has been simulated, the model is solved for positive r-values only.

350

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

Fig. 6. Computational time consumed by the different subroutines in percent of the total computation time. Black bars are for the synthesis gas process and white bars are for the methanol process. Definitions of mathematical operators: Convz: Convection in z-direction, Diffz: Diffusion (dispersion) in z-direction, Diffr: Diffusion (dispersion) in r-direction, Reaction: kinetics, Advecz, Enthalpy (temperature) advection in z-direction, Diffz – erg: heat conduction in z-direction, Diffr – erg: heat of reaction terms, and Q – energy: heating/cooling terms.

5.1. 1D CFD simulations A 1D CFD simulation (Eqs. (35)– (38)) was performed to evaluate the effect of the variations in velocity, pressure and mixture density on the synthesis gas reactor performance. The results are given in Figs. 7 and 8. The profiles obtained from the 1D CFD simulation are similar to the results obtained by the corresponding 1D dispersion model. The total pressure drop profile is almost linear justifying the assumption made regarding this variable for the simpler dispersion models. However, the axial gradients predicted in the velocity and mixture density profiles are significant. The superficial velocity increases from 1.89 m/s at the reactor inlet to about 3.4 m/s at the outlet. This effect is balanced by a drop in the mixture density from 7.89 to about 4.5 kg/m3. The net effect of these changes in superficial velocity and mixture density on the process is that a higher rate of conversion of CH4 and H2O to CO and H2 is obtained. There are no significant differences in the temperature profile.

5.2. 2D CFD simulations A consistent 2D axisymmetric CFD model (Eqs. (39) –(43)) is used to evaluate the influence of a non-

uniform void distribution in the reactor. Results from a 2D simulation with a uniform void fraction distribution are given in Figs. 9 and 10. A vector representation of the two velocity components are given in Fig. 11. The results obtained from the 2D simulations with a uniform void distribution (see Figs. 9–11) do not vary much from the results obtained from the 2D dispersion model. The radial gradients observed at the first 1– 2 m inlet zone for the chemical components in the reactor (Fig. 10) are due to the heat induced at through wall. However, the reactions are endothermic, and not enough energy is transferred through the wall to smooth out the profiles. The radial velocity component is very small, as can be seen from Figs. 9 and 11. The flow pattern can be characterized as plug flow. Results from a 2D simulation with a non-uniform void distribution profile are given in Figs. 12– 14. In Figs. 12–14 the effect of introducing a non-uniform void fraction distribution can be seen. The velocity and pressure profiles are significantly altered. Higher void fraction at the wall induces less friction from the bed, and the gas flow increases while the pressure drop decreases. As the void fraction profile has its maximum value close to the wall, the highest flux of mass are observed close to the wall. The mass

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

Fig. 7. Velocity, density, pressure and temperature profiles from the 1D simulation.

Fig. 8. Mole fraction profiles for components from the 1D simulation.

351

352

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

Fig. 9. Velocity, density, pressure and temperature profiles from the 2D simulation with a uniform void fraction profile. Note, as a 2D axisymmetric tube has been simulated, the model is solved for positive r-values only.

Fig. 10. Mole fraction profiles for the components from a 2D simulation with a uniform void fraction distribution. Note, as a 2D axisymmetric tube has been simulated, the model is solved for positive r-values only.

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

353

Fig. 11. Mole fraction profiles for the components from a 2D simulation with a uniform void fraction distribution. Note, as a 2D axisymmetric tube has been simulated, the model is solved for positive r-values only.

Fig. 12. Velocity, density, pressure and temperature profiles from a 2D simulation with a non-uniform void fraction distribution. Note, as a 2D axisymmetric tube has been simulated, the model is solved for positive r-values only.

354

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

Fig. 13. Mole fraction profiles for the components from a 2D simulation with a non-uniform void fraction distribution. Note, as a 2D axisymmetric tube has been simulated, the model is solved for positive r-values only.

entering the reactor in plug flow migrates from the tube center towards the wall in a very narrow inlet zone. The radial velocity component at the inlet is about ten times larger than the corresponding velocity simulated with a uniform void fraction distribution. Fig. 14 shows the velocity profiles in the reactor. The pressure drop through the reactor is also reduced considerably (only 1/6 of the results obtained for uniform void fraction distributions), as shown in Fig. 12. However, the influences on the temperature and the mole fraction profiles were small but noticeable.

Fig. 14. Vector representation of velocities from a 2D simulation with a non-uniform void fraction distribution. Note, as a 2D axisymmetric tube has been simulated, the model is solved for positive r-values only.

This latter observation was a little surprising as larger humps of the type shown in these figures have been observed performing experimental studies on packed beds (e.g. Froment & Bischoff, 1990). Froment and Bischoff (1990) also stated that changes in the velocity profile are practically never accounted for in previous work, because it immediately complicates the computation in a serious way. In addition, very few data are available to date and no general correlation could be set up for the velocity profile. Therefore, as these simplifying assumptions are embedded in the derivation of the effective dispersion coefficients, the effect of the velocity profile has not been accounted for in a consistent manner. This may introduce a serious deficiency in the present correlations for uer as both solid and fluid are involved in the heat transfer processes at the wall, in contrast with Der. The flux by effective conduction is considered to consist of two contributions, a dynamic part and a static part. The parameterizations for the dynamic part that is dependent on the flow conditions are thus based on the superficial velocity, rather than the local velocity at the wall, which may be beneficial making the correlations more accurate. Together with further experimental studies, this paper provides a 2D CFD model that can be applied performing detailed analysis for optimal reactor operation, scale-up and design. This is especially important for systems where relatively small improvements may have a large impact on the process economy.

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

5.2.1. Computational costs The computational costs are considerably increased comparing the CFD simulations with the dispersion model simulations. Due to stability considerations in the flow part of the model, a very small Dt has to be used. A typical time step is of the order Dt = 5 × 10 − 5. It was found that most of the computational time of the CFD calculations is spend by the Poisson solver (75%) for pressure and by the reaction calculations (15%). Both these operators are solved by implicit methods where iterations are needed to obtain solutions. Parallelization and performance optimization of these model parts are important issues to consider for reducing the computational costs.

6. Heat and mass budgets Global mixture and component mass balance monitors were implemented in the programme. For the mixture mass balance the sum of mass in the system at computation end, plus the sum of mass out of the system during integration, should be equal to the sum of initial mass in the system, plus the sum of mass entering the system during integration. The corresponding component mass balances contain two more terms, one term on the right hand side adding the mass transformed (i.e. mass sink) to other components by reaction, and one term on the left hand side adding the component mass generated (i.e. mass source) from other components by reaction. A similar enthalpy (or temperature) balance monitor was also implemented. For a given grid cell, the temperature difference between the initial value and the final value should be equal to the sum of the individual temperature modifications obtained from the various operators as advection in both axial and radial directions, heat conduction in space, heat of reaction, and cooling/heating. The maximum cell error obtained in the whole calculation domain is used as error monitor. In general, the mass- and enthalpy budget calculations performed in the program show that the transported quantities were balanced with the precision determined by the computer word length.

7. Conclusion A numerical method is tailored for the solution of a reduced set of model equations developed for the description of reactive flows in chemical reactors. The performance and behavior of the operator-split scheme are assessed based on simulations of two industrial chemical processes (i.e. the synthesis gas and methanol production processes) performed in multitube fixed bed reactors.

355

Both 1D and 2D pseudo-homogeneous dispersion models (i.e. heat and mass balances) with prescribed velocity, mixture density and pressure profiles are used describing the reactor performance. The predicted profiles for both the methanol and synthesis gas processes were in accordance with results reported in the literature. The 2D simulations also reveal that the temperature and mole fraction profiles are not sensitive to radial variations in the bed structure of the synthesis gas reactor. Both 1D and 2D CFD models were developed and implemented in FORTRAN to simulate multi-tube fixed bed reactors for the synthesis gas process. A 1D CFD simulation was performed to evaluate the effect of the variations in velocity, pressure and mixture density on the reactor performance. The changes in these variables significantly effects the chemical composition and conversion in the reactor. The 2D CFD model was used for further studies of the synthesis gas process, analyzing the influence of non-uniform void fraction distributions on the flow and chemical conversion. The non-uniform void fraction distribution induces a significant reduction in axial pressure drop and a much higher fluid velocity close to the wall. However, the influence on the temperature and mole fraction profiles was hardly noticeable. Mass- and enthalpy budgets were implemented for the heat and mass quantities in the program. These budgets show that the enthalpy and both the mixture and component mass balances are fulfilled in the simulations. The numerical CFD algorithm developed in this paper has been found to be both computationally stable and accurate. Computational efficiency analysis shows that the Poisson solver requires about 75% of the total computational time. The computational time spend on the chemistry part of the model is about 15% of the total time. About 10% of the cost used on the chemistry part is spend on the chemistry solver. To provide results within more reasonable time limits, parallelization and performance optimization of the model will be considered.

References Almgren, A. S., Bell, J. B., Colella, P., Howell, L. H., & Welcome, M. L. (1998). A conservative adaptive projection method for variable density incompressible Navier – Stokes equations. Journal of Computational Physics, 142, 1 – 46 Article No. CP985890. Amsden, A. A., & Harlow, F. H. (1970a). A simplified MAC technique for incompressible fluid flow calculations. Journal of Computational Physics, 6, 322 – 325. Amsden, A. A. & Harlow, F. H. (1970b). The SMAC method: a numerical technique for calculating incompressible fluid flows. Los Alamos Scientific Lab., Report No. LA-4370. Anderson, D. A., Tannehill, J. C., & Pletcher, R. H. (1984). Computational fluid mechanics and heat transfer. Washington/New York: Hemisphere Publishing Corporation/McGraw-Hill.

356

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357

Beam, R. M., & Warming, R. F. (1978). n Implicit factored scheme for compressible Navier –Stokes equations. AIAA Journal, 16 (4), 393– 402. Bell, J. B., & Marcus, D. L. (1992). A second-order projection method for variable-density flows. Journal of Computational Physics, 101, 334 – 348. Benenati, R. F., & Brosilow, C. B. (1962). Void fraction distribution in beds of spheres. American Institute of Chemical Engineering Journal, 8, 359 – 361. Berge, E., & Jakobsen, H. A. (1998). A regional scale multi-layer model for the calculation of long term transport and deposition of air pollution in Europe. Tellus, 50 (3), 205 –223. Bird, R. B., Stewart, W. E., & Lightfoot, E. N. (1960). Transport phenomena. New York: Wiley. Braza, M., Chassaing, P., & Ha Minh, H. H. (1986). Numerical study and physical analysis of the pressure and velocity fields in the near wake of a circular cylinder. Journal of Fluid Mechanics, 165, 79 – 130. Chorin, A. J. (1968). Numerical solution of the Navier –Stokes equations. Mathematical Computation, 22, 745–762. Dabdub, D., & Seinfeld, J. H. (1995). Extrapolation techniques used in the solution of stiff ODEs associated with chemical kinetics of air quality models. Atmospheric En6ironment, 29 (3), 403 –410. De Groote, A. M., & Froment, G. F. (1995). Reactor modeling and simulation in synthesis gas production. Re6iews in Chemical Engineering, 11 (2), 145 – 183. Delhaye, J. M. & Achard, J. L. (1976). On the use of averaging operators introduced in two-phase flow modeling. Proc. CSNI Specialists’ Meeting, Toronto, 3 –4 August 1976. In: Banerjee, S. and Weaver, K. R. (Eds.) Transient two-phase flow (Vol. 1) (pp. 5– 84) AECL, Canada. Delhaye, J. M. & Achard, J. L. (1977). On the use of averaging operators in two-phase flow modeling. In: O. C., Jr. & Bankoff, S. G. (Eds.) Thermal and hydraulic aspects of nuclear reactor safety, 1: light water reactors. ASME Winter Meeting, 1977. Erlebacher, G., Hussaini, M. Y., Speziale, C. G., & Zang, T. A. (1992). Towards the large-eddy simulation of compressible turbulent flows. Journal of Fluid Mechanics, 238, 155–185. Faires, J. D., & Burden, R. L. (1993). Numerical methods. Boston: PWS Publishing Company. Ferziger, J. H., & Peric´ , M. (1996). Computational methods for fluid dynamics. Berlin: Springer. Froment, G. F., & Bischoff, K. B. (1990). Chemical reactor analysis and design. New York: Wiley. Godunov, S. K. (1959). A difference method for the numerical calculation of discontinuous solutions of hydrodynamic equations. Mat. Sbornik, 47 (3), 271 –306 In Russian. Translated U.S. Joint Publications Research Service 1636 Connecticut Ave., N.W. Washington 25, D.C., JPRS: 7225, 1-61, 1960. Govindarao, V. M. H., & Froment, G. F. (1986). Voidage profiles in packed beds of spheres. Chemical Engineering Science, 41 (3), 533– 539. Harlow, F. H., & Welch, J. E. (1965). Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. The Physics of Fluids, 8 (12), 2182 –2189. Hundsdorfer, W., & Verwer, J. G. (1995). A note on splitting errors for advection-reaction equations. Applied Numerical Mathematics, 18, 191 – 199. Hutchinson, B. R., & Raithby, G. D. (1986). A multigrid method based on the additive correction strategy. Numerical Heat. Transfer, 9, 511 – 537. Hutchinson, B. R., Galpin, P. F., & Raithby, G. D. (1988). Application of additive correction multi-grid to the coupled fluid flow equations. Journal of Numerical Heat. Transfer, 13, 133–147. Jacobson, M. Z. (1996). Application of a sparse-matrix, vectorized gear-type code in a new air pollution modeling system. ZAMM Z. Angew. Mathematical Mechanics, 76 (S2), 333 –336.

Jacobson, M. Z., Tabazadeh, A., & Turco, R. P. (1996). Simulation equilibrium within aerosols and non-equilibrium between gases and aerosols. Journal of Geophysics Research, 101 (D4), 9079 – 9091. Jonson, J. E., Bartnicki, J., Olendrzynski, K., Jakobsen, H. A., & Berge, E. (1998). EMEP Eulerian model for atmospheric transport and deposition of nitrogen species over Europe. En6ironmental Pollution, 102 (S1), 289 – 298. Kelkar, K. M., & Patankar, S. V. (1989). Development of generalized block correction procedures for the solution of discretized Navier – Stokes equations. Computer Physics Communications, 53, 329 – 336. Klemp, J. B., & Wilhelmson, R. (1978). The simulation of three-dimensional convective storm dynamics. J. Atmos. Sci., 35, 1070– 1096. Kim, J., & Moin, P. (1985). Application of a fractional-step method to incompressible Navier – Stokes equations. Journal of Computational Physics, 59, 308 – 323. Knoth, O., & Wolke, R. (1995). Numerical methods for the solution of large kinetic systems. Applied Numerical Mathematics, 18, 211 – 221. Knoth, O., & Wolke, R. (1998). An explicit – implicit numerical approach for atmospheric chemistry — transport modeling. Atmospheric En6ironment, 32 (10), 1785 – 1797. Knio, O. M., Habib, N., & Wyckoff, P. S. (1999). A semi-implicit numerical scheme for reacting flow. Journal of Computational Physics, 154, 428 – 467. Kuo, K. K. (1986). Principles of combustion. New York: Wiley. Lahey, R. T. Jr, & Drew, D. A. (1989). The three-dimensional time and volume averaged conservation equations of two-phase flow. Ad6ances in Nuclear Science and Technology, 20, 1 – 69. Lanser, D., & Verwer, J. G. (1999). Analysis of operator splitting for advection-diffusion-reaction problems from air pollution modeling. Journal of Computational and Applied Mathematics, 111, 201 – 216. Leonard, B. P. (1994). Comparison of truncation error of finite-difference and finite-volume formulations of convection terms. Applied Mathematical Modelling, 18, 46 – 50. Le Veque, R. J. (1990). Numerical methods for conser6ati6e laws. Basel: Birkhauser Chapter 16. Lydersen, A. L. (1979). Fluid Flow and Heat Transfer. Chichester: John Wiley & Sons. MacCormack, R. W. (1969). The effect of viscosity in hypervelocity impact catering. AIAA -69-354, 1 – 7. MacCormack, R. W. (1976). An Efficient Numerical Method for Solving the Time-Dependent Compressible Navier-Stokes Equations at High-Reynolds number. NASA TM 73129, also in: Computing in applied mechanics, Applied Mech. Div., Vol. 18. Majda, A., & Sethian, J. (1985). The derivation and numerical solution of the equations for zero Mach number combustion. Combustion Science and Technology, 42, 185 – 205. Marchuk, G. I. (1974). In A. Arakawa, & Y. Mintz, Numerical methods in weather prediction. New York: Academic Press. Najm, H. N., Wychoff, P. S., & Knio, O. M. (1998). A semi-implicit numerical scheme for reacting flows I. Stiff chemistry. Journal of Computational Physics, 143, 381 – 402. Ogura, Y., & Phillips, N. A. (1962). Scale analysis of deep and shallow convection in the atmosphere. Journal of Atmosphere Science, 19, 171 – 179. Olendrzynski, K., Jonson, J. E., Bartnicki, J., Jakobsen, H. A., & Berge, E. (2000). EMEP Eulerian model for acid deposition over Europe. International Journal of En6ironment and Pollution, 14 (1– 6), 391 – 399. Panton, R. L. (1996). Incompressible flow (second ed.). New York: Wiley. Patankar, S. V. (1980). Numerical heat transfer and fluid flow. New York: Hemisphere Publishing Corporation.

H.A. Jakobsen et al. / Computers and Chemical Engineering 26 (2002) 333–357 Peyret, R., & Taylor, T. D. (1983). Computational methods for fluid flow. New York: Springer. Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. (1992). Numerical recipes in FORTRAN 77. The art of scientific computing, vol. 1 (second ed.). Cambridge University Press. Rehm, R. G. & Baum, H. R. (1978). The equations of motion for thermally driven, buoyant flows. Journal of Research of the National Bureau of Standards, 83 (3) 297 –308. Reid, R. C., Prausnitz, J. M., & Poling, B. E. (1988). The Properties of Gases and Liquids (fourth ed.). New York: McGraw-Hill Book Company. Ridgway, K., & Tarbuck, K. J. (1968). Voidage fluctuations in randomly packed beds of spheres adjacent to a containing wall. Chemical Engineering Science, 23, 1147–1155. Roe, P. L. (1986). Characteristic-based schemes for the Euler equations. Annual Re6iew of Fluid Mechanics, 18, 337–365. Rosner, D. E. (1986). Transport processes in chemically reacting flow systems. Boston: Butterworths. Sandler, S. I. (1999). Chemical and engineering thermodynamics (third ed.). New York: Wiley. Sandu, A., Verweer, J. G., Van Loon, M., et al. (1997a). Benchmarking stiff ODE solvers for atmospheric chemistry problems II: implicit vs. explicit. Atmospheric En6ironment, 31 (19), 3151 – 3166. Sandu, A., Verweer, J. G., Blom, J. G., Spee, E. J., Carmichael, G. R., & Potra, F. A. (1997b). Benchmarking stiff ODE solvers for atmospheric chemistry problems II: Rosenbrock solvers. Atmospheric En6ironment, 31 (20), 3459 –3472. Sathyamurthy, P. S., & Patankar, S. V. (1994). Block-correctionbased multigrid method for fluid flow problems. Num. Heat Transfer, Part B, 25, 375–394. Saylor, R. D., & Ford, G. D. (1995). On the comparison of numerical methods for the integration of kinetic equations in atmospheric chemistry and transport models. Atmospheric En6ironment, 29 (19), 2585 – 2593. Schiesser, W. E. (1991). The numerical method of lines. Integration of partial differential equations. San Diego, USA: Academic Press. Settari, A., & Aziz, K. (1973). A generalization of the additive correction methods for the iterative solution of matrix equations. SIAM Journal of Numerical Analysis, 10 (3), 506 –521. Skaare, S.H. (1993) Reaction and heat transfer in a wall-cooled fixed bed reactor, Dr. ing. Thesis, The Norwegian Institute of Technology. Skamarock, W. C., & Klemp, J. B. (1992). The Stability of time-split numerical methods for the hydrostatic and the non-hydrostatic elastic equations. Monthly Weather Re6iew, 120, 2109–2127. Slattery, J. C. (1990). Interfacial transport phenomena. New York: Springer.

357

Sod, G. A. (1978). A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws. Journal of Computational Physics, 27, 1 – 31. Steger, J. L., & Kutler, P. (1977). Implicit finite-difference procedures for the computation of vortex wakes. AIAA Journal, 15 (4), 581 – 590. Strang, G. (1968). On the construction and comparison of difference schemes. SIAM Journal of Numerical Analysis, 5 (3), 506 –517. Sweby, P. K. (1984). High resolution schemes using flux limiters for hyperbolic conservation laws. SIAM Journal of Numerical Analysis, 21 (5), 995 – 1011. Tapp, M. C., & White, P. W. (1976). A non-hydrostatic mesoscale model. Quarterly Journal of Royal Met. Society, 102, 277–296. Thompson, P. A. (1972). Compressible-fluid dynamics. New York: McGraw-Hill. Toro, E. F. (1999). Riemann sol6ers and numerical methods for fluid dynamics. A practical introduction (second ed.). Berlin: Springer. van Leer, B. (1974). Towards the ultimate conservation difference scheme II. Monotonicity and conservation combined in a second order scheme. Journal of Computational Physics, 14, 361–370. Vanden Bussche, K. M., & Froment, G. F. (1996). A steady-state kinetic model for methanol synthesis and the water gas shift reaction on a commercial Cu/ZnO/Al2O3 catalyst. Journal of Catalysis, 161 (156), 1 – 10 Academic Press Inc. Versteeg, H. K., & Malalasekera, W. (1995). An introduction to computational fluid dynamics. The finite 6olume method. London: Longman. Verwer, J. G., & Simpson, D. (1995). Explicit methods for stiff ODEs from atmospheric chemistry. Applied Numerical Mathematics, 18, 413 – 430. Verwer, J. G., Blom, J. G., & Hundsdorfer, W. (1996a). An implicit – explicit approach for atmospheric transport — chemistry problems. Applied Numerical Mathematics, 20, 191 – 209. Verwer, J. G., Blom, J. G., Van Loon, M., & Spee, E. J. (1996b). A comparison of stiff ODE solvers for atmospheric chemistry problems. Atmospheric En6ironment, 30 (1), 49 – 58. Verwer, J. G., Spee, E. J., Blom, J. G., & Hundsdorfer, W. (1999). A second-order Rosenbrock method applied to photochemical dispersion problems. SIAM Journal of Science Computers, 20 (4), 1456 – 1480. Xu, J., & Froment, G. F. (1989a). Methane steam reforming, methanation and water – gas shift: I. Intrinsic kinetics. American Institute of Chemical Engineering Journal, 35 (1), 88 – 96. Xu, J., & Froment, G. F. (1989b). Methane steam reforming, methanation and water – gas shift: II. Diffusional limitations and reactor simulations. American Institute of Chemical Engineering Journal, 35 (1), 97 – 103.