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

Calculation of a refrigeration cycle using the Peng-Robinson equation

This notebook has been written in  Mathematica by 

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

mjm@nd.edu
http://www.nd.edu/~mjm/


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: 7/6/20

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 for calculation of a refrigeration cycle

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 a component

Critical constants for Methane

Molecular weight is 16

These are constants for Methane

componentname = Methane;
Pc=45.95/1.013;  
Tc=190.55;           
 ω=0.011;

Ideal gas heat capacity coefficients from Elliot and Lira

cpa=19.25;
cpb = .05213;
cpc=1.197 10^-5;
cpd=-1.132 10^-8;

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;

rr = 8.314;

If you clicked on the physical properties PReq is now

PReq

Enthalpy departure

Set up of departure functions

We would like to get the enthalpy departure functions from Peng-Robinson for various calculations.   Here is the definition from Sandler’s book.

HVdeptPR := R T (Z -1)+ Log[ ( Z+ ( 1 + Sqrt[2]) B)/( Z + ( 1 - Sqrt[2])B)]  /2/Sqrt[2]/b       ( - T .45723553 R^2 Tc^2/Pc  κ Sqrt[α/T/Tc] -a)

Some addition constants are needed

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

B:= b P /R/T;

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

Here is the entropy departure formula

SdeptPR := Log[ Z-B] R+ Log[ ( Z+ ( 1 + Sqrt[2]) B)/( Z + ( 1 - Sqrt[2])B)]  /2/Sqrt[2]/b( - .45723553 R^2 Tc^2/Pc  κ Sqrt[α/T/Tc])

idealgascp[T_, aa_,bb_,cc_,dd_] := aa + bb T + cc T^2 + dd T^3

hideal1 =Integrate[idealgascp[T, cpa,cpb,cpc,cpd],T]

sideal1 = Integrate[idealgascp[T, cpa,cpb,cpc,cpd]/T,T]

sidealfunction:=( sideal1/.T->T2 )- (sideal1/.T->T1)- 8.314 Log[p2/p1]

hidealfunction:= ( hideal1/.T->T2 )- (hideal1/.T->T1)

SdeptPR := Log[ Z-B] R+ Log[ ( Z+ ( 1 + Sqrt[2]) B)/( Z + ( 1 - Sqrt[2])B)]  /2/Sqrt[2]/b( - .45723553 R^2 Tc^2/Pc  κ Sqrt[α/T/Tc])

sreal :=(SdeptPR/.{T->T2, Z-> v2 p2/R/T2,P->p2})8.314/R + sidealfunction -  (SdeptPR/.{T->T1, Z-> v1 p1/R/T1,P->p1})8.314/R

HVdeptPR := R T (Z -1)+ Log[ ( Z+ ( 1 + Sqrt[2]) B)/( Z + ( 1 - Sqrt[2])B)]  /2/Sqrt[2]/b       ( - T .45723553 R^2 Tc^2/Pc  κ Sqrt[α/T/Tc] -a)

hreal :=(HVdeptPR/.{T->T2, Z-> v2 p2/R/T2,P->p2})8.314/R + hidealfunction -  (HVdeptPR/.{T->T1, Z-> v1 p1/R/T1,P->p1})8.314/R

Definition of fugacity coefficient to make the flash calculations easier

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

Start of refrigeration cycle calculation

For this calculation we pick the starting point to be the fluid as it leaves the cold side heat exchanger.  It will be chosen close to a saturated vapor.  Then we calculate the compression step.

Pick 120 K.

Tevapfinal = 120;

We should find the saturated conditions at the temperature that we want to have in the evaporators.  This will determine the pressure and serve as a starting point for the calculation

initialconditions = FindRoot[{  ((logfslashP/. {z-> P  vl/R/Tevapfinal})/.T->Tevapfinal) ==
((logfslashP/. {z-> P  vv/R/Tevapfinal})/.T->Tevapfinal) ,
  (P-PReq==0)/.{v->vl,T->Tevapfinal},  
  (P-PReq==0)/.{v->vv,T->Tevapfinal}},
{{P,10},{vl,.03},{vv,1}}]

Pevapfinal = P/.initialconditions

Check the volume and find compressibility at these conditions

satvap = Solve[Evaluate[Pevapfinal-PReq/.{T->Tevapfinal}] == 0,v]

