`***continuum_contact`#

Description#

***continuum_contact is an alternative to ***contact, and is based on a penalty method to enforce the contact constraints. Compared to ***contact that solves a local problem restricted to nodes/degrees of freedom entering contact, in this formulation, contact elements are dynamically added to the global FE equilibrium system, depending on the contact state at the current Newton iteration.

An important implication of this is that since all non-linearities are treated at the same time (i.e. contact, but also material or geometric non-linearities), a larger number of equilibrium iterations are generally needed by this method compared to the previous one. This may become a major cause of divergence in case of so-called contact fluttering (new contact elements entering or leaving the non-linear system at each equilibrium iteration). A special procedure has been devised to alleviate this particular problem, and can be triggered by the **stabilize_behavior command.

On the other hand, large CPU times (mainly involved in the calculation of the flexibility matrix, see page ) may be needed to build the local contact problem with the ***contact algorithm. As a result, the present formulation is in general more efficient, as soon as the number of nodes entering contact becomes large. Also, another major advantage of ***continuum_contact, is that it’s the only contact formulation that supports Z-set’s domain-decomposition parallel solvers, thus extending the field of HPC to large contact problems.

Contact elements and definition of the warning zone#

Definition of contacting parts in the FE mesh is similar to the one used in the ***contact. An arbitrary number of contact zones can be added, each zone being characterized by:

a slave (impactor) nset,
a master (target) bset,
a warning distance, defining a critical value of the distance to the master surface below which contact constraints should be added for nodes belonging to the slave nset.

As shown in the next figure, when the altitude of a given slave node is below the warning distance, a contact element is automatically created that links this slave node to nodes of the closest bset element in the master surface definition.

The contact penalty formulation#

In the penalty method, contact constraints are only enforced approximately and some degree of penetration is allowed, depending on the value used for the $p_n$ penalty parameter. Large values decrease the amount of penetration allowable, but can lead to convergence problems, such that choosing an optimal value for $p_n$ is a tradeoff between precision/efficiency requirements. An automatic choice of this parameter, based on the stiffness of the underlying master surface, can be activated by means of the *automatic_penalty command.

Current limitations and use with quadratic elements#

in the current release (8.6), ***continuum_contact is limited to frictionless contact only.
specific issues related to the handling of contact at middle nodes of quadratic elements are not treated explicitly in the formulation. The only solution currently available for quadratic meshes is to use the *add_mpc option, that automatically adds relationships at middle nodes, to tie them linearly to corner nodes of the corresponding edge of second-order elements.
***continuum_contact is not yet supported for dynamic problems.

Stabilization of the material behavior#

Command *stabilize_behavior greatly improves convergence in the general case where both material and contact non-linearities are included. When defined, convergence of the equilibrium problem for each increment, is split-up in 2 successive phases:

severe non-linear iterations: during this initial phase, material non-linearities are de-activated (i.e. the material behavior is linearized around the state obtained at the end of the last converged increment), and only contact non-linearities are active.
standard equilibrium iterations: once a preliminary stabilized contact state has been obtained at the end of the previous phase, material non-linearities are released, and conventional equilibrium iterations proceed until convergence. Note that the contact non-linearities are included, such that the contact state may indeed change during this phase. However, smaller contact-state variations are expected compared to the case where the first phase is not done.

Experience has shown that convergence problems associated with large contact-state fluctuations interacting in an harmful way with a strongly non-linear material response in the contact zone, are alleviated, and even completely removed, by the use of this *stabilize_behavior procedure.

Use with parallel domain-decomposition solvers#

This contact formulation is now fully supported by the generic_dd (or feti) parallel solvers, meaning that the occurrence of contact between sub-domains is automatically accounted for. In this case, new elements are added to the sub-domain containing the master surface, while the slave nodes of the impactor sub-domains are automatically added to the list of nodes of the master domain, and to the corresponding interfaces.

For example, in the case of the figure hereafter, where contact is defined between slave nodes of Domain 1 and a master surface included in Domain 2, steps involved in the parallel implementation are the following one:

Parallel contact detection: slave node $S$ of Domain 1 is the warning zone of the master surface of Domain 2 (distance $PS < w_d$, where P is the orthogonal projection of $S$ on the master surface, and $w_d$ is the warning distance).
Update of the master subdomain definition: Slave node $S$ is added to list of nodes of Domain 2, while a new contact element $(S,M_1,M_2)$ is also automatically created in this domain.
Management of the interfaces between domains: Node $S$ is now shared by both domains, and corresponding degrees of freedoms are added to the interface definition.

../../_images/parallel_continuum_contact.svg

Some specific issues related to the use of domain-decomposition solvers then need some further description:

as shown in the figure above, if a node ($S$) of the slave nset is within the warning zone of the master bset defined in another domain, new nodes/elements are added to the master domain (Domain 1 in the figure). In the event that no contact at node $S$ is detected during equilibrium iterations (positive altitude, no penetration of the master surface), the stiffness associated to this contact element is zero, such that $S$ is not connected to other elements in Domain 2, giving rise to null pivots during tangent matrix factorization. To avoid this problem, as shown in the figure below, a small stiffness is automatically added, the value of which is defined by the Parallel.Contact.Stabilize global parameter (default value is 1.e-3). The stiffness of those stabilization springs can be redefined by using the ***global_parameter command, or directly at the command line level using the -s switch:
```
$ Zrun -s Parallel.Contact.Stabilize 1.e-1 ...
```

