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 "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 : #include "GhostBoundary.h"
20 :
21 : #include "libmesh/string_to_enum.h"
22 :
23 : registerMooseObject("HeatTransferApp", GapHeatTransfer);
24 :
25 : InputParameters
26 7196 : GapHeatTransfer::validParams()
27 : {
28 7196 : InputParameters params = IntegratedBC::validParams();
29 7196 : params.addClassDescription("Transfers heat across a gap between two "
30 : "surfaces dependent on the gap geometry specified.");
31 14392 : params.addParam<std::string>(
32 : "appended_property_name", "", "Name appended to material properties to make them unique");
33 :
34 : // Common
35 14392 : params.addParam<Real>("min_gap", 1.0e-6, "A minimum gap size");
36 14392 : params.addParam<Real>("max_gap", 1.0e6, "A maximum gap size");
37 21588 : params.addRangeCheckedParam<unsigned int>(
38 14392 : "min_gap_order", 0, "min_gap_order<=1", "Order of the Taylor expansion below min_gap");
39 :
40 : // Deprecated parameter
41 14392 : MooseEnum coord_types("default XYZ cyl", "default");
42 14392 : params.addDeprecatedParam<MooseEnum>(
43 : "coord_type",
44 : coord_types,
45 : "Gap calculation type (default or XYZ).",
46 : "The functionality of this parameter is replaced by 'gap_geometry_type'.");
47 :
48 14392 : MooseEnum gap_geom_types("PLATE CYLINDER SPHERE");
49 14392 : params.addParam<MooseEnum>("gap_geometry_type",
50 : gap_geom_types,
51 7196 : "Gap calculation type. Choices are: " + gap_geom_types.getRawNames());
52 :
53 14392 : params.addParam<RealVectorValue>("cylinder_axis_point_1",
54 : "Start point for line defining cylindrical axis");
55 14392 : params.addParam<RealVectorValue>("cylinder_axis_point_2",
56 : "End point for line defining cylindrical axis");
57 14392 : params.addParam<RealVectorValue>("sphere_origin", "Origin for sphere geometry");
58 :
59 : // Quadrature based
60 14392 : params.addParam<bool>("quadrature",
61 14392 : false,
62 : "Whether or not to do Quadrature point based gap heat "
63 : "transfer. If this is true then gap_distance and "
64 : "gap_temp should NOT be provided (and will be "
65 : "ignored) however paired_boundary IS then required.");
66 14392 : params.addParam<BoundaryName>("paired_boundary", "The boundary to be penetrated");
67 :
68 7196 : MooseEnum orders(AddVariableAction::getNonlinearVariableOrders());
69 14392 : params.addParam<MooseEnum>("order", orders, "The finite element order");
70 :
71 14392 : params.addParam<bool>(
72 14392 : "warnings", false, "Whether to output warning messages concerning nodes not being found");
73 :
74 14392 : params.addCoupledVar(
75 : "displacements",
76 : "The displacements appropriate for the simulation geometry and coordinate system");
77 :
78 : // Node based options
79 14392 : params.addCoupledVar("gap_distance", "Distance across the gap");
80 14392 : params.addCoupledVar("gap_temp", "Temperature on the other side of the gap");
81 :
82 14392 : params.addRelationshipManager(
83 : "GhostBoundary",
84 : Moose::RelationshipManagerType::GEOMETRIC,
85 5730 : [](const InputParameters & obj_params, InputParameters & rm_params)
86 : {
87 11460 : auto & boundary = rm_params.set<std::vector<BoundaryName>>("boundary");
88 5730 : boundary = obj_params.get<std::vector<BoundaryName>>("boundary");
89 5730 : boundary.push_back(obj_params.get<BoundaryName>("paired_boundary"));
90 5730 : });
91 :
92 7196 : return params;
93 7196 : }
94 :
95 941 : GapHeatTransfer::GapHeatTransfer(const InputParameters & parameters)
96 : : IntegratedBC(parameters),
97 941 : _gap_geometry_type(declareRestartableData<GapConductance::GAP_GEOMETRY>("gap_geometry_type",
98 1882 : GapConductance::PLATE)),
99 1882 : _quadrature(getParam<bool>("quadrature")),
100 941 : _secondary_flux(!_quadrature ? &_sys.getVector("secondary_flux") : NULL),
101 2823 : _gap_conductance(getMaterialProperty<Real>("gap_conductance" +
102 : getParam<std::string>("appended_property_name"))),
103 941 : _gap_conductance_dT(getMaterialProperty<Real>(
104 1882 : "gap_conductance" + getParam<std::string>("appended_property_name") + "_dT")),
105 1882 : _min_gap(getParam<Real>("min_gap")),
106 1882 : _min_gap_order(getParam<unsigned int>("min_gap_order")),
107 1882 : _max_gap(getParam<Real>("max_gap")),
108 941 : _gap_temp(0),
109 941 : _gap_distance(std::numeric_limits<Real>::max()),
110 941 : _edge_multiplier(1.0),
111 941 : _has_info(false),
112 941 : _disp_vars(3, libMesh::invalid_uint),
113 941 : _gap_distance_value(_quadrature ? _zero : coupledValue("gap_distance")),
114 941 : _gap_temp_value(_quadrature ? _zero : coupledValue("gap_temp")),
115 941 : _penetration_locator(
116 941 : !_quadrature ? NULL
117 2124 : : &getQuadraturePenetrationLocator(
118 : parameters.get<BoundaryName>("paired_boundary"),
119 2357 : getParam<std::vector<BoundaryName>>("boundary")[0],
120 708 : Utility::string_to_enum<Order>(parameters.get<MooseEnum>("order")))),
121 1882 : _warnings(getParam<bool>("warnings")),
122 1882 : _p1(declareRestartableData<Point>("cylinder_axis_point_1", Point(0, 1, 0))),
123 1882 : _p2(declareRestartableData<Point>("cylinder_axis_point_2", Point(0, 0, 0))),
124 941 : _pinfo(nullptr),
125 941 : _secondary_side_phi(nullptr),
126 941 : _secondary_side(nullptr),
127 941 : _secondary_j(0)
128 : {
129 1882 : if (isParamValid("displacements"))
130 : {
131 : // modern parameter scheme for displacements
132 240 : for (unsigned int i = 0; i < coupledComponents("displacements"); ++i)
133 80 : _disp_vars[i] = coupled("displacements", i);
134 : }
135 :
136 941 : if (_quadrature)
137 : {
138 1416 : if (!parameters.isParamValid("paired_boundary"))
139 0 : mooseError(std::string("No 'paired_boundary' provided for ") + _name);
140 : }
141 : else
142 : {
143 466 : if (!isCoupled("gap_distance"))
144 0 : mooseError(std::string("No 'gap_distance' provided for ") + _name);
145 :
146 466 : if (!isCoupled("gap_temp"))
147 0 : mooseError(std::string("No 'gap_temp' provided for ") + _name);
148 : }
149 941 : }
150 :
151 : void
152 941 : GapHeatTransfer::initialSetup()
153 : {
154 : std::set<SubdomainID> subdomain_ids;
155 :
156 93199 : for (const auto & bnd_elem : *_mesh.getBoundaryElementRange())
157 : {
158 92258 : Elem * elem = bnd_elem->_elem;
159 92258 : subdomain_ids.insert(elem->subdomain_id());
160 : }
161 :
162 941 : if (subdomain_ids.empty())
163 0 : mooseError("No boundary elements found");
164 :
165 941 : Moose::CoordinateSystemType coord_system = _fe_problem.getCoordSystem(*subdomain_ids.begin());
166 :
167 2793 : for (auto sid : subdomain_ids)
168 1852 : if (_fe_problem.getCoordSystem(sid) != coord_system)
169 0 : mooseError("The GapHeatTransfer model requires all boundary elements to have the same "
170 : "coordinate system.");
171 :
172 941 : GapConductance::setGapGeometryParameters(
173 941 : _pars, coord_system, _fe_problem.getAxisymmetricRadialCoord(), _gap_geometry_type, _p1, _p2);
174 941 : }
175 :
176 : Real
177 15109682 : GapHeatTransfer::computeQpResidual()
178 : {
179 15109682 : computeGapValues();
180 :
181 15109682 : if (!_has_info)
182 : return 0.0;
183 :
184 13610190 : Real grad_t = (_u[_qp] - _gap_temp) * _edge_multiplier * _gap_conductance[_qp];
185 :
186 : // This is keeping track of this residual contribution so it can be used as the flux on the other
187 : // side of the gap.
188 13610190 : if (!_quadrature)
189 : {
190 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mutex);
191 6518933 : const Real secondary_flux = computeSecondaryFluxContribution(grad_t);
192 6518933 : _secondary_flux->add(_var.dofIndices()[_i], secondary_flux);
193 : }
194 :
195 13610190 : return _test[_i][_qp] * grad_t;
196 : }
197 :
198 : Real
199 6518933 : GapHeatTransfer::computeSecondaryFluxContribution(Real grad_t)
200 : {
201 6518933 : return _coord[_qp] * _JxW[_qp] * _test[_i][_qp] * grad_t;
202 : }
203 :
204 : void
205 38191 : GapHeatTransfer::computeJacobian()
206 : {
207 38191 : prepareMatrixTag(_assembly, _var.number(), _var.number());
208 :
209 219035 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
210 : {
211 : // compute this up front because it only depends on the quadrature point
212 180844 : computeGapValues();
213 :
214 2324245 : for (_i = 0; _i < _test.size(); _i++)
215 38800362 : for (_j = 0; _j < _phi.size(); _j++)
216 36656961 : _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJacobian();
217 :
218 : // Ok now do the contribution from the secondary side
219 180844 : if (_quadrature && _has_info)
220 : {
221 : std::vector<dof_id_type> secondary_side_dof_indices;
222 :
223 112469 : _sys.dofMap().dof_indices(_secondary_side, secondary_side_dof_indices, _var.number());
224 :
225 112469 : DenseMatrix<Number> K_secondary(_var.dofIndices().size(), secondary_side_dof_indices.size());
226 :
227 : mooseAssert(
228 : _secondary_side_phi->size() == secondary_side_dof_indices.size(),
229 : "The number of shapes does not match the number of dof indices on the secondary elem");
230 :
231 1218569 : for (_i = 0; _i < _test.size(); _i++)
232 7488124 : for (_secondary_j = 0;
233 7488124 : _secondary_j < static_cast<unsigned int>(secondary_side_dof_indices.size());
234 : ++_secondary_j)
235 6382024 : K_secondary(_i, _secondary_j) += _JxW[_qp] * _coord[_qp] * computeSecondaryQpJacobian();
236 :
237 112469 : addJacobian(_subproblem.assembly(_tid, _sys.number()),
238 : K_secondary,
239 : _var.dofIndices(),
240 : secondary_side_dof_indices,
241 112469 : _var.scalingFactor());
242 112469 : }
243 : }
244 :
245 38191 : accumulateTaggedLocalMatrix();
246 :
247 38191 : if (_has_diag_save_in)
248 : {
249 : unsigned int rows = _local_ke.m();
250 0 : DenseVector<Number> diag(rows);
251 0 : for (unsigned int i = 0; i < rows; i++)
252 0 : diag(i) = _local_ke(i, i);
253 :
254 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
255 0 : for (unsigned int i = 0; i < _diag_save_in.size(); i++)
256 0 : _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
257 : }
258 38191 : }
259 :
260 : void
261 6904 : GapHeatTransfer::computeOffDiagJacobian(const unsigned int jvar_num)
262 : {
263 6904 : if (jvar_num == _var.number())
264 : {
265 4600 : computeJacobian();
266 4600 : return;
267 : }
268 :
269 : const auto & jvar = getVariable(jvar_num);
270 :
271 2304 : prepareMatrixTag(_assembly, _var.number(), jvar_num);
272 :
273 2304 : auto phi_size = jvar.dofIndices().size();
274 :
275 6912 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
276 : {
277 : // compute this up front because it only depends on the quadrature point
278 4608 : computeGapValues();
279 :
280 23040 : for (_i = 0; _i < _test.size(); _i++)
281 92160 : for (_j = 0; _j < phi_size; _j++)
282 73728 : _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobian(jvar_num);
283 :
284 : // Ok now do the contribution from the secondary side
285 4608 : if (_quadrature && _has_info)
286 : {
287 : std::vector<dof_id_type> secondary_side_dof_indices;
288 :
289 4608 : _sys.dofMap().dof_indices(_secondary_side, secondary_side_dof_indices, jvar_num);
290 :
291 4608 : DenseMatrix<Number> K_secondary(_var.dofIndices().size(), secondary_side_dof_indices.size());
292 :
293 : mooseAssert(
294 : _secondary_side_phi->size() == secondary_side_dof_indices.size(),
295 : "The number of shapes does not match the number of dof indices on the secondary elem");
296 :
297 23040 : for (_i = 0; _i < _test.size(); _i++)
298 55296 : for (_secondary_j = 0;
299 55296 : _secondary_j < static_cast<unsigned int>(secondary_side_dof_indices.size());
300 : ++_secondary_j)
301 36864 : K_secondary(_i, _secondary_j) +=
302 36864 : _JxW[_qp] * _coord[_qp] * computeSecondaryQpOffDiagJacobian(jvar_num);
303 :
304 4608 : addJacobian(_subproblem.assembly(_tid, _sys.number()),
305 : K_secondary,
306 : _var.dofIndices(),
307 : secondary_side_dof_indices,
308 4608 : _var.scalingFactor());
309 4608 : }
310 : }
311 :
312 2304 : accumulateTaggedLocalMatrix();
313 : }
314 :
315 : Real
316 36656961 : GapHeatTransfer::computeQpJacobian()
317 : {
318 36656961 : if (!_has_info)
319 : return 0.0;
320 :
321 33640281 : return _test[_i][_qp] *
322 33640281 : ((_u[_qp] - _gap_temp) * _edge_multiplier * _gap_conductance_dT[_qp] +
323 33640281 : _edge_multiplier * _gap_conductance[_qp]) *
324 33640281 : _phi[_j][_qp];
325 : }
326 :
327 : Real
328 6382024 : GapHeatTransfer::computeSecondaryQpJacobian()
329 : {
330 6382024 : return _test[_i][_qp] *
331 6382024 : ((_u[_qp] - _gap_temp) * _edge_multiplier * _gap_conductance_dT[_qp] -
332 6382024 : _edge_multiplier * _gap_conductance[_qp]) *
333 6382024 : (*_secondary_side_phi)[_secondary_j][0];
334 : }
335 :
336 : Real
337 73728 : GapHeatTransfer::computeQpOffDiagJacobian(unsigned jvar)
338 : {
339 73728 : if (!_has_info)
340 : return 0.0;
341 :
342 : unsigned int coupled_component;
343 : bool active = false;
344 110592 : for (coupled_component = 0; coupled_component < _disp_vars.size(); ++coupled_component)
345 110592 : if (jvar == _disp_vars[coupled_component])
346 : {
347 : active = true;
348 : break;
349 : }
350 :
351 : Real dRdx = 0.0;
352 73728 : if (active)
353 : {
354 : // Compute dR/du_[xyz]
355 : // Residual is based on
356 : // h_gap = h_conduction() + h_contact() + h_radiation();
357 : // grad_t = (_u[_qp] - _gap_temp) * h_gap;
358 : // So we need
359 : // (_u[_qp] - _gap_temp) * (dh_gap/du_[xyz]);
360 : // Assuming dh_contact/du_[xyz] = dh_radiation/du_[xyz] = 0,
361 : // we need dh_conduction/du_[xyz]
362 : // Given
363 : // h_conduction = gapK / gapLength, then
364 : // dh_conduction/du_[xyz] = -gapK/gapLength^2 * dgapLength/du_[xyz]
365 : // Given
366 : // gapLength = ((u_x-m_x)^2+(u_y-m_y)^2+(u_z-m_z)^2)^1/2
367 : // where m_[xyz] is the primary coordinate, then
368 : // dGapLength/du_[xyz] =
369 : // 1/2*((u_x-m_x)^2+(u_y-m_y)^2+(u_z-m_z)^2)^(-1/2)*2*(u_[xyz]-m_[xyz])
370 : // = (u_[xyz]-m_[xyz])/gapLength
371 : // This is the normal vector.
372 :
373 73728 : const Real gapL = gapLength();
374 :
375 : // THIS IS NOT THE NORMAL WE NEED.
376 : // WE NEED THE NORMAL FROM THE CONSTRAINT, THE NORMAL FROM THE
377 : // PRIMARY SURFACE. HOWEVER, THIS IS TRICKY SINCE THE NORMAL
378 : // FROM THE PRIMARY SURFACE WAS COMPUTED FOR A POINT ASSOCIATED
379 : // WITH A SECONDARY NODE. NOW WE ARE AT A SECONDARY INTEGRATION POINT.
380 : //
381 : // HOW DO WE GET THE NORMAL WE NEED?
382 : //
383 : // Until we have the normal we need,
384 : // we'll hope that the one we have is close to the negative of the one we need.
385 73728 : const Point & normal(_normals[_qp]);
386 :
387 73728 : const Real dgap = dgapLength(-normal(coupled_component));
388 73728 : dRdx = -(_u[_qp] - _gap_temp) * _edge_multiplier * _gap_conductance[_qp] *
389 73728 : GapConductance::gapAttenuation(gapL, _min_gap, _min_gap_order) * dgap;
390 : }
391 73728 : return _test[_i][_qp] * dRdx * _phi[_j][_qp];
392 : }
393 :
394 : Real
395 36864 : GapHeatTransfer::computeSecondaryQpOffDiagJacobian(unsigned jvar)
396 : {
397 36864 : if (!_has_info)
398 : return 0.0;
399 :
400 : unsigned int coupled_component;
401 : bool active = false;
402 55296 : for (coupled_component = 0; coupled_component < _disp_vars.size(); ++coupled_component)
403 55296 : if (jvar == _disp_vars[coupled_component])
404 : {
405 : active = true;
406 : break;
407 : }
408 :
409 : Real dRdx = 0.0;
410 36864 : if (active)
411 : {
412 36864 : const Real gapL = gapLength();
413 :
414 36864 : const Point & normal(_normals[_qp]);
415 :
416 36864 : const Real dgap = dgapLength(-normal(coupled_component));
417 :
418 : // The sign of the secondary side should presumably be opposite that of the primary side
419 36864 : dRdx = (_u[_qp] - _gap_temp) * _edge_multiplier * _gap_conductance[_qp] *
420 36864 : GapConductance::gapAttenuation(gapL, _min_gap, _min_gap_order) * dgap;
421 : }
422 36864 : return _test[_i][_qp] * dRdx * (*_secondary_side_phi)[_secondary_j][0];
423 : }
424 :
425 : Real
426 221184 : GapHeatTransfer::gapLength() const
427 : {
428 221184 : if (_has_info)
429 221184 : return GapConductance::gapLength(_gap_geometry_type, _radius, _r1, _r2, _max_gap);
430 :
431 : return 1.0;
432 : }
433 :
434 : Real
435 110592 : GapHeatTransfer::dgapLength(Real normalComponent) const
436 : {
437 110592 : const Real gap_L = gapLength();
438 : Real dgap = 0.0;
439 :
440 110592 : if (_min_gap <= gap_L && gap_L <= _max_gap)
441 : dgap = normalComponent;
442 :
443 110592 : return dgap;
444 : }
445 :
446 : void
447 15295134 : GapHeatTransfer::computeGapValues()
448 : {
449 15295134 : if (!_quadrature)
450 : {
451 6567517 : _has_info = true;
452 6567517 : _gap_temp = _gap_temp_value[_qp];
453 6567517 : _gap_distance = _gap_distance_value[_qp];
454 : }
455 : else
456 : {
457 8727617 : Node * qnode = _mesh.getQuadratureNode(_current_elem, _current_side, _qp);
458 :
459 8727617 : _gap_temp = 0.0;
460 8727617 : _gap_distance = std::numeric_limits<Real>::max();
461 8727617 : _has_info = false;
462 8727617 : _edge_multiplier = 1.0;
463 :
464 8727617 : if (auto it = _penetration_locator->_penetration_info.find(qnode->id());
465 8727617 : it != _penetration_locator->_penetration_info.end() && it->second)
466 : {
467 7208334 : _pinfo = it->second;
468 7208334 : _gap_distance = _pinfo->_distance;
469 7208334 : _has_info = true;
470 :
471 7208334 : _secondary_side = _pinfo->_side;
472 7208334 : _secondary_side_phi = &_pinfo->_side_phi;
473 7208334 : _gap_temp = _variable->getValue(_secondary_side, *_secondary_side_phi);
474 :
475 7208334 : Real tangential_tolerance = _penetration_locator->getTangentialTolerance();
476 7208334 : if (tangential_tolerance != 0.0)
477 : {
478 0 : _edge_multiplier = 1.0 - _pinfo->_tangential_distance / tangential_tolerance;
479 0 : if (_edge_multiplier < 0.0)
480 0 : _edge_multiplier = 0.0;
481 : }
482 : }
483 : else
484 : {
485 1519283 : _pinfo = nullptr;
486 1519283 : if (_warnings)
487 0 : mooseWarning("No gap value information found for node ",
488 0 : qnode->id(),
489 : " on processor ",
490 0 : processor_id());
491 : }
492 : }
493 :
494 15295134 : GapConductance::computeGapRadii(
495 15295134 : _gap_geometry_type, _q_point[_qp], _p1, _p2, _gap_distance, _normals[_qp], _r1, _r2, _radius);
496 15295134 : }
|
