# The singular dynamic method for constrained second order hyperbolic equations: Application to dynamic contact problems

Yves Renard

► To cite this version:

Yves Renard. The singular dynamic method for constrained second order hyperbolic equations: Application to dynamic contact problems. Journal of Computational and Applied Mathematics, 2010, 234, pp.906-923. 10.1016/j.cam.2010.01.058. hal-01461799

HAL Id: hal-01461799

https://hal.science/hal-01461799

Submitted on 8 Feb 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
# The singular dynamic method for constrained second order hyperbolic equations. Application to dynamic contact problems

Yves RENARD<fnref n="1" />

## Abstract

The purpose of this paper is to present a new family of numerical methods for the approximation of second order hyperbolic partial differential equations submitted to a convex constraint on the solution. The main application is dynamic contact problems. The principle consists in the use of a singular mass matrix obtained by the mean of different discretizations of the solution and of its time derivative. We prove that the semi-discretized problem is well-posed and energy conserving. Numerical experiments show that this is a crucial property to build stable numerical schemes.

Keywords: hyperbolic partial differential equation, constrained equation, finite element methods, variational inequalities.

65M60, 35L87, 74M15, 35Q74.

## Introduction

An interesting class of hyperbolic partial differential equations with constraints on the solution consists in elastodynamic contact problems for which the vast majority of traditional numerical schemes show spurious oscillations on the contact displacement and stress (see for instance [12, 9, 10]). Moreover, these oscillations do not disappear when the time step decreases. Typically, they have instead tended to increase. This is a characteristic of order two hyperbolic equations with unilateral constraints that makes it very difficult to build stable numerical schemes. These difficulties have already led to many research under which a variety of solutions were proposed. Some of them consists in adding damping terms (see [24] for instance), but with a loss of accuracy on the solution, or to implicit the contact stress [7, 6] but with a loss of kinetic energy which could be independent of the discretization parameters (see the numerical experiments). Some energy conserving schemes have also been proposed in [11, 25, 17, 16, 9, 2, 10]. Unfortunately, these schemes, although more satisfactory than the most other schemes, lead to large oscillations on the contact stress. Besides, most of them do not strictly respect the constraint.

In this paper, we propose a new class of methods whose principle is to make different approximations of the solution and of its time derivative. Compared to the classical space semi-discretization, this corresponds to a singular modification of the mass matrix. In this sense, it is in the same class of methods than the mass redistribution method proposed in [12, 13] for elastodynamic contact problems. Indeed, in this latter method, the mass matrix has zero components for all the nodes on the contact boundary (which limits its application to constraints on a boundary). The singular modification of the mass matrix justify the proposed terminology “singular dynamic method”. The main feature is to provide a well-posed space semi-discretization. The numerical tests show that it has a crucial influence

<footnote n="1">¹Université de Lyon, CNRS, INSA-Lyon, ICJ UMR5208, F-69621, Villeurbanne, France. Yves.Renard@insa-lyon.fr</footnote>
on the stability of standard scheme and on the quality of the approximation, especially for the computation of Lagrange multipliers corresponding to the constraints.

The classical semi-discretizations, for example with finite element methods, give a problem in time which is a measure differential inclusion (see [19, 20, 21, 22]). Such a differential inclusion is systematically ill-posed, unless an additional impact law is considered. However, the scheme obtained with the addition of an impact law in [21] leads also to spurious oscillations.

The semi-discretization we propose here leads to a problem which is equivalent to a regular Lipschitz ordinary differential equation. Thus, time integration schemes at least converge for a fixed space discretization when the time step tends to zero. This work generalizes in a sense the methods presented in [13, 8].

The outline of the paper is the following. Section 1 is devoted to the description of the abstract hyperbolic equation with constraints and the equivalent variational inequality. Section 2 presents the new approximation methods and the main results of well-posedness and energy conservation. Then, in section 3, a non-trivial model problem which corresponds for instance to the dynamics of a thin membrane under an obstacle condition is developed. An example of well-posed discretization is also built in this section. Section 4 briefly describes the fully discrete problem obtained with the finite difference midpoint scheme and presents some numerical experiments on this model which shows in particular that the midpoint scheme is stable with well-posed semi-discretizations and unstable otherwise. Finally, in Section 5, the proposed method is applied to an elastodynamic contact problem.

# 1 The abstract hyperbolic equation

Let $\Omega \subset \mathbb{R}^d$ be a Lipschitz domain and $H = L^2(\Omega)$ the standard Hilbert space of square integrable functions on $\Omega$. Let $W$ be a Hilbert space such that

$$W \subset H \subset W',$$

with dense compact and continuous inclusions and let

$$A: W \to W'$$

be a linear self-adjoint elliptic continuous operator, i.e. which satisfies

$$\langle Aw, v \rangle_{W', W} = \langle Av, w \rangle_{W', W}, \quad \forall v, w \in W,$$

$$\exists \alpha > 0, \quad \forall w \in W, \quad \langle Aw, w \rangle_{W', W} \geq \alpha \|w\|_W^2, \quad \exists c > 0, \quad \forall w \in W, \quad \|Aw\|_{W'} \leq c \|w\|_W.$$

We consider the following problem

$$\left\{
\begin{array}{ll}
\text{Find } u: [0, T] \to K \text{ such that} & \\
\displaystyle\frac{\partial^2 u}{\partial t^2}(t) + Au(t) \in f - N_K(u(t)) & \text{for a.e. } t \in (0, T], \\
u(0) = u_0, \quad \displaystyle\frac{\partial u}{\partial t}(0) = v_0,
\end{array}
\right.
\qquad (1)$$

where $K$ is a closed convex nonempty subset of $W$, $f \in W'$, $u_0 \in K$, $v_0 \in H$, $T > 0$ and $N_K(u)$ is the normal cone to $K$ defined by (see for instance [5] for a detailed presentation of differential inclusions)

$$N_K(u) = \begin{cases} \emptyset & \text{if } u \notin K, \\ \{f \in W' : \langle f, w - u \rangle_{W', W} \le 0, \forall w \in K\} & \text{if } u \in K. \end{cases}$$
This means that $u(t)$ satisfies the second order hyperbolic equation and is constrained to remain in the convex $K$. As far as we know, there is no general result of existence and uniqueness for the solution to this kind of equation. Some existence results for a scalar Signorini problem can be found in [18, 15]. Introducing now the linear and bilinear symmetric maps

$$l(v) = \langle f, v \rangle_{W', W}, \quad a(u, v) = \langle Au, v \rangle_{W', W}$$

Problem (1) can be rewritten as the following variational inequality:

