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

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. Copyright © 2013 John Wiley & Sons, Ltd.


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 488 Q.-Z. CHEN ET AL. by vibrations of structural elements such as the blades and the tower coupled with strong nonlinearities 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 continuous 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 experiments, 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, usually 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 spacediscretized 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 timestepping 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][10][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][15][16] is based on a Â-method [17]. Time-stepping methods have been proven to be convergent and robust, and are extensively applied as the solution to the 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 [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 integration schemes for structural and flexible multibody systems, that is, HHT method and generalized-m ethod, 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-˛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 including 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 instabilities 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 modifications 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][24][25][26][27][28][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][15][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 489 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 implemention 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.

EQUATIONS OF MOTION
The equations of a flexible multibody system including bilateral and unilateral constraints can be expressed in the following form: .q/ D 0,

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 with where P e v is an acceleration variable related with the smooth force f, and dw is a differential measure because of the contact reaction and bilateral force. A key idea of the proposed method is to integrate the variable P e v using a second-order accurate generalized-˛method, whereas the differential measure dw 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 where M nC1 D M.q nC1 / and f nC1 D f.q nC1 , v nC1 , t nC1 /. For the other part, we keep the values of the impulse over the time step .t n , t nC1 as primary unknowns, that is, where W nC1 represents a velocity jump over the time step; ƒ b nC1 and ƒ u nC1 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, Z The discretization of the contact law is written as if N g 6 0 then 0 6 g q,nC1 v nC1 C e g q,n v n ? ƒ u nC1 > 0, (19) where N g is defined either as N g D g n or as N g D g p nC1 , with g p nC1 defined as the prediction of the unilateral constraints at time t nC1 . Equation (19) can be rewritten in an equation form as if N g 6 0 then ‰ ƒ u nC1 , g q,nC1 v nC1 C e g q,n v n D 0, with ‰.a, b/ D a max.0, a Ä b/, Ä > 0. Defining the set of active constraints g A as the set of all the unilateral constraints satisfying N g 6 0 and ƒ A as the corresponding Lagrange multipliers, the discretized form of the equations of motion becomes 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-m ethod is obtained by further taking the impulsive terms into account so that it takes the following form: where e q, e v are intermediate variables representing some smooth positions and velocities respectively, whereas q, v stand for the generalized position and velocity in the general sense; a is an intermediate acceleration-like variable defined as in (27) and initialized at the initial time as The form of Equations (28-29) is justified as follows. Suppose that an impact is detected at time t i 2 OEt n , t nC1 /, 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: where W nC1 also satisfies the discrete equations of motion in (21)(22)(23)(24).
In the unconstrained case for the smooth generalized-˛algorithm, the condition of second-order accuracy yields: A small amount of numerical damping, represented by the spectral radius of the algorithm at the high-frequency limit 1 with 1 2 OE0, 1 can be introduced. The optimal parameters for the generalized-˛method are then derived as [2]: For the HHT algorithm,˛m D 0 and˛f 2 OE0, 1=3, whereas for the Newmark algorithm, f D˛m D 0.
The proposed method will be compared with the Moreau-Jean Â-method, which is obtained by replacing Equations (25)(26)(27)(28)(29) by the following equations with Â 2 OE0, 1. The Moreau-Jean method with Â D 1=2 is equivalent to the nonsmooth generalized-˛scheme with˛m D˛f D 0,ˇD 1=4, D 1=2 and a nC1 D P e v nC1 . 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 !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 R q C ! 2 q D 0. The algorithmic parameters are chosen as follows: for the generalized-˛method, the spectral radius at infinity is chosen as 1 D 0.8; for the Moreau-Jean method, Â D 1, which implies the implicit Euler method; for the so-called fully implicit Newmark method, D 1 andˇD 1=4. C 1=2/ 2 D 0.5625. As one can see from Figure 1, the generalized-˛method has the least numerical dissipation. In contrast, the Moreau-Jean method is the most numerically dissipative one.

A generalized Newton method
Formulae (25)(26)(27)(28)(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 e v, v, ƒ b and ƒ A as independent sets of unknowns. The linearized equation for Newton iteration becomes: 2 where the numerical parameters: t D .1 ˛m/=OE.1 ˛f / h andˇt Dˇh= 1=2 h; the tangent stiffness and damping matrices are respectively defined as: K t D @.M P e v h/=@q and C t D @h=@v; ˆq 2 is the derivative ofˆq ,nC1 v nC1 with respect to q nC1 ; ‰ q , ‰ v and ‰ ƒ A are the partial derivatives of the function ‰ ƒ A nC1 , g A q,nC1 v nC1 C e g A q,n v n Á with respect to q nC1 , v nC1 and ƒ A nC1 respectively and are indeed the generalized Jacobian matrices at the point of discontinuity; res q , res con , resˆand res ‰ denote respectively the residuals of the four sets of equations in (21)(22)(23)(24). Using a superscript '*' to denote the value of a variable at Newton iteration k, the residuals are The detailed computation of generalized Jacobians of ‰ 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.

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 M and a vector r such that 0 6 Mz C r ? z > 0. (21)(22)(23) for iteration k C 1 can be recast as The unilateral constraint at the velocity level is also linearized as Therefore, by substituting Equation (38) into (40), the complementarity condition in Equation (24) is obtained in the standard form as To sum up, the flowchart of the predictor-corrector iterative method for the nonsmooth generalizedmethod based on a sequential LCP is shown in Algorithm 1. 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.

Bouncing ball
The first example is a standard bouncing ball on a rigid plane, as depicted in Figure 2 as follows: mass m D 1kg, radius R D 0.2m, gravity acceleration g D 10m=s 2 , initial height h 0 D 1.001m (with an initial height of h 0 D 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-˛, the Moreau-Jean and the fully implicit Newmark schemes. The numerical parameters are set as: nominal time step h D 10 3 s for all the methods; for the nonsmooth generalized-˛method, 1 is chosen as 0.8, for the Moreau-Jean time stepping method, Â D 1 and for the fully implicit Newmark, D 1 andˇD 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-˛algorithm ( 1 D 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-˛second-order method. The fully implicit Newmark scheme is a special case of this method with D 1 and˛f D˛m D 0. We also consider formulations where the constraint at velocity level in Equation (24) is replaced by the constraint at position level: 0 D ƒ nC1 max.0, ƒ nC1 g.q nC1 //. As one can see from Figure 4, the numerical methods based on the position-level constraint are not reliable. The smooth generalized-˛method with 1 D 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 analyticallyexact solution, Appendix A [10]. Figure 5 shows the convergence analysis of these three methods.  The relative error is analyzed on the L 1 norm, which is defined for constant stepsize as where N is the number of time steps, e i D 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 OE0, 4s. 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.

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 D 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 D 1kg and radius R D 0.2m for each ball; the stiffness of the spring k D 10 4 N=m and the unstretched length l D 1m; the initial velocity is zero and the initial height h 0 D 1.001m; the gravity acceleration g D 10m=s 2 .
The nominal time step size is set as h D 5 10 4 s and the other algorithmic parameters are: 1 D 0.8 for the nonsmooth generalized-˛method; Â D 1 for the Moreau-Jean method and D 1 andˇD 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 semianalytical 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 D 0.4s till the occurrence of the third impact t 1 D 1.0s. Figure 7 shows the convergence rate on L 1 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-C      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-˛scheme is smaller than the other two, Figure 7. Figure 8 further suggests that nonsmooth generalized-m ethod 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  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.

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 t D 2L p =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 D 900Pa, density D 1kg=m 3 , undeformed initial length L D 10m, initial height from the ground to the bottom of the bar h 0 D 5m and initial velocity v 0 D 10m=s. In this case, the bar stays in closed contact during a period of t D 2=3s. 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 D 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 2 OE0, 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 D 2 10 3 s. Other algorithmic parameters are as: 1 D 0.8 for the nonsmooth generalized-˛method; Â D 1 for the Moreau-Jean method and D 1 andˇD 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  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 D 2 10 4 s and h D 2 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-m ethod 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.

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 D 1m, width b D 0.01m, height d D 0.01m, cross-sectional area A D bd , crosssectional inertia I D bd 3 =12, reduced cross-sectional area for shear A 2 D 5=6A, Young's modulus E D 2.1 10 11 N=m 2 , Poisson ratio D 0.3 and density D 7800kg=m 3 . The acceleration of gravity is g D 10m=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: p 2=2 6 x ? u > 0. Likewise, the complementarity condition is expressed at the velocity level as in Equation (4). The coefficient of restitution is set as e D 0.8.
The beam is discretized into four finite elements. The model is first simulated with time step size h D 5 10 4 s. Because it is difficult to obtain the exact solution, a finer step size of h D 2 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: 1 D 0.8 for the nonsmooth generalized-m ethod, Â D 1 for the Moreau-Jean method and D 1 andˇD 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 502 Q.-Z. CHEN ET AL.    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-˛scheme is still convergent in Newton iteration, but the other two methods are not. If one chooses Â D 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-˛method has the best energy performance between impacts, where no physical dissipation is expected.

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: q D OEx y Â T . The nonlinear bilateral constraints are expressed as:ˆD OEx l cos Â y C l sin Â T . The pendulum is subjected to the gravity acceleration g, and thus the vector of external forces reads: f ext D OE0 mg 0 T . The mass matrix is M D diag OEm m J. The physical parameters of the model are as: mass m D 1kg, moment of inertia J D 0.1kg m 2 , length l D 1m, gravity acceleration g D 10m=s 2 and initial position q D OE1 0 0 T . The pendulum is released from rest and swings clockwise because of the gravity force. An obstacle is placed 503 such that a unilateral constraint is enforced as p 2=2 6 x ? u > 0, and is then expressed at velocity level. The coefficient of restitution is chosen as e D 0.8.
In the numerical tests, the nominal time step size is 10 3 s; 1 D 0.8 for the nonsmooth generalized-˛method, Â D 1 for the Moreau-Jean time stepping method and D 1 andˇD 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 k k 2 norm is used for the projection method, as where N 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 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 D 1s. 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 D 10 7 s. The order of accuracy remains first order with larger errors with the Moreau-Jean method compared with the other two.

Multiple impacts of a crank-slider mechanism
This section deals with the analysis of multiple impacts of a planar crank-slider 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.

Rigid crank-slider 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 D 0.153m, length of the connecting rod l 2 D 0.306m, width of the slider 2a D 0.05m, length of the slider 2b D 0.1m and distance between the guide surfaces d D 0.052m. Initial configuration parameters are Â 1 D 0rad, Â 2 D 0rad, Â 3 D 0rad and both the upper and lower clearances between the slider and the guide surfaces are c D 0.001m.
The mass and inertia of the crank: m 1 D 0.038kg and J 1 D 7.4 10 5 kg m 2 ; the connecting rod: m 2 D 0.038kg and J 2 D 5.9 10 4 kg m 2 ; the slider: m 3 D 0.076kg and J 3 D 2.7 10 6 kg m 2 .
The crank is released with an initial speed ! 1 D 150rad=s and the initial speed of the connecting rod is set as ! 2 D 75rad=s. The gravity acceleration g D 10m 2 =s and the restitution coefficient The model is first simulated with time step size h D 10 4 s. In order to study the influence of the step size, the Moreau-Jean method with a finer step size of h D 10 5 s is also considered. Other numerical parameters are as: 1 D 0.8 for the nonsmooth generalized-˛method, Â D 1 for the Moreau-Jean method, and D 1 andˇD 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].

Flexible crank-slider 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 connecting rod as in the rigid crank-slider mechanism, the width and the height of the crank are chosen as b 1 D 0.0064m and d 1 D 0.005m; the width and the height of the connecting rod are b 2 D 0.0032m and d 2 D 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 ! 1 D 150rad=s and ! 2 D 75rad=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 D 10 4 s, then a finer step size of h D  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.

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 508 Q.-Z. CHEN ET AL.
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 D 2m=s 2 and e D 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: where n is the dimension of A: n D 4 in this example. The coefficients˛i .t / are as follows: Ái , Ái . 2) works only for the smooth situation. After each discontinuity due to the impacts, Q 0 and t 0 should be updated. Here, the impact time instants are numerically obtained by solving q 1 D 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 D 10 4 s, the exact solution of the system is computed till t D 3.46s, just before the accumulation and the closed contact situation.