Added DQ as input pair; closes #822

CoolProp · Oct 23, 2015 · 37f75dd · 37f75dd
1 parent 108f3c4
commit 37f75dd
Show file tree

Hide file tree

Showing 5 changed files with 63 additions and 1 deletion.
diff --git a/include/DataStructures.h b/include/DataStructures.h
@@ -211,6 +211,8 @@ enum input_pairs{
     QSmass_INPUTS, ///< Molar quality, Entropy in J/kg/K
     HmolarQ_INPUTS, ///< Enthalpy in J/mol, Molar quality
     HmassQ_INPUTS, ///< Enthalpy in J/kg, Molar quality
+    DmolarQ_INPUTS, ///< Density in mol/m^3, Molar quality
+    DmassQ_INPUTS,  ///< Density in kg/m^3, Molar quality
 
     PT_INPUTS, ///< Pressure in Pa, Temperature in K
 
@@ -326,6 +328,12 @@ template<class T> CoolProp::input_pairs generate_update_pair(parameters key1, T
         else if (match_pair(key1, key2, iDmolar, iP, swap)){
             pair = DmolarP_INPUTS; // Molar density in mol/m^3, Pressure in Pa
         }
+        else if (match_pair(key1, key2, iDmass, iQ, swap)){
+            pair = DmassQ_INPUTS; // Mass density in kg/m^3, molar vapor quality
+        }
+        else if (match_pair(key1, key2, iDmolar, iQ, swap)){
+            pair = DmolarQ_INPUTS; // Molar density in mol/m^3, molar vapor quality
+        }
         else if (match_pair(key1, key2, iHmass, iP, swap)){
             pair = HmassP_INPUTS; // Enthalpy in J/kg, Pressure in Pa
         }

diff --git a/src/Backends/Helmholtz/FlashRoutines.cpp b/src/Backends/Helmholtz/FlashRoutines.cpp
@@ -253,6 +253,47 @@ void FlashRoutines::DP_flash(HelmholtzEOSMixtureBackend &HEOS)
         }
     }
 }
+
+class DQ_flash_residual : public FuncWrapper1DWithTwoDerivs
+{
+public:
+    HelmholtzEOSMixtureBackend &HEOS;
+    double rhomolar, Q_target;
+    DQ_flash_residual(HelmholtzEOSMixtureBackend &HEOS, double rhomolar, double Q_target) : HEOS(HEOS), rhomolar(rhomolar), Q_target(Q_target) {};
+    double call(double T){
+        HEOS.update(QT_INPUTS, 0, T); // Doesn't matter whether liquid or vapor, we are just doing a full VLE call for given T
+        double rhoL = HEOS.saturated_liquid_keyed_output(iDmolar);
+        double rhoV = HEOS.saturated_vapor_keyed_output(iDmolar);
+        return (1/rhomolar - 1/rhoL)/(1/rhoV - 1/rhoL) - Q_target;
+    }
+    double deriv(double T){
+        return _HUGE;
+    }
+    double second_deriv(double T){
+        return _HUGE;
+    }
+};
+
+void FlashRoutines::DQ_flash(HelmholtzEOSMixtureBackend &HEOS)
+{
+    SaturationSolvers::saturation_PHSU_pure_options options;
+    options.use_logdelta = false;
+    HEOS.specify_phase(iphase_twophase);
+    if (HEOS.is_pure_or_pseudopure){
+        double Tmax = HEOS.T_critical() - 0.1;
+        double Tmin = HEOS.Tmin();
+        std::string errstr;
+        DQ_flash_residual resid(HEOS, HEOS._rhomolar, HEOS._Q);
+        Brent(resid, Tmin, Tmax, DBL_EPSILON, 1e-10, 100, errstr);
+        HEOS._p = HEOS.SatV->p();
+        HEOS._T = HEOS.SatV->T();
+        HEOS._rhomolar = HEOS.SatV->rhomolar();
+        HEOS._phase = iphase_twophase;
+    }
+    else{
+        throw NotImplementedError("DQ_flash not ready for mixtures");
+    }
+}
 void FlashRoutines::HQ_flash(HelmholtzEOSMixtureBackend &HEOS, CoolPropDbl Tguess)
 {
     SaturationSolvers::saturation_PHSU_pure_options options;
@@ -277,7 +318,7 @@ void FlashRoutines::HQ_flash(HelmholtzEOSMixtureBackend &HEOS, CoolPropDbl Tgues
         HEOS._phase = iphase_twophase;
     }
     else{
-        throw NotImplementedError("QS_flash not ready for mixtures");
+        throw NotImplementedError("HQ_flash not ready for mixtures");
     }
 }
 void FlashRoutines::QS_flash(HelmholtzEOSMixtureBackend &HEOS)

