# A Well-Posed Semi-Discretization of Elastodynamic Contact Problems with Friction

Tomas Ligursky, Yves Renard

► To cite this version:

Tomas Ligursky, Yves Renard. A Well-Posed Semi-Discretization of Elastodynamic Contact Problems with Friction. Quarterly Journal of Mechanics and Applied Mathematics, 2011, 64 (2), pp.215-238. 10.1093/qjmam/hbr004. hal-01577746

HAL Id: hal-01577746

https://hal.science/hal-01577746

Submitted on 28 Aug 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A WELL-POSED SEMI-DISCRETIZATION OF ELASTODYNAMIC CONTACT PROBLEMS WITH FRICTION

by T. LIGURSKÝ†

(Department of Numerical Mathematics, Charles University in Prague, Sokolovská 83, 18675 Prague 8, Czech Republic)

and

Y. RENARD‡

(Université de Lyon, CNRS, INSA-Lyon, ICJ UMR5208, F 69621, Villeurbanne, France)

It has been shown numerically in previous works that the well-posedness of the spatial semi-discretization plays a crucial role in obtaining stable numerical schemes for elastodynamic frictionless contact problems. The purpose of this paper is thus to introduce a mass redistribution method adapted to elastodynamic contact with Coulomb friction that guarantees the well-posedness of the semi-discrete problem. It is shown that a differentiated treatment has to be applied to the friction condition. Some numerical tests illustrating the gain in stability for the midpoint time integration scheme are presented. They suggest also that, although the differentiated treatment is necessary for the well-posedness, it is not always mandatory from the numerical viewpoint.

## 1. Introduction

The aim of this paper is to describe a spatial well-posed semi-discretization for elastodynamic unilateral contact problems with Coulomb friction. This kind of problems is of interest in computational mechanics, where situations with a frictional contact condition are fairly common. Moreover, it is well known that the full discretization of elastodynamic contact problems induces a number of difficulties.

Numerous works have already been dedicated to the construction of numerical schemes that are as much as possible stable, respecting the contact constraint and not leading to spurious oscillations. Among the strategies already proposed in the literature, we refer to (1) for a time integration scheme that is adapted to take into account a restitution coefficient coming from an impact law. Although this approach is better fitting to the case of rigid solids, the addition of an impact law makes the semi-discrete problem also well posed. However, the nature of this problem is a measure differential inclusion in time (see (1 to 4)), which is a very low regular problem. Another proposed strategy is to build energy dissipative schemes. This is the case in (5, 6), where the contact force is implicit. The drawback of this method is that the kinetic energy of the contacting nodes is cancelled at each
impact. A less radical solution is to build energy conserving schemes. Such schemes are introduced in (7 to 9). However, energy conserving schemes either introduce spurious oscillations on the contact boundary or allow a small interpenetration. It is possible to build energy conserving schemes with a penalized contact condition (9 to 11), but this also leads to important oscillations of the normal stress. In this context, it was detected early that a key point is the satisfaction of the complementarity condition between the sliding velocity and the contact pressure, the so-called persistency condition (8, 10, 12). But a compromise has to be made between the satisfaction of this condition and the non-penetration condition.

A common point to all these works is that they are focused on finding a good time integration scheme. However, in (13) and (14), it is shown that this is rather obtaining a well-posed and regular spatial semi-discrete problem that allows for stable schemes (see also (15) and (16) for further developments). The spatial semi-discretizations proposed in (13) and (14) allow the use of any reasonable time integration scheme while almost all time integration schemes are unstable with the standard discretization. However, these works are focused on contact conditions without friction. One might think that the strategies developed there are directly applicable to friction condition. We will see thereafter that this does not provide the well-posedness results and therefore a strategy adapted to the friction condition is needed. This is also reflected in studies presented in (17) and (18), where it was shown that adding a mass on the contact boundary regularizes the tangential friction problem and prevents the occurrence of multiple solutions in elastodynamics. Note that the semi-discrete problem obtained by finite elements naturally adds a mass on the nodes of the contact boundary (but this is also the main difficulty for the unilateral contact condition!).

The method proposed in this paper is to apply the redistribution mass method introduced in (13) only on the unilateral contact condition not on the friction one. We show that in this case, the semi-discrete problem in space reduces to a differential inclusion with a unique Lipschitz continuous solution (not to a measure differential inclusion as in the standard semi-discretization).

For the sake of simplicity, we limit ourselves to the small deformations framework. However, the same kind of difficulty exists for large deformation problems and similar strategies can be applied.

The outline of the paper is the following. In Section 2, we present a classical semi-discretization of an elastodynamic contact problem with friction. In Section 3, we propose an adaptation of the mass redistribution method, namely to apply it only on the normal component. Then, the well-posedness of the obtained semi-discrete problem is proved in Section 4. The unique solution of the semi-discrete problem is proved to be energy decreasing in Section 5. An elementary example is described in Section 6. It shows that the well-posedness of the fully discrete problem cannot be ensured when the mass redistribution method is applied to both contact and friction conditions. Finally, Section 7 is devoted to a numerical test that confirms the advantage of the mass redistribution method for the stability of the midpoint scheme.

## 2. A classical finite-element approximation

In this section, we introduce a classical spatial semi-discretization based on the finite-element method. Since we are mainly interested in the semi-discrete problem, we do not describe the weak formulation of the continuous problem. More details about such a discretization can be found for instance in (19 to 21).

Let $\Omega \subset \mathbb{R}^2$ be a bounded Lipschitz domain representing the reference configuration of a linearly elastic body. The Neumann condition is prescribed on $\Gamma_N$, the Dirichlet one on $\Gamma_D$ and a unilateral contact with the Coulomb friction law with respect to a rigid foundation on $\Gamma_C$ (see Fig. 1).
<figure><img src="image_2.png" /><figcaption>Fig. 1 Elastic body in contact with a rigid foundation</figcaption></figure>

We suppose that $\Gamma_N$, $\Gamma_D$ and $\Gamma_C$ form a partition of $\partial\Omega$, the boundary of $\Omega$. Let also $\rho$, $\sigma(u)$, $\varepsilon(u)$ and $A$ be the mass density, the stress tensor, the linearized strain tensor and the elasticity tensor, respectively. The elastodynamic problem consists in finding the displacement field $u(t, x)$ satisfying

$$ egin{cases} 
\rho\ddot{u} - \operatorname{div}\sigma(u) = f & \text{in } (0, T] \times \Omega, \\
\sigma(u) = A\varepsilon(u) & \text{in } (0, T] \times \Omega, \\
u = 0 & \text{on } (0, T] \times \Gamma_D, \\
\sigma(u)v = g & \text{on } (0, T] \times \Gamma_N, \\
u(0) = u^0, \dot{u}(0) = v^0 & \text{in } \Omega,
\end{cases} 
\quad (1) $$

where, additionally, $T > 0$ determines the time interval of interest, $v \in \mathbb{R}^2$ is the outward unit normal vector to $\Omega$ on $\partial\Omega$ and $f, g$ are some given external loads. Assuming the $C^1$ regularity for $\Gamma_C$, we decompose the displacement and the stress vector into normal and tangential components on $\Gamma_C$ as follows:

$$ egin{aligned} u_N &= u \cdot v, & u_T &= u \cdot \tau, \\ \sigma_N(u) &= (\sigma(u)v) \cdot v, & \sigma_T(u) &= (\sigma(u)v) \cdot \tau, \end{aligned} $$

where $\tau \in \mathbb{R}^2$ is a tangent unit vector orthogonal to $v$. Without real loss of generality, we also assume that there is no initial gap between the solid and the rigid foundation. Denoting by $\mathcal{F}$ the friction coefficient, the unilateral contact condition with Coulomb friction is expressed as follows:

$$ egin{cases} u_N \le 0, \sigma_N(u) \le 0, u_N\sigma_N(u) = 0 \\ |\sigma_T(u)| \le -\mathcal{F}\sigma_N(u) \\ \sigma_T(u) = \mathcal{F}\sigma_N(u) \frac{\dot{u}_T}{|\dot{u}_T|} & \text{if } \dot{u}_T \ne 0 \end{cases} 
\quad (2) $$

Now, we consider a vector Lagrange finite-element method defined on $\Omega$. Let $a_1, \dots, a_n$ be the finite-element nodes and $\varphi_1, \dots, \varphi_{n_p}$ the (vector) shape functions of the finite-element displacement space. We denote by $n_c$ the number of nodes on $\Gamma_C$ and by $n_p$ the number of degrees of
freedom. Let $\mathbf{u}(t)$ be the vector of degrees of freedom of the finite-element displacement field $u^h(t, x)$ such that

$$u^h(t, x) = \sum_{1 \le i \le n_p} u_i(t) \varphi_i(x) \quad \text{and} \quad \mathbf{u}(t) = (u_i(t)) \in \mathbb{R}^{n_p}.$$

Let $a_{\alpha_i}, i = 1, \dots, n_c$, be the $i$-th contact node and $\mathbf{v}_i, \boldsymbol{\tau}_i \in \mathbb{R}^{n_p}, i = 1, \dots, n_c$, be the vectors linking a displacement vector with its normal and tangential displacements at $a_{\alpha_i}$, that is,

