Merge branch 'develop' into fix_SpecialOutput_SonicBoom

pcarruscag · web-flow · commit 5cb3a3c327be · 2021-08-09T14:58:51.000+01:00
diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
@@ -475,7 +475,6 @@ class CConfig {
   unsigned short Kind_Solver,      /*!< \brief Kind of solver Euler, NS, Continuous adjoint, etc.  */
   Kind_FluidModel,                 /*!< \brief Kind of the Fluid Model: Ideal or Van der Walls, ... . */
   Kind_InitOption,                 /*!< \brief Kind of Init option to choose if initializing with Reynolds number or with thermodynamic conditions   */
-  Kind_TransCoeffModel,            /*!< \brief Transport coefficient Model for NEMO solver. */
   Kind_GridMovement,               /*!< \brief Kind of the static mesh movement. */
   *Kind_SurfaceMovement,           /*!< \brief Kind of the static mesh movement. */
   nKind_SurfaceMovement,           /*!< \brief Kind of the dynamic mesh movement. */
@@ -1142,6 +1141,7 @@ class CConfig {
   monoatomic;                               /*!< \brief Flag for monoatomic mixture. */
   string GasModel,                          /*!< \brief Gas Model. */
   *Wall_Catalytic;                          /*!< \brief Pointer to catalytic walls. */
+  TRANSCOEFFMODEL   Kind_TransCoeffModel;   /*!< \brief Transport coefficient Model for NEMO solver. */
 
   /*!
    * \brief Set the default values of config options not set in the config file using another config object.
@@ -3623,7 +3623,7 @@ class CConfig {
    * \brief Get the transport coefficient model.
    * \return Index of transport coefficient model.
    */
-  unsigned short GetKind_TransCoeffModel(void) const { return Kind_TransCoeffModel; }
+  TRANSCOEFFMODEL GetKind_TransCoeffModel(void) const { return Kind_TransCoeffModel; }
 
   /*!
    * \brief Get the total number of heat flux markers.
diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
@@ -540,13 +540,15 @@ MakePair("ONESPECIES", ONESPECIES)
 /*!
  * \brief types of coefficient transport model
  */
-enum ENUM_TRANSCOEFFMODEL {
-  WILKE      = 0,
-  GUPTAYOS   = 1
-};
-static const MapType<std::string, ENUM_TRANSCOEFFMODEL> TransCoeffModel_Map = {
-MakePair("WILKE", WILKE)
-MakePair("GUPTA-YOS", GUPTAYOS)
+enum class TRANSCOEFFMODEL {
+  WILKE,
+  GUPTAYOS,
+  CHAPMANN_ENSKOG
+};
+static const MapType<std::string, TRANSCOEFFMODEL> TransCoeffModel_Map = {
+MakePair("WILKE", TRANSCOEFFMODEL::WILKE)
+MakePair("GUPTA-YOS", TRANSCOEFFMODEL::GUPTAYOS)
+MakePair("CHAPMANN-ENSKOG", TRANSCOEFFMODEL::CHAPMANN_ENSKOG)
 };
 
 /*!
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
@@ -1149,7 +1149,7 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Specify chemical model for multi-species simulations - read by Mutation++ library*/
   addStringOption("GAS_MODEL", GasModel, string("N2"));
   /* DESCRIPTION: Specify transport coefficient model for multi-species simulations */
-  addEnumOption("TRANSPORT_COEFF_MODEL", Kind_TransCoeffModel, TransCoeffModel_Map, WILKE);
+  addEnumOption("TRANSPORT_COEFF_MODEL", Kind_TransCoeffModel, TransCoeffModel_Map, TRANSCOEFFMODEL::WILKE);
   /* DESCRIPTION: Specify mass fraction of each species */
   addDoubleListOption("GAS_COMPOSITION", nSpecies, Gas_Composition);
   /* DESCRIPTION: Specify if mixture is frozen */
@@ -3587,8 +3587,12 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     SU2_MPI::Error("Only STANDARD_AIR fluid model can be used with US Measurement System", CURRENT_FUNCTION);
   }
 
-  if (nemo && Kind_TransCoeffModel != WILKE ) {
-    SU2_MPI::Error("Only WILKE transport model is stable for the NEMO solver.", CURRENT_FUNCTION);
+  if (Kind_FluidModel == SU2_NONEQ && Kind_TransCoeffModel != TRANSCOEFFMODEL::WILKE ) {
+    SU2_MPI::Error("Only WILKE transport model is stable for the NEMO solver using SU2TClib. Use Mutation++ instead.", CURRENT_FUNCTION);
+  }
+
+  if (Kind_FluidModel == MUTATIONPP && (Kind_TransCoeffModel != TRANSCOEFFMODEL::WILKE && Kind_TransCoeffModel != TRANSCOEFFMODEL::CHAPMANN_ENSKOG)) {
+    SU2_MPI::Error("Only WILKE and Chapmann-Enskog transport model can be used with Mutation++ at the moment.", CURRENT_FUNCTION);
   }
 
   if (!ideal_gas && !nemo) {
diff --git a/SU2_CFD/include/fluid/CNEMOGas.hpp b/SU2_CFD/include/fluid/CNEMOGas.hpp
@@ -48,8 +48,8 @@ class CNEMOGas : public CFluidModel {
   nHeavy,                                /*!< \brief Number of heavy particles in gas */
   nEl,                                   /*!< \brief Number of electrons in gas */
   nDim,                                  /*!< \brief Number of dimensions. */
-  nEnergyEq = 2,                         /*!< \brief Number of energy equations for the 2T model. */
-  Kind_TransCoeffModel;                  /*!< \brief Transport coefficients model for NEMO solver. */
+  nEnergyEq = 2;                         /*!< \brief Number of energy equations for the 2T model. */
+  TRANSCOEFFMODEL Kind_TransCoeffModel;  /*!< \brief Transport coefficients model for NEMO solver. */
 
   unsigned iSpecies,                     /*!< \brief Common iteration counter for species */
   jSpecies,                              /*!< \brief Common iteration counter for species */
diff --git a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp
@@ -81,6 +81,7 @@ class CNEMOEulerVariable : public CVariable {
   LAM_VISC_INDEX, EDDY_VISC_INDEX, nSpecies;
 
   su2double Tve_Freestream; /*!< \brief Freestream vib-el temperature. */
+  const bool implicit;      /*!< \brief Implicit flag. */
 
 public:
 
diff --git a/SU2_CFD/src/fluid/CMutationTCLib.cpp b/SU2_CFD/src/fluid/CMutationTCLib.cpp
@@ -42,10 +42,12 @@ CMutationTCLib::CMutationTCLib(const CConfig* config, unsigned short val_nDim):
   /*--- Set up inputs to define type of mixture in the Mutation++ library ---*/
 
   /*--- Define transport model ---*/
-  if(Kind_TransCoeffModel == WILKE)
+  if(Kind_TransCoeffModel == TRANSCOEFFMODEL::WILKE)
     transport_model = "Wilke";
-  else if (Kind_TransCoeffModel == GUPTAYOS)
+  else if (Kind_TransCoeffModel == TRANSCOEFFMODEL::GUPTAYOS)
     transport_model = "Gupta-Yos";
+  else if (Kind_TransCoeffModel == TRANSCOEFFMODEL::CHAPMANN_ENSKOG)
+    transport_model = "Chapmann-Enskog_LDLT";
 
   opt.setStateModel("ChemNonEqTTv");
   if (frozen) opt.setMechanism("none");
diff --git a/SU2_CFD/src/fluid/CSU2TCLib.cpp b/SU2_CFD/src/fluid/CSU2TCLib.cpp
@@ -1004,9 +1004,9 @@ vector<su2double>& CSU2TCLib::ComputeSpeciesEnthalpy(su2double val_T, su2double
 
 vector<su2double>& CSU2TCLib::GetDiffusionCoeff(){
 
-  if(Kind_TransCoeffModel == WILKE)
+  if(Kind_TransCoeffModel == TRANSCOEFFMODEL::WILKE)
    DiffusionCoeffWBE();
-  if(Kind_TransCoeffModel == GUPTAYOS)
+  if(Kind_TransCoeffModel == TRANSCOEFFMODEL::GUPTAYOS)
    DiffusionCoeffGY();
 
   return DiffusionCoeff;
@@ -1015,9 +1015,9 @@ vector<su2double>& CSU2TCLib::GetDiffusionCoeff(){
 
 su2double CSU2TCLib::GetViscosity(){
 
-  if(Kind_TransCoeffModel == WILKE)
+  if(Kind_TransCoeffModel == TRANSCOEFFMODEL::WILKE)
     ViscosityWBE();
-  if(Kind_TransCoeffModel == GUPTAYOS)
+  if(Kind_TransCoeffModel == TRANSCOEFFMODEL::GUPTAYOS)
     ViscosityGY();
 
   return Mu;
@@ -1026,9 +1026,9 @@ su2double CSU2TCLib::GetViscosity(){
 
 vector<su2double>& CSU2TCLib::GetThermalConductivities(){
 
-  if(Kind_TransCoeffModel == WILKE)
+  if(Kind_TransCoeffModel == TRANSCOEFFMODEL::WILKE)
     ThermalConductivitiesWBE();
-  if(Kind_TransCoeffModel == GUPTAYOS)
+  if(Kind_TransCoeffModel == TRANSCOEFFMODEL::GUPTAYOS)
     ThermalConductivitiesGY();
 
   return ThermalConductivities;
diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp
@@ -1286,16 +1286,21 @@ void CNEMOEulerSolver::SetNondimensionalization(CConfig *config, unsigned short
     if (viscous) {
 
       switch(config->GetKind_TransCoeffModel()){
-      case WILKE:
+      case TRANSCOEFFMODEL::WILKE:
         ModelTable << "Wilke-Blottner-Eucken";
         NonDimTable.PrintFooter();
         break;
 
-      case GUPTAYOS:
+      case TRANSCOEFFMODEL::GUPTAYOS:
         ModelTable << "Gupta-Yos";
         NonDimTable.PrintFooter();
         break;
 
+      case TRANSCOEFFMODEL::CHAPMANN_ENSKOG:
+        ModelTable << "CHAPMANN-ENSKOG_LDLT";
+        NonDimTable.PrintFooter();
+        break;
+
       default:
         break;
       }
diff --git a/SU2_CFD/src/variables/CNEMOEulerVariable.cpp b/SU2_CFD/src/variables/CNEMOEulerVariable.cpp
@@ -43,7 +43,8 @@ CNEMOEulerVariable::CNEMOEulerVariable(su2double val_pressure,
                                                                          ndim,
                                                                          nvar,
                                                                          config ),
-                                       Gradient_Reconstruction(config->GetReconstructionGradientRequired() ? Gradient_Aux : Gradient_Primitive) {
+                                       Gradient_Reconstruction(config->GetReconstructionGradientRequired() ? Gradient_Aux : Gradient_Primitive),
+                                       implicit(config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT) {
 
   unsigned short iDim, iSpecies;
 
@@ -267,41 +268,20 @@ bool CNEMOEulerVariable::Cons2PrimVar(su2double *U, su2double *V,
   V[TVE_INDEX] = T[1];
 
   // Determine if the temperature lies within the acceptable range
-  if (V[T_INDEX] <= Tmin)      { nonPhys = true;}
-  else if (V[T_INDEX] >= Tmax) { nonPhys = true;}
-  else if (V[T_INDEX] != V[T_INDEX] || V[TVE_INDEX] != V[TVE_INDEX]){
-    nonPhys = true;}
+  if (V[T_INDEX] <= Tmin)      { nonPhys = true; return nonPhys;}
+  else if (V[T_INDEX] >= Tmax) { nonPhys = true; return nonPhys;}
+  else if (V[T_INDEX] != V[T_INDEX]){ nonPhys = true; return nonPhys;}
 
-  /*--- Vibrational-Electronic Temperature ---*/
-  vector<su2double> eves_min = fluidmodel->ComputeSpeciesEve(Tvemin);
-  vector<su2double> eves_max = fluidmodel->ComputeSpeciesEve(Tvemax);
-
-  // Check for non-physical solutions
   if (!monoatomic){
-    su2double rhoEve_min = 0.0;
-    su2double rhoEve_max = 0.0;
-    for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
-      rhoEve_min += U[iSpecies] * eves_min[iSpecies];
-      rhoEve_max += U[iSpecies] * eves_max[iSpecies];
-    }
-
-    if (rhoEve < rhoEve_min) {
-
-      nonPhys      = true;
-      V[TVE_INDEX] = Tvemin;
-      U[nSpecies+nDim+1] = rhoEve_min;
-    } else if (rhoEve > rhoEve_max) {
-      nonPhys      = true;
-      V[TVE_INDEX] = Tvemax;
-      U[nSpecies+nDim+1] = rhoEve_max;
-    }
-  } else {
-    //TODO: can e-modes/vibe modes be active?
-    V[TVE_INDEX] = Tve_Freestream;
+    if (V[TVE_INDEX] <= Tvemin)      { nonPhys = true; return nonPhys;}
+    else if (V[TVE_INDEX] >= Tvemax) { nonPhys = true; return nonPhys;}
+    else if (V[TVE_INDEX] != V[TVE_INDEX]){ nonPhys = true; return nonPhys;}
   }
+  else {V[TVE_INDEX] = Tve_Freestream;}
 
   // Determine other properties of the mixture at the current state
   fluidmodel->SetTDStateRhosTTv(rhos, V[T_INDEX], V[TVE_INDEX]);
+
   const auto& cvves = fluidmodel->ComputeSpeciesCvVibEle();
   vector<su2double> eves  = fluidmodel->ComputeSpeciesEve(V[TVE_INDEX]);
 
@@ -310,11 +290,8 @@ bool CNEMOEulerVariable::Cons2PrimVar(su2double *U, su2double *V,
     val_Cvves[iSpecies] = cvves[iSpecies];
   }
 
-  su2double rhoCvtr = fluidmodel->ComputerhoCvtr();
-  su2double rhoCvve = fluidmodel->ComputerhoCvve();
-
-  V[RHOCVTR_INDEX] = rhoCvtr;
-  V[RHOCVVE_INDEX] = rhoCvve;
+  V[RHOCVTR_INDEX] = fluidmodel->ComputerhoCvtr();
+  V[RHOCVVE_INDEX] = fluidmodel->ComputerhoCvve();
 
   /*--- Pressure ---*/
   V[P_INDEX] = fluidmodel->ComputePressure();
@@ -325,9 +302,11 @@ bool CNEMOEulerVariable::Cons2PrimVar(su2double *U, su2double *V,
   }
 
   /*--- Partial derivatives of pressure and temperature ---*/
-  fluidmodel->ComputedPdU  (V, eves, val_dPdU  );
-  fluidmodel->ComputedTdU  (V, val_dTdU );
-  fluidmodel->ComputedTvedU(V, eves, val_dTvedU);
+  if(implicit){
+    fluidmodel->ComputedPdU  (V, eves, val_dPdU  );
+    fluidmodel->ComputedTdU  (V, val_dTdU );
+    fluidmodel->ComputedTvedU(V, eves, val_dTvedU);
+  }
 
   /*--- Sound speed ---*/
   V[A_INDEX] = fluidmodel->ComputeSoundSpeed();
diff --git a/config_template.cfg b/config_template.cfg
@@ -1260,6 +1260,12 @@ TIME_DISCRE_ADJTURB= EULER_IMPLICIT
 % Reduction factor of the CFL coefficient in the adjoint turbulent problem
 CFL_REDUCTION_ADJTURB= 0.01
 
+
+% -------------------- NEMO NUMERICAL METHOD DEFINITION -----------------------%
+%
+% Mixture transport properties (WILKE,GUPTA-YOS,CHAPMANN-ENSKOG)
+TRANSPORT_COEFF_MODEL = WILKE
+
 % ----------------------- GEOMETRY EVALUATION PARAMETERS ----------------------%
 %
 % Marker(s) of the surface where geometrical based function will be evaluated
