All Azeotropes for a Ternary System by bifurcation : Acetone - Chloroform - Methanol System​
​Authors: Housam Binous and Ahmed Bellagi

Acknowledgement: An analogous problem and its Matlab solution were developed and refined in the Chemical Engineering department at MIT by several people, including the late Prof. James Swan. Here, we provide the solution using Mathematica for a different ternary system.​Let us consider a ternary mixture of chloroform, acetone, and methanol at
101.325kPa
. The vapor-liquid equilibrium (VLE) behavior is described by the modified Raoult's law with activity coefficients predicted by the Wilson model [1]. ​This mixture presents seven azeotropes, whose locations and boiling points are given below:​​
​​ ​This problem was treated by Fidkowski, Malone and Doherty (Ref. 1) and their results are reproduced below:​​
​​​References:​[1] Z.T. Fidkowski, M.F. Malone, M.F. Doherty, Computing azeotropes in multicomponent mixtures, Computers & Chemical Engineering, Volume 17, Issue 12, 1993, Pages 1141-1155, https://doi.org/10.1016/0098-1354(93)80095-5​​​​Chloroform = component 1Acetone = component 2Methanol = component 3
Total pressure (1 atm) in torr
In[]:=
P=760;
Wilson model
In[]:=
l12=-506.85180;l21=116.11710;l13=-361.79440;​​l31=1694.02400;l23=-124.93280;l32=551.45450;​​lambdas=Transpose@{{0,l12,l13},{l21,0,l23},{l31,l32,0}};(*incal/mol*)​​​​volumes={80.7,74.0,40.7};(*inmL/mol*)​​​​wilson[x_,T_]:=Module[{Lambdas,gammas,i,j},​​For[i=1,i≤3,i++,For[j=1,j≤3,j++,​​Lambdas[i,j]=volumes[[j]]/volumes[[i]]*Exp[-lambdas[[j,i]]/(1.9872041*T)]]];​​For[i=1,i≤3,i++,gammas[i]=1-Log[Sum[Lambdas[i,k]x[[k]],{k,1,3}]];​​For[j=1,j≤3,j++,gammas[i]=gammas[i]-Lambdas[j,i]x[[j]]/Sum[Lambdas[j,k]x[[k]],{k,1,3}]]];​​Table[Exp[gammas[i]],{i,1,3}]]
Vapor pressure obtained using Antoine equation
In[]:=
A2=7.11714;B2=1210.595;C2=229.664;​​A1=6.95465;B1=1170.966;C1=226.232;​​A3=8.08097;B3=1582.271;C3=239.726;​​myPsats={10^(A1-B1/(C1+(T-273.15))),10^(A2-B2/(C2+(T-273.15))),10^(A3-B3/(C3+(T-273.15)))};
System of algebraic equations, which solutions yield the azeotropes compositions and boiling point when λ=1
In[]:=
mysyshom[x_,T_,l_]:=Module[{gammas,psat},gammas=wilson[x,T];psat=myPsats;{x[[1]]-((1-l)+lgammas[[1]])psat[[1]]x[[1]]/P0,x[[2]]-((1-l)+lgammas[[2]])psat[[2]]x[[2]]/P0,x[[3]]-((1-l)+lgammas[[3]])psat[[3]]x[[3]]/P0,x[[1]]+x[[2]]+x[[3]]-10}]​​
Pure components' boiling point versus λ
In[]:=
l=Table[i,{i,0,1,0.01}];data=.;​​guess={1,10^-12,10^-12,300};​​data[1]=Table[{l[[i]],-273.15+T/.FindRoot[mysyshom[{x1,x2,x3},T,l[[i]]],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{T,guess[[4]]}]},{i,1,101}];​​guess={10^-12,1,10^-12,300};​​data[2]=Table[{l[[i]],-273.15+T/.FindRoot[mysyshom[{x1,x2,x3},T,l[[i]]],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{T,guess[[4]]}]},{i,1,101}];​​guess={10^-12,10^-12,1,300};​​data[3]=Table[{l[[i]],-273.15+T/.FindRoot[mysyshom[{x1,x2,x3},T,l[[i]]],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{T,guess[[4]]}]},{i,1,101}];p1=ListPlot[{data[1],data[2],data[3]},GridLinesAutomatic,FrameTrue,AspectRatio1,Joined->True,PlotRange{{0,1},{50,70}},FrameLabel{Style["λ",16],Style["T",16]}]
Out[]=
Compute the Jacobian matrix
In[]:=
x={x1,x2,x3};​​l=.;(*lisλ*)​​gammas=wilson[x,T];psat=myPsats;​​jac=Outer[D,{x[[1]]-((1-l)+lgammas[[1]])psat[[1]]x[[1]]/P,x[[2]]-((1-l)+lgammas[[2]])psat[[2]]x[[2]]/P,x[[3]]-((1-l)+lgammas[[3]])psat[[3]]x[[3]]/P,x[[1]]+x[[2]]+x[[3]]-1},{x1,x2,x3,T}];
Starting with the pure component mixture and track bifurcations

1st branch

