Numerical integration of multibody dynamic systems involving nonholonomic equality constraints

This work considers a class of multibody dynamic systems involving bilateral nonholonomic constraints. An appropriate set of equations of motion is employed first. This set is derived by application of Newton’s second law and appears as a coupled system of strongly nonlinear second-order ordinary differential equations in both the generalized coordinates and the Lagrange multipliers associated with the motion constraints. Next, these equations are manipulated properly and converted to a weak form. Furthermore, the position, velocity and momentum type quantities are subsequently treated as independent. This yields a three-field set of equations of motion, which is then used as a basis for performing a suitable temporal discretization, leading to a complete time integration scheme. In order to test and validate its accuracy and numerical efficiency, this scheme is applied next to challenging mechanical examples, exhibiting rich dynamics. In all cases, the emphasis is put on highlighting the advantages of the new method by direct comparison with existing analytical solutions as well as with results of current state-of-the-art numerical methods. Finally, a comparison is also performed with results available for a benchmark problem.


Introduction
This work is concerned with dynamics of multibody mechanical systems subject to motion constraints. Such constraints arise frequently in modeling the dynamics of a plethora of engineering systems and need a special treatment [1][2][3]. In particular, the class of systems examined involves nonholonomic bilateral constraints. That is, nonintegrable equality constraints on their generalized velocities. In contrast to integrable (or holonomic) constraints, which limit both the velocities and the positions of a system, nonholonomic constraints impose restrictions on the velocity level only, without limiting the possible positions on their configuration space. In classical applications, such constraints arise in the presence of rolling or sliding motions. They lead to quite interesting dynamics, not exhibited by holonomic systems and represent a central subject of Nonholonomic Mechanics [4,5].
The main emphasis of this work is placed on developing and presenting the essential steps of a new systematic methodology, leading to a robust, accurate and efficient numerical integration of the equations of motion of mechanical systems involving nonholonomic constraints. Classical numerical methodologies, treating the equations of motion as a coupled system of differential-algebraic equations (DAEs), lead to severe difficulties and inaccuracies, related to response stability issues [6]. Specifically, small errors in the initial values as well as errors developed through numerical integration are amplified and lead to an erroneous response as time goes on, independently on the accuracy level imposed on the calculations. Most of the times, increasing these accuracy levels helps in just extending the acceptable accuracy of the solution for a little longer time duration but this is not sufficient to prevent an eventual drift from the real solution. This prompted the initiation of many research efforts, focusing toward eliminating or reducing these undesirable consequences. Besides consideration of dissipative time integration schemes [1,3], this included efforts to properly scale the governing equations and constraints or to perform appropriate velocity projections [7][8][9].
In this work, the effort to eliminate the stability problems arising during integration of the equations of motion of constrained systems starts from their theoretical foundation and continues with the development and application of suitable numerical methodologies. More specifically, the essential reason for achieving an accurate and reliable numerical solution is the utilization of an appropriate set of equations, derived by a consistent application of Newton's law of motion [10]. Their derivation is based on classical concepts and tools of Analytical Dynamics and leads directly to a set of second-order ordinary differential equations (ODEs) in both the generalized constraints and the Lagrange multipliers, associated with the motion constraints [11]. This eliminates the wellknown pathologies associated with the classical DAE formulations. Moreover, the robustness of the new numerical scheme is reinforced by first putting these equations in a suitable weak form and then treating the generalized coordinates, velocities and momenta as independent quantities, leading to a convenient threefield weak formulation [12]. This follows earlier mixed multi-field methods in computational mechanics, based on the de Veubeke-Washizu variational principle [13,14]. Finally, the process is completed by developing and applying a suitable temporal discretization scheme.
A similar path was followed in a recent work of the authors [15]. In that work, the emphasis was placed on mechanical systems with holonomic constraints. Instead, the present work focuses on systems involving nonholonomic constraints. This leads to some important deviations in both the analytical formulation and the numerical scheme developed, which are noticed and explained in detail in Sects. 2-4. In addition, the set of unknowns arising in each time step of the present treatment includes both the generalized velocities and the time derivatives of the Lagrange multipliers. In contrast, the method applied in [15] was based on an augmented Lagrangian approach and treated the generalized velocities and the time derivatives of the Lagrange multipliers at different solution levels. Finally, all the numerical examples investigated in that work involved holonomic constraints, while the examples examined in this study involve nonholonomic constraints.
An outline of the present work is as follows. First, the set of equations of motion employed is presented briefly in Sect. 2, giving emphasis on explaining and modeling the effect of the nonholonomic constraints. Then, the essential details needed for putting these equations in a suitable three-field weak form are included in Sect. 3. This provides a solid foundation for developing an accurate and efficient methodology for the temporal discretization of the equations of motion, which is done in Sect. 4. Next, a selected set of numerical results is included in Sect. 5 for three typical mechanical examples. For the first two of them, analytical solutions are available, while the third example is a benchmark problem. In this way, direct comparison of the numerical results provides strong evidence and demonstrates the advantages of the new method, over classical DAE formulations. Finally, a synopsis of the study, together with possible future extensions, is presented in Sect. 6.

