Multiscale Modeling and Solution Multiplicity in Catalytic Pellet

Oct 17, 2008 - To whom correspondence should be addressed. Phone: (312) 413-7743. Fax: (312) 413-7803. E-mail: [email protected]., †. University of ...
2 downloads 0 Views 2MB Size
Download PDF
8572

Ind. Eng. Chem. Res. 2008, 47, 8572–8581

KINETICS, CATALYSIS, AND REACTION ENGINEERING Multiscale Modeling and Solution Multiplicity in Catalytic Pellet Reactors Kedar Kulkarni,† Jeonghwa Moon,† Libin Zhang,† Angelo Lucia,‡ and Andreas A. Linninger*,† Laboratory for Product and Process Design, Department of Chemical and Bioengineering, UniVersity of Illinois, Chicago, Illinois 60607, and Department of Chemical Engineering, UniVersity of Rhode Island, Kingston, Rhode Island 02881

Transport and reaction phenomena in catalytic pellet reactors are often difficult to analyze because of coupling between heat and mass transport occurring at different space and time scales. To calculate the reactor concentrations and temperatures, it is necessary to account for the species reaction and transport occurring in the reactor bulk at the macroscopic level as well as the catalyst pellets at the microscopic level. The resulting approach yields a large system of nonlinear partial differential equations with multiple scales and solutions that are difficult to find numerically. In addition, the catalyst pellets may operate in multiple steady states for identical conditions. Conventional computational methods may entirely miss the multiplicity phenomenon at the catalyst pellet level and, as a result, may not correctly predict overall reactor yields. In this paper, we introduce two numerical techniques to address multiple scales and multiplicity in heterogeneous reaction models. The first method expands existing bisection with “shooting”; the second global method deploys orthogonal collocation oVer finite elements with niche eVolutionary algorithms. We also propose a new multiscale method entitled effectiVeness factor maps to expedite and simplify the numerical effort to solve transport and reaction phenomena at different length scales. 1. Introduction Many industrial reactors involve heterogeneous reaction kinetics of packed catalytic pellets in fixed-bed reactors. To predict the concentration and temperature profiles in a fixedbed catalytic pellet reactor, mass and energy conservation equations for the reactor bulk need to be solved simultaneously with heat and mass transport inside the catalytic pellets. The robust numerical solution of this multiscale problem of reactions in the bulk and inside the pellets is still challenging. One more issue further compounding the difficult multiscale reactor problem is the possibility of multiple steady states in the catalytic pellets. Catalyst pellets are capable of operating in multiple steady states for identical reaction conditions; each steady state corresponds to different concentration and temperature profiles in the catalyst pellet leading to different overall realized conversions.1 Thus, multiple reactor concentrations and temperature profiles are possible. It is thus extremely important to quantify all possible steady states and study their impact on the resulting reactor performance. In this paper, we consider the case of a simplified fixed bed catalytic pellet reactor with coupled reaction and transport. In a catalytic pellet reactor, the reactant is converted into the product as it traverses the reactor length. It is a distributed system because the system states, such as the concentrations and temperatures, vary spatially. In addition, the reactor is packed with a catalyst bed to enhance overall yield and selectivity. The reactant may diffuse deeply into the catalyst pellets, while simultaneously being consumed within the porous microstructure. The reaction and transport phenomena in the reactor thus occur at multiple length scales in the reactor bulk

at the macroscopic leVel, as well as in the catalyst pellets at the microscopic leVel. Discretization of the distributed mass and energy balances for the reactor bulk, as well as the catalyst pellets, yields a large system of nonlinear algebraic equations. There are numerical solution techniques to solve these coupled boundary value problems. However, ill-conditioning caused by the multiple length scales and solution multiplicity make many algorithms unsuitable. This article proposes numerical techniques to overcome the challenges posed by multiple scales and solution multiplicity. The paper is organized as follows. Section 2 briefly introduces the general transport and reaction equations. Section 3 introduces two novel methods to handle steady state multiplicity. Section 4 discusses problems in direct simulation of the reactor bulk with full consideration of the pellet kinetics and mass transfer limitation. A novel interpolation approach entitled effectiVeness factor maps will be presented as a robust and efficient approximation technique. Section 5 closes the paper with conclusions and future work. 2. Reactions and Species Transport in a Packed Bed Reactor This section develops general reaction and transport equations in the packed bed reactor as well as in the catalyst pellets. A simplified packed bed catalytic pellet reactor with a spherical catalyst pellet packing is shown in Figure 1. 2.1. Macroscopic Level. We consider the heterogeneous overall stoichiometry of the reaction in the fixed bed reactor, as shown in eq. (1). catalyst

A f B ∆H

* To whom correspondence should be addressed. Phone: (312) 4137743. Fax: (312) 413-7803. E-mail: [email protected]. † University of Illinois. ‡ University of Rhode Island.

(1)

The species and heat balance in the axial and radial directions for a cylindrical control volume of the catalytic packed bed reactor are given in eqs 2 and 3, respectively.2 The momentum

