# Time discretization of vibro-impact

Laetitia Paoli

► To cite this version:

Laetitia Paoli. Time discretization of vibro-impact. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2001, 359, pp.2405-2428. 10.1098/rsta.2001.0858. hal-01566940

HAL Id: hal-01566940

https://hal.science/hal-01566940

Submitted on 21 Jul 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
# Time discretization of vibro-impact

LAETITIA PAOLI

Equipe d'Analyse Numérique de Saint Etienne et MAPLY UMR 5585 - CNRS,
Université Jean Monnet, 23 rue du Docteur Paul Michelon,
Saint Etienne 42023 Cedex 2, France

We consider vibro-impact problems, i.e. mechanical systems with a finite number of degrees of freedom submitted to perfect unilateral constraints. The dynamics is basically described by a second-order measure differential inclusion for the unknown position completed with a constitutive impact law. Another formulation of the problem as a frictionless sweeping process is possible: the unknown velocity belongs to an appropriate functional space and satisfies a first order measure differential inclusion.

The equivalence of these two formulations is studied. They lead to time-discretizations written in terms of positions or in terms of velocities, respectively. We present these different schemes and we compare them on the simple test-problem of a bouncing ball. We recall the convergence results in the single constraint case. Moreover, an example of implementation of the scheme derived from the basic description of the dynamics is presented. Finally, in the multi-constraint case, we point out some theoretical and computational difficulties.

Keywords: vibro-impact; frictionless sweeping process; time discretization; convergence; single/multi constraints

## 1. Introduction: the basic description of the dynamics

We consider a dynamical system with $d$ degrees of freedom submitted to perfect (i.e. frictionless) unilateral constraints. We denote by $q \in \mathbb{R}^d$ the representative point of the system and we assume that the constraints are described by the following inequalities,

$$ \varphi_{\alpha}(q) \geq 0, \quad \alpha \in \{1, \dots, \ell\}, $$

with smooth functions $\varphi_{\alpha}$ (at least $C^1$) such that $d\varphi_{\alpha}$ does not vanish in a neighbourhood of $\{q : \varphi_{\alpha}(q) = 0\}$.

The set of admissible positions is given by

$$ K = \{q \in \mathbb{R}^d : \varphi_{\alpha}(q) \geq 0 \ \forall \alpha \in \{1, \dots, \ell\}\}. \quad (1.1) $$

For all $q \in \mathbb{R}^d$ we define the set of active constraints by

$$ J(q) = \{\alpha \in \{1, \dots, \ell\} : \varphi_{\alpha}(q) \leq 0\}, $$

and we extend the definition of the tangent cone to $K$ at $q$ by setting

$$ \forall q \in \mathbb{R}^d \quad T_K(q) = \{w \in \mathbb{R}^d : d\varphi_{\alpha}(q)w \geq 0 \ \forall \alpha \in J(q)\}. $$

We assume moreover that the active constraints are functionally independent, i.e.

$$ \forall q \in K \quad (d\varphi_{\alpha}(q))_{\alpha \in J(q)} \text{ are linearly independent.} $$
Whenever $q(t)$ belongs to $\text{Int}(K)$ the motion is described by a second-order ordinary differential equation (ODE),

$$M(q)\ddot{q} = f(t, q, p), \quad p = M(q)\dot{q},$$

where $M(q)$ is the mass matrix of the system. We should note that dissipative terms may be included in this formulation.

If $q(t) \in \text{Int}(K)$ for all $t \in (t_0, t_1) \cup (t_1, t_2)$ with $q(t_1) \in \partial K$, we obtain

$$\dot{q}(t_1 + 0) \in T_K(q(t_1)), \quad -\dot{q}(t_1 - 0) \in T_K(q(t_1)).$$

Hence, if $\dot{q}(t_1 - 0) \notin T_K(q(t_1))$, a discontinuity in the velocity will occur. We denote $\nu$ the reaction of the constraints, we have

$$M(q)\ddot{q} = f(t, q, p) + \nu, \quad p = M(q)\dot{q}, \qquad (1.2)$$

and, since we expect discontinuities in the velocity at impacts, $\nu$ is a measure such that

$$\text{supp}(\nu) \subset \{t: q(t) \in \partial K\}. \qquad (1.3)$$

Moreover, the constraints are perfect and we get (Moreau 1963)

$$\nu = \sum_{\alpha=1}^{\ell} \lambda_{\alpha} \mathrm{d}\varphi_{\alpha}(q) \quad \text{with} \quad \begin{cases} \lambda_{\alpha} \ge 0 & \text{if } \alpha \in J(q) \\ \lambda_{\alpha} = 0 & \text{if } \alpha \notin J(q) \end{cases} \qquad (1.4)$$

are perfect.

These equations do not describe completely the motion, since the coefficients $(\lambda_{\alpha})_{\alpha \in J(q)}$ are unknown and we should add a constitutive law of impact.

In the case of a single active constraint, i.e. $J(q(t)) = \{\alpha\}$ with $\alpha \in \{1, \dots, \ell\}$, we decompose the momenta $p(t \pm 0)$ as

$$p(t \pm 0) = p_N(t \pm 0) + p_T(t \pm 0)$$

with

$$\left. \begin{array}{l} p_N(t \pm 0) \in \mathbb{R} d\varphi_{\alpha}(q(t)), \\ d\varphi_{\alpha}(q(t))(M(q(t))^{-1}p_T(t \pm 0)) = 0. \end{array} \right\} \qquad (1.5)$$

Let us denote by $\langle \cdot, \cdot \rangle_q$ and $\langle \cdot, \cdot \rangle_q^*$ the inner products given by

$$\langle \xi, \eta \rangle_q = {}^t\xi M(q)\eta \quad (\text{kinetic local metric}),$$

$$\langle \xi, \eta \rangle_q^* = {}^t\xi M(q)^{-1}\eta \quad (\text{momentum local metric}).$$

The corresponding norms are $|\cdot|_q$ and $|\cdot|_q^*$. Relations (1.5) define the normal and tangential components of $p(t \pm 0)$ relative to the local momentum metric.

The inequality $\varphi_{\alpha}(q(s)) \ge 0$ for all $s$ and equation (1.4) imply that

$$p_T(t + 0) = p_T(t - 0),$$

$$p_N(t + 0) \in \mathbb{R}^+ d\varphi_{\alpha}(q(t)), \quad p_N(t - 0) \in \mathbb{R}^- d\varphi_{\alpha}(q(t)).$$

If we assume moreover that the kinetic energy does not increase at impacts, we obtain

$$\begin{aligned} 2T(t + 0) &= |p_N(t + 0)|_{q(t)}^{*2} + |p_T(t + 0)|_{q(t)}^{*2} \\ &\le |p_N(t - 0)|_{q(t)}^{*2} + |p_T(t - 0)|_{q(t)}^{*2} = 2T(t - 0) \end{aligned}$$
and thus

$$|p_N(t+0)|_{q(t)}^* \le |p_N(t-0)|_{q(t)}^*.$$

Hence, the constitutive law of impact is given by

$$p_T(t+0) = p_T(t-0), \quad p_N(t+0) = -ep_N(t-0) \quad \text{if } q(t) \in \partial K, \quad (1.6)$$

with a restitution coefficient $e \in [0, 1]$.

The kinetic energy is conserved if and only if $e = 1$ (totally elastic shocks), and for $e = 0$ we find the inelastic shocks defined in Moreau (1983), since (1.6) yields

$$\dot{q}(t+0) = \text{Proj}_{q(t)}(T_K(q(t)), \dot{q}(t-0)),$$

where $\text{Proj}_q(T_K(q), z)$ is the projection of $z$ on $T_K(q)$ relative to the local kinetic metric at $q$. When $e \in (0, 1)$, we obtain partly elastic shocks and the relation between the restitution coefficient and the 'dissipation index' $\delta$ of Moreau (1988) is given by

$$e = \frac{1 - \delta}{1 + \delta}.$$

When the cardinal of $J(q(t))$ is greater than 1 we can decompose $p(t \pm 0)$ as

$$p(t \pm 0) = p_N(t \pm 0) + p_T(t \pm 0)$$

with

$$p_N(t \pm 0) \in \text{Span}\{d\varphi_\alpha(q(t)), \alpha \in J(q(t))\}, \\ d\varphi_\alpha(q(t))(M(q(t))^{-1}p_T(t \pm 0)) = 0 \quad \forall \alpha \in J(q(t)),$$

which is a generalization of (1.5) and we find again

$$p_T(t+0) = p_T(t-0), \quad |p_N(t+0)|_{q(t)}^* \le |p_N(t-0)|_{q(t)}^*. \quad (1.7)$$

But now (1.7) does not necessarily yield an impact law like (1.6) because $p_N(t+0)$ and $p_N(t-0)$ are not necessarily collinear.

We will come back to this case in § 5, and in the next three sections we assume that

$$\forall t \quad \text{Card} \, J(q(t)) \le 1. \quad (1.8)$$

We are interested in the Cauchy problem: for given admissible initial data $(t_0, q_0, v_0) \in \mathbb{R} \times K \times T_K(q_0)$ we look for a solution of (1.2)–(1.4) and (1.6). More precisely, we consider the following problem.

**Problem (P).** To find $\tau > 0$ and a function $q$ from $[t_0, t_0 + \tau]$ to $\mathbb{R}^d$ satisfying the following variational properties,

