Liquid-Liquid Phase Split Calculations for the n-Octane + Toluene + Xylene + Sulfolane System
Housam Binous
Liquid-Liquid Phase Split Calculations for the n-Octane + Toluene + Xylene + Sulfolane System
Housam Binous
Housam Binous
The notebook implements liquid–liquid phase-split (flash) calculations for the quaternary system n‑Octane(1) + Toluene(2) + Xylene(3) + Sulfolane(4) at 298.15 K, using the NRTL (non-random two-liquid) activity coefficient model with binary interaction parameters τ_ij taken from Table 4 of Chen et al. The raw τ_ij values (in K) are divided by T to render them dimensionless and the non-randomness parameters are all equal to 0.2. The function lnGamma[i, x, T] implements the full multicomponent NRTL expression (Ref.1 Eq. 2) to compute the natural log of the activity coefficient of component i in a mixture of composition x, and gamma[i, x, T] exponentiates it to give the activity coefficient directly.
Before solving for a tie line, the notebook defines a TPD[w1, w2, w3] function that evaluates the tangent-plane-distance criterion (Michelsen-style stability analysis) for a trial composition w relative to the overall feed composition z. This function checks physical feasibility (all mole fractions must be positive and sum to 1, returning a large penalty value of 10^10 otherwise) and computes the weighted sum of differences in chemical potential (via log of activity × composition) between the trial phase and the feed. NMinimize with the DifferentialEvolution method is used to search for a global minimum of TPD over the trial composition space; a negative TPDmin confirms that the single-phase feed is thermodynamically unstable and will split into two liquid phases, while also yielding an initial guess (w1, w2, w3) for one of the two phase compositions.
The core eqns[xa1, xa2, xa3, xb1, xb2, xb3, beta] function sets up the full isoactivity (equal fugacity) and Rachford–Rice-style mass-balance system for a two-liquid-phase flash: it computes activity coefficients for trial compositions xa (phase I) and xb (phase II), forms the K-values K = γ_a/γ_b, and assembles seven equations — three isoactivity conditions (x_a·γ_a = x_b·γ_b for components 1–3), three mass-balance relations linking each phase composition back to the overall feed z via the phase-split fraction β, and a Rachford–Rice-type summation residual: rr. For each tie line, the TPD-derived trial composition is used to seed xaGuess, and xbGuess is backed out algebraically from the mass balance; FindRoot is then called (with MaxIterations → 1000, AccuracyGoal/PrecisionGoal → 10) to solve the nonlinear system simultaneously for both phase compositions and β.
Validation Against Experimental Tie Lines: the notebook repeats this full TPD → guess-generation → FindRoot workflow for four separate experimental tie lines (Tie line 1 through 4) drawn from the paper’s Table 5 for this quaternary system. For each tie line, an overall feed z is constructed as the midpoint of the two experimental phase compositions (xPhaseI, xPhaseII), and the solver’s converged result is cross-checked using FindMinimum on the sum of squared residuals of eqns — in every case returning a near-zero objective value (e.g., 10^−21 to 10^−31), confirming that the NRTL-predicted phase split essentially exactly reproduces the algebraic flash equations once a good initial guess is supplied, validating both the τ_ij parameterization and the solver’s convergence to the correct two-phase equilibrium.
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
The notebook implements liquid–liquid phase-split (flash) calculations for the quaternary system n‑Octane(1) + Toluene(2) + Xylene(3) + Sulfolane(4) at 298.15 K, using the NRTL (non-random two-liquid) activity coefficient model with binary interaction parameters τ_ij taken from Table 4 of Chen et al. The raw τ_ij values (in K) are divided by T to render them dimensionless and the non-randomness parameters are all equal to 0.2. The function lnGamma[i, x, T] implements the full multicomponent NRTL expression (Ref.1 Eq. 2) to compute the natural log of the activity coefficient of component i in a mixture of composition x, and gamma[i, x, T] exponentiates it to give the activity coefficient directly.
Before solving for a tie line, the notebook defines a TPD[w1, w2, w3] function that evaluates the tangent-plane-distance criterion (Michelsen-style stability analysis) for a trial composition w relative to the overall feed composition z. This function checks physical feasibility (all mole fractions must be positive and sum to 1, returning a large penalty value of 10^10 otherwise) and computes the weighted sum of differences in chemical potential (via log of activity × composition) between the trial phase and the feed. NMinimize with the DifferentialEvolution method is used to search for a global minimum of TPD over the trial composition space; a negative TPDmin confirms that the single-phase feed is thermodynamically unstable and will split into two liquid phases, while also yielding an initial guess (w1, w2, w3) for one of the two phase compositions.
The core eqns[xa1, xa2, xa3, xb1, xb2, xb3, beta] function sets up the full isoactivity (equal fugacity) and Rachford–Rice-style mass-balance system for a two-liquid-phase flash: it computes activity coefficients for trial compositions xa (phase I) and xb (phase II), forms the K-values K = γ_a/γ_b, and assembles seven equations — three isoactivity conditions (x_a·γ_a = x_b·γ_b for components 1–3), three mass-balance relations linking each phase composition back to the overall feed z via the phase-split fraction β, and a Rachford–Rice-type summation residual: rr. For each tie line, the TPD-derived trial composition is used to seed xaGuess, and xbGuess is backed out algebraically from the mass balance; FindRoot is then called (with MaxIterations → 1000, AccuracyGoal/PrecisionGoal → 10) to solve the nonlinear system simultaneously for both phase compositions and β.
Validation Against Experimental Tie Lines: the notebook repeats this full TPD → guess-generation → FindRoot workflow for four separate experimental tie lines (Tie line 1 through 4) drawn from the paper’s Table 5 for this quaternary system. For each tie line, an overall feed z is constructed as the midpoint of the two experimental phase compositions (xPhaseI, xPhaseII), and the solver’s converged result is cross-checked using FindMinimum on the sum of squared residuals of eqns — in every case returning a near-zero objective value (e.g., 10^−21 to 10^−31), confirming that the NRTL-predicted phase split essentially exactly reproduces the algebraic flash equations once a good initial guess is supplied, validating both the τ_ij parameterization and the solver’s convergence to the correct two-phase equilibrium.
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
T=298.15;
In[]:=
tau={{0.0000,0.1155*10^3,-0.3318*10^3,0.2235*10^4},{-0.8422*10^2,0.0000,0.0000,0.5623*10^3},{0.3645*10^3,0.0000,0.0000,0.4137*10^3},{0.1231*10^4,-0.5359*10^2,0.1451*10^3,0.0000}}/T;
In[]:=
tau//TableForm
Out[]//TableForm=
0. | 0.387389 | -1.11286 | 7.49623 |
-0.282475 | 0. | 0. | 1.88596 |
1.22254 | 0. | 0. | 1.38756 |
4.12879 | -0.179742 | 0.486668 | 0. |
In[]:=
(*---G[i,j]asafunctionoftemperatureT---*)G[i_,j_,T_]:=Exp[-0.2*tau[[i,j]]](*tau[i,i]=0andG[i,i]=1byconvention*)G[i_,i_,T_]:=1(*NRTLactivitycoefficient(Eq.2inthepaper)*)lnGamma[i_,x_List,T_]:=Module[{n=Length[x]},Sum[tau[[j,i]]*G[j,i,T]*x[[j]],{j,n}]/Sum[G[k,i,T]*x[[k]],{k,n}]+Sum[(G[i,j,T]*x[[j]])/Sum[G[k,j,T]*x[[k]],{k,n}]*(tau[[i,j]]-Sum[tau[[l,j]]*G[l,j,T]*x[[l]],{l,n}]/Sum[G[k,j,T]*x[[k]],{k,n}]),{j,n}]]gamma[i_,x_List,T_]:=Exp[lnGamma[i,x,T]]
In[]:=
TPD[w1_?NumericQ,w2_?NumericQ,w3_?NumericQ]:=Module[{w4,w,gammaW},w4=1-w1-w2-w3;If[w1<=0||w2<=0||w3<=0||w4<=0,Return[10^10]];w={w1,w2,w3,w4};gammaW=Table[gamma[i,w,T],{i,1,4}];Sum[w[[i]]*(Log[w[[i]]*gammaW[[i]]]-Log[z[[i]]*gamma[i,z,T]]),{i,1,4}]];
In[]:=
eqns[xa1_,xa2_,xa3_,xb1_,xb2_,xb3_,beta_]:=Module[{xa4,xb4,xa,xb,ga,gb,K,rr},xa4=1-xa1-xa2-xa3;xb4=1-xb1-xb2-xb3;xa={xa1,xa2,xa3,xa4};xb={xb1,xb2,xb3,xb4};ga=Table[gamma[i,xa,T],{i,1,4}];gb=Table[gamma[i,xb,T],{i,1,4}];K=ga/gb;rr=Sum[z[[i]]*(K[[i]]-1)/(1+beta*(K[[i]]-1)),{i,1,4}];{xa[[1]]*ga[[1]]-xb[[1]]*gb[[1]],xa[[2]]*ga[[2]]-xb[[2]]*gb[[2]],xa[[3]]*ga[[3]]-xb[[3]]*gb[[3]],xa[[1]]-z[[1]]/(1+beta*(K[[1]]-1)),xa[[2]]-z[[2]]/(1+beta*(K[[2]]-1)),xa[[3]]-z[[3]]/(1+beta*(K[[3]]-1)),rr}];
Tie line 1
In[]:=
xPhaseI={0.251,0.171,0.507,0.071};xPhaseII={0.019,0.091,0.209,0.681};z=(xPhaseI+xPhaseII)/2
Out[]=
{0.135,0.131,0.358,0.376}
In[]:=
sol=NMinimize[{TPD[w1,w2,w3],w1>=10^-6,w2>=10^-6,w3>=10^-6,1-w1-w2-w3>=10^-6},{w1,w2,w3},Method->{"DifferentialEvolution"}];TPDmin=sol[[1]];Print[" TPDmin = ",NumberForm[TPDmin,{12,8}]]
TPDmin = -0.06848818
In[]:=
sol
Out[]=
{-0.0684882,{w10.00853591,w20.0571731,w30.1315}}
In[]:=
betaGuess=0.5;
In[]:=
xbGuess={w1,w2,w3}/.sol[[2]];(*fromTPD*)(*backoutxafrommassbalanceusingnormalizedz*)xaGuess=(z[[1;;3]]-(1-betaGuess)*xbGuess)/betaGuess;
In[]:=
solRoot=FindRoot[eqns[xa1,xa2,xa3,xb1,xb2,xb3,beta],{{xa1,xaGuess[[1]]},{xa2,xaGuess[[2]]},{xa3,xaGuess[[3]]},{xb1,xbGuess[[1]]},{xb2,xbGuess[[2]]},{xb3,xbGuess[[3]]},{beta,betaGuess}},MaxIterations->1000,AccuracyGoal->10,PrecisionGoal->10]
Out[]=
{xa10.255554,xa20.169616,xa30.504078,xb10.0114728,xb20.0914312,xb30.208319,beta0.493909}
In[]:=
FindMinimum[Sum[eqns[xa1,xa2,xa3,xb1,xb2,xb3,beta][[i]]^2,{i,1,7}],{{xa1,xaGuess[[1]]},{xa2,xaGuess[[2]]},{xa3,xaGuess[[3]]},{xb1,xbGuess[[1]]},{xb2,xbGuess[[2]]},{xb3,xbGuess[[3]]},{beta,betaGuess}},MaxIterations->1000,AccuracyGoal->10,PrecisionGoal->10]
Out[]=
{7.6019×,{xa10.255554,xa20.169616,xa30.504078,xb10.0114728,xb20.0914312,xb30.208319,beta0.493909}}
-23
10
Tie line 2
In[]:=
xPhaseI={0.649,0.154,0.184,0.013};xPhaseII={0.007,0.058,0.041,0.894};z=(xPhaseI+xPhaseII)/2
Out[]=
{0.328,0.106,0.1125,0.4535}
In[]:=
sol=NMinimize[{TPD[w1,w2,w3],w1>=10^-6,w2>=10^-6,w3>=10^-6,1-w1-w2-w3>=10^-6},{w1,w2,w3},Method->{"DifferentialEvolution"}];TPDmin=sol[[1]];Print[" TPDmin = ",NumberForm[TPDmin,{12,8}]]
TPDmin = -0.53066654
In[]:=
sol
Out[]=
{-0.530667,{w10.891235,w20.0410038,w30.067337}}
In[]:=
betaGuess=0.5;
In[]:=
xaGuess={w1,w2,w3}/.sol[[2]];(*fromTPD*)(*backoutxafrommassbalanceusingnormalizedz*)xbGuess=(z[[1;;3]]-(1-betaGuess)*xaGuess)/betaGuess;
In[]:=
solRoot=FindRoot[eqns[xa1,xa2,xa3,xb1,xb2,xb3,beta],{{xa1,xaGuess[[1]]},{xa2,xaGuess[[2]]},{xa3,xaGuess[[3]]},{xb1,xbGuess[[1]]},{xb2,xbGuess[[2]]},{xb3,xbGuess[[3]]},{beta,betaGuess}},MaxIterations->1000,AccuracyGoal->10,PrecisionGoal->10]
Out[]=
{xa10.660708,xa20.150867,xa30.182495,xb10.00617711,xb20.0626012,xb30.0447949,beta0.508315}
In[]:=
FindMinimum[Sum[eqns[xa1,xa2,xa3,xb1,xb2,xb3,beta][[i]]^2,{i,1,7}],{{xa1,xaGuess[[1]]},{xa2,xaGuess[[2]]},{xa3,xaGuess[[3]]},{xb1,xbGuess[[1]]},{xb2,xbGuess[[2]]},{xb3,xbGuess[[3]]},{beta,betaGuess}},MaxIterations->1000,AccuracyGoal->10,PrecisionGoal->10]
Out[]=
{6.30233×,{xa10.660708,xa20.150867,xa30.182495,xb10.00617711,xb20.0626012,xb30.0447949,beta0.508315}}
-24
10
Tie line 3
In[]:=
xPhaseI={0.548,0.244,0.187,0.022};xPhaseII={0.011,0.105,0.057,0.827};z=(xPhaseI+xPhaseII)/2
Out[]=
{0.2795,0.1745,0.122,0.4245}
In[]:=
sol=NMinimize[{TPD[w1,w2,w3],w1>=10^-6,w2>=10^-6,w3>=10^-6,1-w1-w2-w3>=10^-6},{w1,w2,w3},Method->{"DifferentialEvolution"}];TPDmin=sol[[1]];Print[" TPDmin = ",NumberForm[TPDmin,{12,8}]]
Tie line 4