PengRobinson EOS in FPROPS

Documentation

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.

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