Relatório ELV

Prévia do material em texto

UNIVERSIDADE DO ESTADO DO RIO DE JANEIRO
INSTITUTO DE QUÍMICA 
Departamento Físico Química 
Termodinâmica Experimental para Engenharia Química 
Prof. Pedro Alijó
Equilíbrio Líquido-Vapor
Grupo 1: Gustavo Anciens
 Thaiane Nolasco
 Thainá Tussini
 Viviane Borges
Realização do experimento: 08/09/2015 e 22/09/2015
Rio de Janeiro
2015
Sumário
INTRODUÇÃO --------------------------------------------------------------------------3
OBJETIVO ------------------------------------------------------------------------------- 4
METODOLOGIA ----------------------------------------------------------------------- 4
MODELAGEM -------------------------------------------------------------------------- 6
RESULTADOS E DISCUSSÃO ---------------------------------------------------10
CONCLUSÃO --------------------------------------------------------------------------16
REFERÊNCIAS BIBLIOGRÁFICAS ---------------------------------------------16 
APÊNDICE I – ROTINA DE CÁLCULOS ---------------------------------------16
Introdução teórica
O equilíbrio líquido-vapor é caracterizado quando a taxa das moléculas que atravessam a interface em um sentido (do líquido para o vapor) e no outro (do vapor para o líquido) se igualam, não ocorrendo variações das propriedades macroscópicas de um determinado sistema em função do tempo decorrido. Então, é possível dizer que todos os potenciais que podem afetar a condição de equilíbrio do sistema são igualados.
As equações de equilíbrio líquido-vapor para um sistema fechado à temperatura e pressão constantes são dadas pela igualdade dos potenciais químicos para todos os componentes em ambas as fases:
 (1.1.)
O modelo mais simples para prever comportamento de sistemas em equilíbrio líquido-vapor é a lei de Raoult. As principais hipóteses para a redução dos cálculos do equilíbrio líquido-vapor para a Lei de Raoult dizem que a fase vapor é um gás ideal (pressões baixas e moderadas) e a fase líquida é uma solução ideal (componentes quimicamente semelhantes). 
O diagrama representativo do equilíbrio líquido-vapor geralmente é apresentado como se segue.
Figura 1: Diagrama geral do equilíbrio líquido-vapor
No diagrama acima a fase vapor encontra-se acima da curva dos pontos de orvalho e a fase líquida abaixo da curva dos pontos de bolha. O ponto de orvalho é determinado quando a primeira gota de líquido se faz presente no sistema. Já o ponto de bolha é determinado quando a primeira bolha de vapor se faz presente. A coexistência líquido-vapor está presente entre as curvas do ponto de orvalho e de bolha. Para uma determinada composição do componente i na fase líquida tem-se uma composição do componente i na fase vapor, e, consequentemente, uma determinada linha de amarração.
Objetivo
Construir o envelope de fases que representa o equilíbrio liquido-vapor (ELV) do sistema Água/Etanol à pressão constante, fazendo uso da curva de calibração do ELV do mesmo sistema, bem como realizar a modelagem matemática que prevê o envelope usando a abordagem ϕ-ϕ e a equação de estado cúbica de Redlich Kwong.
Metodologia
Materiais
Balança analítica de precisão;
Condensador com água de Refluxo;
Ebuliômetro;
Frascos de capacidade de 20 mL;
Manta de aquecimento;
Picnômetro de vidro; 
Pró-pipete;
Pipetas graduadas;
Provetas.
3.2 Reagentes
Etanol comercial
Água destilada
3.3 Procedimento Experimental
	O procedimento consiste em realização de séries de destilações de variadas composições de misturas Água/Etanol em ebuliômetro com sistema de refluxo, em dois dias diferentes.
