Chapter 8 : Phase Equilibria III Phase Diagrams¶

Example 8.1 Page No : 470¶

In [5]:
from scipy.optimize import fsolve
import math

# Variables
A_C5H12 = 9.2131 ; 			#From table E8.2A
B_C5H12 = 2477.07 ; 			#From table E8.2A
C_C5H12 = -39.94 ; 			#From table E8.2A
A_C6H12 = 9.1325 ; 			#From table E8.2A
B_C6H12 = 2766.63 ; 			#From table E8.2A
C_C6H12 = -50.50 ; 			#From table E8.2A
A_C6H14 = 9.2164 ; 			#From table E8.2A
B_C6H14 = 2697.55 ; 			#From table E8.2A
C_C6H14 = -48.78 ; 			#From table E8.2A
A_C7H16 = 9.2535 ; 			#From table E8.2A
B_C7H16 = 2911.32 ; 			#From table E8.2A
C_C7H16 = -56.51 ; 			#From table E8.2A
x_C5H12 = 0.3 ;
x_C6H12 = 0.3 ;
x_C6H14 = 0.2 ;
x_C7H16 = 0.2 ;

# Calculations
def f82(T):
return -1 + (x_C5H12 * math.exp(A_C5H12 - B_C5H12 / (T + C_C5H12)) + x_C6H12 *\
math.exp(A_C6H12 - B_C6H12 / (T + C_C6H12)) + x_C6H14 * \
math.exp(A_C6H14 - B_C6H14 / (T + C_C6H14)) + x_C5H12 * math.exp(A_C5H12 - B_C5H12 / (T \
+ C_C5H12)) + x_C7H16 * math.exp(A_C7H16 - B_C7H16 / (T + C_C7H16)));

y =fsolve(f82,300)

# Results
print "    The temperature at which the liquid develops the first bubble of vapour = %d K"%(y);

# Note : answer may vary because of fsolve function of python.

    The temperature at which the liquid develops the first bubble of vapour = 317 K


Example 8.2 Page No : 471¶

In [6]:
from scipy.optimize import fsolve
import math

# Variables
A_C5H12 = 9.2131 ; 			#From table E8.2A
B_C5H12 = 2477.07 ; 			#From table E8.2A
C_C5H12 = -39.94 ; 			#From table E8.2A
A_C6H12 = 9.1325 ; 			#From table E8.2A
B_C6H12 = 2766.63 ; 			#From table E8.2A
C_C6H12 = -50.50 ; 			#From table E8.2A
A_C6H14 = 9.2164 ; 			#From table E8.2A
B_C6H14 = 2697.55 ; 			#From table E8.2A
C_C6H14 = -48.78 ; 			#From table E8.2A
A_C7H16 = 9.2535 ; 			#From table E8.2A
B_C7H16 = 2911.32 ; 			#From table E8.2A
C_C7H16 = -56.51 ; 			#From table E8.2A

y_C5H12 = 0.3 ;
y_C6H12 = 0.3 ;
y_C6H14 = 0.2 ;
y_C7H16 = 0.2 ;
P = 1 ; 			#[bar]

# Calculations and Results
def f83(T):
return  -1 + P * ( y_C5H12 / math.exp(A_C5H12 - B_C5H12 / (T +
C_C5H12)) + y_C6H12 / math.exp(A_C6H12 - B_C6H12 / (T + C_C6H12)) + \
y_C6H14 / math.exp(A_C6H14 - B_C6H14 / (T + C_C6H14)) + y_C7H16 / math.exp(A_C7H16 - B_C7H16 / (T + C_C7H16)));

y =fsolve(f83,300)

print "The temperature at which vapour develops the first drop of liquid = %.f K"%(y) ;

T = y ;
P_sat_C5H12 = math.exp(A_C5H12 - B_C5H12 / (T + C_C5H12)) ;
p_sat_C6H12 = math.exp(A_C6H12 - B_C6H12 / (T + C_C6H12)) ;
P_sat_C6H14 = math.exp(A_C6H14 - B_C6H14 / (T + C_C6H14)) ;
P_sat_C7H16 = math.exp(A_C7H16 - B_C7H16 / (T + C_C7H16)) ;

