Detail.cpp File Reference

Implementation of the AGA 8 Part 1 DETAIL equation of state for natural gas mixtures. More...

#include "Detail.h"
#include <math.h>
#include <cmath>
#include <iostream>
#include <cstdlib>

Include dependency graph for Detail.cpp:

Go to the source code of this file.

Functions
static void	xTermsDetail (const std::vector< double > &x)
static void	Alpha0Detail (const double T, const double D, const std::vector< double > &x, double a0[3])
	Calculates the ideal gas Helmholtz energy and its derivatives with respect to T and D.
static void	AlpharDetail (const int itau, const int idel, const double T, const double D, double ar[4][4])
	Calculates derivatives of the residual Helmholtz energy with respect to temperature and density.
double	sq (double x)
void	MolarMassDetail (const std::vector< double > &x, double &Mm)
	Calculate molar mass of a mixture based on composition.
void	PressureDetail (const double T, const double D, const std::vector< double > &x, double &P, double &Z)
	Calculates pressure and compressibility factor as a function of temperature and density.
void	DensityDetail (const double T, const double P, const std::vector< double > &x, double &D, int &ierr, std::string &herr)
	Calculates density as a function of temperature and pressure using an iterative method.
void	PropertiesDetail (const double T, const double D, const std::vector< double > &x, double &P, double &Z, double &dPdD, double &d2PdD2, double &d2PdTD, double &dPdT, double &U, double &H, double &S, double &Cv, double &Cp, double &W, double &G, double &JT, double &Kappa, double &Cf)
	Calculates thermodynamic properties as a function of temperature and density.
void	SetupDetail ()
	The following routine must be called once before any other routine.

Variables
static double	RDetail
static const int	NcDetail = 21
static const int	MaxFlds = 21
static const int	NTerms = 58
static const double	epsilon = 1e-15
static int	fn [NTerms+1]
static int	gn [NTerms+1]
static int	qn [NTerms+1]
static double	an [NTerms+1]
static double	un [NTerms+1]
static int	bn [NTerms+1]
static int	kn [NTerms+1]
static double	Bsnij2 [MaxFlds+1][MaxFlds+1][18+1]
static double	Bs [18+1]
static double	Csn [NTerms+1]
static double	Fi [MaxFlds+1]
static double	Gi [MaxFlds+1]
static double	Qi [MaxFlds+1]
static double	Ki25 [MaxFlds+1]
static double	Ei25 [MaxFlds+1]
static double	Kij5 [MaxFlds+1][MaxFlds+1]
static double	Uij5 [MaxFlds+1][MaxFlds+1]
static double	Gij5 [MaxFlds+1][MaxFlds+1]
static double	Tun [NTerms+1]
static double	Told
static double	n0i [MaxFlds+1][7+1]
static double	th0i [MaxFlds+1][7+1]
static double	MMiDetail [MaxFlds+1]
static double	K3
static double	xold [MaxFlds+1]
static double	dPdDsave

Detailed Description

Implementation of the AGA 8 Part 1 DETAIL equation of state for natural gas mixtures.

This implementation is based on the AGA Report No. 8, Part 1 - DETAIL method for the calculation of thermodynamic properties of natural gas mixtures. The code provides routines for calculating:

• Pressure from temperature and density
• Density from temperature and pressure (iterative)
• Various thermodynamic properties

The mixture compositions are specified using mole fractions for 21 components in the following order:

1. Methane 8. Isopentane 15. Hydrogen
2. Nitrogen 9. n-Pentane 16. Oxygen
3. Carbon dioxide 10. n-Hexane 17. Carbon monoxide
4. Ethane 11. n-Heptane 18. Water
5. Propane 12. n-Octane 19. Hydrogen sulfide
6. Isobutane 13. n-Nonane 20. Helium
7. n-Butane 14. n-Decane 21. Argon

Note

Before using any other functions, SetupDetail() must be called once to initialize constants

Author

Eric W. Lemmon (Original FORTRAN)

Ian H. Bell (C++ Translation)

Volker Heinemann, RMG Messtechnik GmbH (Contributor)

Jason Lu, Thermo Fisher Scientific (Contributor)

Copyright

Public Domain - developed by NIST employees

Version

2.0 (April 2017)

Definition in file Detail.cpp.

Function Documentation

Alpha0Detail()

static void Alpha0Detail ( const double T, const double D, const std::vector< double > & x, double a0[3] )

static

Calculates the ideal gas Helmholtz energy and its derivatives with respect to T and D.

This function computes the ideal gas Helmholtz energy and its first and second temperature derivatives. This calculation is not required when only pressure (P) or compressibility factor (Z) is needed.

Parameters

T	Temperature in Kelvin (K)
D	Density in moles per liter (mol/l)
x	Vector of mole fractions representing composition
a0	Array to store results: a0[0]: Ideal gas Helmholtz energy (J/mol) a0[1]: First temperature derivative ∂(a0)/∂T [J/(mol-K)] a0[2]: Second temperature derivative T*∂²(a0)/∂T² [J/(mol-K)]

Note

The function uses hyperbolic functions and logarithmic terms to compute the ideal gas contributions for each component in the mixture.

Warning

Input array x must be of sufficient size to handle all components (NcDetail)

Output array a0 must be of size 3 to store all computed values

Definition at line 550 of file Detail.cpp.

{
    // Private Sub Alpha0Detail(T, D, x, a0)
 
    // Calculate the ideal gas Helmholtz energy and its derivatives with respect to T and D.
    // This routine is not needed when only P (or Z) is calculated.
 
    // Inputs:
    //      T - Temperature (K)
    //      D - Density (mol/l)
    //    x() - Composition (mole fraction)
 
    // Outputs:
    // a0(0) - Ideal gas Helmholtz energy (J/mol)
    // a0(1) -   partial  (a0)/partial(T) [J/(mol-K)]
    // a0(2) - T*partial^2(a0)/partial(T)^2 [J/(mol-K)]
 
    double LogT, LogD, LogHyp, th0T, LogxD;
    double SumHyp0, SumHyp1, SumHyp2;
    double em, ep, hcn, hsn;
 
    a0[0] = 0;
    a0[1] = 0;
    a0[2] = 0;
    if (D > epsilon)
    {
        LogD = log(D);
    }
    else
    {
        LogD = log(epsilon);
    }
    LogT = log(T);
    for (std::size_t i = 1; i <= NcDetail; ++i)
    {
        if (x[i] > 0)
        {
            LogxD = LogD + log(x[i]);
            SumHyp0 = 0;
            SumHyp1 = 0;
            SumHyp2 = 0;
            for (int j = 4; j <= 7; ++j)
            {
                if (th0i[i][j] > 0)
                {
                    th0T = th0i[i][j] / T;
                    ep = exp(th0T);
                    em = 1 / ep;
                    hsn = (ep - em) / 2;
                    hcn = (ep + em) / 2;
                    if (j == 4 || j == 6)
                    {
                        LogHyp = log(std::abs(hsn));
                        SumHyp0 += n0i[i][j] * LogHyp;
                        SumHyp1 += n0i[i][j] * (LogHyp - th0T * hcn / hsn);
                        SumHyp2 += n0i[i][j] * sq(th0T / hsn);
                    }
                    else
                    {
                        LogHyp = log(std::abs(hcn));
                        SumHyp0 += -n0i[i][j] * LogHyp;
                        SumHyp1 += -n0i[i][j] * (LogHyp - th0T * hsn / hcn);
                        SumHyp2 += +n0i[i][j] * sq(th0T / hcn);
                    }
                }
            }
            a0[0] += x[i] * (LogxD + n0i[i][1] + n0i[i][2] / T - n0i[i][3] * LogT + SumHyp0);
            a0[1] += x[i] * (LogxD + n0i[i][1] - n0i[i][3] * (1 + LogT) + SumHyp1);
            a0[2] += -x[i] * (n0i[i][3] + SumHyp2);
        }
    }
    a0[0] = a0[0] * RDetail * T;
    a0[1] = a0[1] * RDetail;
    a0[2] = a0[2] * RDetail;
}

References epsilon, n0i, NcDetail, RDetail, sq(), and th0i.

Referenced by PropertiesDetail().

AlpharDetail()

static void AlpharDetail ( const int itau, const int idel, const double T, const double D, double ar[4][4] )

static

Calculates derivatives of the residual Helmholtz energy with respect to temperature and density.

This function computes various derivatives of the residual Helmholtz energy based on temperature and density inputs. The xTerms subroutine must be called before this routine if x has changed.

Parameters

