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