x_C5H12 = y_C5H12 * P / P_sat_C5H12 ;
x_C6H12 = y_C6H12 * P / p_sat_C6H12 ;
x_C6H14 = y_C6H14 * P / P_sat_C6H14 ;
x_C7H16 = y_C7H16 * P / P_sat_C7H16 ;

print "x_C5H12 = %.3f     x_C6H12 = %.3f         x_C6H14 = %.3f     x_C7H16 = %.3f"%(x_C5H12,x_C6H12,x_C6H14,x_C7H16) ;

The temperature at which vapour develops the first drop of liquid = 349 K
x_C5H12 = 0.091     x_C6H12 = 0.345         x_C6H14 = 0.159     x_C7H16 = 0.405


Example 8.3 Page No : 471¶

In [8]:
from scipy.integrate import *
from numpy import *
import math

# Variables
P_a_sat = 0.53 ;  			#[bar]
P_b_sat = 0.16 ; 			#[bar]
X = 1./3 ;
Y = 1- X ;
x_a_feed = 0.5 ;
x_b_feed = 0.5 ;

# Calculations
a = Y * -(x_a_feed + x_b_feed) + Y**2 ;
b = X * Y *(P_a_sat + P_b_sat) - (x_a_feed * P_b_sat + x_b_feed * P_a_sat)*X ;
c =  P_a_sat * P_b_sat * X**2;

k=poly1d([a,b,c])#'k',0);
#P = c + b*k**1 + a*k**2 ;
M = roots(k);

# Results
for i in range(2):
ans = sign(M[0]) ;
if ans == 1:
print "     Pressure  = %.2f bar"%(M[i]) ;
Xa = x_a_feed / (P_a_sat / M[0] * X + Y) ;			#....E8.4D
Ya = Xa * P_a_sat / M[0]                			#.....E8.4B
print "       Xa = %.2f        Ya = %.1f"%(Xa,Ya);
break

     Pressure  = 0.31 bar
Xa = 0.40        Ya = 0.7


Example 8.5 Page No : 476¶

In [5]:
from scipy.optimize import fsolve
import math

# Variables
P1_sat = 0.72 ; 			#[bar]
P2_sat = 0.31 ; 			#[bar]
A = 3590. ;
B = -1180. ;
R = 8.314 ;
T = 70. + 273 ;			#[K]

# Calculations
def f85(x1):
return -.48 + ( x1 * math.exp((A + 3*B) * (1 - x1)**2 / (R * T) -
4 * B * (1 - x1)**3  / (R * T)) * P1_sat) / ( x1 * math.exp((A + 3*B) *
(1 - x1)**2 / (R * T) -4 * B * (1 - x1)**3  / (R * T)) * P1_sat +(1 - x1) *
math.exp((A - 3*B) * x1**2 / (R * T) -4 * B * x1**3  / (R * T)) * P2_sat ) ;

y = fsolve(f85,0.1)[0]
x1 = y ;
P = ( x1 * math.exp((A + 3*B) * (1 - x1)**2 / (R * T) - 4 * B *
(1 - x1)**3  / (R * T)) * P1_sat) + (1 - x1) * math.exp((A - 3*B) * x1**2 /
(R * T) -4 * B * x1**3  / (R * T)) * P2_sat ;

# Results
print "The value of x1 = %.3f"%(y) ;
print "Pressure = %.2f bar"%(P) ;

The value of x1 = 0.114
Pressure = 0.55 bar


Example 8.7 Page No : 488¶

In [6]:
# Variables
import math
P = 0.223 ; 			#[bar]
P_a_sat = 0.156 ; 			# [bar]
P_b_sat = 0.124 ; 			#[bar]
R = 8.314 ;
T = 50 + 273 ;
Xa = 0.554 ;
Xb = 1 - Xa ;

