%~Mouliné par MaN_auto v.0.40.4 (bbb487d5) 2026-04-16 16:56:52
\documentclass[CRMECA,Unicode,biblatex,published]{cedram}

\addbibresource{CRMECA_Rey_20251009.bib}

\usepackage[normalem]{ulem}
\usepackage[algo2e,ruled,lined,titlenumbered,commentsnumbered]{algorithm2e}

\newcommand{\disp}{\mathrm{d}}
\newcommand{\cont}{\mathrm{c}}
\newcommand{\rrm}{\mathrm{r}}
\newcommand{\normal}{\mathrm{N}}
\newcommand{\tangent}{\mathrm{T}}
\newcommand{\nnormal}{\mathrm{N}}
\newcommand{\ntangent}{\mathrm{T}}
\newcommand{\Dirichlet}{\mathrm{d}}
\newcommand{\Neumann}{\mathrm{n}}
\newcommand{\Hrm}{\mathrm{H}}
\newcommand{\Prm}{\mathrm{P}}

\newcommand{\dif}{\mathop{}\!{\operatorfont{d}}}

\let\Arctan\arctan

\DeclareMathOperator{\sign}{sign}

\makeatletter
\renewcommand\uuline{\bgroup\markoverwith%
 {%
 \rule[-0.64ex]{2pt}{0.45pt}%
 \llap{\rule[-0.84ex]{2pt}{0.45pt}}%
}%
 \ULon}
\renewcommand\underline{\bgroup\markoverwith%
 {%
 \rule[-0.64ex]{2pt}{0.45pt}%
}%
 \ULon}
\makeatother

