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 "PenaltyFrictionUserObject.h"
11 : #include "MooseVariableFE.h"
12 : #include "SystemBase.h"
13 : #include "MortarUtils.h"
14 : #include "MooseUtils.h"
15 : #include "MathUtils.h"
16 : #include "MortarContactUtils.h"
17 : #include "ADReal.h"
18 :
19 : #include <Eigen/Core>
20 :
21 : registerMooseObject("ContactApp", PenaltyFrictionUserObject);
22 :
23 : InputParameters
24 198 : PenaltyFrictionUserObject::validParams()
25 : {
26 198 : InputParameters params = WeightedVelocitiesUserObject::validParams();
27 198 : params += PenaltyWeightedGapUserObject::validParams();
28 :
29 198 : params.addClassDescription(
30 : "Computes the mortar frictional contact force via a penalty approach.");
31 396 : params.addParam<Real>("penalty_friction",
32 : "The penalty factor for frictional interaction. If not provide, the normal "
33 : "penalty factor is also used for the frictional problem.");
34 396 : params.addRequiredParam<Real>("friction_coefficient",
35 : "The friction coefficient ruling Coulomb friction equations.");
36 396 : params.addRangeCheckedParam<Real>(
37 : "slip_tolerance",
38 : "slip_tolerance > 0",
39 : "Acceptable slip distance at which augmented Lagrange iterations can be stopped");
40 396 : MooseEnum adaptivity_penalty_friction("SIMPLE FRICTION_LIMIT", "FRICTION_LIMIT");
41 396 : adaptivity_penalty_friction.addDocumentation(
42 : "SIMPLE", "Keep multiplying by the frictional penalty multiplier between AL iterations");
43 396 : adaptivity_penalty_friction.addDocumentation(
44 : "FRICTION_LIMIT",
45 : "This strategy will be guided by the Coulomb limit and be less reliant on the initial "
46 : "penalty factor provided by the user.");
47 396 : params.addParam<MooseEnum>(
48 : "adaptivity_penalty_friction",
49 : adaptivity_penalty_friction,
50 : "The augmented Lagrange update strategy used on the frictional penalty coefficient.");
51 594 : params.addRangeCheckedParam<Real>(
52 : "penalty_multiplier_friction",
53 396 : 1.0,
54 : "penalty_multiplier_friction > 0",
55 : "The penalty growth factor between augmented Lagrange "
56 : "iterations for penalizing relative slip distance if the node is under stick conditions.");
57 198 : return params;
58 198 : }
59 :
60 99 : PenaltyFrictionUserObject::PenaltyFrictionUserObject(const InputParameters & parameters)
61 : /*
62 : * We are using virtual inheritance to avoid the "Diamond inheritance" problem. This means that
63 : * that we have to construct WeightedGapUserObject explicitly as it will_not_ be constructed in
64 : * the intermediate base classes PenaltyWeightedGapUserObject and WeightedVelocitiesUserObject.
65 : * Virtual inheritance ensures that only one instance of WeightedGapUserObject is included in this
66 : * class. The inheritance diagram is as follows:
67 : *
68 : * WeightedGapUserObject <----- PenaltyWeightedGapUserObject
69 : * ^ ^
70 : * | |
71 : * WeightedVelocitiesUserObject <----- PenaltyFrictionUserObject
72 : *
73 : */
74 : : WeightedGapUserObject(parameters),
75 : PenaltyWeightedGapUserObject(parameters),
76 : WeightedVelocitiesUserObject(parameters),
77 99 : _penalty(getParam<Real>("penalty")),
78 396 : _penalty_friction(isParamValid("penalty_friction") ? getParam<Real>("penalty_friction")
79 99 : : getParam<Real>("penalty")),
80 270 : _slip_tolerance(isParamValid("slip_tolerance") ? getParam<Real>("slip_tolerance") : 0.0),
81 198 : _friction_coefficient(getParam<Real>("friction_coefficient")),
82 198 : _penalty_multiplier_friction(getParam<Real>("penalty_multiplier_friction")),
83 99 : _adaptivity_friction(
84 99 : getParam<MooseEnum>("adaptivity_penalty_friction").getEnum<AdaptivityFrictionalPenalty>()),
85 99 : _epsilon_tolerance(1.0e-40)
86 :
87 : {
88 198 : if (!_augmented_lagrange_problem == isParamValid("slip_tolerance"))
89 0 : paramError("slip_tolerance",
90 : "This parameter must be supplied if and only if an augmented Lagrange problem "
91 : "object is used.");
92 99 : }
93 :
94 : const VariableTestValue &
95 99 : PenaltyFrictionUserObject::test() const
96 : {
97 99 : return _aux_lm_var ? _aux_lm_var->phiLower() : _disp_x_var->phiLower();
98 : }
99 :
100 : const ADVariableValue &
101 5318272 : PenaltyFrictionUserObject::contactTangentialPressureDirOne() const
102 : {
103 5318272 : return _frictional_contact_traction_one;
104 : }
105 :
106 : const ADVariableValue &
107 3440640 : PenaltyFrictionUserObject::contactTangentialPressureDirTwo() const
108 : {
109 3440640 : return _frictional_contact_traction_two;
110 : }
111 :
112 : void
113 343 : PenaltyFrictionUserObject::timestepSetup()
114 : {
115 : // these functions do not call WeightedGapUserObject::timestepSetup to avoid double initialization
116 343 : WeightedVelocitiesUserObject::selfTimestepSetup();
117 343 : PenaltyWeightedGapUserObject::selfTimestepSetup();
118 :
119 : // instead we call it explicitly here
120 343 : WeightedGapUserObject::timestepSetup();
121 :
122 : // Clear step slip (values used in between AL iterations for penalty adaptivity)
123 413 : for (auto & map_pr : _dof_to_step_slip)
124 : {
125 : auto & [step_slip, old_step_slip] = map_pr.second;
126 : old_step_slip = {0.0, 0.0};
127 : step_slip = {0.0, 0.0};
128 : }
129 :
130 : // save off accumulated slip from the last timestep
131 1629 : for (auto & map_pr : _dof_to_accumulated_slip)
132 : {
133 : auto & [accumulated_slip, old_accumulated_slip] = map_pr.second;
134 : old_accumulated_slip = accumulated_slip;
135 : }
136 :
137 343 : for (auto & dof_lp : _dof_to_local_penalty_friction)
138 0 : dof_lp.second = _penalty_friction;
139 :
140 : // save off tangential traction from the last timestep
141 1629 : for (auto & map_pr : _dof_to_tangential_traction)
142 : {
143 : auto & [tangential_traction, old_tangential_traction] = map_pr.second;
144 : old_tangential_traction = {MetaPhysicL::raw_value(tangential_traction(0)),
145 : MetaPhysicL::raw_value(tangential_traction(1))};
146 2572 : tangential_traction = {0.0, 0.0};
147 : }
148 :
149 343 : for (auto & [dof_object, delta_tangential_lm] : _dof_to_frictional_lagrange_multipliers)
150 : delta_tangential_lm.setZero();
151 343 : }
152 :
153 : void
154 9322 : PenaltyFrictionUserObject::initialize()
155 : {
156 : // these functions do not call WeightedGapUserObject::initialize to avoid double initialization
157 9322 : WeightedVelocitiesUserObject::selfInitialize();
158 9322 : PenaltyWeightedGapUserObject::selfInitialize();
159 :
160 : // instead we call it explicitly here
161 9322 : WeightedGapUserObject::initialize();
162 9322 : }
163 :
164 : Real
165 129696 : PenaltyFrictionUserObject::getFrictionalContactPressure(const Node * const node,
166 : const unsigned int component) const
167 : {
168 129696 : const auto it = _dof_to_tangential_traction.find(_subproblem.mesh().nodePtr(node->id()));
169 :
170 129696 : if (it != _dof_to_tangential_traction.end())
171 4136 : return MetaPhysicL::raw_value(it->second.first(component));
172 : else
173 : return 0.0;
174 : }
175 :
176 : Real
177 129696 : PenaltyFrictionUserObject::getAccumulatedSlip(const Node * const node,
178 : const unsigned int component) const
179 : {
180 129696 : const auto it = _dof_to_accumulated_slip.find(_subproblem.mesh().nodePtr(node->id()));
181 :
182 129696 : if (it != _dof_to_accumulated_slip.end())
183 4136 : return MetaPhysicL::raw_value(it->second.first(component));
184 : else
185 : return 0.0;
186 : }
187 :
188 : Real
189 124056 : PenaltyFrictionUserObject::getTangentialVelocity(const Node * const node,
190 : const unsigned int component) const
191 : {
192 124056 : const auto it = _dof_to_real_tangential_velocity.find(_subproblem.mesh().nodePtr(node->id()));
193 :
194 124056 : if (it != _dof_to_real_tangential_velocity.end())
195 3796 : return MetaPhysicL::raw_value(it->second[component]);
196 : else
197 : return 0.0;
198 : }
199 :
200 : Real
201 0 : PenaltyFrictionUserObject::getDeltaTangentialLagrangeMultiplier(const Node * const node,
202 : const unsigned int component) const
203 : {
204 : const auto it =
205 0 : _dof_to_frictional_lagrange_multipliers.find(_subproblem.mesh().nodePtr(node->id()));
206 :
207 0 : if (it != _dof_to_frictional_lagrange_multipliers.end())
208 0 : return MetaPhysicL::raw_value(it->second[component]);
209 : else
210 : return 0.0;
211 : }
212 :
213 : void
214 153192 : PenaltyFrictionUserObject::reinit()
215 : {
216 : // Normal contact pressure with penalty
217 153192 : PenaltyWeightedGapUserObject::reinit();
218 :
219 : // Reset frictional pressure
220 153192 : _frictional_contact_traction_one.resize(_qrule_msm->n_points());
221 153192 : _frictional_contact_traction_two.resize(_qrule_msm->n_points()); // 3D
222 531256 : for (const auto qp : make_range(_qrule_msm->n_points()))
223 : {
224 378064 : _frictional_contact_traction_one[qp] = 0.0;
225 378064 : _frictional_contact_traction_two[qp] = 0.0;
226 : }
227 :
228 : // zero vector
229 153192 : const static TwoVector zero{0.0, 0.0};
230 :
231 : // iterate over nodes
232 531256 : for (const auto i : make_range(_test->size()))
233 : {
234 : // current node
235 378064 : const Node * const node = _lower_secondary_elem->node_ptr(i);
236 :
237 : const auto penalty_friction = findValue(
238 378064 : _dof_to_local_penalty_friction, static_cast<const DofObject *>(node), _penalty_friction);
239 :
240 : // utilized quantities
241 378064 : const auto & normal_pressure = _dof_to_normal_pressure[node];
242 :
243 : // map the tangential traction and accumulated slip
244 378064 : auto & [tangential_traction, old_tangential_traction] = _dof_to_tangential_traction[node];
245 378064 : auto & [accumulated_slip, old_accumulated_slip] = _dof_to_accumulated_slip[node];
246 :
247 : Real normal_lm = -1;
248 378064 : if (_dof_to_lagrange_multiplier.find(node) != _dof_to_lagrange_multiplier.end())
249 165568 : normal_lm = libmesh_map_find(_dof_to_lagrange_multiplier, node);
250 :
251 : // Keep active set fixed from second Uzawa loop
252 378064 : if (normal_lm < -TOLERANCE && normal_pressure > TOLERANCE)
253 : {
254 : using namespace std;
255 :
256 : const auto & real_tangential_velocity =
257 62048 : libmesh_map_find(_dof_to_real_tangential_velocity, node);
258 62048 : const ADTwoVector slip_distance = {real_tangential_velocity[0] * _dt,
259 62048 : real_tangential_velocity[1] * _dt};
260 :
261 : // frictional lagrange multiplier (delta lambda^(k)_T)
262 : const auto & tangential_lm =
263 62048 : _augmented_lagrange_problem ? _dof_to_frictional_lagrange_multipliers[node] : zero;
264 :
265 : // tangential trial traction (Simo 3.12)
266 : // Modified for implementation in MOOSE: Avoid pingponging on frictional sign (max. 0.4
267 : // capacity)
268 62048 : ADTwoVector inner_iteration_penalty_friction = penalty_friction * slip_distance;
269 :
270 62048 : const auto slip_metric = MetaPhysicL::raw_value(slip_distance).cwiseAbs().norm();
271 :
272 124096 : if (slip_metric > _epsilon_tolerance &&
273 186144 : penalty_friction * slip_distance.norm() >
274 248192 : 0.4 * _friction_coefficient * std::abs(normal_pressure))
275 : {
276 : inner_iteration_penalty_friction =
277 55232 : MetaPhysicL::raw_value(0.4 * _friction_coefficient * std::abs(normal_pressure) /
278 110464 : (penalty_friction * slip_distance.norm())) *
279 55232 : penalty_friction * slip_distance;
280 : }
281 :
282 : ADTwoVector tangential_trial_traction =
283 62048 : old_tangential_traction + tangential_lm + inner_iteration_penalty_friction;
284 :
285 : // Nonlinearity below
286 62048 : ADReal tangential_trial_traction_norm = tangential_trial_traction.norm();
287 124096 : ADReal phi_trial = tangential_trial_traction_norm - _friction_coefficient * normal_pressure;
288 : tangential_traction = tangential_trial_traction;
289 :
290 : // Simo considers this a 'return mapping'; we are just capping friction to the Coulomb limit.
291 62048 : if (phi_trial > 0.0)
292 : // Simo 3.13/3.14 (the penalty formulation has an error in the paper)
293 : tangential_traction -=
294 9104 : phi_trial * tangential_trial_traction / tangential_trial_traction_norm;
295 :
296 : // track accumulated slip for output purposes
297 62048 : accumulated_slip = old_accumulated_slip + MetaPhysicL::raw_value(slip_distance).cwiseAbs();
298 :
299 : // Keep track of slip vector for adaptive penalty
300 62048 : auto & [step_slip, old_step_slip] = _dof_to_step_slip[node];
301 : step_slip = MetaPhysicL::raw_value(slip_distance);
302 : }
303 : else
304 : {
305 : // reset slip and clear traction
306 : accumulated_slip.setZero();
307 316016 : tangential_traction.setZero();
308 : }
309 :
310 : // Now that we have consistent nodal frictional values, create an interpolated frictional
311 : // pressure variable.
312 378064 : const auto & test_i = (*_test)[i];
313 1420912 : for (const auto qp : make_range(_qrule_msm->n_points()))
314 : {
315 2085696 : _frictional_contact_traction_one[qp] += test_i[qp] * tangential_traction(0);
316 1042848 : _frictional_contact_traction_two[qp] += test_i[qp] * tangential_traction(1);
317 : }
318 : }
319 153192 : }
320 :
321 : void
322 9322 : PenaltyFrictionUserObject::finalize()
323 : {
324 9322 : WeightedVelocitiesUserObject::finalize();
325 9322 : PenaltyWeightedGapUserObject::selfFinalize();
326 9322 : }
327 :
328 : bool
329 143 : PenaltyFrictionUserObject::isAugmentedLagrangianConverged()
330 : {
331 : // save off step slip
332 : // This method is called at the beginning of the AL iteration.
333 548 : for (auto & map_pr : _dof_to_step_slip)
334 : {
335 : auto & [step_slip, old_step_slip] = map_pr.second;
336 : old_step_slip = step_slip;
337 : }
338 :
339 143 : std::pair<Real, dof_id_type> max_slip{0.0, 0};
340 :
341 6446 : for (const auto & [dof_object, traction_pair] : _dof_to_tangential_traction)
342 : {
343 6303 : const auto & tangential_traction = traction_pair.first;
344 6303 : auto normal_pressure = _dof_to_normal_pressure[dof_object];
345 :
346 : // We may not find a node in the map because AL machinery gets called at different system
347 : // configurations. That is, we may want to find a key from a node that computed a zero traction
348 : // on the verge of not projecting. Since, when we computed the velocity, the system was at a
349 : // slightly different configuration, that node may not a computed physical weighted velocity.
350 : // This doesn't seem an issue at all as non-projecting nodes should be a corner case not
351 : // affecting the physics.
352 : TwoVector slip_velocity = {0.0, 0.0};
353 6303 : if (_dof_to_real_tangential_velocity.find(dof_object) != _dof_to_real_tangential_velocity.end())
354 : {
355 : const auto & real_tangential_velocity =
356 6303 : libmesh_map_find(_dof_to_real_tangential_velocity, dof_object);
357 : slip_velocity = {MetaPhysicL::raw_value(real_tangential_velocity[0]),
358 : MetaPhysicL::raw_value(real_tangential_velocity[1])};
359 : }
360 :
361 : // Check slip/stick
362 18909 : if (_friction_coefficient * normal_pressure < tangential_traction.norm() * (1 + TOLERANCE))
363 : {
364 : // If it's slipping, any slip distance is physical.
365 : }
366 6223 : else if (slip_velocity.norm() * _dt > _slip_tolerance && normal_pressure > TOLERANCE)
367 : {
368 : const auto new_slip =
369 516 : std::make_pair<Real, dof_id_type>(slip_velocity.norm() * _dt, dof_object->id());
370 : if (new_slip > max_slip)
371 : max_slip = new_slip;
372 : }
373 : }
374 :
375 : // Communicate max_slip here where all ranks get to
376 143 : this->_communicator.max(max_slip);
377 :
378 : // Check normal contact convergence now, to make sure all ranks get here
379 143 : if (!PenaltyWeightedGapUserObject::isAugmentedLagrangianConverged())
380 : return false;
381 :
382 : // Did we observe any above tolerance slip anywhere?
383 127 : if (max_slip.first > _slip_tolerance)
384 : {
385 101 : if (this->_communicator.rank() == 0)
386 : {
387 : mooseInfoRepeated("Stick tolerance fails. Slip distance for sticking node ",
388 : max_slip.second,
389 : " is: ",
390 : max_slip.first,
391 : ", but slip tolerance is chosen to be ",
392 56 : _slip_tolerance);
393 : }
394 101 : return false;
395 : }
396 :
397 : return true;
398 : }
399 :
400 : void
401 117 : PenaltyFrictionUserObject::updateAugmentedLagrangianMultipliers()
402 : {
403 117 : PenaltyWeightedGapUserObject::updateAugmentedLagrangianMultipliers();
404 :
405 432 : for (auto & [dof_object, tangential_lm] : _dof_to_frictional_lagrange_multipliers)
406 : {
407 : auto & penalty_friction = _dof_to_local_penalty_friction[dof_object];
408 315 : if (penalty_friction == 0.0)
409 90 : penalty_friction = _penalty_friction;
410 :
411 : // normal quantities
412 315 : const auto & normal_lm = libmesh_map_find(_dof_to_lagrange_multiplier, dof_object);
413 :
414 : // tangential quantities
415 : const auto & real_tangential_velocity =
416 315 : libmesh_map_find(_dof_to_real_tangential_velocity, dof_object);
417 : const TwoVector slip_velocity = {MetaPhysicL::raw_value(real_tangential_velocity[0]),
418 : MetaPhysicL::raw_value(real_tangential_velocity[1])};
419 :
420 : auto & [tangential_traction, old_tangential_traction] = _dof_to_tangential_traction[dof_object];
421 :
422 : const TwoVector tangential_trial_traction =
423 315 : old_tangential_traction + tangential_lm + penalty_friction * slip_velocity * _dt;
424 : const Real tangential_trial_traction_norm = tangential_trial_traction.norm();
425 :
426 : // Augment
427 315 : if (tangential_trial_traction_norm * (1 + TOLERANCE) <=
428 315 : std::abs(_friction_coefficient * normal_lm))
429 : {
430 251 : tangential_lm += penalty_friction * slip_velocity * _dt;
431 : }
432 : else
433 : {
434 : tangential_lm = -tangential_trial_traction / tangential_trial_traction_norm *
435 64 : penalty_friction * normal_lm -
436 64 : old_tangential_traction;
437 : }
438 :
439 : // Update penalty.
440 315 : if (_adaptivity_friction == AdaptivityFrictionalPenalty::SIMPLE)
441 : {
442 180 : if (_slip_tolerance < _dt * slip_velocity.norm())
443 68 : penalty_friction *= _penalty_multiplier_friction;
444 :
445 : // Provide the user the ability of setting this maximum penalty
446 90 : if (penalty_friction > _penalty_friction * _max_penalty_multiplier)
447 0 : penalty_friction = _penalty_friction * _max_penalty_multiplier;
448 : }
449 225 : else if (_adaptivity_friction == AdaptivityFrictionalPenalty::FRICTION_LIMIT)
450 : {
451 225 : const auto & step_slip = libmesh_map_find(_dof_to_step_slip, dof_object);
452 : // No change of direction: Adjust penalty factor for the frictional problem
453 225 : if (step_slip.first.dot(step_slip.second) > 0.0 && std::abs(normal_lm) > TOLERANCE &&
454 450 : _slip_tolerance < _dt * slip_velocity.norm())
455 : {
456 170 : penalty_friction =
457 170 : (_friction_coefficient * std::abs(normal_lm)) / 2 / (_dt * slip_velocity.norm());
458 : // Alternative: accumulated_slip.norm() - old_accumulated_slip.norm()
459 : }
460 : // Change of direction: Reduce penalty factor to avoid lack of convergence
461 55 : else if (step_slip.first.dot(step_slip.second) < 0.0)
462 0 : penalty_friction /= 2.0;
463 :
464 : // Heuristics to bound the penalty factor
465 225 : if (penalty_friction < _penalty_friction)
466 16 : penalty_friction = _penalty_friction;
467 209 : else if (penalty_friction > _penalty_friction * _max_penalty_multiplier)
468 0 : penalty_friction = _penalty_friction * _max_penalty_multiplier;
469 : }
470 : }
471 117 : }
|
