Chapter 9 : Chemical Equilibria

Example 9.1.A page : 359

In [2]:
			
# Variables
g11 = 178900 			#kJ/kmol
g12 = 207037 			#kJ/kmol
g21 = 211852 			#kJ/kmol
g22 = 228097 			#kJ/kmol
			
# Calculations
dG = g21-g11
			
# Results
print "Standard free energy change  =  %d kJ/kmol"%(dG)

Standard free energy change  =  32952 kJ/kmol

Example 9.1B page : 359

In [3]:
import math 
			
# Variables
m1 = 54.1
m2 = 56.1
m3 = 2.
cp1 = 2.122 			#kJ/kmol K
cp2 = 2.213 			#kJ/kmol K
cp3 = 14.499 			#kJ/kmol K
hf1 = 110200. 			#kJ/kmol
hf2 = -126. 			#kJ/kmol
T = 700.    			#K
Ts = 298. 	    		#K
			
# Calculations
hf = hf1-hf2
cpn = cp1*m1-cp2*m2+cp3*m3
h700 = hf+ cpn*(T-Ts)
s298 = 103.7
s700  =  s298 + cpn*math.log(T/Ts)
G700 = h700-T*s700
			
# Results
print "Change in gibbs energy  =  %d kJ/kmol"%(G700)
print ("The answer is a bit different due to rounding off error in textbook")

Change in gibbs energy  =  33888 kJ/kmol
The answer is a bit different due to rounding off error in textbook

Example 9.2 page : 361

In [4]:
import math 
			
# Variables
g1 = 150670. 			#kJ/kmol
g2 = 71500. 			#kJ/kmol
R = 8.314
Ts = 298. 			#K
T = 700. 			#K

#calculation
G = g1-g2
G2 = 33875 			#kJ/kmol
K1 = math.exp(-G/R/Ts)
K2 = math.exp(-G2/R/T)
			
# Results
print "In case 1, equilibrium constant  =  %.2e"%(K1)
print " In case 2, equilibrium constant  =  %.5f"%(K2)

In case 1, equilibrium constant  =  1.33e-14
 In case 2, equilibrium constant  =  0.00297

Example 9.3 page : 363

In [5]:
import math 
from scipy.integrate import quad
			
# Variables
R = 8.3143
T1 = 1273 			#K
T2 = 2273 			#K
k2 = 0.0018
A = 123.94
B = 7.554
C = 8.552*10**-3
D = -13.25e-6
E = 7.002e-9
F = 13.494e-13
			
# Calculations
def  cp(T):
    return A/T**2 +B/T +C +D*T +E*T**2 -F*T**3


lnk = 1/R * quad(cp,T1,T2)[0]
k1 = k2/ math.exp(lnk)
			
# Results
print "Equilibrium constant  =  %.5f "%(k1)

Equilibrium constant  =  0.00112

Example 9.4 page : 368

In [2]:
import math
from scipy.optimize import fsolve
			
# Variables
G = -30050. 			#kJ/kmol
R = 8.314
T = 573.     			#K
			
# Calculations
lnk = G/(R*T)
k = math.exp(lnk)

def f(x):
    return 4*x**2 - k*(1-x)**2

vec = fsolve(f,0,full_output=1)

x2 = vec[0]
# Results
print "Mole fraction of HCN  =  %.4f"%(x2)

Mole fraction of HCN  =  0.0209

Example 9.4B page : 368

In [3]:
import math 
from numpy import roots
			
# Variables
G = -30050. 			#kJ/kmol
R = 8.314
T = 573. 			#K
phi1 = 0.980
phi2 = 0.915
phi3 = 0.555
			
# Calculations
lnk = G/(R*T)
k = math.exp(lnk)
kexp =  k*phi1*phi2/phi3**2 /4
vec = roots([1-kexp,2*kexp,-kexp])
x2 = vec[1]
			
# Results
print "Mole fraction of HCN  =  %.4f"%(x2)

Mole fraction of HCN  =  0.0351

Example 9.5 page : 370

In [4]:
import math 
from scipy.optimize import fsolve
			
# Variables
kp = 74.
kp2 = kp**2			
# Calculations
def fun(f):
    return f**2 *(100-6*f) - kp**2 *(1-f)**2 *(9-6*f)

vec = fsolve(fun,0,full_output=1)
fn = vec[0]
			
# Results
print "Fractional conversion  =  %.3f"%(fn)

Fractional conversion  =  0.934

Example 9.6 page : 372

In [10]:
			
# Variables
C = 3.
phi = 3.
R = 1.
Sc = 0.

def fun(C,phi,R,Sc):
    return 2+C-phi-R-Sc
			
# Calculations
V = fun(C,phi,R,Sc)
			
# Results
print "Degrees of freedom  =  %d "%(V)

Degrees of freedom  =  1

Example 9.6B page : 373

In [11]:
			
# Variables
C = 3.
phi = 1.
R = 1.
Sc = 1.

def fun(C,phi,R,Sc):
    return 2+C-phi-R-Sc

# Calculations
V = fun(C,phi,R,Sc)
			
# Results
print "Degrees of freedom  =  %d "%(V)

Degrees of freedom  =  2

Example 9.6C page : 373

In [12]:
			
# Variables
C = 6.
phi = 1.
R = 3.
Sc = 0.

def fun(C,phi,R,Sc):
    return 2+C-phi-R-Sc

			
# Calculations
V = fun(C,phi,R,Sc)
			
# Results
print "Degrees of freedom  =  %d "%(V)

Degrees of freedom  =  4

Example 9.7 page : 377

In [13]:
			
# Variables
a1 = 0.956
y = 0.014
x = 0.956
M = 18.
z = 0.475
P = 8.37 			#Mpa
			
# Calculations
m = y/(x*M) *10**3
w = 0.0856
phi1 = -0.04
phi2 = 0.06
phi = 10**(phi1+ w*phi2)
f = z*phi*P
K = m/(f*a1)
			
# Results
print "Equilibrium constant  =  %.3f"%(K)

Equilibrium constant  =  0.232

Example 9.9 page : 385

In [14]:
from numpy import array
			
# Variables
y = 0.18
z = 0.6
			
# Calculations
mole = array([1-y-z, 5-y-2*z, y, 3*y+4*z, z])
s = sum( mole)
molef = mole/s
			
# Results
print "Product composition moles   =  ",
print (mole)
print "Mole fraction  =  ",
print (molef)

Product composition moles   =   [ 0.22  3.62  0.18  2.94  0.6 ]
Mole fraction  =   [ 0.02910053  0.47883598  0.02380952  0.38888889  0.07936508]

Example 9.10 page : 388

In [5]:
from sympy import Symbol,solve
			
# Variables
kp = 1.09
			
# Calculations
x = Symbol('x')
vec = solve(kp/4**4 /4 *(1-x)*(5-2*x)**2 *(6+2*x)**2 -x**5)
x = vec[0]
pro = [1-x, 5-2*x, x, 4*x, 0]
			
# Results
print "Equlibrium composition (moles) =  ",
print (pro)

Equlibrium composition (moles) =   [0.273026192093833, 3.54605238418767, 0.726973807906167, 2.90789523162467, 0]

Example 9.10B page : 390

In [6]:
import math 
from numpy import array
from scipy.optimize import fsolve
			
# Variables
kp = 1.09
kp2 = 0.154
feed = array([ 1, 5, 0, 0, 0 ])
			
# Calculations

def f1(x):
    return kp/4**4 /4 *(1-x)*(5-2*x)**2 *(6+2*x)**2 -x**5

vec = fsolve(f1,0,full_output=1)[0]
x = vec[0]
pro = feed - array([x, 2*x, -x, -4*x, 0])

def f2(y):
    return kp2*(0.273-y)*(0.727-y)*(7.454+2*y)**2 - 4*y**2 *(2.908+2*y)**2 *4

vec2 = fsolve(f2,0,full_output=1)[0]    
y = vec2[0]
pro2 = pro- array([ y, 0, y, -2*y, -2*y])

def f3(z):
    return kp*(0.189-z)*(3.546-2*z)**2 *(7.622+2*z)**2 -(0.643+z)*(3.076+4*z)**4 *4
    
vec3 =  fsolve(f3,0,full_output=1)[0]
z = vec3[0]
pro3 = pro2 - array([z, 2*z, -z, -4*z, 0])

def f4(w):
    return kp2*(0.229-w)*(0.603-w)*(7.542+2*w) - (2.916+2*w)**2 *(0.168+2*w)**2 *4

vec4 = fsolve(f4,0,full_output=1)[0]
w = vec4[0]
w = 0.01
pro4 = pro3 - array([w, 0, w, -2*w, -2*w])
			
# Results
print "feed  =  ",
print (feed)
print "After reactor 1,",
print (pro)
print "After reactor 2,",
print (pro2)
print "After reactor 3,",
print (pro3)
print "After reactor 4",
print (pro4)
print ("The answers are a bit different due to rounding off error in textbook")

feed  =   [1 5 0 0 0]
After reactor 1, [ 0.27302619  3.54605238  0.72697381  2.90789523  0.        ]
After reactor 2, [ 0.18846228  3.54605238  0.6424099   3.07702305  0.16912782]
After reactor 3, [ 0.22284245  3.61481272  0.60802973  2.93950238  0.16912782]
After reactor 4 [ 0.21284245  3.61481272  0.59802973  2.95950238  0.18912782]
The answers are a bit different due to rounding off error in textbook