# Chapter 9 : Chemical reaction Equilibria¶

## Example 9.1 Page No : 570¶

In [1]:
# Variables
n_o_CH3OH = 1. ; 			#[mol]
n_o_H2O = 3.     			#[mol]
S = 0.87 ;

# Calculations
n_CH3OH = 1 - S ;
n_H2O = 2 - S ;
n_CO2 = S ;
n_H2 = 3 * S ;
n_v = n_CH3OH + n_CO2 + n_H2O + n_H2 ;
y_H2 = n_H2 / n_v ;

# Results
print "No of moles of H2 produced for 1mol of CH3OH = %.3f mol" %( n_H2)
print "Mole fraction of H2 = %.2f" %( y_H2) ;

No of moles of H2 produced for 1mol of CH3OH = 2.610 mol
Mole fraction of H2 = 0.55


## Example 9.2 Page No : 574¶

In [6]:
import math

# Variables
del_gf_0_CO2 = -394.36 ; 			#[kJ/mol],From Appendix A.3
del_gf_0_H2 = 0 ; 			#[kJ/mol],From Appendix A.3
del_gf_0_H2O = -228.57 ; 			#[kJ/mol],From Appendix A.3
del_gf_0_CH3OH = -161.96 ; 			#[kJ/mol],From Appendix A.3
n_CO2 = 1 ;
n_H2 = 3 ;
n_CH3OH = 1 ;
n_H2O = 1 ;
T = 298.15 ;			#[K]
R = 8.314 ; 			#[J/molK]

# Calculations
del_g0_rxn = (n_CO2 * del_gf_0_CO2 + n_H2 * del_gf_0_H2 - n_H2O * del_gf_0_H2O - n_CH3OH * del_gf_0_CH3OH) * 10**3 ; 			# [J/mol]
K_298 = math.exp( - del_g0_rxn / (R * T)) ;

# Results
print "The equillibrium constant K298 = %.2f "%(K_298) ;

The equillibrium constant K298 = 4.69


## Example 9.3 page no : 575¶

In [1]:
import math

# Variables
K_298 = 4.69

# calculations
deltaH =  (-393.51) +  3 * 0 - (-241.82) - (-200.66)
lnK333_469 = round((-deltaH*1000/8.314 )*(1./333 - 1./298),2)
lnK333 = lnK333_469 + math.log(K_298)
K333 = math.e**lnK333

# Results
print "Equilibrium constant = %.2f "%K333

# Note: answer is slightly different because of rounding off error.

Equilibrium constant = 37.54


## Example 9.4 Page No : 577¶

In [5]:
import math
# Variables
del_gf_0_CH2O = -110.0 ; 			#[kJ/mol],From Appendix A.2 & A.3
del_gf_0_H2 = 0 ; 			#[kJ/mol],From Appendix A.2 & A.3
del_gf_0_CH4O = -162.0 ; 			#[kJ/mol],From Appendix A.2 & A.3
del_hf_0_CH2O = -116.0 ; 			#[kJ/mol],From Appendix A.2 & A.3
del_hf_0_H2 = 0 ; 			#[kJ/mol],From Appendix A.2 & A.3
del_hf_0_CH4O = -200.7 ; 			#[kJ/mol],From Appendix A.2 & A.3n_CH20 = 1 ;
n_H2 = 1. ;
n_CH4O = 1. ;
n_CH2O = 1. ;
T1 = 298.   			#[K]
T2 = 873. ; 			# [K]
R = 8.314 ; 			#[J/molK]
Del_A = 3.302 ;
Del_B = -4.776 * 10**-3 ;
Del_C = 1.57 * 10**-6 ;
Del_D = 0.083 * 10**5 ;

# Calculations and Results
del_g_rxn_298 = n_CH2O * del_gf_0_CH2O + n_H2 * del_gf_0_H2 - n_CH4O * del_gf_0_CH4O ;
K_298 = math.exp( - del_g_rxn_298 * 10**3 / (R * T1)) ;
print "a) K_298 = %.2e    As the equilibrium constant is very small very little amount of\
formaldehyde will be formed ."%(K_298) ;

