diff --git a/Web/coolprop/snippets/mixture_derivative_table.cxx b/Web/coolprop/snippets/mixture_derivative_table.cxx new file mode 100644 index 00000000..38759cba --- /dev/null +++ b/Web/coolprop/snippets/mixture_derivative_table.cxx @@ -0,0 +1,92 @@ +#include "CoolProp.h" +#include "MixtureDerivatives.h" +#include +using namespace CoolProp; +int main() +{ + // Ethane/Propane mixture, 25/75 molar + std::vector components(2,"Ethane"); components[1] = "Propane"; + std::vector z(2,0.25); z[1] = 0.75; + + shared_ptr HEOS(new HelmholtzEOSMixtureBackend(components)); + HelmholtzEOSMixtureBackend &rHEOS = *(HEOS.get()); + HEOS->set_mole_fractions(z); + HEOS->specify_phase(iphase_gas); // So that we don't do a + HEOS->update(DmolarT_INPUTS, 300, 300); + + std::vector terms; + terms.push_back("p"); + terms.push_back("p2(deriv)"); + terms.push_back("rhor"); + terms.push_back("Tr"); + terms.push_back("dalphar_dDelta"); + terms.push_back("dTr_dxi"); + terms.push_back("drhor_dxi"); + terms.push_back("ndpdV__constT_n"); + terms.push_back("dpdxj__constT_V_xi"); + + terms.push_back("ln_fugacity_coefficient"); + terms.push_back("ndpdni__constT_V_nj"); + terms.push_back("dln_fugacity_coefficient_dxj__constT_p_xi"); + terms.push_back("d2nalphar_dxj_dni__constT_V"); + + terms.push_back("d2alphar_dxi_dDelta"); + + for (std::vector::iterator it = terms.begin(); it != terms.end(); ++it) + { + if (!it->compare("p")){ + printf("p: %0.16Lg\n", HEOS->p()); + } + else if (!it->compare("p2(deriv)")){ + printf("p calculated by rho*R*T*(1+deltadar_dDelta): %0.16Lg\n", HEOS->rhomolar()*HEOS->gas_constant()*HEOS->T()*(1+HEOS->delta()*HEOS->dalphar_dDelta())); + } + else if (!it->compare("dalphar_dDelta")){ + printf("dalphar_dDelta: %0.16Lg\n", HEOS->dalphar_dDelta()); + } + else if (!it->compare("rhor")){ + printf("rhor: %0.16Lg\n", HEOS->get_reducing_state().rhomolar); + } + else if (!it->compare("Tr")){ + printf("Tr: %0.16Lg\n", HEOS->get_reducing_state().T); + } + else if (!it->compare("dTr_dxi")){ + printf("dTr_dxi: %0.16Lg\n", HEOS->Reducing.p->dTrdxi__constxj(rHEOS.get_mole_fractions(), 0, XN_DEPENDENT)); + } + else if (!it->compare("drhor_dxi")){ + printf("drhor_dxi: %0.16Lg\n", HEOS->Reducing.p->drhormolardxi__constxj(rHEOS.get_mole_fractions(), 0, XN_DEPENDENT)); + } + else if(!it->compare("ndpdV__constT_n")){ + printf("ndpdV__constT_n: %0.16Lg\n", MixtureDerivatives::ndpdV__constT_n(rHEOS)); + } + else if(!it->compare("ln_fugacity_coefficient")){ + printf("ln_fugacity_coefficient(0): %0.16Lg\n", MixtureDerivatives::ln_fugacity_coefficient(rHEOS, 0, XN_DEPENDENT)); + printf("ln_fugacity_coefficient(1): %0.16Lg\n", MixtureDerivatives::ln_fugacity_coefficient(rHEOS, 1, XN_DEPENDENT)); + } + else if(!it->compare("dln_fugacity_coefficient_dxj__constT_p_xi")){ + printf("dln_fugacity_coefficient_dxj__constT_p_xi(0,0): %0.16Lg\n", MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rHEOS, 0, 0, XN_DEPENDENT)); + //printf("dln_fugacity_coefficient_dxj__constT_p_xi(0,1): %0.16Lg\n", MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rHEOS, 0, 1, XN_DEPENDENT)); + printf("dln_fugacity_coefficient_dxj__constT_p_xi(1,0): %0.16Lg\n", MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rHEOS, 1, 0, XN_DEPENDENT)); + //printf("dln_fugacity_coefficient_dxj__constT_p_xi(1,1): %0.16Lg\n", MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rHEOS, 1, 1, XN_DEPENDENT)); + } + else if(!it->compare("d2nalphar_dxj_dni__constT_V")){ + printf("d2nalphar_dxj_dni__constT_V(0,0): %0.16Lg\n", MixtureDerivatives::d2nalphar_dxj_dni__constT_V(rHEOS, 0, 0, XN_DEPENDENT)); + //printf("d2nalphar_dxj_dni__constT_V(0,1): %0.16Lg\n", MixtureDerivatives::d2nalphar_dxj_dni__constT_V(rHEOS, 0, 1, XN_DEPENDENT)); + printf("d2nalphar_dxj_dni__constT_V(1,0): %0.16Lg\n", MixtureDerivatives::d2nalphar_dxj_dni__constT_V(rHEOS, 1, 0, XN_DEPENDENT)); + //printf("d2nalphar_dxj_dni__constT_V(1,1): %0.16Lg\n", MixtureDerivatives::d2nalphar_dxj_dni__constT_V(rHEOS, 1, 1, XN_DEPENDENT)); + } + else if(!it->compare("d2alphar_dxi_dDelta")){ + printf("d2alphar_dxi_dDelta(0): %0.16Lg\n", MixtureDerivatives::d2alphar_dxi_dDelta(rHEOS, 0, XN_DEPENDENT)); + //printf("d2alphar_dxi_dDelta(1): %0.16Lg\n", MixtureDerivatives::d2alphar_dxi_dDelta(rHEOS, 1, XN_DEPENDENT)); + } + else if(!it->compare("ndpdni__constT_V_nj")){ + printf("ndpdni__constT_V_nj(0): %0.16Lg\n", MixtureDerivatives::ndpdni__constT_V_nj(rHEOS, 0, XN_DEPENDENT)); + //printf("ndpdni__constT_V_nj(1): %0.16Lg\n", MixtureDerivatives::ndpdni__constT_V_nj(rHEOS, 1, XN_DEPENDENT)); + } + else if(!it->compare("dpdxj__constT_V_xi")){ + printf("dpdxj__constT_V_xi(0): %0.16Lg\n", MixtureDerivatives::dpdxj__constT_V_xi(rHEOS, 0, XN_DEPENDENT)); + //printf("dpdxj__constT_V_xi(1): %0.16Lg\n", MixtureDerivatives::dpdxj__constT_V_xi(rHEOS, 1, XN_DEPENDENT)); + } + } + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/Web/coolprop/snippets/mixture_derivative_table.cxx.output b/Web/coolprop/snippets/mixture_derivative_table.cxx.output new file mode 100644 index 00000000..2191462e --- /dev/null +++ b/Web/coolprop/snippets/mixture_derivative_table.cxx.output @@ -0,0 +1,17 @@ +p: 675616.5260238834 +p calculated by rho*R*T*(1+deltadar_dDelta): 675616.5260238834 +rhor: 5375.073740431392 +Tr: 354.7802371525603 +dalphar_dDelta: -1.740349537597104 +dTr_dxi: -60.84527088778825 +drhor_dxi: 1587.916430678275 +ndpdV__constT_n: -181257322.2040499 +dpdxj__constT_V_xi(0): 48181.54565089135 +ln_fugacity_coefficient(0): -0.04441530290159221 +ln_fugacity_coefficient(1): -1.#IND +ndpdni__constT_V_nj(0): 640327.2332516684 +dln_fugacity_coefficient_dxj__constT_p_xi(0,0): -1.#IND +dln_fugacity_coefficient_dxj__constT_p_xi(1,0): -1.#IND +d2nalphar_dxj_dni__constT_V(0,0): -1.#IND +d2nalphar_dxj_dni__constT_V(1,0): -1.#IND +d2alphar_dxi_dDelta(0): 0.001221607211340736 diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 1b81cb10..b0055bad 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -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 & mole_fractions) diff --git a/src/Backends/Helmholtz/MixtureDerivatives.cpp b/src/Backends/Helmholtz/MixtureDerivatives.cpp index ceefc2dd..a8530a10 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.cpp +++ b/src/Backends/Helmholtz/MixtureDerivatives.cpp @@ -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 &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 &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 &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 names(2); - names[0] = "Methane"; names[1] = "Ethane"; + names[0] = "Ethane"; names[1] = "Propane"; std::vector z(2); - z[0] = 0.5; z[1] = 1-z[0]; + z[0] = 0.25; z[1] = 1-z[0]; shared_ptr 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 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 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 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 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 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 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 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 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); + } + } } diff --git a/src/Backends/Helmholtz/MixtureDerivatives.h b/src/Backends/Helmholtz/MixtureDerivatives.h index 11d3cdb0..a8c36baf 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.h +++ b/src/Backends/Helmholtz/MixtureDerivatives.h @@ -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 */ diff --git a/src/Backends/Helmholtz/ReducingFunctions.cpp b/src/Backends/Helmholtz/ReducingFunctions.cpp index 9cc9d7bb..1546afe6 100644 --- a/src/Backends/Helmholtz/ReducingFunctions.cpp +++ b/src/Backends/Helmholtz/ReducingFunctions.cpp @@ -94,9 +94,9 @@ ReducingFunction *ReducingFunction::factory(const std::vector &c beta_T.resize(N, std::vector(N, 0)); gamma_T.resize(N, std::vector(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 &c return new GERG2008ReducingFunction(components,beta_v, gamma_v, beta_T, gamma_T); } -long double ReducingFunction::d_ndTrdni_dxj__constxi(const std::vector &x, int i, int j) +long double ReducingFunction::d_ndTrdni_dxj__constxi(const std::vector &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 &x, int i, int j) +long double ReducingFunction::d_ndrhorbardni_dxj__constxi(const std::vector &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 &x, int i) +long double ReducingFunction::ndrhorbardni__constnj(const std::vector &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 &x, int i) +long double ReducingFunction::ndTrdni__constnj(const std::vector &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 &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 &x, int i) +long double GERG2008ReducingFunction::dTrdxi__constxj(const std::vector &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 &x, int i) +long double GERG2008ReducingFunction::rhormolar(const std::vector &x) +{ + return 1/Yr(x, beta_v, gamma_v, v_c, Yc_v); +} +long double GERG2008ReducingFunction::drhormolardxi__constxj(const std::vector &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 &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 &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 &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 &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 &x, int i, int j) +long double GERG2008ReducingFunction::d2Trdxidxj(const std::vector &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 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 &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 &x, int i, int j) +long double GERG2008ReducingFunction::d2vrmolardxidxj(const std::vector &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 &x, int i) -{ - return -pow(rhormolar(x),2)*dvrmolardxi__constxj(x,i); -} -long double GERG2008ReducingFunction::d2vrmolardxi2__constxj(const std::vector &x, int i) + +long double GERG2008ReducingFunction::d2vrmolardxi2__constxj(const std::vector &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 &x, int i) + +long double GERG2008ReducingFunction::Yr(const std::vector &x, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector &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 &x, int i, int j) +long double GERG2008ReducingFunction::dYrdxi__constxj(const std::vector &x, std::size_t i, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector &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 &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 &x, int k, int i, std::vector< std::vector< long double> > &beta) +long double GERG2008ReducingFunction::dfYkidxi__constxk(const std::vector &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 &x, int i, int k, std::vector< std::vector< long double> > &beta) +long double GERG2008ReducingFunction::dfYikdxi__constxk(const std::vector &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 &x, int i, int j, std::vector< std::vector< long double> > &beta) +long double GERG2008ReducingFunction::f_Y_ij(const std::vector &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 &x, int i, int k, std::vector< std::vector< long double> > &beta) +long double GERG2008ReducingFunction::d2fYikdxi2__constxk(const std::vector &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 &x, int k, int i, std::vector< std::vector< long double> > &beta) +long double GERG2008ReducingFunction::d2fYkidxi2__constxk(const std::vector &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 &x, int i, int j, std::vector< std::vector< long double> > &beta) +long double GERG2008ReducingFunction::d2fYijdxidxj(const std::vector &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)); } diff --git a/src/Backends/Helmholtz/ReducingFunctions.h b/src/Backends/Helmholtz/ReducingFunctions.h index 4ccc74ed..b8fbb18f 100644 --- a/src/Backends/Helmholtz/ReducingFunctions.h +++ b/src/Backends/Helmholtz/ReducingFunctions.h @@ -8,6 +8,9 @@ namespace CoolProp{ typedef std::vector > 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 &x) = 0; /// The derivative of reduced temperature with respect to component i mole fraction - virtual long double dTrdxi__constxj(const std::vector &x, int i) = 0; + virtual long double dTrdxi__constxj(const std::vector &x, std::size_t i, x_N_dependency_flag xN_flag) = 0; /// The molar reducing density virtual long double rhormolar(const std::vector &x) = 0; ///Derivative of the molar reducing density with respect to component i mole fraction - virtual long double drhormolardxi__constxj(const std::vector &x, int i) = 0; + virtual long double drhormolardxi__constxj(const std::vector &x, std::size_t i, x_N_dependency_flag xN_flag) = 0; - virtual long double d2rhormolardxi2__constxj(const std::vector &x, int i) = 0; - virtual long double d2rhormolardxidxj(const std::vector &x, int i, int j) = 0; - virtual long double d2Trdxi2__constxj(const std::vector &x, int i) = 0; - virtual long double d2Trdxidxj(const std::vector &x, int i, int j) = 0; + virtual long double d2rhormolardxi2__constxj(const std::vector &x, std::size_t i, x_N_dependency_flag xN_flag) = 0; + virtual long double d2rhormolardxidxj(const std::vector &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) = 0; + virtual long double d2Trdxi2__constxj(const std::vector &x, std::size_t i, x_N_dependency_flag xN_flag) = 0; + virtual long double d2Trdxidxj(const std::vector &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 &x, int i, int j); + long double d_ndTrdni_dxj__constxi(const std::vector &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 &x, int i, int j); + long double d_ndrhorbardni_dxj__constxi(const std::vector &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); - long double ndrhorbardni__constnj(const std::vector &x, int i); - long double ndTrdni__constnj(const std::vector &x, int i); + long double ndrhorbardni__constnj(const std::vector &x, std::size_t i, x_N_dependency_flag xN_flag); + long double ndTrdni__constnj(const std::vector &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 pFluids; ///< List of pointers to fluids + std::vector Yc_T; ///< Vector of critical temperatures for all components + std::vector Yc_v; ///< Vector of critical molar volumes for all components + std::vector pFluids; ///< List of postd::size_ters to fluids public: GERG2008ReducingFunction(std::vector 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(N,0)); v_c.resize(N,std::vector(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 &x); /// The derivative of reduced temperature with respect to component i mole fraction - long double dTrdxi__constxj(const std::vector &x, int i); + long double dTrdxi__constxj(const std::vector &x, std::size_t i, x_N_dependency_flag xN_flag); /// The molar reducing density long double rhormolar(const std::vector &x); ///Derivative of the molar reducing density with respect to component i mole fraction - long double drhormolardxi__constxj(const std::vector &x, int i); - long double dvrmolardxi__constxj(const std::vector &x, int i); + long double drhormolardxi__constxj(const std::vector &x, std::size_t i, x_N_dependency_flag xN_flag); + long double dvrmolardxi__constxj(const std::vector &x, std::size_t i, x_N_dependency_flag xN_fla); - long double d2vrmolardxi2__constxj(const std::vector &x, int i); - long double d2rhormolardxi2__constxj(const std::vector &x, int i); - long double d2vrmolardxidxj(const std::vector &x, int i, int j); - long double d2rhormolardxidxj(const std::vector &x, int i, int j); - long double d2Trdxi2__constxj(const std::vector &x, int i); - long double d2Trdxidxj(const std::vector &x, int i, int j); + long double d2vrmolardxi2__constxj(const std::vector &x, std::size_t i, x_N_dependency_flag xN_flag); + long double d2rhormolardxi2__constxj(const std::vector &x, std::size_t i, x_N_dependency_flag xN_flag); + long double d2vrmolardxidxj(const std::vector &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); + long double d2rhormolardxidxj(const std::vector &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); + long double d2Trdxi2__constxj(const std::vector &x, std::size_t i, x_N_dependency_flag xN_flag); + long double d2Trdxidxj(const std::vector &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 &x, int i, int j, std::vector< std::vector< long double> > &beta); + long double Yr(const std::vector &x, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector &Yc); + long double dYrdxi__constxj(const std::vector &x, std::size_t i, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector &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 &x, std::size_t i, std::size_t j, const STLMatrix &beta); - long double dfYkidxi__constxk(const std::vector &x, int k, int i,std::vector< std::vector< long double> > &beta); - long double dfYikdxi__constxk(const std::vector &x, int i, int k, std::vector< std::vector< long double> > &beta); - long double d2fYkidxi2__constxk(const std::vector &x, int k, int i, std::vector< std::vector< long double> > &beta); - long double d2fYikdxi2__constxk(const std::vector &x, int i, int k, std::vector< std::vector< long double> > &beta); - long double d2fYijdxidxj(const std::vector &x, int i, int k, std::vector< std::vector< long double> > &beta); + long double dfYkidxi__constxk(const std::vector &x, std::size_t k, std::size_t i, const STLMatrix &beta); + long double dfYikdxi__constxk(const std::vector &x, std::size_t i, std::size_t k, const STLMatrix &beta); + long double d2fYkidxi2__constxk(const std::vector &x, std::size_t k, std::size_t i, const STLMatrix &beta); + long double d2fYikdxi2__constxk(const std::vector &x, std::size_t i, std::size_t k, const STLMatrix &beta); + long double d2fYijdxidxj(const std::vector &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 &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 &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)); }; }; diff --git a/src/Backends/Helmholtz/VLERoutines.cpp b/src/Backends/Helmholtz/VLERoutines.cpp index 09f81284..a0b1c18c 100644 --- a/src/Backends/Helmholtz/VLERoutines.cpp +++ b/src/Backends/Helmholtz/VLERoutines.cpp @@ -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(N, 0)); -// -// for (std::size_t j = 0; j < N; ++j) -// { -// std::vector KK = K; -// KK[j] += 1e-6; -// build_arrays(HEOS,beta0,T0,rhomolar_liq0,rhomolar_vap0,z,KK); -// std::vector 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 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 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 &z, std::vector &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 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); } }