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 "ComputeDynamicFrictionalForceLMMechanicalContact.h"
11 : #include "DisplacedProblem.h"
12 : #include "Assembly.h"
13 : #include "Function.h"
14 : #include "MortarContactUtils.h"
15 : #include "ContactFrictionUtils.h"
16 : #include "AutomaticMortarGeneration.h"
17 :
18 : #include "metaphysicl/metaphysicl_version.h"
19 : #include "metaphysicl/dualsemidynamicsparsenumberarray.h"
20 : #include "metaphysicl/parallel_dualnumber.h"
21 : #if METAPHYSICL_MAJOR_VERSION < 2
22 : #include "metaphysicl/parallel_dynamic_std_array_wrapper.h"
23 : #else
24 : #include "metaphysicl/parallel_dynamic_array_wrapper.h"
25 : #endif
26 : #include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
27 : #include "timpi/parallel_sync.h"
28 :
29 : #include <limits>
30 :
31 : registerMooseObject("ContactApp", ComputeDynamicFrictionalForceLMMechanicalContact);
32 :
33 : InputParameters
34 208 : ComputeDynamicFrictionalForceLMMechanicalContact::validParams()
35 : {
36 208 : InputParameters params = ComputeDynamicWeightedGapLMMechanicalContact::validParams();
37 208 : params.addClassDescription("Computes the tangential frictional forces for dynamic simulations");
38 416 : params.addRequiredCoupledVar("friction_lm", "The frictional Lagrange's multiplier");
39 416 : params.addCoupledVar("friction_lm_dir",
40 : "The frictional Lagrange's multiplier for an addtional direction.");
41 416 : params.addParam<FunctionName>(
42 : "function_friction",
43 : "Coupled function to evaluate friction with values from contact pressure and relative "
44 : "tangential velocities (from the previous step).");
45 416 : params.addParam<Real>("c_t", 1e0, "Numerical parameter for tangential constraints");
46 624 : params.addRangeCheckedParam<Real>(
47 : "epsilon",
48 416 : 1.0e-7,
49 : "epsilon > 0",
50 : "Minimum value of contact pressure that will trigger frictional enforcement");
51 416 : params.addRangeCheckedParam<Real>(
52 : "mu", "mu >= 0", "The friction coefficient for the Coulomb friction law");
53 416 : params.addParam<MooseEnum>("friction_coefficient_regularization",
54 416 : Moose::Contact::frictionCoefficientRegularizationOptions(),
55 : "The regularization applied to the Coulomb friction coefficient.");
56 624 : params.addRangeCheckedParam<Real>(
57 : "friction_reference_slip",
58 416 : 0.0,
59 : "friction_reference_slip >= 0",
60 : "Reference slip increment used by friction coefficient regularization.");
61 624 : params.addRangeCheckedParam<Real>(
62 : "friction_elastic_slip",
63 416 : 0.0,
64 : "friction_elastic_slip >= 0",
65 : "Tangential elastic slip distance over which the Coulomb friction bound is reached.");
66 208 : return params;
67 0 : }
68 :
69 104 : ComputeDynamicFrictionalForceLMMechanicalContact::ComputeDynamicFrictionalForceLMMechanicalContact(
70 104 : const InputParameters & parameters)
71 : : ComputeDynamicWeightedGapLMMechanicalContact(parameters),
72 208 : _c_t(getParam<Real>("c_t")),
73 104 : _secondary_x_dot(_secondary_var.adUDot()),
74 104 : _primary_x_dot(_primary_var.adUDotNeighbor()),
75 104 : _secondary_y_dot(adCoupledDot("disp_y")),
76 104 : _primary_y_dot(adCoupledNeighborValueDot("disp_y")),
77 104 : _secondary_z_dot(_has_disp_z ? &adCoupledDot("disp_z") : nullptr),
78 104 : _primary_z_dot(_has_disp_z ? &adCoupledNeighborValueDot("disp_z") : nullptr),
79 208 : _epsilon(getParam<Real>("epsilon")),
80 404 : _mu(isParamValid("mu") ? getParam<Real>("mu") : std::numeric_limits<double>::quiet_NaN()),
81 104 : _friction_coefficient_regularization(
82 104 : getParam<MooseEnum>("friction_coefficient_regularization")
83 : .getEnum<Moose::Contact::FrictionCoefficientRegularization>()),
84 208 : _friction_reference_slip(getParam<Real>("friction_reference_slip")),
85 208 : _friction_elastic_slip(getParam<Real>("friction_elastic_slip")),
86 214 : _function_friction(isParamValid("function_friction") ? &getFunction("function_friction")
87 : : nullptr),
88 208 : _has_friction_function(isParamValid("function_friction")),
89 104 : _3d(_has_disp_z)
90 : {
91 398 : if (!_has_friction_function && !isParamValid("mu"))
92 0 : paramError("mu",
93 : "A coefficient of friction needs to be provided as a constant value or via a "
94 : "function.");
95 :
96 116 : if (_has_friction_function && isParamValid("mu"))
97 0 : paramError("mu",
98 : "Either provide a constant coefficient of friction or a function defining the "
99 : "coefficient of friction. Both inputs cannot be provided simultaneously.");
100 :
101 208 : if (!getParam<bool>("use_displaced_mesh"))
102 0 : paramError("use_displaced_mesh",
103 : "'use_displaced_mesh' must be true for the "
104 : "ComputeFrictionalForceLMMechanicalContact object");
105 :
106 182 : if (_3d && !isParamValid("friction_lm_dir"))
107 0 : paramError("friction_lm_dir",
108 : "Three-dimensional mortar frictional contact simulations require an additional "
109 : "frictional Lagrange's multiplier to enforce a second tangential pressure");
110 :
111 104 : if (_friction_coefficient_regularization !=
112 26 : Moose::Contact::FrictionCoefficientRegularization::NONE &&
113 26 : _friction_reference_slip <= 0.0)
114 0 : paramError("friction_reference_slip",
115 : "A positive friction_reference_slip is required when "
116 : "friction_coefficient_regularization is not NONE.");
117 :
118 104 : _friction_vars.push_back(getVar("friction_lm", 0));
119 :
120 104 : if (_3d)
121 52 : _friction_vars.push_back(getVar("friction_lm_dir", 0));
122 :
123 104 : if (!_friction_vars[0]->isNodal())
124 0 : if (_friction_vars[0]->feType().order != static_cast<Order>(0))
125 0 : paramError(
126 : "friction_lm",
127 : "Frictional contact constraints only support elemental variables of CONSTANT order");
128 :
129 : // Request the old solution state in unison
130 104 : _sys.solutionOld();
131 104 : }
132 :
133 : void
134 565044 : ComputeDynamicFrictionalForceLMMechanicalContact::computeQpProperties()
135 : {
136 : // Compute the value of _qp_gap
137 565044 : ComputeDynamicWeightedGapLMMechanicalContact::computeQpProperties();
138 :
139 : // It appears that the relative velocity between weighted gap and this class have a sign
140 : // difference
141 565044 : _qp_tangential_velocity_nodal = -_relative_velocity * (_JxW_msm[_qp] * _coord[_qp]);
142 565044 : _qp_real_tangential_velocity_nodal = -_relative_velocity;
143 565044 : }
144 :
145 : void
146 1986576 : ComputeDynamicFrictionalForceLMMechanicalContact::computeQpIProperties()
147 : {
148 : // Get the _dof_to_weighted_gap map
149 1986576 : ComputeDynamicWeightedGapLMMechanicalContact::computeQpIProperties();
150 :
151 1986576 : const auto & nodal_tangents = amg().getNodalTangents(*_lower_secondary_elem);
152 :
153 : // Get the _dof_to_weighted_tangential_velocity map
154 : const DofObject * const dof =
155 1986576 : _friction_vars[0]->isNodal()
156 1986576 : ? static_cast<const DofObject *>(_lower_secondary_elem->node_ptr(_i))
157 0 : : static_cast<const DofObject *>(_lower_secondary_elem);
158 :
159 : _dof_to_weighted_tangential_velocity[dof][0] +=
160 3973152 : _test[_i][_qp] * _qp_tangential_velocity_nodal * nodal_tangents[0][_i];
161 :
162 : _dof_to_real_tangential_velocity[dof][0] +=
163 3973152 : _test[_i][_qp] * _qp_real_tangential_velocity_nodal * nodal_tangents[0][_i];
164 :
165 : // Get the _dof_to_weighted_tangential_velocity map for a second direction
166 1986576 : if (_3d)
167 : {
168 : _dof_to_weighted_tangential_velocity[dof][1] +=
169 3425952 : _test[_i][_qp] * _qp_tangential_velocity_nodal * nodal_tangents[1][_i];
170 :
171 : _dof_to_real_tangential_velocity[dof][1] +=
172 3425952 : _test[_i][_qp] * _qp_real_tangential_velocity_nodal * nodal_tangents[1][_i];
173 : }
174 1986576 : }
175 :
176 : void
177 6122 : ComputeDynamicFrictionalForceLMMechanicalContact::residualSetup()
178 : {
179 : // Clear both maps
180 6122 : ComputeDynamicWeightedGapLMMechanicalContact::residualSetup();
181 : _dof_to_weighted_tangential_velocity.clear();
182 : _dof_to_real_tangential_velocity.clear();
183 6122 : }
184 :
185 : void
186 317 : ComputeDynamicFrictionalForceLMMechanicalContact::timestepSetup()
187 : {
188 :
189 317 : ComputeDynamicWeightedGapLMMechanicalContact::timestepSetup();
190 :
191 : _dof_to_old_real_tangential_velocity.clear();
192 :
193 2087 : for (const auto & [dof, real_tangential_velocity] : _dof_to_real_tangential_velocity)
194 : _dof_to_old_real_tangential_velocity.emplace(
195 : dof,
196 1770 : std::array<Real, 2>{{MetaPhysicL::raw_value(real_tangential_velocity[0]),
197 : MetaPhysicL::raw_value(real_tangential_velocity[1])}});
198 317 : }
199 :
200 : void
201 451 : ComputeDynamicFrictionalForceLMMechanicalContact::post()
202 : {
203 451 : ComputeDynamicWeightedGapLMMechanicalContact::post();
204 :
205 451 : Moose::Mortar::Contact::communicateVelocities(
206 451 : _dof_to_weighted_tangential_velocity, _mesh, _nodal, _communicator, false);
207 :
208 451 : Moose::Mortar::Contact::communicateVelocities(
209 451 : _dof_to_real_tangential_velocity, _mesh, _nodal, _communicator, false);
210 :
211 : // Enforce frictional complementarity constraints
212 3430 : for (const auto & pr : _dof_to_weighted_tangential_velocity)
213 : {
214 2979 : const DofObject * const dof = pr.first;
215 :
216 2979 : if (dof->processor_id() != this->processor_id())
217 0 : continue;
218 :
219 : // Use always weighted gap for dynamic PDASS. Omit the dynamic weighted gap approach that is
220 : // used in normal contact where the discretized gap velocity is enforced if a node has
221 : // identified to be into contact.
222 :
223 2979 : _weighted_gap_ptr = &_dof_to_weighted_gap[dof].first;
224 2979 : _normalization_ptr = &_dof_to_weighted_gap[dof].second;
225 2979 : _tangential_vel_ptr[0] = &(pr.second[0]);
226 :
227 2979 : if (_3d)
228 : {
229 2979 : _tangential_vel_ptr[1] = &(pr.second[1]);
230 2979 : enforceConstraintOnDof3d(dof);
231 : }
232 : else
233 0 : enforceConstraintOnDof(dof);
234 : }
235 451 : }
236 :
237 : void
238 5671 : ComputeDynamicFrictionalForceLMMechanicalContact::incorrectEdgeDroppingPost(
239 : const std::unordered_set<const Node *> & inactive_lm_nodes)
240 : {
241 5671 : ComputeDynamicWeightedGapLMMechanicalContact::incorrectEdgeDroppingPost(inactive_lm_nodes);
242 :
243 5671 : Moose::Mortar::Contact::communicateVelocities(
244 5671 : _dof_to_weighted_tangential_velocity, _mesh, _nodal, _communicator, false);
245 :
246 5671 : Moose::Mortar::Contact::communicateVelocities(
247 5671 : _dof_to_real_tangential_velocity, _mesh, _nodal, _communicator, false);
248 :
249 : // Enforce frictional complementarity constraints
250 51352 : for (const auto & pr : _dof_to_weighted_tangential_velocity)
251 : {
252 45681 : const DofObject * const dof = pr.first;
253 :
254 : // If node inactive, skip
255 45681 : if ((inactive_lm_nodes.find(static_cast<const Node *>(dof)) != inactive_lm_nodes.end()) ||
256 45327 : (dof->processor_id() != this->processor_id()))
257 726 : continue;
258 :
259 : // Use always weighted gap for dynamic PDASS
260 44955 : _weighted_gap_ptr = &_dof_to_weighted_gap[dof].first;
261 44955 : _normalization_ptr = &_dof_to_weighted_gap[dof].second;
262 44955 : _tangential_vel_ptr[0] = &pr.second[0];
263 :
264 44955 : if (_3d)
265 : {
266 7689 : _tangential_vel_ptr[1] = &pr.second[1];
267 7689 : enforceConstraintOnDof3d(dof);
268 : }
269 : else
270 37266 : enforceConstraintOnDof(dof);
271 : }
272 5671 : }
273 :
274 : void
275 10668 : ComputeDynamicFrictionalForceLMMechanicalContact::enforceConstraintOnDof3d(
276 : const DofObject * const dof)
277 : {
278 : using std::max, std::sqrt;
279 :
280 : // Get normal LM
281 10668 : const auto normal_dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
282 10668 : const ADReal & weighted_gap = *_weighted_gap_ptr;
283 10668 : ADReal contact_pressure = (*_sys.currentSolution())(normal_dof_index);
284 : Moose::derivInsert(contact_pressure.derivatives(), normal_dof_index, 1.);
285 :
286 : // Get friction LMs
287 : std::array<const ADReal *, 2> & tangential_vel = _tangential_vel_ptr;
288 : std::array<dof_id_type, 2> friction_dof_indices;
289 : std::array<ADReal, 2> friction_lm_values;
290 :
291 : const unsigned int num_tangents = 2;
292 32004 : for (const auto i : make_range(num_tangents))
293 : {
294 21336 : friction_dof_indices[i] = dof->dof_number(_sys.number(), _friction_vars[i]->number(), 0);
295 21336 : friction_lm_values[i] = (*_sys.currentSolution())(friction_dof_indices[i]);
296 : Moose::derivInsert(friction_lm_values[i].derivatives(), friction_dof_indices[i], 1.);
297 : }
298 :
299 : // Get normalized c and c_t values (if normalization specified
300 10668 : const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
301 10668 : const Real c_t = _normalize_c ? _c_t / *_normalization_ptr : _c_t;
302 :
303 10668 : const Real contact_pressure_old = _sys.solutionOld()(normal_dof_index);
304 :
305 : // Compute the friction coefficient (constant or function)
306 : const auto & current_real_tangential_velocity =
307 10668 : libmesh_map_find(_dof_to_real_tangential_velocity, dof);
308 : const auto old_real_tangential_velocity_it = _dof_to_old_real_tangential_velocity.find(dof);
309 : const std::array<Real, 2> current_function_real_tangential_velocity{
310 : {MetaPhysicL::raw_value(current_real_tangential_velocity[0]),
311 10668 : MetaPhysicL::raw_value(current_real_tangential_velocity[1])}};
312 : const auto & function_real_tangential_velocity =
313 : old_real_tangential_velocity_it != _dof_to_old_real_tangential_velocity.end()
314 12432 : ? old_real_tangential_velocity_it->second
315 : : current_function_real_tangential_velocity;
316 : const ADReal slip_increment =
317 10668 : sqrt(current_real_tangential_velocity[0] * current_real_tangential_velocity[0] +
318 10668 : current_real_tangential_velocity[1] * current_real_tangential_velocity[1] + 1.0e-24) *
319 10668 : _dt;
320 : ADReal mu_ad = computeFrictionValue(contact_pressure_old,
321 : function_real_tangential_velocity[0],
322 : function_real_tangential_velocity[1],
323 10668 : slip_increment);
324 :
325 : ADReal dof_residual;
326 : ADReal dof_residual_dir;
327 :
328 : // Primal-dual active set strategy (PDASS)
329 10668 : if (contact_pressure < _epsilon)
330 : {
331 2899 : dof_residual = friction_lm_values[0];
332 2899 : dof_residual_dir = friction_lm_values[1];
333 : }
334 : else
335 : {
336 : const Real epsilon_sqrt = 1.0e-48;
337 :
338 7769 : const auto lamdba_plus_cg = contact_pressure + c * weighted_gap;
339 7769 : const auto normal_bound = max(0.0, lamdba_plus_cg);
340 : const auto friction_bound = mu_ad * normal_bound;
341 : const auto tangential_compliance =
342 7769 : _friction_elastic_slip > 0.0 ? _friction_elastic_slip / (friction_bound + _epsilon) : 0.0;
343 :
344 : std::array<ADReal, 2> lambda_t_plus_ctu;
345 : lambda_t_plus_ctu[0] =
346 7769 : friction_lm_values[0] +
347 15538 : c_t * (*tangential_vel[0] * _dt - tangential_compliance * friction_lm_values[0]);
348 : lambda_t_plus_ctu[1] =
349 7769 : friction_lm_values[1] +
350 15538 : c_t * (*tangential_vel[1] * _dt - tangential_compliance * friction_lm_values[1]);
351 :
352 : const auto tangential_trial_norm =
353 7769 : sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
354 7769 : lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt);
355 :
356 7769 : const auto term_1_x = max(friction_bound, tangential_trial_norm) * friction_lm_values[0];
357 :
358 15538 : const auto term_1_y = max(friction_bound, tangential_trial_norm) * friction_lm_values[1];
359 :
360 : const auto term_2_x = friction_bound * lambda_t_plus_ctu[0];
361 :
362 : const auto term_2_y = friction_bound * lambda_t_plus_ctu[1];
363 :
364 7769 : dof_residual = term_1_x - term_2_x;
365 7769 : dof_residual_dir = term_1_y - term_2_y;
366 : }
367 :
368 10668 : addResidualsAndJacobian(_assembly,
369 21336 : std::array<ADReal, 1>{{dof_residual}},
370 21336 : std::array<dof_id_type, 1>{{friction_dof_indices[0]}},
371 : _friction_vars[0]->scalingFactor());
372 10668 : addResidualsAndJacobian(_assembly,
373 21336 : std::array<ADReal, 1>{{dof_residual_dir}},
374 21336 : std::array<dof_id_type, 1>{{friction_dof_indices[1]}},
375 : _friction_vars[1]->scalingFactor());
376 10668 : }
377 :
378 : void
379 37266 : ComputeDynamicFrictionalForceLMMechanicalContact::enforceConstraintOnDof(
380 : const DofObject * const dof)
381 : {
382 : using std::max, std::abs, std::sqrt;
383 :
384 : // Get friction LM
385 37266 : const auto friction_dof_index = dof->dof_number(_sys.number(), _friction_vars[0]->number(), 0);
386 37266 : const ADReal & tangential_vel = *_tangential_vel_ptr[0];
387 37266 : ADReal friction_lm_value = (*_sys.currentSolution())(friction_dof_index);
388 : Moose::derivInsert(friction_lm_value.derivatives(), friction_dof_index, 1.);
389 :
390 : // Get normal LM
391 37266 : const auto normal_dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
392 37266 : const ADReal & weighted_gap = *_weighted_gap_ptr;
393 37266 : ADReal contact_pressure = (*_sys.currentSolution())(normal_dof_index);
394 : Moose::derivInsert(contact_pressure.derivatives(), normal_dof_index, 1.);
395 :
396 37266 : const Real contact_pressure_old = _sys.solutionOld()(normal_dof_index);
397 :
398 : // Get normalized c and c_t values (if normalization specified
399 37266 : const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
400 37266 : const Real c_t = _normalize_c ? _c_t / *_normalization_ptr : _c_t;
401 :
402 : // Compute the friction coefficient (constant or function)
403 : const auto & current_real_tangential_velocity =
404 37266 : libmesh_map_find(_dof_to_real_tangential_velocity, dof);
405 : const auto old_real_tangential_velocity_it = _dof_to_old_real_tangential_velocity.find(dof);
406 : const std::array<Real, 2> current_function_real_tangential_velocity{
407 : {MetaPhysicL::raw_value(current_real_tangential_velocity[0]),
408 37266 : MetaPhysicL::raw_value(current_real_tangential_velocity[1])}};
409 : const auto & function_real_tangential_velocity =
410 : old_real_tangential_velocity_it != _dof_to_old_real_tangential_velocity.end()
411 69190 : ? old_real_tangential_velocity_it->second
412 : : current_function_real_tangential_velocity;
413 : const ADReal slip_increment =
414 37266 : sqrt(current_real_tangential_velocity[0] * current_real_tangential_velocity[0] + 1.0e-24) *
415 37266 : _dt;
416 : ADReal mu_ad = computeFrictionValue(
417 37266 : contact_pressure_old, function_real_tangential_velocity[0], 0.0, slip_increment);
418 :
419 : ADReal dof_residual;
420 : // Primal-dual active set strategy (PDASS)
421 37266 : if (contact_pressure < _epsilon)
422 22933 : dof_residual = friction_lm_value;
423 : else
424 : {
425 14333 : const auto lambda_plus_cg = contact_pressure + c * weighted_gap;
426 14333 : const auto normal_bound = max(0.0, lambda_plus_cg);
427 : const auto friction_bound = mu_ad * normal_bound;
428 : const auto tangential_compliance =
429 18469 : _friction_elastic_slip > 0.0 ? _friction_elastic_slip / (friction_bound + _epsilon) : 0.0;
430 : const auto lambda_t_plus_ctu =
431 : friction_lm_value +
432 14333 : c_t * (tangential_vel * _dt - tangential_compliance * friction_lm_value);
433 :
434 28666 : const auto term_1 = max(friction_bound, abs(lambda_t_plus_ctu)) * friction_lm_value;
435 : const auto term_2 = friction_bound * lambda_t_plus_ctu;
436 :
437 14333 : dof_residual = term_1 - term_2;
438 : }
439 :
440 37266 : addResidualsAndJacobian(_assembly,
441 74532 : std::array<ADReal, 1>{{dof_residual}},
442 74532 : std::array<dof_id_type, 1>{{friction_dof_index}},
443 : _friction_vars[0]->scalingFactor());
444 37266 : }
445 :
446 : ADReal
447 47934 : ComputeDynamicFrictionalForceLMMechanicalContact::computeFrictionValue(
448 : const ADReal & contact_pressure,
449 : const Real & function_tangential_vel,
450 : const Real & function_tangential_vel_dir,
451 : const ADReal & slip_increment)
452 : {
453 : using std::sqrt;
454 :
455 : // TODO: Introduce temperature dependence in the function. Do this when we have an example.
456 : ADReal mu_ad;
457 :
458 47934 : if (!_has_friction_function)
459 46710 : mu_ad = _mu;
460 : else
461 : {
462 : ADReal tangential_vel_magnitude =
463 2448 : sqrt(function_tangential_vel * function_tangential_vel +
464 1224 : function_tangential_vel_dir * function_tangential_vel_dir + 1.0e-24);
465 :
466 1224 : mu_ad = _function_friction->value<ADReal>(0.0, contact_pressure, tangential_vel_magnitude, 0.0);
467 : }
468 :
469 : return Moose::Contact::regularizedFrictionCoefficient(
470 47934 : mu_ad, slip_increment, _friction_coefficient_regularization, _friction_reference_slip);
471 : }
|
