SciELO - Scientific Electronic Library Online

 
vol.77 número163MODELO DE PLANEAMIENTO DE LA EXPANSIÓN DE LA GENERACIÓN CONSIDERANDO RESTRICCIONES DE EMISIONESINFLUENCIA DE LA RESISTENCIA DE FALLA Y LA RESISTIVIDAD DEL SUELO EN LOS MÉTODOS DE LOCALIZACIÓN DE FALLAS BASADOS EN LA ESTIMACIÓN DE LA IMPEDANCIA. UN ANÁLISIS COMPARATIVO índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Servicios Personalizados

Revista

Articulo

Indicadores

Links relacionados

  • En proceso de indezaciónCitado por Google
  • No hay articulos similaresSimilares en SciELO
  • En proceso de indezaciónSimilares en Google

Compartir


DYNA

versión impresa ISSN 0012-7353versión On-line ISSN 2346-2183

Dyna rev.fac.nac.minas v.77 n.163 Medellín jul./sep. 2010

 

MAPPING THE REGION OF INSTABILITY FOR ADIABATIC PACKED BED REACTORS USING A HOMOTOPY CONTINUATION METHOD

DETERMINACIÓN DE LA REGION DE INESTABILIDAD PARA REACTORES ADIABATICOS DE LECHO EMPACADO UTILIZANDO UN METODO DE CONTINUACION POR HOMOTOPIA

 

JUAN PABLO GUTIÉRREZ-HERNANDEZ
Ingeniero Químico, M. Sc., Universidad Nacional de Colombia, Sede Manizales,J.P.Gutierrez-Hernandez@tue.nl

JAVIER FONTALVO ALZATE
Profesor Asociado, Universidad Nacional de Colombia, Sede Manizales, jfontalvoa@unal.edu.co

MIGUEL ÁNGEL GÓMEZ GARCÍA
Profesor Asociado, Universidad Nacional de Colombia, Sede Manizales, magomez@unal.edu.co

 

Received for review September 18th, 2009, accepted March 2th, 2010, final version April, 29th, 2010

 


ABSTRACT: The pioneer schematic ideas of Kimura and Levenspiel (Ind. Eng. Chem. Proc. Des. Dev., 16 (1977) 145 - 148) have been developed to find numerically the region of instability for adiabatic packed bed reactors. Three different cases of special industrial interest and complexity are presented. The highly exothermic gas-phase reactions: ammonia synthesis, methanol production from syn-gas, and SO2 oxidation. Equations were parameterized and solved according to a continuation homotopy numerical method. The results showed that concentration of inerts and total pressure influences the size of the instability region.

KEYWORDS: Instability region, adiabatic packed bed reactors, homotopy continuation method.

RESUMEN: Las pioneras esquemáticas ideas de Kimura y Levenspiel (Ind. Eng. Chem. Proc. Des. Dev., 16 (1977) 145 - 148) han sido desarrolladas para determinar numéricamente, para reactores de lecho empacado adiabáticos, cuándo un punto de operación es o no estable, y para localizar la envolvente de las condiciones a las cuales se debe evitar su operación. Tres reacciones en fase gas, de especial interés industrial y complejidad, altamente exotérmicas, se utilizan para ilustrar el método de análisis propuesto: la síntesis de amoníaco, la producción de metanol a partir de gas de síntesis, y la oxidación de SO2. Las ecuaciones son parametrizadas y resueltas de acuerdo a un método de continuación por homotopía. Los resultados demuestran que la concentración de inertes y la presión tienen marcada influencia sobre el tamaño de la región de inestabilidad.

PALABRAS CLAVE: Región de inestabilidad, Reactores de lecho empacado adiabáticos, Método de continuación por homotopía.


 

1. INTRODUCTION

