Relatório ELLT

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ó
Diagrama de fases - Binodal do Sistema ternário
Grupo 1: Gustavo Anciens
 Thaiane Nolasco
 Thainá Tussini
 Viviane Borges
Realização do experimento: 03/09/2015
Entrega: 24/09/2015
Rio de Janeiro
2015
Sumário
INTRODUÇÃO -------------------------------------------------------------------------- 2
OBJETIVO ------------------------------------------------------------------------------- 4
METODOLOGIA ------------------------------------------------------------------------ 4
MODELAGEM -------------------------------------------------------------------------- 5
RESULTADOS E DISCUSSÃO ---------------------------------------------------- 8
CONCLUSÃO -------------------------------------------------------------------------- 15
REFERÊNCIAS BIBLIOGRÁFICAS ---------------------------------------------- 16
APÊNDICE I – ROTINA DE CÁLCULOS ---------------------------------------- 16
Introdução teórica
Em um sistema com mais de duas substâncias, pode ocorrer a formação de um sistema heterogêneo, ou seja, o aparecimento de mais de uma fase com diferentes composições. Isso ocorre porque o sistema polifásico é mais estável do que o monofásico. Quando as diferentes fases de um sistema encontram-se no equilíbrio termodinâmico, referimo-nos ao equilíbrio líquido-líquido. 
O equilíbrio líquido-líquido pode ser observado em um sistema constituído por dois líquidos parcialmente miscíveis a uma determinada temperatura T. À medida que a temperatura é alterada, a miscibilidade entre os dois componentes muda até o ponto em que o sistema apresenta transição de uma única fase para duas. No caso de um sistema agitado, esta separação pode ser facilmente percebida pela ocorrência de turbidez na mistura. Porém em uma solução de dois líquidos imiscíveis, a adição de um terceiro componente miscível nos dois pode alterar significativamente a solubilidade.
Todos os sistemas buscam o equilíbrio termodinâmico. Um critério de estabilidade deve ser satisfeito, estabelecendo que a uma temperatura e pressão constantes, um estado estável é aquele que apresenta um mínimo na energia livre de Gibbs. Em uma mistura com duas ou mais substâncias, define-se dG como a diferença entre a energia livre de Gibbs da solução e a dos compostos puros. Caso o dG ≤ 0, forma-se uma solução monofásica estável, porém se dG=0, a solução é instável e o sistema é então composto por duas ou mais fases, a fim de minimizar a energia livre de Gibbs. Assim há a formação de sistemas bifásicos ou até multifásicos, e o critério de equilíbrio é representado por:
fi= fi
Onde fi  é a fugacidade da espécie i em solução na fase alfa e fi  é a fugacidade da espécie i em solução na fase beta. Introduzindo os coeficientes de atividade em cada fase no critério de equilibro, e considerando que cada espécie pura existe como líquido na temperatura do sistema, fi= fi= fi, resultamos em:
xi  i  = xi   i 
Em um sistema de equilíbrio líquido-líquido ternário e bifásico, de acordo com a regra das fases, é necessária a especificação prévia de três variáveis
intensivas (pressão, temperatura e uma outra a se determinar).
Para a representação deste sistema, utiliza-se um diagrama triangular no qual os vértices do triângulo representam os componentes puros e os lados correspondem às misturas binárias dos componentes relativos aos respectivos vértices. Para a determinação da composição do sistema em um determinado ponto no diagrama, são traçadas três retas paralelas as faces do triângulo originárias deste ponto, sendo o valor obtido em cada face a fração molar ou mássica de cada componente, como ilustrado na figura 1.
Figura 1: Diagrama de fases triangular
Existem faixas de composição onde se observa a presença de uma única fase, enquanto que há outras faixas em que se observa duas fases. A linha no diagrama triangular que separa essas regiões é denominada de curva binodal ou curva de solubilidade. Como o sistema estudado (ciclohexano-água-etanol) forma um par de líquidos parcialmente miscíveis, a curva binodal usada para representar este sistema é do tipo 1, como ilustrado na figura 2.
Figura 2: diagrama ternário do tipo 1
A curva DLMNE representa a binodal, e em seu interior existem linhas que conectam as composições de equilíbrio de duas fases (linhas de amarração).
Objetivo
Modelar o sistema em equilíbrio de uma solução ternária composta de água, ciclohexano e etanol. Através do modelo UNIFAC, encontrar os coeficientes de atividade e calcular a s composições, além de construir o diagrama de solubilidade do sistema.
Metodologia
3.1 Materiais
1 bureta de 50mL
1 suporte universal
2 pipetas graduada de 1 mL
1 pipeta graduada de 20 mL
2 pêras
2 erlenmeyers de 125 mL
2 pipetas graduada de 5 mL
3 béqueres de 50 mL
3.2 Reagentes
Etanol 95% (álcool comercial)
Água destilada
Ciclohexano
3.3 Procedimento Experimental
	Este experimento foi realizado em dois erlenmeyers:
	Erlnmeyer 1: Esta prática baseia-se em uma titulação que para o primeiro ponto foi adicionado 15 mL de ciclohexano, 5 mL de etanol comercial 95% e nenhuma quantidade de água destilada. Inicia-se a titulação da solução adicionando etanol comercial 95% com uma bureta graduada ao mesmo tempo em que o béquer com a solução é agitado manualmente, até o sistema sair do equilíbrio (solução turva). Anotou-se o volume de etanol titulado. Em seguida foi adicionado diferentes volumes de água destilada, até que seja observado o sistema em equilíbrio (a solução passava de turva para incolor).
	
