measure the heat flux used in BCs instead of computing from the normal gradient

Author: pcarruscag

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -2363,6 +2363,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
   const su2double Beta = config->GetAoS() * PI_NUMBER / 180.0;
   const su2double RefLength = config->GetRefLength();
   const su2double RefHeatFlux = config->GetHeat_Flux_Ref();
+  const su2double RefTemperature = config->GetTemperature_Ref();
   const su2double Gas_Constant = config->GetGas_ConstantND();
   auto Origin = config->GetRefOriginMoment(0);
 
@@ -2396,6 +2397,8 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
     Marker_Tag = config->GetMarker_All_TagBound(iMarker);
     if (!config->GetViscous_Wall(iMarker)) continue;
 
+    const bool py_custom = config->GetMarker_All_PyCustom(iMarker);
+
     /*--- Obtain the origin for the moment computation for a particular marker ---*/
 
     const auto Monitoring = config->GetMarker_All_Monitoring(iMarker);
@@ -2506,26 +2509,55 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
       /*--- Compute total and maximum heat flux on the wall ---*/
 
-      su2double dTdn = -GeometryToolbox::DotProduct(nDim, Grad_Temp, UnitNormal);
-
-      if (!nemo){
-
+      if (!nemo) {
         if (FlowRegime == ENUM_REGIME::COMPRESSIBLE) {
-
           Cp = (Gamma / Gamma_Minus_One) * Gas_Constant;
           thermal_conductivity = Cp * Viscosity / Prandtl_Lam;
         }
         if (FlowRegime == ENUM_REGIME::INCOMPRESSIBLE) {
-          if (!energy) dTdn = 0.0;
           thermal_conductivity = nodes->GetThermalConductivity(iPoint);
         }
-        HeatFlux[iMarker][iVertex] = -thermal_conductivity * dTdn * RefHeatFlux;
 
+        if (config->GetMarker_All_KindBC(iMarker) == BC_TYPE::HEAT_FLUX) {
+          if (py_custom) {
+            HeatFlux[iMarker][iVertex] = geometry->GetCustomBoundaryHeatFlux(iMarker, iVertex);
+          } else {
+            HeatFlux[iMarker][iVertex] = config->GetWall_HeatFlux(Marker_Tag);
+            if (config->GetIntegrated_HeatFlux()) {
+              HeatFlux[iMarker][iVertex] /= geometry->GetSurfaceArea(config, iMarker);
+            }
+          }
+        } else if (config->GetMarker_All_KindBC(iMarker) == BC_TYPE::ISOTHERMAL) {
+          su2double Twall = 0.0;
+          if (py_custom) {
+            Twall = geometry->GetCustomBoundaryTemperature(iMarker, iVertex) / RefTemperature;
+          } else {
+            Twall = config->GetIsothermal_Temperature(Marker_Tag) / RefTemperature;
+          }
+          iPointNormal = geometry->vertex[iMarker][iVertex]->GetNormal_Neighbor();
+          Coord_Normal = geometry->nodes->GetCoord(iPointNormal);
+          su2double Vec_ij[MAXNDIM] = {0.0};
+          GeometryToolbox::Distance(nDim, Coord, Coord_Normal, Vec_ij);
+
+          /*--- Prevent divisions by 0 by limiting the normal projection. ---*/
+          const su2double dist_ij = fmax(
+            fabs(GeometryToolbox::DotProduct(int(MAXNDIM), Vec_ij, UnitNormal)),
+            fmax(0.05 * GeometryToolbox::Norm(int(MAXNDIM), Vec_ij), EPS));
+
+          const su2double There = nodes->GetTemperature(iPointNormal);
+
+          HeatFlux[iMarker][iVertex] = -thermal_conductivity * (There - Twall) / dist_ij * RefHeatFlux;
+        } else {
+          su2double dTdn = -GeometryToolbox::DotProduct(nDim, Grad_Temp, UnitNormal);
+          if (FlowRegime == ENUM_REGIME::INCOMPRESSIBLE && !energy) dTdn = 0.0;
+          HeatFlux[iMarker][iVertex] = -thermal_conductivity * dTdn * RefHeatFlux;
+        }
       } else {
 
         const auto& thermal_conductivity_tr = nodes->GetThermalConductivity(iPoint);
         const auto& thermal_conductivity_ve = nodes->GetThermalConductivity_ve(iPoint);
 
+        const su2double dTdn = -GeometryToolbox::DotProduct(nDim, Grad_Temp, UnitNormal);
         const su2double dTvedn = -GeometryToolbox::DotProduct(nDim, Grad_Temp_ve, UnitNormal);
 
         /*--- Surface energy balance: trans-rot heat flux, vib-el heat flux ---*/

SU2_CFD/src/solvers/CIncNSSolver.cpp

@@ -472,18 +472,21 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con
 
       const auto Coord_i = geometry->nodes->GetCoord(iPoint);
       const auto Coord_j = geometry->nodes->GetCoord(Point_Normal);
-      su2double Edge_Vector[MAXNDIM];
+      su2double Edge_Vector[MAXNDIM] = {0.0};
       GeometryToolbox::Distance(nDim, Coord_j, Coord_i, Edge_Vector);
-      su2double dist_ij_2 = GeometryToolbox::SquaredNorm(nDim, Edge_Vector);
-      su2double dist_ij = sqrt(dist_ij_2);
+
+      /*--- Prevent divisions by 0 by limiting the normal projection. ---*/
+      const su2double dist_ij = fmax(
+        fabs(GeometryToolbox::DotProduct(int(MAXNDIM), Edge_Vector, Normal)) / Area,
+        fmax(0.05 * GeometryToolbox::Norm(int(MAXNDIM), Edge_Vector), EPS));
 
       /*--- Compute the normal gradient in temperature using Twall ---*/
 
-      su2double dTdn = -(nodes->GetTemperature(Point_Normal) - Twall)/dist_ij;
+      const su2double dTdn = -(nodes->GetTemperature(Point_Normal) - Twall)/dist_ij;
 
       /*--- Get thermal conductivity ---*/
 
-      su2double thermal_conductivity = nodes->GetThermalConductivity(iPoint);
+      const su2double thermal_conductivity = nodes->GetThermalConductivity(iPoint);
 
       /*--- Apply a weak boundary condition for the energy equation.
       Compute the residual due to the prescribed heat flux. ---*/
@@ -493,10 +496,7 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con
       /*--- Jacobian contribution for temperature equation. ---*/
 
       if (implicit) {
-        su2double proj_vector_ij = 0.0;
-        if (dist_ij_2 > 0.0)
-          proj_vector_ij = GeometryToolbox::DotProduct(nDim, Edge_Vector, Normal) / dist_ij_2;
-        Jacobian.AddVal2Diag(iPoint, nDim+1, thermal_conductivity*proj_vector_ij);
+        Jacobian.AddVal2Diag(iPoint, nDim+1, thermal_conductivity * Area / dist_ij);
       }
       break;
     } // switch

SU2_CFD/src/solvers/CNSSolver.cpp

@@ -672,22 +672,25 @@ void CNSSolver::BC_Isothermal_Wall_Generic(CGeometry *geometry, CSolver **solver
 
     const auto Coord_i = geometry->nodes->GetCoord(iPoint);
     const auto Coord_j = geometry->nodes->GetCoord(Point_Normal);
+    su2double Vec_ij[MAXNDIM] = {0.0};
+    GeometryToolbox::Distance(nDim, Coord_i, Coord_j, Vec_ij);
 
-    su2double dist_ij = GeometryToolbox::Distance(nDim, Coord_i, Coord_j);
+    /*--- Prevent divisions by 0 by limiting the normal projection. ---*/
+    const su2double dist_ij = fmax(
+      fabs(GeometryToolbox::DotProduct(int(MAXNDIM), Vec_ij, UnitNormal)),
+      fmax(0.05 * GeometryToolbox::Norm(int(MAXNDIM), Vec_ij), EPS));
 
     /*--- Store the corrected velocity at the wall which will
      be zero (v = 0), unless there is grid motion (v = u_wall)---*/
 
     if (dynamic_grid) {
       nodes->SetVelocity_Old(iPoint, geometry->nodes->GetGridVel(iPoint));
-    }
-    else {
+    } else {
       su2double zero[MAXNDIM] = {0.0};
       nodes->SetVelocity_Old(iPoint, zero);
     }
 
-    for (auto iDim = 0u; iDim < nDim; iDim++)
-      LinSysRes(iPoint, iDim+1) = 0.0;
+    for (auto iDim = 0u; iDim < nDim; iDim++) LinSysRes(iPoint, iDim+1) = 0.0;
     nodes->SetVel_ResTruncError_Zero(iPoint);
 
     /*--- Get transport coefficients ---*/
@@ -708,8 +711,7 @@ void CNSSolver::BC_Isothermal_Wall_Generic(CGeometry *geometry, CSolver **solver
     if (cht_mode) {
       Twall = GetCHTWallTemperature(config, val_marker, iVertex, dist_ij,
                                     thermal_conductivity, There, Temperature_Ref);
-    }
-    else if (config->GetMarker_All_PyCustom(val_marker)) {
+    } else if (config->GetMarker_All_PyCustom(val_marker)) {
       Twall = geometry->GetCustomBoundaryTemperature(val_marker, iVertex) / Temperature_Ref;
     }