10.1021/ie8003978 CCC: $40.75  2008 American Chemical Society Published on Web 10/18/2008

Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008 8573

transfer in the packed bed of the reactor can be approximated by Darcy’s law as in eq 4. ∇(u bCA) ) ∇ · (DA ∇ CA) + rA species balance

(2)

bT) ) ∇ · (Kb ∇ T) + rA∆H heat balance FgCp ∇ (u

(3)

-∇P)

µ b u Darcy’s law km

(4)

CA|z)0 ) Ci, T |z)0 ) Ti inlet boundary conditions ∇zT|z)L ) 0, ∇zCA|z)L ) 0 outlet boundary conditions(5a) ∇rCA|r)0 ) 0, ∇rT|r)0 ) 0 reactor-axis symmetry ∇rCA|r)Rr ) 0, ∇rT|r)Rr ) UA(T - Tcool) reactor-mantle: no mass outlet and cooling(5b) The inlet and outlet boundary conditions are summarized in eq 5a. They include the inlet gas concentration,CA|z)0 ) Ci, and temperature, T|z ) 0 ) Ti, as well as negligible diffusive outlet mass,∇zCA|z)L ) 0, and heat flux,∇zT|z)L ) 0, of the gas. The boundary conditions for the reactor axis and reactor mantle are summarized in (5b). They include the symmetry for concentration and temperature in the reactor’s center axis,∇rCA|r ) 0 ) 0, ∇rT|r ) 0 ) 0. There is no mass transfer across the reactor wall,∇rCA|r ) Rr ) 0. Cooling links the reactor mantle temperature to the cooling temperature,∇rT|r ) Rr ) UA(T - Tcool). The reaction in the bulk is coupled with reaction and diffusion in the porous catalyst pellets. Mass transfer resistance at the surface of the catalyst pellets causes a concentration gradient across the boundary layer thus coupling the mass transfer with the reactions inside the catalyst.3 The total reaction rate is the volume integral over the entire pellet. In order to compute the total reaction rate, the concept of effectiveness factor is introduced. The internal effectiVeness factor, η, is the ratio of the total overall reaction rate,rA,T, to a hypothetical homogeneous first order reaction rate at pellet surface conditions, rA(CAs), given by the following Arrhenius expression:

[ (

rA(CAs) ) -Sakref exp -

E RTref

)]

Tref - 1 CAs T

(6)

The total overall reaction rate,rA,T, can then be written as rA,T ) η rA(CAs)

(7)

Reference conditions for the effectiveness factor are local pellet surface temperatures and concentrations; these may vary radially and axially throughout the reactor. The possibility to dynamically adjust these reference conditions will be shown as a critical advantage in the effectiVeness factor maps introduced in section 4. The effectiveness factor can be computed by solving the pellet scale problem described next. 2.2. Reaction and Diffusion in Catalyst Pellets. A single catalyst pellet of radius, R, can be treated as a porous medium. The steady state species and energy balances for diffusive transport and reaction inside the catalyst pellets are given in eqs 8-10.4 De∇2CA + rA ) 0 species balance

(8)

Ke∇2T + rA∆H ) 0 heat balance

(9)

CA |r)R ) CAs ) CA(bulk) surface concentration (10a) T |r)R ) Ts ) T(bulk) surface temperature

(10b)

∇CA |r)0 ) ∇T |r)0 ) 0 symmetry boundary condition (10c) When assuming infinite mass and heat transfer at the catalyst

pellet surface, the surface concentrations and temperatures can be given as Dirichlet boundary conditions in eqs 10a and 10b. It is easy to replace eq 10a and 10b with proper flux Neumann boundary conditions, but then a heat transfer coefficient must be given. Because of symmetry, mass and energy flux at the center of the catalyst pellet are zero as in eq 10c. System 8-10 represents a nonlinear PDE system for coupled heat and mass transfer in a spherical nonisothermal catalyst pellet. Cylindrical or slab pellet geometry models are similar in mathematical structure as described elsewhere.4 Assuming that the pellet properties De and Ke are constants and that the reaction constant k is a function of the temperature only, the system in eqs 8-10 can be written in terms of the dimensionless concentration, y ) CA/CAs, as well as the dimensionless activation energy, γ, the heat of reaction, β, and the Thiele modulus, φ, as follows: γβ(1 - y) d2y 2 dy - φ2y exp + )0 2 x dx 1 + β(1 - y) dx

(11)

dy )0 dx x)0

(12)

y |x)1 ) 1

(13)

(

)

The independent variable x is the dimensionless catalyst pellet radius,x ) r/R. The parameters γ, β, and φ can be expressed in terms of the bulk transport and reaction properties, as well as the pellet surface concentration,CAs, and temperature,Ts, as follows: γ) β)

dimensionless activation energy

-∆HDe CAs Ke Ts



φ)R

E RgTs

dimensionless heat of reaction

[

(

)]

kref E Tref exp -1 De RgTref Ts

(14) (15)

