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