# Simultaneous enforcement of constraints at position and velocity levels in the nonsmooth generalized-α scheme

Olivier Brüls, Vincent Acary, Alberto Cardona

► To cite this version:

Olivier Brüls, Vincent Acary, Alberto Cardona. Simultaneous enforcement of constraints at position and velocity levels in the nonsmooth generalized-α scheme. Computer Methods in Applied Mechanics and Engineering, 2014, 281, pp.131-161. 10.1016/j.cma.2014.07.025. hal-01059823

HAL Id: hal-01059823

https://inria.hal.science/hal-01059823

Submitted on 30 Nov 2015

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
Simultaneous enforcement of constraints at position and
velocity levels in the nonsmooth generalized-α scheme

Olivier Brüls<sup>a,*</sup>, Vincent Acary<sup>b</sup>, Alberto Cardona<sup>c</sup>

<sup>a</sup>University of Liège, Department of Aerospace and Mechanical Engineering (LTAS), Chemin des Chevreuils, 1 (B52), 4000 Liège, Belgium

<sup>b</sup>INRIA Rhône-Alpes, Centre de recherche Grenoble, 655 avenue de l'Europe, Inovallée de Montbonnot, 38334 St Ismier Cedex, France

<sup>c</sup>Universidad Nacional del Litoral - Conicet, CIMEC, Colectora Ruta Nac 168 / Paraje El Pozo, 3000 Santa Fe, Argentina

**Abstract**

This paper presents a formalism for the transient simulation of nonsmooth dynamic mechanical systems composed of rigid and flexible bodies, kinematic joints and frictionless contact conditions. The proposed algorithm guarantees the exact satisfaction of the bilateral and unilateral constraints both at position and velocity levels. Thus, it significantly differs from penalty techniques since no penetration is allowed. The numerical scheme is obtained in two main steps. Firstly, a splitting method is used to isolate the contributions of impacts, which shall be integrated with only first-order accuracy, from smooth contributions which can be integrated using a higher order scheme. Secondly, following the idea of Gear, Gupta and Leimkuhler, the equations of motion are reformulated so that the bilateral and unilateral constraints appear both at position and velocity levels. After time discretization, the equations of motion involve two complementarity conditions and it can be solved at each time step using a monolithic semi-smooth Newton method. The numerical behaviour of the proposed method is studied and compared to other approaches for a number of numerical examples. It is shown that the formulation offers a unified and valid approach for the description of contact conditions between rigid bodies as well as between flexible bodies.

*Keywords:* nonsmooth contact dynamics, flexible multibody system, time integration, generalized-α method, time-stepping schemes, index reduction

**1. Introduction**

This paper studies numerical algorithms for the simulation of mechanical sys-
tems including rigid and flexible bodies, kinematic joints and frictionless contact

<footnote>*Corresponding author, phone: +32.4.366.91.84, fax: +32.366.92.17</footnote>

<footnote>Email addresses: o.bruls@ulg.ac.be (Olivier Brüls), vincent.acary@inria.fr
(Vincent Acary), acardona@intec.unl.edu.ar (Alberto Cardona)</footnote>
conditions. The condition of impenetrability of the bodies in contact is expressed
as a unilateral constraint, with the consequence that impacts and/or instanta-
neous changes in the velocities may arise in the dynamic response. The design of
a consistent and stable time integration scheme for such systems requires great
care.

If dynamic contacts are analyzed between elastic solids and structures, the
gap velocity is necessarily discontinuous otherwise a non-physical penetration
of the bodies would occur. Due to the flexibility of the bodies, the contact
pressure remain nevertheless finite and no impulsive force is observed in this
case. Standard schemes for mechanical systems, such as the Newmark, HHT
or generalized-α methods, are based on the assumption that displacements and
velocities evolve smoothly as a function of time and are thus unable to deal
consistently with velocity jumps. This motivated the development of energy-
consistent schemes or stabilization techniques [1, 2, 3] for the analysis of dynamic
contacts between flexible bodies.

When contacts between rigid bodies are studied, e.g., in multibody systems,
impulsive reaction forces can occur leading to an instantaneous change in the
linear and angular momenta of each body. In this case, the aforementioned
standard schemes completely fail since the numerical response may artificially
generate energy when a contact occurs, see, e.g., [4]. In order to analyze such
impact phenomena, nonsmooth time integration methods have been proposed in
the literature and can be classified into two main groups, namely, event-driven
schemes and time-stepping schemes. Event-driven schemes are based on an ac-
curate event detection and the time step is adapted such that the end of the
step coincides with an event. At this time instant, the event is solved with the
help of an impact law. Such schemes are accurate for the free flight smooth
motions and are especially suitable for small multi-body systems with a limited
number of events, but they become inefficient if frequent transitions occur in a
short time. Contrary to the event-driven schemes, time-stepping schemes do not
adapt their time step size on events but only on some accuracy requirements if
needed. Two main families of time-stepping schemes have been designed up to
now: the Schatzman-Paoli scheme [5, 6] is based on a central difference scheme
and the Moreau-Jean scheme [7, 8, 9] is based on a θ-method. Time-stepping
methods have been proven to be convergent and robust, and are extensively ap-
plied as the solution to nonsmooth system models. In contrast to event-driven
schemes, time-stepping schemes are expected to have order-one accuracy even in
the smooth part of the motion. Thus, the accuracy is less satisfactory unless a
very small step size is applied. However, they remain robust and efficient even
for a large number of events. Some other strategies rely on an implicit detec-
tion of the impact times [10] and thus share properties of both event-driven and
time-stepping schemes.

The algorithm proposed here for systems with rigid and/or flexible bodies is
inspired from the Moreau-Jean time-stepping strategy. A fundamental property
of the Moreau-Jean scheme is that the unilateral constraints are imposed at
velocity level. As recognized by many authors [1, 2, 3, 8, 9], this choice leads to interesting consistency and stability properties of a simulation algorithm for dynamic contact analysis. In the Moreau–Jean scheme, the unilateral constraints at position level are only used to support the decision to activate the constraint but they are not exactly satisfied. The consequence is that some penetration can be observed in the numerical solution, which may not be physically acceptable.

In order to prevent such penetration problems, this paper presents an algo-
rithm which enforces the constraint not only at velocity level, so as to inherit good
consistency and stability properties, but also at position level. For that purpose,
the Gear-Gupta-Leimkuhler (GGL) approach, which was initially developed for
the stabilization of index-3 differential-algebraic equations [11], is generalized to
systems with unilateral constraints following a similar idea as in [12, 13]. After
time discretization, the equations of motion involve one linear complementar-
ity condition at position level and a second one at velocity level. However, the
present work differs from [12, 13] by the choice of the Lagrange multipliers which
appear in these complementarity conditions, as well as by the definition of an im-
plicit algorithm which solves the dynamic equilibrium and the complementarity
conditions in a monolithic way. At each time step, the problem is solved using
a semi-smooth Newton process, which can also be interpreted as an active set
method as discussed in [14, 15]. In order to guarantee the exact satisfaction of
the unilateral constraints by the converged solution without any penetration, the
active set is updated at each Newton iteration.

Another property of the Moreau–Jean algorithm is that the complete system
is integrated in time using a method which is only first-order accurate. This
appears as a disadvantage in applications where the nonsmooth phenomena are
localized in some mechanical parts of the system, while other parts exhibit smooth
motion. For example, the dynamic response of wind turbines systems is charac-
terized by a coupling between smooth motions and vibrations of large structural
components, such as the blades and the tower, and nonsmooth dynamic phenom-
ena in the gearbox and in the transmission line. In this case, the description of
vibration of phenomena in flexible bodies with a first-order method would require
quite small time steps and would thus be inefficient.

In order to improve the accuracy in the smooth part of motion, Chen et
al [16] observed that some terms in the equations of motion of a nonsmooth
dynamic system, such as the elastic forces, are smooth and can be integrated
in the time domain using a higher-order scheme, e.g., the generalized-α method.
All impulsive terms are still treated using a fully implicit integration scheme to
ensure consistency. This means that different integration formulae are used for
the different contributions to the equations of motion. Globally, the order of con-
vergence of the method is still limited to 1. Nevertheless, the advantage of this
approach is that the numerical dissipation of the generalized-α method is signif-
icantly smaller than the numerical dissipation of the Euler implicit scheme, so
that the energy behavior is strongly improved, especially for mechanical systems
exhibiting both impacts and structural vibrations. In some sense, in this scheme,
the high numerical dissipation of the Euler implicit scheme is only acting locally at the contact region and during the contact time, and not on the full system during the whole trajectory. A similar idea is exploited here with two noticeable differences. Firstly, the smooth part of the reaction forces of bilateral constraints are now treated as smooth variables and integrated with a higher-order scheme, whereas in [16] the total reaction forces of bilateral constraints were integrated using the first-order scheme only. Secondly, the treatment of the unilateral and of the bilateral constraint both at velocity and position levels prevents from any drift-off phenomenon.

The paper is organized as follows. Section 2 formulates the equations of motion with unilateral and bilateral constraints as an equality of differential measures. Moreau's sweeping process is described in Section 3. Section 4 presents the splitting strategy which is used to isolate smooth terms in the equations of motion and the GGL formulation of the constraints is described in Section 5. The time discretization method is developed in Section 6 and its properties are studied in Section 7. A semi-smooth Newton algorithm which solves the discretized equations of motion at each time step is presented in Section 8. In order to perform some numerical comparisons, two algorithms from the literature, the Moreau-Jean $\theta$-scheme and the smooth generalized-$\alpha$ method are briefly described in Section 9. Numerical examples with rigid impacts, accumulation phenomena, closed contact situations, and flexible contact conditions are studied in Section 10. Finally, the paper ends with some conclusions in Section 11.

## 2. Equations of motion

After spatial semi-discretization, the equations of motion of a flexible multi-body system including bilateral and unilateral constraints can be expressed in the following form:

$$ \dot{\mathbf{q}}^{+} = \mathbf{v}^{+} \qquad (1a) $$

$$ \mathbf{M}(\mathbf{q}) \mathrm{d}\mathbf{v} - \mathbf{g}_{\mathbf{q}}^{T} \mathrm{d}\mathbf{i} = \mathbf{f}(\mathbf{q}, \mathbf{v}, t) \mathrm{d}t \qquad (1b) $$

$$ \mathbf{g}^{\bar{\mathcal{U}}}(\mathbf{q}) = \mathbf{0} \qquad (1c) $$

$$ \mathbf{0} \leq \mathbf{g}^{\mathcal{U}}(\mathbf{q}) \perp \mathrm{d}\mathbf{i}^{\mathcal{U}} \geq \mathbf{0} \qquad (1d) $$

where

* $\mathbf{q}$ is the vector of coordinates, e.g., absolute nodal coordinates;

* $\dot{\mathbf{q}}^{+}(t) = \lim_{\tau \to t, \tau > t} \dot{\mathbf{q}}(\tau)$ and $\mathbf{v}^{+}(t) = \lim_{\tau \to t, \tau > t} \mathbf{v}(\tau)$ are the right limits of the velocity, which are functions of bounded variations;

* $\mathbf{f} = \mathbf{f}^{\mathrm{ext}}(t) - \mathbf{f}^{\mathrm{damp}}(\mathbf{q}, \mathbf{v}) - \mathbf{f}^{\mathrm{int}}(\mathbf{q})$ collects the external, damping and internal forces;

* $\mathbf{M}$ is the mass matrix which may, in general, depend on the coordinates;
* $dv$ is the differential measure associated with the velocity $v$ assumed to be of bounded variations;

* $t$ is time, and $dt$ is the corresponding standard Lebesgue measure;

* $g$ is the combined set of bilateral and unilateral constraints, and $g_q(q)$ is the corresponding matrix of constraint gradients;

* $di$ is the impulse measure of the contact reactions and of the bilateral forces;

* $U$ denotes the set of indices of the unilateral constraints, $\overline{U}$ is its complementarity set, i.e., the set of bilateral constraints, $C = U \cup \overline{U}$ is the full set of constraints, and we have

$$ g = \begin{bmatrix} g^U \\ g^{\overline{U}} \end{bmatrix}, \quad di = \begin{bmatrix} di^U \\ di^{\overline{U}} \end{bmatrix} \quad (2) $$

For the sake of notation simplicity, the convention $v(t) = v^+(t)$ and $\dot{q}(t) = \dot{q}^+(t)$ shall be used in the remaining part of the paper.

Since the motion might be nonsmooth, jumps in the velocity and impacts are expected yielding the following decomposition of the measures, neglecting singular measures:

$$ dv = \dot{v} dt + \sum_{i} (v(t_i) - v^-(t_i)) \delta_{t_i} \quad (3) $$

$$ di = \lambda dt + \sum_{i} p_i \delta_{t_i} \quad (4) $$

where $\lambda$ is the vector of nonimpulsive Lagrange multipliers associated with the Lebesgue measurable constraint forces; $v^-(t) = \lim_{\tau \to t, \tau < t} v(\tau)$; $(v(t_i) - v^-(t_i))$ is the jump in velocity at the instant $t_i$, $\delta_{t_i}$ is the Dirac delta supported at $t_i$, and $p_i$ is the impulse producing the jump at the instant $t_i$.

Equation (1d) is formulated as a complementarity condition, i.e., $(g^U(q))^T di^U = 0$ with both $g^U(q)$ and $di^U$ being non-negative. To complete the model, an impact law is needed. The Newton impact law defines the normal velocity jump in case of an impact for the constraint $j \in U$ as $g_q^j v(t) = -e^j g_q^j v^-(t)$, where $e^j \in [0, 1]$ is the coefficient of restitution, and $v^-(t)$ and $v(t)$ denote the velocity before and after the impact, respectively. Therefore, the contact condition is expressed at velocity level as [9]

$$ \text{if } g^j(q(t_i)) \le 0 \text{ then } 0 \le g_q^j v(t_i) + e^j g_q^j v^-(t_i) \perp p_i^j \ge 0, \quad \forall j \in U \quad (5) $$

