Merge pull request #2021 from PENGYAN777/feature_coolprop

CoolProp enhancement

Author: PENGYAN777

SU2_CFD/include/fluid/CCoolProp.hpp

@@ -48,6 +48,7 @@ class CCoolProp final : public CFluidModel {
   su2double Temperature_Critical{0.0}; /*!< \brief critical temperature */
   su2double acentric_factor{0.0};      /*!< \brief acentric factor */
   const su2double dp{0.01};            /*!< threshold for pressure */
+  const su2double dt{0.01};            /*!< threshold for temperature */
 #ifdef USE_COOLPROP
   std::unique_ptr<CoolProp::AbstractState> fluid_entity; /*!< \brief fluid entity */
 #endif
@@ -62,6 +63,18 @@ class CCoolProp final : public CFluidModel {
       Pressure = fmin(Pressure, (1 - dp) * Pressure_Critical);
   }
 
+  /*!
+   * \brief Avoid critical temperature
+   * \param[in,out] Temperature: Modified so that it is not too close to critical temperature to avoid issues in
+   * CoolProp.
+   */
+  void CheckTemperature(su2double& Temperature) const {
+    if (Temperature > Temperature_Critical)
+      Temperature = fmax(Temperature, (1 + dt) * Temperature_Critical);
+    else
+      Temperature = fmin(Temperature, (1 + dt) * Temperature_Critical);
+  }
+
  public:
   /*!
    * \brief Constructor of the class.
@@ -155,4 +168,4 @@ class CCoolProp final : public CFluidModel {
    * \return Value of the constant: Gamma
    */
   su2double GetGamma(void) const { return Gamma; }
-};
+};

SU2_CFD/src/fluid/CCoolProp.cpp

@@ -56,25 +56,20 @@ void CCoolProp::SetTDState_rhoe(su2double rho, su2double e) {
   dTdrho_e = fluid_entity->first_partial_deriv(CoolProp::iT, CoolProp::iDmass, CoolProp::iUmass);
   dTde_rho = fluid_entity->first_partial_deriv(CoolProp::iT, CoolProp::iUmass, CoolProp::iDmass);
   if (fluid_entity->phase() == CoolProp::iphase_twophase) {
-    // assume it is pure gas
-    fluid_entity->specify_phase(CoolProp::iphase_gas);
+    // impose gas phase
+    Temperature = Temperature + 0.1;
     CheckPressure(Pressure);
+    CheckTemperature(Temperature);
     fluid_entity->update(CoolProp::PT_INPUTS, Pressure, Temperature);
-    if (abs(fluid_entity->rhomass() / Density - 1) < dp) {
-      // origial phase is near saturation gas, then just compute sound speed
-      SoundSpeed2 = pow(fluid_entity->speed_sound(), 2);
-    } else {
-      // original phase is not near saturation gas, then specify the phase as gas phase
-      fluid_entity->specify_phase(CoolProp::iphase_gas);
-      SetTDState_PT(Pressure, Temperature);
-    }
+    SoundSpeed2 = pow(fluid_entity->speed_sound(), 2);
   } else {
     SoundSpeed2 = pow(fluid_entity->speed_sound(), 2);
   }
 }
 
 void CCoolProp::SetTDState_PT(su2double P, su2double T) {
   CheckPressure(P);
+  CheckTemperature(T);
   fluid_entity->update(CoolProp::PT_INPUTS, P, T);
   su2double rho = fluid_entity->rhomass();
   su2double e = fluid_entity->umass();
@@ -130,4 +125,4 @@ CCoolProp::CCoolProp(const string& fluidname) {
       "or autodiff",
       CURRENT_FUNCTION);
 }
-#endif
+#endif