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