From bdd46a3f4c7f0c776e50ff34434bebc718d51614 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Wed, 7 Jan 2015 20:18:02 -0700 Subject: [PATCH] Some work on derivatives, drho/dp|h doesn't seem to work Signed-off-by: Ian Bell --- .../Helmholtz/HelmholtzEOSMixtureBackend.cpp | 83 ++++++++++++++----- src/Tests/CoolProp-Tests.cpp | 34 ++++++-- 2 files changed, 91 insertions(+), 26 deletions(-) diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 92462855..97624ff4 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -2635,6 +2635,25 @@ long double HelmholtzEOSMixtureBackend::calc_second_two_phase_deriv(parameters O long double d_dvdh_dp__consth = (denominator*dnumerator - numerator*ddenominator)/POW2(denominator); return -POW2(rhomolar())*d_dvdh_dp__consth + dv_dh_constp*(-2*rhomolar())*drhomolar_dp__consth; } + else if (Of == iDmass && ((Wrt1 == iHmass && Constant1 == iP && Wrt2 == iP && Constant2 == iHmass) || (Wrt2 == iHmass && Constant2 == iP && Wrt1 == iP && Constant1 == iHmass))){ + parameters h_key = iHmass, rho_key = iDmass, p_key = iP; + long double rho = keyed_output(rho_key); + // taking the derivative of (drho/dv)*(dv/dh|p) with respect to p with h constant + long double dv_dh_constp = calc_first_two_phase_deriv(rho_key,h_key,p_key)/(-POW2(rho)); + long double drho_dp__consth = first_two_phase_deriv(rho_key, p_key, h_key); + + // Calculate the derivative of dvdh|p with respect to p at constant h + long double dhL_dp_sat = SatL->calc_first_saturation_deriv(h_key, p_key, *SatL, *SatV); + long double dhV_dp_sat = SatV->calc_first_saturation_deriv(h_key, p_key, *SatL, *SatV); + long double drhoL_dp_sat = SatL->calc_first_saturation_deriv(rho_key, p_key, *SatL, *SatV); + long double drhoV_dp_sat = SatV->calc_first_saturation_deriv(rho_key, p_key, *SatL, *SatV); + long double numerator = 1/SatV->keyed_output(rho_key) - 1/SatL->keyed_output(rho_key); + long double denominator = SatV->keyed_output(h_key) - SatL->keyed_output(h_key); + long double dnumerator = -1/POW2(SatV->keyed_output(rho_key))*drhoV_dp_sat + 1/POW2(SatL->keyed_output(rho_key))*drhoL_dp_sat; + long double ddenominator = dhV_dp_sat - dhL_dp_sat; + long double d_dvdh_dp__consth = (denominator*dnumerator - numerator*ddenominator)/POW2(denominator); + return -POW2(rho)*d_dvdh_dp__consth + dv_dh_constp*(-2*rho)*drho_dp__consth; + } else{ throw ValueError(); } @@ -2704,14 +2723,17 @@ long double HelmholtzEOSMixtureBackend::calc_first_two_phase_deriv_splined(param // Faking single-phase long double rho_liq = Liq->keyed_output(rho_key); - long double drho_dh_liq = Liq->first_partial_deriv(rho_key, h_key, p_key); + long double drho_dh_liq__constp = Liq->first_partial_deriv(rho_key, h_key, p_key); - // Spline coordinates - long double a = 1/POW3(Delta_end) * (2*rho_liq - 2*rho_end + Delta_end * (drho_dh_liq + drho_dh_end)); - long double b = 3/POW2(Delta_end) * (-rho_liq + rho_end) - 1/Delta_end * (drho_dh_end + 2 * drho_dh_liq); - long double c = drho_dh_liq; + // Spline coordinates a, b, c, d + long double Abracket = (2*rho_liq - 2*rho_end + Delta_end * (drho_dh_liq__constp + drho_dh_end)); + long double a = 1/POW3(Delta_end) * Abracket; + long double b = 3/POW2(Delta_end) * (-rho_liq + rho_end) - 1/Delta_end * (drho_dh_end + 2 * drho_dh_liq__constp); + long double c = drho_dh_liq__constp; long double d = rho_liq; + std::cout <second_partial_deriv(rho_key, h_key, p_key, p_key, h_key); - long double drho_dp_liq = Liq->first_partial_deriv(rho_key, p_key, h_key); + // At the saturated state + long double rhoL = SatL->keyed_output(rho_key); + long double rhoV = SatV->keyed_output(rho_key); + long double hL = SatL->keyed_output(h_key); + long double hV = SatV->keyed_output(h_key); // Derivatives *along* the saturation curve using the special internal method long double dhL_dp_sat = SatL->calc_first_saturation_deriv(h_key, p_key, *SatL, *SatV); long double dhV_dp_sat = SatV->calc_first_saturation_deriv(h_key, p_key, *SatL, *SatV); - long double drhoL_dp_sat = SatL->calc_first_saturation_deriv(rho_key, p_key, *SatL, *SatV); + long double drhoL_dp_sat = SatL->calc_first_saturation_deriv(rho_key, p_key, *SatL, *SatV); + long double drhoV_dp_sat = SatV->calc_first_saturation_deriv(rho_key, p_key, *SatL, *SatV); + + long double drho_dp_end = POW2(End->keyed_output(rho_key))*(x_end/POW2(rhoV)*drhoV_dp_sat + (1-x_end)/POW2(rhoL)*drhoL_dp_sat); + + // Faking single-phase + long double drho_dp__consth_liq = Liq->first_partial_deriv(rho_key, p_key, h_key); + long double d2rhodhdp_liq = Liq->second_partial_deriv(rho_key, h_key, p_key, p_key, h_key); // ? + + //d2rhodhdp_liq = -5.88936E-08; // Derivatives at the end point - long double drho_dp_end = End->calc_first_two_phase_deriv(rho_key, p_key, h_key); - long double drhodhdp_end = End->calc_second_two_phase_deriv(rho_key, h_key, p_key, p_key, h_key); + long double drho_dp__consth_end = End->calc_first_two_phase_deriv(rho_key, p_key, h_key); + long double d2rhodhdp_end = End->calc_second_two_phase_deriv(rho_key, h_key, p_key, p_key, h_key); - long double dDelta_end_dp = x_end * (dhL_dp_sat - dhV_dp_sat); - long double dDelta_dp = -dhL_dp_sat; + // Reminder: + // Delta = Q()*(hV-hL) = h-hL + // Delta_end = x_end*(hV-hL); + long double d_Delta_dp__consth = -dhL_dp_sat; + long double d_Delta_end_dp__consth = x_end*(dhV_dp_sat - dhL_dp_sat); - // First pressure derivative at constant h of the coefficients - long double da_dp = (-6/POW4(Delta_end)) * dDelta_end_dp * (rho_liq - rho_end) + (2/POW3(Delta_end)) * (drho_dp_liq - drho_dp_end) + (1/POW2(Delta_end)) * (drhodhdp_liq + drhodhdp_end) - (2/POW3(Delta_end)) * (drho_dh_liq + drho_dh_end)*dDelta_end_dp; - long double db_dp = (-6/POW3(Delta_end)) * dDelta_end_dp * (-rho_liq + rho_end) + (3/POW2(Delta_end)) * (-drho_dp_liq + drho_dp_end) + (1/POW2(Delta_end)) * dDelta_end_dp * (drho_dh_end + 2 * drho_dh_liq) - (1/Delta_end) * (drhodhdp_end + 2*drhodhdp_liq); - long double dc_dp = drhodhdp_liq; - long double dd_dp = drho_dp_liq; + // First pressure derivative at constant h of the coefficients a,b,c,d + // long double Abracket = (2*rho_liq - 2*rho_end + Delta_end * (drho_dh_liq__constp + drho_dh_end)); + long double d_Abracket_dp_consth = (2*drhoL_dp_sat - 2*drho_dp__consth_end + Delta_end*(d2rhodhdp_liq + d2rhodhdp_end) + d_Delta_dp__consth*(drho_dh_liq__constp + drho_dh_end)); + long double da_dp = 1/POW3(Delta_end)*d_Abracket_dp_consth + Abracket*(-3/POW4(Delta_end)*d_Delta_end_dp__consth); + long double db_dp = - 6/POW3(Delta_end)*d_Delta_end_dp__consth*(rho_end - rho_liq) + + (3/POW2(Delta_end))*(drho_dp__consth_end - drhoL_dp_sat) + + (1/POW2(Delta_end)*d_Delta_end_dp__consth) * (drho_dh_end + 2*drho_dh_liq__constp) + - (1/Delta_end) * (d2rhodhdp_end + 2*d2rhodhdp_liq); + long double dc_dp = d2rhodhdp_liq; + long double dd_dp = drhoL_dp_sat; - long double d_rho_spline_dp__consth = (3*a*POW2(Delta) + 2*b*Delta + c)*dDelta_dp + POW3(Delta)*da_dp + POW2(Delta)*db_dp + Delta*dc_dp + dd_dp; + std::cout < AS(new CoolProp::HelmholtzEOSBackend("n-Propane")); AS->update(QT_INPUTS, 0.3, 300); long double analytical = AS->second_two_phase_deriv(iDmolar, iHmolar, iP, iP, iHmolar); @@ -1565,16 +1565,34 @@ TEST_CASE("Check the second two-phase derivative", "[second_two_phase_deriv]") long double v2 = AS->first_two_phase_deriv(iDmolar, iHmolar, iP); long double numerical = (v1 - v2)/(pplus - pminus); CAPTURE(numerical); - CHECK(std::abs(numerical/analytical-1) < 1e-40); + CHECK(std::abs(numerical/analytical-1) < 1e-6); + } + SECTION("d2rhodhdp using mass",""){ + shared_ptr AS(new CoolProp::HelmholtzEOSBackend("n-Propane")); + AS->update(QT_INPUTS, 0.3, 300); + long double analytical = AS->second_two_phase_deriv(iDmass, iHmass, iP, iP, iHmass); + CAPTURE(analytical); + long double pplus = AS->p()*1.001, pminus = AS->p()*0.999, h = AS->hmass(); + AS->update(HmassP_INPUTS, h, pplus); + long double v1 = AS->first_two_phase_deriv(iDmass, iHmass, iP); + AS->update(HmassP_INPUTS, h, pminus); + long double v2 = AS->first_two_phase_deriv(iDmass, iHmass, iP); + long double numerical = (v1 - v2)/(pplus - pminus); + CAPTURE(numerical); + CHECK(std::abs(numerical/analytical-1) < 1e-6); } } -TEST_CASE("Check the first two-phase derivative using splines", "[first_two_phase_deriv]") +TEST_CASE("Check the first two-phase derivative using splines", "[first_two_phase_deriv_splined]") { const int number_of_pairs = 4; struct pair {parameters p1, p2, p3;}; - pair pairs[number_of_pairs] = {{iDmass, iP, iHmass}, {iDmolar, iP, iHmolar}, - {iDmolar, iHmolar, iP}, {iDmass, iHmass, iP}}; + pair pairs[number_of_pairs] = { + {iDmass, iP, iHmass}, + {iDmolar, iP, iHmolar}, + {iDmolar, iHmolar, iP}, + {iDmass, iHmass, iP} + }; shared_ptr AS(new CoolProp::HelmholtzEOSBackend("n-Propane")); for (std::size_t i = 0; i < number_of_pairs; ++i) { @@ -1583,7 +1601,7 @@ TEST_CASE("Check the first two-phase derivative using splines", "[first_two_phas ss1 << "for (" << get_parameter_information(pairs[i].p1,"short") << ", " << get_parameter_information(pairs[i].p2,"short") << ", " << get_parameter_information(pairs[i].p3,"short") << ")"; SECTION(ss1.str(),"") { - AS->update(QT_INPUTS, 0.15, 300); + AS->update(QT_INPUTS, 0.2, 300); long double numerical; long double analytical = AS->first_two_phase_deriv_splined(pairs[i].p1, pairs[i].p2, pairs[i].p3, 0.3); CAPTURE(analytical); @@ -1594,16 +1612,18 @@ TEST_CASE("Check the first two-phase derivative using splines", "[first_two_phas v3base = AS->keyed_output(pairs[i].p3); long double v2plus = v2base*1.001; long double v2minus = v2base*0.999; + CoolProp::input_pairs input_pair1 = generate_update_pair(pairs[i].p2, v2plus, pairs[i].p3, v3base, out1, out2); AS->update(input_pair1, out1, out2); long double v1 = AS->first_two_phase_deriv_splined(pairs[i].p1, pairs[i].p1, pairs[i].p1, 0.3); + CoolProp::input_pairs input_pair2 = generate_update_pair(pairs[i].p2, v2minus, pairs[i].p3, v3base, out1, out2); AS->update(input_pair2, out1, out2); long double v2 = AS->first_two_phase_deriv_splined(pairs[i].p1, pairs[i].p1, pairs[i].p1, 0.3); numerical = (v1 - v2)/(v2plus - v2minus); CAPTURE(numerical); - CHECK(std::abs(numerical/analytical-1) < 1e-4); + CHECK(std::abs(numerical/analytical-1) < 1e-10); } } }