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 164 : PenaltyFrictionUserObject::validParams()
25 : {
26 164 : InputParameters params = WeightedVelocitiesUserObject::validParams();
27 164 : params += PenaltyWeightedGapUserObject::validParams();
28 :
29 164 : params.addClassDescription(
30 : "Computes the mortar frictional contact force via a penalty approach.");
31 328 : 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 328 : params.addRequiredParam<Real>("friction_coefficient",
35 : "The friction coefficient ruling Coulomb friction equations.");
36 328 : params.addRangeCheckedParam<Real>(
37 : "slip_tolerance",
38 : "slip_tolerance > 0",
39 : "Acceptable slip distance at which augmented Lagrange iterations can be stopped");
40 328 : MooseEnum adaptivity_penalty_friction("SIMPLE FRICTION_LIMIT", "FRICTION_LIMIT");
41 328 : adaptivity_penalty_friction.addDocumentation(
42 : "SIMPLE", "Keep multiplying by the frictional penalty multiplier between AL iterations");
43 328 : 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 328 : params.addParam<MooseEnum>(
48 : "adaptivity_penalty_friction",
49 : adaptivity_penalty_friction,
50 : "The augmented Lagrange update strategy used on the frictional penalty coefficient.");
51 492 : params.addRangeCheckedParam<Real>(
52 : "penalty_multiplier_friction",
53 328 : 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 164 : return params;
58 164 : }
59 :
60 82 : 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 82 : _penalty(getParam<Real>("penalty")),
78 328 : _penalty_friction(isParamValid("penalty_friction") ? getParam<Real>("penalty_friction")
79 82 : : getParam<Real>("penalty")),
80 198 : _slip_tolerance(isParamValid("slip_tolerance") ? getParam<Real>("slip_tolerance") : 0.0),
81 164 : _friction_coefficient(getParam<Real>("friction_coefficient")),
82 164 : _penalty_multiplier_friction(getParam<Real>("penalty_multiplier_friction")),
83 82 : _adaptivity_friction(
84 82 : getParam<MooseEnum>("adaptivity_penalty_friction").getEnum<AdaptivityFrictionalPenalty>()),
85 82 : _epsilon_tolerance(1.0e-40)
86 :
87 : {
88 164 : 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 82 : }
93 :
94 : const VariableTestValue &
95 82 : PenaltyFrictionUserObject::test() const
96 : {
97 82 : return _aux_lm_var ? _aux_lm_var->phiLower() : _disp_x_var->phiLower();
98 : }
99 :
100 : const ADVariableValue &
101 5049904 : PenaltyFrictionUserObject::contactTangentialPressureDirOne() const
102 : {
103 5049904 : 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 361 : PenaltyFrictionUserObject::timestepSetup()
114 : {
115 : // these functions do not call WeightedGapUserObject::timestepSetup to avoid double initialization
116 361 : WeightedVelocitiesUserObject::selfTimestepSetup();
117 361 : PenaltyWeightedGapUserObject::selfTimestepSetup();
118 :
119 : // instead we call it explicitly here
120 361 : WeightedGapUserObject::timestepSetup();
121 :
122 : // Clear step slip (values used in between AL iterations for penalty adaptivity)
123 452 : 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 1984 : 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 366 : for (auto & dof_lp : _dof_to_local_penalty_friction)
138 5 : dof_lp.second = _penalty_friction;
139 :
140 : // save off tangential traction from the last timestep
141 1984 : 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 3246 : tangential_traction = {0.0, 0.0};
147 : }
148 :
149 366 : for (auto & [dof_object, delta_tangential_lm] : _dof_to_frictional_lagrange_multipliers)
150 : delta_tangential_lm.setZero();
151 361 : }
152 :
153 : void
154 9016 : PenaltyFrictionUserObject::initialize()
155 : {
156 : // these functions do not call WeightedGapUserObject::initialize to avoid double initialization
157 9016 : WeightedVelocitiesUserObject::selfInitialize();
158 9016 : PenaltyWeightedGapUserObject::selfInitialize();
159 :
160 : // instead we call it explicitly here
161 9016 : WeightedGapUserObject::initialize();
162 9016 : }
163 :
164 : Real
165 155970 : PenaltyFrictionUserObject::getFrictionalContactPressure(const Node * const node,
166 : const unsigned int component) const
167 : {
168 155970 : const auto it = _dof_to_tangential_traction.find(_subproblem.mesh().nodePtr(node->id()));
169 :
170 155970 : if (it != _dof_to_tangential_traction.end())
171 4432 : return MetaPhysicL::raw_value(it->second.first(component));
172 : else
173 : return 0.0;
174 : }
175 :
176 : Real
177 155970 : PenaltyFrictionUserObject::getAccumulatedSlip(const Node * const node,
178 : const unsigned int component) const
179 : {
180 155970 : const auto it = _dof_to_accumulated_slip.find(_subproblem.mesh().nodePtr(node->id()));
181 :
182 155970 : if (it != _dof_to_accumulated_slip.end())
183 4432 : return MetaPhysicL::raw_value(it->second.first(component));
184 : else
185 : return 0.0;
186 : }
187 :
188 : Real
189 150330 : PenaltyFrictionUserObject::getTangentialVelocity(const Node * const node,
190 : const unsigned int component) const
191 : {
192 150330 : const auto it = _dof_to_real_tangential_velocity.find(_subproblem.mesh().nodePtr(node->id()));
193 :
194 150330 : if (it != _dof_to_real_tangential_velocity.end())
195 4088 : 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 136419 : PenaltyFrictionUserObject::reinit()
215 : {
216 : // Normal contact pressure with penalty
217 136419 : PenaltyWeightedGapUserObject::reinit();
218 :
219 : // Reset frictional pressure
220 136419 : _frictional_contact_traction_one.resize(_qrule_msm->n_points());
221 136419 : _frictional_contact_traction_two.resize(_qrule_msm->n_points()); // 3D
222 480937 : for (const auto qp : make_range(_qrule_msm->n_points()))
223 : {
224 344518 : _frictional_contact_traction_one[qp] = 0.0;
225 344518 : _frictional_contact_traction_two[qp] = 0.0;
226 : }
227 :
228 : // zero vector
229 136419 : const static TwoVector zero{0.0, 0.0};
230 :
231 : // iterate over nodes
232 480937 : for (const auto i : make_range(_test->size()))
233 : {
234 : // current node
235 344518 : const Node * const node = _lower_secondary_elem->node_ptr(i);
236 :
237 689036 : const auto penalty_friction = findValue(
238 344518 : _dof_to_local_penalty_friction, static_cast<const DofObject *>(node), _penalty_friction);
239 :
240 : // utilized quantities
241 344518 : const auto & normal_pressure = _dof_to_normal_pressure[node];
242 :
243 : // map the tangential traction and accumulated slip
244 344518 : auto & [tangential_traction, old_tangential_traction] = _dof_to_tangential_traction[node];
245 344518 : auto & [accumulated_slip, old_accumulated_slip] = _dof_to_accumulated_slip[node];
246 :
247 : Real normal_lm = -1;
248 344518 : if (_dof_to_lagrange_multiplier.find(node) != _dof_to_lagrange_multiplier.end())
249 117180 : normal_lm = libmesh_map_find(_dof_to_lagrange_multiplier, node);
250 :
251 : // Keep active set fixed from second Uzawa loop
252 344518 : if (normal_lm < -TOLERANCE && normal_pressure > TOLERANCE)
253 : {
254 : using namespace std;
255 :
256 : const auto & real_tangential_velocity =
257 60748 : libmesh_map_find(_dof_to_real_tangential_velocity, node);
258 60748 : const ADTwoVector slip_distance = {real_tangential_velocity[0] * _dt,
259 60748 : real_tangential_velocity[1] * _dt};
260 :
261 : // frictional lagrange multiplier (delta lambda^(k)_T)
262 : const auto & tangential_lm =
263 60748 : _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 60748 : ADTwoVector inner_iteration_penalty_friction = penalty_friction * slip_distance;
269 :
270 60748 : const auto slip_metric = MetaPhysicL::raw_value(slip_distance).cwiseAbs().norm();
271 :
272 121490 : if (slip_metric > _epsilon_tolerance &&
273 182232 : penalty_friction * slip_distance.norm() >
274 242974 : 0.4 * _friction_coefficient * std::abs(normal_pressure))
275 : {
276 : inner_iteration_penalty_friction =
277 55436 : MetaPhysicL::raw_value(0.4 * _friction_coefficient * std::abs(normal_pressure) /
278 110872 : (penalty_friction * slip_distance.norm())) *
279 55436 : penalty_friction * slip_distance;
280 : }
281 :
282 : ADTwoVector tangential_trial_traction =
283 60748 : old_tangential_traction + tangential_lm + inner_iteration_penalty_friction;
284 :
285 : // Nonlinearity below
286 60748 : ADReal tangential_trial_traction_norm = tangential_trial_traction.norm();
287 121496 : 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 60748 : if (phi_trial > 0.0)
292 : // Simo 3.13/3.14 (the penalty formulation has an error in the paper)
293 : tangential_traction -=
294 9544 : phi_trial * tangential_trial_traction / tangential_trial_traction_norm;
295 :
296 : // track accumulated slip for output purposes
297 60748 : accumulated_slip = old_accumulated_slip + MetaPhysicL::raw_value(slip_distance).cwiseAbs();
298 :
299 : // Keep track of slip vector for adaptive penalty
300 60748 : 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 283770 : tangential_traction.setZero();
308 : }
309 :
310 : // Now that we have consistent nodal frictional values, create an interpolated frictional
311 : // pressure variable.
312 344518 : const auto & test_i = (*_test)[i];
313 1320274 : for (const auto qp : make_range(_qrule_msm->n_points()))
314 : {
315 1951512 : _frictional_contact_traction_one[qp] += test_i[qp] * tangential_traction(0);
316 975756 : _frictional_contact_traction_two[qp] += test_i[qp] * tangential_traction(1);
317 : }
318 : }
319 136419 : }
320 :
321 : void
322 9016 : PenaltyFrictionUserObject::finalize()
323 : {
324 9016 : WeightedVelocitiesUserObject::finalize();
325 9016 : PenaltyWeightedGapUserObject::selfFinalize();
326 9016 : }
327 :
328 : bool
329 96 : PenaltyFrictionUserObject::isAugmentedLagrangianConverged()
330 : {
331 : // save off step slip
332 : // This method is called at the beginning of the AL iteration.
333 383 : 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 96 : std::pair<Real, dof_id_type> max_slip{0.0, 0};
340 :
341 4602 : for (const auto & [dof_object, traction_pair] : _dof_to_tangential_traction)
342 : {
343 4506 : const auto & tangential_traction = traction_pair.first;
344 4506 : 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 4506 : if (_dof_to_real_tangential_velocity.find(dof_object) != _dof_to_real_tangential_velocity.end())
354 : {
355 : const auto & real_tangential_velocity =
356 4506 : 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 13518 : if (_friction_coefficient * normal_pressure < tangential_traction.norm() * (1 + TOLERANCE))
363 : {
364 : // If it's slipping, any slip distance is physical.
365 : }
366 4446 : else if (slip_velocity.norm() * _dt > _slip_tolerance && normal_pressure > TOLERANCE)
367 : {
368 : const auto new_slip =
369 354 : 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 96 : this->_communicator.max(max_slip);
377 :
378 : // Check normal contact convergence now, to make sure all ranks get here
379 96 : if (!PenaltyWeightedGapUserObject::isAugmentedLagrangianConverged())
380 : return false;
381 :
382 : // Did we observe any above tolerance slip anywhere?
383 85 : if (max_slip.first > _slip_tolerance)
384 : {
385 67 : 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 40 : _slip_tolerance);
393 : }
394 67 : return false;
395 : }
396 :
397 : return true;
398 : }
399 :
400 : void
401 78 : PenaltyFrictionUserObject::updateAugmentedLagrangianMultipliers()
402 : {
403 78 : PenaltyWeightedGapUserObject::updateAugmentedLagrangianMultipliers();
404 :
405 299 : for (auto & [dof_object, tangential_lm] : _dof_to_frictional_lagrange_multipliers)
406 : {
407 : auto & penalty_friction = _dof_to_local_penalty_friction[dof_object];
408 221 : if (penalty_friction == 0.0)
409 61 : penalty_friction = _penalty_friction;
410 :
411 : // normal quantities
412 221 : const auto & normal_lm = libmesh_map_find(_dof_to_lagrange_multiplier, dof_object);
413 :
414 : // tangential quantities
415 : const auto & real_tangential_velocity =
416 221 : 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 221 : 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 221 : if (tangential_trial_traction_norm * (1 + TOLERANCE) <=
428 221 : std::abs(_friction_coefficient * normal_lm))
429 : {
430 175 : tangential_lm += penalty_friction * slip_velocity * _dt;
431 : }
432 : else
433 : {
434 : tangential_lm = -tangential_trial_traction / tangential_trial_traction_norm *
435 46 : penalty_friction * normal_lm -
436 46 : old_tangential_traction;
437 : }
438 :
439 : // Update penalty.
440 221 : if (_adaptivity_friction == AdaptivityFrictionalPenalty::SIMPLE)
441 : {
442 112 : if (_slip_tolerance < _dt * slip_velocity.norm())
443 42 : penalty_friction *= _penalty_multiplier_friction;
444 :
445 : // Provide the user the ability of setting this maximum penalty
446 56 : if (penalty_friction > _penalty_friction * _max_penalty_multiplier)
447 0 : penalty_friction = _penalty_friction * _max_penalty_multiplier;
448 : }
449 165 : else if (_adaptivity_friction == AdaptivityFrictionalPenalty::FRICTION_LIMIT)
450 : {
451 165 : 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 165 : if (step_slip.first.dot(step_slip.second) > 0.0 && std::abs(normal_lm) > TOLERANCE &&
454 330 : _slip_tolerance < _dt * slip_velocity.norm())
455 : {
456 125 : penalty_friction =
457 125 : (_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 40 : 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 165 : if (penalty_friction < _penalty_friction)
466 12 : penalty_friction = _penalty_friction;
467 153 : else if (penalty_friction > _penalty_friction * _max_penalty_multiplier)
468 0 : penalty_friction = _penalty_friction * _max_penalty_multiplier;
469 : }
470 : }
471 78 : }
|
