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