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

AngPass · 2025-09-16T10:34:22Z

Proposed Changes

Hybrid RANS-LES methods aim to leverage the strengths of both Reynolds-Averaged Navier-Stokes (RANS) approaches and Large Eddy Simulation (LES), applying RANS near solid walls for computational efficiency and LES in outer flow regions for improved fidelity. A persistent challenge in these techniques is the appearance of the so-called grey area—the transitional region between RANS and LES zones—which is characterized by the absence of resolved turbulence. This leads to degraded accuracy in capturing free shear layers. The Stochastic Backscatter Model (Kok, 2017) aims to mitigate this issue. The model introduces stochastic forcing terms into the momentum and total energy equations that emulate the backscatter of energy from unresolved (subgrid) scales to resolved scales. The random forcing is formulated as the curl of a stochastic vector potential, say $\underline\xi$. The latter must be constructed using a three-step procedure:

A vector of three mutually-uncorrelated normally-distributed random variables is generated at each node.
A Laplacian smoothing in space (Passariello, 2025) is performed on each component of this vector.
The smoothed fields are used as source terms in three stochastic differential equations (also known as Langevin equations). The latter are solved for the components of the vector potential $\underline\xi$.
The curl of $\underline\xi$ is added to the momentum and total energy equations as source term. This can be done modifying the definition of the modeled stress tensor.

References:

Kok, J. C. (2017). A stochastic backscatter model for grey-area mitigation in detached eddy simulations. Flow, Turbulence and Combustion, 99(1), 119-150.
Passariello, A. (2025). Stochastic Backscatter Model for Unstructured Solvers. arXiv preprint arXiv:2508.17958.

Roadmap

Modifying the definition of the modeled stress tensor using a stochastic contribution. As preliminary step, the vector of uncorrelated variables (point 1 in the procedure above) is employed.
Performing a smoothing on the uncorrelated variables. This requires the solution of a Laplace equation at the beginning of each step.
Implementing three additional equations for the components of the stochastic vector potential. These equations require the use of symmetry-preserving central scheme for the discretization of the convective terms (Kok, 2017).

Related Work

None.

PR Checklist

I am submitting my contribution to the develop branch.
My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
I have added a test case that demonstrates my contribution, if necessary.
I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary.

The Reynolds stress tensor definition is modified to include a random contribution, which is expressed as the curl of a normally-distributed stochastic vector potential.

SU2_CFD/include/numerics/CNumerics.hpp

Add Langevin equations to Spalart-Allmaras solver (uncorrelated random source term) + add stochastic contribution to Reynolds stress tensor (using conservative variables from Langevin equations).

SU2_CFD/include/numerics/CNumerics.hpp

SU2_CFD/src/solvers/CTurbSASolver.cpp

SU2_CFD/src/variables/CTurbSAVariable.cpp

Fix sanitizer warnings and regression test failures - Address issues reported by clang sanitizers. - Adjust implementation to ensure regression tests pass.

- Enhance readability. - Ensure consistent declaration types.

- Extend implementation to 2D configurations. - Improve robustness.

Enforce boundary conditions in Langevin equations.

- Implement a skew-symmetric scheme for the discretization of the convective terms in the Langevin equations. - Fix the seed for the generation of the random source terms at the beginning of the simulation.

- Add the flag for the simulation of the Decaying Isotropic Homogeneous Turbulence. - Generate the initial velocity field matching the experimental data by Comte-Bellot & Corrsin. - Include a preliminary test case.

- Implement flag to enforce LES in the whole domain. - Add option to set the LES filter width to a user-specified value. - Add velocity divergence to volume outputs. - Remove internal generation of initial velocity field for simulating the Decaying Isotropic Homogeneous Turbulence (DIHT).

- Add consistent evaluation of the turbulent kinetic energy in the random source term appearing in the momentum equations. - Add backscatter intensity coefficient in the configuration file. - Add random initialization of the Langevin variables.

- Add Laplacian smoothing of source terms in Langevin equations. Remark: pseudo-time integration is employed (residuals are printed on screen with fixed frequency, the maximum number of time iterations can be defined by the user).

- Fix redundancy in virtual member definition.

- Add stochastic source term to Spalart-Allmaras turbulence model equation (to ensure exchange of modeled and resolved kinetic energy).

