[WIP] Preconditioning for multicomponent flows

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,16 @@ class CFluidModel {
     return mass_diffusivity;
   }
 
+  /*!
+   * \brief Get heat diffusivity terms.
+   */
+  virtual void GetEnthalpyDiffusivity(su2double* enthalpy_diffusions = nullptr) {}
+
+  /*!
+   * \brief Get gradient heat diffusivity terms.
+   */
+  virtual void GetGradEnthalpyDiffusivity(su2double* grad_enthalpy_diffusions = nullptr) {}
+
   /*!
    * \brief Get fluid pressure partial derivative.
    */
@@ -339,6 +355,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.
+   */
+  virtual void ComputeTempFromEnthalpy(su2double val_enthalpy, su2double* val_temperature,
+                                       const su2double* val_scalars = nullptr) {}
+
   /*!
    * \brief Set fluid eddy viscosity provided by a turbulence model needed for computing effective thermal conductivity.
    */

SU2_CFD/include/fluid/CFluidScalar.hpp

@@ -91,6 +91,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.
    */
@@ -137,6 +142,22 @@ class CFluidScalar final : public CFluidModel {
    */
   inline su2double GetMassDiffusivity(int ivar) override { return massDiffusivity[ivar]; }
 
+  /*!
+   * \brief Get enthalpy diffusivity terms.
+   */
+  void GetEnthalpyDiffusivity(su2double* enthalpy_diffusions) override;
+
+  /*!
+   * \brief Get gradient enthalpy diffusivity terms.
+   */
+  void GetGradEnthalpyDiffusivity(su2double* grad_enthalpy_diffusions) override;
+
+  /*!
+   * \brief Compute Temperature from Enthalpy and scalars.
+   */
+  void ComputeTempFromEnthalpy(const su2double val_enthalpy, su2double* val_temperature,
+                               const su2double* val_scalars) override;
+
   /*!
    * \brief Set the Dimensionless State using Temperature.
    * \param[in] t - Temperature value at the point.

SU2_CFD/include/numerics/CNumerics.hpp

@@ -130,6 +130,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
+  Jac_HeatFluxDiffusion;   /*!< \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. */
@@ -188,6 +192,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:
@@ -751,6 +757,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_jac_heatfluxdiffusion) {
+    Jac_HeatFluxDiffusion = val_jac_heatfluxdiffusion;
+  }
+
   /*!
    * \brief Set the laminar viscosity.
    * \param[in] val_laminar_viscosity_i - Value of the laminar viscosity at point i.

SU2_CFD/include/numerics/flow/flow_diffusion.hpp

@@ -51,6 +51,8 @@ class CAvgGrad_Base : public CNumerics {
   *heat_flux_jac_i = nullptr,             /*!< \brief Jacobian of the molecular + turbulent heat flux vector, projected onto the normal vector. */
   **tau_jacobian_i = nullptr;             /*!< \brief Jacobian of the viscous + turbulent stress tensor, projected onto the normal vector. */
   su2double *Mean_PrimVar = nullptr;      /*!< \brief Mean primitive variables. */
+  su2double Mean_HeatFluxDiffusion;       /*!< \brief Mean heat flux due to enthalpy diffusion for multicomponent flows. */
+  su2double Mean_JacHeatFluxDiffusion;    /*!< \brief Mean Jacobian heat flux due to enthalpy diffusion for multicomponent flows. */
   const su2double
   *PrimVar_i = nullptr,
   *PrimVar_j = nullptr;                   /*!< \brief Primitives variables at point i and j. */
