Ryo Kikuuwe ´bernard Brogliato

B Brogliato
email: bernard.brogliato@inria.fr

A New Representation of Systems with Frictional Unilateral Constraints and Its Baumgarte-Like Relaxation

Keywords: Baumgarte stabilization, frictional unilateral constraints, differential algebraic inclusions

A New Representation of Systems with Frictional Unilateral Constraints and Its Baumgarte-Like Relaxation

Ryo Kikuuwe, Bernard Brogliato

Introduction

For many mechanical systems composed of very stiff objects, the rigid-body formalism is a reasonable approximation suited for the purposes of simulation, analysis and control. It however requires a special care in the mathematical treatment because of its discontinuous and nonsmooth nature. If the system involves unilateral constraints with Coulomb friction, the system's dynamics is usually described by a set of inequalities and equations involving force and kinematic variables, such as position, velocity and acceleration.

One important feature of such a system is that the structure drastically changes at the transitions between the sticking, slipping, and non-contact phases of each contact point. For example, in a system with a unilateral constraint, the velocity must satisfy an inequality condition only if the position satisfies an equality condition. More specifically, a velocity-level condition, the so-called velocity Signorini condition [27, eq.( 5)], is active only when the system is in the contact phase. Thus, the representation usually includes "if" conditionals, which complicate the problem to find a unified, concise representation that describes all conditions of a system with multiple frictional unilateral constraints especially in the continuous-time domain. This feature is in contrast to a system with only an equality condition, which has a standard form of a differential algebraic equation (DAE) with a Lagrange multiplier.