Equations of motion
The equations of motion employed for the class of systems examined in this work have been obtained and presented in a previous publication of the authors [11]. In this section, a brief summary of the essential steps performed and the main results derived in that work are presented in order to set up the notation and enhance the completeness of this study. Specifically, by adopting the classical viewpoint of Analytical Dynamics, the motion is described by a set of generalized coordinates q ¼ ðq 1 . . .q n Þ, representing motion of a fictitious point p on a manifold M, as a function of time t [16][17][18]. Moreover, the tangent vector v to the motion curve at point p, known as generalized velocity, belongs to an n-dimensional vector space T p M [4]. Therefore, if B e ¼ f e 1 . . . e n g is a basis of T p M, vector v can be put in the compact form v ¼ v i e i by using the summation convention on repeated indices [18]. Also, the elements of the corresponding cotangent (or dual) space T Ã p M, known as covectors, represent generalized momenta [17].
In dynamics, the inner product is expressed in terms of the components g ij of the metric tensor at point p, which are determined through consideration of the kinetic energy of the system.
When there are no motion constraints, a minimal and independent set of generalized coordinates can be selected. In general, the geometry of the corresponding configuration manifold M is non-Euclidean. For this reason, the exact solution path is determined by application of Newton's second law in the form The term in the right-hand side is a covariant derivative along a path on M, with tangent vector v [10]. It is expressed in the component form with i; j; k ¼ 1; . . .; n, where the quantities K k ji are known as affinities. They are components of the affine connection r and dictate the transition from a tangent space to any neighboring tangent space of manifold M [19]. Moreover, covectors p represent the generalized momenta and the applied forces, respectively [10]. In fact, direct application of Eq. (1) leads to Then, by introducing the classical matrix notation.
where vector hðq; _ qÞ includes the inertia terms related to the affinities K k ij . The picture changes drastically when the system is subjected to motion constraints. In particular, the class of mechanical systems examined is assumed to be subject to a set of k scleronomic nonholonomic and linearly independent constraints, with general form where A ¼ ½a R i ðqÞ is a known k Â n matrix. The manipulation of these constraints is quite critical for the derivation of the equations of motion. Essentially, they provide the tools to decompose the original configuration manifold M locally into an (n À k)dimensional manifold M A , described by an appropriate set of n À k minimal coordinates, plus k one-dimensional manifolds M R , one for each constraint. This is achieved by using Eq. (6) for defining a suitable set of linear operators, acting between the tangent and dual spaces of these manifolds. First, these operators help in selecting the geometric properties (i.e., metric and connection) of manifolds M A and M R , so that Newton's law of motion is transferred from M to M A and to M R (for R ¼ 1; . . .; k) in an invariant form, similar to that expressed by Eq. (2) [10]. In addition, they help in splitting each element of the cotangent space T Ã p M in a unique manner. More specifically, the elements of the special class of Newton covectors on M, defined by are expressed in the form The Due to the presence of the motion constraints, the original set of coordinates q ¼ ðq 1 . . . q n Þ becomes redundant. Consequently, condition (2) is no longer true and is replaced by instead. Then, by taking Eqs. (9) and (10) into account, Next, substitution of Eqs. (3) and (4) in Eq. (7) yields with components Likewise, each Newton covector h s Ã R on the constraint manifolds M R (for R ¼ 1; . . .; k) takes the form with components In the last two equations and in the sequel, the convention on repeated indices does not apply to index R. Moreover, during the solution process, the terms are determined through a projection along special directions c R on T p M, for each of the constraints. In particular, the components of the n-vector c R are chosen so that Finally, each coefficient c RR is selected so that the corresponding term in Eq. (15) represents a corrective force applied on the figurative point, when a velocity violation tends to develop along direction c R . For instance, if the applied forces depend on the generalized velocity v, a convenient choice is Next, by combining Eqs. (11)- (15), it easily turns out that for i ¼ 1; . . .; n. Therefore, after enhancing the notation introduced by Eq. (5) with and the array f , as determined by Eqs.
Equation (20), or equivalently Eq. (22), corresponds to a set of n second-order ODEs, involving the n þ k unknowns q i and k R . A complete mathematical formulation is obtained by including the k equations of the constraints. Originally, these are expressed by Eq. (6) but can eventually be put in the second-order ODE form through appropriate mappings from M to M R , forcing _ w R to become zero, for R ¼ 1; . . .; k. In contrast to cases involving holonomic constraints, there is neither a w R term in Eq. (23) nor a k term in Eq. (22). Complete details on the formulation presented in this section as well as an extensive comparison with previous formulations can be found in the original publication [11]. In brief, the present theoretical approach brings significant advantages when compared to other approaches applied so far in the field of Multibody Dynamics. These advantages are related to the physically consistent and correct elimination of the singularities associated with the usual sets of highindex DAEs of motion [1][2][3]. More specifically, by treating the motion constraints as an integral part of the overall process of deriving the equations of motion led to a set of dynamic Lagrange multipliers, expressed by Eq. (15). In this way, each constraint introduces effective inertia terms in Eqs. (22) and (23). Consequently, both of these equations appear now in a pure second-order ODE form, from the onset and in a natural way. In addition, the extra terms appearing in these equations are evaluated in a systematic and analytical manner through exact application of Newton's second law on configuration manifolds possessing general geometric properties. As a result, there is no need for an ad hoc selection of their values for numerical stabilization or for scaling the constraints and the equations of motion. Some extra advantages are also realized during the development of the numerical integration process, which is started by first putting the equations of motion in a suitable multifield weak form and is then integrated by performing an appropriate temporal discretization, as explained in Sects. 3 and 4, respectively. All these advantages are further supported and illuminated by the numerical results presented in Sect. 5.
3 Three-field weak form of the equations of motion for arbitrary time instances t 1 and t 2 . Then, by performing lengthy manipulations, including a usual operation involving integration by parts (see ref. [12] for details), the above relation leads eventually to the following weak form of the equations of motion The quantities w i and dv i represent variations. They are evaluated by taking into account that the variation of any scalar function f ðqÞ on manifold M is defined as the derivative of f along an arbitrary vector w of T p M, according to the rule Therefore, for each holonomic coordinate this yields w i ¼ dq i , while a little more involved relation is obtained in case of anholonomic coordinates [12,19]. In this way, the variation of the velocity components is defined by dv i ¼ wðv i Þ. Moreover, the total time derivative Da R i Dt, appearing in the second line of Eq. (25), is determined by with partial derivatives Finally, the terms s i jk represent components of the torsion of the connection selected on manifold M , while the quantities r i jk are arbitrary constants, with r i kj ¼ Àr i jk [12]. A similar treatment of each constraint equation, expressed by Eq. (23), leads first to which, after integration by parts, becomes Next, in order to exploit more advantages of the weak formulation, the position, velocity and momentum variables are considered as independent quantities [13,20,21]. For this, a new velocity field t is introduced on manifold M, which is forced to become identical to the true velocity field v ¼ dq . dt in an average rather than in a pointwise sense through the relations where dp i are arbitrary constants. Likewise, a new quantity l R is introduced for each constraint, which is forced to become identical to _ where dr R are arbitrary constants. In addition, these last sets of equations are enhanced by an accompanying set of conditions, consisting of where p i and r R are appropriate Lagrange multipliers, leading to independent variations dt i , dv i and dl R , d _ k R , respectively. Therefore, by augmenting Eq. (25) with Eq. (32) and performing another round of manipulations, it can be shown (see [12] for details) that this leads to the following variational equation A similar combination of Eqs. (29) and (33) yields By construction, the variations w i (representing dq i or d# i for a true coordinate q i or a pseudo-coordinate # i , respectively), dk R , dt i , dl R , dp i and dr R are independent for all i ¼ 1; . . .; n and R ¼ 1; . . .; k. Consequently, Eqs. (20)- (25) can be employed to construct a dynamically equivalent set of first-order ODEs. This was done in an earlier work, leading to the classical Hamilton's canonical equations [12]. Here, the weak form expressed by these equations is used, instead, as a foundation for performing an appropriate temporal discretization of the equations of motion, as explained in the following section.

Temporal discretization of the equations of motion
First, the presence of the motion constraints causes a non-flatness of the configuration space of the class of systems examined. This, in turn, makes advantageous the utilization of a special type of curves in performing the temporal discretization. In this work, these curves are selected to be the autoparallels, corresponding to the ''straightest'' curves on the configuration space [19,22]. Moreover, following earlier approaches [15,23], all the variations in Eqs. (30)-(35) are assumed to remain constant within each time interval [t m ; t mþ1 ], while the corresponding variables vary linearly in time. Then, it can be shown that these conditions convert Eqs. (30) and (31) into relations with the following general form and where Dt ¼ t mþ1 À t m is the time step. In particular, the function gðt mþ1 ; t m Þ in Eq. (37) is determined in terms of the components of the connection selected.
For instance, in the special case of a Euclidean configuration space [23], where K k ij ¼ 0, it turns out that Equations (37) and (38) represent update formulas for the generalized coordinates and the Lagrange multipliers, respectively. In addition, collecting the terms in Eqs. (34) and (35) multiplied by dt i and dl R yields and respectively. These equations can be used in order to update the values of the generalized momenta p i and r R , respectively, when necessary. Finally, collecting the remaining terms of Eqs. (34) and (35), multiplied by the variations w i and dk R , yields and respectively. In the sequel, a summary is provided for the essential information needed in understanding the basic features of the new numerical scheme. As a result of the three-field formulation performed in this work, the set of unknowns of the resulting mathematical problem is enhanced and includes the generalized coordinates q i and k R (i.e., the Lagrange multipliers are also considered as generalized coordinates from hereon), together with the corresponding generalized velocities t i and l R as well as with the generalized momenta p i and r R . In total, these quantities give rise to a set of 3ðn þ kÞ unknowns. Their determination is accomplished by employing the system of Eqs. (30), (31) and (39)-(42). Here, this is achieved by developing and applying an implicit numerical scheme in order to carry out the required temporal discretization. As usual, the outcome of this process is a system of nonlinear algebraic equations, whose solution provides the values of the unknown quantities at the end t mþ1 of the time step considered in terms of their known values at earlier time steps. In fact, this set of equations is solved by applying a block-type iterative technique within each time step, as explained next.
First, it is assumed that the values of all the unknowns, but the components of the weak velocity vectors, defined by t ¼ ð t 1 Á Á Á t n Þ T and l ¼ ð l 1 Á Á Á l k Þ T , are fixed. Then, an involved set of n þ k nonlinear algebraic equations is obtained through application of Eqs. (41) and (42), which can be put in the general form gðt mþ1 ; l mþ1 ; t m ; l m ; q m ; k m ; p m ; r m Þ ¼ 0: In the last expression, t mþ1 , l mþ1 and t m , l m represent the values of t, l at times t mþ1 and t m , respectively. Moreover, a similar meaning is also given to the vector quantities, p ¼ ð p 1 Á Á Á p n Þ T and r ¼ ð r 1 Á Á Á r k Þ T . Following common practice, Eq. (43) is then solved by applying a Newton-Raphson approach, with respect to the unknown To achieve this, given an estimate x ' mþ1 , a corrected value x 'þ1 mþ1 is obtained, according to where the correction Dx ' is determined by solving the linearized problem resulting by substituting Eq. (45) into Eq. (43), with Jacobian matrix and residual vector ; t m ; l m ; q m ; k m ; p m ; r m Þ: ð48Þ The resulting system of equations involves only n þ k of the original unknowns. Each problem, as expressed by Eq. (46), is solved by employing a direct linear solver. The computations are stopped when the set of weak velocities x mþ1 and the constraint equations are satisfied up to a prespecified accuracy or the iterations exceed a critical number. In the latter case, the time step Dt is reduced and the process is restarted. Next, the values of the generalized coordinates q mþ1 and k mþ1 are determined through a direct update, based on Eqs. (37) and (38), respectively. The iteration process is completed when the residual in the right-hand side of Eq. (46) becomes also sufficiently small. Otherwise, the process is repeated after decreasing the time step. Finally, the new values of the momentum variables p mþ1 and r mþ1 can also be obtained after the iterations are finished by using the subsystem of equations resulting from application of Eqs. (39) and (40), respectively.
For better clarity, the numerical implementation of the algorithm employed for solving the discretized set of equations of motion, given by Eqs. (41) and (42), is presented next. Here, the analysis is restricted to systems possessing rigid bodies only. In this case, through the choices r ' ij ¼ 0 and s ' ij ¼ K ' ij , with the affinities as determined in [10,22], Eq. (41) is replaced by This equation can be written in a convenient matrix form where the terms _ a R i are defined by Eq. (36). Moreover, the action of operator A is equivalent to the identity or the corresponding rotation matrix of a rigid body when it applies to its translational or rotational degrees of freedom, respectively. In a similar manner, Eq. (42) is also put in a matrix form Therefore, application of the classical trapezoidal rule to Eqs. (50) and (52) leads eventually to respectively. Equations (54) and (55) where R s v and R s l is the numerical residual of Eq. (54) and (55), respectively. Then, the square of the norm of the total residual is evaluated in the form Likewise, the square of the norm of the residual of the equation constraints, given by Eq. (6), is obtained by For a given tolerance level, the calculations are assumed to have converged when both jjR s T jj and jjR C jj are smaller than the selected tolerance value.

Numerical results
In this section, a set of characteristic results is presented for three challenging mechanical examples.
Special emphasis is placed on highlighting the advantages of the integration scheme developed. For this, the attention is focused on comparing numerical results with analytical results as well as with similar results arising by employing standard DAE solvers for two classical nonholonomic systems. Finally, comparison with results for a benchmark problem involving nonholonomic constraints is also performed. Specifically, the numerical solutions captured by the method presented in this work (indicated by the label NM-ODE) are compared with existing analytical solutions [5] and results of a well-known benchmark problem [24,25]. Moreover, the results are compared with similar results, obtained by using two state-ofthe-art commercial codes [26,27]. More specifically, these codes set up the equations of motion as a system of high-index DAEs and solve them numerically, by employing classical integration schemes, based on backward differentiation formulas. In particular, the GSTIFF and DASPK method were selected in solving the equations of motion by ADAMS and Motion Solve. Finally, results obtained by applying the new method after setting in Eq. (20), for R ¼ 1; . . .; k, are also included. In this way, the set of equations employed is reduced to the following system of equations of motion which is identical to that employed by multibody dynamics DAE formulations [28][29][30]. Moreover, the value is also selected in Eq. (23). In the sequel, this modified set of equations is referred to as MM-DAE.

Sphere rolling inside a cylinder
In the first example, motion of a sphere of unit mass and radius a, rolling on the rough inside wall of a fixed vertical cylinder with radius b, is investigated. The sphere and the cylinder, together with a fixed coordinate system Oxyz and a rotating coordinate system Ox 1 y 1 z 1 , are shown in Fig. 1. Moreover, the action of gravity is along the negative Oz axis. Then, if f e x e y e z g and f e 1 e 2 e 3 g are appropriate orthonormal bases for the coordinate frames Oxyz and Ox 1 y 1 z 1 , respectively, it easily turns out that where r ¼ r OG is the position vector of the center of mass G of the sphere, I 3 is the 3 Â 3 identity matrix and E ¼ e z e T z : In addition, if C is the contact point of the sphere with the cylinder, the condition of rolling without slipping of the sphere on the cylinder is expressed by where x is the angular velocity of the sphere and r GC ¼ ae 1 . Therefore, after relating the components of the angular velocity of the sphere in the Oxyz and a body frame (fixed in the sphere), through the corresponding rotation matrix R, i.e., the rolling condition is put eventually in the form v G þx r GC ¼ v G Àr GC ðR XÞ ¼ 0: As usual,x represents the 3 Â 3 antisymmetric matrix . Equivalently, after performing some straightforward operations, this leads to This condition belongs to the class of constraints expressed by Eq. (3) with An exact solution for this problem can be found in pp. 95-98 of ref. [5]. Next, a selected set of numerical results is presented in Fig. 2. These results were obtained for the parameter values a ¼ 0:4 m, b ¼ 1 m, z 0 ¼ 1:224 m, x z0 ¼ Fig. 2 Rolling sphere on a fixed vertical cylinder: a vertical displacement z of the sphere center for 0\t\10 sec, b displacement z for 10 sec \t\20 sec, c displacement z for 95 sec \t\100 sec, d mechanical energy for 0\t\5 sec, e norm jjR C jj of the residual of the constraint equations for 0\t\5 sec, f time step for 0\t\5 sec f f x 30 ¼ À5 rad/s and a gravity acceleration g ¼ 9:81 m s 2 . Also, selecting the amplitude A ¼ 1 m and the phase E ¼ 0 in the exact expression of x 2 , presented by Eq. (2.24) in [5], leads to the initial angular velocity values x 20 ¼ 0 and x 10 ¼ À20:265 rad/s. First, in Fig. 2a are presented results for the time history of the vertical displacement z of the sphere center. According to the analytical predictions, the sphere center executes a harmonic oscillation. This is obvious from the continuous curve in Fig. 2a, representing the exact solution. Also, identical results were obtained by the new numerical method (NM-ODE), represented by the full dots. However, the picture arising by employing the DAE solvers is quite different. First, only Motion Solve run and provided results for this example, while attempts to run it with ADAMS were not successful. In any case, the results obtained by Motion Solve, using the higher tolerance value, start deviating from the very beginning. In addition, decreasing the tolerance value causes a temporary improvement, which is quickly lost as time increases. This is verified by the results presented in The gradual worsening in the quality of the solution obtained by the DAE solvers, as time goes on, can also be realized by looking at the mechanical energy of the system. Since there is no slipping event, the system examined is conservative. This means that its mechanical energy should be preserved. This is clearly the case for the results obtained by the new method, as illustrated in Fig. 2d. On the contrary, the results obtained by the DAE solvers exhibit a gradual and noticeable change in the mechanical energy. Similar behavior has also been observed and reported in other applications [8,9,15].
Finally, some more light is shed on the performance of the numerical methods employed by looking at the time variation of the norm jjR C jj of the residual of the constraint equations and the corresponding time steps. In all cases, an initial value of 0.01 s was selected for the time step. First, the results in Fig. 2e and f indicate that the new method converges to the correct solution by keeping the time step constant to the initial value of 0.01 s throughout the calculations, for a prespecified value of the tolerance jjR s T jj. On the other hand, the results in Fig. 2e demonstrate that in order for all the DAE numerical schemes applied to converge to a solution, independently of how good or bad their ultimate accuracy is, the norm jjR C jj should be reduced to extremely small levels eventually. As a consequence, the corresponding time step must also be reduced dramatically (i.e., more than two orders of magnitude), as verified by the results included in Fig. 2f. This demonstrates the numerical efficiency of the new method.

The falling rolling disk on a horizontal plane
In the second example, dynamics of a circular rigid disk, rolling without sliding on a horizontal plane, is examined. The disk has a unit mass and its radius is a. Moreover, the diametral and polar mass moments of inertia of the disk with respect to its center of mass G are known. The position of the disk is completely specified by five generalized coordinates, defined in Fig. 3. Namely, the coordinates x and y of the lowest point A of the disk, which is in contact with the horizontal plane, with respect to a frame Oxyz, together with three orientation angles h, u and w. The frame Oxyz is fixed on the ground, so that the gravity force is along the negative Oz axis, with gravity acceleration g.
If e x e y e z È É and e n e g e f È É are orthonormal bases for the fixed coordinate frame Oxyz and the moving frame Ongf, shown in Fig. 3, respectively, it turns out that where R is the rotation matrix between the fixed frame Oxyz and a frame fixed on the disk, while Then, if the disk center of mass is located at point G r GA ¼ Àae g and the condition of rolling of the disk without slipping on the ground is expressed by where x is the angular velocity of the disk. Therefore, after some manipulation, the rolling condition can eventually be written in the form which belongs to the class of constraints expressed by Eq. (3). In all the subsequent numerical examples, the disk radius was selected as a ¼ 1 m, the diametral and polar mass moment of inertia of the disk with respect to its center G was 0.25 and 0:5 kg m 2 , respectively, while the gravity acceleration was chosen as g ¼ 9:81 m s 2 . Also, without a loss in generality, the initial values of the coordinates x, y, u and w were selected to be zero. Finally, for compatibility in the notation with [5], the angular velocity x of the disk is expressed in the moving frame Ongf as so that the components of the angular velocity in the moving frame are A quite complete analysis of the dynamics of the disk under the action of gravity, together with an investigation of certain stability issues, is presented in pp. 55-60 of ref. [5]. In general, an arbitrary set of initial conditions leads to complex motions, with a variable angle h. However, certain combinations of the system initial parameters lead to special motions. One such category of motions are the cyclic motions, which are characterized by constant h, q and r. Namely, for some special cases with h ¼ h 0 6 ¼ 0, the point of the disk in contact with the ground undergoes a circular orbit. Instead, when h 0 ¼ 0, the disk can follow a uniform rectilinear rolling or a uniform spinning about a fixed vertical diameter [5].
First, in Fig. 4 are presented results obtained for h 0 ¼ 1 , p 0 ¼ q 0 ¼ 0 and three different initial values r 0 . For this set of initial conditions, the results presented indicate that complex motions arise, with time varying h, p, q and r. In all cases examined, the results of the new method are found to be virtually indistinguishable from those obtained by the analytical solution, verifying the accuracy of the new method. Moreover, in terms of dynamics, the results demonstrate that as the absolute value of r 0 increases, the magnitude of the variations of h, p, q and r decreases gradually. Eventually, trajectories and dynamics with h, q and r approaching constant values are reached, for sufficiently large absolute value of r 0 , corresponding to cyclic motions. The dynamics of the motions presented in Fig. 4 is illuminated in a better way by Fig. 5a, depicting the corresponding trajectories of the disk contact point on the horizontal plane. The orbits obtained exhibit a geometric complexity for a small absolute value of r 0 . This complexity is diminished gradually as the absolute value of r 0 increases and the orbit tends to Fig. 3 A falling rolling disk on a horizontal plane g 123 reach a clean circular shape. In addition, the results of Fig. 5b, obtained for the largest absolute value of r 0 , that is r 0 ¼ À2:25 rad=s, with h 0 ¼ 0:5 o , indicate that the radius of the circular trajectory increases as the value of h 0 decreases. In this way, the motions examined reach another special type of motion, known as uniform rectilinear rolling, in the limit value h 0 ¼ 0 [5].
Cyclic motions, with h, q and r constant, can be captured from the start, through an appropriate choice of the initial conditions. For instance, the results presented in Fig. 6 were obtained for h 0 ¼ 20 o , p 0 ¼ 0 and q 0 ¼ 1 rad=s, while the value of r 0 was selected by using Eq. (2.9) of [5]. Once again, the new method (NM-ODE) gives identical results with the analytical solution. In addition, the pathologies of the DAE   Fig. 6a and b, respectively. These results demonstrate a gradual and monotonically increasing loss of accuracy in the computations, as time goes on. Furthermore, decreasing the tolerance value causes an improvement, which is only temporary. Similar behavior was also noticed for the other two angular velocity components p and q. As a consequence of these computational errors, the trajectory of the disk contact point on the horizontal plane starts deviating from the circular orbit, as shown in Fig. 6c, for the higher tolerance value. The result for the lower tolerance value is not included since it leads to similar behavior and just complicates the figure. In any case, the results indicate a gradual decrease in the radius of the orbit. This is justified by looking at the corresponding time histories of the disk mechanical energy, which are included in Fig. 6d. Obviously, the new method predicts a constant energy state, which is correct, while omission of the extra terms in the governing equations and the constraint equations leads to a gradual loss of energy.
Finally, the study focused on motions with h ¼ 0, where the disk executes a uniform rectilinear rolling. In the first set of results, presented in Fig. 7, the initial values h 0 ¼ 0, p 0 ¼ 0 and q 0 ¼ 0 were selected. At the same time, the initial value r 0 was chosen based on results of a stability analysis. Specifically, a critical value r cr was selected, based on Eq. (2.14) of [5]. For the chosen combination of parameters, this numerical value was r cr ¼ 1:808 rad=s. Results were also obtained for a higher value (r 0 ¼ 2r cr ) and a lower value (r 0 ¼ 0:9 r cr ), corresponding to a supercritical and a subcritical case, respectively. For each of these cases, the trajectory of the disk contact point on the horizontal plane and the corresponding history of the  Fig. 7. In all cases, the results of the new method were found to coincide with those obtained by the analytical solution, once again. Specifically, the trajectory of the disk contact point on the horizontal plane was a straight path, while the angle h remained zero for all times. In fact, this was possible even after selecting an arbitrarily large time step and tolerance value, since the specific motion (translation along a line plus rotation about a fixed axis) takes place along an autoparallel curve of the configuration space. However, the picture was quite different for the results obtained by using ADAMS. Namely, significant deviations were observed to occur, especially in the critical case and became even more pronounced in the subcritical case. Moreover, it is apparent that decreasing the tolerance level in the calculations does not prevent initial errors from increasing as time increases.
Some important and rapid changes were observed to occur in the disk motion by imposing a small perturbation in the value of p 0 . For instance, in Fig. 8 is presented a similar set of results, obtained for the same set of parameters as in the previous example, but with p 0 ¼ À0:1 rad=s, instead. In this case, the trajectory exhibits some noticeable deviations from a straight line. At the same time, the inclination angle h Fig. 7 Uniform rectilinear rolling motions of the disk (p 0 ¼ 0). Trajectory of the disk contact point on the horizontal plane and history of angle h for: a and b a supercritical case, with r 0 ¼ 2r cr , c and d the critical case (r 0 ¼ r cr ), e and f a subcritical case, with r 0 ¼ 0:9 r cr also takes nonzero values. These deviations are relatively small in the supercritical range but are amplified quickly with a reduction in the value of parameter r 0 , especially for values smaller than r cr , where the changes are quite rapid and the motion deviates quickly and significantly from the rectilinear motion. Once again, the results of the new method coincide with those obtained by utilizing the analytical solution, in all cases examined. In contrast, all the predictions of the DAE method employed are quite unreliable. Moreover, these predictions are subject to qualitatively similar errors with the case p 0 ¼ 0.

The uncontrolled bicycle
In the final example, results of the new method are compared with results available for a benchmark problem [31]. More specifically, response of an uncontrolled bicycle problem presented first in [24] is investigated. The corresponding mechanical model supercritical case, with r 0 ¼ 2 r cr , c and d the critical case (r 0 ¼ r cr ) and e and f a subcritical case, with r 0 ¼ 0:9 r cr consists of four rigid members (i.e., the rear and front frame assembly, together with the rear and front wheel), as shown in Fig. 9. All the parameters of this model are taken from [31].
Assuming that both wheels roll on the ground without slipping, introduces nonholonomic constraints in the formulation. Eventually, the model examined possesses three degrees of freedom. Namely, the tilt or roll angle, the steering angle and the forward displacement of the bicycle. The bicycle is under gravity forcing, acting along the z axis. Also, at the starting position, the bicycle is in a vertical position and the steering is straight. Then, it is given an initial forward velocity, together with some roll velocity as a perturbation. According to a stability analysis of the reference position [24,25], the subsequent motion of the bicycle is stable only when the forward speed lies within an interval of two critical values, say v c1 and v c2 .
In particular, when the forward speed of the bicycle is smaller than v c1 , the instability is static in nature [24,25]. In accordance to the benchmark problem, selecting a forward speed of 4.0 m/s and a roll angle speed of 0.05 rad/s leads to Figs. 10a-c for the time histories of the forward velocity, the roll angle and the steer angle, respectively, over the first 20 s of the motion. Likewise, in Fig. 10d is shown the corresponding mechanical energy E m of the system. In all cases, the benchmark results are represented by continuous curves, while the results of the new method are indicated by broken lines with full dots.
The virtual coincidence of the results illustrates the accuracy of the new method. Moreover, the corresponding percentage of variation of the mechanical energy E m [31], defined by was calculated to be 5:5 Á 10 À4 , which is much better than the value of 17:0 Á 10 À4 , obtained by the reported benchmark results. Similarly, the results presented in Fig. 11 were obtained for an initial forward velocity value of the bicycle lying between the two critical values v c1 and v c2 , leading to an asymptotically stable solution [24,25]. More specifically, these results were obtained for a forward speed of 4.6 m/s and a roll angle speed of 0.50 rad/s. In addition, the corresponding percentage of variation of the mechanical energy was calculated to be 5:2 Á 10 À4 , which is better than the value of 9:5 Á 10 À4 determined by the reported benchmark results. Once again, the results illustrate the numerical accuracy and stability properties of the new method.
Finally, similar conclusions can be drawn by considering the set of results presented in Fig. 12. These results were obtained for an initial forward velocity value of the bicycle greater than v c2 , leading to an oscillatory unstable solution [24,25]. More specifically, these results were obtained for a forward speed of 8.0 m/s and a roll angle speed of 0.05 rad/s. In this case, the corresponding percentage of variation of the mechanical energy was calculated to be 0:105 Á 10 À4 , which is about the same with the benchmark value of 0:106 Á 10 À4 .
In closing this section, it is mentioned that results obtained by application of all the DAE methodologies applied in the earlier examples were quite different than the benchmark results, in all three cases examined for the bicycle example. Also, qualitatively similar results were obtained for other typical mechanical systems subject to nonholonomic constraints, arising from a sliding rather than a rolling action (e.g., knifeedge, Chaplygin's sleigh [4,5]).

Synopsis and extensions
The basic ingredients of a new method developed for the numerical integration of the equations of motion governing the behavior of multibody mechanical systems, involving nonholonomic motion constraints, were presented in the first part of this work. This   Initially, an appropriate set of equations of motion was employed, consisting of a coupled and strongly nonlinear system of second-order ODEs for both the generalized coordinates and the Lagrange multipliers. These equations include suitable terms, causing an automatic stabilization and scaling of the governing equations and the equations of the constraints. These terms were obtained by a systematic and consistent application of Newton's law of motion, avoiding ad hoc and incorrect selections of parameters. Then, these equations were put eventually in a three-field weak form, by treating the generalized coordinates, velocities and momenta as independent quantities. Finally, a suitable temporal discretization scheme was developed for the class of systems examined, based on this solid theoretical foundation.
The robustness, accuracy and efficiency of the new numerical method were demonstrated in the second part of this study by applying it to typical mechanical examples of nonholonomic mechanics. First, dynamics of the classical example of a sphere rolling inside a cylinder was examined. Then, dynamics of a falling rolling disk on a horizontal plane was investigated in detail. For both of these examples, analytical solutions are available. In all cases examined, the results of the new method were found to be virtually coincident with those obtained by the analytical solution, validating the excellent accuracy and stability properties as well as the numerical efficiency of the new method in the best possible way. Moreover, special emphasis was placed on comparing the results of the new method with those obtained by application of classical DAE approaches. This brought up and illustrated the classical instability problems arising when using such methods. Similar conclusions were also drawn by considering a well-known benchmark problem, referring to response and stability of an uncontrolled bicycle model. Once again, comparison of the numerical results verified the accuracy and stability properties of the new method.
The new method can easily be applied to more complex and challenging mechanical systems, possessing a large number of degrees of freedom [32]. In such cases, it is quite possible to run into cases involving the presence of redundant constraints and/or the appearance of singular configurations during the motion of the system. Then, in order to be able to accommodate such numerically difficult situations, the new method should first be extended to an augmented Lagrangian formulation. Finally, the present formulation can provide a solid basis to develop a numerical scheme for systems with unilateral constraints. In fact, it will be quite helpful in efforts to examine special situations arising in contact problems, where conditions of contact are fulfilled temporarily or change instantly [33][34][35]. In such cases, where the unilateral constraints can be treated as bilateral up to their violation, the constraints cannot be put in an integrable form and should therefore be treated as nonholonomic. Declarations