Thiele modulus (16)

Finally, the desired effectiveness factor, η, accounting for the overall consumption of the reactant can be obtained by solving eq 11 and taking the flux at the surface as shown in eq 17. η)

3 dy x)1 φ2 dx

(17)

3. Multiplicity in Catalytic Pellets The solution to the nonlinear system in eqs 11-13 gives the steady state concentration profile of reactant A as a function of the dimensionless pellet radius. For the same reactor or pellet conditions, multiple reaction rates corresponding to different conversions are physically possible. Steady state multiplicity in the catalyst pellet means that for identical conditions given by the dimensionless parameters γ, β and φ, the effectiveness factor can assume several distinct, yet physically meaningful solutions. Figure 2 displays η-φ relations over various ranges of γ, β, and φ.1 In some ranges of γ, β, and φ, multiple concentration profiles satisfy the same species transport and energy balances. The precise value of the effectiveness factors thus depends on which branch of the multiple concentration profiles is selected. For the values γ ) 30, β ) 0.6, and φ ) 0.2, there are three distinct solutions with three different effectiveness factors, whose values differ by several orders of magnitude (η1 ) 1.0639, η2 ) 9.8791, η3 ) 513.928). Given the dimensionless reaction and transport coefficients, γ, β, and

8574 Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008

Figure 1. Schematic depicting the multiscale nature of the simulation problem in the catalytic pellet reactor. The multiscale model accounts for transport and reaction occurring in the reactor bulk (macroscopic level), as well as in the catalyst pellets (microscopic level).

φ, it is desirable to identify all steady state catalyst pellet profiles physically possible for the boundary value problem in eqs 11-13. While a number of methods, such as finite difference, shooting, collocation, and spectral methods, apply to boundary value problems, few specifically address multiplicity of pellet profile. One example is a recent contribution from Stadtherr’s group that addressed boundary value problem with interval method. 10 Most iterative methods depend strongly on initial

guesses. Only the solution “nearest” to the initial guess is typically found; there is no systematic process to ascertain the existence and values of different solutions. 3.1. Complete Solutions to the Pellet Problem. This section discusses two methods that can identify all solutions to the pellet problem in practice. It is more concise to state that our method is “highly likely” to identify all solutions because kinetic problems are non-polynomial hard and because floating point precision on a digital computer is used. The first method is an interval-based global bisection method modeled on shooting. The second method uses global hybrid niche evolutionary algorithms combined with orthogonal collocation over finite elements (OCFE). A third method based on global terrain methods has been presented elsewhere.5 The section closes with a comparison of the methods’ performance and robustness. 3.2. Shooting-Based Global Bisection Method over Intervals. This subsection introduces a global shooting bisection method to determine all possible effectiveness factors and corresponding concentration profiles. Shooting methods solve the boundary value problem by accurate forward integration.1,6 However, the earlier methods did not address multiplicity. These methods start with an initial guess for the center concentration. In our adapted method, we “shoot” for the Thiele modulus, φ, as a function of the center concentration, y0, with eq 18. φmodel(y0) - φ ) 0

(18)

Here, φmodel(y0), is the value of the calculated Thiele modulus as a function of a center concentration, y0. In our global bisection method, the entire domain of possible center concentrations, y0, is divided into a fixed number of intervals. We then search each of the brackets successively using global bisection. This method uses only function evaluations without requiring gradients. Gradient-based methods that are potentially faster

Figure 2. η-φ curves from Weisz and Hicks1 depicting the existence of multiple steady states leading to multiple effectiveness factor values. For a specific Thiele modulus, φ, more than one overall reaction rate is possible.

Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008 8575

Figure 3. (a) Schematic depicting niche areas in a search space of two variables. (b) A flowchart summarizing the main steps of the hybrid niche algorithm including population learning, solution acceleration and niche tracking.

were tried but often failed because of the steepness of concentration profiles in the reaction regimes with multiplicity. 3.3. Global Hybrid Niche Evolutionary Methods with OCFE Discretization. A different strategy adopts orthogonal collocation over finite elements (OCFE) to solve the equation system (eqs 11-13). The following two subsections introduce OCFE with stochastic niche evolutionary methods to identify all solutions to the pellet problem. Deterministic approaches are described elsewhere.5,10,11 3.3.1. Global Discretization of the Pellet Problem by OCFE. Orthogonal collocation over finite elements (OCFE) is a global polynomial approximation method based on approximation of the actual solution of partial differential equations with piecewise continuous polynomials over the domain of the independent variable.4,9 The OCFE discretization of eqs 11-13, leads to the following nonlinear algebraic system: Ay + G (y) ) 0

(19)

By ) c

(20)

where y is the vector of dimensionless nodal concentrations. Eqs 19 and 20 are the transport and reaction equations suitably discretized at internal nodes of each finite element domain; G (y) is a nonlinear vector function given in eq 21.

(

2 [j] G (y[j] i ) ) -φ yi exp

γβ(1 - y[j] i ) 1 + β(1 - y[j] i )

)