zees = Pevapfinal * v/R/Tevapfinal;

zees/.satvap

Calculate the compressor final temperature for chosen final pressure

Choose a pressure for the exit of the compressor.  We don’t have to worry about the phase changing!

Pcompressorfinal = 10;

Clear[T1]

compressans = FindRoot[{(sreal==0)/.{p1->Pevapfinal,p2->Pcompressorfinal,T1->Tevapfinal},(Pcompressorfinal-PReq==0)/.{T->T2,v->v2},(Pevapfinal-PReq==0)/.{T->Tevapfinal,v->v1}},{{T2,300},{v2,3},{v1,5}}]

Tcompressexit = T2/.compressans

Vcompressexit = v2/.compressans

Here is just a check of the answer.  If we know the final temperature, we can get the volume.

Solve[Evaluate[Pcompressorfinal-PReq/.{T->Tcompressexit}] == 0,v]

The enthalpy change is

compressorwork =(hreal/.{p1->Pevapfinal,p2->Pcompressorfinal,T1->Tevapfinal})/.compressans

Vcompressexit

Tcompressexit

calculate the hotside heat exchanger.  (constant pressure heat removal.)

Check to see the temperature at saturation for the chosen pressure

heatxconditions = FindRoot[{  (logfslashP/. {P->Pcompressorfinal,z-> Pcompressorfinal vl/R/T}) ==
(logfslashP/. {P->Pcompressorfinal,z-> Pcompressorfinal vv/R/T}) ,
  (Pcompressorfinal-PReq==0)/.{v->vl},  
  (Pcompressorfinal-PReq==0)/.{v->vv}},
{{T,Tcompressexit-50 },{vl,Vcompressexit/30},{vv,Vcompressexit .7}}]

Tsatheatexchange = T/.heatxconditions

this is just a constant P calculation,  Find v at chosen T and then get h

Tcondensor =149.3

Thotsidefinal = 149.3

satliqhot = Solve[Evaluate[Pcompressorfinal-PReq/.{T->Thotsidefinal}] == 0,v]

vhotexit = v/.satliqhot[[1]]

The condenser heat load is:

hreal/.{p1->Pcompressorfinal,p2->Pcompressorfinal,T1->Tcompressexit,T2->Thotsidefinal,v1->Vcompressexit,v2->vhotexit}

8400 J/mol (mol/16 /g)

%/.{g-> kg/1000, J->kJ/1000.}

Isenthalpic pressure drop   (Joule-Thompson expansion)

Pick the final pressure for the expansion

Pevapfinal

We need to do a flash calculation that is adiabatic.  
β is the fraction of vapor.

evapans = FindRoot[{  (logfslashP/. {P->Pevapfinal,z-> Pevapfinal vl/R/T}) ==
(logfslashP/. {P->Pevapfinal,z-> Pevapfinal vv/R/T}) ,
  (Pevapfinal-PReq==0)/.{v->vl},  
  (Pevapfinal-PReq==0)/.{v->vv},
  ((8.314/R  HVdeptPR/.Z->P v/R/T)/.{v->vhotexit,T->Thotsidefinal,P->Pcompressorfinal} )-(hidealfunction/.{T1->Thotsidefinal,T2->T}) ==
  (         (1-β)((  8.314/R HVdeptPR/.Z->P vl/R/T)/.{v->vl,P->Pevapfinal} )+ (β ((8.314/R HVdeptPR/.Z->P vv/R/T)/.{v->vv,P->Pevapfinal})))},
{{T,140},{vl,.04},{vv,5},{β,.8}}]

Tevap = T/.evapans

vevapliquid = vl/.evapans

vevapvapor = vv/.evapans

ϕvap = β/.evapans

This is the delta h along an isotherm from the liquid to the vapor volumes.  Hence... the heat of vaporization.  The (1- ϕ) gives the amount of liquid that can be evaporated in the “evaporator” inside your freezer.

evapheat =(1-ϕvap) hreal/.{p1->Pevapfinal,p2->Pevapfinal,T2->Tevapfinal,T1->Tevapfinal,v2->vevapvapor,v1->vevapliquid}

The coefficient of performance the heat removed in the evaporator divided by the work input to the compressor.

evapheat/compressorwork

To match the p-H diagram for methane, you need to turn these Joules/(gram- mole) into kJ/kg.