$$\left. \begin{array}{l} (i) \quad q \text{ is continuous with values in } K, \\ (ii) \quad \dot{q} \text{ is a function of bounded variation,} \end{array} \right\} \quad (1.9)$$

such that (1.2)–(1.4) and (1.6) hold on $(t_0, t_0 + \tau)$ and

$$q(t_0) = q_0, \quad \dot{q}(t_0 + 0) = v_0. \quad (1.10)$$

**Remark 1.1.** With the regularity assumption (1.9)(ii), $\dot{q}(t+0)$ is defined for all $t \in [t_0, t_0 + \tau)$, $\dot{q}(t-0)$ is defined for all $t \in (t_0, t_0 + \tau]$ and $\ddot{q}$ in (1.2) denotes the Stieltjes measure of $\dot{q}$.
Under assumption (1.8) we can restrict ourselves, without loss of generality, to the case of a set K given by

$$K = \{q \in \mathbb{R}^d : \varphi(q) \ge 0\} \quad (\ell = 1). \qquad (1.11)$$

In § 2 we recall briefly the mathematical background (existence and uniqueness results) and the motivations for numerical methods based on time discretization. In § 3 we present a scheme derived directly from problem (P) and we apply it to compute the motion of a guided slender pipe. Then in § 4, we show that the basic description of vibro-impact problems is equivalent to a frictionless sweeping process and we present the time discretizations derived from this second formulation.

## 2. Single constraint case: mathematical background and motivations for time discretization

The existence and uniqueness of solutions for the Cauchy problem (P) under assumption (1.11) is still an open problem in the general case. One of the oldest contributions in this field is due to Bressan (1959) with a counterexample to uniqueness. The first existence results in the functional framework defined by (1.9) had been proved by Schatzman (1977, 1978). She considered only the case of totally elastic shocks and she assumed that the set of admissible positions is convex and the mass matrix is trivial (thus the kinetic local metric is simply the Euclidean metric of $\mathbb{R}^d$). She used a penalty method which allows very weak regularity for $\partial K$.

For totally elastic shocks, more general results had been obtained by Buttazzo & Percivale (1981, 1983) and Percivale (1986), uniqueness had been studied by Carriero & Pascali (1982, 1985) in the one-dimensional case and by Percivale (1985, 1991) in the general case.

For inelastic shocks, Monteiro Marques (1985, 1993) gave a proof of existence using a time discretization under the assumption of a trivial mass matrix.

Since the beginning of the 1990s, all these results have been extended to the case of partly elastic shocks. In Paoli & Schatzman (1993a, c) a generalization of Schatzman (1978) is obtained for any restitution coefficient $e$ belonging to $(0, 1]$. A time discretization of the dynamics is proposed in Paoli & Schatzman (1993b, 1999a) which allows the treatment of the case $e = 0$ as well as the general case of a non-trivial mass matrix and/or a non-convex set of admissible positions. The convergence of this scheme, which yields an existence result, is proved in Paoli (1993) and Paoli & Schatzman (1999b).

Another discretization, inspired by Monteiro Marques (1993), is given in Mabrouk (1998a), the convergence being established only when the mass matrix is trivial in Mabrouk (1998b).

Finally, in the more general case of a restitution coefficient depending on $t$ and $q$, an existence result when all the data are analytic is stated in Ballard (1999). Uniqueness is also studied and a generalization of a previous work of Schatzman (1998) is obtained.

We can consider three different approaches to compute the solution.

1. Event-driven schemes: the idea is to integrate by any method (exact or approximate solution) the ODE of the unconstrained motion, to determine the impact times and to apply the impact law to restart the algorithm.

2. Compliant models: the idea is to relax the constraint in order to replace the measure differential inclusion (1.2)–(1.4) by a (stiff) ODE.
3. Time discretization: the idea is to provide a kind of Euler discretization of the measure differential inclusion which describes the dynamics, including some information about the restitution coefficient and avoiding the precise determination of the impact times.

The first approach is conceptually the most simple and is very convenient when the motion can be decomposed into a finite number of intervals of regular motion. Unfortunately, since the work of Delassus (1917), it is known that more complex dynamical behaviours exist, with impacts accumulation, for instance. In such a case the numerical detection of all the impacts is impossible and a threshold $\delta$ must be chosen: we decide that the constraint is active is $\varphi(q(t)) \le \delta$. In order to keep a good approximation of the global dynamics of the system, $\delta$ should be chosen to be small enough and thus we still have to detect a lot of impacts. Since for each detection we have to solve the equation $\varphi(q(t)) = \delta$ this approach may become very expensive.

The second approach (compliant models) is more sophisticated and is mechanically and mathematically very attractive. From a mechanical point of view, a body is never perfectly rigid and the system (1.2)–(1.4) and (1.6) describes a limit behaviour when the stiffness tends to $+\infty$. Hence it seems very natural to approximate the problem by replacing the rigid boundary of $K$ by an elastic one with a stiffness $k \gg 1$. From a mathematical point of view, this idea is substantiated by the convergence of the penalized solutions which has been established in Schatzman (1978), Paoli & Schatzman (1993c) and Paoli (2000).

But for the implementation of compliant models we have to fix the value of $k$, which is a physical parameter depending on the stiffness of the impacting bodies. The estimation of this stiffness is experimentally quite difficult but we should expect large values: $k = 10^{10} \text{ N m}^{-1}$ for systems with joint clearance (Ravn 1998), $k = 5.5 \times 10^7 \text{ N m}^{-1}$ for an impacting bar (Stoianovici & Hurmuzlu 1996). Since the penalized solutions leave $K$ during a time of order $O(1/\sqrt{k})$ (Paoli & Schatzman 1993c), we should choose a time-step $h \le O(1/\sqrt{k})$. Thus this approach is very expensive too. Moreover, the collision outcome may be quite sensitive with respect to these parameters (Brogliato 1999; Paoli & Schatzman 2000).

The third approach is interesting for theoretical and computational reasons. The time discretizations proposed in Paoli & Schatzman (1993b, 1999a), Monteiro Marques (1985, 1993) and Mabrouk (1998a, b) deal directly with the constraint and the impact law: they do not require a systematic detection of impacts and no refinement near the boundary of $K$ is needed. Thus they avoid the computational drawbacks of the previous approaches. Moreover, their convergence is stated in Paoli (1993), Paoli & Schatzman (1999b), Monteiro Marques (1985, 1993) and Mabrouk (1998b).

In the next two sections we describe successively these two schemes. We do not present in this paper the time discretization proposed by Stewart (1997, 1998), since it only allows the treatment of inelastic impacts.

### 3. Time discretization of problem (P)

#### (a) Description of the scheme

This discretization of equations (1.2)–(1.4) and (1.6) is written in terms of positions only. Thus the relevant metric is the kinetic metric of the system. In order to be able to define, at least locally, the projection $P_K$ on $K$ relative to this metric, we introduce the following technical assumptions:
(H1) the function $\varphi$ is of class $C^3$,

(H2) the mapping $q \mapsto M(q)$ is of class $C^3$ from $\mathbb{R}^d$ to the set of symmetric positive definite matrices.

We assume moreover that

(H3) the function $f$ is continuous from $[0, T] \times \mathbb{R}^d \times \mathbb{R}^d$ ($T > 0$) to $\mathbb{R}^d$, and is locally Lipschitz continuous with respect to its last two arguments.

Let $F$ be a continuous function from $[0, T] \times \mathbb{R}^d \times \mathbb{R}^d \times \mathbb{R}^d \times [0, h^*]$ ($h^* > 0$) to $\mathbb{R}^d$ such that

(H4) $F$ is locally Lipschitz continuous with respect to its second, third and fourth arguments and

$$ \forall (t, q, v) \in [0, T] \times (\mathbb{R}^d)^2 \quad F(t, q, q, v, 0) = M(q)^{-1} f(t, q, M(q)v). $$

For given initial data $(t_0, q_0, v_0) \in [0, T] \times K \times T_K(q_0)$ we define the following approximation of the Cauchy problem (P):

$$ q^0 = q_0, \quad q^1 = q_0 + hv_0 + hz(h) \quad \text{with} \quad \lim_{h \to 0} z(h) = 0, \qquad (3.1) $$

and, for all $i \ge 1$,

$$ q^{i+1} = -eq^{i-1} + (1+e)P_K \left( \frac{2q^i - (1-e)q^{i-1} + h^2 F^i}{1+e} \right), \qquad (3.2) $$

with

$$ F^i = F \left( t_0 + ih, q^i, q^{i-1}, \frac{q^{i+1} - q^{i-1}}{2h}, h \right), \qquad (3.3) $$

provided that $q^{i+1}$ is uniquely defined in a neighbourhood of $q^i$.

At first glance this discretization does not seem much related to problem (P) and it is necessary to explain its construction. This scheme has been proposed first in the case of a trivial mass matrix (i.e. $M(q) \equiv \text{Id}_{\mathbb{R}^d}$) with a convex set of admissible positions (Paoli & Schatzman 1993b). Under these assumptions equations (1.2)–(1.4) can be rewritten as

$$ \ddot{q} + \partial\psi_K(q) \in f(t, q, \dot{q}), $$

where $\psi_K$ denotes the indicatrix function of $K$ and $\partial\psi_K$ is its subdifferential (Rockafellar 1970). The idea is to apply a centred scheme to the ODE $\ddot{q} = f(t, q, \dot{q})$ and to add a discretization of the term $\partial\psi_K(q)$. Thus we propose

