Chapter 14 : Chemical reaction equilibrium

Example 14.1 Page No : 489

In [2]:
import math

# Variables
T = 298.15;			 #temperature in K
del_Gf = [-137.327,-228.600,-394.815,0]; #the standard Gibbs free energy of formation of CO(g),H2O(g),CO2(g) and H2(g) in kJ
n = [1,1,-1,-1]         			 #stoichiometric coefficients of CO(g),H2O(g),CO2(g) and H2(g) respectively (no unit)
R = 8.314;			                 #universal gas constant in J/molK

# Calculations
del_G = (n[0]*del_Gf[0])+(n[1]*del_Gf[1])+(n[2]*del_Gf[2])+(n[3]*del_Gf[3]);
Ka = math.exp((-(del_G*10**3))/(R*T))

# Results
print 'The standard Gibbs free energy of the water gas shift reaction at 298.15K = %0.3f kJ '%(del_G);
print 'The equilibrium constant of the water gas shift reaction at 298.15K = %0.3e '%(Ka);
The standard Gibbs free energy of the water gas shift reaction at 298.15K = 28.888 kJ 
The equilibrium constant of the water gas shift reaction at 298.15K = 8.685e-06 

Example 14.2 Page No : 490

In [3]:
import math
# Variables
T = 298.15;			 #temperature in K
P_s = 0.16716;			 #saturation pressure of CH3OH in bar at T
del_G1 = -166.215
R = 8.314;			 #universal gas constant in J/molK

# Calculations
f_v = 1
f_l = P_s
del_G2 = R*T*math.log(f_v/f_l)*10**-3;			 # Calculations of the value of del_G2 in kJ
del_G = del_G2+del_G1;		

# Results
print 'The standard Gibbs free energy of formation of CH3OHg = %0.3f kJ '%(del_G);
The standard Gibbs free energy of formation of CH3OHg = -161.781 kJ 

Example 14.3 Page No : 491

In [4]:
import math
# Variables
#The water gas shift reaction is given by: CO2(g)+H2(g)--->CO(g)+H2O(g)
T1 = 298.15     			 #initial temperature in K
Ka1 = 8.685*10**-6;			 #equilibrium constant for the water-gas shift reaction at T1 (no unit)
T2 = 1000.		        	 #temperature at which the equilibrium constant has to be determined in K
R = 8.314;			         #universal gas constant in J/molK
del_Hf = [-110.532,-241.997,-393.978,0];			 #the standard enthalpy of formation of CO(g),H2O(g),CO2(g) and H2(g) in kJ
n = [1,1,-1,-1];			 #stoichiometric coefficients of CO(g),H2O(g),CO2(g) and H2(g) respectively (no unit)

# Calculations
#It is assumed that del_H is constant in the temperature range T1 and T2
del_H = (n[0]*del_Hf[0])+(n[1]*del_Hf[1])+(n[2]*del_Hf[2])+(n[3]*del_Hf[3])
Ka2 = Ka1*math.exp(((del_H*10**3)/R)*((1./T1)-(1./T2)));			 

# Results
print 'The equilibrium constant for the water gas shift reaction at 1000K = %f '%(Ka2);
The equilibrium constant for the water gas shift reaction at 1000K = 1.085357 

Example 14.4 Page No : 492

In [22]:
import math
# Variables
P = 0.1     			 #pressure in MPa
T1 = 298.15;			 #initial temperature in K
Ka1 = 8.685*10**-6;			 #equilibrium constant for the water-gas shift reaction at T1 (no unit) (from Example 14.1)
T2 = 1000.  			 #temperature at which the equilibrium constant is to be found, in K
del_H = 41.449;			 #smath.tan(math.radiansard enthalpy of the reaction at T1 in kJ (from Example 14.3)
a = [28.068,28.850,45.369,27.012];
b = [4.631*10**-3,12.055*10**-3,8.688*10**-3,3.509*10**-3]
c = [0,0,0,0];
d = [0,0,0,0];
e = [-0.258*10**5,1.006*10**5,-9.619*10**5,0.690*10**5]
n = [1,1,-1,-1]
R = 8.314;		
Ka2_prev = 1.0855


# Calculations
del_a = (n[0]*a[0])+(n[1]*a[1])+(n[2]*a[2])+(n[3]*a[3]);
del_b = (n[0]*b[0])+(n[1]*b[1])+(n[2]*b[2])+(n[3]*b[3]);
del_c = (n[0]*c[0])+(n[1]*c[1])+(n[2]*c[2])+(n[3]*c[3]);
del_d = (n[0]*d[0])+(n[1]*d[1])+(n[2]*d[2])+(n[3]*d[3]);
del_e = (n[0]*e[0])+(n[1]*e[1])+(n[2]*e[2])+(n[3]*e[3]);
del_H0 = (del_H*10**3)-((del_a*T1)+((del_b/2)*T1**2)+((del_c/3)*T1**3)+((del_d/4)*T1**4)-(del_e/T1));
I = (math.log(Ka1))-((1./R)*((-del_H0/T1)+(del_a*math.log(T1))+((del_b/2)*T1)+((del_c/6)*T1**2)+((del_d/12)*T1**3)+((del_e/(2*T1**2)))));
Ka2 = math.exp(((1./R)*((-del_H0/T2)+(del_a*math.log(T2))+((del_b/2)*T2)+((del_c/6)*T2**2)+((del_d/12)*T2**3)+((del_e/(2*T2**2)))))+I);

			 # Results
print 'The equilibrium constant for the water gas shift reaction at 1000K\
 by taking into account the variation of del_H with temperature = %f '%(Ka2);
print 'The equilibrium constant for the water gas shift reaction at 1000K without\
 considering the variation of del_H with temperature as given by Example14.3 = %0.4f '%(Ka2_prev);
The equilibrium constant for the water gas shift reaction at 1000K by taking into account the variation of del_H with temperature = 0.664976 
The equilibrium constant for the water gas shift reaction at 1000K without considering the variation of del_H with temperature as given by Example14.3 = 1.0855 

Example 14.5 page no : 494

In [103]:
from sympy import Symbol
from sympy.solvers import solve
import math

# variables
deltaGf = -161.781       # kJ from exa. 14.2
deltaG298 = deltaGf - (-137.327)
deltaH298 = -200.660 - (-110.532)
deltaA = 18.382 - 28.068 - 2 * 27.012
deltaB = (101.564 - 4.631 - 2 * 3.509) * 10**-3
deltaC = -28.683 * 10**-6
deltaD = 0
deltaE = (0.258 - 2 * 0.690) * 10**5
T = 298.15

# calculations
deltaHf298 = -238.648 + 37.988
deltaH0 = Symbol('deltaH0')
Eq = deltaH0 + deltaA*T + deltaB*T**2/2 + deltaC*T**3/3 - deltaE/T - deltaH298*1000
deltaH0 = solve(Eq,deltaH0)[0]/1000
I = Symbol('I')

Eq1 = deltaH0*1000 - deltaA*T*math.log(T) - deltaB*T**2/2 -  deltaC*T**3/6 - deltaD*T**4/12 - deltaE/(2*T) - 8.314 * T * I + 24454
I = round(solve(Eq1,I)[0],3)

T = 500
deltaG500 = (deltaH0*1000 - deltaA*T*math.log(T) - deltaB*T**2/2 -  deltaC*T**3/6 - deltaD*T**4/12 - deltaE/(2*T) - 8.314 * T * I)/1000

Ka = math.exp(-22048/(8.314*T))
e = Symbol('e')
Eq3 = ( (e*(3-2*e)**2) / (1 - e)**3) / 0.4973 -1
e = round(solve(Eq3,e)[0],4)
yCH3OH = e/(3-2*e)
yco = (1-e)/(3-2*e)
yh2 = (2*(1-e) / (3-2*e))

# Results
print "E = %.4f"%e
print "YCH30H = %.4f"%yCH3OH
print "YCO = %.4f"%yco
print "YH2 = %.4f"%yh2
E = 0.0506
YCH30H = 0.0175
YCO = 0.3275
YH2 = 0.6550

Example 14.6, Page 496

In [23]:
from sympy.solvers import solve 
from sympy import Symbol 

#Variables
e = Symbol('e')

#Calculations
x = solve((((e/(3-2*e))/(((1-e)/(3-2*e))*(((2*(1-e))/(3-2*e))**2)))-49.73),e)
print "e =",round(x[0],3)
ych3oh = x[0]/(3-2*x[0])
yco = (1-x[0])/(3-2*x[0])
yh2 = ((2*(1-x[0]))/(3-2*x[0]))
print "yCH3OH =",round(ych3oh,4)
print "yCO =",round(yco,4)
print "yH2 =",round(yh2,4)
e = 0.801
yCH3OH = 0.5731
yCO = 0.1423
yH2 = 0.2846

Example 14.7 pageno : 497

