Aplicabilidad De La Mezcla Cloroformo-acetona, Equilibro De Fases Con Ecuaciones De Estado

DIANA MARITZA VASQUEZ MOLINA FRANKY ESTEBAN BEDOYA LORA

SISTEMA A TRABAJAR Y APLICABILIDAD

El sistema a trabajar es la solución Cloroformo-Acetona, adjunto al mensaje se envía la tabla que contiene todos los datos del sistema, la fuente bibliográfica fue: “Vapor liquid equilibrium data collection” de la serie DECHEMA.

APLICABILIDAD DE LA MEZCLA CLOROFORMO-ACETONA

La mezcla cloroformo-acetona tiene aplicabilidad principalmente a nivel de laboratorio según lo consultado, se encontró usos muy específicos en el caso de prácticas de laboratorio en particular para la determinación de impurezas mecánicas a una muestra de propóleos y en el siguiente procedimiento para la determinación de compuestos fenólicos. En el área de investigación esta mezcla es bastante utilizada en análisis cromatográfico empleando el método de capa fina (TLC) de los cuales se encontraron ejemplos como la tipificacion de aflatoxinas (toxinas producidas por un moho que crece en las nueces, en semillas y en las legumbres) mediante este procedimiento utilizando un sistema de solventes cloroformo-acetona, también se utiliza en la separación de antibióticos específicamente para el Pyoluteorin y el Gliotoxin, en estos procedimientos el solvente debe presentar cierta composición de cloroformo y acetona para efectuar una separación más definida, es por esto que el equilibrio líquido-vapor de esta solución es importante ya que en muchos casos se desea mantener cierta concentración de las especies en el estado liquido, no solo mientras se realiza el proceso de separación por TLC si no también en su forma de almacenamiento y preparación, y esto depende de las condiciones a las cuales se encuentre la solución, como la presión y temperatura de la habitación en la cual se realiza en análisis. Adjunto al mensaje se envía una tabla y uno de los artículos como referencia.

Diana Maritza Vásquez Molina 1036.602.284 Franky Esteban Bedoya Lora 1017.143.241 INFORME: Sistema: Acetona-Cloroformo a 760 mmHg a) Para evaluar los coeficientes de fugacidad de cada componente en la mezcla: 1. Al escoger el modelo de cálculo adecuado, se tuvo en cuenta las propiedades de las sustancias trabajadas, en este caso la polaridad fue una de la más importantes ya que ambos compuestos (Acetona – Cloroformo) tenían un comportamiento polar significante; por esta razón en un principio se utilizó la correlación de Tsonopoulos (la cual incluye el efecto de la polaridad). Los cálculos de los coeficientes de fugacidad de cada componente puro para la primera temperatura fue el siguiente:

T1 = 335.7 K Tc (acetona) = 508 K Pc(acetona) = 4.76Mpa Pr = P / Pc

P = 101.325Mpa Tc(cloroformo) = 536.6 K Pc(cloroformo) = 5.47 Mpa Tr = T / Tc

BPc = f 0 + wf 1 + f 2 (1) RTc 0.330 0.1385 0.0121 0.000607 Donde f 0 = 0.1445 − − − − Tr Tr 2 Tr 3 Tr 8 f 1 = 0.0637 +

0.331 0.423 0.008 − − Tr 2 Tr 3 Tr 8

a a (acetona) = - 0.0309 a(cloroformo) = - 0.04329 Tr 6 Para la acetona hallamos B de la ecuación (1) y calculamos el coeficiente de fugacidad puros así: f2 =

ϕ = exp( BP / RT )

φ acetona = 0.957

Realizando el mismo procedimiento para el cloroformo: φ cloroformo = 0.949

Como se puede deducir de los resultados por este método ninguna de las dos sustancias se puede considerar de comportamiento ideal, por lo tanto es necesario calcular los coeficientes de mezcla, para verificar si es una mezcla ideal; para que el procedimiento fuera menos complejo, se pensó en analizar la desviación de los anteriores resultados con respecto a la ecuación virial, de la cual si teníamos reglas de mezclado.

Los cálculos de los coeficientes de fugacidad de cada componente puro para la primera temperatura usando las correlaciones generalizadas fue el siguiente: BPc = B 0 + WB1 (2) RTc Donde

B 0 = 0.083 −

W acetona = 0.3124

0.422 Tr1.6

B = 0.139 − 1

0.172 Tr 4.2

W cloroformo = 0.219

Reemplazando en (2) y utilizando: ϕ = exp( BP / RT )

ϕ acetona = 0.9684

ϕcloroformo = 0.9694

El porcentaje de error para la acetona fue: 1.2% y para el cloroformo: 2.2 %, se puede concluir que no se presenta mucha diferencia, por consiguiente hubiera sido correcto utilizar las reglas de mezclado para las correlaciones generalizadas de la ecuación virial, sin embargo decidimos utilizar otra correlación mas apropiada para este sistema. 2. Para mayor exactitud en los cálculos se decidió utilizar el modelo de O’connell-Prausnitz, muy adecuado para gases polares. Para los componentes puros:

(1)

Donde:

es la presión critica del componente i es la temperatura critica del componente i

f u ( µ R ,TR ) = −5.237220 + 5.665807( Lnµ R ) − 2.133816( Lnµ R ) 2 + 0.2525373( Lnµ R )3 + 1 / Tr [5.769770 − 6.181427( Lnµ R ) + 2.283270( Lnµ R ) 2 − 0.2649074( Lnµ R )3 3

f a (TR ) = exp[6.6(0.7 − TR )] TR = T / Tc

Tc / acetona) = 508 K

PR = P / Pc

Pc (acetona) = 46.798atm

De tablas del artículo: Whacetona = 0.187

Tc / cloroformo) = 536.6 K Pc (cloroformo) = 53.985atm Whcloroformo = 0.187

µ acetona ( Debye) = 2.88 η acetona = 0

µcloroformo ( Debye) = 1.02 ηcloroformo = 0.28

Despejando de la ecuación (1), Bii y a partir de la ecuación virial:

ϕ = exp( BiiP / RT ) Así se calcularon los coeficientes para los componentes puros, donde el subíndice i se reemplazo por acetona y posteriormente por cloroformo (B11 y B22). Obtuvimos los siguientes resultados

ϕ acetona = 0.9569

ϕ cloroformo = 0.9945

Se puede deducir que el comportamiento del cloroformo no difiere mucho del ideal, pues, su coeficiente de fugacidad con un valor de 0.9945 se acerca a 1, por el contrario para la acetona el coeficiente de fugacidad con un valor de 0.9569 no es gas ideal (estos datos fueron tomados para la primera temperatura). Siguiendo las reglas de mezclado del método escogido tenemos:

µ rij =

105 µi µ j Pcij Tcij

2

De la ecuación (1), para hallar el coeficiente de fugacidad de la mezcla cambiamos el subíndice i por ij TR = T / Tcij

δ ij = 2 * Bij − Bi − B j

δ 12 = 2 * B12 − B1 − B2 Luego la regla de mezclado para la ecuación virial, hallamos los coeficientes de fugacidad para cada compuesto en la mezcla: ^  P  2 ln(ϕ acetona ) =  ( B11 + y 2 δ 12 )  RT  ^  P  2 ln(ϕ cloroformo ) =  ( B22 + y1 δ 12 ) RT  

Despejando

^

 P   2 ( B11 + y 2 δ 12 )  RT    P   2 = exp  ( B22 + y1 δ 12 )  RT  

ϕ acetona = exp  ^

ϕ cloroformo

Los coeficientes de fugacidad de los compuestos en estado de saturación ϕisat se hallaron de forma análoga a los puros y utilizando la correlación de O’connell-Prausnitz para el segundo coeficiente virial (B) pero con una presión igual a la presión de saturación del compuesto a la temperatura dada, esto se hizo recurriendo a la ecuación de Antoine.

b) Para los coeficientes de actividad: ^

^

^

fL fv y Pϕ γ i = iid = i = i sat isat xi f i xi Pi ϕ i fi

Asumiendo f i = Pi sat Donde para el modelo real se convierte en f i = Pi sat ϕ isat ^ v

^ L

Además para el sistema en equilibrio µ iv = µ iL , entonces f i = f i y por lo tanto

DESCRIPCIÓN DEL PROGRAMA EN MATLAB

El programa en matlab realiza el procedimiento descrito anteriormente usando el modelo de O’connell- Prausnitz, calculando los coeficientes de fugacidad de los componentes puros y los de la mezcla, para las diferentes temperaturas. Se denota con un subíndice a, cuando nos referimos a la acetona, con el subíndice c cuando nos referimos al cloroformo y m para la mezcla. Fi coeficiente de fugacidad, r coeficientes de actividad. BPc para el cloroformo RTc BP pbrta = c para la acetona RTc BP pbrtm = c para la mezcla RTc Las constantes de Antoine para calcular la presión de saturación fueron tomadas de la serie DECHEMA en la página de nuestro sistema de mezcla pbrtc =

ANALISIS DE RESULTADOS En la primera gráfica vemos como el coeficiente de fugacidad disminuye para la acetona cuando la fracción de vapor aumenta, mientras que para el cloroformo aumenta cuando se incrementa la fracción de vapor en la mezcla. Cuando la fracción de vapor tiende a 0, se puede notar que el coeficiente de fugacidad del cloroformo en la mezcla se acerca a la curva de coeficientes de fugacidad para el cloroformo puro y cuando la fracción de vapor tiende a 1 la curva de coeficientes de fugacidad de la acetona se acerca a la de los coeficientes de fugacidad para la acetona pura, por conceptos debemos tener en cuenta que cuando lo anterior sucede es que se está acercando al comportamiento de una mezcla ideal, es decir, cuando el compuesto esta interactuando cada vez mas consigo mismo. También se puede apreciar como los valores de los coeficientes de fugacidad para el cloroformo a diferentes fracciones de vapor, se mantienen casi constantes, con un valor muy aproximado a 1, lo que nos indica que no sería un error tomarlo como gas ideal en este caso; para la acetona como componente puro la función de los coeficientes de fugacidad va aumentando levemente a medida que la fracción de vapor crece, pero ninguno de estos valores es suficientemente cercano a 1 para poder considerarlo como gas ideal.

Para el siguiente grupo de gráficas: 1. La grafica de la temperatura vs la fracción líquida muestra una variación directa hasta una fracción de 0.35 aproximadamente e indirecta en las siguientes, presentándose un máximo de temperatura en el azeótropo.

2. Donde se observa el intercepto entre las dos gráficas se identifica el azeótropo y también se presenta un máximo de temperatura. 3. En esta gráfica al compararla con una línea de 45º donde los valores son iguales, se analiza como en el azeótropo las composiciones tanto en la fase líquida como en la fase vapor son iguales. 4. Analizando el coeficiente de actividad como medida del grado de divergencia del comportamiento de la sustancia con respecto al ideal se puede deducir como mientras la fracción del liquido de la acetona aumenta, el coeficiente de actividad se acerca a 1 desde abajo para la acetona mientras que el del cloroformo disminuye, entre más cercano a 1 sea el coeficiente de actividad, menor es la desviación del comportamiento ideal de mezcla, por lo tanto esta gráfica nos dice que para el cloroformo con fracción liquida tendiente a 1 (X21), se puede asumir como gas ideal en una mezcla ideal y para la acetona con fracción liquida tendiente a 1 (X11) se podría asumir como mezcla ideal solamente, porque su coeficiente de fugacidad puro es un poco alejado de 1. El modelo de la gráfica es simétrico. BIBLIOGRAFIA: AICHE Journal Chemical Engineering Research and Development, Marzo 1974, “An empirical Correlation of Second Virial Coefficients” Constantine Tsonopoulos pág 263. I&EC Process Desing and Development, Vol 6 No 2 Abril 1967, “Empirical correlation of second virial coefficients for Vapor-Liquid equilibrium calculations”, J. P. O’Connell, y J.M. Prausnitz, pág 245.

