Systematic Determination of Azeotropes Using Relative Volatility Functions​
​Authors: Housam Binous and Ahmed Bellagi

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].
​
We plot the relative volatility distribution for this system on both ternary equilateral and rectangular diagrams. We also determine all four azeotropes.
​
Reference:
[1] G. M. Wilson, "Vapor-Liquid Equilibrium XI: A New Expression for the Excess Free Energy of Mixing," Journal of the American Chemical Society, 86(2), 1964 pp. 127–130.
The coding below is based on the ideas presented by various contributors. See: https://mathematica.stackexchange.com/questions/39733/how-to-plot-ternary-density-plots​
​
The underlying idea behind plotting the density plot on a ternary equilateral/rectangular diagram is as follows:
1- construct a polygon (the ternary equilateral/rectangular diagram)
2- fill it with an image corresponding to the relative volatility distribution density plot obtained using the interpolating function f.
In[]:=
ClearAll[tickGraphics,triangleTicks];​​Module[{ticklen=0.020,(*lengthofticks(constparameter)*)textoffset=0.05},tickGraphics[tickrange:{0.,total_},angle_]:=With[{rotSideXF=RotationTransform[angle,total{1.,Tan[Pi/6]}/2.]},{Line[rotSideXF/@Thread[{tickrange,0}]],With[{rotTickXF=Composition[RotationTransform[π/3.,rotSideXF@{#,0.}],rotSideXF]},{Text[NumberForm[N@#,{3,2}],rotTickXF[total{-ticklen-textoffset,0}+{#,0.}],{0,0.}],{Line[rotTickXF/@{total{-ticklen,0}+{#,0.},{#,0.}}]}}]&/@Select[N@FindDivisions[tickrange,5],0.<=#<=total&]}];​​triangleTicks[g_Graphics,total_:1]:=Show[Graphics[{First@g}],Graphics[{AbsoluteThickness[0.5],Table[tickGraphics[{0.,total},N@angle],{angle,0,4Pi/3,2Pi/3}]}],Axes->False,Frame->None,PlotRange->total{{0,1+ticklen},{0,Sqrt[3]/2+ticklen/2}},PlotRangePadding->0.07total,PlotRangeClipping->False,ImagePadding->{{20,10},{20,10}}]]
In[]:=
texturize[f_,n_,colf_]:=#/.Polygon[{v1_,v2_,v3_}]:>{EdgeForm[],Texture@ImageData@Colorize[Image@f[#1,#2]&[#,Transpose[#]]&@ConstantArray[Range[0,1,1./n],n+1],ColorFunction->colf,ColorFunctionScaling->True],Polygon[{v3,v2,v1},VertexTextureCoordinates->{{1-1.5/(n+3),1-1.5/(n+3)},{1.5/(n+3),1.5/(n+3)},{1.5/(n+3),1-1.5/(n+3)}}]}&;
Gas Constant (expressed in kcal/K*kmol) and Total Pressure (expressed in KPa)
In[]:=
R=1.987;P=101.325;

Chloroform = component 1
Acetone = component 2
Methanol = component 3

Vapor Pressure Using Antoine Equation : (Temperature T expressed in Kelvin)​​
ln
(
sat
P
i
)=
A
i
+
B
i
C
i
+T
+
D
i
ln(T)+
E
i
 
F
i
T
In[]:=
A1=73.7058;B1=-6055.60;C1=0;D1=-8.91890;E1=7.7440710^-6;F1=2.0;​​A2=71.3031;B2=-5952.00;C2=0.0;D2=-8.53128;E2=7.8239310^-6;F2=2.0;​​A3=59.8373;B3=-6282.89;C3=0;D3=-6.37873;E3=4.6174610^-6;F3=2.0;
In[]:=
PS1=Exp[A1+B1/(C1+T)+D1Log[T]+E1T^F1];​​PS2=Exp[A2+B2/(C2+T)+D2Log[T]+E2T^F2];​​PS3=Exp[A3+B3/(C3+T)+D3Log[T]+E3T^F3];
In[]:=
BP1=-273.15+T/.FindRoot[PS1==P,{T,300}](*b.p.ofChloroform*)
Out[]=
61.0994
In[]:=
BP2=-273.15+T/.FindRoot[PS2==P,{T,300}](*b.p.ofAcetone*)
Out[]=
56.0648
In[]:=
BP3=-273.15+T/.FindRoot[PS3==P,{T,300}](*b.p.ofMethanol*)
Out[]=
64.4832
Binary interaction Parameters and Activities for Wilson Model
In[]:=
l12=-506.85180;l21=116.11710;l13=-361.79440;​​l31=1694.02400;l23=-124.93280;l32=551.45450;
In[]:=
V1=80.7;V2=74.0;V3=40.7;(*molarvolumes*)
In[]:=
A12=V2/V1Exp[-l12/(R*(T))];​​A13=V3/V1Exp[-l13/(R*(T))];​​A32=V2/V3Exp[-l32/(R*(T))];​​A21=V1/V2Exp[-l21/(R*(T))];​​A31=V1/V3Exp[-l31/(R*(T))];​​A23=V3/V2Exp[-l23/(R*(T))];​​​​GAM1=Exp[-Log[x1+x2*A12+x3*A13]+1-(x1/(x1+x2*A12+x3*A13)+x2*A21/(x1*A21+x2+x3*A23)+x3*A31/(x1*A31+x2*A32+x3))];​​GAM2=Exp[-Log[x1*A21+x2+x3*A23]+1-(x1*A12/(x1+x2*A12+x3*A13)+x2/(x1*A21+x2+x3*A23)+x3*A32/(x1*A31+x2*A32+x3))];​​GAM3=Exp[-Log[x1*A31+x2*A32+x3]+1-(x1*A13/(x1+x2*A12+x3*A13)+x2*A23/(x1*A21+x2+x3*A23)+x3/(x1*A31+x2*A32+x3))];
Calculation of the relative volatility distributions for this system at different locations of the ternary diagram.​We solve the following bubble temperature equation :​​
P=
N
C
∑
i=1
x
i 
sat
P
i
γ
i
where
N
C
=3
.​The summation rules indicate that the sum of the mole fractions in either the liquid-phase or the vapor-phase are equal to unity: ​​
N
C
∑
i=1
x
k,i
=1
and
N
C
∑
i=1
y
k,i
=1
​​At low pressure, the equilibrium relationship is given by the modified Raoult’s law:​​
 y
i 
P=
x
i 
sat
P
i
γ
i
​where the Wilson model, which allows the prediction of the activity coefficients, is used. This latter model is valid since we have a homogeneous liquid-phase. ​The Wilson model is given by:​​
ln(
γ
k
)=-ln
N
c
∑
j=1
x
j
A
kj
+1-
N
c
∑
i=1
x
i
A
ik
N
c
∑
j=1
x
j
A
ij
​​with
A
ij
=
v
j
v
i
exp-
λ
ij
-
λ
ii
R T
​​The relative volatilities are defined by:​​
α
ij
=
y
i
x
j
x
i
y
j
fori,j=1to3andi≠j
​​
In[]:=
y1=GAM1x1PS1/P;​​y2=GAM2x2PS2/P;​​y3=GAM3x3PS3/P;
In[]:=
sol[a_,b_,c_]:=FindRoot[PGAM1PS1x1+GAM2PS2x2+GAM3PS2x3/.x1a/.x2b/.x3c,{T,100}]
In[]:=
α
13
[x_,y_]:=If0≤x+y≤1,
y1y
xy3
/.sol[x,1-x-y,y]/.x1x/.x2(1-x-y)/.x3y,0
In[]:=
α
12
[x_,y_]:=If0≤x+y≤1,
y1y
xy2
/.sol[x,y,1-x-y]/.x1x/.x2y/.x3(1-x-y),0
In[]:=
α
23
[x_,y_]:=If0≤x+y≤1,
y2y
xy3
/.sol[1-x-y,x,y]/.x2x/.x3y/.x1(1-x-y),0
This problem is suitable for parallelization.
In[]:=
LaunchKernels[]
Out[]=
KernelObject
Local kernel
KernelID: 1
,KernelObject
Local kernel
KernelID: 2
,KernelObject
Local kernel
KernelID: 3
,KernelObject
Local kernel
KernelID: 4

Building an interpolating function labelled f
α
ij
In[]:=
tbl13=Parallelize@Table[{{x,y},
α
13
[x,y]},{x,0.0001,1,0.005},{y,0.0001,1,0.005}];//Chop//Quiet
In[]:=
fα
13
=Interpolation[Flatten[tbl13,1]];
In[]:=
tbl12=Parallelize@Table[{{x,y},
α
12
[x,y]},{x,0.0001,1,0.005},{y,0.0001,1,0.005}];//Chop//Quiet
In[]:=
fα
12
=Interpolation[Flatten[tbl12,1]];
In[]:=
tbl23=Parallelize@Table[{{x,y},
α
23
[x,y]},{x,0.0001,1,0.005},{y,0.0001,1,0.005}];//Chop//Quiet
In[]:=
fα
23
=Interpolation[Flatten[tbl23,1]];
We find three binary azeotropes and a ternary azeotrope using the Relative Volatility Functions defined above (i.e.,
α
12
,
α
13
and
α
23
).​​This mixture presents four azeotropes, whose locations and boiling points are:​​
mole%
BP°C
chloroform
methanol
acetone
A
1
65.52
34.48
-
53.82
A
2
-
20.93
79.07
55.35
A
3
66.13
-
33.21
64.45
A
4
22.90
43.89
33.21
57.32
​​