Información Tecnológica Vol. 21(4), 3-10 (2010) doi:10.1612/inf.tecnol.4341 it.09

 

MATEMATICAS APLICADAS

 

Una Solución de la Ecuación de Rachford-Rice para Sistemas Multifases Aplicando el Método Newton-Raphson, un Parámetro de Broyden y el Flash Negativo

 

A Solution of the Rachford-Rice Equation for Multiphase Systems by Using the Newton-Raphson Method, Broyden Parameter and the Negative Flash

 

Miguel A. Mueses y Fiderman Machuca

Escuela de Ingeniería Química, Grupo de Investigación en Procesos Avanzados de Oxidación GAOX Universidad del Valle, Ciudad Universitaria Meléndez, A. A. 25360, Cali-Colombia. (e-mail: mgmueses@univalle.edu.co, fiderman@univalle.edu.co)

---

Resumen

Se obtuvo una solución de la ecuación de Rachford Rice Generalizada (R-RG) para múltiples fases y múltiples componentes, usando un planteamiento numérico que acopló un método modificado de Newton-Raphson, un parámetro de amortiguamiento tipo Broyden (obtenido a partir de la Norma Euclidiana del vector de funciones residuales R-RG) y el concepto de Flash negativo, el cual condicionó la naturaleza positiva de las fracciones de los componentes en las fases resultantes. Esta solución se probó con sistemas hipotéticos de P fases y N componentes generados aleatoriamente, al igual que los vectores de valor inicial. Se encontró que la solución propuesta es altamente estable y convergente para cualquier tipo de vector de inicio y que el número de iteraciones es afectado críticamente por éste vector y no necesariamente por la cantidad de fases o número de componentes.

Palabras Clave: ecuación Rachford-Rice, parámetro de Broyden, norma euclidiana, flash negativo

---

Abstract

A solution is obtained for the generalized Rachford-Rice equation (GR-R) for multiple phases and multicomponent systems by using a numeric approach that coupled a modified Newton-Raphson’s method, a Broyden-type damping parameter (obtained from the Euclidian vector norm of the GR-R residual functions) and the Negative Flash, which conditioned the positive nature of the component fractions in the resulting phases. This solution was tested on hypothetical systems of P phases and N components randomly generated, and the same was done with the initial value vectors. It was found that the proposed solution is highly stable and converges for any type of the initial value vector and that the number of iterations is critically affected by this vector and not necessarily for the amount of phases or the number of components.

Keywords: Rachford-Rice Equation, Broyden parameter, euclidean norm, negative flash

---

 

INTRODUCCIÓN

Convencionalmente, la solución de la ecuación de Rachford-Rice ha sido abordada bajo dos enfoques: uno por solución no lineal de las ecuaciones de separación (Leibovici y Neoschil, 1995; Nelson, 1987) y otro por minimización de la energía de Gibbs del sistema (Leibovici y Nichita, 2008; Nichita et al, 2006; Vaca et al, 2006; Tiscarreño et al, 1998; Michelsen, 1994). En los dos enfoques el principal problema es la definición del vector de valor inicial dado que si este no se encuentra cerca de la región factible, el sistema no converge. Los esquemas convencionales aplicados en la solución de estos sistemas como Newton-Raphson original (NR), sustitución sucesiva o método de la secante, son fácilmente vulnerables si los valores de arranque no están dentro de la región factible o no son cercanos a la solución puntual del sistema. Para resolver esta limitación se ha implementado, para el caso de dos y tres fases, una restricción de la región de solución por planos de acotamiento generados por análisis de los coeficientes de distribución (Leibovici y Nichita, 2008; Nichita et al, 2006; Nelson, 1987). Para cuatro o más fases, se ha propuesto un parámetro potencial en función del flash negativo (Leibovici y Neoschil, 1995) pero no presenta relación matemática con las ecuaciones R-RG. Existe otra solución por aproximación geométrica que redefine el concepto de flash negativo (Juanes, 2008), no obstante sus restricciones matemáticas la convierten en una aplicación específica.

