# 16. Contact

section{Contact virtual work}

As a starting point for the treatment of contact, its contribution to the virtual work expression can be stated as:

(16.1)\[\int_{S^3} \left( -t_N \delta g_N + t_{T_\alpha} \delta g_T^{\alpha} \right) \text{d} a ,\]

where \(S^3\) is the common surface between two continua, \(t_N\) is the contact normal traction (positive in compression) \(t_{T_\alpha}\) is the contact tangential traction in one of two local (tangent plane) directions \(\alpha\), and $delta g_N:math:` and delta g_T^{alpha}:math: are the directional derivatives of the contact normal gap g_N$ and tangential slip :math:`g_T^{alpha} in the direction of \(\dot{\varphi}\), i.e.:

(16.2)\[\begin{split}\begin{aligned} \delta g_N &:= \left.\frac{\textrm{d}}{\textrm{d}\beta}\right|_{\beta=0} \left[ g_N \left(\varphi +\beta \dot{\varphi}\right) \right] \\, \delta g_T^{\alpha} &:= \left.\frac{\textrm{d}}{\textrm{d}\beta}\right|_{\beta=0} \left[ g_T^{\alpha} \left(\varphi +\beta \dot{\varphi}\right) \right] . \end{aligned}\end{split}\]

In (16.1), the deformation is subject to the following constraints, referred to as the Kuhn-Tucker conditions. The Kuhn-Tucker conditions are a set of constraints to be considered representative of the mechanical contact problem in continuum mechanics, and can be written as:

(16.3)\[\begin{split} \begin{aligned} \delta g_N \geq 0 & \text{ (a) impenetrability constraint},\\ t_N \geq 0 & \text{ (b) no adhesion condition},\\ t_N g_n \geq 0 & \text{ (c) complementary condition}, \\ t_N \dot{g}_N \geq 0 & \text{ (d) persistency condition}, \end{aligned}\end{split}\]

for frictionless response and

(16.4)\[\begin{split} \begin{aligned} \Phi := \| t_T \| - \mu t_N \leq 0 & \text{ (a) slip function}, \\ L_v g_T -\varsigma \frac{t_T}{\|t_T\|} = 0 & \text{ (b) slip rule},\\ \varsigma \geq 0 & \text{ (c) consistency parameter},\\ \Phi \varsigma = 0 & \text{ (d) complementary condition}, \end{aligned}\end{split}\]

for the prescription of a Coulomb friction (where \(\mu\) is the friction coefficient). In (16.3), the gap \(g_N\) is defined with respect to all material points \(\mathbf{Y} \in S^3\) as:

(16.5)\[g_N ( \mathbf{X},t) = \min_{\mathbf{Y} \in S^3} \| \varphi ( \mathbf{X} ) - \varphi ( \mathbf{Y} ) \| *sign* (g_N),\]

where

(16.6)\[\begin{split}\textrm{sign} (g_N) = \begin{cases} -1 & \text{if } \varphi (\mathbf{X}) \text{ lies on the interior of the contacted body}, \\ 1 & \text{otherwise} . \end{cases}\end{split}\]

The tangential gap rate \(L_v g_T\) in (16.4) is defined as follows:

(16.7)\[ L_v g_T = \left( \varphi \left( \mathbf{X} \right) - \varphi \left( \bar{\mathbf{Y}} (\mathbf{X}) \right) \right) \cdot \left( \mathbf{p}^{\alpha} \otimes \mathbf{p}_{\alpha} \right),\]

where \(\mathbf{p}^{\alpha}\) and \(\mathbf{p}_{\alpha}\) are base vectors associated with any appropriate surface coordinate system used to describe \(S^3\), with these base vectors being evaluated at the current contact point \((\bar{\mathbf{Y}}(\mathbf{X}))\) that satisfies the minimization of (16.3). Use of the notation \(L_v\) is meant to imply a Lie derivative, which can be understood to be the time derivative of an object as viewed from an embedded reference frame, in this case the convected frame \(\mathbf{p}_{\alpha}\) frame, that moves along with the point.

## 16.1. Discretized forms of contact constraints

The question is then, how to represent these conditions in discretized form suitable for FE solution methods. A simple example, shown in Fig. 16.1.1, serves to demonstrate the concern that this question embodies. Two discretizations for the interface are evident and, as this simple example indicates, leads to an ambiguous definition of the interface.

