Chapter 9: Distillation¶
Ex9.1: Page 349¶
In [70]:
# Illustration 9.1
# Page: 349
print'Illustration 9.1 - Page: 349\n\n'
# solution
import numpy
#****Data****#
# a:n-heptane b:n-octane
Pt = 760; # [mm Hg]
#*****#
Tempa = 98.4;# [boiling point of A,OC]
Tempb = 125.6;# [boiling point of B,OC]
x = numpy.zeros(6);
y_star = numpy.zeros(6);
alpha = numpy.zeros(6);
# Data = [Temp Pa (mm Hg) Pb(mm Hg)]
Data = [(98.4, 760.0, 333.0),(105.0 ,940.0, 417.0),(110.0, 1050.0, 484.0),(115.0, 1200.0, 561.0),(120.0, 1350.0, 650.0),(125.6 ,1540.0, 760.0)];
for i in range(0,6):
x[i] = (Pt-Data[i][2])/(Data[i][1]-Data[i][2]);# [mole fraction of heptane in liquid]
y_star[i] = (Data[i][1]/Pt)*x[i];
alpha[i] = Data[i][1]/Data[i][2];
print"\t\t T(OC)\t\t\t\t Pa(mm Hg)\t\t\t\t\t\t\t Pb(mm Hg)\t\t\t\t\t\t\t\t x\t\t\t\t\t\t\t\t\t\t\t y*\t\t\t\t\t\t\t\t\t alpha\n"
for i in range(0,6):
print "\t \t ",Data[i][0],"\t \t \t \t",Data[i][1],"\t \t \t \t",Data[i][2],"\t \t \t \t ",round(x[i],3),"\t \t \t \t ",round(y_star[i],3),"\t\t\t\t\t\t\t\t",round(alpha[i],2)
Ex9.2: Page 354¶
In [71]:
# Illustration 9.2
# Page: 354
print'Illustration 9.2 - Page: 354\n\n'
# solution
#****Data****#
# a:water b:ethylaniline
Pt = 760.0; # [mm Hg]
ma1 = 50.0;# [g]
mb1 = 50.0;# [g]
#*******#
# Data = [Temp Pa(mm Hg) Pb(mm Hg)]
Data = [(38.5, 51.1 ,1.0),(64.4 ,199.7, 5.0),(80.6 ,363.9 ,10.0),(96.0, 657.6, 20.0),(99.15 ,737.2 ,22.8),(113.2, 1225.0, 40.0)];
Ma = 18.02;# [kg/kmol]
Mb = 121.1;# [kg/kmol]
for i in range(0,6):
p = Data[i][1]+Data[i][2];
if p==Pt:
pa = Data[4][1];# [mm Hg]
pb = Data[i][2];# [mm Hg]
T = Data[i][0];# [OC]
ya_star = pa/Pt;
yb_star = pb/Pt;
ya1 = ma1/Ma;# [g mol water]
yb1 = mb1/Mb;# [g mol ethylalinine]
Y = ya1*(yb_star/ya_star);# [g mol ethylalinine]
print"The original mixture contained",round(ya1,2),"g mol water and ",round(yb1,2)," g mol ethylanaline\n"
print"The mixture will continue to boil at ",T," degree C, where the equilibrium vapour of the indicated composition,until all the water evaporated together with ",round(Y,3),"g mol ethylaniline\n"
print"The temparature will then rise to 204 degree C, and the equilibrium vapour will be of pure ethylanaline"
Ex9.3: Page 362¶
In [72]:
# Illustration 9.3
# Page: 362
print'Illustration 9.3 - Page: 362\n\n'
import numpy
# solution
#****Data****#
# a:n-C3H8 b:n-C4H10 c:n-C5H12 d:n-C6H14
# Bubble Point Calculation
xa = 0.05;
xb = 0.30;
xc = 0.40;
xd = 0.25;
P = 350;# [kN/square m]
#******#
# Assume:
Temp = 60;# [OC]
x = [0.05, 0.30, 0.40, 0.25];
m = [4.70, 1.70 ,0.62 ,0.25];# [At 60 OC]
# Reference: C5H12
mref = m[3];
Sum = 0;
alpha = numpy.zeros(4)
alpha_x = numpy.zeros(4);
for i in range(0,4):
alpha[i] = m[i]/m[3];
alpha_x[i] = alpha[i]*x[i];
Sum = Sum+alpha_x[i];
# From Eqn. 9.23:
SumF = Sum;
Sum = 0;
mref = 1/SumF;
# Corresponding Temparature from the nomograph:
Temp = 56.8;# [OC]
m = [4.60, 1.60, 0.588, 0.235];# [At 56.8 OC]
for i in range(0,4):
alpha[i] = m[i]/m[2];
alpha_x[i] = alpha[i]*x[i];
Sum = Sum+alpha_x[i];
SumF = Sum;
mref = 1/SumF;
# Corresponding Temparature from the nomograph:
Temp = 56.7;# [OC]
Bt = 56.8;# [OC]
yi = numpy.zeros(4);
for i in range(0,4):
yi[i] = alpha_x[i]/Sum;
print"The Bubble Point is ",Bt," degree C\n"
print"Bubble point vapour composition \n"
print"\t yi\n";
print"\n n-C3\t ",round(yi[0],3)
print"\n n-C4\t ",round(yi[1],3)
print"\n n-C5\t ",round(yi[2],3)
print"\n n-C6\t ",round(yi[3],3)
print"\n \n \n"
# Dew Point Calculation
# Asume:
ya = 0.05;
yb = 0.30;
yc = 0.40;
yd = 0.25;
Temp = 80;# [OC]
y = [0.05, 0.30 ,0.40, 0.25];
m = [6.30 ,2.50 ,0.96 ,0.43];# [At 60 OC]
# Reference: C5H12
mref = m[3];
Sum = 0;
alpha = numpy.zeros(4)
alpha_y = numpy.zeros(4);
for i in range(0,4):
alpha[i] = m[i]/m[3];
alpha_y[i] = y[i]/alpha[i];
Sum = Sum+alpha_y[i];
# From Eqn. 9.29:
SumF = Sum;
Sum = 0;
mref = SumF;
# Corresponding Temparature from the nomograph:
Temp = 83.7;# [OC]
m = [6.60, 2.70, 1.08, 0.47];# [At 56.8 OC]
for i in range(0,4):
alpha[i] = m[i]/m[2];
alpha_y[i] = y[i]/alpha[i];
Sum = Sum+alpha_y[i];
SumF = Sum;
mref = 1.0/SumF;
# Corresponding Temparature from the nomograph:
Temp = 84.0;# [OC]
Dt = 84.0;# [OC]
xi = numpy.zeros(4);
for i in range(0,4):
xi[i] = alpha_y[i]/Sum;
print"The Dew Point is ",Dt," degree C\n"
print"Dew point liquid composition \n"
print"\t xi\n"
print"\n n-C3\t ",round(xi[0],3)
print"\n n-C4\t ",round(xi[1],3)
print"\n n-C5\t ",round(xi[2],3)
print"\n n-C6\t ",round(xi[3],3)
Ex9.4: Page 365¶
In [6]:
# Illustration 9.4
# Page: 365
print'Illustration 9.4 - Page: 365\n\n'
import matplotlib.pyplot as plt
%matplotlib inline
import numpy
from scipy.optimize import fsolve
# solution
#****Data****#
# Basis:
F = 100.0;# [mol feed]
zF = 0.5;
D = 60.0;# [mol]
W = 40.0;# [mol]
#*******#
# From Illustration 9.1, Equilibrium data:
Data = numpy.array([[1.0 ,1.0],[0.655, 0.810],[0.487 ,0.674],[0.312, 0.492],[0.1571 ,0.279],[0, 0]]);
Feed = numpy.array([[0 ,0],[1.0 ,1.0]]);
# The operating line is drawn with a slope -(W/D) to cut the equilibrium line.
def f44(x):
return -((W/D)*(x-zF))+zF
x = numpy.arange(0.2,0.7,0.1);
plt.plot(Data[:,0],Data[:,1],label="Equilibrium Line")
plt.plot(Feed[:,0],Feed[:,1],label="Feed Line")
plt.plot(x,f44(x),label="Operating Line");
plt.grid('on')
ax = pylab.gca()
ax.set_xlabel("Mole fraction of heptane in liquid")
ax.set_ylabel("Mole fraction of heptane in vapour")
plt.legend(loc='lower right');
plt.show()
# The point at which the operating line cuts the equilibrium line has the following composition* temparature:
yd = 0.575;# [mole fraction heptane in vapour phase]
xW = 0.387;# [mole fraction heptane in liquid phase]
Temp = 113;# [OC]
print"mole fraction of heptane in vapour phase ",yd
print"mole fraction of heptane in liquid phase ",xW
print"Temperature is ",Temp," degree C\n"
Ex9.5: Page 366¶
In [74]:
# Illustration 9.5
# Page: 366
print'Illustration 9.5 - Page: 366\n\n'
# solution
import numpy
import pylab
import numpy.linalg as lin
#****Data****#
Pt = 760.0;# [mm Hg]
zFa = 0.5;# [mol fraction benzene]
zFb = 0.25;# [mol fraction toulene]
zFc = 0.25;# [mol fraction o-xylene]
#********#
# Basis:
F = 100.0;# [mol feed]
# For Summtion of Yd_star to be unity, W/D = 2.08
# The Eqn.are
# (1): W+D = F
# (2): W-2.08D = 0
a =numpy.array([[1.0 ,1.0],[1.0 ,-2.08]]);
b = numpy.array([[F*1.0],[0]]);
soln = lin.solve(a,b)
W = soln[0];
D = soln[1];
Sub = ['A','B','C'];
p =numpy.array([1370 ,550, 200]);# [mm Hg]
m = numpy.zeros(3);
zF = [zFa ,zFb, zFc];# [Given]
yd_star = numpy.array([0,0,0]);
xW = numpy.zeros(3);
for i in range(0,3):
m[i] = p[i]/Pt;
yd_star[i]=(zF[i])*((W/D)+1)#/(1+(W/(D*m[i])));
xW[i] = yd_star[i]/m[i];
print"\t \t \t \t \t \t \t \t At W/D = 2.08\n\n\n"
print"Substance \t \t p(mm Hg)\t \t m\t \t \t \t \t \t \t \t \t \t zF\t \t \t \t \t \t \t yd*\t\t\t\t\t\txW\n"
for i in range(0,3):
print "\n",Sub[i]," \t \t \t \t ",p[i],"\t \t \t \t ",m[i],"\t \t \t",m[i],"\t \t \t",zF[i]," \t \t \t",yd_star[i],"\t",xW[i]
Ex9.6: Page 370¶
In [75]:
# Illustration 9.6
# Page: 370
print'Illustration 9.6 - Page: 370\n\n'
# solution
from scipy.optimize import fsolve
import math
#****Data****#
# Basis:
F = 100.0;# [mol]
xF = 0.5;
D = 0.6*100;# [mol]
#******#
W = F-D;# [mol]
# From Illustration 9.1:
alpha = 2.16;# [average value of alpha]
# From Eqn.9.46;
def f45(xW):
return math.log(F*xF/(W*xW))-(alpha*math.log(F*(1-xF)/(W*(1-xW))))
xW = fsolve(f45,0.5);# [mole fraction heptane]
def f46(yD):
return F*xF-((D*yD)+(W*xW))
yD = fsolve(f46,100);# [mole fraction heptane]
print"Mole Fraction of heptane in the distillate is ",round(yD,3),"\n"
print"Mole Fraction of heptane in the residue is ",round(xW,3)," \n"
Ex9.7: Page 371¶
In [76]:
# Illustration 9.7
# Page: 371
from scipy.optimize import fsolve
print'Illustration 9.7 - Page: 371\n\n'
# solution
#****Data****#
# a:benzene b:toulene c:o-xylene
# Assume:
Bt = 100.0;#[OC]
pa = 1370.0;# [mm Hg]
pb = 550.0;# [mm Hg]
pc = 200.0;# [mm Hg]
xFa = 0.5;# [mole fraction]
xFb = 0.25;# [mole fraction]
xFc = 0.25;# [mole fraction]
# Basis:
F = 100.0;# [mol]
D = 32.5;# [mol]
#*******#
ref = pb;
alpha_a = pa/ref;
alpha_b = pb/ref;
alpha_c = pc/ref;
W = F-D;# [mol]
xbW = 0.3;# [mol]
xaW = 0.4;# [mol]
xcW = 0.3;# [mol]
err = 1.0;
while(err>(10**(-1))):
# From Eqn. 9.47:
def f47(xaW):
return math.log(F*xFa/(W*xaW))-(alpha_a*math.log(F*xFb/(W*xbW)))
xaW = fsolve(f47,xbW);
def f48(xcW):
return math.log(F*xFc/(W*xcW))-(alpha_c*math.log(F*xFb/(W*xbW)))
xcW = fsolve(f48,xbW);
xbW_n = 1-(xaW+xcW);
err = abs(xbW-xbW_n);
xbw = xbW_n;
# Material balance:
# for A:
def f49(yaD):
return F*xFa-((D*yaD)+(W*xaW))
yaD = fsolve(f49,100);# [mole fraction benzene]
# For B:
def f50(ybD):
return F*xFb-((D*ybD)+(W*xbW))
ybD = fsolve(f50,100);# [mole fraction toulene]
# For C:
def f51(ycD):
return F*xFc-((D*ycD)+(W*xcW))
ycD = fsolve(f51,100);# [mole fraction o-xylene]
print"The residual compositions are:\n"
print"Benzene:\n",round(xaW,3)
print"Toulene:\n",round(xbW,3)
print"o-xylene:\n",round(xcW,3)
print"\n The composited distillate compositions are:\n"
print"Benzene:\n",round(yaD,3)
print"Toulene:\n",round(ybD,3)
print"o-xylene:\n",round(ycD,3)
#the answers are slightly different in textbook due to approximation while here answers are precise
Ex9.8: Page 388¶
In [77]:
# Illustration 9.8
# Page: 388
print'Illustration 9.8 - Page: 388\n\n'
import numpy.linalg as lin
# solution
#****Data*****#
# a:methanol b:water
Xa = 0.5;# [Wt fraction]
Temp1 = 26.7;# [OC]
Temp2 = 37.8;# [OC]
F1 = 5000.0;# [kg/hr]
#******#
#(a)
Ma = 32.04;# [kg/kmol]
Mb = 18.02;# [kg/kmol]
Xa = 0.5;# [Wt fraction]
Xb = 1-Xa;# [Wt fraction]
Temp1 = 26.7;# [OC]
Temp2 = 37.8;# [OC]
F1 = 5000.0;# [kg/hr];
# Basis: 1hr
F = (F1*Xa/Ma)+(F1*Xb/Mb);# [kmol/hr]
# For feed:
zF = (F1*Xa/Ma)/F;# [mole fracton methanol]
MavF = F1/F;# [kg/kmol]
# For distillate:
xD = (95/Ma)/((95/Ma)+(5/Mb));# [mole fraction methanol]
MavD = 100.0/((95/Ma)+(5/Mb));# [kg/kmol]
# For residue:
xW = (1/Ma)/((1/Ma)+(99/Mb));# [mole fraction methanol]
MavR = 100/((1/Ma)+(99/Mb));# [kg/kmol]
# (1): D+W = F [Eqn.9.75]
# (2): D*xD+W*xW = F*zF [Eqn. 9.76]
# Solvving simultaneously:
a = numpy.array([[1.0 ,1.0],[xD ,xW]]);
b = numpy.array([F,F*zF]);
soln = lin.solve(a,b);
D = soln[0];# [kmol/h]
W = soln[1];# [kmol/h]
print"Quantity of Distillate is", round(D*MavD)," kg/hr\n"
print"Quantity of Residue is ",round(W*MavR)," kg/hr\n"
print"\n"
# (b)
# For the vapour-liquid equilibria:
Tempo = 19.69;# [Base Temp. according to "International Critical Tables"]
BtR = 99.0;# [Bubble point of the residue, OC]
hR = 4179.0;# [J/kg K]
hF = 3852.0;# [J/kg K]
def f52(tF):
return (F1*hF*(tF-Temp1))-((W*MavR)*hR*(BtR-Temp2))
tF = fsolve(f52,Temp1);# [OC]
BtF = 76.0;# [Bubble point of feed, OC]
# For the feed:
delta_Hs = -902.5;# [kJ/kmol]
Hf = ((hF/1000.0)*MavF*(tF-Tempo))+delta_Hs;# [kJ/kmol]
# From Fig 9.27:
HD = 6000.0;# [kJ/kmol]
HLo = 3640.0;# [kJ/kmol]
HW = 6000.0;# [kJ/kmol]
print"The enthalpy of feed is ",round(Hf)," kJ/kmol\n"
print"The enthalpy of the residue is ",round(HW)," kJ/kmol\n"
print"\n"
# (c)
# From Fig.9.27:
# The miium reflux ratio is established by the tie line (x = 0.37 y = 0.71), which extended pass through F,the feed.
# At Dm:
Qm = 62570.0;# [kJ/kmol]
Hg1 = 38610.0;# [kJ/kmol]
# From Eqn. 9.65:
Rm = (Qm-Hg1)/(Hg1-HLo);
print"The minimum reflux ratio is ",round(Rm,4),"\n"
print"\n"
# (d)
# From Fig. 9.28:
Np = 4.9;
# But it include the reboiler.
Nm = Np-1;
print"The minimum number of theoretical trays required is ",round(Nm)," \n"
print"\n"
# (e)
R = 1.5*Rm;
# Eqn. 9.65:
def f53(Q_prime):
return R-((Q_prime-Hg1)/(Hg1-HLo))
Q_prime = fsolve(f53,2);# [kJ/kmol]
def f54(Qc):
return Q_prime-(HD+(Qc/D))
Qc = fsolve(f54,2);# [kJ/hr]
Qc = Qc/3600.0;# [kW]
print"The Condensor heat load is ",round(Qc)," kW\n"
# From Eqn. 9.77:
def f55(Q_dprime):
return (F*Hf)-((D*Q_prime)+(W*Q_dprime))
Q_dprime = fsolve(f55,2);
def f56(Qb):
return Q_dprime-(HW-(Qb/W))
Qb = fsolve(f56,2);# [kJ/hr]
Qb = Qb/3600.0;# [kW]
print"The Reboiler heat load is ",round(Qb)," kW\n"
print"\n"
# (f)
# From Fig: 9.28
Np = 9.0;
# But it is including the reboiler
print"No. of theoretical trays in tower is",Np-1,"\n",
G1 = D*(R+1);# [kmol/hr]
Lo = D*R;# [kmol/hr]
# From Fig. 9.28:
# At the feed tray:
x4 = 0.415;
y5 = 0.676;
x5 = 0.318;
y6 = 0.554;
# From Eqn. 9.64:
def f57(L4):
return (L4/D)-((xD-y5)/(y5-x4))
L4 = fsolve(f57,2);# [kmol/hr]
# From Eqn. 9.62:
def f58(G5):
return (L4/G5)-((xD-y5)/(xD-x4))
G5 = fsolve(f58,2);# [kmol/hr]
# From Eqn. 9.74:
def f59(L5_bar):
return (L5_bar/W)-((y6-xW)/(y6-x5))
L5_bar = fsolve(f59,2);# [kmol/hr]
# From Eqn. 9.72:
def f60(G6_bar):
return (L5_bar/G6_bar)-((y6-xW)/(x5-xW))
G6_bar = fsolve(f60,2);# [kmol/hr]
# At the bottom:
# Material Balance:
# Eqn. 9.66:
# (1): L8_bar-GW_bar = W;
# From Fig. 9.28:
yW = 0.035;
x8 = 0.02;
# From Eqn. 9.72:
L8ByGW_bar = (yW-xW)/(x8-xW);
# (2): L8_bar-(L8ByGW_bar*Gw_bar) = 0
a = numpy.array([[1 ,-1],[1 ,-L8ByGW_bar]]);
b = numpy.array([W,0]);
soln = lin.solve(a,b)
L8_bar = soln[0];# [kmol/h]
GW_bar = soln[1];# [kmol/h]
print"The Liquid quantity inside the tower is ",round(L8_bar)," kmol/hr\n"
print"The vapour quantity inside the tower is ",round(GW_bar)," kmol/hr\n"
# The answers are slightly different in textbook due to approximation while in python the answers are precise
Ex9.9: Page 395¶
In [78]:
# Illustration 9.9
# Page: 395
print'Illustration 9.9 - Page: 395\n\n'
# solution
import scipy
import numpy
import numpy.linalg as lin
#****Data****#
P = 695.0;# [kN/square m]
#********#
# a:methanol b:water
# From Illustration 9.8:
Ma = 32.04;# [kg/kmol]
Mb = 18.02;# [kg/kmol]
F = 216.8;# [kmol/h]
Tempo = 19.7;# [OC]
zF = 0.360;# [mole fraction methanol]
HF = 2533;# [kJ/kmol]
D = 84.4;# [kkmol/h]
zD = 0.915;# [mole fraction methanol]
HD = 3640.0;# [kJ/kmol]
Qc = 5990000.0;# [kJ/h]
# Since the bottom will essentially be pure water:
HW = 6094.0;# [kJ/kmol]
# From Steam tables:
Hs = 2699.0;# [enthalpy of saturated steam, kJ/kg]
hW = 4.2*(Tempo-0);# [enthalpy of water, kJ/kg]
HgNpPlus1 = (Hs-hW)*Mb;# [kJ/kmol]
# (1): GNpPlus1-W = D-F [From Eqn. 9.86]
# (2): (GNpPlus1*HgNpPlus1)-(W*HW) = (D*HD)+Qc-(F*HF) [From Eqn. 9.88]
a = numpy.array([[1 ,-1],[HgNpPlus1 ,-HW]]);
b = numpy.array([[D-F],[(D*HD)+Qc-(F*HF)]]);
soln=lin.solve(a,b)
GNpPlus1 = soln[0];# [kmol/h]
W = soln[1];# [kmol/h]
# From Eqn. 9.87:
def f61(xW):
return (F*zF)-((D*zD)+(W*xW))
xW = fsolve(f61,2);
# The enthalpy of the solution at its bubble point is 6048 kJ/kmol, sufficiently closed to 6094 assumed earlier.
# For delta_w:
xdelta_w = W*xW/(W-GNpPlus1);
Q_dprime = ((W*HW)-(GNpPlus1*HgNpPlus1))/(W-GNpPlus1);# [kJ/kmol]
# From Fig 9.27 ad Fig. 9.28, and for the stripping section:
Np = 9.5;
print"Steam Rate: ",round(GNpPlus1,1),"kmol/h\n"
print"Bottom Composition: xW:",round(xW,5),"\n"
print"Number of theoretical stages: ",Np,"\n"
Ex9.10: Page 412¶
In [79]:
# Illustration 9.10
# Page: 412
print'Illustration 9.10 - Page: 412\n\n'
# solution
# a:methanol b:water
Ma = 32.04;# [kg/kmol]
Mb = 18.02;# [kg/kmol]
# Feed:
F1 = 5000;# [kg/h]
F = 216.8;# [kmol/h]
Tempo = 19.7;# [OC]
zF = 0.360;# [mole fraction methanol]
MavF = 23.1;# [kg/kmol]
Tempf = 58.3;# [OC]
# Distillate:
D1 = 2620;# [kg/h]
D = 84.4;# [kkmol/h]
xD = 0.915;# [mole fraction methanol]
# Residue:
R1 = 2380;# [kg/h]
R = 132.4;# [kmol/h]
xW = 0.00565;# [mole fraction methanol]
# From Fig. 9.42 (Pg 413):
BtF = 76.0;# [Bubble point if the feed, OC]
DtF = 89.7;# [Dew point of the feed, OC]
# Latent heat of vaporisation at 76 OC
lambda_a = 1046.7;# [kJ/kg]
lambda_b = 2284;# [kJ/kg]
ha = 2.721;# [kJ/kg K]
hb = 4.187;# [kJ/kg K]
hF = 3.852;# [kJ/kg K]
# If heats of solution is ignaored:
# Enthalpy of the feed at the bubble point referred to the feed temp.
HF = hF*MavF*(BtF-Tempf);# [kJ/kmol]
# enthalpy of the saturated vapour at dew point referred to the liquid at feed temp.
HL = (zF*((ha*Ma*(DtF-Tempf))+(lambda_a*Ma)))+((1-zF)*((hb*Mb*(DtF-Tempf))+(lambda_b*Mb)));# [kJ/kmol]
q = HL/(HL-HF);
slope = q/(q-1);
# In fig. 9.42: xD,xW & zF are located on the 45 degree diagonal & the q line is drawn with slope = 'slope' .
# The operating line for minimum reflux ratio in this case pass through the intersection of the q line and the equilibrium curve.
ordinate = 0.57;
def f62(Rm):
return ordinate-(xD/(Rm+1))
Rm = fsolve(f62,0);# [mole reflux/mole distillate]
# from fig. 9.42 (Pg 413):
# The minimum number of theoretical trays is determied using the 45 degree diagonal as operating line.
Np = 4.9;# [including the reboiler]
R = 1.5*Rm;# [mole reflux/mole distillate]
# From Eqn. 9.49:
L = R*D;# [kmol/h]
# From Eqn. 9.115:
G = D*(R+1);# [kmol/h]
# From Eqn. 9.126:
L_bar = (q*F)+L;# [kmol/h]
# From Eqn. 9.127:
G_bar = (F*(q-1))+G;# [kmol/h]
ordinateN = xD/(R+1);
# As in Fig. 9.43:
# The y-intercept = ordinateN and enriching and exhausting operating lines are plotted.
# Number of theoretical stages are determined.
NpN = 8.8;# [including the reboiler]
print"Number of theoretical stages is ",NpN-1,"\n"
Ex9.11: Page 423¶
In [83]:
# Illustration 9.11
# Page: 423
print'Illustration 9.11 - Page: 423\n\n'
# solution
import math
#****Data****#
# a:ethanol b:water
zF = 0.3;
xa = 0.3;# [mole fraction of ethanol]
Temp = 78.2;# [OC]
Ao = 0.0462;# [Area of perforations,square m]
t = 0.450;# [m]
#******#
Ma = 46.05;# [kg/kmol]
Mb = 18.02;# [kg/kmol]
xb = 1-xa;# [mole fraction of water]
ma = 0.3*Ma/((0.3*Ma)+(xb*Mb));# [mass fraction of ethanol]
mb = 1-ma;# [mass fraction of water]
# Feed:
F1 = 910.0;# [kg/h]
Xa = F1*ma/Ma;# [moles of ethanol]
Xb = F1*mb/Mb;# [moles of water]
F = Xa+Xb;# [Total moles]
# Distillate:
xD = 0.80;# [mole fraction of ethanol]
# If essentially all the ethanol is removed from the residue:
D = Xa/xD;# [kmol/h]
MavD = (xD*Ma)+((1-xD)*Mb);# [kg/kmol]
D1 = D*MavD;# [kg/h]
Density_G = (MavD/22.41)*(273.0/(273+Temp));# [kg/cubic meter]
Density_L = 744.9;# [kg/cubic meter]
sigma = 0.021;# [N/m]
# From Table 6.2,Pg 169:
alpha = (0.0744*t)+0.01173;
beeta = (0.0304*t)+0.015;
At = math.pi*(0.760**2)/4;# [Tower cross sectional Area, square m]
WByT = 530.0/760;# [Table 6.1, Pg 162]
Ad = 0.0808*At;# [Downspout area,square m]
Aa = At-(2*Ad);# [Active area,square m]
# abcissa = (L/G)*(density_G/Density_L)^0.5
# Assume:
abcissa = 0.1;
# From Eqn.6.30:
Cf = (alpha*math.log10(1/abcissa)+beeta)*(sigma/0.020)**0.2;
# From Eqn. 6.29:
Vf = Cf*((Density_L-Density_G)/Density_G)**(1.0/2);# [m/s]
An = At-Ad;# [square m]
R = 3.0;# [Reflux Ratio]
G = D*(R+1);
G1 = (G*22.41/3600)*((273.0+Temp)/273);# [cubic meter/s]
V = G1/An;# [Vapour velocity,m/s]
percent = (V/Vf)*100;
# Vapour velocity is 58 percent of flooding velocity (amply safe)
L = R*D;# [kmol/h]
L1 = L*MavD;# [kg/h]
abcissa = (L1/(G1*3600.0*Density_G))*(Density_G/Density_L)**0.5;
# Since the value of abcissa is less than0.1, the calculaed value of Cf is correct.
# Since the feed is at the buubble point.
q = 1;
# From Eqn. 9.126:
L_bar = L+(q*F);# [kmol/h]
# From Eqn. 9.127:
G_bar = G+F*(q-1);# [kmol/h]
# The enthalpy of saturated steam,referred to 0 OC,69 kN/square m:
HGNpPlus1 = 2699.0;# [kN m/kg]
# This will be the enthalpy as it enters the tower if expanded adiabatically to the tower pressure
# The enthalpy of steam at 1 std. atm:
HGsat = 2676.0;# [kN m/kg]
Lambda = 2257.0;# [kN m/kg]
# From Eqn. 9.140:
def f63(GNpPlus1_bar):
return G_bar-(GNpPlus1_bar*(1+((HGNpPlus1-HGsat)*Mb/(Lambda*Mb))))
GNpPlus1_bar = fsolve(f63,7);
# From Eqn. 9.141:
LNp_bar = L_bar-(G_bar-GNpPlus1_bar);
# Tray Efficiencies:
# Consider the situation:
x = 0.5;
y_star = 0.962;
Temp = 79.8;# [OC]
# This is in the enriching section.
Density_L = 791;# [kg/cubic meter]
Density_G = 1.253;# [kg/cubic meter]
# From equilibrium data:
m = 0.42;
A = L/(m*G);
# From chapter 2:
ScG = 0.930;
Dl = 2.065*10**(-9);# [square m/s]
# For L = 38.73 kmol/h
q = 4.36*10**(-4);# [cubic meter/s]
# For G = 51.64 kmol/h
Va = 1.046;# [m/s]
# From tray dimensions:
z = 0.647;# [m]
Z = 0.542;# [m]
hW = 0.06;# [m]
# From Eqn. 6.61:
NtG = (0.776+(4.57*hW)-(0.238*Va*Density_G**0.5)+(104.6*q/Z))/(ScG**0.5);
# From Eqn. 6.38
hL = 6.10*10**(-3)+(0.725*hW)-(0.238*hW*Va*(Density_G)**0.5)+(1.225*q/z);# [m]
# From Eqn. 6.64:
thetha_L = hL*z*Z/q;# [s]
# From Eqn. 6.62:
NtL = 40000*(Dl**0.5)*((0.213*Va*Density_G**0.5)+0.15)*thetha_L;
# From Eqn. 6.52:
NtoG = 1/((1/NtG)+(1/(A*NtL)));
# From Eqn. 6.51:
EoG = 1-math.exp(-NtoG);
# From Eqn. 6.63:
DE = ((3.93*10**(-3))+(0.0171*Va)+(3.67*q/Z)+(0.1800*hW))**2;
# From Eqn. 6.59:
Pe = Z**2/(DE*thetha_L);
# From Eqn. 6.58:
eta = (Pe/2)*((1+(4*m*G1*EoG/(L1*Pe)))**0.5-1);
# From Eqn. 6.57:
EMG = EoG*(((1-math.exp(-(eta+Pe)))/((eta+Pe)*(1+(eta+Pe)/eta)))+((math.exp(eta)-1)/(eta*(1+(eta/(eta+Pe))))));
# Entrainment is neglible:
# Similarly for other x
# Value = [x Entrainment]
#Value = [0 0.48;0.1 .543;0.3 0.74;0.5 EMG;0.7 0.72];
# Tray Calculation:
op_intercept = xD/(R+1);
# From Fig. 9.48:
# The exhausting section operating line, on this scale plot, for all practical purposes passes through the origin.
# The broken curve is located so that, at each concentration, vertical distances corresponding to lines BC and AC are in the ratio of EMG.
# This curve is used instead of equilibrium trays to locate the ideal trays.
# The feed tray is thirteenth.
x14 = 0.0150;
alpha = 8.95;
EMG = 0.48;
A_bar = L_bar/(alpha*G_bar);
# From Eqn. 8.16:
Eo = math.log(1+(EMG*((1/A_bar)-1)))/math.log(1/A_bar);
# The 6 real trays corresponds to:
NRp = 6*Eo;
xW = 0.015/((math.exp(NRp*math.log(1/A_bar))-A_bar)/(1-A_bar));# [mole fraction ethanol]
# This corresponds to ethanol loss of 0.5 kg/day.
print"The mole fraction of ethanol in residue is",round(xW,8)
print"The Reflux ratio of ",R," will cause the ethanol loss of 0.5 kg/day\n"
print"Larger reflux ratios would reduce this, but the cost of additional steam will probaby make them not worthwile.\n"
print"Smaller values of R, with corresponding reduced steam cost and larger ethanol loss, should be considered, but care must be taken to ensure vapour velocities above the weeping velocities."
Ex-9.12: Pg- 429¶
In [3]:
# Illustration 9.12
# Page: 429
print'Illustration 9.12 - Page: 429\n\n'
# solution
import math
import numpy
import matplotlib.pyplot as plt
%matplotlib inline
# a:methanol b:water
# Vapour and liquid quantities throughout the tower, as in Illustration 9.8, with the Eqn. 9.62, 9.64, 9.72, 9.74:
# Data = [x tL(OC) y tG(OC) Vapor(kmol/h) Vapor(kg/h) Liquid(kmol/h) Liquid(kg/h)]
Ma = 34.02;# [kg/kmol]
Mb = 18.02;# [kg/kmol]
Temp = 78.7;# [OC]
x = numpy.array([0.915, 0.600 ,0.370, 0.370, 0.200, 0.100, 0.02]);
y = numpy.array([0.915, 0.762, 0.656, 0.656, 0.360 ,0.178, 0.032]);
plt.plot(x,y);
plt.grid('on');
ax = pylab.gca()
ax.set_xlabel("mole fraction of methanol in liquid");
ax.set_ylabel("mole fraction of methanol in vapour");
plt.title("Operating Line curve");
plt.legend(loc="lower right")
plt.show()
#x = 0.370: the dividing point between stripping and enriching section
tL =numpy.array([66, 71, 76, 76, 82, 87, 96.3]);# [Bubble point, OC]
tG = numpy.array([68.2 ,74.3 ,78.7 ,78.7 ,89.7 ,94.7 ,99.3]);# [Dew Point, OC]
Vapor = numpy.array([171.3, 164.0 ,160.9, 168.6, 161.6, 160.6, 127.6]);# [kmol/h]
Vapor1 = numpy.array([5303, 4684, 4378, 4585, 3721, 3296 ,2360]);# [kg/h]
Liquid = numpy.array([86.7 ,79.6 ,76.5 ,301, 294, 293, 260]);# [kmol/h]
Liquid1 = numpy.array([2723, 2104, 1779 ,7000, 6138, 5690 ,4767]);# [kg/h]
Data = numpy.zeros(shape=(7,8));
for j in range(1,7):
Data[j,0]= x[j];
Data[j,1]= tL[j];
Data[j,2]= y[j];
Data[j,3]= tG[j];
Data[j,4]= Vapor[j];
Data[j,5]= Vapor1[j];
Data[j,6]= Liquid[j];
Data[j,7]= Liquid1[j];
# The tower diameter will be set by the conditions at the top of the stripping section because of the large liquid flow at this point.
# From Illustration 9.8:
G = Data[3,5];
L = Data[3,7];
Density_G = (Data[3,5]/(22.41*Data[3,4]))*(273.0/(273+Temp));# [kg/cubic m]
Density_L = 905.0;# [kg/cubic m]
# abcissa = (L/G)*(Density_L/Density_G)^0.5
abcissa = (Data[3,7]/Data[3,5])*(Density_G/Density_L)**0.5;
# From Fig. 6.34, choose a gas pressure drop of 450 N/square m/m
ordinate = 0.0825;
# From Table 6.3 (Pg 196):
Cf = 95;
viscosity_L = 4.5*10**(-4);# [kg/m.s]
sigma = 0.029;# [N/m]
J = 1;
G_prime = (ordinate*Density_G*(Density_L-Density_G)/(Cf*viscosity_L**0.1))**0.5;# [kg/square m.s]
A = G/(3600*G_prime);# [Tower ,cross section area,square m]
L_prime = L/(A*3600);# [kg/square m.s]
# Mass transfer will be computed for the same location:
# From Table 6.4 (Pg 205):
m = 36.4;
n = (0.0498*L_prime)-0.1013;
p = 0.274;
aAW = m*((808*G_prime/Density_G**0.5)**n)*L_prime**p;# [square m/cubic m]
# From Table 6.5 (Pg 206):
dS = 0.0530;# [m]
beeta = 1.508*dS**0.376;
shi_LsW = 2.47*10**(-4)/dS**1.21;
shi_LtW = ((2.09*10**(-6))*(737.5*L_prime)**beeta)/dS**2;
shi_LOW = shi_LtW-shi_LsW;
shi_Ls = (0.0486*viscosity_L**0.02*sigma**0.99)/(dS**1.21*Density_L**0.37);
H = ((975.7*L_prime**0.57*viscosity_L**0.13)/(Density_L**0.84*((2.024*L_prime**0.430)-1)))*(sigma/0.073)**(0.1737-0.262*math.log10(L_prime));# [m]
shi_Lo = shi_LOW*H;
shi_Lt = shi_Lo+shi_Ls;
# From Eqn. 6.73:
aA = aAW*(shi_Lo/shi_LOW);# [square m/cubic m]
# From Table 6.3 (Pg 196):
e = 0.71;
# From Eqn. 6.71:
eLo = e-shi_Lt;
# From Chapter 2:
ScG = 1;
MavG = 0.656*Ma+(1-0.656)*Mb;# [kg/kmol]
G = G_prime/MavG;
viscosity_G = 2.96*10**(-5);# [kg/m.s]
# From Eqn. 6.70:
Fg = (1.195*G/ScG**(2/3))*((dS*G_prime/(viscosity_G*(1-eLo)))**(-0.36));# [kmol/square m s (mole fraction)]
kY_prime = Fg;# [kmol/square m s (mole fraction)]
DL = 4.80*10**(-9);# [square m/s]
ScL = viscosity_L/(Density_L*DL);
# From Eqn. 6.72:
kL = (25.1*DL/dS)*((dS*L_prime/viscosity_L)**0.45)*ScL**0.5;# [kmol/square m s (kmol/cubic m)]
# At 588.33 OC
Density_W = 53.82;# [kg/cubic m]
kx_prime = Density_W*kL;# [kmol/square m s (mole fraction)]
# Value1 = [x G a ky_prime*10^3 kx_prime]
Value1 = numpy.array([[0.915 ,0.0474 ,20.18 ,1.525, 0.01055],[0.6, 0.0454 ,21.56 ,1.542, 0.00865],[0.370 ,0.0444 ,21.92 ,1.545 ,0.00776],[0.370, 0.0466 ,38, 1.640, 0.0143],[0.2 ,0.0447, 32.82 ,1.692 ,0.0149],[0.1 ,0.0443 ,31.99 ,1.766 ,0.0146],[0.02, 0.0352 ,22.25 ,1.586 ,0.0150]]);
# From Fig: 9.50
# At x = 0.2:
y = 0.36;
slope = -(Value1[4,4]/(Value1[4,3]*10**(-3)));
# The operating line drawn from(x,y) with slope. The point where it cuts the eqb. line gives yi.
# K = ky_prime*a(yi-y)
# For the enriching section:
# En = [y yi 1/K Gy]
En = numpy.array([[0.915 ,0.960, 634 ,0.0433],[0.85 ,0.906 ,532.8 ,0.0394],[0.8 ,0.862 ,481.1 ,0.0366],[0.70, 0.760 ,499.1, 0.0314],[0.656, 0.702, 786.9, 0.0292]]);
# For the Stripping section:
# St = [y yi 1/K Gy]
St = numpy.array([[0.656, 0.707, 314.7, 0.0306],[0.50, 0.639, 124.6 ,0.0225],[0.40 ,0.580, 99.6 ,0.01787],[0.3 ,0.5 ,89 ,0.0134],[0.2 ,0.390, 92.6 ,0.00888],[0.10, 0.232, 154.5, 0.00416],[0.032 ,0.091, 481 ,0.00124]])
# Graphical Integration, according to Eqn.9.52::
plt.plot(En[:,3],En[:,2],'g');
plt.grid();
ax = pylab.gca()
ax.set_xlabel("Gy");
ax.set_ylabel("1 / (ky_prime*a*(yi-y))");
plt.title("Graphical Integration for Enriching section");
plt.show()
# From Area under the curve:
Ze = 7.53;# [m]
# Graphical Integration:
plt.plot(St[:,3],St[:,2],'r');
plt.grid('on');
ax = pylab.gca()
ax.set_xlabel("Gy");
ax.set_ylabel("1 / (ky_prime*a*(yi-y))");
plt.title("Graphical Integration for Stripping section");
plt.show()
# From Area under the curve:
Zs = 4.54;# [m]
# Since the equlibrium curve slope varies so greatly that the use of overall mass transfer coeffecient is not recommended:
print"Height of Tower for enriching Section is ",Ze," m\n"
print"Height of Tower for Stripping Section is ",Zs," m\n"
Ex9.13: Page 436¶
In [85]:
# Illustration 9.13:
print'Illustration 9.13\n\n'
#**************************Calculation Of Minimum Reflux ratio************************#
# Page: 436
print'Page: 436\n\n'
import math
import numpy
from scipy import interp
from scipy.optimize import fsolve
import numpy.linalg as lin
#***Data***#
# C1:CH4 C2:C2H6 C3:n-C3H8 C4:n-C4H10 C5:n-C5H12 C6:n-C6H14
# zF = [zF(C1) zF(C2) zF(C3) zF(C4) zF(C5) zF(C6)]
zF = numpy.array([0.03 ,0.07 ,0.15 ,0.33 ,0.30 ,0.12]);# [mole fraction]
LF_By_F = 0.667;
Temp = 82;# [OC]
ylk = 0.98;
yhk = 0.01;
#**********#
# Data = [m HG HL(30 OC);m HG HL(60 OC);m HG HL(90 OC);m HG HL(120 OC);]
Data1 = numpy.array([[16.1 ,12790 ,9770],[19.3 ,13910, 11160],[21.8 ,15000, 12790],[24.0 ,16240, 14370]]);# [For C1]
Data2 = numpy.array([[3.45, 22440, 16280],[4.90 ,24300 ,18140],[6.25 ,26240 ,19890],[8.15 ,28140, 21630]]);# [For C2]
Data3 = numpy.array([[1.10, 31170, 16510],[2.00 ,33000 ,20590],[2.90, 35800 ,25600],[4.00 ,39000, 30900]]);# [For C3]
Data4 = numpy.array([[0.35, 41200 ,20350],[0.70 ,43850 ,25120],[1.16 ,46500, 30000],[1.78 ,50400 ,35400]]);# [For C4]
Data5 = numpy.array([[0.085, 50500, 24200],[0.26, 54000 ,32450],[0.50 ,57800 ,35600],[0.84, 61200 ,41400]]);# [For C5]
Data6 = numpy.array([[0.0300, 58800 ,27700],[0.130, 63500, 34200],[0.239 ,68150 ,40900],[0.448, 72700 ,48150]]);# [For C6]
# T = [Temparature]
T = numpy.array([30,60,0,120]);
# Flash vaporisation of the Feed:
# Basis: 1 kmol feed throughout
# After Several trials, assume:
F = 1.0;# [kmol]
GF_By_F = 0.333;
LF_By_GF = LF_By_F/GF_By_F;
m82 = numpy.zeros(6);
y = numpy.zeros(6);
m82[0] = interp(Temp,T,Data1[:,1]);# [For C1]
m82[1] = interp(Temp,T,Data2[:,0]);# [For C2]
m82[2] = interp(Temp,T,Data3[:,0]);# [For C3]
m82[3] = interp(Temp,T,Data4[:,0]);# [For C4]
m82[4] = interp(Temp,T,Data5[:,0]);# [For C5]
m82[5] = interp(Temp,T,Data6[:,0]);# [For C6]
for i in range (0,6):
y[i] = zF[i]*(LF_By_GF+1)/(1.0+(2/m82[i]));
Sum = sum(y);
# Since Sum is sufficiently close to 1.0, therefore:
q = 0.67;# [LF_By_F]
# Assume:
# C3: light key
# C5: heavy key
zlkF = zF[2];# [mole fraction]
zhkF = zF[4];# [mole fraction]
ylkD = ylk*zF[2];# [kmol]
yhkD = yhk*zF[4];# [kmol]
# Estimate average Temp to be 80 OC
m80 = numpy.zeros(6);
alpha80 = numpy.zeros(6);
m80[0] = interp(Temp,T,Data1[:,0]);# [For C1]
m80[1] = interp(Temp,T,Data2[:,0]);# [For C2]
m80[2] = interp(Temp,T,Data3[:,0]);# [For C3]
m80[3] = interp(Temp,T,Data4[:,0]);# [For C4]
m80[4] = interp(Temp,T,Data5[:,0]);# [For C5]
m80[5] = interp(Temp,T,Data6[:,0]);# [For C6]
for i in range(0,6):
alpha80[i] = m80[i]/m80[4];
# By Eqn. 9.164:
yD_By_zF1 = (((alpha80[0]-1)/(alpha80[2]-1))*(ylkD/zF[2]))+(((alpha80[2]-alpha80[0])/(alpha80[2]-1))*(yhkD/zF[4]));# [For C1]
yD_By_zF2 = (((alpha80[1]-1)/(alpha80[2]-1))*(ylkD/zF[2]))+(((alpha80[2]-alpha80[1])/(alpha80[2]-1))*(yhkD/zF[4]));# [For C2]
yD_By_zF6 = (((alpha80[5]-1)/(alpha80[2]-1))*(ylkD/zF[2]))+(((alpha80[2]-alpha80[5])/(alpha80[2]-1))*(yhkD/zF[4]));# [For C6]
# The distillate contains:
yC1 = 0.03;# [kmol C1]
yC2 = 0.07;# [kmol C2]
yC6 = 0;# [kmol C6]
# By Eqn 9.165:
def g1(phi):
return (((alpha80[0]*zF[0])/(alpha80[0]-phi))+((alpha80[1]*zF[1])/(alpha80[1]-phi))+((alpha80[2]*zF[2])/(alpha80[2]-phi))+((alpha80[3]*zF[3])/(alpha80[3]-phi))+((alpha80[4]*zF[4])/(alpha80[4]-phi))+((alpha80[5]*zF[5])/(alpha80[5]-phi)))-(F*(1-q))
# Between alphaC3 & alphaC4:
phi1 = fsolve(g1,3);
# Between alphaC4 & alphaC5:
phi2 = fsolve(g1,1.5);
# From Eqn. 9.166:
# Val = D*(Rm+1)
# (alpha80(1)*yC1/(alpha80(1)-phi1))+(alpha80(2)*yC2/(alpha80(2)-phi1))+(alpha80(3)*ylkD/(alpha80(3)-phi1))+(alpha80(4)*yD/(alpha80(4)-phi1))+(alpha80(i)*yhkD/(alpha80(5)-phi1))+(alpha80(6)*yC6/(alpha80(6)-phi1)) = Val.....................(1)
# (alpha80(1)*yC1/(alpha80(1)-phi2))+(alpha80(2)*yC2/(alpha80(2)-phi2))+(alpha80(3)*ylkD/(alpha80(3)-phi2))+(alpha80(4)*yD/(alpha80(4)-phi2))+(alpha80(i)*yhkD/(alpha80(5)-phi2))+(alpha80(6)*yC6/(alpha80(6)-phi2)) = Val ....................(2)
# Solving simultaneously:
a =numpy.array([[-alpha80[3]/(alpha80[3]-phi1), 1],[-alpha80[3]/(alpha80[3]-phi2), 1]]);
b =numpy.array([[alpha80[0]*yC1/[alpha80[0]-phi1]]+[alpha80[1]*yC2/[alpha80[1]-phi1]]+[alpha80[2]*ylkD/[alpha80[2]-phi1]]+[alpha80[i]*yhkD/[alpha80[4]-phi1]]+[alpha80[5]*yC6/[alpha80[5]-phi1]],[alpha80[0]*yC1/[alpha80[0]-phi2]]+[alpha80[1]*yC2/[alpha80[1]-phi2]]+[alpha80[2]*ylkD/[alpha80[2]-phi2]]+[alpha80[i]*yhkD/[alpha80[4]-phi2]]+[alpha80[5]*yC6/[alpha80[5]-phi2]]])
soln = lin.solve(a,b);
yC4 =0.1313547 # [kmol C4 in the distillate]
Val =0.617469; # [kmol C4 in the distillate]
# For the distillate, at a dew point of 46 OC
ydD = numpy.array([yC1,yC2 ,ylkD ,yC4 ,yhkD ,yC6]);
D = sum(ydD);
yD = zeros(6);
m46 = zeros(6);
alpha46 = zeros(6);
Ratio1= zeros(6);
m46[0] = interp(Temp,T,Data1[:,0]);# [For C1]
m46[1] = interp(Temp,T,Data2[:,0]);# [For C2]
m46[2] = interp(Temp,T,Data3[:,0]);# [For C3]
m46[3] = interp(Temp,T,Data4[:,0]);# [For C4]
m46[4] = interp(Temp,T,Data5[:,0]);# [For C5]
m46[5] = interp(Temp,T,Data6[:,0]);# [For C6]
yD=numpy.array([0.0786,0.1835,0.3854,0.34,0.007866,0.0])
# mhk = mC5 at 46.6 OC, the assumed 46 OC is satisfactory.
# For the residue, at a dew point of 46 OC
xwW =numpy.array([zF[0]-yC1, zF[1]-yC2 ,zF[2]-ylkD, zF[3]-yC4, zF[4]-yhkD, zF[5]-yC6]);
W = sum(xwW);
xW = zeros(6);
m113 = zeros(6);
alpha113 = zeros(6);
alphalk_av=zeros(6);
alpha_av=zeros(6);
Value=zeros(6);
m113[0] = interp(Temp,T,Data1[:,1]);# [For C1]
m113[1] = interp(Temp,T,Data2[:,1]);# [For C2]
m113[2] = interp(Temp,T,Data3[:,1]);# [For C3]
m113[3] = interp(Temp,T,Data4[:,1]);# [For C4]
m113[4] = interp(Temp,T,Data5[:,1]);# [For C5]
m113[5] = interp(Temp,T,Data6[:,1]);# [For C6]
for i in range(0,6):
alpha113[i] = m113[i]/m113[4];
xW[i] = xwW[i]/W;
# Ratio = yD/alpha46
Value[i] = alpha113[i]*xW[i];
# mhk = mC5 at 114 OC, the assumed 113 OC is satisfactory.
Temp_Avg = (114+46.6)/2;# [OC]
# Temp_avg is very close to the assumed 80 OC
Rm = (Val/D)-1;
print"Minimum Reflux Ratio is ",Rm," mol reflux/mol distillate\n \n"
print"*****************Distillate Composition*********************\n"
print"C1\t \t \t \t:",yD[0]
print"C2\t \t \t \t:",yD[1]
print"C3\t \t \t \t:",yD[2]
print"C4\t \t \t \t:",yD[3]
print"C5\t \t \t \t:",yD[4]
print"C6\t \t \t \t:",yD[5]
print"\n"
print"*****************Residue Composition*********************\n"
print"C1\t \t \t \t: ",xW[0]
print"C2\t \t \t \t: ",xW[1]
print"C3\t \t \t \t: ",xW[2]
print"C4\t \t \t \t: ",xW[3]
print"C5\t \t \t \t: ",xW[4]
print"C6\t \t \t \t: ",xW[5]
print"\n"
#**********************Number of Theoretical stage***********************#
# Page:440
print'Page: 440\n\n'
for i in range(0,6):
alpha_av[i] = (alpha46[i]*alpha113[i])**0.5;
alphalk_av = alpha_av[1];
# By Eqn. 9.167:
xhkW = xwW[3];
xlkW = xwW[1];
Nm = 3.496;
# Ratio = yD/xW
Ratio2= zeros(6)
for i in range(0,6):
Ratio2[i] = (alpha_av[i]**(Nm+1))*yhkD/xhkW;
# For C1:
# yC1D-Ratio(1)*xC1W = 0
# yC1D+xC1W = zF(1)
# Similarly for others
yD2=zeros(6)
xW2=zeros(6)
for i in range(0,6):
a = numpy.array([[1 ,-Ratio2[i]],[1, 1]]);
b = [0,zF[i]];
soln =lin.solve(a,b);
yD2[i] = soln[0];# [kmol]
xW2[i] = soln[1];# [kmol]
D = sum(yD2);# [kmol]
W = sum(xW2);# [kmol]
# The distillate dew point computes to 46.6 OC and the residue bubble point computes to 113 OC, which is significantly close to the assumed.
#***************Product composition at R = 0.8***********************#
# Page:441
print'Page: 441\n\n'
# Since C1 and C2 do not enter in the residue nor C6 in the distillate, appreciably at total reflux or minimum reflux ratio, it will be assumed that they will not enter R = 0.8. C3 and C5 distribution are fixed by specifications. Only that C4 remains to be estimated.
# R = [Infinte 0.8 0.58] [Reflux ratios For C4]
R = [inf ,0.8, 0.58];
# Val = R/(R+1)
val=[ 0 , 2.0 , 2.0]
# ydD = [Inf 0.58]
y4D = [0.1255, 0.1306];
yC4D = 0.1306 ;# by Linear Interpolation
# For Distillate:
Sum1 = sum(Ratio1);
x0 = numpy.array([0.004,0.0444501,0.2495,0.65640,0.0451,0.0])
print"For the reflux ratio of 0.8\n"
print"*****************Distillate Composition*********************\n"
print"\t\t\t Liquid reflux in equilibrium with the distillate vapour\n"
for i in range(0,6):
print"C",i,"\t \t \t \t\t:",x0[i]
# For boiler:
#**********Number Of Theoretical Trays***************#
# Page: 443
print'Page: 443\n\n'
R = 0.8;# [reflux ratio]
# From Eqn. 9.175
intersection = (zlkF-(ylkD/D)*(1-q)/(R+1))/(zhkF-(yhkD/D)*(1-q)/(R+1));
# Enriching Section:
y1 = zeros(5);
L = R*D;# [kmol]
G = L+D;# [kmol]
# Assume: Temp1 = 57 OC
# alpha57 = [C1 C2 C3 C4 C5]
alpha57 = numpy.array([79.1 ,19.6 ,7.50, 2.66, 1]);
# From Eqn. 9.177, n = 0:
Val57=zeros(6)
for i in range(0,5):
y1[i] = (L/G)*x0[i]+((D/G)*yD[i]);
Val57[i] = y1[i]/alpha57[i];
x1 = Val57/sum(Val57);
mC5 = sum(Val57);
Temp1 = 58.4; # [OC]
# Liquid x1's is in equilibrium with y1's.
xlk_By_xhk1 = x1[2]/x1[4];
# Tray 1 is not the feed tray.
# Assume: Temp2 = 63 OC
# alpha63 = [C1 C2 C3 C4 C5]
alpha63 = numpy.array([68.9 ,17.85, 6.95, 2.53, 1.00]);
# From Eqn. 9.177, n = 1:
y2=zeros(6)
Val63=zeros(6)
for i in range(0,5):
y2[i] = (L/G)*x1[i]+((D/G)*yD[i]);
Val63[i] = y1[i]/alpha63[i];
mC5 = sum(Val63);
x2 = Val63/sum(Val63);
xlk_By_xhk2 = x2[2]/x2[4];
# The tray calculation are continued downward in this manner.
# Results for trays 5 & 6 are:
# Temp 75.4 [OC]
# x5 = [C1 C2 C3 C4 C5]
x5 = numpy.array([0.00240, 0.0195, 0.1125, 0.4800, 0.3859]);
xlk_By_xhk5 = x5[2]/x5[4];
# Temp6 = 79.2 OC
# x6 = [C1 C2 C3 C4 C5]
x6 = numpy.array([0.00204 ,0.0187 ,0.1045, 0.4247 ,0.4500]);
xlk_By_xhk6 = x6[2]/x6[4];
# From Eqn. 9.176:
# Tray 6 is the feed tray
Np1 = 6;
# Exhausting section:
# Assume Temp = 110 OC
L_bar = L+(q*F);# [kmol]
G_bar = L_bar-W;# [kmol]
# alpha57 = [C3 C4 C5 C6]
alpha110 = numpy.array([5 ,2.2 ,1, 0.501]);
# From Eqn. 9.178:
xNp = zeros(4);
Val110=zeros(6)
k = 0;
for i in range(2,6):
xNp[k] = ((G_bar/L_bar)*yNpPlus1[i])+((W/L_bar)*xW[i]);
Val110[k] = alpha110[k]*xNp[k];
k = k+1;
yNp = Val110/sum(Val110);
mC5 = 1/sum(Val110);
# yNp is in Eqb. with xNp:
xlk_By_xhkNp = xNp[0]/xNp[3];
# Results for Np-7 to Np-9 trays:
# For Np-7
# Temp = 95.7 OC
# xNpMinus7 = [C3 C4 C5 C6]
xNpMinus7 = numpy.array([0.0790 ,0.3944 ,0.3850, 0.1366]);
xlk_By_xhkNpMinus7 = xNpMinus7[0]/xNpMinus7[2];
# For Np-8
# Temp = 94.1 OC
# xNpMinus8 = [C3 C4 C5 C6]
xNpMinus8 = numpy.array([0.0915, 0.3897 ,0.3826, 0.1362]);
xlk_By_xhkNpMinus8 = xNpMinus8[0]/xNpMinus8[2];
# For Np-9
# Temp = 93.6 OC
# xNpMinus9 = [C3 C4 C5 C6]
xNpMinus9 = numpy.array([0.1032, 0.3812, 0.3801 ,0.1355]);
xlk_By_xhkNpMinus9 = xNpMinus9[0]/xNpMinus9[2];
# From Eqn. 9.176:
# Np-8 is the feed tray.
def g2(Np):
return Np-8-Np1
Np = fsolve(g2,7);
print"Number of theoretical Trays required for R = 0.8: ",Np[0]
print"\n"
#**************Composition Correction*****************#
# Page: 446
print'Page: 446\n\n'
# New Bubble Point:
# Temp = 86.4 OC
x6_new = x6*(1-xNpMinus8[3]);
x6_new[4] = xNpMinus8[3];
# alpha86 = [C1 C2 C3 C4 C5 C6]
alpha86 =numpy.array([46.5, 13.5, 5.87, 2.39, 1.00, 0.467]);
# From Eqn. 9.181:
xhkn = x5[3];
xhknPlus1 = x6_new[3];
xC65 = alpha86[5]*x6_new[4]*xhkn/xhknPlus1;
x5_new = x5*(1-xC65);
x5_new[4] = 1-sum(x5_new);
# Tray 5 has a bubble point of 80 OC
# Similarly , the calculations are continued upward:
# x2_new = [C1 C2 C3 C4 C5 C6]
x2_new = numpy.array([0.0021, 0.0214 ,0.1418, 0.6786, 0.1553, 0.00262]);
# y2_new = [C1 C2 C3 C4 C5 C6]
y2_new = numpy.array([0.0444, 0.111 ,0.2885, 0.5099, 0.0458 ,0.00034]);
# x1_new = [C1 C2 C3 C4 C5 C6]
x1_new = numpy.array([0.00226, 0.0241, 0.1697 ,0.7100, 0.0932, 0.00079]);
# y1_new = [C1 C2 C3 C4 C5 C6]
y1_new = numpy.array([0.0451 ,0.1209 ,0.3259 ,0.4840 ,0.0239 ,0.000090]);
# x0_new = [C1 C2 C3 C4 C5 C6]
x0_new = numpy.array([0.00425 ,0.0425 ,0.2495, 0.6611 ,0.0425 ,0.00015]);
# yD_new = [C1 C2 C3 C4 C5 C6]
yD_new = numpy.array([0.0789 ,0.1842 ,0.3870 ,0.3420 ,0.0079, 0.00001]);
# From Eqn. 9.184:
# For C1 & C2
alphalkm = alpha86[2];
xlkmPlus1 = xNpMinus7[0];
xlkm = x6_new[2];
xC17 = x6_new[0]*alpha86[2]*xlkmPlus1/(alpha86[0]*xlkm);
xC27 = x6_new[1]*alpha86[2]*xlkmPlus1/(alpha86[1]*xlkm);
# Since xC17 = 1-xC27
# The adjusted value above constitute x7's.
# The new bubbl point is 94 OC
# The calculations are continued down in the same fashion.
# The new tray 6 has:
# xC1 = 0.000023 & xC2 = 0.00236
# It is clear that the conc. of these components are reducing so rapidly that there is no need to go an further.
print"******Corrected Composition***********\n"
print"Component\t \tx2\t \t \t y2\t \t \t x1\t \t \t y1\t \t \tx0\t \t \tyD\n"
for i in range(0,6):
print"C",i,"\t \t \t",x2_new[i],"\t \t \t \t ",y2_new[i],"\t \t \t \t ",x1_new[i],"\t \t \t \t",y1_new[i],"\t \t \t \t \t",x0_new[i],"\t \t \t \t",yD_new[i]
print"\n"
#*************Heat Load of Condensor & Boiler & L/G ratio**********#
# Page 448
print'Page: 448\n\n'
# Values of x0, yD & y1 are taken from the corrected concentration.
# HD46 = [C1 C2 C3 C4 C5 C6]
HD46 = numpy.array([13490, 23380, 32100, 42330, 52570, 61480]);# [kJ/kmol]
yDHD= zeros(6)
for i in range(0,6):
yDHD[i] = yD_new[i]*HD46[i];
HD = sum(yDHD);# [kJ]
# HL46 = [C1 C2 C3 C4 C5 C6]
HL46 = numpy.array([10470, 17210, 18610, 22790, 27100, 31050]);# [kJ/kmol]
xHL=zeros(6)
for i in range(0,6):
xHL[i] = x0_new[i]*HL46[i];
HL0 = sum(xHL);# [kJ]
# HG58 = [C1 C2 C3 C4 C5 C6]
HG58 = numpy.array([13960, 24190, 37260, 43500, 53900, 63500]);# [kJ/kmol]
yHG1=zeros(6)
for i in range(0,6):
yHG1[i] = y1_new[i]*HG58[i];
HG1 = sum(yHG1);# [kJ]
# From Eqn. 9.54:
Qc = D*((R+1)*HG1-(R*HL0)-HD);# [kJ/kmol feed]
# Similarly:
HW = 39220;# [kJ]
HF = 34260;# [kJ]
# From Eqn. 9.55:
Qb = (D*HD)+(W*HW)+Qc-(F*HF);# [kJ/kmol feed]
# For tray n = 1
G1 = D*(R+1);# [kmol]
# With x1 & y2 from corrected composition;
# HG66 = [C1 C2 C3 C4 C5 C6]
HG66 = numpy.array([14070, 24610, 33800, 44100, 54780, 64430]);# [kJ/kmol feed]
yHG2=zeros(6)
for i in range(0,6):
yHG2[i] = y2_new[i]*HG66[i];
HG2 = sum(yHG2);# [kJ]
# HL58 = [C1 C2 C3 C4 C5 C6]
HL58 =numpy.array([11610 ,17910 ,20470, 24900, 29500, 33840]);# [kJ/kmol feed]
xHL1=zeros(6)
for i in range(0,6):
xHL1[i] = x1_new[i]*HL58[i];
HL1 = sum(xHL1);# [kJ]
# From Eqn. 9.185:
G2 = (Qc+D*(HD-HL1))/(HG2-HL1);# [kmol]
L2 = G2-D;# [kmol]
L2_By_G2 = L2/G2;
# Similarly, the calculations are made for other trays in enriching section.
# For tray, Np = 14:
# C1 & C2 are absent.
# HG113 = [C3 C4 C5 C6]
HG113 = numpy.array([38260, 49310 ,60240, 71640]);# [kJ/kmol feed]
k = 2;
yHG15=zeros(6)
for i in range(0,4):
yHG15[i] = yNpPlus1[k]*HG113[i];
k = k+1;
HG15 = sum(yHG15);
# HL107 = [C3 C4 C5 C6]
HL107 = numpy.array([29310 ,31870, 37680 ,43500]);# [kJ/kmol feed]
xHL14=zeros(6)
for i in range(0,4):
xHL14[i] = xNp[i]*HL107[i];
HL14 = sum(xHL14);# [kJ]
# Similarly:
HL13 = 36790;# [kJ]
HG14 = 52610;# [kJ]
# From Eqn. 9.186:
G15_bar = (Qb+(W*(HL14-HW)))/(HG15-HL14);# [kmol]
L14_bar = W+G15_bar;# [kmol]
G14_bar = (Qb+(W*(HL13-HW)))/(HG14-HL13);# [kmol]
L14_By_G14 = L14_bar/G14_bar;
print"Condensor Heat Load kJ:\n",HL0
print"Reboiler Heat Load kJ:\n",HG15
# For other Exhausting Section Trays:
# Result = [Tray No. L_By_G Temp(OC)]
# Tray 0: Condensor
# Tray 15: Reboiler
Result = numpy.array([[0,0.80 ,46.6],[1 ,0.432 ,58.4],[2, 0.437, 66],[3, 0.369, 70.4],[4 ,0.305, 74],[5 ,0.310, 80.3],[6, 1.53, 86.4],[7, 4.05 ,94.1],[8 ,3.25 ,96.3],[9, 2.88 ,97.7],[10 ,2.58 ,99],[11, 2.48 ,100],[12 ,2.47 ,102.9],[13 ,2.42 ,104.6],[14 ,2.18 ,107.9],[15, 1.73 ,113.5]]);
print"**************L/G*************\n"
print"Tray No. \t\t L/G\t\t\t\t Temp(OC)\n"
for i in range(0,16):
print Result[i,0],"\t\t \t \t",Result[i,1],"\t \t \t",Result[i,2];
# These values are not final.
# They scatter eratically because they are based on the temp. and conc. computed with the assumption of constant L/G
print"\n"
#**************Thiele Geddes Method******************#
# Page:452
print'Page: 452\n\n'
# Use the tray Temperature to obtain m.
# For C4:
# m = [0(Condensor) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15(Reboiler)]
m = numpy.array([0.50 ,0.66, 0.75 ,0.81 ,0.86 ,0.95 ,1.07 ,1.22 ,1.27 ,1.29 ,1.30, 1.32, 1.40, 1.45, 1.51, 1.65]);
A = numpy.array([1.6,0.65,0.582,0.4555,0.354,0.326,1.42990])
S = numpy.array([0.3012,0.39076,0.4479,0.503875,0.53225,0.56680,0.59917,0.69,0.95375])
# f = Tray No. 6
# From Eqn. 9.196:
# Value1 = Gf*yf/(D*zD)
Sum = 0;
for i in range(0,6):
Val = 1;
for j in range(0,6):
Val = Val*A[j];
Sum = Sum+Val;
Value1 = 1+Sum;
# From Eqn. 9.206:
# Value2 = Lf_bar*xf/(W*xW);
Sum = 0.5316
Value2 = 1+Sum;
# From Eqn. 9.208:
# Value3 = W*xW/(D*zD)
Value3 = A[6]*Value1/Value2;
# From Eqn. 9.210:
DyD = F*zF[3]/(Value3+1);# [kmol,C4]
# From Eqn. 9.209:
WxW = ((F*zF[3]))-(DyD);# [kmol, C4]
# Similarly:
# For [C1; C2; C3; C4; C5; C6]
# Result2 = [Value1 Value2 Value3 DyD WxW]
Result2 = numpy.array([[1.0150, 254*10**6 ,288*10**(-10), 0.03, 0],[1.0567, 8750, 298*10**(-5) ,0.07 ,0],[1.440, 17.241 ,0.0376 ,0.1447, 0.0053],[1.5778 ,1.5306 ,1.475, 0.1335 ,0.1965],[15580, 1.1595, 45.7 ,0.00643 ,0.29357],[1080 ,1.0687 ,7230 ,0.0000166 ,0.1198]]);
D = sum(Result2[:,2]);# [kmol]
W = sum(Result2[:,3]);# [kmol]
# In the Distillate:
DyD_C3 = Result[1,2];# [kmol]
zFC3 = zF[2];# [kmol]
percentC3 = (DyD_C3/zFC3)*100;
DyD_C5 = Result2[3,2];# [kmol]
zFC5 = zF[4];# [kmol]
percentC5 = (DyD_C5/zFC5)*100;
# These do not quite meet the original specification.
# For Tray 2 & C4
# From Eqn. 9.195:
# Value4 = G2*y2/(D*zD)
n = 2;
Sum = 0;
for i in range(0,n):
Val = 1;
for j in range(i,n):
Val = Val*A[j];
Sum = Sum+Val;
Value4 = 1+Sum;
# From The enthalpy Balnce:
G2 = 0.675;
# From Eqn. 9.211:
y2 = Value4*DyD/G2;
# Similarly:
# Value4 = [C1 C2 C3 C4 C5 C6]
Value4 = numpy.array([1.0235, 1.1062, 1.351, 2.705, 10.18 ,46.9]);
y2= numpy.array([0.04548,0.114,0.2896,0.53498,0.09697,0.001153]);
Y2 = sum(y2);
# Since Y2 is not equal to 1. THerefore the original temperature is incorrect. By adjusting y2 to unity.
# The dew point is 77 OC instead of 66 OC
# y2_adjusted = [C1 C2 C3 C4 C5 C6]
y2_adjusted = numpy.array([0.0419 ,0.1059 ,0.2675 ,0.4939, 0.0896, 0.00106]);
print"*****************Composition By Thiele Geddes Method*****************\n"
print"Component\t \t \t y2\n"
for i in range(0,6):
print"C",i,"\t \t \t \t",y2_adjusted[i]
# some values of solution in the textbook are incorrect