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 "GapConductance.h"
11 :
12 : // MOOSE includes
13 : #include "Function.h"
14 : #include "MooseMesh.h"
15 : #include "MooseVariable.h"
16 : #include "PenetrationLocator.h"
17 : #include "SystemBase.h"
18 : #include "AddVariableAction.h"
19 :
20 : #include "libmesh/string_to_enum.h"
21 :
22 : registerMooseObject("HeatTransferApp", GapConductance);
23 :
24 : InputParameters
25 1382 : GapConductance::validParams()
26 : {
27 1382 : InputParameters params = Material::validParams();
28 1382 : params += GapConductance::actionParameters();
29 1382 : params.addClassDescription(
30 : "Material to compute an effective gap conductance based on the thermal conductivity of the "
31 : "gap and diffusive approximation to the radiative heat transfer.");
32 :
33 2764 : params.addRequiredCoupledVar("variable", "Temperature variable");
34 :
35 : // Node based
36 2764 : params.addCoupledVar("gap_distance", "Distance across the gap");
37 2764 : params.addCoupledVar("gap_temp", "Temperature on the other side of the gap");
38 2764 : params.addParam<Real>("gap_conductivity", 1.0, "The thermal conductivity of the gap material");
39 2764 : params.addParam<FunctionName>(
40 : "gap_conductivity_function",
41 : "Thermal conductivity of the gap material as a function. Multiplied by gap_conductivity.");
42 2764 : params.addCoupledVar("gap_conductivity_function_variable",
43 : "Variable to be used in the gap_conductivity_function in place of time");
44 :
45 : // Quadrature based
46 2764 : params.addParam<bool>("quadrature",
47 2764 : false,
48 : "Whether or not to do quadrature point based gap heat "
49 : "transfer. If this is true then gap_distance and "
50 : "gap_temp should NOT be provided (and will be "
51 : "ignored); however, paired_boundary and variable are "
52 : "then required.");
53 2764 : params.addParam<BoundaryName>("paired_boundary", "The boundary to be penetrated");
54 :
55 2764 : params.addParam<Real>("stefan_boltzmann", 5.670374e-8, "The Stefan-Boltzmann constant");
56 :
57 2764 : params.addParam<bool>("use_displaced_mesh",
58 2764 : true,
59 : "Whether or not this object should use the "
60 : "displaced mesh for computation. Note that in "
61 : "the case this is true but no displacements "
62 : "are provided in the Mesh block the "
63 : "undisplaced mesh will still be used.");
64 :
65 2764 : params.addParam<bool>(
66 2764 : "warnings", false, "Whether to output warning messages concerning nodes not being found");
67 :
68 1382 : MooseEnum orders(AddVariableAction::getNonlinearVariableOrders());
69 2764 : params.addParam<MooseEnum>("order", orders, "The finite element order");
70 :
71 : // These parameter groups match the ones in ThermalContactAction
72 2764 : params.addParamNamesToGroup("gap_distance", "Gap size");
73 2764 : params.addParamNamesToGroup("paired_boundary", "Gap surface definition");
74 2764 : params.addParamNamesToGroup("order quadrature", "Integration");
75 2764 : params.addParamNamesToGroup("warnings", "Diagnostics and debug");
76 2764 : params.addParamNamesToGroup(
77 : "gap_conductivity gap_conductivity_function gap_conductivity_function_variable",
78 : "Gap conductivity");
79 :
80 1382 : return params;
81 1382 : }
82 :
83 : InputParameters
84 4578 : GapConductance::actionParameters()
85 : {
86 4578 : InputParameters params = emptyInputParameters();
87 9156 : params.addParam<std::string>(
88 : "appended_property_name", "", "Name appended to material properties to make them unique");
89 9156 : MooseEnum gap_geom_types("PLATE CYLINDER SPHERE");
90 9156 : params.addParam<MooseEnum>("gap_geometry_type", gap_geom_types, "Gap calculation type.");
91 :
92 9156 : params.addParam<RealVectorValue>("cylinder_axis_point_1",
93 : "Start point for line defining cylindrical axis");
94 9156 : params.addParam<RealVectorValue>("cylinder_axis_point_2",
95 : "End point for line defining cylindrical axis");
96 9156 : params.addParam<RealVectorValue>("sphere_origin", "Origin for sphere geometry");
97 :
98 9156 : params.addRangeCheckedParam<Real>("emissivity_primary",
99 : 1,
100 : "emissivity_primary>=0 & emissivity_primary<=1",
101 : "The emissivity of the primary surface");
102 9156 : params.addRangeCheckedParam<Real>("emissivity_secondary",
103 : 1,
104 : "emissivity_secondary>=0 & emissivity_secondary<=1",
105 : "The emissivity of the secondary surface");
106 : // Common
107 13734 : params.addRangeCheckedParam<Real>(
108 9156 : "min_gap", 1e-6, "min_gap>0", "A minimum gap (denominator) size");
109 13734 : params.addRangeCheckedParam<Real>(
110 9156 : "max_gap", 1e6, "max_gap>=0", "A maximum gap (denominator) size");
111 13734 : params.addRangeCheckedParam<unsigned int>(
112 9156 : "min_gap_order", 0, "min_gap_order<=1", "Order of the Taylor expansion below min_gap");
113 :
114 9156 : params.addParamNamesToGroup("appended_property_name", "Material property retrieval");
115 9156 : params.addParamNamesToGroup("gap_geometry_type cylinder_axis_point_1 cylinder_axis_point_2",
116 : "Gap geometry");
117 9156 : params.addParamNamesToGroup("emissivity_primary emissivity_secondary", "Radiative heat transfer");
118 9156 : params.addParamNamesToGroup("min_gap max_gap min_gap_order", "Gap size");
119 :
120 4578 : return params;
121 4578 : }
122 :
123 893 : GapConductance::GapConductance(const InputParameters & parameters)
124 : : Material(parameters),
125 893 : _appended_property_name(getParam<std::string>("appended_property_name")),
126 893 : _temp(coupledValue("variable")),
127 893 : _gap_geometry_type(declareRestartableData<GapConductance::GAP_GEOMETRY>("gap_geometry_type",
128 893 : GapConductance::PLATE)),
129 1786 : _quadrature(getParam<bool>("quadrature")),
130 893 : _gap_temp(0),
131 893 : _gap_distance(88888),
132 893 : _radius(0),
133 893 : _r1(0),
134 893 : _r2(0),
135 893 : _has_info(false),
136 893 : _gap_distance_value(_quadrature ? _zero : coupledValue("gap_distance")),
137 893 : _gap_temp_value(_quadrature ? _zero : coupledValue("gap_temp")),
138 893 : _gap_conductance(declareProperty<Real>("gap_conductance" + _appended_property_name)),
139 1786 : _gap_conductance_dT(declareProperty<Real>("gap_conductance" + _appended_property_name + "_dT")),
140 893 : _gap_thermal_conductivity(declareProperty<Real>("gap_conductivity")),
141 1786 : _gap_conductivity(getParam<Real>("gap_conductivity")),
142 893 : _gap_conductivity_function(isParamValid("gap_conductivity_function")
143 893 : ? &getFunction("gap_conductivity_function")
144 : : nullptr),
145 893 : _gap_conductivity_function_variable(isCoupled("gap_conductivity_function_variable")
146 893 : ? &coupledValue("gap_conductivity_function_variable")
147 : : nullptr),
148 1786 : _stefan_boltzmann(getParam<Real>("stefan_boltzmann")),
149 1786 : _emissivity_primary(getParam<Real>("emissivity_primary")),
150 1786 : _emissivity_secondary(getParam<Real>("emissivity_secondary")),
151 1786 : _min_gap(getParam<Real>("min_gap")),
152 1786 : _min_gap_order(getParam<unsigned int>("min_gap_order")),
153 1786 : _max_gap(getParam<Real>("max_gap")),
154 893 : _temp_var(_quadrature ? getVar("variable", 0) : nullptr),
155 893 : _penetration_locator(nullptr),
156 893 : _serialized_solution(_quadrature ? &_temp_var->sys().currentSolution() : nullptr),
157 893 : _dof_map(_quadrature ? &_temp_var->sys().dofMap() : nullptr),
158 1786 : _warnings(getParam<bool>("warnings")),
159 1786 : _p1(declareRestartableData<Point>("cylinder_axis_point_1", Point(0, 1, 0))),
160 2679 : _p2(declareRestartableData<Point>("cylinder_axis_point_2", Point(0, 0, 0)))
161 : {
162 170 : _emissivity = _emissivity_primary != 0.0 && _emissivity_secondary != 0.0
163 1063 : ? 1.0 / _emissivity_primary + 1.0 / _emissivity_secondary - 1
164 : : 0.0;
165 :
166 893 : if (_quadrature)
167 : {
168 1380 : if (!parameters.isParamValid("paired_boundary"))
169 0 : mooseError("No 'paired_boundary' provided for ", _name);
170 :
171 2070 : _penetration_locator = &_subproblem.geomSearchData().getQuadraturePenetrationLocator(
172 : parameters.get<BoundaryName>("paired_boundary"),
173 1380 : getParam<std::vector<BoundaryName>>("boundary")[0],
174 1380 : Utility::string_to_enum<Order>(parameters.get<MooseEnum>("order")));
175 : }
176 : else
177 : {
178 406 : if (!isCoupled("gap_distance"))
179 0 : paramError("gap_distance", "needed if not using quadrature point based gap heat");
180 :
181 406 : if (!isCoupled("gap_temp"))
182 0 : paramError("gap_temp", "needed if not using quadrature point based gap heat");
183 : }
184 :
185 893 : if (_mesh.uniformRefineLevel() != 0)
186 2 : mooseError("GapConductance does not work with uniform mesh refinement.");
187 891 : }
188 :
189 : void
190 891 : GapConductance::initialSetup()
191 : {
192 : ///set generated from the passed in vector of subdomain names
193 : const auto & check_subdomains =
194 891 : blockRestricted() && !blockIDs().empty() ? blockIDs() : meshBlockIDs();
195 891 : if (check_subdomains.empty())
196 0 : mooseError("No subdomains found");
197 :
198 : // make sure all subdomains are using the same coordinate system
199 891 : Moose::CoordinateSystemType coord_system = _fe_problem.getCoordSystem(*check_subdomains.begin());
200 2683 : for (auto subdomain : check_subdomains)
201 1792 : if (_fe_problem.getCoordSystem(subdomain) != coord_system)
202 0 : mooseError(
203 : "The GapConductance model requires all subdomains to have the same coordinate system.");
204 :
205 : // Select proper coordinate system and geometry (plate, cylinder, spheres)
206 891 : setGapGeometryParameters(
207 891 : _pars, coord_system, _fe_problem.getAxisymmetricRadialCoord(), _gap_geometry_type, _p1, _p2);
208 891 : }
209 :
210 : void
211 1832 : GapConductance::setGapGeometryParameters(const InputParameters & params,
212 : const Moose::CoordinateSystemType coord_sys,
213 : unsigned int axisymmetric_radial_coord,
214 : GAP_GEOMETRY & gap_geometry_type,
215 : Point & p1,
216 : Point & p2)
217 : {
218 3664 : if (params.isParamSetByUser("gap_geometry_type"))
219 : {
220 648 : gap_geometry_type =
221 648 : GapConductance::GAP_GEOMETRY(int(params.get<MooseEnum>("gap_geometry_type")));
222 : }
223 : else
224 : {
225 1184 : if (coord_sys == Moose::COORD_XYZ)
226 892 : gap_geometry_type = GapConductance::PLATE;
227 292 : else if (coord_sys == Moose::COORD_RZ)
228 292 : gap_geometry_type = GapConductance::CYLINDER;
229 0 : else if (coord_sys == Moose::COORD_RSPHERICAL)
230 0 : gap_geometry_type = GapConductance::SPHERE;
231 : }
232 :
233 1832 : if (gap_geometry_type == GapConductance::PLATE)
234 : {
235 964 : if (coord_sys == Moose::COORD_RSPHERICAL)
236 0 : ::mooseError("'gap_geometry_type = PLATE' cannot be used with models having a spherical "
237 : "coordinate system.");
238 : }
239 868 : else if (gap_geometry_type == GapConductance::CYLINDER)
240 : {
241 652 : if (coord_sys == Moose::COORD_XYZ)
242 : {
243 720 : if (!params.isParamValid("cylinder_axis_point_1") ||
244 720 : !params.isParamValid("cylinder_axis_point_2"))
245 0 : ::mooseError("For 'gap_geometry_type = CYLINDER' to be used with a Cartesian model, "
246 : "'cylinder_axis_point_1' and 'cylinder_axis_point_2' must be specified.");
247 240 : p1 = params.get<RealVectorValue>("cylinder_axis_point_1");
248 240 : p2 = params.get<RealVectorValue>("cylinder_axis_point_2");
249 : }
250 412 : else if (coord_sys == Moose::COORD_RZ)
251 : {
252 824 : if (params.isParamValid("cylinder_axis_point_1") ||
253 824 : params.isParamValid("cylinder_axis_point_2"))
254 0 : ::mooseError("The 'cylinder_axis_point_1' and 'cylinder_axis_point_2' cannot be specified "
255 : "with axisymmetric models. The y-axis is used as the cylindrical axis of "
256 : "symmetry.");
257 :
258 412 : if (axisymmetric_radial_coord == 0) // R-Z problem
259 : {
260 392 : p1 = Point(0, 0, 0);
261 392 : p2 = Point(0, 1, 0);
262 : }
263 : else // Z-R problem
264 : {
265 20 : p1 = Point(0, 0, 0);
266 20 : p2 = Point(1, 0, 0);
267 : }
268 : }
269 0 : else if (coord_sys == Moose::COORD_RSPHERICAL)
270 0 : ::mooseError("'gap_geometry_type = CYLINDER' cannot be used with models having a spherical "
271 : "coordinate system.");
272 : }
273 216 : else if (gap_geometry_type == GapConductance::SPHERE)
274 : {
275 216 : if (coord_sys == Moose::COORD_XYZ || coord_sys == Moose::COORD_RZ)
276 : {
277 312 : if (!params.isParamValid("sphere_origin"))
278 0 : ::mooseError("For 'gap_geometry_type = SPHERE' to be used with a Cartesian or axisymmetric "
279 : "model, 'sphere_origin' must be specified.");
280 156 : p1 = params.get<RealVectorValue>("sphere_origin");
281 : }
282 60 : else if (coord_sys == Moose::COORD_RSPHERICAL)
283 : {
284 120 : if (params.isParamValid("sphere_origin"))
285 0 : ::mooseError("The 'sphere_origin' cannot be specified with spherical models. x=0 is used "
286 : "as the spherical origin.");
287 60 : p1 = Point(0, 0, 0);
288 : }
289 : }
290 1832 : }
291 :
292 : void
293 2103646 : GapConductance::computeQpProperties()
294 : {
295 2103646 : computeGapValues();
296 2103646 : computeQpConductance();
297 2103646 : }
298 :
299 : void
300 2103646 : GapConductance::computeQpConductance()
301 : {
302 2103646 : if (_has_info)
303 : {
304 1935352 : _gap_conductance[_qp] = h_conduction() + h_radiation();
305 1935352 : _gap_conductance_dT[_qp] = dh_conduction() + dh_radiation();
306 : }
307 : else
308 : {
309 168294 : _gap_conductance[_qp] = 0;
310 168294 : _gap_conductance_dT[_qp] = 0;
311 : }
312 2103646 : }
313 :
314 : Real
315 2045944 : GapConductance::gapAttenuation(Real adjusted_length, Real min_gap, unsigned int min_gap_order)
316 : {
317 : mooseAssert(min_gap > 0, "min_gap must be larger than zero.");
318 :
319 2045944 : if (adjusted_length > min_gap)
320 2025034 : return 1.0 / adjusted_length;
321 : else
322 20910 : switch (min_gap_order)
323 : {
324 10374 : case 0:
325 10374 : return 1.0 / min_gap;
326 :
327 10536 : case 1:
328 10536 : return 1.0 / min_gap - (adjusted_length - min_gap) / (min_gap * min_gap);
329 :
330 0 : default:
331 0 : ::mooseError("Invalid Taylor expansion order");
332 : }
333 : }
334 :
335 : Real
336 1935352 : GapConductance::h_conduction()
337 : {
338 1935352 : _gap_thermal_conductivity[_qp] = gapK();
339 1935352 : const Real adjusted_length = gapLength(_gap_geometry_type, _radius, _r1, _r2, _max_gap);
340 1935352 : return _gap_thermal_conductivity[_qp] * gapAttenuation(adjusted_length, _min_gap, _min_gap_order);
341 : }
342 :
343 : Real
344 1935352 : GapConductance::dh_conduction()
345 : {
346 1935352 : return 0.0;
347 : }
348 :
349 : Real
350 1935352 : GapConductance::h_radiation()
351 : {
352 : /*
353 : Gap conductance due to radiation is based on the diffusion approximation:
354 :
355 : q12 = sigma*Fe*(T1^4 - T2^4) ~ hr(T1 - T2)
356 : where sigma is the Stefan-Boltzmann constant, Fe is an emissivity
357 : function, T1 and T2 are the temperatures of the two surfaces, and
358 : hr is the radiant gap conductance. Solving for hr,
359 :
360 : hr = sigma*Fe*(T1^4 - T2^4) / (T1 - T2)
361 : which can be factored to give:
362 :
363 : hr = sigma*Fe*(T1^2 + T2^2) * (T1 + T2)
364 :
365 : Assuming the gap is between infinite parallel planes, the emissivity
366 : function is given by:
367 :
368 : Fe = 1 / (1/e1 + 1/e2 - 1)
369 :
370 : For cylinders and spheres, see Fundamentals of Heat and Mass Transfer,
371 : Sixth Edition, John Wiley & Sons, Table 13.3.
372 :
373 : For cylinders:
374 :
375 : Fe = 1 / (1/e1 + (1/e2 - 1) * (r1/r2))
376 :
377 : q21 = -q12 * (r1/r2)
378 :
379 : For spheres:
380 :
381 : Fe = 1 / (1/e1 + (1/e2 - 1) * (r1/r2)^2)
382 :
383 : q21 = -q12 * (r1/r2)^2
384 : */
385 :
386 1935352 : if (_emissivity == 0.0)
387 : return 0.0;
388 :
389 : // We add 'surface_integration_factor' to account for the surface integration of the conductance
390 : // due to radiation.
391 : Real surface_integration_factor = 1.0;
392 :
393 617056 : if (_gap_geometry_type == GapConductance::CYLINDER)
394 : {
395 588768 : _emissivity = 1.0 / _emissivity_primary + (1.0 / _emissivity_secondary - 1) * _r1 / _r2;
396 588768 : if (_r2 == _radius)
397 294384 : surface_integration_factor = _r1 / _r2;
398 : }
399 28288 : else if (_gap_geometry_type == GapConductance::SPHERE)
400 : {
401 23976 : _emissivity =
402 23976 : 1.0 / _emissivity_primary + (1.0 / _emissivity_secondary - 1) * _r1 * _r1 / (_r2 * _r2);
403 23976 : if (_r2 == _radius)
404 8400 : surface_integration_factor = _r1 * _r1 / (_r2 * _r2);
405 : }
406 :
407 : const Real temp_func =
408 617056 : (_temp[_qp] * _temp[_qp] + _gap_temp * _gap_temp) * (_temp[_qp] + _gap_temp);
409 :
410 617056 : return _stefan_boltzmann * temp_func / _emissivity * surface_integration_factor;
411 : }
412 :
413 : Real
414 1935352 : GapConductance::dh_radiation()
415 : {
416 1935352 : if (_emissivity == 0.0)
417 : return 0.0;
418 :
419 : Real surface_integration_factor = 1.0;
420 :
421 617056 : if (_gap_geometry_type == GapConductance::CYLINDER)
422 : {
423 588768 : _emissivity = 1.0 / _emissivity_primary + (1.0 / _emissivity_secondary - 1) * _r1 / _r2;
424 588768 : if (_r2 == _radius)
425 294384 : surface_integration_factor = _r1 / _r2;
426 : }
427 28288 : else if (_gap_geometry_type == GapConductance::SPHERE)
428 : {
429 23976 : _emissivity =
430 23976 : 1.0 / _emissivity_primary + (1.0 / _emissivity_secondary - 1) * _r1 * _r1 / (_r2 * _r2);
431 23976 : if (_r2 == _radius)
432 8400 : surface_integration_factor = _r1 * _r1 / (_r2 * _r2);
433 : }
434 :
435 617056 : const Real temp_func = 3 * _temp[_qp] * _temp[_qp] + _gap_temp * (2 * _temp[_qp] + _gap_temp);
436 :
437 617056 : return _stefan_boltzmann * temp_func / _emissivity * surface_integration_factor;
438 : }
439 :
440 : Real
441 2156536 : GapConductance::gapLength(const GapConductance::GAP_GEOMETRY & gap_geom,
442 : const Real radius,
443 : const Real r1,
444 : const Real r2,
445 : const Real max_gap)
446 : {
447 2156536 : if (gap_geom == GapConductance::CYLINDER)
448 832874 : return gapCyl(radius, r1, r2, max_gap);
449 1323662 : else if (gap_geom == GapConductance::SPHERE)
450 489421 : return gapSphere(radius, r1, r2, max_gap);
451 : else
452 834241 : return gapRect(r2 - r1, max_gap);
453 : }
454 :
455 : Real
456 834241 : GapConductance::gapRect(const Real distance, const Real max_gap)
457 : {
458 834241 : return std::min(distance, max_gap);
459 : }
460 :
461 : Real
462 832874 : GapConductance::gapCyl(const Real radius, const Real r1, const Real r2, const Real max_denom)
463 : {
464 832874 : const Real denominator = radius * std::log(r2 / r1);
465 832874 : return std::min(denominator, max_denom);
466 : }
467 :
468 : Real
469 489421 : GapConductance::gapSphere(const Real radius, const Real r1, const Real r2, const Real max_denom)
470 : {
471 489421 : const Real denominator = radius * radius * ((1.0 / r1) - (1.0 / r2));
472 489421 : return std::min(denominator, max_denom);
473 : }
474 :
475 : Real
476 1935352 : GapConductance::gapK()
477 : {
478 1935352 : Real gap_conductivity = _gap_conductivity;
479 :
480 1935352 : if (_gap_conductivity_function)
481 : {
482 0 : if (_gap_conductivity_function_variable)
483 0 : gap_conductivity *= _gap_conductivity_function->value(
484 0 : (*_gap_conductivity_function_variable)[_qp], _q_point[_qp]);
485 : else
486 0 : gap_conductivity *= _gap_conductivity_function->value(_t, _q_point[_qp]);
487 : }
488 :
489 1935352 : return gap_conductivity;
490 : }
491 :
492 : void
493 2103646 : GapConductance::computeGapValues()
494 : {
495 2103646 : if (!_quadrature)
496 : {
497 729389 : _has_info = true;
498 729389 : _gap_temp = _gap_temp_value[_qp];
499 729389 : _gap_distance = _gap_distance_value[_qp];
500 : }
501 : else
502 : {
503 1374257 : Node * qnode = _mesh.getQuadratureNode(_current_elem, _current_side, _qp);
504 1374257 : PenetrationInfo * pinfo = _penetration_locator->_penetration_info[qnode->id()];
505 :
506 1374257 : _gap_temp = 0.0;
507 1374257 : _gap_distance = 88888;
508 1374257 : _has_info = false;
509 :
510 1374257 : if (pinfo)
511 : {
512 1205963 : _gap_distance = pinfo->_distance;
513 1205963 : _has_info = true;
514 :
515 1205963 : const Elem * secondary_side = pinfo->_side;
516 : std::vector<std::vector<Real>> & secondary_side_phi = pinfo->_side_phi;
517 : std::vector<dof_id_type> secondary_side_dof_indices;
518 :
519 1205963 : _dof_map->dof_indices(secondary_side, secondary_side_dof_indices, _temp_var->number());
520 :
521 5758586 : for (unsigned int i = 0; i < secondary_side_dof_indices.size(); ++i)
522 : {
523 : // The zero index is because we only have one point that the phis are evaluated at
524 4552623 : _gap_temp +=
525 4552623 : secondary_side_phi[i][0] * (*(*_serialized_solution))(secondary_side_dof_indices[i]);
526 : }
527 1205963 : }
528 : else
529 : {
530 168294 : if (_warnings)
531 0 : mooseWarning("No gap value information found for node ",
532 0 : qnode->id(),
533 : " on processor ",
534 0 : processor_id(),
535 : " at coordinate ",
536 0 : Point(*qnode));
537 : }
538 : }
539 :
540 2103646 : Point current_point(_q_point[_qp]);
541 2103646 : computeGapRadii(
542 2103646 : _gap_geometry_type, current_point, _p1, _p2, _gap_distance, _normals[_qp], _r1, _r2, _radius);
543 2103646 : }
544 :
545 : void
546 17398780 : GapConductance::computeGapRadii(const GapConductance::GAP_GEOMETRY gap_geometry_type,
547 : const Point & current_point,
548 : const Point & p1,
549 : const Point & p2,
550 : const Real & gap_distance,
551 : const Point & current_normal,
552 : Real & r1,
553 : Real & r2,
554 : Real & radius)
555 : {
556 17398780 : if (gap_geometry_type == GapConductance::CYLINDER)
557 : {
558 : // The vector _p1 + t*(_p2-_p1) defines the cylindrical axis. The point along this
559 : // axis closest to current_point is found by the following for t:
560 : const Point p2p1(p2 - p1);
561 : const Point p1pc(p1 - current_point);
562 4923842 : const Real t = -(p1pc * p2p1) / p2p1.norm_sq();
563 :
564 : // The nearest point on the cylindrical axis to current_point is p.
565 : const Point p(p1 + t * p2p1);
566 : Point rad_vec(current_point - p);
567 4923842 : Real rad = rad_vec.norm();
568 : rad_vec /= rad;
569 : Real rad_dot_norm = rad_vec * current_normal;
570 :
571 4923842 : if (rad_dot_norm > 0)
572 : {
573 2615484 : r1 = rad;
574 2615484 : r2 = rad - gap_distance; // note, gap_distance is negative
575 2615484 : radius = r1;
576 : }
577 2308358 : else if (rad_dot_norm < 0)
578 : {
579 2308358 : r1 = rad + gap_distance;
580 2308358 : r2 = rad;
581 2308358 : radius = r2;
582 : }
583 : else
584 0 : ::mooseError("Issue with cylindrical flux calc. normals.\n");
585 : }
586 12474938 : else if (gap_geometry_type == GapConductance::SPHERE)
587 : {
588 : const Point origin_to_curr_point(current_point - p1);
589 : const Real normal_dot = origin_to_curr_point * current_normal;
590 4813378 : const Real curr_point_radius = origin_to_curr_point.norm();
591 4813378 : if (normal_dot > 0) // on inside surface
592 : {
593 4325333 : r1 = curr_point_radius;
594 4325333 : r2 = curr_point_radius - gap_distance; // gap_distance is negative
595 4325333 : radius = r1;
596 : }
597 488045 : else if (normal_dot < 0) // on outside surface
598 : {
599 488045 : r1 = curr_point_radius + gap_distance; // gap_distance is negative
600 488045 : r2 = curr_point_radius;
601 488045 : radius = r2;
602 : }
603 : else
604 0 : ::mooseError("Issue with spherical flux calc. normals. \n");
605 : }
606 : else
607 : {
608 7661560 : r2 = -gap_distance;
609 7661560 : r1 = 0;
610 7661560 : radius = 0;
611 : }
612 17398780 : }
|
