[WIP] Implementation of DDES and SAS formulations for SST

Common/include/option_structure.hpp

@@ -990,6 +990,8 @@ enum class SST_OPTIONS {
   V,           /*!< \brief Menter k-w SST model with vorticity production terms. */
   KL,          /*!< \brief Menter k-w SST model with Kato-Launder production terms. */
   UQ,          /*!< \brief Menter k-w SST model with uncertainty quantification modifications. */
+  SAS_SIMPLE,         /*!< \brief Menter k-w SST model with Scale Adaptive Simulations modifications. */
+  SAS_COMPLICATED,         /*!< \brief Menter k-w SST model with Scale Adaptive Simulations modifications. */
 };
 static const MapType<std::string, SST_OPTIONS> SST_Options_Map = {
   MakePair("NONE", SST_OPTIONS::NONE)
@@ -1002,6 +1004,8 @@ static const MapType<std::string, SST_OPTIONS> SST_Options_Map = {
   MakePair("VORTICITY", SST_OPTIONS::V)
   MakePair("KATO-LAUNDER", SST_OPTIONS::KL)
   MakePair("UQ", SST_OPTIONS::UQ)
+  MakePair("SAS_SIMPLE", SST_OPTIONS::SAS_SIMPLE)
+  MakePair("SAS_COMPLICATED", SST_OPTIONS::SAS_COMPLICATED)
 };
 
 /*!
@@ -1010,8 +1014,10 @@ static const MapType<std::string, SST_OPTIONS> SST_Options_Map = {
 struct SST_ParsedOptions {
   SST_OPTIONS version = SST_OPTIONS::V1994;   /*!< \brief Enum SST base model. */
   SST_OPTIONS production = SST_OPTIONS::NONE; /*!< \brief Enum for production corrections/modifiers for SST model. */
+  SST_OPTIONS sasModel = SST_OPTIONS::SAS_SIMPLE;   /*!< \brief Enum SST base model. */
   bool sust = false;                          /*!< \brief Bool for SST model with sustaining terms. */
   bool uq = false;                            /*!< \brief Bool for using uncertainty quantification. */
+  bool sas = false;                           /*!< \brief Bool for using Scale Adaptive Simulations. */
   bool modified = false;                      /*!< \brief Bool for modified (m) SST model. */
 };
 
@@ -1048,6 +1054,16 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
   const bool sst_v = IsPresent(SST_OPTIONS::V);
   const bool sst_kl = IsPresent(SST_OPTIONS::KL);
   const bool sst_uq = IsPresent(SST_OPTIONS::UQ);
+  const bool sst_sas_simple = IsPresent(SST_OPTIONS::SAS_SIMPLE);
+  const bool sst_sas_comp = IsPresent(SST_OPTIONS::SAS_COMPLICATED);
+  const bool sst_sas = sst_sas_simple || sst_sas_comp;
+  if (sst_sas_simple && sst_sas_comp) {
+    SU2_MPI::Error("Two versions (Simple and Complicated) selected for SAS under SST_OPTIONS. Please choose only one.", CURRENT_FUNCTION);
+  } else if (sst_sas_simple) {
+    SSTParsedOptions.sasModel = SST_OPTIONS::SAS_SIMPLE;
+  } else {
+    SSTParsedOptions.sasModel = SST_OPTIONS::SAS_COMPLICATED;
+  } 
 
   if (sst_1994 && sst_2003) {
     SU2_MPI::Error("Two versions (1994 and 2003) selected for SST_OPTIONS. Please choose only one.", CURRENT_FUNCTION);
@@ -1071,6 +1087,7 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
   SSTParsedOptions.sust = sst_sust;
   SSTParsedOptions.modified = sst_m;
   SSTParsedOptions.uq = sst_uq;
+  SSTParsedOptions.sas = sst_sas;
   return SSTParsedOptions;
 }
 
@@ -1380,14 +1397,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)
 };
 
 /*!
@@ -2450,6 +2473,7 @@ enum PERIODIC_QUANTITIES {
   PERIODIC_LIM_PRIM_1 ,  /*!< \brief Primitive limiter communication phase 1 of 2 (periodic only). */
   PERIODIC_LIM_PRIM_2 ,  /*!< \brief Primitive limiter communication phase 2 of 2 (periodic only). */
   PERIODIC_IMPLICIT   ,  /*!< \brief Implicit update communication to ensure consistency across periodic boundaries. */