A tabela a seguir mostra os volumes de água adicionados.
Tabela 3.1 - Resultados do 1º erlenmeyer
	Erlenmeyer 1
	Água
	Ciclohexano
	Etanol comercial
(95% ºGL, 92,8 ºINPM)
	Ponto experimental
	Volume Total (mL)
	Volume Total (mL)
	Volume Inicial (mL)
	Volume Titulado (mL)
	1
	0
	15
	5
	11,3
	2
	0,2
	15
	0
	16
	3
	0,7
	15
	0
	23,7
	4
	3,6
	15
	0
	48,7
	Erlenmeyer 2: Da mesma forma que no Erlenmeyer 1 realizou-se a titulação, porém neste, foi adicionado inicialmente 4 mL de água, 0,3 mL de ciclohexano e 5 mL de etanol. A este sistema, foi adicionado ciclohexano variando o volume e titulando novamente em etanol, anotando-se os resultados, até que seja observado o sistema em equilíbrio como mostrado na tabela abaixo:
Tabela 3.2 - Resultados do 2º erlenmeyer.
	Erlenmeyer 2
	Água
	Ciclohexano
	Etanol comercial
(95% ºGL, 92,8 ºINPM)
	Ponto experimental
	Volume Total (mL)
	Volume Total (mL)
	Volume Inicial (mL)
	Volume Titulado (mL)
	5
	4
	0,3
	5
	1,8
	6
	4
	1,0
	0
	14,4
	7
	4
	3,5
	0
	24
Modelagem
Modelagem UNIFAC
O modelo termodinâmico UNIFAC calcula o coeficiente de atividade das misturas ternárias relacionando as interações entre os grupos estruturais do sistema, ou seja, utiliza a contribuição de grupos onde uma mistura pode ser considerada uma solução de unidades estruturais (grupos fundamentais); além disso, considera que o coeficiente de atividade dos componentes na mistura seja dado por uma contribuição combinatorial (levando em conta as diferenças de tamanho e forma das moléculas), e por uma contribuição residual (levando em conta as interações moleculares). 
As equações referentes ao modelo estão representadas a seguir:
 (5)
 
 (6)
 (7)
 (8)
 
 (9)
Parâmetro de volume para o componente i:(10) 
Onde é o número de grupos do tipo k na molécula i. 
Parâmetro de área superficial para o componente i:
 
 (11) 
 (12) 
 (13)
 (14)
 (15)
 (16)
Algoritmo 
A estratégia escolhida para o cálculo das composições do sistema que compõe a curva binodal foi a do flash bifásico líquido-líquido que irá utilizar as equações descritas abaixo e as equações do modelo UNIFAC para obter um algoritmo de cálculo. A modelagem foi realizada efetuando-se os cálculos do algoritmo no programa Scilab. 
 (17)
 (18) 
 (19) 
 (20)
Calcular 
γ
i
α
, 
γ
i
β
 pelo modelo (UNIFAC).
Calcular K
i
 pela equação 
(17)
.
Encontrar novo 
α
 através da 
equação de flash bifásico líquido-líquido 
(18)
 usando o
 Método de 
