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

Implemented the x_N dependent form of the composition derivatives. Mostly working but derivatives with respect to x_N DO NOT work

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-08-28 01:22:10 +02:00
parent b518684404
commit e7aa3df917
8 changed files with 725 additions and 375 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -1993,7 +1993,8 @@ long double HelmholtzEOSMixtureBackend::calc_gibbsmolar(void)
}
long double HelmholtzEOSMixtureBackend::calc_fugacity_coefficient(int i)
{
return exp(MixtureDerivatives::ln_fugacity_coefficient(*this, i));
x_N_dependency_flag xN_flag = XN_DEPENDENT;
return exp(MixtureDerivatives::ln_fugacity_coefficient(*this, i, xN_flag));
}
SimpleState HelmholtzEOSMixtureBackend::calc_reducing_state_nocache(const std::vector<long double> & mole_fractions)

436
src/Backends/Helmholtz/MixtureDerivatives.cpp
View File

@@ -2,113 +2,191 @@
namespace CoolProp{
long double MixtureDerivatives::dalphar_dxi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i)
long double MixtureDerivatives::dalphar_dxi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
return HEOS.components[i]->pEOS->baser(HEOS._tau, HEOS._delta) + HEOS.Excess.dalphar_dxi(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i);
if (xN_flag == XN_INDEPENDENT){
return HEOS.components[i]->pEOS->baser(HEOS._tau, HEOS._delta) + HEOS.Excess.dalphar_dxi(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i);
}
else if(xN_flag == XN_DEPENDENT){
std::vector<long double> &x = HEOS.mole_fractions;
std::size_t N = x.size();
if (i==N-1) return _HUGE;
double dar_dxi = HEOS.components[i]->pEOS->baser(HEOS._tau, HEOS._delta) - HEOS.components[N-1]->pEOS->baser(HEOS._tau, HEOS._delta);
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->alphar(HEOS._tau, HEOS._delta);
dar_dxi += (1-2*x[i])*FiNariN;
for (std::size_t k = 0; k < N-1; ++k){
if (i==k) continue;
double Fikarik = HEOS.Excess.F[i][k]*HEOS.Excess.DepartureFunctionMatrix[i][k]->alphar(HEOS._tau, HEOS._delta);
double FkNarkN = HEOS.Excess.F[k][N-1]*HEOS.Excess.DepartureFunctionMatrix[k][N-1]->alphar(HEOS._tau, HEOS._delta);
dar_dxi += x[i]*(Fikarik - FiNariN - FkNarkN);
}
return dar_dxi;
}
else{
throw ValueError(format("xN_flag is invalid"));
}
}
long double MixtureDerivatives::d2alphar_dxi_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i)
long double MixtureDerivatives::d2alphar_dxi_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
return HEOS.components[i]->pEOS->dalphar_dTau(HEOS._tau, HEOS._delta) + HEOS.Excess.d2alphar_dxi_dTau(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i);
if (xN_flag == XN_INDEPENDENT){
return HEOS.components[i]->pEOS->dalphar_dTau(HEOS._tau, HEOS._delta) + HEOS.Excess.d2alphar_dxi_dTau(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i);
}
else if(xN_flag == XN_DEPENDENT){
std::vector<long double> &x = HEOS.mole_fractions;
std::size_t N = x.size();
if (i==N-1) return _HUGE;
double d2ar_dxi_dTau = HEOS.components[i]->pEOS->dalphar_dTau(HEOS._tau, HEOS._delta) - HEOS.components[N-1]->pEOS->dalphar_dTau(HEOS._tau, HEOS._delta);
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->dalphar_dTau(HEOS._tau, HEOS._delta);
d2ar_dxi_dTau += (1-2*x[i])*FiNariN;
for (std::size_t k = 0; k < N-1; ++k){
if (i==k) continue;
double Fikarik = HEOS.Excess.F[i][k]*HEOS.Excess.DepartureFunctionMatrix[i][k]->dalphar_dTau(HEOS._tau, HEOS._delta);
double FkNarkN = HEOS.Excess.F[k][N-1]*HEOS.Excess.DepartureFunctionMatrix[k][N-1]->dalphar_dTau(HEOS._tau, HEOS._delta);
d2ar_dxi_dTau += x[i]*(Fikarik - FiNariN - FkNarkN);
}
return d2ar_dxi_dTau;
}
else{
throw ValueError(format("xN_flag is invalid"));
}
}
long double MixtureDerivatives::d2alphar_dxi_dDelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i)
long double MixtureDerivatives::d2alphar_dxi_dDelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
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, std::size_t i, std::size_t j)
{
return 0 + HEOS.Excess.d2alphardxidxj(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i, j);
if (xN_flag == XN_INDEPENDENT){
return HEOS.components[i]->pEOS->dalphar_dDelta(HEOS._tau, HEOS._delta) + HEOS.Excess.d2alphar_dxi_dDelta(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i);
}
else if(xN_flag == XN_DEPENDENT){
std::vector<long double> &x = HEOS.mole_fractions;
std::size_t N = x.size();
if (i==N-1) return _HUGE;
double d2ar_dxi_dDelta = HEOS.components[i]->pEOS->dalphar_dDelta(HEOS._tau, HEOS._delta) - HEOS.components[N-1]->pEOS->dalphar_dDelta(HEOS._tau, HEOS._delta);
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->dalphar_dDelta(HEOS._tau, HEOS._delta);
d2ar_dxi_dDelta += (1-2*x[i])*FiNariN;
for (std::size_t k = 0; k < N-1; ++k){
if (i==k) continue;
double Fikarik = HEOS.Excess.F[i][k]*HEOS.Excess.DepartureFunctionMatrix[i][k]->dalphar_dDelta(HEOS._tau, HEOS._delta);
double FkNarkN = HEOS.Excess.F[k][N-1]*HEOS.Excess.DepartureFunctionMatrix[k][N-1]->dalphar_dDelta(HEOS._tau, HEOS._delta);
d2ar_dxi_dDelta += x[i]*(Fikarik - FiNariN - FkNarkN);
}
return d2ar_dxi_dDelta;
}
else{
throw ValueError(format("xN_flag is invalid"));
}
}
long double MixtureDerivatives::fugacity_i(HelmholtzEOSMixtureBackend &HEOS, std::size_t i){
return HEOS.mole_fractions[i]*HEOS.rhomolar()*HEOS.gas_constant()*HEOS.T()*exp( dnalphar_dni__constT_V_nj(HEOS, i));
}
long double MixtureDerivatives::ln_fugacity_coefficient(HelmholtzEOSMixtureBackend &HEOS, std::size_t i)
long double MixtureDerivatives::d2alphardxidxj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
return HEOS.alphar() + ndalphar_dni__constT_V_nj(HEOS, i)-log(1+HEOS._delta.pt()*HEOS.dalphar_dDelta());
if (xN_flag == XN_INDEPENDENT){
return 0 + HEOS.Excess.d2alphardxidxj(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i, j);
}
else if(xN_flag == XN_DEPENDENT){
std::size_t N = HEOS.mole_fractions.size();
if (i == N-1){ return 0;}
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->alphar(HEOS._tau, HEOS._delta);
if (i == j) { return FiNariN; }
if (j == N-1){ return _HUGE; }
double Fijarij = HEOS.Excess.F[i][j]*HEOS.Excess.DepartureFunctionMatrix[i][j]->alphar(HEOS._tau, HEOS._delta);
double FjNarjN = HEOS.Excess.F[j][N-1]*HEOS.Excess.DepartureFunctionMatrix[j][N-1]->alphar(HEOS._tau, HEOS._delta);
return Fijarij-FiNariN-FjNarjN;
}
else{
throw ValueError(format("xN_flag is invalid"));
}
}
long double MixtureDerivatives::dln_fugacity_coefficient_dT__constrho_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i)
long double MixtureDerivatives::fugacity_i(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag){
return HEOS.mole_fractions[i]*HEOS.rhomolar()*HEOS.gas_constant()*HEOS.T()*exp( dnalphar_dni__constT_V_nj(HEOS, i, xN_flag));
}
long double MixtureDerivatives::ln_fugacity_coefficient(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
return HEOS.alphar() + ndalphar_dni__constT_V_nj(HEOS, i, xN_flag)-log(1+HEOS._delta.pt()*HEOS.dalphar_dDelta());
}
long double MixtureDerivatives::dln_fugacity_coefficient_dT__constrho_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
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;
return (HEOS.dalphar_dTau() + d_ndalphardni_dTau(HEOS, i, xN_flag)-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, std::size_t i)
long double MixtureDerivatives::dln_fugacity_coefficient_drho__constT_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
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;
return (HEOS.dalphar_dDelta() + d_ndalphardni_dDelta(HEOS, i, xN_flag)-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, std::size_t i)
long double MixtureDerivatives::dnalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
// GERG Equation 7.42
return HEOS.alphar() + ndalphar_dni__constT_V_nj(HEOS, i);
return HEOS.alphar() + ndalphar_dni__constT_V_nj(HEOS, i, xN_flag);
}
long double MixtureDerivatives::d2nalphar_dni_dT(HelmholtzEOSMixtureBackend &HEOS, std::size_t i)
long double MixtureDerivatives::d2nalphar_dni_dT(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
return -HEOS._tau.pt()/HEOS._T*(HEOS.dalphar_dTau() + d_ndalphardni_dTau(HEOS, i));
return -HEOS._tau.pt()/HEOS._T*(HEOS.dalphar_dTau() + d_ndalphardni_dTau(HEOS, i, xN_flag));
}
long double MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i)
long double MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
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);
return d2nalphar_dni_dT(HEOS, i, xN_flag) + 1/T-partial_molar_volume(HEOS, i, xN_flag)/(R_u*T)*dpdT__constV_n(HEOS);
}
long double MixtureDerivatives::partial_molar_volume(HelmholtzEOSMixtureBackend &HEOS, std::size_t i)
long double MixtureDerivatives::partial_molar_volume(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
return -ndpdni__constT_V_nj(HEOS, i)/ndpdV__constT_n(HEOS);
return -ndpdni__constT_V_nj(HEOS, i, xN_flag)/ndpdV__constT_n(HEOS);
}
long double MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i)
long double MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
// GERG equation 7.30
long double R_u = HEOS.gas_constant();
double partial_molar_volumeval = partial_molar_volume(HEOS, i); // [m^3/mol]
double partial_molar_volumeval = partial_molar_volume(HEOS, i, xN_flag); // [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, std::size_t i, std::size_t j)
long double MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
// 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_dxj_dni__constT_V(HEOS, j, i) - partial_molar_volume(HEOS, i)/(R_u*HEOS._T)*dpdxj__constT_V_xi(HEOS, j);
return d2nalphar_dxj_dni__constT_V(HEOS, j, i, xN_flag) - partial_molar_volume(HEOS, i, xN_flag)/(R_u*HEOS._T)*dpdxj__constT_V_xi(HEOS, j, xN_flag);
}
long double MixtureDerivatives::dpdxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j)
long double MixtureDerivatives::dpdxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag)
{
// 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));
return HEOS._rhomolar*R_u*HEOS._T*(ddelta_dxj__constT_V_xi(HEOS, j, xN_flag)*HEOS.dalphar_dDelta()+HEOS._delta.pt()*d_dalpharddelta_dxj__constT_V_xi(HEOS, j, xN_flag));
}
long double MixtureDerivatives::d_dalpharddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j)
long double MixtureDerivatives::d_dalpharddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag)
{
// 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);
return HEOS.d2alphar_dDelta2()*ddelta_dxj__constT_V_xi(HEOS, j, xN_flag)
+ HEOS.d2alphar_dDelta_dTau()*dtau_dxj__constT_V_xi(HEOS, j, xN_flag)
+ d2alphar_dxi_dDelta(HEOS, j, xN_flag);
}
long double MixtureDerivatives::dalphar_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j)
long double MixtureDerivatives::dalphar_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag)
{
//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);
return HEOS.dalphar_dDelta()*ddelta_dxj__constT_V_xi(HEOS, j, xN_flag)+HEOS.dalphar_dTau()*dtau_dxj__constT_V_xi(HEOS, j, xN_flag)+dalphar_dxi(HEOS, j, xN_flag);
}
long double MixtureDerivatives::d_ndalphardni_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j)
long double MixtureDerivatives::d_ndalphardni_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
// 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);
return d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i,j, xN_flag)
+ ddelta_dxj__constT_V_xi(HEOS, j, xN_flag)*d_ndalphardni_dDelta(HEOS, i, xN_flag)
+ dtau_dxj__constT_V_xi(HEOS, j, xN_flag)*d_ndalphardni_dTau(HEOS, i, xN_flag);
}
long double MixtureDerivatives::ddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j)
long double MixtureDerivatives::ddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag)
{
// Gernert 3.121 (Catch test provided)
return -HEOS._delta.pt()/HEOS._reducing.rhomolar*HEOS.Reducing.p->drhormolardxi__constxj(HEOS.mole_fractions,j);
return -HEOS._delta.pt()/HEOS._reducing.rhomolar*HEOS.Reducing.p->drhormolardxi__constxj(HEOS.mole_fractions,j, xN_flag);
}
long double MixtureDerivatives::dtau_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j)
long double MixtureDerivatives::dtau_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag)
{
// Gernert 3.122 (Catch test provided)
return 1/HEOS._T*HEOS.Reducing.p->dTrdxi__constxj(HEOS.mole_fractions,j);
return 1/HEOS._T*HEOS.Reducing.p->dTrdxi__constxj(HEOS.mole_fractions,j, xN_flag);
}
long double MixtureDerivatives::dpdT__constV_n(HelmholtzEOSMixtureBackend &HEOS)
@@ -126,103 +204,115 @@ 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, std::size_t i)
long double MixtureDerivatives::ndpdni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
// 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 ndrhorbar_dni__constnj = HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions,i, xN_flag);
double ndTr_dni__constnj = HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions,i, xN_flag);
double summer = 0;
for (unsigned int k = 0; k < HEOS.mole_fractions.size(); ++k)
std::size_t kmax = HEOS.mole_fractions.size();
if (xN_flag == XN_DEPENDENT){ kmax--; }
for (unsigned int k = 0; k < kmax; ++k)
{
summer += HEOS.mole_fractions[k]*d2alphar_dxi_dDelta(HEOS, k);
summer += HEOS.mole_fractions[k]*d2alphar_dxi_dDelta(HEOS, k, xN_flag);
}
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;
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, xN_flag)-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, std::size_t i)
long double MixtureDerivatives::ndalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
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 term1 = HEOS._delta.pt()*HEOS.dalphar_dDelta()*(1-1/HEOS._reducing.rhomolar*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions,i, xN_flag));
double term2 = HEOS._tau.pt()*HEOS.dalphar_dTau()*(1/HEOS._reducing.T)*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions,i, xN_flag);
double s = 0;
for (unsigned int k = 0; k < HEOS.mole_fractions.size(); k++)
std::size_t kmax = HEOS.mole_fractions.size();
if (xN_flag == XN_DEPENDENT){ kmax--; }
for (unsigned int k = 0; k < kmax; k++)
{
s += HEOS.mole_fractions[k]*dalphar_dxi(HEOS, k);
s += HEOS.mole_fractions[k]*dalphar_dxi(HEOS, k, xN_flag);
}
double term3 = dalphar_dxi(HEOS, i);
double term3 = dalphar_dxi(HEOS, i, xN_flag);
return term1 + term2 + term3 - s;
}
long double MixtureDerivatives::ndln_fugacity_coefficient_dnj__constT_p(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j)
long double MixtureDerivatives::ndln_fugacity_coefficient_dnj__constT_p(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
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);
return nd2nalphardnidnj__constT_V(HEOS, j, i, xN_flag) + 1 - partial_molar_volume(HEOS, j, xN_flag)/(R_u*HEOS._T)*ndpdni__constT_V_nj(HEOS, i, xN_flag);
}
long double MixtureDerivatives::nddeltadni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i)
long double MixtureDerivatives::nddeltadni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
return HEOS._delta.pt()-HEOS._delta.pt()/HEOS._reducing.rhomolar*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions, i);
return HEOS._delta.pt()-HEOS._delta.pt()/HEOS._reducing.rhomolar*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag);
}
long double MixtureDerivatives::ndtaudni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i)
long double MixtureDerivatives::ndtaudni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
return HEOS._tau.pt()/HEOS._reducing.T*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i);
return HEOS._tau.pt()/HEOS._reducing.T*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
}
long double MixtureDerivatives::d_ndalphardni_dxj__constdelta_tau_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j)
long double MixtureDerivatives::d_ndalphardni_dxj__constdelta_tau_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
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 line1 = HEOS._delta.pt()*d2alphar_dxi_dDelta(HEOS, j, xN_flag)*(1-1/HEOS._reducing.rhomolar*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag));
double line2 = -HEOS._delta.pt()*HEOS.dalphar_dDelta()*(1/HEOS._reducing.rhomolar)*(HEOS.Reducing.p->d_ndrhorbardni_dxj__constxi(HEOS.mole_fractions, i, j, xN_flag)-1/HEOS._reducing.rhomolar*HEOS.Reducing.p->drhormolardxi__constxj(HEOS.mole_fractions,j, xN_flag)*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions,i, xN_flag));
double line3 = HEOS._tau.pt()*d2alphar_dxi_dTau(HEOS, j, xN_flag)*(1/HEOS._reducing.T)*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
double line4 = HEOS._tau.pt()*HEOS.dalphar_dTau()*(1/HEOS._reducing.T)*(HEOS.Reducing.p->d_ndTrdni_dxj__constxi(HEOS.mole_fractions,i,j, xN_flag)-1/HEOS._reducing.T*HEOS.Reducing.p->dTrdxi__constxj(HEOS.mole_fractions, j, xN_flag)*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag));
double s = 0;
for (unsigned int m = 0; m < HEOS.mole_fractions.size(); m++)
std::size_t mmax = HEOS.mole_fractions.size();
if (xN_flag == XN_DEPENDENT){ mmax--; }
for (unsigned int m = 0; m < mmax; m++)
{
s += HEOS.mole_fractions[m]*d2alphardxidxj(HEOS, j,m);
s += HEOS.mole_fractions[m]*d2alphardxidxj(HEOS, j,m, xN_flag);
}
double line5 = d2alphardxidxj(HEOS, i,j)-dalphar_dxi(HEOS, j)-s;
double line5 = d2alphardxidxj(HEOS, i,j, xN_flag)-dalphar_dxi(HEOS, j, xN_flag)-s;
return line1+line2+line3+line4+line5;
}
long double MixtureDerivatives::nd2nalphardnidnj__constT_V(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j)
long double MixtureDerivatives::nd2nalphardnidnj__constT_V(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
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 line0 = ndalphar_dni__constT_V_nj(HEOS, j, xN_flag); // First term from 7.46
double line1 = d_ndalphardni_dDelta(HEOS, i, xN_flag)*nddeltadni__constT_V_nj(HEOS, j, xN_flag);
double line2 = d_ndalphardni_dTau(HEOS, i, xN_flag)*ndtaudni__constT_V_nj(HEOS, j, xN_flag);
double summer = 0;
for (unsigned int k = 0; k < HEOS.mole_fractions.size(); k++)
std::size_t kmax = HEOS.mole_fractions.size();
if (xN_flag == XN_DEPENDENT){ kmax--; }
for (unsigned int k = 0; k < kmax; k++)
{
summer += HEOS.mole_fractions[k]*d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, k);
summer += HEOS.mole_fractions[k]*d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, k, xN_flag);
}
double line3 = d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, j)-summer;
double line3 = d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, j, xN_flag)-summer;
return line0 + line1 + line2 + line3;
}
long double MixtureDerivatives::d_ndalphardni_dDelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i)
long double MixtureDerivatives::d_ndalphardni_dDelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
// 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));
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, xN_flag));
// 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);
double term2 = HEOS._tau.pt()*HEOS.d2alphar_dDelta_dTau()*(1/HEOS._reducing.T)*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
// The third line
double term3 = d2alphar_dxi_dDelta(HEOS, i);
for (unsigned int k = 0; k < HEOS.mole_fractions.size(); k++)
double term3 = d2alphar_dxi_dDelta(HEOS, i, xN_flag);
std::size_t kmax = HEOS.mole_fractions.size();
if (xN_flag == XN_DEPENDENT){ kmax--; }
for (unsigned int k = 0; k < kmax; k++)
{
term3 -= HEOS.mole_fractions[k]*d2alphar_dxi_dDelta(HEOS, k);
term3 -= HEOS.mole_fractions[k]*d2alphar_dxi_dDelta(HEOS, k, xN_flag);
}
return term1 + term2 + term3;
}
long double MixtureDerivatives::d_ndalphardni_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i)
long double MixtureDerivatives::d_ndalphardni_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
// 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));
double term1 = HEOS._delta.pt()*HEOS.d2alphar_dDelta_dTau()*(1-1/HEOS._reducing.rhomolar*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag));
// 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);
double term2 = (HEOS._tau.pt()*HEOS.d2alphar_dTau2()+HEOS.dalphar_dTau())*(1/HEOS._reducing.T)*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
// The third line
double term3 = d2alphar_dxi_dTau(HEOS, i);
for (unsigned int k = 0; k < HEOS.mole_fractions.size(); k++)
double term3 = d2alphar_dxi_dTau(HEOS, i, xN_flag);
std::size_t kmax = HEOS.mole_fractions.size();
if (xN_flag == XN_DEPENDENT){ kmax--; }
for (unsigned int k = 0; k < kmax; k++)
{
term3 -= HEOS.mole_fractions[k]*d2alphar_dxi_dTau(HEOS, k);
term3 -= HEOS.mole_fractions[k]*d2alphar_dxi_dTau(HEOS, k, xN_flag);
}
return term1 + term2 + term3;
}
@@ -237,13 +327,14 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
{
/* Set up a test class for the mixture tests */
std::vector<std::string> names(2);
names[0] = "Methane"; names[1] = "Ethane";
names[0] = "Ethane"; names[1] = "Propane";
std::vector<long double> z(2);
z[0] = 0.5; z[1] = 1-z[0];
z[0] = 0.25; z[1] = 1-z[0];
shared_ptr<HelmholtzEOSMixtureBackend> HEOS(new HelmholtzEOSMixtureBackend(names));
HelmholtzEOSMixtureBackend &rHEOS = *(HEOS.get());
rHEOS.set_mole_fractions(z);
x_N_dependency_flag xN_flag = XN_DEPENDENT;
// These ones only require the i index
for (std::size_t i = 0; i< z.size();++i)
{
@@ -255,36 +346,37 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
rHEOS.specify_phase(iphase_gas);
rHEOS.update(PT_INPUTS, 101325, T1);
double analytic = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rHEOS, i);
double analytic = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rHEOS, i, xN_flag);
rHEOS.update(PT_INPUTS, 101325, T1 + dT);
double v1 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i);
double v1 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i, xN_flag);
rHEOS.update(PT_INPUTS, 101325, T1 - dT);
double v2 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i);
double v2 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i, xN_flag);
double numeric = (v1 - v2)/(2*dT);
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-6);
CHECK(err < 1e-8);
}
std::ostringstream ss2;
ss2 << "dln_fugacity_coefficient_dp__constT_n, i=" << i;
SECTION(ss2.str(), "")
{
double p1 = 101325, dP = 1;
double p0 = 101325, dP = 1e-4;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(PT_INPUTS, p1, 300);
double analytic = MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(rHEOS, i);
rHEOS.update(PT_INPUTS, p0, 300);
double analytic = MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(rHEOS, i, xN_flag);
double rho1 = rHEOS.rhomolar();
rHEOS.update(PT_INPUTS, p1 + dP, 300);
double v1 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i);
rHEOS.update(PT_INPUTS, p1 - dP, 300);
double v2 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i);
rHEOS.update(DmolarT_INPUTS, rho1 + dP, 300);
double v1 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i, xN_flag), p1 = rHEOS.p();
rHEOS.update(DmolarT_INPUTS, rho1 - dP, 300);
double v2 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i, xN_flag), p2 = rHEOS.p();
double numeric = (v1 - v2)/(2*dP);
double numeric = (v1 - v2)/(p1 - p2);
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-6);
CHECK(err < 1e-8);
}
std::ostringstream ss3;
ss3 << "d_ndalphardni_dDelta, i=" << i;
@@ -294,16 +386,16 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
rHEOS.specify_phase(iphase_gas);
rHEOS.update(PT_INPUTS, p1, 300);
double analytic = MixtureDerivatives::d_ndalphardni_dDelta(rHEOS, i);
double analytic = MixtureDerivatives::d_ndalphardni_dDelta(rHEOS, i, xN_flag);
rHEOS.update(PT_INPUTS, p1 + dP, 300);
double v1 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS, i), delta1 = rHEOS.delta();
double v1 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS, i, xN_flag), delta1 = rHEOS.delta();
rHEOS.update(PT_INPUTS, p1 - dP, 300);
double v2 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS, i), delta2 = rHEOS.delta();
double v2 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS, i, xN_flag), delta2 = rHEOS.delta();
double numeric = (v1 - v2)/(delta1 - delta2);
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-6);
CHECK(err < 1e-8);
}
std::ostringstream ss4;
ss4 << "d_ndalphardni_dTau, i=" << i;
@@ -315,18 +407,79 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
double rho1 = rHEOS.rhomolar();
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double analytic = MixtureDerivatives::d_ndalphardni_dTau(rHEOS, i);
double analytic = MixtureDerivatives::d_ndalphardni_dTau(rHEOS, i, xN_flag);
rHEOS.update(DmolarT_INPUTS, rho1, 300 + dT);
double v1 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS, i), tau1 = rHEOS.tau();
double v1 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS, i, xN_flag), tau1 = rHEOS.tau();
rHEOS.update(DmolarT_INPUTS, rho1, 300 - dT);
double v2 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS, i), tau2 = rHEOS.tau();
double v2 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS, i, xN_flag), tau2 = rHEOS.tau();
double numeric = (v1 - v2)/(tau1 - tau2);
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-6);
CHECK(err < 1e-8);
}
std::ostringstream ss5;
ss5 << "dpdxj__constT_V_xi, i=" << i;
SECTION(ss5.str(), "")
{
double dz = 1e-6;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(DmolarT_INPUTS, 300, 300);
double rho1 = rHEOS.rhomolar();
double analytic = MixtureDerivatives::dpdxj__constT_V_xi(rHEOS, i, xN_flag);
std::vector<long double> zp = z; /// Copy base composition
zp[i] += dz;
rHEOS.set_mole_fractions(zp);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v1 = rHEOS.p();
std::vector<long double> zm = z; /// Copy base composition
zm[i] -= dz;
rHEOS.set_mole_fractions(zm);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v2 = rHEOS.p();
double numeric = (v1 - v2)/(2*dz);
double err = std::abs((numeric-analytic)/analytic);
CAPTURE(numeric);
CAPTURE(analytic);
CHECK(err < 1e-80);
}
std::ostringstream ss6;
ss6 << "d_dalpharddelta_dxj__constT_V_xi, i=" << i;
SECTION(ss6.str(), "")
{
double dz = 1e-6;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(DmolarT_INPUTS, 300, 300);
double rho1 = rHEOS.rhomolar();
double analytic = MixtureDerivatives::d_dalpharddelta_dxj__constT_V_xi(rHEOS, i, xN_flag);
// Increment mole fraction
std::vector<long double> zp = z; /// Copy base composition
zp[i] += dz;
rHEOS.set_mole_fractions(zp);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v1 = rHEOS.dalphar_dDelta();
// Decrement mole fraction
std::vector<long double> zm = z; /// Copy base composition
zm[i] -= dz;
rHEOS.set_mole_fractions(zm);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v2 = rHEOS.dalphar_dDelta();
double numeric = (v1 - v2)/(2*dz);
double err = std::abs((numeric-analytic)/analytic);
CAPTURE(numeric);
CAPTURE(analytic);
CHECK(err < 1e-60);
}
// These derivatives depend on both the i and j indices
for (std::size_t j = 0; j < z.size(); ++j){
std::ostringstream ss1;
@@ -335,24 +488,63 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
{
double dz = 1e-6;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(PT_INPUTS, 101325, 300);
rHEOS.update(DmolarT_INPUTS, 300, 300);
double p = rHEOS.p();
CAPTURE(p);
double rho1 = rHEOS.rhomolar();
double analytic = MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rHEOS, i, j);
double analytic = MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rHEOS, i, j, xN_flag);
std::vector<long double> zp = z; /// Copy base composition
zp[j] += dz;
rHEOS.set_mole_fractions(zp);
rHEOS.update(PT_INPUTS, 101325, 300);
double v1 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i);
rHEOS.update(PT_INPUTS, p, 300);
double v1 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i, xN_flag);
std::vector<long double> zm = z; /// Copy base composition
zm[j] -= dz;
rHEOS.set_mole_fractions(zm);
rHEOS.update(PT_INPUTS, 101325, 300);
double v2 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i);
rHEOS.update(PT_INPUTS, p, 300);
double v2 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i, xN_flag);
double numeric = (v1 - v2)/(2*dz);
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-6);
CAPTURE(numeric);
CAPTURE(analytic);
CHECK(err < 1e-8);
}
std::ostringstream ss3;
ss2 << "d_ndalphardni_dxj__constT_V_xi, i=" << i << ", j=" << j;
SECTION(ss2.str(), "")
{
double dz = 1e-6;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(DmolarT_INPUTS, 300, 300);
double rho1 = rHEOS.rhomolar();
double analytic = MixtureDerivatives::d_ndalphardni_dxj__constT_V_xi(rHEOS, i, j, xN_flag);
// Increment mole fraction
std::vector<long double> zp = z; /// Copy base composition
zp[j] += dz;
rHEOS.set_mole_fractions(zp);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v1 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS, i, xN_flag);
// Decrement mole fraction
std::vector<long double> zm = z; /// Copy base composition
zm[j] -= dz;
rHEOS.set_mole_fractions(zm);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v2 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS, i, xN_flag);
double numeric = (v1 - v2)/(2*dz);
double err = std::abs((numeric-analytic)/analytic);
CAPTURE(numeric);
CAPTURE(analytic);
CHECK(err < 1e-60);
}
}
}

