Heribert Blum

Hannah Frohne

Jörg Frohne

Andreas Rademacher
email: andreas.rademacher@tu-dortmund.de

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

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

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.

Introduction

Frictional, elasto-plastic multi-body contact problems play an important role in mechanical engineering [START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF][START_REF] Siebrecht | Simulation of grinding processes using finite element analysis and geometric simulation of individual grains[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF]. 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 [START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF]. 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 [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth newton method[END_REF] 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 [START_REF] Pang | Newton's method for b-differentiable equations[END_REF][START_REF] Qi | A nonsmooth version of Newton's method[END_REF][START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF]. 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 [START_REF] Hüeber | Discretization techniques and efficient algorithms for contact problems[END_REF][START_REF] Hüeber | A primal-dual active set algorithm for three-dimensional contact problems with Coulomb friction[END_REF][START_REF] Sitzmann | A dual Lagrange method for contact problems with regularized frictional contact conditions: Modelling micro slip[END_REF][START_REF] Koziara | Semismooth Newton method for frictional contact between pseudo-rigid bodies[END_REF]. The papers [START_REF] Popp | A dual mortar approach for 3D finite deformation contact with consistent linearization[END_REF][START_REF] Gitterle | Finite deformation frictional mortar contact using a semi-smooth Newton method with consistent linearization[END_REF] show a dual mortar approach for finite deformation contact problems. Linear and quadratic finite elements have been regarded in [START_REF] Hüeber | A priori error estimates and an inexact primal-dual active set strategy for linear and quadratic finite elements applied to multibody contact problems[END_REF]. Semi-smooth Newton methods for elasto-plastic problems are analyzed involving frictionless contact in [START_REF] Hintermüller | A duality-based path-following semismooth Newton method for elasto-plastic contact problems[END_REF][START_REF] De Los Reyes | A duality based semismooth Newton framework for solving variational inequalities of the second kind[END_REF] and frictional contact in [START_REF] Christensen | A semi-smooth Newton method for elasto-plastic contact problems[END_REF][START_REF] Hager | Nonlinear complementarity functions for plasticity problems with frictional contact[END_REF][START_REF] Seitz | A semi-smooth newton method for orthotropic plasticity and frictional contact at finite strains[END_REF]. Besides Mortar methods there exist a series of strategies for the numerical treatment of frictional contact problems.

One approach forms classical fixpoint methods [START_REF] Dostál | Implementation of the fixed point method in contact problems with coulomb friction based on a dual splitting type technique[END_REF]. General ideas with penalty methods and Lagrange multiplier techniques are described in [START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF]. A combination of both is augmented Lagrange multiplier techniques [START_REF] Simo | An augmented lagrangian treatment of contact problems involving friction[END_REF]. Domain decomposition methods as FETI techniques [START_REF] Dostál | FETI based algorithms for contact problems: scalability, large displacements and 3d coulomb friction[END_REF] or Dirichlet-Neumann algorithms in combination with Mortar discretizations [START_REF] Krause | A Dirichlet-Neumann type algorithm for contact problems with friction[END_REF] are efficient parallel computable solvers. By monotone multigrid constructed global convergent solvers are suggested in [START_REF] Krause | Monotone multigrid methods for Signorini's problem with friction[END_REF]. A cascadic multigrid algorithm for variational inequalities is presented in [START_REF] Blum | A cascadic multigrid algorithm for variational inequalities[END_REF]. Under the use of higher-order DG-discretizations active set strategies have been applied to linear-elastic obstacle problems [START_REF] Banz | hp-finite element and boundary element methods for elliptic, elliptic stochastic[END_REF].

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 [START_REF] Hüeber | Discretization techniques and efficient algorithms for contact problems[END_REF][START_REF] Hüeber | A primal-dual active set algorithm for three-dimensional contact problems with Coulomb friction[END_REF] to mixed finite elements introduced by Haslinger [START_REF] Haslinger | Mixed formulation of elliptic variational inequalities and its approximation[END_REF] and higher-order discretizations presented in [START_REF] Schröder | Mixed finite element methods for two-body contact problems[END_REF][START_REF] Schröder | Mixed FEM of higher order for contact problems with friction[END_REF]. 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 [START_REF] Schröder | Mixed FEM of higher order for contact problems with friction[END_REF][START_REF] Kleemann | Adaptive FEM für Mehrkörperkontaktprobleme[END_REF] a solution scheme for linear-elastic, frictional multi-body contact problems and higherorder 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 [START_REF] Rannacher | A posteriori error estimation and mesh adaptation for finite element models in elasto-plasticity[END_REF]. 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 [START_REF] Benzi | Numerical solution of saddle point problems[END_REF] 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 [START_REF] Kleemann | Adaptive FEM für Mehrkörperkontaktprobleme[END_REF][START_REF] Hüeber | A primal-dual active set strategy for non-linear multibody contact problems[END_REF].

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.

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 Ω m with m = 1, 2 in d = 2 or 3 spatial dimensions, on which volume forces f m ∈ L 2 (Ω m , R d ) are acting. Their boundaries are denoted by Γ m , m = 1, 2. With the outer normal vector n on Γ m one can define the surface stresses σ n (u m ) := σ (u m )n. The normal part of these stresses is given by the scalar value σ nn := n ⊤ σ (u m )n, whereas the tangential stress vectors are calculated by σ nt,i := n ⊤ σ (u m )t i . The matrix t ∈ R d×d-1 contains the tangential vectors on Γ m that build an orthonormal system with the outer normal n. Define the trace operators γ m M :

H 1 (Ω m ) → L 2 (M) for M ⊂ Γ m and H D (Ω m ) :=  v ∈ H 1 (Ω m )   γ Γ m D (v) = 0  . Their d-dimensional cartesian product space is denoted by V :=  H D (Ω 1 )  d ×  H D (Ω 2 )  d .
We are interested in the displacements u = (u 1 , u 2 ) ∈ V , which fulfill the following conditions for m = 1, 2:

ε(u m ) = A m σ (u m ) + ε m,P in Ω m (1) -divσ (u m ) = f m in Ω m (2) ε m,P (τ -σ (u m )) ≥ 0 ∀τ with F m,iso (τ, |ε m,P | F ) ≤ 0 in Ω m (3) u m = 0 on Γ m D (4) σ n (u m ) = p m on Γ m N . (5)
Relation (1) describes the material law, the relation between the linearized strain ε(u m ) = 1 2  ∇u m + ∇u m,⊤  and the stress σ (u m ). The strain is split up into an elastic part A m σ (u m ) with the fourth order compliance tensor A m corresponding to isotropic material and a plastic part ε m,P . Eq. ( 2

:= Γ m \  Γ m D ∪ Γ m C  surface stresses p m ∈ L 2 (Γ m N , R d ) act on each body Ω m by relation (5). Let Φ : Γ 1 C → Γ 2
C be an appropriate, bijective, sufficiently smooth mapping between the contact boundaries of the slave body Ω 1 and the master body Ω 2 . In practical applications the mapping function Φ 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 [START_REF] Laursen | Computational Contact and Impact Mechanics: Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis[END_REF]Chapter 4.2]. In order to model contact conditions we define a generalized normal vector for x ∈ Γ 1 by

n δ (x) :=    Φ(x) -x ∥Φ(x) -x∥ , x ̸ = Φ(x) n 1 (x) = -n 2 (x), x = Φ(x)
and corresponding tangential vectors t δ ∈ L 2 (Γ 1 ) d×(d-1) such that (n δ (x), t δ (x)) form an orthonormal system. We define the normal jump

[v] n δ (x) := γ Γ 1 C (v 1 )(x) • n δ (x) -γ Γ 2 C (v 2 )(Φ(x)) • n δ (x)
and the tangential jump

[v] t δ (x) := t δ (x) ⊤ γ Γ 1 C (v 1 )(x) -t δ (x) ⊤ γ Γ 2 C (v 2 )(Φ(x))
on Γ 1 C . The distance of Ω 1 and Ω 2 is given by the gap function g(x) := |Φ(x) -x|. With this notion the contact conditions, which are defined on the slave body, read as follows:

[u] n δ ≤ g on Γ 1 C (6) σ n δ n δ (u 1 ) ≤ 0 on Γ 1 C (7) σ n δ n δ (u 1 )([u] n δ -g) = 0 on Γ 1 C ( 8
)
σ n δ (u 1 ) = -Θ * σ n δ (u 2 ) on Γ 1 C . (9)
The bodies are not allowed to penetrate each other [START_REF] Pang | Newton's method for b-differentiable equations[END_REF] and only negative or vanishing contact forces are allowed [START_REF] Qi | A nonsmooth version of Newton's method[END_REF]. By the complementarity condition (8) the property that either contact occurs or the normal contact forces vanish is modeled. Eq. ( 9) contains the adjoint Θ * of a transfer operator Θ :

L 2 (Γ 2 C ) → L 2 (Γ 1 C ) which is defined by Θ(v 2 )(x 1 ) := v(Φ(x 1 )
) and ensures equality of the contact forces on Γ 1 C and Γ 2 C . Besides the described normal contact we regard frictional constraints on

Γ 1 C :   σn δ t δ (u 1 )    ≤ s(σ n δ n δ (u 1 )) (10)
  σn δ t δ (u 1 )    < s(σ n δ n δ (u 1 )) ⇒ [u] t δ = 0 (11)   σn δ t δ (u 1 )    = s(σ n δ n δ (u 1 )) ⇒ ∃α ∈ R ≥0 : [u] t δ = ασ n δ t δ (u) (12)
with the euclidean norm ∥ • ∥. 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 [START_REF] Sitzmann | A dual Lagrange method for contact problems with regularized frictional contact conditions: Modelling micro slip[END_REF]. Otherwise the slip condition [START_REF] Koziara | Semismooth Newton method for frictional contact between pseudo-rigid bodies[END_REF] 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 [START_REF] Rannacher | A posteriori error estimation and mesh adaptation for finite element models in elasto-plasticity[END_REF] we introduce a primal-mixed formulation of the frictional elasto-plastic two-body contact problem (1)- [START_REF] Koziara | Semismooth Newton method for frictional contact between pseudo-rigid bodies[END_REF] and project the stresses onto the admissible set by the projector

P Π (τ ) :=    τ, |τ D | F ≤ σ m 0  γ m iso 2µ m + γ m iso +  1 - γ m iso 2µ m + γ m iso  σ m 0 |τ D | F  τ D + 1 d tr(τ )I d×d , |τ D | F > σ m 0
with the shear modulus µ m of the mth body material. This yields a variational inequality in the displacements

a(u)(ϕ -u) -f ext (ϕ -u) + j (ϕ) -j (u) ≥ 0 ∀ϕ ∈ K on the convex set K :=  v ∈ V | [v] n δ ≤ g  ⊂ V
containing the contact conditions and with the functional

j (ϕ) :=  Γ 1 C s∥[ϕ] t δ ∥dσ
describing the frictional constraints. The semi-linearform a is defined by

a : V × V → R, a(v)(w) := 2  m=1  Ω m P Π  (A m ) -1 ε(v m )  : ε(w m )d x,
and the linear form

f ext (v) := 2  m=1 ( f m , v m ) 0 + ( p m , v m ) 0,Γ m N
corresponds to the external energy with the L 2 -scalar product (•, •) 0 . Introducing the Lagrange functional,

L (v, µ n , µ t ) := 1 2 a(v)(v) -f ext (v) + ⟨µ n , [v] n δ -g⟩ + (µ t , s[v] t δ ) 0,Γ 1 C the dual cone Λ n :=  µ ∈ H 1 2 (Γ 1 C )    µ ≤ 0 a.e.  ′ and
Λ t :=  µ ∈ L 2 (Γ 1 C ) d-1    ∥µ∥ ≤ s 
we replace the contact and the frictional constraints and end up in a mixed formulation of the form:

Find (u, λ n , λ t ) ∈ V × Λ n × Λ t such that a(u)(v) + b n (λ n , v) + b t (λ t , v) = f ext (v) ∀v ∈ V (13) b n (µ n -λ n , u) + b t (µ t -λ t , u) ≥ ⟨µ n -λ n , g⟩ ∀µ n ∈ Λ n , ∀µ t ∈ Λ t . (14)
Here ⟨•, •⟩ denotes the dual pairing on Λ n . The bilinear forms that include the contact and frictional conditions are defined by

b n : Λ n × V → R, b n (µ, w) := ⟨µ, [w] n δ ⟩, b t : Λ t × V → R, b t (µ, w) := (µ, s[w] t δ ) 0,Γ C .

Discretization

In this section we introduce a higher-order discretization for the described problem ( 13)-( 14) that was mentioned in [START_REF] Schröder | Mixed finite element methods for two-body contact problems[END_REF]. Denote by T h,m and B H triangulations of Ω m , m = 1, 2 and Γ 1 C , respectively, with corresponding affine transformations

F m,T : T := [-1, 1] d → T ∈ T h,m
and

F E : Ê := [-1, 1] d-1 → E ∈ B H .
Let S r l be the polynomial tensor product space of order r on the reference element [-1, 1] l . We define ansatz functions of polynomial degree p ∈ N on body m = 1, 2 by

V p h,m :=  v ∈ (H D (Ω m )) d    ∀T ∈ T h,m : v| T • F m,T ∈ S p d (T )  =: span{φ m i } n m i=1
and the tensor space

V h = V p 1 h,1 × V p 2 h,2 .
The cumulated dimensions of the discrete spaces V p h,m are named n := n 1 + n 2 . By introduction of the space

M q H :=  v ∈ L 2 (Γ 1 C )    ∀E ∈ B H : v| E • F E ∈ S q d-1 (E)
 with a polynomial degree q ∈ N the admissible set for the Lagrange multipliers concerning geometrical contact is defined as follows:

Λ n,H :=  v ∈ M q H   ∀E ∈ B H : ∀x ∈ C q : v(F E (x)) ≤ 0  =: span  ψ n i  m 1 i=1 . The sign condition v(F E (x)) ≤ 0 is only defined on the finite set C q ⊂ [-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

 Γ 1 C λ n,H ψ H do ≥ 0 for functions ψ H ∈ Λ 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 [START_REF] Dörsek | Adaptive hp-FEM for the contact problem with Tresca friction in linear elasticity[END_REF] and convergence is shown for an elastic two-body problem in [START_REF] Schröder | Mixed finite element methods for two-body contact problems[END_REF]. 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 ∈ B H . Accordingly the discrete space for the Lagrange multiplier concerning friction is defined by the non-conforming ansatz

Λ t,H := {v ∈ (M q H ) d-1   ∀E ∈ B H : ∀x ∈ C q : v(F E (x)) ≤ 1 } =: span{ψ t i } m 2 i=1
and becomes conforming again for q = 0, 1. The discrete space Λ t,H is spanned by the

m 2 := (d -1)m 1 basis vectors  ψ n j e k , j = 1, . . . , m 1 , k = 1, . . . , d -1  with the kth unit vector e k in R d-1 .
In two space dimensions this leads to coinciding bases for Λ n,H and Λ t,H . Remark 1. The advantage of this choice is the given split of the contact pressure σ n (x) on Γ 1 C in its normal and tangential part σ n (x) = σ nn (x)n(x) + σ nt (x)t (x). The normal and tangential components of the contact pressure are directly calculated by

σ nn (x) = λ n (x)n(x) and σ nt (x) =  d-1 i=1 λ 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 [START_REF] Schröder | Mixed finite element methods for two-body contact problems[END_REF]. It is shown that if the inequality

2  m=1  h H -1 max{1, q} 2 p -1  θ m ≤ ε
holds for an ε > 0 sufficiently small and 1 + θ m -regularity of the solution u there exists an α ∈ R ≥0 independent of h, H, p and q such that the Brezzi-Babuska condition

α   (µ n,H , µ t,H )   -1/2 ≤ sup v h ∈V h ,∥v h ∥=1  b n (µ n,H , v h ) + b t (µ t,H , v h ) 
is fulfilled for all (µ n,H , µ t,H ) ∈ Λ n,H × Λ 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. [START_REF] Kleemann | Adaptive FEM für Mehrkörperkontaktprobleme[END_REF]. Eventually the discrete problem is to find

(u h , λ n,H , λ t,H ) ∈ V h × Λ n,H × Λ t,H with a(u h )(v h ) + b n (λ n,H , v h ) + b t (λ t,H , v h ) = f ext (v h ) b n (µ n,H -λ n,H , u h ) + b t (µ t,H -λ t,H , u h ) ≥ ⟨µ n,H -λ n,H , g⟩ (15)
for all v h ∈ V h , µ n,H ∈ Λ n,H and µ t,H ∈ Λ t,H .

Semi-smooth Newton methods

Following the work [START_REF] Hüeber | Discretization techniques and efficient algorithms for contact problems[END_REF] we transfer the presented active-set strategies for Mortar methods to mixed finite elements, which were introduced by Haslinger [START_REF] Haslinger | Mixed formulation of elliptic variational inequalities and its approximation[END_REF] for low order. We generalize the solving methods to higher-order discretizations that are proposed in [START_REF] Schröder | Mixed finite element methods for two-body contact problems[END_REF] and general friction laws depending on the contact forces.

Active-set strategy for contact

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

 Γ 1 C ([u h ] n -g)ψ H do ≤ 0 (16)  Γ 1 C λ n,H ψ H do ≥ 0 (17)   Γ 1 C λ n,H ψ H do    Γ 1 C ([u h ] n -g)ψ H do  = 0 (18)
for all ψ H ∈ Λ 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 ∈ R m 1 ×n m concerning the bases of V h and Λ H,n by (N 1 ) i j :=  Γ 1 C ψ n i (x)γ Γ 1 C (φ 1 j )(x) • n do, (N 2 ) i j :=  Γ 1 C ψ n i (x)γ Γ 2 C (φ 2 j )(Φ(x)) • n do
as well as the mass matrix

M ∈ R m 1 ×m 1 , M i j :=  Γ 1 C ψ n j ψ n i do
and the gap vector

ḡ ∈ R m 1 , ḡi =  Γ 1 C gψ n i 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 ūh -ḡ ≤ 0, M λn,H ≤ 0,  M λn,H  i (N ūh -ḡ) i = 0 ∀i = 1, . . . , m 1 . (19)
To reformulate [START_REF] Hager | Nonlinear complementarity functions for plasticity problems with frictional contact[END_REF] in terms of an equation we define the NCP function

C N : R n × R m 1 → R m 1 C N ( ūh , λn,H ) i := (-M λn,H ) i -max  0, (-M λn,H ) i + c n (N ūh -ḡ) i 
with a positive constant c n . Following the arguments in [9, Chapter 4] the discrete, weak contact conditions [START_REF] Hager | Nonlinear complementarity functions for plasticity problems with frictional contact[END_REF] can equivalently be expressed by

C N ( ūh , λn,H ) = 0. ( 20
)
With variations δu k h , δλ k n,H and characteristical functions χ i defined by

χ i :=  1,  -M λn,H + c n (N ūh -ḡ)  i > 0 0,  -M λn,H + c n (N ūh -ḡ)  i ≤ 0
the generalized derivative of the NCP function reads

C ′ N ( ūh , λn,H )( δu h , δλ n,H ) i = -χ i  -M δλ n,H + c n N δu h  i -  M δλ n,H  i .
We apply a semi-smooth Newton's method

C ′ N ( ūk-1 h , λk-1 n,H )( δu k h , δλ k n,H ) = -C N ( ūk-1 h , λk-1 n,H )
to solve [START_REF] Seitz | A semi-smooth newton method for orthotropic plasticity and frictional contact at finite strains[END_REF]. The new iterates are calculated by ūk h := ūk-1

A k n :=  i ∈ {1, . . . , m 1 }    -M λn,H + c n (N ūh -ḡ)  i > 0  ( 21
)
I k n :=  i ∈ {1, . . . , m 1 }    -M λn,H + c n (N ūh -ḡ)  i ≤ 0  . (22)
Then the iterate ( ūk h , λk n,H ) of Newton's method w.r.t. the contact conditions can equivalently be expressed by

(N ūk h ) i = ḡi , i ∈ A k n (M λk n,H ) i = 0, i ∈ I k n for all i = 1, . . . , m 1 .

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 analogon of the frictional constraints (10)- [START_REF] Koziara | Semismooth Newton method for frictional contact between pseudo-rigid bodies[END_REF] reads

∥σ nt (u 1 h )∥ ≤ s(σ nn (u 1 h )) ∥σ nt (u 1 h )∥ < s(σ nn (u 1 h )) ⇒ [u h ] t = 0 ∥σ nt (u 1 h )∥ = s(σ nn (u 1 h )) ⇒ ∃α ∈ R ≥0 : [u h ] t = ασ nt (u h ). (23)
Instead of these pointwise conditions [START_REF] Dostál | FETI based algorithms for contact problems: scalability, large displacements and 3d coulomb friction[END_REF] we use a weak formulation with test functions ψ n ∈ Λ n respectively ψ t ∈ Λ t :

 Γ 1 C   λ t,H   ψ n do ≤  Γ 1 C s(λ n,H )ψ n do  Γ 1 C   λ t,H   ψ n do <  Γ 1 C s(λ n,H )ψ n do ⇒  Γ 1 C [u h ] t ψ t do = 0  Γ 1 C   λ t,H   ψ n do =  Γ 1 C s(λ n,H )ψ n do ⇒ ∃α ∈ R ≥0 :  Γ 1 C [u h ] t ψ t do = α  Γ 1 C λ t,H ψ t do.
We replace the normal and tangential stresses σ nn (u h ) and σ nt (u h ) by the corresponding Lagrange multipliers λ n and λ t . The integrals

 Γ 1 C   λ t,H   ψ n i do and  Γ 1 C
λ t ψ t i do are approximated respectively integrated exactly by a (q + 1) d-1 -point gaussian quadrature with weights α i . The nodes of the Lagrangian basis of Λ n,H and Λ t,H are placed in these gaussian quadrature points

C q := { x1 , . . . , x(q+1) d-1 } on the reference element [-1, 1] d-1 . Define the matrix T ∈ R m 2 ×n with partition T = [ T1 T2 ] by  T1  i j := 1 ω i  Γ 1 C ψ t i (x) t (x) ⊤ γ Γ 1 C (φ 1 j )(x) do  T2  i j := 1 ω i  Γ 1 C ψ t i (x) t (x) ⊤ γ Γ 2 C (φ 2 j )(Φ(x)) do.
The weight ω i is given by ω

i := α i | det ∇ F E ( xi ) ⊤ ∇ F E ( xi )| ≥ 0 on E = supp(ψ t i ).
Let

s i := 1 ω i  Γ 1 C s(λ - n,H )ψ n i do (24)
with λ - n,H := min{0, λ n,H } be the algebraic representation of the positive friction bound s. The resulting discrete frictional constraints read

 λ t,H,I(i)   ≤ s i  λ t,H,I(i)   < s i ⇒  T ūh  I(i) = 0  λ t,H,I(i)   = s i ⇒ ∃α ∈ R ≥0 :  T ūh  I(i) = α λt,H,I(i) (25)
for the dimension-dependent index map

I(i) :=  {i}, d = 2 {2i -1, 2i}, d = 3.
In this sense for a vector u ∈ R m 2 the selection u I(i) is a scalar value in the case of two space dimensions and a vector in R 2 if d = 3. Following [9, Section 5] we express the conditions (25) by the NCP function

C T : R n × R m 1 ≤0 × R m 2 → R m 2 , C T ( ūh , λn,H , λt,H ) I(i) := max  s i ,    λt,H,I(i) + c t T ( ūh ) I(i)    λt,H,I(i) -s i  λt,H,I(i) + c t T ( ūh ) I(i) 
with c t > 0 and i ∈ {1, . . . , m 1 }. Analog arguments to the proof in [9, Theorem 5.1] lead to the equivalence of the equation C T ( ūh , λn,H , λt,H ) = 0 and the frictional constraints [START_REF] Krause | Monotone multigrid methods for Signorini's problem with friction[END_REF]. Defining the characteristical functions

χ A t,i :=  1,  λ t,H,I(i) + c t T ūh,I(i)   -s i > 0 0,  λ t,H,I(i) + c t T ūh,I(i)   -s i ≤ 0 and χ I t,i :=  1,  λ t,H,I(i) + c t T ūh,I(i)   -s i ≤ 0 0,  λ t,H,I(i) + c t T ūh,I(i)   -s i > 0
the generalized derivative for λn,H,i < 0 is given by

C ′ T ( ūh , λn,H , λt,H )( δu h , δλ n,H , δλ t,H ) I(i) = χ A t,i λt,H,I(i)  λt,H,I(i) + c t ( T ūh ) I(i)  ⊤ ∥ λt,H,I(i) + c t ( T ūh ) I(i) ∥ ( δλ t,H,I(i) + c t ( T δu h ) I(i) ) + δλ t,I(i) max  s i ,   ( λt,H,I(i) + c t ( T ūh ) I(i) )    -s i ( δλ t,I(i) + c t ( T δu h ) I(i) ) -s ′ i,λ n,H ( δλ n,H )( λt,H,I(i) + c t ( T ūh ) I(i) ) + χ I t,i s ′ i,λ n,H ( δλ n,H ) λt,I(i) .
We formulate the semi-smooth Newton method

C ′ T ( ūk-1 h , λk-1 n,H , λk-1 t,H )( δu h , δλ n,H , δλ t,H ) = -C T ( ūk-1 h , λk-1 n,H , λk-1 t,H ) (26)
with increments δu

k h = ūk h -ūk-1 h , δλ k n,H = λk n,H -λk-1 n,H and δλ k t,H = λk t,H -λk-1 t,H
. The Newton equation ( 26) is equivalent to

 Γ 1 C U ( ūk-1 h , λk-1 n,H , λk-1 t,H ) i λk t,H,I(i) ψ do +  Γ 1 C V ( ūk-1 h , λk-1 n,H , λk-1 t,H ) i λk n,H,i ψ do +  Γ 1 C W ( ūk-1 h , λk-1 n,H , λk-1 t,H ) i ( T ūk h ) I(i) ψ do =  Γ 1 C R( ūk-1 h , λk-1 n,H , λk-1 t,H ) i ψ do (27)
for all i ∈ {1, . . . , m 1 } and ψ ∈ Λ n,H . We define the dimension-dependent index map J ( j) =  j d-1



and separate the set of indices {1, . . . , m 2 } into slip indices with contact

A k t :=  j ∈ {1, . . . , m 2 }    i := J ( j) :    λk-1 t,H,I(i) + T ( ūk-1 h ) I(i)    -s k-1 i > 0, λk-1 n,H,i < 0  , (28)
stick indices in contact

I k t :=  j ∈ {1, . . . , m 2 }    i := J ( j) :    λk-1 t,H,I(i) + T ( ūk-1 h ) I(i)    -s k-1 i ≤ 0, λk-1 n,H,i < 0  , (29)
and indices without contact

I k tn :=  j ∈ {1, . . . , m 2 } | λk-1 n,H,J ( j) = 0  . ( 30
)
For j ∈ A k t and i := J ( j) we obtain the functions

U ( ūk-1 h , λk-1 n,H , λk-1 t,H ) i = λk-1 t,H,I(i) ( λk-1 t,H,I(i) + c t ( T ūk-1 h ) I(i) ) ⊤ ∥( λk-1 t,H,I(i) + c t ( T ūk-1 h ) I(i) )∥ + I (d-1)×(d-1)  ∥( λk-1 t,H,I(i) + c t ( T ūk-1 h ) I(i) )∥ -s k-1 i  , V ( ūk-1 h , λk-1 n,H , λk-1 t,H ) i = -(s ′ i,λ n,H ) k-1 ( λk-1 t,H,I(i) + c T ( ūk-1 h ) I(i) ), W ( ūk-1 h , λk-1 n,H , λk-1 t,H ) i = c t λk-1 t,H,I(i) ( λk-1 t,H,I(i) + c t ( T ūk-1 h ) I(i) ) ⊤ ∥ λk-1 t,H,I(i) + c t ( T ūk-1 h ) I(i) ∥ -c t I (d-1)×(d-1) s k-1 i and R( ūk-1 h , λk-1 n,H , λk-1 t,H ) i = c t λk-1 t,H,I(i)  λk-1 t,H,I(i) + c t ( T ūk-1 h ) I(i)  ⊤ ∥ λk-1 t,H,I(i) + c t ( T ūk-1 h ) I(i) ∥ • ( λk-1 t,H,I(i) + c t ( T ūk-1 h ) I(i) ) -(s ′ i,λ n,H ) k-1 ( λk-1 n,i )( λk-1 t,H,I(i) + c t ( T ūk-1 h ) I(i) ).
In the slip case j ∈ I k t the functions W i = 0,

U ( ūk-1 h , λk-1 n,H , λk-1 t,H ) i = 1 ω i I (d-1)×(d-1) , V ( ūk-1 h , λk-1 n,H , λk-1 t,H ) i = ( T ūk-1 ) I(i) (s ′ i,λ n,H ) k-1 s k-1 i and R( ūk-1 h , λk-1 n,H , λk-1 t,H ) i = ( T ūk-1 ) I(i) (s ′ i,λ n,H ) k-1 λk-1 n,H,i s k-1 i
determine the Newton equation ( 26). In the limit case j ∈ I k tn of no contact we choose λt,H,I(i) = 0 and therefore vanishing U i , V i and R i as well as W i = I (d-1)×(d-1) . We define the matrices

U = [U 1 U 2 ] ∈ R m 2 ×n , V ∈ R m 2 ×m 1 , W ∈ R m 2 ×m 2 and the vector r ∈ R m 2 by (U 1 ) ji :=  Γ 1 C U ( ūk-1 h , λk-1 n,H , λk-1 t,H ) J ( j) t δ (x) ⊤ γ Γ 1 C (φ 1 i )(x)ψ t j (x)dσ (U 2 ) ji :=  Γ 1 C U ( ūk-1 h , λk-1 n,H , λk-1 t,H ) J ( j) t δ (x) ⊤ γ Γ 1 C (φ 2 i )(x)ψ t j (Φ(x))dσ V ji :=  Γ 1 C V ( ūk-1 h , λk-1 n,H , λk-1 t,H ) J ( j) ψ n i ψ t j dσ W ji :=  Γ 1 C W ( ūk-1 h , λk-1 n,H , λk-1 t,H ) J ( j) ψ t i ψ t j dσ r j :=  Γ 1 C R( T ūk-1 h , λk-1 n,H , λk-1 t,H ) J ( j) ψ t j dσ.
Then the Newton equation ( 27) is equivalent to the vector-valued equation

U ūk h + V λk n,H + W λk t,H = r .

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 T let the unscaled matrix T = [T 1 T 2 ] ∈ R m 2 ×n be defined by

(T 1 ) i j :=  Γ 1 C ψ t i (x) t (x) ⊤ γ Γ 1 C (φ 1 j )(x) do, (T 2 ) i j :=  Γ 1 C ψ t i (x) t (x) ⊤ γ Γ 2 C (φ 2 j )(Φ(x)) do,
as well as the block-diagonal matrix

K k = diag(K k 1 , K k 2 ) =  a ′ (u k-1 h )(φ m j , φ m i )  i, j=1,...,n
by the Frechét derivative of the semi-linearform in direction u k-1 h . Define the vector

L k (u k-1 h ) = [L k 1 L k 2 ] ⊤ =  a(u k-1 h )(φ m i )  i=1,...,n
and the Newton right hand side of the plastic problem

f k m := K k m ūk-1 h,m -L k m (u k-1 h ) + fm .
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

                       K k 1 0 N ⊤ 1,A k n N ⊤ 1,I k n T ⊤ 1,A k t T ⊤ 1,I k t T ⊤ 1,I k tn 0 K k 2 N ⊤ 2,A k n N ⊤ 2,I k n T ⊤ 2,A k t T ⊤ 2,I k t T ⊤ 2,I k tn N 1,A k n N 2,A k n 0 0 0 0 0 0 0 M A k n I k n M I k n I k n 0 0 0 U 1,A k t U 2,A k t V A k t A k n V A k t I k n W A k t A k t W A k t I k t W A k t I k tn T1,I k t T2,I k t V I k t A k n V I k t I k n 0 0 0 0 0 0 0 0 0 I                                               ūk h,1 ūk h,2 λk n,H,A k n λk n,H,I k n λk t,H,A k t λk t,H,I k t λk t,H,I k tn                        =                        f1 f2 ḡA k n 0 rA k t rI k t 0                        (31)
with the following notation. For a matrix A ∈ R k×l the submatrix

A IJ :=  A i, j  i∈I, j∈J ∈ R |I|×|J |
consists of the rows and columns of A, whose indices belong to the sets I respectively J of indices. If only rows of I are selected the notation is

A I :=  A i,•  i∈I ∈ R |I|×l .
The linear system corresponds to a saddle point problem of the structure

 K k B ⊤ 1 B k 2 -C k   ūk h λk  = b ( 32
)
with partitions

λk :=   λk n,H λk t,H   , B ⊤ 1 :=   N ⊤ 1 T ⊤ 1 N ⊤ 2 T ⊤ 2   , B k 2 :=           N 1,A k n N 2,A k n 0 0 U 1,A k t U 2,A k t T1,I k t T2,I k t 0 0           and -C k :=           0 0 0 0 0 M A k n I k n M I k n I k n 0 0 0 V A k t A k n V A k t I k n W A k t A k t W A k t I k t W A k t I k tn V I k t A k n V I k t I k n 0 0 0 0 0 0 0 I           , b :=                   f1 f2 ḡA k n 0 rA k t rI k t 0                  
.

In general the system matrix of the problem ( 32) is asymmetric. The matrices B 1 and B k 2 are nonquadratic whereas the quadratic matrix C k might be singular. To solve the discrete contact problem [START_REF] Hüeber | A priori error estimates and an inexact primal-dual active set strategy for linear and quadratic finite elements applied to multibody contact problems[END_REF] we have to solve a sequence of saddle point problems [START_REF] Rannacher | A posteriori error estimation and mesh adaptation for finite element models in elasto-plasticity[END_REF]. This semi-smooth Newton method is summarized in Algorithm 1.

res plast := 2  m=1 ∥K k m ūk h,m -L k m ( ūk h,m ) + f m ∥ C k N := C N ( ūk h , λk n,H ) C k T := C T ( ūk h , λk n,H , λk t,H ) 8 if ∥C k N ∥ + ∥C k T ∥ + res plast ≤ tol then stop 9 perform damping if k > 2 : ũ0 m := ūk h,m
, res 0 := res plast , ω 0 := 0.5 for j = 0 to n damp do

ũ j m := ω j ūk h,m + (1 -ω j ) ūk-1 h,m res j :=  2 m=1 ∥K k m ũ j m -L k m ( ũ j m ) + f m ∥ if res j ≤ res j-1 then ω j := 0.5 ω j-1 else ūk h := ũ j-1 m break
In every Newton step the linear system ( 32) is solved to a tolerance tol k . We apply an inexact strategy and start with a relatively high value of tol k 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

P k =  K k B ⊤ 1 0 S k  . (33)
The matrix S k is defined as the generalized Schur complement matrix

S k := -(B k 2 (K k ) -1 B ⊤ 1 + C k ).
This approach was proposed in the work [START_REF] Benzi | Numerical solution of saddle point problems[END_REF]Section 10.1.2]. We choose an approximation Ŝk to S k implicitly by the action of ( Ŝk ) -1 on a given vector v. The Schur complement matrix is asymmetric and not given explicitly, so the approximation Ŝk to S k is defined by the solution with GMRES to a given tolerance. In an analogous manner we define an approximation K k to K k . The approximation ( K k ) -1 has to be calculated once per GMRES step in the solving process with Ŝ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 ( 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

P k S =  diag  -(B k 2  diag(K k )  -1 B ⊤ 1 + C k )  -1 . (34)
The explicit form of

 P k  -1 reads  P k  -1 =  ( K k ) -1 0 0 I   I B ⊤ 1 0 -I   I 0 0 -( Ŝk ) -1  .
Basically the effort of applying P k S 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 k N , C k T and res 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 [START_REF] Nocedal | Numerical Optimization[END_REF] beginning in the step k = 3 of the algorithm.

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.

Motivation

We consider two flat rectangular elastic bodies

Ω 1 = (0, 20) × (0, 2) and Ω 2 = (18, 38) × (-2.2, -0.2) that are clamped at the sides Γ 1 0 = {0} × [0, 2] and Γ 2 0 = {38} × [-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 ν 1 = ν 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 [START_REF] Braess | Finite elements in solid mechanics[END_REF]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. [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF]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).

Two dimensional example

We consider the contact of an elasto-plastic and a linear-elastic body, which are represented by Ω 1 = (0, 1)×(0, 1) and Ω 2 = (0.91, 1.91) × (0, 1), cf. C that fulfills H = 2h. The Young's moduli E 1 = 10 4 , E 2 = 10 3 and Poisson's ratios ν 1 = 0.25, ν 2 = 0.22 define the compliance tensors A 1 and A 2 . We choose hardening parameters γ 1 iso = 0.01, γ 2 iso = 1.0 and yield stresses σ 1 0 = 10 2 , σ 2 0 = 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 |σ (u 1 h ) D | F and the portion of plastified quadrature points are shown in Fig. 4. We consider a friction law after Betten [START_REF] Horwatitsch | Entwicklung eines Reibmodells für hohe Temperaturen und hohe Umformgrade[END_REF]Section 2.4] given by the (a) p = 1, q = 0 and p = 2, q = 1.

(b) p = 1, q = 0 and p = 5, q = 4.

(c) p = 1, q = 0 and p = 2, q = 1.

(d) p = 1, q = 0 and p = 5, q = 4. 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 1 Comparison of preconditioners with p = 1, q = 0 on a uniform refined mesh, n = 288, m 1 = m 2 = 4, tol = 10 -9 .

Step active indices the quantity C k N tends to zero very fast. After some iterations in which the displacements have to be damped (α = 0.5) the residuals C k T and res 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 [START_REF] Benzi | Numerical solution of saddle point problems[END_REF]. In the third version the generalized Schur complement matrix Ŝk within the evaluation of P k is additionally preconditioned by P k S given in [START_REF] Hüeber | A primal-dual active set strategy for non-linear multibody contact problems[END_REF]. The number of outer GMRES steps is named as iter, whereas iter S labels the accumulated number of GMRES iterations that are performed to solve with the Schur complement. Every evaluation of Ŝk includes a solution with the block-diagonal matrix K k . Due to the chosen friction law and the start value λ 0 n,H = 0 the friction bound and therefore λ 1 t,H 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 S 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 necessiates 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 [START_REF] Hüeber | A priori error estimates and an inexact primal-dual active set strategy for linear and quadratic finite elements applied to multibody contact problems[END_REF][START_REF] Wohlmuth | An abstract framework for a priori estimates for contact problems in 3D with quadratic finite elements[END_REF][START_REF] Coorevits | Mixed finite element methods for unilateral problems: Convergence analysis and numerical studies[END_REF] and linear-elastic problems with volumetric inequality constraints [START_REF] Bonfils | Treating volumetric inequality constraint in a continuum media with a coupled X-FEM/level-set strategy[END_REF]. The full benefit of optimal convergence for higher-order methods can only be achieved by using adaptivity [START_REF] Kleemann | Adaptive FEM für Mehrkörperkontaktprobleme[END_REF][START_REF] Dörsek | Adaptive hp-FEM for the contact problem with Tresca friction in linear elasticity[END_REF].

Threedimensional example

As an application from mechanical engineering we consider a two-body contact problem that models a grinding process. The first body Ω 1 represents a cylindrical mounted point with diameter 2.2. Its spindle is clamped at the end. A hexahedral body Ω 2 = (-2, 3) × (1.1, 3.1) × (-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(λ n ) := min{µ|λ n |, Y 0 }
with constants µ = 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 , ν 1 = ν 2 = 0.3, σ 1 0 = 3, σ 1 0 = 2 and γ 1 iso = γ 2 iso = 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 k T appears more unstable. This behavior indicates the sensitivity of the method concerning the choice of the constant c t .

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 [START_REF] Hager | Nonlinear complementarity functions for plasticity problems with frictional contact[END_REF] 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 [START_REF] Dickopf | Efficient simulation of multi-body contact problems on complex geometries: A flexible decomposition approach using constrained minimization[END_REF] and separation of two bodies [START_REF] Wriggers | Computational Contact Mechanics[END_REF] as well as the development of appropriate adaptive strategies are planed.

within the Collaborative Research Center (SFB) 708 "3D-Surface Engineering of Tools for Sheet Metal Forming -Manufacturing, Modeling, Machining".

  ) ensures that outer and inner forces are balanced. The deviatoric part of a tensor τ is denoted by τ D := τ -1 d tr(τ )I d×d and | • | F indicates the Frobenius norm. Defining the yield function F m,iso (τ, η) = |τ D | F -(σ m 0 + γ m iso η) with the yield stress σ m 0 and the isotropic hardening parameter γ m iso 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 Γ m D ⊂ Γ m with positive measure. These Dirichlet boundary conditions are described in Eq. (4). At parts Γ m C ⊂ Γ m the two bodies may come into contact. We assume that the open set Γ m C fulfills Γ m C Γ m \Γ m D . On the remaining part of the boundaries Γ m N

  n,H := λk-1 n,H + δλ k n,H . Let active and inactive indices in Newton step k be determined in the following way

Algorithm 1 : 3 calculate friction bound ( 24 ) 4 calculate 5 assemble

 132445 Semi-smooth Newton method for k = 0 to maxIter do 1 assemble linearization K k and right hand side f k 2 calculate active and inactive sets A k n , I k n for contact (21, 22) active and inactive sets A k t , I k t , I k tn for friction (28, 29, 30) Newton matrices U k , V k , W k and right hand side r k 6 solve linear system (32) to tolerance tol k 7 calculate residuals

Fig. 1 .

 1 Fig. 1. Comparison of deformations for different polynomial degrees and SRI.

Fig. 2 .

 2 Fig. 2. Setting of the 2D two-body contact problem.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Deformations in xand y-direction.

Fig. 2 (

 2 a). The body Ω 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 Γ 1 D = {0} × [0, 1] and Γ 2 D = {1.91} × [0, 1]. The two bodies may come into contact on the contact boundaries Γ 1 C = {1} × [0, 1] respectively Γ 2 C = {0.91} × [0, 1], which is caused by the overlapping of Ω 1 and Ω 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 Γ 1

Fig. 5 .Fig. 6 . 0  p  1 p

 5601 Fig. 5. Lagrange multipliers λ n and λ t for different polynomial degrees.

Fig. 7 .

 7 Fig. 7. Convergence of the error in the energy norm with uniform mesh refinement.

2 .Fig. 8 .

 28 Fig. 8. Overview of the 3D example.

Fig. 9 .

 9 Fig. 9. Plastic variables of the mounted point Ω 1 .

Fig. 10 .Fig. 11 .Fig. 12 .

 101112 Fig. 10. Resulting mechanical variables in the contact zone of workpiece Ω 2 .

L: 1 , m 2 :

 12 Lagrange functional F m,T , F E : Transformation on T ∈ T h,m resp. E ∈ B H C N , C T : NCP functions concerning normal contact and friction χ: Characteristical function I, J : Index maps U , V , W , R: Functions defining the Newton equation Sets Ω m : Deformable body m = 1, 2 Γ m , Γ m D , Γ m N , Γ m C : Boundary of Ω m , Dirichlet, Neumann and contact boundary T h,m , B h,m : Triangulations of Ω m and Γ m C C q : Set of q + 1 Gauss-quadrature points A k n , I k n : Sets of active/inactive indices concerning geometrical contact A k t , I k t , I k tn : 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 δ : Outer normal vector on Γ m and generalized normal vector on Γ 1 C ḡ: Gap vector T, T : Unweighted and weighted tangential coupling matrix S: Schur complement matrix ūh , λn,H , λt,H : Vector-valued representations of u h , λ n,H and λ t,H s: Friction bound U, V, W, r : Vector-valued representations of U , V , W , R P k : Block triangular preconditioner P k S : Preconditioner of Schur complement matrix Ŝ, K : Approximations to S resp. K by solving with GMRES Constants n, n m : Dimensions of V h resp. Dimensions of Λ n,H resp. Λ t,H h, H : Mesh sizes of Ω m and Γ m C 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

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"

Glossary

× (H D (Ω 2 )) d V p h,m : Set of ansatz functions of polynomial degree p on body m with meshsize h V h :

Λ n , Λ t : Continuous normal and tangential Lagrange ansatz spaces S r l : Polynomial tensor product space of order l on [-1, 1] l Λ n,H , Λ t,H : Admissible set for discrete Lagrange multipliers concerning contact and friction