www.mooseframework.org
GapHeatTransfer.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 #include "GapHeatTransfer.h"
8 
9 // MOOSE includes
10 #include "AddVariableAction.h"
11 #include "Assembly.h"
12 #include "MooseMesh.h"
13 #include "MooseVariable.h"
14 #include "PenetrationLocator.h"
15 #include "SystemBase.h"
16 
17 #include "libmesh/string_to_enum.h"
18 
19 template <>
20 InputParameters
22 {
23  InputParameters params = validParams<IntegratedBC>();
24  params.addClassDescription("Transfers heat across a gap between two "
25  "surfaces dependant on the gap geometry specified.");
26  params.addParam<std::string>(
27  "appended_property_name", "", "Name appended to material properties to make them unique");
28 
29  // Common
30  params.addParam<Real>("min_gap", 1.0e-6, "A minimum gap size");
31  params.addParam<Real>("max_gap", 1.0e6, "A maximum gap size");
32 
33  // Deprecated parameter
34  MooseEnum coord_types("default XYZ cyl", "default");
35  params.addDeprecatedParam<MooseEnum>(
36  "coord_type",
37  coord_types,
38  "Gap calculation type (default or XYZ).",
39  "The functionality of this parameter is replaced by 'gap_geometry_type'.");
40 
41  MooseEnum gap_geom_types("PLATE CYLINDER SPHERE");
42  params.addParam<MooseEnum>("gap_geometry_type",
43  gap_geom_types,
44  "Gap calculation type. Choices are: " + gap_geom_types.getRawNames());
45 
46  params.addParam<RealVectorValue>("cylinder_axis_point_1",
47  "Start point for line defining cylindrical axis");
48  params.addParam<RealVectorValue>("cylinder_axis_point_2",
49  "End point for line defining cylindrical axis");
50  params.addParam<RealVectorValue>("sphere_origin", "Origin for sphere geometry");
51 
52  // Quadrature based
53  params.addParam<bool>("quadrature",
54  false,
55  "Whether or not to do Quadrature point based gap heat "
56  "transfer. If this is true then gap_distance and "
57  "gap_temp should NOT be provided (and will be "
58  "ignored) however paired_boundary IS then required.");
59  params.addParam<BoundaryName>("paired_boundary", "The boundary to be penetrated");
60 
61  MooseEnum orders(AddVariableAction::getNonlinearVariableOrders());
62  params.addParam<MooseEnum>("order", orders, "The finite element order");
63 
64  params.addParam<bool>(
65  "warnings", false, "Whether to output warning messages concerning nodes not being found");
66 
67  params.addCoupledVar("disp_x", "The x displacement");
68  params.addCoupledVar("disp_y", "The y displacement");
69  params.addCoupledVar("disp_z", "The z displacement");
70 
71  params.addCoupledVar(
72  "displacements",
73  "The displacements appropriate for the simulation geometry and coordinate system");
74 
75  // Node based options
76  params.addCoupledVar("gap_distance", "Distance across the gap");
77  params.addCoupledVar("gap_temp", "Temperature on the other side of the gap");
78 
79  return params;
80 }
81 
82 GapHeatTransfer::GapHeatTransfer(const InputParameters & parameters)
83  : IntegratedBC(parameters),
84  _gap_geometry_type(declareRestartableData<GapConductance::GAP_GEOMETRY>("gap_geometry_type",
85  GapConductance::PLATE)),
86  _quadrature(getParam<bool>("quadrature")),
87  _slave_flux(!_quadrature ? &_sys.getVector("slave_flux") : NULL),
88  _gap_conductance(getMaterialProperty<Real>("gap_conductance" +
89  getParam<std::string>("appended_property_name"))),
90  _gap_conductance_dT(getMaterialProperty<Real>(
91  "gap_conductance" + getParam<std::string>("appended_property_name") + "_dT")),
92  _min_gap(getParam<Real>("min_gap")),
93  _max_gap(getParam<Real>("max_gap")),
94  _gap_temp(0),
95  _gap_distance(std::numeric_limits<Real>::max()),
96  _edge_multiplier(1.0),
97  _has_info(false),
98  _disp_vars(3, libMesh::invalid_uint),
99  _gap_distance_value(_quadrature ? _zero : coupledValue("gap_distance")),
100  _gap_temp_value(_quadrature ? _zero : coupledValue("gap_temp")),
101  _penetration_locator(
102  !_quadrature ? NULL
103  : &getQuadraturePenetrationLocator(
104  parameters.get<BoundaryName>("paired_boundary"),
105  getParam<std::vector<BoundaryName>>("boundary")[0],
106  Utility::string_to_enum<Order>(parameters.get<MooseEnum>("order")))),
107  _warnings(getParam<bool>("warnings")),
108  _p1(declareRestartableData<Point>("cylinder_axis_point_1", Point(0, 1, 0))),
109  _p2(declareRestartableData<Point>("cylinder_axis_point_2", Point(0, 0, 0)))
110 {
111  if (isParamValid("displacements"))
112  {
113  // modern parameter scheme for displacements
114  for (unsigned int i = 0; i < coupledComponents("displacements"); ++i)
115  _disp_vars[i] = coupled("displacements", i);
116  }
117  else
118  {
119  // Legacy parameter scheme for displacements
120  if (isParamValid("disp_x"))
121  _disp_vars[0] = coupled("disp_x");
122  if (isParamValid("disp_y"))
123  _disp_vars[1] = coupled("disp_y");
124  if (isParamValid("disp_z"))
125  _disp_vars[2] = coupled("disp_z");
126 
127  // TODO: these are only used in one Bison test. Deprecate soon!
128  }
129 
130  if (_quadrature)
131  {
132  if (!parameters.isParamValid("paired_boundary"))
133  mooseError(std::string("No 'paired_boundary' provided for ") + _name);
134  }
135  else
136  {
137  if (!isCoupled("gap_distance"))
138  mooseError(std::string("No 'gap_distance' provided for ") + _name);
139 
140  if (!isCoupled("gap_temp"))
141  mooseError(std::string("No 'gap_temp' provided for ") + _name);
142  }
143 }
144 
145 void
147 {
149  _pars, _assembly.coordSystem(), _gap_geometry_type, _p1, _p2);
150 }
151 
152 Real
154 {
156 
157  if (!_has_info)
158  return 0.0;
159 
160  Real grad_t = (_u[_qp] - _gap_temp) * _edge_multiplier * _gap_conductance[_qp];
161 
162  // This is keeping track of this residual contribution so it can be used as the flux on the other
163  // side of the gap.
164  if (!_quadrature)
165  {
166  Threads::spin_mutex::scoped_lock lock(Threads::spin_mutex);
167  const Real slave_flux = computeSlaveFluxContribution(grad_t);
168  _slave_flux->add(_var.dofIndices()[_i], slave_flux);
169  }
170 
171  return _test[_i][_qp] * grad_t;
172 }
173 
174 Real
176 {
177  return _coord[_qp] * _JxW[_qp] * _test[_i][_qp] * grad_t;
178 }
179 
180 Real
182 {
184 
185  if (!_has_info)
186  return 0.0;
187 
188  return _test[_i][_qp] * ((_u[_qp] - _gap_temp) * _edge_multiplier * _gap_conductance_dT[_qp] +
190  _phi[_j][_qp];
191 }
192 
193 Real
195 {
197 
198  if (!_has_info)
199  return 0.0;
200 
201  unsigned int coupled_component;
202  bool active = false;
203  for (coupled_component = 0; coupled_component < _disp_vars.size(); ++coupled_component)
204  if (jvar == _disp_vars[coupled_component])
205  {
206  active = true;
207  break;
208  }
209 
210  Real dRdx = 0.0;
211  if (active)
212  {
213  // Compute dR/du_[xyz]
214  // Residual is based on
215  // h_gap = h_conduction() + h_contact() + h_radiation();
216  // grad_t = (_u[_qp] - _gap_temp) * h_gap;
217  // So we need
218  // (_u[_qp] - _gap_temp) * (dh_gap/du_[xyz]);
219  // Assuming dh_contact/du_[xyz] = dh_radiation/du_[xyz] = 0,
220  // we need dh_conduction/du_[xyz]
221  // Given
222  // h_conduction = gapK / gapLength, then
223  // dh_conduction/du_[xyz] = -gapK/gapLength^2 * dgapLength/du_[xyz]
224  // Given
225  // gapLength = ((u_x-m_x)^2+(u_y-m_y)^2+(u_z-m_z)^2)^1/2
226  // where m_[xyz] is the master coordinate, then
227  // dGapLength/du_[xyz] = 1/2*((u_x-m_x)^2+(u_y-m_y)^2+(u_z-m_z)^2)^(-1/2)*2*(u_[xyz]-m_[xyz])
228  // = (u_[xyz]-m_[xyz])/gapLength
229  // This is the normal vector.
230 
231  const Real gapL = gapLength();
232 
233  // THIS IS NOT THE NORMAL WE NEED.
234  // WE NEED THE NORMAL FROM THE CONSTRAINT, THE NORMAL FROM THE
235  // MASTER SURFACE. HOWEVER, THIS IS TRICKY SINCE THE NORMAL
236  // FROM THE MASTER SURFACE WAS COMPUTED FOR A POINT ASSOCIATED
237  // WITH A SLAVE NODE. NOW WE ARE AT A SLAVE INTEGRATION POINT.
238  //
239  // HOW DO WE GET THE NORMAL WE NEED?
240  //
241  // Until we have the normal we need,
242  // we'll hope that the one we have is close to the negative of the one we need.
243  const Point & normal(_normals[_qp]);
244 
245  const Real dgap = dgapLength(-normal(coupled_component));
246  dRdx = -(_u[_qp] - _gap_temp) * _edge_multiplier * _gap_conductance[_qp] / gapL * dgap;
247  }
248  return _test[_i][_qp] * dRdx * _phi[_j][_qp];
249 }
250 
251 Real
253 {
254  if (_has_info)
256 
257  return 1.0;
258 }
259 
260 Real
261 GapHeatTransfer::dgapLength(Real normalComponent) const
262 {
263  const Real gap_L = gapLength();
264  Real dgap = 0.0;
265 
266  if (_min_gap <= gap_L && gap_L <= _max_gap)
267  dgap = normalComponent;
268 
269  return dgap;
270 }
271 
272 void
274 {
275  if (!_quadrature)
276  {
277  _has_info = true;
278  _gap_temp = _gap_temp_value[_qp];
280  }
281  else
282  {
283  Node * qnode = _mesh.getQuadratureNode(_current_elem, _current_side, _qp);
284  PenetrationInfo * pinfo = _penetration_locator->_penetration_info[qnode->id()];
285 
286  _gap_temp = 0.0;
287  _gap_distance = std::numeric_limits<Real>::max();
288  _has_info = false;
289  _edge_multiplier = 1.0;
290 
291  if (pinfo)
292  {
293  _gap_distance = pinfo->_distance;
294  _has_info = true;
295 
296  const Elem * slave_side = pinfo->_side;
297  std::vector<std::vector<Real>> & slave_side_phi = pinfo->_side_phi;
298  _gap_temp = _variable->getValue(slave_side, slave_side_phi);
299 
300  Real tangential_tolerance = _penetration_locator->getTangentialTolerance();
301  if (tangential_tolerance != 0.0)
302  {
303  _edge_multiplier = 1.0 - pinfo->_tangential_distance / tangential_tolerance;
304  if (_edge_multiplier < 0.0)
305  _edge_multiplier = 0.0;
306  }
307  }
308  else
309  {
310  if (_warnings)
311  mooseWarning("No gap value information found for node ",
312  qnode->id(),
313  " on processor ",
314  processor_id());
315  }
316  }
317 
319  _gap_geometry_type, _q_point[_qp], _p1, _p2, _gap_distance, _normals[_qp], _r1, _r2, _radius);
320 }
GapConductance::GAP_GEOMETRY & _gap_geometry_type
virtual Real dgapLength(Real normalComponent) const
virtual void computeGapValues()
virtual Real gapLength() const
InputParameters validParams< GapHeatTransfer >()
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)
const MaterialProperty< Real > & _gap_conductance_dT
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
virtual Real computeQpJacobian() override
const MaterialProperty< Real > & _gap_conductance
static Real gapLength(const GAP_GEOMETRY &gap_geom, Real radius, Real r1, Real r2, Real min_gap, Real max_gap)
const Real _min_gap
const Real _max_gap
Generic gap heat transfer model, with h_gap = h_conduction + h_contact + h_radiation.
std::vector< unsigned int > _disp_vars
const bool _warnings
PenetrationLocator * _penetration_locator
const VariableValue & _gap_temp_value
virtual Real computeSlaveFluxContribution(Real grad_t)
GapHeatTransfer(const InputParameters &parameters)
virtual Real computeQpResidual() override
const VariableValue & _gap_distance_value
virtual Real computeQpOffDiagJacobian(unsigned jvar) override
virtual void initialSetup() override
static void setGapGeometryParameters(const InputParameters &params, const Moose::CoordinateSystemType coord_sys, GAP_GEOMETRY &gap_geometry_type, Point &p1, Point &p2)
NumericVector< Number > * _slave_flux