Rectifying and Stripping Cascade Trajectories 
for Isopropyl Acetate Chemistry 
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Reference :

S. B. Gadewar, N. Chadda, M. F. Malone and M. F. Doherty, Feasibility and Process Alternatives for Reactive Distillation, Reactive Distillation, Status and Future Directions, Edited by Kai Sundmacher & Achim Kienle, Wiley VCH 2002.

Activity Coefficients based on NRTL model
Component 1: Acetic Acid
Component 2: Isopropanol
Component 3: Isopropyl Acetate
Component 4: Water

Gam[i_]:=Exp[(((∑_(j=1)^4)( tau[j,i]G[j,i] x[j]))/((∑_(l=1)^4) (G[l,i] x[l])))+(∑_(j=1)^4) (((G[i,j] x[j])/((∑_(l=1)^4)( G[l,j] x[l])))(tau[i,j]-(((∑_(r=1)^4)( tau[r,j]G[r,j] x[r]))/((∑_(l=1)^4) (G[l,j] x[l])))))]/.x[1]->x1/.x [2]->x2/.x[3]->x3/.x[4]->x4

G[i_,j_]:=Exp[-a[i,j] tau[i,j]]

tau[i_,j_]:=g[i,j]/(R  T)

g[1,1]=0;g[2,1]=81.3926;g[3,1]=154.7885;g[4,1]=842.6081;
g[1,2]=-281.4482;g[2,2]=0;g[3,2]=140.0972;g[4,2]=1655.2550;
g[1,3]=141.0082;g[2,3]=269.9609;g[3,3]=0;g[4,3]=1270.2036;
g[1,4]=-219.7238;g[2,4]=39.8541;g[3,4]=1165.709;g[4,4]=0;

a[1,1]=0;a[2,1]=0.3048;a[3,1]=0.3014;a[4,1]=0.2997;
a[1,2]=0.3048;a[2,2]=0;a[3,2]=0.3009;a[4,2]=0.3255;
a[1,3]=0.3014;a[2,3]=0.3009;a[3,3]=0;a[4,3]=0.3300;
a[1,4]=0.2997;a[2,4]=0.3255;a[3,4]=0.3300;a[4,4]=0;

Gas Constant and Total Pressure

R=1.987;P=1.013 10^5;

Vapor Pressure based on Antoine Equation and Normal Boiling Points

Psat[i_]:=Exp[(Aa[i]+Ba[i]/(T+Ca[i]))]

Aa[1]=23.3618;Aa[2]=25.3358;
Aa[3]=21.7798;Aa[4]=23.4776;
Ba[1]=-4457.83;Ba[2]=-4628.96;Ba[3]=-3307.73;Ba[4]=-3984.92;
Ca[1]=-14.699;Ca[2]=-20.514;Ca[3]=-39.485;Ca[4]=-39.724;

FindRoot[Psat[1]== 1.013 10^5,{T,300}]

FindRoot[Psat[2]== 1.013 10^5,{T,300}]

FindRoot[Psat[3]== 1.013 10^5,{T,300}]

FindRoot[Psat[4]== 1.013 10^5,{T,300}]

Equilibrium Constant for Dimerization of Acetic Acid

K=10^(-12.5454+3166/(T));

Equilibrium constant for esterification

Keq=8.7;

Rectifying Cascade Trajectory at Damköhler number Da=3

t=0;Y[1,0]=0.5;Y[2,0]=0.5;y[1,-1]=0.5;y[2,-1]=0.5;y[3,-1]=0.0;y[4,-1]=0.0;x[1,-1]=0.5;x[2,-1]=0.5;x[3,-1]=0.0;x[4,-1]=0.0;

ϕ=0.5; Dd=0.75;

