# Semi-smooth Newton methods for mixed FEM discretizations of higher-order for frictional, elasto-plastic two-body contact problems

Heribert Blum, Hannah Frohne, Jörg Frohne, Andreas Rademacher

► To cite this version:

Heribert Blum, Hannah Frohne, Jörg Frohne, Andreas Rademacher. Semi-smooth Newton methods for mixed FEM discretizations of higher-order for frictional, elasto-plastic two-body contact problems. Computer Methods in Applied Mechanics and Engineering, 2016, 309, pp.131-151. 10.1016/j.cma.2016.06.004. hal-01371133

HAL Id: hal-01371133

https://hal.science/hal-01371133

Submitted on 24 Sep 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Semi-smooth Newton methods for mixed FEM discretizations of higher-order for frictional, elasto-plastic two-body contact problems

Heribert Blum, Hannah Frohne, Jörg Frohne, Andreas Rademacher

Chair of Scientific Computing, Technische Universität Dortmund, Vogelpothsweg 87, 44227 Dortmund, Germany

E-mail address: andreas.rademacher@tu-dortmund.de (A. Rademacher).

## Abstract

In this article a semi-smooth Newton method for frictional two-body contact problems and a solution algorithm for the resulting sequence of linear systems are presented. It is based on a mixed variational formulation of the problem and a discretization by finite elements of higher-order. General friction laws depending on the normal stresses and elasto-plastic material behavior with linear isotropic hardening are considered. Numerical results show the efficiency of the presented algorithm.

**Keywords:** Higher-order FEM; Frictional contact; Semi-smooth Newton method; Complementarity function; Plasticity; Hardening

## 1. Introduction

Frictional, elasto-plastic multi-body contact problems play an important role in mechanical engineering [1–3]. The nonlinearities caused by geometric contact and frictional constraints combined with the nonlinearity in the material law result in challenging numerical problems in forms of variational inequalities and therefore efficient solving methods are needed. In the Mortar context semi-smooth Newton methods have been established as powerful methods for solving contact problems of various kinds [4]. The nonlinear contact and frictional constraints as well as the elasto-plastic material behavior can be reformulated in terms of equations defined by semi-smooth functions and handled in one monolithic semi-smooth Newton method. As shown in [5] this approach is equivalent to active set strategies. The condition number of the resulting linear systems does not increase during the iteration process. A base of these developments forms generalized Newton methods [6–8]. In contact mechanics semi-smooth Newton methods of various kinds have been developed. Active set strategies for geometrical contact and frictional constraints as well as multi-body contact are analyzed in [9–12]. The papers [13,14] show a dual mortar approach for finite deformation contact problems. Linear and quadratic finite elements have been regarded in [15]. Semi-smooth Newton methods for elasto-plastic problems are analyzed involving frictionless contact in [16,17] and frictional contact in [18–20]. Besides Mortar methods there exist a series of strategies for the numerical treatment of frictional contact problems.
One approach forms classical fixpoint methods [21]. General ideas with penalty methods and Lagrange multiplier techniques are described in [1]. A combination of both is augmented Lagrange multiplier techniques [22]. Domain decomposition methods as FETI techniques [23] or Dirichlet–Neumann algorithms in combination with Mortar discretizations [24] are efficient parallel computable solvers. By monotone multigrid constructed global convergent solvers are suggested in [25]. A cascadic multigrid algorithm for variational inequalities is presented in [26]. Under the use of higher-order DG-discretizations active set strategies have been applied to linear-elastic obstacle problems [27].

The main object of our paper is to develop an efficient framework to solve frictional multi-body contact problems using higher-order discretizations which avoid locking effects in solid mechanics for instance. Our approach has compared to Mortar methods the advantage that the higher-order basis functions are given by simple Lagrange basis functions. We transfer the approach of active-set strategies that have been developed for Mortar methods [9,10] to mixed finite elements introduced by Haslinger [28] and higher-order discretizations presented in [29,30]. Lagrange multipliers capture the geometrical contact and frictional constraints. They are discretized on a coarser mesh with ansatz functions of different polynomial degrees than the primal variable. This ansatz is conforming if d-linear or d-quadratic finite elements with $d = 2, 3$ are chosen to discretize the primal variable and becomes nonconforming for higher degrees. In [30,31] a solution scheme for linear-elastic, frictional multi-body contact problems and higher-order discretizations is suggested. It is based on the dual formulation of the discrete mixed variational formulation and leads to an optimization problem in the Lagrange multipliers. This problem is solved by a SQP method, which instantiates the contact constraints as sign- and the frictional constraints as nonlinear constraints. To include nonlinear material behavior inside a Newton iteration a contact problem has to be solved in every step to full accuracy. A more convenient ansatz in the context of nonlinear material behavior is the application of an inexact, monolithic Newton method based on active sets.

In the presented work the discrete weak formulations of the constraint inequalities are reformulated by NCP functions in terms of equations. As requirement for this proceeding we choose a convenient combination of basis functions for the higher-order Lagrange ansatz spaces and quadrature rules. Within the regarded discretization the coupling matrices are nonquadratic and nondiagonal in contrast to the Mortar ansatz. On the other hand the construction of basis functions is easier. Concerning plasticity we use a primal-mixed formulation in the displacements and project the stresses onto the admissible set [32]. The nonlinearity in the material law is processed by linearization and damping strategies. We present efficient inexact strategies to solve the arising full saddle point systems. A block triangular preconditioner is adapted from [33] and a cheap but effective preconditioning method for the Schur complement matrix is suggested. Representative for one time step of a quasi-static process a static problem is regarded. A generalization to quasi-static or dynamic multi-body problems can easily be performed [31,34].

This paper is organized as follows. In Section 2 we present the regarded frictional, elasto-plastic two-body contact problem in its strong as well as mixed variational formulation. A higher-order discretization is given in Section 3 whereas in Section 4 semi-smooth Newton methods for contact and general frictional constraints are developed for the described mixed finite elements. Numerical results of problems in two and three space dimensions substantiating the efficiency of the presented approach are discussed in Section 5.

## 2. Problem formulation

In this section we give the strong and weak formulations of the regarded frictional two-body contact problem with an elasto-plastic material law and linear isotropic hardening. We consider two deformable bodies $\Omega^m$ with $m = 1, 2$ in $d = 2$ or $3$ spatial dimensions, on which volume forces $f^m \in L^2(\Omega^m, \mathbb{R}^d)$ are acting. Their boundaries are denoted by $\Gamma^m$, $m = 1, 2$. With the outer normal vector $n$ on $\Gamma^m$ one can define the surface stresses $\sigma_n(u^m) := \sigma(u^m)n$. The normal part of these stresses is given by the scalar value $\sigma_{nn} := n^\top\sigma(u^m)n$, whereas the tangential stress vectors are calculated by $\sigma_{nt,i} := n^\top\sigma(u^m)_{ti}$. The matrix $t \in \mathbb{R}^{d \times d-1}$ contains the tangential vectors on $\Gamma_m$ that build an orthonormal system with the outer normal $n$. Define the trace operators $\gamma_M^m: H^1(\Omega^m) \to L^2(M)$ for $M \subset \Gamma^m$ and

$$H_D(\Omega^m) := \{ v \in H^1(\Omega^m) \mid \gamma_{\Gamma_D^m}(v) = 0 \}.$$

Their $d$-dimensional cartesian product space is denoted by $V := (H_D(\Omega^1))^d \times (H_D(\Omega^2))^d$. We are interested in the displacements $u = (u^1, u^2) \in V$, which fulfill the following conditions for $m = 1, 2$:

$$\varepsilon(u^m) = A^m \sigma(u^m) + \varepsilon^{m, P} \quad \text{in } \Omega^m \qquad (1)$$
$$ -\mathrm{div}\sigma(u^m) = f^m \qquad \text{in } \Omega^m \qquad (2) $$

$$ \varepsilon^{m,P}(\tau - \sigma(u^m)) \geq 0 \quad \forall \tau \text{ with } \mathcal{F}^{m,\text{iso}}(\tau, |\varepsilon^{m,P}|_F) \leq 0 \qquad \text{in } \Omega^m \qquad (3) $$

$$ u^m = 0 \qquad \text{on } \Gamma_D^m \qquad (4) $$

$$ \sigma_n(u^m) = p^m \qquad \text{on } \Gamma_N^m. \qquad (5) $$

Relation (1) describes the material law, the relation between the linearized strain $\varepsilon(u^m) = \frac{1}{2}(\nabla u^m + \nabla u^{m,\top})$ and the stress $\sigma(u^m)$. The strain is split up into an elastic part $A^m\sigma(u^m)$ with the fourth order compliance tensor $A^m$ corresponding to isotropic material and a plastic part $\varepsilon^{m,P}$. Eq. (2) ensures that outer and inner forces are balanced. The deviatoric part of a tensor $\tau$ is denoted by $\tau^D := \tau - \frac{1}{d}\mathrm{tr}(\tau)I_d \times d$ and $|\cdot|_F$ indicates the Frobenius norm. Defining the yield function $\mathcal{F}^{m,\text{iso}}(\tau, \eta) = |\tau^D|_F - (\sigma_0^m + \gamma_{\text{iso}}^m\eta)$ with the yield stress $\sigma_0^m$ and the isotropic hardening parameter $\gamma_{\text{iso}}^m$ the complementarity condition (3) ensures that plastic strain only may occur if the yield function is zero. The bodies are fixed at some closed subset of their boundary $\Gamma_D^m \subset \Gamma^m$ with positive measure. These Dirichlet boundary conditions are described in Eq. (4). At parts $\Gamma_C^m \subset \Gamma^m$ the two bodies may come into contact. We assume that the open set $\Gamma_C^m$ fulfills $\Gamma_C^m \subsetneq \Gamma^m \setminus \Gamma_D^m$. On the remaining part of the boundaries $\Gamma_N^m := \Gamma^m \setminus (\Gamma_D^m \cup \bar{\Gamma}_C^m)$ surface stresses $p^m \in L^2(\Gamma_N^m, \mathbb{R}^d)$ act on each body $\Omega^m$ by relation (5). Let $\Phi: \Gamma_C^1 \to \Gamma_C^2$ be an appropriate, bijective, sufficiently smooth mapping between the contact boundaries of the slave body $\Omega^1$ and the master body $\Omega^2$. In practical applications the mapping function $\Phi$ may not be bijective. The difficulties come up e.g. for contact boundaries that are not convex. In general the weaker condition of local convexity is sufficient to determine the gap. If this condition is not given suitable corrections have to be performed, see [35, Chapter 4.2]. In order to model contact conditions we define a generalized normal vector for $x \in \Gamma^1$ by

$$ n_{\delta}(x) := \begin{cases} \frac{\Phi(x) - x}{\|\Phi(x) - x\|}, & x \neq \Phi(x) \\ n^1(x) = -n^2(x), & x = \Phi(x) \end{cases} $$

and corresponding tangential vectors $t_{\delta} \in L^2(\Gamma^1)^{d \times (d-1)}$ such that $(n_{\delta}(x), t_{\delta}(x))$ form an orthonormal system. We define the normal jump

$$ [v]_{n_{\delta}}(x) := \gamma_{\Gamma_C^1}(v^1)(x) \cdot n_{\delta}(x) - \gamma_{\Gamma_C^2}(v^2)(\Phi(x)) \cdot n_{\delta}(x) $$

