Contact Boundary Conditions¶

Formulation¶

Rigid-body contact boundary conditions are implemented in Ratel using a Nitsche method approach as described in [Mli18], which we describe in detail for completeness. Currently, a rigid plane (platen) and rigid cylinder are supported.

Deformable-deformable contact is in development, see Deformable-Deformable Contact below.

Large-Deformation Frictional Contact¶

Given a boundary region \(\Gamma_c\), the gap between the rigid shape \(\mathcal{B}\) and a deformed configuration point \(\bm{x}\in\Gamma_c\) is given by the distance between \(\bm{x}\) and the nearest point on the shape \(\bm{y}\in\partial\mathcal{B}\). That is,

\[ \bm{y} = \arg\min_{\bm y\in\partial\mathcal{B}} \| \bm x - \bm y \|_2, \]

and the signed gap function is given by

\[ g(\bm{x}) = \bm{n}_y \cdot \left( \bm{x} - \bm{y} \right). \]

We can formulate the gap in terms of the displacement field \(\bm u\) as

\[ g(\bm u) = \bm{n}_y \cdot(\bm{u} + \bm X - \bm c) = \bm{n}_y \cdot \bm u + \bm{n}_y (\bm X - \bm c) = \bm{n}_y \cdot \bm u + g_0, \]

where \(g_0\) is the gap in the initial configuration.

Denote the contact traction by \({\bm{\sigma}}_N(\bm{u})\), which is computed as the first Piola-Kirchoff traction vector \({\bm{P} \cdot \hat{\bm N}}\), where \(\hat{\bm N}\) is the normal at the material point \(\bm{X} = \phi^{-1}(\bm{x})\). Note that regardless of whether the volume operator is computed using the initial or current configurations, the contact traction is always computed in the initial configuration to simplify the linearization. The contact traction is then decomposed into normal and tangential components, relative to the normal \(\bm{n}_y\), as

\[ {\bm{\sigma}}_{N} = -{\sigma}_n\bm{n}_y + \bm{T}_{\bm{n}_y}{\bm{\sigma}}_{N} = -{\sigma}_n\bm{n}_y + {\bm{\sigma}}_t, \]

where \(\bm{T}_{\bm{v}} = \bm{I} - \bm{v} \otimes \bm{v}\) is the projection operator onto the tangent plane defined by the normal vector \(\bm{v}\). The normal component \({\sigma}_n\) is the contact pressure.

Using these quantities, normal contact requires the following:

(329)¶\[ \begin{aligned} g(\bm u) &\geq 0\\ {\sigma}_n (\bm u) &\leq 0\\ {\sigma}_n(\bm u) g(\bm u) &= 0. \end{aligned} \]

That is, if the point is in contact with the shape, then the contact pressure must be negative; otherwise, the contact pressure must be zero.

For Coulomb friction, we also require:

(330)¶\[ \begin{cases} \|{\bm{\sigma}}_t\| \leq -\mathscr{F} {\sigma}_n & \text{if } \bm{\mathrm{v}} = 0,\\ {\bm{\sigma}}_t = \mathscr{F} {\sigma}_n \frac{\bm{\mathrm{v}}}{\|\bm{\mathrm{v}}\|} & \text{otherwise}, \end{cases} \]

where \(\bm{\mathrm{v}}\) is the relative slip velocity defined by

\[ \bm{\mathrm{v}} = \dot{\bm{x}} - \dot{\bm{y}} - g \dot{\bm{n}}_y, \]

and \(\dot{\bm{x}} = \dot{\bm{u}}\).

Adopting the notation of [Mli18], we define the projection operator \(P_{B(\bm{n}, r)}(\bm{v})\) which first projects the vector \(\bm{v}\) onto the tangent plane defined by the normal vector \(\bm{n}\), then projects the result onto a ball of radius \(r\). Specifically, that is,

\[ P_{B(\bm{n}, r)}(\bm{v}) = \begin{cases} \bm{T}_{\bm{n}} \bm{v} & \text{if } \|\bm{T}_{\bm{n}} \bm{v}\| \leq r, \\ r \frac{\bm{T}_{\bm{n}} \bm{v}}{\|\bm{T}_{\bm{n}} \bm{v}\|} & \text{otherwise}. \end{cases} \]

Additional friction models are also supported, see Friction Models.

Nitsche Method¶

The Nitsche method reformulates the conditions (329) by defining the Nitsche-penalized normal contact pressure as

