Semi{Smooth Newton Method for Solving 2D Contact Problems with Tresca and Coulomb Friction

The contribution deals with contact problems for two elastic bodies with friction. After the description of the problem we present its discretization based on linear or bilinear finite elements. The semi– smooth Newton method is used to find the solution, from which we derive active sets algorithms. Finally, we arrive at the globally convergent dual implementation of the algorithms in terms of the Langrange multipliers for the Tresca problem. Numerical experiments conclude the paper.


Introduction
Nowadays, solving contact problems with friction counts among very challenging tasks in mechanics and is of crucial importance in various practical applications.As an example we can state biological processes, the design of machines and transportation systems, and metal forming.This is a strong motivation for the development of methods allowing a reliable and fast simulation, i.e., to implement robust numerical solvers.
As for the history, in 1933, Signorini was one of the first savants, who considered the frictionless contact of an elastic body with a rigid foundation.These problems are often named after him.
Contact problems with friction are inherently nonlinear.This fact complicates the teoretical analysis and the proposal of an efficient numerical solver.General information about contact problems one can find, e.g., in [17].The main difficulties of contact problems are two nonlinear conditions: the non-penetration of bodies and friction law.Predominant friction laws are called Coulomb or Tresca.
Newton-type methods for solving contact problems with friction were used already in 1991, see [1].In more recent contributions is shown that the performance of generalized Newton-like methods is superior to interior point methods, see [6].Semi-smooth Newton method in finite dimensions appears already in 1993, see [22], and also more recently in infinite dimensions.For that, the concept of the slant differentiability was introduced; see [2] and the references given therein.The proof of the superlinear convergence of the semi-smooth Newton method is stated in [23].
This paper proposes effective solvers for contact problems with friction using the semi-smooth Newton method.Moreover, a new idea of the globalization strategy that is (surprisingly) similar to the MPRGP algorithm, see [3], is presented.Let us introduce conventions that we use throughout the whole text.For any non-empty set of indices I and a matrix M ∈ R m×n , we denote by M I a submatrix of M with the rows given by I.

Formulation
Let us consider two homogeneous isotropic elastic bodies represented by bounded domains Ω k ⊂ R , where we consider three contact conditions: the non-penetration of bodies (the unilateral contact law), the transmission of contact stresses, and the Coulomb friction law.Elastic properties of Ω k are described by the Lame constants λ k , µ k > 0. Finally, assume the volume forces Fig. 1: Geometry of the bodies.
Our aim is to find an equilibrium state of Ω 1 and Ω 2 .By the solution of this problem we mean displacement vector fields , satisfying the equilibrium equations and the Dirichlet and Neumann conditions: where σ k := σ(u k ) is the stress tensor in Ω k and n k stands for the unit outward normal vector to ∂Ω k , k = 1, 2. The tress tensor is related to the linearized strain tensor ε k := ε(u k ) = 1/2(∇u k + ∇ u k ) by the Hooke law for linear isotropic materials: where "tr" denotes the trace of matrices and I ∈ R 2×2 is the identity matrix.Let us note that λ k , µ k are given by the Young modulus E k > 0 and the Poisson ratio ν k ∈ (0, 0.5) as To formulate the contact conditions we introduce a predefined one-to-one transfer mapping χ : by means of which we define the initial distance between the bodies at x ∈ Γ 1 c as d(x) := χ(x) − x and the critical direction ν(x . stands for the Euclidean norm in R 2 .The unilateral contact law reads as follows: where u ν (x) := (u 1 (x) − u 2 (χ(x))) ν(x) stands for the relative normal contact displacement and σ ν (x) := ν(x) σ 1 (x)ν(x) is the normal contact stress at x ∈ Γ 1 c .Let t := t(x) be an unit vector orthogonal to ν := ν(x) so that the pair {ν, t} is a local orthonormal basis in R 2 with the origin at x ∈ Γ 1 c .By F := F(x) we denote given positive coefficient of friction.The Coulomb friction law for x ∈ Γ 1 c reads as follows: where u t (x) := (u 1 (x) − u 2 (χ(x))) t(x) and σ t (x) := t(x) σ 1 (x)ν(x) are the relative tangential contact displacement and the tangential contact stress at x ∈ Γ 1 c , respectively.Finally we require the transmission of the contact stresses: Note that the satisfaction of ( 4) is inherently hidden in the weak form of the problem.Therefore, the discrete version of the problem does not include this condition explicitely.The formulation of the contact problem with Coulomb friction consists in finding a pair u := (u 1 , u 2 ) of the displacement fields u k in Ω k , k = 1, 2, satisfying (1), ( 2), (3), and (4).Let us note that this problem results in an implicit variational inequality of the elliptic type.One of possible ways how to obtain the solution is the fixed-point approach.To this end we replace the slip bounds in (3) by an apriori given non-negative function g ∈ L 2 + (Γ 1 c ), g ≥ 0, so that for x ∈ Γ 1 c : Now, (1), (2), (5), and (4) describes the (auxiliary) contact problem with Tresca friction.It is well-known that this problem can be represented by the variational inequality of the second type for which there is the unique solution for each g ∈ L 2 + (Γ 1 c ).Let us consider the mapping Ψ : where σ ν (g) := σ ν (u(g)) is the normal contact stress corresponding to the solution of the contact problem with Tresca friction with given g.It is easily seen that the fixed point g * of Ψ defines the solution to the contact problem with Coulomb friction, i.e., the point g * such that Ψ(g * ) = g * .The natural way how to compute fixed points is the method of successive approximations.It is based on the following iterative scheme: c 2013 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING These iterations converge when Ψ is contractive that is guaranteed by sufficiently small coefficient of friction F. Finally, note that each evaluation of the mapping Ψ in (6) requires to solve the contact problem with Tresca friction.
When concerning the discrete problem with Tresca friction it is known that there exists a unique solution, see [3].The situation is different for the discrete contact problem with Coulomb friction.In this case, the solution exists, but the uniqueness is guaranteed only for F sufficiently small, see [8], that we consider in this work.The sufficient smallness is mesh dependent.Some of these problems concerning the non-unique solutions are discussed in [9] using continuation in contact problems with Coulomb friction.

Discrete Problems
We describe the finite element approximation based on the linear or bilinear finite elements.Let Let N k cont ⊂ N k be the subset of the contact nodes of Ω k , i.e., the nodes lying on Γ k c \ Γ k u , and let m k be their number.Moreover, we denote n = n 1 + n 2 and m = m 1 .The stiffness matrix and the load vector corresponding to Ω k are denoted by With the contact nodes of Ω 1 we associate two matrices N 1 , T 1 ∈ R m×n 1 so that the ith row of N 1 and T 1 contains at most two non-zero entries given by components of ν(x 1 i ) and t(x 1 i ), respectively, The ith row of N 2 and T 2 contains at most four nonzero entries given by the components of . By d ∈ R m we denote the vector whose entries are given by d i = d(x 1 i ).Finally, F stands for the diagonal matrix cont .Note that the relative tangential displacement at the contact node The discrete problems use two Lagrange multipliers λ ν , λ t ∈ R m that approximates −σ ν and −σ t at the contact nodes, respectively.The discrete contact problem with Coulomb friction reads as follows: for i = 1, . . ., m.
In order to obtain the discrete contact problem with Tresca friction, we replace F i λ ν,i in (9) by the entries g i of g ∈ R m , g ≥ 0, that approximate values of g from ( 5) at the nodes of N 1 cont .Next we use the equivalent formulation of the previous problems as the systems of non-smooth equations.To this end we introduce the projection mappings: respectively, where the inequalities and the absolute value are understood componentwisely and r ∈ R m + .The definitions of the components of P R m + , P Λ(r) are based on the max-function in R 1 : Denote y := (u , λ ν , λ t ) ∈ R 2n+2m and consider a parameter ρ > 0. The discrete contact problem with Coulomb friction is equivalent to the equation: where G : R 2n+2m → R 2n+2m is defined by: The discrete contact problem with Tresca friction is equivalent to the equation: where G g : R 2n+2m → R 2n+2m is defined by: Note that both G and G g are piecewise smooth functions with the non-smooth points determined by the max-function in (10) and (11).

Algorithms
There are two principal possibilities how to derive an algorithm: (i) solving directly the discrete contact problem with Coulomb friction; (ii) solving the discrete contact problem with Tresca friction and combine it with the method of successive approximations (6).In both cases we will use a variant of the Newton method that is based on the iterative scheme: ), (14) for k = 1, 2, . . ., where F : R 2n+2m → R 2n+2m is a sufficiently smooth function and J F (y) approximates the Jacobi matrix to F at y ∈ R 2n+2m .It is well known that the sequence {y (k) } generated by ( 14) converges superlinearly to the solution of F (y) = 0, if F is slantly differentiable and if the initial iteration y (0) ∈ R 2n+2m is sufficiently close to the solution, see [2].The slanting function F o , see [23], is used as J F and ( 14) is called the semi-smooth Newton method (SSNM).As the max-function is slantly differentiable, we can apply the SSNM for solving (12) and ( 13).
We will show how to derive the slanting function The idea is based on introducing the active and inactive sets that determine the smooth pieces of G or G g .These smooth pieces are differentiable so that the standard differential rules can be used.Such implementation of the SSNM is equivalent to an active set algorithm.The distinction between G o and G o g consists in the fact that G o g may be symmetrized.
Before discussing both algorithms we introduce notations.Let M = {1, 2, . . ., m} be the set of all indices and let y = (u , λ ν , λ t ) ∈ R 2n+2m be given.The active set A ν := A ν (y) corresponding to the nonpenetration condition is defined by: The respective inactive set is its complement I ν := I ν (y) = M \ A ν (y).For S ⊆ M we introduce the diagonal matrix D S = diag(s 1 , . . ., s m ) ∈ R m×m with:

Non-Symmetric Case
First we propose the algorithm for direct solving discrete contact problems with Coulomb friction.We introduce two inactive sets I + t := I + t (y), I − t := I − t (y) corresponding to the conditions of Coulomb friction: The respective active set is the complement A t := A t (y) = M \ (I + t (y) ∪ I − t (y)).Due to the definition of the projections (10) and (11), we get the following expression for G at y: After differentiating the right-hand side, we arrive at the slanting function: and We get the active-set implementation of the SSNM.
) ) define the active and inactive sets: (2.) Solve: The most expensive part of each iteration is solving the block linear systems in step (2.).These systems correspond to the finite element approximation of the linear elasticity problem with the following mixed conditions on Γ 1 c : the relative normal contact displacements, the zero relative normal contact stresses, and the zero relative tangential contact displacements at the contact nodes with indices of A ν , I ν , and A t , respectively; the achievement of the slip bounds λ t,i = −F i λ ν,i and λ t,i = F i λ ν,i at the contact nodes with indices of I − t and I + t , respectively.Note that ρ is discarded from the system in step (2.).Indeed, if this system is solved exactly, then ρ would not play any role in the definitions of the active and inactive sets (except of the first iteration).