Diana Maritza Vásquez Molina 1036.602.284 Franky Esteban Bedoya Lora 1017.143.241 INFORME 2º ENTREGA Sistema: Acetona-Cloroformo a 760 mmHg CALCULO DE COEFICIENTES b12 b21 Y α PARA EL MODELO NRTL De la entrega anterior tomamos los coeficientes de actividad calculados con la ecuación virial y correlaciones para el segundo coeficiente. Son los siguientes: T (ºC) 62.7 63.55 63.95 64.4 64.55 64.25 63.55 62.3 61.8 61.27 60.8 60.1 59.3 57.95 57.1 56.7

X1 0.1013 0.1792 0.2585 0.3022 0.3697 0.4418 0.5268 0.6318 0.6683 0.7020 0.7315 0.7605 0.8137 0.8946 0.9433 0.9652

γ1 (experimentales) 0.6078 0.6434 0.6825 0.7277 0.7808 0.8286 0.8837 0.9214 0.9411 0.9499 0.9567 0.9698 0.9755 0.9911 0.9955 0.9980

γ2(experimentales) 0.9820 0.9690 0.9619 0.9316 0.8980 0.8674 0.8191 0.7796 0.7441 0.7328 0.7226 0.7182 0.7061 0.6394 0.6261 0.6119

Calculamos los coeficientes b12, b21 y α mediante dos procedimientos: Iteración normal con minimización de Energia de gibbs en exceso y por el método de Barker. METODO NORMAL DE ITERACION (ENERGIA DE GIBBS EN EXCESO) Utilizamos el programa MATLAB 6.5 con el toolbox para optimización, ya que este contenía las funciones adecuadas para realizar iteraciones. Para la constante de gases ideales R escogimos 1.987 Cal/(K*mol), estas dimensiones satisfacen las unidades de b12 y b21 encontradas en el DECHEMA, además la iteración converge mas rápidamente con esta constante. Por lo tanto b12 y b21 toman las unidades de mol/cal. Se crea un vector z, que contiene las variables a iterar, donde z(1)=b12, z(2)=b21 y z(3)=α, se crea una funci ón “minimo” la cual se optimiza con la herramienta “lsqnonlin”. “minimo” tiene la siguiente forma:

 G E  G E     = minimo −    RT   Experimental  RT  Teorico  Donde  G E   = x1 * ln(γ 1* ) + x 2 * ln(γ 2* )  RT   Experimental 

(1)

(2)

γ 1* y γ 2* son los coeficientes de actividad experimentales  G E   = x1 * ln(γ 1 ) + x 2 * ln(γ 2 )  RT   teorico 

(3)

Donde 2     G21 G12 * τ 12  +  ln(γ 1 ) = x * τ 21  (x2 + x1 * G12 )2    x1 + x 2 * G21  2     G12 G21 * τ 21 2  +  ln(γ 2 ) = x1 * τ 12  (x1 + x2 * G21 )2    x 2 + x1 * G12  2 2

G12 =exp(-α* τ 12 ) G21 =exp(-α* τ 21 ) b b τ 21 = 21 τ 12 = 12 RT RT

(4)

(5)

(6) (7)

Esta función se minimizó con el método Gillmurray, partiendo de valores b12=0, b21=0 y α=0.3, con rangos de -10000 a 10000 para b12 y b21 y de 0.2 a 0.47 para α. La herramienta “lsqnonlin” nos retorna los valores de ‘z’ y una variable llamada ‘fval’ que contiene el numero mínimo alcanzado por la función ‘minimo’, que en nuestro caso da del orden 10e-4, lo cual es aceptable. Con los coeficientes hallados calculamos de nuevo los valores de γ 1 y γ 2 , estos serán los teóricos (‘rteorica1’ y ‘rteorica2’), y calculamos el porcentaje de error con respecto a los experimentales. 2       G21 G12 * τ 12 2  +  (8) γ 1 = exp x 2 * τ 21  2  * x x G + ( ) x x G * +   2 21   2 1 12     1 2       G12 G21 * τ 21 2  +  (9) γ 2 = exp x1 * τ 12  (x1 + x2 * G21 )2     x 2 + x1 * G12   A continuación graficamos los coeficientes de actividad teóricos γ 1 y γ 2 con * * respecto a x1, y en la misma grafica los γ 1 y γ 2 experimentales.

METODO DE BARKER PARA ITERACION (PRESIONES) Para este método se procedió de la misma forma que el descrito anteriormente, solo que la función a minimizar es diferente, las demás variables se nombran igual pero con una ‘b’ al final. Fue necesario calcular las Presiones de saturación para los compuestos 1 y 2 utilizando las constantes de Antoine, inicialmente dadas en mmHg fueron llevadas a atmósferas. La presión total del sistema es 1 atm. La ecuación a minimizar fue la siguiente: Pteorica – Pexperimental = minimo

(10)

Donde : Pteorica = Psistema = 1 atm

(11)

Pexperimental = x1 * P1sat * γ 1 + x 2 * P2sat * γ 2

(12)

γ 1 y γ 2 siguen las ecuaciones (8) y (9) del método anterior que contienen las

variables a iterar b12, b21 y α. Se procede de igual forma para la iteración, solo que en este método tomamos como valor de referencia para los coeficientes b12, b21 y α (‘zB’) los mismos valores hallados en el anterior procedimiento

(‘z’). Calculamos γ 1 y γ 2 teóricos, y los porcentajes de error con respecto al experimental. Graficamos γ 1 y γ 2 teóricos y experimentales versus x1.

ANALISIS DE RESULTADOS Analizando la dos gráficas (Método Normal y Método Barker) tenemos que los valores de γ 1 se acercan bastante a los valores de γ 1 esto lo podemos comprobar gráficamente porque las curvas tienen puntos que están demasiado cerca, también lo podemos verificar analizando los porcentajes de error que son muy pequeños y van disminuyendo conforme la fracción líquida del * componente 1 aumenta, en comparación con las curvas de γ 2 y γ 2 las cuales *

aunque son muy cercanas se desvían en mayor proporción que las de γ 1 , y se sigue el mismo patrón o comportamiento, mientras la x2 aumenta (x1 disminuye) mas exacta es la predicción. Este comportamiento se debe a que los coeficientes de actividad experimentales son menos predecibles conforme la concentración o fracción del componente disminuye. *

Los valores hallados de b12, b21 y α son diferentes en los dos étodos m en comparación con los dados en el DECHEMA, esto porque la iteración se puede realizar desde diferentes puntos iniciales arrojando convergencias distintas, además estas poseen infinitas soluciones. Como podemos ver en los porcentajes de error el método de Barker es más inexacto y presenta mayor desviación con respecto a los datos experimentales, esta observación se ve claramente al comparar las dos graficas, además los

porcentajes de error máximo para Barker están alrededor de 12% en cambio para el método normal el mayor es de 7%. Estas diferencias se deben a que el método Barker utiliza para su iteración valores de Presión de saturación los cuales son hallados utilizando las constantes de antoine que son aproximaciones y por lo tanto contribuyen significativamente al error. BIBLIOGRAFIA SMITH, VAN NESS Y ABBOT; “Introducción a la Termodinámica en la Ingeniería Química”, 4 edición, editorial McGRAW-HILL RIVERA, Martín “Optimización con MATLAB“, Universidad Iberoamericana Ciudad de México.

Diana Maritza Vásquez Molina 1036.602.284 Franky Esteban Bedoya Lora 1017.143.241 INFORME 3º ENTREGA Sistema: Acetona-Cloroformo a 760 mmHg CALCULO DE COEFICIENTES DE ACTIVIDAD POR EL METODO UNIFAC Tomando los datos de temperatura y fracción del componente 1 (acetona) en la fase liquida, calculamos los coeficientes utilizando las ecuaciones del modelo UNIFAC tomadas del libro “Introducción a la Termodinámica en Ingeniería Química” T (ºC) 62.7 63.55 63.95 64.4 64.55 64.25 63.55 62.3 61.8 61.27 60.8 60.1 59.3 57.95 57.1 56.7

X1 0.1013 0.1792 0.2585 0.3022 0.3697 0.4418 0.5268 0.6318 0.6683 0.7020 0.7315 0.7605 0.8137 0.8946 0.9433 0.9652

La ecuaciones fueron las siguientes ln γ i = ln γ iC + ln γ iR  J J  ln γ iC = 1 − J i + ln J i − 5q1 − i + ln i  Li   Li  S S  ln γ iR = qi (1 − ln Li ) − ∑ θ k ki − Gki ln ki  ηk ηk  k  ri Ji = ∑ rj x j j

Li =

qi ∑qjxj

ri = ∑ vk(i ) Rk k

qi = ∑ vk(i )Qk

Gki = vk(i )Qk

θ k = ∑ Gki xi i

k

j

ski = ∑ Gmiτ mk m

η k = ∑ ski xi τ mk = exp i

− amk T

El subíndice i identifica los componentes (1 y 2) y j es el subíndice de unión que opera sobre y todos los componentes. El subíndice k identifica los subgrupos y m es un índice de unión que opera sobre todos los subgrupos. Los valores de los parámetros Rk y Qk y de los parámetros de interacción de grupos amk los leímos del libro “The Properties of Gases and Liquids” de Prausnitz, John M. Los subgrupos derivados de las sustancias fueron los siguientes: Acetona CH3COCH3 ===> CH3CO (1) y CH3 (2) Cloroformo CHCl3 #

Formula

1 2 3

CH3CO CH3 CHCl3

===> CHCl3 (3)

k 18 1 50

Rk 1.6724 0.9011 2.8700

Qk 1.488 0.848 2.410

vk(1) 1 1 0

vk( 2 ) 0 0 1

Parámetros de interacción entre los grupos k 18 1 50

18 0 476.40 552.10

1 26.76 0 36.70

50 -354.60 24.90 0

ANALISIS Como podemos ver en los resultados y comparando con los coeficientes de actividad experimentales hay un alto porcentaje de error para aquellos datos que tienen una fracción menor a 0.5, por ejemplo para la sustancia 1 los errores son apreciables para las composiciones de 1 pequeñas, igual pasa con la sustancia 2, esto se debe en gran parte a que determinar la concentración de una sustancia de composición pequeña es muy difícil experimentalmente debido a la inexactitud técnica. Sin embargo para concentraciones relativamente altas la predicción es bastante exacta y satisfactoria.

Diana Maritza Vásquez Molina 1036.602.284 Franky Esteban Bedoya Lora 1017.143.241 INFORME 4º ENTREGA CALCULO DE LA TEMPERATURA DE BURBUJA Se tomaron los mismos 16 datos trabajados durante todo el proyecto, en este caso solo las fracciones líquidas X1(Real) 0.1013 0.1792 0.2585 0.3022 0.3697 0.4418 0.5268 0.6318 0.6683 0.7020 0.7315 0.7605 0.8137 0.8946 0.9433 0.9652 Los coeficientes a utilizar en la ecuación NRTL fueron los siguientes (Obtenidos del DECHEMA): A12= -481.7574 A21=106.6503 α=0.3030 Nos basamos en el siguiente algoritmo para hallar las respectivas Temperaturas de burbuja y fracciones en el vapor.

Ecuaciones (1) TkSat =

Bk − Ck Ak − log 10( P)

(2) T = ∑ x x TkSat k

  Bk  Antoine (3) PkSat = 10^  Ak − T (º C ) + C k   (4) Los coeficientes de actividad se calcularon con la ecuación NRTL asi: 2       G21 G12 * τ 12 2  +  γ 1 = exp x 2 * τ 21  (x2 + x1 * G12 )2     x1 + x 2 * G21   2       G12 G21 * τ 21 2  +  γ 2 = exp x1 * τ 12  (x1 + x2 * G21 )2     x 2 + x1 * G12   G12 =exp(-α* τ 12 ) G21 =exp(-α* τ 21 ) A A τ 12 = 12 τ 21 = 21 RT RT

