# Chapter 17 : Chemical Reactions Equilibria¶

### Example 17.2 Page Number : 598¶

In [1]:

from scipy.optimize import fsolve
import math

# Variables
P = 1;			#[atm] - Reactor pressure
T = 749;			#[K] - Reactor temperature
K = 74;			# Equlibrium constant

# SO2 + (1/2)*O2 - SO3

# Calculations and Results
Kp = P**(1);
Ky = K/Kp;

#(1)
# Initial number of moles of the components are
n_SO2_1_in = 12;
n_O2_1_in = 9;
n_SO3_1_in = 0;

# Let the reaction coordinate at equilibrium for the reaction be X
# At equilibrium, the moles of the components be
# n_SO2_1_eq = 12 - X
# n_O2_1_eq = 9 - 0.5*X
# n_SO3_1_eq = X
# Total moles = 21 - 0.5*X

# The mole fractions of the components at equilibrium are
# y_SO3 = X/(21-0.5*X)
# y_SO2 = (12-X)/(21-0.5*X)
# y_O2 = (9-0.5*X)/(21-0.5*X)

# Ky = y_SO3/(y_SO2*y_O2**(2))
# Ky = (X*(21-0.5*X)**(1/2))/((12-X)*(9-0.5*X)**(1/2))
def f(X):
return  Ky-(X*(21-0.5*X)**(1./2))/((12-X)*(9-0.5*X)**(1./2))
X_1 = fsolve(f,11)

y_SO3_1 = X_1/(21-0.5*X_1);
y_SO2_1 = (12-X_1)/(21-0.5*X_1);
y_O2_1 = (9-0.5*X_1)/(21-0.5*X_1);

print " 1).The moles of SO3 formed = %f mol"%(X_1);
print "     The mole fractions at equilibrium are y_S03 = %f y_SO2 = %f and y_O2 = %f"%(y_SO3_1,y_SO2_1,y_O2_1);

#(2)
# Initial number of moles of the components are
n_SO2_2_in = 24;
n_O2_2_in = 18;
n_SO3_2_in = 0;

# At equilibrium, the moles of the components be
# n_SO2_1_eq = 24 - X
# n_O2_1_eq = 18 - 0.5*X
# n_SO3_1_eq = X
# Total moles = 42 - 0.5*X

# The mole fractions of the components at equilibrium are
# y_SO3 = X/(42-0.5*X)
# y_SO2 = (24-X)/(42-0.5*X)
# y_O2 = (18-0.5*X)/(42-0.5*X)

# Ky = y_SO3/(y_SO2*y_O2**(2))
# Ky = (X*(42-0.5*X)**(1/2))/((24-X)*(18-0.5*X)**(1/2))
def f1(X):
return  Ky-(X*(42-0.5*X)**(1./2))/((24-X)*(18-0.5*X)**(1./2))
X_2 = fsolve(f1,22)

y_SO3_2 = X_2/(42-0.5*X_2);
y_SO2_2 = (24-X_2)/(42-0.5*X_2);
y_O2_2 = (18-0.5*X_2)/(42-0.5*X_2);
print " 2).The moles of SO3 formed = %f mol"%(X_2);
print "     The mole fractions at equilibrium are y_S03 = %f, y_SO2 = %f and y_O2 = %f"%(y_SO3_2,y_SO2_2,y_O2_2);

#(3)
# Initial number of moles of the components are
n_SO2_3_in = 12;
n_O2_3_in = 9;
n_SO3_3_in = 0;
n_N2 = 79;

# At equilibrium, the moles of the components be
# n_SO2_1_eq = 12 - X
# n_O2_1_eq = 9 - 0.5*X
# n_SO3_1_eq = X
# Total moles = 100 - 0.5*X

# The mole fractions of the components at equilibrium are
# y_SO3 = X/(100-0.5*X)
# y_SO2 = (12-X)/(100-0.5*X)
# y_O2 = (9-0.5*X)/(100-0.5*X)

# Ky = y_SO3/(y_SO2*y_O2**(2))
# Ky = (X*(100-0.5*X)**(1/2))/((12-X)*(9-0.5*X)**(1/2))
def f2(X):
return  Ky-(X*(100-0.5*X)**(1./2))/((12-X)*(9-0.5*X)**(1./2))
X_3 = fsolve(f2,10)

y_SO3_3 = X_3/(100-0.5*X_3);
y_SO2_3 = (12-X_3)/(100-0.5*X_3);
y_O2_3 = (9-0.5*X_3)/(100-0.5*X_3);

print " 3).The moles of SO3 formed = %f mol"%(X_3);
print "     The mole fractions at equilibrium are y_S03 = %f y_SO2 = %f and y_O2 = %f"%(y_SO3_3,y_SO2_3,y_O2_3);

 1).The moles of SO3 formed = 11.655537 mol
The mole fractions at equilibrium are y_S03 = 0.768215 y_SO2 = 0.022704 and y_O2 = 0.209081
2).The moles of SO3 formed = 23.311074 mol
The mole fractions at equilibrium are y_S03 = 0.768215, y_SO2 = 0.022704 and y_O2 = 0.209081
3).The moles of SO3 formed = 11.202213 mol
The mole fractions at equilibrium are y_S03 = 0.118669 y_SO2 = 0.008451 and y_O2 = 0.036006


### Example 17.3 Page Number : 599¶

In [2]:

# Variables
T = 600;			#[K] -  Reactor temperature
P = 300;			#[atm] - Reactor pressure
K = 0.91*10**(-4);			# Equilibrium constant

# The fugacity coefficients of the components are
phi_CO = 1.0;
phi_H2 = 1.2;
phi_CH3OH = 0.47;

# CO + 2*H2 - CH3OH

# For gas phase reactions the standard state is pure ideal gas and thus fi_0 = 1 atm and thus
# ai_cap = fi_cap/fi_0 = yi*P*phi_i_cap/1
# Thus K = Ky*Kp*K_phi
# Calculations and Results
Kp = P**(1-3);
K_phi = phi_CH3OH/(phi_CO*phi_H2**(2));
Ky = K/(Kp*K_phi);

# Let the reaction coordinate at equilibrium for the reaction be X
# At equilibrium ,the moles of the components be
# n_CO = 1 - X
# n_H2 = 3 - 2*X
# n_CH3OH = X
# Total moles = 4 - 2*X

# The mole fractions of the components at equilibrium are
# y_CO = (1-X)/(4-2*X)
# y_H2 = (3-2*X)/(4-2*X)
# y_CH3OH = (X)/(4-2*X)

# Ky = y_CH3OH/(y_CO*y_H2**(2)) = (X/(4-2*X))/(((1-X)/(4-2*X))*((3-2*X)/(4-2*X))**(2))
def f(X):
return Ky-(X/(4-2*X))/(((1-X)/(4-2*X))*((3-2*X)/(4-2*X))**(2))
X = fsolve(f,0.1)

# Therefore at equilibrium
y_CO = (1-X)/(4-2*X);
y_H2 = (3-2*X)/(4-2*X);
y_CH3OH = (X)/(4-2*X);

print " The mole fractions at equilibrium are y_CO = %f y_H2 = %f and y_CH3OH = %f"%(y_CO,y_H2,y_CH3OH);

 The mole fractions at equilibrium are y_CO = 0.051857 y_H2 = 0.551857 and y_CH3OH = 0.396286


### Example 17.4 Page Number : 600¶

In [3]:

# Variables
T = 600;			#[K] -  Reactor temperature
P = 4;			#[atm] - Reactor pressure
K = 1.175;			# Equilibrium constant

# (1/2)*N2 + (3/2)*H_2 - NH3

# Initial number of moles of the components are
n_N2 = 1;
n_H2 = 3;
n_HN3 = 0;

# Calculations and Results
# Let the reaction coordinate at equilibrium for the reaction be X.
# At equilibrium ,the moles of the components be
# n_N2 = 1 - 0.5*X
# n_H2 = 3 - 1.5*X
# n_NH3 = X
# Total moles = 4 - X

# We have, K = Ky*Kp
Kp = P**(1-2);			#[atm**(-1)]
Ky = K/(Kp);

# Ky = y_NH3/(y_N2**(1/2)*y_H2**(3/2)) = (X/(4-X))/(((1-0.5*X)/(4-X))**(1/2)*((3-1.5*X)/(4-X))**(3/2))
# Solving the above equation we get
def f(X):
return Ky - (X/(4-X))/(((1-0.5*X)/(4-X))**(1./2)*((3-1.5*X)/(4-X))**(3./2))
X = fsolve(f,0.1)

y_NH3 = X/(4-X);			# Mole fraction of NH3 at equilibrium

print " The value of Kp = %f and Ky = %f "%(Kp,Ky);
print " The mole fractions of NH3 at equilibrium is %f"%(y_NH3);

# We know that for ideal gas, P*V = n*R*T and thus P is directly proportional to n at constant V and T.
# Let P = k*n
# Initially P = 4 atm and n = 4 moles, thus K = 1 and we get p = n, where P is in atm.
# Thus at equilibrium P = 4 - X

