su2code · Cristopher-Morales · Dec 13, 2025 · Oct 22, 2024 · Oct 25, 2024 · Oct 26, 2024
diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
@@ -919,6 +919,7 @@ 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 */
@@ -1958,6 +1959,12 @@ 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.

diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
@@ -1229,6 +1229,9 @@ 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*/

diff --git a/SU2_CFD/include/fluid/CConstantDensity.hpp b/SU2_CFD/include/fluid/CConstantDensity.hpp
@@ -56,5 +56,16 @@ class CConstantDensity final : public CFluidModel {
        decoupled equation. Hence, we update the value.
        Note Cp = Cv, (gamma = 1).*/
     Temperature = t;
+    Enthalpy = Cp * Temperature;
+  }
+
+  /*!
+   * \brief Set the Dimensionless State using Enthalpy.
+   * \param[in] val_enthalpy - Enthalpy value at the point.
+   * \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;
   }
 };
diff --git a/SU2_CFD/include/fluid/CFluidFlamelet.hpp b/SU2_CFD/include/fluid/CFluidFlamelet.hpp
@@ -118,6 +118,13 @@ class CFluidFlamelet final : public CFluidModel {
    */
   void SetTDState_T(su2double val_temperature, const su2double* val_scalars = nullptr) override;
 
+  /*!
+   * \brief Set the thermodynamic state.
+   * \param[in] val_enthalpy - enthalpy
+   * \param[in] val_scalars - pointer to species mass fractions
+   */
+  void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override;
+
   /*!
    * \brief Evaluate data-set for flamelet simulations.
    * \param[in] input_scalar - controlling variables used to interpolate manifold.

diff --git a/SU2_CFD/include/fluid/CFluidModel.hpp b/SU2_CFD/include/fluid/CFluidModel.hpp
@@ -47,6 +47,7 @@ class CLookUpTable;
 class CFluidModel {
  protected:
   su2double StaticEnergy{0.0};     /*!< \brief Internal Energy. */
+  su2double Enthalpy{0.0};         /*!<*\brief Enthalpy. */
   su2double Entropy{0.0};          /*!< \brief Entropy. */
   su2double Density{0.0};          /*!< \brief Density. */
   su2double Pressure{0.0};         /*!< \brief Pressure. */
