Addition of catalytic isothermal wall boundary conditions in SU2-NEMO

.github/workflows/regression.yml

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

Common/include/CConfig.hpp

@@ -311,7 +311,6 @@ class CConfig {
   su2double *Isothermal_Temperature;         /*!< \brief Specified isothermal wall temperatures (static). */
   su2double *HeatTransfer_Coeff;             /*!< \brief Specified heat transfer coefficients. */
   su2double *HeatTransfer_WallTemp;          /*!< \brief Specified temperatures at infinity alongside heat transfer coefficients. */
-  su2double *Wall_Catalycity;                /*!< \brief Specified wall species mass-fractions for catalytic boundaries. */
   su2double *Heat_Flux;                      /*!< \brief Specified wall heat fluxes. */
   su2double *Roughness_Height;               /*!< \brief Equivalent sand grain roughness for the marker according to config file. */
   su2double *Displ_Value;                    /*!< \brief Specified displacement for displacement boundaries. */
@@ -1169,20 +1168,24 @@ class CConfig {
   unsigned short maxBasisDim,               /*!< \brief Maximum number of POD basis dimensions. */
   rom_save_freq;                            /*!< \brief Frequency of unsteady time steps to save. */
 
-  unsigned short nSpecies,                  /*!< \brief Number of transported species equations (for NEMO and species transport)*/
+  unsigned short nSpecies;                  /*!< \brief Number of transported species equations (for NEMO and species transport)*/
 
   /* other NEMO configure options*/
-  iWall_Catalytic,
-  nWall_Catalytic;                          /*!< \brief No of catalytic walls */
-  su2double *Gas_Composition,               /*!< \brief Initial mass fractions of flow [dimensionless] */
+  unsigned short nSpecies_Cat_Wall,         /*!< \brief No. of species for a catalytic wall. */
+  iWall_Catalytic,                          /*!< \brief Iterator over catalytic walls. */
+  nWall_Catalytic;                          /*!< \brief No. of catalytic walls. */
+  su2double *Gas_Composition,               /*!< \brief Initial mass fractions of flow [dimensionless]. */
+  *Supercatalytic_Wall_Composition,         /*!< \brief Supercatalytic wall mass fractions [dimensionless]. */
   pnorm_heat;                               /*!< \brief pnorm for heat-flux. */
   bool frozen,                              /*!< \brief Flag for determining if mixture is frozen. */
   ionization,                               /*!< \brief Flag for determining if free electron gas is in the mixture. */
   vt_transfer_res_limit,                    /*!< \brief Flag for determining if residual limiting for source term VT-transfer is used. */
-  monoatomic;                               /*!< \brief Flag for monoatomic mixture. */
+  monoatomic,                               /*!< \brief Flag for monoatomic mixture. */
+  Supercatalytic_Wall;                      /*!< \brief Flag for supercatalytic wall. */
   string GasModel,                          /*!< \brief Gas Model. */
   *Wall_Catalytic;                          /*!< \brief Pointer to catalytic walls. */
   TRANSCOEFFMODEL   Kind_TransCoeffModel;   /*!< \brief Transport coefficient Model for NEMO solver. */
+  su2double CatalyticEfficiency;            /*!< \brief Wall catalytic efficiency. */
 
   /*--- Additional species solver options ---*/
   bool Species_Clipping;           /*!< \brief Boolean that activates solution clipping for scalar transport. */
@@ -3029,13 +3032,13 @@ class CConfig {
 
   /*!
    * \brief Retrieves the number of periodic time instances for Harmonic Balance.
-   * \return: Number of periodic time instances for Harmonic Balance.
+   * \return Number of periodic time instances for Harmonic Balance.
    */
   unsigned short GetnTimeInstances(void) const { return nTimeInstances; }
 
   /*!
    * \brief Retrieves the period of oscillations to be used with Harmonic Balance.
-   * \return: Period for Harmonic Balance.
+   * \return Period for Harmonic Balance.
    */
   su2double GetHarmonicBalance_Period(void) const { return HarmonicBalance_Period; }
 
@@ -3724,6 +3727,12 @@ class CConfig {
    */
   string GetWall_Catalytic_TagBound(unsigned short val_marker) const { return Wall_Catalytic[val_marker]; }
 
+  /*!
+   * \brief Get wall catalytic efficiency.
+   * \return wall catalytic efficiency value.
+   */
+  su2double GetCatalytic_Efficiency(void) const { return CatalyticEfficiency; }
+
   /*!
    * \brief Fluid model that we are using.
    * \return Fluid model that we are using.
@@ -5211,22 +5220,22 @@ class CConfig {
   TIME_MARCHING GetTime_Marching() const { return TimeMarching; }
 
   /*!
-   * \brief Provides the number of species present in the plasma
-   * \return: The number of species present in the plasma, read from input file
+   * \brief Provides the number of species present in the gas mixture.
+   * \return The number of species present in the gas mixture.
    */
   unsigned short GetnSpecies() const { return nSpecies; }
 
-   /*!
-   * \brief Get the wall heat flux on a constant heat flux boundary.
-   * \return The heat flux.
+  /*!
+   * \brief Provides the gas mass fractions of the flow.
+   * \return Gas Mass fractions.
    */
-  const su2double *GetWall_Catalycity(void) const { return Wall_Catalycity; }
+  const su2double *GetGas_Composition(void) const { return Gas_Composition; }
 
   /*!
-   * \brief Provides the gas mass fractions of the flow
-   * \return: Gas Mass fractions
+   * \brief Provides the gas mass fractions at the wall for supercat wall.
+   * \return Supercat wall gas mass fractions.
    */
-  const su2double *GetGas_Composition(void) const { return Gas_Composition; }
+  const su2double *GetSupercatalytic_Wall_Composition(void) const { return Supercatalytic_Wall_Composition; }
 
   /*!
    * \brief Provides the restart information.
@@ -5304,6 +5313,11 @@ class CConfig {
    */
   bool GetMonoatomic(void) const { return monoatomic; }
 
+  /*!
+   * \brief Indicates whether supercatalytic wall is used.
+   */
+  bool GetSupercatalytic_Wall(void) const { return Supercatalytic_Wall; }
+
   /*!
    * \brief Information about computing and plotting the equivalent area distribution.
    * \return <code>TRUE</code> or <code>FALSE</code>  depending if we are computing the equivalent area.
@@ -6217,17 +6231,30 @@ class CConfig {
   bool GetViscous_Wall(unsigned short iMarker) const;
 
   /*!
-   * \brief Determines if problem is adjoint
-   * \return true if Adjoint
+   * \brief Determines whether a marker with index iMarker is a catalytic boundary.
+   * \param iMarker
+   * \return <TRUE> it marker with index iMarker is a catalytic boundary.
+   */
+  bool GetCatalytic_Wall(unsigned short iMarker) const;
+
+  /*!
+   * \brief Determines if problem is adjoint.
+   * \return true if Adjoint.
    */
   bool GetContinuous_Adjoint(void) const { return ContinuousAdjoint; }
 
   /*!
-   * \brief Determines if problem is viscous
-   * \return true if Viscous
+   * \brief Determines if problem is viscous.
+   * \return true if Viscous.
    */
   bool GetViscous(void) const { return Viscous; }
 
+  /*!
+   * \brief Determines if problem has catalytic walls.
+   * \return true if catalytic walls are present.
+   */
+  bool GetCatalytic(void) const { return nWall_Catalytic > 0; }
+
   /*!
    * \brief Provides the index of the solution in the container.
    * \param[in] val_eqsystem - Equation that is being solved.

Common/src/CConfig.cpp

@@ -1057,7 +1057,6 @@ void CConfig::SetPointersNull(void) {
 
   Kind_TimeNumScheme = EULER_IMPLICIT;
 
-  Gas_Composition = nullptr;
 }
 
 void CConfig::SetConfig_Options() {
@@ -1183,7 +1182,13 @@ void CConfig::SetConfig_Options() {
   addBoolOption("VT_RESIDUAL_LIMITING", vt_transfer_res_limit, false);
   /* DESCRIPTION: List of catalytic walls */
   addStringListOption("CATALYTIC_WALL", nWall_Catalytic, Wall_Catalytic);
-
+  /* DESCRIPTION: Specfify super-catalytic wall */
+  addBoolOption("SUPERCATALYTIC_WALL", Supercatalytic_Wall, false);
+  /* DESCRIPTION: Wall mass fractions for supercatalytic case */
+  addDoubleListOption("SUPERCATALYTIC_WALL_COMPOSITION", nSpecies_Cat_Wall, Supercatalytic_Wall_Composition);
+  /* DESCRIPTION: Specfify catalytic efficiency of wall if using gamma model */
+  addDoubleOption("CATALYTIC_EFFICIENCY", CatalyticEfficiency, 1.0);
+  /*!\brief MARKER_MONITORING\n DESCRIPTION: Marker(s) of the surface where evaluate the non-dimensional coefficients \ingroup Config*/
 
   /*--- Options related to VAN der WAALS MODEL and PENG ROBINSON ---*/
 
@@ -7796,6 +7801,17 @@ bool CConfig::GetViscous_Wall(unsigned short iMarker) const {
           Marker_All_KindBC[iMarker] == CHT_WALL_INTERFACE);
 }
 
+bool CConfig::GetCatalytic_Wall(unsigned short iMarker) const {
+
+  bool catalytic = false;
+  for (unsigned short iMarker_Catalytic = 0; iMarker_Catalytic < nWall_Catalytic; iMarker_Catalytic++){
+    string Catalytic_Tag = Wall_Catalytic[iMarker_Catalytic];
+    if (Catalytic_Tag == Marker_All_TagBound[iMarker]) { catalytic = true; break; }
+  }
+
+  return catalytic;
+}
+
 bool CConfig::GetSolid_Wall(unsigned short iMarker) const {
 
   return GetViscous_Wall(iMarker) ||

SU2_CFD/include/fluid/CNEMOGas.hpp

@@ -84,6 +84,8 @@ class CNEMOGas : public CFluidModel {
   Enthalpy_Formation,                    /*!< \brief Enthalpy of formation */
   Ref_Temperature;                       /*!< \brief Reference temperature for thermodynamic relations */
 
+  su2matrix<int> CatRecombTable;         /*!< \brief Table for catalytic wall recombination pairs. */
+
 public:
 
   /*!
@@ -258,4 +260,10 @@ class CNEMOGas : public CFluidModel {
    * \brief Get species formation enthalpy.
    */
   virtual const vector<su2double>& GetSpeciesFormationEnthalpy() = 0;
+
+  /*!
+   * \brief Get catalytic wall recombination indices and constants.
+   */
+  inline const su2matrix<int>& GetCatalyticRecombination() const {return CatRecombTable;}
+
 };

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -2437,9 +2437,10 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
   unsigned long iVertex, iPoint, iPointNormal;
   unsigned short iMarker, iMarker_Monitoring, iDim, jDim;
-  su2double Viscosity = 0.0, Area, Density = 0.0, GradTemperature = 0.0, WallDistMod, FrictionVel,
+  su2double Viscosity = 0.0, Area, Density = 0.0, WallDistMod, FrictionVel,
             UnitNormal[3] = {0.0}, TauElem[3] = {0.0}, Tau[3][3] = {{0.0}}, Cp,
-            thermal_conductivity, MaxNorm = 8.0, Grad_Vel[3][3] = {{0.0}}, Grad_Temp[3] = {0.0}, AxiFactor;
+            thermal_conductivity, MaxNorm = 8.0, Grad_Vel[3][3] = {{0.0}}, Grad_Temp[3] = {0.0},
+            Grad_Temp_ve[3] = {0.0}, AxiFactor;
   const su2double *Coord = nullptr, *Coord_Normal = nullptr, *Normal = nullptr;
   const su2double minYPlus = config->GetwallModel_MinYPlus();
 
@@ -2520,6 +2521,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
           Grad_Vel[iDim][jDim] = nodes->GetGradient_Primitive(iPoint, prim_idx.Velocity() + iDim, jDim);
         }
         Grad_Temp[iDim] = nodes->GetGradient_Primitive(iPoint, prim_idx.Temperature(), iDim);
+        if (nemo) Grad_Temp_ve[iDim] = nodes->GetGradient_Primitive(iPoint, prim_idx.Temperature_ve(), iDim);
       }
 
       Viscosity = nodes->GetLaminarViscosity(iPoint);
@@ -2591,35 +2593,53 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
       /*--- Compute total and maximum heat flux on the wall ---*/
 
-      /// TODO: Move these ifs to specialized functions.
+      su2double dTdn = -GeometryToolbox::DotProduct(nDim, Grad_Temp, UnitNormal);
+
       if (!nemo){
 
         if (FlowRegime == ENUM_REGIME::COMPRESSIBLE) {
-          GradTemperature = -GeometryToolbox::DotProduct(nDim, Grad_Temp, UnitNormal);
 
           Cp = (Gamma / Gamma_Minus_One) * Gas_Constant;
           thermal_conductivity = Cp * Viscosity / Prandtl_Lam;
         }
         if (FlowRegime == ENUM_REGIME::INCOMPRESSIBLE) {
-          if (energy)
-            GradTemperature = -GeometryToolbox::DotProduct(nDim, Grad_Temp, UnitNormal);
-
+          if (!energy) dTdn = 0.0;
           thermal_conductivity = nodes->GetThermalConductivity(iPoint);
         }
-        HeatFlux[iMarker][iVertex] = -thermal_conductivity * GradTemperature * RefHeatFlux;
+        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 auto& Grad_PrimVar            = nodes->GetGradient_Primitive(iPoint);
-
-        su2double dTn   = GeometryToolbox::DotProduct(nDim, Grad_PrimVar[prim_idx.Temperature()], UnitNormal);
-        su2double dTven = GeometryToolbox::DotProduct(nDim, Grad_PrimVar[prim_idx.Temperature_ve()], UnitNormal);
+
+        const su2double dTvedn = -GeometryToolbox::DotProduct(nDim, Grad_Temp_ve, UnitNormal);
+
+        /*--- Surface energy balance: trans-rot heat flux, vib-el heat flux ---*/
+        HeatFlux[iMarker][iVertex] = -(thermal_conductivity_tr*dTdn + thermal_conductivity_ve*dTvedn);
+
+        /*--- Compute enthalpy transport to surface due to mass diffusion ---*/
+        bool catalytic = config->GetCatalytic_Wall(iMarker);
+        if (catalytic){
+
+          const auto nSpecies = config->GetnSpecies();
+          const auto& Grad_PrimVar = nodes->GetGradient_Primitive(iPoint);
+          const auto& PrimVar = nodes->GetPrimitive(iPoint);
+          const auto& Ds = nodes->GetDiffusionCoeff(iPoint);
+          const auto& hs = nodes->GetEnthalpys(iPoint);
+          const su2double rho = PrimVar[prim_idx.Density()];
+
+          su2double sumJhs = 0.0;
+          for (auto iSpecies = 0u; iSpecies < nSpecies; iSpecies++) {
+            for (auto iDim = 0u; iDim < nDim; iDim++) {
+              su2double dYdn = 1.0/rho*(Grad_PrimVar[iSpecies][iDim] - PrimVar[iSpecies]*Grad_PrimVar[prim_idx.Density()][iDim]/rho);
+              sumJhs += rho*Ds[iSpecies]*hs[iSpecies]*dYdn*UnitNormal[iDim];
+            }
+          }
+          /*--- Surface energy balance: mass diffusion ---*/
+          HeatFlux[iMarker][iVertex] += sumJhs;
 
-        /*--- Surface energy balance: trans-rot heat flux, vib-el heat flux,
-        enthalpy transport due to mass diffusion ---*/
-        HeatFlux[iMarker][iVertex] = thermal_conductivity_tr*dTn + thermal_conductivity_ve*dTven;
+        }
       }
 
       /*--- Note that heat is computed at the

SU2_CFD/src/fluid/CMutationTCLib.cpp

@@ -38,6 +38,7 @@ CMutationTCLib::CMutationTCLib(const CConfig* config, unsigned short val_nDim):
   Cv_ks.resize(nEnergyEq*nSpecies,0.0);
   es.resize(nEnergyEq*nSpecies,0.0);
   omega_vec.resize(1,0.0);
+  CatRecombTable.resize(nSpecies,2) = 0;
 
   /*--- Set up inputs to define type of mixture in the Mutation++ library ---*/
 
@@ -69,6 +70,33 @@ CMutationTCLib::CMutationTCLib(const CConfig* config, unsigned short val_nDim):
   }
   else { nHeavy = nSpecies;   nEl = 0; }
 
+  /*--- Set up catalytic recombination table. ---*/
+  // Creation/Destruction (+1/-1), Index of monoatomic reactants
+  // Monoatomic species (N,O) recombine into diaatomic (N2, O2)
+  if (gas_model == "N2") {
+    CatRecombTable(0,0) =  1; CatRecombTable(0,1) = 1;
+    CatRecombTable(1,0) = -1; CatRecombTable(1,1) = 1;
+
+  } else if (gas_model == "air_5"){
+    CatRecombTable(0,0) = -1; CatRecombTable(0,1) = 0;
+    CatRecombTable(1,0) = -1; CatRecombTable(1,1) = 1;
+    CatRecombTable(2,0) =  0; CatRecombTable(2,1) = 4;
+    CatRecombTable(3,0) =  1; CatRecombTable(3,1) = 0;
+    CatRecombTable(4,0) =  1; CatRecombTable(4,1) = 1;
+
+  } else if (gas_model == "air_6") {
+    CatRecombTable(0,0) = -1; CatRecombTable(0,1) = 0;
+    CatRecombTable(1,0) = -1; CatRecombTable(1,1) = 1;
+    CatRecombTable(2,0) =  0; CatRecombTable(2,1) = 4;
+    CatRecombTable(3,0) =  1; CatRecombTable(3,1) = 0;
+    CatRecombTable(4,0) =  1; CatRecombTable(4,1) = 1;
+    CatRecombTable(5,0) =  0; CatRecombTable(5,1) = 4;
+
+  } else {
+    if (config->GetCatalytic())
+      SU2_MPI::Error("Cataylic wall recombination not implemented for specified Mutation gas model.", CURRENT_FUNCTION);
+  }
+
 }
 
 CMutationTCLib::~CMutationTCLib(){}