In the design of adiabatic packed bed reactors, there are two types of temperature excursions that must be specially considered: (1) the existence of "hot spot"; temperatures, characteristic of reactors with no recycle of gas; and (2) the unstable "runaway"; temperature, characterized by large recycle of gas. Kimura and Levenspiel [1] schematically have shown that a conversion vs. temperature chart can become a convenient tool for identifying these unstable points. For their recognition, they considered that "hot spot"; and "runaway"; temperatures can be related to two ideal reactors models: plug flow tubular reactor (PFTR) and completely stirred tank reactor (CSTR), respectively. With this background, they have developed a criterion to find stable points and an envelope of operating conditions which must be avoided in operation of adiabatic packed bed reactors. From the CSTR model, they linked the molar and adiabatic energy balances in a conversion vs. temperature plot. For a given feed temperature and operating residence time (t or t´ = t/CAo) the adiabatic energy balance line and the S-shaped material balance intersect each other at either one or more stable or unstable points. At the unstable points it is possible to generalize as follows:

Equation (1) was claimed by Kimura and Levenspiel [1] as a fundamental criterion for mapping the instability operation region. Thus, to map the instability region, it is necessary to establish the points where the slope of the constant t' line (molar balance) equals to the slope of the corresponding adiabatic line (adiabatic energy balance). These points are joined in a t'-conversion-temperature chart to obtain the boundary of unstable region. Those kinds of charts are used in the literature to solve optimization problems: e.g., to calculate reactor sizes, and to compare alternative designs in reaction rate-conversion-temperature charts [2, 3]. To create those kinds of charts, besides simple rate laws (e.g., first order reactions), computational procedures with complex numerical methods are required due to nonlinearity rate laws and the resulting models. Continuation or homotopy methods have been demonstrated to be robust for the location of all real roots of numerous types of chemical engineering problems involving single or systems of nonlinear equations [4]. Homotopy methods obtain the solution of nonlinear equations by tracking a path from a starting point that is the solution of a much simpler function. Both, the original function and the simpler function are embedded into a parameterized initial-value problem in ordinary differential equations, solved by a predictor-corrector algorithm.

This paper presents how numerically, by using a homotopy continuation method, it is possible to find the instability regions of complex reaction rate systems, based on the pioneer schematic ideas of Kimura and Levenspiel [1]. Three different cases of special industrial interest and complexity are presented: ammonia, methanol, SO3 synthesis.

 

2. MATHEMATICAL MODEL OF UNSTABLE REGIONS IN PACKED BED REACTORS

For the CSTR, the adiabatic energy balance can be arranged as follows [5]:

where XA corresponds to the conversion of the limiting reactant. The CSTR molar balance is defined as:

The objective is to find the coordinated points of equation (4), (XA, T), for a given t', with a value of J that satisfy equation (3). It also means to solve the following equation:

It is difficult to evaluate the term d(-rA)/dT due to its nonlineal dependence on conversion and temperature. This implies to solve the following equation:

This is a general expression that corresponds to the instability criteria presented in equation (1). The complexity of applying and solving equation (7) depends on the rate law expression.

To accomplish a numerical solution of equation (1), a second function can be defined as follows:

Thus, equations (5) and (8) have to be solved using a homotopy method (details presented in the appendix). As a result, these two equations are grouped in a new variable named H as follows:

Equation (9) contains XA and T, which are grouped in a new variable Z. As XA = XA(T), orT = T(XA). the derivative of equation (9), with respect to Z is:

Basically, the algorithm consists of finding values of the new variable, Z, for a given t'. This can be made using a multivariable Newton-Raphson method.

 

3. REACTION SYSTEMS

Here the rate law expressions, for three cases considered in this work, are presented. The rate for ammonia synthesis is [2]:

The formation of methanol based on CO with H2 is [6]:

The oxidation of sulfur dioxide to sulfur trioxide, given by [2]:

 

4. RESULTS

The proposed algorithm proved to be robust and no convergence problems were detected for any of the analyzed rate laws. Figure 1 presents the r esults obtained for ammonia, methanol and sulfur trioxide synthesis. Calculation time depended on the complexity of the rate law.


