Properties of liquids, gases and vapors from the Peng-Robinson equation

Calculation of volumetric behavior, two-phase region, departure functions, isentropic path, fugacity and the two-phase region by equal fugacities.

This notebook has been written in  Mathematica by 

Mark J. McCready
Professor of Chemical and Biomolecular Engineering
University of Notre Dame
Notre Dame, IN 46556
USA

mjm@nd.edu
http://www.chemeprof.com


It is copyrighted to the extent allowed by which ever laws pertain to the World Wide Web and the Internet.

I would hope that as a professional courtesy, this notice remain visible to other users. 
There is no charge for copying and dissemination   

Version: 10/25/18

Homework

#1

Background

To design and operate chemical processes, the pressure, temperature and volumetric information of  pure substances and mixtures is often required.  In addition, if there is any energy exchange occurring, properties such as enthalpy and entropy will also be needed.

A common calculation procedure for such information is to base the T,P, V behavior on a cubic equation of state and then to get the energetic information from “departure functions”.  This Notebook illustrates this using the Peng-Robinson equation

Definition of the terms of the P-R equation

When you use this notebook, start here.  Always hit these initializations if you have done calculations in Mathematica using other equations and definition to re-set everything for the calculations below.

Note that the “:=”, tells Mathematica to recalculate “b” using the most recent values of “Tc and Pc”.  So if half way down, you want to change to a different substance, then it will calculate new “a’s. b’s”. ....

Set up of equation and constants

The constants

b:=.07779607 R Tc/Pc;

a:= .45723553 R^2 Tc^2/Pc (1 + κ (1-Sqrt[T/Tc]))^2;

κ:= .37464 + 1.54226 ω - .26992 ω^2

The equation

PReq := R T/( v-b) - a / (v ( v+b) + b ( v-b))

We “call its name” to see the constants substituted in.

PReq

Pick the components

Critical constants for Propane

componentname = hexane;
Pc=30.1/1.013;  
Tc=507.5;           
 ω=0.299;

Critical constants for Propane

componentname = octane;
Pc=24.9/1.013;  
Tc=568.8;           
 ω=0.398;

Gas constant

The volumetric gas constant is used with conversion at the end to Joules for Enthalpy and Entropy.  It is easier to give it a value right here.

R =.082057;

If you clicked on the physical properties PReq is now

PReq

Check of updating of “:=” definitions

And for a check I’ll change Tc and change it back

Tc=500

PReq

You can see that PReq updated.

Tc = 369.9;

PReq

Pressure volume plots using the P-R equation

Let’s look at the behavior of this equation

Above the critical temperature, this is typical behavior

Plot[PReq/.T->530,{v,.12,3}, PlotRange->{-20,100},AxesLabel->{"volume, (l)","pressure (ATM)"}, PlotLabel->"Peng-Robinson Equation for" componentname]

Plot[PReq/.T->490,{v,.13,4}, PlotRange->{-100,100},AxesLabel->{"volume, (l)","pressure (ATM)"}, PlotLabel->"Peng-Robinson Equation for" componentname]

Plot[PReq/.T->440,{v,.15,3}, PlotRange->All,AxesLabel->{"volume, (l)","pressure (ATM)"}, PlotLabel->"Peng-Robinson Equation for" componentname]

Plot[PReq/.T->400,{v,.15,3}, PlotRange->All,AxesLabel->{"volume, (l)","pressure (ATM)"}, PlotLabel->"Peng-Robinson Equation for" componentname]

Close to the critical point you see a flat spot developing

Plot[PReq/.T->370,{v,.13,2},PlotRange->All, AxesLabel->{"volume, (l)","pressure (ATM)"}, PlotLabel->"Peng-Robinson Equation for" componentname]

Below the critical point, there is a definite “dome”

Plot[PReq/.T->350,{v,.13,1},PlotRange->All, AxesLabel->{"volume, (l)","pressure (ATM)"}, PlotLabel->"Peng-Robinson Equation for" componentname]

Since we can’t have negative volume, there must be more to this problem. 

What we realize is that below the critical temperature, there is a two-phase region!

Here is the overall behavior.

Show[%,%%,%%%,%%%%,%%%%%]

From the plots we see a couple of problems.

1.  Below the critical temperature, there are 3 volumes corresponding to a given pressure.
2.  In the two phase region, for a single component, F = C- P +2 = 1, so that if T is chosen, which value P should we pick.