Newton –
Raphson
.
Calcular os novos 
x
i
α
, 
x
i
β
 pelas equações 
(19)
 e 
(20)
.
Comparar os novos 
x
i
α
, 
x
i
β
 encontrados com os antigos.
Dados iniciais
:
 
T, 
z
i
 
Chutes iniciais:
 
x
i
α
, 
x
i
β
 e 
α
Σ| x
i 
- 
x
io
 
| < ε
NÃO
S
IM
Calcular novo 
z
i
 = (
x
i
α
 
+
 
x
i
β
) / 2 
Prosseguir com 
z
i
 adicionando etanol
 e 
usando
 
z
i
 de convergência anterior como estimativa inicial 
|
z
i
 – 
x
i
β
| < ε
NÃO
S
IM	
FIM
Resultados e Discussão
Conforme os resultados obtidos experimentalmente e as respectivas propriedades dos componentes, foi possível calcular as composições molares globais (z).
Tabela 5-1 Propriedaddes do Componentes 
	
	Água
	Ciclohexano
	Etanol
	Massa Molar (g/mol)
	18,015
	84,160
	46,068
	Massa Específica (g/mL)
	1,000
	0,779
	0,789
Tabela 5-2 Composições molares globais no Erlenmeyer 1
	
	ÁGUA
	CICLOHEXANO
	ETANOL
	Solução
	Vol.Total (mL)
	Qtd de Mols
	Fração Molar (z)
	Vol.Total (mL)
	Qtd de Mols
	Fração Molar (z)
	Vol.Total corrigido (mL)
	Qtd de Mols
	Fração Molar (z)
	1
	0,0
	0,03136
	0,08858
	15
	0,13884
	0,39214
	10,735
	0,18386
	0,51928
	2
	0,20
	0,05551
	0,12208
	15
	0,13884
	0,30536
	15,200
	0,26033
	0,57255
	3
	0,70
	0,10464
	0,16633
	15
	0,13884
	0,22070
	22,515
	0,38561
	0,61297
	4
	3,6
	0,33500
	0,26457
	15
	0,13884
	0,10965
	46,265
	0,79237
	0,62578
Tabela 5-3 Composições molares globais no Erlenmeyer 2
	
	ÁGUA
	CICLOHEXANO
	ETANOL
	Solução
	Vol.Total (mL)
	Qtd de Mols
	Fração Molar (z)
	Vol.Total (mL)
	Qtd de Mols
	Fração Molar (z)
	Vol.Total corrigido (mL)
	Qtd de Mols
	Fração Molar (z)
	5
	4,0000
	0,2270
	0,8762
	0,3000
	0,0028
	0,0107
	1,7100
	0,0293
	0,1130
	6
	4,0000
	0,2620
	0,5182
	1,0000
	0,0093
	0,0183
	13,6800
	0,2343
	0,4634
	7
	4,0000
	0,2886
	0,4057
	3,5000
	0,0324
	0,0455
	22,8000
	0,3905
	0,5488
A partir desses dados, é possível plotar um diagrama ternário com as frações molares dos componentes:
Gráfico 1: Diagrama ternário dos pontos experimentais x pontos da literatura
No gráfico 1 obteve-se também a plotagem dos dados do artigo PLACKOV, D., STERN, I. 1992. Liquid-liquid equilibria for ternary systems of cyclohexane-water and C1 to C3 alcohols: data and predictions.FluidPhaseEquilibria, 71: 189-209, Para a comparação dos pontos experimentais com os da literatura. Os dados da literatura utilizados seguem abaixo:
Figura 3- Dados da literatura de Plackov
Observa-se que no Gráfico 1 a curva experimental ficou bem próxima da curva da literatura. 
E para a construção do diagrama (Gráfico 2) com os parâmetros utilizados na modelagem para o cálculo dos gamas foram retirados do livro Van Ness S; INTRODUÇÃO A TERMODINÂMICA DA ENGENHARIA QUÍMICA. Editora LTC, 7ª edição, 2007. 
Tabela 5-4: Parâmetros para o modelo UNIFAC
	Molécula
	ri
	qi
	e1,i
	e2,i
	e15,i
	e17,i
	Água
	0,92
	1,4
	0
	0
	0
	1
	Etanol
	2,5755
	2,588
	0.327666151
	0.208655332
	0.465928946
	0
	Ciclohexano
	4,0464
	3,24
	0
	1
	0
	0
