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 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 temperatura fue el siguiente: siguiente: T 1 = 335.7 K Tc(acetona) = 508K Pc(acetona) = 4.76 Mpa Pr = P / Pc BPc RT c
P = 101.325 Mpa Tc(cloroformo) = 536.6 K
Pc(cloroformo ) = 5.47 Mpa Tr = T / Tc
= f 0 + wf 1 + f 2 (1)
Donde f 0 = 0.1445 −
f 1 = 0.0637 +
f 2 =
a Tr 6
0.330
−
Tr
0.331 Tr 2
−
0.1385 Tr 2
0.423 Tr 3
−
−
0.0121 Tr 3
−
0.000607 Tr 8
0.008 Tr 8
a (acetona) = - 0.0309
a(cloroformo) = - 0.04329
Para la acetona hallamos B de la ecuación (1) y calculamos el coeficiente de fugacidad puros así: ϕ = 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 puede considerar de comportamiento comportamiento ideal, 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 usando las correlaciones correlaciones generalizadas fue el siguiente: BPc RTc
= B 0 + WB1 (2)
Donde
B = 0.083 −
W acetona =
0
0.3124
0.422
1
B = 0.139 −
Tr 1.6 W cloroformo =
0.172 Tr 4.2
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 ,T R ) = −5.237220 + 5.665807( Lnµ R ) − 2.133816( Lnµ R ) 2 + 0.2525373( Lnµ R )3
+ 1 / T r [5.769770 − 6.181427( Lnµ R ) + 2.283270 ( Lnµ R ) 2 − 0.2649074( Lnµ R )3 3
f a (T R ) = exp[ 6.6(0.7 − T R )] T R = T / T c
T c / acetona) = 508K
P R = P / Pc
Pc (acetona) = 46.798atm
De tablas del artículo: W hacetona = 0.187
T c / cloroformo) = 536.6 K
Pc (cloroformo ) = 53.985atm W hcloroformo = 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 (B 11 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 T
2
De la ecuación (1), para hallar el coeficiente de fugacidad de la mezcla cambiamos el subíndice i por ij T R = 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 ( B11 + y 2 δ 12 ) RT ^ P 2 ln(ϕ cloroformo ) = ( B22 + y1 δ 12 ) RT ^
ln(ϕ acetona ) =
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 para el segundo coeficiente coeficiente virial (B) pero con una presión 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: γ i =
^ L i id i
f f
^
=
f i v xi f i
^
=
y i P ϕ i xi Pi sat ϕ isat
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 µ = µ , entonces f i = f i y por lo tanto v i
L i
DESCRIPCIÓN DEL PROGRAMA EN MATLAB 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. 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. pbrtc =
BPc RT c
pbrta = pbrtm =
BPc RT c BPc RT c
para el cloroformo para la acetona para la mezcla
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
ANALISIS DE RESULTADOS En la primera gráfica vemos como el coeficiente de fugacidad disminuye 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 presentándose un máximo de temperatura temperatura en el azeótropo. 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 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 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 (X21), se puede asumir como gas ideal en una mezcla ideal y para la acetona con fracción liquida tendiente a 1 (X11) 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 b 12 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) (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 RT Teorico Experiment al Donde
(1)
G E = x1 * ln(γ 1* ) + x 2 * ln(γ 2* ) RT Experiment al *
γ 1
y
*
γ 2
(2)
son los coeficientes coeficientes de actividad experimentales experimentales
G E = x1 * ln(γ 1 ) + x 2 * ln(γ 2 ) RT teorico
(3)
Donde 2 G G12 * τ 12 2 21 x = + ln(γ 1 ) 2 * τ 21 2 + x x G * ( x 2 + x1 * G12 ) 1 2 21 2 G G21 * τ 21 12 + ln(γ 2 ) = x12 * τ 12 2 x x G * + ( x1 + x 2 * G21 ) 2 1 12
=exp(-α* τ 12 G12 =exp(-α* τ 12
=
b12 RT
)
(4)
(5)
=exp(-α* τ 21 ) G 21 =exp(-α* τ 21
=
(6)
b21
(7)
RT
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 G 21 G12 * τ 12 2 (8) γ 1 = exp x 2 * τ 21 x + x * G + ( x + x * G )2 1 2 21 2 1 12 2 G G 21 * τ 21 2 12 (9) γ 2 = exp x1 * τ 12 x + x * G + ( x + x * G )2 1 12 1 2 21 2 A continuación graficamos los coeficientes de actividad teóricos respecto a x1, y en la misma grafica los
* γ 1
y
* γ 2
experimentales.
γ 1
y
γ 2
con
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)
sat sat Pexperimental = x1 * P1 * γ 1 + x 2 * P2 * γ 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
experimental. Graficamos
teóricos, y los porcentajes de error con respecto al γ 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 experimentales son menos predecibles conforme 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 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 γ i R
ln γ iC = 1 − J i + ln J i − 5q1 −
J i Li
+ ln
J i
Li
S S ln γ i R = qi (1 − ln Li ) − ∑ θ k ki − Gki ln ki η k η k k J i =
r i
∑ r x j
j
j
Li =
qi
r i = ∑ vk i Rk
qi = ∑ vk i Qk
( )
∑ q j x j
( )
k
( )
Gki = vk i Qk
θ k = ∑ Gki xi
k
i
j
ski = ∑ Gmiτ mk
η k =
m
∑s
ki
xi τ mk = exp
− amk
i
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 R
Q
a
parámetros k y k y de los parámetros de interacción de grupos mk los leímos del libro “The Properties of Gases and and Liquids” de Prausnitz, John M. Los subgrupos derivados de las sustancias fueron los siguientes: Acetona CH3COCH3 ===>
CH3CO (1) y CH3 (2)
Cloroformo CHCl 3
CHCl3 (3)
===>
#
Formula
k
Rk
Qk
vk
vk ( 2)
1 2 3
CH3CO CH3 CHCl3
18 1 50
1.6724 0.9011 2.8700
1.488 0.848 2.410
1 1 0
0 0 1
(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 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) T k Sat =
Bk Ak − log 10( P)
− C k
(2) T = ∑ x x T k Sat k
(3) Pk Sat = 10^ Ak −
Antoine T (º C ) + C k Bk
(4) Los coeficientes de actividad se calcularon con la ecuación NRTL asi:
γ 1
γ 2
2 G G12 * τ 12 2 21 + = exp x 2 * τ 21 2 x x G * + x x G * ( ) + 1 2 21 2 1 12 2 G G 21 * τ 21 2 12 + = exp x1 * τ 12 2 * x x G + ( x1 + x 2 * G21 ) 2 1 12
G12 =exp(-α* τ 12 ) G21 =exp(-α* τ 21 ) τ 12
(5) P jSat =
(6) T =
=
A12
τ 21
RT
=
A21 RT
P Sat x k γ k Pk ∑k Φ P Sat k j
B j A j − log 10( P jSat )
− C j
Sat
(7) y k =
x k γ k Pk
Φ 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.
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 T k Sat k Sat j
(10)
P
(11)
x k =
(12)
x k =
Sat P j y k Φ k Sat = P ∑ γ k k Pk y k Φ k P Sat
γ k Pk
x k
∑ x
k
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 coeficiente s 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
con la fórmula correspondiente que parte del equilibrio f i
^ =
L
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 CH 3COCH3 ===>
CH3CO (1) y CH 3 (2)
Cloroformo CHCl 3
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 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 ura=10^5*ua^2*Pca/Tca^2; /Tca^2; Vca=0.082057*(Tca/Pca)*(0.293-0.0 Vca=0.082057*(Tca /Pca)*(0.293-0.08*wa);% 8*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 urc=10^5*uc^2*Pcc/Tcc^2; /Tcc^2; Vcc=0.082057*(Tcc/Pcc)*(0.293-0.0 Vcc=0.082057*(Tcc /Pcc)*(0.293-0.08*wc);% 8*wc);% m3/mol %para la mezcla Tcm=(Tca*Tcc)^0.5;%K Pcm=4*Tcm*((Pca*Vca)/Tca+(Pcc*Vcc)/ Pcm=4*Tcm*((Pca*V ca)/Tca+(Pcc*Vcc)/Tcc)/(Vca^(1/3)+V Tcc)/(Vca^(1/3)+Vcc^(1/3))^3;%atm cc^(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 %--------------------Calculo Temp Burbuja---------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.* alpha.*(A21./(R.* TemK))))./(X1+X2.*(exp(-alpha.*(A (exp(-alpha.*(A21./(R.*TemK))))) 21./(R.*TemK)))))).^2+((exp().^2+((exp(alpha.*(A12./(R.*TemK)))).*(A12./(R alpha.*(A12./(R.* TemK)))).*(A12./(R.*TemK)))./(X2+X1 .*TemK)))./(X2+X1.*(exp(.*(exp(alpha.*(A12./(R.*TemK))))).^2)); r2=exp(X1.^2.*((A12./(R.*TemK)).*((exp(alpha.*(A12./(R.*TemK))))./(X2+X1.* alpha.*(A12./(R.* TemK))))./(X2+X1.*(exp(-alpha.*(A (exp(-alpha.*(A12./(R.*TemK))))) 12./(R.*TemK)))))).^2+((exp().^2+((exp(alpha.*(A21./(R.*TemK)))).*(A21./(R alpha.*(A21./(R.* TemK)))).*(A21./(R.*TemK)))./(X1+X2 .*TemK)))./(X1+X2.*(exp(.*(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) Psatj=P./((X1.*r1./FI1).*(Psat1/Psa t1)+(X2.*r2./FI2).*(Psat2/Psat1)); .*(Psat2/Psat1)); %mmHg Tem=B1./(A1-log10(Psatj))-C1; 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 Psat1=10.^(A1-B1./(Tem+C1)); Psat2=10.^(A2-B2./(Tem+C2)); 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 O'conell-Prausnitz %Para la acetona pura f0a=0.1445-0.330./Tra-0.1385./Tra.^2-0.0121./Tra.^ f0a=0.1445-0.330./Tra-0.1385./Tra. ^2-0.0121./Tra.^3; 3; f1a=0.073+0.46./Tra-0.5./Tra.^2-0 f1a=0.073+0.46./T ra-0.5./Tra.^2-0.097./Tra.^3-0.00 .097./Tra.^3-0.0073./Tra.^8; 73./Tra.^8; f2a=-5.237220+5.665807*(log(ura))2.133816*(log(ura))^2+0.2525373*(lo 2.133816*(log(ura) )^2+0.2525373*(log(ura))^3+(1./Tra g(ura))^3+(1./Tra)*(5.769770)*(5.7697706.181427*(log(ura))+2.28327*(log(ur 6.181427*(log(ura) )+2.28327*(log(ura))^2-0.2649074 a))^2-0.2649074*(log(ura))^3); *(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.^ f0c=0.1445-0.330./Trc-0.1385./Trc. ^2-0.0121./Trc.^3; 3; f1c=0.073+0.46./Trc-0.5./Trc.^2-0 f1c=0.073+0.46./T rc-0.5./Trc.^2-0.097./Trc.^3-0.00 .097./Trc.^3-0.0073./Trc.^8; 73./Trc.^8; f2c=-5.237220+5.665807*(log(urc))2.133816*(log(urc))^2+0.2525373*(lo 2.133816*(log(urc ))^2+0.2525373*(log(urc))^3+(1./Trc g(urc))^3+(1./Trc)*(5.769770)*(5.7697706.181427*(log(urc))+2.28327*(log(ur 6.181427*(log(urc ))+2.28327*(log(urc))^2-0.2649074 c))^2-0.2649074*(log(urc))^3); *(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.^ f0m=0.1445-0.330./Trm-0.1385./Trm. ^2-0.0121./Trm.^3; 3; f1m=0.073+0.46./Trm-0.5./Trm.^2-0 f1m=0.073+0.46./T rm-0.5./Trm.^2-0.097./Trm.^3-0.00 .097./Trm.^3-0.0073./Trm.^8; 73./Trm.^8; f2m=-5.237220+5.665807*(log(urm))2.133816*(log(urm))^2+0.2525373*(lo 2.133816*(log(urm ))^2+0.2525373*(log(urm))^3+(1./Trm g(urm))^3+(1./Trm)*(5.769770)*(5.7697706.181427*(log(urm))+2.28327*(log(ur 6.181427*(log(urm ))+2.28327*(log(urm))^2-0.2649074 m))^2-0.2649074*(log(urm))^3); *(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.*d Fiamezcla=exp((Patm./(0.082057.*Te mK)).*(Ba+Y2.^2.*d12)); 12)); Ficmezcla=exp((Patm./(0.082057.*Te Ficmezcla=exp((Pa tm./(0.082057.*TemK)).*(Bc+Y1.^2.*d mK)).*(Bc+Y1.^2.*d12)); 12)); FI1=Fiamezcla./Fiasat; FI2=Ficmezcla./Ficsat; r1=exp(X2.^2.*((A21./(R.*TemK)).*((exp(alpha.*(A21./(R.*TemK))))./(X1+X2. alpha.*(A21./(R.* TemK))))./(X1+X2.*(exp(-alpha.*(A *(exp(-alpha.*(A21./(R.*TemK))))) 21./(R.*TemK)))))).^2+((exp().^2+((exp(alpha.*(A12./(R.*TemK)))).*(A12./( alpha.*(A12./(R.* TemK)))).*(A12./(R.*TemK)))./(X2+X1 R.*TemK)))./(X2+X1.*(exp(.*(exp(alpha.*(A12./(R.*TemK))))).^2)); r2=exp(X1.^2.*((A12./(R.*TemK)).*((exp(alpha.*(A12./(R.*TemK))))./(X2+X1. alpha.*(A12./(R.* TemK))))./(X2+X1.*(exp(-alpha.*(A *(exp(-alpha.*(A12./(R.*TemK))))) 12./(R.*TemK)))))).^2+((exp().^2+((exp(alpha.*(A21./(R.*TemK)))).*(A21./( alpha.*(A21./(R.* TemK)))).*(A21./(R.*TemK)))./(X1+X2 R.*TemK)))./(X1+X2.*(exp(.*(exp(alpha.*(A21./(R.*TemK))))).^2)); Psatj=P./((X1.*r1./FI1).*(Psat1/Ps Psatj=P./((X1.*r1 ./FI1).*(Psat1/Psat1)+(X2.*r2./FI2) at1)+(X2.*r2./FI2).*(Psat2/Psat1)); .*(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 %--------------------Calculo Temp Rocio-----------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) Psatj=P.*(Y1.*FI1./(r1).*(Psat1/Psa t1)+Y2.*FI2./(r2).*(Psat1/Psat2)); .*(Psat1/Psat2)); Tem=B1./(A1-log10(Psatj))-C1; Tem=B1./(A1-log10 (Psatj))-C1; %ºC TemK=Tem+273; Psat1=10.^(A1-B1./(Tem+C1)); %mmHg Psat1=10.^(A1-B1./(Tem+C1)); Psat2=10.^(A2-B2./(Tem+C2)); 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 O'conell-Prausnitz %Para la acetona pura f0a=0.1445-0.330./Tra-0.1385./Tra.^2-0.0121./Tra.^ f0a=0.1445-0.330./Tra-0.1385./Tra. ^2-0.0121./Tra.^3; 3; f1a=0.073+0.46./Tra-0.5./Tra.^2-0 f1a=0.073+0.46./T ra-0.5./Tra.^2-0.097./Tra.^3-0.00 .097./Tra.^3-0.0073./Tra.^8; 73./Tra.^8; f2a=-5.237220+5.665807*(log(ura))2.133816*(log(ura))^2+0.2525373*(lo 2.133816*(log(ura) )^2+0.2525373*(log(ura))^3+(1./Tra g(ura))^3+(1./Tra)*(5.769770)*(5.7697706.181427*(log(ura))+2.28327*(log(ur 6.181427*(log(ura) )+2.28327*(log(ura))^2-0.2649074 a))^2-0.2649074*(log(ura))^3); *(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.^ f0c=0.1445-0.330./Trc-0.1385./Trc. ^2-0.0121./Trc.^3; 3; f1c=0.073+0.46./Trc-0.5./Trc.^2-0 f1c=0.073+0.46./T rc-0.5./Trc.^2-0.097./Trc.^3-0.00 .097./Trc.^3-0.0073./Trc.^8; 73./Trc.^8; f2c=-5.237220+5.665807*(log(urc))2.133816*(log(urc))^2+0.2525373*(lo 2.133816*(log(urc ))^2+0.2525373*(log(urc))^3+(1./Trc g(urc))^3+(1./Trc)*(5.769770)*(5.7697706.181427*(log(urc))+2.28327*(log(ur 6.181427*(log(urc ))+2.28327*(log(urc))^2-0.2649074 c))^2-0.2649074*(log(urc))^3); *(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.^ f0m=0.1445-0.330./Trm-0.1385./Trm. ^2-0.0121./Trm.^3; 3; f1m=0.073+0.46./Trm-0.5./Trm.^2-0 f1m=0.073+0.46./T rm-0.5./Trm.^2-0.097./Trm.^3-0.00 .097./Trm.^3-0.0073./Trm.^8; 73./Trm.^8; f2m=-5.237220+5.665807*(log(urm))2.133816*(log(urm))^2+0.2525373*(lo 2.133816*(log(urm ))^2+0.2525373*(log(urm))^3+(1./Trm g(urm))^3+(1./Trm)*(5.769770)*(5.7697706.181427*(log(urm))+2.28327*(log(ur 6.181427*(log(urm ))+2.28327*(log(urm))^2-0.2649074 m))^2-0.2649074*(log(urm))^3); *(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.*d Fiamezcla=exp((Patm./(0.082057.*Te mK)).*(Ba+Y2.^2.*d12)); 12)); Ficmezcla=exp((Patm./(0.082057.*Te Ficmezcla=exp((Pa tm./(0.082057.*TemK)).*(Bc+Y1.^2.*d mK)).*(Bc+Y1.^2.*d12)); 12)); 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. alpha.*(A21./(R.* TemK))))./(X1+X2.*(exp(-alpha.*(A *(exp(-alpha.*(A21./(R.*TemK))))) 21./(R.*TemK)))))).^2+((exp().^2+((exp(alpha.*(A12./(R.*TemK)))).*(A12./( alpha.*(A12./(R.* TemK)))).*(A12./(R.*TemK)))./(X2+X1 R.*TemK)))./(X2+X1.*(exp(.*(exp(alpha.*(A12./(R.*TemK))))).^2)); r2=exp(X1.^2.*((A12./(R.*TemK)).*((exp(alpha.*(A12./(R.*TemK))))./(X2+X1. alpha.*(A12./(R.* TemK))))./(X2+X1.*(exp(-alpha.*(A *(exp(-alpha.*(A12./(R.*TemK))))) 12./(R.*TemK)))))).^2+((exp().^2+((exp(alpha.*(A21./(R.*TemK)))).*(A21./( alpha.*(A21./(R.* TemK)))).*(A21./(R.*TemK)))./(X1+X2 R.*TemK)))./(X1+X2.*(exp(.*(exp(alpha.*(A21./(R.*TemK))))).^2)); Psatj=P.*(Y1.*FI1./(r1).*(Psat1/Psat1)+Y2.*FI2./(r2) Psatj=P.*(Y1.*FI1./(r1).*(Psat1/Ps at1)+Y2.*FI2./(r2).*(Psat1/Psat2)); .*(Psat1/Psat2)); Tem=B1./(A1-log10(Psatj))-C1; 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 Psat1=10.^(A1-B1./(Tem+C1)); Psat2=10.^(A2-B2./(Tem+C2)); 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 O'conell-Prausnitz %Para la acetona pura f0a=0.1445-0.330./Tra-0.1385./Tra.^2-0.0121./Tra.^ f0a=0.1445-0.330./Tra-0.1385./Tra. ^2-0.0121./Tra.^3; 3; f1a=0.073+0.46./Tra-0.5./Tra.^2-0 f1a=0.073+0.46./T ra-0.5./Tra.^2-0.097./Tra.^3-0.00 .097./Tra.^3-0.0073./Tra.^8; 73./Tra.^8; f2a=-5.237220+5.665807*(log(ura))2.133816*(log(ura))^2+0.2525373*(lo 2.133816*(log(ura ))^2+0.2525373*(log(ura))^3+(1./Tra g(ura))^3+(1./Tra)*(5.769770)*(5.7697706.181427*(log(ura))+2.28327*(log(ur 6.181427*(log(ura ))+2.28327*(log(ura))^2-0.2649074 a))^2-0.2649074*(log(ura))^3); *(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.^ f0c=0.1445-0.330./Trc-0.1385./Trc. ^2-0.0121./Trc.^3; 3; f1c=0.073+0.46./Trc-0.5./Trc.^2-0 f1c=0.073+0.46./T rc-0.5./Trc.^2-0.097./Trc.^3-0.00 .097./Trc.^3-0.0073./Trc.^8; 73./Trc.^8; f2c=-5.237220+5.665807*(log(urc))2.133816*(log(urc))^2+0.2525373*(lo 2.133816*(log(urc ))^2+0.2525373*(log(urc))^3+(1./Trc g(urc))^3+(1./Trc)*(5.769770)*(5.7697706.181427*(log(urc))+2.28327*(log(ur 6.181427*(log(urc ))+2.28327*(log(urc))^2-0.2649074 c))^2-0.2649074*(log(urc))^3); *(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.^ f0m=0.1445-0.330./Trm-0.1385./Trm. ^2-0.0121./Trm.^3; 3; f1m=0.073+0.46./Trm-0.5./Trm.^2-0 f1m=0.073+0.46./T rm-0.5./Trm.^2-0.097./Trm.^3-0.00 .097./Trm.^3-0.0073./Trm.^8; 73./Trm.^8; f2m=-5.237220+5.665807*(log(urm))2.133816*(log(urm))^2+0.2525373*(lo 2.133816*(log(urm ))^2+0.2525373*(log(urm))^3+(1./Trm g(urm))^3+(1./Trm)*(5.769770)*(5.7697706.181427*(log(urm))+2.28327*(log(ur 6.181427*(log(urm ))+2.28327*(log(urm))^2-0.2649074 m))^2-0.2649074*(log(urm))^3); *(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.*d Fiamezcla=exp((Patm./(0.082057.*Tem K)).*(Ba+Y2.^2.*d12)); 12)); Ficmezcla=exp((Patm./(0.082057.*Tem Ficmezcla=exp((Pa tm./(0.082057.*TemK)).*(Bc+Y1.^2.*d K)).*(Bc+Y1.^2.*d12)); 12)); FI1=Fiamezcla./Fiasat; FI2=Ficmezcla./Ficsat; deltar1=1; deltar2=1; while deltar1>0.000001 & deltar2>0.00000 deltar2>0.000001 1 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.* alpha.*(A21./(R.* TemK))))./(X1+X2.*(exp(-alpha.*(A (exp(-alpha.*(A21./(R.*TemK))))) 21./(R.*TemK)))))).^2+((exp().^2+((exp(alpha.*(A12./(R.*TemK)))).*(A12./(R alpha.*(A12./(R.* TemK)))).*(A12./(R.*TemK)))./(X2+X1 .*TemK)))./(X2+X1.*(exp(.*(exp(alpha.*(A12./(R.*TemK))))).^2)); r2=exp(X1.^2.*((A12./(R.*TemK)).*((exp(alpha.*(A12./(R.*TemK))))./(X2+X1.* alpha.*(A12./(R.* TemK))))./(X2+X1.*(exp(-alpha.*(A (exp(-alpha.*(A12./(R.*TemK))))) 12./(R.*TemK)))))).^2+((exp().^2+((exp(alpha.*(A21./(R.*TemK)))).*(A21./(R alpha.*(A21./(R.* TemK)))).*(A21./(R.*TemK)))./(X1+X2 .*TemK)))./(X1+X2.*(exp(.*(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) Psatj=P.*(Y1.*FI1./(r1).*(Psat1/Psa t1)+Y2.*FI2./(r2).*(Psat1/Psat2)); .*(Psat1/Psat2)); Tem=B1./(A1-log10(Psatj))-C1; 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 legend('Temperatu ra Burbuja','Tempera Burbuja','Temperatura tura 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 legend('Experimen tal','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*Pc ura=10^5*ua^2*Pca/Tca^2; a/Tca^2; Vca=0.082057*(Tca/Pca)*(0.293-0.0 Vca=0.082057*(Tca /Pca)*(0.293-0.08*wa);% 8*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*Pc urc=10^5*uc^2*Pcc/Tcc^2; c/Tcc^2; Vcc=0.082057*(Tcc/Pcc)*(0.293-0.0 Vcc=0.082057*(Tcc /Pcc)*(0.293-0.08*wc);% 8*wc);% m3/mol %para la mezcla Tcm=(Tca*Tcc)^0.5;%K Pcm=4*Tcm*((Pca*Vca)/Tca+(Pcc*Vcc)/ Pcm=4*Tcm*((Pca*V ca)/Tca+(Pcc*Vcc)/Tcc)/(Vca^(1/3)+V Tcc)/(Vca^(1/3)+Vcc^(1/3))^3;%atm cc^(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 O'conell-Prausnitz %Para la acetona pura f0a=0.1445-0.330./Tra-0.1385./Tra.^2-0.0121./Tra.^ f0a=0.1445-0.330./Tra-0.1385./Tra. ^2-0.0121./Tra.^3; 3; f1a=0.073+0.46./Tra-0.5./Tra.^2-0 f1a=0.073+0.46./T ra-0.5./Tra.^2-0.097./Tra.^3-0.00 .097./Tra.^3-0.0073./Tra.^8; 73./Tra.^8; f2a=-5.237220+5.665807*(log(ura))2.133816*(log(ura))^2+0.2525373*(lo 2.133816*(log(ura ))^2+0.2525373*(log(ura))^3+(1./Tra g(ura))^3+(1./Tra)*(5.769770)*(5.7697706.181427*(log(ura))+2.28327*(log(ur 6.181427*(log(ura ))+2.28327*(log(ura))^2-0.2649074 a))^2-0.2649074*(log(ura))^3); *(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.^ f0c=0.1445-0.330./Trc-0.1385./Trc. ^2-0.0121./Trc.^3; 3; f1c=0.073+0.46./Trc-0.5./Trc.^2-0 f1c=0.073+0.46./T rc-0.5./Trc.^2-0.097./Trc.^3-0.00 .097./Trc.^3-0.0073./Trc.^8; 73./Trc.^8; f2c=-5.237220+5.665807*(log(urc))2.133816*(log(urc))^2+0.2525373*(lo 2.133816*(log(urc ))^2+0.2525373*(log(urc))^3+(1./Trc g(urc))^3+(1./Trc)*(5.769770)*(5.7697706.181427*(log(urc))+2.28327*(log(ur 6.181427*(log(urc ))+2.28327*(log(urc))^2-0.2649074 c))^2-0.2649074*(log(urc))^3); *(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.^ f0m=0.1445-0.330./Trm-0.1385./Trm. ^2-0.0121./Trm.^3; 3; f1m=0.073+0.46./Trm-0.5./Trm.^2-0 f1m=0.073+0.46./T rm-0.5./Trm.^2-0.097./Trm.^3-0.00 .097./Trm.^3-0.0073./Trm.^8; 73./Trm.^8; f2m=-5.237220+5.665807*(log(urm))2.133816*(log(urm))^2+0.2525373*(lo 2.133816*(log(urm ))^2+0.2525373*(log(urm))^3+(1./Trm g(urm))^3+(1./Trm)*(5.769770)*(5.7697706.181427*(log(urm))+2.28327*(log(ur 6.181427*(log(urm ))+2.28327*(log(urm))^2-0.2649074 m))^2-0.2649074*(log(urm))^3); *(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.*d Fiamezcla=exp((P./(0.082057.*T)).* (Ba+(1-Y1).^2.*d12)); 12)); Ficmezcla=exp((P./(0.082057.*T)).* Ficmezcla=exp((P. /(0.082057.*T)).*(Bc+Y1.^2.*d12)); (Bc+Y1.^2.*d12)); ra=Y1.*760.*Fiamezcla./(X1.*Psat1.*Fiasat); rc=(1-Y1).*760.*Ficmezcla./((1-X1 rc=(1-Y1).*760.*F icmezcla./((1-X1).*Psat2.*Ficsat) ).*Psat2.*Ficsat); ; datos1=[X1',Fiapuro',Ficpuro',Fiasat',Ficsat',Fiamez datos1=[X1',Fiapuro',Ficpuro',Fias at',Ficsat',Fiamezcla',Ficmezcla']; cla',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 disp('coeficiente s 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,F plot(Y1,Fiamezcla ,Y1,Ficmezcla,Y1,Fiapuro,Y1,Ficpuro iapuro,Y1,Ficpuro) ) ylabel('Fi') xlabel('Y1') title('Coeficientes title('Coeficient es de fugacidad vs Y1') legend('Fi^-Acetona','Fi^-Clorofor legend('Fi^-Aceto na','Fi^-Cloroformo','FiAcetona mo','FiAcetona Puro','FiClorofo Puro','FiCloroformo rmo 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 R(1:16)=1.987;%ca l · 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 minimo=inline('(x 1.*log(r1)+x2.*log(r2))-(x1.*(x2. (r2))-(x1.*(x2.^2.*((z(2)./(R.*T ^2.*((z(2)./(R.*T)).*((exp()).*((exp(z(3).*(z(2)./(R.*T))))./(x1+x2.*(ex z(3).*(z(2)./(R.* T))))./(x1+x2.*(exp(-z(3).*(z(2). p(-z(3).*(z(2)./(R.*T)))))).^2+( /(R.*T)))))).^2+((exp((exp(z(3).*(z(1)./(R.*T)))).*(z(1)./(R.* z(3).*(z(1)./(R.* T)))).*(z(1)./(R.*T)))./(x2+x1.*(ex T)))./(x2+x1.*(exp(p(z(3).*(z(1)./(R.*T))))).^2))+x2.*(x z(3).*(z(1)./(R.* T))))).^2))+x2.*(x1.^2.*((z(1)./(R. 1.^2.*((z(1)./(R.*T)).*((exp(*T)).*((exp(z(3).*(z(1)./(R.*T))))./(x2+x1.*(ex z(3).*(z(1)./(R.* T))))./(x2+x1.*(exp(-z(3).*(z(1). p(-z(3).*(z(1)./(R.*T)))))).^2+( /(R.*T)))))).^2+((exp((exp(-
z(3).*(z(2)./(R.*T)))).*(z(2)./(R.*T)))./(x1+x2.*(ex z(3).*(z(2)./(R.*T)))).*(z(2)./(R.* T)))./(x1+x2.*(exp(p(z(3).*(z(2)./(R.*T))))).^2)))','z', z(3).*(z(2)./(R.* T))))).^2)))','z','R','T','x1','x2' 'R','T','x1','x2','r1','r2'); ,'r1','r2'); options=optimset('LargeScale','on', options=optimset( 'LargeScale','on','Display','off',' 'Display','off','HessUpdate','gill HessUpdate','gillmurray'); murray'); [z,fval]=lsqnonlin(minimo,[0,0,0.3] [z,fval]=lsqnonli n(minimo,[0,0,0.3],[-10000,,[-10000,10000,0.2],[10000,10000,0.47],optio 10000,0.2],[10000 ,10000,0.47],options,R,T,x1,x2,r1,r ns,R,T,x1,x2,r1,r2); 2); rteorica1=exp(x2.^2.*((z(2)./(R.*T)).*((exp(z(3).*(z(2)./(R.*T))))./(x1+x2.*(ex z(3).*(z(2)./(R.* T))))./(x1+x2.*(exp(-z(3).*(z(2). p(-z(3).*(z(2)./(R.*T)))))).^2+( /(R.*T)))))).^2+((exp((exp(z(3).*(z(1)./(R.*T)))).*(z(1)./(R.* z(3).*(z(1)./(R.* T)))).*(z(1)./(R.*T)))./(x2+x1.*(ex T)))./(x2+x1.*(exp(p(z(3).*(z(1)./(R.*T))))).^2)); rteorica2=exp(x1.^2.*((z(1)./(R.*T)).*((exp(z(3).*(z(1)./(R.*T))))./(x2+x1.*(ex z(3).*(z(1)./(R.* T))))./(x2+x1.*(exp(-z(3).*(z(1). p(-z(3).*(z(1)./(R.*T)))))).^2+( /(R.*T)))))).^2+((exp((exp(z(3).*(z(2)./(R.*T)))).*(z(2)./(R.* z(3).*(z(2)./(R.* T)))).*(z(2)./(R.*T)))./(x1+x2.*(ex T)))./(x1+x2.*(exp(p(z(3).*(z(2)./(R.*T))))).^2)); error1=abs((r1-rteorica1)./r1)*100; error2=abs((r2-rteorica2)./r2)*100; datos1=[x1',r1',rteorica1',error1', datos1=[x1',r1',r teorica1',error1',r2',rteorica2',er r2',rteorica2',error2']; ror2']; 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','r1e legend('r1teorica ','r2teorica','r1exp','r2exp',0) xp','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.*ex minimo=inline('P(x1.*Psat1atm.*exp(x2.^2.*((zB(2). p(x2.^2.*((zB(2)./(R.*T)).*((exp(/(R.*T)).*((exp(zB(3).*(zB(2)./(R.*T))))./(x1+x2.*( zB(3).*(zB(2)./(R .*T))))./(x1+x2.*(exp(-zB(3).*(zB exp(-zB(3).*(zB(2)./(R.*T)))))). (2)./(R.*T)))))).^2+((exp(^2+((exp(zB(3).*(zB(1)./(R.*T)))).*(zB(1)./( zB(3).*(zB(1)./(R .*T)))).*(zB(1)./(R.*T)))./(x2+x1.* R.*T)))./(x2+x1.*(exp((exp(zB(3).*(zB(1)./(R.*T))))).^2))+x2.* zB(3).*(zB(1)./(R .*T))))).^2))+x2.*Psat2atm.*exp(x1. Psat2atm.*exp(x1.^2.*((zB(1)./(R.* ^2.*((zB(1)./(R.*T)).*((exp( T)).*((exp( -zB(3).*(zB(1)./(R.*T))))./(x2+x1.* -zB(3).*(zB(1)./( R.*T))))./(x2+x1.*(exp(-zB(3).*(z (exp(-zB(3).*(zB(1)./(R.*T)))))) B(1)./(R.*T)))))).^2+((exp(.^2+((exp(-
zB(3).*(zB(2)./(R.*T)))).*(zB(2)./(R.*T)))./(x1+x2.* zB(3).*(zB(2)./(R.*T)))).*(zB(2)./( R.*T)))./(x1+x2.*(exp((exp(zB(3).*(zB(2)./(R.*T))))).^2)))','z zB(3).*(zB(2)./(R .*T))))).^2)))','zB','R','T','x1',' B','R','T','x1','x2','P','Psat1atm x2','P','Psat1atm','Psat2atm ','Psat2atm '); options=optimset('LargeScale','on', options=optimset( 'LargeScale','on','Display','off',' 'Display','off','HessUpdate','stee HessUpdate','steepdesc'); pdesc'); [zB,fval]=lsqnonlin(minimo,z,[-10000,10000,0.2],[10000,10000,0.47],optio 10000,0.2],[10000 ,10000,0.47],options,R,T,x1,x2,P,Ps ns,R,T,x1,x2,P,Psat1atm,Psat2atm); at1atm,Psat2atm); rteorica1B=exp(x2.^2.*((zB(2)./(R.*T)).*((exp(rteorica1B=exp(x2.^2.*((zB(2)./(R.* T)).*((exp(zB(3).*(zB(2)./(R.*T))))./(x1+x2.*( zB(3).*(zB(2)./(R .*T))))./(x1+x2.*(exp(-zB(3).*(zB exp(-zB(3).*(zB(2)./(R.*T)))))). (2)./(R.*T)))))).^2+((exp(^2+((exp(zB(3).*(zB(1)./(R.*T)))).*(zB(1)./( zB(3).*(zB(1)./(R .*T)))).*(zB(1)./(R.*T)))./(x2+x1.* R.*T)))./(x2+x1.*(exp((exp(zB(3).*(zB(1)./(R.*T))))).^2)); rteorica2B=exp(x1.^2.*((zB(1)./(R.* rteorica2B=exp(x1 .^2.*((zB(1)./(R.*T)).*((exp(T)).*((exp(zB(3).*(zB(1)./(R.*T))))./(x2+x1.*( zB(3).*(zB(1)./(R .*T))))./(x2+x1.*(exp(-zB(3).*(zB exp(-zB(3).*(zB(1)./(R.*T)))))). (1)./(R.*T)))))).^2+((exp(^2+((exp(zB(3).*(zB(2)./(R.*T)))).*(zB(2)./( zB(3).*(zB(2)./(R .*T)))).*(zB(2)./(R.*T)))./(x1+x2.* R.*T)))./(x1+x2.*(exp((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' datos2=[x1',r1',rteorica1B',error1B ',r2',rteorica2B',error2B']; ,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','r1e legend('rteorica1 ','rteorica2','r1exp','r2exp',0) xp','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) VARIABLE(SUSTAN CIA,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,:);% 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,:)) 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)+l lnrC(sst,dat)=1 -JJ(sst,1,dat)+log(JJ(sst,1,dat))og(JJ(sst,1,dat))-5.*q(sst,1).*(1 5.*q(sst,1).*(1JJ(sst,1,dat)./LL(sst,1,dat)+log(JJ JJ(sst,1,dat)./LL (sst,1,dat)+log(JJ(sst,1,dat)./LL(s (sst,1,dat)./LL(sst,1,dat))); st,1,dat))); for k=1:length(Rk) sumatoria(k)=thk(1,k,dat).*(Sk(sst sumatoria(k)=thk( 1,k,dat).*(Sk(sst,k,dat)./nn(1,k,d ,k,dat)./nn(1,k,dat))at))G(sst,k)*log(Sk(sst,k,dat)./nn(1,k,dat));
end lnrR(sst,dat)=q(sst,1).*(1-log( lnrR(sst,dat)=q (sst,1).*(1-log(LL(sst,1,dat)))LL(sst,1,dat)))-sum(sumatoria); sum(sumatoria); runifac(sst,dat)=exp(lnrC(sst,dat runifac(sst,dat )=exp(lnrC(sst,dat)+lnrR(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,: plot(x1,runifac(1, :),x1,runifac(2,:),x1,rexp(1,:),x1 ),x1,rexp(1,:),x1,rexp(2,:)); ,rexp(2,:)); title('rUNIFAC(1) y rUNIFAC(2) vs X1'); xlabel('X1') ylabel('r') xlim([0 1]) grid on legend('rUNIFAC(1)','rUNIFAC(2)','r legend('rUNIFAC(1) ','rUNIFAC(2)','rexp(1)','rexp(2)' exp(1)','rexp(2)',0) ,0) error1=abs((rexp(1,:)-runifac(1,:))./rexp(1,:))*100 error1=abs((rexp(1,:)-runifac(1,:) )./rexp(1,:))*100; ; error2=abs((rexp(2,:)-runifac(2,:) error2=abs((rexp(2 ,:)-runifac(2,:))./rexp(2,:))*100 )./rexp(2,:))*100; ; datos1=[x1',rexp(1,:)',runifac(1,:) datos1=[x1',rexp(1 ,:)',runifac(1,:)',error1',rexp(2, ',error1',rexp(2,:)',runifac(2,:)' :)',runifac(2,:)',error2']; ,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,'ScreenSiz get(0,'ScreenSize'); e'); pos1 = [1/10*scnsize(3)+bdwi [1/10*scnsize(3)+bdwidth,1/7*scnsize(4 dth,1/7*scnsize(4) ) + bdwidth,8/10*scn bdwidth,8/10*scnsize(3) size(3) 2*bdwidth,... 5/7*scnsize(4) - (topbdwidth + bdwidth)];
%---------------Coericientes de actividad experimentales(x %---------------Coericientes experimentales(x1;r1;r2)--------1;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)-----------%---------------Coeficientes -------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)---------%---------------Coeficientes ---------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 %---------------Temperatura 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 temburbujaNRTL(1,:)=[0.1013 0.6683 0.7020 0.7315 0.7605 0.8137 0.8946 0.9433 0.9652]; %---------------Temperatura de Rocio Calculada (Y1;Temperatura de Rocio)-------%---------------Temperatura ------------
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 temrocioNRTL(1,:)=[0.0740 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)---%------Coeficientes 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 %------Coeficiente s 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 %------Temperaturas (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 %------Temperaturas (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 %----------------T emperatura de Burbuja y rocioexperimenta rocioexperimental------------l------------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 %---------------Co mposicion 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','of 'NumberTitle','off','MenuBar','Non f','MenuBar','None', e', 'Name',... 'Coeficientes de actividad con varios modelos','Color',[1,1,1],'position',pos1) subplot(2,2,1)
plot(rexperimentales(1,:),rexperimentales(2,:),rexpe plot(rexperimentales(1,:),rexperime ntales(2,:),rexperimentales(1,:),r rimentales(1,:),rexperimenta experimenta les(3,:),runifac(1,:),runifac(2,:), les(3,:),runifac( 1,:),runifac(2,:),runifac(1,:),runi runifac(1,:),runifac(3,:)) fac(3,:)) title('Coeficientes title('Coeficient es 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,:),rexperime plot(rexperimenta les(1,:),rexperimentales(2,:),rexpe ntales(2,:),rexperimentales(1,:),r rimentales(1,:),rexperimenta experimenta les(3,:),rUNIFACPROII(1,:),rUNIFACP les(3,:),rUNIFACP ROII(1,:),rUNIFACPROII(2,:),rUNIFAC ROII(2,:),rUNIFACPROII(1,:),rUNIFA PROII(1,:),rUNIFACPROII(3,:) CPROII(3,:) ) title('Coeficientes title('Coeficient es 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,:),rexperime plot(rexperimenta les(1,:),rexperimentales(2,:),rexpe ntales(2,:),rexperimentales(1,:),r rimentales(1,:),rexperimenta experimenta les(3,:),rntrl(1,:),rntrl(2,:),rntr les(3,:),rntrl(1, :),rntrl(2,:),rntrl(1,:),rntrl(3,:) l(1,:),rntrl(3,:)) ) title('Coeficientes title('Coeficient es 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,:),rexperime plot(rexperimenta les(1,:),rexperimentales(2,:),rexpe ntales(2,:),rexperimentales(1,:),r rimentales(1,:),rexperimenta experimenta les(3,:),rNRTLPROII(1,:),rNRTLPROII les(3,:),rNRTLPRO II(1,:),rNRTLPROII(2,:),rNRTLPROII( (2,:),rNRTLPROII(1,:),rNRTLPROII(3 1,:),rNRTLPROII(3,:)) ,:)) title('Coeficientes title('Coeficient es 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', experimental', 'r1 NRTL PROII', 'r2 NRTL PROII',0) figure(2) set(gcf, 'NumberTitle','of 'NumberTitle','off','MenuBar','Non f','MenuBar','None', e', 'Name',... 'Puntos de Burbuja y Rocio','Color', Rocio','Color',[1,1,1],'position' [1,1,1],'position',pos1) ,pos1) subplot(2,2,1) plot(temburbujaNRTL(1,:),temburbuja plot(temburbujaNR TL(1,:),temburbujaNRTL(2,:),temroci NRTL(2,:),temrocioNRTL(1,:),temroc oNRTL(1,:),temrocioNRTL(2,:) ioNRTL(2,:) ,X1,TC,Y1,TC) title('Temperatura title('Temperatur a 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 plot(temNRTLPROII (1,:),temNRTLPROII(2,:),temNRTLPROI (2,:),temNRTLPROII(3,:),temNRTLPRO I(3,:),temNRTLPROII(4,:),X1, II(4,:),X1, TC,Y1,TC) title('Temperatura title('Temperatur a 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,:),temUNIFACP plot(temUNIFACPRO II(1,:),temUNIFACPROII(2,:),temUNIF ROII(2,:),temUNIFACPROII(3,:),temU ACPROII(3,:),temUNIFACPROII( NIFACPROII( 4,:),X1,TC,Y1,TC) title('Temperatura title('Temperatur a 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,:),X plot(X1,Y1,X1,tem burbujaNRTL(3,:),XY(1,:),XY(2,:),a, Y(1,:),XY(2,:),a,a) a) title('Grafica Equilibrio L-V (Y1 vs X1)') xlabel('X1') ylabel('T') xlim([0 1]) ylim([0 1]) grid on legend('Experimental', legend('Experimen tal', 'NRTL', 'NRTL PROII',0)
