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 "ContactSlipDamper.h"
      11             : #include "FEProblem.h"
      12             : #include "DisplacedProblem.h"
      13             : #include "AuxiliarySystem.h"
      14             : #include "PenetrationLocator.h"
      15             : #include "NearestNodeLocator.h"
      16             : 
      17             : registerMooseObject("ContactApp", ContactSlipDamper);
      18             : 
      19             : InputParameters
      20         234 : ContactSlipDamper::validParams()
      21             : {
      22         234 :   InputParameters params = GeneralDamper::validParams();
      23         468 :   params.addParam<std::vector<BoundaryName>>(
      24             :       "primary", "IDs of the primary surfaces for which slip reversals should be damped");
      25         468 :   params.addParam<std::vector<BoundaryName>>(
      26             :       "secondary", "IDs of the secondary surfaces for which slip reversals should be damped");
      27         468 :   params.addParam<Real>(
      28         468 :       "max_iterative_slip", std::numeric_limits<Real>::max(), "Maximum iterative slip");
      29         702 :   params.addRangeCheckedParam<Real>("min_damping_factor",
      30         468 :                                     0.0,
      31             :                                     "min_damping_factor < 1.0",
      32             :                                     "Minimum permissible value for damping factor");
      33         468 :   params.addParam<Real>("damping_threshold_factor",
      34         468 :                         1.0e3,
      35             :                         "If previous iterations's slip is below the slip tolerance, "
      36             :                         "only damp a slip reversal if the slip magnitude is greater "
      37             :                         "than than this factor times the old slip.");
      38         468 :   params.addParam<bool>("debug_output", false, "Output detailed debugging information");
      39         234 :   params.addClassDescription(
      40             :       "Damp the iterative solution to minimize oscillations in frictional contact "
      41             :       "constriants between nonlinear iterations");
      42         234 :   return params;
      43           0 : }
      44             : 
      45         117 : ContactSlipDamper::ContactSlipDamper(const InputParameters & parameters)
      46             :   : GeneralDamper(parameters),
      47         117 :     _aux_sys(parameters.get<FEProblemBase *>("_fe_problem_base")->getAuxiliarySystem()),
      48         117 :     _displaced_problem(parameters.get<FEProblemBase *>("_fe_problem_base")->getDisplacedProblem()),
      49         117 :     _num_contact_nodes(0),
      50         117 :     _num_sticking(0),
      51         117 :     _num_slipping(0),
      52         117 :     _num_slipping_friction(0),
      53         117 :     _num_stick_locked(0),
      54         117 :     _num_slip_reversed(0),
      55         117 :     _max_iterative_slip(parameters.get<Real>("max_iterative_slip")),
      56         117 :     _min_damping_factor(parameters.get<Real>("min_damping_factor")),
      57         117 :     _damping_threshold_factor(parameters.get<Real>("damping_threshold_factor")),
      58         351 :     _debug_output(parameters.get<bool>("debug_output"))
      59             : {
      60         117 :   if (!_displaced_problem)
      61           0 :     mooseError("Must have displaced problem to use ContactSlipDamper");
      62             : 
      63         234 :   std::vector<BoundaryName> primary = getParam<std::vector<BoundaryName>>("primary");
      64         351 :   std::vector<BoundaryName> secondary = getParam<std::vector<BoundaryName>>("secondary");
      65             : 
      66         117 :   unsigned int num_interactions = primary.size();
      67         117 :   if (num_interactions != secondary.size())
      68           0 :     mooseError(
      69             :         "Sizes of primary surface and secondary surface lists must match in ContactSlipDamper");
      70         117 :   if (num_interactions == 0)
      71           0 :     mooseError("Must define at least one primary/secondary pair in ContactSlipDamper");
      72             : 
      73         234 :   for (unsigned int i = 0; i < primary.size(); ++i)
      74             :   {
      75         117 :     std::pair<int, int> ms_pair(_subproblem.mesh().getBoundaryID(primary[i]),
      76         117 :                                 _subproblem.mesh().getBoundaryID(secondary[i]));
      77         117 :     _interactions.insert(ms_pair);
      78             :   }
      79         117 : }
      80             : 
      81             : void
      82         540 : ContactSlipDamper::timestepSetup()
      83             : {
      84         540 :   GeometricSearchData & displaced_geom_search_data = _displaced_problem->geomSearchData();
      85             :   const auto & penetration_locators = displaced_geom_search_data._penetration_locators;
      86             : 
      87        1080 :   for (const auto & pl : penetration_locators)
      88             :   {
      89         540 :     PenetrationLocator & pen_loc = *pl.second;
      90             : 
      91         540 :     if (operateOnThisInteraction(pen_loc))
      92             :     {
      93        2940 :       for (const auto & secondary_node_num : pen_loc._nearest_node._secondary_nodes)
      94        2400 :         if (pen_loc._penetration_info[secondary_node_num])
      95             :         {
      96        2360 :           PenetrationInfo & info = *pen_loc._penetration_info[secondary_node_num];
      97        2360 :           const Node & node = _displaced_problem->mesh().nodeRef(secondary_node_num);
      98             : 
      99        2360 :           if (node.processor_id() == processor_id())
     100             :             //              && info.isCaptured())                  //TODO maybe just set this
     101             :             //              everywhere?
     102        1660 :             info._slip_reversed = false;
     103             :         }
     104             :     }
     105             :   }
     106         540 : }
     107             : 
     108             : Real
     109        3547 : ContactSlipDamper::computeDamping(const NumericVector<Number> & solution,
     110             :                                   const NumericVector<Number> & /*update*/)
     111             : {
     112             :   std::map<unsigned int, const NumericVector<Number> *> nl_soln;
     113        3547 :   nl_soln.emplace(_sys.number(), &solution);
     114             : 
     115             :   // Do new contact search to update positions of slipped nodes
     116        3547 :   _displaced_problem->updateMesh(nl_soln, *_aux_sys.currentSolution());
     117             : 
     118        3547 :   Real damping = 1.0;
     119             : 
     120        3547 :   _num_contact_nodes = 0;
     121        3547 :   _num_sticking = 0;
     122        3547 :   _num_slipping = 0;
     123        3547 :   _num_slipping_friction = 0;
     124        3547 :   _num_stick_locked = 0;
     125        3547 :   _num_slip_reversed = 0;
     126             : 
     127        3547 :   GeometricSearchData & displaced_geom_search_data = _displaced_problem->geomSearchData();
     128             :   const auto & penetration_locators = displaced_geom_search_data._penetration_locators;
     129             : 
     130        7094 :   for (const auto & pl : penetration_locators)
     131             :   {
     132        3547 :     PenetrationLocator & pen_loc = *pl.second;
     133             : 
     134        3547 :     if (operateOnThisInteraction(pen_loc))
     135             :     {
     136       20162 :       for (const auto & secondary_node_num : pen_loc._nearest_node._secondary_nodes)
     137       16615 :         if (pen_loc._penetration_info[secondary_node_num])
     138             :         {
     139       16395 :           PenetrationInfo & info = *pen_loc._penetration_info[secondary_node_num];
     140       16395 :           const Node & node = _displaced_problem->mesh().nodeRef(secondary_node_num);
     141             : 
     142       16395 :           if (node.processor_id() == processor_id())
     143             :           {
     144       10590 :             if (info.isCaptured())
     145             :             {
     146        6780 :               _num_contact_nodes++;
     147        6780 :               if (info._mech_status == PenetrationInfo::MS_STICKING)
     148         324 :                 _num_sticking++;
     149        6456 :               else if (info._mech_status == PenetrationInfo::MS_SLIPPING)
     150          66 :                 _num_slipping++;
     151        6390 :               else if (info._mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
     152        6390 :                 _num_slipping_friction++;
     153        6780 :               if (info._stick_locked_this_step >= 2) // TODO get from contact interaction
     154         334 :                 _num_stick_locked++;
     155             : 
     156             :               RealVectorValue tangential_inc_slip_prev_iter =
     157             :                   info._incremental_slip_prev_iter -
     158             :                   (info._incremental_slip_prev_iter * info._normal) * info._normal;
     159             :               RealVectorValue tangential_inc_slip =
     160             :                   info._incremental_slip - (info._incremental_slip * info._normal) * info._normal;
     161             : 
     162             :               RealVectorValue tangential_it_slip =
     163             :                   tangential_inc_slip - tangential_inc_slip_prev_iter;
     164        6780 :               Real node_damping_factor = 1.0;
     165        6780 :               if ((tangential_inc_slip_prev_iter * tangential_inc_slip < 0.0) &&
     166             :                   info._mech_status == PenetrationInfo::MS_SLIPPING_FRICTION)
     167             :               {
     168         132 :                 info._slip_reversed = true;
     169         132 :                 _num_slip_reversed++;
     170         132 :                 Real prev_iter_slip_mag = tangential_inc_slip_prev_iter.norm();
     171             :                 RealVectorValue prev_iter_slip_dir =
     172             :                     tangential_inc_slip_prev_iter / prev_iter_slip_mag;
     173             :                 Real cur_it_slip_in_old_dir = tangential_it_slip * prev_iter_slip_dir;
     174             : 
     175         132 :                 if (prev_iter_slip_mag > info._slip_tol ||
     176          16 :                     cur_it_slip_in_old_dir > -_damping_threshold_factor * prev_iter_slip_mag)
     177         116 :                   node_damping_factor =
     178         116 :                       1.0 - (cur_it_slip_in_old_dir + prev_iter_slip_mag) / cur_it_slip_in_old_dir;
     179             : 
     180         132 :                 if (node_damping_factor < 0.0)
     181           0 :                   mooseError("Damping factor can't be negative");
     182             : 
     183         132 :                 if (node_damping_factor < _min_damping_factor)
     184           0 :                   node_damping_factor = _min_damping_factor;
     185             :               }
     186             : 
     187        6780 :               if (tangential_it_slip.norm() > _max_iterative_slip)
     188           0 :                 node_damping_factor =
     189           0 :                     (tangential_it_slip.norm() - _max_iterative_slip) / tangential_it_slip.norm();
     190             : 
     191        6780 :               if (_debug_output && node_damping_factor < 1.0)
     192           0 :                 _console << "Damping node: " << node.id()
     193           0 :                          << " prev iter slip: " << info._incremental_slip_prev_iter
     194           0 :                          << " curr iter slip: " << info._incremental_slip
     195           0 :                          << " slip_tol: " << info._slip_tol
     196           0 :                          << " damping factor: " << node_damping_factor << std::endl;
     197             : 
     198        6780 :               if (node_damping_factor < damping)
     199         116 :                 damping = node_damping_factor;
     200             :             }
     201             :           }
     202             :         }
     203             :     }
     204             :   }
     205        3547 :   _console << std::flush;
     206        3547 :   _communicator.sum(_num_contact_nodes);
     207        3547 :   _communicator.sum(_num_sticking);
     208        3547 :   _communicator.sum(_num_slipping);
     209        3547 :   _communicator.sum(_num_slipping_friction);
     210        3547 :   _communicator.sum(_num_stick_locked);
     211        3547 :   _communicator.sum(_num_slip_reversed);
     212        3547 :   _communicator.min(damping);
     213             : 
     214             :   _console << "   ContactSlipDamper: Damping     #Cont    #Stick     #Slip #SlipFric #StickLock  "
     215        3547 :               "#SlipRev\n";
     216             : 
     217        3547 :   _console << std::right << std::setw(29) << damping << std::setw(10) << _num_contact_nodes
     218        3547 :            << std::setw(10) << _num_sticking << std::setw(10) << _num_slipping << std::setw(10)
     219        3547 :            << _num_slipping_friction << std::setw(11) << _num_stick_locked << std::setw(10)
     220        3547 :            << _num_slip_reversed << "\n\n";
     221        3547 :   _console << std::flush;
     222             : 
     223        3547 :   return damping;
     224             : }
     225             : 
     226             : bool
     227        4087 : ContactSlipDamper::operateOnThisInteraction(const PenetrationLocator & pen_loc)
     228             : {
     229             :   bool operate_on_this_interaction = false;
     230             :   std::set<std::pair<int, int>>::iterator ipit;
     231             :   std::pair<int, int> ms_pair(pen_loc._primary_boundary, pen_loc._secondary_boundary);
     232             :   ipit = _interactions.find(ms_pair);
     233        4087 :   if (ipit != _interactions.end())
     234             :     operate_on_this_interaction = true;
     235        4087 :   return operate_on_this_interaction;
     236             : }