We propose to find all azeotropes for a quaternary mixture composed of Methanol, Benzene, i-Propanol and n-Propanol at 1 atm using the Wilson model.This problem was first treated by Fidkowski, Malone and Doherty (Ref. 1) and their results are reproduced below:Finally, this problem was also investigated by the late Prof. Floudas and co-workers using GAMS: see Chapter 14 at : http://titan.princeton.edu/TestProblems/
We benchmark our results, summarized in the Table below, against data generated by Aspen-Properties as well as a Matlab bifurcation code written by several people at the CHE department of MIT including the late Prof. Swan.
In[]:=
(*vaporpressureofallcomponentsusingAnotineEquationTisinKelvinandpessureisinmmHg*)aa[1]=8.08097;bb[1]=1582.271;cc[1]=239.726;aa[2]=6.87987;bb[2]=1196.760;cc[2]=219.161;aa[3]=8.87829;bb[3]=2010.330;cc[3]=252.636;aa[4]=8.37895;bb[4]=1788.020;cc[4]=227.438;p1[T_]:=10^(aa[1]-bb[1]/(cc[1]+T));(*Methanol*)p2[T_]:=10^(aa[2]-bb[2]/(cc[2]+T));(*Benzene*)p3[T_]:=10^(aa[3]-bb[3]/(cc[3]+T));(*i-Propanol*)p4[T_]:=10^(aa[4]-bb[4]/(cc[4]+T));(*n-Propanol*)
In[]:=
P=760;(*totalpressureinmmHg*)
In[]:=
(*normalboilingtemperatuesofMethanol,Benzene,i-Propanolandn-Propanol*)
In[]:=
TM=T/.FindRoot[p1[T]P,{T,60}]
Out[]=
64.5477
In[]:=
TB=T/.FindRoot[p2[T]P,{T,60}]
Out[]=
80.0996
In[]:=
TiP=T/.FindRoot[p3[T]P,{T,60}]
Out[]=
82.56
In[]:=
TnP=T/.FindRoot[p4[T]P,{T,60}]
Out[]=
97.7667
In[]:=
R=1.987;(*universalgasconstantinkcal/kmol*)
In[]:=
(*parametersofWilsonmodelarereproducedfromFlouda'stestproblem2insection:LocatingAllHomogeneousAzeotropes*)lamda[i_,i_]:=0;lamda[1,2]=1734.4181;lamda[2,1]=183.0383;lamda[1,3]=-38.4124;lamda[3,1]=176.3185;lamda[1,4]=150.2057;lamda[4,1]=-33.8819;lamda[2,3]=392.7020;lamda[3,2]=813.9198;lamda[2,4]=312.1615;lamda[4,2]=996.2753;lamda[3,4]=818.3291;lamda[4,3]=-481.0590;(*molarvolumes*)v[1]=38.8196;v[2]=94.0952;v[3]=79.6109;v[4]=77.8926;wilson[x_,T_]:=Module[{Lambdas,gammas,i,j},For[i=1,i≤4,i++,For[j=1,j≤4,j++,Lambdas[i,j]=v[j]/v[i]Exp[-lamda[i,j]/(T+273.15)/R]]];For[i=1,i≤4,i++,gammas[i]=1-Log[Sum[Lambdas[i,k]x[[k]],{k,1,4}]];For[j=1,j≤4,j++,gammas[i]=gammas[i]-Lambdas[j,i]x[[j]]/Sum[Lambdas[j,k]x[[k]],{k,1,4}]]];Table[Exp[gammas[i]],{i,1,4}]]
In[]:=
myPsats[T_]:={p1[T],p2[T],p3[T],p4[T]};
In[]:=
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[[4]]-((1-l)+lgammas[[4]])psat[[4]]x[[4]]/P0,x[[1]]+x[[2]]+x[[3]]+x[[4]]-10}]
Approach based on changing eaach time the initial guess in order to find the azeotropes . This approach compares favorably with results found in Ref . 1 .
In[]:=
guess={0.6,0.4,0.0,0.0,60};FindRoot[mysyshom[{x1,x2,x3,x4},T,1],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{x4,guess[[4]]},{T,guess[[5]]}]
Out[]=
{x10.620093,x20.379907,x30.,x40.,T57.9336}
In[]:=
guess={0.0,0.6,0.4,0.0,60};FindRoot[mysyshom[{x1,x2,x3,x4},T,1],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{x4,guess[[4]]},{T,guess[[5]]}]
Out[]=
{x1-8.36736×,x20.586058,x30.413942,x40.,T71.968}
-37
10
In[]:=
guess={0.0,0.4,0.0,0.5,60};FindRoot[mysyshom[{x1,x2,x3,x4},T,1],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{x4,guess[[4]]},{T,guess[[5]]}]
Out[]=
{x17.69569×,x20.782293,x31.28092×,x40.217707,T77.0483}
-31
10
-24
10
Pure components' boiling point versus λ
In[]:=
l=Table[i,{i,0,1,0.01}];data=.;guess={1,10^-12,10^-12,10^-12,70};data[1]=Table[{l[[i]],T/.Quiet@FindRoot[mysyshom[{x1,x2,x3,x4},T,l[[i]]],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{x4,guess[[4]]},{T,guess[[5]]}]},{i,1,101}];guess={10^-12,1,10^-12,10^-12,70};data[2]=Table[{l[[i]],T/.Quiet@FindRoot[mysyshom[{x1,x2,x3,x4},T,l[[i]]],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{x4,guess[[4]]},{T,guess[[5]]}]},{i,1,101}];guess={10^-12,10^-12,1,10^-12,70};data[3]=Table[{l[[i]],T/.Quiet@FindRoot[mysyshom[{x1,x2,x3,x4},T,l[[i]]],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{x4,guess[[4]]},{T,guess[[5]]}]},{i,1,101}];guess={10^-12,10^-12,10^-12,1,70};data[4]=Table[{l[[i]],T/.Quiet@FindRoot[mysyshom[{x1,x2,x3,x4},T,l[[i]]],{x1,guess[[1]]},{x2,guess[[2]]},{x3,guess[[3]]},{x4,guess[[4]]},{T,guess[[5]]}]},{i,1,101}];
In[]:=
plt=ListPlot[{data[1],data[2],data[3],data[4]},GridLinesAutomatic,FrameTrue,AspectRatio1,Joined->True,PlotRange{{0,1},{60,100}},FrameLabel{Style["λ",16],Style["T in Kelvin",16]}]
Out[]=
Bifurcation Approach inspired from Matlab code by CHE Depart. of MIT
Bifurcation Approach inspired from Matlab code by CHE Depart. of MIT
In[]:=
x={x1,x2,x3,x4};l=.;(*lisλ*)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[[4]]-((1-l)+lgammas[[4]])psat[[4]]x[[4]]/P,x[[1]]+x[[2]]+x[[3]]+x[[4]]-1},{x1,x2,x3,x4,T}];
2nd branch
2nd branch
3rd branch
3rd branch
Comparison with output from MIT code (created by several peoples at the CHE Dept. of MIT including the late Prof. Swan)
Comparison with output from MIT code (created by several peoples at the CHE Dept. of MIT including the late Prof. Swan)
MATLAB OUTPUT (Code created by several peoples at the CHE Dept . of MIT including the late Prof . Swan)
All three azeotropes
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
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
All azeotropes of a quaternary system by bifurcation
by Housam Binous and Ahmed Bellagi
Wolfram Community, STAFF PICKS, July 7, 2025
https://community.wolfram.com/groups/-/m/t/3492340
by Housam Binous and Ahmed Bellagi
Wolfram Community, STAFF PICKS, July 7, 2025
https://community.wolfram.com/groups/-/m/t/3492340