diff --git a/src/Backends/Helmholtz/FlashRoutines.h b/src/Backends/Helmholtz/FlashRoutines.h
@@ -50,6 +50,10 @@ class FlashRoutines{
     /// @param HEOS The HelmholtzEOSMixtureBackend to be used
     static void QS_flash(HelmholtzEOSMixtureBackend &HEOS);
 
+    /// Flash for given molar density and (molar) quality
+    /// @param HEOS The HelmholtzEOSMixtureBackend to be used
+    static void FlashRoutines::DQ_flash(HelmholtzEOSMixtureBackend &HEOS);
+
     /// Flash for given molar enthalpy and (molar) quality
     /// @param HEOS The HelmholtzEOSMixtureBackend to be used
     /// @param Tguess (optional) The guess temperature in K to start from, ignored if < 0

diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp
@@ -1033,6 +1033,7 @@ void HelmholtzEOSMixtureBackend::mass_to_molar_inputs(CoolProp::input_pairs &inp
         case SmassT_INPUTS: ///< Entropy in J/kg/K, Temperature in K
         //case TUmass_INPUTS: ///< Temperature in K, Internal energy in J/kg (NOT CURRENTLY IMPLEMENTED)
         case DmassP_INPUTS: ///< Mass density in kg/m^3, Pressure in Pa
+        case DmassQ_INPUTS: ///< Mass density in kg/m^3, molar quality
         case HmassP_INPUTS: ///< Enthalpy in J/kg, Pressure in Pa
         case PSmass_INPUTS: ///< Pressure in Pa, Entropy in J/kg/K
         case PUmass_INPUTS: ///< Pressure in Pa, Internal energy in J/kg
@@ -1055,6 +1056,7 @@ void HelmholtzEOSMixtureBackend::mass_to_molar_inputs(CoolProp::input_pairs &inp
                 case SmassT_INPUTS: input_pair = SmolarT_INPUTS; value1 *= mm;  break;
                 //case TUmass_INPUTS: input_pair = TUmolar_INPUTS; value2 *= mm;  break; (NOT CURRENTLY IMPLEMENTED)
                 case DmassP_INPUTS: input_pair = DmolarP_INPUTS; value1 /= mm;  break;
+                case DmassQ_INPUTS: input_pair = DmolarQ_INPUTS; value1 /= mm;  break;
                 case HmassP_INPUTS: input_pair = HmolarP_INPUTS; value1 *= mm;  break;
                 case PSmass_INPUTS: input_pair = PSmolar_INPUTS; value2 *= mm;  break;
                 case PUmass_INPUTS: input_pair = PUmolar_INPUTS; value2 *= mm;  break;
@@ -1136,6 +1138,8 @@ void HelmholtzEOSMixtureBackend::update(CoolProp::input_pairs input_pair, double
             _Q = value1; _smolar = value2; FlashRoutines::QS_flash(*this); break;
         case HmolarQ_INPUTS:
             _hmolar = value1; _Q = value2; FlashRoutines::HQ_flash(*this); break;
+        case DmolarQ_INPUTS:
+            _rhomolar = value1; _Q = value2; FlashRoutines::DQ_flash(*this); break;
         default:
             throw ValueError(format("This pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str()));
     }

diff --git a/src/DataStructures.cpp b/src/DataStructures.cpp
@@ -406,6 +406,9 @@ const input_pair_info input_pair_list[] = {
     { QSmass_INPUTS,       "QS_INPUTS",           "Molar quality, Entropy in J/kg/K"                   },
     { HmolarQ_INPUTS,      "HQ_INPUTS",           "Enthalpy in J/mol, Molar quality"                   },
     { HmassQ_INPUTS,       "HQ_INPUTS",           "Enthalpy in J/kg, Molar quality"                    },
+    { DmassQ_INPUTS,       "DmassQ_INPUTS",       "Molar density kg/m^3, Molar quality"               },
+    { DmolarQ_INPUTS,      "DmolarQ_INPUTS",      "Molar density in mol/m^3, Molar quality"            },
+
     { PQ_INPUTS,           "PQ_INPUTS",           "Pressure in Pa, Molar quality"                      },
 
     { PT_INPUTS,           "PT_INPUTS",           "Pressure in Pa, Temperature in K"                   },
@@ -500,6 +503,8 @@ void split_input_pair(input_pairs pair, parameters &p1, parameters &p2)
         case TUmolar_INPUTS: p1 = iT; p2 = iUmolar; break;
         case DmassP_INPUTS: p1 = iDmass; p2 = iP; break;
         case DmolarP_INPUTS: p1 = iDmolar; p2 = iP; break;
+        case DmassQ_INPUTS: p1 = iDmass; p2 = iQ; break;
+        case DmolarQ_INPUTS: p1 = iDmolar; p2 = iQ; break;
         case HmassP_INPUTS: p1 = iHmass; p2 = iP; break;
         case HmolarP_INPUTS: p1 = iHmolar; p2 = iP; break;
         case PSmass_INPUTS: p1 = iP; p2 = iSmass; break;