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

Reference :

Doherty and Malone, Conceptual Design of Distillation Systems, McGraw-Hill, 2001.

Relative volatilities

(Subscript[α, AC])=5;(Subscript[α, BC])=3;(Subscript[α, CC])=1;

Molecular Weigths

M[1]=86.18;M[2]=100.21;M[3]=186.39;

Equilibrium constant for reaction

Keq=2;

Rectifying Cascade Trajectory at Damköhler number Da=0.277

t=0;y[1,-1]=0.4;y[2,-1]=0.4;y[3,-1]=0.2;x[1,-1]=0.4;x[2,-1]=0.4;x[3,-1]=0.2;

ϕ=0.5;Da=0.277;

While[t<60,
{
sol=Chop[FindRoot[{
y[1,t-1]== ϕ (((∑_(i=1)^3)y[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3])) y1+(1-ϕ  (((∑_(i=1)^3)y[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3]))) x1-(((∑_(i=1)^3)y[i,t-1] M[i])/(x1 M[1]+x2 M[2]+x3 M[3]))Da  (x1 x2-x3/Keq) (-1+x1),y[2,t-1]== ϕ (((∑_(i=1)^3)y[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3])) y2+(1-ϕ(((∑_(i=1)^3)y[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3]))) x2-(((∑_(i=1)^3)y[i,t-1] M[i])/(x1 M[1]+x2 M[2]+x3 M[3]))Da (x1 x2-x3/Keq) (-1+x2),
y[3,t-1]== ϕ (((∑_(i=1)^3)y[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3]))y3+(1-ϕ(((∑_(i=1)^3)y[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3]))) x3-(((∑_(i=1)^3)y[i,t-1] M[i])/(x1 M[1]+x2 M[2]+x3 M[3]))Da (x1 x2-x3/Keq) (1+x3),
y1==x1/(x1+(Subscript[α, BC])/(Subscript[α, AC]) x2+1/(Subscript[α, AC]) x3),
y2==x2/((Subscript[α, AC])/(Subscript[α, BC])x1+x2+1/(Subscript[α, BC]) x3),
y3==x3/((Subscript[α, AC]) x1+(Subscript[α, BC]) x2+x3)},{x3,x[3,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[3,t]=x3/.sol,x[2,t]=x2/.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

plt1=ListPlot[Table[{y[1,t],y[2,t]},{t,-1,59}],PlotRange-> {{0,1},{0,1}},AspectRatio-> 1,Epilog-> Line[{{1,0},{0,1}}]]

Stripping Cascade Trajectory at Damköhler number Da=0.277

t=0;y[1,-1]=0.4;y[2,-1]=0.4;y[3,-1]=0.2;x[1,-1]=0.4;x[2,-1]=0.4;x[3,-1]=0.2;

ϕ=0.5;Da=0.277;

While[t<60,
{
sol=Chop[FindRoot[{
x[1,t-1]== ϕ (((∑_(i=1)^3)x[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3]))y1+(1-ϕ(((∑_(i=1)^3)x[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3]))) x1-Da (((∑_(i=1)^3)x[i,t-1] M[i])/(x1 M[1]+x2 M[2]+x3 M[3]))(x1 x2-x3/Keq) (-1+x1),x[2,t-1]== ϕ (((∑_(i=1)^3)x[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3]))y2+(1-ϕ(((∑_(i=1)^3)x[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3]))) x2-Da (((∑_(i=1)^3)x[i,t-1] M[i])/(x1 M[1]+x2 M[2]+x3 M[3]))(x1 x2-x3/Keq) (-1+x2),
x1+x2+x3==1,
y1==x1/(x1+(Subscript[α, BC])/(Subscript[α, AC]) x2+1/(Subscript[α, AC]) x3),
y2==x2/((Subscript[α, AC])/(Subscript[α, BC])x1+x2+1/(Subscript[α, BC]) x3),
y3==x3/((Subscript[α, AC]) x1+(Subscript[α, BC]) x2+x3)},{x3,x[3,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->1000000]],
x[3,t]=x3/.sol,x[2,t]=x2/.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

plt2=ListPlot[Table[{x[1,t],x[2,t]},{t,-1,59}],PlotRange-> {{0,1},{0,1}},AspectRatio-> 1,PlotStyle-> Red]

{x[1,59],x[2,59]}

FindRoot[{
x1== ϕ ((x1 M[1]+x2 M[2]+x3 M[3])/(y1 M[1]+y2 M[2]+y3 M[3]))y1+(1-ϕ((x1 M[1]+x2 M[2]+x3 M[3])/(y1 M[1]+y2 M[2]+y3 M[3]))) x1-Da ((x1 M[1]+x2 M[2]+x3 M[3])/(x1 M[1]+x2 M[2]+x3 M[3]))(x1 x2-x3/Keq) (-1+x1),x2== ϕ ((x1 M[1]+x2 M[2]+x3 M[3])/(y1 M[1]+y2 M[2]+y3 M[3]))y2+(1-ϕ((x1 M[1]+x2 M[2]+x3 M[3])/(y1 M[1]+y2 M[2]+y3 M[3]))) x2-Da ((x1 M[1]+x2 M[2]+x3 M[3])/(x1 M[1]+x2 M[2]+x3 M[3]))(x1 x2-x3/Keq) (-1+x2),
x1+x2+x3==1,
y1==x1/(x1+(Subscript[α, BC])/(Subscript[α, AC]) x2+1/(Subscript[α, AC]) x3),
y2==x2/((Subscript[α, AC])/(Subscript[α, BC])x1+x2+1/(Subscript[α, BC]) x3),
y3==x3/((Subscript[α, AC]) x1+(Subscript[α, BC]) x2+x3)},{x3,0.7},{x1,0.1},{x2,0.2},{y1,0.1},{y2,0.3},{y3,0.6},MaxIterations->1000000]

Stripping & Rectifying Profiles at Damköhler number Da=0.277

Show[plt1,plt2]

Rectifying Cascade Trajectory at Damköhler number Da=11.1

t=0;y[1,-1]=0.4;y[2,-1]=0.4;y[3,-1]=0.2;x[1,-1]=0.4;x[2,-1]=0.4;x[3,-1]=0.2;

ϕ=0.5;Da=11.1;

While[t<60,
{
sol=Chop[FindRoot[{
y[1,t-1]== ϕ (((∑_(i=1)^3)y[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3])) y1+(1-ϕ  (((∑_(i=1)^3)y[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3]))) x1-(((∑_(i=1)^3)y[i,t-1] M[i])/(x1 M[1]+x2 M[2]+x3 M[3]))Da  (x1 x2-x3/Keq) (-1+x1),y[2,t-1]== ϕ (((∑_(i=1)^3)y[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3])) y2+(1-ϕ(((∑_(i=1)^3)y[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3]))) x2-(((∑_(i=1)^3)y[i,t-1] M[i])/(x1 M[1]+x2 M[2]+x3 M[3]))Da (x1 x2-x3/Keq) (-1+x2),
y[3,t-1]== ϕ (((∑_(i=1)^3)y[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3]))y3+(1-ϕ(((∑_(i=1)^3)y[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3]))) x3-(((∑_(i=1)^3)y[i,t-1] M[i])/(x1 M[1]+x2 M[2]+x3 M[3]))Da (x1 x2-x3/Keq) (1+x3),
y1==x1/(x1+(Subscript[α, BC])/(Subscript[α, AC]) x2+1/(Subscript[α, AC]) x3),
y2==x2/((Subscript[α, AC])/(Subscript[α, BC])x1+x2+1/(Subscript[α, BC]) x3),
y3==x3/((Subscript[α, AC]) x1+(Subscript[α, BC]) x2+x3)},{x3,x[3,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[3,t]=x3/.sol,x[2,t]=x2/.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

plt1=ListPlot[Table[{y[1,t],y[2,t]},{t,-1,59}],PlotRange-> {{0,1},{0,1}},AspectRatio-> 1,Epilog-> Line[{{1,0},{0,1}}]]

Stripping Cascade Trajectory at Damköhler number Da=11.1

t=0;y[1,-1]=0.4;y[2,-1]=0.4;y[3,-1]=0.2;x[1,-1]=0.4;x[2,-1]=0.4;x[3,-1]=0.2;

ϕ=0.5;Da=11.1;

While[t<200,
{
sol=Chop[FindRoot[{
x[1,t-1]== ϕ (((∑_(i=1)^3)x[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3]))y1+(1-ϕ(((∑_(i=1)^3)x[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3]))) x1-Da (((∑_(i=1)^3)x[i,t-1] M[i])/(x1 M[1]+x2 M[2]+x3 M[3]))(x1 x2-x3/Keq) (-1+x1),x[2,t-1]== ϕ (((∑_(i=1)^3)x[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3]))y2+(1-ϕ(((∑_(i=1)^3)x[i,t-1] M[i])/(y1 M[1]+y2 M[2]+y3 M[3]))) x2-Da (((∑_(i=1)^3)x[i,t-1] M[i])/(x1 M[1]+x2 M[2]+x3 M[3]))(x1 x2-x3/Keq) (-1+x2),
x1+x2+x3==1,
y1==x1/(x1+(Subscript[α, BC])/(Subscript[α, AC]) x2+1/(Subscript[α, AC]) x3),
y2==x2/((Subscript[α, AC])/(Subscript[α, BC])x1+x2+1/(Subscript[α, BC]) x3),
y3==x3/((Subscript[α, AC]) x1+(Subscript[α, BC]) x2+x3)},{x3,x[3,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->1000000]],
x[3,t]=x3/.sol,x[2,t]=x2/.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

plt2=ListPlot[Table[{x[1,t],x[2,t]},{t,-1,199}],PlotRange-> {{0,1},{0,1}},AspectRatio-> 1,PlotStyle->Red]

{x[1,199],x[2,199]}

FindRoot[{
x1== ϕ ((x1 M[1]+x2 M[2]+x3 M[3])/(y1 M[1]+y2 M[2]+y3 M[3]))y1+(1-ϕ((x1 M[1]+x2 M[2]+x3 M[3])/(y1 M[1]+y2 M[2]+y3 M[3]))) x1-Da ((x1 M[1]+x2 M[2]+x3 M[3])/(x1 M[1]+x2 M[2]+x3 M[3]))(x1 x2-x3/Keq) (-1+x1),x2== ϕ ((x1 M[1]+x2 M[2]+x3 M[3])/(y1 M[1]+y2 M[2]+y3 M[3]))y2+(1-ϕ((x1 M[1]+x2 M[2]+x3 M[3])/(y1 M[1]+y2 M[2]+y3 M[3]))) x2-Da ((x1 M[1]+x2 M[2]+x3 M[3])/(x1 M[1]+x2 M[2]+x3 M[3]))(x1 x2-x3/Keq) (-1+x2),
x1+x2+x3==1,
y1==x1/(x1+(Subscript[α, BC])/(Subscript[α, AC]) x2+1/(Subscript[α, AC]) x3),
y2==x2/((Subscript[α, AC])/(Subscript[α, BC])x1+x2+1/(Subscript[α, BC]) x3),
y3==x3/((Subscript[α, AC]) x1+(Subscript[α, BC]) x2+x3)},{x3,0.7},{x1,0.1},{x2,0.2},{y1,0.1},{y2,0.3},{y3,0.6},MaxIterations->1000000]

Chop[x1 x2-x3/Keq/.sol]

Stripping & Rectifying Profiles at Damköhler number Da=11.1

Show[plt1,plt2]