DDES, IDDES and shear-layer-adapted DDES for SST, and fixes to SA EDDES

.github/workflows/regression.yml

@@ -235,7 +235,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:261011-1107
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}}
+          args: -b ${{github.ref}} -t develop -c feature_SST_IDDES -s ${{matrix.testscript}}
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:261011-1107
         with:
@@ -284,7 +284,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:261011-1107
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tapetests"
+          args: -b ${{github.ref}} -t develop -c feature_SST_IDDES -s ${{matrix.testscript}} -a "--tapetests"
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:261011-1107
         with:
@@ -334,7 +334,7 @@ jobs:
           PMIX_MCA_gds: hash
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tsan"
+          args: -b ${{github.ref}} -t develop -c feature_SST_IDDES -s ${{matrix.testscript}} -a "--tsan"
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2-tsan:261011-1107
         with:
@@ -379,7 +379,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2-asan:261011-1107
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--asan"
+          args: -b ${{github.ref}} -t develop -c feature_SST_IDDES -s ${{matrix.testscript}} -a "--asan"
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2-asan:261011-1107
         with:

Common/include/option_structure.hpp

@@ -1196,7 +1196,6 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
   const bool sst_compWilcox = IsPresent(SST_OPTIONS::COMP_Wilcox);
   const bool sst_compSarkar = IsPresent(SST_OPTIONS::COMP_Sarkar);
   const bool sst_dll = IsPresent(SST_OPTIONS::DLL);
-
   if (sst_1994 && sst_2003) {
     SU2_MPI::Error("Two versions (1994 and 2003) selected for SST_OPTIONS. Please choose only one.", CURRENT_FUNCTION);
   } else if (sst_2003) {
@@ -1702,17 +1701,25 @@ static const MapType<std::string, WINDOW_FUNCTION> Window_Map = {
  */
 enum ENUM_HYBRIDRANSLES {
   NO_HYBRIDRANSLES = 0,  /*!< \brief No turbulence model. */
-  SA_DES   = 1,          /*!< \brief Kind of Hybrid RANS/LES (SA - Detached Eddy Simulation (DES)). */
-  SA_DDES  = 2,          /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Delta_max SGS ). */
-  SA_ZDES  = 3,          /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Vorticity based SGS like Zonal DES). */
-  SA_EDDES = 4           /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Shear Layer Adapted SGS: Enhanced DDES). */
+  SA_DES           = 1,  /*!< \brief Kind of Hybrid RANS/LES (SA - Detached Eddy Simulation (DES)). */
+  SA_DDES          = 2,  /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Delta_max SGS ). */
+  SA_ZDES          = 3,  /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Vorticity based SGS like Zonal DES). */
+  SA_EDDES         = 4,  /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Shear Layer Adapted SGS: Enhanced DDES). */
+  SST_DDES         = 5,  /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): DDES). */
+  SST_IDDES        = 6,  /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Improved DDES). */
+  SST_SIDDES       = 7,  /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Simplified Improved DDES). */
+  SST_EDDES        = 8   /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Enhanced (SLA) DDES). */
 };
 static const MapType<std::string, ENUM_HYBRIDRANSLES> HybridRANSLES_Map = {
   MakePair("NONE", NO_HYBRIDRANSLES)
   MakePair("SA_DES", SA_DES)
   MakePair("SA_DDES", SA_DDES)
   MakePair("SA_ZDES", SA_ZDES)
   MakePair("SA_EDDES", SA_EDDES)
+  MakePair("SST_DDES", SST_DDES)
+  MakePair("SST_IDDES", SST_IDDES)
+  MakePair("SST_SIDDES", SST_SIDDES)
+  MakePair("SST_EDDES", SST_EDDES)
 };
 
 /*!

Common/src/CConfig.cpp

@@ -3124,7 +3124,7 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Window (weight) function for the cost-functional in the reverse sweep */
   addEnumOption("WINDOW_FUNCTION", Kind_WindowFct, Window_Map, WINDOW_FUNCTION::SQUARE);
 