and the tangential jump

$$ [v]_{t_{\delta}}(x) := t_{\delta}(x)^{\top} \gamma_{\Gamma_C^1}(v^1)(x) - t_{\delta}(x)^{\top} \gamma_{\Gamma_C^2}(v^2)(\Phi(x)) $$

on $\Gamma_C^1$. The distance of $\Omega^1$ and $\Omega^2$ is given by the gap function $g(x) := |\Phi(x) - x|$. With this notion the contact conditions, which are defined on the slave body, read as follows:

$$ [u]_{n_{\delta}} \leq g \qquad \text{on } \Gamma_C^1 \qquad (6) $$

$$ \sigma_{n_{\delta}n_{\delta}}(u^1) \leq 0 \qquad \text{on } \Gamma_C^1 \qquad (7) $$

$$ \sigma_{n_{\delta}n_{\delta}}(u^1)([u]_{n_{\delta}} - g) = 0 \qquad \text{on } \Gamma_C^1 \qquad (8) $$

$$ \sigma_{n_{\delta}}(u^1) = -\Theta^* \sigma_{n_{\delta}}(u^2) \qquad \text{on } \Gamma_C^1. \qquad (9) $$

The bodies are not allowed to penetrate each other (6) and only negative or vanishing contact forces are allowed (7). By the complementarity condition (8) the property that either contact occurs or the normal contact forces vanish is modeled. Eq. (9) contains the adjoint $\Theta^*$ of a transfer operator $\Theta: L^2(\Gamma_C^2) \to L^2(\Gamma_C^1)$ which is defined by $\Theta(v^2)(x^1) := v(\Phi(x^1))$ and ensures equality of the contact forces on $\Gamma_C^1$ and $\Gamma_C^2$. Besides the described normal contact we regard frictional constraints on $\Gamma_C^1$:

$$ \|\sigma_{n_{\delta}t_{\delta}}(u^1)\| \leq s(\sigma_{n_{\delta}n_{\delta}}(u^1)) \qquad (10) $$

$$ \|\sigma_{n_{\delta}t_{\delta}}(u^1)\| < s(\sigma_{n_{\delta}n_{\delta}}(u^1)) \Rightarrow [u]_{t_{\delta}} = 0 \qquad (11) $$

$$ \|\sigma_{n_{\delta}t_{\delta}}(u^1)\| = s(\sigma_{n_{\delta}n_{\delta}}(u^1)) \Rightarrow \exists \alpha \in \mathbb{R}_{\ge 0}: [u]_{t_{\delta}} = \alpha \sigma_{n_{\delta}t_{\delta}}(u) \qquad (12) $$
with the euclidean norm $|| \cdot ||$. The tangential stresses are bounded by a functional $s$ which represents a general friction law depending on the normal stresses of the slave body. If this threshold is not achieved the bodies stick (11). Otherwise the slip condition (12) holds and a tangential movement occurs, that is proportional to the tangential stresses. Eqs. (1)-(12) generate a strong formulation of the regarded static, frictional, elasto-plastic two-body contact problem.

Following [32] we introduce a primal-mixed formulation of the frictional elasto-plastic two-body contact problem (1)-(12) and project the stresses onto the admissible set by the projector

$$P_{\Pi}(\tau) := \begin{cases} \tau, & |\tau^D|_F \le \sigma_0^m \\ \left( \frac{\gamma_{\text{iso}}^m}{2\mu^m + \gamma_{\text{iso}}^m} + \left( 1 - \frac{\gamma_{\text{iso}}^m}{2\mu^m + \gamma_{\text{iso}}^m} \right) \frac{\sigma_0^m}{|\tau^D|_F} \right) \tau^D + \frac{1}{d} \text{tr}(\tau) I_{d \times d}, & |\tau^D|_F > \sigma_0^m \end{cases}$$

with the shear modulus $\mu^m$ of the mth body material. This yields a variational inequality in the displacements

$$a(u)(\varphi - u) - f_{\text{ext}}(\varphi - u) + j(\varphi) - j(u) \ge 0 \quad \forall \varphi \in K$$

on the convex set

$$K := \{ v \in V | [v]_{n_\delta} \le g \} \subset V$$

containing the contact conditions and with the functional

$$j(\varphi) := \int_{\Gamma_C^1} s \|[\varphi]_{t_\delta}\| d\sigma$$

describing the frictional constraints. The semi-linear form $a$ is defined by

$$a: V \times V \to \mathbb{R}, \quad a(v)(w) := \sum_{m=1}^{2} \int_{\Omega^m} P_{\Pi}((A^m)^{-1}\varepsilon(v^m)) : \varepsilon(w^m)dx,$$

and the linear form

$$f_{\text{ext}}(v) := \sum_{m=1}^{2} (f^m, v^m)_0 + (p^m, v^m)_{0, \Gamma_N^m}$$

corresponds to the external energy with the $L^2$-scalar product $(\cdot, \cdot)_0$. Introducing the Lagrange functional,

$$\mathcal{L}(v, \mu_n, \mu_t) := \frac{1}{2} a(v)(v) - f_{\text{ext}}(v) + \langle \mu_n, [v]_{n_\delta} - g \rangle + (\mu_t, s[v]_{t_\delta})_{0, \Gamma_C^1}$$

the dual cone

$$\Lambda_n := \left( \left\{ \mu \in H^{\frac{1}{2}}(\Gamma_C^1) \mid \mu \le 0 \text{ a.e.} \right\} \right)'$$

and

$$\Lambda_t := \left\{ \mu \in L^2(\Gamma_C^1)^{d-1} \mid \|\mu\| \le s \right\}$$

we replace the contact and the frictional constraints and end up in a mixed formulation of the form: Find $(u, \lambda_n, \lambda_t) \in V \times \Lambda_n \times \Lambda_t$ such that

$$a(u)(v) + b_n(\lambda_n, v) + b_t(\lambda_t, v) = f_{\text{ext}}(v) \quad \forall v \in V \quad (13)$$

$$b_n(\mu_n - \lambda_n, u) + b_t(\mu_t - \lambda_t, u) \ge \langle \mu_n - \lambda_n, g \rangle \quad \forall \mu_n \in \Lambda_n, \forall \mu_t \in \Lambda_t. \quad (14)$$

Here $\langle \cdot, \cdot \rangle$ denotes the dual pairing on $\Lambda_n$. The bilinear forms that include the contact and frictional conditions are defined by

$$\begin{aligned} b_n: \Lambda_n \times V &\to \mathbb{R}, \\ b_t: \Lambda_t \times V &\to \mathbb{R}, \end{aligned} \quad \begin{aligned} b_n(\mu, w) &:= \langle \mu, [w]_{n_\delta} \rangle, \\ b_t(\mu, w) &:= (\mu, s[w]_{t_\delta})_{0, \Gamma_C}. \end{aligned}$$
### 3. Discretization

In this section we introduce a higher-order discretization for the described problem (13)–(14) that was mentioned in [29]. Denote by $\mathcal{T}_{h,m}$ and $\mathcal{B}_H$ triangulations of $\Omega_m$, $m = 1, 2$ and $\Gamma_C^1$, respectively, with corresponding affine transformations

$$F_{m,T} : \hat{T} := [-1, 1]^d \rightarrow T \in \mathcal{T}_{h,m}$$

and

$$F_E : \hat{E} := [-1, 1]^{d-1} \rightarrow E \in \mathcal{B}_H.$$

Let $S_l^p$ be the polynomial tensor product space of order $r$ on the reference element $[-1, 1]^l$. We define ansatz functions of polynomial degree $p \in \mathbb{N}$ on body $m = 1, 2$ by

$$V_{h,m}^p := \{ v \in (H_D(\Omega_m))^d \mid \forall T \in \mathcal{T}_{h,m} : v|_T \circ F_{m,T} \in S_d^p(T) \} := \operatorname{span}\{\phi_i^m\}_{i=1}^{n_m}$$

and the tensor space

$$V_h = V_{h,1}^{p_1} \times V_{h,2}^{p_2}.$$

The cumulated dimensions of the discrete spaces $V_{h,m}^p$ are named $n := n_1 + n_2$. By introduction of the space

$$\mathcal{M}_H^q := \{ v \in L^2(\Gamma_C^1) \mid \forall E \in \mathcal{B}_H : v|_E \circ F_E \in S_{d-1}^q(E) \}$$

with a polynomial degree $q \in \mathbb{N}$ the admissible set for the Lagrange multipliers concerning geometrical contact is defined as follows:

$$\Lambda_{n,H} := \{ v \in \mathcal{M}_H^q \mid \forall E \in \mathcal{B}_H : \forall x \in C_q : v(F_E(x)) \le 0 \} := \operatorname{span}\{\psi_i^n\}_{i=1}^{m_1}.$$

The sign condition $v(F_E(x)) \le 0$ is only defined on the finite set $C_q \subset [-1, 1]^{d-1}$ which consists of the $(q+1)^{d-1}$ Gauss-quadrature points, because for polynomials of higher degree a uniform condition is hard to satisfy. This kind of discrete sign-condition combined with numerical integration using $C_q$ yields the integral sign-condition $\int_{\Gamma_C^1} \lambda_{n,H} \psi_H \, d\lambda \ge 0$ for functions $\psi_H \in \Lambda_{n,H}$. For constant ansatz functions we define $C_0 := \{0\}_{d-1}$, whereas for $q = 1$ the set $C_1$ consists of the corners of the reference element $[-1, 1]^{d-1}$. The non-conforming ansatz was proposed in [36] and convergence is shown for an elastic two-body problem in [29]. In the case of lower polynomial degrees $q = 0, 1$, this ansatz becomes conforming, since the choice of $C_q$ leads to a uniform sign condition on the elements $E \in \mathcal{B}_H$. Accordingly the discrete space for the Lagrange multiplier concerning friction is defined by the non-conforming ansatz

$$\Lambda_{t,H} := \{ v \in (\mathcal{M}_H^q)^{d-1} \mid \forall E \in \mathcal{B}_H : \forall x \in C_q : v(F_E(x)) \le 1 \} := \operatorname{span}\{\psi_i^t\}_{i=1}^{m_2}.$$

and becomes conforming again for $q = 0, 1$. The discrete space $\Lambda_{t,H}$ is spanned by the $m_2 := (d-1)m_1$ basis vectors $\{\psi_j^n e_k, j = 1, \dots, m_1, k = 1, \dots, d-1\}$ with the $k$th unit vector $e_k$ in $\mathbb{R}^{d-1}$. In two space dimensions this leads to coinciding bases for $\Lambda_{n,H}$ and $\Lambda_{t,H}$.

**Remark 1.** The advantage of this choice is the given split of the contact pressure $\sigma_n(x)$ on $\Gamma_C^1$ in its normal and tangential part $\sigma_n(x) = \sigma_{nn}(x)n(x) + \sigma_{nt}(x)t(x)$. The normal and tangential components of the contact pressure are directly calculated by $\sigma_{nn}(x) = \lambda_n(x)n(x)$ and $\sigma_{nt}(x) = \sum_{i=1}^{d-1} \lambda_{t,i}(x)t_i(x)$. Clearly this ansatz does not need to represent constant tangential forces by a constant Lagrange multiplier. Furthermore the tangential system has to be chosen in correspondence to the discretization.

By this definition the solution depends on the choice of tangential vectors. However the splitting of the normal and tangential parts is clear. Stability of the described discretization is proven for the elastic case and balanced $h, H, p$, and $q$ in [29]. It is shown that if the inequality

$$\sum_{m=1}^{2} (hH^{-1} \max\{1, q\}^2 p^{-1})^{\theta_m} \le \varepsilon$$
holds for an $\varepsilon > 0$ sufficiently small and $1 + \theta_m$-regularity of the solution $u$ there exists an $\alpha \in \mathbb{R}_{\ge 0}$ independent of $h, H, p$ and $q$ such that the Brezzi-Babuska condition

$$ \alpha \|(\mu_{n,H}, \mu_{t,H})\|_{-1/2} \le \sup_{v_h \in V_h, \|v_h\|=1} (b_n(\mu_{n,H}, v_h) + b_t(\mu_{t,H}, v_h)) $$

is fulfilled for all $(\mu_{n,H}, \mu_{t,H}) \in \Lambda_{n,H} \times \Lambda_{t,H}$. Practical investigations show that the choice $q = \max\{p^m\} - 1$ and $H = 2 \max\{h^m\}$ leads to a stable discretization, cf. [31]. Eventually the discrete problem is to find $(u_h, \lambda_{n,H}, \lambda_{t,H}) \in V_h \times \Lambda_{n,H} \times \Lambda_{t,H}$ with

$$ \begin{aligned} & a(u_h)(v_h) + b_n(\lambda_{n,H}, v_h) + b_t(\lambda_{t,H}, v_h) = f_{\text{ext}}(v_h) \\ & b_n(\mu_{n,H} - \lambda_{n,H}, u_h) + b_t(\mu_{t,H} - \lambda_{t,H}, u_h) \ge \langle \mu_{n,H} - \lambda_{n,H}, g \rangle \end{aligned} \quad (15) $$

for all $v_h \in V_h$, $\mu_{n,H} \in \Lambda_{n,H}$ and $\mu_{t,H} \in \Lambda_{t,H}$.

# 4. Semi-smooth Newton methods

Following the work [9] we transfer the presented active-set strategies for Mortar methods to mixed finite elements, which were introduced by Haslinger [28] for low order. We generalize the solving methods to higher-order discretizations that are proposed in [29] and general friction laws depending on the contact forces.

## 4.1. Active-set strategy for contact

We begin with presenting a semi-smooth Newton method for geometrical contact constraints. For simplicity we skip the subscript $\delta$ of the generalized normals $n_\delta$ in the two-body contact formulations. The point of departure is formed by the weak version of the contact conditions (6)–(8)

$$ \int_{\Gamma_C^1} ([u_h]_n - g) \psi_H \, do \le 0 \quad (16) $$

$$ \int_{\Gamma_C^1} \lambda_{n,H} \psi_H \, do \ge 0 \quad (17) $$

$$ \left( \int_{\Gamma_C^1} \lambda_{n,H} \psi_H \, do \right) \left( \int_{\Gamma_C^1} ([u_h]_n - g) \psi_H \, do \right) = 0 \quad (18) $$

for all $\psi_H \in \Lambda_{n,H}$. The weak contact conditions (16)–(18) correspond to the discrete mixed formulation (15) if the numerical quadrature rule based on the gaussian quadrature points $C_q$ is applied. Biorthogonal ansatz functions and suitable lumping techniques would also lead to decoupled conditions. The first Eq. (16) decouples by biorthogonality whereas the inequalities (17) and (18) end up in nodewise conditions by lumping. Define the coupling matrix $N = [N_1 \ N_2]$ of size $N_m \in \mathbb{R}^{m_1 \times n_m}$ concerning the bases of $V_h$ and $\Lambda_{H,n}$ by

$$ (N_1)_{ij} := \int_{\Gamma_C^1} \psi_i^n(x) \gamma_{\Gamma_C^1} (\phi_j^1)(x) \cdot n \, do, $$

$$ (N_2)_{ij} := \int_{\Gamma_C^1} \psi_i^n(x) \gamma_{\Gamma_C^2} (\phi_j^2)(\Phi(x)) \cdot n \, do $$

as well as the mass matrix

$$ M \in \mathbb{R}^{m_1 \times m_1}, \quad M_{ij} := \int_{\Gamma_C^1} \psi_j^n \psi_i^n \, do $$

and the gap vector

$$ \bar{g} \in \mathbb{R}^{m_1}, \quad \bar{g}_i = \int_{\Gamma_C^1} g \psi_i^n \, do. $$
In the following a bar generally indicates the vector-valued representation of the corresponding discrete function in the belonging basis. The vector-valued formulation of (16)–(18) reads

$$N\bar{u}_h - \bar{g} \le 0, \quad M\bar{\lambda}_{n,H} \le 0, \quad (M\bar{\lambda}_{n,H})_i (N\bar{u}_h - \bar{g})_i = 0 \quad \forall i = 1, \dots, m_1. \tag{19}$$

To reformulate (19) in terms of an equation we define the NCP function

$$C_N : \mathbb{R}^n \times \mathbb{R}^{m_1} \to \mathbb{R}^{m_1} \\ C_N(\bar{u}_h, \bar{\lambda}_{n,H})_i := (-M\bar{\lambda}_{n,H})_i - \max\{0, (-M\bar{\lambda}_{n,H})_i + c_n(N\bar{u}_h - \bar{g})_i\}$$

with a positive constant $c_n$. Following the arguments in [9, Chapter 4] the discrete, weak contact conditions (19) can equivalently be expressed by

$$C_N(\bar{u}_h, \bar{\lambda}_{n,H}) = 0. \tag{20}$$

With variations $\delta\bar{u}_h^k, \delta\bar{\lambda}_{n,H}^k$ and characteristical functions $\chi_i$ defined by

$$\chi_i := \begin{cases} 1, & (-M\bar{\lambda}_{n,H} + c_n(N\bar{u}_h - \bar{g}))_i > 0 \\ 0, & (-M\bar{\lambda}_{n,H} + c_n(N\bar{u}_h - \bar{g}))_i \le 0 \end{cases}$$

the generalized derivative of the NCP function reads

$$C'_N(\bar{u}_h, \bar{\lambda}_{n,H})(\delta\bar{u}_h, \delta\bar{\lambda}_{n,H})_i = -\chi_i (-M\delta\bar{\lambda}_{n,H} + c_n N\delta\bar{u}_h)_i - (M\delta\bar{\lambda}_{n,H})_i.$$

We apply a semi-smooth Newton's method

$$C'_N(\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1})(\delta\bar{u}_h^k, \delta\bar{\lambda}_{n,H}^k) = -C_N(\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1})$$

to solve (20). The new iterates are calculated by $\bar{u}_h^k := \bar{u}_h^{k-1} + \delta\bar{u}_h^k$ and $\bar{\lambda}_{n,H}^k := \bar{\lambda}_{n,H}^{k-1} + \delta\bar{\lambda}_{n,H}^k$. Let active and inactive indices in Newton step $k$ be determined in the following way

$$A_n^k := \{i \in \{1, \dots, m_1\} \mid (-M\bar{\lambda}_{n,H} + c_n(N\bar{u}_h - \bar{g}))_i > 0\} \tag{21}$$

$$\mathcal{I}_n^k := \{i \in \{1, \dots, m_1\} \mid (-M\bar{\lambda}_{n,H} + c_n(N\bar{u}_h - \bar{g}))_i \le 0\}. \tag{22}$$

Then the iterate $(\bar{u}_h^k, \bar{\lambda}_{n,H}^k)$ of Newton's method w.r.t. the contact conditions can equivalently be expressed by

$$(N\bar{u}_h^k)_i = \bar{g}_i, \quad i \in A_n^k$$

$$(M\bar{\lambda}_{n,H}^k)_i = 0, \quad i \in \mathcal{I}_n^k$$

for all $i = 1, \dots, m_1$.

## 4.2. Active-set strategy for friction

In this subsection we present an active-set strategy for frictional constraints in the context of higher-order finite elements. The discrete analog of the frictional constraints (10)–(12) reads

$$\begin{align}
\|\sigma_{nt}(u_h^1)\| &\le s(\sigma_{nn}(u_h^1)) \nonumber \\
\|\sigma_{nt}(u_h^1)\| &< s(\sigma_{nn}(u_h^1)) \Rightarrow [u_h]_t = 0 \tag{23} \\
\|\sigma_{nt}(u_h^1)\| &= s(\sigma_{nn}(u_h^1)) \Rightarrow \exists \alpha \in \mathbb{R}_{\ge 0} : [u_h]_t = \alpha\sigma_{nt}(u_h). \nonumber
\end{align}$$
Instead of these pointwise conditions (23) we use a weak formulation with test functions $\psi_n \in \Lambda_n$ respectively $\psi_t \in \Lambda_t$:

$$
\begin{align*}
\int_{\Gamma_C^1} \|\lambda_{t,H}\| \psi_n \, do &\le \int_{\Gamma_C^1} s(\lambda_{n,H}) \psi_n \, do \\
\int_{\Gamma_C^1} \|\lambda_{t,H}\| \psi_n \, do &< \int_{\Gamma_C^1} s(\lambda_{n,H}) \psi_n \, do \Rightarrow \int_{\Gamma_C^1} [u_h]_t \psi_t \, do = 0 \\
\int_{\Gamma_C^1} \|\lambda_{t,H}\| \psi_n \, do &= \int_{\Gamma_C^1} s(\lambda_{n,H}) \psi_n \, do \Rightarrow \exists \alpha \in \mathbb{R}_{\ge 0}: \int_{\Gamma_C^1} [u_h]_t \psi_t \, do = \alpha \int_{\Gamma_C^1} \lambda_{t,H} \psi_t \, do.
\end{align*}
$$

We replace the normal and tangential stresses $\sigma_{nn}(u_h)$ and $\sigma_{nt}(u_h)$ by the corresponding Lagrange multipliers $\lambda_n$ and $\lambda_t$. The integrals $\int_{\Gamma_C^1} \|\lambda_{t,H}\| \psi_i^n \, do$ and $\int_{\Gamma_C^1} \lambda_t \psi_i^t \, do$ are approximated respectively integrated exactly by a $(q+1)^{d-1}$-point gaussian quadrature with weights $\alpha_i$. The nodes of the Lagrangian basis of $\Lambda_{n,H}$ and $\Lambda_{t,H}$ are placed in these gaussian quadrature points $C_q := \{\hat{x}_1, \dots, \hat{x}_{(q+1)^{d-1}}\}$ on the reference element $[-1, 1]^{d-1}$. Define the matrix $\bar{T} \in \mathbb{R}^{m_2 \times n}$ with partition $\bar{T} = [\bar{T}_1 \bar{T}_2]$ by

$$
\begin{align*}
(\bar{T}_1)_{ij} &:= \frac{1}{\omega_i} \int_{\Gamma_C^1} \psi_i^t(x) t(x)^{\top} \gamma_{\Gamma_C^1} (\phi_j^1)(x) \, do \\
(\bar{T}_2)_{ij} &:= \frac{1}{\omega_i} \int_{\Gamma_C^1} \psi_i^t(x) t(x)^{\top} \gamma_{\Gamma_C^2} (\phi_j^2)(\Phi(x)) \, do.
\end{align*}
$$

The weight $\omega_i$ is given by $\omega_i := \alpha_i |\det \nabla F_E(\hat{x}_i)^{\top} \nabla F_E(\hat{x}_i)| \ge 0$ on $E = \text{supp}(\psi_i^t)$. Let