# Calculations
gama_a = P / P_a_sat ;
A1 = R * T * math.log(gama_a) / (Xb**2) * 10**-3 ;
gama_b = P / P_b_sat ;
A2 = R * T * math.log(gama_b) / (Xa**2) * 10**-3 ;

A = math.ceil((A1 + A2) / 2) ;

# Results
print "Value of two suffix Marguels parameter = %.1f kJ/mol"%(A);

Value of two suffix Marguels parameter = 5.0 kJ/mol


Example 8.9 Page No : 491¶

In [10]:
# Variables
import math
R = 8.314 ;
T = 10 + 273 ; 			#[K]
A_B = 9.2806 ; 			# From Appendix A , Table A1.1
B_B = 2788.5 ; 			# From Appendix A , Table A1.1
C_B = -52.36 ; 			# From Appendix A , Table A1.1
A_C = 9.1325 ; 			# From Appendix A , Table A1.1
B_C = 2766.63 ; 			# From Appendix A , Table A1.1
C_C = -50.50 ; 			# From Appendix A , Table A1.1

x1 = [0 ,0.0610 ,0.2149 ,0.3187 ,0.4320 ,0.5246 ,0.6117 ,0.7265 ,0.8040 ,0.8830 ,0.8999 ,1] ; 			#From table E8.9A
P_exp = [6344 ,6590 ,6980 ,7140 ,7171 ,7216 ,7140 ,6974 ,6845 ,6617 ,6557 ,6073] ; 			#From table E8.9A

P_1_sat = 6072.15 ; 			#[Pa]
P_2_sat = 6344 ; 			#[Pa]

A = [1390 ,1391 ,1392 ,1393 ,1394 ,1395 ,1396 ,1397 ,1398 ,1399 ,1400 ,1401 ,1402 ,1403 ,1404 ,1405 ,1406 ,1407 ,1408 ,1409 ,1410 ] ;

C = zeros([21,12])
for k in range(21):
y = A[k] ;
P = zeros(12)
for i in range(12):
P[i] = x1[i] * math.exp( y / (R * T ) * (1 - x1[i])**2) * P_1_sat+(1 - x1[i]) * math.exp(y / (R * T ) * x1[i]**2) * P_2_sat ;
C[k,i] = (P[i] - P_exp[i])**2 ;

R = zeros(21)
for k in range(21):
y = 0 ;
for i in range(12):
y = y + C[k,i] ;
R[k] = y ;

k = 100000. ;
for i in range(21):
K = R[i] ;
if K < k:
k = K ;

for i in range(21):
if R[i] == k :
print "The two suffix Margules co-efficient is = %g J/mol" %(A[i]) ;

#Note : answer may vary because of rounding off error.

The two suffix Margules co-efficient is = 1401 J/mol


Example 8.10 Page No : 494¶

In [3]:
# Variables
import math
R = 8.314 ;
T = 10 + 273.15			#[K]
A_B = 9.2806 ; 			# From Appendix A , Table A1.1
B_B = 2788.5 ; 			# From Appendix A , Table A1.1
C_B = -52.36 ; 			# From Appendix A , Table A1.1
A_C = 9.1325 ; 			# From Appendix A , Table A1.1
B_C = 2766.63 ; 			# From Appendix A , Table A1.1
C_C = -50.50 ; 			# From Appendix A , Table A1.1

x1 = [0,0.0610 ,0.2149 ,0.3187 ,0.4320 ,0.5246 ,0.6117 ,0.7265 ,0.8040 ,0.8830 ,0.8999 ,1] ; 			#From table E8.9A
P_exp = [6344 ,6590 ,6980 ,7140 ,7171 ,7216 ,7140 ,6974 ,6845 ,6617 ,6557 ,6073] ; 			#From table E8.9A