Em cada um dos dias foi realizada uma série de destilações que se iniciaram com a destilação de um dos componentes puros. Os compostos são introduzidos no ebuliômetro com o auxílio da proveta graduada, todo o sistema é fechado, sendo que em uma das saídas do equipamento é acoplado um termômetro que registrará a temperatura do vapor formado.
	O sistema é aquecido, e então o vapor é direcionado ao sistema de condensação, e posteriormente capturado no outro compartimento do ebuliômetro. Após o sistema apresentar-se completo, o aquecimento, agora constante, e o refluxo são mantidos por mais 30 minutos. Alcançado o tempo determinado, amostras das duas fases (fase líquida e fase vapor – condensado recuperado) são coletadas com o auxílio de pipetas e transferidas para frascos. Em cada ensaio subseqüente, tanto no primeiro quanto no segundo dia, apenas adicionou-se o componente referenciado na tabela ao sistema do ensaio anterior.
	Após resfriadas, as soluções tiveram suas densidades aferidas através do método do picnômetro.
Tabela 1: Descrição das misturas utilizadas.
Modelagem RK
A lei de Raoult, conforme já foi dito, é o modelo termodinâmico mais simples que descreve sistemas em equilíbrio líquido-vapor. Considera as hipóteses de fase vapor como um gás ideal e fase líquida como uma solução ideal, o que restringe sua aplicabilidade apenas a sistemas a pressões baixas a moderadas, compostos por espécies quimicamente semelhantes. A partir desse modelo é possível calcular a temperatura, pressão e composição das fases líquida e vapor em equilíbrio. 
A expressão da lei de Raoult é mostrada logo abaixo:
yi P = Pisat xi ( i = 1,2, ... , N) (4.1.)
Onde, 
xi = fração molar do componente i na fase líquida;
yi = fração molar do componente i na fase vapor;
Pisat = pressão de vapor da espécie i pura na temperatura do sistema.
	A pressão parcial de um componente é dada pelo produto de yi por P, expresso no lado esquerdo da equação (4.1).
	Os critérios de equilíbrio termodinâmido de um sistema líquido-vapor compreendem a igualdade de temperatura, pressão e fugacidade de todos os componentes nas duas fases. Em uma situação de equilíbrio, vale assim a relação de igualdade de fugacidades dos componentes da mistura nas fases do sistema, conforme mostrado a seguir:
 (i = 1, 2,..., N)
Onde, os expoentes v, l e sat indicam a fugacidade do componente i na mistura, na fase vapor, líquida e líquido/vapor saturado, respectivamente. 
	Em sistemas em equilíbrio a pressões baixas a moderadas, e misturas compostas por espécies quimicamente semelhantes, podem ser adotadas duas hipóteses simplificadores: a fase vapor é um gás ideal e a fase líquida é uma solução ideal. 
	Igualando a fugacidade da fase líquida à da fase vapor, tem-se para a situação de equilíbrio:
yi P = Pisat xi ( i = 1,2, ... , N)
 
	Existem algumas abordagens empregadas na previsão de equilíbrio líquido-vapor, dentre as quais está a abordagem ϕ-ϕ. Para tal, pode-se expressar a fugacidade em função de um novo parâmetro, o chamado coeficiente de fugacidade (), definido como:
	Reescrevendo a relação de igualdade de fugacidade em função do coeficiente de fugacidade, tem-se que:
Assim, numa abordagem , a equação de isofugacidade pode ser escrita como:
Dividindo toda a equação por P, tem-se então:
	O cálculo que se pretende efetuar é chamado de bolha T, ou seja, deseja-se calcular a temperatura do ponto de bolha da mistura líquida (temperatura de ebulição), dada a sua composição na fase líquida (xi) e a pressão total do sistema (P = P atmosférica), que são conhecidas. Para isso, escolhemos o modelo de equilíbrio em que a fase gasosa é um gás ideal e a fase líquida uma mistura ideal, conforme visto anteriormente.
