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         249 : GrayLambertSurfaceRadiationBase::validParams()
      19             : {
      20         249 :   InputParameters params = SideUserObject::validParams();
      21         498 :   params.addParam<Real>(
      22             :       "stefan_boltzmann_constant",
      23         498 :       5.670367e-8,
      24             :       "The Stefan-Boltzmann constant. Default value is in units of [W / m^2 K^4].");
      25         498 :   params.addRequiredCoupledVar("temperature", "The coupled temperature variable.");
      26         498 :   params.addRequiredParam<std::vector<FunctionName>>("emissivity",
      27             :                                                      "Emissivities for each boundary.");
      28         498 :   params.addParam<std::vector<BoundaryName>>(
      29             :       "fixed_temperature_boundary",
      30             :       {},
      31             :       "The list of boundary IDs from the mesh with fixed temperatures.");
      32         498 :   params.addParam<std::vector<FunctionName>>(
      33             :       "fixed_boundary_temperatures", {}, "The temperatures of the fixed boundary.");
      34         498 :   params.addParam<std::vector<BoundaryName>>(
      35             :       "adiabatic_boundary", {}, "The list of boundary IDs from the mesh that are adiabatic.");
      36             : 
      37         249 :   params.addClassDescription(
      38             :       "This object implements the exchange of heat by radiation between sidesets.");
      39         249 :   return params;
      40           0 : }
      41             : 
      42         134 : GrayLambertSurfaceRadiationBase::GrayLambertSurfaceRadiationBase(const InputParameters & parameters)
      43             :   : SideUserObject(parameters),
      44         134 :     _sigma_stefan_boltzmann(getParam<Real>("stefan_boltzmann_constant")),
      45         134 :     _n_sides(boundaryIDs().size()),
      46         134 :     _temperature(coupledValue("temperature")),
      47         134 :     _radiosity(_n_sides),
      48         134 :     _heat_flux_density(_n_sides),
      49         134 :     _side_temperature(_n_sides),
      50         134 :     _side_type(_n_sides),
      51         134 :     _areas(_n_sides),
      52         134 :     _beta(_n_sides),
      53         268 :     _surface_irradiation(_n_sides)
      54             : {
      55             :   // get emissivity functions
      56         134 :   auto & eps_names = getParam<std::vector<FunctionName>>("emissivity");
      57         134 :   _emissivity.resize(eps_names.size());
      58        1012 :   for (unsigned int j = 0; j < _emissivity.size(); ++j)
      59         878 :     _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         268 :     std::vector<BoundaryName> boundary_names = getParam<std::vector<BoundaryName>>("boundary");
      65             : 
      66        1018 :     for (unsigned int j = 0; j < boundary_names.size(); ++j)
      67             :     {
      68         884 :       if (boundary_names[j] == "ANY_BOUNDARY_ID")
      69           0 :         paramError("boundary", "boundary must be explicitly provided.");
      70             : 
      71         884 :       _side_id_index[_mesh.getBoundaryID(boundary_names[j])] = j;
      72         884 :       _side_type[j] = VARIABLE_TEMPERATURE;
      73             :     }
      74             : 
      75             :     // consistency check on emissivity, must be as many entries as boundary
      76         134 :     if (boundary_names.size() != _emissivity.size())
      77           2 :       paramError("emissivity", "The number of entries must match the number of boundary entries.");
      78         132 :   }
      79             : 
      80             :   // get the fixed boundaries of the system if any are provided
      81         264 :   if (isParamValid("fixed_temperature_boundary"))
      82             :   {
      83             :     // if fixed_temperature_boundary is valid, then fixed_side_temperatures must be
      84             :     // valid, too
      85         264 :     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         396 :     auto fst_fn = getParam<std::vector<FunctionName>>("fixed_boundary_temperatures");
      90         288 :     for (auto & fn : fst_fn)
      91         156 :       _fixed_side_temperature.push_back(&getFunctionByName(fn));
      92             : 
      93             :     // get the fixed boundaries and temperatures
      94             :     std::vector<BoundaryName> boundary_names =
      95         396 :         getParam<std::vector<BoundaryName>>("fixed_temperature_boundary");
      96             : 
      97         132 :     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         282 :     for (auto & name : boundary_names)
     104             :     {
     105         152 :       _fixed_side_id_index[_mesh.getBoundaryID(name)] = index;
     106         152 :       index++;
     107             :     }
     108             : 
     109             :     // check that fixed side ids is a subset of boundary ids
     110             :     // and update _side_type info
     111         280 :     for (auto & p : _fixed_side_id_index)
     112             :     {
     113         152 :       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         150 :       _side_type[_side_id_index.find(p.first)->second] = FIXED_TEMPERATURE;
     117             :     }
     118         128 :   }
     119             : 
     120             :   // get the fixed boundaries of the system if any are provided
     121         256 :   if (isParamValid("adiabatic_boundary"))
     122             :   {
     123             :     // get the adiabatic boundaries and temperatures
     124             :     std::vector<BoundaryName> boundary_names =
     125         384 :         getParam<std::vector<BoundaryName>>("adiabatic_boundary");
     126             : 
     127         448 :     for (auto & name : boundary_names)
     128         320 :       _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         444 :     for (auto & id : _adiabatic_side_ids)
     133             :     {
     134         318 :       if (_side_id_index.find(id) == _side_id_index.end())
     135           2 :         paramError("adiabatic_boundary", "adiabatic_boundary must be a subset of boundary.");
     136         316 :       _side_type[_side_id_index.find(id)->second] = ADIABATIC;
     137             :     }
     138             : 
     139             :     // make sure that adiabatic boundaries are not already fixed boundaries
     140         442 :     for (auto & id : _adiabatic_side_ids)
     141         316 :       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         126 :   }
     144         126 : }
     145             : 
     146             : void
     147      130524 : 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      130524 :   unsigned int index = _side_id_index.find(_current_boundary_id)->second;
     152             : 
     153      458004 :   for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
     154             :   {
     155      327480 :     _areas[index] += _JxW[qp] * _coord[qp];
     156             : 
     157             :     Real temp = 0;
     158      327480 :     if (_side_type[index] == ADIABATIC)
     159      128340 :       continue;
     160      199140 :     else if (_side_type[index] == VARIABLE_TEMPERATURE)
     161      182528 :       temp = _temperature[qp];
     162       16612 :     else if (_side_type[index] == FIXED_TEMPERATURE)
     163             :     {
     164       16612 :       unsigned int iso_index = _fixed_side_id_index.find(_current_boundary_id)->second;
     165       16612 :       temp = _fixed_side_temperature[iso_index]->value(_t, _q_point[qp]);
     166             :     }
     167             : 
     168      199140 :     _beta[index] += _JxW[qp] * _coord[qp] * _sigma_stefan_boltzmann *
     169      199140 :                     _emissivity[index]->value(temp, Point()) * MathUtils::pow(temp, 4);
     170      199140 :     _side_temperature[index] += _JxW[qp] * _coord[qp] * temp;
     171             :   }
     172      130524 : }
     173             : 
     174             : void
     175        8759 : 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        8759 :   _view_factors = setViewFactors();
     180             : 
     181             :   // initialize areas, beta, side temps
     182       68754 :   for (unsigned int j = 0; j < _n_sides; ++j)
     183             :   {
     184       59995 :     _areas[j] = 0;
     185       59995 :     _beta[j] = 0;
     186       59995 :     _side_temperature[j] = 0;
     187             :   }
     188        8759 : }
     189             : 
     190             : void
     191        7641 : GrayLambertSurfaceRadiationBase::finalize()
     192             : {
     193             :   // need to do some parallel communiction here
     194        7641 :   gatherSum(_areas);
     195        7641 :   gatherSum(_beta);
     196        7641 :   gatherSum(_side_temperature);
     197             : 
     198             :   // first compute averages from the totals
     199       59746 :   for (unsigned int j = 0; j < _n_sides; ++j)
     200             :   {
     201       52105 :     _beta[j] /= _areas[j];
     202       52105 :     _side_temperature[j] /= _areas[j];
     203             :   }
     204             : 
     205             :   // matrix and rhs vector for the view factor calculation
     206        7641 :   DenseMatrix<Real> matrix(_n_sides, _n_sides);
     207        7641 :   DenseVector<Real> rhs(_n_sides);
     208        7641 :   DenseVector<Real> radiosity(_n_sides);
     209       59746 :   for (unsigned int i = 0; i < _n_sides; ++i)
     210             :   {
     211       52105 :     rhs(i) = _beta[i];
     212       52105 :     matrix(i, i) = 1;
     213      446444 :     for (unsigned int j = 0; j < _n_sides; ++j)
     214             :     {
     215      394339 :       if (_side_type[i] == ADIABATIC)
     216      125148 :         matrix(i, j) -= _view_factors[i][j];
     217             :       else
     218      269191 :         matrix(i, j) -=
     219      269191 :             (1 - _emissivity[i]->value(_side_temperature[i], Point())) * _view_factors[i][j];
     220             :     }
     221             :   }
     222             : 
     223             :   // compute the radiosityes
     224        7641 :   matrix.lu_solve(rhs, radiosity);
     225             : 
     226             :   // store the radiosity, temperatures and heat flux density for each surface
     227       59746 :   for (unsigned int i = 0; i < _n_sides; ++i)
     228             :   {
     229       52105 :     _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       52105 :     _heat_flux_density[i] = radiosity(i);
     235      446444 :     for (unsigned int j = 0; j < _n_sides; ++j)
     236      394339 :       _heat_flux_density[i] -= _view_factors[i][j] * radiosity(j);
     237             : 
     238       52105 :     if (_side_type[i] == ADIABATIC)
     239             :     {
     240       17748 :       Real e = _emissivity[i]->value(_side_temperature[i], Point());
     241       17748 :       _side_temperature[i] = std::pow(
     242       17748 :           (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       52105 :     _surface_irradiation[i] = 0;
     247      446444 :     for (unsigned int j = 0; j < _n_sides; ++j)
     248      394339 :       _surface_irradiation[i] += _view_factors[i][j] * radiosity(j);
     249             :   }
     250        7641 : }
     251             : 
     252             : void
     253        1118 : GrayLambertSurfaceRadiationBase::threadJoin(const UserObject & y)
     254             : {
     255             :   const auto & pps = static_cast<const GrayLambertSurfaceRadiationBase &>(y);
     256             : 
     257        9008 :   for (unsigned int j = 0; j < _n_sides; ++j)
     258             :   {
     259        7890 :     _areas[j] += pps._areas[j];
     260        7890 :     _side_temperature[j] += pps._side_temperature[j];
     261        7890 :     _beta[j] += pps._beta[j];
     262             :   }
     263        1118 : }
     264             : 
     265             : std::set<BoundaryID>
     266           9 : GrayLambertSurfaceRadiationBase::getSurfaceIDs() const
     267             : {
     268             :   std::set<BoundaryID> surface_ids;
     269          45 :   for (auto & p : _side_id_index)
     270          36 :     surface_ids.insert(p.first);
     271           9 :   return surface_ids;
     272             : }
     273             : 
     274             : Real
     275      617728 : GrayLambertSurfaceRadiationBase::getSurfaceIrradiation(BoundaryID id) const
     276             : {
     277      617728 :   if (_side_id_index.find(id) == _side_id_index.end())
     278             :     return 0;
     279      617728 :   return _surface_irradiation[_side_id_index.find(id)->second];
     280             : }
     281             : 
     282             : Real
     283      243963 : GrayLambertSurfaceRadiationBase::getSurfaceHeatFluxDensity(BoundaryID id) const
     284             : {
     285      243963 :   if (_side_id_index.find(id) == _side_id_index.end())
     286             :     return 0;
     287      243963 :   return _heat_flux_density[_side_id_index.find(id)->second];
     288             : }
     289             : 
     290             : Real
     291          82 : GrayLambertSurfaceRadiationBase::getSurfaceTemperature(BoundaryID id) const
     292             : {
     293          82 :   if (_side_id_index.find(id) == _side_id_index.end())
     294             :     return 0;
     295          82 :   return _side_temperature[_side_id_index.find(id)->second];
     296             : }
     297             : 
     298             : Real
     299          73 : GrayLambertSurfaceRadiationBase::getSurfaceRadiosity(BoundaryID id) const
     300             : {
     301          73 :   if (_side_id_index.find(id) == _side_id_index.end())
     302             :     return 0;
     303          73 :   return _radiosity[_side_id_index.find(id)->second];
     304             : }
     305             : 
     306             : Real
     307      979364 : GrayLambertSurfaceRadiationBase::getSurfaceEmissivity(BoundaryID id) const
     308             : {
     309             :   auto p = _side_id_index.find(id);
     310      979364 :   if (p == _side_id_index.end())
     311             :     return 1;
     312      979364 :   return _emissivity[p->second]->value(_side_temperature[p->second], Point());
     313             : }
     314             : 
     315             : Real
     316           0 : GrayLambertSurfaceRadiationBase::getViewFactor(BoundaryID from_id, BoundaryID to_id) const
     317             : {
     318           0 :   if (_side_id_index.find(from_id) == _side_id_index.end())
     319             :     return 0;
     320           0 :   if (_side_id_index.find(to_id) == _side_id_index.end())
     321             :     return 0;
     322           0 :   return _view_factors[_side_id_index.find(from_id)->second][_side_id_index.find(to_id)->second];
     323             : }