github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-02-08 21:05:14 -05:00
Code Issues Packages Projects Releases Wiki Activity

Moved mixture derivatives into their own file and more functions take references rather than pointers

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-08-25 13:32:25 +02:00
parent 6d5b36d26f
commit 9435b2d136
7 changed files with 632 additions and 569 deletions
Download Patch File Download Diff File Expand all files Collapse all files

20
src/Backends/Helmholtz/FlashRoutines.cpp
View File

@@ -75,7 +75,7 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
options.omega = omega;
// Actually call the solver
SaturationSolvers::saturation_T_pure(&HEOS, HEOS._T, options);
SaturationSolvers::saturation_T_pure(HEOS, HEOS._T, options);
// If you get here, there was no error, all is well
break;
@@ -91,7 +91,7 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
catch(std::exception &){
try{
// We may need to polish the solution at low pressure
SaturationSolvers::saturation_T_pure_1D_P(&HEOS, T, options);
SaturationSolvers::saturation_T_pure_1D_P(HEOS, T, options);
}
catch(std::exception &){
throw;
@@ -144,10 +144,10 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
options.Nstep_max = 20;
// Get an extremely rough guess by interpolation of ln(p) v. T curve where the limits are mole-fraction-weighted
long double pguess = SaturationSolvers::saturation_preconditioner(&HEOS, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions);
long double pguess = SaturationSolvers::saturation_preconditioner(HEOS, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions);
// Use Wilson iteration to obtain updated guess for pressure
pguess = SaturationSolvers::saturation_Wilson(&HEOS, HEOS._Q, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions, pguess);
pguess = SaturationSolvers::saturation_Wilson(HEOS, HEOS._Q, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions, pguess);
// Actually call the successive substitution solver
SaturationSolvers::successive_substitution(HEOS, HEOS._Q, HEOS._T, pguess, HEOS.mole_fractions, HEOS.K, options);
@@ -210,7 +210,7 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
options.omega = omega;
// Actually call the solver
SaturationSolvers::saturation_PHSU_pure(&HEOS, HEOS._p, options);
SaturationSolvers::saturation_PHSU_pure(HEOS, HEOS._p, options);
// If you get here, there was no error, all is well
break;
@@ -225,7 +225,7 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
}
catch(std::exception &){
// We may need to polish the solution at low pressure
SaturationSolvers::saturation_P_pure_1D_T(&HEOS, HEOS._p, options);
SaturationSolvers::saturation_P_pure_1D_T(HEOS, HEOS._p, options);
}
// Load the outputs
@@ -250,10 +250,10 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
io.Nstep_max = 20;
// Get an extremely rough guess by interpolation of ln(p) v. T curve where the limits are mole-fraction-weighted
long double Tguess = SaturationSolvers::saturation_preconditioner(&HEOS, HEOS._p, SaturationSolvers::imposed_p, HEOS.mole_fractions);
long double Tguess = SaturationSolvers::saturation_preconditioner(HEOS, HEOS._p, SaturationSolvers::imposed_p, HEOS.mole_fractions);
// Use Wilson iteration to obtain updated guess for temperature
Tguess = SaturationSolvers::saturation_Wilson(&HEOS, HEOS._Q, HEOS._p, SaturationSolvers::imposed_p, HEOS.mole_fractions, Tguess);
Tguess = SaturationSolvers::saturation_Wilson(HEOS, HEOS._Q, HEOS._p, SaturationSolvers::imposed_p, HEOS.mole_fractions, Tguess);
// Actually call the successive substitution solver
SaturationSolvers::successive_substitution(HEOS, HEOS._Q, Tguess, HEOS._p, HEOS.mole_fractions, HEOS.K, io);
@@ -350,14 +350,14 @@ void FlashRoutines::PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
if (HEOS._rhomolar > HEOS._crit.rhomolar)
{
options.imposed_rho = SaturationSolvers::saturation_D_pure_options::IMPOSED_RHOL;
SaturationSolvers::saturation_D_pure(&HEOS, HEOS._rhomolar, options);
SaturationSolvers::saturation_D_pure(HEOS, HEOS._rhomolar, options);
// SatL and SatV have the saturation values
Sat = HEOS.SatL;
}
else
{
options.imposed_rho = SaturationSolvers::saturation_D_pure_options::IMPOSED_RHOV;
SaturationSolvers::saturation_D_pure(&HEOS, HEOS._rhomolar, options);
SaturationSolvers::saturation_D_pure(HEOS, HEOS._rhomolar, options);
// SatL and SatV have the saturation values
Sat = HEOS.SatV;
}

229
src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp
View File

@@ -26,6 +26,7 @@
#include "VLERoutines.h"
#include "FlashRoutines.h"
#include "TransportRoutines.h"
#include "MixtureDerivatives.h"
static int deriv_counter = 0;
@@ -1132,7 +1133,7 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot
// Run the saturation routines to determine the saturation densities and pressures
HelmholtzEOSMixtureBackend HEOS(components);
SaturationSolvers::saturation_T_pure_options options;
SaturationSolvers::saturation_T_pure(&HEOS, _T, options);
SaturationSolvers::saturation_T_pure(HEOS, _T, options);
long double Q;
@@ -1991,7 +1992,7 @@ long double HelmholtzEOSMixtureBackend::calc_gibbsmolar(void)
}
long double HelmholtzEOSMixtureBackend::calc_fugacity_coefficient(int i)
{
return exp(mixderiv_ln_fugacity_coefficient(i));
return exp(MixtureDerivatives::ln_fugacity_coefficient(*this, i));
}
SimpleState HelmholtzEOSMixtureBackend::calc_reducing_state_nocache(const std::vector<long double> & mole_fractions)
@@ -2291,228 +2292,4 @@ long double HelmholtzEOSMixtureBackend::calc_d3alpha0_dTau3(void)
}
long double HelmholtzEOSMixtureBackend::mixderiv_dalphar_dxi(int i)
{
return components[i]->pEOS->baser(_tau, _delta) + Excess.dalphar_dxi(_tau, _delta, mole_fractions, i);
}
long double HelmholtzEOSMixtureBackend::mixderiv_d2alphar_dxi_dTau(int i)
{
return components[i]->pEOS->dalphar_dTau(_tau, _delta) + Excess.d2alphar_dxi_dTau(_tau, _delta, mole_fractions, i);
}
long double HelmholtzEOSMixtureBackend::mixderiv_d2alphar_dxi_dDelta(int i)
{
return components[i]->pEOS->dalphar_dDelta(_tau, _delta) + Excess.d2alphar_dxi_dDelta(_tau, _delta, mole_fractions, i);
}
long double HelmholtzEOSMixtureBackend::mixderiv_d2alphardxidxj(int i, int j)
{
return 0 + Excess.d2alphardxidxj(_tau, _delta, mole_fractions, i, j);
}
long double HelmholtzEOSMixtureBackend::mixderiv_ln_fugacity_coefficient(int i)
{
return alphar() + mixderiv_ndalphar_dni__constT_V_nj(i)-log(1+_delta.pt()*dalphar_dDelta());
}
long double HelmholtzEOSMixtureBackend::mixderiv_dln_fugacity_coefficient_dT__constrho_n(int i)
{
double dtau_dT = -_tau.pt()/_T; //[1/K]
return (dalphar_dTau() + mixderiv_d_ndalphardni_dTau(i)-1/(1+_delta.pt()*dalphar_dDelta())*(_delta.pt()*d2alphar_dDelta_dTau()))*dtau_dT;
}
long double HelmholtzEOSMixtureBackend::mixderiv_dln_fugacity_coefficient_drho__constT_n(int i)
{
double ddelta_drho = 1/_reducing.rhomolar; //[m^3/mol]
return (dalphar_dDelta() + mixderiv_d_ndalphardni_dDelta(i)-1/(1+_delta.pt()*dalphar_dDelta())*(_delta.pt()*d2alphar_dDelta2()+dalphar_dDelta()))*ddelta_drho;
}
long double HelmholtzEOSMixtureBackend::mixderiv_dnalphar_dni__constT_V_nj(int i)
{
// GERG Equation 7.42
return alphar() + mixderiv_ndalphar_dni__constT_V_nj(i);
}
long double HelmholtzEOSMixtureBackend::mixderiv_d2nalphar_dni_dT(int i)
{
return -_tau.pt()/_T*(dalphar_dTau() + mixderiv_d_ndalphardni_dTau(i));
}
long double HelmholtzEOSMixtureBackend::mixderiv_dln_fugacity_coefficient_dT__constp_n(int i)
{
double T = _reducing.T/_tau.pt();
long double R_u = static_cast<long double>(_gas_constant);
return mixderiv_d2nalphar_dni_dT(i) + 1/T-mixderiv_partial_molar_volume(i)/(R_u*T)*mixderiv_dpdT__constV_n();
}
long double HelmholtzEOSMixtureBackend::mixderiv_partial_molar_volume(int i)
{
return -mixderiv_ndpdni__constT_V_nj(i)/mixderiv_ndpdV__constT_n();
}
long double HelmholtzEOSMixtureBackend::mixderiv_dln_fugacity_coefficient_dp__constT_n(int i)
{
// GERG equation 7.30
long double R_u = static_cast<long double>(_gas_constant);
double partial_molar_volume = mixderiv_partial_molar_volume(i); // [m^3/mol]
double term1 = partial_molar_volume/(R_u*_T); // m^3/mol/(N*m)*mol = m^2/N = 1/Pa
double term2 = 1.0/p();
return term1 - term2;
}
long double HelmholtzEOSMixtureBackend::mixderiv_dln_fugacity_coefficient_dxj__constT_p_xi(int i, int j)
{
// Gernert 3.115
long double R_u = static_cast<long double>(_gas_constant);
// partial molar volume is -dpdn/dpdV, so need to flip the sign here
return mixderiv_d2nalphar_dni_dxj__constT_V(i,j) - mixderiv_partial_molar_volume(i)/(R_u*_T)*mixderiv_dpdxj__constT_V_xi(j);
}
long double HelmholtzEOSMixtureBackend::mixderiv_dpdxj__constT_V_xi(int j)
{
// Gernert 3.130
long double R_u = static_cast<long double>(_gas_constant);
return _rhomolar*R_u*_T*(mixderiv_ddelta_dxj__constT_V_xi(j)*dalphar_dDelta()+_delta.pt()*mixderiv_d_dalpharddelta_dxj__constT_V_xi(j));
}
long double HelmholtzEOSMixtureBackend::mixderiv_d_dalpharddelta_dxj__constT_V_xi(int j)
{
// Gernert Equation 3.134 (Catch test provided)
return d2alphar_dDelta2()*mixderiv_ddelta_dxj__constT_V_xi(j)
+ d2alphar_dDelta_dTau()*mixderiv_dtau_dxj__constT_V_xi(j)
+ mixderiv_d2alphar_dxi_dDelta(j);
}
long double HelmholtzEOSMixtureBackend::mixderiv_dalphar_dxj__constT_V_xi(int j)
{
//Gernert 3.119 (Catch test provided)
return dalphar_dDelta()*mixderiv_ddelta_dxj__constT_V_xi(j)+dalphar_dTau()*mixderiv_dtau_dxj__constT_V_xi(j)+mixderiv_dalphar_dxi(j);
}
long double HelmholtzEOSMixtureBackend::mixderiv_d_ndalphardni_dxj__constT_V_xi(int i, int j)
{
// Gernert 3.118
return mixderiv_d_ndalphardni_dxj__constdelta_tau_xi(i,j)
+ mixderiv_ddelta_dxj__constT_V_xi(j)*mixderiv_d_ndalphardni_dDelta(i)
+ mixderiv_dtau_dxj__constT_V_xi(j)*mixderiv_d_ndalphardni_dTau(i);
}
long double HelmholtzEOSMixtureBackend::mixderiv_ddelta_dxj__constT_V_xi(int j)
{
// Gernert 3.121 (Catch test provided)
return -_delta.pt()/_reducing.rhomolar*Reducing.p->drhormolardxi__constxj(mole_fractions,j);
}
long double HelmholtzEOSMixtureBackend::mixderiv_dtau_dxj__constT_V_xi(int j)
{
// Gernert 3.122 (Catch test provided)
return 1/_T*Reducing.p->dTrdxi__constxj(mole_fractions,j);
}
long double HelmholtzEOSMixtureBackend::mixderiv_dpdT__constV_n()
{
long double R_u = static_cast<long double>(_gas_constant);
return _rhomolar*R_u*(1+_delta.pt()*dalphar_dDelta()-_delta.pt()*_tau.pt()*d2alphar_dDelta_dTau());
}
long double HelmholtzEOSMixtureBackend::mixderiv_dpdrho__constT_n()
{
long double R_u = static_cast<long double>(_gas_constant);
return R_u*_T*(1+2*_delta.pt()*dalphar_dDelta()+pow(_delta.pt(),2)*d2alphar_dDelta2());
}
long double HelmholtzEOSMixtureBackend::mixderiv_ndpdV__constT_n()
{
long double R_u = static_cast<long double>(_gas_constant);
return -pow(_rhomolar,2)*R_u*_T*(1+2*_delta.pt()*dalphar_dDelta()+pow(_delta.pt(),2)*d2alphar_dDelta2());
}
long double HelmholtzEOSMixtureBackend::mixderiv_ndpdni__constT_V_nj(int i)
{
// Eqn 7.64 and 7.63
long double R_u = static_cast<long double>(_gas_constant);
double ndrhorbar_dni__constnj = Reducing.p->ndrhorbardni__constnj(mole_fractions,i);
double ndTr_dni__constnj = Reducing.p->ndTrdni__constnj(mole_fractions,i);
double summer = 0;
for (unsigned int k = 0; k < mole_fractions.size(); ++k)
{
summer += mole_fractions[k]*mixderiv_d2alphar_dxi_dDelta(k);
}
double nd2alphar_dni_dDelta = _delta.pt()*d2alphar_dDelta2()*(1-1/_reducing.rhomolar*ndrhorbar_dni__constnj)+_tau.pt()*d2alphar_dDelta_dTau()/_reducing.T*ndTr_dni__constnj+mixderiv_d2alphar_dxi_dDelta(i)-summer;
return _rhomolar*R_u*_T*(1+_delta.pt()*dalphar_dDelta()*(2-1/_reducing.rhomolar*ndrhorbar_dni__constnj)+_delta.pt()*nd2alphar_dni_dDelta);
}
long double HelmholtzEOSMixtureBackend::mixderiv_ndalphar_dni__constT_V_nj(int i)
{
double term1 = _delta.pt()*dalphar_dDelta()*(1-1/_reducing.rhomolar*Reducing.p->ndrhorbardni__constnj(mole_fractions,i));
double term2 = _tau.pt()*dalphar_dTau()*(1/_reducing.T)*Reducing.p->ndTrdni__constnj(mole_fractions,i);
double s = 0;
for (unsigned int k = 0; k < mole_fractions.size(); k++)
{
s += mole_fractions[k]*mixderiv_dalphar_dxi(k);
}
double term3 = mixderiv_dalphar_dxi(i);
return term1 + term2 + term3 - s;
}
long double HelmholtzEOSMixtureBackend::mixderiv_ndln_fugacity_coefficient_dnj__constT_p(int i, int j)
{
long double R_u = static_cast<long double>(_gas_constant);
return mixderiv_nd2nalphardnidnj__constT_V(j, i) + 1 - mixderiv_partial_molar_volume(j)/(R_u*_T)*mixderiv_ndpdni__constT_V_nj(i);
}
long double HelmholtzEOSMixtureBackend::mixderiv_nddeltadni__constT_V_nj(int i)
{
return _delta.pt()-_delta.pt()/_reducing.rhomolar*Reducing.p->ndrhorbardni__constnj(mole_fractions, i);
}
long double HelmholtzEOSMixtureBackend::mixderiv_ndtaudni__constT_V_nj(int i)
{
return _tau.pt()/_reducing.T*Reducing.p->ndTrdni__constnj(mole_fractions, i);
}
long double HelmholtzEOSMixtureBackend::mixderiv_d_ndalphardni_dxj__constdelta_tau_xi(int i, int j)
{
double line1 = _delta.pt()*mixderiv_d2alphar_dxi_dDelta(j)*(1-1/_reducing.rhomolar*Reducing.p->ndrhorbardni__constnj(mole_fractions, i));
double line2 = -_delta.pt()*dalphar_dDelta()*(1/_reducing.rhomolar)*(Reducing.p->d_ndrhorbardni_dxj__constxi(mole_fractions, i, j)-1/_reducing.rhomolar*Reducing.p->drhormolardxi__constxj(mole_fractions,j)*Reducing.p->ndrhorbardni__constnj(mole_fractions,i));
double line3 = _tau.pt()*mixderiv_d2alphar_dxi_dTau(j)*(1/_reducing.T)*Reducing.p->ndTrdni__constnj(mole_fractions, i);
double line4 = _tau.pt()*dalphar_dTau()*(1/_reducing.T)*(Reducing.p->d_ndTrdni_dxj__constxi(mole_fractions,i,j)-1/_reducing.T*Reducing.p->dTrdxi__constxj(mole_fractions,j)*Reducing.p->ndTrdni__constnj(mole_fractions, i));
double s = 0;
for (unsigned int m = 0; m < mole_fractions.size(); m++)
{
s += mole_fractions[m]*mixderiv_d2alphardxidxj(j,m);
}
double line5 = mixderiv_d2alphardxidxj(i,j)-mixderiv_dalphar_dxi(j)-s;
return line1+line2+line3+line4+line5;
}
long double HelmholtzEOSMixtureBackend::mixderiv_nd2nalphardnidnj__constT_V(int i, int j)
{
double line0 = mixderiv_ndalphar_dni__constT_V_nj(j); // First term from 7.46
double line1 = mixderiv_d_ndalphardni_dDelta(i)*mixderiv_nddeltadni__constT_V_nj(j);
double line2 = mixderiv_d_ndalphardni_dTau(i)*mixderiv_ndtaudni__constT_V_nj(j);
double summer = 0;
for (unsigned int k = 0; k < mole_fractions.size(); k++)
{
summer += mole_fractions[k]*mixderiv_d_ndalphardni_dxj__constdelta_tau_xi(i, k);
}
double line3 = mixderiv_d_ndalphardni_dxj__constdelta_tau_xi(i, j)-summer;
return line0 + line1 + line2 + line3;
}
long double HelmholtzEOSMixtureBackend::mixderiv_d_ndalphardni_dDelta(int i)
{
// The first line
double term1 = (_delta.pt()*d2alphar_dDelta2()+dalphar_dDelta())*(1-1/_reducing.rhomolar*Reducing.p->ndrhorbardni__constnj(mole_fractions, i));
// The second line
double term2 = _tau.pt()*d2alphar_dDelta_dTau()*(1/_reducing.T)*Reducing.p->ndTrdni__constnj(mole_fractions, i);
// The third line
double term3 = mixderiv_d2alphar_dxi_dDelta(i);
for (unsigned int k = 0; k < mole_fractions.size(); k++)
{
term3 -= mole_fractions[k]*mixderiv_d2alphar_dxi_dDelta(k);
}
return term1 + term2 + term3;
}
long double HelmholtzEOSMixtureBackend::mixderiv_d_ndalphardni_dTau(int i)
{
// The first line
double term1 = _delta.pt()*d2alphar_dDelta_dTau()*(1-1/_reducing.rhomolar*Reducing.p->ndrhorbardni__constnj(mole_fractions, i));
// The second line
double term2 = (_tau.pt()*d2alphar_dTau2()+dalphar_dTau())*(1/_reducing.T)*Reducing.p->ndTrdni__constnj(mole_fractions, i);
// The third line
double term3 = mixderiv_d2alphar_dxi_dTau(i);
for (unsigned int k = 0; k < mole_fractions.size(); k++)
{
term3 -= mole_fractions[k]*mixderiv_d2alphar_dxi_dTau(k);
}
return term1 + term2 + term3;
}
} /* namespace CoolProp */

