Héctor Armando Durán Peralta1 y Luis Fernando Córdoba C.2
1 M.Sc. Ingeniería Química. Profesor, Departamento de Ingeniería Química, Universidad Nacional de Colombia, Bogotá. haduranp@unal.edu.co
2 M.Sc. Ph.D., Ciencias Químicas. Profesor, Departamento de Ingeniería Química, Universidad Nacional de Colombia, Bogotá. lfcordobac@unal.edu.co
RESUMEN
]]> Abunda la literatura referente al análisis de estabilidad de reactores con parámetros globalizados de concentración y temperatura (por ejemplo el CSTR), en cambio es escasa la literatura sobre la estabilidad de reactores con parámetros distribuidos donde existe distribución espacial de concentración y temperatura, como es el caso del reactor tubular PFTR. Este documento analiza la estabilidad del reactor PFTR isotérmico y no isotérmico para una reacción con cinética de primer orden utilizando la funcional de Lyapunov. Se trabaja con una cinética de primer orden pues un objetivo de este artículo es mostrar cómo se aplica la funcional de Lyapunov al análisis de un reactor de parámetros distribuidos, dado que es casi inexistente la literatura sobre el método de la funcional de Lyapunov aplicada a la estabilidad de reactores (técnica usada en el análisis de estabilidad de sistemas en ingeniería eléctrica). El análisis de estabilidad dio como resultado perfiles de temperatura y concentración asintóticamente estables para los casos PFTR isotérmico, no isotérmico con constante cinética independiente de la temperatura y PFTR no isotérmico adiabático. Para el PFTR con retiro de calor el análisis condujo a una región de estabilidad asintótica y a una región incierta donde puede o no haber oscilaciones.Palabras clave: funcional de Lyapunov, PFTR, estabilidad en reactores
ABSTRACT
The stability of reactors having encompassing concentration and temperature parameters, such as continuous flow stirred tank reactors (CSTR), has been widely explored in the literature; however, there are few papers about the stability of tubular reactor having distributed spatial concentration and temperature parameters such as the plow flow tubular reactor (PFTR). This paper analyses the stability of isothermal and non-isothermal PFTR reactors using the Lyapunov functional method. The first order kinetic reaction was selected because one of this paper’s objectives was to apply Lyapunov functionals to stability analysis of distributed parameter reactors (technique used in electrical engineering systems’ stability analysis). The stability analysis revealed asymptotically stable temperature and concentration profiles for isothermal PFTR, non-isothermal PFTR with kinetic constant independent of temperature and adiabatic non-isothermal PFTR. Analysis revealed an asymptotically stability region for the heat exchange reactor and an uncertain region where it may have oscillations.
Keywords: Lyapunov functional, PFTR, stability reactor.
Recibido: septiembre 19 de 2006
Aceptado: marzo 1 de 2007
Introducción
En la solución de diferentes problemas de ingeniería a veces se necesita acudir a métodos de cálculo como lo son los variacionales, ello implica utilizar funcionales, las cuales intervienen en la solución de ecuaciones diferenciales ordinarias y ecuaciones diferenciales parciales (Burden y Douglas, 1998), estos métodos también son fundamentales en la solución de problemas de optimización (Dale y Watson, 1982). Además, son una herramienta básica en el estudio de la estabilidad de sistemas dinámicos con parámetros distribuidos; cuando se utiliza el método variacional para el estudio de la estabilidad resulta una funcional que se denomina la funcional de Lyapunov. Algunas situaciones de estudios de estabilidad en las cuales interviene la funcional de Lyapunov pueden consultarse, por ejemplo, en: (Kouachi Said, 2001; Kolmanovskii y Shaikhet, 2002; Hatvani, 2001; Shaikhet, 2004). Uno de los objetivos de este documento es mostrar cómo se utiliza la funcional de Lyapunov en el estudio de la estabilidad de un reactor tubular.
Normalmente cuando se opera un reactor en una planta química preferentemente se opera en estado estacionario (E.E). Y es importante que el ingeniero se pregunte si la operación del reactor es estable o inestable, es decir, si el reactor se ha diseñado para operar en un régimen estacionario para unas determinadas condiciones. Puede suceder que una pequeña perturbación en las variables de operación saque al reactor del E.E y a medida que pasa el tiempo el reactor no regrese a las condiciones que se desea. En este caso no se ha previsto si en el diseño el E.E es estable, lo cual puede repercutir negativamente, el desarrollo del proceso. Por ello es importante analizar la estabilidad de los reactores y en general la de los equipos que se proyectan para trabajar en estado estacionario. Este artículo analiza la estabilidad del reactor tubular para una reacción con cinética de primer orden utilizando el método de la funcional de Lyapunov; se trabaja este método dado que en la literatura existe abundante información sobre la estabilidad del reactor continuo de mezcla completa (CSTR) pero es muy escasa la literatura sobre el análisis del reactor tubular PFTR utilizando el método de la funcional de Lyapunov.
]]> En análisis de reactores es casi inexistente, por no decir nula, la información sobre la estabilidad de reactores para una reacción específica, ya que las industrias que trabajan en la síntesis de un determinado producto hace el análisis de estabilidad del reactor con la reacción que produce el compuesto de interés pero no divulgan los resultados. Por lo tanto, en la literatura de análisis de reactores suele encontrarse sistemas de reacción genéricos pero no datos para una reacción en particular. En este documento se analiza la estabilidad del estado estacionario del reactor PFTR tanto isotérmico como no isotérmico, para el sistema químico homogéneo no isotérmico
con cinética de primer orden (-rA=kCA) utilizando la funcional de Lyapunov. Estabilidad del reactor tubular isotérmico (PFTR)
La Ecuación del modelo que simula el PFTR isotérmico sin cambios de densidad de la mezcla es (Rutherford, 1975):
Condiciones: C(z,0)=C0; C(0,t)=C0 t≥0
Antes de iniciar el análisis de estabilidad se expresa el balance de materia (ecuación 1) y la condición inicial y de frontera en forma adimensional:
Condiciones: x(ζ,0)=1; x(0,τ)=0 si τ≥0
]]> Donde la concentración adimensional, temperatura adimensional, tiempo adimensional y longitud adimensional son respectivamente:
Y “P” es un número adimensional 
A continuación puede comenzar el análisis determinando el perfil de concentración de estado estacionario (ecuación 3) utilizando el balance de materia adimensional (ecuación 2),
Para analizar la estabilidad del perfil de estado estacionario, (ecuación 3), en sistemas de parámetros distribuidos, una forma es utilizar la funcional de Lyapunov. Es de recordar en la teoría de la estabilidad que la función de Lyapunov y la funcional de Lyapunov no son únicas, es decir, hay una infinidad de funciones y funcionales de Lyapunov, lo cual es una ventaja pues puede determinarse por ensayo y error. Aquí se procedió formalmente por un proceso al límite para construir la funcional de Lyapunov; este proceso condujo a la siguiente expresión:
A continuación se va a comprobar si la funcional (4) es funcional de Lyapunov. La funcional propuesta debe satisfacer tres condiciones (Zubov, 1957; Hahn,1963, y Wang, 1964):
]]> i- En una vecindad del perfil de E.E, la funcional debe ser positiva:
ii- La funcional ha de cumplir 
en el perfil de E.E y el perfil de E.E debe ser una extremal de la funcional 
. En este caso el perfil de E.E tiene que minimizar la funcional v.
iii- La derivada 
es una vecindad del perfil de E.E.
si se cumplen las tres condiciones anteriores la funcional es una funcional de Lyapunov y ello implica: a) que el perfil de concentración es asintóticamente estable si se cumple la desigualdad de la tercera condición, y b) es estable si se cumple la igualdad de la tercera condición.
A continuación va a verificarse si la funcional propuesta (ecuación 4) es funcional de Lyapunov.
Primera condición:
La funcional debe ser positiva en un entorno cerca de la solución de estado estacionario.
Para la funcional (ecuación 4) en el corchete el perfil transiente 
siempre es mayor que el perfil de E.E 
, esto hace que la perturbación [
] sea positiva pues el reactivo se está consumiendo y además la perturbación está al cuadrado, por lo cual la funcional es positiva en una vecindad de radio finito con centro en el perfil o solución de estado estacionario. Esta aseveración se justifica dado que la funcional es una integral y una propiedad de las integrales dice que al provenir la integral como el límite de una sumatoria de términos positivos (como es el caso de esta funcional) entonces la integral será positiva.
Segunda condición: a) Para ello se evalúa la funcional en el perfil de E.E.
Aquí F es:
Reemplazando F en la ecuación de Euler y efectuando las manipulaciones algebraicas del caso, se llega a:
Entonces
]]>
La cual es el perfil de E.E 
y por lo tanto satisface la condición de la extremal de Euler, cumpliéndose la segunda condición para la funcional de Lyapunov.
Tercera condición:
La derivada de la funcional de Lyapunov con respecto al tiempo es:
Esta funcional, mediante el balance de materia en el reactor, se transforma en:
Se llega a:
El primer corchete del integrando es >0 porque el perfil de concentración en el instante t es menor que el perfil de concentración en estado estacionario, y el segundo corchete es <0 pues este es el término de acumulación de masa, que al consumirse el reactivo el cambio, será negativo. Sabiendo que la suma de términos negativos es negativa, y que la integral proviene de una sumatoria de términos negativos, la integral será entonces negativa, luego en todo el dominio de concentraciones adimensionales.
Como satisface las tres condiciones, esta funcional es una funcional de Lyapunov y por lo tanto el perfil de E.E es asintóticamente estable en todo el rango de concentraciones. Un teorema dice que cuando existe la funcional de Lyapunov en una región “R” el sistema en esa región no tiene soluciones periódicas (Wang, 1964). Finalmente, para este caso del reactor tubular isotérmico el balance de materia es una ecuación lineal, por lo cual el perfil de E.E es único, y no puede haber multiplicidad de perfiles de E.E.
Análisis de la estabilidad del reactor tubular no isotérmico (PFTR). Caso de constante cinética independiente de la temperatura
El modelo para este reactor (Perlmutter, 1972) es:
]]>
Ecuaciones válidas en la región espacial R del reactor, donde R: 0 ≤ z < L.
Condiciones:
En el presente documento, para evitar confusiones en el signo del calor de reacción y no enredar el análisis que se hará posteriormente, se adoptará el siguiente convenio: se trabajará con el valor absoluto del calor de reacción Abs(ΔH)≡ΔH. Así, en la ecuación (6) el término ΔHkapC indica por el balance que la reacción es exotérmica.
A continuación se expresan los balances en forma adimensional:
Balance de materia
Balance de energíaEcuaciones válidas en la región espacial R del reactor, donde R: 0≤ζ<1. En los balances anteriores los números adimensionales y, P, S y H son:
]]>
Y las condiciones del problema pasan a:
Para analizar la estabilidad de este reactor se utilizará el método del plano compuesto de fases (CPP). En dicho método el sistema se globaliza utilizando coordenadas lagrangianas, lo cual implica utilizar el operador derivada material, que es equivalente a lo que en fenómenos de transporte se denomina derivada substancial (Bird, Stewart y Lightfoot 1982), 
. Expresando el balance de materia y energía en términos de la derivada material, los balances se transforman:
Como resultado el sistema original queda globalizado con el plano de fases, lo cual da como resultado el plano de fases compuesto (CPP). Al determinar la solución de estado estacionario del sistema anterior resulta:
]]>
xE=0; yE=ywHomogeneizando el sistema, es decir, trasladando las ecuaciones a la solución de estado estacionario (xE,xE):
Con los valores propios se determina la estabilidad de la solución de E.E. La ecuación característica que da los valores propios es:
Los valores propios son:
y 
dependen de parámetros que siempre son positivos; así, los valores propios nunca serán positivos ni imaginarios. La Ecuación para obtener las trayectorias en el plano compuesto CPP es:
La ecuación de las trayectorias de fase en el plano compuesto es:
Donde Θ0=y0-yw.
Para valores típicos reportados por Wang (1968) en un reactor tubular con un sistema que tiene cinética de primer orden, k0(1.2*1010h-1), energía de activación b = 9.400 ºR y coeficiente de transferencia de calor h = 7.000 Btu/ft3/°F /h. Lo cual permite evaluar kap:
Utilizando los valores de los parámetros anteriores se determinaron las trayectorias de fase en el plano compuesto CPP obtenidas de la ecuación (8). La figura 1 muestra estas trayectorias de fase.
Obsérvese en la Figura 1 que al ubicar la posición inicial (x0,θ0), cualquier posición inicial conduce a (x=0,θ=0), valor de estado estacionario, que en el plano CPP es un punto, el cual viene a corresponder en coordenadas estacionarias al perfil de estado estacionario.
El balance de energía incluye ahora el término de disipación de calor y la constante de la cinética varia con la temperatura. En estas condiciones los balances en forma adimensional son:
Haciendo extensivo el método de la funcional de Lyapunov que se aplicó al reactor tubular isotérmico, la funcional queda para este caso:
En este caso la extremal debe satisfacer:
La funcional propuesta para determinar la estabilidad es:
El siguiente paso es verificar si se satisfacen las tres condiciones para que la funcional sea funcional de Lyapunov. Por analogía con el caso isotérmico se hace el mismo tratamiento para este caso cumpliéndose las dos primeras condiciones para la funcional de Lyapunov. La tercera usualmente es la que presenta problema, para este caso la tercera condición a cumplir es:
Utilizando los balances de materia y energía, la anterior ecuación se transforma:
]]>
Para el caso adiabático no se retira calor, de modo que H=0
La primera integral es idéntica a la del caso isotérmico y sigue el mismo análisis que se hizo para ese caso, por lo tanto la primera integral es negativa en todo el dominio de concentración. En la segunda integral los perfiles transientes de temperatura son menores que el perfil de E.E por ser la reacción exotérmica y el reactor adiabático, por lo cual (y-yE) es <0, y para que el integrando de la segunda integral sea negativo es necesario exigir que
En la expresión anterior 
, ya que el gradiente adimensional de temperatura es positivo, luego la expresión 
es positiva, esto completa la tercera condición y por lo tanto la funcional propuesta es funcional de Lyapunov, y ello garantiza la estabilidad asintótica de los perfiles de concentración y temperatura cuando el reactor opera no isotérmicamente y adiabáticamente en el dominio de operación de concentración y temperatura del reactor tubular. Desde el punto de vista físico al ser la reacción irreversible llega un momento en el cual se agota el reactivo y por lo tanto deja de liberarse calor y la temperatura del sistema no aumenta más, estabilizándose, esto en el caso de un reactor de gran longitud, en el caso ideal de longitud infinita.
Para el caso no isotérmico con retiro de calor:
]]> Antes de continuar con la funcional se quiere identificar regiones donde pueda descartarse que existan oscilaciones en el tiempo en concentración y temperatura. Para ello se trabaja con coordenadas móviles siguiendo el movimiento. Las ecuaciones de balance de materia y energía en términos de la derivada material son:
El estado estacionario en el plano compuesto se obtiene de resolver las ecuaciones (9) y (10):
Siendo la solución: x=0 y=yw
Entonces el estado estacionario en el plano de fase compuesto es: (x,y)=(0,yw)
]]> Para obtener las trayectorias de fase se usa la expresión:
El campo vectorial asociado según la ecuación (11) es:
Para poder descartar trayectorias cerradas en alguna región del plano compuesto CPP se utiliza la divergencia del campo vectorial asociado:
El teorema de Bendixon dice que si en una región del espacio de fases la divergencia del campo vectorial asociado resulta ser de signo constante, entonces en esa región no pueden haber trayectorias cerradas y por lo tanto el sistema no oscila en esa región; en nuestro caso, si la divergencia en una región R resulta ser de signo constante, allí no oscilarán ni la concentración ni la temperatura. A continuación se desea establecer la región donde exista exclusión de oscilaciones.
]]> Para que la divergencia sea de signo constante (por ejemplo negativa), en la ecuación (13) debe cumplirse:
Utilizando las expresiones correspondientes a los números adimensionales de la expresión (14), esta queda:
La desigualdad (15) permitió construir la Figura 2 en términos de parámetros adimensionales y su equivalente (Figura 3) en términos de temperatura y concentración de reactivo, donde se muestra la región de exclusión de trayectorias cerradas (oscilaciones), es decir, se obtuvo un RAS (región de estabilidad asintótica). La condición (15) en términos de concentración y temperatura es:
Con los valores de los parámetros reportados por Wang (1968), se obtiene:
]]>
Tomemos una concentración inicial de C0=0.035 Lbmol/ft3 con T0=600 ºR. Ubicando el punto (0.035, 600) en el gráfico anterior el punto queda situado en la región fuera de exclusión de oscilaciones (región no sombreada), donde no se puede asegurar si hay o no oscilaciones, pero si el reactor comienza a operar desde el punto B (C=0.019) y como la concentración de reactivo C al transcurrir el tiempo disminuye a partir de 0.019 y la reacción es exotérmica, la reacción evoluciona adentrándose más en la región sombreada, en la cual se puede asegurar que no se presentan oscilaciones de temperatura ni de concentración.
En la región del plano χy donde se cumpla la desigualdad x<χy2 para la ecuación (13) se excluye la existencia de trayectorias cerradas en el campo vectorial en el plano de fase compuesto (CPP). En la región anterior donde no hay trayectorias cerradas, y si esta incluye los puntos estacionarios, hay un teorema que dice que esta es una condición de existencia de la funcional de Lyapunov.
Retomando la funcional
como no hay posibilidad de oscilaciones en la región donde se cumple la desigualdad (15), y con retiro de calor fijo siempre que se produzca calor debido a la reacción, entonces la segunda integral de la funcional es negativa, y como la primera integral es negativa tenemos que la tercera condición se cumple, la funcional propuesta es funcional de Lyapunov, y por lo tanto si se opera el reactor con los parámetros de modo que se satisfaga la desigualdad (14) el sistema resultará con perfiles de temperatura y concentración asintóticamente estables.
Conclusiones
El presente documento analizó la estabilidad del reactor tubular tanto isotérmico como no isotérmico utilizando la funcional de Lyapunov. En el caso isotérmico, el método de la funcional de Lyapunov permitió predecir que el PFTR isotérmico tiene perfil de concentración único de estado estacionario y es asintóticamente estable en todo el dominio de concentración si el reactor es semiinfinito.
]]> Para el caso de reactor PFTR no isotérmico con constante cinética independiente de la temperatura se analizó la estabilidad de los perfiles de temperatura y concentración de estado estacionario utilizando el método del plano compuesto CPP y empleando valores propios, al utilizar el plano compuesto CPP fue necesario utilizar coordenadas siguiendo el movimiento, observándose que los perfiles de concentración y temperatura resultaron asintóticamente estables.Para el caso de reactor PFTR no isotérmico adiabático con constante cinética variable con la temperatura se utilizó la funcional de Lyapunov con la finalidad de analizar la estabilidad de los perfiles de estado estacionario, determinándose que existe la funcional de Lyapunov, y el reactor en este caso tiene perfiles de concentración y temperatura de estado estacionario asintóticamente estables en el dominio de operación de concentración y temperatura, ello si el reactor es semiinfinito.
Para el caso de reactor PFTR no isotérmico con retiro de calor y constante cinética variable con la temperatura, también se utilizó el método del plano de fase compuesto CPP para descartar oscilaciones en temperatura y concentración, como resultado se obtuvo que el sistema no oscila en la región sombreada mostrada en las figuras 2 y 3 y por lo tanto en esta región existe la funcional de Lyapunov y los perfiles de temperatura y concentración son asintóticamente estables. La región no sombreada es una región incierta pues no se sabe si el sistema oscila o no. En un futuro trabajo se podría proponer una funcional de Lyapunov que amplíe la región de estabilidad asintótica (RAS) encontrada en el presente artículo.
Nomenclatura
C : Concentración (Lbmol/pie3)
t : Tiempo (h)
z : Coordenada axial (pie)
x : Concentración adimensional
y : Temperatura adimensional
τ : Tiempo adimensional
]]> ζ : Longitud adimensionalL : Longitud del reactor tubular (pie)
u : Velocidad lineal de la corriente fluida en el reactor (pie/h)
T : Temperatura local dentro del reactor tubular (°R)
TL : Temperatura en la corriente de salida (°R)
T0 : Temperatura en la corriente de entrada (°R)
C0 : Concentración en la entrada al reactor
P : Número adimensional
kap : Constante de velocidad de reacción (h-1)
xE : Concentración adimensional de estado estacionario
]]> yE : Temperatura adimensional de estado estacionario : Funcional de Lyapunov
F : Función en la condición de Euler para que exista la Extremal de la funcional
cp : Calor específico (Btu/pie3/°R )
ΔH : Calor de reacción (Btu /mol)
h : Coeficiente de transferencia de calor convectivo (Btu/pie3/°R/h )
Tw : Temperatura del fluido en la longitud infinita del reactor
S : Número adimensional
H : Número adimensional
λ : Valor propio
]]> yw : Temperatura adimensional del fluido en la longitud infinita del reactorθ : Temperatura adimensional trasladada
k0 : Factor de frecuencia en el modelo cinético de la reacción (h-1)
Bibliografía
Bird, R. B.: Stewart, E. and Ligthfood, E., Fenómenos de transporte, Editorial Reverte, 1982. [ Links ]
Burden Richard, L. and Douglas Faires, J., Numerical analysis., 6ª edn, Brooks Cole Publishing, 1998. [ Links ]
Dale, F. Rudd y Watson Charles, C. Estrategia en Ingeniería de Procesos., Ed. Alambra, Madrid, 1982.
Hahn, W., Theory and Application of Lyapunov’s Direct Methods., Englewood Cliffs, N, prentice-hall, 1963. [ Links ]
Hatvani, L., On the Asymptotic Stability by Lyapunov Functionals with Semidefinite Derivatives., Nonlinear Analysis, 47, 2001, pp. 4333-4343. [ Links ]
Kolmanovskii, V. and Shaikhet, L., Some Peculiarities of the General Method of Lyapunov Functionals Construction., Applied Mathematics Letters, Vol 15, 2002, pp. 355-360. [ Links ]
Kouachi Said., Existence of Global Solutions to Reaction-diffusion Systems Via a Lyapunov Functional., Electronic Journal of Diferential Equations, Vol. 21, No. 68, 2001, pp. 1-10. [ Links ]
Perlmutter Ben., Stability of Chemical Reactors., 1ª edn, New Jersey, Prentice-Hall, 1972. [ Links ]
Rutherford Aris., Questions of Uniqueness Stability and Transient Behavior., Clarendon Press Oxford, Great Britain, 1975. [ Links ]
Shaikhet, L., About Lyapunov Functionals Construction for Difference Equations with Continuous Time., Mathematics and Computers in Simulation, 66, 2004, pp. 509–521. [ Links ]
Wang, P. K. C., Advances in Control Systems.,. N, Y: Academic Press, 1964. [ Links ]
Zubov, V. I., Methods of A.M. Lyapunov and their Application., Leningrad, 1957. [ Links ] ]]>