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 "HSCoupler2D3DUserObject.h"
      11             : #include "StoreVariableByElemIDSideUserObject.h"
      12             : #include "MeshAlignment2D3D.h"
      13             : #include "THMUtils.h"
      14             : #include "HeatTransferModels.h"
      15             : #include "Function.h"
      16             : 
      17             : registerMooseObject("ThermalHydraulicsApp", HSCoupler2D3DUserObject);
      18             : 
      19             : InputParameters
      20          73 : HSCoupler2D3DUserObject::validParams()
      21             : {
      22          73 :   InputParameters params = SideUserObject::validParams();
      23             : 
      24         146 :   params.addRequiredCoupledVar("temperature", "Temperature");
      25         146 :   params.addRequiredParam<Real>("radius_2d", "Radius of the 2D heat structure boundary [m]");
      26         146 :   params.addParam<FunctionName>(
      27             :       "emissivity_2d",
      28         146 :       0,
      29             :       "Emissivity of the 2D heat structure boundary as a function of temperature [K]");
      30         146 :   params.addParam<FunctionName>(
      31             :       "emissivity_3d",
      32         146 :       0,
      33             :       "Emissivity of the 3D heat structure boundary as a function of temperature [K]");
      34         146 :   params.addRequiredParam<FunctionName>("gap_thickness",
      35             :                                         "Gap thickness [m] as a function of temperature [K]");
      36         146 :   params.addRequiredParam<FunctionName>(
      37             :       "gap_thermal_conductivity",
      38             :       "Gap thermal conductivity [W/(m-K)] as a function of temperature [K]");
      39         146 :   params.addParam<FunctionName>(
      40         146 :       "gap_htc", 0, "Gap heat transfer coefficient [W/(m^2-K)] as a function of temperature [K]");
      41         146 :   params.addRequiredParam<UserObjectName>(
      42             :       "temperature_2d_uo",
      43             :       "StoreVariableByElemIDSideUserObject containing the temperature values on the 2D boundary");
      44         146 :   params.addRequiredParam<MeshAlignment2D3D *>("mesh_alignment", "Mesh alignment object");
      45         219 :   params.addRangeCheckedParam<Real>("symmetry_factor",
      46         146 :                                     1.0,
      47             :                                     "symmetry_factor>=1.0",
      48             :                                     "Azimuthal symmetry correction factor (>= 1.0). Equal to 2*pi "
      49             :                                     "divided by the azimuthal angle covered by the 3D mesh.");
      50             : 
      51          73 :   params.addClassDescription("Computes heat fluxes for HSCoupler2D3D.");
      52             : 
      53          73 :   return params;
      54           0 : }
      55             : 
      56          39 : HSCoupler2D3DUserObject::HSCoupler2D3DUserObject(const InputParameters & parameters)
      57             :   : SideUserObject(parameters),
      58             : 
      59          39 :     _T_3d(adCoupledValue("temperature")),
      60          78 :     _r_2d(getParam<Real>("radius_2d")),
      61          39 :     _emissivity_2d_fn(getFunction("emissivity_2d")),
      62          39 :     _emissivity_3d_fn(getFunction("emissivity_3d")),
      63          96 :     _include_radiation(isParamSetByUser("emissivity_2d") && isParamSetByUser("emissivity_3d")),
      64          39 :     _gap_thickness_fn(getFunction("gap_thickness")),
      65          39 :     _k_gap_fn(getFunction("gap_thermal_conductivity")),
      66          39 :     _htc_gap_fn(getFunction("gap_htc")),
      67          39 :     _temperature_2d_uo(getUserObject<StoreVariableByElemIDSideUserObject>("temperature_2d_uo")),
      68          78 :     _mesh_alignment(*getParam<MeshAlignment2D3D *>("mesh_alignment")),
      69         117 :     _symmetry_factor(getParam<Real>("symmetry_factor"))
      70             : {
      71          39 :   _mesh_alignment.buildCoupledElemQpIndexMapSecondary(_assembly);
      72          39 : }
      73             : 
      74             : void
      75        1155 : HSCoupler2D3DUserObject::initialize()
      76             : {
      77             :   _elem_id_to_heat_flux.clear();
      78        1155 : }
      79             : 
      80             : void
      81       64253 : HSCoupler2D3DUserObject::execute()
      82             : {
      83       64253 :   const auto elem_id_3d = _current_elem->id();
      84       64253 :   const auto elem_id_2d = _mesh_alignment.getCoupledPrimaryElemID(elem_id_3d);
      85             : 
      86       64253 :   const auto n_qp_2d = _mesh_alignment.getPrimaryNumberOfQuadraturePoints();
      87       64253 :   const auto n_qp_3d = _qrule->n_points();
      88             : 
      89       64253 :   const auto & T_2d_values = _temperature_2d_uo.getVariableValues(elem_id_2d);
      90       64253 :   const auto & area_2d = _mesh_alignment.getPrimaryArea(elem_id_2d);
      91             : 
      92       64253 :   std::vector<ADReal> heat_flux_3d(n_qp_3d, 0.0);
      93       64253 :   std::vector<ADReal> heat_flux_2d(n_qp_2d, 0.0);
      94      321253 :   for (unsigned int qp_3d = 0; qp_3d < n_qp_3d; qp_3d++)
      95             :   {
      96      257003 :     const auto qp_2d = _mesh_alignment.getCoupledPrimaryElemQpIndex(elem_id_3d, qp_3d);
      97             : 
      98      257003 :     const auto & T_2d = T_2d_values[qp_2d];
      99      257003 :     const auto & T_3d = _T_3d[qp_3d];
     100      514006 :     const auto T_gap = 0.5 * (T_2d + T_3d);
     101             : 
     102      257003 :     const auto gap_thickness = evaluateTemperatureFunction(_gap_thickness_fn, T_gap);
     103      257003 :     const auto k_gap = evaluateTemperatureFunction(_k_gap_fn, T_gap);
     104             : 
     105      257003 :     if (!MooseUtils::absoluteFuzzyGreaterThan(gap_thickness, 0))
     106           3 :       mooseError("Gap thickness must be > 0.");
     107             : 
     108      257000 :     const auto r_3d = _r_2d + gap_thickness;
     109             : 
     110             :     const auto heat_flux_cond =
     111      257000 :         HeatTransferModels::cylindricalGapConductionHeatFlux(k_gap, _r_2d, r_3d, T_2d, T_3d);
     112      257000 :     auto heat_flux = heat_flux_cond;
     113             : 
     114      257000 :     const auto htc = evaluateTemperatureFunction(_htc_gap_fn, T_gap);
     115      257000 :     heat_flux += htc * (T_2d - T_3d);
     116             : 
     117      257000 :     if (_include_radiation)
     118             :     {
     119      102000 :       const auto emissivity_2d = evaluateTemperatureFunction(_emissivity_2d_fn, T_2d);
     120      102000 :       const auto emissivity_3d = evaluateTemperatureFunction(_emissivity_3d_fn, T_3d);
     121             : 
     122      102000 :       heat_flux += HeatTransferModels::cylindricalGapRadiationHeatFlux(
     123      102000 :           _r_2d, r_3d, emissivity_2d, emissivity_3d, T_2d, T_3d);
     124             :     }
     125             : 
     126      257000 :     heat_flux_3d[qp_3d] = heat_flux;
     127             :     heat_flux_2d[qp_2d] -=
     128      514000 :         _symmetry_factor * _JxW[qp_3d] * _coord[qp_3d] * heat_flux / area_2d[qp_2d];
     129             :   }
     130             : 
     131             :   // Store values in maps
     132       64250 :   _elem_id_to_heat_flux[elem_id_3d] = heat_flux_3d;
     133             :   auto it = _elem_id_to_heat_flux.find(elem_id_2d);
     134       64250 :   if (it == _elem_id_to_heat_flux.end())
     135        8784 :     _elem_id_to_heat_flux[elem_id_2d] = heat_flux_2d;
     136             :   else
     137             :   {
     138       55466 :     auto & heat_flux_2d_existing = _elem_id_to_heat_flux[elem_id_2d];
     139      166398 :     for (const auto qp_2d : index_range(heat_flux_2d_existing))
     140      110932 :       heat_flux_2d_existing[qp_2d] += heat_flux_2d[qp_2d];
     141             :   }
     142       64250 : }
     143             : 
     144             : void
     145         144 : HSCoupler2D3DUserObject::threadJoin(const UserObject & uo)
     146             : {
     147             :   const auto & other_uo = static_cast<const HSCoupler2D3DUserObject &>(uo);
     148        8904 :   for (auto & it : other_uo._elem_id_to_heat_flux)
     149        8760 :     if (_elem_id_to_heat_flux.find(it.first) == _elem_id_to_heat_flux.end())
     150        8760 :       _elem_id_to_heat_flux[it.first] = it.second;
     151             :     else
     152             :     {
     153           0 :       auto & existing = _elem_id_to_heat_flux[it.first];
     154           0 :       for (const auto qp : index_range(existing))
     155           0 :         existing[qp] += it.second[qp];
     156             :     }
     157         144 : }
     158             : 
     159             : void
     160        1008 : HSCoupler2D3DUserObject::finalize()
     161             : {
     162        1008 :   THM::allGatherADVectorMapSum(comm(), _elem_id_to_heat_flux);
     163        1008 : }
     164             : 
     165             : const std::vector<ADReal> &
     166     2056000 : HSCoupler2D3DUserObject::getHeatFlux(dof_id_type elem_id) const
     167             : {
     168             :   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
     169             : 
     170             :   auto it = _elem_id_to_heat_flux.find(elem_id);
     171     2056000 :   if (it != _elem_id_to_heat_flux.end())
     172     2056000 :     return it->second;
     173             :   else
     174           0 :     mooseError(name(), ": Requested heat flux for element ", elem_id, " was not computed.");
     175             : }
     176             : 
     177             : ADReal
     178      975006 : HSCoupler2D3DUserObject::evaluateTemperatureFunction(const Function & fn, const ADReal & T) const
     179             : {
     180      975006 :   ADReal f = fn.value(T);
     181      975006 :   f.derivatives() = T.derivatives() * fn.timeDerivative(T.value());
     182             : 
     183      975006 :   return f;
     184             : }