../../_images/parallel_continuum_contact2.svg

the generic_dd solver supports many options (choice of primal/dual treatment for user-specified interface degrees of freedom), and preconditioning methods for the iterative resolution of the interface problem. Particular options may or may not be well adapted to the handling of those contact problems, in particular because of a poor conditioning of the global stiffness due to the penalty method. Experience has shown that a full dual treatment of the interface problem associated to the use of a Dirichlet preconditioner is the most robust combination. Typical input commands of generic_dd are included hereafter for reference:
```
***linear_solver generic_dd
 **local_solver mumps *no_option
 **dof_kind
  *dof all dual
 **iterative_solver nncg
   precision 1.e-06
   max_iteration 1000
 **precond full
 **scaling stiffness
 **projector_schur none
 **projector_scaling stiffness
```

Syntax#

***continuum_contact $~\,$ **zone penalty $~\,~\,$ *master mname $~\,~\,$ *slave sname $~\,~\,$ *warning_distance dist $~\,$ [ *penalty_normal pn ] $~\,$ [ *automatic_penalty ] $~\,$ [ *auto_factor fact ] $~\,$ [ *keep_master ] $~\,$ [ *add_mpc [ middle | corner ] [ full | everywhere ] ] $~\,$ [ *stabilize [ everything ] stab ] $~\,$ [ *damp ] $~\,$ [ *remove ctip ] $~\,$ [ **zone penalty $~\,~\,~\,$ … ] [ **stabilize_behavior ]

**zone

this command starts the definition of a contact zone. Note that as many zones as needed may be added to describe contact between various sub-components in the global mesh.

*master

defines the master target bset. mname is a liset for 2D problems, or a faset in 3D. Care must be taken that the normal to this bset is pointing outwards the interior of the associated body.

*slave

defines the name sname of an nset containing impactor slave nodes for this particular zone.

*warning_distance

defines the critical distance between a slave node and the closest point on the master surface, below which a contact element is created for this particular slave. Note that if the slave node is not found in contact during equilibrium iterations (positive altitude, no penetration), null stiffness/nodal forces will be calculated for the corresponding contact element, such that setting a too large value for dist does not change the results. However, the size of the global stiffness matrix is increased by those useless elements, and for efficiency reasons dist should not be taken too big.

*penalty_normal

this optional command allows to define the value of the penalty parameter pn used to enforce non-penetration of slave nodes inside the master surface. As discussed above, setting manually an optimal value to pn may be cumbersome, and using the *automatic_penalty command is a practical alternative, especially when the slave surface is not homogeneous, such that using the exact same value for all the contact zone is not suitable.

*automatic_penalty

with this option a value of the penalty coefficient is automatically computed for each slave node, based on the stiffness of the underlying master surface at this particular node. *auto_factor this optional command allows to modify the heuristic used in the previous *automatic_penalty method, where pn is defined by the following equation:

(58)#\[pn~=~fact~.~stiffness\_bset\]

with $stiffness\_bset$ is the average of the stiffness at nodes of the closest bset element. Experience has shown that the default value of $fact=0.1$ is a good choice, such that using of *auto_factor to redefine fact is scarcely needed.

*keep_master

in the default mode, all bset elements within the warning zone of a given slave node are selected as contacting candidates at the beginning of the increment. Then, during equilibrium iterations, a contact element is created only with the closest one in this list of candidate, a choice that may change from one iteration to the other. If command *keep_master is defined, the closest element at the beginning of the increment is retained during all iterations. This command may be used to speed-up convergence in case of small-sliding.

*add_mpc

this command should be used for contact between quadratic FE meshes. In this case, contact conditions are only enforced for a subset of the slave nodes of quadratic elements, i.e. corner nodes when corner is defined (default) or alternatively middle nodes when middle is used. In all cases, linear relationships are automatically added to middle nodes in order that they stay aligned with corner nodes of the corresponding element edge. Optional subcommand all adds those MPCs at all the middle nodes of the slave nset regardless of the contact state. If everywhere is used instead of all, this procedure is extended to all middle nodes of the master surface.

*remove

this command typically takes as argument the name ctip of an nset containing nodes at the crack tip, and may for example be used in Zcracks simulations when contact between the crack lips is defined. In this case one side of the crack is defined as master and the other one as slave, while the crack tip is indeed common to both sides, such that enforcing contact conditions for those nodes has no sense an could inhibit convergence.

*stabilize

this command is intended to remove rigid-body movements associated with the definition of a contacting component whose motion is not completely prescribed by boundary conditions, and is only restrained by contact with other components. In this case springs with a stiffness corresponding to the stab argument are added:

at the slave node for each contact element (grounding of the slave node)
between the slave and the master nodes of each contact element.

When the option everything is added, not only the slave node, but also the master nodes are grounded. The stiffness damp should be chosen significantly smaller than the stiffness of the contacting bodies, but large enough to remove rigid motions. Since adding those stabilization springs, only change the global stiffness matrix, setting a too large value for damp can degrade convergence.

*damp

when this option is used in conjunction with *stabilize, corresponding internal forces are also added to the problem residual. This strategy improves convergence, but the quality of the solution is now altered for too large stab values.

**stabilize_behavior

as described previously this command may be used to improve convergence when both material and contact non-linearities are included at the same time.

Example#

***continuum_contact
 **zone penalty
  *slave imp
  *master tar
  *warning_distance 1.0
  *automatic_penalty

***continuum_contact

On this page

***continuum_contact#