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);
# 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);
# 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);
# 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);
# 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"
# 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);
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);
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";
# 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);
# 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);
# 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);
# 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);
# 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);
# 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);
# 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);
# 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)
math.log_K2_K1_4 = quad(f3,T_1,T_2)[0]
# 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);
# 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
# 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";
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);
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);
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";