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 "AutomaticMortarGeneration.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 <limits>
29 :
30 : registerMooseObject("ContactApp", ComputeDynamicFrictionalForceLMMechanicalContact);
31 :
32 : InputParameters
33 204 : ComputeDynamicFrictionalForceLMMechanicalContact::validParams()
34 : {
35 204 : InputParameters params = ComputeDynamicWeightedGapLMMechanicalContact::validParams();
36 204 : params.addClassDescription("Computes the tangential frictional forces for dynamic simulations");
37 408 : params.addRequiredCoupledVar("friction_lm", "The frictional Lagrange's multiplier");
38 408 : params.addCoupledVar("friction_lm_dir",
39 : "The frictional Lagrange's multiplier for an addtional direction.");
40 408 : params.addParam<FunctionName>(
41 : "function_friction",
42 : "Coupled function to evaluate friction with values from contact pressure and relative "
43 : "tangential velocities (from the previous step).");
44 408 : params.addParam<Real>("c_t", 1e0, "Numerical parameter for tangential constraints");
45 408 : params.addParam<Real>(
46 : "epsilon",
47 408 : 1.0e-7,
48 : "Minimum value of contact pressure that will trigger frictional enforcement");
49 408 : params.addParam<Real>("mu", "The friction coefficient for the Coulomb friction law");
50 204 : return params;
51 0 : }
52 :
53 102 : ComputeDynamicFrictionalForceLMMechanicalContact::ComputeDynamicFrictionalForceLMMechanicalContact(
54 102 : const InputParameters & parameters)
55 : : ComputeDynamicWeightedGapLMMechanicalContact(parameters),
56 204 : _c_t(getParam<Real>("c_t")),
57 102 : _secondary_x_dot(_secondary_var.adUDot()),
58 102 : _primary_x_dot(_primary_var.adUDotNeighbor()),
59 102 : _secondary_y_dot(adCoupledDot("disp_y")),
60 102 : _primary_y_dot(adCoupledNeighborValueDot("disp_y")),
61 102 : _secondary_z_dot(_has_disp_z ? &adCoupledDot("disp_z") : nullptr),
62 102 : _primary_z_dot(_has_disp_z ? &adCoupledNeighborValueDot("disp_z") : nullptr),
63 204 : _epsilon(getParam<Real>("epsilon")),
64 392 : _mu(isParamValid("mu") ? getParam<Real>("mu") : std::numeric_limits<double>::quiet_NaN()),
65 212 : _function_friction(isParamValid("function_friction") ? &getFunction("function_friction")
66 : : nullptr),
67 204 : _has_friction_function(isParamValid("function_friction")),
68 102 : _3d(_has_disp_z)
69 : {
70 384 : if (!_has_friction_function && !isParamValid("mu"))
71 0 : mooseError(
72 : "A coefficient of friction needs to be provided as a constant value of via a function.");
73 :
74 118 : if (_has_friction_function && isParamValid("mu"))
75 0 : paramError("mu",
76 : "Either provide a constant coefficient of friction or a function defining the "
77 : "coefficient of friction. Both inputs cannot be provided simultaneously.");
78 :
79 204 : if (!getParam<bool>("use_displaced_mesh"))
80 0 : paramError("use_displaced_mesh",
81 : "'use_displaced_mesh' must be true for the "
82 : "ComputeFrictionalForceLMMechanicalContact object");
83 :
84 204 : if (_3d && !isParamValid("friction_lm_dir"))
85 0 : paramError("friction_lm_dir",
86 : "Three-dimensional mortar frictional contact simulations require an additional "
87 : "frictional Lagrange's multiplier to enforce a second tangential pressure");
88 :
89 102 : _friction_vars.push_back(getVar("friction_lm", 0));
90 :
91 102 : if (_3d)
92 68 : _friction_vars.push_back(getVar("friction_lm_dir", 0));
93 :
94 102 : if (!_friction_vars[0]->isNodal())
95 0 : if (_friction_vars[0]->feType().order != static_cast<Order>(0))
96 0 : paramError(
97 : "friction_lm",
98 : "Frictional contact constraints only support elemental variables of CONSTANT order");
99 :
100 : // Request the old solution state in unison
101 102 : _sys.solutionOld();
102 102 : }
103 :
104 : void
105 686640 : ComputeDynamicFrictionalForceLMMechanicalContact::computeQpProperties()
106 : {
107 : // Compute the value of _qp_gap
108 686640 : ComputeDynamicWeightedGapLMMechanicalContact::computeQpProperties();
109 :
110 : // It appears that the relative velocity between weighted gap and this class have a sign
111 : // difference
112 686640 : _qp_tangential_velocity_nodal = -_relative_velocity * (_JxW_msm[_qp] * _coord[_qp]);
113 686640 : _qp_real_tangential_velocity_nodal = -_relative_velocity;
114 686640 : }
115 :
116 : void
117 2474080 : ComputeDynamicFrictionalForceLMMechanicalContact::computeQpIProperties()
118 : {
119 : // Get the _dof_to_weighted_gap map
120 2474080 : ComputeDynamicWeightedGapLMMechanicalContact::computeQpIProperties();
121 :
122 2474080 : const auto & nodal_tangents = amg().getNodalTangents(*_lower_secondary_elem);
123 :
124 : // Get the _dof_to_weighted_tangential_velocity map
125 : const DofObject * const dof =
126 2474080 : _friction_vars[0]->isNodal()
127 2474080 : ? static_cast<const DofObject *>(_lower_secondary_elem->node_ptr(_i))
128 0 : : static_cast<const DofObject *>(_lower_secondary_elem);
129 :
130 : _dof_to_weighted_tangential_velocity[dof][0] +=
131 4948160 : _test[_i][_qp] * _qp_tangential_velocity_nodal * nodal_tangents[0][_i];
132 :
133 2474080 : _dof_to_real_tangential_velocity[dof][0] +=
134 2474080 : _test[_i][_qp] * MetaPhysicL::raw_value(_qp_real_tangential_velocity_nodal) *
135 2474080 : nodal_tangents[0][_i];
136 :
137 : // Get the _dof_to_weighted_tangential_velocity map for a second direction
138 2474080 : if (_3d)
139 : {
140 : _dof_to_weighted_tangential_velocity[dof][1] +=
141 4403200 : _test[_i][_qp] * _qp_tangential_velocity_nodal * nodal_tangents[1][_i];
142 :
143 2201600 : _dof_to_real_tangential_velocity[dof][1] +=
144 2201600 : _test[_i][_qp] * MetaPhysicL::raw_value(_qp_real_tangential_velocity_nodal) *
145 2201600 : nodal_tangents[1][_i];
146 : }
147 2474080 : }
148 :
149 : void
150 6805 : ComputeDynamicFrictionalForceLMMechanicalContact::residualSetup()
151 : {
152 : // Clear both maps
153 6805 : ComputeDynamicWeightedGapLMMechanicalContact::residualSetup();
154 : _dof_to_weighted_tangential_velocity.clear();
155 : _dof_to_real_tangential_velocity.clear();
156 6805 : }
157 :
158 : void
159 339 : ComputeDynamicFrictionalForceLMMechanicalContact::timestepSetup()
160 : {
161 :
162 339 : ComputeDynamicWeightedGapLMMechanicalContact::timestepSetup();
163 :
164 : _dof_to_old_real_tangential_velocity.clear();
165 :
166 2129 : for (auto & map_pr : _dof_to_real_tangential_velocity)
167 : _dof_to_old_real_tangential_velocity.emplace(map_pr);
168 339 : }
169 :
170 : void
171 663 : ComputeDynamicFrictionalForceLMMechanicalContact::post()
172 : {
173 663 : ComputeDynamicWeightedGapLMMechanicalContact::post();
174 :
175 663 : Moose::Mortar::Contact::communicateVelocities(
176 663 : _dof_to_weighted_tangential_velocity, _mesh, _nodal, _communicator, false);
177 :
178 663 : if (_has_friction_function)
179 400 : Moose::Mortar::Contact::communicateVelocities(
180 400 : _dof_to_real_tangential_velocity, _mesh, _nodal, _communicator, false);
181 :
182 : // Enforce frictional complementarity constraints
183 4578 : for (const auto & pr : _dof_to_weighted_tangential_velocity)
184 : {
185 3915 : const DofObject * const dof = pr.first;
186 :
187 3915 : if (dof->processor_id() != this->processor_id())
188 0 : continue;
189 :
190 : // Use always weighted gap for dynamic PDASS. Omit the dynamic weighted gap approach that is
191 : // used in normal contact where the discretized gap velocity is enforced if a node has
192 : // identified to be into contact.
193 :
194 3915 : _weighted_gap_ptr = &_dof_to_weighted_gap[dof].first;
195 3915 : _normalization_ptr = &_dof_to_weighted_gap[dof].second;
196 3915 : _tangential_vel_ptr[0] = &(pr.second[0]);
197 :
198 3915 : if (_3d)
199 : {
200 3915 : _tangential_vel_ptr[1] = &(pr.second[1]);
201 3915 : enforceConstraintOnDof3d(dof);
202 : }
203 : else
204 0 : enforceConstraintOnDof(dof);
205 : }
206 663 : }
207 :
208 : void
209 6142 : ComputeDynamicFrictionalForceLMMechanicalContact::incorrectEdgeDroppingPost(
210 : const std::unordered_set<const Node *> & inactive_lm_nodes)
211 : {
212 6142 : ComputeDynamicWeightedGapLMMechanicalContact::incorrectEdgeDroppingPost(inactive_lm_nodes);
213 :
214 6142 : Moose::Mortar::Contact::communicateVelocities(
215 6142 : _dof_to_weighted_tangential_velocity, _mesh, _nodal, _communicator, false);
216 :
217 6142 : if (_has_friction_function)
218 0 : Moose::Mortar::Contact::communicateVelocities(
219 0 : _dof_to_real_tangential_velocity, _mesh, _nodal, _communicator, false);
220 :
221 : // Enforce frictional complementarity constraints
222 54072 : for (const auto & pr : _dof_to_weighted_tangential_velocity)
223 : {
224 47930 : const DofObject * const dof = pr.first;
225 :
226 : // If node inactive, skip
227 47930 : if ((inactive_lm_nodes.find(static_cast<const Node *>(dof)) != inactive_lm_nodes.end()) ||
228 47508 : (dof->processor_id() != this->processor_id()))
229 1214 : continue;
230 :
231 : // Use always weighted gap for dynamic PDASS
232 46716 : _weighted_gap_ptr = &_dof_to_weighted_gap[dof].first;
233 46716 : _normalization_ptr = &_dof_to_weighted_gap[dof].second;
234 46716 : _tangential_vel_ptr[0] = &pr.second[0];
235 :
236 46716 : if (_3d)
237 : {
238 9672 : _tangential_vel_ptr[1] = &pr.second[1];
239 9672 : enforceConstraintOnDof3d(dof);
240 : }
241 : else
242 37044 : enforceConstraintOnDof(dof);
243 : }
244 6142 : }
245 :
246 : void
247 13587 : ComputeDynamicFrictionalForceLMMechanicalContact::enforceConstraintOnDof3d(
248 : const DofObject * const dof)
249 : {
250 : using std::max, std::sqrt;
251 :
252 : // Get normal LM
253 13587 : const auto normal_dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
254 13587 : const ADReal & weighted_gap = *_weighted_gap_ptr;
255 13587 : ADReal contact_pressure = (*_sys.currentSolution())(normal_dof_index);
256 : Moose::derivInsert(contact_pressure.derivatives(), normal_dof_index, 1.);
257 :
258 : // Get friction LMs
259 : std::array<const ADReal *, 2> & tangential_vel = _tangential_vel_ptr;
260 : std::array<dof_id_type, 2> friction_dof_indices;
261 : std::array<ADReal, 2> friction_lm_values;
262 :
263 : const unsigned int num_tangents = 2;
264 40761 : for (const auto i : make_range(num_tangents))
265 : {
266 27174 : friction_dof_indices[i] = dof->dof_number(_sys.number(), _friction_vars[i]->number(), 0);
267 27174 : friction_lm_values[i] = (*_sys.currentSolution())(friction_dof_indices[i]);
268 : Moose::derivInsert(friction_lm_values[i].derivatives(), friction_dof_indices[i], 1.);
269 : }
270 :
271 : // Get normalized c and c_t values (if normalization specified
272 13587 : const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
273 13587 : const Real c_t = _normalize_c ? _c_t / *_normalization_ptr : _c_t;
274 :
275 13587 : const Real contact_pressure_old = _sys.solutionOld()(normal_dof_index);
276 :
277 : // Compute the friction coefficient (constant or function)
278 : ADReal mu_ad = computeFrictionValue(contact_pressure_old,
279 : _dof_to_old_real_tangential_velocity[dof][0],
280 13587 : _dof_to_old_real_tangential_velocity[dof][1]);
281 :
282 : ADReal dof_residual;
283 : ADReal dof_residual_dir;
284 :
285 : // Primal-dual active set strategy (PDASS)
286 13587 : if (contact_pressure < _epsilon)
287 : {
288 3536 : dof_residual = friction_lm_values[0];
289 3536 : dof_residual_dir = friction_lm_values[1];
290 : }
291 : else
292 : {
293 : const Real epsilon_sqrt = 1.0e-48;
294 :
295 10051 : const auto lamdba_plus_cg = contact_pressure + c * weighted_gap;
296 : std::array<ADReal, 2> lambda_t_plus_ctu;
297 20102 : lambda_t_plus_ctu[0] = friction_lm_values[0] + c_t * *tangential_vel[0] * _dt;
298 20102 : lambda_t_plus_ctu[1] = friction_lm_values[1] + c_t * *tangential_vel[1] * _dt;
299 :
300 10051 : const auto term_1_x = max(mu_ad * lamdba_plus_cg,
301 10051 : sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
302 10051 : lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt)) *
303 : friction_lm_values[0];
304 :
305 10051 : const auto term_1_y = max(mu_ad * lamdba_plus_cg,
306 10051 : sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
307 10051 : lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt)) *
308 : friction_lm_values[1];
309 :
310 10051 : const auto term_2_x = mu_ad * max(0.0, lamdba_plus_cg) * lambda_t_plus_ctu[0];
311 :
312 20102 : const auto term_2_y = mu_ad * max(0.0, lamdba_plus_cg) * lambda_t_plus_ctu[1];
313 :
314 10051 : dof_residual = term_1_x - term_2_x;
315 10051 : dof_residual_dir = term_1_y - term_2_y;
316 : }
317 :
318 13587 : addResidualsAndJacobian(_assembly,
319 27174 : std::array<ADReal, 1>{{dof_residual}},
320 27174 : std::array<dof_id_type, 1>{{friction_dof_indices[0]}},
321 : _friction_vars[0]->scalingFactor());
322 13587 : addResidualsAndJacobian(_assembly,
323 27174 : std::array<ADReal, 1>{{dof_residual_dir}},
324 27174 : std::array<dof_id_type, 1>{{friction_dof_indices[1]}},
325 : _friction_vars[1]->scalingFactor());
326 13587 : }
327 :
328 : void
329 37044 : ComputeDynamicFrictionalForceLMMechanicalContact::enforceConstraintOnDof(
330 : const DofObject * const dof)
331 : {
332 : using std::max, std::abs;
333 :
334 : // Get friction LM
335 37044 : const auto friction_dof_index = dof->dof_number(_sys.number(), _friction_vars[0]->number(), 0);
336 37044 : const ADReal & tangential_vel = *_tangential_vel_ptr[0];
337 37044 : ADReal friction_lm_value = (*_sys.currentSolution())(friction_dof_index);
338 : Moose::derivInsert(friction_lm_value.derivatives(), friction_dof_index, 1.);
339 :
340 : // Get normal LM
341 37044 : const auto normal_dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
342 37044 : const ADReal & weighted_gap = *_weighted_gap_ptr;
343 37044 : ADReal contact_pressure = (*_sys.currentSolution())(normal_dof_index);
344 : Moose::derivInsert(contact_pressure.derivatives(), normal_dof_index, 1.);
345 :
346 37044 : const Real contact_pressure_old = _sys.solutionOld()(normal_dof_index);
347 :
348 : // Get normalized c and c_t values (if normalization specified
349 37044 : const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
350 37044 : const Real c_t = _normalize_c ? _c_t / *_normalization_ptr : _c_t;
351 :
352 : // Compute the friction coefficient (constant or function)
353 : ADReal mu_ad =
354 37044 : computeFrictionValue(contact_pressure_old, _dof_to_old_real_tangential_velocity[dof][0], 0.0);
355 :
356 : ADReal dof_residual;
357 : // Primal-dual active set strategy (PDASS)
358 37044 : if (contact_pressure < _epsilon)
359 22326 : dof_residual = friction_lm_value;
360 : else
361 : {
362 14718 : const auto term_1 = max(mu_ad * (contact_pressure + c * weighted_gap),
363 29436 : abs(friction_lm_value + c_t * tangential_vel * _dt)) *
364 : friction_lm_value;
365 14718 : const auto term_2 = mu_ad * max(0.0, contact_pressure + c * weighted_gap) *
366 29436 : (friction_lm_value + c_t * tangential_vel * _dt);
367 :
368 14718 : dof_residual = term_1 - term_2;
369 : }
370 :
371 37044 : addResidualsAndJacobian(_assembly,
372 74088 : std::array<ADReal, 1>{{dof_residual}},
373 74088 : std::array<dof_id_type, 1>{{friction_dof_index}},
374 : _friction_vars[0]->scalingFactor());
375 37044 : }
376 :
377 : ADReal
378 50631 : ComputeDynamicFrictionalForceLMMechanicalContact::computeFrictionValue(
379 : const ADReal & contact_pressure, const Real & tangential_vel, const Real & tangential_vel_dir)
380 : {
381 : using std::sqrt;
382 :
383 : // TODO: Introduce temperature dependence in the function. Do this when we have an example.
384 : ADReal mu_ad;
385 :
386 50631 : if (!_has_friction_function)
387 48381 : mu_ad = _mu;
388 : else
389 : {
390 : ADReal tangential_vel_magnitude =
391 2250 : sqrt(tangential_vel * tangential_vel + tangential_vel_dir * tangential_vel_dir + 1.0e-24);
392 :
393 2250 : mu_ad = _function_friction->value<ADReal>(0.0, contact_pressure, tangential_vel_magnitude, 0.0);
394 : }
395 :
396 50631 : return mu_ad;
397 : }
|
