# A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints

Qiong Zhong Chen, Vincent Acary, Geoffrey Virlez, Olivier Brüls

► To cite this version:

Qiong Zhong Chen, Vincent Acary, Geoffrey Virlez, Olivier Brüls. A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints. International Journal for Numerical Methods in Engineering, 2013, 96 (8), pp.487-511. 10.1002/nme.4563 . hal-00858066

HAL Id: hal-00858066

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

Submitted on 28 Oct 2017

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints

Qiong-zhong Chen<fnref n="1" />, Vincent Acary<fnref n="2" />, Geoffrey Virlez<fnref n="1" /> and Olivier Brüls<fnref n="1" />,*,†

Mechanical systems are usually subjected not only to bilateral constraints, but also to unilateral constraints. Inspired by the generalized-α time integration method for smooth flexible multibody dynamics, this paper presents a nonsmooth generalized-α method, which allows a consistent treatment of the nonsmooth phenomena induced by unilateral constraints and an accurate description of the structural vibrations during free motions. Both the algorithm and the implementation are illustrated in detail. Numerical example tests are given in the scope of both rigid and flexible body models, taking account for both linear and nonlinear systems and comprising both unilateral and bilateral constraints. The extended nonsmooth generalized-α method is verified through comparison to the traditional Moreau-Jean method and the fully implicit Newmark method. Results show that the nonsmooth generalized-α method benefits from the accuracy and stability property of the classical generalized-α method with controllable numerical damping. In particular, when it comes to the analysis of flexible systems, the nonsmooth generalized-α method shows much better accuracy property than the other two methods.

**KEY WORDS:** nonsmooth dynamics; unilateral constraints; generalized-α method; time integration scheme; flexible multibody systems

## 1. INTRODUCTION

The equations of motion for a mechanical system with bilateral constraints can be characterized by a mixed set of second order differential and algebraic equations. One of the classical research interests is in the efficient numerical solutions. Among the classical numerical solvers, the Hilber-Hughes-Taylor (HHT) method [1] and the generalized-α method [2], which are extended from the original implicit Newmark approach [3], are popular and have been extensively applied in structural dynamics. With an optimal choice of parameters, these solvers yield second-order accuracy and unconditional stability, which are of utmost importance because the structural dynamics equations can be very stiff with a broad range of system eigenfrequencies [4]. Cardona and Géradin [5] proposed to apply these solvers to index-3 DAEs in flexible multibody dynamics. These authors have demonstrated that, in the linear regime, the numerical instability in the Lagrangian multipliers can be removed with a small amount of numerical dissipation at high frequencies [4]. An extended analysis on the global convergence can be found in [6].

In practice, mechanical systems are usually subjected not only to bilateral constraints, but also to unilateral constraints. For instance, the global dynamic behaviour of wind turbines is characterized

<footnote n="1">¹Department of Aerospace and Mechanical Engineering (LTAS), University of Liège Chemin des Chevreuil, 1 (B52/3), 4000 Liège, Belgium</footnote>

<footnote n="2">²INRIA Rhône-Alpes, Centre de recherche Grenoble, 655 avenue de l'Europe, Inovallée de Montbonnot, 38334 St Ismier Cedex, France</footnote>

<footnote>*Correspondence to: Olivier Brüls, Department of Aerospace and Mechanical Engineering (LTAS), University of Liège Chemin des Chevreuil, 1 (B52/3), 4000 Liège, Belgium.</footnote>

<footnote>†E-mail: o.bruls@ulg.ac.be</footnote>
by vibrations of structural elements such as the blades and the tower coupled with strong nonlinear-
ities in the transmission line due to contacts between stiff or rigid bodies with backlash and friction
phenomena. Unilateral constraints can be modelled either by regularization techniques based on
single-valued force laws, or by set-valued force laws that lead to nonsmooth models. Regularization
techniques aim at replacing the discontinuity by a stiff but compliant model, for example, the con-
tinuous impact models proposed in [7]. However, the knowledge of a number of physical contact
parameters such as contact stiffness and damping, which are difficult to be obtained by experi-
ments, can be required [7], and these techniques result in stiff differential equations and inexact
solutions [8]. Nonsmooth models, represented as differential inclusions, complementarity systems
or variational inequalities [9], are used in order to model the unilateral nature of the constraints, usu-
ally modelled by the Signorini condition. Alternatively, nonsmooth models can also express impact
laws, required by the unilateral constraints in finite-dimensional systems (rigid models or space-
discretized systems). In the sequel, Newton's impact law will be used. When nonsmooth models are
involved, one of the main computational challenges remains the numerical time integration of the
dynamics of systems, which may also be nonsmooth with impacts and jumps in velocities.

In the context of rigid multibody systems, occurrences of impacts and/or abrupt changes in
the velocities due to Signorini's condition and Coulomb's friction require great cares in the time
integration method in order to ensure consistency and stability. There exist two main groups of
numerical integration schemes for nonsmooth systems, namely, event-driven schemes and time-
stepping schemes. Event-driven schemes are based on an accurate 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,
but fail to handle frequent transitions in a short time and no convergence proof has been provided
[9–11]. In practice, they are only suitable for small multibody systems with a limited number of
events. 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 schemes
have been designed up to now : the Schatzman–Paoli scheme [12, 13] is based on central difference
scheme and the Moreau–Jean scheme [14–16] is based on a $\theta$-method [17]. Time-stepping methods
have been proven to be convergent and robust, and are extensively applied as the solution to the non-
smooth 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 [10, 18]. Thus, the accuracy is less
satisfactory unless a very small step size is applied. However, it remains robust and efficient even for
a large number of events as it is usual in structural dynamics. Attempts have recently been made to
improve the global order of accuracy of time-stepping schemes with nonsmooth events [10, 19, 20].
But for the moment, the application to nonsmooth mechanical systems of the classical time integra-
tion schemes for structural and flexible multibody systems, that is, HHT method and generalized-$\alpha$
method, has not been successfully carried out. Within the time-stepping strategy, it remains difficult
to expect to increase the order of accuracy of the scheme when an impact occurs. Nevertheless,
the goal of this paper is to bridge this gap in order to benefit from the qualitative properties of the
generalized-$\alpha$ scheme and the HHT scheme, in terms of stability and controlled numerical damping.

In the context of computational contact mechanics of solids and structures, a lot of variants of
the Newmark scheme have been proposed to solve the contact dynamics problems, see [21, 22] for
two recent reviews. Mostly, it is implicitly assumed in these references that the solutions includ-
ing positions, displacements, velocities and contact forces are sufficiently smooth so that Newmark
schemes of order two are applicable. Nevertheless, the lack of robustness and the numerical insta-
bilities of the midpoint, trapezoidal rule and even dissipative HHT schemes, often leading to the
blow-up of the numerical computation, have motivated the development of alternatives and modifi-
cations to standard integration schemes. To cite a few of them, the invocation of energy-conserving
and dissipative schemes was the first attempt to remedy this problem [23–29]. It generally results in
a scheme of lower order where the contact forces or impulses are treated in a fully implicit manner
[26, 27], and where the constraint is expressed at the velocity level [24, 25, 29]. It should be noted
that these requirements exactly meet the principles of design of the Moreau–Jean time-stepping
schemes [14–16], but in the latter case, they are motivated by the nonsmoothness of solutions. By the
way, the reader will find in [30] a comparison between the Moreau–Jean scheme and the Newmark
scheme for space-discretized impacting bars. The convergence of the Newmark scheme with central differences and with an implicit treatment of the contact constraint at the position level has been proven in [31] for vibrating beams between two stops. Fortunately enough, it has been shown in [32] that the question of the coefficient of restitution is not raised on this case. The implicit treatment of the contact forces was also proposed in the context of explicit integrators in a seminal paper [33]. Other attempts have been based on (1) a mass redistribution, which consists in removing the mass from the contact boundaries [34, 35] or (2) a contact stabilization of the relative velocities at contact [36]. In our case, the contact activation between discrete mechanical models induces the nonsmoothness of the solutions and, because rigid bodies are involved, it is difficult to remove masses and to stabilize contact without resorting to an impact law. Therefore, the direct application of higher order schemes in this context may be hazardous.

Within the framework of structural dynamics with rigid and flexible bodies, this paper proposes to take the best of the generalized-α method and the Moreau–Jean scheme by

* keeping the integration rule of the generalized-α scheme for the nonlinear smooth dynamics; the aim is to enjoy the stability, the robustness and the controlled numerical dissipation at high frequencies;

* dealing, as in the Moreau–Jean scheme, with the contact forces and impacts through their associated impulses; in this way, the consistency of the evaluation of the contact reactions is ensured when the time-step vanishes and an impact occurs;

* treating, as in the Moreau–Jean scheme, the unilateral constraint at the velocity level together with the impact law; this yields an energetically consistent treatment of the jumps in velocity together with a correct stabilization of the velocity at contact.

The proposed method aims at providing a general solver for structural dynamics, which can take both smooth and nonsmooth mechanics into account. Even though the global accuracy remains order one because of the nonsmoothness [8, 9] or the presence of constraint forces, it will improve the control of high-frequency numerical damping in comparison with the Moreau–Jean time-stepping scheme based on the θ-method. Stability and convergence will be empirically analyzed and the properties of the proposed method will be illustrated based on academic examples.

This paper is organized as follows. Section 2 describes the theory of dynamical equations of a flexible multibody system with both bilateral and unilateral constraints. Based on the equations of motion, Section 3 presents the formulation of a time stepping method, herewith called the nonsmooth generalized-α method. Numerical implementation of this method is introduced in Section 4. In Section 5, numerical examples, consisting of impacts of both rigid and flexible multibody systems, ranging from single unilateral constraint to multi unilateral constraints, are given for analysis. Finally, conclusions are presented in Section 6.

