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 RT\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

$${\mathrm{\Lambda }}_{xy}^r=(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 RT{\mathrm{\Lambda }}_{01}^r$$

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 RT}=1+{\mathrm{\Lambda }}_{01}^r$$

Internal energy ($u=a+Ts$):

$$\frac{u}{RT}={\mathrm{\Lambda }}_{10}^{\mathrm{t}\mathrm{o}\mathrm{t}}$$

Enthalpy ($h=u+p/\rho$):

$$\frac{h}{RT}=1+{\mathrm{\Lambda }}_{01}^r+{\mathrm{\Lambda }}_{10}^{\mathrm{t}\mathrm{o}\mathrm{t}}$$

Entropy ($s\equiv -(\partial a/\partial T{)}_v$):

$$\frac{s}{R}={\mathrm{\Lambda }}_{10}^{\mathrm{t}\mathrm{o}\mathrm{t}}-{\mathrm{\Lambda }}_{00}^{\mathrm{t}\mathrm{o}\mathrm{t}}$$

Gibbs energy ($g=h-Ts$):

$$\frac{g}{RT}=1+{\mathrm{\Lambda }}_{01}^r+{\mathrm{\Lambda }}_{00}^{\mathrm{t}\mathrm{o}\mathrm{t}}$$

Derivatives of pressure:

$${\left(\frac{\partial p}{\partial \rho }\right)}_T=RT(1+2{\mathrm{\Lambda }}_{01}^r+{\mathrm{\Lambda }}_{02}^r)$$

$${\left(\frac{\partial p}{\partial T}\right)}_{\rho }=R\rho (1+{\mathrm{\Lambda }}_{01}^r-{\mathrm{\Lambda }}_{11}^r)$$

Isochoric specific heat ($c_v\equiv (\partial u/\partial T{)}_v$):

$$\frac{c_v}{R}=-{\mathrm{\Lambda }}_{20}^{\mathrm{t}\mathrm{o}\mathrm{t}}$$

Isobaric specific heat ($c_p\equiv (\partial h/\partial T{)}_p$; see Eq. 3.56 from Span for the derivation):

$$\frac{c_p}{R}=-{\mathrm{\Lambda }}_{20}^{\mathrm{t}\mathrm{o}\mathrm{t}}+\frac{(1+{\mathrm{\Lambda }}_{01}^r-{\mathrm{\Lambda }}_{11}^r{)}^2}{1+2{\mathrm{\Lambda }}_{01}^r+{\mathrm{\Lambda }}_{02}^r}$$

In teqp, these derivatives are obtained from methods like

• get_Arxy()

• get_Ar06n()

where the A in this context indicates the variable $\mathrm{\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 $\mathrm{\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)}=\underset{\rho \to 0}{lim}{\left(\frac{{\partial }^i{\alpha }^r}{\partial {\rho }^i}\right)}_{T,\overrightarrow{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 $\mathrm{\Psi }$ as a function of temperature and molar densities $\overrightarrow{\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=-\mathrm{\Psi }+\sum _{i=1}^N{\rho }_i{\mu }_i$$

with the chemical potential ${\mu }_i$ obtained from

$${\mu }_i={\left(\frac{\partial \mathrm{\Psi }}{\partial {\rho }_i}\right)}_{T,{\rho }_{j\ne 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=$\overrightarrow{\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])