The present formalism is developed for the analysis of contact conditions between rigid or flexible bodies. For rigid bodies, the coefficient of restitution defines the amount of energy dissipated during an impact. For flexible bodies, the physical meaning of a coefficient of restitution is not clear. The spatial discretization of a
flexible body using the finite element method leads to a finite dimensional system with finite masses. An impact law with a coefficient of restitution is thus needed to describe contact conditions. In practice, for flexible bodies, a value $e^j = 0$ may be used so that the condition $g_q^j \mathbf{v}(t_i) = 0$ is imposed when the constraint is active.

Inserting Eqs. (3-4) in (1a-1d) yields the standard equations of motion for almost every time $t$:

$$ \mathbf{M}(\mathbf{q})\dot{\mathbf{v}} - \mathbf{g}_{\mathbf{q}}^{T}\boldsymbol{\lambda} = \mathbf{f}(\mathbf{q}, \mathbf{v}, t) \quad (6a) $$

$$ \mathbf{g}_{\mathbf{q}}^{\overline{\mathcal{U}}}(\mathbf{q}) = \mathbf{0} \quad (6b) $$

$$ \mathbf{0} \le \mathbf{g}^{\mathcal{U}}(\mathbf{q}) \perp \boldsymbol{\lambda}^{\mathcal{U}} \ge \mathbf{0} \quad (6c) $$

and the impact equations at time $t_i$

$$ \mathbf{M}(\mathbf{q})(\mathbf{v}(t_i) - \mathbf{v}^{-}(t_i)) - \mathbf{g}_{\mathbf{q}}^{T}\mathbf{p}_i = \mathbf{0} \quad (6d) $$

$$ \mathbf{g}_{\mathbf{q}}^{\overline{\mathcal{U}}} \mathbf{v}(t_i) = \mathbf{0} \quad (6e) $$

$$ \text{if } g^j(\mathbf{q}(t_i)) \le 0 \text{ then } 0 \le g_{\mathbf{q}}^j \mathbf{v}(t_i) + e^j g_{\mathbf{q}}^j \mathbf{v}^{-}(t_i) \perp p_i^j \ge \mathbf{0}, \quad \forall j \in \mathcal{U} \quad (6f) $$

### 3. Moreau's sweeping process

The Moreau-Jean method is based on a reformulation of the equations of motion with the constraints at velocity level. For unilateral constraints, the constraint at velocity level is provided by the Newton's impact law. This leads to the following set of equations

$$ \dot{\mathbf{q}} = \mathbf{v} \quad (7a) $$

$$ \mathbf{M}(\mathbf{q}) \mathrm{d}\mathbf{v} - \mathbf{g}_{\mathbf{q}}^{T} \mathrm{d}\mathbf{i} = \mathbf{f}(\mathbf{q}, \mathbf{v}, t) \mathrm{d}t \quad (7b) $$

$$ \mathbf{g}_{\mathbf{q}}^{\overline{\mathcal{U}}} \mathbf{v} = \mathbf{0} \quad (7c) $$

$$ \text{if } g^j(\mathbf{q}) \le 0 \text{ then } 0 \le g_{\mathbf{q}}^j \mathbf{v} + e^j g_{\mathbf{q}}^j \mathbf{v}^{-} \perp \mathrm{d}i^j \ge 0, \quad \forall j \in \mathcal{U} \quad (7d) $$

As a consequence of Moreau's viability Lemma [17], the exact solution of these equations also satisfies the complementarity condition at position level in Eq. (1d).

If all unilateral constraints are closed, the original system of equations (1) with the constraints at position level is an index-3 DAE, while the system (7) with the constraints at velocity level is an index-2 DAE. From a numerical point of view, the difficulties encountered when solving an index-2 DAE are less severe than when solving the original index-3 DAE [18]. Therefore, the reformulation of the equations of motion in (7) can be interpreted as an index reduction process for systems with unilateral constraints. A time discretization of Eq. (7) leads to a discrete system in which the constraints are imposed at velocity level. This is an important property to guarantee the stability and the consistency of the numerical solution for nonsmooth dynamic systems with unilateral constraints and impacts. The consequence is that the constraints at velocity level are exactly satisfied but,
due to the presence of round-off and truncation errors, the constraints at position level are only approximately satisfied by the numerical solution. This leads to a drift phenomenon which may be characterized by a progressive increase in the constraint violation at position level during the simulation. Starting from the formulation in Eq. (7), this paper aims at developing a consistent approach which guarantees the exact satisfaction of the constraints both at velocity and position levels, in order to completely eliminate any drift problem.

## 4. Splitting strategy

Following a similar approach as in [16], a splitting of the variables and of the equations of motion is proposed to isolate impulsive terms, which are integrated with only first order accuracy, from the other contributions, which are assumed to be sufficiently smooth to be integrated using a higher order method. In this paper, the splitting will also be exploited to derive the GGL formulation presented in the next section.

For a given time step $(t_n, t_{n+1}]$, the position $\tilde{q}(t)$, velocity $\tilde{v}(t)$ and Lagrange multiplier $\tilde{\lambda}(t)$ are defined as the solution of the following initial value problem coupled with the variables $q(t)$ and $v(t)$

$$ \dot{\tilde{q}} = \tilde{v} \qquad (8a) $$

$$ M(\tilde{q}) \dot{\tilde{v}} - g_{\tilde{q}}^T(\tilde{q}) \tilde{\lambda} = f(\tilde{q}, v, t) \qquad (8b) $$

$$ g_{\tilde{q}}^{\overline{\lambda}}(\tilde{q}) \tilde{v} = 0 \qquad (8c) $$

$$ \tilde{\lambda}^{\overline{\lambda}} = 0 \qquad (8d) $$

with the initial value $\tilde{v}(t_n) = v(t_n)$, $\tilde{q}(t_n) = q(t_n)$. If there is no unilateral constraint, one observes that $\tilde{v}$ exactly satisfies the equations of motion (6a) with the bilateral constraints at velocity level, so that $\tilde{v} = v$ and $\tilde{\lambda} = \lambda$. For unilaterally constrained systems, Eq. (8) thus represent the dynamics of the system if the unilateral constraints and the contact forces are partially ignored on the interval $(t_n, t_{n+1}]$. Indeed, the unilateral constraints still have an indirect influence on the smooth solution $\tilde{q}(t)$, $\tilde{v}(t)$ because Eq. (8) is still coupled with the variables $q(t)$, $v(t)$. In the literature, $\tilde{v}$ is sometimes called the free velocity, i.e., the velocity of the system when the reaction forces of the unilateral constraints are null [19]. If $f(q(t), v(t), t)$ is a sufficiently smooth function of time and if the matrix

$$ \begin{bmatrix} M & g_{\tilde{q}}^{\overline{\lambda},T} \\ g_{\tilde{q}}^{\overline{\lambda}} & 0 \end{bmatrix} \qquad (9) $$

is nonsingular, Eq. (8) is a well-defined and smooth system of index-2 differential-algebraic equations. Let us remark that the mass matrix itself can be singular. Then, the variable $dw$ is defined by the equation

$$ dw = dv - \dot{v} dt \qquad (10) $$
The motivation for this splitting is to express $dv$ as the sum of a smooth contribution which is represented by the variable $\dot{\tilde{v}}$ and of some impulsive contributions which are all collected by the variable $dw$. The variable $dw$ includes all discontinuous contributions to the velocity, which result from the various impacts that occur during the time step. Based on this idea, $\tilde{v}$ can be evaluated by time integration using a second-order method while, to ensure consistency, $dw$ is evaluated using the Moreau-Jean scheme and the Euler implicit method. However, second-order accuracy cannot always be guaranteed for $\tilde{v}$ since $f(q(t), v(t), t)$ can be a discontinuous function of time whenever a velocity jump occurs. An elimination of $dv$ and $\dot{\tilde{v}}$ from Eqs. (7b,8b,10) yields

$$ M(q) dw - g_q^T (di - \tilde{\lambda} dt) = 0 \quad (11) $$

In terms of these variables, the equations of motion become

$$ \dot{q} = v \quad (12a) $$

$$ dv = dw + \dot{v} dt \quad (12b) $$

$$ M(q) \dot{v} - g_q^{\bar{\lambda},T} \tilde{\lambda}^{\bar{\lambda}} = f(q, v, t) \quad (12c) $$

$$ g_q^{\bar{\lambda}} \dot{v} = 0 \quad (12d) $$

$$ \tilde{\lambda}^{\mu} = 0 \quad (12e) $$

$$ M(q) dw - g_q^T (di - \tilde{\lambda} dt) = 0 \quad (12f) $$

$$ g_q^{\bar{\lambda}} v = 0 \quad (12g) $$

$$ \text{if } g^j(q) \le 0 \text{ then } 0 \le g_q^j v + e^j g_q^j v^- \perp di^j \ge 0, \quad \forall j \in \mathcal{U} \quad (12h) $$

The set of 3 equations (12b, 12c, 12f) for the three variables $dv$, $\dot{\tilde{v}}$ and $dw$ is strictly equivalent to the single Eq. (7b) for the variable $dv$.

## 5. Gear-Gupta-Leimkuhler formulation

The equations of motion (12) involve the bilateral and unilateral constraints at velocity level only. As mentioned in Section 3, the exact solution would also satisfy the constraints at position level

$$ g^{\bar{\lambda}}(q) = 0 \quad (13a) $$

$$ g^{\mu}(q) \ge 0 \quad (13b) $$

However, a time integration algorithm based on Eq. (12) that enforces the constraints at velocity level only will not satisfy the constraints at position level because of the drift phenomenon. In this section, we propose to include the constraints at velocity and position levels in a unique set of equations following a method similar to that proposed by Gear, Gupta and Leimkuhler [11]. The advantage of this strategy is that the time discretization procedure will lead to an
algorithm which simultaneously satisfies the constraints at position and velocity levels.

The idea is to supplement Eq. (12) with the constraints at position level. Equation (12a) is multiplied by the mass matrix

$$M(q)\dot{q} = M(q)v \quad (14)$$

and an additional Lagrange multiplier $\mu$ is introduced in order to accommodate for this duplicated set of position and velocity constraints and to have a well-posed problem. Therefore, Eq. (14) is modified as

$$M(q) \dot{q} - g_q^T \mu = M(q) v \quad (15)$$

where it is expected that $\mu$ vanishes identically for the exact solution as a consequence of the redundancy between the constraints at position and velocity levels. As a consequence, the variables $\dot{q}$ and $v$ are not formally equivalent, the variable $\dot{q}$ being related with the position variable $q$ through the equation

$$q(t) = q(t_n) + \int_{(t_n, t]} \dot{q}(\tau) d\tau \quad (16)$$

Following [12], the equilibrium equations are then obtained after replacing Eq. (12) by

$$M(q) \dot{q} - g_q^T \mu = M(q) v \quad (17a)$$

$$g_q^{\overline{\mu}}(q) = 0 \quad (17b)$$

$$0 \le g_q^{\overline{\mu}}(q) \perp \mu^{\overline{\mu}} \ge 0 \quad (17c)$$

$$d\mathbf{v} = d\mathbf{w} + \dot{\mathbf{v}} dt \quad (17d)$$

$$M(q) \dot{\mathbf{v}} - g_q^{\overline{\mu},T} \tilde{\mathbf{\lambda}}^{\overline{\mu}} = \mathbf{f}(q, v, t) \quad (17e)$$

$$g_q^{\overline{\mu}} \tilde{\mathbf{v}} = 0 \quad (17f)$$

$$\tilde{\mathbf{\lambda}}^{\overline{\mu}} = 0 \quad (17g)$$

$$M(q) d\mathbf{w} - g_q^T (d\mathbf{i} - \tilde{\mathbf{\lambda}} dt) = 0 \quad (17h)$$

$$g_q^{\overline{\mu}} \mathbf{v} = 0 \quad (17i)$$

$$\text{if } g^j(q) \le 0 \text{ then } 0 \le g_q^j v + e^j g_q^j v^- \perp d\mathbf{i}^j \ge 0, \quad \forall j \in \mathcal{U} \quad (17j)$$

By construction, any solution of Eq. (12) also satisfies Eq. (17) with $\mu = 0$. Let us also observe that this formulation is still valid when the mass matrix is singular, provided that the matrix in Eq. (9) is nonsingular.

If all unilateral constraints are closed, the GGL form of the equations of motion is an index-2 DAE. In this special case, the formalism presented here boils down to the DAE formalism proposed by Gear, Gupta and Leimkuhler [11]. An important difference compared to the index-2 DAE in Eq. (7) is that the constraints are now imposed both at velocity and position levels.
## 6. Time discretization

By construction, the variables $\tilde{q}(t)$, $\tilde{v}(t)$ and $\tilde{\lambda}(t)$ only capture a sufficiently smooth part of the motion during the time step $(t_n, t_{n+1}]$ and they shall be computed using a second-order time integration scheme. In order to obtain the complete time integration scheme, discrete nonsmooth corrections are now defined so that the total motion satisfies a discrete form of the bilateral constraints in Eqs. (17b,17i) as well as a discrete form of the complementarity conditions in Eqs. (17c,17j), which represent the unilateral constraints.

### 6.1. Discrete nonsmooth variables

The velocity jump and position correction variables $\mathbf{W}(t_n, t)$ and $\mathbf{U}(t_n; t)$ are defined as

$$ \mathbf{W}(t_n; t) = \int_{(t_n, t]} \mathrm{d}\mathbf{w} \qquad (18a) $$

$$ \mathbf{U}(t_n; t) = \int_{t_n}^{t} (\dot{\mathbf{q}} - \tilde{\mathbf{v}}) \mathrm{d}t \qquad (18b) $$

Since $\mathbf{v}(t_n) = \tilde{\mathbf{v}}(t_n)$ and $\mathbf{q}(t_n) = \tilde{\mathbf{q}}(t_n)$, we get

