Simulador de reactores químicos basado en Excel (página 4)

4.5. SOLUCIÓN DE ECUACIONES DE
UNA SOLA VARIABLE EN UN SIMULADOR DE REACTORES DE TANQUE
AGITADO

Es poco común que los métodos
numéricos se apliquen a las ecuaciones algebraicas
manualmente, de hecho, ésta es una práctica que se
lleva a cabo únicamente con fines académicos, en la
mayoría de los casos. El notable avance de los
computadores así como su fácil acceso al
público en general ha originado una enorme difusión
y un gran auge de los métodos
numéricos, puesto que los algoritmos son
programados y el computador se
encarga de realizar los cálculos y de mostrar
después los resultados por medio de tablas o gráficos, lo cual reduce significativamente
el tiempo de
resolución de los problemas y
evita la resolución consecutiva del método
manualmente. Estos algoritmos contienen, en ocasiones, varios
métodos numéricos incorporados y se caracterizan
por poseer varios criterios de decisión, los cuales,
dependiendo del comportamiento
de la ecuación que se desea resolver, eligen el
método más adecuado para garantizar la convergencia
a la solución.

Un algoritmo
computacional que merece especial atención por su eficacia es el
propuesto por Brent (1971) que modifica al propuesto antes por
Dekker (1969). Este algoritmo se describe con detalle a
continuación, puesto que será el empleado por el
simulador del presente trabajo de
grado para determinar la temperatura en
cada paso del intervalo de tiempo que se fije para la simulación.

4.5.1. Algoritmos de Brent y de Dekker. La
solución de la ecuación (4.1) dentro del intervalo
requiere que
(hay que anotar
que si , entonces
alguno de los extremos del intervalo es la raíz buscada).
El objetivo del
algoritmo es determinar la aproximación de la raíz real,
dentro de una tolerancia
, mediante la
evaluación de en un reducido número de puntos.
Según se observa en la expresión (4.7), cuando se
emplea el método de la bisección, la
aproximación de la solución será encontrada
en un número de iteraciones cercano a . Aunque el algoritmo de
Brent no reduce este número de pasos, si se cumple que
sea continua y
diferenciable en las cercanías de la raíz y si se
ignoran los errores de redondeo, puede afirmarse que el algoritmo
posee convergencia superlineal a la raíz, haciendo que la
técnica sea más rápida que la
bisección.

El algoritmo de Brent es similar al de Dekker y sus
variaciones, aunque es importante hacer énfasis en el
hecho de que ninguno de los últimos garantiza la
convergencia en menos de evaluaciones de la función,
mientras que el primero lo hace en cerca de evaluaciones de la
función. Esta diferencia es importante en casos en los
cuales, por ejemplo, es igual a la unidad y el valor de
es , puesto que el
método de Brent requiere 1600 evaluaciones de la
función, al tiempo que el método de Dekker requiere
evaluaciones.
Además, en funciones
polinómicas de grado moderado con raíces
adecuadamente espaciadas, el método en cuestión ha
probado ser igual o incluso más rápido que el de
Dekker, lo cual se suma al hecho de que ambos son, a su vez,
más rápidos que el de la
bisección.

– Algoritmo de Dekker. Este procedimiento
considera que, en un paso típico, se deben tener 3 puntos
, y , tales que es igual a , además y con la condición
de que . El
método asume que es la mejor aproximación a la solución
hecha hasta el paso en cuestión, además, es el valor previo de
y la raíz
debe estar entre
y . Si es equivalente a cero, la
solución ha sido encontrada, mientras que si no lo es,
entonces debe definirse la cantidad como se muestra en
(4.14).

(4.14)

Si el valor absoluto de es menor que la tolerancia, entonces es la aproximación
de la raíz buscada, mientras que si no lo es, debe
interpolarse o extrapolarse a linealmente entre y , generando un nuevo punto , que se calcula como se muestra en
(4.15).

(4.15)

Los parámetros y se definen mediante las expresiones (4.16) y (4.17),
tomando en consideración que la cantidad se define, asimismo, en
(4.18).

(4.16)

(4.17)

(4.18)

Cabe anotar que no siempre es necesario calcular el
valor de ; en
general, esta cantidad solo se calcula cuando se satisface la
condición observada en (4.19).

(4.19)

Antes de calcular la nueva aproximación de la
raíz, con la cual se continúa al siguiente paso, se
debe calcular el valor de , magnitud que se define en (4.20).

(4.20)

Finalmente, la aproximación de la raíz,
, se determina de
acuerdo con las condiciones mostradas en la expresión
(4.21).

(4.21)

De esta manera, la terna se reemplaza por la terna y el algoritmo comienza
de nuevo hasta satisfacer alguno de los criterios de
parada.

– Modificaciones de Brent al algoritmo de Dekker. La
principal modificación de Brent al algoritmo propuesto por
Dekker, consiste en asegurar que una bisección se hace al
menos una vez, cuando se realiza un número de pasos
consecutivos igual a . Dicha modificación parte de definir la magnitud
como se puede
apreciar en (4.22).

(4.22)

Si cumple la condición (4.23) o la
condición (4.24), mostradas a continuación,
entonces se hace una bisección y en otro caso, puede
hacerse una interpolación o una bisección, como en
el método de Dekker.

(4.23)

(4.24)

De acuerdo con lo anterior, se reduce al menos a la mitad de su valor
en cada paso siguiente y cuando dicho valor es menor que
, se realiza una
bisección, luego de la cual el algoritmo impone la
igualdad
(4.25).

(4.26)

Por otra parte, el algoritmo de Brent aprovecha el hecho
de que se cuente en cada paso con tres puntos para realizar una
interpolación cuadrática, la cual cuenta con mayor
precisión en sus aproximaciones que la lineal realizada
por el algoritmo de Dekker. Además de lo mencionado, la
interpolación realizada es inversa en lugar de ser directa
puesto que con ésta última es necesario resolver
una ecuación cuadrática para determinar el valor de
, lo cual aumenta
el número de cálculos.

Otra modificación de Brent al algoritmo de Dekker
consiste en evitar el error de desbordamiento o la
división por cero al calcular el valor de . Dado que es la más reciente
aproximación a la raíz y considerando que es el valor previo de
, es posible
afirmar que el método realiza una bisección cuando
se verifica la condición (4.27).

(4.27)

Si la condición (4.27) no se cumple, entonces es
válido afirmar que la condición (4.28) debe ser
verdadera y, por ende, el valor de debe ser calculado de una manera tal que
se eviten los errores antes mencionados.

(4.28)

Para calcular de manera segura el valor de la variable
, se requiere
calcular inicialmente las constantes , y ,
las cuales se definen como se muestra entre (4.29) y
(4.31).

(4.29) (4.30)

(4.31)

Una vez que las constantes , y ,
han sido calculadas, pueden deducirse los valores de
y a través de las
ecuaciones (4.32) y (4.33). Cabe anotar que estas nuevas
igualdades sustituyen en el algoritmo de Brent a las expresiones
(4.16) y (4.17), dadas por Dekker.

(4.32)
(4.33)

De nuevo, el valor de se calcula a través de la expresión
(4.15), sin olvidar que si la condición (4.19) se cumple,
entonces el método debe hacer una bisección en
lugar de una interpolación.

Sumada a sus demás virtudes, el algoritmo de
Brent posee una importante característica: cuenta con
convergencia superlineal a la raíz. Este concepto se
asocia matemáticamente con la capacidad de un
método para satisfacer la expresión (4.34), toda
vez que sea de
clase 1
(), es decir, que
ella y su primer derivada sean continuas, en los alrededores de
.

(4.34)

La convergencia superlineal del método de Brent
hace que el orden de convergencia aumente, como mínimo, en
un número de veces igual a , comparado con los métodos que emplean
interpolación lineal. Todo lo anterior hace que el
método de Brent sea tanto o más rápido que
el de Dekker, con la ventaja de que el primero solo
requirió para su convergencia, durante las pruebas
reportadas en la literatura, un número
de evaluaciones de la función nunca mayor al triple de las
que se hubieran empleado si se hubiese recurrido al método
de la bisección, mientras que el método de Dekker
requiere, en algunas ocasiones, más de evaluaciones de la
función.

-Tolerancia en el método de Brent. En este punto
es apropiado comentar la forma mediante la cual opera la
tolerancia en el método de Brent. En esta técnica,
la tolerancia ()
es una combinación de una tolerancia relativa () y una tolerancia
absoluta (), de
manera tal que, en cada paso, se tiene la igualdad mostrada en
(4.35).

(4.35)

En la expresión (4.35) se supone que es el valor de la
aproximación en un paso determinado, al tiempo que
es una
tolerancia absoluta positiva cuyo valor es constante.
Además,
se define como la precisión relativa de la máquina
en la cual se programa el
método, la cual puede calcularse como se observa en
(4.36).

(4.36)

En la ecuación (4.36) se ha establecido que
es el
número de dígitos en que trunca la
aritmética de punto flotante cuya base es , aunque si la
aritmética empleada es de redondeo, se emplea la mitad del
valor obtenido a través de la ecuación
(4.36).

Como puede observarse en la igualdad (4.35), depende del valor de la
aproximación en cada paso, la cual puede estar en
cualquier lugar del intervalo de solución, razón
por la que puede afirmarse que el método de Brent posee
una tolerancia variable que asume un nuevo valor en cada paso. Si
se desea que esta técnica numérica le asigne un
único valor a la tolerancia, debe reemplazarse a por su mínimo
valor positivo dentro del intervalo, sin olvidar que dicho valor
no puede excesivamente pequeño, puesto que ello
obligaría a realizar un elevado número de
evaluaciones de la función, lo cual no es conveniente al
tener en cuenta que este hecho es, precisamente, uno de los
criterios de parada de la técnica.

– Modificación al algoritmo de Brent. El
algoritmo desarrollado por Brent se conoce como zero y
puede ser programado con relativa facilidad* . Sin embargo, los autores del presente trabajo
de grado modificaron este algoritmo sustancialmente, de manera
tal que el tiempo de cálculo se
redujo considerablemente en las pruebas realizadas con el
simulador. La modificación es simple: como se verá
en el capítulo 5, el sistema de
ecuaciones conformado por los balances de materia y
energía y las relaciones adicionales, posee un carácter algebraico-diferencial que
requiere, para su solución, un conjunto de condiciones
iniciales. Dentro de este conjunto se encuentra un valor inicial
de la temperatura en el interior del reactor, la cual se
suministra al método de Brent con la intención de
que, en lugar de partir de la igualdad de con , emplee siempre tal valor
inicial como la primera aproximación de la raíz.
Con esta modificación se logra que el intervalo en el cual
se busca la solución (definido entre 1 y 800) se reduzca
rápidamente en cada paso del intervalo de tiempo, haciendo
que la simulación se concluya en un tiempo menor que con
el método original y cerca de 3 veces más
rápido que con la técnica de Newton-Raphson
antes expuesta, con la ventaja de haberse evitado el
cálculo de la derivada numérica de la
función, lo cual suele introducir errores significativos
en el cálculo de la aproximación.

5.
MODELO GENERAL
PARA LA SIMULACIÓN DE REACTORES
DE TANQUE AGITADO

1. DEFINICIONES PREVIAS

1. 

   Figura 2. Volumen de control seleccionado
   para el reactor que será simulado.

   Fuente: SMITH, C. A. y CORRIPIO, A. B.
   Principles and practice of automatic process control. 2nd
   Ed. New York: John Wiley & Sons, 1997. p.
   129

2. Volumen de
   control.
   En este caso es el volumen en el interior del reactor que se
   encuentra ocupado por el líquido. Esto implica que el
   volumen de control está lleno de líquido en
   todo momento, que la frontera
   del sistema es móvil y por tanto la superficie de
   control es variable. En la Figura 2 se observa el volumen de
   control seleccionado, se observa también que el
   reactor cuenta con una corriente de enfriamiento o de
   calentamiento (según el caso), con un agitador y con
   una válvula en la corriente de salida.
3. Marco de referencia. El marco de referencia
   empleado es inercial y su sistema coordenado tiene su origen
   en el centro de la circunferencia que conforma la base del
   reactor. Se emplean, debido a la forma del reactor,
   coordenadas cilíndricas suponiendo que el equipo
   está orientado verticalmente en el eje Z.

En la Figura 3 se muestra la ubicación del origen
del sistema de coordenadas.

Figura 3. Ubicación del origen de
coordenadas dentro del volumen de control.

Fuente: SMITH, C. A. y CORRIPIO, A. B. Principles
and practice of automatic process control. 2nd Ed. New York: John
Wiley & Sons, 1997. p. 458.

5.2. ECUACIONES DE LAS LEYES DE
CONSERVACIÓN,

5.2.1. Balances de materia por componente. El balance de
materia para un componente en un volumen de control en el cual ocurre un
número determinado reacciones químicas se presenta
a continuación en la ecuación (5.1).

(5.1)

Donde:

:
Flujo molar de entrada al reactor del componente
i.

:
Flujo molar de salida del reactor del componente
i.

:
Volumen de líquido en el interior del
reactor.

:
Coeficiente estequiométrico del componente i en
la reacción j.

:
Velocidad
global de la reacción j.

:
Número total de moles del componente i en el
interior del reactor.

:
Número de componentes presentes en la corriente de
alimentación o en el interior del
reactor.

:
Número de reacciones que ocurren en el interior del
reactor.

– Suposiciones de los balances de materia:

• El estado de la
  masa dentro del volumen de control puede cambiar con el tiempo,
  pero en cualquier instante el estado es
  uniforme en todo el volumen de control.
