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 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 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}}\]

In teqp, these derivatives are obtained from methods like

• get_Arxy()

• get_Ar06n()

where the A in this context indicates the variable \(\Lambda\) above. This naming is perhaps not ideal, since \(A\) is sometimes the total Helmholtz energy, but it was a close visual mnemonic to the character \(\Lambda\).

import teqp
teqp.__version__

'0.22.0'

import numpy as np

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

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

-0.06836660379313926

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

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

array([-6.96613834e-02, -6.83666038e-02,  2.53578225e-03, -1.57011622e-04,
        1.68186288e-05, -2.23059409e-06,  3.82592585e-07])

But more derivatives are slower than fewer:

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

567 ns ± 0.747 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)
1.11 μs ± 4.12 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

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

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:

model.get_B2vir(300, z)

-0.00023661263734465424

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

model.get_Bnvir(7, 300, z)

{2: -0.0002366126373446542,
 3: 3.001768410777936e-08,
 4: -3.2409760373816364e-12,
 5: 3.961781646633723e-16,
 6: -4.5529239838367004e-20,
 7: 5.375927851118494e-24}

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

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

1.0095625628421257e-10

Isochoric Thermodynamics Derivatives

In the isochoric thermodynamics formalism, the EOS is expressed in the Helmholtz energy density \(\Psi\) as a function of temperature and molar densities \(\vec\rho\). This formalism is handy because it allows for a concise mathematical structure, well suited to implementation in teqp. For instance the pressure is obtained from (see https://doi.org/10.1002/aic.16074):

\[p=-\Psi + \sum_{i=1}^N\rho_i\mu_i\]

with the chemical potential \(\mu_i\) obtained from

\[\mu_i = \left(\frac{\partial \Psi}{\partial \rho_i}\right)_{T,\rho_{j\neq i}}\]

The molar densities \(\rho_i\) are related to the total density and the mole fractions:

\[\rho_i = x_i\rho\]

In teqp, the isochoric derivative functions like get_fugacity_coefficients, get_partial_molar_volumes take as arguments the temperature T and the vector of molar concentrations rhovec=\(\vec\rho\), which are obtained by multiplying the mole fractions by the total density.

Example:

model = teqp.build_multifluid_model(["CO2","Argon"], teqp.get_datapath())
T, rhovec = 300, np.array([0.3,0.4])*300 # K, mol/m^3
display(model.get_fugacity_coefficients(T, rhovec))
display(model.get_partial_molar_volumes(T, rhovec))

array([0.97884567, 0.99866748])

array([0.00470644, 0.00480351])