En este trabajo se presenta un esquema de solución para el sistema de funciones residuales no lineales de Rachford-Rice, para múltiples fases en equilibrio y múltiples componentes, mediante un análisis numérico acoplando un método modificado de Newton-Raphson con un parámetro de amortiguamiento tipo Broyden (Broyden, 1965), obtenido por aproximación cuadrática de la Norma Euclidiana del vector de funciones residuales R-RG y el Flash negativo (Whitson y Michelsen, 1989; Leibovici y Neoschil, 1995), que condiciona la naturaleza positiva de las fracciones de los componentes en las fases reales o hipotéticas resultantes del proceso de separación instantánea.

ECUACIÓN DE RACHFORD-RICE GENERALIZADA (R-RG)

La ecuación de Rachford-Rice Generalizada, R-RG, puede ser obtenida utilizando la metodología de las ecuaciones MES (mass equilibrium stoichiometry) (Henley y Seader, 1988), en un separador instantáneo (flash) multifásico hipotético (con, k = 1,..., P fases), para una mezcla de N componentes con composición global {z_i} (con i = 1... N).

Al aplicar el balance de materia (masa o molar) total al sistema se obtiene:

Donde Ψ_k = p/F, es la fracción de fase y p es su flujo de masa. Análogamente, el balance de materia por componente, se expresa como una función adimensional, en términos de la composición x del componente i en la fase k, así:

De la termodinámica de equilibrio de fases sea K_k,i el coeficiente de distribución (Henley y Seader, 1988) de la fase k para el componente i y tomando una fase R no nula, e independiente del alimento F, como fase de referencia (para 1 ≤ R ≤ P), se obtiene:

Expandiendo y acoplando las ecuaciones (1) y (2) para la fase R y utilizando la ecuación (3), se obtiene la composición de cada fase:

Por condición estequiométrica para cada fase P la suma de las fracciones de masa debe ser igual a 1. Restando la suma de cada fase a la fase de referencia se obtiene un grupo de (P-1) ecuaciones residuales o de discrepancia, correspondientes a la ecuación de Rachford-Rice Generalizada para P fases y N componentes (Leibovici y Neoschil, 1995):

Esta ecuación es una función generatriz de las ecuaciones flash, donde, E, es una función residual de la fase k. Por tanto, las ecuaciones que describen la separación "flash" comprenden un sistema de (P-1) ecuaciones, una por cada fracción de fase Ψ. Las fracciones de los componentes de fase, χ, se calculan utilizando la ecuación (4).

RESTRICCIONES DE LA ECUACIÓN R-RG

El denominador de la ecuación R-RG, debe cumplir con las siguientes restricciones matemáticas:

La igualdad a cero de la ecuación (6) corresponde a una indeterminación de la ecuación R-RG, lo cual implica que el espacio de solución factible es confinado en hipercubos para el grupo de valores de Ψ_k. Para el caso de P = 2 y 3 fases, la igualdad a cero, permite calcular los coeficientes de distribución (K_max y K_min) que acotan la solución (Leibovici y Nichita, 2008; Nichita et al, 2006, Vaca et al, 2006; Tiscarreño et al, 1998; Michelsen, 1994). Para P fases, la nulidad de la ecuación permite establecer el grupo de valores de la fracción de fase que delimitan la solución factible (Leibovici y Neoschil, 1995). La desigualdad positiva (>0) corresponde a la restricción impuesta por la naturaleza física de las fracciones de los componentes en cada fase. Esta condición es conocida como la condición de no negatividad o flash negativo, X_k¡ ≥ 0 (Whitson y Michelsen, 1989).

SOLUCIÓN DE LA ECUACIÓN R-RG

La naturaleza no lineal de la ecuación R-RG implica alta vulnerabilidad al utilizar métodos convencionales de solución iterativa. Sin embargo es posible utilizar un método modificado de Newton-Raphson multivariable, que para el sistema de (P-1) funciones residuales se exprese para cada iteración (m), como un sistema de ecuaciones lineales de la forma:

Donde ΔΨ_k^(m) es la solución del sistema lineal, E_k(Ψ_k^(m)) es la evaluación del las ecuaciones R-RG en la iteración m y la sumatoria corresponde al Jacobiano de las funciones residuales. La solución de la ecuación R-RG, está dada por la forma modificada de ΔΨ_k^(m)(Leibovici y Neoschil, 1995):

Donde λ* corresponde a un coeficiente de amortiguamiento. En este trabajo se ha considerado un coeficiente tipo Broyden (Broyden, 1965) que surge por minimización de la Norma Euclidiana del vector de funciones de discrepancia. El valor numérico de λ* se restringe a dos condiciones fundamentales:

La propuesta de la inclusión del parámetro λ* tipo Broyden permite un control permanente del camino de convergencia y el cumplimiento de la restricción de no negatividad (flash negativo), además de tener una funcionalidad matemática directa con las ecuaciones residuales R-RG (a diferencia del parámetro de amortiguamiento propuesto por Leibovici y Neoschil (1995), el cual no posee ninguna relación matemática con estas ecuaciones), reflejándose en alta estabilidad numérica y convergencia efectiva a la solución del sistema. Para lograr tal estabilidad fue necesario expresar los elementos de la matriz Jacobiana de forma analítica, así:

Para acelerar la convergencia de la solución del sistema y ayudar a la estabilidad de método de Newton-Raphson, se utilizó un algoritmo de reducción matricial de Gauss-Jordan con pivoteo parcial y re-escalado de columna.

RESULTADOS Y DISCUSIÓN

El objetivo de ésta solución fue obtener el conjunto de valores de {Ψ_k} que anulen las funciones residuales (ec. 5) con tolerancias numéricas r menores a 10^-10. Para evaluar la efectividad de la solución propuesta y su capacidad de convergencia se consideraron sistemas hipotéticos con 2 ≤ P (fases) ≤ 3 y 2 ≤ N (componentes) ≤ 10; y 3≤ P ≤ 8 con 10 ≤ N ≤ 50.

Tanto las composiciones {z_i} del alimento y { χ_k,j} de las P fases, fueron generadas aleatoriamente de forma similar a lo propuesto por Leibovici y Neoschil (1995). A diferencia de estos autores, en este trabajo, los vectores de valor inicial para las fracciones de fase y la composición de referencia fueron también generados aleatoriamente por invocación de una subrutina numérica que consideró el tiempo del sistema (en segundos) como parámetro de arranque (para superar la limitación de repetición de grupos de números). Lo anterior garantiza que la aplicación del método sea generalizable a cualquier sistema de equilibrio de fases, independiente del valor inicial.

Debido a la naturaleza hipotética de los sistemas evaluados, los coeficientes de distribución fueron inicialmente generados a partir de los valores de las fases utilizando la ecuación (3) y se mantuvieron constantes a través del procedimiento iterativo. Esta es una ventaja relevante en la generalización del método propuesto, dado que al ser los sistemas aleatoriamente generados y obtener convergencia absoluta en cada uno de ellos; si se cuenta con sistemas reales con funcionalidad termodinámica, se puede garantizar la capacidad y efectividad de la convergencia de la solución (como un caso particular de los sistemas aleatorios). En el caso de sistemas reales, el método propuesto puede implementarse como una subrutina numérica en algoritmos robustos de descripción de equilibrio de fases, donde se corrigen en ciclos externos los coeficientes de distribución con cálculos de fugacidad y/o actividad.