-  /* DESCRIPTION: DES Constant */
+  /* DESCRIPTION: DES Constant of the SA hybrid models. The SST hybrid models blend their own constants with F1. */
   addDoubleOption("DES_CONST", Const_DES, 0.65);
 
   /* DESCRIPTION: SBS lengthscale coefficient */
@@ -3145,7 +3145,7 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Specify if the Stochastic Backscatter Model must be activated */
   addBoolOption("STOCHASTIC_BACKSCATTER", SBSParam.StochasticBackscatter, false);
 
-  /* DESCRIPTION: Specify if the LES mode must be enforced */
+  /* DESCRIPTION: Specify if the LES mode must be enforced (SA hybrid models only) */
   addBoolOption("ENFORCE_LES", enforceLES, false);
 
   /* DESCRIPTION: Specify if the stochastic source term must be included in the turbulence model equation */
@@ -3166,7 +3166,7 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Shielding function lower threshold for application of Stochastic Backscatter Model */
   addDoubleOption("SBS_FD_LOWER_THRESHOLD", SBSParam.stochFdThreshold, 0.9);
 
-  /* DESCRIPTION: Filter width for LES (if negative, it is computed based on the local cell size) */
+  /* DESCRIPTION: Filter width for LES (if negative, it is computed based on the local cell size, SA hybrid models only) */
   addDoubleOption("LES_FILTER_WIDTH", LES_FilterWidth, -1.0);
 
   /* DESCRIPTION: Roe with low dissipation for unsteady flows */
@@ -4149,6 +4149,52 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
                    "TIME_STEPPING, DUAL_TIME_STEPPING-1ST_ORDER or DUAL_TIME_STEPPING-2ND_ORDER", CURRENT_FUNCTION);
   }
 
+  if (Kind_HybridRANSLES != NO_HYBRIDRANSLES && !Time_Domain) {
+    SU2_MPI::Error("Hybrid RANS/LES models (HYBRID_RANSLES) require an unsteady simulation, set TIME_DOMAIN= YES.",
+                   CURRENT_FUNCTION);
+  }
+
+  /*--- The SA_* hybrid models are built on the SA model and the SST_* ones on the SST model. ---*/
+  const bool hybridSA = Kind_HybridRANSLES == SA_DES || Kind_HybridRANSLES == SA_DDES ||
+                        Kind_HybridRANSLES == SA_ZDES || Kind_HybridRANSLES == SA_EDDES;
+  const bool hybridSST = Kind_HybridRANSLES == SST_DDES || Kind_HybridRANSLES == SST_IDDES ||
+                         Kind_HybridRANSLES == SST_SIDDES || Kind_HybridRANSLES == SST_EDDES;
+  if ((hybridSA && Kind_Turb_Model != TURB_MODEL::SA) || (hybridSST && Kind_Turb_Model != TURB_MODEL::SST)) {
+    SU2_MPI::Error(string("The hybrid RANS/LES model (HYBRID_RANSLES) must match the turbulence model (KIND_TURB_MODEL): ") +
+                   "SA_* options require SA, SST_* options require SST.", CURRENT_FUNCTION);
+  }
+
+  /*--- The stochastic backscatter model adds the Langevin equations to the SA model (turbulence variables 1 to 3). ---*/
+  if (SBSParam.StochasticBackscatter && hybridSST) {
+    SU2_MPI::Error(string("The stochastic backscatter model (STOCHASTIC_BACKSCATTER= YES) is only available with the SA ") +
+                   "hybrid RANS/LES models (HYBRID_RANSLES= SA_DES, SA_DDES, SA_ZDES or SA_EDDES).", CURRENT_FUNCTION);
+  }
+
+  /*--- How the SST compressibility corrections should act in the LES region of the hybrid models is not validated yet
+   (e.g. Forsythe et al. 2002 keep them only in the RANS length scale). Stop rather than give unvalidated results. ---*/
+  if (hybridSST && (sstParsedOptions.compWilcox || sstParsedOptions.compSarkar)) {
+    SU2_MPI::Error(string("The SST compressibility corrections (COMPRESSIBILITY-WILCOX, COMPRESSIBILITY-SARKAR) are not ") +
+                   "supported yet with the SST hybrid RANS/LES models (HYBRID_RANSLES= SST_*).", CURRENT_FUNCTION);
+  }
+
+  /*--- The SST hybrid models compute their length scale from the grid with the constant C_DES = 0.78 F1 + 0.61 (1 - F1)
+   (Gritskevich et al. 2012), they do not use these options of the SA hybrid models. Stop rather than ignore them. ---*/
+  if (hybridSST) {
+    if (enforceLES) {
+      SU2_MPI::Error("ENFORCE_LES= YES is only available with the SA hybrid RANS/LES models (HYBRID_RANSLES= SA_*).",
+                     CURRENT_FUNCTION);
+    }
+    if (LES_FilterWidth > 0.0) {
+      SU2_MPI::Error(string("A user-specified LES_FILTER_WIDTH is only available with the SA hybrid RANS/LES models ") +
+                     "(HYBRID_RANSLES= SA_*). The SST hybrid models compute the filter width from the grid.", CURRENT_FUNCTION);
+    }
+    if (OptionIsSet("DES_CONST")) {
+      SU2_MPI::Error(string("DES_CONST is not used by the SST hybrid RANS/LES models (HYBRID_RANSLES= SST_*), which blend ") +
+                     "their own constants: C_DES = 0.78 F1 + 0.61 (1 - F1). Remove DES_CONST from the configuration file.",
+                     CURRENT_FUNCTION);
+    }
+  }
+
   if (Time_Domain){
     Delta_UnstTime = Time_Step;
 
@@ -6847,13 +6893,22 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
         cout << "Hybrid RANS/LES: ";
         switch (Kind_HybridRANSLES) {
           case NO_HYBRIDRANSLES: cout << "No Hybrid RANS/LES" << endl; break;
-          case SA_DES:   cout << "Detached Eddy Simulation (DES97) " << endl; break;
-          case SA_DDES:  cout << "Delayed Detached Eddy Simulation (DDES) with Standard SGS" << endl; break;
-          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;
+          case SA_DES:         cout << "Detached Eddy Simulation (DES97) " << endl; break;
+          case SA_DDES:        cout << "Delayed Detached Eddy Simulation (DDES) with Standard SGS" << endl; break;
+          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;
+          case SST_DDES:       cout << "Delayed Detached Eddy Simulation (DDES)" << endl; break;
+          case SST_IDDES:      cout << "Improved Delayed Detached Eddy Simulation (IDDES)" << endl; break;
+          case SST_SIDDES:     cout << "Simplified Improved Delayed Detached Eddy Simulation (SIDDES)" << endl; break;
+          case SST_EDDES:      cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break;
         }
         if (Kind_HybridRANSLES != NO_HYBRIDRANSLES) {
           if (LES_FilterWidth > 0.0) cout << "User-specified LES filter width: " << LES_FilterWidth << endl;
+          if (Kind_HybridRANSLES == SST_DDES || Kind_HybridRANSLES == SST_IDDES ||
+              Kind_HybridRANSLES == SST_SIDDES || Kind_HybridRANSLES == SST_EDDES)
+            cout << "DES constant blended with F1: C_DES = 0.78 F1 + 0.61 (1 - F1)." << endl;
+          else
+            cout << "DES constant: " << Const_DES << endl;
           cout << "Stochastic Backscatter: ";
           if (SBSParam.StochasticBackscatter) {
             cout << "ON" << endl;

SU2_CFD/include/numerics/CNumerics.hpp

@@ -142,6 +142,7 @@ class CNumerics {
   const su2double
   *TurbPsi_i,  /*!< \brief Vector of adjoint turbulent variables at point i. */
   *TurbPsi_j;  /*!< \brief Vector of adjoint turbulent variables at point j. */
+  su2double lengthScale_i, lengthScale_j; /*!< \brief length scale for SST */
   CMatrixView<const su2double>
   ConsVar_Grad_i,  /*!< \brief Gradient of conservative variables at point i. */
   ConsVar_Grad_j,  /*!< \brief Gradient of conservative variables at point j. */
@@ -964,6 +965,16 @@ class CNumerics {
     dist_j = val_dist_j;
   }
 
+  /*!
+   * \brief Set the value of the length scale for SST.
+   * \param[in] val_lengthScale_i - Value of of the length scale for SST from point i.
+   * \param[in] val_lengthScale_j - Value of of the length scale for SST from point j.
+   */
+  void SetLengthScale(su2double val_lengthScale_i, su2double val_lengthScale_j) {
+    lengthScale_i = val_lengthScale_i;
+    lengthScale_j = val_lengthScale_j;
+  }
+
   /*!
    * \brief Set the stochastic source term for the Langevin equations (Backscatter Model).
    * \param[in] val_stoch_source - Value of stochastic source term.

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

@@ -882,6 +882,11 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
     AD::SetPreaccIn(V_i[idx.Velocity() + 1]);
     AD::SetPreaccIn(V_i[idx.SoundSpeed()]);
 
+    /*--- The hybrid RANS/LES length scale depends on the flow and turbulence solution, it must be a preaccumulation
+     input or its sensitivities are lost in the discrete adjoint. ---*/
+    const bool hybridRANSLES = config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES;
+    if (hybridRANSLES) AD::SetPreaccIn(lengthScale_i);
+
     Density_i = V_i[idx.Density()];
     Laminar_Viscosity_i = V_i[idx.LaminarViscosity()];
     Eddy_Viscosity_i = V_i[idx.EddyViscosity()];
@@ -996,6 +1001,9 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       /*--- Dissipation ---*/
 
       su2double dk = beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * (1.0 + zetaFMt);
+      if (hybridRANSLES)
+        dk = Density_i * sqrt(pow(ScalarVar_i[0], 3)) / lengthScale_i;
+
       su2double dw = beta_blended * Density_i * ScalarVar_i[1] * ScalarVar_i[1] * (1.0 - 0.09/beta_blended * zetaFMt);
 
       /*--- LM model coupling with production and dissipation term for k transport equation---*/
@@ -1024,6 +1032,13 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
 
       Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume * (1.0 + zetaFMt);
       Jacobian_i[0][1] = -beta_star * ScalarVar_i[0] * Volume * (1.0 + zetaFMt);
+      if (hybridRANSLES) {
+        /*--- Hybrid k destruction dk = rho k^(3/2) / l_hyb, with the length scale l_hyb frozen (it is computed in the
+         preprocessing and treated as a constant here): d(dk)/d(rho k) = 1.5 dk / (rho k) and it does not depend on omega.
+         The RANS form above would use the wrong derivative where l_hyb differs from l_RANS = sqrt(k) / (beta* omega). ---*/
+        Jacobian_i[0][0] = -1.5 * sqrt(ScalarVar_i[0]) / lengthScale_i * Volume;
+        Jacobian_i[0][1] = 0.0;
+      }
       Jacobian_i[1][0] = 0.0;
       Jacobian_i[1][1] = -2.0 * beta_blended * ScalarVar_i[1] * Volume * (1.0 - 0.09/beta_blended * zetaFMt);
     }

SU2_CFD/include/solvers/CTurbSSTSolver.hpp

@@ -69,6 +69,17 @@ class CTurbSSTSolver final : public CTurbSolver {
   void BoundaryFlux(CGeometry* geometry, CSolver** solver_container, CConfig* config, const ScalarFluxOptions& opt,
                     unsigned short val_marker);
 
+
+   /*!
+   * \brief A virtual member.
+   * \param[in] solver - Solver container
+   * \param[in] geometry - Geometrical definition.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetDES_LengthScale(CSolver** solver,
+                          CGeometry *geometry,
+                          CConfig *config);
+
 public:
   /*!
    * \brief Constructor.

SU2_CFD/include/solvers/CTurbSolver.hpp

@@ -42,6 +42,20 @@ class CTurbSolver : public CScalarSolver<CTurbVariable> {
 
   vector<su2activematrix> Inlet_TurbVars;  /*!< \brief Turbulence variables at inlet profiles */
 
+  /*!
+   * \brief Shear-layer-adapted subgrid length-scale of the EDDES models, Delta_SLA = Delta_omega * F_KH(<VTM>).
+   * \note Shur et al., An Enhanced Version of DES with Rapid Transition from RANS to LES in Separated Flows,
+   *       Flow Turbulence Combust 95, 2015, Eqs. (1)-(6). Delta_omega is computed on the vertices of the median-dual
+   *       control volume of the point, so it applies to any type of grid.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] iPoint - Point index.
+   * \param[in] vorticity - Vorticity vector at the point.
+   * \param[in] shielded - True in the RANS region (f_d < 0.99), where F_KH is set to 1 (Eq. 6).
+   * \return Delta_SLA, or the maximum edge length where the vorticity vanishes.
+   */
+  su2double ShearLayerAdaptedLengthScale(const CGeometry* geometry, unsigned long iPoint, const su2double* vorticity,
+                                         bool shielded) const;
+
 public:
   /*!
    * \brief Destructor of the class.

SU2_CFD/include/variables/CTurbSAVariable.hpp

@@ -39,12 +39,9 @@
 class CTurbSAVariable final : public CTurbVariable {
 
 private:
-  VectorType DES_LengthScale;
-  VectorType lesMode;
   su2vector<int8_t> sbsInBox;
   MatrixType stochSource;
   MatrixType stochSourceOld;
-  VectorType Vortex_Tilting;
   VectorType besselIntegral;
 
 public:
@@ -65,19 +62,6 @@ class CTurbSAVariable final : public CTurbVariable {
    */
   ~CTurbSAVariable() override = default;
 
-  /*!
-   * \brief Get the DES length scale
-   * \param[in] iPoint - Point index.
-   * \return Value of the DES length Scale.
-   */
-  inline su2double GetDES_LengthScale(unsigned long iPoint) const override { return DES_LengthScale(iPoint); }
-
-  /*!
-   * \brief Set the DES Length Scale.
-   * \param[in] iPoint - Point index.
-   */
-  inline void SetDES_LengthScale(unsigned long iPoint, su2double val_des_lengthscale) override { DES_LengthScale(iPoint) = val_des_lengthscale; }
-
   /*!
    * \brief Get the source terms for the stochastic equations.
    * \param[in] iPoint - Point index.
@@ -124,31 +108,6 @@ class CTurbSAVariable final : public CTurbVariable {
    */
   inline void SetSbsInBox(unsigned long iPoint, int8_t val_sbsInBox) override { sbsInBox(iPoint) = val_sbsInBox; }
 
-  /*!
-   * \brief Set the LES sensor.
-   */
-  inline void SetLES_Mode(unsigned long iPoint, su2double val_les_mode) override { lesMode(iPoint) = val_les_mode; }
-
-  /*!
-   * \brief Get the LES sensor.
-   * \return Value of the LES sensor.
-   */
-  inline su2double GetLES_Mode(unsigned long iPoint) const override { return lesMode(iPoint); }
-
-  /*!
-   * \brief Set the vortex tilting measure for computation of the EDDES length scale
-   * \param[in] iPoint - Point index.
-   */
-  void SetVortex_Tilting(unsigned long iPoint, CMatrixView<const su2double>,
-                         const su2double* Vorticity, su2double LaminarViscosity) override;
-
-  /*!
-   * \brief Get the vortex tilting measure for computation of the EDDES length scale
-   * \param[in] iPoint - Point index.
-   * \return Value of the DES length Scale
-   */
-  inline su2double GetVortex_Tilting(unsigned long iPoint) const override { return Vortex_Tilting(iPoint); }
-
   /*!
    * \brief Set the integral of the product of three Bessel functions appearing in Laplacian smoothing.
    * \param[in] iPoint - Point index.