Many of previous analytical studies investigated the system behavior under each particular equality condition separately. For example, Génot and Brogliato [START_REF] Génot | New results on Painlevé paradoxes[END_REF] and Song et al. [START_REF] Song | Analysis of rigid-body dynamic models for simulation of systems with frictional contacts[END_REF] analyzed systems with a single frictional unilateral constraint in the slipping phase, at which the contact force stays on the boundary of the friction cone. In such approaches, a transition from slipping to sticking causes non-existence of the solution. Event driven schemes [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF]Chap. 8] [26] also depend on representations at different kinematic levels, i.e., of position, velocity, and acceleration levels. In the studies on the development of time-stepping schemes [START_REF] Acary | Higher order event capturing time-stepping schemes for nonsmooth multibody systems with unilateral constraints and impacts[END_REF][START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb's friction[END_REF][START_REF] Anitescu | Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems[END_REF][START_REF] Stewart | An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and Coulomb friction[END_REF], derivations are performed mostly in the discrete-time domain. A continuous-time representation of the system behavior in all conditions should provide better insights to the development of computational techniques.

This paper presents two main contributions. One is a unified representation of mechanical systems subject to multiple frictional unilateral constraints. The new representation is of the form of a differential algebraic inclusion (DAI), and is expressed as a straightforward generalization of conventional representations involving Lagrange multipliers for equality constraints, which are represented as DAEs. In the proposed DAI representation, the set of the constraints is represented as an inclusion in a normal cone of the set of admissible values of the Lagrange multiplier. The normal cone used here is defined with a singular norm metric, defined as a limit of a family of normal cones with parameterized non-Euclidean norm metric.

Another contribution of this paper is a scheme with which the aforementioned DAI is relaxed into another DAI that can be equivalently converted into an ordinary differential equation (ODE). This relaxation scheme can be seen as a generalization of the Baumgarte stabilization [START_REF] Baumgarte | Stabilization of constraints and integrals of motion in dynamical systems[END_REF], which is originally for stabilizing numerical solutions of DAEs. Numerical integration of the obtained ODE provides approximate solutions of the original DAI preserving its nonsmooth nature, such as static friction. Besides that, although the relaxed representation resides in the continuous-time domain, it has a certain analogy to implicitly-discretized forms of the original DAI and includes representations at all kinematic levels as its asymptotic limits. Due to these features, the presented framework allows for understanding previous studies on time-stepping schemes [START_REF] Acary | Higher order event capturing time-stepping schemes for nonsmooth multibody systems with unilateral constraints and impacts[END_REF][START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb's friction[END_REF][START_REF] Anitescu | Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems[END_REF][START_REF] Stewart | An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and Coulomb friction[END_REF] and eventdriven schemes [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF]Chap. 7] [START_REF] Haddouni | Comparison of index-3, index-2 and index-1 DAE solvers for nonsmooth multibody systems with unilateral and bilateral constraints[END_REF] from the continuous-time point of view. As an example, Section 5.3 of this paper shows a geometric interpretation of the solution non-uniqueness of Painlevé paradox [START_REF] Génot | New results on Painlevé paradoxes[END_REF] in the continuous-time domain. This paper is organized as follows. Section 2 provides mathematical preliminaries. Section 3 overviews the treatment of DAEs with the Baumgarte method, and also discusses its relation to implicit Euler integration. Section 4 presents a new representation of mechanical systems with frictional unilateral constraints and its Baumgarte-like relaxation. Section 5 shows illustrative example applications of the presented method. Section 6 provides concluding remarks.

Mathematical Preliminaries

Throughout this paper, 0 is the zero vector or zero matrix of appropriate dimensions, In denotes the n-dimensional identity matrix, R + denotes the set of non-negative real values, and ∂X and IntX denote the boundary and the interior, respectively, of a set X .

Cones and Projections

Let us recall the definitions of cones, normal cones, and projections.

Definition 1 (Cone

) A set C ⊂ R n is said to be a cone if and only if x ∈ C ⇒ x/ε ∈ C for all ε > 0.
Definition 2 (Normal Cone) Let a closed non-empty convex set X ∈ R n and a symmetric and positive-definite matrix U ∈ R n×n be given. The normal cone to X at the point x ∈ R n with the metric U is defined as follows:

N U X (x) ∆ = ȷ {ξ ∈ R n | ξ T U (x * -x) ≤ 0 ∀x * ∈ X } if x ∈ X ∅ if x ̸ ∈ X . ( 1
)
Definition 3 (Projection) Let a closed non-empty convex set X ∈ R n and a symmetric and positive-definite matrix U ∈ R n×n be given. The projection of x ∈ R n on the set X with the metric U is defined as follows:

proj U X (x) ∆ = argmin ξ∈X (x -ξ) T U (x -ξ). ( 2
)
For notational brevity, we write

N In X (x) = N X (x) and proj In X (x) = proj X (x). It should be noted that, if x ∈ IntX , N U X (x) = {0} and proj U X (x) =
x are satisfied. The following relations are useful for the derivation in this paper:

N U X (u) = U -1 N X (u) (3) x -u ∈ -N U X (x) ⇐⇒ x = proj U X (u). (4)
The relation (3) can be obtained through a straightforward derivation from the definition (1) of N U X (u). The relation ( 4) is explained in Section A.3 of [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF].

x y N X (r b ) N X (r a ) N X (r c ) (a) (b) X r a r b r c X r a r b r c
Fig. 1 Normal cones (a) with the Euclidean metric and (b) with the singular metric defined by (9).

Limit of a Family of Sets

Let us recall that a limit of a function is formally defined as follows:

lim ε→+0 f (ε) ∆ = f 0 s.t. `∀ϵ 1 > 0, ∃δ 1 > 0 s.t. 0 < ε < δ 1 ⇒ |f (ε) -f 0 | < ϵ 1 ´, (5)
which is usually referred to as the ε-δ definition. In the same way, a limit of a family of sets Xε ∈ R n parameterized by a positive scalar ϵ can be defined as follows:

lim ε→+0 Xε ∆ = ˘ξ ∈ R n ˛∀ϵ 1 > 0, ∃δ 1 > 0 s.t. 0 < ε < δ 1 ⇒ `∃ξ 1 ∈ Xε s.t. |ξ -ξ 1 | < ϵ 1 ´¯. (6)
The following example illustrates this definition:

lim ε→+0 ȷ » y/ε y - ∈ R 2 ˛y ≥ 0 ff = ȷ » x 0 - ∈ R 2 ˛x ≥ 0 ff . ( 7
)
Here, note that the y component (the second entry) of the vector in the left-hand side never reaches zero with ε → 0, but the limit of the left-hand side is the set with y = 0. In this sense, this example appears paradoxical but it matches the definition (6).

Normal Cones with a Singular Metric

This paper uses limits of families of normal cones in parameterized metrics of the form

lim ε→+0 N Uε X (r) (8)
where X ⊂ R n , r ∈ R n and Uε ∈ R n×n is a symmetric and positive-definite matrix for all ε > 0. When Uε includes the division by ε, it can be referred to as a singular metric.

As an illustrative example, let us consider a 2-dimensional case with

Uε ∆ = diag[1, 1/ε]. (9)
In this definition of the metric Uε, the normal cones of the form (8) can be said to be singular-metric normal cones. Fig. 1 illustrates some instances of Euclideanmetric and singular-metric normal cones. Now, let us focus on N X (ra) in Fig. 1(a). It can be written as follows:

N X (ra) = ȷ » x αx - ∈ R 2 ˛α ∈ [α 0 , α 1 ], x ≤ 0 ff ( 10
)
where α 0 and α 1 are scalars satisfying α 0 < α 1 < 0. By using the rule (3), one can easily obtain

lim ε→+0 N Uε X (ra) = lim ε→+0 ȷ » x εαx - ∈ R 2 ˛α ∈ [α 0 , α 1 ], x ≤ 0 ff = ȷ » x 0 - ∈ R 2 ˛x ≤ 0 ff , ( 11
)
which is illustrated in Fig. 1(b). The cones at the point r b can also be obtained by modifying (10)(11) by x ≥ 0 and setting α 0 = α 1 > 0.

As another example, let us consider N X (rc) in Fig. 1(a), which can be written as

N X (rc) = ȷ » αy y - ∈ R 2 ˛α ∈ [α 0 , α 1 ], y ≤ 0 ff (12)
where α 0 < 0 < α 1 . This leads to the following:

lim ε→+0 N Uε X (rc) = lim ε→+0 ȷ » αỹ εỹ - ∈ R 2 ˛α ∈ [α 0 , α 1 ], ỹ ≤ 0 ff = ȷ » x y - ∈ R 2 ˛x ∈ R, y ≤ 0 ff . ( 13
)
The derivation in (13) can be explained by the fact that, for any x ∈ R, y ≤ 0, α 0 < 0 and α 1 > 0, one can choose an arbitrarily small ε > 0 that satisfies εx/y ∈ [α 0 , α 1 ], and thus can also choose ỹ = y/ε and α = εx/y, which make [αỹ, εỹ] T belong to N Uε X (rc). In conclusions, one can see that the singular-metric normal cones in Fig 1(b) can be obtained by infinitely elongating the Euclidean-metric normal cones in Fig 1(a) in the x direction. This paper uses this type of singular-metric normal cones to describe the relation between the force and kinematic variables such as velocity.

Overview on Equality Constraints

This section overviews DAEs that represent mechanical systems involving positionlevel and velocity-level constraints. Such a DAE is represented as a pair of a constraint equation and a differential equation with a Lagrange multiplier, which represents the constraint force. The Baumgarte stabilization is also overviewed as a method to convert DAEs to ODEs. Its similarity to implicit discretization and its difference from the penalty-based relaxation are also discussed here.

Equality Constraints and Lagrange Multiplier

Let us consider the class of n-dimensional constrained mechanical systems described as follows:

M (q)q + Φ(q, q, t) = J(q) T λ ( 14
) » Jv(q) q Ψp(q) - = 0 ( 15
)
where

J(q) ∆ = » Jv(q) Jp(q) - ∈ R m×n (16) Jp(q) ∆ = ∂Ψp(q)/∂q ∈ R mp×n ( 17
)
and mp ≤ m ≤ n. Here, M (q) ∈ R n×n is a symmetric and positive-definite matrix that represents the inertia of the system, q ∈ R n is the generalized coordinate of the system, and Φ : R n × R n × R → R n represents the sum of nonlinear inertial forces and external forces. Equation ( 15) represents a constraint imposed on q and q, and λ ∈ R m is the Lagrange multiplier, which represents the force caused by the constraint. The constraint Ψp(q) = 0 can be said to be a position-level (or holonomic) constraint. The constraint Jv(q) q = 0 is a non-holonomic constraint that can be said to be a velocity-level constraint, which appears in, for example, slipless wheeled systems [START_REF] Kim | Dynamic analysis of a nonholonomic two-wheeled inverted pendulum robot[END_REF] and other classes of vehicles [START_REF] Kobilarov | Lie group integrators for animation and control of vehicles[END_REF].

One way to numerically solve the problem (14)( 15) is to take the derivatives of (15) to isolate q to construct a set of simultaneous equations with [START_REF] Bertails-Descoubes | A nonsmooth Newton solver for capturing exact Coulomb friction in fiber assemblies[END_REF]. Specifically, by taking the first-and second-derivatives of the first and the second lines of [START_REF] Bonnefon | Quartic formulation of Coulomb 3D frictional contact[END_REF], respectively, one obtains the following constraint: (Jv(q) q)

d 2 Ψp(q) dt 2 3 7 5 = 0, (18)
or equivalently, J(q)q + H(q, q) q = 0 (

where

H(q, q) ∆ = » Hv(q, q) Hp(q, q) - ∆ = J(q). ( 20
)
Eliminating q from ( 14) and [START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF] provides the following:

J(q)M (q) -1 (J(q) T λ -Φ(q, q, t)) + H(q, q) q = 0, [START_REF] De Saxcé | The bipotential method: A constructive approach to design the complete contact law with friction and improved numerical algorithms[END_REF] and when the Delassus operator

A(q) ∆ = J(q)M (q) -1 J(q) T ( 22
)
is invertible, it can be solved with respect to λ as follows:

λ = A(q) -1 " J(q)M (q) -1 Φ(q, q, t) -H(q, q) q" . ( 23
)
By substituting [START_REF] García De Jalón | Multibody dynamics with redundant constraints and singular mass matrix: existence, uniqueness, and determination of solutions for accelerations and constraint forces[END_REF] in [START_REF] Bertails-Descoubes | A nonsmooth Newton solver for capturing exact Coulomb friction in fiber assemblies[END_REF] to eliminate λ, one can view the DAE (14)(15) as an ODE, which can be integrated easily, with the invariant set {(q, q) | Jv(q) q = 0, Ψp(q) = 0}.

Baumgarte Stabilization

One problem caused by the use of the differentiated constraint ( 19) is that it causes numerical drift from the original constraint [START_REF] Bonnefon | Quartic formulation of Coulomb 3D frictional contact[END_REF]. The Baumgarte stabilization method [START_REF] Baumgarte | Stabilization of constraints and integrals of motion in dynamical systems[END_REF] is developed as a computational technique for avoiding this artifact. It is fully developed in the continuous-time domain, being in contrast to discrete-time integration schemes [START_REF] Arnold | Convergence of the generalized-α scheme for constrained mechanical systems[END_REF][START_REF] Lunk | Solving constrained mechanical systems by the family of Newmark and α-methods[END_REF] elaborated to avoid drifting. Although some practical weaknesses have been pointed out [START_REF] Ascher | Stabilization of constrained mechanical systems with DAEs and invariant manifolds[END_REF][START_REF] Flores | A parametric study on the Baumgarte stabilization method for forward dynamics of constrained multibody systems[END_REF], it has been recognized as a useful tool for some applications [START_REF] Joli | A modular modeling approach to simulate interactively multibody systems with a Baumgarte/Uzawa formulation[END_REF][START_REF] Marquis-Favre | A planar mechanical library in the AMESim simulation software. part I: Formulation of dynamics equations[END_REF]. In addition, its continuous-time representation may give some insights for the analytical purposes, being independent of discretization schemes.

In this method, the constraint ( 15) is replaced by

2 6 4 d dt (Jv(q) q) d 2 Ψp(q) dt 2 3 7 5 + 2ζ δ 2 4 Jv(q) q dΨp(q) dt 3 5 + 1 δ 2 » 0 Ψp(q) - = 0, (24)
or equivalently,

J(q)q + H(q, q) q + 2ζ δ J(q) q + 1 δ 2 » 0 Ψp(q) - = 0 (25)
where ζ and δ are appropriately chosen positive scalars. The initial condition must be Ψp(q) = 0 and J(q) q = 0 at t = 0. One can say that the modified constraint ( 25) is analytically equivalent to the original constraint [START_REF] Bonnefon | Quartic formulation of Coulomb 3D frictional contact[END_REF] or the differentiated constraint [START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF] since, for any functions x of time, x(t) = 0 ∀t ≥ 0 is analytically equivalent to ẋ(t) = 0 ∀t ≥ 0 with the initial condition x(t) = 0 at t = 0. Therefore, one can say that the modification from [START_REF] Bonnefon | Quartic formulation of Coulomb 3D frictional contact[END_REF] or [START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF] to [START_REF] Glocker | Complementarity problems in multibody systems with planar friction[END_REF] is not an approximation in the analytical sense. In numerical computation, however, the last two terms of the left-hand side of (25) may deviate from zero, and thus this modification can be said to be an approximation.

Based on the constraint (25), λ is obtained as follows:

λ = A(q) -1 b(δ, q, q, t) (26)
where b(δ, q, q, t)

∆ = J(q)M (q) -1 Φ(q, q, t) -H(q, q) q - 2ζ δ J(q) q - 1 δ 2 » 0 Ψp(q) - . ( 27
)
By substituting ( 26) into ( 14) and using the initial values of Ψp(q) = 0 and J(q) q = 0, one can integrate the system (14)(25) as an ODE, of which the numerical solution with respect to q satisfies (14) [START_REF] Bonnefon | Quartic formulation of Coulomb 3D frictional contact[END_REF] approximately.

It may need to be pointed out that there is no physical (more precisely, rheological) interpretation for [START_REF] Glocker | Complementarity problems in multibody systems with planar friction[END_REF]. The position-level constraint in [START_REF] Bonnefon | Quartic formulation of Coulomb 3D frictional contact[END_REF] is replaced by a second-order differential equation that may be seen analogous to a massspring-dashpot system. The system is however an autonomous system, which is not influenced by any other factors, although the evolution of q is influenced by the force term Φ(q, q). Considering this, it is not obvious to construct a consistent interpretation in which the mass M (q) is connected to an imaginary spring-dashpot system. This fact differentiates the Baumgarte stabilization from the penalty-based method, which will be discussed in section 3.4.

Relations to Implicit Integration

Another approach for the time integration of the system (14)( 15) is to use implicit integration schemes. We here show an analogic relation between the implicit Euler integration and the Baumgarte stabilization. In a fully-implicit Euler integration of ( 14) [START_REF] Bonnefon | Quartic formulation of Coulomb 3D frictional contact[END_REF], one should solve the following nonlinear equations with respect to q k+1 , u k+1 and λ k+1 at every timestep:

M (q k+1 ) u k+1 -u k h + Φ(q k+1 , u k+1 , (k + 1)h) = J(q k+1 ) T λ k+1 ( 28
) » Jv(q k+1 )u k+1 Ψp(q k+1 ) - = 0 ( 29
)
q k+1 = q k + hu k+1 . ( 30
)
Here, u ∆ = q, the subscripts are the discrete-time indices, and h > 0 is the timestep size. Assuming that h is small enough, we can consider the following approximations:

M (q k+1 ) ≈ M (q k ) (31)
J(q k+1 ) T λ k+1 ≈ J(q k ) T λ k+1 (32) Φ(q k+1 , u k+1 , (k + 1)h) ≈ Φ(q k , u k , kh) ( 33
)
Jv(q k+1 )u k+1 ≈ Jv(q k )u k+1 + hHv(q k , u k )u k (34) Ψp(q k+1 ) ≈ Ψp(q k ) + hJp(q k )u k+1 + h 2 Hp(q k , u k )u k . ( 35
)
By using these approximations, the nonlinear problem (28)(29)(30) can be linearly approximated and λ k+1 can be analytically obtained as follows:

λ k+1 = A(q k ) -1 " J(q k )M (q k ) -1 Φ(q k , u k , kh) -H(q k , u k )u k - 1 h J(q k )u k - 1 h 2 » 0 Ψp(q k ) -« . ( 36
)
This is equivalent to

λ k+1 = A(q k ) -1 b(h, q k , u k , kh) with ζ = 1/2 according to the definition (27) of b : R × R n × R n × R → R m .
Notice the similarity between (26)( 27) and ( 36).

There could be some variations in the discretization scheme. For example, the terms involving H(q k , u k ) may be omitted from ( 34) and ( 35) if the temporal change in J(q) is slow enough. Moreover, if the temporal change in Φ(q, q, t) is not slow enough, the Jacobians of Φ(q, q, t) may need to be included in [START_REF] Klarbring | A mathematical programming approach to three-dimensional contact problems with friction[END_REF]. Although λ k+1 may not be exactly the same form as that of [START_REF] Jean | The non-smooth contact dynamics method[END_REF], one can see that λ k+1 with h → +0 can be close to that of [START_REF] Jean | The non-smooth contact dynamics method[END_REF] with δ → +0. This fact may imply the existence of some underlying relations between the Baumgarte stabilization and the implicit Euler discretization, which may need further theoretical investigation. More practically, it suggests that some numerical methods in the previous studies based on implicit discretization can be applied to the numerical integration of Baumgarte-stabilized systems [START_REF] Bertails-Descoubes | A nonsmooth Newton solver for capturing exact Coulomb friction in fiber assemblies[END_REF] [START_REF] Glocker | Complementarity problems in multibody systems with planar friction[END_REF].

Relation to Penalty-based Relaxation

It should be noted that the Baumgarte stabilization is intrinsically different from the penalty-based relaxation. When the penalty-based method is applied to deal with the constraint [START_REF] Bonnefon | Quartic formulation of Coulomb 3D frictional contact[END_REF], a force λ is usually determined as a spring-dashpot force penalizing the deviation from the constraint [START_REF] Bonnefon | Quartic formulation of Coulomb 3D frictional contact[END_REF]. Although there are several variations [START_REF] Bauchau | Review of contemporary approaches for constraint enforcement in multibody systems[END_REF] to realize this, the following scheme can be thought of as an example:

λ = K "» ev ep - + β » ėv ėp -« (37)
»

Jv(q) q + ėv Ψp(q) + ep - = 0. ( 38
)
Here, K is a symmetric and positive-definite (and often diagonal) matrix representing the stiffness of the spring, and Kβ can be seen as a damping coefficient matrix. The vector [e T v , e T p ] T represents the deviation from the constraint [START_REF] Bonnefon | Quartic formulation of Coulomb 3D frictional contact[END_REF]. An important difference between (37)( 38) and ( 25) is that (37)( 38) is different from the original constraint (15) even in the continuous-time domain. In the Baumgarte stabilization (25), the ideal constraint ( 15) can be satisfied with t → ∞, but in the penalty-based relaxation (37) [START_REF] Lunk | Solving constrained mechanical systems by the family of Newmark and α-methods[END_REF], it can never be satisfied unless Φ(q, q, t) ≡ 0. The deviation from the constraint may be reduced by setting higher values to the entries of K, but it results in the high stiffness of the equations, threatening the numerical stability especially with explicit integration schemes. An advantage of the penalty-based relaxation (37) [START_REF] Lunk | Solving constrained mechanical systems by the family of Newmark and α-methods[END_REF] is that the force λ is independent of the inertia matrix M (q), and thus it is easier in computation. Thus, one can say that the penalty-based relaxation is a more relaxed method than the Baumgarte stabilization to realize easier computation by sacrificing the accuracy.

New Approach for Inequality Constraints

This section considers the system [START_REF] Bertails-Descoubes | A nonsmooth Newton solver for capturing exact Coulomb friction in fiber assemblies[END_REF] with the force λ being restricted to a closed set C ⊂ R m . In such a case, ( 15) is satisfied when λ ∈ IntC, but can be violated when λ ∈ ∂C. Mechanical systems with unilateral constraints and/or three-dimensional Coulomb friction are included in this class of systems. This section provides a new DAI representation of this class of systems with the use of a singular-metric normal cone. Moreover, its ODE approximation inspired by the Baumgarte stabilization is also presented.

New Representation of Inequality Constraints

This paper presents the following form as a general representation of rigid-body systems with friction and unilateral constraints:

M (q)q + Φ(q, q, t) = J(q) T λ ( 39
) » Jv(q) q Ψp(q) - ∈ -lim ε→+0 N Uε C (λ). ( 40
)
The matrix Uε is a symmetric and positive-definite matrix for all ε > 0, which defines the system behavior under the condition λ ∈ ∂C. It should be noted that the problem (39)(40) includes the problem (14)(15) as a special case of C = R m .

Here we show some examples of systems represented by the form (39) [START_REF] Mazhar | Chrono: a parallel multi-physics library for rigid-body, flexible-body, and fluid dynamics[END_REF]. First, let us consider m frictionless unilateral constraints:

Ψp(q) ≥ 0 (41)
where the inequality is component-wise. In such a system, the contact force λ satisfies 0 ≤ λ ⊥ Ψp(q) ≥ 0, and thus the system can be written as follows:

M (q)q + Φ(q, q, t) = (∂Ψp(q)/∂q) T λ (42)

Ψp(q) ∈ -N R m + (λ). ( 43
)
It should be noted that ( 43) could be equivalently rewritten as λ ∈ -N R m + (Ψp(q)). Nevertheless, the representation (42)( 43) is convenient to show that it is a special case of (39) [START_REF] Mazhar | Chrono: a parallel multi-physics library for rigid-body, flexible-body, and fluid dynamics[END_REF] with m = mp, Uε = Im, J(q) = Jp(q) = ∂Ψp(q)/∂q, and C = R m + . As another example, let us consider a point mass M > 0 subject to an external force fe ∈ R 2 and a Coulomb friction force λ ∈ R 2 . The system can be represented as follows:

M q = fe + λ (44) λ ∈ -F sgn( q) (45)
where sgn(x) is a two-dimensional set-valued signum function sgn(x)

∆ = ȷ x/∥x∥ if x ̸ = 0 B 2 otherwise, (46)
B 2 denotes the 2-dimensional unit ball and F > 0. (Here we do not consider the effect of normal force or the maximum static friction force being larger than the kinetic friction force.) This system can be rewritten as follows:

M q -fe = λ (47) q ∈ -N F B2 (λ), (48)
which is a special case of (39)(40) with m = 2, mp = 0, J(q) = Jv(q) = Uε = I 2 , and

C = F B 2 .
The main advantage of the new representation (40) appears in cases with three-dimensional frictional unilateral constraints. Let us denote the position of the contact point by [p T t , pn] T = Ψ (q) ∈ R 3 where Ψ is a map from R n to R 3 . Let us assume that the unilateral constraint pn ≥ 0 (pn ∈ R) is imposed by a contact surface, and that ṗt ∈ R 2 is the tangential velocity. Then, the contact force λ = [λ T t , λn] T ∈ R 3 can be related to pn and ṗt as follows:

λ t ∈ -µλnsgn( ṗt ) ( 49
) λn ∈ ȷ 0 if pn > 0 R + if pn = 0, (50)
or equivalently,

ṗt ∈ -N µλnB2 (λ t ) ( 51
) pn ∈ -N R + (λn). ( 52
)
Carefully observing the relation between the vector λ and the vector [ ṗt , pn] T , one can see that the relation (51)(52) can be rewritten as follows:

» ṗt pn - ∈ -lim ε→+0 N Uε C (λ) (53)
where

C ∆ = n [λ T t , λn] T ∈ R 3 ˛∥λ t ∥ ≤ µλn o ( 54
)
Uε ∆ = diag[1, 1, 1/ε]. ( 55
)
The singular-metric normal cone in the right-hand side of ( 53) is illustrated in Fig. 2. Note that ( 53) is a special case of ( 40) with m = 3, mp = 1, ( 54), (55) and

Jv(q) = » 1 0 0 0 1 0 - ∂Ψ (q) ∂q ( 56
)
Ψp(q) = ˆ0 0 1 ˜Ψ (q). ( 57
)
The same idea can be applied to more complicated systems with many frictional contacts. Let us consider an n-dimensional system with mp potential threedimensional frictional unilateral constraints. Let us denote the mp unilateral constraints by Ψp(q) ≥ 0 with Ψp : R n → R mp and the tangential velocities at the contact points by Jv(q) q ∈ R 2mp with Jv : R n → R 2mp . Let us also denote the contact forces by λ ∈ R 3mp , of which λ {2mp+i} and (λ {2i-1} , λ {2i} ) are the normal and tangential components, respectively, of the i-th contact force, where the subscript {i} denotes the i-th component of the vector. Then, such a system can be represented by [START_REF] Mazhar | Chrono: a parallel multi-physics library for rigid-body, flexible-body, and fluid dynamics[END_REF] specialized by

C ∆ = n λ ∈ R 3mp ˛q(λ {2i-1} ) 2 + (λ {2i} ) 2 ≤ µ i λ {2mp+i} , ∀i ∈ {1, • • • , mp} o (58) Uε ∆ = » I 2mp 0 0 Im p /ε - (59)
where µ i is the friction coefficient of the i-th potential contact.

The use of normal cones has been known to be useful for representing nonsmooth dynamical systems [START_REF] Anitescu | An iterative approach for cone complementarity problems for nonsmooth dynamics[END_REF][START_REF] Kanno | Three-dimensional quasi-static frictional contact by using second-order cone linear complementarity problem[END_REF][START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Studer | Representation of normal cone inclusion problems in dynamics via non-linear equations[END_REF], but the form [START_REF] Mazhar | Chrono: a parallel multi-physics library for rigid-body, flexible-body, and fluid dynamics[END_REF], involving a singularmetric normal cone for the representation of frictional unilateral constraints, is not found in the previous studies. One important point here is to use the normal cone defined in the constraint force space R m (the space of the Lagrange multipliers λ) to preserve the structure of the left-hand side of the conventional form [START_REF] Bonnefon | Quartic formulation of Coulomb 3D frictional contact[END_REF] for the equality constraint. By using the representation [START_REF] Mazhar | Chrono: a parallel multi-physics library for rigid-body, flexible-body, and fluid dynamics[END_REF], inequality constraints caused by frictional unilateral constraints can be seen as a straightforward generalization of equality constraints of the form (15).

Baumgarte-Like Relaxation

In order to solve the DAI (39) [START_REF] Mazhar | Chrono: a parallel multi-physics library for rigid-body, flexible-body, and fluid dynamics[END_REF], one needs to isolate q from (40), as can be done with the DAE (14) [START_REF] Bonnefon | Quartic formulation of Coulomb 3D frictional contact[END_REF]. Recall that the equality constraint (15) can be replaced by its differentiated constraint [START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF]. This approach, however, cannot be applied to Fig. 2 The friction cone C and its singular-metric normal cones with Uε defined by (55). Note that the normal cone at λ = 0 is the half-space z ≤ 0.

the inclusion constraint (40) because, with a function x of time and a set X , x ∈ X has in general nothing to do with ẋ ∈ X , while x ≡ 0 implies ẋ = 0.

A key observation here is that, with a cone C and a twice-differentiable function x of time t, if x satisfies x ∈ C for all t > 0,

lim δ→+0 " ẍ + 2ζ ẋ δ + x δ 2 « ∈ C ( 60
)
is satisfied for any ζ > 0. Note that here we exploit the property of cones indicated by Definition 1. This can be viewed as a generalization of the lexicographic nonnegativity [START_REF] Van Der Schaft | Complementarity modeling of hybrid systems[END_REF] [16, p.214], which is written as

(x, ẋ, ẍ) ≽ 0 (61)
and can also be written as

lim δ→+0 " ẍ + 2ζ ẋ δ + x δ 2 « ≥ 0. ( 62
)
These expressions mean that all of {x, ẋ, ẍ} are zero or the first nonzero entry of {x, ẋ, ẍ} is positive. The expression (62) is a special case of (60) with C = R + .

In the same light as (60), the following representation can be seen as an approximation of (40): (Jv(q) q)

d 2 Ψp(q) dt 2 3 7 5 + 2ζ δ 2 4 Jv(q) q dΨp(q) dt 3 5 + 1 δ 2 » 0 Ψp(q) - ∈ -lim ε→+0 N Uε C (λ), (63)
or equivalently, J(q)q + J(q) q + 2ζ δ J(q) q + 1

δ 2 » 0 Ψp(q) - ∈ -lim ε→+0 N Uε C (λ). ( 64
)
This representation can be seen as a relaxation of the original constraint [START_REF] Mazhar | Chrono: a parallel multi-physics library for rigid-body, flexible-body, and fluid dynamics[END_REF], and the left-hand side has exactly the same form as the Baumgarte stabilization [START_REF] Glocker | Complementarity problems in multibody systems with planar friction[END_REF]. Thus, we hereafter refer to it as a Baumgarte-like relaxation. By eliminating q from (39) and (64), one can obtain the following expression:

A(q)λ -b(δ, q, q, t) ∈ -lim ε→+0 N Uε C (λ) (65)
where A(q) and b(δ, q, q, t) are defined in ( 22) and [START_REF] Jean | The non-smooth contact dynamics method[END_REF], respectively. This expression includes (26) as a special case of C = R m . By substituting the solution λ of (65) into (39) and using the initial values of Ψp(q) = 0 and J(q) q = 0, one can integrate the relaxed DAI (39)(64) to obtain q that satisfies the original DAI (39) [START_REF] Mazhar | Chrono: a parallel multi-physics library for rigid-body, flexible-body, and fluid dynamics[END_REF] approximately. In other words, under the assumptions that the inclusion (65) is uniquely solvable with respect to λ and that the solution λ is a Lipschitz function of q and q, the relaxed DAI (39)(64) can be considered equivalent to an ODE. The solution of (65), however, is not always obvious to obtain. From ( 3) and (4), if A(q)Uε is a symmetric and positive-definite matrix (e.g., if either or both of Uε and A(q) are scalar multiples of the identity matrix), the closed-form solution of (65) can be obtained as follows:

λ = lim ε→+0 proj A(q)Uε C (A(q) -1 b(δ, q, q, t)).
(66)

Otherwise, the operator "proj" cannot be used and the existence and the uniqueness of the solution of the algebraic inclusion (65) are not guaranteed. As will be shown in section 5.3, the problems known as Painlevé paradoxes [START_REF] Génot | New results on Painlevé paradoxes[END_REF] are raised when (65) has non-unique solutions. Some previous studies [START_REF] Acary | A formulation of the linear discrete Coulomb friction problem via convex optimization[END_REF][START_REF] Alart | A mixed formulation for frictional contact problems prone to newton like solution methods[END_REF][START_REF] Bertails-Descoubes | A nonsmooth Newton solver for capturing exact Coulomb friction in fiber assemblies[END_REF][START_REF] Bonnefon | Quartic formulation of Coulomb 3D frictional contact[END_REF][START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF] provide numerical methods to solve problems equivalent to (65) as will be discussed in section 4.5.

One should note that the relaxed constraint (64) is analytically different from the original constraint [START_REF] Mazhar | Chrono: a parallel multi-physics library for rigid-body, flexible-body, and fluid dynamics[END_REF]. This fact is in contrast to the fact that the original Baumgarte stabilization for equality constraints does not alter the constraint analytically. It should also be noted that, as will be shown in section 5.2, the free parameter ζ influences the behavior of the solution both numerically and analytically.

Relations to Previous Representations

To represent the system with inequality constraints, some previous studies employ different representations at different kinematic levels, i.e., acceleration, velocity, and position levels [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF]Sec.8.3]. For example, let us assume that q ∈ R is an absolutely continuous function of time t and that there exists a position-level inequality constraint q ≥ 0. Then, one should consider a velocity-level representation q ≥ 0 of the constraint, which is valid only when q = 0. In the same light, there is also an acceleration-level representation q ≥ 0, which is valid only if q = 0 and q = 0.

The expression (64) includes all conditions in a unified form, and can be reduced into each kinematic level by exploiting the property of cones indicated by Definition 1 and setting δ → +0. For example, when Ψp(q) > 0, setting δ → +0 in (64) results in

» 0 Ψp(q) - ∈ -lim ε→+0 N Uε C (λ) (67)
because the left-hand side can be multiplied by δ 2 due to Definition 1 of cones. This representation (67) can be seen as a position-level representation. In the same way, when Ψp(q) = 0 and when Ψp(q) = 0 and J(q) q = 0, respectively, setting δ → +0 with (64) yields the following:

J(q) q ∈ -lim ε→+0 N Uε C (λ) (68) J(q)q + J(q) q ∈ -lim ε→+0 N Uε C (λ). ( 69
)
These representations can be named velocity-level and acceleration-level representations, respectively. The acceleration-level representation is especially important for event-driven schemes [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF]Chap. 8] [26] to detect discontinuous events, as have been exemplified in [START_REF] Glocker | Complementarity problems in multibody systems with planar friction[END_REF][START_REF] Leine | Nonlinear dynamics and modeling of various wooden toys with impact and friction[END_REF][START_REF] Pfeiffer | Unilateral problems of dynamics[END_REF][START_REF] Studer | Representation of normal cone inclusion problems in dynamics via non-linear equations[END_REF].

It should be noted that the presented representation does not explicitly involve impulse because, due to the relaxation using the non-zero value δ, the acceleration is always bounded and discontinuous velocity jumps do not happen. It poses a difficulty to relate the presented representation to measure-differential inclusions [START_REF] Acary | Higher order Moreau's sweeping process: mathematical formulation and numerical simulation[END_REF][START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: Definition and outlook[END_REF]. It is also a non-trivial issue to find what kind of impact laws can be obtained as the result of this scheme of relaxation. Its relation to the simplest one-dimensional impact law will be discussed in section 5.2.

Relations to Implicit Integration

In the same way as the discussion in section 3.3, we can see that the presented scheme has analogical relation to the implicit integration scheme. By using the approximations [START_REF] Kikuuwe | Admittance and impedance representations of friction based on implicit Euler integration[END_REF] to [START_REF] Klarbring | A mathematical programming approach to contact problems with friction and varying contact surface[END_REF] and setting ζ = 1/2, if the timestep size h is small enough, the original DAI (39)(40) can be approximated by the following discretetime representation:

A(q k )λ k+1 -b(h, q k , u k , kh) ∈ -lim ε→+0 N Uε C (λ k+1 ). ( 70
)
As mentioned in the last paragraph of section 3.3, there may be some variations in b(h, q k , u k , kh) according to possible choice of approximations, but the form of the algebraic problem (70), which is exactly of the same as (65), is preserved.