itau	Set to 1 to calculate derivatives with respect to T [ar(1,0), ar(1,1), ar(2,0)], 0 otherwise
idel	Currently not used, reserved for future use in specifying highest density derivative
T	Temperature in Kelvin
D	Density in mol/l
ar[4][4]	Output array containing the following derivatives: ar[0][0]: Residual Helmholtz energy (J/mol) ar[0][1]: D*∂(ar)/∂D (J/mol) ar[0][2]: D²*∂²(ar)/∂D² (J/mol) ar[0][3]: D³*∂³(ar)/∂D³ (J/mol) ar[1][0]: ∂(ar)/∂T [J/(mol-K)] ar[1][1]: D*∂²(ar)/∂D∂T [J/(mol-K)] ar[2][0]: T*∂²(ar)/∂T² [J/(mol-K)]

Note

This function is part of the GERG-2008 equation of state calculations

The function assumes that global variables Told, Tun, K3, RDetail, un, Bs, Csn, bn, kn are properly initialized

Definition at line 649 of file Detail.cpp.

{
    // Private Sub AlpharDetail(itau, idel, T, D, ar)
 
    // Calculate the derivatives of the residual Helmholtz energy (ar) with respect to T and D.
    // itau and idel are inputs that contain the highest derivatives needed.
    // Outputs are returned in the array ar.
    // Subroutine xTerms must be called before this routine if x has changed
 
    // Inputs:
    //  itau - Set this to 1 to calculate "ar" derivatives with respect to T [i.e., ar(1,0), ar(1,1), and ar(2,0)], otherwise set it to 0.
    //  idel - Currently not used, but kept as an input for future use in specifing the highest density derivative needed.
    //     T - Temperature (K)
    //     D - Density (mol/l)
 
    // Outputs:
    // ar(0,0) - Residual Helmholtz energy (J/mol)
    // ar(0,1) -   D*partial  (ar)/partial(D) (J/mol)
    // ar(0,2) - D^2*partial^2(ar)/partial(D)^2 (J/mol)
    // ar(0,3) - D^3*partial^3(ar)/partial(D)^3 (J/mol)
    // ar(1,0) -     partial  (ar)/partial(T) [J/(mol-K)]
    // ar(1,1) -   D*partial^2(ar)/partial(D)/partial(T) [J/(mol-K)]
    // ar(2,0) -   T*partial^2(ar)/partial(T)^2 [J/(mol-K)]
 
    double ckd, bkd, Dred;
    double Sum, s0, s1, s2, s3, RT;
    double Sum0[NTerms + 1], SumB[NTerms + 1], Dknn[9 + 1], Expn[4 + 1];
    double CoefD1[NTerms + 1], CoefD2[NTerms + 1], CoefD3[NTerms + 1];
    double CoefT1[NTerms + 1], CoefT2[NTerms + 1];
 
    for (int i = 0; i <= 3; ++i)
    {
        for (int j = 0; j <= 3; ++j)
        {
            ar[i][j] = 0;
        }
    }
    if (std::abs(T - Told) > 0.0000001)
    {
        for (int n = 1; n <= 58; ++n)
        {
            Tun[n] = pow(T, -un[n]);
        }
    }
    Told = T;
 
    // Precalculation of common powers and exponents of density
    Dred = K3 * D;
    Dknn[0] = 1;
    for (int n = 1; n <= 9; ++n)
    {
        Dknn[n] = Dred * Dknn[n - 1];
    }
    Expn[0] = 1;
    for (int n = 1; n <= 4; ++n)
    {
        Expn[n] = exp(-Dknn[n]);
    }
    RT = RDetail * T;
 
    for (int n = 1; n <= 58; ++n)
    {
        // Contributions to the Helmholtz energy and its derivatives with respect to temperature
        CoefT1[n] = RDetail * (un[n] - 1);
        CoefT2[n] = CoefT1[n] * un[n];
        // Contributions to the virial coefficients
        SumB[n] = 0;
        Sum0[n] = 0;
        if (n <= 18)
        {
            Sum = Bs[n] * D;
            if (n >= 13)
            {
                Sum += -Csn[n] * Dred;
            }
            SumB[n] = Sum * Tun[n];
        }
        if (n >= 13)
        {
            // Contributions to the residual part of the Helmholtz energy
            Sum0[n] = Csn[n] * Dknn[bn[n]] * Tun[n] * Expn[kn[n]];
            // Contributions to the derivatives of the Helmholtz energy with respect to density
            bkd = bn[n] - kn[n] * Dknn[kn[n]];
            ckd = kn[n] * kn[n] * Dknn[kn[n]];
            CoefD1[n] = bkd;
            CoefD2[n] = bkd * (bkd - 1) - ckd;
            CoefD3[n] = (bkd - 2) * CoefD2[n] + ckd * (1 - kn[n] - 2 * bkd);
        }
        else
        {
            CoefD1[n] = 0;
            CoefD2[n] = 0;
            CoefD3[n] = 0;
        }
    }
 
    for (int n = 1; n <= 58; ++n)
    {
        // Density derivatives
        s0 = Sum0[n] + SumB[n];
        s1 = Sum0[n] * CoefD1[n] + SumB[n];
        s2 = Sum0[n] * CoefD2[n];
        s3 = Sum0[n] * CoefD3[n];
        ar[0][0] = ar[0][0] + RT * s0;
        ar[0][1] = ar[0][1] + RT * s1;
        ar[0][2] = ar[0][2] + RT * s2;
        ar[0][3] = ar[0][3] + RT * s3;
        // Temperature derivatives
        if (itau > 0)
        {
            ar[1][0] = ar[1][0] - CoefT1[n] * s0;
            ar[1][1] = ar[1][1] - CoefT1[n] * s1;
            ar[2][0] = ar[2][0] + CoefT2[n] * s0;
            // The following are not used, but fully functional
            // ar(1, 2) = ar(1, 2) - CoefT1(n) * s2;
            // ar(1, 3) = ar(1, 3) - CoefT1(n) * s3;
            // ar(2, 1) = ar(2, 1) + CoefT2(n) * s1;
            // ar(2, 2) = ar(2, 2) + CoefT2(n) * s2;
            // ar(2, 3) = ar(2, 3) + CoefT2(n) * s3;
        }
    }
}

References bn, Bs, Csn, K3, kn, NTerms, RDetail, Told, Tun, and un.

Referenced by PressureDetail(), and PropertiesDetail().

DensityDetail()

void DensityDetail	(	const double	T,
		const double	P,
		const std::vector< double > &	x,
		double &	D,
		int &	ierr,
		std::string &	herr )

Calculates density as a function of temperature and pressure using an iterative method.

This function uses an iterative Newton's method that calls PressureDetail to find the correct state point. Generally only 6 iterations at most are required. If the iteration fails to converge, the ideal gas density and an error message are returned. No checks are made to determine the phase boundary, which would have guaranteed that the output is in the gas phase. It is up to the user to locate the phase boundary, and thus identify the phase of the T and P inputs. If the state point is 2-phase, the output density will represent a metastable state.

Parameters

T	Temperature in Kelvin (K)
P	Pressure in kiloPascals (kPa)
x	Vector of mole fractions representing composition
D	Density in mol/l (can be negative to use as initial guess)
ierr	Error number (0 indicates no error)
herr	Error message string (empty if no error)

Note

If D is passed as negative, its absolute value will be used as initial estimate

If P is zero (or very close to zero), D will be set to 0

If calculation fails to converge, ideal gas density is returned with error message

See also

DensityDetail_wrapper for the Emscripten wrapped version of this function

Definition at line 246 of file Detail.cpp.

{
    // Sub DensityDetail(T, P, x, D, ierr, herr)
 
    // Calculate density as a function of temperature and pressure.  This is an iterative routine that calls PressureDetail
    // to find the correct state point.  Generally only 6 iterations at most are required.
    // If the iteration fails to converge, the ideal gas density and an error message are returned.
    // No checks are made to determine the phase boundary, which would have guaranteed that the output is in the gas phase.
    // It is up to the user to locate the phase boundary, and thus identify the phase of the T and P inputs.
    // If the state point is 2-phase, the output density will represent a metastable state.
 
    // Inputs:
    //      T - Temperature (K)
    //      P - Pressure (kPa)
    //    x() - Composition (mole fraction)
 
    // Outputs:
    //      D - Density (mol/l) (make D negative and send as an input to use as an initial guess)
    //   ierr - Error number (0 indicates no error)
    //   herr - Error message if ierr is not equal to zero
 
    double plog, vlog, P2, Z, dpdlv, vdiff, tolr;
 
    ierr = 0;
    herr = "";
    if (std::abs(P) < epsilon)
    {
        D = 0;
        return;
    }
    tolr = 0.0000001;
    if (D > -epsilon)
    {
        D = P / RDetail / T; // Ideal gas estimate
    }
    else
    {
        D = std::abs(D); // If D<0, then use as initial estimate
    }
    plog = log(P);
    vlog = -log(D);
    for (int it = 1; it <= 20; ++it)
    {
        if (vlog < -7 || vlog > 100)
        {
            ierr = 1;
            herr = "Calculation failed to converge in DETAIL method, ideal gas density returned.";
            D = P / RDetail / T;
            return;
        }
        D = exp(-vlog);
        PressureDetail(T, D, x, P2, Z);
        if (dPdDsave < epsilon || P2 < epsilon)
        {
            vlog += 0.1;
        }
        else
        {
            // Find the next density with a first order Newton's type iterative scheme, with
            // log(P) as the known variable and log(v) as the unknown property.
            // See AGA 8 publication for further information.
            dpdlv = -D * dPdDsave; // d(p)/d[log(v)]
            vdiff = (log(P2) - plog) * P2 / dpdlv;
            vlog = vlog - vdiff;
            if (std::abs(vdiff) < tolr)
            {
                D = exp(-vlog);
                return; // Iteration converged
            }
        }
    }
    ierr = 1;
    herr = "Calculation failed to converge in DETAIL method, ideal gas density returned.";
    D = P / RDetail / T;
    return;
}