• Las propiedades en el interior del reactor son
  iguales a las propiedades de la corriente de salida (estado
  uniforme entre el interior del rector y la salida).
• El reactor y sus corrientes están conformados
  por una sola fase líquida.
• El reactor es un cilindro circunferencial
  recto.
• El reactor se encuentra abierto a la atmósfera.
• La tubería de salida está conectada al
  fondo del reactor.
• La corriente de alimentación está
  ubicada en la parte superior del reactor.

2. Balance de energía. La forma general del
   balance de energía (primera ley de la
   termodinámica) para un volumen de
   control en estado transitorio es la indicada en
   (5.2).

(5.2)

Donde:

: Flujo
molar total de entrada al reactor.

:
Flujo molar total de salida del reactor.

:
Energía molar del fluido que entra al volumen de
control.

:
Energía molar del fluido que sale del volumen de
control.

:
Energía total del fluido que se encuentra dentro del
volumen de control.

:
Calor neto
que entra al volumen de control desde los
alrededores.

:
Trabajo neto que entra al volumen de control desde los
alrededores.

– Suposiciones del balance de energía:

• El trabajo de eje posee un valor despreciable frente
  a las demás formas de energía.
• El sistema carece de efectos eléctricos,
  magnéticos, gravitacionales, de movimiento y
  de tensión superficial (sistema compresible
  simple).
• La energía cinética y la energía
  potencial poseen valores
  despreciables frente a las demás formas de
  energía.

El balance de energía, teniendo en cuenta las
simplificaciones que introducen las suposiciones, queda de la
manera mostrada en (5.3).

(5.3)

Donde:

:
Energía interna molar del fluido que entra al volumen de
control.

:
Energía interna molar del fluido que sale del volumen de
control.

:
Energía interna total del fluido que se encuentra
dentro del volumen de control.

:
Trabajo de flujo neto que entra al volumen de control desde
los alrededores.

Reemplazando el trabajo de
flujo en la ecuación anterior, queda la expresión
(5.4).

(5.4)

Donde:

:
Presión
de la corriente de alimentación.

:
Presión de la corriente de salida.

:
Volumen molar de la corriente de
alimentación.

:
Volumen molar de la corriente de salida.

El balance de energía, finalmente, de acuerdo con
la definición de entalpía, tiene la estructura que
se muestra a continuación en (5.5).

(5.5)

Donde:

:
Entalpía molar del fluido que entra al volumen de
control.

:
Entalpía molar del fluido que entra al volumen de
control.

El flujo másico de entrada y su
composición son constantes, se supondrá
además que el flujo másico a la salida es
constante, sin embargo, los balances tienen unidades molares, por
lo cual es necesario establecer una relación entre ambos
flujos, como se muestra en la ecuación (5.6).

(5.6)

5.3. RELACIONES ADICIONALES,

Las relaciones adicionales son expresiones algebraicas
que permiten relacionar algunas cantidades presentes en los
balances entre sí y con ello se reducen los grados de
libertad del
modelo. Estas relaciones se listan a continuación y sus
representaciones matemáticas son las ecuaciones (5.7) a
(5.16).

– Definiciones de los flujos globales:

(5.7) (5.8)

– El flujo de salida está controlado por una
válvula:

(5.9)

Es preciso mencionar que, dentro del simular
desarrollado, el flujo másico de salida se asume como un
dato fijo a lo largo de la simulación, lo cual equivale a
decir que puede
mantenerse constante y que la válvula establece un flujo
determinado que no varía; esto es cierto siempre y cuando
se garantice un caudal mínimo de operación que
depende, entre otros factores, del tamaño de la
tubería.

– Las expresiones para la velocidad de reacción
poseen la siguiente forma general (aunque su estructura matemática
será expuesta en la sección 5.4):

(5.10)

– La concentración se relaciona con los flujos
mediante el volumen molar:

(5.11)

– Las moles a la salida dependen del volumen en el
reactor y la concentración a la salida (la cual es la
misma que en el interior del reactor):

(5.12)

– La energía interna se relaciona con la
entalpía mediante la definición
siguiente:

(5.13)

– De acuerdo con la termodinámica, si se conocen
la temperatura, la presión y la composición de una
mezcla se pueden hallar las demás propiedades
termodinámicas, por lo cual las entalpías y el
volumen molar a la salida están dadas por relaciones entre
estas cantidades:

(5.14)

(5.15)

(5.16)

Donde:

:
Concentración molar de salida al reactor del componente
i.

:
Volumen molar a la salida del reactor

:
Número total de moles en el interior del
reactor.

5.4. ESTRUCTURA MATEMÁTICA DE LAS LEYES DE
VELOCIDAD