Numerical Methods to Solve (65)

The problem (65) is generally not easy to solve, but many previous studies have been devoted to algebraic problems of this form in different notations. One of major approaches is to facetize the friction cones that comprise the cone C, that is, to approximate the friction cones by polygonal cones. The purpose of this approximation is to reformulate the problem as a Linear Complementarity Problem (LCP) for the ease of application of established numerical solvers. This approach has been introduced for quasi-static problems with elastic continua [START_REF] Klarbring | A mathematical programming approach to three-dimensional contact problems with friction[END_REF][START_REF] Klarbring | A mathematical programming approach to contact problems with friction and varying contact surface[END_REF], and have been employed in Stewart and Trinkle's method [START_REF] Stewart | An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and Coulomb friction[END_REF] and Anitescu and Potra's method [START_REF] Anitescu | Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems[END_REF] for simulating dynamical systems. Another class of approaches is to reformulate (65) into a cone complementarity problem [START_REF] Anitescu | An iterative approach for cone complementarity problems for nonsmooth dynamics[END_REF], which is of the form

C * g ∋ Agxg -bg ⊥ xg ∈ Cg, (71)
or equivalently,

Agxg -bg ∈ -N Cg (xg) ( 72
)
where xg is an unknown vector, Ag is a given symmetric and positive-definite matrix, bg is a given vector, Cg is a given cone, and C * g is the dual cone of Cg. The problem (65) does not fit to this form in general. In some previous studies [START_REF] Acary | A formulation of the linear discrete Coulomb friction problem via convex optimization[END_REF][START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF][START_REF] De Saxcé | New inequality and functional for contact with friction: The implicit standard material approach[END_REF][START_REF] De Saxcé | The bipotential method: A constructive approach to design the complete contact law with friction and improved numerical algorithms[END_REF], however, the left-hand side of (65) is appropriately modified so that it fits to the form (71) or (72). Such an approach was initially introduced by De Saxcé as detailed in [3, Sec.3.9.2].

Here we discuss the relation between the proposed approach and De Saxcé's approach. Let us again consider a system including mp potential three-dimensional frictional unilateral constraints, which is described by [START_REF] Marquis-Favre | A planar mechanical library in the AMESim simulation software. part I: Formulation of dynamics equations[END_REF] [START_REF] Mazhar | Chrono: a parallel multi-physics library for rigid-body, flexible-body, and fluid dynamics[END_REF] specialized by (58)(59). Based on the proposed relaxation scheme, one arrives in the problem (65) to obtain λ. De Saxcé's approach is to rewrite (65) of this special case into the following form:

A(q)λ -b(δ, q, q, t) + Es(η(δ, q, q, λ, t)) ∈ -N C (λ), ( 73
)
or equivalently,

C * ∋ A(q)λ -b(δ, q, q, t) + Es(η(δ, q, q, λ, t)) ⊥ λ ∈ C ( 74
)
where η(δ, q, q, λ, t)

∆ = A(q)λ -b(δ, q, q, t) ∈ R 3mp (75) s(η) ∆ = 2 6 6 6 4 q (η {1} ) 2 + (η {2} ) 2 . . . q (η {2mp-1} ) 2 + (η {2mp} ) 2 3 7 7 7 5 ∈ R mp (76) E ∆ = » 0 2mp×mp diag[µ 1 , • • • , µm p ] - ∈ R 3mp×mp .
(77) Fig. 3 illustrates (73) of the case of mp = 1. Here, the left-hand side of ( 73) is orthogonal to the reaction force λ. Some methods have been presented [START_REF] Acary | A formulation of the linear discrete Coulomb friction problem via convex optimization[END_REF][START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF] to solve the simultaneous equations of the form (73)(75)(76) with respect to η, s, and λ. Kanno et al. [29] also treated a problem equivalent to (65) by reformulating it into another form of cone complementarity problem by introducing an auxiliary variable.

Anitescu and Tasora [START_REF] Anitescu | An iterative approach for cone complementarity problems for nonsmooth dynamics[END_REF] proposed a method in which the term +Es in ( 73) is simply neglected1 to fit the problem in the standard form of a cone complementarity problem. This modification corresponds to neglecting Uε in (65), resulting in the solution of the following form:

λ = proj A(q) C (A(q) -1 b(δ, q, q, t)), (78)
or equivalently,

A(q)λ -b(δ, q, q, t) ∈ -N C (λ). (79)
Their method has been employed in an open-source software package [START_REF] Mazhar | Chrono: a parallel multi-physics library for rigid-body, flexible-body, and fluid dynamics[END_REF]. As has been recognized, and as is obvious from its physical interpretation presented in [7, Fig. 1], their method causes small motion in the normal direction even during pure tangential slipping. In contrast to this is Nakaoka et al.'s [START_REF] Nakaoka | Constraint-based dynamics simulator for humanoid robots with shock absorbing mechanisms[END_REF] approach, in which A(q) as a factor of norm metric in (65) is neglected. That is, their method can be described as follows:

λ = lim ε→+0 proj Uε C (A(q) -1 b(δ, q, q, t)), (80)
or equivalently,

λ -A(q) -1 b(δ, q, q, t) ∈ -lim ε→+0 N Uε C (λ). (81)
They applied a version of projected Gauss-Seidel method to obtain λ of (80). The method presented by Silcowitz et al. [START_REF] Silcowitz | Interactive rigid body dynamics using a projected Gauss-Seidel subspace minimization method[END_REF] is similar to (80), but they decouple each contact into two one-dimensional friction problems. This modification is equivalent to approximating the friction cones, which comprise C, by square pyramids.

The problem (65) can also be rewritten as an equality of the form Γ AC (A(q)λ -b(δ, q, q, t), λ) = 0 (82) by using an appropriate nonlinear function Γ AC : R m × R m → R m , as proposed by Alart and Curnier [START_REF] Alart | A mixed formulation for frictional contact problems prone to newton like solution methods[END_REF] and detailed in [START_REF] Bertails-Descoubes | A nonsmooth Newton solver for capturing exact Coulomb friction in fiber assemblies[END_REF]. For the case with a single frictional unilateral constraint (i.e., m = 3 and mp = 1), it has been shown that the problem (65), or equivalently (73), is solvable based on particular quartic polynomials [START_REF] Bonnefon | Quartic formulation of Coulomb 3D frictional contact[END_REF] [19, Appendix B].

Relations to Penalty-based Relaxation

Here we show the difference between the presented relaxation scheme and the penalty-based relaxation scheme. Let us again consider the original constraint [START_REF] Mazhar | Chrono: a parallel multi-physics library for rigid-body, flexible-body, and fluid dynamics[END_REF] imposed on the dynamical system [START_REF] Marquis-Favre | A planar mechanical library in the AMESim simulation software. part I: Formulation of dynamics equations[END_REF]. Considering an imaginary spring-dashpot system of which the force is λ and the displacement is [e T v , e T p ] T , and considering that [e T v , e T p ] T can be determined by the deviation from the constraint [START_REF] Mazhar | Chrono: a parallel multi-physics library for rigid-body, flexible-body, and fluid dynamics[END_REF], one can obtain the following relation:

λ = K "» ev ep - + β » ėv ėp -« (83)
»

Jv(q) q + ėv Ψp(q) + ep - ∈ -lim ε→+0 N Uε C (λ). (84)
Here, as have been defined in section 3.4, K is the stiffness matrix and Kβ is a damping coefficient matrix. It should be noted that the constraint ( 40) is approximated by this DAI (83)(84), by introducing additional state vector [e T v , e T p ] T , of which the rate of change is determined by the DAI (83)(84).

An approach for penalty-based relaxation of frictional unilateral constraint has been proposed by Xiong et al. [START_REF] Xiong | A differential-algebraic method to approximate nonsmooth mechanical systems by ordinary differential equations[END_REF]. Special cases of (83)(84) appear in Xiong et al. [START_REF] Xiong | A differential-algebraic method to approximate nonsmooth mechanical systems by ordinary differential equations[END_REF], which are [54, eq.( 21)] for unilateral constraint and [54, eq.( 17)] for Coulomb friction. They also have pointed out that, in the case of unilateral constraint, the DAI cannot be solved in the continuous time because ėp is not included in the second row of the left-hand side of (84). To circumvent this problem, Xiong et al. [START_REF] Xiong | A differential-algebraic method to approximate nonsmooth mechanical systems by ordinary differential equations[END_REF] have introduced an artificial, non-physical term +η ėp where η > 0, which allows for the conversion from a DAI [54, eqs.(22) [START_REF] García De Jalón | Multibody dynamics with redundant constraints and singular mass matrix: existence, uniqueness, and determination of solutions for accelerations and constraint forces[END_REF]] into an ODE [54, eqs.(24) [START_REF] Glocker | Complementarity problems in multibody systems with planar friction[END_REF]]. In contrast, the approach presented by Kikuuwe and Fujimoto [START_REF] Kikuuwe | Incorporating geometric algorithms in impedance-and admittance-type haptic rendering[END_REF] is based on the implicit Euler discretization of (83)(84), of which the analytical solution can be obtained by using the relation [START_REF] Acary | Higher order Moreau's sweeping process: mathematical formulation and numerical simulation[END_REF].

A major difference between the penalty-based relaxation (83)(84) and the proposed relaxation scheme (64) is that, in (83)(84), λ is independent of M (q) and A(q). Thus, the solution is always unique in the penalty-based method. A downside of the penalty-based scheme is, as discussed in Section 3.4, its low accuracy; the original constraint (40) is never satisfied unless Φ(q, q, t) ≡ 0.

Case Studies

This section shows some examples of the system described in the form [START_REF] Marquis-Favre | A planar mechanical library in the AMESim simulation software. part I: Formulation of dynamics equations[END_REF] [START_REF] Mazhar | Chrono: a parallel multi-physics library for rigid-body, flexible-body, and fluid dynamics[END_REF] and their relaxed representations of the form (39)(64).

Planar Friction

Now we analyze the influence of the proposed relaxation scheme on the behavior of a simple system with Coulomb friction. Let us consider again the system (47) [START_REF] Silcowitz | Interactive rigid body dynamics using a projected Gauss-Seidel subspace minimization method[END_REF], which represents a point mass M > 0 subject to an external force fe ∈ R 2 and a Coulomb friction force λ ∈ R 2 . Based on the approach presented in Section 4, (47)(48) can be relaxed as follows:

M q -fe = λ (85) q + q δ ∈ -N F B2 (λ) (86)
where δ > 0 is a small positive value. Note that (86) is equivalent to the following:

λ ∈ -F sgn(q + q/δ). ( 87)

When q = 0, (87) reduces to the acceleration-level spatial friction constraint, which was discussed in [START_REF] Studer | Representation of normal cone inclusion problems in dynamics via non-linear equations[END_REF]Sec. 3.2]. By eliminating q from (85) and ( 86), one obtains the following:

λ M + fe M + q δ ∈ -N F B2 (λ), (88)
which is, by the relation ( 4), equivalent to

λ = -proj F B2 (fe + M q/δ). (89)
This means that the velocity q evolves along the following ordinary differential equation:

M q = fe -proj F B2 (fe + M q/δ). ( 90
)
One important feature of the original system (47) [START_REF] Silcowitz | Interactive rigid body dynamics using a projected Gauss-Seidel subspace minimization method[END_REF] is that q converges to 0 in finite time as long as the external force fe satisfies ∥fe∥ < F . Now, it is shown that, in the relaxed form (85)(86), or equivalently, (90), q → 0 holds asymptotically as long as ∥fe∥ < F . It can be seen by the observation of the time derivative ∥ q∥ 2 /2, which is obtained as follows:

d dt " ∥ q∥ 2 2 « = qT M " fe -proj F B2 " fe + M q δ «« = 8 < : -∥ q∥ 2 /δ if ∥fe + M q/δ∥ ≤ F F qT M " fe F - fe + M q/δ ∥fe + M q/δ∥ « otherwise. ( 91
)
It is easy to prove that, when ∥fe∥ < F , d(∥ q∥ 2 /2)/dt < 0 is always satisfied for all q ̸ = 0. Thus, q → 0 is achieved as t → ∞ as long as ∥fe∥ ≤ F , even under the non-zero external force fe. This implies that the relaxation from ( 48) to (86) preserves the characteristics of having the static friction state.

It should be noted that (90) is exactly of the same form as the "model A" of Kikuuwe et al. [START_REF] Kikuuwe | Admittance and impedance representations of friction based on implicit Euler integration[END_REF], which is obtained by the implicit discretization of [START_REF] Payr | Oblique frictional impact of a bar: Analysis and comparison of different impact laws[END_REF]. Replacing δ by the timestep size h exactly provides [31, eq.(26)]. This corresponds to the relation between the Baumgarte-like relaxation and implicit discretization, discussed in Section 4.4. In contrast to this is the penalty-based relaxation of ( 48), which can be described as follows:

λ = K(e t + β ėt ) (92) q + ėt ∈ -N F B2 (λ). (93)
This representation exactly corresponds to [54, eq.( 17)] and the "model C" of [31, eqs.(39) [START_REF] Mazhar | Chrono: a parallel multi-physics library for rigid-body, flexible-body, and fluid dynamics[END_REF]]. It should be noted that, unlike in the proposed relaxation (89), λ does not depend on fe or M in (92)(93).

Relations to Impact Laws

Let us provide a relation between the parameter ζ and the coefficient of restitution of Newton's impact law. Let us consider a one-dimensional system, in which a point mass M whose position is q ∈ R in contact with a rigid surface at q = 0. The equation of motion of this system can be described as follows:

M q = λ + fe (94) q ∈ -N R+ (λ). (95)
This is the special case of (39) [START_REF] Mazhar | Chrono: a parallel multi-physics library for rigid-body, flexible-body, and fluid dynamics[END_REF] with n = m = mp = 1, Uε = 1, Ψp(q) = q and C = R + . Based on the proposed approach, the constraint (95) can be relaxed into the following form:

q + 2ζ q δ + q δ 2 ∈ -N R+ (λ). ( 96
)
Let us analyze the behavior of the system (94)(96) that is initialized by q = 0 and q = -V at t = 0 where V > 0. This initial condition can be seen as the beginning of a collision, the end of the collision is the time at which λ = 0. The velocity Va at the end of collision is needed to obtain the coefficient of restitution. Under the condition λ > 0, the right-hand side of (96) is {0}, and thus the solution of (96) can be obtained as follows:

q δ (t) = -δV Γ ζ (t/δ) (97)
where

Γ ζ (τ ) ∆ = exp(-ζτ ) sin(τ p 1 -ζ 2 ) p 1 -ζ 2 . ( 98
)
Here, it should be noted that the function Γ ζ (τ ) does not depend of δ, and is a continuous real function of τ ≥ 0 and ζ ≥ 0, which is real even if ζ > 1 and is continuous even at ζ = 1 by the limit of ζ → 1. By using the first and second derivatives of Γ ζ (τ ), one can write the velocity and the acceleration as follows:

qδ (t) = -V Γ ′ ζ (t/δ) (99) qδ (t) = -V Γ ′′ ζ (t/δ)/δ. ( 100
)
By using them and considering the relation (94) between λ and q, one can see that the collision terminates at the time:

t = t δ ∆ = {t s.t. M qδ (t) -fe = 0} = δ(Γ ′′ ζ ) -1 " -δfe M V « (101)
where (Γ ′′ ζ ) -1 is the inverse mapping of Γ ′′ ζ . Therefore, the velocity Va at the time

t δ is Va = -V Γ ′ ζ " (Γ ′′ ζ ) -1 " -δfe M V «« . ( 102
)
Although the analytical form of the inverse map (Γ ′′ ζ ) -1 is difficult to obtain, its value at zero can be obtained as follows:

(Γ ′′ ζ ) -1 (0) = cos -1 (2ζ 2 -1) p 1 -ζ 2 . ( 103
)
By using this, when δfe is sufficiently small, one can relate ζ and the coefficient of restitution Va/V as follows:

Va V = exp - ζ cos -1 (2ζ 2 -1) p 1 -ζ 2 ! , ( 104
)
which is shown in Fig. 4. This relation has already been obtained in previous studies [18, p.89] [47] where the contact point is assumed to have linear viscoelasticity. In contrast, the present analysis is derived from the non-physical, numerical relaxation, but it provides similar results due to the similarity in the form of equations. Impact laws of two or more dimensional collisions, such as those investigated in [START_REF] Payr | Oblique frictional impact of a bar: Analysis and comparison of different impact laws[END_REF][START_REF] Stronge | Smooth dynamics of oblique impact with friction[END_REF], are more complicated, involving another one or more parameters that define the behaviors in the tangential direction. The relation between two-or threedimensional impact laws, which may involve tangential impact, and the presented relaxation scheme is left as a subject for future study.

Painlevé's Paradox

Here we show the application of the proposed relaxation scheme to the Painlevé's Paradox [START_REF] Génot | New results on Painlevé paradoxes[END_REF], in which a rigid rod is in frictional contact with a rigid horizontal surface, as illustrated in Fig. 5. The main point of this problem is that, in the rigid-body formalism, the system sometimes does not provide solutions unless discontinuous jumps in the velocity are permitted, and even when a solution exists, it may not be unique. Here we show how the presented relaxation scheme deals with this point.

As illustrated in Fig. 5, let us denote the coordinate of the gravity center of the rod by q = [q t , qn, θ] T . The position p of the contact point is written as

p ∆ = » p t pn - ∆ = » q t -ℓ cos θ qn -ℓ sin θ - ( 105
)
where ℓ is the half of the length of the rod. Without loss of generality, we assume θ ∈ [0, π/2]. We also assume that pn = 0, ṗn = 0, ṗt < 0, with which the uniqueness and existence of the solution are questioned [START_REF] Génot | New results on Painlevé paradoxes[END_REF]. The equation of motion of the system is described as follows:

M q = J(q) T λ + f (106) » ṗt pn - ∈ -lim ε→+0 N Uε F (λ) (107)
where

M ∆ = diag[m, m, mℓ 2 /3] (108) J(q) ∆ = » 1 0 ℓ sin θ 0 1 -ℓ cos θ - ( 109
) f ∆ = [0, -mg, 0] T (110) Uε ∆ = diag[1, 1/ε] (111) F ∆ = n [λ t , λn] T ∈ R 2 ˛|λ t | ≤ µλn o . ( 112
)
Here, note that (107) is the same as [START_REF] Van Der Schaft | Complementarity modeling of hybrid systems[END_REF], which was derived from (51) [START_REF] Studer | Representation of normal cone inclusion problems in dynamics via non-linear equations[END_REF]. By using the presented relaxation method, (107) is relaxed as follows:

p + 2ζ δ ṗ + 1 δ 2 » 0 pn - ∈ -lim ε→+0 N Uε F (λ) (113)
Thus, we obtain the form (65), or equivalently, r(δ, q, q) ∈ C(λ, q) (114)

where

C(λ, q) ∆ = λ + A(q) -1 " lim ε→+0 N Uε F (λ) « (115) r(δ, q, q) ∆ = ra(q, q) - 2ζ ṗt δ r b (q) - " 2ζ ṗn δ + pn δ 2 « A(q) -1 » 0 1 - (116) A(q) ∆ = J(q) T M -1 J(q) (117) ra(q, q) ∆ = mg 4 » 3 sin θ cos θ 4 -3 cos 2 θ - -mℓ θ2 » cos θ sin θ - (118) r b (q) ∆ = A(q) -1 » 1 0 - . (119)
Let us now discuss the problem (114) through a geometric consideration on C(λ, q). From (115), it is easy to see that C(λ, q) = {λ} if λ ∈ IntF. If λ ∈ ∂F\{0}, C(λ, q) is a half line that projects from λ and parallels r b (q), as illustrated in Fig. 6. If λ = 0, C(λ, q) is the half-space S(q) in Fig. 6, of which the boundary parallels r b (q). In addition, because we are assuming δ → +0, ṗt < 0, ṗn = 0 and pn = 0, r(δ, q, q) is a infinitely far point on the half-line extending from ra(q) in the direction of r b (q). This observation implies that the problem (114) has a unique solution if S(q) does not overlap F . The necessary and sufficient condition for this is

G B (q) ∆ = [1, -µ]r b (q) > 0.
(120)

Moreover, ra(q, q) ∈ S(q) is satisfied if and only if G A (q, q) ∆ = [0, 1]A(q)ra(q, q) < 0 (121) is satisfied. That is, under the condition of G B (q) > 0, the solution of ( 114) is λ = 0 if G A (q, q) < 0, and is a particular value on ∂F if G A (q, q) > 0. Note that G A (q, q) and G B (q) defined here correspond to A and B defined in [START_REF] Génot | New results on Painlevé paradoxes[END_REF], respectively. The problem is complicated when G B (q) < 0, i.e., S(q) and F overlap each other, as illustrated in Fig. 7. In such a case, r(δ, q, q) is in IntF as δ → +0, and thus λ = r(δ, q, q) ∈ IntF is a solution. This contradicts the assumption that the rod is slipping with ṗt < 0. This contradiction can be resolved by supposing that the rod instantaneously sticks due to the infinitely large λ, which constitutes a tangential impact [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control[END_REF]Sec.5.5]. This instantaneous state transition is sometimes referred to as a dynamic jamming [START_REF] Or | Investigation of Painlevé's paradox and dynamic jamming during mechanism sliding motion[END_REF]. Another important observation is that, when ra(q, q) ∈ S, i.e., G A (q, q) < 0, the problem has three solutions of λ, which are λ 1 = r(δ, q, q), a λ 2 on ∂F , and λ 3 = 0, as illustrated in Fig. 7. These solutions r a (q,q) F S(q) C(λ,q) r b (q) λ r(δ,q,q) Fig. 6 Geometric interpretation of (114) of the case of G A (q, q) > 0 and G B (q) > 0, where the unique solution λ exists.

F S(q)

λ 3 =0 r(δ,q,q) =λ 1 r a (q,q) r b (q) F ∩ S(q) λ 2 C(λ 2 ,q) C(λ 3 ,q) Fig. 7 Geometric interpretation of (114) of the case of G A (q, q) < 0 and G B (q) < 0, where (114) has three solutions: λ 1 , λ 2 and λ 3 .

correspond to sticking, slipping, and detachment, respectively. These observations are consistent with those of [START_REF] Génot | New results on Painlevé paradoxes[END_REF].

One important feature of the analysis in this section, based on the presented relaxation scheme, is that it includes both slipping and sticking by dealing with the effect of ṗt as a directional infinity. This allows geometric consideration on the problem as illustrated in Fig. 6 and Fig. 7, and both the interior and the boundary of the friction cone F can be taken into account. This point is in contrast to the analysis of [START_REF] Génot | New results on Painlevé paradoxes[END_REF][START_REF] Song | Analysis of rigid-body dynamic models for simulation of systems with frictional contacts[END_REF], which focuses the slipping state, in which the force λ is assumed to stay on the boundary of the friction cone F .

Conclusions

This article has proposed a new representation of mechanical systems involving frictional unilateral constraints. The representation is of the form of a differential algebraic inclusion (DAI), and it includes the conventional DAE representation employing Lagrange multipliers as its special case. The new representation is based on a normal cone with a singular norm metric, which is a limit of a family of normal cones with parameterized non-Euclidean norm metrics.

The paper has also presented a relaxed version of the aforementioned DAI representation. The relaxed approximation is also a DAI, but it can be algebraically equivalent to an ODE if a particular algebraic problem has a unique solution. It can be viewed as a generalization of the Baumgarte stabilization method, which is originally for equality constraints. A main benefit of the presented representation is that it can be subject to various analysis in the continuous-time domain. From the presented representation, one can derive some mathematical expressions that are similar to those in previous studies on time-stepping and event-driven schemes. That is, previous mathematical expressions in the discrete time or at different kinematic levels can now be interpreted in the continuous time in a unified framework.

The whole discussion of the paper has been based on the assumption that the Delassus operator A(q) and the mass matrix M (q) are non-singular. It has been known that singular Delassus operators [START_REF] García De Jalón | Multibody dynamics with redundant constraints and singular mass matrix: existence, uniqueness, and determination of solutions for accelerations and constraint forces[END_REF] can be caused by redundant constraints and that singular mass matrices [START_REF] Brogliato | Singular mass matrix and redundant constraints in unilaterally constrained Lagrangian and Hamiltonian systems[END_REF][START_REF] García De Jalón | Multibody dynamics with redundant constraints and singular mass matrix: existence, uniqueness, and determination of solutions for accelerations and constraint forces[END_REF] can be resulted from redundant coordinates. Further analysis will be needed to extend the scheme to be able to deal with such systems. Another open problem is to relate the presented representations to impact laws. In addition, extensions to more sophisticated models, such as those including rotational frictional contacts and area contacts with pressure distribution [START_REF] Leine | A set-valued force law for spatial Coulomb-Contensou friction[END_REF], may need to be considered in future studies.

Fig. 3

 3 Fig. 3 Relation (73) of the case of mp = 1.

Fig. 4

 4 Fig. 4 Relation (104) between the coefficient of restitution Va/V and ζ.

Fig. 5

 5 Fig. 5 Painlevé's example.

Removing +Es from (73) results in[9, equ.(16)].