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
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
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)]
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
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)
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
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
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
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