# Ky = K/Kp = 1.175*P = 1.175*(4 - X)
# (X/(4-X))/(((1-0.5*X)/(4-X))**(1/2)*((3-1.5*X)/(4-X))**(3/2)) = 1.175*(4 - X)
# Solving the above equation we get
def f1(X):
return (X/(4-X))/(((1-0.5*X)/(4-X))**(1./2)*((3-1.5*X)/(4-X))**(3./2))-1.175*(4-X)
X_prime = fsolve(f1,1)

# Therefore at equilibrium
P_prime = 4 - X_prime;
y_NH3_prime = X_prime/(4-X_prime);

print " If reaction is carried out at constant temperature and volumethen"
print " The equilibrium pressure is %f atm"%(P_prime);
print " The equilibrium mole fractions of NH3 in the reactor is %f"%(y_NH3_prime);

 The value of Kp = 0.250000 and Ky = 4.700000
The mole fractions of NH3 at equilibrium is 0.454388
If reaction is carried out at constant temperature and volumethen
The equilibrium pressure is 2.863057 atm
The equilibrium mole fractions of NH3 in the reactor is 0.397108


### Example 17.5 Page Number : 601¶

In [4]:

# Variables
T = 400;			#[K] -  Reactor temperature
P = 1;			#[atm] - Reactor pressure
K = 1.52;			# Equilibrium constant
y_H2 = 0.4;			# Equilibrium mole fraction of hydrogen

# Calculations
# CO(g) + 2*H_2(g) - CH3OH(g)

# K = y_CH3OH/(y_CO*y_H2**(2)*P**(2))
# Let total number of moles at equilibrium be 1
# y_CH3OH = 0.6 - y_CO;
# (0.6 - y_CO)/y_CO = K*P**(2)*y_H2**(2)

y_CO = 0.6/(1 + K*P**(2)*y_H2**(2));
y_CH3OH = 0.6 - y_CO;

# Results
print " The mole fractions are y_CO = %f and y_CH3OH = %f "%(y_CO,y_CH3OH);

 The mole fractions are y_CO = 0.482625 and y_CH3OH = 0.117375


### Example 17.6 Page Number : 602¶

In [5]:


# Variables
T = 749.;			#[K] - Reactor temperature
P = 1.;			#[atm] - Reactor pressure
K = 74.;

# Calculations and Results
Kp = P**(-1./2);			#[atm**(-1/2)]
Ky = K/Kp;

# SO2 + (1/2)*O2 - SO3

# Initial number of moles of the components are
n_SO2_1_in = 10;
n_O2_1_in = 8;
n_SO3_1_in = 0;

# Let the reaction coordinate at equilibrium for the reaction be X
# At equilibrium, the moles of the components be
# n_SO2_1_eq = 10 - X
# n_O2_1_eq = 8 - 0.5*X
# SO3_1_eq = X
# Total moles = 18 - 0.5*X

# The mole fractions of the components at equilibrium are
# y_SO3 = X/(18-0.5*X)
# y_SO2 = (10-X)/(18-0.5*X)
# y_O2 = (8-0.5*X)/(18-0.5*X)

# Ky = y_SO3/(y_SO2*y_O2**(2))
# Ky = (X*(18-0.5*X)**(1/2))/((10-X)*(8-0.5*X)**(1/2))
def f(X):
return  Ky-(X*(18-0.5*X)**(1./2))/((10-X)*(8-0.5*X)**(1./2))
X_1 = fsolve(f,11)

n_SO3 = X_1;
n_SO2 = 10 - X_1;
n_O2 = 8 - 0.5*X_1;

print " 1).The moles of the components at equilibrium are n_SO3 = %f mol n_SO2 = %f mol and n_O2 = %f mol"%(n_SO3,n_SO2,n_O2);

# Now for the reaction
# 2*SO2 + O2 - 2*SO3

# The equilibrium constant for this reaction is KP**(2)
Ky_prime = Ky**(2);

# At equilibrium, the moles of the components be
# n_SO2_1_eq = 10 - 2*X
# n_O2_1_eq = 8 - X
# SO3_1_eq = 2*X
# Total moles = 18 - X

# The mole fractions of the components at equilibrium are
# y_SO3 = 2*X/(18-X)
# y_SO2 = (10-2*X)/(18-X)
# y_O2 = (8- X)/(18-X)

# Ky_prime = y_SO3**(2)/(y_SO2**(2)*y_O2)
# Ky_prime = ((2*X)**(2)*(18-X))/((10-2*X)**(2)*(8-X))
def f1(X):
return  Ky_prime-((2*X)**(2)*(18-X))/(((10-2*X)**(2))*(8-X))
X_2 = fsolve(f1,6)

n_SO3_prime = 2*X_2;
n_SO2_prime = 10 - 2*X_2;
n_O2_prime = 8 - X_2;

print " 2).The moles of the components at equilibrium are n_SO3 = %f mol n_SO2 = %f mol and n_O2 = %f mol"%(n_SO3_prime,n_SO2_prime,n_O2_prime);
print "     Thus the number of moles remains the same irrespective of the stoichoimetry of the reaction"

 1).The moles of the components at equilibrium are n_SO3 = 9.730824 mol n_SO2 = 0.269176 mol and n_O2 = 3.134588 mol
2).The moles of the components at equilibrium are n_SO3 = 9.730824 mol n_SO2 = 0.269176 mol and n_O2 = 3.134588 mol
Thus the number of moles remains the same irrespective of the stoichoimetry of the reaction


### Example 17.7 Page Number : 603¶

In [7]:


# Variables
T = 500;			#[K]
# For the reaction,  0.5*A2 + 0.5*B2 - AB
delta_G = -4200;			#[J/mol]
R = 8.314;			#[J/mol*K] - Universal gas constant

#(1)
# A2 + B2 - 2*AB

# Calculations and Results
# We know delta_G_rkn_0 = -R*T*math.log(K)
delta_G_1 = 2*delta_G;
K_1 = math.exp(-delta_G_1/(R*T));			# Equilibrium constant at 500 K for the above reaction
# As can be seen the reaction is not affected by pressure and therefore K = Ky as Kp = 1
Ky = K_1;

# Initial number of moles of the components are
n_A2_1_in = 0.5;
n_B2_1_in = 0.5;
n_AB_1_in = 0;

# Let the reaction coordinate at equilibrium for the reaction be X
# At equilibrium, the moles of the components be
# n_A2_1_eq = 0.5 - X
# n_B2_1_eq = 0.5- X
# n_AB_1_eq = 2*X
# Total moles = 1

# Ky = (2*X)**(2)/(0.5-X)**(2)
def f(X):
return  Ky-(2*X)**(2)/(0.5-X)**(2)
X_1 = fsolve(f,0.2)

# The mole fractions of the components at equilibrium are
y_A2_1 = 0.5 - X_1;
y_B2_1 = 0.5- X_1;
y_AB_1 = 2*X_1;

print " 1).The mole fractions at equilibrium are y_A2 = %f y_B2 = %f and y_AB = %f"%(y_A2_1,y_B2_1,y_AB_1);

#(2)
# 0.5*A2 + 0.5*B2 - AB

# We know delta_G_rkn_0 = -R*T*math.log(K)
delta_G_2 = delta_G;
K_2 = math.exp(-delta_G_2/(R*T));			# Equilibrium constant at 500 K for the above reaction

# As can be seen the reaction is not affected by pressure and therefore K = Ky as Kp = 1
Ky_2 = K_2;

# Initial number of moles of the components are
n_A2_2_in = 0.5;
n_B2_2_in = 0.5;
n_AB_2_in = 0;

# Let the reaction coordinate at equilibrium for the reaction be X
# At equilibrium, the moles of the components be
# n_A2_2_eq = 0.5 - 0.5*X
# n_B2_2_eq = 0.5- 0.5*X
# n_AB_2_eq = X
# Total moles = 1

# Ky = y_AB/(y_A2**(1/2)*y_B2**(1/2))
# Ky = X/(0.5 - 0.5*X)
X_2 = 0.5*Ky_2/(1+0.5*Ky_2);

# The mole fractions of the components at equilibrium are
y_A2_2 = 0.5 - 0.5*X_2;
y_B2_2 = 0.5- 0.5*X_2;
y_AB_2 = X_2;

print " 2).The mole fractions at equilibrium are y_A2 = %f y_B2 = %f and y_AB = %f"%(y_A2_2,y_B2_2,y_AB_2);

#(3)
# 2*AB -  A2 + B2

K_3 = 1/K_1;			# Equilibrium constant at 500 K for the above reaction
# As can be seen the reaction is not affected by pressure and therefore K = Ky as Kp = 1
Ky_3 = K_3;

# Initial number of moles of the components are
n_AB_3_in = 1;
n_A2_3_in = 0;
n_B2_3_in = 0;

# Let the reaction coordinate at equilibrium for the reaction be X
# At equilibrium, the moles of the components be
# n_AB_3_eq = 1 - X
# n_A2_3_eq = X/2
# n_B2_3_eq = X/2
# Total moles = 1