$$ \mathbf{W}(t_n; t) = \mathbf{v}(t) - \tilde{\mathbf{v}}(t) \qquad (19a) $$

$$ \mathbf{U}(t_n; t) = \mathbf{q}(t) - \tilde{\mathbf{q}}(t) \qquad (19b) $$

as well as $\mathbf{W}(t_n; t_n) = \mathbf{0}$ and $\mathbf{U}(t_n; t_n) = \mathbf{0}$. The multipliers $\mathbf{\Lambda}(t_n; t)$ and $\mathbf{\nu}(t_n; t)$ are defined as

$$ \mathbf{\Lambda}(t_n; t) = \int_{(t_n, t]} (\mathrm{d}\mathbf{i} - \tilde{\mathbf{\lambda}}(\tau)) \mathrm{d}\tau \qquad (20a) $$

$$ \mathbf{\nu}(t_n; t) = \int_{t_n}^{t} (\mathbf{\mu}(\tau) + \mathbf{\Lambda}(t_n; \tau)) \mathrm{d}\tau \qquad (20b) $$

with $\mathbf{\Lambda}(t_n; t_n) = \mathbf{\nu}(t_n; t_n) = \mathbf{0}$.

The discrete approximation of Eqs. (17a) and (17b) is based on the following result.

**Theorem 1.** The velocity jump and position correction variables, which are caused by the unilateral constraints, satisfy

$$ \mathbf{M}(\mathbf{q}(t_{n+1})) \mathbf{W}(t_n; t_{n+1}) - \mathbf{g}_{\mathbf{q}}^{T}(\mathbf{q}(t_{n+1})) \mathbf{\Lambda}(t_n; t_{n+1}) = \mathcal{O}(h) \qquad (21a) $$

$$ \mathbf{M}(\mathbf{q}(t_{n+1})) \mathbf{U}(t_n; t_{n+1}) - \mathbf{g}_{\mathbf{q}}^{T}(\mathbf{q}(t_{n+1})) \mathbf{\nu}(t_n; t_{n+1}) = \mathcal{O}(h^2) \qquad (21b) $$

in which $h = t_{n+1} - t_n$ is the step size.
Proof. The proof is obtained by integrating Eqs. (17a) and (17h) over the time step $(t_n, t_{n+1}]$.

Firstly, we study the integrated form of Eq. (17h). Since $\mathbf{q}(t)$ is continuous, we have $\mathbf{M}(\mathbf{q}(\tau)) = \mathbf{M}(\mathbf{q}(t)) + \mathcal{O}(h)$ and $\mathbf{g}_{\mathbf{q}}^{T}(\mathbf{q}(\tau)) = \mathbf{g}_{\mathbf{q}}^{T}(\mathbf{q}(t)) + \mathcal{O}(h)$, $\forall \tau \in (t_n, t]$. Then, Eq. (21a) is obtained by observing that, $\forall t \in (t_n, t_{n+1}]$,

$$ 
\begin{aligned}
\mathbf{0} &= \int_{(t_n, t]} \left( \mathbf{M}(\mathbf{q})\mathrm{d}\mathbf{w} - \mathbf{g}_{\mathbf{q}}^{T}(\mathrm{d}\mathbf{i} - \tilde{\mathbf{\lambda}}\mathrm{d}\tau) \right) \\
&= \mathbf{M}(\mathbf{q}(t))\mathbf{W}(t_n; t) - \mathbf{g}_{\mathbf{q}}^{T}(\mathbf{q}(t))\mathbf{\Lambda}(t_n; t) + \mathcal{O}(h)
\end{aligned} 
\quad (22) $$

Secondly, we study the integrated form of Eq. (17a). Using $\mathbf{M}(\mathbf{q}(t)) = \mathbf{M}(\mathbf{q}(t_{n+1})) + \mathcal{O}(h)$, $\mathbf{g}_{\mathbf{q}}^{T}(\mathbf{q}(t)) = \mathbf{g}_{\mathbf{q}}^{T}(\mathbf{q}(t_{n+1})) + \mathcal{O}(h)$, $\forall t \in (t_n, t_{n+1}]$, and $\mathbf{q}_n = \tilde{\mathbf{q}}_n$, the integrated form of Eq. (17a) is obtained as

$$ 
\begin{aligned}
\mathbf{0} &= \int_{t_n}^{t_{n+1}} (\mathbf{M}(\mathbf{q}(t))(\dot{\mathbf{q}}(t) - \mathbf{v}(t)) - \mathbf{g}_{\mathbf{q}}^{T}(\mathbf{q}(t))\boldsymbol{\mu}(t)) \mathrm{d}t \\
&= \int_{t_n}^{t_{n+1}} (\mathbf{M}(\mathbf{q}(t))(\dot{\mathbf{q}}(t) - \tilde{\mathbf{v}}(t) - \mathbf{W}(t_n; t)) - \mathbf{g}_{\mathbf{q}}^{T}(\mathbf{q}(t))\boldsymbol{\mu}(t)) \mathrm{d}t \\
&= \int_{t_n}^{t_{n+1}} (\mathbf{M}(\mathbf{q}(t))(\dot{\mathbf{q}}(t) - \tilde{\mathbf{v}}(t)) - \mathbf{g}_{\mathbf{q}}^{T}(\mathbf{q}(t)) (\mathbf{\Lambda}(t_n; t) + \boldsymbol{\mu}(t)) + \mathcal{O}(h)) \mathrm{d}t \\
&= \mathbf{M}(\mathbf{q}(t_{n+1})) (\mathbf{q}(t_{n+1}) - \mathbf{q}(t_n) - (\tilde{\mathbf{q}}(t_{n+1}) - \tilde{\mathbf{q}}(t_n))) - \mathbf{g}_{\mathbf{q}}^{T}(\mathbf{q}(t_{n+1})) \boldsymbol{\nu}(t_n; t_{n+1}) + \mathcal{O}(h^2) \\
&= \mathbf{M}(\mathbf{q}(t_{n+1}))\mathbf{U}(t_n; t_{n+1}) - \mathbf{g}_{\mathbf{q}}^{T}(\mathbf{q}(t_{n+1})) \boldsymbol{\nu}(t_n; t_{n+1}) + \mathcal{O}(h^2)
\end{aligned}
$$

where Eqs. (19a), (22), (20b) and (19b) have been successively used. □

Let us remark that, if impacts occur and $h \to 0$, $\mathbf{W}(t_n; t_{n+1}) = \mathcal{O}(1)$ and $\mathbf{\Lambda}(t_n; t_{n+1}) = \mathcal{O}(1)$ since the velocities may exhibit finite jumps, whereas $\mathbf{U}(t_n; t_{n+1}) = \mathcal{O}(h)$ and $\boldsymbol{\nu}(t_n; t_{n+1}) = \mathcal{O}(h)$ since the positions remain continuous. The position corrections eliminate the drift of the constraints, and the velocity jumps are approximated over the time step.

## 6.2. Discrete complementarity conditions

The discrete complementarity conditions at velocity and position levels can be obtained from the inequalities $\mathrm{d}\mathbf{i}^{\mathcal{U}} \ge \mathbf{0}$ in Eq. (17j), and $\boldsymbol{\mu}^{\mathcal{U}} \ge \mathbf{0}$ in Eq. (17c). At velocity level, Eqs. (17g) and (20a) imply that

$$ \int_{(t_n, t_{n+1}]} \mathrm{d}\mathbf{i}^{\mathcal{U}} = \mathbf{\Lambda}^{\mathcal{U}}(t_n; t_{n+1}) \quad (23) $$

so that, following a similar idea as in the Moreau–Jean method, the complementarity condition at velocity level is expressed as

$$ \text{if } g^j(\mathbf{q}_{n+1}^*) \le 0 \text{ then } 0 \le g_{\mathbf{q},n+1}^j \mathbf{v}_{n+1} + e^j g_{\mathbf{q},n}^j \mathbf{v}_n \perp \mathbf{\Lambda}_{n+1}^j \ge 0, \quad \forall j \in \mathcal{U} \quad (24) $$
At position level, a fully implicit complementarity expression between the collocated value of the constraint $g^U(q)$ at time $t_{n+1}$ and the multiplier $\mu^U(t)$ at time $t_{n+1}$ would lead to

$$0 \le g^U(q_{n+1}) \perp \mu_{n+1}^U \ge 0 \qquad (25)$$

This relation is written in terms of the impulse $\mu^U$, which is not directly accessible in the algorithm (integrated values along the time step are instead used, i.e., $\nu^U(t_n; t_{n+1})$ and $\Lambda^U(t_n; t_{n+1})$).

Starting from Eq. (20b), we can get by using the $\theta$-method

$$\nu^U(t_n; t_{n+1}) = h \left( \theta \left( \mu^U(t_{n+1}) + \Lambda^U(t_n; t_{n+1}) \right) + (1-\theta) \mu^U(t_n) \right) + \mathcal{O}(h). \quad (26)$$

This result can be obtained from Lemma 1 in [20], assuming that $\mu^U(t)+\Lambda^U(t_n; t)$ is a bounded variation function. The discontinuous behaviour of $\Lambda^U(t_n; t)$ in the interval $t \in (t_n, t_{n+1}]$ makes it difficult to get a better numerical approximation, because the contact time instant is not known. Only in the special case where no impact occurs in the interval $(t_n, t_{n+1}]$, the variable $\Lambda^U(t_n; t)$ is continuous and the local truncation error in Eq. (26) drops to $\mathcal{O}(h^2)$.

For instance, a fully implicit Euler approximation with $\theta = 1$ leads to

$$h\mu^U(t_{n+1}) \approx \nu^U(t_n; t_{n+1}) - h\Lambda^U(t_n; t_{n+1}) \qquad (27)$$

However, both terms $\nu^U(t_n; t_{n+1})$ and $h\Lambda^U(t_n; t_{n+1})$ are $\mathcal{O}(h)$, so that the relative error on this estimate is at best $\mathcal{O}(1)$, i.e., it is dominated by numerical artifacts. Experiments using this approximation in Eq. (25) evidenced spurious chattering phenomena, as shown in Section 10.1. For this reason, an alternative approach is proposed.

Note that, if contact phenomena occur during the time step $(t_n; t_{n+1}]$, then $\mu^j(t_n; t_{n+1}) \ge 0$, $\Lambda^j(t_n; t_{n+1}) > 0$ and therefore $\nu^j(t_n; t_{n+1}) > 0$. Based on this observation, we make the assumption that the numerical solution should satisfy the complementarity condition

$$0 \le g^U(q_{n+1}) \perp \nu_{n+1}^U \ge 0 \qquad (28)$$

From Eq. (21b), the situation $\nu^j(t_n; t_{n+1}) > 0$ implies that $U(t_n, t_{n+1}) \ne 0$, i.e., that the trajectory departs from the smooth trajectory and contact phenomena have occurred needing a position correction, imposed by the constraint $g^j(q_{n+1}) = 0$. Conversely, if $g^j(q_{n+1}) > 0$, then there is no impact at time $t_{n+1}$, no contact interaction has occurred and the condition $\nu_{n+1}^j = 0$ needs to be satisfied. This assumption does not affect the order of accuracy of the algorithm since it specifies the influence of an $\mathcal{O}(h)$ correction to the position variable, while $\mathcal{O}(h)$ local errors are anyway accepted when impacts occur. However, this assumption is essential because it significantly affects the qualitative behaviour of the numerical solution. It will be shown below that it results in a robust and chattering-free numerical scheme.
### 6.3. Time stepping scheme

The complete time integration scheme is obtained by combining a first-order approximation of the nonsmooth variables with a one-step and second-order time integration scheme for the smooth variables. In this work, the generalized-α method is used for the smooth part, however, other DAE time integration schemes could also be considered. According to [21], the generalized-α method combines the equations of motion at time step n + 1 with some difference equations. This hybrid time integration scheme permitting to advance the solution at each step is formulated as follows

$$M(\mathbf{q}_{n+1})\dot{\mathbf{v}}_{n+1} - f(\mathbf{q}_{n+1}, \mathbf{v}_{n+1}, t_{n+1}) - \mathbf{g}_{\mathbf{q},n+1}^{\bar{\mathcal{U}},T} \tilde{\boldsymbol{\lambda}}_{n+1}^{\bar{\mathcal{U}}} = \mathbf{0} \quad (29a)$$

$$\mathbf{g}_{\mathbf{q},n+1}^{\bar{\mathcal{U}}} \tilde{\mathbf{v}}_{n+1} = \mathbf{0} \quad (29b)$$

$$M(\mathbf{q}_{n+1})\mathbf{U}_{n+1} - \mathbf{g}_{\mathbf{q},n+1}^{T} \boldsymbol{\nu}_{n+1} = \mathbf{0} \quad (29c)$$

$$\mathbf{g}^{\bar{\mathcal{U}}}(\mathbf{q}_{n+1}) = \mathbf{0} \quad (29d)$$

$$\mathbf{0} \leq \mathbf{g}^{\mathcal{U}}(\mathbf{q}_{n+1}) \perp \boldsymbol{\nu}_{n+1}^{\mathcal{U}} \geq \mathbf{0} \quad (29e)$$

$$M(\mathbf{q}_{n+1})\mathbf{W}_{n+1} - \mathbf{g}_{\mathbf{q},n+1}^{T} \mathbf{\Lambda}_{n+1} = \mathbf{0} \quad (29f)$$

$$\mathbf{g}_{\mathbf{q},n+1}^{\bar{\mathcal{U}}} \mathbf{v}_{n+1} = \mathbf{0} \quad (29g)$$

$$\text{if } g^j(\mathbf{q}_{n+1}^*) \le 0 \text{ then } 0 \le g_{\mathbf{q},n+1}^j \mathbf{v}_{n+1} + e^j g_{\mathbf{q},n}^j \mathbf{v}_n \perp \mathbf{\Lambda}_{n+1}^j \ge 0, \quad (29h) \\ \forall j \in \mathcal{U}$$

