Merge branch 'develop' into fix_cornernode_comp

su2code · Sep 23, 2024 · 56e1ce5 · 56e1ce5
2 parents bb45348 + d0e98e1
commit 56e1ce5
Show file tree

Hide file tree

Showing 4 changed files with 31 additions and 37 deletions.
diff --git a/SU2_CFD/include/fluid/CFluidScalar.hpp b/SU2_CFD/include/fluid/CFluidScalar.hpp
@@ -41,7 +41,6 @@ class CFluidScalar final : public CFluidModel {
  private:
   const int n_species_mixture;            /*!< \brief Number of species in mixture. */
   su2double Gas_Constant;           /*!< \brief Specific gas constant. */
-  const su2double Gamma;                  /*!< \brief Ratio of specific heats of the gas. */
   const su2double Pressure_Thermodynamic; /*!< \brief Constant pressure thermodynamic. */
   const su2double GasConstant_Ref;        /*!< \brief Gas constant reference needed for Nondimensional problems. */
   const su2double Prandtl_Number;         /*!< \brief Prandlt number.*/
@@ -106,7 +105,7 @@ class CFluidScalar final : public CFluidModel {
   /*!
    * \brief Constructor of the class.
    */
-  CFluidScalar(su2double val_Cp, su2double val_gas_constant, su2double val_operating_pressure, const CConfig* config);
+  CFluidScalar(su2double val_operating_pressure, const CConfig* config);
 
   /*!
    * \brief Set viscosity model.

diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
@@ -48,11 +48,11 @@ struct CSAVariables {
   const su2double cb2_sigma = cb2 / sigma;
   const su2double cw1 = cb1 / k2 + (1 + cb2) / sigma;
   const su2double cr1 = 0.5;
-  const su2double CRot = 1.0;
+  const su2double CRot = 2.0;
   const su2double c2 = 0.7, c3 = 0.9;
 
   /*--- List of auxiliary functions ---*/
-  su2double ft2, d_ft2, r, d_r, g, d_g, glim, fw, d_fw, Ji, d_Ji, S, Shat, d_Shat, fv1, d_fv1, fv2, d_fv2;
+  su2double ft2, d_ft2, r, d_r, g, d_g, glim, fw, d_fw, Ji, d_Ji, Shat, d_Shat, fv1, d_fv1, fv2, d_fv2, Prod;
 
   /*--- List of helpers ---*/
   su2double Omega, dist_i_2, inv_k2_d2, inv_Shat, g_6, norm2_Grad;
@@ -151,20 +151,7 @@ class CSourceBase_TurbSA : public CNumerics {
     Residual = 0.0;
     Jacobian_i[0] = 0.0;
 
-    /*--- Evaluate Omega with a rotational correction term. ---*/
-
-    Omega::get(Vorticity_i, nDim, PrimVar_Grad_i + idx.Velocity(), var);
-
-    /*--- Dacles-Mariani et. al. rotation correction ("-R"). ---*/
-    if (options.rot) {
-      var.Omega += var.CRot * min(0.0, StrainMag_i - var.Omega);
-      /*--- Do not allow negative production for SA-neg. ---*/
-      if (ScalarVar_i[0] < 0) var.Omega = abs(var.Omega);
-    }
-
     if (dist_i > 1e-10) {
-      /*--- Vorticity ---*/
-      var.S = var.Omega;
 
       var.dist_i_2 = pow(dist_i, 2);
       const su2double nu = laminar_viscosity / density;
@@ -188,12 +175,24 @@ class CSourceBase_TurbSA : public CNumerics {
       var.fv2 = 1 - ScalarVar_i[0] / (nu + ScalarVar_i[0] * var.fv1);
       var.d_fv2 = -(1 / nu - Ji_2 * var.d_fv1) / pow(1 + var.Ji * var.fv1, 2);
 
-      /*--- Compute ft2 term ---*/
-      ft2::get(var);
+      /*--- Evaluate Omega with a rotational correction term. ---*/
+
+      Omega::get(Vorticity_i, nDim, PrimVar_Grad_i + idx.Velocity(), var);
 
       /*--- Compute modified vorticity ---*/
       ModVort::get(ScalarVar_i[0], nu, var);
       var.inv_Shat = 1.0 / var.Shat;
+      var.Prod = var.Shat;
+
+      /*--- Dacles-Mariani et. al. rotation correction ("-R"). ---*/
+      if (options.rot) {
+        var.Prod += var.CRot * min(0.0, StrainMag_i - var.Omega);
+        /*--- Do not allow negative production for SA-neg. ---*/
+        if (ScalarVar_i[0] < 0) var.Prod = abs(var.Prod);
+      }
+
+      /*--- Compute ft2 term ---*/
+      ft2::get(var);
 
       /*--- Compute auxiliary function r ---*/
       rFunc::get(ScalarVar_i[0], var);
@@ -234,7 +233,6 @@ class CSourceBase_TurbSA : public CNumerics {
       } else if (transition_LM){
 
         var.intermittency = intermittency_eff_i;
-        //var.intermittency = 1.0;
         // Is wrong the reference from NASA?
         // Original max(min(gamma, 0.5), 1.0) always gives 1 as result.
         var.interDestrFactor = min(max(intermittency_i, 0.5), 1.0);
@@ -347,12 +345,12 @@ struct Bsl {
 
     /*--- Limiting of \hat{S} based on "Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model"
      * Note 1 option c in https://turbmodels.larc.nasa.gov/spalart.html ---*/
-    if (Sbar >= - c2 * var.S) {
-      var.Shat = var.S + Sbar;
+    if (Sbar >= - c2 * var.Omega) {
+      var.Shat = var.Omega + Sbar;
     } else {
-      const su2double Num = var.S * (c2 * c2 * var.S + c3 * Sbar);
-      const su2double Den = (c3 - 2 * c2) * var.S - Sbar;
-      var.Shat = var.S + Num / Den;
+      const su2double Num = var.Omega * (c2 * c2 * var.Omega + c3 * Sbar);
+      const su2double Den = (c3 - 2 * c2) * var.Omega - Sbar;
+      var.Shat = var.Omega + Num / Den;
     }
     if (var.Shat <= 1e-10) {
       var.Shat = 1e-10;
@@ -366,12 +364,12 @@ struct Bsl {
 /*! \brief Edward. */
 struct Edw {
   static void get(const su2double& nue, const su2double& nu, CSAVariables& var) {
-    var.Shat = max(var.S * ((1.0 / max(var.Ji, 1.0e-16)) + var.fv1), 1.0e-16);
+    var.Shat = max(var.Omega * ((1.0 / max(var.Ji, 1.0e-16)) + var.fv1), 1.0e-16);
     var.Shat = max(var.Shat, 1.0e-10);
     if (var.Shat <= 1.0e-10) {
       var.d_Shat = 0.0;
     } else {
-      var.d_Shat = -var.S * pow(var.Ji, -2) / nu + var.S * var.d_fv1;
+      var.d_Shat = -var.Omega * pow(var.Ji, -2) / nu + var.Omega * var.d_fv1;
     }
   }
 };
@@ -383,7 +381,7 @@ struct Neg {
       // Baseline solution
       Bsl::get(nue, nu, var);
     } else {
-      var.Shat = 1.0e-10;
+      var.Shat = var.Omega;
       var.d_Shat = 0.0;
     }
     /*--- Don't check whether Sbar <>= -cv2*S.
@@ -447,8 +445,8 @@ struct Bsl {
   static void ComputeProduction(const su2double& nue, const CSAVariables& var, su2double& production,
                                 su2double& jacobian) {
     const su2double factor = var.intermittency * var.cb1;
-    production = factor * (1.0 - var.ft2) * var.Shat * nue;
-    jacobian += factor * (-var.Shat * nue * var.d_ft2 + (1.0 - var.ft2) * (nue * var.d_Shat + var.Shat));
+    production = factor * (1.0 - var.ft2) * var.Prod * nue;
+    jacobian += factor * (-var.Prod * nue * var.d_ft2 + (1.0 - var.ft2) * (nue * var.d_Shat + var.Prod));
   }
 
   static void ComputeDestruction(const su2double& nue, const CSAVariables& var, su2double& destruction,
@@ -481,7 +479,7 @@ struct Neg {
 
   static void ComputeProduction(const su2double& nue, const CSAVariables& var, su2double& production,
                                 su2double& jacobian) {
-    const su2double dP_dnu = var.intermittency * var.cb1 * (1.0 - var.ct3) * var.S;
+    const su2double dP_dnu = var.intermittency * var.cb1 * (1.0 - var.ct3) * var.Prod;
     production = dP_dnu * nue;
     jacobian += dP_dnu;
   }

diff --git a/SU2_CFD/src/fluid/CFluidScalar.cpp b/SU2_CFD/src/fluid/CFluidScalar.cpp
@@ -41,12 +41,9 @@
 #include "../../include/fluid/CPolynomialViscosity.hpp"
 #include "../../include/fluid/CSutherland.hpp"
 
-CFluidScalar::CFluidScalar(su2double val_Cp, su2double val_gas_constant, su2double value_pressure_operating,
-                           const CConfig* config)
+CFluidScalar::CFluidScalar(su2double value_pressure_operating, const CConfig* config)
     : CFluidModel(),
       n_species_mixture(config->GetnSpecies() + 1),
-      Gas_Constant(val_gas_constant),
-      Gamma(config->GetGamma()),
       Pressure_Thermodynamic(value_pressure_operating),
       GasConstant_Ref(config->GetGas_Constant_Ref()),
       Prandtl_Number(config->GetPrandtl_Turb()),

diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
@@ -325,7 +325,7 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i
 
       config->SetGas_Constant(UNIVERSAL_GAS_CONSTANT / (config->GetMolecular_Weight() / 1000.0));
       Pressure_Thermodynamic = config->GetPressure_Thermodynamic();
-      auxFluidModel = new CFluidScalar(config->GetSpecific_Heat_Cp(), config->GetGas_Constant(), Pressure_Thermodynamic, config);
+      auxFluidModel = new CFluidScalar(Pressure_Thermodynamic, config);
       auxFluidModel->SetTDState_T(Temperature_FreeStream, config->GetSpecies_Init());
       break;
 
@@ -482,7 +482,7 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i
         break;
 
       case FLUID_MIXTURE:
-        fluidModel = new CFluidScalar(Specific_Heat_CpND, Gas_ConstantND, Pressure_ThermodynamicND, config);
+        fluidModel = new CFluidScalar(Pressure_ThermodynamicND, config);
         fluidModel->SetTDState_T(Temperature_FreeStreamND, config->GetSpecies_Init());
         break;