All Azeotropes of a System Composed of Ethanol - Methylcyclopentane - Benzene - Hexane at 1 atm​
​Authors: Housam Binous and Ahmed Bellagi

We propose to find all azeotropes for a quaternary mixture composed of Ethanol - Methylcyclopentane - Benzene - Hexane at 1 atm using the UNIQUAC model.​​This problem was also investigated by the late Prof. Floudas and co-workers using GAMS: see Test problem 7 of Section Locating All Homogeneous Azeotropes in 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.11220;bb[1]=1592.864;cc[1]=226.184;​​aa[3]=6.8664;bb[3]=1188.05;cc[3]=226.276;​​aa[2]=6.87987;bb[2]=1196.760;cc[2]=219.161;​​aa[4]=6.91058;bb[4]=1189.64;cc[4]=226.28;​​​​p1[T_]:=10^(aa[1]-bb[1]/(cc[1]+T-273.15));(*Ethanol*)​​p2[T_]:=10^(aa[2]-bb[2]/(cc[2]+T-273.15));(*Benzene*)​​p3[T_]:=10^(aa[3]-bb[3]/(cc[3]+T-273.15));(*Methylcyclopentane*)​​p4[T_]:=10^(aa[4]-bb[4]/(cc[4]+T-273.15));(*Hexane*)
In[]:=
P=760;(*totalpressureinmmHg*)
normal boiling temperatures of Ethanol - Methylcyclopentane - Benzene - Hexane
In[]:=
TE=T/.FindRoot[p1[T]P,{T,360}]
Out[]=
351.448
In[]:=
TB=T/.FindRoot[p2[T]P,{T,360}]
Out[]=
353.25
In[]:=
TMcp=T/.FindRoot[p3[T]P,{T,360}]
Out[]=
344.961
In[]:=
TH=T/.FindRoot[p4[T]P,{T,360}]
Out[]=
342.083
In[]:=
(*parametersofUNIQUACmodelarereproducedfromAspen-Properties*)​​​​DATA1={{-0.0464,0.,0.,0.,1.2281,0.},{0.4665,0.,0.,0.,-1.1078,0.},{76.5759,89.8657,65.0854,-9.9976,-385.704,121.71},{-556.475,-598.883,-534.892,-44.1416,294.108,-137.655}};​​​​c1=DATA1[[1]];​​c2=DATA1[[2]];​​d1=DATA1[[3]];​​d2=DATA1[[4]];​​​​A[1,2]=c1[[1]];​​A[1,3]=c1[[2]];​​A[1,4]=c1[[3]];​​A[2,3]=c1[[4]];​​A[2,4]=c1[[5]];​​A[3,4]=c1[[6]];​​A[2,1]=c2[[1]];​​A[3,1]=c2[[2]];​​A[4,1]=c2[[3]];​​A[3,2]=c2[[4]];​​A[4,2]=c2[[5]];​​A[4,3]=c2[[6]];​​​​B[1,2]=d1[[1]];​​B[1,3]=d1[[2]];​​B[1,4]=d1[[3]];​​B[2,3]=d1[[4]];​​B[2,4]=d1[[5]];​​B[3,4]=d1[[6]];​​B[2,1]=d2[[1]];​​B[3,1]=d2[[2]];​​B[4,1]=d2[[3]];​​B[3,2]=d2[[4]];​​B[4,2]=d2[[5]];​​B[4,3]=d2[[6]];​​​​A[i_,i_]:=0;B[i_,i_]:=0;​​​​r[1]=2.11;​​r[3]=3.97;​​r[2]=3.19;​​r[4]=4.50;​​​​q[1]=1.97;​​q[3]=3.01;​​q[2]=2.40;​​q[4]=3.86;​​​​qp[1]=1.97;​​qp[3]=3.01;​​qp[2]=2.40;​​qp[4]=3.86;​​​​uniquac[x_,T_]:=Module{tau,vp,v,GamC,GamR,z,l,phi},​​tau[i_,j_]:=Exp[A[i,j]+(B[i,j])/(T)];​​GamC[i_]:=ExpLog[phi[i]/x[[i]]]+z2q[i]Log[v[i]/phi[i]]+l[i]-phi[i]x[[i]]
4
∑
j=1
x[[j]]l[j];l[i_]:=z2(r[i]-q[i])-(r[i]-1);​​​​z=10;​​GamR[i_]:=Expqp[i]1-Log
4
∑
j=1
vp[j]tau[j,i]-
4
∑
j=1
vp[j]tau[i,j]
4
∑
k=1
vp[k]tau[k,j]
;​​phi[i_]:=r[i]x[[i]]
4
∑
j=1
r[j]x[[j]];​​v[i_]:=q[i]x[[i]]
4
∑
j=1
q[j]x[[j]];​​vp[i_]:=qp[i]x[[i]]
4
∑
j=1
qp[j]x[[j]];​​G[1]=GamC[1]GamR[1];​​G[2]=GamC[2]GamR[2];​​G[3]=GamC[3]GamR[3];​​G[4]=GamC[4]GamR[4];​​Table[G[i],{i,1,4}];
In[]:=
myPsats[T_]:={p1[T],p2[T],p3[T],p4[T]};​​
In[]:=
mysyshom[x_,T_,l_]:=Module[{gammas,psat},gammas=uniquac[x,T];psat=myPsats[T];{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[[4]]-((1-l)+lgammas[[4]])psat[[4]]x[[4]]/P0,x[[1]]+x[[2]]+x[[3]]+x[[4]]-10}]​​
Approach based on changing each time the initial guess in order to find the azeotropes. This approach compares favorably with results given below and found using Aspen-Properties.
Pure components' boiling point versus λ

Bifurcation Approach inspired from Matlab code by CHE Depart. of MIT

Comparison with output from MIT code (created by several peoples at the CHE Dept. of MIT including the late Prof. Swan)

All seven azeotropes