From ea1485bb6fba79aedf875e4e8b96c09957e71da3 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Fri, 29 Jan 2016 10:06:35 -0700 Subject: [PATCH] Allow gibbs as input to first_partial_deriv; closes #951 --- src/AbstractState.cpp | 105 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 105 insertions(+) diff --git a/src/AbstractState.cpp b/src/AbstractState.cpp index 9755b085..dccf7f36 100644 --- a/src/AbstractState.cpp +++ b/src/AbstractState.cpp @@ -579,6 +579,22 @@ void get_dT_drho(AbstractState &AS, parameters index, CoolPropDbl &dT, CoolPropD } break; } + case iGmass: + case iGmolar: + { + // dg/dT|rho + double dTau_dT = 1/dT_dtau; + dT = R*AS.T()*(AS.dalpha0_dTau()+AS.dalphar_dTau()+AS.delta()*AS.d2alphar_dDelta_dTau())*dTau_dT + R*(1+AS.alpha0() + AS.alphar() + AS.delta()*AS.dalphar_dDelta()); + // dg/drho|T + double dDelta_drho = 1/rhor; + drho = AS.T()*R*(AS.dalpha0_dDelta()+AS.dalphar_dDelta()+AS.delta()*AS.d2alphar_dDelta2() + AS.dalphar_dDelta())*dDelta_drho; + if (index == iGmass){ + // dg/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass + drho /= AS.molar_mass(); + dT /= AS.molar_mass(); + } + break; + } case iTau: dT = 1/dT_dtau; drho = 0; break; case iDelta: @@ -765,6 +781,95 @@ TEST_CASE("Check AbstractState","[AbstractState]") { CHECK_NOTHROW(shared_ptr Water(CoolProp::AbstractState::factory("REFPROP", "Water"))); } +} + +TEST_CASE("Check derivatives in first_partial_deriv","[derivs_in_first_partial_deriv]") +{ + shared_ptr Water(CoolProp::AbstractState::factory("HEOS", "Water")); + shared_ptr WaterplusT(CoolProp::AbstractState::factory("HEOS", "Water")); + shared_ptr WaterminusT(CoolProp::AbstractState::factory("HEOS", "Water")); + shared_ptr Waterplusrho(CoolProp::AbstractState::factory("HEOS", "Water")); + shared_ptr Waterminusrho(CoolProp::AbstractState::factory("HEOS", "Water")); + + double dT = 1e-3, drho = 1; + Water->update(CoolProp::PT_INPUTS, 101325, 300); + WaterplusT->update(CoolProp::DmolarT_INPUTS, Water->rhomolar(), 300+dT); + WaterminusT->update(CoolProp::DmolarT_INPUTS, Water->rhomolar(), 300-dT); + Waterplusrho->update(CoolProp::DmolarT_INPUTS, Water->rhomolar()+drho, 300); + Waterminusrho->update(CoolProp::DmolarT_INPUTS, Water->rhomolar()-drho, 300); + + // Numerical derivatives + CoolPropDbl dP_dT_num = (WaterplusT->p() - WaterminusT->p())/(2*dT); + CoolPropDbl dP_drho_num = (Waterplusrho->p() - Waterminusrho->p())/(2*drho); + + CoolPropDbl dHmolar_dT_num = (WaterplusT->hmolar() - WaterminusT->hmolar())/(2*dT); + CoolPropDbl dHmolar_drho_num = (Waterplusrho->hmolar() - Waterminusrho->hmolar())/(2*drho); + CoolPropDbl dHmass_dT_num = (WaterplusT->hmass() - WaterminusT->hmass())/(2*dT); + CoolPropDbl dHmass_drho_num = (Waterplusrho->hmass() - Waterminusrho->hmass())/(2*drho); + + CoolPropDbl dSmolar_dT_num = (WaterplusT->smolar() - WaterminusT->smolar())/(2*dT); + CoolPropDbl dSmolar_drho_num = (Waterplusrho->smolar() - Waterminusrho->smolar())/(2*drho); + CoolPropDbl dSmass_dT_num = (WaterplusT->smass() - WaterminusT->smass())/(2*dT); + CoolPropDbl dSmass_drho_num = (Waterplusrho->smass() - Waterminusrho->smass())/(2*drho); + + CoolPropDbl dUmolar_dT_num = (WaterplusT->umolar() - WaterminusT->umolar())/(2*dT); + CoolPropDbl dUmolar_drho_num = (Waterplusrho->umolar() - Waterminusrho->umolar())/(2*drho); + CoolPropDbl dUmass_dT_num = (WaterplusT->umass() - WaterminusT->umass())/(2*dT); + CoolPropDbl dUmass_drho_num = (Waterplusrho->umass() - Waterminusrho->umass())/(2*drho); + + CoolPropDbl dGmolar_dT_num = (WaterplusT->gibbsmolar() - WaterminusT->gibbsmolar())/(2*dT); + CoolPropDbl dGmolar_drho_num = (Waterplusrho->gibbsmolar() - Waterminusrho->gibbsmolar())/(2*drho); + CoolPropDbl dGmass_dT_num = (WaterplusT->gibbsmass() - WaterminusT->gibbsmass())/(2*dT); + CoolPropDbl dGmass_drho_num = (Waterplusrho->gibbsmass() - Waterminusrho->gibbsmass())/(2*drho); + + // Pressure + CoolPropDbl dP_dT_analyt, dP_drho_analyt; + CoolProp::get_dT_drho(*Water, CoolProp::iP, dP_dT_analyt, dP_drho_analyt); + // Enthalpy + CoolPropDbl dHmolar_dT_analyt, dHmolar_drho_analyt; + CoolProp::get_dT_drho(*Water, CoolProp::iHmolar, dHmolar_dT_analyt, dHmolar_drho_analyt); + CoolPropDbl dHmass_dT_analyt, dHmass_drho_analyt; + CoolProp::get_dT_drho(*Water, CoolProp::iHmass, dHmass_dT_analyt, dHmass_drho_analyt); + // Entropy + CoolPropDbl dSmolar_dT_analyt, dSmolar_drho_analyt; + CoolProp::get_dT_drho(*Water, CoolProp::iSmolar, dSmolar_dT_analyt, dSmolar_drho_analyt); + CoolPropDbl dSmass_dT_analyt, dSmass_drho_analyt; + CoolProp::get_dT_drho(*Water, CoolProp::iSmass, dSmass_dT_analyt, dSmass_drho_analyt); + // Internal energy + CoolPropDbl dUmolar_dT_analyt, dUmolar_drho_analyt; + CoolProp::get_dT_drho(*Water, CoolProp::iUmolar, dUmolar_dT_analyt, dUmolar_drho_analyt); + CoolPropDbl dUmass_dT_analyt, dUmass_drho_analyt; + CoolProp::get_dT_drho(*Water, CoolProp::iUmass, dUmass_dT_analyt, dUmass_drho_analyt); + // Gibbs + CoolPropDbl dGmolar_dT_analyt, dGmolar_drho_analyt; + CoolProp::get_dT_drho(*Water, CoolProp::iGmolar, dGmolar_dT_analyt, dGmolar_drho_analyt); + CoolPropDbl dGmass_dT_analyt, dGmass_drho_analyt; + CoolProp::get_dT_drho(*Water, CoolProp::iGmass, dGmass_dT_analyt, dGmass_drho_analyt); + + double eps = 1e-3; + + CHECK( std::abs(dP_dT_analyt/dP_dT_num-1) < eps); + CHECK( std::abs(dP_drho_analyt/dP_drho_num-1) < eps); + + CHECK( std::abs(dHmolar_dT_analyt/dHmolar_dT_num-1) < eps); + CHECK( std::abs(dHmolar_drho_analyt/dHmolar_drho_num-1) < eps); + CHECK( std::abs(dHmass_dT_analyt/dHmass_dT_num-1) < eps); + CHECK( std::abs(dHmass_drho_analyt/dHmass_drho_num-1) < eps); + + CHECK( std::abs(dSmolar_dT_analyt/dSmolar_dT_num-1) < eps); + CHECK( std::abs(dSmolar_drho_analyt/dSmolar_drho_num-1) < eps); + CHECK( std::abs(dSmass_dT_analyt/dSmass_dT_num-1) < eps); + CHECK( std::abs(dSmass_drho_analyt/dSmass_drho_num-1) < eps); + + CHECK( std::abs(dUmolar_dT_analyt/dUmolar_dT_num-1) < eps); + CHECK( std::abs(dUmolar_drho_analyt/dUmolar_drho_num-1) < eps); + CHECK( std::abs(dUmass_dT_analyt/dUmass_dT_num-1) < eps); + CHECK( std::abs(dUmass_drho_analyt/dUmass_drho_num-1) < eps); + + CHECK( std::abs(dGmolar_dT_analyt/dGmolar_dT_num-1) < eps); + CHECK( std::abs(dGmolar_drho_analyt/dGmolar_drho_num-1) < eps); + CHECK( std::abs(dGmass_dT_analyt/dGmass_dT_num-1) < eps); + CHECK( std::abs(dGmass_drho_analyt/dGmass_drho_num-1) < eps); }