References dPdDsave, epsilon, PressureDetail(), and RDetail.

Referenced by DensityDetail_wrapper().

MolarMassDetail()

void MolarMassDetail	(	const std::vector< double > &	x,
		double &	Mm )

Calculate molar mass of a mixture based on composition.

This function calculates the molar mass of a mixture using the mole fractions of each component provided in the composition array.

Parameters

	x	Vector of mole fractions for each component in the mixture. Must sum to 1.0. Use mole fractions only (not mass fractions or mole percents). Components must be ordered according to the fluid order defined in the implementation.
[out]	Mm	Calculated molar mass of the mixture in g/mol

Note

The composition array x must contain valid mole fractions that sum to 1.0

See also

MolarMassDetail_wrapper for the Emscripten wrapped version of this function

Definition at line 158 of file Detail.cpp.

{
    // Calculate molar mass of the mixture with the compositions contained in the x() input array
 
    // Inputs:
    //    x() - Composition (mole fraction)
    //          Do not send mole percents or mass fractions in the x() array, otherwise the output will be incorrect.
    //          The sum of the compositions in the x() array must be equal to one.
    //          The order of the fluids in this array is given at the top of this code.
 
    // Outputs:
    //     Mm - Molar mass (g/mol)
 
    Mm = 0;
    for (std::size_t i = 1; i <= NcDetail; ++i)
    {
        Mm += x[i] * MMiDetail[i];
    }
}

References MMiDetail, and NcDetail.

Referenced by MolarMassDetail_wrapper(), and PropertiesDetail().

PressureDetail()

void PressureDetail	(	const double	T,
		const double	D,
		const std::vector< double > &	x,
		double &	P,
		double &	Z )

Calculates pressure and compressibility factor as a function of temperature and density.

This function computes the pressure and compressibility factor for a gas mixture using the GERG-2008 equation of state. It also calculates d(P)/d(D) which is cached for use in iterative density calculations.

Parameters

	T	Temperature in Kelvin (K)
	D	Density in moles per liter (mol/l)
	x	Vector of composition mole fractions (must sum to 1.0)
[out]	P	Pressure in kiloPascals (kPa)
[out]	Z	Compressibility factor (dimensionless)

Note

The composition vector x must contain mole fractions, not mole percents or mass fractions

The sum of all mole fractions in x must equal 1.0

The derivative d(P)/d(D) is cached internally but not returned as an argument

See also

PressureDetail_wrapper for the Emscripten wrapped version of this function

Definition at line 196 of file Detail.cpp.

{
    // Sub Pressure(T, D, x, P, Z)
 
    // Calculate pressure as a function of temperature and density.  The derivative d(P)/d(D) is also calculated
    // for use in the iterative DensityDetail subroutine (and is only returned as a common variable).
 
    // Inputs:
    //      T - Temperature (K)
    //      D - Density (mol/l)
    //    x() - Composition (mole fraction)
    //          Do not send mole percents or mass fractions in the x() array, otherwise the output will be incorrect.
    //          The sum of the compositions in the x() array must be equal to one.
 
    // Outputs:
    //      P - Pressure (kPa)
    //      Z - Compressibility factor
    //   dPdDsave - d(P)/d(D) [kPa/(mol/l)] (at constant temperature)
    //            - This variable is cached in the common variables for use in the iterative density solver, but not returned as an argument.
 
    double ar[3 + 1][3 + 1];
    xTermsDetail(x);
    AlpharDetail(0, 2, T, D, ar);
    Z = 1 + ar[0][1] / RDetail / T; // ar(0,1) is the first derivative of alpha(r) with respect to density
    P = D * RDetail * T * Z;
    dPdDsave = RDetail * T + 2 * ar[0][1] + ar[0][2]; // d(P)/d(D) for use in density iteration
}

References AlpharDetail(), dPdDsave, RDetail, and xTermsDetail().

Referenced by DensityDetail(), and PressureDetail_wrapper().

PropertiesDetail()

void PropertiesDetail	(	const double	T,
		const double	D,
		const std::vector< double > &	x,
		double &	P,
		double &	Z,
		double &	dPdD,
		double &	d2PdD2,
		double &	d2PdTD,
		double &	dPdT,
		double &	U,
		double &	H,
		double &	S,
		double &	Cv,
		double &	Cp,
		double &	W,
		double &	G,
		double &	JT,
		double &	Kappa,
		double &	Cf )

Calculates thermodynamic properties as a function of temperature and density.

If density is unknown, call DensityDetail first with known pressure and temperature values. Makes calls to Molarmass, Alpha0Detail, and AlpharDetail subroutines.

Parameters

	T	Temperature in Kelvin (K)
	D	Density in mol/l
	x	Vector of mole fractions representing composition
[out]	P	Pressure in kPa
[out]	Z	Compressibility factor
[out]	dPdD	First derivative of pressure with respect to density at constant temperature [kPa/(mol/l)]
[out]	d2PdD2	Second derivative of pressure with respect to density at constant temperature [kPa/(mol/l)^2]
[out]	d2PdTD	Second derivative of pressure with respect to temperature and density [kPa/(mol/l)/K]
[out]	dPdT	First derivative of pressure with respect to temperature at constant density (kPa/K)
[out]	U	Internal energy in J/mol
[out]	H	Enthalpy in J/mol
[out]	S	Entropy in J/(mol-K)
[out]	Cv	Isochoric heat capacity in J/(mol-K)
[out]	Cp	Isobaric heat capacity in J/(mol-K)
[out]	W	Speed of sound in m/s
[out]	G	Gibbs energy in J/mol
[out]	JT	Joule-Thomson coefficient in K/kPa
[out]	Kappa	Isentropic Exponent
[out]	Cf	Critical Flow Factor (dimensionless)

See also

PropertiesDetail_wrapper for the Emscripten wrapped version of this function

Definition at line 350 of file Detail.cpp.

{
    // Sub Properties(T, D, x, P, Z, dPdD, d2PdD2, d2PdTD, dPdT, U, H, S, Cv, Cp, W, G, JT, Kappa)
 
    // Calculate thermodynamic properties as a function of temperature and density.  Calls are made to the subroutines
    // Molarmass, Alpha0Detail, and AlpharDetail.  If the density is not known, call subroutine DensityDetail first
    // with the known values of pressure and temperature.
 
    // Inputs:
    //      T - Temperature (K)
    //      D - Density (mol/l)
    //    x() - Composition (mole fraction)
 
    // Outputs:
    //      P - Pressure (kPa)
    //      Z - Compressibility factor
    //   dPdD - First derivative of pressure with respect to density at constant temperature [kPa/(mol/l)]
    // d2PdD2 - Second derivative of pressure with respect to density at constant temperature [kPa/(mol/l)^2]
    // d2PdTD - Second derivative of pressure with respect to temperature and density [kPa/(mol/l)/K] (currently not calculated)
    //   dPdT - First derivative of pressure with respect to temperature at constant density (kPa/K)
    //      U - Internal energy (J/mol)
    //      H - Enthalpy (J/mol)
    //      S - Entropy [J/(mol-K)]
    //     Cv - Isochoric heat capacity [J/(mol-K)]
    //     Cp - Isobaric heat capacity [J/(mol-K)]
    //      W - Speed of sound (m/s)
    //      G - Gibbs energy (J/mol)
    //     JT - Joule-Thomson coefficient (K/kPa)
    //  Kappa - Isentropic Exponent
    //     Cf - Critical Flow Factor (dimensionless)
 
    double a0[2 + 1], ar[3 + 1][3 + 1], Mm, A, R, RT;
 
    MolarMassDetail(x, Mm);
    xTermsDetail(x);
 
    // Calculate the ideal gas Helmholtz energy, and its first and second derivatives with respect to temperature.
    Alpha0Detail(T, D, x, a0);
 
    // Calculate the real gas Helmholtz energy, and its derivatives with respect to temperature and/or density.
    AlpharDetail(2, 3, T, D, ar);
 
    R = RDetail;
    RT = R * T;
    Z = 1 + ar[0][1] / RT;
    P = D * RT * Z;
    dPdD = RT + 2 * ar[0][1] + ar[0][2];
    dPdT = D * R + D * ar[1][1];
    A = a0[0] + ar[0][0];
    S = -a0[1] - ar[1][0];
    U = A + T * S;
    Cv = -(a0[2] + ar[2][0]);
    if (D > epsilon)
    {
        H = U + P / D;
        G = A + P / D;
        Cp = Cv + T * sq(dPdT / D) / dPdD;
        d2PdD2 = (2 * ar[0][1] + 4 * ar[0][2] + ar[0][3]) / D;
        JT = (T / D * dPdT / dPdD - 1) / Cp / D;
    }
    else
    {
        H = U + RT;
        G = A + RT;
        Cp = Cv + R;
        d2PdD2 = 0;
        JT = 1E+20; //=(dB/dT*T-B)/Cp for an ideal gas, but dB/dT is not calculated here
    }
    W = 1000 * Cp / Cv * dPdD / Mm;
    if (W < 0)
    {
        W = 0;
    }
    W = sqrt(W);
    Kappa = W * W * Mm / (RT * 1000 * Z);
    Cf = sqrt(Kappa * pow( (2 / (Kappa + 1)), ((Kappa + 1) / (Kappa - 1))));
    d2PdTD = 0;
}

