diff --git a/CHANGELOG.md b/CHANGELOG.md index d8740a33f..87df11d7f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [Unreleased] +### Changed +- Changed the internal implementation of the association contribution to accomodate more general association schemes. [#150](https://github.com/feos-org/feos/pull/150) +- To comply with the new association implementation, the default values of `na` and `nb` are now `0` rather than `1`. Parameter files have been adapted accordingly. [#150](https://github.com/feos-org/feos/pull/150) ## [0.4.2] - 2023-03-20 - Python only: Release the changes introduced in `feos-core` 0.4.1 and `feos-dft` 0.4.1. diff --git a/docs/theory/models/FeOs_Association.png b/docs/theory/models/FeOs_Association.png new file mode 100644 index 000000000..dd43fd2f9 Binary files /dev/null and b/docs/theory/models/FeOs_Association.png differ diff --git a/docs/theory/models/association.md b/docs/theory/models/association.md new file mode 100644 index 000000000..9313bbe46 --- /dev/null +++ b/docs/theory/models/association.md @@ -0,0 +1,83 @@ +# 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 + +$$\frac{1}{X_\alpha}=1+\sum_\beta\rho_\beta X_\beta\Delta^{\alpha\beta}$$ (eqn:x_assoc) + +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}$. + +```{image} FeOs_Association.png +:alt: Association strength matrix +:class: bg-primary +:width: 300px +:align: center +``` + +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](https://pubs.acs.org/doi/10.1021/ie00107a014) 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](https://pubs.acs.org/doi/full/10.1021/ie060029x). 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](https://pubs.acs.org/doi/full/10.1021/ie060029x) 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. {eq}`eqn:x_assoc` 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. + +## Helmholtz energy functional +The Helmholtz energy contribution proposed by [Yu and Wu 2002](https://aip.scitation.org/doi/abs/10.1063/1.1463435) is used to model association in inhomogeneous systems. It uses the same weighted densities that are used in [Fundamental Measure Theory](hard_spheres) (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}$$ \ No newline at end of file diff --git a/docs/theory/models/hard_spheres.md b/docs/theory/models/hard_spheres.md new file mode 100644 index 000000000..26c2a278f --- /dev/null +++ b/docs/theory/models/hard_spheres.md @@ -0,0 +1,57 @@ +# Hard spheres + +$\text{FeO}_\text{s}$ provides an implementation of the Boublík-Mansoori-Carnahan-Starling-Leland (BMCSL) equation of state ([Boublík, 1970](https://doi.org/10.1063/1.1673824), [Mansoori et al., 1971](https://doi.org/10.1063/1.1675048)) for hard-sphere mixtures which is often used as reference contribution in SAFT equations of state. The implementation is generalized to allow the description of non-sperical or fused-sphere reference fluids. + +The reduced Helmholtz energy density is calculated according to + +$$\frac{\beta A}{V}=\frac{6}{\pi}\left(\frac{3\zeta_1\zeta_2}{1-\zeta_3}+\frac{\zeta_2^3}{\zeta_3\left(1-\zeta_3\right)^2}+\left(\frac{\zeta_2^3}{\zeta_3^2}-\zeta_0\right)\ln\left(1-\zeta_3\right)\right)$$ + +with the packing fractions + +$$\zeta_k=\frac{\pi}{6}\sum_\alpha C_{k,\alpha}\rho_\alpha d_\alpha^k,~~~~~~~~k=0\ldots 3$$ + +The geometry coefficients $C_{k,\alpha}$ and the segment diameters $d_\alpha$ depend on the context in which the model is used. The following table shows how the expression can be reused in various models. For details on the fused-sphere chain model, check the [repository](https://github.com/feos-org/feos-fused-chains) or the [publication](https://doi.org/10.1103/PhysRevE.105.034110). + +||Hard spheres|PC-SAFT|Fused-sphere chains| +|-|:-:|:-:|:-:| +|$d_\alpha$|$\sigma_\alpha$|$\sigma_\alpha\left(1-0.12e^{\frac{-3\varepsilon_\alpha}{k_\mathrm{B}T}}\right)$|$\sigma_\alpha$| +|$C_{0,\alpha}$|$1$|$m_\alpha$|$1$| +|$C_{1,\alpha}$|$1$|$m_\alpha$|$A_\alpha^*$| +|$C_{1,\alpha}$|$1$|$m_\alpha$|$A_\alpha^*$| +|$C_{1,\alpha}$|$1$|$m_\alpha$|$V_\alpha^*$| + +## Fundamental measure theory + +An model for inhomogeneous mixtues of hard spheres is provided by fundamental measure theory (FMT, [Rosenfeld, 1989](https://doi.org/10.1103/PhysRevLett.63.980)). Different variants have been proposed of which only those that are consistent with the BMCSL equation of state in the homogeneous limit are currently considered in $\text{FeO}_\text{s}$ (exluding, e.g., the original Rosenfeld and White-Bear II variants). + +The Helmholtz energy density is calculated according to + +$$\beta f=-n_0\ln\left(1-n_3\right)+\frac{n_{12}}{1-n_3}+\frac{1}{36\pi}n_2n_{22}f_3(n_3)$$ + +The expressions for $n_{12}$ and $n_{22}$ depend on the FMT version. + +|version|$n_{12}$|$n_{22}$|references| +|-|:-:|:-:|:-:| +|WhiteBear|$n_1n_2-\vec n_1\cdot\vec n_2$|$n_2^2-3\vec n_2\cdot\vec n_2$|[Roth et al., 2002](https://doi.org/10.1088/0953-8984/14/46/313), [Yu and Wu, 2002](https://doi.org/10.1063/1.1520530)| +|KierlikRosinberg|$n_1n_2$|$n_2^2$|[Kierlik and Rosinberg, 1990](https://doi.org/10.1103/PhysRevA.42.3382)| +|AntiSymWhiteBear|$n_1n_2-\vec n_1\cdot\vec n_2$|$n_2^2\left(1-\frac{\vec n_2\cdot\vec n_2}{n_2^2}\right)^3$|[Rosenfeld et al., 1997](https://doi.org/10.1103/PhysRevE.55.4245), [Kessler et al., 2021](https://doi.org/10.1016/j.micromeso.2021.111263)| + +For small $n_3$, the value of $f(n_3)$ numerically diverges. Therefore, it is approximated with a Taylor expansion. + +$$f_3=\begin{cases}\frac{n_3+\left(1-n_3\right)^2\ln\left(1-n_3\right)}{n_3^2\left(1-n_3\right)^2}&\text{if }n_3>10^{-5}\\ +\frac{3}{2}+\frac{8}{3}n_3+\frac{15}{4}n_3^2+\frac{24}{5}n_3^3+\frac{35}{6}n_3^4&\text{else}\end{cases}$$ + +The weighted densities $n_k(\mathbf{r})$ are calculated by convolving the density profiles $\rho_\alpha(\mathbf{r})$ with weight functions $\omega_k^\alpha(\mathbf{r})$ + +$$n_k(\mathbf{r})=\sum_\alpha n_k^\alpha(\mathbf{r})=\sum_\alpha\int\rho_\alpha(\mathbf{r}')\omega_k^\alpha(\mathbf{r}-\mathbf{r}')\mathrm{d}\mathbf{r}'$$ + +which differ between the different FMT versions. + +||WhiteBear/AntiSymWhiteBear|KierlikRosinberg| +|-|:-:|:-:| +|$\omega_0^\alpha(\mathbf{r})$|$\frac{C_{0,\alpha}}{\pi\sigma_\alpha^2}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|$C_{0,\alpha}\left(-\frac{1}{8\pi}\,\delta''\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)+\frac{1}{2\pi\|\mathbf{r}\|}\,\delta'\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)\right)$| +|$\omega_1^\alpha(\mathbf{r})$|$\frac{C_{1,\alpha}}{2\pi\sigma_\alpha}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|$\frac{C_{1,\alpha}}{8\pi}\,\delta'\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$| +|$\omega_2^\alpha(\mathbf{r})$|$C_{2,\alpha}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|$C_{2,\alpha}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$| +|$\omega_3^\alpha(\mathbf{r})$|$C_{3,\alpha}\,\Theta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|$C_{3,\alpha}\,\Theta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$| +|$\vec\omega_1^\alpha(\mathbf{r})$|$C_{3,\alpha}\frac{\mathbf{r}}{2\pi\sigma_\alpha\|\mathbf{r}\|}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|-| +|$\vec\omega_2^\alpha(\mathbf{r})$|$C_{3,\alpha}\frac{\mathbf{r}}{\|\mathbf{r}\|}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|-| \ No newline at end of file diff --git a/docs/theory/models/index.md b/docs/theory/models/index.md index 148f360b4..efccaf1cd 100644 --- a/docs/theory/models/index.md +++ b/docs/theory/models/index.md @@ -6,4 +6,7 @@ It is currently still under construction. You can help by [contributing](https:/ ```{eval-rst} .. toctree:: :maxdepth: 1 + + hard_spheres + association ``` \ No newline at end of file diff --git a/parameters/pcsaft/gross2002.json b/parameters/pcsaft/gross2002.json index b888201d8..3b8737154 100644 --- a/parameters/pcsaft/gross2002.json +++ b/parameters/pcsaft/gross2002.json @@ -13,7 +13,9 @@ "sigma": 3.23, "epsilon_k": 188.9, "kappa_ab": 0.035176, - "epsilon_k_ab": 2899.5 + "epsilon_k_ab": 2899.5, + "na": 1.0, + "nb": 1.0 }, "molarweight": 32.042 }, @@ -31,7 +33,9 @@ "sigma": 3.1771, "epsilon_k": 198.24, "kappa_ab": 0.032384, - "epsilon_k_ab": 2653.4 + "epsilon_k_ab": 2653.4, + "na": 1.0, + "nb": 1.0 }, "molarweight": 46.069 }, @@ -49,7 +53,9 @@ "sigma": 3.2522, "epsilon_k": 233.4, "kappa_ab": 0.015268, - "epsilon_k_ab": 2276.8 + "epsilon_k_ab": 2276.8, + "na": 1.0, + "nb": 1.0 }, "molarweight": 60.096 }, @@ -67,7 +73,9 @@ "sigma": 3.6139, "epsilon_k": 259.59, "kappa_ab": 0.006692, - "epsilon_k_ab": 2544.6 + "epsilon_k_ab": 2544.6, + "na": 1.0, + "nb": 1.0 }, "molarweight": 74.123 }, @@ -85,7 +93,9 @@ "sigma": 3.4508, "epsilon_k": 247.28, "kappa_ab": 0.010319, - "epsilon_k_ab": 2252.1 + "epsilon_k_ab": 2252.1, + "na": 1.0, + "nb": 1.0 }, "molarweight": 88.15 }, @@ -103,7 +113,9 @@ "sigma": 3.6735, "epsilon_k": 262.32, "kappa_ab": 0.005747, - "epsilon_k_ab": 2538.9 + "epsilon_k_ab": 2538.9, + "na": 1.0, + "nb": 1.0 }, "molarweight": 102.177 }, @@ -121,7 +133,9 @@ "sigma": 3.545, "epsilon_k": 253.46, "kappa_ab": 0.001155, - "epsilon_k_ab": 2878.5 + "epsilon_k_ab": 2878.5, + "na": 1.0, + "nb": 1.0 }, "molarweight": 116.203 }, @@ -139,7 +153,9 @@ "sigma": 3.7145, "epsilon_k": 262.74, "kappa_ab": 0.002197, - "epsilon_k_ab": 2754.8 + "epsilon_k_ab": 2754.8, + "na": 1.0, + "nb": 1.0 }, "molarweight": 130.23 }, @@ -157,7 +173,9 @@ "sigma": 3.7292, "epsilon_k": 263.64, "kappa_ab": 0.001427, - "epsilon_k_ab": 2941.9 + "epsilon_k_ab": 2941.9, + "na": 1.0, + "nb": 1.0 }, "molarweight": 144.257 }, @@ -175,7 +193,9 @@ "sigma": 3.2085, "epsilon_k": 208.42, "kappa_ab": 0.024675, - "epsilon_k_ab": 2253.9 + "epsilon_k_ab": 2253.9, + "na": 1.0, + "nb": 1.0 }, "molarweight": 60.096 }, @@ -193,7 +213,9 @@ "sigma": 3.9053, "epsilon_k": 266.01, "kappa_ab": 0.001863, - "epsilon_k_ab": 2618.8 + "epsilon_k_ab": 2618.8, + "na": 1.0, + "nb": 1.0 }, "molarweight": 88.15 }, @@ -211,7 +233,9 @@ "sigma": 3.0007, "epsilon_k": 366.51, "kappa_ab": 0.034868, - "epsilon_k_ab": 2500.7 + "epsilon_k_ab": 2500.7, + "na": 1.0, + "nb": 1.0 }, "molarweight": 18.015 }, @@ -229,7 +253,9 @@ "sigma": 2.8906, "epsilon_k": 214.94, "kappa_ab": 0.095103, - "epsilon_k_ab": 684.3 + "epsilon_k_ab": 684.3, + "na": 1.0, + "nb": 1.0 }, "molarweight": 31.06 }, @@ -247,7 +273,9 @@ "sigma": 3.1343, "epsilon_k": 221.53, "kappa_ab": 0.017275, - "epsilon_k_ab": 854.7 + "epsilon_k_ab": 854.7, + "na": 1.0, + "nb": 1.0 }, "molarweight": 45.09 }, @@ -265,7 +293,9 @@ "sigma": 3.5347, "epsilon_k": 250.52, "kappa_ab": 0.022674, - "epsilon_k_ab": 1028.1 + "epsilon_k_ab": 1028.1, + "na": 1.0, + "nb": 1.0 }, "molarweight": 59.11 }, @@ -283,7 +313,9 @@ "sigma": 3.4777, "epsilon_k": 231.8, "kappa_ab": 0.02134, - "epsilon_k_ab": 932.2 + "epsilon_k_ab": 932.2, + "na": 1.0, + "nb": 1.0 }, "molarweight": 59.11 }, @@ -301,7 +333,9 @@ "sigma": 3.7021, "epsilon_k": 335.47, "kappa_ab": 0.074883, - "epsilon_k_ab": 1351.6 + "epsilon_k_ab": 1351.6, + "na": 1.0, + "nb": 1.0 }, "molarweight": 93.13 }, @@ -319,7 +353,9 @@ "sigma": 3.8582, "epsilon_k": 211.59, "kappa_ab": 0.07555, - "epsilon_k_ab": 3044.4 + "epsilon_k_ab": 3044.4, + "na": 1.0, + "nb": 1.0 }, "molarweight": 60.053 } diff --git a/parameters/pcsaft/loetgeringlin2015_homo.json b/parameters/pcsaft/loetgeringlin2015_homo.json index ff9df6b16..030332969 100644 --- a/parameters/pcsaft/loetgeringlin2015_homo.json +++ b/parameters/pcsaft/loetgeringlin2015_homo.json @@ -277,6 +277,8 @@ "epsilon_k": 488.66, "epsilon_k_ab": 2517.0, "kappa_ab": 0.006825, + "na": 1.0, + "nb": 1.0, "viscosity": [ -15.7583e-3, -2.5654e-1, @@ -294,6 +296,8 @@ "epsilon_k": 467.59, "epsilon_k_ab": 1064.6, "kappa_ab": 0.026662, + "na": 1.0, + "nb": 1.0, "viscosity": [ -4.4048e-3, -0.6089e-1, diff --git a/parameters/pcsaft/loetgeringlin2018.json b/parameters/pcsaft/loetgeringlin2018.json index 989e3439c..153e40ed3 100644 --- a/parameters/pcsaft/loetgeringlin2018.json +++ b/parameters/pcsaft/loetgeringlin2018.json @@ -896,6 +896,8 @@ "mu": 0.9204, "kappa_ab": 0.03, "epsilon_k_ab": 1547.4422, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.6824, -2.129, @@ -921,6 +923,8 @@ "mu": 1.0313, "kappa_ab": 0.03, "epsilon_k_ab": 959.5053, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.3807, 2.1585, @@ -946,6 +950,8 @@ "mu": 1.07, "kappa_ab": 0.03, "epsilon_k_ab": 39.2881, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.8987, -2.1083, @@ -1040,6 +1046,8 @@ "mu": 1.391, "kappa_ab": 0.03, "epsilon_k_ab": 1497.5173, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.8691, -2.3216, @@ -1065,6 +1073,8 @@ "mu": 1.6698, "kappa_ab": 0.03, "epsilon_k_ab": 1846.2029, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.0065, -1.88906, @@ -1090,6 +1100,8 @@ "mu": 1.6608, "kappa_ab": 0.03, "epsilon_k_ab": 1718.0488, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.9864, -1.3818, @@ -1363,6 +1375,8 @@ "mu": 1.409, "kappa_ab": 0.03, "epsilon_k_ab": 2644.5966, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.3992, -3.2065, @@ -1388,6 +1402,8 @@ "mu": 1.6189, "kappa_ab": 0.03, "epsilon_k_ab": 1967.8312, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.5371, -2.97397, @@ -1526,6 +1542,8 @@ "mu": 1.6908, "kappa_ab": 0.03, "epsilon_k_ab": 2002.9613, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.6684, -3.9471, @@ -1639,6 +1657,8 @@ "mu": 2.4103, "kappa_ab": 0.03, "epsilon_k_ab": 2711.6591, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.7689, -0.9258, @@ -1664,6 +1684,8 @@ "mu": 1.6908, "kappa_ab": 0.03, "epsilon_k_ab": 2514.0609, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.08153, -1.1998, @@ -1870,6 +1892,8 @@ "mu": 1.6608, "kappa_ab": 0.03, "epsilon_k_ab": 2220.2939, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.1698, -3.4917, @@ -1895,6 +1919,8 @@ "mu": 1.7388, "kappa_ab": 0.03, "epsilon_k_ab": 1956.8424, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.3075, -2.24414, @@ -1920,6 +1946,8 @@ "mu": 1.6488, "kappa_ab": 0.03, "epsilon_k_ab": 1878.5922, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.3001, -1.7492, @@ -2080,6 +2108,8 @@ "mu": 1.6698, "kappa_ab": 0.03, "epsilon_k_ab": 2140.2545, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.8993, -3.0227, @@ -2127,6 +2157,8 @@ "mu": 1.5889, "kappa_ab": 0.03, "epsilon_k_ab": 82.5108, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.0803, -2.3348, @@ -2152,6 +2184,8 @@ "mu": 1.6488, "kappa_ab": 0.03, "epsilon_k_ab": 1930.6231, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.218, -2.0634, @@ -2177,6 +2211,8 @@ "mu": 1.6578, "kappa_ab": 0.03, "epsilon_k_ab": 1793.3615, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.2075, -1.7111, @@ -2291,6 +2327,8 @@ "mu": 1.3101, "kappa_ab": 0.03, "epsilon_k_ab": 1103.6088, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.369, -1.2753, @@ -2338,6 +2376,8 @@ "mu": 1.6998, "kappa_ab": 0.03, "epsilon_k_ab": 2519.7116, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.577, -0.44059, @@ -2454,6 +2494,8 @@ "mu": 0.8994, "kappa_ab": 0.03, "epsilon_k_ab": 2731.5672, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.3279, -3.1517, @@ -2479,6 +2521,8 @@ "mu": 1.6099, "kappa_ab": 0.03, "epsilon_k_ab": 2073.3555, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.4657, -2.71167, @@ -2592,6 +2636,8 @@ "mu": 1.6608, "kappa_ab": 0.03, "epsilon_k_ab": 2150.7482, + "na": 1.0, + "nb": 1.0, "viscosity": [ -2.0036, -3.7363, @@ -2639,6 +2685,8 @@ "mu": 1.421, "kappa_ab": 0.03, "epsilon_k_ab": 2476.9879, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.2518, -2.9995, @@ -2664,6 +2712,8 @@ "mu": 1.6488, "kappa_ab": 0.03, "epsilon_k_ab": 1869.8446, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.3896, -2.58867, @@ -2689,6 +2739,8 @@ "mu": 1.6488, "kappa_ab": 0.03, "epsilon_k_ab": 1721.3869, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.3845, -1.1782, @@ -2781,6 +2833,8 @@ "mu": 1.5499, "kappa_ab": 0.03, "epsilon_k_ab": 1923.3059, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.9812, -2.6533, @@ -2806,6 +2860,8 @@ "mu": 1.7, "kappa_ab": 0.03, "epsilon_k_ab": 1835.0315, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.1188, -2.22062, @@ -2831,6 +2887,8 @@ "mu": 1.6668, "kappa_ab": 0.03, "epsilon_k_ab": 1713.2429, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.1043, -1.0389, @@ -2856,6 +2914,8 @@ "mu": 1.6399, "kappa_ab": 0.03, "epsilon_k_ab": 1806.0011, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.1043, -0.9823, @@ -2971,6 +3031,8 @@ "mu": 1.1692, "kappa_ab": 0.03, "epsilon_k_ab": 86.0572, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.7383, -1.7723, @@ -2996,6 +3058,8 @@ "mu": 1.6788, "kappa_ab": 0.03, "epsilon_k_ab": 2044.5298, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.8756, -1.47721, @@ -3021,6 +3085,8 @@ "mu": 1.6608, "kappa_ab": 0.03, "epsilon_k_ab": 1871.8788, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.847, -1.0735, @@ -3181,6 +3247,8 @@ "mu": 1.55, "kappa_ab": 0.03, "epsilon_k_ab": 2238.1838, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.7882, -2.8133, @@ -3294,6 +3362,8 @@ "mu": 1.6488, "kappa_ab": 0.03, "epsilon_k_ab": 2212.4364, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.7295, -2.8003, @@ -3341,6 +3411,8 @@ "mu": 1.6698, "kappa_ab": 0.03, "epsilon_k_ab": 1888.8681, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.6044, -3.44231, diff --git a/parameters/pcsaft/rehner2020.json b/parameters/pcsaft/rehner2020.json index 87d9d42e6..df289c35b 100644 --- a/parameters/pcsaft/rehner2020.json +++ b/parameters/pcsaft/rehner2020.json @@ -14,7 +14,9 @@ "sigma": 2.937523956051823, "epsilon_k": 272.02757407828676, "kappa_ab": 0.04448012165716923, - "epsilon_k_ab": 3125.3202766200056 + "epsilon_k_ab": 3125.3202766200056, + "na": 1.0, + "nb": 1.0 } }, { @@ -33,7 +35,8 @@ "epsilon_k": 238.31794591948915, "kappa_ab": 0.037807044304069184, "epsilon_k_ab": 2749.004567423258, - "na": 2 + "na": 2.0, + "nb": 1.0 } }, { @@ -52,8 +55,8 @@ "epsilon_k": 169.7751872692369, "kappa_ab": 0.13373757842938191, "epsilon_k_ab": 1772.0393059972052, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -72,7 +75,9 @@ "epsilon_k": 166.6147951235982, "mu": 1.6152087869692175, "kappa_ab": 0.09819448826630345, - "epsilon_k_ab": 2667.2518268470913 + "epsilon_k_ab": 2667.2518268470913, + "na": 1.0, + "nb": 1.0 } }, { @@ -92,7 +97,8 @@ "mu": 1.9373715367486852, "kappa_ab": 0.038235772849856756, "epsilon_k_ab": 2377.871413812318, - "na": 2 + "na": 2.0, + "nb": 1.0 } }, { @@ -112,8 +118,8 @@ "mu": 1.5050185176397637, "kappa_ab": 0.082906656369032, "epsilon_k_ab": 1784.1459137848506, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -132,8 +138,8 @@ "epsilon_k": 101.0845669390466, "kappa_ab": 0.11953491759977186, "epsilon_k_ab": 1834.845540500216, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -152,8 +158,8 @@ "epsilon_k": 124.58467252505255, "kappa_ab": 0.10067627754508424, "epsilon_k_ab": 1810.4338396651838, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -172,8 +178,8 @@ "epsilon_k": 144.37733584205245, "kappa_ab": 0.05810650285845608, "epsilon_k_ab": 1959.8294504365563, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -192,8 +198,8 @@ "epsilon_k": 163.6648110273199, "kappa_ab": 0.07216895211562263, "epsilon_k_ab": 1862.0561323915151, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -212,8 +218,8 @@ "epsilon_k": 179.754303082289, "kappa_ab": 0.0762520371571717, "epsilon_k_ab": 1824.6154575116475, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -231,7 +237,9 @@ "sigma": 3.7482827652225086, "epsilon_k": 270.86189939448803, "kappa_ab": 0.002566318648210027, - "epsilon_k_ab": 2778.9261790720657 + "epsilon_k_ab": 2778.9261790720657, + "na": 1.0, + "nb": 1.0 } }, { @@ -249,7 +257,9 @@ "sigma": 3.620930465025628, "epsilon_k": 252.7630330258869, "kappa_ab": 0.003060843493411822, - "epsilon_k_ab": 2700.937487295063 + "epsilon_k_ab": 2700.937487295063, + "na": 1.0, + "nb": 1.0 } }, { @@ -267,7 +277,9 @@ "sigma": 3.8281719743325215, "epsilon_k": 268.72917896983984, "kappa_ab": 0.0031104598923094376, - "epsilon_k_ab": 2804.2762044751535 + "epsilon_k_ab": 2804.2762044751535, + "na": 1.0, + "nb": 1.0 } }, { @@ -285,7 +297,9 @@ "sigma": 3.942380472821992, "epsilon_k": 278.68277806898686, "kappa_ab": 0.0014907242131781597, - "epsilon_k_ab": 3103.6530564736395 + "epsilon_k_ab": 3103.6530564736395, + "na": 1.0, + "nb": 1.0 } }, { @@ -303,7 +317,9 @@ "sigma": 4.042094734564122, "epsilon_k": 281.2912655243266, "kappa_ab": 0.002493036173908781, - "epsilon_k_ab": 3023.609567972825 + "epsilon_k_ab": 3023.609567972825, + "na": 1.0, + "nb": 1.0 } }, { @@ -321,7 +337,9 @@ "sigma": 3.977088533521953, "epsilon_k": 274.5543622578429, "kappa_ab": 0.0012672291257527345, - "epsilon_k_ab": 3223.361822388517 + "epsilon_k_ab": 3223.361822388517, + "na": 1.0, + "nb": 1.0 } }, { @@ -340,8 +358,8 @@ "epsilon_k": 136.16495641117777, "kappa_ab": 0.09762722816926367, "epsilon_k_ab": 1718.8359393009746, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -360,8 +378,8 @@ "epsilon_k": 149.97423748959227, "kappa_ab": 0.16261464675774262, "epsilon_k_ab": 1469.6277485887485, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -380,8 +398,8 @@ "epsilon_k": 159.24901703474268, "kappa_ab": 0.0787091926637352, "epsilon_k_ab": 1771.398524849797, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -400,8 +418,8 @@ "epsilon_k": 162.73433804276743, "kappa_ab": 0.12199296128478557, "epsilon_k_ab": 1516.0966359637514, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -420,8 +438,8 @@ "epsilon_k": 211.90069259634961, "kappa_ab": 0.041207808637705026, "epsilon_k_ab": 2550.2519514335218, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -440,8 +458,8 @@ "epsilon_k": 183.09717953416344, "kappa_ab": 0.08499573407848357, "epsilon_k_ab": 2268.8453113675146, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -460,8 +478,8 @@ "epsilon_k": 299.7826857254255, "kappa_ab": 0.019561736415810844, "epsilon_k_ab": 3069.7332849899476, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } } ] \ No newline at end of file diff --git a/parameters/pcsaft/sauer2014_hetero.json b/parameters/pcsaft/sauer2014_hetero.json index da9755997..fc5485f57 100644 --- a/parameters/pcsaft/sauer2014_hetero.json +++ b/parameters/pcsaft/sauer2014_hetero.json @@ -192,7 +192,9 @@ "sigma": 2.7702, "epsilon_k": 334.29, "epsilon_k_ab": 2575.9, - "kappa_ab": 0.009583 + "kappa_ab": 0.009583, + "na": 1.0, + "nb": 1.0 }, "molarweight": 17.00734 }, @@ -203,7 +205,9 @@ "sigma": 3.1129, "epsilon_k": 309.93, "epsilon_k_ab": 1471.5, - "kappa_ab": 0.005769 + "kappa_ab": 0.005769, + "na": 1.0, + "nb": 1.0 }, "molarweight": 16.02238 } diff --git a/parameters/pcsaft/sauer2014_hetero_joback.json b/parameters/pcsaft/sauer2014_hetero_joback.json index 9198b27bf..2ddb5dd7e 100644 --- a/parameters/pcsaft/sauer2014_hetero_joback.json +++ b/parameters/pcsaft/sauer2014_hetero_joback.json @@ -332,7 +332,9 @@ "sigma": 2.7702, "epsilon_k": 334.29, "epsilon_k_ab": 2575.9, - "kappa_ab": 0.009583 + "kappa_ab": 0.009583, + "na": 1.0, + "nb": 1.0 }, "molarweight": 17.00734, "ideal_gas_record": { @@ -350,7 +352,9 @@ "sigma": 3.1129, "epsilon_k": 309.93, "epsilon_k_ab": 1471.5, - "kappa_ab": 0.005769 + "kappa_ab": 0.005769, + "na": 1.0, + "nb": 1.0 }, "molarweight": 16.02238, "ideal_gas_record": { diff --git a/parameters/pcsaft/sauer2014_homo.json b/parameters/pcsaft/sauer2014_homo.json index 67aa45680..cb92a9fb0 100644 --- a/parameters/pcsaft/sauer2014_homo.json +++ b/parameters/pcsaft/sauer2014_homo.json @@ -17,7 +17,6 @@ }, "molarweight": 14.02658 }, - { "identifier": ">CH", "model_record": { @@ -57,7 +56,7 @@ { "identifier": "=C<", "model_record": { - "m": 0.86367, + "m": 0.86367, "sigma": 3.1815, "epsilon_k": 156.31 }, @@ -66,7 +65,7 @@ { "identifier": "C≡CH", "model_record": { - "m": 1.3279, + "m": 1.3279, "sigma": 2.9421, "epsilon_k": 223.05 }, @@ -75,7 +74,7 @@ { "identifier": "CH2_hex", "model_record": { - "m": 0.39496, + "m": 0.39496, "sigma": 3.9126, "epsilon_k": 289.03 }, @@ -84,7 +83,7 @@ { "identifier": "CH_hex", "model_record": { - "m": 0.02880, + "m": 0.02880, "sigma": 8.9779, "epsilon_k": 1306.7 }, @@ -93,7 +92,7 @@ { "identifier": "CH2_pent", "model_record": { - "m": 0.46742, + "m": 0.46742, "sigma": 3.7272, "epsilon_k": 267.16 }, @@ -102,7 +101,7 @@ { "identifier": "CH_pent", "model_record": { - "m": 0.03314, + "m": 0.03314, "sigma": 7.7190, "epsilon_k": 1297.7 }, @@ -111,7 +110,7 @@ { "identifier": "CH_arom", "model_record": { - "m": 0.42335, + "m": 0.42335, "sigma": 3.7270, "epsilon_k": 274.41 }, @@ -169,7 +168,7 @@ { "identifier": "HCOO", "model_record": { - "m": 1.7525, + "m": 1.7525, "sigma": 2.9043, "epsilon_k": 229.63, "mu": 2.7916 @@ -193,7 +192,9 @@ "sigma": 3.2859, "epsilon_k": 488.66, "epsilon_k_ab": 2517.0, - "kappa_ab": 0.006825 + "kappa_ab": 0.006825, + "na": 1.0, + "nb": 1.0 }, "molarweight": 17.00734 }, @@ -204,9 +205,10 @@ "sigma": 3.6456, "epsilon_k": 467.59, "epsilon_k_ab": 1064.6, - "kappa_ab": 0.026662 + "kappa_ab": 0.026662, + "na": 1.0, + "nb": 1.0 }, "molarweight": 16.02238 } -] - +] \ No newline at end of file diff --git a/parameters/pcsaft/sauer2014_homo_joback.json b/parameters/pcsaft/sauer2014_homo_joback.json index 534c9ae75..c56e8632b 100644 --- a/parameters/pcsaft/sauer2014_homo_joback.json +++ b/parameters/pcsaft/sauer2014_homo_joback.json @@ -332,7 +332,9 @@ "sigma": 3.2859, "epsilon_k": 488.66, "epsilon_k_ab": 2517.0, - "kappa_ab": 0.006825 + "kappa_ab": 0.006825, + "na": 1.0, + "nb": 1.0 }, "molarweight": 17.00734, "ideal_gas_record": { @@ -350,7 +352,9 @@ "sigma": 3.6456, "epsilon_k": 467.59, "epsilon_k_ab": 1064.6, - "kappa_ab": 0.026662 + "kappa_ab": 0.026662, + "na": 1.0, + "nb": 1.0 }, "molarweight": 16.02238, "ideal_gas_record": { diff --git a/src/association/dft.rs b/src/association/dft.rs index e824efb24..05c92b378 100644 --- a/src/association/dft.rs +++ b/src/association/dft.rs @@ -52,9 +52,6 @@ where // number of segments let n = self.association_parameters.component_index.len(); - // number of associating segments - let nassoc = self.association_parameters.assoc_comp.len(); - // number of dimensions let dim = (weighted_densities.shape()[0] - 1) / n - 1; @@ -70,7 +67,7 @@ where .collect(); let n3 = weighted_densities.index_axis(Axis(0), n * (dim + 1)); - // calculate rho0 (only associating segments) + // calculate rho0 let [_, _, c2, _] = p.geometry_coefficients(temperature); let diameter = p.hs_diameter(temperature); let mut n2i = n0i.to_owned(); @@ -88,9 +85,6 @@ where *rho0 = n0i; } }); - let rho0 = Array2::from_shape_fn((nassoc, n3.len()), |(i, j)| { - rho0[(self.association_parameters.assoc_comp[i], j)] - }); // calculate xi let n2v: Vec<_> = n2vi.iter().map(|n2vi| n2vi.sum_axis(Axis(0))).collect(); @@ -111,56 +105,143 @@ where // auxiliary variables let n3i = n3.mapv(|n3| (-n3 + 1.0).recip()); - // only one associating component - if nassoc == 1 { - // association strength - let k = &n2 * &n3i * diameter[self.association_parameters.assoc_comp[0]] * 0.5; - let deltarho = (((&k / 18.0 + 0.5) * &k * xi + 1.0) * n3i) - * ((temperature.recip() * self.association_parameters.epsilon_k_aibj[(0, 0)]) - .exp_m1() - * self.association_parameters.sigma3_kappa_aibj[(0, 0)]) - * rho0.index_axis(Axis(0), 0); - - let na = self.association_parameters.na[0]; - let nb = self.association_parameters.nb[0]; - let f = |x: N| x.ln() - x * 0.5 + 0.5; - if nb > 0.0 { - // no cross association, two association sites - let xa = deltarho.mapv(|d| Self::assoc_site_frac_ab(d, na, nb)); - let xb = (&xa - 1.0) * (na / nb) + 1.0; - Ok((xa.mapv(f) * na + xb.mapv(f) * nb) * rho0.index_axis(Axis(0), 0)) - } else { - // no cross association, one association site - let xa = deltarho.mapv(|d| Self::assoc_site_frac_a(d, na)); - - Ok(xa.mapv(f) * na * rho0.index_axis(Axis(0), 0)) + self.calculate_helmholtz_energy_density(temperature, &rho0, &n2, &n3i, &xi) + } +} + +impl Association

