# Chapter 3 - Chemical Engineering Thermodynamics¶

## Exa 3.1 Page 32¶

In [1]:
from __future__ import division
from sympy import symbols
vold=symbols('vold')
#from scipy.misc import derivative
print "all values in m3/mol"
T=373.15
P=101325
Tc=647.1
Pc=220.55*10**5
w=.345
R=8.314
Tr=T/Tc # reduced temperature
Pr=P/Pc #reduced pressure
b0=.083-.422*Tr**-1.6
b1=.139-.172*Tr**-4.2
B=(b0+w*b1)*R*Tc/Pc
vnew=1
e1=1
vold_=R*T/P+B
print "the soln by virial gas eqn. of volume in m3/mol by Z(T,P) is",vold
while e1>1e-6:
vold_=vnew
#def Fh(vold):
y=P*vold/(R*T)-1-B/vold
#    return y
#ydash=derivative(Fh,vold)#
ydash=y.diff().subs(vold,vold_)
#vnew=vold-Fh(vold)/ydash#
vnew=vold_-y.subs(vold,vold_)/ydash
e1=abs(vold_-vnew)

print "the soln by virial gas eqn. of volume in m3/mol by Z(T,V) is",vold_
#by peng robinson method
k=.37464+1.54226*w-.26992*w**2
s=1+k*(1-Tr**.5)
lpha=s**2
a=.45724*R**2*Tc**2*lpha/Pc
b=.0778*R*Tc/Pc
vnew=b
e2=1
vol=b
fe=0
fd=0
print "the volume of saturated liq. and saturated vapour by peng-robinson method respectively is"
for i in range(0,3):
vol=vnew
y2=vol**3*P+vol**2*(P*b-R*T)-vol*(3*P*b**2+2*b*R*T-a)+P*b**3+b**2*R*T-a*b
ydash2=3*P*vol**2+(P*b-R*T)*2*vol-(3*P*b**2+2*b*R*T-a)
vnew=vol-y2/ydash2
e2=abs(vol-vnew)
if i==1 and e2<1e-6:
fd=vnew
vnew=R*T/P

print vol,fd

#by redlich-kwong method
i=0
a=.42748*R**2*Tc**2.5/Pc
b=.08664*R*Tc/Pc
Vnew=b ; e3=1
print "the vol of saturated liq. and vapour by redlich kwong method respectively are"

from math import sqrt
VV=symbols('VV')
for i in range(0,3):
V=Vnew
y3=VV**3*P-R*T*VV**2-VV*(P*b**2+b*R*T-a/sqrt(T))-a*b/sqrt(T)
deriv=y3.diff().subs(VV,V)
Vnew=V-y3.subs(VV,V)/deriv
e3=abs(Vnew-V)
if i==1 and e3<1e-6:
de=Vnew
Vnew=R*T/P          #for saturated liq.

print V,de
#vander waals method
i=0
a=27*R**2*Tc**2/(64*Pc)
b=R*Tc/(8*Pc)
vnew=b; v=b ; e=1
vv=symbols('vv')
for i in range(0,3):
v=vnew
v3=vv**3*P-vv**2*(P*b+R*T)+a*vv-a*b
deriva=v3.diff().subs(vv,v)
vnew=v-v3.subs(vv,v)/deriva
e4=abs(v-vnew)
if i==1 and e4<10**-6:
sol=vnew
vnew=R*T/P

print "the vol of saturated liq. and vapour by vander waal method respectively are"
print vnew,sol

all values in m3/mol
the soln by virial gas eqn. of volume in m3/mol by Z(T,P) is vold
the soln by virial gas eqn. of volume in m3/mol by Z(T,V) is 0.0302510618121734
the volume of saturated liq. and saturated vapour by peng-robinson method respectively is
0.0306180024673 2.2509412835e-05
the vol of saturated liq. and vapour by redlich kwong method respectively are
2.62463386857767e-5 2.64047733063277e-5
the vol of saturated liq. and vapour by vander waal method respectively are
0.0304707748423797 3.90181676996140e-5


