- Try changes to SA-EDDES

Author: rois1995

Common/include/option_structure.hpp

@@ -1494,14 +1494,16 @@ enum ENUM_HYBRIDRANSLES {
   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_EDDES = 8           /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Enhanced (SLA) DDES). */
+  SST_EDDES = 8,           /*!< \brief Kind of Hybrid RANS/LES (SST - Delayed DES (DDES): Enhanced (SLA) DDES). */
+  SA_EDDES_MOD = 9           /*!< \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_MOD", SA_EDDES_MOD)
   MakePair("SST_DDES", SST_DDES)
   MakePair("SST_IDDES", SST_IDDES)
   MakePair("SST_SIDDES", SST_SIDDES)

SU2_CFD/src/solvers/CTurbSASolver.cpp

@@ -184,7 +184,7 @@ void CTurbSASolver::Preprocessing(CGeometry *geometry, CSolver **solver_containe
 
     /*--- Set the vortex tilting coefficient at every node if required ---*/
 
-    if (kind_hybridRANSLES == SA_EDDES){
+    if (kind_hybridRANSLES == SA_EDDES || kind_hybridRANSLES == SA_EDDES_MOD){
       auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
 
       SU2_OMP_FOR_STAT(omp_chunk_size)
@@ -1511,6 +1511,62 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC
         const su2double distDES = constDES * maxDelta;
         lengthScale = wallDistance-f_d*max(0.0,(wallDistance-distDES));
 
+        break;
+      }
+      case SA_EDDES_MOD: {
+        /*--- An Enhanced Version of DES with Rapid Transition from RANS to LES in Separated Flows.
+         Shur et al.
+         Flow Turbulence Combust - 2015
+         ---*/
+
+        su2double vortexTiltingMeasure = nodes->GetVortex_Tilting(iPoint);
+
+        const su2double omega = GeometryToolbox::Norm(3, vorticity);
+
+        su2double ratioOmega[MAXNDIM] = {};
+
+        for (auto iDim = 0; iDim < 3; iDim++){
+          ratioOmega[iDim] = vorticity[iDim]/omega;
+        }
+
+        const su2double deltaDDES = geometry->nodes->GetMaxLength(iPoint);
+
+        su2double ln_max = -1;
+        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[MAXNDIM] = {};
+            // This should only be performed on 3D cases anyway
+            for (auto iDim = 0u; iDim < 3; 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, ratioOmega, l_n_minus_m);
+            ln_max = max(ln_max, GeometryToolbox::Norm(nDim, l_n_minus_m));
+          }
+          vortexTiltingMeasure += nodes->GetVortex_Tilting(jPoint);
+        }
+        
+        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 = (kinematicViscosityTurb+kinematicViscosity)/(uijuij*k2*pow(wallDistance, 2.0));
+        const su2double f_d = 1.0-tanh(pow(8.0*r_d,3.0));
+
+        su2double maxDelta = (ln_max/sqrt(3.0)) * f_kh;
+        if (f_d < 0.999){
+          maxDelta = deltaDDES;
+        }
+
+        const su2double distDES = constDES * maxDelta;
+        lengthScale = wallDistance-f_d*max(0.0,(wallDistance-distDES));
+
         break;
       }
     }