@@ -287,6 +289,7 @@ class CAvgGrad_Flow final : public CAvgGrad_Base {
 class CAvgGradInc_Flow final : public CAvgGrad_Base {
 private:
   su2double Mean_Thermal_Conductivity; /*!< \brief Mean value of the effective thermal conductivity. */
+  su2double Mean_Heat_Capacity;        /*!< \brief Mean value of the heat capacity. */
   bool energy;                         /*!< \brief computation with the energy equation. */
 
   /*!
@@ -298,10 +301,10 @@ class CAvgGradInc_Flow final : public CAvgGrad_Base {
    * \param[in] val_gradprimvar - Gradient of the primitive variables.
    * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face.
    * \param[in] val_thermal_conductivity - Thermal conductivity.
+   * \param[in] val_heatDiffusion - Heat diffusion
    */
-  void GetViscousIncProjFlux(const su2double* const *val_gradprimvar,
-                             const su2double *val_normal,
-                             su2double val_thermal_conductivity);
+  void GetViscousIncProjFlux(const su2double* const* val_gradprimvar, const su2double* val_normal,
+                             su2double val_thermal_conductivity, su2double val_heatFluxDiffusion);
 
   /*!
    * \brief Compute the projection of the viscous Jacobian matrices.

SU2_CFD/include/numerics/flow/flow_sources.hpp

@@ -126,7 +126,8 @@ class CSourceGeneralAxisymmetric_Flow final : public CSourceAxisymmetric_Flow {
 class CSourceIncAxisymmetric_Flow final : public CSourceBase_Flow {
   bool implicit, /*!< \brief Implicit calculation. */
   viscous,       /*!< \brief Viscous incompressible flows. */
-  energy;        /*!< \brief computation with the energy equation. */
+  energy,        /*!< \brief computation with the energy equation. */
+  multicomponent; /*!< \brief multicomponent incompressible flows. */
 
 public:
   /*!

SU2_CFD/include/output/CFlowIncOutput.hpp

@@ -40,6 +40,7 @@ class CFlowIncOutput final: public CFlowOutput {
 private:
   TURB_MODEL turb_model;     /*!< \brief The kind of turbulence model*/
   bool heat;                 /*!< \brief Boolean indicating whether have a heat problem*/
+  bool multicomponent;       /*!< \brief Boolean indicating whether have a multicomponent problem*/
   bool weakly_coupled_heat;  /*!< \brief Boolean indicating whether have a weakly coupled heat equation*/
   bool flamelet;  /*!< \brief Boolean indicating whether we solve the flamelet equations */
   unsigned short streamwisePeriodic;   /*!< \brief Boolean indicating whether it is a streamwise periodic simulation. */

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -69,7 +69,8 @@ class CFVMFlowSolverBase : public CSolver {
   su2double Mach_Inf = 0.0;          /*!< \brief Mach number at the infinity. */
   su2double Density_Inf = 0.0;       /*!< \brief Density at the infinity. */
   su2double Energy_Inf = 0.0;        /*!< \brief Energy at the infinity. */
-  su2double Temperature_Inf = 0.0;   /*!< \brief Energy at the infinity. */
+  su2double Temperature_Inf = 0.0;   /*!< \brief Temperature at the infinity. */
+  su2double Enthalpy_Inf = 0.0;      /*!< \brief Enthalpy at the infinity. */
   su2double Pressure_Inf = 0.0;      /*!< \brief Pressure at the infinity. */
   su2double* Velocity_Inf = nullptr; /*!< \brief Flow Velocity vector at the infinity. */
 

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -434,7 +434,7 @@ void CFVMFlowSolverBase<V, R>::Viscous_Residual_impl(unsigned long iEdge, CGeome
   const bool implicit  = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
   const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST);
 
-  CVariable* turbNodes = nullptr;
+  CVariable* turbNodes = nullptr; 
   if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes();
 
   /*--- Points, coordinates and normal vector in edge ---*/

SU2_CFD/include/solvers/CIncEulerSolver.hpp

@@ -191,6 +191,19 @@ class CIncEulerSolver : public CFVMFlowSolverBase<CIncEulerVariable, ENUM_REGIME
                       CConfig *config,
                       unsigned short iMesh) final;
 
+  /*!
+   * \brief Recompute the extrapolated quantities, after MUSCL reconstruction,
+   *        in a more thermodynamically consistent way.
+   * \note This method is static to improve the chances of it being used in a
+   *       thread-safe manner.
+   * \param[in,out] fluidModel - The fluid model.
+   * \param[in] nDim - Number of physical dimensions.
+   * \param[in,out] primitive - Primitive variables.
+   * \param[in] scalar - scalar variable.
+   */
+  static void ComputeConsistentExtrapolation(CFluidModel* fluidModel, unsigned short nDim, su2double* primitive,
+                                             const su2double* scalar);
+
   /*!
    * \brief Source term integration.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/include/variables/CIncEulerVariable.hpp

@@ -40,7 +40,7 @@
  */
 class CIncEulerVariable : public CFlowVariable {
 public:
-  static constexpr size_t MAXNVAR = 12;
+  static constexpr size_t MAXNVAR = 13;
 
   template <class IndexType>
   struct CIndices {
@@ -59,11 +59,11 @@ class CIncEulerVariable : public CFlowVariable {
     inline IndexType ThermalConductivity() const { return nDim+6; }
     inline IndexType CpTotal() const { return nDim+7; }
     inline IndexType CvTotal() const { return nDim+8; }
+    inline IndexType Enthalpy() const { return nDim + 9; }
 
     /*--- For compatible interface with NEMO. ---*/
     inline IndexType SpeciesDensities() const { return std::numeric_limits<IndexType>::max(); }
     inline IndexType Temperature_ve() const { return std::numeric_limits<IndexType>::max(); }
-    inline IndexType Enthalpy() const { return std::numeric_limits<IndexType>::max(); }
   };
 
  protected:
@@ -121,6 +121,14 @@ class CIncEulerVariable : public CFlowVariable {
     return val_temperature <= 0.0;
   }
 
+  /*!
+   * \brief Set the value of the enthalpy for multicomponent incompressible flows with energy equation.
+   * \param[in] iPoint - Point index.
+   */
+  inline void SetEnthalpy(unsigned long iPoint, su2double val_enthalpy) {
+    Primitive(iPoint, indices.Enthalpy()) = val_enthalpy;
+  }
+
   /*!
    * \brief Set the value of the beta coeffient for incompressible flows.
    * \param[in] iPoint - Point index.
@@ -153,6 +161,12 @@ class CIncEulerVariable : public CFlowVariable {
    */
   inline su2double GetTemperature(unsigned long iPoint) const final { return Primitive(iPoint, indices.Temperature()); }
 
+  /*!
+   * \brief Get the enthalpy of the flow.
+   * \return Value of the enthalpy of the flow.
+   */
+  inline su2double GetEnthalpy(unsigned long iPoint) const final { return Primitive(iPoint, indices.Enthalpy()); }
+
   /*!
    * \brief Get the velocity of the flow.
    * \param[in] iDim - Index of the dimension.

SU2_CFD/include/variables/CIncNSVariable.hpp

@@ -40,6 +40,7 @@ class CIncNSVariable final : public CIncEulerVariable {
 private:
   VectorType Tau_Wall;        /*!< \brief Magnitude of the wall shear stress from a wall function. */
   VectorType DES_LengthScale;
+  bool Energy_Multicomponent = false;
 
 public:
   /*!

SU2_CFD/src/fluid/CFluidScalar.cpp

@@ -212,13 +212,43 @@ su2double CFluidScalar::ComputeMeanSpecificHeatCp(const su2double* val_scalars)
   return mean_cp;
 }
 
+su2double CFluidScalar::ComputeEnthalpyFromT(const su2double val_temperature, const su2double* val_scalars){
+  su2double val_Enthalpy = Cp * (val_temperature - 298.15);
+  return val_Enthalpy;
+}
+
+void CFluidScalar::ComputeTempFromEnthalpy(const su2double val_enthalpy, su2double* val_temperature,
+                                           const su2double* val_scalars) {
+  MassToMoleFractions(val_scalars);
+  su2double val_cp = ComputeMeanSpecificHeatCp(val_scalars);
+  *val_temperature = val_enthalpy / val_cp + 298.15;
+}
+
+void CFluidScalar::GetEnthalpyDiffusivity(su2double* enthalpy_diffusions) {
+  for (int iVar = 0; iVar < n_species_mixture - 1; iVar++) {
+    enthalpy_diffusions[iVar] = Density *
+                               (specificHeat[iVar] * massDiffusivity[iVar] -
+                                specificHeat[n_species_mixture - 1] * massDiffusivity[n_species_mixture - 1]) *
+                               (Temperature - 298.15);
+  }
+}
+
+void CFluidScalar::GetGradEnthalpyDiffusivity(su2double* grad_enthalpy_diffusions){
+  for (int iVar = 0; iVar < n_species_mixture - 1; iVar++) {
+    grad_enthalpy_diffusions[iVar] = Density *
+                               (specificHeat[iVar] * massDiffusivity[iVar] -
+                                specificHeat[n_species_mixture - 1] * massDiffusivity[n_species_mixture - 1]);
+  }
+}
+
 void CFluidScalar::SetTDState_T(const su2double val_temperature, const su2double* val_scalars) {
   MassToMoleFractions(val_scalars);
   ComputeGasConstant();
   Temperature = val_temperature;
   Density = Pressure_Thermodynamic / (Temperature * Gas_Constant);
   Cp = ComputeMeanSpecificHeatCp(val_scalars);
   Cv = Cp - Gas_Constant;
+  Enthalpy = ComputeEnthalpyFromT(Temperature, val_scalars);
 
   if (wilke) {
     Mu = WilkeViscosity(val_scalars);