Platen Contact Boundary Conditions

Formulation

Platen-specific contact boundary conditions are implemented in Ratel using a Nitsche method approach as described in [Mli18], which we describe in detail for completeness.

Large-Deformation Frictional Contact

The platen is modeled as an abstract plane 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. Given a boundary region \(\Gamma_c\), the gap between the platen and a deformed configuration point \(\bm{x}\in\Gamma_c\) is given by projection onto the plane of the platen, that is,

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

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

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

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

Denote the contact traction by \(\hat{\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 platen, as

\[ \hat{\bm{\sigma}}_{N}(\bm{u}) = -\hat{\sigma}_n (u) \bm{n}_p + \bm{T}_{\bm{n}_p}\hat{\bm{\sigma}}_{N}(\bm{u}) = -\hat{\sigma}_n (u) \bm{n}_p + \hat{\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 \(\hat{\sigma}_n\) is the contact pressure.

Using these quantities, normal contact requires the following:

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

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

For Coulomb friction, we also require:

(156)\[ \begin{cases} \|\hat{\bm{\sigma}}_t\| \leq -\mathscr{F} \hat{\sigma}_n(\bm u) & \text{if } \dot{\bm{u}} = 0,\\ \hat{\bm{\sigma}}_t = \mathscr{F} \hat{\sigma}_n(\bm u) \frac{\dot{\bm{u}}}{\|\dot{\bm{u}}\|} & \text{otherwise}. \end{cases} \]

Adopting the notation of [Mli18], we define the projection operator \(P_{B(\bm{n}, \delta)}(\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 \(\delta\). Specifically, that is,

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

Nitsche Method

The Nitsche method reformulates the conditions (155) as

\[ \hat{\sigma}_n (\bm{u}) = \left[ \hat{\sigma}_n(\bm u) + \gamma g(\bm u) \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. In general, \(\gamma\) need only be a positive function defined on \(\Gamma_c\).

The friction conditions (156) are similarly reformulated as

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

The resulting formulation for frictional contact is then

\[ \hat{\bm{\sigma}}_N = -\left[ -\hat{\bm{\sigma}}_N \cdot \bm{n}_p + \gamma g \right]_{\mathbb{R}^-}\bm{n}_p + \left[\hat{\bm{\sigma}}_N - \gamma \dot{\bm{u}}\right]_{-\mathscr{F}\left[ -\hat{\bm{\sigma}}_N\cdot\bm{n}_p + \gamma g \right]_{\mathbb{R}^-}}. \]

Weak Form

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

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:

(157)\[ \int_{\Gamma_c} \bm{v}\cdot\left(\left[ -\hat{\bm{\sigma}}_N \cdot \bm{n}_p + \gamma g \right]_{\mathbb{R}^-}\bm{n}_p\right) \ dS - \int_{\Gamma_c} \bm{v}\cdot P_{B\left(\bm{n}_p, -\mathscr{F}\left[ -\hat{\bm{\sigma}}_N\cdot\bm{n}_p + \gamma g \right]_{\mathbb{R}^-}\right)}(\hat{\bm{\sigma}}_N - \gamma \dot{\bm{u}}) \ dS \]

Newton Linearization

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

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

(158)\[ \mathrm{d}\left(\left[ \hat{\sigma}_n(\bm u) + \gamma g(\bm u) \right]_{\mathbb{R}^-}\right) = H\left(-\hat{\sigma}_n(\bm u) - \gamma g(\bm u)\right) \left( -\mathrm{d}\hat{\bm{\sigma}}_{N} \cdot \bm{n}_p + \gamma \mathrm{d} \bm{u}\right). \]

The case with friction requires some additional consideration. First, denote \(\bm{q} = \hat{\bm{\sigma}}_N - \gamma \dot{\bm{u}}\), \(\bm{q}_t = \bm{T}_{n_p} \bm{q}\), and \(\delta = -\mathscr{F}\left[ \hat{\sigma}_n(\bm u) + \gamma g(\bm u) \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}_p,\ \delta \right)}(\bm q) \ dS \]

The linearization of the friction term is then

\[ \mathrm{d} \left(P_{B(\bm{n}_p, \delta)}(\bm q)\right) = \begin{cases} 0, & \text{if } \delta \leq 0 \text{ or } \|\bm{q}_t\| = 0 \\ \mathrm{d}\bm{q}_t, &\text{if } \| \bm{q}_t \| \leq \delta \\ \frac{\delta}{\|\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}\delta, & \text{otherwise}. \end{cases} \]

Since \(\mathrm{d}\delta\) is given in (158) and \(\mathrm{d} \bm{q}\) can be written as

\[ \mathrm{d} \bm{q}_t = \mathrm{d} \left( \hat{\bm{\sigma}}_t - \gamma \dot{\bm{u}}_t \right) = \bm{T}_{\bm{n}_p} \left(\mathrm{d} \hat{\bm{\sigma}}_N - \gamma \mathrm{d}\dot{\bm{u}}\right), \]

the full linearization for contact with friction is given as

\[ \int_{\Gamma_c} \bm{v}\cdot \mathrm{d}\left(\left[ \hat{\sigma}_n(\bm u) + \gamma g(\bm u) \right]_{\mathbb{R}^-}\right) \bm{n}_p \ dS - \int_{\Gamma_c} \bm{v} \cdot \mathrm{d} \left(P_{B\left(\bm{n}_p, -\mathscr{F}\left[ -\hat{\bm{\sigma}}_N\cdot\bm{n}_p + \gamma g \right]_{\mathbb{R}^-}\right)}(\hat{\bm{\sigma}}_N - \gamma \dot{\bm{u}})\right)\ dS. \]