feat:Implement Dynamic Smagorinsky SGS model for LES (Germano-Lilly, local approximation)

Common/include/option_structure.hpp

@@ -1479,14 +1479,16 @@ enum class TURB_SGS_MODEL {
   IMPLICIT_LES, /*!< \brief Implicit LES, i.e. no explicit SGS model. */
   SMAGORINSKY , /*!< \brief Smagorinsky SGS model. */
   WALE        , /*!< \brief Wall-Adapting Local Eddy-viscosity SGS model. */
-  VREMAN        /*!< \brief Vreman SGS model. */
+  VREMAN,        /*!< \brief Vreman SGS model. */
+  DYNAMIC_SMAGORINSKY /*!< \brief Dynamic Smagorinsky SGS model. */
 };
 static const MapType<std::string, TURB_SGS_MODEL> SGS_Model_Map = {
   MakePair("NONE",         TURB_SGS_MODEL::NONE)
   MakePair("IMPLICIT_LES", TURB_SGS_MODEL::IMPLICIT_LES)
   MakePair("SMAGORINSKY",  TURB_SGS_MODEL::SMAGORINSKY)
   MakePair("WALE",         TURB_SGS_MODEL::WALE)
   MakePair("VREMAN",       TURB_SGS_MODEL::VREMAN)
+  MakePair("DYNAMIC_SMAGORINSKY", TURB_SGS_MODEL::DYNAMIC_SMAGORINSKY)
 };
 
 /*!

Common/src/CConfig.cpp

@@ -6576,6 +6576,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
           case TURB_SGS_MODEL::SMAGORINSKY:  cout << "Smagorinsky " << endl; break;
           case TURB_SGS_MODEL::WALE:         cout << "WALE"         << endl; break;
           case TURB_SGS_MODEL::VREMAN:       cout << "VREMAN"         << endl; break;
+          case TURB_SGS_MODEL::DYNAMIC_SMAGORINSKY: cout << "Dynamic Smagorinsky" << endl; break;
           default:
             SU2_MPI::Error("Subgrid Scale model not specified.", CURRENT_FUNCTION);
 

SU2_CFD/include/sgs_model.hpp

@@ -784,4 +784,47 @@ class CVremanModel : public CSGSModel {
                                    su2double &dMuTdy,
                                    su2double &dMuTdz) override;
 };
+
+class CDynamicSmagorinskyModel : public CSGSModel {
+  public:
+    su2double filter_ratio; /*!< \brief Multiplier to get filter width from grid length scale. */
+    su2double Cs2_clip_max; /*!< \brief Maximum value of the clipped Cs^2. */
+
+    /*!
+     * \brief Constructor of the class.
+     */
+    CDynamicSmagorinskyModel(void);
+
+    /*!
+     * \brief Destructor of the class.
+     */
+    ~CDynamicSmagorinskyModel(void) override;
+
+    /*!
+     * \brief Compute eddy viscosity for 2D (not implemented).
+     */
+    su2double ComputeEddyViscosity_2D(const su2double rho,
+                                      const su2double dudx,
+                                      const su2double dudy,
+                                      const su2double dvdx,
+                                      const su2double dvdy,
+                                      const su2double lenScale,
+                                      const su2double distToWall) override;
+
+    /*!
+     * \brief Compute eddy viscosity for 3D using Germano-Lilly dynamic procedure.
+     */
+    su2double ComputeEddyViscosity_3D(const su2double rho,
+                                      const su2double dudx,
+                                      const su2double dudy,
+                                      const su2double dudz,
+                                      const su2double dvdx,
+                                      const su2double dvdy,
+                                      const su2double dvdz,
+                                      const su2double dwdx,
+                                      const su2double dwdy,
+                                      const su2double dwdz,
+                                      const su2double lenScale,
+                                      const su2double distToWall) override;
+};
 #include "sgs_model.inl"

SU2_CFD/include/sgs_model.inl

@@ -532,3 +532,94 @@ inline void CVremanModel::ComputeGradEddyViscosity_3D(const su2double rho,
                                                     su2double &dMuTdz) {
   SU2_MPI::Error("Not implemented yet", CURRENT_FUNCTION);
 }
+
+inline CDynamicSmagorinskyModel::CDynamicSmagorinskyModel(void) : CSGSModel() {
+  filter_ratio = 2.0;
+  Cs2_clip_max = 0.04;
+}
+
+inline CDynamicSmagorinskyModel::~CDynamicSmagorinskyModel(void){}
+
+inline su2double CDynamicSmagorinskyModel::ComputeEddyViscosity_2D(const su2double rho,
+                                                    const su2double dudx,
+                                                    const su2double dudy,
+                                                    const su2double dvdx,
+                                                    const su2double dvdy,
+                                                    const su2double lenScale,
+                                                    const su2double distToWall) {
+  SU2_MPI::Error("Not implemented yet", CURRENT_FUNCTION);
+  return 0.0;
+}
+
+inline su2double CDynamicSmagorinskyModel::ComputeEddyViscosity_3D(const su2double rho,
+                                                    const su2double dudx,
+                                                    const su2double dudy,
+                                                    const su2double dudz,
+                                                    const su2double dvdx,
+                                                    const su2double dvdy,
+                                                    const su2double dvdz,
+                                                    const su2double dwdx,
+                                                    const su2double dwdy,
+                                                    const su2double dwdz,
+                                                    const su2double lenScale,
+                                                    const su2double distToWall) {
+  /* Germano lilly Dynamic Smagorinsky (Local approximation)
+
+  Grid filter width is calculated as: Delta = lenScale
+  Test filter width is calculated as: L_t = filter_ratio * lenScale
+   * Step 1 — grid-filter strain rate tensor S_ij and its magnitude |S|
+   * Step 2 — test-filter strain rate hatS_ij = S_ij (local approx, no stencil)
+   * Step 3 — Germano identity residual tensor L_ij (Bardina approximation)
+   * Step 4 — M_ij = 2 Delta^2 |S| S_ij  -  2 hatDelta^2 |hatS| hatS_ij
+   * Step 5 — Lilly least-squares:  C_s^2 = <L_ij M_ij> / <M_ij M_ij>
+   * Step 6 — Clip C_s^2 >= 0 (no backscatter), return rho * C_s^2 * Delta^2 * |S|
+
+
+
+  */
+  const su2double S11 = dudx, S22 = dvdy, S33 = dwdz;
+  const su2double S12 = 0.5*(dudy + dvdx);
+  const su2double S13 = 0.5*(dudz + dwdx);
+  const su2double S23 = 0.5*(dvdz + dwdy);
+  // const su2double S21 = S12, S31 = S13, S32 = S23;
+
+  const su2double S_mag = sqrt(2.0*(S11*S11 + S22*S22 + S33*S33
+                          + 2.0*(S12*S12 + S13*S13 + S23*S23)));
+  // Step 2 
+  const su2double hatS11 = S11, hatS22 = S22, hatS33 = S33;
+  const su2double hatS12 = S12, hatS13 = S13, hatS23 = S23;
+  // const su2double hatS21 = hatS12, hatS31 = hatS13, hatS32 = hatS23;
+
+  // Step 3 
+  // Lij = (hatDelta^2 - Delta^2) * 2 * S_ik S_kj 
+  const su2double Delta = lenScale; 
+  const su2double hatDelta = filter_ratio * lenScale;
+  const su2double Delta2 = Delta * Delta;
+  const su2double hatDelta2 = hatDelta * hatDelta;
+  const su2double scaleL = hatDelta2 - Delta2;
+
+   /* L_ij = scaleL * 2 * S_ik S_kj */
+   const su2double L11 = scaleL * 2.0 * (S11*S11 + S12*S12 + S13*S13);
+   const su2double L22 = scaleL * 2.0 * (S12*S12 + S22*S22 + S23*S23);
+   const su2double L33 = scaleL * 2.0 * (S13*S13 + S23*S23 + S33*S33);
+   const su2double L12 = scaleL * 2.0 * (S11*S12 + S12*S22 + S13*S23);
+   const su2double L13 = scaleL * 2.0 * (S11*S13 + S12*S23 + S13*S33);
+   const su2double L23 = scaleL * 2.0 * (S12*S13 + S22*S23 + S23*S33);
+
+   /* --- Step 4: M_ij = 2 Delta^2 |S| S_ij  -  2 hatDelta^2 |hatS| hatS_ij --- */
+   const su2double M11 = 2.0 * Delta2 * S_mag * S11 - 2.0 * hatDelta2 * S_mag * hatS11;
+   const su2double M22 = 2.0 * Delta2 * S_mag * S22 - 2.0 * hatDelta2 * S_mag * hatS22;
+   const su2double M33 = 2.0 * Delta2 * S_mag * S33 - 2.0 * hatDelta2 * S_mag * hatS33;
+   const su2double M12 = 2.0 * Delta2 * S_mag * S12 - 2.0 * hatDelta2 * S_mag * hatS12;
+   const su2double M13 = 2.0 * Delta2 * S_mag * S13 - 2.0 * hatDelta2 * S_mag * hatS13;
+   const su2double M23 = 2.0 * Delta2 * S_mag * S23 - 2.0 * hatDelta2 * S_mag * hatS23;
+
+   /* --- Step 5: Lilly least-squares:  C_s^2 = <L_ij M_ij> / <M_ij M_ij> --- */
+   const su2double LdotM = L11*M11 + L22*M22 + L33*M33 + 2.0*(L12*M12 + L13*M13 + L23*M23);
+   const su2double MdotM = M11*M11 + M22*M22 + M33*M33 + 2.0*(M12*M12 + M13*M13 + M23*M23);
+
+   const su2double Cs2 = LdotM / (MdotM + 1.0E-20);
+   return rho * max(Cs2, 0.0) * Delta2 * S_mag;
+
+
+}

