Chapter 10 : Stability and phase transition in thermodynamic systems

Example 10.2 Page No : 369

In [1]:
# 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"
 The number of degrees of freedom  =  2 
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

Example 10.3 Page No : 370

In [34]:
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.
 The vapour pressure of n-octane at 427.85K  =  0.21196  MPa