# Chapter 10 : Residual Properties by Equations of State¶

### Example 10.2 Page Number : 334¶

In :
import math

# Variables
T = 40 + 273.15;			#[C] - Temperature
P_1 = 0;			#[bar]
P_2 = 10;			#[bar]
V_liq = 90.45;			#[cm**(3)/mol]
V_liq = V_liq*10**(-6);			#[m**(3)/mol]
P_sat = 4.287;			#[bar]

T_c = 425.0;			#[K] - Critical temperature
P_c = 43.3;			#[bar] - Critical pressure
P_c = P_c*10**(5);			#[N/m**(2)]
w = 0.195;			# Acentric factor
R = 8.314;			#[J/mol*K] - Universal gas constant

# Calculations and Results
# Let us calculate second virial coefficient at 40 C
Tr = T/T_c;			# Reduced temperature
B_0 = 0.083-(0.422/(Tr)**(1.6));
B_1 = 0.139-(0.172/(Tr)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0+(w*B_1)
B = ((B_0 + (w*B_1))*(R*T_c))/P_c;			#[m**(3)/mol] - Second virial coefficient

# math.log(f/P) = (B*P)/(R*T)
# f = P*exp((B*P)/(R*T))

print " The table is as follows"
print " Pbar \t\t fbar \t\t                 phi";
from numpy import zeros
P = [1,2,3,4,4.287,5,6,8,10];
f = zeros(9);
phi = zeros(9);
for i in range(5):
f[i]=P[i]*(math.exp((B*P[i]*10**(5))/(R*T)));			#[bar]  # Pressure inside the exponential term has to be in N/m**(2)
phi[i]= (f[i]/P[i]);
print " %f \t %f \t\t\t %f"%(P[i],f[i],phi[i])

f_sat = f;

# From pressure of 4.287 bar onwards the valid equation to compute fugacity of compressed liquid is given by
# f = f_sat*exp[V_liq*(P-P_sat)/(R*T)]

for j in range(5,9):
f[j] = f_sat*math.exp((V_liq*(P[j]-P_sat)*10**(5))/(R*T));	#[bar]  # Pressure inside the exponential term has to be in N/m**(2)
phi[j] = f[j]/P[j];
print " %f \t %f \t\t\t %f"%(P[j],f[j],phi[j]);

 The table is as follows
Pbar 		 fbar 		                 phi
1.000000 	 0.978336 			 0.978336
2.000000 	 1.914284 			 0.957142
3.000000 	 2.809221 			 0.936407
4.000000 	 3.664484 			 0.916121
4.287000 	 3.902801 			 0.910380
5.000000 	 3.912481 			 0.782496
6.000000 	 3.926097 			 0.654349
8.000000 	 3.953471 			 0.494184
10.000000 	 3.981037 			 0.398104


### Example 10.3 Page Number : 334¶

In :


import math

# Variables
n = 100.;			#[mol] - No of moles
T_1 = 600.;			#[K] - Initial temperature
T_2 = 300.;			#[K] - Final temperature
P_1 = 10.;			#[atm] - Initial pressure
P_1 = P_1*101325;			#[Pa]
P_2 = 5.;			#[atm] - Final presssure
P_2 = P_2*101325;			#[Pa]
Tc = 369.8;			#[K] - Critical temperature
Pc = 42.48;			#[bar] - Critical pressure
Pc = Pc*10**(5);			#[Pa]
w = 0.152;
R = 8.314;			#[J/mol*K] - Universal gas constant

# Calculations
# At 600 K
Tr = T_1/Tc;			# Reduced temperature
B_0 = 0.083-(0.422/(Tr)**(1.6));
B_1 = 0.139-(0.172/(Tr)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B = ((B_0 + (w*B_1))*(R*Tc))/Pc;			#[m**(3)/mol] - Second virial coefficient
dB0_dT = 0.422*1.6*Tc**(1.6)*T_1**(-2.6);			# (dB_0/dT)
dB1_dT = 0.172*4.2*Tc**(4.2)*T_1**(-5.2);			# (dB_1/dT)
dB_dT = ((R*Tc)/(Pc))*(dB0_dT + w*dB1_dT);			# dB/dT

# Now let us calculate B and dB/dT at 300 K
Tr_prime = T_2/Tc;			# Reduced temperature
B_0_prime = 0.083-(0.422/(Tr_prime)**(1.6));
B_1_prime = 0.139-(0.172/(Tr_prime)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B_prime = ((B_0_prime + (w*B_1_prime))*(R*Tc))/Pc;			#[m**(3)/mol] - Second virial coefficient
dB0_dT_prime = 0.422*1.6*Tc**(1.6)*T_2**(-2.6);			# (dB_0/dT)
dB1_dT_prime = 0.172*4.2*Tc**(4.2)*T_2**(-5.2);			# (dB_1/dT)
dB_dT_prime = ((R*Tc)/(Pc))*(dB0_dT_prime + w*dB1_dT_prime);			# dB/dT

# The change in enthalpy for ideal gas is given by

def f16(T):
return -0.966+7.279*10**(-2)*T-3.755*10**(-5)*T**(2)+7.58*10**(-9)*T**(3)

delta_H_ig = delta_H_ig*4.184;			#[J/mol]

# We know that delta_H_ig = delta_U_ig + R*delta_T. Therefore change in internal energy is given by
delta_U_ig = delta_H_ig - R*(T_2 - T_1);			#[J/mol]

# The change in entropy of ideal gas is given by

def f17(T):
return Cp_0/T

def f18(T):
return (-0.966+7.279*10**(-2)*T-3.755*10**(-5)*T**(2)+7.58*10**(-9)*T**(3))/T

# Now let us calculate the change in enthalpy of gas. We know that
# delta_H = delta_H_ig + delta_H_R
# delta_H_R = H_2_R - H_1_R
H_2_R = B_prime*P_2 - P_2*T_2*dB_dT_prime;			# [J/mol]
H_1_R = B*P_1 - P_1*T_1*dB_dT;			# [J/mol]
delta_H_R = H_2_R - H_1_R;			# [J/mol]
delta_H = delta_H_ig + delta_H_R;			# [J/mol]

# Let us calculate the residual entropy of gas
S_2_R = -P_2*dB_dT_prime;			#[J/mol-K]
S_1_R = -P_1*dB_dT;			#[J/mol-K]
delta_S = delta_S_ig + (S_2_R - S_1_R);			#[J/mol-K]

# Let us calculate the residual internal energy of gas
U_2_R = -P_2*T_2*dB_dT_prime;			#[J/mol-K]
U_1_R = -P_1*T_1*dB_dT;			#[J/mol-K]
delta_U = delta_U_ig + (U_2_R - U_1_R);			#[J/mol-K]

# For 100 mol sample,
delta_H_ig = delta_H_ig*n*10**(-3);			#[kJ/mol]
delta_H = delta_H*n*10**(-3);			#[kJ/mol]

delta_U_ig = delta_U_ig*n*10**(-3);			#[kJ/mol]
delta_U = delta_U*n*10**(-3);			#[kJ/mol]

delta_S_ig = delta_S_ig*n*10**(-3);			#[kJ/mol]
delta_S = delta_S*n*10**(-3);			#[kJ/mol]

# Results
print " The value of delta_H = %f kJ/mol"%(delta_H);
print " The value of delta_H_ig ideal gas)= %f J/mol"%(delta_H_ig);
print " The value of delta_U = %f kJ/mol"%(delta_U);
print " The value of delta_U_ig ideal gas) = %f kJ/mol"%(delta_U_ig);
print " The value of delta_S = %f kJ/mol"%(delta_S);
print " The value of delta_S_ig ideal gas) = %f kJ/mol"%(delta_S_ig);

 The value of delta_H = -3130.406247 kJ/mol
The value of delta_H_ig ideal gas)= -3096.763542 J/mol
The value of delta_U = -2867.753839 kJ/mol
The value of delta_U_ig ideal gas) = -2847.343542 kJ/mol
The value of delta_S = -6.466831 kJ/mol
The value of delta_S_ig ideal gas) = -6.358994 kJ/mol


### Example 10.4 Page Number : 337¶

In :

# Variables
T = 35 + 273.15;			#[K] - Temperature
P = 10;			#[atm] - Pressure
P = P*101325;			#[Pa]
# Methane obeys the equation of state
# Z = 1 + (P*B)/(R*T)

# Calculations
# At 35 C,
B = -50;			#[cm**(3)/mol]
dB_dT = 1.0;			#[cm**(3)/mol-K] - dB/dT
dB_dT = dB_dT*10**(-6);			#[m**(3)/mol-K]
d2B_dT2 = -0.01;			#[cm**(3)/mol-K**(2)] - d**2(B)/d(T**2)
d2B_dT2 = d2B_dT2*10**(-6);			#[m**(3)/mol-K**(2)]

# Ideal gas molar heat capacity of methane is given by
# Cp_0 = 4.75 + 1.2*10**(-2)*T + 0.303*10**(-5)*T**(2) - 2.63*10**(-9)*T**(3)

# The molar heat capacity is given by
# Cp = Cp_0 + Cp_R
# For virial gas equation of state
Cp_R = -P*T*d2B_dT2;			#[J/mol-K]

# thus heat capacity is given by
# Cp = a + b*T + c*T**(2) + d*T**(3) - P*T*d2B_dT2
# Putting the values, we get
Cp = (4.75 + 1.2*10**(-2)*T + 0.303*10**(-5)*T**(2) - 2.63*10**(-9)*T**(3))*4.184 - P*T*d2B_dT2;			#[J/mol-K]

# Results
print " The molar heat capacity of methane is %f J/mol-K"%(Cp);

 The molar heat capacity of methane is 39.349753 J/mol-K


### Example 10.5 Page Number : 338¶

In :

from scipy.optimize import fsolve

# Variables
T_1 = 360;			#[K] - Initial temperature
P_1 = 10;			#[bar] - Initial pressure
P_1 = P_1*10**(5);			#[Pa]
Tc = 408.1;			#[K] - Critical temperature
Pc = 36.48;			#[bar] - Critical pressure
Pc = Pc*10**(5);			#[Pa]
w = 0.181;
R = 8.314;			#[J/mol*K] - Universal gas constant
Cv_0 = 106.0;			#[J/mol-K]

