[WIP] Composition dependent fluid properties for multicomponent flow

.github/workflows/regression.yml

@@ -9,6 +9,7 @@ on:
     branches: 
       - 'develop'
       - 'master'
+      - 'feature_flamelet'
 
 jobs:
   build:
@@ -60,14 +61,16 @@ jobs:
     strategy:
       fail-fast: false
       matrix: 
-        testscript: ['tutorials.py', 'parallel_regression.py', 'flamelet_regression.py', 'parallel_regression_AD.py', 'serial_regression.py', 'serial_regression_AD.py', 'hybrid_regression.py']
+        testscript: ['tutorials.py', 'parallel_regression.py', 'flamelet_regression.py', 'multicomp_regression.py', 'parallel_regression_AD.py', 'serial_regression.py', 'serial_regression_AD.py', 'hybrid_regression.py']
         include:
           - testscript: 'tutorials.py'
             tag: MPI
           - testscript: 'parallel_regression.py'
             tag: MPI
           - testscript: 'flamelet_regression.py'
             tag: MPI
+          - testscript: 'multicomp_regression.py'
+            tag: MPI
           - testscript: 'parallel_regression_AD.py'
             tag: MPI
           - testscript: 'serial_regression.py'
@@ -92,7 +95,7 @@ jobs:
         uses: docker://su2code/test-su2:20200303
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c feature_flamelet -s ${{matrix.testscript}}
+          args: -b ${{github.ref}} -t develop -c feature_multicomp -s ${{matrix.testscript}}
   unit_tests:
     runs-on: ubuntu-latest
     name: Unit Tests

Common/include/option_structure.hpp