$$
s_i := \frac{1}{\omega_i} \int_{\Gamma_C^1} s(\lambda_{n,H}^-) \psi_i^n \, do \quad (24)
$$

with

$$
\lambda_{n,H}^{-} := \min\{0, \lambda_{n,H}\}
$$

be the algebraic representation of the positive friction bound $s$. The resulting discrete frictional constraints read

$$
\begin{align}
\| \bar{\lambda}_{t,H,\mathcal{I}(i)} \| &\le s_i \notag \\
\| \bar{\lambda}_{t,H,\mathcal{I}(i)} \| &< s_i \Rightarrow (\bar{T} \bar{u}_h)_{\mathcal{I}(i)} = 0 \tag{25} \\
\| \bar{\lambda}_{t,H,\mathcal{I}(i)} \| &= s_i \Rightarrow \exists \alpha \in \mathbb{R}_{\ge 0} : (\bar{T} \bar{u}_h)_{\mathcal{I}(i)} = \alpha \bar{\lambda}_{t,H,\mathcal{I}(i)} \notag
\end{align}
$$

for the dimension-dependent index map

$$
\mathcal{I}(i) := \begin{cases} \{i\}, & d=2 \\ \{2i-1, 2i\}, & d=3. \end{cases}
$$

In this sense for a vector $u \in \mathbb{R}^{m_2}$ the selection $u_{\mathcal{I}(i)}$ is a scalar value in the case of two space dimensions and a vector in $\mathbb{R}^2$ if $d=3$. Following [9, Section 5] we express the conditions (25) by the NCP function

$$
C_T : \mathbb{R}^n \times \mathbb{R}_{\le 0}^{m_1} \times \mathbb{R}^{m_2} \to \mathbb{R}^{m_2}, \\
C_T(\bar{u}_h, \bar{\lambda}_{n,H}, \bar{\lambda}_{t,H})_{\mathcal{I}(i)} := \max \{s_i, \|(\bar{\lambda}_{t,H,\mathcal{I}(i)} + c_t \bar{T}(\bar{u}_h)_{\mathcal{I}(i)})\| \} \bar{\lambda}_{t,H,\mathcal{I}(i)} - s_i (\bar{\lambda}_{t,H,\mathcal{I}(i)} + c_t \bar{T}(\bar{u}_h)_{\mathcal{I}(i)})
$$

with $c_t > 0$ and $i \in \{1, \dots, m_1\}$. Analog arguments to the proof in [9, Theorem 5.1] lead to the equivalence of the equation $C_T(\bar{u}_h, \bar{\lambda}_{n,H}, \bar{\lambda}_{t,H}) = 0$ and the frictional constraints (25). Defining the characteristical functions

$$
\chi_{A_{t,i}} := \begin{cases} 1, & \| \bar{\lambda}_{t,H,\mathcal{I}(i)} + c_t \bar{T} \bar{u}_{h,\mathcal{I}(i)} \| - s_i > 0 \\ 0, & \| \bar{\lambda}_{t,H,\mathcal{I}(i)} + c_t \bar{T} \bar{u}_{h,\mathcal{I}(i)} \| - s_i \le 0 \end{cases}
$$
and

$$ \chi_{\mathcal{I}_{t,i}} := \begin{cases} 1, & \| \bar{\lambda}_{t,H,\mathcal{I}(i)} + c_t \bar{T} \bar{u}_h \mathcal{I}(i) \| - s_i \le 0 \\ 0, & \| \bar{\lambda}_{t,H,\mathcal{I}(i)} + c_t \bar{T} \bar{u}_h \mathcal{I}(i) \| - s_i > 0 \end{cases} $$

the generalized derivative for $\bar{\lambda}_{n,H,i} < 0$ is given by

$$ \begin{aligned} C'_T(\bar{u}_h, \bar{\lambda}_{n,H}, \bar{\lambda}_{t,H}) (\delta\bar{u}_h, \delta\bar{\lambda}_{n,H}, \delta\bar{\lambda}_{t,H})_{\mathcal{I}(i)} & = \chi_{\mathcal{A}_{t,i}} \frac{\bar{\lambda}_{t,H,\mathcal{I}(i)} (\bar{\lambda}_{t,H,\mathcal{I}(i)} + c_t(\bar{T}\bar{u}_h)_{\mathcal{I}(i)})^\top (\delta\bar{\lambda}_{t,H,\mathcal{I}(i)} + c_t(\bar{T}\delta\bar{u}_h)_{\mathcal{I}(i)})}{\|\bar{\lambda}_{t,H,\mathcal{I}(i)} + c_t(\bar{T}\bar{u}_h)_{\mathcal{I}(i)}\|} \\ & \quad + \delta\bar{\lambda}_{t,\mathcal{I}(i)} \max\{s_i, \|(\bar{\lambda}_{t,H,\mathcal{I}(i)} + c_t(\bar{T}\bar{u}_h)_{\mathcal{I}(i)})\| - s_i(\delta\bar{\lambda}_{t,\mathcal{I}(i)} + c_t(\bar{T}\delta\bar{u}_h)_{\mathcal{I}(i)}) \\ & \quad - s'_{i,\lambda_{n,H}}(\delta\bar{\lambda}_{n,H})(\bar{\lambda}_{t,H,\mathcal{I}(i)} + c_t(\bar{T}\bar{u}_h)_{\mathcal{I}(i)}) + \chi_{\mathcal{I}_{t,i}} s'_{i,\lambda_{n,H}}(\delta\bar{\lambda}_{n,H}) \bar{\lambda}_{t,\mathcal{I}(i)}. \end{aligned} $$

We formulate the semi-smooth Newton method

$$ C'_{T}(\bar{u}_{h}^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1})(\delta\bar{u}_{h}, \delta\bar{\lambda}_{n,H}, \delta\bar{\lambda}_{t,H}) = -C_{T}(\bar{u}_{h}^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1}) \quad (26) $$

with increments $\delta\bar{u}_h^k = \bar{u}_h^k - \bar{u}_h^{k-1}$, $\delta\bar{\lambda}_{n,H}^k = \bar{\lambda}_{n,H}^k - \bar{\lambda}_{n,H}^{k-1}$ and $\delta\bar{\lambda}_{t,H}^k = \bar{\lambda}_{t,H}^k - \bar{\lambda}_{t,H}^{k-1}$. The Newton equation (26) is equivalent to

$$ \int_{\Gamma_C^1} \mathcal{W}(\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1})_i \bar{\lambda}_{t,H,\mathcal{I}(i)}^k \psi \, do + \int_{\Gamma_C^1} \mathcal{V}(\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1})_i \bar{\lambda}_{n,H,i}^k \psi \, do \\ + \int_{\Gamma_C^1} \mathcal{W}(\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1})_i (\bar{T}\bar{u}_h^k)_{\mathcal{I}(i)} \psi \, do = \int_{\Gamma_C^1} \mathcal{R}(\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1})_i \psi \, do \quad (27) $$

for all $i \in \{1, \dots, m_1\}$ and $\psi \in \Lambda_{n,H}$. We define the dimension-dependent index map $\mathcal{J}(j) = \lceil \frac{j}{d-1} \rceil$ and separate the set of indices $\{1, \dots, m_2\}$ into slip indices with contact

$$ \mathcal{A}_t^k := \{ j \in \{1, \dots, m_2\} \mid i := \mathcal{J}(j) : \| \bar{\lambda}_{t,H,\mathcal{I}(i)}^{k-1} + \bar{T}(\bar{u}_h^{k-1})_{\mathcal{I}(i)} \| - s_i^{k-1} > 0, \bar{\lambda}_{n,H,i}^{k-1} < 0 \}, \quad (28) $$

stick indices in contact

$$ \mathcal{I}_t^k := \{ j \in \{1, \dots, m_2\} \mid i := \mathcal{J}(j) : \| \bar{\lambda}_{t,H,\mathcal{I}(i)}^{k-1} + \bar{T}(\bar{u}_h^{k-1})_{\mathcal{I}(i)} \| - s_i^{k-1} \le 0, \bar{\lambda}_{n,H,i}^{k-1} < 0 \}, \quad (29) $$

and indices without contact

$$ \mathcal{I}_{tn}^k := \{ j \in \{1, \dots, m_2\} \mid \bar{\lambda}_{n,H,\mathcal{J}(j)}^{k-1} = 0 \}. \quad (30) $$

For $j \in \mathcal{A}_t^k$ and $i := \mathcal{J}(j)$ we obtain the functions

$$ \mathcal{W}(\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1})_i = \frac{\bar{\lambda}_{t,H,\mathcal{I}(i)}^{k-1} (\bar{\lambda}_{t,H,\mathcal{I}(i)}^{k-1} + c_t(\bar{T}\bar{u}_h^{k-1})_{\mathcal{I}(i)})^\top}{\| (\bar{\lambda}_{t,H,\mathcal{I}(i)}^{k-1} + c_t(\bar{T}\bar{u}_h^{k-1})_{\mathcal{I}(i)}) \|} \\ + I_{(d-1) \times (d-1)} \left( \| (\bar{\lambda}_{t,H,\mathcal{I}(i)}^{k-1} + c_t(\bar{T}\bar{u}_h^{k-1})_{\mathcal{I}(i)}) \| - s_i^{k-1} \right), $$

$$ \mathcal{V}(\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1})_i = -(s'_{i,\lambda_{n,H}})^{k-1} (\bar{\lambda}_{t,H,\mathcal{I}(i)}^{k-1} + c\bar{T}(\bar{u}_h^{k-1})_{\mathcal{I}(i)}), $$

$$ \mathcal{W}(\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1})_i = c_t \frac{\bar{\lambda}_{t,H,\mathcal{I}(i)}^{k-1} (\bar{\lambda}_{t,H,\mathcal{I}(i)}^{k-1} + c_t(\bar{T}\bar{u}_h^{k-1})_{\mathcal{I}(i)})^\top}{\| \bar{\lambda}_{t,H,\mathcal{I}(i)}^{k-1} + c_t(\bar{T}\bar{u}_h^{k-1})_{\mathcal{I}(i)} \|} - c_t I_{(d-1) \times (d-1)} s_i^{k-1} $$
and

$$
\begin{equation}
\begin{split}
\mathcal{R}(\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1})_i = c_t \frac{\bar{\lambda}_{t,H,\mathcal{I}(i)}^{k-1} (\bar{\lambda}_{t,H,\mathcal{I}(i)}^{k-1} + c_t(\bar{T}\bar{u}_h^{k-1})_{\mathcal{I}(i)})^{\top}}{\|\bar{\lambda}_{t,H,\mathcal{I}(i)}^{k-1} + c_t(\bar{T}\bar{u}_h^{k-1})_{\mathcal{I}(i)}\|} \cdot & \\
& (\bar{\lambda}_{t,H,\mathcal{I}(i)}^{k-1} + c_t(\bar{T}\bar{u}_h^{k-1})_{\mathcal{I}(i)}) \\
& - (s'_{i,\lambda_{n,H}})^{k-1} (\bar{\lambda}_{n,i}^{k-1}) (\bar{\lambda}_{t,H,\mathcal{I}(i)}^{k-1} + c_t(\bar{T}\bar{u}_h^{k-1})_{\mathcal{I}(i)}).
\end{split}
\end{equation}
$$

In the slip case $j \in \mathcal{I}_i^k$ the functions $\mathcal{W}_i = 0$,

$$
\mathcal{W}(\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1})_i = \frac{1}{\omega_i} I_{(d-1)\times(d-1)},
$$

