# Chapter 13 : Chemical Reaction Equilibria¶

## Example 13.1 page no : 216¶

In [1]:
# Variables
#CH4 + H2O --> CO + 3H2
n_CH4 = 2;			#Moles of CH4
n_H2O = 1;			#Moles of H20
n_CO = 1;			#Moles of CO
n_H2 = 4;			#Moles of H2

v_CH4 = -1;
v_H2O = -1
v_CO = 1;
v_H2 = 3;

# Calculations and Results
v = v_CH4+v_H2O+v_CO+v_H2;
n = n_CH4+n_H2O+n_CO+n_H2;

#y_CH4 = (n_CH4+(v_CH4e)/n+(v*e))
#y_H2O = (n_H2O+(v_H2Oe)/n+(v*e))
#y_CO = (n_CO+(v_CO*e)/n+(v*e))
#y_H2 = (n_H2+(v_H2*e)/n+(v*e))

y_CH4 = '(n_CH4+(v_CH4e)/n+(v*e))';
y_H2O = '(n_H2O+(v_H2Oe)/n+(v*e))';
y_CO = '(n_CO+(v_CO*e)/n+(v*e))';
y_H2 = '(n_H2+(v_H2*e)/n+(v*e))';

#Hence

y_CH4 = '(2-e/8+2e)';
y_H2O = '(1-e/8+2e)';
y_CO = '(1+e/8+2e)';
y_H2 = '(4+3e/8+2e)';

print 'y_CH4  =  ',y_CH4
print 'y_H2O  =  ',y_H2O
print 'y_CO  =  ',y_CO
print 'y_H2  =  ',y_H2

y_CH4  =   (2-e/8+2e)
y_H2O  =   (1-e/8+2e)
y_CO  =   (1+e/8+2e)
y_H2  =   (4+3e/8+2e)


## Example 13.3 page no : 218¶

In [1]:
# Program to Determine the Expression for yi for two reactions

#CH4 + H2O  --> CO + 3H2    (A)

#CH4 + 2H2O --> CO2 + 4H2   (B)

# Variables
Species = ['CH4','H2O','CO','CO2','H2','sum'];
v_A = ['-1','-1','1','0','3','2'];
v_B = ['-1','-2','0','1','4','2'];
m_CH4 = 2;
m_H2O = 3;

# Calculations
mt = m_CH4+m_H2O;
y = ['(m_CH4+vie1+vje2)/(mt+vie1+vje2)','(m_H2O+vie1+vje2)/(mt+vie1+vje2)','(vie1)/(mt+vie1+vje2)','(vje2)/(mt+vie1+vje2)','(vie1+vje2)/(mt+vie1+vje2)',' '];

#Hence
yf = ['(2-e1-e2)/(5+2e1+2e2)','(3-e1-2e2)/(5+2e1+3e2)','e1/(5+2e1+2e2)','e2/(5+2e1+2e2)','(3e1+4e2)/(5+2e1+2e2)',' '];
Ans = [Species,v_A,v_B,y,yf];

# Results
print ' i   v_A  v_B      y(Before substitution)         y_species'
print Ans

 i   v_A  v_B      y(Before substitution)         y_species
[['CH4', 'H2O', 'CO', 'CO2', 'H2', 'sum'], ['-1', '-1', '1', '0', '3', '2'], ['-1', '-2', '0', '1', '4', '2'], ['(m_CH4+vie1+vje2)/(mt+vie1+vje2)', '(m_H2O+vie1+vje2)/(mt+vie1+vje2)', '(vie1)/(mt+vie1+vje2)', '(vje2)/(mt+vie1+vje2)', '(vie1+vje2)/(mt+vie1+vje2)', ' '], ['(2-e1-e2)/(5+2e1+2e2)', '(3-e1-2e2)/(5+2e1+3e2)', 'e1/(5+2e1+2e2)', 'e2/(5+2e1+2e2)', '(3e1+4e2)/(5+2e1+2e2)', ' ']]


## Example 13.4 page no : 219¶

In [3]:
import math

def IDCPH(T0,T,dA,dB,dC,dD):
t = T/T0;
return (dA+((dB/2)*T0*(t+1))+((dC/3)*T0*T0*((t**2)+t+1))+(dD/(t*T0*T0)))*(T-T0)

def IDCPS(T0,T,dA,dB,dC,dD):
t = T/T0;
return ((dA)*math.log(t))+(((dB*T0)+(((dC*T0*T0)+(dD/(t*t*T0*T0)))*(t+1)/2))*(t-1))

# Variables
T0 = 298.16;			#[K]
T1 = 418.15;			#[K]
T2 = 593.15;			#[K]
R = 8.314;

#C2H4(g) + H2O(g) -->  C2H5OH(g)
#Values From Table C.1 At T = 298.15K

A_ethanol = 3.518;
A_ethene = 1.424;
A_water = 3.470;

# Calculations
B_ethanol = 20.001*10**-3;
B_ethene = 14.394*10**-3;
B_water = 1.450*10**-3;

C_ethanol = -6.002*10**-6;
C_ethene = -4.392*10**-6;
C_water = 0;

D_ethanol = 0;
D_ethene = 0;
D_water = 0.121*10**5;

dA = A_ethanol-A_ethene-A_water
dB = B_ethanol-B_ethene-B_water
dC = C_ethanol-C_ethene-C_water
dD = D_ethanol-D_ethene-D_water

# Values from Table C.4 at T = 298.15K
H_ethanol = -235100;			#[J/mol]
H_ethene = 52510;			#[J/mol]
H_water = -241572;			#[J/mol]

G_ethanol = -168490;			#[J/mol]
G_ethene = 68460;			#[J/mol]
G_water = -228572;			#[J/mol]

dHo = H_ethanol-H_ethene-H_water
dGo = G_ethanol-G_ethene-G_water

I1 = round(IDCPH(T0,T1,dA,dB,dC,dD),3)
I2 = round(IDCPS(T0,T1,dA,dB,dC,dD),5)

#Using Eqn 13.18
#dG_418/RT = ((dGo - dHo)/RTo)+(dHo/RT)+((1/T)*I1)-I2  c1 = dG_418/RT

c1 = round(((dGo-dHo)/(R*T0))+(dHo/(R*T1))+((1/T1)*I1)-I2,4)
I3 = round(IDCPH(T0,T2,dA,dB,dC,dD),3)
I4 = round(IDCPS(T0,T2,dA,dB,dC,dD),5)

#Using Eqn 13.18
#dG_593/RT = ((dGo - dHo)/RTo)+(dHo/RT)+((1/T)*I1)-I2  c2 = dG_593/RT

c2 = round(((dGo-dHo)/(R*T0))+(dHo/(R*T2))+((1/T2)*I3)-I4,4)

K_413 = round(math.exp(-c1),4);
K_593 = math.exp(-c2);

# Results
print 'Equilibrium Consmath.tant at T  =  413.15K is ',K_413*10,'X 10**-1'
print 'Equilibrium Constant at T  =  593.15K is ',round(K_593*1000,3),'X 10**-3'

Equilibrium Consmath.tant at T  =  413.15K is  1.404 X 10**-1
Equilibrium Constant at T  =  593.15K is  2.802 X 10**-3


## Example 13.4 page no : 220¶

In [8]:
from numpy import array,round,ones,exp,log
import math

#Example 13.4  Alternate
# Alternate Program to 13.4

# Variables
T = array([298.15,418.15,593.15]);
t = round(T/T[0],4);
R = 8.314;

#C2H4(g) + H2O(g) -->  C2H5OH(g)
#Values From Table C.1 At T = 298.15K

A_ethanol = 3.518;
A_ethene = 1.424;
A_water = 3.470;

# Calculations
B_ethanol = 20.001*10**-3;
B_ethene = 14.394*10**-3;
B_water = 1.450*10**-3;

C_ethanol = -6.002*10**-6;
C_ethene = -4.392*10**-6;
C_water = 0;

D_ethanol = 0;
D_ethene = 0;
D_water = 0.121*10**5;

dA = A_ethanol-A_ethene-A_water
dB = B_ethanol-B_ethene-B_water
dC = C_ethanol-C_ethene-C_water
dD = D_ethanol-D_ethene-D_water

# Values from Table C.4 at T = 298.15K
H_ethanol = -235100;			#[J/mol]
H_ethene = 52510;			#[J/mol]
H_water = -241572;			#[J/mol]

G_ethanol = -168490;			#[J/mol]
G_ethene = 68460;			#[J/mol]
G_water = -228572;			#[J/mol]

dHo = H_ethanol-H_ethene-H_water
dGo = G_ethanol-G_ethene-G_water

#Umath.sing Eqn(13.21)
Ko = round(exp(-dGo/(R*T[0])),3)
K0 = Ko*ones((1,3));
#Umath.sing Eqn(13.22)
K1 = exp((dHo/(R*T0))*(1-(T[0]/T)));
#Umath.sing Eqn(13.24)
K2 = round(exp((dA*(log(t)-((t-1)/t)))+(0.5*dB*T[0]*((t-1)**2)/t)+((1/6.)*dC*T[0]*T[0]*((t-1)**2)*(t+2)/t)+(0.5*dD*((t-1)**2)/((T[0]**2)*(t)**2))),4);

K = K0*K1*K2;

Ans = [T,t,K0,K1,K2,K];

# Results
print '    T/K       t          K0          K1         K2         K',Ans

    T/K       t          K0          K1         K2         K [array([ 298.15,  418.15,  593.15]), array([ 1.    ,  1.4025,  1.9894]), array([[ 29.366,  29.366,  29.366]]), array([  1.00000000e+00,   4.84531951e-03,   9.74032575e-05]), array([ 1.    ,  0.9862,  0.9794]), array([[  2.93660000e+01,   1.40324083e-01,   2.80142097e-03]])]


## Example 13.5 page no : 222¶

In [12]:
import math
from numpy import array,round,ones,exp,log,poly1d,roots

# Variables
#CO(g) + H2O(g) --> CO2(g) + H2(g)
v_CO = -1;
v_H2O = -1;
v_CO2 = 1;
v_H2 = 1;
v = v_CO+v_H2O+v_CO2+v_H2;

#Calculate e(Fraction of Stream) in each case
#(a)
n_H2O_a = 1;			#mol
n_CO_a = 1;			#mol
T_a = 1100;			#[K]
P_a = 1;			#[bar]