O coeficiente de fugacidade de um componente de uma mistura são funções da sua composição, temperatura e pressão. Para a determinação dos coeficientes de fugacidade serão empregadas as equações a seguir:
	O parâmetro Zi é o fator de compressibilidade do componente i na mistura, nas fases líquida e vapor. O fator de compressibilidade pode ser determinado para dada T e P por equações de estado de terceiro grau, chamadas por isso de equações de estado cúbicas. 
	Apesar da grande quantidade de equações disponíveis, as equações cúbicas de estado são,ainda, as mais utilizadas para predições de equilíbrio de fases e cálculo de propriedades em aplicações de engenharia como projeto, simulação e otimização de processos. Trata-se de equações de ordem 3 com relação ao volume, que podem ser facilmente estendidas de substâncias puras para sistemas multicomponente, quando associadas a regras de mistura. 
Nas equações de estado cúbicas são utilizadas as propriedades críticas das substâncias, sendo elas a temperatura crítica (TC), pressão crítica (PC) e o fator acêntrico (ω). Isso as torna genéricas, podendo ser aplicadas para os mais variados tipos de sistemas (polares, apolares, assimétricos e outros). 
	Muitas modificações foram propostas à equação cúbica de estado de Van der Walls (VdW), a primeira grande contribuição à termodinâmica na representação do comportamento PVT de sistemas. Dentre as equações de estado cúbicas, tem-se Redlich e Kwong (RK), Soave-Redlich-Kwong (SRK) e Peng-Robinson (PR). 
	Uma equação cúbica de estado genérica é a fonte da equação de van derWaals e de muitas modificações posteriores disponíveis na literatura, conforme segue:
	Nesta equação, os parâmetros e dependem apenas da equação de estado empregada. Já o termo atrativo (a(T)) e o termo de covolume (b) dependem da substância, e são determinados pelas expressões que seguem:
Onde, 
Tr = Temperatura reduzida, dada por: Tr = TTc
Os demais símbolos que aparecem nas equações podem ser vistos na tabela abaixo:
Tabela 2 – Parâmetros característicos de diferentes equações de estado cúbicas
	No presente trabalho, a equação empregada nos cálculos do ELV será a equação de Redlich e Kwong (RK), em função de Z. Desta forma, a equação cúbica de estado apresenta a seguinte forma:
	As equações complementares necessárias seguem a abaixo:
	Além disso, para o cálculo do fator de compressibilidade a partir da pressão, emprega-se a seguinte equação:
Finalmente, 
	A partir do conjunto de equações exposto e dos parâmetros para a equação RK (tabela 2), é possível implementar o algoritmo do cálculo bolha T, um processo iterativo, com o auxílio do programa Scilab.
Resultados e discussão 
	A partir dos valores de densidade de cada fase, obtidos pelo método do pcinômetro, foram previstos os valores de frações molares referentes à composição de cada amostra recolhida. Para tal, empregou-se a curva de calibração obtida na prática anterior, com um ajuste de terceira ordem. Tal curva é reproduzida no presente relatório, conforme segue:
Gráfico 1: Curva de calibração do ELV com ajuste polinomial de terceiro grau
y = -57,424x3 + 166,52x2 - 164,18x + 55,071
y se refere a fração molar de etanol e x, a densidade.
	Os valores de frações molares de etanol em cada fase foram obtidos a partir das densidades experimentais com o uso da equação polinomial de terceiro grau acima. Os valores de frações molares de água puderam ser determinados em seguida por diferença, já que o somatório das frações molares dos dois componentes deve ser igual a um. Os resultados obtidos, bem como as temperaturas de ebulição alcançadas por cada mistura, estão resumidos na tabela abaixo:
	
	Amostra
	Volume água (mL)
	Volume Etanol (mL)
	Massa Picnômetro Vazio (g)
	Massa Picnômetro com Água (g)
	Massa Picnômetro com Amostra (g)
	Densidade
	Fração Molar de Etanol
(1)
	Fração Molar de Água
