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 "PenaltySimpleCohesiveZoneModel.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             : registerMooseObject("ContactApp", PenaltySimpleCohesiveZoneModel);
      27             : 
      28             : InputParameters
      29         144 : PenaltySimpleCohesiveZoneModel::validParams()
      30             : {
      31         144 :   InputParameters params = WeightedVelocitiesUserObject::validParams();
      32         144 :   params += PenaltyWeightedGapUserObject::validParams();
      33             : 
      34         144 :   params.addClassDescription(
      35             :       "Computes the mortar frictional contact force via a penalty approach.");
      36         288 :   params.addParam<Real>("penalty_friction",
      37             :                         "The penalty factor for frictional interaction. If not provide, the normal "
      38             :                         "penalty factor is also used for the frictional problem.");
      39         288 :   params.addRequiredParam<Real>("friction_coefficient",
      40             :                                 "The friction coefficient ruling Coulomb friction equations.");
      41             : 
      42         144 :   return params;
      43           0 : }
      44             : 
      45          72 : PenaltySimpleCohesiveZoneModel::PenaltySimpleCohesiveZoneModel(const InputParameters & parameters)
      46             :   /*
      47             :    * We are using virtual inheritance to avoid the "Diamond inheritance" problem. This means that
      48             :    * that we have to construct WeightedGapUserObject explicitly as it will _not_ be constructed in
      49             :    * the intermediate base classes PenaltyWeightedGapUserObject and WeightedVelocitiesUserObject.
      50             :    * Virtual inheritance ensures that only one instance of WeightedGapUserObject is included in this
      51             :    * class. The inheritance diagram is as follows:
      52             :    *
      53             :    * WeightedGapUserObject         <-----   PenaltyWeightedGapUserObject
      54             :    *          ^                                         ^
      55             :    *          |                                         |
      56             :    * WeightedVelocitiesUserObject  <-----   PenaltySimpleCohesiveZoneModel
      57             :    *
      58             :    */
      59             :   : WeightedGapUserObject(parameters),
      60             :     PenaltyWeightedGapUserObject(parameters),
      61             :     WeightedVelocitiesUserObject(parameters),
      62          72 :     _penalty(getParam<Real>("penalty")),
      63         288 :     _penalty_friction(isParamValid("penalty_friction") ? getParam<Real>("penalty_friction")
      64          72 :                                                        : getParam<Real>("penalty")),
      65         144 :     _friction_coefficient(getParam<Real>("friction_coefficient")),
      66          72 :     _epsilon_tolerance(1.0e-40)
      67             : 
      68             : {
      69          72 :   if (_augmented_lagrange_problem)
      70           0 :     mooseError("PenaltySimpleCohesiveZoneModel constraints cannot be enforced with an augmented "
      71             :                "Lagrange approach.");
      72          72 : }
      73             : 
      74             : const VariableTestValue &
      75          72 : PenaltySimpleCohesiveZoneModel::test() const
      76             : {
      77          72 :   return _aux_lm_var ? _aux_lm_var->phiLower() : _disp_x_var->phiLower();
      78             : }
      79             : 
      80             : void
      81         585 : PenaltySimpleCohesiveZoneModel::timestepSetup()
      82             : {
      83             :   // these functions do not call WeightedGapUserObject::timestepSetup to avoid double initialization
      84         585 :   WeightedVelocitiesUserObject::selfTimestepSetup();
      85         585 :   PenaltyWeightedGapUserObject::selfTimestepSetup();
      86             : 
      87             :   // instead we call it explicitly here
      88         585 :   WeightedGapUserObject::timestepSetup();
      89             : 
      90             :   // save off accumulated slip from the last timestep
      91        1996 :   for (auto & map_pr : _dof_to_accumulated_slip)
      92             :   {
      93             :     auto & [accumulated_slip, old_accumulated_slip] = map_pr.second;
      94             :     old_accumulated_slip = accumulated_slip;
      95             :   }
      96             : 
      97         585 :   for (auto & dof_lp : _dof_to_local_penalty_friction)
      98           0 :     dof_lp.second = _penalty_friction;
      99             : 
     100             :   // save off tangential traction from the last timestep
     101        1996 :   for (auto & map_pr : _dof_to_tangential_traction)
     102             :   {
     103             :     auto & [tangential_traction, old_tangential_traction] = map_pr.second;
     104             :     old_tangential_traction = {MetaPhysicL::raw_value(tangential_traction(0)),
     105             :                                MetaPhysicL::raw_value(tangential_traction(1))};
     106        2822 :     tangential_traction = {0.0, 0.0};
     107             :   }
     108             : 
     109         585 :   for (auto & [dof_object, delta_tangential_lm] : _dof_to_frictional_lagrange_multipliers)
     110             :     delta_tangential_lm.setZero();
     111         585 : }
     112             : 
     113             : void
     114        7091 : PenaltySimpleCohesiveZoneModel::initialize()
     115             : {
     116             :   // these functions do not call WeightedGapUserObject::initialize to avoid double initialization
     117        7091 :   WeightedVelocitiesUserObject::selfInitialize();
     118        7091 :   PenaltyWeightedGapUserObject::selfInitialize();
     119             : 
     120             :   // instead we call it explicitly here
     121        7091 :   WeightedGapUserObject::initialize();
     122        7091 : }
     123             : 
     124             : void
     125       15034 : PenaltySimpleCohesiveZoneModel::reinit()
     126             : {
     127             :   // Normal contact pressure with penalty
     128       15034 :   PenaltyWeightedGapUserObject::reinit();
     129             : 
     130             :   // Compute all rotations that were created as material properties in CZMComputeDisplacementJump
     131       15034 :   prepareJumpKinematicQuantities();
     132             : 
     133             :   // Reset frictional pressure
     134       15034 :   _frictional_contact_traction_one.resize(_qrule_msm->n_points());
     135       15034 :   _frictional_contact_traction_two.resize(_qrule_msm->n_points()); // 3D
     136             : 
     137       45102 :   for (const auto qp : make_range(_qrule_msm->n_points()))
     138             :   {
     139       30068 :     _frictional_contact_traction_one[qp] = 0.0;
     140       30068 :     _frictional_contact_traction_two[qp] = 0.0;
     141             :   }
     142             : 
     143             :   // zero vector
     144       15034 :   const static TwoVector zero{0.0, 0.0};
     145             : 
     146             :   // iterate over nodes
     147       45102 :   for (const auto i : make_range(_test->size()))
     148             :   {
     149             :     // current node
     150       30068 :     const Node * const node = _lower_secondary_elem->node_ptr(i);
     151             : 
     152             :     // Compute the weighted nodal deformation gradient and rotation tensors.
     153       30068 :     computeFandR(node);
     154             : 
     155             :     // The call below is a 'macro' call. Create a utility function or user object for it.
     156       30068 :     computeBilinearMixedModeTraction(node);
     157             : 
     158             :     // Build final traction vector
     159       30068 :     computeGlobalTraction(node);
     160             : 
     161             :     // Compute mechanical contact until end of method.
     162       30068 :     const auto penalty_friction = findValue(
     163       30068 :         _dof_to_local_penalty_friction, static_cast<const DofObject *>(node), _penalty_friction);
     164             : 
     165             :     // utilized quantities
     166       30068 :     const auto & normal_pressure = _dof_to_normal_pressure[node];
     167             : 
     168             :     // map the tangential traction and accumulated slip
     169       30068 :     auto & [tangential_traction, old_tangential_traction] = _dof_to_tangential_traction[node];
     170       30068 :     auto & [accumulated_slip, old_accumulated_slip] = _dof_to_accumulated_slip[node];
     171             : 
     172             :     Real normal_lm = -1;
     173       30068 :     if (_dof_to_lagrange_multiplier.find(node) != _dof_to_lagrange_multiplier.end())
     174           0 :       normal_lm = libmesh_map_find(_dof_to_lagrange_multiplier, node);
     175             : 
     176             :     // Keep active set fixed from second Uzawa loop
     177       30068 :     if (normal_lm < -TOLERANCE && normal_pressure > TOLERANCE)
     178             :     {
     179             :       using namespace std;
     180             : 
     181             :       const auto & real_tangential_velocity =
     182        1960 :           libmesh_map_find(_dof_to_real_tangential_velocity, node);
     183        1960 :       const ADTwoVector slip_distance = {real_tangential_velocity[0] * _dt,
     184        1960 :                                          real_tangential_velocity[1] * _dt};
     185             : 
     186             :       // frictional lagrange multiplier (delta lambda^(k)_T)
     187             :       const auto & tangential_lm =
     188        1960 :           _augmented_lagrange_problem ? _dof_to_frictional_lagrange_multipliers[node] : zero;
     189             : 
     190             :       // tangential trial traction (Simo 3.12)
     191             :       // Modified for implementation in MOOSE: Avoid pingponging on frictional sign (max. 0.4
     192             :       // capacity)
     193        1960 :       ADTwoVector inner_iteration_penalty_friction = penalty_friction * slip_distance;
     194             : 
     195             :       const auto slip_metric = std::abs(MetaPhysicL::raw_value(slip_distance).cwiseAbs()(0)) +
     196        1960 :                                std::abs(MetaPhysicL::raw_value(slip_distance).cwiseAbs()(1));
     197             : 
     198        3776 :       if (slip_metric > _epsilon_tolerance &&
     199        5592 :           penalty_friction * slip_distance.norm() >
     200        7408 :               0.4 * _friction_coefficient * std::abs(normal_pressure))
     201             :       {
     202             :         inner_iteration_penalty_friction =
     203           0 :             MetaPhysicL::raw_value(0.4 * _friction_coefficient * std::abs(normal_pressure) /
     204           0 :                                    (penalty_friction * slip_distance.norm())) *
     205           0 :             penalty_friction * slip_distance;
     206             :       }
     207             : 
     208             :       ADTwoVector tangential_trial_traction =
     209        1960 :           old_tangential_traction + tangential_lm + inner_iteration_penalty_friction;
     210             : 
     211             :       // Nonlinearity below
     212        1960 :       ADReal tangential_trial_traction_norm = 0;
     213        3920 :       if (std::abs(tangential_trial_traction(0)) + std::abs(tangential_trial_traction(1)) >
     214             :           _epsilon_tolerance)
     215           0 :         tangential_trial_traction_norm = tangential_trial_traction.norm();
     216             : 
     217        3920 :       ADReal phi_trial = tangential_trial_traction_norm - _friction_coefficient * normal_pressure;
     218             :       tangential_traction = tangential_trial_traction;
     219             : 
     220             :       // Simo considers this a 'return mapping'; we are just capping friction to the Coulomb limit.
     221        1960 :       if (phi_trial > 0.0)
     222             :         // Simo 3.13/3.14 (the penalty formulation has an error in the paper)
     223             :         tangential_traction -=
     224           0 :             phi_trial * tangential_trial_traction / tangential_trial_traction_norm;
     225             : 
     226             :       // track accumulated slip for output purposes
     227        1960 :       accumulated_slip = old_accumulated_slip + MetaPhysicL::raw_value(slip_distance).cwiseAbs();
     228             :     }
     229             :     else
     230             :     {
     231             :       // reset slip and clear traction
     232             :       accumulated_slip.setZero();
     233       28108 :       tangential_traction.setZero();
     234             :     }
     235             : 
     236             :     // Now that we have consistent nodal frictional values, create an interpolated frictional
     237             :     // pressure variable.
     238       30068 :     const auto & test_i = (*_test)[i];
     239       90204 :     for (const auto qp : make_range(_qrule_msm->n_points()))
     240             :     {
     241      120272 :       _frictional_contact_traction_one[qp] += test_i[qp] * tangential_traction(0);
     242       60136 :       _frictional_contact_traction_two[qp] += test_i[qp] * tangential_traction(1);
     243             :     }
     244             :   }
     245       15034 : }
     246             : 
     247             : void
     248        7091 : PenaltySimpleCohesiveZoneModel::finalize()
     249             : {
     250        7091 :   WeightedVelocitiesUserObject::finalize();
     251        7091 :   PenaltyWeightedGapUserObject::selfFinalize();
     252        7091 : }
     253             : 
     254             : const ADVariableValue &
     255           0 : PenaltySimpleCohesiveZoneModel::contactTangentialPressureDirOne() const
     256             : {
     257           0 :   return _frictional_contact_traction_one;
     258             : }
     259             : 
     260             : const ADVariableValue &
     261           0 : PenaltySimpleCohesiveZoneModel::contactTangentialPressureDirTwo() const
     262             : {
     263           0 :   return _frictional_contact_traction_two;
     264             : }