- Added SLA DDES formulation to SST

Author: rois1995

Common/include/option_structure.hpp

@@ -1450,7 +1450,8 @@ enum ENUM_HYBRIDRANSLES {
   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_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)
@@ -1461,6 +1462,7 @@ static const MapType<std::string, ENUM_HYBRIDRANSLES> HybridRANSLES_Map = {
   MakePair("SST_DDES", SST_DDES)
   MakePair("SST_IDDES", SST_IDDES)
   MakePair("SST_SIDDES", SST_SIDDES)
+  MakePair("SST_EDDES", SST_EDDES)
 };
 
 /*!

SU2_CFD/include/variables/CTurbSAVariable.hpp

@@ -39,7 +39,6 @@
 class CTurbSAVariable final : public CTurbVariable {
 
 private:
-  VectorType Vortex_Tilting;
 
   VectorType k, Omega;  /*!< \brief SST variables as computed through SA solution. */
 
@@ -61,18 +60,4 @@ class CTurbSAVariable final : public CTurbVariable {
    */
   ~CTurbSAVariable() override = default;
 
-  /*!
-   * \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

@@ -45,6 +45,7 @@ class CTurbVariable : public CScalarVariable {
   VectorType turb_index;            /*!< \brief Value of the turbulence index for transition simulations. */
   VectorType intermittency;         /*!< \brief Value of the intermittency for the transition model. */
   VectorType SRSGridSize;           /*!< \brief alue of the desired grid size for Scale Resolving Simulations. */
+  VectorType Vortex_Tilting;
 
   /*!
    * \brief Constructor of the class.
@@ -129,5 +130,18 @@ class CTurbVariable : public CScalarVariable {
    */
   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/src/solvers/CTurbSSTSolver.cpp

@@ -205,6 +205,21 @@ void CTurbSSTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contain
 
   if (kind_hybridRANSLES != NO_HYBRIDRANSLES) {
 
+    /*--- Set the vortex tilting coefficient at every node if required ---*/
+
+    if (kind_hybridRANSLES == SST_EDDES){
+      auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
+
+      SU2_OMP_FOR_STAT(omp_chunk_size)
+      for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){
+        auto Vorticity = flowNodes->GetVorticity(iPoint);
+        auto PrimGrad_Flow = flowNodes->GetGradient_Primitive(iPoint);
+        auto Laminar_Viscosity = flowNodes->GetLaminarViscosity(iPoint);
+        nodes->SetVortex_Tilting(iPoint, PrimGrad_Flow, Vorticity, Laminar_Viscosity);
+      }
+      END_SU2_OMP_FOR
+    }
+
     /*--- Compute the DES length scale ---*/
 
     SetDES_LengthScale(solver_container, geometry, config);
@@ -1154,7 +1169,8 @@ void CTurbSSTSolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, C
   for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){
 
     const su2double StrainMag = max(flowNodes->GetStrainMag(iPoint), 1e-12);
-    const su2double VortMag = max(GeometryToolbox::Norm(3, flowNodes->GetVorticity(iPoint)), 1e-12);
+    const auto Vorticity      = flowNodes->GetVorticity(iPoint);
+    const su2double VortMag   = max(GeometryToolbox::Norm(3, Vorticity), 1e-12);
 
     const su2double KolmConst2 = 0.41*0.41;
     const su2double wallDist2 = geometry->nodes->GetWall_Distance(iPoint)*geometry->nodes->GetWall_Distance(iPoint); 
@@ -1242,6 +1258,65 @@ void CTurbSSTSolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, C
 
         break;
       }
