In  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.  and . 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 . 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.
Let be closed and . Then
is called the (Bouligand) tangent cone to at ;
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
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 ()
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
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
3 The used model
The fundamental results concerning unilateral contact problems with Coulomb friction have been established in . 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 .
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 vectorin 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 , 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  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 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
In the Newton step we put
is an identity matrix and block diagonal matricesattain 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 .
Stopping rule It is possible to show (even for more general Coulomb friction model ) 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 .
|of K (sec)||Schur (sec)||time (sec)||iters|
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
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
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.
Authors are grateful to Petr Beremlijski (TU Ostrava) for providing original Matlab codes of  and discussions leading to various improvements of our implementation.
-  P. Beremlijski, J. Haslinger, M. Kočvara, R. Kučera, J.V. Outrata, Shape optimization in 3D contact problems with Coulomb friction, SIAM J. Optimization 20(2009), 416–444.
-  M. Čermák, S. Sysala, J. Valdman, Efficient and flexible MATLAB implementation of 2D and 3D elastoplastic problems, Applied Mathematics and Computation 355 (2019), 595-614.
-  H. Gfrerer, J.V. Outrata, On a semismooth* Newton method for solving generalized equations, accepted in SIOPT, arXiv:1904.09167.
-  H. Gfrerer, J.V. Outrata, J. Valdman, On the application of the semismooth* Newton method to variational inequalities of the second kind, arXiv:1904.09167.
-  J. Haslinger, I. Hlaváček, J. Nečas, Numerical methods for unilateral problems in solid mechanics, in P. G. Ciarlet and J. L. Lions, editors, Handbook of numerical analysis, pages 313–485. Elsevier Science, 1996.
-  J. Nečas, J. Jarušek, J. Haslinger, On the solution of the variational inequality to the Signorini problem with small friction, Boll. Unione Mat. Ital. V. Ser. B17 (1980), 796–811.
-  P. Neittaanmäki, S. Repin, J. Valdman, Estimates of deviations from exact solutions of elasticity problems with nonlinear boundary conditions, Russ. J. Numer. Anal. Math. Model. 28 (2013), No.6, 597–630.
-  J. V. Outrata, M. Kočvara, J. Zowe, Nonsmooth Approach to Optimization Problems with Equilibrium Constraints, Kluwer Academic Publishers, Dordrecht, 1998.