## 2. EQUATIONS OF MOTION

The equations of a flexible multibody system including bilateral and unilateral constraints can be expressed in the following form:

$$ M(\mathbf{q}) \, d\mathbf{v} + \mathbf{\Phi}_{\mathbf{q}}^{T}(\mathbf{q}) \, d\mathbf{i}^{b} = \mathbf{f}(\mathbf{q}, \mathbf{v}, t) \, d\mathbf{t} + \mathbf{g}_{\mathbf{q}}^{T}(\mathbf{q}) \, d\mathbf{i}^{u}, \quad (1) $$

$$ \mathbf{\Phi}(\mathbf{q}) = \mathbf{0}, \quad (2) $$

$$ \mathbf{0} \le \mathbf{g}(\mathbf{q}) \perp \mathbf{d}\mathbf{i}^{u} \ge \mathbf{0}, \quad (3) $$

where

* $\mathbf{q}$ is the vector of generalized coordinates, for example, absolute nodal coordinates;

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

* $M$ is the mass matrix, which is assumed to be non-singular;
* $\Phi$ is the vector of bilateral holonomic constraints, and $\Phi_q$ is the matrix of constraint gradients;

* $g$ is the vector of unilateral holonomic constraints, and $g_q$ is the corresponding matrix of constraint gradients;

* $dv$ is the differential measure associated with the velocity $v$ assumed to be of bounded variation;

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

* $di^u$ and $di^b$ are the impulse measure of the contact reaction and the bilateral force, respectively.