Fig. 16.1.1 Concerns in constraints choices for contact problems.

A historical treatment of contact has focused on applying the Kuhn-Tucker condition directly to the discretized form, leading to what we will be referring to here as a node-face treatment of contact, or node-face contact. As we will review here in this chapter, this treatment of contact is relatively straightforward from a conceptual standpoint, however it does have several issues - even to the point of the overall approach being pathological in some applications.

Alternatively, more recent investigations have focused on addressing these issues, leading to what we will be referring to here as a face-face treatment of contact, or face-face contact. These methods consider the weak form more directly, thus leading to a variationally consistent approach (e.g., mortar methods are an example of this approach and are, at the moment, prevalent in the literature).

### 16.1.1. Node-Face contact

For node-face contact the Kuhn-Tucker conditions are assumed to apply to one side of the contacting surfaces.Thus the gap \(g_N\) is defined with respect to all nodal points \(\mathbf{Y}_I\) as:

(16.1.1)\[ g_N ( \mathbf{X}_I,t) = \min_{\mathbf{Y}_I \in S^3} \| \varphi ( \mathbf{X}_I ) - \varphi ( \mathbf{Y} ) \| *sign* (g_N),\]

where \(I\) refers to a nodal point on one side of the interface, whose coordinates are \(\mathbf{X}_I,t\) at time \(t\) of interest. The right-hand side of (16.1.1) is the discrete form of (16.5) but is more commonly called the closest-point projection, which will be discussed in some detail in Section 16.2.3.

As mentioned, there are issues associated with node-face constraints. They stem from the application of contact constraints directly to the discretized problem. As shown in Fig. 16.1.1, the potential to over constrain the interface is avoided by applying the impenetrability constraint only at selected points along the interface. In the Solid Mechanics module these points coincide with the nodes, as it makes it convenient to obtain contact results (normal and tangential tractions, stick/slip results, etc.) and interpret them in post-processing.

However, this approach does not truly alleviate over-constraining. This is easily demonstrated with an enlightening example (we will make use of this example for the discussion of node-face contact and face-face contact, so making a proper introduction is worthwhile). Fig. 16.1.2, shows a beam bending problem that is being modeled with continuum elements, in this case hex8 elements. The beam is cantilevered at its left end appropriately, i.e., fixed at the neutral axis and constrained from motion only in the x-direction elsewhere.

Fig. 16.1.2 A continuum beam subjected to pure bending.

The analytic solution to this problem is one where the neutral axis should take the displacement corresponding to an arc of a circle. When the moment is prescribed to be \(M^*\) the beam should deform into a perfect circle.

When the beam is meshed with either an all coarse mesh (4 elements through its thickness) or an all fine mesh (16 elements through its thickness), the Finite Element results appear to be quite acceptable, producing the pure bending solution.

However, lets now combine coarse and fine discretizations to solve the problem. In this case a mesh tying constraint is required to obtain the solution, which is seen to be fundamentally a contact problem with adhesion and infinite frictional capacity. The combining of coarse and fine discretizations can be done is a couple of canonical ways, as shown in Fig. 16.1.3; one where the interface between the discretizations is vertical and the other where is along the neutral axis.

Fig. 16.1.3 A continuum beam that includes mesh tying subjected to pure bending.

In both cases, we apply the standard rule of thumb: given the same material on both sides of the interface, apply the contact constraints on the finer discretization. Subjecting the beam to the prescribed moment reveals at once the issue: kinematically enforcing a zero gap condition at each node is exactly correct in one case, where the interface is through the depth of the beam, and severely over constraining in the other, where the interface is along the neutral axis. As Fig. 16.1.4 shows, the over constraint can be severe and may produce spurious stress distributions in the fine mesh, particularly near the neutral axis.

Fig. 16.1.4 Results for a continuum beam that includes mesh tying subjected to pure bending.

## 16.2. Contact Search

The contact search algorithm is a logical component of the overall contact capability. Much of the reason to consider search as a separate component is due to a need to revisit and replace algorithms as they demonstrate better performance. As problem sizes grow there is an increasing computational cost of this aspect of computational solid mechanics, particularly on distributed memory (parallel) computers.