# Ky = (X/2)**(2)/(1-X)**(2)
def f1(X):
return  Ky_3-(X/2)**(2)/(1-X)**(2)
X_3 = fsolve(f1,0.4)

# The mole fractions of the components at equilibrium are
y_A2_3 = X_3/2;
y_B2_3 = X_3/2;
y_AB_3 = 1-X_3;

print " 3).The mole fractions at equilibrium are y_A2 = %f y_B2 = %f and y_AB = %f"%(y_A2_3,y_B2_3,y_AB_3);

 1).The mole fractions at equilibrium are y_A2 = 0.210680 y_B2 = 0.210680 and y_AB = 0.578641
2).The mole fractions at equilibrium are y_A2 = 0.210680 y_B2 = 0.210680 and y_AB = 0.578641
3).The mole fractions at equilibrium are y_A2 = 0.210680 y_B2 = 0.210680 and y_AB = 0.578641


### Example 17.9 Page Number : 606¶

In [8]:


from scipy.optimize import fsolve
import math
from scipy.integrate import quad

# Variables
R = 1.987;			#[cal/mol-K]

delta_H_SO2_298 = -70.96;			#[kcal/mol] - Enthalpy of formation of S02 at 298.15 K
delta_H_SO3_298 = -94.45;			#[kcal/mol] - Enthalpy of formation of S03 at 298.15 K
delta_G_SO2_298 = -71.79;			#[kcal/mol] - Gibbs free energy change for formation of SO2 at 298.15 K
delta_G_SO3_298 = -88.52;			#[kcal/mol] - Gibbs free energy change for formation of SO3 at 298.15 K

# Cp_0 = a + b*T + c*T**(2) + d*T**(3)

a_SO2 = 6.157;
a_SO3 = 3.918;
a_O2 = 6.085;
b_SO2 = 1.384*10**(-2);
b_SO3 = 3.483*10**(-2);
b_O2 = 0.3631*10**(-2);
c_SO2 = -0.9103*10**(-5);
c_SO3 = -2.675*10**(-5);
c_O2 = -0.1709*10**(-5);
d_SO2 = 2.057*10**(-9);
d_SO3 = 7.744*10**(-9);
d_O2 = 0.3133*10**(-9);

#(1)
T_1 = 298.15;			#[K]

# Calculations and Results
delta_H_rkn_298 = delta_H_SO3_298 - delta_H_SO2_298;			#[kcal]
delta_H_rkn_298 = delta_H_rkn_298*10**(3);			#[cal]
delta_G_rkn_298 = delta_G_SO3_298 - delta_G_SO2_298;			#[kcal]
delta_G_rkn_298 = delta_G_rkn_298*10**(3);			#[cal]

delta_a = a_SO3 - a_SO2 - (a_O2/2);
delta_b = b_SO3 - b_SO2 - (b_O2/2);
delta_c = c_SO3 - c_SO2 - (c_O2/2);
delta_d = d_SO3 - d_SO2 - (d_O2/2);

def f60(T):
return delta_a+(delta_b*T)+(delta_c*T**(2))+(delta_d*T**(3))

# delta_H_rkn_T = delta_H_rkn_298 +  quad(f60,T_1,T)[0]

# On simplification we get
# delta_H_rkn_T = -22630.14 - 5.2815*T + 0.9587*10**(-2)*T**(2) - 0.5598*10**(-5)*T**(3) + 1.3826*10**(-9)*T**(4)

print " 1.The math.expression for delta_H_rkn_T as a function of T is given by";
print "     delta_H_rkn_T = -22630.14 - 5.2815*T + 0.9587*10**-2*T**2 - 0.5598*10**-5*T**3 + 1.3826*10**-9*T**4";

#(2)

#def f61(T):
#	 return K_T/K_298) = integrate(delta_H_rkn_T/T**(2)

# R*math.log(K_T/K_298) =  quad(f61,T_1,T)[0]

# First let us calculate K_298.
# delta_G_rkn_T = - R*T*math.log(K)
K_298 = math.exp(-delta_G_rkn_298/(R*T_1));

# On substituting the values and simplifying we get the math.expression
# math.log(K) = 3.87 + 11380.10/T - 2.6580*math.log(T) + 0.4825*10**(-2)*T - 0.1409*10**(-5)*T**(2) + 0.2320*10**(-9)*T**(3)

print " 2.The math.expression for math.logK as a function of T is given by";
print "     math.logK = 3.87 + 11380.10/T - 2.6580*math.logT + 0.4825*10**-2*T - 0.1409*10**-5*T**2 + 0.2320*10**-9*T**3";

#(3)
P = 1;			#[atm]
T = 880;			#[K]
K = math.exp(3.87 + 11380.10/T - 2.6580*math.log(T) + 0.4825*10**(-2)*T - 0.1409*10**(-5)*T**(2) + 0.2320*10**(-9)*T**(3));
Kp = P**(-1./2);			#[atm**(-1/2)]
Ky = K/Kp;

# Let the reaction coordinate at equilibrium for the reaction be X
# At equilibrium, the moles of the components be
# n_SO2_eq = 1 - X
# n_O2_eq = 0.5- 0.5*X
# n_SO3_1_eq = X
# Total moles = 1.5-0.5*X

# Ky = (X*(1.5-0.5*X)**(1/2))/((1-X)*(0.5-0.5*X)**(1/2))
def f(X):
return  Ky - (X*(1.5-0.5*X)**(1./2))/((1-X)*(0.5-0.5*X)**(1./2))
X = fsolve(f,0.8)

# The mole fraction of SO3 at equilibrium is given by
y_SO3 = X/(1.5-0.5*X);

print " 3).The mole fraction of SO3 at equilibrium is given by y_SO3 = %f"%(y_SO3);

 1.The math.expression for delta_H_rkn_T as a function of T is given by
delta_H_rkn_T = -22630.14 - 5.2815*T + 0.9587*10**-2*T**2 - 0.5598*10**-5*T**3 + 1.3826*10**-9*T**4
2.The math.expression for math.logK as a function of T is given by
math.logK = 3.87 + 11380.10/T - 2.6580*math.logT + 0.4825*10**-2*T - 0.1409*10**-5*T**2 + 0.2320*10**-9*T**3
3).The mole fraction of SO3 at equilibrium is given by y_SO3 = 0.649167


### Example 17.10 Page Number : 609¶

In [9]:


from numpy import *

# Variables
# (1/2)*N2 + (1/2)*O2 - NO

R = 1.987;			#[cal/mol-K]

delta_H_NO_298 = 21.600;			#[kcal/mol] - Enthalpy of formation of S02 at 298.15 K
delta_G_NO_298 = 20.719;			#[kcal/mol] - Gibbs free energy change for formation of SO2 at 298.15 K

# Cp_0 = a + b*T + c*T**(2) + d*T**(3)

a_N2 = 6.157;
a_O2 = 6.085;
a_NO = 6.461;
b_N2 = -0.03753*10**(-2);
b_O2 = 0.3631*10**(-2);
b_NO = 0.2358*10**(-2);
c_N2 = 0.1930*10**(-5);
c_O2 = -0.1709*10**(-5);
c_NO = -0.07705*10**(-5);
d_N2 = -0.6861*10**(-9);
d_O2 = 0.3133*10**(-9);
d_NO = 0.08729*10**(-9);

#(1)
T_1 = 298.15;			#[K]

# Calculations and Results
delta_H_rkn_298 = delta_H_NO_298;			#[kcal]
delta_H_rkn_298 = delta_H_rkn_298*10**(3);			#[cal]
delta_G_rkn_298 = delta_G_NO_298;			#[kcal]
delta_G_rkn_298 = delta_G_rkn_298*10**(3);			#[cal]

delta_a = a_NO - (a_N2/2) - (a_O2/2);
delta_b = b_NO - (b_N2/2) - (b_O2/2);
delta_c = c_NO - (c_N2/2) - (c_O2/2);
delta_d = d_NO - (d_N2/2) - (d_O2/2);

def f49(T):
return delta_a+(delta_b*T)+(delta_c*T**(2))+(delta_d*T**(3))

# delta_H_rkn_T = delta_H_rkn_298 +  quad(f49,T_1,T)[0]

# On simplification we get
# delta_H_rkn_T = 21584.63 - 0.033*T + 0.0365*10**(-2)*T**(2) - 0.0293*10**(-5)*T**(3) + 0.0685*10**(-9)*T**(4)

print " The math.expression for delta_H_rkn_T as a function of T is given by";
print " delta_H_rkn_T = 21584.63 - 0.033*T + 0.0365*10**-2*T**2 - 0.0293*10**-5*T**3 + 0.0685*10**-9*T**4";

# Now let us calculate K_298 (at 298 K)
# We know delta_G_rkn_298 = -R*T*math.log(K_298)
K_298 = math.exp(-delta_G_rkn_298/(R*T_1));			# Equilibrium constant at 298.15 K

#def f50(T):
#	 return K_2/K_1) = integrate(delta_H_rkn_298/(R*T**(2))