References Alpha0Detail(), AlpharDetail(), epsilon, MolarMassDetail(), RDetail, sq(), and xTermsDetail().

Referenced by PropertiesDetail_wrapper().

SetupDetail()

void SetupDetail

(

)

The following routine must be called once before any other routine.

Initializes all constants and parameters for the DETAIL model equations of state

This function sets up all the necessary parameters for the AGA8 equation of state in its "DETAIL" version. It initializes:

• Gas constant (RDetail)
• Molar masses for 21 components (MMiDetail)
• Equation of state coefficients (an)
• Density exponents (bn)
• Temperature exponents (un)
• Various flags (fn, gn, qn, sn, wn)
• Binary interaction parameters:
  • Energy parameters (Eij)
  • Size parameters (Kij)
  • Orientation parameters (Gij)
  • Conformal energy parameters (Uij)
• Component-specific parameters:
  • Energy parameters (Ei)
  • Size parameters (Ki)
  • Orientation parameters (Gi)
  • Quadrupole parameters (Qi)
  • Dipole parameters (Si)
  • Association parameters (Wi)
• Ideal gas parameters (n0i, th0i)

The function also performs necessary precalculations of constants used in the equation of state calculations to optimize performance.

This implementation corresponds to the reference equations for natural gas mixtures as published in the DETAIL formulation.

Note

this function is directly bind to SetupDetail() in the Emscripten wrapper

Definition at line 806 of file Detail.cpp.