In [2]:
from sympy import Symbol
from sympy.solvers import solve

# variables
Ka = 4.973*10**-3          # from example 14.5
Ky = 25 * Ka               # from example 14.5

# calculations
e = Symbol('e')
Eq1 = ( (e*(4-e)**2)/(1-e)**3 ) / .1243 - 1
e = round(solve(Eq1,e)[0],5)
yCH3OH = e/(8-2*e)
yco = (1-e)/(8-2*e)
yh2 = (2*(1-e) / (8-2*e))
yA = 5/(8-2*e)

# Results
print "E = %.5f"%e
print "YCH30H = %.5f"%yCH3OH
print "YCO = %.5f"%yco
print "YH2 = %.5f"%yh2
print "YA = %.5f"%yA
E = 0.00762
YCH30H = 0.00095
YCO = 0.12428
YH2 = 0.24857
YA = 0.62619

Example 14.8 pageno : 198

In [134]:
from sympy import Symbol
from sympy.solvers import solve

# variables
Ka = 4.973*10**-3           # example 14.5
Ky = 25 * Ka                # example 14.5

# calculations and results
e = Symbol('e')
Eq1= ((e*(1-e))/(1-2*e)**2)/.03108 - 1
e = solve(Eq1,e)
print "part a"
print "e = %.5f and e = %.5f"%(e[0],e[1])
print "The admissible value is e = %.5f"%e[0]

print "part b"
print "e = 0.0506 ( from example 14.5 )"

x = 4
e = Symbol('e')
Eq1=  ((e/(1-e) ) * ( (5 - 2*e)/ (2-e) )**2)/ 0.4973 - 1
e = solve(Eq1,e)
print "part c"
print "e = %.4f" %(e[0])
part a
e = 0.02845 and e = 0.97155
The admissible value is e = 0.02845
part b
e = 0.0506 ( from example 14.5 )
part c
e = 0.0727

Example 14.9 pageno : 499

In [138]:
from sympy import Symbol
from sympy.solvers import solve

# variables
Ka = 4.973*10**-3           # example 14.5
Ky = 25 * Ka                # example 14.5

# calculations
e = Symbol('e')
Eq1 = (0.02 + e)*(1.51-e)**2 / (1-e)**3 / Ky - 1
e = solve(Eq1,e)[0]

# Results
print "E = %.5f"%e
E = 0.03165

Example 14.10 pageno : 500

In [7]:
%matplotlib inline

from matplotlib.pyplot import *
from sympy import Symbol
from sympy.solvers import solve
from scipy.integrate import quad
import math

# variables
deltaA = round(28.850 - 27.012 - 0.5 * 30.255,2)
deltaB = (12.055 - 3.509 - 0.5 * 4.207) * 10**-3
deltaC = 0
deltaD = 0
deltaE = round((1.006 - 0.690 + 0.5 * 1.887),3) * 10**5
deltaH298 =  -241.997
Ts = [2000,2500,3000,3500,3800]
deltaHt = []
deltaGt = []
Ka = []
Ea = []
Eb = []

# Calculations
T = 298.15
deltaH0 = Symbol('deltaH0')
Eq = deltaH0 + deltaA*T + deltaB*T**2/2 + deltaC*T**3/3 + deltaD*T**4/4 - deltaE/T - deltaH298*1000
deltaH0 = round(solve(Eq,deltaH0)[0]/1000,3)
I = Symbol('I')
Eq1 = deltaH0*1000 - deltaA*T*math.log(T) - deltaB*T**2/2 -  deltaC*T**3/6 - deltaD*T**4/12 - deltaE/(2*T) - 8.314 * T * I + 228600
I = round(solve(Eq1,I)[0],3)

def fun1(T1):
    #return 42.1395*(T1 - 298.15) + 5.613*10**-3/2*(T1**2 - 298.15**2) + .2535 * 10**5 * (1./T1 - 1/298.15)
    return 42.1395*T1 + 5.613*10**-3/2 * T1**2 + 0.2535*10**5/T1 - 12898.37

for T in Ts:
    Ht = deltaH0*1000 + deltaA*T + deltaB*T**2/2 + deltaC*T**3/3 + deltaD*T**4/4 - deltaE/T 
    Gt = deltaH0*1000 - deltaA*T*math.log(T) - deltaB*T**2/2 -  deltaC*T**3/6 - deltaD*T**4/12 - deltaE/(2*T) - 8.314 * T * I 
    ka = math.exp(-Gt/(8.314*T))
    e = Symbol('e')
    Eq2 = ( e * (3-e)**(1./2) ) / ( 1-e)**(3./2) / ka - 1
    b = round(solve(Eq2)[0],4)
    a = (quad(fun1,298.15,T)[0]/1000) / -Ht
    deltaHt.append(Ht)
    deltaGt.append(Gt)
    Ka.append(ka)
    Ea.append(a)
    Eb.append(b)


