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 "ComputeFrictionalForceLMMechanicalContact.h"
11 : #include "DisplacedProblem.h"
12 : #include "Assembly.h"
13 : #include "MortarContactUtils.h"
14 : #include "ContactFrictionUtils.h"
15 : #include "WeightedVelocitiesUserObject.h"
16 :
17 : #include "metaphysicl/metaphysicl_version.h"
18 : #include "metaphysicl/dualsemidynamicsparsenumberarray.h"
19 : #include "metaphysicl/parallel_dualnumber.h"
20 : #if METAPHYSICL_MAJOR_VERSION < 2
21 : #include "metaphysicl/parallel_dynamic_std_array_wrapper.h"
22 : #else
23 : #include "metaphysicl/parallel_dynamic_array_wrapper.h"
24 : #endif
25 : #include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
26 : #include "timpi/parallel_sync.h"
27 :
28 : #include <cmath>
29 :
30 : registerMooseObject("ContactApp", ComputeFrictionalForceLMMechanicalContact);
31 :
32 : InputParameters
33 354 : ComputeFrictionalForceLMMechanicalContact::validParams()
34 : {
35 354 : InputParameters params = ComputeWeightedGapLMMechanicalContact::validParams();
36 354 : params.addClassDescription("Computes the tangential frictional forces");
37 708 : params.addRequiredCoupledVar("friction_lm", "The frictional Lagrange's multiplier");
38 708 : params.addCoupledVar("friction_lm_dir",
39 : "The frictional Lagrange's multiplier for an addtional direction.");
40 708 : params.addParam<FunctionName>(
41 : "function_friction",
42 : "Coupled function to evaluate friction with values from contact pressure and relative "
43 : "tangential velocities");
44 708 : params.addParam<Real>("c_t", 1e0, "Numerical parameter for tangential constraints");
45 1062 : params.addRangeCheckedParam<Real>(
46 : "epsilon",
47 708 : 1.0e-7,
48 : "epsilon > 0",
49 : "Minimum value of contact pressure that will trigger frictional enforcement");
50 708 : params.addRangeCheckedParam<Real>(
51 : "mu", "mu >= 0", "The friction coefficient for the Coulomb friction law");
52 708 : params.addParam<MooseEnum>("friction_coefficient_regularization",
53 708 : Moose::Contact::frictionCoefficientRegularizationOptions(),
54 : "The regularization applied to the Coulomb friction coefficient.");
55 1062 : params.addRangeCheckedParam<Real>(
56 : "friction_reference_slip",
57 708 : 0.0,
58 : "friction_reference_slip >= 0",
59 : "Reference slip increment used by friction coefficient regularization.");
60 1062 : params.addRangeCheckedParam<Real>(
61 : "friction_elastic_slip",
62 708 : 0.0,
63 : "friction_elastic_slip >= 0",
64 : "Tangential elastic slip distance over which the Coulomb friction bound is reached.");
65 708 : params.addRequiredParam<UserObjectName>("weighted_velocities_uo",
66 : "The weighted tangential velocities user object.");
67 :
68 354 : return params;
69 0 : }
70 :
71 177 : ComputeFrictionalForceLMMechanicalContact::ComputeFrictionalForceLMMechanicalContact(
72 177 : const InputParameters & parameters)
73 : : ComputeWeightedGapLMMechanicalContact(parameters),
74 177 : _weighted_velocities_uo(getUserObject<WeightedVelocitiesUserObject>("weighted_velocities_uo")),
75 354 : _c_t(getParam<Real>("c_t")),
76 177 : _secondary_x_dot(_secondary_var.adUDot()),
77 177 : _primary_x_dot(_primary_var.adUDotNeighbor()),
78 177 : _secondary_y_dot(adCoupledDot("disp_y")),
79 177 : _primary_y_dot(adCoupledNeighborValueDot("disp_y")),
80 177 : _secondary_z_dot(_has_disp_z ? &adCoupledDot("disp_z") : nullptr),
81 177 : _primary_z_dot(_has_disp_z ? &adCoupledNeighborValueDot("disp_z") : nullptr),
82 354 : _epsilon(getParam<Real>("epsilon")),
83 524 : _mu(isParamValid("function_friction") ? std::numeric_limits<double>::quiet_NaN()
84 517 : : getParam<Real>("mu")),
85 177 : _friction_coefficient_regularization(
86 177 : getParam<MooseEnum>("friction_coefficient_regularization")
87 : .getEnum<Moose::Contact::FrictionCoefficientRegularization>()),
88 354 : _friction_reference_slip(getParam<Real>("friction_reference_slip")),
89 354 : _friction_elastic_slip(getParam<Real>("friction_elastic_slip")),
90 361 : _function_friction(isParamValid("function_friction") ? &getFunction("function_friction")
91 : : nullptr),
92 354 : _has_friction_function(isParamValid("function_friction")),
93 177 : _3d(_has_disp_z)
94 :
95 : {
96 354 : if (parameters.isParamSetByUser("mu") && _has_friction_function)
97 0 : paramError(
98 : "mu",
99 : "Please only provide friction either as a function or as a constant value, but not both.");
100 354 : else if (!parameters.isParamSetByUser("mu") && !_has_friction_function)
101 0 : paramError("mu", "Please provide a value or a function for the coefficient of friction.");
102 :
103 354 : if (!getParam<bool>("use_displaced_mesh"))
104 0 : paramError("use_displaced_mesh",
105 : "'use_displaced_mesh' must be true for the "
106 : "ComputeFrictionalForceLMMechanicalContact object");
107 :
108 318 : if (_3d && !isParamValid("friction_lm_dir"))
109 0 : paramError("friction_lm_dir",
110 : "Three-dimensional mortar frictional contact simulations require an additional "
111 : "frictional Lagrange's multiplier to enforce a second tangential pressure");
112 :
113 177 : if (_friction_coefficient_regularization !=
114 16 : Moose::Contact::FrictionCoefficientRegularization::NONE &&
115 16 : _friction_reference_slip <= 0.0)
116 2 : paramError("friction_reference_slip",
117 : "A positive friction_reference_slip is required when "
118 : "friction_coefficient_regularization is not NONE.");
119 :
120 175 : _friction_vars.push_back(getVar("friction_lm", 0));
121 :
122 175 : if (_3d)
123 94 : _friction_vars.push_back(getVar("friction_lm_dir", 0));
124 :
125 175 : if (!_friction_vars[0]->isNodal())
126 0 : if (_friction_vars[0]->feType().order != static_cast<Order>(0))
127 0 : paramError(
128 : "friction_lm",
129 : "Frictional contact constraints only support elemental variables of CONSTANT order");
130 175 : }
131 :
132 : void
133 0 : ComputeFrictionalForceLMMechanicalContact::computeQpProperties()
134 : {
135 0 : }
136 :
137 : void
138 0 : ComputeFrictionalForceLMMechanicalContact::computeQpIProperties()
139 : {
140 0 : }
141 :
142 : void
143 42077 : ComputeFrictionalForceLMMechanicalContact::residualSetup()
144 : {
145 42077 : }
146 :
147 : void
148 20000 : ComputeFrictionalForceLMMechanicalContact::post()
149 : {
150 : const auto & dof_to_weighted_tangential_velocity =
151 20000 : _weighted_velocities_uo.dofToWeightedVelocities();
152 :
153 : const std::unordered_map<const DofObject *, std::pair<ADReal, Real>> & dof_to_weighted_gap =
154 20000 : _weighted_gap_uo.dofToWeightedGap();
155 :
156 : // Enforce frictional constraints
157 :
158 93878 : for (const auto & [dof_object, weighted_velocities_pr] : dof_to_weighted_tangential_velocity)
159 : {
160 73878 : if (dof_object->processor_id() != this->processor_id())
161 1557 : continue;
162 :
163 : const auto & [weighted_gap_pr, normalization] =
164 72321 : libmesh_map_find(dof_to_weighted_gap, dof_object);
165 72321 : _weighted_gap_ptr = &weighted_gap_pr;
166 72321 : _normalization_ptr = &normalization;
167 72321 : _tangential_vel_ptr[0] = &(weighted_velocities_pr[0]);
168 :
169 72321 : if (_3d)
170 : {
171 0 : _tangential_vel_ptr[1] = &(weighted_velocities_pr[1]);
172 0 : enforceConstraintOnDof3d(dof_object);
173 : }
174 : else
175 72321 : enforceConstraintOnDof(dof_object);
176 : }
177 20000 : }
178 :
179 : void
180 31824 : ComputeFrictionalForceLMMechanicalContact::incorrectEdgeDroppingPost(
181 : const std::unordered_set<const Node *> & inactive_lm_nodes)
182 : {
183 : const auto & dof_to_weighted_tangential_velocity =
184 31824 : _weighted_velocities_uo.dofToWeightedVelocities();
185 31824 : const auto & dof_to_weighted_gap = _weighted_gap_uo.dofToWeightedGap();
186 : // Enforce frictional complementarity constraints
187 151909 : for (const auto & [dof_object, weighted_velocities_pr] : dof_to_weighted_tangential_velocity)
188 : {
189 : // If node inactive, skip
190 237950 : if ((inactive_lm_nodes.find(static_cast<const Node *>(dof_object)) !=
191 120085 : inactive_lm_nodes.end()) ||
192 120069 : (dof_object->processor_id() != this->processor_id()))
193 2220 : continue;
194 :
195 117865 : _weighted_gap_ptr = &dof_to_weighted_gap.at(dof_object).first;
196 117865 : _normalization_ptr = &dof_to_weighted_gap.at(dof_object).second;
197 117865 : _tangential_vel_ptr[0] = &weighted_velocities_pr[0];
198 :
199 117865 : if (_3d)
200 : {
201 39890 : _tangential_vel_ptr[1] = &weighted_velocities_pr[1];
202 39890 : enforceConstraintOnDof3d(dof_object);
203 : }
204 : else
205 77975 : enforceConstraintOnDof(dof_object);
206 : }
207 31824 : }
208 :
209 : void
210 39890 : ComputeFrictionalForceLMMechanicalContact::enforceConstraintOnDof3d(const DofObject * const dof)
211 : {
212 : using std::max, std::sqrt;
213 :
214 39890 : ComputeWeightedGapLMMechanicalContact::enforceConstraintOnDof(dof);
215 :
216 : // Get normal LM
217 39890 : const auto normal_dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
218 39890 : const ADReal & weighted_gap = *_weighted_gap_ptr;
219 39890 : ADReal contact_pressure = (*_sys.currentSolution())(normal_dof_index);
220 : Moose::derivInsert(contact_pressure.derivatives(), normal_dof_index, 1.);
221 :
222 : // Get friction LMs
223 : std::array<const ADReal *, 2> & tangential_vel = _tangential_vel_ptr;
224 : std::array<dof_id_type, 2> friction_dof_indices;
225 : std::array<ADReal, 2> friction_lm_values;
226 :
227 : const unsigned int num_tangents = 2;
228 119670 : for (const auto i : make_range(num_tangents))
229 : {
230 79780 : friction_dof_indices[i] = dof->dof_number(_sys.number(), _friction_vars[i]->number(), 0);
231 79780 : friction_lm_values[i] = (*_sys.currentSolution())(friction_dof_indices[i]);
232 : Moose::derivInsert(friction_lm_values[i].derivatives(), friction_dof_indices[i], 1.);
233 : }
234 :
235 : // Get normalized c and c_t values (if normalization specified
236 39890 : const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
237 39890 : const Real c_t = _normalize_c ? _c_t / *_normalization_ptr : _c_t;
238 :
239 : // Compute the friction coefficient (constant or function)
240 : const auto & real_tangential_velocity =
241 39890 : libmesh_map_find(_weighted_velocities_uo.dofToRealVelocities(), dof);
242 : ADReal mu_ad = computeFrictionValue(
243 39890 : contact_pressure, real_tangential_velocity[0], real_tangential_velocity[1]);
244 :
245 : ADReal dof_residual;
246 : ADReal dof_residual_dir;
247 :
248 : // Primal-dual active set strategy (PDASS)
249 39890 : if (contact_pressure < _epsilon)
250 : {
251 11818 : dof_residual = friction_lm_values[0];
252 11818 : dof_residual_dir = friction_lm_values[1];
253 : }
254 : else
255 : {
256 : // Epsilon to avoid automatic differentiation singularity
257 : const Real epsilon_sqrt = 1.0e-48;
258 :
259 28072 : const auto lamdba_plus_cg = contact_pressure + c * weighted_gap;
260 28072 : const auto normal_bound = max(0.0, lamdba_plus_cg);
261 : const auto friction_bound = mu_ad * normal_bound;
262 : const auto tangential_compliance =
263 28072 : _friction_elastic_slip > 0.0 ? _friction_elastic_slip / (friction_bound + _epsilon) : 0.0;
264 :
265 : std::array<ADReal, 2> lambda_t_plus_ctu;
266 : lambda_t_plus_ctu[0] =
267 28072 : friction_lm_values[0] +
268 56144 : c_t * (*tangential_vel[0] * _dt - tangential_compliance * friction_lm_values[0]);
269 : lambda_t_plus_ctu[1] =
270 28072 : friction_lm_values[1] +
271 56144 : c_t * (*tangential_vel[1] * _dt - tangential_compliance * friction_lm_values[1]);
272 :
273 : const auto tangential_trial_norm =
274 28072 : sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
275 28072 : lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt);
276 :
277 28072 : const auto term_1_x = max(friction_bound, tangential_trial_norm) * friction_lm_values[0];
278 :
279 56144 : const auto term_1_y = max(friction_bound, tangential_trial_norm) * friction_lm_values[1];
280 :
281 : const auto term_2_x = friction_bound * lambda_t_plus_ctu[0];
282 :
283 : const auto term_2_y = friction_bound * lambda_t_plus_ctu[1];
284 :
285 28072 : dof_residual = term_1_x - term_2_x;
286 28072 : dof_residual_dir = term_1_y - term_2_y;
287 : }
288 :
289 39890 : addResidualsAndJacobian(_assembly,
290 79780 : std::array<ADReal, 1>{{dof_residual}},
291 79780 : std::array<dof_id_type, 1>{{friction_dof_indices[0]}},
292 : _friction_vars[0]->scalingFactor());
293 39890 : addResidualsAndJacobian(_assembly,
294 79780 : std::array<ADReal, 1>{{dof_residual_dir}},
295 79780 : std::array<dof_id_type, 1>{{friction_dof_indices[1]}},
296 : _friction_vars[1]->scalingFactor());
297 39890 : }
298 :
299 : void
300 150296 : ComputeFrictionalForceLMMechanicalContact::enforceConstraintOnDof(const DofObject * const dof)
301 : {
302 : using std::abs, std::max;
303 :
304 150296 : ComputeWeightedGapLMMechanicalContact::enforceConstraintOnDof(dof);
305 :
306 : // Get friction LM
307 150296 : const auto friction_dof_index = dof->dof_number(_sys.number(), _friction_vars[0]->number(), 0);
308 150296 : const ADReal & tangential_vel = *_tangential_vel_ptr[0];
309 150296 : ADReal friction_lm_value = (*_sys.currentSolution())(friction_dof_index);
310 : Moose::derivInsert(friction_lm_value.derivatives(), friction_dof_index, 1.);
311 :
312 : // Get normal LM
313 150296 : const auto normal_dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
314 150296 : const ADReal & weighted_gap = *_weighted_gap_ptr;
315 150296 : ADReal contact_pressure = (*_sys.currentSolution())(normal_dof_index);
316 : Moose::derivInsert(contact_pressure.derivatives(), normal_dof_index, 1.);
317 :
318 : // Get normalized c and c_t values (if normalization specified
319 150296 : const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
320 150296 : const Real c_t = _normalize_c ? _c_t / *_normalization_ptr : _c_t;
321 :
322 : // Compute the friction coefficient (constant or function)
323 : const auto & real_tangential_velocity =
324 150296 : libmesh_map_find(_weighted_velocities_uo.dofToRealVelocities(), dof);
325 150296 : ADReal mu_ad = computeFrictionValue(contact_pressure, real_tangential_velocity[0], 0.0);
326 :
327 : ADReal dof_residual;
328 : // Primal-dual active set strategy (PDASS)
329 150296 : if (contact_pressure < _epsilon)
330 42743 : dof_residual = friction_lm_value;
331 : else
332 : {
333 107553 : const auto lambda_plus_cg = contact_pressure + c * weighted_gap;
334 107553 : const auto normal_bound = max(0.0, lambda_plus_cg);
335 : const auto friction_bound = mu_ad * normal_bound;
336 : const auto tangential_compliance =
337 115001 : _friction_elastic_slip > 0.0 ? _friction_elastic_slip / (friction_bound + _epsilon) : 0.0;
338 : const auto lambda_t_plus_ctu =
339 : friction_lm_value +
340 107553 : c_t * (tangential_vel * _dt - tangential_compliance * friction_lm_value);
341 :
342 215106 : const auto term_1 = max(friction_bound, abs(lambda_t_plus_ctu)) * friction_lm_value;
343 : const auto term_2 = friction_bound * lambda_t_plus_ctu;
344 :
345 107553 : dof_residual = term_1 - term_2;
346 : }
347 :
348 150296 : addResidualsAndJacobian(_assembly,
349 300592 : std::array<ADReal, 1>{{dof_residual}},
350 300592 : std::array<dof_id_type, 1>{{friction_dof_index}},
351 : _friction_vars[0]->scalingFactor());
352 150296 : }
353 :
354 : ADReal
355 190186 : ComputeFrictionalForceLMMechanicalContact::computeFrictionValue(const ADReal & contact_pressure,
356 : const ADReal & tangential_vel,
357 : const ADReal & tangential_vel_dir)
358 : {
359 : using std::sqrt;
360 :
361 : // TODO: Introduce temperature dependence in the function. Do this when we have an example.
362 : ADReal mu_ad;
363 :
364 190186 : if (!_has_friction_function)
365 183578 : mu_ad = _mu;
366 : else
367 : {
368 : ADReal tangential_vel_magnitude =
369 6608 : sqrt(tangential_vel * tangential_vel + tangential_vel_dir * tangential_vel_dir + 1.0e-24);
370 6608 : mu_ad = _function_friction->value<ADReal>(0.0, contact_pressure, tangential_vel_magnitude, 0.0);
371 : }
372 :
373 : const auto slip_increment =
374 190186 : sqrt(tangential_vel * tangential_vel + tangential_vel_dir * tangential_vel_dir + 1.0e-24) *
375 190186 : _dt;
376 : return Moose::Contact::regularizedFrictionCoefficient(
377 190186 : mu_ad, slip_increment, _friction_coefficient_regularization, _friction_reference_slip);
378 : }
|
