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 "WeightedVelocitiesUserObject.h"
15 :
16 : #include "metaphysicl/metaphysicl_version.h"
17 : #include "metaphysicl/dualsemidynamicsparsenumberarray.h"
18 : #include "metaphysicl/parallel_dualnumber.h"
19 : #if METAPHYSICL_MAJOR_VERSION < 2
20 : #include "metaphysicl/parallel_dynamic_std_array_wrapper.h"
21 : #else
22 : #include "metaphysicl/parallel_dynamic_array_wrapper.h"
23 : #endif
24 : #include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
25 : #include "timpi/parallel_sync.h"
26 :
27 : #include <cmath>
28 :
29 : registerMooseObject("ContactApp", ComputeFrictionalForceLMMechanicalContact);
30 :
31 : InputParameters
32 322 : ComputeFrictionalForceLMMechanicalContact::validParams()
33 : {
34 322 : InputParameters params = ComputeWeightedGapLMMechanicalContact::validParams();
35 322 : params.addClassDescription("Computes the tangential frictional forces");
36 644 : params.addRequiredCoupledVar("friction_lm", "The frictional Lagrange's multiplier");
37 644 : params.addCoupledVar("friction_lm_dir",
38 : "The frictional Lagrange's multiplier for an addtional direction.");
39 644 : params.addParam<FunctionName>(
40 : "function_friction",
41 : "Coupled function to evaluate friction with values from contact pressure and relative "
42 : "tangential velocities");
43 644 : params.addParam<Real>("c_t", 1e0, "Numerical parameter for tangential constraints");
44 644 : params.addParam<Real>(
45 : "epsilon",
46 644 : 1.0e-7,
47 : "Minimum value of contact pressure that will trigger frictional enforcement");
48 644 : params.addRangeCheckedParam<Real>(
49 : "mu", "mu > 0", "The friction coefficient for the Coulomb friction law");
50 644 : params.addRequiredParam<UserObjectName>("weighted_velocities_uo",
51 : "The weighted tangential velocities user object.");
52 :
53 322 : return params;
54 0 : }
55 :
56 161 : ComputeFrictionalForceLMMechanicalContact::ComputeFrictionalForceLMMechanicalContact(
57 161 : const InputParameters & parameters)
58 : : ComputeWeightedGapLMMechanicalContact(parameters),
59 161 : _weighted_velocities_uo(getUserObject<WeightedVelocitiesUserObject>("weighted_velocities_uo")),
60 322 : _c_t(getParam<Real>("c_t")),
61 161 : _secondary_x_dot(_secondary_var.adUDot()),
62 161 : _primary_x_dot(_primary_var.adUDotNeighbor()),
63 161 : _secondary_y_dot(adCoupledDot("disp_y")),
64 161 : _primary_y_dot(adCoupledNeighborValueDot("disp_y")),
65 161 : _secondary_z_dot(_has_disp_z ? &adCoupledDot("disp_z") : nullptr),
66 161 : _primary_z_dot(_has_disp_z ? &adCoupledNeighborValueDot("disp_z") : nullptr),
67 322 : _epsilon(getParam<Real>("epsilon")),
68 476 : _mu(isParamValid("function_friction") ? std::numeric_limits<double>::quiet_NaN()
69 469 : : getParam<Real>("mu")),
70 329 : _function_friction(isParamValid("function_friction") ? &getFunction("function_friction")
71 : : nullptr),
72 322 : _has_friction_function(isParamValid("function_friction")),
73 161 : _3d(_has_disp_z)
74 :
75 : {
76 161 : _weighted_velocities_uo.includeNodalNormalDerivatives();
77 :
78 322 : if (parameters.isParamSetByUser("mu") && _has_friction_function)
79 0 : paramError(
80 : "mu",
81 : "Please only provide friction either as a function or as a constant value, but not both.");
82 322 : else if (!parameters.isParamSetByUser("mu") && !_has_friction_function)
83 0 : paramError("mu", "Please provide a value or a function for the coefficient of friction.");
84 :
85 322 : if (!getParam<bool>("use_displaced_mesh"))
86 0 : paramError("use_displaced_mesh",
87 : "'use_displaced_mesh' must be true for the "
88 : "ComputeFrictionalForceLMMechanicalContact object");
89 :
90 302 : if (_3d && !isParamValid("friction_lm_dir"))
91 0 : paramError("friction_lm_dir",
92 : "Three-dimensional mortar frictional contact simulations require an additional "
93 : "frictional Lagrange's multiplier to enforce a second tangential pressure");
94 :
95 161 : _friction_vars.push_back(getVar("friction_lm", 0));
96 :
97 161 : if (_3d)
98 94 : _friction_vars.push_back(getVar("friction_lm_dir", 0));
99 :
100 161 : if (!_friction_vars[0]->isNodal())
101 0 : if (_friction_vars[0]->feType().order != static_cast<Order>(0))
102 0 : paramError(
103 : "friction_lm",
104 : "Frictional contact constraints only support elemental variables of CONSTANT order");
105 161 : }
106 :
107 : void
108 0 : ComputeFrictionalForceLMMechanicalContact::computeQpProperties()
109 : {
110 0 : }
111 :
112 : void
113 0 : ComputeFrictionalForceLMMechanicalContact::computeQpIProperties()
114 : {
115 0 : }
116 :
117 : void
118 36698 : ComputeFrictionalForceLMMechanicalContact::residualSetup()
119 : {
120 36698 : }
121 :
122 : void
123 18533 : ComputeFrictionalForceLMMechanicalContact::post()
124 : {
125 : const auto & dof_to_weighted_tangential_velocity =
126 18533 : _weighted_velocities_uo.dofToWeightedVelocities();
127 :
128 : const std::unordered_map<const DofObject *, std::pair<ADReal, Real>> & dof_to_weighted_gap =
129 18533 : _weighted_gap_uo.dofToWeightedGap();
130 :
131 : // Enforce frictional constraints
132 :
133 89478 : for (const auto & [dof_object, weighted_velocities_pr] : dof_to_weighted_tangential_velocity)
134 : {
135 70945 : if (dof_object->processor_id() != this->processor_id())
136 1422 : continue;
137 :
138 : const auto & [weighted_gap_pr, normalization] =
139 69523 : libmesh_map_find(dof_to_weighted_gap, dof_object);
140 69523 : _weighted_gap_ptr = &weighted_gap_pr;
141 69523 : _normalization_ptr = &normalization;
142 69523 : _tangential_vel_ptr[0] = &(weighted_velocities_pr[0]);
143 :
144 69523 : if (_3d)
145 : {
146 0 : _tangential_vel_ptr[1] = &(weighted_velocities_pr[1]);
147 0 : enforceConstraintOnDof3d(dof_object);
148 : }
149 : else
150 69523 : enforceConstraintOnDof(dof_object);
151 : }
152 18533 : }
153 :
154 : void
155 27085 : ComputeFrictionalForceLMMechanicalContact::incorrectEdgeDroppingPost(
156 : const std::unordered_set<const Node *> & inactive_lm_nodes)
157 : {
158 : const auto & dof_to_weighted_tangential_velocity =
159 27085 : _weighted_velocities_uo.dofToWeightedVelocities();
160 27085 : const auto & dof_to_weighted_gap = _weighted_gap_uo.dofToWeightedGap();
161 : // Enforce frictional complementarity constraints
162 135314 : for (const auto & [dof_object, weighted_velocities_pr] : dof_to_weighted_tangential_velocity)
163 : {
164 : // If node inactive, skip
165 214308 : if ((inactive_lm_nodes.find(static_cast<const Node *>(dof_object)) !=
166 108229 : inactive_lm_nodes.end()) ||
167 108221 : (dof_object->processor_id() != this->processor_id()))
168 2150 : continue;
169 :
170 106079 : _weighted_gap_ptr = &dof_to_weighted_gap.at(dof_object).first;
171 106079 : _normalization_ptr = &dof_to_weighted_gap.at(dof_object).second;
172 106079 : _tangential_vel_ptr[0] = &weighted_velocities_pr[0];
173 :
174 106079 : if (_3d)
175 : {
176 38454 : _tangential_vel_ptr[1] = &weighted_velocities_pr[1];
177 38454 : enforceConstraintOnDof3d(dof_object);
178 : }
179 : else
180 67625 : enforceConstraintOnDof(dof_object);
181 : }
182 27085 : }
183 :
184 : void
185 38454 : ComputeFrictionalForceLMMechanicalContact::enforceConstraintOnDof3d(const DofObject * const dof)
186 : {
187 : using std::max, std::sqrt;
188 :
189 38454 : ComputeWeightedGapLMMechanicalContact::enforceConstraintOnDof(dof);
190 :
191 : // Get normal LM
192 38454 : const auto normal_dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
193 38454 : const ADReal & weighted_gap = *_weighted_gap_ptr;
194 38454 : ADReal contact_pressure = (*_sys.currentSolution())(normal_dof_index);
195 : Moose::derivInsert(contact_pressure.derivatives(), normal_dof_index, 1.);
196 :
197 : // Get friction LMs
198 : std::array<const ADReal *, 2> & tangential_vel = _tangential_vel_ptr;
199 : std::array<dof_id_type, 2> friction_dof_indices;
200 : std::array<ADReal, 2> friction_lm_values;
201 :
202 : const unsigned int num_tangents = 2;
203 115362 : for (const auto i : make_range(num_tangents))
204 : {
205 76908 : friction_dof_indices[i] = dof->dof_number(_sys.number(), _friction_vars[i]->number(), 0);
206 76908 : friction_lm_values[i] = (*_sys.currentSolution())(friction_dof_indices[i]);
207 : Moose::derivInsert(friction_lm_values[i].derivatives(), friction_dof_indices[i], 1.);
208 : }
209 :
210 : // Get normalized c and c_t values (if normalization specified
211 38454 : const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
212 38454 : const Real c_t = _normalize_c ? _c_t / *_normalization_ptr : _c_t;
213 :
214 : // Compute the friction coefficient (constant or function)
215 : ADReal mu_ad = computeFrictionValue(contact_pressure,
216 : _dof_to_real_tangential_velocity[dof][0],
217 38454 : _dof_to_real_tangential_velocity[dof][1]);
218 :
219 : ADReal dof_residual;
220 : ADReal dof_residual_dir;
221 :
222 : // Primal-dual active set strategy (PDASS)
223 38454 : if (contact_pressure < _epsilon)
224 : {
225 11290 : dof_residual = friction_lm_values[0];
226 11290 : dof_residual_dir = friction_lm_values[1];
227 : }
228 : else
229 : {
230 : // Espilon to avoid automatic differentiation singularity
231 : const Real epsilon_sqrt = 1.0e-48;
232 :
233 27164 : const auto lamdba_plus_cg = contact_pressure + c * weighted_gap;
234 : std::array<ADReal, 2> lambda_t_plus_ctu;
235 54328 : lambda_t_plus_ctu[0] = friction_lm_values[0] + c_t * *tangential_vel[0] * _dt;
236 54328 : lambda_t_plus_ctu[1] = friction_lm_values[1] + c_t * *tangential_vel[1] * _dt;
237 :
238 27164 : const auto term_1_x = max(mu_ad * lamdba_plus_cg,
239 27164 : sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
240 27164 : lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt)) *
241 : friction_lm_values[0];
242 :
243 27164 : const auto term_1_y = max(mu_ad * lamdba_plus_cg,
244 27164 : sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
245 27164 : lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt)) *
246 : friction_lm_values[1];
247 :
248 27164 : const auto term_2_x = mu_ad * max(0.0, lamdba_plus_cg) * lambda_t_plus_ctu[0];
249 :
250 54328 : const auto term_2_y = mu_ad * max(0.0, lamdba_plus_cg) * lambda_t_plus_ctu[1];
251 :
252 27164 : dof_residual = term_1_x - term_2_x;
253 27164 : dof_residual_dir = term_1_y - term_2_y;
254 : }
255 :
256 38454 : addResidualsAndJacobian(_assembly,
257 76908 : std::array<ADReal, 1>{{dof_residual}},
258 76908 : std::array<dof_id_type, 1>{{friction_dof_indices[0]}},
259 : _friction_vars[0]->scalingFactor());
260 38454 : addResidualsAndJacobian(_assembly,
261 76908 : std::array<ADReal, 1>{{dof_residual_dir}},
262 76908 : std::array<dof_id_type, 1>{{friction_dof_indices[1]}},
263 : _friction_vars[1]->scalingFactor());
264 38454 : }
265 :
266 : void
267 137148 : ComputeFrictionalForceLMMechanicalContact::enforceConstraintOnDof(const DofObject * const dof)
268 : {
269 : using std::abs, std::max;
270 :
271 137148 : ComputeWeightedGapLMMechanicalContact::enforceConstraintOnDof(dof);
272 :
273 : // Get friction LM
274 137148 : const auto friction_dof_index = dof->dof_number(_sys.number(), _friction_vars[0]->number(), 0);
275 137148 : const ADReal & tangential_vel = *_tangential_vel_ptr[0];
276 137148 : ADReal friction_lm_value = (*_sys.currentSolution())(friction_dof_index);
277 : Moose::derivInsert(friction_lm_value.derivatives(), friction_dof_index, 1.);
278 :
279 : // Get normal LM
280 137148 : const auto normal_dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
281 137148 : const ADReal & weighted_gap = *_weighted_gap_ptr;
282 137148 : ADReal contact_pressure = (*_sys.currentSolution())(normal_dof_index);
283 : Moose::derivInsert(contact_pressure.derivatives(), normal_dof_index, 1.);
284 :
285 : // Get normalized c and c_t values (if normalization specified
286 137148 : const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
287 137148 : const Real c_t = _normalize_c ? _c_t / *_normalization_ptr : _c_t;
288 :
289 : // Compute the friction coefficient (constant or function)
290 : ADReal mu_ad =
291 274296 : computeFrictionValue(contact_pressure, _dof_to_real_tangential_velocity[dof][0], 0.0);
292 :
293 : ADReal dof_residual;
294 : // Primal-dual active set strategy (PDASS)
295 137148 : if (contact_pressure < _epsilon)
296 41147 : dof_residual = friction_lm_value;
297 : else
298 : {
299 96001 : const auto term_1 = max(mu_ad * (contact_pressure + c * weighted_gap),
300 192002 : abs(friction_lm_value + c_t * tangential_vel * _dt)) *
301 : friction_lm_value;
302 96001 : const auto term_2 = mu_ad * max(0.0, contact_pressure + c * weighted_gap) *
303 192002 : (friction_lm_value + c_t * tangential_vel * _dt);
304 :
305 96001 : dof_residual = term_1 - term_2;
306 : }
307 :
308 137148 : addResidualsAndJacobian(_assembly,
309 274296 : std::array<ADReal, 1>{{dof_residual}},
310 274296 : std::array<dof_id_type, 1>{{friction_dof_index}},
311 : _friction_vars[0]->scalingFactor());
312 137148 : }
313 :
314 : ADReal
315 175602 : ComputeFrictionalForceLMMechanicalContact::computeFrictionValue(const ADReal & contact_pressure,
316 : const ADReal & tangential_vel,
317 : const ADReal & tangential_vel_dir)
318 : {
319 : using std::sqrt;
320 :
321 : // TODO: Introduce temperature dependence in the function. Do this when we have an example.
322 : ADReal mu_ad;
323 :
324 175602 : if (!_has_friction_function)
325 170034 : mu_ad = _mu;
326 : else
327 : {
328 : ADReal tangential_vel_magnitude =
329 5568 : sqrt(tangential_vel * tangential_vel + tangential_vel_dir * tangential_vel_dir + 1.0e-24);
330 5568 : mu_ad = _function_friction->value<ADReal>(0.0, contact_pressure, tangential_vel_magnitude, 0.0);
331 : }
332 :
333 175602 : return mu_ad;
334 : }
|
