Association

The Helmholtz contribution due to short range attractive interaction (“association”) in SAFT models can be conveniently expressed as

\[\frac{A^\mathrm{assoc}}{k_\mathrm{B}TV}=\sum_\alpha\rho_\alpha\left(\ln X_\alpha-\frac{X_\alpha}{2}+\frac{1}{2}\right)\]

Here, \(\alpha\) is the index of all distinguishable association sites in the system. The site density \(\rho_\alpha\) is the density of the segment or molecule on which the association site is located times the multiplicity of the site. The fraction of non-bonded association sites \(X_\alpha\) is calculated implicitly from

(1)\[\frac{1}{X_\alpha}=1+\sum_\beta\rho_\beta X_\beta\Delta^{\alpha\beta}\]

where \(\Delta^{\alpha\beta}\) is the association strength between sites \(\alpha\) and \(\beta\). The exact expression for the association strength varies between different SAFT versions. We implement the expressions used in PC-SAFT but a generalization to other models is straightforward.

The number of degrees of freedom in the association strength matrix \(\Delta\) is vast and for practical purposes it is useful to reduce the number of non-zero elements in \(\Delta\). In \(\text{FeO}_\text{s}\) we use 3 kinds of association sites: \(A\), \(B\) and \(C\). Association sites of kind \(A\) only associate with \(B\) and vice versa, whereas association sites of kind \(C\) only associate with other sites of kind \(C\). By sorting the sites by their kind, the entire \(\Delta\) matrix can be reduced to two smaller matrices: \(\Delta^{AB}\) and \(\Delta^{CC}\).

In practice, the \(AB\) association can be linked to hydrogen bonding sites. The \(CC\) association is less widely used but implemented to ensure that all the association schemes defined in Huang and Radosz 1990 are covered.

Calculation of the fraction of non-bonded association sites

The algorithm to solve the fraction of non-bonded sites for general association schemes follows Michelsen 2006. The problem is rewritten as an optimization problem with gradients

\[g_\alpha=\frac{1}{X_\alpha}-1-\sum_\beta\rho_\beta X_\beta\Delta^{\alpha\beta}\]

The analytic Hessian is

\[H_{\alpha\beta}=-\frac{\delta_{\alpha\beta}}{X_\alpha^2}-\rho_\beta\Delta^{\alpha\beta}\]

with the Kronecker delta \(\delta_{\alpha\beta}\). However, Michelsen 2006 shows that it is beneficial to use

\[\hat{H}_{\alpha\beta}=-\frac{\delta_{\alpha\beta}}{X_\alpha}\left(1+\sum_\beta\rho_\beta X_\beta\Delta^{\alpha\beta}\right)-\rho_\beta\Delta^{\alpha\beta}\]

instead. \(X_\alpha\) can then be solved robustly using a Newton algorithm with a check that ensures that the values of \(X_\alpha\) remain positive. With the split in \(AB\) and \(CC\) association the two kinds different versions could be solved separately from each other. This is currently not implemented, because use cases are rare to nonexistent and the benefit is small.

A drastic improvement in performance, however, can be achieved by solving eq. (1) analytically for simple cases. If there is only one \(A\) and one \(B\) site the corresponding fractions of non-bonded association sites can be calculated from

\[X_A=\frac{2}{\sqrt{\left(1+\left(\rho_A-\rho_B\right)\Delta^{AB}\right)^2+4\rho_B\Delta^{AB}}+1+\left(\rho_B-\rho_A\right)\Delta^{AB}}\]

\[X_B=\frac{2}{\sqrt{\left(1+\left(\rho_A-\rho_B\right)\Delta^{AB}\right)^2+4\rho_B\Delta^{AB}}+1+\left(\rho_A-\rho_B\right)\Delta^{AB}}\]

The specific form of the expressions (with the square root in the denominator) is particularly robust for cases where \(\Delta\) and/or \(\rho\) are small or even 0.

Analogously, for a single \(C\) site, the expression becomes

\[X_C=\frac{2}{\sqrt{1+4\rho_C\Delta^{CC}}+1}\]

PC-SAFT expression for the association strength

In \(\text{FeO}_\text{s}\) every association site \(\alpha\) is parametrized with the dimensionless association volume \(\kappa^\alpha\) and the association energy \(\varepsilon^\alpha\). The association strength between sites \(\alpha\) and \(\beta\) is then calculated from

\[\Delta^{\alpha\beta}=\left(\frac{1}{1-\zeta_3}+\frac{\frac{3}{2}d_{ij}\zeta_2}{\left(1-\zeta_3\right)^2}+\frac{\frac{1}{2}d_{ij}^2\zeta_2^2}{\left(1-\zeta_3\right)^3}\right)\sqrt{\sigma_i^3\kappa^\alpha\sigma_j^3\kappa^\beta}\left(e^{\frac{\varepsilon^\alpha+\varepsilon^\beta}{2k_\mathrm{B}T}}-1\right)\]

with