#Solution(b)
del_h_rxn_298 = (n_CH2O * del_hf_0_CH2O + n_H2 * del_hf_0_H2 - n_CH4O * del_hf_0_CH4O) * 10**3 ;			#[J/mol]
K_873 = K_298 * math.exp((-del_h_rxn_298  * (1/T2 - 1/T1)) / R) ;
print "b)            i)  K_873 = %.2f "%(K_873) ;

#Solution(c)
x =  ( -del_h_rxn_298 / R + Del_A * T1 + Del_B / 2 * T1**2 + Del_C /3 * T1**3 - Del_D / T1 ) *(1/T2 - 1/T1) + Del_A * math.log(T2 / T1) + Del_B / 2 * (T2 -T1) + Del_C / 6 * (T2**2 -T1**2) + Del_D / 2 * (1/(T2**2) -1/(T1**2)) ;
K_873 = K_298 * math.exp(x) ;
print "             ii)  K_873 = %.2f "%(K_873) ;

# Note : answer are slightly different because of rounding off error.

a) K_298 = 7.67e-10    As the equilibrium constant is very small very little amount of formaldehyde will be formed .
b)            i)  K_873 = 4.61
ii)  K_873 = 8.71


## Example 9.6 Page No : 581¶

In [1]:
# Variables
import math
del_g0_f_C6H6 = -32.84 ; 			#[kJ/mol] , From Table E9.6
del_g0_f_C2H4 = 68.15 ; 			#[kJ/mol] , From Table E9.6
del_g0_f_H2 = 0 ;  			        #[kJ/mol] , From Table E9.6
del_h0_f_C6H6 = -84.68 ; 			#[kJ/mol] , From Table E9.6
del_h0_f_C2H4 = 52.26 ; 			#[kJ/mol] , From Table E9.6
del_h0_f_H2 = 0 ;  			        #[kJ/mol] , From Table E9.6
T1 = 298.2 ;			            #[K]
P = 1.                  			#[bar]
R = 8.31 ;
T2 = 1273.               			# [K]

# Calculations
del_g0_f_rxn = del_g0_f_C2H4 + del_g0_f_H2 - del_g0_f_C6H6 ;
K_298 = math.exp ( - (del_g0_f_rxn * 10**3) / (R * T1)) ;
del_h0_f_rxn = (del_h0_f_C2H4 + del_h0_f_H2 - del_h0_f_C6H6) * 10**3 ;
K_1273 = K_298 * math.exp( - del_h0_f_rxn / R * (1/T2 - 1/T1)) ;
x = math.sqrt( K_1273 / ( K_1273 + P)) ;

# Results
print "n_C2H6 = %.2f mol     n_C2H4 = %.2f mol     n_H2 = %.2f mol"%(1-x,x,x) ;

n_C2H6 = 0.09 mol     n_C2H4 = 0.91 mol     n_H2 = 0.91 mol


## Example 9.7 Page No : 583¶

In [10]:
from matplotlib.pyplot import *
from numpy import *
from scipy.integrate import *
import math

# Variables
del_h0_f_NH3 = -46.11 ; 			# [kJ/mol],From table E9.7
del_h0_f_N2 = 0          			# [kJ/mol],From table E9.7
del_h0_f_H2 = 0 ; 		        	# [kJ/mol],From table E9.7
del_g0_f_NH3 = -16.45 ; 			# [kJ/mol],From table E9.7
del_g0_f_N2 = 0 ; 		        	# [kJ/mol],From table E9.7
del_g0_f_H2 = 0 ; 			        # [kJ/mol],From table E9.7
n_NH3 = 2. ;
n_N2 = -1. ;
n_H2 = -3. ;
A_NH3 = 3.578
B_NH3 = 3.02 * 10**-3
D_NH3 = -0.186 * 10**5 ;
A_N2 = 3.280
B_N2 = 0.593 * 10**-3
D_N2 = 0.040 * 10**5 ;
A_H2 = 3.249
B_H2 = 0.422 * 10**-3
D_H2 = 0.083 * 10**5 ;
R = 8.314 ;
T = 298. ;
T2 = 773. ;
P = 1. ; 			#[bat]

