• Procedures to calculate the parameters, the fugacities and solve the EOS.

Calculation of the Bubble Point Pressure for a mixture using the Peng Robinson Equation of State and simple mixing rules.

2000 John Kitchin
All Rights Reserved
University of Delaware
Department of Chemical Engineering

> restart:

This worksheet calculates the bubble point pressure of a mixture using the complete generalized Peng-Robinson equation of state shown below. The parameters are calculated using the Van der Waal mixing rules from the pure component properties.

The bubble point pressure of a mixture is defined as the pressure where the first bubble of vapor is formed at a given temperature. The strategy to calculate this pressure is to assume that the very first bubble will not change the composition of the liquid, and will not have the same composition as the liquid (otherwise it would be the liquid phase). We will make a guess for what the pressure and vapor composition should be, and then iterate on both of them until all of the species fugacities (calculated using the Peng Robinson EOS for a mixture) are equal in each phase. The end result will be the bubble point pressure and vapor composition in equilibrium with the liquid mixture.

The algorithm is based on the one presented in Chemical and Engineering Thermodynamics, 3rd Ed., S.I. Sandler on page523. All of the data used is obtained from tables in this book.

> R:=0.00008314: gas constant
Property matrix of critical properties of pure components in form of [[Tc1, Pc1, omega1], [Tc2, Pc2, omega2]] Tc in K, Pc in bar. In this case, CO2 and n-butane.
PropMatrix:=[[304.2,73.76,.225],[425.2,38.00,.193]];
k:=[[0,0.130],[0.130, 0]]; binary interaction parameter if known (set all elements to 0 if not known).

Procedures to calculate the parameters, the fugacities and solve the EOS.

These are the equations used to calculate the coefficients of the PR EOS for the pure components.

> PRparamsPure:=proc()
local q; global Tc,Pc,omega,kap,b,ac,a;
for q from 1 to 2 do
Tc[q]:=PropMatrix[q,1]: assign the critical properties of the pure components.
Pc[q]:=PropMatrix[q,2]:
omega[q]:=PropMatrix[q,3];
kap[q]:=0.37464+1.54226*omega[q]-0.26992*(omega[q])^2;
b[q]:=0.07780*R*Tc[q]/Pc[q]:
ac[q]:=0.45724*R^2*Tc[q]^2/Pc[q];
a[q] := ac[q]*(1+kap[q]*(1-sqrt(T/Tc[q])))^2
od:
end:

>

These are the mixing rules used to calculate the parameters for the mixture and the EOS of the mixture in terms of the compressibility. They are known as the Van der Waal one-fluid mixing rues. The equations for the coefficients are

where refers to the pure component parameter, and refers to the binary interaction parameter.

For simplicity, I created two procedures for calculating the parameters for the EOS, one for the vapor phase which uses the vapor composition, and one for the liquid that uses the liquid composition.

> PRmixParamsL:=proc()
global am,bm,A,B,EOS;
PRparamsPure();
am:=sum(sum('x[n]*x[m]*sqrt(a[m]*a[n])*(1-k[m,n])','m' = 1 .. 2),'n' = 1 .. 2);
bm:=sum('x[i]*b[i]','i'=1..2);
A:=am*P/(R*T)^2;
B:=P*bm/R/T;
EOS:=Z^3+(B-1)*Z^2+(A-3*B^2-2*B)*Z+(-A*B+B^2+B^3)=0: the EOS in terms of the compressibility.
end:

> PRmixParamsV:=proc()
global am,bm,A,B,EOS;
PRparamsPure();
am:=sum(sum('y[n]*y[m]*sqrt(a[m]*a[n])*(1-k[m,n])','m' = 1 .. 2),'n' = 1 .. 2);
bm:=sum('y[i]*b[i]','i'=1..2);
A:=am*P/(R*T)^2;
B:=P*bm/R/T;
EOS:=Z^3+(B-1)*Z^2+(A-3*B^2-2*B)*Z+(-A*B+B^2+B^3)=0:
end:

This procedure solves the EOS for the compressibility roots, sets all imaginary roots to zero, sorts them and takes the largest for ZV and smallest for ZL.

> PRsolvEOS:=proc()
global ZZ,ZL,ZV,Z1,Z2,Z3,zz;
ZZ:=solve(EOS,Z);
Z1:=ZZ[1];Z2:=ZZ[2];Z3:=ZZ[3];
if abs(Im(Z1))>0 then Z1:=0 fi;
if abs(Im(Z2))>0 then Z2:=0 fi;
if abs(Im(Z3))>0 then Z3:=0 fi;
zz:=sort([Z1,Z2,Z3]);
if zz[1]=0 then zz[1]:=zz[3] fi; These two steps are necessary to ensure two real roots, one for vapor and one for liquid.
if zz[3]=0 then zz[3]:=zz[1] fi;
ZL:=zz[1]; set the smallest root to the liquid compressibility
ZV:=zz[3]; set the largest root to the vapor compressibility
end:

Procedure to calculate fugacity in the liquid phase using the liquid composition and compressibility.

> PRfugl:=proc()
local i;
global fl,phil,fugl;
for i from 1 to 2 do
fl[i]:=(ZL-1)*b[i]/bm-ln(ZL-B)-am/(2*sqrt(2)*bm*R*T)*(2*sum(x[j]*sqrt(a[i]*a[j])*(1-k[i,j]),j=1..2)/am-b[i]/bm)*ln((ZL+(1+sqrt(2))*B)/(ZL+(1-sqrt(2))*B)):
phil[i]:=exp(fl[i]):
fugl[i]:=evalf(P*x[i]*phil[i]);
od:
end:

Procedure to calculate fugacity of species in the vapor phase using the vapor composition and vapor compressibility.

> PRfugv:=proc()
local i;
global fv,phiv,fugv;
for i from 1 to 2 do
fv[i]:=(ZV-1)*b[i]/bm-ln(ZV-B)-am/(2*sqrt(2)*bm*R*T)*(2*sum(y[j]*sqrt(a[i]*a[j])*(1-k[i,j]),j=1..2)/am-b[i]/bm)*ln((ZV+(1+sqrt(2))*B)/(ZV+(1-sqrt(2))*B)):
phiv[i]:=exp(fv[i]):
fugv[i]:=evalf(P*y[i]*phiv[i]);
od:
end:

This procedure iterates on the vapor phase composition until the vapor phase species fugacity is equal to the liquid phase species fugacity.

> InnerIterate:=proc()
global fugv,y; local i,test1,test2;
for i from 1 to 50 do
PRmixParamsV(); Calculate the parameters for the vapor phase EOS
PRsolvEOS(); Solve the EOS
PRfugv(); Calculate the vapor phase fugacities.
These two if/then statements check to see if the species fugacities are the same within the specified tolerance and adjust the mole fractions accordingly.
if (abs(1-fugl[1]/fugv[1])>.00001) then y[1]:=y[1]*fugl[1]/fugv[1]; else test1:=1; fi;
if (abs(1-fugl[2]/fugv[2])>.00001) then y[2]:=y[2]*fugl[2]/fugv[2]; else test2:=1 fi;
if test1=1 and test2=1 then break fi; If both mole fractions have converged then break the iteration loop.
od:
end:

Once the vapor phase fugacity and liquid phase fugacity are equal, we check to make sure the vapor phase mole fractions sum to 1, if not, we iterate on the pressure guess until they do. When the iteration is done, the final pressure and vapor compositions are the Bubble point pressure and the vapor composition.

> BubblePressure:=proc()
local test;
global P,i;
for i from 1 to 50 do
PRmixParamsL(); The structure of the procedures is to calculate the parameters for the EOS for the conditions desired,
PRsolvEOS(); then use a separate procedure to solve it.
PRfugl(); and finally, a procedure to calculate the species fugacities.
InnerIterate();
if ((abs((y[1]+y[2])-1)>0.00001)) then P:=P*(y[1]+y[2]) else test:=1 fi;
if test=1 then break fi;
od:
end:

> T:=273.15+50: Temperature of the feed in Kelvin
x[1]:=.3585: mole fraction of componenet 1 in the feed
x[2]:=1-x[1]:
P:=40; estimated bubble point pressure, perhaps where Pvap could be calculated from the pure component properties and the vapor pressure worksheet. This value will be iterated on.
y[1]:=.8; estimated mole fraction of component 1 in vapor phase perhaps . If you don't know, guess a high number like 0.99. It seems to work better if the guess is above the actual value. These values will also be iterated on. Note that y[i] refers to the mole fraction in the vapor phase and x[i] is the mole fraction in the liquid phase.
y[2]:=1-y[1];

> BubblePressure(); Call the procedure that does the iterations.

> "Bubble Point Pressure (bar): ",P;"y[1] ",y[1];"y[2] ",y[2];

If your initial guesses are not very good, the program may not converge, or can converge to a trivial solution (that is where the vapor phase has the same composition as the liquid phase and so there is only one phase).

In conclusion, we have successfully applied Maple to a calculation that would be extremely tedious to carry out by hand. Furthermore, by creating procedures for each step of the calculation, it is very easy to perform a new calculation for different conditions or different mixtures. We need only change the pure component properties and the desired conditions and rerun a single command, BubblePressure(). This makes the worksheet very amenable to calculating the bubble pressure of a mixture across a broad range of temperatures and compositions with very simple do loops.