{
    // Initialize all the constants and parameters in the DETAIL model.
    // Some values are modified for calculations that do not depend on T, D, and x in order to speed up the program.
 
    int sn[NTerms + 1], wn[NTerms + 1];
    double Ei[MaxFlds + 1], Ki[MaxFlds + 1], Si[MaxFlds + 1], Wi[MaxFlds + 1], Bsnij;
    double Kij[MaxFlds + 1][MaxFlds + 1], Gij[MaxFlds + 1][MaxFlds + 1], Eij[MaxFlds + 1][MaxFlds + 1], Uij[MaxFlds + 1][MaxFlds + 1];
    double d0;
 
    RDetail = 8.31451;
 
    // Molar masses (g/mol)
    MMiDetail[1] = 16.043;   // Methane
    MMiDetail[2] = 28.0135;  // Nitrogen
    MMiDetail[3] = 44.01;    // Carbon dioxide
    MMiDetail[4] = 30.07;    // Ethane
    MMiDetail[5] = 44.097;   // Propane
    MMiDetail[6] = 58.123;   // Isobutane
    MMiDetail[7] = 58.123;   // n-Butane
    MMiDetail[8] = 72.15;    // Isopentane
    MMiDetail[9] = 72.15;    // n-Pentane
    MMiDetail[10] = 86.177;  // Hexane
    MMiDetail[11] = 100.204; // Heptane
    MMiDetail[12] = 114.231; // Octane
    MMiDetail[13] = 128.258; // Nonane
    MMiDetail[14] = 142.285; // Decane
    MMiDetail[15] = 2.0159;  // Hydrogen
    MMiDetail[16] = 31.9988; // Oxygen
    MMiDetail[17] = 28.01;   // Carbon monoxide
    MMiDetail[18] = 18.0153; // Water
    MMiDetail[19] = 34.082;  // Hydrogen sulfide
    MMiDetail[20] = 4.0026;  // Helium
    MMiDetail[21] = 39.948;  // Argon
 
    // Initialize constants
    Told = 0;
    for (int i = 1; i <= NTerms; ++i)
    {
        an[i] = 0;
        bn[i] = 0;
        gn[i] = 0;
        fn[i] = 0;
        kn[i] = 0;
        qn[i] = 0;
        sn[i] = 0;
        un[i] = 0;
        wn[i] = 0;
    }
 
    for (int i = 1; i <= MaxFlds; ++i)
    {
        Ei[i] = 0;
        Fi[i] = 0;
        Gi[i] = 0;
        Ki[i] = 0;
        Qi[i] = 0;
        Si[i] = 0;
        Wi[i] = 0;
        xold[i] = 0;
        for (int j = 1; j <= MaxFlds; ++j)
        {
            Eij[i][j] = 1;
            Gij[i][j] = 1;
            Kij[i][j] = 1;
            Uij[i][j] = 1;
        }
    }
 
    // Coefficients of the equation of state
    an[1] = 0.1538326;
    an[2] = 1.341953;
    an[3] = -2.998583;
    an[4] = -0.04831228;
    an[5] = 0.3757965;
    an[6] = -1.589575;
    an[7] = -0.05358847;
    an[8] = 0.88659463;
    an[9] = -0.71023704;
    an[10] = -1.471722;
    an[11] = 1.32185035;
    an[12] = -0.78665925;
    an[13] = 0.00000000229129;
    an[14] = 0.1576724;
    an[15] = -0.4363864;
    an[16] = -0.04408159;
    an[17] = -0.003433888;
    an[18] = 0.03205905;
    an[19] = 0.02487355;
    an[20] = 0.07332279;
    an[21] = -0.001600573;
    an[22] = 0.6424706;
    an[23] = -0.4162601;
    an[24] = -0.06689957;
    an[25] = 0.2791795;
    an[26] = -0.6966051;
    an[27] = -0.002860589;
    an[28] = -0.008098836;
    an[29] = 3.150547;
    an[30] = 0.007224479;
    an[31] = -0.7057529;
    an[32] = 0.5349792;
    an[33] = -0.07931491;
    an[34] = -1.418465;
    an[35] = -5.99905E-17;
    an[36] = 0.1058402;
    an[37] = 0.03431729;
    an[38] = -0.007022847;
    an[39] = 0.02495587;
    an[40] = 0.04296818;
    an[41] = 0.7465453;
    an[42] = -0.2919613;
    an[43] = 7.294616;
    an[44] = -9.936757;
    an[45] = -0.005399808;
    an[46] = -0.2432567;
    an[47] = 0.04987016;
    an[48] = 0.003733797;
    an[49] = 1.874951;
    an[50] = 0.002168144;
    an[51] = -0.6587164;
    an[52] = 0.000205518;
    an[53] = 0.009776195;
    an[54] = -0.02048708;
    an[55] = 0.01557322;
    an[56] = 0.006862415;
    an[57] = -0.001226752;
    an[58] = 0.002850908;
 
    // Density exponents
    bn[1] = 1;
    bn[2] = 1;
    bn[3] = 1;
    bn[4] = 1;
    bn[5] = 1;
    bn[6] = 1;
    bn[7] = 1;
    bn[8] = 1;
    bn[9] = 1;
    bn[10] = 1;
    bn[11] = 1;
    bn[12] = 1;
    bn[13] = 1;
    bn[14] = 1;
    bn[15] = 1;
    bn[16] = 1;
    bn[17] = 1;
    bn[18] = 1;
    bn[19] = 2;
    bn[20] = 2;
    bn[21] = 2;
    bn[22] = 2;
    bn[23] = 2;
    bn[24] = 2;
    bn[25] = 2;
    bn[26] = 2;
    bn[27] = 2;
    bn[28] = 3;
    bn[29] = 3;
    bn[30] = 3;
    bn[31] = 3;
    bn[32] = 3;
    bn[33] = 3;
    bn[34] = 3;
    bn[35] = 3;
    bn[36] = 3;
    bn[37] = 3;
    bn[38] = 4;
    bn[39] = 4;
    bn[40] = 4;
    bn[41] = 4;
    bn[42] = 4;
    bn[43] = 4;
    bn[44] = 4;
    bn[45] = 5;
    bn[46] = 5;
    bn[47] = 5;
    bn[48] = 5;
    bn[49] = 5;
    bn[50] = 6;
    bn[51] = 6;
    bn[52] = 7;
    bn[53] = 7;
    bn[54] = 8;
    bn[55] = 8;
    bn[56] = 8;
    bn[57] = 9;
    bn[58] = 9;
 
    // Exponents on density in EXP[-cn*D^kn] part
    // The cn part in this term is not included in this program since it is 1 when kn<>0][and 0 otherwise
    kn[13] = 3;
    kn[14] = 2;
    kn[15] = 2;
    kn[16] = 2;
    kn[17] = 4;
    kn[18] = 4;
    kn[21] = 2;
    kn[22] = 2;
    kn[23] = 2;
    kn[24] = 4;
    kn[25] = 4;
    kn[26] = 4;
    kn[27] = 4;
    kn[29] = 1;
    kn[30] = 1;
    kn[31] = 2;
    kn[32] = 2;
    kn[33] = 3;
    kn[34] = 3;
    kn[35] = 4;
    kn[36] = 4;
    kn[37] = 4;
    kn[40] = 2;
    kn[41] = 2;
    kn[42] = 2;
    kn[43] = 4;
    kn[44] = 4;
    kn[46] = 2;
    kn[47] = 2;
    kn[48] = 4;
    kn[49] = 4;
    kn[51] = 2;
    kn[53] = 2;
    kn[54] = 1;
    kn[55] = 2;
    kn[56] = 2;
    kn[57] = 2;
    kn[58] = 2;
 
    // Temperature exponents
    un[1] = 0;
    un[2] = 0.5;
    un[3] = 1;
    un[4] = 3.5;
    un[5] = -0.5;
    un[6] = 4.5;
    un[7] = 0.5;
    un[8] = 7.5;
    un[9] = 9.5;
    un[10] = 6;
    un[11] = 12;
    un[12] = 12.5;
    un[13] = -6;
    un[14] = 2;
    un[15] = 3;
    un[16] = 2;
    un[17] = 2;
    un[18] = 11;
    un[19] = -0.5;
    un[20] = 0.5;
    un[21] = 0;
    un[22] = 4;
    un[23] = 6;
    un[24] = 21;
    un[25] = 23;
    un[26] = 22;
    un[27] = -1;
    un[28] = -0.5;
    un[29] = 7;
    un[30] = -1;
    un[31] = 6;
    un[32] = 4;
    un[33] = 1;
    un[34] = 9;
    un[35] = -13;
    un[36] = 21;
    un[37] = 8;
    un[38] = -0.5;
    un[39] = 0;
    un[40] = 2;
    un[41] = 7;
    un[42] = 9;
    un[43] = 22;
    un[44] = 23;
    un[45] = 1;
    un[46] = 9;
    un[47] = 3;
    un[48] = 8;
    un[49] = 23;
    un[50] = 1.5;
    un[51] = 5;
    un[52] = -0.5;
    un[53] = 4;
    un[54] = 7;
    un[55] = 3;
    un[56] = 0;
    un[57] = 1;
    un[58] = 0;
 
    // Flags
    fn[13] = 1;
    fn[27] = 1;
    fn[30] = 1;
    fn[35] = 1;
    gn[5] = 1;
    gn[6] = 1;
    gn[25] = 1;
    gn[29] = 1;
    gn[32] = 1;
    gn[33] = 1;
    gn[34] = 1;
    gn[51] = 1;
    gn[54] = 1;
    gn[56] = 1;
    qn[7] = 1;
    qn[16] = 1;
    qn[26] = 1;
    qn[28] = 1;
    qn[37] = 1;
    qn[42] = 1;
    qn[47] = 1;
    qn[49] = 1;
    qn[52] = 1;
    qn[58] = 1;
    sn[8] = 1;
    sn[9] = 1;
    wn[10] = 1;
    wn[11] = 1;
    wn[12] = 1;
 
    // Energy parameters
    Ei[1] = 151.3183;
    Ei[2] = 99.73778;
    Ei[3] = 241.9606;
    Ei[4] = 244.1667;
    Ei[5] = 298.1183;
    Ei[6] = 324.0689;
    Ei[7] = 337.6389;
    Ei[8] = 365.5999;
    Ei[9] = 370.6823;
    Ei[10] = 402.636293;
    Ei[11] = 427.72263;
    Ei[12] = 450.325022;
    Ei[13] = 470.840891;
    Ei[14] = 489.558373;
    Ei[15] = 26.95794;
    Ei[16] = 122.7667;
    Ei[17] = 105.5348;
    Ei[18] = 514.0156;
    Ei[19] = 296.355;
    Ei[20] = 2.610111;
    Ei[21] = 119.6299;
 
    // Size parameters
    Ki[1] = 0.4619255;
    Ki[2] = 0.4479153;
    Ki[3] = 0.4557489;
    Ki[4] = 0.5279209;
    Ki[5] = 0.583749;
    Ki[6] = 0.6406937;
    Ki[7] = 0.6341423;
    Ki[8] = 0.6738577;
    Ki[9] = 0.6798307;
    Ki[10] = 0.7175118;
    Ki[11] = 0.7525189;
    Ki[12] = 0.784955;
    Ki[13] = 0.8152731;
    Ki[14] = 0.8437826;
    Ki[15] = 0.3514916;
    Ki[16] = 0.4186954;
    Ki[17] = 0.4533894;
    Ki[18] = 0.3825868;
    Ki[19] = 0.4618263;
    Ki[20] = 0.3589888;
    Ki[21] = 0.4216551;
 
    // Orientation parameters
    Gi[2] = 0.027815;
    Gi[3] = 0.189065;
    Gi[4] = 0.0793;
    Gi[5] = 0.141239;
    Gi[6] = 0.256692;
    Gi[7] = 0.281835;
    Gi[8] = 0.332267;
    Gi[9] = 0.366911;
    Gi[10] = 0.289731;
    Gi[11] = 0.337542;
    Gi[12] = 0.383381;
    Gi[13] = 0.427354;
    Gi[14] = 0.469659;
    Gi[15] = 0.034369;
    Gi[16] = 0.021;
    Gi[17] = 0.038953;
    Gi[18] = 0.3325;
    Gi[19] = 0.0885;
 
    // Quadrupole parameters
    Qi[3] = 0.69;
    Qi[18] = 1.06775;
    Qi[19] = 0.633276;
    Fi[15] = 1;      // High temperature parameter
    Si[18] = 1.5822; // Dipole parameter
    Si[19] = 0.39;   // Dipole parameter
    Wi[18] = 1;      // Association parameter
 
    // Energy parameters
    Eij[1][2] = 0.97164;
    Eij[1][3] = 0.960644;
    Eij[1][5] = 0.994635;
    Eij[1][6] = 1.01953;
    Eij[1][7] = 0.989844;
    Eij[1][8] = 1.00235;
    Eij[1][9] = 0.999268;
    Eij[1][10] = 1.107274;
    Eij[1][11] = 0.88088;
    Eij[1][12] = 0.880973;
    Eij[1][13] = 0.881067;
    Eij[1][14] = 0.881161;
    Eij[1][15] = 1.17052;
    Eij[1][17] = 0.990126;
    Eij[1][18] = 0.708218;
    Eij[1][19] = 0.931484;
    Eij[2][3] = 1.02274;
    Eij[2][4] = 0.97012;
    Eij[2][5] = 0.945939;
    Eij[2][6] = 0.946914;
    Eij[2][7] = 0.973384;
    Eij[2][8] = 0.95934;
    Eij[2][9] = 0.94552;
    Eij[2][15] = 1.08632;
    Eij[2][16] = 1.021;
    Eij[2][17] = 1.00571;
    Eij[2][18] = 0.746954;
    Eij[2][19] = 0.902271;
    Eij[3][4] = 0.925053;
    Eij[3][5] = 0.960237;
    Eij[3][6] = 0.906849;
    Eij[3][7] = 0.897362;
    Eij[3][8] = 0.726255;
    Eij[3][9] = 0.859764;
    Eij[3][10] = 0.855134;
    Eij[3][11] = 0.831229;
    Eij[3][12] = 0.80831;
    Eij[3][13] = 0.786323;
    Eij[3][14] = 0.765171;
    Eij[3][15] = 1.28179;
    Eij[3][17] = 1.5;
    Eij[3][18] = 0.849408;
    Eij[3][19] = 0.955052;
    Eij[4][5] = 1.02256;
    Eij[4][7] = 1.01306;
    Eij[4][9] = 1.00532;
    Eij[4][15] = 1.16446;
    Eij[4][18] = 0.693168;
    Eij[4][19] = 0.946871;
    Eij[5][7] = 1.0049;
    Eij[5][15] = 1.034787;
    Eij[6][15] = 1.3;
    Eij[7][15] = 1.3;
    Eij[10][19] = 1.008692;
    Eij[11][19] = 1.010126;
    Eij[12][19] = 1.011501;
    Eij[13][19] = 1.012821;
    Eij[14][19] = 1.014089;
    Eij[15][17] = 1.1;
 
    // Conformal energy parameters
    Uij[1][2] = 0.886106;
    Uij[1][3] = 0.963827;
    Uij[1][5] = 0.990877;
    Uij[1][7] = 0.992291;
    Uij[1][9] = 1.00367;
    Uij[1][10] = 1.302576;
    Uij[1][11] = 1.191904;
    Uij[1][12] = 1.205769;
    Uij[1][13] = 1.219634;
    Uij[1][14] = 1.233498;
    Uij[1][15] = 1.15639;
    Uij[1][19] = 0.736833;
    Uij[2][3] = 0.835058;
    Uij[2][4] = 0.816431;
    Uij[2][5] = 0.915502;
    Uij[2][7] = 0.993556;
    Uij[2][15] = 0.408838;
    Uij[2][19] = 0.993476;
    Uij[3][4] = 0.96987;
    Uij[3][10] = 1.066638;
    Uij[3][11] = 1.077634;
    Uij[3][12] = 1.088178;
    Uij[3][13] = 1.098291;
    Uij[3][14] = 1.108021;
    Uij[3][17] = 0.9;
    Uij[3][19] = 1.04529;
    Uij[4][5] = 1.065173;
    Uij[4][6] = 1.25;
    Uij[4][7] = 1.25;
    Uij[4][8] = 1.25;
    Uij[4][9] = 1.25;
    Uij[4][15] = 1.61666;
    Uij[4][19] = 0.971926;
    Uij[10][19] = 1.028973;
    Uij[11][19] = 1.033754;
    Uij[12][19] = 1.038338;
    Uij[13][19] = 1.042735;
    Uij[14][19] = 1.046966;
 
    // Size parameters
    Kij[1][2] = 1.00363;
    Kij[1][3] = 0.995933;
    Kij[1][5] = 1.007619;
    Kij[1][7] = 0.997596;
    Kij[1][9] = 1.002529;
    Kij[1][10] = 0.982962;
    Kij[1][11] = 0.983565;
    Kij[1][12] = 0.982707;
    Kij[1][13] = 0.981849;
    Kij[1][14] = 0.980991;
    Kij[1][15] = 1.02326;
    Kij[1][19] = 1.00008;
    Kij[2][3] = 0.982361;
    Kij[2][4] = 1.00796;
    Kij[2][15] = 1.03227;
    Kij[2][19] = 0.942596;
    Kij[3][4] = 1.00851;
    Kij[3][10] = 0.910183;
    Kij[3][11] = 0.895362;
    Kij[3][12] = 0.881152;
    Kij[3][13] = 0.86752;
    Kij[3][14] = 0.854406;
    Kij[3][19] = 1.00779;
    Kij[4][5] = 0.986893;
    Kij[4][15] = 1.02034;
    Kij[4][19] = 0.999969;
    Kij[10][19] = 0.96813;
    Kij[11][19] = 0.96287;
    Kij[12][19] = 0.957828;
    Kij[13][19] = 0.952441;
    Kij[14][19] = 0.948338;
 
    // Orientation parameters
    Gij[1][3] = 0.807653;
    Gij[1][15] = 1.95731;
    Gij[2][3] = 0.982746;
    Gij[3][4] = 0.370296;
    Gij[3][18] = 1.67309;
 
    // Ideal gas parameters
    n0i[1][3] = 4.00088;
    n0i[1][4] = 0.76315;
    n0i[1][5] = 0.0046;
    n0i[1][6] = 8.74432;
    n0i[1][7] = -4.46921;
    n0i[1][1] = 29.83843397;
    n0i[1][2] = -15999.69151;
    n0i[2][3] = 3.50031;
    n0i[2][4] = 0.13732;
    n0i[2][5] = -0.1466;
    n0i[2][6] = 0.90066;
    n0i[2][7] = 0;
    n0i[2][1] = 17.56770785;
    n0i[2][2] = -2801.729072;
    n0i[3][3] = 3.50002;
    n0i[3][4] = 2.04452;
    n0i[3][5] = -1.06044;
    n0i[3][6] = 2.03366;
    n0i[3][7] = 0.01393;
    n0i[3][1] = 20.65844696;
    n0i[3][2] = -4902.171516;
    n0i[4][3] = 4.00263;
    n0i[4][4] = 4.33939;
    n0i[4][5] = 1.23722;
    n0i[4][6] = 13.1974;
    n0i[4][7] = -6.01989;
    n0i[4][1] = 36.73005938;
    n0i[4][2] = -23639.65301;
    n0i[5][3] = 4.02939;
    n0i[5][4] = 6.60569;
    n0i[5][5] = 3.197;
    n0i[5][6] = 19.1921;
    n0i[5][7] = -8.37267;
    n0i[5][1] = 44.70909619;
    n0i[5][2] = -31236.63551;
    n0i[6][3] = 4.06714;
    n0i[6][4] = 8.97575;
    n0i[6][5] = 5.25156;
    n0i[6][6] = 25.1423;
    n0i[6][7] = 16.1388;
    n0i[6][1] = 34.30180349;
    n0i[6][2] = -38525.50276;
    n0i[7][3] = 4.33944;
    n0i[7][4] = 9.44893;
    n0i[7][5] = 6.89406;
    n0i[7][6] = 24.4618;
    n0i[7][7] = 14.7824;
    n0i[7][1] = 36.53237783;
    n0i[7][2] = -38957.80933;
    n0i[8][3] = 4;
    n0i[8][4] = 11.7618;
    n0i[8][5] = 20.1101;
    n0i[8][6] = 33.1688;
    n0i[8][7] = 0;
    n0i[8][1] = 43.17218626;
    n0i[8][2] = -51198.30946;
    n0i[9][3] = 4;
    n0i[9][4] = 8.95043;
    n0i[9][5] = 21.836;
    n0i[9][6] = 33.4032;
    n0i[9][7] = 0;
    n0i[9][1] = 42.67837089;
    n0i[9][2] = -45215.83;
    n0i[10][3] = 4;
    n0i[10][4] = 11.6977;
    n0i[10][5] = 26.8142;
    n0i[10][6] = 38.6164;
    n0i[10][7] = 0;
    n0i[10][1] = 46.99717188;
    n0i[10][2] = -52746.83318;
    n0i[11][3] = 4;
    n0i[11][4] = 13.7266;
    n0i[11][5] = 30.4707;
    n0i[11][6] = 43.5561;
    n0i[11][7] = 0;
    n0i[11][1] = 52.07631631;
    n0i[11][2] = -57104.81056;
    n0i[12][3] = 4;
    n0i[12][4] = 15.6865;
    n0i[12][5] = 33.8029;
    n0i[12][6] = 48.1731;
    n0i[12][7] = 0;
    n0i[12][1] = 57.25830934;
    n0i[12][2] = -60546.76385;
    n0i[13][3] = 4;
    n0i[13][4] = 18.0241;
    n0i[13][5] = 38.1235;
    n0i[13][6] = 53.3415;
    n0i[13][7] = 0;
    n0i[13][1] = 62.09646901;
    n0i[13][2] = -66600.12837;
    n0i[14][3] = 4;
    n0i[14][4] = 21.0069;
    n0i[14][5] = 43.4931;
    n0i[14][6] = 58.3657;
    n0i[14][7] = 0;
    n0i[14][1] = 65.93909154;
    n0i[14][2] = -74131.45483;
    n0i[15][3] = 2.47906;
    n0i[15][4] = 0.95806;
    n0i[15][5] = 0.45444;
    n0i[15][6] = 1.56039;
    n0i[15][7] = -1.3756;
    n0i[15][1] = 13.07520288;
    n0i[15][2] = -5836.943696;
    n0i[16][3] = 3.50146;
    n0i[16][4] = 1.07558;
    n0i[16][5] = 1.01334;
    n0i[16][6] = 0;
    n0i[16][7] = 0;
    n0i[16][1] = 16.8017173;
    n0i[16][2] = -2318.32269;
    n0i[17][3] = 3.50055;
    n0i[17][4] = 1.02865;
    n0i[17][5] = 0.00493;
    n0i[17][6] = 0;
    n0i[17][7] = 0;
    n0i[17][1] = 17.45786899;
    n0i[17][2] = -2635.244116;
    n0i[18][3] = 4.00392;
    n0i[18][4] = 0.01059;
    n0i[18][5] = 0.98763;
    n0i[18][6] = 3.06904;
    n0i[18][7] = 0;
    n0i[18][1] = 21.57882705;
    n0i[18][2] = -7766.733078;
    n0i[19][3] = 4;
    n0i[19][4] = 3.11942;
    n0i[19][5] = 1.00243;
    n0i[19][6] = 0;
    n0i[19][7] = 0;
    n0i[19][1] = 21.5830944;
    n0i[19][2] = -6069.035869;
    n0i[20][3] = 2.5;
    n0i[20][4] = 0;
    n0i[20][5] = 0;
    n0i[20][6] = 0;
    n0i[20][7] = 0;
    n0i[20][1] = 10.04639507;
    n0i[20][2] = -745.375;
    n0i[21][3] = 2.5;
    n0i[21][4] = 0;
    n0i[21][5] = 0;
    n0i[21][6] = 0;
    n0i[21][7] = 0;
    n0i[21][1] = 10.04639507;
    n0i[21][2] = -745.375;
    th0i[1][4] = 820.659;
    th0i[1][5] = 178.41;
    th0i[1][6] = 1062.82;
    th0i[1][7] = 1090.53;
    th0i[2][4] = 662.738;
    th0i[2][5] = 680.562;
    th0i[2][6] = 1740.06;
    th0i[2][7] = 0;
    th0i[3][4] = 919.306;
    th0i[3][5] = 865.07;
    th0i[3][6] = 483.553;
    th0i[3][7] = 341.109;
    th0i[4][4] = 559.314;
    th0i[4][5] = 223.284;
    th0i[4][6] = 1031.38;
    th0i[4][7] = 1071.29;
    th0i[5][4] = 479.856;
    th0i[5][5] = 200.893;
    th0i[5][6] = 955.312;
    th0i[5][7] = 1027.29;
    th0i[6][4] = 438.27;
    th0i[6][5] = 198.018;
    th0i[6][6] = 1905.02;
    th0i[6][7] = 893.765;
    th0i[7][4] = 468.27;
    th0i[7][5] = 183.636;
    th0i[7][6] = 1914.1;
    th0i[7][7] = 903.185;
    th0i[8][4] = 292.503;
    th0i[8][5] = 910.237;
    th0i[8][6] = 1919.37;
    th0i[8][7] = 0;
    th0i[9][4] = 178.67;
    th0i[9][5] = 840.538;
    th0i[9][6] = 1774.25;
    th0i[9][7] = 0;
    th0i[10][4] = 182.326;
    th0i[10][5] = 859.207;
    th0i[10][6] = 1826.59;
    th0i[10][7] = 0;
    th0i[11][4] = 169.789;
    th0i[11][5] = 836.195;
    th0i[11][6] = 1760.46;
    th0i[11][7] = 0;
    th0i[12][4] = 158.922;
    th0i[12][5] = 815.064;
    th0i[12][6] = 1693.07;
    th0i[12][7] = 0;
    th0i[13][4] = 156.854;
    th0i[13][5] = 814.882;
    th0i[13][6] = 1693.79;
    th0i[13][7] = 0;
    th0i[14][4] = 164.947;
    th0i[14][5] = 836.264;
    th0i[14][6] = 1750.24;
    th0i[14][7] = 0;
    th0i[15][4] = 228.734;
    th0i[15][5] = 326.843;
    th0i[15][6] = 1651.71;
    th0i[15][7] = 1671.69;
    th0i[16][4] = 2235.71;
    th0i[16][5] = 1116.69;
    th0i[16][6] = 0;
    th0i[16][7] = 0;
    th0i[17][4] = 1550.45;
    th0i[17][5] = 704.525;
    th0i[17][6] = 0;
    th0i[17][7] = 0;
    th0i[18][4] = 268.795;
    th0i[18][5] = 1141.41;
    th0i[18][6] = 2507.37;
    th0i[18][7] = 0;
    th0i[19][4] = 1833.63;
    th0i[19][5] = 847.181;
    th0i[19][6] = 0;
    th0i[19][7] = 0;
    th0i[20][4] = 0;
    th0i[20][5] = 0;
    th0i[20][6] = 0;
    th0i[20][7] = 0;
    th0i[21][4] = 0;
    th0i[21][5] = 0;
    th0i[21][6] = 0;
    th0i[21][7] = 0;
 
    // Precalculations of constants
    for (int i = 1; i <= MaxFlds; ++i)
    {
        Ki25[i] = pow(Ki[i], 2.5);
        Ei25[i] = pow(Ei[i], 2.5);
    }
    for (int i = 1; i <= MaxFlds; ++i)
    {
        for (int j = i; j <= MaxFlds; ++j)
        {
            for (int n = 1; n <= 18; ++n)
            {
                Bsnij = 1;
                if (gn[n] == 1)
                {
                    Bsnij = Gij[i][j] * (Gi[i] + Gi[j]) / 2;
                }
                if (qn[n] == 1)
                {
                    Bsnij = Bsnij * Qi[i] * Qi[j];
                }
                if (fn[n] == 1)
                {
                    Bsnij = Bsnij * Fi[i] * Fi[j];
                }
                if (sn[n] == 1)
                {
                    Bsnij = Bsnij * Si[i] * Si[j];
                }
                if (wn[n] == 1)
                {
                    Bsnij = Bsnij * Wi[i] * Wi[j];
                }
                Bsnij2[i][j][n] = an[n] * pow(Eij[i][j] * sqrt(Ei[i] * Ei[j]), un[n]) * pow(Ki[i] * Ki[j], 1.5) * Bsnij;
            }
            Kij5[i][j] = (pow(Kij[i][j], 5) - 1) * Ki25[i] * Ki25[j];
            Uij5[i][j] = (pow(Uij[i][j], 5) - 1) * Ei25[i] * Ei25[j];
            Gij5[i][j] = (Gij[i][j] - 1) * (Gi[i] + Gi[j]) / 2;
        }
    }
    // Ideal gas terms
    d0 = 101.325 / RDetail / 298.15;
    for (int i = 1; i <= MaxFlds; ++i)
    {
        n0i[i][3] = n0i[i][3] - 1;
        n0i[i][1] = n0i[i][1] - log(d0);
    }
    return;
 
    // Code to produce nearly exact values for n0[1] and n0[2]
    // This is not called in the current code, but included below to show how the values were calculated.  The return above can be removed to call this code.
    // T0 = 298.15;
    // d0 = 101.325 / RDetail / T0;
    // for (int i=1; i <= MaxFlds; ++i){
    //   n1 = 0; n2 = 0;
    //   if (th0i[i][4] > 0) {n2 = n2 - n0i[i][4] * th0i[i][4] / tanh(th0i[i][4] / T0); n1 = n1 - n0i[i][4] * log(sinh(th0i[i][4] / T0));}
    //   if (th0i[i][5] > 0) {n2 = n2 + n0i[i][5] * th0i[i][5] * tanh(th0i[i][5] / T0); n1 = n1 + n0i[i][5] * log(cosh(th0i[i][5] / T0));}
    //   if (th0i[i][6] > 0) {n2 = n2 - n0i[i][6] * th0i[i][6] / tanh(th0i[i][6] / T0); n1 = n1 - n0i[i][6] * log(sinh(th0i[i][6] / T0));}
    //   if (th0i[i][7] > 0) {n2 = n2 + n0i[i][7] * th0i[i][7] * tanh(th0i[i][7] / T0); n1 = n1 + n0i[i][7] * log(cosh(th0i[i][7] / T0));}
    //   n0i[i][2] = n2 - n0i[i][3] * T0;
    //   n0i[i][3] = n0i[i][3] - 1;
    //   n0i[i][1] = n1 - n2 / T0 + n0i[i][3] * (1 + log(T0)) - log(d0);
    // }
}