(5) PjSat =

(6) T =

P x γ ∑k  Φk k  k Bj

 PkSat   P Sat  j

A j − log 10( PjSat )

   

−Cj

x k γ k PkSat Φk P (8) Los coeficientes de fugacidad los hallamos utilizando la correlación O'conellPrausnitz para el segundo de la ecuación virial. Este método esta descrito en el informe 1 de este proyecto. (7) y k =

CALCULO DE LA TEMPERATURA DE ROCIO Se tomaron 16 datos que corresponden a los mismos puntos en el vapor dados en el DECHEMA Y1 (Real) 0.0740 0.1428 0.2221 0.2814 0.3724 0.4695 0.5862 0.7070 0.7526 0.7852 0.8123 0.8376 0.8793 0.9411 0.9699 0.9822 Nos basamos en el siguiente algoritmo para hallar las respectivas Temperaturas de rocio y fracciones en el liquido (9) T = ∑ y x TkSat k

(10) (11)

(12)

y Φ PjSat = P ∑  k k k  γk y Φk P x k = k Sat γ k Pk xk =

xk ∑ xk

Sat  Pj  Sat   Pk

   

ANALISIS DE RESULTADOS Analizando las temperaturas de burbuja observamos que estás presentan un máximo en su punto azeotrópico, esto se debe a que los coeficientes de actividad calculados anteriormente eran menores a uno, por tal motivo los componentes prefieren estar juntos (desviaciones negativas), esto se manifiesta como un máximo en la temperatura de ebullición ya que los compuestos como prefieren estar juntos en su fase líquida presentan mayor resistencia a formar vapor, este comportamiento para esta mezcla en particular se debe a que los cloros del cloroformo forman puentes de hidrógeno con la acetona produciendo fuertes enlaces de unión y presentando una mejor afinidad en mezcla que estando puros. Al igual que en las entregas anteriores podemos ver como el mayor porcentaje de error (6.7%) se presenta cuando las composiciones son pequeñas, asumiríamos que se debe a un error de lectura experimental, ya que las composiciones pequeñas son difíciles de leer. El modelo NRTL es muy exacto, las curvas de equilibrio son demasiado parecidas y solo presentan porcentajes de error alrededor de 1% para composiciones moderadas. Para la temperatura de rocío hallada, la cual es un poco mayor a la temperatura de burbuja, también se presentan desviaciones bajas para composiciones liquidas moderadas, cuando estas son pequeñas las desviaciones crecen (5.9% para la mayor desviación) siguiendo el mismo comportamiento para las composiciones en la fase vapor halladas para el punto de burbuja. En el punto azeotrópico ambas curvas, tanto la de burbuja como la de rocío alcanzan su máximo y se intersecan, esto porque al formarse el azeótropo las composiciones de la fase vapor son iguales a la fase liquida y este valor es único para ambas. Las curva de equilibrio se trazó con la composiciones en el vapor halladas según el modelo NRTL y con base en las composiciones en el líquido reales, estas composiciones (Y1) fueron las obtenidas por el procedimiento para el cálculo de las temperaturas de burbuja; podemos apreciar que esta curva de equilibrio no difiere mucho de la experimental, tal y como vimos en la 2 entrega, el modelo NRTL es muy funcional para este sistema Acetona-Cloroformo y es aún más fiel que el modelo UNIFAC.

BIBLIOGRAFIA SMITH, VAN NESS Y ABBOT; “Introducción a la Termodinámica en la Ingeniería Química”, 4 edición, editorial McGRAW-HILL

I&EC Process Desing and Development, Vol 6 No 2 Abril 1967, “Empirical correlation of second virial coefficients for Vapor-Liquid equilibrium calculations”, J. P. O’Connell, y J.M. Prausnitz, pág 245.

PROYECTO FINAL (S06-1) - EQUILIBRIO DE FASES SISTEMA ACETONA-CLOROFORMO

DIANA MARITZA VASQUEZ MOLINA 1036.602.284 FRANKY ESTEBAN BEDOYA LORA 1017.143.241

TERMODINAMICA II

PROFESOR FELIPE BUSTAMANTE

UNIVERSIDAD DE ANTIOQUIA INGENIERIA QUIMICA FACULTAD DE INGENIERIA MEDELLIN 2006

INFORME FINAL PROYECTO TERMODINAMICA II

1. INTRODUCCION

En el presente proyecto se trabajó el equilibrio liquido – vapor para la mezcla Acetona – Cloroformo, desde ahora (1) y (2) respectivamente, con el fin de establecer un modelo eficaz que evalúe sus propiedades a diferentes temperaturas y presión constante, además se utilizo el software PRO/II para comparar los datos obtenidos con los experimentales y los calculados, los algoritmos fueron realizados utilizando MATLAB. El proyecto se divide en 4 partes a saber:

1. Evaluación de los coeficientes de fugacidad y cálculo de coeficientes de actividad experimentales. 2. Calculo de los coeficientes para el modelo NRTL y hallar coeficientes de actividad con estos parámetros. 3. Calculo de coeficientes de actividad según el modelo UNIFAC 4. Determinación de las temperaturas de rocío y burbuja utilizando los parámetros para el modelo NRTL dados en el DECHEMA.

El presente informe final, contiene el análisis y comparación de los distintos modelos usados, y la apreciación del mejor modelo que describa el comportamiento de los datos experimentales. Además se presenta una breve descripción de los procesos llevados a cabo para realizar cada parte del proyecto, enfocándonos en las dificultades presentadas.

2. DESCRIPCION Y DIFICULTADES

Para realizar la primera entrega, (Ver archivo “LiqVap.m”)se nos pidió a partir de los datos experimentales, evaluar el coeficiente de fugacidad de cada componente en la fase vapor, y de cada componente como vapor saturado, para realizarlo se pensó primero en utilizar la ecuación de estado de virial con sus reglas de mezclado, para facilitar un poco los cálculos, pero en nuestro caso, para las sustancias acetona y cloroformo, era necesario tener en cuenta las propiedades moleculares debido a que ambas sustancias poseen un carácter polar que no se puede despreciar fácilmente.

Para tomar en cuenta este efecto, se decidió utilizar la ecuación de Tsonopoulos, pero después se encontró el artículo indicado que abarcaba las ecuaciones necesarias para realizar los cálculos de coeficientes de fugacidad para los componentes puros los cuales eran necesarios compararlos con un valor de 1 el cual identificaba el comportamiento ideal; como resultado obtenido, la diferencia no era muy grande, pero tampoco era lo suficientemente pequeña como para aproximar las sustancias a este comportamiento ( se puede hacer una salvedad en el caso del cloroformo que dio un valor de 0.9945 que a diferencia de la acetona tiende a tener el ϕ cloroformo = 1 ), por lo que posteriormente se debió calcular los coeficientes de cada componente, ahora, en la mezcla, las reglas de mezclado fueron obtenidas de la correlación de O´conell-Prausnitz para el segundo coeficiente virial (B), este procedimiento se hizo para analizar si la mezcla era o no ideal, se llegaba a que era ideal si los resultados de los valores de los coeficientes de fugacidad para cada componente puro era igual al valor del coeficiente de fugacidad pero en la mezcla. La correlación de O´conell-Prausnitz era adecuada para los componentes trabajados, incluso en el artículo se obtuvieron valores tabulados para estas sustancias. Las reglas de mezclado para esta correlación no diferían mucho con respecto a las utilizadas por Tsonopoulos, sin embargo eran un poco más complejas.

Para realizar los anteriores cálculos solo fue necesario buscar de tablas datos como: Pc, Tc y w tabulados para cada sustancia respectivamente.

Los valores de los coeficientes de fugacidad de cada sustancia tanto pura como en mezcla era necesarios para conocer los valores de los coeficientes de actividad ^ v

^ L

con la fórmula correspondiente que parte del equilibrio f i = f i .

Para la segunda entrega (ver archivo “nrtl.m”) calculamos los parámetros de la ecuación NRTL, siguiendo dos métodos de iteración el primero era el método normal (energía de Gibbs en exceso) y el segundo el método de Barker, donde lo que se deseaba encontrar eran los valores de b12, b21 y α. Como primera dificultad que presentó fue encontrar la función correcta para realizar el programa en Matlab como solución se buscó en Internet y se encontró ejemplos parecidos que minimizaban funciones en Matlab, haciendo una analogía con nuestro problema, se escogió la herramienta “lsqnonlin”. También se encontró un método de iteración en el cual nos basamos: Gillmurray. Estos métodos convergían en valores un poco diferentes a los registrados en el DECHEMA, sin embargo, se sabía que existen infinitas soluciones para estos coeficientes.

Para la tercera entrega (ver archivo “unifac.m”), se debió hacer un programa en Matlab que calculara los coeficientes de actividad por el método UNIFAC.

El método UNIFAC tiene ecuaciones bastante complejas. Lo primero que se hizo fue escoger los grupos en los cuales se iban a dividir nuestras sustancias, para analizar su contribución, fue sencillo en nuestro caso pues el cloroformo, por ejemplo, constituía un solo grupo y en el caso de la acetona fueron solo dos grupos, metil y metilcetona:

Acetona CH3COCH3 ===>

CH3CO (1) y CH3 (2)

Cloroformo CHCl3

CHCl3 (3)

===>

En esta entrega nos enfrentamos a una gran cantidad de datos que debían ser procesados, para esto elegimos un método basado en vectores y de esta forma podíamos realizar una gran cantidad de operaciones con estos, el esquema fue el siguiente VARIABLE(SUSTANCIA,SUBGRUPO,DATO), así, si necesitábamos referirnos a las sustancias simplemente cambiábamos la primera coordenada del vector y así mismo con los subgrupos y los datos que correspondían a los 16 puntos tomados del DECHEMA.