La ecuación (5.10) únicamente enuncia las
variables de
las cuales dependen las leyes de velocidad que se llevan a cabo
en el interior del reactor, aunque no muestra la manera
explícita como se estructura matemáticamente esta
dependencia. En su caso más general, las reacciones
químicas que tengan lugar dentro del reactor simulado
deben seguir la ley de velocidad de que se muestra en la
ecuación (5.17), la cual se conoce como ley de velocidad
catalítica heterogénea reversible. Dado que este es
el caso más general de ley de velocidad que puede
simularse, es posible simplificar dicha ecuación de manera
que se ajuste la ley de velocidad según la
situación.

(5.17)

Los términos presentes en la ley de velocidad
(5.17) y siguientes se definen a continuación:

• ,
  y : Factores de
  frecuencia de la reacción directa, de la inversa y de
  los pasos intermedios de la reacción,
  respectivamente.
• ,
  y : Energías de
  activación de la reacción directa, de la inversa
  y de los pasos intermedios de la reacción,
  respectivamente.
• ,
  y : Funciones de la
  concentración de las especies que reaccionan en la
  reacción directa, en la inversa y en los pasos
  intermedios de la reacción, respectivamente.
• :
  Temperatura en el interior del reactor.
• :
  Constante universal de los gases
  ().
• 

: Constante que depende del tipo de mecanismo que siga
la reacción.

La primera simplificación que puede realizarse
resulta de afirmar que la reacción cuya ley de velocidad
se está representando no es reversible, con lo cual se
obtiene la ley de velocidad mostrada en (5.18).

(5.18)

Por otra parte, si se establece que la reacción
es reversible aunque no catalítica heterogénea, el
denominador de la ecuación (5.17) se hace igual a la
unidad y se llega a la expresión (5.19).

(5.19)

Finalmente, si la reacción no es
catalítica heterogénea ni reversible, la
ecuación (5.17) se simplifica a su versión
más simple, la cual se esquematiza en la expresión
(5.20).

(5.20)

Es de especial importancia mencionar la estructura
matemática común de los términos presentes
en las anteriores leyes de velocidad, , y ,
puesto que es indispensable especificarla correctamente en el
simulador Retromezcla, como se verá en el capítulo
6. En este simulador, tales funciones están conformadas
por una productoria de las composiciones molares de las
sustancias, elevadas a algún exponente real (), como se aprecia en la
ecuación (5.21).

(5.21)

Como se pudo observar en esta sección, el
simulador permite ingresar un conjunto amplio de reacciones e
incluso varias reacciones al tiempo, siempre que se tengan los
datos
cinéticos suficientes en lo que respecta a factores de
frecuencia, energías de activación, constantes y
relaciones de la concentración de las especies con la
velocidad de reacción. Más adelante, en el
capítulo 6, se verá la forma mediante la cual se
ingresan los datos cinéticos en el simulador, así
como los demás datos necesarios para realizar la
simulación.

6.
RETROMEZCLA: SIMULADOR DE REACTORES DE TANQUE
AGITADO

6.1. DESCRIPCIÓN DEL PROBLEMA

Como pudo observarse en el capítulo 5, la
naturaleza
dinámica del objeto de estudio le confiere
un carácter diferencial a las ecuaciones que representan
los balances de materia y energía y su solución
requiere un método numérico para ecuaciones
diferenciales. Por otro lado, los métodos
termodinámicos son relaciones complicadas que en caso de
resolverse para la temperatura, precisan un método
numérico para ecuaciones algebraicas. Cada estrategia
numérica puede ser implantada por separado de forma
relativamente sencilla para lograr un objetivo específico,
sin embargo, cuando un sistema de ecuaciones es
simultáneamente algebraico y diferencial, la
solución tiene un mayor grado de dificultad y puede variar
de un problema a otro. En general, existen dos formas
prácticas de encontrar sistemas
algebraico-diferenciales: cuando se aplican métodos
numéricos implícitos multipaso del tipo
predictor-corrector para la resolución de un sistema de
ecuaciones diferenciales o toda vez que se pretenda resolver el
sistema de ecuaciones resultante del análisis de un sistema
dinámico.

