Merge pull request idaholab#19347 from GiudGiud/PR_functorthermalBC

Add functor versions of some FV BCs

framework/include/fvbcs/FVFluxBC.h

@@ -71,6 +71,16 @@ class FVFluxBC : public FVBoundaryCondition,
    */
   std::pair<SubdomainID, SubdomainID> faceArgSubdomains() const;
 
+  /**
+   * Determine the single sided face argument when evaluating a functor on a face.
+   * This is used to perform evluations of material properties with the actual face values of
+   * their dependences, rather than interpolate the material property to the boundary.
+   * @param limiter_type the limiter type, to be specified if more than the default average
+   *        interpolation is required for the parameters of the functor
+   */
+  std::tuple<const FaceInfo *, Moose::FV::LimiterType, bool, SubdomainID> singleSidedFaceArg(
+      Moose::FV::LimiterType limiter_type = Moose::FV::LimiterType::CentralDifference) const;
+
 private:
   /**
    * This creates a tuple of an element, \p FaceInfo, and subdomain ID. The element returned will

framework/src/fvbcs/FVFluxBC.C

@@ -262,3 +262,15 @@ FVFluxBC::faceArgSubdomains() const
 {
   return std::make_pair(std::get<2>(elemFromFace()), std::get<2>(neighborFromFace()));
 }
+
+std::tuple<const FaceInfo *, Moose::FV::LimiterType, bool, SubdomainID>
+FVFluxBC::singleSidedFaceArg(Moose::FV::LimiterType limiter_type) const
+{
+  const bool use_elem = _face_info->faceType(_var.name()) == FaceInfo::VarFaceNeighbors::ELEM;
+
+  if (use_elem)
+    return std::make_tuple(_face_info, limiter_type, true, _face_info->elem().subdomain_id());
+  else
+    return std::make_tuple(
+        _face_info, limiter_type, true, _face_info->neighborPtr()->subdomain_id());
+}

modules/heat_conduction/doc/content/source/fvbcs/FunctorThermalResistanceBC.md

@@ -0,0 +1,12 @@
+# FunctorThermalResistanceBC
+
+!syntax description /FVBCs/FunctorThermalResistanceBC
+
+The FunctorThermalResistanceBC is simply the [FVThermalResistanceBC.md] using functor
+material properties. See its documentation for more details.
+
+!syntax parameters /FVBCs/FunctorThermalResistanceBC
+
+!syntax inputs /FVBCs/FunctorThermalResistanceBC
+
+!syntax children /FVBCs/FunctorThermalResistanceBC

modules/heat_conduction/include/fvbcs/FVThermalResistanceBC.h

@@ -37,8 +37,10 @@ class FVThermalResistanceBC : public FVFluxBC
 protected:
   virtual ADReal computeQpResidual() override;
 
+  /// Computes the parallel heat flux resistance for a combined radiation-convection boundary
   void computeParallelResistance();
 
+  /// Computes the serial resistance of multiple conductive layers
   void computeConductionResistance();
 
   /// Whether to use a `cylindrical` or `cartesian` form for the thermal resistances.

modules/heat_conduction/include/fvbcs/FunctorThermalResistanceBC.h

@@ -0,0 +1,92 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "FVFluxBC.h"
+
+/**
+ * This BC applies a heat flux to a boundary, where the heat flux is
+ * determined using series conduction resistances, and parallel convection
+ * and radiation resistances at the surface. It is assumed that this BC is
+ * applied to an equation where `_u` represents a temperature. Note that this
+ * BC is only valid for pseudo steady-state problems and with no heat sources in
+ * the conduction layers. It is also assumed that the surface is isothermal - if
+ * this is not the case, this BC essentially assumes parallel heat transfer
+ * between each quadrature point and the surface, with no interaction between
+ * adjacent quadrature points. This will underpredict the heat transfer if the
+ * surface is not isothermal.
+ *
+ * When a radiation resistance at the surface is included, the solution requires
+ * an iterative procedure because the effective heat transfer coefficient for
+ * radiative heat transfer depends on the surface temperature. This iteration is
+ * performed with a simple underrelaxation method.
+ */
+class FunctorThermalResistanceBC : public FVFluxBC
+{
+public:
+  FunctorThermalResistanceBC(const InputParameters & parameters);
+  static InputParameters validParams();
+
+protected:
+  virtual ADReal computeQpResidual() override;
+
+  /// Computes the parallel heat flux resistance for a combined radiation-convection boundary
+  void computeParallelResistance();
+
+  /// Computes the serial resistance of multiple conductive layers
+  void computeConductionResistance();
+
+  /// Whether to use a `cylindrical` or `cartesian` form for the thermal resistances.
+  /// When using a cylindrical geometry, the inner-most radius must be provided.
+  const Moose::CoordinateSystemType _geometry;
+
+  /// Radius corresponding to the cylindrical surface (when using a cylindrical geometry)
+  const Real _inner_radius;
+
+  /// temperature variable
+  const Moose::Functor<ADReal> & _T;
+
+  /// ambient temperature for convection and radiation heat transfer
+  const Real _T_ambient;
+
+  /// thermal conductivities for each conduction layer, listed in order closest to
+  /// the boundary
+  const std::vector<Real> & _k;
+
+  /// thicknesses for each conduction layer, listed in order closest to the boundary
+  const std::vector<Real> & _dx;
+
+  /// convective heat transfer coefficient
+  const Moose::Functor<ADReal> & _h;
+
+  /// boundary emissivity
+  const Real & _emissivity;
+
+  /// maximum number of iterations (when radiative heat transfer is included)
+  const unsigned int & _max_iterations;
+
+  /// tolerance of iterations (when radiative heat transfer is included)
+  const Real & _tolerance;
+
+  /// underrelaxation factor (when radiative heat transfer is included)
+  const Real & _alpha;
+
+  /// surface temperature
+  ADReal _T_surface;
+
+  /// outer radius of surface
+  ADReal _outer_radius;
+
+  /// conduction thermal resistance
+  ADReal _conduction_resistance;
+
+  /// parallel convection and radiation thermal resistance
+  ADReal _parallel_resistance;
+};