Figure 1.
The instability region of the adiabatic packed bed reactors for ammonia, methanol, and SO3 synthesis. Included curves at t' constant and, for the ammonia case, also at rate constant

For the three instability regions of Figure 1 computation time were of 115,35; 101,45; and 49,2 seconds for ammonia, methanol, and SO3, respectively using a laptop with 1,8 Mhz processor and 776 Mb of RAM. In general, a region of instability is far from the equilibrium (-rA = 0). It is always to the left of the optimal progression of temperature (maximum reaction rate line in the reaction rate - conversion - temperature chart).

The size of the instability region highly depends on the inert concentration as presented for ammonia synthesis (Figure 2). However, for a high concentration of inerts the instability region disappears. Since the instability region has its origins in the difference between generated and removed energy of the system, it can be expected that a decrease in the reactant concentration diminishes the extension of the instability region.


Figure 2.
Effect of concentration of inerts on the size of the instability region of an adiabatic packed bed reactor for ammonia synthesis

The trajectories of the instability region have a maximum temperature (Figure 2). The maximum temperature represents the temperature limit for the instability region, and it could be correlated to the adiabatic operation line of a packed bed reactor. The curve of the locus of maxima is also shown in Figure 2. Above this locus curve the derivative dXA/dT is always negative thus, reactor conversion decreases as temperature increases. Below the locus curve reactor conversions always increases. Consequently, the minimum feed temperature can be found from the adiabatic reactor energy balance (equation (3)). The locus of maxima curve is controlled by operation parameters, such as: type and concentration of inerts (its heat capacity), total pressure, molar feed rate, heat transfer area, temperature of fluid inside the heat-exchange tube; and reactor construction material (its global heat-transfer coefficient) etc... Two of them are presented here: the effect of inert concentration in Figure 2 and the influence of the total pressure on the instability region that is presented in Figure 3, both for ammonia synthesis. In Figure 3, it is possible to observe that upon a decrease on total pressure (at 186 atm, light gray line), the locus curve moves downwards, and slightly to lower temperatures, comparing to those at higher pressures (e.g., 296 atm).


Figure 3.
Effect of total pressure on instability region size of an adiabatic packed bed reactor for ammonia synthesis

An increase of total pressure leads to an increase of the extent of the instability region. This can be expected since for an equilibrium reaction, where the number of moles of products is lower than that of reactants, high pressures shifts the equilibrium to the right. Thus, at higher pressures the reaction extent is higher as well as the heat evolution. The locus of maxima of those trajectories is adjusted to straight line, parallel to the temperature axis.

 

5. CONCLUSIONS

In this paper, a numerical approach, along the original schematic lines of Kimura and Levenspiel [1], for mapping the region of instability of adiabatic packed bed reactors has been shown for three reaction rate laws of high complexity and industrial interest. It was found that the instability region is far from the equilibrium (-rA = 0), and to the left of the optimal progression of temperatures. The shape and size of the instability region of an adiabatic packed bed reactor are controlled by type and concentration of inerts (its heat capacity), total pressure, molar feed rate, heat transfer area, temperature of fluid inside the heat-exchange tube; and reactor construction material. For a high concentration of inerts, the instability region disappears and the increase in total pressure leads to an increase of the extent of the instability region. The proposed homotopy continuation method appears to be a reasonable simple and accurate numerical technique, for which a simple personal computer amply suffices.

Nomenclature

CAO = Concentration of reactant in feed, mol/m3.
CPi = Specific heat of the fluid per mole of entering reactant, cal/(mol.K).
-DHRxn = Heat of reaction per mole of reactant A, cal/mol.
k i = Reaction rate constant.
Ki = Adsorption equilibrium constant.
Kp = Partial pressure equilibrium constant.
pi = Partial pressure of i.
R = Ideal gas law constant, 1.9859 cal/(mol.K).
-rA = Rate of reaction.
T = Temperature, K.
To = Entering temperature, K.
XA = Fraction of reactant A converted into product.
t = Space time based on unit mass or volume of catalyst.
qi = Ratio of the molar flow of species i entering to the molar flow of A entering

 