with the notation $\mathbf{g}_{\mathbf{q},n+1}^T = \mathbf{g}_{\mathbf{q}}^T(\mathbf{q}_{n+1})$ and the difference equations

$$\tilde{\mathbf{q}}_{n+1} = \mathbf{q}_n + h\mathbf{v}_n + h^2(0.5 - \beta)\mathbf{a}_n + h^2\beta\mathbf{a}_{n+1} \quad (30a)$$

$$\mathbf{q}_{n+1} = \tilde{\mathbf{q}}_{n+1} + \mathbf{U}_{n+1} \quad (30b)$$

$$\tilde{\mathbf{v}}_{n+1} = \mathbf{v}_n + h(1 - \gamma)\mathbf{a}_n + h\gamma\mathbf{a}_{n+1} \quad (30c)$$

$$\mathbf{v}_{n+1} = \tilde{\mathbf{v}}_{n+1} + \mathbf{W}_{n+1} \quad (30d)$$

$$(1 - \alpha_m)\mathbf{a}_{n+1} + \alpha_m \mathbf{a}_n = (1 - \alpha_f)\dot{\tilde{\mathbf{v}}}_{n+1} + \alpha_f \dot{\tilde{\mathbf{v}}}_n \quad (30e)$$

The numerical coefficients $\gamma$, $\beta$, $\alpha_m$, and $\alpha_f$ can be selected from a desired value of the spectral radius at the infinity $\rho_\infty \in [0, 1]$, see [22]

$$\alpha_m = \frac{2\rho_\infty - 1}{\rho_\infty + 1}, \quad \alpha_f = \frac{\rho_\infty}{\rho_\infty + 1}, \quad \gamma = 0.5 + \alpha_f - \alpha_m, \quad \beta = 0.25(\gamma + 0.5)^2 \quad (31)$$

The following remarks can be formulated.

1. The variable $\mathbf{a}_{n+1}$ is a pseudo-acceleration variable which is initialized as $\mathbf{a}_0 = \dot{\tilde{\mathbf{v}}}_0$. This variable can be interpreted as an approximation of the true acceleration $\dot{\tilde{\mathbf{v}}}$ at a shifted time, i.e., $\mathbf{a}_{n+1} \simeq \dot{\tilde{\mathbf{v}}}(t_{n+1} + (\alpha_m - \alpha_f)h)$. Ref. [21] includes a more detailed discussion about the generalized-α method for the simulation of constrained mechanical systems.
2. As shown in [23], the generalized-$\alpha$ method is mathematically equivalent to a two-step method for the velocity variable $\tilde{\mathbf{v}}$ and to a three-step method for the position variable $\tilde{\mathbf{q}}$. Indeed, the value of the acceleration-like variable $\mathbf{a}_n$ depends on the history of $\tilde{\mathbf{v}}$. Therefore, the time integration of $\tilde{\mathbf{v}}$ over the time step $(t_n, t_{n+1}]$ implicitly exploits some values of $\tilde{\mathbf{v}}(t)$ for $t < t_n$. When a jump in the velocity $\mathbf{v}$ occurs, the forces $\mathbf{f}(\mathbf{q}, \mathbf{v}, t)$ in Eq. (8b) are also discontinuous and so is the variable $\tilde{\mathbf{v}}(t)$. Moreover, the variables $\tilde{\mathbf{q}}$ and $\tilde{\mathbf{v}}$ are reinitialized at the beginning of each time step with the values $\mathbf{q}_n$ and $\mathbf{v}_n$, i.e., $\tilde{\mathbf{q}}$ and $\tilde{\mathbf{v}}$ are discontinuous at $t_n$ if an impact has occurred over the time step $(t_{n-1}, t_n]$. In this case, the algorithm, which relies on past values of $\tilde{\mathbf{v}}(t)$, loses its accuracy properties and leads to $O(h)$ errors in the velocity variable $\tilde{\mathbf{v}}(t)$. The proposed algorithm is nevertheless acceptable in the sense that $O(h)$ errors are tolerated anyway for time steps with impacts and its numerical damping can efficiently stabilize the numerical but transient oscillations induced by these disturbances.

3. For numerical reasons, a prediction value $\mathbf{q}_{n+1}^*$ is used in Eq. (29h) instead of the actual value $\mathbf{q}_{n+1}$, similar to what is usually done in the Moreau-Jean scheme. In the following, the prediction $\mathbf{q}_{n+1}^* = \tilde{\mathbf{q}}_{n+1}$ shall be used. The reasons for this choice are discussed in Section 7.2.

4. Equations (29f) and (29h) for the velocity jump $\mathbf{W}_{n+1}$ and the multiplier $\mathbf{\Lambda}_{n+1}$ are similar to the conditions imposed in the classical Moreau-Jean method. The multiplier $\mathbf{\Lambda}_{n+1}$ is thus interpreted as the integral of the contact reaction forces over the time step $(t_n, t_{n+1}]$. Equation (29h) expresses the complementarity between the unilateral constraint at velocity level and the Lagrange multiplier $\mathbf{\Lambda}_{n+1}^j$, which controls the amplitude of the velocity jump.

# 7. Properties and consistency of the scheme

## 7.1. Estimates of the Lagrange multipliers

We can derive an estimate of the multipliers $\nu_{n+1}$ and $\mathbf{\Lambda}_{n+1}$. Let us define $\mathcal{A} \subset \mathcal{C}$ the set of constraints active at position level, which means that $\mathbf{g}^A(\mathbf{q}_{n+1}) = \mathbf{0}$, and $\mathcal{B} \subset \mathcal{C}$ the set of constraints active at velocity level, which means that $\mathbf{g}_{\mathbf{q},n+1}^\mathcal{B} \mathbf{v}_{n+1} + \mathbf{E}^\mathcal{B} \mathbf{g}_{\mathbf{q},n}^\mathcal{B} \mathbf{v}_n = \mathbf{0}$, where $\mathbf{E}^\mathcal{B} = \text{diag}\{e^j\}$ is the diagonal matrix containing the coefficients of restitution of the constraints in $\mathcal{B}$. In this expression, some coefficients of restitution have been artificially introduced for the bilateral constraints and are defined as $e^j = 0, \forall j \in \bar{\mathcal{U}}$. By construction, the bilateral constraints are in both $\mathcal{A}$ and $\mathcal{B}$.

**Property 1.** If the mass matrix is nonsingular, we have

$$ \nu_{n+1}^A = -(\mathbf{g}_{\mathbf{q},n+1}^A \mathbf{M}_{n+1}^{-1} \mathbf{g}_{\mathbf{q},n+1}^{A,T})^{-1} \mathbf{g}^A(\tilde{\mathbf{q}}_{n+1}) + \mathcal{O}(h^2) \quad (32a) $$

$$ \mathbf{\Lambda}_{n+1}^\mathcal{B} = -(\mathbf{g}_{\mathbf{q},n+1}^\mathcal{B} \mathbf{M}_{n+1}^{-1} \mathbf{g}_{\mathbf{q},n+1}^{\mathcal{B},T})^{-1} (\mathbf{g}_{\mathbf{q},n+1}^\mathcal{B} \tilde{\mathbf{v}}_{n+1} + \mathbf{E}^\mathcal{B} \mathbf{g}_{\mathbf{q},n}^\mathcal{B} \mathbf{v}_n) \quad (32b) $$
Proof. We have $g_{q,n+1}^T \Lambda_{n+1} = g_{q,n+1}^{\beta,T} \Lambda_{n+1}^{\beta}$ and $g_{q,n+1}^T \nu_{n+1} = g_{q,n+1}^{A,T} \nu_{n+1}^A$. The proof follows (i) from the approximated form of the active constraints at position level $g^A(q_{n+1}) = g^A(\tilde{q}_{n+1}) + g_{q,n+1}^A U_{n+1} + \mathcal{O}(h^2)$, which results from $U_{n+1} = \mathcal{O}(h)$, with Eq. (29c) and (ii) from the exact form of the active constraints at velocity level $g_{q,n+1}^B v_{n+1} + E^B g_{q,n}^B v_n = g_{q,n+1}^B (\tilde{v}_{n+1} + W_{n+1}) + E^B g_{q,n}^B v_n$ with Eq. (29f). $\square$

Therefore, $\nu_{n+1}$ can be seen as the set of multipliers that appears to compensate the violation of the constraints at position level by the free displacement; this quantity tends to converge to zero when the time step decreases and it is independent of the coefficient of restitution. Similarly, $\Lambda_{n+1}$ can be seen as the set of multipliers that appears to compensate the violation of the constraints at velocity level, i.e., the Newton impact law, by the free velocity; this quantity depends on the value of the coefficient of restitution.

Based on Eq. (20a), the impulse of the constraint reaction forces can be estimated as

$$ \int_{(t_n, t]} \mathrm{d}\boldsymbol{i} \approx \Lambda_{n+1} + h\tilde{\boldsymbol{\lambda}}_{n+1} \quad (33) $$

If impacts occur, $\Lambda_{n+1} = \mathcal{O}(1)$ and the reaction forces are expected to be unbounded. Theorem 1 indicates that $\Lambda_{n+1}$ is affected by an $\mathcal{O}(h)$ local error.

When all constraints remain either closed or open, i.e., no impact occurs, $\Lambda_{n+1} = \mathcal{O}(h)$ and a more refined analysis reveals that this variable is only affected by an $\mathcal{O}(h^2)$ local error. The variable $(1/h)\Lambda_{n+1} + \tilde{\boldsymbol{\lambda}}_{n+1}$ is thus a meaningful approximation the constraint reaction forces. Since $\tilde{\boldsymbol{\lambda}}^U = \mathbf{0}$, the unilateral contact forces are represented by $(1/h)\Lambda_{n+1}^U$. Also, the generalized reaction forces of the bilateral constraints can be estimated as $(1/h)\Lambda_{n+1}^{\overline{U}} + \tilde{\boldsymbol{\lambda}}_{n+1}^{\overline{U}}$. Illustrations of this estimate will be given in Section 10.2.

## 7.2. Choice of the prediction value for the activation of the velocity constraint

Different choices for the prediction $q_{n+1}^*$ in Eq. (29h) could be considered such as $q_{n+1}^* = \tilde{q}_{n+1}$, $q_{n+1}^* = q_n$, $q_{n+1}^* = q_n + h\nu_n$ or $q_{n+1}^* = q_{n+1}$, which would change the behaviour of the algorithm. It may seem attractive to use $q_{n+1}^* = q_{n+1}$ in order to prevent the activation of a constraint at velocity level when it is not active at position level. However, a tolerance on the criterion $g^j(q_{n+1}^*) \le tol$ is then needed where tol is defined as an upper bound on the numerical errors for the satisfaction of the constraint at position level. This tolerance would be particularly important for systems with nonlinear unilateral constraints, since any iterative Newton-type algorithm would only converge to the position level constraint up to certain tolerance. However, if the tolerance is well-chosen, the behavior of the algorithm would be similar as if $q_{n+1}^* = \tilde{q}_{n+1}$, because Eq. (32a) with the condition $\nu_{n+1}^A \ge \mathbf{0}$ imply that $g^A(\tilde{q}_{n+1}) \le \mathbf{0}$ for a sufficiently small time step. For that reason, the choice $q_{n+1}^* = \tilde{q}_{n+1}$ shall be used in the following since it does not require the selection of a tolerance.
### 7.3. Fundamental properties in a one degree-of-freedom test problem

In this section, the scheme is studied for the simple example of a bouncing ball, which can be seen as a single degree-of-freedom model equation for the contact problem.

The equations of motion of the bouncing ball take the form

$$ \dot{q} = v \qquad (34a) $$

$$ mdv + di = -ma_g dt \qquad (34b) $$

$$ 0 \le q \perp di \ge 0 \qquad (34c) $$

where $m$ is the mass of the ball and $a_g$ is the gravity acceleration.

According to the proposed time stepping scheme, the numerical solution is such that $\dot{v}_{n+1} = a_{n+1} = -a_g$ and

$$ mU_{n+1} - \nu_{n+1} = 0 \qquad (35a) $$

$$ 0 \le q_{n+1} \perp \nu_{n+1} \ge 0 \qquad (35b) $$

$$ mW_{n+1} - \Lambda_{n+1} = 0 \qquad (35c) $$

$$ \text{if } \tilde{q}_{n+1} \le 0 \text{ then } 0 \le v_{n+1} + e v_n \perp \Lambda_{n+1} \ge 0 \qquad (35d) $$

$$ \tilde{q}_{n+1} = q_n + hv_n - h^2 a_g / 2 \qquad (35e) $$

$$ q_{n+1} = \tilde{q}_{n+1} + U_{n+1} \qquad (35f) $$

$$ \tilde{v}_{n+1} = v_n - ha_g \qquad (35g) $$

$$ v_{n+1} = \tilde{v}_{n+1} + W_{n+1} \qquad (35h) $$

The next property is related with the equilibrium points of the numerical solution. Let us recall that a state $q_n$, $v_n$ of the discrete system is an equilibrium point (or a fixed point) if we have $q_n = q_e$ and $v_n = 0 \forall n$.

**Property 2.** For the system represented by Eq. (35), the state $q_n = 0, v_n = 0$ is the unique equilibrium point.

*Proof.* Let us assume we are at an equilibrium point with $q_n = q_e$, $v_n = 0$. Then, we have

$$ \tilde{q}_{n+1} = q_e - h^2 a_g / 2 $$

$$ \tilde{v}_{n+1} = -ha_g $$

If $q_e > h^2 a_g / 2$, the complementarity condition at velocity level is not activated so that $W_{n+1} = 0, v_{n+1} = -ha_g$ and the state is not an equilibrium point.

If $q_e \in [0, h^2 a_g / 2)$, after some manipulations, the complementarity condition at velocity level becomes

$$ 0 \le -ha_g + W_{n+1} \perp mW_{n+1} \ge 0 $$
which yields $W_{n+1} = ha_g$ so that $\Lambda_{n+1} = mha_g$ and $v_{n+1} = 0$, i.e., the constraint is active at velocity level. The complementarity condition at position level can be reformulated as

