Residue Curve Map for Ternary Systems : Case of Acetone - Chloroform - Methanol System​
​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]. ​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
​​We plot a few residue curves for this system on a ternary equilateral diagram.​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.[2] H. Binous, M. A. Al-Harthi and A. Bellagi, Experimental and Theoretical Study of Multicomponent Batch Distillation, Computer Applications in Engineering Education Journal, 23(6), pp. 813-823, 2015.

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

Gas Constant (expressed in kcal/K*kmol) and Total Pressure (expressed in KPa)
In[]:=
R=1.987;P=101.325;
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[]:=
PSat1[t_]:=Exp[A1+B1/(C1+T[t])+D1Log[T[t]]+E1T[t]^F1];​​PSat2[t_]:=Exp[A2+B2/(C2+T[t])+D2Log[T[t]]+E2T[t]^F2];​​PSat3[t_]:=Exp[A3+B3/(C3+T[t])+D3Log[T[t]]+E3T[t]^F3];
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[t]))];​​A13=V3/V1Exp[-l13/(R*(T[t]))];​​A32=V2/V3Exp[-l32/(R*(T[t]))];​​A21=V1/V2Exp[-l21/(R*(T[t]))];​​A31=V1/V3Exp[-l31/(R*(T[t]))];​​A23=V3/V2Exp[-l23/(R*(T[t]))];
In[]:=
GAM1[i_]:=Exp[-Log[x1[i]+x2[i]*A12+x3[i]*A13]+1-(x1[i]/(x1[i]+x2[i]*A12+x3[i]*A13)+x2[i]*A21/(x1[i]*A21+x2[i]+x3[i]*A23)+x3[i]*A31/(x1[i]*A31+x2[i]*A32+x3[i]))];​​GAM2[i_]:=Exp[-Log[x1[i]*A21+x2[i]+x3[i]*A23]+1-(x1[i]*A12/(x1[i]+x2[i]*A12+x3[i]*A13)+x2[i]/(x1[i]*A21+x2[i]+x3[i]*A23)+x3[i]*A32/(x1[i]*A31+x2[i]*A32+x3[i]))];​​GAM3[i_]:=Exp[-Log[x1[i]*A31+x2[i]*A32+x3[i]]+1-(x1[i]*A13/(x1[i]+x2[i]*A12+x3[i]*A13)+x2[i]*A23/(x1[i]*A21+x2[i]+x3[i]*A23)+x3[i]/(x1[i]*A31+x2[i]*A32+x3[i]))];
Equilibrium relations
In[]:=
y1[t_]=x1[t]PSat1[t]GAM1[t]/P;​​y2[t_]=x2[t]PSat2[t]GAM2[t]/P;​​y3[t_]=x3[t]PSat3[t]GAM3[t]/P;
Summation rule
In[]:=
x3[t_]:=1-x1[t]-x2[t];
Calculation of Residue Curve Map for different initial conditions.​We solve the following DAE system :​​
d
x
i
dξ
=
x
i
-
y
i
with i = 1..3 and ξ is the warped time (i.e., a deformed time).​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
​​where
N
C
=3
.​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
​​For more information, see:​M. F. Doherty and M. F. Malone, Conceptual Design of Distillation Systems, New York: McGraw–Hill, 2001.​
In[]:=
sol1=Quiet@NDSolve[{D[x1[t],{t,1}]y1[t]-x1[t],D[x2[t],{t,1}]y2[t]-x2[t],(x1[t]PSat1[t]GAM1[t]+x2[t]PSat2[t]GAM2[t]+x3[t]PSat3[t]GAM3[t]-P)0,T[0]340,x1[0]0.2,x2[0]0.2},{x1,x2,T},{t,-80,80}];
In[]:=
sol2=Quiet@NDSolve[{D[x1[t],{t,1}]y1[t]-x1[t],D[x2[t],{t,1}]y2[t]-x2[t],(x1[t]PSat1[t]GAM1[t]+x2[t]PSat2[t]GAM2[t]+x3[t]PSat3[t]GAM3[t]-P)0,T[0]340,x1[0]0.5,x2[0]0.1},{x1,x2,T},{t,-80,80}];
In[]:=
sol3=Quiet@NDSolve[{D[x1[t],{t,1}]y1[t]-x1[t],D[x2[t],{t,1}]y2[t]-x2[t],(x1[t]PSat1[t]GAM1[t]+x2[t]PSat2[t]GAM2[t]+x3[t]PSat3[t]GAM3[t]-P)0,T[0]340,x1[0]0.05,x2[0]0.6},{x1,x2,T},{t,-80,80}];
In[]:=
sol4=Quiet@NDSolve[{D[x1[t],{t,1}]y1[t]-x1[t],D[x2[t],{t,1}]y2[t]-x2[t],(x1[t]PSat1[t]GAM1[t]+x2[t]PSat2[t]GAM2[t]+x3[t]PSat3[t]GAM3[t]-P)0,T[0]340,x1[0]0.05,x2[0]0.05},{x1,x2,T},{t,-80,80}];
In[]:=
sol5=Quiet@NDSolve[{D[x1[t],{t,1}]y1[t]-x1[t],D[x2[t],{t,1}]y2[t]-x2[t],(x1[t]PSat1[t]GAM1[t]+x2[t]PSat2[t]GAM2[t]+x3[t]PSat3[t]GAM3[t]-P)0,T[0]340,x1[0]0.75,x2[0]0.15},{x1,x2,T},{t,-80,80}];
In[]:=
sol6=Quiet@NDSolve[{D[x1[t],{t,1}]y1[t]-x1[t],D[x2[t],{t,1}]y2[t]-x2[t],(x1[t]PSat1[t]GAM1[t]+x2[t]PSat2[t]GAM2[t]+x3[t]PSat3[t]GAM3[t]-P)0,T[0]340,x1[0]0.65,x2[0]0.25},{x1,x2,T},{t,-80,80}];
In[]:=
sol7=Quiet@NDSolve[{D[x1[t],{t,1}]y1[t]-x1[t],D[x2[t],{t,1}]y2[t]-x2[t],(x1[t]PSat1[t]GAM1[t]+x2[t]PSat2[t]GAM2[t]+x3[t]PSat3[t]GAM3[t]-P)0,T[0]340,x1[0]0.5,x2[0]0.05},{x1,x2,T},{t,-80,80}];
In[]:=
sol8=Quiet@NDSolve[{D[x1[t],{t,1}]y1[t]-x1[t],D[x2[t],{t,1}]y2[t]-x2[t],(x1[t]PSat1[t]GAM1[t]+x2[t]PSat2[t]GAM2[t]+x3[t]PSat3[t]GAM3[t]-P)0,T[0]340,x1[0]0.35,x2[0]0.45},{x1,x2,T},{t,-80,80}];
Plot of Residue Curves
Azeotropes are show as green dots
In[]:=
DATA1=Flatten[Table[{x1[t],x2[t],x3[t]}/.sol1,{t,-80,80,0.4}],1];
In[]:=
DATA2=Flatten[Table[{x1[t],x2[t],x3[t]}/.sol2,{t,-80,80,0.4}],1];