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