Bifurcation Analysis of Simple Distillation:
Case of a Hypothetical Ternary Mixture with Chemical Reaction
Author's Data: Housam Binous 
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Reference :

G. Venimadhavan, M. F. Malone and M. F. Doherty, Bifurcation Study of Kinetic Effects in Recative Distillation, AIChE Journal, March 1999, Vol. 45, No. 3, pp 546-556.

Vapor Pressure of components A, B and C

a1=23.1962;a2=22.6867;a3=21.5872;
b1=-3502.77;b2=-3502.77;b3=-3502.77;
c1=0;c2=0;c3=0;

Psat[1]=Exp[a1+b1/(c1+T)];
Psat[2]=Exp[a2+b2/(c2+T)];
Psat[3]=Exp[a3+b3/(c3+T)];

Total Pressure

P=1.01325 10^5;

y1=Psat[1] x1/P;
y2=Psat[2] x2/P;
y3=Psat[3] x3/P;

Summation rule

x3=1-x1-x2;

Equilibrium constant for the reaction A + B ⇄ C

Keq=5;

Computing the Jacobian Matrix

Grad[fon_,var_List]:=D[fon,#]&/@var

F[x1_,x2_,T_]:={(1-Dd) (x1-y1)+Dd (-1+x1) (x1 x2-x3/Keq),
(1-Dd) (x2-y2)+Dd (-1+x2) (x1 x2-x3/Keq),1-Psat[1] x1/P-Psat[2] x2/P-Psat[3] x3/P}

M=Grad[F[x1,x2,T],{x1,x2,T}]

Computing Branch 4

Dd=0.;i=0

sol[0]=FindRoot[{(1-Dd) (x1-y1)+Dd (-1+x1) (x1 x2-x3/Keq)==0,
(1-Dd) (x2-y2)+Dd (-1+x2) (x1 x2-x3/Keq)==0,Psat[1] x1/P+Psat[2] x2/P+Psat[3] x3/P== 1},{x1,0.0},{x2,0.1},{T,320.0}]//Chop

Eigenvalues[M]/.sol[0]

While[i< 100,i=i+1;{Dd=0.01*i,sol[i]=FindRoot[{(1-Dd) (x1-y1)+Dd (-1+x1) (x1 x2-x3/Keq)==0,
(1-Dd) (x2-y2)+Dd (-1+x2) (x1 x2-x3/Keq)==0,Psat[1] x1/P+Psat[2] x2/P+Psat[3] x3/P== 1},{x1,x1/.sol[i-1]},{x2,x2/.sol[i-1]},{T,T/.sol[i-1]},AccuracyGoal-> 5],j=0,For[k=1,k≤ 3,k++,If[Re[(Eigenvalues[M][[k]]/.sol[i])]<0,j=j+1]],test=If[j==3,"stable",If[j== 2,"saddle","unstable"]], Print[i," ",sol[i],Dd," ",Eigenvalues[M]/.sol[i]," ",test," ",j]}]

FindRoot[{(x1-y1) (-1+x2)== (x2-y2) (-1+x1),x1 x2-x3/Keq==0,y1+y2+y3==1},{x1,0.2},{x2,0.5},{T,330}]

plt1=ListPlot[Table[{0.01*i,(T/.sol[i])-273.15},{i,0,100,5}],PlotRange->{40,80},PlotMarkers->{"●"}]

plt2=ListPlot[Table[{x1/.sol[i],x2/.sol[i]},{i,0,100,5}],PlotRange->{{0,1},{0,1}},AspectRatio->1,PlotMarkers->{"●",14}]

Computing Branch 3

Dd=0.;i=0

sol[0]=FindRoot[{(1-Dd) (x1-y1)+Dd (-1+x1) (x1 x2-x3/Keq)==0,
(1-Dd) (x2-y2)+Dd (-1+x2) (x1 x2-x3/Keq)==0,Psat[1] x1/P+Psat[2] x2/P+Psat[3] x3/P== 1},{x1,0.0},{x2,0.1},{T,273.0}]//Chop

Eigenvalues[M]/.sol[0]

While[i< 100,i=i+1;{Dd=0.01*i,sol[i]=FindRoot[{(1-Dd) (x1-y1)+Dd (-1+x1) (x1 x2-x3/Keq)==0,
(1-Dd) (x2-y2)+Dd (-1+x2) (x1 x2-x3/Keq)==0,Psat[1] x1/P+Psat[2] x2/P+Psat[3] x3/P== 1},{x1,x1/.sol[i-1]},{x2,x2/.sol[i-1]},{T,T/.sol[i-1]},AccuracyGoal->5],j=0,For[k=1,k≤ 3,k++,If[Re[(Eigenvalues[M][[k]]/.sol[i])]<0,j=j+1]],test=If[j==3,"stable",If[j== 2,"saddle","unstable"]], Print[i," ",sol[i],Dd," ",Eigenvalues[M]/.sol[i]," ",test," ",j]}]

plt3=ListPlot[Table[{0.01*i,(T/.sol[i])-273.15},{i,0,100,5}],PlotRange->{30,50},PlotStyle->Red,PlotMarkers->{"□"}]

plt4=ListPlot[Table[{x1/.sol[i],x2/.sol[i]},{i,0,100,5}],PlotRange->{{-0.2,1},{0,1.2}},AspectRatio->1,PlotStyle->Red,PlotMarkers->{"□",18}]

Computing Branch 1

Dd=0;i=0

sol[0]=FindRoot[{(1-Dd) (x1-y1)+Dd (-1+x1) (x1 x2-x3/Keq)==0,
(1-Dd) (x2-y2)+Dd (-1+x2) (x1 x2-x3/Keq)==0,y1+y2+y3== 1},{x1,1.2,1.4},{x2,-0.5,-1},{T,275.0,300.0}]//Chop

Eigenvalues[M]/.sol[0]

While[i< 101,i=i+1;{Dd=0.01*i,sol[i]=Chop[FindRoot[{(1-Dd) (x1-y1)+Dd (-1+x1) (x1 x2-x3/Keq)==0,
(1-Dd) (x2-y2)+Dd (-1+x2) (x1 x2-x3/Keq)==0,Psat[1] x1/P+Psat[2] x2/P+Psat[3] x3/P== 1},{x1,x1/.sol[i-1]},{x2,x2/.sol[i-1]},{T,T/.sol[i-1]},MaxIterations-> 100000]],j=0,For[k=1,k≤ 3,k++,If[(Re[Eigenvalues[M][[k]]/.sol[i]])<0,j=j+1]],test=If[j==3,"stable",If[j== 2,"saddle","unstable"]], Print[i," ",sol[i],Dd," ",Eigenvalues[M]/.sol[i]," ",test," ",j]}]

plt5=ListPlot[{Table[{0.01*i,(T/.sol[i])-273.15},{i,0,26,5}],Table[{0.01*i,(T/.sol[i])-273.15},{i,27,100,5}]},PlotRange->{20,40},PlotMarkers->{"○","□"},PlotStyle-> Green]

plt6=ListPlot[{Table[{x1/.sol[i],x2/.sol[i]},{i,0,26,5}],Table[{x1/.sol[i],x2/.sol[i]},{i,27,100,5}]},PlotRange->{{0,1.2},{-0.2,1}},AspectRatio->1,PlotMarkers->{{"○",18},{"□",18}},PlotStyle-> Green]

Computing Branch 2

Dd=0.2;i=0

sol[0]=FindRoot[{(1-Dd) (x1-y1)+Dd (-1+x1) (x1 x2-x3/Keq)==0,
(1-Dd) (x2-y2)+Dd (-1+x2) (x1 x2-x3/Keq)==0,y1+y2+y3== 1},{x1,1.2,1.4},{x2,-0.5,-1},{T,275.0,300.0}]//Chop

Eigenvalues[M]/.sol[0]

While[i< 19,i=i+1;{Dd=0.2-0.01*i,sol[i]=Chop[FindRoot[{(1-Dd) (x1-y1)+Dd (-1+x1) (x1 x2-x3/Keq)==0,
(1-Dd) (x2-y2)+Dd (-1+x2) (x1 x2-x3/Keq)==0,Psat[1] x1/P+Psat[2] x2/P+Psat[3] x3/P== 1},{x1,x1/.sol[i-1]},{x2,x2/.sol[i-1]},{T,T/.sol[i-1]},MaxIterations-> 100000]],j=0,For[k=1,k≤ 3,k++,If[(Re[Eigenvalues[M][[k]]/.sol[i]])<0,j=j+1]],test=If[j==3,"stable",If[j== 2,"saddle","unstable"]], Print[i," ",sol[i],Dd," ",Eigenvalues[M]/.sol[i]," ",test," ",j]}]

Dd=0.2;i=0

solu[0]=FindRoot[{(1-Dd) (x1-y1)+Dd (-1+x1) (x1 x2-x3/Keq)==0,
(1-Dd) (x2-y2)+Dd (-1+x2) (x1 x2-x3/Keq)==0,y1+y2+y3== 1},{x1,1.2,1.4},{x2,-0.5,-1},{T,275.0,300.0}]//Chop

Eigenvalues[M]/.solu[0]

While[i<7 ,i=i+1;{Dd=0.2+0.01*i,solu[i]=Chop[FindRoot[{(1-Dd) (x1-y1)+Dd (-1+x1) (x1 x2-x3/Keq)==0,
(1-Dd) (x2-y2)+Dd (-1+x2) (x1 x2-x3/Keq)==0,Psat[1] x1/P+Psat[2] x2/P+Psat[3] x3/P== 1},{x1,x1/.solu[i-1]},{x2,x2/.solu[i-1]},{T,T/.solu[i-1]},MaxIterations-> 100000]],j=0,For[k=1,k≤ 3,k++,If[(Re[Eigenvalues[M][[k]]/.solu[i]])<0,j=j+1]],test=If[j==3,"stable",If[j== 2,"saddle","unstable"]], Print[i," ",solu[i],Dd," ",Eigenvalues[M]/.solu[i]," ",test," ",j]}]

Dd=0.28;i=0

sol[28]=FindRoot[{(1-Dd) (x1-y1)+Dd (-1+x1) (x1 x2-x3/Keq)==0,
(1-Dd) (x2-y2)+Dd (-1+x2) (x1 x2-x3/Keq)==0,y1+y2+y3== 1},{x1,0.8,0.9},{x2,0.1,0.2},{T,275.0,300.0}]//Chop

Eigenvalues[M]/.sol[28]

i=28;

While[i<100 ,i=i+1;{Dd=0.01*i,sol[i]=Chop[FindRoot[{(1-Dd) (x1-y1)+Dd (-1+x1) (x1 x2-x3/Keq)==0,
(1-Dd) (x2-y2)+Dd (-1+x2) (x1 x2-x3/Keq)==0,Psat[1] x1/P+Psat[2] x2/P+Psat[3] x3/P== 1},{x1,x1/.sol[i-1]},{x2,x2/.sol[i-1]},{T,T/.sol[i-1]},AccuracyGoal->5]],j=0,For[k=1,k≤ 3,k++,If[(Re[Eigenvalues[M][[k]]/.sol[i]])<0,j=j+1]],test=If[j==3,"stable",If[j== 2,"saddle","unstable"]], Print[i," ",sol[i],Dd," ",Eigenvalues[M]/.sol[i]," ",test," ",j]}]

plt7=ListPlot[{Table[{0.2-0.01*i,(T/.sol[i])-273.15},{i,1,19,5}],Table[{0.2+0.01*i,(T/.solu[i])-273.15},{i,1,7,5}],Table[{0.01*i,(T/.sol[i])-273.15},{i,28,100,5}]},PlotRange->{0,40},PlotMarkers->{"□","□","○"},PlotStyle-> Magenta]

plt8=ListPlot[{Table[{x1/.sol[i],x2/.sol[i]},{i,1,15,5}],Table[{x1/.solu[i],x2/.solu[i]},{i,1,7,5}],Table[{x1/.sol[i],x2/.sol[i]},{i,27,100,5}]},PlotRange->All,AspectRatio->1,PlotMarkers->{{"□",16},{"□",16},{"○",16}},PlotStyle-> Magenta]

Bifurcation diagram for the hypothetical ternary mixture of A, B and C with chemical reaction A+B⇄C, (Subscript[K, eq])=5 at 1atm

Show[plt1,plt3,plt5,plt7,PlotRange->{0,80},AspectRatio->1,Frame->True]

Chemical equilibrium curve

plt9=ContourPlot[x1 x2-x3/Keq==0,{x1,0,1},{x2,0,1},ContourStyle-> Red]

Solution branches on a phase plan for the hypothetical ternary mixture of A, B and C with chemical reaction A+B⇄C, (Subscript[K, eq])=5 at 1atm

Show[plt2,plt4,plt6,plt8,plt9,AspectRatio->1,Frame->True,Epilog-> Line[{{0,1},{1,0}}],PlotRange->{{-0.2,1.5},{-0.5,1.2}}]