Symmetric case
To propose the algorithm for solving the discrete contact problem with Tresca friction we define the inactive sets I + t := I + t (y), I − t := I − t (y) as follows: The respective active set is again the complement.Using (10) and (11), we get: The slanting function G o g to G g reads as follows: The implementation of the SSNM uses the fact that components of the solutions of the inner linear systems corresponding to the indices of I ν , I + t and I − t may be found without any complication.Eliminating them, we arrive at the reduced system with the symmetric saddle-point matrix for the remaining components.

• Algorithm ActiveSetTresca
This algorithm differs from the previous one Ac-tiveSetCoulomb only in step (2.): (2.) Solve: Here, we use the index set as the multiindex.In the case of a matrix, it represents its submatrix composed by rows with indices of the index set.The system in step (2.) corresponds again to the finite element approximation of a mixed problem on Γ 1 c .Comparing with that one in the previous section, the only distinction consists in the interpretation of conditions corresponding to the indices of I − t and I + t .

Implementation
In this section, we interpret step (2.) of algorithm ActiveSetTresca as the minimization of a strictly quadratic objective function in terms of λ.Therefore, the conjugate gradient method can be used as the inner solver.First of all we introduce notation.
Let us combine the active and inactive sets defined by ( 15), ( 18) and ( 19) in A and I for y = (u , λ ν , λ t ) as follows: A = {i| i ∈ A ν }∪{i+m| i ∈ A t } and I = {1, 2, . . ., 2m} \ A. Denote: where The gradient to q at λ reads as: The linear system in step (2.) of algorithm Ac-tiveSetTresca takes the form: where . The remaining components of λ (k) are given by: Lemma 1: 1) .Then λ (k) is fully determined by (21) and A is the minimizer to the problem: where The remaining components of λ (k) are given by (21) and Proof : The statement (i) is straightforward.Let us assume A = ∅.The first block equation in (20) gives A ) that enables us to eliminate the unknown u (k) .The resulting linear system in terms of λ (k) reads as: Note that A AA is the symmetric, positive definite so that (23) and the minimization problem (22) result in the same solution.The lemma is proved.
The next lemma shows how to define the active and inactive sets without the knowledge of u (k) .

Globalization Strategy
The standard result on convergence of the semi-smooth Newton method [2,23] requires an initial iteration "sufficiently close" to the solution.Unfortunately, it is impractical from the point of view of computations.To overcome this drawback, it is proposed by Jungho Lee in [19] to find an initial guess of the solution by projecting an auxiliary solution of an appropriate unconstrained problem, i.e., the problem without the nonpenetration and friction conditions.
In this section, we comment briefly an other globalization strategy.The main idea consists in inexact solving of inner linear systems using several steps of the conjugate gradient method.Indeed, our implementation of the algorithm ActiveSetTresca may be interpreted as the restarted conjugate gradient method, in which all inner iterations connected to one sequence {λ (l) } give the decreasing sequence {q(λ (l) )}.Moreover, we adopt from [3,22] the criterion for the adaptive (inner) precision control.Therefore, we arrive at the final implementation of the inexact semi-smooth Newton method that is closely related to the MPRGP algorithm of [3].The main distinction consists in definitions of the active and inactive sets that are given in our case by (15), (18) and (19).Nevertheless, the same convergence result may be proved as for the MPRGP that is valid for each initial iteration λ 0 from an apriori known feasible set R m + × Λ(g).Note that the MPRGP uses the reduced gradient that is defined by the steplength α whose value is bounded by a multiple of the inverse to the largest eigenvalue of A. In our implementation, the role of α plays the parameter ρ from the definitions of the active and inactive sets (15), ( 18) and (19).

Numerical Experiments
Numerical experiments have been computed by Matlab R2009a using computer with RAM 4GB, Intel Core i5 (2.53GHz) 460M.
Let us consider two plane elastic bodies: and Ω 2 = (0, 3) × (0, 1), made of an isotropic, homogeneous material characterized by the Young modulus 2.119 × 10 11 and the Poisson ratio 0, 277 (steel).The decompositions of ∂Ω 1 and ∂Ω 2 are as follows: The given function g representing the slip bound is equal to 1, 7 × 10 7 .The volume forces vanish for both bodies.The non-vanishing surface tractions 2 ) act on Γ 1 p so that: where p 1 2,L = −7 × 10 7 , p 1 2,R = −1/3 × 10 7 , p 1 2,B = 4 × 10 6 , and p 1 2,U = 1, 8 × 10 7 .In tables below, we compare the performance of various implementations of the algorithm ActiveSetTresca with the MPRGP for different numbers of dofs n and m.We report the number n K of matrix-vector multiplications by K −1 and the number k max of the outer (Newton) iterations in the case of the algorithm Ac-tiveSetTresca.The final terminating precision (ε λ ) is set 10 −6 in all cases.
In Tab. 1, we present results of numerical tests without the globalization strategy.The column labeled MPRGP is computed (by the MPRGP algorithm) with α = 2α −1 max , where α max denotes the largest eigenvalue of A. The column labeled SNM1 is obtained by exact solving the inner linear systems using the fixed inner precision 10 −12 and with ρ = 1.
In the column SNM2, we show the results, when the inner linear systems were solved inexactly using the heuristic adaptive (inner) precision control as proposed in [21] and with ρ = 1 again.
Finally in the column SNM3, we use the same adaptive precision control but we set ρ = α −1 max .One can see that the last choice leads to the most efficient behavior of the algorithm ActiveSetTresca.In the column SNM4 of Tab. 2 we report the results of computations with the globalization strategy for ρ = 20α −1 max .One can observe that the efficiency in this case is certain compromise between SSNM4 and MPRGP.

Conclusion
In the present contribution, we have investigated a two-body 2D contact problem with Coulomb and Tresca friction.Primarily, we have designed and implement an effective solver for these problems in Matlab that is based on the semi-smooth Newton method and the active set strategy.Moreover, we have arrived at the globally convergent dual implementation of the algorithm.Our numerical experiments indicate better results than MPRGP.

2
with sufficiently smooth boundaries ∂Ω k , k = 1, 2. Each boundary consists of three disjoint parts Γ k u , Γ k p , and Γ k c open in ∂Ω k so that ∂Ω k = Γ k u ∪Γ k p ∪Γ k c and Γ k u = ∅; see Fig. 1.The zero displacements are prescribed on Γ k u while surface tractions p k ∈ (L 2 (Γ k p )) 2 act on Γ k p .The bodies may get into contact on the contact interface given by Γ 1 c and Γ 2 c and λ t := λ (k) t .(4.) Set k := k + 1 and go to step (1).