$$\left\{
\begin{array}{l}
\text{Find } u : [0, T] \to K \text{ such that for a.e. } t \in (0, T], \\
\langle \frac{\partial^2 u}{\partial t^2}(t), w - u(t) \rangle_{W', W} + a(u(t), w - u(t)) \ge l(w - u(t)) \quad \forall w \in K, \\
u(0) = u_0, \quad \frac{\partial u}{\partial t}(0) = v_0.
\end{array}
\right.
\qquad (2)$$

Note that the terminology “variational inequality” is used here in the sense that Problem (1) derives from the conservation of the energy functional

$$J(t) = \frac{1}{2} \int_{\Omega} \left( \frac{\partial u}{\partial t}(t) \right)^2 dx + \frac{1}{2} a(u(t), u(t)) - l(u(t)) + I_K(u(t)),$$

where $I_K(u(t))$ is the convex indicator function of $K$. However, it is generally not possible to prove that each solution to Problem (2) is energy conserving, due to the weak regularity involved.

## 2 Approximation and well-posedness result

The goal of this section is to present well-posed space semi-discretizations of Problem (2). The strategy adopted is to use a Galerkin method with different approximations of $u$ and of $v = \frac{\partial u}{\partial t}$. Let $W^h$ and $H^h$ be two finite dimensional vector subspaces of $W$ and $H$ respectively. Let $K^h \subset W^h$ be a closed convex nonempty approximation of $K$. The proposed approximation of Problem (2) is the following mixed approximation:

$$\left\{
\begin{array}{l}
\text{Find } u^h : [0, T] \to K^h \text{ and } v^h : [0, T] \to H^h \text{ such that} \\
\displaystyle \int_{\Omega} \frac{\partial v^h}{\partial t} (w^h - u^h) dx + a(u^h, w^h - u^h) \ge l(w^h - u^h) \quad \forall w^h \in K^h, \quad \forall t \in (0, T], \\
\displaystyle \int_{\Omega} (v^h - \frac{\partial u^h}{\partial t}) q^h dx = 0 \quad \forall q^h \in H^h, \quad \forall t \in (0, T], \\
u^h(0) = u_0^h, \quad v^h(0) = v_0^h,
\end{array}
\right.
\qquad (3)$$

where $u_0^h \in K^h$ and $v_0^h \in H^h$ are some approximations of $u_0$ and $v_0$ respectively. Of course, when $H^h = W^h$ this corresponds to a standard Galerkin approximation of Problem (2).

Let $\varphi_i$, $1 \le i \le N_W$ and $\psi_i$, $1 \le i \le N_H$ be some basis of $W^h$ and $H^h$ respectively, and let the matrices $A, B$ and $C$, of sizes $N_W \times N_W, N_H \times N_W$ and $N_H \times N_H$ respectively, and the vectors $L, U$ and $V$, of size $N_W, N_W$ and $N_H$ respectively, be defined by

$$A_{i,j} = a(\varphi_i, \varphi_j), \quad B_{i,j} = \int_{\Omega} \psi_i \varphi_j dx, \quad C_{i,j} = \int_{\Omega} \psi_i \psi_j dx,$$

$$L_i = l(\varphi_i), \quad u^h = \sum_{i=1}^{N_W} U_i \varphi_i, \quad v^h = \sum_{i=1}^{N_H} V_i \psi_i.$$
Then, the algebraic expression of Problem (3) is the following:

$$ egin{cases} 
\text{Find } U : [0, T] \to \bar{K}^h \text{ and } V : [0, T] \to \mathbb{R}^{NH} \text{ such that } \forall t \in (0, T], \\ 
(W - U(t))^T (B^T \dot{V}(t) + AU(t)) \ge (W - U(t))^T L, \quad \forall W \in \bar{K}^h, \\ 
CV(t) = B\dot{U}(t), \\ 
U(0) = U_0, \quad V(0) = V_0. 
\end{cases} 
\quad (4) $$

where $\dot{U}(t)$ and $\dot{V}(t)$ denote the derivative with respect to $t$ of $U(t)$ and $V(t)$ respectively and $\bar{K}^h$ is defined by

$$ \bar{K}^h = \{W \in \mathbb{R}^{Nw} : \sum_{i=1}^{Nw} W_i \varphi_i \in K^h\}. $$

At this stage, the unknown $V$ can be eliminated since $C$ is always invertible which leads to the relation $V(t) = C^{-1}B\dot{U}(t)$. Thus denoting

$$ M = B^T C^{-1} B, $$

Problem (4) can be rewritten

$$ egin{cases} 
\text{Find } U : [0, T] \to \bar{K}^h \text{ such that} \\ 
(W - U(t))^T (M\ddot{U}(t) + AU(t)) \ge (W - U(t))^T L, \quad \forall W \in \bar{K}^h, \forall t \in (0, T], \\ 
U(0) = U_0, \quad B\dot{U}(0) = CV_0. 
\end{cases} 
\quad (5) $$

**Remark 1** If the couple of discretization spaces $(H^h, W^h)$ satisfies a classical inf-sup condition then the matrix $B$ is surjective and the initial condition $B\dot{U}(0) = CV_0$ is always admissible. Conversely, if $B$ is not surjective then the initial condition $V_0$ has to satisfy the following condition

$$ V_0 \in \text{Im}(C^{-1}B). \quad (6) $$

This condition is also implicitly contained in Problem (4).

In comparison with the standard approximation where $H^h = W^h$ the only difference introduced by the presented method is replacing the standard mass matrix $\left(\int_{\Omega} \varphi_i \varphi_j dx\right)_{i,j}$ by $M = B^T C^{-1} B$. In the interesting cases where $\dim(H^h) < \dim(W^h)$ this corresponds to replace the standard invertible mass matrix by a singular one. We propose to call this kind of method a singular dynamic method. Of course, the numerical implementation will be facilitated when the matrix $C$ is diagonal. This is the case for instance when $H^h$ is defined with $P_0$ finite element method or with a more general finite element method using an adapted sub-integration (lumped mass matrix). We will see how, rather surprisingly, the introduction of a singular mass matrix allows to recover the well-posedness of the approximation.

The goal now is to give a sufficient condition for Problem (5) (or equivalently Problem (3) or (4)) to be well-posed. To this end, we will define a more restrictive framework (see the concluding remarks for a possibility to extend this framework). We will suppose that $K^h$ is defined by

$$ K^h = \{w^h \in W^h : g^i(w^h) \le \alpha^i, 1 \le i \le N_g\}, $$
where $\alpha^i \in \mathbb{R}$ and $g^i : W^h \to \mathbb{R}$, $1 \le i \le N_g$ are some linearly independent linear maps. Of course, this restricts the possibilities concerning the convex $K$ since $K^h$ is supposed to be an approximation of $K$. With vector notations this leads to

$$ \bar{K}^h = \{W \in \mathbb{R}^{N_w} : (G^i)^T W \le \alpha_i, 1 \le i \le N_g\}, $$

where $G^i \in \mathbb{R}^{N_w}$ are such that $g^i(w^h) = (G^i)^T W$, $1 \le i \le N_g$. We will also denote by $G$ the $N_W \times N_g$ matrix whose components are

$$ G_{ij} = (G^i)_j. $$

Let us consider the subspace $F^h$ of $W^h$ defined by

$$ F^h = \{w^h \in W^h : \int_{\Omega} w^h q^h = 0 \quad \forall q^h \in H^h\}. $$

Then, the corresponding set $F = \{W \in \mathbb{R}^{N_w} : \sum_{i=1}^{N_W} W_i \varphi_i \in F^h\}$, is such that

$$ F = \text{Ker}(B). $$

In this framework, we will prove that the following condition is sufficient for the well-posedness of the discrete problem (5):

$$ \inf_{\substack{Q \in \mathbb{R}^{N_g} \\ Q \neq 0}} \sup_{\substack{W \in F \\ W \neq 0}} \frac{Q^T GW}{|Q| |W|} > 0, \quad (7) $$

where $|Q|$ and $|W|$ stands for the Euclidean norm of $Q$ in $\mathbb{R}^{N_g}$ and $W$ in $\mathbb{R}^{N_w}$ respectively. This condition is equivalent to the fact that the linear maps $g^i$ are independent on $F^h$ and also to the fact that $G$ is surjective on $F$. A direct consequence is that it implies $\dim(F^h) \ge N_g$, and consequently

$$ \dim(H^h) \le \dim(W^h) - N_g. $$

This again prescribed some conditions on the approximation made which links $W^h$, $H^h$ and also $K^h$. We will see in Section 3 that this condition can be satisfied for interesting practical situations. We can now prove the following result:

**Theorem 1** If $W^h$, $H^h$ and $K^h$ satisfy the condition (7) then Problem (5) admits a unique solution. Moreover, this solution is Lipschitz-continuous with respect to $t$.

First, let us establish the following intermediary result:

**Lemma 1** If $W^h$, $H^h$ and $K^h$ satisfy the condition (7) then there exists $F^c$ a sub-space of $\mathbb{R}^{N_w}$ such that $F^c \subset \text{Ker}(G)$ and such that $F$ and $F^c$ are complementary sub-spaces.

**Proof.** For $W \in \mathbb{R}^{N_w}$ let $X_F \in F$ be such that

$$ G(X_F) = G(W). $$

Such an $X_F$ exists since a consequence of condition (7) is that the matrix $G$ defines a surjective linear map from $F$ to $\mathbb{R}^{N_g}$. Thus

$$ W = (W - X_F) + X_F, $$
is a decomposition of W with $W - X_F \in \text{Ker}(G)$ and $X_F \in F$. This proves that $\mathbb{R}^{N_W} = F + \text{Ker}(G)$. The result of the lemma is then a consequence of the basis extension theorem. $\square$

**Proof of Theorem 1.** Now, using the result of Lemma 1, let us decompose $U, W \in \mathbb{R}^{N_W}$ as

$$U = U_F + U_{Fc}, \quad W = W_F + W_{Fc},$$

with $U_F, W_F \in F$ and $U_{Fc}, W_{Fc} \in F^c$. The inequality of (5) can be written for all $t \in (0, T]$

$$\begin{aligned} & (W_{Fc} - U_{Fc})^T (M \ddot{U}_{Fc} + AU_{Fc} + AU_F) + (W_F - U_F)^T (AU_{Fc} + AU_F) \\ & \geq (W_{Fc} - U_{Fc})^T L + (W_F - U_F)^T L, \quad \forall W_F \in \bar{K}^h \cap F, \forall W_{Fc} \in F^c. \end{aligned} \quad (8)$$

Taking now $W_{Fc} = U_{Fc}$ one obtains

$$(W_F - U_F)^T A U_F \geq (W_F - U_F)^T (L - A U_{Fc}), \quad \forall W_F \in \bar{K}^h \cap F. \quad (9)$$

This is a variational inequality for the unknown $U_F$. The solution to this variational inequality minimizes the quadratic functional

$$J_F(W_F) = \frac{1}{2} W_F^T A W_F - W_F^T (L - A U_{Fc})$$

over the closed convex $\bar{K}^h \cap F$. The ellipticity assumption implies that the matrix A is coercive. This leads to the existence and uniqueness of the solution $U_F$ to this variational inequality due to the Stampacchia theorem. Moreover, this is a classical result that $U_F$ depends Lipschitz-continuously on $U_{Fc}$. Indeed, let $U_F^1$ and $U_F^2$ be two solutions for $U_{Fc}^1$ and $U_{Fc}^2$ respectively. Then it can be straightforwardly deduced from the variational inequality that

$$(U_F^2 - U_F^1)^T A (U_F^2 - U_F^1) \leq (U_F^2 - U_F^1)^T A (U_{Fc}^2 - U_{Fc}^1).$$

Thus due to the coercivity of the matrix A one obtains for $c > 0$ a generic constant

$$|U_F^2 - U_F^1| \leq c|U_{Fc}^2 - U_{Fc}^1|.$$

We will thus use the notation $U_F(U_{Fc})$. Now, since inequality (8) has to be satisfied for all $W_{Fc} \in F^c$, this implies that $U_{Fc}(t)$ verifies for all $t \in (0, T]$

$$W_{Fc}^T M \ddot{U}_{Fc} = W_{Fc}^T (L - A U_{Fc} - A U_F (U_{Fc})) \quad \forall W_{Fc} \in F^c, \quad (10)$$

which represents an ordinary differential equation with Lipschitz-continuous right-hand side. Since the matrix M is nonsingular on $F^c$ (because $F^c$ is complementary to $F = \text{Ker}(B)$ and $M = B^T C^{-1} B$) there exists a unique solution to the associated initial value problem with the initial conditions $U_{Fc}(0) = (U_0)_{Fc}$ and $B \dot{U}_{Fc}(0) = CV_0$ assuming condition (6) when B is not surjective.

Since $U_{Fc}(t)$ is the solution to a second order autonomous ordinary differential equation with Lipschitz-continuous right-hand side, it has at least the regularity $U_{Fc} \in W^{3,\infty}(0, T; F^c)$. Finally, the whole $U(t)$ is Lipschitz-continuous with respect to $t$ due to the fact that $U_F$ depends Lipschitz-continuously on $U_{Fc}$. $\square$

Now, an interesting property is that the solution to Problem (5) satisfies the so-called persistency condition (see [16, 17]). This is a condition between $\dot{U}(t)$ and the Lagrange multipliers corresponding to the constraints. In a sense, this is a stronger condition than the
complementary condition which links $U(t)$ and the Lagrange multipliers. In fact, Problem (5) can be re-written

$$ egin{cases} 
\text{Find } U : [0, T] \to \bar{K}^h \text{ such that} \\ 
M\ddot{U}(t) + AU(t) \in L - N_{\bar{K}^h}(U(t)) \quad \forall t \in (0, T], \\ 
U(0) = U_0, \quad B\dot{U}(0) = CV_0, 
\end{cases} 
\quad (11) $$

where $N_{\bar{K}^h}(U(t))$ is the normal cone to $\bar{K}^h$. A straightforward computation leads to the following result:

$$ N_{\bar{K}^h}(U(t)) = \left\{ \begin{array}{ll} 
\emptyset & \text{if } U(t) \notin \bar{K}^h, \\ 
\left\{ \sum_{\substack{1 \le i \le N_g \\ (G^i)^T U(t) < \alpha_i}} \mu^i G^i : \mu_i \ge 0 \right\} & \text{if } U(t) \in \bar{K}^h. 
\end{array} \right. $$

Thus, introducing Lagrange multipliers, the discrete problem is also equivalent to the following one:

$$ egin{cases} 
\text{Find } U : [0, T] \to \bar{K}^h \text{ and } \lambda^i : [0, T] \to \mathbb{R}, 1 \le i \le N_g \text{ such that } \forall t \in (0, T] \\ 
M\ddot{U}(t) + AU(t) = L + \sum_{i=1}^{N_g} \lambda^i(t)G^i, \\ 
\lambda^i(t) \le 0, \quad (G^i)^T U(t) - \alpha_i \le 0, \quad \lambda^i(t)((G^i)^T U(t) - \alpha_i) = 0, \quad 1 \le i \le N_g, \\ 
U(0) = U_0, \quad B\dot{U}(0) = CV_0, 
\end{cases} 
\quad (12) $$

**Proposition 1** If $W^h$, $H^h$ and $K^h$ satisfy the condition (7) then the solution $U(t)$ to Problem (5) verifies the following persistency condition

$$ \lambda^i(t)(G^i)^T \dot{U}(t) = 0 \quad \forall t \in (0, T], \quad 1 \le i \le N_g. $$

**Proof.** With still the same decomposition as the one in Theorem 1 we deduce from (12) that $\lambda^i$, $1 \le i \le N_g$ satisfy

$$ W_F^T A U_F = W_F^T (L - A U_{Fc} + \sum_{i=1}^{N_g} \lambda^i G^i), \quad \forall W_F \in F. $$

Since $U_F$ depends Lipschitz-continuously on $U_{Fc}$, this equation implies that each $\lambda^i$ depends also Lipschitz-continuously on $U_{Fc}$. Thus each $\lambda^i(t)$ is Lipschitz-continuous with respect to $t$. But

$$ \lambda^i = 0 \text{ on } \text{Supp}((G^i)^T U - \alpha_i) = \omega^i \subset [0, T], \quad 1 \le i \le N_g, $$

where $\text{Supp}(f)$ denotes the support of the function $f(t)$. The continuity of $\lambda^i(t)$ implies

$$ \lambda^i = 0 \text{ on } \bar{\omega}^i. $$

Since $(G^i)^T U - \alpha_i = 0$ on the complement of $\omega_i$, then its derivative $(G^i)^T \dot{U}$ vanishes also on the interior of complementary of $\omega_i$ which proves the result of the proposition. $\square$

Now, we can prove that the persistency condition implies the energy conservation.
**Theorem 2** If $W^h$, $H^h$ and $K^h$ satisfy the condition (7) then the solution $U(t)$ to Problem (5) is energy conserving in the sense that the discrete energy

$$J^h(t) = \frac{1}{2} \dot{U}^T(t) M \dot{U}(t) + \frac{1}{2} U^T(t) A U(t) - U^T(t) L,$$

is constant with respect to $t$.

**Proof.** The first equation of (12) implies

$$\dot{U}^T M \ddot{U}(t) + \dot{U}^T A U(t) = \dot{U}^T L + \sum_{i=1}^{N_g} \dot{U}^T \lambda^i(t) G^i, \quad \text{on } [0, T].$$

Integrating from 0 to $t$ and using Proposition 1 one can conclude that

$$\frac{1}{2} \dot{U}^T(t) M \dot{U}(t) + \frac{1}{2} U^T(t) A U(t) - U^T(t) L = J^h(0).$$

□

## 3 A model problem

The goal of this section is to provide a simple but interesting situation for which some consistent approximations satisfy the condition (7). With $W = H^1(\Omega)$ and $K = \{w \in W : w \ge 0 \text{ a.e. on } \Omega\}$ we consider the following problem:

$$\left\{
\begin{array}{ll}
\text{Find } u : [0, T] \to K \text{ such that} & \\
\frac{\partial^2 u}{\partial t^2}(t) - \Delta u(t) \in f - N_K(u(t)) \quad \text{in } \Omega, \quad \text{for a.e. } t \in (0, T], & \\
\frac{\partial u}{\partial n} = 0 \quad \text{on } \Gamma_N, & \\
u = 0 \quad \text{on } \Gamma_D, & \\
u(0) = u_0, \quad \frac{\partial u}{\partial t}(0) = v_0, &
\end{array}
\right.$$

where $\Gamma_N$ and $\Gamma_D$ is a partition of $\partial\Omega$, $\Gamma_D$ being of non zero measure in $\partial\Omega$. This models for instance the contact between an antiplane elastic structure with a rigid foundation or a stretched drum membrane under an obstacle condition. In this situation, the mass redistribution method presented in [13] is not usable since the area subjected to potential contact is the whole domain. Consequently, this method would lead to suppress the mass on the whole domain which is a non consistent drastic change of the problem.

We build now the approximation spaces thanks to finite element methods. Let $\mathcal{T}^h$ a regular triangular mesh of $\Omega$ (in the sense of Ciarlet [4], $h$ being the diameter of the largest element) and $W^h$ be the $P_1$+ finite element space

$$W^h = \{w^h \in \mathcal{C}^0(\Omega) : w^h = \sum_{a_i \in \mathcal{A}} w_i \varphi_i + \sum_{T \in \mathcal{T}^h} w_T \varphi_T\},$$

where $\mathcal{A}$ is the set of the vertices of the mesh which do not lie on $\Gamma_D$, $\varphi_i$, $i \in \mathcal{A}$ are the piecewise affine function satisfying

$$\varphi_i(a_j) = \delta_{ij},$$
i.e. the shape functions of a $P_1$ Lagrange finite element method on $\mathcal{T}^h$. Each functions $\varphi_T$, $T \in \mathcal{T}^h$ is a cubic bubble function whose support is $T$. Let $H^h$ be the $P_0$ finite element space

$$H^h = \{v^h \in L^2(\Omega) : v^h = \sum_{T \in \mathcal{T}^h} v_T \mathbb{I}_T\},$$

and finally, let $K^h$ be defined as

$$K^h = \{w^h \in W^h : w^h(a_i) \ge 0 \text{ for all } a_i \text{ vertex of } \mathcal{T}^h\}, \quad (13)$$

which means that the constraints are only prescribed at the vertices of the mesh.

**Lemma 2** *This choice of $W^h$, $H^h$ and $K^h$ satisfies condition (7).*

**Proof.** The computation of $F^h$ gives

$$\begin{aligned} F^h &= \{w^h \in W^h : \int_T w^h dx = 0 \quad \forall T \in \mathcal{T}^h\} \\ &= \{w^h = \sum_{a_i \in \mathcal{A}} w_i \varphi_i + \sum_{T \in \mathcal{T}^h} w_T \varphi_T : w_T = - \int_T \sum_{a_i \in \mathcal{A}} w_i \varphi_i dx\}, \end{aligned}$$

while the functions $g^i$ are defined by

$$g^i(w^h) = -w^h(a_i), \quad a_i \in \mathcal{A}.$$

Thus one has $g^i(w^h) = -w_i$ and the surjectivity of $G$ on $F$ is obvious. $\square$

<figure><img src="image_2.png" /><figcaption>Figure 1: $h = 0.05$.</figcaption></figure>
## 4 Numerical experiments

<figure><img src="image_3.png" /><figcaption>Figure 2: Evolution of the energy, the displacement at the center point (0.5, 0.5) and the contact stress at the center point for a $P_1$/+P$_0$ method, a midpoint scheme and with $h = 0.1$.</figcaption></figure>

We present now some numerical experiments done on the problem described in the previous section, with

$$ \Omega = (0,1) \times (0,1), \quad \Gamma_D = \partial\Omega, \quad \Gamma_N = \emptyset, \quad f = -0.6. $$

The initial condition is

$$ u(0,x) = 0.02, \quad \dot{u}(0,x) = 0, \quad x \in \Omega, $$

and we consider a non-homogeneous Dirichlet condition

$$ u(t,x) = 0.02, \quad x \in \partial\Omega. $$

The structured mesh used can be viewed on Figure 1 where the solution is represented during the first impact on the obstacle. The numerical experiments are performed with our finite element library Getfem++ [23] (the program itself is available on Getfem++ web site). A semi-smooth Newton method is used to solve the discrete problem (see [1, 3, 14]). All the numerical experiments use the same definition of $K^h$ done by (13).

We mainly use a midpoint scheme for the time discretization of the problem. The midpoint scheme is interesting since it is energy conserving on the linear part (equation without constraint) but of course any other reasonably stable scheme can be applied. A midpoint scheme on Problem (5) has the following expression ($dt$ is the time step):
$$
\left\{
\begin{array}{l}
U^0 \text{ and } \bar{V}^0 \text{ be given. For } n \ge 0, \text{ find } U^{n+\frac{1}{2}} \in \bar{K}^h \\
(W - U^{n+\frac{1}{2}})^T (MZ^{n+\frac{1}{2}} + AU^{n+\frac{1}{2}}) \ge (W - U^{n+\frac{1}{2}})^T L, \quad \forall W \in \bar{K}^h, \\
U^{n+\frac{1}{2}} = \frac{U^n + U^{n+1}}{2}, \quad \bar{V}^{n+\frac{1}{2}} = \frac{\bar{V}^n + \bar{V}^{n+1}}{2}, \\
U^{n+1} = U^n + dt\bar{V}^{n+\frac{1}{2}}, \quad \bar{V}^{n+1} = \bar{V}^n + dtZ^{n+\frac{1}{2}},
\end{array}
\right.
\qquad (14)
$$

where $\bar{V}^n$ is here an approximation of $\dot{U}(ndt)$ and does not refer to the velocity $V(t)$ of Problem (4). Note that another equivalent possibility would be to apply a midpoint scheme directly to Problem (4) in order for instance to avoid the explicit computation of $M = B^T C^{-1} B$.

In the system (14) for $U^n$ and $V^n$ be given, $U^{n+\frac{1}{2}}$ is the solution to the inequation

$$
U^{n+\frac{1}{2}} \in \bar{K}^h, \quad (W - U^{n+\frac{1}{2}})^T \left( \frac{4}{dt^2} M U^{n+\frac{1}{2}} + A U^{n+\frac{1}{2}} \right) \geq (W - U^{n+\frac{1}{2}}) \tilde{L}, \quad \forall W \in \bar{K}^h,
$$

where $\tilde{L}$ depends on $U^n$ and $V^n$. Due to the coercivity of the matrix $A$, this variational inequality admits a unique solution, whatever the choice of $W^h$, $H^h$ and $K^h$ (even with a standard discretization, i.e. in the case $W^h = H^h$). Of course, this well-posedness of Problem (14) augurs nothing of the stability of the whole scheme.

The first numerical test is made with the midpoint scheme and the approximation presented in Section 3, that is a $P_1+/P_0$ method ($P_1+$ for $u^h$ and $P_0$ for $\dot{u}^h$).

In good accordance with the theoretical results, the curves on Figure 2 show that the energy tends to be conserved when the time step decreases (an experiment with $dt = 10^{-4}$ has been performed but the difference with the one for $dt = 10^{-3}$ is not visible). Moreover, both the displacement and the contact stress taken at the point $(0.5, 0.5)$ are smooth and converge satisfactorily when the time step diminish.

<figure><img src="image_4.png" /><figcaption>Figure 3: Evolution of the energy, the displacement at the center point (0.5, 0.5) and the contact stress at the center point for a $P_1/P_0$ method, a midpoint scheme and with $h = 0.1$.</figcaption></figure>
Conversely, the curves on Figures 3 and 4 obtained for a $P_1/P_0$ method and a $P_1/P_1$ method respectively are unstable. The energy is growing very fast after the first impact. The displacement and the contact stress are very oscillating and do not converge. Moreover, the instabilities are more important for the smallest time step. These two methods do not satisfy the condition (7) since $\dim(H^h) \ge \dim(W^h)$. Note that the $P_1/P_1$ method corresponds to a standard Galerkin approximation of the problem.

Even thought we do not have a proof that condition (7) is satisfied for $P_1+/P_1$ and $P_2/P_1$ methods (still with the same $K^h$), Figures 5 and 6 show that the midpoint scheme is stable and converging for these two methods (here also, some experiments with $dt = 10^{-4}$ have been performed with no visible differences with the one for $dt = 10^{-3}$).

An interesting situation is also presented in Figures 7, 8 and 9 where a backward Euler scheme is used. This time integration scheme is unconditionally stable because it is possible to prove that the discrete energy decreases from an iteration to another (see [12] for instance). This is the case for any choice of $W^h$ and $H^h$. Consequently, this method presents some smooth results for the displacement and the contact stress. However, the energy decreases rapidly for large time steps. Figure 7 shows that for a well-posed method, the energy tends to be conserved for small time steps, but Figures 8 and 9 show that with an ill-posed method (such as classical discretizations) there is an energy loss at the impact which do not vanish when the time step and the mesh size decrease. This means that with an ill-posed method, we do not approximate a physical solution of the problem whenever one expects energy conservation to be satisfied at the limit.

<figure><img src="image_5.png" /><figcaption>Figure 4: Evolution of the energy, the displacement at the center point (0.5, 0.5) and the contact stress at the center point for a $P_1/P_1$ method, a midpoint scheme and with $h = 0.1$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Evolution of the energy, the displacement at the center point (0.5, 0.5) and the contact stress at the center point for a $P_1$/ $P_1$ method, a midpoint scheme and with $h = 0.1$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Evolution of the energy, the displacement at the center point (0.5, 0.5) and the contact stress at the center point for a $P_2$/ $P_1$ method, a midpoint scheme and with $h = 0.1$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Evolution of the energy, the displacement at the center point (0.5, 0.5) and the contact stress at the center point for a $P_1$/$P_0$ method, a backward Euler scheme and with $h = 0.1$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Evolution of the energy, the displacement at the center point (0.5, 0.5) and the contact stress at the center point for a $P_1$/$P_0$ method, a backward Euler scheme and with $h = 0.1$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: *Evolution of the energy for a $P_1/P_0$ method, a backward Euler scheme and $dt = 0.001$ for different values of the mesh size.*</figcaption></figure>

A convergence test for both $h$ and $dt$ decreasing is shown in Fig. 10. The midpoint scheme is used for the time integration. The curves represent the displacement at the center point (the curve are vertically shifted for better visibility) for decreasing values of the mesh size and the time step. The ratio $dt/h$ is kept constant. The test on the finest mesh ($dt = 3.125 \times 10^{-4}$ and $h = 6.25 \times 10^{-2}$) gives some similar results for the two methods ($P_1 + /P_0$ and $P_2/P_1$). The convergence may seem rather slow. However, the evolution of the displacement is not smooth due to impacts. This of course limits the convergence rate of the scheme (see also the test for linear elasticity in section 5). Moreover, this is only a representation of the displacement on a single point.

<figure><img src="image_11.png" /><figcaption>Figure 10: *Evolution of the displacement at the center point (0.5, 0.5) for both $h$ and $dt$ decreasing. The curves are vertically shifted for better visibility.*</figcaption></figure>

The evolution of the stability of the method for decreasing mesh sizes is illustrated by Fig. 11. The error plotted is the relative difference (in %) between the energy at the final time and the initial energy of the system. The midpoint scheme is still used. Each curve in this figure corresponds to a fixed mesh size and a decreasing time step $dt$. Once again, it can be noted that for a fixed space discretization, the energy tends to be conserved when the time step decreases for both $P_1 + /P_0$ and $P_2/P_1$ methods. For better readability of results, the abscissa is $dt/h$. It can be seen that for a constant $dt/h$, there is not a degradation of
the stability when the discretization parameters decrease.

<figure><img src="image_12.png" /><figcaption>Figure 11: Relative energy error on the final time with respect to the ratio dt/h for different mesh sizes.</figcaption></figure>

## 5 Extension to an elastodynamic contact problem

In this section, we consider the dynamic evolution of an elastic solid $\Omega$ submitted to a frictionless contact condition with a flat rigid foundation on a part of its boundary $\Gamma_C$ (see figure 12). Note that the consideration of a friction condition is not an additional difficulty concerning the stability of the numerical scheme because, for instance, the midpoint scheme applied to the friction term is stable (see for instance [12]). On the rest of the boundary, a Dirichlet condition is prescribed on $\Gamma_D$ and a Neumann one on $\Gamma_N$. With a linearized elasticity constitutive law, the problem reads as follows:

<figure><img src="image_13.png" /><figcaption>Figure 12: Elastic solid in contact with a flat rigid foundation</figcaption></figure>

$$ \rho\ddot{u} - \operatorname{div} \boldsymbol{\sigma}(u) = f, \quad \text{in } (0, T] \times \Omega, \quad (15) $$

$$ \boldsymbol{\sigma}(u) = \mathcal{A}\boldsymbol{\varepsilon}(u), \quad \text{in } (0, T] \times \Omega, \quad (16) $$

$$ \boldsymbol{\sigma}(u)\mathbf{n} = g, \quad \text{on } (0, T] \times \Gamma_N, \quad (17) $$

$$ u = 0, \quad \text{on } (0, T] \times \Gamma_D, \quad (18) $$

$$ u_N \le 0, \quad \boldsymbol{\sigma}_N(u) \le 0, \quad u_N \boldsymbol{\sigma}_N(u) = 0, \quad \boldsymbol{\sigma}_T(u) = 0, \quad \text{on } (0, T] \times \Gamma_C, \quad (19) $$

$$ u(0, x) = u_0(x), \quad \dot{u}(0, x) = u_1(x), \quad \text{in } \Omega, \quad (20) $$

where $u$ is the displacement, $\boldsymbol{\sigma}$ is the stress tensor, $\mathcal{A}$ is the fourth order tensor of elasticity coefficients, $u_N = u.n$ is the normal displacement on $\Gamma_C$, $\boldsymbol{\sigma}_N = (\boldsymbol{\sigma}n).n$ is the normal stress
(contact pressure) and $\sigma_T = \sigma_n - \sigma_N$ is the tangential stress on $\Gamma_C$. The weak form of this problem can be written using a Lagrange multiplier $\lambda_N$ being the contact stress as follows (hybrid formulation of the contact condition):

$$ egin{cases} 
\text{Find } u : [0, T] \to V \text{ and } \lambda_N : [0, T] \to X'_N \text{ satisfying} \\ 
\langle \rho \ddot{u}, v \rangle_{V',V} + a(u, v) = l(v) + \langle \lambda_N, v_N \rangle_{X'_N, X_N} \quad \forall v \in V, \\ 
\lambda_N \in \Lambda_N, \quad \langle \mu_N - \lambda_N, u_N \rangle_{X'_N, X_N} \ge 0, \quad \forall \mu_N \in \Lambda_N, \\ 
u(0) = u_0, \quad \dot{u}(0) = u_1. 
\end{cases} $$

where

$$ V = \{v \in H^1(\Omega; \mathbb{R}^3) : v = 0 \text{ sur } \Gamma_D\}, $$

$$ a(u, v) = \int_{\Omega} \mathcal{A}\varepsilon(u) : \varepsilon(v), \quad l(v) = \int_{\Omega} f.vdx + \int_{\Gamma_N} g.vd\Gamma, $$

$$ X_N = \{v_N|_{\Gamma_C} : v \in V\}, \quad \Lambda_N = \{\mu_N \in X'_N : \langle \mu_N, v_N \rangle_{X'_N, X_N} \ge 0, \forall v_N \in X_N, v_N \le 0\}. $$

Note that a primal formulation of contact could also be used as in Section 2. The goal here is to present the application of the singular dynamic method to the hybrid formulation of contact. The mixed approximation in displacement and velocity can be written as follows:

$$ egin{cases} 
\text{Find } u^h : [0, T] \to V^h, v^h : [0, T] \to H^h \text{ and } \lambda_N^h : [0, T] \to X_N^{h} \text{ satisfying for all } t \in (0, T] \\ 
\int_{\Omega} \rho \frac{\partial v^h}{\partial t} \cdot w^h dx + a(u^h, w^h) = l(w^h) + \langle \lambda_N^h, w_N^h \rangle_{X'_N, X_N} \quad \forall w^h \in V^h, \\ 
\int_{\Omega} (v^h - \frac{\partial u^h}{\partial t}) \cdot q^h dx = 0 \quad \forall q^h \in H^h, \\ 
\lambda_N^h \in \Lambda_N^h, \quad \langle \mu_N^h - \lambda_N^h, u_N^h \rangle_{X'_N, X_N} \ge 0, \quad \forall \mu_N^h \in \Lambda_N^h, \\ 
u^h(0) = u_0^h, \quad v^h(0) = v_0^h, 
\end{cases} 
\tag{21} $$

where $V^h, H^h$ and $X_N^h$ are some given finite element spaces approximating $V, L^2(\Omega)$ and $X'_N$, respectively, and $\Lambda_N^h \subset X_N^h$ is a given approximation of $\Lambda_N$. Note that the well-posedness and the quality of the approximation is submitted to the verification of an inf-sup condition between $X_N^h$ and $V^h$ (see [14] for instance). The matrix expression is the following:

$$ egin{cases} 
B^T \dot{V}(t) + AU(t) = L + B_N^T \lambda_N(t), \quad \forall t \in (0, T], \\ 
CV(t) = B \dot{U}(t), \quad \forall t \in (0, T], \\ 
\lambda_N(t) \in \bar{\Lambda}_N^h, \quad (\lambda_N(t) - \mu_N)^T B_N U \ge 0, \quad \forall \mu_N \in \bar{\Lambda}_N^h, \quad \forall t \in (0, T], \\ 
U(0) = U_0, \quad V(0) = V_0, 
\end{cases} 
\tag{22} $$

where the matrices $A, B$ and $C$ have the same meaning as in Section 2 and matrix $B_N$ represent the discrete trace operator on $\Gamma_C$. The set $\bar{\Lambda}_N^h$ is the set corresponding to $\Lambda_N^h$ in the matrix expression. The unknown $V(t)$ can still be eliminated using $V(t) = C^{-1}B\dot{U}(t)$. It follows that a midpoint scheme for the time integration of Problem (22) can be written
with $M = B^T C^{-1} B$:

$$
\left\{
\begin{array}{l}
U^0 \text{ and } V^0 \text{ be given. For } n \ge 0, \text{ find } U^{n+\frac{1}{2}} \text{ such that} \\
MZ^{n+\frac{1}{2}} + AU^{n+\frac{1}{2}} = L + B_N^T \lambda_N^{n+\frac{1}{2}}, \\
\lambda_N^{n+\frac{1}{2}} \in \bar{\Lambda}_N^h, \quad (\lambda_N^{n+\frac{1}{2}} - \mu_N)^T B_N U^{n+\frac{1}{2}} \ge 0, \quad \forall \mu_N \in \bar{\Lambda}_N^h, \quad \forall t \in (0, T], \\
U^{n+\frac{1}{2}} = \frac{U^n + U^{n+1}}{2}, \quad V^{n+\frac{1}{2}} = \frac{V^n + V^{n+1}}{2}, \\
U^{n+1} = U^n + dtV^{n+\frac{1}{2}}, \quad V^{n+1} = V^n + dtZ^{n+\frac{1}{2}}.
\end{array}
\right.
\qquad (23)
$$

<figure><img src="image_14.png" /><figcaption>Figure 13: *Mesh of the disc (with quadratic elements).*</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 14: *Evolution of the energy, the displacement and the contact stress at the point A for a standard finite element semi-discretization and a midpoint scheme.*</figcaption></figure>

We present now some numerical experiments performed on a disc (plane strain approximation of a cylinder) represented on Fig. 13 with the mesh used for the first experiment. This mesh has curved (quadratic) elements on the boundary. This situation is similar to the one presented in [13]. The main parameters of the problem and of the discretization are summarized in Table 1. We denote by $A$ the lowest point of the disc (the point which is the first to come into contact with the foundation).

<table><caption>Table 1: *Characteristics of the elastic disc and the resolution method.*</caption><tr><td>ρ</td><td>10<sup>3</sup>kg/m<sup>3</sup>,</td></tr><tr><td>diameter</td><td>2cm,</td></tr><tr><td>Lamé coefficients</td><td>λ = 300 MPa, μ = 150 MPa</td></tr><tr><td>u<sup>0</sup>, v<sup>0</sup></td><td>0.01 m, -1 m/s</td></tr><tr><td>Simulation time</td><td>0.02 s</td></tr><tr><td>Mesh size</td><td>≈ 2 × 10<sup>-3</sup> m</td></tr></table>

In Fig. 14 the experiment correspond to a classical $P_2$ Lagrange finite element method (a standard mass matrix) with a contact condition on each finite element node lying on the boundary of the mesh. Once again, even though the midpoint scheme is stable (and energy conserving) for the linearized elastodynamic problems, it is not stable for the contact problem. It is remarkable that for the treated case, the discrete problem is close to be energy conserving (first graphic of Fig. 14) for a time step equal to $10^{-4}$s but completely
instable for a smaller time step equal to $10^{-5}s$ (the energy is blowing up very rapidly at the beginning of the impact). The second graphic of Fig. 14 show the vertical displacement at the point A. For the readability of the graphic, the curve for $dt = 10^{-5}s$ is translated with an increment of 1. The instabilities clearly appears for $dt = 10^{-5}s$.

<figure><img src="image_16.png" /><figcaption>Figure 15: Evolution of the energy, the displacement and the contact stress at the point A for a $P_2/P_1$ method, a midpoint scheme and contact condition at the vertices.</figcaption></figure>

The third graphic represents the contact stress at point A. The contact stress for $dt = 10^{-5}s$ is not represented because it is too much noisy. It can be seen that for the contact stress the instabilities already appear for $dt = 10^{-4}s$.

In Fig. 15 the experiment corresponds to the proposed method with a $P_2$ Lagrange finite element for the displacement, a $P_1$ Lagrange finite element method for the velocity and a contact condition applied on each vertex of the boundary of the mesh. The scheme is now perfectly stable and the convergence toward an energy conserving solution is numerically verified when the time step is decreasing. Note that on the third graphic the contact stress at the point A seems also to converge (for the readability of the graphic, curves for $dt = 10^{-4}s$ and $dt = 10^{-5}s$ are translated). We can conclude that the proposed method for elastodynamic contact problems present the same advantages than the mass redistribution method (see the numerical experiments in [13]).

A difference with the mass redistribution method is that here, compared to a classical approximation, the modification of the mass matrix occurs also in the interior of the domain. However, in the framework of the proposed method, it should be possible to make an approximation of the velocity which is different of the one of the displacement only on the elements lying on the boundary of the domain. For instance, this could be achieved by considering the same approximation for the velocity and the displacement except on the
contact boundary either by removing a selection of velocity degrees of freedom or by adding some displacement degrees of freedom (stabilization by bubble functions). This would also lead to a mass matrix which is modified only on the boundary of the domain.

<table><caption>Table 2: Discretization parameters for the convergence test. The relative $H^1$ error is computed on the final time step in comparison with the numerical solution on the finest mesh.</caption><thead><tr><th>dt</th><th>h</th><th>number of d.o.f.</th><th>relative H<sup>1</sup> error</th></tr></thead><tbody><tr><td>8 × 10<sup>-4</sup></td><td>2mm</td><td>622</td><td>5.86%</td></tr><tr><td>4 × 10<sup>-4</sup></td><td>1mm</td><td>2782</td><td>4.41%</td></tr><tr><td>2 × 10<sup>-4</sup></td><td>0.5mm</td><td>11306</td><td>2.57%</td></tr><tr><td>10<sup>-4</sup></td><td>0.25mm</td><td>45706</td><td>1.42%</td></tr><tr><td>5 × 10<sup>-5</sup></td><td>0.14mm</td><td>126892</td><td>0.43%</td></tr><tr><td>2.5 × 10<sup>-5</sup></td><td>0.07mm</td><td>511386</td><td>—</td></tr></tbody></table>

Finally, Fig. 16 presents a convergence test, still with the $P_2/P_1$ method for $h$ (the mesh size) and $dt$ decreasing together. The discretization parameters are listed in Table 2. The relative $H^1$ error is computed on the solution at the last time step (for $t = 0.02s$). In the absence of an exact solution, the comparison is done with the numerical solution on the finest mesh ($h = 0.07mm$) taken to be the reference solution. The first graphic of Fig. 16 represents the evolution of energy for three numerical solutions. There is a numerical convergence toward an energy conserving solution. This is confirmed in the second graphic which represents the evolution of the stress at the point A (still with a vertical shift for $dt = 10^{-4}$ and $2.5 \times 10^{-5}$). There is a fairly good convergence in spite of the fact that we consider a constraint on a single point. The third graphics represent the relative $H^1$ error. The solution for $h = 0.14mm$ is omitted because it is to close to the reference solution ($h = 0.07mm$). Given the low regularity of the solution, a convergence rate of 0.69 is rather satisfactory.
<figure><img src="image_17.png" /><figcaption>Figure 16: A convergence test for the $P_2/P_1$ method, a midpoint scheme and contact condition at the vertices. The mesh size $h$ and $dt$ are decreasing together.</figcaption></figure>

# 6 Concluding remarks

The classical space semi-discretizations of second order hyperbolic problems with constraints lead to ill-posed problems. Moreover, among the solutions of the standard space semi-discretization problem, none of them are really satisfactory. On the one hand, the solutions corresponding to a non-elastic contact (which can be approximated for instance with a backward Euler scheme) may not converge toward an energy conserving solution. On the other hand, the energy conserving solutions have a very oscillatory normal displacement on the contact area which is non-physical and mesh dependent. This is why the great majority of time integration schemes are unstable when they are applied to such semi-discretizations, and why, even with stable time integration schemes, classical space semi-discretizations can lead to non physical solutions. This analysis clearly shows that the problem cannot be resolved at the level of the definition of the time integration scheme. This is indeed the space discretization which has to be adapted. The proposed strategy allows to have well-posed semi-discretizations and ensure that the standard time integration schemes converge toward an energy conserving solution at least for a fixed space semi-discretization. This study identifies the concept of an inf-sup condition linking the approximation of the velocity, the one of the displacement and the one of the convex of constraints. A perspective of this work is to extend the mathematical analysis to more general sets of constraints for instance by replacing the condition (7) by a condition on the tangent cone at each instant. Note that we did not discuss the overall stability of the full discretization, which is also a perspective of this work.
An advantage compared to the mass redistribution method presented in [13, 8] is that here no artificial modification of the mass matrix is necessary. Moreover arbitrary order methods can be obtained and a mass is still present on the contact boundary. Finally, contrary to the mass redistribution method the presented strategy can be applied to thin structures.

## References

[1] P. ALART, A. CURNIER. *A mixed formulation for frictional contact problems prone to Newton like solution methods*. Comp. Meth. Appl. Mech. Engng., 92 (1991), pp 353-375.

[2] M. BARBOTEU. *An efficient algorithm to solve large non linear elastodynamic problems with contact friction*. Bull. Math. Soc. Sc. Math. Roumanie, 48:96:2(2005), pp119-137.

[3] P.W. CHRISTENSEN, J.S. PANG, *Frictional contact algorithms based on semi-smooth Newton methods*. In Reformulation: nonsmooth, piecewise smooth, semi-smooth and smoothing methods. Appl. Optim, vol. 22 (1999), pp 81-116, Kluwer.

[4] P.G. CIARLET, *The finite element method for elliptic problems*, Studies in Mathematics and its Applications No 4, North Holland, 1978.

[5] K. DEIMLING, *Multivalued Differential Equations*, de Gruyter, 1992.

[6] P. DEUFLHARD, R. KRAUSE, S. ERTEL. *A Contact-Stabilized Newmark Method for dynamical contact problems*. Int. J. Num. Meth. Eng, 73(2007), 9:1274-1290.

[7] N. J. CARPENTER. *Lagrange constraints for transient finite element surface contact*, Int. J. Num. Meth. Eng, 32 (1991), pp 103-128.

[8] C. HAGER, B. WOHLMUTH. *Analysis of a modified mass lumping method for the stabilization of frictional contact problems*, to appear in SIAM J. Numer. Anal.

[9] P. HAURET. *Numerical methods for the dynamic analysis of two-scale incompressible nonlinear structures*, Thèse de Doctorat, Ecole Polytechnique, France, 2004.

[10] P. HAURET, P. LE TALLEC. *Energy controlling time integration methods for nonlinear elastodynamics and low-velocity impact*, Comput. Meth. Appl. Mech. Engrg., 195 (2006), pp. 4890-4916.

[11] T. J.R. HUGUES, R. L. TAYLOR, J. L. SACKMAN, A. CURNIER, W. KANOK-NUKULCHAI. *A finite method for a class of contact-impact problems*, 8 (1976), 249-276.

[12] H.B. KHENOUS. *Problèmes de contact unilatéral avec frottement de Coulomb en élastostatique et élastodynamique. Etude mathématique et résolution numérique*, PhD thesis, INSA de Toulouse, France, 2005.

[13] H.B. KHENOUS, P. LABORDE, Y. RENARD. *Mass redistribution method for finite element contact problems in elastodynamics*. Eur. J. Mech., A/Solids, 27(5) (2008), 918-932.

[14] H.B. KHENOUS, J. POMMIER, Y. RENARD. *Hybrid discretization of the Signorini problem with Coulomb friction. Theoretical aspects and comparison of some numerical solvers*. Applied Numerical Mathematics, 56:2 (2006), pp 163-192.
[15] J.U. KIM, *A boundary thin obstacle problem for a wave equation*, Com. part. diff. eqs., 14 (8&9) (1989), 1011–1026.

[16] T.A. LAURSEN, V. CHAWLA, *Design of energy conserving algorithms for frictionless dynamic contact problems*. Int. J. Num. Meth. Engrg, 40 (1997), 863–886.

[17] T.A. LAURSEN, G.R. LOVE. *Improved implicit integrators for transient impact problems-geometric admissibility within the conserving framework*. Int. J. Num. Meth. Eng, 53 (2002), 245–274.

[18] G. LEBEAU, M. SCHATZMAN. *A wave problem in a half-space with a unilateral constraint at the boundary*. J. diff. eqs., 55 (1984), 309–361.

[19] J.J. MOREAU. *Liaisons unilatérales sans frottement et chocs inélastiques*. C.R.A.S. série II, 296 pp 1473-1476, 1983.

[20] J.J. MOREAU, *Numerical aspects of the sweeping process*. Comp. Meth. Appl. Mech. Engrg., 177 (1999), 329–349.

[21] L. PAOLI. *Time discretization of vibro-impact*. Phil. Trans. R. Soc. Lond. A., 359 (2001), 2405–2428.

[22] L. PAOLI, M. SCHATZMAN. *Approximation et existence en vibro-impact*. C. R. Acad. Sci. Paris, Sér. I, 329 (1999), 1103–1107.

[23] Y. RENARD, J. POMMIER. *Getfem++*. An Open Source generic C++ library for finite element methods, http://home.gna.org/getfem.

[24] K. SCHWEIZERHOF, J.O. HALLQUIST, D. STILLMAN. *Efficiency Refinements of Contact Strategies and Algorithms in Explicit Finite Element Programming*. In Computational Plasticity, eds. Owen, Onate, Hinton, Pineridge Press, pp 457–482, 1992.

[25] R. L. TAYLOR, P. PAPADOPOULOS. *On a finite element method for dynamic contact-impact problems*. Int. J. for Num. Meth. Eng., Vol 36 (1993), pp 2123–2140.