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

Common/include/CConfig.hpp

@@ -1081,6 +1081,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 */
@@ -9524,6 +9525,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/include/toolboxes/random_toolbox.hpp

@@ -0,0 +1,85 @@
+/*!
+ * \file random_toolbox.hpp
+ * \brief Collection of utility functions for random number generation.
+ * \version 8.3.0 "Harrier"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md)
+ *
+ * SU2 is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * SU2 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include <random>
+#include <functional>
+
+namespace RandomToolbox {
+/// \addtogroup RandomToolbox
+/// @{
+
+/*!
+ * \brief Combine two 64-bit integers into a single hash value.
+ * \param[in] v1 Current hash value.
+ * \param[in] v2 Value to mix in.
+ * \return Combined hash value.
+ */
+inline uint64_t HashCombine(uint64_t v1, uint64_t v2) {
+  const uint64_t prime = 1099511628211ULL;
+  v1 ^= v2;
+  v1 *= prime;
+  return v1;
+}
+
+/*!
+ * \brief Convert a double to a 64-bit integer suitable for hashing.
+ * \param[in] x Double to integer.
+ * \return Hash value of the double.
+ */
+inline uint64_t ToUInt64(double x) {
+  return std::hash<double>{}(x); 
+}
+
+/*!
+ * \brief Build a deterministic seed from physical time.
+ * \param[in] time Physical time.
+ * \return 64-bit seed value.
+ */
+inline uint64_t GetSeed(double x) {
+  uint64_t h = 1469598103934665603ULL; // Offset
+  h = HashCombine(h, ToUInt64(x));
+  return h;
+}
+
+/*!
+ * \brief Generate a standard normally-distributed random number.
+ * \param[in] seed Seed for the random number generator.
+ * \param[in] mean Mean of the normal distribution (default 0).
+ * \param[in] stddev Standard deviation of the normal distribution (default 1).
+ * \return Normally-distributed random number.
+ */
+inline double GetRandomNormal(uint64_t seed,
+                                      double mean = 0.0,
+                                      double stddev = 1.0) {
+  std::mt19937 gen(static_cast<unsigned int>(seed)); 
+  std::normal_distribution<double> rnd(mean, stddev);
+  return rnd(gen);
+}
+
+/// @}
+}  // namespace RandomToolbox

Common/src/CConfig.cpp

@@ -2914,6 +2914,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);
 
@@ -6447,6 +6450,13 @@ 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;
+        if (StochasticBackscatter && Kind_HybridRANSLES == NO_HYBRIDRANSLES)
+          SU2_MPI::Error("Stochastic Backscatter can only be activated with Hybrid RANS/LES.", CURRENT_FUNCTION);
         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"
@@ -179,6 +180,15 @@
   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. */
+  su2double
+  lesSensor_i,        /*!< \brief LES sensor at point i. */
+  lesSensor_j;        /*!< \brief LES sensor at point j. */
+  su2double lastTime; /*!< \brief Physical time of unsteady simulation. */
+  su2double stochSource[3]; /*!< \brief Source term for Langevin equations in Stochastic Backscatter Model. */
+  const su2double
+  *stochVar_i, /*!< \brief Stochastic variables at point i for Stochastic Backscatter Model. */
+  *stochVar_j; /*!< \brief Stochastic variables at point j for Stochastic 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 */
@@ -633,6 +643,50 @@
       }
     }
   }
+
+  /*!
+   * \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, const su2double *rndVec, 
+                                                 Mat2& stochReynStress) {
+
+    /* --- 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.
@@ -817,6 +871,16 @@
     turb_ke_j = val_turb_ke_j;
   }
 
+  /*!
+   * \brief Set the stochastic variables from Langevin equations (Stochastic Backscatter Model).
+   * \param[in] val_stochvar_i - Value of the stochastic variable at point i.
+   * \param[in] val_stochvar_j - Value of the stochastic variable at point j.
+   */
+  inline void SetStochVar(su2double *val_stochvar_i, su2double *val_stochvar_j) {
+    stochVar_i = val_stochvar_i;
+    stochVar_j = val_stochvar_j;
+  }
+
   /*!
    * \brief Set the value of the distance from the nearest wall.
    * \param[in] val_dist_i - Value of of the distance from point i to the nearest wall.
@@ -827,6 +891,39 @@
     dist_j = val_dist_j;
   }
 
+  /*!
+   * \brief Set the value of the LES sensor.
+   * \param[in] val_les_i - Value of the LES sensor at point point i.
+   * \param[in] val_les_j - Value of the LES sensor at point point j.
+   */
+  void SetLESSensor(su2double val_les_i, su2double val_les_j) {
+    lesSensor_i = val_les_i;
+    lesSensor_j = val_les_j;
+  }
+
+  /*!
+   * \brief Set the value of physical time.
+   * \param[in] val_last_time - Value of physical time.
+   */
+  void SetLastTime(su2double val_last_time) {
+    lastTime = val_last_time;
+  }
+
+  /*!
+   * \brief Get the value of physical time.
+   * \param[out] lastTime - Value of physical time.
+   */
+  inline su2double GetLastTime() const { return lastTime; }
+
+  /*!
+   * \brief Set the stochastic source term for the Langevin equations (Backscatter Model).
+   * \param[in] val_stoch_source - Value of stochastic source term.
+   * \param[in] iDim - Index of Langevin equation.
+   */
+  void SetStochSource(su2double val_stoch_source, unsigned short iDim) {
+    stochSource[iDim] = val_stoch_source;
+  }
+
   /*!
    * \brief Set the value of the roughness from the nearest wall.
    * \param[in] val_dist_i - Value of of the roughness of the nearest wall from point i

SU2_CFD/include/numerics/flow/flow_diffusion.hpp

@@ -61,6 +61,7 @@ class CAvgGrad_Base : public CNumerics {
   Mean_Cp,                                /*!< \brief Mean value of the specific heat capacity at constant pressure. */
   Mean_turb_ke,                           /*!< \brief Mean value of the turbulent kinetic energy. */
   Mean_TauWall,                           /*!< \brief Mean wall shear stress (wall functions). */
+  *Mean_StochVar,                         /*!< \brief Mean stochastic variables (Stochastic Backscatter Model). */
   TauWall_i, TauWall_j,                   /*!< \brief Wall shear stress at point i and j (wall functions). */
   dist_ij_2,                              /*!< \brief Length of the edge and face, squared */
   Edge_Vector[MAXNDIM] = {0.0},           /*!< \brief Vector from point i to point j. */
@@ -200,12 +201,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/include/numerics/turbulent/turb_convection.hpp

@@ -59,9 +59,20 @@ class CUpwSca_TurbSA final : public CUpwScalar<FlowIndices> {
    * \param[in] config - Definition of the particular problem.
    */
   void FinishResidualCalc(const CConfig* config) override {
-    Flux[0] = a0*ScalarVar_i[0] + a1*ScalarVar_j[0];
-    Jacobian_i[0][0] = a0;
-    Jacobian_j[0][0] = a1;
+    bool backscatter = config->GetStochastic_Backscatter();
+    if (!backscatter) {
+      Flux[0] = a0*ScalarVar_i[0] + a1*ScalarVar_j[0];
+      Jacobian_i[0][0] = a0;
+      Jacobian_j[0][0] = a1;
+    } else {
+      for (unsigned short iVar = 0; iVar < 4; iVar++) {
+        Flux[iVar] = a0*ScalarVar_i[iVar] + a1*ScalarVar_j[iVar];
+        for (unsigned short jVar = 0; jVar < 4; jVar++) {
+          Jacobian_i[iVar][jVar] = (iVar == jVar) ? a0 : 0.0;
+          Jacobian_j[iVar][jVar] = (iVar == jVar) ? a1 : 0.0;
+        }
+      }
+    }
   }
 
 public: