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 :
275 26614938 : if (_gap_geometry_type == GapGeometry::CYLINDER)
276 : {
277 179622 : const auto denominator = _radius * std::log(_r2 / _r1);
278 89811 : return std::min(denominator, _max_gap);
279 : }
280 26525127 : else if (_gap_geometry_type == GapGeometry::SPHERE)
281 : {
282 7572420 : const auto denominator = _radius * _radius * ((1.0 / _r1) - (1.0 / _r2));
283 2524140 : return std::min(denominator, _max_gap);
284 : }
285 : else
286 48001974 : return std::min(_r2 - _r1, _max_gap);
287 : }
288 :
289 : void
290 26614938 : ModularGapConductanceConstraint::computeGapRadii(const ADReal & gap_length)
291 : {
292 26614938 : const Point & current_point = _q_point[_qp];
293 :
294 26614938 : if (_gap_geometry_type == GapGeometry::CYLINDER)
295 : {
296 : // The vector _p1 + t*(_p2-_p1) defines the cylindrical axis. The point along this
297 : // axis closest to current_point is found by the following for t:
298 89811 : const Point p2p1(_p2 - _p1);
299 : const Point p1pc(_p1 - current_point);
300 89811 : const Real t = -(p1pc * p2p1) / p2p1.norm_sq();
301 :
302 : // The nearest point on the cylindrical axis to current_point is p.
303 : const Point p(_p1 + t * p2p1);
304 : Point rad_vec(current_point - p);
305 89811 : Real rad = rad_vec.norm();
306 : rad_vec /= rad;
307 : Real rad_dot_norm = rad_vec * _normals[_qp];
308 :
309 89811 : if (rad_dot_norm > 0)
310 : {
311 0 : _r1 = rad;
312 0 : _r2 = rad + gap_length;
313 0 : _radius = _r1;
314 : }
315 89811 : else if (rad_dot_norm < 0)
316 : {
317 89811 : _r1 = rad - gap_length;
318 89811 : _r2 = rad;
319 89811 : _radius = _r2;
320 : }
321 : else
322 0 : mooseError("Issue with cylindrical flux calculation normals.\n");
323 : }
324 26525127 : else if (_gap_geometry_type == GapGeometry::SPHERE)
325 : {
326 2524140 : const Point origin_to_curr_point(current_point - _p1);
327 : const Real normal_dot = origin_to_curr_point * _normals[_qp];
328 2524140 : const Real curr_point_radius = origin_to_curr_point.norm();
329 2524140 : if (normal_dot > 0) // on inside surface
330 : {
331 2440620 : _r1 = curr_point_radius;
332 2440620 : _r2 = curr_point_radius + gap_length;
333 2440620 : _radius = _r1;
334 : }
335 83520 : else if (normal_dot < 0) // on outside surface
336 : {
337 83520 : _r1 = curr_point_radius - gap_length;
338 83520 : _r2 = curr_point_radius;
339 83520 : _radius = _r2;
340 : }
341 : else
342 0 : mooseError("Issue with spherical flux calculation normals. \n");
343 : }
344 : else
345 : {
346 24000987 : _r2 = gap_length;
347 24000987 : _r1 = 0;
348 24000987 : _radius = 0;
349 : }
350 26614938 : }
351 :
352 : ADReal
353 182406006 : ModularGapConductanceConstraint::computeQpResidual(Moose::MortarType mortar_type)
354 : {
355 182406006 : switch (mortar_type)
356 : {
357 77895534 : case Moose::MortarType::Primary:
358 77895534 : return _lambda[_qp] * _test_primary[_i][_qp];
359 :
360 77895534 : case Moose::MortarType::Secondary:
361 155791068 : return -_lambda[_qp] * _test_secondary[_i][_qp];
362 :
363 26614938 : case Moose::MortarType::Lower:
364 : {
365 : // we are creating an AD version of phys points primary and secondary here...
366 26614938 : ADRealVectorValue ad_phys_points_primary = _phys_points_primary[_qp];
367 26614938 : ADRealVectorValue ad_phys_points_secondary = _phys_points_secondary[_qp];
368 :
369 : // ...which uses the derivative vector of the primary and secondary displacements as
370 : // an approximation of the true phys points derivatives when the mesh is displacing
371 26614938 : if (_displaced)
372 : {
373 : // Trim interior node variable derivatives
374 : const auto & primary_ip_lowerd_map = amg().getPrimaryIpToLowerElementMap(
375 41976 : *_lower_primary_elem, *_lower_primary_elem->interior_parent(), *_lower_secondary_elem);
376 : const auto & secondary_ip_lowerd_map =
377 41976 : amg().getSecondaryIpToLowerElementMap(*_lower_secondary_elem);
378 :
379 41976 : std::array<const MooseVariable *, 3> var_array{{_disp_x_var, _disp_y_var, _disp_z_var}};
380 : std::array<ADReal, 3> primary_disp;
381 : std::array<ADReal, 3> secondary_disp;
382 :
383 125928 : for (unsigned int i = 0; i < _n_disp; ++i)
384 : {
385 83952 : primary_disp[i] = (*_disp_primary[i])[_qp];
386 83952 : secondary_disp[i] = (*_disp_secondary[i])[_qp];
387 : }
388 :
389 41976 : trimInteriorNodeDerivatives(primary_ip_lowerd_map, var_array, primary_disp, false);
390 41976 : trimInteriorNodeDerivatives(secondary_ip_lowerd_map, var_array, secondary_disp, true);
391 :
392 : // Populate quantities with trimmed derivatives
393 125928 : for (unsigned int i = 0; i < _n_disp; ++i)
394 : {
395 83952 : ad_phys_points_primary(i).derivatives() = primary_disp[i].derivatives();
396 : ad_phys_points_secondary(i).derivatives() = secondary_disp[i].derivatives();
397 : }
398 : }
399 :
400 : // compute an ADReal gap width to pass to each gap flux model
401 26614938 : _gap_width = (ad_phys_points_primary - ad_phys_points_secondary) * _normals[_qp];
402 :
403 : // First, compute radii with _gap_width
404 26614938 : computeGapRadii(_gap_width);
405 :
406 : // With radii, compute adjusted gap length for conduction
407 26614938 : _adjusted_length = computeGapLength();
408 :
409 : // Ensure energy balance for non-flat (non-PLATE) general geometries when using radiation
410 26614938 : _surface_integration_factor = computeSurfaceIntegrationFactor();
411 :
412 : // Sum up all flux contributions from all supplied gap flux models
413 26614938 : ADReal flux = 0.0;
414 79742902 : for (const auto & model : _gap_flux_models)
415 53127964 : flux += model->computeFluxInternal(*this);
416 :
417 : // The Lagrange multiplier _is_ the gap flux
418 53229876 : return (_lambda[_qp] - flux) * _test[_i][_qp];
419 : }
420 :
421 0 : default:
422 0 : return 0;
423 : }
424 : }
425 :
426 : void
427 38 : ModularGapConductanceConstraint::deduceGeometryParameters()
428 : {
429 : Point position;
430 38 : Real area = 0.0;
431 : const auto my_pid = processor_id();
432 :
433 : // build side element list as (element, side, id) tuples
434 38 : const auto bnd = _mesh.buildActiveSideList();
435 :
436 38 : std::unique_ptr<const Elem> side_ptr;
437 4598 : for (auto [elem_id, side, id] : bnd)
438 4560 : if (id == _primary_id)
439 : {
440 1520 : const auto * elem = _mesh.elemPtr(elem_id);
441 1520 : if (elem->processor_id() != my_pid)
442 560 : continue;
443 :
444 : // update side_ptr
445 960 : elem->side_ptr(side_ptr, side);
446 :
447 : // area of the (linearized) side
448 960 : const auto side_area = side_ptr->volume();
449 :
450 : // position of the side
451 960 : const auto side_position = side_ptr->true_centroid();
452 :
453 : // sum up
454 : position += side_position * side_area;
455 960 : area += side_area;
456 : }
457 :
458 : // parallel communication
459 38 : _communicator.sum(position);
460 38 : _communicator.sum(area);
461 :
462 : // set axis
463 38 : if (area == 0.0)
464 : {
465 0 : if (_gap_geometry_type == GapGeometry::CYLINDER)
466 0 : paramError("gap_geometry_type",
467 : "Unable to decuce cylinder axis automatically, please specify "
468 : "'cylinder_axis_point_1' and 'cylinder_axis_point_2'.");
469 0 : else if (_gap_geometry_type == GapGeometry::SPHERE)
470 0 : paramError("gap_geometry_type",
471 : "Unable to decuce sphere origin automatically, please specify "
472 : "'sphere_origin'.");
473 : else
474 0 : mooseError("Internal error.");
475 : }
476 :
477 38 : _p1 = position / area;
478 38 : _p2 = _p1 + Point(0, 0, 1);
479 38 : }
|