## Exa 3.2 Page 33¶

In [1]:
from sympy import symbols,log,exp
z1=.5 ; P=101.325                                      #given data
a1=14.3916 ; b1=2795.32 ;c1=230.002
a2=16.262 ; b2=3799.887 ;c2=226.346
x1=z1 ; x2=1-z1
T1sat=b1/(a1-log(P))-c1
T2sat=b2/(a2-log(P))-c2
Told_=273 ; e=1
Tnew=x1*T1sat+x2*T2sat
Told = symbols('Told')
while e>1e-6:
Told_=Tnew

P1sat=exp(a1-b1/(Told+c1))
P2sat=exp(a2-b2/(Told+c2))
y1=P-x1*P1sat-x2*P2sat

ydashd=y1.diff().subs(Told,Told_)
Tnew=Told_-y1.subs(Told,Told_)/ydashd
e=abs(Told_-Tnew)

print "the bubble point temp. in Celsius is",Tnew

#calc of dew point
y1=z1 ; y2=1-z1 ; e=1
Tnew=y1*T1sat+y2*T2sat
Told = symbols('Told')
while e>1e-6:
Told_=Tnew
P1sat=exp(a1-b1/(Told+c1))
P2sat=exp(a2-b2/(Told+c2))
y11=1/P-y1/P1sat-y2/P2sat

ydashd=y11.diff().subs(Told,Told_)
Tnew=Told_-y11.subs(Told,Told_)/ydashd
e=abs(Told_-Tnew)

print "the dew point temp. in Celsius is",Tnew

the bubble point temp. in Celsius is 71.8054048538555
the dew point temp. in Celsius is 87.0799975600241


## Exa 3.3 Page 33¶

In [1]:
from sympy import symbols,log,exp
from numpy import zeros
psat=[0,195.75,97.84,50.32]            #given data
z=[0,.3,.3,.4]
bpp=z[(1)]*psat[(1)]+z[(2)]*psat[(2)]+z[(3)]*psat[(3)]            #Bubble point pressure
dpp=1/(z[(1)]/psat[(1)]+z[(2)]/psat[(2)]+z[(3)]/psat[(3)])        #dew point pressure
print "the bubble point pressure and dew point pressure respectively are",bpp,dpp
P=100 ; v=1 ; vnew=1 ; e=1 ; y2=0 ; der=0                #given pressure is between BPP & DPP
k=zeros(5)
for i in range(1,4):
k[(i)]=psat[(i)]/P

while e>1e-6:
v=vnew
for i in range(1,4):
t1=(1-v+v*k[(i)])
y2=y2+z[(i)]*(k[(i)]-1)/t1
der=der-z[(i)]*(k[(i)]-1)**2/t1**2

vnew=v-y2/der
e=abs(vnew-v)
x=zeros(5)
y=zeros(5)
for i in range(1,4):
x[(i)]=z[(i)]/(1-v+v*k[(i)])
y[i]=x[i]*k[i]

print "mol frctn of acetone in liq. phase is",x[(1)]
print "mol frctn of acetone in vapour phase is",y[(1)]
print "mol frctn of acetonitrile in liq. phase is",x[(2)]
print "mol frctn of acetonitrile in vapour. phase is",y[(2)]
print "mol frctn of nitromethane in liq. phase is",x[(3)]
print "mol frctn of nitromethane in vapour phase is",y[(3)]

the bubble point pressure and dew point pressure respectively are 108.205 79.6944627367
mol frctn of acetone in liq. phase is 0.226211512882
mol frctn of acetone in vapour phase is 0.442809036467
mol frctn of acetonitrile in liq. phase is 0.30222391249
mol frctn of acetonitrile in vapour. phase is 0.29569587598
mol frctn of nitromethane in liq. phase is 0.481489923277
mol frctn of nitromethane in vapour phase is 0.242285729393


## Exa 3.4 Page 35¶