Por un lado, la elección de un método
numérico es, en esencia, un problema matemático que
puede resolverse según el criterio de aquel que lo emplea,
puesto que la aplicación de un método
numérico explícito requiere menos cálculos
por paso que el implícito aunque su tamaño de paso
se encuentra limitado por consideraciones acerca de la
estabilidad y por ello es recomendable emplear un método
implícito cuando el tamaño de paso es
pequeño. Sin embargo, cuando se elige un método
predictor-corrector, se aplica un método explícito
para generar una predicción de la aproximación que
es utilizada como valor inicial en la iteración con el
método implícito con el cual se determina la
verdadera aproximación de la solución. Este tipo de
métodos tiene un alto grado de desarrollo
pues incorporan estimaciones automáticas del error que
optimizan el orden de aproximación de las derivadas, el
tamaño de paso y la frecuencia con que se calcula la
matriz
jacobiana en métodos como el de Newton, con la
única desventaja de tener que resolver en cada paso un
sistema algebraico de ecuaciones, cuyo procedimiento de
solución puede resultar complicado debido al
carácter implícito de las variables. De cualquier
manera, la aplicación de los métodos
predictor-corrector no es obligatoria y según la
dificultad del problema al que puedan aplicarse y teniendo en
cuenta el grado de precisión requerido, pueden sustituirse
por métodos más sencillos.

Por otro lado, en el caso de los sistemas
dinámicos, la solución conjunta de ecuaciones
algebraicas y diferenciales no es una cuestión opcional
puesto que el carácter dinámico de un sistema se
asocia matemáticamente con la derivada, mientras que las
relaciones adicionales entre las variables suelen ser de
naturaleza algebraica. En particular, como se observó en
el capítulo 5, para el caso del diseño
de reactores en estado transitorio, los balances de masa y
energía son ecuaciones
diferenciales ordinarias cuya variable independiente es el
tiempo, mientras que las relaciones volumétricas, las
expresiones cinéticas y las relaciones entre propiedades
físicas son expresiones algebraicas entre los valores
instantáneos de las variables.

Para resolver un sistema de ecuaciones de este estilo,
Reklaitis et al resaltan la necesidad de eliminar de los balances
las variables cuyas derivadas no aparecen en las ecuaciones
diferenciales (llamadas variables algebraicas), para que el
sistema resultante solo se encuentre en términos de las
variables cuyas derivadas sí aparecen en ellas (llamadas
variables diferenciales). Esta eliminación puede ser
simbólica (manipulando las ecuaciones), aunque resulta
equivalente hacerla numéricamente mediante la
ejecución de un algoritmo. Otra alternativa propuesta por
Reklaitis et al consiste en derivar el sistema de ecuaciones
algebraico con respecto al tiempo para obtener un sistema
completamente diferencial que puede ser resuelto por alguno de
los métodos descritos en el capítulo 3, suponiendo
que el sistema generado tenga solución, lo cual no siempre
es cierto, como mostraron los autores citados; aunque, en
ocasiones, pueden existir soluciones
locales a pesar de que la matriz jacobiana del sistema sea
singular cuando una o varias ecuaciones algebraicas son
derivadas. En cuanto a este tema, Reklaitis et al afirman que
cuando existe al menos una variable algebraica que solo aparece
en las ecuaciones diferenciales, el sistema generalmente no tiene
solución, aunque es posible que al derivar ciertas
ecuaciones y realizar ciertas manipulaciones matemáticas,
se obtengan nuevas expresiones algebraicas que permitan llegar a
la solución del sistema algebraico-diferencial resultante.
Este proceso
derivativo orientado a la obtención de nuevas relaciones
independientes entre las variables algebraicas se conoce como
algoritmo de reducción de índice de Gear-Petzold,
en el cual se establece que el índice del sistema
corresponde al número de veces que es necesario derivar
las ecuaciones para alcanzar una solución.

6.2. ALGORITMO DE SOLUCIÓN
DESARROLLADO

El problema descrito en la sección 6.1 requiere
para resolverse, un algoritmo que muestre a través de una
sucesión de pasos, el proceso mediante el cual se
administran la información y las herramientas
numéricas disponibles, de manera que se consiga,
finalmente, una aproximación de la solución en cada
paso del intervalo de tiempo especificado. Los pasos del
algoritmo mostrado a continuación permiten encontrar, a
través del tiempo, la temperatura y la composición
en el interior del reactor, así como el volumen de
reacción y las moles que lo componen.

a) Se ingresan los valores que permanecerán
constantes a lo largo del proceso, los cuales se muestran en la
Tabla 2, de acuerdo con la notación empleada en el
capítulo 5. Cabe mencionar que la suposición de que
estos valores permanezcan constantes es inherente a este
algoritmo y no propiamente al modelo expuesto en el
capítulo 5. Es necesario ingresar también el
intervalo de tiempo y el tamaño del paso.

Tabla 2. Datos constantes a través del
proceso.

b) Se fijan los valores iniciales de las condiciones de
salida que son variables a través del tiempo, así
como el volumen inicial de reacción, como se muestra a
continuación en la Tabla 3.