# Calculations
# At 360 K
Tr = T_1/Tc;			# Reduced temperature
B_0 = 0.083-(0.422/(Tr)**(1.6));
B_1 = 0.139-(0.172/(Tr)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B = ((B_0 + (w*B_1))*(R*Tc))/Pc;			#[m**(3)/mol] - Second virial coefficient
dB0_dT = 0.422*1.6*Tc**(1.6)*T_1**(-2.6);			# (dB_0/dT)
dB1_dT = 0.172*4.2*Tc**(4.2)*T_1**(-5.2);			# (dB_1/dT)
dB_dT = ((R*Tc)/(Pc))*(dB0_dT + w*dB1_dT);			# dB/dT

# At state 1
V_1 = B + (R*T_1)/P_1;			#[m**(3)/mol] - Molar volume
# At state 1
V_2 = 10*V_1;			#[m**(3)/mol] - Molar volume

T_2 = -(P_1*T_1*(dB_dT))/Cv_0 + T_1;			#[K]

T = 350;

fault = 10;
def f(T):
return  106*(T-T_1)+972.72-((R*T**(2))/(V_2-B_prime))*dB_dT_prime
while(fault>0.007):
Tr_prime = T/Tc;			# Reduced temperature
B_0_prime = 0.083-(0.422/(Tr_prime)**(1.6));
B_1_prime = 0.139-(0.172/(Tr_prime)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B_prime = ((B_0_prime + (w*B_1_prime))*(R*Tc))/Pc;			#[m**(3)/mol] - Second virial coefficient
dB0_dT_prime = 0.422*1.6*Tc**(1.6)*T_2**(-2.6);			# (dB_0/dT)
dB1_dT_prime = 0.172*4.2*Tc**(4.2)*T_2**(-5.2);			# (dB_1/dT)
dB_dT_prime = ((R*Tc)/(Pc))*(dB0_dT_prime + w*dB1_dT_prime);			# dB/dT
T_prime = fsolve(f,0.15)
fault=abs(T-T_prime);
T = T + 0.001;

# Results
print " The final temperature is %f K"%(T);

 The final temperature is 351.910000 K


### Example 10.6 Page Number : 339¶

In :

# Variables
T = 220 + 273.15;			#[K] - Temperature
Tc = 562.2;			#[K] - Critical temperature
Pc = 48.98;			#[bar] - Critical pressure
Pc = Pc*10**(5);			#[Pa]
w = 0.210;
R = 8.314;			#[J/mol*K] - Universal gas constant
P_sat = 1912.86;			#[kPa] - Saturation pressure at 220 C
P_sat = P_sat*10**(3);			#[Pa]
Mol_wt = 78.114;			#[g/mol] - Molecular weight of benzene

# Calculations and Results
#(1)

# At 220 C
Tr = T/Tc;			# Reduced temperature
B_0 = 0.083-(0.422/(Tr)**(1.6));
B_1 = 0.139-(0.172/(Tr)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B = ((B_0 + (w*B_1))*(R*Tc))/Pc;			#[m**(3)/mol] - Second virial coefficient

# We know that math.log(f/P) = (B*P)/(R*T)
# Thus at saturated conditions
# math.log(f_sat/P_sat) = B*P_sat/(R*T)
f_sat = P_sat*(math.exp((B*P_sat)/(R*T)));			#[Pa]
f_sat = f_sat*10**(-3);			#[kPa]

print " 1).The fugacity of liquid benzene is %f kPa"%(f_sat);

#(2)
P = 2014.7;			# [psia] - Total gauge pressure
P = 138.94;			# [bar]
P = P*10**(5);			# [Pa]
den = 0.63;			# [g/cm**(3)] - density of benzene
den = den*10**(3);			# [kg/m**(3)]

# Therefore specific volume is
V = 1/den;			#[m/**(3)/kg]
# Molar volume is given by
V = V*Mol_wt*10**(-3);			#[m**(3)/mol]

# Thus fugacity at 220 C and pressure P is given by
f = f_sat*(math.exp((V*(P-P_sat))/(R*T)));

print " 2).The fugacity of liquid benzene is %f kPa"%(f);

 1).The fugacity of liquid benzene is 1551.083216 kPa
2).The fugacity of liquid benzene is 2228.386532 kPa


### Example 10.7 Page Number : 341¶

In :

# Variables
# C = -0.067 + 30.7/T
# D = 0.0012 - 0.416/T

T = 80 + 273.15;			#[K]
P = 30;			#[bar]
#P = P;			#[N/m**(2)]
R = 8.314;			#[J/mol*K] - Universal gas constant

# Calculations
# We have the relation derived in the book
# d(G/(R*T)) = (V/(R*T))dP - (H/(R*T**(2)))dT
# Writing the same equation for ideal gas and subtracting it from the above equation we get
# d(G_R/(R*T)) = (V_R/(R*T))dP - (H_R/(R*T**(2)))dT
# Therefore, H_R/(R*T**(2)) = -[del((G_R)/(R*T))/del(T)]_P

# Substituting the relation G_R/(R*T) = math.log(f/P), we get
# H_R/(R*T**(2)) = -[del(math.log(f/P))/del(T)]_P = -[del(-C*P - D*P**(2))/del(T)]_P
# H_R/R = - 30.7*P + 0.416*P**(2)

# Substituting the given conditions we get
H_R = R*(-30.7*P + 0.416*P**(2));			#[J/mol]

# Results
print " The molar enthalpy of the gas relative to that of the ideal gas at 80 C and 30 bar pressure is H_R = %f J/mol"%(H_R);

 The molar enthalpy of the gas relative to that of the ideal gas at 80 C and 30 bar pressure is H_R = -4544.432400 J/mol


### Example 10.8 Page Number : 341¶

In :

# Variables
# (1)
T = 311;			#[K] - Temperature
R = 8.314;			#[J/mol*K] - Universal gas constant
# Pressure in 'bar' is given below
P = [0.690,1.380,2.760,5.520,8.280,11.034,13.800,16.550];
# Molar volume in 'm**(3)/mol' is given below
V = [0.0373,0.0186,0.00923,0.00455,0.00298,0.00220,0.00175,0.00144];

# Z = 1 + (B/V) + (C/V**(2))
# (Z-1)*V = B + (C/V)

from numpy import zeros
Z=zeros(8);
k=zeros(8);
t=zeros(8);

# Calculations and Results
for i in range(8):
Z[i]=(P[i]*10**(5)*V[i])/(R*T);
k[i]=(Z[i]-1)*V[i];
t[i]=1./V[i];

from scipy.stats import linregress
from numpy import array
t = array(t)
k = array(k)
C,B,c,d,e = linregress(t.T,k.T)
#[C,B,sig]=reglin(t',k');

#From the regression, we get intercept = B and slope = C,and thus,

print " 1).The second virial coefficient of CO2 is given by B = %e m**3)/mol"%(B);
print "     The thied virial coefficient of CO2 is given by C = %e m**6)/mol**2)"%(C);

# (2)
P_final = 13.8;			#[bar]

def f40(P):
return V-(R*T)/P

# We know that R*T*math.log(f/P) =  quad(f40,0,P)

# Therefore we have to plot V - (R*T)/P versus P and calculate the area beneath the curve from 0 to 13.8 bar
# For this we need the value of the term V - (R*T)/P at P = 0. At low pressure the virial equation becomes
# Z = 1 + (B/V)
#and V - (R*T)/P = (Z*R*T)/P - (R*T)/P = (1 + (B/V))*((R*T)/P) - (R*T)/P = (B*R*T)/(P*V) = (B/Z)
# Thus lim P tending to zero (V - (R*T)/P) = B    ( as P tend to zero, Z tend to 1 )

P_prime = [0.000,0.690,1.380,2.760,5.520,8.280,11.034,13.800];
V_prime = [0.000,0.0373,0.0186,0.00923,0.00455,0.00298,0.00220,0.00175];
summation = 0;
x=zeros(8);
y=zeros(8);
z=zeros(8);
for j in range(2,8):
x[j]=V_prime[j]-(R*T)/(P_prime[j]*10**(5));			#[m**(3)/mol]
y[j]=(x[j] + x[j-1])/2;
z[j]=y[j]*((P_prime[j]-P_prime[j-1]))*10**(5);
summation = summation + z[j] ;			#[J/mol]

summation = summation + 2*z - z;			# Because in the above calculation,in order to calculate the average a summation of z(2) is not included,only half of it gets added

def f41(P):
return V - (R*T)/P

# Now we have, area =  quad(f41,0,13.8*10**(5))

# R*T*math.log(f/P) = summation
f = P_final*(math.exp(summation/(R*T)));			#[bar]

print " 2).The fugacity of steam at 311 K and 13.8 bar pressure is %f bar"%(f);

 1).The second virial coefficient of CO2 is given by B = -1.515264e-04 m**3)/mol
The thied virial coefficient of CO2 is given by C = 5.618543e-08 m**6)/mol**2)
2).The fugacity of steam at 311 K and 13.8 bar pressure is 12.892780 bar


### Example 10.9 Page Number : 344¶

In :

# Variables
T = 0 + 273.15;			#[K] - Temperature
R = 8.314;			#[J/mol*K] - Universal gas constant
# Pressure in 'atm' is given below
P = [100,200,300,400,500,600,700,800,900,1000];
# The compressibility factor values are
Z = [1.069,1.138,1.209,1.283,1.356,1.431,1.504,1.577,1.649,1.720];

# Z = 1 + (B/V) + (C/V**(2))
# (Z-1)*V = B + (C/V)
from numpy import zeros

V = zeros(10);
k = zeros(10);
t = zeros(10);

# Calculations and Results
for i in range(10):
V[i]=Z[i]*R*T/(P[i]*101325);			#[m**(3)/mol]
k[i]=(Z[i]-1)*V[i];
t[i]=1/V[i];

from scipy.stats import linregress
C,B,c,d,e = linregress(t,k)
#[C,B,sig]=reglin(t,k);

#From the regression, we get intercept = B and slope = C,and thus,

print " 1).The second virial coefficient of H2 is given by B = %e m**3)/mol"%(B);
print "     The thied virial coefficient of H2 is given by C = %e m**6)/mol**2)"%(C);

# (2)
# We know that, limit P tending to zero (V-(R*T)/P) = B, therfore P = 0,  V-(R*T)/P = B

def f42(P):
return V-(R*T)/P 			#and determine the integral integrate((V-(R*T)/P)

# Now let us tabulate V-(R*T)/P and determine the integral  quad(f42,0,1000)

P_prime = [0,100,200,300,400,500,600,700,800,900,1000];
Z_prime = [0,1.069,1.138,1.209,1.283,1.356,1.431,1.504,1.577,1.649,1.720];

summation = 0;
V_prime = zeros(11);
x = zeros(11);
y = zeros(11);
z = zeros(11);
for j in range(1,11):
V_prime[j]=Z_prime[j]*R*T/(P_prime[j]*101325);			#[m**(3)/mol]
x[j]=V_prime[j]-(R*T)/(P_prime[j]*101325);
y[j]=(x[j] + x[j-1])/2;
z[j]=y[j]*((P_prime[j]-P_prime[j-1]))*101325;
summation = summation + z[j] ;			#[J/mol]

summation = summation + 2*z - z;

# Now we have
# R*T*math.log(f/P) = summation
P_dash = 1000;			#[atm] - Pressure at which fugacity is to be calculated
T_dash = 273.15;			#[K] - Temperature at which fugacity is to be calculated
f = P_dash*math.exp(summation/(R*T_dash));			#[atm]

print " 2).The fugacity of H2 at 0 C and 1000 atm pressure is  f = %f atm"%(f);

 1).The second virial coefficient of H2 is given by B = 1.345242e-05 m**3)/mol
The thied virial coefficient of H2 is given by C = 5.286525e-10 m**6)/mol**2)
2).The fugacity of H2 at 0 C and 1000 atm pressure is  f = 2030.305169 atm


### Example 10.10 Page Number : 345¶

In :

from scipy.optimize import fsolve
import math

# Variables
P_1 = 1*10.**(6.);			#[Pa] - Initial pressure
T_1 = 200. + 273.15;		#[K] - Initial temperature
P_2 = 8*10.**(6);			#[Pa] - Final pressure
Tc = 647.1;			        #[K] - Critical temperature of water
Pc = 220.55;    			#[bar] - Critical pressure of water
Pc = Pc*10**(5);			#[Pa]
w = 0.345;
R = 8.314;		    	#[J/mol*K] - Universal gas constant

# For the virial gas the following are the relations for residual enthalpy and entropy
# H_R = B*P - P*T*(dB/dT)
# S_R = -P*(dB/dT)
# Where, (dB/dT) = ((R*Tc)/Pc)*((dB_0/dT) + w*(dB_1/dT))
# dB0_dT = 0.422*1.6*Tc**(1.6)*T**(-2.6);			# (dB_0/dT)
# dB1_dT = 0.172*4.2*Tc**(4.2)*T**(-5.2);			# (dB_1/dT)

# (1)
Cp_0 = 29.114;			#[J/mol-K] - Specific heat capacity at constant pressure
# For the isentropic process entropy change is zero, thus
# delta_S = Cp_0*math.log(T_2/T_1) - P_2*(dB/dT)_2 + P_1*(dB/dT)_1 = 0

# Calculations
# At state 1,
Tr_1 = T_1/Tc;
B0_1 = 0.083 - 0.422/(Tr_1**(1.6));
B1_1 = 0.139 - 0.172/(Tr_1**(4.2));
# (B*Pc)/(R*Tc) = B0 + w*B1
B_1 = ((B0_1 + w*B1_1)*(R*Tc))/Pc;			# [m**(3)/mol] - Second virial coefficient at state 1
dB0_dT_1 = 0.422*1.6*Tc**(1.6)*T_1**(-2.6);			# (dB_0/dT)
dB1_dT_1 = 0.172*4.2*Tc**(4.2)*T_1**(-5.2);			# (dB_1/dT)
dB_dT_1 = ((R*Tc)/Pc)*((dB0_dT_1) + w*(dB1_dT_1));			# (dB/dT)_1