P_1_sat = 6073. ; 			#[Pa]
P_2_sat = 6344. ; 			#[Pa]
A = linspace(1390,1410,21) ;
B = linspace(60,80,21)
w = 1 / (R * T) ;
S = zeros([21,21])
for k in range(21):
y = A[k] ;
for i in range(21):
z = B[i] ;
P = zeros(12)
R = zeros(12)
for j in range(12):
P[j] = x1[j] * math.exp((y + 3 * z) * (1 - (x1[j]))**2 *w-4*z*(1-x1[j])**3* w )* P_1_sat + (1-x1[j])*exp((y -3*z)*(x1[j])**2 * w + 4 * z * (x1[j]**3) * w )*P_2_sat ;
R[j] =(P[j] - P_exp[j])**2 ;
m = 0 ;
for l in range(12):
m = m + R[l]
S[k,i] = m ;

D = zeros(21)
for i in range(21):
k = S[i,0]
for l in range(1,21):
if S[i,l] < k:
k = S[i,l]
D[i] = k ;

a = D[0] ;
for i in range(2,21):
if D[i] < a:
a = D[i] ;

# Results
for i in range(21):
if D[i] == a :
for l in range(21):
if S[i,l] == a:
print "     A = %g J/mol"%(A[i]);
print "     B = %g J/mol"%(B[l]) ;

     A = 1397 J/mol
B = 69 J/mol


Example 8.11 Page No : 494¶

In [1]:
%matplotlib inline

from matplotlib.pyplot import *

# Given
R = 8.314 ;
T = 10 + 273.15 ; 			#[K]
x1 = [0 ,0.0610 ,0.2149 ,0.3187 ,0.4320 ,0.5246 ,0.6117 ,0.7265 ,0.8040 ,0.8830 ,0.8999 ,1] ; 			#From table E8.9A
P_exp = [6344 ,6590 ,6980 ,7140 ,7171 ,7216 ,7140 ,6974 ,6845 ,6617 ,6557 ,6073 ,6073] ; 			#From table E8.9A
y1 = [ 1 ,0.0953 ,0.2710 ,0.3600 ,0.4453,0.5106 ,0.5735 ,0.6626 ,0.7312 ,0.8200 ,0.8382, 0 ] ;			#From table E8.9A
P_1_sat = 6073. ; 			#[Pa]
P_2_sat = 6344. ; 			#[Pa]

n = 0 ;
x2 = zeros(11)
y2 = zeros(11)
g_E = zeros(11)
ydata = zeros(11)
xdata = zeros(11)
for i in range(1,11):
x2[i] = 1 - x1[i] ;
y2[i] = 1 - y1[i] ;
g_E[i] = R * T *( x1[i] * math.log (( y1[i] * P_exp[i]) / (x1[i]* P_1_sat)) + x2[i] * math.log((y2[i] * P_exp[i]) / (x2[i] * P_2_sat)) ) ;
n = n + g_E[i] / ((x1[i] * x2[i]) * 10) ;
ydata[i-1] = (g_E[i]/(x1[i]*x2[i]));
xdata[i-1] = x1[i] - x2[i] ;

m= 0
n=0
o = 0
p= 0
N = 10
for i in range(1,11):
m = m + g_E[i] * (2 * x1[i] - 1) / ( x1[i] *  x2[i]) ;
n = n + g_E[i] / ( x1[i] *  x2[i]) ;
o = o + (2 * x1[i] - 1)  ;
p = p + (2 * x1[i] - 1)**2 ;

x_bar = o / N ;
y_bar = n / N ;
a1 = (N * m - n * o)/(N * p - o**2) ;
a0 = y_bar - a1 * x_bar ;

ydata2 = zeros(11)
for i in range(11):
ydata2[i] = a0 + a1*xdata[i] ;

# Results
plot(xdata,ydata,"+") ;
plot(xdata,ydata2) ;
xlabel("x1-x2")
ylabel("g**E/x1x2 [J/mol]")
#suptitle("Fit of g E /x 1 x 2 vs. x 1 2 x 2 for benzene (1)â€“cyclohexane (2) at 10Â°C to straight line")
ylim(1320,1500)

print "      From average,  the value of A = %d J/mol"%(n/10) ;
print "      From linear regression best fit line the values of A and B are \
%.1f J/mol    &     %.1f J/mol   respectively ."%(a0, a1) ;
#Readers can refer figure E8.11 .
show()

Populating the interactive namespace from numpy and matplotlib
From average,  the value of A = 1408 J/mol
From linear regression best fit line the values of A and B are     1401.3 J/mol    &     76.1 J/mol   respectively .


Example 8.12 Page No : 500¶

In [15]:
# Variables
H_O2 = 44253.9			#[bar] , From table 8.1
p_O2 = 0.21  			#[bar]

# Calculations
x_O2 = p_O2 / H_O2 ;
v_H2O = 1/(1/0.001 * 1/0.018 * 0.001 );
_O2_ = x_O2 / v_H2O ; 			#[M]

# Results
print "Mole fraction of O2 = %.2e"%(x_O2 ) ;
print "Concentration of O2 = %.2e M "%(_O2_) ;

Mole fraction of O2 = 4.75e-06
Concentration of O2 = 2.64e-04 M


Example 8.13 Page No : 500¶

In [13]:
# Variables
import math
P = 300. ; 			#[bar]
V_bar_inf_N2 = 3.3 * 10**-5 ;
R = 8.314 ;
T = 298. ; 			#[K]
y_N2 = 1. ; 			# At 25*C vapour pressure of water is small
H_N2_1 = 87365. ; 			#[bar]
P_c = 33.8 ; 			#[bar]
T_c = 126.2 ;			# [K]
w = 0.039 ; 			# From Appendix A.1
log_w_0 = 0.013 ;
log_w_1 = 0.210 ;

# Calculations
H_N2_300 = H_N2_1 * math.exp((V_bar_inf_N2 * (P -1) * 10**5 )/ (R * T)) ;
k = log_w_0 + w * log_w_1 ;
sai_N2 = 10**k ;
x_N2 = y_N2 * sai_N2 * P / H_N2_300 ;

# Results
print "Solubility of N2 in water = %.5f"%(x_N2) ;

Solubility of N2 in water = 0.00242


Example 8.20 Page No : 518¶

In [2]:
%matplotlib inline

import math
from matplotlib.pyplot import *
from numpy import *

# Variables
R = 8.314 ;
T = 20. + 273 ;			#[K]
A = 6000. ; 			#[J/mol]
B = -384. ; 			#[J/mol]
x_a = [0.001 ,0.03 ,0.05 ,0.06 ,0.075 ,0.1 ,0.12 , 0.13 ,0.15 ,0.2 ,0.25 ,0.3 ,0.35 ,0.4 ,0.45,0.475 ,0.5 ,0.55 ,0.6 ,0.65  ,0.7 ,0.75 ,0.8 ,0.8475 ,0.85  ,0.9 ,0.925 ,0.95 ,0.975 ,0.999] ;

y_data = zeros(30)
y_data2 = zeros(30)

# Calculations
for i in range(30):
y_data[i] = R * T * ( x_a[i] * math.log(x_a[i]) + (1 - x_a[i]) * math.log(1- x_a[i])) + x_a[i] * (1 - x_a[i]) * (A + B * (2*x_a[i] - 1 )) ;
y_data2[i] = - 82 * x_a[i]- 185.6 ;

m = min(y_data) ;
for i in range(30):
if y_data[i] == m:
a = x_a[i] ;

for i in range(30):
y_data2[i] = -(R * T *( math.log(a)  - math.log(1 - a)) + A * (1 - 2*a) + B * (6 * a - 1 - 6 * a**2)) * (x_a[i] - a) + m ;

y_data3 = zeros(20)
for i in range(20):
y_data3[i] = y_data[i] - y_data2[i] ;
n = min(y_data3) ;

for i in range(20):
if y_data3[i] == n:
b = x_a[i] ;

# Results
plot(x_a ,y_data) ;
plot(x_a ,y_data2) ;
xlabel("Xa")
ylabel("g-XaGa-XbGb")
suptitle("the binary system described ")
print "The equilibrium composition can be found by drawing a line tangent to\
the minima .  In this case the answer is %.2f and %.1f   ." %(a ,b)
show()