i ) 1, ..., m j ) 1, ..., n

all solutions, and at the same time, is accelerated by the quadratically convergent Newton method. Niche evolutionary algorithms are stochastic search methods based on the principles of genetics and natural selection used to obtain multiple solutions to multimodal unconstrained optimization problems.12,13 One local solution is identified first using the principles of stochastic search. Once the first solution and its surrounding area are detected, a niche is defined excluding future iterations from that niche as shown schematically in Figure 3a. By tracking niches containing different local solutions, it is possible to identify all local extrema contained in the search space. 3.3.3. Population Learning Acceleration. Figure 3b summarizes the main steps of the hybrid niche algorithm. The unknown concentration profiles of the initial population are improved by crossover and mutation. Some candidate solutions in each generation are subjected to a population learning process, which performs only a few iterations with a local quasi-Newton method. The learning procedure thus accelerates the convergence to the local minimum of the candidate population. Once the cloud of candidate solutions coalesces close to a local minimum as determined by the criterion in eq 22, the local quasi-Newton method is launched to complete convergence starting from the currently best candidate solution. More implementation details can be found in Moon et al.13 ms g ε · N

(21) Equation 20 implements boundary conditions at x ) 0 and x ) 1, as well as the C0 and C1 continuity conditions at finite element boundaries. More details on the OCFE discretization of the spherical catalyst problem are discussed in Appendix. 3.3.2. Global Hybrid Niche Evolutionary Methods. Unfortunately, the discretized system (eqs 19 and 20) is difficult to solve, and most methods do not render multiple solutions. We next introduce a method capable of rendering all solutions where they exist, reliably and in reasonable computational time. We propose a hybrid algorithm to combine the advantages of robustness and performance. We use globally optimal multimodal niche algorithms for multiplicity and robustness.7 A population learning method uses a Newton-based algorithm to accelerate refined candidate solutions or solution clusters. The combined hybrid niche algorithm is thus globally optimal, tracks

N

with ms )

∑ sh(d, j) and sh(d, j) ) j)1

{

(22) 1-

( ) d σsh 0

2

if d < σsh otherwise

Here, ms sums the sharing functions sh(d, j) measuring a penalty for each individual candidate solution j as a function of its Euclidean distance, d, from the best candidate solution and a given niche radius, σsh. The parameter N is the population size, and ε is a scalar between 0 and 1, typically set to ε ) 0.75. Individuals in subsequent populations are directed away from the current niche through fitness sharing.7,12 Thus, multiple solutions are sequentially obtained using the global hybrid evolutionary methods with population learning. 3.3.4. Results. The interVal-based bisection method, as well as the global hybrid niche eVolutionary method, was tested over different sets of γ, β, and φ parameter ranges with multiple

8576 Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008 Table 1. Comparison of Performance in CPU Time, as well as the Effectiveness Factor Values for Global Hybrid Niche Evolutionary Methods and Interval Based Bisection Methodsa performance comparison

a

hybrid niche

bisection method

MATLAB fsolve

γ

β

φ

CPU time

η

CPU time

η

CPU time

20

0.8

0.2

13.95 s

53.98 s

0.3

0.5

16.95 s

30

0.6

0.1

15.64 s

40

0.3

0.3

15.72 s

149.38 26.67 1.04 22.39 5.43 1.192 874.63 76.12 1.01 122.1 8.7 1.07

1.18 s

30

155.95 27 1.04 23.66 5.44 1.192 818.32 67.33 1.01 113.6 8.47 1.07

62.14 s 67.85 s 51.87 s

η

1.04 1.39 s 1.192 1.18 s 1.01 1.32 s 1.07

All results were calculated on a Pentium IV computer (1.60 GHz with 1 GB RAM).

Figure 4. Comparison of pellet concentration profiles evaluated by OCFE with niche algorithms (bold lines) and interval based bisection methods (dashed lines). The solution profiles are in close agreement.

concentration profiles. Their performance was evaluated on the basis of the CPU time needed to obtain the concentration profiles and evaluate the effectiveness factors. Table 1 shows that both methods computed multiple effectiveness factors accurately and robustly. In contrast, the MATLAB equation solver suite “fsolve” converged only when initialized close to the solution and is not apt to identify all solutions. We also compared the calculated pellet concentration profiles and effectiveness factors with literature values and found a good match as shown in Figure 4. So far, our efforts demonstrated that the pellet problem alone was problematic because of ill-conditioning and solution multiplicity. The next section discusses further difficulties when trying to proceed from the pellet kinetics to the reactor level.

4. Multiscale Simulation Using Effectiveness Factor Maps The reactor bulk species and heat balance in eqs 2-5 are strongly coupled with the species and energy transport in the catalyst pellet given by eqs 8-10. To determine the profiles in the entire reactor, the coupled nonlinear PDE system has to be solved. The most direct approach to solve the reactor and pellet multiscale problem is to solve the coupled nonlinear equation system resulting from discretization of mass and energy balances of reactor bulk and the catalyst pellets either simultaneously or using a nested iterative algorithm. The simultaneous solution is typically impractical because good initial guesses for pellets

Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008 8577

