Merge pull request #2152 from su2code/feature_solid-solid_cht

Solid to solid conjugate heat transfer with contact resistance

.github/workflows/regression.yml

@@ -209,7 +209,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}}
+          args: -b ${{github.ref}} -t feature_solid-solid_cht -c develop -s ${{matrix.testscript}}
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536
         with:

Common/include/CConfig.hpp

@@ -192,6 +192,7 @@ class CConfig {
   nMarker_Fluid_Load,             /*!< \brief Number of markers in which the flow load is computed/employed. */
   nMarker_Fluid_InterfaceBound,   /*!< \brief Number of fluid interface markers. */
   nMarker_CHTInterface,           /*!< \brief Number of conjugate heat transfer interface markers. */
+  nMarker_ContactResistance,      /*!< \brief Number of CHT interfaces with contact resistance. */
   nMarker_Inlet,                  /*!< \brief Number of inlet flow markers. */
   nMarker_Inlet_Species,          /*!< \brief Number of inlet species markers. */
   nSpecies_per_Inlet,             /*!< \brief Number of species defined per inlet markers. */
@@ -396,6 +397,7 @@ class CConfig {
   su2double **Periodic_RotCenter;            /*!< \brief Rotational center for each periodic boundary. */
   su2double **Periodic_RotAngles;            /*!< \brief Rotation angles for each periodic boundary. */
   su2double **Periodic_Translation;          /*!< \brief Translation vector for each periodic boundary. */
+  su2double *CHT_ContactResistance;          /*!< \brief Contact resistance values for each solid-solid CHT interface. */
   string *Marker_CfgFile_TagBound;           /*!< \brief Global index for markers using config file. */
   unsigned short *Marker_All_KindBC,         /*!< \brief Global index for boundaries using grid information. */
   *Marker_CfgFile_KindBC;                    /*!< \brief Global index for boundaries using config file. */
@@ -589,6 +591,7 @@ class CConfig {
   bool EulerPersson;       /*!< \brief Boolean to determine whether this is an Euler simulation with Persson shock capturing. */
   bool FSI_Problem = false,/*!< \brief Boolean to determine whether the simulation is FSI or not. */
   Multizone_Problem;       /*!< \brief Boolean to determine whether we are solving a multizone problem. */
+  //bool ContactResistance = false; /*!< \brief Apply contact resistance for conjugate heat transfer. */
   unsigned short nID_DV;   /*!< \brief ID for the region of FEM when computed using direct differentiation. */
 
   bool AD_Mode;             /*!< \brief Algorithmic Differentiation support. */
@@ -3662,6 +3665,13 @@ class CConfig {
    */
   unsigned short GetMarker_n_ZoneInterface(void) const { return nMarker_ZoneInterface; }
 
+  /*!
+   * \brief Get the contact resistance value of a specified interface.
+   * \param[in] val_interface interface index.
+   * \return Contact resistance value (zero by default).
+   */
+  su2double GetContactResistance(unsigned short val_interface) const { return (nMarker_ContactResistance > 0) ? CHT_ContactResistance[val_interface] : 0.0; }
+
   /*!
    * \brief Get the DV information for a marker <i>val_marker</i>.
    * \param[in] val_marker - 0 or 1 depending if the the marker is going to be affected by design variables.

Common/include/option_structure.hpp

@@ -409,6 +409,7 @@ enum ENUM_TRANSFER {
   CONJUGATE_HEAT_WEAKLY_FS          = 17,   /*!< \brief Conjugate heat transfer (between incompressible fluids and solids). */
   CONJUGATE_HEAT_SF                 = 18,   /*!< \brief Conjugate heat transfer (between solids and compressible fluids). */
   CONJUGATE_HEAT_WEAKLY_SF          = 19,   /*!< \brief Conjugate heat transfer (between solids and incompressible fluids). */
+  CONJUGATE_HEAT_SS                 = 20,   /*!< \brief Conjugate heat transfer (between two solids). */
 };
 
 /*!

Common/src/CConfig.cpp

@@ -1517,6 +1517,8 @@ void CConfig::SetConfig_Options() {
   addStringListOption("MARKER_ZONE_INTERFACE", nMarker_ZoneInterface, Marker_ZoneInterface);
   /*!\brief MARKER_CHT_INTERFACE \n DESCRIPTION: CHT interface boundary marker(s) \ingroup Config*/
   addStringListOption("MARKER_CHT_INTERFACE", nMarker_CHTInterface, Marker_CHTInterface);
+  /*!\brief CHT_INTERFACE_CONTACT_RESISTANCE: Thermal contact resistance values for each CHT inerface. \ingroup Config*/
+  addDoubleListOption("CHT_INTERFACE_CONTACT_RESISTANCE", nMarker_ContactResistance, CHT_ContactResistance);
   /* DESCRIPTION: Internal boundary marker(s) */
   addStringListOption("MARKER_INTERNAL", nMarker_Internal, Marker_Internal);
   /* DESCRIPTION: Custom boundary marker(s) */
@@ -3560,6 +3562,25 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
       SU2_MPI::Error(string("You probably want to set INC_ENERGY_EQUATION= YES for the fluid solver. \n"), CURRENT_FUNCTION);
   }
 
+  /*--- Check correctness and consistency of contact resistance options. ---*/
+  if (nMarker_ContactResistance > 0) {
+
+    /*--- Set constant contact resistance across CHT interfaces if a single value is provided. ---*/
+    if (nMarker_ContactResistance == 1) {
+      auto val_CHTInterface = CHT_ContactResistance[0];
+      delete [] CHT_ContactResistance;
+      CHT_ContactResistance = new su2double[nMarker_CHTInterface];
+      for (auto iCHTMarker=0u; iCHTMarker < nMarker_CHTInterface; iCHTMarker++)
+        CHT_ContactResistance[iCHTMarker] = val_CHTInterface;
+    }else if((nMarker_CHTInterface/2) != nMarker_ContactResistance){
+      SU2_MPI::Error("Number of CHT interfaces does not match number of contact resistances.", CURRENT_FUNCTION);
+    }
+    for (auto iCHTMarker=0u; iCHTMarker < nMarker_ContactResistance; iCHTMarker++){
+      if (CHT_ContactResistance[iCHTMarker] < 0)
+        SU2_MPI::Error("Contact resistance value should be positive.", CURRENT_FUNCTION);
+    }
+  }
+
   /*--- By default, in 2D we should use TWOD_AIRFOIL (independenly from the input file) ---*/
 
   if (val_nDim == 2) Geo_Description = TWOD_AIRFOIL;

SU2_CFD/include/interfaces/CInterface.hpp

@@ -219,4 +219,9 @@ class CInterface {
    */
   void GatherAverageValues(CSolver *donor_solution, CSolver *target_solution, unsigned short donorZone);
 
+  /*!
+   * \brief Set the contact resistance value for the solid-to-solid heat transfer interface.
+   * \param[in] val_contact_resistance - Contact resistance value in m^2/W
+   */
+  inline virtual void SetContactResistance(su2double val_contact_resistance) {};
 };

SU2_CFD/include/interfaces/cht/CConjugateHeatInterface.hpp

@@ -35,6 +35,7 @@
  * \ingroup Interfaces
  */
 class CConjugateHeatInterface : public CInterface {
+  su2double ContactResistance = 0; /*!<\brief Contact resistance value of the current inerface. */
 public:
   /*!
    * \brief Constructor of the class.
@@ -70,4 +71,10 @@ class CConjugateHeatInterface : public CInterface {
    */
   void SetTarget_Variable(CSolver *target_solution, CGeometry *target_geometry, const CConfig *target_config,
                           unsigned long Marker_Target, unsigned long Vertex_Target, unsigned long Point_Target) override;
+
+  /*!
+   * \brief Set the contact resistance value for the solid-to-solid heat transfer interface.
+   * \param[in] val_contact_resistance - Contact resistance value in m^2/W
+   */
+  void SetContactResistance(su2double val_contact_resistance) override { ContactResistance = val_contact_resistance; }
 };

SU2_CFD/src/drivers/CDriver.cpp

@@ -2494,19 +2494,21 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet
           }
         }
         else if (heat_donor || heat_target) {
-          if (heat_donor && heat_target)
-            SU2_MPI::Error("Conjugate heat transfer between solids is not implemented.", CURRENT_FUNCTION);
+          if (heat_donor && heat_target){
+            interface_type = CONJUGATE_HEAT_SS;
 
-          const auto fluidZone = heat_target? donor : target;
-
-          if (config[fluidZone]->GetEnergy_Equation() || (config[fluidZone]->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE)
-              || (config[fluidZone]->GetKind_FluidModel() == ENUM_FLUIDMODEL::FLUID_FLAMELET))
-            interface_type = heat_target? CONJUGATE_HEAT_FS : CONJUGATE_HEAT_SF;
-          else if (config[fluidZone]->GetWeakly_Coupled_Heat())
-            interface_type = heat_target? CONJUGATE_HEAT_WEAKLY_FS : CONJUGATE_HEAT_WEAKLY_SF;
-          else
-            interface_type = NO_TRANSFER;
+          } else {
 
+            const auto fluidZone = heat_target? donor : target;
+            if (config[fluidZone]->GetEnergy_Equation() || (config[fluidZone]->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE)
+                || (config[fluidZone]->GetKind_FluidModel() == ENUM_FLUIDMODEL::FLUID_FLAMELET))
+              interface_type = heat_target? CONJUGATE_HEAT_FS : CONJUGATE_HEAT_SF;
+            else if (config[fluidZone]->GetWeakly_Coupled_Heat())
+              interface_type = heat_target? CONJUGATE_HEAT_WEAKLY_FS : CONJUGATE_HEAT_WEAKLY_SF;
+            else
+              interface_type = NO_TRANSFER;
+          }
+
           if (interface_type != NO_TRANSFER) {
             auto nVar = 4;
             interface[donor][target] = new CConjugateHeatInterface(nVar, 0);

SU2_CFD/src/drivers/CMultizoneDriver.cpp

@@ -563,6 +563,9 @@ bool CMultizoneDriver::TransferData(unsigned short donorZone, unsigned short tar
         BroadcastData(SPECIES_SOL, SPECIES_SOL);
       }
       break;
+    case CONJUGATE_HEAT_SS:
+      BroadcastData(HEAT_SOL, HEAT_SOL);
+      break;
     case CONJUGATE_HEAT_FS:
       BroadcastData(FLOW_SOL, HEAT_SOL);
       break;

SU2_CFD/src/interfaces/CInterface.cpp

@@ -108,6 +108,11 @@ void CInterface::BroadcastData(const CInterpolator& interpolator,
     su2activematrix sendDonorVar(nLocalVertexDonor, nVar);
 
     if (markDonor >= 0) {
+
+      /*--- Apply contact resistance if specified. ---*/
+
+      SetContactResistance(donor_config->GetContactResistance(iMarkerInt));
+
       for (auto iVertex = 0ul, iSend = 0ul; iVertex < donor_geometry->GetnVertex(markDonor); iVertex++) {
         const auto iPoint = donor_geometry->vertex[markDonor][iVertex]->GetNode();
 

SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp

@@ -132,8 +132,10 @@ void CConjugateHeatInterface::GetDonor_Variable(CSolver *donor_solution, CGeomet
     if ((donor_config->GetKind_CHT_Coupling() == CHT_COUPLING::DIRECT_TEMPERATURE_ROBIN_HEATFLUX) ||
         (donor_config->GetKind_CHT_Coupling() == CHT_COUPLING::AVERAGED_TEMPERATURE_ROBIN_HEATFLUX)) {
 
+      /*--- Apply contact resistance to solid-to-solid heat transfer boundary ---*/
       const su2double rho_cp_solid = donor_config->GetSpecific_Heat_Cp()*donor_config->GetMaterialDensity(0);
-      conductivity_over_dist  = thermal_diffusivity*rho_cp_solid/dist;
+      thermal_conductivity = thermal_diffusivity * rho_cp_solid;
+      conductivity_over_dist  = thermal_conductivity/(dist + thermal_conductivity * ContactResistance);
     }
   }
 

TestCases/parallel_regression.py

@@ -1339,6 +1339,7 @@ def main():
     sp_pinArray_3d_cht_mf_hf_tp.multizone = True
     test_list.append(sp_pinArray_3d_cht_mf_hf_tp)
 
+
     ##########################
     ###   Python wrapper   ###
     ##########################

TestCases/tutorials.py

@@ -64,6 +64,15 @@ def main():
     cht_incompressible.multizone = True
     test_list.append(cht_incompressible)
 
+    # Solid-to-solid and solid-to-fluid CHT with contact resistance
+    cht_CR           = TestCase('cht_solid_solid')
+    cht_CR.cfg_dir   = "../Tutorials/multiphysics/contact_resistance_cht"
+    cht_CR.cfg_file  = "master.cfg"
+    cht_CR.test_iter = 80
+    cht_CR.test_vals = [ -8.857438, -9.377593, -10.097769, -2.122358]
+    cht_CR.multizone = True
+    test_list.append(cht_CR)
+
     ### Incompressible Flow
 
     # 2D pin case massflow periodic with heatflux BC and prescribed extracted outlet heat