SU2_CFD/src/solvers/CFEM_DG_NSSolver.cpp

@@ -140,6 +140,11 @@ CFEM_DG_NSSolver::CFEM_DG_NSSolver(CGeometry *geometry, CConfig *config, unsigne
         SGSModel     = new CVremanModel;
         SGSModelUsed = true;
         break;
+
+      case TURB_SGS_MODEL::DYNAMIC_SMAGORINSKY:
+        SGSModel     = new CDynamicSmagorinskyModel;
+        SGSModelUsed = true;
+        break;
 
       default:
         SU2_MPI::Error("Unknown SGS model encountered", CURRENT_FUNCTION);

TestCases/incomp_navierstokes/dynamic_smagorinsky_sc20012/sc20012_les.cfg

@@ -0,0 +1,26 @@
+SOLVER= INC_NAVIER_STOKES
+KIND_TURB_MODEL= NONE
+KIND_SGS_MODEL= DYNAMIC_SMAGORINSKY
+FLUID_MODEL= CONSTANT_DENSITY
+INC_DENSITY_MODEL= CONSTANT
+INC_DENSITY_INIT= 1.225
+VISCOSITY_MODEL= CONSTANT_VISCOSITY
+MU_CONSTANT= 1.789e-5
+INC_VELOCITY_INIT= ( 102.0, 0.0, 0.0 )
+MARKER_HEATFLUX= ( airfoil, 0.0 )
+MARKER_FAR= ( farfield )
+MARKER_SYM= ( symmetry, wake_axis )
+MARKER_OUTLET= ( outflow, 0.0, )
+INC_OUTLET_TYPE= PRESSURE_OUTLET
+NUM_METHOD_GRAD= GREEN_GAUSS
+CONV_NUM_METHOD_FLOW= FDS
+TIME_DOMAIN= YES
+TIME_MARCHING= TIME_STEPPING
+TIME_STEP= 1e-5
+MAX_TIME= 0.01
+TIME_ITER= 10
+INNER_ITER= 5
+MESH_FILENAME= sc20012.su2
+MESH_FORMAT= SU2
+SCREEN_OUTPUT= TIME_ITER, RMS_PRESSURE, RMS_VELOCITY-X
+OUTPUT_FILES= NONE

TestCases/serial_regression.py

@@ -811,6 +811,15 @@ def main():
     ddes_flatplate.unsteady  = True
     test_list.append(ddes_flatplate)
 
+    # Dynamic Smagorinsky SGS model - SC20012 airfoil
+    dynsmag_sc20012 = TestCase('dynsmag_sc20012')
+    dynsmag_sc20012.cfg_dir = "incomp_navierstokes/dynamic_smagorinsky_sc20012"
+    dynsmag_sc20012.cfg_file = "sc20012_les.cfg"
+    dynsmag_sc20012.test_iter = 10
+    dynsmag_sc20012.test_vals = [-2.985351, -2.986340]
+    dynsmag_sc20012.unsteady = True
+    test_list.append(dynsmag_sc20012)
+
     # unsteady pitching NACA0015, SA
     unst_inc_turb_naca0015_sa           = TestCase('unst_inc_turb_naca0015_sa')
     unst_inc_turb_naca0015_sa.cfg_dir   = "unsteady/pitching_naca0015_rans_inc"