$$0 \le q_e - h^2 a_g / 2 + U_{n+1} \perp mU_{n+1} \ge 0$$

$\forall q_e \in [0, h^2 a_g / 2)$, the solution is then obtained as

$$q_{n+1} = 0, \quad U_{n+1} \ge 0 \qquad (36)$$

Any $q_e \neq 0$ yields $q_{n+1} = 0 \neq q_e$ and is thus not an equilibrium point. The value $q_e = 0$ yields $q_{n+1} = 0$ and is thus the unique equilibrium point. □

**Property 3.** The numerical solution of the bouncing ball problem in Eq. (34) is such that the total energy decays monotonically

$$E_{n+1} \le E_n$$

where $E_n = mv_n^2/2 + ma_gq_n$.

*Proof.* The proof relies on the analysis of the different activation scenarios of the unilateral constraint.

The first situation is the free flight phase. Since the smooth motion is solved exactly, the energy is exactly preserved by the algorithm.

The second situation is a simultaneous activation of the constraint at position and velocity levels, which implies $v_{n+1} = -e v_n$ and $q_{n+1} = 0$. The variation of energy is then equal to

$$E_{n+1} - E_n = -q_n ma_g - \frac{1}{2} m(1 - e^2)v_n^2$$

Using $q_n \ge 0$ and $e \in [0, 1]$, we obtain $E_{n+1} - E_n \le 0$.

The third situation is such that the constraint is not active at velocity level because $\tilde{v}_{n+1} + e v_n > 0$ but is active at position level because $\tilde{q}_{n+1} < 0$. We would have simultaneously

$$\tilde{q}_{n+1} = q_n + h v_n - \frac{1}{2} h^2 a_g < 0 \\ \tilde{v}_{n+1} + e v_n = v_n - h a_g + e v_n > 0$$

which would imply that

$$\frac{h^2 a_g}{1 + e} < h v_n < \frac{h^2 a_g}{2} - q_n$$

Since $e \in [0, 1]$ and $q_n \ge 0$, these conditions can never be satisfied so that this situation cannot be encountered.

The fourth situation is such that the constraint is active at velocity level and not at position level. This situation is forbidden in the proposed algorithm. □
## 8. Newton semi-smooth method

This section addresses the solution of the discrete nonlinear system (29,30) for a given time step $n+1$, using a Newton procedure with an active set strategy. The algorithm is described in the special case $q_{n+1}^* = \tilde{q}_{n+1}$.

The complementarity conditions in Eq. (29) can be reformulated to get a nonlinear semi-smooth problem as follows

$$ \mathbf{M}(\mathbf{q}_{n+1})\dot{\tilde{\mathbf{v}}}_{n+1} - \mathbf{f}(\mathbf{q}_{n+1}, \mathbf{v}_{n+1}, t_{n+1}) - \mathbf{g}_{\mathbf{q},n+1}^{\overline{\mathcal{U}},T} \tilde{\mathbf{\lambda}}_{n+1}^{\overline{\mathcal{U}}} = \mathbf{0} \quad (37a) $$

$$ \mathbf{g}_{\mathbf{q},n+1}^{\overline{\mathcal{U}}} \tilde{\mathbf{v}}_{n+1} = \mathbf{0} \quad (37b) $$

$$ \mathbf{M}(\mathbf{q}_{n+1})\mathbf{U}_{n+1} - \mathbf{g}_{\mathbf{q},n+1}^{T} \boldsymbol{\nu}_{n+1} = \mathbf{0} \quad (37c) $$

$$ \mathbf{g}^{\overline{\mathcal{U}}}(\mathbf{q}_{n+1}) = \mathbf{0} \quad (37d) $$

$$ \boldsymbol{\nu}_{n+1}^{\mathcal{U}} - \max(\mathbf{0}, \boldsymbol{\nu}_{n+1}^{\mathcal{U}} - r \mathbf{g}^{\mathcal{U}}(\mathbf{q}_{n+1})) = \mathbf{0} \quad (37e) $$

$$ \mathbf{M}(\mathbf{q}_{n+1})\mathbf{W}_{n+1} - \mathbf{g}_{\mathbf{q},n+1}^{T} \mathbf{\Lambda}_{n+1} = \mathbf{0} \quad (37f) $$

$$ \mathbf{g}_{\mathbf{q},n+1}^{\overline{\mathcal{U}}} \mathbf{v}_{n+1} = \mathbf{0} \quad (37g) $$

$$ \text{if } g^j(\tilde{\mathbf{q}}_{n+1}) \le 0 \text{ then } \mathbf{\Lambda}^j - \max(0, \mathbf{\Lambda}^j - r(g_{\mathbf{q},n+1}^j \mathbf{v}_{n+1} + e^j g_{\mathbf{q},n}^j \mathbf{v}_n)) = 0, \quad (37h) \\ \forall j \in \mathcal{U} $$

with the numerical coefficient $r > 0$. This nonlinear semi-smooth problem leads to a monolithic scheme to iteratively compute the solution and contact conditions. The solution of this semi-smooth problem is used to advance the computation at each time step.

The active sets $\mathcal{A}_{n+1}$ and $\mathcal{B}_{n+1}$ and their complementarity sets $\overline{\mathcal{A}}_{n+1}$ and $\overline{\mathcal{B}}_{n+1}$ are defined as the sets of constraints

$$ \mathcal{A}_{n+1} = \overline{\mathcal{U}} \cup \{j \in \mathcal{U} : \nu_{n+1}^j - r g^j(\mathbf{q}_{n+1}) \ge 0\} \quad (38a) $$

$$ \overline{\mathcal{A}}_{n+1} = \mathcal{C} \setminus \mathcal{A}_{n+1} \quad (38b) $$

$$ \mathcal{B}_{n+1} = \overline{\mathcal{U}} \cup \{j \in \mathcal{U} : g^j(\tilde{\mathbf{q}}_{n+1}) \le 0 \text{ and} \\ \qquad \Lambda_{n+1}^j - r(g_{\mathbf{q},n+1}^j \mathbf{v}_{n+1} + e^j g_{\mathbf{q},n}^j \mathbf{v}_n) \ge 0\} \quad (38c) $$

$$ \overline{\mathcal{B}}_{n+1} = \mathcal{C} \setminus \mathcal{B}_{n+1} \quad (38d) $$

One observes that these two activation criteria are completely implicit and have a similar structure as the activation criterion proposed by Alart and Curnier [24]
for quasi-static problems. The scaled vectors of residuals are defined as follows:

$$ \mathbf{r}^s = \left[ \begin{array}{c} \mathbf{M}(\mathbf{q}_{n+1})h\dot{\tilde{\mathbf{v}}}_{n+1} - \mathbf{g}_{\mathbf{q},n+1}^{\bar{\mathcal{U}},T} h\tilde{\lambda}_{n+1}^{\bar{\mathcal{U}}} - h\mathbf{f}(\mathbf{q}_{n+1}, \mathbf{v}_{n+1}, t_{n+1}) \\ \mathbf{g}_{\mathbf{q},n+1}^{\bar{\mathcal{U}}} \tilde{\mathbf{v}}_{n+1} \end{array} \right] \quad (39a) $$

$$ \mathbf{r}^p = \frac{1}{h} \left[ \begin{array}{c} \mathbf{M}(\mathbf{q}_{n+1})\mathbf{U}_{n+1} - \mathbf{g}_{\mathbf{q},n+1}^{A,T} \boldsymbol{\nu}_{n+1}^A \\ \mathbf{g}^A(\mathbf{q}_{n+1}) \\ \boldsymbol{\nu}_{n+1}^{\bar{\mathcal{A}}} \end{array} \right] \quad (39b) $$

$$ \mathbf{r}^v = \left[ \begin{array}{c} \mathbf{M}(\mathbf{q}_{n+1})\mathbf{W}_{n+1} - \mathbf{g}_{\mathbf{q},n+1}^{\beta,T} \mathbf{\Lambda}_{n+1}^{\beta} \\ \mathbf{g}_{\mathbf{q},n+1}^{\beta} \mathbf{v}_{n+1} + \mathbf{E}^{\beta} \mathbf{g}_{\mathbf{q},n}^{\beta} \mathbf{v}_{n} \\ \mathbf{\Lambda}^{\bar{\mathcal{B}}} \end{array} \right] \quad (39c) $$

$$ \mathbf{r} = \left[ \begin{array}{c} \mathbf{r}^s \\ \mathbf{r}^p \\ \mathbf{r}^v \end{array} \right] \quad (39d) $$

In the above expressions, $\mathcal{A}_{n+1}$ and $\mathcal{B}_{n+1}$ are denoted by $\mathcal{A}$ and $\mathcal{B}$ for the sake of readability, whereas the matrix of restitution coefficients $\mathbf{E}^{\beta}$ was introduced in Section 7.1. It can be observed that the numerical coefficient $r$ only appears in the activation criteria but not in the residuals (39), which is an important difference compared to the Augmented Lagrangian formulation described in [24]. Let us remark that the residuals could be reorganized in a form similar to an Augmented Lagrangian formulation by the addition of penalty terms, i.e., by replacing the terms

$$ \mathbf{g}_{\mathbf{q},n+1}^{\bar{\mathcal{U}},T} h \tilde{\lambda}_{n+1}^{\bar{\mathcal{U}}}, \mathbf{g}_{\mathbf{q},n+1}^{A,T} \boldsymbol{\nu}_{n+1}^{A}, \mathbf{g}_{\mathbf{q},n+1}^{\beta,T} \mathbf{\Lambda}_{n+1}^{\beta} $$

by

$$ \mathbf{g}_{\mathbf{q},n+1}^{\bar{\mathcal{U}},T} h(\tilde{\lambda}_{n+1}^{\bar{\mathcal{U}}} - r \mathbf{g}_{\mathbf{q},n+1}^{\bar{\mathcal{U}}} \tilde{\mathbf{v}}_{n+1}), \mathbf{g}_{\mathbf{q},n+1}^{A,T} (\boldsymbol{\nu}_{n+1}^{A} - r \mathbf{g}^{A}(\mathbf{q}_{n+1})), \\ \mathbf{g}_{\mathbf{q},n+1}^{\beta,T} (\mathbf{\Lambda}_{n+1}^{\beta} - r(\mathbf{g}_{\mathbf{q},n+1}^{\beta} \mathbf{v}_{n+1} + \mathbf{E}^{\beta} \mathbf{g}_{\mathbf{q},n}^{\beta} \mathbf{v}_{n})), $$

respectively. These penalty terms would not modify the numerical solution but they could tend to bring some convexity to the problem, thereby improving the convergence of the Newton semi-smooth iterations. However, this reformulation was not necessary in the test cases described below.

The prediction is based on the initial guesses $\tilde{\mathbf{v}}_{n+1}^0 = \mathbf{0}$, $\mathbf{W}_{n+1}^0 = \mathbf{0}$, $\mathbf{U}_{n+1}^0 = \mathbf{0}$,
and $\dot{\mathbf{q}}_n = \mathbf{v}_n$, which implies:

$$ \mathbf{q}_{n+1}^0 = \mathbf{q}_n + h\mathbf{v}_n + h^2(0.5 - \beta)\mathbf{a}_n + h^2\beta(\alpha_f \dot{\mathbf{v}}_n - \alpha_m \mathbf{a}_n)/(1 - \alpha_m) \quad (40a) $$

$$ \tilde{\mathbf{v}}_{n+1}^0 = \mathbf{v}_n + h(1 - \gamma)\mathbf{a}_n + h\gamma(\alpha_f \dot{\mathbf{v}}_n - \alpha_m \mathbf{a}_n)/(1 - \alpha_m) \quad (40b) $$

$$ \mathbf{a}_{n+1}^0 = (\alpha_f \dot{\mathbf{v}}_n - \alpha_m \mathbf{a}_n)/(1 - \alpha_m) \quad (40c) $$

$$ \mathbf{v}_{n+1}^0 = \tilde{\mathbf{v}}_{n+1}^0 \quad (40d) $$

$$ \mathbf{U}_{n+1}^0 = \mathbf{0} \quad (40e) $$

$$ \mathbf{W}_{n+1}^0 = \mathbf{0} \quad (40f) $$

$$ \tilde{\mathbf{\lambda}}_{n+1}^0 = \mathbf{0} \quad (40g) $$

$$ \nu_{n+1}^0 = \mathbf{0} \quad (40h) $$

$$ \Lambda_{n+1}^0 = \mathbf{0} \quad (40i) $$

Then, the correction terms should satisfy the difference formulae, so that the corrections $\Delta\tilde{\mathbf{v}}_{n+1}$ and $\Delta\mathbf{q}_n$ can be eliminated in terms of the corrections $\Delta\tilde{\mathbf{v}}_{n+1}$, $\Delta\mathbf{W}_{n+1}$ and $\Delta\mathbf{U}_{n+1}$.

$$ \Delta\mathbf{v}_{n+1} = \Delta\tilde{\mathbf{v}}_{n+1} + \Delta\mathbf{W}_{n+1} \quad (41a) $$

$$ \Delta\tilde{\mathbf{v}}_{n+1} = (1 - \alpha_m)/((1 - \alpha_f)\gamma h) \Delta\tilde{\mathbf{v}}_{n+1} \quad (41b) $$

$$ \Delta\mathbf{q}_{n+1} = h\beta/\gamma \Delta\tilde{\mathbf{v}}_{n+1} + \Delta\mathbf{U}_{n+1} \quad (41c) $$

Introducing the vector of independent corrections