# math.log(K_2/K_1) =  quad(f50,T_1,T)[0]

# On substituting the values and simplifying we get the math.expression
# math.log(K) = 1.5103 - 10862.92/T - 0.0166*math.log(T) + 1.84*10**(-4)*T - 7.35*10**(-8)*T**(2) + 1.15*10**(-11)*T**(3)

print " The math.expression for math.logK as a function of T is given by";
print " math.logK = 1.5103 - 10862.92/T - 0.0166*math.logT + 1.84*10**-4*T - 7.35*10**-8*T**2 + 1.15*10**-11*T**3"

T = [500,1000,1500,2000,2500];
K = zeros(5);

print " T K  \t\t\t  K  ";

for i in range(5):
K[i] = math.exp(1.5103-10862.92/T[i] - 0.0166*math.log(T[i]) + 1.84*10**(-4)*T[i] - 7.35*10**(-8)*T[i]**(2) + 1.15*10**(-11)*T[i]**(3));

print " %f  \t\t  %e  "%(T[i],K[i]);

# It can be seen that for endothermic reactions equilibrium constant increases with temperature.
print " It can be seen that for endothermic reactions equilibrium constant increases with temperature";

 The math.expression for delta_H_rkn_T as a function of T is given by
delta_H_rkn_T = 21584.63 - 0.033*T + 0.0365*10**-2*T**2 - 0.0293*10**-5*T**3 + 0.0685*10**-9*T**4
The math.expression for math.logK as a function of T is given by
math.logK = 1.5103 - 10862.92/T - 0.0166*math.logT + 1.84*10**-4*T - 7.35*10**-8*T**2 + 1.15*10**-11*T**3
T K  			  K
500.000000  		  1.615470e-09
1000.000000  		  8.737610e-05
1500.000000  		  3.333913e-03
2000.000000  		  2.062328e-02
2500.000000  		  6.176400e-02
It can be seen that for endothermic reactions equilibrium constant increases with temperature


### Example 17.11 Page Number : 611¶

In [2]:


# Variables
# SO2 + (1/2)*O2 - SO3
R = 8.314;			#[J/mol-K]

K_800 = 0.0319;			# Equilibrium constant at 800 K
K_900 = 0.153;			# Equilibrium constant at 900 K
T_1 = 800.;			#[K]
T_2 = 900.;			#[K]

# Calculations
# We have the relation
# math.log(K_2/K_1) = -(delta_H_rkn/R)*(1/T_2 - 1/T_1)
# math.log(K_900/K_800) = -(delta_H_rkn_850/R)*(1/T_2 - 1/T_1)
delta_H_rkn_850 = -R*math.log(K_900/K_800)/(1/T_2 - 1/T_1);			#[J]
delta_H_rkn_850 = delta_H_rkn_850*10**(-3);			#[kJ]

# Results
print " The mean standard enthalpy change of reaction in the region 800 t0 900 is given by delta_H_rkn_850 = %f KJ"%(delta_H_rkn_850);

 The mean standard enthalpy change of reaction in the region 800 t0 900 is given by delta_H_rkn_850 = 93.851672 KJ


### Example 17.13 Page Number : 611¶

In [12]:


# Variables
T_1 = 298.15;			#[K]
T = 2600;			#[K]
R = 1.987;			#[cal/mol-K] - Universal gas constant

# Cp_0 = a + b*T + c*T**(2) + d*T**(3)
delta_H_CO_298 = -26.416;			#[kcal/mol] - Enthalpy of formation of CO at 298.15 K
delta_G_CO_298 = -32.808;			#[kcal/mol] - Gibbs free energy change for formation of CO at 298.15 K
delta_H_CO2_298 = -94.052;			#[kcal/mol] - Enthalpy of formation of C02 at 298.15 K
delta_G_CO2_298 = -94.260;			#[kcal/mol] - Gibbs free energy change for formation of CO2 at 298.15 K

# CO + (1/2)*O2 - CO2

a_CO = 6.726;
a_O2 = 6.0685;
a_CO2 = 5.316;
b_CO = 0.04001*10**(-2);
b_O2 = 0.3631*10**(-2);
b_CO2 = 1.4285*10**(-2);
c_CO = 0.1283*10**(-5);
c_O2 = -0.1709*10**(-5);
c_CO2 = -0.8362*10**(-5);
d_CO = -0.5307*10**(-9);
d_O2 = 0.3133*10**(-9);
d_CO2 = 1.784*10**(-9);

# Calculations and Results
delta_H_rkn_298 = delta_H_CO2_298 - delta_H_CO_298;			#[kcal]
delta_H_rkn_298 = delta_H_rkn_298*10**(3);			#[cal]
delta_G_rkn_298 = delta_G_CO2_298 - delta_G_CO_298;			#[kcal]
delta_G_rkn_298 = delta_G_rkn_298*10**(3);			#[cal]

delta_a = a_CO2 - (a_CO) - (a_O2/2);
delta_b = b_CO2 - (b_CO) - (b_O2/2);
delta_c = c_CO2 - (c_CO) - (c_O2/2);
delta_d = d_CO2 - (d_CO) - (d_O2/2);

def f47(T):
return delta_a+(delta_b*T)+(delta_c*T**(2))+(delta_d*T**(3))

# delta_H_rkn_T = delta_H_rkn_298 +  quad(f47,T_1,T)[0]

# On simplification we get
delta_H_rkn_T = -66773.56 - 4.45*T + 0.605*10**(-2)*T**(2) - 0.29*10**(-5)*T**(3) + 0.54*10**(-9)*T**(4);

#def f48(T):
#	 return K/K_298) = integrate(delta_H_rkn_T/(R*T**(2))

# math.log(K/K_298) =  quad(f48,T_1,T)[0]

# We know that  delta_G_rkn_T = -R*T*math.log(K)
# At 298.15 K
K_298 = math.exp(-delta_G_rkn_298/(R*T_1) );

# Therfore on simplification we get
#math.log(K) = 2.94 + 33605.2/T - 2.24*math.log(T) + 0.304*10(-2)*T - 0.073*10**(-5)*T**(2) + 0.09*10**(-9)*T**(3)
K = math.exp(2.94 + 33605.2/T - 2.24*math.log(T) + 0.304*10**(-2)*T - 0.073*10**(-5)*T**(2) + 0.09*10**(-9)*T**(3));

print "  The value of equilibrium constant at 2600 K is given by K_298 = %f"%(K);

#(a)
P_1 = 1;			#[atm]
Kp_1 = P_1**(-1./2);
Ky_1 = K/Kp_1;

# Let the reaction coordinate at equilibrium for the reaction be X
# At equilibrium, the moles of the components be
# n_CO_1_eq = 1 - X
# n_02_1_eq = 0.5- 0.5X
# n_CO2_1_eq = X
# Total moles = 1.5 - 0.5*X

# Ky = y_CO2/(y_CO**(1/2)*y_CO)
#ky = (X*(1.5-0.5*X)**(1/2))/((1-X)*(0.5-0.5*X)**(1/2))

def f(X):
return  Ky_1-(X*(1.5-0.5*X)**(1./2))/((1-X)*(0.5-0.5*X)**(1./2))
X_1 = fsolve(f,0.9)

y_CO2_1 = X_1/(1.5-0.5*X_1);
y_CO_1 = (1-X_1)/(1.5-0.5*X_1);
y_O2_1 = (0.5-0.5*X_1)/(1.5-0.5*X_1);

print " a).The equilibrium composition at 1 atm) is given by y_CO2 = %f y_CO = %f and y_O2 = %f"%(y_CO2_1,y_CO_1,y_O2_1);

#(b)
P_2 = 10;			#[atm]
Kp_2 = P_2**(-1./2);
Ky_2 = K/Kp_2;

# Ky = y_CO2/(y_CO**(1/2)*y_CO)
#ky = (X*(1.5-0.5*X)**(1/2))/((1-X)*(0.5-0.5*X)**(1/2))

def f1(X):
return  Ky_2-(X*(1.5-0.5*X)**(1./2))/((1-X)*(0.5-0.5*X)**(1./2))
X_2 = fsolve(f1,0.9)

y_CO2_2 = X_2/(1.5-0.5*X_2);
y_CO_2 = (1-X_2)/(1.5-0.5*X_2);
y_O2_2 = (0.5-0.5*X_2)/(1.5-0.5*X_2);

print " b).The equilibrium composition at 10 atm) is given by y_CO2 = %f y_CO = %f and y_O2 = %f"%(y_CO2_2,y_CO_2,y_O2_2);

#(c)
P_3 = 1;			#[atm]
Kp_3 = P_3**(-1./2);
Ky_3 = K/Kp_3;

# Ky = y_CO2/(y_CO**(1/2)*y_CO)
#ky = (X*(1.5-0.5*X)**(1/2))/((1-X)*(0.5-0.5*X)**(1/2))

# At equilibrium, the moles of the components be
# n_CO_3_eq = 1 - X
# n_02_3_eq = 0.5- 0.5X
# n_CO2_3_eq = X
# n_N2_eq = 1;
# Total moles = 2.5 - 0.5*X

