www.mooseframework.org
ContactMaster.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
8 // MOOSE includes
9 #include "ContactMaster.h"
10 #include "NodalArea.h"
11 #include "SystemBase.h"
12 #include "PenetrationInfo.h"
13 #include "MooseMesh.h"
14 
15 #include "libmesh/sparse_matrix.h"
16 #include "libmesh/string_to_enum.h"
17 
18 template <>
19 InputParameters
21 {
22  MooseEnum orders("CONSTANT FIRST SECOND THIRD FOURTH", "FIRST");
23 
24  InputParameters params = validParams<DiracKernel>();
25  params.addRequiredParam<BoundaryName>("boundary", "The master boundary");
26  params.addRequiredParam<BoundaryName>("slave", "The slave boundary");
27  params.addRequiredParam<unsigned int>("component",
28  "An integer corresponding to the direction "
29  "the variable this kernel acts in. (0 for x, "
30  "1 for y, 2 for z)");
31 
32  params.addCoupledVar("disp_x", "The x displacement");
33  params.addCoupledVar("disp_y", "The y displacement");
34  params.addCoupledVar("disp_z", "The z displacement");
35 
36  params.addCoupledVar(
37  "displacements",
38  "The displacements appropriate for the simulation geometry and coordinate system");
39 
40  params.addRequiredCoupledVar("nodal_area", "The nodal area");
41  params.addParam<std::string>("model", "frictionless", "The contact model to use");
42 
43  params.set<bool>("use_displaced_mesh") = true;
44  params.addParam<Real>(
45  "penalty",
46  1e8,
47  "The penalty to apply. This can vary depending on the stiffness of your materials");
48  params.addParam<Real>("friction_coefficient", 0, "The friction coefficient");
49  params.addParam<Real>("tangential_tolerance",
50  "Tangential distance to extend edges of contact surfaces");
51  params.addParam<Real>(
52  "capture_tolerance", 0, "Normal distance from surface within which nodes are captured");
53  params.addParam<Real>(
54  "normal_smoothing_distance",
55  "Distance from edge in parametric coordinates over which to smooth contact normal");
56  params.addParam<std::string>("normal_smoothing_method",
57  "Method to use to smooth normals (edge_based|nodal_normal_based)");
58  params.addParam<MooseEnum>("order", orders, "The finite element order");
59 
60  params.addParam<Real>("tension_release",
61  0.0,
62  "Tension release threshold. A node in contact "
63  "will not be released if its tensile load is below "
64  "this value. No tension release if negative.");
65 
66  params.addParam<std::string>("formulation", "default", "The contact formulation");
67  params.addParam<bool>(
68  "normalize_penalty",
69  false,
70  "Whether to normalize the penalty parameter with the nodal area for penalty contact.");
71  return params;
72 }
73 
74 ContactMaster::ContactMaster(const InputParameters & parameters)
75  : DiracKernel(parameters),
76  _component(getParam<unsigned int>("component")),
77  _model(contactModel(getParam<std::string>("model"))),
78  _formulation(contactFormulation(getParam<std::string>("formulation"))),
79  _normalize_penalty(getParam<bool>("normalize_penalty")),
80  _penetration_locator(
81  getPenetrationLocator(getParam<BoundaryName>("boundary"),
82  getParam<BoundaryName>("slave"),
83  Utility::string_to_enum<Order>(getParam<MooseEnum>("order")))),
84  _penalty(getParam<Real>("penalty")),
85  _friction_coefficient(getParam<Real>("friction_coefficient")),
86  _tension_release(getParam<Real>("tension_release")),
87  _capture_tolerance(getParam<Real>("capture_tolerance")),
88  _updateContactSet(true),
89  _residual_copy(_sys.residualGhosted()),
90  _mesh_dimension(_mesh.dimension()),
91  _vars(3, libMesh::invalid_uint),
92  _nodal_area_var(getVar("nodal_area", 0)),
93  _aux_system(_nodal_area_var->sys()),
94  _aux_solution(_aux_system.currentSolution())
95 {
96  if (isParamValid("displacements"))
97  {
98  // modern parameter scheme for displacements
99  for (unsigned int i = 0; i < coupledComponents("displacements"); ++i)
100  _vars[i] = coupled("displacements", i);
101  }
102  else
103  {
104  // Legacy parameter scheme for displacements
105  if (isParamValid("disp_x"))
106  _vars[0] = coupled("disp_x");
107  if (isParamValid("disp_y"))
108  _vars[1] = coupled("disp_y");
109  if (isParamValid("disp_z"))
110  _vars[2] = coupled("disp_z");
111 
112  mooseDeprecated("use the `displacements` parameter rather than the `disp_*` parameters (those "
113  "will go away with the deprecation of the Solid Mechanics module).");
114  }
115 
116  if (parameters.isParamValid("tangential_tolerance"))
117  _penetration_locator.setTangentialTolerance(getParam<Real>("tangential_tolerance"));
118 
119  if (parameters.isParamValid("normal_smoothing_distance"))
120  _penetration_locator.setNormalSmoothingDistance(getParam<Real>("normal_smoothing_distance"));
121 
122  if (parameters.isParamValid("normal_smoothing_method"))
123  _penetration_locator.setNormalSmoothingMethod(
124  parameters.get<std::string>("normal_smoothing_method"));
125 
127  _penetration_locator.setUpdate(false);
128 
129  if (_friction_coefficient < 0)
130  mooseError("The friction coefficient must be nonnegative");
131 }
132 
133 void
135 {
136  if (_component == 0)
137  {
138  if (_updateContactSet)
140 
141  _updateContactSet = true;
142  }
143 }
144 
145 void
147 {
148  if (_component == 0)
149  {
150  updateContactSet(true);
151  _updateContactSet = false;
152  }
153 }
154 
155 void
156 ContactMaster::updateContactSet(bool beginning_of_step)
157 {
158  std::map<dof_id_type, PenetrationInfo *>::iterator
159  it = _penetration_locator._penetration_info.begin(),
160  end = _penetration_locator._penetration_info.end();
161  for (; it != end; ++it)
162  {
163  const dof_id_type slave_node_num = it->first;
164  PenetrationInfo * pinfo = it->second;
165 
166  // Skip this pinfo if there are no DOFs on this node.
167  if (!pinfo || pinfo->_node->n_comp(_sys.number(), _vars[_component]) < 1)
168  continue;
169 
170  if (beginning_of_step)
171  {
172  pinfo->_locked_this_step = 0;
173  pinfo->_starting_elem = it->second->_elem;
174  pinfo->_starting_side_num = it->second->_side_num;
175  pinfo->_starting_closest_point_ref = it->second->_closest_point_ref;
176  pinfo->_contact_force_old = pinfo->_contact_force;
177  pinfo->_accumulated_slip_old = pinfo->_accumulated_slip;
178  pinfo->_frictional_energy_old = pinfo->_frictional_energy;
179  }
180 
181  const Real contact_pressure = -(pinfo->_normal * pinfo->_contact_force) / nodalArea(*pinfo);
182  const Real distance = pinfo->_normal * (pinfo->_closest_point - _mesh.nodeRef(slave_node_num));
183 
184  // Capture
185  if (!pinfo->isCaptured() &&
186  MooseUtils::absoluteFuzzyGreaterEqual(distance, 0, _capture_tolerance))
187  {
188  pinfo->capture();
189 
190  // Increment the lock count every time the node comes back into contact from not being in
191  // contact.
192  if (_formulation == CF_KINEMATIC)
193  ++pinfo->_locked_this_step;
194  }
195  // Release
196  else if (_model != CM_GLUED && pinfo->isCaptured() && _tension_release >= 0 &&
197  -contact_pressure >= _tension_release && pinfo->_locked_this_step < 2)
198  {
199  pinfo->release();
200  pinfo->_contact_force.zero();
201  }
202 
203  if (_formulation == CF_AUGMENTED_LAGRANGE && pinfo->isCaptured())
204  pinfo->_lagrange_multiplier -= getPenalty(*pinfo) * distance;
205  }
206 }
207 
208 void
210 {
211  _point_to_info.clear();
212 
213  std::map<dof_id_type, PenetrationInfo *>::iterator
214  it = _penetration_locator._penetration_info.begin(),
215  end = _penetration_locator._penetration_info.end();
216 
217  for (; it != end; ++it)
218  {
219  PenetrationInfo * pinfo = it->second;
220 
221  // Skip this pinfo if there are no DOFs on this node.
222  if (!pinfo || pinfo->_node->n_comp(_sys.number(), _vars[_component]) < 1)
223  continue;
224 
225  if (pinfo->isCaptured())
226  {
227  addPoint(pinfo->_elem, pinfo->_closest_point);
228  _point_to_info[pinfo->_closest_point] = pinfo;
229  computeContactForce(pinfo);
230  }
231  }
232 }
233 
234 void
235 ContactMaster::computeContactForce(PenetrationInfo * pinfo)
236 {
237  const Node * node = pinfo->_node;
238 
239  RealVectorValue res_vec;
240  // Build up residual vector
241  for (unsigned int i = 0; i < _mesh_dimension; ++i)
242  {
243  dof_id_type dof_number = node->dof_number(0, _vars[i], 0);
244  res_vec(i) = _residual_copy(dof_number);
245  }
246 
247  const Real area = nodalArea(*pinfo);
248 
249  RealVectorValue distance_vec(_mesh.nodeRef(node->id()) - pinfo->_closest_point);
250  RealVectorValue pen_force(_penalty * distance_vec);
251  if (_normalize_penalty)
252  pen_force *= area;
253 
254  if (_model == CM_FRICTIONLESS)
255  {
256  switch (_formulation)
257  {
258  case CF_DEFAULT:
259  pinfo->_contact_force = -pinfo->_normal * (pinfo->_normal * res_vec);
260  break;
261  case CF_PENALTY:
262  pinfo->_contact_force = pinfo->_normal * (pinfo->_normal * pen_force);
263  break;
265  pinfo->_contact_force =
266  (pinfo->_normal *
267  (pinfo->_normal * (pen_force + pinfo->_lagrange_multiplier * pinfo->_normal)));
268  break;
269  default:
270  mooseError("Invalid contact formulation");
271  break;
272  }
273  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
274  }
275  else if (_model == CM_COULOMB && _formulation == CF_PENALTY)
276  {
277  distance_vec =
278  pinfo->_incremental_slip +
279  (pinfo->_normal * (_mesh.nodeRef(node->id()) - pinfo->_closest_point)) * pinfo->_normal;
280  pen_force = _penalty * distance_vec;
281  if (_normalize_penalty)
282  pen_force *= area;
283 
284  // Frictional capacity
285  // const Real capacity( _friction_coefficient * (pen_force * pinfo->_normal < 0 ? -pen_force *
286  // pinfo->_normal : 0) );
287  const Real capacity(_friction_coefficient *
288  (res_vec * pinfo->_normal > 0 ? res_vec * pinfo->_normal : 0));
289 
290  // Elastic predictor
291  pinfo->_contact_force =
292  pen_force +
293  (pinfo->_contact_force_old - pinfo->_normal * (pinfo->_normal * pinfo->_contact_force_old));
294  RealVectorValue contact_force_normal((pinfo->_contact_force * pinfo->_normal) * pinfo->_normal);
295  RealVectorValue contact_force_tangential(pinfo->_contact_force - contact_force_normal);
296 
297  // Tangential magnitude of elastic predictor
298  const Real tan_mag(contact_force_tangential.norm());
299 
300  if (tan_mag > capacity)
301  {
302  pinfo->_contact_force = contact_force_normal + capacity * contact_force_tangential / tan_mag;
303  pinfo->_mech_status = PenetrationInfo::MS_SLIPPING;
304  }
305  else
306  {
307  pinfo->_mech_status = PenetrationInfo::MS_STICKING;
308  }
309  }
310  else if (_model == CM_GLUED || (_model == CM_COULOMB && _formulation == CF_DEFAULT))
311  {
312  switch (_formulation)
313  {
314  case CF_DEFAULT:
315  pinfo->_contact_force = -res_vec;
316  break;
317  case CF_PENALTY:
318  pinfo->_contact_force = pen_force;
319  break;
321  pinfo->_contact_force =
322  pen_force + pinfo->_lagrange_multiplier * distance_vec / distance_vec.norm();
323  break;
324  default:
325  mooseError("Invalid contact formulation");
326  break;
327  }
328  pinfo->_mech_status = PenetrationInfo::MS_STICKING;
329  }
330  else
331  {
332  mooseError("Invalid or unavailable contact model");
333  }
334 }
335 
336 Real
338 {
339  PenetrationInfo * pinfo = _point_to_info[_current_point];
340  Real resid = -pinfo->_contact_force(_component);
341  return _test[_i][_qp] * resid;
342 }
343 
344 Real
346 {
347 
348  PenetrationInfo * pinfo = _point_to_info[_current_point];
349  Real penalty = _penalty;
350  if (_normalize_penalty)
351  penalty *= nodalArea(*pinfo);
352 
353  switch (_model)
354  {
355  case CM_FRICTIONLESS:
356  switch (_formulation)
357  {
358  case CF_DEFAULT:
359  return 0;
360  break;
361  case CF_PENALTY:
363  return _test[_i][_qp] * penalty * _phi[_j][_qp] * pinfo->_normal(_component) *
364  pinfo->_normal(_component);
365  break;
366  default:
367  mooseError("Invalid contact formulation");
368  break;
369  }
370  break;
371  case CM_GLUED:
372  case CM_COULOMB:
373  switch (_formulation)
374  {
375  case CF_DEFAULT:
376  return 0;
377  break;
378  case CF_PENALTY:
380  return _test[_i][_qp] * penalty * _phi[_j][_qp];
381  break;
382  default:
383  mooseError("Invalid contact formulation");
384  break;
385  }
386  break;
387  default:
388  mooseError("Invalid or unavailable contact model");
389  break;
390  }
391 
392  return 0;
393  /*
394  if (_i != _j)
395  return 0;
396 
397  Node * node = pinfo->_node;
398 
399  RealVectorValue jac_vec;
400 
401  // Build up jac vector
402  for (unsigned int i=0; i<_dim; i++)
403  {
404  long int dof_number = node->dof_number(0, _vars[i], 0);
405  jac_vec(i) = _jacobian_copy(dof_number, dof_number);
406  }
407 
408  Real jac_mag = pinfo->_normal * jac_vec;
409 
410  return _test[_i][_qp]*pinfo->_normal(_component)*jac_mag;
411  */
412 }
413 
416 {
417  ContactModel model(CM_INVALID);
418  std::transform(name.begin(), name.end(), name.begin(), ::tolower);
419 
420  if ("frictionless" == name)
421  model = CM_FRICTIONLESS;
422  else if ("glued" == name)
423  model = CM_GLUED;
424  else if ("coulomb" == name)
425  model = CM_COULOMB;
426  else
427  ::mooseError("Invalid contact model found: ", name);
428 
429  return model;
430 }
431 
434 {
435  ContactFormulation formulation(CF_INVALID);
436  std::transform(name.begin(), name.end(), name.begin(), ::tolower);
437 
438  if ("default" == name || "kinematic" == name)
439  formulation = CF_DEFAULT;
440 
441  else if ("penalty" == name)
442  formulation = CF_PENALTY;
443 
444  else if ("augmented_lagrange" == name)
445  formulation = CF_AUGMENTED_LAGRANGE;
446 
447  else if ("tangential_penalty" == name)
448  formulation = CF_TANGENTIAL_PENALTY;
449 
450  if (formulation == CF_INVALID)
451  ::mooseError("Invalid formulation found: ", name);
452 
453  return formulation;
454 }
455 
456 Real
457 ContactMaster::nodalArea(PenetrationInfo & pinfo)
458 {
459  const Node * node = pinfo._node;
460 
461  dof_id_type dof = node->dof_number(_aux_system.number(), _nodal_area_var->number(), 0);
462 
463  Real area = (*_aux_solution)(dof);
464  if (area == 0.0)
465  {
466  if (_t_step > 1)
467  mooseError("Zero nodal area found");
468  else
469  area = 1.0; // Avoid divide by zero during initialization
470  }
471  return area;
472 }
473 
474 Real
475 ContactMaster::getPenalty(PenetrationInfo & pinfo)
476 {
477  Real penalty = _penalty;
478  if (_normalize_penalty)
479  penalty *= nodalArea(pinfo);
480 
481  return penalty;
482 }
MooseVariable * _nodal_area_var
Definition: ContactMaster.h:74
const Real _tension_release
Definition: ContactMaster.h:62
virtual void addPoints()
ContactModel
Definition: ContactMaster.h:14
Real getPenalty(PenetrationInfo &pinfo)
SystemBase & _aux_system
Definition: ContactMaster.h:75
const Real _friction_coefficient
Definition: ContactMaster.h:61
ContactFormulation
Definition: ContactMaster.h:22
ContactFormulation contactFormulation(const std::string &the_name)
virtual void jacobianSetup()
virtual void timestepSetup()
const ContactFormulation _formulation
Definition: ContactMaster.h:56
void computeContactForce(PenetrationInfo *pinfo)
const unsigned int _mesh_dimension
Definition: ContactMaster.h:70
const unsigned int _component
Definition: ContactMaster.h:54
NumericVector< Number > & _residual_copy
Definition: ContactMaster.h:66
ContactMaster(const InputParameters &parameters)
Definition: ContactMaster.C:74
const Real _penalty
Definition: ContactMaster.h:60
virtual void updateContactSet(bool beginning_of_step=false)
bool _updateContactSet
Definition: ContactMaster.h:64
const Real _capture_tolerance
Definition: ContactMaster.h:63
std::map< Point, PenetrationInfo * > _point_to_info
Definition: ContactMaster.h:68
const bool _normalize_penalty
Definition: ContactMaster.h:57
PenetrationLocator & _penetration_locator
Definition: ContactMaster.h:58
static ContactFormulation contactFormulation(std::string name)
virtual Real computeQpJacobian()
ContactModel contactModel(const std::string &the_name)
static ContactModel contactModel(std::string name)
virtual Real computeQpResidual()
const ContactModel _model
Definition: ContactMaster.h:55
std::vector< unsigned int > _vars
Definition: ContactMaster.h:72
Real nodalArea(PenetrationInfo &pinfo)
InputParameters validParams< ContactMaster >()
Definition: ContactMaster.C:20