+  PERIODIC_VEL_LAPLACIAN  ,  /*!< \brief Velocity Laplacian communication for SAS (periodic only). */
 };
 
 /*!
@@ -2481,6 +2505,7 @@ enum MPI_QUANTITIES {
   MESH_DISPLACEMENTS   ,  /*!< \brief Mesh displacements at the interface. */
   SOLUTION_TIME_N      ,  /*!< \brief Solution at time n. */
   SOLUTION_TIME_N1     ,  /*!< \brief Solution at time n-1. */
+  VELOCITY_LAPLACIAN   ,  /*!< \brief Velocity Laplacian communication. */
 };
 
 /*!

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,10 @@ class CNumerics {
 
   bool bounded_scalar = false;    /*!< \brief Flag for bounded scalar problem */
 
+  su2double lengthScale_i, lengthScale_j;
+  su2double FTrans;   /*!< \brief SAS function */
+  su2double VelLapl_X, VelLapl_Y, VelLapl_Z;
+
 public:
   /*!
    * \brief Return type used in some "ComputeResidual" overloads to give a
@@ -707,6 +711,26 @@ class CNumerics {
    */
   virtual void SetCrossDiff(su2double val_CDkw_i) {/* empty */};
 
+  /*!
+   * \brief Get the value of the value of FTrans.
+   */
+  inline virtual su2double GetFTrans() const { return 0.0; }
+
+  /*!
+   * \brief Get the value of the value of FTrans.
+   */
+  inline void SetVelLapl(su2double val_VelLapl_X, su2double val_VelLapl_Y) {
+    VelLapl_X = val_VelLapl_X;
+    VelLapl_Y = val_VelLapl_Y;
+  }
+
+  /*!
+   * \brief Get the value of the value of FTrans.
+   */
+  inline void SetVelLapl_Z(su2double val_VelLapl_Z) {
+    VelLapl_Z = val_VelLapl_Z;
+  }
+
   /*!
    * \brief Set the value of the effective intermittency for the LM model.
    * \param[in] intermittency_eff_i - Value of the effective intermittency at point i.
@@ -828,6 +852,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

@@ -605,6 +605,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
   /*--- Closure constants ---*/
   const su2double sigma_k_1, sigma_k_2, sigma_w_1, sigma_w_2, beta_1, beta_2, beta_star, a1, alfa_1, alfa_2;
   const su2double prod_lim_const;
+  const su2double cTrans;
 
   /*--- Ambient values for SST-SUST. ---*/
   const su2double kAmb, omegaAmb;
@@ -689,6 +690,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
         alfa_1(constants[8]),
         alfa_2(constants[9]),
         prod_lim_const(constants[10]),