# Now let us assume the exit temperature to be 870 K, at this temperature
# T_2 = 870;			#[K] -
# At this temperature
# delta_S = Cp_0*math.log(T_2/T_1) - P_2*(dB/dT)_2 + P_1*(dB/dT)_1 =

T_2 = 860.;			#[K] - Exit temperature
# Therefore at state 2, we have
Tr_2 = T_2/Tc;
B0_2 = 0.083 - 0.422/(Tr_2**(1.6));
B1_2 = 0.139 - 0.172/(Tr_2**(4.2));
# (B*Pc)/(R*Tc) = B0 + w*B1
B_2 = ((B0_2 + w*B1_2)*(R*Tc))/Pc;			# [m**(3)/mol] - Second virial coefficient at state 2
dB0_dT_2 = 0.422*1.6*Tc**(1.6)*T_2**(-2.6);			# (dB_0/dT)
dB1_dT_2 = 0.172*4.2*Tc**(4.2)*T_2**(-5.2);			# (dB_1/dT)
dB_dT_2 = ((R*Tc)/Pc)*((dB0_dT_2) + w*(dB1_dT_2));			# (dB/dT)_2

delta_H_s = Cp_0*(T_2 - T_1) + B_2*P_2 -P_2*T_2*(dB_dT_2) - B_1*P_1 + P_1*T_1*(dB_dT_1);			#[J/mol] - Enthalpy change

# As no heat exchange is assumed to take place with the surroundings,work transfer is given by
W_1 = - delta_H_s;			# [J/mol]

# Results
print " 1).The exit temperature is %f K"%(T_2);
print "     The required amount of work is %f J/mol"%(W_1);

# (2)
eff = 0.8;			# Adiabatic efficiency
delta_H_a = delta_H_s/0.8;			# Actual enthalpy change

# Now for calculating the value of T_exit
# delta_H_a = Cp_0*(T_exit - T_1) + B*P_2 -P_2*T_exit*(dB_dT) - B_1*P_1 + P_1*T_1*(dB_dT_1)
# On simplification we get
# 29.114*(T_2 - T_1)*B_2*8*10**(6)-8*10**(6)*T_2*(dB/dT)_2 = 12643.77

# Let us assume a temperature of say
T = 900.;			#[K]
fault=10.;

def f(T_exit):
return  delta_H_a - Cp_0*(T_exit - T_1) + B*P_2 -P_2*T_exit*(dB_dT) - B_1*P_1 + P_1*T_1*(dB_dT_1)

while(fault>0.3):
Tr = T/Tc;
B0 = 0.083 - 0.422/(Tr**(1.6));
B1 = 0.139 - 0.172/(Tr**(4.2));
# (B*Pc)/(R*Tc) = B0 + w*B1
B = ((B0 + w*B1)*(R*Tc))/Pc;			# [m**(3)/mol] - Second virial coefficient at state 2
dB0_dT = 0.422*1.6*Tc**(1.6)*T**(-2.6);			# (dB_0/dT)
dB1_dT = 0.172*4.2*Tc**(4.2)*T**(-5.2);			# (dB_1/dT)
dB_dT = ((R*Tc)/Pc)*((dB0_dT) + w*(dB1_dT));			# (dB/dT)_1

T_exit = fsolve(f,900)
fault=abs(T-T_exit);
T = T + 0.2;

Texit = T;

W_2 = - delta_H_a;			# [J/mol]

print " 2).The exit temperature is %f K"%(Texit);
print "     The required amount of work is %f J/mol"%(W_2);

def f20(T):
return Cp_0/T

T_prime = 700.;			#[K]
fault1=10.;

def f1(T_out):
return  32.2168*math.log(T_out) + 0.1922*10**(-2)*T_out + 0.5274*10**(-5) \
*T_2**(2) - 1.1976*10**(-9)*T_out**(3)-8*10**(6)*dB_dT_prime -216.64

while(fault1>0.5):
Tr_prime = T_prime/Tc;
B0_prime = 0.083 - 0.422/(Tr_prime**(1.6));
B1_prime = 0.139 - 0.172/(Tr_prime**(4.2));
B_prime = ((B0_prime + w*B1_prime)*(R*Tc))/Pc;			# [m**(3)/mol] - Second virial coefficient at state 2
dB0_dT_prime = 0.422*1.6*Tc**(1.6)*T_prime**(-2.6);			# (dB_0/dT)
dB1_dT_prime = 0.172*4.2*Tc**(4.2)*T_prime**(-5.2);			# (dB_1/dT)
dB_dT_prime = ((R*Tc)/Pc)*((dB0_dT_prime) + w*(dB1_dT_prime));			# (dB/dT)_1
T_out = fsolve(f1,10)
fault1=abs(T_prime-T_out);
T_prime = T_prime + 0.5;

T_out = T_prime;

# Now we have to calculate enthalpy change as W = -delta_H

def f21(T):
return (7.7 + 0.04594*10**(-2.)*T + 0.2521*10**(-5.)*T**(2.) - 0.8587*10.**(-9.)*T**(3.))*4.184

delta_H_3 =  quad(f21,T_1,T_out) + B_prime*P_2 - P_2*T_out*dB_dT_prime - B_1*P_1 + P_1*T_1*dB_dT_1
print T_1,T_out

W_3 = - delta_H_3;			# [J/mol]

print " 3).The exit temperature is %f K"%(T_out);
print "     The required amount of work is %f J/mol"%(W_3);

#(4)
n = 0.8;			# Adiabatic efficiency
delta_H_a_4 = delta_H_3/n;			#[J/mol]
W_4 = -delta_H_a_4;			#[J/mol]

T_prime1 = 700.;			#[K]
fault2=10.;

def f2(T_2):
return  7.7*4.184*(T_2-T_1) + ((0.04594*4.184*10**(-2))/2)* \
(T_2**(2)-T_1**(2)) + ((0.2521*4.184*10**(-5))/3)*(T_2**(3)-T_1**(3)) - \
((0.8587*4.184*10**(-9))/4)*(T_2**(4)-T_1**(4)) + B_prime1*8*10**(6) \
- 8*10**(6)*T_2*dB_dT_prime1 + 191.7 + 496.81 - delta_H_a_4

while(fault2>0.5):
Tr_prime1 = T_prime1/Tc;
B0_prime1 = 0.083 - 0.422/(Tr_prime1**(1.6));
B1_prime1 = 0.139 - 0.172/(Tr_prime1**(4.2));
# (B*Pc)/(R*Tc) = B0 + w*B1
B_prime1 = ((B0_prime1 + w*B1_prime1)*(R*Tc))/Pc;			# [m**(3)/mol] - Second virial coefficient at state 2
dB0_dT_prime1 = 0.422*1.6*Tc**(1.6)*T_prime1**(-2.6);			# (dB_0/dT)
dB1_dT_prime1 = 0.172*4.2*Tc**(4.2)*T_prime1**(-5.2);			# (dB_1/dT)
dB_dT_prime1 = ((R*Tc)/Pc)*((dB0_dT_prime1) + w*(dB1_dT_prime1));			# (dB/dT)_1
T_out1 = fsolve(f2,100)
fault2=abs(T_prime1-T_out1);
T_prime1 = T_prime1 + 0.5;

T_out1 = T_prime1;

print " 4).The exit temperature is %f K"%(T_out1);
print "     The required amount of work is %f J/mol"%(W_4);

# answers are varies because of rounding error. please check it manually.

 1).The exit temperature is 860.000000 K
The required amount of work is -10667.724150 J/mol
2).The exit temperature is 916.600000 K
The required amount of work is -13334.655188 J/mol
473.15 755.5
3).The exit temperature is 755.500000 K
The required amount of work is -9282.701079 J/mol
4).The exit temperature is 809.500000 K
The required amount of work is -11603.376349 J/mol


### Example 10.11 Page Number : 348¶

In :

# Variables
Vol = 0.15;			    #[m**(3)] - Volume of the cylinder
P_1 = 100.;			    #[bar] - Initial pressure
P_1 = P_1*10**(5);		#[Pa]
T_1 = 170.;			    #[K] - Initial temperature
n_withdrawn = 500.;		#[mol] - Withdrawn number of moles
R = 8.314;		    	#[J/mol*K] - Universal gas constant

# Calculations and Results
#(1)
Y = 1.4;			# Coefficient of adiabatic expansion
n_total = (P_1*Vol)/(R*T_1);			#[mol] - Total number of moles
n_2 = n_total - n_withdrawn;			#[mol] - Left number of moles
V_1 = Vol/n_total;			#[m**(3)/mol] - Molar volume at initial state.
# At final state
V_2 = Vol/n_2;			#[m**(3)/mol] - Molar volume at final state

# During duscharging  P_1*V_1**(Y) = P_2*V_2**(Y), therefore
P_2_1 = P_1*((V_1/V_2)**(Y));			#[Pa] - Final pressure
P_2_1 = P_2_1*10**(-5);			#[bar]
T_2_1 = ((P_2_1*10**(5))*V_2)/R;			#[K] - Final temperature

print " 1).The final temperature %f K"%(T_2_1);
print "     The final pressure %f bar"%(P_2_1);

def f53(T):
return Cp_0/T

T_prime = 150;			#[K]
error = 10;
while(error>1):
f_T = 18.886*math.log(T_prime) + 4.2*10**(-3)*T_prime - 92.4;
f_T_dash = 18.886/T_prime + 4.2*10**(-3);
T_new = T_prime - (f_T/f_T_dash);
error=abs(T_prime - T_new);
T_prime = T_new;

T_2_2 = T_prime;			#[K] - Final temperature
P_2_2 = ((n_2*R*T_2_2)/Vol)*10**(-5);			#[bar] - Final pressure

print " 2).The final temperature %f K"%(T_2_2);
print "     The final pressure %f bar"%(P_2_2);

#(3)
Tc = 126.2;			#[K] - Critical temperature of nitrogen
Pc = 34.0;			#[bar] - Critical pressure of nitrogen
Pc = Pc*10**(5);			#[Pa]
w = 0.038;			# Acentric factor

dB0_dT = 0.422*1.6*Tc**(1.6)*T_1**(-2.6);			# (dB_0/dT) at state 1
dB1_dT = 0.172*4.2*Tc**(4.2)*T_1**(-5.2);			# (dB_1/dT) at state 1
dB_dT = ((R*Tc)/Pc)*((dB0_dT) + w*(dB1_dT));			# (dB/dT) at state 1
# The residual entropy at the initial state is given by
S_R_1 = -P_1*(dB_dT);			#[J/mol-K]

# Now let us calculate molar volume at initial state
Tr = T_1/Tc;			# Reduced temperature
B_0 = 0.083-(0.422/(Tr)**(1.6));
B_1 = 0.139-(0.172/(Tr)**(4.2));

#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B = ((B_0+(w*B_1))*(R*Tc))/Pc;			#[m**(3)/mol]

V_1_3 = B + (R*T_1)/P_1;			#[m**(3)/mol]
# Therefore number of moles in the initial state is
n_1_3 = Vol/V_1_3;			#[mol]
# Therefore final number of moles is
n_2_3 = n_1_3 - n_withdrawn;

# Therefore molar volume at final state is
V_2_3 = Vol/n_2_3;			#[m**(3)/mol]

# Now let us determine the relation between pressure and temperature in the final state
# P_2_3 = (R*T_2_3)/(V_2_3 - B_2)
#delta_S = 0, thus delta_S_ig + delta_S_R = 0
delta_S_R = - S_R_1;

def f54(T):
return Cp_0/T

T_2_prime = 135;			#[K]
delta = 0.1;
error = 10;
while(error>0.01):
T_r = T_2_prime/Tc;			# Reduced temperature
B_0_3 = 0.083-(0.422/(T_r)**(1.6));
B_1_3 = 0.139-(0.172/(T_r)**(4.2));
B_3 = ((B_0_3+(w*B_1_3))*(R*Tc))/Pc;			#[m**(3)/mol]
dB0_dT_3 = 0.422*1.6*Tc**(1.6)*T_2_prime**(-2.6);			# (dB_0/dT)
dB1_dT_3 = 0.172*4.2*Tc**(4.2)*T_2_prime**(-5.2);			# (dB_1/dT)
dB_dT_3 = ((R*Tc)/Pc)*((dB0_dT_3) + w*(dB1_dT_3));			# (dB/dT)
P_2_3 = (R*T_2_prime)/(V_2_3 - B_3);
delta_S = 27.2*(math.log(T_2_prime/T_1)) + 4.2*10**(-3)*(T_2_prime - T_1) - R*(math.log(P_2_3/P_1)) - P_2_3*(dB_dT_3) + delta_S_R;
T_new = T_2_prime + delta;
error=abs(delta_S);
T_2_prime = T_new;