$$ \frac{q^{i+1} - 2q^i + q^{i-1}}{h^2} + \partial\psi_K \left( \frac{q^{i+1} + eq^{i-1}}{1+e} \right) \ni F^i, $$

which we rewrite as equation (3.2) by using the classical properties of convex analysis. The constraint is satisfied at each time-step by the average position

$$ \bar{q}^i = \frac{q^{i+1} + eq^{i-1}}{1+e}, $$
and whenever $\bar{q}^i$ belongs to Int($K$) equation (3.2) reduces to

$$q^{i+1} = 2q^i - q^{i-1} + h^2 F^i,$$

which is a second-order discretization of the unconstrained dynamics. Moreover, the weights involved in this average position yield a correct reflection of the velocities at impacts. Let us illustrate this property with the test-problem of a bouncing ball: $d = 1$, $K = \mathbb{R}^+$, $M(q) = 1$, $f = 0$ and $(t_0, q_0, v_0) = (0, 1, -1)$. The solution is defined by

$$q(t) = \begin{cases} 1-t & \text{for all } t \in [0, 1], \\ e(t-1) & \text{for all } t \ge 1. \end{cases}$$

We choose $F \equiv 0$, $q^0 = 1$ and $q^1 = 1 - h$. The algorithm is given by

$$\forall i \ge 1 \qquad q^{i+1} = -eq^{i-1} + \max(2q^i - (1-e)q^{i-1}, 0). \qquad (3.4)$$

Let $h < (1 + e)/2$; there exists a greatest integer $p \ge 2$ such that

$$\forall i \in \{1, \dots, p-1\} \qquad 2q^i - (1-e)q^{i-1} \ge 0.$$

With (3.4) we obtain

$$\forall i \in \{1, \dots, p-1\} \qquad q^{i+1} - q^i = q^i - q^{i-1} = -h,$$

and, by definition of $p$, $q^{p+1} = -eq^{p-1}$. It follows that

$$2q^{p+1} - (1-e)q^p = -2eq^{p-1} - (1-e)(2q^{p-1} - q^{p-2}) \\ = -(2q^{p-1} - (1-e)q^{p-2}) \le 0,$$

and thus $q^{p+2} = -eq^p$. Moreover,

$$2q^{p+2} - (1-e)q^{p+1} = -e(2q^p - (1-e)q^{p-1}) \ge 0,$$

and we obtain by induction that

$$\forall i \ge p+2 \qquad q^i = q^{p+2} + eh(i-p-2).$$

Hence the velocity is reversed in two time-steps and the approximate impact time $t_p = ph$ satisfies

$$1 - \frac{1-e}{1+e}h < t_p \le 1 + \frac{2eh}{1+e}.$$

It follows that the constraints are violated during at most two time-steps.

### (b) Convergence result

The convergence of this scheme has been proved first in the case of a convex set of constraints with a trivial mass matrix in Paoli (1993), then an extension to non-convex constraints (but still trivial mass matrix) has been proposed in Paoli & Schatzman (1998) and finally the general case has been considered in Paoli & Schatzman (1999b). Let us describe the sketch of the proof.

First we show that there exists $\tau > 0$ such that the discrete velocities and accelerations remain uniformly bounded on $[t_0, t_0 + \tau]$. We define the approximate solutions
$q_h$ by linear interpolation of the $q^i$. The Ascoli–Arzela theorem and Helly’s theorem
imply that there exists a subsequence, still denoted $q_h$, such that

$$
q_h \to q \quad \text{strongly in} \quad C^0([t_0, t_0 + \tau]; \mathbb{R}^d),
$$

$$
\dot{q}_h \to \dot{q} \quad \text{pointwise in} \quad [t_0, t_0 + \tau],
$$

$$
\ddot{q}_h \to \ddot{q} \quad \text{weakly* in} \quad M^1((t_0, t_0 + \tau); \mathbb{R}^d).
$$

We use here the classical notation of the ‘weak*’ topology for measures, i.e. conver-
gence relative to the duality product with continuous functions. Then we establish
that the function $q$ takes its values in $K$ and that equations (1.2)–(1.4) and (1.6)
and (1.10) are satisfied.

The main difficulties are the proof of uniform estimates for the discrete velocities
and the study of transmission of energy at impacts. In both cases we use systemat-
ically local coordinates which transform K into a half-space. This is a very natural
geometrical idea. Unfortunately, curvature effects create quadratic terms which inter-
act with the unilateral constraint.

More precisely, for all $\bar{q} \in \partial K$, we define locally a diffeomorphism $\Phi$ such that

$$
P_K(q) = \Phi^{-1} \left( \begin{array}{c} S(q) \\ \max(Y(q), 0) \end{array} \right) \quad \text{with} \quad \Phi(q) = \left( \begin{array}{c} S(q) \\ Y(q) \end{array} \right) \in \mathbb{R}^{d-1} \times \mathbb{R}
$$

and we apply $\Phi$ to equations (3.2)–(3.3). We obtain

$$
s^{i+1} = 2s^i - s^{i-1} + h^2 \kappa^i,
$$

$$
y^{i+1} = -ey^{i-1} + \max(2y^i - (1-e)y^{i-1}, 0) + h^2 \lambda^i,
$$

with

$$
\Phi(q^i) = \begin{pmatrix} s^i \\ y^i \end{pmatrix}
$$

and we have the following lemma.

**Lemma 3.1 (the heart of the estimates).** Let $\eta^i = (y^{i+1} - y^i)/h$. We have

$$
|\eta^i| \le \max(|\eta^{i-1}|, e^{|\eta^{i-2}|}) + h^2|\lambda^i| + h^2|\lambda^{i-1}|. \quad (3.5)
$$

We prove that ($|\kappa^i|$)$_{0\le i\le \lfloor\tau/h\rfloor-1}$ and ($|\lambda^i|$)$_{0\le i\le \lfloor\tau/h\rfloor-1}$ remain uniformly bounded
for some $\tau > 0$: these estimates are rather technical since $\kappa^i$ and $\lambda^i$ depend on $F^i$,
$q^{i+1}$, $q^i$, $q^{i-1}$ and $\Phi$.

Lemma 3.1 also plays an essential role in the study of the velocity reflection at impacts. Let $\bar{t}$ be such that $q(\bar{t}) \in \partial K$ and $p_N(\bar{t}-0) \neq p_N(\bar{t}+0)$. By using the local coordinates and lemma 3.1 we obtain that $p_N(\bar{t}-0) \neq 0$. Then, applying the same ideas as in the example of the bouncing ball, we prove that the discrete velocities are correctly reversed in two time-steps.

(c) *Implementation*

In this subsection we present an example of computation. We consider a clamped
free tube excited by a shaker and guided (see figure 1). Let us denote by $\rho$ the mass
density, $E$ the Young’s modulus of elasticity, $S$ the cross-sectional area, $L$ the length
and $I$ the cross-sectional moment of inertia of the pipe. We assume that the external
excitement $g(x, t)$ is a concentrated force parallel to the $y$-axis.
<figure><img src="image_2.png" /><figcaption>Figure 1. The mechanical system.</figcaption></figure>

Let $u : [0, T] \times [0, L] \to \mathbb{R}$ ($T > 0$) be the deflection along the y-axis. The unconstrained motion is described by the following partial differential equation (Landau & Lifschitz 1959)

$$ \rho S \frac{\partial^2 u}{\partial t^2} + EI \frac{\partial^4 u}{\partial x^4} = g(x, t) \qquad (3.6) $$

with

$$ \forall t \in [0, T] \quad \begin{cases} u(0, t) = \frac{\partial u}{\partial x}(0, t) = 0 & (\text{clamped end}), \\ \frac{\partial^2 u}{\partial x^2}(L, t) = \frac{\partial^3 u}{\partial x^3}(L, t) = 0 & (\text{free end}). \end{cases} \qquad (3.7) $$

If we assume that the guide abscissa is $x_0$, we have the constraint

$$ \forall t \in [0, T] \quad u(x_0, t) \in [\phi^{-}, \phi^{+}]. \qquad (3.8) $$

First we discretize in space the system (3.6), (3.7) by using finite differences. Let $k^2 = EI/\rho S$, $d = \lfloor L/\Delta x \rfloor$ and

$$ A_{\Delta x} = \frac{1}{\Delta x^4} \begin{pmatrix} 6 & -4 & 1 & 0 & \cdots & \cdots & \cdots & 0 \\ -4 & 6 & -4 & 1 & 0 & \cdots & \cdots & 0 \\ 1 & -4 & 6 & -4 & 1 & 0 & \cdots & 0 \\ 0 & 1 & -4 & 6 & -4 & 1 & 0 & \cdots & 0 \\ \vdots & & & & & & & & \vdots \\ 0 & \cdots & 0 & 1 & -4 & 6 & -4 & 1 & 0 \\ 0 & \cdots & \cdots & 0 & 1 & -4 & 6 & -4 & 1 \\ 0 & \cdots & \cdots & \cdots & 0 & 1 & -4 & 5 & -2 \\ 0 & \cdots & \cdots & \cdots & \cdots & 0 & 1 & -2 & 1 \end{pmatrix} \in \mathcal{M}_{d,d}(\mathbb{R}). $$

The system (3.6), (3.7) yields the following ODE,

$$ \frac{d^2 U_{\Delta x}}{dt^2} + k^2 A_{\Delta x} U_{\Delta x} = G_{\Delta x}, \qquad (3.9) $$

where $G_{\Delta x}$ is an approximation of $g/\rho S$ at the spatial nodes. We define $K$ by

$$ K = \left( \prod_{j=1, j \neq j_0}^{d} \mathbb{R} \right) \times [\phi^{-}, \phi^{+}] \qquad (3.10) $$
with $|j_0\Delta x - x_0| = \min_{1 \le j \le d} |j\Delta x - x_0|$ and the constraint (3.8) is approximated by

$$ \forall t \in [0, T] \quad U_{\Delta x}(t) \in K. \qquad (3.11) $$

We discretize in time the vibro-impact problem (3.9)-(3.11) by using the scheme. We choose $F$ as follows: for all $(t, u_1, u_2, v, h) \in \mathbb{R} \times (\mathbb{R}^d)^3 \times [0, 1]$

$$ F(t, u_1, u_2, v, h) = \frac{1}{4}(G_{\Delta x}(t-h) + 2G_{\Delta x}(t) + G_{\Delta x}(t+h)) - \frac{1}{2}k^2 A_{\Delta x}(u_1 + hv + u_2). $$

With this choice, equation (3.2) reduces to a stable, second-order, Newmark's scheme for (3.6) whenever

$$ \frac{U_{\Delta x}^{i+1} + eU_{\Delta x}^{i-1}}{1 + e} $$

belongs to $\text{Int}(K)$. Moreover, the usual properties of projections yield an explicit formulation of the algorithm: for all $i \ge 1$

$$ (\text{Id}_{\mathbb{R}^d} + \frac{1}{4}h^2k^2A_{\Delta x})(U_{\Delta x}^{i+1} + eU_{\Delta x}^{i-1}) + h^2\partial\psi_K\left(\frac{U_{\Delta x}^{i+1} + eU_{\Delta x}^{i-1}}{1+e}\right) \ni X^i \quad (3.12) $$

with

$$ X^i = \frac{1}{4}h^2(G_{\Delta x}((i+1)h) + 2G_{\Delta x}(ih) + G_{\Delta x}((i-1)h)) \\ + 2U_{\Delta x}^i - \frac{1}{2}k^2h^2A_{\Delta x}U_{\Delta x}^i - (1-e)(\text{Id}_{\mathbb{R}^d} + \frac{1}{4}h^2k^2A_{\Delta x})U_{\Delta x}^{i-1}. \quad (3.13) $$

In order to solve (3.12) we use the following lemma.

**Lemma 3.2.** Let $M$ be a symmetric definite matrix, $X \in \mathbb{R}^d$ and $\delta > 0$. The inclusion

$$ Mu + \delta \partial\psi_K(u) \ni X \qquad (3.14) $$

has a unique solution $u$ given by

$$ \left. \begin{array}{l} u_{j_0} = P_{[\phi^{-}, \phi^{+}]}((M^{-1}X)_{j_0}), \\ (Mu)_j = X_j \quad \forall j \in \{1, \dots, d\} \setminus \{j_0\}. \end{array} \right\} \qquad (3.15) $$

*Proof.* We denote by $(\cdot, \cdot)$ the Euclidean inner product of $\mathbb{R}^d$. By definition of $\partial\psi_K$ we have

$$ \partial\psi_K(u) = \{z \in \mathbb{R}^d : (z, v - u) \le 0 \ \forall v \in K\} \quad \text{if } u \in K, \quad \partial\psi_K(u) = \emptyset \quad \text{if } u \notin K. $$

Thus we infer that $u$ solves (3.14) if and only if $u \in K$ and $(Mu - X, v - u) \ge 0$ for all $v \in K$, i.e. if and only if $u$ is the unique solution of the problem

$$ \left. \begin{array}{l} \text{find } u \in K \text{ such that} \\ \frac{1}{2}(Mu, u) - (X, u) = \min_{v \in K} \frac{1}{2}(Mv, v) - (X, v). \end{array} \right\} \qquad (3.16) $$

Then we just have to check that the vector $u$ given by (3.15) satisfies (3.16) (Paoli 1993).

Consequently, we have to solve at most two linear systems at each time-step of the algorithm (3.12), (3.13).
<figure><img src="image_3.png" /><figcaption>Figure 2. The mapping $t \mapsto u(x_{\text{ob}}, t)$: numerical results with $h = 10^{-4}$.</figcaption></figure>

The numerical results have been compared with experimental results obtained by the AMV Department (Département Acoustique et Mécanique Vibratoire) of Électricité De France (EDF), Clamart (see also Jacquart & Gay (1992) for the description of a numerical code based on compliant models for impact and friction). The numerical and experimental values of the mapping $t \mapsto u(x_{\text{ob}}, t)$ ($x_{\text{ob}}$ is the observation abscissa) are plotted in figures 2 and 3. We have also compared the reaction force $\nu$ given by

$$ \nu = -g(x_0, \cdot) + EI \frac{\partial^4 u}{\partial x^4}(x_0, \cdot) + \rho S \frac{\partial^2 u}{\partial t^2}(x_0, \cdot). $$

The value of $\nu$ at $t_i = ih$ is approximated by

$$ (\rho S)^{-1} \nu_{\Delta x}^{i} = -\left(\frac{1}{4}G_{\Delta x, j_0}(t_{i+1}) + 2G_{\Delta x, j_0}(t_i) + G_{\Delta x, j_0}(t_{i-1})\right) \\ + \frac{U_{\Delta x, j_0}^{i+1} - 2U_{\Delta x, j_0}^{i} + U_{\Delta x, j_0}^{i-1}}{h^2} + \frac{1}{4}k^2 A_{\Delta x}(U_{\Delta x, j_0}^{i+1} + 2U_{\Delta x, j_0}^{i} + U_{\Delta x, j_0}^{i-1}) \\ \text{for } 1 \le i \le \lfloor T/h \rfloor - 1. $$

The numerical and experimental results are plotted in figures 4 and 5.

In both cases we obtain good agreement. We should note that very close impacts occur: for this example the approximation of the motion by event-driven schemes or compliant models would be very expensive.

This time discretization has been implemented on other examples. In Panet et al. (1995) and Paoli & Schatzman (1995, 2000) a model of a tight joint with one degree of freedom is studied, and Paoli & Schatzman (1999c) consider a slender bar in large deformations (the mass matrix is no longer trivial). The results are compared with those obtained by other numerical methods (event-driven scheme
<figure><img src="image_4.png" /><figcaption>Figure 3. The mapping $t \mapsto u(x_{\text{ob}}, t)$: experimental results (Reproduced with permission from Electricité De France, Clamart, France).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. The reaction force $\nu$: numerical results with $h = 10^{-4}$.</figcaption></figure>

in Panet *et al.* (1995) and Paoli & Schatzman (1995, 1999c), compliant models in Paoli & Schatzman (2000) and also with experimental results in the last example (Stoianovici & Hurmuzlu 1996).

**4. Time discretization of the frictionless sweeping process**

The time discretizations presented in this section rely on another formulation of the problem, first introduced by Moreau to describe inelastic shocks (Moreau 1983, 1985). The general case, with a restitution coefficient $e \in [0, 1]$, has been considered in Moreau (1988) and Jean & Moreau (1992).
<figure><img src="image_6.png" /><figcaption>Figure 5. The reaction force $\nu$: experimental results (Reproduced with permission from the AMV Department, Electricité De France, Clamart, France).</figcaption></figure>

In the next subsection we recall Moreau’s formulation and we establish that it is equivalent to the basic description of the dynamics when assumption (1.11) holds. Then we focus on the time discretization derived by Monteiro Marques and Mabrouk, the convergence of which has been proved in Monteiro Marques (1985, 1993) and Mabrouk (1998b).

A more sophisticated time discretization, with a mid-point approximation of the velocity, has been proposed by Moreau (1986, 1988, 1999) and Jean & Moreau (1992). Moreover, with some appropriate modifications, frictional contact can be considered too. These algorithms have been successfully used to simulate granular matter problems (Moreau 1994, 1995), but their convergence has never been completely studied. Thus we will not describe them in this paper.

### (a) Moreau’s sweeping process

Let $\tau$ be a non-negative number. We define the functional space $bvem([t_0, t_0 + \tau])$ by the following definition.

**Definition 4.1.** A function $u : [t_0, t_0 + \tau] \to \mathbb{R}^d$ belongs to $bvem([t_0, t_0 + \tau])$ if and only if $u$ is a function of bounded variation and

$$u(t) = \frac{u(t+0) + eu(t-0)}{1+e}$$

for all $t \in (t_0, t_0 + \tau)$.

Let us denote $\psi_{T_K(q)}$ the indicatrix function of the convex set $T_K(q)$ and $\partial\psi_{T_K(q)}$ its subdifferential and let $(t_0, q_0, v_0)$ be admissible initial data. We consider the following problem.

**Problem (P') (Moreau's sweeping process).** To find $\tau > 0$ and a function $u \in bvem([t_0, t_0 + \tau])$ with $u(t_0) = u(t_0 + 0) = v_0$ such that, putting

$$q(t) = q_0 + \int_{t_0}^{t} u(s) \, ds$$
for all $t \in [t_0, t_0 + \tau]$ we have

$$ \forall t \in [t_0, t_0 + \tau] \quad u(t) \in T_K(q(t)) \qquad (4.1) $$

and

$$ f(t, q(t), M(q(t))u(t)) - M(q(t))\dot{u}(t) \in \partial\psi_{T_K(q(t))}(u(t)) $$

in the following sense: there exists a non-negative measure $\mu$ on $[t_0, t_0 + \tau]$ such that the Lebesgue's measure and the Stieltjes measure $\dot{u}$ are absolutely continuous with respect to $\mu$ and the respective densities $t'_\mu \in L^1([t_0, t_0 + \tau]; \mathbb{R})$ and $u'_\mu \in L^1([t_0, t_0 + \tau]; \mathbb{R}^d; \mathrm{d}\mu)$ satisfy the differential inclusion

$$ f(t, q(t), M(q(t))u(t))t'_{\mu}(t) - M(q(t))u'_{\mu}(t) \in \partial\psi_{T_K(q(t))}(u(t)) \quad \mathrm{d}\mu \text{ a.e.} \qquad (4.2) $$

In order to compare problems (P) and (P') let us come back to the basic description of the dynamics. Let $q$ be a solution of problem (P) and assume that (1.11) holds. Then equations (1.2)–(1.4) can be rewritten as

$$ M(q)\ddot{q} - f(t, q, p) = \nu = \lambda \mathrm{d}\varphi(q), \qquad (4.3) $$

where $\lambda$ is a non-negative measure satisfying

$$ \mathrm{supp}(\lambda) \subset \{t: q(t) \in \partial K\}. \qquad (4.4) $$

Moreover, whenever $q(t) \in \partial K$, (1.6) implies that

$$ \begin{aligned} \mathrm{d}\varphi(q(t))(\dot{q}(t+0) + e\dot{q}(t-0)) &= \langle p(t+0) + ep(t-0), \mathrm{d}\varphi(q(t)) \rangle_{q(t)}^* \\ &= \langle (1+e)p_T(t-0), \mathrm{d}\varphi(q(t)) \rangle_{q(t)}^* = 0 \end{aligned} $$

and we infer that

$$ \frac{\dot{q}(t+0) + e\dot{q}(t-0)}{1+e} \in \partial T_K(q(t)). \qquad (4.5) $$

But, for all $q \in \partial K$, $T_K(q)$ is a half-space and it follows that

$$ \forall q \in \partial K, \quad \forall z \in \partial T_K(q) \quad \partial\psi_{T_K(q)}(z) = \mathbb{R}^{-} \mathrm{d}\varphi(q). \qquad (4.6) $$

If $q$ belongs to $\mathrm{Int}(K)$, $T_K(q) = \mathbb{R}^d$ and we have

$$ \forall z \in \mathbb{R}^d \quad \partial\psi_{T_K(q)}(z) = \{0\}. \qquad (4.7) $$

Taking into account (4.5)–(4.7) we obtain that

$$ -M(q)\ddot{q} + f(t, q, p) \in \partial\psi_{T_K(q)}\left(\frac{\dot{q}(t+0) + e\dot{q}(t-0)}{1+e}\right). $$

Let us define $u$ by

$$ \left. \begin{array}{l} \forall t \in (t_0, t_0 + \tau) \quad u(t) = \frac{\dot{q}(t+0) + e\dot{q}(t-0)}{1+e}, \\ u(t_0) = \dot{q}(t_0+0), \quad u(t_0 + \tau) = \dot{q}(t_0 + \tau - 0). \end{array} \right\} \qquad (4.8) $$

Since $\dot{q}$ is a function of bounded variation, we have $u(t) = \dot{q}(t)$ except perhaps on a countable subset of $[t_0, t_0 + \tau]$ and thus $u \in \mathrm{bvem}([t_0, t_0 + \tau])$ and $\dot{u} = \ddot{q}$. It follows
that $u$ is a solution of problem (P'): the differential inclusion (4.2) is satisfied with the choice $d\mu = |du| + dt$ for instance.

Inversely, let $u$ be a solution of problem (P') and let $q$ be defined by

$$ \forall t \in [t_0, t_0 + \tau] \qquad q(t) = q_0 + \int_{t_0}^{t} u(s) \, ds. \qquad (4.9) $$

It is clear that $q$ is continuous and $\dot{q}$ is a function of bounded variation such that $\dot{q}(t) = u(t)$ except perhaps on a countable subset of $[t_0, t_0 + \tau]$. Moreover, (4.1) implies that

$$ d\varphi(q(t))u(t) \ge 0 \quad \text{if } q(t) \notin \text{Int}(K). $$

Since $\varphi(q_0) = \varphi(q(t_0)) \ge 0$ we infer that $\varphi(q(t)) \ge 0$, i.e. $q(t) \in K$ for all $t \in [t_0, t_0 + \tau]$. We define the measure $\nu$ by

$$ \nu = M(q)\ddot{q} - f(t, q, M(q)\dot{q}). $$

Relation (4.2) implies that $\nu$ admits the following density with respect to $\mu$

$$ \nu'_{\mu}(t) = M(q(t))u'_{\mu}(t) - f(t, q(t), M(q(t)))t'_{\mu}(t) \in -\partial\psi_{T_K(q(t))}(u(t)) \quad d\mu \text{ a.e.} $$

Then, using (4.6) and (4.7) we obtain

$$ \text{supp}(\nu) \subset \{t: q(t) \in \partial K\} $$

and we infer that there exists a non-negative measure $\lambda$ such that

$$ \nu = \lambda d\varphi(q). $$

Furthermore, we have the following lemma.

**Lemma 4.2.** Let $t \in (t_0, t_0 + \tau)$ be a discontinuity point of $\dot{q}$. Then $q(t) \in \partial K$ and the impact law (1.6) is satisfied.

*Proof.* Let $t \in (t_0, t_0 + \tau)$ be such that $\dot{q}(t-0) \ne \dot{q}(t+0)$. Then (4.9) implies that $u(t-0) = \dot{q}(t-0)$, $u(t+0) = \dot{q}(t+0)$ and $\dot{u}$ has an atom at $t$. Since $\dot{u}$ is absolutely continuous with respect to $\mu$ we deduce that $\mu$ has also an atom at $t$ and

$$ 0 \ne u(t+0) - u(t-0) = u'_{\mu}(t)\mu(\{t\}) \quad \text{with } \mu(\{t\}) > 0. $$

From relation (4.2) we infer that

$$ 0 \ne M(q(t)) \frac{u(t-0) - u(t+0)}{\mu(\{t\})} \in \partial\psi_{T_K(q(t))}(u(t)), $$

which implies that $T_K(q(t)) \ne \mathbb{R}^d$ and $u(t) \in \partial T_K(q(t))$. It follows that $q(t)$ belongs to $\partial K$ and with (4.6) we obtain

$$ p(t-0) - p(t+0) \in \mathbb{R}^{-} d\varphi(q(t)), \\ d\varphi(q(t))u(t) = d\varphi(q(t))M(q(t))^{-1} \left( \frac{p(t+0) + ep(t-0)}{1+e} \right) = 0, $$

and (1.6) follows.
In order to conclude let us remark that, for all $t \in (t_0, t_0 + \tau)$ such that $q(t) \in \partial K$ with $\dot{q}(t+0) = \dot{q}(t-0)$, we have

$$ \dot{q}(t+0) = \dot{q}(t-0) \in -T_K(q(t)) \cap T_K(q(t)) = \partial T_K(q(t)) $$

and thus

$$ \begin{aligned} p_N(t \pm 0) &= \frac{\langle p(t \pm 0), \mathrm{d}\varphi(q(t)) \rangle_{q(t)}^*}{\langle \mathrm{d}\varphi(q(t)), \mathrm{d}\varphi(q(t)) \rangle_{q(t)}^*} \mathrm{d}\varphi(q(t)) \\ &= \frac{\mathrm{d}\varphi(q(t))\dot{q}(t \pm 0)}{\langle \mathrm{d}\varphi(q(t)), \mathrm{d}\varphi(q(t)) \rangle_{q(t)}^*} \mathrm{d}\varphi(q(t)) = 0. \end{aligned} $$

Hence the impact law (1.6) holds for all $t$ such that $q(t) \in \partial K$. Consequently, we have proved the following proposition.

**Proposition 4.3.** Let us assume (1.11). For all admissible initial data $(t_0, q_0, v_0)$ the following properties are equivalent:

(i) $(\tau, q)$ is a solution of problem (P),

(ii) $(\tau, u)$ is a solution of problem (P'),

with $u$ and $q$ given by (4.8) and (4.9).

### (b) Monteiro Marques and Mabrouk time discretization

On the contrary of Paoli and Schatzman scheme, this time discretization is written in terms of velocities. For given $(t_0, q_0, v_0)$ the algorithm is defined by

$$ q^0 = q_0, \quad hf(t_0, q_0, M(q_0)v^0) - M(q_0)(v^0 - v_0) \in \partial\psi_{T_K(q_0)}\left(\frac{v^0 + ev_0}{1+e}\right) \quad (4.10) $$

and, for all $i \ge 0$,

$$ q^{i+1} = q^i + hv^i, \quad (4.11) $$

$$ \begin{aligned} hf(t_0 + (i+1)h, q^{i+1}, M(q^{i+1})v^{i+1}) - M(q^{i+1})(v^{i+1} - v^i) \\ \in \partial\psi_{T_K(q^{i+1})}\left(\frac{v^{i+1} + ev^i}{1+e}\right). \end{aligned} \quad (4.12) $$

If we denote $q_h$ the linear interpolation of the $q^i$, then $v^i$ and $v^{i+1}$ are respectively the right and left velocities at $t_{i+1} = t_0 + (i+1)h$ and (4.12) is a very natural discretization of the differential inclusion (4.2). Moreover, we could remark that, whenever $q^i$ belongs to $\mathrm{Int}(K)$, the algorithm reduces to a centred scheme of order 2 for the ODE $M(q)\ddot{q} = f(t, q, M(q)\dot{q})$. More precisely, if $q^i \in \mathrm{Int}(K)$, we have $T_K(q^i) = \mathbb{R}^d$, thus (4.12) implies that

$$ v^i = v^{i-1} + hM(q^i)^{-1}f(t_i, q^i, M(q^i)v^i) $$

and using (4.11) we get

$$ \frac{q^{i+1} - 2q^i + q^{i-1}}{h^2} = M(q^i)^{-1}f\left(t_i, q^i, M(q^i)\frac{q^{i+1} - q^i}{h}\right). \quad (4.13) $$
**Remark 4.4.** The choice

$$F(t, u_1, u_2, v, h) = M(u_1)^{-1}f(t, u_1, M(u_1)(2v + (u_2 - u_1)/h))$$

in (3.3) yields also (4.13) whenever $(q^{i+1} + eq^{i-1})/(1 + e)$ belongs to $\text{Int}(K)$.

Since, for all $q \in \mathbb{R}^d$, $T_K(q)$ is a closed, non-empty, convex set we can rewrite (4.10)–(4.12) as follows:

$$\left. \begin{aligned} q^0 &= q_0, \\ v^0 &= -ev_0 + (1+e) \operatorname{Proj}_{q_0} \left( T_K(q_0), v_0 + \frac{h}{1+e} M(q_0)^{-1} f(t_0, q_0, M(q_0)v^0) \right) \end{aligned} \right\} \quad (4.14)$$

and, for all $i \ge 0$,

$$\left. \begin{aligned} q^{i+1} &= q^i + hv^i, \\ v^{i+1} &= -ev^i + (1+e) \operatorname{Proj}_{q^{i+1}} \left( T_K(q^{i+1}), v^i + \frac{h}{1+e} M(q^{i+1})^{-1} f(t_{i+1}, q^{i+1}, M(q^{i+1})v^{i+1}) \right). \end{aligned} \right\} \quad (4.15)$$

**Remark 4.5.** We recall that $\operatorname{Proj}_q(T_K(q), \cdot)$ denotes the projection on $T_K(q)$ for the kinetic metric frozen at $q$, i.e. for the distance defined by

$$\forall (z, z') \in (\mathbb{R}^d)^2 \quad \operatorname{dist}(z, z') = [^t(z - z')M(q)(z - z')]^{1/2}.$$

In order to understand how the constraint on positions and the impact law are taken into account by this time discretization, let us consider once again the very simple test-problem introduced in §3a: $d = 1$, $K = \mathbb{R}^+$, $M(q) = 1$, $f = 0$, $(t_0, q_0, v_0) = (0, 1, -1)$. We obtain

$$q^0 = 1, \quad v^0 = -1,$$

and, for all $i \ge 0$,

$$q^{i+1} = q^i + hv^i, \quad v^{i+1} = -ev^i + (1+e) \operatorname{Proj}_{q^{i+1}} (T_K(q^{i+1}), v^i),$$

with

$$\forall z \in \mathbb{R}^d \quad \operatorname{Proj}_{q^{i+1}}(T_K(q^{i+1}), z) = \begin{cases} z & \text{if } q^{i+1} > 0, \\ \max(z, 0) & \text{if } q^{i+1} \le 0. \end{cases}$$

Let us assume that $h \in (0, 1)$ and define $p = \max\{i : q^k > 0 \ \forall k \in \{0, \dots, i\}\}$. We have $p \ge 1$ and

$$\forall i \in \{0, \dots, p\} \quad q^{i+1} = 1 - (i+1)h, \quad v^i = -1.$$

By definition of $p$ we have $q^{p+1} \le 0$ and thus

$$v^{p+1} = -ev^p + (1+e) \max(v^p, 0) = e.$$

By an immediate induction we infer that

$$\forall n \ge 1 \quad q^{p+1+n} = 1 - (p+1)h + neh, \quad v^{p+n} = e.$$
<figure><img src="image_7.png" /><figcaption>Figure 6. Comparison between the exact solution and the time discretization algorithms: the bouncing ball test-problem with $e = 0.3$.</figcaption></figure>

The velocity is reversed and multiplied by $e$ instantaneously at $t_{p+1} \in [1, 1+h]$, but the number of approximate positions which do not satisfy the constraint depends on $e$. The results obtained with this algorithm are compared with those obtained with the scheme (3.1)–(3.3) in figure 6.

**Remark 4.6.** With the scheme (3.1)–(3.3) the discrete velocity is not reversed instantaneously (two time-steps are needed) but at most two approximate positions are outside $K$. Although (3.1)–(3.3) discretize directly the constraints by means of the average position $\bar{q}^i = (q^{i+1} + eq^{i-1})/(1+e)$, the scheme (4.14), (4.15) discretizes the constraints at the velocity level since

$$ \forall i \ge 0 \quad \frac{v^{i+1} + ev^i}{1+e} \in T_K(q^{i+1}), $$

which implies

$$ d\varphi(q^{i+1}) \left( \frac{v^{i+1} + ev^i}{1+e} \right) \ge 0 \quad \text{if } q^{i+1} \notin \text{Int}(K). $$

### (c) Convergence result

The convergence of the scheme (4.14), (4.15) is proved in Mabrouk (1998b) under the following assumptions:

(H'1) $\varphi \in C^{1,\beta}$ with $\beta > 1/2$;

(H'2) the kinetic metric can be locally identified to the usual Euclidean metric of $\mathbb{R}^d$;

(H'3) the function $f$ is globally continuous and uniformly Lipschitzian in its third variable.
The approximate velocity $v_h$ is defined by

$$ \forall t \in [t_0 + ih, t_0 + (i + 1)h) \quad v_h(t) = v^i. $$

We have the following theorem.

**Theorem 4.7 (Mabrouk 1998b).** Assume $(H'1)-(H'3)$. There exists $\tau > 0$ such that

(i) the sequence $(q_h, v_h)$ possesses a subsequence which converges uniformly $\times$ pointwise on $[t_0, t_0 + \tau]$ to a limit $(q, v)$, where $v$ is a function of bounded variation and

$$ q(t) = q_0 + \int_{t_0}^{t} v(s) \, ds \quad \text{for all } t \in [t_0, t_0 + \tau], $$

(ii) the couple $(\tau, u)$ with

$$ \begin{array}{ll} \forall t \in (t_0, t_0 + \tau) & u(t) = \frac{v(t+0) + ev(t-0)}{1+e}, \\ u(t_0) = v(t_0+0), & u(t_0 + \tau) = v(t_0 + \tau - 0), \end{array} $$

is a solution of problem (P').

This is a generalization of the convergence result obtained by Monteiro Marques in the case of inelastic shocks ($e = 0$) and a bounded function $f$ depending only on $t$ and $q$ (Monteiro Marques 1985, 1993).

The sketch of the proof is the same as in Paoli & Schatzman (1999b). First, uniform estimates for the discrete velocities and accelerations are obtained and the Ascoli-Arzela theorem and Helly's theorem imply (i). Then the differential inclusion (4.2) is established with $d\mu = |du| + dt$. Contrary to Paoli & Schatzman (1999b), general properties of convex analysis play an essential role throughout the proof, in particular, the two following lemmas.

**Lemma 4.8 ('lemma of the two cones' (Moreau 1962)).** Let $V$ and $N$ be a pair of mutually polar closed convex cones in a real Hilbert space $E$. Then, for all $(x, y, z) \in E^3$ we have

$$ \begin{cases} x = \operatorname{Proj}(V, z), \\ y = \operatorname{Proj}(N, z), \end{cases} \quad \text{if and only if} \quad \begin{cases} z = x + y, \\ x \in V, y \in N, x \cdot y = 0. \end{cases} $$

**Lemma 4.9 (Moreau 1978).** Let $C$ be a closed convex cone of a Hilbert space $E$, $C$ containing a ball $\bar{B}(a, r)$. Then

$$ \forall x \in E \quad |x - \operatorname{Proj}(C, x)| \le \frac{1}{2r} (|x - a|^2 - |\operatorname{Proj}(C, x) - a|^2). $$

Note that the mapping $q \mapsto T_K(q)$ is lower semi-continuous, which ensures that there exists $\delta > 0$ and a fixed ball $\bar{B}(a, r)$ such that

$$ \forall q \in \bar{B}(q_0, \delta) \quad \bar{B}(a, r) \subset T_K(q). $$

For a more detailed presentation of this proof, see the review paper of Kunze & Monteiro Marques (2000).
5. The multi-constraint case

Let us assume now that assumption (1.8) does not hold any more, i.e.

$$
\exists t \text{ s.t. } \mathrm{Card} \, J(q(t)) \ge 2. \tag{5.1}
$$

As we have already explained at § 1, the decomposition of the momenta $p(t \pm 0)$ as
$p(t \pm 0) = p_N(t \pm 0) + p_T(t \pm 0)$ with

$$
p_N(t \pm 0) \in \operatorname{Span}\{\mathrm{d}\varphi_\alpha(q(t)), \alpha \in J(q(t))\}, \\
\mathrm{d}\varphi_\alpha(q(t))(M(q(t))^{-1}p_T(t \pm 0)) = 0 \quad \forall \alpha \in J(q(t)),
$$

does not necessarily yield a simple relation between $p(t+0)$ and $p_N(t-0)$, $p_T(t-0)$ when (5.1) holds. In order to generalize easily the impact law (1.6) we rewrite it in terms of velocities. When $J(q(t)) = \{\alpha\}$ the impact law is given by

$$
p(t + 0) = p(t - 0) - (1 + e)p_N(t - 0)
$$

with

$$
p_N(t-0) = \frac{\langle p(t-0), \mathrm{d}\varphi_\alpha(q(t)) \rangle_{q(t)}^*}{\langle \mathrm{d}\varphi_\alpha(q(t)), \mathrm{d}\varphi_\alpha(q(t)) \rangle_{q(t)}^*} \mathrm{d}\varphi_\alpha(q(t)),
$$

which is equivalent to

$$
\dot{q}(t+0) = \dot{q}(t-0) - (1+e) \frac{\mathrm{d}\varphi_{\alpha}(q(t))\dot{q}(t-0)}{|\mathrm{d}\varphi_{\alpha}(q(t))|_{q(t)}^{*2}} M(q(t))^{-1} \mathrm{d}\varphi_{\alpha}(q(t)). \quad (5.2)
$$

For all $q \in \mathbb{R}^d$, let us denote by $N(q)$ the polar cone of $T_K(q)$ for the kinetic metric,
i.e.

$$
N(q) = \{w \in \mathbb{R}^d : \langle w, z \rangle_q = {}^t w M(q) z \le 0 \ \forall z \in T_K(q)\}.
$$

Then relation (5.2) can be rewritten as follows,

$$
\dot{q}(t+0) = \dot{q}(t-0) - (1+e) \mathrm{Proj}_{q(t)}(N(q(t)), \dot{q}(t-0)), \quad (5.3)
$$

where $\mathrm{Proj}_q(N(q), z)$ is the projection of $z$ on $N(q)$ relative to the kinetic metric
frozen at $q$. Whenever $\mathrm{Card}\, J(q(t)) = 1$, the impact law (1.6) is equivalent to (5.3)
and we assume that (5.3) still holds if $\mathrm{Card}\, J(q(t)) \ge 2$.

We should note that (5.3) is energetically consistent: if $q(t) \in \partial K$ let us denote

$$
\dot{q}_{\text{T}}^{-} = \text{Proj}_{q(t)}(T_{K}(q(t)), \dot{q}(t-0)), \quad \dot{q}_{\text{N}}^{-} = \text{Proj}_{q(t)}(N(q(t)), \dot{q}(t-0)).
$$

By using lemma 4.8 we get

$$
\begin{align*}
2T(t-0) = {}^t\dot{q}(t-0)M(q(t))\dot{q}(t-0) &= |\dot{q}_T^-|_{q(t)}^2 + |\dot{q}_N^-|_{q(t)}^2 \\
&\geq |\dot{q}_T^-|_{q(t)}^2 + e^2|\dot{q}_N^-|_{q(t)}^2 = 2T(t+0).
\end{align*}
$$

Moreover, if $q(t) \in \partial K$, we have

$$
\nu(\{t\}) = M(q(t))(\dot{q}(t+0) - \dot{q}(t-0)) = -(1+e)M(q(t))\dot{q}_N^{-} \\
\in \left\{ \sum_{\alpha \in J(q(t))} \lambda_{\alpha} \, \mathrm{d}\varphi_{\alpha}(q(t)), \lambda_{\alpha} \ge 0 \right\}
$$

and relation (1.4) still holds.

Thus, for given initial data $(t_0, q_0, v_0)$, the problem (P) is generalized as follows.
**Problem (P1).** To find $\tau > 0$ and a function $q$ from $[t_0, t_0 + \tau]$ to $\mathbb{R}^d$ satisfying the variational properties (1.9), the initial data (1.10) and such that (1.2)–(1.4) and (5.3) hold on $(t_0, t_0 + \tau)$.

Assuming that all the data are analytic, Ballard has performed a complete theoretical study of problem (P1) (Ballard 1999): an existence and uniqueness theorem is proved, including also the more general case of a restitution coefficient depending on $t$ and $q$.

Furthermore, Moreau’s formulation encompass the multi-constraint case and we can prove that problem (P1) is equivalent to problem (P’). This result relies on the same arguments as proposition 4.3 completed with the following.

**Lemma 5.1.** Let $t \in (t_0, t_0 + \tau)$ such that $q(t) \in \partial K$. Then (5.3) is equivalent to

$$M(q(t))(\dot{q}(t-0) - \dot{q}(t+0)) = -\nu(\{t\}) \in \partial\psi_{T_K(q(t))}\left(\frac{\dot{q}(t+0) + e\dot{q}(t-0)}{1+e}\right).$$

*Proof.* We just have to remark that (5.3) is equivalent to

$$\dot{q}_{\mathrm{N}}^{-} = \frac{\dot{q}(t-0) - \dot{q}(t+0)}{1+e}, \quad \dot{q}_{\mathrm{T}}^{-} = \frac{\dot{q}(t+0) + e\dot{q}(t-0)}{1+e},$$

and that $N(q(t))$ and $T_K(q(t))$ are mutually polar convex cones for the kinetic metric at $q(t)$. ■

As a consequence, the time discretizations of frictionless sweeping process allow to treat the multi-constraint case too (Moreau 1988, 1999; Jean & Moreau 1992) but the convergence is not proved: the assumption $(H'1)$ on the regularity of $\partial K$ plays an essential role in the proof of theorem 4.7 and it seems difficult to weaken it.

Let us point out that equations (1.2)–(1.4) and the decay of kinetic energy at impacts do not necessarily yield the impact law (5.3). Consider for instance the following example: $d = 2$, $K = \{(x, y) \in \mathbb{R}^2 : x \le 0, y \ge 0\}$, $M(q) = \mathrm{Id}_{\mathbb{R}^2}$, $f = 0$, $q_0 = (-1, 0)$ and $v_0 = (1, 0)$. The motions given by

$$q_1(t) = \begin{cases} (-1+t, 0) & \forall t \in [0, 1] \\ (0, t-1) & \forall t \ge 1 \end{cases} \quad \text{and} \quad q_2(t) = \begin{cases} (-1+t, 0) & \forall t \in [0, 1] \\ (1-t, 0) & \forall t \ge 1 \end{cases}$$

are two energy-conserving solutions of (1.2)–(1.4). But only $q_2$ satisfies (5.3) (with $e = 1$ of course).

Consequently, when the kinetic energy is conserved at impacts we should replace (5.3) by the more general impact law

$$|\dot{q}(t+0)| = |\dot{q}(t-0)| \quad \text{if } q(t) \in \partial K.$$

An existence result has been proved in Paoli (2000) by using a penalty method under the following assumptions: $K$ is a convex set and the mass matrix is trivial. Although the time discretization (3.1)–(3.3) still makes sense, its convergence does not seem easy to prove due to the lack of regularity of $\partial K$.

Moreover, in the multi-constraint case, a new difficulty occurs in the computation of approximate solutions: the motion is not continuous with respect to initial data in general. Even when the convergence of approximate solutions is proved the round-up errors may lead to a kind of unpredictibility. We refer the reader to Paoli (2000) for an illustration of this bad property.
# 6. Conclusion

In this article we have presented the time discretizations of vibro-impact problems derived from the basic description of the dynamics and from Moreau's sweeping process. We have recalled the associated convergence results and a qualitative comparison has been performed on a very simple test-problem. But nothing is said about the convergence order of these methods, which could seem surprising! Two reasons can be recalled. Firstly, the study of the test-problem shows that all the methods presented here are at most of order 1 but the determination of the convergence order in the general case seems very difficult and there is no theoretical result. Secondly, the dynamics of vibro-impact problems is often complex: sensitivity to initial data and chaotic behaviour may occur even for low-dimensional systems (see, for example, Guckenheimer & Holmes 1983; Shaw & Rand 1989). In this setting the convergence order is not an essential point since any precision in the computation of individual motions is lost in finite time. On the contrary, it is crucial to check that the mean values and the dynamical invariants (like attractors) are correctly approximated by the schemes. Consequently, the comparison of numerical methods for vibro-impact problems (event-driven scheme, compliant methods as well as time discretization algorithms) would require suitable benchmarks corresponding to various dynamical situations.

# References

Ballard, P. 1999 Dynamique des systèmes mécaniques discrets avec liaisons unilatérales parfaites. *C. R. Acad. Sci. Paris Sér. IIb* **327**, 953-958.

Bressan, A. 1959 Questioni di regolarita e di unicita del moto in presenza di vincoli olonomi unilaterali. *Rend. Sem. Mat. Univ. Padova* **29**, 271-315.

Brogliato, B. 1999 *Non-smooth mechanics*, 2nd edn. Springer Communications and Control Engineering Series.

Buttazzo, G. & Percivale, D. 1981 The bounce problem on *n*-dimensional Riemannian manifolds. *Atti Accad. Naz. Lincei Cl. Sci. Fis. Mat. Natur.* **70**, 246-250.

Buttazzo, G. & Percivale, D. 1983 On the approximation of the elastic bounce problem on Riemannian manifolds. *J. Diff. Eqns* **47**, 227-245.

Carriero, M. & Pascali, E. 1982 Uniqueness of the one-dimensional bounce problem as a generic property in $L^1([0, T]; \mathbb{R})$. *Boll. UMI* **5**, 87-91.

Carriero, M. & Pascali, E. 1985 Uniqueness for the elastic one-dimensional bounce problem. In *Nonlinear variational problems* (ed. A. Marino, L. Modica, S. Spagnolo & M. Degiovanni) pp. 88-91. (Research Notes in Mathematics, no. 127.) Pitman.

Delassus, E. 1917 Mémoire sur la théorie des liaisons finies unilatérales. Memorandum no. 34, Ecole Normale Supérieure, France.

Guckenheimer, J. & Holmes, P. 1983 *Nonlinear oscillations, dynamical systems and bifurcations of vector fields*. Springer.

Jacquart, G. & Gay, N. 1992 Computation of impact-friction interaction between a vibrating tube and its loose supports. In *Proc. ASME 92 Annual Winter Meeting, Int. Symp. on Flow Induced Vibration and Noise*. EDF Report no HP-61/92.133.

Jean, M. & Moreau, J. J. 1992 Unilaterality and dry friction in the dynamics of rigid body collections. In *Proc. Contact Mechanics Int. Symp., Lausanne, Switzerland, 7-9 October 1992*, pp. 31-48.

Kunze, M. & Monteiro Marques, M. D. P. 2000 An introduction to Moreau's sweeping process. In *Impacts in mechanical systems, analysis and modeling* (ed. B. Brogliato), pp. 1-60. Springer Lecture Notes in Physics.
Landau, L. D. & Lifschitz, E. M. 1959 *Theory of elasticity*. Pergamon.

Mabrouk, M. 1998a Liaisons unilatérales et chocs élastiques quelconques: un résultat d'existence. *C. R. Acad. Sci. Paris Sér. I* **326**, 1353-1357.

Mabrouk, M. 1998b A unified variational model for the dynamics of perfect unilateral constraints. *Eur. J. Mech. A* **17**, 819-842.

Monteiro Marques, M. P. D. 1985 Chocs inélastiques standards: un résultat d'existence. Semi-naire d'Analyse Convexe 15, exp. 4, Université des Sciences et Techniques du Languedoc.

Monteiro Marques, M. P. D. 1993 *Differential inclusions in non-smooth mechanical problems: shocks and dry friction*. In *Progress in nonlinear differential equations*, vol. 9. Birkhäuser.

Moreau, J. J. 1962 Décomposition orthogonale d'un espace hilbertien selon deux cônes mutuelle-ment polaires. *C. R. Séanc. Hebd. Acad. Sci. Paris* **255**, 238-240.

Moreau, J. J. 1963 Les liaisons unilatérales et le principe de Gauss. *C. R. Séanc. Hebd. Acad. Sci. Paris Séri. A* **256**, 871-874.

Moreau, J. J. 1978 Un cas de convergence des itérés d'une contraction d'un espace hilbertien. *C. R. Hebd. Séanc. Acad. Sci. Paris* **286**, 143-144.

Moreau, J. J. 1983 Liaisons unilatérales sans frottement et chocs inélastiques. *C. R. Séanc. Acad. Sci. Paris Sér. II* **296**, 1473-1476.

Moreau, J. J. 1985 Standard inelastic shocks and the dynamics of unilateral constraints. In *Unilateral problems in structural analysis* (ed. G. Del Piero & F. Maceri), pp. 173-221. Springer CISM Courses and Lectures, no. 288.

Moreau, J. J. 1986 Dynamique des systèmes à liaisons unilatérales avec frottement sec éventuel; essais numériques. Technical note no. 85-1, LMGC, Université des Sciences et Techniques du Languedoc.

Moreau, J. J. 1988 Unilateral contact and dry friction in finite freedom dynamics. In *Nonsmooth mechanics and applications* (ed. J. J. Moreau & P. D. Panagiotopoulos), pp. 1-82. Springer CISM Courses and Lectures, no. 302.

Moreau, J. J. 1994 Some numerical methods in multibody dynamics: application to granular materials. *Eur. J. Mech. A* **13**, 93-114.

Moreau, J. J. 1995 Numerical experiments in granular dynamics: vibration-induced size segregation. In *Proc. Contact Mechanics Int. Symp. 2, Carry-le-Rouet, France, 19-23 September 1994*, pp. 347-358.

Moreau, J. J. 1999 Numerical aspects of sweeping process. *Comput. Meth. Appl. Mech. Engng* **177**, 329-349.

Panet, M., Paoli, L. & Schatzman, M. 1995 Theoretical and numerical study for a model of vibrations with unilateral constraints. In *Proc. Contact Mechanics Int. Symp. 2, Carry-le-Rouet, France, 19-23 September 1994*, pp. 457-464.

Paoli, L. 1993 Analyse numérique de vibrations avec contraintes unilatérales. PhD thesis, Université Lyon I.

Paoli, L. 2000 An existence result for vibrations with unilateral constraints: case of a non-smooth set of constraints. *Math. Models Meth. Appl. Sci.* **10**, 815-831.

Paoli, L. & Schatzman, M. 1993a Vibrations avec contraintes unilatérales et perte d'énergie aux impacts, en dimension finie. *C. R. Acad. Sci. Paris Sér. I* **317**, 97-101.

Paoli, L. & Schatzman, M. 1993b Schéma numérique pour un modèle de vibrations avec contraintes unilatérales et perte d'énergie aux impacts, en dimension finie. *C. R. Acad. Sci. Paris Sér. I* **317**, 211-215.

Paoli, L. & Schatzman, M. 1993c Mouvement à un nombre fini de degrés de liberté avec contraintes unilatérales: cas avec perte d'énergie. *Modél. Math. Analyse Num.* **27**, 673-717.

Paoli, L. & Schatzman, M. 1995 Dynamical impact problem with loss of energy: numerical study for large times. Application to a model of tight joints. In *Proc. Contact Mechanics* 95, Ferrara, Italy, 11-13 July 1995, pp. 233-240.
Paoli, L. & Schatzman, M. 1998 Theoretical and numerical study for vibrations with unilateral constraints: case of a non-convex set of constraints. In *Proc. 4th World Congress on Computational Mechanics, New Trends and Applications, Buenos Aires, Argentina, 29 June-2 July 1998*. CD ROM.

Paoli, L. & Schatzman, M. 1999a Approximation et existence en vibro-impact, *C. R. Acad. Sci. Paris Sér. I* **329**, 1103-1107.

Paoli, L. & Schatzman, M. 1999b A numerical scheme for impact problems. Preprint, UMR 5585 MAPLY, CNRS. (Available at http://maply.univ-lyon1.fr.)

Paoli, L. & Schatzman, M. 1999c Dynamics of an impacting bar: a model and numerics. In *Proc. Eur. Conf. on Computational Mechanics, Munich, Germany, 31 August-3 September 1999*. CD ROM.

Paoli, L. & Schatzman, M. 2000 Ill-posedness in vibro-impact and its numerical consequences. In *Proc. Eur. Congress on Computational Mechanics, Barcelona, Spain, 11-14 September 2000*. CD ROM.

Percivale, D. 1985 Uniqueness in the elastic bounce problem. *J. Diff. Eqns* **56**, 206-215.

Percivale, D. 1986 Bounce problem with weak hypotheses of regularity. *Ann. Mat. Pura Appl.* **4**, 259-274.

Percivale, D. 1991 Uniqueness in the elastic bounce problem. II. *J. Diff. Eqns* **90**, 304-315.

Ravn, P. 1998 A continuous analysis method for planar multibody systems with joint clearance. *Multibody Syst. Dynam.* **2.1**, 1-24.

Rockafellar, R. T. 1970 *Convex analysis*. Princeton University Press.

Schatzman, M. 1977 Le système différentiel $(d^2u/dt^2) + \partial\varphi(u) \ni f$ avec conditions initiales. *C. R. Acad. Sci. Paris Sér. I* **284**, 603-606.

Schatzman, M. 1978 A class of nonlinear differential equations of second order in time. *Nonlin. Analysis* **2**, 355-373.

Schatzman, M. 1998 Uniqueness and continuous dependence on data for one-dimensional impact problems. *Math. Comput. Modelling* **28**, 1-18.

Shaw, S. W. & Rand, R. H. 1989 The transition to chaos in simple mechanical systems. *Int. J. Nonlin. Mech.* **24**, 41-56.

Stewart, D. E. 1997 Existence of solutions to rigid body dynamics and the Painlevé paradoxes. *C. R. Acad. Sci. Paris Sér. I* **325**, 689-693.

Stewart, D. E. 1998 Convergence of a time-stepping scheme for rigid body dynamics and resolution of Painlevé's paradoxes. *Arch. Ration. Mech. Analysis* **145**, 215-260.

Stoianovici, D. & Hurmuzlu, Y. 1996 A critical study of the applicability of rigid body collision theory. *ASME J. Appl. Mech.* **63**, 307-316.