From class we have learned that the answer to the first question is that the equation is giving a liquid volume (smallest), a vapor volume (largest) and an intermediate volume that corresponds to an unstable state.  So for the two-region calculations, we choose the smallest and largest roots as needed.

The second question is a bit more subtle.  The value of P for a given T in the two-phase region is determined by the smallest value of the Gibbs Energy.  If this criteria is applied then if a horizontal line is drawn at the correct pressure, the area of the P-V curve above this line will be the same as the area below this line.  Thus we can applied an integrated area criteria and find the pressure by iteration.

Numerical behavior of equation

Root behavior

Let’s make sure we know how the roots come out.

For T=380 > Tc here are the values of roots. p=10ATM.  The p- PReq ==0 is the Peng-Robinson equation and I pick T and P to substitute into it.

Solve[(p-PReq==0)/.{T->520, p->10},v]

The result is a “gas” with a volume V = 3.02 l and two complex conjugate roots that don’t correspond to physically occurring volumes.

Now look at lower T.  T = 320, P = 20

Solve[(p-PReq==0)/.{T->440, p->10},v]

The result is a liquid with a volume V = .0949 and a vapor with a volume, V = .825 l.

Correct pressure for a given isotherm

Here is the correct pressure, P= 24.11 ATM,  for T = 340K.  It is hard to “eyeball” these.

Plot[{PReq/.T->460,14.11},{v,.13,3},PlotRange->{-20,60},AxesLabel->{"volume, (l)","pressure (ATM)"}, PlotLabel->"Peng-Robinson Equation for" componentname]

