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);
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);
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);
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);
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
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)
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
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])
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
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
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)'
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
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
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);