Statistical Thermodynamic Models: Non-random Two-liquid (NRTL) models

Currently the Cebule engine supports the following models:

• COSMO-RS & COSMO-SAC (Link)
• Non-random Two-liquid (NRTL, eNRTL), Wilson model and COSMO-NRTL (this page)
• Group contribution (GC) models including UNIQUAC, Original UNIFAC, UNIFAC (Dortmund), UNIFAC (Lyngby), UNIFAC 2.0 (Kaiserslautern) and Predictive Soave-Redlich Kwong (PSRK) (Link)
• Hybrid models: UNIFAC-VISCO, UNIFAC-IL (Link)

Non-random Two-liquid (NRTL) models

Local composition models represent a class of thermodynamic models that account for the non-uniform distribution of molecules in liquid mixtures. Unlike ideal solution theory which assumes random molecular mixing, these models recognize that molecular interactions lead to preferential local arrangements. The concept was pioneered by G.M. Wilson in 1964¹ and subsequently extended by Renon and Prausnitz (1968)² with the NRTL model, and later by Chen and co-workers for electrolyte systems (eNRTL)³.

The fundamental assumption underlying all local composition models is that the local concentration of molecules around a central molecule differs from the bulk concentration due to differences in interaction energies between molecular pairs. This local composition approach provides a more physically realistic description of liquid mixtures, particularly for systems exhibiting significant non-ideality.

Wilson

The Wilson model was the first local composition model for activity coefficients. It is based on the concept that the local mole fraction of component \(j\) around a central molecule of component \(i\) differs from the bulk mole fraction.

Wilson model equations

The excess Gibbs energy for a multicomponent system according to the Wilson model is expressed as:

\[ \frac{G^E}{RT} = -\sum_{i=1}^{n} x_i \ln\left(\sum_{j=1}^{n} x_j \Lambda_{ij}\right) \]

where \(x_i\) and \(x_j\) are the mole fractions of components \(i\) and \(j\), respectively, and \(\Lambda_{ij}\) are the Wilson parameters.

The binary interaction parameters \(\Lambda_{ij}\) for the Wilson model are defined as:

\[ \Lambda_{ij} = \frac{V_j}{V_i} \exp\left(-\frac{\lambda_{ij} - \lambda_{ii}}{RT}\right) = \frac{V_j}{V_i} \exp\left(-\frac{\Delta\lambda_{ij}}{RT}\right) \]

where \(V_j\) and \(V_i\) are the molar volumes of component \(j\) and \(i\) in the liquid phase, \(\lambda_{ij}\) and \(\lambda_{ii}\) are the binary energy parameters representing the interaction energy between molecules \(i\) and \(j\), and \(\Delta\lambda_{ij} = \lambda_{ij} - \lambda_{ii}\) is the energy difference parameter. Note that \(\Lambda_{ii} = \Lambda_{jj} = 1\) by definition.

The activity coefficient for component \(i\) in a multicomponent mixture is given by:

\[ \ln \gamma_i = 1 - \ln\left(\sum_{j=1}^{n} x_j \Lambda_{ij}\right) - \sum_{k=1}^{n} \frac{x_k \Lambda_{ki}}{\sum_{j=1}^{n} x_j \Lambda_{kj}} \]

For a binary system, the Wilson activity coefficient equations simplify to:

\[ \ln \gamma_1 = -\ln(x_1 + \Lambda_{12}x_2) + x_2\left(\frac{\Lambda_{12}}{x_1 + \Lambda_{12}x_2} - \frac{\Lambda_{21}}{x_2 + \Lambda_{21}x_1}\right) \]

\[ \ln \gamma_2 = -\ln(x_2 + \Lambda_{21}x_1) - x_1\left(\frac{\Lambda_{12}}{x_1 + \Lambda_{12}x_2} - \frac{\Lambda_{21}}{x_2 + \Lambda_{21}x_1}\right) \]

The activity coefficients at infinite dilution are:

\[ \ln \gamma_1^\infty = -\ln(\Lambda_{12}) + 1 - \Lambda_{21} \]

\[ \ln \gamma_2^\infty = -\ln(\Lambda_{21}) + 1 - \Lambda_{12} \]

Limitations of the Wilson model

The Wilson model has one significant limitation: it cannot predict liquid-liquid phase separation (liquid-liquid equilibrium, LLE). This is because the mathematical form of the Wilson equation always predicts a single liquid phase. For systems that exhibit liquid-liquid immiscibility, the NRTL or UNIQUAC models should be used instead.

NRTL

The Non-Random Two-Liquid (NRTL) model was developed as an extension of Wilson's local composition concept. The key innovation of NRTL is the introduction of a non-randomness parameter (\(\alpha\)) that characterizes the tendency of molecules to be distributed in a non-random fashion.

NRTL model equations

The excess Gibbs energy for a multicomponent system is expressed as:

\[ \frac{G^E}{RT} = \sum_{i=1}^{n} x_i \frac{\sum_{j=1}^{n} \tau_{ji} G_{ji} x_j}{\sum_{k=1}^{n} G_{ki} x_k} \]

where the dimensionless interaction parameter \(\tau_{ij}\) and the non-randomness factor \(G_{ij}\) are defined as:

\[ \tau_{ij} = \frac{g_{ij} - g_{jj}}{RT} = \frac{\Delta g_{ij}}{RT} \]

\[ G_{ij} = \exp(-\alpha_{ij} \tau_{ij}) \]

Here, \(g_{ij}\) represents the interaction energy between molecules \(i\) and \(j\), \(\Delta g_{ij} = g_{ij} - g_{jj}\) is the energy difference parameter, and \(\alpha_{ij}\) is the non-randomness parameter (typically \(\alpha_{ij} = \alpha_{ji}\)).

The activity coefficient for component \(i\) in a multicomponent mixture is:

\[ \ln \gamma_i = \frac{\sum_{j=1}^{n} \tau_{ji} G_{ji} x_j}{\sum_{k=1}^{n} G_{ki} x_k} + \sum_{j=1}^{n} \frac{x_j G_{ij}}{\sum_{k=1}^{n} G_{kj} x_k} \left(\tau_{ij} - \frac{\sum_{m=1}^{n} x_m \tau_{mj} G_{mj}}{\sum_{k=1}^{n} G_{kj} x_k}\right) \]

For a binary system, the NRTL activity coefficient equations become:

\[ \ln \gamma_1 = x_2^2 \left[\tau_{21}\left(\frac{G_{21}}{x_1 + x_2 G_{21}}\right)^2 + \frac{\tau_{12} G_{12}}{(x_2 + x_1 G_{12})^2}\right] \]

\[ \ln \gamma_2 = x_1^2 \left[\tau_{12}\left(\frac{G_{12}}{x_2 + x_1 G_{12}}\right)^2 + \frac{\tau_{21} G_{21}}{(x_1 + x_2 G_{21})^2}\right] \]

where:

\[ G_{12} = \exp(-\alpha_{12} \tau_{12}), \quad G_{21} = \exp(-\alpha_{12} \tau_{21}) \]

The activity coefficients at infinite dilution for binary systems are:

\[ \ln \gamma_1^\infty = \tau_{21} + \tau_{12} \exp(-\alpha_{12} \tau_{12}) \]

\[ \ln \gamma_2^\infty = \tau_{12} + \tau_{21} \exp(-\alpha_{12} \tau_{21}) \]

Temperature-dependent NRTL parameters

For improved accuracy over wide temperature ranges, the interaction parameters \(\tau_{ij}\) can be expressed as temperature-dependent functions. A common form used in process simulators (e.g., Aspen Plus, AVEVA PRO/II, DWSIM) is:

\[ \tau_{ij} = a_{ij} + \frac{b_{ij}}{T} + e_{ij} \ln(T) + f_{ij} T \]

\[ \alpha_{ij} = c_{ij} + d_{ij}(T - 273.15) \]

where \(a_{ij}\), \(b_{ij}\), \(c_{ij}\), \(d_{ij}\), \(e_{ij}\), and \(f_{ij}\) are adjustable parameters.

The non-randomness parameter \(\alpha\)

The non-randomness parameter \(\alpha_{ij}\) characterizes the tendency of molecules \(i\) and \(j\) to be distributed in a non-random fashion. When \(\alpha_{ij} = 0\), the model reduces to the random mixing case (similar to the Margules two-suffix equation).

Typical values of \(\alpha\) depend on the molecular nature of the mixture:

System type	Recommended \(\alpha\)
Nonpolar mixtures	0.20
Saturated hydrocarbons + polar non-associated	0.30
Nonpolar + polar non-associated	0.30
Water + polar non-associated	0.40
Aqueous organic solutions	0.20
Strongly self-associated + nonpolar	0.47

In practice, \(\alpha\) is often fixed at 0.30 for vapor-liquid equilibrium (VLE) calculations and 0.20 for liquid-liquid equilibrium (LLE) calculations to reduce the number of adjustable parameters.

eNRTL

The electrolyte Non-Random Two-Liquid (eNRTL) model extends the NRTL framework to describe the thermodynamic properties of electrolyte solutions. The model accounts for both long-range electrostatic interactions between ions and short-range interactions between all species (ions and molecules).

Theoretical framework

The eNRTL model expresses the excess Gibbs energy as the sum of two contributions:

\[ G^{E,*} = G^{E,*}_{PDH} + G^{E,*}_{NRTL} \]

where \(G^{E,*}_{PDH}\) represents the long-range electrostatic contribution (Pitzer-Debye-Hückel) and \(G^{E,*}_{NRTL}\) represents the short-range local composition contribution. The asterisk (\(*\)) denotes the unsymmetric convention (infinite dilution reference state for ions).

Long-range interaction contribution: Pitzer-Debye-Hückel

The long-range electrostatic interactions between ions are described by the Pitzer-Debye-Hückel equation:

\[ \frac{G^{E,*}_{PDH}}{RT} = -\left(\sum_i x_i\right) \frac{4 A_\phi I_x}{\rho} \ln\left(1 + \rho I_x^{1/2}\right) \]

where \(A_\phi\) is the Debye-Hückel parameter, \(I_x\) is the ionic strength on a mole fraction basis, and \(\rho\) is the closest approach parameter (typically set to 14.9 for aqueous systems).

The Debye-Hückel parameter is defined as:

\[ A_\phi = \frac{1}{3}\left(\frac{2\pi N_A d_s}{1000}\right)^{1/2} \left(\frac{e^2}{\varepsilon_s k_B T}\right)^{3/2} \]

where \(N_A\) is Avogadro's number, \(d_s\) is the solvent density, \(e\) is the electronic charge, \(\varepsilon_s\) is the dielectric constant of the solvent, and \(k_B\) is Boltzmann's constant.

Short-range interaction contribution: Local composition NRTL

The short-range contribution follows the NRTL formalism with two key assumptions for electrolyte systems:

1. Like-ion repulsion assumption: Due to strong electrostatic repulsion, the local composition of cations around cations (and anions around anions) is zero.

2. Local electroneutrality assumption: The local composition of cations and anions around a central solvent molecule satisfies local electroneutrality.

The short-range excess Gibbs energy is expressed as:

\[ \frac{G^{E,*}_{NRTL}}{RT} = \sum_m x_m \frac{\sum_j X_j G_{jm} \tau_{jm}}{\sum_k X_k G_{km}} + \sum_c x_c \sum_{a'} \frac{X_{a'}}{\sum_{a''} X_{a''}} \frac{\sum_j X_j G_{jc,a'c} \tau_{jc,a'c}}{\sum_k X_k G_{kc,a'c}} + \sum_a x_a \sum_{c'} \frac{X_{c'}}{\sum_{c''} X_{c''}} \frac{\sum_j X_j G_{ja,c'a} \tau_{ja,c'a}}{\sum_k X_k G_{ka,c'a}} \]

where the subscripts \(m\), \(c\), and \(a\) refer to molecular species, cations, and anions, respectively.

Binary interaction parameters

The eNRTL model requires the following binary interaction parameters:

Interaction type	Parameters
Molecule-molecule	\(\tau_{mm'}\), \(\alpha_{mm'}\) (standard NRTL)
Molecule-electrolyte	\(\tau_{m,ca}\), \(\tau_{ca,m}\), \(\alpha_{m,ca}\)
Electrolyte-electrolyte	\(\tau_{ca,c'a}\), \(\tau_{ca,ca'}\), \(\alpha_{ca,c'a}\)

The non-randomness factor for electrolyte interactions is typically set to \(\alpha = 0.2\).

eNRTL parameter table

Nonrandomness Factor of \(\alpha\) = 0.2

Molecule (1)	Water	Hexane	Methanol
Electrolyte (2)	NaCl	NaCl	NaCl
\(\tau_{12}\)	8.885 (a)	15.000 (b)	3.624 (c)
\(\tau_{21}\)	-4.549 (a)	5.000 (b)	-0.789 (c)

(a) Chen et al., 1982³ (b) Chen and Song, 2004⁴ (c) Yang and Lee, 1998⁵

Activity coefficient expressions

The activity coefficient for a molecular species \(m\) in the eNRTL model is:

\[ \ln \gamma_m^* = \ln \gamma_m^{*,PDH} + \ln \gamma_m^{*,NRTL} \]

For cation \(c\) and anion \(a\), similar expressions apply with appropriate modifications for the reference state convention.

COSMO-NRTL: Deriving NRTL parameters from COSMO-RS

A very useful application of COSMO-RS and COSMO-SAC models is their ability to generate NRTL binary interaction parameters for systems where experimental data is unavailable. This hybrid approach, often referred to as COSMO-NRTL, leverages the predictive capability of quantum chemistry-based models to extend the applicability of semi-empirical models commonly used in process simulation software.

Motivation for COSMO-NRTL

While COSMO-RS and COSMO-SAC provide excellent predictive capabilities directly, many industrial process simulation tools (such as Aspen Plus, HYSYS, AVEVA PRO/II or DWSIM) are built around the NRTL, UNIQUAC/UNIFAC, or Wilson activity coefficient models. The COSMO-NRTL approach bridges this gap by:

1. Using COSMO-RS/SAC to generate pseudo-experimental activity coefficient data
2. Regressing NRTL parameters against this generated data
3. Enabling seamless integration with existing process simulation infrastructure

This methodology is particularly valuable for biobased molecules, pharmaceuticals, and novel chemical systems where experimental phase equilibrium data is limited or unavailable.

Regression methodology

The NRTL parameter regression from COSMO-RS data follows these steps:

Step 1: Generate reference activity coefficient data

Using COSMO-RS or COSMO-SAC, calculate activity coefficients \(\gamma_1\) and \(\gamma_2\) for a range of compositions at the temperature(s) of interest:

\[ \gamma_i^{COSMO}(x_1, T) \quad \text{for} \quad x_1 \in [0, 1] \]

The composition grid should include sufficient points to capture the non-linear behavior of the activity coefficients, including the limiting cases at infinite dilution (\(x_1 \to 0\) and \(x_1 \to 1\)).

Step 2: Define the objective function

The NRTL parameters are obtained by minimizing the difference between COSMO-predicted and NRTL-calculated activity coefficients:

\[ \text{min} \sum_k \left[ w_1 \left(\ln \gamma_1^{COSMO}(x_1^k) - \ln \gamma_1^{NRTL}(x_1^k)\right)^2 + w_2 \left(\ln \gamma_2^{COSMO}(x_1^k) - \ln \gamma_2^{NRTL}(x_1^k)\right)^2 \right] \]

where \(w_1\) and \(w_2\) are weighting factors that can be adjusted to emphasize infinite dilution values.

Step 3: Optimization

Common optimization algorithms for NRTL parameter regression include:

• Differential evolution: A global optimization method that is robust for the non-convex NRTL parameter space
• L-BFGS-B: A local optimization method suitable when good initial guesses are available
• Hybrid approaches: Global search followed by local refinement

Step 4: Validation

The regressed parameters should be validated by comparing NRTL predictions with the original COSMO data and, when available, experimental VLE or LLE data.

Implementation with Cebule SDK

The Cebule engine provides integrated support for COSMO-NRTL parameter derivation. The workflow involves:

1. Create COSMO calculations for all components
2. Generate sigma profiles using the SIGMA task
3. Calculate activity coefficients at various compositions using the ACTIVITY_COEFFICIENT task
4. Regress NRTL parameters against the COSMO-derived data

Example workflow:

from cebule import Session, TaskType

session = Session()

# Step 1: COSMO calculations for each component
cosmo_comp1 = session.cebule.create_task("COSMO Component 1",
    TaskType.COSMO,
    cosmo_input_comp1,
    nprocs=16)

cosmo_comp2 = session.cebule.create_task("COSMO Component 2",
    TaskType.COSMO,
    cosmo_input_comp2,
    nprocs=16)

# Step 2: Generate sigma profiles
sigma_comp1 = session.cebule.create_task("Sigma Profile 1",
    TaskType.SIGMA,
    connected_task_id=cosmo_comp1.id,
    cosmo_method="cosmo-rs")

sigma_comp2 = session.cebule.create_task("Sigma Profile 2",
    TaskType.SIGMA,
    connected_task_id=cosmo_comp2.id,
    cosmo_method="cosmo-rs")

# Step 3: Calculate activity coefficients at various compositions
activity_data = []
for x1 in [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]:
    task = session.cebule.create_task("Activity Coefficient",
        TaskType.ACTIVITY_COEFFICIENT,
        connected_task_id=[sigma_comp1.id, sigma_comp2.id],
        calculation_method="cosmo-rs",
        mole_fractions=[x1, 1.0 - x1],
        temperature=298.15)
    activity_data.append(task)

# Step 4: Collect results and perform NRTL regression
# (Use scipy.optimize or similar for parameter fitting)

Best practices for COSMO-NRTL parameter regression

1. Composition range: Include compositions from 0 to 1 with emphasis on the infinite dilution regions (e.g., \(x_1 = 0.001, 0.01, 0.99, 0.999\))

2. Temperature dependence: For temperature-dependent NRTL parameters, generate COSMO data at multiple temperatures and fit the extended parameter form

3. Non-randomness factor: The value of \(\alpha\) can either be:

4. Fixed at a typical value (0.2 or 0.3)

5. Included as an additional fitted parameter

6. Weighting: Apply higher weights to infinite dilution activity coefficients if accurate limiting behavior is important for the application

7. Validation: Always compare the regressed NRTL model with the original COSMO predictions to ensure the regression quality

Advantages of COSMO-NRTL

Advantage	Description
Predictive	No experimental VLE data required
Consistent	Same COSMO model underlies all parameter derivations
Extensible	Works for any molecule that can be calculated with DFT/COSMO
Compatible	Output parameters directly usable in commercial process simulators
Efficient	Once sigma profiles are computed, activity coefficients are fast to calculate

Limitations

• The accuracy of derived NRTL parameters is limited by the accuracy of the underlying COSMO-RS/SAC model
• For highly associating systems or electrolytes, additional model refinements may be needed
• The NRTL functional form may not perfectly capture complex COSMO-RS predictions, leading to some residual error

---

1. Grant M. Wilson. Vapor-liquid equilibrium. xi. a new expression for the excess free energy of mixing. Journal of the American Chemical Society, 86:127–130, 1964. URL: https://api.semanticscholar.org/CorpusID:102136993. ↩

2. Henri Renon and J. M. Prausnitz. Local compositions in thermodynamic excess functions for liquid mixtures. AIChE Journal, 14(1):135–144, 1968. doi:https://doi.org/10.1002/aic.690140124. ↩

3. Chau-Chyun Chen, H. I. Britt, J. F. Boston, and L. B. Evans. Local composition model for excess gibbs energy of electrolyte systems. part i: single solvent, single completely dissociated electrolyte systems. AIChE Journal, 28(4):588–596, 1982. doi:https://doi.org/10.1002/aic.690280410. ↩↩

4. Chau-Chyun Chen and Yuhua Song. Generalized electrolyte-nrtl model for mixed-solvent electrolyte systems. AIChE Journal, 50(8):1928–1941, 2004. doi:https://doi.org/10.1002/aic.10151. ↩

5. Sung-Oh Yang and Chul Soo Lee. Vapor−liquid equilibria of water + methanol in the presence of mixed salts. Journal of Chemical & Engineering Data, 43(4):558–561, 1998. doi:10.1021/je970286w. ↩