Merge pull request #1360 from bigfooted/fix_wallfunctions

Cleanup wall functions and add user controls

Author: TobiKattmann

Common/include/CConfig.hpp

@@ -821,10 +821,13 @@ class CConfig {
   Pressure_Thermodynamic,          /*!< \brief Thermodynamic pressure of the fluid. */
   Temperature_FreeStream,          /*!< \brief Total temperature of the fluid.  */
   Temperature_ve_FreeStream;       /*!< \brief Total vibrational-electronic temperature of the fluid.  */
-  su2double Prandtl_Lam,           /*!< \brief Laminar Prandtl number for the gas.  */
-  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 */
+  unsigned short wallModel_MaxIter; /*!< \brief maximum number of iterations for the Newton method for the wall model */
+  su2double wallModel_Kappa,        /*!< \brief von Karman constant kappa for turbulence wall modeling */
+  wallModel_B,                      /*!< \brief constant B for turbulence wall modeling */
+  wallModel_RelFac,                 /*!< \brief relaxation factor for the Newton method used in the wall model */
+  wallModel_MinYplus;               /*!< \brief minimum Y+ value, below which the wall model is not used anymore */
+  su2double Prandtl_Lam,      /*!< \brief Laminar Prandtl number for the gas.  */
+  Prandtl_Turb,               /*!< \brief Turbulent Prandtl number for the gas.  */
   Length_Ref,                 /*!< \brief Reference length for non-dimensionalization. */
   Pressure_Ref,               /*!< \brief Reference pressure for non-dimensionalization.  */
   Temperature_Ref,            /*!< \brief Reference temperature for non-dimensionalization.*/
@@ -833,7 +836,7 @@ class CConfig {
   Velocity_Ref,               /*!< \brief Reference velocity for non-dimensionalization.*/
   Time_Ref,                   /*!< \brief Reference time for non-dimensionalization. */
   Viscosity_Ref,              /*!< \brief Reference viscosity for non-dimensionalization. */
-  Conductivity_Ref,           /*!< \brief Reference conductivity for non-dimensionalization. */
+  Thermal_Conductivity_Ref,   /*!< \brief Reference conductivity for non-dimensionalization. */
   Energy_Ref,                 /*!< \brief Reference viscosity for non-dimensionalization. */
   Wall_Temperature,           /*!< \brief Temperature at an isotropic wall in Kelvin. */
   Omega_Ref,                  /*!< \brief Reference angular velocity for non-dimensionalization. */
@@ -1671,13 +1674,31 @@ class CConfig {
    * \brief Get the value of the von Karman constant kappa for turbulence wall modeling.
    * \return von Karman constant.
    */
-  su2double GetwallModelKappa(void) const { return wallModelKappa; }
+  su2double GetwallModel_Kappa() const { return wallModel_Kappa; }
 
   /*!
-   * \brief Get the value of the von Karman constant kappa for turbulence wall modeling.
-   * \return von Karman constant.
+   * \brief Get the value of the max. number of Newton iterations for turbulence wall modeling.
+   * \return Max number of iterations.
+   */
+  unsigned short GetwallModel_MaxIter() const { return wallModel_MaxIter; }
+
+  /*!
+   * \brief Get the value of the relaxation factor for turbulence wall modeling.
+   * \return Relaxation factor.
+   */
+  su2double GetwallModel_RelFac() const { return wallModel_RelFac; }
+
+  /*!
+   * \brief Get the value of the minimum Y+ value below which the wall function is deactivated.
+   * \return Minimum Y+ value.
+   */
+  su2double GetwallModel_MinYPlus() const { return wallModel_MinYplus; }
+
+  /*!
+   * \brief Get the value of the wall model constant B for turbulence wall modeling.
+   * \return Wall model constant B.
    */
-  su2double GetwallModelB(void) const { return wallModelB; }
+  su2double GetwallModel_B() const { return wallModel_B; }
 
   /*!
    * \brief Get the value of the thermal diffusivity for solids.
@@ -1753,10 +1774,10 @@ class CConfig {
   su2double GetFan_Poly_Eff(void) const { return Fan_Poly_Eff; }
 
   /*!
-   * \brief Get the value of the reference conductivity for non-dimensionalization.
-   * \return Reference conductivity for non-dimensionalization.
+   * \brief Get the value of the reference thermal conductivity for non-dimensionalization.
+   * \return Reference thermal conductivity for non-dimensionalization.
    */
-  su2double GetConductivity_Ref(void) const { return Conductivity_Ref; }
+  su2double GetThermal_Conductivity_Ref(void) const { return Thermal_Conductivity_Ref; }
 
   /*!
    * \brief Get the value of the reference angular velocity for non-dimensionalization.
@@ -2376,7 +2397,7 @@ class CConfig {
    * \brief Set the reference conductivity for nondimensionalization.
    * \param[in] val_conductivity_ref - Value of the reference conductivity.
    */
-  void SetConductivity_Ref(su2double val_conductivity_ref) { Conductivity_Ref = val_conductivity_ref; }
+  void SetConductivity_Ref(su2double val_conductivity_ref) { Thermal_Conductivity_Ref = val_conductivity_ref; }
 
   /*!
    * \brief Set the nondimensionalized freestream pressure.

Common/src/CConfig.cpp

@@ -1223,10 +1223,20 @@ 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);
+
+  /*--- Options related to wall models. ---*/
+
+  /*!\brief WALLMODEL_KAPPA \n DESCRIPTION: von Karman constant used for the wall model \n DEFAULT 0.41 \ingroup Config*/
+  addDoubleOption("WALLMODEL_KAPPA", wallModel_Kappa, 0.41);
+  /*!\brief WALLMODEL_MAXITER \n DESCRIPTION: Max iterations used for the wall model \n DEFAULT 200 \ingroup Config*/
+  addUnsignedShortOption("WALLMODEL_MAXITER", wallModel_MaxIter, 200);
+  /*!\brief WALLMODEL_RELFAC \n DESCRIPTION: Relaxation factor used for the wall model \n DEFAULT 0.5 \ingroup Config*/
+  addDoubleOption("WALLMODEL_RELFAC", wallModel_RelFac, 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", wallModel_MinYplus, 5.0);
+  /*!\brief WALLMODEL_B \n DESCRIPTION: constant B used for the wall model \n DEFAULT 5.5 \ingroup Config*/
+  addDoubleOption("WALLMODEL_B", wallModel_B, 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.  */
@@ -3295,11 +3305,18 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
       if (Kind_WallFunctions[iMarker] != WALL_FUNCTIONS::NONE)
         Wall_Functions = true;
 
-      if ((Kind_WallFunctions[iMarker] == WALL_FUNCTIONS::ADAPTIVE_FUNCTION) || (Kind_WallFunctions[iMarker] == WALL_FUNCTIONS::SCALABLE_FUNCTION)
-        || (Kind_WallFunctions[iMarker] == WALL_FUNCTIONS::NONEQUILIBRIUM_MODEL))
-
+      if ((Kind_WallFunctions[iMarker] == WALL_FUNCTIONS::ADAPTIVE_FUNCTION) ||
+          (Kind_WallFunctions[iMarker] == WALL_FUNCTIONS::SCALABLE_FUNCTION) ||
+          (Kind_WallFunctions[iMarker] == WALL_FUNCTIONS::NONEQUILIBRIUM_MODEL))
         SU2_MPI::Error(string("For RANS problems, use NONE, STANDARD_WALL_FUNCTION or EQUILIBRIUM_WALL_MODEL.\n"), CURRENT_FUNCTION);
 
+      if (Kind_WallFunctions[iMarker] == WALL_FUNCTIONS::STANDARD_FUNCTION) {
+        if (!((Kind_Solver == RANS) || (Kind_Solver == INC_RANS)))
+          SU2_MPI::Error(string("Wall model STANDARD_FUNCTION only available for RANS or INC_RANS.\n"), CURRENT_FUNCTION);
+        if (nRough_Wall != 0)
+          SU2_MPI::Error(string("Wall model STANDARD_FUNCTION and WALL_ROUGHNESS migh not be compatible. Checking required!\n"), CURRENT_FUNCTION);
+      }
+
     }
   }
 

Common/src/wall_model.cpp

@@ -41,7 +41,7 @@ CWallModel::CWallModel(CConfig *config) {
   Pr_lam  = config->GetPrandtl_Lam();
   Pr_turb = config->GetPrandtl_Turb();
 
-  karman = config->GetwallModelKappa(); // von Karman constant -> k = 0.41; or 0.38;
+  karman = config->GetwallModel_Kappa(); // von Karman constant -> k = 0.41; or 0.38;
 }
 
 void CWallModel::WallShearStressAndHeatFlux(const su2double rhoExchange,

SU2_CFD/include/numerics/CNumerics.hpp

@@ -1554,7 +1554,7 @@ class CNumerics {
    * \param[in] val_tauwall_i - Tauwall at point i
    * \param[in] val_tauwall_j - Tauwall at point j
    */
-  inline virtual void SetTauWall(su2double val_tauwall_i, su2double val_tauwall_j) { }
+  inline virtual void SetTau_Wall(su2double val_tauwall_i, su2double val_tauwall_j) { }
 
   /*!
    * \brief - Calculate the central/upwind blending function for a face

SU2_CFD/include/numerics/flow/flow_diffusion.hpp

@@ -186,7 +186,7 @@ class CAvgGrad_Base : public CNumerics {
    * \param[in] val_tauwall_i - Value of the wall shear stress at point i.
    * \param[in] val_tauwall_j - Value of the wall shear stress at point j.
    */
-  inline void SetTauWall(su2double val_tauwall_i, su2double val_tauwall_j) override {
+  inline void SetTau_Wall(su2double val_tauwall_i, su2double val_tauwall_j) override {
     TauWall_i = val_tauwall_i;
     TauWall_j = val_tauwall_j;
   }

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,7 +159,8 @@ 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);
+
+    if(wallFun) addTauWall(iPoint, jPoint, solution.GetTau_Wall(), unitNormal, tau);
 
     Double cond = derived->thermalConductivity(avgV);
     VectorDbl<nDim> heatFlux;

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -443,8 +443,8 @@ void CFVMFlowSolverBase<V, R>::Viscous_Residual_impl(unsigned long iEdge, CGeome
 
   /*--- Wall shear stress values (wall functions) ---*/
 
-  numerics->SetTauWall(nodes->GetTauWall(iPoint),
-                       nodes->GetTauWall(jPoint));
+  numerics->SetTau_Wall(nodes->GetTau_Wall(iPoint),
+                        nodes->GetTau_Wall(jPoint));
 
   /*--- Compute and update residual ---*/
 
@@ -2393,10 +2393,11 @@ 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, Density = 0.0, GradTemperature = 0.0,
-            UnitNormal[3] = {0.0}, TauElem[3] = {0.0}, TauTangent[3] = {0.0}, Tau[3][3] = {{0.0}}, Cp,
+  su2double Viscosity = 0.0, Area, Density = 0.0, GradTemperature = 0.0, WallDistMod, FrictionVel,
+            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->GetwallModel_MinYPlus();
 
   string Marker_Tag, Monitoring_Tag;
 
@@ -2526,30 +2527,37 @@ 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 ---*/
 
-      su2double TauNormal = GeometryToolbox::DotProduct(nDim, TauElem, UnitNormal);
+      su2double TauTangent[MAXNDIM] = {0.0};
+      GeometryToolbox::TangentProjection(nDim, Tau, UnitNormal, TauTangent);
 
-      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] = GeometryToolbox::Norm(int(MAXNDIM), TauTangent);
+
+      /*--- For wall functions, the wall stresses need to be scaled by the wallfunction stress Tau_Wall---*/
+      su2double Tau_Wall, scale;
+      if (wallfunctions && (YPlus[iMarker][iVertex] > minYPlus)){
+        Tau_Wall = nodes->GetTau_Wall(iPoint);
+        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] = GeometryToolbox::Norm(nDim, TauTangent);
 
-      for (iDim = 0; iDim < nDim; iDim++) WallDist[iDim] = (Coord[iDim] - Coord_Normal[iDim]);
-      
-      su2double WallDistMod = GeometryToolbox::Norm(nDim, WallDist);
+      for (iDim = 0; iDim < nDim; iDim++) {
+        CSkinFriction[iMarker](iVertex,iDim) = TauTangent[iDim] * factorFric;
+      }
 
-      /*--- Compute y+ and non-dimensional velocity ---*/
+      WallDistMod = GeometryToolbox::Distance(nDim, Coord, Coord_Normal);
 
-      su2double FrictionVel = sqrt(fabs(WallShearStress[iMarker][iVertex]) / Density);
+      /*--- Compute non-dimensional velocity and y+ ---*/
 
-      /* --- 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)
+      FrictionVel = sqrt(fabs(WallShearStress[iMarker][iVertex]) / Density);
+      
+      if (!wallfunctions) {
         YPlus[iMarker][iVertex] = WallDistMod * FrictionVel / (Viscosity / Density);
+      }
 
       /*--- Compute total and maximum heat flux on the wall ---*/
 

SU2_CFD/include/solvers/CIncNSSolver.hpp

@@ -71,7 +71,7 @@ class CIncNSSolver final : public CIncEulerSolver {
    * \param[in] solver_container - Container vector with all the solutions.
    * \param[in] config - Definition of the particular problem.
    */
-  void SetTauWall_WF(CGeometry *geometry,
+  void SetTau_Wall_WF(CGeometry *geometry,
                      CSolver** solver_container,
                      const CConfig* config);
 

SU2_CFD/include/solvers/CNEMONSSolver.hpp

@@ -70,7 +70,7 @@ class CNEMONSSolver final : public CNEMOEulerSolver {
    * \param[in] solver_container - Container vector with all the solutions.
    * \param[in] config - Definition of the particular problem.
    */
-  void SetTauWall_WF(CGeometry *geometry, CSolver** solver_container, const CConfig* config);
+  void SetTau_Wall_WF(CGeometry *geometry, CSolver** solver_container, const CConfig* config);
 
 public: