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 "HSCoupler2D2DRadiationUserObject.h"
      11             : #include "StoreVariableByElemIDSideUserObject.h"
      12             : #include "MeshAlignment2D2D.h"
      13             : #include "THMUtils.h"
      14             : #include "HeatConductionNames.h"
      15             : 
      16             : registerMooseObject("ThermalHydraulicsApp", HSCoupler2D2DRadiationUserObject);
      17             : 
      18             : InputParameters
      19         117 : HSCoupler2D2DRadiationUserObject::validParams()
      20             : {
      21         117 :   InputParameters params = SideUserObject::validParams();
      22             : 
      23         234 :   params.addRequiredParam<std::vector<Real>>("emissivities", "Emissivities of each surface");
      24         234 :   params.addRequiredParam<std::vector<std::vector<Real>>>(
      25             :       "view_factors", "View factors between each surface, as a matrix");
      26         234 :   params.addRequiredParam<bool>(
      27             :       "include_environment",
      28             :       "Whether or not to include an environment surrounding all of the surfaces.");
      29         234 :   params.addParam<Real>("T_environment", 0.0, "Environment temperature.");
      30         234 :   params.addRequiredParam<UserObjectName>(
      31             :       "temperature_uo",
      32             :       "StoreVariableByElemIDSideUserObjects containing the temperature values for each element");
      33         234 :   params.addRequiredParam<MeshAlignment2D2D *>("mesh_alignment", "Mesh alignment object");
      34             : 
      35         117 :   params.addClassDescription("Computes heat fluxes for HSCoupler2D2D.");
      36             : 
      37         117 :   return params;
      38           0 : }
      39             : 
      40          63 : HSCoupler2D2DRadiationUserObject::HSCoupler2D2DRadiationUserObject(
      41          63 :     const InputParameters & parameters)
      42             :   : SideUserObject(parameters),
      43             : 
      44          63 :     _emissivities(getParam<std::vector<Real>>("emissivities")),
      45         126 :     _view_factors(getParam<std::vector<std::vector<Real>>>("view_factors")),
      46         126 :     _include_environment(getParam<bool>("include_environment")),
      47         126 :     _T_environment(getParam<Real>("T_environment")),
      48          63 :     _temperature_uo(getUserObject<StoreVariableByElemIDSideUserObject>("temperature_uo")),
      49         126 :     _mesh_alignment(*getParam<MeshAlignment2D2D *>("mesh_alignment")),
      50          63 :     _n_hs(_emissivities.size()),
      51         126 :     _n_surfaces(_include_environment ? _n_hs + 1 : _n_hs)
      52             : {
      53          63 : }
      54             : 
      55             : void
      56        2736 : HSCoupler2D2DRadiationUserObject::initialize()
      57             : {
      58             :   _elem_id_to_heat_flux.clear();
      59        2736 : }
      60             : 
      61             : void
      62       13680 : HSCoupler2D2DRadiationUserObject::execute()
      63             : {
      64             :   // store element IDs corresponding to this side
      65             :   std::vector<dof_id_type> elem_ids;
      66       13680 :   elem_ids.push_back(_current_elem->id());
      67       13680 :   const auto coupled_elem_ids = _mesh_alignment.getCoupledSecondaryElemIDs(elem_ids[0]);
      68       13680 :   elem_ids.insert(elem_ids.end(), coupled_elem_ids.begin(), coupled_elem_ids.end());
      69             : 
      70             :   // store temperatures
      71             :   std::vector<std::vector<ADReal>> temperatures;
      72       50480 :   for (const auto & elem_id : elem_ids)
      73       36800 :     temperatures.push_back(_temperature_uo.getVariableValues(elem_id));
      74             : 
      75             :   // loop over the quadrature points and compute the heat fluxes
      76       13680 :   const auto n_qp = _qrule->n_points();
      77       13680 :   std::vector<std::vector<ADReal>> heat_fluxes(_n_hs, std::vector<ADReal>(n_qp, 0.0));
      78       41040 :   for (unsigned int qp = 0; qp < n_qp; qp++)
      79             :   {
      80             :     // form matrix and RHS
      81       27360 :     DenseMatrix<ADReal> matrix(_n_surfaces, _n_surfaces);
      82       27360 :     DenseVector<ADReal> emittances(_n_surfaces);
      83      110720 :     for (unsigned int i = 0; i < _n_surfaces; ++i)
      84             :     {
      85       83360 :       if (_include_environment && i == _n_surfaces - 1) // environment surface
      86             :       {
      87        9760 :         emittances(i) = HeatConduction::Constants::sigma * std::pow(_T_environment, 4);
      88        9760 :         matrix(i, i) = 1.0;
      89             :       }
      90             :       else // heat structure surface
      91             :       {
      92             :         emittances(i) =
      93      147200 :             _emissivities[i] * HeatConduction::Constants::sigma * std::pow(temperatures[i][qp], 4);
      94       73600 :         matrix(i, i) = 1.0;
      95      315840 :         for (unsigned int j = 0; j < _n_surfaces; ++j)
      96      242240 :           matrix(i, j) -= (1 - _emissivities[i]) * _view_factors[i][j];
      97             :       }
      98             :     }
      99             : 
     100             :     // solve for radiosities
     101       27360 :     DenseVector<ADReal> radiosities(_n_surfaces);
     102       27360 :     matrix.lu_solve(emittances, radiosities);
     103             : 
     104             :     // compute heat fluxes
     105      100960 :     for (unsigned int i = 0; i < _n_hs; ++i)
     106       73600 :       heat_fluxes[i][qp] =
     107      220800 :           1.0 / (1.0 - _emissivities[i]) * (emittances(i) - _emissivities[i] * radiosities(i));
     108       27360 :   }
     109             : 
     110             :   // store the heat fluxes by element ID
     111       50480 :   for (unsigned int i = 0; i < _n_hs; ++i)
     112       36800 :     _elem_id_to_heat_flux[elem_ids[i]] = heat_fluxes[i];
     113       27360 : }
     114             : 
     115             : void
     116         513 : HSCoupler2D2DRadiationUserObject::threadJoin(const UserObject & uo)
     117             : {
     118             :   const auto & other_uo = static_cast<const HSCoupler2D2DRadiationUserObject &>(uo);
     119         513 :   for (auto & it : other_uo._elem_id_to_heat_flux)
     120           0 :     if (_elem_id_to_heat_flux.find(it.first) == _elem_id_to_heat_flux.end())
     121           0 :       _elem_id_to_heat_flux[it.first] = it.second;
     122             :     else
     123             :     {
     124           0 :       auto & existing = _elem_id_to_heat_flux[it.first];
     125           0 :       for (const auto qp : index_range(existing))
     126           0 :         existing[qp] += it.second[qp];
     127             :     }
     128         513 : }
     129             : 
     130             : void
     131        2223 : HSCoupler2D2DRadiationUserObject::finalize()
     132             : {
     133        2223 :   THM::allGatherADVectorMapSum(comm(), _elem_id_to_heat_flux);
     134        2223 : }
     135             : 
     136             : const std::vector<ADReal> &
     137      289280 : HSCoupler2D2DRadiationUserObject::getHeatFlux(dof_id_type elem_id) const
     138             : {
     139             :   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
     140             : 
     141             :   auto it = _elem_id_to_heat_flux.find(elem_id);
     142      289280 :   if (it != _elem_id_to_heat_flux.end())
     143      289280 :     return it->second;
     144             :   else
     145           0 :     mooseError(name(), ": Requested heat flux for element ", elem_id, " was not computed.");
     146             : }