References an, bn, Bsnij2, Ei25, Fi, fn, Gi, Gij5, gn, Ki25, Kij5, kn, MaxFlds, MMiDetail, n0i, NTerms, Qi, qn, RDetail, th0i, Told, Uij5, un, and xold.

Referenced by EMSCRIPTEN_BINDINGS().

sq()

double sq ( double x )

inline

Definition at line 142 of file Detail.cpp.

{ return x * x; }

Referenced by Alpha0Detail(), PropertiesDetail(), and xTermsDetail().

xTermsDetail()

static void xTermsDetail ( const std::vector< double > & x )

static

Definition at line 430 of file Detail.cpp.

{
    // Calculate terms dependent only on composition
    //
    // Inputs:
    //    x() - Composition (mole fraction)
 
    double G, Q, F, U, Q2, xij, xi2;
    int icheck;
 
    // Check to see if a component fraction has changed.  If x is the same as the previous call, then exit.
    icheck = 0;
    for (std::size_t i = 1; i <= NcDetail; ++i)
    {
        if (std::abs(x[i] - xold[i]) > 0.0000001)
        {
            icheck = 1;
        }
        xold[i] = x[i];
    }
    if (icheck == 0)
    {
        return;
    }
 
    K3 = 0;
    U = 0;
    G = 0;
    Q = 0;
    F = 0;
    for (int n = 1; n <= 18; ++n)
    {
        Bs[n] = 0;
    }
 
    // Calculate pure fluid contributions
    for (std::size_t i = 1; i <= NcDetail; ++i)
    {
        if (x[i] > 0)
        {
            xi2 = sq(x[i]);
            K3 += x[i] * Ki25[i]; // K, U, and G are the sums of a pure fluid contribution and a
            U += x[i] * Ei25[i];  // binary pair contribution
            G += x[i] * Gi[i];
            Q += x[i] * Qi[i]; // Q and F depend only on the pure fluid parts
            F += xi2 * Fi[i];
            for (int n = 1; n <= 18; ++n)
            {
                Bs[n] = Bs[n] + xi2 * Bsnij2[i][i][n]; // Pure fluid contributions to second virial coefficient
            }
        }
    }
    K3 = sq(K3);
    U = sq(U);
 
    // Binary pair contributions
    for (std::size_t i = 1; i <= NcDetail - 1; ++i)
    {
        if (x[i] > 0)
        {
            for (std::size_t j = i + 1; j <= NcDetail; ++j)
            {
                if (x[j] > 0)
                {
                    xij = 2 * x[i] * x[j];
                    K3 = K3 + xij * Kij5[i][j];
                    U = U + xij * Uij5[i][j];
                    G = G + xij * Gij5[i][j];
                    for (int n = 1; n <= 18; ++n)
                    {
                        Bs[n] = Bs[n] + xij * Bsnij2[i][j][n]; // Second virial coefficients of mixture
                    }
                }
            }
        }
    }
    K3 = pow(K3, 0.6);
    U = pow(U, 0.2);
 
    // Third virial and higher coefficients
    Q2 = sq(Q);
    for (int n = 13; n <= 58; ++n)
    {
        Csn[n] = an[n] * pow(U, un[n]);
        if (gn[n] == 1)
        {
            Csn[n] = Csn[n] * G;
        }
        if (qn[n] == 1)
        {
            Csn[n] = Csn[n] * Q2;
        }
        if (fn[n] == 1)
        {
            Csn[n] = Csn[n] * F;
        }
    }
}

References an, Bs, Bsnij2, Csn, Ei25, Fi, fn, Gi, Gij5, gn, K3, Ki25, Kij5, NcDetail, Qi, qn, sq(), Uij5, un, and xold.

Referenced by PressureDetail(), and PropertiesDetail().

Variable Documentation

an

double an[NTerms+1]

static

Definition at line 131 of file Detail.cpp.

Referenced by SetupDetail(), and xTermsDetail().

bn

int bn[NTerms+1]

static

Definition at line 132 of file Detail.cpp.

Referenced by AlpharDetail(), and SetupDetail().

Bs

double Bs[18+1]

static

Definition at line 133 of file Detail.cpp.

Referenced by AlpharDetail(), and xTermsDetail().

Bsnij2

double Bsnij2[MaxFlds+1][MaxFlds+1][18+1]

static

Definition at line 133 of file Detail.cpp.

Referenced by SetupDetail(), and xTermsDetail().

Csn

double Csn[NTerms+1]

static

Definition at line 133 of file Detail.cpp.

Referenced by AlpharDetail(), and xTermsDetail().

dPdDsave

double dPdDsave

static

Definition at line 140 of file Detail.cpp.

Referenced by DensityDetail(), DensityGERG(), DensityGross(), PressureDetail(), PressureGERG(), and PressureGross().

