[WIP] Implementation of eN transition model

Common/include/CConfig.hpp

@@ -866,6 +866,9 @@ class CConfig {
   Tke_FreeStream,                  /*!< \brief Total turbulent kinetic energy of the fluid.  */
   Intermittency_FreeStream,        /*!< \brief Freestream intermittency (for sagt transition model) of the fluid.  */
   ReThetaT_FreeStream,             /*!< \brief Freestream Transition Momentum Thickness Reynolds Number (for LM transition model) of the fluid.  */
+  TurbulenceIntensity_FreeStream,  /*!< \brief Freestream turbulent intensity (for sagt transition model) of the fluid.  */
+  AmplificationFactor_FreeStream,  /*!< \brief Freestream amplifictation factor for the eN 1 equation transition model.  */
+  Turb2LamViscRatio_FreeStream,    /*!< \brief Ratio of turbulent to laminar viscosity. */
   NuFactor_FreeStream,             /*!< \brief Ratio of turbulent to laminar viscosity. */
   NuFactor_Engine,                 /*!< \brief Ratio of turbulent to laminar viscosity at the engine. */
   SecondaryFlow_ActDisk,           /*!< \brief Ratio of turbulent to laminar viscosity at the actuator disk. */
@@ -1959,6 +1962,12 @@ class CConfig {
    * \return Freestream momentum thickness Reynolds number.
    */
   su2double GetReThetaT_FreeStream() const { return ReThetaT_FreeStream; }
+
+  /*!
+   * \brief Get the value of the non-dimensionalized freestream amplification factor.
+   * \return Non-dimensionalized freestream amplification factor.
+   */
+  su2double GetAmplificationFactor_FreeStream(void) const { return AmplificationFactor_FreeStream; }
 
   /*!
    * \brief Get the value of the non-dimensionalized freestream turbulence intensity.

Common/include/option_structure.hpp

@@ -1161,10 +1161,12 @@ inline SA_ParsedOptions ParseSAOptions(const SA_OPTIONS *SA_Options, unsigned sh
 enum class TURB_TRANS_MODEL {
   NONE,  /*!< \brief No transition model. */
   LM,    /*!< \brief Kind of transition model (Langtry-Menter (LM) for SST and Spalart-Allmaras). */
+  EN,	 /*!< \brief Kind of transition model using Amplification factor purely for Spalart-Allmaras*/
 };
 static const MapType<std::string, TURB_TRANS_MODEL> Trans_Model_Map = {
   MakePair("NONE", TURB_TRANS_MODEL::NONE)
   MakePair("LM", TURB_TRANS_MODEL::LM)
+  MakePair("EN", TURB_TRANS_MODEL::EN)
 };
 
 /*!

Common/src/CConfig.cpp

@@ -4797,7 +4797,17 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     for (int i=0; i<7; ++i) eng_cyl[i] /= 12.0;
   }
 
-  if(Turb_Fixed_Values && !OptionIsSet("TURB_FIXED_VALUES_DOMAIN")){
+  if (Kind_Trans_Model == TURB_TRANS_MODEL::EN) {
+    if (saParsedOptions.bc){
+      SU2_MPI::Error("Please select only one transition model for the SA turbulence model (Choose either eN or BCM).", CURRENT_FUNCTION);
+    } else if (!saParsedOptions.ft2) {
+      SU2_MPI::Error("Please select SA_OPTIONS = WITHFT2 when using SA-Ft2-eN transtion model.", CURRENT_FUNCTION);
+    } else if ((Kind_Regime == ENUM_REGIME::COMPRESSIBLE) && (Ref_NonDim == 0)) {
+      SU2_MPI::Error("Please select a non-dimensionalization option other than 'REF_DIMENSIONALIZATION = DIMENSIONAL' when using SA-Ft2-eN for a compressible case.", CURRENT_FUNCTION);
+    }
+  }
+
+  if (Turb_Fixed_Values && !OptionIsSet("TURB_FIXED_VALUES_DOMAIN")){
     SU2_MPI::Error("TURB_FIXED_VALUES activated, but no domain set with TURB_FIXED_VALUES_DOMAIN.", CURRENT_FUNCTION);
   }
 
@@ -6073,6 +6083,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
         }
         switch (Kind_Trans_Model) {
           case TURB_TRANS_MODEL::NONE:  break;
+          case TURB_TRANS_MODEL::EN:    cout << "Low-turbulence Transition model: eN 1 equation model (2014)" << endl; break;
           case TURB_TRANS_MODEL::LM: {
             cout << "Transition model: Langtry and Menter's 4 equation model";
             if (lmParsedOptions.LM2015) {
@@ -6103,6 +6114,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
               break;
           }
         }
+
         cout << "Hybrid RANS/LES: ";
         switch (Kind_HybridRANSLES) {
           case NO_HYBRIDRANSLES: cout << "No Hybrid RANS/LES" << endl; break;

SU2_CFD/include/drivers/CDriver.hpp

@@ -233,6 +233,7 @@ class CDriver : public CDriverBase {
   template <class FlowIndices>
   void InstantiateTransitionNumerics(unsigned short nVar_Trans, int offset, const CConfig* config,
                                      const CSolver* trans_solver, CNumerics****& numerics) const;
+
   /*!
    * \brief Helper to instantiate species transport numerics specialized for different flow solvers.
    */

SU2_CFD/include/numerics/CNumerics.hpp

@@ -80,6 +80,8 @@ class CNumerics {
   turb_ke_i,  /*!< \brief Turbulent kinetic energy at point i. */
   turb_ke_j;  /*!< \brief Turbulent kinetic energy at point j. */
   su2double
+  amplification_factor_i;  /*!< \brief amplification factor at point i. */
+  su2double
   intermittency_eff_i, /*!< \brief effective intermittency at point i. */
   intermittency_i; /*!< \brief intermittency at point i. */
   su2double
@@ -709,6 +711,19 @@ class CNumerics {
    */
   virtual void SetCrossDiff(su2double val_CDkw_i) {/* empty */};
 
+  /*!
+  * \brief Get the Amplification factor for the e^N model.
+  */
+  inline su2double GetAmplificationFactor() const { return amplification_factor_i; }
+
+  /*!
+   * \brief Set the value of the amplification factor for the e^N model.
+   * \param[in] amplification_factor_i - Value of the amplification factor at point i.
+   */
+  void SetAmplificationFactor(su2double val_amplification_factor_i) {
+    amplification_factor_i = val_amplification_factor_i;
+  };
+
   /*!
    * \brief Set the value of the effective intermittency for the LM model.
    * \param[in] intermittency_eff_i - Value of the effective intermittency at point i.

SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp

@@ -2,7 +2,7 @@
  * \file trans_convection.hpp
  * \brief Delarations of numerics classes for discretization of
  *        convective fluxes in transition problems.
- * \author S. Kang
+ * \author S. Kang, R. Roos
  * \version 7.5.1 "Blackbird"
  *
  * SU2 Project Website: https://su2code.github.io
@@ -27,8 +27,8 @@
  */
 
 #pragma once
-
 #include "../turb_convection.hpp"
+#include "../../scalar/scalar_convection.hpp"
 
 /*!
  * \class CUpwSca_TransLM
@@ -38,3 +38,11 @@
 template <class FlowIndices>
 using CUpwSca_TransLM  = CUpwSca_TurbSST<FlowIndices>;
 
+/*!
+ * \class CUpwSca_TransLM
+ * \brief Re-use the SA convective fluxes for the scalar upwind discretization of eN transition model equations.
+ * \ingroup ConvDiscr
+ */
+template <class FlowIndices>
+using CUpwSca_TransEN  = CUpwSca_TurbSA<FlowIndices>;
+

SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp

@@ -109,7 +109,7 @@ class TransLMCorrelations {
         break;
       }
       case TURB_TRANS_CORRELATION::DEFAULT:
-        SU2_MPI::Error("Transition correlation is set to DEFAULT but no default value has ben set in the code.",
+        SU2_MPI::Error("Transition correlation is set to DEFAULT but no default value has been set in the code.",
                        CURRENT_FUNCTION);
         break;
     }
@@ -182,7 +182,7 @@ class TransLMCorrelations {
         break;
       }
       case TURB_TRANS_CORRELATION::DEFAULT:
-        SU2_MPI::Error("Transition correlation is set to DEFAULT but no default value has ben set in the code.",
+        SU2_MPI::Error("Transition correlation is set to DEFAULT but no default value has been set in the code.",
                        CURRENT_FUNCTION);
         break;
     }

SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp

@@ -2,7 +2,7 @@
  * \file trans_diffusion.hpp
  * \brief Declarations of numerics classes for discretization of
  *        viscous fluxes in transition problems.
- * \author S. Kang
+ * \author S. Kang, R. Roos
  * \version 7.5.1 "Blackbird"
  *
  * SU2 Project Website: https://su2code.github.io
@@ -103,3 +103,72 @@ class CAvgGrad_TransLM final : public CAvgGrad_Scalar<FlowIndices> {
   }
 
 };
+
+/*!
+ * \class CAvgGrad_TransEN
+ * \brief Class for computing viscous term using average of gradient with correction (e^N transition model).
+ * \ingroup ViscDiscr
+ * \author R. Roos
+ */
+template <class FlowIndices>
+class CAvgGrad_TransEN final : public CAvgGrad_Scalar<FlowIndices> {
+private:
+  using Base = CAvgGrad_Scalar<FlowIndices>;
+  using Base::Laminar_Viscosity_i;
+  using Base::Laminar_Viscosity_j;
+  using Base::Eddy_Viscosity_i;
+  using Base::Eddy_Viscosity_j;
+  using Base::Density_i;
+  using Base::Density_j;
+  using Base::ScalarVar_i;
+  using Base::ScalarVar_j;
+  using Base::Proj_Mean_GradScalarVar;
+  using Base::proj_vector_ij;
+  using Base::Flux;
+  using Base::Jacobian_i;
+  using Base::Jacobian_j;
+
+  const su2double sigma_n = 1.0;
+
+  /*!
+   * \brief Adds any extra variables to AD
+   */
+  void ExtraADPreaccIn() override {}
+
+  /*!
+   * \brief SA specific steps in the ComputeResidual method
+   * \param[in] config - Definition of the particular problem.
+   */
+  void FinishResidualCalc(const CConfig* config) override {
+	const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT;
+
+	/*--- Compute mean effective dynamic viscosity ---*/
+	const su2double diff_i_amplification = (Laminar_Viscosity_i + Eddy_Viscosity_i)/sigma_n;
+	const su2double diff_j_amplification = (Laminar_Viscosity_j + Eddy_Viscosity_j)/sigma_n;
+
+	const su2double diff_amplification = 0.5*(diff_i_amplification + diff_j_amplification);
+
+    Flux[0] = diff_amplification*Proj_Mean_GradScalarVar[0];
+
+    /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/
+
+    if (implicit) {
+      Jacobian_i[0][0] = (0.5*Proj_Mean_GradScalarVar[0]-diff_amplification*proj_vector_ij);
+      Jacobian_j[0][0] = (0.5*Proj_Mean_GradScalarVar[0]+diff_amplification*proj_vector_ij);
+    }
+  }
+
+public:
+  /*!
+   * \brief Constructor of the class.
+   * \param[in] val_nDim - Number of dimensions of the problem.
+   * \param[in] val_nVar - Number of variables of the problem.
+   * \param[in] correct_grad - Whether to correct gradient for skewness.
+   * \param[in] config - Definition of the particular problem.
+   */
+  CAvgGrad_TransEN(unsigned short val_nDim, unsigned short val_nVar,
+                  bool correct_grad, const CConfig* config)
+    : CAvgGrad_Scalar<FlowIndices>(val_nDim, val_nVar, correct_grad, config) {
+  }
+
+};

SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp

@@ -323,3 +323,142 @@ class CSourcePieceWise_TransLM final : public CNumerics {
     return ResidualType<>(Residual, Jacobian_i, nullptr);
   }
 };
+
+/*!
+ * \class CSourcePieceWise_TranEN
+ * \brief Class for integrating the source terms of the e^N transition model equations.
+ * \ingroup SourceDiscr
+ * \author R. Roos
+ */
+template <class FlowIndices>
+class CSourcePieceWise_TransEN final : public CNumerics {
+ private:
+  const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */
+
+  su2double Residual, *Jacobian_i;
+  su2double Jacobian_Buffer; /*!< \brief Static storage for the Jacobian (which needs to be pointer for return type). */
+
+ public:
+  /*!
+   * \brief Constructor of the class.
+   * \param[in] val_nDim - Number of dimensions of the problem.
+   * \param[in] val_nVar - Number of variables of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  CSourcePieceWise_TransEN(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config)
+      : CNumerics(val_nDim, 1, config),
+        idx(val_nDim, config->GetnSpecies()) {
+
+    /*--- "Allocate" the Jacobian using the static buffer. ---*/
+	Jacobian_i = &Jacobian_Buffer;
+  }
+
+  /*!
+   * \brief Residual for source term integration.
+   * \param[in] config - Definition of the particular problem.
+   * \return A lightweight const-view (read-only) of the residual/flux and Jacobians.
+   */
+  ResidualType<> ComputeResidual(const CConfig* config) override {
+
+    AD::StartPreacc();
+    AD::SetPreaccIn(V_i[idx.Density()], V_i[idx.Pressure()], V_i[idx.LaminarViscosity()], StrainMag_i, ScalarVar_i[0], Volume, dist_i);
+    AD::SetPreaccIn(&V_i[idx.Velocity()], nDim);
+    AD::SetPreaccIn(Vorticity_i, 3);
+    AD::SetPreaccIn(PrimVar_Grad_i + idx.Velocity(), nDim, nDim);
+    AD::SetPreaccIn(ScalarVar_Grad_i[0], nDim);
+
+    su2double rho 			= V_i[idx.Density()];
+    su2double p 			= V_i[idx.Pressure()];
+    su2double muLam 		= V_i[idx.LaminarViscosity()];
+
+    const su2double VorticityMag = GeometryToolbox::Norm(3, Vorticity_i);
+
+    const su2double vel_u 	= V_i[idx.Velocity()];
+    const su2double vel_v 	= V_i[1+idx.Velocity()];
+    const su2double vel_w 	= (nDim ==3) ? V_i[2+idx.Velocity()] : 0.0;
+
+    const su2double vel_mag = sqrt(vel_u*vel_u + vel_v*vel_v + vel_w*vel_w);
+
+    su2double rhoInf		= config->GetDensity_FreeStreamND();
+    su2double pInf 			= config->GetPressure_FreeStreamND();
+    const su2double *velInf = config->GetVelocity_FreeStreamND();
+
+    su2double velInf2 = 0.0;
+    for(unsigned short iDim = 0; iDim < nDim; ++iDim) {
+    	velInf2 += velInf[iDim]*velInf[iDim];
+    }
+
+    Residual = 0.0;
+    Jacobian_i[0] = 0.0;
+
+    if (dist_i > 1e-10) {
+
+      su2double u_e;
+      if (config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) {
+
+  		/*--- Estimate of the equivalent flow velocity at the edge of the boundary layer based on compressible Bernoulli's equation ---*/
+  		const su2double Gamma = config->GetGamma();
+  		const su2double G_over_G = Gamma/(Gamma-1);
+
+  		const su2double rho_e = pow(((pow(rhoInf,Gamma)/pInf)*p),(1/Gamma));
+
+  		/*--- Abs value taken to account of negative difference between pInf and local p ---*/
+  		//u_e = sqrt(2*(G_over_G*( (pInf/rhoInf) -(p/rho_e) )) + velInf2 );
+  		u_e = sqrt(2*(G_over_G*(pInf/rhoInf) + (velInf2/2) - G_over_G*(p/rho_e)));
+
+      } else {
+
+    	/*--- Inviscid edge velocity based on incompressible Bernoulli's equation---*/
+    	//u_e = sqrt((rhoInf*velInf2 + 2*(pInf-p))/rho);
+    	u_e = sqrt((rhoInf*velInf2 + 2*(p-pInf))/rho);
+
+      }
+
+      /*--- Local pressure-gradient parameter for the boundary layer shape factor. Minimum value of 0.328 for stability ---*/
+      const su2double H_L 		= max(((StrainMag_i*dist_i)/u_e),0.328);
+
+      /*--- Integral shape factor ---*/
+      const su2double H_12 		= 13.9766*pow(H_L,4) - 22.9166*pow(H_L,3) + 13.7227*pow(H_L,2) - 1.0023*H_L + 1.6778;
+
+      /*--- F growth parameters ---*/
+      const su2double DH_12 	= (0.0616*pow(H_12,2) + 0.2339*H_12 + 3.4298)/
+    		  	  	  	  	  	     (0.0047*pow(H_12,3) - 0.1056*pow(H_12,2) + 0.9350*H_12 - 1.2071);
+
+      const su2double lH_12 	= (6.54*H_12 - 14.07)/pow(H_12,2);
+      const su2double mH_12 	= (0.058*(pow((H_12 - 4),2)/(H_12 - 1)) - 0.068)*(1/lH_12);
+
+      const su2double F_growth 	= DH_12*((1 + mH_12)*lH_12)/2;
+
+      /*--- F crit parameters ---*/
+      const su2double Re_y 		= (rho*vel_mag*dist_i)/muLam;
+      const su2double k_y 		= -0.00315*pow(H_12,3) + 0.0986*pow(H_12,2) - 0.242*H_12 + 3.739;
+      const su2double Re_d2_0 	= pow(10,(0.7*tanh((14/(H_12 - 1)) - 9.24) + 2.492/pow((H_12 - 1),0.43) + 0.62));
+      const su2double Re_y_0	= k_y * Re_d2_0;
+
+      short int F_crit = 0;
+      if (Re_y > Re_y_0){
+        F_crit = 1;
+      } else {
+        F_crit = 0;
+      }
+
+      /*--- Source term expresses stream wise growth of Tollmien_schlichting instabilities ---*/
+      const su2double dn_over_dRe_d2 = 0.028*(H_12 - 1) - 0.0345*exp(-pow((3.87/(H_12 - 1) - 2.52),2));
+
+      /*--- Production term ---*/
+      const su2double P_amplification = rho*VorticityMag*F_crit*F_growth*dn_over_dRe_d2;
+
+      /*--- Add Production to residual ---*/
+      Residual += P_amplification * Volume;
+
+      /*--- Implicit part ---*/
+	  Jacobian_i[0] = (rho*VorticityMag*F_crit*F_growth) * Volume;
+
+    }
+
+    AD::SetPreaccOut(Residual);
+    AD::EndPreacc();
+
+    return ResidualType<>(&Residual, &Jacobian_i, nullptr);
+  }
+};