modules/heat_conduction/src/fvbcs/FVThermalResistanceBC.C

@@ -70,8 +70,9 @@ FVThermalResistanceBC::FVThermalResistanceBC(const InputParameters & parameters)
     _parallel_resistance(0.0)
 {
   if (_k.size() != _dx.size())
-    mooseError("Number of specified thermal conductivities must match the number "
-               "of conduction layers!");
+    paramError("conduction_thicknesses",
+               "Number of specified thermal conductivities must match "
+               "the number of conduction layers!");
 
   if (_geometry == Moose::CoordinateSystemType::COORD_RZ)
     for (const auto & d : _dx)
@@ -87,7 +88,7 @@ FVThermalResistanceBC::computeConductionResistance()
 {
   Real r = _inner_radius;
 
-  for (std::size_t i = 0; i < _k.size(); ++i)
+  for (const auto i : index_range(_k))
   {
     switch (_geometry)
     {
@@ -123,7 +124,7 @@ FVThermalResistanceBC::computeQpResidual()
   computeParallelResistance();
 
   // other iteration requirements
-  Real iteration = 0;
+  unsigned int iteration = 0;
   ADReal norm = 2 * _tolerance;
   ADReal T_surface_previous;
 

modules/heat_conduction/src/fvbcs/FunctorThermalResistanceBC.C

@@ -0,0 +1,183 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "FunctorThermalResistanceBC.h"
+#include "HeatConductionNames.h"
+
+registerMooseObject("HeatConductionApp", FunctorThermalResistanceBC);
+
+InputParameters
+FunctorThermalResistanceBC::validParams()
+{
+  auto params = FVFluxBC::validParams();
+  params.addParam<MooseFunctorName>("temperature", "temperature variable");
+  params.addRequiredParam<Real>(HeatConduction::T_ambient, "constant ambient temperature");
+  params.addRequiredParam<MooseFunctorName>("htc", "heat transfer coefficient");
+
+  params.addRequiredRangeCheckedParam<Real>(HeatConduction::emissivity,
+                                            HeatConduction::emissivity + " >= 0.0 & " +
+                                                HeatConduction::emissivity + " <= 1.0",
+                                            "emissivity of the surface");
+
+  params.addRequiredParam<std::vector<Real>>(
+      "thermal_conductivities",
+      "vector of thermal conductivity values used for the conduction layers");
+  params.addRequiredParam<std::vector<Real>>("conduction_thicknesses",
+                                             "vector of conduction layer thicknesses");
+
+  MooseEnum geometry("cartesian cylindrical", "cartesian");
+  params.addParam<MooseEnum>("geometry", geometry, "type of geometry");
+  params.addRangeCheckedParam<Real>("inner_radius",
+                                    "inner_radius > 0.0",
+                                    "coordinate corresponding to the first resistance layer");
+
+  params.addRangeCheckedParam<Real>(
+      "step_size", 0.1, "step_size > 0.0", "underrelaxation step size");
+
+  params.addRangeCheckedParam<unsigned int>(
+      "max_iterations", 100, "max_iterations >= 0", "maximum iterations");
+
+  params.addRangeCheckedParam<Real>(
+      "tolerance", 1E-3, "tolerance > 0.0", "tolerance to converge iterations");
+  params.addClassDescription("Thermal resistance heat flux boundary condition for the "
+                             "fluid and solid energy equations");
+  return params;
+}
+
+FunctorThermalResistanceBC::FunctorThermalResistanceBC(const InputParameters & parameters)
+  : FVFluxBC(parameters),
+    _geometry(getParam<MooseEnum>("geometry").getEnum<Moose::CoordinateSystemType>()),
+    _inner_radius(
+        _geometry == Moose::CoordinateSystemType::COORD_RZ ? getParam<Real>("inner_radius") : 1.0),
+    _T(isParamValid("temperature") ? getFunctor<ADReal>("temperature")
+                                   : getFunctor<ADReal>("variable")),
+    _T_ambient(getParam<Real>(HeatConduction::T_ambient)),
+    _k(getParam<std::vector<Real>>("thermal_conductivities")),
+    _dx(getParam<std::vector<Real>>("conduction_thicknesses")),
+    _h(getFunctor<ADReal>(getParam<MooseFunctorName>("htc"))),
+    _emissivity(getParam<Real>(HeatConduction::emissivity)),
+    _max_iterations(getParam<unsigned int>("max_iterations")),
+    _tolerance(getParam<Real>("tolerance")),
+    _alpha(getParam<Real>("step_size")),
+    _T_surface(0.0),
+    _outer_radius(_inner_radius),
+    _conduction_resistance(0.0),
+    _parallel_resistance(0.0)
+{
+  if (_k.size() != _dx.size())
+    paramError("conduction_thicknesses",
+               "Number of specified thermal conductivities must match "
+               "the number of conduction layers!");
+
+  if (_geometry == Moose::CoordinateSystemType::COORD_RZ)
+    for (const auto & d : _dx)
+      _outer_radius += d;
+
+  // because the thermal conductivities are constant, we only need to compute
+  // the conduction resistance one time
+  computeConductionResistance();
+}
+
+void
+FunctorThermalResistanceBC::computeConductionResistance()
+{
+  Real r = _inner_radius;
+
+  for (const auto i : index_range(_k))
+  {
+    switch (_geometry)
+    {
+      case Moose::CoordinateSystemType::COORD_XYZ:
+        _conduction_resistance += _dx[i] / _k[i];
+        break;
+      case Moose::CoordinateSystemType::COORD_RZ:
+      {
+        _conduction_resistance += std::log((_dx[i] + r) / r) / _k[i];
+        r += _dx[i];
+        break;
+      }
+      default:
+        mooseError("Unhandled 'GeometryEnum' in 'FunctorThermalResistanceBC'!");
+    }
+  }
+}
+
+ADReal
+FunctorThermalResistanceBC::computeQpResidual()
+{
+  // Evaluate material properties on the face
+  const auto face_arg = singleSidedFaceArg();
+
+  // radiation resistance has to be solved iteratively, since we don't know the
+  // surface temperature. We do know that the heat flux in the conduction layers
+  // must match the heat flux in the parallel convection-radiation segment. For a
+  // first guess, take the surface temperature as the average of _T and T_ambient.
+  _T_surface = 0.5 * (_T(face_arg) + _T_ambient);
+
+  // total flux perpendicular to boundary
+  ADReal flux;
+
+  // resistance component representing the sum of the convection and radiation
+  // resistances, in parallel
+  computeParallelResistance();
+
+  // other iteration requirements
+  unsigned int iteration = 0;
+  ADReal norm = 2 * _tolerance;
+  ADReal T_surface_previous;
+
+  // iterate to find the approximate surface temperature needed for evaluating the
+  // radiation resistance. We only do this iteration if we have radiation transfer.
+  if (_emissivity > 1e-8)
+    while (norm > (_tolerance * _alpha))
+    {
+      T_surface_previous = _T_surface;
+
+      // compute the flux based on the conduction part of the circuit
+      flux = (_T(face_arg) - _T_surface) / _conduction_resistance;
+
+      computeParallelResistance();
+
+      // use the flux computed from the conduction half to update T_surface
+      _T_surface = flux * _parallel_resistance + _T_ambient;
+      _T_surface = _alpha * _T_surface + (1 - _alpha) * T_surface_previous;
+      norm = std::abs(_T_surface - T_surface_previous) / std::abs(T_surface_previous);
+
+      if (iteration == _max_iterations)
+      {
+        mooseWarning("Maximum number of iterations reached in 'FunctorThermalResistanceBC'!");
+        break;
+      }
+      else
+        iteration += 1;
+    }
+
+  // once we have determined T_surface, we can finally evaluate the complete
+  // resistance to find the overall heat flux. For Cartesian, dividing by the
+  // 'inner_radius' has no effect, but it is required for correct normalization
+  // for cylindrical geometries.
+  flux =
+      (_T(face_arg) - _T_ambient) / (_conduction_resistance + _parallel_resistance) / _inner_radius;
+  return flux;
+}
+
+void
+FunctorThermalResistanceBC::computeParallelResistance()
+{
+  const auto face_arg = singleSidedFaceArg();
+
+  // compute the parallel convection and radiation resistances, assuming they
+  // act on the same surface area size
+  ADReal hr = _emissivity * HeatConduction::Constants::sigma *
+              (_T_surface * _T_surface + _T_ambient * _T_ambient) * (_T_surface + _T_ambient);
+
+  // for Cartesian, dividing by the 'outer_radius' has no effect, but it is
+  // required for correct normalization for cylindrical geometries
+  _parallel_resistance = 1.0 / (hr + _h(face_arg)) / _outer_radius;
+}