Cleanup wall functions and add user controls

Common/include/CConfig.hpp

@@ -766,6 +766,7 @@ class CConfig {
   Wrt_Projected_Sensitivity, /*!< \brief Write projected sensitivities (dJ/dx) on surfaces to ASCII file. */
   Plot_Section_Forces;       /*!< \brief Write sectional forces for specified markers. */
   unsigned short
+  wallModelMaxIter,     /*!< \brief maximum number of iterations for the Newton method for the wall model */
   Console_Output_Verb,  /*!< \brief Level of verbosity for console output */
   Kind_Average;         /*!< \brief Particular average for the marker analyze. */
   su2double Gamma,      /*!< \brief Ratio of specific heats of the gas. */
@@ -825,6 +826,8 @@ class CConfig {
   Prandtl_Turb,                    /*!< \brief Turbulent Prandtl number for the gas.  */
   wallModelKappa,                  /*!< \brief von Karman constant kappa for turbulence wall modeling */
   wallModelB,                  /*!< \brief constant B for turbulence wall modeling */
+  wallModelRelFac,             /*!< \brief relaxation factor for the Newton method used in the wall model */
+  wallModelMinYplus,          /*!< \brief minimum Y+ value, below which the wall model is not used anymore */
   Length_Ref,                 /*!< \brief Reference length for non-dimensionalization. */
   Pressure_Ref,               /*!< \brief Reference pressure for non-dimensionalization.  */
   Temperature_Ref,            /*!< \brief Reference temperature for non-dimensionalization.*/
@@ -1673,6 +1676,24 @@ class CConfig {
    */
   su2double GetwallModelKappa(void) const { return wallModelKappa; }
 
+  /*!
+   * \brief Get the value of the max. number of Newton iterations for turbulence wall modeling.
+   * \return max number of iterations.
+   */
+  unsigned short GetwallModelMaxIter() const { return wallModelMaxIter; }
+
+  /*!
+   * \brief Get the value of the relaxation factor for turbulence wall modeling.
+   * \return relaxation factor.
+   */
+  su2double GetwallModelRelFac() const { return wallModelRelFac; }
+
+  /*!
+   * \brief Get the value of the minimum Y+ value below which the wall function is deactivated
+   * \return minimum Y+ value.
+   */
+  su2double GetwallModelMinYPlus() const { return wallModelMinYplus; }
+
   /*!
    * \brief Get the value of the von Karman constant kappa for turbulence wall modeling.
    * \return von Karman constant.

Common/src/CConfig.cpp

@@ -1223,10 +1223,16 @@ void CConfig::SetConfig_Options() {
   addDoubleOption("PRANDTL_LAM", Prandtl_Lam, 0.72);
   /*!\brief PRANDTL_TURB \n DESCRIPTION: Turbulent Prandtl number (0.9 (air), only for compressible flows) \n DEFAULT 0.90 \ingroup Config*/
   addDoubleOption("PRANDTL_TURB", Prandtl_Turb, 0.90);
-  /*!\brief WALLMODELKAPPA \n DESCRIPTION: von Karman constant used for the wall model \n DEFAULT 0.41 \ingroup Config*/
-  addDoubleOption("WALLMODELKAPPA", wallModelKappa, 0.41);
-  /*!\brief WALLMODELB \n DESCRIPTION: constant B used for the wall model \n DEFAULT 5.0 \ingroup Config*/
-  addDoubleOption("WALLMODELB", wallModelB, 5.5);
+  /*!\brief WALLMODEL_KAPPA \n DESCRIPTION: von Karman constant used for the wall model \n DEFAULT 0.41 \ingroup Config*/
+  addDoubleOption("WALLMODEL_KAPPA", wallModelKappa, 0.41);
+  /*!\brief WALLMODEL_MAXITER \n DESCRIPTION: Max iterations used for the wall model \n DEFAULT 200 \ingroup Config*/
+  addUnsignedShortOption("WALLMODEL_MAXITER", wallModelMaxIter, 200);
+  /*!\brief WALLMODEL_RELFAC \n DESCRIPTION: Relaxation factor used for the wall model \n DEFAULT 0.5 \ingroup Config*/
+  addDoubleOption("WALLMODEL_RELFAC", wallModelRelFac, 0.5);
+  /*!\brief WALLMODEL_MINYPLUS \n DESCRIPTION: lower limit for Y+ used for the wall model \n DEFAULT 5.0 \ingroup Config*/
+  addDoubleOption("WALLMODEL_MINYPLUS", wallModelMinYplus, 5.0);
+  /*!\brief WALLMODEL_B \n DESCRIPTION: constant B used for the wall model \n DEFAULT 5.0 \ingroup Config*/
+  addDoubleOption("WALLMODEL_B", wallModelB, 5.5);
   /*!\brief BULK_MODULUS \n DESCRIPTION: Value of the Bulk Modulus  \n DEFAULT 1.42E5 \ingroup Config*/
   addDoubleOption("BULK_MODULUS", Bulk_Modulus, 1.42E5);
   /* DESCRIPTION: Epsilon^2 multipier in Beta calculation for incompressible preconditioner.  */

SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp

@@ -151,7 +151,7 @@ FORCEINLINE void addQCR(const MatrixType& grad, MatrixDbl<nDim>& tau) {
 
 /*!
  * \brief Scale the stress tensor according to the target (from a
- *        wall function) magnitute in the tangential direction.
+ *        wall function) magnitude in the tangential direction.
  */
 template<class Container, size_t nDim>
 FORCEINLINE void addTauWall(Int iPoint, Int jPoint,

SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp

@@ -159,6 +159,7 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD {
       addPerturbedRSM(avgV, avgGrad, turb_ke, tau,
                       uq_eigval_comp, uq_permute, uq_delta_b, uq_urlx);
     }
+
     if(wallFun) addTauWall(iPoint, jPoint, solution.GetTauWall(), unitNormal, tau);
 
     Double cond = derived->thermalConductivity(avgV);

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -2393,11 +2393,12 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
   unsigned long iVertex, iPoint, iPointNormal;
   unsigned short iMarker, iMarker_Monitoring, iDim, jDim;
   unsigned short T_INDEX = 0, TVE_INDEX = 0, VEL_INDEX = 0;
-  su2double Viscosity = 0.0, WallDist[3] = {0.0}, Area, TauNormal, dTn, dTven,
+  su2double Viscosity = 0.0, Area, dTn, dTven,
             GradTemperature, Density = 0.0, WallDistMod, FrictionVel,
-            UnitNormal[3] = {0.0}, TauElem[3] = {0.0}, TauTangent[3] = {0.0}, Tau[3][3] = {{0.0}}, Cp,
+            UnitNormal[3] = {0.0}, TauElem[3] = {0.0}, Tau[3][3] = {{0.0}}, Cp,
             thermal_conductivity, MaxNorm = 8.0, Grad_Vel[3][3] = {{0.0}}, Grad_Temp[3] = {0.0}, AxiFactor;
   const su2double *Coord = nullptr, *Coord_Normal = nullptr, *Normal = nullptr;
+  const su2double minYPlus = config->GetwallModelMinYPlus();
 
   string Marker_Tag, Monitoring_Tag;
 
@@ -2527,35 +2528,33 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
       /*--- Compute wall shear stress (using the stress tensor). Compute wall skin friction coefficient, and heat flux
        * on the wall ---*/
 
-      TauNormal = 0.0;
-      for (iDim = 0; iDim < nDim; iDim++) TauNormal += TauElem[iDim] * UnitNormal[iDim];
+      su2double TauTangent[MAXNDIM] = {0.0};
+      GeometryToolbox::TangentProjection(nDim, Tau, UnitNormal, TauTangent);
+
+      WallShearStress[iMarker][iVertex] = GeometryToolbox::Norm(int(MAXNDIM), TauTangent);
+
+      /*--- For wall functions, the wall stresses need to be scaled by the wallfunction stress Tau_Wall---*/
+      if (wallfunctions && (YPlus[iMarker][iVertex] > minYPlus)){
+        const su2double Tau_Wall = nodes->GetTauWall(iPoint); 
+        const su2double scale = Tau_Wall / WallShearStress[iMarker][iVertex];
+        for (iDim = 0; iDim < nDim; iDim++) {
+          TauTangent[iDim] *= scale;
+          TauElem[iDim] *= scale;
+        }
+        WallShearStress[iMarker][iVertex] = Tau_Wall;
+      }
 
-      WallShearStress[iMarker][iVertex] = 0.0;
       for (iDim = 0; iDim < nDim; iDim++) {
-        TauTangent[iDim] = TauElem[iDim] - TauNormal * UnitNormal[iDim];
-        /* --- in case of wall functions, we have computed the skin friction in the turbulence solver --- */
-        /* --- Note that in the wall model, we switch off the computation when the computed y+ < 5    --- */
-        /* --- We put YPlus to 1.0 so we have to compute skinfriction and the actual y+ in that case as well --- */
-        if (!wallfunctions || (wallfunctions && YPlus[iMarker][iVertex] < 5.0))
-          CSkinFriction[iMarker](iVertex,iDim) = TauTangent[iDim] * factorFric;
-
-        WallShearStress[iMarker][iVertex] += TauTangent[iDim] * TauTangent[iDim];
+        CSkinFriction[iMarker](iVertex,iDim) = TauTangent[iDim] * factorFric;
       }
-      WallShearStress[iMarker][iVertex] = sqrt(WallShearStress[iMarker][iVertex]);
 
-      for (iDim = 0; iDim < nDim; iDim++) WallDist[iDim] = (Coord[iDim] - Coord_Normal[iDim]);
-      WallDistMod = 0.0;
-      for (iDim = 0; iDim < nDim; iDim++) WallDistMod += WallDist[iDim] * WallDist[iDim];
-      WallDistMod = sqrt(WallDistMod);
+      WallDistMod = GeometryToolbox::Distance(nDim, Coord, Coord_Normal);
 
-      /*--- Compute y+ and non-dimensional velocity ---*/
+      /*--- Compute non-dimensional velocity and y+ ---*/
 
       FrictionVel = sqrt(fabs(WallShearStress[iMarker][iVertex]) / Density);
 
-      /* --- in case of wall functions, we have computed YPlus in the turbulence class --- */
-      /* --- Note that we do not recompute y+ when y+<5 because y+ can become > 5 again --- */
-      if (!wallfunctions)
-        YPlus[iMarker][iVertex] = WallDistMod * FrictionVel / (Viscosity / Density);
+      YPlus[iMarker][iVertex] = WallDistMod * FrictionVel / (Viscosity / Density);
 
       /*--- Compute total and maximum heat flux on the wall ---*/
 

SU2_CFD/src/solvers/CIncNSSolver.cpp

@@ -653,13 +653,13 @@ void CIncNSSolver::SetTauWall_WF(CGeometry *geometry, CSolver **solver_container
 
   const su2double Gas_Constant = config->GetGas_ConstantND();
   const su2double Cp = (Gamma / Gamma_Minus_One) * Gas_Constant;
-  constexpr unsigned short max_iter =200; /*--- maximum number of iterations for the Newton Solver---*/
-  const su2double tol = 1e-12;            /*--- convergence criterium for the Newton solver, note that 1e-10 is too large ---*/
-  const su2double relax = 0.5;            /*--- relaxation factor for the Newton solver ---*/
+  const su2double tol = 1e-12;  /*!< \brief convergence criterium for the Newton solver, note that 1e-10 is too large */
+  const unsigned short max_iter =config->GetwallModelMaxIter();
+  const su2double relax = config->GetwallModelRelFac(); 
 
-  /*--- Compute the recovery factor ---*/
+  /*--- Compute the recovery factor 
+   * use Molecular (Laminar) Prandtl number (see Nichols & Nelson, nomenclature ) ---*/
 
-  // Molecular (Laminar) Prandtl number (see Nichols & Nelson, nomenclature )
   const su2double Recovery = pow(config->GetPrandtl_Lam(), (1.0/3.0));
 
   /*--- Typical constants from boundary layer theory ---*/
@@ -679,15 +679,6 @@ void CIncNSSolver::SetTauWall_WF(CGeometry *geometry, CSolver **solver_container
 
     if (config->GetWallFunction_Treatment(Marker_Tag) != WALL_FUNCTIONS::STANDARD_FUNCTION) continue;
 
-    /*--- Get the specified wall heat flux from config ---*/
-    // note that we can get the heat flux from the temperature gradient
-    su2double q_w = 0.0;
-
-    if (config->GetMarker_All_KindBC(iMarker) == HEAT_FLUX)
-      q_w = config->GetWall_HeatFlux(Marker_Tag);
-
-    // heat flux from temperature: q_w = h*(T_wall - T_fluid)
-
     /*--- Loop over all of the vertices on this boundary marker ---*/
 
     SU2_OMP_FOR_DYN(OMP_MIN_SIZE)
@@ -744,13 +735,22 @@ void CIncNSSolver::SetTauWall_WF(CGeometry *geometry, CSolver **solver_container
        * In incompressible flows, we can assume that there is no velocity-related
        * temperature change Prandtl: T+ = Pr*y+ ---*/
       su2double T_Wall = nodes->GetTemperature(iPoint);
+      su2double Conductivity_Wall = nodes->GetThermalConductivity(iPoint);
+
+      /*--- if a wall temperature was given, we compute the local heat flux using k*dT/dn ---*/
 
-      /*--- Extrapolate the pressure from the interior & compute the
-      wall density using the equation of state ---*/
+      su2double q_w = 0.0;
 
+      if (config->GetMarker_All_KindBC(iMarker) == ISOTHERMAL) {
+        const su2double T_n = nodes->GetTemperature(Point_Normal);
+        q_w = Conductivity_Wall * (T_Wall - T_n) / WallDistMod;
+      }  
+      else if (config->GetMarker_All_KindBC(iMarker) == HEAT_FLUX) {
+        q_w = config->GetWall_HeatFlux(Marker_Tag);  
+      }
+
       /*--- incompressible formulation ---*/
       su2double Density_Wall = nodes->GetDensity(iPoint);
-      su2double Conductivity_Wall = nodes->GetThermalConductivity(iPoint);
       su2double Lam_Visc_Normal = nodes->GetLaminarViscosity(Point_Normal);
 
       /*--- Compute the shear stress at the wall in the regular fashion
@@ -759,8 +759,7 @@ void CIncNSSolver::SetTauWall_WF(CGeometry *geometry, CSolver **solver_container
       su2double tau[MAXNDIM][MAXNDIM] = {{0.0}};
       su2double Lam_Visc_Wall = nodes->GetLaminarViscosity(iPoint);
       su2double Eddy_Visc_Wall = nodes->GetEddyViscosity(iPoint);
-      // do we need the total viscosity for the stress tensor?
-      //su2double total_viscosity = (Lam_Visc_Wall + Eddy_Visc_Wall);
+
       CNumerics::ComputeStressTensor(nDim, tau, nodes->GetGradient_Primitive(iPoint)+1, Lam_Visc_Wall);
 
       su2double TauTangent[MAXNDIM] = {0.0};
@@ -775,72 +774,68 @@ void CIncNSSolver::SetTauWall_WF(CGeometry *geometry, CSolver **solver_container
       unsigned long counter = 0;
       su2double diff = 1.0;
       su2double U_Tau = max(1.0e-6,sqrt(WallShearStress/Density_Wall));
-      su2double Y_Plus = 5.0; // clipping value
+      su2double Y_Plus = 0.99*config->GetwallModelMinYPlus(); // clipping value
 
       su2double Y_Plus_Start = Density_Wall * U_Tau * WallDistMod / Lam_Visc_Wall;
 
-      /*--- Automatic switch off when y+ < 5 according to Nichols & Nelson (2004) ---*/
+      /*--- Automatic switch off when y+ < "limit" according to Nichols & Nelson (2004) ---*/
 
-      if (Y_Plus_Start < Y_Plus) {
-        /*--- impose a minimum y+ for stability reasons---*/
+      if (Y_Plus_Start < config->GetwallModelMinYPlus()) {
         smallYPlusCounter++;
       }
-      else {
-        while (fabs(diff) > tol) {
+      else while (fabs(diff) > tol) {
 
-          /*--- Friction velocity and u+ ---*/
+        /*--- Friction velocity and u+ ---*/
 
-          su2double U_Plus = VelTangMod/U_Tau;
+        su2double U_Plus = VelTangMod/U_Tau;
 
-          /*--- Gamma, Beta, Q, and Phi, defined by Nichols & Nelson (2004) page 1110 ---*/
+        /*--- Gamma, Beta, Q, and Phi, defined by Nichols & Nelson (2004) page 1110 ---*/
 
-          su2double Gam  = Recovery*U_Tau*U_Tau/(2.0*Cp*T_Wall);
-          /*--- nijso: heated wall needs validation testcase! ---*/
-          su2double Beta = q_w*Lam_Visc_Wall/(Density_Wall*T_Wall*Conductivity_Wall*U_Tau); // TODO: nonzero heatflux needs validation case
-          su2double Q    = sqrt(Beta*Beta + 4.0*Gam);
-          su2double Phi  = asin(-1.0*Beta/Q);
+        su2double Gam  = Recovery*U_Tau*U_Tau/(2.0*Cp*T_Wall);
+        su2double Beta = q_w*Lam_Visc_Wall/(Density_Wall*T_Wall*Conductivity_Wall*U_Tau); 
+        su2double Q    = sqrt(Beta*Beta + 4.0*Gam);
+        su2double Phi  = asin(-1.0*Beta/Q);
 
-          /*--- Y+ defined by White & Christoph (compressibility and heat transfer) negative value for (2.0*Gam*U_Plus - Beta)/Q ---*/
+        /*--- Y+ defined by White & Christoph (compressibility and heat transfer) negative value for (2.0*Gam*U_Plus - Beta)/Q ---*/
 
-          su2double Y_Plus_White = exp((kappa/sqrt(Gam))*(asin((2.0*Gam*U_Plus - Beta)/Q) - Phi))*exp(-1.0*kappa*B);
+        su2double Y_Plus_White = exp((kappa/sqrt(Gam))*(asin((2.0*Gam*U_Plus - Beta)/Q) - Phi))*exp(-1.0*kappa*B);
 
-          /*--- Spalding's universal form for the BL velocity with the
-           outer velocity form of White & Christoph above. ---*/
+        /*--- Spalding's universal form for the BL velocity with the
+         outer velocity form of White & Christoph above. ---*/
 
-          su2double kUp = kappa*U_Plus;
-          Y_Plus = U_Plus + Y_Plus_White - (exp(-1.0*kappa*B)* (1.0 + kUp + 0.5*kUp*kUp + kUp*kUp*kUp/6.0));
+        const su2double kUp = kappa*U_Plus;
+        Y_Plus = U_Plus + Y_Plus_White - (exp(-1.0*kappa*B)* (1.0 + kUp + 0.5*kUp*kUp + kUp*kUp*kUp/6.0));
 
-          su2double dypw_dyp = 2.0*Y_Plus_White*(kappa*sqrt(Gam)/Q)*sqrt(1.0 - pow(2.0*Gam*U_Plus - Beta,2.0)/(Q*Q));
+        su2double dypw_dyp = 2.0*Y_Plus_White*(kappa*sqrt(Gam)/Q)*sqrt(1.0 - pow(2.0*Gam*U_Plus - Beta,2.0)/(Q*Q));
 
-          Eddy_Visc_Wall = Lam_Visc_Wall*(1.0 + dypw_dyp - kappa*exp(-1.0*kappa*B)*
-                                           (1.0 + kappa*U_Plus + kappa*kappa*U_Plus*U_Plus/2.0)
-                                           - Lam_Visc_Normal/Lam_Visc_Wall);
-          Eddy_Visc_Wall = max(1.0e-6, Eddy_Visc_Wall);
+        Eddy_Visc_Wall = Lam_Visc_Wall*(1.0 + dypw_dyp - kappa*exp(-1.0*kappa*B)*
+                                         (1.0 + kappa*U_Plus + kappa*kappa*U_Plus*U_Plus/2.0)
+                                         - Lam_Visc_Normal/Lam_Visc_Wall);
+        Eddy_Visc_Wall = max(1.0e-6, Eddy_Visc_Wall);
 
-          /* --- Define function for Newton method to zero --- */
+        /* --- Define function for Newton method to zero --- */
 
-          diff = (Density_Wall * U_Tau * WallDistMod / Lam_Visc_Wall) - Y_Plus;
+        diff = (Density_Wall * U_Tau * WallDistMod / Lam_Visc_Wall) - Y_Plus;
 
-          /* --- Gradient of function defined above --- */
+        /* --- Gradient of function defined above --- */
 
-          su2double grad_diff = Density_Wall * WallDistMod / Lam_Visc_Wall + VelTangMod / (U_Tau * U_Tau) +
-                    kappa /(U_Tau * sqrt(Gam)) * asin(U_Plus * sqrt(Gam)) * Y_Plus_White -
-                    exp(-1.0 * B * kappa) * (0.5 * pow(VelTangMod * kappa / U_Tau, 3) +
-                    pow(VelTangMod * kappa / U_Tau, 2) + VelTangMod * kappa / U_Tau) / U_Tau;
+        su2double grad_diff = Density_Wall * WallDistMod / Lam_Visc_Wall + VelTangMod / (U_Tau * U_Tau) +
+                  kappa /(U_Tau * sqrt(Gam)) * asin(U_Plus * sqrt(Gam)) * Y_Plus_White -
+                  exp(-1.0 * B * kappa) * (0.5 * pow(VelTangMod * kappa / U_Tau, 3) +
+                  pow(VelTangMod * kappa / U_Tau, 2) + VelTangMod * kappa / U_Tau) / U_Tau;
 
-          /* --- Newton Step --- */
+        /* --- Newton Step --- */
 
-          U_Tau = U_Tau - relax*(diff / grad_diff);
+        U_Tau = U_Tau - relax*(diff / grad_diff);
 
-          counter++;
-          if (counter > max_iter) {
-            notConvergedCounter++;
-            // use some safe values for convergence
-            Y_Plus = 30.0;
-            Eddy_Visc_Wall = 1.0;
-            U_Tau = 1.0;
-            break;
-          }
+        counter++;
+        if (counter > max_iter) {
+          notConvergedCounter++;
+          // use some safe values for convergence
+          Y_Plus = 30.0;
+          Eddy_Visc_Wall = 1.0;
+          U_Tau = 1.0;
+          break;
         }
       }
 
@@ -852,21 +847,9 @@ void CIncNSSolver::SetTauWall_WF(CGeometry *geometry, CSolver **solver_container
       EddyViscWall[iMarker][iVertex] = Eddy_Visc_Wall;
       UTau[iMarker][iVertex] = U_Tau;
 
-      // wall model value
       su2double Tau_Wall = (1.0/Density_Wall)*pow(Y_Plus*Lam_Visc_Wall/WallDistMod,2.0);
 
-      // nijso: skinfriction for wall functions gives opposite sign?
-
-      for (auto iDim = 0u; iDim < nDim; iDim++)
-        CSkinFriction[iMarker](iVertex,iDim) = (Tau_Wall/WallShearStress)*TauTangent[iDim] / DynamicPressureRef;
-
       nodes->SetTauWall(iPoint, Tau_Wall);
-      // for compressible flow:
-      //nodes->SetTemperature(iPoint,T_Wall);
-      //nodes->SetSolution(iPoint, 0, Density_Wall);
-      //nodes->SetPrimitive(iPoint, nDim + 1, P_Wall);
-      // for incompressible flow:
-      // ...?
 
     }
     END_SU2_OMP_FOR