Ei25

double Ei25[MaxFlds+1]

static

Definition at line 135 of file Detail.cpp.

Referenced by SetupDetail(), and xTermsDetail().

epsilon

const double epsilon = 1e-15

static

Definition at line 129 of file Detail.cpp.

Referenced by Alpha0Detail(), Alpha0GERG(), AlpharGERG(), DensityDetail(), DensityGERG(), DensityGross(), GrossMethod1(), GrossMethod2(), PropertiesDetail(), PropertiesGERG(), PseudoCriticalPointGERG(), ReducingParametersGERG(), and tTermsGERG().

Fi

double Fi[MaxFlds+1]

static

Definition at line 134 of file Detail.cpp.

Referenced by SetupDetail(), and xTermsDetail().

fn

int fn[NTerms+1]

static

Definition at line 130 of file Detail.cpp.

Referenced by SetupDetail(), and xTermsDetail().

Gi

double Gi[MaxFlds+1]

static

Definition at line 134 of file Detail.cpp.

Referenced by SetupDetail(), and xTermsDetail().

Gij5

double Gij5[MaxFlds+1][MaxFlds+1]

static

Definition at line 136 of file Detail.cpp.

Referenced by SetupDetail(), and xTermsDetail().

gn

int gn[NTerms+1]

static

Definition at line 130 of file Detail.cpp.

Referenced by SetupDetail(), and xTermsDetail().

K3

double K3

static

Definition at line 139 of file Detail.cpp.

Referenced by AlpharDetail(), and xTermsDetail().

Ki25

double Ki25[MaxFlds+1]

static

Definition at line 135 of file Detail.cpp.

Referenced by SetupDetail(), and xTermsDetail().

Kij5

double Kij5[MaxFlds+1][MaxFlds+1]

static

Definition at line 136 of file Detail.cpp.

Referenced by SetupDetail(), and xTermsDetail().

kn

int kn[NTerms+1]

static

Definition at line 132 of file Detail.cpp.

Referenced by AlpharDetail(), and SetupDetail().

MaxFlds

const int MaxFlds = 21

static

Definition at line 128 of file Detail.cpp.

Referenced by SetupDetail(), and SetupGERG().

MMiDetail

double MMiDetail[MaxFlds+1]

static

Definition at line 139 of file Detail.cpp.

Referenced by MolarMassDetail(), and SetupDetail().

n0i

double n0i[MaxFlds+1][7+1]

static

Definition at line 138 of file Detail.cpp.

Referenced by Alpha0Detail(), Alpha0GERG(), SetupDetail(), and SetupGERG().

NcDetail

const int NcDetail = 21

static

Definition at line 128 of file Detail.cpp.

Referenced by Alpha0Detail(), MolarMassDetail(), and xTermsDetail().

NTerms

const int NTerms = 58

static

Definition at line 128 of file Detail.cpp.

Referenced by AlpharDetail(), and SetupDetail().

Qi

double Qi[MaxFlds+1]

static

Definition at line 134 of file Detail.cpp.

Referenced by SetupDetail(), and xTermsDetail().

qn

int qn[NTerms+1]

static

Definition at line 130 of file Detail.cpp.

Referenced by SetupDetail(), and xTermsDetail().

RDetail

double RDetail

static

Definition at line 127 of file Detail.cpp.

Referenced by Alpha0Detail(), AlpharDetail(), DensityDetail(), PressureDetail(), PropertiesDetail(), and SetupDetail().

th0i

double th0i[MaxFlds+1][7+1]

static

Definition at line 138 of file Detail.cpp.

Referenced by Alpha0Detail(), Alpha0GERG(), SetupDetail(), and SetupGERG().

Told

double Told

static

Definition at line 137 of file Detail.cpp.

Referenced by AlpharDetail(), AlpharGERG(), ReducingParametersGERG(), SetupDetail(), and SetupGERG().

Tun

double Tun[NTerms+1]

static

Definition at line 137 of file Detail.cpp.

Referenced by AlpharDetail().

Uij5

double Uij5[MaxFlds+1][MaxFlds+1]

static

Definition at line 136 of file Detail.cpp.

Referenced by SetupDetail(), and xTermsDetail().

un

double un[NTerms+1]

static

Definition at line 131 of file Detail.cpp.

Referenced by AlpharDetail(), SetupDetail(), and xTermsDetail().

xold

double xold[MaxFlds+1]

static

Definition at line 139 of file Detail.cpp.

Referenced by ReducingParametersGERG(), SetupDetail(), SetupGERG(), and xTermsDetail().