$$\left\{ \begin{array}{ll} u_{\mathrm{N}}^{h}(t, a_{\alpha_i}) = \mathbf{v}_{i}^{T} \mathbf{u}(t), & i = 1, \dots, n_c, \\ u_{\mathrm{T}}^{h}(t, a_{\alpha_i}) = \boldsymbol{\tau}_{i}^{T} \mathbf{u}(t), & \end{array} \right. \qquad (3)$$

and satisfying

$$\left\{ \begin{array}{l} \|\mathbf{v}_i\| = 1, \quad \|\boldsymbol{\tau}_i\| = 1, \quad \mathbf{v}_i^T \boldsymbol{\tau}_j = 0, \quad \forall i, j = 1, \dots, n_c, \\ \mathbf{v}_i^T \mathbf{v}_j = 0, \quad \boldsymbol{\tau}_i^T \boldsymbol{\tau}_j = 0, \quad \forall i, j = 1, \dots, n_c, \quad i \neq j. \end{array} \right. \qquad (4)$$

We denote by $\|\cdot\|$ the Euclidean norm of vectors in $\mathbb{R}^n$ as well as the matrix norm in $\mathcal{M}_{n,m}(\mathbb{R})$ generated by the Euclidean vector norm. Using a nodal approximation of the contact condition, the spatial semi-discretization of Problem (1) and (2) can be written as follows:

$$\left\{ \begin{array}{l} \text{Find } \mathbf{u}: [0, T] \to \mathbb{R}^{n_p}, \lambda_{\nu}, \lambda_{\tau}: [0, T] \to \mathbb{R}^{n_c} \text{ such that} \\ \lambda_{\nu}(t) \in \Lambda_{\nu}, \lambda_{\tau}(t) \in \Lambda_{\tau}(\mathcal{F}\lambda_{\nu}(t)) \quad \text{a.e. in } (0, T), \\ M\ddot{\mathbf{u}}(t) + A\mathbf{u}(t) = \mathbf{f} + \mathbf{B}_{\nu}^{T}\boldsymbol{\lambda}_{\nu}(t) + \mathbf{B}_{\tau}^{T}\boldsymbol{\lambda}_{\tau}(t) \quad \text{a.e. in } (0, T), \\ (\boldsymbol{\mu}_{\nu} - \boldsymbol{\lambda}_{\nu}(t))^{T}\mathbf{B}_{\nu}\mathbf{u}(t) \ge 0, \quad \forall \boldsymbol{\mu}_{\nu} \in \Lambda_{\nu} \quad \text{a.e. in } (0, T), \\ (\boldsymbol{\mu}_{\tau} - \boldsymbol{\lambda}_{\tau}(t))^{T}\mathbf{B}_{\tau}\dot{\mathbf{u}}(t) \ge 0, \quad \forall \boldsymbol{\mu}_{\tau} \in \Lambda_{\tau}(\mathcal{F}\boldsymbol{\lambda}_{\nu}(t)) \quad \text{a.e. in } (0, T), \\ \mathbf{u}(0) = \mathbf{u}^{0}, \quad \dot{\mathbf{u}}(0) = \mathbf{v}^{0}, \end{array} \right. \qquad (5)$$

where

$$A_{ij} = \int_{\Omega} \mathcal{A}\varepsilon(\varphi_i): \varepsilon(\varphi_j) \, dx \quad \text{and} \quad M_{ij} = \int_{\Omega} \rho \varphi_i \cdot \varphi_j \, dx \quad (1 \le i, j \le n_p)$$

are the components of the stiffness matrix $\mathbf{A} \in \mathcal{M}_{n_p}(\mathbb{R})$ and of the mass matrix $\mathbf{M} \in \mathcal{M}_{n_p}(\mathbb{R})$, respectively. We assume that the tensor $\mathcal{A}$ of elasticity coefficients obey the usual symmetry and uniform ellipticity conditions, the density $\rho$ is bounded from below by a positive constant and $\Gamma_D$ is of non-zero measure on $\partial\Omega$. As a consequence, $\mathbf{A}$ and $\mathbf{M}$ are both symmetric positive definite matrices. The components of the load vector $\mathbf{f} \in \mathbb{R}^{n_p}$ are given by

$$f_i = \int_{\Omega} f \cdot \varphi_i \, dx + \int_{\Gamma_C} g \cdot \varphi_i \, d\Gamma.$$

We assume for simplicity that the load vector is time independent. Finally, $\boldsymbol{\lambda}_{\nu} = (\lambda_{\nu,1}, \dots, \lambda_{\nu,n_c})^T$ and $\boldsymbol{\lambda}_{\tau} = (\lambda_{\tau,1}, \dots, \lambda_{\tau,n_c})^T$ are the normal and tangential Lagrange multipliers, respectively, $\mathbf{B}_{\nu} = (\mathbf{v}_1, \dots, \mathbf{v}_{n_c})^T$, $\mathbf{B}_{\tau} = (\boldsymbol{\tau}_1, \dots, \boldsymbol{\tau}_{n_c})^T \in \mathcal{M}_{n_c,n_p}(\mathbb{R})$ and

$$\mathbf{\Lambda}_{\nu} = \mathbb{R}^{n_c},$$

$$\mathbf{\Lambda}_{\tau}(\mathcal{F}\mathbf{\mu}_{\nu}) = \{\boldsymbol{\mu}_{\tau} \in \mathbb{R}^{n_c}: |\boldsymbol{\mu}_{\tau,i}| \le -\mathcal{F}\boldsymbol{\mu}_{\nu,i}, \forall i = 1, \dots, n_c\}, \quad \mathbf{\mu}_{\nu} \in \mathbf{\Lambda}_{\nu},$$

stand for the Lagrange multiplier sets, $\mathbb{R}^{n_c}$ being the cone of all non-positive vectors in $\mathbb{R}^{n_c}$.
Problem (5) can be viewed as a measure differential inclusion (see (1, 3)). It is ill-posed unless an impact law is added on each contact node. Even in this case, the solutions have a very low regularity.

### 3. The mass redistribution method

The analysis presented in (13) highlights the fact that the main cause of ill-posedness is due to the inertia of the nodes on the contact boundary. It is proposed a method that consists in the redistribution of the mass near the contact boundary. This technique allows to recover the well-posedness of the semi-discrete problem and ensures the solution to be energy conserving. Moreover, it transforms the measure differential inclusion corresponding to (5) into a regular Lipschitz continuous ordinary differential equation, which can be approximated by any reasonable time integration scheme.

The singular dynamic method introduced in (14) for unilateral conditions is similar and more general than the mass redistribution method since, for instance, it can be applied to thin structures. However, we use here the mass redistribution method. The reason is that we need a differentiated treatment of unilateral and friction conditions, which would be more difficult to obtain with the singular dynamic method. In Section 6, an elementary example illustrates the fact that an undifferentiated treatment leads to a potential multiplicity of solutions.

Let $\mathcal{N} := \text{span}\{\boldsymbol{v}_1, \dots, \boldsymbol{v}_{n_c}\}$ and $\mathcal{N}^{\perp}$ denote the space spanned by $\boldsymbol{v}_i$ and its orthogonal complement, respectively. We shall consider the redistributed mass matrix $\mathbf{M}_r \in \mathcal{M}_{n_p}(\mathbb{R})$ satisfying:

$$ egin{cases} 
(i) \quad \mathbf{M}_r = \mathbf{M}_r^T, \\ 
(ii) \quad \text{Ker}(\mathbf{M}_r) = \mathcal{N}, \\ 
(iii) \quad \boldsymbol{w}^T \mathbf{M}_r \boldsymbol{w} > 0, \quad \forall \boldsymbol{w} \in \mathcal{N}^{\perp}, \quad \boldsymbol{w} \neq \mathbf{0}, 
\end{cases} 
\quad (6) $$

that is, being symmetric positive semi-definite with the kernel equal to $\mathcal{N}$. In (13), a simple algorithm is proposed to build the redistributed mass matrix preserving the main characteristics of the mass matrix (total mass, center of gravity and moments of inertia).

Using the decomposition $\mathbf{u}(t) = \mathbf{u}_{\mathcal{N}^{\perp}}(t) + \mathbf{u}_{\mathcal{N}}(t)$, $\mathbf{u}_{\mathcal{N}^{\perp}}(t) \in \mathcal{N}^{\perp}$, $\mathbf{u}_{\mathcal{N}}(t) \in \mathcal{N}$, of the displacement vector for any time $t$ and replacing $\mathbf{M}$ with $\mathbf{M}_r$, Problem (5) becomes

$$ egin{cases} 
\text{Find } \mathbf{u}_{\mathcal{N}^{\perp}}: [0, T] \to \mathcal{N}^{\perp}, \quad \mathbf{u}_{\mathcal{N}}: [0, T] \to \mathcal{N}, \quad \boldsymbol{\lambda}_{\nu}, \boldsymbol{\lambda}_{\tau}: [0, T] \to \mathbb{R}^{n_c} \text{ such that} \\ 
\qquad \boldsymbol{\lambda}_{\nu}(t) \in \mathbf{\Lambda}_{\nu}, \quad \boldsymbol{\lambda}_{\tau}(t) \in \mathbf{\Lambda}_{\tau}(\mathcal{F}\boldsymbol{\lambda}_{\nu}(t)) \quad \text{a.e. in } (0, T), \\ 
\qquad \mathbf{M}_r \ddot{\boldsymbol{u}}_{\mathcal{N}^{\perp}}(t) + \mathbf{A}(\mathbf{u}_{\mathcal{N}^{\perp}}(t) + \mathbf{u}_{\mathcal{N}}(t)) = \mathbf{f} + \mathbf{B}_{\nu}^{T} \boldsymbol{\lambda}_{\nu}(t) + \mathbf{B}_{\tau}^{T} \boldsymbol{\lambda}_{\tau}(t) \quad \text{a.e. in } (0, T), \\ 
\qquad (\boldsymbol{\mu}_{\nu} - \boldsymbol{\lambda}_{\nu}(t))^{\mathrm{T}} \mathbf{B}_{\nu} \mathbf{u}_{\mathcal{N}}(t) \ge 0, \quad \forall \boldsymbol{\mu}_{\nu} \in \mathbf{\Lambda}_{\nu} \quad \text{a.e. in } (0, T), \\ 
\qquad (\boldsymbol{\mu}_{\tau} - \boldsymbol{\lambda}_{\tau}(t))^{\mathrm{T}} \mathbf{B}_{\tau} \dot{\boldsymbol{u}}_{\mathcal{N}^{\perp}}(t) \ge 0, \quad \forall \boldsymbol{\mu}_{\tau} \in \mathbf{\Lambda}_{\tau}(\mathcal{F}\boldsymbol{\lambda}_{\nu}(t)) \quad \text{a.e. in } (0, T), \\ 
\qquad \mathbf{u}_{\mathcal{N}^{\perp}}(0) = \mathbf{u}_{\mathcal{N}^{\perp}}^{0}, \quad \dot{\mathbf{u}}_{\mathcal{N}^{\perp}}(0) = \mathbf{v}_{\mathcal{N}^{\perp}}^{0}, 
\end{cases} 
\quad (7) $$

where $\mathbf{u}_{\mathcal{N}^{\perp}}^0$ and $\mathbf{v}_{\mathcal{N}^{\perp}}^0$ are the projections of the initial values of the displacement and velocity vectors into $\mathcal{N}^{\perp}$, respectively. Moreover, it is possible to express the unilateral contact and friction
conditions in an equivalent way (see (20) for instance) and rewrite the problem as follows:

$$
\left\{
\begin{array}{l}
\text{Find } \mathbf{u}_{\mathcal{N}^{\perp}}: [0, T] \to \mathcal{N}^{\perp}, \mathbf{u}_{\mathcal{N}}: [0, T] \to \mathcal{N}, \lambda_{\nu}, \lambda_{\tau}: [0, T] \to \mathbb{R}^{n_c} \text{ such that} \\
\mathbf{M}_{\tau} \ddot{\mathbf{u}}_{\mathcal{N}^{\perp}}(t) + \mathbf{A}(\mathbf{u}_{\mathcal{N}^{\perp}}(t) + \mathbf{u}_{\mathcal{N}}(t)) = \mathbf{f} + \sum_{i=1}^{n_c} \lambda_{\nu,i}(t) \mathbf{v}_i + \sum_{i=1}^{n_c} \lambda_{\tau,i}(t) \mathbf{\tau}_i \quad \text{a.e. in } (0, T), \\
-\lambda_{\nu,i}(t) \in N_{\mathbb{R}_-}(\mathbf{v}_i^T \mathbf{u}_{\mathcal{N}}(t)), \quad \forall i = 1, \dots, n_c \quad \text{a.e. in } (0, T), \\
\lambda_{\tau,i}(t) \in \mathcal{F} \lambda_{\nu,i}(t) \text{Sgn}(\mathbf{\tau}_i^T \dot{\mathbf{u}}_{\mathcal{N}^{\perp}}(t)), \quad \forall i = 1, \dots, n_c \quad \text{a.e. in } (0, T), \\
\mathbf{u}_{\mathcal{N}^{\perp}}(0) = \mathbf{u}_{\mathcal{N}^{\perp}}^0, \quad \dot{\mathbf{u}}_{\mathcal{N}^{\perp}}(0) = \mathbf{v}_{\mathcal{N}^{\perp}}^0,
\end{array}
\right.
\tag{8}
$$

where $N_{\mathbb{R}_-}$ denotes the normal cone of $\mathbb{R}_-$ and the multifunction Sgn: $\mathbb{R} \Rightarrow \mathbb{R}$ is the sub-differential of the function $a \mapsto |a|$, that is,

$$
\operatorname{Sgn}(a) = \begin{cases} \frac{a}{|a|} & \text{if } a \neq 0, \\ [-1, 1] & \text{if } a = 0, \end{cases} \quad a \in \mathbb{R}.
$$

REMARK 3.1. *From* (4), it immediately follows that there exists $\beta > 0$ such that

$$
\sup_{0 \neq w \in \mathbb{R}^{n_p}} \frac{\mathbf{w}^T \mathbf{B}_v^T \boldsymbol{\mu}_v + \mathbf{w}^T \mathbf{B}_{\tau}^T \boldsymbol{\mu}_{\tau}}{\|\mathbf{w}\|} \geq \beta (\|\boldsymbol{\mu}_v\| + \|\boldsymbol{\mu}_{\tau}\|), \quad \forall \boldsymbol{\mu}_v, \boldsymbol{\mu}_{\tau} \in \mathbb{R}^{n_c}. \quad (9)
$$

4. Well-posedness result

In this section, we shall establish the well-posedness of Problem (7). First, owing to (4) and (6), the first three variables of any $(\mathbf{u}_{\mathcal{N}^{\perp}}, \mathbf{u}_{\mathcal{N}}, \boldsymbol{\lambda}_{\nu}, \boldsymbol{\lambda}_{\tau})$ solving (7) have to satisfy

$$
\left\{
\begin{array}{l}
\mathbf{u}_{\mathcal{N}^{\perp}}(t) \in \mathcal{N}^{\perp}, \quad \mathbf{u}_{\mathcal{N}}(t) \in \mathcal{N}, \quad \boldsymbol{\lambda}_{\nu}(t) \in \mathbf{\Lambda}_{\nu}, \\
\mathbf{w}^T \mathbf{A}(\mathbf{u}_{\mathcal{N}^{\perp}}(t) + \mathbf{u}_{\mathcal{N}}(t)) = \mathbf{w}^T \mathbf{f} + \mathbf{w}^T \mathbf{B}_{\nu}^T \boldsymbol{\lambda}_{\nu}(t), \quad \forall \mathbf{w} \in \mathcal{N}, \\
(\boldsymbol{\mu}_{\nu} - \boldsymbol{\lambda}_{\nu}(t))^T \mathbf{B}_{\nu} \mathbf{u}_{\mathcal{N}}(t) \ge 0, \quad \forall \boldsymbol{\mu}_{\nu} \in \mathbf{\Lambda}_{\nu}
\end{array}
\right.
\tag{10}
$$

for almost all $t \in (0, T)$. From here, $\mathbf{u}_{\mathcal{N}}$ and $\boldsymbol{\lambda}_{\nu}$ are uniquely determined by $\mathbf{u}_{\mathcal{N}^{\perp}}$ as states the following assertion.

LEMMA 4.1. Let (4) be satisfied and $f \in \mathbb{R}^{n_p}$ be arbitrary. Then there exist unique functions $g_1: \mathcal{N}^{\perp} \to \mathcal{N}$ and $g_2: \mathcal{N}^{\perp} \to \mathbf{\Lambda}_{\nu}$ such that the triplet $(\mathbf{u}_{\mathcal{N}^{\perp}}(t), \mathbf{u}_{\mathcal{N}}(t), \boldsymbol{\lambda}_{\nu}(t))$ with $\mathbf{u}_{\mathcal{N}}(t) := g_1(\mathbf{u}_{\mathcal{N}^{\perp}}(t)), \boldsymbol{\lambda}_{\nu} := g_2(\mathbf{u}_{\mathcal{N}^{\perp}}(t))$, satisfies (10) for any $\mathbf{u}_{\mathcal{N}^{\perp}}(t) \in \mathcal{N}^{\perp}$ and any $t \in [0, T]$. Moreover, the functions $g_1$ and $g_2$ are Lipschitz continuous:

$$
\exists L_1, L_2 > 0: \|g_i(\mathbf{w}^1) - g_i(\mathbf{w}^2)\| \le L_i \|\mathbf{w}^1 - \mathbf{w}^2\|, \quad \forall \mathbf{w}^1, \mathbf{w}^2 \in \mathcal{N}^\perp, i = 1, 2. \quad (11)
$$

*Proof.* In fact, it suffices to analyse the static problem:

$$
\left\{
\begin{array}{l}
\text{Find } (\tilde{\mathbf{u}}_{\mathcal{N}}, \tilde{\boldsymbol{\lambda}}_{\nu}) := (\tilde{\mathbf{u}}_{\mathcal{N}}(\tilde{\mathbf{u}}_{\mathcal{N}^{\perp}}), \tilde{\boldsymbol{\lambda}}_{\nu}(\tilde{\mathbf{u}}_{\mathcal{N}^{\perp}})) \in \mathcal{N} \times \mathbf{\Lambda}_{\nu} \text{ such that} \\
\mathbf{w}^T \mathbf{A} \tilde{\mathbf{u}}_{\mathcal{N}} = \mathbf{w}^T (\mathbf{f} - \mathbf{A} \tilde{\mathbf{u}}_{\mathcal{N}^{\perp}}) + \mathbf{w}^T \mathbf{B}_{\nu}^T \tilde{\boldsymbol{\lambda}}_{\nu}, \quad \forall \mathbf{w} \in \mathcal{N}, \\
(\boldsymbol{\mu}_{\nu} - \tilde{\boldsymbol{\lambda}}_{\nu})^T \mathbf{B}_{\nu} \tilde{\mathbf{u}}_{\mathcal{N}} \ge 0, \quad \forall \boldsymbol{\mu}_{\nu} \in \mathbf{\Lambda}_{\nu},
\end{array}
\right.
\tag{12}
$$
for $\tilde{u}_{\mathcal{N}^{\perp}} \in \mathcal{N}^{\perp}$ given. It is readily seen that this problem is equivalent to finding a saddle-point $(\tilde{u}_{\mathcal{N}}, \lambda_{\nu})$ of the Lagrangian

$$ \mathcal{L}(w, \mu_v) := \frac{1}{2} w^T A w - w^T (f - A\tilde{u}_{\mathcal{N}^{\perp}}) - w^T B_v^T \mu_v, \quad (w, \mu_v) \in \mathbb{R}^{n_p} \times \mathbb{R}^{n_c}, $$

on $\mathcal{N} \times \Lambda_v$. Since $A$ is positive definite and

$$ \beta \|\mu_v\| \le \sup_{0 \ne w \in \mathbb{R}^{n_p}} \frac{w^T B_v^T \mu_v}{\|w\|} = \sup_{0 \ne w \in \mathcal{N}} \frac{w^T B_v^T \mu_v}{\|w\|}, \quad \forall \mu_v \in \mathbb{R}^{n_c} \quad (13) $$

due to (4), where $\beta$ is the constant from (9), Problem (12) possesses a unique solution for any $\tilde{u}_{\mathcal{N}^{\perp}} \in \mathcal{N}^{\perp}$, which depends Lipschitz continuously on the data $\tilde{u}_{\mathcal{N}^{\perp}}$ (see (22) for instance). This yields the existence, the uniqueness and the Lipschitz continuity of the functions $g_1$ and $g_2$.

From the other side, if $(u_{\mathcal{N}^{\perp}}, u_{\mathcal{N}}, \lambda_{\nu}, \lambda_{\tau})$ solves (8), then

$$ \left\{ \begin{array}{l} w^T M_r \ddot{u}_{\mathcal{N}^{\perp}}(t) + w^T A(u_{\mathcal{N}^{\perp}}(t) + u_{\mathcal{N}}(t)) = w^T f + w^T \left( \sum_{i=1}^{n_c} \lambda_{\tau,i}(t) \tau_i \right), \\ \qquad \forall w \in \mathcal{N}^{\perp} \text{ a.e. in } (0, T), \\[1.5ex] \lambda_{\tau,i}(t) \in \mathcal{F} \lambda_{\nu,i}(t) \text{Sgn}(\tau_i^T \dot{u}_{\mathcal{N}^{\perp}}(t)), \quad \forall i = 1, \dots, n_c \text{ a.e. in } (0, T), \\[1.5ex] u_{\mathcal{N}^{\perp}}(0) = u_{\mathcal{N}^{\perp}}^0, \quad \dot{u}_{\mathcal{N}^{\perp}}(0) = v_{\mathcal{N}^{\perp}}^0. \end{array} \right. \quad (14) $$

Substituting the inclusion for $\lambda_{\tau,i}(t)$ into the equality and taking $u_{\mathcal{N}}(t) = g_1(u_{\mathcal{N}^{\perp}}(t))$, $\lambda_{\nu,i}(t) = g_{2,i}(u_{\mathcal{N}^{\perp}}(t))$ according to Lemma 4.1, this becomes

$$ \left\{ \begin{array}{l} w^T M_r \ddot{u}_{\mathcal{N}^{\perp}}(t) \in w^T (f - A u_{\mathcal{N}^{\perp}}(t) - A g_1(u_{\mathcal{N}^{\perp}}(t))) \\ \qquad + w^T \left( \sum_{i=1}^{n_c} \mathcal{F} g_{2,i} (u_{\mathcal{N}^{\perp}}(t)) \text{Sgn}(\tau_i^T \dot{u}_{\mathcal{N}^{\perp}}(t)) \tau_i \right), \quad \forall w \in \mathcal{N}^{\perp} \text{ a.e. in } (0, T), \\ u_{\mathcal{N}^{\perp}}(0) = u_{\mathcal{N}^{\perp}}^0, \quad \dot{u}_{\mathcal{N}^{\perp}}(0) = v_{\mathcal{N}^{\perp}}^0. \end{array} \right. \quad (15) $$

**LEMMA 4.2.** Let (4) and (6) be fulfilled and $f \in \mathbb{R}^{n_p}$, $u_{\mathcal{N}^{\perp}}^0$, $v_{\mathcal{N}^{\perp}}^0 \in \mathcal{N}^{\perp}$ be arbitrary. Then there exists a unique Lipschitz continuous function $u_{\mathcal{N}^{\perp}}: [0, T] \to \mathcal{N}^{\perp}$ with $\ddot{u}_{\mathcal{N}^{\perp}} \in L^1(0, T; \mathbb{R}^{n_p})$ solving (15).

*Proof.* Introducing the matrix $P \in \mathcal{M}_{n_p, \bar{n}}(\mathbb{R})$, $\bar{n} := \dim \mathcal{N}^{\perp}$, columns of which form an orthonormal basis of $\mathcal{N}^{\perp}$, any vector $w \in \mathcal{N}^{\perp}$ can be represented by $\bar{w} \in \mathbb{R}^{\bar{n}}$ with

$$ \bar{w} = P^T w, \quad w = P P^T w = P \bar{w}, $$

and (15) is equivalent to

$$ \left\{ \begin{array}{l} \bar{w}^T \bar{M}_r \ddot{u}(t) \in \bar{w}^T (\bar{f} - \bar{A}\bar{u}(t) - \bar{g}_1(\bar{u}(t))) + \bar{w}^T \left( \sum_{i=1}^{n_c} \mathcal{F} \bar{g}_{2,i} (\bar{u}(t)) \text{Sgn}(\bar{\tau}_i^T \dot{\bar{u}}(t)) \bar{\tau}_i \right), \\ \qquad \forall \bar{w} \in \mathbb{R}^{\bar{n}} \text{ a.e. in } (0, T), \\[1.5ex] \bar{u}(0) = \bar{u}^0, \quad \dot{\bar{u}}(0) = P^T v_{\mathcal{N}^{\perp}}^0, \end{array} \right. $$

where

$$ \begin{align*} \bar{M}_r &= P^T M_r P, \quad \bar{A} = P^T A P, \quad \bar{g}_1(\bar{u}(t)) = P^T A g_1(P\bar{u}(t)), \\ \bar{g}_2(\bar{u}(t)) &= g_2(P\bar{u}(t)), \quad \bar{g}_2(\bar{u}(t)) = (\bar{g}_{2,1}(\bar{u}(t)), \dots, \bar{g}_{2,n_c}(\bar{u}(t)))^T, \\ \bar{u} &= P^T u_{\mathcal{N}^{\perp}}, \quad \bar{u}^0 = P^T u_{\mathcal{N}^{\perp}}^0, \quad \bar{f} = P^T f, \quad \bar{\tau}_i = P^T \tau_i, \quad i = 1, \dots, n_c. \end{align*} $$
Having in mind (6), this can be written as

$$
\left\{
\begin{array}{l}
\ddot{\bar{u}}(t) \in \bar{M}_r^{-1} \left[ \bar{f} - \bar{A}\bar{u}(t) - \bar{g}_1(\bar{u}(t)) + \sum_{i=1}^{n_c} \mathcal{F}\bar{g}_{2,i}(\bar{u}(t)) \text{Sgn}(\bar{\tau}_i^T \dot{\bar{u}}(t)) \bar{\tau}_i \right] \quad \text{a.e. in } (0, T), \\
\dot{\bar{u}}(0) = \bar{u}^0, \quad \dot{\bar{u}}(0) = \mathbf{P}^T \mathbf{v}_{\mathcal{N}^{\perp}}^0,
\end{array}
\right.
$$

and denoting $\bar{v} = \bar{M}_r^{1/2} \dot{\bar{u}}$, $\bar{v}^0 = \bar{M}_r^{1/2} \mathbf{P}^T \mathbf{v}_{\mathcal{N}^{\perp}}^0$, this leads to the following differential inclusion of the first order:

$$
\left\{
\begin{array}{l}
\left( \begin{pmatrix} \dot{\bar{u}}(t) \\ \dot{\bar{v}}(t) \end{pmatrix} \right) \in \left( \begin{pmatrix} \bar{M}_r^{-1/2} \bar{v}(t) \\ \bar{M}_r^{-1/2} \left[ \bar{f} - \bar{A}\bar{u}(t) - \bar{g}_1(\bar{u}(t)) + \sum_{i=1}^{n_c} \mathcal{F}\bar{g}_{2,i}(\bar{u}(t)) \text{Sgn} \left( \bar{\tau}_i^T \bar{M}_r^{-1/2} \bar{v}(t) \right) \bar{\tau}_i \right] \end{pmatrix} \right) \\
\text{a.e. in } (0, T), \\
\begin{pmatrix} \bar{u}(0) \\ \bar{v}(0) \end{pmatrix} = \begin{pmatrix} \bar{u}^0 \\ \bar{v}^0 \end{pmatrix}.
\end{array}
\right.
$$

Thus, we have to solve

$$
\left\{
\begin{array}{ll}
\dot{\mathbf{x}}(t) \in \mathbf{h}(\mathbf{x}(t)) & \text{a.e. in } (0, T), \\
\mathbf{x}(0) = \mathbf{x}^0 &
\end{array}
\right.
\qquad (16)
$$

with the multifunction $h: \mathbb{R}^{2\bar{n}} \Rightarrow \mathbb{R}^{2\bar{n}}$ defined by

$$
\mathbf{h}(\mathbf{y}) = \left( \bar{M}_r^{-1/2} \left[ \bar{f} - \bar{A}\mathbf{y}_1 - \bar{g}_1(\mathbf{y}_1) + \sum_{i=1}^{n_c} \mathcal{F}\bar{g}_{2,i}(\mathbf{y}_1) \text{Sgn} \left( \bar{\tau}_i^T \bar{M}_r^{-1/2} \mathbf{y}_2 \right) \bar{\tau}_i \right] \right),
$$
$$
\mathbf{y} = \begin{pmatrix} \mathbf{y}_1 \\ \mathbf{y}_2 \end{pmatrix} \in \mathbb{R}^{2\bar{n}},
$$

and $\mathbf{x}^0 = ((\bar{\mathbf{u}}^0)^T, (\bar{\mathbf{v}}^0)^T)^T$.

Obviously, $h$ is upper semi-continuous, that is, $h^{-1}(C)$ is closed whenever $C \subset \mathbb{R}^{2\bar{n}}$ is closed and has closed convex values. Furthermore, there exists $c > 0$ such that

$$
\|h(\mathbf{y})\| \equiv \sup\{\|z\| \mid z \in h(\mathbf{y})\} \le c(1 + \|\mathbf{y}\|), \quad \forall \mathbf{y} \in \mathbb{R}^{2\bar{n}}. \quad (17)
$$

Indeed,

$$
\begin{align*}
\|\mathbf{h}(\mathbf{y})\| &\le \|\bar{M}_r^{-1/2}\| \left[ \|\mathbf{y}_2\|^2 + \left\| \bar{f} - \bar{A}\mathbf{y}_1 - \bar{g}_1(\mathbf{y}_1) + \sum_{i=1}^{n_c} \mathcal{F}\bar{g}_{2,i}(\mathbf{y}_1) \text{Sgn} \left( \bar{\tau}_i^T \bar{M}_r^{-1/2} \mathbf{y}_2 \right) \bar{\tau}_i \right\|^2 \right]^{1/2} \\
&\le \|\bar{M}_r^{-1/2}\| \times \left[ \|\mathbf{y}_2\|^2 + (\|\bar{f}\| + \|\bar{A}\|\|\mathbf{y}_1\| + \|\bar{g}_1(\mathbf{y}_1)\| \right. \\
&\qquad \left. + \left\| \sum_{i=1}^{n_c} \mathcal{F}\bar{g}_{2,i}(\mathbf{y}_1) \text{Sgn} \left( \bar{\tau}_i^T \bar{M}_r^{-1/2} \mathbf{y}_2 \right) \bar{\tau}_i \right\|^2 \right]^{1/2}.
\end{align*}
$$
First,

$$ \left\| \sum_{i=1}^{n_c} \mathcal{F} \bar{g}_{2,i}(\mathbf{y}_1) \operatorname{Sgn} \left( \bar{\tau}_i^T \bar{M}_r^{-1/2} \mathbf{y}_2 \right) \bar{\tau}_i \right\| \le \left( \sum_{i=1}^{n_c} (\mathcal{F} \bar{g}_{2,i}(\mathbf{y}_1))^2 \right)^{1/2} = \mathcal{F} \| \bar{g}_2(\mathbf{y}_1) \| $$

in virtue of the orthonormality of $\bar{\tau}_i$ and the definition of the mapping Sgn. Second, making use of (11) and of the form of $\mathbf{P}$, we have

$$ \begin{aligned} \|\bar{g}_1(\mathbf{y}_1)\| &= \|\mathbf{P}^T \mathbf{A} g_1(\mathbf{P}\mathbf{y}_1)\| \le \|\mathbf{A}\| \|g_1(\mathbf{P}\mathbf{y}_1)\|, \\ \|g_1(\mathbf{P}\mathbf{y}_1)\| - \|g_1(\mathbf{P}\mathbf{0})\| &\le L_1 \|\mathbf{P}(\mathbf{y}_1 - \mathbf{0})\| = L_1 \|\mathbf{y}_1\|, \end{aligned} $$

consequently

$$ \|\bar{g}_1(\mathbf{y}_1)\| \le \|\mathbf{A}\|(\|g_1(\mathbf{0})\| + L_1\|\mathbf{y}_1\|), $$

and in a similar way, one can show that

$$ \|\bar{g}_2(\mathbf{y}_1)\| \le \|g_2(\mathbf{0})\| + L_2\|\mathbf{y}_1\|. $$

Hence,

$$ \|\mathbf{h}(\mathbf{y})\| \le \|\bar{\mathbf{M}}_r^{-1/2}\| \times \left[ \|\mathbf{y}_2\|^2 + (\|\bar{\mathbf{f}}\| + \|\bar{\mathbf{A}}\|\|\mathbf{y}_1\| + \|\mathbf{A}\|(\|g_1(\mathbf{0})\| + L_1\|\mathbf{y}_1\|) \right. \\ \left. + \mathcal{F}(\|g_2(\mathbf{0})\| + L_2\|\mathbf{y}_1\|) \right]^{1/2} $$

from which the expression for the constant $c$ in (17) follows. Therefore, Theorem 5.1 in (23) guarantees that (16) has an absolutely continuous solution $\mathbf{x}$ in $[0, T]$ for any $\mathbf{x}^0 \in \mathbb{R}^{2\bar{n}}$, that is, a function $\mathbf{x}: [0, T] \to \mathbb{R}^{2\bar{n}}$ with $\dot{\mathbf{x}} \in L^1(0, T; \mathbb{R}^{2\bar{n}})$ satisfying

$$ \mathbf{x}(t) = \mathbf{x}^0 + \int_0^t \dot{\mathbf{x}}(s) \, ds \quad \text{for all } t \in [0, T] \quad \text{and} \quad \dot{\mathbf{x}}(t) \in \mathbf{h}(\mathbf{x}(t)) \quad \text{a.e. in } (0, T). $$

This gives the existence part of the assertion. To prove the uniqueness, it suffices to show that $\mathbf{h}$ is one-sided Lipschitz (see for instance Theorem 10.4 in (23)), that is,

$$ \exists k > 0: (\mathbf{y}^1 - \mathbf{y}^2)^T (\mathbf{h}(\mathbf{y}^1) - \mathbf{h}(\mathbf{y}^2)) \le k \|\mathbf{y}^1 - \mathbf{y}^2\|^2, \quad \forall \mathbf{y}^1, \mathbf{y}^2 \in \mathbb{R}^{2\bar{n}}. $$

From the definition of $\mathbf{h}$

$$ \begin{aligned} (\mathbf{y}^1 - \mathbf{y}^2)^T (\mathbf{h}(\mathbf{y}^1) - \mathbf{h}(\mathbf{y}^2)) = & (\mathbf{y}_1^1 - \mathbf{y}_1^2)^T \bar{\mathbf{M}}_r^{-1/2} (\mathbf{y}_2^1 - \mathbf{y}_2^2) + (\mathbf{y}_2^1 - \mathbf{y}_2^2)^T \bar{\mathbf{M}}_r^{-1/2} \bar{\mathbf{A}} (\mathbf{y}_1^2 - \mathbf{y}_1^2) \\ & + (\mathbf{y}_2^1 - \mathbf{y}_2^2)^T \bar{\mathbf{M}}_r^{-1/2} (\bar{\mathbf{g}}_1 (\mathbf{y}_1^2) - \bar{\mathbf{g}}_1 (\mathbf{y}_1^1)) \\ & + (\mathbf{y}_2^1 - \mathbf{y}_2^2)^T \bar{\mathbf{M}}_r^{-1/2} \left( \sum_{i=1}^{n_c} \mathcal{F} (\bar{g}_{2,i} (\mathbf{y}_1^1) \operatorname{Sgn} (\bar{\tau}_i^T \bar{\mathbf{M}}_r^{-1/2} \mathbf{y}_2^1) \right. \\ & \left. \qquad - \bar{g}_{2,i} (\mathbf{y}_1^2) \operatorname{Sgn} (\bar{\tau}_i^T \bar{\mathbf{M}}_r^{-1/2} \mathbf{y}_2^2) \right) \bar{\tau}_i \right) \\ =: & S_1 + S_2 + S_3 + S_4. \end{aligned} $$
Clearly,

$$S_1 \le \|\bar{\mathbf{M}}_r^{-1/2}\| \|\mathbf{y}^1 - \mathbf{y}^2\|^2, \quad S_2 \le \|\bar{\mathbf{M}}_r^{-1/2} \bar{\mathbf{A}}\| \|\mathbf{y}^1 - \mathbf{y}^2\|^2$$

and

$$S_3 \le \|\bar{\mathbf{M}}_r^{-1/2}\| \|\mathbf{y}^1 - \mathbf{y}^2\| \|\bar{\mathbf{g}}_1(\mathbf{y}_1^2) - \bar{\mathbf{g}}_1(\mathbf{y}_1^1)\| \\ \le \|\mathbf{A}\| \|\bar{\mathbf{M}}_r^{-1/2}\| \|\mathbf{y}^1 - \mathbf{y}^2\| \|\mathbf{g}_1(\mathbf{P}\mathbf{y}_1^2) - \mathbf{g}_1(\mathbf{P}\mathbf{y}_1^1)\| \le L_1 \|\mathbf{A}\| \|\bar{\mathbf{M}}_r^{-1/2}\| \|\mathbf{y}^1 - \mathbf{y}^2\|^2$$

by (11). Furthermore,

$$S_4 = \sum_{i=1}^{n_c} \mathcal{F} \left( \bar{g}_{2,i}(\mathbf{y}_1^1) \operatorname{Sgn} \left( \bar{\tau}_i^T \bar{\mathbf{M}}_r^{-1/2} \mathbf{y}_2^1 \right) - \bar{g}_{2,i}(\mathbf{y}_1^2) \operatorname{Sgn} \left( \bar{\tau}_i^T \bar{\mathbf{M}}_r^{-1/2} \mathbf{y}_2^2 \right) \right) \\ \times \left( \mathbf{y}_2^1 T \bar{\mathbf{M}}_r^{-1/2} \bar{\tau}_i - \mathbf{y}_2^2 T \bar{\mathbf{M}}_r^{-1/2} \bar{\tau}_i \right).$$

Hence, fixing $i$ and setting

$$a^1 = \bar{g}_{2,i}(\mathbf{y}_1^1), \quad a^2 = \bar{g}_{2,i}(\mathbf{y}_1^2), \quad b^1 = \bar{\tau}_i^T \bar{\mathbf{M}}_r^{-1/2} \mathbf{y}_2^1, \quad b^2 = \bar{\tau}_i^T \bar{\mathbf{M}}_r^{-1/2} \mathbf{y}_2^2,$$

the $i$th summand of $S_4$ takes the form

$$\mathcal{F}(a^1 \operatorname{Sgn}(b^1) - a^2 \operatorname{Sgn}(b^2))(b^1 - b^2).$$

Note that $a^1, a^2 \le 0$. We claim that in this case,

$$(a^1 \operatorname{Sgn}(b^1) - a^2 \operatorname{Sgn}(b^2))(b^1 - b^2) \le |a^1 - a^2| |b^1 - b^2|. \quad (18)$$

Indeed, for $\zeta \in \operatorname{Sgn}(b^1)$ and $\xi \in \operatorname{Sgn}(b^2)$, we get

$$(a^1 \zeta - a^2 \xi)(b^1 - b^2) = (a^1 \zeta - a^1 \xi + a^1 \xi - a^2 \xi)(b^1 - b^2) \le (a^1 - a^2) \xi(b^1 - b^2)$$

due to the monotonicity of the multifunction Sgn. And of course (18) can be deduced from

$$(a^1 - a^2)\xi(b^1 - b^2) \le |a^1 - a^2| |b^1 - b^2|.$$

Applying this together with the Cauchy–Schwarz inequality and (11) we get

$$S_4 \le \mathcal{F} \sum_{i=1}^{n_c} |\bar{g}_{2,i}(\mathbf{y}_1^1) - \bar{g}_{2,i}(\mathbf{y}_1^2)| \|\bar{\tau}_i^T \bar{\mathbf{M}}_r^{-1/2} \mathbf{y}_2^1 - \bar{\tau}_i^T \bar{\mathbf{M}}_r^{-1/2} \mathbf{y}_2^2\| \\ \le \mathcal{F} \|\bar{g}_2(\mathbf{y}_1^1) - \bar{g}_2(\mathbf{y}_1^2)\| \|\mathbf{B}_\tau \bar{\mathbf{M}}_r^{-1/2}(\mathbf{y}_2^1 - \mathbf{y}_2^2)\| \le \mathcal{F} L_2 \|\bar{\mathbf{M}}_r^{-1/2}\| \|\mathbf{y}^1 - \mathbf{y}^2\|^2.$$

All in all, the one-sided Lipschitz property of $h$ is verified.

On the basis of the previous two lemmas, we arrive at the announced well-posedness result.

**THEOREM 4.3.** Let $f \in \mathbb{R}^{n_p}$, $u_{\mathcal{N}^{\perp}}^0$, $v_{\mathcal{N}^{\perp}}^0 \in \mathcal{N}^{\perp}$ be arbitrary. If (4) and (6) are satisfied, then there exist a unique Lipschitz continuous function $u_{\mathcal{N}^{\perp}}: [0, T] \to \mathcal{N}^{\perp}$ with $\ddot{u}_{\mathcal{N}^{\perp}} \in L^1(0, T; \mathbb{R}^{n_p})$ and unique functions $u_{\mathcal{N}}: [0, T] \to \mathcal{N}$ and $\lambda_{\nu}, \lambda_{\tau}: [0, T] \to \mathbb{R}^{n_c}$ such that the quadruplet $(u_{\mathcal{N}^{\perp}}, u_{\mathcal{N}}, \lambda_{\nu}, \lambda_{\tau})$ solves (7). In addition, $u_{\mathcal{N}}, \lambda_{\nu}$ are Lipschitz continuous in $[0, T]$ and $\lambda_{\tau} \in L^{\infty}(0, T; \mathbb{R}^{n_c})$.
*Proof.* The existence and uniqueness as well as the Lipschitz continuity of $u_{\mathcal{N}^{\perp}}$ and $u_{\mathcal{N}}$, $\lambda_{\nu}$ are ensured by Lemmas 4.2 and 4.1, respectively. Consequently, the existence of $\lambda_{\tau}$ is readily seen from the relation between (14) and (15). If $(u_{\mathcal{N}^{\perp}}, u_{\mathcal{N}}, \lambda_{\nu}, \lambda_{\tau}^{1})$ and $(u_{\mathcal{N}^{\perp}}, u_{\mathcal{N}}, \lambda_{\nu}, \lambda_{\tau}^{2})$ were two solutions to (7), then

$$ w^T B_{\tau}^T (\lambda_{\tau}^1(t) - \lambda_{\tau}^2(t)) = 0, \quad \forall w \in \mathbb{R}^{n_p}, \text{ a.e. in } (0, T) $$

by the first equation in (7) and

$$ \beta \|\lambda_{\tau}^1(t) - \lambda_{\tau}^2(t)\| \le \sup_{0 \ne w \in \mathbb{R}^{n_p}} \frac{w^T B_{\tau}^T (\lambda_{\tau}^1(t) - \lambda_{\tau}^2(t))}{\|w\|} = 0 \quad \text{a.e. in } (0, T) $$

due to (9). In a similar way, one also obtains that $\lambda_{\tau} \in L^{\infty}(0, T; \mathbb{R}^{n_c})$ from the Lipschitz continuity of $\lambda_{\nu}$ and the second inclusion of (8).

**REMARK 4.4.** The well-posedness result can easily be extended to three-dimensional (3D) problems since the key point of the proof is the monotonicity of the multifunction Sgn. For the 3D problems, the friction condition can be expressed by means of the subdifferential of the function $a \mapsto \|a\|$, which is also monotonic. This allows to obtain an equivalent relation to (18). An extension of the well-posedness result can also be obtained for a load vector which is a Lipschitz continuous function of time.

**REMARK 4.5.** The spatial semi-discrete problem (7) being equivalent to the one-sided Lipschitz regular differential inclusion (16), most of the classical time integration schemes will be convergent (for a fixed mesh) due to, for instance, the result obtained in (24). Moreover, the fully discrete problem is also ensured to be well-posed for a sufficiently small time step (because of its monotonicity).

## 5. Energy decreasing result

First, note that the result of Proposition 1 in (13) is still valid, which means that the so-called persistency condition holds

$$ \lambda_{\nu,i}(t)(v_i^T \dot{u}(t)) = 0 \quad \text{a.e. in } (0, T), \quad i = 1, \dots, n_c. $$

This allows to prove the following result.

**PROPOSITION 5.1** Still assuming the load vector $f$ to be constant in time and denoting $u = u_{\mathcal{N}} + u_{\mathcal{N}^{\perp}}$ the solution to (7), the energy

$$ E(t) = \frac{1}{2} \dot{u}^T(t) M_r \dot{u}(t) + \frac{1}{2} u^T(t) Au(t) - u^T(t) f \quad (19) $$

is decreasing in time.

*Proof.* The first equation in (8) implies

$$ \dot{u}^T(t) M_r \ddot{u}(t) + \dot{u}^T(t) Au(t) = \dot{u}^T(t) f + \sum_{i=1}^{n_c} \lambda_{\nu,i}(t) \dot{u}^T(t) v_i + \sum_{i=1}^{n_c} \lambda_{\tau,i}(t) \dot{u}^T(t) \tau_i. $$

Integrating from $t_0$ to $t_1$, it follows:

$$ E(t_1) = E(t_0) + \sum_{i=1}^{n_c} \int_{t_0}^{t_1} \lambda_{\nu,i}(t) \dot{u}^T(t) v_i \, dt + \sum_{i=1}^{n_c} \int_{t_0}^{t_1} \lambda_{\tau,i}(t) \dot{u}^T(t) \tau_i \, dt. $$
Due to the fact that $\lambda_{\tau,i}(t) \in \mathcal{F}\lambda_{\nu,i}(t)\text{Sgn}(\dot{u}^T(t)\tau_i)$ a.e. in $(0, T)$, one has

$$ \int_{t_0}^{t_1} \lambda_{\tau,i}(t) \dot{u}^T(t) \tau_i dt \le 0, \quad i = 1, \dots, n_c. $$

Together with the persistency condition, this gives the result.

## 6. An elementary example

This section presents the mass redistribution method for an elementary contact problem involving a single linear triangular finite-element depicted in Fig. 2. The aim is to show that an undifferentiated treatment of the contact and friction condition may lead to an ill-posedness of the fully discrete problem whatever is the length of the time step. Using the time discretization by the midpoint scheme, we shall compare different possibilities of the redistribution of the mass. The studied contact problem is in fact a dynamic case of the elementary example studied, for example, in (25) and (26).

Denoting the lengths of the sides of the triangle by $\ell$, $\ell$, $\sqrt{2}\ell$ and employing Hooke's constitutive law for homogeneous, isotropic material, the formulation of the problem in inclusions reads as follows:

$$ \left\{ \begin{array}{l} \text{Find } \mathbf{u}: [0, T] \to \mathbb{R}^2, \lambda_{\nu}, \lambda_{\tau}: [0, T] \to \mathbb{R} \text{ such that} \\ \mathbf{M}\ddot{\mathbf{u}}(t) + \mathbf{A}\mathbf{u}(t) = \mathbf{f}(t) + \lambda_{\nu}(t)\mathbf{v} + \lambda_{\tau}(t)\mathbf{\tau} \quad \text{a.e. in } (0, T), \\ -\lambda_{\nu}(t) \in N_{\mathbb{R}_-}(\mathbf{v}^T\mathbf{u}(t)) \quad \text{a.e. in } (0, T), \\ \lambda_{\tau}(t) \in \mathcal{F}\lambda_{\nu}(t)\text{Sgn}(\mathbf{\tau}^T\dot{\mathbf{u}}(t)) \quad \text{a.e. in } (0, T), \\ \mathbf{u}(0) = \mathbf{u}^0, \quad \dot{\mathbf{u}}(0) = \mathbf{v}^0, \end{array} \right. $$

where

$$ M = \begin{pmatrix} \frac{\rho \ell^2}{12} & 0 \\ 0 & \frac{\rho \ell^2}{12} \end{pmatrix}, \quad A = \begin{pmatrix} \frac{\lambda + 3\mu}{2} & -\frac{\lambda + \mu}{2} \\ -\frac{\lambda + \mu}{2} & \frac{\lambda + 3\mu}{2} \end{pmatrix}, \quad \mathbf{v} = \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \quad \mathbf{\tau} = \begin{pmatrix} 0 \\ 1 \end{pmatrix}. $$

Here, $\rho > 0$ is constant, $\lambda \ge 0$, $\mu > 0$ are the Lamé coefficients and $f$ is assumed to be dependent on $t$.

Obviously, the mass redistribution method consists in replacing the matrix $M$ by $M_r = \begin{pmatrix} m_1 & 0 \\ 0 & m_2 \end{pmatrix}$ with $m_1, m_2 \ge 0$. Further, to do the time discretization by the midpoint method, we divide the

<figure><img src="image_3.png" /><figcaption>Fig. 2 Geometry of the elementary example</figcaption></figure>
interval $[0, T]$ uniformly into $n_t$ subintervals and set $\Delta t = T/n_t$ and $t_k = k\Delta t$ for $k = 0, \dots, n_t$. We seek the approximations $\mathbf{u}^{k+1}$ and $\mathbf{v}^{k+1}$ of $\mathbf{u}(t_{k+1})$ and $\dot{\mathbf{u}}(t_{k+1})$, respectively, for $k = 0, \dots, n_t - 1$ so that

$$ \begin{cases} \mathbf{u}^{k+1} = \mathbf{u}^k + \Delta t \mathbf{v}^{k+1/2}, & \mathbf{v}^{k+1} = \mathbf{v}^k + \Delta t \mathbf{a}^{k+1/2}, \\ \mathbf{u}^{k+1/2} = \frac{\mathbf{u}^{k+1} + \mathbf{u}^k}{2}, & \mathbf{v}^{k+1/2} = \frac{\mathbf{v}^{k+1} + \mathbf{v}^k}{2} \end{cases} \qquad (20) $$

and

$$ \left\{ \begin{array}{l} \mathbf{M}_r \mathbf{a}^{k+1/2} + \mathbf{A} \mathbf{u}^{k+1/2} = \mathbf{f}\left(\frac{t_k+t_{k+1}}{2}\right) + \lambda_{\nu}^{k+1/2} \mathbf{v} + \lambda_{\tau}^{k+1/2} \mathbf{\tau}, \\ -\lambda_{\nu}^{k+1/2} \in N_{\mathbb{R}_-}(\mathbf{v}^T \mathbf{u}^{k+1/2}), \\ \lambda_{\tau}^{k+1/2} \in \mathcal{F} \lambda_{\nu}^{k+1/2} \text{Sgn}(\mathbf{\tau}^T \mathbf{v}^{k+1/2}). \end{array} \right. \qquad (21) $$

From (20) one can express $\mathbf{v}^{k+1/2}$ and $\mathbf{a}^{k+1/2}$ as

$$ \mathbf{v}^{k+1/2} = \frac{2}{\Delta t} \mathbf{u}^{k+1/2} - \frac{2}{\Delta t} \mathbf{u}^k, \quad \mathbf{a}^{k+1/2} = -\frac{4}{\Delta t^2} \mathbf{u}^{k+1/2} - \frac{4}{\Delta t^2} \mathbf{u}^k - \frac{2}{\Delta t} \mathbf{v}^k, \qquad (22) $$

which inserted into (21) leads to

$$ \left\{ \begin{array}{l} \left( \frac{4}{\Delta t^2} \mathbf{M}_r + \mathbf{A} \right) \mathbf{u}^{k+1/2} = \hat{\mathbf{f}}^{k+1/2} + \lambda_{\nu}^{k+1/2} \mathbf{v} + \lambda_{\tau}^{k+1/2} \mathbf{\tau}, \\ -\lambda_{\nu}^{k+1/2} \in N_{\mathbb{R}_-} (\mathbf{v}^T \mathbf{u}^{k+1/2}), \\ \lambda_{\tau}^{k+1/2} \in \mathcal{F} \lambda_{\nu}^{k+1/2} \text{Sgn} \left( \frac{2}{\Delta t} \mathbf{\tau}^T (\mathbf{u}^{k+1/2} - \mathbf{u}^k) \right) \end{array} \right. $$

with

$$ \hat{\mathbf{f}}^{k+1/2} = \mathbf{f}\left(\frac{t_k + t_{k+1}}{2}\right) + \frac{4}{\Delta t^2} \mathbf{M}_r \mathbf{u}^k + \frac{2}{\Delta t} \mathbf{M}_r \mathbf{v}^k. $$

Finally, we consider the decomposition

$$ \mathbf{u}^i = (u_{\nu}^i, u_{\tau}^i), \mathbf{v}^i = (v_{\nu}^i, v_{\tau}^i), \mathbf{a}^i = (a_{\nu}^i, a_{\tau}^i), \hat{\mathbf{f}}^i = (\hat{f}_{\nu}^i, \hat{f}_{\tau}^i) $$

and denote

$$ a := \left( \frac{4}{\Delta t^2} m_1 + \frac{\lambda + 3\mu}{2} \right), \quad b := \frac{\lambda + \mu}{2}, \quad c := \left( \frac{4}{\Delta t^2} m_2 + \frac{\lambda + 3\mu}{2} \right). $$

In each time step, we then obtain the following problem:

$$ \left\{ \begin{array}{l} \text{Find } (u_{\nu}^{k+1/2}, u_{\tau}^{k+1/2}, \lambda_{\nu}^{k+1/2}, \lambda_{\tau}^{k+1/2}) \in \mathbb{R}^4 \text{ such that} \\ au_{\nu}^{k+1/2} - bu_{\tau}^{k+1/2} = \hat{f}_{\nu}^{k+1/2} + \lambda_{\nu}^{k+1/2}, \\ -bu_{\nu}^{k+1/2} + cu_{\tau}^{k+1/2} = \hat{f}_{\tau}^{k+1/2} + \lambda_{\tau}^{k+1/2}, \\ -\lambda_{\nu}^{k+1/2} \in N_{\mathbb{R}_-}(u_{\nu}^{k+1/2}), \\ \lambda_{\tau}^{k+1/2} \in \mathcal{F} \lambda_{\nu}^{k+1/2} \text{Sgn}(u_{\tau}^{k+1/2} - u_{\tau}^k), \end{array} \right. \qquad (23) $$
after resolution of which the values of $(u_v^{k+1}, u_\tau^{k+1})$ and $(v_v^{k+1}, v_\tau^{k+1})$ are determined by (20) and (22).

We shall derive exact solutions of problem (23) for an arbitrary $k \in \{0, \dots, n_t - 1\}$ by considering all possible situations occurring in the inclusions $(23)_4$ and $(23)_5$.

(i) Let $\lambda_v^{k+1/2} = 0$. From (23)$_5$, it follows that $\lambda_\tau^{k+1/2} = 0$ and solving the equations (23)$_2$ and (23)$_3$ we get

$$u_v^{k+1/2} = \frac{c \hat{f}_v^{k+1/2} + b \hat{f}_\tau^{k+1/2}}{ac - b^2}, \quad u_\tau^{k+1/2} = \frac{a \hat{f}_\tau^{k+1/2} + b \hat{f}_v^{k+1/2}}{ac - b^2}.$$

Since $u_v^{k+1/2} \le 0$ by (23)$_4$, this solution is valid under the following constraint:

$$c \hat{f}_v^{k+1/2} + b \hat{f}_\tau^{k+1/2} \le 0.$$

(ii) If $\lambda_v^{k+1/2} < 0$ and $u_\tau^{k+1/2} = u_\tau^k$, then $u_v^{k+1/2} = 0$ according to (23)$_4$ and (23)$_{2,3}$ yield

$$\lambda_v^{k+1/2} = -bu_\tau^k - \hat{f}_v^{k+1/2}, \quad \lambda_\tau^{k+1/2} = cu_\tau^k - \hat{f}_\tau^{k+1/2}.$$

Our assumption $\lambda_v^{k+1/2} < 0$ and the condition $\mathcal{F}\lambda_v^{k+1/2} \le \lambda_\tau^{k+1/2} \le -\mathcal{F}\lambda_\tau^{k+1/2}$ implied by (23)$_5$ give the restrictions:

$$\hat{f}_v^{k+1/2} > -bu_\tau^k, \quad (c - b\mathcal{F})u_\tau^k - \mathcal{F}\hat{f}_v^{k+1/2} \le \hat{f}_\tau^{k+1/2} \le (c + b\mathcal{F})u_\tau^k + \mathcal{F}\hat{f}_v^{k+1/2}.$$

(iii) Let us pose $\lambda_v^{k+1/2} < 0$, $u_\tau^{k+1/2} > u_\tau^k$. Making use of (23)$_{4,5}$, we have $u_v^{k+1/2} = 0$, $\lambda_\tau^{k+1/2} = \mathcal{F}\lambda_v^{k+1/2}$ and (23)$_{2,3}$ lead to

$$u_\tau^{k+1/2} = \frac{\hat{f}_\tau^{k+1/2} - \mathcal{F}\hat{f}_v^{k+1/2}}{c + b\mathcal{F}}, \quad \lambda_v^{k+1/2} = -\frac{c\hat{f}_v^{k+1/2} + b\hat{f}_\tau^{k+1/2}}{c + b\mathcal{F}}.$$

From $\lambda_v^{k+1/2} < 0$ and $u_\tau^{k+1/2} > u_\tau^k$, one can see that

$$c \hat{f}_v^{k+1/2} + b \hat{f}_\tau^{k+1/2} > 0, \quad \hat{f}_\tau^{k+1/2} > (c + b\mathcal{F})u_\tau^k + \mathcal{F}\hat{f}_v^{k+1/2}.$$

(iv) Suppose that $\lambda_v^{k+1/2} < 0$ and $u_\tau^{k+1/2} < u_\tau^k$. Consequently, $u_v^{k+1/2} = 0$ and $\lambda_\tau^{k+1/2} = -\mathcal{F}\lambda_v^{k+1/2}$. If $\mathcal{F} \ne c/b$, then we obtain

$$u_\tau^{k+1/2} = \frac{\hat{f}_\tau^{k+1/2} + \mathcal{F}\hat{f}_v^{k+1/2}}{c - b\mathcal{F}}, \quad \lambda_v^{k+1/2} = -\frac{c\hat{f}_v^{k+1/2} + b\hat{f}_\tau^{k+1/2}}{c - b\mathcal{F}}.$$

under the condition

$$\mathcal{F} < \frac{c}{b}, \quad c \hat{f}_v^{k+1/2} + b \hat{f}_\tau^{k+1/2} > 0, \quad \hat{f}_\tau^{k+1/2} < (c - b\mathcal{F})u_\tau^k - \mathcal{F}\hat{f}_v^{k+1/2}$$

or

$$\mathcal{F} > \frac{c}{b}, \quad c \hat{f}_v^{k+1/2} + b \hat{f}_\tau^{k+1/2} < 0, \quad \hat{f}_\tau^{k+1/2} > (c - b\mathcal{F})u_\tau^k - \mathcal{F}\hat{f}_v^{k+1/2}.$$

In the case of $\mathcal{F} = c/b$, there exists the whole solution set

$$\left\{ (u_\tau^{k+1/2}, \lambda_v^{k+1/2}) \in \mathbb{R}^2 \mid \lambda_v^{k+1/2} = -bu_\tau^{k+1/2} - \hat{f}_v^{k+1/2} \right\}$$

with the following restrictions:

$$\hat{f}_\tau^{k+1/2} + \mathcal{F}\hat{f}_v^{k+1/2} = 0, \quad -\frac{\hat{f}_v^{k+1/2}}{b} < u_\tau^{k+1/2} < u_\tau^k, \quad -bu_\tau^k < \hat{f}_v^{k+1/2}.$$
To summarize the results, introduce the linear functions $S_i^{k+1/2}: \mathbb{R}^2 \times (0, +\infty) \to \mathbb{R}^4$, $1 \le i \le 4$, and the multi-valued function $S_5^{k+1/2}: \mathbb{R}^2 \times (0, +\infty) \to \mathbb{R}^4$ by

$$S_1^{k+1/2}(\hat{f}, \mathcal{F}) = \left( \frac{c\hat{f}_v + b\hat{f}_\tau}{ac - b^2}, \frac{a\hat{f}_\tau + b\hat{f}_v}{ac - b^2}, 0, 0 \right)^T, \quad \hat{f} \in \mathbb{R}^2, \mathcal{F} > 0,$$

$$S_2^{k+1/2}(\hat{f}, \mathcal{F}) = (0, u_\tau^k, -(\hat{f}_v + bu_\tau^k), cu_\tau^k - \hat{f}_\tau)^T, \quad \hat{f} \in \mathbb{R}^2, \mathcal{F} > 0,$$

$$S_3^{k+1/2}(\hat{f}, \mathcal{F}) = \left( 0, \frac{\hat{f}_\tau - \mathcal{F}\hat{f}_v}{c + b\mathcal{F}}, -\frac{c\hat{f}_v + b\hat{f}_\tau}{c + b\mathcal{F}}, -\mathcal{F}\frac{c\hat{f}_v + b\hat{f}_\tau}{c + b\mathcal{F}} \right)^T, \quad \hat{f} \in \mathbb{R}^2, \mathcal{F} > 0,$$

$$S_4^{k+1/2}(\hat{f}, \mathcal{F}) = \left( 0, \frac{\hat{f}_\tau + \mathcal{F}\hat{f}_v}{c - b\mathcal{F}}, -\frac{c\hat{f}_v + b\hat{f}_\tau}{c - b\mathcal{F}}, \mathcal{F}\frac{c\hat{f}_v + b\hat{f}_\tau}{c - b\mathcal{F}} \right)^T, \quad \hat{f} \in \mathbb{R}^2, \mathcal{F} \in (0, +\infty) \setminus \left\{ \frac{c}{b} \right\},$$

$$S_5^{k+1/2}(\hat{f}, \mathcal{F}) = \left\{ (u_\nu, u_\tau, \lambda_\nu, \lambda_\tau)^T \in \mathbb{R}^4 \mid u_\nu = 0, -\frac{\hat{f}_\nu}{b} \le u_\tau \le u_\tau^k, \right. \\ \left. \qquad \lambda_\nu = -(\hat{f}_\nu + bu_\tau), \lambda_\tau = \mathcal{F}(\hat{f}_\nu + bu_\tau) \right\}, \quad \hat{f} \in \mathbb{R}^2, \mathcal{F} = \frac{c}{b},$$

and for $\mathcal{F} > 0$ define the sets:

$$\sigma_1^{k+1/2}(\mathcal{F}) = \{\hat{f} \in \mathbb{R}^2 \mid c\hat{f}_v + b\hat{f}_\tau \le 0\},$$

$$\sigma_2^{k+1/2}(\mathcal{F}) = \{\hat{f} \in \mathbb{R}^2 \mid \hat{f}_v \ge -bu_\tau^k, (c - b\mathcal{F})u_\tau^k - \mathcal{F}\hat{f}_v \le \hat{f}_\tau \le (c + b\mathcal{F})u_\tau^k + \mathcal{F}\hat{f}_v\},$$

$$\sigma_3^{k+1/2}(\mathcal{F}) = \{\hat{f} \in \mathbb{R}^2 \mid c\hat{f}_v + b\hat{f}_\tau \ge 0, \hat{f}_\tau \ge (c + b\mathcal{F})u_\tau^k + \mathcal{F}\hat{f}_v\},$$

$$\sigma_4^{k+1/2}(\mathcal{F}) = \begin{cases} \{\hat{f} \in \mathbb{R}^2 \mid \hat{f}_v \ge -bu_\tau^k, c\hat{f}_v + b\hat{f}_\tau \ge 0, \hat{f}_\tau \le (c - b\mathcal{F})u_\tau^k - \mathcal{F}\hat{f}_v\} & \text{if } \mathcal{F} \in (0, \frac{c}{b}), \\ \{\hat{f} \in \mathbb{R}^2 \mid \hat{f}_v \ge -bu_\tau^k, c\hat{f}_v + b\hat{f}_\tau \le 0, \hat{f}_\tau \ge (c - b\mathcal{F})u_\tau^k - \mathcal{F}\hat{f}_v\} & \text{if } \mathcal{F} \in (\frac{c}{b}, +\infty), \end{cases}$$

$$\sigma_5^{k+1/2}(\mathcal{F}) = \{\hat{f} \in \mathbb{R}^2 \mid \hat{f}_v \ge -bu_\tau^k, c\hat{f}_v + b\hat{f}_\tau = 0\}, \quad \mathcal{F} = \frac{c}{b}.$$

One can easily verify that $S_i^{k+1/2}(\hat{f}^{k+1/2}, \mathcal{F})$ solves (23) for $\hat{f}^{k+1/2} \in \sigma_i^{k+1/2}(\mathcal{F}), \mathcal{F} > 0$, $1 \le i \le 4$, and $S_5^{k+1/2}(\hat{f}^{k+1/2}, \mathcal{F})$ is the set of solutions to (23) for $\hat{f}^{k+1/2} \in \sigma_5^{k+1/2}(\mathcal{F}), \mathcal{F} = c/b$. Hence, the structure of the solution set to (23) depends on the mutual position of $\sigma_i^{k+1/2}(\mathcal{F})$, which depends on the magnitude of $\mathcal{F}$.

If $\mathcal{F} \in (0, c/b)$, then the interiors of $\sigma_i^{k+1/2}(\mathcal{F})$ are mutually disjoint for $1 \le i \le 4$ and

$$S_i^{k+1/2}(\hat{f}, \mathcal{F}) = S_j^{k+1/2}(\hat{f}, \mathcal{F}), \quad \forall \hat{f} \in \partial\sigma_i^{k+1/2}(\mathcal{F}) \cap \partial\sigma_j^{k+1/2}(\mathcal{F}), \quad \forall i, j \in \{1, \dots, 4\}.$$
Consequently, (23) has a unique solution for any $\hat{f}^{k+1/2} \in \mathbb{R}^2$ (see Fig. 3; note that $u_\nu^{k+1/2}, u_\tau^{k+1/2}$ and $\lambda_\tau^{k+1/2}$ are uniquely determined by the values of $\lambda_\nu^{k+1/2}$). If $\mathcal{F} > c/b$, then $\sigma_4^{k+1/2}(\mathcal{F}) = \sigma_1^{k+1/2}(\mathcal{F}) \cap \sigma_2^{k+1/2}(\mathcal{F})$ and its interior is non-empty. In this case, there exists a unique solution to (23) if $\hat{f}^{k+1/2} \in (\mathbb{R}^2 \setminus \sigma_4^{k+1/2}(\mathcal{F})) \cup \{(-bu_\tau^k, cu_\tau^k)\}$, there are two solutions on $\partial\sigma_4^{k+1/2}(\mathcal{F}) \setminus \{(-bu_\tau^k, cu_\tau^k)\}$ and three solutions in $\text{Int} \sigma_4^{k+1/2}(\mathcal{F})$ (see Fig. 4). Finally, if $\mathcal{F} = c/b$, $\sigma_1^{k+1/2}(\mathcal{F}) \cap \sigma_2^{k+1/2}(\mathcal{F}) = \sigma_5^{k+1/2}(\mathcal{F})$ is a half-line and there exists a unique solution to (23) for $\hat{f}^{k+1/2} \in (\mathbb{R}^2 \setminus \sigma_5^{k+1/2}(\mathcal{F})) \cup \{(-bu_\tau^k, cu_\tau^k)\}$, whereas the continuous branch $S_5^{k+1/2}(\hat{f}^{k+1/2}, \mathcal{F})$ of solutions connects $S_1^{k+1/2}(\hat{f}^{k+1/2}, \mathcal{F})$ and $S_2^{k+1/2}(\hat{f}^{k+1/2}, \mathcal{F})$ for $\hat{f}^{k+1/2} \in \sigma_5^{k+1/2}(\mathcal{F}) \setminus \{(-bu_\tau^k, cu_\tau^k)\}$ (as depicted in Fig. 5).

<figure><img src="image_4.png" /><figcaption>Fig. 3 Structure of the solution for $\mathcal{F} \in (0, c/b)$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4 Structure of the solution for $\mathcal{F} \in (c/b, +\infty)$</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 5 Structure of the solution for $\mathcal{F} = c/b$</figcaption></figure>

Now take the redistributed mass matrix $M_r$ such that $m_1 = 0$ and $m_2 > 0$, that is, (6) is fulfilled. Then, for any $\mathcal{F} > 0$ given, one can find $\Delta t_0 > 0$ satisfying

$$ \frac{c}{b} = \frac{\frac{4}{\Delta t^2}m_2 + \frac{\lambda+3\mu}{2}}{\frac{\lambda+\mu}{2}} > \mathcal{F}, \quad \forall \Delta t \in (0, \Delta t_0) $$

and the analysis above ensures the unique solvability of (23) for any $\int_{\tau}^{k+1/2} \in \mathbb{R}^2$ and any $\Delta t \in (0, \Delta t_0)$. Observe that this is in good accordance with the well-posedness result established in Section 4.

On the contrary, consider $M_r$ with $m_1 > 0$, $m_2 = 0$ or $m_1 = m_2 = 0$, which corresponds to the elimination of the mass in the tangential direction and the total elimination of the mass on the contact zone, respectively. If the coefficient $\mathcal{F}$ is larger than $(\lambda + 3\mu)/(\lambda + \mu) = c/b$, one can always find $\int_{\tau}^{k+1/2}$ such that (23) possesses multiple solutions whatever small $\Delta t$ is. This suggests that the well-posedness is not reached in such cases.

## 7. Numerical tests

The numerical simulations presented in (13) show the effectiveness of the mass redistribution method to remove the spurious oscillations caused by the contact condition in discretized dynamical problems. Here, we chose a test case where the sliding is much more present with still a partial loss of contact (see Fig. 6). An elastic body whose reference configuration is the square $\{(x_1, x_2) \in ]0, 10 \text{ cm} \times ]0, 10 \text{ cm}]\}$ is in contact with a rigid foundation at $x_2 = 0$. This rigid foundation is moving horizontally at a constant speed of 20 m/s. At the top of the structure ($x_2 = 10 \text{ cm}$) the Dirichlet condition $u_1 = 0$, $u_2 = -2.5 \times 10^{-3} \text{ cm}$ ensures an initial compression of the body. The friction coefficient is larger than 1 ($\mathcal{F} = 1.2$). The other characteristics of the simulation are summarized in Table 1. Additionally, Problem (7) is approximated in time with the midpoint scheme (see Section 6 or (13) for more details). The C++ program that performs the tests is available along with Getfem++ (27).
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Structured mesh of the square $\{(x_1, x_2) \in ]0, 10 \text{ cm}[ \times ]0, 10 \text{ cm}]$ and an example of deformation (at $t = 0.01$ s). The rigid foundation lies at $x_2 = 0$ and has a constant horizontal velocity of 20 m/s. The structure is clamped on its top ($x_2 = 10 \text{ cm}$)</figcaption></figure>

<table><caption>**Table 1** Main characteristics of the simulation</caption><tr><td>Density ρ</td><td>10<sup>3</sup> kg/m<sup>3</sup></td></tr><tr><td>Domain</td><td>]0, 10 cm[×]0, 10 cm[</td></tr><tr><td>Lamé coefficients</td><td>λ = 300 MPa, μ = 150 MPa</td></tr><tr><td>Simulation time</td><td>0.02 s</td></tr><tr><td>Friction coefficient F</td><td>1.2</td></tr><tr><td>Horizontal velocity of the rigid foundation</td><td>20 m/s</td></tr></table>

Due to friction and the fact that the rigid foundation has a constant horizontal velocity, there is a source of energy in the system. The evolution of the total energy (given by (19)) for an element size $h = 0.5$ cm and various time steps is shown in Fig. 7 for four situations.

The first graph of Fig. 7 corresponds to the standard semi-discretization. It clearly shows the instability of the midpoint scheme applied to the standard semi-discretization. A rather unique feature, reserved to the discretization of dynamic contact problems, is that the smaller the time step is, the most the scheme is unstable. One can also remark that the scheme is reasonably stable in the first half of the simulation period. It is probably due to the fact that, at the beginning of the simulation, the body is pressed down to the foundation and is submitted to a relatively monotone loading due to the friction force. Consequently, the number of transitions between contact and noncontact is low.

The second and third graphs correspond to the mass redistribution applied only on the normal component of the displacement and on both components, respectively. The stability of the midpoint scheme is recovered.

An interesting situation is described on the fourth graph of this figure, which corresponds to the standard semi-discretization with an adaptation of the midpoint scheme where the contact forces have been implicated (see (5, 6)). Even though this scheme is proven to be stable (energy dissipative in fact) and this is also the case here, one can see that the result is rather different than the two previous graphs.
<figure><img src="image_8.png" /><figcaption>Fig. 7 Evolution of the energy for the following four situations: standard discretization, partial mass redistribution, whole-mass redistribution and standard discretization with implicit contact forces</figcaption></figure>

This is more clear in Fig. 8, where the graphs correspond to the difference between the total energy and the energy transferred (provided and dissipated) by friction. The latter is given at the time step $l\Delta t$ by

$$ \xi_l = \sum_{k=1}^{l} \sum_{i=1}^{n_c} \lambda_{\tau,i}^k (\dot{u}^k)^T \boldsymbol{\tau}_i \Delta t, $$

where $\lambda_{\tau,i}^k$ and $\dot{u}^k$ are the corresponding quantities to $\lambda_{\tau,i}(k\Delta t)$ and $\dot{u}(k\Delta t)$ at the $k$th time step of the midpoint scheme. These graphs reflect better the stability of each scheme. Of course, this does not change the conclusion for the midpoint scheme with the standard semi-discretization, which still appears to be unstable. However, from the two graphs corresponding to the mass redistribution method, we can see that for both cases, the energy conservation is obtained asymptotically for a time step going to zero. This is not the case for the fourth graph, corresponding to the standard semi-discretization with implicit contact forces. The scheme is stable, but it does not converge towards an energy conserving solution. This is due to the fact that a certain amount of energy is lost at each impact of each node, independently of the length of the time step. Since the number of
<figure><img src="image_9.png" /><figcaption>**Fig. 8** Evolution of the difference between the total energy and the energy transferred by friction</figcaption></figure>

impacts grows in this simulation for decreasing time steps, the smallest the time step is, the most dissipative the scheme is. In conclusion, we can say that this leads to a non-physical solution.

A convergence test on the time interval [0, 0.005 s] has been performed. The results are shown in Figs 9 and 10. A reference solution has been computed on a refined mesh ($h = 0.0625$ cm) for a small time step ($\Delta t = 5 \times 10^{-7}$ s) using the partial mass redistribution. Then the differences between this reference solution and four experiments whose characteristics are presented in Table 2 are computed. The curves in Fig. 9 present the maxima of the $L^2(\Omega)$-norm and $H^1(\Omega)$-semi-norm in [0, 0.005 s] for the four experiments. Numerical convergence with an order <1 is found. Due to the weak regularity of the solution, a faster rate of convergence cannot reasonably be expected. Note that a mathematical result of convergence of numerical solutions towards a solution of the continuous problem is an open problem. Until now, neither existence nor uniqueness results have been established on this model (unless a certain number of regularizations).

The convergence is also illustrated in Fig. 10, where the evolution of the density of friction force for different experiments is shown. Both the cases with a redistribution of mass only on the vertical component and on both components of the displacement are shown. In the two cases, the numerical solution seems to converge towards the same solution. This means that at least for the presented test case, the differentiated treatment of the two conditions is not strictly mandatory to obtain reliable
<figure><img src="image_10.png" /><figcaption>**Fig. 9** Convergence test for the partial mass redistribution. Maxima of the $L^2(\Omega)$-norm and $H^1(\Omega)$-semi-norm in $[0, 0.005 \text{ s}]$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10** Comparison of the density of friction force at the point $(0, 0)$ for the convergence test</figcaption></figure>

<table><caption>**Table 2** Element sizes and time steps for the convergence test</caption><thead><tr><th></th><th>Experiment 1</th><th>Experiment 2</th><th>Experiment 3</th><th>Experiment 4</th><th>Reference solution</th></tr></thead><tbody><tr><td>Element size <i>h</i> (cm)</td><td>2</td><td>1</td><td>0.5</td><td>0.25</td><td>0.0625</td></tr><tr><td>Time step Δ<i>t</i> (s)</td><td>1.6 × 10<sup>-5</sup></td><td>8 × 10<sup>-6</sup></td><td>4 × 10<sup>-6</sup></td><td>2 × 10<sup>-6</sup></td><td>5 × 10<sup>-7</sup></td></tr></tbody></table>
numerical results. Of course, with the redistribution of mass on both components, the well-posedness of the semi-discretization is not guaranteed (see the discussion in Section 6). A real difference may occur if there is a dynamical bifurcation. But the exhibition of such a dynamical bifurcation is also still an open problem.

## 8. Concluding remarks

We adapted the mass redistribution method for the elastodynamic contact problem with friction. The proposed strategy, which is to apply the mass redistribution only on the normal component corresponding to the contact condition, allows to transform the semi-discrete problem into a regular one-sided Lipschitz differential inclusion. The advantage is that any reasonable time integration scheme is then convergent (see (24)) at least for a fixed mesh. Moreover, the fully discrete problem is also well-posed for a sufficiently small time step. The simple example described in Section 6 shows that this is not the case when the mass redistribution is applied on both the contact and friction conditions.

The test case presented in Section 7 confirms that the stability is gained by the mass redistribution method. However, for this simple test case, there is no significant difference between the two strategies how to apply the mass redistribution. One may think that it makes no difference numerically.

For the moment, we do not have any dynamical bifurcation example for the Coulomb friction law with a constant friction coefficient that would test the difference. A perspective of this work would be to consider a coefficient of friction depending on the sliding velocity. In (17, 28), such a friction coefficient has been considered and multi-solutions are given in a one-dimensional case. It is proven that an additional mass on the contact boundary allows to recover the uniqueness of the solution. Moreover, it selects a particular solution which is related to the perfect delay criterion introduced in (29) for contact problems with friction. In this context, obtaining a well-posed semi-discrete problem would be more crucial and it would be interesting to see if the same solution is selected.

## Acknowledgements

T. Ligurský has been supported by the grant no. 18008 of the Charles University Grant Agency and under grant no. 201/07/0294 of the Grant Agency of the Czech Republic. He acknowledges also the support of the Nečas Center for Mathematical Modeling.

## References

1. L. Paoli, Time discretization of vibro-impact, *Philos. Trans. R. Soc. Lond. A* **359** (2001) 2405–2428.

2. J. J. Moreau, Liaisons unilatérales sans frottement et chocs inélastiques, *C. R. Acad. Sci. Ser. II* **296** (1983) 1473–1476.

3. J. J. Moreau, Numerical aspects of the sweeping process, *Comput. Methods Appl. Mech. Engng* **177** (1999) 329–349.

4. L. Paoli and M. Schatzman, Approximation et existence en vibro-impact, *C. R. Acad. Sci. Paris Sér. I* **329** (1999) 1103–1107.

5. J. O. Hallquist, K. Schweizerhof and D. Stillman, Efficiency refinements of contact strategies and algorithms in explicit finite element programming, *Proceedings of the Third
*International Conference on Computational Plasticity: Fundamentals and Applications*
(ed. D. R. J. Owen, E. Onate & E. Hinton; Pineridge Press, Swansea 1992) 457–481.

6. N. J. Carpenter, Lagrange constraints for transient finite element surface contact, *Int. J. Numer. Methods Engng* **32** (1991) 103–128.

7. T. A. Laursen and G. R. Love, Improved implicit integrators for transient impact problems- geometric admissibility within the conserving framework, *ibid.* **53** (2002) 245–274.

8. T. A. Laursen and V. Chawla, Design of energy conserving algorithms for frictionless dynamic contact problems, *ibid.* **40** (1997) 863–886.

9. P. Hauret and P. Le Tallec, Energy controlling time integration methods for nonlinear elastodynamics and low-velocity impact, *Comput. Methods Appl. Mech. Engng* **195** (2006) 4890–4916.

10. F. Armero and E. Petocz, Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems, *ibid.* **158** (1998) 269–300.

11. P. Hauret, Numerical methods for the dynamic analysis of two-scale incompressible nonlinear structures. Thèse de Doctorat (2004, Ecole Polytechnique, France).

12. T. J. R. Hughes, R. L. Taylor, J. L. Sackman, A. Curnier and W. Kanok-Nukulchai, A finite element method for a class of contact-impact problems, *Comput. Methods Appl. Mech. Engng.* **8** (1976) 249–276.

13. H. B. Khenous, P. Laborde and Y. Renard, Mass redistribution method for finite element contact problems in elastodynamics, *Eur. J. Mech. A Solids* **27** (2008) 918–932.

14. Y. Renard, The singular dynamic method for constrained second order hyperbolic equations: application to dynamic contact problems, *J. Comput. Appl. Math.* **234** (2010) 906–923.

15. P. Hauret, Mixed interpretation and extensions of the equivalent mass matrix approach for elastodynamics with contacts, *Comput. Methods Appl. Mech. Engng* **199** (2010) 2941–2957.

16. C. Hager, S. Hüeber and B. I. Wohlmuth, A stable energy conserving approach for frictional contact problems based on quadrature formulas, *Int. J. Numer. Methods Engng* **73** (2008) 205–225.

17. Y. Renard, Modélisation des instabilités liées au frottement sec des solides élastiques, aspects théoriques et numériques. Ph.D. Thesis (1998, LMC-IMAG Grenoble, France).

18. J.-C. Paumier and Y. Renard, Surface perturbation of an elastodynamic contact problem with friction, *Eur. J. Appl. Math.* **14** (2003) 465–483.

19. N. Kikuchi and J. T. Oden, *Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods* (SIAM, Philadelphia 1988).

20. H. B. Khenous, J. Pommier and Y. Renard, Hybrid discretization of the Signorini problem with Coulomb friction. Theoretical aspects and comparison of some numerical solvers, *Appl. Numer. Math.* **56** (2006) 163–192.

21. P. Laborde and Y. Renard, Fixed points strategies for elastostatic frictional contact problems, *Math. Methods Appl. Sci.* **31** (2008) 415–441.

22. I. Ekeland and R. Temam, *Analyse Convexe et Problèmes Variationnels* (Études Mathématiques, Dunod, Gauthier-Villars, Paris 1974).

23. K. Deimling, *Multivalued Differential Equations*, De Gruyter Studies in Nonlinear Analysis and Applications 1 (Walter de Gruyter, Berlin 1992).

24. A. Dontchev and F. Lempio, *Difference Methods for Differential Inclusions: A Survey*, SIAM Rev **34** (1992) 263–294.

25. P. Hild and Y. Renard, Local uniqueness and continuation of solutions for the discrete Coulomb friction problem in elastostatics, *Q. Appl. Math.* **63** (2005) 553–573.
26. T. Ligurský, Theoretical analysis of discrete contact problems with Coulomb friction, *Appl. Math*. Submitted (2009).

27. Y. Renard and J. Pommier, *Getfem++*, An open source generic C++ library for finite element methods, http://home.gna.org/getfem.

28. Y. Renard, Singular perturbation approach to an elastic dry friction problem with a non-monotone friction coefficient, *Q. Appl. Math. LVIII* (2000) 303–324.

29. I. R. Ionescu and J.-C. Paumier, On the contact problem with slip rate dependent friction in elastodynamic, *Eur. J. Mech. A Solids* **4** (1994) 555–568.