diff --git a/include/CoolPropFluid.h b/include/CoolPropFluid.h index 9c03a29f..e37b3c52 100644 --- a/include/CoolPropFluid.h +++ b/include/CoolPropFluid.h @@ -305,16 +305,12 @@ public: SaturationAncillaryFunction(rapidjson::Value &json_code) { - std::string type = cpjson::get_string(json_code,"type"); if (!type.compare("rational_polynomial")) { - std::vector A = cpjson::get_double_array(json_code["A"]); - std::vector B = cpjson::get_double_array(json_code["B"]); - std::reverse(A.begin(), A.end()); - std::reverse(B.begin(), B.end()); - num_coeffs = vec_to_eigen(A); - den_coeffs = vec_to_eigen(B); + num_coeffs = vec_to_eigen(cpjson::get_double_array(json_code["A"])); + den_coeffs = vec_to_eigen(cpjson::get_double_array(json_code["B"])); + max_abs_error = cpjson::get_double(json_code,"max_abs_error"); } else { @@ -500,7 +496,7 @@ public: struct Ancillaries { - SaturationAncillaryFunction pL, pV, rhoL, rhoV, hL, hV, sL, sV; + SaturationAncillaryFunction pL, pV, rhoL, rhoV, hL, hLV, sL, sLV; MeltingLineVariables melting_line; SurfaceTensionCorrelation surface_tension; }; diff --git a/src/Backends/Helmholtz/FlashRoutines.cpp b/src/Backends/Helmholtz/FlashRoutines.cpp index a918ee57..1ec5d8da 100644 --- a/src/Backends/Helmholtz/FlashRoutines.cpp +++ b/src/Backends/Helmholtz/FlashRoutines.cpp @@ -327,6 +327,95 @@ void FlashRoutines::PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, int other) } void FlashRoutines::HSU_P_flash_singlephase(HelmholtzEOSMixtureBackend &HEOS, int other, long double T0, long double rhomolar0) { + double A[2][2], B[2][2]; + long double y; + HelmholtzEOSMixtureBackend _HEOS(HEOS.get_components()); + _HEOS.update(DmolarT_INPUTS, rhomolar0, T0); + long double Tc = HEOS.calc_T_critical(); + long double rhoc = HEOS.calc_rhomolar_critical(); + long double R = HEOS.gas_constant(); + long double p = HEOS.p(); + switch (other) + { + case iHmolar: y = HEOS.hmolar(); break; + case iSmolar: y = HEOS.smolar(); break; + } + + long double worst_error = 999; + int iter = 0; + bool failed = false; + long double omega = 1.0, f2, df2_dtau, df2_ddelta; + long double tau = _HEOS.tau(), delta = _HEOS.delta(); + while (worst_error>1e-6 && failed == false) + { + + // All the required partial derivatives + long double a0 = _HEOS.calc_alpha0_deriv_nocache(0,0,HEOS.mole_fractions, tau, delta,Tc,rhoc); + long double da0_ddelta = _HEOS.calc_alpha0_deriv_nocache(0,1,HEOS.mole_fractions, tau, delta,Tc,rhoc); + long double da0_dtau = _HEOS.calc_alpha0_deriv_nocache(1,0,HEOS.mole_fractions, tau, delta,Tc,rhoc); + long double d2a0_dtau2 = _HEOS.calc_alpha0_deriv_nocache(2,0,HEOS.mole_fractions, tau, delta,Tc,rhoc); + long double d2a0_ddelta_dtau = 0.0; + + long double ar = _HEOS.calc_alphar_deriv_nocache(0,0,HEOS.mole_fractions, tau, delta); + long double dar_dtau = _HEOS.calc_alphar_deriv_nocache(1,0,HEOS.mole_fractions, tau, delta); + long double dar_ddelta = _HEOS.calc_alphar_deriv_nocache(0,1,HEOS.mole_fractions, tau, delta); + long double d2ar_ddelta_dtau = _HEOS.calc_alphar_deriv_nocache(1,1,HEOS.mole_fractions, tau, delta); + long double d2ar_ddelta2 = _HEOS.calc_alphar_deriv_nocache(0,2,HEOS.mole_fractions, tau, delta); + long double d2ar_dtau2 = _HEOS.calc_alphar_deriv_nocache(2,0,HEOS.mole_fractions, tau, delta); + + long double f1 = delta/tau*(1+delta*dar_ddelta)-p/(rhoc*R*Tc); + long double df1_dtau = (1+delta*dar_ddelta)*(-delta/tau/tau)+delta/tau*(delta*d2ar_ddelta_dtau); + long double df1_ddelta = (1.0/tau)*(1+2.0*delta*dar_ddelta+delta*delta*d2ar_ddelta2); + switch (other) + { + case iHmolar: + { + f2 = (1+delta*dar_ddelta)+tau*(da0_dtau+dar_dtau) - tau*y/(R*Tc); + df2_dtau = delta*d2ar_ddelta_dtau+da0_dtau+dar_dtau+tau*(d2a0_dtau2+d2ar_dtau2) - y/(R*Tc); + df2_ddelta = (dar_ddelta+delta*d2ar_ddelta2)+tau*(d2a0_ddelta_dtau+d2ar_ddelta_dtau); + break; + } + case iSmolar: + { + f2 = tau*(da0_dtau+dar_dtau)-ar-a0-y/R; + df2_dtau = tau*(d2a0_dtau2 + d2ar_dtau2)+(da0_dtau+dar_dtau)-dar_dtau-da0_dtau; + df2_ddelta = tau*(d2a0_ddelta_dtau+d2ar_ddelta_dtau)-dar_ddelta-da0_ddelta; + break; + } + default: + break; + } + + //First index is the row, second index is the column + A[0][0]=df1_dtau; + A[0][1]=df1_ddelta; + A[1][0]=df2_dtau; + A[1][1]=df2_ddelta; + + //double det = A[0][0]*A[1][1]-A[1][0]*A[0][1]; + + MatInv_2(A,B); + tau -= omega*(B[0][0]*f1+B[0][1]*f2); + delta -= omega*(B[1][0]*f1+B[1][1]*f2); + + if (fabs(f1)>fabs(f2)) + worst_error=fabs(f1); + else + worst_error=fabs(f2); + + if (!ValidNumber(f1) || !ValidNumber(f2)) + { + throw SolutionError(format("Invalid values for inputs p=%g y=%g for fluid %s in HSU_P_flash_singlephase", p, y, _HEOS.name().c_str())); + } + + iter += 1; + if (iter>100) + { + throw SolutionError(format("HSU_P_flash_singlephase did not converge with inputs p=%g h=%g for fluid %s", p, y,_HEOS.name().c_str())); + } + } + + HEOS.update(DmolarT_INPUTS, rhoc*delta, Tc/tau); } void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, int other) { @@ -354,28 +443,15 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, int other) } else { - HEOS._phase = iphase_gas; throw NotImplementedError("HSU_P_flash does not support mixtures (yet)"); - // Find the phase, while updating all internal variables possible } } - if (HEOS.isHomogeneousPhase() && !ValidNumber(HEOS._p)) + if (HEOS.isHomogeneousPhase() && !ValidNumber(HEOS._T)) { - switch (other) - { - case iDmolar: - break; - case iHmolar: - HEOS._rhomolar = HEOS.solver_for_rho_given_T_oneof_HSU(HEOS._T, HEOS._hmolar, iHmolar); break; - case iSmolar: - HEOS._rhomolar = HEOS.solver_for_rho_given_T_oneof_HSU(HEOS._T, HEOS._smolar, iSmolar); break; - case iUmolar: - HEOS._rhomolar = HEOS.solver_for_rho_given_T_oneof_HSU(HEOS._T, HEOS._umolar, iUmolar); break; - default: - break; - } - HEOS.calc_pressure(); + long double T0 = 300, rho0 = 100; + HSU_P_flash_singlephase(HEOS, other, T0, rho0); + HEOS._Q = -1; } } diff --git a/src/Backends/Helmholtz/Fluids/FluidLibrary.h b/src/Backends/Helmholtz/Fluids/FluidLibrary.h index c3f97cf3..440bd128 100644 --- a/src/Backends/Helmholtz/Fluids/FluidLibrary.h +++ b/src/Backends/Helmholtz/Fluids/FluidLibrary.h @@ -927,14 +927,31 @@ protected: fluid.ancillaries.rhoL = SaturationAncillaryFunction(ancillaries["rhoL"]); fluid.ancillaries.rhoV = SaturationAncillaryFunction(ancillaries["rhoV"]); - if (ancillaries.HasMember("hL")) - { + if (ancillaries.HasMember("hL")){ fluid.ancillaries.hL = SaturationAncillaryFunction(ancillaries["hL"]); } - else - { + else{ if (get_debug_level() > 0){ std::cout << "Missing hL ancillary for fluid " << fluid.name; } } + if (ancillaries.HasMember("hLV")){ + fluid.ancillaries.hLV = SaturationAncillaryFunction(ancillaries["hLV"]); + } + else{ + if (get_debug_level() > 0){ std::cout << "Missing hLV ancillary for fluid " << fluid.name; } + } + + if (ancillaries.HasMember("sL")){ + fluid.ancillaries.sL = SaturationAncillaryFunction(ancillaries["sL"]); + } + else{ + if (get_debug_level() > 0){ std::cout << "Missing sL ancillary for fluid " << fluid.name; } + } + if (ancillaries.HasMember("sLV")){ + fluid.ancillaries.sLV = SaturationAncillaryFunction(ancillaries["sLV"]); + } + else{ + if (get_debug_level() > 0){ std::cout << "Missing sLV ancillary for fluid " << fluid.name; } + } }; diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index e17ccaeb..098937d1 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -642,20 +642,34 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot case iHmolar: { long double h_liq = components[0]->ancillaries.hL.evaluate(_TLanc); - long double h_vap = components[0]->ancillaries.hV.evaluate(_T); + long double h_liq_error_band = components[0]->ancillaries.hL.get_max_abs_error(); + long double h_vap = h_liq + components[0]->ancillaries.hLV.evaluate(_TLanc); + long double h_vap_error_band = h_liq_error_band + components[0]->ancillaries.hLV.get_max_abs_error(); // Check if in range given the accuracy of the fit - if (value > h_vap + components[0]->ancillaries.hV.get_max_abs_error()){ + if (value > h_vap + h_vap_error_band){ this->_phase = iphase_gas; _Q = -1000; return; } - else if (value < h_liq - components[0]->ancillaries.hL.get_max_abs_error()){ + else if (value < h_liq - h_liq_error_band){ this->_phase = iphase_liquid; _Q = 1000; return; } break; } case iSmolar: { - // Add entropy ancillary code here + long double s_liq = components[0]->ancillaries.sL.evaluate(_TLanc); + long double s_liq_error_band = components[0]->ancillaries.sL.get_max_abs_error(); + long double s_vap = s_liq + components[0]->ancillaries.sLV.evaluate(_TLanc); + long double s_vap_error_band = s_liq_error_band + components[0]->ancillaries.sLV.get_max_abs_error(); + + // Check if in range given the accuracy of the fit + if (value > s_vap + s_vap_error_band){ + this->_phase = iphase_gas; _Q = -1000; return; + } + else if (value < s_liq - s_liq_error_band){ + this->_phase = iphase_liquid; _Q = 1000; return; + } + break; } case iUmolar: {