# Results    
plot(Ts,Ea,"g")
plot(Ts,Eb,"b")
xlabel("T(k)")
ylabel("E")
suptitle("Plot of e versus adiabatic reaction temperature")
show()
print "from plot, it can be seen that both lines are simultaneously satisfied at the point \
\nintersection where e = 0.68 and T = 3440 K"
e = 0.68
T = 3440
yh2 = (1- e)/(1.5 - 0.5*e)
yo2 = 0.5*(1-e) / (1.5-0.5*e)
yh2o = e/(1.5- 0.5*e)

print "yH2 = %.4f"%yh2
print "yO2 = %.4f"%yo2
print "yH2O = %.4f"%yh2o
Populating the interactive namespace from numpy and matplotlib
WARNING: pylab import has clobbered these variables: ['solve', 'draw_if_interactive', 'e']
`%pylab --no-import-all` prevents importing * from pylab and numpy
from plot, it can be seen that both lines are simultaneously satisfied at the point 
intersection where e = 0.68 and T = 3440 K
yH2 = 0.2759
yO2 = 0.1379
yH2O = 0.5862

Example 14.11 Page No : 506

In [169]:
from numpy import *

# Variables
stoichio_matrix = matrix([[-1,-1, 1, 3, 0, 0],[0, -1, -1, 1, 1 ,0],[-1, -2, 0 ,4, 1, 0], \
                [0,0, 1, 0, -1, 0.5],[-1, 2 ,0, 0, 1, -2],[-1, 1, 1, 1, 0, -1]])			 #Framing the stoichiometric coefficient matrix


# Calculations

r = linalg.matrix_rank(stoichio_matrix)

stoichio_matrix[0] = -stoichio_matrix[0];
stoichio_matrix[2] = stoichio_matrix[2]+stoichio_matrix[0];
stoichio_matrix[5:] = stoichio_matrix[5:]+stoichio_matrix[0];
stoichio_matrix[6:] = stoichio_matrix[6:]+stoichio_matrix[0];
stoichio_matrix[1] = -stoichio_matrix[1];
stoichio_matrix[2] = stoichio_matrix[2]+stoichio_matrix[1];
stoichio_matrix[5:] = stoichio_matrix[5:]-(3*stoichio_matrix[1]);
stoichio_matrix[6:] = stoichio_matrix[6:]-(2*stoichio_matrix[1]);
x = stoichio_matrix[:3]
y = stoichio_matrix[:4];
stoichio_matrix[:4] = y;
stoichio_matrix[:3] = x;
stoichio_matrix[5:] = stoichio_matrix[5:]+(4*stoichio_matrix[3]);
stoichio_matrix[6:] = stoichio_matrix[6:]+(2*stoichio_matrix[3]);



# Results
print ' The stoichiometric coefficient matrix after performing the elementary row operations = ';
print (stoichio_matrix);
print ' The number of primary reactions = %d'%(r);
print ' The non zero rows are 1,2,4'
print ' The primary reactions are: CH4g)+H2Og)--->COg)+3H2g)%( COg)+H2Og)--->CO2g)+H2g),( CO2g)--->COg)+1./2)O2g)'
 The stoichiometric coefficient matrix after performing the elementary row operations = 
[[ 1.   1.  -1.  -3.  -0.  -0. ]
 [-0.   1.   1.  -1.  -1.  -0. ]
 [ 0.   0.   0.   0.   0.   0. ]
 [ 0.   0.   1.   0.  -1.   0.5]
 [-1.   2.   0.   0.   1.  -2. ]
 [ 0.  -1.   1.   1.  -1.   1. ]]
 The number of primary reactions = 3
 The non zero rows are 1,2,4
 The primary reactions are: CH4g)+H2Og)--->COg)+3H2g)%( COg)+H2Og)--->CO2g)+H2g),( CO2g)--->COg)+1./2)O2g)

Example 14.13 pageno : 510

In [214]:
from sympy import Symbol
from sympy.solvers import solve

# variables
Ka1 = 0.1429
Ka2 = 2

# calculations and Results