En la última entrega (ver archivo “burbuja.m”, la cual fue relativamente más fácil que las demás, debido a que se tenía el algoritmo y las formulas que se debían utilizar (SMITH, VAN NESS Y ABBOT; “Introducción a la Termodinámica en la Ingeniería Química”, 4 edición ) solo se hizo necesario la función while para realizar las iteraciones; además los métodos para hallar los coeficientes de fugacidad y de actividad estaban dados en otros programas de entregas anteriores, así que tan solo se tuvo que realizar algunos cambios de variable y copiarlos en el nuevo programa. Al realizar la grafica de las temperaturas de rocío burbuja ambas en función de la composición del componente 1 liquido (X1) obtuvimos la misma curva, sin embargo se tenía un error y era que la Temperatura de roció se grafica con respecto a Y1 y no con X1.

3. ANALÍSIS DE RESULTADOS

Nota: Para observar los resultados ejecutar el programa “resultados.m” en MATLAB

La figura

(1) contiene en la primera fila los coeficientes de actividad para la

sustancia 1 y 2 calculados con el modelo UNIFAC, el primero hallado con el programa realizado en Matlab y el segundo basado en los datos aportados por PRO/II, de las graficas se puede decir que son exactamente las mismas, presentan las mismas desviaciones respecto a la experimental, por lo que

concluimos que el programa realizado fue satisfactorio y semejante al usado por el programa PRO/II.

En la segunda fila encontramos las graficas con los coeficientes de actividad hallados, el primero, con el programa de la 2 entrega, el cual se basaba en la iteración de los parámetros del modelo NRTL, este es el que presenta las menores desviación respecto a los demás, sin embargo, como es básicamente la representación de unos datos ajustados es natural que presente los menores porcentajes de error respecto al experimental; al comparar este último con los datos calculados según el modelo NRTL obtenidos de PRO/II podemos ver, que estos últimos se alejan más de los experimentales al compararlos con el mismo modelo calculado en MATLAB, y aún más al compararlos con el modelo UNIFAC. Esta precisión del modelo UNIFAC se debe a que tiene presente las interacciones de cada subgrupo incluidos en las moléculas de cada compuesto, lo que conlleva a tener que procesar una mayor cantidad de parámetros, y a calcular propiedades del sistema con mayor exactitud. Este modelo de contribución de grupo se basa en estos parámetros intermoleculares, a diferencia del modelo NRTL que evalúa de una forma macro el comportamiento de las moléculas y del cual solo se derivan 3 parámetros que dependen del compuesto en sí, y no de sus grupos funcionales. Sin embargo, notamos que para sistemas con pocos subgrupos, como el presente, es lógico pensar que la precisión es mucho mayor que si existiese una gran variedad de subgrupos en las moléculas.

Evaluar los parámetros del modelo NRTL por métodos iterativos resulta ser muy útil para sistemas en los cuales se conoce el estado a una temperatura o presión dada, los parámetros hallados se pueden interpolar para diferentes presiones y temperaturas pero solo de forma moderada ya que este método basado en el concepto de composición local tiene una flexibilidad limitada para el ajuste de datos.

Analizando el sistema en sí, vemos que las desviaciones respecto a la idealidad son negativas, es decir, los compuestos presentan cierta afinidad y “les gusta” estar más en solución que puros, esto debido a los puentes de hidrogeno formados entre la acetona y el cloroformo, fenómeno que conlleva a tener un máximo de temperatura de ebullición en el punto azeotrópico, y es por esto también que los coeficientes de actividad son menores a 1 para ambos compuestos. Por otra parte, en los coeficientes de fugacidad, hallados con la correlación de O´conell-Prausnitz para el segundo coeficiente virial, vemos como en estado puro el vapor del cloroformo presenta un mayor acercamiento al estado ideal que la acetona. En la mezcla se observa claramente que para sistemas diluidos los coeficientes de fugacidad se acercan más a uno que en su estado puro, por ejemplo, para la acetona cuando está en bajas concentraciones su coeficiente de fugacidad en mezcla se acerca a 1, de igual forma para el cloroformo; se evidencia entonces, de igual manera, su afinidad y su gusto por permanecer juntos, aunque en fase vapor estas interacciones deberían actuar en sentido contrario, es decir, alejarlas de la idealidad, vemos como para esta particular solución se presenta el fenómeno contrario.

En todas las graficas es evidente como mientras la composición de un componente baja el error de los datos calculados respecto a los datos experimentales es mayor, esto aplica para cualquier modelo. Una explicación sería la dificultad técnica en la determinación de las composiciones de las sustancias, tarea que se hace mucho más difícil e imprecisa si estas son bajas, lo que trae consigo problemas de consistencia termodinámica, propiedad que no fue evaluada en este proyecto.

Pasando al análisis de la temperatura de de rocío y burbuja, Figura (2), se ve claramente como las composiciones y las temperaturas calculadas con el modelo NRTL se acercan mucho a las temperaturas experimentales, la diferencia entre las temperaturas máximas en el punto azeotrópico es tan solo de algunas decimas, desafortunadamente PRO/II al utilizar el modelo NRTL no reproduce estos datos

de manera satisfactoria y nos brinda una gráfica con temperaturas de rocío y burbuja mucho más altas que las calculadas, esto puede deberse a que el programa PRO/II utiliza el modelo Soave-Redlich-Kwong para gases reales mientras que en el presente trabajo se utilizó una correlación que es mucho más acertada para dichos compuestos, sin embargo como los compuestos no difieren mucho del comportamiento ideal, se espera que el error en los datos de PRO/II se deban a alguna causa de mayor peso, como puede ser la utilización de otros valores en los parámetros del modelo NRTL, esta razonamiento se ve favorecido por el hecho de que cuando PRO/II utiliza el modelo UNIFAC las curvas para temperatura de rocío y burbuja obtenidas se acercan mucho a las experimentales, pero no tanto como las del modelo NRTL usado en nuestro algoritmo, concluimos entonces que para predecir las temperaturas de burbuja y rocío es mejor el modelo NRTL si

se realiza el algoritmo de forma manual y no utilizando el

programa PRO/II que tiene errores en este tipo de modelo, u otra forma es utilizar un modelo más sofisticado como UNIFAC el cual presenta pequeñas desviaciones si se utiliza PRO/II, sin embargo este presenta mayores desviaciones que el modelo NRTL propuesto en nuestro algoritmo utilizando los parámetros suministrados por el DECHEMA.

La exactitud en la predicción de las composiciones en equilibrio líquido-vapor es muy buena para el algoritmo planteado usando el modelo NRTL, esta exactitud se debe en parte a que el sistema trabajado no presenta mucha desviación de la idealidad (línea 45º). El modelo NRTL que como vimos anteriormente tiene ciertas fallas presenta una desviación significativa respecto a los datos experimentales.

4. CONCLUSIONES

Existen muchas formas de evaluar las propiedades de un sistema liquido-vapor, depende de las características de los compuestos y la condiciones del sistema el modelo a escoger, es así como para condiciones de presiones bajas (1 atm) se

puede suponer comportamiento de gas ideal en la fase vapor, para nuestro sistema aunque esta consideración pudo ser acertada, decidimos utilizar una correlación de la ecuación virial para gases reales y así calcular datos más cercanos a los experimentales. Entre los modelos NRTL y UNIFAC definitivamente este último presente un mayor acercamiento en el cálculo de coeficientes de actividad, sin embargo, tal cual como anotamos en los análisis, si se tienen los suficientes datos para evaluar los parámetros de NRTL una ajustamiento de datos experimentales por este modelo brindaría del mismo modo predicciones muy acertadas.

Para el cálculo de la temperatura de rocío y burbuja, de forma cualitativa (basados en la grafica) se observa una mayor exactitud en el modelo NRTL evaluando las propiedades del sistema según los algoritmos presentados en el SMITH – VAN – NESS. U otra forma seria evaluar las composiciones y estas temperaturas utilizando el modelo UNIFAC del programa PRO/II, porque definitivamente el modelo NRTL de dicho programa presenta serios problemas en la predicción de estas propiedades.

Para la predicción de la composición en la fase vapor, el modelo NRTL es suficientemente exacto, y debido a su sencillez resulta más adecuado utilizar que otros modelos más complejos como el UNIFAC, que por supuesto arrojaría datos más exactos. Cabe anotar que el modelo NRTL de PRO/II también presenta fallas en la predicción de estas composiciones.

5. BIBLIOGRAFÌA

AICHE Journal Chemical Engineering Research and Development, Marzo 1974, “An empirical Correlation of Second Virial Coefficients” Constantine Tsonopoulos pág 263.

I&EC Process Desing and Development, Vol 6 No 2 Abril 1967, “Empirical correlation of second virial coefficients for Vapor-Liquid equilibrium calculations”, J. P. O’Connell, y J.M. Prausnitz, pág 245.

SMITH, VAN NESS Y ABBOT; “Introducción a la Termodinámica en la Ingeniería Química”, 4 edición, editorial McGRAW-HILL

RIVERA, Martín “Optimización con MATLAB“, Universidad Iberoamericana Ciudad de México.

TREYBAL, Robert E. “Mass-Transfer Operations” McGraw-Hill International Editions, 2 edition.

PROGRAMAS PARA USA EN EL SOFTWARE MATLAB

burbuja.m clear all clc

P=760;%760mmHg Patm=1;%atm A1=7.11714; B1=1210.595; C1=229.664; %Constantes de antoine A2=6.95465; B2=1170.966; C2=226.232; %Constantes de antoine TC=[62.7 63.55 63.95 64.4 64.55 64.25 63.55 62.3 61.8 61.27 57.1 56.7];%°C TK=TC+273.15;%K X1=[0.1013 0.1792 0.2585 0.3022 0.3697 0.4418 0.5268 0.6318 0.7605 0.8137 0.8946 0.9433 0.9652]; X1Real=X1; X2=1-X1; X2Real=X2; R=1.987; A12=-481.7574; A21=106.6503; alpha=0.3030; nB=0; nR=0;

Acetona Cloroformo 60.8 60.1 59.3 57.95

0.6683 0.7020 0.7315

%Para la acetona Tca=508; Pca=46.978; wa=0.3124; wha=0.187; ua=2.88; na=0;% K, atm Pra=P./Pca; ura=10^5*ua^2*Pca/Tca^2; Vca=0.082057*(Tca/Pca)*(0.293-0.08*wa);% m3/mol %para el cloroformo Tcc=536.6; Pcc=53.985; wc=0.219; whc=0.187; uc=1.02; nc=0.28;% K, atm Prc=P./Pcc; urc=10^5*uc^2*Pcc/Tcc^2; Vcc=0.082057*(Tcc/Pcc)*(0.293-0.08*wc);% m3/mol %para la mezcla Tcm=(Tca*Tcc)^0.5;%K Pcm=4*Tcm*((Pca*Vca)/Tca+(Pcc*Vcc)/Tcc)/(Vca^(1/3)+Vcc^(1/3))^3;%atm whm=0.5*(wha+whc); urm=10^5*ua*uc*Pcm/Tcm^2; nm=0.5*(na+nc); Prm=P./Pcm;

%---------------------Calculo Temp Burbuja-------------------

FI1=1; FI2=1; Tsat1=B1/(A1-log10(P))-C1; Tsat2=B2/(A2-log10(P))-C2; Tem=X1.*Tsat1+X2.*Tsat2;%ºC TemK=273.15+Tem;%K. r1=exp(X2.^2.*((A21./(R.*TemK)).*((exp(alpha.*(A21./(R.*TemK))))./(X1+X2.*(exp(-alpha.*(A21./(R.*TemK)))))).^2+((exp(alpha.*(A12./(R.*TemK)))).*(A12./(R.*TemK)))./(X2+X1.*(exp(alpha.*(A12./(R.*TemK))))).^2)); r2=exp(X1.^2.*((A12./(R.*TemK)).*((exp(alpha.*(A12./(R.*TemK))))./(X2+X1.*(exp(-alpha.*(A12./(R.*TemK)))))).^2+((exp(alpha.*(A21./(R.*TemK)))).*(A21./(R.*TemK)))./(X1+X2.*(exp(alpha.*(A21./(R.*TemK))))).^2)); Psat1=10.^(A1-B1./(Tem+C1));%mmHg Psat2=10.^(A2-B2./(Tem+C2));%mmHg Psatj=P./((X1.*r1./FI1).*(Psat1/Psat1)+(X2.*r2./FI2).*(Psat2/Psat1)); %mmHg Tem=B1./(A1-log10(Psatj))-C1; %ºC TemK=Tem+273; delta=1; while delta>0.0000001 TemKComp=TemK; Psat1=10.^(A1-B1./(Tem+C1)); %mmHg Psat2=10.^(A2-B2./(Tem+C2)); %mmHg Psat1atm=Psat1./760;%Atmosferas Psat2atm=Psat2./760;%Atmosferas Y1=X1.*r1.*Psat1./(FI1.*P); Y2=X2.*r2.*Psat2./(FI2.*P); Tra=TemK./Tca; Trc=TemK./Tcc; Trm=TemK./Tcm; %Trabajando con la correlacion O'conell-Prausnitz %Para la acetona pura f0a=0.1445-0.330./Tra-0.1385./Tra.^2-0.0121./Tra.^3; f1a=0.073+0.46./Tra-0.5./Tra.^2-0.097./Tra.^3-0.0073./Tra.^8; f2a=-5.237220+5.665807*(log(ura))2.133816*(log(ura))^2+0.2525373*(log(ura))^3+(1./Tra)*(5.7697706.181427*(log(ura))+2.28327*(log(ura))^2-0.2649074*(log(ura))^3); f3a=exp(6.6*(0.7-Tra)); pbrta=f0a+wha*f1a+f2a+na*f3a; %Para el cloroformo puro

f0c=0.1445-0.330./Trc-0.1385./Trc.^2-0.0121./Trc.^3; f1c=0.073+0.46./Trc-0.5./Trc.^2-0.097./Trc.^3-0.0073./Trc.^8; f2c=-5.237220+5.665807*(log(urc))2.133816*(log(urc))^2+0.2525373*(log(urc))^3+(1./Trc)*(5.7697706.181427*(log(urc))+2.28327*(log(urc))^2-0.2649074*(log(urc))^3); f3c=exp(6.6*(0.7-Trc)); pbrtc=f0c+whc*f1c+f2c+nc*f3c; %Para la mezcla f0m=0.1445-0.330./Trm-0.1385./Trm.^2-0.0121./Trm.^3; f1m=0.073+0.46./Trm-0.5./Trm.^2-0.097./Trm.^3-0.0073./Trm.^8; f2m=-5.237220+5.665807*(log(urm))2.133816*(log(urm))^2+0.2525373*(log(urm))^3+(1./Trm)*(5.7697706.181427*(log(urm))+2.28327*(log(urm))^2-0.2649074*(log(urm))^3); f3m=exp(6.6.*(0.7-Trm)); pbrtm=f0m+whm*f1m+f2m+nm*f3m; Ba=0.082057*(Tca/Pca).*pbrta; Bc=0.082057*(Tcc/Pcc).*pbrtc; Bm=0.082057*(Tcm/Pcm).*pbrtm;

%B11 %B22 %B12

d12=2.*Bm-Ba-Bc;

Fiasat=exp(Ba.*Psat1atm./(0.082057.*TemK)); Ficsat=exp(Bc.*Psat2atm./(0.082057.*TemK)); Fiamezcla=exp((Patm./(0.082057.*TemK)).*(Ba+Y2.^2.*d12)); Ficmezcla=exp((Patm./(0.082057.*TemK)).*(Bc+Y1.^2.*d12)); FI1=Fiamezcla./Fiasat; FI2=Ficmezcla./Ficsat; r1=exp(X2.^2.*((A21./(R.*TemK)).*((exp(alpha.*(A21./(R.*TemK))))./(X1+X2.*(exp(-alpha.*(A21./(R.*TemK)))))).^2+((exp(alpha.*(A12./(R.*TemK)))).*(A12./(R.*TemK)))./(X2+X1.*(exp(alpha.*(A12./(R.*TemK))))).^2)); r2=exp(X1.^2.*((A12./(R.*TemK)).*((exp(alpha.*(A12./(R.*TemK))))./(X2+X1.*(exp(-alpha.*(A12./(R.*TemK)))))).^2+((exp(alpha.*(A21./(R.*TemK)))).*(A21./(R.*TemK)))./(X1+X2.*(exp(alpha.*(A21./(R.*TemK))))).^2)); Psatj=P./((X1.*r1./FI1).*(Psat1/Psat1)+(X2.*r2./FI2).*(Psat2/Psat1)); Tem=B1./(A1-log10(Psatj))-C1; TemK=Tem+273.15; delta=abs(TemKComp-TemK); nB=nB+1; end

TemBur=Tem; TemBurK=TemK; Y1Bur=Y1; Y2Bur=Y2; %---------------------Calculo Temp Rocio------------------Y1=[0.0740 0.1428 0.2221 0.2814 0.3724 0.4695 0.5862 0.7070 0.7526 0.7852 0.8123 0.8376 0.8793 0.9411 0.9699 0.9822]; Y1Real=Y1; Y2=1-Y1; Y2Real=Y2; FI1=1; FI2=1; r1=1; r2=1; Tsat1=B1/(A1-log10(P))-C1; Tsat2=B2/(A2-log10(P))-C2; Tem=Y1.*Tsat1+Y2.*Tsat2;%ºC TemK=273.15+Tem;%K. Psat1=10.^(A1-B1./(Tem+C1));%mmHg Psat2=10.^(A2-B2./(Tem+C2));%mmHg %j=1 Psatj=P.*(Y1.*FI1./(r1).*(Psat1/Psat1)+Y2.*FI2./(r2).*(Psat1/Psat2)); Tem=B1./(A1-log10(Psatj))-C1; %ºC TemK=Tem+273; Psat1=10.^(A1-B1./(Tem+C1)); %mmHg Psat2=10.^(A2-B2./(Tem+C2)); %mmHg Psat1atm=Psat1./760;%Atmosferas Psat2atm=Psat2./760;%Atmosferas Tra=TemK./Tca; Trc=TemK./Tcc; Trm=TemK./Tcm; %Trabajando con la correlacion O'conell-Prausnitz %Para la acetona pura f0a=0.1445-0.330./Tra-0.1385./Tra.^2-0.0121./Tra.^3; f1a=0.073+0.46./Tra-0.5./Tra.^2-0.097./Tra.^3-0.0073./Tra.^8; f2a=-5.237220+5.665807*(log(ura))2.133816*(log(ura))^2+0.2525373*(log(ura))^3+(1./Tra)*(5.7697706.181427*(log(ura))+2.28327*(log(ura))^2-0.2649074*(log(ura))^3); f3a=exp(6.6*(0.7-Tra)); pbrta=f0a+wha*f1a+f2a+na*f3a; %Para el cloroformo puro

f0c=0.1445-0.330./Trc-0.1385./Trc.^2-0.0121./Trc.^3; f1c=0.073+0.46./Trc-0.5./Trc.^2-0.097./Trc.^3-0.0073./Trc.^8; f2c=-5.237220+5.665807*(log(urc))2.133816*(log(urc))^2+0.2525373*(log(urc))^3+(1./Trc)*(5.7697706.181427*(log(urc))+2.28327*(log(urc))^2-0.2649074*(log(urc))^3); f3c=exp(6.6*(0.7-Trc)); pbrtc=f0c+whc*f1c+f2c+nc*f3c; %Para la mezcla f0m=0.1445-0.330./Trm-0.1385./Trm.^2-0.0121./Trm.^3; f1m=0.073+0.46./Trm-0.5./Trm.^2-0.097./Trm.^3-0.0073./Trm.^8; f2m=-5.237220+5.665807*(log(urm))2.133816*(log(urm))^2+0.2525373*(log(urm))^3+(1./Trm)*(5.7697706.181427*(log(urm))+2.28327*(log(urm))^2-0.2649074*(log(urm))^3); f3m=exp(6.6.*(0.7-Trm)); pbrtm=f0m+whm*f1m+f2m+nm*f3m; Ba=0.082057*(Tca/Pca).*pbrta; Bc=0.082057*(Tcc/Pcc).*pbrtc; Bm=0.082057*(Tcm/Pcm).*pbrtm;

%B11 %B22 %B12

d12=2.*Bm-Ba-Bc;

Fiasat=exp(Ba.*Psat1atm./(0.082057.*TemK)); Ficsat=exp(Bc.*Psat2atm./(0.082057.*TemK)); Fiamezcla=exp((Patm./(0.082057.*TemK)).*(Ba+Y2.^2.*d12)); Ficmezcla=exp((Patm./(0.082057.*TemK)).*(Bc+Y1.^2.*d12)); FI1=Fiamezcla./Fiasat; FI2=Ficmezcla./Ficsat; X1=Y1.*FI1.*P./(r1.*Psat1); X2=Y2.*FI2.*P./(r2.*Psat2); r1=exp(X2.^2.*((A21./(R.*TemK)).*((exp(alpha.*(A21./(R.*TemK))))./(X1+X2.*(exp(-alpha.*(A21./(R.*TemK)))))).^2+((exp(alpha.*(A12./(R.*TemK)))).*(A12./(R.*TemK)))./(X2+X1.*(exp(alpha.*(A12./(R.*TemK))))).^2)); r2=exp(X1.^2.*((A12./(R.*TemK)).*((exp(alpha.*(A12./(R.*TemK))))./(X2+X1.*(exp(-alpha.*(A12./(R.*TemK)))))).^2+((exp(alpha.*(A21./(R.*TemK)))).*(A21./(R.*TemK)))./(X1+X2.*(exp(alpha.*(A21./(R.*TemK))))).^2)); Psatj=P.*(Y1.*FI1./(r1).*(Psat1/Psat1)+Y2.*FI2./(r2).*(Psat1/Psat2)); Tem=B1./(A1-log10(Psatj))-C1; %ºC TemK=Tem+273; deltaT=1;

while deltaT>0.0000001 TemKComp=TemK; Psat1=10.^(A1-B1./(Tem+C1)); %mmHg Psat2=10.^(A2-B2./(Tem+C2)); %mmHg Psat1atm=Psat1./760;%Atmosferas Psat2atm=Psat2./760;%Atmosferas Tra=TemK./Tca; Trc=TemK./Tcc; Trm=TemK./Tcm; %Trabajando con la correlacion O'conell-Prausnitz %Para la acetona pura f0a=0.1445-0.330./Tra-0.1385./Tra.^2-0.0121./Tra.^3; f1a=0.073+0.46./Tra-0.5./Tra.^2-0.097./Tra.^3-0.0073./Tra.^8; f2a=-5.237220+5.665807*(log(ura))2.133816*(log(ura))^2+0.2525373*(log(ura))^3+(1./Tra)*(5.7697706.181427*(log(ura))+2.28327*(log(ura))^2-0.2649074*(log(ura))^3); f3a=exp(6.6*(0.7-Tra)); pbrta=f0a+wha*f1a+f2a+na*f3a; %Para el cloroformo puro f0c=0.1445-0.330./Trc-0.1385./Trc.^2-0.0121./Trc.^3; f1c=0.073+0.46./Trc-0.5./Trc.^2-0.097./Trc.^3-0.0073./Trc.^8; f2c=-5.237220+5.665807*(log(urc))2.133816*(log(urc))^2+0.2525373*(log(urc))^3+(1./Trc)*(5.7697706.181427*(log(urc))+2.28327*(log(urc))^2-0.2649074*(log(urc))^3); f3c=exp(6.6*(0.7-Trc)); pbrtc=f0c+whc*f1c+f2c+nc*f3c; %Para la mezcla f0m=0.1445-0.330./Trm-0.1385./Trm.^2-0.0121./Trm.^3; f1m=0.073+0.46./Trm-0.5./Trm.^2-0.097./Trm.^3-0.0073./Trm.^8; f2m=-5.237220+5.665807*(log(urm))2.133816*(log(urm))^2+0.2525373*(log(urm))^3+(1./Trm)*(5.7697706.181427*(log(urm))+2.28327*(log(urm))^2-0.2649074*(log(urm))^3); f3m=exp(6.6.*(0.7-Trm)); pbrtm=f0m+whm*f1m+f2m+nm*f3m; Ba=0.082057*(Tca/Pca).*pbrta; Bc=0.082057*(Tcc/Pcc).*pbrtc; Bm=0.082057*(Tcm/Pcm).*pbrtm; d12=2.*Bm-Ba-Bc;

%B11 %B22 %B12

Fiasat=exp(Ba.*Psat1atm./(0.082057.*TemK)); Ficsat=exp(Bc.*Psat2atm./(0.082057.*TemK)); Fiamezcla=exp((Patm./(0.082057.*TemK)).*(Ba+Y2.^2.*d12)); Ficmezcla=exp((Patm./(0.082057.*TemK)).*(Bc+Y1.^2.*d12)); FI1=Fiamezcla./Fiasat; FI2=Ficmezcla./Ficsat; deltar1=1; deltar2=1; while deltar1>0.000001 & deltar2>0.000001 r1Comp=r1; r2Comp=r2; X1=Y1.*FI1.*P./(r1.*Psat1); X2=Y2.*FI2.*P./(r2.*Psat2); X1=X1./(X1+X2); X2=X2./(X1+X2); r1=exp(X2.^2.*((A21./(R.*TemK)).*((exp(alpha.*(A21./(R.*TemK))))./(X1+X2.*(exp(-alpha.*(A21./(R.*TemK)))))).^2+((exp(alpha.*(A12./(R.*TemK)))).*(A12./(R.*TemK)))./(X2+X1.*(exp(alpha.*(A12./(R.*TemK))))).^2)); r2=exp(X1.^2.*((A12./(R.*TemK)).*((exp(alpha.*(A12./(R.*TemK))))./(X2+X1.*(exp(-alpha.*(A12./(R.*TemK)))))).^2+((exp(alpha.*(A21./(R.*TemK)))).*(A21./(R.*TemK)))./(X1+X2.*(exp(alpha.*(A21./(R.*TemK))))).^2)); deltar1=abs(r1-r1Comp); deltar2=abs(r2-r2Comp); end Psatj=P.*(Y1.*FI1./(r1).*(Psat1/Psat1)+Y2.*FI2./(r2).*(Psat1/Psat2)); Tem=B1./(A1-log10(Psatj))-C1; %ºC TemK=Tem+273; deltaT=abs((TemKComp-TemK)); nR=nR+1; end TemRoc=Tem; TemRocK=TemK; X1Roc=X1; X2Roc=X2; X1Error=abs((X1Real-X1Roc)./X1Real.*100); Y1Error=abs((Y1Real-Y1Bur)./Y1Real.*100); datos=[Y1Real' Y1Bur' Y1Error' TemBur' X1Real' X1Roc' X1Error' TemRoc']; disp('Los Datos obtenidos para la temperatura de Burbuja son los siguientes') disp(' ')

disp(' Y1 Y1Burbuja Y1%Error Temperatura Burbuja (ºC) ') disp(datos(:,1:4)) disp('Los Datos obtenidos para la temperatura de Rocio son los siguientes') disp(' ') disp(' X1 X1Rocio X1%Error Temperatura Rocio (ºC)') disp(datos(:,5:8)) figure(1) plot(X1Real,TemBur,Y1Real,TemRoc) title('X1 real Vs T'); xlabel('X1') ylabel('T ºC') xlim([0 1]) grid on legend('Temperatura Burbuja','Temperatura Rocio',0)

figure(2) plot(X1Real,Y1Real,X1Real,Y1Bur) title('X1 real Vs Y1 Calculado (NRTL)'); xlabel('X1') ylabel('Y1') xlim([0 1]) ylim([0 1]) grid on legend('Experimental','Calculado con NRTL',0) LiqVap.m clear all clc P=1;%atm o 760mmHg A1=7.11714; B1=1210.595; C1=229.664; %Constantes de antoine A2=6.95465; B2=1170.966; C2=226.232; %Constantes de antoine TC=[62.7 63.55 63.95 64.4 64.55 64.25 63.55 62.3 61.8 61.27 57.1 56.7];%°C X1=[0.1013 0.1792 0.2585 0.3022 0.3697 0.4418 0.5268 0.6318 0.7605 0.8137 0.8946 0.9433 0.9652]; Y1=[0.0740 0.1428 0.2221 0.2814 0.3724 0.4695 0.5862 0.7070 0.8376 0.8793 0.9411 0.9699 0.9822]; T=TC+273.15;%K

Acetona Cloroformo 60.8 60.1 59.3 57.95 0.6683 0.7020 0.7315 0.7526 0.7852 0.8123

Psat1=10.^(A1-B1./(TC+C1));%mmHg Psat2=10.^(A2-B2./(TC+C2));%mmHg Psat1atm=Psat1./760;%Atmosferas Psat2atm=Psat2./760;%Atmosferas %Para la acetona Tca=508; Pca=46.978; wa=0.3124; wha=0.187; ua=2.88; na=0;% K, atm Tra=T./Tca; Pra=P./Pca; ura=10^5*ua^2*Pca/Tca^2; Vca=0.082057*(Tca/Pca)*(0.293-0.08*wa);% m3/mol

%para el cloroformo Tcc=536.6; Pcc=53.985; wc=0.219; whc=0.187; uc=1.02; nc=0.28;% K, atm Trc=T./Tcc; Prc=P./Pcc; urc=10^5*uc^2*Pcc/Tcc^2; Vcc=0.082057*(Tcc/Pcc)*(0.293-0.08*wc);% m3/mol %para la mezcla Tcm=(Tca*Tcc)^0.5;%K Pcm=4*Tcm*((Pca*Vca)/Tca+(Pcc*Vcc)/Tcc)/(Vca^(1/3)+Vcc^(1/3))^3;%atm whm=0.5*(wha+whc); urm=10^5*ua*uc*Pcm/Tcm^2; nm=0.5*(na+nc); Trm=T./Tcm; Prm=P./Pcm; %Trabajando con la correlacion O'conell-Prausnitz %Para la acetona pura f0a=0.1445-0.330./Tra-0.1385./Tra.^2-0.0121./Tra.^3; f1a=0.073+0.46./Tra-0.5./Tra.^2-0.097./Tra.^3-0.0073./Tra.^8; f2a=-5.237220+5.665807*(log(ura))2.133816*(log(ura))^2+0.2525373*(log(ura))^3+(1./Tra)*(5.7697706.181427*(log(ura))+2.28327*(log(ura))^2-0.2649074*(log(ura))^3); f3a=exp(6.6*(0.7-Tra)); pbrta=f0a+wha*f1a+f2a+na*f3a; %Para el cloroformo puro f0c=0.1445-0.330./Trc-0.1385./Trc.^2-0.0121./Trc.^3; f1c=0.073+0.46./Trc-0.5./Trc.^2-0.097./Trc.^3-0.0073./Trc.^8; f2c=-5.237220+5.665807*(log(urc))2.133816*(log(urc))^2+0.2525373*(log(urc))^3+(1./Trc)*(5.7697706.181427*(log(urc))+2.28327*(log(urc))^2-0.2649074*(log(urc))^3); f3c=exp(6.6*(0.7-Trc)); pbrtc=f0c+whc*f1c+f2c+nc*f3c; %Para la mezcla f0m=0.1445-0.330./Trm-0.1385./Trm.^2-0.0121./Trm.^3; f1m=0.073+0.46./Trm-0.5./Trm.^2-0.097./Trm.^3-0.0073./Trm.^8; f2m=-5.237220+5.665807*(log(urm))2.133816*(log(urm))^2+0.2525373*(log(urm))^3+(1./Trm)*(5.7697706.181427*(log(urm))+2.28327*(log(urm))^2-0.2649074*(log(urm))^3); f3m=exp(6.6.*(0.7-Trm)); pbrtm=f0m+whm*f1m+f2m+nm*f3m; Ba=0.082057*(Tca/Pca).*pbrta; Bc=0.082057*(Tcc/Pcc).*pbrtc; Bm=0.082057*(Tcm/Pcm).*pbrtm;

%B11 %B22 %B12

d12=2.*Bm-Ba-Bc; Fiapuro=exp(Ba.*P./(0.082057.*T)); Ficpuro=exp(Bc.*P./(0.082057.*T)); Fiasat=exp(Ba.*Psat1atm./(0.082057.*T)); Ficsat=exp(Bc.*Psat2atm./(0.082057.*T)); Fiamezcla=exp((P./(0.082057.*T)).*(Ba+(1-Y1).^2.*d12)); Ficmezcla=exp((P./(0.082057.*T)).*(Bc+Y1.^2.*d12)); ra=Y1.*760.*Fiamezcla./(X1.*Psat1.*Fiasat); rc=(1-Y1).*760.*Ficmezcla./((1-X1).*Psat2.*Ficsat); datos1=[X1',Fiapuro',Ficpuro',Fiasat',Ficsat',Fiamezcla',Ficmezcla']; disp('Para el sistema de equilibrio de fases Acetona(1)-Cloroformo(2)') disp('a 760 mmHg se obtuvieron los siguientes datos utilizando la') disp('correlacion de OConnell - Prausnitz para el segundo termino') disp('de la ecuacion Virial:') disp(' ') disp(' X1 Fi1puro Fi2Puro Fi1Sat Fi2Sat Fi1Mezc Fi2Mezc') disp(datos1) datos2=[X1',ra',rc']; disp(' X1 r1 r2') disp(datos2) disp('A continuacion se presenta un grafico donde se muestran los') disp('coeficientes de fugacidad (fi) para los componentes puros y en la') disp('mezcla, en funcion de la composicion de (1) en la fase vapor (Y1)') plot(Y1,Fiamezcla,Y1,Ficmezcla,Y1,Fiapuro,Y1,Ficpuro) ylabel('Fi') xlabel('Y1') title('Coeficientes de fugacidad vs Y1') legend('Fi^-Acetona','Fi^-Cloroformo','FiAcetona Puro','FiCloroformo Puro','Location','Best') grid on disp('Presione enter para continuar') pause

disp('A continuacion se presentan los diferentes graficos requeridos para') disp('la realizacion del trabajo y su analisis') subplot(2,2,1) hold on plot(X1,T,'r') plot(Y1,T) title('T vs Y') xlabel('X y Y') ylabel('T') xlim([0 1]) legend('X','Y',0) grid on

subplot(2,2,3) hold on a=[0 1]; plot(X1,Y1,a,a) title('Y vs X') xlabel('X') ylabel('Y') xlim([0 1]) grid on subplot(2,2,4) plot(X1,ra,X1,rc); title('r1 y r2 vs X1'); xlabel('X1') ylabel('r') xlim([0 1]) grid on legend('r1','r2',0) ntrl.m %De los datos obtenidos anteriormente clear all clc TC=[62.7 63.55 63.95 64.4 64.55 64.25 63.55 62.3 61.8 61.27 60.8 60.1 59.3 57.95 57.1 56.7];%°C T=TC+273.15;%K x1=[0.1013 0.1792 0.2585 0.3022 0.3697 0.4418 0.5268 0.6318 0.6683 0.7020 0.7315 0.7605 0.8137 0.8946 0.9433 0.9652]; r1=[0.6078 0.6434 0.6825 0.7277 0.7808 0.8286 0.8837 0.9214 0.9411 0.9499 0.9567 0.9698 0.9755 0.9911 0.9955 0.9980]; r2=[0.9820 0.9690 0.9619 0.9316 0.8980 0.8674 0.8191 0.7796 0.7441 0.7328 0.7226 0.7182 0.7061 0.6394 0.6261 0.6119]; x2=1-x1; P=1;%atm R(1:16)=1.987;%cal · K-1 · mol-1, Estas dimensiones satisfancen las unidades de B12 y B21 encontradas en el DECHEMA, ademas la iteracion converge mas rapidamente con esta constante.

disp('METODO NORMAL DE ITERACION (ENERGIA DE GIBBS EN EXCESO)') disp(' ') %z(1)=b12, z(2)=b21, z(3)=alfa minimo=inline('(x1.*log(r1)+x2.*log(r2))-(x1.*(x2.^2.*((z(2)./(R.*T)).*((exp(z(3).*(z(2)./(R.*T))))./(x1+x2.*(exp(-z(3).*(z(2)./(R.*T)))))).^2+((exp(z(3).*(z(1)./(R.*T)))).*(z(1)./(R.*T)))./(x2+x1.*(exp(z(3).*(z(1)./(R.*T))))).^2))+x2.*(x1.^2.*((z(1)./(R.*T)).*((exp(z(3).*(z(1)./(R.*T))))./(x2+x1.*(exp(-z(3).*(z(1)./(R.*T)))))).^2+((exp(-

z(3).*(z(2)./(R.*T)))).*(z(2)./(R.*T)))./(x1+x2.*(exp(z(3).*(z(2)./(R.*T))))).^2)))','z','R','T','x1','x2','r1','r2'); options=optimset('LargeScale','on','Display','off','HessUpdate','gillmurray'); [z,fval]=lsqnonlin(minimo,[0,0,0.3],[-10000,10000,0.2],[10000,10000,0.47],options,R,T,x1,x2,r1,r2); rteorica1=exp(x2.^2.*((z(2)./(R.*T)).*((exp(z(3).*(z(2)./(R.*T))))./(x1+x2.*(exp(-z(3).*(z(2)./(R.*T)))))).^2+((exp(z(3).*(z(1)./(R.*T)))).*(z(1)./(R.*T)))./(x2+x1.*(exp(z(3).*(z(1)./(R.*T))))).^2)); rteorica2=exp(x1.^2.*((z(1)./(R.*T)).*((exp(z(3).*(z(1)./(R.*T))))./(x2+x1.*(exp(-z(3).*(z(1)./(R.*T)))))).^2+((exp(z(3).*(z(2)./(R.*T)))).*(z(2)./(R.*T)))./(x1+x2.*(exp(z(3).*(z(2)./(R.*T))))).^2)); error1=abs((r1-rteorica1)./r1)*100; error2=abs((r2-rteorica2)./r2)*100; datos1=[x1',r1',rteorica1',error1',r2',rteorica2',error2']; disp('Despues de la iteracion obtenemos los siguientes coeficientes') disp('B12 = ');disp(z(1)) disp('B21 = ');disp(z(2)); disp('alfa = ');disp(z(3)); disp('El valor minimo hallado de tal ecuacion fue = ');disp(fval) disp('A continuacion se presentan los datos teoricos, hallados con los coeficientes B12, B21 y alfa ') disp(' X1 r1(exp) r1(teo) %Error r2(exp) r2(teo) %Error') disp(datos1(:,:)) subplot(1,2,1) plot(x1,rteorica1,x1,rteorica2,x1,r1,x1,r2); title('rteorica1 y rteorica2 vs X1'); xlabel('X1') ylabel('r') xlim([0 1]) grid on legend('r1teorica','r2teorica','r1exp','r2exp',0) disp('METODO DE BARKER PARA ITERACION (PRESIONES)') disp(' ') A1=7.11714; B1=1210.595; C1=229.664; %Constantes de antoine Acetona A2=6.95465; B2=1170.966; C2=226.232; %Constantes de antoine Cloroformo Psat1=10.^(A1-B1./(TC+C1));%mmHg Psat2=10.^(A2-B2./(TC+C2));%mmHg Psat1atm=Psat1./760;%Atmosferas Psat2atm=Psat2./760;%Atmosferas %zB(1)=b12, zB(2)=b21, zB(3)=alfa minimo=inline('P-(x1.*Psat1atm.*exp(x2.^2.*((zB(2)./(R.*T)).*((exp(zB(3).*(zB(2)./(R.*T))))./(x1+x2.*(exp(-zB(3).*(zB(2)./(R.*T)))))).^2+((exp(zB(3).*(zB(1)./(R.*T)))).*(zB(1)./(R.*T)))./(x2+x1.*(exp(zB(3).*(zB(1)./(R.*T))))).^2))+x2.*Psat2atm.*exp(x1.^2.*((zB(1)./(R.*T)).*((exp( -zB(3).*(zB(1)./(R.*T))))./(x2+x1.*(exp(-zB(3).*(zB(1)./(R.*T)))))).^2+((exp(-

zB(3).*(zB(2)./(R.*T)))).*(zB(2)./(R.*T)))./(x1+x2.*(exp(zB(3).*(zB(2)./(R.*T))))).^2)))','zB','R','T','x1','x2','P','Psat1atm','Psat2atm '); options=optimset('LargeScale','on','Display','off','HessUpdate','steepdesc'); [zB,fval]=lsqnonlin(minimo,z,[-10000,10000,0.2],[10000,10000,0.47],options,R,T,x1,x2,P,Psat1atm,Psat2atm); rteorica1B=exp(x2.^2.*((zB(2)./(R.*T)).*((exp(zB(3).*(zB(2)./(R.*T))))./(x1+x2.*(exp(-zB(3).*(zB(2)./(R.*T)))))).^2+((exp(zB(3).*(zB(1)./(R.*T)))).*(zB(1)./(R.*T)))./(x2+x1.*(exp(zB(3).*(zB(1)./(R.*T))))).^2)); rteorica2B=exp(x1.^2.*((zB(1)./(R.*T)).*((exp(zB(3).*(zB(1)./(R.*T))))./(x2+x1.*(exp(-zB(3).*(zB(1)./(R.*T)))))).^2+((exp(zB(3).*(zB(2)./(R.*T)))).*(zB(2)./(R.*T)))./(x1+x2.*(exp(zB(3).*(zB(2)./(R.*T))))).^2)); error1B=abs((r1-rteorica1B)./r1)*100; error2B=abs((r2-rteorica2B)./r2)*100; datos2=[x1',r1',rteorica1B',error1B',r2',rteorica2B',error2B']; disp('Despues de la iteracion por metodo Barker obtenemos los siguientes coeficientes') disp('B12 = ');disp(zB(1)) disp('B21 = ');disp(zB(2)); disp('alfa = ');disp(zB(3)); disp('El valor minimo hallado de tal ecuacion fue = ');disp(fval) disp('A continuacion se presentan los datos teoricos, hallados con los coeficientes B12, B21 y alfa ') disp(' X1 r1(exp) r1(teo) %Error r2(exp) r2(teo) %Error') disp(datos2(:,:)) subplot(1,2,2) plot(x1,rteorica1B,x1,rteorica2B,x1,r1,x1,r2); title('rteorica1 y rteorica2 vs X1 (Barker)'); xlabel('X1') ylabel('r') xlim([0 1]) grid on legend('rteorica1','rteorica2','r1exp','r2exp',0) unifac.m

