View source for PengRobinson EOS in FPROPS

Work on this is still going on as a part of GSoC 2010, Project [[User:Ankitml|Ankit]]

Comments and suggestions are welcome
==Algorithm==
===Algorithm From Literature===
*Guess a Pressure.  
*Find roots of volume using PR-EOS, if roots<3 guess P again.
*Min V is the molar liquid volume, and the max V is vapour molar volume.
*get Zl and Zv.
*Get fugacities of both phases.
*Pnew = Pold*fl/fv

===The current Algorithm===
*Calculate PV data using PR_EOS from 0.001 to 6666 cm3/mol.
*Get the range of Pressure in which volume has 3 roots.
*Give this range as lower and upper bounds to a golden search algorithm, minimizing the difference in fugacities.
*To calculate fugacities a non iterative analytical expression was derived for PR-EOS.

===Issues That Forced to Deviate From Literature===
*Find a non iterative way to solve a cubic Equation for volume. The normal iterative one makes the algorithm very slow, as it is required to solve the volume cubic many number of times to arrive at the solution. So right now generating PV data for a range with a constant increment in V. 
*This also helps in guessing the first P in the range of 3 roots, which is required as the first step.
*The golden search algorithm is just experimental over the standard Pnew = Pold*Fl/Fv iteration, which requires a good guess whereas Golden search requires a bracket in which solution lies. It is yet to be tested which one is better.

==Limitations==
===Errors===
*Not working at some temperatures. [It should work on entire range from triple point to Critical Temperature]

