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