Add fields to VOLUME_OUTPUT

- Include time-averaged skin friction coefficient and instantaneous energy backscatter ratio into volume output fields.
- Fix boundary conditions for the implicit smoothing of the stochastic source term in Langevin equations.
- Fix computation of the source term in turbulence model equation.
- Fix computation of the subgrid kinetic energy (divide eddy viscosity by flow density).
- Diagonalize Jacobian in turbulence equations for numerical robustness.

Author: Angelo Passariello

Common/include/toolboxes/random_toolbox.hpp

@@ -91,32 +91,26 @@ inline double GetRandomUniform(std::mt19937 gen, double xmin = 0.0, double xmax
  * \return Value of Bessel function.
  */
 inline double GetBesselZero(double x) {
-    double abx = fabs(x);
-
-    if (abx < 3.75) {
-        double t = x / 3.75;
-        double p = 1.0 +
-            t*t*(3.5156229 +
-            t*t*(3.0899424 +
-            t*t*(1.2067492 +
-            t*t*(0.2659732 +
-            t*t*(0.0360768 +
-            t*t*0.0045813)))));
-        return log(p);
-    } else {
-        double t = 3.75 / abx;
-        double poly =
-            0.39894228 +
-            t*(0.01328592 +
-            t*(0.00225319 +
-            t*(-0.00157565 +
-            t*(0.00916281 +
-            t*(-0.02057706 +
-            t*(0.02635537 +
-            t*(-0.01647633 +
-            t*0.00392377)))))));
-        return abx - 0.5*log(abx) + log(poly);
-    }
+  double abx = fabs(x);
+
+  if (abx < 3.75) {
+    double t = x / 3.75;
+    double p =
+        1.0 +
+        t * t *
+            (3.5156229 +
+             t * t * (3.0899424 + t * t * (1.2067492 + t * t * (0.2659732 + t * t * (0.0360768 + t * t * 0.0045813)))));
+    return log(p);
+  } else {
+    double t = 3.75 / abx;
+    double poly =
+        0.39894228 +
+        t * (0.01328592 +
+             t * (0.00225319 +
+                  t * (-0.00157565 +
+                       t * (0.00916281 + t * (-0.02057706 + t * (0.02635537 + t * (-0.01647633 + t * 0.00392377)))))));
+    return abx - 0.5 * log(abx) + log(poly);
+  }
 }
 
 /*!
@@ -127,34 +121,33 @@ inline double GetBesselZero(double x) {
  * \return Value of the integral.
  */
 inline double GetBesselIntegral(double beta_x, double beta_y, double beta_z) {
+  const double A = 1.0 + 2.0 * (beta_x + beta_y + beta_z);
+  const double Bx = 2.0 * beta_x;
+  const double By = 2.0 * beta_y;
+  const double Bz = 2.0 * beta_z;
 
-    const double A  = 1.0 + 2.0*(beta_x + beta_y + beta_z);
-    const double Bx = 2.0*beta_x;
-    const double By = 2.0*beta_y;
-    const double Bz = 2.0*beta_z;
+  const int N = 4000;
+  const double t_max = 40.0;
+  const double dt = t_max / N;
 
-    const int    N      = 4000;
-    const double t_max  = 40.0;
-    const double dt     = t_max / N;
+  double sum = 0.0;
 
-    double sum = 0.0;
+  for (int i = 1; i < N; i++) {
+    double t = i * dt;
 
-    for (int i = 1; i < N; i++) {
-        double t = i * dt;
+    double e = exp(-A * t);
 
-        double e = exp(-A*t);
+    double lx = GetBesselZero(Bx * t);
+    double ly = GetBesselZero(By * t);
+    double lz = GetBesselZero(Bz * t);
 
-        double lx = GetBesselZero(Bx*t);
-        double ly = GetBesselZero(By*t);
-        double lz = GetBesselZero(Bz*t);
+    double lin = log(t) - A * t + lx + ly + lz;
 
-        double lin = log(t) - A*t + lx + ly + lz;
+    double integrand = exp(lin);
+    sum += integrand;
+  }
 
-        double integrand = exp(lin);
-        sum += integrand;
-    }
-
-    return sum * dt;
+  return sum * dt;
 }
 
 /// @}

Common/src/CConfig.cpp

@@ -2946,7 +2946,7 @@ void CConfig::SetConfig_Options() {
   addBoolOption("ENFORCE_LES", enforceLES, false);
 
   /* DESCRIPTION: Specify if the stochastic source term must be included in the turbulence model equation */
-  addBoolOption("STOCH_SOURCE_NU", stochSourceNu, false);
+  addBoolOption("STOCH_SOURCE_NU", stochSourceNu, true);
 
   /* DESCRIPTION: Filter width for LES (if negative, it is computed based on the local cell size) */
   addDoubleOption("LES_FILTER_WIDTH", LES_FilterWidth, -1.0);

SU2_CFD/include/numerics/CNumerics.hpp

@@ -666,23 +666,15 @@ class CNumerics {
 
     /* --- Calculate stochastic tensor --- */
 
-    if (lesSensor > 0.999) {
-      stochReynStress[1][0] = - Cmag * density * turbKE * rndVec[2];
-      stochReynStress[2][0] =   Cmag * density * turbKE * rndVec[1];
-      stochReynStress[2][1] = - Cmag * density * turbKE * rndVec[0];
-      for (size_t iDim = 0; iDim < nDim; iDim++) {
-        for (size_t jDim = 0; jDim <= iDim; jDim++) {
-          if (iDim==jDim) {
-            stochReynStress[iDim][jDim] = 0.0;
-          } else {
-            stochReynStress[jDim][iDim] = - stochReynStress[iDim][jDim];
-          }
-        }
-      }
-    } else {
-      for (size_t iDim = 0; iDim < nDim; iDim++) {
-        for (size_t jDim = 0; jDim < nDim; jDim++) {
+    stochReynStress[1][0] = - std::nearbyint(lesSensor) * Cmag * density * turbKE * rndVec[2];
+    stochReynStress[2][0] =   std::nearbyint(lesSensor) * Cmag * density * turbKE * rndVec[1];
+    stochReynStress[2][1] = - std::nearbyint(lesSensor) * Cmag * density * turbKE * rndVec[0];
+    for (size_t iDim = 0; iDim < nDim; iDim++) {
+      for (size_t jDim = 0; jDim <= iDim; jDim++) {
+        if (iDim==jDim) {
           stochReynStress[iDim][jDim] = 0.0;
+        } else {
+          stochReynStress[jDim][iDim] = - stochReynStress[iDim][jDim];
         }
       }
     }

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

@@ -122,18 +122,14 @@ class CSourceBase_TurbSA : public CNumerics {
                                      const CSAVariables& var, TIME_MARCHING time_marching) {
 
     const su2double& nue = ScalarVar_i[0];
-
     const auto& density = V_i[idx.Density()];
 
-    const su2double nut_small = 1.0e-10;
-
-    su2double Ji_2 = pow(var.Ji,2);
-    su2double Ji_3 = Ji_2 * var.Ji;
+    const su2double nut_small = 1.0e-12;
 
-    const su2double nut = nue * var.fv1;
+    const su2double nut = max(nue * var.fv1, nut_small);
 
-    su2double tLES = ct * pow(lengthScale/DES_const, 2) / max(nut, nut_small);
-    su2double tRANS = pow(lengthScale, 2) / max(nut, nut_small);
+    su2double tLES = ct * pow(lengthScale/DES_const, 2) / nut;
+    su2double tRANS = pow(lengthScale, 2) / nut;
     su2double tLR = tLES / tRANS;
     su2double tTurb = tLES / (lesMode_i + (1.0-lesMode_i)*tLR);
     su2double tRat = timeStep / tTurb;
@@ -146,7 +142,7 @@ class CSourceBase_TurbSA : public CNumerics {
     }
 
     su2double scaleFactor = 0.0;
-    if (lesMode_i > 0.999) scaleFactor = 1.0/tTurb * sqrt(2.0/tRat) * density * corrFac;
+    if (std::nearbyint(lesMode_i) == 1) scaleFactor = 1.0/tTurb * sqrt(2.0/tRat) * density * corrFac;
 
     ResidSB[0] = Residual;
     for (unsigned short iVar = 1; iVar < nVar; iVar++ ) {
@@ -157,22 +153,15 @@ class CSourceBase_TurbSA : public CNumerics {
     for (unsigned short iVar = 0; iVar < nVar; iVar++ ) {
       for (unsigned short jVar = 0; jVar < nVar; jVar++ ) {
         JacobianSB_i[iVar][jVar] = 0.0;
-        if (iVar == jVar) JacobianSB_i[iVar][jVar] = -1.0/tTurb * density * Volume;
       }
     }
-    JacobianSB_i[0][0] = Jacobian_i[0];
-
-    su2double dnut_dnue = var.fv1 + 3.0 * var.cv1_3 * Ji_3 / pow(Ji_3 + var.cv1_3, 2);
-    su2double dtTurb_dnut = - ct * pow(lengthScale,2) / (max(nut, nut_small)*max(nut, nut_small));
 
     for (unsigned short iVar = 1; iVar < nVar; iVar++ ) {
-      JacobianSB_i[iVar][0] = - 1.0/tTurb * scaleFactor * stochSource[iVar-1]
-                              + 1.0/(tTurb*tTurb) * ScalarVar_i[iVar]
-                              + density * corrFac * stochSource[iVar-1] / 
-                                (tTurb * sqrt(2.0*tTurb*timeStep));
-      JacobianSB_i[iVar][0] *= dtTurb_dnut * dnut_dnue * Volume;
+      JacobianSB_i[iVar][iVar] = -1.0/tTurb * density * Volume;
     }
 
+    JacobianSB_i[0][0] = Jacobian_i[0];
+
   }
 
   /*!
@@ -182,23 +171,21 @@ class CSourceBase_TurbSA : public CNumerics {
   inline void AddStochSource(const CSAVariables& var, const MatrixType& velGrad, const su2double Cmag) {
 
     su2double dist2 = dist_i * dist_i;
-    const su2double eps = 1.0e-10;
+    const su2double eps = 1.0e-12;
     su2double xi3 = pow(var.Ji, 3);
 
     su2double factor = dist2 / (2.0 * var.fv1 * ScalarVar_i[0] + eps);
     factor /= (3.0 * xi3 * var.cv1_3 / pow(xi3 + var.cv1_3, 2) + var.fv1 + eps);
 
-    const auto& density = V_i[idx.Density()];
-
     su2double tke = pow(var.fv1*ScalarVar_i[0]/dist_i, 2); 
 
-    su2double R12 =   Cmag * density * tke * ScalarVar_i[2];
-    su2double R13 = - Cmag * density * tke * ScalarVar_i[1];
-    su2double R23 =   Cmag * density * tke * ScalarVar_i[0];
+    su2double R12 = (nDim==3 ?   Cmag * tke * ScalarVar_i[3] : 0.0);
+    su2double R13 =            - Cmag * tke * ScalarVar_i[2];
+    su2double R23 =              Cmag * tke * ScalarVar_i[1];
 
-    su2double RGradU = R12 * (velGrad[0][1] - velGrad[1][0]) +
-                       R13 * (velGrad[0][2] - velGrad[2][0]) +
-                       R23 * (velGrad[1][2] - velGrad[2][1]);
+    su2double RGradU =            R12 * (velGrad[0][1] - velGrad[1][0]) +
+                       (nDim==3 ? R13 * (velGrad[0][2] - velGrad[2][0]) +
+                                  R23 * (velGrad[1][2] - velGrad[2][1]) : 0.0);
 
     su2double source_k = - RGradU;
     su2double source_nu = factor * source_k;
@@ -379,7 +366,7 @@ class CSourceBase_TurbSA : public CNumerics {
       /*--- Compute residual for Langevin equations (Stochastic Backscatter Model). ---*/
 
       if (backscatter) {
-        if (lesMode_i > 0.999 && config->GetStochSourceNu())
+        if (std::nearbyint(lesMode_i) == 1 && config->GetStochSourceNu())
           AddStochSource(var, PrimVar_Grad_i + idx.Velocity(), config->GetSBS_Cmag());
         const su2double DES_const = config->GetConst_DES();
         const su2double ctau = config->GetSBS_Ctau();

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -475,10 +475,11 @@ void CFVMFlowSolverBase<V, R>::Viscous_Residual_impl(unsigned long iEdge, CGeome
       numerics->SetStochVar(turbNodes->GetSolution(iPoint, 1 + iDim),
                             turbNodes->GetSolution(jPoint, 1 + iDim), iDim);
     }
-    su2double eddy_visc_i, eddy_visc_j, DES_length_i,
+    su2double rho, eddy_visc_i, eddy_visc_j, DES_length_i,
               DES_length_j, tke_i, tke_j, lesMode_i, lesMode_j;
-    eddy_visc_i = turbNodes->GetmuT(iPoint);
-    eddy_visc_j = turbNodes->GetmuT(jPoint);
+    rho = nodes->GetDensity(iPoint);
+    eddy_visc_i = turbNodes->GetmuT(iPoint) / rho;
+    eddy_visc_j = turbNodes->GetmuT(jPoint) / rho;
     DES_length_i = turbNodes->GetDES_LengthScale(iPoint);
     DES_length_j = turbNodes->GetDES_LengthScale(jPoint);
     lesMode_i = turbNodes->GetLES_Mode(iPoint);