T_2_3 = T_2_prime;			#[K] - Final temperature
# Therefore at T_2_3
P_2_3 = P_2_3*10**(-5);			#[bar] - Final pressure

print " 3).The final temperature %f K"%(T_2_3);
print "     The final pressure %f bar"%(P_2_3);

 1).The final temperature 131.761821 K
The final pressure 40.991361 bar
2).The final temperature 129.504151 K
The final pressure 40.288995 bar
3).The final temperature 141.300000 K
The final pressure 57.079997 bar


### Example 10.12 Page Number : 351¶

In :

# Variables
P_1 = 80.;			        #[bar] - Initial pressure
P_1 = P_1*10.**(5);			#[Pa]
T_1 = 300. + 273.15;		#[T] - Initial temperature
P_2 = 40.;			        #[bar] - Final pressure
P_2 = P_2*10**(5);			#[Pa]
T_2 = 300. + 273.15;		#[K] - Final temperature
T_0 = 25. + 273.15;			#[K] - Surrounding temperature
P_0 = 1.;			        #[atm] - Surrounding pressure
P_0 = P_0*101325;			#[Pa]
Tc = 647.1;			        #[K]
Pc = 220.55;			    #[bar]
Pc = Pc*10.**(5);			#[Pa]
R = 8.314;			        #[J/mol*K] - Universal gas constant

# Calculations
# For van der Walls equation of state
a = (27*R**(2)*Tc**(2))/(64*Pc);			#[Pa-m**(6)/mol**(2)]
b = (R*Tc)/(8*Pc);			#[m**(3)/mol]

def f(V):
return V**(3)-(b+(R*T_1)/P_1)*V**(2)+(a/P_1)*V-(a*b)/P_1
V_1_1=fsolve(f,0.1)
V_1_2=fsolve(f,10)
V_1_2=fsolve(f,100)
# The largest root is considered because of vapour
V_1 = V_1_1;

U_R_1 = -a/V_1;			#[J/mol] - Internal energy
H_R_1 = P_1*V_1 - R*T_1 - a/V_1;			#[J/mol] - Enthalpy
S_R_1 = R*math.log((P_1*(V_1-b))/(R*T_1));

def f1(V):
return V**(3)-(b+(R*T_2)/P_2)*V**(2)+(a/P_2)*V-(a*b)/P_2
V_2_1 = fsolve(f1,0.1)
V_2_2 = fsolve(f1,10)
V_2_3 = fsolve(f1,100)
V_2 = V_2_1;

U_R_2 = -a/V_2;			#[J/mol] - Internal energy
H_R_2 = P_2*V_2 - R*T_2 - a/V_2;			#[J/mol] - Enthalpy
S_R_2 = R*math.log((P_2*(V_2-b))/(R*T_2));

delta_U_R = U_R_2 - U_R_1;			#
delta_H_R = H_R_2 - H_R_1;			#
delta_S_R = S_R_2 - S_R_1;			#

delta_U_ig = 0;			#[J/mol] - As temperature is constant
delta_H_ig = 0;			#[J/mol] - As temperature is constant
delta_S_ig = - R*math.log(P_2/P_1);			# [J/mol-K]
delta_U = delta_U_R + delta_U_ig;			#[J/mol]
delta_H = delta_H_R + delta_H_ig;			#[J/mol]
delta_S = delta_S_R + delta_S_ig;			#[J/mol-K]

# Change in exergy is given by
# delta_phi = phi_1 - phi_2 = U_1 - U_2 + P_0*(V_1 - _V_2) - T_0*(S_1 - S_2)
delta_phi = - delta_U + P_0*(V_1 - V_2) - T_0*(-delta_S);			#[J/mol]

# Results
print " The change in internal energy is %f J/mol"%(delta_U);
print " The change in enthalpy is %f J/mol"%(delta_H);
print " The change in entropy is %f J/mol-K"%(delta_S);
print " The change in exergy is %f J/mol"%(delta_phi);

 The change in internal energy is 615.070900 J/mol
The change in enthalpy is 1053.220316 J/mol
The change in entropy is 6.930264 J/mol-K
The change in exergy is 1389.940735 J/mol


### Example 10.13 Page Number : 353¶

In :

# Variables
T_1 = 500.;			#[K] - Initial temperature
P_1 = 30.;			#[atm] - Initial pressure
P_1 = P_1*101325;			#[Pa]
P_2 = 1;			#[atm] - Final pressure
P_2 = P_2*101325;			#[Pa]
R = 8.314;			#[J/mol*K] - Universal gas constant
# For chlorine
Tc = 417.2;			#[K] - Critical temperature
Pc = 77.10;			#[bar] - Critical pressure
Pc = Pc*10**(5);			#[Pa]

#Redlich Kwong equation of state,
a = (0.42748*(R**(2))*(Tc**(2.5)))/Pc;			# [Pa*m**(6)*K**(1/2)/mol]
b = (0.08664*R*Tc)/Pc;			# [m**(3)/mol]

# Calculations
def f1(V):
return V**(3)-((R*T_1)/P_1)*V**(2)-((b**(2))+((b*R*T_1)/P_1)-(a/(T_1**(1./2)*P_1)))*V-(a*b)/(T_1**(1./2)*P_1)
V_1=fsolve(f1,1)
V_2=fsolve(f1,10)
V_3=fsolve(f1,100)

V = V_1;			#[m**(3)/mol]

Z = (P_1*V_1)/(R*T_1);			#compressibility factor

# The residual enthalpy at state 1 is given by
H_R_1 = (Z-1)*R*T_1 + ((3*a)/(2*b*T_1**(1./2)))*(math.log(V/(V+b)));			#[J/mol]

H_R_2 = 0;			# Residual enthalpy at state 2
delta_H_R = H_R_2 - H_R_1;			#[J/mol] - Residual enthalpy change
# and since isothermal conditions are maintained, therfore
delta_H_ig = 0;			# Enthalpy change under ideal condition
delta_H = delta_H_R + delta_H_ig;			#[J/mol]

# Results
print " The change in enthalpy is given by delta_H = %f J/mol"%(delta_H);

 The change in enthalpy is given by delta_H = 1053.558471 J/mol


### Example 10.14 Page Number : 353¶

In :


# Variables
Vol_1 = 0.1;			#[m**(3)] - Initial volume of each compartment
n_1 = 400;			#[mol] - Initial number of moles in compartment 1
V_1 = Vol_1/n_1;			#[m**(3)/mol] - Molar volume at state 1
T_1 = 294;			#[K]
Vol_2 = 0.2;			#[m**(3)] - Final volume of the compartment after removing the partition.
n_2 = n_1;			#[mol] - Number of moles remains the same
V_2 = Vol_2/n_2;			#[m**(3)/mol] - Molar volume at state 2

a = 0.1362;			#[Pa-m**(6)/mol**(2)]
b = 3.215*10**(-5);			#[m**(3)/mol]
Cv_0 = 12.56;			#[J/mol-K] - Heat capacity in ideal gas state

# Calculations
# For overall system q = 0, and no work is done, therefore delta_U = 0
# Therfore from the relation proved in part (1), we have
T_2 = T_1 + (a/Cv_0)*(1/V_2 - 1/V_1);			#[K]

# Results
print " 2).The final temperatutre is %f K"%(T_2)

 2).The final temperatutre is 272.312102 K


### Example 10.17 Page Number : 356¶

In :

# Variables
P_1 = 1*10**(6);			#[Pa] - Initial pressure
T_1 = 200 + 273.15;			#[K] - Initial temperature
P_2 = 8*10**(6);			#[Pa]
R = 8.314;			#[J/mol*K] - Universal gas constant
Y = 1.4;			# Index of expansion
Cp_0 = 29.114;			#[J/mol-K]

a = 0.55366;			#[Pa-m**(6)/mol**(2)]
b = 3.049*10**(-5);			#[m**(3)/mol]

# Calculations
V_1 = 3.816*10**(-3);			#[m**(3)/mol]
Z_1 = (P_1*V_1)/(R*T_1);

T_2 = T_1*(P_2/P_1)**((Y-1)/Y);			#[K]

# At 8 MPa and T_2,
# The molar volume of steam following van der Walls equation of state (as reported in the book) is
V_2 = 8.41*10**(-4);			#[m**(3)/mol]
# And the compressibility factor is
Z_2 = (P_2*V_2)/(R*T_2);

# For van der Walls equation of state we know that
# delta_S_R/R = math.log(Z_2/Z_1) + math.log((V_2 - b)/V_2) - math.log((V_1 - b)/V_1)
delta_S_R = R*(math.log(Z_2/Z_1) + math.log((V_2 - b)/V_2) - math.log((V_1 - b)/V_1));			#[J/mol]

# delta_S_ig = Cp_0*math.log(T_2/T_1) - R*math.log(P_2/P_1)
# The entropy change is therefore
# delta_S = delta_S_ig + delta_S_R

# Let us assume a temperature, say T = 870 K
# At 870 K the molar volume of steam following van der Walls equation of state (as reported in the book) is
#  V_3 = 8.57*10**(-4);			# [m**(3)/mol]
# Therefore
# Z_3 = (P_2*V_3)/(R*T_2);
# At this temperature,
# delta_S = Cp_0*math.log(T/T_1) - R*math.log(P_2/P_1) + R*(math.log(Z/Z_1) + R*math.log((V - b)/V) - R*math.log((V_1 - b)/V_1))

T = 800;			#[K]
fault=10;

def f1(V):
return V**(3)-(b+(R*T)/P_2)*V**(2)+(a/P_2)*V-(a*b)/P_2

def f2(T):
return  Cp_0*math.log(T/T_1) - R*math.log(P_2/P_1) + R*(math.log(Z/Z_1) + R*(math.log((V - b)/V)) - R*(math.log((V_1 - b)/V_1)))

while(fault>0.3):
# At T and 8 MPa
V = fsolve(f1,1)
Z = (P_2*V)/(R*T);

T_exit = fsolve(f2,0.1)
fault=abs(T-T_exit);
T = T + 0.5;

Texit = T;

# W = - delta_H

# For van der Walls gas the enthalpy change is given by
delta_H_s = Cp_0*(T_exit - T_1) + (Z - 1)*R*T_exit - a/V - (Z_1-1)*R*T_1 + a/V_1;			#[J/mol]
W = - delta_H_s;			#[J/mol]

# Results
print " 1).The exit temperature is %f K"%(Texit);
print "     The work required is given by W = %f J/mol"%(W);

#(2)
eff = 0.8;			# Adiabatic efficiency
delta_H_a = eff*delta_H_s;			#[J/mol] - Actual enthalpy change
W_2 = - delta_H_a;

# Let us assume a temperature, say
T_prime= 900;			#[K]
fault1=10;
def f22(V):
return V**(3)-(b+(R*T_prime)/P_2)*V**(2)+(a/P_2)*V-(a*b)/P_2

def f3(T_prime):
return  Cp_0*(T_prime - T_1) + (Z_prime - 1)*R*T_prime - a/V_prime - 13230.49

while(fault1>0.3):
# At T_prime and 8 MPa
V_prime=fsolve(f22,1)
Z_prime = (P_2*V_prime)/(R*T_prime);
T_exit1 = fsolve(f3,100)
fault1=abs(T_prime-T_exit1);
T_prime = T_prime + 0.2;

Texit1 = T_prime;

print " 2).The exit temperature is %f K"%(Texit1);
print "     The work required is given by W = %f J/mol"%(W_2);

 1).The exit temperature is 916.500000 K
The work required is given by W = -12195.093996 J/mol
2).The exit temperature is 958.400000 K
The work required is given by W = -9756.075197 J/mol


### Example 10.19 Page Number : 359¶

In :

# Variables
T = 100 + 273.15;			#[K] - Temperature
Tc = 647.1;			#[K] - Critical temperature of water
Pc = 220.55;			#[bar] - Critical pressure of water
Pc = Pc*10**(5);			#[Pa]
R = 8.314;			#[J/mol*K] - Universal gas constant