% % % % % % % %

Nomenclatura sst: Contador para SUSTANCIA de 1(acetona) a 2(cloroformo); m y k: Contador para SUBGRUPO dat: Contador para DATO de 1 a 16; Casi todas las variables siguen la convencion: VARIABLE(SUSTANCIA,SUBGRUPO,DATO); Convencion subgrupos (18)=(1), (1)=(2), (50)=(3) donde subgrupo (1)=acetona, (2)=metil, (3)=Cloroformo

clc; clear all

TC=[62.7 63.55 63.95 64.4 64.55 64.25 63.55 62.3 61.8 61.27 60.8 60.1 59.3 57.95 57.1 56.7];% °C T=TC+273.15;% K x(1,1,:)=[0.1013 0.1792 0.2585 0.3022 0.3697 0.4418 0.5268 0.6318 0.6683 0.7020 0.7315 0.7605 0.8137 0.8946 0.9433 0.9652]; x(2,1,:)=1-x(1,1,:); x1(1,:)=x(1,1,:);% Variabla de solo dos dimensiones para poder hacer la grafica con plot rexp=[0.6078 0.6434 0.6825 0.7277 0.7808 0.8286 0.8837 0.9214 0.9411 0.9499 0.9567 0.9698 0.9755 0.9911 0.9955 0.9980 0.9820 0.9690 0.9619 0.9316 0.8980 0.8674 0.8191 0.7796 0.7441 0.7328 0.7226 0.7182 0.7061 0.6394 0.6261 0.6119]; Rk=[1.6724,0.9011,2.8700]; Qk=[1.488,0.848,2.410]; Vk=[1 1 0 0 0 1]; a=[ 0 476.40 552.10

26.76 -354.60 0 24.90 36.70 0];

