Add in Sutherland's law for NEMO problems + some clean up to viscous functions

Author: WallyMaier

Common/include/option_structure.hpp

@@ -613,11 +613,13 @@ MakePair("ONESPECIES", ONESPECIES)
  * \brief types of coefficient transport model
  */
 enum class TRANSCOEFFMODEL {
+  SUTHERLAND,
   WILKE,
   GUPTAYOS,
   CHAPMANN_ENSKOG
 };
 static const MapType<std::string, TRANSCOEFFMODEL> TransCoeffModel_Map = {
+MakePair("SUTHERLAND", TRANSCOEFFMODEL::SUTHERLAND)
 MakePair("WILKE", TRANSCOEFFMODEL::WILKE)
 MakePair("GUPTA-YOS", TRANSCOEFFMODEL::GUPTAYOS)
 MakePair("CHAPMANN-ENSKOG", TRANSCOEFFMODEL::CHAPMANN_ENSKOG)

Common/src/CConfig.cpp

@@ -3792,8 +3792,8 @@ 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 ) {
-    SU2_MPI::Error("Only WILKE transport model is stable for the NEMO solver using SU2TClib. Use Mutation++ instead.", 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);
   }
 
   if (Kind_FluidModel == MUTATIONPP && (Kind_TransCoeffModel != TRANSCOEFFMODEL::WILKE && Kind_TransCoeffModel != TRANSCOEFFMODEL::CHAPMANN_ENSKOG)) {

SU2_CFD/include/fluid/CSU2TCLib.hpp

@@ -64,6 +64,12 @@ class CSU2TCLib : public CNEMOGas {
   phis, mus,                      /*!< \brief Auxiliary vectors to be used in Wilke/Blottner/Eucken model */
   A;                              /*!< \brief Auxiliary vector to be used in net production rate computation */
 
+  vector<su2double>
+  mu_ref,                         /*!< \brief Vector containing reference viscosity for Sutherland's law */
+  k_ref,                          /*!< \brief Vector containing reference thermal conducivities for Sutherland's law */
+  Sm_ref,                         /*!< \brief Vector containing Sutherland's constant for viscosity */
+  Sk_ref;                         /*!< \brief Vector containing Sutherland's constant for thermal conductivities */
+
   su2activematrix CharElTemp,    /*!< \brief Characteristic temperature of electron states. */
   ElDegeneracy,                  /*!< \brief Degeneracy of electron states. */
   RxnConstantTable,              /*!< \brief Table of chemical equiibrium reaction constants */
@@ -232,6 +238,16 @@ class CSU2TCLib : public CNEMOGas {
    */
   void ThermalConductivitiesGY();
 
+  /*!
+   * \brief Get viscosity with Sutherland's transport model.
+   */
+  void ViscositySuth();
+
+  /*!
+   * \brief Get T-R and V-E thermal conductivities vector with Sutherland's transport model.
+   */
+  void ThermalConductivitiesSuth();
+
   /*!
    * \brief Get reference temperature.
    */

SU2_CFD/src/fluid/CSU2TCLib.cpp

@@ -52,10 +52,17 @@ CSU2TCLib::CSU2TCLib(const CConfig* config, unsigned short val_nDim, bool viscou
   eve_eq.resize(nSpecies,0.0);
   eve.resize(nSpecies,0.0);
 
-  if(viscous){
+  if (viscous) {
     MolarFracWBE.resize(nSpecies,0.0);
     phis.resize(nSpecies,0.0);
     mus.resize(nSpecies,0.0);
+
+    /*--- Reference values for Sutherland's law viscosity ---*/
+    //Note: Can be extended with better coefficients of individual species,.
+    k_ref.resize(1,0.0);
+    mu_ref.resize(1,0.0);
+    Sm_ref.resize(1,0.0);
+    Sk_ref.resize(1,0.0);
   }
 
   if (gas_model =="ARGON"){
@@ -110,6 +117,13 @@ CSU2TCLib::CSU2TCLib(const CConfig* config, unsigned short val_nDim, bool viscou
     ElDegeneracy(0,5) = 5;
     ElDegeneracy(0,6) = 15;
 
+    if (viscous) {
+      mu_ref[0] = 2.125E-5;
+      k_ref[0] = 0.0163;
+      Sm_ref[0] = 114.0;
+      Sk_ref[0] = 170;
+    }
+
   } else if (gas_model == "N2"){
     /*--- Check for errors in the initialization ---*/
     if (nSpecies != 2) {
@@ -252,6 +266,13 @@ CSU2TCLib::CSU2TCLib(const CConfig* config, unsigned short val_nDim, bool viscou
     Omega11(1,0,0) = -8.3493693E-03;  Omega11(1,0,1) = 1.7808911E-01;   Omega11(1,0,2) = -1.4466155E+00;  Omega11(1,0,3) = 1.9324210E+03;
     Omega11(1,1,0) = -7.7439615E-03;  Omega11(1,1,1) = 1.7129007E-01;   Omega11(1,1,2) = -1.4809088E+00;  Omega11(1,1,3) = 2.1284951E+03;
 
+    if (viscous) {
+      k_ref[0] = 0.0242;
+      mu_ref[0] = 1.663E-5;
+      Sm_ref[0] = 107.0;
+      Sk_ref[0] = 150.0;
+    }
+
   } else if (gas_model == "AIR-5"){
 
     /*--- Check for errors in the initialization ---*/
@@ -598,6 +619,13 @@ CSU2TCLib::CSU2TCLib(const CConfig* config, unsigned short val_nDim, bool viscou
     Omega11(4,3,0) = -5.0478143E-03;  Omega11(4,3,1) = 1.0236186E-01;   Omega11(4,3,2) = -9.0058935E-01;  Omega11(4,3,3) = 4.4472565E+02;
     Omega11(4,4,0) = -4.2451096E-03;  Omega11(4,4,1) = 9.6820337E-02;   Omega11(4,4,2) = -9.9770795E-01;  Omega11(4,4,3) = 8.3320644E+02;
 
+    if (viscous) {
+      k_ref[0] = 0.0241;
+      mu_ref[0] = 1.716E-5;
+      Sm_ref[0] = 111.0;
+      Sk_ref[0] = 194.0;
+    }
+
   } else if (gas_model == "AIR-7"){
 
     /*--- Check for errors in the initialization ---*/
@@ -1027,6 +1055,13 @@ CSU2TCLib::CSU2TCLib(const CConfig* config, unsigned short val_nDim, bool viscou
     Omega11(4,2,0) = -1.0066279E-03;  Omega11(4,2,1) = 1.1029264E-02;   Omega11(4,2,2) = -2.0671266E-01;  Omega11(4,2,3) = 8.2644384E+01;
     Omega11(4,3,0) = -5.0478143E-03;  Omega11(4,3,1) = 1.0236186E-01;   Omega11(4,3,2) = -9.0058935E-01;  Omega11(4,3,3) = 4.4472565E+02;
     Omega11(4,4,0) = -4.2451096E-03;  Omega11(4,4,1) = 9.6820337E-02;   Omega11(4,4,2) = -9.9770795E-01;  Omega11(4,4,3) = 8.3320644E+02;
+
+    if (viscous) {
+      k_ref[0] = 0.0241;
+      mu_ref[0] = 1.716E-5;
+      Sm_ref[0] = 111.0;
+      Sk_ref[0] = 194.0;
+    }
   }
 
   if (ionization) { nHeavy = nSpecies-1; nEl = 1; }
@@ -1598,6 +1633,8 @@ vector<su2double>& CSU2TCLib::GetDiffusionCoeff(){
    DiffusionCoeffWBE();
   if(Kind_TransCoeffModel == TRANSCOEFFMODEL::GUPTAYOS)
    DiffusionCoeffGY();
+  if(Kind_TransCoeffModel == TRANSCOEFFMODEL::SUTHERLAND)
+   DiffusionCoeffWBE();
 
   return DiffusionCoeff;
 
@@ -1609,6 +1646,8 @@ su2double CSU2TCLib::GetViscosity(){
     ViscosityWBE();
   if(Kind_TransCoeffModel == TRANSCOEFFMODEL::GUPTAYOS)
     ViscosityGY();
+  if(Kind_TransCoeffModel == TRANSCOEFFMODEL::SUTHERLAND)
+    ViscositySuth();
 
   return Mu;
 
@@ -1620,6 +1659,8 @@ vector<su2double>& CSU2TCLib::GetThermalConductivities(){
     ThermalConductivitiesWBE();
   if(Kind_TransCoeffModel == TRANSCOEFFMODEL::GUPTAYOS)
     ThermalConductivitiesGY();
+  if(Kind_TransCoeffModel == TRANSCOEFFMODEL::SUTHERLAND)
+    ThermalConductivitiesSuth();
 
   return ThermalConductivities;
 
@@ -1802,7 +1843,8 @@ void CSU2TCLib::DiffusionCoeffGY(){
     //}
 
     /*--- Assign species diffusion coefficient ---*/
-    DiffusionCoeff[iSpecies] = gam_t*gam_t*Mi*(1-Mi*gam_i) / denom;
+    if (denom == 0) { DiffusionCoeff[iSpecies] = 0.0; }
+    else { DiffusionCoeff[iSpecies] = gam_t*gam_t*Mi*(1-Mi*gam_i) / denom; }
   }
   // if (ionization) {
   //TODO: Update correct iElectron....
@@ -1971,18 +2013,50 @@ void CSU2TCLib::ThermalConductivitiesGY(){
     }
 
     /*--- Translational contribution to thermal conductivity ---*/
-    ThermalCond_tr    += (15.0/4.0)*kb*gam_i/denom_t;
+    if (denom_t != 0) ThermalCond_tr += (15.0/4.0)*kb*gam_i/denom_t;
 
     /*--- Translational contribution to thermal conductivity ---*/
-    if (RotationModes[iSpecies] != 0.0)
-      ThermalCond_tr  += kb*gam_i/denom_r;
+    if (RotationModes[iSpecies] != 0.0 && denom_r != 0) ThermalCond_tr  += kb*gam_i/denom_r;
 
     /*--- Vibrational-electronic contribution to thermal conductivity ---*/
-    ThermalCond_ve += kb*Cvve/R*gam_i / denom_r;
+    if (denom_r != 0) ThermalCond_ve += kb*Cvve/R*gam_i / denom_r;
   }
 
   ThermalConductivities[0] = ThermalCond_tr;
   ThermalConductivities[1] = ThermalCond_ve;
+}
+
+void CSU2TCLib::ViscositySuth(){
+
+  su2double T_ref = 273.15;
+  su2double T_nd = T / T_ref;
+
+  /*--- Calculate mixture laminar viscosity ---*/
+  Mu = mu_ref[0] * T_nd * sqrt(T_nd) * ((T_ref + Sm_ref[0]) / (T + Sm_ref[0]));
+
+}
+
+void CSU2TCLib::ThermalConductivitiesSuth(){
+
+  /*--- Compute mixture quantities ---*/
+  su2double mass = 0.0, rho = 0.0;
+  for (unsigned short ii=0; ii<nSpecies; ii++) rho  += rhos[ii];
+  for (unsigned short ii=0; ii<nSpecies; ii++) mass += rhos[ii]/rho*MolarMass[ii];
+
+  su2double Cvtr = ComputerhoCvtr()/rho;
+  su2double Cvve = ComputerhoCvve()/rho;
+
+  /*--- Compute simple Kve scaling factor ---*/
+  su2double scl  = Cvve/Cvtr;
+
+  /*--- Compute k's using Sutherland's law ---*/
+  su2double T_ref = 273.15;
+  su2double T_nd = T / T_ref;
+  su2double k = k_ref[0] * T_nd * sqrt(T_nd) * ((T_ref + Sk_ref[0]) / (T + Sk_ref[0]));
+  su2double kve = scl*k;
+
+  ThermalConductivities[0] = k;
+  ThermalConductivities[1] = kve;
 
 }
 

SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp

@@ -228,61 +228,58 @@ void CNEMONumerics::GetViscousProjFlux(const su2double *val_primvar,
                                        su2double val_therm_conductivity_ve,
                                        const CConfig *config) {
 
+  if (ionization) {
+    SU2_MPI::Error("NEED TO IMPLEMENT IONIZED FUNCTIONALITY!!!",CURRENT_FUNCTION);
+  }
+
   // Requires a slightly non-standard primitive vector:
   // Assumes -     V = [Y1, ... , Yn, T, Tve, ... ]
   // and gradient GV = [GY1, ... , GYn, GT, GTve, ... ]
   // rather than the standard V = [r1, ... , rn, T, Tve, ... ]
 
-  unsigned short iSpecies, iVar, iDim, jDim;
-  su2double mu, ktr, kve, rho, T, Tve, RuSI, Ru;
-  auto& Ms = fluidmodel->GetSpeciesMolarMass();
-
   su2activematrix Flux_Tensor(nVar,nDim);
 
   /*--- Initialize ---*/
-  for (iVar = 0; iVar < nVar; iVar++) {
+  for (auto iVar = 0; iVar < nVar; iVar++) {
     Proj_Flux_Tensor[iVar] = 0.0;
-    for (iDim = 0; iDim < nDim; iDim++)
+    for (auto iDim = 0; iDim < nDim; iDim++)
       Flux_Tensor[iVar][iDim] = 0.0;
   }
 
-  /*--- Rename for convenience ---*/
+  /*--- Rename variables for convenience ---*/
+  const auto& Ms = fluidmodel->GetSpeciesMolarMass();
   const auto& Ds  = val_diffusioncoeff;
-  mu  = val_lam_viscosity+val_eddy_viscosity;
-  ktr = val_therm_conductivity;
-  kve = val_therm_conductivity_ve;
-  rho = val_primvar[RHO_INDEX];
-  T   = val_primvar[T_INDEX];
-  Tve = val_primvar[TVE_INDEX];
+  const su2double mu  = val_lam_viscosity+val_eddy_viscosity;
+  su2double ktr = val_therm_conductivity;
+  su2double kve = val_therm_conductivity_ve;
+  const su2double rho = val_primvar[RHO_INDEX];
+  const su2double T = val_primvar[T_INDEX];
+  const su2double Tve = val_primvar[TVE_INDEX];
   const auto& V   = val_primvar;
   const auto& GV  = val_gradprimvar;
-  RuSI= UNIVERSAL_GAS_CONSTANT;
-  Ru  = 1000.0*RuSI;
-
+  const su2double Ru = 1000.0*UNIVERSAL_GAS_CONSTANT;
   const auto& hs = fluidmodel->ComputeSpeciesEnthalpy(T, Tve, val_eve);
 
   /*--- Scale thermal conductivity with turb visc ---*/
   // TODO: Need to determine proper way to incorporate eddy viscosity
   // This is only scaling Kve by same factor as ktr
   // NOTE: V[iSpecies] is == Ys.
   su2double Mass = 0.0;
-  su2double tmp1, scl, Cptr;
-  for (iSpecies=0;iSpecies<nSpecies;iSpecies++)
+  for (auto iSpecies = 0;iSpecies<nSpecies;iSpecies++)
     Mass += V[iSpecies]*Ms[iSpecies];
-  Cptr = V[RHOCVTR_INDEX]/V[RHO_INDEX]+Ru/Mass;
-  tmp1 = Cptr*(val_eddy_viscosity/Prandtl_Turb);
-  scl  = tmp1/ktr;
+
+  su2double Cptr = V[RHOCVTR_INDEX]/V[RHO_INDEX]+Ru/Mass;
+  su2double tmp1 = Cptr*(val_eddy_viscosity/Prandtl_Turb);
+  su2double scl = tmp1/ktr;
   ktr += Cptr*(val_eddy_viscosity/Prandtl_Turb);
   kve  = kve*(1.0+scl);
   //Cpve = V[RHOCVVE_INDEX]+Ru/Mass;
   //kve += Cpve*(val_eddy_viscosity/Prandtl_Turb);
 
-  /*--- Pre-compute mixture quantities ---*/
-
+  /*--- Pre-compute mixture quantities ---*/  //TODO
   su2double Vector[MAXNDIM] = {0.0};
-
-  for (iDim = 0; iDim < nDim; iDim++) {
-    for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
+  for (auto iDim = 0; iDim < nDim; iDim++) {
+    for (auto iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
       Vector[iDim] += rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][iDim];
     }
   }
@@ -291,37 +288,34 @@ void CNEMONumerics::GetViscousProjFlux(const su2double *val_primvar,
   ComputeStressTensor(nDim,tau,val_gradprimvar+VEL_INDEX, mu);
 
   /*--- Populate entries in the viscous flux vector ---*/
-  for (iDim = 0; iDim < nDim; iDim++) {
+  for (auto iDim = 0; iDim < nDim; iDim++) {
+
     /*--- Species diffusion velocity ---*/
-    for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
+    for (auto iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
       Flux_Tensor[iSpecies][iDim] = rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][iDim]
           - V[RHOS_INDEX+iSpecies]*Vector[iDim];
     }
-    if (ionization) {
-      SU2_MPI::Error("NEED TO IMPLEMENT IONIZED FUNCTIONALITY!!!",CURRENT_FUNCTION);
-    }
 
-    /*--- Shear stress related terms ---*/
+    /*--- Shear-stress/momentum related terms ---*/
     Flux_Tensor[nSpecies+nDim][iDim] = 0.0;
-    for (jDim = 0; jDim < nDim; jDim++) {
+    for (auto jDim = 0; jDim < nDim; jDim++) {
       Flux_Tensor[nSpecies+jDim][iDim]  = tau[iDim][jDim];
       Flux_Tensor[nSpecies+nDim][iDim] += tau[iDim][jDim]*val_primvar[VEL_INDEX+jDim];
     }
 
     /*--- Diffusion terms ---*/
-    for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
+    for (auto iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
       Flux_Tensor[nSpecies+nDim][iDim]   += Flux_Tensor[iSpecies][iDim] * hs[iSpecies];
       Flux_Tensor[nSpecies+nDim+1][iDim] += Flux_Tensor[iSpecies][iDim] * val_eve[iSpecies];
     }
 
     /*--- Heat transfer terms ---*/
-    Flux_Tensor[nSpecies+nDim][iDim]   += ktr*GV[T_INDEX][iDim] +
-        kve*GV[TVE_INDEX][iDim];
+    Flux_Tensor[nSpecies+nDim][iDim] += ktr*GV[T_INDEX][iDim] + kve*GV[TVE_INDEX][iDim];
     Flux_Tensor[nSpecies+nDim+1][iDim] += kve*GV[TVE_INDEX][iDim];
   }
 
-  for (iVar = 0; iVar < nVar; iVar++) {
-    for (iDim = 0; iDim < nDim; iDim++) {
+  for (auto iVar = 0; iVar < nVar; iVar++) {
+    for (auto iDim = 0; iDim < nDim; iDim++) {
       Proj_Flux_Tensor[iVar] += Flux_Tensor[iVar][iDim]*val_normal[iDim];
     }
   }

SU2_CFD/src/numerics/NEMO/NEMO_diffusion.cpp

@@ -128,14 +128,10 @@ CNumerics::ResidualType<> CAvgGrad_NEMO::ComputeResidual(const CConfig *config)
     PrimVar_j[iSpecies] = V_j[iSpecies]/V_j[RHO_INDEX];
     Mean_PrimVar[iSpecies] = 0.5*(PrimVar_i[iSpecies] + PrimVar_j[iSpecies]);
     for (auto iDim = 0; iDim < nDim; iDim++) {
-      Mean_GradPrimVar[iSpecies][iDim] = 0.5*(1.0/V_i[RHO_INDEX] *
-                                              (PrimVar_Grad_i[iSpecies][iDim] -
-                                               PrimVar_i[iSpecies] *
-                                               PrimVar_Grad_i[RHO_INDEX][iDim]) +
-                                              1.0/V_j[RHO_INDEX] *
-                                              (PrimVar_Grad_j[iSpecies][iDim] -
-                                               PrimVar_j[iSpecies] *
-                                               PrimVar_Grad_j[RHO_INDEX][iDim]));
+      Mean_GradPrimVar[iSpecies][iDim] = 0.5*(1.0/V_i[RHO_INDEX] * (PrimVar_Grad_i[iSpecies][iDim] -
+                                              PrimVar_i[iSpecies] * PrimVar_Grad_i[RHO_INDEX][iDim]) +
+                                              1.0/V_j[RHO_INDEX] * (PrimVar_Grad_j[iSpecies][iDim] -
+                                              PrimVar_j[iSpecies] * PrimVar_Grad_j[RHO_INDEX][iDim]));
     }
   }
 

SU2_CFD/src/solvers/CNEMOEulerSolver.cpp

@@ -1312,6 +1312,11 @@ void CNEMOEulerSolver::SetNondimensionalization(CConfig *config, unsigned short
     if (viscous) {
 
       switch(config->GetKind_TransCoeffModel()){
+      case TRANSCOEFFMODEL::SUTHERLAND:
+        ModelTable << "Sutherland + WBE";
+        NonDimTable.PrintFooter();
+        break;
+
       case TRANSCOEFFMODEL::WILKE:
         ModelTable << "Wilke-Blottner-Eucken";
         NonDimTable.PrintFooter();

config_template.cfg

@@ -1346,7 +1346,7 @@ CFL_REDUCTION_ADJTURB= 0.01
 
 % -------------------- NEMO NUMERICAL METHOD DEFINITION -----------------------%
 %
-% Mixture transport properties (WILKE,GUPTA-YOS,CHAPMANN-ENSKOG)
+% Mixture transport properties (WILKE,GUPTA-YOS,CHAPMANN-ENSKOG, SUTHERLAND)
 TRANSPORT_COEFF_MODEL = WILKE
 
 % ----------------------- GEOMETRY EVALUATION PARAMETERS ----------------------%