# assume e1 = .3
e2 = Symbol('e2')
Eq1 = ((0.3 - e2)*(0.3+e2) /((0.7 - e2) * 0.7)) / 0.1429 -1
e2 =  solve(Eq1)
print "E2 = %.1f and %.1f "%(e2[0],e2[1])

# since negative value is inadmissible. e2 = 0.2
e2 = round(e2[1],1)
e1 = Symbol('e1')
Eq2 = ((e1 + 0.2)*(0.2) /((0.8 - e1) *(e1 - 0.2) ) ) / 2 -1
e1 =  solve(Eq2)
print "E1 = %.1f and %.1f "%(e1[0],e1[1])
e1 = round(e1[0],1)
print "Satisfied results when e1 = %.1f and e2 = %.1f"%(e1,e2)

yA = (1-e1-e2)/2
yB = (1-e1)/2
yC = (e1-e2)/2
yD = (e1+e2)/2
yE = e2/2

print "yA = %.2f"%yA
print "yB = %.2f"%yB
print "yC = %.2f"%yC
print "yD = %.2f"%yD
print "yE = %.2f"%yE
E2 = -0.1 and 0.2 
E1 = 0.3 and 0.6 
Satisfied results when e1 = 0.3 and e2 = 0.2
yA = 0.25
yB = 0.35
yC = 0.05
yD = 0.25
yE = 0.10

Example 14.14 page no : 515

In [220]:
from sympy import Symbol
from sympy.solvers import solve

# variables
Kc = 2.92

# Calculations and results
e = Symbol('e')
Eq1 = (e*(10+e)/((5-e)*(10-e)))/Kc - 1
e = round(solve(Eq1)[0],4)
print "E = %.4f"%e

cA = 5-e
cB = 10-e
cC = e
cD = 10+e

print "CA = %.4f kmol/m*3"%cA
print "CB = %.4f kmol/m*3"%cB
print "CC = %.4f kmol/m*3"%cC
print "CD = %.4f kmol/m*3"%cD
E = 3.0446
CA = 1.9554 kmol/m*3
CB = 6.9554 kmol/m*3
CC = 3.0446 kmol/m*3
CD = 13.0446 kmol/m*3

Example 14.15, Page 517

In [6]:
import math
# Variables
T = 1200.   			 #temperature in K
T0 = 298.15 			 #reference temperature in K
a = [41.84,45.369,82.34];
b = [20.25*10**-3,8.688*10**-3,49.75*10**-3];
c = [0,0,0];
d = [0,0,0];
e = [-4.51*10**5,-9.619*10**5,-12.87*10**5];
del_Gf = [-604.574,-394.815,-1129.515]
del_Hf = [-635.975,-393.978,-1207.683]
n = [1,1,-1]
R = 8.314

# Calculations
del_G = (n[0]*del_Gf[0])+(n[1]*del_Gf[1])+(n[2]*del_Gf[2]);
del_H = (n[0]*del_Hf[0])+(n[1]*del_Hf[1])+(n[2]*del_Hf[2]);
del_a = (n[0]*a[0])+(n[1]*a[1])+(n[2]*a[2]);
del_b = (n[0]*b[0])+(n[1]*b[1])+(n[2]*b[2]);
del_c = (n[0]*c[0])+(n[1]*c[1])+(n[2]*c[2]);
del_d = (n[0]*d[0])+(n[1]*d[1])+(n[2]*d[2]);
del_e = (n[0]*e[0])+(n[1]*e[1])+(n[2]*e[2]);
del_H0 = ((del_H*10**3)-((del_a*T0)+((del_b/2)*T0**2)+((del_c/3)*T0**3)+((del_d/4)*T0**4)-(del_e/T0)))*10**-3;
IR = (1./(T0))*((del_H0*10**3)-(del_a*T0*math.log(T0))-((del_b/2)*T0**2)-((del_c/6)*T0**3)-((del_d/12)*T0**4)-((del_e/2)*(1./T0))-(del_G*10**3));
del_G_T = ((del_H0*10**3)-(del_a*T*math.log(T))-((del_b/2)*T**2)-((del_c/6)*T**3)-((del_d/12)*T**4)-((del_e/2)*(1./T))-(IR*T))*10**-3;
Ka = math.exp((-del_G_T*10**3)/(R*T))
y = 1;
phi = 1;
f0 = 1;
P = (Ka*f0)/(phi*y)

# Results
print ' The decomposition pressure,P = %f bar '%(P);
 The decomposition pressure,P = 2.384295 bar