Merge pull request #1887 from su2code/feature_NEMO_visc_ion

Adding ionization capabilities for viscous simulations (NEMO)

Author: WallyMaier

Common/include/option_structure.hpp

@@ -94,6 +94,7 @@ const su2double STANDARD_GRAVITY = 9.80665;           /*!< \brief Acceleration d
 const su2double UNIVERSAL_GAS_CONSTANT = 8.3144598;   /*!< \brief Universal gas constant in J/(mol*K) */
 const su2double BOLTZMANN_CONSTANT = 1.3806503E-23;   /*!< \brief Boltzmann's constant [J K^-1] */
 const su2double AVOGAD_CONSTANT = 6.0221415E26; /*!< \brief Avogadro's constant, number of particles in one kmole. */
+const su2double FUND_ELEC_CHARGE_CGS = 4.8032047E-10; /*!< \brief Fundamental electric charge in CGS units, cm^(3/2) g^(1/2) s^(-1). */
 
 const su2double EPS = 1.0E-16;        /*!< \brief Error scale. */
 const su2double TURB_EPS = 1.0E-16;   /*!< \brief Turbulent Error scale. */

Common/src/CConfig.cpp

@@ -3931,16 +3931,27 @@ 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 (Kind_FluidModel == SU2_NONEQ && (Kind_TransCoeffModel != TRANSCOEFFMODEL::WILKE && Kind_TransCoeffModel != TRANSCOEFFMODEL::SUTHERLAND) ) {
-      SU2_MPI::Error("Only WILKE and SUTHERLAND transport models are stable for the NEMO solver using SU2TClib. Use Mutation++ instead.", CURRENT_FUNCTION);
+    /* --- Check for NEMO compatibility issues ---*/
+    if (Kind_FluidModel == SU2_NONEQ && (Kind_TransCoeffModel != TRANSCOEFFMODEL::WILKE && Kind_TransCoeffModel != TRANSCOEFFMODEL::SUTHERLAND && Kind_TransCoeffModel != TRANSCOEFFMODEL::GUPTAYOS) ) {
+      SU2_MPI::Error("Transport model not available for NEMO solver using SU2TCLIB. Please use the WILKE, SUTHERLAND or GUPTAYOS transport model instead.", CURRENT_FUNCTION);
+    }
+
+    if (Kind_Solver == MAIN_SOLVER::NEMO_NAVIER_STOKES) {
+      if (Kind_FluidModel == SU2_NONEQ && GasModel == "AIR-7" && Kind_TransCoeffModel != TRANSCOEFFMODEL::GUPTAYOS) {
+        SU2_MPI::Error("Only Gupta-Yos transport model available for ionized flows using SU2TCLIB.", 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.",
+      SU2_MPI::Error("Transport model not available for NEMO solver using MUTATIONPP. Please use the WILKE or CHAPMANN_ENSKOG transport model instead..",
                      CURRENT_FUNCTION);
     }
 
+    if (Kind_FluidModel == SU2_NONEQ && GasModel == "AIR-7" && nWall_Catalytic != 0) {
+      SU2_MPI::Error("Catalytic wall recombination is not yet available for ionized flows in SU2_NEMO.", CURRENT_FUNCTION);
+    }
+
     if (!ideal_gas && !nemo) {
       if (Kind_Upwind_Flow != UPWIND::ROE && Kind_Upwind_Flow != UPWIND::HLLC && Kind_Centered_Flow != CENTERED::JST) {
         SU2_MPI::Error("Only ROE Upwind, HLLC Upwind scheme, and JST scheme can be used for Non-Ideal Compressible Fluids", CURRENT_FUNCTION);
@@ -6086,14 +6097,14 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
         if (Kind_Regime == ENUM_REGIME::COMPRESSIBLE) cout << "Compressible two-temperature thermochemical non-equilibrium Euler equations." << endl;
         if (Kind_FluidModel == SU2_NONEQ){
           if ((GasModel != "N2") && (GasModel != "AIR-5") && (GasModel != "AIR-7") && (GasModel != "ARGON"))
-          SU2_MPI::Error("The GAS_MODEL given as input is not valid. Choose one of the options: N2, AIR-5, AIR-7, ARGON.", CURRENT_FUNCTION);
+            SU2_MPI::Error("The GAS_MODEL given is unavailable using CSU2TCLIB. Choose one of the options: N2, AIR-5, AIR-7, or ARGON.", CURRENT_FUNCTION);
         }
         break;
       case MAIN_SOLVER::NEMO_NAVIER_STOKES:
         if (Kind_Regime == ENUM_REGIME::COMPRESSIBLE) cout << "Compressible two-temperature thermochemical non-equilibrium Navier-Stokes equations." << endl;
         if (Kind_FluidModel == SU2_NONEQ){
-          if ((GasModel != "N2") && (GasModel != "AIR-5") && (GasModel != "ARGON"))
-          SU2_MPI::Error("The GAS_MODEL given as input is not valid. Choose one of the options: N2, AIR-5, ARGON.", CURRENT_FUNCTION);
+          if ((GasModel != "N2") && (GasModel != "AIR-5") && (GasModel != "AIR-7") && (GasModel != "ARGON"))
+          SU2_MPI::Error("The GAS_MODEL given is unavailable using CSU2TCLIB. Choose one of the options: N2, AIR-5, AIR-7, or ARGON.", CURRENT_FUNCTION);
         }
         break;
       case MAIN_SOLVER::FEM_LES:

SU2_CFD/include/fluid/CSU2TCLib.hpp

@@ -70,15 +70,33 @@ class CSU2TCLib : public CNEMOGas {
   std::array<su2double,1> Sk_ref; /*!< \brief Vector containing Sutherland's constant for thermal conductivities */
 
   const su2double T_ref_suth = 273.15; /*!<\brief Reference temperature for Sutherland's model [K] */
+  
+  // Coulomb potential constant values source: Scalabrin, NUMERICAL SIMULATION OF WEAKLY IONIZED HYPERSONIC
+  // FLOW OVER REENTRY CAPSULES, 2007.
+  /*--- Attractive Coulombic potential constants ---*/
+  const su2double D1_a = 0.784;
+  const su2double C1_a = -0.476;
+  const su2double c1_a = 0.0313;
+  const su2double D2_a = 1.262;
+  const su2double C2_a = -0.146;
+  const su2double c2_a = 0.0377;
+
+  /*--- Repulsive Coulombic potential constants ---*/
+  const su2double D1_r = 0.765;
+  const su2double C1_r = 0.138;
+  const su2double c1_r = 0.0106;
+  const su2double D2_r = 1.235;
+  const su2double C2_r = 0.157;
+  const su2double c2_r = 0.0274;
 
   su2activematrix CharElTemp,    /*!< \brief Characteristic temperature of electron states. */
   ElDegeneracy,                  /*!< \brief Degeneracy of electron states. */
   RxnConstantTable,              /*!< \brief Table of chemical equiibrium reaction constants */
   Blottner,                      /*!< \brief Blottner viscosity coefficients */
   Dij;                           /*!< \brief Binary diffusion coefficients. */
 
-  C3DDoubleMatrix Omega00,       /*!< \brief Collision integrals (Omega(0,0)) */
-  Omega11;                       /*!< \brief Collision integrals (Omega(1,1)) */
+  C3DDoubleMatrix Omega11,       /*!< \brief Collision integrals (Omega^(1,1)) */
+  Omega22;                       /*!< \brief Collision integrals (Omega^(2,2)) */
 
   /*--- Implicit variables ---*/
   su2double                     /*!< \brief Derivatives w.r.t. conservative variables */
@@ -249,6 +267,29 @@ class CSU2TCLib : public CNEMOGas {
    */
   void ThermalConductivitiesSuth();
 
+  /*!
+   *\brief Compute the collision terms (deltas)
+   *\param[in] iSpecies - Species of gas
+   *\param[in] jSpecies - Collision partner of species i
+   *\param[in] Mi - Molar mass of species i
+   *\param[in] Mj - Molar mass of species j
+   *\param[in] T - Temperature of gas
+   *\param[in] d1 - Whether delta_1 or delta_2 is computed
+   *\param[out] Delta_1 or Delta_2 - Collision term between species i and j
+   */
+  su2double ComputeCollisionDelta(unsigned iSpecies, unsigned jSpecies, su2double Mi, su2double Mj, su2double T, bool d1);
+  
+  /*!
+   *\brief Compute the collision cross section
+   *\param[in] iSpecies - Species of gas
+   *\param[in] jSpecies - Collision partner of species i
+   *\param[in] T - Temperature of gas
+   *\param[in] d1 - Whether omega^(1,1) or omega^(2,2)
+   *\param[in] coulomb - Whether to use Coulomb potential
+   *\param[out] Omega_ij - collision cross section
+   */
+  su2double ComputeCollisionCrossSection(unsigned iSpecies, unsigned jSpecies, su2double T, bool d1, bool coulomb);
+
   /*!
    * \brief Get reference temperature.
    */

SU2_CFD/src/fluid/CMutationTCLib.cpp

@@ -58,43 +58,38 @@ CMutationTCLib::CMutationTCLib(const CConfig* config, unsigned short val_nDim):
   /* Initialize mixture object */
   mix.reset(new Mutation::Mixture(opt));
 
-  for(iSpecies = 0; iSpecies < nSpecies; iSpecies++) MolarMass[iSpecies] = 1000* mix->speciesMw(iSpecies); // x1000 to have Molar Mass in kg/kmol
-
-  if (mix->hasElectrons()) {
-    if (config->GetViscous()) {
-      SU2_MPI::Error("Ionization is not yet operational for a viscous flow in the NEMO solver.", CURRENT_FUNCTION);
-    } else {
-      nHeavy = nSpecies-1;
-      nEl = 1;
-    }
-  }
-  else { nHeavy = nSpecies;   nEl = 0; }
+  // x1000 to have Molar Mass in kg/kmol
+  for(iSpecies = 0; iSpecies < nSpecies; iSpecies++)
+    MolarMass[iSpecies] = 1000* mix->speciesMw(iSpecies);
+
+  if (mix->hasElectrons()) { nHeavy = nSpecies-1; nEl = 1; }
+  else { nHeavy = nSpecies; nEl = 0; }
 
   /*--- Set up catalytic recombination table. ---*/
   // Creation/Destruction (+1/-1), Index of monoatomic reactants
-  // Monoatomic species (N,O) recombine into diaatomic (N2, O2)
+  // Monoatomic species (N,O) recombine into diaatomic (N2, O2) species
   if (gas_model == "N2") {
-    CatRecombTable(0,0) =  1; CatRecombTable(0,1) = 1;
-    CatRecombTable(1,0) = -1; CatRecombTable(1,1) = 1;
+    CatRecombTable(0,0) =  1; CatRecombTable(0,1) = 1; // N2
+    CatRecombTable(1,0) = -1; CatRecombTable(1,1) = 1; // N
 
   } else if (gas_model == "air_5"){
-    CatRecombTable(0,0) = -1; CatRecombTable(0,1) = 0;
-    CatRecombTable(1,0) = -1; CatRecombTable(1,1) = 1;
-    CatRecombTable(2,0) =  0; CatRecombTable(2,1) = 4;
-    CatRecombTable(3,0) =  1; CatRecombTable(3,1) = 0;
-    CatRecombTable(4,0) =  1; CatRecombTable(4,1) = 1;
+    CatRecombTable(0,0) = -1; CatRecombTable(0,1) = 0; // N
+    CatRecombTable(1,0) = -1; CatRecombTable(1,1) = 1; // O
+    CatRecombTable(2,0) =  0; CatRecombTable(2,1) = 4; // NO
+    CatRecombTable(3,0) =  1; CatRecombTable(3,1) = 0; // N2
+    CatRecombTable(4,0) =  1; CatRecombTable(4,1) = 1; // O2
 
   } else if (gas_model == "air_6") {
-    CatRecombTable(0,0) = -1; CatRecombTable(0,1) = 0;
-    CatRecombTable(1,0) = -1; CatRecombTable(1,1) = 1;
-    CatRecombTable(2,0) =  0; CatRecombTable(2,1) = 4;
-    CatRecombTable(3,0) =  1; CatRecombTable(3,1) = 0;
-    CatRecombTable(4,0) =  1; CatRecombTable(4,1) = 1;
-    CatRecombTable(5,0) =  0; CatRecombTable(5,1) = 4;
+    CatRecombTable(0,0) = -1; CatRecombTable(0,1) = 0; // N
+    CatRecombTable(1,0) = -1; CatRecombTable(1,1) = 1; // O
+    CatRecombTable(2,0) =  0; CatRecombTable(2,1) = 4; // NO
+    CatRecombTable(3,0) =  1; CatRecombTable(3,1) = 0; // N2
+    CatRecombTable(4,0) =  1; CatRecombTable(4,1) = 1; // O2
+    CatRecombTable(5,0) =  0; CatRecombTable(5,1) = 4; // Ar
 
   } else {
     if (config->GetCatalytic())
-      SU2_MPI::Error("Cataylic wall recombination not implemented for specified Mutation gas model.", CURRENT_FUNCTION);
+      SU2_MPI::Error("Catalytic wall recombination not implemented for specified Mutation gas model.", CURRENT_FUNCTION);
   }
 
 }