$$
\mathcal{V}(\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1})_i = (\bar{T}\bar{u}^{k-1})_{\mathcal{I}(i)} \frac{(s'_{i,\lambda_{n,H}})^{k-1}}{s_i^{k-1}}
$$

and

$$
\mathcal{R}(\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1})_i = (\bar{T}\bar{u}^{k-1})_{\mathcal{I}(i)} \frac{(s'_{i,\lambda_{n,H}})^{k-1} \bar{\lambda}_{n,H,i}^{k-1}}{s_i^{k-1}}
$$

determine the Newton equation (26). In the limit case $j \in \mathcal{I}_n^k$ of no contact we choose $\bar{\lambda}_{t,H,\mathcal{I}(i)} = 0$ and therefore vanishing $\mathcal{W}_i$, $\mathcal{V}_i$ and $\mathcal{R}_i$ as well as $\mathcal{W}_i = I_{(d-1)\times(d-1)}$. We define the matrices $U = [U_1 \ U_2] \in \mathbb{R}^{m_2 \times n}$, $V \in \mathbb{R}^{m_1 \times m_2}$, $W \in \mathbb{R}^{m_1 \times m_2}$ and the vector $\bar{r} \in \mathbb{R}^{m_2}$ by

$$
(U_1)_{ji} := \int_{\Gamma_C^1} \mathcal{W}(\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1}) \mathcal{J}_{(j)} t_\delta(x)^\top \gamma_{\Gamma_C^1}(\phi_i^1)(x) \psi_j^t(x) d\sigma
$$

$$
(U_2)_{ji} := \int_{\Gamma_C^1} \mathcal{W}(\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1}) \mathcal{J}_{(j)} t_\delta(x)^\top \gamma_{\Gamma_C^1}(\phi_i^2)(x) \psi_j^t(\Phi(x)) d\sigma
$$

$$
V_{ji} := \int_{\Gamma_C^1} \mathcal{V}(\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1}) \mathcal{J}_{(j)} \psi_i^n \psi_j^t d\sigma
$$

$$
W_{ji} := \int_{\Gamma_C^1} \mathcal{W}(\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1}) \mathcal{J}_{(j)} \psi_i^t \psi_j^t d\sigma
$$

$$
\bar{r}_j := \int_{\Gamma_C^1} \mathcal{R}(\bar{T}\bar{u}_h^{k-1}, \bar{\lambda}_{n,H}^{k-1}, \bar{\lambda}_{t,H}^{k-1}) \mathcal{J}_{(j)} \psi_j^t d\sigma.
$$

Then the Newton equation (27) is equivalent to the vector-valued equation

$$
U \bar{u}_h^k + V \bar{\lambda}_{n,H}^k + W \bar{\lambda}_{t,H}^k = \bar{r}.
$$

4.3. Algebraic representation of the saddle-point system

In this subsection we give the linear systems that have to be solved during the Newton iteration. Furthermore we present the resulting algorithm and the used solving techniques. According to $\bar{T}$ let the unscaled matrix $T = [T_1 \ T_2] \in \mathbb{R}^{m_2 \times n}$ be defined by

$$
(T_1)_{ij} := \int_{\Gamma_C^1} \psi_i^t(x) t(x)^\top \gamma_{\Gamma_C^1}(\phi_j^1)(x) \, do,
$$

$$
(T_2)_{ij} := \int_{\Gamma_C^1} \psi_i^t(x) t(x)^\top \gamma_{\Gamma_C^2}(\phi_j^2)(\Phi(x)) \, do,
$$

as well as the block-diagonal matrix

