www.mooseframework.org
GapConductance.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
8 #include "GapConductance.h"
9 
10 // MOOSE includes
11 #include "Function.h"
12 #include "MooseMesh.h"
13 #include "MooseVariable.h"
14 #include "PenetrationLocator.h"
15 #include "SystemBase.h"
16 #include "AddVariableAction.h"
17 
18 #include "libmesh/string_to_enum.h"
19 
20 template <>
21 InputParameters
23 {
24  InputParameters params = validParams<Material>();
26 
27  params.addRequiredCoupledVar("variable", "Temperature variable");
28 
29  // Node based
30  params.addCoupledVar("gap_distance", "Distance across the gap");
31  params.addCoupledVar("gap_temp", "Temperature on the other side of the gap");
32  params.addParam<Real>("gap_conductivity", 1.0, "The thermal conductivity of the gap material");
33  params.addParam<FunctionName>(
34  "gap_conductivity_function",
35  "Thermal conductivity of the gap material as a function. Multiplied by gap_conductivity.");
36  params.addCoupledVar("gap_conductivity_function_variable",
37  "Variable to be used in the gap_conductivity_function in place of time");
38 
39  // Quadrature based
40  params.addParam<bool>("quadrature",
41  false,
42  "Whether or not to do quadrature point based gap heat "
43  "transfer. If this is true then gap_distance and "
44  "gap_temp should NOT be provided (and will be "
45  "ignored); however, paired_boundary and variable are "
46  "then required.");
47  params.addParam<BoundaryName>("paired_boundary", "The boundary to be penetrated");
48 
49  params.addParam<Real>("stefan_boltzmann", 5.669e-8, "The Stefan-Boltzmann constant");
50 
51  params.addParam<bool>("use_displaced_mesh",
52  true,
53  "Whether or not this object should use the "
54  "displaced mesh for computation. Note that in "
55  "the case this is true but no displacements "
56  "are provided in the Mesh block the "
57  "undisplaced mesh will still be used.");
58 
59  return params;
60 }
61 
62 InputParameters
64 {
65  InputParameters params = emptyInputParameters();
66  params.addParam<std::string>(
67  "appended_property_name", "", "Name appended to material properties to make them unique");
68  MooseEnum gap_geom_types("PLATE CYLINDER SPHERE");
69  params.addParam<MooseEnum>("gap_geometry_type", gap_geom_types, "Gap calculation type.");
70 
71  params.addParam<RealVectorValue>("cylinder_axis_point_1",
72  "Start point for line defining cylindrical axis");
73  params.addParam<RealVectorValue>("cylinder_axis_point_2",
74  "End point for line defining cylindrical axis");
75  params.addParam<RealVectorValue>("sphere_origin", "Origin for sphere geometry");
76 
77  params.addRangeCheckedParam<Real>("emissivity_1",
78  0.0,
79  "emissivity_1>=0 & emissivity_1<=1",
80  "The emissivity of the fuel surface");
81  params.addRangeCheckedParam<Real>("emissivity_2",
82  0.0,
83  "emissivity_2>=0 & emissivity_2<=1",
84  "The emissivity of the cladding surface");
85 
86  params.addParam<bool>(
87  "warnings", false, "Whether to output warning messages concerning nodes not being found");
88 
89  MooseEnum orders(AddVariableAction::getNonlinearVariableOrders());
90  params.addParam<MooseEnum>("order", orders, "The finite element order");
91 
92  // Common
93  params.addRangeCheckedParam<Real>(
94  "min_gap", 1e-6, "min_gap>=0", "A minimum gap (denominator) size");
95  params.addRangeCheckedParam<Real>(
96  "max_gap", 1e6, "max_gap>=0", "A maximum gap (denominator) size");
97 
98  return params;
99 }
100 
101 GapConductance::GapConductance(const InputParameters & parameters)
102  : Material(parameters),
103  _appended_property_name(getParam<std::string>("appended_property_name")),
104  _temp(coupledValue("variable")),
105  _gap_geometry_type(declareRestartableData<GapConductance::GAP_GEOMETRY>("gap_geometry_type",
107  _quadrature(getParam<bool>("quadrature")),
108  _gap_temp(0),
109  _gap_distance(88888),
110  _radius(0),
111  _r1(0),
112  _r2(0),
113  _has_info(false),
114  _gap_distance_value(_quadrature ? _zero : coupledValue("gap_distance")),
115  _gap_temp_value(_quadrature ? _zero : coupledValue("gap_temp")),
116  _gap_conductance(declareProperty<Real>("gap_conductance" + _appended_property_name)),
117  _gap_conductance_dT(declareProperty<Real>("gap_conductance" + _appended_property_name + "_dT")),
118  _gap_thermal_conductivity(declareProperty<Real>("gap_conductivity")),
119  _gap_conductivity(getParam<Real>("gap_conductivity")),
120  _gap_conductivity_function(isParamValid("gap_conductivity_function")
121  ? &getFunction("gap_conductivity_function")
122  : NULL),
123  _gap_conductivity_function_variable(isCoupled("gap_conductivity_function_variable")
124  ? &coupledValue("gap_conductivity_function_variable")
125  : NULL),
126  _stefan_boltzmann(getParam<Real>("stefan_boltzmann")),
127  _emissivity(getParam<Real>("emissivity_1") != 0.0 && getParam<Real>("emissivity_2") != 0.0
128  ? 1.0 / getParam<Real>("emissivity_1") + 1.0 / getParam<Real>("emissivity_2") -
129  1
130  : 0.0),
131  _min_gap(getParam<Real>("min_gap")),
132  _max_gap(getParam<Real>("max_gap")),
133  _temp_var(_quadrature ? getVar("variable", 0) : NULL),
134  _penetration_locator(NULL),
135  _serialized_solution(_quadrature ? &_temp_var->sys().currentSolution() : NULL),
136  _dof_map(_quadrature ? &_temp_var->sys().dofMap() : NULL),
137  _warnings(getParam<bool>("warnings")),
138  _p1(declareRestartableData<Point>("cylinder_axis_point_1", Point(0, 1, 0))),
139  _p2(declareRestartableData<Point>("cylinder_axis_point_2", Point(0, 0, 0)))
140 {
141  if (_quadrature)
142  {
143  if (!parameters.isParamValid("paired_boundary"))
144  mooseError(std::string("No 'paired_boundary' provided for ") + _name);
145  }
146  else
147  {
148  if (!isCoupled("gap_distance"))
149  mooseError(std::string("No 'gap_distance' provided for ") + _name);
150 
151  if (!isCoupled("gap_temp"))
152  mooseError(std::string("No 'gap_temp' provided for ") + _name);
153  }
154 
155  if (_quadrature)
156  {
157  _penetration_locator = &_subproblem.geomSearchData().getQuadraturePenetrationLocator(
158  parameters.get<BoundaryName>("paired_boundary"),
159  getParam<std::vector<BoundaryName>>("boundary")[0],
160  Utility::string_to_enum<Order>(parameters.get<MooseEnum>("order")));
161  }
162 }
163 
164 void
166 {
167  setGapGeometryParameters(_pars, _coord_sys, _gap_geometry_type, _p1, _p2);
168 }
169 
170 void
171 GapConductance::setGapGeometryParameters(const InputParameters & params,
172  const Moose::CoordinateSystemType coord_sys,
173  GAP_GEOMETRY & gap_geometry_type,
174  Point & p1,
175  Point & p2)
176 {
177  if (params.isParamSetByUser("gap_geometry_type"))
178  {
179  gap_geometry_type =
180  GapConductance::GAP_GEOMETRY(int(params.get<MooseEnum>("gap_geometry_type")));
181  }
182  else
183  {
184  if (coord_sys == Moose::COORD_XYZ)
185  gap_geometry_type = GapConductance::PLATE;
186  else if (coord_sys == Moose::COORD_RZ)
187  gap_geometry_type = GapConductance::CYLINDER;
188  else if (coord_sys == Moose::COORD_RSPHERICAL)
189  gap_geometry_type = GapConductance::SPHERE;
190  }
191 
192  if (gap_geometry_type == GapConductance::PLATE)
193  {
194  if (coord_sys == Moose::COORD_RSPHERICAL)
195  ::mooseError("'gap_geometry_type = PLATE' cannot be used with models having a spherical "
196  "coordinate system.");
197  }
198  else if (gap_geometry_type == GapConductance::CYLINDER)
199  {
200  if (coord_sys == Moose::COORD_XYZ)
201  {
202  if (!params.isParamValid("cylinder_axis_point_1") ||
203  !params.isParamValid("cylinder_axis_point_2"))
204  ::mooseError("For 'gap_geometry_type = CYLINDER' to be used with a Cartesian model, "
205  "'cylinder_axis_point_1' and 'cylinder_axis_point_2' must be specified.");
206  p1 = params.get<RealVectorValue>("cylinder_axis_point_1");
207  p2 = params.get<RealVectorValue>("cylinder_axis_point_2");
208  }
209  else if (coord_sys == Moose::COORD_RZ)
210  {
211  if (params.isParamValid("cylinder_axis_point_1") ||
212  params.isParamValid("cylinder_axis_point_2"))
213  ::mooseError("The 'cylinder_axis_point_1' and 'cylinder_axis_point_2' cannot be specified "
214  "with axisymmetric models. The y-axis is used as the cylindrical axis of "
215  "symmetry.");
216  p1 = Point(0, 0, 0);
217  p2 = Point(0, 1, 0);
218  }
219  else if (coord_sys == Moose::COORD_RSPHERICAL)
220  ::mooseError("'gap_geometry_type = CYLINDER' cannot be used with models having a spherical "
221  "coordinate system.");
222  }
223  else if (gap_geometry_type == GapConductance::SPHERE)
224  {
225  if (coord_sys == Moose::COORD_XYZ || coord_sys == Moose::COORD_RZ)
226  {
227  if (!params.isParamValid("sphere_origin"))
228  ::mooseError("For 'gap_geometry_type = SPHERE' to be used with a Cartesian or axisymmetric "
229  "model, 'sphere_origin' must be specified.");
230  p1 = params.get<RealVectorValue>("sphere_origin");
231  }
232  else if (coord_sys == Moose::COORD_RSPHERICAL)
233  {
234  if (params.isParamValid("sphere_origin"))
235  ::mooseError("The 'sphere_origin' cannot be specified with spherical models. x=0 is used "
236  "as the spherical origin.");
237  p1 = Point(0, 0, 0);
238  }
239  }
240 }
241 
242 void
244 {
247 }
248 
249 void
251 {
252  if (_has_info)
253  {
256  }
257  else
258  {
259  _gap_conductance[_qp] = 0;
260  _gap_conductance_dT[_qp] = 0;
261  }
262 }
263 
264 Real
266 {
268  return _gap_thermal_conductivity[_qp] /
270 }
271 
272 Real
274 {
275  return 0.0;
276 }
277 
278 Real
280 {
281  /*
282  Gap conductance due to radiation is based on the diffusion approximation:
283 
284  qr = sigma*Fe*(Tf^4 - Tc^4) ~ hr(Tf - Tc)
285  where sigma is the Stefan-Boltztmann constant, Fe is an emissivity function, Tf and Tc
286  are the fuel and clad absolute temperatures, respectively, and hr is the radiant gap
287  conductance. Solving for hr,
288 
289  hr = sigma*Fe*(Tf^4 - Tc^4) / (Tf - Tc)
290  which can be factored to give:
291 
292  hr = sigma*Fe*(Tf^2 + Tc^2) * (Tf + Tc)
293 
294  Approximating the fuel-clad gap as infinite parallel planes, the emissivity function is given by:
295 
296  Fe = 1 / (1/ef + 1/ec - 1)
297  */
298 
299  if (_emissivity == 0.0)
300  return 0.0;
301 
302  const Real temp_func =
303  (_temp[_qp] * _temp[_qp] + _gap_temp * _gap_temp) * (_temp[_qp] + _gap_temp);
304  return _stefan_boltzmann * temp_func / _emissivity;
305 }
306 
307 Real
309 {
310  if (_emissivity == 0.0)
311  return 0.0;
312 
313  const Real temp_func = 3 * _temp[_qp] * _temp[_qp] + _gap_temp * (2 * _temp[_qp] + _gap_temp);
314  return _stefan_boltzmann * temp_func / _emissivity;
315 }
316 
317 Real
319  Real radius,
320  Real r1,
321  Real r2,
322  Real min_gap,
323  Real max_gap)
324 {
325  if (gap_geom == GapConductance::CYLINDER)
326  return gapCyl(radius, r1, r2, min_gap, max_gap);
327  else if (gap_geom == GapConductance::SPHERE)
328  return gapSphere(radius, r1, r2, min_gap, max_gap);
329  else
330  return gapRect(r2 - r1, min_gap, max_gap);
331 }
332 
333 Real
334 GapConductance::gapRect(Real distance, Real min_gap, Real max_gap)
335 {
336  return std::max(min_gap, std::min(distance, max_gap));
337 }
338 
339 Real
340 GapConductance::gapCyl(Real radius, Real r1, Real r2, Real min_denom, Real max_denom)
341 {
342  Real denominator = radius * std::log(r2 / r1);
343  return std::max(min_denom, std::min(denominator, max_denom));
344 }
345 
346 Real
347 GapConductance::gapSphere(Real radius, Real r1, Real r2, Real min_denom, Real max_denom)
348 {
349  Real denominator = std::pow(radius, 2.0) * ((1.0 / r1) - (1.0 / r2));
350  return std::max(min_denom, std::min(denominator, max_denom));
351 }
352 
353 Real
355 {
356  Real gap_conductivity = _gap_conductivity;
357 
359  {
361  gap_conductivity *= _gap_conductivity_function->value(
362  (*_gap_conductivity_function_variable)[_qp], _q_point[_qp]);
363  else
364  gap_conductivity *= _gap_conductivity_function->value(_t, _q_point[_qp]);
365  }
366 
367  return gap_conductivity;
368 }
369 
370 void
372 {
373  if (!_quadrature)
374  {
375  _has_info = true;
376  _gap_temp = _gap_temp_value[_qp];
378  }
379  else
380  {
381  Node * qnode = _mesh.getQuadratureNode(_current_elem, _current_side, _qp);
382  PenetrationInfo * pinfo = _penetration_locator->_penetration_info[qnode->id()];
383 
384  _gap_temp = 0.0;
385  _gap_distance = 88888;
386  _has_info = false;
387 
388  if (pinfo)
389  {
390  _gap_distance = pinfo->_distance;
391  _has_info = true;
392 
393  const Elem * slave_side = pinfo->_side;
394  std::vector<std::vector<Real>> & slave_side_phi = pinfo->_side_phi;
395  std::vector<dof_id_type> slave_side_dof_indices;
396 
397  _dof_map->dof_indices(slave_side, slave_side_dof_indices, _temp_var->number());
398 
399  for (unsigned int i = 0; i < slave_side_dof_indices.size(); ++i)
400  {
401  // The zero index is because we only have one point that the phis are evaluated at
402  _gap_temp += slave_side_phi[i][0] * (*(*_serialized_solution))(slave_side_dof_indices[i]);
403  }
404  }
405  else
406  {
407  if (_warnings)
408  mooseWarning("No gap value information found for node ",
409  qnode->id(),
410  " on processor ",
411  processor_id(),
412  " at coordinate ",
413  Point(*qnode));
414  }
415  }
416 
417  Point current_point(_q_point[_qp]);
419  _gap_geometry_type, current_point, _p1, _p2, _gap_distance, _normals[_qp], _r1, _r2, _radius);
420 }
421 
422 void
424  const Point & current_point,
425  const Point & p1,
426  const Point & p2,
427  const Real & gap_distance,
428  const Point & current_normal,
429  Real & r1,
430  Real & r2,
431  Real & radius)
432 {
433  if (gap_geometry_type == GapConductance::CYLINDER)
434  {
435  // The vector _p1 + t*(_p2-_p1) defines the cylindrical axis. The point along this
436  // axis closest to current_point is found by the following for t:
437  const Point p2p1(p2 - p1);
438  const Point p1pc(p1 - current_point);
439  const Real t = -(p1pc * p2p1) / p2p1.norm_sq();
440 
441  // The nearest point on the cylindrical axis to current_point is p.
442  const Point p(p1 + t * p2p1);
443  Point rad_vec(current_point - p);
444  Real rad = rad_vec.norm();
445  rad_vec /= rad;
446  Real rad_dot_norm = rad_vec * current_normal;
447 
448  if (rad_dot_norm > 0)
449  {
450  r1 = rad;
451  r2 = rad - gap_distance; // note, gap_distance is negative
452  radius = r1;
453  }
454  else if (rad_dot_norm < 0)
455  {
456  r1 = rad + gap_distance;
457  r2 = rad;
458  radius = r2;
459  }
460  else
461  ::mooseError("Issue with cylindrical flux calc. normals.\n");
462  }
463  else if (gap_geometry_type == GapConductance::SPHERE)
464  {
465  const Point origin_to_curr_point(current_point - p1);
466  const Real normal_dot = origin_to_curr_point * current_normal;
467  const Real curr_point_radius = origin_to_curr_point.norm();
468  if (normal_dot > 0) // on inside surface
469  {
470  r1 = curr_point_radius;
471  r2 = curr_point_radius - gap_distance; // gap_distance is negative
472  radius = r1;
473  }
474  else if (normal_dot < 0) // on outside surface
475  {
476  r1 = curr_point_radius + gap_distance; // gap_distance is negative
477  r2 = curr_point_radius;
478  radius = r2;
479  }
480  else
481  ::mooseError("Issue with spherical flux calc. normals. \n");
482  }
483  else
484  {
485  r2 = -gap_distance;
486  r1 = 0;
487  radius = 0;
488  }
489 }
static Real gapRect(Real distance, Real min_gap, Real max_gap)
MaterialProperty< Real > & _gap_conductance_dT
const VariableValue * _gap_conductivity_function_variable
const Real _stefan_boltzmann
static InputParameters actionParameters()
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)
MaterialProperty< Real > & _gap_thermal_conductivity
Function *const _gap_conductivity_function
const Real _gap_conductivity
const std::string _appended_property_name
InputParameters validParams< GapConductance >()
virtual Real h_radiation()
virtual Real dh_radiation()
virtual void initialSetup() override
virtual Real h_conduction()
const VariableValue & _gap_temp_value
static Real gapLength(const GAP_GEOMETRY &gap_geom, Real radius, Real r1, Real r2, Real min_gap, Real max_gap)
virtual Real dh_conduction()
static Real gapCyl(Real radius, Real r1, Real r2, Real min_denom, Real max_denom)
const NumericVector< Number > ** _serialized_solution
Generic gap heat transfer model, with h_gap = h_conduction + h_contact + h_radiation.
virtual void computeGapValues()
const bool _warnings
const VariableValue & _gap_distance_value
virtual void computeQpProperties() override
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
virtual void computeQpConductance()
Override this to compute the conductance at _qp.
GAP_GEOMETRY & _gap_geometry_type
virtual Real gapK()
const VariableValue & _temp
static Real gapSphere(Real radius, Real r1, Real r2, Real min_denom, Real max_denom)
PenetrationLocator * _penetration_locator
MooseVariable * _temp_var
GapConductance(const InputParameters &parameters)
static void setGapGeometryParameters(const InputParameters &params, const Moose::CoordinateSystemType coord_sys, GAP_GEOMETRY &gap_geometry_type, Point &p1, Point &p2)
MaterialProperty< Real > & _gap_conductance