Tabla 3. Valores iniciales de la salida y del
interior del reactor.

c) A partir de los valores ingresados anteriormente en
el paso a), se calculan las magnitudes dependientes de dichos
valores, las cuales se listan en la Tabla 4, mediante las
fórmulas que se indican en la segunda columna de dicha
Tabla. Cabe anotar que, debido a que estas magnitudes dependen de
valores constantes, de la misma manera su comportamiento
será constante a través del tiempo.

Tabla 4. Valores constantes
calculados.

d) Con base en los valores especificados en a) y b), se
calculan los valores iniciales de las demás variables
necesarias para resolver numéricamente los balances de
masa y energía. Las variables calculadas junto con sus
ecuaciones de cálculo pueden observarse en la Tabla
5.

Tabla 5. Valores calculados de las
variables.

e) Con todos los valores iniciales especificados o
calculados anteriormente, se resuelve el siguiente sistema de
ecuaciones diferenciales presentado en la Tabla 6 para las
variables y
, mediante alguno
de los métodos numéricos para ecuaciones
diferenciales que fueron expuestos en el capítulo
3.

Tabla 6. Sistema de ecuaciones diferenciales del
volumen de control.

Es importante mencionar que para la
implementación en MS Excel® de
este algoritmo, se empleó el método de Runge-Kutta
de cuarto orden (RK4). Dado que el método RK4 requiere el
cálculo de 4 constantes intermedias antes de generar la
aproximación real de la solución, después de
calcular cada constante, deben realizarse los cálculos
descritos a continuación en el paso f) y recalcularse los
valores mencionados en la Tabla 5.

f) Empleando los valores obtenidos en el paso e), se
calculan los valores de las variables a través del tiempo,
después del tiempo inicial. Inicialmente se calculan las
composiciones de los componentes mediante la ecuación
(6.1).

(6.1)

La temperatura del paso 1) se determina aplicando un
método numérico de una variable a la
ecuación (5.12), puesto que en tal ecuación la
temperatura es una variable implícita.

(5.12)

Para el desarrollo del simulador en MS Excel®, el
método numérico algebraico seleccionado fue el de
Brent, el cual fue expuesto en el capítulo 4.

Una vez que se ha determinado la temperatura a la
salida, se calcula el volumen molar de la mezcla en el interior
del reactor y con estos datos puede calcularse el volumen de
reacción reacomodando la expresión (5.11) de manera
que dicha variable sea explícita, obteniéndose la
ecuación (6.2).

(6.2)

g) Se repite el algoritmo desde el paso d), para cada
uno de los valores del tiempo hasta el final del intervalo
definido, teniendo en cuenta el tamaño del paso fijado en
el paso a).

En cualquier paso del algoritmo presentado es posible
que se encuentren valores negativos para algunas de las variables
calculadas. Dado que la entidad en cuestión es la masa y
sabiendo que las temperaturas y las presiones son absolutas, no
es coherente obtener algún valor menor que cero para
ninguna variable, de manera que cuando esto sucede, se debe
interrumpir el algoritmo y revisar los cálculos, las
especificaciones o las leyes de velocidad, a fin de detectar el
error. En ocasiones, la obtención de valores negativos al
calcular el volumen de reacción está asociada con
el consumo total
del reactivo límite.

El anterior algoritmo fue programado en VBA (Visual Basic para
Aplicaciones) a manera de complemento en MS Excel® por los
autores del presente trabajo de grado, dando origen al simulador
dinámico de reactores de tanque agitado Retromezcla, en el
cual pueden simularse reactores continuos (CSTR), con alguna
entrada o alguna salida (semibatch) o discontinuos
(Batch).

6.3. GUÍA PARA UTILIZAR EL SIMULADOR
RETROMEZCLA

Ahora se mostrará la operación del
simulador Retromezcla a través de la simulación de
un sistema hipotético de 2 reacciones en serie formado por
una dimerización irreversible seguida de una
isomerización irreversible, ambas en fase líquida,
cuya estequiometría se muestra en (6.3) y (6.4)
(las masas moleculares son 46, 92 y 92 para , y
,
respectivamente).

(6.3)

(6.4)

La velocidad de reacción de la primera
reacción sigue una cinética de segundo orden que da
origen a la ecuación (6.5), mientras que la segunda
reacción sigue la cinética de primer orden descrita
por (6.6).

(6.5)

(6.6)

