The liquid–liquid equilibrium (LLE) behavior of the ternary system consisting of isobutyl acetate, acetic acid, and water at 283.15 K was modeled using the Non-Random Two-Liquid (NRTL) activity coefficient model. The binary interaction parameters and non-randomness factors reported by Cháfer et al. were used to calculate activity coefficients and the excess Gibbs free energy of mixing. The dimensionless Gibbs energy surface, (g/RT), was formulated as a function of composition and served as the thermodynamic basis for determining phase stability and phase splitting. The NRTL model captures the strong non-ideal interactions present in partially miscible liquid systems and is therefore well suited for describing liquid–liquid equilibria in this ternary mixture.
​
The equilibrium compositions of the two coexisting liquid phases were determined using the common tangent plane criterion. This approach requires that the Gibbs free energy surface possess a common tangent plane touching the surface at two distinct compositions corresponding to the equilibrium phases. Mathematically, the method involves solving a set of nonlinear equations that enforce equality of the tangent plane slopes and satisfy the Gibbs energy balance between the two phase compositions. An additional distance constraint was introduced to prevent convergence toward the trivial solution where both phases have identical compositions. The resulting system of equations defines a point on the binodal curve and provides the basis for tracing the complete liquid–liquid equilibrium boundary.
​
To generate the entire coexistence curve, a continuation technique based on arc-length parameterization was employed. Starting from an initial equilibrium solution obtained by FindRoot, the common tangent plane equations were transformed into a differential-algebraic system and solved using NDSolve. The continuation variable enabled smooth progression along the equilibrium manifold, allowing the binodal curve to be traced from one end of the two-phase region to the other. The method simultaneously produced the compositions of both equilibrium phases and identified the plait point, corresponding to the critical composition where the two liquid phases become indistinguishable.
​
The last Figure presents the calculated ternary phase diagram together with experimental data reported by Cháfer et al. The blue and yellow curve represents the predicted binodal boundary separating the single-phase and two-phase regions, while the red point indicates the calculated plait point. Red tie-lines connect the compositions of the two coexisting phases predicted by the NRTL model, whereas the blue dashed tie-lines correspond to experimental measurements reported in Ref. 1. The overall agreement between calculated and experimental results demonstrates the ability of the NRTL model combined with the common tangent plane continuation method to accurately reproduce the liquid–liquid equilibrium behavior of the isobutyl acetate–acetic acid–water system. Small deviations observed near the plait point are expected due to the increased sensitivity of phase compositions in the critical region.
Method based on common tangent plane and arc-length continuation
In[]:=
ClearAll["Global`*"];
NRTL Model-Cháfer et al.,Fluid Phase Equilibria 271 (2008)
Parameters at T=283.15 K,Table 4
R=8.314;(*J/(molK)*)​​T=283.15;(*K*)
alpha matrix (symmetric)
α={{0,0.2,0.2},{0.2,0,0.3},{0.2,0.3,0}};
tau matrix:tau[[i,j]]=Aij (energy parameter,J/mol),A1[[i,i]]=0
τ={{0,-3215.26,2744.08},{1829.92,0,-3695.34},{14809.44,6050.47,0}}/R/T;
In[]:=
α//MatrixForm
Out[]//MatrixForm=
0
0.2
0.2
0.2
0
0.3
0.2
0.3
0
In[]:=
τ//MatrixForm
Out[]//MatrixForm=
0.
-1.36581
1.16566
0.77733
0.
-1.56974
6.29089
2.57017
0.
In[]:=
G[i_,j_]:=Exp[-α[[i,j]]τ[[i,j]]];
In[]:=
TableForm@Table[G[i,j],{i,1,3},{j,1,3}]
Out[]//TableForm=
1.
1.31411
0.79205
0.856016
1.
1.60147
0.284171
0.462526
1.
Composition vector
xvec[x1_,x2_]:={x1,x2,1-x1-x2};
NRTL activity coefficients
lnGamma[i_,x1_,x2_]:=Module[{x,term1,term2},x=xvec[x1,x2];​​term1=Sum[x[[j]]τ[[j,i]]G[j,i],{j,1,3}]/Sum[x[[k]]G[k,i],{k,1,3}];​​term2=Sum[x[[j]]G[i,j]/Sum[x[[k]]G[k,j],{k,1,3}](τ[[i,j]]-Sum[x[[m]]τ[[m,j]]G[m,j],{m,1,3}]/Sum[x[[m]]G[m,j],{m,1,3}]),{j,1,3}];​​term1+term2];​​​​gamma[i_,x1_,x2_]:=Exp[lnGamma[i,x1,x2]];
Dimensionless Gibbs energy of mixing: g/RT
gmix[x1_,x2_]:=Module[{x},x=xvec[x1,x2];​​Sum[x[[i]]Log[x[[i]]],{i,1,3}]+Sum[x[[i]]lnGamma[i,x1,x2],{i,1,3}]];
Derivatives
dgdx1[x1_,x2_]:=Evaluate[D[gmix[u,v],u]/.{u->x1,v->x2}];​​dgdx2[x1_,x2_]:=Evaluate[D[gmix[u,v],v]/.{u->x1,v->x2}];
Unknowns
vars={x1a,x2a,x1b,x2b,m1,m2};
Common tangent equations
eq1=dgdx1[x1a,x2a]==m1;​​eq2=dgdx1[x1b,x2b]==m1;​​eq3=dgdx2[x1a,x2a]==m2;​​eq4=dgdx2[x1b,x2b]==m2;​​​​eq5=gmix[x1b,x2b]-gmix[x1a,x2a]==m1(x1b-x1a)+m2(x2b-x2a);
Extra constraint - prevent trivial solution xα=xβ
eq6=(x1b-x1a)^2+(x2b-x2a)^2==d;
In[]:=
(x1b-x1a)^2+(x2b-x2a)^2/.x1b->0.04/.x1a->0.2/.x2b->0.22/.x2a->0.3
Out[]=
0.032
First value for d --- very small d=0.001 yield approximately the plait point
d0t=0.8;
Solution
solct=FindRoot[{eq1,eq2,eq3,eq4,eq5,eq6/.d->d0t},{{x1a,0.925},{x2a,0.004},{x1b,0.0004},{x2b,0.0004},{m1,-10.0},{m2,40}},MaxIterations->10000]
Out[]=
{x1a0.895148,x2a0.0255223,x1b0.000881603,x2b0.00858798,m1-0.0788415,m2-4.74014}
In[]:=
init={x1a[0]==(x1a/.solct),x2a[0]==(x2a/.solct),x1b[0]==(x1b/.solct),x2b[0]==(x2b/.solct),m1[0]==(m1/.solct),m2[0]==(m2/.solct),d[0]==d0t}
Out[]=
{x1a[0]0.895148,x2a[0]0.0255223,x1b[0]0.000881603,x2b[0]0.00858798,m1[0]-0.0788415,m2[0]-4.74014,d[0]0.8}
In[]:=
initdev={x1a'[0]==-0.0001,x2a'[0]==0.0001,x1b'[0]==0.0001,x2b'[0]==0.0001,m1'[0]==0.01,m2'[0]==0.01,d'[0]==-0.00001}
Out[]=
{
′
x1a
[0]-0.0001,
′
x2a
[0]0.0001,
′
x1b
[0]0.0001,
′
x2b
[0]0.0001,
′
m1
[0]0.01,
′
m2
[0]0.01,
′
d
[0]-0.00001}
In[]:=
SYS={eq1,eq2,eq3,eq4,eq5,eq6}/.x1a->x1a[s]/.x1b->x1b[s]/.x2a->x2a[s]/.x2b->x2b[s]/.m1->m1[s]/.m2->m2[s]/.d->d[s];
In[]:=
auxeq=x1a'[s]^2+x1b'[s]^2+x2a'[s]^2+x2b'[s]^2+m1'[s]^2+m2'[s]^2+d'[s]^2==1
Out[]=
2
′
d
[s]
+
2
′
m1
[s]
+
2
′
m2
[s]
+
2
′
x1a
[s]
+
2
′
x1b
[s]
+
2
′
x2a
[s]
+
2
′
x2b
[s]
1
In[]:=
SYSFinal=Join[SYS,{auxeq},init,initdev];
In[]:=
SOL=Quiet@NDSolve[SYSFinal,{x1a,x2a,x1b,x2b,m1,m2,d},{s,-10,10},AccuracyGoal->5];
In[]:=
DATA1=Table[{x1a[s],x2a[s],1-x1a[s]-x2a[s]}/.SOL,{s,-8.064,3.435,0.001}];​​DATA2=Table[{x1b[s],x2b[s],1-x1b[s]-x2b[s]}/.SOL,{s,-8.064,3.435,0.001}];
In[]:=
Pt={{2,1},{0,Sqrt[3.]}}/2;
Type I liquid-liquid equilibrium (LLE) diagram - Plait point indicated in red
In[]:=
p0=Show[TernaryListPlot[{Flatten[DATA1,1],Flatten[DATA2,1]},Joined->True,PlotStyle->Thick,GridLines->Automatic,FrameLabel->{"Isobutyl acetate","Acetic acid","Water"},PlotRange->All,Method->{"Backdrop"->LightOrange}],Graphics[{Red,PointSize[0.02],Point[Pt.{x1a[s],x2a[s]}/.SOL/.s->3.435]}]]
Out[]=
Data
Experimental LLE data-Table 2
​
Isobutyl acetate (1)+Acetic acid (2)+Water (3)
​
T=283.15 K
​
Columns: x1,x2 for each phase (x3=1-x1-x2)
Isobutyl acetate-rich phase: {x1,x2}
acetateRichPhase283={{0.945,0.000},{0.892,0.020},{0.847,0.054},{0.792,0.092},{0.744,0.120},{0.704,0.140},{0.664,0.159},{0.615,0.187},{0.572,0.206},{0.522,0.227},{0.483,0.248},{0.447,0.264},{0.412,0.275},{0.363,0.289},{0.320,0.295},{0.265,0.309},{0.205,0.314}};
Water-rich phase: {x1,x2}
waterRichPhase283={{0.001,0.000},{0.002,0.011},{0.001,0.023},{0.002,0.038},{0.002,0.049},{0.002,0.057},{0.003,0.068},{0.003,0.079},{0.004,0.089},{0.005,0.100},{0.005,0.116},{0.006,0.128},{0.008,0.139},{0.011,0.158},{0.015,0.173},{0.024,0.199},{0.040,0.227}};
Add x3 column (water/acetate,by difference) and full triples
acetateRichFull283=Map[Append[#,1-Total[#]]&,acetateRichPhase283];​​waterRichFull283=Map[Append[#,1-Total[#]]&,waterRichPhase283];
Display as tables with headers
TableForm[acetateRichFull283,TableHeadings->{None,{"x1 (acetate)","x2 (acetic acid)","x3 (water)"}}]​​​​TableForm[waterRichFull283,TableHeadings->{None,{"x1 (acetate)","x2 (acetic acid)","x3 (water)"}}]
Out[]//TableForm=
x1 (acetate)
x2 (acetic acid)
x3 (water)
0.945
0.
0.055
0.892
0.02
0.088
0.847
0.054
0.099
0.792
0.092
0.116
0.744
0.12
0.136
0.704
0.14
0.156
0.664
0.159
0.177
0.615
0.187
0.198
0.572
0.206
0.222
0.522
0.227
0.251
0.483
0.248
0.269
0.447
0.264
0.289
0.412
0.275
0.313
0.363
0.289
0.348
0.32
0.295
0.385
0.265
0.309
0.426
0.205
0.314
0.481
Out[]//TableForm=
x1 (acetate)
x2 (acetic acid)
x3 (water)
0.001
0.
0.999
0.002
0.011
0.987
0.001
0.023
0.976
0.002
0.038
0.96
0.002
0.049
0.949
0.002
0.057
0.941
0.003
0.068
0.929
0.003
0.079
0.918
0.004
0.089
0.907
0.005
0.1
0.895
0.005
0.116
0.879
0.006
0.128
0.866
0.008
0.139
0.853
0.011
0.158
0.831
0.015
0.173
0.812
0.024
0.199
0.777
0.04
0.227
0.733
Combine into tie-line pairs:{phaseI,phaseII}

References

A. Cháfer, E. Lladosa, J. de la Torre, M.C. Burguet, Study of liquid–liquid equilibrium of the systems isobutyl acetate+acetic acid+water and isobutyl alcohol+acetic acid+water at different temperatures, Fluid Phase Equilibria, Volume 271, Issues 1–2, 2008, Pages 76-81, https://doi.org/10.1016/j.fluid.2008.07.001

CITE THIS NOTEBOOK

Arc-length continuation of common tangent plane solutions in NRTL-based liquid–liquid equilibrium​
by Housam Binous​
Wolfram Community, STAFF PICKS, June 25, 2026
​https://community.wolfram.com/groups/-/m/t/3739509