(331)¶\[ \hat{\sigma}_n = \left[ \sigma_n + \gamma g \right]_{\mathbb{R}^-}, \]

where \([v]_{\mathbb{R}^-} = \frac{1}{2}(v-|v|)\) is the projection onto the negative real numbers and \(\gamma > 0\) is the Nitsche method parameter. The Nitsche method parameter \(\gamma\) serves a similar role to a penalty method, and larger values of \(\gamma\) result in more robust non-penetration, but may hinder convergence. A reasonable choice of \(\gamma\) is \(200\cdot E\) where \(E\) is the Young’s modulus of the elastic material. In general, \(\gamma\) need only be a positive function defined on \(\Gamma_c\).

The friction conditions (330) are similarly reformulated as

\[ \hat{\bm{\sigma}}_t = P_{B\left(\bm{n}_y, -\mathscr{F\left[ {\sigma}_n + \gamma g \right]_{\mathbb{R}^-}}\right)}({\bm{\sigma}}_N - \gamma \bm{\mathrm{v}}). \]

Weak Form¶

For simplicity, we only describe the boundary part of the weak form, which replaces the final (boundary) term in, e.g. (93).

There are several variants of the Nitsche method formulation depending on an auxiliary parameter \(\theta\in\mathbb{R}\). Typically, \(\theta=-1,0,1\) are the values considered. Currently, Ratel implements boundary conditions with \(\theta=0\), which results in the simplest Jacobian. The weak form of the boundary conditions on \(\Gamma_c\) can be written as:

(332)¶\[ \int_{\Gamma_c} \bm{v}\cdot\left(\left[ {\sigma}_n + \gamma g \right]_{\mathbb{R}^-}\bm{n}_y - P_{B\left(\bm{n}_y, -\mathscr{F}\left[ {\sigma}_n + \gamma g \right]_{\mathbb{R}^-}\right)}({\bm{\sigma}}_N - \gamma \bm{\mathrm{v}})\right) \ \diff S \]

Newton Linearization¶

In order to solve using Newton-Krylov solvers, we require the linearization of (332).

Each shape provides its own function to compute \(\diff \bm{y}\) and \(\diff \bm{n}_y\) for a given \(\bm{x}\) and variation \(\diff\bm{x} = \diff\bm{u}\), which depending on the shape may be quite complex. Given these, the linearization of the frictionless case is straightforward to compute. Let \(H\) be the Heaviside function, defined as

\[ H(x) = \begin{cases} 0, & x < 0;\\ 1, & x \geq 0. \end{cases} \]

Then, we can compute the variation in the contact pressure as

\[ \begin{aligned} \diff \hat{\sigma}_n = -\diff \left( \bm{P} \cdot \hat{\bm{N}} \cdot \bm{n}_y\right) &= -\diff \bm P \cdot \hat{\bm{N}} \cdot \bm{n}_y - \bm P \cdot \hat{\bm{N}} \cdot \diff \bm{n}_y\\ &= -\diff {\bm{\sigma}}_N \cdot \bm{n}_y - {\bm{\sigma}}_N \cdot \diff \bm{n}_y, \end{aligned} \]

where \(\diff \bm{P}\) is the variation in the first Piola-Kirchhoff stress tensor computed via the constitutive model. The variation in the signed normal gap \(g\) is given as

\[ \diff g = \diff \left( \bm n_y \cdot(\bm{x} - \bm{y}) \right) = \bm n_y \cdot (\diff \bm x - \diff \bm y) + \diff \bm{n}_y \cdot (\bm{x} - \bm{y}). \]

Combining these, the variation in the magnitude of the Nitsche contact force is

(333)¶\[ \diff\left(\left[ {\sigma}_n + \gamma g \right]_{\mathbb{R}^-}\right) = H\left(-{\sigma}_n - \gamma g\right) ((\diff \bm x - \diff \bm y - \diff {\bm{\sigma}}_N) \cdot \bm{n}_y + (x - y - {\bm{\sigma}}_N) \cdot \diff \bm{n}_y). \]

Finally, the linearization of the frictionless contact traction is

\[ \diff\left(\left[ {\sigma}_n + \gamma g \right]_{\mathbb{R}^-}\bm{n}_y\right) = \diff\left(\left[ {\sigma}_n + \gamma g \right]_{\mathbb{R}^-}\right) \bm{n}_y + \left[ {\sigma}_n + \gamma g \right]_{\mathbb{R}^-}\diff\bm{n}_y. \]