Equation (3) is formulated as a complementarity condition, that is, $g^T(q)di^u = 0$ with both $g(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 as $g_qv(t^+) = -e g_qv(t^-)$, where $e$ is the coefficient of restitution, and $v(t^-)$ and $v(t^+)$ denote the velocity before and after the impact, respectively. The complementarity condition can then be expressed at the velocity level in the following form:

$$ \text{if } g(q) \le 0 \text{ then } 0 \le g_q v(t^+) + e g_q v(t^-) \perp di^u \ge 0. \quad (4) $$

This reformulation means that if a unilateral holonomic constraint is active, the complementarity condition is expressed with respect to the unilateral constraint at the velocity level.

Because the motion might be nonsmooth, jumps in the velocity and the corresponding impacts are expected yielding the following decomposition of the measures$^{‡}$:

$$ dv = \dot{v} dt + \sum_i (v(t_i^+) - v(t_i^-)) \delta_{t_i}, \quad (5) $$

$$ di^u = \lambda^u dt + \sum_i p_i^u \delta_{t_i}, \quad (6) $$

$$ di^b = \lambda^b dt + \sum_i p_i^b \delta_{t_i}, \quad (7) $$

where $\lambda^u$ and $\lambda^b$ are the vectors of Lagrange multipliers associated with Lebesgue measurable uni- lateral and bilateral forces, respectively; $(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$; $p_i^u$ and $p_i^b$ are impulses. Inserting Equations (5-7) in (1-3) yields the standard equations of motion for almost every time $t$:

$$ M\dot{v} + \Phi_q^T \lambda^b = f + g_q^T \lambda^u, \quad (8) $$

$$ \Phi(q) = 0, \quad (9) $$

$$ 0 \le g(q) \perp \lambda^u \ge 0, \quad (10) $$

and the impact equations at time $t_i$

$$ M(v(t_i^+) - v(t_i^-)) + \Phi_q^T p_i^b = g_q^T p_i^u, \quad (11) $$

$$ \Phi_q(q)v(t_i^+) = 0, \quad (12) $$

$$ 0 \le g_q v(t_i^+) + e g_q v(t_i^-) \perp p_i^u \ge 0. \quad (13) $$

<footnote>‡we neglect singular measures</footnote>
### 3. TIME INTEGRATION METHOD

In order to present the time integration method, we propose to split the equation in a slightly different way by separating the contribution of the contact reaction and bilateral force from the contribution of smooth forces such that

$$d\mathbf{v} = \tilde{\mathbf{v}}dt + d\mathbf{w}, \qquad (14)$$

with

$$\begin{aligned} \mathbf{M}\tilde{\mathbf{v}} &= \mathbf{f}, \\ \mathbf{M}d\mathbf{w} + \mathbf{\Phi}_q^T d\mathbf{i}^b &= \mathbf{g}_q^T d\mathbf{i}^u, \end{aligned} \qquad (15)$$

where $\tilde{\mathbf{v}}$ is an acceleration variable related with the smooth force $\mathbf{f}$, and $d\mathbf{w}$ is a differential measure because of the contact reaction and bilateral force. A key idea of the proposed method is to integrate the variable $\tilde{\mathbf{v}}$ using a second-order accurate generalized-$\alpha$ method, whereas the differential measure $d\mathbf{w}$ is treated using a time-stepping method based on the first-order Euler implicit scheme. If the constraints are furthermore expressed at the velocity level, this combination ensures a consistent treatment of nonsmooth dynamic effects and an accurate description of structural vibrations during free motions. The discrete evaluation of the first equation in (15) will be denoted as

$$\mathbf{M}_{n+1}\tilde{\mathbf{v}}_{n+1} = \mathbf{f}_{n+1}, \qquad (16)$$

where $\mathbf{M}_{n+1} = \mathbf{M}(\mathbf{q}_{n+1})$ and $\mathbf{f}_{n+1} = \mathbf{f}(\mathbf{q}_{n+1}, \mathbf{v}_{n+1}, t_{n+1})$. For the other part, we keep the values of the impulse over the time step $(t_n, t_{n+1}]$ as primary unknowns, that is,

$$\mathbf{M}_{n+1}\mathbf{W}_{n+1} + \mathbf{\Phi}_{q,n+1}^T \mathbf{\Lambda}_{n+1}^b = \mathbf{g}_{q,n+1}^T \mathbf{\Lambda}_{n+1}^u, \qquad (17)$$

where $\mathbf{W}_{n+1}$ represents a velocity jump over the time step; $\mathbf{\Lambda}_{n+1}^b$ and $\mathbf{\Lambda}_{n+1}^u$ are the Lagrange multipliers linked to bilateral force and unilateral contact reaction respectively. Equation (17) is actually a first-order approximation of the integral of the second equation in (15), that is,

$$\int_{[t_n, t_{n+1}]} (\mathbf{M} d\mathbf{w} + \mathbf{\Phi}_q^T d\mathbf{i}^b) = \int_{[t_n, t_{n+1}]} \mathbf{g}_q^T d\mathbf{i}^u. \qquad (18)$$

The discretization of the contact law is written as

$$\text{if } \bar{\mathbf{g}} \le \mathbf{0} \text{ then } \mathbf{0} \le \mathbf{g}_{q,n+1} \mathbf{v}_{n+1} + e \mathbf{g}_{q,n} \mathbf{v}_n \perp \mathbf{\Lambda}_{n+1}^u \ge \mathbf{0}, \qquad (19)$$

where $\bar{\mathbf{g}}$ is defined either as $\bar{\mathbf{g}} = \mathbf{g}_n$ or as $\bar{\mathbf{g}} = \mathbf{g}_{n+1}^p$, with $\mathbf{g}_{n+1}^p$ defined as the prediction of the unilateral constraints at time $t_{n+1}$. Equation (19) can be rewritten in an equation form as

$$\text{if } \bar{\mathbf{g}} \le \mathbf{0} \text{ then } \Psi (\mathbf{\Lambda}_{n+1}^u, \mathbf{g}_{q,n+1} \mathbf{v}_{n+1} + e \mathbf{g}_{q,n} \mathbf{v}_n) = \mathbf{0}, \qquad (20)$$

with $\Psi(a, b) = a - \max(0, a - \kappa b)$, $\kappa > 0$. Defining the set of active constraints $\mathbf{g}^A$ as the set of all the unilateral constraints satisfying $\bar{\mathbf{g}} \le \mathbf{0}$ and $\mathbf{\Lambda}^A$ as the corresponding Lagrange multipliers, the discretized form of the equations of motion becomes

$$\mathbf{M}_{n+1}\tilde{\mathbf{v}}_{n+1} = \mathbf{f}_{n+1}, \qquad (21)$$

$$\mathbf{M}_{n+1}\mathbf{W}_{n+1} + \mathbf{\Phi}_{q,n+1}^T \mathbf{\Lambda}_{n+1}^b = (\mathbf{g}_{q,n+1}^A)^T \mathbf{\Lambda}_{n+1}^A, \qquad (22)$$

$$\mathbf{0} = \mathbf{\Phi}_{q,n+1} \mathbf{v}_{n+1}, \qquad (23)$$

$$\mathbf{0} = \Psi (\mathbf{\Lambda}_{n+1}^A, \mathbf{g}_{q,n+1}^A \mathbf{v}_{n+1} + e \mathbf{g}_{q,n}^A \mathbf{v}_n). \qquad (24)$$
In this expression, the bilateral constraints are formulated at the velocity level. A consequence of this index-2 formulation is the drift-off phenomena. To deal with the drift of the constraints, projection methods can be applied if needed, see [17, 37].

For smooth dynamics, the classical generalized-α algorithm can be obtained by introducing a vector **a** of acceleration-like variables in the Newmark formulae [2]. The nonsmooth generalized-α method is obtained by further taking the impulsive terms into account so that it takes the following form:

$$ \tilde{\mathbf{q}}_{n+1} = \mathbf{q}_n + h\mathbf{v}_n + h^2\left(\frac{1}{2} - \beta\right)\mathbf{a}_n + h^2\beta\mathbf{a}_{n+1}, \quad (25) $$

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

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

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

$$ \mathbf{q}_{n+1} = \tilde{\mathbf{q}}_{n+1} + \frac{1}{2}h\mathbf{W}_{n+1}, \quad (29) $$

where $\tilde{\mathbf{q}}, \tilde{\mathbf{v}}$ are intermediate variables representing some smooth positions and velocities respectively, whereas $\mathbf{q}, \mathbf{v}$ stand for the generalized position and velocity in the general sense; $\mathbf{a}$ is an intermediate acceleration-like variable defined as in (27) and initialized at the initial time as $\mathbf{a}_0 = \tilde{\mathbf{v}}_0$; the time step size is $h = t_{n+1} - t_n$.

The form of Equations (28–29) is justified as follows. Suppose that an impact is detected at time $t_i \in [t_n, t_{n+1})$, the impact reaction then occurs only in the atomic measure at the instant $t_i$. The influence of the impact is then considered only at the current time step and not at any other time steps. The impulsive increments to the velocity and position are then derived as:

$$ \mathbf{v}_{n+1} - \tilde{\mathbf{v}}_{n+1} = \int_{[t_n, t_{n+1})} \mathrm{d}\mathbf{w} = \mathbf{W}_{n+1}, \quad (30) $$

$$ \mathbf{q}_{n+1} - \tilde{\mathbf{q}}_{n+1} = \int_{[t_n, t_{n+1})} \int_{[t_n, t)} \mathrm{d}\mathbf{w} \mathrm{d}t = \frac{1}{2} h \mathbf{W}_{n+1}, \quad (31) $$

where $\mathbf{W}_{n+1}$ also satisfies the discrete equations of motion in (21–24).

In the unconstrained case for the smooth generalized-α algorithm, the condition of second-order accuracy yields:

$$ \gamma = \frac{1}{2} + \alpha_f - \alpha_m. $$

A small amount of numerical damping, represented by the spectral radius of the algorithm at the high-frequency limit $\rho_\infty$ with $\rho_\infty \in [0, 1]$ can be introduced. The optimal parameters for the generalized-α method are then derived as [2]:

$$ \alpha_m = \frac{2\rho_\infty - 1}{\rho_\infty + 1}, \quad \alpha_f = \frac{\rho_\infty}{\rho_\infty + 1}, \quad \beta = \frac{1}{4} \left( \gamma + \frac{1}{2} \right)^2. $$

For the HHT algorithm, $\alpha_m = 0$ and $\alpha_f \in [0, 1/3]$, whereas for the Newmark algorithm, $\alpha_f = \alpha_m = 0$.

The proposed method will be compared with the Moreau–Jean θ-method, which is obtained by replacing Equations (25–29) by the following equations

$$ \tilde{\mathbf{q}}_{n+1} = \mathbf{q}_n + h(1 - \theta)\mathbf{v}_n + h\theta\tilde{\mathbf{v}}_{n+1}, \quad (32) $$
<figure><img src="image_2.png" /><figcaption>Figure 1. Spectral properties of time integration schemes.</figcaption></figure>

$$ \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 (33) $$

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

$$ \mathbf{q}_{n+1} = \tilde{\mathbf{q}}_{n+1} + h\theta\mathbf{W}_{n+1}. \quad (35) $$

with $\theta \in [0, 1]$. The Moreau-Jean method with $\theta = 1/2$ is equivalent to the nonsmooth generalized-$\alpha$ scheme with $\alpha_m = \alpha_f = 0$, $\beta = 1/4$, $\gamma = 1/2$ and $\mathbf{a}_{n+1} = \tilde{\mathbf{v}}_{n+1}$.

In order to compare the numerical damping of the proposed nonsmooth time integration algorithms, which will be tested in Section 5, Figure 1 shows the spectral radius as a function of the frequency parameter $\omega h$. It is worth noting that these spectral properties are obtained in case of unconstrained motion. The spectral analysis is made for a linear test problem $\ddot{q} + \omega^2 q = 0$. The algorithmic parameters are chosen as follows: for the generalized-$\alpha$ method, the spectral radius at infinity is chosen as $\rho_\infty = 0.8$; for the Moreau-Jean method, $\theta = 1$, which implies the implicit Euler method; for the so-called fully implicit Newmark method, $\gamma = 1$ and $\beta = 1/4(\gamma + 1/2)^2 = 0.5625$. As one can see from Figure 1, the generalized-$\alpha$ method has the least numerical dissipation. In contrast, the Moreau-Jean method is the most numerically dissipative one.

# 4. IMPLEMENTATION OF NONSMOOTH GENERALIZED-$\alpha$ ALGORITHM

## 4.1. A generalized Newton method

Formulae (25–29) can be solved together with the discretized equations of motion (21–24) using a predictor-corrector scheme coupled with a solver for linear complementarity problems. Let us treat $\tilde{\mathbf{v}}$, $\mathbf{v}$, $\mathbf{\Lambda}^b$ and $\mathbf{\Lambda}^A$ as independent sets of unknowns. The linearized equation for Newton iteration becomes:

$$ \begin{bmatrix} \gamma_t \mathbf{M} + \beta_t \mathbf{K}_t & \mathbf{C}_t + \frac{1}{2} h \mathbf{K}_t & \mathbf{0} \\ -\mathbf{M} + \beta_t \mathbf{K}_t^{\text{con}} & \mathbf{M} + \frac{1}{2} h \mathbf{K}_t^{\text{con}} & \mathbf{\Phi}_q^T \\ \beta_t \mathbf{\Phi}_{q2} & \mathbf{\Phi}_q + \frac{1}{2} h \mathbf{\Phi}_{q2} & \mathbf{0} \\ \hline \beta_t \mathbf{\Psi}_q & \mathbf{\Psi}_v + \frac{1}{2} h \mathbf{\Psi}_q & \mathbf{0} \end{bmatrix} \begin{bmatrix} \Delta \tilde{\mathbf{v}} \\ \Delta \mathbf{v} \\ \frac{\Delta \mathbf{\Lambda}^b}{\Delta \mathbf{\Lambda}^A} \end{bmatrix} = \begin{bmatrix} -\mathrm{res}^q \\ -\mathrm{res}^w \\ -\mathrm{res}^\Phi \\ -\mathrm{res}^\Psi \end{bmatrix}, \quad (36) $$

where

* the numerical parameters: $\gamma_t = (1 - \alpha_m)/[(1 - \alpha_f)\gamma h]$ and $\beta_t = \beta h/\gamma - 1/2 h$;

* the tangent stiffness and damping matrices are respectively defined as: $\mathbf{K}_t = \partial(\mathbf{M}\tilde{\mathbf{v}} - \mathbf{h})/\partial\mathbf{q}$ and $\mathbf{C}_t = -\partial\mathbf{h}/\partial\mathbf{v}$;

* $\mathbf{K}_t^{\text{con}}$ is defined as $\mathbf{K}_t^{\text{con}} = \partial(\mathbf{M}\mathbf{W} + \mathbf{\Phi}_q^T \mathbf{\Lambda}^b - (\mathbf{g}_q^A)^T \mathbf{\Lambda}^A)/\partial\mathbf{q}$;
• $\Phi_{q2}$ is the derivative of $\Phi_{q,n+1}v_{n+1}$ with respect to $q_{n+1}$;

• $\Psi_q, \Psi_v$ and $\Psi_{\Lambda^A}$ are the partial derivatives of the function $\Psi(\Lambda_{n+1}^A, g_{q,n+1}^A v_{n+1} + e g_{q,n}^A v_n)$ with respect to $q_{n+1}, v_{n+1}$ and $\Lambda_{n+1}^A$ respectively and are indeed the generalized Jacobian matrices at the point of discontinuity;

• $\text{res}^q, \text{res}^{\text{con}}, \text{res}^{\Phi}$ and $\text{res}^{\Psi}$ denote respectively the residuals of the four sets of equations in (21–24). Using a superscript '*' to denote the value of a variable at Newton iteration $k$, the residuals are

$$ \begin{bmatrix} \text{res}^q \\ \text{res}^W \\ \text{res}^\Phi \\ \text{res}^\Psi \end{bmatrix} = \begin{bmatrix} M_{n+1}^* \tilde{v}_{n+1}^* - h_{n+1}^* \\ M_{n+1}^* W_{n+1}^* + (\Phi_{q,n+1}^*)^T \Lambda_{n+1}^{b*} - (g_{q,n+1}^{A*})^T \Lambda_{n+1}^{A*} \\ \Phi_{q,n+1}^* v_{n+1}^* \\ \Psi_{n+1}^* \end{bmatrix}. \quad (37) $$

The detailed computation of generalized Jacobians of $\Psi$ can be found in [38, 39] together with some generalized Newton method that allows to solve the contact problem directly by iteratively solving Equation (36). The extension to the frictional case is then straightforward. In the next subsection, we propose a simpler approach to the frictionless case based on a sequence of linear complementarity problems (LCPs) that are solved at each Newton iteration.

## 4.2. Reduction to a sequential linear complementarity problems

In this section, we propose to reformulate the iteration in (36) as a LCP, see [9, 40, 41], which replaces the linearized system that we have to solve at each Newton’s step. A LCP solver is able to find a vector $z$ given a matrix $\bar{M}$ and a vector $\bar{r}$ such that $0 \le \bar{M}z + \bar{r} \perp z \ge 0$.

If we define $\Delta\mathbf{x}^T = [\Delta\tilde{\mathbf{v}}^T \Delta\mathbf{v}^T (\Delta\mathbf{\Lambda}^b)^T]$, the linearized form of equations (21–23) for iteration $k+1$ can be recast as

$$ S_t \Delta \mathbf{x} = -\mathbf{res} + \mathbf{b} (\mathbf{\Lambda}_{n+1}^A - \mathbf{\Lambda}_{n+1}^{A*}) \quad (38) $$

with

$$ S_t = \begin{bmatrix} \gamma_t \mathbf{M}^* + \beta_t \mathbf{K}_t^* & \mathbf{C}_t^* + 1/2 h \mathbf{K}_t^* & \mathbf{0} \\ -\mathbf{M}^* + \beta_t \mathbf{K}_t^{\text{con}*} & \mathbf{M}^* + 1/2 h \mathbf{K}_t^{\text{con}*} & \mathbf{\Phi}_q^{*T} \\ \beta_t \mathbf{\Phi}_{q2}^* & \mathbf{\Phi}_q^* + 1/2 h \mathbf{\Phi}_{q2}^* & \mathbf{0} \end{bmatrix}, \quad (39) $$

$$ \text{res} = \begin{bmatrix} \text{res}^q \\ \text{res}^W \\ \text{res}^\Phi \end{bmatrix} \quad \text{and} \quad \mathbf{b} = \begin{bmatrix} \mathbf{0} \\ (\mathbf{g}_q^{A*})^T \\ \mathbf{0} \end{bmatrix}. $$

The unilateral constraint at the velocity level is also linearized as

$$ \mathbf{0} \le \mathbf{g}_{q,n+1}^{A*} \mathbf{v}_{n+1}^* + e \mathbf{g}_{q,n}^A \mathbf{v}_n + \mathbf{c}^T \Delta \mathbf{x} \quad (40) $$

with

$$ \mathbf{c}^T = [\beta_t \mathbf{g}_{q2,n+1}^{A*} \quad \mathbf{g}_{q,n+1}^{A*} + 1/2 h \mathbf{g}_{q2,n+1}^{A*} \quad \mathbf{0}], $$

where $\mathbf{g}_{q2,n+1}^A$ is defined as $\mathbf{g}_{q2,n+1}^A = \partial(\mathbf{g}_{q,n+1}^A \mathbf{v}_{n+1} + e \mathbf{g}_{q,n}^A \mathbf{v}_n) / \partial \mathbf{q}_{n+1}$.

Therefore, by substituting Equation (38) into (40), the complementarity condition in Equation (24) is obtained in the standard form as

$$ \mathbf{0} \le \bar{\mathbf{M}} \mathbf{\Lambda}_{n+1}^A + \bar{\mathbf{r}} \perp \mathbf{\Lambda}_{n+1}^A \ge \mathbf{0} \quad (41) $$
with

$$
\begin{cases}
\bar{M} = c^T S_t^{-1} b \\
\bar{r} = g_{q,n+1}^{A*} v_{n+1}^* - c^T S_t^{-1} \text{res} - c^T S_t^{-1} b \Lambda_{n+1}^{A*} + e g_{q,n}^A v_n
\end{cases}
$$

To sum up, the flowchart of the predictor-corrector iterative method for the nonsmooth generalized-α method based on a sequential LCP is shown in Algorithm 1.

<caption>**Algorithm 1** Nonsmooth generalized-α time integration scheme</caption>

Initial values $\mathbf{q}_0$ and $\mathbf{v}_0$
$\tilde{\mathbf{v}}_0 := \mathbf{M}^{-1}\mathbf{h}(\mathbf{q}_0, \mathbf{v}_0, t_0)$
$\mathbf{a}_0 := \tilde{\mathbf{v}}_0$
**for** $n = 0$ to $n_{\text{final}} - 1$ **do**
&emsp; $\tilde{\mathbf{v}}_{n+1} := \mathbf{0}$
&emsp; $\mathbf{a}_{n+1} := 1/(1-\alpha_m) (\alpha_f \tilde{\mathbf{v}}_n - \alpha_m \mathbf{a}_n)$
&emsp; $\tilde{\mathbf{v}}_{n+1} := \mathbf{v}_n + h(1-\gamma)\mathbf{a}_n + \gamma\mathbf{a}_{n+1}$
&emsp; $\mathbf{v}_{n+1} := \tilde{\mathbf{v}}_{n+1}$
&emsp; $\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}$
&emsp; $\mathbf{\Lambda}_{n+1}^b := \mathbf{0}$
&emsp; Compute $\bar{\mathbf{g}}_{n+1}$, find out the index of the active set: $A = \text{find } (\bar{\mathbf{g}}_{n+1} \le 0)$
&emsp; $\mathbf{\Lambda}_{n+1}^A := \mathbf{0}$
&emsp; **for** $i = 1$ to $i_{\text{max}}$ **do**
&emsp; Compute the residual $\mathbf{res}$ and $\mathbf{\Psi}$
&emsp; **if** $\sqrt{\|\mathbf{res}\|^2 + \|\mathbf{\Psi}\|^2} < \text{tol}$ **then**
&emsp;&emsp; break
&emsp; **end if**
&emsp; $\mathbf{\Lambda}_{n+1}^{A*} := \mathbf{\Lambda}_{n+1}^A$
&emsp; Compute $S_t^{-1}\mathbf{b}$ and $S_t^{-1}\mathbf{res}$
&emsp; Compute $\bar{\mathbf{M}}$ and $\bar{\mathbf{r}}$
&emsp; $\mathbf{\Lambda}_{n+1}^A := \text{LCP solver}(\bar{\mathbf{M}}, \bar{\mathbf{r}})$
&emsp; $\Delta\mathbf{x} := -S_t^{-1}\mathbf{res} + S_t^{-1}\mathbf{b} (\mathbf{\Lambda}_{n+1}^A - \mathbf{\Lambda}_{n+1}^{A*})$
&emsp; Extract the increments $\Delta\tilde{\mathbf{v}}$, $\Delta\mathbf{v}$ and $\Delta\mathbf{\Lambda}^b$ from $\Delta\mathbf{x}$
&emsp; $\tilde{\mathbf{v}}_{n+1} := \tilde{\mathbf{v}}_{n+1} + \Delta\tilde{\mathbf{v}}$
&emsp; $\mathbf{v}_{n+1} := \mathbf{v}_{n+1} + \Delta\mathbf{v}$
&emsp; $\mathbf{q}_{n+1} := \mathbf{q}_{n+1} + \beta_t \Delta\tilde{\mathbf{v}} + 1/2 h \Delta\mathbf{v}$
&emsp; $\tilde{\mathbf{v}}_{n+1} := \tilde{\mathbf{v}}_{n+1} + \gamma_t \Delta\tilde{\mathbf{v}}$
&emsp; $\mathbf{\Lambda}_{n+1}^b := \mathbf{\Lambda}_{n+1}^b + \Delta\mathbf{\Lambda}^b$
&emsp; **end for**
&emsp; $\mathbf{a}_{n+1} := \mathbf{a}_{n+1} + (1-\alpha_f)/(1-\alpha_m)\tilde{\mathbf{v}}_{n+1}$
&emsp; **end for**

5. NUMERICAL EXAMPLES

Examples are given in this section to test the numerical methods. In order to highlight the properties
of the presented time integration method, the examples comprise of both rigid-body systems and
discretized flexible systems, as shown in Figure 2. The first three examples, namely, a bouncing
ball, a bouncing oscillator and a bouncing elastic bar, have a linear behaviour in free motion and can
be compared with exact solutions. Besides that, two pendulum examples are also studied: one is the
<figure><img src="image_3.png" /><figcaption>Figure 2. Examples: (a) bouncing ball; (b) linear vertical oscillator; (c) bouncing of an elastic bar; (d) impact of a flexible rotating beam; (e) impact of a rigid pendulum; (f) multiple impacts of a crank-slider mechanism.</figcaption></figure>

impact of a flexible rotating beam with geometric nonlinearities; the other is the impact of a simple rigid pendulum with the simultaneous presence of unilateral and bilateral constraints. The application is finally extended to an example of a crank-slider mechanism with a translational clearance joint. Multiple impacts encountered in the crank-slider mechanism composed of only rigid bodies or a mix of rigid and flexible bodies are then studied.

### 5.1. Bouncing ball

The first example is a standard bouncing ball on a rigid plane, as depicted in Figure 2(a). The bouncing ball is a well-known benchmark with one degree of freedom (DOF) in the field of nonsmooth mechanics. As usual, both the ball and the plane are rigid. It starts from rest, falls because of the gravity and rebounds against the rigid plane with a restitution coefficient of $e = 0.8$. Thus, it introduces a unilateral constraint on the vertical position. The dynamic response includes an accumulation phenomenon with an infinite number of rebounds over a finite time interval. An event-capturing technique would fail to simulate this phenomenon because it is not possible to stop and restart the numerical integration an infinite number of times. Other physical parameters of this model are
as follows: mass $m = 1kg$, radius $R = 0.2m$, gravity acceleration $g = 10m/s^2$, initial height $h_0 = 1.001m$ (with an initial height of $h_0 = 1m$, the first impact would be synchronized with a time step; this special case would be numerically favourable but not representative of a general situation).

The simulation is first carried out using three valid time integration algorithms, namely, the nonsmooth generalized-$\alpha$, the Moreau-Jean and the fully implicit Newmark schemes. The numerical parameters are set as: nominal time step $h = 10^{-3}$s for all the methods; for the nonsmooth generalized-$\alpha$ method, $\rho_\infty$ is chosen as 0.8, for the Moreau-Jean time stepping method, $\theta = 1$ and for the fully implicit Newmark, $\gamma = 1$ and $\beta = 0.5625$. Figure 3 shows the position and velocity of the ball.

Then, Figure 4 demonstrates the inconsistency of some alternative schemes for this simple example. The smooth generalized-$\alpha$ algorithm ($\rho_\infty = 0.8$) is based on the assumption that both bilateral and unilateral forces are treated as non-impulsive quantities and are integrated according to the generalized-$\alpha$ second-order method. The fully implicit Newmark scheme is a special case of this method with $\gamma = 1$ and $\alpha_f = \alpha_m = 0$. We also consider formulations where the constraint at velocity level in Equation (24) is replaced by the constraint at position level: $\mathbf{0} = \mathbf{\Lambda}_{n+1} - \max(\mathbf{0}, \mathbf{\Lambda}_{n+1} - \mathbf{g}(\mathbf{q}_{n+1}))$. As one can see from Figure 4, the numerical methods based on the position-level constraint are not reliable. The smooth generalized-$\alpha$ method with $\rho_\infty = 0.8$ is not consistent either, even when the constraint is imposed at velocity level.

For the three consistent methods, the errors are computed by comparison with an analytically-exact solution, Appendix A [10]. Figure 5 shows the convergence analysis of these three methods.

<figure><img src="image_4.png" /><figcaption>Figure 3. Position and velocity of the bouncing ball vs. time: numerical methods with velocity-level constraint.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Position and velocity of the bouncing ball vs. time: numerical methods with non-consistent response.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5. Relative errors on L¹ norm for the bouncing ball: (a) position, (b) velocity.</figcaption></figure>

The relative error is analyzed on the L¹ norm, which is defined for constant stepsize as

$$ \|e\|_1 = \sum_{i=0}^{N} |e_i| / \sum_{i=0}^{N} |f(t_i)|, \qquad (42) $$

where N is the number of time steps, $e_i = f_i - f(t_i)$, $f_i$ is the numerical solution and $f(t_i)$ is the exact solution.

The convergence analysis is made on the interval [0, 4]s. As one can see from the figure, all the three methods remain first order accurate in the overall range. However, the nonsmooth generalized-α and the fully implicit Newmark methods have a slightly better accuracy.

## 5.2. Bouncing of a linear oscillator

In this example, the bouncing of a vertical linear oscillator model is studied, Figure 2(b). The oscillator consists of two masses connected by a spring. It is subject to the gravity force and has two DOFs in the vertical direction. After the lower mass impacts against the plane, it bounces back with a restitution coefficient of $e = 0.8$. In the meanwhile, it is also subjected to a force by the compressed spring. Thus, a second impact or multiple impacts can occur right after the first impact. In the free-flight mode after impacts, the system is oscillating with its natural frequency. Physical parameters used in this model are: mass $m = 1$kg and radius $R = 0.2$m for each ball; the stiffness of the spring $k = 10^4$N/m and the unstretched length $l = 1$m; the initial velocity is zero and the initial height $h_0 = 1.001$m; the gravity acceleration $g = 10$m/s².

The nominal time step size is set as $h = 5 \cdot 10^{-4}$s and the other algorithmic parameters are: $\rho_\infty = 0.8$ for the nonsmooth generalized-α method; $\theta = 1$ for the Moreau-Jean method and $\gamma = 1$ and $\beta = 0.5625$ for the fully implicit Newmark method. Results in Figure 6 show respectively the position and velocity of the bottom mass, and the numerical results are compared with a semi-analytical solution, Appendix B. In the analytical analysis, the exact time instants can be obtained using an iterative method. Thus, here in this example, accuracy is compared for the results including only a couple of impacts: from the first impact $t_0 = 0.4$s till the occurrence of the third impact $t_1 = 1.0$s. Figure 7 shows the convergence rate on L¹ norm. Figure 8 shows the total energy of the oscillator system.

As one can see from Figure 6, for all the three methods, double impacts are well accounted for. During the free flight after the first double impacts, the system is oscillating with slight numerical energy dissipation using the nonsmooth generalized-α method; however, the numerical energy dissipation tends to be comparatively larger for the fully implicit Newmark and Moreau-Jean methods. These findings can also be referred to Figure 1.

From the comparison to the exact solution, the nonsmooth generalized-α has better accuracy than the other two methods with the same time step size. In [42], the authors demonstrated that decreasing the step size for Moreau-Jean method, the energy dissipation gets smaller and the results tend to be closer to the nonsmooth HHT method — a particular case of the nonsmooth generalized-α
<figure><img src="image_7.png" /><figcaption>Figure 6. Numerical results of the bottom ball of the bouncing oscillator: (a) position, (b) velocity.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Relative errors on $L^1$ norm for the bouncing oscillator: (a) position, (b) velocity.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Numerical results for the total energy of the bouncing oscillator.</figcaption></figure>

method. Similar analysis will also be shown in the following examples. Though the order of accuracy remains first order for all the three methods, the error of the nonsmooth generalized-$\alpha$ scheme is smaller than the other two, Figure 7. Figure 8 further suggests that nonsmooth generalized-$\alpha$ method has better quality of global energy behaviour compared with the other two methods.

With the same computer hardware configuration, the relative central processing unit time needed to compute the bouncing oscillator system model over the interval $[0,4]$s is compared in Table I. For a given time step size, the computational effort is similar for the three methods. In Figure 7, it can be
<table><caption>Table I. Relative CPU time of the time integration schemes for the bouncing oscillator.</caption><thead><tr><th></th><th>NS generalized-α<br>(h = 5·10<sup>-4</sup>s)</th><th>Moreau–Jean<br>(h = 5·10<sup>-4</sup>s)</th><th>Fully implicit Newmark<br>(h = 5·10<sup>-4</sup>s)</th><th>Moreau–Jean<br>(h = 5·10<sup>-5</sup>s)</th></tr></thead><tbody><tr><td>Rel. CPU time</td><td>1.02</td><td>1</td><td>1.01</td><td>10.42</td></tr></tbody></table>

observed that the nonsmooth generalized-α is almost 10 times more accurate than the Moreau–Jean scheme. As it is difficult to set a specific relative error for each scheme, the Moreau–Jean scheme with one-tenth of the step size is presented for comparison.

### 5.3. Bouncing of an elastic bar

This benchmark illustrates the behaviour of the algorithm when a closed contact occurs. The elastic bar is placed above a rigid ground. It is dropped from an initial height of $h_0$, with a uniform initial velocity of $v_0$ and bounces back against the ground. Before the impact, the bar is undeformed. Once the bar reaches the ground, a shock wave is caused by the impact and travels from the bottom to the top, and then again, travels back to the bottom. During the travelling of the shock wave, the bar stays in contact with the ground because of its compression during a period of $\Delta t = 2L\sqrt{\rho/E}$. As soon as the wave reaches back to the bottom, the bar takes off and the contact is released.

An analytically-exact solution for this benchmark is detailed in [22]. For comparison, the same parameters as in [22] are applied in this test example: Youngs Modulus $E = 900$Pa, density $\rho = 1$kg/m$^3$, undeformed initial length $L = 10$m, initial height from the ground to the bottom of the bar $h_0 = 5$m and initial velocity $v_0 = 10$m/s. In this case, the bar stays in closed contact during a period of $\Delta t = 2/3$s. The gravity acceleration $g$ is set to 0 so that only one close impact will occur. It is worth noting that the restitution coefficient for the impact is chosen as $e = 0$. However, this benchmark does not involve any impact between rigid bodies, and we observe in preliminary tests that the numerical response tends to be insensitive with respect to the value of $e \in [0, 1)$ when the spatial and temporal discretization is refined.

This bouncing bar example is a one-dimensional problem. The bar is uniformly discretized in space by 200 linear finite elements. Time step size can be chosen based on the evaluation of the Courant number, a relevant ratio, which links the mesh size and the step size [22]. The step size with this mesh discretization is then chosen as $h = 2·10^{-3}$s. Other algorithmic parameters are as: $\rho_\infty = 0.8$ for the nonsmooth generalized-α method; $\theta = 1$ for the Moreau–Jean method and $\gamma = 1$ and $\beta = 0.5625$ for the fully implicit Newmark method.

Figure 9 shows the position and the pressure on the bottom of the bar. Also, the total energy of the bouncing elastic bar is analyzed, as shown in Figure 10. The numerical results of the position response and the pressure are compared with the exact solution. As one can tell from the

<figure><img src="image_10.png" /><figcaption>Figure 9. Comparison of the numerical results for the bouncing elastic bar: (a) position, (b) pressure.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10. Numerical results for the total energy of the bouncing elastic bar.</figcaption></figure>

figures, closed contact analysis is stable for all the three methods. Compared with Moreau–Jean and fully implicit Newmark methods, the nonsmooth generalized-α method has better accuracy for the position response and the pressure, in particular for the period near/after the take-off. An energy performance analysis is carried out and finer step sizes ($h = 2 \cdot 10^{-4}$s and $h = 2 \cdot 10^{-5}$s) are adjusted for the Moreau–Jean method in order to compare the energy dissipation. The analysis of the Moreau–Jean scheme shows that the numerical solution tends to converge to a behaviour without energy dissipation after impact. As can be seen from Figure 10, the nonsmooth generalized-α method has the least energy dissipation after the impact. One can thus conclude that the nonsmooth generalized-α scheme has better quality of global energy behaviour, as already shown in Section 5.2.

### 5.4. Impact of a flexible rotating beam

The following example studies the impact of a flexible rotating beam. The beam works as a flexible simple pendulum, that is, it is subjected to the gravity force and swings around a pivot point, as shown in Figure 2(d). The beam is modelled using a geometrically exact, two-dimensional finite element formulation, based on Timoshenko’s theory. The planar beam finite element is capable of handling arbitrarily large finite rotations. It is a special case of the 3D model presented in [4]. It is of interest to verify the nonsmooth time integration methods for such a highly nonlinear application.

The beam studied in this example has a square cross section. The parameters of the beam are as follows: length $l = 1$m, width $b = 0.01$m, height $d = 0.01$m, cross-sectional area $A = bd$, cross-sectional inertia $I = bd^3/12$, reduced cross-sectional area for shear $A_2 = 5/6A$, Young’s modulus $E = 2.1 \cdot 10^{11}$N/m$^2$, Poisson ratio $\nu = 0.3$ and density $\rho = 7800$kg/m$^3$. The acceleration of gravity is $g = 10$m/s$^2$.

In the beginning, the beam is placed horizontally with the tip pointing to the right. It is released from rest and is restricted to planar motion under gravity. An obstacle is placed such that a unilateral constraint is applied as follows: $\sqrt{2}/2 \le x \perp \lambda^u \ge 0$. Likewise, the complementarity condition is expressed at the velocity level as in Equation (4). The coefficient of restitution is set as $e = 0.8$.

The beam is discretized into four finite elements. The model is first simulated with time step size $h = 5 \cdot 10^{-4}$s. Because it is difficult to obtain the exact solution, a finer step size of $h = 2 \cdot 10^{-5}$s is then chosen for the Moreau–Jean method in order to compare the accuracy. This makes sense as the Moreau–Jean method is commonly known as a validated method and has been proven to be convergent. Other numerical parameters are as: $\rho_\infty = 0.8$ for the nonsmooth generalized-α method, $\theta = 1$ for the Moreau–Jean method and $\gamma = 1$ and $\beta = 0.5625$ for the fully implicit Newmark method.

Figure 11 shows the response of position and velocity at the tip of the beam in x-direction. Comparatively, the nonsmooth generalized-α scheme tends to have results much closer to the Moreau–Jean method with a smaller time step size. Higher numerical damping in the Moreau–Jean and fully implicit Newmark methods leads to higher energy dissipation even during the smooth
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11. Numerical results for the impact of a flexible rotating beam: (a) position, (b) velocity.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12. Numerical results for the total energy of the bouncing beam.</figcaption></figure>

motion, as can be seen in the velocity response in Figure 11(b). However, higher numerical damping does not necessarily result in better stability performance. As a matter of fact, if one further increase the time step size to $10^{-3}$s, the nonsmooth generalized-$\alpha$ scheme is still convergent in Newton iteration, but the other two methods are not. If one chooses $\theta = 1/2$ for the Moreau-Jean method, which is equivalent to the trapezoidal method in the free flight phases (i.e. no numerical damping is present), the results become divergent. It indicates that controllable numerical damping is necessary in this highly nonlinear case. Comparison in the overall scale can be seen in the energy performance, as shown in Figure 12. The nonsmooth generalized-$\alpha$ method has the best energy performance between impacts, where no physical dissipation is expected.

### 5.5. Impact of a rigid pendulum

This example studies the impact of a simple rigid pendulum with both unilateral and bilateral constraints. The purpose is to validate the proposed time integration methods in the regime of nonlinear, bilaterally-constrained problems. 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 Cartesian coordinates by the DOFs: $\mathbf{q} = [x \ y \ \theta]^T$. The nonlinear bilateral constraints are expressed as: $\mathbf{\Phi} = [x - l \cos \theta \ y + l \sin \theta]^T$. The pendulum is subjected to the gravity acceleration $g$, and thus the vector of external forces reads: $\mathbf{f}^{\text{ext}} = [0 \ -mg \ 0]^T$. The mass matrix is $\mathbf{M} = \text{diag}[m \ m \ J]$.

The physical parameters of the model are as: mass $m = 1$kg, moment of inertia $J = 0.1$kg·m$^2$, length $l = 1$m, gravity acceleration $g = 10$m/s$^2$ and initial position $\mathbf{q} = [1 \ 0 \ 0]^T$. The pendulum is released from rest and swings clockwise because of the gravity force. An obstacle is placed
such that a unilateral constraint is enforced as $\sqrt{2}/2 \le x \perp \lambda^u \ge 0$, and is then expressed at velocity level. The coefficient of restitution is chosen as $e = 0.8$.

In the numerical tests, the nominal time step size is $10^{-3}$s; $\rho_\infty = 0.8$ for the nonsmooth generalized-$\alpha$ method, $\theta = 1$ for the Moreau-Jean time stepping method and $\gamma = 1$ and $\beta = 0.5625$ for the fully implicit Newmark method. The model is first simulated without coordinate projection to deal with the drifting phenomena, and then, projection method is applied with the Moreau-Jean method at every time step for comparison. The $\|\cdot\|_2$ norm is used for the projection method, as

$$ \min_{q_n} \|q_n - \bar{q}_n\|_2 \quad \text{s.t.} \quad \Phi(q_n) = 0, \qquad (43) $$

where $\bar{q}_n$ and $q_n$ denote the vector of generalized coordinates before and after projection, respectively.

Figure 13 shows the results of position and velocity at the tip of the pendulum in x-direction, and also the drifting of the bilateral constraints. All the three methods have consistent results. However, Moreau-Jean method without projection has larger drift-off effect than the other two methods, for which the drifting of the bilateral constraints is less than $5 \cdot 10^{-6}$m. Figure 14 shows the results of the bilateral and unilateral constraints, which are expressed in terms of impulse. Figure 15 shows the convergence rate on $L^1$ norm. The time interval for the relative error analysis is chosen from the beginning till $t = 1$s. Because the analytical solution for the bouncing pendulum example is difficult to obtain, it is reasonable to use as a reference solution the numerical results obtained using the well-validated Moreau-Jean method with a tiny stepsize $h = 10^{-7}$s. The order of accuracy remains first order with larger errors with the Moreau-Jean method compared with the other two.

<figure><img src="image_15.png" /><figcaption>Figure 13. Comparison of the numerical results for the impact of a simple rigid pendulum: (a) position, (b) pressure, (c) drifting of the bilateral constraints.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14. Lagrangian multipliers of the nonsmooth generalized-$\alpha$ method: (a) bilateral, (b) unilateral.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 15. Relative errors on $L^1$ norm for the bouncing pendulum: (a) position, (b) velocity.</figcaption></figure>

## 5.6. Multiple impacts of a crankasiDER mechanism

This section deals with the analysis of multiple impacts of a planar crankasiDER mechanism with a
translational clearance joint. This mechanism consists of four bodies: the crank, the connecting rod,
the slider and the ground. These bodies are connected via three revolution joints and one translation
joint. The revolution joints are modelled as ideal constraints, whereas clearance is considered in
the translation joint, which implies four unilateral constraints on the four vertices and consequently
results in the impacts of the slider with the guide surfaces during the rotation of the crank. The
slider and the ground are modelled as rigid bodies. For the sake of comparison, the crank and the
connecting rod are firstly modelled as rigid bodies and then are modelled with flexibility using beam
elements as described in Section 5.4.

**5.6.1. Rigid crankasiDER mechanism.** In this study, the four bodies of the crank-slider mechanism are modelled as rigid bodies. Similar physical parameters are taken from [17,43], where only rigid multibody systems are taken into account. The geometrical parameters are: length of the crank $l_1 = 0.153\text{m}$, length of the connecting rod $l_2 = 0.306\text{m}$, width of the slider $2a = 0.05\text{m}$, length of the slider $2b = 0.1\text{m}$ and distance between the guide surfaces $d = 0.052\text{m}$. Initial configuration parameters are $\theta_1 = 0\text{rad}, \theta_2 = 0\text{rad}, \theta_3 = 0\text{rad}$ and both the upper and lower clearances between the slider and the guide surfaces are $c = 0.001\text{m}$.

The mass and inertia of the crank: $m_1 = 0.038\text{kg}$ and $J_1 = 7.4 \cdot 10^{-5}\text{kg m}^2$; the connecting rod: $m_2 = 0.038\text{kg}$ and $J_2 = 5.9 \cdot 10^{-4}\text{kg m}^2$; the slider: $m_3 = 0.076\text{kg}$ and $J_3 = 2.7 \cdot 10^{-6}\text{kg m}^2$.

The crank is released with an initial speed $\omega_1 = 150\text{rad/s}$ and the initial speed of the connecting rod is set as $\omega_2 = -75\text{rad/s}$. The gravity acceleration $g = 10\text{m}^2/\text{s}$ and the restitution coefficient $e = 0.4$.

The model is first simulated with time step size $h = 10^{-4}\text{s}$. In order to study the influence of the step size, the Moreau-Jean method with a finer step size of $h = 10^{-5}\text{s}$ is also considered. Other numerical parameters are as: $\rho_\infty = 0.8$ for the nonsmooth generalized-$\alpha$ method, $\theta = 1$ for the Moreau-Jean method, and $\gamma = 1$ and $\beta = 0.5625$ for the fully implicit Newmark method.

Figure 16 shows the position of the corners of the slider over two crank revolutions. Multiple impacts and closed contacts occur in all the four vertices with the guide surfaces, and are correctly represented by the three algorithms. Figure 17 shows the crank speed and the position of the centre of the crank. One can observe that the simulation results show a good match with [17,43].

**5.6.2. Flexible crankasiDER mechanism.** From an optimization point of view, structural flexibility is usually taken into account in the real application of a crank-slider mechanism. In this study, the configuration of the mechanism is the same as the previous example except that the crank and the connecting rod are replaced by flexible beams instead of rigid bodies. The beam applied in this example has the same material parameters as described in section 5.4. Both the crank and the connecting rod have a square cross section. In order to apply the same mass for the crank and the
<table><tr><td>Nonsmooth generalized-α (h=10<sup>-4</sup>s)</td><td>Moreau-Jean (h=10<sup>-4</sup>s)</td><td>Fully implicit Newmark (h=10<sup>-4</sup>s)</td><td>Moreau-Jean (h=10<sup>-5</sup>s)</td></tr></table>

<figure><img src="image_18.png" /><figcaption>Figure 16. Slider corners' position of the rigid crank-slider mechanism with different time integration schemes.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 17. Simulation results with the nonsmooth generalized-α method: (a) crank speed, (b) position of the center of the slider.</figcaption></figure>

connecting rod as in the rigid crank-slider mechanism, the width and the height of the crank are chosen as $b_1 = 0.0064m$ and $d_1 = 0.005m$; the width and the height of the connecting rod are $b_2 = 0.0032m$ and $d_2 = 0.005m$.

The mechanism is released from rest in the beginning and is constrained in a planar motion under the gravity force. Initial rotating speed of the crank and the connecting rod is enforced as $\omega_1 = 150\text{rad/s}$ and $\omega_2 = -75\text{rad/s}$ respectively and a full initial velocity vector of each node is calculated in a way assuming that it is a rigid mechanism. After the mechanism is released, impacts occur between the four vertices of the crank and the guide surfaces, and therefore, cause vibration on the flexible parts. The purpose of this analysis is to compare the numerical dissipation of different time integration methods with a preference that the low-frequency vibration will be maintained so that accuracy may be satisfactory.

Both the flexible crank and connecting rod are discretized into two finite elements. The model is first simulated with time step size $h = 10^{-4}\text{s}$, then a finer step size of $h = 10^{-6}\text{s}$ is chosen for the Moreau-Jean method for comparison. Other numerical parameters are as: $\rho_\infty = 0.8$ for the nonsmooth generalized-α method, $\theta = 1$ for the Moreau-Jean method and $\gamma = 1$ and $\beta = 0.5625$ for the fully implicit Newmark method.
<figure><img src="image_20.png" /><figcaption>Figure 18. Slider corners’ position of the flexible crank-slider mechanism with different time integration schemes : (a) slider corners’ position vs. crank revolutions and (b) slider corners’ position versus time.</figcaption></figure>

Figure 18 shows the corner positions of the slider. Compared with the other two schemes with the same step size, the results obtained from the nonsmooth generalized-α method match much better with those from the Moreau–Jean scheme with a 100 times finer step size, especially for the time-domain response, Figure 18(b). The main cause is that the nonsmooth generalized-α scheme dissipates less of the vibration of the flexible crank and connecting rod after impacts, and consequently attains better accuracy. Similar observation can be found in Figure 19 in terms of the crank speed and energy behaviour. It further proves that the nonsmooth generalized-α method saves much computation cost as the scale of step size is critical in the simulation time, in particular when it comes to a complex flexible multibody system. The relative central processing unit time for the computation of the flexible crank-slider mechanism is compared in Table II. For a given time step size, the computational effort is similar for the three methods.
<figure><img src="image_21.png" /><figcaption>Figure 19. Comparison of different time integration methods for the flexible crank-prüfer mechanism : (a) crank speed and (b) total energy.</figcaption></figure>

<table><caption>Table II. Relative central processing unit time of the time integration schemes for the flexible crank-prüfer mechanism.</caption><thead><tr><th></th><th>NS generalized-α<br>(h = 10<sup>-4</sup>s)</th><th>Moreau-Jean<br>(h = 10<sup>-4</sup>s)</th><th>Fully implicit Newmark<br>(h = 10<sup>-4</sup>s)</th><th>Moreau-Jean<br>(h = 10<sup>-6</sup>s)</th></tr></thead><tbody><tr><td>Rel. CPU time</td><td>1.07</td><td>1</td><td>1.05</td><td>88.24</td></tr></tbody></table>

## 6. CONCLUSIONS

Considering that many mechanical systems are subjected to both bilateral constraints and unilateral constraints, an extension of the generalized-α method is proposed to account for the nonsmooth phenomena induced by unilateral constraints. The method presented here is denoted as the nonsmooth generalized-α method, and is expected to benefit from good qualitative properties of the classical generalized-α and HHT methods in terms of accuracy, stability and controllable numerical damping.

Numerical tests are studied and results show that the Moreau-Jean, the fully implicit Newmark and the nonsmooth generalized-α methods are valid when unilateral constraints are expressed at the velocity level. Unlike event-capturing schemes, these methods can be used in the presence of accumulation phenomena. Based on the analysis of test examples including linear and nonlinear problems, rigid and flexible system models, unilateral and bilateral constraints, the three methods are stable and the presented nonsmooth generalized-α method has considerably better quality of accuracy, especially for the description of vibrations in flexible system dynamics. The nonsmooth generalized-α method allows the introduction of high-frequency numerical dissipation, which is important for the solution of stiff numerical models while preserving a good energy behaviour in the low-frequency range. Moreover, in the presence of holonomic bilateral constraints, one finds out
that the nonsmooth generalized-α method has less drifting phenomena compared with the Moreau-Jean method. Therefore, large time steps can be used to obtain a physically acceptable numerical solution. Hence, we conclude that the nonsmooth generalized-α time integration method can substantially decrease the computational cost for high-fidelity models of flexible systems, such as wind turbine systems, in which unilateral contacts might occur.

Because of the first order formulation of constraint forces as impulses, the accuracy is brought back to first order in case of contact or in the presence of constraints. Future efforts can be put in obtaining improved accuracy of the nonsmooth generalized-α method, for example, to maintain second-order accuracy for bilaterally constrained systems in case of no contact. The extension of the method to constrained systems with singular mass matrix may also be investigated.

## APPENDIX A: ANALYTICAL SOLUTION OF BOUNCING BALL

Reference [10] presents an analytical solution of the bouncing ball example, which is a particular case with the parameters $g = 2m/s^2$ and $e = 0.5$. A general analytical solution of the bouncing ball example is however derived in this section. With the aforementioned notations, the solution can be written as:

$$ 
\left\{ 
\begin{aligned}
q(t) &= -\frac{1}{2}gt^2 + h_0 \\
v(t) &= -gt
\end{aligned}
\right. 
\qquad t \in [0, \sqrt{\frac{2(h_0-R)}{g}}), \\
\left\{ 
\begin{aligned}
q(t) &= -\frac{1}{2}g(t-t_n)^2 - e^n v_1^{-}(t-t_n) + R \\
v(t) &= -g(t-t_n) - e^n v_1^{-}
\end{aligned}
\right. 
\qquad t \in [t_n, t_{n+1}), \\
\left\{ 
\begin{aligned}
q(t) &= R \\
v(t) &= 0
\end{aligned}
\right. 
\qquad t \in [\frac{1+e}{1-e}\sqrt{\frac{2(h_0-R)}{g}}, +\infty),
$$

with

$$ 
\begin{aligned}
t_n &= \sqrt{\frac{2(h_0 - R)}{g}} \left( 2e^{\frac{1-e^{n-1}}{1-e}} + 1 \right), \\
v_1^{-} &= -\sqrt{2g(h_0 - R)},
\end{aligned}
$$

where $n$ denotes the $n$-th impact.

## APPENDIX B: EXACT SOLUTION OF BOUNCING LINEAR OSCILLATOR

The solution of the oscillator example can be derived in the following way. Denote

$$ Q = ( q_1 \quad \dot{q}_1 \quad q_2 \quad \dot{q}_2 )^T, $$

then the equations of motion can be written in the form of state equations as:

$$ M\dot{Q} + KQ = F, \qquad (B.1) $$

where

$$ M = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & m & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & m \end{pmatrix}, \quad K = \begin{pmatrix} 0 & -1 & 0 & 0 \\ k & 0 & -k & 0 \\ 0 & 0 & 0 & -1 \\ -k & 0 & k & 0 \end{pmatrix}, \quad F = \begin{pmatrix} 0 \\ -mg - kl \\ 0 \\ -mg + kl \end{pmatrix}. $$

The analytical solution can then be derived by the following equation:

$$ Q = e^{A(t-t_0)} Q_0 + \int_{t_0}^{t} e^{A(t-s)} F(s) ds, \qquad (B.2) $$
where $\mathbf{A} = -\mathbf{M}^{-1}\mathbf{K}$ and $\mathbf{Q}_0$ denotes the initial condition at the instant $t_0$. According to the Cayley-Hamilton method, the matrix exponential $e^{\mathbf{A}t}$ can be expressed in a standard polynomial form as

$$e^{\mathbf{A}t} = \sum_{i=0}^{n-1} \alpha_i(t) \mathbf{A}^i, \qquad (\text{B.3})$$

where $n$ is the dimension of $\mathbf{A}$: $n = 4$ in this example. The coefficients $\alpha_i(t)$ are as follows:

$$\begin{align*}
\alpha_0(t) &= 1, \\
\alpha_1(t) &= t, \\
\alpha_2(t) &= m/(2k) \left[ 1 - \cos \left( t \sqrt{2k/m} \right) \right], \tag{B.4} \\
\alpha_3(t) &= \left( \sqrt{m/(2k)} \right)^3 \left[ t \sqrt{2k/m} - \sin \left( t \sqrt{2k/m} \right) \right].
\end{align*}$$

As both $\mathbf{A}$ and $\mathbf{F}$ are constant matrices, integration of the second term in Equation (B.2) is actually the integral of $\alpha_i(t)$, which are expressions in simple polynomials and trigonometric functions. Therefore, Equation (B.2) can be analytically expressed of simple elementary functions. However, the exact solution of $\mathbf{Q}$ computed from Equation (B.2) works only for the smooth situation. After each discontinuity due to the impacts, $\mathbf{Q}_0$ and $t_0$ should be updated. Here, the impact time instants are numerically obtained by solving $q_1 = R$ using the Newton-Raphson recursive algorithm. With the Newton impact law, the initial condition after impact can then be updated.

This method allows to compute the solution when the number of impacts is finite, but it cannot be used to describe an accumulation phenomenon. In the example presented in Section 5.2, with step size as $h = 10^{-4}s$, the exact solution of the system is computed till $t = 3.46s$, just before the accumulation and the closed contact situation.

## ACKNOWLEDGEMENTS

This research work was carried out under grant number 850533 (DYNAWIND) from the Walloon Region (Belgium) which is gratefully acknowledged. The second author has been supported by the French National Research Agency (ANR) through COSINUS program (project SALADYN ANR-08-COSI-014). The third author would like to acknowledge the Belgian National Fund for Scientific research (FRIA) for its financial support.

## REFERENCES

1. Hilber H, Hughes T, Taylor R. Improved numerical dissipation for time integration algorithms in structural dynamics. *Earthquake Engineering and Structural Dynamics* 1977; 5(3):283-292.

2. Chung J, Hulbert G. A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-$\alpha$ method. *ASME Journal of Applied Mechanics* 1993; 60(2):371-375.

3. Newmark NM. A method for computation for structural dynamics. *ASCE Journal of the Engineering Mechanics Division* 1959; 85(3):67-94.

4. Géradin M, Cardona A. *Flexible Multibody Dynamics: A Finite Element Approach*. John Wiley & Sons: New York, 2001.

5. Cardona A, Géradin M. Time integration of the equations of motion in mechanism analysis. *Computers and Structures* 1989; 33(3):801-820.

6. Arnold M, Brüls O. Convergence of the generalized-$\alpha$ scheme for constrained mechanical systems. *Multibody System Dynamics* 2007; 18(2):185-202.

7. Lankarani H, Nikravesh P. Continuous contact force models for impact analysis in multibody analysis. *Nonlinear Dynamics* 1994; 5(2):193-207.

8. Pfeiffer F. On non-smooth dynamics. *Meccanica* 2008; 43(5):533-554.

9. Acary V, Brogliato B. *Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics*, Lecture Notes in Applied and Computational Mechanics, Vol. 35. Springer Verlag: Berlin Heidelberg, 2008.

10. Acary V. Higher order event capturing time-stepping schemes for nonsmooth multibody systems with unilateral constraints and impacts. *Applied Numerical Mathematics* 2012; 62(10):1259-1275.
11. Acary V. Numerical methods for the time integration of nonsmooth mechanical systems. Part II. Nonsmooth mechanical systems with impacts. published in 6th Workshop SDS2010 Structural Dynamical System: Computational Aspects, Monopoli, Italy, June 8-11, 2010.

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

13. Paoli L, Schatzman M. A numerical scheme for impact problems II: The multi-dimensional case. *SIAM Journal of Numerical Analysis* 2002; **40**(2):734–768.

14. Moreau JJ. Unilateral contact and dry friction in finite freedom dynamics. In *J.J. Moreau and P.D. Panagiotopoulos, Editors, Nonsmooth Mechanics and Applications, CISM Courses and Lectures, Wien - New York: Springer Verlag* 1988; **302**:1–82.

15. Jean M, Moreau JJ. Unilaterality and dry friction in the dynamics of rigid bodies collections. In *Proc. of Contact Mech. Int. Symp.*, Vol. 1, Curnier A (ed.). Presses Polytechniques et Universitaires Romandes: Lausanne, Switzerland, 1992; 31–48.

16. Jean M. The non-smooth contact dynamics method. *Computer Methods in Applied Mechanics and Engineering* 1999; **177**(3-4):235–257.

17. Acary V. Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb's friction. *Computer Methods in Applied Mechanics and Engineering* 2013; **256**(10):224–250.

18. Studer C. *Numerics of Unilateral Contacts and Friction: Modeling and Numerical Time Integration in Non-smooth Dynamics*, Lecture Notes in Applied and Computational Mechanics, Vol. 47. Springer Verlag: Berlin Heidelberg, 2009.

19. Studer C, Leine RI, Glocker C. Step size adjustment and extrapolation for time stepping schemes in non-smooth dynamics. *International Journal for Numerical Methods in Engineering* 2008; **76**(11):1747–1781.

20. Schindler T, Acary V. Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: definition and outlook. *Mathematics and Computers in Simulation* 2013. In press.

21. Krause R, Walloth M. Presentation and comparison of selected algorithms for dynamic contact based on the Newmark scheme. *Applied Numerical Mathematics* 2012; **62**(10):1393–1410.

22. Doyen D, Ern A, Piperno S. Time-integration schemes for the finite element dynamic Signorini problem. *SIAM Journal on Scientific Computing* 2011; **33**(1):223–249.

23. Simo JC, Tarnow N. The discrete energy-momentum method. Conservering algorithms for nonlinear elastodynamics. *Zangnew Math Phys* 1992; **43**(5):757–793.

24. Laursen T, Chawla V. Design of energy conserving algorithms for frictionless dynamic contact problems. *International Journal for Numerical Methods in Engineering* 1997; **40**(5):863–886.

25. Chawla V, Laursen T. Energy consistent algorithms for frictional contact problem. *International Journal for Numerical Methods in Engineering* 1998; **42**(5):799–827.

26. Laursen T, Love G. Improved implicit integrators for transient impact problems-geometric admissibility within the conserving framework. *International Journal for Numerical Methods in Engineering* 2002; **53**(2): 245–274.

27. Laursen T, Love G. Improved implicit integrators for transient impact problems: Dynamic frictional dissipation within an admissible conserving framework. *Computer Methods in Applied Mechanics and Engineering* 2003; **192**(19):2223–2248.

28. Hauret P, Le Tallec P. Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact. *Computer Methods in Applied Mechanics and Engineering* 2006; **195**(37-40):4890–4916.

29. Zolghadr Jahromi H, Izzuddin B. Energy conserving algorithms for dynamic contact analysis using Newmark methods. *Computers and Structures* 2013; **118**:74–89.

30. Vola D, Pratt E, Jean M, Raous M. Consistent time discretization for dynamical frictional contact problems and complementarity techniques. *European Journal of Computational Mechanics / Revue européenne des éléments finis* 1998; **7**(1-2-3):149–162.

31. Dumont Y, Paoli L. Vibrations of a beam between obstacles. Convergence of a fully discretized approximation. *ESAIM, Mathematical Models and Numerical Analysis* 2006; **40**(4):705–734.

32. Schatzman M, Bercovier M. Numerical approximation of a wave equation with unilateral constraints. *Mathematics of Computations* 1989; **53**(187):55–79.

33. Carpenter NJ, Taylor RL, Katona MG. Lagrange constraints for transient finite element surface contact. *International Journal for Numerical Methods in Engineering* 1991; **32**(1):103–128.

34. Khenous H, Laborde P, Renard Y. On the discretization of contact problems in elastodynamics. In *Analysis and simulation of contact problems*, Vol. 27, Wriggers P, et al. (ed.), Lecture Notes in Applied and Computational Mechanics. Springer Verlag: Berlin Heidelberg, 2006; 31–38.

35. Khenous HB, Laborde P, Renard Y. Mass redistribution method for finite element contact problems in elastodynamics. *European Journal of Mechanics, A, Solids* 2008; **27**(5):918–932.

36. Deufhard P, Krause R, Ertel S. A contact-stabilized Newmark method for dynamical contact problems. *International Journal for Numerical Methods in Engineering* 2007; **73**(9):1274–1290.

37. Eich E. Convergence results for a coordinate projection method applied to mechanical systems with algebraic constraints. *SIAM Journal on Numerical Analysis* 1993; **30**(5):1467–1482.

38. Alart P, Curnier A. A mixed formulation for frictional contact problems prone to Newton like solution method. *Computer Methods in Applied Mechanics and Engineering* 1991; **92**(3):353–375.
39. Christensen P, Klarbring A, Pang J, Stromberg N. Formulation and comparison of algorithms for frictional contact problems. *International Journal for Numerical Methods in Engineering* 1998; **42**(1):145–172.

40. Cottle R, Pang J, Stone R. *The Linear Complementarity Problem*. Academic Press: Boston, MA, 1992.

41. Klarbring A. A mathematical programming approach to three-dimensional contact problems with friction. *Computer Methods in Applied Mechanics and Engineering* 1986; **58**(2):175–200.

42. Chen Q, Acary V, Virlez G, Brüls O. A Newmark-type integrator for flexible systems considering nonsmooth uni-lateral constraints. *Proceedings of the 2nd Joint International Conference on Multibody System Dynamics*, Stuttgart, Germany, 2012.

43. Flores P, Leine R, Glocker C. Modeling and analysis of planar rigid multibody systems with translational clearance joints based on the non-smooth dynamics approach. *Multibody System Dynamics* 2010; **23**(2):165–190.