[WIP] Stochastic Backscatter Model for Grey Area Mitigation in Hybrid RANS-LES Simulations

Common/include/CConfig.hpp

@@ -1077,6 +1077,7 @@ class CConfig {
   su2double Const_DES;                 /*!< \brief Detached Eddy Simulation Constant. */
   WINDOW_FUNCTION Kind_WindowFct;      /*!< \brief Type of window (weight) function for objective functional. */
   unsigned short Kind_HybridRANSLES;   /*!< \brief Kind of Hybrid RANS/LES. */
+  bool StochasticBackscatter;          /*!< \brief Option to include Stochastic Backscatter Model. */
   unsigned short Kind_RoeLowDiss;      /*!< \brief Kind of Roe scheme with low dissipation for unsteady flows. */
 
   unsigned short nSpanWiseSections; /*!< \brief number of span-wise sections */
@@ -9466,6 +9467,12 @@ class CConfig {
    */
   unsigned short GetKind_HybridRANSLES(void) const { return Kind_HybridRANSLES; }
 
+  /*!
+   * \brief Get if the Stochastic Backscatter Model must be activated.
+   * \return TRUE if the Stochastic Backscatter Model is activated.
+   */
+  bool GetStochastic_Backscatter(void) const { return StochasticBackscatter; }
+
   /*!
    * \brief Get the Kind of Roe Low Dissipation Scheme for Unsteady flows.
    * \return Value of Low dissipation approach.

Common/src/CConfig.cpp

@@ -2911,6 +2911,9 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Specify Hybrid RANS/LES model */
   addEnumOption("HYBRID_RANSLES", Kind_HybridRANSLES, HybridRANSLES_Map, NO_HYBRIDRANSLES);
 
+  /* DESCRIPTION: Specify if the Stochastic Backscatter Model must be activated */
+  addBoolOption("STOCHASTIC_BACKSCATTER", StochasticBackscatter, false);
+
   /* DESCRIPTION: Roe with low dissipation for unsteady flows */
   addEnumOption("ROE_LOW_DISSIPATION", Kind_RoeLowDiss, RoeLowDiss_Map, NO_ROELOWDISS);
 
@@ -6448,6 +6451,11 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
           case SA_ZDES:  cout << "Delayed Detached Eddy Simulation (DDES) with Vorticity-based SGS" << endl; break;
           case SA_EDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break;
         }
+        cout << "Stochastic Backscatter: ";
+        if (StochasticBackscatter)
+          cout << "ON" << endl;
+        else
+          cout << "OFF" << endl;
         break;
       case MAIN_SOLVER::NEMO_EULER:
         if (Kind_Regime == ENUM_REGIME::COMPRESSIBLE) cout << "Compressible two-temperature thermochemical non-equilibrium Euler equations." << endl;

SU2_CFD/include/numerics/CNumerics.hpp

@@ -32,6 +32,7 @@
 #include <iostream>
 #include <limits>
 #include <cstdlib>
+#include <random>
 
 #include "../../../Common/include/CConfig.hpp"
 #include "../../../Common/include/linear_algebra/blas_structure.hpp"
@@ -180,6 +181,7 @@ class CNumerics {
   roughness_j = 0.0;                       /*!< \brief Roughness of the wall nearest to point j. */
 
   su2double MeanPerturbedRSM[3][3];   /*!< \brief Perturbed Reynolds stress tensor  */
+  su2double stochReynStress[3][3];    /*!< \brief Stochastic contribution to Reynolds stress tensor for Backscatter Model. */
   SST_ParsedOptions sstParsedOptions; /*!< \brief additional options for the SST turbulence model */
   unsigned short Eig_Val_Comp;    /*!< \brief Component towards which perturbation is perfromed */
   su2double uq_delta_b;           /*!< \brief Magnitude of perturbation */
@@ -635,6 +637,62 @@ class CNumerics {
     }
   }
 
+  /*!
+   * \brief Compute a random contribution to the Reynolds stress tensor (Stochastic Backscatter Model).
+   * \details See: Kok, Johan C. "A stochastic backscatter model for grey-area mitigation in detached
+   * eddy simulations." Flow, Turbulence and Combustion 99.1 (2017): 119-150.
+   * \param[in] nDim - Dimension of the flow problem, 2 or 3.
+   * \param[in] density - Density.
+   * \param[in] eddyVis - Eddy viscosity.
+   * \param[in] velGrad - Velocity gradient matrix.
+   * \param[out] stochReynStress - Stochastic tensor (to be added to the Reynolds stress tensor).
+   */
+  template<class Mat1, class Mat2, class Scalar>
+  NEVERINLINE static void ComputeStochReynStress(size_t nDim, Scalar density, Scalar eddyVis,
+                                                 const Mat1& velGrad, Mat2& stochReynStress) {
+
+    /* --- Initialize seed ---*/
+
+    static std::default_random_engine gen;
+    std::random_device rd;
+    gen.seed(rd());
+
+    /* --- Generate a vector of three independent normally-distributed samples ---*/                                              
+
+    std::normal_distribution<Scalar> rnd(0.0,1.0);
+    Scalar rndVec [3] = {0.0};
+    for (size_t iDim = 0; iDim < nDim; iDim++) 
+      rndVec[iDim] = rnd(gen);
+
+    /* --- Estimate turbulent kinetic energy --- */
+
+    Scalar turbKE = 0.0, strainMag = 0.0;
+    for (size_t iDim = 0; iDim < nDim; iDim++) {
+      for (size_t jDim = 0; jDim < nDim; jDim++) {
+        strainMag += pow(0.5 * (velGrad[iDim][jDim] + velGrad[jDim][iDim]), 2);
+      }
+    }
+    strainMag = sqrt(2.0 * strainMag);
+    turbKE = eddyVis * strainMag;
+    turbKE = max(turbKE, 1E-10);
+
+    /* --- Calculate stochastic tensor --- */
+
+    stochReynStress[1][0] = - density * turbKE * rndVec[2];
+    stochReynStress[2][0] =   density * turbKE * rndVec[1];
+    stochReynStress[2][1] = - density * turbKE * rndVec[0];
+    for (size_t iDim = 0; iDim < nDim; iDim++) {
+      for (size_t jDim = 0; jDim <= iDim; jDim++) {
+        if (iDim==jDim) {
+          stochReynStress[iDim][jDim] = 0.0;
+        } else {
+          stochReynStress[jDim][iDim] = - stochReynStress[iDim][jDim];
+        }
+      }
+    }
+
+  }
+
   /*!
    * \brief Project average gradient onto normal (with or w/o correction) for viscous fluxes of scalar quantities.
    * \param[in] nDim - Dimension of the flow problem, 2 or 3.

SU2_CFD/include/numerics/flow/flow_diffusion.hpp

@@ -198,12 +198,14 @@ class CAvgGrad_Base : public CNumerics {
    * \param[in] val_turb_ke - Turbulent kinetic energy
    * \param[in] val_laminar_viscosity - Laminar viscosity.
    * \param[in] val_eddy_viscosity - Eddy viscosity.
+   * \param[in] config - Definition of the particular problem.
    */
   void SetStressTensor(const su2double *val_primvar,
                        const su2double* const *val_gradprimvar,
                        su2double val_turb_ke,
                        su2double val_laminar_viscosity,
-                       su2double val_eddy_viscosity);
+                       su2double val_eddy_viscosity,
+                       const CConfig* config);
 
   /*!
    * \brief Get a component of the viscous stress tensor.

SU2_CFD/src/numerics/flow/flow_diffusion.cpp

@@ -122,7 +122,8 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar,
                            const su2double* const *val_gradprimvar,
                            const su2double val_turb_ke,
                            const su2double val_laminar_viscosity,
-                           const su2double val_eddy_viscosity) {
+                           const su2double val_eddy_viscosity,
+                           const CConfig* config) {
 
   const su2double Density = val_primvar[nDim+2];
 
@@ -142,6 +143,15 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar,
     // turb_ke is not considered in the stress tensor, see #797
     ComputeStressTensor(nDim, tau, val_gradprimvar+1, total_viscosity, Density, su2double(0.0));
   }
+
+  /* --- If the Stochastic Backscatter Model is active, add random contribution to stress tensor ---*/
+
+  if (config->GetStochastic_Backscatter()) {
+    for (unsigned short iDim = 0 ; iDim < nDim; iDim++)
+      for (unsigned short jDim = 0 ; jDim < nDim; jDim++)
+        tau[iDim][jDim] += stochReynStress[iDim][jDim];
+  }
+
 }
 
 void CAvgGrad_Base::AddTauWall(const su2double *UnitNormal,
@@ -432,10 +442,17 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config)
                         Mean_turb_ke, MeanPerturbedRSM);
   }
 