x = 10**4/T_a;
#at this x the value of ln K = 0 form Graph
y = 0;
nt = n_H2O_a+n_CO_a;
K = exp(y);

#y_H2O = (n_H2O+(v_H2O*e))/nt
#y_CO = (n_CO+(v_CO*e))/nt
#y_H2 = (v_H2*e)/nt
#y_CO2 = (v_CO2*e)/nt

y_H2O = '(1-e)/2'
y_CO = '(1-e)/2'
y_H2 = 'e/2'
y_CO2 = 'e/2'

#K = (y_H2*y_CO2)/(y_CO*y_H2O)

K = 'e**2/((1-e)**2) = 1';

#Solving
e_a = 0.5;

#(b) same as in (a) P_b = 10bar

#Pressure has no effect on fraction of stream
e_b = e_a;

#(c) same as in (a) N2 = 2mols included

#Since N2 just act as diluent
#It only changes the total moles fraction remains the same
e_c = e_a;

#(d)
n_H2O_d = 2;			#mol
n_CO_d = 1;			#mol
T_d = 1100;			#[K]
P_d = 1;			#[bar]

x = 10**4/T_d;
#at this x the value of ln K = 0 form Graph
y = 0;
nt = n_H2O_d+n_CO_d;
K = math.exp(y);

#y_H2O = (n_H2O+(v_H2O*e))/nt
#y_CO = (n_CO+(v_CO*e))/nt
#y_H2 = (v_H2*e)/nt
#y_CO2 = (v_CO2*e)/nt

y_H2O = '(1-e)/3'
y_CO = '(2-e)/3'
y_H2 = 'e/3'
y_CO2 = 'e/3'

#K = (y_H2*y_CO2)/(y_CO*y_H2O)

K = 'e**2/((1-e)(2-e)) = 1';

#Solving
e = round(2./3,3);
e_d = round(e/2.,3);

#(e)
#Here the y_CO and y_H2O are interchanged

#No change in Fraction of stream
e_e = e_d;

#(f)
n_H2O_f = 1;			#mol
n_CO_f = 1;			#mol
n_CO2_f = 1;			#[mol]
T_f = 1100;			#[K]
P_f = 1;			#[bar]

x = 10**4/T_f;
#at this x the value of ln K = 0 form Graph
y = 0;
nt = n_H2O_f+n_CO_f+n_CO2_f;
K = math.exp(y);

#y_H2O = (n_H2O+(v_H2O*e))/nt
#y_CO = (n_CO+(v_CO*e))/nt
#y_CO2 = (n_CO2+(v_CO2*e))/nt
#y_H2 = (v_H2*e)/nt

y_H2O = '(1-e)/3'
y_CO = '(1-e)/3'
y_H2 = '(1+e)/3'
y_CO2 = 'e/2'

#K = (y_H2*y_CO2)/(y_CO*y_H2O)

K = '(e*(e+1))/((1-e)**2) = 1';

#Solving
e_f = round(1./3,3);

#(g)
n_H2O_g = 1;			#mol
n_CO_g = 1;			#mol
T_g = 1650;			#[K]
P_g = 1;			#[bar]

x = 10**4/T_g;
#at this x the value of ln K = 0 form Graph
y = -1.15;
nt = n_H2O_g+n_CO_g;
K = math.exp(y);

#y_H2O = (n_H2O+(v_H2O*e))/nt
#y_CO = (n_CO+(v_CO*e))/nt
#y_H2 = (v_H2*e)/nt
#y_CO2 = (v_CO2*e)/nt

y_H2O = '(1-e)/2'
y_CO = '(1-e)/2'
y_H2 = 'e/2'
y_CO2 = 'e/2'

#K = (y_H2*y_CO2)/(y_CO*y_H2O)

Exp = 'e**2/((1-e)**2) = 0.316';

#Solving
p = poly1d([K, -2*K, K-1],'e','c');

root = roots(p);
e_g = round(root[0],2);

#Other Root is negative and the Fraction of steam cannot be negative

# Results
print ('(a)1mol H20 and 1 mol CO T = 1100K  P = 1bar')
print ('(b)1mol H20 and 1 mol CO T = 1100K  P = 10bar')
print ('(c)1mol H20 and 1 mol CO 2mol N2 T = 1100K  P = 1bar')
print ('(d)2mol H20 and 1 mol CO T = 1100K  P = 1bar')
print ('(e)1mol H20 and 2 mol CO T = 1100K  P = 1bar')
print ('(f)1mol H20 and 1 mol CO 1mol CO2 T = 1100K  P = 1bar')
print ('(g)1mol H20 and 1 mol CO T = 1650K  P = 1bar')
e = [e_a, e_b, e_c, e_d, e_e, e_f, e_g];
print 'Fraction Of Steam Reacted in each case',e

