Distillation Curves for Ternary Systems
​Case of Acetone - Chloroform - Benzene System
​Authors: Housam Binous and Ahmed Bellagi

Reference: L. T. Biegler, I. E. Grossmann and A. W. Westerberg, Systematic Methods of Chemical Process Design, Prentice Hall, 1999.

Acetone = component 1
Chloroform = component 2
Benzene = component 3

In[]:=
Off[General::"spell1"]
Gas Constant and Total Pressure
In[]:=
R=1.987;P=760;
Vapor Pressure Using Antoine Equation
In[]:=
A1=7.11714;B1=1210.595;C1=229.664;​​A2=6.95465;B2=1170.966;C2=226.232;​​A3=6.87987;B3=1196.760;C3=219.161;​​
In[]:=
PS1[i_]:=10^(A1-B1/(C1+T));​​PS2[i_]:=10^(A2-B2/(C2+T));​​PS3[i_]:=10^(A3-B3/(C3+T));
Binary interaction Parameters and Activities for Wilson Model
In[]:=
l12=116.1171;l21=-506.8519;l13=682.40610;​​l31=-243.96510;l23=-71.81089;l32=-11.82310;​​V1=74.05;V2=80.67;V3=89.41;​​
In[]:=
A12=V2/V1Exp[-l12/(R*(273.15+T))];​​A13=V3/V1Exp[-l13/(R*(273.15+T))];​​A32=V2/V3Exp[-l32/(R*(273.15+T))];​​A21=V1/V2Exp[-l21/(R*(273.15+T))];​​A31=V1/V3Exp[-l31/(R*(273.15+T))];​​A23=V3/V2Exp[-l23/(R*(273.15+T))];
In[]:=
GAM1[i_]:=Exp[-Log[x1[i]+x2[i]*A12+x3[i]*A13]+1-(x1[i]/(x1[i]+x2[i]*A12+x3[i]*A13)+x2[i]*A21/(x1[i]*A21+x2[i]+x3[i]*A23)+x3[i]*A31/(x1[i]*A31+x2[i]*A32+x3[i]))];​​GAM2[i_]:=Exp[-Log[x1[i]*A21+x2[i]+x3[i]*A23]+1-(x1[i]*A12/(x1[i]+x2[i]*A12+x3[i]*A13)+x2[i]/(x1[i]*A21+x2[i]+x3[i]*A23)+x3[i]*A32/(x1[i]*A31+x2[i]*A32+x3[i]))];​​GAM3[i_]:=Exp[-Log[x1[i]*A31+x2[i]*A32+x3[i]]+1-(x1[i]*A13/(x1[i]+x2[i]*A12+x3[i]*A13)+x2[i]*A23/(x1[i]*A21+x2[i]+x3[i]*A23)+x3[i]/(x1[i]*A31+x2[i]*A32+x3[i]))];
In[]:=
GAMMA1[i_]:=GAM1[i]/.x1[i]x1/.x2[i]x2/.x3[i]x3;
In[]:=
GAMMA2[i_]:=GAM2[i]/.x1[i]x1/.x2[i]x2/.x3[i]x3;
In[]:=
GAMMA3[i_]:=GAM3[i]/.x1[i]x1/.x2[i]x2/.x3[i]x3
Calculation of Distillation Curves for different initial conditions. We get the same figure as figure 14.8 page 466 of the above reference.
Green point is the azeotrope.
In[]:=
x1[0]=0.2;x2[0]=0.2;x3[0]=0.6;i=0;
In[]:=
While[i<200,{T=Quiet@FindRoot[PEvaluate[x1[i]PS1[i]GAM1[i]+x2[i]PS2[i]GAM2[i]+x3[i]PS3[i]GAM3[i]],{T,100}][[1,2]],y1[i]=x1[i]PS1[i]GAM1[i]/P,​​y2[i]=x2[i]PS2[i]GAM2[i]/P,​​y3[i]=x3[i]PS3[i]GAM3[i]/P,x1[i+1]=y1[i],x2[i+1]=y2[i],x3[i+1]=y3[i],i++}];
In[]:=
t1=Table[{x2[i],x3[i],1-x2[i]-x3[i]},{i,0,200}];
In[]:=
T=.
In[]:=
i=0;x1[0]=0.2;x2[0]=0.2;x3[0]=0.6;y1[0]=0.2;y2[0]=0.2;y3[0]=0.6;
In[]:=
While[i<200,{sol=Quiet@FindRoot[{x1+x2+x31,y1[i]x1PS1[i]GAMMA1[i]/P,​​y2[i]x2PS2[i]GAMMA2[i]/P,​​y3[i]x3PS3[i]GAMMA3[i]/P},{T,60,90},{x1,0.01,0.2},{x2,0.01,0.3},{x3,0.5,0.9},MaxIterations1000],x1[i+1]=Chop[sol[[2,2]]],x2[i+1]=Chop[sol[[3,2]]],x3[i+1]=Chop[sol[[4,2]]],y1[i+1]=x1[i],y2[i+1]=x2[i],y3[i+1]=x3[i],i++}];
In[]:=
t2=Table[{x2[i],x3[i],1-x2[i]-x3[i]},{i,0,200}];
In[]:=
plt2=TernaryListPlot[{t1,t2},PlotStyle->{{Darker@Blue,Thick},{Darker@Blue,Thick}},GridLines->True,PlotRangeAll,PlotMarkers->{Darker@Blue,Darker@Blue},Joined->True,Method->{"Backdrop"LightBlue}];
In[]:=
x1[0]=0.8;x2[0]=0.1;x3[0]=0.1;i=0;
In[]:=
While[i<200,{T=Quiet@FindRoot[PEvaluate[x1[i]PS1[i]GAM1[i]+x2[i]PS2[i]GAM2[i]+x3[i]PS3[i]GAM3[i]],{T,100}][[1,2]],y1[i]=x1[i]PS1[i]GAM1[i]/P,​​y2[i]=x2[i]PS2[i]GAM2[i]/P,​​y3[i]=x3[i]PS3[i]GAM3[i]/P,x1[i+1]=y1[i],x2[i+1]=y2[i],x3[i+1]=y3[i],i++}];
In[]:=
t3=Table[{x2[i],x3[i],1-x2[i]-x3[i]},{i,0,200}];
In[]:=
T=.
In[]:=
i=0;x1[0]=0.8;x2[0]=0.1;x3[0]=0.1;y1[0]=0.8;y2[0]=0.1;y3[0]=0.1;
In[]:=
While[i<200,{sol=Quiet@FindRoot[{x1+x2+x31,y1[i]x1PS1[i]GAMMA1[i]/P,​​y2[i]x2PS2[i]GAMMA2[i]/P,​​y3[i]x3PS3[i]GAMMA3[i]/P},{T,60,90},{x1,0.01,0.2},{x2,0.01,0.3},{x3,0.5,0.9},MaxIterations1000],x1[i+1]=Chop[sol[[2,2]]],x2[i+1]=Chop[sol[[3,2]]],x3[i+1]=Chop[sol[[4,2]]],y1[i+1]=x1[i],y2[i+1]=x2[i],y3[i+1]=x3[i],i++}];
In[]:=
t4=Table[{x2[i],x3[i],1-x2[i]-x3[i]},{i,0,200}];
In[]:=
plt4=TernaryListPlot[{t3,t4},PlotStyle->{{Darker@Blue,Thick},{Darker@Blue,Thick}},PlotMarkers->{Darker@Blue,Darker@Blue},GridLines->True,PlotRangeAll,Joined->True];
In[]:=
x1[0]=0.1;x2[0]=0.8;x3[0]=0.1;i=0;
In[]:=
While[i<200,{T=Quiet@FindRoot[PEvaluate[x1[i]PS1[i]GAM1[i]+x2[i]PS2[i]GAM2[i]+x3[i]PS3[i]GAM3[i]],{T,100}][[1,2]],y1[i]=x1[i]PS1[i]GAM1[i]/P,​​y2[i]=x2[i]PS2[i]GAM2[i]/P,​​y3[i]=x3[i]PS3[i]GAM3[i]/P,x1[i+1]=y1[i],x2[i+1]=y2[i],x3[i+1]=y3[i],i++}];
In[]:=
t5=Table[{x2[i],x3[i],1-x2[i]-x3[i]},{i,0,200}];
In[]:=
T=.
In[]:=
i=0;x1[0]=0.1;x2[0]=0.8;x3[0]=0.1;y1[0]=0.1;y2[0]=0.8;y3[0]=0.1;
In[]:=
While[i<200,{sol=Quiet@FindRoot[{x1+x2+x31,y1[i]x1PS1[i]GAMMA1[i]/P,​​y2[i]x2PS2[i]GAMMA2[i]/P,​​y3[i]x3PS3[i]GAMMA3[i]/P},{T,60,90},{x1,0.01,0.2},{x2,0.01,0.3},{x3,0.5,0.9},MaxIterations1000],x1[i+1]=Chop[sol[[2,2]]],x2[i+1]=Chop[sol[[3,2]]],x3[i+1]=Chop[sol[[4,2]]],y1[i+1]=x1[i],y2[i+1]=x2[i],y3[i+1]=x3[i],i++}];
In[]:=
t6=Table[{x2[i],x3[i],1-x2[i]-x3[i]},{i,0,200}];
In[]:=
plt6=TernaryListPlot[{t5,t6},PlotStyle->{{Darker@Blue,Thick},{Darker@Blue,Thick}},GridLines->True,PlotRangeAll,Joined->True,PlotMarkers->{Darker@Blue,Darker@Blue}];
In[]:=
x1[0]=0.33;x2[0]=0.65;x3[0]=0.02;i=0;
In[]:=
While[i<200,{T=Quiet@FindRoot[PEvaluate[x1[i]PS1[i]GAM1[i]+x2[i]PS2[i]GAM2[i]+x3[i]PS3[i]GAM3[i]],{T,100}][[1,2]],y1[i]=x1[i]PS1[i]GAM1[i]/P,​​y2[i]=x2[i]PS2[i]GAM2[i]/P,​​y3[i]=x3[i]PS3[i]GAM3[i]/P,x1[i+1]=y1[i],x2[i+1]=y2[i],x3[i+1]=y3[i],i++}];
In[]:=
t7=Table[{x2[i],x3[i],1-x2[i]-x3[i]},{i,0,200}];
In[]:=
T=.
In[]:=
i=0;x1[0]=0.33;x2[0]=0.65;x3[0]=0.02;y1[0]=0.33;y2[0]=0.65;y3[0]=0.02;
In[]:=
While[i<200,{sol=Quiet@FindRoot[{x1+x2+x31,y1[i]x1PS1[i]GAMMA1[i]/P,​​y2[i]x2PS2[i]GAMMA2[i]/P,​​y3[i]x3PS3[i]GAMMA3[i]/P},{T,60,90},{x1,0.01,0.2},{x2,0.01,0.3},{x3,0.5,0.9},MaxIterations1000],x1[i+1]=Chop[sol[[2,2]]],x2[i+1]=Chop[sol[[3,2]]],x3[i+1]=Chop[sol[[4,2]]],y1[i+1]=x1[i],y2[i+1]=x2[i],y3[i+1]=x3[i],i++}];
In[]:=
t8=Table[{x2[i],x3[i],1-x2[i]-x3[i]},{i,0,200}];
In[]:=
plt8=TernaryListPlot[{t7,t8},PlotMarkers->{Darker@Blue,Darker@Blue},PlotStyle->{{Darker@Blue,Thick},{Darker@Blue,Thick}},GridLines->True,PlotRangeAll,Joined->True];
In[]:=
x1[0]=0.31;x2[0]=0.67;x3[0]=0.02;i=0;
In[]:=
While[i<200,{T=Quiet@FindRoot[PEvaluate[x1[i]PS1[i]GAM1[i]+x2[i]PS2[i]GAM2[i]+x3[i]PS3[i]GAM3[i]],{T,100}][[1,2]],y1[i]=x1[i]PS1[i]GAM1[i]/P,​​y2[i]=x2[i]PS2[i]GAM2[i]/P,​​y3[i]=x3[i]PS3[i]GAM3[i]/P,x1[i+1]=y1[i],x2[i+1]=y2[i],x3[i+1]=y3[i],i++}];
In[]:=
t9=Table[{x2[i],x3[i],1-x2[i]-x3[i]},{i,0,200}];
In[]:=
T=.
In[]:=
i=0;x1[0]=0.31;x2[0]=0.67;x3[0]=0.02;y1[0]=0.31;y2[0]=0.67;y3[0]=0.02;
In[]:=
While[i<200,{sol=Quiet@FindRoot[{x1+x2+x31,y1[i]x1PS1[i]GAMMA1[i]/P,​​y2[i]x2PS2[i]GAMMA2[i]/P,​​y3[i]x3PS3[i]GAMMA3[i]/P},{T,60,90},{x1,0.01,0.2},{x2,0.01,0.3},{x3,0.5,0.9},MaxIterations1000],x1[i+1]=Chop[sol[[2,2]]],x2[i+1]=Chop[sol[[3,2]]],x3[i+1]=Chop[sol[[4,2]]],y1[i+1]=x1[i],y2[i+1]=x2[i],y3[i+1]=x3[i],i++}];