Tabela 5-5: Parâmetros de interação de grupos para T=298,15K
	k
	a1,k
	a2,k
	a15,k
	a17,k
	1
	0
	0
	156,4
	300
	2
	0
	0
	331,69
	300
	15
	461,56
	461,56
	0
	-779,24
	17
	1318
	1318
	-540,2
	0
Gráfico 2: Diagrama ternário calculado com parâmetros originais de Smith, Van Ness e Abbott, 2007
A modelagem utilizando os parâmetros do livro não foi satisfatória, pois representa o equilíbrio líquido-líquido somente entre a água e o ciclohexano. 
Então, com o objetivo de encontrar o melhor ajuste para representar os dados calculados tentou-se um ajuste empírico na tabela (5-6) sobre os parâmetros de interação de grupos “aij” da tabela (5-5). 
Tabela 5-6: Parâmetros de interação de grupos modificados para T=298,15K
	k
	a1,k
	a2,k
	a15,k
	a17,k
	1
	0
	0
	331,69
	300
	2
	0
	0
	331,69
	300
	15
	1200
	461,56
	0
	90
	17
	800
	800
	30
	0
O diagrama calculado a partir desses parâmetros pode ser observado a seguir:
Gráfico 3: Diagrama calculado com parâmetros de Plackov
A curva binodal, na região com menor fração de água, não fechou pelo fato de ocorrer uma divisão por zero no cálculo do Flash bifásico, e com o objetivo de representar os dados experimentais, foram realizadas diversas modificações para conseguir ajustar, na melhor maneira, a binodal aos resultados calculados. Assim utilizou-se para o cálculo final os parâmetros de interação descritos na tabela abaixo. 
Tabela 5-7: Parâmetros de interação de grupos finais modificados para T=298,15K
	k
	a1,k
	a2,k
	a15,k
	a17,k
	1
	0
	-10
	250
	950
	2
	-5
	0
	450
	1100
	15
	990
	461,56
	0
	150
	17
	1300
	1300
	-540
	0
	
 Corrigindo o critério de parada para o fechamento da curva, pode-se observar no Gráfico 4:
Gráfico 4: Diagrama calculado com parâmetros ajustados
Alguns parâmetros apresentaram uma sensibilidade maior ao ajuste do que outros. Entretanto, fazer uma análise de sensibilidade de cada parâmetro é uma tarefa complexa.
Com os ajustes finais, chegou-se nas composições de cada componente em fase alfa e beta e seus respectivos gamas descritos nas tabelas (5.8) e (5.9); onde x1, x2, e x3 correspondem respectivamente às composições da Água, Etanol e Ciclohexano. 
Com esses dados obtidos foi possível plotar o gráfico ternário com a curva binodal calculada em comparativo com os dados experimentais e os da literatura. Segue abaixo o diagrama de comparação entre o modelado, literatura e o experimental;
Gráfico 5: Comparação entre modelado, literatura e experimental
Verifica-se que o modelo apresentou resultados bem satisfatórios para a curva binodal calculada, ficando bastante semelhante à curva dos pontos da literatura. 
As tabelas a seguir mostram alguns pontos calculados na modelagem e seus respectivos gamas
Tabela 5-8 Composições e gamas calculados para a fase alfa
	Fase alfa
	x1
	x2
	x3
	1
	2
	3
	3.69148E-4
	0.0
	0.999631
	2708.93808
	5.194837
	1.0000063
	0.0011105
	0.0075555
	0.991333
	805.39862
	4.925017
	1.0010208
	0.0033401
	0.0273851
	0.969274
	212.32050
	4.480472
	1.0054526
	0.00801343
	0.0671496
	0.9248368
	64.5648593
	3.861263
	1.0197571
	0.00993050
	0.08393840.9061310
	46.719344
	3.6504923
	1.0277325
	0.01134813
	0.0966778
	0.8919740
	37.851007
	3.5041955
	1.03453360
	0.023109859
	0.2170350
	0.75985505
	10.315738
	2.4981579
	1.13254358
	0.02591
	0.24971
	0.724372
	8.074243
	2.304423
	1.17109
	0.026084
	0.251802
	0.722112
	7.956230
	2.292899
	1.17375
	0.028926
	0.285876
	0.685196
	6.332782
	2.116677
	1.221007
	0.029133
	0.28836
	0.682498
	6.233324
	2.104649
	1.224755
	0.052136
	0.49267
	0.455191
	2.165678
	1.386964
	1.759948
	0.055983
	0.51183
	0.432186
	1.989674
	1.338241
	1.851685
	0.060761
	0.53503
	0.404200
	1.800278
	1.284199
	1.977156
