[WIP] Implementation of DDES and SAS formulations for SST

Common/include/option_structure.hpp

@@ -1380,14 +1380,20 @@ enum ENUM_HYBRIDRANSLES {
   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_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). */
 };
 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)
 };
 
 /*!

Common/src/CConfig.cpp

@@ -6157,6 +6157,9 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
           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) with Shear-layer Adapted SGS" << endl; break;
+          case SST_IDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break;
+          case SST_SIDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break;
         }
         break;
       case MAIN_SOLVER::NEMO_EULER:

SU2_CFD/include/numerics/CNumerics.hpp

@@ -190,6 +190,8 @@ class CNumerics {
 
   bool bounded_scalar = false;    /*!< \brief Flag for bounded scalar problem */
 
+  su2double lengthScale_i, lengthScale_j;
+
 public:
   /*!
    * \brief Return type used in some "ComputeResidual" overloads to give a
@@ -828,6 +830,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 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/turbulent/turb_sources.hpp

@@ -834,6 +834,9 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       /*--- Dissipation ---*/
 
       su2double dk = beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0];
+      if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES)
+        dk = Density_i * sqrt(ScalarVar_i[0]*ScalarVar_i[0]*ScalarVar_i[0]) / lengthScale_i;
+
       su2double dw = beta_blended * Density_i * ScalarVar_i[1] * ScalarVar_i[1];
 
       /*--- LM model coupling with production and dissipation term for k transport equation---*/

SU2_CFD/include/solvers/CTurbSSTSolver.hpp

