Partes: 1, 2, 3, 4, 5

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.

Datos constantes
	Presión a la entrada
	Temperatura a la entrada
	Flujo másico total a la entrada
	Fracciones molares a la entrada ()
	Presión a la salida
	Flujo másico total a la salida
	Corriente de energía a través del reactor

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.

Valores iniciales
	Fracciones molares a la salida ()
	Volumen de reacción
	Temperatura a la salida

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.

Constantes calculadas
	Flujo molar a la entrada. es la masa molecular del componente ()
	Flujos por componente a la entrada ()
	Entalpía molar a la entrada ()

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.

Variables y ecuaciones de cálculo
	Flujo molar a la salida ()
	Flujo por componente a la salida ()
	Volumen molar total a la salida ()
	Entalpía molar a la salida ()
Variables y ecuaciones de cálculo
	Concentración molar a la salida ()
	Velocidad molar de reacción (, )
	Número de moles en el interior del reactor ()
	Energía interna de la masa en el interior del reactor

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.

Balances de materia y energía
	Balances de masa por componente (, )
	Balance de energía

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.