diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 73266c7dacd..b8d4fad1e05 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1479,7 +1479,8 @@ 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 SGS_Model_Map = { MakePair("NONE", TURB_SGS_MODEL::NONE) @@ -1487,6 +1488,7 @@ static const MapType SGS_Model_Map = { 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) }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 0f7a5c6821f..17140bc88b6 100644 --- a/Common/src/CConfig.cpp +++ b/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); diff --git a/SU2_CFD/include/sgs_model.hpp b/SU2_CFD/include/sgs_model.hpp index 962b81ffc42..16a82cb7eb8 100644 --- a/SU2_CFD/include/sgs_model.hpp +++ b/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" diff --git a/SU2_CFD/include/sgs_model.inl b/SU2_CFD/include/sgs_model.inl index fdcf1796345..9c45af2ee8b 100644 --- a/SU2_CFD/include/sgs_model.inl +++ b/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 = / + * 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 = / --- */ + 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; + + +} \ No newline at end of file diff --git a/SU2_CFD/src/solvers/CFEM_DG_NSSolver.cpp b/SU2_CFD/src/solvers/CFEM_DG_NSSolver.cpp index 6290a9b54f0..0e1324ca631 100644 --- a/SU2_CFD/src/solvers/CFEM_DG_NSSolver.cpp +++ b/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); diff --git a/TestCases/incomp_navierstokes/dynamic_smagorinsky_sc20012/sc20012_les.cfg b/TestCases/incomp_navierstokes/dynamic_smagorinsky_sc20012/sc20012_les.cfg new file mode 100644 index 00000000000..ab95a20108f --- /dev/null +++ b/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 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 279923612d6..b829f29a2f3 100755 --- a/TestCases/serial_regression.py +++ b/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"