Add example of calling superancillary for PC-SAFT to docs

d13a4641 · Ian Bell · c4494a3a · d13a4641
Commit d13a4641 authored 2 years ago by Ian Bell
--- a/doc/source/models/PCSAFT.ipynb
+++ b/doc/source/models/PCSAFT.ipynb
@@ -104,7 +104,51 @@
   "source": [
    "## Superancillary\n",
    "\n",
-    "The superancillary equation for PC-SAFT has been developed, and is much more involved than that of the cubic EOS. As a consequence, the superancillary equation has been provided as a separate package rather than integrating it into to teqp to minimize the binary size of teqp. It can be installed from PYPI with: ``pip install PCSAFTsuperanc``"
+    "The superancillary equation for PC-SAFT has been developed, and is much more involved than that of the cubic EOS. As a consequence, the superancillary equation has been provided as a separate package rather than integrating it into to teqp to minimize the binary size of teqp. It can be installed from PYPI with: ``pip install PCSAFTsuperanc``\n",
+    "\n",
+    "The scaling in the superancillaries uses reduced variables:\n",
+    "\n",
+    "$$ \\tilde T = T/(\\epsilon/k_{\\rm B}) $$\n",
+    "$$ \\tilde\\rho = \\rho_{\\rm N}\\sigma^3 $$\n",
+    "\n",
+    "where $\\rho_{\\rm N}$ is the number density, and the other parameters are from the PC-SAFT model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f6e3b8d2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import PCSAFTsuperanc\n",
+    "\n",
+    "sigma_m = 3e-10 # [meter]\n",
+    "e_over_k = 150.0 # [K]\n",
+    "m = 5\n",
+    "\n",
+    "# The saturation temperature\n",
+    "T = 300\n",
+    "\n",
+    "[Ttilde_crit, Ttilde_min] = PCSAFTsuperanc.get_Ttilde_crit_min(m=m)\n",
+    "print('Ttilde crit:', Ttilde_crit)\n",
+    "\n",
+    "# Get the scaled densities for liquid and vapor phases\n",
+    "[tilderhoL, tilderhoV] = PCSAFTsuperanc.PCSAFTsuperanc_rhoLV(Ttilde=T/e_over_k, m=m)\n",
+    "# Convert back to molar densities\n",
+    "N_A = PCSAFTsuperanc.N_A # The value of Avogadro's constant used in superancillaries\n",
+    "rhoL, rhoV = [tilderho/(N_A*sigma_m**3) for tilderho in [tilderhoL, tilderhoV]]\n",
+    "\n",
+    "# As a sanity check, confirm that we got the same pressure in both phases\n",
+    "c = teqp.SAFTCoeffs()\n",
+    "c.sigma_Angstrom = sigma_m*1e10\n",
+    "c.epsilon_over_k = e_over_k \n",
+    "c.m = m\n",
+    "model = teqp.PCSAFTEOS([c])\n",
+    "z = np.array([1.0])\n",
+    "pL = rhoL*model.get_R(z)*T*(1+model.get_Ar01(T, rhoL, z))\n",
+    "pV = rhoV*model.get_R(z)*T*(1+model.get_Ar01(T, rhoV, z))\n",
+    "print('Pressures are:', pL, pV, 'Pa')"
   ]
  },
  {

 %% Cell type:markdown id:f406bbb5 tags:

 # PC-SAFT

 The PC-SAFT implementation in teqp is based on the implementation of Gross and Sadowski (https://doi.org/10.1021/ie0003887), with the typo from their paper fixed.  It does NOT include the association contribution, only the dispersive contributions.

 The model in teqp requires the user to specify the values of ``sigma``, ``epsilon/kB``, and ``m`` for each substance.  A very few substances are hardcoded in teqp, for testing purposes.

 %% Cell type:raw id:d9efd027 tags:

 The Python class is here: :py:class:`PCSAFTEOS <teqp.teqp.PCSAFTEOS>`

 %% Cell type:code id:984925ce tags:

 ``` python
 import teqp
 import numpy as np
 teqp.__version__
 ```

 %% Cell type:code id:7bbd7129 tags:

 ``` python
 TeXkey = 'Gross-IECR-2001'
 ms = [1.0, 1.6069, 2.0020]
 eoverk = [150.03, 191.42, 208.11]
 sigmas = [3.7039, 3.5206, 3.6184]

 coeffs = []
 for i in range(len(ms)):
    c = teqp.SAFTCoeffs()
    c.m = ms[i]
    c.epsilon_over_k = eoverk[i]
    c.sigma_Angstrom = sigmas[i]
    coeffs.append(c)

 model = teqp.PCSAFTEOS(coeffs)
 ```

 %% Cell type:markdown id:578630c8 tags:

 The model parameters can be queried:

 %% Cell type:code id:d4e47e54 tags:

 ``` python
 model.get_m(), model.get_epsilon_over_k_K(), model.get_sigma_Angstrom()
 ```

 %% Cell type:markdown id:5cbb382d tags:

 ## Adjusting k_ij

 Fine-tuned values of $k_{ij}$ can be provided when instantiating the model.  A complete matrix of all the $k_{ij}$ values must be provided. This allows for asymmetric mixing models in which $k_{ij}\neq k_{ji}$.

 %% Cell type:code id:a32c41b5 tags:

 ``` python
 k_01 = 0.01; k_10 = k_01
 kmat = [[0,k_01,0],[k_10,0,0],[0,0,0]]
 teqp.PCSAFTEOS(coeffs, kmat)
 ```

 %% Cell type:markdown id:ca52e844 tags:

 ## Superancillary

 The superancillary equation for PC-SAFT has been developed, and is much more involved than that of the cubic EOS. As a consequence, the superancillary equation has been provided as a separate package rather than integrating it into to teqp to minimize the binary size of teqp. It can be installed from PYPI with: ``pip install PCSAFTsuperanc``

+The scaling in the superancillaries uses reduced variables:
+
+$$ \tilde T = T/(\epsilon/k_{\rm B}) $$
+$$ \tilde\rho = \rho_{\rm N}\sigma^3 $$
+
+where $\rho_{\rm N}$ is the number density, and the other parameters are from the PC-SAFT model
+
+%% Cell type:code id:f6e3b8d2 tags:
+
+``` python
+import PCSAFTsuperanc
+
+sigma_m = 3e-10 # [meter]
+e_over_k = 150.0 # [K]
+m = 5
+
+# The saturation temperature
+T = 300
+
+[Ttilde_crit, Ttilde_min] = PCSAFTsuperanc.get_Ttilde_crit_min(m=m)
+print('Ttilde crit:', Ttilde_crit)
+
+# Get the scaled densities for liquid and vapor phases
+[tilderhoL, tilderhoV] = PCSAFTsuperanc.PCSAFTsuperanc_rhoLV(Ttilde=T/e_over_k, m=m)
+# Convert back to molar densities
+N_A = PCSAFTsuperanc.N_A # The value of Avogadro's constant used in superancillaries
+rhoL, rhoV = [tilderho/(N_A*sigma_m**3) for tilderho in [tilderhoL, tilderhoV]]
+
+# As a sanity check, confirm that we got the same pressure in both phases
+c = teqp.SAFTCoeffs()
+c.sigma_Angstrom = sigma_m*1e10
+c.epsilon_over_k = e_over_k
+c.m = m
+model = teqp.PCSAFTEOS([c])
+z = np.array([1.0])
+pL = rhoL*model.get_R(z)*T*(1+model.get_Ar01(T, rhoL, z))
+pV = rhoV*model.get_R(z)*T*(1+model.get_Ar01(T, rhoV, z))
+print('Pressures are:', pL, pV, 'Pa')
+```
+
 %% Cell type:markdown id:0bdf568f tags:

 ## Maximum density

 The maximum number density allowed by the EOS is defined based on the packing fraction. To get a molar density, divide by Avogadro's number. The function is conveniently exposed in Python:

 %% Cell type:code id:3c8491a9 tags:

 ``` python
 max_rhoN = teqp.PCSAFTEOS(coeffs).max_rhoN(130.0, np.array([0.3, 0.3, 0.4]))
 display(max_rhoN)
 max_rhoN/6.022e23 # the maximum molar density in mol/m^3
 ```