Composition derivatives

For other mixture calculations composition derivatives of the form

\[\Lambda^r_{xyz_i} = (1/T)^x(\rho)^y\left(\frac{\partial^{x+y+z_i}(\alpha^r)}{\partial (1/T)^x\partial \rho^y \partial \mathbf{Z}_i^{z_i}} \right)\]

are needed. This function is exposed in teqp (as of version 0.19) as the function get_ATrhoXi. In order to limit the binary size and compilation time, x has a max of 2 and y does as well. z_i can be up to 3, and must be at least 1, otherwise you can use the other derivative functions that do not require any composition derivatives.

The mixed composition derivative of the form

\[\Lambda^r_{xyz_i z_j } = (1/T)^x(\rho)^y\left(\frac{\partial^{x+y+z_i+z_j}(\alpha^r)}{\partial (1/T)^x\partial \rho^y \partial \mathbf{Z}_i^{z_i} \partial \mathbf{Z}_j^{z_j} }\right)\]

supports \(x\) and \(y\) either 0 or 1, with at most two composition derivatives. The triple composition derivative of the form

\[\Lambda^r_{xyz_i z_j z_k} = (1/T)^x(\rho)^y\left(\frac{\partial^{x+y+z_i+z_j+z_k}(\alpha^r)}{\partial (1/T)^x\partial \rho^y \partial \mathbf{Z}_i^{z_i} \partial \mathbf{Z}_j^{z_j} \partial \mathbf{Z}_k^{z_k} }\right)\]

supports \(x\) and \(y\) either 0 or 1, with up to first derivatives in each composition variable. If this is not enough derivatives, open a feature request here : https://github.com/usnistgov/teqp/issues

\(\tau\) and \(\delta\) derivatives

In the multi-fluid modeling approach used in NIST REFPROP and the GERG-2004 GERG-2008 models, the derivatives are in the form

\[\Lambda^r_{xyz_i} = \tau^x\delta^y\left(\frac{\partial^{x+y+z_i}(\alpha^r)}{\partial \tau^x\partial \delta^y \partial \mathbf{Z}_i^{z_i}} \right)\]

with \(\tau=T_{\rm red}(\mathbf{Z})/T\) and \(\delta=\rho/\rho_{\rm red}(\mathbf{Z})\). The higher derivatives are similarly equal to

\[\Lambda^r_{xyz_i z_j } = \tau^x\delta^y\left(\frac{\partial^{x+y+z_i+z_j}(\alpha^r)}{\partial \tau^x\partial \delta^y \partial \mathbf{Z}_i^{z_i} \partial \mathbf{Z}_j^{z_j} }\right)\]
\[\Lambda^r_{xyz_i z_j z_k} = \tau^x\delta^y\left(\frac{\partial^{x+y+z_i+z_j+z_k}(\alpha^r)}{\partial \tau^x\partial \delta^y \partial \mathbf{Z}_i^{z_i} \partial \mathbf{Z}_j^{z_j} \partial \mathbf{Z}_k^{z_k} }\right)\]

The same limitations on the numbers of derivatives are used for the derivatives with \((1/T)\) and \(\rho\) as independent variables.

The Python methods are documented here:

xN (in)dependent