In [1]:
from __future__ import division
from math import exp,log
a12=437.98*4.186
a21=1238*4.186
v1=76.92
v2=18.07                                                               #calc of BPP
t=100
x1=.5 ; R=8.314
a1=16.678 ; b1=3640.2 ; c1=219.61
a2=16.2887 ; b2=3816.44 ; c2=227.02
x2=1-x1
p1sat=exp(a1-b1/(c1+t))
p2sat=exp(a2-b2/(c2+t))
h12=v2*exp(-a12/(R*(t+273.15)))/v1
h21=v1*exp(-a21/(R*(t+273.15)))/v2
m=h12/(x1+x2*h12)-h21/(x2+x1*h21)
g1=exp(-log(x1+x2*h12)+x2*m)
g2=exp(-log(x2+x1*h21)-x1*m)
p=x1*g1*p1sat+x2*g2*p2sat
print "boiling point pressure in kPa is",p
#calc of BPT
p=101.325 ; x1=.5 ;  e=1
x2=1-x1
t1sat=b1/(a1-log(p))-c1
t2sat=b2/(a2-log(p))-c2
tnew=x1*t1sat+x2*t2sat
while e>10**-4:
told=tnew
p1sat=exp(a1-b1/(c1+told))
p2sat=exp(a2-b2/(c2+told))
p1sat=p/(g1*x1+g2*x2*(p2sat/p1sat))
tnew=b1/(a1-log(p1sat))-c1
e=abs(tnew-told)

print "boiling point temperature in Celsius is",tnew
#calc of dpp
e1=1 ;  e2=1 ;  e3=1 ; pold=1
t=100 ; y1=.5
y2=1-y1
p1sat=exp(a1-b1/(c1+t))
p2sat=exp(a2-b2/(c2+t))
g1=1 ; g2=1 ; g11=1 ; g22=1
pnew=1/(y1/(g1*p1sat)+y2/(g2*p2sat))
while e1>.0001:
pold=pnew
while e2>.0001 and e3>.0001:
g1=g11 ; g2=g22
x1=y1*pold/(g1*p1sat)
x2=y2*pold/(g2*p2sat)
x1=x1/(x1+x2)
x2=1-x1
h12=v2*exp(-a12/(R*(t+273.15)))/v1
h21=v1*exp(-a21/(R*(t+273.15)))/v2
m=h12/(x1+x2*h12)-h21/(x2+x1*h21)
g11=exp(-log(x1+x2*h12)+x2*m)
g22=exp(-log(x2+x1*h21)-x1*m)
e2=abs(g11-g1) ; e3=abs(g22-g2)

pnew=1/(y1/(g1*p1sat)+y2/(g2*p2sat))
e1=abs(pnew-pold)

print "dew point pressure  in kPa is",pnew
#calc dpt
p=101.325 ; y1=.5 ;  e4=1 ;  e5=1 ; e6=1
y2=1-y1
t1sat=b1/(a1-log(p))-c1
t2sat=b2/(a2-log(p))-c2
tnew=y1*t1sat+y2*t2sat
g11=1 ; g22=1
while e4>.0001:
told=tnew
p1sat=exp(a1-b1/(c1+told))
p2sat=exp(a2-b2/(c2+told))
while e5>.0001 and e6>.0001:
g1=g11
g2=g22
x1=y1*p/(g1*p1sat)
x2=y2*p/(g2*p2sat)
x1=x1/(x1+x2)
x2=1-x1
h12=v2*exp(-a12/(R*(t+273.15)))/v1
h21=v1*exp(-a21/(R*(t+273.15)))/v2
m=h12/(x1+x2*h12)-h21/(x2+x1*h21)
g11=exp(-log(x1+x2*h12)+x2*m)
g22=exp(-log(x2+x1*h21)-x1*m)
e5=abs(g11-g1)
e6=abs(g22-g2)

p1sat=p*(y1/g1+y2*p1sat/(g2*p2sat))
tnew=b1/(a1-log(p1sat))-c1
e4=abs(tnew-told)

print "dew point temperature  in Celsius is",tnew