La concentración se encuentra en unidades de
, mientras que
las velocidades de reacción tienen unidades de . Las reacciones se
efectuarán en un CSTR adiabático cuyos flujos de
entrada y salida equivalen a 36. Las presiones a la entrada y a la salida son
120; la
temperatura de la alimentación, así como la inicial
en el interior del reactor, es 350. La alimentación contiene a puro y el reactor cuenta
con un volumen inicial de 0.5, en el cual se encuentra únicamente . La energía de
activación de la primera reacción tiene un valor de
7000, al tiempo
que la segunda es 6000.

En este caso se supondrá que las propiedades del
reactivo pueden
aproximarse a las del etanol; se supondrá además
que las propiedades de la sustancia pueden aproximarse a las del 1-butanol y
que las del producto
pueden
aproximarse a las del 1,3-butadieno. Por último, se
supondrá que la solución formada por esta terna de
sustancias es ideal.

6.3.1. Presentación. Para comenzar, la Figura 4
muestra la presentación del simulador Retromezcla, el cual
funciona como complemento de MS Excel® y fue programado en
VBA (Visual Basic for
Applications). En ella se pueden apreciar los nombres de los
autores del presente trabajo de grado, así como los
directores del mismo. A continuación, se explicarán
brevemente las partes del simulador Retromezcla, así como
la secuencia de pasos necesaria para llevar a cabo la
simulación deseada.

Figura 4. Presentación del
simulador Retromezcla.

Es preciso señalar que cada sección
contiene una ventana de ayuda que se encuentra siempre a la vista
del usuario y es posible abrirla para su lectura
haciendo clic con el Mouse sobre
ella. Cada ventana de ayuda explica brevemente lo que debe
hacerse en la sección a la que pertenece. Un ejemplo de
dicha ventana se muestra en la Figura 5.

Figura 5. Ventana de ayuda del
simulador Retromezcla.

6.3.2. Sustancias. Una vez pasada la ventana de
presentación, el simulador se sitúa en la ventana
"Sustancias", mostrada en la Figura 6, a través de la cual
es posible seleccionar de la base de datos
el conjunto de sustancias que intervendrán en la
simulación, creándose una lista de trabajo (en este
caso se adicionan las sustancias hipotéticas , y ).

Figura 6. Sección
"Sustancias" del simulador Retromezcla.

La adición de sustancias a la lista de trabajo
puede realizarse a través del botón "Agregar" o
haciendo un doble clic con el Mouse sobre el nombre de la
sustancia en la base de datos. La eliminación de
sustancias de la lista de trabajo puede realizarse
análogamente a través del botón "Eliminar" o
haciendo un doble clic sobre el nombre de la sustancia en la
lista de trabajo. Si se desea verificar si una sustancia se
encuentra dentro de la base de datos, puede emplearse el buscador
de sustancias, que es la casilla bajo el mensaje "Buscar una
sustancia en la base de datos"; en esta casilla se escribe el
nombre sistemático o común de la sustancia
requerida y si está en ella, aparece seleccionada en ella
inmediatamente. En la ventana "Sustancias" existen también
otros 3 botones: "Adicionar compuesto", "Coeficientes binarios" y
"Editar sustancia". Las funciones de cada uno de ellos,
así como la forma de emplearlos se explican a
continuación.

– Adicionar compuesto. Cuando se presiona este
botón se accede a la ventana que se muestra en la Figura
7. Ella permite adicionar una nueva sustancia y sus propiedades a
la base de datos del programa, de manera que pueda emplearse en
la simulación, luego de ingresar las propiedades
requeridas y de presionar el botón "Guardar".

Figura 7. Ventana "Adicionar
compuesto" del simulador Retromezcla.

Siempre que se quiera adicionar una nueva sustancia debe
ingresarse obligatoriamente su nombre, su peso molecular en
, su temperatura
crítica
en , su
presión crítica en , su constante A de Bureš en , su constante B de
Bureš en ,
su constante C de Bureš en . Por otra parte, puede escogerse entre la
determinación de los parámetros empíricos de
la ecuación PTWS, y ,
o ingresar el factor acéntrico de la sustancia, aunque
esta última opción solo se recomienda para
sustancias no polares o ligeramente polares, como se
sugirió en el capítulo 2.

– Coeficientes binarios. Al apretar este botón se
accede a la ventana que se observa en la Figura 8. En ella deben
seleccionarse dos sustancias (Sustancia 1 y Sustancia 2),
presentes en la base de datos, para conocer o modificar sus
coeficientes binarios de 3 tipos: coeficientes de interacción binarios de la regla de
Wong-Sandler, coeficientes "aij" del modelo de
actividad NRTL y coeficientes "ij", también del modelo NRTL.
Si no se adiciona ningún valor para estas cantidades, el
programa asume que son cero.

Figura 8. Ventana "Coeficientes
binarios" del simulador Retromezcla.