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, only a rigid plane (platen) is supported.

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:

(307)\[ \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:

(308)\[ \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 (307) by defining the Nitsche-penalized normal contact pressure as

(309)\[ \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 (308) 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:

(310)\[ \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 (310).

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

(311)\[ \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 (311) 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 prescibes 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 (312). 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 (317). 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} \]