@@ -54,6 +54,16 @@ class CTurbSSTSolver final : public CTurbSolver {
                      const CConfig *config,
                      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/variables/CTurbSAVariable.hpp

@@ -39,7 +39,6 @@
 class CTurbSAVariable final : public CTurbVariable {
 
 private:
-  VectorType DES_LengthScale;
   VectorType Vortex_Tilting;
 
 public:
@@ -60,19 +59,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 Set the vortex tilting measure for computation of the EDDES length scale
    * \param[in] iPoint - Point index.

SU2_CFD/include/variables/CTurbSSTVariable.hpp

@@ -45,6 +45,13 @@ class CTurbSSTVariable final : public CTurbVariable {
   VectorType F2;    /*!< \brief Menter blending function for blending of k-w and k-eps. */
   VectorType CDkw;  /*!< \brief Cross-diffusion. */
   SST_ParsedOptions sstParsedOptions;
+
+  VectorType ftilda_d;
+  VectorType l_RANS;
+  VectorType l_LES;
+  VectorType r_dl;
+  VectorType r_dt;
+  VectorType r_d;
 public:
   /*!
    * \brief Constructor of the class.
@@ -87,4 +94,37 @@ class CTurbSSTVariable final : public CTurbVariable {
    * \brief Get the value of the cross diffusion of tke and omega.
    */
   inline su2double GetCrossDiff(unsigned long iPoint) const override { return CDkw(iPoint); }
+
+  /*!
+   * \brief Set the DES Length Scale.
+   * \param[in] iPoint - Point index.
+   */
+  inline void SetDebug_Quantities(CConfig *config, unsigned long iPoint, su2double val_ftilda_d, su2double val_l_RANS, su2double val_l_LES, su2double val_r_d) override { 
+    ftilda_d(iPoint) = val_ftilda_d;
+    l_RANS(iPoint) = val_l_RANS;
+    l_LES(iPoint) = val_l_LES;
+
+    if ( config->GetKind_HybridRANSLES() == SST_DDES)
+      r_d(iPoint) = val_r_d;
+    else
+      r_dt(iPoint) = val_r_d;
+  }
+
+  inline void SetDebug_Quantities(CConfig *config, unsigned long iPoint, su2double val_ftilda_d, su2double val_l_RANS, su2double val_l_LES, su2double val_r_dl, su2double val_r_dt) override { 
+    ftilda_d(iPoint) = val_ftilda_d;
+    l_RANS(iPoint) = val_l_RANS;
+    l_LES(iPoint) = val_l_LES;
+    r_dt(iPoint) = val_r_dt;
+    r_dl(iPoint) = val_r_dl;
+  }
+
+  inline su2double Get_L_RANS(unsigned long iPoint) const override { return l_RANS(iPoint); }
+  inline su2double Get_L_LES(unsigned long iPoint) const override { return l_LES(iPoint); }
+  inline su2double Get_ftilda_d(unsigned long iPoint) const override { return ftilda_d(iPoint); }
+  inline su2double Get_r_dt(unsigned long iPoint) const override { return r_dt(iPoint); }
+  inline su2double Get_r_dl(unsigned long iPoint) const override { return r_dl(iPoint); }
+  inline su2double Get_r_d(unsigned long iPoint) const override { return r_d(iPoint); }
+
+
+
 };

SU2_CFD/include/variables/CTurbVariable.hpp

@@ -38,6 +38,7 @@
 class CTurbVariable : public CScalarVariable {
 protected:
   VectorType muT; /*!< \brief Eddy viscosity. */
+  VectorType DES_LengthScale;
 
 public:
   static constexpr size_t MAXNVAR = 2;
@@ -100,5 +101,18 @@ class CTurbVariable : public CScalarVariable {
    */
   inline void SetIntermittency(unsigned long iPoint, su2double val_intermittency) final { intermittency(iPoint) = val_intermittency; }
 
+  /*!
+   * \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; }
+
 };
 

SU2_CFD/include/variables/CVariable.hpp

@@ -381,6 +381,25 @@ class CVariable {
    */
   inline virtual void SetDES_LengthScale(unsigned long iPoint, su2double val_des_lengthscale) {}
 
+  /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   */
+  inline virtual void SetDebug_Quantities(CConfig *config, unsigned long iPoint, su2double val_ftilda_d, su2double val_l_RANS, su2double val_l_LES, su2double val_r_d) {}
+
+  inline virtual su2double Get_L_RANS(unsigned long iPoint) const { return 0.0; }
+  inline virtual su2double Get_L_LES(unsigned long iPoint) const { return 0.0; }
+  inline virtual su2double Get_ftilda_d(unsigned long iPoint) const { return 0.0; }
+  inline virtual su2double Get_r_dt(unsigned long iPoint) const { return 0.0; }
+  inline virtual su2double Get_r_dl(unsigned long iPoint) const { return 0.0; }
+  inline virtual su2double Get_r_d(unsigned long iPoint) const { return 0.0; }
+
+  /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   */
+  inline virtual void SetDebug_Quantities(CConfig *config, unsigned long iPoint, su2double val_ftilda_d, su2double val_l_RANS, su2double val_l_LES, su2double val_r_dl, su2double val_r_dt) {}
+
   /*!
    * \brief A virtual member.
    * \param[in] iPoint - Point index.

SU2_CFD/src/output/CFlowOutput.cpp

@@ -1464,6 +1464,24 @@ void CFlowOutput::SetVolumeOutputFieldsScalarMisc(const CConfig* config) {
   if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES) {
     AddVolumeOutput("DES_LENGTHSCALE", "DES_LengthScale", "DDES", "DES length scale value");
     AddVolumeOutput("WALL_DISTANCE", "Wall_Distance", "DDES", "Wall distance value");
+    if ( config->GetKind_HybridRANSLES() == SST_DDES){
+      AddVolumeOutput("F_D", "f_d", "DDES", "DES length scale value");
+      AddVolumeOutput("L_RANS", "l_RANS", "DDES", "DES length scale value");
+      AddVolumeOutput("L_LES", "l_LES", "DDES", "DES length scale value");
+      AddVolumeOutput("R_D", "r_d", "DDES", "DES length scale value");
+    } else if ( config->GetKind_HybridRANSLES() == SST_IDDES){
+      AddVolumeOutput("F_D", "f_d", "DDES", "DES length scale value");
+      AddVolumeOutput("L_RANS", "l_RANS", "DDES", "DES length scale value");
+      AddVolumeOutput("L_LES", "l_LES", "DDES", "DES length scale value");
+      AddVolumeOutput("R_DT", "r_dt", "DDES", "DES length scale value");
+      AddVolumeOutput("R_DL", "r_dl", "DDES", "DES length scale value");
+    } else if ( config->GetKind_HybridRANSLES() == SST_SIDDES){
+      AddVolumeOutput("F_D", "f_d", "DDES", "DES length scale value");
+      AddVolumeOutput("L_RANS", "l_RANS", "DDES", "DES length scale value");
+      AddVolumeOutput("L_LES", "l_LES", "DDES", "DES length scale value");
+      AddVolumeOutput("R_DT", "r_dt", "DDES", "DES length scale value");
+    }
+    AddVolumeOutput("LESIQ", "LESIQ", "DDES", "LESIQ index for SRS simulations");
   }
 
   if (config->GetViscous()) {
@@ -1563,6 +1581,27 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
   if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES) {
     SetVolumeOutputValue("DES_LENGTHSCALE", iPoint, Node_Flow->GetDES_LengthScale(iPoint));
     SetVolumeOutputValue("WALL_DISTANCE", iPoint, Node_Geo->GetWall_Distance(iPoint));
+    if ( config->GetKind_HybridRANSLES() == SST_DDES){
+      SetVolumeOutputValue("F_D", iPoint, Node_Turb->Get_ftilda_d(iPoint));
+      SetVolumeOutputValue("L_RANS", iPoint, Node_Turb->Get_L_RANS(iPoint));
+      SetVolumeOutputValue("L_LES", iPoint, Node_Turb->Get_L_LES(iPoint));
+      SetVolumeOutputValue("R_D", iPoint, Node_Turb->Get_r_d(iPoint));
+    } else if ( config->GetKind_HybridRANSLES() == SST_IDDES){
+      SetVolumeOutputValue("F_D", iPoint, Node_Turb->Get_ftilda_d(iPoint));
+      SetVolumeOutputValue("L_RANS", iPoint, Node_Turb->Get_L_RANS(iPoint));
+      SetVolumeOutputValue("L_LES", iPoint, Node_Turb->Get_L_LES(iPoint));
+      SetVolumeOutputValue("R_DT", iPoint, Node_Turb->Get_r_dt(iPoint));
+      SetVolumeOutputValue("R_DL", iPoint, Node_Turb->Get_r_dl(iPoint));
+    } else if ( config->GetKind_HybridRANSLES() == SST_SIDDES){
+      SetVolumeOutputValue("F_D", iPoint, Node_Turb->Get_ftilda_d(iPoint));
+      SetVolumeOutputValue("L_RANS", iPoint, Node_Turb->Get_L_RANS(iPoint));
+      SetVolumeOutputValue("L_LES", iPoint, Node_Turb->Get_L_LES(iPoint));
+      SetVolumeOutputValue("R_DT", iPoint, Node_Turb->Get_r_dt(iPoint));
+    }
+    const su2double mut = Node_Flow->GetEddyViscosity(iPoint);
+    const su2double mu = Node_Flow->GetLaminarViscosity(iPoint); 
+    const su2double LESIQ = 1.0/(1.0+0.05*pow((mut+mu)/mu, 0.53));
+    SetVolumeOutputValue("LESIQ", iPoint, LESIQ);
   }
 
   switch (config->GetKind_Species_Model()) {