def f2(X):
return  Ky_3-(X*(2.5-0.5*X)**(1./2))/((1-X)*(0.5-0.5*X)**(1./2))
X_3 = fsolve(f2,0.9)

y_CO2_3 = X_3/(2.5-0.5*X_3);
y_CO_3 = (1-X_3)/(2.5-0.5*X_3);
y_O2_3 = (0.5-0.5*X_3)/(2.5-0.5*X_3);
y_N2 = 1./(2.5-0.5*X_3);

print " c).The equilibrium composition at 1 atm and 1 mol N2) is given by y_CO2 = %f y_CO = %f y_O2 = %f and y_N2 = %f"%(y_CO2_3,y_CO_3,y_O2_3,y_N2);

  The value of equilibrium constant at 2600 K is given by K_298 = 16.484152
a).The equilibrium composition at 1 atm) is given by y_CO2 = 0.757531 y_CO = 0.161646 and y_O2 = 0.080823
b).The equilibrium composition at 10 atm) is given by y_CO2 = 0.876008 y_CO = 0.082662 and y_O2 = 0.041331
c).The equilibrium composition at 1 atm and 1 mol N2) is given by y_CO2 = 0.373822 y_CO = 0.100943 y_O2 = 0.050471 and y_N2 = 0.474764


### Example 17.14 Page Number : 614¶

In [13]:


# Variables
T = 25 + 298.15;			#[K] - Temperature
R = 8.314;			#[J/mol-K]
delta_G_298 = -1000;			#[J] - Gibbs free energy change at 298 K

# G_E/(R*T) = x_1*x_2
# Calculations and Results
# We know that  delta_G_rkn = - R*T*math.log(K), therefore
K = math.exp(-delta_G_298/(R*T));

#(1)
# Let x_1 is the mole fraction of A and x_2 be the mole fraction of B
# If A and B form an ideal mixture then,
Ky = 1;
# and K = Kx = x_2/x_1
# and at equilibrium x_2/x_1 = K
# (1-x_1)/x_1 = K
x_1 = 1/(1+K);

print " 1).The equilibrium composition for ideal behaviour) is given by x_1 = %f"%(x_1);

#(2)
# The activity coefficients are given by
# math.log(Y1) = [del(n*G_E/(R*T))/del(n_1)]_T,P,n_2 = x_2**(2)
# similarly, math.log(Y2) = x_1**(2)

# The equilibrium constant for the liquid phase reaction is given by
# K = Kx*Ky = (x_2*Y2)/(x_1*Y1) = (x_2*math.exp(x_1**(2)))/(x_1*math.exp(x_2**(2)))
# Solving the above equation we get

def f2(x_1):
return  K -((1-x_1)*math.exp(x_1**(2)))/(x_1*math.exp((1-x_1)**(2)))
x_1_prime = fsolve(f2,0.9)

print " 2).The equilibrium composition for non-ideal behaviour) is given by x_1 = %f"%(x_1_prime);

 1).The equilibrium composition for ideal behaviour) is given by x_1 = 0.408008
2).The equilibrium composition for non-ideal behaviour) is given by x_1 = 0.328409


### Example 17.15 Page Number : 615¶

In [14]:


# Variables
T_1 = 298.15;			#[K] - Smath.tan(math.radiansard reaction temperature
T_2 = 373.;			#[K]  - Reaction temperature
P = 1.;			#[atm]
R = 1.987;			#[cal/mol-K] - Universal gas constant

# CH3COOH (l) + C2H5OH (l) - CH3COOC2H5 (l) + H2O (l)

delta_H_CH3COOH_298 = -116.2*10**(3.);			# [cal/mol]
delta_H_C2H5OH_298 = -66.35*10**(3.);			# [cal/mol]
delta_H_CH3COOC2H5_298 = -110.72*10**(3.);			# [cal/mol]
delta_H_H2O_298 = -68.3174*10**(3.);			# [cal/mol]

delta_G_CH3COOH_298 = -93.56*10**(3.);			# [cal/mol]
delta_G_C2H5OH_298 = -41.76*10**(3.);			# [cal/mol]
delta_G_CH3COOC2H5_298 = -76.11*10**(3.);			# [cal/mol]
delta_G_H2O_298 = -56.6899*10**(3.);			# [cal/mol]

# Calculations and Results
delta_H_rkn_298 = delta_H_CH3COOC2H5_298 + delta_H_H2O_298 - delta_H_CH3COOH_298 - delta_H_C2H5OH_298;			#[cal/mol]
delta_G_rkn_298 = delta_G_CH3COOC2H5_298 + delta_G_H2O_298 - delta_G_CH3COOH_298 - delta_G_C2H5OH_298;			#[cal/mol]

# We know that  delta_G_rkn_T = -R*T*math.log(K)
# At 298.15 K
K_298 = math.exp(-delta_G_rkn_298/(R*T_1) );

# We know that dmath.log(K)/dT = delta_H_rkn/(R*T**(2))
# If delta_H_rkn is assumed constant we get
# math.log(K_2/K_1) = (-delta_H_rkn/R)*(1/T_2 - 1/T_1)
# math.log(K_373/K_298) = (-delta_H_rkn_298/R)*(1/T_2 - 1/T_1)

K_373 = math.exp(math.log(K_298) + (-delta_H_rkn_298/R)*(1./T_2 - 1./T_1));
# Note that the equilibrium constant value rises becauses the reaction is endothermic

print " The value of equilibrium constant at 373 K is  K_373 = %f"%(K_373);

# Let the reaction coordinate at equilibrium for the reaction be X
# At equilibrium, the moles of the components be
# n_CH3COOH = 1 - X
# n_C2H5OH = 1 - X
# n_CH3COOC2H5 = X
# n_H20 = X
# Total moles = 2

# Kx = (x_CH3COOH*x_C2H5OH)/(x_CH3COOC2H5*x_H2O)
# Assuming the liquid mixture to be ideal,that is Ky = 1, therefore K_x = K
K_x = K_373;
# X**(2)/(1-X)**(2) = K_x
X = (K_x)**(1./2)/(1+(K_x)**(1./2));

# The mole fraction of ethyl acetate is given by
x_CH3COOC2H5 = X/2;

print " The mole fraction of ethyl acetate in the equilibrium reaction mixture is given by x_CH3COOC2H5 = %f"%(x_CH3COOC2H5);

 The value of equilibrium constant at 373 K is  K_373 = 0.046697
The mole fraction of ethyl acetate in the equilibrium reaction mixture is given by x_CH3COOC2H5 = 0.088848


### Example 17.16 Page Number : 617¶

In [15]:


# Variables
# CaCO3 (s1) - CaO (s2) + CO2 (g)
T_1 = 898 + 273.15;			#[K]
T_2 = 700 + 273.15;			#[K]
R = 8.314;			#[J/mol-K] - Universal gas constant

P_CO2_T_1 = 1;			#[atm] - Decomposition pressure of  CO2 over CaCO3 at 898 C
P_CO2_T_2 = 0.0333;			#[atm] - Decomposition pressure of  CO2 over CaCO3 at 700 C

# The equilibrium constant of the reaction is given by
# K = (a_CO2*a_CaO)/a_CaCO3

# Since the pressure is small therefore carbon dioxide can be assumed to behave as ideal gas and thus
# a_CO2 = y_CO2*P/1 = P_CO2

# The activity of CaO is (CaO is pure)
# a_CaO = f_CaO/f_0_CaO = exp[V_CaO*(P - P_0)/(R*T)] = 1  (since pressure is low)

# The activity of CaCO3 is (CaCO3 is pure)
# a_CaCO3 = f_CaCO3/f_0_CaCO3 = exp[V_CaCO3*(P - P_0)/(R*T)] = 1  (since pressure is low)

#Since pressures are around 1 atm,therefore Poynting factors are also near 1, and thus activity of CaO and CaCO3 is unity and equilibrium constant is given by
#K = P_CO2 , therefore
# Calculations
# At 898 C
K_T_1 = P_CO2_T_1;
delta_G_T_1 = -R*T_1*math.log(K_T_1);

# At 700 C
K_T_2 = P_CO2_T_2;
delta_G_T_2 = -R*T_2*math.log(K_T_2);

# Results
print " The value of delta_G_rkn at 700 C is %f J"%(delta_G_T_1);
print " The value of delta_G_rkn at 898 C is %f J"%(delta_G_T_2);

 The value of delta_G_rkn at 700 C is -0.000000 J
The value of delta_G_rkn at 898 C is 27526.397496 J


### Example 17.17 Page Number : 618¶

In [3]:


# Variables
T = 700 + 273.15;			#[K]
K = 7.403;			# Equilibrium constant for the reaction at 700 C

# CH4 - C (s) + 2*H2

# The equilibrium constant of the reaction is given by
# K = (a_C*a_H2**(2))/a_CH4