$$ \Delta \mathbf{x} = \begin{bmatrix} \Delta \tilde{\mathbf{v}} \\ -h \Delta \tilde{\mathbf{\lambda}}^{\bar{\mathcal{U}}} \\ \Delta \mathbf{U}/h \\ -\Delta \mathbf{\nu}^{\mathcal{A}}/h \\ -\Delta \mathbf{\nu}^{\bar{\mathcal{A}}}/h \\ \Delta \mathbf{W} \\ -\Delta \mathbf{\Lambda}^{\mathcal{B}} \\ -\Delta \mathbf{\Lambda}^{\bar{\mathcal{B}}} \end{bmatrix} \quad (42) $$

the correction equation is obtained as

$$ \mathbf{S}_t \Delta \mathbf{x} = -\mathbf{r} \quad (43) $$

In this expression, the iteration matrix is defined as

$$ \mathbf{S}_t = \begin{bmatrix} \mathbf{S}_t^* & \mathbf{g}_{\mathbf{q}}^{\bar{\mathcal{U}},T} & h^2 \mathbf{K}_t & 0 & 0 & h \mathbf{C}_t & 0 & 0 \\ \mathbf{g}_{\mathbf{q}}^{\bar{\mathcal{U}}} + \frac{h\beta}{\gamma} \mathbf{G}^{*s} & \mathbf{0} & h \mathbf{G}^{*s} & 0 & 0 & \mathbf{0} & 0 & 0 \\ \frac{\beta}{\gamma} \mathbf{G}^p & \mathbf{0} & \mathbf{M} + \mathbf{G}^p & \mathbf{g}_{\mathbf{q}}^{\mathcal{A},T} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \frac{\beta}{\gamma} \mathbf{g}_{\mathbf{q}}^{\mathcal{A}} & \mathbf{0} & \mathbf{g}_{\mathbf{q}}^{\mathcal{A}} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{I}^{\bar{\mathcal{A}}} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \frac{h\beta}{\gamma} \mathbf{G}^v & \mathbf{0} & h \mathbf{G}^v & \mathbf{0} & \mathbf{0} & \mathbf{M} & \mathbf{g}_{\mathbf{q}}^{\mathcal{B},T} & \mathbf{0} \\ \mathbf{g}_{\mathbf{q}}^{\mathcal{B}} + \frac{h\beta}{\gamma} \mathbf{G}^{*v} & \mathbf{0} & h \mathbf{G}^{*v} & \mathbf{0} & \mathbf{0} & \mathbf{g}_{\mathbf{q}}^{\mathcal{B}} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{I}^{\bar{\mathcal{B}}} \end{bmatrix} \quad (44) $$
where $\overline{I^A}$ and $\overline{I^B}$ are identity matrices and

$$ 
\begin{aligned}
\mathbf{S}_t^* &= \frac{1 - \alpha_m}{(1 - \alpha_f)\gamma} \mathbf{M} + h\mathbf{C}_t + \frac{\beta h^2}{\gamma} \mathbf{K}_t, \quad \mathbf{K}_t = \frac{\partial(\mathbf{M}(\mathbf{q})\dot{\tilde{\mathbf{v}}} - \mathbf{f}(\mathbf{q}, \mathbf{v}, t) - \mathbf{g}_{\mathbf{q}}^{\overline{\mathcal{U}},T}\tilde{\mathbf{\lambda}}^{\overline{\mathcal{U}}})}{\partial \mathbf{q}}, \\
\mathbf{C}_t &= -\frac{\partial \mathbf{f}(\mathbf{q}, \mathbf{v}, t)}{\partial \mathbf{v}}, \quad \mathbf{G}^p = \frac{\partial(\mathbf{M}(\mathbf{q})\mathbf{U} - \mathbf{g}_{\mathbf{q}}^{A,T}\boldsymbol{\nu}^A)}{\partial \mathbf{q}}, \quad \mathbf{G}^v = \frac{\partial(\mathbf{M}(\mathbf{q})\mathbf{W} - \mathbf{g}_{\mathbf{q}}^{\mathcal{B},T}\boldsymbol{\Lambda}^\mathcal{B})}{\partial \mathbf{q}}, \\
\mathbf{G}^{*s} &= \frac{\partial(\mathbf{g}_{\mathbf{q}}^{\overline{\mathcal{U}}\tilde{\mathbf{v}}})}{\partial \mathbf{q}}, \quad \mathbf{G}^{*v} = \frac{\partial(\mathbf{g}_{\mathbf{q}}^{\mathcal{B}\mathbf{v}})}{\partial \mathbf{q}}
\end{aligned}
\quad (45) $$

with the tangent stiffness matrix $\mathbf{K}_t$ and the tangent damping matrix $\mathbf{C}_t$. The matrices $\mathbf{G}^p$, $\mathbf{G}^v$, $\mathbf{G}^{*s}$ and $\mathbf{G}^{*v}$ are related with the second derivatives of the constraints and they vanish for linear constraints. The scaling of the multipliers in Eq. (42) by negative numbers has been selected in order to obtain symmetric contributions in the iteration matrix $\mathbf{S}_t$ and to ensure that this matrix does not become ill-conditioned for small time steps, see, e.g., [25] for a more detailed discussion on DAE scaling strategies. Indeed, since $\mathbf{U} = \mathcal{O}(h)$ and $\boldsymbol{\nu} = \mathcal{O}(h)$, we have $\mathbf{G}^p = \mathcal{O}(h)$ and one can observe that