# Calculations
Del_h0_rxn = (n_NH3 * del_h0_f_NH3 + n_N2 * del_h0_f_N2 + n_H2 * del_h0_f_H2) * 10**3 ;
Del_g0_rxn = (n_NH3 * del_g0_f_NH3 + n_N2 * del_g0_f_N2 + n_H2 * del_g0_f_H2) * 10**3 ;
del_A = n_NH3 * A_NH3 + n_N2 * A_N2 + n_H2 * A_H2 ;
del_B = n_NH3 * B_NH3 + n_N2 * B_N2 + n_H2 * B_H2 ;
del_D = n_NH3 * D_NH3 + n_N2 * D_N2 + n_H2 * D_H2 ;

K_298 = math.exp( - Del_g0_rxn / ( R * T)) ;
K_T = K_298 * math.exp( - Del_h0_rxn / R * (1 / T2 - 1 / T)) ;
A = K_T * P**2 *27 -16 ;
B = 64 - K_T * P**2 * 108 ;
C = -64 + K_T * P**2 * 162 ;
D = -108 * K_T * P**2 ;
E = 27 * K_T * P**2 ;

mycoeff =[E , D ,C , B ,A];
M = roots([A,B,C,D,E]);

for i in range(3):
ans = isreal(M[i]) ;
if ans == True:
y = M[i] / M[i+1] - 1 ;
ans = sign(y) ;

if ans == 1 :
x = M[i]
else:
x = M[i+1]

# Results
print "a)Extent of reaction = %.3f "%(x);

#(b)
X = (-Del_h0_rxn  / R + del_A * T + del_B / 2 * T**2 - del_D / T) * (1/T2 - 1/T) + del_A * math.log(T2 / T) + del_B / 2 * (T2 - T) + del_D / 2 * (1/(T2**2) - 1/(T**2) );
K_T = K_298 * math.exp(X) ;

A = K_T * P**2 *27 -16 ;
B = 64 - K_T * P**2 * 108 ;
C = -64 + K_T * P**2 * 162 ;
D = -108 * K_T * P**2 ;
E = 27 * K_T * P**2 ;

mycoeff =[E , D ,C , B ,A];
M1 = roots([A,B,C,D,E]);

for i in range(3):
ans = isreal(M1[i])
if ans == True :
y = M1[i] / M1[i+1] - 1 ;
ans = sign(y) ;
if ans == True:
x1 = M1[i]
else:
x1 = M1[i+1]

print "b)  Extent of reaction = %.3f"%(x1);
print "Under these conditions we do not expect to produce an appreciable amount of ammonia ."

a)Extent of reaction = 0.005
b)  Extent of reaction = 0.003
Under these conditions we do not expect to produce an appreciable amount of ammonia .


## Example 9.9 Page No : 585¶

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

# Variables
K_T = 1.51 * 10**-5 ;
P = 300.  ; 			#[bar]
T = 500 + 273.2 ; 			#[K]
R = 8.314 ;

# Calculations and Results
def f991(k):
return ((2 * k)**2 * (4 - 2 * k)**2 / ((1 - k) * (3 - 3*k)**3)) * P**-2 - K_T

z1 = fsolve(f991,0.3)[0]

print "a) Extent of reaction = %.2f "%(z1);

P_c = [111.3 * 101325 , 33.5 * 101325 , 12.8 * 101325] ;
T_c = [405.5 , 126.2 , 33.3] ;

