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