(a)1mol H20 and 1 mol CO T = 1100K  P = 1bar
(b)1mol H20 and 1 mol CO T = 1100K  P = 10bar
(c)1mol H20 and 1 mol CO 2mol N2 T = 1100K  P = 1bar
(d)2mol H20 and 1 mol CO T = 1100K  P = 1bar
(e)1mol H20 and 2 mol CO T = 1100K  P = 1bar
(f)1mol H20 and 1 mol CO 1mol CO2 T = 1100K  P = 1bar
(g)1mol H20 and 1 mol CO T = 1650K  P = 1bar
Fraction Of Steam Reacted in each case [0.5, 0.5, 0.5, 0.33400000000000002, 0.33400000000000002, 0.33300000000000002, -0.68000000000000005]


## Example 13.6 page no : 224¶

In [13]:
import math
from numpy import array,round,ones,exp,log

def IDCPH(T0,T,dA,dB,dC,dD):
t = T/T0;
return (dA+((dB/2)*T0*(t+1))+((dC/3)*T0*T0*((t**2)+t+1))+(dD/(t*T0*T0)))*(T-T0)

def IDCPS(T0,T,dA,dB,dC,dD):
t = T/T0;
return  ((dA)*math.log(t))+(((dB*T0)+(((dC*T0*T0)+(dD/(t*t*T0*T0)))*(t+1)/2))*(t-1))

def PHIB(Tr,Pr,omega):
B0 = 0.083-(0.422/(Tr**1.6));
B1 = 0.139-(0.172/(Tr**4.2));
return exp((Pr/Tr)*(B0+(omega*B1)));

# Variables
T0 = 298.16;			#[K]
T1 = 523.15;			#[K]
P = 35.;    			#[bar]
R = 8.314;

#C2H4(g) + H2O(g) -->  C2H5OH(g)
#Values From Table C.1 At T = 298.15K

A_ethanol = 3.518;
A_ethene = 1.424;
A_water = 3.470;

B_ethanol = 20.001*10**-3;
B_ethene = 14.394*10**-3;
B_water = 1.450*10**-3;

C_ethanol = -6.002*10**-6;
C_ethene = -4.392*10**-6;
C_water = 0;

D_ethanol = 0;
D_ethene = 0;
D_water = 0.121*10**5;

# Calculations and Results
dA = A_ethanol-A_ethene-A_water
dB = B_ethanol-B_ethene-B_water
dC = C_ethanol-C_ethene-C_water
dD = D_ethanol-D_ethene-D_water

# Values from Table C.4 at T = 298.15K
H_ethanol = -235100;			#[J/mol]
H_ethene = 52510;			#[J/mol]
H_water = -241572;			#[J/mol]

G_ethanol = -168490;			#[J/mol]
G_ethene = 68460;			#[J/mol]
G_water = -228572;			#[J/mol]

dHo = H_ethanol-H_ethene-H_water
dGo = G_ethanol-G_ethene-G_water

I1 = round(IDCPH(T0,T1,dA,dB,dC,dD),3)
I2 = round(IDCPS(T0,T1,dA,dB,dC,dD),5)

#Umath.sing Eqn 13.18
#dG_418/RT = ((dGo - dHo)/RTo)+(dHo/RT)+((1/T)*I1)-I2  c1 = dG_418/RT

c1 = round(((dGo-dHo)/(R*T0))+(dHo/(R*T1))+((1/T1)*I1)-I2,4)
K_523 = round(exp(-c1),4);
print 'Equilibrium Consmath.tant at T  =  523.15K is ',round(K_523*1000,3),'X 10**-3'

#Values Frm App B
Tc = array([282.3,647.1,513.9]);
Pc = array([50.4,220.55,61.48]);
omega = array([0.087,0.345,0.645]);

Tr = round(T1/Tc,3);
Pr = round(P/Pc,3);
si = round(PHIB(Tr,Pr,omega),3);

#Umath.sing eqn
#(y_ETOH*si_ETOH)/(y_C2H4*si_C2H4*y_H20*si_H2O) = (P/Po)K
#y_ETOH/(y_C2H4*y_H20) = c = ((si_C2H4*si_H2O)/si_ETOH)(P/Po)K
c = round(((si[0]*si[1])/si[2])*(P*K_523),3)

#y_C2H4  =  (1-e)/(6-e)
#y_ETOH  =  (5-e)/(6-e)
#y_H2O  =  (e)/(6-e)

#Solving we get a Eqn
#poly([1.342 -6 1],'e','c')
root = round(roots(poly1d([1.342, -6, 1],'e','c')),3)

r = root[0]*100;
#Since e > 1 not possible so e = 0.233

print 'The Maximum Conversion of ethylene to ethanol by Vapor-Phase Hydration is ',r,'%'

Equilibrium Consmath.tant at T  =  523.15K is  9.8 X 10**-3
The Maximum Conversion of ethylene to ethanol by Vapor-Phase Hydration is  -600.0 %


## Example 13.7 page no : 226¶

In [14]:
from numpy import roots, poly1d
import math

# Variables
T = 1393.15;			#K
P = 1.;			#[bar]
x = 10.**4/T;
#C2H2 --> 2C + H2  (I)
#2C + 2H2 --> C2H4 (II)

# Calculations
#Values Of ln K (at  1000/T )for Reactions I and II from Graph
K_I = math.exp(12.9);
K_II = math.exp(-12.9);
K = K_I*K_II;

#Application in Eqn (13.5)

#y_C2H4/(y_C2H2*y_H2) = c = (P/Po)K
c = P*K;
#y_H2 = y_C2H2 = (1-e)/(2-e)
#y_C2H4 = e/(2-e)

