Consider a ternary mixture of Ethanol-MEK-Water at . 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 four azeotropes (3 pure components, 3 binary mixtures, and 1 ternary mixture), whose locations and boiling points are as follows :We propose to find all the azeotropes of this ternary system by starting with the pure component mixture and detecting bifurcations.
101.325kPa
This problem was treated by Fidkowski, Malone and Doherty (Ref. 1) and their results are reproduced below:
Wilson model
lambdas={{0,-149.7978,926.2630},{694.0825,0,1888.8509},{393.1971,6811.3433,0}};(*incal/mol*)volumes={56.68,90.17,18.07};(*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
myPsats[T_]:={10.^(8.11220-1592.864/(226.184+(T-273.15))),10.^(7.06356-1261.339/(221.969+(T-273.15))),10.^(8.07131-1730.630/(233.426+(T-273.15)))}
System of algebraic equations, which solutions yield the azeotropes compositions and boiling point
mysys[x_,T_]:=Module[{gammas,psat},gammas=wilson[x,T];psat=myPsats[T];{x[[1]]-gammas[[1]]psat[[1]]x[[1]]/P0,x[[2]]-gammas[[2]]psat[[2]]x[[2]]/P0,x[[3]]-gammas[[3]]psat[[3]]x[[3]]/P0,x[[1]]+x[[2]]+x[[3]]-10}]
System of algebraic equations, which solutions yield the azeotropes compositions and boiling point
mysyshom[x_,T_,l_]:=Module[{gammas,psat},gammas=wilson[x,T];psat=myPsats[T];{x[[1]]-((1-l)+lgammas[[1]])psat[[1]]x[[1]]/P0,x[[2]]-((1-l)+lgammas[[2]])psat[[2]]x[[2]]/P0,x[[3]]-((1-l)+lgammas[[3]])psat[[3]]x[[3]]/P0,x[[1]]+x[[2]]+x[[3]]-10}]
Total pressure (1 atm) in torr
P=760;
Pure components' boiling point versus λ
In[]:=
l=Table[i,{i,0,1,0.01}];data=.;guess=Log@{1,10^-8,10^-8,350};data[1]=Quiet@Table[{l[[i]],-273.15+Exp@T/.Quiet@FindRoot[mysyshom[Exp@{x1,x2,x3},Exp@T,l[[i]]],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{T,guess[[4]]}]},{i,1,101}];guess=Log@{10^-8,1,10^-8,350};data[2]=Quiet@Table[{l[[i]],-273.15+Exp@T/.Quiet@FindRoot[mysyshom[Exp@{x1,x2,x3},Exp@T,l[[i]]],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{T,guess[[4]]}]},{i,1,101}];guess=Log@{10^-8,10^-8,1,350};data[3]=Quiet@Table[{l[[i]],-273.15+Exp@T/.Quiet@FindRoot[mysyshom[Exp@{x1,x2,x3},Exp@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]},GridLinesAutomatic,FrameTrue,AspectRatio1,Joined->True,PlotRange{{0,1},{70,105}},FrameLabel{Style["λ",16],Style["T",16]},PlotStyle->Thick]
Out[]=
In[]:=
p1=ListPlot[{data[1],data[2],data[3]},GridLinesAutomatic,FrameTrue,AspectRatio1,Joined->True,PlotRange{{0,1},{70,100}},FrameLabel{Style["λ",16],Style["T",16]}];
Compute the Jacobian matrix
x={x1,x2,x3};l=.;gammas=wilson[x,T];psat=myPsats[T];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}];
Initial guesses
In[]:=
guessT={{1,10^-12,10^-12,350},{1,10^-12,10^-12,350},{10^-12,1.0,10^-12,350},{1.0,10^-12,10^-12,350}};
Labels
In[]:=
AZ={"AZ1","AZ2","AZ3","AZ4"}
Out[]=
{AZ1,AZ2,AZ3,AZ4}
Starting with the pure component mixture and track bifurcations
In[]:=
i=1;While[i<=4,j=1;guessI=guessT[[i]];While[j≤501,k=1;a=If[i==2,25,2];guess=If[j1,Log@guessI,Log@guess];s=Quiet@FindRoot[mysyshom[Exp@{x1,x2,x3},Exp@T,(j)0.002],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{T,guess[[4]]},MaxIterations10000];sol[j,1,1]=Exp@x1/.s;sol[j,2,1]=Exp@x2/.s;sol[j,3,1]=Exp@x3/.s;sol[j,4,1]=Exp@T/.s;jacs[j]=jac/.x1sol[j,1,1]/.x2sol[j,2,1]/.x3sol[j,3,1]/.Tsol[j,4,1]/.l(j)0.002;If[j>a&&Sign[Det[jacs[j]]Det[jacs[j-1]]]<0,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.1;guess={x1,x2,x3,T}/.x1sol[j,1,1]/.x2sol[j,2,1]/.x3sol[j,3,1]/.Tsol[j,4,1];k=1;While[Abs[diff[[1]]]<10^-3&&Abs[diff[[2]]]<10^-3&&Abs[diff[[3]]]<10^-3,{guess=Re@Log@(guess+(incevect));s=Quiet@FindRoot[mysyshom[Exp@{x1,x2,x3},Exp@T,(j+1)0.002],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{T,guess[[4]]},MaxIterations10000];nbw[1]=Exp@x1/.s;nbw[2]=Exp@x2/.s;nbw[3]=Exp@x3/.s;nbw[4]=Exp@T/.s;inc=-1.5inc;diff=Table[sol[j,i,1]-nbw[i],{i,1,3}],k++}];If[i==4,Table[sol[j+1,i,1]=nbw[i],{i,1,4}],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}]];If[i==4,If[236≥j≥2&&k≥2,Table[guess[[i]]=nbw[i],{i,1,4}],Table[guess[[i]]=sol[j,i,1],{i,1,4}]]];j++];dataaz[i]=Table[{j0.002,-273.15+sol[j,4,1]},{j,1,500,20}];j=500;d[i]={AZ[[i]],sol[j,1,1],sol[j,2,1],sol[j,3,1],sol[j,4,1]};paz1[i]=ListPlot[dataaz[i],GridLinesAutomatic,FrameTrue,AspectRatio1,PlotMarkersStyle["◼",Hue[i/5],12],PlotRange{{0,1},{70,100}}];Print[i];i++];
1
2
3
4
Pure components' boiling points
In[]:=
TE=T/.FindRoot[myPsats[T][[1]]P,{T,350}]
Out[]=
351.448
In[]:=
TMEK=T/.FindRoot[myPsats[T][[2]]P,{T,350}]
Out[]=
352.739
In[]:=
TW=T/.FindRoot[myPsats[T][[3]]P,{T,350}]
Out[]=
373.147
In[]:=
Eth={"E",1.0,0.0,0.0,TE};MEK={"M",0.0,1.0,0.0,TMEK};Water={"W",0.0,0.0,1.0,TW};
All seven azeotropes
In[]:=
Text@Style[Grid[{Eth,MEK,Water,d[1],d[2],d[3],d[4]}//Chop,Frame->All,Background->{None,{Lighter[Blend[{Blue,Red}],.8],Lighter[Blend[{Blue,Red}],.8],Lighter[Blend[{Blue,Red}],.8],{Lighter[Blend[{Yellow,Darker@Green}],.2],Lighter[Green,.5]}}},ItemSize->{1.2{5,4.2,4.2,4.2,5},2},Frame->Darker[Gray,.6],ItemStyle->12,Spacings->{1.,.5}]]
Comparison with Matlab output
Acknowledgement
Acknowledgement
This 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. Problem statement can be found in HW#3 Problem 1 at : https://ocw.mit.edu/courses/10-34-numerical-methods-applied-to-chemical-engineering-fall-2015/pages/assignments/
References
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
Volume 17, Issue 12, 1993, Pages 1141-1155, https://doi.org/10.1016/0098-1354(93)80095-5
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
All azeotropes for ternary system: Ethanol - MEK - water system using bifurcation tracking
by Housam Binous and Ahmed Bellagi
Wolfram Community, STAFF PICKS, June 18, 2025
https://community.wolfram.com/groups/-/m/t/3481435
by Housam Binous and Ahmed Bellagi
Wolfram Community, STAFF PICKS, June 18, 2025
https://community.wolfram.com/groups/-/m/t/3481435