The case with friction requires some additional consideration.

The linearization of the relative slip velocity is

\[ \diff \bm{\mathrm{v}} = \diff \dot{\bm u} - g (\diff \dot{\bm{n}}_y) - (\diff g) \dot{\bm{n}}_y, \]

as since the shape is rigid, \(\diff \dot{\bm{y}} = 0\).

Denote \(\bm{q} = {\bm{\sigma}}_N - \gamma \bm{\mathrm{v}}\), \(\bm{q}_t = \bm{T}_{\bm{n}_y} \bm{q}\), and \(r = -\mathscr{F}\left[ {\sigma}_n + \gamma g \right]_{\mathbb{R}^-}\). Then, the frictional part of the weak form is rewritten as

\[ -\int_{\Gamma_c} \bm{v}\cdot P_{B\left(\bm{n}_y,\ r \right)}(\bm q) \ dS \]

The linearization of the friction term is then

\[ \mathrm{d} \left(P_{B(\bm{n}_y, r)}(\bm q)\right) = \begin{cases} 0, & \text{if } r \leq 0 \text{ or } \|\bm{q}_t\| = 0 \\ \mathrm{d}\bm{q}_t, &\text{if } \| \bm{q}_t \| \leq r \\ \frac{r}{\|\bm{q}_t\|} \left(I - \frac{\bm{q}_t}{\|\bm{q}_t\|} \otimes \frac{\bm{q}_t}{\|\bm{q}_t\|} \right) \mathrm{d}\bm{q}_t + \frac{\bm{q}_t}{\|\bm{q}_t\|}\mathrm{d}r, & \text{otherwise}. \end{cases} \]

Note that \(\diff r\) is given in (333) and \(\diff \bm{q}\) can be written as

\[ \diff \bm q_t = \diff \bm q - (\bm q \cdot \bm n_y) \diff \bm n_y - (\diff \bm q \cdot \bm n_y) \bm n_y - (\bm q \cdot \diff \bm n_y) \bm n_y, \]

where \(\diff \bm q = \diff {\bm{\sigma}}_N - \gamma \diff \bm{\mathrm{v}}\).

The full linearization for contact with friction is then given as

\[ \begin{aligned} &\int_{\Gamma_c} \bm{v}\cdot \left(\diff\left(\left[ {\sigma}_n + \gamma g \right]_{\mathbb{R}^-}\right) \bm{n}_y + \left(\left[ {\sigma}_n + \gamma g \right]_{\mathbb{R}^-}\right) \diff \bm{n}_y\right) \ \diff S \\ &\quad- \int_{\Gamma_c} \bm{v} \cdot \diff \left(P_{B\left(\bm{n}_y, -\mathscr{F}\left[ {\sigma}_n + \gamma g \right]_{\mathbb{R}^-}\right)}({\bm{\sigma}}_N - \gamma \bm{\mathrm{v}})\right)\ \diff S. \end{aligned} \]

Shapes¶

Each shape requires an implementation of RatelBCContactProject, which computes the projection of a given point onto the nearest point \(\bm{y}\) on the contact surface and the normal to the contact surface \(\bm{n}_y\). Additionally, all shapes support linearly-interpolated rigid-body translation through space, specified via the -bc_contact_[face]_translate option (see Contact for more details on the command-line interface). To support these rigid-body motions, each shape requires an implementation of RatelBCContactUpdateModel. This function takes a given time and returns an updated shape-specific struct with the geometric properties of the shape at that time.

Platen¶

The platen is modeled as an abstract half-space with a point \(\bm{c}\) and a normal vector \(\bm{n}_p\), which is displaced by a distance \(d > 0\) along \(\bm{n}_p\) linearly throughout the simulation. For the platen, finding the nearest point \(\bm{y}\in\partial\mathcal{B}\) to a given current configuration point \(\bm{x}\) is quite simple. Note that since the platen has the same normal everywhere, \(\bm{n}_y = \bm{n}_p\) and the normal gap can be computed trivially as

\[ g(\bm{x}) = \bm{n}_p \cdot \left(\bm{x} - \bm{c}\right). \]

Then, the nearest point is given by subtracting the normal gap from \(x\) as:

\[ \bm{y} = \bm{x} - g(\bm{x}) \bm{n}_p. \]

The linearization \(\diff \bm{y}\) is then