Figure 5. Main steps of the nested iterative algorithm used to obtain bulk concentrations and temperatures given the bulk and pellet transport, and reaction parameters.

and reactor cannot easily be attained. In each major iteration of the two stage approach, the pellet subproblem is solved locally with OCFE using the current iterates of the reactor concentrations and temperatures. The solutions to the pellet problem are further used to generate new reactor state iterates. The resulting procedure takes a nested structure as depicted in Figure 5. The nonlinear system of equations in each stage can be solved using a trust region dogleg method.8,14 The brute force simultaneous or nested computation of the reactor and pellet scale transport and reaction phenomena is slow, but more importantly, convergence failure is common with state-of-the-art algorithms. While the two numerical techniques developed in earlier sections solve for all possible steady states of a single catalyst pellet, their naı¨ve integration with the reactor level is impractical. One solution approach to overcome these challenges is based on so-called effectiVeness factor maps. The basic idea is simple. We extend the method of global polynomial approximation already used in the pellet to parametrize the pellet solution space in terms of γ, β, and φ. By solving for parameter ranges of interest with the methods of section 3, a complete threedimensional map of the effectiveness factors can be constructed. When executing the calculations a priori, the challenging pellet problem can be decoupled completely from the bulk reactor level. Solving the pellet kinetics reduces to retrieving the effectiveness factor values from these precalculated solution maps with suitable polynomial interpolation, if necessary. The next section describes how we generate the effectiveness factor maps for the heterogeneous kinetics. While the method of effectiveness factor maps is introduced here to solve a heterogeneous packed bed reactor problem, the notion of introducing solution maps as a nonlinear surrogate might be generalizable for other hard to solve nonlinear multiscale problems in process systems engineering. 4.1. Generation of Effectiveness Factor Maps. Effectiveness factor maps are generated within anticipated operational ranges of the dimensionless parameters γ, β, and φ; these

possible operational ranges are estimated on the basis of the maximum expected conversion and temperature change in the reactor. Initially, all solutions for effectiveness factors for operating condition sets (γ, β, φ) are calculated by solving the system of eqs 11-13 using any of the two global numerical approaches. The effectiveness factors corresponding to intermediates γ, β, and φ are approximated through cubic spline interpolation over the solution map ηˆ ) f(γ,β,φ), as in eq 23. 2

η(γ, β, φ) =

2

2

∑∑∑η

ijkΦi(γ)Φj(β)Φk(φ)

(23)

i)1 j)1 k)1

where ηijk is the precalculated effectiveness factor value at the grid-point (γi,βj,φk) and Φi(γ), Φj(β), and Φk(φ) are suitable base functions like Lagrange polynomials written as follows: 2

Φi(γ) )

(γ - γl)

∏ (γ - γ ) , l)1 l*i

i

l

2

Φj(β) )

(β - βl)

∏ (β - β ) , l)1

j

l

l*j 2

Φk(β) )

(φ - φl)

∏ (φ - φ ) (24) l)1

k

l

l*k

Effectiveness factor values for parameters (γ, β, φ) not corresponding to the nodal values can be approximated with the interpolation formula of eq 23. The solution map can be interpreted as a nonlinear parametrization of the pellet transport and reaction mechanisms. Figure 6 depicts the lattice of sample points for the effectiveness factor solution map. 4.2. Multiscale Reactor Simulation. The bulk reaction and transport properties of a catalytic pellet reactor with 20 m length and 0.5 m diameter are displayed in Table 2. The physical domain of the reactor was discretized into 50 finite volumes and the catalyst pellet was discretized using three collocation nodes in a single finite element. The problem was solved with the direct nested approach, as well as with effectiVeness factor maps. The effectiveness factor map was generated using 500

8578 Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008

Figure 6. For each vertex, the effectiveness factor values are precalculated with the methods of section 3. The effectiveness factors for any state (γ, β, φ) inside the cubes are approximated by cubic spline interpolation over the eight vertex values. Table 2. Bulk and Pellet Properties for the Catalytic Pellet Reactor with a Single Final Reactor Concentration and Temperature Profile properties of the catalyst pellets

properties of the bed

De ) 8.1 × 10-12 m2/s Ke ) 2.63 × 10-9 W/m.K R ) 0.003 m  ) 0.4 Fc ) 2.8 × 103 kg/m3

length ) 20 m, diameter ) 0.5 m DA ) 6.55 × 10-5 m2/s Kb ) 3.84 × 10-3 W/m.K b ) 0.4 Fb ) Fc(1 - b)

reaction properties

properties of inlet gas

E ) 50 000 J/mol ∆H ) -96 650 J/mol

Ci ) 2 mol/m3 Ti ) 700 K

