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 14712 : GapHeatTransfer::validParams()
27 : {
28 14712 : InputParameters params = IntegratedBC::validParams();
29 14712 : params.addClassDescription("Transfers heat across a gap between two "
30 : "surfaces dependent on the gap geometry specified.");
31 29424 : params.addParam<std::string>(
32 : "appended_property_name", "", "Name appended to material properties to make them unique");
33 :
34 : // Common
35 29424 : params.addParam<Real>("min_gap", 1.0e-6, "A minimum gap size");
36 29424 : params.addParam<Real>("max_gap", 1.0e6, "A maximum gap size");
37 44136 : params.addRangeCheckedParam<unsigned int>(
38 29424 : "min_gap_order", 0, "min_gap_order<=1", "Order of the Taylor expansion below min_gap");
39 :
40 : // Deprecated parameter
41 29424 : MooseEnum coord_types("default XYZ cyl", "default");
42 29424 : 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 29424 : MooseEnum gap_geom_types("PLATE CYLINDER SPHERE");
49 29424 : params.addParam<MooseEnum>("gap_geometry_type",
50 : gap_geom_types,
51 14712 : "Gap calculation type. Choices are: " + gap_geom_types.getRawNames());
52 :
53 29424 : params.addParam<RealVectorValue>("cylinder_axis_point_1",
54 : "Start point for line defining cylindrical axis");
55 29424 : params.addParam<RealVectorValue>("cylinder_axis_point_2",
56 : "End point for line defining cylindrical axis");
57 29424 : params.addParam<RealVectorValue>("sphere_origin", "Origin for sphere geometry");
58 :
59 : // Quadrature based
60 29424 : params.addParam<bool>("quadrature",
61 29424 : 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 29424 : params.addParam<BoundaryName>("paired_boundary", "The boundary to be penetrated");
67 :
68 14712 : MooseEnum orders(AddVariableAction::getNonlinearVariableOrders());
69 29424 : params.addParam<MooseEnum>("order", orders, "The finite element order");
70 :
71 29424 : params.addParam<bool>(
72 29424 : "warnings", false, "Whether to output warning messages concerning nodes not being found");
73 :
74 29424 : params.addCoupledVar(
75 : "displacements",
76 : "The displacements appropriate for the simulation geometry and coordinate system");
77 :
78 : // Node based options
79 29424 : params.addCoupledVar("gap_distance", "Distance across the gap");
80 29424 : params.addCoupledVar("gap_temp", "Temperature on the other side of the gap");
81 :
82 29424 : params.addRelationshipManager(
83 : "GhostBoundary",
84 : Moose::RelationshipManagerType::GEOMETRIC,
85 11646 : [](const InputParameters & obj_params, InputParameters & rm_params)
86 : {
87 23292 : auto & boundary = rm_params.set<std::vector<BoundaryName>>("boundary");
88 11646 : boundary = obj_params.get<std::vector<BoundaryName>>("boundary");
89 11646 : boundary.push_back(obj_params.get<BoundaryName>("paired_boundary"));
90 11646 : });
91 :
92 14712 : return params;
93 14712 : }
94 :
95 1991 : GapHeatTransfer::GapHeatTransfer(const InputParameters & parameters)
96 : : IntegratedBC(parameters),
97 1991 : _gap_geometry_type(declareRestartableData<GapConductance::GAP_GEOMETRY>("gap_geometry_type",
98 3982 : GapConductance::PLATE)),
99 3982 : _quadrature(getParam<bool>("quadrature")),
100 1991 : _secondary_flux(!_quadrature ? &_sys.getVector("secondary_flux") : NULL),
101 5973 : _gap_conductance(getMaterialProperty<Real>("gap_conductance" +
102 : getParam<std::string>("appended_property_name"))),
103 1991 : _gap_conductance_dT(getMaterialProperty<Real>(
104 3982 : "gap_conductance" + getParam<std::string>("appended_property_name") + "_dT")),
105 3982 : _min_gap(getParam<Real>("min_gap")),
106 3982 : _min_gap_order(getParam<unsigned int>("min_gap_order")),
107 3982 : _max_gap(getParam<Real>("max_gap")),
108 1991 : _gap_temp(0),
109 1991 : _gap_distance(std::numeric_limits<Real>::max()),
110 1991 : _edge_multiplier(1.0),
111 1991 : _has_info(false),
112 1991 : _disp_vars(3, libMesh::invalid_uint),
113 1991 : _gap_distance_value(_quadrature ? _zero : coupledValue("gap_distance")),
114 1991 : _gap_temp_value(_quadrature ? _zero : coupledValue("gap_temp")),
115 1991 : _penetration_locator(
116 1991 : !_quadrature ? NULL
117 4482 : : &getQuadraturePenetrationLocator(
118 : parameters.get<BoundaryName>("paired_boundary"),
119 4979 : getParam<std::vector<BoundaryName>>("boundary")[0],
120 1494 : Utility::string_to_enum<Order>(parameters.get<MooseEnum>("order")))),
121 3982 : _warnings(getParam<bool>("warnings")),
122 3982 : _p1(declareRestartableData<Point>("cylinder_axis_point_1", Point(0, 1, 0))),
123 3982 : _p2(declareRestartableData<Point>("cylinder_axis_point_2", Point(0, 0, 0))),
124 1991 : _pinfo(nullptr),
125 1991 : _secondary_side_phi(nullptr),
126 1991 : _secondary_side(nullptr),
127 1991 : _secondary_j(0)
128 : {
129 3982 : if (isParamValid("displacements"))
130 : {
131 : // modern parameter scheme for displacements
132 528 : for (unsigned int i = 0; i < coupledComponents("displacements"); ++i)
133 176 : _disp_vars[i] = coupled("displacements", i);
134 : }
135 :
136 1991 : if (_quadrature)
137 : {
138 2988 : if (!parameters.isParamValid("paired_boundary"))
139 0 : mooseError(std::string("No 'paired_boundary' provided for ") + _name);
140 : }
141 : else
142 : {
143 994 : if (!isCoupled("gap_distance"))
144 0 : mooseError(std::string("No 'gap_distance' provided for ") + _name);
145 :
146 994 : if (!isCoupled("gap_temp"))
147 0 : mooseError(std::string("No 'gap_temp' provided for ") + _name);
148 : }
149 1991 : }
150 :
151 : void
152 1991 : GapHeatTransfer::initialSetup()
153 : {
154 : std::set<SubdomainID> subdomain_ids;
155 :
156 182597 : for (const auto & bnd_elem : *_mesh.getBoundaryElementRange())
157 : {
158 180606 : Elem * elem = bnd_elem->_elem;
159 180606 : subdomain_ids.insert(elem->subdomain_id());
160 : }
161 :
162 1991 : if (subdomain_ids.empty())
163 0 : mooseError("No boundary elements found");
164 :
165 1991 : Moose::CoordinateSystemType coord_system = _fe_problem.getCoordSystem(*subdomain_ids.begin());
166 :
167 5907 : for (auto sid : subdomain_ids)
168 3916 : 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 1991 : GapConductance::setGapGeometryParameters(
173 1991 : _pars, coord_system, _fe_problem.getAxisymmetricRadialCoord(), _gap_geometry_type, _p1, _p2);
174 1991 : }
175 :
176 : Real
177 20070876 : GapHeatTransfer::computeQpResidual()
178 : {
179 20070876 : computeGapValues();
180 :
181 20070876 : if (!_has_info)
182 : return 0.0;
183 :
184 17811732 : 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 17811732 : if (!_quadrature)
189 : {
190 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mutex);
191 7734225 : const Real secondary_flux = computeSecondaryFluxContribution(grad_t);
192 7734225 : _secondary_flux->add(_var.dofIndices()[_i], secondary_flux);
193 : }
194 :
195 17811732 : return _test[_i][_qp] * grad_t;
196 : }
197 :
198 : Real
199 7734225 : GapHeatTransfer::computeSecondaryFluxContribution(Real grad_t)
200 : {
201 7734225 : return _coord[_qp] * _JxW[_qp] * _test[_i][_qp] * grad_t;
202 : }
203 :
204 : void
205 55400 : GapHeatTransfer::computeJacobian()
206 : {
207 55400 : prepareMatrixTag(_assembly, _var.number(), _var.number());
208 :
209 305372 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
210 : {
211 : // compute this up front because it only depends on the quadrature point
212 249972 : computeGapValues();
213 :
214 3033088 : for (_i = 0; _i < _test.size(); _i++)
215 48530184 : for (_j = 0; _j < _phi.size(); _j++)
216 45747068 : _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJacobian();
217 :
218 : // Ok now do the contribution from the secondary side
219 249972 : if (_quadrature && _has_info)
220 : {
221 : std::vector<dof_id_type> secondary_side_dof_indices;
222 :
223 159030 : _sys.dofMap().dof_indices(_secondary_side, secondary_side_dof_indices, _var.number());
224 :
225 159030 : 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 1641592 : for (_i = 0; _i < _test.size(); _i++)
232 9657292 : for (_secondary_j = 0;
233 9657292 : _secondary_j < static_cast<unsigned int>(secondary_side_dof_indices.size());
234 : ++_secondary_j)
235 8174730 : K_secondary(_i, _secondary_j) += _JxW[_qp] * _coord[_qp] * computeSecondaryQpJacobian();
236 :
237 159030 : addJacobian(_subproblem.assembly(_tid, _sys.number()),
238 : K_secondary,
239 : _var.dofIndices(),
240 : secondary_side_dof_indices,
241 159030 : _var.scalingFactor());
242 159030 : }
243 : }
244 :
245 55400 : accumulateTaggedLocalMatrix();
246 :
247 55400 : 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 55400 : }
259 :
260 : void
261 10384 : GapHeatTransfer::computeOffDiagJacobian(const unsigned int jvar_num)
262 : {
263 10384 : if (jvar_num == _var.number())
264 : {
265 6928 : computeJacobian();
266 6928 : return;
267 : }
268 :
269 : const auto & jvar = getVariable(jvar_num);
270 :
271 3456 : prepareMatrixTag(_assembly, _var.number(), jvar_num);
272 :
273 3456 : auto phi_size = jvar.dofIndices().size();
274 :
275 10368 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
276 : {
277 : // compute this up front because it only depends on the quadrature point
278 6912 : computeGapValues();
279 :
280 34560 : for (_i = 0; _i < _test.size(); _i++)
281 138240 : for (_j = 0; _j < phi_size; _j++)
282 110592 : _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobian(jvar_num);
283 :
284 : // Ok now do the contribution from the secondary side
285 6912 : if (_quadrature && _has_info)
286 : {
287 : std::vector<dof_id_type> secondary_side_dof_indices;
288 :
289 6912 : _sys.dofMap().dof_indices(_secondary_side, secondary_side_dof_indices, jvar_num);
290 :
291 6912 : 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 34560 : for (_i = 0; _i < _test.size(); _i++)
298 82944 : for (_secondary_j = 0;
299 82944 : _secondary_j < static_cast<unsigned int>(secondary_side_dof_indices.size());
300 : ++_secondary_j)
301 55296 : K_secondary(_i, _secondary_j) +=
302 55296 : _JxW[_qp] * _coord[_qp] * computeSecondaryQpOffDiagJacobian(jvar_num);
303 :
304 6912 : addJacobian(_subproblem.assembly(_tid, _sys.number()),
305 : K_secondary,
306 : _var.dofIndices(),
307 : secondary_side_dof_indices,
308 6912 : _var.scalingFactor());
309 6912 : }
310 : }
311 :
312 3456 : accumulateTaggedLocalMatrix();
313 : }
314 :
315 : Real
316 45747068 : GapHeatTransfer::computeQpJacobian()
317 : {
318 45747068 : if (!_has_info)
319 : return 0.0;
320 :
321 41227848 : return _test[_i][_qp] *
322 41227848 : ((_u[_qp] - _gap_temp) * _edge_multiplier * _gap_conductance_dT[_qp] +
323 41227848 : _edge_multiplier * _gap_conductance[_qp]) *
324 41227848 : _phi[_j][_qp];
325 : }
326 :
327 : Real
328 8174730 : GapHeatTransfer::computeSecondaryQpJacobian()
329 : {
330 8174730 : return _test[_i][_qp] *
331 8174730 : ((_u[_qp] - _gap_temp) * _edge_multiplier * _gap_conductance_dT[_qp] -
332 8174730 : _edge_multiplier * _gap_conductance[_qp]) *
333 8174730 : (*_secondary_side_phi)[_secondary_j][0];
334 : }
335 :
336 : Real
337 110592 : GapHeatTransfer::computeQpOffDiagJacobian(unsigned jvar)
338 : {
339 110592 : if (!_has_info)
340 : return 0.0;
341 :
342 : unsigned int coupled_component;
343 : bool active = false;
344 165888 : for (coupled_component = 0; coupled_component < _disp_vars.size(); ++coupled_component)
345 165888 : if (jvar == _disp_vars[coupled_component])
346 : {
347 : active = true;
348 : break;
349 : }
350 :
351 : Real dRdx = 0.0;
352 110592 : 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 110592 : 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 110592 : const Point & normal(_normals[_qp]);
386 :
387 110592 : const Real dgap = dgapLength(-normal(coupled_component));
388 110592 : dRdx = -(_u[_qp] - _gap_temp) * _edge_multiplier * _gap_conductance[_qp] *
389 110592 : GapConductance::gapAttenuation(gapL, _min_gap, _min_gap_order) * dgap;
390 : }
391 110592 : return _test[_i][_qp] * dRdx * _phi[_j][_qp];
392 : }
393 :
394 : Real
395 55296 : GapHeatTransfer::computeSecondaryQpOffDiagJacobian(unsigned jvar)
396 : {
397 55296 : if (!_has_info)
398 : return 0.0;
399 :
400 : unsigned int coupled_component;
401 : bool active = false;
402 82944 : for (coupled_component = 0; coupled_component < _disp_vars.size(); ++coupled_component)
403 82944 : if (jvar == _disp_vars[coupled_component])
404 : {
405 : active = true;
406 : break;
407 : }
408 :
409 : Real dRdx = 0.0;
410 55296 : if (active)
411 : {
412 55296 : const Real gapL = gapLength();
413 :
414 55296 : const Point & normal(_normals[_qp]);
415 :
416 55296 : 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 55296 : dRdx = (_u[_qp] - _gap_temp) * _edge_multiplier * _gap_conductance[_qp] *
420 55296 : GapConductance::gapAttenuation(gapL, _min_gap, _min_gap_order) * dgap;
421 : }
422 55296 : return _test[_i][_qp] * dRdx * (*_secondary_side_phi)[_secondary_j][0];
423 : }
424 :
425 : Real
426 331776 : GapHeatTransfer::gapLength() const
427 : {
428 331776 : if (_has_info)
429 331776 : return GapConductance::gapLength(_gap_geometry_type, _radius, _r1, _r2, _max_gap);
430 :
431 : return 1.0;
432 : }
433 :
434 : Real
435 165888 : GapHeatTransfer::dgapLength(Real normalComponent) const
436 : {
437 165888 : const Real gap_L = gapLength();
438 : Real dgap = 0.0;
439 :
440 165888 : if (_min_gap <= gap_L && gap_L <= _max_gap)
441 : dgap = normalComponent;
442 :
443 165888 : return dgap;
444 : }
445 :
446 : void
447 20327760 : GapHeatTransfer::computeGapValues()
448 : {
449 20327760 : if (!_quadrature)
450 : {
451 7795579 : _has_info = true;
452 7795579 : _gap_temp = _gap_temp_value[_qp];
453 7795579 : _gap_distance = _gap_distance_value[_qp];
454 : }
455 : else
456 : {
457 12532181 : Node * qnode = _mesh.getQuadratureNode(_current_elem, _current_side, _qp);
458 12532181 : _pinfo = _penetration_locator->_penetration_info[qnode->id()];
459 :
460 12532181 : _gap_temp = 0.0;
461 12532181 : _gap_distance = std::numeric_limits<Real>::max();
462 12532181 : _has_info = false;
463 12532181 : _edge_multiplier = 1.0;
464 :
465 12532181 : if (_pinfo)
466 : {
467 10243449 : _gap_distance = _pinfo->_distance;
468 10243449 : _has_info = true;
469 :
470 10243449 : _secondary_side = _pinfo->_side;
471 10243449 : _secondary_side_phi = &_pinfo->_side_phi;
472 10243449 : _gap_temp = _variable->getValue(_secondary_side, *_secondary_side_phi);
473 :
474 10243449 : Real tangential_tolerance = _penetration_locator->getTangentialTolerance();
475 10243449 : if (tangential_tolerance != 0.0)
476 : {
477 0 : _edge_multiplier = 1.0 - _pinfo->_tangential_distance / tangential_tolerance;
478 0 : if (_edge_multiplier < 0.0)
479 0 : _edge_multiplier = 0.0;
480 : }
481 : }
482 : else
483 : {
484 2288732 : if (_warnings)
485 0 : mooseWarning("No gap value information found for node ",
486 0 : qnode->id(),
487 : " on processor ",
488 0 : processor_id());
489 : }
490 : }
491 :
492 20327760 : GapConductance::computeGapRadii(
493 20327760 : _gap_geometry_type, _q_point[_qp], _p1, _p2, _gap_distance, _normals[_qp], _r1, _r2, _radius);
494 20327760 : }
|