As a way of introducing the concepts inherent in contact search algorithms, we recognize the similarity of contact search to many other other problems in the simulation domain (e.g., the video gaming industry). In this more abstract sense, a significant part of the contact search algorithm is a proximity determination of one object with respect to another. Collision detection, as it is also referred to, is computationally intensive but also studied thoroughly to obtain the best performance possible. Thus the contact algorithm in Sierra/SolidMechanics is comprised of the more general proximity search followed by the much more specific detailed detection of contact in the context of a discrete Finite Element method.

### 16.2.1. Proximity search algorithms

Although various proximity search algorithms have been developed over the years, those that have been used in Sierra/SolidMechanics are discussed. Proximity search algorithms deal with bounding boxes. Construction of bounding boxes are straightforwardly computed as the vector of min/max coordinates of the volume swept by the predicted motion of either nodes or faces over a time step. Specifically, over the time step \(( t \rightarrow t+\Delta t )\), the axis-aligned bounding box for node \(I\) is computed as follows:

(16.2.1)\[ aabb(N_I) = \min( \varphi(\mathbf{X}_I , t) , \varphi(\mathbf{X}_I , t+\Delta t) ) \text{ , } \max( \varphi(\mathbf{X}_I , t) , \varphi(\mathbf{X}_I , t+\Delta t) ).\]

The resulting data of an axis aligned bounding box calculation is the absolute minimum amount of data (min and max \(x,y,z\) coordinates) representing a contact entity. The example expressed in (16.2.1) is for a node; with the extension to a face being straightforwardly computed as:

(16.2.2)\[ aabb(F_M) = \min( aabb(N_J), \text{ }J=\text{nodes of face } M )\text{ , } \max( aabb(N_J), \text{ }J=\text{nodes of face } M ).\]

This allows a simplification or reduction of the contact problem to the more general proximity detection based solely on axis-aligned bounding boxes of contact entities. Many of the structural modeling capabilities within Sierra/SM are converted to these primitives (i.e., nodes and faces). Beams and shells, for example, do not explicitly represent volume, but a volume is inferred with ancillary data such as thickness. These elements are converted to contact primitives by lofting the finite element geometry explicitly to volumetric discretizations.

### 16.2.2. Parallel search algorithms

The strategy for computational simulation of contact on distributed memory architectures (parallel computing) is to decompose the contact problem in a distributed manner among the compute nodes. Typically known as domain decomposition methods, there are several that are directly applicable to the contact problem. Inertial decomposition and recursive coordinate bisection (RCB) decomposition are two geometric-based algorithms that are examples. This geometry-based decomposition approach is depicted graphically in Fig. 16.2.1.

Fig. 16.2.1 Simple illustration of the domain decomposition for contact problems.

The proximity search algorithm thus performs with parallel scalability and serial efficiently on each processor.

### 16.2.3. Contact kinematics

Recall that the output of the proximity search is a collection of bounding box pairs whose volumes overlap. The contact entities associated with the bounding boxes are then considered for contact in what we call a detailed search. Detailed searching is a term applied to computing the contact kinematic quantities associated with either the node-face or face-face algorithm.

Closest point projection for node-face contact

The underpinning of the node-face contact approach is the choice of a set of points at which to apply the Kuhn-Tucker conditions. Once this choice is made, it follows that a contact point must be determined for each node. The closest point projection is the name given to this calculation, which is simply the point on the opposing surface that minimizes the gap, i.e.,

(16.2.3)\[ g_N ( \mathbf{X_I},t) = \min_{\mathbf{Y} \in \{ *facets on side A surface* \} } \| \varphi ( \mathbf{X_I} ) - \varphi ( \mathbf{Y} ) \| *sign* (g_N).\]

This can be seen simply as the discrete form of the gap function expressed in (16.3). Immediately with the discrete surface not possessing \(C_1\) continuity, the closest point projection is no longer unique. At and around edges and corners are the regions on the side A surface where the issues arise. Fig. 16.2.2 depicts a simple illustration of the non-uniqueness encountered.

Fig. 16.2.2 Simple illustration of the non-uniqueness in the closest point projection.

Minimum volume overlap for face-face contact

The face-face contact approach seeks to avoid these issues by considering the volume overlap between the discrete sides of the interface.