points within γ, β, and φ ranges of γ ∈ [8,9], β ∈ [0,0.8], φ ∈ [1,2] and yielded effectiveness factors in the range η ∈ [0.806,4.468] as shown in Figure 7. Figure 8 shows that the simulation results using both methods match really well. The reactant concentration profiles show a monotonically decreasing trend along the length of the reactor because of its conversion into the product. because the reaction that occurs inside the reactor is exothermic, the temperature profile along the length of the reactor shows a hotspot. Figure 9 summarizing the performances of both methods indicates that the effectiveness factor map approach was 3-5 times faster than the nested collocation method in our experiments. This comparison also does not take into account that the nested approach often failed. These results demonstrate the superiority of the solution maps over simultaneous methods even for single solution problems. The next subsection deploys effectiveness factor maps to address

Figure 8. Comparison of simulation results for the bulk concentrations and temperatures in a catalytic pellet reactor (absolute error mean for the dimensionless concentration C ) 0.03; absolute error mean for the dimensionless temperature T ) 0.0072).

Figure 9. Performance comparison of nested collocation and 3D interpolation methods in terms of CPU time and number of finite volumes in the reactor. All results were computed on a Pentium IV computer (1.60 GHz with 1GB RAM).

Figure 7. Effectiveness factor map in the reactor’s expected γ, β, and φ ranges. The ranges chosen in this study were γ (8-9), β (0-0.8), and φ (1-2).

the problem of multiplicity, which is not tractable with traditional methods. 4.3. Solution Maps Make Multiplicity Tractable. The effectiveness factor maps developed in the previous section were also used to simulate the concentrations and temperatures in a catalytic pellet reactor for parameter ranges of (γ, β, φ) with multiple steady states. In these cases, the overall conversion in the reactor is not unique. For rigorously computing the reactor

Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008 8579

Figure 10. Effectiveness factor maps corresponding to regions of multiplicity for the ranges γ ∈ [20,22], β ∈ [0.1,0.9], and φ ∈ [0.4,0.7]. (a) All 500 solution points in this cube are located on the upper effectiveness factor solution branch (branch A). (b) The solutions encoded in this alternative effectiveness factor map correspond to the lower branch (branch B). Some regions in which the solution for η is unique are common in both maps.

Figure 11. Multiscale simulation using effectiveness factor maps corresponding to the upper and lower effectiveness factor branches (a) Simulation of reactor concentrations the upper effectiveness factor branch led to an overall conversion of almost 100%. (b) Simulation of reactor concentrations with the lower effectiveness factor branch led to an overall conversion of only 75% for the same reactor conditions.

performance in the presence of multiplicity, there are again at least two options. The direct method would solve simultaneously for the reactor bulk level, while identifying exhaustively all pellet solutions. Unfortunately, the direct simultaneous approach often diverges, even when solutions are unique as reported above. The EffectiVeness factor solution map addresses this tough multiscale problem. We generate a priori multiple solution maps tracking alternative multiplicity ranges in the dimensionless parameter space. To calculate pellet kinetics, solutions were retrieved from different maps as needed. Figure 10a and b shows two alternative effectiveness factor maps for heterogeneous kinetics in the parameter ranges γ ∈ [20,22], β ∈ [0.1,0.9], and φ ∈ [0.4,0.7]. The effectiveness factor cubes are different only in the ranges of multiplicity, but they are identical in parameter ranges where there is a single solution. This makes it convenient to apply as shown in these figures. The results demonstrated in the figures also show that the solution map approach is able to tackle this very tough multistate multiple length scale problem. The overall achievable conversions for the reactor corresponding to multiple states of the pellet kinetics are shown in the concentration profiles of Figure 11a and b The bulk reaction and transport properties of this catalytic pellet reactor with 6 m length and 0.5 m diameter are displayed in Table 3. Accordingly, the conversion was almost complete when pellets follow branch A, with η ∈ [1.01,378.28]. On the other hand, only 75%

Table 3. Bulk and Pellet Properties for the Catalytic Pellet Reactor with Multiple Final Concentration and Temperature Profiles properties of the catalyst pellets

properties of the reactor/bed

De ) 5.2 × 10-6 m2/s Ke ) 8.39 × 10-4 W/m K R ) 0.003 m  ) 0.4 Fc ) 2.8 × 103 kg/m3

length ) 6 m, diameter ) 0.5 m DA ) 6.55 × 10-5 m2/s Kb ) 3.84 × 10-3 W/m K b ) 0.4 Fb ) Fc(1 - b)

reaction properties E ) 125 000 J/mol ∆H ) -96 650 J/mol

properties of inlet gas Ci ) 2 mol/m3 Ti ) 700 K

conversion was achieved in the same reactor, for the exact same reaction conditions when pellets followed branch B, with η ∈ [1.01,277.68]. This difference in outcome would certainly be an unpleasant surprise in a real reactor. The present case study assumed that all pellets chose to follow consistently the same effectiveness factor branch. However, it is possible to imagine that different pellets adopt different steady states in the same reactor. So, the pellets in one section could operate on branch A, while pellets in a different section could operate on branch B. If this situation occurred in practice, the outcomes of the reactors would be chaotic, despite deterministic kinetic rate laws. Although we have no experimental evidence for such phenomena at this