+      case SST_EDDES: {
+
+        // Improved DDES version with the Shear-Layer-Adapted augmentation 
+        // found in Detached Eddy Simulation: Recent Development and Application to Compressor Tip Leakage Flow, Xiao He, Fanzhou Zhao, Mehdi Vahdati
+        // originally from Application of SST-Based SLA-DDES Formulation to Turbomachinery Flows, Guoping Xia, Zifei Yin and Gorazd Medic
+        // I could be naming it either as SST_EDDES to follow the same notation as for the SA model or as SST_SLA_DDES to follow the paper notation
+
+        const su2double f_max = 1.0, f_min = 0.1, a1 = 0.15, a2 = 0.3;
+
+        const auto nNeigh = geometry->nodes->GetnPoint(iPoint);
+
+        su2double vortexTiltingMeasure = nodes->GetVortex_Tilting(iPoint);
+
+        su2double deltaOmega = -1;
+        su2double vorticityDir[MAXNDIM] = {};
+
+        for (auto iDim = 0; iDim < 3; iDim++){
+          vorticityDir[iDim] = Vorticity[iDim]/VortMag;
+        }
+        
+        for (const auto jPoint : geometry->nodes->GetPoints(iPoint)){
+          const auto coord_j = geometry->nodes->GetCoord(jPoint);
+
+          for (const auto kPoint : geometry->nodes->GetPoints(iPoint)){
+            const auto coord_k = geometry->nodes->GetCoord(kPoint);
+
+            su2double delta[3];
+            for (auto iDim = 0u; iDim < nDim; iDim++){
+              delta[iDim] = (coord_j[iDim] - coord_k[iDim])/2.0; // Should I divide by 2 as I am interested in the dual volume?
+            }
+            su2double l_n_minus_m[3];
+            GeometryToolbox::CrossProduct(delta, vorticityDir, l_n_minus_m);
+            deltaOmega = max(deltaOmega, GeometryToolbox::Norm(nDim, l_n_minus_m));
+          }
+
+          // Add to VTM(iPoint) to perform the average
+          vortexTiltingMeasure += nodes->GetVortex_Tilting(jPoint);
+        }
+        deltaOmega /= (1/sqrt(3));
+        vortexTiltingMeasure /= (nNeigh+1);
+
+        const su2double f_kh = max(f_min,
+                                   min(f_max,
+                                       f_min + ((f_max - f_min)/(a2 - a1)) * (vortexTiltingMeasure - a1)));
+
+        const su2double r_d = (eddyVisc + lamVisc) / max((KolmConst2*wallDist2 * sqrt(0.5 * (StrainMag*StrainMag + VortMag*VortMag))), 1e-10);
+        const su2double C_d1 = 20.0;
+        const su2double C_d2 = 3.0;
+
+        const su2double f_d = 1 - tanh(pow(C_d1 * r_d, C_d2));
+
+        su2double delta = deltaOmega * f_kh;
+        if (f_d < 0.99){
+          delta = h_max;
+        }
+
+        const su2double l_LES = C_DES * delta;
+        DES_lengthScale = l_RANS - f_d * max(0.0, l_RANS - l_LES);
+      }
     }
 
     nodes->SetDES_LengthScale(iPoint, DES_lengthScale);

SU2_CFD/src/variables/CTurbSAVariable.cpp

@@ -46,55 +46,7 @@ CTurbSAVariable::CTurbSAVariable(su2double val_nu_tilde, su2double val_muT, unsi
     Solution_time_n1 = Solution;
   }
 
-  Vortex_Tilting.resize(nPoint);
 
   k.resize(nPoint) = su2double(0.0);
   Omega.resize(nPoint) = su2double(0.0);
-}
-
-void CTurbSAVariable::SetVortex_Tilting(unsigned long iPoint, CMatrixView<const su2double> PrimGrad_Flow,
-                                        const su2double* Vorticity, su2double LaminarViscosity) {
-
-  su2double Strain[3][3] = {{0,0,0}, {0,0,0}, {0,0,0}}, Omega, StrainDotVort[3], numVecVort[3];
-  su2double numerator, trace0, trace1, denominator;
-
-  AD::StartPreacc();
-  AD::SetPreaccIn(PrimGrad_Flow, nDim+1, nDim);
-  AD::SetPreaccIn(Vorticity, 3);
-  /*--- Eddy viscosity ---*/
-  AD::SetPreaccIn(muT(iPoint));
-  /*--- Laminar viscosity --- */
-  AD::SetPreaccIn(LaminarViscosity);
-
-  Strain[0][0] = PrimGrad_Flow[1][0];
-  Strain[1][0] = 0.5*(PrimGrad_Flow[2][0] + PrimGrad_Flow[1][1]);
-  Strain[0][1] = 0.5*(PrimGrad_Flow[1][1] + PrimGrad_Flow[2][0]);
-  Strain[1][1] = PrimGrad_Flow[2][1];
-  if (nDim == 3){
-    Strain[0][2] = 0.5*(PrimGrad_Flow[3][0] + PrimGrad_Flow[1][2]);
-    Strain[1][2] = 0.5*(PrimGrad_Flow[3][1] + PrimGrad_Flow[2][2]);
-    Strain[2][0] = 0.5*(PrimGrad_Flow[1][2] + PrimGrad_Flow[3][0]);
-    Strain[2][1] = 0.5*(PrimGrad_Flow[2][2] + PrimGrad_Flow[3][1]);
-    Strain[2][2] = PrimGrad_Flow[3][2];
-  }
-
-  Omega = sqrt(Vorticity[0]*Vorticity[0] + Vorticity[1]*Vorticity[1]+ Vorticity[2]*Vorticity[2]);
-
-  StrainDotVort[0] = Strain[0][0]*Vorticity[0]+Strain[0][1]*Vorticity[1]+Strain[0][2]*Vorticity[2];
-  StrainDotVort[1] = Strain[1][0]*Vorticity[0]+Strain[1][1]*Vorticity[1]+Strain[1][2]*Vorticity[2];
-  StrainDotVort[2] = Strain[2][0]*Vorticity[0]+Strain[2][1]*Vorticity[1]+Strain[2][2]*Vorticity[2];
-
-  numVecVort[0] = StrainDotVort[1]*Vorticity[2] - StrainDotVort[2]*Vorticity[1];
-  numVecVort[1] = StrainDotVort[2]*Vorticity[0] - StrainDotVort[0]*Vorticity[2];
-  numVecVort[2] = StrainDotVort[0]*Vorticity[1] - StrainDotVort[1]*Vorticity[0];
-
-  numerator = sqrt(6.0) * sqrt(numVecVort[0]*numVecVort[0] + numVecVort[1]*numVecVort[1] + numVecVort[2]*numVecVort[2]);
-  trace0 = 3.0*(pow(Strain[0][0],2.0) + pow(Strain[1][1],2.0) + pow(Strain[2][2],2.0));
-  trace1 = pow(Strain[0][0] + Strain[1][1] + Strain[2][2],2.0);
-  denominator = pow(Omega, 2.0) * sqrt(trace0-trace1);
-
-  Vortex_Tilting(iPoint) = (numerator/denominator) * max(1.0,0.2*LaminarViscosity/muT(iPoint));
-
-  AD::SetPreaccOut(Vortex_Tilting(iPoint));
-  AD::EndPreacc();
-}
+}