@@ -113,6 +114,11 @@ class CFluidModel {
    */
   su2double GetStaticEnergy() const { return StaticEnergy; }
 
+  /*!
+   * \brief Get fluid enthalpy.
+   */
+  su2double GetEnthalpy() const { return Enthalpy; }
+
   /*!
    * \brief Get fluid density.
    */
@@ -186,6 +192,30 @@ class CFluidModel {
     return mass_diffusivity;
   }
 
+  /*!
+   * \brief Get the enthalpy diffusivity terms for all species being solved.
+   *
+   * This function computes and retrieves the enthalpy diffusion terms required in the energy equation
+   * for multicomponent flows.
+   *
+   * \param[in,out] enthalpy_diffusions - Array containing the enthalpy diffusion terms for all
+   * species to be solved. The size of \p enthalpy_diffusions must be at least (n_species_mixture - 1),
+   * corresponding to the number of species transport equations in the system.
+   */
+  virtual void GetEnthalpyDiffusivity(su2double* enthalpy_diffusions = nullptr) const {}
+
+  /*!
+   * \brief Get the gradient of enthalpy diffusivity terms for all species being solved.
+   *
+   * This function computes and retrieves the gradient of the enthalpy diffusion terms with respect to temperature.
+   * These terms are required for implicit computations when solving the energy equation for multicomponent flows.
+   *
+   * \param[in,out] grad_enthalpy_diffusions - Array containing the gradient of enthalpy diffusion terms for all
+   * species to be solved. The size of \p grad_enthalpy_diffusions must be at least (n_species_mixture - 1),
+   * corresponding to the number of species transport equations in the system.
+   */
+  virtual void GetGradEnthalpyDiffusivity(su2double* grad_enthalpy_diffusions = nullptr) const {}
+
   /*!
    * \brief Get fluid pressure partial derivative.
    */
@@ -339,6 +369,13 @@ class CFluidModel {
    */
   virtual void SetTDState_T(su2double val_Temperature, const su2double* val_scalars = nullptr) {}
 
+  /*!
+   * \brief Virtual member.
+   * \param[in] val_enthalpy - Enthalpy value at the point.
+   * \param[in] val_scalars - Scalar mass fractions. 
+   */
+  virtual void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) {}
+
   /*!
    * \brief Set fluid eddy viscosity provided by a turbulence model needed for computing effective thermal conductivity.
    */

diff --git a/SU2_CFD/include/fluid/CFluidScalar.hpp b/SU2_CFD/include/fluid/CFluidScalar.hpp
@@ -40,10 +40,12 @@
 class CFluidScalar final : public CFluidModel {
  private:
   const int n_species_mixture;            /*!< \brief Number of species in mixture. */
-  su2double Gas_Constant;           /*!< \brief Specific gas constant. */
+  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 Prandtl_Number;         /*!< \brief Prandlt number.*/
+  const su2double Prandtl_Turb_Number;    /*!< \brief Prandlt turbulent number.*/
+  const su2double Schmidt_Turb_Number;    /*!< \brief Schmidt turbulent number.*/
 
   const bool wilke;
   const bool davidson;
@@ -91,6 +93,11 @@ class CFluidScalar final : public CFluidModel {
    */
   su2double ComputeMeanSpecificHeatCp(const su2double* val_scalars);
 
+  /*!
+   * \brief Compute Enthalpy given the temperature and scalars.
+   */
+  su2double ComputeEnthalpyFromT(const su2double val_temperature, const su2double* val_scalars);
+
   /*!
    * \brief Compute gas constant for mixture.
    */
@@ -130,16 +137,47 @@ class CFluidScalar final : public CFluidModel {
   /*!
    * \brief Get fluid thermal conductivity.
    */
-  inline su2double GetThermalConductivity() override { return Kt + Mu_Turb * Cp / Prandtl_Number; }
+  inline su2double GetThermalConductivity() override { return Kt + Mu_Turb * Cp / Prandtl_Turb_Number; }
 
   /*!
    * \brief Get fluid mass diffusivity.
    */
   inline su2double GetMassDiffusivity(int ivar) override { return massDiffusivity[ivar]; }
 
+  /*!
+   * \brief Get the enthalpy diffusivity terms for all species being solved.
+   *
+   * This function computes and retrieves the enthalpy diffusion terms required in the energy equation
+   * for multicomponent flows.
+   *
+   * \param[in,out] enthalpy_diffusions - Array containing the enthalpy diffusion terms for all
+   * species to be solved. The size of \p enthalpy_diffusions must be at least (n_species_mixture - 1),
+   * corresponding to the number of species transport equations in the system.
+   */
+  void GetEnthalpyDiffusivity(su2double* enthalpy_diffusions) const override;
+
+  /*!
+   * \brief Get the gradient of enthalpy diffusivity terms for all species being solved.
+   *
+   * This function computes and retrieves the gradient of the enthalpy diffusion terms with respect to temperature.
+   * These terms are required for implicit computations when solving the energy equation for multicomponent flows.
+   *
+   * \param[in,out] grad_enthalpy_diffusions - Array containing the gradient of enthalpy diffusion terms for all
+   * species to be solved. The size of \p grad_enthalpy_diffusions must be at least (n_species_mixture - 1),
+   * corresponding to the number of species transport equations in the system.
+   */
+  void GetGradEnthalpyDiffusivity(su2double* grad_enthalpy_diffusions) const override;
+
   /*!
    * \brief Set the Dimensionless State using Temperature.
    * \param[in] t - Temperature value at the point.
    */
   void SetTDState_T(su2double val_temperature, const su2double* val_scalars) override;
+
+  /*!
+   * \brief Virtual member.
+   * \param[in] val_enthalpy - Enthalpy value at the point.
+   * \param[in] val_scalars - Scalar mass fractions. 
+   */
+  void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override;
 };
diff --git a/SU2_CFD/include/fluid/CIncIdealGas.hpp b/SU2_CFD/include/fluid/CIncIdealGas.hpp
@@ -58,6 +58,17 @@ class CIncIdealGas final : public CFluidModel {
     /*--- The EoS only depends upon temperature. ---*/
     Temperature = t;
     Density = Pressure / (Temperature * Gas_Constant);
+    Enthalpy = Cp * Temperature;
+  }
+
+  /*!
+   * \brief Set the Dimensionless State using Enthalpy.
+   * \param[in] val_enthalpy - Enthalpy value at the point.
+   */
+  void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override {
+    Enthalpy = val_enthalpy;
+    Temperature = Enthalpy / Cp;
+    Density = Pressure / (Temperature * Gas_Constant);
   }
 
  private:

diff --git a/SU2_CFD/include/fluid/CIncIdealGasPolynomial.hpp b/SU2_CFD/include/fluid/CIncIdealGasPolynomial.hpp
@@ -58,9 +58,15 @@ class CIncIdealGasPolynomial final : public CFluidModel {
    * \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;
     for (int i = 0; i < N; ++i) {
+      t_i *= t_ref;
       coeffs_[i] = config->GetCp_PolyCoeffND(i);
+      Enthalpy_Ref += coeffs_[i] * t_i / (i + 1);
     }
+    Temperature_Min = config->GetTemperatureLimits(0);
   }
 
   /*!
@@ -72,18 +78,70 @@ class CIncIdealGasPolynomial final : public CFluidModel {
     Temperature = t;
     Density = Pressure / (Temperature * Gas_Constant);
 
-    /* Evaluate the new Cp from the coefficients and temperature. */
+    /* Evaluate the new Cp and enthalpy from the coefficients and temperature. */
     Cp = coeffs_[0];
+    Enthalpy = coeffs_[0] * t - Enthalpy_Ref;
     su2double t_i = 1.0;
     for (int i = 1; i < N; ++i) {
       t_i *= t;
       Cp += coeffs_[i] * t_i;
+      Enthalpy += coeffs_[i] * t_i * t / (i + 1);
+    }
+    Cv = Cp / Gamma;
+  }
+
+  /*!
+   * \brief Set the Dimensionless State using enthalpy.
+   * \param[in] val_enthalpy - Enthalpy value at the point.
+   */
+  void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override {
+    Enthalpy = val_enthalpy;
+    /*--- convergence criterion for temperature in [K], high accuracy needed for restarts. ---*/
+    const su2double toll = 1e-5;
+    su2double temp_iter = 300.0;
+    su2double Cp_iter = 0.0;
+    su2double delta_temp_iter = 1e10;
+    su2double delta_enthalpy_iter;
+    const int counter_limit = 20;
+
+    int counter = 0;
+
+    /*--- Computing 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 t_i = 1.0;
+      for (int i = 1; i < N; ++i) {
+        t_i *= temp_iter;
+        Cp_iter += coeffs_[i] * t_i;
+        Enthalpy_iter += coeffs_[i] * t_i * temp_iter / (i + 1);
+      }
+
+      delta_enthalpy_iter = Enthalpy - Enthalpy_iter;
+
+      delta_temp_iter = delta_enthalpy_iter / Cp_iter;
+
+      temp_iter += delta_temp_iter;
+      if (temp_iter < Temperature_Min) {
+        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;
     }
+    Density = Pressure / (Temperature * Gas_Constant);
     Cv = Cp / Gamma;
   }
 
  private:
   su2double Gas_Constant{0.0}; /*!< \brief Specific Gas Constant. */
   su2double Gamma{0.0};        /*!< \brief Ratio of specific heats. */
-  array<su2double, N> coeffs_; /*!< \brief Polynomial coefficients for conductivity as a function of temperature. */
+  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. */
 };
diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp
@@ -129,6 +129,10 @@ class CNumerics {
   const su2double
   *ScalarVar_i,   /*!< \brief Vector of scalar variables at point i. */
   *ScalarVar_j;   /*!< \brief Vector of scalar variables at point j. */
+  su2double
+  HeatFluxDiffusion;   /*!< \brief Heat flux due to enthalpy diffusion for multicomponent. */
+  su2double
+  JacHeatFluxDiffusion;   /*!< \brief Heat flux jacobian due to enthalpy diffusion for multicomponent. */
   const su2double
   *TransVar_i,  /*!< \brief Vector of turbulent variables at point i. */
   *TransVar_j;  /*!< \brief Vector of turbulent variables at point j. */
@@ -187,6 +191,8 @@ class CNumerics {
 
   bool nemo;                      /*!< \brief Flag for NEMO problems  */
 
+  bool energy_multicomponent = false; /*!< \brief Flag for multicomponent and reacting flow  */
+
   bool bounded_scalar = false;    /*!< \brief Flag for bounded scalar problem */
 
 public:
@@ -750,6 +756,20 @@ class CNumerics {
     Diffusion_Coeff_j = val_diffusioncoeff_j;
   }
 
+  /*!
+   * \brief Set the heat flux due to enthalpy diffusion
+   * \param[in] val_heatfluxdiffusion - Value of the heat flux due to enthalpy diffusion.
+   */
+  inline void SetHeatFluxDiffusion(su2double val_heatfluxdiffusion) { HeatFluxDiffusion = val_heatfluxdiffusion; }
+
+  /*!
+   * \brief Set Jacobian of the heat flux due to enthalpy diffusion
+   * \param[in] val_jacheatfluxdiffusion - Value of the heat flux jacobian due to enthalpy diffusion.
+   */
+  inline void SetJacHeatFluxDiffusion(su2double val_jacheatfluxdiffusion) {
+    JacHeatFluxDiffusion = val_jacheatfluxdiffusion;
+  }
+
   /*!
    * \brief Set the laminar viscosity.
    * \param[in] val_laminar_viscosity_i - Value of the laminar viscosity at point i.
@@ -1038,8 +1058,7 @@ class CNumerics {
    * \param[in] val_density - Value of the density.
    * \param[in] val_velocity - Pointer to the velocity.
    * \param[in] val_betainc2 - Value of the artificial compresibility factor.
-   * \param[in] val_cp - Value of the specific heat at constant pressure.
-   * \param[in] val_temperature - Value of the temperature.
+   * \param[in] val_enthalpy - Value of the enthalpy.
    * \param[in] val_dRhodT - Value of the derivative of density w.r.t. temperature.
    * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face.
    * \param[in] val_scale - Scale of the projection.
@@ -1048,8 +1067,7 @@ class CNumerics {
   void GetInviscidIncProjJac(const su2double *val_density,
                              const su2double *val_velocity,
                              const su2double *val_betainc2,
-                             const su2double *val_cp,
-                             const su2double *val_temperature,
+                             const su2double *val_enthalpy,
                              const su2double *val_dRhodT,
                              const su2double *val_normal,
                              su2double val_scale,
@@ -1060,17 +1078,15 @@ class CNumerics {
    * \param[in] val_density - Value of the density.
    * \param[in] val_velocity - Pointer to the velocity.
    * \param[in] val_betainc2 - Value of the artificial compresibility factor.
-   * \param[in] val_cp - Value of the specific heat at constant pressure.
-   * \param[in] val_temperature - Value of the temperature.
-   * \param[in] val_dRhodT - Value of the derivative of density w.r.t. temperature.
+   * \param[in] val_enthalpy - Value of the enthalpy.
+   * \param[in] val_dRhodh - Value of the derivative of density w.r.t. enthalpy.
    * \param[out] val_Precon - Pointer to the preconditioning matrix.
    */
   void GetPreconditioner(const su2double *val_density,
                          const su2double *val_velocity,
                          const su2double *val_betainc2,
-                         const su2double *val_cp,
-                         const su2double *val_temperature,
-                         const su2double *val_drhodt,
+                         const su2double *val_enthalpy,
+                         const su2double *val_drhodh,
                          su2double **val_Precon) const;
 
   /*!