[WIP] Implementation of DDES and SAS formulations for SST

Common/include/option_structure.hpp

@@ -994,6 +994,8 @@ enum class SST_OPTIONS {
   COMP_Wilcox, /*!< \brief Menter k-w SST model with Compressibility correction of Wilcox. */
   COMP_Sarkar, /*!< \brief Menter k-w SST model with Compressibility correction of Sarkar. */
   DLL,         /*!< \brief Menter k-w SST model with dimensionless lower limit clipping of turbulence variables. */
+  SAS_TRAVIS,         /*!< \brief Menter k-w SST model with Scale Adaptive Simulations modifications. */
+  SAS_BABU,         /*!< \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)
@@ -1009,6 +1011,8 @@ static const MapType<std::string, SST_OPTIONS> SST_Options_Map = {
   MakePair("COMPRESSIBILITY-WILCOX", SST_OPTIONS::COMP_Wilcox)
   MakePair("COMPRESSIBILITY-SARKAR", SST_OPTIONS::COMP_Sarkar)
   MakePair("DIMENSIONLESS_LIMIT", SST_OPTIONS::DLL)
+  MakePair("SAS_TRAVIS", SST_OPTIONS::SAS_TRAVIS)
+  MakePair("SAS_BABU", SST_OPTIONS::SAS_BABU)
 };
 
 /*!
@@ -1017,6 +1021,7 @@ 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::NONE;   /*!< \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 modified = false;                      /*!< \brief Bool for modified (m) SST model. */
@@ -1061,6 +1066,15 @@ 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);
+  const bool sst_sas_simple = IsPresent(SST_OPTIONS::SAS_TRAVIS);
+  const bool sst_sas_comp = IsPresent(SST_OPTIONS::SAS_BABU);
+  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_TRAVIS;
+  } else {
+    SSTParsedOptions.sasModel = SST_OPTIONS::SAS_BABU;
+  } 
 
   if (sst_1994 && sst_2003) {
     SU2_MPI::Error("Two versions (1994 and 2003) selected for SST_OPTIONS. Please choose only one.", CURRENT_FUNCTION);
@@ -1433,14 +1447,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)
 };
 
 /*!
@@ -2519,6 +2539,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). */
 };
 
 /*!
@@ -2550,6 +2571,7 @@ enum class 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

@@ -6225,6 +6225,17 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
             else cout << "\nusing default hard coded lower limit clipping";
 
             cout << "." << endl;
+            if(sstParsedOptions.sasModel != SST_OPTIONS::NONE) cout << "Scale Adaptive Simulation model: ";
+            switch (sstParsedOptions.sasModel) {
+              case SST_OPTIONS::SAS_TRAVIS:
+                cout << "Travis et al. (2004)";
+                break;
+              case SST_OPTIONS::SAS_BABU:
+                cout << "Babu et al. (2016)";
+                break;
+              default:
+                break;
+            }
             break;
         }
         switch (Kind_Trans_Model) {
@@ -6266,6 +6277,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)" << 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

@@ -136,6 +136,9 @@ 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 */
+  su2double FTrans;   /*!< \brief SAS function */
+  su2double *VelLapl; /*!< \brief Laplacian of the velocity */
   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. */
@@ -707,6 +710,18 @@ 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) {
+    VelLapl = val_VelLapl;
+  }
+
   /*!
    * \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 +843,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

@@ -654,6 +654,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;
@@ -738,6 +739,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. ---*/
@@ -905,6 +907,9 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       /*--- Dissipation ---*/
 
       su2double dk = beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * (1.0 + zetaFMt);
+      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] * (1.0 - 0.09/beta_blended * zetaFMt);
 
       /*--- LM model coupling with production and dissipation term for k transport equation---*/
@@ -918,9 +923,43 @@ 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.sasModel == SST_OPTIONS::SAS_BABU) {
+
+        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
+        const su2double VelLaplMag = GeometryToolbox::SquaredNorm(nDim, VelLapl);
+        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);
+
+        const su2double gradTKE = GeometryToolbox::SquaredNorm(nDim, ScalarVar_Grad_i[0])/ScalarVar_i[0];
+        const su2double gradOmega = GeometryToolbox::SquaredNorm(nDim, ScalarVar_Grad_i[1])/ScalarVar_i[1];
+
+        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.sasModel == SST_OPTIONS::SAS_TRAVIS) {
+        FTrans = max(1.0, pow(StrainMag_i / (cTrans * VorticityMag), 2.0));
+      }
+      Residual[0] -= dk * Volume * FTrans;
       Residual[1] -= dw * Volume;
 
       /*--- Cross diffusion ---*/
@@ -933,15 +972,26 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
 
       /*--- Implicit part ---*/
 
-      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);
+      Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume * (1.0 + zetaFMt) * FTrans;
+      Jacobian_i[0][1] = -beta_star * ScalarVar_i[0] * Volume * (1.0 + zetaFMt) * FTrans;
       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);
+
+      if (sstParsedOptions.sasModel == SST_OPTIONS::SAS_BABU) {
+        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.