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 "ExplicitDynamicsContactAction.h"
11 :
12 : #include "Factory.h"
13 : #include "FEProblem.h"
14 : #include "Conversion.h"
15 : #include "AddVariableAction.h"
16 : #include "NonlinearSystemBase.h"
17 : #include "Parser.h"
18 :
19 : #include "NanoflannMeshAdaptor.h"
20 : #include "PointListAdaptor.h"
21 :
22 : #include <set>
23 : #include <algorithm>
24 : #include <unordered_map>
25 : #include <limits>
26 :
27 : #include "libmesh/petsc_nonlinear_solver.h"
28 : #include "libmesh/string_to_enum.h"
29 :
30 : // Counter for distinct contact action objects
31 : static unsigned int ed_contact_action_counter = 0;
32 :
33 : registerMooseAction("ContactApp", ExplicitDynamicsContactAction, "add_aux_variable");
34 : registerMooseAction("ContactApp", ExplicitDynamicsContactAction, "add_contact_aux_variable");
35 : registerMooseAction("ContactApp", ExplicitDynamicsContactAction, "add_aux_kernel");
36 :
37 : registerMooseAction("ContactApp", ExplicitDynamicsContactAction, "add_constraint");
38 : registerMooseAction("ContactApp", ExplicitDynamicsContactAction, "add_user_object");
39 :
40 : InputParameters
41 90 : ExplicitDynamicsContactAction::validParams()
42 : {
43 90 : InputParameters params = Action::validParams();
44 90 : params += ExplicitDynamicsContactAction::commonParameters();
45 :
46 180 : params.addParam<std::vector<BoundaryName>>(
47 : "primary", "The list of boundary IDs referring to primary sidesets");
48 180 : params.addParam<std::vector<BoundaryName>>(
49 : "secondary", "The list of boundary IDs referring to secondary sidesets");
50 180 : params.addParam<std::vector<VariableName>>(
51 : "displacements",
52 : "The displacements appropriate for the simulation geometry and coordinate system");
53 180 : params.addParam<Real>("friction_coefficient", 0, "The friction coefficient");
54 180 : params.addParam<MooseEnum>(
55 180 : "model", ExplicitDynamicsContactAction::getModelEnum(), "The contact model to use");
56 180 : params.addParam<Real>("tangential_tolerance",
57 : "Tangential distance to extend edges of contact surfaces");
58 180 : params.addParam<Real>("penalty", 1e8, "Penalty factor for normal contact.");
59 90 : params.addClassDescription("Sets up all objects needed for mechanical contact enforcement in "
60 : "explicit dynamics simulations.");
61 180 : params.addParam<std::vector<TagName>>(
62 : "extra_vector_tags",
63 : "The tag names for extra vectors that residual data should be saved into");
64 180 : params.addParam<std::vector<TagName>>(
65 : "absolute_value_vector_tags",
66 : "The tags for the vectors this residual object should fill with the "
67 : "absolute value of the residual contribution");
68 :
69 180 : params.addParam<VariableName>("secondary_gap_offset",
70 : "Offset to gap distance from secondary side");
71 180 : params.addParam<VariableName>("mapped_primary_gap_offset",
72 : "Offset to gap distance mapped from primary side");
73 180 : params.addParam<bool>("verbose", false, "Verbose output. May increase runtime.");
74 :
75 90 : return params;
76 0 : }
77 :
78 90 : ExplicitDynamicsContactAction::ExplicitDynamicsContactAction(const InputParameters & params)
79 : : Action(params),
80 180 : _boundary_pairs(getParam<BoundaryName, BoundaryName>("primary", "secondary")),
81 180 : _model(getParam<MooseEnum>("model").getEnum<ExplicitDynamicsContactModel>()),
82 270 : _verbose(getParam<bool>("verbose"))
83 : {
84 : // The resulting velocity of the contact algorithm is applied, as the code stands, by modifying
85 : // the old position. This causes artifacts in the internal forces as old position, new position,
86 : // and velocities are not fully consistent. Results improve with shorter time steps, but it would
87 : // be best to modify all the solution arrays for consistency or replace the way the explicit
88 : // system is solved to obtain accelerations (a = M^{-1}F) such that only velocities on the
89 : // contacting boundary need to be updated with the contact algorithm output.
90 90 : mooseWarning("Verification of explicit dynamics capabilities is an ongoing effort.");
91 90 : }
92 :
93 : void
94 450 : ExplicitDynamicsContactAction::act()
95 : {
96 : // proform problem checks/corrections once during the first feasible task
97 450 : if (_current_task == "add_contact_aux_variable")
98 : {
99 180 : if (!_problem->getDisplacedProblem())
100 0 : mooseError(
101 : "Contact requires updated coordinates. Use the 'displacements = ...' parameter in the "
102 : "Mesh block.");
103 :
104 : // It is risky to apply this optimization to contact problems
105 : // since the problem configuration may be changed during Jacobian
106 : // evaluation. We therefore turn it off for all contact problems so that
107 : // PETSc-3.8.4 or higher will have the same behavior as PETSc-3.8.3.
108 90 : if (!_problem->isSNESMFReuseBaseSetbyUser())
109 : _problem->setSNESMFReuseBase(false, false);
110 : }
111 :
112 450 : addNodeFaceContact();
113 :
114 450 : if (_current_task == "add_aux_kernel" && _verbose)
115 : { // Add ContactPenetrationAuxAction.
116 144 : if (!_problem->getDisplacedProblem())
117 0 : mooseError("Contact requires updated coordinates. Use the 'displacements = ...' line in the "
118 : "Mesh block.");
119 : unsigned int pair_number(0);
120 : // Create auxiliary kernels for each contact pairs
121 144 : for (const auto & contact_pair : _boundary_pairs)
122 : {
123 : {
124 72 : InputParameters params = _factory.getValidParams("PenetrationAux");
125 288 : params.applyParameters(parameters(),
126 : {"secondary_gap_offset", "mapped_primary_gap_offset", "order"});
127 :
128 : std::vector<VariableName> displacements =
129 72 : getParam<std::vector<VariableName>>("displacements");
130 72 : const auto order = _problem->systemBaseNonlinear(/*nl_sys_num=*/0)
131 72 : .system()
132 72 : .variable_type(displacements[0])
133 : .order.get_order();
134 :
135 144 : params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
136 288 : params.set<ExecFlagEnum>("execute_on") = {EXEC_INITIAL, EXEC_LINEAR};
137 216 : params.set<std::vector<BoundaryName>>("boundary") = {contact_pair.second};
138 144 : params.set<BoundaryName>("paired_boundary") = contact_pair.first;
139 144 : params.set<AuxVariableName>("variable") = "penetration";
140 144 : if (isParamValid("secondary_gap_offset"))
141 0 : params.set<std::vector<VariableName>>("secondary_gap_offset") = {
142 0 : getParam<VariableName>("secondary_gap_offset")};
143 144 : if (isParamValid("mapped_primary_gap_offset"))
144 0 : params.set<std::vector<VariableName>>("mapped_primary_gap_offset") = {
145 0 : getParam<VariableName>("mapped_primary_gap_offset")};
146 72 : params.set<bool>("use_displaced_mesh") = true;
147 144 : std::string name = _name + "_contact_" + Moose::stringify(pair_number);
148 72 : pair_number++;
149 144 : _problem->addAuxKernel("PenetrationAux", name, params);
150 72 : }
151 : }
152 :
153 72 : addContactPressureAuxKernel();
154 : }
155 :
156 450 : if (_current_task == "add_contact_aux_variable")
157 : {
158 180 : std::vector<VariableName> displacements = getParam<std::vector<VariableName>>("displacements");
159 90 : const auto order = _problem->systemBaseNonlinear(/*nl_sys_num=*/0)
160 90 : .system()
161 90 : .variable_type(displacements[0])
162 : .order.get_order();
163 :
164 : // Add gap rate for output
165 : {
166 90 : auto var_params = _factory.getValidParams("MooseVariable");
167 180 : var_params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
168 180 : var_params.set<MooseEnum>("family") = "LAGRANGE";
169 180 : _problem->addAuxVariable("MooseVariable", "gap_rate", var_params);
170 90 : }
171 90 : if (_verbose)
172 : {
173 : // Add ContactPenetrationVarAction
174 : {
175 72 : auto var_params = _factory.getValidParams("MooseVariable");
176 144 : var_params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
177 144 : var_params.set<MooseEnum>("family") = "LAGRANGE";
178 :
179 144 : _problem->addAuxVariable("MooseVariable", "penetration", var_params);
180 72 : }
181 : {
182 72 : auto var_params = _factory.getValidParams("MooseVariable");
183 144 : var_params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
184 144 : var_params.set<MooseEnum>("family") = "LAGRANGE";
185 :
186 144 : _problem->addAuxVariable("MooseVariable", "contact_pressure", var_params);
187 72 : }
188 : // Add nodal area contact variable
189 : {
190 72 : auto var_params = _factory.getValidParams("MooseVariable");
191 144 : var_params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
192 144 : var_params.set<MooseEnum>("family") = "LAGRANGE";
193 :
194 144 : _problem->addAuxVariable("MooseVariable", "nodal_area", var_params);
195 72 : }
196 : // Add nodal density variable
197 : {
198 72 : auto var_params = _factory.getValidParams("MooseVariable");
199 144 : var_params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
200 144 : var_params.set<MooseEnum>("family") = "LAGRANGE";
201 :
202 144 : _problem->addAuxVariable("MooseVariable", "nodal_density", var_params);
203 72 : }
204 : // Add nodal wave speed
205 : {
206 72 : auto var_params = _factory.getValidParams("MooseVariable");
207 144 : var_params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
208 144 : var_params.set<MooseEnum>("family") = "LAGRANGE";
209 :
210 144 : _problem->addAuxVariable("MooseVariable", "nodal_wave_speed", var_params);
211 72 : }
212 : }
213 90 : }
214 :
215 450 : if (_current_task == "add_user_object" && _verbose)
216 : {
217 : {
218 72 : auto var_params = _factory.getValidParams("NodalArea");
219 :
220 : // Get secondary_boundary_vector from possibly updated set from the
221 : // ContactAction constructor cleanup
222 72 : const auto actions = _awh.getActions<ExplicitDynamicsContactAction>();
223 :
224 : std::vector<BoundaryName> secondary_boundary_vector;
225 144 : for (const auto * const action : actions)
226 144 : for (const auto j : index_range(action->_boundary_pairs))
227 72 : secondary_boundary_vector.push_back(action->_boundary_pairs[j].second);
228 :
229 144 : var_params.set<std::vector<BoundaryName>>("boundary") = secondary_boundary_vector;
230 216 : var_params.set<std::vector<VariableName>>("variable") = {"nodal_area"};
231 :
232 : mooseAssert(_problem, "Problem pointer is NULL");
233 288 : var_params.set<ExecFlagEnum>("execute_on", true) = {EXEC_INITIAL, EXEC_TIMESTEP_BEGIN};
234 72 : var_params.set<bool>("use_displaced_mesh") = true;
235 :
236 216 : _problem->addUserObject(
237 72 : "NodalArea", "nodal_area_object_" + Moose::stringify(name()), var_params);
238 144 : }
239 : // Add nodal density and nodal wave speed user
240 : {
241 : // Add nodal density
242 72 : auto var_params = _factory.getValidParams("NodalDensity");
243 :
244 : // Get secondary_boundary_vector from possibly updated set from the
245 : // ContactAction constructor cleanup
246 72 : const auto actions = _awh.getActions<ExplicitDynamicsContactAction>();
247 :
248 : std::vector<BoundaryName> secondary_boundary_vector;
249 144 : for (const auto * const action : actions)
250 144 : for (const auto j : index_range(action->_boundary_pairs))
251 72 : secondary_boundary_vector.push_back(action->_boundary_pairs[j].second);
252 :
253 144 : var_params.set<std::vector<BoundaryName>>("boundary") = secondary_boundary_vector;
254 216 : var_params.set<std::vector<VariableName>>("variable") = {"nodal_density"};
255 :
256 : mooseAssert(_problem, "Problem pointer is NULL");
257 288 : var_params.set<ExecFlagEnum>("execute_on", true) = {EXEC_INITIAL, EXEC_TIMESTEP_BEGIN};
258 72 : var_params.set<bool>("use_displaced_mesh") = true;
259 :
260 216 : _problem->addUserObject(
261 72 : "NodalDensity", "nodal_density_object_" + Moose::stringify(name()), var_params);
262 144 : }
263 : {
264 : // Add wave speed
265 72 : auto var_params = _factory.getValidParams("NodalWaveSpeed");
266 :
267 : // Get secondary_boundary_vector from possibly updated set from the
268 : // ContactAction constructor cleanup
269 72 : const auto actions = _awh.getActions<ExplicitDynamicsContactAction>();
270 :
271 : std::vector<BoundaryName> secondary_boundary_vector;
272 144 : for (const auto * const action : actions)
273 144 : for (const auto j : index_range(action->_boundary_pairs))
274 72 : secondary_boundary_vector.push_back(action->_boundary_pairs[j].second);
275 :
276 144 : var_params.set<std::vector<BoundaryName>>("boundary") = secondary_boundary_vector;
277 216 : var_params.set<std::vector<VariableName>>("variable") = {"nodal_wave_speed"};
278 :
279 : mooseAssert(_problem, "Problem pointer is NULL");
280 288 : var_params.set<ExecFlagEnum>("execute_on", true) = {EXEC_INITIAL, EXEC_TIMESTEP_BEGIN};
281 72 : var_params.set<bool>("use_displaced_mesh") = true;
282 :
283 216 : _problem->addUserObject(
284 72 : "NodalWaveSpeed", "nodal_wavespeed_object_" + Moose::stringify(name()), var_params);
285 144 : }
286 : }
287 882 : }
288 :
289 : void
290 72 : ExplicitDynamicsContactAction::addContactPressureAuxKernel()
291 : {
292 : // Add ContactPressureAux: Only one object for all contact pairs
293 72 : const auto actions = _awh.getActions<ExplicitDynamicsContactAction>();
294 :
295 : // Increment counter for contact action objects
296 72 : ed_contact_action_counter++;
297 :
298 : // Add auxiliary kernel if we are the last contact action object.
299 72 : if (ed_contact_action_counter == actions.size() && _verbose)
300 : {
301 : std::vector<BoundaryName> boundary_vector;
302 : std::vector<BoundaryName> pair_boundary_vector;
303 :
304 144 : for (const auto * const action : actions)
305 144 : for (const auto j : index_range(action->_boundary_pairs))
306 : {
307 72 : boundary_vector.push_back(action->_boundary_pairs[j].second);
308 72 : pair_boundary_vector.push_back(action->_boundary_pairs[j].first);
309 : }
310 :
311 72 : InputParameters params = _factory.getValidParams("ContactPressureAux");
312 144 : params.applyParameters(parameters(), {"order"});
313 :
314 72 : std::vector<VariableName> displacements = getParam<std::vector<VariableName>>("displacements");
315 72 : const auto order = _problem->systemBaseNonlinear(/*nl_sys_num=*/0)
316 72 : .system()
317 72 : .variable_type(displacements[0])
318 : .order.get_order();
319 :
320 144 : params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
321 72 : params.set<std::vector<BoundaryName>>("boundary") = boundary_vector;
322 144 : params.set<std::vector<BoundaryName>>("paired_boundary") = pair_boundary_vector;
323 144 : params.set<AuxVariableName>("variable") = "contact_pressure";
324 144 : params.addRequiredCoupledVar("nodal_area", "The nodal area");
325 216 : params.set<std::vector<VariableName>>("nodal_area") = {"nodal_area"};
326 72 : params.set<bool>("use_displaced_mesh") = true;
327 :
328 72 : std::string name = _name + "_contact_pressure";
329 216 : params.set<ExecFlagEnum>("execute_on", true) = {EXEC_TIMESTEP_END};
330 144 : _problem->addAuxKernel("ContactPressureAux", name, params);
331 72 : }
332 288 : }
333 :
334 : void
335 450 : ExplicitDynamicsContactAction::addNodeFaceContact()
336 : {
337 450 : if (_current_task != "add_constraint")
338 360 : return;
339 :
340 90 : std::string action_name = MooseUtils::shortName(name());
341 270 : std::vector<VariableName> displacements = getParam<std::vector<VariableName>>("displacements");
342 90 : const unsigned int ndisp = displacements.size();
343 :
344 : std::string constraint_type;
345 :
346 : constraint_type = "ExplicitDynamicsContactConstraint";
347 :
348 90 : InputParameters params = _factory.getValidParams(constraint_type);
349 :
350 540 : params.applyParameters(parameters(),
351 : {"displacements",
352 : "secondary_gap_offset",
353 : "mapped_primary_gap_offset",
354 : "primary",
355 : "secondary"});
356 :
357 90 : const auto order = _problem->systemBaseNonlinear(/*nl_sys_num=*/0)
358 90 : .system()
359 90 : .variable_type(displacements[0])
360 : .order.get_order();
361 :
362 90 : params.set<std::vector<VariableName>>("displacements") = displacements;
363 90 : params.set<bool>("use_displaced_mesh") = true;
364 180 : params.set<Real>("penalty") = getParam<Real>("penalty");
365 :
366 180 : params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
367 :
368 180 : for (const auto & contact_pair : _boundary_pairs)
369 : {
370 :
371 180 : if (isParamValid("vel_x"))
372 216 : params.set<std::vector<VariableName>>("vel_x") = getParam<std::vector<VariableName>>("vel_x");
373 180 : if (isParamValid("vel_y"))
374 216 : params.set<std::vector<VariableName>>("vel_y") = getParam<std::vector<VariableName>>("vel_y");
375 180 : if (isParamValid("vel_z"))
376 216 : params.set<std::vector<VariableName>>("vel_z") = getParam<std::vector<VariableName>>("vel_z");
377 :
378 270 : params.set<std::vector<VariableName>>("gap_rate") = {"gap_rate"};
379 :
380 90 : params.set<BoundaryName>("boundary") = contact_pair.first;
381 180 : if (isParamValid("secondary_gap_offset"))
382 0 : params.set<std::vector<VariableName>>("secondary_gap_offset") = {
383 0 : getParam<VariableName>("secondary_gap_offset")};
384 180 : if (isParamValid("mapped_primary_gap_offset"))
385 0 : params.set<std::vector<VariableName>>("mapped_primary_gap_offset") = {
386 0 : getParam<VariableName>("mapped_primary_gap_offset")};
387 :
388 360 : for (unsigned int i = 0; i < ndisp; ++i)
389 : {
390 810 : std::string name = action_name + "_constraint_" + Moose::stringify(contact_pair, "_") + "_" +
391 270 : Moose::stringify(i);
392 :
393 270 : params.set<unsigned int>("component") = i;
394 :
395 270 : params.set<BoundaryName>("primary") = contact_pair.first;
396 270 : params.set<BoundaryName>("secondary") = contact_pair.second;
397 540 : params.set<NonlinearVariableName>("variable") = displacements[i];
398 810 : params.set<std::vector<VariableName>>("primary_variable") = {displacements[i]};
399 810 : params.applySpecificParameters(parameters(),
400 : {"extra_vector_tags", "absolute_value_vector_tags"});
401 270 : _problem->addConstraint(constraint_type, name, params);
402 : }
403 : }
404 900 : }
405 :
406 : MooseEnum
407 540 : ExplicitDynamicsContactAction::getModelEnum()
408 : {
409 1080 : return MooseEnum("frictionless frictionless_balance", "frictionless");
410 : }
411 :
412 : InputParameters
413 450 : ExplicitDynamicsContactAction::commonParameters()
414 : {
415 450 : InputParameters params = emptyInputParameters();
416 :
417 900 : params.addParam<MooseEnum>(
418 900 : "model", ExplicitDynamicsContactAction::getModelEnum(), "The contact model to use");
419 :
420 : // Gap rate input
421 900 : params.addCoupledVar("vel_x", "x-component of velocity.");
422 900 : params.addCoupledVar("vel_y", "y-component of velocity.");
423 900 : params.addCoupledVar("vel_z", "z-component of velocity.");
424 : // Gap rate output
425 900 : params.addCoupledVar("gap_rate", "Gap rate for output.");
426 :
427 450 : return params;
428 0 : }
|