for sst=1:2 r(sst,1)=sum(Vk(sst,:).*Rk); q(sst,1)=sum(Vk(sst,:).*Qk); G(sst,:)=Vk(sst,:).*Qk; end for dat=1:length(x(1,1,:)) JJ(:,1,dat)=r./sum(r.*x(:,1,dat)); LL(:,1,dat)=q./sum(q.*x(:,1,dat)); tao(:,:,dat)=exp(-a./T(dat)); for sst=1:2 for k=1:length(Rk) for m=1:length(Rk) Sksum(m)=G(sst,m).*tao(m,k,dat); end Sk(sst,k,dat)=sum(Sksum); end end for k=1:length(Rk) thk(1,k,dat)=sum(G(:,k).*x(:,1,dat)); nn(1,k,dat)=sum(Sk(:,k,dat).*x(:,1,dat)); end for sst=1:2 lnrC(sst,dat)=1-JJ(sst,1,dat)+log(JJ(sst,1,dat))-5.*q(sst,1).*(1JJ(sst,1,dat)./LL(sst,1,dat)+log(JJ(sst,1,dat)./LL(sst,1,dat))); for k=1:length(Rk) sumatoria(k)=thk(1,k,dat).*(Sk(sst,k,dat)./nn(1,k,dat))G(sst,k)*log(Sk(sst,k,dat)./nn(1,k,dat));