Below we will write an iteration loop, but we can try and adjust to find a value for P (if say the integration does not converge or if we want to make sure how to do it.
First get the three roots,
pick P and T.

ptest = 15;
Ttest = 460;

rootscheck= Solve[(p-PReq==0)/.{T->Ttest, p->ptest},v]

Now we can get the integral.  We know from below that for an accurate answer, the integral will need to be less that .01 or so.

ans1 =NIntegrate[(p-PReq)/.{T->Ttest,p->ptest},{v,v/.rootscheck[[1]],v/.rootscheck[[3]]}]

Try another pressure.

ptest = 14;

rootscheck= Solve[(p-PReq==0)/.{T->Ttest, p->ptest},v]

Now we can get the integral.  We know from below that for an accurate answer, the integral will need to be less that .01 or so.

ans1 =NIntegrate[(p-PReq)/.{T->Ttest,p->ptest},{v,v/.rootscheck[[1]],v/.rootscheck[[3]]}]

This is much better.  If we try p = 24.115

ptest = 14.425;

rootscheck= Solve[(p-PReq==0)/.{T->Ttest, p->ptest},v]

Now we can get the integral.  We know from below that for an accurate answer, the integral will need to be less that .01 or so.

ans1 =NIntegrate[(p-PReq)/.{T->Ttest,p->ptest},{v,v/.rootscheck[[1]],v/.rootscheck[[3]]}]

We have a good result!

So if you need to do this “by hand” use the equations above.

Two phase region calculation

Here I initialize and run a double loop that computes the P-T values for propane from T = 315 up to close to the critical point.  If you want higher and lower values maybe you can adjust the initial pressure test or “delta-P” numbers or something else.  This was all I could get for propane in an automatic loop.

Note that I just increment T.  The pressure is adjusted in the inner loop with a discrete Newton-Raphson correction procedure (which in case you are interested is just the first term of a Taylor series expansion of a function!).

T0 = 465;
deltp=.1;
p1test = 15;
ptest = 15;

For[i =0,i<7,i++;tt[i] = T0+ (i-1)*5;
p1test =Re[ ptest];
xtest =1;
xtestold=2;
roots1= Solve[p1test-(PReq/.T->tt[i]) ==0,v];
ans1 =NIntegrate[p1test-(PReq/.T->tt[i]),{v,v/.roots1[[1]],v/.roots1[[3]]}];
p2test = p1test + deltp;
roots =Solve[p2test-(PReq/.T->tt[i]) ==0,v];
ans2=NIntegrate[p2test-(PReq/.T->tt[i]),{v,v/.roots[[1]],v/.roots[[3]]}];
ptest=p1test-ans2(ans2-ans1)/deltp;
ansold = ans2;
poldtest = p2test;
While[Abs[xtest]-(1/(10^2))≥0,
roots =Solve[ptest-(PReq/.T->tt[i]) ==0,v];
ans1= NIntegrate[p1test-(PReq/.T->tt[i]),{v,v/.roots[[1]],v/.roots[[3]]}];
xtest = ans1;
ptest = Re[p1test -xtest(ans1-ansold)/(p1test -poldtest)];
ansold = ans1; poldtest = p1test;p1test = ptest;
Print["roots =  ",roots,"  xtest=",xtest,"  p = ",ptest]];
xtest = 1;
Print["p=  ",ptest,"  T= ",tt[i], "  i = ",i]]

Fugacity calculation and use of this to determine the two-phase region

We have established that requiring the Gibbs energy for a component for each phase must be equal at equilibrium.  This is the same as saying that the fugacity of a component is the same in each phase.  It turns out that the fugacity of a component can be easily calculated from the equation of state using a “fugacity coefficient” ϕ = f/P.

		(f^L)=(f^V)
		(ϕ^L)P = (ϕ^V)P
		
but of course the pressure is the same in each phase so that phase equilibrium occurs, at const T and P, if the fugacity coefficients are equal.

		(ϕ^L) = (ϕ^V)
		
The fugacity coefficient for the Peng - Robinson equation is conveniently given as ln[ϕ] =((A ln(((((Sqrt[2])+1) B+z)/((1-(Sqrt[2])) B+z))))/((2 (Sqrt[2])) B))-ln(z-B)+(z-1)  and the value for the liquid or vapor is obtained by merely substituting the correct value of z for either the liquid and vapor.

The natural log of the fugacity coefficient for the Peng Robinson equation is given by:

Some addition constants are needed

A:=a P /R^2 /T^2;

B:= b P /R/T;

α := ( 1 + κ ( 1- Sqrt[T/Tc]))^2

logfslashP := (z-1) - Log[(z-B)] - A/2/Sqrt[2]/B Log[( z + (1 + Sqrt[2])B)/( z + (1 - Sqrt[2]) B)]

logfslashP

In this equation, T, P and z are explicit, but recall that in the two-phase region, which is where we would want to calculate equilibrium for  two-phase mixture, picking T specifies P.  From the P-R equation, you can then calculate liquid and vapor z’s.

T1=440;

p1-PReq/.T->T1

Check of equation and premise from our two-phase region calculation above.

Solve[(p1-PReq==0)/.T->T1]

T1 = 460, p1= 16 ATM, does this work?

T1=460;
p1=16.2;

FindRoot[(p1-PReq ==0)/.{T->T1,P->p1},{v,1}]

FindRoot[(p1-PReq ==0)/.{T->T1,P->p1},{v,.1}]

This might have given us the correct roots, but it is easier to just solve the PR equation so we don’t make a mistake.

Solve[(p1-PReq ==0)/.{T->T1,P->p1},v]

Now is this an equilibrium state?  Will the fugacities for the liquid and vapor match?

logfslashP/.{P->p1,T->T1,z-> p1 1.45/R/T1}

logfslashP/.{P->p1,T->T1,z-> p1 0.1976/R/T1}

So these do not match, at T= 315, p = 10 ATM is not the correct pressure.  

Let’s guess p = 14.12  !

p1=14.12;

FindRoot[(p1-PReq ==0)/.{T->T1,P->p1},{v,1}]

FindRoot[(p1-PReq ==0)/.{T->T1,P->p1},{v,.1}]

This might have given us the correct roots, but it is easier to just solve the PR equation so we don’t make a mistake.

Solve[(p1-PReq ==0)/.{T->T1,P->p1},v]

Check the fugacity coefficients.

logfslashP/.{P->p1,T->T1,z-> p1 1.388224/R/T1}

logfslashP/.{P->p1,T->T1,z-> p1 0.093064425/R/T1}

Now this does work as expected from above.

Calculation of the two-phase region (hexane)

We can re-do the two phase region calculation that was done above, but now using fugacity equality of the two phases to pick P, if T is incremented.

The calculation is easily done from 430 up to the critical temperature.  You could probably go lower by adjusting the “guess” for “vv” the vapor volume. (now at 2 l)

Clear[p,vvl,vvv]

T0=430;

For[i =0,i<17,i++;tt[i] = T0+ (i-1)*5;
T1=tt[i];
ans = FindRoot[{(logfslashP/. {P->p,z-> p vl/R/T1,T->T1}) ==(logfslashP/. {P->p,z-> p vv/R/T1,T->T1}) ,(p-PReq==0)/.{v->vl,T->T1},(p-PReq==0)/.{T->T1,v->vv}},{{p,10},{vl,.2},{vv,1.3}}];
pp[i] =pstep = p/. ans;
vvl[i]=vliquid = vl/. ans;
vvv[i] = vvapor = vv/. ans;
zvapor = pstep vvapor/R/tt[i];
fugacitycoeff = Exp[logfslashP/. {P->pstep,z->  pstep vvapor/R/tt[i],T->tt[i],vv->vvapor}];
Print["temperature =  ", tt[i],"  p = ",pstep " ,  vliquid =  ",vliquid,"  vvapor=  ",vvapor,"  zvapor=  ",zvapor,"  ϕ = ",fugacitycoeff];
]

satliqdata =Table[{vvl[i],pp[i]},{i,1,17}]

satvapdata =Table[{vvv[i],pp[i]},{i,1,17}]

ListPlot[{satliqdata,satvapdata},AxesLabel->{"volume","pressure"},PlotLabel->"Hexane"]

vaporpressuredata =Table[{tt[i],pp[i]},{i,1,17}]

vaporpressuredata ={{430,8.463178886822217`},{435,9.2969092403946`},{440,10.190945359163154`},{445,11.148098477292745`},{450,12.171251702520605`},{455,13.26336348533547`},{460,14.427471526058188`},{465,15.666697147213075`},{470,16.98425016255361`},{475,18.38343427866704`},{480,19.867653070298992`},{485,21.440416576509637`},{490,23.105348571606687`},{495,24.86619457263638`},{500,26.72683065421613`},{505,28.6912731518667`},{510,30.703898897890003`}}

ListPlot[{vaporpressuredata},AxesLabel->{"temperature","pressure"},PlotLabel->"Hexane"]

Plot[10^(4 - (1171.53/(T-48.78))),{T,430,507}, AxesLabel->{"T, K","P, ATM"},PlotLabel->"Antoine Vapor pressure prediction"]

Show[%,%%]

Now let’s do Octane

Critical constants for Octane

componentname = octane;
Pc=24.9/1.013;  
Tc=568.8;           
 ω=0.398;

Gas constant

The volumetric gas constant is used with conversion at the end to Joules for Enthalpy and Entropy.  It is easier to give it a value right here.

R =.082057;

Fugacity calculation and use of this to determine the two-phase region

We have established that requiring the Gibbs energy for a component for each phase must be equal at equilibrium.  This is the same as saying that the fugacity of a component is the same in each phase.  It turns out that the fugacity of a component can be easily calculated from the equation of state using a “fugacity coefficient” ϕ = f/P.

		(f^L)=(f^V)
		(ϕ^L)P = (ϕ^V)P
		
but of course the pressure is the same in each phase so that phase equilibrium occurs, at const T and P, if the fugacity coefficients are equal.

		(ϕ^L) = (ϕ^V)
		
The fugacity coefficient for the Peng - Robinson equation is conveniently given as ln[ϕ] =((A ln(((((Sqrt[2])+1) B+z)/((1-(Sqrt[2])) B+z))))/((2 (Sqrt[2])) B))-ln(z-B)+(z-1)  and the value for the liquid or vapor is obtained by merely substituting the correct value of z for either the liquid and vapor.

The natural log of the fugacity coefficient for the Peng Robinson equation is given by:

Some addition constants are needed

A:=a P /R^2 /T^2;

B:= b P /R/T;

α := ( 1 + κ ( 1- Sqrt[T/Tc]))^2

logfslashP := (z-1) - Log[(z-B)] - A/2/Sqrt[2]/B Log[( z + (1 + Sqrt[2])B)/( z + (1 - Sqrt[2]) B)]

logfslashP

In this equation, T, P and z are explicit, but recall that in the two-phase region, which is where we would want to calculate equilibrium for  two-phase mixture, picking T specifies P.  From the P-R equation, you can then calculate liquid and vapor z’s.

T1=490;
p1=10;

p1-PReq/.T->T1

Check of equation and premise from our two-phase region calculation above.

Solve[(p1-PReq==0)/.T->T1]

T1 = 460, p1= 16 ATM, does this work?

T1=460;
p1=6;

FindRoot[(p1-PReq ==0)/.{T->T1,P->p1},{v,1}]

FindRoot[(p1-PReq ==0)/.{T->T1,P->p1},{v,.1}]

This might have given us the correct roots, but it is easier to just solve the PR equation so we don’t make a mistake.

Solve[(p1-PReq ==0)/.{T->T1,P->p1},v]

Now is this an equilibrium state?  Will the fugacities for the liquid and vapor match?

logfslashP/.{P->p1,T->T1,z-> p1 1.45/R/T1}

logfslashP/.{P->p1,T->T1,z-> p1 0.1976/R/T1}

So these do not match, at T= 315, p = 10 ATM is not the correct pressure.  

Let’s guess p = 14.12  !

p1=4;

FindRoot[(p1-PReq ==0)/.{T->T1,P->p1},{v,1}]

FindRoot[(p1-PReq ==0)/.{T->T1,P->p1},{v,.1}]

This might have given us the correct roots, but it is easier to just solve the PR equation so we don’t make a mistake.

Solve[(p1-PReq ==0)/.{T->T1,P->p1},v]

Check the fugacity coefficients.

logfslashP/.{P->p1,T->T1,z-> p1 1.388224/R/T1}

logfslashP/.{P->p1,T->T1,z-> p1 0.093064425/R/T1}

Now this does work as expected from above.

Calculation of the two-phase region (hexane)

We can re-do the two phase region calculation that was done above, but now using fugacity equality of the two phases to pick P, if T is incremented.

The calculation is easily done from 430 up to the critical temperature.  You could probably go lower by adjusting the “guess” for “vv” the vapor volume. (now at 2 l)

Clear[p,vvl,vvv]

T0=470;

For[i =0,i<22,i++;tt[i] = T0+ (i-1)*5;
T1=tt[i];
ans = FindRoot[{(logfslashP/. {P->p,z-> p vl/R/T1,T->T1}) ==(logfslashP/. {P->p,z-> p vv/R/T1,T->T1}) ,(p-PReq==0)/.{v->vl,T->T1},(p-PReq==0)/.{T->T1,v->vv}},{{p,4.3},{vl,.2},{vv,1.3}}];
pp[i] =pstep = p/. ans;
vvl[i]=vliquid = vl/. ans;
vvv[i] = vvapor = vv/. ans;
zvapor = pstep vvapor/R/tt[i];
fugacitycoeff = Exp[logfslashP/. {P->pstep,z->  pstep vvapor/R/tt[i],T->tt[i],vv->vvapor}];
Print["temperature =  ", tt[i],"  p = ",pstep " ,  vliquid =  ",vliquid,"  vvapor=  ",vvapor,"  zvapor=  ",zvapor,"  ϕ = ",fugacitycoeff];
]

satliqdata =Table[{vvl[i],pp[i]},{i,1,22}]

satvapdata =Table[{vvv[i],pp[i]},{i,1,22}]

ListPlot[{satliqdata,satvapdata},AxesLabel->{"volume","pressure"},PlotLabel->"Octane"]

vaporpressuredata =Table[{tt[i],pp[i]},{i,1,22}]

ListPlot[{vaporpressuredata},AxesLabel->{"temperature","pressure"},PlotLabel->"Hexane"]

Plot[10^(4 - (1171.53/(T-48.78))),{T,430,507}]

Plot[{10^(4 - (1171.53/(T-48.78))),10^(4.0487 - (1355.126/(T-63.633)))},{T,350,420},PlotRange->{0,3}]

Show[%,%107]

Now do a mixture.

Critical constants for Propane

componentname = hexane;
Pc=30.1/1.013;  
Tc=507.5;           
 ω=0.299;

Critical constants for Octane

componentname = octane;
Pc=24.9/1.013;  
Tc=568.8;           
 ω=0.398;

The constants

b:=.07779607 R Tc/Pc;

a:= .45723553 R^2 Tc^2/Pc (1 + κ (1-Sqrt[T/Tc]))^2;

κ:= .37464 + 1.54226 ω - .26992 ω^2

A:=a P /R^2 /T^2;

B:= b P /R/T;

α := ( 1 + κ ( 1- Sqrt[T/Tc]))^2

logfslashP := (z-1) - Log[(z-B)] - A/2/Sqrt[2]/B Log[( z + (1 + Sqrt[2])B)/( z + (1 - Sqrt[2]) B)]