Transcript
9. ESTABILIDAD DE SISTEMAS Para el análisis del comportamiento dinámico de sistemas, especialmente en el caso de los no lineales de difícil solución matemática, la naturaleza estable o inestable de su respuesta se puede determinar mediante algunos procedimientos como la construcción de un Diagrama de Fase. En este capítulo se plantea el modelo de un reactor de mezcla completa adiabático cuya estabilidad se analiza por los métodos numéricos empleados en los casos anteriores y mediante la construcción de su diagrama de fase. Se incluye, además, el análisis del comportamiento del reactor en estado estacionario.
9.1 DIAGRAMAS DE FASE Para un sistema con dos variables de estado, el diagrama de fase es el conjunto de curvas que representan a una variable de estado en función de la otra, cada una de ellas para una condición inicial diferente. El diagrama de fase permite determinar para un conjunto diferente de condiciones iniciales si el sistema es estable o inestable, es decir, si converge a unas condiciones finales en estado estacionario.
Estabilidad de sistemas con dos variables de estado El modelo matemático de un sistema con dos variables de estado y cuyas ecuaciones diferenciales no incluyen la variable tiempo en sus miembros derechos se puede representar de la siguiente manera:
x&1 (t ) = f1 ( x1 , x2 ) x& 2 (t ) = f 2 ( x1 , x2 )
(9.1a) (9.1b)
Se entiende por Punto de Equilibrio, todos aquellos puntos del sistema para los cuales
x&1 (t ) = x& 2 = 0
(9.2)
Las dos funciones f1 ( x1 , x2 ) = 0 y f 2 ( x1 , x2 ) = 0 se llaman Líneas de Fase y dividen el plano cartesiano en cuatro partes
148
Una trayectoria del sistema de ecuaciones (9.1a) y (9.1b) es el conjunto de puntos
{(x (t), x (t)) ∈ R 1
2
2
}
| x1 (t ) ∧ x2 (t ) son solución del sistema de ecuaciones para t ∈[to , t1 ]
Una trayectoria se llama estable si, para t → ∞ , x1 (t ) y x 2 (t ) convergen en el punto de equilibrio
Estabilidad de sistemas lineales Para un sistema lineal con variaciones nulas en sus variables de entrada el modelo en la forma del espacio de los estados se reduce a la ecuación matricial X& = AX
(9.1)
Las características de estabilidad de un sistema lineal se pueden analizar mediante la determinación de los valores propios de la matriz A, que se obtienen resolviendo la ecuación que resulta al igualar a cero el determinante de la matriz λI − A , es decir que:
det(λI − A) = 0
(9.2)
Siendo λ , los valores propios de la matriz A. Para un sistema con dos variables de salida, se tiene que
λ − A11 λI − A = − A21
− A12 λ − A22
det(λI − A) = λ2 − ( A11 + A22 )λ + A11 A22 − A12 A21 Los valores propios son las raíces de la ecuación de segundo grado
Mach
(9.2a)
(9.2b)
149
ó
λ2 − ( A11 + A22 )λ + A11 A22 − A12 A21 = 0
(9.2c)
λ2 − (Traza A)λ + det( A) = 0
(9.2d)
Si los valores propios de la matriz A son números reales el diagrama de fase es nodal y si dichos valores son números complejos conjugados el diagrama de fase es en espiral. En un diagrama de fase nodal las trayectorias de fase se acercan o se alejan al punto de equilibrio tangencialmente a las líneas de fase y este se observa como un receptor (nodo estable) o una fuente (nodo inestable) de trayectorias de fase. Si las trayectorias de fase se acercan y se alejan sin converger en el punto de equilibrio (nodo inestable), este se considera como un punto de silla. En un diagrama de fase en espiral las trayectorias de fase lo hacen describiendo curvas en espiral. Un resumen que muestre la relación entre la naturaleza de los valores propios de la matriz A y el comportamiento del sistema es el siguiente: 1. Si λ1 y λ 2 son números reales Si λ1 < 0 y λ2 < 0 Si λ1 < 0 y λ2 > 0 Si λ1 > 0 y λ 2 > 0
⇒ ⇒ ⇒
⇒
Nodo Estable (Receptor) Nodo Inestable (Silla) Nodo Inestable (Fuente)
2. Si λ1 y λ2 son complejas conjugadas ⇒ Si Re(λ1 ) < 0 Si Re(λ1 ) > 0
⇒ ⇒
Comportamiento Nodal
Comportamiento en Espiral
Nodo Estable (Receptor) Nodo Inestable (Fuente)
9.2 REACTOR DE MEZCLA COMPLETA ADIABATICO Las reacciones químicas son exotérmicas o endotérmicas y, por lo tanto, requieren que esa energía se remueva o se añada al reactor para mantener una temperatura constante. Las reacciones exotérmicas son las más interesantes para el estudio debido a los posibles problemas de seguridad causados por sus aumentos rápidos de temperatura y la posibilidad de un comportamiento singular tales como múltiples estados estacionarios (para un mismo valor de la variable de entrada pueden obtenerse varios posibles valores para la variable de salida)
Mach
150
En este módulo se considera un reactor de mezcla completa como el que se muestra en la figura 9.1 y se estudia una reacción exotérmica irreversible de una cinética de K1
primer orden del tipo A → B y se analiza un comportamiento muy interesante que puede presentarse con esta sencilla reacción
F, CAf(t)
CA T FJ, Tj(t)
FJ, Tji(t) F, CA(t)
Figura 9.1 Reactor de Mezcla Completa Adiabático En la figura 9.1 se ve que se alimenta y se remueve una corriente de un fluido en forma continua a través del reactor. Como el reactor es de mezcla completa, se considera que la corriente de salida presenta la misma concentración y temperatura que la masa reaccionante dentro del reactor. En forma similar, a través de la camisa que cubre al reactor se alimenta y descarga una corriente fluida. Se asume que el fluido que circula por la camisa se encuentra perfectamente mezclado y a una temperatura menor que la del reactor. Entonces, la energía se transfiere desde el reactor hasta la camisa a través de sus paredes, removiéndose el calor generado en la reacción.
Modelo No Lineal – Dominio Tiempo Por simplificación, se asume que la temperatura de enfriamiento en la camisa puede manipularse directamente de tal manera que no se requiera el balance de energía en la camisa. Además, se incluyen las siguientes consideraciones: mezclado perfecto, volumen constante y parámetros constantes. El modelo se plantea con un balance global de materia, un balance del componente A y un balance de energía a través del reactor así: Un balance global de materia a través del reactor:
Mach
151
d (Vρ ) = Fin ρ in − Fout ρ dt
(9.3)
Si se asume que la densidad es constante, entonces: Fin = Fout = F
(9.4)
Un balance de materia del componente A en el reactor:
V
d (C A ) = FC Af − FC A − rV dt
(9.5)
Siendo r, la velocidad de reacción por unidad de volumen Un balance de energía a través del reactor:
VρC p
dT = FρC p (T f − T ) + (− ∆H )rV − UA(T − T j ) dt
(9.6)
Siendo (−∆H )rV la rapidez de energía liberada en la reacción El modelo dinámico se resume a las ecuaciones (9.5) y (9.6) y, además, la ecuación de velocidad de reacción por unidad de volumen las que se resumen así:
Mach
dC A F = f 1 (C A , T ) = (C Af − C A ) − r dt V
(9.7)
dT F UA − ∆H = f 2 (C A , T ) = (T f − T ) + ( )r − (T − T j ) ρC p dt V VρC p
(9.8)
E r = k o exp − C A RT
(9.9)
152
Las ecuaciones (9.7), (9.8) y (9.9) son no lineales porque incluyen términos no lineales como el producto del flujo por la concentración y el exponencial de la temperatura.
Análisis en estado estacionario Las ecuaciones (9.7) a (9.9) se pueden combinar y reducir el sistema a dos ecuaciones algebraicas no lineales en función de la concentración de A y la temperatura en el reactor, para calcular estas variables en condiciones de estado estacionario así:
f 1 (C A , T ) = 0 =
E F C As (C Af − C As ) − k o exp − V RTs
(9.10)
f 2 (C A , T ) = 0 =
E F − ∆H UA C As − (T f − Ts ) + ( )k o exp − (Ts − T j ) ρC p V VρC p RTs
(9.11)
La solución analítica del sistema de ecuaciones (9.10) y (9.11) es complejo debido a su no linealidad. Para el cálculo de la concentración de A y la temperatura en el reactor en estado estacionario, el sistema de ecuaciones se define en el archivo reactor1.m y se aplic el comando fsolve de Matlab al archivo solreactor1.m con los siguientes parámetros: F/V, h-1 ko, h-1 (− ∆H ) , kcal/kgmol E, kcal/kgmol ρC p , kcal/(m3-ºC)
1 9703*3600 5960 11843 500
Tf, ºC CAf, kgmol/m3 UA/V, kcal/(m3-ºC-h) Tj, ºC R, kcal/(kgmol-K)
25 10 150 25 1.987
Al escoger valores iniciales para la solución numérica de un sistema de ecuaciones es importante tener una visión del posible rango de soluciones. Por ejemplo, como la concentración de A es 10 kgmol/m3 y la reacción consume A, el rango posible para la concentración de A es 0 < CA < 10. También es fácil mostrar un límite inferior para la temperatura de 298 K, que ocurriría si no hubiera reacción, porque las temperaturas del alimento y la camisa son 298 K. Se observa, también, que deberían correlacionarse la concentración de A con la temperatura. Si la concentración de A
Mach
153
es alta, significa que no ha ocurrido mucha reacción de manera que se ha liberado muy poca energía por la reacción y por lo tanto la temperatura no será muy diferente con respecto a la de la camisa y el reactor. Se observa en la Tabla 9.1, que las condiciones iniciales que se fijen para la solución del sistema muestran diferentes resultados para la concentración de A y la temperatura en estado estacionario. El caso número uno es de alta concentración y baja temperatura, el número dos de concentración y temperatura intermedia y el caso tres de baja concentración y alta temperatura Tabla 9.1 Condiciones iniciales y Soluciones del Sistema de Ecuaciones
NUMERO 1 2 3
CONDICIONES INICALES CA, kgmol/m3 9 5 1
T, K 300 350 450
SOLUCION CA, kgmol/m3 8.564 5.518 2.359
T, K 311.2 339.1 368.1
Análisis en estado dinámico – Modelo No Lineal – Dominio Tiempo La solución del modelo no lineal en el dominio del tiempo se presenta a continuación para los tres casos ensayados según lo muestra la tabla 1 y posteriormente se muestra el análisis de la respuesta mediante el diagrama de fase correspondiente.
Ensayo 1 (Alta Concentración – Baja Temperatura) La figura (9.2) muestra la respuesta dinámica del sistema para la concentración de A y la temperatura en el reactor para el ensayo 1 de la Tabla 9.1. En este, las condiciones iniciales son muy próximas al estado estacionario de baja temperatura y las gráficas muestran que las variables de estado convergen al estado estacionario de baja temperatura, es decir, 8.564 kgmol/m3 y 311.2 K
Ensayo 2 (Concentración y Temperatura Intermedia) Las figuras (9.3) y (9.4) muestran las respuestas dinámicas del sistema para la concentración de A y la temperatura en el reactor para el ensayo 2 de la Tabla 9.1. Si Mach
154
las condiciones iniciales son muy próximas al estado estacionario de temperatura intermedia. Se observan en dichas gráficas que cuando las condiciones iniciales son 5 kgmol/m3 y 325 K el sistema converge en el estado de baja temperatura y alta concentración, es decir, en 8.564 kgmol/m3 y 311.2 K, mientras que cuando las condiciones iniciales son de 5 kgmol/m3 y 350 K el sistema converge en el estado de alta temperatura y baja concentración, es decir, en 2.359 kgmol/m3 y 368.1 K Si se desarrollan muchas simulaciones con condiciones iniciales próximas al estado estacionario de temperatura intermedia, se encuentra que la temperatura siempre converge o al estado estacionario de baja temperatura o al de alta, pero no al intermedio. Esto indica que el estado intermediario de temperatura es inestable
Figura 9.2 Respuesta del Reactor Adiabático (Caso 1)
Ensayo 3 ( Baja Concentración y Alta Temperatura) La figura (9.5) muestra la respuesta dinámica del sistema para la concentración de A y la temperatura en el reactor para el ensayo 3 de la Tabla 9.1. En este, las condiciones iniciales son muy próximas al estado estacionario de alta temperatura y las gráficas muestran que las variables de estado convergen al estado estacionario de alta temperatura, es decir, 2.359 kgmol/m3 y 368.1 K
Mach
155
Figura 9.3 Respuesta de un Reactor Adiabático (Caso 2)
Figura 9.4 Respuesta de un Reactor Adiabático (Caso 2)
Mach
156
Figura 9.5 Respuesta del Reactor Adiabático (Caso 3)
Análisis en estado dinámico – Diagrama de Fase La Figura (9.6) muestra el diagrama de fase generado con el archivo solreactor2.m, que a su vez resuelve el sistema de ecuaciones diferenciales codificado con Matlab en el archivo reactor2.m. Se muestran los tres estados estacionarios; dos estables (alta y baja temperatura) señalados con pequeños círculos y uno inestable (temperatura intermedia), señalado con una pequeña cruz. Se observa que a condiciones iniciales de bajas concentraciones (0.5 kgmol/m3) y temperaturas relativamente bajas o intermedias (300 – 365 K) todas las trayectorias de fase convergen al estado estacionario de baja temperatura. Cuando la temperatura inicial se aumenta por encima de 365 K, las trayectorias de fase convergen al estado estacionario de alta temperatura. Por otra parte, a condiciones iniciales con alta concentración (10 kgmol/m3) y baja temperatura (300 – 325 K), las trayectorias de fase convergen al estado estacionario de baja temperatura. Cuando la temperatura inicial se aumenta por encima de 325 K, las trayectorias de fase convergen al estado estacionario de alta temperatura. También se observa que cuando la temperatura inicial se aumenta, aproximadamente, a valores mayores que 340 K, ocurre un sobresalto de temperatura a mas de 425 K antes de que el sistema converja al estado estacionario de alta temperatura. No se observa en la Figura (9.6), que a temperaturas iniciales Mach
157
mayores pueden ocurrir sobresaltos a temperaturas mayores que 500 K antes de la convergencia al estado estacionario de alta temperatura. Esto podría ser la causa de potenciales problemas de seguridad si, por ejemplo, ocurren algunas reacciones de descomposición a temperaturas altas. El análisis del comportamiento de un sistema mediante un diagrama de fase nos permite señalar condiciones iniciales problemas.
Figura 9.6 Diagrama de fase del Reactor Adiabático Se observa, además, que para ninguna condición inicial las trayectorias de fase no convergen en el estado estacionario de temperatura intermedia debido a que es un estado inestable. En la siguiente sección se determinan los valores propios de la matriz A del modelo linearizado y se verificarán con la naturaleza de sus valores el carácter estable de las condiciones estacionarias de baja y alta temperatura y el nodo de silla que corresponde a las condiciones de temperatura intermedia. Debería observarse que puede utilizarse un control por retroalimentación para operar el reactor en unas condiciones estacionarias de temperatura intermedia. El controlador mediría la temperatura en el reactor y manipularía la temperatura en el fluído de enfriamiento que se mueve por la camisa (o el flujo) para mantener el estado estacionario en la temperatura intermedia. También, podría utilizarse un controlador por retroalimentación para asegurarse que no ocurra un gran sobresalto de temperatura a partir de ciertas condiciones iniciales. Mach
158
Modelo Lineal – Dominio Tiempo – Espacio de los Estados Los miembros derechos de las ecuaciones (9.7) y (9.8) son funciones no lineales de las variables concentración de A en el reactor y en el alimento, la temperatura del alimento, en el reactor y en la camisa. Por lo tanto, se pueden expresar así: dC A = A11C A + A12T + B11C Af + B12T f + B13T j dt
(9.12)
dT = A21C A + A22T + B21C Af + B22T f + B23T j dt
(9.13)
Siendo,
A11 =
∂f 1 F = − − ks ; ∂C A V
A12 =
∂f 1 = −C As k s' ; ∂T
A21 =
∂f 2 ( − ∆H ) = ks ; ∂C A ρC p
A22 =
∂f 2 (− ∆H ) F UA =− − + C As k s' ∂C B ρC p V VρC p
B11 =
∂f 1 F = ; ∂C Af V
B12 =
∂f 1 = 0; ∂T f
B13 =
∂f 1 = 0; ∂T j
B21 =
∂f 2 = 0; ∂C Af
B22 =
∂f 2 F = ; ∂T f V
B23 =
∂f 2 UA ; = ∂T j VρC p
E k s = k o exp − RT s k s' =
E ∂k s = k o 2 ∂T RTs
F − V − ks A= (−∆H ) ρC k s p
Mach
E E exp − = k s RT 2 RT s s
F UA (−∆H ) ' − − + C As k C V VρC p ρC p − C As k s'
F V B= 0
0
F V
UA VρC P 0
159 0 0 0 D = 0 0 0
1 0 C= 0 1
Análisis de la estabilidad del Reactor en el Espacio de los Estados Para el análisis de estabilidad en el espacio de los estados se calculan los valores propios de la matriz A corriendo el archivo reactor3.m para cada uno de los casos ensayados anteriormente así:
Caso 1 Para una condición inicial en estado estacionario con una concentración de 8.564 kgmol/m3 y una temperatura de 38.2 ºC, la matriz A y los valores propios son: A = -1.1680 -0.0886 2.0030 -0.2443 Lambda = -0.8957 -0.5166 Como ambos valores propios son negativos, la condición estacionaria de baja temperatura es estable, lo que está de acuerdo con los resultados anteriores
Caso 2 Para una condición inicial en estado estacionario con una concentración de 5.518 kgmol/m3 y una temperatura de 66.1 ºC, la matriz A y los valores propios son: A = -1.8124 -0.2324 9.6837 1.4697 Mach
160
lambda = -0.8369 0.4942 Como un valor propio es negativo y el otro es positivo, entonces se verifica que la condición estacionaria de temperatura intermedia es inestable
Caso 3 Para una condición inicial en estado estacionario con una concentración de 2.359 kgmol/m3 y una temperatura de 95.1 ºC, la matriz A y los valores propios son: A = -4.2445 -0.3367 38.6748 2.7132 Lambda = -0.7657 + 0.9584i -0.7657 - 0.9584i Como los valores propios son números complejos conjugados con parte real negativa, se verifica que la condición estacionaria de alta temperatura es estable.
9.3 COMPORTAMIENTO MULTIPLE EN ESTADO ESTACIONARIO El objetivo de esta sección es determinar cómo pueden surgir los múltiples estados estacionarios en un modelo no lineal como el que se estudia en esta lección. También se muestra cómo se generan la curvas de entradas y salidas en estado estacionario que muestren, por ejemplo, cómo varía la temperatura en el reactor en estado estacionario con la variación de la temperatura en la camisa en estado estacionario.
Mach
161
Curvas de Calor removido y de Calor generado La ecuación (9.10) permite deducir una ecuación para calcular la concentración de A en el reactor en estado estacionario en función de la temperatura del reactor que es:
C As
F C Afs V = E F + k o exp − V RT s
(9.14)
La ecuación (9.11) puede arreglarse en la siguiente forma: − ∆H F UA k o exp − E C As (Ts − T fs ) + (Ts − T js ) = ρC V VρC p p RTs
(9.15)
Si ambos miembros se multiplican por VρC p , se halla que: E C As FρC p (Ts − T fs ) + UA(Ts − T js ) = −∆HVk o exp − RTs
(9.16)
El miembro de la izquierda expresa la energía total removida en el reactor mediante el flujo y el intercambio de calor y el miembro derecho expresa la energía generada en la reacción. Se puede escribir como ecuación para el calor total removido la siguiente: Qremovido = [−UAT js − FρC p T fs ] + [UA + FρC p ]Ts
(9.17)
La ecuación (9.17) es lineal con respecto a la temperatura en el reactor en estado estacionario, siendo la pendiente de la línea la suma [UA + FρC p ] y el intercepto con el eje de las ordenadas la suma [−UAT js − FρC p T fs ] . Los cambios que se produzcan en la temperatura de la camisa o del alimento desviarán el intersecto pero no a la pendiente. Los cambios que ocurran en UA o F afectarán tanto a la pendiente como al intersecto. Mach
162
Al considerar el término que expresa el calor liberado podemos escribir la siguiente ecuación: E C As Q generado = (−∆H )Vk o exp − RTs
(9.18)
De tal manera que combinando la ecuación (9.18) con la (9.14) se obtiene la siguiente ecuación que al graficarla muestra una curva en forma de “S” para el calor generado en función de la temperatura del reactor.
Q generado
E C As k o exp − RTs = (−∆H )V E F + k o exp − V RTs
(9.19)
A partir de la ecuación (9.16) se ve que una solución en estado estacionario existe cuando se tiene una intersección entre las curvas de calor generado y calor removido
Efecto de los Parámetros de Diseño Si la pendiente de la curva de calor removido es mayor que la máxima pendiente de la curva de generación de calor, solamente es posible una intersección. Cuando se modifica la temperatura de la camisa o del alimento, las líneas de calor removido se desvían a la izquierda o a la derecha, de tal manera que la intersección puede ser a alta o baja temperatura dependiendo del valor de la temperatura de la camisa o alimentación. Mientras la pendiente de la curva de calor removido es menor que la máxima pendiente de la curva de calor generado, siempre habrá la posibilidad de tres intersecciones con el ajuste adecuado de la temperatura de la camisa o del alimento. (Intersecto). Si se modifica la temperatura de la camisa o del alimento, la línea de calor removido se desvía a la derecha o a la izquierda, donde solamente ocurrirá una intersección (o a baja o a alta temperatura)
Mach
163
Análisis de la estabilidad en el Reactor En la Figura 9.7 se superponen varias gráficas lineales posibles de calor removido con la curva de calor generado en forma de S obtenidas con el archivo reactor4.m codificado con Matlab
Figura 9.7 Calor Removido o Generado versus Temperatura en el Reactor La línea A intersecta la curva de calor generado a baja temperatura; la línea B la intersecta a una baja temperatura y es tangente a una alta temperatura; la línea C intersecta a baja, intermedia y alta temperatura; la línea D es tangente a una baja temperatura e intersecta a una alta temperatura; y la línea E tiene solamente una intersección a alta temperatura. Las líneas A, B, C, D y E se basan todas en los mismos parámetros del sistema con excepción de la temperatura de alimentación que se aumenta al desplazarse desde la línea A hasta la E, con valores respectivamente de 0, 5, 15, 21 y 30 ºC y el flujo de alimento con un valor de 60 m3/h. Este aumento ocasiona cambios en el intersecto sin modificar la pendiente de la línea de calor removido. Se observa en la Figura 9.7 que el balance de energía en estado estacionario se cumple en las condiciones de operación 3, 5 y 7, es decir, hay tres estados
Mach
164
estacionarios. Mediante un simple razonamiento físico se puede analizar cada uno de los tres estados estacionarios de la siguiente manera: Para la operación estacionaria a baja temperatura (Punto 3), si se perturban las condiciones de operación a una temperatura mas fría, es decir, T3 − δT , se genera mas calor que el que se remueve en el reactor y, por lo tanto, la temperatura comienza a aumentar hacia el valor de la temperatura en el punto 3, es decir, T3 . Si la perturbación es a una condición mas caliente, es decir, T3 + δT , se remueve mas calor que el se genera en el reactor y, por lo tanto, la temperatura comienza a disminuir hacia la temperatura en el punto 3. En resumen, la intersección de baja temperatura T3 es una condición de operación estable Para la operación estacionaria a temperatura intermedia (Punto 5), si se perturban las condiciones de operación a una temperatura mas fría, es decir, T5 − δT , se genera menos calor que el que se remueve en el reactor y, por lo tanto, la temperatura comienza a disminuir hacia el valor de la temperatura en el punto 3, es decir, T3 . Si la perturbación es a una condición mas caliente, es decir, T5 + δT , se genera mas calor que el se genera en el reactor y, por lo tanto, la temperatura comienza a aumentar hacia la temperatura en el punto 7, es decir, T7 . En resumen, la intersección de temperatura intermedia T5 es una condición de operación inestable Para la operación estacionaria a alta temperatura (Punto 7), el razonamiento físico es similar al de la operación estacionaria a baja temperatura. Es decir, la intersección de temperatura intermedia T7 es una condición de operación estable
Curvas de entradas y salidas en estado estacionario Se puede utilizar la Figura 9.7 para construir el diagrama de entradas y salidas en estado estacionario que se muestra en la Figura 9.8, siendo la variable de entrada la temperatura del alimento y la variable de salida la temperatura en el reactor. Se observa que la Figura 9.8 exhibe un comportamiento de histéresis. El término histéresis se utiliza para indicar que el comportamiento es diferente dependiendo del sentido en que las variables de entrada cambian. Por ejemplo, si se comienza con una temperatura del alimento baja, el reactor opera a baja temperatura (Punto 1). Cuando se aumenta la temperatura del alimento la temperatura del reactor se aumenta (Puntos 2 y 3) hasta llegar a un límite de temperatura baja (Punto 4). Si la temperatura del alimento es ligeramente aumentada aun mas, sucede un salto en la temperatura del reactor (Ignición) a una alta temperatura (Punto 8) y si a partir de
Mach
165
este punto se aumenta aun mas la temperatura de la camisa ocurre un ligero aumento en la temperatura del reactor.
Figura 9.8. Temperatura en el Reactor versus Temperatura del Alimento Contrasta el comportamiento descrito en el párrafo anterior (comenzando a una baja temperatura) con el que se observa en el caso de comenzar a una alta temperatura en el alimento. Si se comienza a una alta temperatura en el alimento (Punto 9) la temperatura del reactor es alta y disminuye con la disminución de la temperatura del alimento. Al moverse a una temperatura ligeramente inferior al límite de alta temperatura (Punto 6), la temperatura del reactor cae (Extinción) a una baja temperatura (Punto 2). Más disminución en la temperatura del alimento ocasiona pequeñas disminuciones en la temperatura del reactor. . Este comportamiento de histéresis discutido anteriormente es también conocido como comportamiento de Ignición+Extinción. Se observa que la región entre los puntos 4 y 6 son de un comportamiento inestable porque el reactor no parece operar en esta región (al menos en estado estacionario). La Figura 9.8 se construye con los archivos reactor5.m y solreactor5.m, que resuelven el sistema de ecuaciones en su forma no lineal y que aparecen al final del capítulo. Las Figuras (9.9) y (9.10) también muestran la variación de la temperatura en el reactor determinadas en una forma mas sencilla con el archivo reactor6.m. Para construir la figura (9.9) se fija la
Mach
166
temperatura en la camisa y se calcula la temperatura en el alimento para un intervalo de temperaturas en el reactor de la siguiente forma: E C As UA(Ts − Tcs ) − (−∆H )Vk o exp − RTs T fs = Ts + FρC p
(9.20)
Figura 9.9 Temperatura en el Reactor versus Temperatura en el Alimento Para construir la figura (9.10) se fija la temperatura en el alimento y se calcula la temperatura en la camisa para un intervalo de temperaturas en el reactor de la siguiente forma: E C As FρC p (Ts − T fs ) − (−∆H )Vk o exp − RTs Tcs = Ts + UA
Mach
(9.21)
167
Figura 9.10 Temperatura en el Reactor versus Temperatura en la Camisa La Figura (9.11) muestra la variación de la temperatura en el reactor con el cambio en la temperatura de la camisa para diferentes valores de la velocidad espacio
Figura 9.11. Respuesta del Reactor en función de la Velocidad Espacio
Mach
168
Se observa, en este caso, que el comportamiento del reactor en cuanto a la variación de su temperatura con la variación de la temperatura del alimento o de la camisa es un ejemplo de un cúspide catastrófica porque si se analiza esta variación para un cambio en la velocidad espacio se encuentra que cuando esta velocidad aumenta, el comportamiento en estado estacionario del reactor cambia de una respuesta monotónica a un comportamiento de múltiples estados estacionarios. La Figura (9.11) se construye con el archivo solreactor7.m y muestra el comportamiento del reactor para velocidades espacios de 0.01, 0.05, 0.15, 0.20, 0.40, 1.0 y 1.5 min-1. Para bajas velocidades espacio la respuesta del reactor es monotónica y para velocidades espacio mayores o iguales que 0.40 min-1 la respuesta es bifurcada
Análisis de entradas y salidas en el dominio de Laplace El código del archivo reactor3.m desarrolla la función de transferencia del reactor en forma matricial conociendo las matrices A, B y C del modelo expresado en la forma del espacio de los estados de acuerdo a la expresión (8.21).
G ( s ) = C ( sI − A) −1 B
(8.21)
Los resultados para la función de transferencia de la concentración y temperatura en el reactor con respecto a la temperatura en la camisa son las siguientes C A ( s) − 0.02657 − 0.0575 = 2 = T j ( s ) s + 1.412 s + 0.4627 (1.1165s + 1)(1.9357 s + 1) 0.3s + 0.3504 0.7573(0.856 s + 1) T (s) = 2 = T j ( s ) s + 1.412s + 0.4627 (1.1165s + 1)(1.9357 s + 1)
Se halla que la función de transferencia para la concentración es la de un sistema de segundo orden, mientras que la de la temperatura en el reactor contiene un numerador de primer orden y un denominador de segundo orden. Lo anterior indica que para un cambio en la temperatura en la camisa el atraso dinámico es mayor con respecto a la concentración que con respecto a la temperatura en el reactor. Esto tiene un significado físico porque un cambio en la temperatura de la camisa debe afectar primero a la temperatura antes de afectar la concentración en el reactor. En forma similar, se pueden deducir las dinámicas en las respuestas del reactor con respecto a variaciones en la concentración y temperatura en el alimento Mach
169
9.4 MATLAB: COMANDOS UTILIZADOS Comando fsolve Para el análisis en estado estacionario del reactor adiabático se desarrollan los archivos reactor1.m y solreactor1.m codificados con Matlab con los cuales se comprueban que las condiciones estacionarias determinadas dependen de las condiciones iniciales asignadas para la corrida de dichos programas. Se utiliza el comando fsolve cuya sintaxis es: x = fsolve(@reactor1, Inicio) Se observa en la sintaxis que el nombre del archivo se precede de la letra @ y adicionalmente se requiere como argumento las condiciones iniciales de las variables de estado en el reactor. Para el análisis del comportamiento del reactor mediante el diagrama de fase se definen las ecuaciones diferenciales en el archivo reactor2.m y la solución se resuelve con el archivo solreactor2.m. Las trayectorias de fase son las gráficas de temperatura en función de concentración. El programa representa las trayectorias de fase para condiciones iniciales con valores de concentraciones de 0.5 y 10 kgmol/m3 y temperaturas desde 20 ºC hasta 180 ºC para cada una de las concentraciones
Comando eig La estabilidad del rector mediante la representación del modelo en la forma del espacio de los estados requiere de la determinación de la matriz A y de sus correspondientes valores propios que se calculan con el comando eig de Matlab. Este comando se utilizó en el archivo reactor3.m y su sintaxis solo requiere como argumento la matriz A y es: Lambda = eig(A)
Comando fzero Para el análisis del comportamiento reactor adiabático en cuanto a la variación de la temperatura con el cambio en la temperatura de la alimentación en estado estacionario, se desarrollan los archivos reactor5.m y solreactor5.m codificados con
Mach
170
Matlab con los cuales se calculan los valores de las temperaturas en el reactor para temperaturas de alimento desde 1 hasta 30 ºC, asignadas de uno en uno. Se utiliza el comando fzero cuya sintaxis es: x = fzero(@reactor1, Inicio) Se observa en la sintaxis que el nombre del archivo se precede de la letra @ y adicionalmente se requiere como argumento la condición inicial de la temperatura en el reactor. El programa se plantea a sabiendas que según sea la temperatura inicial la solución de la ecuación de balance de energía en el reactor se satisface con una dos o tres respuestas.
9.5 MATLAB - PROGRAMAS CODIFICADOS Archivo reactor1.m function f = reactor1(x) global Fvs xf ko E R yf yj DH rhocp UAV f = [Fvs*(xf - x(1)) - ko*(exp(-E/(R*(x(2) + 273))))*x(1); Fvs*(yf - x(2)) + (DH/rhocp)*ko*(exp(-E/(R*(x(2) + 273))))*x(1) - (UAV/rhocp)*(x(2) - yj)];
Archivo solreactor1.m function f = solreactor1(x) global Fvs xf ko E R yf yj DH rhocp UAV Inicio Fvs = input('Valor de la velocidad espacio = '); ko = input('Valor de la constante de velocidad de reaccion = '); DH = input('Cambio de entalpia en la reaccion = '); E = input('Energia de activacion = '); rhocp = input('Densidad por calor especifico = '); yf = input('Temperatura del alimento = '); xf = input('Valor de la concentracion de A en el alimento = '); UAV = input('Coeficiente por area de transferencia de calor = '); yj = input('Temperatura de la camisa = '); R = input('Constante de los gases = '); Inicio = input('Condiciones iniciales para concentracion y temperatura = '); x = fsolve(@reactor1, Inicio)
Mach
171
Archivo reactor2.m function dx = reactor2(t,x) global Fvs Fv xf ko E R yf yj DH rhocp UAV dx = [Fv*(xf - x(1)) - ko*(exp(-E/(R*(x(2) + 273))))*x(1); Fv*(yf - x(2)) + (DH/rhocp)*ko*(exp(-E/(R*(x(2) + 273))))*x(1) - (UAV/rhocp)*(x(2) - yj)];
Archivo solreactor2.m function f = solreactor2(x) global Fv xf ko E R yf yj DH rhocp UAV Rango Inicio C Fvs = input('Valor de la velocidad espacio en estado estacionario = '); DFv = input('Cambio paso en la velocidad espacio = '); xf = input('Valor de la concentracion de A en el alimento = '); ko = input('Valor de la constante de velocidad de reaccion = '); E = input('Energia de activacion = '); R = input('Constante de los gases = '); yf = input('Temperatura del alimento = '); yj = input('Temperatura de la camisa = '); DH = input('Cambio de entalpia en la reaccion = '); rhocp = input('Densidad por calor especifico = '); UAV = input('Coeficiente por area de transferencia de calor = '); Fv = Fvs + DFv; Rango = input('Intervalo de tiempo = '); C = [0.5 10]; T = [20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180]; for i = 1:2 for j = 1:9 Inicio = [C(i) T(j)]; [t,x] = ode45('reactor2', Rango, Inicio); plot(x(:,1),x(:,2)+273) hold on end end
Archivo reactor3.m global Fv xf ko E R DH rhocp UAV x y Fv = input('Valor de la velocidad espacio en estado estacionario = '); Mach
172
ko = input('Valor de la constante de velocidad de reaccion = '); DH = input('Cambio de entalpia en la reaccion = '); E = input('Energia de activacion = '); rhocp = input('Densidad por calor especifico = '); xf = input('Valor de la concentracion de A en el alimento = '); UAV = input('Coeficiente por area de transferencia de calor = '); R = input('Constante de los gases = '); x = input('Valor de la concentracion de A en estado estacionario = '); y = input('Valor de la temperatura en el reactor en estado estacionario = '); ks = ko*exp(-E/(R*(y + 273))); ksp = ks*(E/(R*(y + 273)*(y + 273))); A(1,1) = -Fv - ks; A(1,2) = -x*ksp; A(2,1) = ks*(-DH)/rhocp; A(2,2) = -Fv -(UAV/rhocp) + (-DH/rhocp)*x*ksp; A Lambda = eig(A)
Archivo reactor4.m F = 60; V = F; ko = 9703*3600; DH = -5960; E = 11843; rhocp = 500; Tjs = 25; Caf = 10; UA = 250; Tfs = [0 5 15 21 30]; T = [0 5 10 15 16 17 27 47 67 87 97 107 117 127 137 147]; for i = 1:5 for j = 1:16 Qr(j) = (-UA*Tjs - F*rhocp*Tfs(i)) + (UA + F*rhocp)*T(j); Cas(j) = (F/V)*Caf/((F/V) + ko*(exp(-E/(1.987*(T(j)+273))))); Qg(j) = (-DH)*V*ko*(exp(-E/(1.987*(T(j)+273))))*Cas(j); end plot(T,Qr,T,Qg) hold on end
Mach
173
Archivo reactor5.m function f = reactor5(x) global F xfs xcs cafs ko E R DH rhocp UA V Inicio1 Inicio2 Inicio3 f = F*rhocp*(x - xfs) + UA*(x - xcs) - (-DH*V)*ko*(exp(-E./(R*(x + 273))))*((F/V)*cafs/((F/V) + ko*(exp(-E./(R*(x + 273))))));
Archivo solreactor5.m function f = solreactor5(x) global F V xfs xcs cafs ko E R DH rhocp UA Inicio1 Inicio2 Inicio3 a b F = 60; V = F; ko = 9703*3600; DH = -5960; E = 11843; rhocp = 500; xcs = 25; cafs = 10; UA = 250; R = 1.987; a = [-5:1:21]; for i = 1:27 xfs = a(i); if xfs <= 21 x(i) = fzero(@reactor5,30); elseif xfs ==5 x(i) = fzero(@reactor5,60); else end end plot(a,x,'k') hold on clear a x b = [5:21]; for j = 1:17 xfs = b(j); Mach
174
x(j) = fzero(@reactor5,90); end plot(b,x,'k') clear b x c = [5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21]; for j = 1:17 xfs = c(j); if xfs == 5 x(j) = fzero(@reactor5,90); else x(j) = fzero(@reactor5,120); end end plot(c,x,'k') clear c x d = [21:1:32]; for j = 1:12 xfs = d(j); x(j) = fzero(@reactor5,177); end end plot(d,x,'k')
Archivo reactor6.m global F V xfs xcs cafs ko E R DH rhocp UA F = 1; V = F; ko = 9703*3600; DH = -5960; E = 11843; rhocp = 500; cafs = 10; UA = 150; R = 1.987; x = [27:1:107]; Mach
175
c = input('Caso a estudiar x vs xcs = 1 o x vs xfs = 2 '); switch c case 1 xfs = 25; for i = 1:81 xcs(i) = x(i) + (F*rhocp*(x(i) - xfs) - (-DH*V)*ko*(exp(-E/(R*(x(i) + 273))))*((F/V)*cafs/((F/V) + ko*(exp(-E/(R*(x(i) + 273)))))))/UA; end figure(1) plot(xcs + 273,x + 273) case 2 xcs = 25; for i = 1:81 xfs(i) = x(i) + (UA*(x(i) - xcs) - (-DH*V)*ko*(exp(-E/(R*(x(i) + 273))))*((F/V)*cafs/((F/V) + ko*(exp(-E/(R*(x(i) + 273)))))))/(F*rhocp); end figure(2) plot(xfs + 273,x + 273) end
Archivo reactor7.m global F V xfs xcs cafs ko E R DH rhocp UA a x y V = 1; ko = 9703*3600; DH = -5960; E = 11843; rhocp = 500; cafs = 10; UA = 150; R = 1.987; x = [27:1:107]; a = [0.01 0.05 0.15 0.2 0.4 1 1.5]*V; xfs = 25; for j = 1:7 F = a(j);
Mach
176
for i = 1:81 xcs(i) = x(i) + (F*rhocp*(x(i) - xfs) - (-DH*V)*ko*(exp(-E/(R*(x(i) + 273))))*((F/V)*cafs/((F/V) + ko*(exp(-E/(R*(x(i) + 273)))))))/UA; end plot(xcs + 273, x + 273) hold on end
Mach