Populating the interactive namespace from numpy and matplotlib
The equilibrium composition can be found by drawing a line tangent to the minima .  In this case the answer is 0.85 and 0.1   .

WARNING: pylab import has clobbered these variables: ['draw_if_interactive']
%pylab --no-import-all prevents importing * from pylab and numpy


Example 8.21 Page No : 519¶

In [15]:
from scipy.optimize import fsolve
import math

# Variables
A = 6349 ; 			#[J/mol]
B = -384. ; 			#[J/mol]
R = 8.314 ;
T = 20 + 273 ; 			#[K]
k = 0.000001 ;

# Calculations
def f816(x_a):
return R * T * (1/x_a + 1/(1 - x_a)) - 2 * A +6 * B * (1 - 2 * x_a) + k

ans1 = fsolve(f816,0.1)[0]
ans2 = fsolve(f816,0.5)[0]

print "%.3f < x_a < %.3f "%(ans1,ans2)

0.225 < x_a < 0.706


Example 8.22 Page No : 522¶

In [1]:
from scipy.optimize import fsolve

# Variables
import math
T = 300.    			#[K]
A = 6235.  	    		#[J/mol]
P_a_sat = 100. * 10**3 ; 			#[Pa]
P_b_sat = 50. * 10**3 ; 			#{Pa}
R = 8.314 ;
w = 1./(R * T)

def f817(R):
x_a_a = R[0]
x_a_b = R[1]
Z817 = [0,0]
Z817[0] = x_a_b * math.exp(A * (1 - x_a_b) ** 2 * w) - x_a_a * math.exp(A*(1-x_a_a)**2*w)
Z817[1] = (1 - x_a_b) * math.exp(A * ( x_a_b) ** 2 * w) - (1 - x_a_a) * math.exp(A * (x_a_a) ** 2 * w )
return Z817
x0 = [0.75 , 0.1] ;
z = fsolve(f817,x0)

# Results
print "The compositions are : x_a_a = %.3f and x_a_b = %.3f"%(z[0],z[1]) ;
P = z[0] * math.exp(A * z[1] ** 2 * w) * P_a_sat + z[1] * math.exp(A * z[0] ** 2 * w) * P_b_sat ;
print "Total pressure = %d kPa"%(P * 10**-3) ;
y_a = z[0] * math.exp(A * z[1] ** 2 * w) * P_a_sat / P ;
print "y_a = %.3f" %( y_a ) ;

# Note : incorrect answer in textbook

The compositions are : x_a_a = 0.272 and x_a_b = 0.272
Total pressure = 49 kPa
y_a = 0.667


Example 8.23 Page No : 535¶

In [17]:
# Variables
T_b = 373.15 ; 	    		#[K]
del_h_vap = 2257. ; 		#[J/g]
MW_salt = 58.5 ; 			#[g/mol]
MW_water = 18. ; 			#[g/mol]
w_salt = 3.5 ;
w_water = 100 - w_salt ;
R = 8.314 ;

# Calculations
x_salt = (w_salt / MW_salt) / (w_salt / MW_salt + w_water / MW_water) ;
x_b = 2 * x_salt
del_T = R * T_b**2 / (del_h_vap * MW_water)* x_b ;

# Results
print "The temperature that sea water boils is = %.2f degreeC"%(100 + del_T);

The temperature that sea water boils is = 100.63 degreeC


Example 8.24 Page No : 538¶

In [18]:
# Variables
rho_w = 1000. 	        		# [kg/m**3]
g = 9.8 ; 	        	    	# [m/s**2]
h = 0.0071 ;			    #[m]
m_b = 1.93 * 10**-3 ; 			# [kg]
V = 520. * 10**-6 ; 			#[m**3]
R = 8.314
T = 298.

# Calculations
PI = rho_w * g * h ;
C_b = m_b / V ;
MW_b = R * T * C_b / PI ;

# Results
print "The molecular weight of the protein = %d kg/mol"%( MW_b );

The molecular weight of the protein = 132 kg/mol