More derivatives

bee14d51 · Ian Bell · 54e94f4b · bee14d51
Commit bee14d51 authored 2 years ago by Ian Bell
--- a/doc/source/derivs/derivs.ipynb
+++ b/doc/source/derivs/derivs.ipynb
@@ -24,7 +24,42 @@
    "\n",
    "$$ p^r = \\rho R T \\Lambda^r_{01}$$\n",
    "\n",
-    "Similar definitions apply for all the other residual thermodynamic properties"
+    "Similar definitions apply for all the other thermodynamic properties, with the tot superscript indicating it is the sum of the residual and ideal-gas (not included in teqp) contributions:\n",
+    "\n",
+    "$$\n",
+    "\\frac{p}{\\rho R T}  = 1 + \\Lambda^r_{01}\n",
+    "$$\n",
+    "Internal energy ($u= a+Ts$):\n",
+    "$$\n",
+    "\t\\frac{u}{RT} = \\Lambda^{\\rm tot}_{10}\n",
+    "$$\n",
+    "Enthalpy ($h= u+p/\\rho$):\n",
+    "$$\n",
+    "\\frac{h}{RT} = 1+\\Lambda^r_{01} + \\Lambda^{\\rm tot}_{10}\n",
+    "$$\n",
+    "Entropy ($s\\equiv -(\\partial a/\\partial T)_v$):\n",
+    "$$\n",
+    "\\frac{s}{R} = \\Lambda^{\\rm tot}_{10}-\\Lambda^{\\rm tot}_{00}\n",
+    "$$\n",
+    "Gibbs energy ($g= h-Ts$):\n",
+    "$$\n",
+    "\t\\frac{g}{RT} = 1+\\Lambda^r_{01}+\\Lambda^{\\rm tot}_{00}\n",
+    "$$\n",
+    "Derivatives of pressure:\n",
+    "$$\n",
+    "\t\\left(\\frac{\\partial p}{\\partial \\rho}\\right)_{T} = RT\\left(1+2\\Lambda^r_{01}+\\Lambda^r_{02}\\right)\n",
+    "$$\n",
+    "$$\n",
+    "\t\\left(\\frac{\\partial p}{\\partial T}\\right)_{\\rho} = R\\rho\\left(1+\\Lambda^r_{01}-\\Lambda^r_{11}\\right)\n",
+    "$$\n",
+    "Isochoric specific heat ($c_v\\equiv (\\partial u/\\partial T)_v$):\n",
+    "$$\n",
+    "\\frac{c_v}{R} = -\\Lambda^{\\rm tot}_{20}\n",
+    "$$\n",
+    "Isobaric specific heat ($c_p\\equiv (\\partial h/\\partial T)_p$; see Eq. 3.56 from Span \\cite{Span-BOOK-2000} for the derivation):\n",
+    "$$\n",
+    "\\frac{c_p}{R} = -\\Lambda^{\\rm tot}_{20}+\\frac{(1+\\Lambda^r_{01}-\\Lambda^r_{11})^2}{1+2\\Lambda^r_{01}+\\Lambda^r_{02}}\n",
+    "$$"
   ]
  },
  {

 %% Cell type:markdown id:1afc4517 tags:

 # Thermodynamic Derivatives

 ## Helmholtz energy derivatives
 Thermodynamic derivatives are at the very heart of teqp.  All models are defined in the form $\alpha^r(T, \rho, z)$, where $\rho$ is the molar density, and z are mole fractions.  There are exceptions for models for which the independent variables are in simulation units (Lennard-Jones and its ilk).

 Thereofore, to obtain the residual pressure, it is obtained as a derivative:

 $$ p^r = \rho R T \left( \rho \left(\frac{\partial \alpha^r}{\partial \rho}\right)_{T}\right)$$

 and other residual thermodynamic properties are defined likewise.

 We can define the concise derivative

 $$ \Lambda^r_{xy} = (1/T)^x(\rho)^y\left(\frac{\partial^{x+y}(\alpha^r)}{\partial (1/T)^x\partial \rho^y}\right)$$

 so we can re-write the derivative above as

 $$ p^r = \rho R T \Lambda^r_{01}$$

-Similar definitions apply for all the other residual thermodynamic properties
+Similar definitions apply for all the other thermodynamic properties, with the tot superscript indicating it is the sum of the residual and ideal-gas (not included in teqp) contributions:
+
+$$
+\frac{p}{\rho R T}  = 1 + \Lambda^r_{01}
+$$
+Internal energy ($u= a+Ts$):
+$$
+	\frac{u}{RT} = \Lambda^{\rm tot}_{10}
+$$
+Enthalpy ($h= u+p/\rho$):
+$$
+\frac{h}{RT} = 1+\Lambda^r_{01} + \Lambda^{\rm tot}_{10}
+$$
+Entropy ($s\equiv -(\partial a/\partial T)_v$):
+$$
+\frac{s}{R} = \Lambda^{\rm tot}_{10}-\Lambda^{\rm tot}_{00}
+$$
+Gibbs energy ($g= h-Ts$):
+$$
+	\frac{g}{RT} = 1+\Lambda^r_{01}+\Lambda^{\rm tot}_{00}
+$$
+Derivatives of pressure:
+$$
+	\left(\frac{\partial p}{\partial \rho}\right)_{T} = RT\left(1+2\Lambda^r_{01}+\Lambda^r_{02}\right)
+$$
+$$
+	\left(\frac{\partial p}{\partial T}\right)_{\rho} = R\rho\left(1+\Lambda^r_{01}-\Lambda^r_{11}\right)
+$$
+Isochoric specific heat ($c_v\equiv (\partial u/\partial T)_v$):
+$$
+\frac{c_v}{R} = -\Lambda^{\rm tot}_{20}
+$$
+Isobaric specific heat ($c_p\equiv (\partial h/\partial T)_p$; see Eq. 3.56 from Span \cite{Span-BOOK-2000} for the derivation):
+$$
+\frac{c_p}{R} = -\Lambda^{\rm tot}_{20}+\frac{(1+\Lambda^r_{01}-\Lambda^r_{11})^2}{1+2\Lambda^r_{01}+\Lambda^r_{02}}
+$$

 %% Cell type:code id:ca35cc29 tags:

 ``` python
 import teqp
 import numpy as np
 ```

 %% Cell type:code id:8d254c40 tags:

 ``` python
 Tc_K = [300]
 pc_Pa = [4e6]
 acentric = [0.01]
 model = teqp.canonical_PR(Tc_K, pc_Pa, acentric)
 ```

 %% Cell type:code id:49ee3453 tags:

 ``` python
 z = np.array([1.0])
 model.get_Ar01(300,300,z)
 ```

 %% Output

    -0.06836660379313926

 %% Cell type:markdown id:498a2350 tags:

 And there are additional methods to obtain all the derivatives up to a given order:

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

 ``` python
 model.get_Ar06n(300,300,z) # derivatives 00, 01, 02, ... 06
 ```

 %% Output

    [-0.06966138343515413,
     -0.06836660379313926,
     0.0025357822532378147,
     -0.00015701162203571184,
     1.6818628788290574e-05,
     -2.2305940927885907e-06,
     3.8259258513417917e-07]

 %% Cell type:markdown id:cd644a3e tags:

 But more derivatives are slower than fewer:

 %% Cell type:code id:4e0609ae tags:

 ``` python
 %timeit model.get_Ar01(300,300,z)
 %timeit model.get_Ar04n(300,300,z)
 ```

 %% Output

    1.9 µs ± 38.5 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
    2.07 µs ± 24.8 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

 %% Cell type:markdown id:d05e9070 tags:

 Note: calling overhead is usually on the order of 1 microsecond

 %% Cell type:markdown id:502e9830 tags:

 ## Virial coefficients

 Virial coefficients represent the thermodynamics of the interaction of two-, three-, ... bodies interacting with each other.  They can be obtained rigorously if the potential energy surface of interaction is fully known.  In general, such a surface can only be constructed for small rigid molecules.  Many simple thermodynamic models do a poor job of predicting the thermodynamics captured by the virial coefficients.

 The i-th virial coefficient is defined by
 $$
 	B_i = \frac{(\alpha^r)^{(i-1)}}{(i-2)!}
 $$
 with the concise derivative term
 $$
 	(\alpha^r)^{(i)} = \lim_{\rho \to 0} \left(\frac{\partial^{i}\alpha^r}{\partial \rho^{i}}\right)_{T,\vec{x}}
 $$

 teqp supports the virial coefficient directly, there is the ``get_B2vir`` method for the second virial coefficient:

 %% Cell type:code id:720e120c tags:

 ``` python
 model.get_B2vir(300, z)
 ```

 %% Output

    -0.00023661263734465424

 %% Cell type:markdown id:2fba4369 tags:

 And the ``get_Bnvir`` method that allows for the calculation of higher virial coefficients:

 %% Cell type:code id:92d01992 tags:

 ``` python
 model.get_Bnvir(7, 300, z)
 ```

 %% Output

    {2: -0.0002366126373446542,
     3: 3.001768410777936e-08,
     4: -3.2409760373816355e-12,
     5: 3.9617816466337214e-16,
     6: -4.552923983836698e-20,
     7: 5.3759278511184914e-24}

 %% Cell type:markdown id:d50998b1 tags:

 The ``get_Bnvir`` method was implemented because when doing automatic differentiation, all the intermediate derivatives are also retained.

 There is also a method to calculate temperature derivatives of a given virial coefficient

 %% Cell type:code id:fa4ff386 tags:

 ``` python
 model.get_dmBnvirdTm(2, 3, 300, z) # third temperature derivative of the second virial coefficient
 ```

 %% Output

    1.0095625628421257e-10