a = zeros(3)
b = zeros(3)
V = zeros(3)
sai = zeros(3)
for i in range(3):
a[i] = 27./ 64 * (R * T_c[i])**2 / P_c[i] ;
b[i] = (R * T_c[i]) / (8 * P_c[i]) ;

def f992(v):
return  (R * T) / (v - b[i]) - a[i] / (v**2) - P * 100000 ;
V[i] = fsolve(f992,.0002)
sai[i] = math.exp( - math.log((V[i] - b[i]) * P * 10**5/ ( R * T)) + b[i] / (V[i] - b[i]) - 2 * a[i] / (R * T * V[i])) ;

def f993(k):
return ((2 * k)**2 * sai[0]**2 * (4 - 2 * k)**2 * 3 / ((1 - k) * sai[1]* (3 - 3*k)**3 * sai[2]**3 ))* P**-2  - K_T

z2 = fsolve(f993,0.3)

x = (z1 - z2) / z1 * 100 ;

print "b)        Extent of reaction = %.2f "%(z2);
print " A correction of about   %d%% is observed from accounting for nonideal behaviour . "%(x)

a) Extent of reaction = 0.37
b)        Extent of reaction = 0.33
A correction of about   11% is observed from accounting for nonideal behaviour .


## Example 9.10 Page No : 586¶

In [3]:
# Variables
import math
del_g0_f_1 = 31.72 ; 			#[kJ/mol]
del_g0_f_2 = 26.89 ; 			#[kJ/mol]
R = 8.314 ;
T = 298 ;			#[K]

# Calculations
del_g0_rxn = del_g0_f_2 - del_g0_f_1 ;
K = math.exp( - del_g0_rxn * 10**3 / (R * T) ) ;
x = K / (1 + K) ;

# Results
print "x = %.3f  \nAt equilibrium %.1f %% of the liquid exists as cyclohexane."%(x ,x * 100) ;

x = 0.875
At equilibrium 87.5 % of the liquid exists as cyclohexane.


## Example 9.11 Page No : 587¶

In [10]:
# Variables
import math
del_g0_f_CaCO3 = -951.25 ;
del_g0_f_CaO = -531.09 ;
del_g0_f_CO2 = -395.81 ;
R = 8.314 ;
T = 1000. ;			# [K]

# Calculations
del_g0_rxn = del_g0_f_CaO + del_g0_f_CO2 - del_g0_f_CaCO3 ;
K = math.exp(-del_g0_rxn * 10**3 / (R * T)) ;
p_CO2 = K ;

# Results
print "Equilibrium pressure = %.3f bar "%(p_CO2) ;

Equilibrium pressure = 0.053 bar


## Example 9.12 Page No : 588¶

In [11]:
# Variables
import math
del_g0_f_B = 124.3 ; 			#[kJ/mol] , From Appendix A.3
del_g0_f_Ac = 209.2 ; 			#[kJ/mol] , From Appendix A.3
R = 8.314 ;
T = 298. ; 			# [K]
A = 9.2806 ;
B = 2788.51 ;
C = -52.36 ;

# Calculations
del_g0_rxn = del_g0_f_B - 3 * del_g0_f_Ac ;
K = math.exp( - del_g0_rxn * 10**3 / (R * T)) ;
P = 1 / K**(1./3) ;
X = A - B / (T + C) ;
P_b = math.exp(X) ;

# Results
print ("At equilibrium , the cylinder is almost completely filled with Benzene .")
print "System pressure = %.3f bar "%(P_b)

At equilibrium , the cylinder is almost completely filled with Benzene .
System pressure = 0.126 bar


## Example 9.13 Page No : 596¶

In [5]:
# Variables
import math
E_0_c = 0.153 ; 			#[V]
E_0_a = -0.521 ; 			# [v]
T = 298. ; 			#[K]
z = 1. ;
F =96485 ; 			#[C/mol e-]
R =8.314 ; 			#[J/mol K ]