boiling point pressure in kPa is 204.264012744
boiling point temperature in Celsius is 81.4675225814
dew point pressure  in kPa is 184.104082046
dew point temperature  in Celsius is 84.0803584712


## Exa 3.5 Page 36¶

In [1]:
from sympy import symbols,log,exp
from numpy import zeros
a12=292.66*4.18
a21=1445.26*4.18
v1=74.05*10**-6
v2=18.07*10**-6
R=8.314
t=100
z1=.3
z2=1-z1
a1=14.39155 ; a2=16.262 ; b1=2795.82 ; b2=3799.89 ; c1=230.002 ; c2=226.35
e1=1 ;e2=1 ;e3=1 ;e4=1 ;e5=1 ;e6=1 ;vnew=0
#calc of BPP
x1=z1; x2=z2
p1sat=exp(a1-(b1/(t+c1)))
p2sat=exp(a2-(b2/(t+c2)))
h12=v2*exp(-a12/(R*(t+273.15)))/v1
h21=v1*exp(-a21/(R*(t+273.15)))/v2
m=(h12/(x1+x2*h12))-(h21/(x2+x1*h21))
g1=exp(-log(x1+x2*h12)+x2*m)
g2=exp(-log(x2+x1*h21)-x1*m)
p=x1*g1*p1sat+x2*g2*p2sat
print"the bubble point pressure is",p
bpp=p ;  gb1=g1 ; gb2=g2                          #g1 & g2 are activity co-efficients
#calc of DPP
y1=z1 ; y2=z2
g1=1 ;  g2=1
pnew=1/(y1/(g1*p1sat)+y2/(g2*p2sat))
g1n=g1 ; g2n=g2
while e1>.0001:
pold=pnew
while e2>.0001 and e3>.0001:
g1=g1n ; g2=g2n
x1=y1*pold/(g1*p1sat)
x2=y2*pold/(g2*p2sat)
x1=x1/(x1+x2)
x2=1-x1
g1n=exp(-log(x1+x2*h12)+x2*m)
g2n=exp(-log(x2+x1*h21)-x1*m)
e2=abs(g1n-g1) ; e3=abs(g2n-g2)

pnew=1/(y1/(g1n*p1sat)+y2/(g2n*p2sat))
e1=abs(pnew-pold)

print "the dew point pressure is",pnew
dpp=pnew ; gd1=g1n ; gd2=g2n
p=200
v=(bpp-p)/(bpp-dpp)
g1=((p-dpp)*(gb1-gd1))/(bpp-dpp)+gd1
g2=((p-dpp)*(gb2-gd2))/(bpp-dpp)+gd2

#calc of distribution co-efficients
from sympy import symbols, exp as exp_, log as log_

v=symbols('v')
while e4>.0001 and e5>.0001:
g1n=g1 ; g2n=g2
k1=g1n*p1sat/p
k2=g2n*p2sat/p
while e6>.00001:
vv=vnew
y1=(k1*z1)/(1-v+v*k1)
y2=(k2*z2)/(1-v+v*k2)
x1_=y1/k1
x2_=y2/k2
f=y1-x1_+y2-x2_
derv=f.diff().subs(v,vv)
vnew=vv-f.subs(v,vv)/derv
e6=abs(vnew-vv)

h12=v2*exp(-a12/(R*(t+273.15)))/v1
h21=v1*exp(-a21/(R*(t+273.15)))/v2
m=(h12/(x1+x2*h12))-(h21/(x2+x1*h21))
g1=exp_(-log_(x1+x2*h12)+x2*m)
g2=exp_(-log_(x2+x1*h21)-x1*m)
e4=abs(g1-g1n)
e5=abs(g2-g2n)

print "the no. of moles in vapour phase is",vv
print "x1 and y1 respectively are",x1.subs(v,v1),y1.subs(v,v1)

the bubble point pressure is 310.299302687235
the dew point pressure is 143.432965325576
the no. of moles in vapour phase is 0.540556419677788
x1 and y1 respectively are 0.0283277058050475 1.88246875380055


## Exa 3.6 Page 37¶