@@ -564,6 +564,7 @@ enum ENUM_FLUIDMODEL {
   MUTATIONPP           = 7, /*!< \brief Mutation++ gas model for nonequilibrium flow. */
   SU2_NONEQ            = 8, /*!< \brief User defined gas model for nonequilibrium flow. */
   FLAMELET_FLUID_MODEL = 9, /*!, \brief Flamelet model */
+  MIXTURE_FLUID_MODEL  = 10,/*!, \brief Mixture model */
 };
 static const MapType<std::string, ENUM_FLUIDMODEL> FluidModel_Map = {
   MakePair("STANDARD_AIR", STANDARD_AIR)
@@ -575,6 +576,7 @@ static const MapType<std::string, ENUM_FLUIDMODEL> FluidModel_Map = {
   MakePair("INC_IDEAL_GAS_POLY", INC_IDEAL_GAS_POLY)
   MakePair("MUTATIONPP", MUTATIONPP)
   MakePair("FLAMELET_FLUID_MODEL", FLAMELET_FLUID_MODEL)
+  MakePair("MIXTURE_FLUID_MODEL", MIXTURE_FLUID_MODEL)
   MakePair("SU2_NONEQ", SU2_NONEQ)
 };
 
@@ -712,12 +714,14 @@ enum class DIFFUSIVITYMODEL {
   CONSTANT_DIFFUSIVITY, /*!< \brief Constant mass diffusivity for scalar transport. */
   CONSTANT_SCHMIDT,     /*!< \brief Constant Schmidt number for mass diffusion in scalar transport. */
   FLAMELET,             /*!< \brief flamelet model */
+  UNITY_LEWIS,             /*!< \brief Unity Lewis model */
 };
 
 static const MapType<std::string, DIFFUSIVITYMODEL> DiffusivityModel_Map = {
   MakePair("CONSTANT_DIFFUSIVITY", DIFFUSIVITYMODEL::CONSTANT_DIFFUSIVITY)
   MakePair("CONSTANT_SCHMIDT", DIFFUSIVITYMODEL::CONSTANT_SCHMIDT)
   MakePair("FLAMELET", DIFFUSIVITYMODEL::FLAMELET)
+  MakePair("UNITY_LEWIS", DIFFUSIVITYMODEL::UNITY_LEWIS)
 };
 
 /*!

SU2_CFD/include/fluid/CConstantDiffusivity.hpp

@@ -50,7 +50,8 @@ class CConstantDiffusivity final : public CDiffusivityModel {
                       su2double rho,
                       su2double mu_lam,
                       su2double mu_turb,
-                      su2double cp) override { }
+                      su2double cp,
+                      su2double kt) override { }
 
   private:
     su2double diff_{0.0};  

SU2_CFD/include/fluid/CConstantSchmidt.hpp

@@ -62,7 +62,7 @@ class CConstantSchmidt final: public CDiffusivityModel {
   /*!
    * \brief Set diffusivity.
    */
-  void SetDiffusivity(su2double T, su2double rho, su2double mu_lam, su2double mu_turb, su2double cp) override {
+  void SetDiffusivity(su2double T, su2double rho, su2double mu_lam, su2double mu_turb, su2double cp, su2double kt) override {
     diff_ = mu_lam / (rho*sc_lam_);
   }
 

SU2_CFD/include/fluid/CConstantSchmidtRANS.hpp

@@ -57,7 +57,7 @@ class CConstantSchmidtRANS : public CDiffusivityModel {
   /*!
    * \brief Set mass diffusivity.
    */
-  void SetDiffusivity(su2double t, su2double rho, su2double mu_lam, su2double mu_turb, su2double cp) override {
+  void SetDiffusivity(su2double t, su2double rho, su2double mu_lam, su2double mu_turb, su2double cp, su2double kt) override {
     diff_ = diff_lam_const_ + mu_turb/(rho*sc_turb_);
   }
 

SU2_CFD/include/fluid/CDiffusivityModel.hpp

@@ -47,5 +47,5 @@ class CDiffusivityModel {
 /*!
  * \brief Set mass diffusivity
  */
-  virtual void SetDiffusivity(su2double T, su2double rho, su2double mu_lam, su2double mu_turb, su2double cp) = 0;
+  virtual void SetDiffusivity(su2double T, su2double rho, su2double mu_lam, su2double mu_turb, su2double cp, su2double kt) = 0;
 };

SU2_CFD/include/fluid/CFluidModel.hpp

@@ -125,6 +125,36 @@ class CFluidModel {
    */
   su2double GetCv() const { return Cv; }
 
+  /*!
+   * \brief Get fluid mean specific heat capacity at constant pressure.
+   */
+  su2double ComputeMeanSpecificHeatCp(const std::vector<su2double> specific_heat_cp, const su2double * const val_scalars) {
+    su2double val_scalars_sum = 0.0;
+    unsigned short n_scalars = specific_heat_cp.size() - 1;
+
+    for (int i_scalar = 0; i_scalar < n_scalars; i_scalar++){
+      Cp += specific_heat_cp[i_scalar] * val_scalars[i_scalar];
+      val_scalars_sum += val_scalars[i_scalar];
+    }
+    return Cp += specific_heat_cp[n_scalars]*(1 - val_scalars_sum);
+  }
+
+  /*!
+   * \brief Get fluid mean molecular weight.
+   */
+  static su2double ComputeMeanMolecularWeight(const std::vector<su2double> molar_masses, const su2double * const val_scalars) {
+    su2double OneOverMeanMolecularWeight = 0.0;
+    su2double val_scalars_sum = 0.0;
+    unsigned short n_scalars = molar_masses.size() - 1;
+
+    for (int i_scalar = 0; i_scalar < n_scalars; 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[n_scalars]/1000);
+    return 1/OneOverMeanMolecularWeight;
+  }
+
   /*!
    * \brief Get the source term of the transported scalar
    */
@@ -164,11 +194,11 @@ class CFluidModel {
    * \brief Get fluid dynamic viscosity.
    */
   virtual inline su2double GetLaminarViscosity() {
-    LaminarViscosity->SetViscosity(Temperature, Density);
-    Mu = LaminarViscosity->GetViscosity();
-    LaminarViscosity->SetDerViscosity(Temperature, Density);
-    dmudrho_T = LaminarViscosity->Getdmudrho_T();
-    dmudT_rho = LaminarViscosity->GetdmudT_rho();
+      LaminarViscosity->SetViscosity(Temperature, Density);
+      Mu = LaminarViscosity->GetViscosity();
+      LaminarViscosity->SetDerViscosity(Temperature, Density);
+      dmudrho_T = LaminarViscosity->Getdmudrho_T();
+      dmudT_rho = LaminarViscosity->GetdmudT_rho();
     return Mu;
   }
 
@@ -188,7 +218,7 @@ class CFluidModel {
    * \brief Get fluid mass diffusivity.
    */
   virtual inline su2double GetMassDiffusivity() {
-    MassDiffusivity->SetDiffusivity(Temperature, Density, Mu, Mu_Turb, Cp);
+    MassDiffusivity->SetDiffusivity(Temperature, Density, Mu, Mu_Turb, Cp, Kt);
     mass_diffusivity = MassDiffusivity->GetDiffusivity();
     //ThermalConductivity->SetDerConductivity(Temperature, Density, dmudrho_T, dmudT_rho, Cp);
     //dktdrho_T = ThermalConductivity->Getdktdrho_T();
@@ -264,12 +294,12 @@ class CFluidModel {
   /*!
    * \brief Set viscosity model.
    */
-  void SetLaminarViscosityModel(const CConfig* config);
+  virtual void SetLaminarViscosityModel(const CConfig* config);
 
   /*!
    * \brief Set thermal conductivity model.
    */
-  void SetThermalConductivityModel(const CConfig* config);
+  virtual void SetThermalConductivityModel(const CConfig* config);
 
   /*!
    * \brief Set mass diffusivity model.

SU2_CFD/include/fluid/CFluidScalar.hpp

@@ -0,0 +1,86 @@
+#pragma once
+
+#include <vector>
+#include <memory>
+
+#include "CFluidModel.hpp"
+
+class CFluidScalar final : public CFluidModel {
+private:
+  unsigned short n_scalars = 0;                       /*!< \brief Number of transported scalars. */
+  unsigned short n_species_mixture = 0;               /*!< \brief Number of species in mixture. */
+  su2double Gas_Constant = 0.0;                       /*!< \brief Specific gas constant. */
+  su2double Gamma = 0.0;                              /*!< \brief Ratio of specific heats of the gas. */
+  su2double Pressure_Thermodynamic = 0.0;             /*!< \brief Constant pressure thermodynamic. */
+
+  bool wilke;
+  bool davidson;
+
+  std::vector<su2double> massFractions;               /*!< \brief Mass fractions of all species. */
+  std::vector<su2double> moleFractions;               /*!< \brief Mole fractions of all species. */
+  std::vector<su2double> molarMasses;                 /*!< \brief Molar masses of all species. */
+  std::vector<su2double> laminarViscosity;            /*!< \brief Laminar viscosity of all species. */
+  std::vector<su2double> specificHeat;                /*!< \brief Specific heat of all species. */
+  std::vector<su2double> laminarthermalConductivity;  /*!< \brief Laminar thermal conductivity of all species. */
+
+  static const int ARRAYSIZE = 100;
+  std::unique_ptr<CViscosityModel> LaminarViscosityPointers[ARRAYSIZE];
+  std::unique_ptr<CConductivityModel> ThermalConductivityPointers[ARRAYSIZE];
+
+  /*!
+   * \brief Convert mass fractions to mole fractions.
+   * \param[in] val_scalars - Scalar mass fraction.
+   */
+  std::vector<su2double>& massToMoleFractions(const su2double * const val_scalars);
+
+  /*!
+   * \brief Wilke mixing law for mixture viscosity.
+   * \param[in] val_scalars - Scalar mass fraction.
+   */
+  su2double wilkeViscosity(const su2double * const val_scalars);
+
+  /*!
+   * \brief Davidson mixing law for mixture viscosity.
+   * \param[in] val_scalars - Scalar mass fraction.
+   */
+  su2double davidsonViscosity(const su2double * const val_scalars);
+
+  /*!
+   * \brief Wilke mixing law for mixture thermal conductivity.
+   * \param[in] val_scalars - Scalar mass fraction.
+   */
+  su2double wilkeConductivity(const su2double * const val_scalars);
+
+public:
+  /*!
+   * \brief Constructor of the class.
+   */
+  CFluidScalar(CConfig *config, const su2double value_pressure_operating);
+
+  /*!
+   * \brief Set viscosity model.
+   */
+  void SetLaminarViscosityModel(const CConfig* config) override;
+
+  /*!
+   * \brief Set thermal conductivity model.
+   */
+  void SetThermalConductivityModel(const CConfig* config) override;
+
+  /*!
+   * \brief Set the Dimensionless State using Temperature.
+   * \param[in] val_temperature - Temperature value at the point.
+   * \param[in] val_scalars - Scalar mass fraction.
+   */
+  unsigned long SetTDState_T(const su2double val_temperature, su2double * const val_scalars) override;
+
+  /*!
+   * \brief Get fluid dynamic viscosity.
+   */
+  inline su2double GetLaminarViscosity() override {return Mu; }
+
+  /*!
+   * \brief Get fluid thermal conductivity.
+   */
+  inline su2double GetThermalConductivity() override { return Kt; }
+};

SU2_CFD/include/fluid/CUnityLewisDiffusivity.hpp

@@ -0,0 +1,57 @@
+/*!
+ * \file CUnityLewisDiffusivity.hpp
+ * \brief Defines unity Lewis mass diffusivity.
+ * \author M.Heimgartner
+ * \version 7.0.6 "Blackbird"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2020, 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 "CDiffusivityModel.hpp"
+
+/*!
+ * \class CUnityLewisDiffusivity
+ * \brief Defines a unity Lewis mass diffusivity model for species equations.
+ * \author M.Heimgartner 
+ */
+class CUnityLewisDiffusivity final : public CDiffusivityModel {
+public:
+  /*!
+   * \brief Constructor of the class.
+   */
+  CUnityLewisDiffusivity() {}
+
+  su2double GetDiffusivity() const override {return diff_;}
+
+  /*!
+   * \brief Set diffusivity.
+   */
+  void SetDiffusivity(su2double T, su2double rho, su2double mu_lam, su2double mu_turb, su2double cp, su2double kt) override { 
+    diff_ = kt / (Lewis * rho * cp);
+  }
+
+  private:
+    su2double diff_{0.0}; 
+    su2double kt_{0.0}; 
+    su2double Lewis{1};  
+};

SU2_CFD/include/numerics/flamelet/scalarLegacy_sources.hpp

@@ -104,7 +104,7 @@ class CSourcePieceWise_transportedScalar_general final : public CNumerics {
           if (flame)
             Residual[iVar] += yinv*Volume*(Diffusion_Coeff_i[iVar]+Mass_Diffusivity_Tur)*ScalarVar_Grad_i[iVar][1];
           else
-            Residual[iVar] += yinv*Volume*Density_i*(Diffusion_Coeff_i[iVar]+Mass_Diffusivity_Tur)*ScalarVar_Grad_i[iVar][1];
+            Residual[iVar] += yinv*Volume*(Density_i*Diffusion_Coeff_i[iVar]+Mass_Diffusivity_Tur)*ScalarVar_Grad_i[iVar][1];
         } 
       }
 

SU2_CFD/include/solvers/CPassiveScalarSolver.hpp

@@ -88,6 +88,15 @@ class CPassiveScalarSolver: public CScalarLegacySolver {
                       CConfig *config,
                       unsigned short iMesh);
 
+   /*!
+   * \brief Compute the primitive variables (diffusivities)
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] Output - boolean to determine whether to print output.
+   * \return - The number of non-physical points.
+   */
+  unsigned long SetPrimitive_Variables(CSolver **solver_container, CConfig *config, bool Output);
+
 
   /*!
    * \brief Set the initial condition for the scalar transport problem.

SU2_CFD/src/fluid/CFluidFlamelet.cpp

@@ -98,7 +98,7 @@ unsigned long CFluidFlamelet::SetScalarSources(su2double *val_scalars){
 
   delete [] table_sources;
 
-  /*--- not used at the moment ---*/ 
+  /*--- not used at the moment ---*/
   return exit_code;
 }
 

SU2_CFD/src/fluid/CFluidModel.cpp

@@ -37,9 +37,11 @@
 #include "../../include/fluid/CPolynomialConductivityRANS.hpp"
 #include "../../include/fluid/CPolynomialViscosity.hpp"
 #include "../../include/fluid/CSutherland.hpp"
+#include "../../include/fluid/CFluidScalar.hpp"
 #include "../../include/fluid/CConstantDiffusivity.hpp"
 #include "../../include/fluid/CConstantSchmidtRANS.hpp"
 #include "../../include/fluid/CConstantSchmidt.hpp"
+#include "../../include/fluid/CUnityLewisDiffusivity.hpp"
 
 void CFluidModel::SetLaminarViscosityModel(const CConfig* config) {
   switch (config->GetKind_ViscosityModel()) {
@@ -114,6 +116,9 @@ void CFluidModel::SetMassDiffusivityModel (const CConfig* config) {
     case DIFFUSIVITYMODEL::FLAMELET:
       /* do nothing. Diffusivity is obtained from the table and set in setTDState_T */
       break;
+    case DIFFUSIVITYMODEL::UNITY_LEWIS:
+      MassDiffusivity = unique_ptr<CUnityLewisDiffusivity>(new CUnityLewisDiffusivity());
+      break;
     default:
       SU2_MPI::Error("Diffusivity model not available.", CURRENT_FUNCTION);
       break;