end lnrR(sst,dat)=q(sst,1).*(1-log(LL(sst,1,dat)))-sum(sumatoria); runifac(sst,dat)=exp(lnrC(sst,dat)+lnrR(sst,dat)); end end figure(1) plot(x1,runifac(1,:),x1,runifac(2,:)); title('rUNIFAC(1) y rUNIFAC(2) vs X1'); xlabel('X1') ylabel('r') xlim([0 1]) grid on legend('rUNIFAC(1)','rUNIFAC(2)',0) figure(2) plot(x1,runifac(1,:),x1,runifac(2,:),x1,rexp(1,:),x1,rexp(2,:)); title('rUNIFAC(1) y rUNIFAC(2) vs X1'); xlabel('X1') ylabel('r') xlim([0 1]) grid on legend('rUNIFAC(1)','rUNIFAC(2)','rexp(1)','rexp(2)',0) error1=abs((rexp(1,:)-runifac(1,:))./rexp(1,:))*100; error2=abs((rexp(2,:)-runifac(2,:))./rexp(2,:))*100; datos1=[x1',rexp(1,:)',runifac(1,:)',error1',rexp(2,:)',runifac(2,:)',error2']; disp('A continuacion se presentan los datos teoricos, hallados con la ecuacion UNIFAC') disp(' X1 r1(exp) r1(teo) %Error r2(exp) r2(teo) %Error') disp(datos1) disp('Las graficas 1 y 2, muestras estos resultados') resultados.m

%En este programa fueron copiados todos los datos calculados anteriormente %y graficados junto con los hallados utilizando PROII, todos los algoritmos %no se pudieron juntar en uno solo ya que muchos de ellos utilizaban las %mismas variables en procesos diferentes produciendo datos erroneos, sin %embargo los datos aqui copiados fueron los mismos arrojados en MATLAB por %dichos programas close all clear all clc bdwidth = 5; topbdwidth = 25; set(0,'Units','pixels') scnsize = get(0,'ScreenSize'); pos1 = [1/10*scnsize(3)+bdwidth,1/7*scnsize(4) + bdwidth,8/10*scnsize(3) 2*bdwidth,... 5/7*scnsize(4) - (topbdwidth + bdwidth)];

%---------------Coericientes de actividad experimentales(x1;r1;r2)------------------rexperimentales=[ 0.1013 0.6078 0.9820 0.1792 0.6434 0.9690 0.2585 0.6825 0.9619 0.3022 0.7277 0.9316 0.3697 0.7808 0.8980 0.4418 0.8286 0.8674 0.5268 0.8837 0.8191 0.6318 0.9214 0.7796 0.6683 0.9411 0.7441 0.7020 0.9499 0.7328 0.7315 0.9567 0.7226 0.7605 0.9698 0.7182 0.8137 0.9755 0.7061 0.8946 0.9911 0.6394 0.9433 0.9955 0.6261 0.9652 0.9980 0.6119]'; %---------------Coeficientes de actividad NRTL Calculados (x1;r1;r2)------------------rntrl=[ 0.1013 0.5626 0.9899 0.1792 0.6365 0.9706 0.2585 0.7062 0.9431 0.3022 0.7421 0.9254 0.3697 0.7928 0.8951 0.4418 0.8407 0.8598 0.5268 0.8887 0.8156 0.6318 0.9353 0.7592 0.6683 0.9483 0.7395 0.7020 0.9589 0.7213 0.7315 0.9671 0.7055 0.7605 0.9741 0.6899 0.8137 0.9848 0.6620 0.8946 0.9953 0.6207 0.9433 0.9987 0.5967 0.9652 0.9995 0.5861]'; %---------------Coeficientes de actividad UNIFAC Calculados (x1;r1;r2)------------------runifac=[ 0.1013 0.5488 0.9948 0.1792 0.6056 0.9796 0.2585 0.6684 0.9533 0.3022 0.7039 0.9350 0.3697 0.7564 0.9019 0.4418 0.8084 0.8613 0.5268 0.8629 0.8086 0.6318 0.9182 0.7396 0.6683 0.9341 0.7151 0.7020 0.9472 0.6925

0.7315 0.7605 0.8137 0.8946 0.9433 0.9652

0.9575 0.9664 0.9800 0.9938 0.9982 0.9993

0.6727 0.6529 0.6181 0.5666 0.5367 0.5235]';

%---------------Temperatura de burbuja Calculada (X1;Temperatura de burbuja, Y1 Calculada)-------------------temburbujaNRTL(2,:)=[ 62.5794 63.4637 64.0931 64.2945 64.3810 64.1664 63.5302 62.2692 61.7365 61.2125 60.7334 60.2479 59.3316 57.9127 57.0652 56.6893]'; temburbujaNRTL(3,:)=[ 0.0690 0.1413 0.2305 0.2850 0.3741 0.4723 0.5859 0.7141 0.7542 0.7889 0.8173 0.8436 0.8874 0.9435 0.9718 0.9832]'; temburbujaNRTL(1,:)=[0.1013 0.1792 0.2585 0.3022 0.3697 0.4418 0.5268 0.6318 0.6683 0.7020 0.7315 0.7605 0.8137 0.8946 0.9433 0.9652]; %---------------Temperatura de Rocio Calculada (Y1;Temperatura de Rocio)-------------------

temrocioNRTL(2,:)=[

62.6563 63.4801 64.0540 64.2883 64.3855 64.1805 63.5314 62.3590 61.7617 61.2734 60.8240 60.3650 59.5135 57.9819 57.1255 56.7248]'; temrocioNRTL(1,:)=[0.0740 0.1428 0.2221 0.2814 0.3724 0.4695 0.5862 0.7070 0.7526 0.7852 0.8123 0.8376 0.8793 0.9411 0.9699 0.9822]; %------Coeficientes de actividad Calculados con PROII usando NRTL (X1;r1;r2)---rNRTLPROII=[ 0 0.38436374 1 0.052631579 0.44144422 0.99640322 0.10526316 0.49876854 0.98633307 0.15789473 0.55526263 0.97083545 0.21052632 0.61000097 0.95087707 0.2631579 0.66222703 0.9273234 0.31578946 0.71135277 0.90093589 0.36842105 0.75694317 0.87237853 0.42105263 0.79869729 0.84223086 0.47368422 0.83642501 0.81099927 0.52631581 0.87002784 0.77912778 0.57894737 0.89948303 0.74700528 0.63157892 0.92483103 0.71496987 0.68421054 0.9461658 0.68331212 0.7368421 0.96362692 0.65227658 0.78947371 0.97739214 0.62206411 0.84210527 0.98767054 0.59283394 0.89473683 0.9946956 0.56470674 0.94736844 0.9987182 0.53776848 1 1 0.51207358]'; %------Coeficientes de actividad Calculados con PROII usando UNIFAC (X1;r1;r2)---rUNIFACPROII=[ 0 0.50041395 1 0.052631579 0.51954639 0.99895895 0.10526316 0.55082732 0.99423611 0.15789473 0.58876866 0.98475116 0.21052632 0.63007325 0.97041172

0.2631579 0.31578946 0.36842105 0.42105263 0.47368422 0.52631581 0.57894737 0.63157892 0.68421054 0.7368421 0.78947371 0.84210527 0.89473683 0.94736844 1

0.67256278 0.71472561 0.75549024 0.79409289 0.82999158 0.8628059 0.89227492 0.91822881 0.94056982 0.95926064 0.97431576 0.98579651 0.99380583 0.9984833 1

0.95157266 0.92876858 0.90258479 0.87359965 0.84236395 0.80939597 0.77518094 0.74017191 0.70478833 0.66941249 0.63438582 0.60000503 0.56651956 0.53413117 0.50299323]';

%------Temperaturas de Burbuja y Rocio Calculados con PROII usando NTRL (X1;Temperatura de Burbuja;Y1;Temperatura de Rocio)----temNRTLPROII=[ 0 61.167778 0 61.167778 0.052631579 62.087181 0.02823332 62.087181 0.10526316 62.982292 0.065623745 62.982292 0.15789473 63.788261 0.11238481 63.788261 0.21052632 64.451042 0.16805036 64.451042 0.2631579 64.930443 0.23146211 64.930443 0.31578946 65.201439 0.30086958 65.201439 0.36842105 65.253258 0.3741098 65.253258 0.42105263 65.090111 0.44887289 65.090111 0.47368422 64.727104 0.52292991 64.727104 0.52631581 64.187767 0.59434128 64.187767 0.57894737 63.50103 0.6615867 63.50103 0.63157892 62.698082 0.72361571 62.698082 0.68421054 61.809898 0.77983093 61.809898 0.7368421 60.865471 0.8300252 60.865471 0.78947371 59.890434 0.87426972 59.890434 0.84210527 58.906487 0.91295588 58.906487 0.89473683 57.931389 0.94645095 57.931389 0.94736844 56.978905 0.97529614 56.978905 1 56.058556 1 56.058556]'; %------Temperaturas de Burbuja y Rocio Calculados con PROII usando UNIFAC (X1;Temperatura de Burbuja;Y1;Temperatura de Rocio)----temUNIFACPROII=[ 0 61.167839 0 61.167839 0.052631579 61.825127 0.033163663 61.825127 0.10526316 62.475853 0.071826845 62.475853 0.15789473 63.089043 0.11747831 63.089043 0.21052632 63.624138 0.17055924 63.624138 0.2631579 64.042198 0.23067702 64.042198 0.31578946 64.311913 0.29674646 64.311913 0.36842105 64.412407 0.36715332 64.412407 0.42105263 64.334129 0.43995443 64.334129

0.47368422 0.52631581 0.57894737 0.63157892 0.68421054 0.7368421 0.78947371 0.84210527 0.89473683 0.94736844 1

64.078331 0.51309609 63.65551 0.58462113 63.083458 0.65283602 62.385033 0.71642017 61.586021 0.77447224 60.713005 0.82650083 59.791862 0.87237346 58.846397 0.91224325 57.897667 0.9464674 56.963432 0.97553188 56.058556 1

64.078331 63.65551 63.083458 62.385033 61.586021 60.713005 59.791862 58.846397 57.897667 56.963432 56.058556]';