In [1]:
from __future__ import division
from scipy import sqrt,log,exp
t=373.15 ; tc=647.1  ; pc=220.55*10**5 ; w=.345 ; R=8.314            #given
f1=1 ; e1=1 ; e2=1 ; vnew=1 ; pnew=1             #assumed values
k=.37464+1.54226*w-.26992*w*2
s=(1+k*(1-(t/tc)**.5))**2
a=.45724*R*R*tc*tc*s/pc
b=.0778*R*tc/pc
#calc of vol. of liq.
while f1>10**-4:
p=pnew ; vnew=b ;
t1=(p*b-R*t)                       #co-efficients of vol. in the eqn.
t2=3*p*b**2+2*b*R*t-a
t3=p*b**3+b**2*R*t-a*b
while e1>1e-6 and vnew>0:
vold=vnew
y1=vold**3*p+vold**2*t1-vold*t2+t3
der=3*vold**2*p+2*vold*t1-t2
vnew=vold-y1/der
e1=abs(vnew-vold)

vliq=vold
#fugacity of liq.
zliq=p*vliq/(R*t)
c=(a/(2*1.414*b*R*t))*(log((vliq+(1+sqrt(2))*b)/(vliq+(1-sqrt(2))*b)))
t5=zliq-p*b/(R*t)
fl2=p*exp(zliq-1-log(t5)-c)
vvnew=R*t/p                 #assumed value close to the ideal value
#calc of vol. of vapour
while e2>1e-6:
vvold=vvnew
x2=vvold**3*p+vvold**2*t1-vvold*t2+t3
der1=3*vvold**2*p+2*vvold*t1-t2
vvnew=vvold-x2/der1
e2=abs(vvnew-vvold)

#fugacity of vapour
vvap=vvold
zvap=p*vvap/(R*t)
d=(a/(2*sqrt(2)*b*R*t))*(log((zvap+(1+sqrt(2))*b*p/(R*t))/(zvap+(1-sqrt(2))*p*b/(R*t))))
t6=zvap-p*b/(R*t)
fv2=p*exp(zvap-1-log(t6)-d)
pnew=p*fl2/fv2                      #updating the value of P
f1=abs((fl2-fv2)/fv2)

print "the vapour pressure of water in Pa is",p
# Note value of fv2 is infinite & it is the reason for warning

the vapour pressure of water in Pa is 154846.968055

/home/asif/.local/lib/python2.7/site-packages/ipykernel/__main__.py:42: RuntimeWarning: overflow encountered in exp
/home/asif/.local/lib/python2.7/site-packages/ipykernel/__main__.py:44: RuntimeWarning: invalid value encountered in double_scalars


## Exa 3.7 Page 39¶