Tabela 5-9 Composições e gamas calculados para a fase beta
	Fase beta
	x1
	x2
	x3
	1
	2
	3
	0.999999
	0.0
	2.71662E-8
	1.00000
	0.827486
	3.679705
	0.894019
	0.105877
	1.03353E-4
	1.000484
	0.351454
	9601.4346
	0.768993
	0.229439
	0.001567
	0.922218
	0.534776
	621.57534
	0.624808
	0.366532
	0.008659
	0.833402
	0.699732
	108.99315
	0.621903
	0.369199
	0.008897
	0.831939
	0.702284
	105.99951
	0.572966
	0.413378
	0.013655
	0.809727
	0.7412494
	68.198147
	0.538863
	0.443226
	0.017910
	0.797118
	0.764346
	51.52338
	0.2990359
	0.617347
	0.083616
	0.797213
	0.878254
	10.291853
	0.2532023
	0.638370
	0.108427
	0.826203
	0.901449
	7.8237093
	0.2074508
	0.651734
	0.140815
	0.875384
	0.931228
	5.93610993
	0.1858527
	0.654338
	0.1598084
	0.909578
	0.949277
	5.19350882
	0.0794608
	0.585063
	0.3354761
	1.439989
	1.165162
	2.3904640
	0.0747177
	0.574361
	0.350920
	1.511164
	1.189702
	2.28288957
	0.0684868
	0.557792
	0.3737208
	1.6264244
	1.2279693
	2.14136863
Conclusão
Comparando-se os dados obtidos na literatura, composições globais do experimento de Plackov, os dados experimentais obtiveram um considerável desvio, este que pode estar associado a erros sistemáticos inerentes ao experimento, porém também podem estar associados às incertezas dos materiais utilizados como pipetas e buretas. 
Com a utilização do modelo de contribuição de grupos UNIFAC, foi possível modelar o comportamento do nosso sistema, mas a modelagem feita com os parâmetros de interação de grupos retirados do livro Smith Van Ness não estava adequado aos pontos experimentais, assim como os parâmetros utilizados em Plackov.
Os parâmetros que melhor descreveram o modelo, foram os ajustados empiricamente, conforme apresentado no gráfico 4.
A modelagem termodinâmica mostra-se coerente ao se analisar as composições de um componente em uma fase e seu respectivo gama calculado. De modo a verificar que à medida que as composições de cada componente, tanto na fase alfa quanto na fase beta, aumentaram (caminhando para se tornarem puros) os gamas tenderam a diminuir até chegar a 1; e à medida que as composições diminuíram (caminhando para a diluição infinita) os gamas tenderam a aumentar de 1 até infinito. Tais valores estão representados na modelagem no apêndice 1.
Referências Bibliográficas
[1] Plackov,D; LIQUID-LIQUID EQUILIBRIA FOR TERNARY SYSTEMS OF CYCLOHEXANE-WATER AND C1 TO C3 ALCHOOLS; 1992;
[2] SMITH, J. M.; VAN NESS, H. C.; ABBOTT, M. M. Introdução a termodinâmica da engenharia química ed.5 Rio de Janeiro: LTC - Livros Técnicos e Científicos Editora S. A., 2000;
 [3] Kackbart, L.M.: EQUILÍBRIO LÍQUIDO-LÍQUIDO DE SISTEMAS CONTENDO FENOL-ÁGUA-SOLVENTE: OBTENÇÃO E MODELAGEM TERMODINÂMICA DE LUCIANA MEIRELES. UNIVERSIDADE FEDERAL DO PARANÁ; 2007.