+        cTrans(1.25),
         kAmb(val_kine_Inf),
         omegaAmb(val_omega_Inf) {
     /*--- "Allocate" the Jacobian using the static buffer. ---*/
@@ -834,6 +836,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---*/
@@ -847,9 +852,53 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       Residual[0] += pk * Volume;
       Residual[1] += pw * Volume;
 
-      /*--- Add the dissipation  terms to the residuals.---*/
+      /*--- Compute SAS source term. ---*/
+      su2double Q_SAS = 0.0;
+
+      if (sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_COMPLICATED) {
+
+        const su2double KolmConst = 0.41;
+        const su2double csi2 = 3.51;
+        const su2double sigma_phi = 2.0/3.0;
+        const su2double C = 2.0;
+        const su2double C_s = 0.5; // Honestly I do not know if it is the right one.
+        const su2double gridSize = pow(Volume, 1.0/3.0);
+        // Scale of the modeled turbulence
+        const su2double L = sqrt(ScalarVar_i[0]) / (pow(beta_star, 0.25) * ScalarVar_i[1]);
+        // Von Karman Length Scale
+        su2double VelLaplMag = VelLapl_X*VelLapl_X + VelLapl_Y*VelLapl_Y;
+        if (nDim == 3) VelLaplMag += VelLapl_Z*VelLapl_Z;
+        const su2double L_vK_1 = KolmConst * StrainMag_i / sqrt(VelLaplMag);
+        const su2double L_vK_2 = C_s * sqrt(KolmConst * csi2 / (beta_blended/beta_star - alfa_blended)) * gridSize;
+        const su2double L_vK = max(L_vK_1, L_vK_2);
+
+        su2double gradOmega = 0.0;
+        for (unsigned short iDim = 0; iDim < nDim; iDim++)
+          gradOmega += ScalarVar_Grad_i[1][iDim]*ScalarVar_Grad_i[1][iDim];
+
+        gradOmega /= ScalarVar_i[1];
+
+        su2double gradTKE = 0.0;
+        for (unsigned short iDim = 0; iDim < nDim; iDim++)
+          gradTKE += ScalarVar_Grad_i[0][iDim]*ScalarVar_Grad_i[0][iDim];
+
+        gradTKE /= ScalarVar_i[0];
+
+        const su2double Q_SAS_1 = csi2 * KolmConst * StrainMag_i * StrainMag_i * (L/L_vK) * (L/L_vK);
+        const su2double Q_SAS_2 = C * (2*ScalarVar_i[0] / sigma_phi) * max(gradOmega, gradTKE);
+
+        Q_SAS = max(Q_SAS_1 - Q_SAS_2, 0.0);
+
+        Residual[1] += Density_i * Q_SAS * Volume;
+
+      }
 
-      Residual[0] -= dk * Volume;
+      /*--- Add the dissipation  terms to the residuals.---*/
+      FTrans = 1.0;
+      if (sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_SIMPLE) {
+        FTrans = max(1.0, pow(StrainMag_i / (cTrans * VorticityMag), 2.0));
+      }
+      Residual[0] -= dk * Volume * FTrans;
       Residual[1] -= dw * Volume;
 
       /*--- Cross diffusion ---*/
@@ -862,15 +911,26 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
 
       /*--- Implicit part ---*/
 
-      Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume;
-      Jacobian_i[0][1] = -beta_star * ScalarVar_i[0] * Volume;
+      Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume * FTrans;
+      Jacobian_i[0][1] = -beta_star * ScalarVar_i[0] * Volume * FTrans;
       Jacobian_i[1][0] = 0.0;
       Jacobian_i[1][1] = -2.0 * beta_blended * ScalarVar_i[1] * Volume;
+
+      if (sstParsedOptions.sas && sstParsedOptions.sasModel == SST_OPTIONS::SAS_COMPLICATED) {
+        Jacobian_i[0][0] += Q_SAS * Volume / ScalarVar_i[0];
+      }
     }
 
     AD::SetPreaccOut(Residual, nVar);
     AD::EndPreacc();
 
     return ResidualType<>(Residual, Jacobian_i, nullptr);
   }
+
+  /*!
+   * \brief Get the value of the FTrans.
+   */
+  inline su2double GetFTrans() const override { return FTrans; }
+
+
 };

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,9 +39,10 @@
 class CTurbSAVariable final : public CTurbVariable {
 
 private:
-  VectorType DES_LengthScale;
   VectorType Vortex_Tilting;
 
+  VectorType k, Omega;  /*!< \brief SST variables as computed through SA solution. */
+
 public:
   /*!
    * \brief Constructor of the class.
@@ -60,19 +61,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.
@@ -87,4 +75,26 @@ class CTurbSAVariable final : public CTurbVariable {
    */
   inline su2double GetVortex_Tilting(unsigned long iPoint) const override { return Vortex_Tilting(iPoint); }
 
+  /*!
+   * \brief Get the value of the turbulence kinetic energy.
+   * \return the value of the turbulence kinetic energy.
+   */
+  inline su2double GetSSTVariables_k(unsigned long iPoint) const { return k(iPoint); }
+
+  /*!
+   * \brief Get the value of the turbulence frequency Omega.
+   * \return the value of the turbulence frequency Omega.
+   */
+  inline su2double GetSSTVariables_omega(unsigned long iPoint) const { return Omega(iPoint); }
+
+  /*!
+   * \brief Set the value of the SST variables computed with SA solution.
+   * \param[in] val_k
+   * \param[in] val_Omega
+   */
+  void SetSSTVariables(unsigned long iPoint, su2double val_k, su2double val_Omega) {
+    k(iPoint) = val_k;
+    Omega(iPoint) = val_Omega;
+  }
+
 };