$$
K^k = \operatorname{diag}(K_1^k, K_2^k) = \left( a'(u_h^{k-1})(\phi_j^m, \phi_i^m) \right)_{i,j=1,...,n}
$$
by the Frechét derivative of the semi-linearform in direction $u_h^{k-1}$. Define the vector

$$L^k(u_h^{k-1}) = [L_1^k L_2^k]^\top = (a(u_h^{k-1})(\phi_i^m))_{i=1,...,n}$$

and the Newton right hand side of the plastic problem

$$\tilde{f}_m^k := K_m^k \bar{u}_{h,m}^{k-1} - L_m^k(u_h^{k-1}) + \bar{f}_m.$$

With this linearization of the plasticity and the described introduction of active-sets for contact and friction we approximate the solution of the discrete mixed formulation (15) by a semi-smooth Newton method in which step k corresponds to solving the linear system

$$\begin{bmatrix} K_1^k & 0 \\ 0 & K_2^k \end{bmatrix} \begin{bmatrix} N_{1,\mathcal{A}_n^k}^\top & N_{1,\mathcal{I}_n^k}^\top \\ N_{2,\mathcal{A}_n^k}^\top & N_{2,\mathcal{I}_n^k}^\top \end{bmatrix} \begin{bmatrix} T_{1,\mathcal{A}_i^k}^\top & T_{1,\mathcal{I}_i^k}^\top & T_{1,\mathcal{I}_{in}^k}^\top \\ T_{2,\mathcal{A}_i^k}^\top & T_{2,\mathcal{I}_i^k}^\top & T_{2,\mathcal{I}_{in}^k}^\top \end{bmatrix} \begin{bmatrix} \bar{u}_{h,1}^k \\ \bar{u}_{h,2}^k \\ \bar{\lambda}_{n,H,\mathcal{A}_n^k}^k \\ \bar{\lambda}_{t,H,\mathcal{A}_t^k}^k \\ \bar{\lambda}_{t,H,\mathcal{I}_t^k}^k \\ \bar{\lambda}_{t,H,\mathcal{I}_{in}^k}^k \end{bmatrix} = \begin{bmatrix} \tilde{f}_1 \\ \tilde{f}_2 \\ \bar{g}_{\mathcal{A}_n^k} \\ 0 \\ \bar{r}_{\mathcal{A}_t^k} \\ \bar{r}_{\mathcal{I}_t^k} \\ 0 \end{bmatrix} \quad (31)$$

with the following notation. For a matrix $A \in \mathbb{R}^{k \times l}$ the submatrix

$$A_{IJ} := [A_{i,j}]_{i \in \mathcal{I}, j \in \mathcal{J}} \in \mathbb{R}^{|\mathcal{I}| \times |\mathcal{J}|}$$

consists of the rows and columns of A, whose indices belong to the sets $\mathcal{I}$ respectively $\mathcal{J}$ of indices. If only rows of $\mathcal{I}$ are selected the notation is

$$A_{\mathcal{I}} := [A_{i,\cdot}]_{i \in \mathcal{I}} \in \mathbb{R}^{|\mathcal{I}| \times l}.$$

The linear system corresponds to a saddle point problem of the structure

$$\begin{bmatrix} K^k & B_1^\top \\ B_2^k & -C^k \end{bmatrix} \begin{bmatrix} \bar{u}_h^k \\ \bar{\lambda}^k \end{bmatrix} = b \quad (32)$$

with partitions

$$\bar{\lambda}^k := \begin{bmatrix} \bar{\lambda}_{n,H}^k \\ \bar{\lambda}_{t,H}^k \end{bmatrix}, \quad B_1^\top := \begin{bmatrix} N_1^\top & T_1^\top \\ N_2^\top & T_2^\top \end{bmatrix}, \quad B_2^k := \begin{bmatrix} N_{1,\mathcal{A}_n^k} & N_{2,\mathcal{A}_n^k} \\ 0 & 0 \\ U_{1,\mathcal{A}_i^k} & U_{2,\mathcal{A}_i^k} \\ \bar{T}_{1,\mathcal{I}_i^k} & \bar{T}_{2,\mathcal{I}_i^k} \\ 0 & 0 \end{bmatrix}$$
and

$$ -C^k := \begin{bmatrix} 0 & 0 & 0 & 0 & 0 \\ M_{A_n^k \mathcal{I}_n^k} & M_{\mathcal{I}_n^k \mathcal{I}_n^k} & 0 & 0 & 0 \\ V_{A_t^k A_n^k} & V_{A_t^k \mathcal{I}_n^k} & W_{A_t^k A_t^k} & W_{A_t^k \mathcal{I}_t^k} & W_{A_t^k \mathcal{I}_{tn}^k} \\ V_{\mathcal{I}_t^k A_n^k} & V_{\mathcal{I}_t^k \mathcal{I}_n^k} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & I \end{bmatrix}, \quad b := \begin{bmatrix} \tilde{f}_1 \\ \tilde{f}_2 \\ \hline \bar{s}_{A_n^k} \\ 0 \\ \hline \bar{r}_{A_t^k} \\ \bar{r}_{\mathcal{I}_t^k} \\ 0 \end{bmatrix} $$

In general the system matrix of the problem (32) is asymmetric. The matrices $B_1$ and $B_2$ are nonquadratic whereas the quadratic matrix $C^k$ might be singular. To solve the discrete contact problem (15) we have to solve a sequence of saddle point problems (32). This semi-smooth Newton method is summarized in Algorithm 1.

<table><caption>**Algorithm 1: Semi-smooth Newton method**</caption>
  <thead>
    <tr>
      <th></th>
      <th>for k = 0 to maxIter do</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>1</td>
      <td>assemble linearization K<sup>k</sup> and right hand side f&#773;<sup>k</sup></td>
    </tr>
    <tr>
      <td>2</td>
      <td>calculate active and inactive sets A<sub>n</sub><sup>k</sup>, I<sub>n</sub><sup>k</sup> for contact (21, 22)</td>
    </tr>
    <tr>
      <td>3</td>
      <td>calculate friction bound (24)</td>
    </tr>
    <tr>
      <td>4</td>
      <td>calculate active and inactive sets A<sub>t</sub><sup>k</sup>, I<sub>t</sub><sup>k</sup>, I<sub>tn</sub><sup>k</sup> for friction (28, 29, 30)</td>
    </tr>
    <tr>
      <td>5</td>
      <td>assemble Newton matrices U<sup>k</sup>, V<sup>k</sup>, W<sup>k</sup> and right hand side r&#773;<sup>k</sup></td>
    </tr>
    <tr>
      <td>6</td>
      <td>solve linear system (32) to tolerance tol<sub>k</sub></td>
    </tr>
    <tr>
      <td>7</td>
      <td>calculate residuals</td>
    </tr>
    <tr>
      <td></td>
      <td>res<sub>plast</sub> := &sum;<sub>m=1</sub><sup>2</sup> ||K<sub>m</sub><sup>k</sup>u&#773;<sub>h,m</sub> - L<sub>m</sub><sup>k</sup>(u&#773;<sub>h,m</sub>) + f&#773;<sup>m</sup>||</td>
    </tr>
    <tr>
      <td></td>
      <td>C<sub>N</sub><sup>k</sup> := C<sub>N</sub>(u&#773;<sub>h</sub><sup>k</sup>, &lambda;&#773;<sub>n,H</sub><sup>k</sup>)</td>
    </tr>
    <tr>
      <td></td>
      <td>C<sub>T</sub><sup>k</sup> := C<sub>T</sub>(u&#773;<sub>h</sub><sup>k</sup>, &lambda;&#773;<sub>n,H</sub><sup>k</sup>, &lambda;&#773;<sub>t,H</sub><sup>k</sup>)</td>
    </tr>
    <tr>
      <td>8</td>
      <td>if ||C<sub>N</sub><sup>k</sup>|| + ||C<sub>T</sub><sup>k</sup>|| + res<sub>plast</sub> &le; tol then</td>
    </tr>
    <tr>
      <td></td>
      <td>&lfloor; stop &rfloor;</td>
    </tr>
    <tr>
      <td>9</td>
      <td>perform damping if k &gt; 2:</td>
    </tr>
    <tr>
      <td></td>
      <td>u&#773;<sub>m</sub><sup>0</sup> := u&#773;<sub>h,m</sub><sup>k</sup>, res<sub>0</sub> := res<sub>plast</sub>, &omega;<sub>0</sub> := 0.5</td>
    </tr>
    <tr>
      <td></td>
      <td>for j = 0 to n<sub>damp</sub> do</td>
    </tr>
    <tr>
      <td></td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;u&#773;<sub>m</sub><sup>j</sup> := &omega;<sub>j</sub>u&#773;<sub>h,m</sub><sup>k</sup> + (1 - &omega;<sub>j</sub>)u&#773;<sub>h,m</sub><sup>k-1</sup></td>
    </tr>
    <tr>
      <td></td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;res<sub>j</sub> := &sum;<sub>m=1</sub><sup>2</sup> ||K<sub>m</sub><sup>k</sup>u&#773;<sub>m</sub><sup>j</sup> - L<sub>m</sub><sup>k</sup>(u&#773;<sub>m</sub><sup>j</sup>) + f&#773;<sup>m</sup>||</td>
    </tr>
    <tr>
      <td></td>
      <td>if res<sub>j</sub> &le; res<sub>j-1</sub> then</td>
    </tr>
    <tr>
      <td></td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;w<sub>j</sub> := 0.5 w<sub>j-1</sub></td>
    </tr>
    <tr>
      <td></td>
      <td>else</td>
    </tr>
    <tr>
      <td></td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;u&#773;<sub>h</sub><sup>k</sup> := u&#773;<sub>m</sub><sup>j-1</sup></td>
    </tr>
    <tr>
      <td></td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;break</td>
    </tr>
  </tbody>
</table>

In every Newton step the linear system (32) is solved to a tolerance tol<sub>k</sub>. We apply an inexact strategy and start with a relatively high value of tol<sub>k</sub> and successively reduce it. Because the system matrix is asymmetric we use a GMRES
solver. Furthermore we specify a block triangular preconditioner of the form

$$ \mathcal{P}^k = \begin{bmatrix} K^k & B_1^\top \\ 0 & S^k \end{bmatrix}. \qquad (33) $$

The matrix $S^k$ is defined as the generalized Schur complement matrix

$$ S^k := -(B_2^k (K^k)^{-1} B_1^\top + C^k). $$

This approach was proposed in the work [33, Section 10.1.2]. We choose an approximation $\hat{S}^k$ to $S^k$ implicitly by the action of $(\hat{S}^k)^{-1}$ on a given vector $v$. The Schur complement matrix is asymmetric and not given explicitly, so the approximation $\hat{S}^k$ to $S^k$ is defined by the solution with GMRES to a given tolerance. In an analogous manner we define an approximation $\hat{K}^k$ to $K^k$. The approximation $(\hat{K}^k)^{-1}$ has to be calculated once per GMRES step in the solving process with $\hat{S}^k$. To keep the effort low we realize this by a conjugated gradient solver preconditioned by an algebraic multigrid method. We perform one V-cycle with symmetric Gauss Seidel pre- and post-smoothing. For systems of moderate size we determine a LU-decomposition of $K^k$ instead and define $(\hat{K}^k)^{-1}v$ as a direct solution. Because the Matrix $K^k$ is constant during an active set step the multigrid structure resp. the LU-decomposition has to be calculated only once per active set step. The computations can be performed in parallel on the bodies because the system is not coupled. Due to the bad condition of the Schur complement matrix we precondition it with

$$ \mathcal{P}_S^k = \left[ \operatorname{diag} \left( -(B_2^k (\operatorname{diag}(K^k))^{-1} B_1^\top + C^k) \right) \right]^{-1}. \qquad (34) $$

The explicit form of $(\mathcal{P}^k)^{-1}$ reads

$$ (\mathcal{P}^k)^{-1} = \begin{bmatrix} (\hat{K}^k)^{-1} & 0 \\ 0 & I \end{bmatrix} \begin{bmatrix} I & B_1^\top \\ 0 & -I \end{bmatrix} \begin{bmatrix} I & 0 \\ 0 & -(\hat{S}^k)^{-1} \end{bmatrix}. $$

Basically the effort of applying $\mathcal{P}_S^k$ equates to one solution with the Schur complement matrix and one solution with $K^k$. The solution process of (32) is terminated if the residuals concerning contact, friction and plasticity $C_N^k$, $C_T^k$ and $\text{res}_{\text{plast}}$ as defined in Algorithm 1 underrun a given tolerance tol. Due to the plastic material we damp the displacements by a line search strategy [37] beginning in the step $k = 3$ of the algorithm.

# 5. Numerical examples

In this section we motivate higher-order methods and apply the presented algorithms to concrete selected numerical examples in two and three space dimensions. We can validate the successful operation of the method and the effectivity of the proposed preconditioning methods.

## 5.1. Motivation

We consider two flat rectangular elastic bodies $\Omega^1 = (0, 20) \times (0, 2)$ and $\Omega^2 = (18, 38) \times (-2.2, -0.2)$ that are clamped at the sides $\Gamma_0^1 = \{0\} \times [0, 2]$ and $\Gamma_0^2 = \{38\} \times [-2.2, -0.2]$. Volume forces $f^1 = (0, -0.0054)$ resp. $f^2 = (0, 0.0054)$ are acting and the material is chosen to be nearly incompressible by Poisson's ratios $\nu^1 = \nu^2 = 0.4999$. The Young's moduli are set to $E^1 = E^2 = 100$. In Fig. 1 the magnitude of the resulting deformations is depicted. We opposed bilinear and biquadratic ansatz functions. The choice of a bilinear ansatz results in a too low deformation and therefore bad approximation which represents locking effects, cf. Fig. 1(a). The contact situation is not resolved, because actually no contact occurs. The effect of locking occurs in the FE simulation of nearly incompressible materials or very thin media by the usage of conform elements with low-order polynomials [38, Chapter 6]. It results in much lower deformations than the true ones. Higher-order ansatz functions avoid these effects as shown in Fig. 1(b). The difference in the norm of deformations is given by a factor of nearly 2.3 and therefore the contact situation can be reproduced accurately. Another strategy to overcome these difficulties is Selective Reduced Integration (SRI), cf. [39, Section 5.3]. It is based on a split of the material tensor into its parts of compression and thrust. The problematic part of compression is integrated with a low-order 1-point Gauss quadrature rule. This method yields very similar results to the usage of polynomial degree 2, cf. Fig. 1(c).
<figure><img src="image_2.png" /><figcaption>Fig. 1. Comparison of deformations for different polynomial degrees and SRI.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Setting of the 2D two-body contact problem.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Deformations in x- and y-direction.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Plastic variables in $\Omega_1$.</figcaption></figure>

## 5.2. Two dimensional example

We consider the contact of an elasto-plastic and a linear-elastic body, which are represented by $\Omega^1 = (0, 1) \times (0, 1)$ and $\Omega^2 = (0.91, 1.91) \times (0, 1)$, cf. Fig. 2(a). The body $\Omega^1$ is subjected to Neumann boundary conditions given by $p^1 = (0, 30)$ on its lower side and by $p^1 = (30, 30)$ on its upper side. Homogeneous Dirichlet boundary conditions hold on $\Gamma_D^1 = \{0\} \times [0, 1]$ and $\Gamma_D^2 = \{1.91\} \times [0, 1]$. The two bodies may come into contact on the contact boundaries $\Gamma_C^1 = \{1\} \times [0, 1]$ respectively $\Gamma_C^2 = \{0.91\} \times [0, 1]$, which is caused by the overlapping of $\Omega_1$ and $\Omega_2$. Resulting displacements calculated on a globally refined mesh of size $n = 197.632$ and $m_1 = m_2 = 128$ with polynomial degrees $p=2, q=1$ are depicted in Fig. 3. For stability reasons we choose a meshsize of $\Gamma_C^1$ that fulfills $H = 2h$. The Young's moduli $E^1 = 10^4$, $E^2 = 10^3$ and Poisson's ratios $\nu^1 = 0.25$, $\nu^2 = 0.22$ define the compliance tensors $A^1$ and $A^2$. We choose hardening parameters $\gamma_{iso}^1 = 0.01$, $\gamma_{iso}^2 = 1.0$ and yield stresses $\sigma_0^1 = 10^2$, $\sigma_0^2 = 10^{10}$. These material parameters cause an ill-conditioned system and are selected to demonstrate the behavior of the algorithm and the effect of the preconditioners. The norm of the deviatoric part of the stresses $|\sigma(u_h^1)^D|_F$ and the portion of plastified quadrature points are shown in Fig. 4. We consider a friction law after Betten [40, Section 2.4] given by the
<figure><img src="image_6.png" /><figcaption>(a) $p = 1, q = 0$ and $p = 2, q = 1$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) $p = 1, q = 0$ and $p = 2, q = 1$.</figcaption></figure>

<caption>(d) $p = 1, q = 0$ and $p = 5, q = 4$.</caption>

<figure><img src="image_8.png" /><figcaption>(a) $p = 3, q = 2, n = 236.480, m_1 = m_2 = 384$.</figcaption></figure>

<caption>(b) $p = 4, q = 3, n = 378.112, m_1 = m_2 = 512$.</caption>

<caption>Fig. 6. Residuals for different polynomial degrees.</caption>

functional

$$s(\lambda_n) := Y_0 \left( \tanh \left( \left[ \frac{\mu |\lambda_n|}{Y_0} \right]^p \right) \right)^{\frac{1}{p}}$$

with parameters $p = 2, \mu = 0.1$ and $Y_0 = 10$. Furthermore the constants in the NCP functions are chosen by $c_n = 1$ and $c_t = E^1$, whereas the tolerance of Algorithm 1 is set to tol = $10^{-10}$. Calculated Lagrange multipliers $\lambda_n$ and $\lambda_t$ for selected polynomial degrees on a locally refined mesh (cf. Fig. 2(b)) are depicted in Fig. 5. The nonoscillating course of these functions indicates the stability of the underlying mixed schemes.

The development of the plastic, frictional and geometrical contact residuals for different polynomial degrees as defined in Algorithm 1 is presented in Fig. 6. For both choices of discretizations $p = 3, q = 2$ and $p = 4, q = 3$ the frictional residual dominates the others permanently. Due to the simple geometric contact situation of exclusively
<table><caption>Table 1
Comparison of preconditioners with $p = 1$, $q = 0$ on a uniform refined mesh, $n = 288$, $m_1 = m_2 = 4$, tol = $10^{-9}$.</caption><thead><tr><th></th><th>Step</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th></tr></thead><tbody><tr><td>No prec.</td><td>#iter</td><td>589</td><td>5401</td><td>10,569</td><td>16,838</td><td>15,463</td><td>5694</td><td>1595</td><td>316</td><td>267</td></tr><tr><td>P<sup>k</sup></td><td>#iter</td><td>5</td><td>4</td><td>4</td><td>4</td><td>4</td><td>3</td><td>3</td><td>2</td><td>/</td></tr><tr><td>P<sup>k</sup>, P<sub>S</sub><sup>k</sup></td><td>#iter</td><td>19</td><td>382</td><td>1,036</td><td>54</td><td>36</td><td>27</td><td>105</td><td>1.016</td><td>/</td></tr><tr><td></td><td>#iter<sub>S</sub></td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>2</td><td>/</td></tr><tr><td></td><td>#iter<sub>S</sub></td><td>13</td><td>27</td><td>27</td><td>27</td><td>27</td><td>27</td><td>27</td><td>18</td><td>/</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Convergence of the error in the energy norm with uniform mesh refinement.</figcaption></figure>

active indices the quantity $C_N^k$ tends to zero very fast. After some iterations in which the displacements have to be damped ($\alpha = 0.5$) the residuals $C_T^k$ and $\text{res}_{\text{plast}}$ also decrease rapidly. Overall the general behavior appears to be similar for both discretizations.

Table 1 shows a comparison of different choices for preconditioning methods applied at the solution process of a problem of very small size. The number of iterations of a Restarted GMRES solver without any preconditioning is opposed to those of a supplemental application of a block triangular preconditioner $P^k$ defined in (33). In the third version the generalized Schur complement matrix $\hat{S}^k$ within the evaluation of $P^k$ is additionally preconditioned by $P_S^k$ given in (34). The number of outer GMRES steps is named as iter, whereas iter<sub>S</sub> labels the accumulated number of GMRES iterations that are performed to solve with the Schur complement. Every evaluation of $\hat{S}^k$ includes a solution with the block-diagonal matrix $\hat{K}^k$. Due to the chosen friction law and the start value $\lambda_{n,H}^0 = 0$ the friction bound and therefore $\lambda_{t,H}^1$ vanishes. This results in a well-conditioned system matrix and a moderate amount of solving iterations in the first semi-smooth Newton step. In the following steps the number of outer iterations rises rapidly if no preconditioning is performed. This indicates the bad condition aroused by the Newton system concerning frictional constraints. In the later Newton steps due to the revising starting value the amount of solving iterations decreases again. The preconditioned versions reduce the number of Newton steps by one. Furthermore preconditioning the Schur complement matrix decreases iter<sub>S</sub> and therefore the costly amount of direct solutions with $K^k$ considerably. The condition of the system matrix is very bad in this application basically due to the frictional constraints, which necessitates preconditioning. An acceptable effort is achieved by the suggested combination of preconditioners.

Eventually we investigate the convergence order in the 2D problem on a uniformly refined mesh. Fig. 7 shows the $H^1$-error to a reference solution calculated on a fine mesh of meshsize $h = 0.5^8$ ($p = 1$) respectively $h = 0.5^7$ ($p = 2$) for bilinear and biquadratic ansatz functions. In both cases we observe $O(h^{1/2})$-convergence which is not optimal. This is caused by the low regularity of the solution due to the elasto-plastic material behavior and the contact situation including friction. Typically the convergence orders are limited by $O(h)$ respectively $O(h^{3/2})$ for frictionless contact problems in linear elasticity [15,41,42] and linear-elastic problems with volumetric inequality
<figure><img src="image_10.png" /><figcaption>Fig. 8. Overview of the 3D example.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Plastic variables of the mounted point $\Omega_1$.</figcaption></figure>

constraints [43]. The full benefit of optimal convergence for higher-order methods can only be achieved by using adaptivity [31,36].

### 5.3. Threedimensional example

As an application from mechanical engineering we consider a two-body contact problem that models a grinding process. The first body $\Omega_1$ represents a cylindrical mounted point with diameter 2.2. Its spindle is clamped at the end. A hexahedral body $\Omega_2 = (-2, 3) \times (1.1, 3.1) \times (-2.5, 2.5)$ represents a workpiece, which is machined by the mounted point. An overview of the contact situation, the corresponding FE meshes, and the resulting deformations of both bodies are depicted in Fig. 8. Around the contact zones and the clamped part of the shaft the meshes are locally refined. Dirichlet boundary conditions hold on its side surfaces and the bottom. We consider a Coulomb–Orowan friction law [3, Section 4.2.5]

$$s(\lambda_n) := \min\{\mu|\lambda_n|, Y_0\}$$

with constants $\mu = 0.5$ and $Y_0 = 1.1$. In contrast to the two dimensional example the present meshes are nonmatching. The elasto-plastic material parameters are chosen by $E^1 = 10^2$, $E^2 = 10^3$, $\nu^1 = \nu^2 = 0.3$, $\sigma_0^1 = 3$, $\sigma_0^2 = 2$ and $\gamma_{iso}^1 = \gamma_{iso}^2 = 0.1$. Resulting deviators and portions of plastified quadrature points are pictured in Figs. 9 and 10. As might be expected the workpiece plastifies at the contact zone. In contrast at the mounted point, the deviatoric part of the stresses gets maximal around the clamped part of the spindle. Resulting Lagrange multipliers for frictional and geometrical contact constraints are shown in Fig. 11. The discretization appears to be stable because no checkerboard patterns are observed.

A comparison of different choices for the constant $c_t$ inside the NCP function is shown in Fig. 12. The development of the tangential and plastic residual is pictured. As values of $c_t$ the Elasticity moduli $E^1 = 10^2$ and $E^2 = 10^3$ are chosen. It is conspicuous that the number of Newton iterations is significantly larger for $c_t = 10^3$. Furthermore for this constant the development of the residual $C_T^k$ appears more unstable. This behavior indicates the sensitivity of the method concerning the choice of the constant $c_t$.
<figure><img src="image_12.png" /><figcaption>Fig. 10. Resulting mechanical variables in the contact zone of workpiece $\Omega_2$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Lagrange multipliers $\lambda_n$ and $\lambda_t$ for $p = 1, q = 0, n_2 = 65.895, n_2 = 20.358, m_1 = 256, m_2 = 512$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Development of $C_T^k$ and $C_N^k$ for $c_t = 10^2$ (a) and $c_t = 10^3$ (b) ($p = 1, q = 0, n_2 = 65.895, n_2 = 20.358, m_1 = 256, m_2 = 512$).</figcaption></figure>

## 6. Conclusions and outlook

In this paper an active-set solution strategy for mixed higher-order discretizations of frictional two-body contact problems is presented. For this ansatz a construction of higher-order dual basis functions is not necessarily compared to Mortar discretization. We suggest a block triangular preconditioner for the full saddle point system and an appropriate preconditioner for the generalized Schur complement matrix. Numerical results show the successful operation of the semi-smooth Newton method and a considerable reduction of the condition by the preconditioners.

In further work following [19] we will extend this solving method by reformulating plastic material behavior in terms of two additional NCP functions. Furthermore an extension of the contact model by including effects like adhesion [44] and separation of two bodies [3] as well as the development of appropriate adaptive strategies are planed.

## Acknowledgments

The authors gratefully acknowledge the support by the German Research Foundation (DFG) through subproject A5 "Simulation supported NC-shape grinding as a finishing operation of thermally coated deep drawing tools"
within the Collaborative Research Center (SFB) 708 "3D-Surface Engineering of Tools for Sheet Metal Forming - Manufacturing, Modeling, Machining".

## References

[1] N. Kikuchi, J.T. Oden, Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods, SIAM, Philadelphia, 1988. http://dx.doi.org/10.1137/1.9781611970845.

[2] T. Siebrecht, D. Biermann, H. Ludwig, S. Rausch, P. Kersting, H. Blum, A. Rademacher, Simulation of grinding processes using finite element analysis and geometric simulation of individual grains, Prod. Eng. Res. Dev. 8 (2014) 345-353. http://dx.doi.org/10.1007/s11740-013-0524-9.

[3] P. Wriggers, Computational Contact Mechanics, second ed., Springer, Berlin and New York, 2006. http://dx.doi.org/10.1007/978-3-540-32609-0.

[4] B. Wohlmuth, Variationally consistent discretization schemes and numerical algorithms for contact problems, Acta Numer. 20 (2011) 569-734. http://dx.doi.org/10.1017/S0962492911000079.

[5] M. Hintermüller, K. Ito, K. Kunisch, The primal-dual active set strategy as a semismooth newton method, SIAM J. Optim. 13 (3) (2003) 865-888. http://dx.doi.org/10.1137/S1052623401383558.

[6] J.-S. Pang, Newton's method for b-differentiable equations, Math. Oper. Res. 15 (2) (1990) 311-341. http://dx.doi.org/10.1287/moor.15.2.311.

[7] L. Qi, J. Sun, A nonsmooth version of Newton's method, Math. Program. 58 (1) (1993) 353-367. http://dx.doi.org/10.1007/BF01581275.

[8] P. Alart, A. Curnier, A mixed formulation for frictional contact problems prone to Newton like solution methods, Comput. Methods Appl. Mech. Engrg. 92 (3) (1991) 353-375. http://dx.doi.org/10.1016/0045-7825(91)90022-X.

[9] S. Hüeber, Discretization techniques and efficient algorithms for contact problems (Ph.D. thesis), Universität Stuttgart, 2008. URL http://elib.uni-stuttgart.de/opus/volltexte/2008/3608.

[10] S. Hüeber, G. Stadler, B. Wohlmuth, A primal-dual active set algorithm for three-dimensional contact problems with Coulomb friction, SIAM J. Sci. Comput. 30 (2) (2008) 572-596. http://dx.doi.org/10.1137/060671061.

[11] S. Sitzmann, K. Willner, B. Wohlmuth, A dual Lagrange method for contact problems with regularized frictional contact conditions: Modelling micro slip, Comput. Methods Appl. Mech. Engrg. 285 (2015) 468-487. http://dx.doi.org/10.1002/nme.4683.

[12] T. Koziara, N. Biéanić, Semismooth Newton method for frictional contact between pseudo-rigid bodies, Comput. Methods Appl. Mech. Engrg. 197 (33-40) (2008) 2763-2777. http://dx.doi.org/10.1016/j.cma.2008.01.006.

[13] A. Popp, M. Gitterle, M.W. Gee, W.A. Wall, A dual mortar approach for 3D finite deformation contact with consistent linearization, Internat. J. Numer. Methods Engrg. 83 (11) (2010) 1428-1465. http://dx.doi.org/10.1002/nme.2866.

[14] M. Gitterle, A. Popp, M. Gee, W. Wall, Finite deformation frictional mortar contact using a semi-smooth Newton method with consistent linearization, Internat. J. Numer. Methods Engrg. 84 (5) (2010) 543-571. http://dx.doi.org/10.1002/nme.2907.

[15] S. Hüeber, M. Mair, B. Wohlmuth, A priori error estimates and an inexact primal-dual active set strategy for linear and quadratic finite elements applied to multibody contact problems, Appl. Numer. Math. 54 (2005) 555-576. http://dx.doi.org/10.1016/j.apnum.2004.09.019.

[16] M. Hintermüller, S. Rösel, A duality-based path-following semismooth Newton method for elasto-plastic contact problems, J. Comput. Appl. Math. 292 (2016) 150-173. http://dx.doi.org/10.1016/j.cam.2015.06.010.

[17] J.C. De los Reyes, M. Hintermüller, A duality based semismooth Newton framework for solving variational inequalities of the second kind, Interfaces Free Bound. 13 (4) (2011) 437-462.

[18] P.W. Christensen, A semi-smooth Newton method for elasto-plastic contact problems, Int. J. Solids Struct. 39 (8) (2002) 2323-2341. http://dx.doi.org/10.1016/S0020-7683(02)00086-0.

[19] C. Hager, B. Wohlmuth, Nonlinear complementarity functions for plasticity problems with frictional contact, Comput. Methods Appl. Mech. Engrg. 198 (41-44) (2009) 3411-3427. http://dx.doi.org/10.1016/j.cma.2009.06.021.

[20] A. Seitz, A. Popp, W.A. Wall, A semi-smooth newton method for orthotropic plasticity and frictional contact at finite strains, Comput. Methods Appl. Mech. Engrg. 285 (2015) 228-254. http://dx.doi.org/10.1016/j.cma.2014.11.003.

[21] Z. Dostál, J. Haslinger, R. Kučera, Implementation of the fixed point method in contact problems with coulomb friction based on a dual splitting type technique, J. Comput. Appl. Math. 140 (1-2) (2002) 245-256. http://dx.doi.org/10.1016/S0377-0427(01)00405-8.

[22] J. Simo, T. Laursen, An augmented lagrangian treatment of contact problems involving friction, Comput. Struct. 42 (1) (1992) 97-116. http://dx.doi.org/10.1016/0045-7949(92)90540-G.

[23] Z. Dostál, D. Horák, R. Kučera, V. Vondrák, J. Haslinger, J. Dobiáš, S. Pták, FETI based algorithms for contact problems: scalability, large displacements and 3d coulomb friction, Comput. Methods Appl. Mech. Engrg. 194 (2-5) (2005) 395-409. http://dx.doi.org/10.1016/j.cma.2004.05.015.

[24] R. Krause, B. Wohlmuth, A Dirichlet-Neumann type algorithm for contact problems with friction, Comput. Vis. Sci. 5 (2002) 139-148. http://dx.doi.org/10.1007/s00791-002-0096-2.

[25] R. Krause, Monotone multigrid methods for Signorini's problem with friction (Ph.D. thesis), Freie Universität Berlin, 2001. URL http://www.diss.fu-berlin.de/diss/receive/FUDISS_thesis_000000000469.

[26] H. Blum, D. Braess, F. Suttmeier, A cascadic multigrid algorithm for variational inequalities, Comput. Vis. Sci. 7 (3-4) (2004) 153-157. http://dx.doi.org/10.1007/s00791-004-0134-3.

[27] L. Banz, hp-finite element and boundary element methods for elliptic, elliptic stochastic, parabolic and hyperbolic obstacle and contact problems (Ph.D. thesis), Gottfried Wilhelm Leibniz Universität Hannover, 2011. URL http://edok01.tib.uni-hannover.de/edoks/e01dh12/686068297.pdf.

[28] J. Haslinger, Mixed formulation of elliptic variational inequalities and its approximation, Appl. Math. 26 (1981) 462-475. URL http://eudml.org/doc/15217.
[29] A. Schröder, H. Blum, H. Kleemann, Mixed finite element methods for two-body contact problems, J. Comput. Appl. Math. 283 (2015) 58–70. http://dx.doi.org/10.1016/j.cam.2014.11.040.

[30] A. Schröder, H. Blum, A. Rademacher, H. Kleemann, Mixed FEM of higher order for contact problems with friction, Int. J. Numer. Anal. Model. 8 (2) (2011) 302–323.

[31] H. Kleemann, Adaptive FEM für Mehrkörperkontaktprobleme (Ph.D. thesis), Technische Universität Dortmund, 2011. URL http://hdl.handle.net/2003/29117.

[32] R. Rannacher, F. Suttmeier, A posteriori error estimation and mesh adaptation for finite element models in elasto-plasticity, Comput. Methods Appl. Mech. Engrg. 176 (1999) 333–361. http://dx.doi.org/10.1016/S0045-7825(98)00344-2.

[33] M. Benzi, G. Golub, J. Liesen, Numerical solution of saddle point problems, Acta Numer. 14 (2005) 1–137. http://dx.doi.org/10.1017/S0962492904000212.

[34] S. Hüeber, B. Wohlmuth, A primal-dual active set strategy for non-linear multibody contact problems, Comput. Methods Appl. Mech. Engrg. 194 (27–29) (2005) 3147–3166. http://dx.doi.org/10.1016/j.cma.2004.08.006.

[35] T. Laursen, Computational Contact and Impact Mechanics: Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis, Springer, Heidelberg, 2002.

[36] P. Dörsek, J. Melenk, Adaptive hp-FEM for the contact problem with Tresca friction in linear elasticity, Appl. Numer. Math. 60 (7) (2010) 689–704. http://dx.doi.org/10.1016/j.apnum.2010.03.011.

[37] J. Nocedal, S. Wright, Numerical Optimization, second ed., Springer, New York, 2006. http://dx.doi.org/10.1007/978-0-387-40065-5.

[38] D. Braess, Finite elements in solid mechanics, in: Finite Elements, third ed., Cambridge University Press, 2007, Cambridge Books Online. URL http://dx.doi.org/10.1017/CBO9780511618635.009.

[39] T. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover Publications, 2000.

[40] D. Horwatitsch, Entwicklung eines Reibmodells für hohe Temperaturen und hohe Umformgrade (Ph.D. thesis), Universität Kassel, 2012. URL http://www.upress.uni-kassel.de/katalog/abstract.php?978-3-86219-646-3.

[41] B. Wohlmuth, A. Popp, M. Gee, W. Wall, An abstract framework for a priori estimates for contact problems in 3D with quadratic finite elements, Comput. Mech. 49 (6) (2012) 735–747. http://dx.doi.org/10.1007/s00466-012-0704-z.

[42] P. Coorevits, P. Hild, K. Lhalouani, T. Sassi, Mixed finite element methods for unilateral problems: Convergence analysis and numerical studies, Math. Comp. 71 (237) (2002) 1–25. http://dx.doi.org/10.1090/S0025-5718-01-01318-7.

[43] N. Bonfils, N. Chevaugon, N. Moës, Treating volumetric inequality constraint in a continuum media with a coupled X-FEM/level-set strategy, Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 16–28. http://dx.doi.org/10.1016/j.cma.2011.02.012.

[44] T. Dickopf, R. Krause, Efficient simulation of multi-body contact problems on complex geometries: A flexible decomposition approach using constrained minimization, Internat. J. Numer. Methods Engrg. 77 (13) (2009) 1834–1862. http://dx.doi.org/10.1002/nme.2481.

## Glossary

### Spaces

$V$: Continuous ansatz space $(H_D(\Omega^1))^d \times (H_D(\Omega^2))^d$

$V_{h,m}^p$: Set of ansatz functions of polynomial degree $p$ on body $m$ with meshsize $h$

$V_h$: Tensor space $V_h = V_{h,1}^{p_1} \times V_{h,2}^{p_2}$

$\Lambda_n, \Lambda_t$: Continuous normal and tangential Lagrange ansatz spaces

$S_f$: Polynomial tensor product space of order $l$ on $[-1, 1]^l$

$\Lambda_{n,H}, \Lambda_{t,H}$: Admissible set for discrete Lagrange multipliers concerning contact and friction

### Functions

$\sigma$: Stress

$\sigma_n, \sigma_{nn}, \sigma_{nt}$: Surface stresses, normal part of $\sigma_n$ and tangential stress vectors

$\gamma_M^m$: Trace operator $\gamma_M: H^1(\Omega^m) \to L^2(M)$ for $M \subset \Gamma^m$

$\varepsilon$: Linearized strain

$A^m$: Fourth order compliance tensor on body $m$

$\varepsilon^{m,p}$: Plastic part of strain on body $m$

$\tau^D$: Deviatoric part of tensor $\tau$

$\mathcal{F}^{m,\text{iso}}$: Yield function

$\sigma_0^m$: Yield stress on body $m$

$p^m$: Surface stresses

$[\cdot]_{n_\delta}, [\cdot]_{t_\delta}$: Normal and tangential jumps

$g$: Gap function

$\theta$: Transfer operator $\theta: L^2(\Gamma_C^2) \to L^2(\Gamma_C^1)$

$P_{\Pi}$: Projector that maps stresses on the admissible set

$\mu^m$: Shear modulus of $m$th body material

$j$: Functional describing frictional constraints

$a(\cdot)(\cdot)$: Semi-linear form $a: V \times V \to \mathbb{R}$

$f_{\text{ext}}$: External energy

$s$: Friction law
L: Lagrange functional

$F_m, T, F_E$: Transformation on $T \in \mathcal{T}_{h,m}$ resp. $E \in \mathcal{B}_H$

$C_N, C_T$: NCP functions concerning normal contact and friction

$\chi$: Characteristical function

$\mathcal{I}, \mathcal{J}$: Index maps

$\mathcal{U}, \mathcal{V}, \mathcal{W}, \mathcal{R}$: Functions defining the Newton equation

Sets

$\Omega^m$: Deformable body $m = 1, 2$

$\Gamma^m, \Gamma_D^m, \Gamma_N^m, \Gamma_C^m$: Boundary of $\Omega^m$, Dirichlet, Neumann and contact boundary

$\mathcal{T}_{h,m}, \mathcal{B}_{h,m}$: Triangulations of $\Omega^m$ and $\Gamma_C^m$

$C_q$: Set of $q + 1$ Gauss-quadrature points

$\mathcal{A}_n^k, \mathcal{T}_n^k$: Sets of active/inactive indices concerning geometrical contact

$\mathcal{A}_t^k, \mathcal{T}_t^k, \mathcal{T}_{tn}^k$: Sets of active/inactive indices concerning friction resp. indices without contact

Vectors and matrices

$K$: Stiffness matrix

$N$: Normal coupling matrix

$M$: Mass matrix

$n, n_\delta$: Outer normal vector on $\Gamma^m$ and generalized normal vector on $\Gamma_C^{-1}$

$\tilde{g}$: Gap vector

$T, \tilde{T}$: Unweighted and weighted tangential coupling matrix

$S$: Schur complement matrix

$\tilde{u}_h, \tilde{\lambda}_{n,H}, \tilde{\lambda}_{t,H}$: Vector-valued representations of $u_h, \lambda_{n,H}$ and $\lambda_{t,H}$

$s$: Friction bound

$U, V, W, \tilde{r}$: Vector-valued representations of $\mathcal{U}, \mathcal{V}, \mathcal{W}, \mathcal{R}$

$\mathcal{P}^k$: Block triangular preconditioner

$\mathcal{P}_S^k$: Preconditioner of Schur complement matrix

$\hat{S}, \hat{K}$: Approximations to $S$ resp. $K$ by solving with GMRES

Constants

$n, n_m$: Dimensions of $V_h$ resp. $V_{h,m}^p$

$m_1, m_2$: Dimensions of $\Lambda_{n,H}$ resp. $\Lambda_{t,H}$

$h, H$: Mesh sizes of $\Omega^m$ and $\Gamma_C^m$

$p$: Polynomial degree of primal ansatz functions

$q$: Polynomial degree of Lagrange ansatz functions

$k$: Iteration index of Newton method

$c_n, c_t$: Constants in NCP functions $C_N$ and $C_T$