Apêndice I – Rotina de Cálculos
exec('functions.sce'); 
//***********************************dados**********************************//
T=298.15; x1a=0.005; x2a=0.000; x3a=1-x1a-x2a;
x1b=0.98; x2b=0.000; x3b=1-x1b-x2b; x=[x1a x2a x3a;x1b x2b x3b];
z1=0.5; z2=0.000; z3=1-z2-z1; z=[z1 z2 z3];
//Ajuste Final
a=[0 -10 200 950; //a1,1 a2,1 a15,1 a17,1 
-5 0 450 1100; //a1,2 a2,2 a15,2 a17,2 
990 461.56 0 150; //a1,15 a2,15 a15,15 a17,15 
1300 1300 -540 0]; //a1,17 a2,17 a15,17 a17,17 
//Parâmetros do livro
//a=[0 0 156.4 300; //a1,1 a2,1 a15,1 a17,1
//0 0 156.4 300; //a1,2 a2,2 a15,2 a17,2
//986.5 986.5 0 -229.1; //a1,15 a2,15 a15,15 a17,15
//1318 1318 353.5 0]; //a1,17 a2,17 a15,17 a17,17
//Parâmetros da literatura
//a=[0 0 331.69 300; //a1,1 a2,1 a15,1 a17,1
//0 0 331.69 300; //a1,2 a2,2 a15,2 a17,2
//461.56 461.56 0 -779.24; //a1,15 a2,15 a15,15 a17,15
//1318 1318 -540.2 0]; //a1,17 a2,17 a15,17 a17,17
r=[0.92 2.5755 4.0464]; //r1 r2 r3
q=[1.4 2.588 3.24]; //q1 q2 q3
e=[0 0 0 1; //e1,1 e2,1 e15,1 e17,1
0.327666151 0.208655332 0.463678516 0; //e1,2 e2,2 e15,2 e17,2
0 1 0 0]; //e1,3 e2,3 e15,3 e17,3
//********************************gráfico************************************//
//Base do gráfico
exec('grafico_ternário.sce');
//Pontos da literatura
exec('pontos_literatura.sce');
//Pontos experimentais
exec('pontos_experimentais.sce');
//*******************Cálculo do ELLT usando UNIFAC***************************//
cont3=0; //contador do 3 loop
w=0.002;//avanço do chute da composição global.
e3=1
while e3>0.0017;//criterio de parada do 3°loop
cont3=cont3+1;
e2=1;
while e2>10^(-5);
y = gama(a, r, q, T, x, e);
K1=y(1,2)/y(1,1);
K2=y(2,2)/y(2,1);
K3=y(3,2)/y(3,1);
K=[K1 K2 K3];
A= NR(z,K)
[xa0,xb0]= compos(K,z,A)
e2= erro(xa0,xb0,x) //criterio de parada do 2 loop)
x=[xa0;xb0];
end
z10=(x(1,1)+x(2,1))/2;
z20=(x(1,2)+x(2,2))/2;
z30=(x(1,3)+x(2,3))/2;
z=[z10-w/2 z20+w z30-w/2];//passo de z
e3=norm(z-[x(2,1) x(2,2) x(2,3)]);
resultados=[resultados;[z(1),z(2),z(3),x(1,1),x(1,2),x(1,3),x(2,1),x(2,2),x(2,3),y(1,1),y(2,1),y(3,1),y(1,2),y(2,2),y(3,2),cont3,A]];
//[z(1),z(2),z(3),x1a,x2a,x3a,x1b,x2b,x3b,y1a,y2a,y3a,y1b,y2b,y3b, cont3,A];
//disp(resultados);
x1alfa=resultados(:,4);
x2alfa=resultados(:,5);
x1beta=resultados(:,7);
x2beta=resultados(:,8);
plot2d(x1alfa+1/2*x2alfa,x2alfa,6)/
plot2d(x1beta+1/2*x2beta,x2beta,3)
end
Functions
//----- gR-----//
 function F=gR(a, T, x, q, e) 
 for j=1:4 //τ1,1 τ2,1 τ15,1 τ17,1 
 for i=1:4 //τ1,2 τ2,2 τ15,2 τ17,2 
 p=exp(-a(j,i)/T); //τ1,15 τ2,15 τ15,15 τ17,15 
 tal(j,i)=p; //τ1,17 τ2,17 τ15,17 τ17,17 
 end 
 end 
for j=1:3 
 for i=1:4 
 w=e(j,1)*tal(i,1)+e(j,2)*tal(i,2)+e(j,3)*tal(i,3)+e(j,4)*tal(i,4); 
 Beta(j,i)=w; //β1,1 β1,2 β1,15 β1,17 
 end //β2,1β2,2 β2,15 β2,17 