# Calculations
E_0_rxn = E_0_c + E_0_a ;
del_g_0_rxn = - z * F * E_0_rxn ;

K = math.exp( - del_g_0_rxn / ( R * T )) ;

# Results
print "The equilibrium constant = %.3g "%(K)
print "The equilibrium constant is small . So the etching will not proceed\
spontaneously . However if we apply work through application of an electrical potential\
, we can etch the copper ."

# Note: answer is slightly different because of rounding off error.

The equilibrium constant = 5.97e-07
The equilibrium constant is small . So the etching will not proceed spontaneously . However if we apply work through application of an electrical potential  , we can etch the copper .


## Example 9.14 Page No : 596¶

In [12]:
# Variables
import math
E_0_c = 0.34 ; 			#[V]
E_0_a = -1.23 			#[V]
T = 298. ;   			# [K]
pH = 1. ;
z = 2. ;
Cu2 = 0.07
F = 96485. 	    		#[C/mol e-]
R = 8.314 ;

# Calculations
E_0_rxn = E_0_c + E_0_a ;
E = E_0_rxn + 2.303 * R * T * 2 * pH / (z * F) + R * T * math.log(Cu2) / (z * F) ;

# Results
print "Del_E_0_rxn = %.2f "%(E_0_rxn ) ;
print "We have to apply potential greater than %.1f V"%(-E) ;

Del_E_0_rxn = -0.89
We have to apply potential greater than 0.9 V


## Example 9.17 Page No : 602¶

In [14]:
# Variables
m = 4. ;
T = 2. ;
Pai = 1. ;
S = 1. ;

# Calculations
R = m - T + 2 - Pai - S ;

# Results
print "We must specify %g independent equations ."%(R)

We must specify 2 independent equations .


## Example 9.19 Page No : 603¶

In [2]:
%matplotlib inline

from numpy import *
from scipy.optimize import fsolve
from matplotlib.pyplot import *
import math

# Variables
del_g_f_CH4 = -50.72 ;
del_g_f_H2 = 0 ;
del_g_f_H2O = -228.57 ;
del_g_f_CO = -137.17 ;
del_g_f_CO2 = -394.36 ;
del_h_f_CH4 = -74.81 ;
del_h_f_H2 = 0 ;
del_h_f_H2O = -241.82 ;
del_h_f_CO = -110.53 ;
del_h_f_CO2 = -393.51 ;

v1_CH4 = -1. ;
v1_H2 = 3. ;
v1_H2O = -1. ;
v1_CO = 1. ;
v1_CO2 = 0. ;
v2_CH4 = -1. ;
v2_H2 = 4. ;
v2_H2O = -2. ;
v2_CO = 0. ;
v2_CO2 = 1. ;

A_CH4 = 1.702 ;
B_CH4 = 9.08 * 10**-3 ;
C_CH4 = -2.16 * 10**-6 ;
D_CH4 = 0 ;
A_H2 = 3.249 ;
B_H2 = 4.22 * 10**-4 ;
C_H2 = 0 ;
D_H2 = 8.30 * 10**3 ;
A_H2O = 3.47 ;
B_H2O = 1.45 * 10**-3 ;
C_H2O = 0 ;
D_H2O = 1.21 * 10**4 ;
A_CO = 3.376 ;
B_CO = 5.57 * 10**-4 ;
C_CO = 0 ;
D_CO = -3.10 * 10**3 ;
A_CO2 = 5.457 ;
B_CO2 = 1.05 * 10**-3 ;
C_CO2 = 0 ;
D_CO2 = -1.16 * 10**5 ;
M = zeros((12,10))
M[:,0] = linspace(600,1150,12)

R = 8.314 ;
P = 1  ; 			#[bar]
T_ref = 298.15 ; 			#[K]

