Thermodynamic Modeling of Ternary LLE (n-Hexane–Methylcyclopentane–Aniline at 25°C and 1 atm)
​
​Housam Binous and Abdelmalek Hasseine
​The ternary system n-hexane–methylcyclopentane–aniline has long served as a benchmark for liquid-liquid extraction studies, since n-hexane (a paraffin) and methylcyclopentane (a naphthene) are notoriously difficult to separate by distillation due to their close boiling points, while aniline acts as a selective solvent that preferentially dissolves the naphthene. The experimental dataset comes from Darwent and Winkler (1943), who measured the mutual solubility and tie-line data for this system at 25.0°C. This system is Type II (both n-hexane–aniline and methylcyclopentane–aniline exhibit partial miscibility, producing two separate solubility lobes). It should be reminded that this ternary LLE systems exhibits a transition between Type I and Type II phase behavior at 34.5°C.​Building on the 1943 experimental data, Binous and Bellagi (2021), in “Calculation of Ternary Liquid-Liquid Equilibrium Data Using Arc-Length Continuation” (Engineering Reports, 3(2), e12296), developed a numerical framework to compute the binodal curve and tie lines for this system using modern activity-coefficient models — specifically NRTL and UNIQUAC — with parameters estimated via particle swarm optimization (PSO). Their explicit goal was to reproduce the Darwent and Winkler experimental solubility and tie-line data as a validation case, alongside six other ternary systems spanning both Type I and Type II behavior.​​Part A — NRTL parameter determination via the Hybrid method​​The first part of the uploaded notebook implements NRTL model determination using a “Hybrid method” consisting of a combination of global search methaheurisitc algorithm (PSO) and local refinement (interior-point of FindMinimum).​​Part B — LLE diagram computation via common tangent plane and arc-length continuation.​​The second part of the notebook implements a new numerical method: given the fitted NRTL parameters from Part A, it computes the full ternary LLE diagram (binodal curve plus tie lines) using the common tangent plane criterion for phase equilibrium, combined with arc-length continuation to trace the solution branch robustly through the curved, sometimes near-singular regions of the phase diagram (particularly near plait points, where standard Newton-Raphson-type solvers can fail to converge).​The Mathematica code sets up the isoactivity/equal-fugacity equations for each component across both phases, parametrizes the solution path by arc length rather than composition directly (avoiding the ill-conditioning that occurs near the plait point where dx/ds → ∞ in ordinary parametrizations), and interpolates the resulting discrete solution points into continuous functions for the binodal curve and tie-line endpoints. We conclude by overlaying the calculated curve against the digitized experimental tie-line data from Darwent and Winkler for direct visual and numerical comparison, essentially reproducing Figure 12 of the Binous and Bellagi paper.​Together, Parts A and B represent a methodology applied to the exact experimental system reported by Darwent and Winkler (1943): first regressing thermodynamically consistent NRTL parameters from binary solubility data, then using those parameters to predict the full ternary phase diagram via a numerically robust continuation method, and finally validating the prediction against the experimental tie-line data. ​​References:​[1] Darwent, B. de B., and Winkler, C. A. “The System n-Hexane–Methylcyclopentane–Aniline.” The Journal of Physical Chemistry, vol. 47, no. 6, 1943, pp. 442–454. https://pubs.acs.org/doi/abs/10.1021/j150429a005[2] Binous, Housam, and Ahmed Bellagi. “Calculation of Ternary Liquid-Liquid Equilibrium Data Using Arc-Length Continuation.” Engineering Reports, vol. 3, no. 2, 2021, e12296. https://doi.org/10.1002/eng2.12296

Part A: NRTL model determination using Hybrid method

Determination of NRTL parameters is based on a Hybrid PSO - InteriorPoint approach
tie lines adapted from data of Ref. 1 (data in mole %) Aniline-Rich phase Hydrocarbon-Rich phase​MCP Aniline Hexane MCP Aniline Hexane​0. 0.921 0.079 0. 0.067 0.933​0.01 0.916 0.074 0.109 0.071 0.82​0.008 0.916 0.076 0.23 0.079 0.691​0.035 0.894 0.071 0.3 0.081 0.619​0.065 0.878 0.057 0.479 0.1 0.421​0.131 0.834 0.035 0.698 0.111 0.191​0.191 0.792 0.017 0.796 0.132 0.072​0.242 0.758 0. 0.857 0.143 0.​​​
In[]:=
ClearAll[NC,T,NTIL,xxexp,yyexp,P,z];​​​​T=298.15;​​​​data={{0.,0.921,0.079,0.,0.067,0.933},{0.01,0.916,0.074,0.109,0.071,0.82},{0.008,0.916,0.076,0.23,0.079,0.691},{0.035,0.894,0.071,0.3,0.081,0.619},{0.065,0.878,0.057,0.479,0.1,0.421},{0.131,0.834,0.035,0.698,0.111,0.191},{0.191,0.792,0.017,0.796,0.132,0.072},{0.242,0.758,0.,0.857,0.143,2.77556*10^-17}};
In[]:=
{m,n}=Dimensions[data];​​​​NC=n/2;​​​​NTIL=m;​​​​x1e=data[[All,1]];x2e=data[[All,2]];x3e=data[[All,3]];​​x1r=data[[All,4]];x2r=data[[All,5]];x3r=data[[All,6]];​​​​(*PHASERAFFINATE/PHASEEXTRACTforMCP(1)+Aniline(2)+Hexane(3)ternarysystem*)​​​​xxexp={x1r,x2r,x3r};​​yyexp={x1e,x2e,x3e};​​
In[]:=
NRTL[x_]:=Module[{R=1.987,alp,tau,G,gam,G1,tau1,G2,s,act,aa,bb,zz,i,j},​​alp=Table[0.,{NC},{NC}];​​alp[[1,2]]=0.2;alp[[1,3]]=0.2;alp[[2,3]]=0.2;​​alp[[2,1]]=alp[[1,2]];alp[[3,1]]=alp[[1,3]];alp[[3,2]]=alp[[2,3]];​​tau=Table[P[[i,j]]/T,{i,NC},{j,NC}];​​G=Table[Exp[-alp[[i,j]]tau[[i,j]]],{i,NC},{j,NC}];​​gam=Table[0.,{NC}];​​G1=Table[0.,{NC},{NC}];​​tau1=Table[0.,{NC},{NC}];​​G2=Table[0.,{NC},{NC}];​​s=Table[0.,{NC},{NC}];​​Do[​​aa=0.;bb=0.;​​Do[​​zz=G[[j,i]]x[[j]];​​aa+=zz;​​bb+=zztau[[j,i]];​​,{j,NC}];​​gam[[i]]=bb/aa;​​Do[​​G1[[j,i]]=G[[j,i]]/aa;​​tau1[[j,i]]=tau[[j,i]]-gam[[i]];​​G2[[j,i]]=G1[[j,i]]tau1[[j,i]];​​s[[j,i]]=x[[i]]G2[[j,i]];​​,{j,NC}];​​,{i,NC}];​​Do[​​Do[gam[[i]]+=s[[i,j]],{j,NC}];​​gam[[i]]=Exp[gam[[i]]];​​,{i,NC}];​​act=Table[x[[i]]gam[[i]],{i,NC}];​​act]
In[]:=
​​EQUILIBRE[X_]:=Module[{x,y,act1,act2,f,Xt,Yt,i},​​x=Take[X,NC];​​y=Take[X,{NC+1,2NC}];​​​​act1=NRTL[x];act2=NRTL[y];​​​​f=Table[0.,{2NC+1}];​​Do[f[[i]]=act1[[i]]-act2[[i]],{i,NC}];​​Xt=Sum[X[[i]],{i,NC}];​​Yt=Sum[X[[i]],{i,NC+1,2NC}];​​f[[NC+1]]=1-Xt;​​f[[NC+2]]=1-Yt;​​Do[​​f[[NC+2+i]]=z[[i]]-X[[i]]X[[2NC+1]]-X[[i+NC]](1-X[[2NC+1]]),​​{i,NC-1}];​​f]

Part B: LLE diagram computation using common tangent plane and arc - length continuation

Method based on common tangent plane and arc - length continuation
Parameters at T=298.15 K from my above hybrid optimization code
tau_ij for MCP(1) + Aniline(2) + Hexane (3)
Alpha matrix (symmetric)
Composition vector
NRTL activity coefficients
Dimensionless Gibbs energy of mixing: g/RT
Derivatives
Unknowns
Common tangent equations
Extra constraint - Prevent trivial solution xα=xβ
First value for d=0.73
Solution is near {0.035,0.88,0.30,0.08}
aniline-rich phase {x1,x2,x3}, hydrocarbon-rich phase {x1,x2,x3}
Display as tables with headers
Combine into tie-line pairs:{phaseI,phaseII}
Show all nine experimental tie lines
Let us compute the RMSD