On the solution of contact problems with Tresca friction by the semismooth* Newton method

An equilibrium of a linear elastic body subject to loading and satisfying the friction and contact conditions can be described by a variational inequality of the second kind and the respective discrete model attains the form of a generalized equation. To its numerical solution we apply the semismooth* Newton method by Gfrerer and Outrata (2019) in which, in contrast to most available Newton-type methods for inclusions, one approximates not only the single-valued but also the multi-valued part. This is performed on the basis of limiting (Morduchovich) coderivative. In our case of the Tresca friction, the multi-valued part amounts to the subdifferential of a convex function generated by the friction and contact conditions. The full 3D discrete problem is then reduced to the contact boundary. Implementation details of the semismooth* Newton method are provided and numerical tests demonstrate its superlinear convergence and mesh independence.



page 4


Lifted contact dynamics for efficient direct optimal control of rigid body systems with contacts

We propose a novel and efficient lifting approach for the direct optimal...

Convergence of a damped Newton's method for discrete Monge-Ampere functions with a prescribed asymptotic cone

We prove the convergence of a damped Newton's method for the nonlinear s...

Newton Type Methods for solving a Hasegawa-Mima Plasma Model

In [1], the non-linear space-time Hasegawa-Mima plasma equation is formu...

Rapid multi-component phase-split calculations using volume functions and reduction methods

We present a new family of fast and robust methods for the calculation o...

Total Variation Diminishing (TVD) method for Elastohydrodynamic Lubrication (EHL) problem on Parallel Computers

In this article, we offer a novel numerical approach for the solution of...

A Reverse Augmented Constraint preconditioner for Lagrange multiplier methods in contact mechanics

Frictional contact is one of the most challenging problems in computatio...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In [3] the authors developed a new, so-called semismooth Newton-type method for the numerical solution of an inclusion


where is a closed-graph multifunction. In contrast to existing Newton-type method is approximated on the basis of the limiting (Mordukhovich) normal cone to the graph of , computed at the respective point. Under appropriate assumptions, this method exhibits local superlinear convergence and, so far, it has been successfully implemented to the solution of a class of variational inequalities (VIs) of the first and second kind, cf. [3] and [4]. This contribution is devoted to the application of the semismooth method to the discrete 3D contact problem with Tresca friction which is modelled as a VI of the second kind. Therefore the implementation can be conducted along the lines of [4]. The paper has the following structure: In Section 2 we describe briefly the main conceptual iterative scheme of the method. Section 3 deals with the considered discrete contact problem and Section 4 concerns the suggested implementation. The results of numerical tests are then collected in Section 5.

We employ the following notation. For a cone , stands for its (negative) polar and for a multifunction , and denote its domain and its graph, respectively. The symbol “” means the convergence within the set , denotes the Frobenius norm of a matrix and signifies the ball around .

2 The semismooth Newton method

For the reader’s convenience we recall fist the definition of the tangent cone and the limiting (Mordukhovich) normal cone.

Definition 1

Let be closed and . Then

  • the cone

    is called the (Bouligand) tangent cone to at ;

  • the cone

    is called the limiting (Mordukhovich) normal cone to at .

The latter cone will be extensively used in the sequel. Let us assign to a pair two matrices such that their i-th rows, say , fulfill the condition


Moreover, let be the set of matrices satisfying (2) and

The general conceptual iterative scheme of the semismooth Newton method is stated in Algorithm 1 below.

   1 Semismooth Newton-type method for generalized equations


1:Choose a starting point , set the iteration counter .
2:If , stop the algorithm.
3:Approximation step: compute
close to such that .
4:Newton step: select and compute the new iterate
5:Set and go to 2.


Let be a (local) solution of (1). Since need not belong to or need not be close to even if is close to ; one performs in step 3 an approximate projection of onto . Therefore the step 3 is called the approximation step. The Newton step 4 is related to the following fundamental property, according to which the method has been named.

Definition 2 ([3])

Let . We say that is semismooth at provided that for every there is some such that the inequality


is valid for all and for all

If we assume that is semismooth at , then it follows from (3) that for every there is some such that for every and every pair one has

cf. [3, Proposition 4.3]. This is the background for the Newton step in Algorithm 1.

Finally, concerning the convergence, assume that is semismooth at and there are positive reals such that for every sufficiently close to the set of quadruples satisfying the conditions


is nonempty. Then it follows from [3, Theorem 4.4], that Algorithm 1 either terminates at after a finite number of steps or converges superlinearly to whenever is sufficiently close to .

The application of the semismooth Newton methods to a concrete problem of type (1) requires thus the construction of an approximation step and the Newton step which fulfill conditions (4)-(6).

3 The used model

The fundamental results concerning unilateral contact problems with Coulomb friction have been established in [6]. The infinite-dimensional model of the contact problem with Tresca friction in form of a variational inequality of the second kind can be found, e.g., in [5, 8]. Other related friction-type contact problems are described, e.g., in [7].

