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 "GrayLambertSurfaceRadiationBase.h"
      11             : #include "MathUtils.h"
      12             : #include "Function.h"
      13             : #include "libmesh/quadrature.h"
      14             : 
      15             : #include <cmath>
      16             : 
      17             : InputParameters
      18         417 : GrayLambertSurfaceRadiationBase::validParams()
      19             : {
      20         417 :   InputParameters params = SideUserObject::validParams();
      21         834 :   params.addParam<Real>(
      22             :       "stefan_boltzmann_constant",
      23         834 :       5.670367e-8,
      24             :       "The Stefan-Boltzmann constant. Default value is in units of [W / m^2 K^4].");
      25         834 :   params.addRequiredCoupledVar("temperature", "The coupled temperature variable.");
      26         834 :   params.addRequiredParam<std::vector<FunctionName>>("emissivity",
      27             :                                                      "Emissivities for each boundary.");
      28         834 :   params.addParam<std::vector<BoundaryName>>(
      29             :       "fixed_temperature_boundary",
      30             :       {},
      31             :       "The list of boundary IDs from the mesh with fixed temperatures.");
      32         834 :   params.addParam<std::vector<FunctionName>>(
      33             :       "fixed_boundary_temperatures", {}, "The temperatures of the fixed boundary.");
      34         834 :   params.addParam<std::vector<BoundaryName>>(
      35             :       "adiabatic_boundary", {}, "The list of boundary IDs from the mesh that are adiabatic.");
      36             : 
      37         417 :   params.addClassDescription(
      38             :       "This object implements the exchange of heat by radiation between sidesets.");
      39         417 :   return params;
      40           0 : }
      41             : 
      42         224 : GrayLambertSurfaceRadiationBase::GrayLambertSurfaceRadiationBase(const InputParameters & parameters)
      43             :   : SideUserObject(parameters),
      44         224 :     _sigma_stefan_boltzmann(getParam<Real>("stefan_boltzmann_constant")),
      45         224 :     _n_sides(boundaryIDs().size()),
      46         224 :     _temperature(coupledValue("temperature")),
      47         224 :     _radiosity(_n_sides),
      48         224 :     _heat_flux_density(_n_sides),
      49         224 :     _side_temperature(_n_sides),
      50         224 :     _side_type(_n_sides),
      51         224 :     _areas(_n_sides),
      52         224 :     _beta(_n_sides),
      53         448 :     _surface_irradiation(_n_sides)
      54             : {
      55             :   // get emissivity functions
      56         224 :   auto & eps_names = getParam<std::vector<FunctionName>>("emissivity");
      57         224 :   _emissivity.resize(eps_names.size());
      58        1636 :   for (unsigned int j = 0; j < _emissivity.size(); ++j)
      59        1412 :     _emissivity[j] = &getFunctionByName(eps_names[j]);
      60             : 
      61             :   // set up the map from the side id to the local index & check
      62             :   // note that boundaryIDs is not in the right order anymore!
      63             :   {
      64         448 :     std::vector<BoundaryName> boundary_names = getParam<std::vector<BoundaryName>>("boundary");
      65             : 
      66        1642 :     for (unsigned int j = 0; j < boundary_names.size(); ++j)
      67             :     {
      68        1418 :       if (boundary_names[j] == "ANY_BOUNDARY_ID")
      69           0 :         paramError("boundary", "boundary must be explicitly provided.");
      70             : 
      71        1418 :       _side_id_index[_mesh.getBoundaryID(boundary_names[j])] = j;
      72        1418 :       _side_type[j] = VARIABLE_TEMPERATURE;
      73             :     }
      74             : 
      75             :     // consistency check on emissivity, must be as many entries as boundary
      76         224 :     if (boundary_names.size() != _emissivity.size())
      77           2 :       paramError("emissivity", "The number of entries must match the number of boundary entries.");
      78         222 :   }
      79             : 
      80             :   // get the fixed boundaries of the system if any are provided
      81         444 :   if (isParamValid("fixed_temperature_boundary"))
      82             :   {
      83             :     // if fixed_temperature_boundary is valid, then fixed_side_temperatures must be
      84             :     // valid, too
      85         444 :     if (!isParamValid("fixed_boundary_temperatures"))
      86           0 :       paramError("fixed_boundary_temperatures",
      87             :                  "fixed_temperature_boundary is provided, but fixed_boundary_temperatures is not.");
      88             : 
      89         666 :     auto fst_fn = getParam<std::vector<FunctionName>>("fixed_boundary_temperatures");
      90         475 :     for (auto & fn : fst_fn)
      91         253 :       _fixed_side_temperature.push_back(&getFunctionByName(fn));
      92             : 
      93             :     // get the fixed boundaries and temperatures
      94             :     std::vector<BoundaryName> boundary_names =
      95         666 :         getParam<std::vector<BoundaryName>>("fixed_temperature_boundary");
      96             : 
      97         222 :     if (boundary_names.size() != _fixed_side_temperature.size())
      98           2 :       paramError(
      99             :           "fixed_boundary_temperatures",
     100             :           "fixed_boundary_temperatures and fixed_temperature_boundary must have the same length.");
     101             : 
     102             :     unsigned int index = 0;
     103         469 :     for (auto & name : boundary_names)
     104             :     {
     105         249 :       _fixed_side_id_index[_mesh.getBoundaryID(name)] = index;
     106         249 :       index++;
     107             :     }
     108             : 
     109             :     // check that fixed side ids is a subset of boundary ids
     110             :     // and update _side_type info
     111         467 :     for (auto & p : _fixed_side_id_index)
     112             :     {
     113         249 :       if (_side_id_index.find(p.first) == _side_id_index.end())
     114           2 :         paramError("fixed_temperature_boundary",
     115             :                    "fixed_temperature_boundary must be a subset of boundary.");
     116         247 :       _side_type[_side_id_index.find(p.first)->second] = FIXED_TEMPERATURE;
     117             :     }
     118         218 :   }
     119             : 
     120             :   // get the fixed boundaries of the system if any are provided
     121         436 :   if (isParamValid("adiabatic_boundary"))
     122             :   {
     123             :     // get the adiabatic boundaries and temperatures
     124             :     std::vector<BoundaryName> boundary_names =
     125         654 :         getParam<std::vector<BoundaryName>>("adiabatic_boundary");
     126             : 
     127         747 :     for (auto & name : boundary_names)
     128         529 :       _adiabatic_side_ids.insert(_mesh.getBoundaryID(name));
     129             : 
     130             :     // check that adiabatic side ids is a subset of boundary ids
     131             :     // and update _side_type info
     132         743 :     for (auto & id : _adiabatic_side_ids)
     133             :     {
     134         527 :       if (_side_id_index.find(id) == _side_id_index.end())
     135           2 :         paramError("adiabatic_boundary", "adiabatic_boundary must be a subset of boundary.");
     136         525 :       _side_type[_side_id_index.find(id)->second] = ADIABATIC;
     137             :     }
     138             : 
     139             :     // make sure that adiabatic boundaries are not already fixed boundaries
     140         741 :     for (auto & id : _adiabatic_side_ids)
     141         525 :       if (_fixed_side_id_index.find(id) != _fixed_side_id_index.end())
     142           0 :         paramError("adiabatic_boundary", "Isothermal boundary cannot also be adiabatic boundary.");
     143         216 :   }
     144         216 : }
     145             : 
     146             : void
     147      172776 : GrayLambertSurfaceRadiationBase::execute()
     148             : {
     149             :   mooseAssert(_side_id_index.find(_current_boundary_id) != _side_id_index.end(),
     150             :               "Current boundary id not in _side_id_index.");
     151      172776 :   unsigned int index = _side_id_index.find(_current_boundary_id)->second;
     152             : 
     153      597192 :   for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
     154             :   {
     155      424416 :     _areas[index] += _JxW[qp] * _coord[qp];
     156             : 
     157             :     Real temp = 0;
     158      424416 :     if (_side_type[index] == ADIABATIC)
     159      158484 :       continue;
     160      265932 :     else if (_side_type[index] == VARIABLE_TEMPERATURE)
     161      245380 :       temp = _temperature[qp];
     162       20552 :     else if (_side_type[index] == FIXED_TEMPERATURE)
     163             :     {
     164       20552 :       unsigned int iso_index = _fixed_side_id_index.find(_current_boundary_id)->second;
     165       20552 :       temp = _fixed_side_temperature[iso_index]->value(_t, _q_point[qp]);
     166             :     }
     167             : 
     168      265932 :     _beta[index] += _JxW[qp] * _coord[qp] * _sigma_stefan_boltzmann *
     169      265932 :                     _emissivity[index]->value(temp, Point()) * MathUtils::pow(temp, 4);
     170      265932 :     _side_temperature[index] += _JxW[qp] * _coord[qp] * temp;
     171             :   }
     172      172776 : }
     173             : 
     174             : void
     175       13250 : GrayLambertSurfaceRadiationBase::initialize()
     176             : {
     177             :   // view factors are obtained here to make sure that another object had
     178             :   // time to compute them on exec initial
     179       13250 :   _view_factors = setViewFactors();
     180             : 
     181             :   // initialize areas, beta, side temps
     182       99311 :   for (unsigned int j = 0; j < _n_sides; ++j)
     183             :   {
     184       86065 :     _areas[j] = 0;
     185       86065 :     _beta[j] = 0;
     186       86065 :     _side_temperature[j] = 0;
     187             :   }
     188       13246 : }
     189             : 
     190             : void
     191       11490 : GrayLambertSurfaceRadiationBase::finalize()
     192             : {
     193             :   // need to do some parallel communiction here
     194       11490 :   gatherSum(_areas);
     195       11490 :   gatherSum(_beta);
     196       11490 :   gatherSum(_side_temperature);
     197             : 
     198             :   // first compute averages from the totals
     199       87113 :   for (unsigned int j = 0; j < _n_sides; ++j)
     200             :   {
     201       75623 :     _beta[j] /= _areas[j];
     202       75623 :     _side_temperature[j] /= _areas[j];
     203             :   }
     204             : 
     205             :   // matrix and rhs vector for the view factor calculation
     206       11490 :   DenseMatrix<Real> matrix(_n_sides, _n_sides);
     207       11490 :   DenseVector<Real> rhs(_n_sides);
     208       11490 :   DenseVector<Real> radiosity(_n_sides);
     209       87113 :   for (unsigned int i = 0; i < _n_sides; ++i)
     210             :   {
     211       75623 :     rhs(i) = _beta[i];
     212       75623 :     matrix(i, i) = 1;
     213      634488 :     for (unsigned int j = 0; j < _n_sides; ++j)
     214             :     {
     215      558865 :       if (_side_type[i] == ADIABATIC)
     216      180388 :         matrix(i, j) -= _view_factors[i][j];
     217             :       else
     218      378477 :         matrix(i, j) -=
     219      378477 :             (1 - _emissivity[i]->value(_side_temperature[i], Point())) * _view_factors[i][j];
     220             :     }
     221             :   }
     222             : 
     223             :   // compute the radiosityes
     224       11490 :   matrix.lu_solve(rhs, radiosity);
     225             : 
     226             :   // store the radiosity, temperatures and heat flux density for each surface
     227       87113 :   for (unsigned int i = 0; i < _n_sides; ++i)
     228             :   {
     229       75623 :     _radiosity[i] = radiosity(i);
     230             : 
     231             :     // _heat_flux_density is obtained from a somewhat cumbersome relation
     232             :     // but it has the advantage that we do not divide by 1 - emissivity
     233             :     // which blows up for black bodies
     234       75623 :     _heat_flux_density[i] = radiosity(i);
     235      634488 :     for (unsigned int j = 0; j < _n_sides; ++j)
     236      558865 :       _heat_flux_density[i] -= _view_factors[i][j] * radiosity(j);
     237             : 
     238       75623 :     if (_side_type[i] == ADIABATIC)
     239             :     {
     240       26394 :       Real e = _emissivity[i]->value(_side_temperature[i], Point());
     241       26394 :       _side_temperature[i] = std::pow(
     242       26394 :           (radiosity(i) + (1 - e) / e * _heat_flux_density[i]) / _sigma_stefan_boltzmann, 0.25);
     243             :     }
     244             : 
     245             :     // compute the surface irradiation into i from the radiosities
     246       75623 :     _surface_irradiation[i] = 0;
     247      634488 :     for (unsigned int j = 0; j < _n_sides; ++j)
     248      558865 :       _surface_irradiation[i] += _view_factors[i][j] * radiosity(j);
     249             :   }
     250       11490 : }
     251             : 
     252             : void
     253        1756 : GrayLambertSurfaceRadiationBase::threadJoin(const UserObject & y)
     254             : {
     255             :   const auto & pps = static_cast<const GrayLambertSurfaceRadiationBase &>(y);
     256             : 
     257       12198 :   for (unsigned int j = 0; j < _n_sides; ++j)
     258             :   {
     259       10442 :     _areas[j] += pps._areas[j];
     260       10442 :     _side_temperature[j] += pps._side_temperature[j];
     261       10442 :     _beta[j] += pps._beta[j];
     262             :   }
     263        1756 : }
     264             : 
     265             : std::set<BoundaryID>
     266          30 : GrayLambertSurfaceRadiationBase::getSurfaceIDs() const
     267             : {
     268             :   std::set<BoundaryID> surface_ids;
     269         150 :   for (auto & p : _side_id_index)
     270         120 :     surface_ids.insert(p.first);
     271          30 :   return surface_ids;
     272             : }
     273             : 
     274             : Real
     275      768800 : GrayLambertSurfaceRadiationBase::getSurfaceIrradiation(BoundaryID id) const
     276             : {
     277      768800 :   if (_side_id_index.find(id) == _side_id_index.end())
     278             :     return 0;
     279      768800 :   return _surface_irradiation[_side_id_index.find(id)->second];
     280             : }
     281             : 
     282             : Real
     283      367725 : GrayLambertSurfaceRadiationBase::getSurfaceHeatFluxDensity(BoundaryID id) const
     284             : {
     285      367725 :   if (_side_id_index.find(id) == _side_id_index.end())
     286             :     return 0;
     287      367725 :   return _heat_flux_density[_side_id_index.find(id)->second];
     288             : }
     289             : 
     290             : Real
     291         126 : GrayLambertSurfaceRadiationBase::getSurfaceTemperature(BoundaryID id) const
     292             : {
     293         126 :   if (_side_id_index.find(id) == _side_id_index.end())
     294             :     return 0;
     295         126 :   return _side_temperature[_side_id_index.find(id)->second];
     296             : }
     297             : 
     298             : Real
     299         111 : GrayLambertSurfaceRadiationBase::getSurfaceRadiosity(BoundaryID id) const
     300             : {
     301         111 :   if (_side_id_index.find(id) == _side_id_index.end())
     302             :     return 0;
     303         111 :   return _radiosity[_side_id_index.find(id)->second];
     304             : }
     305             : 
     306             : Real
     307     1220060 : GrayLambertSurfaceRadiationBase::getSurfaceEmissivity(BoundaryID id) const
     308             : {
     309             :   auto p = _side_id_index.find(id);
     310     1220060 :   if (p == _side_id_index.end())
     311             :     return 1;
     312     1220060 :   return _emissivity[p->second]->value(_side_temperature[p->second], Point());
     313             : }
     314             : 
     315             : Real
     316         240 : GrayLambertSurfaceRadiationBase::getViewFactor(BoundaryID from_id, BoundaryID to_id) const
     317             : {
     318         240 :   if (_side_id_index.find(from_id) == _side_id_index.end())
     319             :     return 0;
     320         240 :   if (_side_id_index.find(to_id) == _side_id_index.end())
     321             :     return 0;
     322         240 :   return _view_factors[_side_id_index.find(from_id)->second][_side_id_index.find(to_id)->second];
     323             : }