82
src/Backends/Helmholtz/MixtureDerivatives.h
View File

@@ -29,11 +29,16 @@ and not pollute the HelmholtzEOSMixtureBackend namespace
class MixtureDerivatives{
public:
static long double dalphar_dxi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
static long double d2alphar_dxi_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
static long double d2alphar_dxi_dDelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
static long double d2alphardxidxj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j);
/** \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.61
*
* The derivative term
@@ -45,16 +50,13 @@ class MixtureDerivatives{
static long double dpdT__constV_n(HelmholtzEOSMixtureBackend &HEOS);
static long double dpdrho__constT_n(HelmholtzEOSMixtureBackend &HEOS);
static long double dalphar_dxi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
static long double d2alphar_dxi_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
static long double d2alphar_dxi_dDelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
static long double d2alphardxidxj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/** \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
*
@@ -64,7 +66,7 @@ class MixtureDerivatives{
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double ndpdni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
static long double ndpdni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 monograph Eqn. 7.32
*
@@ -74,7 +76,7 @@ class MixtureDerivatives{
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double partial_molar_volume(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
static long double partial_molar_volume(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief Fugacity of the i-th component
*
@@ -83,25 +85,25 @@ class MixtureDerivatives{
* f_i(\delta, \tau, \bar x) = x_i\rho R T \exp\left(\frac{\partial n\alpha^r}{\partial n_i}\right)_{T,V,n_{j \neq i}}
* \f]
*/
static long double fugacity_i(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
static long double fugacity_i(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief Natural logarithm of the fugacity coefficient
*
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double ln_fugacity_coefficient(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
static long double ln_fugacity_coefficient(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \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, std::size_t i);
static long double dln_fugacity_coefficient_dT__constrho_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \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, std::size_t i);
static long double dln_fugacity_coefficient_drho__constT_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph Eqn. 7.29
*
@@ -112,7 +114,7 @@ class MixtureDerivatives{
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double dln_fugacity_coefficient_dT__constp_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
static long double dln_fugacity_coefficient_dT__constp_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief Table B4, Kunz, JCED, 2012 for the original term and the subsequent substitutions
*
@@ -128,10 +130,10 @@ class MixtureDerivatives{
* \f}
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double ndalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
static long double ndalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/// GERG Equation 7.42
static long double dnalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
static long double dnalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph Eqn. 7.30
*
@@ -141,7 +143,7 @@ class MixtureDerivatives{
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double dln_fugacity_coefficient_dp__constT_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
static long double dln_fugacity_coefficient_dp__constT_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph Equation 7.31
*
@@ -155,7 +157,7 @@ class MixtureDerivatives{
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double ndln_fugacity_coefficient_dnj__constT_p(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j);
static long double ndln_fugacity_coefficient_dnj__constT_p(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief Gernert Equation 3.115
*
@@ -165,7 +167,7 @@ class MixtureDerivatives{
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double dln_fugacity_coefficient_dxj__constT_p_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j);
static long double dln_fugacity_coefficient_dxj__constT_p_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief Gernert Equation 3.130
*
@@ -175,7 +177,7 @@ class MixtureDerivatives{
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double dpdxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j);
static long double dpdxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief Gernert Equation 3.117
*
@@ -185,27 +187,27 @@ class MixtureDerivatives{
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double d2nalphar_dxj_dni__constT_V(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, std::size_t i){ return MixtureDerivatives::d_ndalphardni_dxj__constT_V_xi(HEOS, i, j) + MixtureDerivatives::dalphar_dxj__constT_V_xi(HEOS, j);};
static long double d2nalphar_dxj_dni__constT_V(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, std::size_t i, x_N_dependency_flag xN_flag){ return MixtureDerivatives::d_ndalphardni_dxj__constT_V_xi(HEOS, i, j, xN_flag) + MixtureDerivatives::dalphar_dxj__constT_V_xi(HEOS, j, xN_flag);};
/// Gernert Equation 3.119
/// Catch test provided
static long double dalphar_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j);
static long double dalphar_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag);
/// Gernert Equation 3.118
/// Catch test provided
static long double d_ndalphardni_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j);
static long double d_ndalphardni_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/// Gernert Equation 3.134
/// Catch test provided
static long double d_dalpharddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j);
static long double d_dalpharddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag);
/// Gernert Equation 3.121
/// Catch test provided
static long double ddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j);
static long double ddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag);
/// Gernert Equation 3.122
/// Catch test provided
static long double dtau_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j);
static long double dtau_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph, equations 7.44 and 7.51
*
@@ -218,7 +220,7 @@ class MixtureDerivatives{
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double d2nalphar_dni_dT(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
static long double d2nalphar_dni_dT(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph Equation 7.51 and Table B4, Kunz, JCED, 2012
*
@@ -230,7 +232,7 @@ class MixtureDerivatives{
* \f}
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double d_ndalphardni_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
static long double d_ndalphardni_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph Equation 7.50 and Table B4, Kunz, JCED, 2012
*
@@ -242,7 +244,7 @@ class MixtureDerivatives{
* \f}
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double d_ndalphardni_dDelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
static long double d_ndalphardni_dDelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph equation 7.41
*
@@ -263,7 +265,7 @@ class MixtureDerivatives{
* \f}
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double nd2nalphardnidnj__constT_V(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j);
static long double nd2nalphardnidnj__constT_V(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph equation 7.48
*
@@ -273,7 +275,7 @@ class MixtureDerivatives{
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double nddeltadni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
static long double nddeltadni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph equation 7.49
*
@@ -283,7 +285,7 @@ class MixtureDerivatives{
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double ndtaudni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
static long double ndtaudni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph equation 7.52
*
@@ -297,7 +299,7 @@ class MixtureDerivatives{
* \f}
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double d_ndalphardni_dxj__constdelta_tau_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j);
static long double d_ndalphardni_dxj__constdelta_tau_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
}; /* class MixtureDerivatives */

301
src/Backends/Helmholtz/ReducingFunctions.cpp
View File

@@ -94,9 +94,9 @@ ReducingFunction *ReducingFunction::factory(const std::vector<CoolPropFluid*> &c
beta_T.resize(N, std::vector<long double>(N, 0));
gamma_T.resize(N, std::vector<long double>(N, 0));
for (unsigned int i = 0; i < N; ++i)
for (std::size_t i = 0; i < N; ++i)
{
for (unsigned int j = 0; j < N; ++j)
for (std::size_t j = 0; j < N; ++j)
{
if (i == j){ continue; }
@@ -144,96 +144,124 @@ ReducingFunction *ReducingFunction::factory(const std::vector<CoolPropFluid*> &c
return new GERG2008ReducingFunction(components,beta_v, gamma_v, beta_T, gamma_T);
}
long double ReducingFunction::d_ndTrdni_dxj__constxi(const std::vector<long double> &x, int i, int j)
long double ReducingFunction::d_ndTrdni_dxj__constxi(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
long double s = 0;
for (unsigned int k = 0; k < N; k++)
for (std::size_t k = 0; k < N; k++)
{
s += x[k]*d2Trdxidxj(x,j,k);
s += x[k]*d2Trdxidxj(x,j,k, xN_flag);
}
return d2Trdxidxj(x,i,j)-dTrdxi__constxj(x,j)-s;
return d2Trdxidxj(x,i,j, xN_flag)-dTrdxi__constxj(x,j, xN_flag)-s;
}
long double ReducingFunction::d_ndrhorbardni_dxj__constxi(const std::vector<long double> &x, int i, int j)
long double ReducingFunction::d_ndrhorbardni_dxj__constxi(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
long double s = 0;
for (unsigned int k = 0; k < N; k++)
for (std::size_t k = 0; k < N; k++)
{
s += x[k]*d2rhormolardxidxj(x,j,k);
s += x[k]*d2rhormolardxidxj(x,j,k, xN_flag);
}
return d2rhormolardxidxj(x,j,i)-drhormolardxi__constxj(x,j)-s;
return d2rhormolardxidxj(x,j,i, xN_flag)-drhormolardxi__constxj(x,j, xN_flag)-s;
}
long double ReducingFunction::ndrhorbardni__constnj(const std::vector<long double> &x, int i)
long double ReducingFunction::ndrhorbardni__constnj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
long double summer_term1 = 0;
for (unsigned int j = 0; j < N; j++)
{
summer_term1 += x[j]*drhormolardxi__constxj(x,j);
}
return drhormolardxi__constxj(x,i)-summer_term1;
if (xN_flag == XN_INDEPENDENT){
long double summer_term1 = 0;
for (std::size_t j = 0; j < N; j++)
{
summer_term1 += x[j]*drhormolardxi__constxj(x,j, xN_flag);
}
return drhormolardxi__constxj(x,i, xN_flag)-summer_term1;
}
else if (xN_flag == XN_DEPENDENT){
long double summer_term1 = 0;
for (std::size_t k = 0; k < N-1; ++k)
{
summer_term1 += x[k]*drhormolardxi__constxj(x, k, xN_flag);
}
return drhormolardxi__constxj(x, i, xN_flag)-summer_term1;
}
else{
throw ValueError(format("xN dependency flag invalid"));
}
}
long double ReducingFunction::ndTrdni__constnj(const std::vector<long double> &x, int i)
long double ReducingFunction::ndTrdni__constnj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
// GERG Equation 7.54
long double summer_term1 = 0;
for (unsigned int j = 0; j < N; j++)
{
summer_term1 += x[j]*dTrdxi__constxj(x,j);
}
return dTrdxi__constxj(x,i)-summer_term1;
if (xN_flag == XN_INDEPENDENT){
// GERG Equation 7.54
long double summer_term1 = 0;
for (std::size_t j = 0; j < N; j++)
{
summer_term1 += x[j]*dTrdxi__constxj(x,j, xN_flag);
}
return dTrdxi__constxj(x,i, xN_flag)-summer_term1;
}
else if (xN_flag == XN_DEPENDENT){
long double summer_term1 = 0;
for (std::size_t k = 0; k < N-1; ++k)
{
summer_term1 += x[k]*dTrdxi__constxj(x, k, xN_flag);
}
return dTrdxi__constxj(x, i, xN_flag)-summer_term1;
}
else{
throw ValueError(format("xN dependency flag invalid"));
}
}
long double GERG2008ReducingFunction::Tr(const std::vector<long double> &x)
{
long double Tr = 0;
for (unsigned int i = 0; i < N; i++)
{
double xi = x[i], Tci = pFluids[i]->pEOS->reduce.T;
Tr += xi*xi*Tci;
// The last term is only used for the pure component, as it is sum_{i=1}^{N-1}sum_{j=1}^{N}
if (i==N-1){ break; }
for (unsigned int j = i+1; j < N; j++)
{
Tr += c_Y_ij(i, j, beta_T, gamma_T, T_c)*f_Y_ij(x, i, j, beta_T);
}
}
return Tr;
return Yr(x, beta_T, gamma_T, T_c, Yc_T);
}
long double GERG2008ReducingFunction::dTrdxi__constxj(const std::vector<long double> &x, int i)
long double GERG2008ReducingFunction::dTrdxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
// See Table B9 from Kunz Wagner 2012 (GERG 2008)
long double xi = x[i];
long double dTr_dxi = 2*xi*pFluids[i]->pEOS->reduce.T;
for (int k = 0; k < i; k++)
{
dTr_dxi += c_Y_ji(k,i,beta_T,gamma_T,T_c)*dfYkidxi__constxk(x,k,i,beta_T);
}
for (unsigned int k = i+1; k < N; k++)
{
dTr_dxi += c_Y_ij(i,k,beta_T,gamma_T,T_c)*dfYikdxi__constxk(x,i,k,beta_T);
}
return dTr_dxi;
return dYrdxi__constxj(x, i, beta_T, gamma_T, T_c, Yc_T, xN_flag);
}
long double GERG2008ReducingFunction::d2Trdxi2__constxj(const std::vector<long double> &x, int i)
long double GERG2008ReducingFunction::rhormolar(const std::vector<long double> &x)
{
return 1/Yr(x, beta_v, gamma_v, v_c, Yc_v);
}
long double GERG2008ReducingFunction::drhormolardxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
return -pow(rhormolar(x),2)*dvrmolardxi__constxj(x, i, xN_flag);
}
long double GERG2008ReducingFunction::dvrmolardxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
return dYrdxi__constxj(x, i, beta_v, gamma_v, v_c, Yc_v, xN_flag);
}
long double GERG2008ReducingFunction::d2rhormolardxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
long double rhor = this->rhormolar(x);
long double dvrbardxi = this->dvrmolardxi__constxj(x,i, xN_flag);
return 2*pow(rhor,(std::size_t)3)*pow(dvrbardxi,(std::size_t)2)-pow(rhor,(std::size_t)2)*this->d2vrmolardxi2__constxj(x,i, xN_flag);
}
long double GERG2008ReducingFunction::d2rhormolardxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
double rhor = this->rhormolar(x);
double dvrbardxi = this->dvrmolardxi__constxj(x,i, xN_flag);
double dvrbardxj = this->dvrmolardxi__constxj(x,j, xN_flag);
return 2*pow(rhor,(std::size_t)3)*dvrbardxi*dvrbardxj-pow(rhor,(std::size_t)2)*this->d2vrmolardxidxj(x,i,j, xN_flag);
}
long double GERG2008ReducingFunction::d2Trdxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
// See Table B9 from Kunz Wagner 2012 (GERG 2008)
long double d2Tr_dxi2 = 2*pFluids[i]->pEOS->reduce.T;
for (int k = 0; k < i; k++)
for (std::size_t k = 0; k < i; k++)
{
d2Tr_dxi2 += c_Y_ij(k,i,beta_T,gamma_T,T_c)*d2fYkidxi2__constxk(x,k,i,beta_T);
}
for (unsigned int k = i+1; k < N; k++)
for (std::size_t k = i+1; k < N; k++)
{
d2Tr_dxi2 += c_Y_ij(i,k,beta_T,gamma_T,T_c)*d2fYikdxi2__constxk(x,i,k,beta_T);
}
return d2Tr_dxi2;
}
long double GERG2008ReducingFunction::d2Trdxidxj(const std::vector<long double> &x, int i, int j)
long double GERG2008ReducingFunction::d2Trdxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
if (i == j)
{
return d2Trdxi2__constxj(x,i);
return d2Trdxi2__constxj(x, i, xN_flag);
}
else
{
@@ -241,122 +269,161 @@ long double GERG2008ReducingFunction::d2Trdxidxj(const std::vector<long double>
return c_Y_ij(i, j, beta_T, gamma_T, T_c)*d2fYijdxidxj(x, i, j, beta_T);
}
}
long double GERG2008ReducingFunction::dvrmolardxi__constxj(const std::vector<long double> &x, int i)
{
// See Table B9 from Kunz Wagner 2012 (GERG 2008)
long double xi = x[i];
long double dvrbar_dxi = 2*xi/pFluids[i]->pEOS->reduce.rhomolar;
for (int k = 0; k < i; k++)
{
dvrbar_dxi += c_Y_ij(k, i, beta_v, gamma_v, v_c)*dfYkidxi__constxk(x, k, i, beta_v);
}
for (unsigned int k = i+1; k < N; k++)
{
dvrbar_dxi += c_Y_ij(i, k, beta_v, gamma_v, v_c)*dfYikdxi__constxk(x, i, k, beta_v);
}
return dvrbar_dxi;
}
long double GERG2008ReducingFunction::d2vrmolardxidxj(const std::vector<long double> &x, int i, int j)
long double GERG2008ReducingFunction::d2vrmolardxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
if (i == j)
{
return d2vrmolardxi2__constxj(x, i);
return d2vrmolardxi2__constxj(x, i, xN_flag);
}
else
{
return c_Y_ij(i, j, beta_v, gamma_v, v_c)*d2fYijdxidxj(x, i, j, beta_v);
}
}
long double GERG2008ReducingFunction::drhormolardxi__constxj(const std::vector<long double> &x, int i)
{
return -pow(rhormolar(x),2)*dvrmolardxi__constxj(x,i);
}
long double GERG2008ReducingFunction::d2vrmolardxi2__constxj(const std::vector<long double> &x, int i)
long double GERG2008ReducingFunction::d2vrmolardxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
// See Table B9 from Kunz Wagner 2012 (GERG 2008)
double d2vrbardxi2 = 2/pFluids[i]->pEOS->reduce.rhomolar;
for (int k = 0; k < i; k++)
for (std::size_t k = 0; k < i; k++)
{
d2vrbardxi2 += c_Y_ij(k, i, beta_v, gamma_v, v_c)*d2fYkidxi2__constxk(x, k, i, beta_v);
}
for (unsigned int k = i+1; k < N; k++)
for (std::size_t k = i+1; k < N; k++)
{
d2vrbardxi2 += c_Y_ij(i, k, beta_v, gamma_v, v_c)*d2fYikdxi2__constxk(x, i, k, beta_v);
}
return d2vrbardxi2;
}
long double GERG2008ReducingFunction::d2rhormolardxi2__constxj(const std::vector<long double> &x, int i)
long double GERG2008ReducingFunction::Yr(const std::vector<long double> &x, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector<long double> &Yc)
{
long double rhor = this->rhormolar(x);
long double dvrbardxi = this->dvrmolardxi__constxj(x,i);
return 2*pow(rhor,(int)3)*pow(dvrbardxi,(int)2)-pow(rhor,(int)2)*this->d2vrmolardxi2__constxj(x,i);
long double Yr = 0;
// if (xN_flag == XN_INDEPENDENT)
// {
for (std::size_t i = 0; i < N; i++)
{
double xi = x[i];
Yr += xi*xi*Yc[i];
// The last term is only used for the pure component, as it is sum_{i=1}^{N-1}sum_{j=1}^{N}
if (i==N-1){ break; }
for (std::size_t j = i+1; j < N; j++)
{
Yr += c_Y_ij(i, j, beta, gamma, Y_c_ij)*f_Y_ij(x, i, j, beta);
}
}
// }
// else if (xN_flag == XN_DEPENDENT){
//
// // A.43 from Gernert, 2014, supplemental information
// for (std::size_t k = 0; k < N-1; ++k)
// {
// double xk = x[k];
// Yr += xk*xk*Yc[k];
// }
// for (std::size_t k = 0; k < N-2; ++k)
// {
// for (std::size_t m = k+1; m < N-1; ++m)
// {
// Yr += c_Y_ij(k, m, beta, gamma, Y_c_ij)*f_Y_ij(x, k, m, beta);
// }
// }
// for (std::size_t k = 0; k < N-1; ++k)
// {
// Yr += c_Y_ij(k, N-1, beta, gamma, Y_c_ij)*f_Y_ij(x, k, N-1, beta);
// }
// double xN = x[N-1];
// Yr += xN*xN*Yc[N-1];
// }
// else{
// throw ValueError(format("xN dependency flag invalid"));
// }
return Yr;
}
long double GERG2008ReducingFunction::d2rhormolardxidxj(const std::vector<long double> &x, int i, int j)
long double GERG2008ReducingFunction::dYrdxi__constxj(const std::vector<long double> &x, std::size_t i, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector<long double> &Yc, x_N_dependency_flag xN_flag)
{
double rhor = this->rhormolar(x);
double dvrbardxi = this->dvrmolardxi__constxj(x,i);
double dvrbardxj = this->dvrmolardxi__constxj(x,j);
return 2*pow(rhor,(int)3)*dvrbardxi*dvrbardxj-pow(rhor,(int)2)*this->d2vrmolardxidxj(x,i,j);
if (xN_flag == XN_INDEPENDENT){
// See Table B9 from Kunz Wagner 2012 (GERG 2008)
long double xi = x[i];
long double dYr_dxi = 2*xi*Yc[i];
for (std::size_t k = 0; k < i; k++)
{
dYr_dxi += c_Y_ji(k,i,beta,gamma,Y_c_ij)*dfYkidxi__constxk(x,k,i,beta);
}
for (std::size_t k = i+1; k < N; k++)
{
dYr_dxi += c_Y_ij(i,k,beta,gamma,Y_c_ij)*dfYikdxi__constxk(x,i,k,beta);
}
return dYr_dxi;
}
else if (xN_flag == XN_DEPENDENT){
// Table S1 from Gernert, 2014, supplemental information
long double dYr_dxi = 2*x[i]*Yc[i] - 2*x[N-1]*Yc[N-1];
for (std::size_t k = 0; k < i; k++)
{
dYr_dxi += c_Y_ji(k, i, beta, gamma, Y_c_ij)*dfYkidxi__constxk(x,k,i,beta);
}
for (std::size_t k = i+1; k < N-1; k++)
{
dYr_dxi += c_Y_ji(i, k, beta, gamma, Y_c_ij)*dfYkidxi__constxk(x,i,k,beta);
}
double beta_Y_iN = beta[i][N-1], xN = x[N-1];
dYr_dxi += c_Y_ji(i, N-1, beta, gamma, Y_c_ij)*(xN*(x[i]+xN)/(pow(beta_Y_iN,2)*x[i]+xN)+(1-beta_Y_iN*beta_Y_iN)*x[i]*xN*xN/pow(beta_Y_iN*beta_Y_iN*x[i]+xN, 2));
for (std::size_t k = 0; k < N-1; ++k)
{
double beta_Y_kN = beta[k][N-1];
dYr_dxi += c_Y_ji(k, N-1, beta, gamma, Y_c_ij)*(-x[k]*(x[k]+xN)/(pow(beta_Y_kN,2)*x[k]+xN)+(1-beta_Y_kN*beta_Y_kN)*xN*x[k]*x[k]/pow(beta_Y_kN*beta_Y_kN*x[k]+xN, 2));
}
return dYr_dxi;
}
else{
throw ValueError(format("xN dependency flag invalid"));
}
}
long double GERG2008ReducingFunction::rhormolar(const std::vector<long double> &x)
{
double vrbar = 0;
for (unsigned int i = 0; i < N; i++)
{
double xi = x[i];
vrbar += xi*xi/pFluids[i]->pEOS->reduce.rhomolar;
if (i == N-1){ break; }
for (unsigned int j = i+1; j < N; j++)
{
vrbar += c_Y_ij(i, j, beta_v, gamma_v, v_c)*f_Y_ij(x, i, j, beta_v);
}
}
return 1/vrbar;
}
long double GERG2008ReducingFunction::dfYkidxi__constxk(const std::vector<long double> &x, int k, int i, std::vector< std::vector< long double> > &beta)
long double GERG2008ReducingFunction::dfYkidxi__constxk(const std::vector<long double> &x, std::size_t k, std::size_t i, const STLMatrix &beta)
{
double xk = x[k], xi = x[i], beta_Y = beta[k][i];
return xk*(xk+xi)/(beta_Y*beta_Y*xk+xi)+xk*xi/(beta_Y*beta_Y*xk+xi)*(1-(xk+xi)/(beta_Y*beta_Y*xk+xi));
}
long double GERG2008ReducingFunction::dfYikdxi__constxk(const std::vector<long double> &x, int i, int k, std::vector< std::vector< long double> > &beta)
long double GERG2008ReducingFunction::dfYikdxi__constxk(const std::vector<long double> &x, std::size_t i, std::size_t k, const STLMatrix &beta)
{
double xk = x[k], xi = x[i], beta_Y = beta[i][k];
return xk*(xi+xk)/(beta_Y*beta_Y*xi+xk)+xi*xk/(beta_Y*beta_Y*xi+xk)*(1-beta_Y*beta_Y*(xi+xk)/(beta_Y*beta_Y*xi+xk));
}
long double GERG2008ReducingFunction::c_Y_ij(int i, int j, std::vector< std::vector< long double> > &beta, std::vector< std::vector< long double> > &gamma, std::vector< std::vector< long double> > &Y_c)
long double GERG2008ReducingFunction::c_Y_ij(std::size_t i, std::size_t j, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c)
{
return 2*beta[i][j]*gamma[i][j]*Y_c[i][j];
}
long double GERG2008ReducingFunction::c_Y_ji(int j, int i, std::vector< std::vector< long double> > &beta, std::vector< std::vector< long double> > &gamma, std::vector< std::vector< long double> > &Y_c)
long double GERG2008ReducingFunction::c_Y_ji(std::size_t j, std::size_t i, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c)
{
return 2/beta[i][j]*gamma[i][j]*Y_c[i][j];
}
long double GERG2008ReducingFunction::f_Y_ij(const std::vector<long double> &x, int i, int j, std::vector< std::vector< long double> > &beta)
long double GERG2008ReducingFunction::f_Y_ij(const std::vector<long double> &x, std::size_t i, std::size_t j, const STLMatrix &beta)
{
double xi = x[i], xj = x[j], beta_Y = beta[i][j];
return xi*xj*(xi+xj)/(beta_Y*beta_Y*xi+xj);
}
long double GERG2008ReducingFunction::d2fYikdxi2__constxk(const std::vector<long double> &x, int i, int k, std::vector< std::vector< long double> > &beta)
long double GERG2008ReducingFunction::d2fYikdxi2__constxk(const std::vector<long double> &x, std::size_t i, std::size_t k, const STLMatrix &beta)
{
double xi = x[i], xk = x[k], beta_Y = beta[i][k];
return 1/(beta_Y*beta_Y*xi+xk)*(1-beta_Y*beta_Y*(xi+xk)/(beta_Y*beta_Y*xi+xk))*(2*xk-xi*xk*2*beta_Y*beta_Y/(beta_Y*beta_Y*xi+xk));
}
long double GERG2008ReducingFunction::d2fYkidxi2__constxk(const std::vector<long double> &x, int k, int i, std::vector< std::vector< long double> > &beta)
long double GERG2008ReducingFunction::d2fYkidxi2__constxk(const std::vector<long double> &x, std::size_t k, std::size_t i, const STLMatrix &beta)
{
double xi = x[i], xk = x[k], beta_Y = beta[k][i];
return 1/(beta_Y*beta_Y*xk+xi)*(1-(xk+xi)/(beta_Y*beta_Y*xk+xi))*(2*xk-xk*xi*2/(beta_Y*beta_Y*xk+xi));
}
long double GERG2008ReducingFunction::d2fYijdxidxj(const std::vector<long double> &x, int i, int j, std::vector< std::vector< long double> > &beta)
long double GERG2008ReducingFunction::d2fYijdxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, const STLMatrix &beta)
{
double xi = x[i], xj = x[j], beta_Y = beta[i][j], beta_Y2 = beta_Y*beta_Y;
return (xi+xj)/(beta_Y2*xi+xj) + xj/(beta_Y2*xi+xj)*(1-(xi+xj)/(beta_Y2*xi+xj))
+xi/(beta_Y2*xi+xj)*(1-beta_Y2*(xi+xj)/(beta_Y2*xi+xj))
-xi*xj/pow(beta_Y2*xi+xj,(int)2)*(1+beta_Y2-2*beta_Y2*(xi+xj)/(beta_Y2*xi+xj));
-xi*xj/pow(beta_Y2*xi+xj,(std::size_t)2)*(1+beta_Y2-2*beta_Y2*(xi+xj)/(beta_Y2*xi+xj));
}