While[t<60,
{
sol=Chop[FindRoot[{
x1+x2+x3+x4== 1,y[1,t-1]== ϕ y1+(1-ϕ) x1+Dd/(1-Dd) (Gam[1] Gam[2]x1 x2-Gam[3] Gam[4]x3 x4/Keq),y[2,t-1]== ϕ y2+(1-ϕ) x2+Dd/(1-Dd) (Gam[1] Gam[2]x1 x2-Gam[3] Gam[4]x3 x4/Keq),
y[3,t-1]== ϕ y3+(1-ϕ) x3-Dd/(1-Dd) (Gam[1] Gam[2]x1 x2-Gam[3] Gam[4]x3 x4/Keq),
Psat[1] Gam[1]x1/((1+Sqrt[(1+4 K Psat[1])])/(1+Sqrt[(1+4 K P y1 (2-y1))]))+Psat[2] Gam[2]x2/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))])))+Psat[3] Gam[3]x3/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))])))+Psat[4] Gam[4]x4/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))])))== P,
y1== Psat[1] Gam[1]x1/((1+Sqrt[(1+4 K Psat[1])])/(1+Sqrt[(1+4 K P y1 (2-y1))])) 1/P,
y2== Psat[2] Gam[2]x2/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))]))) 1/P,
y3== Psat[3] Gam[3]x3/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))])))1/P},{T,350},{x3,x[3,t-1]},{x4,x[4,t-1]},{x1,x[1,t-1]},{x2,x[2,t-1]},{y1,y[1,t-1]},{y2,y[2,t-1]},{y3,y[3,t-1]},MaxIterations->100000]],
Y[1,t+1]=y1+y3/.sol,
Y[2,t+1]=y2+y3/.sol,
x[3,t]=x3/.sol,x[2,t]=x2/.sol,x[4,t]=x4/.sol,x[1,t]=x1/.sol,y[3,t]=y3/.sol,y[2,t]=y2/.sol,y[1,t]=y1/.sol,
t++,Print[t]}]//Timing

ListPlot[Table[{Y[1,t],Y[2,t]},{t,-1,60}],PlotRange-> {{0,1},{0,1}},AspectRatio-> 1]

Rectifying Cascade Trajectory at Damköhler number Da=0.33

t=0;Y[1,0]=0.5;Y[2,0]=0.5;y[1,-1]=0.5;y[2,-1]=0.5;y[3,-1]=0.0;y[4,-1]=0.0;x[1,-1]=0.5;x[2,-1]=0.5;x[3,-1]=0.0;x[4,-1]=0.0;

ϕ=0.5; Dd=0.25;

While[t<60,
{
sol=Chop[FindRoot[{
x1+x2+x3+x4== 1,y[1,t-1]== ϕ y1+(1-ϕ) x1+Dd/(1-Dd) (Gam[1] Gam[2]x1 x2-Gam[3] Gam[4]x3 x4/Keq),y[2,t-1]== ϕ y2+(1-ϕ) x2+Dd/(1-Dd) (Gam[1] Gam[2]x1 x2-Gam[3] Gam[4]x3 x4/Keq),
y[3,t-1]== ϕ y3+(1-ϕ) x3-Dd/(1-Dd) (Gam[1] Gam[2]x1 x2-Gam[3] Gam[4]x3 x4/Keq),
Psat[1] Gam[1]x1/((1+Sqrt[(1+4 K Psat[1])])/(1+Sqrt[(1+4 K P y1 (2-y1))]))+Psat[2] Gam[2]x2/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))])))+Psat[3] Gam[3]x3/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))])))+Psat[4] Gam[4]x4/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))])))== P,
y1== Psat[1] Gam[1]x1/((1+Sqrt[(1+4 K Psat[1])])/(1+Sqrt[(1+4 K P y1 (2-y1))])) 1/P,
y2== Psat[2] Gam[2]x2/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))]))) 1/P,
y3== Psat[3] Gam[3]x3/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))])))1/P},{T,350},{x3,x[3,t-1]},{x4,x[4,t-1]},{x1,x[1,t-1]},{x2,x[2,t-1]},{y1,y[1,t-1]},{y2,y[2,t-1]},{y3,y[3,t-1]},MaxIterations->100000]],
Y[1,t+1]=y1+y3/.sol,
Y[2,t+1]=y2+y3/.sol,
x[3,t]=x3/.sol,x[2,t]=x2/.sol,x[4,t]=x4/.sol,x[1,t]=x1/.sol,y[3,t]=y3/.sol,y[2,t]=y2/.sol,y[1,t]=y1/.sol,
t++,Print[t]}]//Timing

ListPlot[Table[{Y[1,t],Y[2,t]},{t,-1,60}],PlotRange-> {{0,1},{0,1}},AspectRatio-> 1]

Stripping Cascade Trajectory at Damköhler number Da=3

t=0;X[1,0]=0.5;X[2,0]=0.5;y[1,-1]=0.5;y[2,-1]=0.5;y[3,-1]=0.0;y[4,-1]=0.0;x[1,-1]=0.5;x[2,-1]=0.5;x[3,-1]=0.0;x[4,-1]=0.0;

ϕ=0.5; Dd=0.75;

While[t<60,
{
sol=Chop[FindRoot[{
x1+x2+x3+x4== 1,x[1,t-1]== ϕ y1+(1-ϕ) x1+Dd/(1-Dd) (Gam[1] Gam[2]x1 x2-Gam[3] Gam[4]x3 x4/Keq),x[2,t-1]== ϕ y2+(1-ϕ) x2+Dd/(1-Dd) (Gam[1] Gam[2]x1 x2-Gam[3] Gam[4]x3 x4/Keq),
x[3,t-1]== ϕ y3+(1-ϕ) x3-Dd/(1-Dd) (Gam[1] Gam[2]x1 x2-Gam[3] Gam[4]x3 x4/Keq),
Psat[1] Gam[1]x1/((1+Sqrt[(1+4 K Psat[1])])/(1+Sqrt[(1+4 K P y1 (2-y1))]))+Psat[2] Gam[2]x2/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))])))+Psat[3] Gam[3]x3/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))])))+Psat[4] Gam[4]x4/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))])))== P,
y1== Psat[1] Gam[1]x1/((1+Sqrt[(1+4 K Psat[1])])/(1+Sqrt[(1+4 K P y1 (2-y1))])) 1/P,
y2== Psat[2] Gam[2]x2/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))]))) 1/P,
y3== Psat[3] Gam[3]x3/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))])))1/P},{T,350},{x3,x[3,t-1]},{x4,x[4,t-1]},{x1,x[1,t-1]},{x2,x[2,t-1]},{y1,y[1,t-1]},{y2,y[2,t-1]},{y3,y[3,t-1]},MaxIterations->100000]],
X[1,t+1]=x1+x3/.sol,
X[2,t+1]=x2+x3/.sol,
x[3,t]=x3/.sol,x[2,t]=x2/.sol,x[4,t]=x4/.sol,x[1,t]=x1/.sol,y[3,t]=y3/.sol,y[2,t]=y2/.sol,y[1,t]=y1/.sol,
t++,Print[t]}]//Timing

ListPlot[Table[{X[1,t],X[2,t]},{t,-1,60}],PlotRange-> {{0,1},{0,1}},AspectRatio-> 1]

Stripping Cascade Trajectory at Damköhler number Da=0.33

t=0;X[1,0]=0.5;X[2,0]=0.5;y[1,-1]=0.5;y[2,-1]=0.5;y[3,-1]=0.0;y[4,-1]=0.0;x[1,-1]=0.5;x[2,-1]=0.5;x[3,-1]=0.0;x[4,-1]=0.0;

ϕ=0.5; Dd=0.25;

While[t<60,
{
sol=Chop[FindRoot[{
x1+x2+x3+x4== 1,x[1,t-1]== ϕ y1+(1-ϕ) x1+Dd/(1-Dd) (Gam[1] Gam[2]x1 x2-Gam[3] Gam[4]x3 x4/Keq),x[2,t-1]== ϕ y2+(1-ϕ) x2+Dd/(1-Dd) (Gam[1] Gam[2]x1 x2-Gam[3] Gam[4]x3 x4/Keq),
x[3,t-1]== ϕ y3+(1-ϕ) x3-Dd/(1-Dd) (Gam[1] Gam[2]x1 x2-Gam[3] Gam[4]x3 x4/Keq),
Psat[1] Gam[1]x1/((1+Sqrt[(1+4 K Psat[1])])/(1+Sqrt[(1+4 K P y1 (2-y1))]))+Psat[2] Gam[2]x2/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))])))+Psat[3] Gam[3]x3/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))])))+Psat[4] Gam[4]x4/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))])))== P,
y1== Psat[1] Gam[1]x1/((1+Sqrt[(1+4 K Psat[1])])/(1+Sqrt[(1+4 K P y1 (2-y1))])) 1/P,
y2== Psat[2] Gam[2]x2/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))]))) 1/P,
y3== Psat[3] Gam[3]x3/((2 (1-y1+Sqrt[(1+4 K P y1 (2-y1))]))/((2-y1) (1+Sqrt[(1+4 K P y1 (2-y1))])))1/P},{T,350},{x3,x[3,t-1]},{x4,x[4,t-1]},{x1,x[1,t-1]},{x2,x[2,t-1]},{y1,y[1,t-1]},{y2,y[2,t-1]},{y3,y[3,t-1]},MaxIterations->100000]],
X[1,t+1]=x1+x3/.sol,
X[2,t+1]=x2+x3/.sol,
x[3,t]=x3/.sol,x[2,t]=x2/.sol,x[4,t]=x4/.sol,x[1,t]=x1/.sol,y[3,t]=y3/.sol,y[2,t]=y2/.sol,y[1,t]=y1/.sol,
t++,Print[t]}]//Timing

ListPlot[Table[{X[1,t],X[2,t]},{t,-1,60}],PlotRange-> {{0,1},{0,1}},AspectRatio-> 1]