#The Eqn comes out to be
#poly([1 -4 2],'e','c')

root = round(roots(poly1d([1, -4, 2],'e','c')),3)
e = root[0];
#Since e > 1 not possible so e = 0.293
y_C2H2 = round((1-e)/(2-e),3);
y_H2 = y_C2H2;
y_C2H4 = round(e/(2-e),3);

print 'Equilibrium Composition of H2 C2H2 and C2H4 Respectively ',y_C2H4,y_C2H2,y_H2

Equilibrium Composition of H2 C2H2 and C2H4 Respectively  -0.667 0.833 0.833


## Example 13.8 page no : 227¶

In [16]:
from numpy import roots,poly1d
import math

# Variables
T0 = 298.15;
T = 373.15;			#[K]
R = 8.314;
#CH3COOH(l)+C2H5OH(l) --> CH3COOC2H5(l) + H2O(l)

#From Table C.4
dHo_EtAc = -480000.;			#[J]
dHo_H2O = -285830.;			#[J]
dHo_EtOH = -277690.;			#[J]
dHo_AcH = -484500.;			#[J]

dGo_EtAc = -332200.;			#[J]
dGo_H2O = -237130.;			#[J]
dGo_EtOH = -174780.;			#[J]
dGo_AcH = -389900.;			#[J]

# Calculations
dHo_298 = dHo_EtAc+dHo_H2O-dHo_EtOH-dHo_AcH;
dGo_298 = dGo_EtAc+dGo_H2O-dGo_EtOH-dGo_AcH;

K_298 = round(math.exp(-dGo_298/(R*T0)),4);

#Umath.sing Eqn(13.15)
#ln(K_373/K_298) = c = -(dHo_298/R)*((1/373.15)-(1/298.15))
c = round(-(dHo_298/R)*((1/373.15)-(1/298.15)),4);
K_373 = round(K_298*math.exp(c),4);

#x_AcH = x_EtOH = (1-e)/2  and x_EtAc = x_H2O = e/2
#K = (x_EtAc*x_H2O)/(x_AcH*x_EtOH)

#Hence The Eqn is
q = poly1d([K_373, -2*K_373, K_373-1],'e','c')
root = round(roots(q),4)
e = root[0];
#Since Other Root is > 1 hence e = 0.6879
x_EtAc = round(e/2,3);

# Results
print 'Composition of Ethyl Acetate in the Reacting Mixture',x_EtAc

Composition of Ethyl Acetate in the Reacting Mixture -4.859


## Example 13.9 page no : 228¶

In [19]:
import math
from numpy import isreal

def MCPH(T0,T,A,B,C,D):
t = T/T0;
return (A+((B/2)*T0*(t+1))+((C/3)*T0*T0*((t**2)+t+1))+(D/(t*T0*T0)))

def IDCPH(T0,T,dA,dB,dC,dD):
t = T/T0;
return (dA+((dB/2)*T0*(t+1))+((dC/3)*T0*T0*((t**2)+t+1))+(dD/(t*T0*T0)))*(T-T0)

def IDCPS(T0,T,dA,dB,dC,dD):
t = T/T0;
return ((dA)*math.log(t))+(((dB*T0)+(((dC*T0*T0)+(dD/(t*t*T0*T0)))*(t+1)/2))*(t-1))

# Variables
P = 1.;			#[bar]
T0 = 298.15;			#[K]
R = 8.314;

#SO2 + 0.5O2 --> SO3
dHo_298 = -98890.;			#[J/mol]
dGo_298 = -70866.;			#[J/mol]

# Calculations
n_O2_i = 0.5*1.2;			#  Moles O2 Entering
n_N2_i = n_O2_i*(79./21);			#Moles N2 Entering

#n_SO2 = 1-e
#n_O2 = 0.6-(0.5*e)
#n_SO3 = e
#n_N2 = 2.257

#By Energy Balance
#(dHo_298*e)+dHo_P  =  dH  =  0                                   (A)

#dHo_P = Cp*(T-298.15)  Cp = E nCp                                (B)

#Cp_SO2 = R*MCPH(T0,T,5.699,0.801E-3,0,-1.015E+5)
#Cp_O2 = R*MCPH(T0,T,3.639,0.506E-3,0,-0.227E+5)
#Cp_SO3 = R*MCPH(T0,T,8.06,1.056E-3,0,-2.028E+5)
#Cp_N2 = R*MCPH(T0,T,3.28,0.593E-3,0,0.04E+5)

#T = (-(dHo_298*e)/Cp)+T0                                       (C)

#K = (e/(1-e))*((3.857-(0.5*e))/(0.6-(0.5*e)))**0.5              (D)

#ln K  =  ((dHo_298-dGo_298)/(R*T0))-(dHo_298/(R*T))+I1-(I2/T)  (E)

#I1 = IDCPS(T0,T,0.5415,0.002E-3,0,-0.8995E+5)
#I2 = IDCPH(T0,T,0.5415,0.002E-3,0,-0.8995E+5)

#Iteration
A1 = 300.;			#Initial
i = -1.;

