# Variables
P = 2; #number of phases (no unit)
C = 2; #number of components (no unit)
# Calculations
F = C+2-P
# Results
print " The number of degrees of freedom = %d "%(F);
print "Two intensive properties are required to be specified to\
describe the thermodynamic state of the system,and the fundamental relation in\
the Gibbs free energy representation for this system is of the type, G = GT,P,N1,N2"
import cmath
import math
# Variables
T = 427.85; #temperature of n-octane vapour in K
R = 8.314; #universal gas constant in J/molK
Tc = 569.4; #critical temperature of n-octane in K
Pc = 24.97; #critical pressure of n-octane in bar
w = 0.398; #acentric factor (no unit)
# Calculations
Pguess = 0.215
Tr = T/Tc
Pr = (Pguess*10**6)/(Pc*10**5)
S = 0.37464+(1.54226*w)-(0.26992*w**2)
alpha1 = (1+(S*(1-math.sqrt(Tr))))**2;
a = (0.45724*R**2*Tc**2*alpha1)/(Pc*10**5)
b = (0.07780*R*Tc)/(Pc*10**5);
A = (a*Pguess*10**6)/(R*T)**2;
B = (b*Pguess*10**6)/(R*T);
alpha = -1+B
beeta = A-(2*B)-(3*B**2)
gaamma = -(A*B)+(B**2)+(B**3)
p = beeta-(alpha**2)/3;
q = ((2*alpha**3)/27)-((alpha*beeta)/3)+gaamma
D = (((q)**2)/4)+(((p)**3)/27);
if D>0:
Z = ((-q/2)+(math.sqrt(D)))**(1./3)+((-q/2)-(math.sqrt(D)))**(1./3)-(alpha/3);
Z_l = Z;
Z_v = Z;
elif D == 0:
Z1 = ((-2*(q/2))**(1./3))-(alpha/3);
Z2 = ((q/2)**(1./3))-(alpha/3);
Z3 = ((q/2)**(1./3))-(alpha/3);
Z = [Z1 ,Z2, Z3];
Z_l = min(Z);
Z_v = max(Z);
else:
r = math.sqrt((-(p**3)/27));
theta = math.acos((-(q)/2)*(1./r));
Z1 = (2*(r**(1./3))*math.cos(theta/3))-(alpha/3);
Z2 = (2*(r**(1./3))*math.cos(((2*math.pi)+theta)/3))-(alpha/3)
Z3 = (2*(r**(1./3))*math.cos(((4*math.pi)+theta)/3))-(alpha/3)
Z = [Z1, Z2, Z3];
Z_l = min(Z)
Z_v = max(Z)
phi_l = math.exp(Z_l-1-math.log(Z_l-B)-((a/(2*math.sqrt(2)*b*R*T))*math.log((Z_l+(B*(1+math.sqrt(2))))/(Z_l+(B*(1-math.sqrt(2)))))));
phi_v = math.exp(Z_v-1-math.log(Z_v-B)-((a/(2*math.sqrt(2)*b*R*T))*math.log((Z_v+(B*(1+math.sqrt(2))))/(Z_v+(B*(1-math.sqrt(2)))))));
fl = Pguess*phi_l; # Calculations of the fugacity of the liquid in MPa
fv = Pguess*phi_v; # Calculations of the fugacity of the vapour in MPa
tolerance = 1e-3; #defining the tolerance to compare fl and fv
if abs(fl-fv)<tolerance:
P = Pguess;
else:
Prevised = Pguess*(fl/fv)
while abs(fl-fv)>tolerance:
Tr = T/Tc
Pr = (Prevised*10**6)/(Pc*10**5);
S = 0.37464+(1.54226*w)-(0.26992*w**2)
alpha1 = (1+(S*(1-math.sqrt(Tr))))**2;
a = (0.45724*R**2*Tc**2*alpha1)/(Pc*10**5)
b = (0.07780*R*Tc)/(Pc*10**5);
A = (a*Prevised*10**6)/(R*T)**2;
B = (b*Prevised*10**6)/(R*T);
alpha = -1+B;
beeta = A-(2*B)-(3*B**2);
gaamma = -(A*B)+(B**2)+(B**3)
p = beeta-(alpha**2)/3;
q = ((2*alpha**3)/27)-((alpha*beeta)/3)+gaamma
D = (((q)**2)/4)+(((p)**3)/27);
if D > 0:
Z=((-q/2)+(math.sqrt(D)))**(1./3)+((-q/2)-(math.sqrt(D)))**(1./3)-(alpha/3); #One real root given by Eq.(3.32)
Z_l=Z;
Z_v=Z;
elif D==0:
Z1=((-2*(q/2))**(1./3))-(alpha/3); #Three real roots and two equal given by Eq.(3.33)
Z2=((q/2)**(1./3))-(alpha/3);
Z3=((q/2)**(1./3))-(alpha/3);
Z=[Z1, Z2, Z3];
Z_l=min(Z);
Z_v=max(Z);
else:
r = math.sqrt((-(p**3)/27));
theta = math.acos((-(q)/2)*(1./r));
Z1 = (2*(r**(1./3))*math.cos(theta/3))-(alpha/3);
Z2 = (2*(r**(1./3))*math.cos(((2*math.pi)+theta)/3))-(alpha/3)
Z3 = (2*(r**(1./3))*math.cos(((4*math.pi)+theta)/3))-(alpha/3);
Z = [Z1, Z2, Z3];
Z_l = Z[0];
Z_v = Z[1];
phi_l = math.exp(Z_l-1-math.log(Z_l-B)-((a/(2*math.sqrt(2)*b*R*T))*math.log((Z_l+(B*(1+math.sqrt(2))))/(Z_l+(B*(1-math.sqrt(2)))))));
phi_v = math.exp(Z_v-1-math.log(Z_v-B)-((a/(2*math.sqrt(2)*b*R*T))*math.log((Z_v+(B*(1+math.sqrt(2))))/(Z_v+(B*(1-math.sqrt(2)))))));
fl = Prevised*phi_l;
fv = Prevised*phi_v;
Prevised = Prevised*fl/fv
P = Prevised;
# Results
print " The vapour pressure of n-octane at 427.85K = %.5f"%P," MPa"
# Note: answer is slightly differnt because of rounding off error.