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 "CohesiveZoneModelBase.h"
11 : #include "MooseVariableFE.h"
12 : #include "SystemBase.h"
13 : #include "MortarUtils.h"
14 : #include "MooseUtils.h"
15 : #include "MathUtils.h"
16 :
17 : #include "MortarContactUtils.h"
18 : #include "FactorizedRankTwoTensor.h"
19 :
20 : #include "ADReal.h"
21 :
22 : #include "CohesiveZoneModelTools.h"
23 :
24 : #include <Eigen/Core>
25 :
26 : InputParameters
27 90 : CohesiveZoneModelBase::validParams()
28 : {
29 90 : InputParameters params = WeightedVelocitiesUserObject::validParams();
30 90 : params += PenaltyWeightedGapUserObject::validParams();
31 180 : params.addRequiredCoupledVar("displacements",
32 : "The string of displacements suitable for the problem statement");
33 90 : params.addClassDescription(
34 : "Base class for mortar-based cohesive zone model. To handle frictional cohesive interfaces, "
35 : "it computes the mortar frictional contact forces using a penalty approach.");
36 180 : params.addParam<Real>(
37 : "penalty_friction",
38 : "The penalty factor for frictional interaction. If not provided, the normal "
39 : "penalty factor is also used for the frictional problem.");
40 180 : params.addParam<Real>(
41 180 : "friction_coefficient", 0.0, "The friction coefficient ruling Coulomb friction equations.");
42 : // Suppress augmented Lagrange parameters. AL implementation for CZM remains to be done.
43 90 : params.suppressParameter<Real>("max_penalty_multiplier");
44 90 : params.suppressParameter<Real>("penalty_multiplier");
45 90 : params.suppressParameter<Real>("penetration_tolerance");
46 90 : params.setParameters("penalty", 0.0);
47 90 : params.set<bool>("allow_nodal_normal_derivatives") = false;
48 90 : params.set<bool>("use_nodal_normal_derivatives") = false;
49 90 : params.set<bool>("ghost_point_neighbors") = false;
50 90 : return params;
51 0 : }
52 :
53 45 : CohesiveZoneModelBase::CohesiveZoneModelBase(const InputParameters & parameters)
54 : /*
55 : * We are using virtual inheritance to avoid the "Diamond inheritance" problem. This means that
56 : * we have to construct WeightedGapUserObject explicitly as it will _not_ be constructed in
57 : * the intermediate base classes PenaltyWeightedGapUserObject and WeightedVelocitiesUserObject.
58 : * Virtual inheritance ensures that only one instance of WeightedGapUserObject is included in this
59 : * class. The inheritance diagram is as follows:
60 : *
61 : * WeightedGapUserObject <----- PenaltyWeightedGapUserObject
62 : * ^ ^
63 : * | |
64 : * WeightedVelocitiesUserObject <----- CohesiveZoneModelBase
65 : *
66 : */
67 : : WeightedGapUserObject(parameters),
68 : PenaltyWeightedGapUserObject(parameters),
69 : WeightedVelocitiesUserObject(parameters),
70 45 : _ndisp(coupledComponents("displacements")),
71 180 : _penalty_friction(isParamValid("penalty_friction") ? getParam<Real>("penalty_friction")
72 45 : : getParam<Real>("penalty")),
73 90 : _friction_coefficient(getParam<Real>("friction_coefficient")),
74 45 : _dof_to_accumulated_slip(
75 45 : declareRestartableData<std::unordered_map<dof_id_type, std::pair<ADTwoVector, TwoVector>>>(
76 : "dof_to_accumulated_slip",
77 45 : std::unordered_map<dof_id_type, std::pair<ADTwoVector, TwoVector>>{})),
78 45 : _dof_to_tangential_traction(
79 45 : declareRestartableData<std::unordered_map<dof_id_type, std::pair<ADTwoVector, TwoVector>>>(
80 : "dof_to_tangential_traction",
81 90 : std::unordered_map<dof_id_type, std::pair<ADTwoVector, TwoVector>>{})),
82 45 : _epsilon_tolerance(1.0e-40),
83 45 : _dof_to_damage(declareRestartableData<std::unordered_map<dof_id_type, std::pair<ADReal, Real>>>(
84 135 : "dof_do_damage", std::unordered_map<dof_id_type, std::pair<ADReal, Real>>{}))
85 : {
86 45 : _czm_interpolated_traction.resize(_ndisp);
87 :
88 135 : for (unsigned int i = 0; i < _ndisp; ++i)
89 : {
90 90 : _grad_disp.push_back(&adCoupledGradient("displacements", i));
91 180 : _grad_disp_neighbor.push_back(&adCoupledGradient("displacements", i));
92 : }
93 :
94 : // Set non-intervening components to zero
95 135 : for ([[maybe_unused]] const auto i : make_range(_ndisp))
96 : {
97 90 : _grad_disp.push_back(&_ad_grad_zero);
98 90 : _grad_disp_neighbor.push_back(&_ad_grad_zero);
99 : }
100 :
101 45 : if (_augmented_lagrange_problem)
102 0 : mooseError("CohesiveZoneModelBase constraints cannot be enforced with an augmented "
103 : "Lagrange approach.");
104 45 : }
105 :
106 : const VariableTestValue &
107 45 : CohesiveZoneModelBase::test() const
108 : {
109 45 : return _aux_lm_var ? _aux_lm_var->phiLower() : _disp_x_var->phiLower();
110 : }
111 :
112 : void
113 26886 : CohesiveZoneModelBase::computeQpProperties()
114 : {
115 26886 : WeightedVelocitiesUserObject::computeQpProperties();
116 :
117 : // Compute F and R.
118 26886 : const auto F = (ADRankTwoTensor::Identity() +
119 53772 : ADRankTwoTensor::initializeFromRows(
120 53772 : (*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]));
121 26886 : const auto F_neighbor = (ADRankTwoTensor::Identity() +
122 53772 : ADRankTwoTensor::initializeFromRows((*_grad_disp_neighbor[0])[_qp],
123 26886 : (*_grad_disp_neighbor[1])[_qp],
124 53772 : (*_grad_disp_neighbor[2])[_qp]));
125 :
126 : // TODO in follow-on PRs: Trim interior node variable derivatives
127 26886 : _F_interpolation = F * (_JxW_msm[_qp] * _coord[_qp]);
128 26886 : _F_neighbor_interpolation = F_neighbor * (_JxW_msm[_qp] * _coord[_qp]);
129 26886 : }
130 :
131 : void
132 53772 : CohesiveZoneModelBase::computeQpIProperties()
133 : {
134 53772 : WeightedVelocitiesUserObject::computeQpIProperties();
135 : // Get the _dof_to_weighted_gap map
136 53772 : const auto * const dof = static_cast<const DofObject *>(_lower_secondary_elem->node_ptr(_i));
137 :
138 : // TODO: Probably better to interpolate the deformation gradients.
139 107544 : _dof_to_F[dof] += (*_test)[_i][_qp] * _F_interpolation;
140 107544 : _dof_to_F_neighbor[dof] += (*_test)[_i][_qp] * _F_neighbor_interpolation;
141 53772 : }
142 :
143 : void
144 26886 : CohesiveZoneModelBase::computeFandR(const Node * const node)
145 : {
146 : // First call does not have maps available
147 26886 : const bool return_boolean = _dof_to_F.find(node) == _dof_to_F.end();
148 26886 : if (return_boolean)
149 0 : return;
150 :
151 26886 : const auto normalized_F = normalizeQuantity(_dof_to_F, node);
152 26886 : const auto normalized_F_neighbor = normalizeQuantity(_dof_to_F_neighbor, node);
153 :
154 : // This 'averaging' assumption below can probably be improved upon.
155 80658 : _dof_to_interface_F[node] = 0.5 * (normalized_F + normalized_F_neighbor);
156 :
157 107544 : for (const auto i : make_range(Moose::dim))
158 322632 : for (const auto j : make_range(Moose::dim))
159 241974 : if (!std::isfinite(MetaPhysicL::raw_value(normalized_F(i, j))))
160 : throw MooseException("The deformation gradient on the secondary surface is not finite in "
161 0 : "CohesiveZoneModelBase. MOOSE needs to cut the time step size.");
162 :
163 26886 : const auto dof_to_interface_F_node = libmesh_map_find(_dof_to_interface_F, node);
164 :
165 26886 : const ADFactorizedRankTwoTensor C = dof_to_interface_F_node.transpose() * dof_to_interface_F_node;
166 26886 : const auto Uinv = MathUtils::sqrt(C).inverse().get();
167 53772 : _dof_to_interface_R[node] = dof_to_interface_F_node * Uinv;
168 :
169 : // Transform interface jump according to two rotation matrices
170 26886 : const auto global_interface_displacement = _dof_to_interface_displacement_jump[node];
171 53772 : _dof_to_interface_displacement_jump[node] =
172 26886 : (_dof_to_interface_R[node] * _dof_to_rotation_matrix[node]).transpose() *
173 26886 : global_interface_displacement;
174 : }
175 :
176 : void
177 493 : CohesiveZoneModelBase::timestepSetup()
178 : {
179 : // these functions do not call WeightedGapUserObject::timestepSetup to avoid double initialization
180 493 : WeightedVelocitiesUserObject::selfTimestepSetup();
181 493 : PenaltyWeightedGapUserObject::selfTimestepSetup();
182 :
183 : // instead we call it explicitly here
184 493 : WeightedGapUserObject::timestepSetup();
185 :
186 : // Clear step slip (values used in between AL iterations for penalty adaptivity)
187 493 : for (auto & map_pr : _dof_to_step_slip)
188 : {
189 : auto & [step_slip, old_step_slip] = map_pr.second;
190 : old_step_slip = {0.0, 0.0};
191 0 : step_slip = {0.0, 0.0};
192 : }
193 :
194 : // save off accumulated slip from the last timestep
195 1758 : for (auto & map_pr : _dof_to_accumulated_slip)
196 : {
197 : auto & [accumulated_slip, old_accumulated_slip] = map_pr.second;
198 : old_accumulated_slip = MetaPhysicL::raw_value(accumulated_slip);
199 : }
200 :
201 493 : for (auto & dof_lp : _dof_to_local_penalty_friction)
202 0 : dof_lp.second = _penalty_friction;
203 :
204 : // save off tangential traction from the last timestep
205 1758 : for (auto & map_pr : _dof_to_tangential_traction)
206 : {
207 : auto & [tangential_traction, old_tangential_traction] = map_pr.second;
208 : old_tangential_traction = {MetaPhysicL::raw_value(tangential_traction(0)),
209 : MetaPhysicL::raw_value(tangential_traction(1))};
210 2530 : tangential_traction = {0.0, 0.0};
211 : }
212 :
213 493 : for (auto & [dof_object, delta_tangential_lm] : _dof_to_frictional_lagrange_multipliers)
214 : delta_tangential_lm.setZero();
215 :
216 : // save off damage from the last timestep
217 1758 : for (auto & map_pr : _dof_to_damage)
218 : {
219 : auto & [damage, old_damage] = map_pr.second;
220 1265 : old_damage = {MetaPhysicL::raw_value(damage)};
221 1265 : damage = {0.0};
222 : }
223 493 : }
224 :
225 : void
226 7087 : CohesiveZoneModelBase::initialize()
227 : {
228 : // these functions do not call WeightedGapUserObject::initialize to avoid double initialization
229 7087 : WeightedVelocitiesUserObject::selfInitialize();
230 7087 : PenaltyWeightedGapUserObject::selfInitialize();
231 :
232 : // instead we call it explicitly here
233 7087 : WeightedGapUserObject::initialize();
234 : _dof_to_F.clear();
235 : _dof_to_F_neighbor.clear();
236 : _dof_to_interface_displacement_jump.clear();
237 : _dof_to_interface_F.clear();
238 : _dof_to_interface_R.clear();
239 : _dof_to_czm_traction.clear();
240 21918 : for (auto & map_pr : _dof_to_rotation_matrix)
241 14831 : map_pr.second.setToIdentity();
242 7087 : }
243 :
244 : void
245 13443 : CohesiveZoneModelBase::reinit()
246 : {
247 : // Normal contact pressure with penalty
248 13443 : PenaltyWeightedGapUserObject::reinit();
249 :
250 : // Compute all rotations that were created as material properties in CZMComputeDisplacementJump
251 13443 : prepareJumpKinematicQuantities();
252 :
253 : // Reset frictional pressure
254 13443 : _frictional_contact_traction_one.resize(_qrule_msm->n_points());
255 13443 : _frictional_contact_traction_two.resize(_qrule_msm->n_points()); // 3D
256 :
257 40329 : for (const auto qp : make_range(_qrule_msm->n_points()))
258 : {
259 26886 : _frictional_contact_traction_one[qp] = 0.0;
260 26886 : _frictional_contact_traction_two[qp] = 0.0;
261 : }
262 :
263 : // zero vector
264 13443 : const static TwoVector zero{0.0, 0.0};
265 :
266 : // iterate over nodes
267 40329 : for (const auto i : make_range(_test->size()))
268 : {
269 : // current node
270 26886 : const Node * const node = _lower_secondary_elem->node_ptr(i);
271 :
272 : // Compute the weighted nodal deformation gradient and rotation tensors.
273 26886 : computeFandR(node);
274 :
275 : // The call below is a 'macro' call. Create a utility function or user object for it.
276 26886 : computeCZMTraction(node);
277 :
278 : // Build final traction vector
279 26886 : computeGlobalTraction(node);
280 :
281 : // Compute mechanical contact until end of method.
282 53772 : const auto penalty_friction = findValue(
283 26886 : _dof_to_local_penalty_friction, static_cast<const DofObject *>(node), _penalty_friction);
284 :
285 : // utilized quantities
286 26886 : const auto & normal_pressure = _dof_to_normal_pressure[node];
287 :
288 : // map the tangential traction and accumulated slip
289 26886 : auto & [tangential_traction, old_tangential_traction] = _dof_to_tangential_traction[node->id()];
290 26886 : auto & [accumulated_slip, old_accumulated_slip] = _dof_to_accumulated_slip[node->id()];
291 :
292 : Real normal_lm = -1;
293 26886 : if (auto it = _dof_to_lagrange_multiplier.find(node); it != _dof_to_lagrange_multiplier.end())
294 0 : normal_lm = it->second;
295 :
296 : // Keep active set fixed from second Uzawa loop
297 26886 : if (normal_lm < -TOLERANCE && normal_pressure > TOLERANCE)
298 : {
299 : using std::abs;
300 :
301 5358 : const auto & damage = _dof_to_damage[node->id()].first;
302 :
303 : const auto & real_tangential_velocity =
304 5358 : libmesh_map_find(_dof_to_real_tangential_velocity, node);
305 5358 : const ADTwoVector slip_distance = {real_tangential_velocity[0] * _dt,
306 5358 : real_tangential_velocity[1] * _dt};
307 :
308 : // frictional lagrange multiplier (delta lambda^(k)_T)
309 : const auto & tangential_lm =
310 5358 : _augmented_lagrange_problem ? _dof_to_frictional_lagrange_multipliers[node] : zero;
311 :
312 : // tangential trial traction (Simo 3.12)
313 : // Modified for implementation in MOOSE: Avoid pingponging on frictional sign (max. 0.4
314 : // capacity)
315 5358 : ADTwoVector inner_iteration_penalty_friction = penalty_friction * slip_distance;
316 :
317 : const auto slip_metric = std::abs(MetaPhysicL::raw_value(slip_distance).cwiseAbs()(0)) +
318 5358 : std::abs(MetaPhysicL::raw_value(slip_distance).cwiseAbs()(1));
319 :
320 0 : const auto slip_norm = (MetaPhysicL::raw_value(slip_distance)).norm();
321 5358 : const auto friction_limit = 0.4 * _friction_coefficient * damage * abs(normal_pressure);
322 :
323 5358 : if (slip_metric > _epsilon_tolerance && penalty_friction * slip_norm > friction_limit)
324 : {
325 : inner_iteration_penalty_friction =
326 5588 : MetaPhysicL::raw_value(friction_limit / (penalty_friction * slip_norm)) *
327 2794 : penalty_friction * slip_distance;
328 : }
329 :
330 : ADTwoVector tangential_trial_traction =
331 5358 : old_tangential_traction + tangential_lm + inner_iteration_penalty_friction;
332 :
333 : // Nonlinearity below
334 : ADReal tangential_trial_traction_norm =
335 5358 : (MetaPhysicL::raw_value(tangential_trial_traction)).norm();
336 10716 : ADReal phi_trial = tangential_trial_traction_norm - _friction_coefficient * normal_pressure;
337 : tangential_traction = tangential_trial_traction;
338 :
339 : // Simo considers this a 'return mapping'; we are just capping friction to the Coulomb limit.
340 5358 : if (phi_trial > 0.0)
341 : // Simo 3.13/3.14 (the penalty formulation has an error in the paper)
342 : tangential_traction -=
343 2439 : phi_trial * tangential_trial_traction / tangential_trial_traction_norm;
344 :
345 : // track accumulated slip for output purposes
346 5358 : accumulated_slip = old_accumulated_slip + slip_distance.cwiseAbs();
347 : }
348 : else
349 : {
350 : // reset slip and clear traction
351 21528 : accumulated_slip.setZero();
352 21528 : tangential_traction.setZero();
353 : }
354 :
355 : // Now that we have consistent nodal frictional values, create an interpolated frictional
356 : // pressure variable.
357 26886 : const auto test_i = (*_test)[i];
358 80658 : for (const auto qp : make_range(_qrule_msm->n_points()))
359 : {
360 107544 : _frictional_contact_traction_one[qp] += test_i[qp] * tangential_traction(0);
361 53772 : _frictional_contact_traction_two[qp] += test_i[qp] * tangential_traction(1);
362 : }
363 26886 : }
364 :
365 40329 : for (const auto i : index_range(_czm_interpolated_traction))
366 : {
367 26886 : _czm_interpolated_traction[i].resize(_qrule_msm->n_points());
368 :
369 80658 : for (const auto qp : make_range(_qrule_msm->n_points()))
370 53772 : _czm_interpolated_traction[i][qp] = 0.0;
371 : }
372 :
373 : // iterate over nodes
374 40329 : for (const auto i : make_range(_test->size()))
375 : {
376 : // current node
377 26886 : const Node * const node = _lower_secondary_elem->node_ptr(i);
378 :
379 : // End of CZM bilinear computations
380 26886 : auto it = _dof_to_czm_traction.find(node);
381 26886 : if (it == _dof_to_czm_traction.end())
382 : return;
383 :
384 26886 : const auto & test_i = (*_test)[i];
385 80658 : for (const auto qp : make_range(_qrule_msm->n_points()))
386 161316 : for (const auto idx : index_range(_czm_interpolated_traction))
387 215088 : _czm_interpolated_traction[idx][qp] += test_i[qp] * (it->second)(idx);
388 : }
389 : }
390 :
391 : void
392 13443 : CohesiveZoneModelBase::prepareJumpKinematicQuantities()
393 : {
394 : // Compute rotation matrix from secondary surface
395 : // Rotation matrices and local interface displacement jump.
396 40329 : for (const auto i : make_range(_test->size()))
397 : {
398 26886 : const Node * const node = _lower_secondary_elem->node_ptr(i);
399 :
400 : // First call does not have maps available
401 26886 : const bool return_boolean = _dof_to_weighted_gap.find(node) == _dof_to_weighted_gap.end();
402 26886 : if (return_boolean)
403 : return;
404 :
405 53772 : _dof_to_rotation_matrix[node] = CohesiveZoneModelTools::computeReferenceRotation<true>(
406 26886 : _normals[i], _subproblem.mesh().dimension());
407 :
408 : // Every time we used quantities from a map we "denormalize" it from the mortar integral.
409 : // See normalizing member functions.
410 26886 : if (auto it = _dof_to_weighted_displacements.find(node);
411 : it != _dof_to_weighted_displacements.end())
412 53772 : _dof_to_interface_displacement_jump[node] = it->second;
413 : }
414 : }
415 :
416 : void
417 7087 : CohesiveZoneModelBase::finalize()
418 : {
419 7087 : WeightedVelocitiesUserObject::finalize();
420 7087 : PenaltyWeightedGapUserObject::selfFinalize();
421 :
422 7087 : const bool send_data_back = !constrainedByOwner();
423 :
424 14174 : Moose::Mortar::Contact::communicateR2T(
425 7087 : _dof_to_F, _subproblem.mesh(), _nodal, _communicator, send_data_back);
426 :
427 14174 : Moose::Mortar::Contact::communicateR2T(
428 7087 : _dof_to_F_neighbor, _subproblem.mesh(), _nodal, _communicator, send_data_back);
429 :
430 14174 : Moose::Mortar::Contact::communicateRealObject(
431 7087 : _dof_to_weighted_displacements, _subproblem.mesh(), _nodal, _communicator, send_data_back);
432 7087 : }
433 :
434 : const ADVariableValue &
435 31856 : CohesiveZoneModelBase::contactTangentialPressureDirOne() const
436 : {
437 31856 : return _frictional_contact_traction_one;
438 : }
439 :
440 : const ADVariableValue &
441 0 : CohesiveZoneModelBase::contactTangentialPressureDirTwo() const
442 : {
443 0 : return _frictional_contact_traction_two;
444 : }
445 :
446 : void
447 26886 : CohesiveZoneModelBase::computeGlobalTraction(const Node * const node)
448 : {
449 : // First call does not have maps available
450 26886 : const auto it = _dof_to_czm_traction.find(node);
451 26886 : if (it == _dof_to_czm_traction.end())
452 0 : return;
453 :
454 : const auto local_traction_vector = it->second;
455 26886 : const auto rotation_matrix = libmesh_map_find(_dof_to_rotation_matrix, node);
456 :
457 53772 : _dof_to_czm_traction[node] = rotation_matrix * local_traction_vector;
458 : }
459 :
460 : const ADVariableValue &
461 215088 : CohesiveZoneModelBase::czmGlobalTraction(const unsigned int i) const
462 : {
463 : mooseAssert(i <= 3,
464 : "Internal error in czmGlobalTraction. Index exceeds the traction vector size.");
465 :
466 215088 : return _czm_interpolated_traction[i];
467 : }
|