\[d_{ij}=\frac{2d_id_j}{d_i+d_j},~~~~d_k=\sigma_k\left(1-0.12e^{\frac{-3\varepsilon_k}{k_\mathrm{B}T}}\right)\]

and

\[\zeta_2=\frac{\pi}{6}\sum_k\rho_km_kd_k^2,~~~~\zeta_3=\frac{\pi}{6}\sum_k\rho_km_kd_k^3\]

The indices \(i\), \(j\) and \(k\) are used to index molecules or segments rather than sites. \(i\) and \(j\) in particular refer to the molecule or segment that contain site \(\alpha\) or \(\beta\) respectively.

On their own the parameters \(\kappa^\alpha\) and \(\varepsilon^\beta\) have no physical meaning. For a pure system of self-associating molecules it is therefore natural to define \(\kappa^A=\kappa^B\equiv\kappa^{A_iB_i}\) and \(\varepsilon^A=\varepsilon^B\equiv\varepsilon^{A_iB_i}\) where \(\kappa^{A_iB_i}\) and \(\varepsilon^{A_iB_i}\) are now parameters describing the molecule rather than individual association sites. However, for systems that are not self-associating, the more generic parametrization is required.

SAFT-VR Mie expression for the association strength

We provide an implementation of the association strength as published by Lafitte et al. (2013). Every association site \(\alpha\) is parametrized with two distances \(r_c^\alpha\) and \(r_d^\alpha\), and the association energy \(\varepsilon^\alpha\). Whereas \(r_c^\alpha\) is a parameter that is adjusted for each substance, \(r_d^\alpha\) is kept constant as \(r_d^\alpha = 0.4 \sigma\). Note that the parameter \(r_c^\alpha\) has to be provided as dimensionless quantity in the input, i.e. divided by the segment’s \(\sigma\) value.

The association strength between sites \(\alpha\) and \(\beta\) is then calculated from

\[\Delta^{\alpha\beta}=\left(\frac{1}{1-\zeta_3}+\frac{\frac{3}{2}\tilde{d}_{ij}\zeta_2}{\left(1-\zeta_3\right)^2}+\frac{\frac{1}{2}\tilde{d}_{ij}^2\zeta_2^2}{\left(1-\zeta_3\right)^3}\right)K_{ij}^{\alpha\beta}\sigma_{ij}^3\left(e^{\frac{\sqrt{\varepsilon^\alpha\varepsilon^\beta}}{k_\mathrm{B}T}}-1\right)\]

with

\[\tilde{d}_{ij}=\frac{2d_id_j}{d_i+d_j},~~~~\sigma_{ij} = \frac{\sigma_i + \sigma_j}{2}\]

and

\[\zeta_2=\frac{\pi}{6}\sum_k\rho_km_kd_k^2,~~~~\zeta_3=\frac{\pi}{6}\sum_k\rho_km_kd_k^3\]

where

\[d_i = \int_{l}^{\sigma_i} \left[1 - e^{-\beta u_i^\text{Mie}(r)}\right]\mathrm{d}r.\]

The integral is solved using a Gauss-Legendre quadrature of fifth order where the lower limit, \(l\), is determined using the method of Aasen et al. The dimensionless bonding volume \(K^{\alpha\beta}_{ij}\) (see A43) utilizes arithmetic combining rules for \(r_c^{\alpha\beta}\), \(r_d^{\alpha\beta}\) and \(d_{ij}\) of unlike sites and segments, respectively.

Helmholtz energy functional

The Helmholtz energy contribution proposed by Yu and Wu 2002 is used to model association in inhomogeneous systems. It uses the same weighted densities that are used in Fundamental Measure Theory (the White-Bear version). The Helmholtz energy density is then calculated from

\[\beta f=\sum_\alpha\frac{n_0^\alpha\xi_i}{m_i}\left(\ln X_\alpha-\frac{X_\alpha}{2}+\frac{1}{2}\right)\]

with the equations for the fraction of non-bonded association sites adapted to

\[\frac{1}{X_\alpha}=1+\sum_\beta\frac{n_0^\beta\xi_j}{m_j}X_\beta\Delta^{\alpha\beta}\]

The association strength, again using the PC-SAFT parametrization, is

\[\Delta^{\alpha\beta}=\left(\frac{1}{1-n_3}+\frac{\frac{1}{4}d_{ij}n_2\xi}{\left(1-n_3\right)^2}+\frac{\frac{1}{72}d_{ij}^2n_2^2\xi}{\left(1-n_3\right)^3}\right)\sqrt{\sigma_i^3\kappa^\alpha\sigma_j^3\kappa^\beta}\left(e^{\frac{\varepsilon^\alpha+\varepsilon^\beta}{2k_\mathrm{B}T}}-1\right)\]

The quantities \(\xi\) and \(\xi_i\) were introduced to account for the effect of inhomogeneity and are defined as

\[\xi=1-\frac{\vec{n}_2\cdot\vec{n}_2}{n_2^2},~~~~\xi_i=1-\frac{\vec{n}_2^i\cdot\vec{n}_2^i}{{n_2^i}^2}\]