+  /* --- If the Stochastic Backscatter Model is active, add random contribution to stress tensor ---*/
+
+  if (config->GetStochastic_Backscatter()) {
+    ComputeStochReynStress(nDim, Mean_PrimVar[nDim+2], Mean_Eddy_Viscosity, Mean_GradPrimVar+1,
+                           stochReynStress);
+  }
+
   /*--- Get projected flux tensor (viscous residual) ---*/
 
   SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke,
-                  Mean_Laminar_Viscosity, Mean_Eddy_Viscosity);
+                  Mean_Laminar_Viscosity, Mean_Eddy_Viscosity,config);
   if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau);
   if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall);
 
@@ -617,9 +634,16 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi
                         Mean_turb_ke, MeanPerturbedRSM);
   }
 
+  /* --- If the Stochastic Backscatter Model is active, add random contribution to stress tensor ---*/
+
+  if (config->GetStochastic_Backscatter()) {
+    ComputeStochReynStress(nDim, Mean_PrimVar[nDim+2], Mean_Eddy_Viscosity, Mean_GradPrimVar+1,
+                           stochReynStress);
+  }
+
   /*--- Get projected flux tensor (viscous residual) ---*/
   SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke,
-                  Mean_Laminar_Viscosity, Mean_Eddy_Viscosity);
+                  Mean_Laminar_Viscosity, Mean_Eddy_Viscosity,config);
   if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau);
   if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall);
 