In [1]:
from __future__ import division
from scipy import sqrt,log,exp
et=1 ; er=1 ; sold=0 ; snew=0             #assumed values
R=8.314 ; t=100 ; x1=.958 ; a12=107.38*4.18 ; a21=469.55*4.18
tc1=512.6 ; tc2=647.1 ; pc1=80.97*10**5 ; pc2=220.55*10**5 ; w1=.564 ; w2=.345 ; zc1=.224 ; zc2=.229 ; v1=40.73*10**-6 ; v2=18.07*10**-6                     #given data
x2=1-x1
a1=16.5938 ; a2=16.262 ; b1=3644.3 ; b2=3799.89 ; c1=239.76 ; c2=226.35
p1sat=exp(a1-b1/(c1+t))*1000                         #Saturation Pressure
p2sat=exp(a2-b2/(c2+t))*1000
t=t+273.15                                           #Temp in Kelvin req.
h12=(v2/v1)*exp(-a12/(R*t))
h21=(v1/v2)*exp(-a21/(R*t))
z=h12/(x1+x2*h12)-h21/(x2+x1*h21)
g1=exp(-log(x1+x2*h12)+x2*z)                          #Activity Co-efficients
g2=exp(-log(x2+x1*h21)-x1*z)
tr1=t/tc1                                 #Reduced Temp.
b0=.083-.422*tr1**-1.6
b1=.139-.172*tr1**-4.2
b11=(R*tc1/pc1)*(b0+b1*w1)
tr2=t/tc2
b0=.083-.422*tr2**-1.6
b1=.139-.172*tr2**-4.2
b22=(R*tc2/pc2)*(b0+b1*w2)
w12=(w1+w2)**.5
tc12=(tc1*tc2)**.5
zc12=(zc1+zc2)**.5
vc1=zc1*R*tc1/pc1 ; vc2=zc2*R*tc2/pc2
vc12=((vc1**.33+vc2**.33)/2)**3
pc12=zc12*R*tc12/vc12
tr12=t/tc12
b0=.083-.422*tr12**-1.6
b1=.139-.172*tr12**-4.2
b12=(R*tc12/pc12)*(b0+b1*w12)
d12=2*b12-b11-b22
p=x1*g1*p1sat+x2*p2sat*g2
y1=x1*g1*p1sat/p ; y2=x2*g2*p2sat/p
pnew=p
#calc of Pressure
while et>1e-6:
p=pnew
f1=p1sat*(exp(b11*p1sat/(R*t)))*(exp((v1*(p-p1sat)/(R*t))))
f2=p2sat*(exp(b22*p2sat/(R*t)))*(exp((v2*(p-p2sat)/(R*t))))
while er>1e-6:
sold=snew
fc1=exp((p/(R*t))*(b11+y2**2*d12))
fc2=exp((p/(R*t))*(b22+y1**2*d12))
k1=g1*f1/(fc1*p)
k2=g2*f2/(fc2*p)
snew=x1*k1+x2*k2
y1=x1*k1/snew
y2=x2*k2/snew
er=abs(snew-sold)

pnew=(x1*g1*f1/fc1)+(x2*g2*f2/fc2)
y1=x1*g1*f1/(fc1*pnew)
y2=x2*g2*f2/(fc2*pnew)
et=abs(pnew-p)

print "the amt. of methanol in vapour phase  and system pressure in Pa respectively are",y1,p

the amt. of methanol in vapour phase  and system pressure in Pa respectively are 0.981547559133 344637.196487


## Exa 3.9 Page 41¶

In [1]:
from __future__ import division
from scipy import sqrt,log,exp
#let x1 and x2 be the reaction co-ordinate for 1st and 2nd reactions
x1new=.9 ; x2new=.6 ; r1=1 ; r2=1               #assumed values
Kp=1             #since P=1 atm
K1=.574 ; K2=2.21               #given
Kye1=K1 ; Kye2=K2               #at eqm.
while r1>1e-6 and r2>1e-6:
x1=x1new ; x2=x2new
m_CH4=1-x1
m_H2O=5-x1-x2
m_CO=x1-x2
m_H2=3*x1+x2
m_CO2=x2      #moles of reactants and products at eqm.
total=m_CO2+m_H2+m_CO+m_H2O+m_CH4
Ky1=m_CO*m_H2**3/(m_CH4*m_H2O*total**2)
Ky2=m_CO2*m_H2/(m_CO*m_H2O)
f1=Ky1-.574                             #1st function in x1 and x2
f2=Ky2-2.21                             #2nd function in x1 and x2
d3=((3*x1+x2)**2*(12*x1-8*x2))/((1-x1)*(5-x1-x2)*(6+2*x1)**2)
d4=(3*x1+x2)**3*(x1-x2)*(8*x1**2+6*x1*x2-24*x1+2*x2-16)
d5=((1-x1)**2)*((5-x1-x2)**2)*((6+2*x1)**3)
df1_dx1=d3-(d4/d5)                         #df1/dx1- partial derivative of f1 wrt to x1
d6=3*(x1-x2)*((3*x1+x2)**2)-(3*x1+x2)**3
d7=(1-x1)*(5-x1-x2)*((6+x1*2)**2)
d8=((x1-x2)*(3*x1+x2)**3)/((1-x1)*((5-x1-x2)**2)*(6+2*x1)**2)
df1_dx2=(d6/d7)+d8                         #df1/dx2- partial derivative of f1 wrt to x2
d9=(x1-x2)*(5-x1-x2)
df2_dx1=3*x2/d9-(x2*(3*x1+x2)*(5-2*x1))/(d9**2)        #df1/dx2- partial derivative of f1 wrt to x2
d10=(3*x1+2*x2)/d9
d11=x2*(3*x1+x2)*(2*x2-5)/(d9**2)
df2_dx2=d10-d11                             #df1/dx2- partial derivative of f1 wrt to x2
dm=df1_dx1*df2_dx2-df1_dx2*df2_dx1
delta_x1=(f2*df1_dx2-f1*df2_dx2)/dm
delta_x2=(f1*df2_dx1-f2*df1_dx1)/dm
x1new=x1+delta_x1                        #updating the values of x1 & x2
x2new=x2+delta_x2
r1=abs(x1-x1new)
r2=abs(x2new-x2)