while(i == -1):
I1 = IDCPS(T0,A1,0.5415,0.002E-3,0,-0.8995E+5);
I2 = IDCPH(T0,A1,0.5415,0.002E-3,0,-0.8995E+5);
#Applying in Eqn (E)
K  =  math.exp(((dHo_298-dGo_298)/(R*T0))-(dHo_298/(R*A1))+I1-(I2/A1));
#Applying in Eqn (D)
if(isreal(K)):
x = 0;
# p = poly([-0.6*(K**2) 1.7*(K**2) 3.857-(1.6*(K**2)) 0.5*((K**2)-1)],'e','c')
# (0.5*((K**2)-1)*(x**3))+((3.857-(1.6*(K**2)))*(x**2))+(1.7*(K**2)*x)+(-0.6*(K**2))
F_x = (0.5*((K**2)-1)*(x**3))+((3.857-(1.6*(K**2)))*(x**2))+(1.7*(K**2)*x)+(-0.6*(K**2));
F_a = F_x;

x = 1;
F_x = (0.5*((K**2)-1)*(x**3))+((3.857-(1.6*(K**2)))*(x**2))+(1.7*(K**2)*x)+(-0.6*(K**2));
#F_x = (x**3)-(4*x)+1;
F_b = F_x;
root = -100;
A = 0;
B = 1;
i = 1;
while(i == 1):
a = A;
F_a = (0.5*((K**2)-1)*(a**3))+((3.857-(1.6*(K**2)))*(a**2))+(1.7*(K**2)*a)+(-0.6*(K**2));
#F_a = (a**3)-(4*a)+1;
b = B;
F_b = (0.5*((K**2)-1)*(b**3))+((3.857-(1.6*(K**2)))*(b**2))+(1.7*(K**2)*b)+(-0.6*(K**2));
#F_b = (b**3)-(4*b)+1;
x1 = ((a*F_b)-(b*F_a))/(F_b-F_a);
F_x1 = (0.5*((K**2)-1)*(x1**3))+((3.857-(1.6*(K**2)))*(x1**2))+(1.7*(K**2)*x1)+(-0.6*(K**2));
#F_x1 = (x1**3)-(4*x1)+1;

if((F_a*F_x1)<0):
flag = 1;
A = a;
B = x1;
elif ((F_x1*F_b)<0):
flag = 2;
A = x1;
B = b;
x1_a = round(x1,4);
b_a = round(b,4);
a_a = round(a,4);
if(x1_a == b_a):
root = round(x1,5);
i = 0;
break;
elif(x1_a == a_a):
root = round(x1,5);
i = 0;
break;
e = root;
Cp_SO2 = R*MCPH(T0,A1,5.699,0.801E-3,0,-1.015E+5);
Cp_O2 = R*MCPH(T0,A1,3.639,0.506E-3,0,-0.227E+5);
Cp_SO3 = R*MCPH(T0,A1,8.06,1.056E-3,0,-2.028E+5);
Cp_N2 = R*MCPH(T0,A1,3.28,0.593E-3,0,0.04E+5);

n_SO2 = 1-e;
n_O2 = 0.6-(0.5*e);
n_SO3 = e;
n_N2 = 2.257;
if(n_SO2<0 or n_O2<0 or n_SO3<0):
e = 0;

Cp = (n_SO2*Cp_SO2)+(n_O2*Cp_O2)+(n_SO3*Cp_SO3)+(n_N2*Cp_N2);
#Applying in Eqn (C)
B = (-(dHo_298*e)/Cp)+T0;
m = (A1+B)/2;
dT = round(abs(m-A1),2);
if(dT<0.1):
i = 0;
T = round(A1,1);
e = round(e,2);
break;
A1 = m;
i = -1;
else:
i = -1;
A1 = A1+1;

print 'Fraction',e

n_SO2 = 1-e
n_O2 = 0.6-(0.5*e)
n_SO3 = e
n_N2 = 2.257

nt = n_SO2+n_O2+n_SO3+n_N2;

y_SO2 = round(n_SO2/nt,4);
y_O2 = round(n_O2/nt,4);
y_SO3 = round(n_SO3/nt,4);
y_N2 = round(n_N2/nt,4);

# Results
print 'Final Temperature',T
print 'Composition of SO2',y_SO2
print 'Composition of O2',y_O2
print 'Composition of SO3',y_SO3
print 'Composition of N2',y_N2

Fraction 0.77
Final Temperature 855.8
Composition of SO2 0.0662
Composition of O2 0.0619
Composition of SO3 0.2218
Composition of N2 0.6501


## Example 13.12 page no : 229¶

In [21]:
import math

def IDCPH(T0,T,dA,dB,dC,dD):
t = T/T0;
return (dA+((dB/2)*T0*(t+1))+((dC/3)*T0*T0*((t**2)+t+1))+(dD/(t*T0*T0)))*(T-T0)

def IDCPS(T0,T,dA,dB,dC,dD):
t = T/T0;
return ((dA)*math.log(t))+(((dB*T0)+(((dC*T0*T0)+(dD/(t*t*T0*T0)))*(t+1)/2))*(t-1))

# Variables
T0 = 298.16;			#[K]
T1 = 750;			#[K]
R = 8.314;
P = 1.2;			#[bar]

#C4H10 --> C2H4 + C2H6  (I)
#C4H10 --> C3H6 + CH4   (II)