SU2_CFD/src/variables/CTurbVariable.cpp

@@ -37,5 +37,55 @@ CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned
 
     DES_LengthScale.resize(nPoint) = su2double(0.0);
     SRSGridSize.resize(nPoint) = su2double(0.0);
+    
+    Vortex_Tilting.resize(nPoint);
 
    }
+
+
+void CTurbVariable::SetVortex_Tilting(unsigned long iPoint, CMatrixView<const su2double> PrimGrad_Flow,
+                                        const su2double* Vorticity, su2double LaminarViscosity) {
+
+  su2double Strain[3][3] = {{0,0,0}, {0,0,0}, {0,0,0}}, Omega, StrainDotVort[3], numVecVort[3];
+  su2double numerator, trace0, trace1, denominator;
+
+  AD::StartPreacc();
+  AD::SetPreaccIn(PrimGrad_Flow, nDim+1, nDim);
+  AD::SetPreaccIn(Vorticity, 3);
+  /*--- Eddy viscosity ---*/
+  AD::SetPreaccIn(muT(iPoint));
+  /*--- Laminar viscosity --- */
+  AD::SetPreaccIn(LaminarViscosity);
+
+  Strain[0][0] = PrimGrad_Flow[1][0];
+  Strain[1][0] = 0.5*(PrimGrad_Flow[2][0] + PrimGrad_Flow[1][1]);
+  Strain[0][1] = 0.5*(PrimGrad_Flow[1][1] + PrimGrad_Flow[2][0]);
+  Strain[1][1] = PrimGrad_Flow[2][1];
+  if (nDim == 3){
+    Strain[0][2] = 0.5*(PrimGrad_Flow[3][0] + PrimGrad_Flow[1][2]);
+    Strain[1][2] = 0.5*(PrimGrad_Flow[3][1] + PrimGrad_Flow[2][2]);
+    Strain[2][0] = 0.5*(PrimGrad_Flow[1][2] + PrimGrad_Flow[3][0]);
+    Strain[2][1] = 0.5*(PrimGrad_Flow[2][2] + PrimGrad_Flow[3][1]);
+    Strain[2][2] = PrimGrad_Flow[3][2];
+  }
+
+  Omega = sqrt(Vorticity[0]*Vorticity[0] + Vorticity[1]*Vorticity[1]+ Vorticity[2]*Vorticity[2]);
+
+  StrainDotVort[0] = Strain[0][0]*Vorticity[0]+Strain[0][1]*Vorticity[1]+Strain[0][2]*Vorticity[2];
+  StrainDotVort[1] = Strain[1][0]*Vorticity[0]+Strain[1][1]*Vorticity[1]+Strain[1][2]*Vorticity[2];
+  StrainDotVort[2] = Strain[2][0]*Vorticity[0]+Strain[2][1]*Vorticity[1]+Strain[2][2]*Vorticity[2];
+
+  numVecVort[0] = StrainDotVort[1]*Vorticity[2] - StrainDotVort[2]*Vorticity[1];
+  numVecVort[1] = StrainDotVort[2]*Vorticity[0] - StrainDotVort[0]*Vorticity[2];
+  numVecVort[2] = StrainDotVort[0]*Vorticity[1] - StrainDotVort[1]*Vorticity[0];
+
+  numerator = sqrt(6.0) * sqrt(numVecVort[0]*numVecVort[0] + numVecVort[1]*numVecVort[1] + numVecVort[2]*numVecVort[2]);
+  trace0 = 3.0*(pow(Strain[0][0],2.0) + pow(Strain[1][1],2.0) + pow(Strain[2][2],2.0));
+  trace1 = pow(Strain[0][0] + Strain[1][1] + Strain[2][2],2.0);
+  denominator = pow(Omega, 2.0) * sqrt(trace0-trace1);
+
+  Vortex_Tilting(iPoint) = (numerator/denominator) * max(1.0,0.2*LaminarViscosity/muT(iPoint));
+
+  AD::SetPreaccOut(Vortex_Tilting(iPoint));
+  AD::EndPreacc();
+}