(2)
	Temperaturas (°C)
	1º dia
	V1
	150
	0
	17,8413
	23,1416
	23,1022
	0,9926
	0,0121
	0,9879
	99
	
	L1
	150
	0
	17,8413
	23,1416
	23,0768
	0,9878
	0,0280
	0,9720
	99
	
	V2
	150
	30
	17,8413
	23,1416
	22,7797
	0,9317
	0,2114
	0,7886
	92
	
	L2
	150
	30
	17,8413
	23,1416
	23,0879
	0,9899
	0,0210
	0,9790
	92
	
	V3
	150
	60
	17,8413
	23,1416
	22,6095
	0,8996
	0,3297
	0,6703
	86
	
	L3
	150
	60
	17,8413
	23,1416
	23,0635
	0,9853
	0,0362
	0,9638
	86
	2º dia
	V1
	0
	150
	18,9687
	24,4484
	23,4098
	0,8105
	0,8181
	0,1819
	78
	
	L1
	0
	150
	18,9687
	24,4484
	23,4867
	0,8245
	0,7191
	0,2809
	78
	
	V2
	10
	150
	18,9687
	24,4484
	23,5904
	0,8434
	0,6007
	0,3993
	80
	
	L2
	10
	150
	18,9687
	24,4484
	23,7493
	0,8724
	0,4479
	0,5521
	80
	
	V3
	15
	150
	18,9687
	24,4484
	23,7984
	0,8814
	0,4066
	0,5934
	81
	
	L3
	15
	150
	18,9687
	24,4484
	24,1782
	0,9507
	0,1481
	0,8519
	81
Tabela 1: Dados experimentais obtidos na prática
	A partir da determinação das composições das fases em equilíbrio termodinâmico e das temperaturas, foi possível construir o envelope de fases do equilíbrio líquido-vapor do sistema etanol (1) /água (2), conforme segue.
Gráfico 2: Envelope de fases experimental
	Os resultados experimentais foram comparados com dados retirados da literatura. Para tal, utilizaram-se os dados de RIEOER (1949) para o sistema em estudo. Os dados encontrados seguem abaixo:
	T (°C)
	x1
	y1
	99,3
	0,0028
	0,032
	96,9
	0,0118
	0,113
	96
	0,0137
	0,157
	96
	0,01104
	0,135
	95,6
	0,0176
	0,156
	94,8
	0,0222
	0,186
	93,8
	0,02106
	0,212
	93,5
	0,0302
	0,231
	92,9
	0,0331
	0,248
	90,5
	0,0519
	0,318
	90,5
	0,053
	0,314
	89,4
	0,0625
	0,339
	88,4
	0,0673
	0,37
	88,6
	0,0115
	0,362
	87,2
	0,0871
	0,406
	85,4
	0,126
	0,468
	84,5
	0,143
	0,487
	84
	0,172
	0,505
	83,4
	0,206
	0,53
	83
	0,21
	0,527
	82,3
	0,255
	0,552
	82
	0,284
	0,567
	81,4
	0,321
	0,586
	81,5
	0,324
	0,586
	81,2
	0,345
	0,591
	80,9
	0,405
	0,614
	80,5
	0,43
	0,626
	80,2
	0,449
	0,633
	80
	0,506
	0,661
	79,5
	0,545
	0,673
	78,8
	0,663
	0,733
	78,5
	0,735
	0,776
	78,4
	0,804
	0,815
	78,3
	0,917
	0,906
Tabela 2: Dados da literatura
	A partir dos dados obtidos, pôde-se gerar o envelope de fases correspondente, conforme segue:
Gráfico 3: Envelope de fases para dados da literatura
	O gráfico a seguir mostra uma comparação entre as curvas obtidas a partir de dados experimentais e de dados da literatura. 