# Since carbon is pure therefore its activity is given by
# a_C = f/f_0 = 1 as pressure is 1 atm.
# Since the pressure is low,therefore the gas phase can be taken to be ideal,therefore
# K = (y_H2**(2)*P**(2))/(y_CH4*P) = y_H2**(2)/y_CH4     (as P = 1 atm)
Ky = K;			# (Kp = 1 atm)

# Initial moles of the species are
n_CH4 = 1;
n_H2 = 0;
n_C = 0;

# Let the reaction coordinate at equilibrium for the reaction be X
# Moles at equilibrium be
# n_CH4_eq = 1 -X
# n_H2_eq = 2*x
# n_C_eq = X

# Moles present in gas phase
# n_CH4_gas = 1 -X
# n_H2_gas = 2*x
# Total moles = 1 + X

# gas phase mole fraction
# y_CH4_gas = (1 -X)/(1+X)
# y_H2_gas = 2*x/(1+X)

# Calculations and Results
# Ky = y_H2_gas**(2)/y_CH4_gaS
X = (K/(4+K))**(1./2);

print " The number of moles of carbon black formed from one mole of methane  is %f mol"%(X);

# Therefore mole fraction of components in the gas phase at equilibrium is
y_CH4 = (1-X)/(1+X);
y_H2 = 2*X/(1+X);

print " The mole fraction of components in the gas phase at equilibrium is given by y_CH4 = %f and y_H2 = %f "%(y_CH4,y_H2);

 The number of moles of carbon black formed from one mole of methane  is 0.805739 mol
The mole fraction of components in the gas phase at equilibrium is given by y_CH4 = 0.107580 and y_H2 = 0.892420


### Example 17.18 Page Number : 619¶

In [17]:

# Variables
T_1 = 298.15;			#[K] - Smath.tan(math.radiansard reaction temperature
T_2 = 1042;			#[K] - Reaction temperature
R = 1.987;			#[cal/mol-K] - Universal gas constant

# CaCO3 (s1) - CaO (s2) + CO2 (g)

delta_H_CaCO3_298 = -289.5;			#[kcal/mol] - Enthalpy of formation of CaCO3 at 298.15 K
delta_H_CaO_298 = -151.7;			#[kcal/mol] - Enthalpy of formation of CaO at 298.15 K
delta_H_CO2_298 = -94.052;			#[kcal/mol] - Enthalpy of formation of CO2 at 298.15 K
delta_G_CaCO3_298 = -270.8;			#[kcal/mol] - Gibbs free energy change for formation of CaCO3 at 298.15 K
delta_G_CaO_298 = -144.3;			#[kcal/mol] - Gibbs free energy change for formation of CaO at 298.15 K
delta_G_CO2_298 = -94.260;			#[kcal/mol] - Gibbs free energy change for formation of CO2 at 298.15 K

# The standard heat capacity data as a function of temperature are given below
# Cp_CO2 = 5.316 + 1.4285*10**(2)*T - 0.8362*10**(-5)*T**(2) + 1.784*10**(-9)*T**(3)
# Cp_CaO = 12.129 + 0.88*10**(-3)*T - 2.08*10**(5)*T**(-2)
# Cp_CaCO3 = 24.98 + 5.240*10**(-3)*T - 6.199*10**(5)*T**(-2)
# Therefore Cp_0 is given by
# Cp_0 = Cp_CO2 + Cp_CaO - Cp_CaCO3
# Cp_0 = -7.535 + 9.925*10**(-3)*T - 0.8362*10**(-5)*T**(2) + 1.784*10**(-9)*T**(3) + 4.119*10**(5)*T**(-2)
# Calculations and Results
# The smath.tan(math.radiansard enthalpy change of the reaction at 298.15 K is given by
delta_H_rkn_298 = delta_H_CO2_298 + delta_H_CaO_298 - delta_H_CaCO3_298;			#[kcal]
delta_H_rkn_298 = delta_H_rkn_298*10**(3);			#[cal]
delta_G_rkn_298 = delta_G_CO2_298 + delta_G_CaO_298 - delta_G_CaCO3_298;			#[kcal]
delta_G_rkn_298 = delta_G_rkn_298*10**(3);			#[cal]

# The smath.tan(math.radiansard enthalpy change of the reaction at temperature T is given by

def f7(T):
return -7.535 + 9.925*10**(-3)*T - 0.8362*10**(-5)*T**(2) + 1.784*10**(-9)*T**(3) + 4.119*10**(5)*T**(-2)

# delta_H_rkn_T = delta_H_rkn_298 +  quad(f7,T_1,T)[0]

# On simplification we get
# delta_H_rkn_T = 47005.3 - 7.535*T + 4.9625*10**(-3)*T**(2) - 0.2787*10**(-5)*T**(3) + 0.446*10**(-9)*T**(4) - 4.119*10**(5)*T**(-1)

#def f8(T):
#	 return K_2/K_1) = integrate(delta_H_rkn_T/(R*T**(2))

# math.log(K_2/K_1) =  quad(f8,T_1,T)[0]

def f9(T):
return (47005.3-7.535*T+4.9625*10**(-3)*T**(2)-0.2787*10**(-5)*T**(3)+0.446*10**(-9)*T**(4) - 4.119*10**(5)*T**(-1))/T**(2)

math.log_K2_K1 =  quad(f9,T_1,T_2)[0];			# math.log(K_2/K_1)[0]

# We know that  delta_G_rkn_T = -R*T*math.log(K)
# At 298.15 K
K_298 = math.exp(-delta_G_rkn_298/(R*T_1) );

# Putting the values in the above math.expression we get
# math.log(K_1042/K_298) = math.log_K2_K1/R
K_1042 = K_298*math.exp(math.log_K2_K1/R);

print " The value of equilibrium constant at 1042 K is K_1042 = %f"%(K_1042);

# Since for the given reaction K = P_CO2,where P is in atm, therefore,
P_CO2 = K_1042;
# and thus decomposition takes place till the partial pressure of CO2 reaches 0.095 atm
# After that the decomposition in the closed vessel stops as equilibrium is achieved.

print " The equilibrium pressure of CO2 is P_CO2 = %f atm "%(P_CO2);

 The value of equilibrium constant at 1042 K is K_1042 = 0.095017
The equilibrium pressure of CO2 is P_CO2 = 0.095017 atm


### Example 17.19 Page Number : 620¶

In [18]:


# Variables
T_1 = 298.15;			#[k] - Smath.tan(math.radiansard reaction temperature
T_2 = 1200;			#[K] - Reaction temperature
R = 1.987;			#[cal/mol-K] - Universal gas consatnt

# C (s) + CO2 (g) - 2*CO2 (g)  			# Reaction 1
# CO2 + H2 - CO + H2O  			# Reacction 2

K_1 = 63;			# Equilibrium constant for the first reaction
K_2 = 1.4;			# Equilibrium constant for the secind reaction

delta_G_H2O_298 = -54640;			#[cal/mol] - Smath.tan(math.radiansard Gibbs free energy of formation of water vapour
delta_H_H2O_298 = -57800;			#[cal/mol] - Smath.tan(math.radiansard enthalpy change of formation of water vapour
delta_G_rkn_298 = delta_G_H2O_298;

# The smath.tan(math.radiansard heat capacity data of the components in cal/mol-K are given below
# Cp_H2 = 6.947 - 0.2*10**(-3)*T + 4.8*10**(-7)*T**(2)
# Cp_O2 = 6.148 + 3.1*10**(-3)*T - 9.2*10**(-7)*T**(2)
# Cp_H2O = 7.256 + 2.3*10**(-3)*T + 2.8*10**(-7)*T**(2)

# Let us consider two more reactions
# C (s) + (1/2)*O2 - CO 			# Reaction 3
# H2 + (1/2)*O2 - H2O  			# Reaction 4

# Considering the above 4 reactions, it can be shown that reaction (3) = reaction (1) + reaction (4) - reaction (2)
# Thus,  delta_G_rkn_3 = delta_G_rkn_1 + delta_G_rkn_4 - delta_G_rkn_2
# or, -R*T*math.log(K_3) = -R*T*math.log(K_1) + -R*T*math.log(K_4) - -R*T*math.log(K_2)
# K_3 = (K_1*K_4/K_2)

# Now we have to determine K_4 at 1200 K.
# The smath.tan(math.radiansard enthalpy change of reaction (4) as a fuction of temperature is given by
# delta_H_rkn_T = -57020 - 2.765*T + 0.475*10**(-3)*T**(2) + 8.67*10**(-8)*T**(3);

#def f2(T):
#	 return K_4_2/K_4_1) = integrate(delta_H_rkn_T/(R*T**(2))

# math.log(K_4_2/K_4_1) =  quad(f2,T_1,T)[0]

# Calculations
def f3(T):
return (-57020-2.765*T+0.475*10**(-3)*T**(2)+8.67*10**(-8)*T**(3))/T**(2)

# We know that  delta_G_rkn_T = -R*T*math.log(K)
# At 298.15 K
K_4_298 = math.exp(-delta_G_rkn_298/(R*T_1) );