{ + pub fn calculate_helmholtz_energy_density< + N: DualNum + ScalarOperand, + S: Data, + >( + &self, + temperature: N, + rho0: &Array2, + n2: &ArrayBase, + n3i: &Array1, + xi: &Array1, + ) -> EosResult> { + let a = &self.association_parameters; + + let d = self.parameters.hs_diameter(temperature); + + match ( + a.sites_a.len() * a.sites_b.len(), + a.sites_c.len(), + self.force_cross_association, + ) { + (0, 0, _) => Ok(Array::zeros(n3i.len())), + (1, 0, false) => { + Ok(self.helmholtz_energy_density_ab_analytic(temperature, rho0, &d, n2, n3i, xi)) + } + (0, 1, false) => { + Ok(self.helmholtz_energy_density_cc_analytic(temperature, rho0, &d, n2, n3i, xi)) + } + (1, 1, false) => { + Ok( + self.helmholtz_energy_density_ab_analytic(temperature, rho0, &d, n2, n3i, xi) + + self.helmholtz_energy_density_cc_analytic( + temperature, + rho0, + &d, + n2, + n3i, + xi, + ), + ) + } + _ => { + let mut x: Array1 = Array::from_elem(a.sites_a.len() + a.sites_b.len(), 0.2); + let (assoc_comp_ab, n_ab): (Vec<_>, Vec<_>) = a + .sites_a + .iter() + .chain(a.sites_b.iter()) + .map(|s| (s.assoc_comp, s.n)) + .unzip(); + let rhoab = Array2::from_shape_fn((x.len(), n3i.len()), |(i, j)| { + rho0[(assoc_comp_ab[i], j)] * n_ab[i] + }); + rhoab + .axis_iter(Axis(1)) + .zip(n2.iter()) + .zip(n3i.iter()) + .zip(xi.iter()) + .map(|(((rho, &n2), &n3i), &xi)| { + let [delta_ab, delta_cc] = + self.association_strength(temperature, &d, n2, n3i, xi); + Self::helmholtz_energy_density_cross_association( + &rho, + &delta_ab, + &delta_cc, + self.max_iter, + self.tol, + Some(&mut x), + ) + }) + .collect() } - } else { - let mut x: Array1 = Array::from_elem(2 * nassoc, 0.2); - Ok(rho0 - .view() - .into_shape([nassoc, rho0.len() / nassoc]) - .unwrap() - .axis_iter(Axis(1)) - .zip(n2.iter()) - .zip(n3i.iter()) - .zip(xi.iter()) - .map(|(((rho0, &n2), &n3i), &xi)| { - self.helmholtz_energy_density_cross_association( - temperature, - &rho0, - &diameter, - n2, - n3i, - xi, - self.max_iter, - self.tol, - Some(&mut x), - ) - }) - .collect::, _>>()? - .into_shape(n2.raw_dim()) - .unwrap()) } } + + fn helmholtz_energy_density_ab_analytic + ScalarOperand, S: Data>( + &self, + temperature: N, + rho0: &Array2, + diameter: &Array1, + n2: &ArrayBase, + n3i: &Array1, + xi: &Array1, + ) -> Array1 { + let a = &self.association_parameters; + + // site densities + let i = a.sites_a[0].assoc_comp; + let j = a.sites_b[0].assoc_comp; + let rhoa = &rho0.index_axis(Axis(0), i) * a.sites_a[0].n; + let rhob = &rho0.index_axis(Axis(0), j) * a.sites_b[0].n; + + // association strength + let di = diameter[i]; + let dj = diameter[j]; + let k = n2 * n3i * (di * dj / (di + dj)); + let delta = (((&k / 18.0 + 0.5) * &k * xi + 1.0) * n3i) + * ((temperature.recip() * a.epsilon_k_ab[(0, 0)]).exp_m1() * a.sigma3_kappa_ab[(0, 0)]); + + // no cross association, two association sites + let aux = &delta * (&rhob - &rhoa) + 1.0; + let xa = ((&aux * &aux + &delta * &rhoa * 4.0).map(N::sqrt) + &aux).map(N::recip) * 2.0; + let aux = -aux + 2.0; + let xb = ((&aux * &aux + delta * &rhob * 4.0).map(N::sqrt) + aux).map(N::recip) * 2.0; + + let f = |x: N| x.ln() - x * 0.5 + 0.5; + rhoa * xa.mapv(f) + rhob * xb.mapv(f) + } + + fn helmholtz_energy_density_cc_analytic + ScalarOperand, S: Data>( + &self, + temperature: N, + rho0: &Array2, + diameter: &Array1, + n2: &ArrayBase, + n3i: &Array1, + xi: &Array1, + ) -> Array1 { + let a = &self.association_parameters; + + // site densities + let i = a.sites_c[0].assoc_comp; + let rhoc = &rho0.index_axis(Axis(0), i) * a.sites_c[0].n; + + // association strength + let di = diameter[i]; + let k = n2 * n3i * (di * 0.5); + let delta = (((&k / 18.0 + 0.5) * &k * xi + 1.0) * n3i) + * ((temperature.recip() * a.epsilon_k_cc[(0, 0)]).exp_m1() * a.sigma3_kappa_cc[(0, 0)]); + + // no cross association, two association sites + let xc = ((delta * 4.0 * &rhoc + 1.0).map(N::sqrt) + 1.0).map(N::recip) * 2.0; + + let f = |x: N| x.ln() - x * 0.5 + 0.5; + rhoc * xc.mapv(f) + } } diff --git a/src/association/mod.rs b/src/association/mod.rs index aa1b74bb7..8dfcbc9de 100644 --- a/src/association/mod.rs +++ b/src/association/mod.rs @@ -1,13 +1,13 @@ //! Generic implementation of the SAFT association contribution //! that can be used across models. use crate::hard_sphere::HardSphereProperties; -use feos_core::{EosError, HelmholtzEnergyDual, StateHD}; +use feos_core::{EosError, EosResult, HelmholtzEnergyDual, StateHD}; use ndarray::*; use num_dual::linalg::{norm, LU}; use num_dual::*; +use num_traits::Zero; use serde::{Deserialize, Serialize}; use std::fmt; -use std::ops::SubAssign; use std::sync::Arc; #[cfg(feature = "dft")] @@ -17,6 +17,25 @@ mod python; #[cfg(feature = "python")] pub use python::PyAssociationRecord; +#[derive(Clone, Copy, Debug)] +struct AssociationSite { + assoc_comp: usize, + n: f64, + kappa_ab: f64, + epsilon_k_ab: f64, +} + +impl AssociationSite { + fn new(assoc_comp: usize, n: f64, kappa_ab: f64, epsilon_k_ab: f64) -> Self { + Self { + assoc_comp, + n, + kappa_ab, + epsilon_k_ab, + } + } +} + /// Pure component association parameters. #[derive(Serialize, Deserialize, Clone, Copy, Default)] pub struct AssociationRecord { @@ -25,20 +44,27 @@ pub struct AssociationRecord { /// Association energy parameter in units of Kelvin pub epsilon_k_ab: f64, /// \# of association sites of type A - #[serde(skip_serializing_if = "Option::is_none")] - pub na: Option, + #[serde(skip_serializing_if = "f64::is_zero")] + #[serde(default)] + pub na: f64, /// \# of association sites of type B - #[serde(skip_serializing_if = "Option::is_none")] - pub nb: Option, + #[serde(skip_serializing_if = "f64::is_zero")] + #[serde(default)] + pub nb: f64, + /// \# of association sites of type C + #[serde(skip_serializing_if = "f64::is_zero")] + #[serde(default)] + pub nc: f64, } impl AssociationRecord { - pub fn new(kappa_ab: f64, epsilon_k_ab: f64, na: Option, nb: Option) -> Self { + pub fn new(kappa_ab: f64, epsilon_k_ab: f64, na: f64, nb: f64, nc: f64) -> Self { Self { kappa_ab, epsilon_k_ab, na, nb, + nc, } } } @@ -47,8 +73,16 @@ impl fmt::Display for AssociationRecord { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { write!(f, "AssociationRecord(kappa_ab={}", self.kappa_ab)?; write!(f, ", epsilon_k_ab={}", self.epsilon_k_ab)?; - write!(f, ", na={}", self.na.unwrap_or(1.0))?; - write!(f, ", nb={})", self.nb.unwrap_or(1.0)) + if self.na > 0.0 { + write!(f, ", na={}", self.na)?; + } + if self.nb > 0.0 { + write!(f, ", nb={}", self.nb)?; + } + if self.nc > 0.0 { + write!(f, ", nc={}", self.nc)?; + } + write!(f, ")") } } @@ -57,61 +91,88 @@ impl fmt::Display for AssociationRecord { #[derive(Clone)] pub struct AssociationParameters { component_index: Array1, - pub assoc_comp: Array1, - pub kappa_ab: Array1, - pub epsilon_k_ab: Array1, - pub sigma3_kappa_aibj: Array2, - pub epsilon_k_aibj: Array2, - pub na: Array1, - pub nb: Array1, + sites_a: Array1, + sites_b: Array1, + sites_c: Array1, + pub sigma3_kappa_ab: Array2, + pub sigma3_kappa_cc: Array2, + pub epsilon_k_ab: Array2, + pub epsilon_k_cc: Array2, } impl AssociationParameters { pub fn new( - records: &[Option], + records: &[Vec], sigma: &Array1, component_index: Option<&Array1>, ) -> Self { - let mut assoc_comp = Vec::new(); - let mut sigma_assoc = Vec::new(); - let mut kappa_ab = Vec::new(); - let mut epsilon_k_ab = Vec::new(); - let mut na = Vec::new(); - let mut nb = Vec::new(); + let mut sites_a = Vec::new(); + let mut sites_b = Vec::new(); + let mut sites_c = Vec::new(); for (i, record) in records.iter().enumerate() { - if let Some(record) = record.as_ref() { - if record.kappa_ab > 0.0 && record.epsilon_k_ab > 0.0 { - assoc_comp.push(i); - sigma_assoc.push(sigma[i]); - kappa_ab.push(record.kappa_ab); - epsilon_k_ab.push(record.epsilon_k_ab); - na.push(record.na.unwrap_or(1.0)); - nb.push(record.nb.unwrap_or(1.0)); + for site in record { + if site.kappa_ab > 0.0 && site.epsilon_k_ab > 0.0 { + if site.na > 0.0 { + sites_a.push(AssociationSite::new( + i, + site.na, + site.kappa_ab, + site.epsilon_k_ab, + )); + } + if site.nb > 0.0 { + sites_b.push(AssociationSite::new( + i, + site.nb, + site.kappa_ab, + site.epsilon_k_ab, + )); + } + if site.nc > 0.0 { + sites_c.push(AssociationSite::new( + i, + site.nc, + site.kappa_ab, + site.epsilon_k_ab, + )); + } } } } - let sigma3_kappa_aibj = Array2::from_shape_fn([kappa_ab.len(); 2], |(i, j)| { - (sigma_assoc[i] * sigma_assoc[j]).powf(1.5) * (kappa_ab[i] * kappa_ab[j]).sqrt() + let sigma3_kappa_ab = Array2::from_shape_fn([sites_a.len(), sites_b.len()], |(i, j)| { + (sigma[sites_a[i].assoc_comp] * sigma[sites_b[j].assoc_comp]).powf(1.5) + * (sites_a[i].kappa_ab * sites_b[j].kappa_ab).sqrt() + }); + let sigma3_kappa_cc = Array2::from_shape_fn([sites_c.len(); 2], |(i, j)| { + (sigma[sites_c[i].assoc_comp] * sigma[sites_c[j].assoc_comp]).powf(1.5) + * (sites_c[i].kappa_ab * sites_c[j].kappa_ab).sqrt() + }); + let epsilon_k_ab = Array2::from_shape_fn([sites_a.len(), sites_b.len()], |(i, j)| { + 0.5 * (sites_a[i].epsilon_k_ab + sites_b[j].epsilon_k_ab) }); - let epsilon_k_aibj = Array2::from_shape_fn([epsilon_k_ab.len(); 2], |(i, j)| { - 0.5 * (epsilon_k_ab[i] + epsilon_k_ab[j]) + let epsilon_k_cc = Array2::from_shape_fn([sites_c.len(); 2], |(i, j)| { + 0.5 * (sites_c[i].epsilon_k_ab + sites_c[j].epsilon_k_ab) }); Self { component_index: component_index .cloned() .unwrap_or_else(|| Array1::from_shape_fn(records.len(), |i| i)), - assoc_comp: Array1::from_vec(assoc_comp), - kappa_ab: Array1::from_vec(kappa_ab), - epsilon_k_ab: Array1::from_vec(epsilon_k_ab), - sigma3_kappa_aibj, - epsilon_k_aibj, - na: Array1::from_vec(na), - nb: Array1::from_vec(nb), + sites_a: Array1::from_vec(sites_a), + sites_b: Array1::from_vec(sites_b), + sites_c: Array1::from_vec(sites_c), + sigma3_kappa_ab, + sigma3_kappa_cc, + epsilon_k_ab, + epsilon_k_cc, } } + + pub fn is_empty(&self) -> bool { + (self.sites_a.is_empty() | self.sites_b.is_empty()) & self.sites_c.is_empty() + } } /// Implementation of the SAFT association Helmholtz energy @@ -158,17 +219,25 @@ impl Association

{ n2: D, n3i: D, xi: D, - ) -> Array2 { - // Calculate association strength - let ac = &self.association_parameters.assoc_comp; - Array2::from_shape_fn([ac.len(); 2], |(i, j)| { - let k = diameter[ac[i]] * diameter[ac[j]] / (diameter[ac[i]] + diameter[ac[j]]) - * (n2 * n3i); + ) -> [Array2; 2] { + let p = &self.association_parameters; + let delta_ab = Array2::from_shape_fn([p.sites_a.len(), p.sites_b.len()], |(i, j)| { + let di = diameter[p.sites_a[i].assoc_comp]; + let dj = diameter[p.sites_b[j].assoc_comp]; + let k = di * dj / (di + dj) * (n2 * n3i); n3i * (k * xi * (k / 18.0 + 0.5) + 1.0) - * self.association_parameters.sigma3_kappa_aibj[(i, j)] - * (temperature.recip() * self.association_parameters.epsilon_k_aibj[(i, j)]) - .exp_m1() - }) + * p.sigma3_kappa_ab[(i, j)] + * (temperature.recip() * p.epsilon_k_ab[(i, j)]).exp_m1() + }); + let delta_cc = Array2::from_shape_fn([p.sites_c.len(); 2], |(i, j)| { + let di = diameter[p.sites_c[i].assoc_comp]; + let dj = diameter[p.sites_c[j].assoc_comp]; + let k = di * dj / (di + dj) * (n2 * n3i); + n3i * (k * xi * (k / 18.0 + 0.5) + 1.0) + * p.sigma3_kappa_cc[(i, j)] + * (temperature.recip() * p.epsilon_k_cc[(i, j)]).exp_m1() + }); + [delta_ab, delta_cc] } } @@ -177,6 +246,7 @@ impl + ScalarOperand, P: HardSphereProperties> HelmholtzEnergyDu { fn helmholtz_energy(&self, state: &StateHD) -> D { let p: &P = &self.parameters; + let a = &self.association_parameters; // temperature dependent segment diameter let diameter = p.hs_diameter(state.temperature); @@ -186,48 +256,43 @@ impl + ScalarOperand, P: HardSphereProperties> HelmholtzEnergyDu let n2 = zeta2 * 6.0; let n3i = (-n3 + 1.0).recip(); - if self.association_parameters.assoc_comp.len() > 1 || self.force_cross_association { - // extract densities of associating segments - let rho_assoc = self - .association_parameters - .assoc_comp - .mapv(|a| state.partial_density[self.association_parameters.component_index[a]]); - - // Helmholtz energy - self.helmholtz_energy_density_cross_association( - state.temperature, - &rho_assoc, - &diameter, - n2, - n3i, - D::one(), - self.max_iter, - self.tol, - None, - ) - .unwrap_or_else(|_| D::from(std::f64::NAN)) - * state.volume - } else { - // association strength - let c = self.association_parameters.component_index - [self.association_parameters.assoc_comp[0]]; - let deltarho = - self.association_strength(state.temperature, &diameter, n2, n3i, D::one())[(0, 0)] - * state.partial_density[c]; - - let na = self.association_parameters.na[0]; - let nb = self.association_parameters.nb[0]; - if nb > 0.0 { - // no cross association, two association sites - let xa = Self::assoc_site_frac_ab(deltarho, na, nb); - let xb = (xa - 1.0) * (na / nb) + 1.0; - - state.moles[c] * ((xa.ln() - xa * 0.5 + 0.5) * na + (xb.ln() - xb * 0.5 + 0.5) * nb) - } else { - // no cross association, one association site - let xa = Self::assoc_site_frac_a(deltarho, na); - - state.moles[c] * (xa.ln() - xa * 0.5 + 0.5) * na + // association strength + let [delta_ab, delta_cc] = + self.association_strength(state.temperature, &diameter, n2, n3i, D::one()); + + match ( + a.sites_a.len() * a.sites_b.len(), + a.sites_c.len(), + self.force_cross_association, + ) { + (0, 0, _) => D::zero(), + (1, 0, false) => self.helmholtz_energy_ab_analytic(state, delta_ab[(0, 0)]), + (0, 1, false) => self.helmholtz_energy_cc_analytic(state, delta_cc[(0, 0)]), + (1, 1, false) => { + self.helmholtz_energy_ab_analytic(state, delta_ab[(0, 0)]) + + self.helmholtz_energy_cc_analytic(state, delta_cc[(0, 0)]) + } + _ => { + // extract site densities of associating segments + let rho: Array1<_> = a + .sites_a + .iter() + .chain(a.sites_b.iter()) + .chain(a.sites_c.iter()) + .map(|s| state.partial_density[a.component_index[s.assoc_comp]] * s.n) + .collect(); + + // Helmholtz energy + Self::helmholtz_energy_density_cross_association( + &rho, + &delta_ab, + &delta_cc, + self.max_iter, + self.tol, + None, + ) + .unwrap_or_else(|_| D::from(std::f64::NAN)) + * state.volume } } } @@ -240,66 +305,72 @@ impl

fmt::Display for Association

{ } impl Association

{ - pub fn assoc_site_frac_ab>(deltarho: D, na: f64, nb: f64) -> D { - (((deltarho * (na - nb) + 1.0).powi(2) + deltarho * nb * 4.0).sqrt() - + (deltarho * (nb - na) + 1.0)) - .recip() - * 2.0 + fn helmholtz_energy_ab_analytic>(&self, state: &StateHD, delta: D) -> D { + let a = &self.association_parameters; + + // site densities + let rhoa = + state.partial_density[a.component_index[a.sites_a[0].assoc_comp]] * a.sites_a[0].n; + let rhob = + state.partial_density[a.component_index[a.sites_b[0].assoc_comp]] * a.sites_b[0].n; + + // fraction of non-bonded association sites + let sqrt = ((delta * (rhoa - rhob) + 1.0).powi(2) + delta * rhob * 4.0).sqrt(); + let xa = (sqrt + (delta * (rhob - rhoa) + 1.0)).recip() * 2.0; + let xb = (sqrt + (delta * (rhoa - rhob) + 1.0)).recip() * 2.0; + + (rhoa * (xa.ln() - xa * 0.5 + 0.5) + rhob * (xb.ln() - xb * 0.5 + 0.5)) * state.volume } - pub fn assoc_site_frac_a>(deltarho: D, na: f64) -> D { - ((deltarho * 4.0 * na + 1.0).sqrt() + 1.0).recip() * 2.0 + fn helmholtz_energy_cc_analytic>(&self, state: &StateHD, delta: D) -> D { + let a = &self.association_parameters; + + // site density + let rhoc = + state.partial_density[a.component_index[a.sites_c[0].assoc_comp]] * a.sites_c[0].n; + + // fraction of non-bonded association sites + let xc = ((delta * 4.0 * rhoc + 1.0).sqrt() + 1.0).recip() * 2.0; + + rhoc * (xc.ln() - xc * 0.5 + 0.5) * state.volume } #[allow(clippy::too_many_arguments)] fn helmholtz_energy_density_cross_association< - S: Data, D: DualNum + ScalarOperand, + S: Data, >( - &self, - temperature: D, - density: &ArrayBase, - diameter: &Array1, - n2: D, - n3i: D, - xi: D, + rho: &ArrayBase, + delta_ab: &Array2, + delta_cc: &Array2, max_iter: usize, tol: f64, x0: Option<&mut Array1>, - ) -> Result { + ) -> EosResult { // check if density is close to 0 - if density.sum().re() < f64::EPSILON { + if rho.sum().re() < f64::EPSILON { if let Some(x0) = x0 { x0.fill(1.0); } return Ok(D::zero()); } - let assoc_comp = &self.association_parameters.assoc_comp; - let nassoc = assoc_comp.len(); - - // association strength - let delta = self.association_strength(temperature, diameter, n2, n3i, xi); - - // extract parameters of associating components - let na = &self.association_parameters.na; - let nb = &self.association_parameters.nb; - // cross-association according to Michelsen2006 // initialize monomer fraction let mut x = match &x0 { Some(x0) => (*x0).clone(), - None => Array::from_elem(2 * nassoc, 0.2), + None => Array::from_elem(rho.len(), 0.2), }; + let delta_ab_re = delta_ab.map(D::re); + let delta_cc_re = delta_cc.map(D::re); + let rho_re = rho.map(D::re); for k in 0..max_iter { - if Self::newton_step_cross_association::<_, f64>( - nassoc, + if Self::newton_step_cross_association( &mut x, - &delta.map(D::re), - na, - nb, - &density.map(D::re), + &delta_ab_re, + &delta_cc_re, + &rho_re, tol, )? { break; @@ -312,7 +383,7 @@ impl Association

{ // calculate derivatives let mut x_dual = x.mapv(D::from); for _ in 0..D::NDERIV { - Self::newton_step_cross_association(nassoc, &mut x_dual, &delta, na, nb, density, tol)?; + Self::newton_step_cross_association(&mut x_dual, delta_ab, delta_cc, rho, tol)?; } // save monomer fraction @@ -321,48 +392,72 @@ impl Association

{ } // Helmholtz energy density - let xa = x_dual.slice(s![..nassoc]); - let xb = x_dual.slice(s![nassoc..]); let f = |x: D| x.ln() - x * 0.5 + 0.5; - Ok((density * (xa.mapv(f) * na + xb.mapv(f) * nb)).sum()) + Ok((rho * x_dual.mapv(f)).sum()) } - fn newton_step_cross_association, D: DualNum + ScalarOperand>( - nassoc: usize, + fn newton_step_cross_association + ScalarOperand, S: Data>( x: &mut Array1, - delta: &Array2, - na: &Array1, - nb: &Array1, + delta_ab: &Array2, + delta_cc: &Array2, rho: &ArrayBase, tol: f64, - ) -> Result { + ) -> EosResult { + let nassoc = x.len(); // gradient let mut g = x.map(D::recip); // Hessian - let mut h: Array2 = Array::zeros((2 * nassoc, 2 * nassoc)); + let mut h: Array2 = Array::zeros([nassoc; 2]); - // split x array - let (xa, xb) = x.view().split_at(Axis(0), nassoc); + // split arrays + let &[a, b] = delta_ab.shape() else { panic!("wrong shape!") }; + let c = delta_cc.shape()[0]; + let (xa, xc) = x.view().split_at(Axis(0), a + b); + let (xa, xb) = xa.split_at(Axis(0), a); + let (rhoa, rhoc) = rho.view().split_at(Axis(0), a + b); + let (rhoa, rhob) = rhoa.split_at(Axis(0), a); - // calculate gradients and approximate Hessian for i in 0..nassoc { - let d = &delta.index_axis(Axis(0), i) * rho; - - let dnx = (&xb * nb * &d).sum() + 1.0; + // calculate gradients + let (d, dnx) = if i < a { + let d = delta_ab.index_axis(Axis(0), i); + (d, (&xb * &rhob * d).sum() + 1.0) + } else if i < a + b { + let d = delta_ab.index_axis(Axis(1), i - a); + (d, (&xa * &rhoa * d).sum() + 1.0) + } else { + let d = delta_cc.index_axis(Axis(0), i - a - b); + (d, (&xc * &rhoc * d).sum() + 1.0) + }; g[i] -= dnx; - for j in 0..nassoc { - h[(i, nassoc + j)] = -d[j] * nb[j]; - h[(nassoc + i, j)] = -d[j] * na[j]; - } - h[(i, i)] = -dnx / xa[i]; - let dnx = (&xa * na * &d).sum() + 1.0; - g[nassoc + i] -= dnx; - h[(nassoc + i, nassoc + i)] = -dnx / xb[i]; + // approximate hessian + h[(i, i)] = -dnx / x[i]; + if i < a { + for j in 0..b { + h[(i, a + j)] = -d[j] * rhob[j]; + } + } else if i < a + b { + for j in 0..a { + h[(i, j)] = -d[j] * rhoa[j]; + } + } else { + for j in 0..c { + h[(i, a + b + j)] -= d[j] * rhoc[j]; + } + } } // Newton step - x.sub_assign(&LU::new(h)?.solve(&g)); + // avoid stepping to negative values for x (see Michelsen 2006) + let delta_x = LU::new(h)?.solve(&g); + Zip::from(x).and(&delta_x).for_each(|x, &delta_x| { + if delta_x.re() < x.re() * 0.8 { + *x -= delta_x + } else { + *x *= 0.2 + } + }); // check convergence Ok(norm(&g.map(D::re)) < tol) @@ -407,7 +502,7 @@ mod tests_pcsaft { let mut params = water_parameters(); let mut record = params.pure_records.pop().unwrap(); let mut association_record = record.model_record.association_record.unwrap(); - association_record.na = Some(2.0); + association_record.na = 2.0; record.model_record.association_record = Some(association_record); let params = Arc::new(PcSaftParameters::new_pure(record)); let assoc = Association::new(¶ms, ¶ms.association, 50, 1e-10); diff --git a/src/association/python.rs b/src/association/python.rs index cbc8ed680..9744f5552 100644 --- a/src/association/python.rs +++ b/src/association/python.rs @@ -4,19 +4,16 @@ use feos_core::parameter::ParameterError; use pyo3::prelude::*; /// Pure component association parameters -#[pyclass( - name = "AssociationRecord", - text_signature = "(kappa_ab, epsilon_k_ab, na=None, nb=None)" -)] +#[pyclass(name = "AssociationRecord")] #[derive(Clone)] pub struct PyAssociationRecord(pub AssociationRecord); #[pymethods] impl PyAssociationRecord { - #[pyo3(signature = (kappa_ab, epsilon_k_ab, na=None, nb=None))] #[new] - fn new(kappa_ab: f64, epsilon_k_ab: f64, na: Option, nb: Option) -> Self { - Self(AssociationRecord::new(kappa_ab, epsilon_k_ab, na, nb)) + #[pyo3(signature = (kappa_ab, epsilon_k_ab, na=0.0, nb=0.0, nc=0.0))] + fn new(kappa_ab: f64, epsilon_k_ab: f64, na: f64, nb: f64, nc: f64) -> Self { + Self(AssociationRecord::new(kappa_ab, epsilon_k_ab, na, nb, nc)) } #[getter] @@ -30,15 +27,20 @@ impl PyAssociationRecord { } #[getter] - fn get_na(&self) -> Option { + fn get_na(&self) -> f64 { self.0.na } #[getter] - fn get_nb(&self) -> Option { + fn get_nb(&self) -> f64 { self.0.nb } + #[getter] + fn get_nc(&self) -> f64 { + self.0.nc + } + fn __repr__(&self) -> PyResult { Ok(self.0.to_string()) } diff --git a/src/gc_pcsaft/dft/mod.rs b/src/gc_pcsaft/dft/mod.rs index 438af9ba2..e8e502d65 100644 --- a/src/gc_pcsaft/dft/mod.rs +++ b/src/gc_pcsaft/dft/mod.rs @@ -56,7 +56,7 @@ impl GcPcSaftFunctional { contributions.push(Box::new(att)); // Association - if !parameters.association.assoc_comp.is_empty() { + if !parameters.association.is_empty() { let assoc = Association::new( ¶meters, ¶meters.association, diff --git a/src/gc_pcsaft/dft/parameter.rs b/src/gc_pcsaft/dft/parameter.rs index d6c27d9df..d0931dfac 100644 --- a/src/gc_pcsaft/dft/parameter.rs +++ b/src/gc_pcsaft/dft/parameter.rs @@ -80,7 +80,13 @@ impl ParameterHetero for GcPcSaftFunctionalParameters { sigma.push(segment.model_record.sigma); epsilon_k.push(segment.model_record.epsilon_k); - association_records.push(segment.model_record.association_record); + association_records.push( + segment + .model_record + .association_record + .into_iter() + .collect(), + ); psi_dft.push(segment.model_record.psi_dft.unwrap_or(PSI_GC_DFT)); diff --git a/src/gc_pcsaft/eos/mod.rs b/src/gc_pcsaft/eos/mod.rs index c6e0f9906..eec50fea4 100644 --- a/src/gc_pcsaft/eos/mod.rs +++ b/src/gc_pcsaft/eos/mod.rs @@ -60,7 +60,7 @@ impl GcPcSaft { contributions.push(Box::new(Dispersion { parameters: parameters.clone(), })); - if !parameters.association.assoc_comp.is_empty() { + if !parameters.association.is_empty() { contributions.push(Box::new(Association::new( ¶meters, ¶meters.association, diff --git a/src/gc_pcsaft/eos/parameter.rs b/src/gc_pcsaft/eos/parameter.rs index c25add9ad..795ea2792 100644 --- a/src/gc_pcsaft/eos/parameter.rs +++ b/src/gc_pcsaft/eos/parameter.rs @@ -140,10 +140,10 @@ impl ParameterHetero for GcPcSaftEosParameters { let mut assoc = segment.model_record.association_record; if let Some(mut assoc) = assoc.as_mut() { - assoc.na = Some(assoc.na.unwrap_or(1.0) * count); - assoc.nb = Some(assoc.nb.unwrap_or(1.0) * count); + assoc.na *= count; + assoc.nb *= count; }; - association_records.push(assoc); + association_records.push(assoc.into_iter().collect()); m_i += segment.model_record.m * count; sigma_i += segment.model_record.m * segment.model_record.sigma.powi(3) * count; @@ -287,11 +287,17 @@ impl HardSphereProperties for GcPcSaftEosParameters { impl GcPcSaftEosParameters { pub fn to_markdown(&self) -> String { + let gorup_dict: HashMap<&String, &GcPcSaftRecord> = self + .segment_records + .iter() + .map(|r| (&r.identifier, &r.model_record)) + .collect(); + let mut output = String::new(); let o = &mut output; write!( o, - "|component|molarweight|dipole moment|group|$m$|$\\sigma$|$\\varepsilon$|$\\kappa_{{AB}}$|$\\varepsilon_{{AB}}$|$N_A$|$N_B$|\n|-|-|-|-|-|-|-|-|-|-|-|" + "|component|molarweight|dipole moment|group|$m$|$\\sigma$|$\\varepsilon$|$\\kappa_{{AB}}$|$\\varepsilon_{{AB}}$|$N_A$|$N_B$|$N_C$|\n|-|-|-|-|-|-|-|-|-|-|-|-|" ) .unwrap(); for i in 0..self.m.len() { @@ -305,37 +311,39 @@ impl GcPcSaftEosParameters { .as_ref() .unwrap_or(&format!("Component {}", self.component_index[i] + 1)), self.molarweight[self.component_index[i]], - if let Some(d) = self.dipole_comp.iter().position(|&d| d == i) { + if let Some(d) = self + .dipole_comp + .iter() + .position(|&d| d == self.component_index[i]) + { format!("{}", self.mu[d]) } else { "".into() } ) }; - let association = - if let Some(a) = self.association.assoc_comp.iter().position(|&a| a == i) { - format!( - "{}|{}|{}|{}", - self.association.kappa_ab[a], - self.association.epsilon_k_ab[a], - self.association.na[a], - self.association.nb[a] - ) - } else { - "|||".to_string() - }; + let record = gorup_dict[&self.identifiers[i]]; + let association = if let Some(a) = record.association_record { + format!( + "{}|{}|{}|{}|{}", + a.kappa_ab, a.epsilon_k_ab, a.na, a.nb, a.nc + ) + } else { + "||||".to_string() + }; write!( o, "\n|{}|{}|{}|{}|{}|{}|", component, self.identifiers[i], - self.m[i], - self.sigma[i], - self.epsilon_k[i], + record.m, + record.sigma, + record.epsilon_k, association ) .unwrap(); } + write!(o, "\n\n|component|group 1|group 2|bonds|\n|-|-|-|-|").unwrap(); let mut last_component = None; @@ -426,7 +434,7 @@ pub mod test { 2.7702, 334.29, None, - Some(AssociationRecord::new(0.009583, 2575.9, None, None)), + Some(AssociationRecord::new(0.009583, 2575.9, 1.0, 1.0, 0.0)), None, ), None, diff --git a/src/gc_pcsaft/python/mod.rs b/src/gc_pcsaft/python/mod.rs index 6ad10432a..171254a1e 100644 --- a/src/gc_pcsaft/python/mod.rs +++ b/src/gc_pcsaft/python/mod.rs @@ -67,7 +67,7 @@ impl PyGcPcSaftRecord { #[getter] fn get_association_record(&self) -> Option { - self.0.association_record.clone().map(PyAssociationRecord) + self.0.association_record.map(PyAssociationRecord) } fn __repr__(&self) -> PyResult { diff --git a/src/pcsaft/dft/mod.rs b/src/pcsaft/dft/mod.rs index fab80baad..deefc121d 100644 --- a/src/pcsaft/dft/mod.rs +++ b/src/pcsaft/dft/mod.rs @@ -76,7 +76,7 @@ impl PcSaftFunctional { contributions.push(Box::new(att)); // Association - if !parameters.association.assoc_comp.is_empty() { + if !parameters.association.is_empty() { let assoc = Association::new( ¶meters, ¶meters.association, diff --git a/src/pcsaft/dft/pure_saft_functional.rs b/src/pcsaft/dft/pure_saft_functional.rs index 66279ef3d..801a079a8 100644 --- a/src/pcsaft/dft/pure_saft_functional.rs +++ b/src/pcsaft/dft/pure_saft_functional.rs @@ -18,16 +18,18 @@ const PI36M1: f64 = 1.0 / (36.0 * PI); const N3_CUTOFF: f64 = 1e-5; const N0_CUTOFF: f64 = 1e-9; -#[derive(Clone)] pub struct PureFMTAssocFunctional { parameters: Arc, + association: Association, version: FMTVersion, } impl PureFMTAssocFunctional { pub fn new(parameters: Arc, version: FMTVersion) -> Self { + let association = Association::new(¶meters, ¶meters.association, 50, 1e-10); Self { parameters, + association, version, } } @@ -113,33 +115,22 @@ impl + ScalarOperand> FunctionalContributionDual for PureFMTA // association let a = &p.association; - if a.assoc_comp.len() == 1 { + if !a.is_empty() { let mut xi = -(&n2v * &n2v).sum_axis(Axis(0)) / (&n2 * &n2) + 1.0; xi.iter_mut().zip(&n2).for_each(|(xi, &n2)| { if n2.re() < N0_CUTOFF * 4.0 * PI * p.m[0] * r.re().powi(2) { *xi = N::one(); } }); - - let k = &n2 * &n3m1rec * r; - let deltarho = (((&k / 18.0 + 0.5) * &k * &xi + 1.0) * n3m1rec) - * ((temperature.recip() * a.epsilon_k_aibj[(0, 0)]).exp_m1() - * a.sigma3_kappa_aibj[(0, 0)]) - * (&n0 / p.m[0] * &xi); - - let f = |x: N| x.ln() - x * 0.5 + 0.5; - phi = phi - + if a.nb[0] > 0.0 { - let xa = deltarho.mapv(|d| { - Association::::assoc_site_frac_ab(d, a.na[0], a.nb[0]) - }); - let xb = (xa.clone() - 1.0) * a.na[0] / a.nb[0] + 1.0; - (n0 / p.m[0] * xi) * (xa.mapv(f) * a.na[0] + xb.mapv(f) * a.nb[0]) - } else { - let xa = deltarho - .mapv(|d| Association::::assoc_site_frac_a(d, a.na[0])); - n0 / p.m[0] * xi * (xa.mapv(f) * a.na[0]) - }; + let rho0 = (&n0 / p.m[0] * &xi).insert_axis(Axis(0)); + + phi += &(self.association.calculate_helmholtz_energy_density( + temperature, + &rho0, + &n2, + &n3m1rec, + &xi, + ))?; } Ok(phi) diff --git a/src/pcsaft/eos/mod.rs b/src/pcsaft/eos/mod.rs index 48f0b05fa..c11a7714a 100644 --- a/src/pcsaft/eos/mod.rs +++ b/src/pcsaft/eos/mod.rs @@ -86,7 +86,7 @@ impl PcSaft { variant: options.dq_variant, })); }; - if !parameters.association.assoc_comp.is_empty() { + if !parameters.association.is_empty() { contributions.push(Box::new(Association::new( ¶meters, ¶meters.association, diff --git a/src/pcsaft/eos/qspr.rs b/src/pcsaft/eos/qspr.rs index baa2e6436..0f2ce9b4b 100644 --- a/src/pcsaft/eos/qspr.rs +++ b/src/pcsaft/eos/qspr.rs @@ -70,12 +70,13 @@ pub struct QSPR { impl> IdealGasContributionDual for QSPR { fn de_broglie_wavelength(&self, temperature: D, components: usize) -> Array1 { - let (c_300, c_400) = match self.parameters.association.assoc_comp.len() { - 0 => match self.parameters.ndipole + self.parameters.nquadpole { + let (c_300, c_400) = if self.parameters.association.is_empty() { + match self.parameters.ndipole + self.parameters.nquadpole { 0 => (NA_NP_300, NA_NP_400), _ => (NA_P_300, NA_P_400), - }, - _ => (AP_300, AP_400), + } + } else { + (AP_300, AP_400) }; Array1::from_shape_fn(components, |i| { diff --git a/src/pcsaft/parameters.rs b/src/pcsaft/parameters.rs index 7053aa6a6..5cfef1564 100644 --- a/src/pcsaft/parameters.rs +++ b/src/pcsaft/parameters.rs @@ -7,7 +7,6 @@ use feos_core::parameter::{ }; use ndarray::{Array, Array1, Array2}; use num_dual::DualNum; -use num_traits::Zero; use quantity::si::{JOULE, KB, KELVIN}; use serde::{Deserialize, Serialize}; use std::collections::HashMap; @@ -70,14 +69,23 @@ impl FromSegments for PcSaftRecord { [ record.kappa_ab * n, record.epsilon_k_ab * n, - record.na.unwrap_or(1.0) * n, - record.nb.unwrap_or(1.0) * n, + record.na * n, + record.nb * n, + record.nc * n, ] }) }) - .reduce(|a, b| [a[0] + b[0], a[1] + b[1], a[2] + b[2], a[3] + b[3]]) - .map(|[kappa_ab, epsilon_k_ab, na, nb]| { - AssociationRecord::new(kappa_ab, epsilon_k_ab, Some(na), Some(nb)) + .reduce(|a, b| { + [ + a[0] + b[0], + a[1] + b[1], + a[2] + b[2], + a[3] + b[3], + a[4] + b[4], + ] + }) + .map(|[kappa_ab, epsilon_k_ab, na, nb, nc]| { + AssociationRecord::new(kappa_ab, epsilon_k_ab, na, nb, nc) }); // entropy scaling @@ -151,13 +159,23 @@ impl FromSegments for PcSaftRecord { impl FromSegments for PcSaftRecord { fn from_segments(segments: &[(Self, usize)]) -> Result { // We do not allow more than a single segment for q, mu, kappa_ab, epsilon_k_ab + let polar_segments: usize = segments + .iter() + .filter_map(|(s, n)| { + if s.q.is_some() || s.mu.is_some() || s.association_record.is_some() { + Some(n) + } else { + None + } + }) + .sum(); let quadpole_segments: usize = segments.iter().filter_map(|(s, n)| s.q.map(|_| n)).sum(); let dipole_segments: usize = segments.iter().filter_map(|(s, n)| s.mu.map(|_| n)).sum(); let assoc_segments: usize = segments .iter() .filter_map(|(s, n)| s.association_record.map(|_| n)) .sum(); - if quadpole_segments + dipole_segments + assoc_segments > 1 { + if polar_segments > 1 { return Err(ParameterError::IncompatibleParameters(format!( "Too many polar/associating segments (dipolar: {dipole_segments}, quadrupolar {quadpole_segments}, associating: {assoc_segments})." ))); @@ -209,14 +227,19 @@ impl PcSaftRecord { epsilon_k_ab: Option, na: Option, nb: Option, + nc: Option, viscosity: Option<[f64; 4]>, diffusion: Option<[f64; 5]>, thermal_conductivity: Option<[f64; 4]>, ) -> PcSaftRecord { let association_record = match (kappa_ab, epsilon_k_ab) { - (Some(kappa_ab), Some(epsilon_k_ab)) => { - Some(AssociationRecord::new(kappa_ab, epsilon_k_ab, na, nb)) - } + (Some(kappa_ab), Some(epsilon_k_ab)) => Some(AssociationRecord::new( + kappa_ab, + epsilon_k_ab, + na.unwrap_or(0.0), + nb.unwrap_or(0.0), + nc.unwrap_or(0.0), + )), (None, None) => None, _ => { panic!("To model association, both kappa_ab and epsilon_k_ab need to be specified.") @@ -328,7 +351,7 @@ impl Parameter for PcSaftParameters { epsilon_k[i] = r.epsilon_k; mu[i] = r.mu.unwrap_or(0.0); q[i] = r.q.unwrap_or(0.0); - association_records.push(r.association_record); + association_records.push(r.association_record.into_iter().collect()); viscosity.push(r.viscosity); diffusion.push(r.diffusion); thermal_conductivity.push(r.thermal_conductivity); @@ -460,7 +483,7 @@ impl PcSaftParameters { let o = &mut output; write!( o, - "|component|molarweight|$m$|$\\sigma$|$\\varepsilon$|$\\mu$|$Q$|$\\kappa_{{AB}}$|$\\varepsilon_{{AB}}$|$N_A$|$N_B$|\n|-|-|-|-|-|-|-|-|-|-|-|" + "|component|molarweight|$m$|$\\sigma$|$\\varepsilon$|$\\mu$|$Q$|$\\kappa_{{AB}}$|$\\varepsilon_{{AB}}$|$N_A$|$N_B$|$N_C$|\n|-|-|-|-|-|-|-|-|-|-|-|-|" ) .unwrap(); for (i, record) in self.pure_records.iter().enumerate() { @@ -469,10 +492,10 @@ impl PcSaftParameters { let association = record .model_record .association_record - .unwrap_or_else(|| AssociationRecord::new(0.0, 0.0, None, None)); + .unwrap_or_else(|| AssociationRecord::new(0.0, 0.0, 0.0, 0.0, 0.0)); write!( o, - "\n|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|", + "\n|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|", component, record.molarweight, record.model_record.m, @@ -482,8 +505,9 @@ impl PcSaftParameters { record.model_record.q.unwrap_or(0.0), association.kappa_ab, association.epsilon_k_ab, - association.na.unwrap_or(1.0), - association.nb.unwrap_or(1.0) + association.na, + association.nb, + association.nc ) .unwrap(); } @@ -492,33 +516,6 @@ impl PcSaftParameters { } } -impl std::fmt::Display for PcSaftParameters { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - write!(f, "PcSaftParameters(")?; - write!(f, "\n\tmolarweight={}", self.molarweight)?; - write!(f, "\n\tm={}", self.m)?; - write!(f, "\n\tsigma={}", self.sigma)?; - write!(f, "\n\tepsilon_k={}", self.epsilon_k)?; - if !self.dipole_comp.is_empty() { - write!(f, "\n\tmu={}", self.mu)?; - } - if !self.quadpole_comp.is_empty() { - write!(f, "\n\tq={}", self.q)?; - } - if !self.association.assoc_comp.is_empty() { - write!(f, "\n\tassociating={}", self.association.assoc_comp)?; - write!(f, "\n\tkappa_ab={}", self.association.kappa_ab)?; - write!(f, "\n\tepsilon_k_ab={}", self.association.epsilon_k_ab)?; - write!(f, "\n\tna={}", self.association.na)?; - write!(f, "\n\tnb={}", self.association.nb)?; - } - if !self.k_ij.iter().all(|k| k.is_zero()) { - write!(f, "\n\tk_ij=\n{}", self.k_ij)?; - } - write!(f, "\n)") - } -} - #[cfg(test)] pub mod utils { use super::*; @@ -639,7 +636,9 @@ pub mod utils { "sigma": 3.000683, "epsilon_k": 366.5121, "kappa_ab": 0.034867983, - "epsilon_k_ab": 2500.6706 + "epsilon_k_ab": 2500.6706, + "na": 1.0, + "nb": 1.0 }, "molarweight": 18.0152 }"#; diff --git a/src/pcsaft/python.rs b/src/pcsaft/python.rs index 658722422..077e4d7e6 100644 --- a/src/pcsaft/python.rs +++ b/src/pcsaft/python.rs @@ -36,6 +36,7 @@ impl PyPcSaftRecord { epsilon_k_ab: Option, na: Option, nb: Option, + nc: Option, viscosity: Option<[f64; 4]>, diffusion: Option<[f64; 5]>, thermal_conductivity: Option<[f64; 4]>, @@ -50,6 +51,7 @@ impl PyPcSaftRecord { epsilon_k_ab, na, nb, + nc, viscosity, diffusion, thermal_conductivity, @@ -93,12 +95,17 @@ impl PyPcSaftRecord { #[getter] fn get_na(&self) -> Option { - self.0.association_record.and_then(|a| a.na) + self.0.association_record.map(|a| a.na) } #[getter] fn get_nb(&self) -> Option { - self.0.association_record.and_then(|a| a.nb) + self.0.association_record.map(|a| a.nb) + } + + #[getter] + fn get_nc(&self) -> Option { + self.0.association_record.map(|a| a.nc) } #[getter] @@ -172,10 +179,6 @@ impl PyPcSaftParameters { fn _repr_markdown_(&self) -> String { self.0.to_markdown() } - - fn __repr__(&self) -> PyResult { - Ok(self.0.to_string()) - } } #[pymodule] diff --git a/tests/gc_pcsaft/dft.rs b/tests/gc_pcsaft/dft.rs index b777405b5..cc5a76827 100644 --- a/tests/gc_pcsaft/dft.rs +++ b/tests/gc_pcsaft/dft.rs @@ -120,13 +120,13 @@ fn test_bulk_association() -> Result<(), Box> { segment_records.clone(), None, )?); - let eos = Arc::new(GcPcSaft::new(eos_parameters.clone())); + let eos = Arc::new(GcPcSaft::new(eos_parameters)); let func_parameters = Arc::new(GcPcSaftFunctionalParameters::from_segments( vec![ethylene_glycol], segment_records, None, )?); - let func = Arc::new(GcPcSaftFunctional::new(func_parameters.clone())); + let func = Arc::new(GcPcSaftFunctional::new(func_parameters)); let t = 200.0 * KELVIN; let v = 0.002 * METER.powi(3); @@ -135,35 +135,9 @@ fn test_bulk_association() -> Result<(), Box> { let state_func = State::new_nvt(&func, t, v, &n)?; let p_eos = state_eos.pressure_contributions(); let p_func = state_func.pressure_contributions(); - println!( - "Equation of state: - \tcomps: {} - \tkappa_ab: {} - \tepsilon_k_ab: {} - \tna: {} - \tnb: {}", - eos_parameters.association.assoc_comp, - eos_parameters.association.kappa_ab, - eos_parameters.association.epsilon_k_ab, - eos_parameters.association.na, - eos_parameters.association.nb, - ); for (s, x) in &p_eos { println!("{s:18}: {x:21.16}"); } - println!( - "\nHelmholtz energy functional: - \tcomps: {} - \tkappa_ab: {} - \tepsilon_k_ab: {} - \tna: {} - \tnb: {}", - func_parameters.association.assoc_comp, - func_parameters.association.kappa_ab, - func_parameters.association.epsilon_k_ab, - func_parameters.association.na, - func_parameters.association.nb, - ); for (s, x) in &p_func { println!("{s:26}: {x:21.16}"); } diff --git a/tests/pcsaft/test_parameters.json b/tests/pcsaft/test_parameters.json index 988fcb3b8..c49a5de31 100644 --- a/tests/pcsaft/test_parameters.json +++ b/tests/pcsaft/test_parameters.json @@ -93,7 +93,9 @@ "sigma": 3.000683, "epsilon_k": 366.5121, "kappa_ab": 0.034867983, - "epsilon_k_ab": 2500.6706 + "epsilon_k_ab": 2500.6706, + "na": 1.0, + "nb": 1.0 }, "molarweight": 18.0152 },