8580 Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008

point, these observations about the unpredictability as a consequence of state multiplicity could have implications for chaotic behavior and accidents. These results show that three-dimensional effectiveness factor maps are preferable over direct solution for the following reasons: (i) The effectiveness factor map completely decouples the multiple length scales of the reactor and the pellet level. (ii) It eliminates the convergence-failure problems associated with simultaneous methods. (iii) The proposed multiscale approach may be the only method to reliably calculate concentration and temperature profiles in a heterogeneous catalytic pellet reactor with operating conditions favorable for multiple catalyst pellet steady states. 5. Conclusions and Future Work This paper presents numerical techniques to successfully incorporate a microscopic catalytic pellet model into macroscopic bulk reaction and transport for the catalytic pellet reactor. The proposed numerical techniques capture phenomena occurring at multiple length scales. We introduced two techniques to address state multiplicity: global hybrid niche evolutionary methods with orthogonal collocation over finite elements (OCFE) and an interval-based bisection method. The interval based bisection method is typically 3-4 times slower than global hybrid niche method. However, global hybrid niche evolutionary methods are very robust and failed in none of our experiments. We extended the polynomial approximation used for solving the pellet PDE to create global effectiVeness factor solution maps. Intermediate effectiveness factors are approximated by cubic spline interpolation in a precalculated effectiveness factor map constructed on known γ, β, and φ ranges. The effectiveness factor maps decouple the pellet length scale from the reactor level. Finally, it should be noted that discretization techniques for partial differential equations depend on the nodal choices. An issue that we did not address in this article is the proper choice of finite element spacing. However, we refer the interested reader to work by Cuthrell and Biegler.9 For completeness, we should find out other techniques to solve multiscale problems by parametrization. For example, in the context of optimization, parametrization is one of the basic ideas used in dynamic programming, and it is also the basis of the parametric programming approach developed by Pistikopoulos.15 Appendix I The OCFE approximation to the boundary value problem in system 11-13 divides the dimensionless radius domain, x ∈ [0,1], into n finite elements. Within each element, suitable orthogonal base functions, such as Lagrangian polynomials, li[j](x), depicted in Figure A-1approximate the true profiles, such that the integral approximation error over the element is minimized.16,17 For details on node selection, see Villadsen and Michelsen.4 The approximate polynomial concentration profile in each element, j, can be written as a sum of base polynomials with unknown coefficients, yi[j], as shown in eq 25. m

y[j](x) )

∑l

[j] [j] i (x)yi

i)1

m

j ) 1, ..., n where l[j] i (x) )

x - x[j] k

∏x k)1

[j] [j] i - xk

k*i

(25)

Figure A-1. Division of the spherical catalyst pellet (x ) 0-1) into n finite elements. Each element is further divided into m nodes corresponding to the roots of the mth order Jacobi polynomial. The x-axis denotes the dimensionless radius (x) and the y-axis stands for the dimensionless concentration (y).

Insertion of the polynomial concentration profile for the j-the element given by eq 25 into eq 11 gives eq 26 d2 [j] 2 d [y (x)] + [j] [y[j](x)] - φ2y[j](x) × 2 dx xi dx

(

exp

)

γβ(1 - y[j](x)) ) 0, i ) 2, ..., m - 1, j ) 1, ..., n(26) 1 + β(1 - y[j](x))

After insertion of the first and second derivatives of the base polynomials, eq 27 holds for the interior of each element. Equations 28 and 29 enforce zero and first order continuity between elements to ensure smooth transitions. m

∑l

m

[j] [j] [j] i ′′(xi )yi +

i)1

(

exp

∑

2 [j] [j] 2 [j] l[j] i ′ (xi )yi - φ yi × x[j] i)1 i

γβ(1 - y[j] i ) 1 + β(1 - y[j] i )

)

) 0, i ) 2, ..., m - 1, j ) 1, ..., n(27)

), j ) 1, ..., n - 1 C0 continuity (28) y[j](xm[j]) ) y[j+1](x[j+1] 1 d [j] d [y (x)]x)xm[j] ) [y[j+1](x)]x)x1[j+1], j ) 1, ..., n - 1 dx dx C1continuity(29) Finally, the boundary condition in eq 12 for the first node in the first element, x ) 0, can now be written as shown in eq 30. Similarly, the last node of the last element of the boundary condition in eq 13 becomes eq 31 for the final node, x ) 1: m

∑l

[1] [1] i ′ (x) |x)0yi ) 0

(30)

i)1

m

∑l

[n] [n] i (x) |x)1yi ) 1

(31)

i)1

Thus, we finally have a square system of n(m - 2) + 1 + 1 + (n - 1) + (n - 1) ) nm equations for the nm unknown coefficients of the base polynomials. To identify all sets of Values, yi[j], that solve the equation sets 27-31 simultaneously is the objective of the hybrid niche evolutionary methods.

Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008 8581

