www.mooseframework.org
GapHeatTransfer.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 "GapHeatTransfer.h"
11 
12 // MOOSE includes
13 #include "AddVariableAction.h"
14 #include "Assembly.h"
15 #include "MooseMesh.h"
16 #include "MooseVariable.h"
17 #include "PenetrationLocator.h"
18 #include "SystemBase.h"
19 
20 #include "libmesh/string_to_enum.h"
21 
22 registerMooseObject("HeatTransferApp", GapHeatTransfer);
23 
26 {
28  params.addClassDescription("Transfers heat across a gap between two "
29  "surfaces dependent on the gap geometry specified.");
30  params.addParam<std::string>(
31  "appended_property_name", "", "Name appended to material properties to make them unique");
32 
33  // Common
34  params.addParam<Real>("min_gap", 1.0e-6, "A minimum gap size");
35  params.addParam<Real>("max_gap", 1.0e6, "A maximum gap size");
36  params.addRangeCheckedParam<unsigned int>(
37  "min_gap_order", 0, "min_gap_order<=1", "Order of the Taylor expansion below min_gap");
38 
39  // Deprecated parameter
40  MooseEnum coord_types("default XYZ cyl", "default");
42  "coord_type",
43  coord_types,
44  "Gap calculation type (default or XYZ).",
45  "The functionality of this parameter is replaced by 'gap_geometry_type'.");
46 
47  MooseEnum gap_geom_types("PLATE CYLINDER SPHERE");
48  params.addParam<MooseEnum>("gap_geometry_type",
49  gap_geom_types,
50  "Gap calculation type. Choices are: " + gap_geom_types.getRawNames());
51 
52  params.addParam<RealVectorValue>("cylinder_axis_point_1",
53  "Start point for line defining cylindrical axis");
54  params.addParam<RealVectorValue>("cylinder_axis_point_2",
55  "End point for line defining cylindrical axis");
56  params.addParam<RealVectorValue>("sphere_origin", "Origin for sphere geometry");
57 
58  // Quadrature based
59  params.addParam<bool>("quadrature",
60  false,
61  "Whether or not to do Quadrature point based gap heat "
62  "transfer. If this is true then gap_distance and "
63  "gap_temp should NOT be provided (and will be "
64  "ignored) however paired_boundary IS then required.");
65  params.addParam<BoundaryName>("paired_boundary", "The boundary to be penetrated");
66 
68  params.addParam<MooseEnum>("order", orders, "The finite element order");
69 
70  params.addParam<bool>(
71  "warnings", false, "Whether to output warning messages concerning nodes not being found");
72 
73  params.addCoupledVar(
74  "displacements",
75  "The displacements appropriate for the simulation geometry and coordinate system");
76 
77  // Node based options
78  params.addCoupledVar("gap_distance", "Distance across the gap");
79  params.addCoupledVar("gap_temp", "Temperature on the other side of the gap");
80 
81  return params;
82 }
83 
85  : IntegratedBC(parameters),
86  _gap_geometry_type(declareRestartableData<GapConductance::GAP_GEOMETRY>("gap_geometry_type",
87  GapConductance::PLATE)),
88  _quadrature(getParam<bool>("quadrature")),
89  _secondary_flux(!_quadrature ? &_sys.getVector("secondary_flux") : NULL),
90  _gap_conductance(getMaterialProperty<Real>("gap_conductance" +
91  getParam<std::string>("appended_property_name"))),
92  _gap_conductance_dT(getMaterialProperty<Real>(
93  "gap_conductance" + getParam<std::string>("appended_property_name") + "_dT")),
94  _min_gap(getParam<Real>("min_gap")),
95  _min_gap_order(getParam<unsigned int>("min_gap_order")),
96  _max_gap(getParam<Real>("max_gap")),
97  _gap_temp(0),
98  _gap_distance(std::numeric_limits<Real>::max()),
99  _edge_multiplier(1.0),
100  _has_info(false),
101  _disp_vars(3, libMesh::invalid_uint),
102  _gap_distance_value(_quadrature ? _zero : coupledValue("gap_distance")),
103  _gap_temp_value(_quadrature ? _zero : coupledValue("gap_temp")),
104  _penetration_locator(
105  !_quadrature ? NULL
106  : &getQuadraturePenetrationLocator(
107  parameters.get<BoundaryName>("paired_boundary"),
108  getParam<std::vector<BoundaryName>>("boundary")[0],
109  Utility::string_to_enum<Order>(parameters.get<MooseEnum>("order")))),
110  _warnings(getParam<bool>("warnings")),
111  _p1(declareRestartableData<Point>("cylinder_axis_point_1", Point(0, 1, 0))),
112  _p2(declareRestartableData<Point>("cylinder_axis_point_2", Point(0, 0, 0))),
113  _pinfo(nullptr),
114  _secondary_side_phi(nullptr),
115  _secondary_side(nullptr),
116  _secondary_j(0)
117 {
118  if (isParamValid("displacements"))
119  {
120  // modern parameter scheme for displacements
121  for (unsigned int i = 0; i < coupledComponents("displacements"); ++i)
122  _disp_vars[i] = coupled("displacements", i);
123  }
124 
125  if (_quadrature)
126  {
127  if (!parameters.isParamValid("paired_boundary"))
128  mooseError(std::string("No 'paired_boundary' provided for ") + _name);
129  }
130  else
131  {
132  if (!isCoupled("gap_distance"))
133  mooseError(std::string("No 'gap_distance' provided for ") + _name);
134 
135  if (!isCoupled("gap_temp"))
136  mooseError(std::string("No 'gap_temp' provided for ") + _name);
137  }
138 }
139 
140 void
142 {
143  std::set<SubdomainID> subdomain_ids;
144 
145  for (const auto & bnd_elem : *_mesh.getBoundaryElementRange())
146  {
147  Elem * elem = bnd_elem->_elem;
148  subdomain_ids.insert(elem->subdomain_id());
149  }
150 
151  if (subdomain_ids.empty())
152  mooseError("No boundary elements found");
153 
154  Moose::CoordinateSystemType coord_system = _fe_problem.getCoordSystem(*subdomain_ids.begin());
155 
156  for (auto sid : subdomain_ids)
157  if (_fe_problem.getCoordSystem(sid) != coord_system)
158  mooseError("The GapHeatTransfer model requires all boundary elements to have the same "
159  "coordinate system.");
160 
163 }
164 
165 Real
167 {
169 
170  if (!_has_info)
171  return 0.0;
172 
174 
175  // This is keeping track of this residual contribution so it can be used as the flux on the other
176  // side of the gap.
177  if (!_quadrature)
178  {
179  Threads::spin_mutex::scoped_lock lock(Threads::spin_mutex);
180  const Real secondary_flux = computeSecondaryFluxContribution(grad_t);
181  _secondary_flux->add(_var.dofIndices()[_i], secondary_flux);
182  }
183 
184  return _test[_i][_qp] * grad_t;
185 }
186 
187 Real
189 {
190  return _coord[_qp] * _JxW[_qp] * _test[_i][_qp] * grad_t;
191 }
192 
193 void
195 {
197 
198  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
199  {
200  // compute this up front because it only depends on the quadrature point
202 
203  for (_i = 0; _i < _test.size(); _i++)
204  for (_j = 0; _j < _phi.size(); _j++)
206 
207  // Ok now do the contribution from the secondary side
208  if (_quadrature && _has_info)
209  {
210  std::vector<dof_id_type> secondary_side_dof_indices;
211 
212  _sys.dofMap().dof_indices(_secondary_side, secondary_side_dof_indices, _var.number());
213 
214  DenseMatrix<Number> K_secondary(_var.dofIndices().size(), secondary_side_dof_indices.size());
215 
216  mooseAssert(
217  _secondary_side_phi->size() == secondary_side_dof_indices.size(),
218  "The number of shapes does not match the number of dof indices on the secondary elem");
219 
220  for (_i = 0; _i < _test.size(); _i++)
221  for (_secondary_j = 0;
222  _secondary_j < static_cast<unsigned int>(secondary_side_dof_indices.size());
223  ++_secondary_j)
224  K_secondary(_i, _secondary_j) += _JxW[_qp] * _coord[_qp] * computeSecondaryQpJacobian();
225 
227  K_secondary,
228  _var.dofIndices(),
229  secondary_side_dof_indices,
230  _var.scalingFactor());
231  }
232  }
233 
235 
236  if (_has_diag_save_in)
237  {
238  unsigned int rows = _local_ke.m();
239  DenseVector<Number> diag(rows);
240  for (unsigned int i = 0; i < rows; i++)
241  diag(i) = _local_ke(i, i);
242 
243  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
244  for (unsigned int i = 0; i < _diag_save_in.size(); i++)
245  _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
246  }
247 }
248 
249 void
250 GapHeatTransfer::computeOffDiagJacobian(const unsigned int jvar_num)
251 {
252  if (jvar_num == _var.number())
253  {
254  computeJacobian();
255  return;
256  }
257 
258  const auto & jvar = getVariable(jvar_num);
259 
260  prepareMatrixTag(_assembly, _var.number(), jvar_num);
261 
262  auto phi_size = jvar.dofIndices().size();
263 
264  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
265  {
266  // compute this up front because it only depends on the quadrature point
268 
269  for (_i = 0; _i < _test.size(); _i++)
270  for (_j = 0; _j < phi_size; _j++)
271  _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobian(jvar_num);
272 
273  // Ok now do the contribution from the secondary side
274  if (_quadrature && _has_info)
275  {
276  std::vector<dof_id_type> secondary_side_dof_indices;
277 
278  _sys.dofMap().dof_indices(_secondary_side, secondary_side_dof_indices, jvar_num);
279 
280  DenseMatrix<Number> K_secondary(_var.dofIndices().size(), secondary_side_dof_indices.size());
281 
282  mooseAssert(
283  _secondary_side_phi->size() == secondary_side_dof_indices.size(),
284  "The number of shapes does not match the number of dof indices on the secondary elem");
285 
286  for (_i = 0; _i < _test.size(); _i++)
287  for (_secondary_j = 0;
288  _secondary_j < static_cast<unsigned int>(secondary_side_dof_indices.size());
289  ++_secondary_j)
290  K_secondary(_i, _secondary_j) +=
292 
294  K_secondary,
295  _var.dofIndices(),
296  secondary_side_dof_indices,
297  _var.scalingFactor());
298  }
299  }
300 
302 }
303 
304 Real
306 {
307  if (!_has_info)
308  return 0.0;
309 
310  return _test[_i][_qp] *
313  _phi[_j][_qp];
314 }
315 
316 Real
318 {
319  return _test[_i][_qp] *
322  (*_secondary_side_phi)[_secondary_j][0];
323 }
324 
325 Real
327 {
328  if (!_has_info)
329  return 0.0;
330 
331  unsigned int coupled_component;
332  bool active = false;
333  for (coupled_component = 0; coupled_component < _disp_vars.size(); ++coupled_component)
334  if (jvar == _disp_vars[coupled_component])
335  {
336  active = true;
337  break;
338  }
339 
340  Real dRdx = 0.0;
341  if (active)
342  {
343  // Compute dR/du_[xyz]
344  // Residual is based on
345  // h_gap = h_conduction() + h_contact() + h_radiation();
346  // grad_t = (_u[_qp] - _gap_temp) * h_gap;
347  // So we need
348  // (_u[_qp] - _gap_temp) * (dh_gap/du_[xyz]);
349  // Assuming dh_contact/du_[xyz] = dh_radiation/du_[xyz] = 0,
350  // we need dh_conduction/du_[xyz]
351  // Given
352  // h_conduction = gapK / gapLength, then
353  // dh_conduction/du_[xyz] = -gapK/gapLength^2 * dgapLength/du_[xyz]
354  // Given
355  // gapLength = ((u_x-m_x)^2+(u_y-m_y)^2+(u_z-m_z)^2)^1/2
356  // where m_[xyz] is the primary coordinate, then
357  // dGapLength/du_[xyz] =
358  // 1/2*((u_x-m_x)^2+(u_y-m_y)^2+(u_z-m_z)^2)^(-1/2)*2*(u_[xyz]-m_[xyz])
359  // = (u_[xyz]-m_[xyz])/gapLength
360  // This is the normal vector.
361 
362  const Real gapL = gapLength();
363 
364  // THIS IS NOT THE NORMAL WE NEED.
365  // WE NEED THE NORMAL FROM THE CONSTRAINT, THE NORMAL FROM THE
366  // PRIMARY SURFACE. HOWEVER, THIS IS TRICKY SINCE THE NORMAL
367  // FROM THE PRIMARY SURFACE WAS COMPUTED FOR A POINT ASSOCIATED
368  // WITH A SECONDARY NODE. NOW WE ARE AT A SECONDARY INTEGRATION POINT.
369  //
370  // HOW DO WE GET THE NORMAL WE NEED?
371  //
372  // Until we have the normal we need,
373  // we'll hope that the one we have is close to the negative of the one we need.
374  const Point & normal(_normals[_qp]);
375 
376  const Real dgap = dgapLength(-normal(coupled_component));
379  }
380  return _test[_i][_qp] * dRdx * _phi[_j][_qp];
381 }
382 
383 Real
385 {
386  if (!_has_info)
387  return 0.0;
388 
389  unsigned int coupled_component;
390  bool active = false;
391  for (coupled_component = 0; coupled_component < _disp_vars.size(); ++coupled_component)
392  if (jvar == _disp_vars[coupled_component])
393  {
394  active = true;
395  break;
396  }
397 
398  Real dRdx = 0.0;
399  if (active)
400  {
401  const Real gapL = gapLength();
402 
403  const Point & normal(_normals[_qp]);
404 
405  const Real dgap = dgapLength(-normal(coupled_component));
406 
407  // The sign of the secondary side should presumably be opposite that of the primary side
410  }
411  return _test[_i][_qp] * dRdx * (*_secondary_side_phi)[_secondary_j][0];
412 }
413 
414 Real
416 {
417  if (_has_info)
419 
420  return 1.0;
421 }
422 
423 Real
424 GapHeatTransfer::dgapLength(Real normalComponent) const
425 {
426  const Real gap_L = gapLength();
427  Real dgap = 0.0;
428 
429  if (_min_gap <= gap_L && gap_L <= _max_gap)
430  dgap = normalComponent;
431 
432  return dgap;
433 }
434 
435 void
437 {
438  if (!_quadrature)
439  {
440  _has_info = true;
443  }
444  else
445  {
448 
449  _gap_temp = 0.0;
450  _gap_distance = std::numeric_limits<Real>::max();
451  _has_info = false;
452  _edge_multiplier = 1.0;
453 
454  if (_pinfo)
455  {
457  _has_info = true;
458 
462 
463  Real tangential_tolerance = _penetration_locator->getTangentialTolerance();
464  if (tangential_tolerance != 0.0)
465  {
466  _edge_multiplier = 1.0 - _pinfo->_tangential_distance / tangential_tolerance;
467  if (_edge_multiplier < 0.0)
468  _edge_multiplier = 0.0;
469  }
470  }
471  else
472  {
473  if (_warnings)
474  mooseWarning("No gap value information found for node ",
475  qnode->id(),
476  " on processor ",
477  processor_id());
478  }
479  }
480 
483 }
MooseMesh & _mesh
const VariableTestValue & _test
unsigned int getAxisymmetricRadialCoord() const
virtual bool isCoupled(const std::string &var_name, unsigned int i=0) const
GapConductance::GAP_GEOMETRY & _gap_geometry_type
virtual void computeGapValues()
virtual unsigned int coupled(const std::string &var_name, unsigned int comp=0) const
unsigned int _j
Order
void addDeprecatedParam(const std::string &name, const T &value, const std::string &doc_string, const std::string &deprecation_message)
virtual Real computeSecondaryFluxContribution(Real grad_t)
std::vector< MooseVariableFEBase *> _diag_save_in
const unsigned int invalid_uint
const MooseArray< Point > & _normals
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
void computeJacobian() override
Real getValue(const Elem *elem, const std::vector< std::vector< OutputShape >> &phi) const
unsigned int number() const
const Elem *const & _current_elem
Real computeSecondaryQpJacobian()
compute the Jacobian contributions from the secondary side degrees of freedom
Real getTangentialTolerance()
static void computeGapRadii(const GAP_GEOMETRY gap_geometry_type, const Point &current_point, const Point &p1, const Point &p2, const Real &gap_distance, const Point &current_normal, Real &r1, Real &r2, Real &radius)
Compute current gap radii for surface integration of gas conductance.
static InputParameters validParams()
static Real gapLength(const GAP_GEOMETRY &gap_geom, const Real radius, const Real r1, const Real r2, const Real max_gap)
unsigned int _i
unsigned int m() const
Generic gap heat transfer model, with h_gap = h_conduction + h_contact + h_radiation.
static Real gapAttenuation(Real adjusted_length, Real min_gap, unsigned int min_gap_order)
The following methods are specializations for using the Parallel::packed_range_* routines for a vecto...
std::vector< std::vector< Real > > _side_phi
const MaterialProperty< Real > & _gap_conductance_dT
std::string getRawNames() const
std::map< dof_id_type, PenetrationInfo *> & _penetration_info
THREAD_ID _tid
const VariablePhiValue & _phi
DenseMatrix< Number > _local_ke
const MooseVariableFieldBase & getVariable(unsigned int jvar_num) const
void mooseWarning(Args &&... args) const
unsigned int _qp
Real _edge_multiplier
This is a factor that is used to gradually taper down the conductance if the contact point is off the...
const bool _quadrature
auto max(const L &left, const R &right)
bool isParamValid(const std::string &name) const
virtual Real computeQpJacobian() override
SystemBase & _sys
const MooseArray< Point > & _q_point
void addJacobian(Assembly &assembly, const Residuals &residuals, const Indices &dof_indices, Real scaling_factor)
Real computeSecondaryQpOffDiagJacobian(unsigned int jvar)
compute the displacement Jacobian contributions from the secondary side degrees of freedom ...
virtual DofMap & dofMap()
const MaterialProperty< Real > & _gap_conductance
const std::vector< std::vector< Real > > * _secondary_side_phi
The phi values on the secondary side.
const Real _min_gap
static MooseEnum getNonlinearVariableOrders()
const std::vector< dof_id_type > & dofIndices() const final
static void setGapGeometryParameters(const InputParameters &params, const Moose::CoordinateSystemType coord_sys, unsigned int axisymmetric_radial_coord, GAP_GEOMETRY &gap_geometry_type, Point &p1, Point &p2)
SubProblem & _subproblem
T string_to_enum(const std::string &s)
virtual Real dgapLength(Real normalComponent) const
const Real _max_gap
FEProblemBase & _fe_problem
Generic gap heat transfer model, with h_gap = h_conduction + h_contact + h_radiation.
const Elem * _secondary_side
The secondary side element (this really is a side element, e.g.
virtual Real gapLength() const
std::vector< unsigned int > _disp_vars
unsigned int number() const
const Elem * _side
MooseVariableFE< Real > * _variable
const PenetrationInfo * _pinfo
The current PenetratationInfo.
void accumulateTaggedLocalMatrix()
const MooseArray< Real > & _coord
const std::string _name
NumericVector< Number > * _secondary_flux
Node * getQuadratureNode(const Elem *elem, const unsigned short int side, const unsigned int qp)
void computeOffDiagJacobian(unsigned int jvar) override
const bool _warnings
Assembly & _assembly
void addCoupledVar(const std::string &name, const std::string &doc_string)
PenetrationLocator * _penetration_locator
MooseVariable & _var
const VariableValue & _gap_temp_value
unsigned int coupledComponents(const std::string &var_name) const
const QBase *const & _qrule
Real _tangential_distance
const unsigned int & _current_side
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
GapHeatTransfer(const InputParameters &parameters)
virtual Real computeQpResidual() override
const VariableValue & _gap_distance_value
CoordinateSystemType
registerMooseObject("HeatTransferApp", GapHeatTransfer)
void mooseError(Args &&... args) const
const InputParameters & _pars
void addClassDescription(const std::string &doc_string)
const unsigned int _min_gap_order
const InputParameters & parameters() const
Moose::CoordinateSystemType getCoordSystem(SubdomainID sid) const
void addRangeCheckedParam(const std::string &name, const T &value, const std::string &parsed_function, const std::string &doc_string)
StoredRange< MooseMesh::const_bnd_elem_iterator, const BndElement *> * getBoundaryElementRange()
void prepareMatrixTag(Assembly &assembly, unsigned int ivar, unsigned int jvar)
bool _has_diag_save_in
virtual void add(const numeric_index_type i, const Number value)=0
processor_id_type processor_id() const
virtual Real computeQpOffDiagJacobian(unsigned jvar) override
virtual Assembly & assembly(const THREAD_ID tid, const unsigned int nl_sys_num)=0
static InputParameters validParams()
virtual void initialSetup() override
void ErrorVector unsigned int
const VariableValue & _u
const Elem & get(const ElemType type_in)
unsigned int _secondary_j
The secondary side shape index.
const MooseArray< Real > & _JxW
void scalingFactor(const std::vector< Real > &factor)
bool isParamValid(const std::string &name) const