- Replace dual-time integration with Successive Over-Relaxation for Laplacian smoothing. - Initialize stochastic vector potential as equal to the stochastic source terms in Langevin equations. - Add random source term to main balance equations in LES zones exclusively. - Blend RANS and LES turbulence timescales using the LES sensor.

- Scale source terms in Langevin equations using Bessel functions to preserve the variance. - Compute Bessel integral at the beginning of the simulation for optimization. - Add variance monitoring to screen output.

- Add LD2 discretization for the inviscid terms in the main balance equations. Remarks: -- Valid only for the incompressible flow solver, with constant density fluid model. -- JST scheme must be selected in the config file. The novel option LD2_SCHEME must be set to YES.

- Retrieve classic second-order central scheme at periodic boundaries.

Add option to include stochastic contribution to turbulence model equation.

- Use logarithmic approximation for Bessel function. - Compute relative variance of the stochastic field for verification.

- Add option to print time-averaged skin friction coefficient in volume output.

- Minor syntax fix in CConfig.cpp

- Include time-averaged skin friction coefficient and instantaneous energy backscatter ratio into volume output fields. - Fix boundary conditions for the implicit smoothing of the stochastic source term in Langevin equations. - Fix computation of the source term in turbulence model equation. - Fix computation of the subgrid kinetic energy (divide eddy viscosity by flow density). - Diagonalize Jacobian in turbulence equations for numerical robustness.

- Allow LD2 option only for JST scheme.

- Add option to include diagnostics of the stochastic source term smoothing. - Add relaxation factor for the random forcing term in the main balanca equations. - Add modeled and stochastic Reynolds stresses to volume output files.

- Random generators removed from static variables. - Random generators re-initialized at every time step using deterministic seeds.

- Account for zero DES length scale in ghost cells. - Fix threshold (=0.9) for the activation of the backscatter model.

SU2_CFD/src/variables/CTurbSAVariable.cpp