# Calculations
del_g_f_1 = (v1_CO * del_g_f_CO + v1_H2 *del_g_f_H2 + v1_CH4 * del_g_f_CH4 + v1_H2O * del_g_f_H2O) * 1000  ;
del_h_f_1 = (v1_CO * del_h_f_CO + v1_H2 *del_h_f_H2 + v1_CH4 * del_h_f_CH4 + v1_H2O * del_h_f_H2O) * 1000 ;
del_g_f_2 = (v2_CO2 * del_g_f_CO2 + v2_H2 *del_g_f_H2 + v2_CH4 * del_g_f_CH4 + v2_H2O * del_g_f_H2O) * 1000 ;
del_h_f_2 = (v2_CO2 * del_h_f_CO2 + v2_H2 *del_h_f_H2 + v2_CH4 * del_h_f_CH4 + v2_H2O * del_h_f_H2O) * 1000;
Del_A_1 =  v1_CO * A_CO + v1_H2 * A_H2 + v1_CH4 * A_CH4 + v1_H2O * A_H2O  ;
Del_B_1 =  v1_CO * B_CO + v1_H2 * B_H2 + v1_CH4 * B_CH4 + v1_H2O * B_H2O  ;
Del_C_1 =  v1_CO * C_CO + v1_H2 * C_H2 + v1_CH4 * C_CH4 + v1_H2O * C_H2O  ;
Del_D_1 =  v1_CO * D_CO + v1_H2 * D_H2 + v1_CH4 * D_CH4 + v1_H2O * D_H2O  ;
Del_A_2 = v2_CO2 * A_CO2 + v2_H2 * A_H2 + v2_CH4 * A_CH4 + v2_H2O * A_H2O  ;
Del_B_2 = v2_CO2 * B_CO2 + v2_H2 * B_H2 + v2_CH4 * B_CH4 + v2_H2O * B_H2O  ;
Del_C_2 = v2_CO2 * C_CO2 + v2_H2 * C_H2 + v2_CH4 * C_CH4 + v2_H2O * C_H2O  ;
Del_D_2 = v2_CO2 * D_CO2 + v2_H2 * D_H2 + v2_CH4 * D_CH4 + v2_H2O * D_H2O  ;

K_298_1 = math.exp( - del_g_f_1 / (R * T_ref)) ;
K_298_2 = math.exp( - del_g_f_2 / (R * T_ref)) ;

for i in range(12):
X = (-del_h_f_1 / R + Del_A_1 * T_ref + Del_B_1 / 2 * T_ref**2 + \
Del_C_1 /3* T_ref**3- Del_D_1 / T_ref) * (1./M[i,0] - 1./T_ref) + \
Del_A_1*math.log(M[i,0] / T_ref)+ Del_B_1 / 2 * (M[i,0] - T_ref) + \
Del_C_1 / 6 *(M[i,0]**2 - T_ref**2) + Del_D_1 / 2 * (1./(M[i,0]**2) - 1./(T_ref**2))

M[i,1] = K_298_1 * math.exp(X) ;

Y = (-del_h_f_2 / R + Del_A_2 * T_ref + Del_B_2 / 2 * T_ref**2 + Del_C_2/3* T_ref**3- Del_D_2 / T_ref) * \
(1/M[i,0] - 1/T_ref) + Del_A_2 * math.log(M[i,0] / T_ref)+ Del_B_2 / 2 * (M[i,0] - T_ref) + \
Del_C_2 / 6 *(M[i,0]**2 - T_ref**2) + Del_D_2 / 2* (1/(M[i,0]**2) - 1/(T_ref**2));

