PengRobinson EOS in FPROPS

Documentation

Algorithm

The general Algorithm

• Guess a Pressure.
• Find roots of volume using PR-EOS, if roots<3 guess P again.
• Min V is the liquid density at that P, and the max V is vapour.
• 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.

Limitations

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

#include <stdio.h>
#include <math.h>
#include <iostream>
using namespace std;



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.
//Global Variables
double Tc,Pc,Omega,R,Tre,b,a,RT,T,P[20000],V[20000],Phil,Phiv,Vl,Vv;
int j;

int main()
{
        Tc = 408.2;  // Kelvin
        Pc = 36.5;  // bar
        Omega = 0.183; // acentric factor
        R = 83.14; //bar[cm3]/k/mol
	T = 300.1;
        Tre = T/Tc;
	RT = R*T;
	
	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==1)
                        {
                                pEnd = pp - deltaP;
                                pp = 59.0;
                        }
                }
                deltaP = deltaP/10;
        }
	cout<<"At Temp = "<<T<<" K"<<'\n';
	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);
	cout<<"Psat = "<<sol<<" bar"<<'\n';
	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';

	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);
}