244
src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h
View File

@@ -41,9 +41,7 @@ public:
friend class FlashRoutines; // Allows the static methods in the FlashRoutines class to have access to all the protected members and methods of this class
friend class TransportRoutines; // Allows the static methods in the TransportRoutines class to have access to all the protected members and methods of this class
//friend class MixtureDerivatives; //
friend class MixtureDerivatives; // Allows the static methods in the MixtureDerivatives class to have access to all the protected members and methods of this class
// Helmholtz EOS backend uses mole fractions
bool using_mole_fractions(){return true;}
@@ -241,246 +239,6 @@ public:
long double solver_rho_Tp_SRK(long double T, long double p, int phase);
long double solver_for_rho_given_T_oneof_HSU(long double T, long double value, int other);
// ***************************************************************
// ***************************************************************
// ***************** MIXTURE DERIVATIVES ***********************
// ***************************************************************
// ***************************************************************
long double mixderiv_dalphar_dxi(int i);
long double mixderiv_d2alphar_dxi_dTau(int i);
long double mixderiv_d2alphar_dxi_dDelta(int i);
long double mixderiv_d2alphardxidxj(int i, int j);
/*! \brief GERG 2004 Monograph equation 7.61
*
* The derivative term
\f[
\left(\frac{\partial p}{\partial T} \right)_{V,\bar n} = \rho R(1+\delta \alpha_{\delta}^r-\delta \tau \alpha^r_{\delta\tau})
\f]
*/
long double mixderiv_dpdT__constV_n();
long double mixderiv_dpdrho__constT_n();
/** \brief GERG 2004 Monograph equation 7.62
*
* The derivative term
\f[
n\left(\frac{\partial p}{\partial V} \right)_{T,\bar n} = -\rho^2 RT(1+2\delta \alpha_{\delta}^r+\delta^2\alpha^r_{\delta\delta})
\f]
*/
long double mixderiv_ndpdV__constT_n();
/** \brief GERG 2004 Monograph equation 7.63
*
* The derivative term
\f[
n\left(\frac{\partial p}{\partial n_i} \right)_{T,V,n_j} = \rho RT\left[1+\delta\alpha_{\delta}^r\left[2- \frac{1}{\rho_r}\cdot n\left( \frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] +\delta\cdot n\left(\frac{\partial\alpha_{\delta}^r}{\partial n_i}\right)_{T,V,n_j}\right]
\f]
*/
long double mixderiv_ndpdni__constT_V_nj(int i);
/** GERG 2004 monograph Eqn. 7.32
* The partial molar volume
\f[
\hat v_i = \left( \frac{\partial V}{\partial n_i}\right)_{T,p,n_j} = \frac{-\left(\dfrac{\partial p}{\partial n_i}\right)_{T,V,n_j}}{\left(\dfrac{\partial p}{\partial V}\right)_{T,\bar n}}
\f]
*/
long double mixderiv_partial_molar_volume(int i);
/*!
Natural logarithm of the fugacity coefficient
*/
long double mixderiv_ln_fugacity_coefficient(int i);
/*!
Derivative of the natural logarithm of the fugacity coefficient with respect to T
*/
long double mixderiv_dln_fugacity_coefficient_dT__constrho_n(int i);
/*!
Derivative of the natural logarithm of the fugacity coefficient with respect to T
*/
long double mixderiv_dln_fugacity_coefficient_drho__constT_n(int i);
/// GERG 2004 Monograph Eqn. 7.29
/** The derivative term
\f[
\left(\frac{\partial \ln \phi_i}{\partial T} \right)_{p,\bar n} = \left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} + \frac{1}{T}-\frac{\hat v}{RT}\left(\frac{\partial p}{\partial T}\right)_{V,\bar n}
\f]
*/
long double mixderiv_dln_fugacity_coefficient_dT__constp_n(int i);
/// Table B4, Kunz, JCED, 2012 for the original term and the subsequent substitutions
/*! The derivative term
\f[
n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j}
\f]
which is equal to
\f{eqnarray*}{
n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} &=& \delta \phi^r_{\delta}\left[ 1-\frac{1}{\rho_r}\left[\left(\frac{\partial \rho_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial \rho_r}{\partial x_k}\right)_{x_j} \right]\right]\\
&& +\tau \phi^r_{\tau}\frac{1}{T_r}\left[\left(\frac{\partial T_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial T_r}{\partial x_k}\right)_{x_j} \right]\\
&& +\phi^r_{x_i}-\sum_{k=1}^{N}x_k\phi^r_{x_k}
\f}
*/
long double mixderiv_ndalphar_dni__constT_V_nj(int i);
/// GERG Equation 7.42
long double mixderiv_dnalphar_dni__constT_V_nj(int i);
/// GERG 2004 Monograph Eqn. 7.30
/**
The derivative term
\f[
\left(\frac{\partial \ln \phi_i}{\partial p} \right)_{T,\bar n} = \frac{\hat v_i}{RT}-\frac{1}{p}
\f]
*/
long double mixderiv_dln_fugacity_coefficient_dp__constT_n(int i);
///GERG 2004 Monograph Equation 7.31
/** The derivative term
\f[
n\left(\frac{\partial \ln \phi_i}{\partial n_j}\right)_{T,p} = n\left(\frac{\partial^2n\alpha^r}{\partial n_j \partial n_i} \right)_{T,V}+1+\frac{n}{RT}\frac{\left(\frac{\partial p}{\partial n_j}\right)_{T,V,n_i}\left(\frac{\partial p}{\partial n_i}\right)_{T,V,n_j}}{\left(\frac{\partial p}{\partial V}\right)_{T,\bar n}}
\f]
which is also equal to
\f[
n\left(\frac{\partial \ln \phi_i}{\partial n_j}\right)_{T,p} = n\left(\frac{\partial^2n\alpha^r}{\partial n_j \partial n_i} \right)_{T,V}+1-\frac{\hat v_i}{RT}\left[n\left(\frac{\partial p}{\partial n_j}\right)_{T,V,n_i}\right]
\f]
*/
long double mixderiv_ndln_fugacity_coefficient_dnj__constT_p(int i, int j);
/// Gernert Equation 3.115
/** The derivative term
\f[
\left(\frac{\partial \ln \phi_i}{\partial x_j}\right)_{T,p,x_{k\neq j}} = \left(\frac{\partial^2n\alpha^r}{\partial x_j \partial n_i} \right)_{T,V}+\frac{1}{RT}\frac{\left(\frac{\partial p}{\partial n_i}\right)_{T,V,n_{k\neq i}}\left(\frac{\partial p}{\partial x_j}\right)_{T,V,x_{k\neq j}}}{\left(\frac{\partial p}{\partial V}\right)_{T,\bar n}}
\f]
*/
long double mixderiv_dln_fugacity_coefficient_dxj__constT_p_xi(int i, int j);
/// Gernert Equation 3.130
/** The derivative term
\f[
\left(\frac{\partial p}{\partial x_j} \right)_{T,V,x_{k\neq j}} = \rho RT\left(-\frac{1}{\rho_r}\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_{k\neq j}} \delta\alpha_{\delta}^r + \delta\left(\frac{\partial}{\partial x_j}\left(\left( \frac{\partial \alpha^r}{\partial \delta}\right)_{\tau,\bar x}\right)\right)_{T,V,x_{k\neq j}}\right)
\f]
*/
long double mixderiv_dpdxj__constT_V_xi(int j);
/// Gernert Equation 3.117
/** The derivative term
\f[
\left(\frac{\partial^2n\alpha^r}{\partial x_i\partial n_j} \right)_{T,V} = \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_{j\neq i}}\right)\right)_{T,V,x_{k\neq j}} +\left(\frac{\partial \alpha^r}{\partial x_j}\right)_{T,V,x_{k\neq j}}
\f]
*/
long double mixderiv_d2nalphar_dxi_dxj__constT_V(int i, int j){ return mixderiv_d_ndalphardni_dxj__constT_V_xi(i, j) + mixderiv_dalphar_dxj__constT_V_xi(j);};
/// Gernert Equation 3.119
/// Catch test provided
long double mixderiv_dalphar_dxj__constT_V_xi(int j);
/// Gernert Equation 3.118
/// Catch test provided
long double mixderiv_d_ndalphardni_dxj__constT_V_xi(int i, int j);
/// Gernert Equation 3.134
/// Catch test provided
long double mixderiv_d_dalpharddelta_dxj__constT_V_xi(int j);
/// Gernert Equation 3.121
/// Catch test provided
long double mixderiv_ddelta_dxj__constT_V_xi(int j);
/// Gernert Equation 3.122
/// Catch test provided
long double mixderiv_dtau_dxj__constT_V_xi(int j);
/// GERG 2004 Monograph, equations 7.44 and 7.51
/** The derivative term
\f[
\left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} = \left( \frac{\partial}{\partial T}\left(\frac{\partial n \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right)_{V,\bar n}
\f]
\f[
\left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} = -\frac{\tau}{T}\left[\alpha_{\tau}^r +\left( \frac{\partial}{\partial \tau}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\bar x}\right]
\f]
*/
long double mixderiv_d2nalphar_dni_dT(int i);
/// GERG 2004 Monograph Equation 7.51 and Table B4, Kunz, JCED, 2012
/** The derivative term
\f{eqnarray*}{
\frac{\partial }{\partial \tau} \left( n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} \right) &=& \delta \phi^r_{\delta\tau}\left[ 1-\frac{1}{\rho_r}\left[\left(\frac{\partial \rho_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial \rho_r}{\partial x_k}\right)_{x_j} \right]\right]\\
&& +(\tau \phi^r_{\tau\tau}+\phi^r_{\tau})\frac{1}{T_r}\left[\left(\frac{\partial T_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial T_r}{\partial x_k}\right)_{x_j} \right]\\
&& +\phi^r_{x_i\tau}-\sum_{k=1}^{N}x_k\phi^r_{x_k\tau}
\f}
*/
long double mixderiv_d_ndalphardni_dTau(int i);
/** \brief GERG 2004 Monograph Equation 7.50 and Table B4, Kunz, JCED, 2012
*
The derivative term
\f{eqnarray*}{
\left(\frac{\partial }{\partial \delta} \left( n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} \right)\right)_{\tau,\bar x} &=& (\alpha_{\delta}^r+\delta\alpha_{\delta\delta}^r)\left[1-\frac{1}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j} \right] \\
&+&\tau\alpha^r_{\delta\tau}\frac{1}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\\
&+&\phi^r_{\delta x_i}-\sum_{k=1}^{N}x_k\phi^r_{\delta x_k}
\f}
*/
long double mixderiv_d_ndalphardni_dDelta(int i);
/** \brief GERG 2004 Monograph equation 7.41
* The derivative term
\f[
n\left(\frac{\partial^2n\alpha^r}{\partial n_i \partial n_j} \right)_{T,V} = n\left( \frac{\partial}{\partial n_j}\left(\frac{\partial n\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)_{T,V,n_i}
\f]
and
GERG 2004 Monograph equation 7.46:
\f[
n\left( \frac{\partial}{\partial n_j}\left(\frac{\partial n\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)_{T,V,n_i} = n\left( \frac{\partial \alpha^r}{\partial n_j}\right)_{T,V,n_i} + n\left( \frac{\partial}{\partial n_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right) \right)_{T,V,n_i}
\f]
GERG 2004 Monograph equation 7.47:
\f{eqnarray*}{
n\left( \frac{\partial}{\partial n_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right) \right)_{T,V,n_i} &=& \left( \frac{\partial}{\partial \delta}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\tau,\bar x}\cdot n\left(\frac{\partial\delta}{\partial n_j}\right)_{T,V,n_i}\\
&+& \left( \frac{\partial}{\partial \tau}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\tau,\bar x}\cdot n\left(\frac{\partial\tau}{\partial n_j}\right)_{T,V,n_i}\\
&+& \left( \frac{\partial}{\partial x_j}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i}-\sum_{k=1}^{N}x_k \left( \frac{\partial}{\partial x_k}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i}\\
\f}
*/
long double mixderiv_nd2nalphardnidnj__constT_V(int i, int j);
/** \brief GERG 2004 Monograph equation 7.48
The derivative term
\f[
n\left(\frac{\partial \delta}{\partial n_i} \right)_{T,V,n_j} = \delta - \frac{\delta}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i} \right)_{n_j}
\f]
*/
long double mixderiv_nddeltadni__constT_V_nj(int i);
/** \brief GERG 2004 Monograph equation 7.49
*
The derivative term
\f[
n\left(\frac{\partial \tau}{\partial n_i} \right)_{T,V,n_j} = \frac{\tau}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i} \right)_{n_j}
\f]
*/
long double mixderiv_ndtaudni__constT_V_nj(int i);
/** \brief GERG 2004 Monograph equation 7.52
*
The derivative term
\f{eqnarray*}{
\left( \frac{\partial}{\partial x_j}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i} &=& \delta\alpha_{\delta x_j}^{r}\left[ 1-\frac{1}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] \\
&-& \delta\alpha_{\delta}^{r}\frac{1}{\rho_r}\left[ \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right)\right)_{x_i}-\frac{1}{\rho_r}\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_i}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] \\
&+& \tau\alpha_{\tau x_j}^r\frac{1}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\\
&+& \tau\alpha_{\tau}^{r}\frac{1}{T_r}\left[ \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\right)\right)_{x_i}-\frac{1}{T_r}\left(\frac{\partial T_r}{\partial x_j}\right)_{x_i}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\right] \\
&+& \alpha_{x_ix_j}^r-\alpha_{x_j}^r-\sum_{m=1}^Nx_m\alpha_{x_jx_m}^r
\f}
*/
long double mixderiv_d_ndalphardni_dxj__constdelta_tau_xi(int i, int j);
};
} /* namespace CoolProp */