In[]:=
j=0;guessI={10^-12,1.0,10^-12,350};While[j≤501,k=1;guess=If[j0,guessI,guess];s=Quiet@FindRoot[mysyshom[{x1,x2,x3},T,(j)0.002],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{T,guess[[4]]},MaxIterations10000];sol[j,1,1]=x1/.s;sol[j,2,1]=x2/.s;sol[j,3,1]=x3/.s;sol[j,4,1]=T/.s;​​jacs[j]=jac/.x1sol[j,1,1]/.x2sol[j,2,1]/.x3sol[j,3,1]/.Tsol[j,4,1]/.l(j)0.002;​​If[j>2&&Sign[Det[jacs[j]]Det[jacs[j-1]]]<0,Print[j];eval=Eigenvalues[jacs[j]];​​pos=First@Position[Abs@eval,Min[Abs@eval]];evect=Flatten[Re[Eigenvectors[jacs[j]][[pos]]]];​​diff={0,0,0};inc=Exp@0.01;guess={x1,x2,x3,T}/.x1sol[j,1,1]/.x2sol[j,2,1]/.x3sol[j,3,1]/.Tsol[j,4,1];Print[Sign[Det[jacs[j]]Det[jacs[j-1]]]," ",pos," ",guess];​​k=1;​​While[Abs[diff[[1]]]<10^-3&&Abs[diff[[2]]]<10^-3&&Abs[diff[[3]]]<10^-3,{guess=Re@(guess+(incevect));​​s=Quiet@FindRoot[mysyshom[{x1,x2,x3},T,(j+1)0.002],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{T,guess[[4]]},MaxIterations10000];​​nbw[1]=x1/.s;nbw[2]=x2/.s;nbw[3]=x3/.s;nbw[4]=T/.s;inc=-1.5inc;​​diff=Table[sol[j,i,1]-nbw[i],{i,1,3}],k++}];Table[sol[j,i,1]=nbw[i],{i,1,4}];];If[j≥2&&k≥2,Table[guess[[i]]=nbw[i],{i,1,4}],Table[guess[[i]]=sol[j,i,1],{i,1,4}]];j++]
260
-1 {4} {1.17502×
-25
10
,1.,3.88797×
-27
10
,329.251}
261
-1 {4} {-8.62817×
-32
10
,0.998016,0.00198375,329.251}
In[]:=
dataaz1=Table[{j0.002,-273.15+sol[j,4,1]},{j,0,500,20}];
In[]:=
j=500;
In[]:=
d1={"AZ1",sol[j,1,1],sol[j,2,1],sol[j,3,1],sol[j,4,1]-273.15}//Chop
Out[]=
{AZ1,0,0.791777,0.208223,55.3891}
In[]:=
paz1=ListPlot[dataaz1,GridLinesAutomatic,FrameTrue,AspectRatio1,PlotMarkersStyle["■",Darker@Orange,12],PlotRange{{0,1},{70,100}}];

2nd branch

In[]:=
Clear[j,guessI,s,sol,jacs,eval,evect,pos,diff,inc,nbw,k]
In[]:=
j=0;guessI={1.0,10^-12,10^-12,350};While[j≤501,k=1;guess=If[j0,guessI,guess];s=Quiet@FindRoot[mysyshom[{x1,x2,x3},T,(j)0.002],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{T,guess[[4]]},MaxIterations10000];sol[j,1,1]=x1/.s;sol[j,2,1]=x2/.s;sol[j,3,1]=x3/.s;sol[j,4,1]=T/.s;​​jacs[j]=jac/.x1sol[j,1,1]/.x2sol[j,2,1]/.x3sol[j,3,1]/.Tsol[j,4,1]/.l(j)0.002;​​If[j>=2&&Sign[Det[jacs[j]]Det[jacs[j-1]]]<0,Print[j];eval=Eigenvalues[jacs[j]];​​pos=First@Position[Abs@eval,Min[Abs@eval]];evect=Flatten[Re[Eigenvectors[jacs[j]][[pos]]]];​​diff={0,0,0};inc=0.01;guess={x1,x2,x3,T}/.x1sol[j,1,1]/.x2sol[j,2,1]/.x3sol[j,3,1]/.Tsol[j,4,1];Print[Sign[Det[jacs[j]]Det[jacs[j-1]]]," ",pos," ",guess];​​k=1;​​While[Abs[diff[[1]]]<10^-3&&Abs[diff[[2]]]<10^-3&&Abs[diff[[3]]]<10^-3,{guess=Re@(guess+(incevect));​​s=Quiet@FindRoot[mysyshom[{x1,x2,x3},T,(j+1)0.002],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{T,guess[[4]]},MaxIterations10000];​​nbw[1]=x1/.s;nbw[2]=x2/.s;nbw[3]=x3/.s;nbw[4]=T/.s;inc=-1.5inc;​​diff=Table[sol[j,i,1]-nbw[i],{i,1,3}],k++}];Table[sol[j,i,1]=nbw[i],{i,1,4}];];If[j≥2&&k≥2,Table[guess[[i]]=nbw[i],{i,1,4}],Table[guess[[i]]=sol[j,i,1],{i,1,4}]];j++]
12
-1 {4} {1.,8.27748×
-197
10
,-4.71551×
-187
10
,334.354}
13
-1 {4} {0.987925,4.12234×
-228
10
,0.0120753,334.353}
14
-1 {4} {1.,-1.32239×
-59
10
,1.35041×
-34
10
,334.354}
15
-1 {4} {0.975451,7.05682×
-90
10
,0.024549,334.351}
107
-1 {4} {0.747981,1.96152×
-27
10
,0.252019,333.282}
108
-1 {4} {0.738587,0.00770417,0.253709,333.266}

3rd branch

4th branch

Pure components' boiling points
All seven azeotropes
Comparison with Matlab output