print "the value of X1 and X2 respectively is",x1,x2
m_CH4=1-x1 ; m_H2O=5-x1-x2 ; m_CO=x1-x2 ; m_H2=3*x1+x2 ; m_CO2=x2      #moles of reactants and products at eqm.
total=m_CO2+m_H2+m_CO+m_H2O+m_CH4
print "the moles at eqm of CH4,H2O,CO,H2,CO2 are",m_CH4,m_H2O,m_CO,m_H2,m_CO2
print "total number of moles at eqm. is",total

the value of X1 and X2 respectively is 0.912056278984 0.632839210308
the moles at eqm of CH4,H2O,CO,H2,CO2 are 0.0879437210164 3.45510451071 0.279217068676 3.36900804726 0.632839210308
total number of moles at eqm. is 7.82411255797


## Exa 3.10 Page 45¶

In [1]:
from __future__ import division
from scipy import sqrt,log,exp
from sympy.abc import t

u1=1 ; u2=3.5 ; u3=2 ; u4=3                      #moles given 1-C2H6, 2-O2, 3-CO2, 4-H2O
a1=1.648 ; a2=6.085 ; a3=5.316 ; a4=7.7
b1=4.124e-2 ; b2=.3631e-2 ; b3=1.4285e-2 ; b4=.04594e-2
c1=-1.53e-5 ; c2=-.1709e-5 ; c3=-.8362e-5 ; c4=.2521e-5
d1=1.74e-9 ; d2=.3133e-9 ; d3=1.784e-9 ; d4=-.8587e-9
n1=1 ; n2=4 ; n3=10 ; n4=0 ; t0=298.15 ; t1=25 ; e1=1
t1=t1+273.15
#calc. of sum of co-efficients of heat capacity of the rxn.
sa=n1*a1+n2*a2+n3*a3+n4*a4
sb=n1*b1+n2*b2+n3*b3+n4*b4
sc=n1*c1+n2*c2+n3*c3+n4*c4
sd=n1*d1+n2*d2+n3*d3+n4*d4
da=u4*a4+u3*a3-u2*a2-u1*a1
db=u4*b4+u3*b3-u2*b2-u1*b1
dc=u4*c4+u3*c3-u2*c2-u1*c1
dd=u4*d4+u3*d3-u2*d2-u1*d1
h0=(u4*(-57.798)+u3*(-94.05)-u2*0-u1*(-20.236))*1000                #enthalpy of the rxn.
tnew=1000

while e1>1e-6:
tt=tnew
ft=sa*(t-t1)+(sb/2)*(t**2-t1**2)+(sc/3)*(t**3-t1**3)+(sd/4)*(t**4-t1**4)+h0+da*(t-t0)+(db/2)*(t**2-t0**2)+(dc/3)*(t**3-t0**3)+(dd/4)*(t**4-t0**4)
dr=ft.diff().subs(t,tt)
tnew=tt-ft.subs(t,tt)/dr
e1=abs((tnew-tt)/tnew)

print "the adiabatic flame temp in K is",tnew

the adiabatic flame temp in K is 2090.37132921579