Fix rough wall model in SST turbulence model

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -145,6 +145,7 @@ class CFVMFlowSolverBase : public CSolver {
   su2double** HeatFluxTarget = nullptr;    /*!< \brief Heat transfer coefficient for each boundary and vertex. */
   su2double*** CharacPrimVar = nullptr;    /*!< \brief Value of the characteristic variables at each boundary. */
   su2double*** CSkinFriction = nullptr;    /*!< \brief Skin friction coefficient for each boundary and vertex. */
+  su2double** WallShearStress = nullptr;   /*!< \brief Wall Shear Stress for each boundary and vertex. */
   su2double*** HeatConjugateVar = nullptr; /*!< \brief CHT variables for each boundary and vertex. */
   su2double** CPressure = nullptr;         /*!< \brief Pressure coefficient for each boundary and vertex. */
   su2double** CPressureTarget = nullptr;   /*!< \brief Target Pressure coefficient for each boundary and vertex. */
@@ -1553,6 +1554,16 @@ class CFVMFlowSolverBase : public CSolver {
     return CSkinFriction[val_marker][val_dim][val_vertex];
   }
 
+  /*!
+   * \brief Get the wall shear stress.
+   * \param[in] val_marker - Surface marker where the wall shear stress is computed.
+   * \param[in] val_vertex - Vertex of the marker <i>val_marker</i> where the wall shear stress is evaluated.
+   * \return Value of the wall shear stress.
+   */
+  inline su2double GetWallShearStress(unsigned short val_marker, unsigned long val_vertex) const final {
+    return WallShearStress[val_marker][val_vertex];
+  }
+
   /*!
    * \brief Get the skin friction coefficient.
    * \param[in] val_marker - Surface marker where the coefficient is computed.

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -208,6 +208,10 @@ void CFVMFlowSolverBase<V, R>::Allocate(const CConfig& config) {
     }
   }
 
+  /*--- Wall Shear Stress in all the markers ---*/
+
+  Alloc2D(nMarker, nVertex, WallShearStress);
+
   /*--- Store the values of the temperature and the heat flux density at the boundaries,
    used for coupling with a solid donor cell ---*/
   constexpr auto nHeatConjugateVar = 4u;
@@ -473,6 +477,13 @@ CFVMFlowSolverBase<V, R>::~CFVMFlowSolverBase() {
     delete[] CSkinFriction;
   }
 
+  if (WallShearStress != nullptr) {
+    for (iMarker = 0; iMarker < nMarker; iMarker++) {
+      delete[] WallShearStress[iMarker];
+    }
+    delete[] WallShearStress;
+  }
+
   if (HeatConjugateVar != nullptr) {
     for (iMarker = 0; iMarker < nMarker; iMarker++) {
       for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
@@ -2019,7 +2030,7 @@ 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, div_vel, WallDist[3] = {0.0}, Area, WallShearStress, TauNormal, RefTemp, RefVel2 = 0.0,
+  su2double Viscosity = 0.0, div_vel, WallDist[3] = {0.0}, Area, TauNormal, RefTemp, RefVel2 = 0.0,
             RefDensity = 0.0, GradTemperature, Density = 0.0, WallDistMod, FrictionVel, Mach2Vel, Mach_Motion,
             UnitNormal[3] = {0.0}, TauElem[3] = {0.0}, TauTangent[3] = {0.0}, Tau[3][3] = {{0.0}}, Cp,
             thermal_conductivity, thermal_conductivity_tr, thermal_conductivity_ve = 0.0,
@@ -2236,13 +2247,13 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
       TauNormal = 0.0;
       for (iDim = 0; iDim < nDim; iDim++) TauNormal += TauElem[iDim] * UnitNormal[iDim];
 
-      WallShearStress = 0.0;
+      WallShearStress[iMarker][iVertex] = 0.0;
       for (iDim = 0; iDim < nDim; iDim++) {
         TauTangent[iDim] = TauElem[iDim] - TauNormal * UnitNormal[iDim];
         CSkinFriction[iMarker][iDim][iVertex] = TauTangent[iDim] / (0.5 * RefDensity * RefVel2);
-        WallShearStress += TauTangent[iDim] * TauTangent[iDim];
+        WallShearStress[iMarker][iVertex] += TauTangent[iDim] * TauTangent[iDim];
       }
-      WallShearStress = sqrt(WallShearStress);
+      WallShearStress[iMarker][iVertex] = sqrt(WallShearStress[iMarker][iVertex]);
 
       for (iDim = 0; iDim < nDim; iDim++) WallDist[iDim] = (Coord[iDim] - Coord_Normal[iDim]);
       WallDistMod = 0.0;
@@ -2251,7 +2262,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
       /*--- Compute y+ and non-dimensional velocity ---*/
 
-      FrictionVel = sqrt(fabs(WallShearStress) / Density);
+      FrictionVel = sqrt(fabs(WallShearStress[iMarker][iVertex]) / Density);
       YPlus[iMarker][iVertex] = WallDistMod * FrictionVel / (Viscosity / Density);
 
       /*--- Compute total and maximum heat flux on the wall ---*/

SU2_CFD/include/solvers/CSolver.hpp

@@ -3115,6 +3115,15 @@ class CSolver {
                                             unsigned long val_vertex,
                                             unsigned short val_dim) const { return 0; }
 
+  /*!
+   * \brief A virtual member.
+   * \param[in] val_marker - Surface marker where the wall shear stress is computed.
+   * \param[in] val_vertex - Vertex of the marker <i>val_marker</i> where the wall shear stress is evaluated.
+   * \return Value of the wall shear stress.
+   */
+  inline virtual su2double GetWallShearStress(unsigned short val_marker,
+                                              unsigned long val_vertex) const { return 0; }
+
   /*!
    * \brief A virtual member.
    * \param[in] val_marker - Surface marker where the coefficient is computed.

SU2_CFD/src/solvers/CTurbSSTSolver.cpp

@@ -411,12 +411,8 @@ void CTurbSSTSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_cont
         /*--- Set wall values ---*/
         su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint);
         su2double laminar_viscosity = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint);
+        su2double WallShearStress = solver_container[FLOW_SOL]->GetWallShearStress(val_marker, iVertex);
 
-        su2double WallShearStress = 0.0;
-        for (auto iDim = 0; iDim < nDim; iDim++)
-          WallShearStress += pow(solver_container[FLOW_SOL]->GetCSkinFriction(val_marker, iVertex, iDim),2.0);
-
-        WallShearStress = sqrt(WallShearStress);
         /*--- Compute non-dimensional velocity ---*/
         su2double FrictionVel = sqrt(fabs(WallShearStress)/density);