Thermodynamic Modeling of Ternary LLE Systems : NRTL Parameter Fitting by Hybrid PSO
and Phase Diagram Construction via Continuation Methods
Housam Binous and Abdelmalek Hasseine
and Phase Diagram Construction via Continuation Methods
Housam Binous and Abdelmalek Hasseine
Part A The code begins by fixing the system temperature at 298.15 K. Experimental liquid–liquid equilibrium (LLE) tie-line data for the ternary system (n‑Octane (1) + Xylene (2) + Sulfolane (3)) are also defined. Then the notebook automatically determines the number of components (NC) and the number of experimental tie-lines (NTIL) directly from the dimensions of the dataset. The raffinate-phase and extract-phase mole fractions for each component are then extracted into separate variable arrays and packaged into the xxexp and yyexp matrices, which serve as the experimental reference data that the NRTL model will be fitted against.The core thermodynamic model is implemented in the NRTL module, which computes the activity of each component in a liquid mixture from the NRTL (Non-Random Two-Liquid) equation. Given a composition vector, the module builds the binary interaction parameter matrix tau from the temperature-dependent parameters P, computes the corresponding G matrix using a fixed non-randomness parameter (α = 0.2 for all binary pairs), and then evaluates the local composition terms needed to construct each component’s activity coefficient. The output is the vector of thermodynamic activities (mole fraction × activity coefficient) for every component in the mixture, which is the key quantity used to test whether two liquid phases are in equilibrium.Phase equilibrium is enforced through the EQUILIBRE module, which sets up the classic isoactivity conditions (equal component activities in both phases), together with mass-balance and phase-fraction constraints, as a system of nonlinear equations. For each experimental tie-line, MAIN calls FindRoot to numerically solve this system and predict the equilibrium raffinate and extract compositions consistent with a given set of trial NRTL parameters. The Funct module then assembles the full binary interaction parameter matrix P from a flat parameter vector, runs MAIN over all tie-lines, and returns the sum of squared deviations between the model-predicted and experimentally measured compositions — this sum-of-squares residual is the objective function that must be minimized to identify the best-fit NRTL parameters.Because the parameter-estimation problem is highly nonlinear and prone to multiple local minima, the fitting is carried out with a hybrid global–local strategy. Particle Swarm Optimization (psoOptimize) first explores the six-dimensional parameter space (the independent NRTL binary interaction parameters) broadly, using inertia weighting and cognitive/social acceleration coefficients to iteratively update particle positions and velocities toward improving regions of parameter space over a fixed number of iterations. The best particle position found by PSO (gbest) then serves as the initial guess for a local interior-point-style refinement based on FindMinimum, sharpening the fit down to a more precise minimum. The final output includes the optimized P matrix of NRTL binary interaction parameters, the minimized objective function value (finalLoss), and the root-mean-square deviation (rsmdFinal) between predicted and experimental compositions, providing a quantitative measure of how well the fitted NRTL model reproduces the observed ternary LLE behavior.Part B Then, we compute the liquid - liquid equilibrium (LLE) for this ternary system n‑Octane (1) + Xylene (2) + Sulfolane (3) at 298.15 K using the common tangent plane / arc - length continuation method . It defines the NRTL model (τ_ij matrix, α = 0.2 nonrandomness parameters, G_ij, and the dimensionless Gibbs energy of mixing gmix), then derives the partial derivatives dgdx1, dgdx2 symbolically using D operator and Evaluate . The core idea is to find two compositions (x1a, x2a) and (x1b, x2b) on the Gibbs - energy surface that share a common tangent line — the classic equal - chemical - potential condition for phase equilibrium — solved via a system of common - tangent equations (eqs . 1–4) plus normalization (eq5) and an arc - length parameter d (eq6) that' s swept from a small initial value to trace out the full binodal curve . Starting from an initial solution (e . g ., for d = 0.0001) found by FindRoot, the notebook uses NDSolve to continue this solution as a function of the arc - length parameter d (s), producing interpolating functions for x1a (s), x2a (s), x1b (s), x2b (s) and the slopes m1 (s), m2 (s) . This traces the entire two - phase boundary (binodal) in one continuous sweep rather than requiring separate root - finds for each tie line . The plait point (critical point where the two phases merge) is obtained for d (s) = 0 and s = -0.00486519. The continuous binodal curve is plotted on a ternary diagram, then overlaid with the 9 experimental tie lines (shown in dashed blue) from Table 2 of Chen et al . (2000) for direct comparison . For each experimental tie line, the notebook computes the arc - length, s, (see mytielinesarclength list) and calculate the corresponding tie line (shown in red) . The RMSD between calculated and experimental mole fractions across all tie lines and components is computed; RMSD ≈ 0.00506625 confirming good agreement between the continuation - method predictions and the published experimental LLE data .Reference:[1] J . Chen, L . -P . Duan, J . -G . Mi, W . -Y . Fei, Z . -C . Li, "Liquid–liquid equilibria of multi-component systems including n-hexane, n-octane, benzene, toluene, xylene and sulfolane at 298.15 K and atmospheric pressure," Fluid Phase Equilibria, Vol . 173, 2000, pp . 109–119. DOI/PII : S0378 - 3812 (00) 00398 - 8
Part A: NRTL model determination using Hybrid method
Part A: NRTL model determination using Hybrid method
Determination of NRTL parameters is based on a Hybrid PSO - InteriorPoint approach
In[]:=
ClearAll[NC,T,NTIL,xxexp,yyexp,P,z];T=298.15;data={{0.831,0.165,0.005,0.007,0.048,0.944},{0.752,0.242,0.006,0.007,0.063,0.93},{0.677,0.318,0.006,0.008,0.087,0.904},{0.609,0.381,0.01,0.009,0.106,0.885},{0.547,0.436,0.017,0.009,0.125,0.866},{0.486,0.494,0.02,0.01,0.145,0.845},{0.429,0.541,0.03,0.013,0.199,0.788},{0.312,0.64,0.048,0.014,0.231,0.754},{0.262,0.672,0.066,0.015,0.267,0.717}};{m,n}=Dimensions[data];NC=n/2;NTIL=m;x1r=data[[All,1]];x2r=data[[All,2]];x3r=data[[All,3]];x1e=data[[All,4]];x2e=data[[All,5]];x3e=data[[All,6]];(*PHASERAFFINATE/PHASEEXTRACT,n-Octane(1)+Xylene(2)+Sulfolane(3)*)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
Part B: LLE diagram computation using common tangent plane and arc - length continuation
Method based on common tangent plane and arc - length continuation
NRTL Model-Chen et al., Fluid Phase Equilibria 173 (2000)
Parameters at T=298.15 K from my above hybrid optimization code
Parameters at T=298.15 K from my above hybrid optimization code
tau_ij for n-Octane(1)+Xylene(2)+Sulfolane(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 --- very small d=0.0001 yield approximately the plait point
First value for d --- very small d=0.0001 yield approximately the plait point
Solution is near {0.033,0.61}
Type I liquid-liquid equilibrium (LLE) diagram - Plait point indicated in red
Octane-rich phase {x1,x2,x3}, Sulfolane-rich phase {x1,x2,x3}
Octane-rich phase {x1,x2}, Sulfolane-rich phase {x1,x2}
Octane-rich phase {x1,x2}, Sulfolane-rich phase {x1,x2}
Add x3 column (water/acetate,by difference) and full triples
Display as tables with headers
Combine into tie-line pairs:{phaseI,phaseII}
Show all nine experimental tie lines
Let us compute the RMSD