APPENDIX

Continuation or homotopy methods have been served as useful tools in modern mathematics and engineering [4]. Stated briefly, a homotopy method consists of the following. Suppose the situation in which very little of a priori knowledge concerning zero points of a function F(x) is available. A normal iteration method will often fail to determine the values that make zero the function F(x) because poor starting values are likely to be chosen. As a possible solution, a homotopy or deformation function H(x,l) is defined. It consists of a linear combination of two functions:

(A1)

(A2)

where G(x) is a function for which a zero is known, and l, the continuation or homotopy parameter, which allows tracking a solution path that connects an arbitrary starting point to the solution of F(x) = 0. A standard deformation which is often used is the global homotopy function:

(A3)

The choice of G(x) is arbitrary. As the parameter of continuation is in the range [1 0], a series of solutions to H(x,l) = 0 traces a path to the solution of F(x) = 0. H(x,l) is recursively solved at each value of l, using an appropriate numerical method (e.g., a newton type method).

In some instances, even if the function is parametrizable with respect to l, it may be necessary to choose an extremely small increment for the algorithm to succeed. To remedy this failure or poor performance it can be consider a differential arclength form of equation (A3). Thus, defining a new group of variables:

(A4)

In this way:

(A5)

Differentiating the equation (5) with respect to arclength , along the solution path:

(A6)

Defining:

(A7)

Where l(A) is the tangent vector induced by A. Now, numerical methods for solving initial value problems may be applied to

(A6). It should guaranty that the solution curve Z consists of zero (or near zero) points (ui) of H. Hence, one leads to integrate numerically the equation (A6) very coarsely and then locally use an iterative method for solving the equation (A5) as a stabilizer.

To describe how the points are generated along the curve, it is necessary to suppose that a point has been accepted so that, where is a numerical value describing the tolerance of the numerical method (). This accepted point can be presented as follows:

(A8)

where xo is a vector of length defined and with any numerical value.

Such a definition of the equation (A8) is possible since the following statement (A9) is always true:

(A9)

After accepting this point, it proceeds to integrate the following initial value problem:

(A10)

To obtain a new point along the curve , it is necessary to perform a predictor step due to the difficulty to obtain.

Normally, this predictor step is obtained by a simply numerical integration step of the initial value problem (equation A10). Thus, with an Euler predictor method:

(A11)

The point is just an approximation to the point . To obtain the correct point that satisfies , it is necessary to apply a corrector step, such as Newton-type method. This homotopy method applies until . This procedure is applied now in two steps for mapping the instability region for tubular packed bed reactors.

 

REFERENCES

[1] KIMURA, S.; LEVENSPIEL, O. Temperature Excursions in Adiabatic Packed Bed Reactors. Ind. Eng. Chem. Proc. Des. Dev., 16, 145 - 148, 1977.         [ Links ]
[2] FROMENT, G. B.; BISCHOFF, K. B. Chemical Reactor Analysis and Design. 2nd Ed. John Wiley & Sons. 1990.         [ Links ]
[3] LEVENSPIEL, O. Chemical Reaction Engineering. 3rd Ed. John Wiley & Sons. 1999.         [ Links ]
[4] GRITTON, K. S.; SEADER, J. D.; LIN, W-J. Global homotopy continuation procedures for seeking all roots of a nonlinear equation. Comp. Chem. Eng., 25, 1003 - 1019, 2001.         [ Links ]
[5] FOGLER, H. S. Elements of Chemical Reaction Engineering. 4th Ed. Prentice Hall. 2006.         [ Links ]
[6] BELFIORE, L. A. Transport Phenomena for Chemical Reactor Design. John Wiley & Sons. 2003.
        [ Links ]

Creative Commons License Todo el contenido de esta revista, excepto dónde está identificado, está bajo una Licencia Creative Commons