# Calculations
a = (27*R**(2)*Tc**(2))/(64*Pc);			#[Pa-m**(6)/mol**(2)]
b = (R*Tc)/(8*Pc);			#[m**(3)/mol]

# The cubic form of van der Walls equation of state is given by,
# V**(3) - (b + (R*T)/P)*V**(2) + (a/P)*V - (a*b)/P = 0

# For water vapour at 100 C under saturated conditions pressure is 1 atm, therefore
P = 1;			#[atm]
P = P*101325;			#[Pa]

# At 100 C and 1 atm
def f(V):
return V**(3)-(b+(R*T)/P)*V**(2)+(a/P)*V-(a*b)/P
V_1 = fsolve(f,0.1)
V_1 = fsolve(f,10)
V_1 = fsolve(f,100)
V = V_1;			#[m**(3)/mol]

f = P*(math.exp(math.log((R*T)/(P*(V-b))) + b/(V-b) - (2*a)/(R*T*V)));			#[Pa]
f = f/101325;			#[atm]

# Results
print " The molar volume is %f m**3)/mol"%(V);
print " The fugacity is %f atm"%(f);

 The molar volume is 0.030469 m**3)/mol
The fugacity is 0.995168 atm


### Example 10.20 Page Number : 359¶

In :

P_1 = 6;			#[bar] - Initial pressure
P_1 = P_1*10**(5);			#[Pa]
T_1 = 100 + 273.15;			#[T] - Initial temperature
P_2 = 12;			#[bar] - Final pressure
P_2 = P_2*10**(5);			#[Pa]
T_2 = 500 + 273.15;			#[K] - Final temperature
R = 8.314;			#[J/mol*K] - Universal gas constant
Y = 1.126;			# Index of expansion
Cp_0 = (R*Y)/(Y-1);			#[J/mol-K]

# For propane
Tc = 369.8;			#[K]
Pc = 42.48;			#[bar]
Pc = Pc*10**(5);
w = 0.152;

# Calculations and Results
#(1)
# For van der Walls equation of state
a = (27*R**(2)*Tc**(2))/(64*Pc);			#[Pa-m**(6)/mol**(2)]
b = (R*Tc)/(8*Pc);			#[m**(3)/mol]

# The cubic form of van der Walls equation of state is given by,
# V**(3) - (b + (R*T)/P)*V**(2) + (a/P)*V - (a*b)/P = 0

# At state 1 (100 C and 6 bar)
def f(V):
return V**(3)-(b+(R*T_1)/P_1)*V**(2)+(a/P_1)*V-(a*b)/P_1
V_1_1 = fsolve(f,1)
V_1_2 = fsolve(f,10)
V_1_3 = fsolve(f,100)

V_1 = V_1_1;			#[m**(3)/mol]
# Thus compressibility factor is
Z_1 = (P_1*V_1)/(R*T_1);			#compressibility factor

H_R_1 = (Z_1 - 1)*R*T_1 - (a/V_1);			# [J/mol]
S_R_1 = R*math.log((P_1*(V_1-b))/(R*T_1));			# [J/mol-K]

# At state 2 (500 C and 12 bar)
def f1(V):
return V**(3)-(b+(R*T_2)/P_2)*V**(2)+(a/P_2)*V-(a*b)/P_2
V_2_1 = fsolve(f1,1)
V_2_2 = fsolve(f1,10)
V_2_3 = fsolve(f1,100)
# The largest root is considered because of molar volume of vapour phase is to determined
V_2 = V_2_1;			#[m**(3)/mol]
# Thus compressibility factor is
Z_2 = (P_2*V_2)/(R*T_2);			#compressibility factor

H_R_2 = (Z_2 - 1)*R*T_2 - (a/V_2);			# [J/mol]
S_R_2 = R*math.log((P_2*(V_2-b))/(R*T_2));			# [J/mol-K]

# Ideal gas entropy change is given by
delta_S_ig = Cp_0*math.log(T_2/T_1) - R*math.log(P_2/P_1);			#[J/mol-K]
# Entropy change is given by
delta_S = delta_S_ig + (S_R_2 - S_R_1);			#[J/mol-k]

# Ideal gas enthalpy change is given by
delta_H_ig = Cp_0*(T_2 - T_1);			#[J/mol]
# Enthalpy change is given by
delta_H = delta_H_ig + (H_R_2 - H_R_1);			#[J/mol]

print "1).The change in enthalpy is %f J/mol"%(delta_H);
print "    The change in entropy is %f J/mol-K"%(delta_S);

#(2)
# Virial equation of state

# At state 1 (372.15 K, 6 bar) let us calculate B and dB/dT
Tr = T_1/Tc;			# Reduced temperature
B_0 = 0.083-(0.422/(Tr)**(1.6));
B_1 = 0.139-(0.172/(Tr)**(4.2));

#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B = ((B_0+(w*B_1))*(R*Tc))/Pc;			#[m**(3)/mol]
dB0_dT = 0.422*1.6*Tc**(1.6)*T_1**(-2.6);			# (dB_0/dT) at state 1
dB1_dT = 0.172*4.2*Tc**(4.2)*T_1**(-5.2);			# (dB_1/dT) at state 1
dB_dT = ((R*Tc)/Pc)*((dB0_dT) + w*(dB1_dT));			# (dB/dT) at state 1

H_R_1_2 = B*P_1 - P_1*T_1*dB_dT;			#[J/mol] - Residual enthalpy at state 1
S_R_1_2 = -P_1*(dB_dT);			#[J/mol-K] - Residual entropy at state 1

# At state 2 (773.15 K, 12 bar)
Tr_2 = T_2/Tc;			# Reduced temperature
B_0_2 = 0.083-(0.422/(Tr_2)**(1.6));
B_1_2 = 0.139-(0.172/(Tr_2)**(4.2));

#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B_2 = ((B_0_2+(w*B_1_2))*(R*Tc))/Pc;			#[m**(3)/mol]
dB0_dT_2 = 0.422*1.6*Tc**(1.6)*T_2**(-2.6);			# (dB_0/dT) at state 1
dB1_dT_2 = 0.172*4.2*Tc**(4.2)*T_2**(-5.2);			# (dB_1/dT) at state 1
dB_dT_2 = ((R*Tc)/Pc)*((dB0_dT_2) + w*(dB1_dT_2));			# (dB/dT) at state 1

H_R_2_2 = B_2*P_2 - P_2*T_2*dB_dT_2;			#[J/mol] - Residual enthalpy at state 1
S_R_2_2 = -P_2*(dB_dT_2);			#[J/mol-K] - Residual entropy at state 1

delta_H_2 = delta_H_ig + (H_R_2_2 - H_R_1_2);			#[J/mol]
delta_S_2 = delta_S_ig + (S_R_2_2 - S_R_1_2);			#[J/mol]

print "2).The change in enthalpy is %f J/mol"%(delta_H_2);
print "    The change in entropy is %f J/mol-K"%(delta_S_2);

1).The change in enthalpy is 29798.211799 J/mol
The change in entropy is 48.649067 J/mol-K
2).The change in enthalpy is 29992.807365 J/mol
The change in entropy is 49.021724 J/mol-K


### Example 10.21 Page Number : 362¶

In :

# Variables
P = 2.76*10**(6);			#[N/m**(2)] - Pressure
T = 310.93;			#[K] - Temperature
R = 8.314;			#[J/mol*K] - Universal gas constant

# For n-butane
Tc = 425.18;			#[K] - Critical temperature
Pc = 37.97;			#[bar] - Critical pressure
Pc = Pc*10**(5);			#[Pa]
w = 0.193;
den = 0.61;			#[g/cm**(3)]
mol_wt = 58;			#[g/mol] - Molecular weight of butane

# Calculations and Results
# math.log(P_sat) = 15.7374 - 2151.63/(T-36.24)
P_sat = math.exp(15.7374 - 2151.63/(T-36.24));			#[mm Hg]
P_sat = (P_sat/760)*101325;			#[N/m**(2)]

#(1)
# Let us determine the second virial coefficient at 310.93 K
Tr = T/Tc;			# Reduced temperature
B_0 = 0.083-(0.422/(Tr)**(1.6));
B_1 = 0.139-(0.172/(Tr)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B = ((B_0+(w*B_1))*(R*Tc))/Pc;			#[m**(3)/mol]

# Fugacity under saturated conditions is given by
# math.log(f_sat/P_sat) = (B*P_sat)/(R*T)
f_sat = P_sat*(math.exp((B*P_sat)/(R*T)));			#[N/m**(2)]

# The molar volume is given by
V_liq = (1/(den*1000))*(mol_wt/1000);			#[m**(3)/mol]

f = f_sat*math.exp(V_liq*(P-P_sat)/(R*T));

print " 1).The fugacity of n-butane is %e N/m**2)"%(f);

#(2)
# For van der Walls equation of state
a = (27*R**(2)*Tc**(2))/(64*Pc);			#[Pa-m**(6)/mol**(2)]
b = (R*Tc)/(8*Pc);			#[m**(3)/mol]

# The cubic form of van der Walls equation of state is given by,
# V**(3) - (b + (R*T)/P)*V**(2) + (a/P)*V - (a*b)/P = 0

# At 100 C and 1 atm
def f(V):
return V**(3)-(b+(R*T)/P)*V**(2)+(a/P)*V-(a*b)/P
V_1 = fsolve(f,0.1)
V_1 = fsolve(f,10)
V_1 = fsolve(f,100)
# The above equation has only 1 real root, other two roots are imaginary
V = V_1;			#[m**(3)/mol]

# math.log(f/P) = math.log((R*T)/(P*(V-b))) + b/(V-b) -(2*a)/(R*T*V)
f_2 = P*(math.exp(math.log((R*T)/(P*(V-b))) + b/(V-b) -(2*a)/(R*T*V)));

print " 2).The fugacity of n-butane is %e N/m**2)"%(f_2);

 1).The fugacity of n-butane is 3.293476e+05 N/m**2)
2).The fugacity of n-butane is 8.997352e+05 N/m**2)


### Example 10.22 Page Number : 363¶

In :

# Variables
T = 50+273.15;			#[K] - Temperature
P = 25.*10**(3);			#[Pa] - Pressure
y1 = 0.5;			#[mol] - mole fraction of equimolar mixture
y2 = 0.5;
R = 8.314;			#[J/mol*K] - Universal gas constant

#For component 1 (methyl ethyl ketone)
Tc_1  = 535.5;			#[K] - Critical temperature
Pc_1 = 41.5*10**(5);			#[N/m**(2)] - Critical pressure
Vc_1 = 267.;			#[cm**(3)/mol] - Critical volume
Zc_1 = 0.249;			# Critical compressibility factor
w_1 = 0.323;			# acentric factor

#For component 2 (toluene)
Tc_2 = 591.8;			#[K]
Pc_2 = 41.06*10**(5);			#[N/m**(2)]
Vc_2 = 316.;			#[cm**(3)/mol]
Zc_2 = 0.264;
w_2 = 0.262;