Figure 1: The left picture depicts an undeformed elastic prism occupying domain with the left (blue) face attached (Dirichlet condition) and some surface tractions applied to the right and top faces (depicted in green). They press the contact face against the (red) rigid plane foundation. Example of the resulting deformed body is depicted in the right picture. Front faces are not visualized.

We assume that an elastic prism occupying domain is pressed against a rigid plane foundation (cf. Figure 1). A full three-dimensional domain is discretized by a mesh of brick elements and consists of

nodes (vertices). The finite element method using trilinear basis functions is then applied to approximate a displacement field vector

in each mesh node. Entries of are ordered in such a way that and the j-th node is associated with the pair of its tangential and normal displacements, respectively.

A sparse stiffness matrix and the loading (column) vector are first assembled and then both condensed to incorporate zero displacements in Dirichlet nodes corresponding to the (blue) Dirichlet boundary. Secondly, all nodes not lying in the (bottom) contact face are eliminated by the Schur complement technique and the Cholesky factorization resulting in a dense matrix and a vector , where is the number of nodes excluding Dirichlet boundary nodes.

At last, all local blocks of and all blocks of are expanded to blocks and blocks, respectively in order to incorporate the non-penetrability condition

where is the Lagrange multiplier associated with non-penetrability constraint. Here and in the following, we assume that . In this way, we obtain a dense regular matrix and a (column) vector .

Finally, let us simplify the notation via

to define a vector of unknowns .

Following the development in [1], our model attains the form of generalized equation (GE)


where the single-valued function is given by

and the multifunction by

with being the friction coefficient. GEs of the type (7) have been studied in [4] and so all theoretical results derived there are applicable. For our approach it is also important that the Jacobian is positive definite.

4 Implementation of the semismooth method

In order to facilitate the approximation step we will solve, instead of GE (7), the enhanced system


in variables Clearly is a solution of (7), if and only if is a solution of (8).

In the approximation step we suggest to solve for all consecutively the next three low-dimensional strictly convex optimization problems:


obtaining thus their unique solutions , respectively. The solutions can be ordered in vectors and all together in a vector

Thereafter we compute the outcome of the approximation step via

Clearly and, using the theory [4, Section 4], it is possible to show that condition (4) is fulfilled.

In the Newton step we put


In (9),

is an identity matrix and block diagonal matrices

attain the form


where the diagonal blocks have the structure

and submatrices and scalar entries are computed in dependence on values and as follows:

  • If (sticking), we put , otherwise we put

  • If (no contact or weak contact), we put , otherwise we put .

This choice ensures that matrices with given by (9) fulfill conditions (5),(6) with replaced by .

Stopping rule It is possible to show (even for more general Coulomb friction model [1]) that there is a Lipschitz constant such that,

whenever the output of the approximation step lies in a sufficiently small neighborhood of . It follows that, with a sufficiently small positive , the condition


tested after the approximation step, may serve as a simple yet efficient stopping rule.

4.0.1 Computational benchmark

We assume that the domain

is described by elastic parameters (Young’s modulus), (Poisson’s ratio) and subject to surface tractions


and the friction coefficient .

level nodes assembly Cholesky nodes semismooth solver
of K (sec) Schur (sec) time (sec) iters
1 225 0.078 0.003 40 0.017 5
2 225 0.031 0.003 40 0.017 6
3 637 0.094 0.021 84 0.047 6
4 1377 0.141 0.092 144 0.101 6
5 4056 0.516 0.701 299 0.507 7
6 9537 1.297 3.968 544 1.928 7
7 27072 3.156 32.110 1104 9.734 7
8 70785 18.672 1242.211 2112 48.275 8
Table 1: Performance of the MATLAB solver.

The domain is uniformly divided into hexahedra (bricks), where

are numbers of hexahedra along with coordinate axis, denotes the mesh level of refinement and the ceiling function. Consequently the number of nodes and the number of nodes read


Figure 2: The left picture depicts the deformed contact boundary and the right figure shows the corresponding deformed elastic prism, both pictures together with the (red) rigid plane foundation.

Table 1 reports on the performance of the whole method for various meshes assuming zero initial approximation and the stopping criterion . We can clearly see that the number of iterations of the semismooth Newton method (displayed in the last column) only slightly increase with the mesh size. This behaviour shows that the method is mesh-independent.

Figure 2 visualizes a deformed contact boundary together with a deformation of the full domain obtained by post-processing. Displacements of non-contact boundary nodes are then obtained from a linear system of equations with the matrix and the vector .

All pictures and running times were produced by our MATLAB code available for download and testing at

https://www.mathworks.com/matlabcentral/fileexchange/70255 .

It is based on original codes of [1] and its performance is further enhanced by a vectorized assembly of using [2].

4.0.2 Concluding remarks and further perspectives

The choice (9) of matrices in the Newton step of the method is not unique and may be used to simplify the linear system in the Newton step. The convergence may be further accelerated by an appropriate scaling in the approximation step.

4.0.3 Acknowledgment

Authors are grateful to Petr Beremlijski (TU Ostrava) for providing original Matlab codes of [1] and discussions leading to various improvements of our implementation.