Gráfico 4: Comparação entre o envelope de fases experimental e o da literatura
	A observação do gráfico acima revela que há uma grande discrepância entre o envelope experimental e o da literatura, para valores de fração molar de etanol inferiores a 0,4. A diferença pode ser justificada por erros experimentais, tais como desvios nos valores de temperatura aferidos pelo termômetro ou ainda nos valores da massa fornecidos pela balança. 
	Os valores de composição da fase vapor e temperatura de ebulição obtidos por meio da modelagem empregando a abordagem ϕ-ϕ e a equação de estado cúbica de Redlich e Kwong (RK) são apresentados na tabela abaixo:
	Temperatura (K)
	y1
	Temperatura (K)
	y1
	341.28587
	0.
	310.31378
	0.7799786
	329.30004
	0.3826821
	309.96697
	0.7885513
	322.46767
	0.5417511
	309.61803
	0.7972819
	318.8553
	0.6119170
	309.26761
	0.8061500
	316.97117
	0.6453199
	308.91625
	0.8151368
	315.97865
	0.6622314
	308.56444
	0.8242249
	315.44209
	0.6712737
	308.2126
	0.8333980
	315.13448
	0.6764838
	307.86109
	0.8426407
	314.93578
	0.6798971
	307.51024
	0.8519391
	314.78175
	0.6825940
	307.16032
	0.8612801
	314.6381
	0.6851647
	306.81155
	0.8706515
	314.48689
	0.6879336
	306.46416
	0.880042
	314.31927
	0.6910753
	306.11831
	0.8894410
	314.13152
	0.694677
	305.77415
	0.898839
	313.92275
	0.6987748
	305.43181
	0.9082270
	313.6936
	0.7033743
	305.0914
	0.9175965
	313.44554
	0.7084627
	304.75302
	0.9269399
	313.18038
	0.7140173
	304.41674
	0.9362501
	312.90006
	0.7200095
	304.08263
	0.9455205
	312.60652
	0.7264074
	303.75074
	0.9547449
	312.3015
	0.7331803
	303.42111
	0.9639180
	311.98671
	0.7402958
	303.09378
	0.9730345
	311.66365
	0.7477228
	302.76878
	0.9820898
	311.33369
	0.7554318
	302.44612
	0.9910796
	310.998030.7633948
	302.12582
	1.
	310.65774
	0.7715854
	
	
Tabela 3: Resultados obtidos a partir da modelagem usando RK
	
	A partir dos resultados obtidos, foi possível construir o envelope de fases correspondente ao procedimento de modelagem, mostrado a seguir:
Gráfico 5: Envelope de fases segundo a modelagem RK
	Nota-se que o envelope de fases obtido via modelagem difere tanto daquele referente aos dados da literatura quanto do obtido experimentalmente, não sendo possível ainda observar o azeótropo. Isso indica que pode ter havido alguma falha em relação ao procedimento de modelagem ou ainda que a abordagem ϕ-ϕ usando a equação de RK pode ser inadequada na predição das composições. 
Conclusão
Os dados experimentais obtidos nas práticas oferecem resultados satisfatórios uma vez que o objetivo de montar o envelope a partir desses dados foi atingido. Na comparação do gráfico do envelope de fases experimental com o da literatura também se obteve resultados bem satisfatórios visto que ambos os gráficos ficaram semelhantes para o sistema em questão. Já se referindo aos resultados obtidos pela modelagem RK não foi possível montar um diagrama de fases com o mesmo comportamento que o obtido experimentalmente e na literatura sendo que não foi possível observar o azeótropo. De modo a concluir que a rotina de programação não foi realizada da melhor forma ou ainda que o modelo de RK apresenta certa restrição para o sistema em estudo.
Referências Bibliográficas
SMITH, J.M; VAN NESS. Introdução a Termodinâmica da Engenharia Química. Editora LTC, 7 ed, 2007.
RIEOER R.H.;THOMPSON A.R. IND.ENG.CHEM. 41, 2905, 1949.
Apêndice I – Rotina de Cálculos
//EQUILÍBRIO LÍQUIDO-VAPOR (Redlich/Kwong - RK) 
//Componentes: 1 = etanol 2 = água
clear;clc
exec('func RK.sci');
//DADOS DE ENTRADA
R= 83.14472; 
P= 1.01325; 
//Parâmetros do modelo RK
sigma=1;
e=0;
omega=0.08664;
psi=0.42748;
//Propriedades das espécies puras
Teb1=351.52; 
Teb2=373.15; 
Tc1=513.9; 
Tc2=647.1; 
Pc1=61.48;
Pc2=220.55;
 Kij=-0.12;
// Agua
A2 = 16.8958;
B2 = 3795.17;
C2 = 230.918;
// EtOH
A1 = 16.3872;
B1 = 3885.70;
C1 = 230.170;
//Parâmetros a e b 
b1=omega*(R*Tc1)/Pc1;
b2=omega*(R*Tc2)/Pc2;
ac1=psi*(R*Tc1)^2/Pc1;
ac2=psi*(R*Tc2)^2/Pc2;
ac12=sqrt(ac1*ac2)*(1-Kij);
//i=1;
x=[0.0:1/50:1];
for i=1:51; 
x1=x(i);
x2=1-x1; //Fração de água na fase líquida
T=x1*Teb1+x2*Teb2; //Chute inicial (Kelvin)
y1=0.2; //Chute inicial
y2=1-y1;
Teste1=0; // Condição inicial para verificar se o somatório de Kixi é igual a 1
bL=x1*b1+x2*b2
Zll=1.05*(bL*P)/(R*T)
Zl=1e-3; //Chute inicial
Zv=0.98;
Zvv=0.98; //Chute inicial
tol=1e-8; // Tolerância
//disp(Zll,'chute Zliq')
//*****Início do loop
 //Cálculos para a fração líquida
 k = 0
 
 while abs(Teste1-1)>tol do
 k = k + 1
 
 if k == 200 then
 disp('Não converge o loop externo');
 abort
 end 
 
 [Zl, betaliq, bliq, a1, a2, a12, qliq, aliq] = Zliq(x1, x2, T, P, Zll);
 
 [philiq1,philiq2]=philiq(x1, x2, Zl, betaliq, a1, a2, a12, qliq, aliq); 
// disp(philiq1,'philiq1');
 //Cálculos para a fração vapor
 
 [Zv, betavap, bvap, a1, a2, a12, qvap, avap] = Zvap(y1, y2, T, P, Zvv);
 //disp(Zvap);
 [phivap1,phivap2]=phivap(y1, y2, Zv, betavap, a1, a2, a12, qvap, avap);
 
 
 
 K1=philiq1/phivap1; 
 K2=philiq2/phivap2;
 // disp(philiq1,philiq2,'phiL')
 // disp(phivap1,phivap2,'phivap')
// disp(K1,K2,'ctes. eq.') 
// pause
 Teste2=K1*x1+K2*x2; //Teste 2 e Teste 3 determinam se ouve variação com o novo cálculo de y
 Teste3=1;
 
 kk = 0
 while abs(Teste2-Teste3)>tol do
 kk = kk + 1
 if kk == 100 then
 disp('Não converge loop interno')
 abort
 end 
 Teste3=Teste2; 
 
 y1=K1*x1/(K1*x1+K2*x2); 
 y2=1-y1;
 
 //Cálculos para da nova fração vapor
 
 [Zv, betavap, bvap, a1, a2, a12, qvap, avap] = Zvap(y1, y2, T, P, Zvv);
 //disp(Zvap);
 [phivap1,phivap2]=phivap(y1, y2, Zv, betavap, a1, a2, a12, qvap, avap);
 
 K1=philiq1/phivap1; 
 K2=philiq2/phivap2;
 
 Teste2=K1*x1+K2*x2;
 
end
Teste1=Teste2;
 
 
Tvelho=T;
//y1=K1*x1;
//y2=1-y1;
T=Tvelho-10*(Teste2-1);
 
 // disp([Tvelho, T, Teste2])
 
end
 
 y(i)=y1;
 Tcalc(i)=T;
end
 disp(Tcalc,'Temperatura')
 disp(y,'y1')
 disp(x1,'x1')