Let’s suppose that some quantity \(\Upsilon\) depends on mole fractions. If all the mole fractions are considered to be independent, the total differential is obtained from \(\renewcommand{\deriv}[3]{\left(\frac{\partial #1}{\partial #2}\right)_{#3}}\)

\[{\rm d} \Upsilon = \sum_j\deriv{\Upsilon}{x_j}{x_{i\neq j}} {\rm d}x_j\]

If instead the last mole fraction is defined to be dependent on the others via

\[x_N = 1-\sum_{i=1}^{N-1}x_i\]

then the total differential is obtained from

\[{\rm d} \Upsilon = \sum_{j=1}^{N-1}\deriv{\Upsilon}{x_j}{x_{i\neq j}} {\rm d}x_j + \deriv{\Upsilon}{x_N}{x_{i\neq j}} {\rm d}{x_N}\]

where \(x_N\) is considered to be an independent variable in the derivative \(\deriv{\Upsilon}{x_N}{x_{i\neq j}}\). Thus derivatives with respect to one of the dependent mole fractions (\(x_k\) with \(k < N\)) would be equal to

\[\deriv{\Upsilon}{x_k}{x_{j\neq k}} =\deriv{\Upsilon}{x_k}{x_{i\neq k}} - \deriv{\Upsilon}{x_N}{x_{i\neq N}}\]

because

\[\deriv{x_N}{x_i}{} = -1\]

So if the library (e.g., CoolProp and TREND) allows for the fractions to be dependent (either option is allowed in CoolProp, TREND uses N-1 independent mole fractions), you can use the molar composition derivatives with the mole fractions treated as being independent to obtain derivatives with one of the mole fractions dependent on the other \(N-1\) fractions.

[1]:

import teqp, numpy as np
teqp.__version__

[1]:

'0.22.0'

[2]:

j = {
    'components': ["Methane", "Nitrogen", "Oxygen"],
    'root': teqp.get_datapath(),
    'BIP': '',
    'departure': ''
}
model = teqp.make_model({'kind':'multifluid', 'model': j})


T = 300 # K
rhomolar = 3000 # mol/m^3
z = np.array([0.3, 0.5, 0.2]) # mole fractions

Tr = model.get_Tr(z)
rhor = model.get_rhor(z)
tau = Tr/T
delta = rhomolar/rhor

Ntau = 0
Ndelta = 0
Nxi = 1
print(model.get_AtaudeltaXi(tau, Ntau, delta, Ndelta, z, 0, Nxi))

Ntau = 1
Ndelta = 0
Nxi = 1
print(model.get_AtaudeltaXi(tau, Ntau, delta, Ndelta, z, 0, Nxi))

Ntau = 0
Ndelta = 1
Nxi = 1
print(model.get_AtaudeltaXi(tau, Ntau, delta, Ndelta, z, 0, Nxi))

Ntau = 2
Ndelta = 0
Nxi = 1
print(model.get_AtaudeltaXi(tau, Ntau, delta, Ndelta, z, 0, Nxi))

Ntau = 1
Ndelta = 1
Nxi = 1
print(model.get_AtaudeltaXi(tau, Ntau, delta, Ndelta, z, 0, Nxi))

Ntau = 0
Ndelta = 2
Nxi = 1
print(model.get_AtaudeltaXi(tau, Ntau, delta, Ndelta, z, 0, Nxi))

Ntau = 1
Ndelta = 0
Nxi = 1
Nxj = 1
print(model.get_AtaudeltaXiXj(tau, Ntau, delta, Ndelta, z, 0, Nxi, 1, Nxj))

Ntau = 0
Ndelta = 1
Nxi = 1
Nxj = 1
print(model.get_AtaudeltaXiXj(tau, Ntau, delta, Ndelta, z, 0, Nxi, 1, Nxj))

Ntau = 0
Ndelta = 0
Nxi = 1
Nxj = 1
Nxk = 1
print(model.get_AtaudeltaXiXjXk(tau, Ntau, delta, Ndelta, z, 0, Nxi, 1, Nxj, 2, Nxk))

-0.043587384253511226
-0.2118857998812584
-0.03650566667904927
-0.07488856488580686
-0.2069389009652925
0.014468933385218782
-0.005978809921279949
-0.00279185550001082
0.0

With CoolProp, version 6.6.0, the following script in C++:

#include "AbstractState.h"
#include "Backends/Helmholtz/MixtureDerivatives.h"

int main(){
    std::shared_ptr<CoolProp::AbstractState> AS(
        CoolProp::AbstractState::factory("HEOS","Methane&Nitrogen&Oxygen")
    );
    AS->set_mole_fractions({0.3, 0.5, 0.2});
    AS->specify_phase(CoolProp::iphase_gas);
    AS->update(CoolProp::DmolarT_INPUTS, 3000, 300);
    auto& HEOS = *dynamic_cast<CoolProp::HelmholtzEOSMixtureBackend*>(AS.get());
    auto xN  = CoolProp::x_N_dependency_flag::XN_INDEPENDENT;
    using md = CoolProp::MixtureDerivatives;
    std::cout << md::dalphar_dxi(HEOS, 0, xN) << std::endl;

    std::cout << md::d2alphar_dxi_dTau(HEOS, 0, xN)*AS->tau() << std::endl;
    std::cout << md::d2alphar_dxi_dDelta(HEOS, 0, xN)*AS->delta() << std::endl;
    std::cout << md::d3alphar_dxi_dTau2(HEOS, 0, xN)*pow(AS->tau(), 2) << std::endl;
    std::cout << md::d3alphar_dxi_dDelta_dTau(HEOS, 0, xN)*AS->tau()*AS->delta() << std::endl;
    std::cout << md::d3alphar_dxi_dDelta2(HEOS, 0, xN)*pow(AS->delta(), 2) << std::endl;

    std::cout << md::d3alphar_dxi_dxj_dTau(HEOS, 0, 1, xN)*AS->tau() << std::endl;
    std::cout << md::d3alphar_dxi_dxj_dDelta(HEOS, 0, 1, xN)*AS->delta() << std::endl;
    std::cout << md::d3alphardxidxjdxk(HEOS, 0, 1, 2, xN) << std::endl;
}

yields the output:

-0.0435874
-0.211886
-0.0365057
-0.0748886
-0.206939
0.0144689
-0.00597881
-0.00279186
0

which is the same as the above