Fix enthalpy related issues

.github/workflows/regression.yml

@@ -211,7 +211,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:250717-1402
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}}
+          args: -b ${{github.ref}} -t fix_enthalpy -c fix_enthalpy -s ${{matrix.testscript}}
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:250717-1402
         with:

Common/include/CConfig.hpp

@@ -927,7 +927,6 @@ class CConfig {
   Initial_BCThrust,                /*!< \brief Ratio of turbulent to laminar viscosity at the actuator disk. */
   Pressure_FreeStream,             /*!< \brief Total pressure of the fluid. */
   Pressure_Thermodynamic,          /*!< \brief Thermodynamic pressure of the fluid. */
-  Standard_Ref_Temperature,        /*!< \brief Standard reference temperature for multicomponent flows. */
   Temperature_FreeStream,          /*!< \brief Total temperature of the fluid.  */
   Temperature_ve_FreeStream;       /*!< \brief Total vibrational-electronic temperature of the fluid.  */
   unsigned short wallModel_MaxIter; /*!< \brief maximum number of iterations for the Newton method for the wall model */
@@ -1972,12 +1971,6 @@ class CConfig {
    */
   su2double GetPressure_Thermodynamic(void) const { return Pressure_Thermodynamic; }
 
-  /*!
-   * \brief Get the value of the standard reference temperature for multicomponent flows.
-   * \return Standard reference temperature, Non-dimensionalized if it is needed for Non-Dimensional problems.
-   */
-  su2double GetStandard_RefTemperatureND(void) const { return Standard_Ref_Temperature / Temperature_Ref; }
-
   /*!
    * \brief Get the value of the non-dimensionalized thermodynamic pressure.
    * \return Non-dimensionalized thermodynamic pressure.

Common/include/option_structure.hpp

@@ -96,7 +96,7 @@ const su2double UNIVERSAL_GAS_CONSTANT = 8.3144598;   /*!< \brief Universal gas
 const su2double BOLTZMANN_CONSTANT = 1.3806503E-23;   /*!< \brief Boltzmann's constant [J K^-1] */
 const su2double AVOGAD_CONSTANT = 6.0221415E26; /*!< \brief Avogadro's constant, number of particles in one kmole. */
 const su2double FUND_ELEC_CHARGE_CGS = 4.8032047E-10; /*!< \brief Fundamental electric charge in CGS units, cm^(3/2) g^(1/2) s^(-1). */
-
+const su2double STD_REF_TEMP = 298.15;  /*!< \brief Standard reference temperature for enthalpy in Kelvin. */
 const su2double EPS = 1.0E-16;        /*!< \brief Error scale. */
 const su2double TURB_EPS = 1.0E-16;   /*!< \brief Turbulent Error scale. */
 

Common/src/CConfig.cpp

@@ -1241,9 +1241,6 @@ void CConfig::SetConfig_Options() {
   addDoubleOption("GAMMA_VALUE", Gamma, 1.4);
   /*!\brief THERMODYNAMIC_PRESSURE  \n DESCRIPTION: Thermodynamics(operating) Pressure (101325 Pa), only for incompressible flows) \ingroup Config*/
   addDoubleOption("THERMODYNAMIC_PRESSURE", Pressure_Thermodynamic, 101325.0);
-  /*!\brief STANDARD_REFERENCE_TEMPERATURE  \n DESCRIPTION: Standard reference temperature (298.15K), only for
-   * multicomponent incompressible flows) \ingroup Config*/
-  addDoubleOption("STANDARD_REFERENCE_TEMPERATURE", Standard_Ref_Temperature, 298.15);
   /*!\brief CP_VALUE  \n DESCRIPTION: Specific heat at constant pressure, Cp (1004.703 J/kg*K (air), constant density incompressible fluids only) \ingroup Config*/
   addDoubleListOption("SPECIFIC_HEAT_CP", nSpecific_Heat_Cp, Specific_Heat_Cp);
   /*!\brief THERMAL_EXPANSION_COEFF  \n DESCRIPTION: Thermal expansion coefficient (0.00347 K^-1 (air), used for Boussinesq approximation for liquids/non-ideal gases) \ingroup Config*/

SU2_CFD/include/fluid/CConstantDensity.hpp

@@ -39,10 +39,11 @@ class CConstantDensity final : public CFluidModel {
   /*!
    * \brief Constructor of the class.
    */
-  CConstantDensity(su2double val_Density, su2double val_Cp) {
+  CConstantDensity(su2double val_Density, su2double val_Cp, su2double val_Temperature_Ref) {
     Density = val_Density;
     Cp = val_Cp;
     Cv = val_Cp;
+    Std_Ref_Temp_ND = val_Temperature_Ref;
   }
 
   /*!
@@ -56,16 +57,19 @@ class CConstantDensity final : public CFluidModel {
        decoupled equation. Hence, we update the value.
        Note Cp = Cv, (gamma = 1).*/
     Temperature = t;
-    Enthalpy = Cp * Temperature;
+    Enthalpy = Cp * (Temperature - Std_Ref_Temp_ND);  // Sensible enthalpy relative to STD_REF_TEMP
   }
 
   /*!
    * \brief Set the Dimensionless State using Enthalpy.
    * \param[in] val_enthalpy - Enthalpy value at the point.
-   * \param[in] val_scalars - not used here. 
+   * \param[in] val_scalars - not used here.
    */
   void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override {
     Enthalpy = val_enthalpy;
-    Temperature = Enthalpy / Cp;
+    Temperature = Enthalpy / Cp + Std_Ref_Temp_ND;  // Temperature from sensible enthalpy
   }
+
+ private:
+  su2double Std_Ref_Temp_ND{0.0}; /*!< \brief Nondimensional standard reference temperature for enthalpy. */
 };

SU2_CFD/include/fluid/CFluidModel.hpp

@@ -279,7 +279,7 @@ class CFluidModel {
   /*!
    * \brief Set specific heat Cp model.
    */
-  virtual void SetCpModel(const CConfig* config) {}
+  virtual void SetCpModel(const CConfig* config, su2double val_Temperature_Ref) {}
 
   /*!
    * \brief Set viscosity model.
@@ -372,7 +372,7 @@ class CFluidModel {
   /*!
    * \brief Virtual member.
    * \param[in] val_enthalpy - Enthalpy value at the point.
-   * \param[in] val_scalars - Scalar mass fractions. 
+   * \param[in] val_scalars - Scalar mass fractions.
    */
   virtual void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) {}
 

SU2_CFD/include/fluid/CFluidScalar.hpp

@@ -42,8 +42,8 @@ class CFluidScalar final : public CFluidModel {
   const int n_species_mixture;            /*!< \brief Number of species in mixture. */
   su2double Gas_Constant;                 /*!< \brief Specific gas constant. */
   const su2double Pressure_Thermodynamic; /*!< \brief Constant pressure thermodynamic. */
-  const su2double Ref_Temperature;        /*!< \brief Standard Reference temperature, usually set to 298.15 K. */
   const su2double GasConstant_Ref;        /*!< \brief Gas constant reference needed for Nondimensional problems. */
+  const su2double Std_Ref_Temp_ND;        /*!< \brief Nondimensional standard reference temperature for enthalpy. */
   const su2double Prandtl_Turb_Number;    /*!< \brief Prandlt turbulent number.*/
   const su2double Schmidt_Turb_Number;    /*!< \brief Schmidt turbulent number.*/
 
@@ -177,7 +177,7 @@ class CFluidScalar final : public CFluidModel {
   /*!
    * \brief Virtual member.
    * \param[in] val_enthalpy - Enthalpy value at the point.
-   * \param[in] val_scalars - Scalar mass fractions. 
+   * \param[in] val_scalars - Scalar mass fractions.
    */
   void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override;
 };

SU2_CFD/include/fluid/CIncIdealGas.hpp

@@ -39,7 +39,7 @@ class CIncIdealGas final : public CFluidModel {
   /*!
    * \brief Constructor of the class.
    */
-  CIncIdealGas(su2double val_Cp, su2double val_gas_constant, su2double val_operating_pressure) {
+  CIncIdealGas(su2double val_Cp, su2double val_gas_constant, su2double val_operating_pressure, su2double val_Temperature_Ref) {
     /*--- In the incompressible ideal gas model, the thermodynamic pressure
   is decoupled from the governing equations and held constant. The
   density is therefore only a function of temperature variations. ---*/
@@ -48,6 +48,7 @@ class CIncIdealGas final : public CFluidModel {
     Gamma = 1.0;
     Cp = val_Cp;
     Cv = Cp;
+    Std_Ref_Temp_ND = STD_REF_TEMP / val_Temperature_Ref;
   }
 
   /*!
@@ -58,7 +59,7 @@ class CIncIdealGas final : public CFluidModel {
     /*--- The EoS only depends upon temperature. ---*/
     Temperature = t;
     Density = Pressure / (Temperature * Gas_Constant);
-    Enthalpy = Cp * Temperature;
+    Enthalpy = Cp * (Temperature - Std_Ref_Temp_ND);  // Sensible enthalpy relative to REF_TEMP
   }
 
   /*!
@@ -67,11 +68,12 @@ class CIncIdealGas final : public CFluidModel {
    */
   void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override {
     Enthalpy = val_enthalpy;
-    Temperature = Enthalpy / Cp;
+    Temperature = Enthalpy / Cp + Std_Ref_Temp_ND;  // Temperature from sensible enthalpy
     Density = Pressure / (Temperature * Gas_Constant);
   }
 
  private:
   su2double Gas_Constant{0.0}; /*!< \brief Gas Constant. */
   su2double Gamma{0.0};        /*!< \brief Heat Capacity Ratio. */
+  su2double Std_Ref_Temp_ND{0.0}; /*!< \brief Nondimensional standard reference temperature for enthalpy. */
 };

SU2_CFD/include/fluid/CIncIdealGasPolynomial.hpp

@@ -42,7 +42,7 @@ class CIncIdealGasPolynomial final : public CFluidModel {
   /*!
    * \brief Constructor of the class.
    */
-  CIncIdealGasPolynomial(su2double val_gas_constant, su2double val_operating_pressure) {
+  CIncIdealGasPolynomial(su2double val_gas_constant, su2double val_operating_pressure, su2double val_Temperature_Ref) {
     /* In the incompressible ideal gas model, the thermodynamic pressure is decoupled
     from the governing equations and held constant. The density is therefore only a
     function of temperature variations. We also use a molecular weight (g/mol) and the
@@ -51,20 +51,17 @@ class CIncIdealGasPolynomial final : public CFluidModel {
     Gas_Constant = val_gas_constant;
     Pressure = val_operating_pressure;
     Gamma = 1.0;
+    Std_Ref_Temp_ND = val_Temperature_Ref;
   }
 
   /*!
    * \brief Set the temperature polynomial coefficients for variable Cp.
    * \param[in] config - configuration container for the problem.
    */
-  void SetCpModel(const CConfig* config) override {
-    const su2double t_ref = config->GetStandard_RefTemperatureND();
-    Enthalpy_Ref = 0.0;
-    su2double t_i = 1.0;
+  void SetCpModel(const CConfig* config, su2double val_Temperature_Ref) override {
     for (int i = 0; i < N; ++i) {
-      t_i *= t_ref;
+      /*--- Note that we use cp = dh/dt here. ---*/
       coeffs_[i] = config->GetCp_PolyCoeffND(i);
-      Enthalpy_Ref += coeffs_[i] * t_i / (i + 1);
     }
     Temperature_Min = config->GetTemperatureLimits(0);
   }
@@ -80,12 +77,14 @@ class CIncIdealGasPolynomial final : public CFluidModel {
 
     /* Evaluate the new Cp and enthalpy from the coefficients and temperature. */
     Cp = coeffs_[0];
-    Enthalpy = coeffs_[0] * t - Enthalpy_Ref;
+    Enthalpy = coeffs_[0] * (t - Std_Ref_Temp_ND);
     su2double t_i = 1.0;
+    su2double tref_i = 1.0;
     for (int i = 1; i < N; ++i) {
       t_i *= t;
+      tref_i *= Std_Ref_Temp_ND;
       Cp += coeffs_[i] * t_i;
-      Enthalpy += coeffs_[i] * t_i * t / (i + 1);
+      Enthalpy += coeffs_[i] * (t_i * t - tref_i * Std_Ref_Temp_ND) / (i + 1);
     }
     Cv = Cp / Gamma;
   }
@@ -106,16 +105,18 @@ class CIncIdealGasPolynomial final : public CFluidModel {
 
     int counter = 0;
 
-    /*--- Computing temperature given enthalpy using Newton-Raphson. ---*/
+    /*--- Compute temperature given enthalpy using Newton-Raphson. ---*/
     while ((abs(delta_temp_iter) > toll) && (counter++ < counter_limit)) {
       /* Evaluate the new Cp and enthalpy from the coefficients and temperature. */
       Cp_iter = coeffs_[0];
-      su2double Enthalpy_iter = coeffs_[0] * temp_iter - Enthalpy_Ref;
+      su2double Enthalpy_iter = coeffs_[0] * (temp_iter - Std_Ref_Temp_ND);
       su2double t_i = 1.0;
+      su2double tref_i = 1.0;
       for (int i = 1; i < N; ++i) {
         t_i *= temp_iter;
+        tref_i *= Std_Ref_Temp_ND;
         Cp_iter += coeffs_[i] * t_i;
-        Enthalpy_iter += coeffs_[i] * t_i * temp_iter / (i + 1);
+        Enthalpy_iter += coeffs_[i] * (t_i * temp_iter - tref_i * Std_Ref_Temp_ND) / (i + 1);
       }
 
       delta_enthalpy_iter = Enthalpy - Enthalpy_iter;
@@ -124,15 +125,16 @@ class CIncIdealGasPolynomial final : public CFluidModel {
 
       temp_iter += delta_temp_iter;
       if (temp_iter < Temperature_Min) {
-        cout << "Warning: Negative temperature has been found during Newton-Raphson" << endl;
+        cout << "Warning: Negative temperature has been found during Newton-Raphson." << endl;
         temp_iter = Temperature_Min;
         break;
       }
     }
+
     Temperature = temp_iter;
     Cp = Cp_iter;
     if (counter == counter_limit) {
-      cout << "Warning Newton-Raphson exceed number of max iteration in temperature computation" << endl;
+      cout << "Warning: Newton-Raphson exceeds max. iterations in temperature computation." << endl;
     }
     Density = Pressure / (Temperature * Gas_Constant);
     Cv = Cp / Gamma;
@@ -141,7 +143,7 @@ class CIncIdealGasPolynomial final : public CFluidModel {
  private:
   su2double Gas_Constant{0.0}; /*!< \brief Specific Gas Constant. */
   su2double Gamma{0.0};        /*!< \brief Ratio of specific heats. */
+  su2double Std_Ref_Temp_ND{0.0}; /*!< \brief Nondimensional standard reference temperature for enthalpy. */
   array<su2double, N> coeffs_; /*!< \brief Polynomial coefficients for heat capacity as a function of temperature. */
-  su2double Enthalpy_Ref;      /*!< \brief Enthalpy computed at the reference temperature. */
   su2double Temperature_Min;   /*!< \brief Minimum temperature value allowed in Newton-Raphson iterations. */
 };

SU2_CFD/include/variables/CIncEulerVariable.hpp

@@ -72,7 +72,6 @@ class CIncEulerVariable : public CFlowVariable {
   VectorType Streamwise_Periodic_RecoveredPressure,    /*!< \brief Recovered/Physical pressure [Pa] for streamwise periodic flow. */
              Streamwise_Periodic_RecoveredTemperature; /*!< \brief Recovered/Physical temperature [K] for streamwise periodic flow. */
   su2double TemperatureLimits[2];                      /*!< \brief Temperature limits [K]. */
-  su2double TemperatureInc = 0.0;                      /*!< \brief Temperature [K] imposed when energy equation is switch off. */
  public:
   /*!
    * \brief Constructor of the class.

SU2_CFD/src/fluid/CFluidScalar.cpp

@@ -45,8 +45,8 @@ CFluidScalar::CFluidScalar(su2double value_pressure_operating, const CConfig* co
     : CFluidModel(),
       n_species_mixture(config->GetnSpecies() + 1),
       Pressure_Thermodynamic(value_pressure_operating),
-      Ref_Temperature(config->GetStandard_RefTemperatureND()),
       GasConstant_Ref(config->GetGas_Constant_Ref()),
+      Std_Ref_Temp_ND(STD_REF_TEMP / config->GetTemperature_Ref()),
       Prandtl_Turb_Number(config->GetPrandtl_Turb()),
       Schmidt_Turb_Number(config->GetSchmidt_Number_Turbulent()),
       wilke(config->GetKind_MixingViscosityModel() == MIXINGVISCOSITYMODEL::WILKE),
@@ -219,14 +219,14 @@ su2double CFluidScalar::ComputeEnthalpyFromT(const su2double val_temperature, co
    * depend on temperature, but it does depend on mixture composition, enthalpy is directly computed from
    * the expression h_s = Cp(T - T_ref).
    */
-  su2double val_Enthalpy = Cp * (val_temperature - Ref_Temperature);
+  su2double val_Enthalpy = Cp * (val_temperature - Std_Ref_Temp_ND);
   return val_Enthalpy;
 }
 
 void CFluidScalar::GetEnthalpyDiffusivity(su2double* enthalpy_diffusions) const {
-  const su2double enthalpy_species_N = specificHeat[n_species_mixture - 1] * (Temperature - Ref_Temperature);
+  const su2double enthalpy_species_N = specificHeat[n_species_mixture - 1] * (Temperature - Std_Ref_Temp_ND);
   for (int iVar = 0; iVar < n_species_mixture - 1; iVar++) {
-    const su2double enthalpy_species_i = specificHeat[iVar] * (Temperature - Ref_Temperature);
+    const su2double enthalpy_species_i = specificHeat[iVar] * (Temperature - Std_Ref_Temp_ND);
     enthalpy_diffusions[iVar] = Density * (enthalpy_species_i * massDiffusivity[iVar] -
                                            enthalpy_species_N * massDiffusivity[n_species_mixture - 1]);
     enthalpy_diffusions[iVar] += Mu_Turb * (enthalpy_species_i - enthalpy_species_N) / Schmidt_Turb_Number;
@@ -271,7 +271,7 @@ void CFluidScalar::SetTDState_h(const su2double val_enthalpy, const su2double* v
    * depend on temperature, but it does depend on mixture composition, temperature is directly solved from the
    * expression h_s = Cp(T - T_ref).
    */
-  Temperature = val_enthalpy / Cp + Ref_Temperature;
+  Temperature = val_enthalpy / Cp + Std_Ref_Temp_ND;
   Density = Pressure_Thermodynamic / (Temperature * Gas_Constant);
   Cv = Cp - Gas_Constant;
 
@@ -282,4 +282,5 @@ void CFluidScalar::SetTDState_h(const su2double val_enthalpy, const su2double* v
   }
 
   Kt = WilkeConductivity(val_scalars);
+  ComputeMassDiffusivity();
 }

SU2_CFD/src/output/CFlowIncOutput.cpp

@@ -330,7 +330,8 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){
     AddVolumeOutput("LAMINAR_VISCOSITY", "Laminar_Viscosity", "PRIMITIVE", "Laminar viscosity");
     AddVolumeOutput("HEAT_CAPACITY", "Heat_Capacity", "PRIMITIVE", "Heat capacity");
     AddVolumeOutput("THERMAL_CONDUCTIVITY", "Thermal_Conductivity", "PRIMITIVE", "Thermal conductivity");
-    if (heat || flamelet) AddVolumeOutput("TEMPERATURE", "Temperature", "PRIMITIVE", "Temperature");
+    if (!weakly_coupled_heat)
+      AddVolumeOutput("TEMPERATURE", "Temperature", "PRIMITIVE", "Temperature");
 
     AddVolumeOutput("SKIN_FRICTION-X", "Skin_Friction_Coefficient_x", "PRIMITIVE", "x-component of the skin friction vector");
     AddVolumeOutput("SKIN_FRICTION-Y", "Skin_Friction_Coefficient_y", "PRIMITIVE", "y-component of the skin friction vector");
@@ -354,7 +355,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){
   AddVolumeOutput("RES_VELOCITY-Y", "Residual_Velocity_y", "RESIDUAL", "Residual of the y-velocity component");
   if (nDim == 3)
     AddVolumeOutput("RES_VELOCITY-Z", "Residual_Velocity_z", "RESIDUAL", "Residual of the z-velocity component");
-  if (config->GetEnergy_Equation()){
+  if (heat){
     AddVolumeOutput("RES_ENTHALPY", "Residual_Enthalpy", "RESIDUAL", "Residual of the enthalpy");
   }
   SetVolumeOutputFieldsScalarResidual(config);
@@ -440,7 +441,8 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve
     SetVolumeOutputValue("LAMINAR_VISCOSITY", iPoint, Node_Flow->GetLaminarViscosity(iPoint));
     SetVolumeOutputValue("HEAT_CAPACITY", iPoint, Node_Flow->GetSpecificHeatCp(iPoint));
     SetVolumeOutputValue("THERMAL_CONDUCTIVITY", iPoint, Node_Flow->GetThermalConductivity(iPoint));
-    if (heat || flamelet) SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Flow->GetTemperature(iPoint));
+    if (!weakly_coupled_heat)
+      SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Flow->GetTemperature(iPoint));
   }
 
   SetVolumeOutputValue("RES_PRESSURE", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, 0));