#Values From Table C.1 At T = 298.15K

A_butane = 1.935;
A_ethene = 1.424;
A_ethane = 1.131;
A_propene = 1.637;
A_methane = 1.702;

B_butane = 36.915*10**-3;
B_ethene = 14.394*10**-3;
B_ethane = 19.225*10**-3;
B_propene = 22.706*10**-3;
B_methane = 9.081*10**-3;

C_butane = -11.402*10**-6;
C_ethene = -4.392*10**-6;
C_ethane = -5.561*10**-6;
C_propene = -6.915*10**-6;
C_methane = -2.164*10**-6;

D_butane = 0;
D_ethene = 0;
D_ethane = 0;
D_propene = 0;
D_methane = 0;

dA_I = A_ethene+A_ethane-A_butane;
dA_II = A_methane+A_propene-A_butane;

dB_I = B_ethene+B_ethane-B_butane;
dB_II = B_methane+B_propene-B_butane;

dC_I = C_ethene+C_ethane-C_butane;
dC_II = C_methane+C_propene-C_butane;

dD_I = D_ethene+D_ethane-D_butane;
dD_II = D_methane+D_propene-D_butane;

# Values from Table C.4 at T = 298.15K
H_butane = -125790.;			#[J/mol]
H_ethene = 52510.;			#[J/mol]
H_ethane = -83820.;			#[J/mol]
H_propene = 19710.;			#[J/mol]
H_methane = -74520.;			#[J/mol]

G_butane = -16570;			#[J/mol]
G_ethene = 68460.;			#[J/mol]
G_ethane = -31855.;			#[J/mol]
G_propene = 62205.;			#[J/mol]
G_methane = -50460.;			#[J/mol]

# Calculations
dHo_I = H_ethene+H_ethane-H_butane
dHo_II = H_methane+H_propene-H_butane

dGo_I = G_ethene+G_ethane-G_butane
dGo_II = G_methane+G_propene-G_butane

I1_I = round(IDCPH(T0,T1,dA_I,dB_I,dC_I,dD_I),3)
I1_II = round(IDCPH(T0,T1,dA_II,dB_II,dC_II,dD_II),3)
I2_I = round(IDCPS(T0,T1,dA_I,dB_I,dC_I,dD_I),5)
I2_II = round(IDCPS(T0,T1,dA_II,dB_II,dC_II,dD_II),5)

#Using Eqn 13.18
#dG_418/RT = ((dGo - dHo)/RTo)+(dHo/RT)+((1/T)*I1)-I2  c1 = dG_418/RT

c1_I = round(((dGo_I-dHo_I)/(R*T0))+(dHo_I/(R*T1))+((1/T1)*I1_I)-I2_I,4)
c1_II = round(((dGo_II-dHo_II)/(R*T0))+(dHo_II/(R*T1))+((1/T1)*I1_II)-I2_II,4)

K_I = round(math.exp(-c1_I),4)
K_II = round(math.exp(-c1_II),4)

k = (K_II/K_I)**0.5;
e_I = round(((K_I/P)/(1+(K_I*(1/P)*(1+k)*(1+k))))**0.5,4);

e_II = round(k*e_I,4);

n_C4H10 = 1-e_I-e_II;
n_C2H4 = e_I;
n_C2H6 = e_I;
n_C3H6 = e_II;
n_CH4 = e_II;
nt = n_C4H10+n_C2H4+n_C2H6+n_C3H6+n_CH4;

y_C4H10 = round(n_C4H10/nt,4);
y_C2H4 = round(n_C2H4/nt,4);
y_C2H6 = round(n_C2H6/nt,4);
y_C3H6 = round(n_C3H6/nt,4);
y_CH4 = round(n_CH4/nt,4);

y = [y_C4H10, y_C2H4, y_C2H6, y_C3H6, y_CH4];

# Results
print '  Y_C4H10   Y_C2H4    Y_C2H6    Y_C3H6   Y_CH4'
print y

  Y_C4H10   Y_C2H4    Y_C2H6    Y_C3H6   Y_CH4
[0.0014, 0.049399999999999999, 0.049399999999999999, 0.45000000000000001, 0.45000000000000001]


## Example 13.13 page no : 231¶

In [37]:
%matplotlib inline
import math
from matplotlib.pyplot import plot,suptitle,xlabel,ylabel,legend
from numpy import array,linspace,exp,linalg,zeros

# Variables
n_air = 2.381			#[mol]
n_O2 = 0.21*n_air;
n_N2 = 0.79*n_air;
R = 8.314;

P = 20.;    			#[bar]
T = array([1000, 1100, 1200, 1300, 1400, 1500]);
dG_H2O = array([-192420, -187000, -181380, -175720, -170020, -164310]);
dG_CO = array([-200240, -209110, -217830, -226530, -235130, -243740]);
dG_CO2 = array([-395790, -395960, -396020, -396080, -396130, -396160]);

KI = 'y_H2O/((y_O2)**0.5*y_H2)(P/Po)**-0.5'
KII = 'y_CO/((y_O2)**0.5)(P/Po)**0.5'
KIII = 'y_CO2/y_O2'