+  } else {
+    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+      Solution_Old(iPoint, 0) = Solution(iPoint, 0) = val_nu_tilde;
+      for (unsigned short iVar = 1; iVar < nVar; iVar++) {


SU2_CFD/src/solvers/CTurbSASolver.cpp

+
+    if (backscatter && innerIter==0) {
+      SetLangevinSourceTerms(config, geometry);
+      const unsigned short maxIter = config->GetSBS_maxIterSmooth();


SU2_CFD/src/solvers/CTurbSASolver.cpp

+  const su2double LES_FilterWidth = config->GetLES_FilterWidth();
+  const su2double constDES = config->GetConst_DES();
+  const su2double cDelta = config->GetSBS_Cdelta();
+  const unsigned short maxIter = config->GetSBS_maxIterSmooth();


pcarruscag

Some comments to start, I'll look more into the TurbSolver changes later

pcarruscag · 2026-01-10T22:21:28Z

Common/src/CConfig.cpp

+  /* DESCRIPTION: Specify if the LD2 scheme must be employed (incompressible flows, combined with JST discretization). */
+  addBoolOption("LD2_OPTION", LD2_Scheme, false);


Make this part of CONV_NUM_METHOD enum, ok to still re-use JST in the code, but the option to use the scheme should be consistent with all other schemes.

pcarruscag · 2026-01-10T22:24:03Z

SU2_CFD/include/numerics/CNumerics.hpp

+   * \param[in] velGrad - Velocity gradient matrix.
+   * \param[out] stochReynStress - Stochastic tensor (to be added to the Reynolds stress tensor).
+   */
+  template<class Mat, class Scalar, class Vector>
+  NEVERINLINE static void ComputeStochReynStress(size_t nDim, Scalar density, Scalar tke,
+                                                 Vector rndVec, Mat& stochReynStress, Scalar Cmag) {


Cmag before the output and add to the function doc

pcarruscag · 2026-01-10T22:25:12Z

SU2_CFD/include/numerics/CNumerics.hpp

+  inline void SetStochVar(su2double val_stochvar_i, su2double val_stochvar_j, unsigned short iDim) {
+    stochVar_i[iDim] = val_stochvar_i;
+    stochVar_j[iDim] = val_stochvar_j;
+  }


iDim first in these functions to be more consistent with others

pcarruscag · 2026-01-10T22:26:21Z

SU2_CFD/include/numerics/scalar/scalar_convection.hpp

  const FlowIndices idx;            /*!< \brief Object to manage the access to the flow primitives. */
  su2double a0 = 0.0;               /*!< \brief The maximum of the face-normal velocity and 0. */
  su2double a1 = 0.0;               /*!< \brief The minimum of the face-normal velocity and 0. */
+  su2double qij = 0.0;              /*!< \brief The face-normal velocity (Langevin equations). */


it's not a qdot it's an mdot

pcarruscag · 2026-01-10T22:27:18Z

SU2_CFD/include/numerics/scalar/scalar_convection.hpp

+    if (config->GetStochastic_Backscatter()) {
+      qij = 0.0;
+      for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+        qij += 0.5 * (V_i[idx.Density()]*V_i[iDim + idx.Velocity()] + V_j[idx.Density()]*V_j[iDim + idx.Velocity()]) * Normal[iDim];
+      }
+    }


Can't you use a0 and a1 for this?

pcarruscag · 2026-01-10T22:47:44Z

SU2_CFD/include/variables/CVariable.hpp

+  /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   */
+  inline virtual su2double GetLES_Mode(unsigned long iPoint) const { return 0.0; }
+
+  /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   * \param[in] val_les_mode - Value of the LES sensor.
+   */
+  inline virtual void SetLES_Mode(unsigned long iPoint, su2double val_les_mode) {}
+
+  /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   * \param[in] iDim - Dimension index.
+   */
+  inline virtual su2double GetLangevinSourceTerms(unsigned long iPoint, unsigned short iDim) const { return 0.0; }
+
+  /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   * \param[in] iDim - Dimension index.
+   * \param[in] val_stochSource - Source term in Langevin equations.
+   */
+  inline virtual void SetLangevinSourceTerms(unsigned long iPoint, unsigned short iDim, su2double val_stochSource) {}
+
+  /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   * \param[in] iDim - Dimension index.
+   */
+  inline virtual su2double GetLangevinSourceTermsOld(unsigned long iPoint, unsigned short iDim) const { return 0.0; }
+
+  /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   * \param[in] iDim - Dimension index.
+   * \param[in] val_stochSource_old - Old value of source term in Langevin equations.
+   */
+  inline virtual void SetLangevinSourceTermsOld(unsigned long iPoint, unsigned short iDim, su2double val_stochSource_old) {}


I don't think you need all of these to be virtual, try to remove what you can from CVariable because the turbulence solvers cast to CTurbVariable

pcarruscag · 2026-01-10T22:50:18Z

SU2_CFD/src/numerics/flow/convection/centered.cpp

+    su2double velGrad_i[3][3] = {{0.0}};
+    su2double velGrad_j[3][3] = {{0.0}};
+    for (unsigned short jDim = 0; jDim < nDim; jDim++) {
+      for (iDim = 0; iDim < nDim; iDim++) {
+        velGrad_i[iDim][jDim] = PrimVar_Grad_i[iDim+1][jDim];
+        velGrad_j[iDim][jDim] = PrimVar_Grad_j[iDim+1][jDim];
+      }
+    }
+    for (iDim = 0; iDim < nDim; iDim++) {
+      Velocity_i[iDim] += alpha_LD2 * (velGrad_i[iDim][0] * d_ij[0] +
+                                       velGrad_i[iDim][1] * d_ij[1] +
+                                       velGrad_i[iDim][2] * d_ij[2]);
+      Velocity_j[iDim] -= alpha_LD2 * (velGrad_j[iDim][0] * d_ij[0] +
+                                       velGrad_j[iDim][1] * d_ij[1] +
+                                       velGrad_j[iDim][2] * d_ij[2]);


Is this enough to make it LD2? what about pressure and the dissipation terms?

pcarruscag · 2026-01-10T22:51:20Z

SU2_CFD/src/numerics/flow/flow_diffusion.cpp

+  if (config->GetStochastic_Backscatter()) {
+    for (iVar = 0; iVar < nDim; iVar++)
+      Mean_StochVar[iVar] = 0.5*(stochVar_i[iVar] + stochVar_j[iVar]);
+    su2double SBS_Cmag = config->GetSBS_Cmag();
+    su2double SBS_RelaxFactor = config->GetStochSourceRelax();
+    su2double intensityCoeff = SBS_Cmag;
+    if (SBS_RelaxFactor > 0.0) {
+      su2double FS_Vel = config->GetModVel_FreeStream();
+      su2double ReynoldsLength = config->GetLength_Reynolds();
+      su2double timeScale = ReynoldsLength / FS_Vel;
+      unsigned long timeIter = config->GetTimeIter();
+      unsigned long restartIter = config->GetRestart_Iter();
+      su2double timeStep = config->GetTime_Step();
+      su2double currentTime = (timeIter - restartIter) * timeStep;
+      intensityCoeff = SBS_Cmag * (1.0 - exp(- currentTime / (timeScale*SBS_RelaxFactor)));
+    }
+    ComputeStochReynStress(nDim, Mean_PrimVar[nDim+2], Mean_turb_ke,
+                           Mean_StochVar, stochReynStress, intensityCoeff);


Make a helper function that you can share across numerics

pcarruscag · 2026-01-10T22:52:19Z

SU2_CFD/src/output/CFlowCompOutput.cpp

+  AddVolumeOutput("VELOCITY_DIVERGENCE", "Velocity_Divergence", "DERIVED", "Divergence of the velocity field");
+


after primitives please

pcarruscag · 2026-01-10T22:53:28Z

SU2_CFD/src/output/CFlowCompOutput.cpp

+  const auto VelocityGradient = Node_Flow->GetVelocityGradient(iPoint);
+  su2double divVel = 0.0;
+  for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+    divVel += VelocityGradient[iDim][iDim];
+  }
+  SetVolumeOutputValue("VELOCITY_DIVERGENCE", iPoint, divVel);


Move to common flow outputs

AngPass and others added 2 commits September 15, 2025 18:26

Preliminary implementation: uncorrelated stochastic variables

36b6f7a

The Reynolds stress tensor definition is modified to include a random contribution, which is expressed as the curl of a normally-distributed stochastic vector potential.

Merge branch 'su2code:develop' into develop

d6d6e39

pcarruscag reviewed Sep 18, 2025

View reviewed changes

SU2_CFD/include/numerics/CNumerics.hpp Outdated Show resolved Hide resolved

AngPass and others added 2 commits September 23, 2025 10:12

Add solution of Langevin equations

cfa0cd8

Add Langevin equations to Spalart-Allmaras solver (uncorrelated random source term) + add stochastic contribution to Reynolds stress tensor (using conservative variables from Langevin equations).

Merge branch 'develop' into develop

fd800d3

pcarruscag changed the title ~~Stochastic Backscatter Model for Grey Area Mitigation in Hybrid RANS-LES Simulations~~ [WIP] Stochastic Backscatter Model for Grey Area Mitigation in Hybrid RANS-LES Simulations Sep 23, 2025

bigfooted added the changelog:feature label Sep 23, 2025

github-advanced-security bot found potential problems Sep 23, 2025

View reviewed changes

SU2_CFD/include/numerics/CNumerics.hpp Fixed Show fixed Hide fixed

SU2_CFD/include/numerics/CNumerics.hpp Fixed Show fixed Hide fixed

SU2_CFD/src/solvers/CTurbSASolver.cpp Fixed Show fixed Hide fixed

SU2_CFD/src/variables/CTurbSAVariable.cpp Fixed Show fixed Hide fixed

AngPass and others added 22 commits September 24, 2025 11:05

Fix sanitizer warnings and regression failures

34f676f

Fix sanitizer warnings and regression test failures - Address issues reported by clang sanitizers. - Adjust implementation to ensure regression tests pass.

Minor fixes

66293fd

- Enhance readability. - Ensure consistent declaration types.

Merge branch 'develop' into develop

4a11c8b

Ensure consistency with DES

3481df3

- Extend implementation to 2D configurations. - Improve robustness.

Merge branch 'develop' into develop

26bcfda

Boundary conditions for Langevin equations

85b8aac

Enforce boundary conditions in Langevin equations.

Use symmetry-preserving scheme for Langevin equations

6293042

- Implement a skew-symmetric scheme for the discretization of the convective terms in the Langevin equations. - Fix the seed for the generation of the random source terms at the beginning of the simulation.

Merge branch 'develop' into develop

b03f630

Add option for simulating DIHT

1e8e4e1

- Add the flag for the simulation of the Decaying Isotropic Homogeneous Turbulence. - Generate the initial velocity field matching the experimental data by Comte-Bellot & Corrsin. - Include a preliminary test case.

Merge branch 'develop' into develop

25a32f7

Merge branch 'develop' into develop

e1181d3

Add Laplacian smoothing (Langevin equations)

968d8f9

- Add Laplacian smoothing of source terms in Langevin equations. Remark: pseudo-time integration is employed (residuals are printed on screen with fixed frequency, the maximum number of time iterations can be defined by the user).

Fix redundancy in CTurbSAVariable.hpp

2f8245e

- Fix redundancy in virtual member definition.

Add stochastic source to turbulence model equation

bdbfa05

- Add stochastic source term to Spalart-Allmaras turbulence model equation (to ensure exchange of modeled and resolved kinetic energy).

Merge branch 'develop' into develop

7457d95

Merge branch 'su2code:develop' into develop

ddb99bf

Correct scaling of stochastic source terms in Langevin eqs.

15dbff3

- Scale source terms in Langevin equations using Bessel functions to preserve the variance. - Compute Bessel integral at the beginning of the simulation for optimization. - Add variance monitoring to screen output.

Merge branch 'su2code:develop' into develop

4dd0f90

Merge branch 'su2code:develop' into develop

d93ed49

AngPass and others added 22 commits November 26, 2025 14:56

Merge branch 'develop' of https://github.com/AngPass/SU2 into develop

02dbb54

Merge branch 'develop' into develop

2f654c6

Merge branch 'develop' into develop

65e531e

Fix LD2 for periodic boundaries

5c3e35c

- Retrieve classic second-order central scheme at periodic boundaries.

Merge branch 'develop' into develop

71dca9c

Include stochastic source term in turb. equation

6389023

Add option to include stochastic contribution to turbulence model equation.

Minor fixes

d382a2d

- Use logarithmic approximation for Bessel function. - Compute relative variance of the stochastic field for verification.

Add time-averaged skin friction coefficient

d1a1ae8

- Add option to print time-averaged skin friction coefficient in volume output.

Minor fix

babfb58

- Minor syntax fix in CConfig.cpp

Minor fix

562f311

- Minor syntax fix in CConfig.cpp

Merge branch 'develop' into develop

14323bb

Merge branch 'develop' into develop

92a016b

Fix LD2 option in config file

a5b3aa0

- Allow LD2 option only for JST scheme.

Merge branch 'develop' into develop

ebb6806

Merge branch 'develop' into develop

3e44cd6

Merge branch 'develop' into develop

3530ea8

Enhance numerical robustness

0df34ac

- Add option to include diagnostics of the stochastic source term smoothing. - Add relaxation factor for the random forcing term in the main balanca equations. - Add modeled and stochastic Reynolds stresses to volume output files.

Merge branch 'develop' into develop

343ed9b

Merge branch 'develop' into develop

f539511

Enhance numerical efficiency

443838f

- Random generators removed from static variables. - Random generators re-initialized at every time step using deterministic seeds.

AngPass requested a review from pcarruscag January 8, 2026 17:06

AngPass and others added 2 commits January 8, 2026 18:07

Merge branch 'develop' into develop

86ba8d8

Fix bug in viscous flux computation

fb3f8ec

- Account for zero DES length scale in ghost cells. - Fix threshold (=0.9) for the activation of the backscatter model.

github-advanced-security bot found potential problems Jan 10, 2026

View reviewed changes

pcarruscag reviewed Jan 10, 2026

View reviewed changes

		/* DESCRIPTION: Specify if the LD2 scheme must be employed (incompressible flows, combined with JST discretization). */
		addBoolOption("LD2_OPTION", LD2_Scheme, false);

		AddVolumeOutput("VELOCITY_DIVERGENCE", "Velocity_Divergence", "DERIVED", "Divergence of the velocity field");

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

Are you sure you want to change the base?

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

Conversation

AngPass commented Sep 16, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Proposed Changes

Roadmap

Related Work

PR Checklist

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Check failure

Uh oh!

Uh oh!

Check warning

Check warning

pcarruscag left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

AngPass commented Sep 16, 2025 •

edited

Loading