Implementation of DDES and SAS formulations for SST

Common/include/option_structure.hpp

@@ -1062,7 +1062,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) {
@@ -1473,17 +1472,29 @@ 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). */
+  SA_EDDES_UNSTR   = 5,  /*!< \brief Kind of Hybrid RANS/LES (SA - Delayed DES (DDES) with Shear Layer Adapted SGS: Enhanced DDES and for Unstructured grids). */
+  SST_DDES         = 6,  /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): DDES). */
+  SST_IDDES        = 7,  /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Improved DDES). */
+  SST_SIDDES       = 8,  /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Simplified Improved DDES). */
+  SST_EDDES        = 9,  /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Enhanced (SLA) DDES). */
+  SST_EDDES_UNSTR  = 10  /*!< \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("SA_EDDES_UNSTR", SA_EDDES_UNSTR)
+  MakePair("SST_DDES", SST_DDES)
+  MakePair("SST_IDDES", SST_IDDES)
+  MakePair("SST_SIDDES", SST_SIDDES)
+  MakePair("SST_EDDES", SST_EDDES)
+  MakePair("SST_EDDES_UNSTR", SST_EDDES_UNSTR)
 };
 
 /*!

Common/src/CConfig.cpp

@@ -6442,10 +6442,14 @@ 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 SA_EDDES_UNSTR: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS (modified for Unstructured grids)" << 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;
         }
         break;
       case MAIN_SOLVER::NEMO_EULER:

SU2_CFD/include/numerics/CNumerics.hpp

@@ -135,6 +135,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. */
@@ -827,6 +828,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

@@ -885,14 +885,17 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       }
 
       if (sstParsedOptions.production == SST_OPTIONS::COMP_Sarkar) {
-        const su2double Dilatation_Sarkar = -0.15 * pk * Mt + 0.2 * beta_star * (1.0 +zetaFMt) * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * Mt * Mt;
+        const su2double Dilatation_Sarkar = -0.15 * pk * Mt + 0.2 * beta_star * (1.0 +zetaFMt) * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * pow(Mt, 2);
         pk += Dilatation_Sarkar;
       }
 
       /*--- Dissipation ---*/
 
       su2double dk = beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * (1.0 + zetaFMt);
-      su2double dw = beta_blended * Density_i * ScalarVar_i[1] * ScalarVar_i[1] * (1.0 - 0.09/beta_blended * zetaFMt);
+      if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES)
+        dk = Density_i * sqrt(pow(ScalarVar_i[0], 3)) / lengthScale_i;
+
+      su2double dw = beta_blended * Density_i * pow(ScalarVar_i[1], 2) * (1.0 - 0.09/beta_blended * zetaFMt);
 
       /*--- LM model coupling with production and dissipation term for k transport equation---*/
       if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) {

SU2_CFD/include/solvers/CTurbSSTSolver.hpp

@@ -61,6 +61,17 @@ class CTurbSSTSolver final : public CTurbSolver {
    */
   void ComputeUnderRelaxationFactor(const CConfig *config);
 
+
+   /*!
+   * \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,8 +39,6 @@
 class CTurbSAVariable final : public CTurbVariable {
 
 private:
-  VectorType DES_LengthScale;
-  VectorType Vortex_Tilting;
 
 public:
   /*!
@@ -60,31 +58,4 @@ 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.
-   */
-  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); }
-
 };

SU2_CFD/include/variables/CTurbVariable.hpp

