remove frictional wall heating from incompressible wall function model (#1741)

* remove frictional wall heating from incompressible wall function model
* update regression test for incompressible wall functions

Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com>

Author: bigfooted

SU2_CFD/src/solvers/CIncNSSolver.cpp

@@ -619,15 +619,9 @@ void CIncNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_containe
   unsigned long notConvergedCounter = 0;  /*--- Counts the number of wall cells that are not converged ---*/
   unsigned long smallYPlusCounter = 0;    /*--- Counts the number of wall cells where y+ < 5 ---*/
 
-  const su2double Gas_Constant = config->GetGas_ConstantND();
-  const su2double Cp = (Gamma / Gamma_Minus_One) * Gas_Constant;
   const auto max_iter = config->GetwallModel_MaxIter();
   const su2double relax = config->GetwallModel_RelFac();
 
-  /*--- Compute the recovery factor use 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 ---*/
 
   const su2double kappa = config->GetwallModel_Kappa();
@@ -698,34 +692,7 @@ void CIncNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_containe
 
       su2double WallDistMod = GeometryToolbox::Norm(int(MAXNDIM), WallDist);
 
-      /*--- Compute the wall temperature using the Crocco-Buseman equation.
-       * In incompressible flows, we can assume that there is no velocity-related
-       * temperature change Prandtl: T+ = Pr*y+ ---*/
-      su2double T_Wall = nodes->GetTemperature(iPoint);
-      const su2double Conductivity_Wall = nodes->GetThermalConductivity(iPoint);
-
-      /*--- If a wall temperature was given, we compute the local heat flux using k*dT/dn ---*/
-
-      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) / config->GetHeat_Flux_Ref();
-      }
-      else if (config->GetMarker_All_KindBC(iMarker) == HEAT_TRANSFER) {
-        const su2double Transfer_Coefficient = config->GetWall_HeatTransfer_Coefficient(Marker_Tag) *
-                                               config->GetTemperature_Ref()/config->GetHeat_Flux_Ref();
-        const su2double Tinfinity = config->GetWall_HeatTransfer_Temperature(Marker_Tag) / config->GetTemperature_Ref();
-        q_w = Transfer_Coefficient * (Tinfinity - T_Wall);
-      }
-
-      /*--- Incompressible formulation ---*/
-
       su2double Density_Wall = nodes->GetDensity(iPoint);
-      const su2double Lam_Visc_Normal = nodes->GetLaminarViscosity(Point_Normal);
 
       /*--- Compute the shear stress at the wall in the regular fashion
        *    by using the stress tensor on the surface ---*/
@@ -766,41 +733,33 @@ void CIncNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_containe
 
         /*--- Friction velocity and u+ ---*/
 
-        const su2double U_Plus = VelTangMod/U_Tau;
-
-        /*--- Gamma, Beta, Q, and Phi, defined by Nichols & Nelson (2004) page 1110 ---*/
-
-        const su2double Gam  = Recovery*U_Tau*U_Tau/(2.0*Cp*T_Wall);
-        const su2double Beta = q_w*Lam_Visc_Wall/(Density_Wall*T_Wall*Conductivity_Wall*U_Tau);
-        const su2double Q    = sqrt(Beta*Beta + 4.0*Gam);
-        const su2double Phi  = asin(-1.0*Beta/Q);
+        const su2double U_Plus = VelTangMod / U_Tau;
 
-        /*--- Y+ defined by White & Christoph (compressibility and heat transfer) negative value for (2.0*Gam*U_Plus - Beta)/Q ---*/
+        /*--- Y+ defined by White & Christoph ---*/
+ 
+        const su2double kUp = kappa * U_Plus;
 
-        const su2double Y_Plus_White = exp((kappa/sqrt(Gam))*(asin((2.0*Gam*U_Plus - Beta)/Q) - Phi))*exp(-1.0*kappa*B);
+        // incompressible adiabatic result
+        const su2double Y_Plus_White = exp(kUp) * exp(-kappa * B);
 
         /*--- Spalding's universal form for the BL velocity with the
          *    outer velocity form of White & Christoph above. ---*/
-        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));
+        Y_Plus = U_Plus + Y_Plus_White + (exp(-kappa * B)* (1.0 - kUp - 0.5 * kUp * kUp - kUp * kUp * kUp / 6.0));
 
-        const 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));
+        /*--- incompressible formulation ---*/
+        Eddy_Visc_Wall = Lam_Visc_Wall * kappa*exp(-kappa*B) * (exp(kUp) -1.0 - kUp - kUp * kUp / 2.0);
 
-        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 --- */
 
         diff = (Density_Wall * U_Tau * WallDistMod / Lam_Visc_Wall) - Y_Plus;
 
-        /* --- Gradient of function defined above --- */
+        /* --- Gradient of function defined above wrt U_Tau --- */
 