plot(x,Tcalc), plot(y,Tcalc)
Functions utilizadas 
function [Zl, betaliq, bliq, a1, a2, a12, qliq, aliq]=Zliq(x1, x2, T, P, chuteZll)
 Tr1=T/Tc1; 
 Tr2=T/Tc2; 
 alfa1=(Tr1)^(-0.5);
 alfa2=(Tr2)^(-0.5);
 a1=ac1*alfa1;
 a2=ac2*alfa2;
 a12 = ac12*sqrt(alfa1*alfa2)
 Kij=-0.09;
 aliq=x1^2*a1+2*x1*x2*a12+x2^2*a2;
 bliq=x1*b1+x2*b2; 
 qliq=aliq/(bliq*R*T); 
 betaliq=(bliq*P)/(R*T); 
 Zll = chuteZll; Zl = Zll + 10*tol;
 ki=0
 while abs(Zl-Zll)>tol do
 ki == ki + 1
 if ki == 100 then
 disp('Problema de convergencia de Zl')
 abort
 end
 
 Zll=Zl
 Zl=betaliq+(Zll+e*betaliq)*(Zll+betaliq*sigma)*((1+betaliq-Zll)/(qliq*betaliq)); 
 end
endfunction
function [Zv, betavap, bvap, a1, a2, a12, qvap, avap]=Zvap(x1, x2, T, P, chuteZvv)
 Tr1=T/Tc1;
 Tr2=T/Tc2;
 alfa1=(Tr1)^(-0.5);
 alfa2=(Tr2)^(-0.5);
 a1=ac1*alfa1;
 a2=ac2*alfa2;
 a12=ac12*sqrt(alfa1*alfa2);
 avap=x1^2*a1+2*x1*x2*a12+x2^2*a2; 
 bvap=x1*b1+x2*b2;
 qvap=avap/(bvap*R*T);
 betavap=(bvap*P)/(R*T); 
 Zvv = chuteZvv; Zv = Zvv + 10*tol;
 ki=0
 while abs(Zv-Zvv)>tol do
 ki == ki + 1
 if ki == 100 then
 disp('Problema de convergencia de Zl')
 abort
 end
 Zvv=Zv
 Zv= 1+betavap-qvap*betavap*((Zvv-betavap)/((Zvv+e*betavap)*(Zvv+sigma*betavap))) 
 end
endfunction
function [philiq1, philiq2]=philiq(x1, x2, Zl, betaliq, a1, a2, a12, qliq, aliq)
 
 Il=(1/(sigma-e))*log((Zl+sigma*betaliq)/(Zl+e*betaliq)); 
 
 //Cálculo de phi do líquido
 qliqmedio1=qliq*[(2*x1*a1+2*x2*a12)/aliq-b1/bliq]; 
 qliqmedio2=qliq*[(2*x2*a2+2*x1*a12)/aliq-b2/bliq]; 
 
 philiq1=exp((b1/bliq)*(Zl-1)-log(Zl-betaliq)-qliqmedio1*Il); 
 philiq2=exp((b2/bliq)*(Zl-1)-log(Zl-betaliq)-qliqmedio2*Il); 
endfunction
function [phivap1, phivap2]=phivap(y1, y2, Zv, betavap, a1, a2, a12, qvap, avap)
 
 Iv=(1/(sigma-e))*log((Zv+sigma*betavap)/(Zv+e*betavap)); 
 
 //Cálculo de phi do líquido
 qvapmedio1=qvap*[(2*y1*a1+2*y2*a12)/avap-b1/bvap]; 
 qvapmedio2=qvap*[(2*y2*a2+2*y1*a12)/avap-b2/bvap]; 
 
 phivap1=exp((b1/bvap)*(Zv-1)-log(Zv-betavap)-qvapmedio1*Iv); 
 phivap2=exp((b2/bvap)*(Zv-1)-log(Zv-betavap)-qvapmedio2*Iv); 
endfunction