Menu
▾
▴
[Dev-C++] Help
|
From: Nancy A. <nai...@ho...> - 2003-08-22 03:25:39
|
/***************************************************************************
CALCULO DEWT
Sistema ternario: acetona(1)/metano(2)/agua(3)
****************************************************************************
Data: 16/08/2003
Calculo DEWT com o metodo de Wilson para os coeficientes de atividade em um
sistema ternario usando o metodo de Newton-Raphson para a solucao aproximada
do sistema de equacoes estabelecido.
ENTRADAS:
Numero de equacoes, n, e de desconhecidos, z=(z1,...,zn)t, tolerancia TOL;
numero
maximo de iteracoes Max; equacoes para F(z)=0 na funcao F() e para J(z) em
jaco-
biana().
Dados do sistema desde o arquivo cpp-dados.txt: pressao (P); fracao molar
dos com-
ponentes na entrada do sistema, (y(0,..,n-1)); volume molar de cada
componente da
mistura, v(0,..,n); constantes de Antoine, A(0,..,n), B(0,..,n) e C(0,..n),
para o
calculo das presoes de saturacao de cada componente; constantes da equacao
de Wil-
son, a(ij), para cada par de componentes.
SAIDA: solucao aproximada z = (z1,...,zn)t ou mensagen de que o numero de
iteracoes
foi excedido.
******************************************************************************/
#include<fstream.h>
#include <iostream.h>
#include <math.h>
#include <stdlib.h>
/*** Prototipo das funcoes usadas ***/
/************************************/
double pressat(double , double *, double* , double *, double *);
double atividade(int ,double **, double **, double *);
double wilson(double *, double **, double &, double **) ;
void msg1();
double valor(int , int , double *, double **, double *, double *, double *);
double derGx(int, int, double **, double *, double **, double**,double **);
double derF (int , int , double *, double **, double **, double *, double
**);
double Jacobiana(double **, double *, double **, int, double**);
double Inversa(double** , int, double** );
/* Declara variaveis globais */
/****************************/
const double R=1.987, /* constante dos gases*/
ff=1e-4;
int n, c; /*número de equacoes e numero de
especies quimicas*/
/***************************************************************************/
int main () { //INICIO DO MAIN
/***************************************************************************/
/**** Declara variaveis locais ****/
int i,j,k,l,m,kk,jj; /* índices */
int Max; /* número máximo de iterações */
double TOL, /* tolerancia */
P, /* Presao do sistema*/
soma, /* auxiliar de calculo;*/
erro; /* erro de aproximacao */
double *y, /* fracao molar de cada especie no
vapor*/
*v, /* volume molar de cada especie */
*A,*B,*C, /* constantes de Antoine*/
*Gamma, /* coeficiente de atividade*/
*Psat, /* presao de saturacao*/
*F, /* vetor de funcoes*/
**Gder, /* termo a direita para a derivada
numerica*/
**Fder, /* termo a direita paara a derivada
numerica*/
**z, /* vetor de incognitas*/
**vji, /* radio vj/vi na correlacao de Wilson*/
**a, /* constante de interacao binaria de
Wilson*/
**LAMBDA, /* termo de interacao molecular da
equacao de Wilson*/
**J, /* matriz jacobiana */
**JINV; /* jacobiana inversa */
ifstream OpenFile; OpenFile.open("cpp_home.txt");
ofstream SaveFile; SaveFile.open("cpp_DEWT.txt");
SaveFile<<endl<<"Sistema ternario: acetona(1)/metano(2)/agua(3)"<<endl;
SaveFile<<endl<<" Calcula T, x1, x2 e x3 no ponto de orvalho"<<endl;
/***************************/
/* ENTRADAS */
/***************************/
cout<<"Entrar numero maximo de iteracoes, Max: ";
cin>>Max;
cout<<endl<<"Entrar a Tolerancia. TOL (1.0e-5)?: ";
cin>>TOL;
cout<<"Lendo numero de equacoes, n: "<<endl;
OpenFile>>n; cout<<n<<endl;
if (n <= 0)
cout<<"ERROR - n tem que ser maior que zero.\n"<<endl;
SaveFile<<"O sistema tem "<<n<<" equacoes e "<<n<<" desconhecidos."<<endl;
/******** Alocacao dinamica da memoria *******/
c=n-1;
y = new double [c]; if (y==NULL)msg1();
v = new double [c]; if (v==NULL)msg1();
A = new double [c]; if (A==NULL)msg1();
B = new double [c]; if (B==NULL)msg1();
C = new double [c]; if (C==NULL)msg1();
Psat = new double [c]; if (Psat==NULL)msg1();
Gamma = new double [c]; if (Gamma==NULL)msg1();
F = new double [n]; if (F==NULL)msg1();
a = new double *[c]; if (a==NULL)msg1(); for (i=0; i<c; i++){
a[i] = new double [c]; if(a[i]==NULL)msg1();}
z = new double *[n]; if (z==NULL)msg1(); for (i=0; i<n; i++){
z[i] = new double [Max]; if(z[i]==NULL)msg1();}
Gder = new double *[c]; if(Gder==NULL)msg1(); for (i=0; i<c; i++){
Gder[i] = new double [n]; if(Gder[i]==NULL)msg1();}
LAMBDA = new double *[c]; if (LAMBDA==NULL)msg1(); for ( i=0; i<c;
i++){
LAMBDA[i] = new double [c]; if (LAMBDA[i]==NULL)msg1();}
Fder = new double *[n]; if (Fder==NULL)msg1(); for (int i=0; i<n;
i++){
Fder[i] = new double [n]; if(Fder[i]==NULL)msg1();}
J = new double *[n]; if (J==NULL)msg1(); for (int i=0; i<n;
i++){
J[i] = new double [n]; if(J[i]==NULL)msg1();}
JINV = new double *[n]; if (JINV==NULL)msg1(); for (int i=0; i<n;
i++){
JINV[i] = new double [n]; if(JINV[i]==NULL)msg1();}
OpenFile >> P;
for (i=0; i<c-1; i++)OpenFile >>y[i];
for (i=0; i<c; i++)OpenFile >>v[i];
for (i=0; i<c; i++)OpenFile >>A[i];
for (i=0; i<c; i++)OpenFile >>B[i];
for (i=0; i<c; i++)OpenFile >>C[i];
OpenFile >> a[0][1]>> a[0][2];
OpenFile >> a[1][0]>> a[1][2];
OpenFile >> a[2][0]>> a[2][1];
//SaveFile<<"pressao do sistema: "<<P<<endl;
//SaveFile<<"fracao molar na fase vapor: "<<endl;
soma=0;
for(i=0; i<c-1; i++)soma+=y[i];
y[c-1]=1-soma;
for(i=0; i<c; i++)
//SaveFile<<" y"<<i<<"="<<y[i]<<endl;
//SaveFile<<"volume molar das especies:"<<endl;
//for (i=0; i<c; i++) SaveFile<<" v"<<i<<"="<<v[i]<<endl;
//SaveFile<<"coeficientes de Antoine:"<<endl;
for (i=0; i<c; i++)
//SaveFile<<" A"<<i<<"="<<A[i]<<", B"<<i<<"="<<B[i]<<",
C"<<i<<"="<<C[i]<<endl;
//SaveFile<<"parametros de ajuste de Wilson"<<endl;
//SaveFile<<"a12="<<a[0][1]<<" a13="<<a[0][2]<<endl;
//SaveFile<<"a21="<<a[1][0]<<" a23="<<a[1][2]<<endl;
//SaveFile<<"a31="<<a[2][0]<<" a32="<<a[2][1]<<endl;
for(i=0; i<c; i++)a[i][i]=0;
k = 0;
cout<<"Entrar as aproximacoes iniciais:"<<endl; /*
z[i][0]*/
for (i=0;i<c;i++) {//1
cout<<"x"<<i<<"0: ";
cin>>z[i][0];
}//1
cout<<"T0:";
cin>>z[c][0];
for (i=0; i<n; i++) {//2
SaveFile<<"z["<<i<<"]="<<z[i][0]<<endl; /*Imprime as aproximacoes
iniciais*/
cout<<endl;
}//2
for(i=0; i<c; i++)a[i][i]=0;
/*******************************************************/
/* PROCEDIMENTO ITERATIVO */
/******************************************************/
while ( k<Max) { //Inicio do while
SaveFile<<"k"<<k+1<<endl;
/* Calculo dos coeficientes de atividade com a correlacao de Wilson */
wilson(v,a,z[c][k],LAMBDA); /* obtem os parametros adjustaveis de
Wilson*/
/*for (i=0; i<c; i++){//3
for (j=0; j<c; j++)
SaveFile<<"L["<<i<<j<<"]="<<LAMBDA[i][j]<<endl;
}//3 */
// SaveFile<<"Coeficientes de Atividade calculadas"<<endl;
atividade(k,z,LAMBDA,Gamma); /* Obtem os coeficiente de atividade*/
//for(i=0; i<c; i++)SaveFile<<"Gamma["<<i<<"]="<<Gamma[i]<<endl;
//SaveFile<<"Presões de saturação calculadas"<<endl;
pressat(z[c][k],A,B,C,Psat); /* Calcula as presoes de saturacao */
//for (i=0; i<c; i++) SaveFile<<" Psat"<<i<<"="<<Psat[i]<<endl;
//SaveFile<<"Vetor F[z]:"<<endl ;
valor(0, k, y, z, Gamma, Psat, F);
valor(1,k,y,z,Gamma, Psat,F);
valor (2,k,y,z, Gamma, Psat, F);
valor(3,k,y,z,Gamma, Psat,F); /* Calcula o vetor F[i] */
// SaveFile<<"F"<<i<<"="<<F[i]<<endl;
/*************************************************************/
/* Calculo dos elementos da Jacobiana */
/*************************************************************/
for(i=0; i<c; i++) /* Calculo valores no ponto
(1+ff)*z */
derGx(i, k, z, v, a, LAMBDA, Gder);
for(i=0; i<c; i++)
derF(i, k, y, z, Gder, Psat, Fder);
derF(c, k, y, z, Gder, Psat, Fder);
Jacobiana(Fder, F, z, k, J); /* Calcula dFz/dz */
/* for(i=0; i<n; i++)
for(j=0; j<n; j++)
SaveFile<<"J"<<i<<j<<"="<<J[i][j]<<endl; */
/*************************************************************/
/* Solucao do sistema linear nxn J(X)*DeltaZ = -F(X) */
/*************************************************************/
Inversa(J, n, JINV) ; /* Inverte a Jacobiana */
/* for(i=0; i<n; i++)
for(j=0; j<n; j++)
SaveFile<<"JINV"<<i<<j<<"="<<JINV[i][j]<<endl;*/
for (i=0; i<n; i++){//5
soma=0;
for(j=0; j<n; j++){soma+=-JINV[i][j]*F[j];
// SaveFile<<soma<<endl;
}
z[i][k+1]=z[i][k]+soma; /*Novas estimativas*/
SaveFile<<"z["<<i<<"]["<<k+1<<"]="<<z[i][k+1]<<endl; /*
imprime a solucao na
k+1
iteracao no arquivo de saida */
} //5
/***********
* SAIDAS *
***********/
erro=0;
for (i=0; i<n; i++){//6
erro+=fabs((z[i][k+1]-z[i][k])/z[i][k]);
// SaveFile<<erro<<endl;
}//6
//SaveFile<<"erro"<<erro<<endl;
if (k+1>=1){//7a
if (erro<= TOL) {//7 /* Checa precisao da nova estimativa */
/* Libera memoria */
delete []y; delete []v; delete []A; delete []B; delete []C; delete
[]Gamma;
for (int i=0; i<c; i++) delete[] a[i]; delete [] a ;
for (int i=0; i<n; i++) delete[] Gder[i]; delete [] Gder ;
for (int i=0; i<n; i++) delete[] z[i]; delete []z;
for (int i=0; i<c; i++) delete []LAMBDA[i]; delete []LAMBDA;
for (int i=0; i<c; i++) delete[] Fder[i]; delete [] Fder ;
for (int i=0; i<n; i++) delete []J[i]; delete []J;
for (int i=0; i<n; i++) delete []JINV[i]; delete []JINV;
exit (1); /* Procedimento completado com
sucesso */
}//7
}//7a
k++; /* k = k + 1 */
}
SaveFile<<endl<<" Numero Maximo de iteracoes foi excedido."<<endl;
/* Libera a memoria */
delete []y; delete []v; delete []A; delete []B; delete []C; delete
[]Gamma;
delete []Psat; delete []F;
for (int i=0; i<c; i++) delete[] a[i]; delete [] a ;
for (int i=0; i<n; i++) delete[] Gder[i]; delete [] Gder ;
for (int i=0; i<n; i++) delete[] z[i]; delete []z;
for (int i=0; i<c; i++) delete []LAMBDA[i]; delete []LAMBDA;
for (int i=0; i<c; i++) delete[] Fder[i]; delete [] Fder ;
for (int i=0; i<n; i++) delete []J[i]; delete []J;
for (int i=0; i<n; i++) delete []JINV[i]; delete []JINV;
system("PAUSE");
return 0;
} /* Parar */
/*******************/
void msg1()
/******************/
{
cout<<"Nao consegui alocar memoria"<<endl;
exit(1);
}
/***********************************************************************/
double wilson(double *v, double **a, double &z, double **LAMBDA)
/***********************************************************************/
{ //INICIA
int i, j, c=3;
double **vji;
vji=new double *[c];
for(i=0; i<c; i++)vji[i]=new double[c];
for (i=0; i<c; i++){//1
for (j=0; j<c; j++){//2
vji[i][j]= v[j]/v[i];
a[i][j]=a[i][j]/1.987;
}//2
} //1
for (i=0; i<c; i++){//3
for (j=0; j<c; j++){//4
if(j==i)LAMBDA[i][j]=1.;
else
LAMBDA[i][j]=vji[i][j]*exp(-a[i][j]/z);
}//4
}//3
}//FINAL
/*******************************************************************/
double atividade(int k, double **z, double **LMBDA, double *Gamma)
/*******************************************************************/
{//INICIO
int i, kk,j, c=3;
double somakkjj, somak, somaj;
for (i=0; i<c; i++){/*1*/
somak=0;
for (kk=0; kk<c; kk++){ /*2*/
somakkjj=0;
for (j=0; j<c; j++)somakkjj+=z[j][k]*LMBDA[kk][j];
somak+= z[kk][k]*LMBDA[kk][i]/somakkjj;
}/*2*/
somaj=0;
for (j=0; j<c; j++) somaj+=z[j][k]*LMBDA[i][j];
Gamma[i]=exp(1-log(somaj)-somak);
}/*1*/
}//FINAL
/********************************************************************************/
double pressat(double z, double *A, double *B, double *C,double *Psat )
/********************************************************************************/
{ //INICIO
int i, n=3;
for(i=0; i<n; i++){//1
Psat[i]=exp(A[i]-B[i]/(z+C[i]));
}//1
}//FINAL
/***********************************************************************************/
double valor(int f, int k, double *y, double **z, double *G, double *Ps,
double *F)
/***********************************************************************************/
{//INICIO
if(f==0) F[0]=y[0]-z[0][k]*G[0]*(Ps[0]/101.303);
if(f==1) F[1]=y[1]-z[1][k]*G[1]*(Ps[1]/101.303);
if(f==2) F[2]=y[2]-z[2][k]*G[2]*(Ps[2]/101.303);
if(f==c) F[c]=z[0][k]+z[1][k]+z[2][k]-1;
}//FINAL
/***********************************************************************************/
double derGx(int f,int k,double **z,double *v,double **a,double **L, double
**Gder)
/***********************************************************************************/
{ //INICIO
int i,j, c=3;
double *Gd, **zp;
Gd = new double[n];
zp=new double*[n];
for(i=0; i<n; i++)zp[i]=new double[k];
for(i=0; i<n; i++)
zp[i][k]=z[i][k];
if(f==0){//0
zp[0][k]=(1+ff)*zp[0][k]; /* cout<<"gder"<<zp[0][k]<<endl;
cout<<"gder"<<zp[1][k]<<endl; cout<<"gder"<<zp[2][k]<<endl;
cout<<"gder"<<zp[3][k]<<endl; */
atividade(k, zp, L, Gd);
for(i=0; i<c; i++)Gder[i][0]=Gd[i];
}//0
//cout<<endl;
if(f==1){//1 // cout<<f<<endl;
zp[1][k]= (1+ff)*zp[1][k]; /* cout<<"gder"<<zp[0][k]<<endl;
cout<<"gder"<<zp[1][k]<<endl; cout<<"gder"<<zp[2][k]<<endl;
cout<<"gder"<<zp[3][k]<<endl;*/
atividade(k, z, L, Gd);
for(i=0; i<c; i++)Gder[i][1]=Gd[i];
}//1
//cout<<endl;
if(f==2){//2 // cout<<f<<endl;
zp[2][k]= (1+ ff)*zp[2][k]; /* cout<<"gder"<<zp[0][k]<<endl;
cout<<"gder"<<zp[1][k]<<endl; cout<<"gder"<<zp[2][k]<<endl;
cout<<"gder"<<zp[3][k]<<endl;*/
atividade(k, z, L, Gd);
for(i=0; i<c; i++)Gder[i][2]=Gd[i];
}//2
// cout<<endl;
if(f==3){//3 //cout<<f<<endl;
zp[3][k]=(1+ff)*zp[3][k]; /*cout<<"gder"<<zp[0][k]<<endl;
cout<<"gder"<<zp[1][k]<<endl; cout<<"gder"<<zp[2][k]<<endl;
cout<<"gder"<<zp[3][k]<<endl;*/
for(i=0; i<c; i++){//a
for(j=0; j<c;j++){//b
L[i][j]=exp((v[j]/v[i])*exp(-a[i][j]/z[3][k]));
}//b
}//a
atividade(k, z, L, Gd);
for(i=0; i<c; i++)Gder[i][3]=Gd[i];
}//3 /*cout<<"gfin"<<z[0][k]<<endl;
/*cout<<"gfin"<<z[1][k]<<endl;
cout<<"gfin"<<z[2][k]<<endl;
cout<<"gfin"<<z[3][k]<<endl;*/
delete[]Gd;
}//FINAL
/**********************************************************************************************/
double derF(int f, int k, double *y, double **z, double **Gder, double *Ps,
double **Fder)
/**********************************************************************************************/
{//INICIO
int i,j, c=3, n=4;
double *zp;
zp=new double[n];
for (i=0; i<c; i++){//1
zp[i]=(1+ff)*z[i][k];
for (j=0; j<n; j++){ //1a
if(i==j)Fder[i][j]=y[i]-Gder[i][i]*(Ps[i]/101.33)*zp[i];
else
Fder[i][j]=y[i]-Gder[i][j]*(Ps[i]/101.33)*z[i][k];
}//1a
}//1
} //FINAL* */
/***************************************************************************/
double Jacobiana(double **Fder, double *F, double **z, int k, double **J)
/***************************************************************************/
{//INICIO
int i,j, n=4;
double a, d;
for (i=0; i<c; i++){//1
a=z[i][k];
d=(1+ff)*z[i][k]; /*cout<<"a="<<a<<"
d="<<d<<endl;*/
cout<<endl;
for (j=0; j<n; j++){//1a
J[i][j]=(Fder[i][j]-F[i])/(d-a);
/*cout<<endl<<"(Fder"<<i<<j<<"="<<Fder[i][j]<<"-F"<<i<<"="<<F[i]<<
"/(d="<<d<<"-a="<<a<<")=
J"<<i<<j<<"="<<J[i][j];endl; */
}//1a
cout<<endl;
}//1
for(i=0; i<n; i++)J[c][i]=1;
J[c][c]=0; // for (i=0; i<n; i++){//1
for (j=0; j<n; j++) cout<<"J"<<i<<j<<"="<<J[i][j]<<endl;}//1
}//FINAL
/********************************************************************/
double Inversa(double **J, int n, double **JINV)
/*******************************************************************/
{//INICIO
double Temp, **One;
int i,j,m,l;
One= new double *[n];
for(i=0; i<n; i++)One[i]=new double[n];
/*** Criando a matriz identidade ***/
for (i=0; i<n; i++){//1
for (j=0; j<n; j++){//2
if(i==j)One[i][j]=1;
else One[i][j]=0;
// cout<<"one="<<One[i][j]<<endl;
}//2
}//1
/*** Inversao da Jacobiana ***/
for (l=0; l<n; l++){//7
cout<<l<<endl;
Temp=J[l][l];
for (j=0; j<n; j++) {//8
J[l][j]=J[l][j]/Temp;
One[l][j]=One[l][j]/Temp;
}//8
for (m=l+1; m<n; m++) {//9
Temp=J[m][l];
for (j=0; j<n; j++){ //10
J[m][j]=-1*Temp*J[l][j]+J[m][j];
One[m][j]=-1*Temp*One[l][j]+One[m][j];
} //10
} //9
} //7
for (l=c; l>=0; l--) { //11
for (m=l-1; m>=0; m--){ //12
Temp=J[m][l];
for (j=0; j<n; j++) {//13
J[m][j]=-1*Temp*J[l][j]+J[m][j];
cout<<endl;
One[m][j]=-1*Temp*One[l][j]+One[m][j];
} //13
} //12
} //11
for (i=0; i<n; i++){
for (j=0; j<n; j++){
JINV[i][j]=One[i][j];
}
}
}//FINAL
|
