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