# Calculations and Results
# For equation of state Z = 1 + B/V
#For component 1, let us calculate B and dB/dT
Tr_1 = T/Tc_1;			#Reduced temperature
#At reduced temperature
B1_0 = 0.083-(0.422/(Tr_1)**(1.6));
B1_1 = 0.139-(0.172/(Tr_1)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0+(w*B_1)
B_11 = ((B1_0+(w_1*B1_1))*(R*Tc_1))/Pc_1;			# [m**(3)/mol-K]
dB0_dT_1 = 0.422*1.6*Tc_1**(1.6)*T**(-2.6);			# [m**(3)/mol-K] -  (dB_0/dT)
dB1_dT_1 = 0.172*4.2*Tc_1**(4.2)*T**(-5.2);			# [m**(3)/mol-K] -  (dB_1/dT)
dB_dT_1 = ((R*Tc_1)/Pc_1)*((dB0_dT_1) + w_1*(dB1_dT_1));			#[m**(3)/mol-K] -  (dB/dT)_

#Similarly for component 2
Tr_2 = T/Tc_2;			#Reduced temperature
#At reduced temperature Tr_2,
B2_0 = 0.083 - (0.422/(Tr_2)**(1.6));
B2_1 = 0.139 - (0.172/(Tr_2)**(4.2));
B_22 = ((B2_0+(w_2*B2_1))*(R*Tc_2))/Pc_2;			#[m**(3)/mol]
dB0_dT_2 = 0.422*1.6*Tc_2**(1.6)*T**(-2.6);			# [m**(3)/mol-K] -  (dB_0/dT)
dB1_dT_2 = 0.172*4.2*Tc_2**(4.2)*T**(-5.2);			# [m**(3)/mol-K] -  (dB_1/dT)
dB_dT_2 = ((R*Tc_2)/Pc_2)*((dB0_dT_2) + w_2*(dB1_dT_2));			#[m**(3)/mol-K] -  (dB/dT)_

#For cross coeffcient, let us calculate B and dB/dT
Tc_12 = (Tc_1*Tc_2)**(1./2);			#[K]
w_12 = (w_1 + w_2)/2.;
Zc_12 = (Zc_1 + Zc_2)/2;
Vc_12 = (((Vc_1)**(1./3) + (Vc_2)**(1./3))/2)**(3);			#[cm**(3)/mol]
Vc_12 = Vc_12*10**(-6);			#[m**(3)/mol]
Pc_12 = (Zc_12*R*Tc_12)/Vc_12;			#[N/m**(2)]

#Now we have,(B_12*Pc_12)/(R*Tc_12) = B_0+(w_12*B_1)
#where B_0 and B_1 are to be evaluated at Tr_12
Tr_12 = T/Tc_12;
#At reduced temperature Tr_12
B_0 = 0.083 - (0.422/(Tr_12)**(1.6));
B_1 = 0.139 - (0.172/(Tr_12)**(4.2));
B_12 = ((B_0 + (w_12*B_1))*R*Tc_12)/Pc_12;			#[m**(3)/mol]
dB0_dT_12 = 0.422*1.6*Tc_12**(1.6)*T**(-2.6);			# [m**(3)/mol-K] -  (dB_0/dT)
dB1_dT_12 = 0.172*4.2*Tc_12**(4.2)*T**(-5.2);			# [m**(3)/mol-K] -  (dB_1/dT)
dB_dT_12 = ((R*Tc_12)/Pc_12)*((dB0_dT_12) + w_12*(dB1_dT_12));			#[m**(3)/mol-K] -  (dB/dT)_12

#For the mixture
B = y1**(2)*B_11 + 2*y1*y2*B_12 + y2**(2)*B_22;			#[m**(3)/moL]

# The equation of state can be written as
# V**(2) - ((R*T)/P) - (B*R*T)/P = 0
# V**(2) - 0.1075*V + 1.737*10**(-4) = 0
def f(V):
return V**(2) - 0.1075*V + 1.737*10**(-4)
V1 = fsolve(f,0.1)
V2 = fsolve(f,1)
# We will consider the root which is near to R*T/P
V = V1;
# dB/dT = y_1**(2)*dB_11/dT + y_2**(2)*dB_22/dT + 2*y_1*y_2*dB_12/dT
dB_dT = y1**(2)*dB_dT_1 + y2**(2)*dB_dT_2 + 2*y1*y2*dB_dT_12;			#[m**(3)/mol-K]

# For equation of state Z = 1 + B/V
H_R = (B*R*T)/V - ((R*T**(2))/V)*dB_dT;			#[J/mol]

print " 1).The value of H_R for the mixture using virial equation of state is %f J/mol"%(H_R);

#(2)
#  For van der Walls equation of state
a_11 = (27*R**(2)*Tc_1**(2))/(64*Pc_1);			#[Pa-m**(6)/mol**(2)]
a_22 = (27*R**(2)*Tc_2**(2))/(64*Pc_2);			#[Pa-m**(6)/mol**(2)]
a_12 = (a_11*a_22)**(1./2);
b_1 = (R*Tc_1)/(8*Pc_1);			#[m**(3)/mol]
b_2 = (R*Tc_2)/(8*Pc_2);			#[m**(3)/mol]

# For the mixture
a = y1**(2)*a_11 + y2**(2)*a_22 + 2*y1*y2*a_12;			#[Pa-m**(6)/mol**(2)]
b = y1*b_1 + y2*b_2;			#[m**(3)/mol]

# From the cubic form of van der Walls equation of state
def f1(V):
return V**(3)-(b+(R*T)/P)*V**(2)+(a/P)*V-(a*b)/P
V2_1 = fsolve(f1,0.1)
V2_2 = fsolve(f1,10)
V2_3 = fsolve(f1,100)
# The largest root is considered
V_2 = V2_1;

# The residual enthalpy is given by
H_R_2 = P*V_2 - R*T -a/V_2;			#[J/mol]

print " 2).The value of H_R for the mixture using van der Walls equation of state is %f J/mol"%(H_R_2);

 1).The value of H_R for the mixture using virial equation of state is -151.289795 J/mol
2).The value of H_R for the mixture using van der Walls equation of state is -38.476127 J/mol


### Example 10.23 Page Number : 366¶

In :
from scipy.optimize import fsolve
import math

# Variables
T = 320 + 273.15;			#[K]
R = 8.314;			#[J/mol*K] - Universal gas constant

# For water
Tc = 647.1;			#[K]
Pc = 220.55;			#[bar]
Pc = Pc*10**(5);			#[Pa]

# The cubic form of Redlich Kwong equation of state is given by,
# V**(3) - ((R*T)/P)*V**(2) - ((b_1**(2)) + ((b_1*R*T)/P) - (a/(T**(1/2)*P))*V - (a*b)/(T**(1/2)*P) = 0

# At 320 C and 70 bar pressure
P_1 = 70;			#[bar]
P_1 = P_1*10**(5);			#[Pa]

# Calculations and Results
a = (0.42748*(R**(2))*(Tc**(2.5)))/Pc;			#[Pa*m**(6)*K**(1/2)/mol]
b = (0.08664*R*Tc)/Pc;			#[m**(3)/mol]
# Solving the cubic equation
def f1(V):
return V**(3)-((R*T)/P_1)*V**(2)-((b**(2))+((b*R*T)/P_1)-(a/(T**(1./2)*P_1)))*V-(a*b)/(T**(1./2)*P_1)
V1=fsolve(f1,1)
V2=fsolve(f1,10)
V3=fsolve(f1,100)
# The largest root is considered because at 320 C and 70 bar vapour phase exists.
V_1 = V1;			#[m**(3)/mol]
# Thus compressibility factor is
Z_1 = (P_1*V_1)/(R*T);

# For Redlich-Kwong equation of state
# math.log(f/P) = Z - 1 - math.log(V_1/(V_1-b)) + (a/(b*R*(T**(3/2))))*math.log(V/(V+b))
f_1 = P_1*(math.exp(Z_1-1-math.log(Z_1)+math.log(V_1/(V_1-b))+(a/(b*R*(T**(3./2))))*math.log(V_1/(V_1+b))));			#[Pa]
f_1 = f_1*10**(-5);			#[bar]

print " The fugacity of water vapour at 320 C and 70 bar pressure is %f bar"%(f_1);

# At 320 C and 170 bar pressure, we have
P_2 = 170;			#[bar]
P_2 = P_2*10**(5);			#[Pa]

# Solving the cubic equation
def f2(V):
return V**(3)-((R*T)/P_2)*V**(2)-((b**(2))+((b*R*T)/P_2)-(a/(T**(1./2)*P_2)))*V-(a*b)/(T**(1./2)*P_2)
V4 = fsolve(f2,1)
V5 = fsolve(f2,10)
V6 = fsolve(f2,100)
# The above equation has only 1 real root,other two roots are imaginary. Therefore,
V_2 = V6;			#[m**(3)/mol]
# Thus compressibility factor is
Z_2 = (P_2*V_2)/(R*T);

# For Redlich-Kwong equation of state
# math.log(f/P) = Z - 1 - math.log(V_1/(V_1-b)) + (a/(b*R*(T**(3/2))))*math.log(V/(V+b))
f_2 = P_2*(math.exp(Z_2-1-math.log(Z_2)+math.log(V_2/(V_2-b))+(a/(b*R*(T**(3./2))))*math.log(V_2/(V_2+b))));			#[Pa]
f_2 = f_2*10**(-5);			#[bar]

print " The fugacity of water vapour at 320 C and 170 bar pressure is %f bar"%(f_2);

 The fugacity of water vapour at 320 C and 70 bar pressure is 60.456239 bar
The fugacity of water vapour at 320 C and 170 bar pressure is 101.590989 bar

C:\Anaconda\lib\site-packages\scipy\optimize\minpack.py:227: RuntimeWarning: The iteration is not making good progress, as measured by the
improvement from the last ten iterations.
warnings.warn(msg, RuntimeWarning)


### Example 10.24 Page Number : 367¶

In :
from scipy.optimize import fsolve
import math
# Variables
Vol = 0.057;			#[m**(3)] - Volume of car tyre
P_1 = 300.;			#[kPa] - Initial pressure
P_1 = P_1*10**(3);			#[Pa]
T_1 = 300.;			#[K] - Initial temperature
P_2 = 330.;			#[kPa] - Finnal pressure
P_2 = P_2*10**(3);			#[Pa]
R = 8.314;			#[J/mol*K] - Universal gas constant
Cv_0 = 21.;			#[J/mol-K] - Heat capacity for air

# For oxygen
Tc_O2 = 154.6;			#[K] - Critical temperature
Pc_O2 = 50.43;			#[bar] - Critical pressure
Pc_O2 = Pc_O2*10**(5);			#[Pa]
y1 = 0.21;			# - Mole fraction of oxygen
# For nitrogen
Tc_N2 = 126.2;			#[K] - Critical temperature
Pc_N2 = 34.00;			#[bar] - Critical pressure
Pc_N2 = Pc_N2*10**(5);			#[Pa]
y2 = 0.79;			# - Mole fraction of nitrogen

# Calculations and Results
# (1)
# Assuming ideal gas behaviour. The volume remains the same,therefore,we get
# P_1/T_1 = P_2/T_2
T_2 = P_2*(T_1/P_1);			#[K]

n = (P_1*Vol)/(R*T_1);			#[mol] - Number of moles
delta_U = n*Cv_0*(T_2-T_1);			#[J]

print " 1).The change in internal energy for ideal gas behaviour) is %f J"%(delta_U);

#(2)
#  For van der Walls equation of state
a_O2 = (27*R**(2)*Tc_O2**(2))/(64*Pc_O2);			#[Pa-m**(6)/mol**(2)]
a_N2 = (27*R**(2)*Tc_N2**(2))/(64*Pc_N2);			#[Pa-m**(6)/mol**(2)]
a_12 = (a_O2*a_N2)**(1./2);
b_O2 = (R*Tc_O2)/(8*Pc_O2);			#[m**(3)/mol]
b_N2 = (R*Tc_N2)/(8*Pc_N2);			#[m**(3)/mol]

# For the mixture
a = y1**(2)*a_O2 + y2**(2)*a_N2 + 2*y1*y2*a_12;			#[Pa-m**(6)/mol**(2)]
b = y1*b_O2 + y2*b_N2;			#[m**(3)/mol]

# From the cubic form of van der Walls equation of state
# At 300 K and 300 kPa,
def f1(V):
return V**(3)-(b+(R*T_1)/P_1)*V**(2)+(a/P_1)*V-(a*b)/P_1
V_1 = fsolve(f1,0.1)
V_2 = fsolve(f1,10)
V_3 = fsolve(f1,100)
# The above equation has only 1 real root, other two roots are imaginary
V = V_1;			#[m**(3)/mol]

# Now at P = 330 kPa and at molar volume V
# The van der Walls equation of state is
# (P + a/V**(2))*(V - b) = R*T
T_2_prime = ((P_2 + a/V**(2))*(V - b))/R;			#[K]
n_prime = Vol/V;			#[mol]

# Total change in internal energy is given by
# delta_U_prime = n_prime*delta_U_ig + n_prime*(U_R_2 - U_R_1)
# delta_U_prime = n_prime*Cv_0*(T_2_prime - T_1) + n_prime*a(1/V_2 - 1/V_1);
# Since V_1 = V_2 = V, therefore
delta_U_prime = n_prime*Cv_0*(T_2_prime - T_1);