Nomeclature CA ) concentration ofreactant A inside the catalyst pellet CAs ) concentration of reactantA at the surface of catalyst pellet Cp ) specific heat of the gas d ) Euclidean distance betweencandidate and best solution DA ) bulk diffusivity of thereactant A De ) effective diffusivity inside the catalyst pellet E ) activation energy ∆H ) heat of reaction km ) permeability of reactant gas kref )reference reaction constant Kb ) bulk thermal conductivity in the reactor Ke ) effective thermalconductivity inside the catalyst pellet l ) Lagrange polynomial m ) number of nodes in each finite element ms ) sum of sharing functions n ) number of finite elements N ) population size P ) pressure in the reactor rA ) Arrhenius reaction rate rA,T ) Arrhenius total reaction rate per unit volume R ) catalyst pellet radius Rg ) universal gas constant Rr )reactor radius sh ) sharing function Sa )available surface area per kilogram of the catalyst T ) temperature inside the catalyst pellet Tref ) reference temperature Ts ) temperature at the surface of catalyst pellet x ) dimensionless radius of the spherical catalyst pellet xi[j] ) ith collocation node of the jth finite element y ) dimensionless concentration along radial direction of catalyst pellet yi[j] ) dimensionless concentration value ith collocation node of the jth finite element b u ) superficial velocity of the inlet gas Greek Symbols β ) dimensionless heat of reaction γ ) dimensionless activation energy  ) scalar value between 0 and 1 φ ) Thiele modulus Φ ) Lagrange polynomial to interpolate for effectiveness factor η ) effectiveness factor µ ) viscosity of reactant gas

Fg ) density of the gas σsh ) niche radius

Literature Cited (1) Weisz, P. B.; Hicks, J. S. The behavior of porous catalyst particles in view of internal mass and heat diffusion effects. Chem. Eng. Sci. 1962, 17, 265. (2) Fogler, H. S. Elements of Chemical Reaction Engineering; PrenticeHall: Upper Saddle River, NJ, 1999. (3) Froment, G. F.; Bischoff, K. B. Chemical Reactor Analysis and Design; John Wiley and Sons: Singapore, 1990. (4) Villadsen, J. V.; Michelsen M. L. Solution of Differential Equations by Polynomial Approximation; Prentice-Hall: New York, 1978. (5) Lucia, A.; Gattupalli, R.; Kulkarni, K.; Linninger, A. A barrier terrain methodology for global optimization. Ind. Eng. Chem. Res. 2008, 47, 2666. (6) Lee, J.; Kim, D. H. An improved shooting method for computation of effectiveness factors in porous catalysts. Chem. Eng. Sci. 2005, 60, 5569. (7) Beasley, D.; Bull, D. R.; Martin, R. R. Sequential niche technique for multimodal function optimization. EVol. Comput. 1990, 1, 101. (8) Biegler, L.; Grossmann, I. E.; Westerberg, A. W. Systematic Methods of Chemical Process Design; Prentice Hall PTR: Upper Saddle River, NJ, 1997. (9) Cuthrell,; J, E.; Biegler, L. T. On the optimization of differentialalgebraic process systems. AIChE J. 1987, 33, 1257. (10) Lin, Y.; Enszer.; J, A.; Stadtherr, M. A. Enclosing all solutions of two-point boundary value problems for ODEs. Comput. Chem. Eng. 2007, 32 (8), 1714–1725. (11) Piela, L.; Kostrowicki, J.; Scheraga, H. A. The multiple-minima problem in conformational analysis of molecules. Deformation of the potential energy hypersurface by the diffusion equation method. J. Chem. Phys. 1989, 93, 3339. (12) GoldbergD. E. RichardsonJ. Genetic algorithms with sharing for multimodal function optimization Proceedings of the Second International Conference on Genetic Algorithms; Lawrence Erlbaum Associates: Philadelphia, 1987; p 41. (13) Moon, J.; Zhang, L.; Linninger A. A. A hybrid sequential niche genetic algorithm for multimodal objective functions. Presented at the INFORMS Midwest Regional Conference, Evanston, IL, 2007. (14) Tang, W.; Zhang, L.; Linninger, A. A.; Tranter, R.; Brezinsky, K. Solving Kinetic Inversion Problems via a Physically Bounded Gauss-Newton (PGN) Method. Ind. Eng. Chem. Res. 2005, 44, 3626. (15) Pistikopoulos, E. N.; Georgiadis M. C.; Dua V., Multi-Parametric Programming: Theory, Algorithms, and Applications; Wiley-VCH: Weinheim, Germany, 2007. (16) Zhang, L.; Linninger, A. A. Temperature collocation algorithm for fast and robust distillation design. Ind. Eng. Chem. Res. 2004, 43, 3163. (17) Kulkarni, K. PhD dissertation, University of Illinois at Chicago, Chicago, IL, 2007.

ReceiVed for reView March 10, 2008 ReVised manuscript receiVed June 26, 2008 Accepted July 16, 2008 IE8003978