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