\newcommand{\un}[1]{\ensuremath{\smash{\underline{#1}}}}
\newcommand{\dep}{\ensuremath{\smash{\un{u}}}}
\newcommand{\uu}[1]{\ensuremath{\smash{\uuline{#1}}}}
\newcommand{\sig}{\ensuremath{\smash{\uu{\sigma}}}}
\newcommand\strain[1]{\uu{\varepsilon}(#1)}

\definecolor{rose}{rgb}{1.0, 0.33, 0.64}
\definecolor{vert}{rgb}{0.09, 0.45, 0.27}
\newcommand{\dc}{\textcolor{vert}{free}}
\newcommand{\adh}{\textcolor{rose}{stick}}
\newcommand{\slip}{\textcolor{blue}{slip}}

\renewcommand{\boldmath}[1]{#1}

%%%---------------------------------------------------------------------------

%%% DELIMITERS

\usepackage{mathtools}
\DeclarePairedDelimiter{\parens}{\lparen}{\rparen}
\DeclarePairedDelimiter{\braces}{\{}{\}}
\DeclarePairedDelimiter{\abs}{\lvert}{\rvert}
\DeclarePairedDelimiter{\norm}{\lVert}{\rVert}
\DeclarePairedDelimiter{\angles}{\langle}{\rangle}

%%%---------------------------------------------------------------------------

%%% ENVIRONEMENTS

%%% Env. {problem}
\newcounter{problem}
\renewcommand{\theproblem}{\Alph{problem}}
\newenvironment{problem}{\refstepcounter{problem}\begin{enonce*}[remark]{Problem~\theproblem}}{\end{enonce*}}

%%% Env. {system} (mimicks {cases} but only supports a single column)
\makeatletter
\newenvironment{system}{
  \let\@ifnextchar\new@ifnextchar
  \left\lbrace
  \def\arraystretch{2}
  \array {@{}l@{}}
}{\endarray\right.}
\makeatother

%%%---------------------------------------------------------------------------

%%% TABLES

\usepackage{booktabs}

\renewcommand{\arraystretch}{1.2}
\setlength{\cmidrulewidth}{.2pt}

%%%---------------------------------------------------------------------------

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\let\oldtilde\tilde
\renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}}
\let\oldhat\hat
\renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}}
\let\oldforall\forall
\renewcommand*{\forall}{\mathrel{\oldforall}}

\title{On the influence of resolution algorithm's parameters on the converged solution for a Coulomb friction contact problem exhibiting multiple solutions}
\alttitle{Influence des paramètres d’algorithmes de résolution de contact frottant sur la solution obtenue à convergence : étude sur un problème avec frottement de Coulomb à solutions multiples}

\author{\firstname{Valentine} \lastname{Rey}\CDRorcid{0000-0003-1019-1819}}
\address{Nantes Universit\'e, Ecole Centrale Nantes, CNRS, GeM, UMR 6183, F-44000 Nantes, France}
\email{valentine.rey@univ-nantes.fr}
\thanks{The author is supported by R\'egion Pays de la Loire through grant \'Etoiles montantes}

\keywords{\kwd{Quasi-static contact} \kwd{Coulomb friction} \kwd{algorithms}}
\altkeywords{\kwd{Contact frottant quasi-statique} \kwd{frottement de Coulomb} \kwd{algorithmes}}

\begin{abstract}
In this paper, we compare four algorithms solving the quasi-static contact between one elastic body and one rigid obstacle in two dimensions. With the Coulomb's law of friction, this example exhibits multiple solutions if the friction coefficient is larger than 3. After describing the numerical methods tested in this paper, we study the influence of the parameters of the algorithms on the nature of the obtained solution at convergence. On this academic example, we also compute two existing criteria on the uniqueness of the solution. Finally, for friction coefficient larger than 3, we compute a new sliding solution and observe that not all approaches are able to find it.
\end{abstract}

\begin{altabstract}
Ce papier compare quatre méthodes de résolution du problème bi-dimensionnel de contact quasi-statique entre un solide élastique et un obstacle rigide. Sur l'exemple étudié, avec une loi de frottement de Coulomb, plusieurs solutions existent lorsque le coefficient de frottement est plus grand que 3. Après avoir décrit les quatre méthodes numériques comparées, nous les mettons en œuvre sur l'exemple et étudions l'influence des paramètres de ces algorithmes sur la solution obtenue à convergence. Sur cet exemple académique, nous calculons aussi deux critères d'existence de solutions multiples. Enfin, pour un coefficient de frottement strictement supérieur à 3, nous calculons une troisième solution de type glissement et observons que toutes les approches numériques étudiées ne sont pas capables de la trouver.
\end{altabstract}

\COI{The author does not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and has declared no affiliations other than their research organizations.}

\dateposted{2026-05-04}
\begin{document}
%\input{CR-pagedemetas}
%\end{document}
\maketitle

\vskip 4\baselineskip

\section{Introduction}

In computational mechanics, the numerical simulation of structures involving frictional contact is challenging. The contact problem is non-linear and requires dedicated numerical strategies. Indeed, the non-interpenetration condition between solids introduces inequalities in the normal direction. Moreover, the Coulomb's law of friction is used to model the behaviour on the tangential directions at the contact interface~\cite{Wriggers06}. This law cannot be derived from a potential but can be derived from a bi-potential~\cite{de1998bipotential}: the Coulomb's law is not an associated friction law. As a consequence, in the weak formulation, only a variational inequality can be written. Thus, the mathematical results of existence and uniqueness are limited, which is not the case for the Tresca's law of friction. This lack or regularity may lead to the absence of solution, such as the Painlev\'e problem, or the existence of multiple solutions (possibly infinite number of solutions). Examples of such cases are available in the literature: one may cite~\cite{acary2011formulation} for dynamic frictional contact with one degree of freedom or~\cite{agwa2021existence,agwa2021existenceII,hild2005multiple,hild2004non} for multiple solutions on problems with several degrees of freedom. In the case of quasi-static frictional contact with Coulomb's law, the existence of a solution to the continuous problem has been proven for small coefficients of friction~\cite{eck2005unilateral}. The solution of the finite element discretized problem always exists but the uniqueness is proven only for small coefficients of friction~\cite{REF}. However, only few estimations of the critical friction coefficient above which multiple solutions may exist are available~\cite{agwa2021existence,andersson1999quasistatic}. Moreover, the computational cost of the estimation increases exponentially with the number of nodes on the potential contact interface. Amongst all numerical methods available to solve the quasi-static contact problem with Coulomb's law of friction, only few are dedicated to the computation of all the solutions. Some approaches are based on continuation techniques. In~\cite{ligursky2015bifurcations}, the authors propose numerical continuation based on a predictor-corrector and a bifurcation criterion. However, this criterion is based on probable branching scenarios and potential solutions may thus be missed. In~\cite{janovsky2012computing}, the authors use a parametrization of the loading path or friction coefficient and apply continuation techniques on examples with a small number of contact nodes. In~\cite{haslinger2018nonsmooth}, this method is simplified and applied to more complex examples with finer mesh. In~\cite{hassani2003mixed}, the authors propose to solve non linear eigenvalue problems to detect multiple solutions but this method fails in finding isolated solutions. Finally, considering the discretized friction contact problem as a linear complementarity problem with cone constraints, it is possible to apply existing techniques to compute all the solutions~\cite{de1992generalized,el2015linear} but the computational cost increases rapidly with the number of nodes on the contact interface. Thus, they are not suitable for large finite element computations on fine meshes.

In practice, more standard iterative numerical approaches are used. For a panorama of existing methods to solve quasi-static friction contact problems, the reader can refer to~\cite{Wriggers06,acary2018solving}. Once the algorithm has converged, the solution is used for post-process analysis. Indeed, in engineering problems, the quantities of interest such as maximum local stresses, Von~Mises or Tresca equivalent stress or stress intensity factors are computed from the nodal displacement solution. Associated to criterion and material parameters such as yield limit, those numerical results are exploited to validate the design. Yet, the quality of the numerical solution is rarely assessed. More specifically, in case of large friction coefficients, the possibility of multiple solutions is almost never studied.

In this article, we propose to study the behaviour of four iterative algorithms solving~quasi-static~frictional contact on one academic example exhibiting multiple solutions~\cite{hild2005multiple}. The first numerical method considered in this article is the formulation of contact based on augmented Lagrangian presented in~\cite{simo1992augmented}. For given contact forces (Lagrange multipliers), a non-linear problem is solved. From the obtained displacements, the contact forces are updated. These two steps are repeated until convergence. This method is popular and widely implemented in finite element codes. The second method considered in this article is the Nitsche method~\cite{chouly2017overview}. The main point of this method is that the contact forces are written in terms of displacement thanks to the constitutive relation (elasticity in this specific case). It has the advantage of introducing only one scalar parameter and to avoid the two iterative loops (non-linear problem and update of the Lagrange multipliers) existing in the first method. The third method is the one proposed in~\cite{alart1991mixed,alart1997methode}. It is also based on the optimality of augmented Lagrangian but considers both primal and dual variables at the same time. As a consequence, there is no update of the contact forces and the problem can be written only in terms of the contact forces. This method has shown to be effective for finite element computations~\cite{acary2018solving}. Finally, the fourth method considered in this article is the one proposed in~\cite{kanno2022primal}. It is based on the formulation of the mechanical problem as a linear complementarity problem with cone constraints. This problem is then solved with a primal-dual approach based on accelerated projection method coming from convex optimization community. Indeed, recently, numerical methods traditionally dedicated to optimization problems have been applied in the computational mechanics field. We can cite~\cite{ELBOUSTANI} in which second-order cone programming is applied to a frictional contact problem with an associated version of Coulomb's law (more regular and not exhibiting multiple solutions) or~\cite{acary2011formulation} for dynamic Coulomb friction problem, \cite{acary2024second} for multi-body systems made of rigid part or~\cite{zhang2013second} which requires a regularization of the Coulomb's law.

The four methods considered here lead to algorithms with one or two iterative loops. Many parameters are to be chosen by the user. The influence of these parameters on the speed of convergence is studied and well-known and some guidelines are provided to reduce the number of iterations. However, the influence of the initialization of all the iterative processes is rarely studied. The objective of this article is to study to which extent the values of the scalar parameters and the initialization values affect the nature of the solution obtained at convergence on the academic problem with multiple solutions.

In Section~\ref{sec:notations}, we present the quasi-static contact problem with Coulomb's law of friction between one elastic body and one rigid obstacle. We give the governing equations in the continuous framework and in the discretized framework. Then, in Section~\ref{sec:algos}, we briefly introduce the four methods and the associated algorithm that are compared in this paper. In Section~\ref{sec:crit_multi}, we give two criteria that enable to compute an estimation of the coefficient of friction above which multiple solutions may arise. Those two criteria come from the literature and we also comment on a third criteria proposed in~\cite{agwa2021existence}. In Section~\ref{sec:numerical_example}, after defining the academic problem, we study the influence of the algorithms' parameters on the converged solutions. For this problem, we compute the bounds on the coefficient of friction below which the solution is unique and compare them to the exact value derived in~\cite{hild2005multiple}. We consider one or multiple nodes on the potential contact boundary. We also vary the friction coefficient $\mu$ and exhibit a third solution involving sliding for $\mu>3$. Finally, Section~\ref{sec:ccl} concludes the paper.

\section{Description of the mechanical problem}\label{sec:notations}

In this section, we introduce the quasi-static Coulomb friction contact problem. We choose to consider the contact between one elastic domain and one rigid obstacle because this is the case of the numerical example studied in Section~\ref{sec:numerical_example}. Of course, the methods considered in this paper are not limited to one elastic body and are able to solve contact problems between two elastic bodies. In Section~\ref{first_subsection} we give the continuous formulation of the problem. In Section~\ref{second_subsection} we present the formulation obtained after discretization.

\subsection{Continuous problem}\label{first_subsection}

We consider the elastic body $\Omega$ in frictional contact with one rigid obstacle. The linear behaviour of $\Omega$ is characterized by the Hooke tensor $\mathbb{H}$. We note $\sig$ the Cauchy stress tensor, $\dep$ the displacement field and $\strain{\dep}$ the symmetric part of the gradient of displacements. We do the small strain and quasi-static hypotheses. On $\partial_{\Dirichlet}\Omega$, the displacement $\dep_{\Dirichlet}$ is imposed, which corresponds to Dirichlet boundary condition. On $\partial_{\Neumann}\Omega$, the pressure $\un{F}$ is imposed, which corresponds to Neumann boundary condition. We note $\Gamma$ the potential contact zone between the two solids. Finally, $\un{f}$ is the body force applied to $\Omega$. Let $\un{n}_\Gamma$ be the outward normal to $\Gamma$. We also introduce $\sigma_{\nnormal}= (\sig \un{n}_\Gamma) \un{n}_\Gamma $ and $\un{\sigma}_{\ntangent}=\sig \un{n}_\Gamma-\sigma_{\nnormal} \un{n}_\Gamma$, as well as $u_{\nnormal}=\dep\cdot \un{n}_\Gamma $ and $\un{u}_{\ntangent}=\dep-u_{\nnormal} \un{n}_\Gamma$. The gap $g$ is defined as the initial distance between $\Omega$ and the obstacle. The Hertz--Moreau--Signorini conditions in the normal direction read:
\begin{equation}\label{eq:contact_norm}
\begin{alignedat}{2}
	u_{\nnormal}-g & \leq 0
	& \quad & \text{on $\Gamma$},
\\	\sigma_{\nnormal} & \leq 0
	& \quad & \text{on $\Gamma$},
\\	\sigma_{\nnormal} (u_{\nnormal}-g) & = 0
	& \quad & \text{on $\Gamma$}.
\end{alignedat}
\end{equation}
$\mu$ is the friction coefficient and the quasi-static Coulomb law at contact distinguishes two cases: stick contact with no tangent displacement or slip contact with tangent displacement:
\begin{equation}\label{eq:contact_tgt_vitesse}
\begin{alignedat}{2}
	\norm{\un{\sigma}_{\ntangent}} \leq -\mu \sigma_{\nnormal} \quad
	& \text{and} \quad \norm{\un{\dot{u}}_{\ntangent}}=0
	& \quad & \text{on $\Gamma$ : stick},
\\	\norm{\un{\sigma}_{\ntangent}} = -\mu \sigma_{\nnormal} \frac{\un{\dot{u}}_{\ntangent}}{\norm{\un{\dot{u}}_{\ntangent}}} \quad
	& \text{and} \quad \norm{\un{\dot{u}}_{\ntangent}}\neq0
	& \quad & \text{on $\Gamma$ : slip},
\end{alignedat}
\end{equation}
where $\un{\dot{u}}_{\ntangent}$ is the tangential velocity. In the context of computational contact in dynamic, the time interval is usually discretized and the velocities are approximated by ratios of increment of displacements over time step $\un{\dot{u}}_{\ntangent}\simeq \frac{\Delta \un{u}_{\ntangent}}{\Delta t}$ which leads to the friction law written in terms of increments of displacement:
\begin{equation}\label{eq:contact_tgt_delta}
\begin{alignedat}{2}
	\norm{\un{\sigma}_{\ntangent}} \leq -\mu \sigma_{\nnormal} \quad
	& \text{and} \quad \norm{\Delta\un{u}_{\ntangent}}=0 
	& \quad & \text{on $\Gamma$ : stick},
\\	\norm{\un{\sigma}_{\ntangent}} = -\mu \sigma_{\nnormal} \frac{\Delta\un{u}_{\ntangent}}{\norm{\Delta\un{u}_{\ntangent}}} \quad
	& \text{and} \quad \norm{\Delta\un{u}_{\ntangent}}\neq0
	& \quad & \text{on $\Gamma$ : slip}.
\end{alignedat}
\end{equation}
In the quasi-static framework and considering initial zero condition for load and displacements, the friction law reads:
\begin{equation}\label{eq:contact_tgt}
\begin{alignedat}{2}
	\norm{\un{\sigma}_{\ntangent}} \leq -\mu \sigma_{\nnormal} \quad
	& \text{and} \quad \norm{\un{u}_{\ntangent}}=0
	& \quad & \text{on $\Gamma$ : stick},
\\	\norm{\un{\sigma}_{\ntangent}} = -\mu \sigma_{\nnormal} \frac{\un{u}_{\ntangent}}{\norm{\un{u}_{\ntangent}}} \quad
	& \text{and} \quad \norm{\un{u}_{\ntangent}}\neq0
	& \quad & \text{on $\Gamma$ : slip}.
\end{alignedat}
\end{equation}
In addition to the equations of contact on $\Gamma$, the boundary conditions, the constitutive relation and the equilibrium equations read:
\begin{equation}\label{eq:equilibre}
\begin{alignedat}{2}
	\dep & = \dep_{\Dirichlet}
	& \quad & \text{on $\partial_{\Dirichlet}\Omega$},
\\	\sig \un{n} & = \un{F}
	& \quad & \text{on $\partial_{\Neumann}\Omega$},
\\	\sig & = \mathbb{H} \strain{\dep}
	& \quad & \text{in $\Omega$},
\\	\un{\operatorname{div}}(\sig)+\un{f} & = \un{0}
	& \quad & \text{in $\Omega$},
\end{alignedat}
\end{equation}
where $\norm{\un{x}}$ is the Euclidean norm of vector $\un{x}$. The problem \eqref{eq:contact_norm}-\eqref{eq:contact_tgt}-\eqref{eq:equilibre} may be written as weak inequality~\cite{renard2013generalized}. In computational mechanics, continuous problems are rarely solved directly and discretization techniques are employed to obtain an approximation of the solution. In this article, we consider the Finite Element Method (FEM) as the chosen discretization technique.

\subsection{Discretization with finite elements}\label{second_subsection}

We introduce a discretization of the domain $\Omega$. We note $\Omega_{\Hrm}$ the discretized body. Linear shape functions and isoparametric elements are considered. Two functional spaces are defined:
\begin{equation}
\mathcal{V}_{\Hrm}=\braces[\big]{\un{v}\in\mathcal{H}^1(\Omega_{\Hrm}) \ \text{and} \ \un{v}= \dep_{\Dirichlet} \ \text{on $\partial_{\Dirichlet}\Omega_{\Hrm}$}}
\end{equation}
and
\begin{equation}
\mathcal{V}_{\Hrm}^0=\braces[\big]{\un{v}\in\mathcal{H}^1(\Omega_{\Hrm}) \ \text{and} \ \un{v}= \un{0} \ \text{on $\partial_{\Dirichlet}\Omega_{\Hrm}$}}.
\end{equation}

$\mathbf{K}$ is the stiffness matrix, $\mathbf{f}$ is the vector of generalized forces gathering the bodyforce and Neumann condition contributions. The imposed nodal displacement is $\mathbf{u}_{\disp}$. We note $\mathbf{r}_{\cont}$ the vector of reaction due to contact at the potential contact nodes, that is to say on $\Gamma$. Note that $\mathbf{r}_{\cont}=\mathbf{0}$ corresponds to no contact. The nodal displacement at potential contact nodes is $\mathbf{u}_{\cont}$ and $\mathbf{u}_{\rrm}$ is the nodal displacement of nodes neither in $\Gamma$ nor in $\partial_{\Dirichlet}\Omega$. We note $n_{\cont}$ the number of (potential contact) nodes in $\Gamma$. Using block decomposition of the stiffness matrix, the equilibrium reads:
\begin{equation}\label{eq:discr_equi}
{\begin{pmatrix}
\mathbf{Id} & \mathbf{0} & \mathbf{0} \\
\mathbf{0}& \mathbf{K}_{\rrm\rrm} & \mathbf{K}_{\rrm\cont} \\
\mathbf{0}& \mathbf{K}_{\cont\rrm}& \mathbf{K}_{\cont\cont}
\end{pmatrix}
\begin{pmatrix}
\mathbf{u}_{\disp} \\ \mathbf{u}_{\rrm} \\ \mathbf{u}_{\cont} \\
\end{pmatrix}=
\begin{pmatrix}
\mathbf{u}_{\disp} \\ \mathbf{f}_{\rrm}-\mathbf{K}_{\rrm\disp}\mathbf{u}_{\disp} \\ \mathbf{f}_{\cont}-\mathbf{K}_{\cont\disp}\mathbf{u}_{\disp}
\end{pmatrix}+
\begin{pmatrix}
\mathbf{0} \\\mathbf{0} \\\mathbf{r}_{\cont}
\end{pmatrix}}.
\end{equation}
We note $\mathbf{u}_{\normal}$ the normal nodal displacement for potential contact nodes and $\mathbf{u}_{\tangent}$ the tangential nodal displacement. We use the same notations to decompose the contact reaction $\mathbf{r}_{\cont}$ into $\mathbf{r}_{\normal}$ and~$\mathbf{r}_{\tangent}$. The nodal discretization of the continuous contact law with friction presented in~\eqref{eq:contact_norm} and~\eqref{eq:contact_tgt} reads
\begin{equation}\label{eq:discr_contact}
\begin{aligned}
	& \mathbf{0}\geq (\mathbf{u}_{\normal}-\mathbf{g})\perp \mathbf{r}_{\normal}\leq\mathbf{0},
\\	& \text{$\forall$ nodes $i$ such that $u_{\nnormal,i}-g_i=0$, $\abs{r_{\ntangent,i}}\leq -\mu r_{\nnormal,i}$},
\end{aligned}
\end{equation}
where $\perp$ is the orthogonality between two vectors which can be characterized by the scalar product being null. Note that this nodal discretization is not the only possibility. One can perform an integral approximation of the contact condition using some finite element spaces for the approximation of contact forces. The choice of the discretization of the contact conditions is given for the four compared methods in each dedicated section.

This first equality of~\eqref{eq:discr_contact} expresses the Hertz--Moreau--Signorini conditions in normal direction written for each node of $\Gamma$. The use of linear elements ensures that the verification of the conditions at the nodes leads to the verification of the conditions on $\Gamma$. The second equality of~\eqref{eq:discr_contact} expresses the tangent behaviour without telling the stick case from the slip case. We clearly see that the frictional contact brings inequality constraints to the classical FE linear system~\eqref{eq:discr_equi}.

This problem can also be rewritten as a mixed linear complementarity problem~\cite{agwa2021existence}, as follows.

\begin{problem}
Find $\mathbf{r}_{\normal}$, $\mathbf{r}_{\tangent}$, $\mathbf{u}_{\normal}$, $\mathbf{u}_{\tangent}$ such that
\begin{equation}\label{eq:glp}
\begin{pmatrix}
-\mathbf{Id}& \mathbf{0}& \mathbf{0}& \mathbf{F}_{\normal\normal}-\mu\mathbf{F}_{\normal\tangent}& \mathbf{F}_{\normal\tangent}& \mathbf{0} & \mathbf{g}-\mathbf{u}_{\normal}^* \\
\mathbf{0}& -\mathbf{Id}& \mathbf{0}& \mathbf{F}_{\tangent\normal}-\mu\mathbf{F}_{\tangent\tangent}& \mathbf{F}_{\tangent\tangent}& \mathbf{Id}&-\mathbf{u}_{\tangent}^* \\
\mathbf{0}& \mathbf{0}& -\mathbf{Id}& 2\mu\mathbf{Id}& -\mathbf{Id}& \mathbf{0} & \mathbf{0} \\
\end{pmatrix}
\begin{pmatrix}
	\mathbf{g}-\mathbf{u}_{\normal}
\\	\max(-\mathbf{u}_{\tangent};\mathbf{0})
\\	\mathbf{r}_{\tangent}-\mu \mathbf{r}_{\normal}
\\	-\mathbf{r}_{\normal} \\-\mathbf{r}_{\tangent}-\mu \mathbf{r}_{\normal}
\\	\max(\mathbf{u}_{\tangent};\mathbf{0})
\\	1
\end{pmatrix}=
\begin{pmatrix}\mathbf{0} \\\mathbf{0} \\\mathbf{0}
\end{pmatrix}
\end{equation}
and
\[
\begin{pmatrix}\mathbf{0} \\\mathbf{0} \\\mathbf{0}
\end{pmatrix} \leq
\begin{pmatrix}\mathbf{g}-\mathbf{u}_{\normal} \\\max(-\mathbf{u}_{\tangent};\mathbf{0})\\\mathbf{r}_{\tangent}-\mu \mathbf{r}_{\normal}
\end{pmatrix} \perp
\begin{pmatrix}-\mathbf{r}_{\normal} \\-\mathbf{r}_{\tangent}-\mu \mathbf{r}_{\normal} \\\max(\mathbf{u}_{\tangent};\mathbf{0})
\end{pmatrix}\geq
\begin{pmatrix}\mathbf{0} \\\mathbf{0} \\\mathbf{0}
\end{pmatrix} \\
\]
where $\mathbf{F}_{\normal\normal}$, $\mathbf{F}_{\normal\tangent}$, $\mathbf{F}_{\tangent\normal}$, $\mathbf{F}_{\tangent\tangent}$ are the submatrices of the Schur complement $\mathbf{F}_{\cont\cont}$ defined by
\[
\mathbf{F}_{\cont\cont}=\parens{\mathbf{K}_{\cont\cont}-\mathbf{K}_{\cont\rrm}\mathbf{K}_{\rrm\rrm}^{-1}\mathbf{K}_{\rrm\cont}}^{-1}
\]
and $\mathbf{u}_{\normal}^*$ and $\mathbf{u}_{\tangent}^*$ are the subvectors corresponding respectively to normal and tangent directions coming from the vector $\mathbf{u_c^*}$ defined by
\[
\mathbf{u_c^*}=\mathbf{F}_{\cont\cont}\parens*{\mathbf{f}_{\cont}-\mathbf{K}_{\cont\disp}\mathbf{u}_{\disp}}-\mathbf{F}_{\cont\cont}\mathbf{K}_{\cont\rrm}\mathbf{K}_{\rrm\rrm}^{-1}\parens*{\mathbf{f}_{\rrm}-\mathbf{K}_{\rrm\disp}\mathbf{u}_{\disp}}
\]
\end{problem}

Algorithms exist to compute all the solutions, see~\cite{el2015linear,de1992generalized}, but the computational cost becomes prohibitive when the number of contact nodes increases since the complexity is exponential.

\section{Presentation of considered algorithms}\label{sec:algos}

In this section, we introduce the four methods compared in Section~\ref{sec:numerical_example}. For each method, we shortly present the principles before describing the formulation and the associated algorithms. We insist on the parameters that are left to the user to define. In each subsection, we give references for readers who would obtain more details on the methods.

\subsection{Augmented Lagrangian solved with Uzawa}\label{sec:uzawa}

The augmented Lagrangian enables to formulate one minimization problem with constraints into a minimization problem with no constraint. In frictional contact mechanics, the optimality of the augmented Lagrangian is equivalent to the solution of the initial problem for Tresca modelization of friction but not for Coulomb's law. This is due to the fact that Coulomb friction cannot be derived from on potential but only from a bi-potential \cite{de1998bipotential}. However, in practice, the optimality is extended and used for Coulomb's law. The Uzawa algorithm is popular and easy to implement. It consists in defining Lagrange multipliers $\lambda_{\nnormal}$, $\un{\lambda}_{\ntangent}$ interpreted as contact reactions. Then, the following non-linear problem is solved.

\begin{problem}
Find $\dep\in\mathcal{V}_{\Hrm}$ such that $\forall \un{v}\in\mathcal{V}_{\Hrm}^0$
\begin{multline}\label{eq:NL_solve}
\int_{\Omega_{\Hrm}} \strain{\un{v}}:\mathbb{H}: \strain{\dep} \dif \Omega_{\Hrm} -\int_{\Omega_{\Hrm}} \un{f}\cdot\un{v} \dif \Omega_{\Hrm}-\int_{\partial_{\Neumann}\Omega} \un{F}\cdot\un{v} \dif s
\\	= \int_\Gamma\angles[\big]{\lambda_{\nnormal}+\eta_{\nnormal}(u_{\nnormal}-g)} {v}_{\nnormal} \dif s+\int_\Gamma\parens*{-\un{\sigma}_{\ntangent}+ \un{\lambda}_{\ntangent}+\eta_{\ntangent}\parens*{\dep_{\ntangent}-\xi\frac{\un{s}}{\norm{\un{s}}}}} \un{v}_{\ntangent} \dif s
\end{multline}
with $\langle x\rangle=\frac{x+\abs{x}}{2}$, $\un{s}=-\un{\sigma}_{\ntangent}+\un{\lambda}_{\ntangent}+\eta_{\ntangent}\dep_{\ntangent}$ and
\[
\xi=
\begin{cases}
	0 & \text{if $\norm{\un{s}}-\mu\angles[\big]{\lambda_{\nnormal}+\eta_{\nnormal}(u_{\nnormal}-g)}\leq0$},
\\	\frac{\norm{\un{s}}-\mu\angles*{\lambda_{\nnormal}+\eta_{\nnormal}(u_{\nnormal}-g)}}{\eta_{\ntangent}} & \text{otherwise}.
\end{cases}
\]
\end{problem}

This is the continuous formulation. In the FEM, the computation of the left-hand side of the equality is classical. Regarding the right-hand side of the equality, in terms of finite element discretization, we use the integral approximation of the contact condition and not the nodal approximation. In the nodal approximation, the terms on the boundary $\Gamma$ are computed as a sum of nodal values:
\begin{equation}
\int_\Gamma\angles[\big]{\lambda_{\nnormal}+\eta_{\nnormal}(u_{\nnormal}-g)} {v}_{\nnormal} \dif s= \sum_{i=1}^{n_{\cont}} \angles[\big]{\lambda_{\nnormal,i}+\eta_{\nnormal}(u_{\nnormal,i}-g_i)} {v}_{\nnormal,i}.
\end{equation}

With the integral approximation, the contact forces are discretized using shape functions. We actually choose the same shape functions as for the discretization of the displacement, which is not the optimal choice but the most common one~\cite{desmeure2016strategie}.
The non-linear problem is solved using Newton--Raphson scheme which is usually called generalized Newton method as the function is not continuously differentiable. This lack of regularity may lead to an arbitrary choice of sub-gradient in case the tangent matrix has to be computed at one iteration with the iterate being on the edge of the friction cone. In that case, the chosen sub-gradient is the adjacent value corresponding to the outside of the cone.

Once equation~\eqref{eq:NL_solve} is solved, the Lagrange multipliers are updated. Those two steps are repeated until convergence is reached. The resulting algorithm from~\cite{simo1992augmented} is Algorithm~\ref{alg:Uzawa}.

\begin{algorithm2e}[th]\caption{Augmented Lagrangian solved with Uzawa}\label{alg:Uzawa}
Initialization of $\lambda_{\nnormal}^0$, $\un{\lambda}_{\ntangent}^0$ and $k=0$\;
Find the solution $\dep^k$ of equation~\eqref{eq:NL_solve}\;
\eIf{$u^k_{\nnormal}-g\leq\rho_1$ for all nodes of\/ $\Gamma$ and $\norm{\dep_{\ntangent}^k}\leq \rho_2$ for all nodes of\/ $\Gamma$ such that $\norm{\un{\sigma}_{\ntangent}}<\mu\angles[\big]{\lambda_{\nnormal}+\eta_{\nnormal}(u_{\nnormal}-g)}$}
{convergence\;}
{$\lambda_{\nnormal}^{k+1}=\angles[\big]{\lambda_{\nnormal}^k+\eta_{\nnormal}(u_{\nnormal}^k-g)}$\;
$\un{\lambda}_{\ntangent}^{k+1}=
\begin{cases}
	\un{\lambda}_{\ntangent}^{k}+\eta_{\ntangent} \un{u}_{\ntangent}^k
	& \text{if $\norm{-\un{\sigma}_{\ntangent}+\un{\lambda}_{\ntangent}^{k}+\eta_{\ntangent} \un{u}_{\ntangent}^k}\leq\mu\lambda_{\nnormal}^{k+1}$}
\\	\frac{-\un{\sigma}_{\ntangent}+\lambda_{\ntangent}^{k}+\eta_{\ntangent} \un{u}_{\ntangent}^k}{\norm{-\un{\sigma}_{\ntangent}+\un{\lambda}_{\ntangent}^{k}+\eta_{\ntangent} \un{u}_{\ntangent}^k}}\mu\lambda_{\nnormal}^{k+1}-\un{\sigma}_{\ntangent}
	& \text{otherwise}
\end{cases}$\;
$k\leftarrow k+1$\;}
\end{algorithm2e}

The two tolerances $\rho_1$ and $\rho_2$, the value of the coefficients $\eta_{\nnormal}$ and $\eta_{\ntangent}$ and the initialization of $\lambda_{\nnormal}^0$ and $\un{\lambda}_{\ntangent}^0$, as well as the initialization and the stopping criterion of the Newton--Raphson algorithm are the parameters defined by the user. It is well-known that the penalty parameters influence the speed of convergence so the total number of iterations and thus the total computational time. Large penalty parameters speed up the convergence but may lead to ill-conditioning.

\subsection{Nitsche formulation solved with generalized Newton}

The Nitsche method was originally created to take into account Dirichlet boundary conditions. Later, it was extended to unilateral contact without or with friction, see~\cite{chouly2017overview} for a review. The main difference with the augmented Lagrangian is that no multipliers are introduced: they are replaced by stress vectors that can be obtained from primal unknowns using the constitutive behaviour. Symmetric and unsymmetric versions of the Nitsche's method exist. In this paper, the unsymmetric version was chosen as it was observed to behave better than the symmetric one~\cite{renard2013generalized}. The unsymmetric version of the formulation of the problem reads as follows.

\begin{problem}
Find $\dep\in\mathcal{V}_{\Hrm}$ such that $\forall \un{v}\in\mathcal{V}_{\Hrm}^0$
\begin{multline*}
\int_{\Omega_{\Hrm}} \strain{\un{v}}:\mathbb{H}: \strain{\dep} \dif \Omega_{\Hrm} -\int_{\Omega_{\Hrm}} \un{f}\cdot\un{v} \dif \Omega_{\Hrm}-\int_{\partial_{\Neumann}\Omega} \un{F}\cdot\un{v} \dif s
\\	= -\int_\Gamma\parens*{\sigma_{\nnormal}(\dep)-\kappa(u_{\nnormal}-g)}_-v_{\nnormal} \dif s +\int_\Gamma J\parens[\big]{\un{\sigma}_{\ntangent}(\dep)-\kappa\dep_{\ntangent}}\cdot\un{v}_{\ntangent} \dif s
\end{multline*}
with
\begin{multline}\label{eq:J}
J\parens[\big]{\un{\sigma}_{\ntangent}(\dep)-\kappa\dep_{\ntangent}}
\\	=
\begin{cases}
	\un{\sigma}_{\ntangent}(\dep)-\kappa\dep_{\ntangent}
	& \text{if $\norm[\big]{\un{\sigma}_{\ntangent}(\dep)-\kappa\dep_{\ntangent}}\leq \mu\parens*{\sigma_{\nnormal}(\dep)-\kappa(u_{\nnormal}-g)}_-$},
\\	\mu\parens*{\sigma_{\nnormal}(\dep)-\kappa(u_{\nnormal}-g)}_- \frac{\un{\sigma}_{\ntangent}(\dep)-\kappa\dep_{H,T}}{\norm{\un{\sigma}_{\ntangent}(\dep)-\kappa\dep_{\ntangent}}}
	& \text{otherwise},
\end{cases}
\end{multline}
and $(x)_-=0$ if $x>0$ and $(x)_-=-x$ if $x\leq0$.
\end{problem}

Once again, the continuous problem is discretized. The left-hand side of the equation is classical in elastic linear finite element analysis and we use the integral approximation of the contact conditions (right-hand side of the equation).

This problem is solved using a generalized Newton--Raphson algorithm whose parameters are the initialization and the stopping criterion. The function is not continuously differentiable and the associated consequences are the one explained in the previous section. The main parameter of this method is the value of the scalar $\kappa>0$ whose positiveness ensures the non interpenetration. Yet, in~\cite{renard2013generalized}, the authors explain that $\kappa$ has to be chosen sufficiently large in order to keep the coercivity of the tangent matrix.

\subsection{Alart--Curnier formulation solved with generalized Newton}

This approach was originally proposed in~\cite{alart1997methode,alart1991mixed}. It is also based on the augmented Lagrangian whose optimality is solved using a Newton--Raphson scheme. The unsymmetric version of optimality of the augmented Lagrangian writes as follows.

\begin{problem}
Find $\dep\in\mathcal{V}_{\Hrm}$ such that $\forall \un{v}\in\mathcal{V}_{\Hrm}^0$
\begin{equation}\label{eq:alart}
\begin{system}
	\displaystyle \int_{\Omega_{\Hrm}} \strain{\un{v}}:\mathbb{H}: \strain{\dep} \dif \Omega_{\Hrm} =\int_{\Omega_{\Hrm}} \un{f}\cdot\un{v} \dif \Omega_{\Hrm}+\int_{\partial_{\Neumann}\Omega} \un{F}\cdot\un{v} \dif s+\int_\Gamma \lambda_{\nnormal} v_{\nnormal} \dif s +\int_\Gamma \un{\lambda}_{\ntangent} \un{v}_{\ntangent} \dif s,
\\	\displaystyle  \frac{1}{\eta} \int_\Gamma \parens[\big]{\lambda_{\nnormal} + \parens[\big]{\lambda_{\nnormal}-\eta(u_{\nnormal}-g)}_-} w_{\nnormal} \dif s=0
	\quad \forall \un{w}\in\mathcal{W}_{\Hrm},
\\	\displaystyle \frac{1}{\eta} \int_\Gamma \parens[\big]{\un{\lambda}_{\ntangent}-\mathcal{P}(\un{\lambda}_{\ntangent}-\eta \dep_{\ntangent})} \un{w}_{\ntangent} \dif s=0
	\quad \forall \un{w}\in\mathcal{W}_{\Hrm},
\end{system}
\end{equation}
with
\begin{equation}
	\mathcal{P}(\un{\lambda}_{\ntangent}-\eta \dep_{\ntangent})=
\begin{cases}
\un{\lambda}_{\ntangent}-\eta \dep_{\ntangent}
	& \text{if $\norm{\un{\lambda}_{\ntangent}-\eta \dep_{\ntangent}} \leq \mu\parens*{\lambda_{\nnormal}-\eta(u_{\nnormal}-g)}_-$},
\\	\mu\parens*{\lambda_{\nnormal}-\eta(u_{\nnormal}-g)}_- \frac{\un{\lambda}_{\ntangent}-\eta \dep_{\ntangent}}{\norm{\un{\lambda}_{\ntangent}-\eta \dep_{\ntangent}}}
	& \text{otherwise}.
\end{cases}
\end{equation}
\end{problem}

The test function $w_{\nnormal}$ and $w_{\ntangent}$ live in the functional space $\mathcal{W}_{\Hrm}$ composed of usual finite element shape functions. The definition of $\mathcal{P}$ is very close to the one of $J$ defined in equation~\eqref{eq:J}. They both correspond to the projection on the ball of center 0 and radius $\mu$ multiplied by the normal contact force. Contrary to the Uzawa algorithm presented in Section~\ref{sec:uzawa}, the resolution does not separate the primal variables (displacements) from the dual variables (contact forces) being updated on a second step. Indeed, the system~\eqref{eq:alart} is coupled. However, the first equation is linear with respect to the displacement so that the problem can be written as a unique functional of $\lambda$. As for Nitsche and Uzawa, this functional is not continuously differentiable and the same treatment of the lack of regularity is done. The main parameter of this method is the value of the scalar $\eta>0$. One also needs to define the initialization of the Newton's scheme and its stopping criterion.

\subsection{Primal dual accelerated gradient method from convex optimization}

The last method has been proposed in~\cite{kanno2022primal}. It relies on the formulation of contact constraints as cone constraints and on an iterative algorithm coming from convex optimization community. As the Dirichlet boundary condition can be eliminated (see~\eqref{eq:discr_equi}), the problem can be written as follows.

\begin{problem}\label{problem_E}
Find $\mathbf{u}_{\rrm}$, $\mathbf{u}_{\cont}$, $\mathbf{r}_{\cont}$, such that
\begin{equation}\label{convex}
\begin{aligned}
	&	{\begin{pmatrix}
			\mathbf{K}_{\rrm\rrm} & \mathbf{K}_{\rrm\cont}
		\\	\mathbf{K}_{\cont\rrm}& \mathbf{K}_{\cont\cont}
		\end{pmatrix}
		\begin{pmatrix}
			\mathbf{u}_{\rrm} \\ \mathbf{u}_{\cont}
		\end{pmatrix} =
		\begin{pmatrix}
			\mathbf{f}_{\rrm}-\mathbf{K}_{\rrm\disp}\mathbf{u}_{\disp}
		\\	\mathbf{f}_{\cont}-\mathbf{K}_{\cont\disp}\mathbf{u}_{\disp}
		\end{pmatrix} +
		\begin{pmatrix}
			\mathbf{0}
		\\ \mathbf{r}_{\cont}
		\end{pmatrix}},
\\	& \mathcal{F}^*\ni
\begin{pmatrix}
{-g_j-\mu u_{\ntangent,j}+u_{\nnormal,j}} \\
{u_{\ntangent,j}}
\end{pmatrix} \perp
\begin{pmatrix}
{r_{\nnormal,j}} \\{r_{\ntangent,j}}
\end{pmatrix} \in \mathcal{F}
\quad \text{for $j=1, \dotsc, n_{\cont}$}.
\end{aligned}
\end{equation} 
\end{problem}

The friction cone is defined by $\mathcal{F}=\braces[\big]{({r_{\nnormal}},{r}_{\ntangent})\in\mathbb{R}\times\mathbb{R}^2, \ \text{such that $-\mu {r_{\nnormal}}\geq\norm{{r}_{\ntangent}}$}}$. The projector on~$\mathcal{F}$ is $\Pi$ and the dual cone is $\mathcal{F}^*$. Here, the nodal approximation of contact is used. This is a non-linear complementarity problem with cone constraints. Once again, because of the Coulomb model for friction, solving the optimality condition does not guarantee to obtain the solution of Problem~\ref{problem_E}. However, this is what is done in practice

The algorithm proposed in~\cite{kanno2022primal} is given in Algorithm~\ref{alg:primaldual}. It consists in working with both primal and dual variables and is based on primal-dual hybrid gradient from~\cite{chambolle2016introduction} with an acceleration scheme.
In Algorithm~\ref{alg:primaldual}, $\gamma$ is the largest singular value of the matrix of size $(d,3\times n_{\cont})$ defined by
\begin{equation}
\mathbf{T}=
\begin{pmatrix}T_{1n}& T_{2n} & \cdots &T_{n_{\cont}n} \\T_{1t}& T_{2t} & \cdots &T_{n_{\cont}t}
\end{pmatrix}
\end{equation}
where $u_{t,j}=T_{jt}^t\mathbf{u}$ and $u_{n,j}=T_{jn}^t\mathbf{u}$. The matrix $\mathbf{T}$ is composed of $0$ and $1$.
The smallest eigenvalue of matrix $
\parens[\Big]{\begin{smallmatrix}
\mathbf{K}_{\rrm\rrm} & \mathbf{K}_{\rrm\cont} \\
\mathbf{K}_{\cont\rrm}& \mathbf{K}_{\cont\cont}
\end{smallmatrix}}$ is noted $\omega$. The computation of projection $\Pi$ is easily done using spectral factorization, see~\cite{kanno2022primal}. The parameters of this method are the step length $\alpha_0$ and the initialization of $\mathbf{u}^0$ and $\mathbf{r}_{\cont}$ and the stopping criterion $\zeta$.

Note that a cheaper version of the algorithm, not requiring the resolution of one linear system at each iteration, is available~\cite{kanno2023quasi}. In this paper, we used Algorithm~\ref{alg:primaldual}.

\begin{algorithm2e}[th]\caption{Primal dual accelerated gradient algorithm}\label{alg:primaldual}
Initialization of $\mathbf{u}^0=
\parens[\Big]{\begin{smallmatrix}
\mathbf{u}_{\rrm}^0 \\ \mathbf{u}_{\cont}^0
\end{smallmatrix}}$, $\mathbf{r}_{\cont}^0$, $\alpha_0>0$\;
$\hat{\mathbf{u}}^0=\mathbf{u}^0$\; 
$\beta_0=\frac{1}{\alpha_0\gamma^2}$\;
\While{$\norm{\mathbf{u}^{k+1}-\mathbf{u}^k}>\zeta$}{
\For{$j=1, \dotsc, n_{\cont}$}{$\tilde{g}_j^k=g_j^k+\mu\norm{u_{t,j}^k}$\;
${s}_{n,j}={r}_{n,j}^k+\alpha_k(\tilde{g}_j^k-u_{t,j}^k)$\;
${s}_{t,j}={r}_{t,j}^k-\alpha_k(u_{t,j}^k)$\;
$({r}_{n,j}^{k+1},{r}_{t,j}^{k+1})=\Pi({s}_{n,j},{s}_{t,j})$\;}
Find $\mathbf{u}^{k+1}$ solution of $\parens*{\beta_k
\begin{pmatrix}
\mathbf{K}_{\rrm\rrm} & \mathbf{K}_{\rrm\cont} \\
\mathbf{K}_{\cont\rrm}& \mathbf{K}_{\cont\cont}
\end{pmatrix}+\mathbf{Id}}
\begin{pmatrix}
\mathbf{u}_{\rrm}^{k+1} \\ \mathbf{u}_{\cont}^{k+1}
\end{pmatrix}=
\begin{pmatrix}
\mathbf{u}_{\rrm}^k \\ \mathbf{u}_{\cont}^k
\end{pmatrix}+\beta_k
\begin{pmatrix}
\mathbf{f}_{\rrm}-\mathbf{K}_{\rrm\disp}\mathbf{u}_{\disp} \\ \mathbf{f}_{\cont}-\mathbf{K}_{\cont\disp}\mathbf{u}_{\disp}
\end{pmatrix}+
\begin{pmatrix}
\mathbf{0} \\ \mathbf{r}_{\cont}
\end{pmatrix}$\;
$\theta_k=\frac{1}{\sqrt{1+\omega \beta_k}}$\;
$\alpha_{k+1}=\frac{\alpha_k}{\theta_k}$\;
$\beta_{k+1}=\theta_k\beta_k$\;
$\hat{\mathbf{u}}^{k+1}=\mathbf{u}^{k+1}+\theta_k(\mathbf{u}^{k+1}-\mathbf{u}^k)$\;}
\end{algorithm2e}

\section{Two criteria for existence and multiplicity of solution}\label{sec:crit_multi}

In~\cite{agwa2021existence}, the authors review criteria to estimate the onset of multiplicity solutions. Indeed, these criteria enable to determine a critical friction coefficient. If the friction coefficient is strictly smaller than the critical friction coefficient, only one solution exists. Otherwise, multiple solutions may exist. One criterion is based on the property of P-matrices~\cite{trinkle1995dynamic} and states that the quasi-static frictional contact problem has only one solution if the following interval matrix is a P-matrix:
\begin{equation}\label{eq:crit_p}
\begin{pmatrix}
\mathbf{F}_{\normal\normal} & \mathbf{F}_{\normal\tangent} \\
\mathbf{F}_{\tangent\normal} & \mathbf{F}_{\tangent\tangent}
\end{pmatrix}
{\begin{pmatrix}
\mathbf{Id}& \mathbf{0} \\
[-\mu;\mu] & \mathbf{Id}
\end{pmatrix}}.
\end{equation}

The verification of the property for the interval matrix requires the verification of the property for $2^{n_{\cont}}$ matrices of size $(2n_{\cont};2n_{\cont})$ in two dimensions (and $(3n_{\cont};3n_{\cont})$ in three dimensions). A P-matrix is characterized by the fact that all its real eigenvalues are strictly positive and all the real eigenvalues of all its principal submatrices are also strictly positive~\cite{agwa2021existence,cottle2009linear}.

Another criterion is the one proposed in~\cite{andersson1999quasistatic}. It requires to compute a parameter $c$ as the solution of an optimization problem. Then, the critical friction coefficient is defined by
\begin{equation}
\mu_{\mathrm{crit}}=\frac{c}{\sqrt{1-c^2}}.
\end{equation}
Finally, in~\cite{agwa2021existence}, the authors propose to adapt the criterion proposed in~\cite{andersson1999quasistatic} by adding the constraint $(\mathbf{F}_{\cont\cont}\mathbf{r})^p_{\nnormal}=0$ to the optimization problem in order to reduce the computational costs. However, as this constraint is already present in the article~\cite{andersson1999quasistatic} (and Equation~58 in~\cite{agwa2021existence} is not correct), both criteria should be identical. The definition of parameter $c$ is:
\begin{equation}
c= \inf_{\substack{\mathbf{r}_{\cont} \\ \norm{\mathbf{r}_{\cont}}=1}}
\max_{\substack{p \\ \mathbf{r}^p\neq\mathbf{0} \\ (\mathbf{F}_{\cont\cont}\mathbf{r})^p_{\nnormal}=0 \\ (\mathbf{F}_{\cont\cont}\mathbf{r})^p_{\ntangent}\neq\mathbf{0}}}
\frac{r^p_{\ntangent} (\mathbf{F}_{\cont\cont}\mathbf{r})^p_{\ntangent}}{\norm{\mathbf{r}^p}\times \norm[\big]{(\mathbf{F}_{\cont\cont}\mathbf{r})^p_{\ntangent}}}
\end{equation}
where $p$ is the index for one contact node.

The computation of the critical friction coefficient, regardless the choice of the criterion, is expensive and the cost increases exponentially with the number of potential contact nodes $n_{\cont}$. In Section~\ref{sec:2criteria}, we compute these two critical coefficients of friction for the illustrative example.
\section{Application on the academic example}\label{sec:numerical_example}
\subsection{Problem with two distinct solutions}\label{sec:def_pb}

This example comes from~\cite{hild2005multiple}. Let $(\un{e}_x,\un{e}_y)$ be one orthonormal basis. We consider the triangle~$\mathit{ABC}$ with $A(0;0)$, $B(1;0)$, $C\parens*{\frac{3}{4}; \frac{1}{4}}$. This triangle may be in contact with the rigid obstacle located in $y\leq0$. The Young modulus is $E=1$~Pa and the Poisson ratio is $\nu=\frac{1}{5}$. On $[\mathit{AC}]$, the force $\un{F}_{\mathrm{Neumann}}=\frac{-35}{48}\sqrt{10}\delta\un{e}_y$ is imposed with $\delta\in {\bigl]}0;\frac{4}{3}{\bigr[}$. On $[\mathit{BC}]$, the imposed displacement is $\un{u}_{D}=6\delta(x-1)\un{e}_x+\frac{3}{4}\delta (x-1)\un{e}_y$. There is no bodyforce. The edge $[\mathit{AB}]=\Gamma$ is the potential contact zone. If $\mu\geq3$, then the following displacement fields are solutions:
\begin{itemize}
\item $\un{u}_1=\delta(7x+y-7) \un{e}_x+\delta\parens[\big]{1-x-\frac{7}{4}y}\un{e}_y$, which corresponds to no contact;
\item $\un{u}_2=\delta(-6y) \un{e}_x+\delta\parens[\big]{-\frac{3}{4}y}\un{e}_y$, which corresponds to stick contact.
\end{itemize}
If $\mu<3$, the no contact displacement is the only solution. We highlight that the two solutions exist with $\delta$ arbitrarily small. Moreover, the two displacement solutions are linear with respect to $x$ and $y$ which means that the discretization with one linear triangle is sufficient to describe these solutions. As a consequence, the mesh does not introduce discretization error on these two displacement fields. In the computations, we set $\delta=1$.

To begin with, we consider the mixed linear complementarity problem formulation~\eqref{eq:glp} and use the algorithm proposed in~\cite{el2015linear} to compute all the solutions. We do this for different discretizations ($n_{\cont}=1, \dotsc, 7$) and different friction coefficients (between $0$ and $5$ with $0.1$ increments). For $\mu\in \mathopen{]}0;3\mathclose{[}$, in agreement with the derivations in~\cite{hild2005multiple}, the only solution obtained is the no-contact solution (\dc). For $\mu=3$ exactly, as expected, we always obtain only two solutions: stick (\adh) and no contact (\dc). For $\mu\in \mathopen{]}3;5]$, we obtain three solutions: the no-contact solution, the stick solution and also a third solution (\slip) for which some nodes are in slip contact, represented in Figure~\ref{fig:slip}.

\begin{figure}[h!]
\includegraphics[scale=0.42]{slip3_5}
\caption{Third solution with partial slip and stick contact obtained for $\mu=3.5$ and $n_{\cont}=25$.}\label{fig:slip}
\end{figure}

This solution was not mentioned in~\cite{hild2005multiple}. We observe that this solution is not linear with respect to $x$ and $y$ so that the discretization chosen has an influence of the estimation of this solution.

\subsection{Estimation of the critical friction coefficient}\label{sec:2criteria}

We compute the estimate of the two critical coefficients of friction, for several meshes, as described in Section~\ref{sec:crit_multi}. We note $\hat{\mu}_{\mathrm{P}}$ the value obtained with the criterion based on the P-matrices property~\eqref{eq:crit_p}. We note $\hat{\mu}_c$ the value obtained after solving the optimization problem~\eqref{eq:crit_c}. To compute $\hat{\mu}_{\mathrm{P}}$, we start from $\mu=0$ and increase its value by 0.01 increments until the property is not verified, which leads to the possible loss of uniqueness. To solve the optimization problem~\eqref{eq:crit_c}, we use JuMP (domain-specific modeling language for mathematical optimization embedded in Julia~\cite{Lubin2023}) with Gurobi optimizer~\cite{gurobi} except for the case of one contact node ($n_{\cont}=1$) for which the optimization problem can be solved analytically. Indeed, for $n_{\cont}=1$ in the problem defined in Section~\ref{sec:def_pb} with $\delta=1$, we can write:
\begin{equation}
\mathbf{F}_{\cont\cont}=
{\begin{pmatrix}
6.6&-3\\-3&6.6
\end{pmatrix}}.
\end{equation}

Equation~58 in~\cite{agwa2021existence} reads:
\begin{equation}
c = \inf_{\substack{\mathbf{r}_{\cont} \\ \norm{\mathbf{r}_{\cont}}=1}}
\max_{\substack{p \\ \mathbf{r}^p\neq\mathbf{0} \\ (\mathbf{F}_{\cont\cont}\mathbf{r})^p_t\neq\mathbf{0}}}
\frac{r^p_t (\mathbf{F}_{\cont\cont}\mathbf{r})^p_t}{\norm{\mathbf{r}^p}\times \norm[\big]{(\mathbf{F}_{\cont\cont}\mathbf{r})^p_t}}.
\end{equation}
By writing $\mathbf{r}_{\cont}=\parens[\big]{\cos(\theta) \sin(\theta)}^T=\parens{r_n^p r_t^p}^T$ which is a unit vector, noticing that since $n_{\cont}=1$, then $\mathbf{r}^p=\mathbf{r}_{\cont}$, then we have:
\begin{equation}
c = \inf_{\substack{\theta \in[0;2\pi] \\ -3\cos(\theta) +6.6\sin(\theta)\neq0}} \sin(\theta)\sign\parens[\big]{-3\cos(\theta) +6.6\sin(\theta)}
\end{equation}
which leads to
\begin{equation}
c=\inf_{\substack{\theta \in[0;2\pi] \\ \theta\neq\theta_0}} \sin(\theta)\sign\parens[\big]{-3\cos(\theta) +6.6\sin(\theta)}
\end{equation}
with $\theta_0=\Arctan\parens*{\frac{1}{2.2}}$. We obtain analytically $c\simeq-0.4694<0$. It shows that the constraint $(\mathbf{F}_{\cont\cont}\mathbf{r})^p_n=0$ is necessary. If we apply this constraint, we solve:
\begin{equation}
c = \inf_{\substack{\mathbf{r}_{\cont} \\ \norm{\mathbf{r}_{\cont}}=1}}
\max_{\substack{p \\ \mathbf{r}^p\neq\mathbf{0} \\ (\mathbf{F}_{\cont\cont}\mathbf{r})^p_{\nnormal}=0 \\ (\mathbf{F}_{\cont\cont}\mathbf{r})^p_{\ntangent}\neq\mathbf{0}}} \frac{r^p_{\ntangent} (\mathbf{F}_{\cont\cont}\mathbf{r})^p_{\ntangent}}{\norm{\mathbf{r}^p} \times \norm[\big]{(\mathbf{F}_{\cont\cont}\mathbf{r})^p_{\ntangent}}}
\end{equation}
which writes
\begin{equation}\label{eq:crit_c}
c=\inf_{\substack{\theta \in[0;2\pi] \\ \theta\neq\theta_0 \\ \theta=\theta*}} \sin(\theta)\sign\parens[\big]{-3\cos(\theta) +6.6\sin(\theta)}
\end{equation}
with $\theta* = \arctan(2.2)$. We finally obtain $c\simeq0.9103>0$ and then $\hat{\mu}_c=2.2$ which is in agreement with the other critical friction coefficient.

\begin{table}[ht]
\caption{Critical friction coefficients obtained with two criteria and several mesh sizes.}\label{tab:mu_crit}
\small
\begin{tabular}{@{}c@{\hspace{2\tabcolsep}}ccccccc@{}}
\toprule
& \multicolumn{7}{c}{\boldmath{$n_{\cont}$}} \\
\cmidrule{2-8}
& 1 &2 &3&4 &5 &6 &7 \\
\midrule
\boldmath{$\hat{\mu}_{\mathrm{P}}$} & 2.20 & 0.65& 1.41& 1.79 & 1.63 &1.69& 1.78 \\
\boldmath{$\hat{\mu}_c$} & 2.20 & 2.10 & 2.12 & 2.18 & 2.15 &2.16& 2.16 \\
\bottomrule
\end{tabular}
\end{table}

The results for different mesh sizes are gathered in Table~\ref{tab:mu_crit}. We observe that the critical friction coefficient obtained with the P-matrices property strongly depends on the discretization. Moreover, the value obtained is not close to the exact value which is $\mu=3$. Finally, larger mesh sizes were not considered since the computational time was prohibitive.

\subsection{Comparison of the three algorithms for \texorpdfstring{$\mu=3$}{mu=3}}

In this section, we study the influence of the parameters of each algorithm on the result obtained as convergence. Since the two identified solutions are linear, we first start with only one finite element in the domain, that is to say $n_{\cont}=1$. It enables to do a parametric study of the initialization. Then, we consider a finer mesh and focus on the scalar parameters of each method.

\subsubsection{With one triangular element}

We set the stopping criterion of all the Newton--Raphson algorithms to $10^{-10}$ and we also set the stagnation stopping criterion to $\zeta=10^{-10}$ in the primal-dual approach. Moreover, in the Uzawa algorithm, the verification of the contact inequalities is done with the accuracy $\rho_1=\rho_2=10^{-10}$. The converged solution is compared with the two reference solutions (non contact $\un{u}_1$ and stick contact $\un{u}_2$) given in Section~\ref{sec:def_pb}. More precisely, the energetic norm of the error between the converged solution and the exact solutions are computed. In the numerical results, if this error norm between the converged solution and $\un{u}_1$ is smaller than $10^{-10}$, then the converged solution obtained is denoted by \dc. If this error norm between the converged solution and $\un{u}_2$ is smaller than $10^{-10}$, then the converged solution obtained is denoted by \adh.

We first verified that all methods converge at the first iteration to the initial solution when the initial solution is one of the exact solution (\adh\ or \dc), whatever the parameters defined by the user. The only exception is the Uzawa method with initialization with stick contact solution that converges to stick contact for large values of $\eta$.

Then, we study the influence of the initialization by choosing different values. Note that \textcolor{rose}{$\mathbf{u}^0=(0;0)$} (and/or \textcolor{rose}{$\boldsymbol{\lambda}=\mathbf{r}^0= (0.42;1.25)$}) corresponds to the initialization with the stick solution and that \textcolor{vert}{$\mathbf{u}^0=(-7\delta, \delta)$} (and/or \textcolor{vert}{$\boldsymbol{\lambda}=\mathbf{r}^0= (0.;0)$}) corresponds to the initialization with the no-contact solution. We also consider initializations that may not respect the contact inequality (traction for dual variable or interpenetration) in order to test the robustness of the algorithms. The number of total iterations is not limited. Finally, for the Uzawa algorithm, the penalty parameters vary from $10^4$ to $10^{13}$. For the Nitsche method, the scalar parameter $\kappa$ varies from $1$ to $10^7$. For the Alart--Curnier method, the parameter $\eta$ varies from $10^{-7}$ to $10^7$. Finally, the step length $\alpha_0$ in the primal-dual method varies from $10^{-5}$ to $2$.

The results for Uzawa algorithm are available in Tables~\ref{tab:uzawa_u} and~\ref{tab:uzawa_lam}. They correspond to various initializations of displacements in Table~\ref{tab:uzawa_u} and of nodal multipliers in Table~\ref{tab:uzawa_lam}. We first observe that the algorithm always converges to one of the two exact solutions. The stick solution is obtained only for large values of penalty parameters $\eta_{\nnormal}$ and $\eta_{\ntangent}$ if the initial multiplier is $0$. It is not the case for nonzero initializations of the Lagrange multipliers $\lambda$, see Table~\ref{tab:uzawa_lam}.

\begin{table}[ht]
\caption{Uzawa, $\mu=3.0$: type of the converged solution with various initializations $\mathbf{u}^0$ and values of $\eta=\eta_{\nnormal}=\eta_{\ntangent}$ and zero initial Lagrange multipliers.}\label{tab:uzawa_u}
\renewcommand{\arraystretch}{1.3}
{\scriptsize
\begin{tabular}{@{}l@{\hspace{3\tabcolsep}}ccccccccc@{}}
\toprule
& \multicolumn{9}{c}{$\mathbf{u}^0$} \\
\cmidrule{2-10}
\boldmath{$\eta$} & $(0.5;0)$ & $(0.5;0.5)$ & $(0;0.5)$ & $(0;0)$ & $(-7\delta, \delta)$ & $(-0.5;0)$ & $(-0.5;-0.5)$ & $(0;-0.5)$ & $(-1.84\delta;0)$\\
\midrule
$10^4$		& \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc \\
$10^5$		& \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc \\
$10^6$		& \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc \\
$10^7$		& \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc \\
$10^9$		& \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc \\
$10^{11}$	& \dc & \dc & \adh & \adh & \dc & \dc & \adh & \adh & \dc \\
$10^{13}$	& \dc & \dc & \adh & \adh & \dc & \dc & \adh & \adh & \dc \\
\bottomrule
\end{tabular}}
\end{table}

\begin{table}[ht]
\caption{Uzawa, $\mu=3.0$: type of the converged solution with various initializations $\boldsymbol{\lambda}^0=(\lambda_{\nnormal},\Delta\lambda_{\ntangent})$ and values of $\eta=\eta_{\nnormal}=\eta_{\ntangent}$ and zero as initial displacement ($\mathbf{u}^0=(0;0)$).}\label{tab:uzawa_lam}
\renewcommand{\arraystretch}{1.3}
{\scriptsize
\begin{tabular}{@{}l@{\hspace{3\tabcolsep}}cccccccc@{}}
\toprule
& \multicolumn{8}{c}{$\boldsymbol{\lambda}^0$} \\
\cmidrule{2-9}
\boldmath{$\eta$} & $(0.5;0.)$ & $(0.5;0.5)$ & $(0.5;-0.5)$ & $(0.;0.)$ & $(0.;0.5)$ & $(0.; -0.5)$ & $(0.42;1.25)$ & $(0.26\mu;0.26)$\\
\midrule
$10^4$		& \adh & \adh & \adh & \dc & \dc & \dc & \dc & \dc \\
$10^5$		& \adh & \adh & \adh & \dc & \dc & \dc & \dc & \dc \\
$10^6$		& \adh & \adh & \adh & \dc & \dc & \dc & \dc & \dc \\
$10^7$		& \adh & \adh & \adh & \dc & \dc & \dc & \dc & \dc \\
$10^9$		& \adh & \adh & \adh & \dc & \dc & \dc & \dc & \dc \\
$10^{11}$	& \adh & \adh & \adh & \adh & \dc & \dc & \adh & \adh \\
$10^{13}$	& \adh & \adh & \adh & \adh & \dc & \dc & \adh & \adh \\
\bottomrule
\end{tabular}}
\end{table}

The results for the Nitsche method are available in Table~\ref{tab:nitsche}. We first observe that the algorithm always converges to one of the two exact solutions. When interpenetration is forced at the initialization, the obtained solution is stick. Whatever the value of $\kappa$, the initialization with the exact \dc\ or \adh\ solution converges to the same exact solution. Note that for some initializations, the nature of the obtained solution depends on the parameter $\kappa$.

\begin{table}[h]
\caption{Nitsche, $\mu=3.0$: type of the converged solution with various initializations $\mathbf{u}^0$ and values of $\kappa$.}\label{tab:nitsche}
\renewcommand{\arraystretch}{1.3}
{\scriptsize
\begin{tabular}{@{}l@{\hspace{3\tabcolsep}}ccccccccc@{}}
\toprule
& \multicolumn{9}{c}{$\mathbf{u}^0$} \\
\cmidrule{2-10}
\boldmath{$\kappa$} & $(0.5;0)$ & $(0.5;0.5)$ & $(0;0.5)$ & $(0;0)$ & $(-7\delta, \delta)$ & $(-0.5;0)$ & $(-0.5;-0.5)$ & $(0;-0.5)$ & $(-1.84\delta;0)$\\
\midrule
1		& \adh & \adh & \adh & \adh & \dc & \adh & \adh & \adh & \adh \\
2		& \adh & \adh & \adh & \adh & \dc & \adh & \adh & \adh & \adh \\
5		& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \adh \\
10		& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \adh \\
50		& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \adh \\
100		& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \adh \\
500		& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \adh \\
$10^3$	& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \adh \\
$10^5$	& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \adh \\
$10^6$	& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \adh \\
$10^7$	& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \adh \\
\bottomrule
\end{tabular}}
\end{table}

\pagebreak

The results for the Alart--Curnier generalized Newton method are available in Table~\ref{tab:alart}. Once again, the algorithm always converges to one of the two exact solutions. The choice of the parameter $\kappa$ does have an influence on the solution except if the initialization corresponds to one of the exact solutions ($(0;0)$ for no contact or $(0.42;1.25)$ for stick contact).

\begin{table}[ht]
\caption{Alart--Curnier, $\mu=3.0$: type of the converged solution with various initializations $\boldsymbol{\lambda}^0$ and values of $\eta$.}\label{tab:alart}
\renewcommand{\arraystretch}{1.3}
{\scriptsize
\begin{tabular}{@{}l@{\hspace{3\tabcolsep}}cccccccc@{}}
\toprule
& \multicolumn{8}{c}{$\boldsymbol{\lambda}^0$} \\
\cmidrule{2-9}
\boldmath{$\eta$} & $(0.5;0.)$ & $(0.5;0.5)$ & $(0.5;-0.5)$ & $(0.;0.)$ & $(0.;0.5)$ & $(0.; -0.5)$ & $(0.42;1.25)$ & $(0.26\mu;0.26)$\\
\midrule
$10^{-7}$	& \dc & \adh & \dc & \dc & \adh & \dc & \adh & \adh \\
$10^{-5}$	& \dc & \adh & \dc & \dc & \adh & \dc & \adh & \adh \\
$10^{-3}$	& \dc & \adh & \dc & \dc & \adh & \dc & \adh & \adh \\
0.1			& \dc & \adh & \dc & \dc & \adh & \dc & \adh & \adh \\
1			& \dc & \dc & \adh & \dc & \dc & \dc & \adh & \adh \\
10			& \dc & \dc & \adh & \dc & \dc & \adh & \adh & \adh \\
$10^3$		& \dc & \dc & \adh & \dc & \dc & \adh & \adh & \adh \\
$10^5$		& \dc & \dc & \adh & \dc & \dc & \adh & \adh & \adh \\
$10^7$		& \dc & \dc & \adh & \dc & \dc & \adh & \adh & \adh \\
\bottomrule
\end{tabular}}
\end{table}

The results for the primal-dual algorithm are available in Tables~\ref{tab:pd_u} and~\ref{tab:pd_r}. We first observe that the algorithm always converges to one of the two exact solutions. Most of the time, the no-contact solution is obtained.

\begin{table}[!h]
\caption{Primal-dual, $\mu=3.0$: type of the converged solution with various initializations $\mathbf{u}^0$ and values of $\alpha_0$ for $\mathbf{r}^0=(0;0)$.}\label{tab:pd_u}
\renewcommand{\arraystretch}{1.3}
{\scriptsize
\begin{tabular}{@{}l@{\hspace{3\tabcolsep}}ccccccccc@{}}
\toprule
& \multicolumn{9}{c}{$\mathbf{u}^0$} \\
\cmidrule{2-10}
\boldmath{$\alpha_0$} & $(0.5;0)$ & $(0.5;0.5)$ & $(0;0.5)$ & $(0;0)$ & $(-7\delta, \delta)$ & $(-0.5;0)$ & $(-0.5;-0.5)$ & $(0;-0.5)$ & $(-1.84\delta;0)$ \\
\midrule
$10^{-5}$	& \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc \\
$10^{-3}$	& \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc \\
0.01			& \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc \\
0.05			& \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc \\
0.1			& \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc \\
0.5			& \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc \\
1			& \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc \\
2			& \adh & \dc & \dc & \dc & \dc & \dc & \adh & \adh & \dc \\
\bottomrule
\end{tabular}}
\end{table}

\begin{table}[!h]
\caption{Primal-dual, $\mu=3.0$: type of the converged solution with various initializations $\mathbf{r}^0$ and values of $\alpha_0$ for $\mathbf{u}^0=(0;0)$.}\label{tab:pd_r}
\renewcommand{\arraystretch}{1.3}
{\scriptsize
\begin{tabular}{@{}l@{\hspace{3\tabcolsep}}cccccccc@{}}
\toprule
& \multicolumn{8}{c}{$\mathbf{r}^0$} \\
\cmidrule{2-9}
$\alpha_0$ & $(0.5;0.)$ & $(0.5;0.5)$ & $(0.5;-0.5)$ & $(0.;0.)$ & $(0.;0.5)$ & $(0.; -0.5)$ & $(0.42;1.25)$ & $(0.26;0.26\mu)$ \\
\midrule
$10^{-5}$	& \dc & \dc & \dc & \dc & \dc & \dc & \adh & \dc \\
$10^{-3}$	& \dc & \dc & \dc & \dc & \dc & \dc & \adh & \dc \\
0.01			& \dc & \dc & \dc & \dc & \dc & \dc & \adh & \dc \\
0.05			& \dc & \dc & \dc & \dc & \dc & \dc & \adh & \dc \\
0.1			& \dc & \dc & \dc & \dc & \dc & \dc & \adh & \dc \\
0.5			& \dc & \dc & \dc & \dc & \dc & \dc & \adh & \dc \\
1			& \dc & \dc & \dc & \dc & \dc & \dc & \adh & \dc \\
2			& \dc & \dc & \dc & \dc& \dc & \dc & \adh & \dc \\
\bottomrule
\end{tabular}}
\end{table}

\subsubsection{With several elements}\label{sec:sv_3}

In this section, we increase the number of degrees of freedom such that $n_{\cont}=50$. A parameter study of the initialization is no longer possible. Thus, the initial fields are randomly generated. Then, 100 computations are performed and the nature of the converged solution is studied. We highlight that 100 computations may not be sufficient to obtain statistically converged results. The objective of this section is to perform a first coarse analysis on the probability of obtaining contact or no-contact solution.

For Uzawa algorithm, the initial multipliers are such that $\lambda_{\nnormal}$ is generated using a uniform law on the interval $[0;0.5]$. Then, $\lambda_{\ntangent}$ and the initial displacement are generated using a uniform law on the interval $[-0.5;0.5]$. The converged solution is always the no-contact solution, as illustrated in Table~\ref{tab:uzawa_multidof}. We consider $\eta_{\nnormal}=\eta_{\ntangent}=\eta\leq100$ otherwise the algorithm diverges. For the primal-dual accelerated gradient approach, the initial displacements and reactions at contact nodes are generated using a uniform law on $[-0.5;0.5]$. Only the no-contact solution is found by the algorithm, as illustrated in Table~\ref{tab:pd_multidof}. This may be explained by the fact that vanishing the constraints is easier than enforcing them while verifying equilibrium.

\begin{table}[ht]
\caption{Uzawa, $\mu=3.0$: type of the converged solution for 100 random initializations and different values of $\eta_{\nnormal}=\eta_{\ntangent}$.}\label{tab:uzawa_multidof}
\small
\begin{tabular}{@{}c@{\hspace{3\tabcolsep}}cccc@{}}
\toprule
& \multicolumn{3}{c}{$\eta_{\nnormal}=\eta_{\ntangent}$} \\
\cmidrule{2-4}
\textbf{Type of solution} & 10&50&100 \\
\midrule
\dc		& 100& 100&100 \\
\adh		&0 & 0& 0 \\
\bottomrule
\end{tabular}
\end{table}

\begin{table}[ht]
\caption{Primal dual, $\mu=3.0$: type of the converged solution for 100 random initializations and different values of $\alpha_0$.}\label{tab:pd_multidof}
\small
\begin{tabular}{@{}c@{\hspace{3\tabcolsep}}ccccccc@{}}
\toprule
& \multicolumn{7}{c}{$\alpha_0$} \\
\cmidrule{2-8}
\textbf{Type of solution} & $10^{-3}$&0.01&0.05&0.1&0.5&1&2 \\
\midrule
\dc		& 100& 100&100 & 100&100 & 100&100 \\
\adh		& 0 & 0& 0& 0&0 & 0& 0 \\
\bottomrule
\end{tabular}
\end{table}

For the Nitsche method, only the initial displacement is required and it follows a uniform law on $[-0.5;0.5]$. The results obtained are given in Table~\ref{tab:nitsche_multidof}. We observe that with small values of~$\kappa$, only the no-contact solution is found. When $\kappa$ increases, the stick solution is more likely to happen.

\begin{table}[ht]
\caption{Nitsche, $\mu=3.0$: type of the converged solution for 100 random initializations and different values of~$\kappa$.}\label{tab:nitsche_multidof}
\small
\begin{tabular}{@{}c@{\hspace{3\tabcolsep}}ccccccc@{}}
\toprule
& \multicolumn{7}{c}{$\kappa$} \\
\cmidrule{2-8}
\textbf{Type of solution} & 50&100&500& $10^3$ & $10^5$ & $10^6$ & $10^7$ \\
\midrule
\dc		& 100 &100 &100 &100 &19 & 0& 24 \\
\adh		& 0 & 0& 0&0 & 81&100 & 76 \\
\bottomrule
\end{tabular}
\end{table}

For the Alart--Curnier method, only the initial contact forces are required and it follows a uniform law on $[-0.5;0.5]$. The results obtained are given in Table~\ref{tab:alart_multidof} and show that only the no-contact solution is found.

\begin{table}[ht]
\caption{Alart Curnier, $\mu=3.0$: type of the converged solution for 100 random initializations and different values of $\eta$.}\label{tab:alart_multidof}
\small
\begin{tabular}{@{}c@{\hspace{3\tabcolsep}}cccccc@{}}
\toprule
& \multicolumn{6}{c}{$\eta$} \\
\cmidrule{2-7}
\textbf{Type of solution} & $10^{-5}$&$10^{-3}$&0.1&1&10& $10^3$ \\
\midrule
\dc		& 100 &100 &100 &100 &100 & 100 \\
\adh		& 0 & 0& 0&0 & 0&0 \\
\bottomrule
\end{tabular}
\end{table}

\subsection{Considering \texorpdfstring{$\mu=3.5$}{mu=3.5}}

In this section, we perform the same analysis as previously by setting $\mu=3.5$ in order to see if the algorithms are able to capture the sliding solution computed in Section~\ref{sec:def_pb}.

\subsubsection{With only one element}

All methods converge to the initial solution when the initial solution is one of the exact solution both dual and primal variables. As for the case $\mu=3.0$, it is true for the Uzawa method if the penalty parameters $\eta$ are large enough for the stick contact case. Note that initialization \textcolor{blue}{$\mathbf{u}^0=(-1.84\delta;0)$} or \textcolor{blue}{$\mathbf{r}^0=(0.26\mu;0.26)$} corresponds to the slip solution for the discretization $n_{\cont}=1$.

The results for Uzawa algorithm are available in Tables~\ref{tab:uzawa_u_35} and~\ref{tab:uzawa_lam_35}. We first observe that the algorithm always converges to stick or non contact but never to the slip solutions. Contrary to what was obtained for $\mu=3.0$, the stick solution is found more often for different initial displacements.

\begin{table}[!h]
\caption{Uzawa, $\mu=3.5$: type of the converged solution with various initializations $\mathbf{u}^0$ and values of $\eta=\eta_{\nnormal}=\eta_{\ntangent}$ and zero initial Lagrange multipliers.}\label{tab:uzawa_u_35}
\renewcommand{\arraystretch}{1.3}
{\scriptsize
\begin{tabular}{@{}l@{\hspace{3\tabcolsep}}ccccccccc@{}}
\toprule
& \multicolumn{9}{c}{$\mathbf{u}^0$} \\
\cmidrule{2-10}
\boldmath{$\eta$} & $(0.5;0)$ & $(0.5;0.5)$ & $(0;0.5)$ & $(0;0)$ & $(-7\delta, \delta)$ & $(-0.5;0)$ & $(-0.5;-0.5)$ & $(0;-0.5)$ & $(-1.84\delta;0)$ \\
\midrule
$10^4$		& \dc & \dc & \adh & \adh & \dc & \dc & \adh & \adh & \dc \\
$10^5$		& \dc & \dc & \adh & \adh & \dc & \dc & \adh & \adh & \dc \\
$10^6$		& \dc & \dc & \adh & \adh & \dc & \dc & \adh & \adh & \dc \\
$10^7$		& \dc & \dc & \adh & \adh & \dc & \dc & \adh & \adh & \dc \\
$10^9$		& \dc & \dc & \adh & \adh & \dc & \dc & \adh & \adh & \dc \\
$10^{11}$	& \dc & \dc & \adh & \adh & \dc & \dc & \adh & \adh & \dc \\
$10^{13}$	& \dc & \dc & \adh & \adh & \dc & \dc & \adh & \adh & \dc \\
\bottomrule
\end{tabular}}
\end{table}

\begin{table}[ht]
\caption{Uzawa, $\mu=3.5$: type of the converged solution with various initializations $\boldsymbol{\lambda}^0=(\lambda_{\nnormal},\Delta\lambda_{\ntangent})$ and values of $\eta=\eta_{\nnormal}=\eta_{\ntangent}$ and zero as initial displacement ($\mathbf{u}^0=(0;0)$).}\label{tab:uzawa_lam_35}
\renewcommand{\arraystretch}{1.3}
{\scriptsize
\begin{tabular}{@{}l@{\hspace{3\tabcolsep}}cccccccc@{}}
\toprule
& \multicolumn{8}{c}{$\boldsymbol{\lambda}^0$} \\
\cmidrule{2-9}
\boldmath{$\eta$} & $(0.5;0.)$ & $(0.5;0.5)$ & $(0.5;-0.5)$ & $(0.;0.)$ & $(0.;0.5)$ & $(0.; -0.5)$ & $(0.42;1.25)$ & $(0.26\mu;0.26)$ \\
\midrule
$10^4$		& \adh & \adh & \adh & \adh & \dc & \dc & \dc & \dc \\
$10^5$		& \adh & \adh & \adh & \adh & \dc & \dc & \dc & \dc \\
$10^6$		& \adh & \adh & \adh & \adh & \dc & \dc & \dc & \dc \\
$10^7$		& \adh & \adh & \adh & \adh & \dc & \dc & \dc & \dc \\
$10^9$		& \adh & \adh & \adh & \adh & \dc & \dc & \dc & \dc \\
$10^{11}$	& \adh & \adh & \adh & \adh & \dc & \dc & \adh & \adh \\
$10^{13}$	& \adh & \adh & \adh & \adh & \dc & \dc & \adh & \adh \\
\bottomrule
\end{tabular}}
\end{table}

The results for Nitsche method are available in Table~\ref{tab:nitsche_35}. We first observe that the algorithm always converges to one of the three exact solutions. In the last column, we give the energetic error norm between the converged displacement and the exact slip solution. The slip solution is obtained only if the initialization corresponds to the slip solution and for large values of $\kappa$. The stick solution is found more often than previously.

\begin{table}[t]
\caption{Nitsche, $\mu=3.5$: type of the converged solution with various initializations $\mathbf{u}^0$ and values of $\kappa$.}\label{tab:nitsche_35}
\renewcommand{\arraystretch}{1.3}
{\scriptsize
\begin{tabular}{@{}l@{\hspace{3\tabcolsep}}ccccccccc@{}}
\toprule
& \multicolumn{9}{c}{$\mathbf{u}^0$} \\
\cmidrule{2-10}
\boldmath{$\kappa$} & $(0.5;0)$ & $(0.5;0.5)$ & $(0;0.5)$ & $(0;0)$ & $(-7\delta, \delta)$ & $(-0.5;0)$ & $(-0.5;-0.5)$ & $(0;-0.5)$ & $(-1.84\delta;0)$ \\
\midrule
1 		& \adh & \adh & \adh & \adh & \dc & \adh & \adh & \adh & \adh \\
2		& \adh & \adh & \adh & \adh & \dc & \adh & \adh & \adh & \adh \\
5		& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \adh \\
10		& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \textcolor{blue}{3.8 $10^{-3}$} \\
50		& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \textcolor{blue}{1.4 $10^{-4}$} \\
100		& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \textcolor{blue}{3.6 $10^{-5}$} \\
500		& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \textcolor{blue}{1.4 $10^{-6}$} \\
$10^3$	& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \textcolor{blue}{3.6 $10^{-7}$} \\
$10^5$	& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \slip \\
$10^6$	& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \slip \\
$10^7$	& \adh & \dc & \dc & \adh & \dc & \adh & \adh & \adh & \slip \\
\bottomrule
\end{tabular}}
\end{table}

\begin{table}[t]
\caption{Alart--Curnier, $\mu=3.5$: type of the converged solution with various initializations $\boldsymbol{\lambda}^0$ and values of $r$.}\label{tab:alart_35}
\renewcommand{\arraystretch}{1.3}
{\scriptsize
\begin{tabular}{@{}l@{\hspace{3\tabcolsep}}cccccccc@{}}
\toprule
& \multicolumn{8}{c}{$\boldsymbol{\lambda}^0$} \\
\cmidrule{2-9}
\boldmath{$r$} & $(0.5;0.)$ & $(0.5;0.5)$ & $(0.5;-0.5)$ & $(0.;0.)$ & $(0.;0.5)$ & $(0.; -0.5)$ & $(0.42;1.25)$ & $(0.26\mu;0.26)$ \\
\midrule
$10^{-7}$	& \dc & \adh & \dc & \dc & \adh & \dc & \adh & \slip \\
$10^{-5}$	& \dc & \adh & \dc & \dc & \adh & \dc & \adh & \slip \\
$10^{-3}$	& \dc & \adh & \dc & \dc & \adh & \dc & \adh & \slip \\
0.1			& \dc & \adh & \dc & \dc & \adh & \dc & \adh & \slip \\
1			& \dc & \dc & \slip & \dc & \dc & \dc & \adh & \slip \\
10			& \dc & \dc & \slip & \dc & \dc & \slip & \adh & \slip \\
$10^3$		& \dc & \dc & \slip & \dc & \dc & \slip & \adh & \slip \\
$10^5$		& \dc & \dc & \slip & \dc & \dc & \slip & \adh & \slip \\
$10^7$		& \dc & \dc & \slip & \dc & \dc & \slip & \adh & \slip \\
\bottomrule
\end{tabular}}
\end{table}

The results for the Alart--Curnier generalized Newton method are available in Table~\ref{tab:alart_35}. Once again, if the initialization corresponds to one of the exact solutions, the algorithm converges in 0 iteration and $\kappa$ has no influence. Otherwise, both the initialization and $\kappa$ have an influence on the converged solution.

The results for the primal-dual algorithm are available in Tables~\ref{tab:pd_u_35} and~\ref{tab:pd_r_35}. We first observe that the algorithm always converges to stick or non contact but never to the slip solutions. Contrary to what was obtained for $\mu=3.0$, the stick solution is found more often for different initial displacements and reactions.

\begin{table}[ht]
\caption{Primal-dual, $\mu=3.5$: type of the converged solution with various initializations $\mathbf{u}^0$ and values of $\alpha_0$ for $\mathbf{r}^0=(0;0)$.}\label{tab:pd_u_35}
\renewcommand{\arraystretch}{1.3}
{\scriptsize
\begin{tabular}{@{}l@{\hspace{3\tabcolsep}}ccccccccc@{}}
\toprule
& \multicolumn{9}{c}{$\mathbf{u}^0$} \\
\cmidrule{2-10}
\boldmath{$\alpha_0$} & $(0.5;0)$ & $(0.5;0.5)$ & $(0;0.5)$ & $(0;0)$ & $(-7\delta, \delta)$ & $(-0.5;0)$ & $(-0.5;-0.5)$ & $(0;-0.5)$ & $(-1.84\delta;0)$ \\
\midrule
$10^{-5}$	& \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc \\
$10^{-3}$	& \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc \\
0.01			& \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc \\
0.05			& \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc \\
0.1			& \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc & \dc \\
0.5			& \adh & \adh & \adh & \adh & \dc & \adh & \adh & \adh & \dc \\
1			& \adh & \adh & \adh & \adh & \dc & \adh & \adh & \adh & \dc \\
2			& \adh & \adh & \adh & \adh & \dc & \adh & \adh & \adh & \dc \\
\bottomrule
\end{tabular}}
\end{table}

\begin{table}[ht]
\caption{Primal-dual, $\mu=3.5$: type of the converged solution with various initializations $\mathbf{r}^0$ and values of $\alpha_0$ for $\mathbf{u}^0=(0;0)$.}\label{tab:pd_r_35}
\renewcommand{\arraystretch}{1.3}
{\scriptsize
\begin{tabular}{@{}l@{\hspace{3\tabcolsep}}cccccccc@{}}
\toprule
& \multicolumn{8}{c}{$\mathbf{r}^0$} \\
\cmidrule{2-9}
\boldmath{$\alpha_0$} & $(0.5;0.)$ & $(0.5;0.5)$ & $(0.5;-0.5)$ & $(0.;0.)$ & $(0.;0.5)$ & $(0.; -0.5)$ & $(0.42;1.25)$ & $(0.26\mu;0.26)$ \\
\midrule
$10^{-5}$	& \dc & \dc & \dc & \dc & \dc & \dc & \adh & \adh \\
$10^{-3}$	& \dc & \dc & \dc & \dc & \dc & \dc & \adh & \adh \\
0.01			& \dc & \dc & \dc & \dc & \dc & \dc & \adh & \adh \\
0.05			& \dc & \dc & \dc & \dc & \dc & \dc & \adh & \adh \\
0.1			& \dc & \dc & \dc & \dc & \dc & \dc & \adh & \adh \\
0.5			& \adh & \adh & \adh & \adh& \adh & \adh & \adh & \adh \\
1			& \adh & \adh & \adh & \adh& \adh & \adh & \adh & \adh \\
2			& \adh & \adh & \adh & \adh& \adh & \adh & \adh & \adh \\
\bottomrule
\end{tabular}}
\end{table}

\subsubsection{With several elements}

In this section, we increase the number of degrees of freedom such that $n_{\cont}=50$. We generate the initial fields as it is described in Section~\ref{sec:sv_3}. With the Uzawa algorithm, as for $\mu=3.0$, the only solution obtained is no contact. For the Alart--Curnier approach, the only solution obtained is no contact. The results obtained with the Nitsche method are in Table~\ref{tab:nitsche_multidof_35}. The slip solution is obtained only once and the parameter $\kappa$ has to be large enough. In this table, we also indicate that some computations fail. Indeed, the initialization may be very wrong (with many points with interpenetration in the normal direction for instance). These failed computations correspond to increments being larger at each iteration leading to the divergence of the algorithm.

\begin{table}[ht]
\caption{Nitsche, $\mu=3.5$: type of the converged solution for 100 random initializations and different values of~$\kappa$.}\label{tab:nitsche_multidof_35}
\small
\begin{tabular}{@{}c@{\hspace{3\tabcolsep}}ccccccccc@{}}
\toprule
& \multicolumn{7}{c}{$\kappa$} \\
\cmidrule{2-8}
\textbf{Type of solution} & 50&100&500& $10^3$ & $10^5$ & $10^6$ & $10^7$ \\
\midrule
\dc		& 97 &94 &98 &97 &100& 100& 4 \\
\adh		& 0 & 0& 0&0 & 0&0 & 95 \\
\slip	& 0&0&0&0&0&0&1 \\
fail		&3&6&2&3&0&0&0\\
\bottomrule
\end{tabular}
\end{table}

\pagebreak

For the primal-dual accelerated gradient method, the initial displacements and reactions at contact nodes are generated using a uniform law on $[-0.5;0.5]$. The results are given in Table~\ref{tab:pd_multidof_35} and illustrate that the nature of the converged solution only depends on the scalar parameter~$\alpha_0$ and not on the initialization (or weakly and thus a sample of 100 computations is not large enough to be representative in terms of statistics).

\begin{table}[ht]
\caption{Primal dual, $\mu=3.5$: type of the converged solution for 100 random initializations and different values of $\alpha_0$.}\label{tab:pd_multidof_35}
\small
\begin{tabular}{@{}c@{\hspace{3\tabcolsep}}ccccccc@{}}
\toprule
& \multicolumn{7}{c}{$\alpha_0$} \\
\cmidrule{2-8}
\textbf{Type of solution} & $10^{-3}$&0.01&0.05&0.1&0.5&1&2 \\
\midrule
\dc		& 100& 0&0 & 0&0 & 0&0 \\
\adh		& 0 & 100& 100& 100&100 & 100& 100 \\
\bottomrule
\end{tabular}
\end{table}

\section{Conclusion}\label{sec:ccl}

In this paper, we compared four contact algorithms on one specific contact with Coulomb friction example: this quasi-static contact problem between one elastic body and one rigid obstacle exhibits multiple solutions if $\mu\geq3$. We studied the influence of the initialization and the scalar parameters of each approach. We observed that both initialization and parameters have an influence on the nature of the solution. When $\mu=3$, it seems easier for the algorithms to go towards the no-contact solution. In addition, we computed two existing criteria on the presence of multiple solutions for this problem. We conclude that these criteria are very expensive and not informative enough. Indeed, only a lower bound is provided and this bound may strongly depend on the discretization. Finally, considering friction coefficient $\mu>3$ we found a third solution with some sliding nodes. It appears that only the Nitsche and Alart--Curnier approaches are able to obtain this solution which is yet very unlikely. This numerical study highlights the fact that giving guidelines to users regarding the definition of the algorithms' parameters is crucial. More generally, tools to estimate the quality of the approximated solution or the presence of multiple solutions are not sufficient and yet critically needed.

\printCOI

\printbibliography

\end{document}