16 Contact
16.1 Mechanical Contact
- Use ArborX to find the collision candidates. For each triangle, we look for the triangles that satisfy the following conditions:
Must be within the activation distance of the triangle
Is not directly connected to the triangle based on whether vertex IDs are shared
Is not filtered out by the collision filters (See Collision Detection).
- We now compute a double integral of the force between the two surfaces:
- Iterate over all the triangles:
- Iterate over all the quadrature points of that triangle (We use a Gauss-Lobatto rule by just choosing the vertices):
- Iterate over the candidate triangles:
Find the barycentric coordinates of the nearest point on the current candidate triangle
Calculate the position of the quadrature point relative to the nearest point
Calculate the inner integral weight. This is just the area of the candidate triangle within a ball about the quadrature point whose radius is the activation distance. This helps to normalize the contributions from different sized triangles
Calculate effective distance.
Compute the potential gradient by passing the effective distance to the 1d potential derivative.
Compute the outer quadrature weight as one-third the area of the triangle that the quadrature point comes from
Compute the force between the quadrature point and the candidate triangle as the product of the potential gradient, the normal of the candidate triangle, and the two quadrature weights.
If signed-distance force scaling is turned on, the force is then scaled based on how far the point is from the closest point on any of the candidate triangles.
Assemble this vector into the residual for the outer triangle and then negate and scale it by the barycentric coordinates on the candidate triangle to find the contributions to the candidate triangle
The triangle residual is then projected into the spline function space through an extraction matrix
A couple of caveats here.
The actual looping structure moves the outer triangle loop down to the bottom of the looping because the quadrature rule is a Gauss-Lobatto style rule where the quadrature points are on the boundaries between triangles. This means that a lot of calculations are shared across triangles.
There is no actual potential associated with this method. It was inspired by a double integral of a potential, but the inputs to the potential are extensively manipulated and then the results from the potential are further modified. We maintain force balance like a potential, but there is no potential to match.
The aforementioned manipulations are directly made into the normal contact traction shown in the Euler-Lagrange equation derived from the defined double-integral contact potential along with the standard solid mechanics potential. With the manipulations, there is no potential that explicitly matches the weak form any more. This decision was made to maintain a continuous contact force field as well as easy implementation.
If we successfully develop a true or smoother potential (the multiplicative approach or something along that line is a promising candidate), meaning that it’s smooth throughout the 3-D space, we will be able to seamlessly derive the variational formulation in the context of potential minimization. Working with a potential also directly implies force balance. Force balance is an important property of any contact method so that one side does not experience excessive acceleration compared to the other.
16.1.1 Kinematics
Given a motion of a group of solid bodies and two points and in the reference configuration of the bodies, we denote the motion evaluated at those two points by and . Likewise, given a displacement field we denote the displacement field evaluated at those two points by and . The boundary of all bodies in the reference configuration is denoted by . The motion can be written in terms of the displacement as The signed distance between and is where is the surface normal at .
The variation of the signed distance is
16.1.2 Equilibrium
Given a one-dimensional potential, , we define the energy associated with contact as where is the portion of the surface that is influenced by the point . More specifically, is chosen to be the intersection where denotes the activation ball of radius , which is called activation distance, centered at . satisfies the collision detection requirements enumerated in Collision Detection. Note those requirements are simplified to the case of a triangular mesh.
Ignoring the variation of in equation , the first variation of the contact energy is approximately where
In the penalty formulation, the one-dimensional potential is defined as with the contact stiffness serving as a penalty parameter. Note that the expression includes to normalize the inner integral in equation so that the penalty parameter of the contact formulation has the unit consistent with the standard penalty Dirichlet BC formulation.
The simplest case is the quadratic potential:
In order to avoid the discontinuity at in the 1st derivative and improve the performance of nonlinear convergence, we can consider higher order potential, e.g. a cubic potential: or alternatively, a locally smoothed quadratic (LS-Q) potential: where denotes the distance smoothing factor. This potential is designed such that the discontinuous second derivative of the quadratic potential is regularized by a cubic Hermite polynomial within the range as follows: The corresponding first derivative is The above potential reduces to the quadratic potential when is set to zero.
16.1.3 Frictionless residual construction
Based on equation , the contact residual is defined as
Discretizing the solution, the residual then becomes where and are the spline basis functions on the surface at points and , respectively.
In order to take advantage of the speed of nearest point searches on linear triangles, we approximate the spline surfaces with linear tessellations. This approximation is done through extraction. The spline basis functions can be approximated with a linear operator known as the extraction operator. where is the linear basis function with index on a triangle containing point , and similar for . This approximates the basis functions, but we can also approximate the bounds of integration by applying the extraction operator to the spline coefficients. Note that, by the linearity of the integral, we can pull the extraction operator outside. Denoting the tessellation of and by and , respectively, the force then becomes where For simplicity, is just a subset of the triangles in .
We approximate each of the two integrals with quadrature rules. The outer integral is approximated with a simple "off-the-shelf" triangle quadrature rule. We approximate the inner integral by finding the nearest point on the triangle to the quadrature point on the triangle . This forms a single point quadrature rule for each point in the outer quadrature rule. The corresponding quadrature weight is just the area of the intersection , denoted . The benefit of this quadrature weight is that the force smoothly goes back to zero once we have penetrated too far through the surface.
The residuals then reduce to where is the area of the triangle and is the quadrature weight of the outer integral quadrature rule. Choosing a three point Gauss-Lobatto quadrature rule this reduces to where is the location of vertex on triangle .
The summation over the ordered pairs of triangles is greatly reduced through the use of a spatial search tree to eliminate triangles that are too far apart (See Collision Detection).
16.1.4 Frictional residual construction
We introduce the nominally tangential contribution, corresponding to the nominally normal contact traction contribution defined in equation , where the tangential relative velocity is and the unit tangent vector is The function is a regularization of Coulomb friction that determines the magnitude of frictional traction as a function of the magnitude of tangential velocity and normal pressure. A smooth transition of the frictional traction is considered to improve convergence of implicit time stepping. An example is where denotes a regularization velocity and is interpreted as the maximum slip velocity allowed during what should nominally be static friction.
The tangential traction contribution is added to its associated normal contribution to obtain the full residual vector for frictional contact.
16.1.5 Collision Detection
Spurious self-contact filter: There are surface configurations that do not involve surface interpenetration, but the contact algorithm mistakenly considers them to be in a penetrating state. For example, in the case of a thin wall with a thickness smaller than the activation distance, the opposite surface may be incorrectly treated as a penetrating surface, resulting in a large, spurious contact force. Another example is the case of a sharp edge where two surfaces meet, and their normals produce a negative dot product. To resolve this, throughout the entire simulation, we filter out candidate triangles on the same body as the owned triangle that are part of the initial non-filtered collision set and initially produce a negative signed distance as defined in equation .
Angle filter: For computational efficiency, we filter out triangles whose normal vectors have an angular difference of less than a certain value, such as , with the normal of the owned triangle. This filter helps avoid unnecessary preallocation in the tangent matrices. Note that only triangle pairs with a negative dot product between their normals generate contact forces. Therefore, we exclude triangles that are not expected to generate contact forces in the next few steps. When this filter is used, the need to repeat collision detection is also determined by the angular changes of the triangles on the contact surfaces. This allows us to re-include force-generating triangles back into the collision sets.
Area filter (beta feature): Contact can become a very memory-intensive procedure when too many candidate triangles are included in the collision sets, leading to excessive preallocation in the tangent matrix. Diverging solution results in large contact triangles that may grab an excessive amount of candidate triangles. In the current workflow, memory blowups have been observed frequently in large problems during the tangent preallocation stage because the collision search and preallocation are triggered before PETSc reduces the time step size in response to a diverged residual norm. When this filter is used, the collision search and preallocation are skipped if the area of any triangle on the contact surface has significantly increased compared to the last collision search step. This is a beta feature, intended for use until a more effective workflow is developed, but it has already helped some large problems avoid memory blowups.
16.1.6 Contact force modifications
16.1.6.1 Effective signed distance modification
While simple, the signed distance defined in equation can produce multiple issues in practical applications.
First, we scale the signed distance to exclude the contribution from the surface segments on that are not facing towards . where The kink in affects the convergence only near of angle difference, but a smoother version of can be considered.
The next consideration is the contact force transition near the triangle edges. In the current approach of additive potential, we keep the force contribution from point to a triangle, even when the closest point falls on an edge of the triangle. This contribution smoothly fades out as moves further away from the triangle, because the domain of potential’s inner integral in equation becomes smaller. In spite of it, the signed distance defined in equation can generate spurious contact forces in many cases, as long as the triangle is intersected by the activation ball of . An example is that a curved surface is in contact with a plane, where a point on the plane lies outside of the opposite body but generates the signed distance equation with respect to an opposite triangle is negative. To resolve this issue, the following modified signed distance can be considered. where with being the point on the plane of whose orthogonal projection onto the opposite triangle is .
Note that this modification does help resolve the spurious contact forces from curved surface-to-plane type but does not fix all the edge cases, for which the force scaling is available, which is more expensive but effectively removes spurious contact forces. One advantage of this fix over the force scaling is that this does not increase the bandwidth of the system while the the force scaling does.
The aforementioned manipulations are directly made into the normal contact traction shown in equation for simple implementation.
16.1.6.2 Contact force scaling
The motivation of the treatment described in this section is to resolve the spurious contact force generated on arbitrary geometry. The contact force scaling is a multiplicative correction to the normal contact traction in equation . The advantage of this approach is that any spurious contact forces are completely removed as long as the contact surface is not terminated within the current set of candidates. Note that if there is termination in the contact surface within the current set of candidates, the inside/outside determination could be incorrect. A downside of this approach is that the bandwidth of the tangent matrix increases, as there are test vertex-trial vertex connections between candidate vertices.
The normal contact traction contribution in equation is scaled as follows: where is the force scaling factor that is a continuous functional in terms of the closest distance , the signed distance to the closest point on the opposite surface that is the union of all integration cells within the current set of candidates. Note that the closest distance itself is a continuous function, so is a continuous function. However, with faceted triangular representation of the contact surface, is only . Therefore, additional smoothing could be helpful for nonlinear convergence.
The force scaling factor is a regularized step function and, in the current implementation, utilizes a cubic Hermite smoothing in equation as follows: where a small factor that defines the range over which the force is scaled.
16.1.7 Interference fits
The resolution of interference fit relies on the regular contact algorithm to find the equilibrium configuration by removing the initial penetration under the quai-static condition. In using the regular contact algorithm, there are two main issues:
Two inter-penetrating surfaces cannot see each other when the initial inter-penetration is larger than the desired activation distance with which the remaining simulation runs after resolving the interference
The initial contact force can be intractably large with an initially large inter-penetration, which likely leads to failure in the nonlinear convergence.
The first issue is addressed by setting a large initial activation distance that is large enough for the penetrating surfaces to see each other and ramping it down to the desired value over a specified time period. Likewise, the second issue is addressed by setting a small initial contact stiffness to avoid divergence of the nonlinear solver and ramping it up to the desired value over a specified time period. Currently, the following ramping functions are implemented.
where , , , and denote the activation distance at time , initial activation distance, final (desired) activation distance, and ramp time interval.
where , , and denote the contact stiffness at time , initial contact stiffness, and final (desired) contact stiffness.
16.2 Tied Contact
Let be the inversion on the master geometry of a point on the slave geometry denoted by . The displacement of the slave point at the master point can be approximated with We then enforce the constaint
16.2.1 Residual
The residual at a quadrature point is