$$ \lim_{h \to 0} \mathbf{S}_t = \begin{bmatrix} \frac{1 - \alpha_m}{(1 - \alpha_f)\gamma} \mathbf{M} & \mathbf{g}_{\mathbf{q}}^{\overline{\mathcal{U}},T} & 0 & 0 & 0 & 0 & 0 & 0 \\ \mathbf{g}_{\mathbf{q}}^{\overline{\mathcal{U}}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \frac{\beta}{\gamma} \mathbf{G}^p & 0 & \mathbf{M} & \mathbf{g}_{\mathbf{q}}^{A,T} & 0 & 0 & 0 & 0 \\ \frac{\beta}{\gamma} \mathbf{g}_{\mathbf{q}}^{A} & 0 & \mathbf{g}_{\mathbf{q}}^{A} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \overline{\mathbf{I}^A} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \mathbf{M} & \mathbf{g}_{\mathbf{q}}^{\mathcal{B},T} & 0 \\ \mathbf{g}_{\mathbf{q}}^{\mathcal{B}} & 0 & 0 & 0 & 0 & \mathbf{g}_{\mathbf{q}}^{\mathcal{B}} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \overline{\mathbf{I}^B} \end{bmatrix} \quad (46) $$

is a non-singular and well-conditioned matrix if we assume that the matrix in Eq (9) is also well-conditioned.

If the tangent damping and stiffness matrices are equal to 0 (which means that the smooth force vector $\mathbf{f}$ does not depend on $\mathbf{q}$ or $\mathbf{v}$) and if the bilateral constraints are linear (which implies $\mathbf{G}^{*s} = 0$), a lower-triangular block matrix is obtained and the system can be solved based on a sequential solution of the three subproblems. In the general case, this sequential approach also converges in many practical situations, even if the coupling terms are neglected at each iteration. This sequential approach can simplify the implementation in a general purpose code and improve the efficiency of the simulation. The complete procedure is
described in Algorithm 1 and involves the following matrices

$$ S_t^s = \begin{bmatrix} S_t^* & g_{\bar{\mathcal{U}},T} \\ g_{\bar{\mathcal{U}}} + \frac{h\beta}{\gamma} G^{*s} & 0 \end{bmatrix} \quad (47) $$

$$ S_t^p = \begin{bmatrix} M + G^p & g_{\bar{\mathcal{U}}}^{A,T} & 0 \\ g_{\bar{\mathcal{U}}}^A & 0 & 0 \\ 0 & 0 & I^{\bar{\mathcal{A}}} \end{bmatrix} \quad (48) $$

$$ S_t^v = \begin{bmatrix} M & g_{\bar{\mathcal{U}}}^{\beta,T} & 0 \\ g_{\bar{\mathcal{U}}}^{\beta} & 0 & 0 \\ 0 & 0 & I^{\bar{\mathcal{B}}} \end{bmatrix} \quad (49) $$

It can be summarized as follows. At each iteration of the Newton procedure, Algorithm 1 involves 3 steps:

1. solve the linearized equation of the smooth motion (i.e., all contributions of unilateral constraints are disregarded) using the classical procedure obtained from the generalized-α scheme for smooth systems with the bilateral constraints expressed at velocity level;

2. project on the position constraints in $\mathcal{A}$ using the matrix $S_t^p$;

3. project on the velocity constraints in $\mathcal{B}$ using the matrix $S_t^v$.

If convergence difficulties are encountered, one may consider either to reduce the time step, which tends to reduce the influence of the coupling terms $h^2\mathbf{K}_t$, $h\mathbf{G}^{*s}$, and $h\mathbf{C}_t$, or to solve the coupled problem based on the full matrix $\mathbf{S}_t$ in Eq. (44).

## 9. Alternative schemes

In the numerical examples below, the proposed nonsmooth generalized-α GGL method is compared with other two algorithms from the literature, which are briefly described in this section.

### 9.1. Moreau-Jean $\theta$-method

Assuming that the mass matrix is non-singular, the Moreau-Jean $\theta$-method can be formulated based on the following discrete set of equations

$$ M(\mathbf{q}_{n+1})\dot{\mathbf{v}}_{n+1} - \mathbf{f}(\mathbf{q}_{n+1}, \mathbf{v}_{n+1}, t_{n+1}) = \mathbf{0} \quad (50a) $$

$$ M(\mathbf{q}_{n+1})\mathbf{W}_{n+1} - \mathbf{g}_{\mathbf{q},n+1}^T \mathbf{\Lambda}_{n+1} = \mathbf{0} \quad (50b) $$

$$ \mathbf{g}_{\mathbf{q},n+1}^{\bar{\mathcal{U}}} \mathbf{v}_{n+1} = \mathbf{0} \quad (50c) $$

$$ \text{if } g^j(\mathbf{q}_{n+1}^*) \le 0 \text{ then } \Lambda^j - \max(0, \Lambda^j - r(g_{\mathbf{q},n+1}^j \mathbf{v}_{n+1} + e^j g_{\mathbf{q},n}^j \mathbf{v}_n)) = 0, \quad (50d) \\ \forall j \in \mathcal{U} $$

with the integration formulae

$$ \tilde{\mathbf{q}}_{n+1} = \mathbf{q}_n + h(1-\theta)\mathbf{v}_n + h\theta\tilde{\mathbf{v}}_{n+1} \quad (50e) $$

$$ \tilde{\mathbf{v}}_{n+1} = \mathbf{v}_n + h(1-\theta)\dot{\tilde{\mathbf{v}}}_n + h\theta\dot{\tilde{\mathbf{v}}}_{n+1} \quad (50f) $$

$$ \mathbf{v}_{n+1} = \tilde{\mathbf{v}}_{n+1} + \mathbf{W}_{n+1} \quad (50g) $$

$$ \mathbf{q}_{n+1} = \tilde{\mathbf{q}}_{n+1} + h\theta\mathbf{W}_{n+1} \quad (50h) $$
**Algorithm 1** Nonsmooth GGL generalized-α time integration scheme

Inputs: initial values $\mathbf{q}_0$ and $\mathbf{v}_0$

Compute the consistent value of $\dot{\tilde{\mathbf{v}}}_0$

$$ \mathbf{a}_0 := \dot{\tilde{\mathbf{v}}}_0 $$

for $n = 0$ to $n_{\text{final}} - 1$ do

$$ \dot{\tilde{\mathbf{v}}}_{n+1} := \mathbf{0}, \tilde{\mathbf{\lambda}}_{n+1}^{\mathbf{U}} := \mathbf{0}, \mathbf{\nu}_{n+1} := \mathbf{0}, \mathbf{\Lambda}_{n+1} := \mathbf{0} $$

$$ \mathbf{a}_{n+1} := \frac{1}{1 - \alpha_m}(\alpha_f \dot{\tilde{\mathbf{v}}}_n - \alpha_m \mathbf{a}_n) $$

$$ \dot{\tilde{\mathbf{v}}}_{n+1} := \mathbf{v}_n + h(1 - \gamma)\mathbf{a}_n + h\gamma\mathbf{a}_{n+1} $$

$$ \mathbf{v}_{n+1} := \tilde{\mathbf{v}}_{n+1} $$

$$ \mathbf{q}_{n+1} := \mathbf{q}_n + h\mathbf{v}_n + h^2(1/2 - \beta)\mathbf{a}_n + h^2\beta\mathbf{a}_{n+1} $$

for $i = 1$ to $i_{\text{max}}$ do

Compute $\mathbf{r}^s$, $\mathbf{r}^p$ and $\mathbf{r}^v$ using Eqs. (39a-39c)

if $\|\mathbf{r}^s\| < \text{tol}^s$ and $\|\mathbf{r}^p\| < \text{tol}^p$ and $\|\mathbf{r}^v\| < \text{tol}^v$ then
break

end if

Step 1 (smooth motion):

Compute the iteration matrix $\mathbf{S}_t^s$ using Eq. (47)

$$ \begin{bmatrix} \Delta \tilde{\mathbf{v}} \\ -h \Delta \tilde{\mathbf{\lambda}}^{\mathbf{U}} \end{bmatrix} := -(\mathbf{S}_t^s)^{-1} \mathbf{r}^s $$

$$ \tilde{\mathbf{v}}_{n+1} := \tilde{\mathbf{v}}_{n+1} + \Delta \tilde{\mathbf{v}} $$

$$ \dot{\tilde{\mathbf{v}}}_{n+1} := \dot{\tilde{\mathbf{v}}}_{n+1} + \frac{1 - \alpha_m}{(1 - \alpha_f)\gamma h} \Delta \tilde{\mathbf{v}} $$

$$ \mathbf{v}_{n+1} := \tilde{\mathbf{v}}_{n+1} + \mathbf{W}_{n+1} $$

$$ \mathbf{q}_{n+1} := \mathbf{q}_{n+1} + h\beta/\gamma \Delta \tilde{\mathbf{v}} $$

$$ \tilde{\mathbf{\lambda}}_{n+1}^{\mathbf{U}} := \tilde{\mathbf{\lambda}}_{n+1}^{\mathbf{U}} + \Delta \tilde{\mathbf{\lambda}}^{\mathbf{U}} $$

Step 2 (projection on position constraints):

Compute $\mathbf{r}^p$ using Eq. (39b)

Compute $\mathbf{S}_t^p$ using Eq. (48)

$$ \begin{bmatrix} \Delta \mathbf{U} \\ -\Delta \mathbf{\nu}/h \end{bmatrix} := -(\mathbf{S}_t^p)^{-1} \mathbf{r}^p $$

$$ \mathbf{U}_{n+1} := \mathbf{U}_{n+1} + \Delta \mathbf{U} $$

$$ \mathbf{q}_{n+1} := \mathbf{q}_{n+1} + \Delta \mathbf{U} $$

$$ \mathbf{\nu}_{n+1} := \mathbf{\nu}_{n+1} + \Delta \mathbf{\nu} $$

Step 3 (projection on velocity constraints):

Compute $\mathbf{r}^v$ using Eq. (39c)

Compute $\mathbf{S}_t^v$ using Eq. (49)

$$ \begin{bmatrix} \Delta \mathbf{W} \\ -\Delta \mathbf{\Lambda} \end{bmatrix} := -(\mathbf{S}_t^v)^{-1} \mathbf{r}^v $$

$$ \mathbf{W}_{n+1} := \mathbf{W}_{n+1} + \Delta \mathbf{W} $$

$$ \mathbf{v}_{n+1} := \tilde{\mathbf{v}}_{n+1} + \mathbf{W}_{n+1} $$

$$ \mathbf{\Lambda}_{n+1} := \mathbf{\Lambda}_{n+1} + \Delta \mathbf{\Lambda} $$

end for

$$ \mathbf{a}_{n+1} := \mathbf{a}_{n+1} + \frac{1 - \alpha_f}{1 - \alpha_m} \dot{\tilde{\mathbf{v}}}_{n+1} $$

end for
This system of equations can be solved using a Newton method, solving the complementarity condition in Eq. (50d) at every iteration. It is important that the activation of the complementarity condition in Eq. (50d) relies on an evaluation of the constraint $g^j$ at the position $q_{n+1}^*$, where $q_{n+1}^*$ is an explicit prediction of the position which is not updated during the Newton iteration process. The reason is that, if the activation was based on the constraint at position level at time $t_{n+1}$, i.e., on the sign of $g^j(q_{n+1})$, the discrete problem may not always be solvable and the Newton iterations may fail to converge. As a consequence, the Moreau-Jean scheme does not provide any guarantee on the sign of the unilateral constraint at position-level for the final response.

## 9.2. Smooth generalized-α index-3 method

This algorithm relies on a formulation of the constraints at position level only. It assumes that there is no jump in the velocities and is therefore only consistent if no impact occurs.

$$M(q_{n+1})\dot{v}_{n+1} - g_{q,n+1}^T \lambda_{n+1} = f(q_{n+1}, v_{n+1}, t_{n+1}) \quad (51a)$$

$$g_{q,n+1}^{\overline{U}}(q_{n+1}) = 0 \quad (51b)$$

$$\lambda_{n+1}^U - \max(0, \lambda_{n+1}^U - rg^U(q_{n+1})) = 0 \quad (51c)$$

which are combined with the integration formulae

$$q_{n+1} = q_n + h\mathbf{v}_n + h^2(0.5 - \beta)\mathbf{a}_n + h^2\beta\mathbf{a}_{n+1} \quad (52a)$$

$$\tilde{\mathbf{v}}_{n+1} = \mathbf{v}_n + h(1 - \gamma)\mathbf{a}_n + h\gamma\mathbf{a}_{n+1} \quad (52b)$$

$$(1 - \alpha_m)\mathbf{a}_{n+1} + \alpha_m \mathbf{a}_n = (1 - \alpha_f)\dot{\mathbf{v}}_{n+1} + \alpha_f \dot{\mathbf{v}}_n \quad (52c)$$

It is noticeable that this algorithm does not make use of the impact law, i.e., the value of the coefficient of restitution is not needed to obtain a numerical solution.

# 10. Numerical examples

In order to highlight the properties of the proposed time integration method, several examples are studied as shown in Figure 1. The bouncing ball and the elastic bar have a linear behaviour in free motion and can be compared to exact solutions. The impact of a simple rigid pendulum illustrates the behaviour of the scheme when unilateral and bilateral constraints are simultaneously present. The rocking block example is modelled using two nonlinear unilateral constraints and its motion exhibits multiple impacts.

Several algorithms are compared:

* Smooth-α index-3: the generalized-α scheme with constraints at position level only, i.e., no coefficient of restitution is used, the numerical coefficients are defined according to [22] in order to have a spectral radius $\rho_\infty = 0.8$,
<figure><img src="image_2.png" /></figure>

(a)

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

(b)

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

(c)

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

(d)

<caption>Figure 1: Examples: (a) bouncing ball; (b) rigid pendulum; (c) elastic bar; (d) rocking block.</caption>

• Nonsmooth-α, modified LCP: a variant of the proposed scheme when the assumption represented by Eq. (28) is relaxed and the complementarity condition (29e) is replaced by Eq. (25) with the approximation of $\mu^\mu$ in Eq. (27), the numerical coefficients are defined according to [22] in order to have a spectral radius $\rho_\infty = 0.8$,

• Moreau: the Moreau $\theta$-method, where $\theta$ is selected in order to have a spectral radius $\rho_\infty = 0.8$ using the formula $\theta = 1/(\rho_\infty + 1)$,

• Nonsmooth-α GGL: the proposed GGL approach with constraints at velocity and position level, the numerical coefficients are defined according to [22] in order to have a spectral radius $\rho_\infty = 0.8$,

• Reference: a reference solution computed either with Nonsmooth-α GGL scheme and a smaller time step, or analytically.

### 10.1. The bouncing ball

The first example is a standard bouncing rigid ball on a rigid plane, as depicted in Figure 1(a). The ball starts from rest, falls due to the gravity and rebounds against the rigid plane with a restitution coefficient $e = 0.8$. Thus, the model includes a unilateral constraint on the vertical position: $0 \le g^\mu = y - R \perp di^\mu \ge 0$, where $R$ is the radius of the ball. The dynamic response displays an accumulation phenomena with an infinite number of rebounds over a finite time interval. An event-driven technique would fail to simulate this phenomenon since it is not possible to stop and restart the numerical integration an infinite
<figure><img src="image_6.png" /><figcaption>Figure 2: Bouncing ball: gap distance and zoom.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 3: Bouncing ball: velocity and zoom.</figcaption></figure>

number of times. The physical parameters of this model are as follows: mass $m = 1$ kg, radius $R = 0.2$ m, gravity acceleration $a_g = 10$ m/s$^2$, initial height $h_0 = 1.001$ m. With an initial height of $h_0 = 1$ m, the first impact would be synchronized with a time step; this special case would be numerically favorable but not representative of a general situation; this is the reason why the value $h_0 = 1.001$ m is selected. The numerical solution is computed with a time step $h = 2 \cdot 10^{-3}$ s. The reference solution has been obtained analytically and the accumulation instant is $t = 3.602$ s.

In Figure 2, one observes that the smooth algorithm leads to an erroneous dynamic response characterized by increases in the energy at the impact times. Therefore, this algorithm is not consistent in the presence of impacts and it will not be studied in more detail for this example. The variant of the Nonsmooth-α scheme with the modified LCP condition based on Eq. (27) is then tested. The numerical solution behaves well until the accumulation phenomenon. Then, the constraints at position level (Figure 2) and at velocity level (Figure 3) oscillate between active and inactive states. A similar chattering behaviour was already investigated in [12]. This observation demonstrates that, when the assumption represented by Eq. (28) is not used, inconsistent results can be obtained.
<figure><img src="image_8.png" /><figcaption>Figure 4: Bouncing ball: gap distance and zoom.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 5: Bouncing ball: velocity and energy.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 6: Bouncing ball: impulse of the contact reaction force $\Lambda$ and zoom.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 7: Bouncing ball: Lagrange multiplier $\nu$ and zoom.</figcaption></figure>

As can be seen in Figure 4, the Moreau-Jean scheme and the Nonsmooth-α GGL scheme based on Eq. (28) provide a consistent evaluation of the dynamic response and of the accumulation phenomenon. However, the correct enforcement of the constraint at position level is only ensured by the Nonsmooth-α GGL scheme. Figure 5 shows that the Moreau and Nonsmooth-α GGL algorithms can represent jumps in velocities. The energy of the numerical response is a kind of staircase function, which is almost constant during the free-flight phases (it is actually strictly constant for the Nonsmooth-α GGL responses, but slowly decreasing for the Moreau response) and undergoes sudden drops at the impact times.

The impulse of the contact reaction force $\Lambda$ is depicted in Figure 6. One observes some impulsive contributions at each impact instant and a continuous contribution after the accumulation when the contact is closed. The Lagrange multiplier $\nu$, which is shown in Figure 7, can only take positive values. After the accumulation, $\Lambda$ converges to $ma_g h$, so that $\Lambda/h$ represents the reaction force, in agreement with the discussion in Section 7.1.
<figure><img src="image_12.png" /><figcaption>Figure 8: Bouncing pendulum: gap distance and zoom.</figcaption></figure>

10.2. *The bouncing pendulum*

This example studies the impact of a simple rigid pendulum with both uni-
lateral and bilateral constraints. The pendulum is constrained to swing around
a pivot in the x-y plane. It consists of a massless rod and a concentrated mass
at the tip. The kinematics of the simple pendulum is described in a redundant
manner in terms of three Cartesian coordinates **q** = [x y θ]ᵀ. Since the sys-
tem has one degree-of-freedom, two nonlinear bilateral constraints are expressed
as g̅**U** = [x - l cos θ y - l sin θ]ᵀ = 0 with U = {1, 2}. The pendulum is sub-
jected to the gravity acceleration a₉, and thus the vector of external forces reads:
**f** ext = [0 - ma₉ 0]ᵀ. The mass matrix is **M** = diag [m m J].

As in [16], the physical parameters of the model are the mass *m* = 1 kg, the
moment of inertia *J* = 0.1 kgm², the length *l* = 1 m, the gravity acceleration
*a*<sub>g</sub> = 10 m/s², the initial position **q**<sub>0</sub> = [*l* cos θ<sub>0</sub> *l* sin θ<sub>0</sub> θ<sub>0</sub>]<sup>T</sup>, with θ<sub>0</sub> = π/12
rad. The pendulum is released from rest and swings clockwise due to the gravity
force. An obstacle is placed such that a unilateral constraint is enforced as 0 ≤
g<sup>U</sup> = *x* - √2/2 ⊥ d*i*<sup>U</sup> ≥ 0 with U = {3}. The coefficient of restitution is chosen
as *e* = 0.8. In the numerical tests, the nominal time step size is *h* = 10<sup>-3</sup> s and
the reference solution is obtained using the Nonsmooth-α GGL method and a
smaller time step size *h* = 10<sup>-4</sup> s.

Since impacts occur in the dynamic response, the Smooth-α index-3 algo-
rithm is not considered here. As in the bouncing ball example, we observe in
Figures 8 and 9 that both the Moreau and the Nonsmooth-α GGL method lead
to a consistent representation of jumps in velocities and of the accumulation
phenomenon in the dynamic response. Also, the unilateral constraint at posi-
tion level is exactly satisfied by the Nonsmooth-α GGL scheme, but not by the
Moreau method. In Figure 10, the impulse of the contact reaction force includes
impulsive and continuous contributions. When the contact is closed, after the
accumulation phenomenon, the contact force can be estimated as Λ³/h. This is
consistent with the observation that, in that phase, Λ³ is 10 times smaller for the
reference solution, since it has been obtained using a 10 times smaller time step.
<figure><img src="image_13.png" /><figcaption>Figure 9: Bouncing pendulum: velocity and energy.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 10: Bouncing pendulum: unilateral multiplier $\Lambda^3$ and zoom.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 11: Bouncing pendulum: unilateral multiplier $\nu^3$.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 12: Bouncing pendulum: bilateral constraint at position level.</figcaption></figure>

The Lagrange multiplier $\nu^3$ is shown in Figure 11. The value of $\nu^3$ is smaller for the reference solution, which is consistent with the fact that this multiplier is expected to converge to zero when the time step decreases.

Let us study the bilateral constraints. Figure 12 represents the first bilateral constraint at position level. A drift phenomenon is observed for the Moreau method. In contrast, the bilateral constraint is satisfied exactly at position level (as well as at velocity level) by the Nonsmooth-α GGL method. The Lagrange multipliers of the bilateral constraints are shown in Figure 13. In the Nonsmooth-α GGL scheme, $\Lambda^1$ differs from zero only if the unilateral constraint is active, whereas $\tilde{\lambda}^1$ has permanently a non-zero value. Some spurious transient oscillations are observed for $\tilde{\lambda}^1$ after each impact but they are quickly damped out. These oscillations can be explained by the fact that, between impacts, the dynamic response is computed by solving an index-2 GGL formulation of a DAE with disturbed initial values. Indeed, after the impact, the algorithm only guarantees that the constraints are satisfied at position and velocity levels, but not at acceleration-level. The violation of the hidden constraints at acceleration-level after an impact leads to transient numerical oscillations which are damped out thanks to the presence of numerical dissipation in the generalized-α algorithm, see [26] for a more detailed analysis of this algorithm for DAEs. The importance of these disturbances decreases when the time step is decreased or when the spectral radius is decreased, see Figure 14. For a spectral radius $\rho = 0$, the oscillations are completely eliminated from the numerical solution. Even though these oscillations are not desirable, in our opinion, they can be tolerated in practical applications. Let us analyze the multiplier $\Lambda^1$ obtained using the Moreau-Jean scheme (there is no smooth multiplier $\tilde{\lambda}^1$ in this case). Between the impacts and when the contact is closed, this Lagrange multiplier has a continuous component, but impulses occur at every impact time. It appears that these impulses match quite well the value of $\Lambda^1$ obtained from the Nonsmooth-α GGL scheme, and that the continuous part has some similarities with the value $h\tilde{\lambda}^1$. According to Section 7.1, when the contact is either closed or inactive, it is more appropriate to
<figure><img src="image_17.png" /><figcaption>Figure 13: Bouncing pendulum: bilateral multipliers (first column: $\Lambda^1$, $\tilde{\lambda}^1$ and zoom on $\tilde{\lambda}^1$; second column: $\Lambda^1$, zoom 1 on $\Lambda^1$ and zoom 2 on $\Lambda^1$).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 14: Bouncing pendulum: bilateral multiplier $\tilde{\lambda}^1$ for three different values of the spectral radius in the Nonsmooth-$\alpha$ GGL method.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 15: Bouncing pendulum: bilateral reaction force $\Lambda^1/h + \tilde{\lambda}^1$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 16: Bouncing pendulum: bilateral multiplier $\nu^1$.</figcaption></figure>
compare the reaction force $\Lambda^1/h + \tilde{\lambda}^1$ as shown in Figure 15, which indicates that the Moreau–Jean and the Nonsmooth-α GGL algorithms lead to similar results. Finally, the multiplier associated with the bilateral constraint at position level $\nu^1$ is depicted in Figure 16 and similar observations can be made as in Figure 11.

### 10.3. Horizontal impact of an elastic bar

The impact of an elastic bar on a rigid wall is now analyzed. The bar starts its motion from a distance $d_0$ in an undeformed configuration with a uniform initial velocity field $v_0$, and bounces back against the wall. The impact is horizontal so that no gravity needs to be considered. An analytically-exact solution for this benchmark is detailed in [27]. Before the impact, the bar stays undeformed. Once the bar reaches the wall, a shock wave is caused by the impact and travels with a speed $\sqrt{E/\rho}$ to the other extremity, and then travels back to the contact point. During the travelling of the shock wave, the bar stays in closed contact with the wall due to its compression during a period of time $\Delta t = 2L\sqrt{\rho/E}$. As soon as the wave reaches the contact point, the bar takes off and the contact is released.

The same parameters as in [27] are applied in this test example: Young Modulus $E = 900$ Pa, density $\rho = 1$ kg/m$^3$, undeformed initial length $L = 10$ m, initial distance from the wall $d_0 = 5.005$ m, and initial velocity $v_0 = 10$ m/s. In this case, the bar stays in closed contact during a period of time $\Delta t = 2/3$ s. The numerical responses are compared with an analytical solution.

A restitution coefficient between the flexible rod and the ground needs to be defined. For flexible bodies, the restitution coefficient does not have a clear physical meaning but it should rather be considered as a numerical damping coefficient which dissipates a small amount of energy at the contact point when an impact occurs. For this reason, its value is chosen as $e = 0$.

This one-dimensional problem is discretized using 200 finite elements. The time step size is chosen as $h = 2 \cdot 10^{-3}$ s, so that the Courant number is 1.2.

Let us analyze the numerical response. Figure 17 shows the position at the bottom of the bar. The three methods lead to reasonable results and are able to represent to closed contact during a finite time interval. However, one observes some oscillations for the Smooth-α index-3 algorithm and some penetration for the Moreau method. The velocity and the total energy of the bouncing elastic bar are analyzed in Figure 18. For the Smooth-α index-3 algorithm, the velocity constraint is not exactly satisfied when the contact is closed and the energy evolves in a non-monotonous way, which indicates a lack of consistency. In contrast, the constraint at velocity level is satisfied by the Moreau and the Nonsmooth-α GGL methods when the contact is closed. The Nonsmooth-α GGL method leads to less energy dissipation after the impact than the Moreau method, so that one can conclude that the Nonsmooth-α GGL scheme leads to a better energy representation than the Moreau–Jean scheme. The unilateral Lagrange multipliers $\lambda$ (Smooth-α scheme) and $\Lambda$ (Moreau–Jean and Nonsmooth-α schemes) are depicted in Figure 19 and the contact reaction force can then be derived either as $\lambda$ or as $\Lambda/h$. Here, the Smooth-α method leads to significant spurious numerical
<figure><img src="image_21.png" /><figcaption>Figure 17: Bar impact: gap distance and zoom.</figcaption></figure>

oscillations, while the numerical solution is much more stable for the Moreau–Jean and the Nonsmooth-α GGL method. Therefore, the Moreau–Jean and the Nonsmooth-α schemes lead to a significantly better prediction of the contact force than the Smooth-α method.

### 10.4. Rocking block

The rocking block example shown in Fig. 1(d) is also analyzed in [12]. The configuration is represented using three generalized coordinates $\mathbf{q} = [x \ y \ \theta]^T$, i.e., the horizontal and vertical positions of the center of mass and the orientation of the block. The equations of motion include two nonlinear unilateral constraints

$$ \mathbf{M} = \begin{bmatrix} m & 0 & 0 \\ 0 & m & 0 \\ 0 & 0 & J \end{bmatrix}, \quad \mathbf{f} = \begin{bmatrix} 0 \\ -ma_g \\ 0 \end{bmatrix}, \quad \mathbf{g} = \begin{bmatrix} y - \frac{l}{2} \cos \theta + \frac{L}{2} \sin \theta \\ y - \frac{l}{2} \cos \theta - \frac{L}{2} \sin \theta \end{bmatrix} \quad (53) $$

and the constraint gradient is computed as

$$ \mathbf{g}_q = \begin{bmatrix} 0 & 1 & \frac{l}{2} \sin \theta + \frac{L}{2} \cos \theta \\ 0 & 1 & \frac{l}{2} \sin \theta - \frac{L}{2} \cos \theta \end{bmatrix} \quad (54) $$

The system undergoes multiple impacts and an accumulation phenomenon is observed. The numerical parameters are as in [12], with $m = 1$, $L = 1$, $l = 1.5$,
<figure><img src="image_22.png" /><figcaption>Figure 18: Bar impact: velocity and energy.</figcaption></figure>

$J = m(L^2 + l^2)/12$, $a_g = 9.81$, the coefficient of restitution is $e = 0.5$ for both constraints and the initial values are $q_0 = [0 \ 0 \ 1 \ 0.2]^T$, $v_0 = 0$. The nominal time step is $h = 10^{-2}$ s and the reference solution is obtained using the Nonsmooth-α GGL scheme and a time step $h = 2 \cdot 10^{-4}$ s.

The Moreau scheme and the Nonsmooth-α GGL method provide a satisfactory dynamic response, which is shown in Figures 20-24. Similar observations are made as in the preceding examples. The Nonsmooth-α GGL scheme satisfies the constraints both at position and velocity level and it leads to an improved description of the energy behaviour, with less numerical dissipation than in the Moreau scheme.

## 11. Conclusion

This paper proposes a unified and valid formulation for frictionless contact conditions between rigid bodies as well as between flexible bodies for the dynamic analysis of multibody systems. It can naturally represent a wide spectrum of situations such as finite contact durations, multiple contacts, contact force impulses, velocity jumps and accumulation phenomena. The stability and the robustness of the scheme with respect to the choice of numerical parameters have been illustrated in a number of numerical examples.

Based on a GGL formulation, the algorithm guarantees the exact satisfaction
<figure><img src="image_23.png" /><figcaption>Figure 19: Bar impact: unilateral Lagrange multiplier $\Lambda$ and zoom.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 20: Rocking block: vertical position.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure 21: Rocking block: unilateral constraints and zoom.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 22: Rocking block: time derivative of the unilateral constraints.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Figure 23: Rocking block: energy and zoom.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 24: Rocking block: unilateral Lagrange multiplier $\Lambda^1$ and zoom.</figcaption></figure>
of the complementarity condition at position level, which means that no penetra-
tion is tolerated, and at velocity level. After time discretization, the equations of
motion involve two complementarity conditions and the problem can be solved
efficiently at each time step using a monolithic Newton semi-smooth method.
The complementarity conditions are imposed using an active set method, i.e.,
the activation criteria are evaluated at each Newton iteration in a fully implicit
way. The numerical examples demonstrate the fast and robust convergence of
the Newton semi-smooth procedure.

The algorithm also relies on a splitting of the equations to isolate smooth
terms which are integrated with a higher-order accuracy than the impulsive terms.
Using a generalized-α method for the smooth parts of motion, an accurate descrip-
tion of vibrations phenomena in nonsmooth dynamic systems can be obtained
with a controllable numerical dissipation. For a given spectral radius at infinite
frequencies, the numerical dissipation of the proposed Nonsmooth-α GGL scheme
in the low- and mid-frequency ranges is significantly lower than the numerical dis-
sipation of the Moreau–Jean θ scheme. Therefore, the proposed method leads to
a dynamic response with a qualitatively improved energy behaviour.

As a perspective to this work, the formulation could be tested for large scale
finite element models and it could be extended to frictional contact problems.

References

[1] T. Laursen, V. Chawla, Design of energy conserving algorithms for frictionless dynamic contact problems, International Journal for Numerical Methods in Engineering 40 (1997) 863–886.

[2] V. Chawla, T. Laursen, Energy consistent algorithms for frictional contact problem, International Journal for Numerical Methods in Engineering 42 (1998) 799–827.

[3] Y. Ayyad, M. Barboteu, Formulation and analysis of two energy-consistent methods for nonlinear elastodynamic frictional contact problems, Journal of Computational and Applied Mathematics 228 (2009) 254–269.

[4] D. Stewart, Rigid-body dynamics with friction and impact, SIAM Review 42 (2000) 3–39.

[5] L. Paoli, M. Schatzman, A numerical scheme for impact problems I: The one-dimensional case, SIAM Journal of Numerical Analysis 40 (2002) 702–733.

[6] L. Paoli, M. Schatzman, A numerical scheme for impact problems II: The multi-dimensional case, SIAM Journal of Numerical Analysis 40 (2002) 734–768.
[7] M. Jean, J.-J. Moreau, Dynamics in the presence of unilateral contacts and dry friction: a numerical approach, in: G. Del Pietro, F. Maceri (Eds.), Unilateral problems in structural analysis. II, CISM 304, Springer Verlag, 1987, pp. 151–196.

[8] J.-J. Moreau, Unilateral contact and dry friction in finite freedom dynamics, in: J.-J. Moreau, P. Panagiotopoulos (Eds.), Non-smooth mechanics and applications, volume 302, Springer-Verlag, Wien - New York, 1988, pp. 1–82.

[9] M. Jean, The non-smooth contact dynamics method, Computer Methods in Applied Mechanics and Engineering 177 (1999) 235–257.

[10] R. Krause, M. Walloth, A family of space-time connecting discretization schemes for elastodynamic contact problems with local impact detection, International Journal for Numerical Methods in Engineering 200 (2011) 3425–3438.

[11] C. Gear, B. Leimkuhler, G. Gupta, Automatic integration of Euler-Lagrange equations with constraints, Journal of Computational and Applied Mathematics 12–13 (1985) 77–90.

[12] V. Acary, Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb’s friction, Computer Methods in Applied Mechanics and Engineering 256 (2013) 224–250.

[13] S. Schoeder, H. Ulbrich, T. Schindler, Discussion on the Gear-Gupta-Leimkuhler method for impacting mechanical systems, Multibody Systems Dynamics 31 (2013) 477–495.

[14] M. Hintermüller, K. Ito, K. Kunish, The primal-dual active set strategy as a semismooth Newton method, SIAM Journal on Optimization 13 (2003) 865–888.

[15] I. Ben Gharbia, J. Gilbert, Nonconvergence of the plain Newton-min algorithm for linear complementarity problems with a P-matrix, Mathematical Programming, Series A 134 (2012) 349–364.

[16] Q.-z. Chen, V. Acary, G. Virlez, O. Brüls, A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints, International Journal for Numerical Methods in Engineering 96 (2013) 487–511.

[17] J.-J. Moreau, Numerical aspects of the sweeping process, Computer Methods in Applied Mechanics and Engineering 177 (1999) 329–349.

[18] K. Brenan, S. Campbell, L. Petzold, Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, 2nd ed., SIAM, Philadelphia, 1996.
[19] V. Acary, B. Brogliato, Numerical methods for nonsmooth dynamical systems - Applications in Mechanics and Electronics, volume 35 of Lecture Notes in Applied and Computational Mechanics, Springer-Verlag, Berlin, 2008.

[20] V. Acary, Higher order event capturing time-stepping schemes for nonsmooth multibody systems with unilateral constraints and impacts, Applied Numerical Mathematics 62 (2012) 1259–1275.

[21] M. Arnold, O. Brüls, Convergence of the generalized-α scheme for constrained mechanical systems, Multibody System Dynamics 18 (2007) 185–202.

[22] J. Chung, G. Hulbert, A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-α method, ASME Journal of Applied Mechanics 60 (1993) 371–375.

[23] S. Erlicher, L. Bonaventura, O. Bursi, The analysis of the generalized-α method for non-linear dynamic problems, Computational Mechanics 28 (2002) 83–104.

[24] P. Alart, A. Curnier, A mixed formulation for frictional contact problems prone to Newton like solution methods, Computer Methods in Applied Mechanics and Engineering 92 (1991) 353–375.

[25] C. Bottasso, O. Bauchau, A. Cardona, Time-step-size-independent conditioning and sensitivity to perturbations in the numerical solution of index three differential algebraic equations, SIAM Journal on Scientific Computing 29 (2007) 397–414.

[26] M. Arnold, O. Brüls, A. Cardona, Error analysis of generalized-α Lie group time integration methods for constrained mechanical systems, Numerische Mathematik in press (2014).

[27] D. Doyen, A. Ern, S. Piperno, Time-integration schemes for the finite element dynamic Signorini problem, SIAM Journal on Scientific Computing 33 (2011) 223–249.