79
src/Backends/Helmholtz/ReducingFunctions.h
View File

@@ -8,6 +8,9 @@ namespace CoolProp{
typedef std::vector<std::vector<long double> > STLMatrix;
enum x_N_dependency_flag{XN_INDEPENDENT, ///< x_N is an independent variable, and not calculated by \f$ x_N = 1-\sum_i x_i\f$
XN_DEPENDENT}; ///< x_N is an dependent variable, calculated by \f$ x_N = 1-\sum_i x_i\f$
/// A container for the mixing parameters for CoolProp mixtures
/**
@@ -60,7 +63,7 @@ reducing parameters \f$ \rho_r \f$ and \f$ T_r \f$
class ReducingFunction
{
protected:
unsigned int N;
std::size_t N;
public:
ReducingFunction(){};
virtual ~ReducingFunction(){};
@@ -71,32 +74,32 @@ public:
/// The reduced temperature
virtual long double Tr(const std::vector<long double> &x) = 0;
/// The derivative of reduced temperature with respect to component i mole fraction
virtual long double dTrdxi__constxj(const std::vector<long double> &x, int i) = 0;
virtual long double dTrdxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag) = 0;
/// The molar reducing density
virtual long double rhormolar(const std::vector<long double> &x) = 0;
///Derivative of the molar reducing density with respect to component i mole fraction
virtual long double drhormolardxi__constxj(const std::vector<long double> &x, int i) = 0;
virtual long double drhormolardxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag) = 0;
virtual long double d2rhormolardxi2__constxj(const std::vector<long double> &x, int i) = 0;
virtual long double d2rhormolardxidxj(const std::vector<long double> &x, int i, int j) = 0;
virtual long double d2Trdxi2__constxj(const std::vector<long double> &x, int i) = 0;
virtual long double d2Trdxidxj(const std::vector<long double> &x, int i, int j) = 0;
virtual long double d2rhormolardxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag) = 0;
virtual long double d2rhormolardxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) = 0;
virtual long double d2Trdxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag) = 0;
virtual long double d2Trdxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) = 0;
/*! GERG 2004 Monograph equation 7.56:
\f[
\left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial T_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2T_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial T_r}{\partial x_j}\right)_{x_i}-\sum_{k=1}^Nx_k\left(\frac{\partial^2T_r}{\partial x_j \partial x_k}\right)
\f]
*/
long double d_ndTrdni_dxj__constxi(const std::vector<long double> &x, int i, int j);
long double d_ndTrdni_dxj__constxi(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/*! GERG 2004 Monograph equation 7.55:
\f[
\left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \rho_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2\rho_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_i}-\sum_{k=1}^Nx_k\left(\frac{\partial^2\rho_r}{\partial x_j \partial x_k}\right)
\f]
*/
long double d_ndrhorbardni_dxj__constxi(const std::vector<long double> &x, int i, int j);
long double d_ndrhorbardni_dxj__constxi(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
long double ndrhorbardni__constnj(const std::vector<long double> &x, int i);
long double ndTrdni__constnj(const std::vector<long double> &x, int i);
long double ndrhorbardni__constnj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag);
long double ndTrdni__constnj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag);
};
/*!
@@ -114,7 +117,9 @@ protected:
STLMatrix gamma_v; ///< \f$ \gamma_{v,ij} \f$ from GERG-2008
STLMatrix beta_T; ///< \f$ \beta_{T,ij} \f$ from GERG-2008
STLMatrix gamma_T; ///< \f$ \gamma_{T,ij} \f$ from GERG-2008
std::vector<CoolPropFluid *> pFluids; ///< List of pointers to fluids
std::vector<long double> Yc_T; ///< Vector of critical temperatures for all components
std::vector<long double> Yc_v; ///< Vector of critical molar volumes for all components
std::vector<CoolPropFluid *> pFluids; ///< List of postd::size_ters to fluids
public:
GERG2008ReducingFunction(std::vector<CoolPropFluid *> pFluids, STLMatrix beta_v, STLMatrix gamma_v, STLMatrix beta_T, STLMatrix gamma_T)
@@ -127,13 +132,17 @@ public:
this->N = pFluids.size();
T_c.resize(N,std::vector<long double>(N,0));
v_c.resize(N,std::vector<long double>(N,0));
for (unsigned int i = 0; i < N; ++i)
Yc_T.resize(N);
Yc_v.resize(N);
for (std::size_t i = 0; i < N; ++i)
{
for (unsigned int j = 0; j < N; j++)
for (std::size_t j = 0; j < N; j++)
{
T_c[i][j] = sqrt(pFluids[i]->pEOS->reduce.T*pFluids[j]->pEOS->reduce.T);
v_c[i][j] = 1.0/8.0*pow(pow(pFluids[i]->pEOS->reduce.rhomolar, -1.0/3.0)+pow(pFluids[j]->pEOS->reduce.rhomolar, -1.0/3.0),(int)3);
v_c[i][j] = 1.0/8.0*pow(pow(pFluids[i]->pEOS->reduce.rhomolar, -1.0/3.0)+pow(pFluids[j]->pEOS->reduce.rhomolar, -1.0/3.0),(std::size_t)3);
}
Yc_T[i] = pFluids[i]->pEOS->reduce.T;
Yc_v[i] = 1/pFluids[i]->pEOS->reduce.rhomolar;
}
};
@@ -142,29 +151,31 @@ public:
/// The reduced temperature
long double Tr(const std::vector<long double> &x);
/// The derivative of reduced temperature with respect to component i mole fraction
long double dTrdxi__constxj(const std::vector<long double> &x, int i);
long double dTrdxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag);
/// The molar reducing density
long double rhormolar(const std::vector<long double> &x);
///Derivative of the molar reducing density with respect to component i mole fraction
long double drhormolardxi__constxj(const std::vector<long double> &x, int i);
long double dvrmolardxi__constxj(const std::vector<long double> &x, int i);
long double drhormolardxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag);
long double dvrmolardxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_fla);
long double d2vrmolardxi2__constxj(const std::vector<long double> &x, int i);
long double d2rhormolardxi2__constxj(const std::vector<long double> &x, int i);
long double d2vrmolardxidxj(const std::vector<long double> &x, int i, int j);
long double d2rhormolardxidxj(const std::vector<long double> &x, int i, int j);
long double d2Trdxi2__constxj(const std::vector<long double> &x, int i);
long double d2Trdxidxj(const std::vector<long double> &x, int i, int j);
long double d2vrmolardxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag);
long double d2rhormolardxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag);
long double d2vrmolardxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
long double d2rhormolardxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
long double d2Trdxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag);
long double d2Trdxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
long double c_Y_ij(int i, int j, std::vector< std::vector< long double> > &beta, std::vector< std::vector< long double> > &gamma, std::vector< std::vector< long double> > &Y_c);
long double c_Y_ji(int j, int i, std::vector< std::vector< long double> > &beta, std::vector< std::vector< long double> > &gamma, std::vector< std::vector< long double> > &Y_c);
long double f_Y_ij(const std::vector<long double> &x, int i, int j, std::vector< std::vector< long double> > &beta);
long double Yr(const std::vector<long double> &x, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector<long double> &Yc);
long double dYrdxi__constxj(const std::vector<long double> &x, std::size_t i, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector<long double> &Yc, x_N_dependency_flag xN_flag);
long double c_Y_ij(std::size_t i, std::size_t j, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c);
long double c_Y_ji(std::size_t j, std::size_t i, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c);
long double f_Y_ij(const std::vector<long double> &x, std::size_t i, std::size_t j, const STLMatrix &beta);
long double dfYkidxi__constxk(const std::vector<long double> &x, int k, int i,std::vector< std::vector< long double> > &beta);
long double dfYikdxi__constxk(const std::vector<long double> &x, int i, int k, std::vector< std::vector< long double> > &beta);
long double d2fYkidxi2__constxk(const std::vector<long double> &x, int k, int i, std::vector< std::vector< long double> > &beta);
long double d2fYikdxi2__constxk(const std::vector<long double> &x, int i, int k, std::vector< std::vector< long double> > &beta);
long double d2fYijdxidxj(const std::vector<long double> &x, int i, int k, std::vector< std::vector< long double> > &beta);
long double dfYkidxi__constxk(const std::vector<long double> &x, std::size_t k, std::size_t i, const STLMatrix &beta);
long double dfYikdxi__constxk(const std::vector<long double> &x, std::size_t i, std::size_t k, const STLMatrix &beta);
long double d2fYkidxi2__constxk(const std::vector<long double> &x, std::size_t k, std::size_t i, const STLMatrix &beta);
long double d2fYikdxi2__constxk(const std::vector<long double> &x, std::size_t i, std::size_t k, const STLMatrix &beta);
long double d2fYijdxidxj(const std::vector<long double> &x, std::size_t i, std::size_t k, const STLMatrix &beta);
};
/*! From Lemmon, JPCRD, 2000 for the properties of Dry Air, and also from Lemmon, JPCRD, 2004 for the properties of R404A, R410A, etc.
@@ -194,7 +205,7 @@ protected:
LemmonAirHFCReducingFunction(const LemmonAirHFCReducingFunction &);
public:
/// Set the coefficients based on reducing parameters loaded from JSON
static void convert_to_GERG(const std::vector<CoolPropFluid*> &pFluids, int i, int j, Dictionary d, long double &beta_T, long double &beta_v, long double &gamma_T, long double &gamma_v)
static void convert_to_GERG(const std::vector<CoolPropFluid*> &pFluids, std::size_t i, std::size_t j, Dictionary d, long double &beta_T, long double &beta_v, long double &gamma_T, long double &gamma_v)
{
long double xi_ij = d.get_number("xi");
long double zeta_ij = d.get_number("zeta");
@@ -204,7 +215,7 @@ public:
long double v_i = 1/pFluids[i]->pEOS->reduce.rhomolar;
long double v_j = 1/pFluids[j]->pEOS->reduce.rhomolar;
long double one_third = 1.0/3.0;
gamma_v = (v_i + v_j + zeta_ij)/(0.25*pow(pow(v_i, one_third)+pow(v_j, one_third),(int)3));
gamma_v = (v_i + v_j + zeta_ij)/(0.25*pow(pow(v_i, one_third)+pow(v_j, one_third),(std::size_t)3));
};
};

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

@@ -737,18 +737,19 @@ void SaturationSolvers::successive_substitution(HelmholtzEOSMixtureBackend &HEOS
f = 0;
df = 0;
x_N_dependency_flag xN_flag = XN_INDEPENDENT;
for (std::size_t i = 0; i < N; ++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);
ln_phi_liq[i] = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatL.get()), i, xN_flag);
ln_phi_vap[i] = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatV.get()), i, xN_flag);
if (options.sstype == imposed_p){
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);
deriv_liq = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(*(HEOS.SatL.get()), i, xN_flag);
deriv_vap = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(*(HEOS.SatV.get()), i, xN_flag);
}
else if (options.sstype == imposed_T){
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);
deriv_liq = MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(*(HEOS.SatL.get()), i, xN_flag);
deriv_vap = MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(*(HEOS.SatV.get()), i, xN_flag);
}
else {throw ValueError();}
@@ -960,25 +961,26 @@ void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureB
// -------
// Build the residual vector and the Jacobian matrix
x_N_dependency_flag xN_flag = XN_INDEPENDENT;
// For the residuals F_i
for (unsigned int i = 0; i < 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);
long double ln_phi_liq = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatL.get()), i, xN_flag);
long double phi_iT_liq = MixtureDerivatives::dln_fugacity_coefficient_dT__constrho_n(*(HEOS.SatL.get()), i, xN_flag);
dlnphi_drho_liq[i] = MixtureDerivatives::dln_fugacity_coefficient_drho__constT_n(*(HEOS.SatL.get()), i, xN_flag);
for (unsigned int j = 0; j < N; ++j)
{
// I think this is wrong.
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
phi_ij_liq[j] = MixtureDerivatives::ndln_fugacity_coefficient_dnj__constT_p(rSatL, i, j, xN_flag) + (MixtureDerivatives::partial_molar_volume(rSatL, i, xN_flag)/(SatL->gas_constant()*T)-1/p)*MixtureDerivatives::ndpdni__constT_V_nj(rSatL, i, xN_flag); // 7.126 from GERG monograph
}
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);
long double ln_phi_vap = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatV.get()), i, xN_flag);
long double phi_iT_vap = MixtureDerivatives::dln_fugacity_coefficient_dT__constrho_n(*(HEOS.SatV.get()), i, xN_flag);
dlnphi_drho_vap[i] = MixtureDerivatives::dln_fugacity_coefficient_drho__constT_n(*(HEOS.SatV.get()), i, xN_flag);
for (unsigned int j = 0; j < N; ++j)
{
// I think this is wrong.
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
phi_ij_vap[j] = MixtureDerivatives::ndln_fugacity_coefficient_dnj__constT_p(rSatV, i,j, xN_flag) + (MixtureDerivatives::partial_molar_volume(rSatV, i, xN_flag)/(SatV->gas_constant()*T)-1/p)*MixtureDerivatives::ndpdni__constT_V_nj(rSatV, i, xN_flag); ; // 7.126 from GERG monograph
}
r[i] = log(K[i]) + ln_phi_vap - ln_phi_liq;
@@ -1082,7 +1084,7 @@ void SaturationSolvers::newton_raphson_saturation::check_Jacobian()
rSatL.update_TP_guessrho(T, p, rhomolar_liq); rhomolar_liq = rSatL.rhomolar();
build_arrays(); r1 = r;
x = x0; x[0] -= dx;// x[1] += dx;
x = x0; x[0] -= dx; x[1] += dx;
rSatL.set_mole_fractions(x);
rSatL.update_TP_guessrho(T, p, rhomolar_liq); rhomolar_liq = rSatL.rhomolar();
build_arrays(); r2 = r;
@@ -1094,41 +1096,6 @@ void SaturationSolvers::newton_raphson_saturation::check_Jacobian()
std::cout << "analytic: " << vec_to_string(get_col(J0, 0), "%0.11Lg") << std::endl;
int rrr = 0;
// Jnum.resize(N, std::vector<long double>(N, 0));
//
// for (std::size_t j = 0; j < N; ++j)
// {
// std::vector<long double> KK = K;
// KK[j] += 1e-6;
// build_arrays(HEOS,beta0,T0,rhomolar_liq0,rhomolar_vap0,z,KK);
// std::vector<long double> r1 = r;
// for (std::size_t i = 0; i < N+2; ++i)
// {
// Jnum[i][j] = (r1[i]-r0[i])/(log(KK[j])-log(K[j]));
// }
// std::cout << vec_to_string(get_col(Jnum,j),"%12.11f") << std::endl;
// std::cout << vec_to_string(get_col(J,j),"%12.11f") << std::endl;
// }
//
// build_arrays(HEOS,beta0,T0+1e-6,rhomolar_liq0,rhomolar_vap0,z,K);
// std::vector<long double> r1 = r, JN = r;
// for (std::size_t i = 0; i < r.size(); ++i)
// {
// Jnum[i][N] = (r1[i]-r0[i])/(log(T0+1e-6)-log(T0));
// }
// std::cout << vec_to_string(get_col(Jnum,N),"%12.11f") << std::endl;
// std::cout << vec_to_string(get_col(J,N),"%12.11f") << std::endl;
//
// // Build the Jacobian and residual vectors for given inputs of K_i,T,p
// build_arrays(HEOS,beta0,T0,rhomolar_liq0+1e-3,rhomolar_vap0,z,K);
// std::vector<long double> r2 = r, JNp1 = r;
// for (std::size_t i = 0; i < r.size(); ++i)
// {
// Jnum[i][N+1] = (r2[i]-r0[i])/(log(rhomolar_liq0+1e-3)-log(rhomolar_liq0));
// }
// std::cout << vec_to_string(get_col(Jnum, N+1),"%12.11f") << std::endl;
// std::cout << vec_to_string(get_col(J,N+1),"%12.11f") << std::endl;
}
void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBackend &HEOS, const std::vector<long double> &z, std::vector<long double> &z_incipient, newton_raphson_saturation_options &IO)
{
@@ -1147,17 +1114,17 @@ void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBacke
// Hold a pointer to the backend
this->HEOS = &HEOS;
check_Jacobian();
do
{
// Build the Jacobian and residual vectors
build_arrays();
check_Jacobian();
// Solve for the step; v is the step with the contents
// [delta(x_0), delta(x_1), ..., delta(x_{N-2}), delta(spec)]
std::cout << vec_to_string(negative_r, "%0.12Lg") << std::endl;
std::cout << vec_to_string(J, "%0.12Lg") << std::endl;
std::cout << "-r: " << vec_to_string(negative_r, "%0.12Lg") << std::endl;
std::cout << "J: " << vec_to_string(J, "%0.12Lg") << std::endl;
std::vector<long double> v = linsolve(J, negative_r);
max_rel_change = max_abs_value(v);
@@ -1209,30 +1176,31 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
// -------
// Build the residual vector and the Jacobian matrix
x_N_dependency_flag xN_flag = XN_DEPENDENT;
// For the residuals F_i (equality of fugacities)
for (std::size_t i = 0; i < N; ++i)
{
// Equate the liquid and vapor fugacities
long double ln_f_liq = log(MixtureDerivatives::fugacity_i(rSatL, i));
long double ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i));
long double ln_f_liq = log(MixtureDerivatives::fugacity_i(rSatL, i, xN_flag));
long double ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i, xN_flag));
r[i] = ln_f_liq - ln_f_vap;
if (bubble_point){
throw NotImplementedError();
}
else{
else{ // its a dewpoint calculation
for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2
if (i == N-1){
J[i][j] = (-1/x[i] + MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rSatL,i,j));
J[N-1][j] = (-1/x[N-1] + MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rSatL,N-1,j, xN_flag));
}
else if (i != j){
J[i][j] = MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rSatL,i,j);
J[i][j] = MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rSatL,i,j, xN_flag);
}
else{
J[i][i] = (1/x[i] + MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rSatL,i,i));
J[i][i] = (1/x[i] + MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rSatL,i,i, xN_flag));
}
}
J[i][N-1] = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatL, i) - MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatV, i);
J[i][N-1] = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatV, i, xN_flag);
}
}