\[ \diff\bm{y} = \diff\bm{x} - (\diff\bm{x} \cdot \bm{n}_p) \bm{n}_p. \]

For a given time \(t\) with prescribed translation \(\bm{T}\), the geometric properties are updated as

\[ \begin{aligned} \bm{n}_p(t) &= \bm{n}_p\\ \bm{c}(t) &= \bm{c} + \bm{T}, \end{aligned} \]

that is, the center \(\bm{c}\) is translated. Platens also support a simplified translation interface via the -bc_contact_[face]_distance [float list] option, which prescribes the translation as a distance along the normal vector (see Platen for more details).

Cylinder¶

The cylinder is modelled as an infinite cylinder with radius \(r\) and linear axis defined by a normalized vector \(\bm{a}\) and a point on the center line \(\bm{c}\). The direction of the normal vector for the cylinder, i.e. inward or outward, can be controlled via -bc_contact_[face]_inside [bool] flag.

The projection onto the surface of the cylinder requires first computing the normal vector from the center axis of the cylinder to the point \(\bm{x}\in\Gamma_c\). Let \(\tilde{\bm{n}}_c = \bm{x} - \bm{c}\). Note, since \(\bm{c}\) is an arbitrary point on the center line, \(\tilde{\bm{n}}_c\) is not normal to the cylinder in general. Thus, the unnormalized normal vector can be computed by projecting \(\tilde{\bm{n}}_c\) onto the tangent plane to \(\bm{a}\), that is,

\[ \bm{n}_c = \tilde{\bm{n}}_c - (\tilde{\bm{n}}_c \cdot \bm{a}) \bm{a}. \]

If the point \(\bm{x}\) exactly coincides with the center line of the cylinder, \(\bm{n}_c\) will be the zero vector. This possibility can be safely ignored, however, as this would require a far larger timestep than could converge to begin with. Additionally, if a zero vector is used, the residual will simply be zero at that particular quadrature point, and the residual at the surrounding quadrature points will results in sufficient displacement to prevent this from happening in subsequent nonlinear iterations. Finally, the normal vector at the mapped point \(\bm{y}\) is given as

\[ \bm{n}_y = \begin{cases} \sgn(\bm{n}_c),& \quad \text{outside},\\ -\sgn(\bm{n}_c),& \quad \text{inside}, \end{cases} \]

where \(\sgn(\bm{v})\) is defined in (334). Then, the projected point \(\bm{y}\) is given by

\[ \bm{y} = \bm{c} + (\tilde{\bm{n}}_c \cdot \bm{a}) \bm{a} + r \bm{n_y} \]

The linearization \(\diff \bm{n}_y\) can be found by:

\[ \begin{aligned} \diff\tilde{\bm{n}}_c &= \diff \bm{x} \\ \diff\bm{n}_c &= \diff \tilde{\bm{n}}_c - (\diff \tilde{\bm{n}}_c \cdot \bm{a}) \bm{a}\\ \diff\bm{n}_y &= \begin{cases} \diff\sgn(\bm{n}_c),& \quad \text{outside},\\ -\diff\sgn(\bm{n}_c),& \quad \text{inside} \end{cases} \end{aligned} \]

where \(\diff \sgn(\bm{v})\) is given in (339). The linearization \(\diff\bm{y}\) is then

\[ \diff \bm{y} = (\diff \tilde{\bm{n}}_c \cdot \bm{a}) \bm{a} + r \diff \bm{n}_y. \]

Translation of the cylinder is defined in terms of the point \(\bm{c}\). For a given time \(t\) with prescribed translation \(\bm{T}\), the geometric properties are updated as

\[ \begin{aligned} r(t) &= r\\ \bm{a}(t) &= \bm{a}\\ \bm{c}(t) &= \bm{c} + \bm{T}. \end{aligned} \]

Deformable-Deformable Contact¶

Full contact support requires the ability to evaluate contact constraints between two (or more) deformable, meshed surfaces. Work to support mesh-to-mesh contact in Ratel is ongoing.

Segmentation-Based Contact¶

The only method for integrating contact forces between meshes in a robust, accurate, and theoretically rigorous way requires so-called segmentation-based contact [FPW15]. The essential steps of defining the integration cells for segmentation-based contact are illustrated in Figure Segmentation Process:

Define an auxiliary plane based on center and normal of non-mortar face (top left).
Project both the mortar and non-mortar face onto the auxiliary plane (top right).
Compute the intersection of the projected faces (bottom left).
Triangulate the resulting clipping polygon into integration cells (bottom right).

Steps (3) and (4) require specialized algorithms to

../../../../../_images/segmentation.svg — Fig. 9 Segmentation Process¶

Tip

Nearly all geometric operations for testing intersections of points, lines, and polygons require computing the Orient2D function. This function takes three points in 2D, \(A\), \(B\), \(C\), and returns twice the signed area of the triangle they enclose; by convention, the sign is positive if the points are ordered counter-clockwise and negative if they are oriented clockwise.

Naively, this can be implemented via the determinant as:

double Orient2D(double A[2], double B[2], double C[2]) {
    // det([A - C, B - C])
    return (A[0] - C[0]) * (B[1] - C[1]) - (A[1] - C[1]) * (B[0] - C[0])
}

Unfortunately, this naive implementation is numerically unstable, particularly for nearly-collinear points. There exist more expensive mixed-precision approaches to improving the stability of the algorithm, see e.g. [RS97], but such methods require significant branching, making them unsuitable for GPU implementations.

Instead, we can simply apply an error bound derived by [RS97], clamping any values which exceed the bound to zero. The error bound is defined as \(\theta = (1.5 + 4.0 \epsilon_m) * \epsilon_m\), where \(\epsilon_m\) is the machine epsilon for the floating point precision used. Note that \(\theta\) should be precomputed using exact-precision arithmetic. For single and double precision, \(\theta\) has the values:

theta64 = 3.3306690738754716e-16
theta32 = 1.7881399116959074e-07

Then, the more robust Orient2D function can be written as:

CeedScalar Orient2DRobustSingle(CeedScalar A[2], CeedScalar B[2], CeedScalar C[2]) {
    static const CeedScalar theta = sizeof(CeedScalar) == 8 ? 3.3306690738754716e-16 : 1.7881399116959074e-07;
    const CeedScalar det_left = (A[0] - C[0]) * (B[1] - C[1]);
    const CeedScalar det_right = (A[1] - C[1]) * (B[0] - C[0]);
    const CeedScalar det = det_left - det_right;

    return fabs(det) >= THETA * fabs(det_left + det_right) ? det : 0.0
}

Note that while Orient2DRobustSingle prevents the output from having the incorrect sign, it does so by expanding the range of inputs which are considered collinear. The impact of this decision can be somewhat mitigated by computing Orient2DRobustSingle three times and rotating the order of the input points, returning the first accepted value. This is the implementation used in Ratel.

CeedScalar Orient2DRobust(CeedScalar A[2], CeedScalar B[2], CeedScalar C[2]) {
    const CeedScalar det1 = RatelOrient2DRobustSingle(A, B, C);
    const CeedScalar det2 = RatelOrient2DRobustSingle(C, A, B);
    const CeedScalar det3 = RatelOrient2DRobustSingle(B, C, A);

    return det1 != 0.0 ? det1 : (det2 != 0.0 ? det2 : det3);
}

For a visualization of the instability of Orient2D, see this interactive notebook.

Polygon Clipping¶

For Step (3), we must compute the geometric intersection of two polygons. This problem commonly arises in the context of computer graphics, where it is referred to as polygon clipping. We utilize a robust clipping algorithm proposed by [FHP19], which extends upon the classical Greiner-Hormann algorithm [GH98].

Triangulation¶

The original segmentation approaches by [PL04] used a simple method for triangulation: place a point in the center of the clipped polygon, then form triangles with consecutive vertices of the clipping polygon and the center point. While simple to implement, this method has several issues:

The triangulation depends on the current configuration center of the face, so it must be linearized.
The clipped polygon must be convex, which is not guaranteed in general even for linear elements.
The number of triangles is equal to the number of vertices \(n\) on the clipping polygon, even if the clipping polygon is, for example, a triangle. For the most common use-cases of triangle and convex quadrilateral intersections, which have a maximum \(n\) of \(6\) and \(8\), respectively, these unnecessary triangles can result in 20-66% more work than optimal triangulation algorithms.

The authors of [FPW15] suggest a constrained Delaunay triangulation, which ensures a minimal number of triangles (\(n-2\)).

We instead utilize a deterministic, linear-time earcut algorithm by [LCSA22], which is simpler than Delaunay triangulation, but still ensures the optimal \(n-2\) triangles, albeit with potentially more distorted triangles.