%----------------Temperatura de Burbuja y rocioexperimental------------TC=[62.7 63.55 63.95 64.4 64.55 64.25 63.55 62.3 61.8 61.27 60.8 60.1 59.3 57.95 57.1 56.7];%°C X1=[0.1013 0.1792 0.2585 0.3022 0.3697 0.4418 0.5268 0.6318 0.6683 0.7020 0.7315 0.7605 0.8137 0.8946 0.9433 0.9652]; Y1=[0.0740 0.1428 0.2221 0.2814 0.3724 0.4695 0.5862 0.7070 0.7526 0.7852 0.8123 0.8376 0.8793 0.9411 0.9699 0.9822]; %---------------Composicion en el equilibrio NRTL PROII XY(X1,Y1)------XY=[ 0 0 0.052631579 0.02823332 0.10526316 0.065623745 0.15789473 0.11238481 0.21052632 0.16805036 0.2631579 0.23146211 0.31578946 0.30086958 0.36842105 0.3741098 0.42105263 0.44887289 0.47368422 0.52292991 0.52631581 0.59434128 0.57894737 0.6615867 0.63157892 0.72361571 0.68421054 0.77983093 0.7368421 0.8300252 0.78947371 0.87426972 0.84210527 0.91295588 0.89473683 0.94645095 0.94736844 0.97529614 1 1]';

figure(1) set(gcf, 'NumberTitle','off','MenuBar','None', 'Name',... 'Coeficientes de actividad con varios modelos','Color',[1,1,1],'position',pos1) subplot(2,2,1)

plot(rexperimentales(1,:),rexperimentales(2,:),rexperimentales(1,:),rexperimenta les(3,:),runifac(1,:),runifac(2,:),runifac(1,:),runifac(3,:)) title('Coeficientes de Actividad Experimentales y UNIFAC') xlabel('X1') ylabel('r') xlim([0 1]) ylim([0.3 1]) grid on legend('r1 experimental', 'r2 experimental', 'r1 UNIFAC', 'r2 UNIFAC',0) subplot(2,2,2) plot(rexperimentales(1,:),rexperimentales(2,:),rexperimentales(1,:),rexperimenta les(3,:),rUNIFACPROII(1,:),rUNIFACPROII(2,:),rUNIFACPROII(1,:),rUNIFACPROII(3,:) ) title('Coeficientes de Actividad Experimentales y UNIFAC con PROII') xlabel('X1') ylabel('r') xlim([0 1]) ylim([0.3 1]) grid on legend('r1 experimental', 'r2 experimental', 'r1 UNIFAC PROII', 'r2 UNIFAC PROII',0) subplot(2,2,3) plot(rexperimentales(1,:),rexperimentales(2,:),rexperimentales(1,:),rexperimenta les(3,:),rntrl(1,:),rntrl(2,:),rntrl(1,:),rntrl(3,:)) title('Coeficientes de Actividad Experimentales y NRTL') xlabel('X1') ylabel('r') xlim([0 1]) ylim([0.3 1]) grid on legend('r1 experimental', 'r2 experimental', 'r1 NRTL', 'r2 NRTL',0) subplot(2,2,4) plot(rexperimentales(1,:),rexperimentales(2,:),rexperimentales(1,:),rexperimenta les(3,:),rNRTLPROII(1,:),rNRTLPROII(2,:),rNRTLPROII(1,:),rNRTLPROII(3,:)) title('Coeficientes de Actividad Experimentales y NRTL con PROII') xlabel('X1') ylabel('r') xlim([0 1]) ylim([0.3 1]) grid on legend('r1 experimental', 'r2 experimental', 'r1 NRTL PROII', 'r2 NRTL PROII',0) figure(2) set(gcf, 'NumberTitle','off','MenuBar','None', 'Name',... 'Puntos de Burbuja y Rocio','Color',[1,1,1],'position',pos1) subplot(2,2,1) plot(temburbujaNRTL(1,:),temburbujaNRTL(2,:),temrocioNRTL(1,:),temrocioNRTL(2,:) ,X1,TC,Y1,TC) title('Temperatura de Burbuja y rocio (NRTL y Experimental)') xlabel('X1') ylabel('T') xlim([0 1]) ylim([56 66]) grid on

legend('T Burbuja Calculados', 'T Rocio Calculados', 'T Bubuja Experimental', 'T Rocio Experimental',0) subplot(2,2,2) plot(temNRTLPROII(1,:),temNRTLPROII(2,:),temNRTLPROII(3,:),temNRTLPROII(4,:),X1, TC,Y1,TC) title('Temperatura de Burbuja y rocio (NRTL PROII y Experimental)') xlabel('X1') ylabel('T') xlim([0 1]) ylim([56 66]) grid on legend('T Burbuja PROII', 'T Rocio PROII', 'T Bubuja Experimental', 'T Rocio Experimental',0) subplot(2,2,3) plot(temUNIFACPROII(1,:),temUNIFACPROII(2,:),temUNIFACPROII(3,:),temUNIFACPROII( 4,:),X1,TC,Y1,TC) title('Temperatura de Burbuja y rocio (UNIFAC PROII y Experimental)') xlabel('X1') ylabel('T') xlim([0 1]) ylim([56 66]) grid on legend('T Burbuja PROII', 'T Rocio PROII', 'T Bubuja Experimental', 'T Rocio Experimental',0) a=[0 1]; subplot(2,2,4) plot(X1,Y1,X1,temburbujaNRTL(3,:),XY(1,:),XY(2,:),a,a) title('Grafica Equilibrio L-V (Y1 vs X1)') xlabel('X1') ylabel('T') xlim([0 1]) ylim([0 1]) grid on legend('Experimental', 'NRTL', 'NRTL PROII',0)