-        const 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;
+        const su2double dyp_dup = 1.0 + exp(-kappa * B) * (kappa * exp(kUp) - kappa - kUp - 0.5 * kUp * kUp);
+        const su2double dup_dutau = - U_Plus / U_Tau;
+        const su2double grad_diff = Density_Wall * WallDistMod / Lam_Visc_Wall - dyp_dup * dup_dutau;
 
         /* --- Newton Step --- */
 

SU2_CFD/src/solvers/CNSSolver.cpp

@@ -857,9 +857,6 @@ void CNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_container,
       for (auto iDim = 0u; iDim < nDim; iDim++)
         Vel[iDim] = nodes->GetVelocity(Point_Normal,iDim);
 
-      const su2double P_Normal = nodes->GetPressure(Point_Normal);
-      const su2double T_Normal = nodes->GetTemperature(Point_Normal);
-
       /*--- Compute the wall-parallel velocity at first point off the wall ---*/
 
       const su2double VelNormal = GeometryToolbox::DotProduct(int(MAXNDIM), Vel, UnitNormal);
@@ -877,7 +874,10 @@ void CNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_container,
 
       const su2double WallDistMod = GeometryToolbox::Norm(int(MAXNDIM), WallDist);
 
-      /*--- Initial value for wall temperature ---*/
+      su2double T_Wall = nodes->GetTemperature(iPoint);
+      const su2double Conductivity_Wall = nodes->GetThermalConductivity(iPoint);
+
+      /*--- If a wall temperature was given, we compute the local heat flux using k*dT/dn ---*/
 
       su2double q_w = 0.0;
 
@@ -888,14 +888,14 @@ void CNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_container,
       /*--- Extrapolate the pressure from the interior & compute the
        wall density using the equation of state ---*/
 
+      const su2double P_Normal = nodes->GetPressure(Point_Normal);
+      const su2double T_Normal = nodes->GetTemperature(Point_Normal);
+      const su2double P_Wall = P_Normal;
+
       /*--- Compressible formulation ---*/
 
-      su2double T_Wall = nodes->GetTemperature(iPoint);
-      const su2double P_Wall = P_Normal;
       su2double Density_Wall = P_Wall / (Gas_Constant * T_Wall);
       const su2double Lam_Visc_Normal = nodes->GetLaminarViscosity(Point_Normal);
-      const su2double Conductivity_Wall = nodes->GetThermalConductivity(iPoint);
-
 
       /*--- Compute the shear stress at the wall in the regular fashion
        *    by using the stress tensor on the surface ---*/
@@ -938,7 +938,7 @@ void CNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_container,
 
         const 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 1108 ---*/
 
         const su2double Gam  = Recovery*U_Tau*U_Tau/(2.0*Cp*T_Wall);
         const su2double Beta = q_w*Lam_Visc_Wall/(Density_Wall*T_Wall*Conductivity_Wall*U_Tau);

TestCases/serial_regression.py

@@ -592,7 +592,7 @@ def main():
     inc_turb_wallfunction_flatplate_sst.cfg_dir   = "wallfunctions/flatplate/incompressible_SST"
     inc_turb_wallfunction_flatplate_sst.cfg_file  = "turb_SST_flatplate.cfg"
     inc_turb_wallfunction_flatplate_sst.test_iter = 10
-    inc_turb_wallfunction_flatplate_sst.test_vals = [-6.561293, -5.730190, -6.305916, -4.230297, -7.164083, -2.047256, 10, -2.941818, 0.000974, 0.003195, 0.205750]
+    inc_turb_wallfunction_flatplate_sst.test_vals = [-6.561370, -5.730257, -6.305996, -4.230300, -7.164395, -2.047273, 10, -2.941859, 0.000974, 0.003195, 0.205820]
     inc_turb_wallfunction_flatplate_sst.su2_exec  = "SU2_CFD"
     inc_turb_wallfunction_flatplate_sst.timeout   = 1600
     inc_turb_wallfunction_flatplate_sst.tol       = 0.00001
@@ -603,7 +603,7 @@ def main():
     inc_turb_wallfunction_flatplate_sa.cfg_dir   = "wallfunctions/flatplate/incompressible_SA"
     inc_turb_wallfunction_flatplate_sa.cfg_file  = "turb_SA_flatplate.cfg"
     inc_turb_wallfunction_flatplate_sa.test_iter = 10
-    inc_turb_wallfunction_flatplate_sa.test_vals = [-6.561221, -5.718427, -6.306321, -4.230222, -9.586624, 10, -2.927239, 0.000858, 0.003794 ] #last 9 columns
+    inc_turb_wallfunction_flatplate_sa.test_vals = [-6.561303, -5.718581, -6.306403, -4.230225, -9.586904, 10, -2.927392, 0.000858, 0.003792 ] #last 9 columns
     inc_turb_wallfunction_flatplate_sa.su2_exec  = "SU2_CFD"
     inc_turb_wallfunction_flatplate_sa.timeout   = 1600
     inc_turb_wallfunction_flatplate_sa.tol       = 0.00001