print " 2).The change in internal energy for van der Walls equation of state) is %f J"%(delta_U_prime);

 1).The change in internal energy for ideal gas behaviour) is 4319.220592 J
2).The change in internal energy for van der Walls equation of state) is 4299.872282 J


### Example 10.25 Page Number : 369¶

In :
from scipy.optimize import fsolve
import math
# Variables
T_1 = 150 + 273.15;			#[K] - Initial emperature
T_2 = T_1;			# Isothermal process
P_1 = 100.*10**(3);			#[Pa] - Initial pressure
P_2 = 450.*10**(3);			#[Pa] - Final pressure
R = 8.314;			#[J/mol*K] - Universal gas constant
# For water
Tc = 647.1;			#[K] - Critical temperature
Pc = 220.55;			#[bar] - Critical pressure
Pc = Pc*10.**(5);
w = 0.345;
Mol_wt = 18.015;			#[g/mol] - Molecular weight of water
Cp_0 = 4.18;			#[J/mol-K] - Smath.tan(math.radiansard heat capacity of water

# Calculations

# In Peng-Robinson equation of state
m = 0.37464 + 1.54226*w - 0.26992*w**(2);
# At T_1 and P_1, we have
Tr = T_1/Tc;
alpha = (1 + m*(1 - Tr**(1./2)))**(2);
a = ((0.45724*(R*Tc)**(2))/Pc)*alpha;			#[Pa*m**(6)/mol**(2)]
b = (0.07780*R*Tc)/Pc;			#[m**(3)/mol]

# Cubuc form of Peng-Robinson equation of stste is given by
# V**(3)+(b-(R*T)/P)*V**(2)-((3*b**(2))+((2*R*T*b)/P)-(a/P))*V+b**(3)+((R*T*(b**(2))/P)-((a*b)/P)=0;
# Solving the cubic equation
def f(V):
return V**(3)+(b-(R*T_1)/P_1)*V**(2)-((3*b**(2))+((2*R*T_1*b)/P_1)-(a/P_1))*V+b**(3)+((R*T_1*(b**(2)))/P_1)-((a*b)/P_1)
V1 = fsolve(f,-1)
V2 = fsolve(f,0)
V3 = fsolve(f,1)
#The largest root is for vapour phase,
#The largest root is only considered as the systemis gas
V_1 = V3;			#[m**(3)/mol]
# Thus compressibility factor is
Z_1 = (P_1*V_1)/(R*T_1);

# At T_2 and P_2, we have
# Cubuc form of Peng-Robinson equation of stste is given by
# V**(3)+(b-(R*T)/P)*V**(2)-((3*b**(2))+((2*R*T*b)/P)-(a/P))*V+b**(3)+((R*T*(b**(2))/P)-((a*b)/P)=0;
# Solving the cubic equation
def f1(V):
return V**(3)+(b-(R*T_2)/P_2)*V**(2)-((3*b**(2))+((2*R*T_2*b)/P_2)-(a/P_2))*V+b**(3)+((R*T_2*(b**(2)))/P_2)-((a*b)/P_2)
V4 = fsolve(f1,-1)
V5 = fsolve(f1,0)
V6 = fsolve(f1,1)

V_2 = V6;			#[m**(3)/mol]
# Thus compressibility factor is
Z_2 = (P_2*V_2)/(R*T_2);

# In the Peng-Robinson equation of stste
# da/dT = -(a*m)/((alpha*T*Tc)**(1/2))
# The residual enthalpy is given by
# H_R = R*T*(Z-1) + (((T*(da_dT))-a)/(2*2**(1/2)*b))*math.log((Z+(1+2**(1/2)*((P*b)/(R*T))))/(Z+(1-2**(1/2)*((P*b)/(R*T)))))

# At state 1
da_dT_1 = -(a*m)/((alpha*T_1*Tc)**(1./2));			#[Pa*m**(6)/mol**(2)]
H_R_1 = R*T_1*(Z_1-1) + (((T_1*(da_dT_1))-a)/(2*(2**(1./2))*b))* \
math.log((Z_1+(1+2**(1./2))*((P_1*b)/(R*T_1)))/(Z_1+(1-2**(1./2))*((P_1*b)/(R*T_1))));

# At state 2
da_dT_2 = -(a*m)/((alpha*T_2*Tc)**(1/2.));			#[Pa*m**(6)/mol**(2)]
H_R_2 = R*T_2*(Z_2-1) + (((T_2*(da_dT_2))-a)/(2*2**(1./2)* \
b))*math.log((Z_2+(1+2**(1./2))*((P_2*b)/(R*T_1)))/(Z_2+(1-2**(1./2))*((P_2*b)/(R*T_1))));

delta_H = H_R_2 - H_R_1;			#[J/mol]
delta_H = delta_H/Mol_wt;			#[kJ/kg]

# The residual entropy relation for a substance following Peng - Robinson equation of state ia
# S_R = R*math.log(Z - (P*b)/(R*T)) + (da_dT/(2*2**(1/2)*b))*math.log((Z+(1+2**(1/2))*((P*b)/(R*T)))/(Z+(1-2**(1/2))*((P*b)/(R*T))))

# The residual entropy at state 1 is
S_R_1 = R*math.log(Z_1 - (P_1*b)/(R*T_1)) + (da_dT_1/(2*2**(1./2)*b))* \
math.log((Z_1+(1+2**(1./2))*((P_1*b)/(R*T_1)))/(Z_1+(1-2**(1./2))*((P_1*b)/(R*T_1))));

# The residual entropy at state 2 is
S_R_2 = R*math.log(Z_2 - (P_2*b)/(R*T_2)) + (da_dT_2/(2*2**(1./2)*b)) \
*math.log((Z_2+(1+2**(1./2))*((P_2*b)/(R*T_2)))/(Z_2+(1-2**(1./2))*((P_2*b)/(R*T_2))));

delta_S_R = S_R_2 - S_R_1;			#[J/mol-K]

# The ideal gas change in entropy is
delta_S_ig = Cp_0*math.log(T_2/T_1) - R*math.log(P_2/P_1);			#[J/mol-K]

# Therefore
delta_S = delta_S_R + delta_S_ig;			#[J/mol-K]

# Results
print " The enthalpy change is given by delta_H = %f kJ/mol"%(delta_H);
print " The entropy change is given by  delta_S = %f J/mol-K"%(delta_S);

 The enthalpy change is given by delta_H = -11.747722 kJ/mol
The entropy change is given by  delta_S = -12.826050 J/mol-K


### Example 10.26 Page Number : 370¶

In :
from scipy.optimize import fsolve
import math
# Variables
Vol = 0.15;			#[m**(3)]
T_1 = 170;			#[K] - Initial emperature
P_1 = 100;			#[bar] - Initial pressure
P_1 = P_1*10**(5);			#[Pa]
R = 8.314;			#[J/mol*K] - Universal gas constant
# For nitrogen
Tc = 126.2;			#[K] - Critical tempeature
Pc = 34;			#[bar] - Critical pressure
Pc = Pc*10**(5);			#[Pa]
w = 0.038;
# Cp_0 = 27.2+4.2*10**(-3)*T

# Calculations and Results
#(1)
# For van der Walls equation of state
a = (27*R**(2)*Tc**(2))/(64*Pc);			#[Pa-m**(6)/mol**(2)]
b = (R*Tc)/(8*Pc);			#[m**(3)/mol]

# The cubic form of van der Walls equation of state is given by,
# V**(3) - (b + (R*T)/P)*V**(2) + (a/P)*V - (a*b)/P = 0
# On simplification the equation changes to
# V**(3) - 1.799*10**(4)*V**(2) + 1.366*10**(-8)*V - 5.269*10**(-13) = 0

# Solving the cubic equation
def f(V):
return V**(3)-1.799*10**(-4)*V**(2) + 1.366*10**(-8)*V - 5.269*10**(-13)
V1 = fsolve(f,1)
V2 = fsolve(f,10)
V3 = fsolve(f,100)
# The above equation has only 1 real root, other two roots are imagimnary
V_1 = V1;			#[m**(3)/mol]
# Thus total number of moles is given by
n_1 = Vol/V_1;			#[mol]

# After 500 mol are withdrawn, the final number of moles is given by
n_2 = n_1 - 500;			#[mol]
# Thus molar volume at final state is
V_2 = Vol/n_2;			#[m**(3)/mol]

# The ideal entropy change is guven by

def f24(T):
return 27.2+4.2*10**(-3)*T

# delta_S_ig =  quad(f24,T_1,T_2) - R*math.log(P_2/P_1)

# The residual entropy change is given by
# delta_S_R = R*math.log((P_2*(V_2-b))/(R*T_2)) - R*math.log((P_1*(V_1-b))/(R*T_1))
# delta_S = delta_S_ig = delta_S_R

def f25(T):
return 27.2+4.2*10**(-3)*T

# delta_S =  quad(f25,T_1,T_2) + R*math.log((V_2-b)/(V_1-b))

# During discharging delta_S = 0, thus on simplification we get
# 18.886*math.log(T_2) + 4.2*10**(-3)*T_2 - 92.937 = 0
# Solving the above equation we get
def f1(T_2):
return 18.886*math.log(T_2) + 4.2*10**(-3)*T_2 - 92.937
T_2 = fsolve(f1,1)

# Thus at T_2,
P_2 = (R*T_2)/(V_2-b) - a/V_2**(2);			#[N/m**(2)]
P_2 = P_2*10**(-5);			#[bar]

print " 1).The final temperature is %f K"%(T_2);
print "     The final pressure is %f bar"%(P_2);

#(2)
# In Peng-Robinson equation of state
m = 0.37464 + 1.54226*w - 0.26992*w**(2);
# At T_1 and P_1, we have
Tr = T_1/Tc;
alpha = (1 + m*(1 - Tr**(1./2)))**(2);
a_2 = ((0.45724*(R*Tc)**(2))/Pc)*alpha;			#[Pa*m**(6)/mol**(2)]
b_2 = (0.07780*R*Tc)/Pc;			#[m**(3)/mol]

# Cubuc form of Peng-Robinson equation of stste is given by
# V**(3)+(b-(R*T)/P)*V**(2)-((3*b**(2))+((2*R*T*b)/P)-(a/P))*V+b**(3)+((R*T*(b**(2))/P)-((a*b)/P)=0;
# Solving the cubic equation
def f2(V):
return V**(3)+(b_2-(R*T_1)/P_1)*V**(2)-((3*b_2**(2))+((2*R*T_1*b_2)/P_1)-(a_2/P_1))*V+b_2**(3)+((R*T_1*(b_2**(2)))/P_1)-((a_2*b_2)/P_1)
V4 = fsolve(f2,-1)
V5 = fsolve(f2,0)
V6 = fsolve(f2,0.006)
#The above equation has only 1 real root,the other two roots are imaginary
V_1_2 = V6;			#[m**(3)/mol]

# The number of moles in the initial state is given by
n_1_2 = Vol/V_1_2;			#[mol]
# After 500 mol are withdrawn, the final number of moles is given by
n_2_2 = n_1_2 - 500;			#[mol]
# Thus molar volume at final state is
V_2_2 = Vol/n_2_2;			#[m**(3)/mol]

# At the final state the relation between pressure and temperature is
# P_2_2 = (R*T_2_2)/(V_2_2-b_2) - a_2/V_2_2**(2)
# P_2_2 = 7.23*10**(4)*T_2 - 3.93*10**(7)*a_2

# Now let us calculate the residual entropy at initial state
Z_1 = (P_1*V_1_2)/(R*T_1);
da_dT_1 = -(a*m)/((alpha*T_1*Tc)**(1./2));			#[Pa*m**(6)/mol**(2)] - da/dT

# The residual entropy change for Peng-Robinson equatiob of state is given by
# S_R = R*math.log(Z-(P*b)/(R*T)) + (da_dT/(2*2**(1/2)*b))*math.log((V+(1+2**(1/2))*b))/((V+(1-2**(1/2)*b))));
S_R_1 = R*(math.log(Z_1-(P_1*b_2)/(R*T_1))) + (da_dT_1/(2*2**(1.2)*b_2))*(math.log((V_1_2+(1+2**(1./2))*b_2)/(V_1_2+(1-2**(1./2))*b_2)));

# The total entropy change is given by
# delta_S = delta_S_ig + delta_S_R

def f26(T):
return 27.2+4.2*10**(-3)*T

# where, delta_S_ig =  quad(f26,T_1,T_2_2) - R*math.log(P_2_2/P_1)

# and, P_2_2 = (R*T_2_2)/(V_2_2-b_2) - a_2/V_2_2**(2)
# On simplification we get
# delta_S = 27.2*math.log(T_2_2-T_1) + 4.2*10**(-3)*(T_2_2-T_1) - R*math.log(P_2_2/P_1) + R*math.log(Z_2-(P_2_2*b)/(R*T_2_2)) + 6226*(da_dT_2) + 9.22

# Now we have the determine the value of T_2_2 such that delta_S = 0
# Starting with a temperature of 150 K
T_prime = 100.;			#[K]
error = 10.;
while(error>0.1):
Tr_prime = T_prime/Tc;
alpha_prime = (1 + m*(1 - Tr_prime**(1./2)))**(2);
a_prime = ((0.45724*(R*Tc)**(2))/Pc)*alpha_prime;
P_prime = 7.23*10**(4)*T_prime - 3.93*10**(7)*a_prime;
Z_prime = (P_prime*V_2_2)/(R*T_prime);
da_dT_prime = -(a_prime*m)/((alpha_prime*T_prime*Tc)**(1./2));
delta_S = 27.2*math.log(T_prime/T_1) + 4.2*10**(-3)*(T_prime-T_1) - R*math.log(P_prime/P_1) + R*math.log(Z_prime-((P_prime*b_2)/(R*T_prime))) + 6226*(da_dT_prime) + 9.22;
error=abs(delta_S);
T_prime = T_prime + 0.3;

T_2_2 = T_prime;			#[K] - Final temperature
P_2_2 = P_prime*10**(-5);			#[bar] - Final pressure

print " 2).The final temperature is %f K"%(T_2_2);
print "     The final pressure is %f bar"%(P_2_2);

 1).The final temperature is 133.131844 K
The final pressure is 39.636659 bar
2).The final temperature is 134.200000 K
The final pressure is 40.130674 bar


### Example 10.27 Page Number : 374¶

In :
from scipy.optimize import fsolve
import math
# Variables
T = 373.15;			#[K]
Tc = 562.16;			#[K]
Pc = 48.98;			#[bar]
Pc = Pc*10**(5);			#[Pa]
R = 8.314;			#[J/mol-K] - Universal gas constant

# The cubic form of Redlich Kwong equation of state is given by,
# V**(3) - ((R*T)/P)*V**(2) - ((b_1**(2)) + ((b_1*R*T)/P) - (a/(T**(1/2)*P))*V - (a*b)/(T**(1/2)*P) = 0

# Calculations
a = (0.42748*(R**(2))*(Tc**(2.5)))/Pc;			#[Pa*m**(6)*K**(1/2)/mol]
b = (0.08664*R*Tc)/Pc;			#[m**(3)/mol]

# At 373.15 K, let us assume the pressure to be 2.5 bar and under these conditions
P_1 = 2.5;			#[bar]
P_1 = P_1*10**(5);			#[bar]

# Putting the values in Redlich Kwong equation of state, the equation becomes
# V**(3) - 0.0124*V**(2) + 8.326*10**(-6)*V - 7.74*10**(-10) = 0
# Solving the cubic equation

def f(V):
return V**(3) - 0.0124*V**(2) + 8.326*10**(-6)*V - 7.74*10**(-10)
V1=fsolve(f,-9)
V2=fsolve(f,10)
V3=fsolve(f,0.1)
# The largest root and the smallest root is considered for liquid phase and vapour phase respectively.
V_liq = V1;			#[m**(3)/mol] - Molar volume in liquid phase
V_vap = V3;			#[m**(3)/mol] - Molar volume in vapour phase

# Let us calculate the fugacity of vapour phase
# math.log(f_vap/P) = b/(V-b) + math.log((R*T)/(P*(V-b))) - (a/(R*T**(1.5)))*(1/(V+b) - (1/b)*math.log(V/(V+b)))
f_vap = P_1*math.exp(b/(V_vap-b) + math.log((R*T)/(P_1*(V_vap-b))) - (a/(R*T**(1.5)))*(1/(V_vap+b) - (1/b)*math.log(V_vap/(V_vap+b))));			#[Pa]

# Let us calculate the fugacity of the liquid phase
f_liq = P_1*math.exp(b/(V_liq-b) + math.log((R*T)/(P_1*(V_liq-b))) - (a/(R*T**(1.5)))*(1/(V_liq+b) - (1/b)*math.log(V_liq/(V_liq+b))));

# The two fugacities are not same; therefore another pressure is to be assumed. The new pressure is
P_new = P_1*(f_liq/f_vap);			#[Pa]

# At P_new
def f1(V):
return V**(3) - ((R*T)/P_new)*V**(2) - (b**(2) + ((b*R*T)/P_new) - a/(T**(1/2)*P_new))*V - (a*b)/(T**(1/2)*P_new)
V4=fsolve(f1,-9)
V5=fsolve(f1,10)
V6=fsolve(f1,0.1)
# The largest root and the smallest root is considered for liquid phase and vapour phase respectively.
V_liq_2 = V4;			#[m**(3)/mol] - Molar volume in liquid phase
V_vap_2 = V6;			#[m**(3)/mol] - Molar volume in vapour phase

f_vap_prime = P_new*math.exp(b/(V_vap_2-b) + math.log((R*T)/(P_new*(V_vap_2-b))) - (a/(R*T**(1.5)))*(1/(V_vap_2+b) - (1/b)*math.log(V_vap_2/(V_vap_2+b))));			#[Pa]
f_liq_prime = P_new*math.exp(b/(V_liq_2-b) + math.log((R*T)/(P_new*(V_liq_2-b))) - (a/(R*T**(1.5)))*(1/(V_liq_2+b) - (1/b)*math.log(V_liq_2/(V_liq_2+b))));

# Since the fugacities of liquid and vapour phasesare almost same the assumed pressure may be taken as vapour pressure at 373.15 K
P_new = P_new*10**(-5);			#[bar]

# Results
print " The vapour pressure of benzene using Redlich Kwong equation of state is %f bar"%(P_new);

 The vapour pressure of benzene using Redlich Kwong equation of state is 2.666995 bar


### Example 10.28 Page Number : 374¶

In :
from scipy.optimize import fsolve
import math
# Variables
T = 150 + 273.15;			#[K]
Tc = 647.1;			#[K]
Pc = 220.55;			#[bar]
Pc = Pc*10**(5);			#[Pa]
w = 0.345;
R = 8.314;			#[J/mol-K] - Universal gas constant

# Let us assume a pressure of 100 kPa.
P_1 = 100.*10**(3);			#[Pa]

# Calculations
# At 100 kPa and 423.15 K, from Peng-Robinson equation of stste
m = 0.37464 + 1.54226*w - 0.26992*w**(2);
Tr = T/Tc;
alpha = (1 + m*(1 - Tr**(1./2)))**(2);
a = ((0.45724*(R*Tc)**(2))/Pc)*alpha;			#[Pa*m**(6)/mol**(2)]
b = (0.07780*R*Tc)/Pc;			#[m**(3)/mol]
# Cubic form of Peng-Robinson equation of stste is given by
# V**(3)+(b-(R*T)/P)*V**(2)-((3*b**(2))+((2*R*T*b)/P)-(a/P))*V+b**(3)+((R*T*(b**(2))/P)-((a*b)/P)=0;
# Solving the cubic equation
def f(V):
return V**(3)+(b-(R*T)/P_1)*V**(2)-((3*b**(2))+((2*R*T*b)/P_1)-(a/P_1))*V+b**(3)+((R*T*(b**(2)))/P_1)-((a*b)/P_1)
V1 = fsolve(f,-1)
V2 = fsolve(f,0)
V3 = fsolve(f,1)
# The largest root and the smallest root is considered for liquid phase and vapour phase respectively.
V_liq = V1;			#[m**(3)/mol] - Molar volume in liquid phase
V_vap = V3;			#[m**(3)/mol] - Molar volume in vapour phase

# The compressibility factor is given by
Z_vap = (P_1*V_vap)/(R*T);			# For liquid phase
Z_liq = (P_1*V_liq)/(R*T);			# For vapour phase

# The math.expression for fugacity of Peng Robinson equation is
# math.log(f/P) = (Z-1) - math.log(Z-((P*b)/(R*T))) - (a/(2*2**(1/2)*b*R*T))*math.log((Z+(1+2**(1/2))*((P*b)/(R*T)))/((Z+(1-2**(1/2))*((P*b)/(R*T)))
# For vapour phase
f_P_vap = math.exp((Z_vap-1) - math.log(Z_vap-((P_1*b)/(R*T))) - (a/(2*2**(1./2)*b*R*T))*math.log((Z_vap+(1+2**(1./2))*((P_1*b)/(R*T)))/(Z_vap+(1-2**(1./2))*((P_1*b)/(R*T)))));
# For liquid phase
f_P_liq = math.exp((Z_liq-1) - math.log(Z_liq-((P_1*b)/(R*T))) - (a/(2*2**(1./2)*b*R*T))*math.log((Z_liq+(1+2**(1./2))*((P_1*b)/(R*T)))/(Z_liq+(1-2**(1./2))*((P_1*b)/(R*T)))));

# Therefore f_liq/f_vap can be calculated as
fL_fV = (f_P_liq/f_P_vap);

# The two values (f/P)_vap and (f/P)_vap are not same [ (f_P_liq/f_P_vap) >1 ]; therefore another pressure is to be assumed. The new pressure be
P_new = P_1*(f_P_liq/f_P_vap);			#[Pa]

# At P_new and 423.15 K, from Peng-Robinson equation of stste

# V**(3)+(b-(R*T)/P)*V**(2)-((3*b**(2))+((2*R*T*b)/P)-(a/P))*V+b**(3)+((R*T*(b**(2))/P)-((a*b)/P)=0;
# Solving the cubic equation
def f(V):
return V**(3)+(b-(R*T)/P_new)*V**(2)-((3*b**(2))+((2*R*T*b)/P_new)-(a/P_new))*V+b**(3)+((R*T*(b**(2)))/P_new)-((a*b)/P_new)
V4 = fsolve(f,-1)
V5 = fsolve(f,0)
V6 = fsolve(f,1)
# The largest root and the smallest root is considered for liquid phase and vapour phase respectively.
V_liq_2 = V4;			#[m**(3)/mol] - Molar volume in liquid phase
V_vap_2 = V6;			#[m**(3)/mol] - Molar volume in vapour phase

# The compressibility factor is given by
Z_vap_2 = (P_new*V_vap_2)/(R*T);			# For liquid phase
Z_liq_2 = (P_new*V_liq_2)/(R*T);			# For vapour phase

# For vapour phase
f_P_vap_2 = math.exp((Z_vap_2-1) - math.log(Z_vap_2-((P_new*b)/(R*T))) - (a/(2*2**(1./2)*b*R*T))*math.log((Z_vap_2+(1+2**(1./2))*((P_new*b)/(R*T)))/(Z_vap_2+(1-2**(1./2))*((P_new*b)/(R*T)))));
# For liquid phase
f_P_liq_2 = math.exp((Z_liq_2-1) - math.log(Z_liq_2-((P_new*b)/(R*T))) - (a/(2*2**(1./2)*b*R*T))*math.log((Z_liq_2+(1+2**(1./2))*((P_new*b)/(R*T)))/(Z_liq_2+(1-2**(1./2))*((P_new*b)/(R*T)))));

# Therefore f_liq/f_vap can be calculated as
fL_fV_2 = (f_P_liq_2/f_P_vap_2);

# And new pressure is given by
P_new_prime = P_new*(f_P_liq_2/f_P_vap_2);			#[Pa]
P_new_prime = P_new_prime*10**(-5);

# Since the change in pressure is small, so we can take this to be the vapour pressure at 150 C

# Results
print " The vapour pressure of water using Peng-Robinson equation of stste  is %f bar"%(P_new_prime);

 The vapour pressure of water using Peng-Robinson equation of stste  is 4.675976 bar