M[i,2] = K_298_2 * math.exp(Y) ;
def f918(R):
s1 = R[0] ;
s2 = R[1] ;
y = [0,0]
y[0] = (s1 * (3 * s1 + 4 * s2)**3) / ((5 + 2 * s1 + 2 * s2)**2 * (1 - s1 -s2) * (4 - s1 - 2 * s2)) * P**2 - M[i,1] ;
y[1] = (s2 * (3 * s1 + 4 * s2)**4) / ((5 + 2 * s1 + 2 * s2)**2 * (1 - s1 -s2) * (4 - s1 - 2 * s2)**2) * P**2 - M[i,2] ;
return y
z = fsolve(f918,[0.0001,0.0001])
M[i,3] = z[0] ;
M[i,4] = z[1] ;
M[i,5] = (1 - M[i,3] - M[i,4]) / (5 + 2 * M[i,3] + 2 * M[i,4]) ;
M[i,6] = (4 - M[i,3] - 2 * M[i,4]) / (5 + 2 * M[i,3] + 2 * M[i,4]) ;
M[i,7] = (3 * M[i,3] + 4 * M[i,4]) / (5 + 2 * M[i,3] + 2 * M[i,4]) ;
M[i,8] = M[i,3] / (5 + 2 * M[i,3] + 2 * M[i,4]) ;
M[i,9] = M[i,4] / (5 + 2 * M[i,3] + 2 * M[i,4]) ;

n1 = zeros([12,7])
for i in range(12):
for j in range(7):
n1[i,j] = M[i,j] ;
n2 = zeros([12,3])
for i in range(12):
for j in range(3):
n2[i,j] = M[i,j+7]
print "     T              K1               K2          S1          S2       y_CH4     y_H2"
for row in n1:
print "   %5d       %7.2e       %7.2e      %7.3f     %7.3f    %7.3f    %7.3f "%(row[0],row[1],row[2],row[3],row[4],row[5],row[6])
#print (n1) ;
print (" y_H20          y_CO          y_CO2  ") ;
for row in n2:
print " %7.3f    %7.3f    %7.3f "%(row[0],row[1],row[2])

#print (n2) ;
N = zeros([10,10])
for i in range(10):
for j in range(10):
N[i,j] = M[i,j]
'''
plot(N[3],N[0],"+") ;
plot(N[4], N[0],".") ;
plot(N[5] , N[0], "o-") ;
plot(N[6] , N[0], "s-");
plot(N[7] , N[0], "*-") ;
plot(N[8] , N[0], "x-") ;
plot(N[9] , N[0], ".-") ;
show()
'''

suptitle("Figure E9.18    Extent of reaxn vs temp")
xlabel("Temperature(K)")
ylabel("S") ;
#legend("S1","S2") ;
plot(N[:,0] , N[:,5], "o-") ;
plot(N[:,0] , N[:,6], "s-");
plot(N[:,0] , N[:,7], "^-") ;
plot(N[:,0] , N[:,8], "x-") ;
plot(N[:,0] , N[:,9], ".-") ;

# Note : some answers are different because of rounding off error.

Populating the interactive namespace from numpy and matplotlib
T              K1               K2          S1          S2       y_CH4     y_H2
600       4.91e-07       1.40e-05        0.000       0.113      0.170      0.722
650       1.42e-05       2.24e-04        0.003       0.191      0.150      0.671
700       2.59e-04       2.47e-03        0.011       0.294      0.124      0.606
750       3.24e-03       2.01e-02        0.037       0.410      0.094      0.534
800       3.00e-02       1.29e-01        0.099       0.515      0.062      0.461
850       2.15e-01       6.69e-01        0.207       0.578      0.033      0.401
900       1.25e+00       2.93e+00        0.332       0.584      0.012      0.366
950       6.06e+00       1.11e+01        0.424       0.551      0.004      0.356
1000       2.52e+01       3.70e+01        0.485       0.509      0.001      0.358
1050       9.20e+01       1.11e+02        0.531       0.468      0.000      0.362
1100       2.99e+02       3.02e+02        0.319       0.322      0.057      0.484
1150       8.78e+02       7.56e+02        0.343       0.295      0.058      0.488
y_H20          y_CO          y_CO2
0.087      0.000      0.022
0.144      0.000      0.036
0.215      0.002      0.052
0.297      0.006      0.070
0.378      0.016      0.083
0.447      0.032      0.088
0.488      0.049      0.085
0.500      0.061      0.079
0.499      0.069      0.073
0.495      0.076      0.067
0.357      0.051      0.051
0.352      0.055      0.047

