su2code · Cristopher-Morales · Jun 16, 2022 · Apr 13, 2022 · Apr 14, 2022 · Apr 21, 2022
diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
@@ -802,8 +802,9 @@ class CConfig {
   Beta_Factor,          /*!< \brief Value of the epsilon^2 multiplier for Beta for the incompressible preconditioner. */
   Gas_Constant,         /*!< \brief Specific gas constant. */
   Gas_ConstantND,       /*!< \brief Non-dimensional specific gas constant. */
-  Molecular_Weight,     /*!< \brief Molecular weight of an incompressible ideal gas (g/mol). */
-  Specific_Heat_Cp,           /*!< \brief Specific heat at constant pressure. */
+  *Molecular_Weight;    /*!< \brief Molecular weight of an incompressible ideal gas (g/mol). */
+  unsigned short nMolecular_Weight; /*!< \brief Number of species molecular weights. */
+  su2double Specific_Heat_Cp, /*!< \brief Specific heat at constant pressure. */
   Specific_Heat_CpND,         /*!< \brief Non-dimensional specific heat at constant pressure. */
   Specific_Heat_Cv,           /*!< \brief Specific heat at constant volume. */
   Specific_Heat_CvND,         /*!< \brief Non-dimensional specific heat at constant volume. */
@@ -1164,8 +1165,9 @@ class CConfig {
   unsigned short maxBasisDim,               /*!< \brief Maximum number of POD basis dimensions. */
   rom_save_freq;                            /*!< \brief Frequency of unsteady time steps to save. */
 
+  unsigned short nSpecies,                  /*!< \brief Number of transported species equations (for NEMO and species transport)*/
+
   /* other NEMO configure options*/
-  unsigned short nSpecies,                  /*!< \brief No of species present in flow */
   iWall_Catalytic,
   nWall_Catalytic;                          /*!< \brief No of catalytic walls */
   su2double *Gas_Composition,               /*!< \brief Initial mass fractions of flow [dimensionless] */
@@ -1186,7 +1188,6 @@ class CConfig {
   bool Species_StrongBC;           /*!< \brief Boolean whether strong BC's are used for in- outlet of the species solver. */
   su2double* Species_Init;         /*!< \brief Initial uniform value for scalar transport. */
   unsigned short nSpecies_Init;    /*!< \brief Number of entries of SPECIES_INIT */
-  unsigned short n_species;        /*!< \brief Number of species equations! Not species itself which would be 1 more. */
 
   /*!
    * \brief Set the default values of config options not set in the config file using another config object.
@@ -1608,7 +1609,7 @@ class CConfig {
    * \brief Get the value of the molecular weight for an incompressible ideal gas (g/mol).
    * \return Value of the molecular weight for an incompressible ideal gas (g/mol).
    */
-  su2double GetMolecular_Weight(void) const { return Molecular_Weight; }
+  su2double GetMolecular_Weight(unsigned short val_index = 0) const { return Molecular_Weight[val_index]; }
 
   /*!
    * \brief Get the value of specific heat at constant pressure.

diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
@@ -567,7 +567,8 @@ enum ENUM_FLUIDMODEL {
   INC_IDEAL_GAS = 5,      /*!< \brief Incompressible ideal gas model. */
   INC_IDEAL_GAS_POLY = 6, /*!< \brief Inc. ideal gas, polynomial gas model. */
   MUTATIONPP = 7,         /*!< \brief Mutation++ gas model for nonequilibrium flow. */
-  SU2_NONEQ = 8           /*!< \brief User defined gas model for nonequilibrium flow. */
+  SU2_NONEQ = 8,          /*!< \brief User defined gas model for nonequilibrium flow. */
+  FLUID_MIXTURE = 9,      /*!< \brief Species mixture model. */
 };
 static const MapType<std::string, ENUM_FLUIDMODEL> FluidModel_Map = {
   MakePair("STANDARD_AIR", STANDARD_AIR)
@@ -579,6 +580,7 @@ static const MapType<std::string, ENUM_FLUIDMODEL> FluidModel_Map = {
   MakePair("INC_IDEAL_GAS_POLY", INC_IDEAL_GAS_POLY)
   MakePair("MUTATIONPP", MUTATIONPP)
   MakePair("SU2_NONEQ", SU2_NONEQ)
+  MakePair("FLUID_MIXTURE", FLUID_MIXTURE)
 };
 
 /*!

diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
@@ -1050,15 +1050,14 @@ void CConfig::SetPointersNull(void) {
 
   Time_Ref = 1.0;
 
-  Delta_UnstTime   = 0.0;
+  Delta_UnstTime = 0.0;
   Delta_UnstTimeND = 0.0;
-  Total_UnstTime   = 0.0;
+  Total_UnstTime = 0.0;
   Total_UnstTimeND = 0.0;
 
   Kind_TimeNumScheme = EULER_IMPLICIT;
 
   Gas_Composition = nullptr;
-
 }
 
 void CConfig::SetConfig_Options() {
@@ -1168,7 +1167,7 @@ void CConfig::SetConfig_Options() {
   /*!\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*/
   addDoubleOption("THERMAL_EXPANSION_COEFF", Thermal_Expansion_Coeff, 0.00347);
   /*!\brief MOLECULAR_WEIGHT \n DESCRIPTION: Molecular weight for an incompressible ideal gas (28.96 g/mol (air) default) \ingroup Config*/
-  addDoubleOption("MOLECULAR_WEIGHT", Molecular_Weight, 28.96);
+  addDoubleListOption("MOLECULAR_WEIGHT", nMolecular_Weight, Molecular_Weight);
 
   ///* DESCRIPTION: Specify if Mutation++ library is used */
   /*--- Reading gas model as string or integer depending on TC library used. ---*/
@@ -3308,6 +3307,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
   bool ideal_gas = ((Kind_FluidModel == STANDARD_AIR) ||
                     (Kind_FluidModel == IDEAL_GAS) ||
                     (Kind_FluidModel == INC_IDEAL_GAS) ||
+                    (Kind_FluidModel == FLUID_MIXTURE) ||
                     (Kind_FluidModel == INC_IDEAL_GAS_POLY) ||
                     (Kind_FluidModel == CONSTANT_DENSITY));
   bool noneq_gas = ((Kind_FluidModel == MUTATIONPP) ||
@@ -3741,6 +3741,35 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     }
   }
 
+/*--- Set default values for various fluid properties. ---*/
+
+  const su2double Molecular_Weight_Default = 28.96;
+
+  if (Molecular_Weight == nullptr) {
+    Molecular_Weight = new su2double[1];
+    Molecular_Weight[0] = Molecular_Weight_Default;
+    nMolecular_Weight = 1;
+  }
+
+  /*--- Check whether inputs for FLUID_MIXTURE are correctly specified. ---*/
+
+  if (Kind_FluidModel == FLUID_MIXTURE) {
+    /*--- Check whether the number of entries of each specified fluid property equals the number of transported scalar
+     equations solved + 1. nMolecular_Weight is used because it is required for the fluid mixing models. --- */
+    if (nMolecular_Weight != nSpecies_Init + 1) {
+      SU2_MPI::Error(
+          "The use of FLUID_MIXTURE requires the number of entries for MOLECULAR_WEIGHT\n"
+          "to be equal to the number of entries of SCALAR_INIT + 1",
+          CURRENT_FUNCTION);
+    }
+    /*--- Check whether the density model used is correct, in the case of FLUID_MIXTURE the density model must be
+     VARIABLE. Otherwise, if the density model is CONSTANT, the scalars will not have influence the mixture density and
+     it will remain constant through the complete domain. --- */
+    if (Kind_DensityModel != INC_DENSITYMODEL::VARIABLE) {
+      SU2_MPI::Error("The use of FLUID_MIXTURE requires the INC_DENSITY_MODEL option to be VARIABLE", CURRENT_FUNCTION);
+    }
+  }
+
   /*--- Overrule the default values for viscosity if the US measurement system is used. ---*/
 
   if (SystemMeasurements == US) {
@@ -4725,7 +4754,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
   }
 
   if (Kind_DensityModel == INC_DENSITYMODEL::VARIABLE) {
-    if (Kind_FluidModel != INC_IDEAL_GAS && Kind_FluidModel != INC_IDEAL_GAS_POLY) {
+    if (Kind_FluidModel != INC_IDEAL_GAS && Kind_FluidModel != INC_IDEAL_GAS_POLY && Kind_FluidModel != FLUID_MIXTURE) {
       SU2_MPI::Error("Variable density incompressible solver limited to ideal gases.\n Check the fluid model options (use INC_IDEAL_GAS, INC_IDEAL_GAS_POLY).", CURRENT_FUNCTION);
     }
   }
@@ -4738,7 +4767,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
 
   if (Kind_Solver == MAIN_SOLVER::INC_NAVIER_STOKES || Kind_Solver == MAIN_SOLVER::INC_RANS) {
     if (Kind_ViscosityModel == VISCOSITYMODEL::SUTHERLAND) {
-      if ((Kind_FluidModel != INC_IDEAL_GAS) && (Kind_FluidModel != INC_IDEAL_GAS_POLY)) {
+      if ((Kind_FluidModel != INC_IDEAL_GAS) && (Kind_FluidModel != INC_IDEAL_GAS_POLY) && (Kind_FluidModel != FLUID_MIXTURE)) {
         SU2_MPI::Error("Sutherland's law only valid for ideal gases in incompressible flows.\n Must use VISCOSITY_MODEL=CONSTANT_VISCOSITY and set viscosity with\n MU_CONSTANT, or use DENSITY_MODEL= VARIABLE with FLUID_MODEL= INC_IDEAL_GAS or INC_IDEAL_GAS_POLY for VISCOSITY_MODEL=SUTHERLAND.\n NOTE: FREESTREAM_VISCOSITY is no longer used for incompressible flows!", CURRENT_FUNCTION);
       }
     }

diff --git a/SU2_CFD/include/fluid/CConstantDensity.hpp b/SU2_CFD/include/fluid/CConstantDensity.hpp
@@ -49,7 +49,7 @@ class CConstantDensity final : public CFluidModel {
    * \brief Set the Dimensionless State using Temperature.
    * \param[in] t - Temperature value at the point.
    */
-  void SetTDState_T(su2double t) override {
+  void SetTDState_T(su2double t, const su2double *val_scalars = nullptr) override {
     /* Density is constant and thermodynamic pressure is
        not required for incompressible, constant density flows,
        but the energy equation can still be computed as a

diff --git a/SU2_CFD/include/fluid/CFluidModel.hpp b/SU2_CFD/include/fluid/CFluidModel.hpp
@@ -119,10 +119,26 @@ class CFluidModel {
    */
   su2double GetCv() const { return Cv; }
 
+  /*!
+   * \brief Compute and return fluid mean molecular weight in kg/mol.
+   */
+  template <class Vector_t>
+  static su2double ComputeMeanMolecularWeight(const Vector_t& molar_masses, const su2double* val_scalars) {
+    su2double OneOverMeanMolecularWeight = 0.0;
+    su2double val_scalars_sum = 0.0;
+
+    for (size_t i_scalar = 0; i_scalar < molar_masses.size() - 1; i_scalar++) {
+      OneOverMeanMolecularWeight += val_scalars[i_scalar] / (molar_masses[i_scalar] / 1000);
+      val_scalars_sum += val_scalars[i_scalar];
+    }
+    OneOverMeanMolecularWeight += (1 - val_scalars_sum) / (molar_masses[molar_masses.size() - 1] / 1000);
+    return 1 / OneOverMeanMolecularWeight;
+  }
+
   /*!
    * \brief Get fluid dynamic viscosity.
    */
-  su2double GetLaminarViscosity() {
+  inline su2double GetLaminarViscosity() {
     LaminarViscosity->SetViscosity(Temperature, Density);
     Mu = LaminarViscosity->GetViscosity();
     LaminarViscosity->SetDerViscosity(Temperature, Density);
@@ -135,7 +151,7 @@ class CFluidModel {
    * \brief Get fluid thermal conductivity.
    */
 
-  su2double GetThermalConductivity() {
+  inline su2double GetThermalConductivity() {
     ThermalConductivity->SetConductivity(Temperature, Density, Mu, Mu_Turb, Cp);
     Kt = ThermalConductivity->GetConductivity();
     ThermalConductivity->SetDerConductivity(Temperature, Density, dmudrho_T, dmudT_rho, Cp);
@@ -291,11 +307,10 @@ class CFluidModel {
    * \brief Virtual member.
    * \param[in] T - Temperature value at the point.
    */
-  virtual void SetTDState_T(su2double val_Temperature) {}
+  virtual void SetTDState_T(su2double val_Temperature, const su2double* val_scalars = nullptr) {}
 
   /*!
    * \brief Set fluid eddy viscosity provided by a turbulence model needed for computing effective thermal conductivity.
    */
   void SetEddyViscosity(su2double val_Mu_Turb) { Mu_Turb = val_Mu_Turb; }
-
 };
diff --git a/SU2_CFD/include/fluid/CFluidScalar.hpp b/SU2_CFD/include/fluid/CFluidScalar.hpp
@@ -0,0 +1,59 @@
+/*!
+ * \file CFluidScalar.hpp
+ * \brief  Defines the multicomponent incompressible Ideal Gas model for mixtures.
+ * \author T. Economon, Mark Heimgartner, Cristopher Morales Ubal
+ * \version 7.3.1 "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/>.
+ */
+
+#pragma once
+#include <memory>
+#include <vector>
+
+#include "CFluidModel.hpp"
+
+/*!
+ * \class CFluidScalar
+ * \brief Child class for defining an incompressible ideal gas model.
+ * \author: T. Economon
+ */
+class CFluidScalar final : public CFluidModel {
+ private:
+  unsigned short n_species_mixture; /*!< \brief Number of species in mixture. */
+  su2double Gas_Constant;           /*!< \brief Specific gas constant. */
+  su2double Gamma;                  /*!< \brief Ratio of specific heats of the gas. */
+  su2double Pressure_Thermodynamic; /*!< \brief Constant pressure thermodynamic. */
+
+  std::vector<su2double> molarMasses;  /*!< \brief Molar masses of all species. */
+
+ public:
+  /*!
+   * \brief Constructor of the class.
+   */
+  CFluidScalar(su2double val_Cp, su2double val_gas_constant, su2double val_operating_pressure, CConfig* config);
+
+  /*!
+   * \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;
+};
diff --git a/SU2_CFD/include/fluid/CIncIdealGas.hpp b/SU2_CFD/include/fluid/CIncIdealGas.hpp
@@ -54,7 +54,7 @@ class CIncIdealGas final : public CFluidModel {
    * \brief Set the Dimensionless State using Temperature.
    * \param[in] t - Temperature value at the point.
    */
-  void SetTDState_T(su2double t) override {
+  void SetTDState_T(su2double t, const su2double *val_scalars = nullptr) override {
     /*--- The EoS only depends upon temperature. ---*/
     Temperature = t;
     Density = Pressure / (Temperature * Gas_Constant);

diff --git a/SU2_CFD/include/fluid/CIncIdealGasPolynomial.hpp b/SU2_CFD/include/fluid/CIncIdealGasPolynomial.hpp
@@ -67,7 +67,7 @@ class CIncIdealGasPolynomial final : public CFluidModel {
    * \brief Set the Dimensionless State using Temperature.
    * \param[in] t - Temperature value at the point.
    */
-  void SetTDState_T(su2double t) override {
+  void SetTDState_T(su2double t, const su2double *val_scalars = nullptr) override {
     /* The EoS only depends upon temperature. */
     Temperature = t;
     Density = Pressure / (Temperature * Gas_Constant);

diff --git a/SU2_CFD/include/variables/CIncNSVariable.hpp b/SU2_CFD/include/variables/CIncNSVariable.hpp
@@ -104,7 +104,7 @@ class CIncNSVariable final : public CIncEulerVariable {
   /*!
    * \brief Set all the primitive variables for incompressible flows
    */
-  bool SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2double turb_ke, CFluidModel *FluidModel) override;
+  bool SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2double turb_ke, CFluidModel *FluidModel, const su2double *scalar = nullptr);
   using CVariable::SetPrimVar;
 
   /*!

diff --git a/SU2_CFD/obj/Makefile.am b/SU2_CFD/obj/Makefile.am
@@ -51,6 +51,7 @@ endif
 libSU2Core_sources = ../src/definition_structure.cpp \
   ../src/fluid/CFluidModel.cpp \
   ../src/fluid/CIdealGas.cpp \
+  ../src/fluid/CFluidScalar.cpp \
   ../src/fluid/CPengRobinson.cpp \
   ../src/fluid/CVanDerWaalsGas.cpp \
   ../src/fluid/CNEMOGas.cpp \

diff --git a/SU2_CFD/src/fluid/CFluidModel.cpp b/SU2_CFD/src/fluid/CFluidModel.cpp
@@ -25,9 +25,10 @@
  * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include "../../include/fluid/CFluidModel.hpp"
+
 #include <utility>
 
-#include "../../include/fluid/CFluidModel.hpp"
 #include "../../include/fluid/CConstantConductivity.hpp"
 #include "../../include/fluid/CConstantConductivityRANS.hpp"
 #include "../../include/fluid/CConstantPrandtl.hpp"
@@ -61,7 +62,8 @@ void CFluidModel::SetThermalConductivityModel(const CConfig* config) {
         ThermalConductivity = unique_ptr<CConstantConductivityRANS>(
             new CConstantConductivityRANS(config->GetThermal_Conductivity_ConstantND(), config->GetPrandtl_Turb()));
       } else {
-        ThermalConductivity = unique_ptr<CConstantConductivity>(new CConstantConductivity(config->GetThermal_Conductivity_ConstantND()));
+        ThermalConductivity =
+            unique_ptr<CConstantConductivity>(new CConstantConductivity(config->GetThermal_Conductivity_ConstantND()));
       }
       break;
     case CONDUCTIVITYMODEL::CONSTANT_PRANDTL:

diff --git a/SU2_CFD/src/fluid/CFluidScalar.cpp b/SU2_CFD/src/fluid/CFluidScalar.cpp
@@ -0,0 +1,60 @@
+/*!
+ * \file CFluidScalar.hpp
+ * \brief Defines the multicomponent incompressible Ideal Gas model for mixtures.
+ * \author T. Economon, Mark Heimgartner, Cristopher Morales Ubal
+ * \version 7.3.1 "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/CFluidScalar.hpp"
+
+#include <math.h>
+
+#include <cmath>
+#include <iostream>
+#include <numeric>
+#include <vector>
+
+CFluidScalar::CFluidScalar(su2double val_Cp, su2double val_gas_constant, const su2double value_pressure_operating,
+                           CConfig* config)
+    : CFluidModel() {
+  n_species_mixture = config->GetnSpecies() + 1;
+
+  molarMasses.resize(n_species_mixture);
+
+  for (int iVar = 0; iVar < n_species_mixture; iVar++) {
+    molarMasses[iVar] = config->GetMolecular_Weight(iVar);
+  }
+
+  Pressure_Thermodynamic = value_pressure_operating;
+  Gas_Constant = val_gas_constant;
+  Gamma = 1.0;
+  Cp = val_Cp;
+  Cv = Cp;
+  SetLaminarViscosityModel(config);
+  SetThermalConductivityModel(config);
+}
+
+void CFluidScalar::SetTDState_T(su2double val_temperature, const su2double* val_scalars) {
+  const su2double MeanMolecularWeight = ComputeMeanMolecularWeight(molarMasses, val_scalars);
+  Density = Pressure_Thermodynamic / ((val_temperature * UNIVERSAL_GAS_CONSTANT / MeanMolecularWeight));
+}
diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build
@@ -5,6 +5,7 @@ su2_cfd_src = files(['definition_structure.cpp',
 
 su2_cfd_src += files(['fluid/CFluidModel.cpp',
                       'fluid/CIdealGas.cpp',
+                      'fluid/CFluidScalar.cpp',
                       'fluid/CPengRobinson.cpp',
                       'fluid/CVanDerWaalsGas.cpp',
                       'fluid/CNEMOGas.cpp',