@@ -947,10 +971,17 @@ CNumerics::ResidualType<> CGeneralAvgGrad_Flow::ComputeResidual(const CConfig* c
                         Mean_turb_ke, MeanPerturbedRSM);
   }
 
+  /* --- If the Stochastic Backscatter Model is active, add random contribution to stress tensor ---*/
+
+  if (config->GetStochastic_Backscatter()) {
+    ComputeStochReynStress(nDim, Mean_PrimVar[nDim+2], Mean_Eddy_Viscosity, Mean_GradPrimVar+1,
+                           stochReynStress);
+  }
+
   /*--- Get projected flux tensor (viscous residual) ---*/
 
   SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke,
-                  Mean_Laminar_Viscosity, Mean_Eddy_Viscosity);
+                  Mean_Laminar_Viscosity, Mean_Eddy_Viscosity,config);
   if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau);
   if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall);
 

TestCases/ddes/flatplate/ddes_flatplate.cfg

@@ -13,9 +13,11 @@
 %
 SOLVER= RANS
 KIND_TURB_MODEL= SA
-HYBRID_RANSLES= SA_EDDES
+HYBRID_RANSLES= SA_DES
+STOCHASTIC_BACKSCATTER= YES
 MATH_PROBLEM= DIRECT
 RESTART_SOL= NO
+RESTART_ITER= 100
 
 % ----------- COMPRESSIBLE AND INCOMPRESSIBLE FREE-STREAM DEFINITION ----------%
 %
@@ -41,7 +43,7 @@ TIME_MARCHING= DUAL_TIME_STEPPING-2ND_ORDER
 %
 % U_inf = 69.4448 - dt*=0.02 - dt=0.000288
 TIME_STEP= 0.000288
-MAX_TIME= 20.0
+MAX_TIME= 5
 UNST_CFL_NUMBER= 0.0
 INNER_ITER= 20
 
@@ -61,7 +63,7 @@ CFL_NUMBER= 10.0
 CFL_ADAPT= NO
 CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 )
 RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 )
-TIME_ITER= 10
+TIME_ITER= 500
 
 % ----------------------- SLOPE LIMITER DEFINITION ----------------------------%
 %
@@ -103,5 +105,5 @@ VOLUME_ADJ_FILENAME= adjoint
 GRAD_OBJFUNC_FILENAME= of_grad
 SURFACE_FILENAME= surface_flow
 SURFACE_ADJ_FILENAME= surface_adjoint
-OUTPUT_WRT_FREQ= 1000
+OUTPUT_WRT_FREQ= 1000000
 SCREEN_OUTPUT= (TIME_ITER, INNER_ITER, RMS_DENSITY, RMS_NU_TILDE, LIFT, DRAG, TOTAL_HEATFLUX)

config_template.cfg

@@ -180,6 +180,9 @@ HYBRID_RANSLES= SA_DDES
 %
 % DES Constant (0.65)
 DES_CONST= 0.65
+%
+% Stochastic Backscatter Model (NO, YES)
+STOCHASTIC_BACKSCATTER= NO
 
 % -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------%
 %