# Putting the values in the above math.expression we get
# math.log(K_4_1200/K_4_298) = math.log_K2_K1_R/R
K_4_1200 = K_4_298*math.exp(math.log_K2_K1_4/R);
K_4 = K_4_1200;

# Therefore the equilibrium constant for reaction (3) at 1200 K is given by
K_3 = (K_1*K_4)/K_2;

# Results
print " The equilibrium constant at 1200 K for the given reaction is K = %e"%(K_3);

 The equilibrium constant at 1200 K for the given reaction is K = 3.622432e+09


### Example 17.20 Page Number : 622¶

In [5]:

# Variables
delta_G_H2O_298 = -237.034;			#[kJ/mol] - Smath.tan(math.radiansard Gibbs free energy of formation of H2O (l) at 298 K
F = 96485;			#[C/mol] - Faraday constant

#(1)
# For the reaction
# H2 (g) + (1/2)*O2 (g) - H2O (l)
n = 2;			# Number of electrons transferred in the reaction

# Calculations and Results
# The standard Gibbs free energy change of the reaction is given by
# delta_G_rkn = delta_G_for_H2O(l) - delta_G_for_H2(g) - (1/2/)*delta_G_for_O2(g)
# Since delta_G_for_H2  = 0 and delta_G_for_O2 = 0  (pure components)
delta_G_rkn = delta_G_H2O_298;			#[kJ]
delta_G_rkn = delta_G_rkn*10**(3);			#[J]

# delta_G_rkn = -n*F*E_0
# Thus standard equilibrium cell voltage is given by
E_0 = - delta_G_rkn/(n*F);			#/[V]

print " 1).The standard equilibrium cell voltage is %f V"%(E_0);

#(2)
# For the reaction
# 2*H2 (g) + O2 (g) - 2*H2O (l)
n_prime = 4;			# Number of electrons transferred in the reaction

# The standard Gibbs free energy change of the reaction is given by
# delta_G_rkn = 2*delta_G_for_H2O(l) - 2*delta_G_for_H2(g) - delta_G_for_O2(g)
# Since delta_G_for_H2  = 0 and delta_G_for_O2 = 0  (pure components)
delta_G_rkn_prime = 2*delta_G_H2O_298;			#[kJ]
delta_G_rkn_prime = delta_G_rkn_prime*10**(3);			#[J]

# delta_G_rkn = -n*F*E_0
# Thus standard equilibrium cell voltage is given by
E_0_prime = - delta_G_rkn_prime/(n_prime*F);			#/[V]

print " 2).The standard equilibrium cell voltage is %f V"%(E_0_prime);

# Thus the result obtained is same for both the reactions

 1).The standard equilibrium cell voltage is 1.228346 V
2).The standard equilibrium cell voltage is 1.228346 V


### Example 17.21 Page Number : 624¶

In [20]:

# Variables
P = 2;			# Number of phases
C = 5;			# Number of components

# Calculations
# C + 2*H2 = CH4   			# (reaction 1)
# H2 + (1/2)*O2 = H20  			# (reaction 2)
# C + (1/2)*O2 = CO  			# (reaction 3)
# C + O2 = CO2  			# (reaction 4)

# We do not have C in the equilibrium reaction mixture, therefore we have to eliminate it.
# Substituting for C from reaction (1) into reactions (3) and (4) we get the following set of reactions
# H2 + (1/2)*O2 = H20
# CH4 - 2*H2 + (1/2)*O2 = CO
# CH4 - 2*H2 + O2 = CO2

# We do not have O2 in the equilibrium reaction mixture,therefore we have to eliminateit
# Substituting for O2 from the first reaction of the above set into seecond and third reactions of the above set we get the following set of reactions.
# CH4 + H2O - H2 = CO + 2*H2
# CH4 + 2*H20 - 2*H2 = CO2 + 2*H2

# Therefore one set of independent reactions is
# CH4 + H20 = CO + 3*H2
# CH4 + 2*H2O = CO2 + 4*H2

# Another set of independent reactions can be obtained by substituting for CH4 from the first reaction into second and we get
# CH4 + H2O = CO + 3*H2
# CO + 3*H2 - H2O + 2*H2O = CO2 4*H2

# Or,
# CH4 + H2O = CO + 3*H2
# CO + H2O = CO2 + H2

r = 2;

# and the number of degree of freedom becomes,
F = r - P + C;

# Results
print " The number of independent chemical reactions are %f "%(r);

print " The first set of independent reactions are given below";
print " CH4 + H20 = CO + 3*H2";
print " CH4 + 2*H2O = CO2 + 4*H2";

print " The second set of independent reactions are given below";
print " CH4 + H20 = CO + 3*H2";
print " CO + H2O = CO2 + H2";

 The number of independent chemical reactions are 2.000000
The first set of independent reactions are given below
CH4 + H20 = CO + 3*H2
CH4 + 2*H2O = CO2 + 4*H2
The second set of independent reactions are given below
CH4 + H20 = CO + 3*H2
CO + H2O = CO2 + H2


### Example 17.22 Page Number : 626¶

In [3]:
import math
from scipy.optimize import fsolve
# Variables
T = 400.;			#[K] - Temperature
P = 1.;			#[atm] - Pressure
R = 1.987;			#[cal/mol-K] - Universal gas consatnt

delta_G_n_pentane_400 = 9600.;			#[cal/mol] - standard Gibbs free energy of formation of n-pentane at 400 K
delta_G_iso_pentane_400 = 8200.;			#[cal/mol] - standard Gibbs free energy of formation of iso-pentane at 400 K
delta_G_neo_pentane_400 = 9000.;			#[cal/mol] - standard Gibbs free energy of formation of neo-pentane at 400 K

# The three reactions for the formation of these isomers can be written as follows
# 5*C + 6*H2 = n-pentane
# 5*C + 6*H2 = iso-pentane
# 5*C + 6*H2 = neo-pentane

# Calculations and Results
# Let us take the following set of independent reactions
# iso-pentane = n-pentane   			#     (reaction 1)
delta_G_rkn_1 = delta_G_n_pentane_400 - delta_G_iso_pentane_400;			#[cal]
K_1 = math.exp(-delta_G_rkn_1/(R*T));
# iso-pentane = neo-pentane 			#     (reaction 2)
delta_G_rkn_2 = delta_G_neo_pentane_400 - delta_G_iso_pentane_400;			#[cal]
K_2 = math.exp(-delta_G_rkn_2/(R*T));

# Let the initial number of moles be
# n_iso_pentane = 1
# n_n_pentane = 0
# n_neo_pentane = 0

# Let the reaction coordinate at equilibrium for the two reaction be X_1 and X_2 respectively
# At equilibrium, the moles of the components be
# n_iso_pentane_eq = 1 - X_1 - X_2
# n_n_pentane_eq = X_1
# n_neo_pentane_eq = X_2
# Total moles = 1

# Pressure has no effect on these reactions ( P = 1 atm) and therefore
Ky_1 = K_1;
Ky_2 = K_2;

# From reaction (1), we get
# Ky_1 = X_1/(1-X_1-X_2)

# From reaction (2), we get
# Ky_2 = X_2/(1-X_1-X_2)

# Putting the value of X_1 from first equation into the second we get
# X_1 = (Ky_1*(1-X_2))/(1+Ky_1)
# Now putting it into the second equation we grt
# Ky_2 = X_2/(1-((Ky_1*(1-X_2))/(1+Ky_1))-X_2)
# Now solving for X_2
def f(X_2):
return  Ky_2 - X_2/(1-((Ky_1*(1-X_2))/(1+Ky_1))-X_2)
X_2 = fsolve(f,0.8)

# Therefore X_1 can be calculated as
X_1 = (Ky_1*(1-X_2))/(1+Ky_1);

# Finally the mole fractions of the components at equilibrium
y_n_pentane = X_1;
y_neo_pentane = X_2;
y_iso_pentane = 1 -X_1 - X_2;

print " The equilibrium composition is given by y_n_pentane = %f y_neo_pentane = %f and y_iso_pentane = %f"%(y_n_pentane,y_neo_pentane,y_iso_pentane);

# Let us consider another set of independent reactions

# n-pentane = iso-pentane   			#     (reaction 3)
delta_G_rkn_3 = delta_G_iso_pentane_400 - delta_G_n_pentane_400;			#[cal]
K_3 = math.exp(-delta_G_rkn_3/(R*T));
# n-pentane = neo-pentane 			#     (reaction 4)
delta_G_rkn_4 = delta_G_neo_pentane_400 - delta_G_n_pentane_400;			#[cal]
K_4 = math.exp(-delta_G_rkn_4/(R*T));

# Let the initial number of moles be
# n_n_pentane = 1
# n_iso_pentane = 0
# n_neo_pentane = 0

# Let the reaction coordinate at equilibrium for the two reaction be X_3 and X_4 respectively
# At equilibrium, the moles of the components be
# n_n_pentane_eq = 1 - X_3 - X_4
# n_iso_pentane_eq = X_4
# n_neo_pentane_eq = X_4
# Total moles = 1

# Pressure has no effect on these reactions (P = 1 atm) and therefore
Ky_3 = K_3;
Ky_4 = K_4;

