Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "FVThermalResistanceBC.h"
      11             : #include "HeatConductionNames.h"
      12             : 
      13             : registerMooseObject("HeatTransferApp", FVThermalResistanceBC);
      14             : 
      15             : InputParameters
      16         164 : FVThermalResistanceBC::validParams()
      17             : {
      18         164 :   auto params = FVFluxBC::validParams();
      19         328 :   params.addCoupledVar("temperature", "temperature variable");
      20         164 :   params.addRequiredParam<Real>(HeatConduction::T_ambient, "constant ambient temperature");
      21         328 :   params.addRequiredParam<MaterialPropertyName>("htc", "heat transfer coefficient");
      22             : 
      23         328 :   params.addRequiredRangeCheckedParam<Real>(HeatConduction::emissivity,
      24         164 :                                             HeatConduction::emissivity + " >= 0.0 & " +
      25         164 :                                                 HeatConduction::emissivity + " <= 1.0",
      26             :                                             "emissivity of the surface");
      27             : 
      28         328 :   params.addRequiredParam<std::vector<Real>>(
      29             :       "thermal_conductivities",
      30             :       "vector of thermal conductivity values used for the conduction layers");
      31         328 :   params.addRequiredParam<std::vector<Real>>("conduction_thicknesses",
      32             :                                              "vector of conduction layer thicknesses");
      33             : 
      34         328 :   MooseEnum geometry("cartesian cylindrical", "cartesian");
      35         328 :   params.addParam<MooseEnum>("geometry", geometry, "type of geometry");
      36         328 :   params.addRangeCheckedParam<Real>("inner_radius",
      37             :                                     "inner_radius > 0.0",
      38             :                                     "coordinate corresponding to the first resistance layer");
      39             : 
      40         492 :   params.addRangeCheckedParam<Real>(
      41         328 :       "step_size", 0.1, "step_size > 0.0", "underrelaxation step size");
      42             : 
      43         492 :   params.addRangeCheckedParam<unsigned int>(
      44         328 :       "max_iterations", 100, "max_iterations >= 0", "maximum iterations");
      45             : 
      46         492 :   params.addRangeCheckedParam<Real>(
      47         328 :       "tolerance", 1E-3, "tolerance > 0.0", "tolerance to converge iterations");
      48         164 :   params.addClassDescription("Thermal resistance Heat flux boundary condition for the "
      49             :                              "fluid and solid energy equations");
      50         164 :   return params;
      51         164 : }
      52             : 
      53          88 : FVThermalResistanceBC::FVThermalResistanceBC(const InputParameters & parameters)
      54             :   : FVFluxBC(parameters),
      55          88 :     _geometry(getParam<MooseEnum>("geometry").getEnum<Moose::CoordinateSystemType>()),
      56          88 :     _inner_radius(
      57         110 :         _geometry == Moose::CoordinateSystemType::COORD_RZ ? getParam<Real>("inner_radius") : 1.0),
      58         220 :     _T(isParamValid("temperature") ? adCoupledValue("temperature") : _u),
      59          88 :     _T_ambient(getParam<Real>(HeatConduction::T_ambient)),
      60         176 :     _k(getParam<std::vector<Real>>("thermal_conductivities")),
      61         176 :     _dx(getParam<std::vector<Real>>("conduction_thicknesses")),
      62         176 :     _h(getADMaterialPropertyByName<Real>(getParam<MaterialPropertyName>("htc"))),
      63          88 :     _emissivity(getParam<Real>(HeatConduction::emissivity)),
      64         176 :     _max_iterations(getParam<unsigned int>("max_iterations")),
      65         176 :     _tolerance(getParam<Real>("tolerance")),
      66         176 :     _alpha(getParam<Real>("step_size")),
      67          88 :     _T_surface(0.0),
      68          88 :     _outer_radius(_inner_radius),
      69          88 :     _conduction_resistance(0.0),
      70         176 :     _parallel_resistance(0.0)
      71             : {
      72          88 :   if (_k.size() != _dx.size())
      73           0 :     paramError("conduction_thicknesses",
      74             :                "Number of specified thermal conductivities must match "
      75             :                "the number of conduction layers!");
      76             : 
      77          88 :   if (_geometry == Moose::CoordinateSystemType::COORD_RZ)
      78          88 :     for (const auto & d : _dx)
      79             :       _outer_radius += d;
      80             : 
      81             :   // because the thermal conductivities are constant, we only need to compute
      82             :   // the conduction resistance one time
      83          88 :   computeConductionResistance();
      84          88 : }
      85             : 
      86             : void
      87          88 : FVThermalResistanceBC::computeConductionResistance()
      88             : {
      89          88 :   Real r = _inner_radius;
      90             : 
      91         352 :   for (const auto i : index_range(_k))
      92             :   {
      93         264 :     switch (_geometry)
      94             :     {
      95         198 :       case Moose::CoordinateSystemType::COORD_XYZ:
      96         198 :         _conduction_resistance += _dx[i] / _k[i];
      97         198 :         break;
      98          66 :       case Moose::CoordinateSystemType::COORD_RZ:
      99             :       {
     100          66 :         _conduction_resistance += std::log((_dx[i] + r) / r) / _k[i];
     101          66 :         r += _dx[i];
     102          66 :         break;
     103             :       }
     104           0 :       default:
     105           0 :         mooseError("Unhandled 'GeometryEnum' in 'FVThermalResistanceBC'!");
     106             :     }
     107             :   }
     108          88 : }
     109             : 
     110             : ADReal
     111         900 : FVThermalResistanceBC::computeQpResidual()
     112             : {
     113             :   // radiation resistance has to be solved iteratively, since we don't know the
     114             :   // surface temperature. We do know that the heat flux in the conduction layers
     115             :   // must match the heat flux in the parallel convection-radiation segment. For a
     116             :   // first guess, take the surface temperature as the average of _T and T_ambient.
     117        1800 :   _T_surface = 0.5 * (_T[_qp] + _T_ambient);
     118             : 
     119             :   // total flux perpendicular to boundary
     120             :   ADReal flux;
     121             : 
     122             :   // resistance component representing the sum of the convection and radiation
     123             :   // resistances, in parallel
     124         900 :   computeParallelResistance();
     125             : 
     126             :   // other iteration requirements
     127             :   unsigned int iteration = 0;
     128         900 :   ADReal norm = 2 * _tolerance;
     129             :   ADReal T_surface_previous;
     130             : 
     131             :   // iterate to find the approximate surface temperature needed for evaluating the
     132             :   // radiation resistance. We only do this iteration if we have radiation transfer.
     133         900 :   if (_emissivity > 1e-8)
     134       54450 :     while (norm > (_tolerance * _alpha))
     135             :     {
     136       54450 :       T_surface_previous = _T_surface;
     137             : 
     138             :       // compute the flux based on the conduction part of the circuit
     139      108900 :       flux = (_T[_qp] - _T_surface) / _conduction_resistance;
     140             : 
     141       54450 :       computeParallelResistance();
     142             : 
     143             :       // use the flux computed from the conduction half to update T_surface
     144      108900 :       _T_surface = flux * _parallel_resistance + _T_ambient;
     145      163350 :       _T_surface = _alpha * _T_surface + (1 - _alpha) * T_surface_previous;
     146      163350 :       norm = std::abs(_T_surface - T_surface_previous) / std::abs(T_surface_previous);
     147             : 
     148       54450 :       if (iteration == _max_iterations)
     149             :       {
     150         450 :         mooseWarning("Maximum number of iterations reached in 'FVThermalResistanceBC'!");
     151             :         break;
     152             :       }
     153             :       else
     154       54000 :         iteration += 1;
     155             :     }
     156             : 
     157             :   // once we have determined T_surface, we can finally evaluate the complete
     158             :   // resistance to find the overall heat flux. For Cartesian, dividing by the
     159             :   // 'inner_radius' has no effect, but it is required for correct normalization
     160             :   // for cylindrical geometries.
     161        2700 :   flux = (_T[_qp] - _T_ambient) / (_conduction_resistance + _parallel_resistance) / _inner_radius;
     162         900 :   return flux;
     163             : }
     164             : 
     165             : void
     166       55350 : FVThermalResistanceBC::computeParallelResistance()
     167             : {
     168             :   // compute the parallel convection and radiation resistances, assuming they
     169             :   // act on the same surface area size
     170       55350 :   ADReal hr = _emissivity * HeatConduction::Constants::sigma *
     171       55350 :               (_T_surface * _T_surface + _T_ambient * _T_ambient) * (_T_surface + _T_ambient);
     172             : 
     173             :   // for Cartesian, dividing by the 'outer_radius' has no effect, but it is
     174             :   // required for correct normalization for cylindrical geometries
     175      166050 :   _parallel_resistance = 1.0 / (hr + _h[_qp]) / _outer_radius;
     176       55350 : }