229
src/Backends/Helmholtz/MixtureDerivatives.cpp Normal file
View File

@@ -0,0 +1,229 @@
#include "MixtureDerivatives.h"
namespace CoolProp{
long double MixtureDerivatives::dalphar_dxi(HelmholtzEOSMixtureBackend &HEOS, int i)
{
return HEOS.components[i]->pEOS->baser(HEOS._tau, HEOS._delta) + HEOS.Excess.dalphar_dxi(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i);
}
long double MixtureDerivatives::d2alphar_dxi_dTau(HelmholtzEOSMixtureBackend &HEOS, int i)
{
return HEOS.components[i]->pEOS->dalphar_dTau(HEOS._tau, HEOS._delta) + HEOS.Excess.d2alphar_dxi_dTau(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i);
}
long double MixtureDerivatives::d2alphar_dxi_dDelta(HelmholtzEOSMixtureBackend &HEOS, int i)
{
return HEOS.components[i]->pEOS->dalphar_dDelta(HEOS._tau, HEOS._delta) + HEOS.Excess.d2alphar_dxi_dDelta(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i);
}
long double MixtureDerivatives::d2alphardxidxj(HelmholtzEOSMixtureBackend &HEOS, int i, int j)
{
return 0 + HEOS.Excess.d2alphardxidxj(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i, j);
}
long double MixtureDerivatives::ln_fugacity_coefficient(HelmholtzEOSMixtureBackend &HEOS, int i)
{
return HEOS.alphar() + ndalphar_dni__constT_V_nj(HEOS, i)-log(1+HEOS._delta.pt()*HEOS.dalphar_dDelta());
}
long double MixtureDerivatives::dln_fugacity_coefficient_dT__constrho_n(HelmholtzEOSMixtureBackend &HEOS, int i)
{
double dtau_dT = -HEOS._tau.pt()/HEOS._T; //[1/K]
return (HEOS.dalphar_dTau() + d_ndalphardni_dTau(HEOS, i)-1/(1+HEOS._delta.pt()*HEOS.dalphar_dDelta())*(HEOS._delta.pt()*HEOS.d2alphar_dDelta_dTau()))*dtau_dT;
}
long double MixtureDerivatives::dln_fugacity_coefficient_drho__constT_n(HelmholtzEOSMixtureBackend &HEOS, int i)
{
double ddelta_drho = 1/HEOS._reducing.rhomolar; //[m^3/mol]
return (HEOS.dalphar_dDelta() + d_ndalphardni_dDelta(HEOS, i)-1/(1+HEOS._delta.pt()*HEOS.dalphar_dDelta())*(HEOS._delta.pt()*HEOS.d2alphar_dDelta2()+HEOS.dalphar_dDelta()))*ddelta_drho;
}
long double MixtureDerivatives::dnalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i)
{
// GERG Equation 7.42
return HEOS.alphar() + ndalphar_dni__constT_V_nj(HEOS, i);
}
long double MixtureDerivatives::d2nalphar_dni_dT(HelmholtzEOSMixtureBackend &HEOS, int i)
{
return -HEOS._tau.pt()/HEOS._T*(HEOS.dalphar_dTau() + d_ndalphardni_dTau(HEOS, i));
}
long double MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(HelmholtzEOSMixtureBackend &HEOS, int i)
{
double T = HEOS._reducing.T/HEOS._tau.pt();
long double R_u = HEOS.gas_constant();
return d2nalphar_dni_dT(HEOS, i) + 1/T-partial_molar_volume(HEOS, i)/(R_u*T)*dpdT__constV_n(HEOS);
}
long double MixtureDerivatives::partial_molar_volume(HelmholtzEOSMixtureBackend &HEOS, int i)
{
return -ndpdni__constT_V_nj(HEOS, i)/ndpdV__constT_n(HEOS);
}
long double MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(HelmholtzEOSMixtureBackend &HEOS, int i)
{
// GERG equation 7.30
long double R_u = HEOS.gas_constant();
double partial_molar_volumeval = partial_molar_volume(HEOS, i); // [m^3/mol]
double term1 = partial_molar_volumeval/(R_u*HEOS._T); // m^3/mol/(N*m)*mol = m^2/N = 1/Pa
double term2 = 1.0/HEOS.p();
return term1 - term2;
}
long double MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(HelmholtzEOSMixtureBackend &HEOS, int i, int j)
{
// Gernert 3.115
long double R_u = HEOS.gas_constant();
// partial molar volume is -dpdn/dpdV, so need to flip the sign here
return d2nalphar_dxi_dnj__constT_V(HEOS, i, j) - partial_molar_volume(HEOS, i)/(R_u*HEOS._T)*dpdxj__constT_V_xi(HEOS, j);
}
long double MixtureDerivatives::dpdxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j)
{
// Gernert 3.130
long double R_u = HEOS.gas_constant();
return HEOS._rhomolar*R_u*HEOS._T*(ddelta_dxj__constT_V_xi(HEOS, j)*HEOS.dalphar_dDelta()+HEOS._delta.pt()*d_dalpharddelta_dxj__constT_V_xi(HEOS, j));
}
long double MixtureDerivatives::d_dalpharddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j)
{
// Gernert Equation 3.134 (Catch test provided)
return HEOS.d2alphar_dDelta2()*ddelta_dxj__constT_V_xi(HEOS, j)
+ HEOS.d2alphar_dDelta_dTau()*dtau_dxj__constT_V_xi(HEOS, j)
+ d2alphar_dxi_dDelta(HEOS, j);
}
long double MixtureDerivatives::dalphar_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j)
{
//Gernert 3.119 (Catch test provided)
return HEOS.dalphar_dDelta()*ddelta_dxj__constT_V_xi(HEOS, j)+HEOS.dalphar_dTau()*dtau_dxj__constT_V_xi(HEOS, j)+dalphar_dxi(HEOS, j);
}
long double MixtureDerivatives::d_ndalphardni_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int i, int j)
{
// Gernert 3.118
return d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i,j)
+ ddelta_dxj__constT_V_xi(HEOS, j)*d_ndalphardni_dDelta(HEOS, i)
+ dtau_dxj__constT_V_xi(HEOS, j)*d_ndalphardni_dTau(HEOS, i);
}
long double MixtureDerivatives::ddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j)
{
// Gernert 3.121 (Catch test provided)
return -HEOS._delta.pt()/HEOS._reducing.rhomolar*HEOS.Reducing.p->drhormolardxi__constxj(HEOS.mole_fractions,j);
}
long double MixtureDerivatives::dtau_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j)
{
// Gernert 3.122 (Catch test provided)
return 1/HEOS._T*HEOS.Reducing.p->dTrdxi__constxj(HEOS.mole_fractions,j);
}
long double MixtureDerivatives::dpdT__constV_n(HelmholtzEOSMixtureBackend &HEOS)
{
long double R_u = HEOS.gas_constant();
return HEOS._rhomolar*R_u*(1+HEOS._delta.pt()*HEOS.dalphar_dDelta()-HEOS._delta.pt()*HEOS._tau.pt()*HEOS.d2alphar_dDelta_dTau());
}
long double MixtureDerivatives::dpdrho__constT_n(HelmholtzEOSMixtureBackend &HEOS)
{
long double R_u = HEOS.gas_constant();
return R_u*HEOS._T*(1+2*HEOS._delta.pt()*HEOS.dalphar_dDelta()+pow(HEOS._delta.pt(),2)*HEOS.d2alphar_dDelta2());
}
long double MixtureDerivatives::ndpdV__constT_n(HelmholtzEOSMixtureBackend &HEOS)
{
long double R_u = HEOS.gas_constant();
return -pow(HEOS._rhomolar,2)*R_u*HEOS._T*(1+2*HEOS._delta.pt()*HEOS.dalphar_dDelta()+pow(HEOS._delta.pt(),2)*HEOS.d2alphar_dDelta2());
}
long double MixtureDerivatives::ndpdni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i)
{
// Eqn 7.64 and 7.63
long double R_u = HEOS.gas_constant();
double ndrhorbar_dni__constnj = HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions,i);
double ndTr_dni__constnj = HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions,i);
double summer = 0;
for (unsigned int k = 0; k < HEOS.mole_fractions.size(); ++k)
{
summer += HEOS.mole_fractions[k]*d2alphar_dxi_dDelta(HEOS, k);
}
double nd2alphar_dni_dDelta = HEOS._delta.pt()*HEOS.d2alphar_dDelta2()*(1-1/HEOS._reducing.rhomolar*ndrhorbar_dni__constnj)+HEOS._tau.pt()*HEOS.d2alphar_dDelta_dTau()/HEOS._reducing.T*ndTr_dni__constnj+d2alphar_dxi_dDelta(HEOS, i)-summer;
return HEOS._rhomolar*R_u*HEOS._T*(1+HEOS._delta.pt()*HEOS.dalphar_dDelta()*(2-1/HEOS._reducing.rhomolar*ndrhorbar_dni__constnj)+HEOS._delta.pt()*nd2alphar_dni_dDelta);
}
long double MixtureDerivatives::ndalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i)
{
double term1 = HEOS._delta.pt()*HEOS.dalphar_dDelta()*(1-1/HEOS._reducing.rhomolar*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions,i));
double term2 = HEOS._tau.pt()*HEOS.dalphar_dTau()*(1/HEOS._reducing.T)*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions,i);
double s = 0;
for (unsigned int k = 0; k < HEOS.mole_fractions.size(); k++)
{
s += HEOS.mole_fractions[k]*dalphar_dxi(HEOS, k);
}
double term3 = dalphar_dxi(HEOS, i);
return term1 + term2 + term3 - s;
}
long double MixtureDerivatives::ndln_fugacity_coefficient_dnj__constT_p(HelmholtzEOSMixtureBackend &HEOS, int i, int j)
{
long double R_u = HEOS.gas_constant();
return nd2nalphardnidnj__constT_V(HEOS, j, i) + 1 - partial_molar_volume(HEOS, j)/(R_u*HEOS._T)*ndpdni__constT_V_nj(HEOS, i);
}
long double MixtureDerivatives::nddeltadni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i)
{
return HEOS._delta.pt()-HEOS._delta.pt()/HEOS._reducing.rhomolar*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions, i);
}
long double MixtureDerivatives::ndtaudni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i)
{
return HEOS._tau.pt()/HEOS._reducing.T*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i);
}
long double MixtureDerivatives::d_ndalphardni_dxj__constdelta_tau_xi(HelmholtzEOSMixtureBackend &HEOS, int i, int j)
{
double line1 = HEOS._delta.pt()*d2alphar_dxi_dDelta(HEOS, j)*(1-1/HEOS._reducing.rhomolar*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions, i));
double line2 = -HEOS._delta.pt()*HEOS.dalphar_dDelta()*(1/HEOS._reducing.rhomolar)*(HEOS.Reducing.p->d_ndrhorbardni_dxj__constxi(HEOS.mole_fractions, i, j)-1/HEOS._reducing.rhomolar*HEOS.Reducing.p->drhormolardxi__constxj(HEOS.mole_fractions,j)*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions,i));
double line3 = HEOS._tau.pt()*d2alphar_dxi_dTau(HEOS, j)*(1/HEOS._reducing.T)*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i);
double line4 = HEOS._tau.pt()*HEOS.dalphar_dTau()*(1/HEOS._reducing.T)*(HEOS.Reducing.p->d_ndTrdni_dxj__constxi(HEOS.mole_fractions,i,j)-1/HEOS._reducing.T*HEOS.Reducing.p->dTrdxi__constxj(HEOS.mole_fractions, j)*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i));
double s = 0;
for (unsigned int m = 0; m < HEOS.mole_fractions.size(); m++)
{
s += HEOS.mole_fractions[m]*d2alphardxidxj(HEOS, j,m);
}
double line5 = d2alphardxidxj(HEOS, i,j)-dalphar_dxi(HEOS, j)-s;
return line1+line2+line3+line4+line5;
}
long double MixtureDerivatives::nd2nalphardnidnj__constT_V(HelmholtzEOSMixtureBackend &HEOS, int i, int j)
{
double line0 = ndalphar_dni__constT_V_nj(HEOS, j); // First term from 7.46
double line1 = d_ndalphardni_dDelta(HEOS, i)*nddeltadni__constT_V_nj(HEOS, j);
double line2 = d_ndalphardni_dTau(HEOS, i)*ndtaudni__constT_V_nj(HEOS, j);
double summer = 0;
for (unsigned int k = 0; k < HEOS.mole_fractions.size(); k++)
{
summer += HEOS.mole_fractions[k]*d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, k);
}
double line3 = d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, j)-summer;
return line0 + line1 + line2 + line3;
}
long double MixtureDerivatives::d_ndalphardni_dDelta(HelmholtzEOSMixtureBackend &HEOS, int i)
{
// The first line
double term1 = (HEOS._delta.pt()*HEOS.d2alphar_dDelta2()+HEOS.dalphar_dDelta())*(1-1/HEOS._reducing.rhomolar*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions, i));
// The second line
double term2 = HEOS._tau.pt()*HEOS.d2alphar_dDelta_dTau()*(1/HEOS._reducing.T)*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i);
// The third line
double term3 = d2alphar_dxi_dDelta(HEOS, i);
for (unsigned int k = 0; k < HEOS.mole_fractions.size(); k++)
{
term3 -= HEOS.mole_fractions[k]*d2alphar_dxi_dDelta(HEOS, k);
}
return term1 + term2 + term3;
}
long double MixtureDerivatives::d_ndalphardni_dTau(HelmholtzEOSMixtureBackend &HEOS, int i)
{
// The first line
double term1 = HEOS._delta.pt()*HEOS.d2alphar_dDelta_dTau()*(1-1/HEOS._reducing.rhomolar*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions, i));
// The second line
double term2 = (HEOS._tau.pt()*HEOS.d2alphar_dTau2()+HEOS.dalphar_dTau())*(1/HEOS._reducing.T)*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i);
// The third line
double term3 = d2alphar_dxi_dTau(HEOS, i);
for (unsigned int k = 0; k < HEOS.mole_fractions.size(); k++)
{
term3 -= HEOS.mole_fractions[k]*d2alphar_dxi_dTau(HEOS, k);
}
return term1 + term2 + term3;
}
} /* namespace CoolProp */