# From reaction (3), we get
# Ky_3 = X_3/(1-X_3-X_4)

# From reaction (4), we get
# Ky_4 = X_4/(1-X_3-X_4)

# Putting the value of X_3 from first equation into the second we get
# X_3 = (Ky_3*(1-X_4))/(1+Ky_3)
# Now putting it into the second equation we grt
# Ky_4 = X_4/(1-((Ky_1*(1-X_4))/(1+Ky_3))-X_4)
# Now solving for X_4
def f1(X_4):
return  Ky_4 - X_4/(1-((Ky_3*(1-X_4))/(1+Ky_3))-X_4)
X_4 = fsolve(f1,0.8)

# Therefore X_3 can be calculated as
X_3 = (Ky_3*(1-X_4))/(1+Ky_3);

# Finally the mole fractions of the components at equilibrium
y_n_pentane1 = 1 - X_3 - X_4;
y_neo_pentane1 = X_4;
y_iso_pentane1 = X_3;

# The final composition does not depend on the set of reactions considered.

print " For another set of independent reactions considered";
print " The equilibrium composition is given by y_n_pentane = %f y_neo_pentane = %f and y_iso_pentane = %f"%(y_n_pentane1,y_neo_pentane1,y_iso_pentane1);

 The equilibrium composition is given by y_n_pentane = 0.111753 y_neo_pentane = 0.237745 and y_iso_pentane = 0.650501
For another set of independent reactions considered
The equilibrium composition is given by y_n_pentane = 0.111753 y_neo_pentane = 0.237745 and y_iso_pentane = 0.650501


### Example 17.23 Page Number : 628¶

In [6]:

T = 600;			#[K] - Temperature
P = 1;			#[atm] - Pressure
R = 1.987;			#[cal/mol-K] - Universal gas consatnt

# CH4 + H2O = CO + 3*H2 			#    (Reaction 1)
# CO + H2O = CO2 + H2  			#     (Reaction 2)

K_1 = 0.574;			# Equilibrium constant of first reaction
K_2 = 2.21;			# Equilibrium constant of second reaction

# Initial number of moles of the components are
# n_CH4 = 1
# n_H2O = 5
# n_CO = 0
# n_H2 = O
# n_CO2 = 0

# Let the reaction coordinate at equilibrium for the two reaction be X_1 and X_2 respectively
# At equilibrium, the moles of the components be
# n_CH4_eq = 1 - X_1
# n_H20_eq = 5 - X_1 - X_2
# n_CO_eq = X_1 - X_2
# n_H2_eq = 3*X_1 + X_2
# n_CO2_eq = X_2
# Total moles = 6 + 2*X_1

# Since the pressure is 1 atm, K = Kp, Ky = K
Ky_1 = K_1;
Ky_2 = K_2;

# From reaction (1), we get
# Ky_1 = ((X_1-X_2)*(3*X_1+X_2)**(3))/((1-X_1)*(5-X_1-X_2)*(6+2*X_1)**(2))

# From reaction (2), we get
# Ky_2 = (X_2*(3*X_1+X_2))/((X_1-X_2)*(5-X_1-X_2))

# Calculations
# Let us assume a value of X_1
X_1 = 0.1;
def f(X_2):
return  Ky_1-((X_1-X_2)*(3*X_1+X_2)**(3))/((1-X_1)*(5-X_1-X_2)*(6+2*X_1)**(2))

def f1(X_1_prime):
return  Ky_2-(X_2_prime*(3*X_1_prime+X_2_prime))/((X_1_prime-X_2_prime)*(5-X_1_prime-X_2_prime))

fault = 10;
while(fault>0.05):
X_2 = fsolve(f,0.8)
X_2_prime = X_2;
X_1_prime = fsolve(f1,0.8)
fault=abs(X_1 - X_1_prime);
X_1 = X_1 + 0.001;

n_CH4 = 1 - X_1;
n_H20 = 5 - X_1 - X_2;
n_CO = X_1 - X_2;
n_H2 = 3*X_1 + X_2;
n_CO2 = X_2;
Total_moles = 6 + 2*X_1;

# Results
print " The equilibrium composition of the resulting mixture is given by";
print " n_CH4 = % f mol n_H2O = %f mol n_CO = %f mol n_H2 = %f mol and n_CO2 = %f mol"%(n_CH4,n_H20,n_CO,n_H2,n_CO2);
print " The total number of moles is %f mol"%(Total_moles);

 The equilibrium composition of the resulting mixture is given by
n_CH4 =  0.089000 mol n_H2O = 3.471420 mol n_CO = 0.293420 mol n_H2 = 3.350580 mol and n_CO2 = 0.617580 mol
The total number of moles is 7.822000 mol

/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)
/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 five Jacobian evaluations.
warnings.warn(msg, RuntimeWarning)


### Example 17.24 Page Number : 631¶

In [23]:

T = 600 + 273.15;			#[K] - Reaction temperature
P = 1;			#[atm] - Reaction pressure

# The Gibbs free energy of formation of various species at 873.15 K are
delta_G_CH4_RT = -2.82;			# delta_G_CH4/(R*T)
delta_G_H2O_RT = -29.73;			# delta_G_CH4/(R*T)
delta_G_CO_RT = -27.51;			# delta_G_CH4/(R*T)
delta_G_H2_RT = -1.46;			# delta_G_CH4/(R*T)
delta_G_CO2_RT = -56.68;			# delta_G_CH4/(R*T)

# Calculations
# Initial number of moles of the components are
# n_CH4 = 1
# n_H2O = 5
# n_CO = 0
# n_H2 = O
# n_CO2 = 0

# The del(F)/del(n_i) = 0 equations for CH4 (1), H2O (2), CO (3), H2 (4) and CO2 (5) becomes
# delta_G_1_T + R*T*math.log((n_1*P)/n) + lambda_C + 4*lambda_H = 0
# delta_G_2_T + R*T*math.log((n_2*P)/n) + 2*lambda_C + lambda_O = 0
# delta_G_3_T + R*T*math.log((n_3*P)/n) + lambda_c + lambda_O = 0
# delta_G_4_T + R*T*math.log((n_4*P)/n) + 2*lambda_H = 0
# delta_G_5_T + R*T*math.log((n_5*P)/n) + lambda_C + 2*lambda_O = 0

# Where n is the total number of moles in the equilibrium reaction mixture. Dividing the above equations by R*T we get
# delta_G_5_T/(R*T) + math.log((n_5*P)/n) + lambda_C/(R*T) + 2*lambda_O/(R*T) = 0

# Substituting the values of delta_G_i_T/(R*T) in the above equations,the full set of equations (including elemental balance) becomes
# -56.68 + math.log(n_5/n) + lambda_C/(R*T) + 2*lambda_O/(R*T) = 0

# The constraint equations are
# n_1 + n_3 + n_5  = 1 			# (moles of C in the reaction mixture = 1)
# 4*n_1 + 2*n_2 + 2*n_4 = 14  			# (moles of H in the reaction mixture = 14)
# n_2 + n_3 + 2*n_5 = 5  			# (moles of O in the raection mixture = 5)

# The total moles are given by
# n = n_1 + n_2 + n_3 + n_4 + n_5

def solution(x):
f = [0,0,0,0,0,0,0,0,0]
f[0]= -2.82 + math.log(x[0]/x[5]) + x[6] + 4*x[7];
f[1] = -29.73 + math.log(x[1]/x[5]) + 2*x[7] + x[8];
f[2] = -27.51 + math.log(x[2]/x[5]) + x[6] + x[8];
f[3] = -1.46 + math.log(x[3]/x[5]) + 2*x[7];
f[4] = -56.68 + math.log(x[4]/x[5]) + x[6] + 2*x[8];
f[5] = x[0] + x[2] +x[4] - 1;
f[6] = 4*x[0] + 2*x[1] + 2*x[3] - 14;
f[7] = x[1] + x[2] +2*x[4] - 5;
f[8] = x[0]+ x[1] + x[2] + x[3] + x[4] - x[5];
return f

x = [0.01,3.5,0.2,3.0,0.5,5.0,2.0,1.0,25.0];
y = fsolve(solution,x)

# Results
print " n_1 = %f mol"%(y[0]);
print " n_2 = %f mol"%(y[1]);
print " n_3 = %f mol"%(y[2]);
print " n_4 = %f mol"%(y[3]);
print " n_5 = %f mol"%(y[4]);
print " n = %f mol"%(y[5]);
print " lambda_C/RT = %f"%(y[6]);
print " lambda_H/RT = %f"%(y[7]);
print " lambda_O/RT = %f"%(y[8]);

print " The Lagrange multiplier values do not have any physical significance";

 n_1 = 0.091556 mol
n_2 = 3.442034 mol
n_3 = 0.258922 mol
n_4 = 3.374855 mol
n_5 = 0.649522 mol
n = 7.816888 mol
lambda_C/RT = 2.667225
lambda_H/RT = 1.149967
lambda_O/RT = 28.250290
The Lagrange multiplier values do not have any physical significance