/home/jovina/virtualenvs/scipy/local/lib/python2.7/site-packages/scipy/optimize/minpack.py:236: RuntimeWarning: The iteration is not making good progress, as measured by the
improvement from the last ten iterations.
warnings.warn(msg, RuntimeWarning)

Out[2]:
[<matplotlib.lines.Line2D at 0x34b06d0>]

## Example 9.21 Page No : 607¶

In [69]:
# Variables
import math
del_g_0_f_SiCl2 = - 216012. ;
del_g_0_f_SiCl4 = - 492536. ;
del_g_0_f_SiCl3H = -356537. ;
del_g_0_f_SiCl2H2 = -199368. ;
del_g_0_f_SiClH3 = -28482. ;
del_g_0_f_SiH4 = -176152. ;
del_g_0_f_HCl = -102644. ;
del_g_0_f_H2 = 0. ;
del_g_0_f_Si = 0. ;
R = 8.314 ;
T = 1300. ; 			#[K]

# Calculations
Del_g_rxn_1 = del_g_0_f_SiCl2 + 2 * del_g_0_f_HCl - del_g_0_f_SiCl4 - del_g_0_f_H2 ;
Del_g_rxn_2 = del_g_0_f_SiCl3H + del_g_0_f_HCl -  del_g_0_f_SiCl4 - del_g_0_f_H2 ;
Del_g_rxn_3 = del_g_0_f_SiCl2H2 + del_g_0_f_HCl -  del_g_0_f_SiCl3H - del_g_0_f_H2 ;
Del_g_rxn_4 = del_g_0_f_SiClH3 + del_g_0_f_HCl -  del_g_0_f_SiCl2H2 - del_g_0_f_H2 ;
Del_g_rxn_5 = del_g_0_f_SiH4 + del_g_0_f_HCl -  del_g_0_f_SiCl3H - del_g_0_f_H2 ;
Del_g_rxn_6 = del_g_0_f_Si + 4 * del_g_0_f_HCl -  del_g_0_f_SiCl4 - 2 * del_g_0_f_H2 ;

M = zeros([6,4])

M[0,0] = math.exp( - Del_g_rxn_1 / (R * T)) ;
M[1,0] = math.exp( - Del_g_rxn_2 / (R * T)) ;
M[2,0] = math.exp( - Del_g_rxn_3 / (R * T)) ;
M[3,0] = math.exp( - Del_g_rxn_4 / (R * T)) ;
M[4,0] = math.exp( - Del_g_rxn_5 / (R * T)) ;
M[5,0] = math.exp( - Del_g_rxn_6 / (R * T)) ;

S = [0.0763,0.1979,0.0067,0.0001,0.0000,-0.0512] ;
K_cal = [.00137,0.0457,0.00644,0.00181,0.000752,0.000509] ;

for i in range(6):
M[i,1] = S[i] ;
M[i,2] = K_cal[i] ;
M[i,3] = M[i,0] - M[i,2] ;
# Results
print ("       K_i           S      K_i_cal    K_i - K_i_cal") ;
for row in M:
print "    %5.2e     %7.4f    %7.2e    %7.2e"%(row[0],row[1],row[2],row[3])

# Readers can refer figure E9.19 .
# Note : answers are different because of rounding off error.

       K_i           S      K_i_cal    K_i - K_i_cal
1.37e-03      0.0763    1.37e-03    2.77e-06
4.57e-02      0.1979    4.57e-02    -1.95e-05
6.44e-03      0.0067    6.44e-03    2.87e-06
1.81e-03      0.0001    1.81e-03    9.39e-07
7.52e-04      0.0000    7.52e-04    -4.00e-09
5.09e-04     -0.0512    5.09e-04    -3.50e-08