n = '3.38+((e2-e1)/2)'
y_H2 = '-e1/n'
y_CO = 'e2/n'
y_O2 = '((0.5(1-e1-e2))-e3)/n'
y_H2O = '(1+e1)/n'
y_CO2 = 'e3/n'
y_N2 = '1.88/n'

KI = '(1+e1)(2n)**0.5*(P/Po)**-0.5'
KII = '(e3*(P/Po)**0.5)/(1-e1-e2-2e3)**0.5*(n/2)**0.5'
KIII = '2e3/(1-e1-e2-2e3)'

# Calculations
K_I = round(exp(-dG_H2O/(R*T)),1)
K_II = round(exp(-dG_CO/(R*T)),1)
K_III = round(exp(-dG_CO2/(R*T)),1)

#Now math.since the values of KI KII KIII valyes are so High the mole fraction of O2 must be very small
#Hence We eleminate O2,Hence 2 Eqns are,

#C + CO2 --> 2CO      (a)
#H2O + C --> H2 + CO  (b)

Ka = '(y_CO**2/y_CO2)*(P/Po)'
Kb = '((y_H2*y_CO)/y_H2O)*(P/Po)'

n = '3.38+(e_a+e_b)'
y_H2 = 'e_b/n'
y_CO = '(2e_a+e_b)/n'
y_H2O = '(1-e_b)/n'
y_CO2 = '(0.5-e_a)/n'
y_N2 = '1.88/n'

Ka = '(2e_a+e_b)**2/((0.5-e_a)*n)*(P/Po)'
Kb = 'e_b(2e_a+e_b)/((1-e_b)*n)*(P/Po)'

dG_new_a = (2*dG_CO)-dG_CO2;
dG_new_b = dG_CO-dG_H2O;

Ka = round(exp((-dG_new_a/(R*T))),3);
Kb = round(exp((-dG_new_b/(R*T))),3);

#Calculation of e_a and e_b

a = 0.1;			#Initial Value

b = 0.7;			#Initial Value

C1 = Ka/20;
C2 = Kb/20;
e_a = zeros(6)
e_b = zeros(6)
for i in range(6):
c = -1;
while(c == -1):
fa = round((((a**2)*(4+C1[i]))+(b**2)+((4+C1[i])*(a*b))+(2.88*C1[i]*a)-(0.5*C1[i]*b)-(1.69*C1[i])),4);
dfax = round(((2*a*(4+C1[i]))+((4+C1[i])*b)+(2.88*C1[i])),4);
dfay = round((2*b)+((4+C1[i])*a)-(0.5*C1[i]),4);

fb = round(((b**2*(1+C2[i]))+((2+C2[i])*a*b)-(C2[i]*a)+(2.38*C2[i]*b)-(3.38*C2[i])),4);
dfbx = round((((2+C2[i])*b)-C2[i]),4);
dfby = round(((2*b*(1+C2[i]))+((2+C2[i])*a)+(2.38*C2[i])),4);

A = [[dfax ,dfay],[dfbx ,dfby]];
B = [-fa,-fb];
Ans = round(linalg.solve(A,B),4);
da = Ans[0];
db = Ans[1];

if(da == 0 and db == 0):
c = 0;
e_a[i] = a;
e_b[i] = b;
break;

a = a+da;
b = b+db;

n = 3.38+(e_a+e_b);
y_H2 = round(e_b/n,3);
y_CO = round(((2*e_a)+e_b)/n,3);
y_H2O = round((1-e_b)/n,3);
y_CO2 = round((0.5-e_a)/n,3);
y_N2 = round(1.88/n,3);

Ans = [T,Ka,Kb,e_a,e_b];
Ans1 = [T,y_H2,y_CO,y_H2O,y_CO2,y_N2];

# Results
plot(T,y_H2,'r-')
plot(T,y_CO,'b-')
plot(T,y_H2O,'g-')
plot(T,y_CO2,'m-')
plot(T,y_N2,'y-')

suptitle('Equllibrium Compositions')
xlabel('T/K')
ylabel('yi')

print '    T/K       Ka         Kb         e_a       e_b'
print Ans
print '    T/K      y_H2     y_CO    y_H2O    y_CO2    y_N2'
print Ans1

    T/K       Ka         Kb         e_a       e_b
[array([1000, 1100, 1200, 1300, 1400, 1500]), array([    1.758,    11.405,    53.155,   194.79 ,   583.343,  1514.118]), array([   2.561,   11.219,   38.609,  110.064,  268.764,  583.577]), array([-0.0506,  0.121 ,  0.3168,  0.4302,  0.4738,  0.4895]), array([ 0.5336,  0.7124,  0.8551,  0.9357,  0.9713,  0.9863])]
T/K      y_H2     y_CO    y_H2O    y_CO2    y_N2
[array([1000, 1100, 1200, 1300, 1400, 1500]), array([ 0.138,  0.169,  0.188,  0.197,  0.201,  0.203]), array([ 0.112,  0.227,  0.327,  0.378,  0.398,  0.405]), array([ 0.121,  0.068,  0.032,  0.014,  0.006,  0.003]), array([ 0.143,  0.09 ,  0.04 ,  0.015,  0.005,  0.002]), array([ 0.487,  0.446,  0.413,  0.396,  0.39 ,  0.387])]