Compressible CHT fixes

Common/include/config_structure.hpp

@@ -532,7 +532,8 @@ class CConfig {
   Kind_Matrix_Coloring,         /*!< \brief Type of matrix coloring for sparse Jacobian computation. */
   Kind_Solver_Fluid_FSI,		/*!< \brief Kind of solver for the fluid in FSI applications. */
   Kind_Solver_Struc_FSI,		/*!< \brief Kind of solver for the structure in FSI applications. */
-  Kind_BGS_RelaxMethod;				/*!< \brief Kind of relaxation method for Block Gauss Seidel method in FSI problems. */
+  Kind_BGS_RelaxMethod,     /*!< \brief Kind of relaxation method for Block Gauss Seidel method in FSI problems. */
+  Kind_CHT_Coupling;        /*!< \brief Kind of coupling method used at CHT interfaces. */
   bool ReconstructionGradientRequired; /*!< \brief Enable or disable a second gradient calculation for upwind reconstruction only. */
   bool LeastSquaresRequired;  /*!< \brief Enable or disable memory allocation for least-squares gradient methods. */
   bool Energy_Equation;         /*!< \brief Solve the energy equation for incompressible flows. */
@@ -903,7 +904,6 @@ class CConfig {
   bool Sine_Load;                 /*!< \brief option for sine load */
   su2double *SineLoad_Coeff;      /*!< \brief Stores the load coefficient */
   su2double Thermal_Diffusivity;  /*!< \brief Thermal diffusivity used in the heat solver. */
-  bool CHT_Robin;                 /*!< \brief Option for boundary condition method at CHT interfaces. */
   su2double Cyclic_Pitch,         /*!< \brief Cyclic pitch for rotorcraft simulations. */
   Collective_Pitch;               /*!< \brief Collective pitch for rotorcraft simulations. */
   su2double Mach_Motion;          /*!< \brief Mach number based on mesh velocity and freestream quantities. */
@@ -8985,10 +8985,10 @@ class CConfig {
   bool GetWeakly_Coupled_Heat(void);
 
   /*!
-   * \brief Get the boundary condition method for CHT.
-   * \return YES if Robin BC is used.
+   * \brief Get the CHT couling method.
+   * \return Kind of the method.
    */
-  bool GetCHT_Robin(void);
+  unsigned short GetKind_CHT_Coupling(void) const;
 
   /*!
    * \brief Check if values passed to the BC_HeatFlux-Routine are already integrated.

Common/include/config_structure.inl

@@ -1917,7 +1917,7 @@ inline bool CConfig::GetJacobian_Spatial_Discretization_Only(void) {return Jacob
 
 inline bool CConfig::GetWeakly_Coupled_Heat(void) { return Weakly_Coupled_Heat; }
 
-inline bool CConfig::GetCHT_Robin(void) { return CHT_Robin; }
+inline unsigned short CConfig::GetKind_CHT_Coupling(void) const { return Kind_CHT_Coupling; }
 
 inline bool CConfig::GetIntegrated_HeatFlux(void) { return Integrated_HeatFlux; }
 

Common/include/option_structure.hpp

@@ -1112,6 +1112,23 @@ static const map<string, ENUM_DVFEA> DVFEA_Map = CCreateMap<string, ENUM_DVFEA>
 ("DEAD_WEIGHT", DEAD_WEIGHT)
 ("ELECTRIC_FIELD", ELECTRIC_FIELD);
 
+/*!
+ * \brief Kinds of coupling methods at CHT interfaces.
+ * The first (temperature) part determines the BC method on the fluid side, the second (heatflux) part determines
+ * the BC method on the solid side of the CHT interface.
+ */
+enum ENUM_CHT_COUPLING {
+  DIRECT_TEMPERATURE_NEUMANN_HEATFLUX = 0,
+  AVERAGED_TEMPERATURE_NEUMANN_HEATFLUX = 1,
+  DIRECT_TEMPERATURE_ROBIN_HEATFLUX = 2,
+  AVERAGED_TEMPERATURE_ROBIN_HEATFLUX = 3
+};
+static const map<string, ENUM_CHT_COUPLING> CHT_Coupling_Map = CCreateMap<string, ENUM_CHT_COUPLING>
+("DIRECT_TEMPERATURE_NEUMANN_HEATFLUX", DIRECT_TEMPERATURE_NEUMANN_HEATFLUX)
+("AVERAGED_TEMPERATURE_NEUMANN_HEATFLUX", AVERAGED_TEMPERATURE_NEUMANN_HEATFLUX)
+("DIRECT_TEMPERATURE_ROBIN_HEATFLUX", DIRECT_TEMPERATURE_ROBIN_HEATFLUX)
+("AVERAGED_TEMPERATURE_ROBIN_HEATFLUX", AVERAGED_TEMPERATURE_ROBIN_HEATFLUX);
+
 /*!
  * \brief types Riemann boundary treatments
  */

Common/src/config_structure.cpp

@@ -2268,9 +2268,9 @@ void CConfig::SetConfig_Options() {
   /*!\par CONFIG_CATEGORY: Heat solver \ingroup Config*/
   /*--- options related to the heat solver ---*/
 
-  /* DESCRIPTION: Use Robin (default) or Neumann BC at CHT interface. */
+  /* DESCRIPTION: CHT interface coupling methods */
   /*  Options: NO, YES \ingroup Config */
-  addBoolOption("CHT_ROBIN", CHT_Robin, true);
+  addEnumOption("CHT_COUPLING_METHOD", Kind_CHT_Coupling, CHT_Coupling_Map, DIRECT_TEMPERATURE_ROBIN_HEATFLUX);
 
   /* DESCRIPTION: Thermal diffusivity constant */
   addDoubleOption("THERMAL_DIFFUSIVITY", Thermal_Diffusivity, 1.172E-5);

SU2_CFD/src/drivers/CDriver.cpp

@@ -3547,7 +3547,7 @@ void CDriver::Interface_Preprocessing(CConfig **config, CSolver***** solver, CGe
       else if (fluid_donor && heat_target) {
         nVarTransfer = 0;
         nVar = 4;
-        if(config[donorZone]->GetEnergy_Equation())
+        if(config[donorZone]->GetEnergy_Equation() || (config[donorZone]->GetKind_Regime() == COMPRESSIBLE))
           interface_types[donorZone][targetZone] = CONJUGATE_HEAT_FS;
         else if (config[donorZone]->GetWeakly_Coupled_Heat())
           interface_types[donorZone][targetZone] = CONJUGATE_HEAT_WEAKLY_FS;
@@ -3558,7 +3558,7 @@ void CDriver::Interface_Preprocessing(CConfig **config, CSolver***** solver, CGe
       else if (heat_donor && fluid_target) {
         nVarTransfer = 0;
         nVar = 4;
-        if(config[targetZone]->GetEnergy_Equation())
+        if(config[targetZone]->GetEnergy_Equation() || (config[targetZone]->GetKind_Regime() == COMPRESSIBLE))
           interface_types[donorZone][targetZone] = CONJUGATE_HEAT_SF;
         else if (config[targetZone]->GetWeakly_Coupled_Heat())
           interface_types[donorZone][targetZone] = CONJUGATE_HEAT_WEAKLY_SF;

SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp

@@ -53,19 +53,17 @@ void CConjugateHeatInterface::GetDonor_Variable(CSolver *donor_solution, CGeomet
 
   /*--- Check whether the current zone is a solid zone or a fluid zone ---*/
 
-  bool flow = ((donor_config->GetKind_Solver() == NAVIER_STOKES)
+  bool compressible_flow = ((donor_config->GetKind_Solver() == NAVIER_STOKES)
                || (donor_config->GetKind_Solver() == RANS)
                || (donor_config->GetKind_Solver() == DISC_ADJ_NAVIER_STOKES)
-               || (donor_config->GetKind_Solver() == DISC_ADJ_RANS)
-               || (donor_config->GetKind_Solver() == INC_NAVIER_STOKES)
+               || (donor_config->GetKind_Solver() == DISC_ADJ_RANS));
+  bool incompressible_flow = ((donor_config->GetKind_Solver() == INC_NAVIER_STOKES)
                || (donor_config->GetKind_Solver() == INC_RANS)
                || (donor_config->GetKind_Solver() == DISC_ADJ_INC_NAVIER_STOKES)
-               || (donor_config->GetKind_Solver() == DISC_ADJ_INC_RANS));
-
-  bool compressible_flow    = (donor_config->GetKind_Regime() == COMPRESSIBLE) && flow;
-  bool incompressible_flow  = (donor_config->GetEnergy_Equation()) && flow;
-  bool heat_equation        = (donor_config->GetKind_Solver() == HEAT_EQUATION_FVM
-                               || donor_config->GetKind_Solver() == DISC_ADJ_HEAT);
+               || (donor_config->GetKind_Solver() == DISC_ADJ_INC_RANS)
+               && (donor_config->GetEnergy_Equation()));
+  bool heat_equation = (donor_config->GetKind_Solver() == HEAT_EQUATION_FVM
+               || donor_config->GetKind_Solver() == DISC_ADJ_HEAT);
 
   Coord         = donor_geometry->node[Point_Donor]->GetCoord();
 
@@ -130,7 +128,8 @@ void CConjugateHeatInterface::GetDonor_Variable(CSolver *donor_solution, CGeomet
     thermal_conductivityND  = Cp*(laminar_viscosity/Prandtl_Lam);
     heat_flux_density       = thermal_conductivityND*dTdn;
 
-    if (donor_config->GetCHT_Robin()) {
+    if ((donor_config->GetKind_CHT_Coupling() == DIRECT_TEMPERATURE_ROBIN_HEATFLUX) ||
+        (donor_config->GetKind_CHT_Coupling() == AVERAGED_TEMPERATURE_ROBIN_HEATFLUX)) {
 
       thermal_conductivity    = thermal_conductivityND*donor_config->GetViscosity_Ref();
       conductivity_over_dist  = thermal_conductivity/dist;
@@ -143,7 +142,8 @@ void CConjugateHeatInterface::GetDonor_Variable(CSolver *donor_solution, CGeomet
     thermal_conductivityND  = donor_solution->GetNodes()->GetThermalConductivity(iPoint);
     heat_flux_density       = thermal_conductivityND*dTdn;
 
-    if (donor_config->GetCHT_Robin()) {
+    if ((donor_config->GetKind_CHT_Coupling() == DIRECT_TEMPERATURE_ROBIN_HEATFLUX) ||
+        (donor_config->GetKind_CHT_Coupling() == AVERAGED_TEMPERATURE_ROBIN_HEATFLUX)) {
 
       switch (donor_config->GetKind_ConductivityModel()) {
 
@@ -167,7 +167,9 @@ void CConjugateHeatInterface::GetDonor_Variable(CSolver *donor_solution, CGeomet
     thermal_diffusivity     = donor_config->GetThermalDiffusivity_Solid();
     heat_flux_density       = thermal_diffusivity*dTdn;
 
-    if (donor_config->GetCHT_Robin()) {
+    if ((donor_config->GetKind_CHT_Coupling() == DIRECT_TEMPERATURE_ROBIN_HEATFLUX) ||
+        (donor_config->GetKind_CHT_Coupling() == AVERAGED_TEMPERATURE_ROBIN_HEATFLUX)) {
+
       rho_cp_solid            = donor_config->GetSpecific_Heat_Cp()*donor_config->GetDensity_Solid();
       conductivity_over_dist  = thermal_diffusivity*rho_cp_solid/dist;
     }
@@ -180,7 +182,8 @@ void CConjugateHeatInterface::GetDonor_Variable(CSolver *donor_solution, CGeomet
 
   /*--- We only need these for the Robin BC option ---*/
 
-  if (donor_config->GetCHT_Robin()) {
+  if ((donor_config->GetKind_CHT_Coupling() == DIRECT_TEMPERATURE_ROBIN_HEATFLUX) ||
+      (donor_config->GetKind_CHT_Coupling() == AVERAGED_TEMPERATURE_ROBIN_HEATFLUX)) {
 
     Donor_Variable[2] = conductivity_over_dist;
     Donor_Variable[3] = Tnormal*donor_config->GetTemperature_Ref();
@@ -201,7 +204,8 @@ void CConjugateHeatInterface::SetTarget_Variable(CSolver *target_solution, CGeom
   target_solution->SetConjugateHeatVariable(Marker_Target, Vertex_Target, 1,
                                             target_config->GetRelaxation_Factor_CHT(), Target_Variable[1]);
 
-  if (target_config->GetCHT_Robin()) {
+  if ((target_config->GetKind_CHT_Coupling() == DIRECT_TEMPERATURE_ROBIN_HEATFLUX) ||
+      (target_config->GetKind_CHT_Coupling() == AVERAGED_TEMPERATURE_ROBIN_HEATFLUX)) {
 
     target_solution->SetConjugateHeatVariable(Marker_Target, Vertex_Target, 2,
                                               target_config->GetRelaxation_Factor_CHT(), Target_Variable[2]);

SU2_CFD/src/solver_direct_heat.cpp

@@ -1113,78 +1113,69 @@ void CHeatSolverFVM::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **s
 
   if (flow) {
 
-    for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
+    for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) {
 
-      if (config->GetMarker_All_KindBC(iMarker) == CHT_WALL_INTERFACE) {
+      iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
 
-        for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) {
-          iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
-
-          if (geometry->node[iPoint]->GetDomain()) {
+      if (geometry->node[iPoint]->GetDomain()) {
 
-            Normal = geometry->vertex[iMarker][iVertex]->GetNormal();
-            Area = 0.0;
-            for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim];
-            Area = sqrt (Area);
+        Normal = geometry->vertex[val_marker][iVertex]->GetNormal();
+        Area = 0.0;
+        for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim];
+        Area = sqrt (Area);
 
-            T_Conjugate = GetConjugateHeatVariable(iMarker, iVertex, 0)/Temperature_Ref;
+        T_Conjugate = GetConjugateHeatVariable(val_marker, iVertex, 0)/Temperature_Ref;
 
-            nodes->SetSolution_Old(iPoint,&T_Conjugate);
-            LinSysRes.SetBlock_Zero(iPoint, 0);
-            nodes->SetRes_TruncErrorZero(iPoint);
+        nodes->SetSolution_Old(iPoint,&T_Conjugate);
+        LinSysRes.SetBlock_Zero(iPoint, 0);
+        nodes->SetRes_TruncErrorZero(iPoint);
 
-            if (implicit) {
-              for (iVar = 0; iVar < nVar; iVar++) {
-                total_index = iPoint*nVar+iVar;
-                Jacobian.DeleteValsRowi(total_index);
-              }
-            }
+        if (implicit) {
+          for (iVar = 0; iVar < nVar; iVar++) {
+            total_index = iPoint*nVar+iVar;
+            Jacobian.DeleteValsRowi(total_index);
           }
         }
       }
     }
   }
   else {
 
-    for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
+    for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) {
 
-      if (config->GetMarker_All_KindBC(iMarker) == CHT_WALL_INTERFACE) {
+      iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
 
-        for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) {
-          iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
-
-          if (geometry->node[iPoint]->GetDomain()) {
+      if (geometry->node[iPoint]->GetDomain()) {
 
-            Normal = geometry->vertex[iMarker][iVertex]->GetNormal();
-            Area = 0.0;
-            for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim];
-            Area = sqrt(Area);
+        Normal = geometry->vertex[val_marker][iVertex]->GetNormal();
+        Area = 0.0;
+        for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim];
+        Area = sqrt(Area);
 
-            thermal_diffusivity = GetConjugateHeatVariable(iMarker, iVertex, 2)/rho_cp_solid;
+        thermal_diffusivity = GetConjugateHeatVariable(val_marker, iVertex, 2)/rho_cp_solid;
 
-            if (config->GetCHT_Robin()) {
+        if ((config->GetKind_CHT_Coupling() == DIRECT_TEMPERATURE_ROBIN_HEATFLUX) ||
+            (config->GetKind_CHT_Coupling() == AVERAGED_TEMPERATURE_ROBIN_HEATFLUX)) {
 
-              Tinterface        = nodes->GetSolution(iPoint,0);
-              Tnormal_Conjugate = GetConjugateHeatVariable(iMarker, iVertex, 3)/Temperature_Ref;
+          Tinterface        = nodes->GetSolution(iPoint,0);
+          Tnormal_Conjugate = GetConjugateHeatVariable(val_marker, iVertex, 3)/Temperature_Ref;
 
-              HeatFluxDensity   = thermal_diffusivity*(Tinterface - Tnormal_Conjugate);
-              HeatFlux          = HeatFluxDensity * Area;
-            }
-            else {
+          HeatFluxDensity   = thermal_diffusivity*(Tinterface - Tnormal_Conjugate);
+          HeatFlux          = HeatFluxDensity * Area;
+        }
+        else {
 
-              HeatFluxDensity = GetConjugateHeatVariable(iMarker, iVertex, 1)/config->GetHeat_Flux_Ref();
-              HeatFlux        = HeatFluxDensity*Area;
-            }
+          HeatFluxDensity = GetConjugateHeatVariable(val_marker, iVertex, 1)/config->GetHeat_Flux_Ref();
+          HeatFlux        = HeatFluxDensity*Area;
+        }
 
-            Res_Visc[0] = -HeatFlux;
-            LinSysRes.SubtractBlock(iPoint, Res_Visc);
+        Res_Visc[0] = -HeatFlux;
+        LinSysRes.SubtractBlock(iPoint, Res_Visc);
 
-            if (implicit) {
+        if (implicit) {
 
-              Jacobian_i[0][0] = thermal_diffusivity*Area;
-              Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i);
-            }
-          }
+          Jacobian_i[0][0] = thermal_diffusivity*Area;
+          Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i);
         }
       }
     }

SU2_CFD/src/solver_direct_mean.cpp

@@ -16003,6 +16003,8 @@ void CNSSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **solver
   su2double Gas_Constant = config->GetGas_ConstantND();
   su2double Cp = (Gamma / Gamma_Minus_One) * Gas_Constant;
 
+  su2double Temperature_Ref = config->GetTemperature_Ref();
+
   bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
 
   /*--- Identify the boundary ---*/
@@ -16092,17 +16094,31 @@ void CNSSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **solver
       /*--- Compute the normal gradient in temperature using Twall ---*/
 
       There = nodes->GetTemperature(Point_Normal);
-      Tconjugate = GetConjugateHeatVariable(val_marker, iVertex, 0);
+      Tconjugate = GetConjugateHeatVariable(val_marker, iVertex, 0)/Temperature_Ref;
+
+      if ((config->GetKind_CHT_Coupling() == AVERAGED_TEMPERATURE_NEUMANN_HEATFLUX) ||
+          (config->GetKind_CHT_Coupling() == AVERAGED_TEMPERATURE_ROBIN_HEATFLUX)) {
 
-      HF_FactorHere = thermal_conductivity*config->GetViscosity_Ref()/dist_ij;
-      HF_FactorConjugate = GetConjugateHeatVariable(val_marker, iVertex, 2);
+        /*--- Compute wall temperature from both temperatures ---*/
 
-      Twall = (There*HF_FactorHere + Tconjugate*HF_FactorConjugate)/(HF_FactorHere + HF_FactorConjugate);
+        HF_FactorHere = thermal_conductivity*config->GetViscosity_Ref()/dist_ij;
+        HF_FactorConjugate = GetConjugateHeatVariable(val_marker, iVertex, 2);
 
-      // this will be changed soon...
-      Twall = GetConjugateHeatVariable(val_marker, iVertex, 0);
+        Twall = (There*HF_FactorHere + Tconjugate*HF_FactorConjugate)/(HF_FactorHere + HF_FactorConjugate);
+        dTdn = -(There - Twall)/dist_ij;
+      }
+      else if ((config->GetKind_CHT_Coupling() == DIRECT_TEMPERATURE_NEUMANN_HEATFLUX) ||
+              (config->GetKind_CHT_Coupling() == DIRECT_TEMPERATURE_ROBIN_HEATFLUX)) {
 
-      dTdn = -(There - Twall)/dist_ij;
+        /*--- (Directly) Set wall temperature to conjugate temperature. ---*/
+
+        Twall = Tconjugate;
+        dTdn = -(There - Twall)/dist_ij;
+      }
+      else {
+
+        SU2_MPI::Error(string("Unknown CHT coupling method."), CURRENT_FUNCTION);
+      }
 
       /*--- Apply a weak boundary condition for the energy equation.
        Compute the residual due to the prescribed heat flux. ---*/
@@ -16281,7 +16297,6 @@ void CNSSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **solver
           Jacobian.DeleteValsRowi(total_index);
         }
       }
-
     }
   }
 }