@@ -37,12 +37,14 @@
  */
 class CTurbVariable : public CScalarVariable {
 protected:
-  VectorType muT; /*!< \brief Eddy viscosity. */
+  VectorType muT;             /*!< \brief Eddy viscosity. */
+  VectorType DES_LengthScale; /*!< \brief DES Length scale. */
 
 public:
   static constexpr size_t MAXNVAR = 2;
-  VectorType turb_index;
-  VectorType intermittency;         /*!< \brief Value of the intermittency for the trans. model. */
+  VectorType turb_index;            /*!< \brief Value of the turbulence index for transition simulations. */
+  VectorType intermittency;         /*!< \brief Value of the intermittency for the transition model. */
+  VectorType Vortex_Tilting;        /*!< \brief Value of the vortex tilting measure for EDDES models. */
 
   /*!
    * \brief Constructor of the class.
@@ -100,6 +102,32 @@ 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; }
+
+  /*!
+   * \brief Set the vortex tilting measure for computation of the EDDES length scale
+   * \param[in] iPoint - Point index.
+   */
+  void SetVortex_Tilting(unsigned long iPoint, su2double **Strain,
+                         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 Diffusion Coefficients of TKE and omega equations.
    * \param[in] iPoint - Point index.

SU2_CFD/include/variables/CVariable.hpp

@@ -2235,7 +2235,7 @@ class CVariable {
 
   inline virtual su2double GetTau_Wall(unsigned long iPoint) const { return 0.0; }
 
-  inline virtual void SetVortex_Tilting(unsigned long iPoint, CMatrixView<const su2double> PrimGrad_Flow,
+  inline virtual void SetVortex_Tilting(unsigned long iPoint, su2double **Strain,
                                         const su2double* Vorticity, su2double LaminarViscosity) {}
 
   inline virtual su2double GetVortex_Tilting(unsigned long iPoint) const { return 0.0; }

SU2_CFD/src/output/CFlowOutput.cpp

@@ -1483,6 +1483,9 @@ 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");
+    AddVolumeOutput("LESIQ", "LESIQ", "DDES", "LESIQ index for SRS simulations");
+    if (config->GetKind_Turb_Model() == TURB_MODEL::SST)
+      AddVolumeOutput("SRS_GRID_SIZE", "Srs_grid_size", "DDES", "desired grid size for Scale Resolving Simulations");
   }
 
   if (config->GetViscous()) {
@@ -1494,6 +1497,27 @@ void CFlowOutput::SetVolumeOutputFieldsScalarMisc(const CConfig* config) {
       AddVolumeOutput("VORTICITY", "Vorticity", "VORTEX_IDENTIFICATION", "Value of the vorticity");
     }
     AddVolumeOutput("Q_CRITERION", "Q_Criterion", "VORTEX_IDENTIFICATION", "Value of the Q-Criterion");
+    if (config->GetKind_HybridRANSLES() == NO_HYBRIDRANSLES) {
+      AddVolumeOutput("WALL_DISTANCE", "Wall_Distance", "DEBUG", "Wall distance value");
+    }
+    AddVolumeOutput("GRAD_VEL_XX", "Grad_Vel_xx", "VELOCITY_GRADIENT", "X-Gradient of U");
+    AddVolumeOutput("GRAD_VEL_XY", "Grad_Vel_xy", "VELOCITY_GRADIENT", "Y-Gradient of U");
+    AddVolumeOutput("GRAD_VEL_YX", "Grad_Vel_yx", "VELOCITY_GRADIENT", "X-Gradient of U");
+    AddVolumeOutput("GRAD_VEL_YY", "Grad_Vel_yy", "VELOCITY_GRADIENT", "Y-Gradient of V");
+    if (nDim == 3) {
+      AddVolumeOutput("GRAD_VEL_XZ", "Grad_Vel_xz", "VELOCITY_GRADIENT", "Z-Gradient of U");
+      AddVolumeOutput("GRAD_VEL_YZ", "Grad_Vel_yz", "VELOCITY_GRADIENT", "Z-Gradient of V");
+      AddVolumeOutput("GRAD_VEL_ZX", "Grad_Vel_zx", "VELOCITY_GRADIENT", "X-Gradient of W");
+      AddVolumeOutput("GRAD_VEL_ZY", "Grad_Vel_zy", "VELOCITY_GRADIENT", "Y-Gradient of W");
+      AddVolumeOutput("GRAD_VEL_ZZ", "Grad_Vel_zz", "VELOCITY_GRADIENT", "Z-Gradient of W");
+    }
+
+    if (config->GetKind_Turb_Model() == TURB_MODEL::SST){
+      AddVolumeOutput("CDkw", "CDkw", "SST_QUANTITIES", "Cross-Diffusion term");
+      AddVolumeOutput("F1", "F1", "SST_QUANTITIES", "F1 blending function");
+      AddVolumeOutput("F2", "F2", "SST_QUANTITIES", "F2 blending function");
+    }
+
   }
 
   // Timestep info
@@ -1519,14 +1543,35 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
   SetVolumeOutputValue("CFL", iPoint, Node_Flow->GetLocalCFL(iPoint));
 
   if (config->GetViscous()) {
+    const auto VelGrad = Node_Flow->GetVelocityGradient(iPoint);
     if (nDim == 3){
       SetVolumeOutputValue("VORTICITY_X", iPoint, Node_Flow->GetVorticity(iPoint)[0]);
       SetVolumeOutputValue("VORTICITY_Y", iPoint, Node_Flow->GetVorticity(iPoint)[1]);
       SetVolumeOutputValue("VORTICITY_Z", iPoint, Node_Flow->GetVorticity(iPoint)[2]);
     } else {
       SetVolumeOutputValue("VORTICITY", iPoint, Node_Flow->GetVorticity(iPoint)[2]);
     }
-    SetVolumeOutputValue("Q_CRITERION", iPoint, GetQCriterion(Node_Flow->GetVelocityGradient(iPoint)));
+    SetVolumeOutputValue("Q_CRITERION", iPoint, GetQCriterion(VelGrad));
+    SetVolumeOutputValue("WALL_DISTANCE", iPoint, Node_Geo->GetWall_Distance(iPoint));
+
+    SetVolumeOutputValue("GRAD_VEL_XX", iPoint, VelGrad(0,0));
+    SetVolumeOutputValue("GRAD_VEL_XY", iPoint, VelGrad(0,1));
+    SetVolumeOutputValue("GRAD_VEL_YX", iPoint, VelGrad(1,0));
+    SetVolumeOutputValue("GRAD_VEL_YY", iPoint, VelGrad(1,1));
+    if (nDim == 3) {
+      SetVolumeOutputValue("GRAD_VEL_XZ", iPoint, VelGrad(0,2));
+      SetVolumeOutputValue("GRAD_VEL_YZ", iPoint, VelGrad(1,2));
+      SetVolumeOutputValue("GRAD_VEL_ZX", iPoint, VelGrad(2,0));
+      SetVolumeOutputValue("GRAD_VEL_ZY", iPoint, VelGrad(2,1));
+      SetVolumeOutputValue("GRAD_VEL_ZZ", iPoint, VelGrad(2,2));
+    }
+
+    if (config->GetKind_Turb_Model() == TURB_MODEL::SST){
+      SetVolumeOutputValue("CDkw", iPoint, Node_Turb->GetCrossDiff(iPoint));
+      SetVolumeOutputValue("F1", iPoint, Node_Turb->GetF1blending(iPoint));
+      SetVolumeOutputValue("F2", iPoint, Node_Turb->GetF2blending(iPoint));
+    }
+
   }
 
   const bool limiter = (config->GetKind_SlopeLimit_Turb() != LIMITER::NONE);
@@ -1582,6 +1627,17 @@ 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));
+    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);
+    if (config->GetKind_Turb_Model() == TURB_MODEL::SST) {
+      const su2double betaStar = 0.09; // constants[6]
+      const su2double RANSLength = sqrt(Node_Flow->GetSolution(iPoint, 0)) / max(1e-20, (betaStar * Node_Flow->GetSolution(iPoint, 1)));
+      const su2double RatioL = 0.1;  // TODO:: it should be less or equal than 0.2 - 0.1. Should be taken as input from config?
+      const su2double SRSGridSize = RANSLength * RatioL;
+      SetVolumeOutputValue("SRS_GRID_SIZE", iPoint, SRSGridSize);
+    }
   }
 
   switch (config->GetKind_Species_Model()) {

SU2_CFD/src/solvers/CSolver.cpp

@@ -968,7 +968,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
             }
 
             break;
-
+            
           default:
             SU2_MPI::Error("Unrecognized quantity for periodic communication.",
                            CURRENT_FUNCTION);
@@ -1290,7 +1290,7 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry,
                 limiter(iPoint, iVar) = min(limiter(iPoint, iVar), bufDRecv[buf_offset+iVar]);
 
               break;
-
+              
             default:
 
               SU2_MPI::Error("Unrecognized quantity for periodic communication.",