De otra parte, la inclusión del parámetro de Broyden como factor controlante del tamaño de paso en el método NR, garantizó una inspección efectiva de la convergencia del sistema y del cumplimiento de la restricción general de no negatividad de las fases (Ec. 6). Lo anterior es consecuencia de la fuerte funcionalidad del parámetro λ* con la Norma Euclidiana, que a su vez depende directamente de las funciones de discrepancia. Este resultado es de alta relevancia ya que en investigaciones previas de Leibovici y Neoschil (1995), el parámetro de control de la restricción del flash negativo, no posee ninguna relación matemática con la ecuación R-RG. En la Figura 1, se presenta el comportamiento de la Norma Euclidiana en función del número de iteraciones.

Fig. 1: Comportamiento de la Norma Euclidiana en función del número de iteraciones para sistemas hipotéticos con 4 ≤ P ≤ 8 fases y N = 50 componentes.

Este comportamiento muestra la efectividad de la convergencia de la solución propuesta aplicada en cualquier sistema hipotético con vectores de valor inicial aleatorios. La primera observación substancial es que no existe una tendencia predecible del número de iteraciones en función del número de fases, lo cual es consistente debido a la naturaleza aleatoria e hipotética de los sistemas, además de las características altamente no lineales de las funciones de discrepancia.

Por su parte, la velocidad de convergencia fue afectada de forma crítica como consecuencia de los vectores aleatorios de valor inicial. Si su valor no se encuentra dentro de la región factible, el número de iteraciones se incrementa, pero siempre se obtiene una solución. Este resultado se corrobora en las Tabla 1 y Tabla 2, las cuales muestran el número de iteraciones para simulaciones de varios sistemas aleatorios.

Tabla 1. Número de iteraciones en función del número de fases y componentes.

Tabla 2. Número de iteraciones para P = 2 en función del número componentes.

El resultado más importante, de este análisis es que, a pesar de no existir una relación directa entre el número de fases y componentes, con el número de iteraciones, la Norma Euclidiana tiende a cero, lo que garantiza que la convergencia es absolutamente efectiva. El caso particular de P = 2 fases y N componentes, ha sido el de mayor interés debido a las aplicaciones inmediatas en procesos de separación (Leibovici y Nichita, 2008; Juanes, 2008; Nichita et al, 2006; Vaca et al, 2006; Tiscarreño et al, 1998; Michelsen, 1994). Todas estas investigaciones se han centrado en evaluar la restricción de la ecuación de Rachford-Rice, con el enfoque de minimización de la energía de Gibbs para el flash negativo (Michelsen, 1989) y la estabilidad y convergencia a través de la búsqueda del intervalo óptimo de factibilidad en función del coeficiente distribución. Lo anterior implica que los vectores de valor inicial van a depender de la eficiencia de búsqueda del intervalo de convergencia, de lo contrario no es posible obtener una solución (éste enfoque es válido solo para P = 2). El método propuesto fue probado para sistemas hipotéticos desde 2 hasta P fases, generando aleatoriamente los vectores de valor inicial tanto de las fracciones de cada fase como su correspondiente composición. Este resultado muestra la no vulnerabilidad del método en el proceso iterativo por efecto del valor inicial, independiente del número de fases Ψ y componentes N (la única implicación es el número de iteraciones), sin necesidad de otro tipo de análisis matemático adicional para el problema del valor inicial, tal como se reporta convencionalmente para 2 y 3 fases.

Otro de los resultados importantes de esta investigación es la generalización del método propuesto, dado que, si para un sistema hipotético con un número inicial de P fases independiente de su cantidad (2, 3,..., P) y su número de componentes (2, 3,..., N) existe una solución factible; cuando se evalúe un caso particular definido como un sistema real termodinámicamente descrito, la efectividad del método propuesto va ha ser igualmente válida.

De otra parte, todos los sistemas simulados presentaron estabilidad numérica indicando robustez del algoritmo desarrollado, lo cual es consecuencia de evaluar analíticamente el Jacobiano de la ecuación R-RG y de utilizar un método de Gauss-Jordan con pivoteo parcial y re-escalado de columna para aceleración del método de NR. Por lo tanto, el método propuesto se centra en la contribución acoplada de tres métodos, bien conocidos, en la solución efectiva de las ecuaciones flash de Rachford-Rice, los cuales de forma individual convencionalmente requieren de análisis matemáticos adicionales.

