single phase and away from Pc

Author: PENGYAN777

SU2_CFD/src/fluid/CCoolProp.cpp

@@ -1,29 +1,29 @@
 /*!
- * \file CCoolProp.cpp
- * \brief Source of the fluid model from CoolProp.
- * \author P. Yan, G. Gori, A. Guardone
- * \version 7.5.1 "Blackbird"
- *
- * SU2 Project Website: https://su2code.github.io
- *
- * The SU2 Project is maintained by the SU2 Foundation
- * (http://su2foundation.org)
- *
- * Copyright 2012-2023, SU2 Contributors (cf. AUTHORS.md)
- *
- * SU2 is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public
- * License as published by the Free Software Foundation; either
- * version 2.1 of the License, or (at your option) any later version.
- *
- * SU2 is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
- */
+* \file CCoolProp.cpp
+* \brief Source of the fluid model from CoolProp.
+* \author P. Yan, G. Gori, A. Guardone
+* \version 7.4.0 "Blackbird"
+*
+* SU2 Project Website: https://su2code.github.io
+*
+* The SU2 Project is maintained by the SU2 Foundation
+* (http://su2foundation.org)
+*
+* Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md)
+*
+* SU2 is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version.
+*
+* SU2 is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+* Lesser General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public
+* License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+*/
 
 #include "../../include/fluid/CCoolProp.hpp"
 
@@ -56,28 +56,51 @@ 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() == 6) {
+    // assume it is pure gas
+    fluid_entity->specify_phase(CoolProp::iphase_gas);
+    if(abs(Pressure/Pressure_Critical-1)<0.01){
+      // if P is very cloase to Pc, the CoolProp has issue
+      Pressure = Pressure_Critical*0.99;
+    }
+    fluid_entity->update(CoolProp::PT_INPUTS, Pressure, Temperature);
+    if(abs(fluid_entity->rhomass()/Density-1)*100 < 1){
+      // 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);
+    }
   }
   else{
-      SoundSpeed2 = pow(fluid_entity->speed_sound(), 2);
+    SoundSpeed2 = pow(fluid_entity->speed_sound(), 2);
   }
 }
 
 void CCoolProp::SetTDState_PT(su2double P, su2double T) {
+  if(abs(P/Pressure_Critical-1)<0.01) {
+    P = Pressure_Critical * 0.99;
+  }
   fluid_entity->update(CoolProp::PT_INPUTS, P, T);
   su2double rho = fluid_entity->rhomass();
   su2double e = fluid_entity->umass();
   SetTDState_rhoe(rho, e);
 }
 
 void CCoolProp::SetTDState_Prho(su2double P, su2double rho) {
+  if(abs(P/Pressure_Critical-1)<0.01){
+    P = Pressure_Critical*0.99;
+  }
   fluid_entity->update(CoolProp::DmassP_INPUTS, rho, P);
   su2double e = fluid_entity->umass();
   SetTDState_rhoe(rho, e);
 }
 
-void CCoolProp::SetEnergy_Prho(su2double P, su2double rho) {
+void CCoolProp::SetEnergy_Prho(su2double P, su2double rho){
+  if(abs(P/Pressure_Critical-1)<0.01){
+    P = Pressure_Critical*0.99;
+  }
   fluid_entity->update(CoolProp::DmassP_INPUTS, rho, P);
   StaticEnergy = fluid_entity->umass();
 }
@@ -90,17 +113,19 @@ void CCoolProp::SetTDState_hs(su2double h, su2double s) {
 }
 
 void  CCoolProp::SetTDState_Ps(su2double P, su2double s) {
+  if(abs(P/Pressure_Critical-1)<0.01){
+    P = Pressure_Critical*0.99;
+  }
   fluid_entity->update(CoolProp::PSmass_INPUTS, P, s);
-  su2double Rho = fluid_entity->rhomass();
+  su2double rho = fluid_entity->rhomass();
   su2double e = fluid_entity->umass();
-  SetTDState_rhoe(Rho, e);
+  SetTDState_rhoe(rho, e);
 }
 
 void CCoolProp::SetTDState_rhoT(su2double rho, su2double T) {
   fluid_entity->update(CoolProp::DmassT_INPUTS, rho, T);
-  su2double Rho = fluid_entity->rhomass();
   su2double e = fluid_entity->umass();
-  SetTDState_rhoe(Rho, e);
+  SetTDState_rhoe(rho, e);
 }
 
 void  CCoolProp::ComputeDerivativeNRBC_Prho(su2double P, su2double rho) {