[Solution] - Currently trying got find out the reason for this behaviour. - [http://ascendwiki.cheme.cmu.edu/User:Ankitml Ankit ]

*Not predicting Saturation Volumes for Liquids correctly.
[Solution] - Its a problem with Peng Robinson EOS, a modification to PR_EOS was proposed by G. Heyen in the paper, "Liquid and Vapor Properties from a cubic Equation of state", in 2nd International Conference on Phase Equilibria and Fluid Properties, 1980.

===Changes to be done to integrate with ASCEND===
*Data (Pc,Tc,Omega) to be taken from data files and only component name should be required in the code.
*No main function
*ASCEND wrapper
*Code should be in C, remove C++ functions.

=Code=
<source lang="c">
#include <stdio.h>
#include <math.h>
#include <iostream>
using namespace std;

//Global Variables
double Tc,Pc,Omega,R,Tre,b,a,RT,T,P[20000],V[20000],Phil,Phiv,Vl,Vv;
int j;

int noRoots(double Psat); //gives out number of roots (of V) at a particular Pressure
double Pcalc(double T, double V);  // Calculates Pressure at Given T and V using PR-EOS
int PhiCalc(double PsatGuess); //Calculates fugacities at a given Pressure and Stores it in global variables
double phiDiff(double Pguess); //Finds fugacity difference at a particular P, uses function PhiCalc and is required for golden search algorithm
double goldSer(double a,double c,double b); //Golden search algorithm, function to be minimized is phiDiff and a is the lower bound, b is upper bound, c is the center.

int main()
{

	Tc = 304.1282;  // Kelvin
        Pc = 73.773;  // bar
        Omega = 0.22394; // acentric factor
        R = 83.14; //bar[cm3]/k/mol
	T = 260;
	RT = R*T;
        Tre = T/Tc;
	
	double p[20000],v[20000];
	//following for loop calculates P V data till v=6500cm3/mol and stores it in P,V arrays only if P is within bounds of 60 and -160
	j = 0;
        for (int i=0;i<20000;i++)
        {
                v[i] = (i + 0.001)/3;
                p[i] = Pcalc(T,v[i]);
                if(p[i]<60.0 && p[i]>-160.0)
                {
                        P[j] = p[i];
                        V[j] = v[i];
                        j++;
                }
        }

//finding the bounds for P
        double pStart,deltaP;
        pStart = 0.01;
        deltaP = 1.0;
        for(int n = 0;n<3;n++)
        {
                for(double pp = pStart;pp<60.0;pp=pp+deltaP)
                {
                        int nor;
                        nor = noRoots(pp);
                        if (nor==3)
                        {
                                pStart = pp - deltaP;
                                pp = 59.0;
                        }
                }
                deltaP = deltaP/10;
        }
        pStart = pStart + 0.01;

        double pEnd;
        pEnd = pStart;
        deltaP = 1.0;
       for(int n = 0;n<2;n++)
       {
                for(double pp = pEnd;pp<60.0;pp=pp+deltaP)
                {
                        int nor;
                        nor = noRoots(pp);
                        if (nor<3)
                        {
                                pEnd = pp - deltaP;
                                pp = 5000;
                        }
                }
                deltaP = deltaP/10;
        }



	int sandbox;
	double Pa,Pb,Pc,sol,phi,resphi;
	phi = (1 + sqrt(5))/2;
        resphi = 2 - phi;
	Pa = pStart;
	Pb = pEnd;
	Pc = (Pb - Pa)*resphi + Pa;
	sol = goldSer(Pa,Pc,Pb);
	sandbox=PhiCalc(sol);
	cout<<"Molar Volume (L) = "<<Vl<<" cm3/mol"<<'\n';
	cout<<"Molar Volume (V) = "<<Vv<<" cm3/mol"<<'\n';
	cout<<"Difference in fugacities of vapour and liquid = "<<fabs(Phil-Phiv)<<'\n';
	cout<<T<<'\t'<<sol<<'\t'<<Vl<<'\t'<<Vv<<'\t'<<fabs(Phil-Phiv)<<'\n';

	return 0;
}

int noRoots(double Psat)
{
        double Vroot[3];
        Vroot[0] = -99; Vroot[1] = -99; Vroot[2] = -99;
        int rootCount = 0;
        for (int i = 0; i < j-1;i++)
        {
                if((Psat<=P[i+1] && Psat>=P[i]) || (Psat>=P[i+1] && Psat<=P[i]))
                {
                        Vroot[rootCount] = V[i]+(Psat-P[i])*(V[i+1]-V[i])/(P[i+1]-P[i]);
                        rootCount++;
                }
        }
        return rootCount;
}

double Pcalc(double T, double V)
{
        double P,m,tmp;
        b = 0.0778*R*Tc/Pc;
        m = 0.37464+1.54226*Omega-0.26992*pow(Omega,2);
        tmp = 1+m*(1-sqrt(Tre));
        a = 0.45724*R*R*Tc*(Tc/Pc)*pow(tmp,2);
        P = R*T/(V-b) - a/(V*(V+b)+b*(V-b));
        return P;
}

int PhiCalc(double PsatGuess)
{       double Vroot[3],A,B;
        B = PsatGuess*b/RT;
        A = a*PsatGuess/pow(RT,2);
        Vroot[0] = -99; Vroot[1] = -99; Vroot[2] = -99;
        int rootCount = 0;
        for (int i = 0; i < j-1;i++)
        {
                if((PsatGuess<=P[i+1] && PsatGuess>=P[i]) || (PsatGuess>=P[i+1] && PsatGuess<=P[i]))
                {
                        Vroot[rootCount] =V[i]+(PsatGuess-P[i])*(V[i+1]-V[i])/(P[i+1]-P[i]);
                        rootCount++;
                }
        }

        double Zv,Zl;
        double term_a,term_b,term_c;

        Vl = Vroot[0];   //liquid volume is the smaller root
        Vv = Vroot[2];   // vapour volume is the largest one

        Zl = PsatGuess*Vl/(R*T);  //the corresponding Z
        Zv = PsatGuess*Vv/(R*T);

        term_a = (1 + sqrt(2))*B;
        term_b = (1 - sqrt(2))*B;
        term_c = 2*sqrt(2)*B;

        Phil = (Zl-1) - log(Zl-B) - A/term_c*log((Zl+term_a)/(Zl+term_b));
        Phiv = (Zv-1) - log(Zv-B) - A/term_c*log((Zv+term_a)/(Zv+term_b));
	return 0;
}

double phiDiff(double Pguess)
{
	double phiD;
	int sandbox;
	sandbox = PhiCalc(Pguess);
	phiD = fabs(Phil - Phiv);
	return phiD;
}

double goldSer(double a,double c,double b)
{
        double d,phi,resphi;
	phi = (1 + sqrt(5))/2;
        resphi = 2 - phi;
        if (fabs(a-b)<0.000001)
                return ((a+b)/2);

        d = c + resphi*(b - c);

        if (phiDiff(d)<phiDiff(c))
                return goldSer(c,d,b);
        else
                return goldSer(d,c,a);
}

</source>
==Possible Improvements in the code==
==Comparisons==