Finalmente, para todos los casos simulados se cumplió la validez de la ecuación (6); en la mayoría de ellos los valores obtenidos para las fases fueron positivos y menores a la unidad, indicando convergencia sobre la zona factible de separación en equilibrio, no obstante, se obtuvieron también, valores de las fracciones de fase negativas ó mayores a la unidad (tanto en el proceso iterativo como en la solución final, mostrando aparición y desaparición de fases), pero todos los casos, fueron consecuentes con la naturaleza positiva de las fracciones de los componentes, es decir, consistentes con el flash negativo (ec. 6) y con la definición matemática de la ecuación R-RG.

CONCLUSIONES

Se planteó una solución numérica para la ecuación de Rachford-Rice generalizada, acoplando el método de Newton-Raphson modificado con un parámetro de amortiguamiento tipo Broyden y el concepto de flash negativo. La solución fue probada con sistemas hipotéticos de P fases y N componentes, generando aleatoriamente los vectores de valor inicial.

Se encontró que la solución propuesta es altamente estable y convergente para cualquier tipo de vector de inicio y que el número de iteraciones es afectado críticamente por éste vector y no necesariamente por la cantidad de fases o número de componentes.

La introducción del parámetro de Broyden garantizó un control permanente en la Norma Euclidiana y de la convergencia del sistema, además del cumplimiento de la condición de no negatividad de las composiciones de fase. El método obtenido puede ser aplicable y generalizable a cualquier número de fases y componentes.

AGRADECIMIENTOS

Los autores agradecen a la Escuela de Ingeniería Química y a la Vicerrectoria de Investigaciones de la Universidad del Valle por el apoyo financiero mediante el proyecto "Modelamiento y Simulación de la Cinética Intrínseca de Reacciones Fotocatalíticas Heterogéneas" (código interno 2520) y a Colciencias (Programa de Doctorados Nacionales) por el apoyo en la realización de estudios doctorales.

REFERENCIAS

Broyden C., A Class of Methods for Solving Nonlinear Simultaneous Equations, Mathematical of Computations, 19 (92), 577-593 (1965).

Henley, E. y J. Seader, Operaciones de Separación por Etapas de Equilibrio en Ingeniería Química, 611-621. ISBN 84-291-7908-9, Reverté, Barcelona, España (1988).

Juanes, R., A Robust Negative Flash Base don a Parameterization of the Tie-Line Field, Fluid Phase Equilibria, 267, 6-17 (2008).

Leibovici C. y J. Neoschil, A Solution of Rachford-Rice Equations for Multiphase Systems, Fluid Phase Equilibria, 112 217-221 (1995).

Leibovici, C. y D. V. Nichita, A New Look at Multiphase Rachford-Rice Equation for Flashes Negatives, Fluid Phase Equilibria, 267 127-132 (2008).

Michelsen, M., Calculation of Multiphase Equilibrium, Computer and Chemical Engineering, 18 (7), 545-550 (1994).

Nelson, P.A., Rapid Phase Determination in Multiple-Phase Flash Calculations, Computers & Chemical Engineering, 11 (6), 581-591 (1987).

Nichita, D. V., Broseta, D. y C. Leibovici, Consistent Delumping of Multiphase Flash Results, Computer and Chemical Engineering, 30, 1026-1037 (2006).

Tiscarreño, F., Gómez, A., Jiménez, A. y R. Chávez, Multiplicity of the Solutions of the Flash Equations, Chemical Engineering Science, 53 (4), 671-677 (1998).

Vaca, M., Jiménez, A. y R. Monroy-Loperena, On the Multiple Solutions of the Flash Equations, Chemical Engineering Science, 61, 3850-3857 (2006).

Whitson C. y M. Michelsen, The Negative Flash, Fluid Phase Equilibria, 53, 57-71 (1989).