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 "ModularGapConductanceConstraint.h"
11 : #include "GapFluxModelBase.h"
12 : #include "GapFluxModelPressureDependentConduction.h"
13 : #include "AutomaticMortarGeneration.h"
14 :
15 : #include "MooseError.h"
16 :
17 : #include "libmesh/parallel_algebra.h"
18 :
19 : registerMooseObject("HeatTransferApp", ModularGapConductanceConstraint);
20 :
21 : InputParameters
22 1582 : ModularGapConductanceConstraint::validParams()
23 : {
24 1582 : InputParameters params = ADMortarConstraint::validParams();
25 1582 : params.addClassDescription(
26 : "Computes the residual and Jacobian contributions for the 'Lagrange Multiplier' "
27 : "implementation of the thermal contact problem. For more information, see the "
28 : "detailed description here: http://tinyurl.com/gmmhbe9");
29 3164 : params.addParam<std::vector<UserObjectName>>("gap_flux_models",
30 : "List of GapFluxModel user objects");
31 3164 : params.addCoupledVar("displacements", "Displacement variables");
32 :
33 3164 : MooseEnum gap_geometry_type("AUTO PLATE CYLINDER SPHERE", "AUTO");
34 3164 : params.addParam<MooseEnum>(
35 : "gap_geometry_type",
36 : gap_geometry_type,
37 : "Gap calculation type. The geometry type is used to compute "
38 : "gap distances and scale fluxes to ensure energy balance. If AUTO is selected, the gap "
39 : "geometry is automatically set via the mesh coordinate system.");
40 3164 : params.addRangeCheckedParam<Real>("max_gap", 1.0e6, "max_gap>=0", "A maximum gap size");
41 3164 : params.addParam<RealVectorValue>("cylinder_axis_point_1",
42 : "Start point for line defining cylindrical axis");
43 3164 : params.addParam<RealVectorValue>("cylinder_axis_point_2",
44 : "End point for line defining cylindrical axis");
45 3164 : params.addParam<RealVectorValue>("sphere_origin", "Origin for sphere geometry");
46 :
47 : // We should default use_displaced_mesh to true. If no displaced mesh exists
48 : // FEProblemBase::addConstraint will automatically correct it to false. However,
49 : // this will still prompt a call from AugmentSparsityOnInterface to get a displaced
50 : // mortar interface since object._use_displaced_mesh = true.
51 :
52 : // These parameter groups have to match the MortarGapHeatTransferAction's
53 3164 : params.addParamNamesToGroup("gap_flux_models", "Gap flux models");
54 3164 : params.addParamNamesToGroup(
55 : "gap_geometry_type max_gap cylinder_axis_point_1 cylinder_axis_point_2 sphere_origin",
56 : "Gap geometry");
57 :
58 1582 : return params;
59 1582 : }
60 :
61 637 : ModularGapConductanceConstraint::ModularGapConductanceConstraint(const InputParameters & parameters)
62 : : ADMortarConstraint(parameters),
63 1274 : _gap_flux_model_names(getParam<std::vector<UserObjectName>>("gap_flux_models")),
64 1274 : _disp_name(parameters.getVecMooseType("displacements")),
65 637 : _n_disp(_disp_name.size()),
66 637 : _disp_secondary(_n_disp),
67 637 : _disp_primary(_n_disp),
68 1274 : _gap_geometry_type(getParam<MooseEnum>("gap_geometry_type").getEnum<GapGeometry>()),
69 637 : _gap_width(0.0),
70 637 : _surface_integration_factor(1.0),
71 1274 : _p1(declareRestartableData<Point>("cylinder_axis_point_1", Point(0, 1, 0))),
72 1274 : _p2(declareRestartableData<Point>("cylinder_axis_point_2", Point(0, 0, 0))),
73 637 : _r1(0),
74 637 : _r2(0),
75 1274 : _max_gap(getParam<Real>("max_gap")),
76 637 : _adjusted_length(0.0),
77 637 : _disp_x_var(getVar("displacements", 0)),
78 637 : _disp_y_var(getVar("displacements", 1)),
79 1274 : _disp_z_var(_n_disp == 3 ? getVar("displacements", 2) : nullptr)
80 : {
81 1321 : if (_n_disp && !getParam<bool>("use_displaced_mesh"))
82 0 : paramWarning("displacements",
83 : "You are coupling displacement variables but are evaluating the gap width on the "
84 : "undisplaced mesh. This is probably a mistake.");
85 :
86 1093 : for (unsigned int i = 0; i < _n_disp; ++i)
87 : {
88 456 : auto & disp_var = _subproblem.getStandardVariable(_tid, _disp_name[i]);
89 456 : _disp_secondary[i] = &disp_var.adSln();
90 456 : _disp_primary[i] = &disp_var.adSlnNeighbor();
91 : }
92 :
93 1274 : const auto use_displaced_mesh = getParam<bool>("use_displaced_mesh");
94 :
95 1588 : for (const auto & name : _gap_flux_model_names)
96 : {
97 951 : const auto & gap_model = getUserObjectByName<GapFluxModelBase>(name);
98 :
99 : // This constraint explicitly calls the gap flux model user objects to
100 : // obtain contributions to its residuals. It therefore depends on all
101 : // variables and material properties, that these gap flux models use, to be
102 : // computed and up to date. To ensure that we collect all variable and
103 : // material property dependencies of these models and declare them as
104 : // dependencies of this constraint object. This turns an implicit, hidden
105 : // dependency into an explicit dependency that MOOSE will automatically fulfill.
106 :
107 : // pass variable dependencies through
108 : const auto & var_dependencies = gap_model.getMooseVariableDependencies();
109 2035 : for (const auto & var : var_dependencies)
110 2168 : addMooseVariableDependency(var);
111 :
112 : // pass material property dependencies through
113 951 : const auto & mat_dependencies = gap_model.getMatPropDependencies();
114 951 : _material_property_dependencies.insert(mat_dependencies.begin(), mat_dependencies.end());
115 :
116 : // ensure that the constraint and the flux models operate on the same mesh
117 951 : if (gap_model.parameters().get<bool>("use_displaced_mesh") != use_displaced_mesh)
118 0 : paramError(
119 : "use_displaced_mesh",
120 : "The gap flux model '",
121 : name,
122 : "' should operate on the same mesh (displaced/undisplaced) as the constraint object");
123 :
124 : // add gap model to list
125 951 : _gap_flux_models.push_back(&gap_model);
126 : }
127 637 : }
128 :
129 : void
130 637 : ModularGapConductanceConstraint::initialSetup()
131 : {
132 : ///set generated from the passed in vector of subdomain names
133 637 : const auto & subdomainIDs = _mesh.meshSubdomains();
134 :
135 : // make sure all subdomains are using the same coordinate system
136 637 : Moose::CoordinateSystemType coord_system = feProblem().getCoordSystem(*subdomainIDs.begin());
137 3470 : for (auto subdomain : subdomainIDs)
138 2833 : if (feProblem().getCoordSystem(subdomain) != coord_system)
139 0 : mooseError("ModularGapConductanceConstraint requires all subdomains to have the same "
140 : "coordinate system.");
141 :
142 : // Select proper coordinate system and geometry (plate, cylinder, sphere)
143 637 : setGapGeometryParameters(
144 637 : _pars, coord_system, feProblem().getAxisymmetricRadialCoord(), _gap_geometry_type);
145 635 : }
146 :
147 : void
148 637 : ModularGapConductanceConstraint::setGapGeometryParameters(
149 : const InputParameters & params,
150 : const Moose::CoordinateSystemType coord_sys,
151 : unsigned int axisymmetric_radial_coord,
152 : GapGeometry & gap_geometry_type)
153 : {
154 :
155 : // Determine what type of gap geometry we are dealing with
156 : // Either user input or from system's coordinate systems
157 637 : if (gap_geometry_type == GapGeometry::AUTO)
158 : {
159 494 : if (coord_sys == Moose::COORD_XYZ)
160 475 : gap_geometry_type = GapGeometry::PLATE;
161 19 : else if (coord_sys == Moose::COORD_RZ)
162 19 : gap_geometry_type = GapGeometry::CYLINDER;
163 0 : else if (coord_sys == Moose::COORD_RSPHERICAL)
164 0 : gap_geometry_type = GapGeometry::SPHERE;
165 : else
166 0 : mooseError("Internal Error");
167 : }
168 :
169 1274 : if (params.isParamValid("cylinder_axis_point_1") != params.isParamValid("cylinder_axis_point_2"))
170 0 : paramError(
171 : "cylinder_axis_point_1",
172 : "Either specify both `cylinder_axis_point_1` and `cylinder_axis_point_2` or neither.");
173 :
174 : // Check consistency of geometry information
175 : // Inform the user of needed input according to gap geometry (if not PLATE)
176 637 : if (gap_geometry_type == GapGeometry::PLATE)
177 : {
178 491 : if (coord_sys == Moose::COORD_RSPHERICAL)
179 0 : paramError("gap_geometry_type",
180 : "'gap_geometry_type = PLATE' cannot be used with models having a spherical "
181 : "coordinate system.");
182 : }
183 146 : else if (gap_geometry_type == GapGeometry::CYLINDER)
184 : {
185 76 : if (coord_sys == Moose::COORD_XYZ)
186 : {
187 76 : if (params.isParamValid("cylinder_axis_point_1") &&
188 57 : params.isParamValid("cylinder_axis_point_2"))
189 : {
190 19 : _p1 = params.get<RealVectorValue>("cylinder_axis_point_1");
191 19 : _p2 = params.get<RealVectorValue>("cylinder_axis_point_2");
192 : }
193 19 : else if (_mesh.dimension() == 3)
194 0 : paramError("gap_geometry_type",
195 : "For 'gap_geometry_type = CYLINDER' to be used with a Cartesian model in 3D, "
196 : "'cylinder_axis_point_1' and 'cylinder_axis_point_2' must be specified.");
197 : else
198 : {
199 19 : deduceGeometryParameters();
200 19 : mooseInfoRepeated(
201 : "ModularGapConductanceConstraint '", name(), "' deduced cylinder axis as ", _p1, _p2);
202 : }
203 : }
204 38 : else if (coord_sys == Moose::COORD_RZ)
205 : {
206 76 : if (params.isParamValid("cylinder_axis_point_1") ||
207 76 : params.isParamValid("cylinder_axis_point_2"))
208 0 : paramError("cylinder_axis_point_1",
209 : "The 'cylinder_axis_point_1' and 'cylinder_axis_point_2' cannot be specified "
210 : "with axisymmetric models. The y-axis is used as the cylindrical axis of "
211 : "symmetry.");
212 :
213 38 : if (axisymmetric_radial_coord == 0) // R-Z problem
214 : {
215 38 : _p1 = Point(0, 0, 0);
216 38 : _p2 = Point(0, 1, 0);
217 : }
218 : else // Z-R problem
219 : {
220 0 : _p1 = Point(0, 0, 0);
221 0 : _p2 = Point(1, 0, 0);
222 : }
223 : }
224 0 : else if (coord_sys == Moose::COORD_RSPHERICAL)
225 0 : paramError("gap_geometry_type",
226 : "'gap_geometry_type = CYLINDER' cannot be used with models having a spherical "
227 : "coordinate system.");
228 : }
229 70 : else if (gap_geometry_type == GapGeometry::SPHERE)
230 : {
231 70 : if (coord_sys == Moose::COORD_XYZ || coord_sys == Moose::COORD_RZ)
232 : {
233 140 : if (params.isParamValid("sphere_origin"))
234 49 : _p1 = params.get<RealVectorValue>("sphere_origin");
235 21 : else if (coord_sys == Moose::COORD_XYZ)
236 : {
237 19 : deduceGeometryParameters();
238 19 : mooseInfoRepeated(
239 : "ModularGapConductanceConstraint '", name(), "' deduced sphere origin as ", _p1);
240 : }
241 : else
242 2 : paramError("gap_geometry_type",
243 : "For 'gap_geometry_type = SPHERE' to be used with an axisymmetric "
244 : "model, 'sphere_origin' must be specified.");
245 : }
246 0 : else if (coord_sys == Moose::COORD_RSPHERICAL)
247 : {
248 0 : if (params.isParamValid("sphere_origin"))
249 0 : paramError("sphere_origin",
250 : "The 'sphere_origin' cannot be specified with spherical models. x=0 is used "
251 : "as the spherical origin.");
252 0 : _p1 = Point(0, 0, 0);
253 : }
254 : }
255 635 : }
256 :
257 : ADReal
258 26614938 : ModularGapConductanceConstraint::computeSurfaceIntegrationFactor() const
259 : {
260 :
261 26614938 : ADReal surface_integration_factor = 1.0;
262 :
263 26614938 : if (_gap_geometry_type == GapGeometry::CYLINDER)
264 269433 : surface_integration_factor = 0.5 * (_r1 + _r2) / _radius;
265 26525127 : else if (_gap_geometry_type == GapGeometry::SPHERE)
266 10096560 : surface_integration_factor = 0.25 * (_r1 + _r2) * (_r1 + _r2) / (_radius * _radius);
267 :
268 26614938 : return surface_integration_factor;
269 : }
270 :
271 : ADReal
272 26614938 : ModularGapConductanceConstraint::computeGapLength() const
273 : {
274 : using std::log, std::min;
275 :
276 26614938 : if (_gap_geometry_type == GapGeometry::CYLINDER)
277 : {
278 179622 : const auto denominator = _radius * log(_r2 / _r1);
279 89811 : return min(denominator, _max_gap);
280 : }
281 26525127 : else if (_gap_geometry_type == GapGeometry::SPHERE)
282 : {
283 7572420 : const auto denominator = _radius * _radius * ((1.0 / _r1) - (1.0 / _r2));
284 2524140 : return min(denominator, _max_gap);
285 : }
286 : else
287 48001974 : return min(_r2 - _r1, _max_gap);
288 : }
289 :
290 : void
291 26614938 : ModularGapConductanceConstraint::computeGapRadii(const ADReal & gap_length)
292 : {
293 26614938 : const Point & current_point = _q_point[_qp];
294 :
295 26614938 : if (_gap_geometry_type == GapGeometry::CYLINDER)
296 : {
297 : // The vector _p1 + t*(_p2-_p1) defines the cylindrical axis. The point along this
298 : // axis closest to current_point is found by the following for t:
299 89811 : const Point p2p1(_p2 - _p1);
300 : const Point p1pc(_p1 - current_point);
301 89811 : const Real t = -(p1pc * p2p1) / p2p1.norm_sq();
302 :
303 : // The nearest point on the cylindrical axis to current_point is p.
304 : const Point p(_p1 + t * p2p1);
305 : Point rad_vec(current_point - p);
306 89811 : Real rad = rad_vec.norm();
307 : rad_vec /= rad;
308 : Real rad_dot_norm = rad_vec * _normals[_qp];
309 :
310 89811 : if (rad_dot_norm > 0)
311 : {
312 0 : _r1 = rad;
313 0 : _r2 = rad + gap_length;
314 0 : _radius = _r1;
315 : }
316 89811 : else if (rad_dot_norm < 0)
317 : {
318 89811 : _r1 = rad - gap_length;
319 89811 : _r2 = rad;
320 89811 : _radius = _r2;
321 : }
322 : else
323 0 : mooseError("Issue with cylindrical flux calculation normals.\n");
324 : }
325 26525127 : else if (_gap_geometry_type == GapGeometry::SPHERE)
326 : {
327 2524140 : const Point origin_to_curr_point(current_point - _p1);
328 : const Real normal_dot = origin_to_curr_point * _normals[_qp];
329 2524140 : const Real curr_point_radius = origin_to_curr_point.norm();
330 2524140 : if (normal_dot > 0) // on inside surface
331 : {
332 2440620 : _r1 = curr_point_radius;
333 2440620 : _r2 = curr_point_radius + gap_length;
334 2440620 : _radius = _r1;
335 : }
336 83520 : else if (normal_dot < 0) // on outside surface
337 : {
338 83520 : _r1 = curr_point_radius - gap_length;
339 83520 : _r2 = curr_point_radius;
340 83520 : _radius = _r2;
341 : }
342 : else
343 0 : mooseError("Issue with spherical flux calculation normals. \n");
344 : }
345 : else
346 : {
347 24000987 : _r2 = gap_length;
348 24000987 : _r1 = 0;
349 24000987 : _radius = 0;
350 : }
351 26614938 : }
352 :
353 : ADReal
354 182406006 : ModularGapConductanceConstraint::computeQpResidual(Moose::MortarType mortar_type)
355 : {
356 182406006 : switch (mortar_type)
357 : {
358 77895534 : case Moose::MortarType::Primary:
359 77895534 : return _lambda[_qp] * _test_primary[_i][_qp];
360 :
361 77895534 : case Moose::MortarType::Secondary:
362 155791068 : return -_lambda[_qp] * _test_secondary[_i][_qp];
363 :
364 26614938 : case Moose::MortarType::Lower:
365 : {
366 : // we are creating an AD version of phys points primary and secondary here...
367 26614938 : ADRealVectorValue ad_phys_points_primary = _phys_points_primary[_qp];
368 26614938 : ADRealVectorValue ad_phys_points_secondary = _phys_points_secondary[_qp];
369 :
370 : // ...which uses the derivative vector of the primary and secondary displacements as
371 : // an approximation of the true phys points derivatives when the mesh is displacing
372 26614938 : if (_displaced)
373 : {
374 : // Trim interior node variable derivatives
375 : const auto & primary_ip_lowerd_map = amg().getPrimaryIpToLowerElementMap(
376 41976 : *_lower_primary_elem, *_lower_primary_elem->interior_parent(), *_lower_secondary_elem);
377 : const auto & secondary_ip_lowerd_map =
378 41976 : amg().getSecondaryIpToLowerElementMap(*_lower_secondary_elem);
379 :
380 41976 : std::array<const MooseVariable *, 3> var_array{{_disp_x_var, _disp_y_var, _disp_z_var}};
381 : std::array<ADReal, 3> primary_disp;
382 : std::array<ADReal, 3> secondary_disp;
383 :
384 125928 : for (unsigned int i = 0; i < _n_disp; ++i)
385 : {
386 83952 : primary_disp[i] = (*_disp_primary[i])[_qp];
387 83952 : secondary_disp[i] = (*_disp_secondary[i])[_qp];
388 : }
389 :
390 41976 : trimInteriorNodeDerivatives(primary_ip_lowerd_map, var_array, primary_disp, false);
391 41976 : trimInteriorNodeDerivatives(secondary_ip_lowerd_map, var_array, secondary_disp, true);
392 :
393 : // Populate quantities with trimmed derivatives
394 125928 : for (unsigned int i = 0; i < _n_disp; ++i)
395 : {
396 83952 : ad_phys_points_primary(i).derivatives() = primary_disp[i].derivatives();
397 : ad_phys_points_secondary(i).derivatives() = secondary_disp[i].derivatives();
398 : }
399 : }
400 :
401 : // compute an ADReal gap width to pass to each gap flux model
402 26614938 : _gap_width = (ad_phys_points_primary - ad_phys_points_secondary) * _normals[_qp];
403 :
404 : // First, compute radii with _gap_width
405 26614938 : computeGapRadii(_gap_width);
406 :
407 : // With radii, compute adjusted gap length for conduction
408 26614938 : _adjusted_length = computeGapLength();
409 :
410 : // Ensure energy balance for non-flat (non-PLATE) general geometries when using radiation
411 26614938 : _surface_integration_factor = computeSurfaceIntegrationFactor();
412 :
413 : // Sum up all flux contributions from all supplied gap flux models
414 26614938 : ADReal flux = 0.0;
415 79742902 : for (const auto & model : _gap_flux_models)
416 53127964 : flux += model->computeFluxInternal(*this);
417 :
418 : // The Lagrange multiplier _is_ the gap flux
419 53229876 : return (_lambda[_qp] - flux) * _test[_i][_qp];
420 : }
421 :
422 0 : default:
423 0 : return 0;
424 : }
425 : }
426 :
427 : void
428 38 : ModularGapConductanceConstraint::deduceGeometryParameters()
429 : {
430 : Point position;
431 38 : Real area = 0.0;
432 : const auto my_pid = processor_id();
433 :
434 : // build side element list as (element, side, id) tuples
435 38 : const auto bnd = _mesh.buildActiveSideList();
436 :
437 38 : std::unique_ptr<const Elem> side_ptr;
438 4598 : for (auto [elem_id, side, id] : bnd)
439 4560 : if (id == _primary_id)
440 : {
441 1520 : const auto * elem = _mesh.elemPtr(elem_id);
442 1520 : if (elem->processor_id() != my_pid)
443 560 : continue;
444 :
445 : // update side_ptr
446 960 : elem->side_ptr(side_ptr, side);
447 :
448 : // area of the (linearized) side
449 960 : const auto side_area = side_ptr->volume();
450 :
451 : // position of the side
452 960 : const auto side_position = side_ptr->true_centroid();
453 :
454 : // sum up
455 : position += side_position * side_area;
456 960 : area += side_area;
457 : }
458 :
459 : // parallel communication
460 38 : _communicator.sum(position);
461 38 : _communicator.sum(area);
462 :
463 : // set axis
464 38 : if (area == 0.0)
465 : {
466 0 : if (_gap_geometry_type == GapGeometry::CYLINDER)
467 0 : paramError("gap_geometry_type",
468 : "Unable to decuce cylinder axis automatically, please specify "
469 : "'cylinder_axis_point_1' and 'cylinder_axis_point_2'.");
470 0 : else if (_gap_geometry_type == GapGeometry::SPHERE)
471 0 : paramError("gap_geometry_type",
472 : "Unable to decuce sphere origin automatically, please specify "
473 : "'sphere_origin'.");
474 : else
475 0 : mooseError("Internal error.");
476 : }
477 :
478 38 : _p1 = position / area;
479 38 : _p2 = _p1 + Point(0, 0, 1);
480 38 : }
|