end //β3,1 β3,2 β3,15 β3,17 
for j=1:2 //anda na linha de x 
 for i=1:4 //anda na coluna do e 
 t=(x(j,1)*q(1)*e(1,i)+x(j,2)*q(2)*e(2,i)+x(j,3)*q(3)*e(3,i))/(x(j,1)*q(1)+x(j,2)*q(2)+x(j,3)*q(3)); 
 teta(i,j)=t; //teta1a teta1b 
 end //teta2a teta2b 
end //teta15a teta15b
 //teta17a teta17b 
for j=1:2 //anda na coluna de teta 
 for i=1:4 //anda na linha de tal 
 S=teta(1,j)*tal(i,1)+teta(2,j)*tal(i,2)+teta(3,j)*tal(i,3)+teta(4,j)*tal(i,4); 
 s(i,j)=S; //s1a s1b 
 end //s2a s2b 
end //s15a s15b
 //s17a s17b 
for j=1:2 
 for i=1:3 
 d=exp(q(i)*(1-((((teta(1,j)*Beta(i,1))/s(1,j))-(e(i,1)*(log(Beta(i,1)/s(1,j))))))-((((teta(2,j)*Beta(i,2))/s(2,j))-(e(i,2)*(log(Beta(i,2)/s(2,j))))))-((((teta(3,j)*Beta(i,3))/s(3,j))-(e(i,3)*(log(Beta(i,3)/s(3,j))))))-((((teta(4,j)*Beta(i,4))/s(4,j))-(e(i,4)*(log(Beta(i,4)/s(4,j))))))));
 yR(i,j)=d; //y1aR y1bR
end //y2aR y2bR 
end //y3aR y3bR 
F = yR
endfunction 
//----gC-----//function F=gC(x, r, q) 
for j=1:2 
 for i=1:3 
 f=r(i)/(r(1)*x(j,1)+r(2)*x(j,2)+r(3)*x(j,3)); 
 J(i,j)=f; // J1a J1b 
 end // J2a J2b 
end // J3a J3b 
for j=1:2 
 for i=1:3 
 l=q(i)/(q(1)*x(j,1)+q(2)*x(j,2)+q(3)*x(j,3)); 
 L(i,j)=l; // L1a L1b 
 end // L2a L2b
end // L3a L3b 
for j=1:2 
 for i=1:3 
 c=exp(1-J(i,j)+log(J(i,j))-5*q(i)*(1-(J(i,j)/L(i,j))+log(J(i,j)/L(i,j)))); 
 yC(i,j)=c; // y1aC y1bC 
 end // y2aC y2bC 
end // y3aC y3bC 
F=yC 
endfunction 
//------gama----//
function y=gama(a, r, q, T, x, e) 
 yR=gR(a, T, x, q, e); 
 yC=gC(x, r, q); 
 y=exp(log(yC)+log(yR)); 
 endfunction 
//----Newton Raphson----//
function A=NR(z, K) 
 A=0.5// chute inicial 
 e1=1; 
 while e1>10^(-5); 
 cont1=cont1+1; 
 f=(z(1)*(K1-1)/(1+A*(K1-1)))+(z(2)*(K2-1)/(1+A*(K2-1)))+(z(3)*(K3-1)/(1+A*(K3-1))); /// Função para calculo do Flash em f(a) 
 df=-((z(1)*(K1-1)^2)/(1+A*(K1-1))^2+(z(2)*(K2-1)^2)/(1+A*(K2-1))^2+(z(3)*(K3-1)^2)/(1+A*(K3-1))^2);//Derivada da função f 
 A1=A-(f/df);// função de iteração de Newton Raphson 
 e1=abs(A1-A); 
 A=A1; 
 end 
 endfunction 
//-----critério de parado do segundo loop---//
function e2=erro(xa0, xb0, x) 
 e2=(abs(xa0(1,1)-x(1,1))+abs(xa0(1,2)-x(1,2))+abs(xa0(1,3)-x(1,3))+abs(xb0(1,1)-x(2,1))+ abs(xb0(1,2)-x(2,2))+abs(xb0(1,3)-x(2,3)));
endfunction 
//----composições alfa e beta---//
function [xa0, xb0]=compos(K, z, A) 
 for i=1:3 
 xa=K(i)*z(i)/(1+A*(K(i)-1)); 
 xb=z(i)/(1+A*(K(i)-1)); 
 xa0(1,i)=xa; 
 xb0(1,i)=xb; 
 end 
 endfunction