296
src/Backends/Helmholtz/MixtureDerivatives.h Normal file
View File

@@ -0,0 +1,296 @@
/**
* This file contains derivatives needed in the mixture model. The derivatives are quite nasty, and there
* are a lot of them, so they are put in this file for cleanness. The MixtureDerivatives class is a friend class
* of the HelmholtzEOSMixtureBackend, so it can access all the private members of the HelmholtzEOSMixtureBackend
* class
*/
// ***************************************************************
// ***************************************************************
// ***************** MIXTURE DERIVATIVES ***********************
// ***************************************************************
// ***************************************************************
#ifndef MIXTURE_DERIVATIVES_H
#define MIXTURE_DERIVATIVES_H
#include "HelmholtzEOSMixtureBackend.h"
namespace CoolProp{
/**
This class is a friend class of HelmholtzEOSMixtureBackend, therefore the
static methods contained in it have access to the private and
protected variables in the HelmholtzEOSMixtureBackend instance.
In this way the derivative routines can be kept in their own separate file
and not pollute the HelmholtzEOSMixtureBackend namespace
*/
class MixtureDerivatives{
public:
static long double dalphar_dxi(HelmholtzEOSMixtureBackend &HEOS, int i);
static long double d2alphar_dxi_dTau(HelmholtzEOSMixtureBackend &HEOS, int i);
static long double d2alphar_dxi_dDelta(HelmholtzEOSMixtureBackend &HEOS, int i);
static long double d2alphardxidxj(HelmholtzEOSMixtureBackend &HEOS, int i, int j);
/** \brief GERG 2004 Monograph equation 7.61
*
* The derivative term
* \f[
* \left(\frac{\partial p}{\partial T} \right)_{V,\bar n} = \rho R(1+\delta \alpha_{\delta}^r-\delta \tau \alpha^r_{\delta\tau})
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double dpdT__constV_n(HelmholtzEOSMixtureBackend &HEOS);
static long double dpdrho__constT_n(HelmholtzEOSMixtureBackend &HEOS);
/** \brief GERG 2004 Monograph equation 7.62
*
* The derivative term
* \f[
* n\left(\frac{\partial p}{\partial V} \right)_{T,\bar n} = -\rho^2 RT(1+2\delta \alpha_{\delta}^r+\delta^2\alpha^r_{\delta\delta})
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double ndpdV__constT_n(HelmholtzEOSMixtureBackend &HEOS);
/** \brief GERG 2004 Monograph equation 7.63
*
* The derivative term
* \f[
* n\left(\frac{\partial p}{\partial n_i} \right)_{T,V,n_j} = \rho RT\left[1+\delta\alpha_{\delta}^r\left[2- \frac{1}{\rho_r}\cdot n\left( \frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] +\delta\cdot n\left(\frac{\partial\alpha_{\delta}^r}{\partial n_i}\right)_{T,V,n_j}\right]
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double ndpdni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i);
/** \brief GERG 2004 monograph Eqn. 7.32
*
* The partial molar volume
* \f[
* \hat v_i = \left( \frac{\partial V}{\partial n_i}\right)_{T,p,n_j} = \frac{-\left(\dfrac{\partial p}{\partial n_i}\right)_{T,V,n_j}}{\left(\dfrac{\partial p}{\partial V}\right)_{T,\bar n}}
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double partial_molar_volume(HelmholtzEOSMixtureBackend &HEOS, int i);
/** \brief Natural logarithm of the fugacity coefficient
*
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double ln_fugacity_coefficient(HelmholtzEOSMixtureBackend &HEOS, int i);
/** \brief Derivative of the natural logarithm of the fugacity coefficient with respect to T
*
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double dln_fugacity_coefficient_dT__constrho_n(HelmholtzEOSMixtureBackend &HEOS, int i);
/** \brief Derivative of the natural logarithm of the fugacity coefficient with respect to T
*
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double dln_fugacity_coefficient_drho__constT_n(HelmholtzEOSMixtureBackend &HEOS, int i);
/** \brief GERG 2004 Monograph Eqn. 7.29
*
* The derivative term
*
* \f[
* \left(\frac{\partial \ln \phi_i}{\partial T} \right)_{p,\bar n} = \left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} + \frac{1}{T}-\frac{\hat v}{RT}\left(\frac{\partial p}{\partial T}\right)_{V,\bar n}
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double dln_fugacity_coefficient_dT__constp_n(HelmholtzEOSMixtureBackend &HEOS, int i);
/** \brief Table B4, Kunz, JCED, 2012 for the original term and the subsequent substitutions
*
* The derivative term
* \f[
* n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j}
* \f]
* which is equal to
* \f{eqnarray*}{
* n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} &=& \delta \phi^r_{\delta}\left[ 1-\frac{1}{\rho_r}\left[\left(\frac{\partial \rho_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial \rho_r}{\partial x_k}\right)_{x_j} \right]\right]\\
* && +\tau \phi^r_{\tau}\frac{1}{T_r}\left[\left(\frac{\partial T_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial T_r}{\partial x_k}\right)_{x_j} \right]\\
* && +\phi^r_{x_i}-\sum_{k=1}^{N}x_k\phi^r_{x_k}
* \f}
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double ndalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i);
/// GERG Equation 7.42
static long double dnalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i);
/** \brief GERG 2004 Monograph Eqn. 7.30
*
* The derivative term
* \f[
* \left(\frac{\partial \ln \phi_i}{\partial p} \right)_{T,\bar n} = \frac{\hat v_i}{RT}-\frac{1}{p}
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double dln_fugacity_coefficient_dp__constT_n(HelmholtzEOSMixtureBackend &HEOS, int i);
/** \brief GERG 2004 Monograph Equation 7.31
*
* The derivative term
* \f[
* n\left(\frac{\partial \ln \phi_i}{\partial n_j}\right)_{T,p} = n\left(\frac{\partial^2n\alpha^r}{\partial n_j \partial n_i} \right)_{T,V}+1+\frac{n}{RT}\frac{\left(\frac{\partial p}{\partial n_j}\right)_{T,V,n_i}\left(\frac{\partial p}{\partial n_i}\right)_{T,V,n_j}}{\left(\frac{\partial p}{\partial V}\right)_{T,\bar n}}
* \f]
* which is also equal to
* \f[
* n\left(\frac{\partial \ln \phi_i}{\partial n_j}\right)_{T,p} = n\left(\frac{\partial^2n\alpha^r}{\partial n_j \partial n_i} \right)_{T,V}+1-\frac{\hat v_i}{RT}\left[n\left(\frac{\partial p}{\partial n_j}\right)_{T,V,n_i}\right]
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double ndln_fugacity_coefficient_dnj__constT_p(HelmholtzEOSMixtureBackend &HEOS, int i, int j);
/** \brief Gernert Equation 3.115
*
* The derivative term
* \f[
* \left(\frac{\partial \ln \phi_i}{\partial x_j}\right)_{T,p,x_{k\neq j}} = \left(\frac{\partial^2n\alpha^r}{\partial x_j \partial n_i} \right)_{T,V}+\frac{1}{RT}\frac{\left(\frac{\partial p}{\partial n_i}\right)_{T,V,n_{k\neq i}}\left(\frac{\partial p}{\partial x_j}\right)_{T,V,x_{k\neq j}}}{\left(\frac{\partial p}{\partial V}\right)_{T,\bar n}}
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double dln_fugacity_coefficient_dxj__constT_p_xi(HelmholtzEOSMixtureBackend &HEOS, int i, int j);
/** \brief Gernert Equation 3.130
*
* The derivative term
* \f[
* \left(\frac{\partial p}{\partial x_j} \right)_{T,V,x_{k\neq j}} = \rho RT\left(-\frac{1}{\rho_r}\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_{k\neq j}} \delta\alpha_{\delta}^r + \delta\left(\frac{\partial}{\partial x_j}\left(\left( \frac{\partial \alpha^r}{\partial \delta}\right)_{\tau,\bar x}\right)\right)_{T,V,x_{k\neq j}}\right)
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double dpdxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j);
/** \brief Gernert Equation 3.117
*
* The derivative term
* \f[
* \left(\frac{\partial^2n\alpha^r}{\partial x_i\partial n_j} \right)_{T,V} = \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_{j\neq i}}\right)\right)_{T,V,x_{k\neq j}} +\left(\frac{\partial \alpha^r}{\partial x_j}\right)_{T,V,x_{k\neq j}}
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double d2nalphar_dxi_dnj__constT_V(HelmholtzEOSMixtureBackend &HEOS, int i, int j){ return MixtureDerivatives::d_ndalphardni_dxj__constT_V_xi(HEOS, i, j) + MixtureDerivatives::dalphar_dxj__constT_V_xi(HEOS, j);};
/// Gernert Equation 3.119
/// Catch test provided
static long double dalphar_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j);
/// Gernert Equation 3.118
/// Catch test provided
static long double d_ndalphardni_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int i, int j);
/// Gernert Equation 3.134
/// Catch test provided
static long double d_dalpharddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j);
/// Gernert Equation 3.121
/// Catch test provided
static long double ddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j);
/// Gernert Equation 3.122
/// Catch test provided
static long double dtau_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j);
/** \brief GERG 2004 Monograph, equations 7.44 and 7.51
*
* The derivative term
* \f[
* \left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} = \left( \frac{\partial}{\partial T}\left(\frac{\partial n \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right)_{V,\bar n}
* \f]
* \f[
* \left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} = -\frac{\tau}{T}\left[\alpha_{\tau}^r +\left( \frac{\partial}{\partial \tau}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\bar x}\right]
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double d2nalphar_dni_dT(HelmholtzEOSMixtureBackend &HEOS, int i);
/** \brief GERG 2004 Monograph Equation 7.51 and Table B4, Kunz, JCED, 2012
*
* The derivative term
* \f{eqnarray*}{
* \frac{\partial }{\partial \tau} \left( n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} \right) &=& \delta \phi^r_{\delta\tau}\left[ 1-\frac{1}{\rho_r}\left[\left(\frac{\partial \rho_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial \rho_r}{\partial x_k}\right)_{x_j} \right]\right]\\
* && +(\tau \phi^r_{\tau\tau}+\phi^r_{\tau})\frac{1}{T_r}\left[\left(\frac{\partial T_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial T_r}{\partial x_k}\right)_{x_j} \right]\\
* && +\phi^r_{x_i\tau}-\sum_{k=1}^{N}x_k\phi^r_{x_k\tau}
* \f}
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double d_ndalphardni_dTau(HelmholtzEOSMixtureBackend &HEOS, int i);
/** \brief GERG 2004 Monograph Equation 7.50 and Table B4, Kunz, JCED, 2012
*
* The derivative term
* \f{eqnarray*}{
* \left(\frac{\partial }{\partial \delta} \left( n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} \right)\right)_{\tau,\bar x} &=& (\alpha_{\delta}^r+\delta\alpha_{\delta\delta}^r)\left[1-\frac{1}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j} \right] \\
* &+&\tau\alpha^r_{\delta\tau}\frac{1}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\\
* &+&\phi^r_{\delta x_i}-\sum_{k=1}^{N}x_k\phi^r_{\delta x_k}
* \f}
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double d_ndalphardni_dDelta(HelmholtzEOSMixtureBackend &HEOS, int i);
/** \brief GERG 2004 Monograph equation 7.41
*
* The derivative term
* \f[
* n\left(\frac{\partial^2n\alpha^r}{\partial n_i \partial n_j} \right)_{T,V} = n\left( \frac{\partial}{\partial n_j}\left(\frac{\partial n\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)_{T,V,n_i}
* \f]
* and
* GERG 2004 Monograph equation 7.46:
* \f[
* n\left( \frac{\partial}{\partial n_j}\left(\frac{\partial n\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)_{T,V,n_i} = n\left( \frac{\partial \alpha^r}{\partial n_j}\right)_{T,V,n_i} + n\left( \frac{\partial}{\partial n_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right) \right)_{T,V,n_i}
* \f]
* GERG 2004 Monograph equation 7.47:
* \f{eqnarray*}{
* n\left( \frac{\partial}{\partial n_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right) \right)_{T,V,n_i} &=& \left( \frac{\partial}{\partial \delta}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\tau,\bar x}\cdot n\left(\frac{\partial\delta}{\partial n_j}\right)_{T,V,n_i}\\
* &+& \left( \frac{\partial}{\partial \tau}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\tau,\bar x}\cdot n\left(\frac{\partial\tau}{\partial n_j}\right)_{T,V,n_i}\\
* &+& \left( \frac{\partial}{\partial x_j}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i}-\sum_{k=1}^{N}x_k \left( \frac{\partial}{\partial x_k}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i}\\
* \f}
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double nd2nalphardnidnj__constT_V(HelmholtzEOSMixtureBackend &HEOS, int i, int j);
/** \brief GERG 2004 Monograph equation 7.48
*
* The derivative term
* \f[
* n\left(\frac{\partial \delta}{\partial n_i} \right)_{T,V,n_j} = \delta - \frac{\delta}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i} \right)_{n_j}
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double nddeltadni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i);
/** \brief GERG 2004 Monograph equation 7.49
*
* The derivative term
* \f[
* n\left(\frac{\partial \tau}{\partial n_i} \right)_{T,V,n_j} = \frac{\tau}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i} \right)_{n_j}
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double ndtaudni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i);
/** \brief GERG 2004 Monograph equation 7.52
*
* The derivative term
* \f{eqnarray*}{
* \left( \frac{\partial}{\partial x_j}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i} &=& \delta\alpha_{\delta x_j}^{r}\left[ 1-\frac{1}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] \\
* &-& \delta\alpha_{\delta}^{r}\frac{1}{\rho_r}\left[ \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right)\right)_{x_i}-\frac{1}{\rho_r}\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_i}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] \\
* &+& \tau\alpha_{\tau x_j}^r\frac{1}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\\
* &+& \tau\alpha_{\tau}^{r}\frac{1}{T_r}\left[ \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\right)\right)_{x_i}-\frac{1}{T_r}\left(\frac{\partial T_r}{\partial x_j}\right)_{x_i}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\right] \\
* &+& \alpha_{x_ix_j}^r-\alpha_{x_j}^r-\sum_{m=1}^Nx_m\alpha_{x_jx_m}^r
* \f}
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double d_ndalphardni_dxj__constdelta_tau_xi(HelmholtzEOSMixtureBackend &HEOS, int i, int j);
}; /* class MixtureDerivatives */
} /* namepsace CoolProp*/
#endif

143
src/Backends/Helmholtz/VLERoutines.cpp
View File

@@ -2,10 +2,11 @@
#include "HelmholtzEOSMixtureBackend.h"
#include "VLERoutines.h"
#include "MatrixMath.h"
#include "MixtureDerivatives.h"
namespace CoolProp {
void SaturationSolvers::saturation_critical(HelmholtzEOSMixtureBackend *HEOS, parameters ykey, long double y){
void SaturationSolvers::saturation_critical(HelmholtzEOSMixtureBackend &HEOS, parameters ykey, long double y){
class inner_resid : public FuncWrapper1D{
public:
@@ -35,9 +36,9 @@ void SaturationSolvers::saturation_critical(HelmholtzEOSMixtureBackend *HEOS, pa
long double r, T, rhomolar_liq, rhomolar_vap, value, p, gL, gV, rhomolar_crit;
int other;
outer_resid(HelmholtzEOSMixtureBackend *HEOS, CoolProp::parameters ykey, long double y)
: HEOS(HEOS), ykey(ykey), y(y){
rhomolar_crit = HEOS->rhomolar_critical();
outer_resid(HelmholtzEOSMixtureBackend &HEOS, CoolProp::parameters ykey, long double y)
: HEOS(&HEOS), ykey(ykey), y(y){
rhomolar_crit = HEOS.rhomolar_critical();
};
double call(double rhomolar_vap){
this->y = y;
@@ -67,13 +68,13 @@ void SaturationSolvers::saturation_critical(HelmholtzEOSMixtureBackend *HEOS, pa
};
outer_resid resid(HEOS, iT, y);
double rhomolar_crit = HEOS->rhomolar_critical();
double rhomolar_crit = HEOS.rhomolar_critical();
std::string errstr;
Brent(&resid, rhomolar_crit*(1-1e-8), rhomolar_crit*0.5, DBL_EPSILON, 1e-9, 20, errstr);
}
void SaturationSolvers::saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_options &options)
void SaturationSolvers::saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend &HEOS, long double T, saturation_T_pure_options &options)
{
// Define the residual to be driven to zero
@@ -85,8 +86,8 @@ void SaturationSolvers::saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend *HEOS,
long double r, T, rhomolar_liq, rhomolar_vap, value, p, gL, gV;
int other;
solver_resid(HelmholtzEOSMixtureBackend *HEOS, long double T, long double rhomolar_liq_guess, long double rhomolar_vap_guess)
: HEOS(HEOS), T(T), rhomolar_liq(rhomolar_liq_guess), rhomolar_vap(rhomolar_vap_guess){};
solver_resid(HelmholtzEOSMixtureBackend &HEOS, long double T, long double rhomolar_liq_guess, long double rhomolar_vap_guess)
: HEOS(&HEOS), T(T), rhomolar_liq(rhomolar_liq_guess), rhomolar_vap(rhomolar_vap_guess){};
double call(double p){
this->p = p;
// Recalculate the densities using the current guess values
@@ -118,13 +119,13 @@ void SaturationSolvers::saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend *HEOS,
Secant(resid, options.p, options.p*1.1, 1e-10, 100, errstr);
}
catch(std::exception &){
long double pmax = std::min(options.p*1.03, static_cast<long double>(HEOS->p_critical()+1e-6));
long double pmin = std::max(options.p*0.97, static_cast<long double>(HEOS->p_triple()-1e-6));
long double pmax = std::min(options.p*1.03, static_cast<long double>(HEOS.p_critical()+1e-6));
long double pmin = std::max(options.p*0.97, static_cast<long double>(HEOS.p_triple()-1e-6));
Brent(resid, pmin, pmax, LDBL_EPSILON, 1e-8, 100, errstr);
}
}
void SaturationSolvers::saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend *HEOS, long double p, saturation_PHSU_pure_options &options){
void SaturationSolvers::saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend &HEOS, long double p, saturation_PHSU_pure_options &options){
// Define the residual to be driven to zero
class solver_resid : public FuncWrapper1D
@@ -135,8 +136,8 @@ void SaturationSolvers::saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend *HEOS,
long double r, p, rhomolar_liq, rhomolar_vap, value, T, gL, gV;
int other;
solver_resid(HelmholtzEOSMixtureBackend *HEOS, long double p, long double rhomolar_liq_guess, long double rhomolar_vap_guess)
: HEOS(HEOS), p(p), rhomolar_liq(rhomolar_liq_guess), rhomolar_vap(rhomolar_vap_guess){};
solver_resid(HelmholtzEOSMixtureBackend &HEOS, long double p, long double rhomolar_liq_guess, long double rhomolar_vap_guess)
: HEOS(&HEOS), p(p), rhomolar_liq(rhomolar_liq_guess), rhomolar_vap(rhomolar_vap_guess){};
double call(double T){
this->T = T;
// Recalculate the densities using the current guess values
@@ -160,12 +161,12 @@ void SaturationSolvers::saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend *HEOS,
if (!ValidNumber(options.rhoV)){throw ValueError("options.rhoV is not valid in saturation_P_pure_1D_T");};
std::string errstr;
long double Tmax = std::min(options.T + 2, static_cast<long double>(HEOS->T_critical()-1e-6));
long double Tmin = std::max(options.T - 2, static_cast<long double>(HEOS->Ttriple()+1e-6));
long double Tmax = std::min(options.T + 2, static_cast<long double>(HEOS.T_critical()-1e-6));
long double Tmin = std::max(options.T - 2, static_cast<long double>(HEOS.Ttriple()+1e-6));
Brent(resid, Tmin, Tmax, LDBL_EPSILON, 1e-11, 100, errstr);
}
void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, long double specified_value, saturation_PHSU_pure_options &options)
void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, long double specified_value, saturation_PHSU_pure_options &options)
{
/*
This function is inspired by the method of Akasaka:
@@ -179,10 +180,10 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l
std::vector<long double> negativer(3,_HUGE), v;
std::vector<std::vector<long double> > J(3, std::vector<long double>(3,_HUGE));
HEOS->calc_reducing_state();
const SimpleState & reduce = HEOS->get_reducing_state();
shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS->SatL,
SatV = HEOS->SatV;
HEOS.calc_reducing_state();
const SimpleState & reduce = HEOS.get_reducing_state();
shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS.SatL,
SatV = HEOS.SatV;
long double T, rhoL, rhoV, pL, pV;
long double deltaL=0, deltaV=0, tau=0, error;
@@ -196,7 +197,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l
// Invert liquid density ancillary to get temperature
// TODO: fit inverse ancillaries too
try{
T = HEOS->get_components()[0]->ancillaries.pL.invert(specified_value);
T = HEOS.get_components()[0]->ancillaries.pL.invert(specified_value);
}
catch(std::exception &e)
{
@@ -207,11 +208,11 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l
{
throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid",options.specified_variable));
}
if (T > HEOS->T_critical()-1){ T -= 1; }
if (T > HEOS.T_critical()-1){ T -= 1; }
// Evaluate densities from the ancillary equations
rhoV = HEOS->get_components()[0]->ancillaries.rhoV.evaluate(T);
rhoL = HEOS->get_components()[0]->ancillaries.rhoL.evaluate(T);
rhoV = HEOS.get_components()[0]->ancillaries.rhoV.evaluate(T);
rhoL = HEOS.get_components()[0]->ancillaries.rhoL.evaluate(T);
// Apply a single step of Newton's method to improve guess value for liquid
// based on the error between the gas pressure (which is usually very close already)
@@ -372,7 +373,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l
}
while (error > 1e-9);
}
void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long double rhomolar, saturation_D_pure_options &options)
void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend &HEOS, long double rhomolar, saturation_D_pure_options &options)
{
/*
This function is inspired by the method of Akasaka:
@@ -386,10 +387,10 @@ void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long
std::vector<long double> r(2,_HUGE), v;
std::vector<std::vector<long double> > J(2, std::vector<long double>(2,_HUGE));
HEOS->calc_reducing_state();
const SimpleState & reduce = HEOS->get_reducing_state();
shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS->SatL,
SatV = HEOS->SatV;
HEOS.calc_reducing_state();
const SimpleState & reduce = HEOS.get_reducing_state();
shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS.SatL,
SatV = HEOS.SatV;
long double T, rhoL,rhoV;
long double deltaL=0, deltaV=0, tau=0, error, p_error;
@@ -402,16 +403,16 @@ void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long
{
// Invert liquid density ancillary to get temperature
// TODO: fit inverse ancillaries too
T = HEOS->get_components()[0]->ancillaries.rhoL.invert(rhomolar);
rhoV = HEOS->get_components()[0]->ancillaries.rhoV.evaluate(T);
T = HEOS.get_components()[0]->ancillaries.rhoL.invert(rhomolar);
rhoV = HEOS.get_components()[0]->ancillaries.rhoV.evaluate(T);
rhoL = rhomolar;
}
else if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOV)
{
// Invert vapor density ancillary to get temperature
// TODO: fit inverse ancillaries too
T = HEOS->get_components()[0]->ancillaries.rhoV.invert(rhomolar);
rhoL = HEOS->get_components()[0]->ancillaries.rhoL.evaluate(T);
T = HEOS.get_components()[0]->ancillaries.rhoV.invert(rhomolar);
rhoL = HEOS.get_components()[0]->ancillaries.rhoL.evaluate(T);
rhoV = rhomolar;
}
else
@@ -532,7 +533,7 @@ void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long
throw SolutionError(format("saturation_D_pure solver abs error on p [%Lg] > limit [%Lg]", p_error, p_error_limit));
}
}
void SaturationSolvers::saturation_T_pure(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_options &options)
void SaturationSolvers::saturation_T_pure(HelmholtzEOSMixtureBackend &HEOS, long double T, saturation_T_pure_options &options)
{
// Set some imput options
SaturationSolvers::saturation_T_pure_Akasaka_options _options;
@@ -552,7 +553,7 @@ void SaturationSolvers::saturation_T_pure(HelmholtzEOSMixtureBackend *HEOS, long
SaturationSolvers::saturation_T_pure_1D_P(HEOS, T, options);
}
}
void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_Akasaka_options &options)
void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend &HEOS, long double T, saturation_T_pure_Akasaka_options &options)
{
// Start with the method of Akasaka
@@ -566,11 +567,11 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HE
Ancillary equations are used to get a sensible starting point
*/
HEOS->calc_reducing_state();
const SimpleState & reduce = HEOS->get_reducing_state();
long double R_u = HEOS->calc_gas_constant();
shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS->SatL,
SatV = HEOS->SatV;
HEOS.calc_reducing_state();
const SimpleState & reduce = HEOS.get_reducing_state();
long double R_u = HEOS.calc_gas_constant();
shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS.SatL,
SatV = HEOS.SatV;
long double rhoL,rhoV,JL,JV,KL,KV,dJL,dJV,dKL,dKV;
long double DELTA, deltaL=0, deltaV=0, error, PL, PV, stepL, stepV;
@@ -589,13 +590,13 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HE
// Use the density ancillary function as the starting point for the solver
// If very close to the critical temp, evaluate the ancillaries for a slightly lower temperature
if (T > 0.99*HEOS->get_reducing_state().T){
rhoL = HEOS->get_components()[0]->ancillaries.rhoL.evaluate(T-0.1);
rhoV = HEOS->get_components()[0]->ancillaries.rhoV.evaluate(T-0.1);
if (T > 0.99*HEOS.get_reducing_state().T){
rhoL = HEOS.get_components()[0]->ancillaries.rhoL.evaluate(T-0.1);
rhoV = HEOS.get_components()[0]->ancillaries.rhoV.evaluate(T-0.1);
}
else{
rhoL = HEOS->get_components()[0]->ancillaries.rhoL.evaluate(T);
rhoV = HEOS->get_components()[0]->ancillaries.rhoV.evaluate(T);
rhoL = HEOS.get_components()[0]->ancillaries.rhoL.evaluate(T);
rhoV = HEOS.get_components()[0]->ancillaries.rhoV.evaluate(T);
// Apply a single step of Newton's method to improve guess value for liquid
// based on the error between the gas pressure (which is usually very close already)
@@ -738,16 +739,16 @@ void SaturationSolvers::successive_substitution(HelmholtzEOSMixtureBackend &HEOS
for (std::size_t i = 0; i < N; ++i)
{
ln_phi_liq[i] = HEOS.SatL->mixderiv_ln_fugacity_coefficient(i);
ln_phi_vap[i] = HEOS.SatV->mixderiv_ln_fugacity_coefficient(i);
ln_phi_liq[i] = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatL.get()), i);
ln_phi_vap[i] = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatV.get()), i);
if (options.sstype == imposed_p){
deriv_liq = HEOS.SatL->mixderiv_dln_fugacity_coefficient_dT__constp_n(i);
deriv_vap = HEOS.SatV->mixderiv_dln_fugacity_coefficient_dT__constp_n(i);
deriv_liq = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(*(HEOS.SatL.get()), i);
deriv_vap = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(*(HEOS.SatV.get()), i);
}
else if (options.sstype == imposed_T){
deriv_liq = HEOS.SatL->mixderiv_dln_fugacity_coefficient_dp__constT_n(i);
deriv_vap = HEOS.SatV->mixderiv_dln_fugacity_coefficient_dp__constT_n(i);
deriv_liq = MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(*(HEOS.SatL.get()), i);
deriv_vap = MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(*(HEOS.SatV.get()), i);
}
else {throw ValueError();}
@@ -811,15 +812,15 @@ void SaturationSolvers::newton_raphson_VLE_GV::resize(unsigned int N)
// Last entry is 1
neg_dFdS[N+1] = 1.0;
}
void SaturationSolvers::newton_raphson_VLE_GV::check_Jacobian(HelmholtzEOSMixtureBackend *HEOS, const std::vector<long double> &z, std::vector<long double> &K, mixture_VLE_IO &IO)
void SaturationSolvers::newton_raphson_VLE_GV::check_Jacobian(HelmholtzEOSMixtureBackend &HEOS, const std::vector<long double> &z, std::vector<long double> &K, mixture_VLE_IO &IO)
{
// Reset all the variables and resize
pre_call();
std::size_t N = K.size();
resize(N);
shared_ptr<HelmholtzEOSMixtureBackend> SatL(new HelmholtzEOSMixtureBackend(HEOS->get_components())),
SatV(new HelmholtzEOSMixtureBackend(HEOS->get_components()));
shared_ptr<HelmholtzEOSMixtureBackend> SatL(new HelmholtzEOSMixtureBackend(HEOS.get_components())),
SatV(new HelmholtzEOSMixtureBackend(HEOS.get_components()));
SatL->specify_phase(iphase_liquid);
SatV->specify_phase(iphase_gas);
@@ -878,7 +879,7 @@ void SaturationSolvers::newton_raphson_VLE_GV::check_Jacobian(HelmholtzEOSMixtur
std::cout << vec_to_string(get_col(J,N+1),"%12.11f") << std::endl;
}
void SaturationSolvers::newton_raphson_VLE_GV::call(HelmholtzEOSMixtureBackend *HEOS, const std::vector<long double> &z, std::vector<long double> &K, mixture_VLE_IO &IO)
void SaturationSolvers::newton_raphson_VLE_GV::call(HelmholtzEOSMixtureBackend &HEOS, const std::vector<long double> &z, std::vector<long double> &K, mixture_VLE_IO &IO)
{
int iter = 0;
@@ -886,8 +887,8 @@ void SaturationSolvers::newton_raphson_VLE_GV::call(HelmholtzEOSMixtureBackend *
pre_call();
resize(K.size());
shared_ptr<HelmholtzEOSMixtureBackend> SatL(new HelmholtzEOSMixtureBackend(HEOS->get_components())),
SatV(new HelmholtzEOSMixtureBackend(HEOS->get_components()));
shared_ptr<HelmholtzEOSMixtureBackend> SatL(new HelmholtzEOSMixtureBackend(HEOS.get_components())),
SatV(new HelmholtzEOSMixtureBackend(HEOS.get_components()));
SatL->specify_phase(iphase_liquid); // So it will always just use single-phase solution
SatV->specify_phase(iphase_gas); // So it will always just use single-phase solution
@@ -931,7 +932,7 @@ void SaturationSolvers::newton_raphson_VLE_GV::call(HelmholtzEOSMixtureBackend *
IO.y = &y; // Mole fractions in vapor
}
void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureBackend *HEOS, long double beta, long double T, long double rhomolar_liq, const long double rhomolar_vap, const std::vector<long double> &z, std::vector<long double> &K)
void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureBackend &HEOS, long double beta, long double T, long double rhomolar_liq, const long double rhomolar_vap, const std::vector<long double> &z, std::vector<long double> &K)
{
// Step 0:
// --------
@@ -941,6 +942,8 @@ void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureB
// Set the mole fractions in the classes
SatL->set_mole_fractions(x);
SatV->set_mole_fractions(y);
HelmholtzEOSMixtureBackend &rSatV = *SatV, &rSatL = *SatL;
// Update the liquid and vapor classes
SatL->update(DmolarT_INPUTS, rhomolar_liq, T);
@@ -958,22 +961,22 @@ void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureB
// For the residuals F_i
for (unsigned int i = 0; i < N; ++i)
{
long double ln_phi_liq = SatL->mixderiv_ln_fugacity_coefficient(i);
long double phi_iT_liq = SatL->mixderiv_dln_fugacity_coefficient_dT__constrho_n(i);
dlnphi_drho_liq[i] = SatL->mixderiv_dln_fugacity_coefficient_drho__constT_n(i);
long double ln_phi_liq = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatL.get()), i);
long double phi_iT_liq = MixtureDerivatives::dln_fugacity_coefficient_dT__constrho_n(*(HEOS.SatL.get()), i);
dlnphi_drho_liq[i] = MixtureDerivatives::dln_fugacity_coefficient_drho__constT_n(*(HEOS.SatL.get()), i);
for (unsigned int j = 0; j < N; ++j)
{
// I think this is wrong.
phi_ij_liq[j] = SatL->mixderiv_ndln_fugacity_coefficient_dnj__constT_p(i,j) + (SatL->mixderiv_partial_molar_volume(i)/(SatL->gas_constant()*T)-1/p)*SatL->mixderiv_ndpdni__constT_V_nj(i); // 7.126 from GERG monograph
phi_ij_liq[j] = MixtureDerivatives::ndln_fugacity_coefficient_dnj__constT_p(rSatL, i, j) + (MixtureDerivatives::partial_molar_volume(rSatL, i)/(SatL->gas_constant()*T)-1/p)*MixtureDerivatives::ndpdni__constT_V_nj(rSatL, i); // 7.126 from GERG monograph
}
long double ln_phi_vap = SatV->mixderiv_ln_fugacity_coefficient(i);
long double phi_iT_vap = SatV->mixderiv_dln_fugacity_coefficient_dT__constrho_n(i);
dlnphi_drho_vap[i] = SatV->mixderiv_dln_fugacity_coefficient_drho__constT_n(i);
long double ln_phi_vap = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatV.get()), i);
long double phi_iT_vap = MixtureDerivatives::dln_fugacity_coefficient_dT__constrho_n(*(HEOS.SatV.get()), i);
dlnphi_drho_vap[i] = MixtureDerivatives::dln_fugacity_coefficient_drho__constT_n(*(HEOS.SatV.get()), i);
for (unsigned int j = 0; j < N; ++j)
{
// I think this is wrong.
phi_ij_vap[j] = SatV->mixderiv_ndln_fugacity_coefficient_dnj__constT_p(i,j) + (SatV->mixderiv_partial_molar_volume(i)/(SatV->gas_constant()*T)-1/p)*SatV->mixderiv_ndpdni__constT_V_nj(i); ; // 7.126 from GERG monograph
phi_ij_vap[j] = MixtureDerivatives::ndln_fugacity_coefficient_dnj__constT_p(rSatV, i,j) + (MixtureDerivatives::partial_molar_volume(rSatV, i)/(SatV->gas_constant()*T)-1/p)*MixtureDerivatives::ndpdni__constT_V_nj(rSatV, i); ; // 7.126 from GERG monograph
}
r[i] = log(K[i]) + ln_phi_vap - ln_phi_liq;
@@ -1009,12 +1012,12 @@ void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureB
r[N+1] = p_liq-p_vap;
for (unsigned int j = 0; j < N; ++j)
{
J[N+1][j] = HEOS->gas_constant()*T*K[j]*z[j]/pow(1-beta+beta*K[j],(int)2)*((1-beta)*dlnphi_drho_vap[j]+beta*dlnphi_drho_liq[j]);
J[N+1][j] = HEOS.gas_constant()*T*K[j]*z[j]/pow(1-beta+beta*K[j],(int)2)*((1-beta)*dlnphi_drho_vap[j]+beta*dlnphi_drho_liq[j]);
}
// dF_{N+1}/d(ln(T))
J[N+1][N] = T*(SatL->mixderiv_dpdT__constV_n() - SatV->mixderiv_dpdT__constV_n());
J[N+1][N] = T*(MixtureDerivatives::dpdT__constV_n(*(HEOS.SatL.get())) - MixtureDerivatives::dpdT__constV_n(*(HEOS.SatV.get())));
// dF_{N+1}/d(ln(rho'))
J[N+1][N+1] = rhomolar_liq*SatL->mixderiv_dpdrho__constT_n();
J[N+1][N+1] = rhomolar_liq*MixtureDerivatives::dpdrho__constT_n(*(HEOS.SatL.get()));
// Flip all the signs of the entries in the residual vector since we are solving Jv = -r, not Jv=r
// Also calculate the rms error of the residual vector at this step
@@ -1027,7 +1030,7 @@ void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureB
error_rms = sqrt(error_rms); // Square-root (The R in RMS)
}
void PhaseEnvelope::PhaseEnvelope_GV::build(HelmholtzEOSMixtureBackend *HEOS, const std::vector<long double> &z, std::vector<long double> &K, SaturationSolvers::mixture_VLE_IO &IO)
void PhaseEnvelope::PhaseEnvelope_GV::build(HelmholtzEOSMixtureBackend &HEOS, const std::vector<long double> &z, std::vector<long double> &K, SaturationSolvers::mixture_VLE_IO &IO)
{
// Use the residual function based on ln(K_i), ln(T) and ln(rho') as independent variables. rho'' is specified
SaturationSolvers::newton_raphson_VLE_GV NRVLE;

40
src/Backends/Helmholtz/VLERoutines.h
View File

@@ -45,14 +45,14 @@ namespace SaturationSolvers
@param p Pressure [Pa]
@param i Index of component [-]
*/
static long double Wilson_lnK_factor(HelmholtzEOSMixtureBackend *HEOS, long double T, long double p, int i){
EquationOfState *EOS = (HEOS->get_components())[i]->pEOS;
static long double Wilson_lnK_factor(HelmholtzEOSMixtureBackend &HEOS, long double T, long double p, int i){
EquationOfState *EOS = (HEOS.get_components())[i]->pEOS;
return log(EOS->reduce.p/p)+5.373*(1 + EOS->accentric)*(1-EOS->reduce.T/T);
};
void saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long double rhomolar, saturation_D_pure_options &options);
void saturation_T_pure(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_options &options);
void saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_Akasaka_options &options);
void saturation_D_pure(HelmholtzEOSMixtureBackend &HEOS, long double rhomolar, saturation_D_pure_options &options);
void saturation_T_pure(HelmholtzEOSMixtureBackend &HEOS, long double T, saturation_T_pure_options &options);
void saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend &HEOS, long double T, saturation_T_pure_Akasaka_options &options);
/**
*/
@@ -67,7 +67,7 @@ namespace SaturationSolvers
/**
*/
void saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, long double specified_value, saturation_PHSU_pure_options &options);
void saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, long double specified_value, saturation_PHSU_pure_options &options);
/* \brief This is a backup saturation_p solver for the case where the Newton solver cannot approach closely enough the solution
*
@@ -77,7 +77,7 @@ namespace SaturationSolvers
* @param p Imposed pressure in kPa
* @param options Options to be passed to the function (at least T, rhoL and rhoV must be provided)
*/
void saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend *HEOS, long double p, saturation_PHSU_pure_options &options);
void saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend &HEOS, long double p, saturation_PHSU_pure_options &options);
/* \brief This is a backup saturation_T solver for the case where the Newton solver cannot approach closely enough the solution
*
@@ -87,7 +87,7 @@ namespace SaturationSolvers
* @param T Imposed temperature in K
* @param options Options to be passed to the function (at least p, rhoL and rhoV must be provided)
*/
void saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_options &options);
void saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend &HEOS, long double T, saturation_T_pure_options &options);
/* \brief A robust but slow solver in the very-near-critical region
*
@@ -102,7 +102,7 @@ namespace SaturationSolvers
* @param ykey The CoolProp::parameters key to be imposed - one of iT or iP
* @param y The value for the imposed variable
*/
void saturation_critical(HelmholtzEOSMixtureBackend *HEOS, CoolProp::parameters ykey, long double y);
void saturation_critical(HelmholtzEOSMixtureBackend &HEOS, CoolProp::parameters ykey, long double y);
void successive_substitution(HelmholtzEOSMixtureBackend &HEOS,
const long double beta,
@@ -127,8 +127,8 @@ namespace SaturationSolvers
std::vector<long double> *K;
HelmholtzEOSMixtureBackend *HEOS;
WilsonK_resid(HelmholtzEOSMixtureBackend *HEOS, double beta, double imposed_value, int input_type, const std::vector<long double> &z, std::vector<long double> &K){
this->z = &z; this->K = &K; this->HEOS = HEOS; this->beta = beta; this->input_type = input_type;
WilsonK_resid(HelmholtzEOSMixtureBackend &HEOS, double beta, double imposed_value, int input_type, const std::vector<long double> &z, std::vector<long double> &K){
this->z = &z; this->K = &K; this->HEOS = &HEOS; this->beta = beta; this->input_type = input_type;
if (input_type == imposed_T){
this->T = imposed_value;
}
@@ -145,19 +145,19 @@ namespace SaturationSolvers
T = input_value; // Iterate on temperature, pressure imposed
}
for (unsigned int i = 0; i< (*z).size(); i++) {
(*K)[i] = exp(Wilson_lnK_factor(HEOS,T,p,i));
(*K)[i] = exp(Wilson_lnK_factor(*HEOS,T,p,i));
summer += (*z)[i]*((*K)[i]-1)/(1-beta+beta*(*K)[i]);
}
return summer;
};
};
inline double saturation_preconditioner(HelmholtzEOSMixtureBackend *HEOS, double input_value, int input_type, const std::vector<long double> &z)
inline double saturation_preconditioner(HelmholtzEOSMixtureBackend &HEOS, double input_value, int input_type, const std::vector<long double> &z)
{
double ptriple = 0, pcrit = 0, Ttriple = 0, Tcrit = 0;
for (unsigned int i = 0; i < z.size(); i++)
{
EquationOfState *EOS = (HEOS->get_components())[i]->pEOS;
EquationOfState *EOS = (HEOS.get_components())[i]->pEOS;
ptriple += EOS->sat_min_liquid.p*z[i];
pcrit += EOS->reduce.p*z[i];
@@ -175,14 +175,14 @@ namespace SaturationSolvers
}
else{ throw ValueError();}
}
inline double saturation_Wilson(HelmholtzEOSMixtureBackend *HEOS, double beta, double input_value, int input_type, const std::vector<long double> &z, double guess)
inline double saturation_Wilson(HelmholtzEOSMixtureBackend &HEOS, double beta, double input_value, int input_type, const std::vector<long double> &z, double guess)
{
double T;
std::string errstr;
// Find first guess for T using Wilson K-factors
WilsonK_resid Resid(HEOS, beta, input_value, input_type, z, HEOS->get_K());
WilsonK_resid Resid(HEOS, beta, input_value, input_type, z, HEOS.get_K());
T = Secant(Resid, guess, 0.001, 1e-10, 100, errstr);
if (!ValidNumber(T)){throw ValueError("saturation_p_Wilson failed to get good T");}
@@ -234,7 +234,7 @@ namespace SaturationSolvers
@param z Bulk mole fractions [-]
@param K Array of K-factors [-]
*/
void call(HelmholtzEOSMixtureBackend *HEOS, const std::vector<long double> &z, std::vector<long double> &K, mixture_VLE_IO &IO);
void call(HelmholtzEOSMixtureBackend &HEOS, const std::vector<long double> &z, std::vector<long double> &K, mixture_VLE_IO &IO);
/*! Build the arrays for the Newton-Raphson solve
@@ -246,11 +246,11 @@ namespace SaturationSolvers
@param z Bulk mole fractions [-]
@param K Array of K-factors [-]
*/
void build_arrays(HelmholtzEOSMixtureBackend *HEOS, long double beta, long double T, long double rhomolar_liq, const long double rho_vapor, const std::vector<long double> &z, std::vector<long double> & K);
void build_arrays(HelmholtzEOSMixtureBackend &HEOS, long double beta, long double T, long double rhomolar_liq, const long double rho_vapor, const std::vector<long double> &z, std::vector<long double> & K);
/** Check the derivatives in the Jacobian using numerical derivatives.
*/
void check_Jacobian(HelmholtzEOSMixtureBackend *HEOS, const std::vector<long double> &z, std::vector<long double> &K, mixture_VLE_IO &IO);
void check_Jacobian(HelmholtzEOSMixtureBackend &HEOS, const std::vector<long double> &z, std::vector<long double> &K, mixture_VLE_IO &IO);
};
};
@@ -299,7 +299,7 @@ namespace PhaseEnvelope
public:
PhaseEnvelopeLog bubble, dew;
void build(HelmholtzEOSMixtureBackend *HEOS, const std::vector<long double> &z, std::vector<long double> &K, SaturationSolvers::mixture_VLE_IO &IO);
void build(HelmholtzEOSMixtureBackend &HEOS, const std::vector<long double> &z, std::vector<long double> &K, SaturationSolvers::mixture_VLE_IO &IO);
};
};