Bernard Brogliato
email: bernard.brogliato@inria.fr

Daniel Goeleven
email: goeleven@hotmail.fr

Singular mass matrix and redundant constraints in unilaterally constrained Lagrangian and Hamiltonian systems

Keywords: Singular Lagrangian system, Singular Hamiltonian system, Unilateral constraint, Bilateral holonomic constraint, Complementarity conditions, Mixed linear complementarity problem, Recession cone, Convex analysis, Redundant constraints, Legendre-Fenchel transformation

This article deals with the analysis of the contact complementarity problem for Lagrangian systems subjected to unilateral constraints, and with a singular mass matrix and redundant constraints. Previous results by the authors on existence and uniqueness of solutions of some classes of variational inequalities are used to characterize the well-posedness of the contact problem. Criteria involving conditions on the tangent cone and the constraints gradient are given. It is shown that the proposed criteria easily extend to the case where the system is also subjected to a set of bilateral holonomic constraints, in addition to the unilateral ones. In the second part, it is shown how basic convex analysis may be used to show the equivalence between the Lagrangian and the Hamiltonian formalisms when the mass matrix is singular.

as vectors with real components). The Lagrange dynamics is written when constraints are present as M(q) q + F (q, q, t) = ∇h(q)λ [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF] where M(q) = M(q) T is the mass matrix, F : R n × R n × R + → R n contains centripetal, Coriolis forces, forces that derive from a potential, and possible external forces like control inputs or perturbations, q is the vector of generalized accelerations. The system is free, or subjected to constraints (bilateral h(q) = 0, or unilateral h(q) ≥ 0) represented by a differentiable function h : R n → R m with gradient = ∇h(q), and λ ∈ R m is the vector of Lagrange multipliers associated with the constraints. The mass matrix M(q) ∈ R n×n is often supposed to be positive definite. If in addition the constraints are functionally independent, then both the acceleration and the multipliers can be calculated uniquely. In the case of bilateral constraints, this boils down to solving a linear equation obtained with the Karush-Kuhn-Tucker matrix, stemming from the so-called index-1 DAE formulation [START_REF] Garcia De Jalon | Multibody dynamics with redundant constraints and singular mass matrix: existence, uniqueness, and determination of solutions for accelerations and constraints forces[END_REF][START_REF] Laulusa | Review of classical approaches for constraint enforcement in multibody systems[END_REF][START_REF] Simeon | Computational Flexible Multibody Dynamics. A Differential-Algebraic Approach[END_REF]:

M(q) -∇h(q) ∇h(q) T 0 q λ = -F (q, q, t) -( d dt ∇h(q) T ) q .

(2)

If both M(q) and ∇h(q) are full-rank matrices, one calculates λ using the Schur complement of the KKT matrix in (2): λ = A(q) -1 H 0 (q, q, t), where H 0 (q, q, t) = ∇h(q) T M(q) -1 F (q, q, t) -d dt ∇h(q) T q and A(q) = ∇h(q) T M(q) -1 ∇h(q) ∈ R m×m . When M(q) and/or ∇h(q) become singular, uniqueness of q and of ∇h(q)λ is guaranteed for any F (q, q, t) if and only if [START_REF] Garcia De Jalon | Multibody dynamics with redundant constraints and singular mass matrix: existence, uniqueness, and determination of solutions for accelerations and constraints forces[END_REF]:

ker M(q) ∩ ker ∇h(q) T = {0},

that is equivalent to rank[M(q)∇h(q)] = n (the result is stated assuming that the system lives on the constraint submanifold {(q, q) : h(q) = 0, ∇h(q) T q = 0}, so that the second equality in (2) is satisfied). Uniqueness of λ requires that ∇h(q) be full rank. The wellposedness of (2), together with the initial position and velocity compatibility with the constraint manifold, guarantees the existence and uniqueness of solutions to the constrained Lagrangian dynamics. Numerical methods also have to solve discretized versions of (2).

In the case of unilateral constraints, inducing complementarity conditions 0 ≤ λ ⊥ h(q) ≥ 0, and under the assumption that M(q) has rank n, one computes the so-called contact linear complementarity problem (LCP) [START_REF] Brogliato | Inertial couplings between unilateral and bilateral holonomic constraints in frictionless Lagrangian systems[END_REF][START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF]: 0 ≤ λ ⊥ A(q)λ + H 0 (q, q, t) ≥ 0 [START_REF] Agrawal | Inertia matrix singularity of planar series-chain manipulators[END_REF] where A(q) is called the Delassus' matrix in the literature on unilaterally constrained systems. The process followed to obtain the contact LCP may be seen as an extension of the index reduction of the bilateral DAE case. The conditions for the solvability, the feasibility and the uniqueness of solutions of the LCP in (4) are given in [START_REF] Brogliato | Inertial couplings between unilateral and bilateral holonomic constraints in frictionless Lagrangian systems[END_REF]Propositions 1,[START_REF] Addi | A qualitative mathematical analysis of a class of linear variational inequalities via semi-complementarity problems[END_REF][START_REF] Agrawal | Inertia matrix singularity of series-chain spatial manipulators with point masses[END_REF], building on previous results by Moreau [START_REF] Moreau | Les liaisons unilatérales et le principe de Gauss[END_REF][START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF] and Lötstedt [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF]. It is seen from the results in [21, §1, 2] and [START_REF] Brogliato | Inertial couplings between unilateral and bilateral holonomic constraints in frictionless Lagrangian systems[END_REF]Propositions 1,[START_REF] Addi | A qualitative mathematical analysis of a class of linear variational inequalities via semi-complementarity problems[END_REF][START_REF] Agrawal | Inertia matrix singularity of series-chain spatial manipulators with point masses[END_REF] that despite the linear equation in [START_REF] Addi | A qualitative mathematical analysis of a class of linear variational inequalities via semi-complementarity problems[END_REF] and the contact LCP in (4) involving the same terms (in particular the Delassus's matrix is present in both), their well-posedness conditions differ. For instance, solvability of (2) relies on the above rank condition from a basic result in linear algebra (see, e.g. [START_REF] Bernstein | Matrix, Mathematics. Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF]Proposition 6.1.7]), while the solvability of (4) relies on complementarity theory [START_REF] Brogliato | Inertial couplings between unilateral and bilateral holonomic constraints in frictionless Lagrangian systems[END_REF]Proposition 3].

The above results on the contact LCP all rely on the assumption that M(q) is positive definite. However, it is known that the rank of M(q) may drop to r < n, hence it becomes positive semi-definite only [START_REF] Agrawal | Inertia matrix singularity of series-chain spatial manipulators with point masses[END_REF][START_REF] Agrawal | Inertia matrix singularity of planar series-chain manipulators[END_REF][START_REF] Balakrishnan | Vibrating systems with singular mass-inertia matrices[END_REF][START_REF] Bhat | Second-order systems with singular mass matrix and an extension of Guyan reduction[END_REF][START_REF] Havelkova | A geometric analysis of dynamical systems with singular Lagrangians[END_REF][START_REF] Hurtado | Lagrangian mechanics of overparameterized systems[END_REF][START_REF] Schutte | New approach to the modeling of complex multibody dynamical systems[END_REF][START_REF] Simeon | Computational Flexible Multibody Dynamics. A Differential-Algebraic Approach[END_REF][START_REF] Udwadia | Explicit equations of motion for constrained mechanical systems with singular mass matrices and applications to multi-body dynamics[END_REF][START_REF] Udwadia | Equations of motion for general constrained systems in Lagrangian mechanics[END_REF]. The analysis of bilaterally constrained systems when the mass matrix has low rank has been made in [20-22, 27, 37, 42, 43]. The analysis of mixed bilaterally/unilaterally constrained systems is made in [START_REF] Brogliato | Inertial couplings between unilateral and bilateral holonomic constraints in frictionless Lagrangian systems[END_REF]. In this case one obtains a contact LCP with an augmented Delassus' matrix A c (q), which is distorted by the bilateral constraints effects. It is shown that even if M(q) has full rank, the contact LCP may be ill-posed, i.e. A c (q) may be positive semi-definite only. Criteria are provided in [START_REF] Brogliato | Inertial couplings between unilateral and bilateral holonomic constraints in frictionless Lagrangian systems[END_REF] which guarantee that A c (q) is positive definite.

The case when unilateral constraints act on the system, and when r < n, however, has received no attention, to the best of the authors' knowledge. This is the topic of this article, which is dedicated to the analysis of the contact complementarity problem when low rank conditions on both M(q) and ∇h(q) are present. The article is organized as follows: Some tools from convex analysis are recalled in the next paragraph. The contact complementarity problem, which takes the form of a mixed linear complementarity problem (MLCP), is constructed in Sect. 2. Its well-posedness is analyzed in Sect. 3. We show in Sect. 4 that the presence of holonomic bilateral constraints, in addition to the unilateral ones, may be treated within the foregoing framework. Section 5 shows how the MLCP may be transformed into an LCP, and paves the way towards Sect. 6 which is dedicated to proving the equivalence between degenerated Lagrangian and Hamiltonian systems with unilateral constraints. Two simple examples are treated in Sect. 7. Conclusions end the article in Sect. 8, and some mathematical results are given in Appendices A and B.

Notation and definitions.

Let h : R n → R m be differentiable. Its gradient at q is ∇h(q) ∈ R n×m , and its Jacobian is ∂h ∂q (q) = ∇h(q) T ∈ R m×n , where the superscript T means transpose. Given a vector x ∈ R n , x ≥ 0 means x i ≥ 0 for all its entries x i (same for > 0). Given a matrix A ∈ R n×m , ker(A) is its kernel and Im(A) is its range, or image. Given two vectors x and y, the product x, y = x T y.

Mathematical tools in convex analysis. The next mathematical tools are mainly taken from classical references [START_REF] Bertsekas | Convex Optimization Theory[END_REF][START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF][START_REF] Hiriart-Urruty | Fundamentals of Convex Analysis[END_REF][START_REF] Rockafellar | Convex Analysis[END_REF]. Let f : R n → R ∪ {+∞} be a proper convex and lower semi-continuous function, we denote by dom(f

) := {x ∈ R n : f (x) < +∞} the domain of the function f (•). The epigraph of f (•) is the set epi(f ) := {(x, α) ∈ R n × R : α ≥ f (x)}.
Recall that the Fenchel transform f (•) of f (•) is the proper, convex and lower semicontinuous function defined by

∀z ∈ R n : f (z) = sup x∈ dom(f ) x, z -f (x) .
(

) 5
The subdifferential ∂f (x

) of f (•) at x ∈ R n is defined by ∂f (x) = ω ∈ R n : f (v) -f (x) ≥ ω, v -x , ∀v ∈ R n .
We denote by Dom(∂f

) := {x ∈ R n : ∂f (x) = ∅}
Dom(∂f ) ⊂ dom(f ) ⊂ Dom(∂f ). ( 6
)
Let x 0 be any element in the domain dom(f The function f ∞ : R n → R ∪ {+∞} is a proper convex and lower semi-continuous function which describes the asymptotic behavior of f (•).

) of f (•), the recession function f ∞ (•) of f (•) is defined by ∀x ∈ R n : f ∞ (x) = lim λ→+∞ 1 λ f (x 0 + λx).
Let K ⊂ R n be a nonempty closed convex set. Let x 0 be any element in K. The recession cone of K is defined by

K ∞ = λ>0 1 λ (K -x 0 ) = u ∈ R n : x + λu ∈ K, ∀λ ≥ 0, ∀x ∈ K .
The set K ∞ is a nonempty closed convex cone that is described in terms of the directions which recede from K. The indicator function of

K is Ψ K (x) = 0 if x ∈ K, Ψ K (x) = +∞ if x / ∈ K. Example 1 Let f (x) = x 2 if x<0, 2x if x≥0. . Then f ∞ (x) = +∞ if x<0, x if x≥0. . Indeed, if x < 0 then lim λ→+∞ 1 λ f (λx) = lim λ→+∞ λx 2 = +∞ while if x ≥ 0 then lim λ→+∞ 1 λ f (λx) = lim λ→+∞ x = x.
Example 2 Let M ∈ R n×n be a positive semi-definite and symmetric matrix, q ∈ R n and f

(x) = 1 2 x T Mx + q T x. Then f ∞ (x) = q T x if x∈ker(M), +∞ if x / ∈ker(M). . Example 3 Let K = {(x 1 , x 2 ) ∈ R 2 : x 1 > 0, x 1 x 2 ≥ 1}. Then K ∞ = {(x 1 , x 2 ) ∈ R 2 : x 1 ≥ 0, x 2 ≥ 0}
. This is depicted in Fig. 1a.

Example 4 Let K = {(x 1 , x 2 ) ∈ R 2 : x 1 ≥ 2, -x 1 + 4 ≤ x 2 ≤ x 1 + 4}. Then K ∞ = {(x 1 , x 2 ) ∈ R 2 : x 1 ≥ 0, -x 1 ≤ x 2 ≤ x 1 }
. This is depicted in Fig. 1b.

Let us here recall some important properties of the recession function and recession cone (see, e.g. [START_REF] Rockafellar | Variational Analysis[END_REF]Chaps. 8 and 9], or [START_REF] Bertsekas | Convex Optimization Theory[END_REF]Proposition 1.4.8]): Proposition 1 (a) Let f 1 : R n → R ∪ {+∞} and f 2 : R n → R ∪ {+∞} be two proper, convex and lower semi-continuous functions. Suppose that

f 1 + f 2 is proper. Then ∀x ∈ R n : (f 1 + f 2 ) ∞ (x) = (f 1 ) ∞ (x) + (f 2 ) ∞ (x). ( 7
)
Fig. 2 Equilibrium positions (b) Let f : R n → R ∪ {+∞} be a proper, convex and lower semi-continuous function and let K be a nonempty closed convex set such that f

+ Ψ K is proper (equivalently, dom(f ) ∩ K is nonempty). Then ∀x ∈ R n : (f + Ψ K ) ∞ (x) = f ∞ (x) + (Ψ K ) ∞ (x). (8)
(c) Let K ⊂ R n be a nonempty, closed and convex set. Then

∀x ∈ R n : (Ψ K ) ∞ (x) = Ψ K∞ (x), (9)
(∀x ∈ K, e ∈ K ∞ ) : x + e ∈ K. ( 10
) (d) If K ⊂ R n is a nonempty closed and convex cone, then K ∞ = K. Example 5 Let K = [a, b], (a, b ∈ R, a < b). The set K is bounded and thus K ∞ = {0}. It results that (Ψ K ) ∞ (x) = Ψ {0} = +∞ if x = 0, 0 ifx = 0.
Example 6 Let us consider a ball M subjected to its weight G and constrained to stay in a box K. In case of Fig. 2a, all directions of K ∞ make an obtuse angle with the direction of the applied force, i.e. (∀x ∈ K ∞ ) : G T x < 0. In that case, we note that the ball admits a unique equilibrium position. The condition (∀x ∈ K ∞ ) : G T x < 0 may also be written as

∀x ∈ R 2 : G T x + Ψ K∞ (x) < 0. ( 11
)
Let us now set U(x) = G T x + Ψ K (x). The mapping x → G T x is the potential energy corresponding to the gravity force while Ψ K (x) is a penalty energy associated with the constraints. Using (8), we see that U ∞ (x) = G T x + Ψ K∞ (x), and thus condition [START_REF] Bertsekas | Convex Optimization Theory[END_REF] is equivalent to

(∀x ∈ R 2 ) : U ∞ (x) < 0.
In case of Fig. 2b, there exists a direction v ∞ ∈ K ∞ making an acute angle with the direction of the applied force, i.e. G T v ∞ > 0, and we see clearly that there is no possible equilibrium position for the ball. This simple example shows the importance of the tools of recession to determine equilibrium conditions in the case of systems subjected to unilateral constraints.

Example 7 A polyhedral set is defined as the set of solutions to a finite set of linear inequalities, i.e.

P (A, b) = x ∈ R n : Ax ≥ b , ( 12
)
where A ∈ R m×n and b ∈ R m . It is clear that P (A, b) is a closed convex set. A function f : R n → R ∪ {+∞} whose epigraph is polyhedral is called a polyhedral function [35, p. 172]. Suppose that P (A, b) = ∅. Then there exists x 0 ∈ R n : Ax 0 ≥ b and

P (A, b) ∞ = λ>0 1 λ P (A, b) -x 0 = P (A, 0). ( 13
)
Indeed,

x ∈ P (A, b) ∞ ⇔ (∀λ > 0) : A(λx + x 0 ) ≥ b ⇔ (∀λ > 0) : λAx + Ax 0 -b ≥ 0 ⇔ Ax ≥ 0,
where one checks that (10) is satisfied.

Let K be a subset of R n . The polar cone of K is defined by

K + = w ∈ R n : w, v ≤ 0, ∀v ∈ K .
The dual cone of K is defined by the relation

K * = -K + .
Let K be a nonempty subset of R n . To every point x ∈ K, we can associate the set

T K (x) = y ∈ R n : ∃{x n } ⊂ K, {λ n } ⊂ R + , λ n 0, x n → x, x n -x λ n → y
called the tangent cone of K at x. The normal cone to K at x is defined as the polar cone of T K (x), i.e.

N K (x) = T K (x) + .
If in addition K is convex then

N K (x) = z ∈ R n : z, v -x ≤ 0 for all v ∈ K ( 14
)
and N K (x) + = T K (x). It represents the cone generated by the normals to the set K at x, and is an outwards cone. If K is a nonempty closed convex set, then one has

∂Ψ K (•) = N K (•)
and both the tangent and normal cones are convex closed sets.

into account only the active constraints). We further assume that the acceleration is rightcontinuous, i.e. q(t) = q(t + ). The same holds for the multiplier λ(•). Developing a bit more, one obtains

0 ≤ λ(t) ⊥ ∇h(q) T q(t) + d dt ∇h q(t) T q(t) ≥ 0. ( 15
)
Let us denote x = q(t), λ = λ(t) and H (q, q) = d dt (∇h(q(t)) T ) q(t). Then ( 1) and ( 15) give rise to

M(q)x + F (q, q, t) -∇h(q)λ = 0, 0 ≤ λ ⊥ ∇h(q) T x + H (q, q) ≥ 0. ( 16
)
The set of equations in ( 16) is a mixed LCP (MLCP) [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF]Definition 12.22], which is a combination between an LCP and linear equations. If M(q) has rank n it reduces to the LCP in ( 4). This MLCP is the counterpart of (2) for unilaterally constrained systems. In order to analyze the well-posedness of this MLCP, let us transform it into a suitable form. First, it follows from a basic convex analysis result that

0 ≤ λ ⊥ ∇h(q) T x + H (q, q) ≥ 0 ⇔ λ ∈ -∂Ψ K ∇h(q) T x + H (q, q) (17)
with

K = R m + .
Thus the MLCP is rewritten as the inclusion

M(q)x + F (q, q, t) ∈ -∇h(q)∂Ψ K ∇h(q) T x + H (q, q) . ( 18
)
Here q and q are considered as parameters, the argument of the subdifferential being x. We set

K = z ∈ R m : z + H (q, q) ∈ K = K -H (q, q).
Let us make the following:

Assumption 1 We suppose that Im ∇h(q) T ∩ K = ∅.
Let us also set

Φ = x ∈ R n : ∇h(q) T x + H (q, q) ∈ K . ( 19
)
Assumption 1 entails that Φ is nonempty. It is also clear that Φ is closed and convex. This set (as well as K) should be denoted with arguments as Φ(q, q); however, to lighten the notation, we drop the arguments. We have

Ψ Φ (x) = Ψ K ∇h(q) T x + H (q, q) = Ψ K ∇h(q) T x .
Let us first consider Ψ Φ (•) as the composition mapping

Ψ Φ (x) = Ψ K ∇h(q) T x .
The set K is a polyhedral set in R m (see [START_REF] Bhat | Second-order systems with singular mass matrix and an extension of Guyan reduction[END_REF]). Indeed,

K = P I m×m , -H (q, q) .
It results that the indicator function Ψ K (•) is a polyhedral function. Then, by Theorem 1 in Appendix B and Assumption 1, it follows that (here

h(•) = Ψ K (•), A = ∇h(q) T and f (•) = Ψ Φ (•)) ∂Ψ Φ (x) = ∇h(q)∂Ψ K ∇h(q) T x .
We have also

∂Ψ K ∇h(q) T x = ∂Ψ K-H (q, q) ∇h(q) T x = ∂Ψ K ∇h(q) T x + H (q, q) ,
and thus

∇h(q) ∂Ψ K ∇h(q) T x + H (q, q) = ∂Ψ Φ (x). ( 20
)
Let us also notice that a necessary condition for the MLCP in ( 16) to possess at least one solution is that Φ is nonempty. This is guaranteed by Assumption 1. 1 Thus we have shown that the MLCP in ( 16) is equivalent, under the constraint qualification of Assumption 1, to the generalized equation

M(q)x + F (q, q, t) ∈ -∂Ψ Φ (x) ⇔ Find x ∈ Φ : M(q)x + F (q, q, t), v -x ≥ 0 for all v ∈ Φ, ( 21
)
where the second formalism is a variational inequality.

Remark 1 Usually the contact problem is studied in terms of the multiplier from the LCP in ( 4). Here we address the problem by studying the existence and uniqueness of the acceleration x in ( 21).

Well-posedness of the MLCP

If rank(M(q)) = n then Theorem 6 in [START_REF] Addi | A qualitative mathematical analysis of a class of linear variational inequalities via semi-complementarity problems[END_REF] applies to show the existence and uniqueness of x (hence providing an alternative proof of the result shown by Moreau in [START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF]). By elementary convex analysis, one finds that the generalized equation ( 21) has a unique solution given by an orthogonal projection in the kinetic metric (see [15, §4.4] for further developments):

x = proj M(q) Φ; -M(q) -1 F (q, q, t) . ( 22
)
However, the case of interest here is when r < n. Using results from [START_REF] Addi | A qualitative mathematical analysis of a class of linear variational inequalities via semi-complementarity problems[END_REF], we are going to characterize the existence (solvability) and uniqueness of solutions to [START_REF] Garcia De Jalon | Multibody dynamics with redundant constraints and singular mass matrix: existence, uniqueness, and determination of solutions for accelerations and constraints forces[END_REF]. Let us first recall the following result, where we concatenate Theorem 3 and Corollary 4 in [START_REF] Addi | A qualitative mathematical analysis of a class of linear variational inequalities via semi-complementarity problems[END_REF] (Corollary 5 in [START_REF] Addi | A qualitative mathematical analysis of a class of linear variational inequalities via semi-complementarity problems[END_REF] cannot be applied unless H (q, q) = 0 in which case Φ is a closed convex cone). They concern variational inequalities of the form: Find u ∈ R n such that

Mu + q, v -u + ϕ(v) -ϕ(u) ≥ 0, ∀v ∈ R n ( 23
)
where M ∈ R n×n is a real matrix, q ∈ R n a vector and ϕ : R n → R ∪ {+∞} a proper convex and lower semicontinuous function. The analogy with the generalized equation in ( 21) is 1 Another constraint qualification guaranteeing nonemptyness of K, hence of Φ, is range (∇h(q) T ) -R m + = R m , which does not imply that ∇h(q) T has rank m [START_REF] Brogliato | Well-posedness results for non-autonomous complementarity systems[END_REF].

clear taking ϕ(•) = Ψ Φ (•), thus restricting the variation of v to Φ, q = F (q, q, t) and M = M(q). The problem in ( 23) is denoted as VI(M, q, ϕ) in the next proposition. We also set

K(M, ϕ) = x ∈ R n |Mx ∈ dom(ϕ ∞ ) * . ( 24
)
Note that (dom(ϕ ∞ )) * is the dual cone of the domain of the recession function ϕ ∞ while (dom(ϕ)) ∞ is the recession cone of dom(ϕ).

Proposition 2 [START_REF] Addi | A qualitative mathematical analysis of a class of linear variational inequalities via semi-complementarity problems[END_REF] Let ϕ : R n → R ∪ {+∞} be a proper, convex and lower semicontinuous function with closed domain, and suppose that M ∈ R n×n is positive semi-definite and symmetric.

(a) If (dom(ϕ)) ∞ ∩ ker{M} ∩ K(M, ϕ) = {0} then for each q ∈ R n , problem VI(M, q, ϕ) has at least one solution. (b) Suppose that (dom(ϕ)) ∞ ∩ ker{M} ∩ K(M, ϕ) = {0}. If q, v + ϕ ∞ (v) > 0, ∀v ∈ dom(ϕ) ∞ ∩ ker{M} ∩ K(M, ϕ), v = 0, ( 25
)
then problem VI(M, q, ϕ) has at least one solution. (c) If u 1 and u 2 denote two solutions of problem VI(M, q, ϕ) then

u 1 -u 2 ∈ ker{M} ( 26
)
and

q, u 1 -u 2 = ϕ(u 2 ) -ϕ(u 1 ). ( 27
) (d) If ϕ(x + z) = ϕ(x), ∀x ∈ dom(ϕ), z ∈ ker{M}
and q, e = 0, ∀e ∈ ker{M}, e = 0, then problem VI(M, q, ϕ) has at most one solution.

Notice that the function ϕ(•) will never be strictly convex in our case (it is an indicator function) so that the strict convexity argument of Corollary 4 in [START_REF] Addi | A qualitative mathematical analysis of a class of linear variational inequalities via semi-complementarity problems[END_REF] does not apply and is not recalled. Let us now apply Proposition 2 to [START_REF] Garcia De Jalon | Multibody dynamics with redundant constraints and singular mass matrix: existence, uniqueness, and determination of solutions for accelerations and constraints forces[END_REF]. We need to perform some calculations on the recession and dual cones first to prepare the result statement.

We have

dom(ϕ) = dom(Ψ Φ ) = Φ = P ∇h(q) T , -H (q, q) .
Thus using [START_REF] Brezis | Analyse Fonctionnelle. Théorie et Applications[END_REF],

dom(ϕ) ∞ = Φ ∞ = P ∇h(q) T , 0 = x ∈ R n : ∇h(q) T x ≥ 0 . ( 28
)
Let S = {q ∈ R n : h(q) ≥ 0} be the admissible domain of configurations. Suppose that all m constraints are active, i.e. h(q(t)) = 0 for all t ∈ I where I is an interval (not reduced to a singleton), and suppose that the Mangasarian-Fromovitz constraint qualification (MFCQ) holds (see Appendix A). This constraint qualification reduces here to the existence of v ∈ R n such that ∇h(q) T v > 0. Then we may write

x ∈ R n : ∇h(q) T x ≥ 0 = T S (q). ( 29
)
From ( 9) one has ϕ ∞ (•) = Ψ Φ∞ (•). Thus the dual cone in ( 24) is given by

dom(ϕ ∞ ) * = Φ * ∞ .
Using the well known Farkas' lemma, we get

Φ * ∞ = x ∈ R n : ∃α ≥ 0 such that x = ∇h(q)α = x ∈ R n : ∃α 1 , . . . , α m ≥ 0 such that x = m i=1 ∇h i (q)α i . ( 30
)
Still assuming that all m constraints are active and that the MFCQ holds, one obtains

Φ * ∞ = dom(ϕ ∞ ) * = -N S (q).
Consequently, using (24), we have

K(M, ϕ) = x ∈ R n : ∃ α ≥ 0 such that M(q)x = ∇h(q)α . ( 31
)
In the next propositions, q(t), q(t) and t are considered as given parameters. Recall that by a solution of the MLCP one means a vector x ∈ R n (and not a multiplier λ), since we are studying the variational inequality [START_REF] Garcia De Jalon | Multibody dynamics with redundant constraints and singular mass matrix: existence, uniqueness, and determination of solutions for accelerations and constraints forces[END_REF].

Remark 2 When the constraints are bilateral (i.e. h(q) = 0), the tangent cone T S (q) and the normal N S (q) cone reduce to the usual tangent and normal subspaces to the submanifold Σ = {q ∈ R n : h(q) = 0} = S, respectively. Then T S (q) = ker(∇h(q) T ) and N S (q) = (T S (q)) ⊥ = Im(∇h(q)).

Proposition 3 (MLCP solvability) Assume that the m unilateral constraints are active, that Assumption 1 holds and that one of the following conditions is satisfied:

(i) P (∇h(q) T , 0) ∩ ker(M(q)) = {0}, (ii) P (∇h(q) T , 0) ∩ ker(M(q)) = {0} and F (q, q, t) T v > 0 for all v ∈ P (∇h(q) T , 0) ∩ ker(M(q)), v = 0.

Then the MLCP in [START_REF] Brogliato | Kinetic quasi-velocities in unilaterally constrained Lagrangian mechanics with impacts and friction[END_REF] has at least one solution.

Proof It follows from Proposition 2(a) and (b). Notice that

Φ ∞ ∩ ker M(q) ∩ K M(q), Ψ Φ = x ∈ R n : ∇h(q) T
x ≥ 0, M(q)x = 0 and M(q)x = ∇h(q)α for some α ≥ 0 = x ∈ R n : ∇h(q) T x ≥ 0 and M(q)x = 0 and 0 = ∇h(q)α for some α ≥ 0

= x ∈ R n : ∇h(q) T x ≥ 0 and M(q)x = 0 = Φ ∞ ∩ ker M(q) .
The results follow using [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF].

Fig. 3 A system satisfying Proposition 4(ii);

S = {q ∈ R n : h(q) ≥ 0}
If the MFCQ is satisfied and all m constraints are active then using [START_REF] Mcclamroch | Feedback stabilization and tracking of constrained robots[END_REF] we may rewrite the previous result as follows:

Proposition 4 (MLCP solvability) Suppose that all m unilateral constraints are active, i.e. h(q) = 0, and suppose that the MFCQ holds. Under Assumption 1 and if one of the following conditions is satisfied:

(i) T S (q) ∩ ker(M(q)) = {0}, (ii) T S (q) ∩ ker(M(q)) = {0} and F (q, q, t) T v > 0 for all v ∈ T S (q) ∩ ker(M(q)), v = 0, then the MLCP in [START_REF] Brogliato | Kinetic quasi-velocities in unilaterally constrained Lagrangian mechanics with impacts and friction[END_REF] has at least one solution.

Figure 3 depicts a case when Proposition 4(ii) is applicable, where the tangent and normal cones are the dashed areas. Exchanging ker(M(q)) with Im(M(q)) and vice versa, one obtains T S (q) ∩ ker(M(q)) = {0}.

Remark 3 (i) Condition (i) of Proposition 3 holds if M(q) is invertible since in this case ker(M(q)) = {0}. (ii) If M(q) is not invertible and m = 1 then P (∇h(q) T , 0) ∩ ker(M(q)) = ∅. Indeed, dim(ker(M(q))) = ν ≥ 1; let ξ 1 , . . . , ξ ν be a basis of ker(M). Setting α i = 1 when ∇h(q) T ξ i ≥ 0 and α i = -1 when ∇h(q) T ξ i < 0, we see that

x = α 1 ξ 1 + • • • + α ν ξ ν ∈ ker M(q) and ∇h(q) T x = α 1 ∇h(q) T ξ 1 + • • • + α ν ∇h(q) T ξ ν ≥ 0.

Thus

x ∈ P ∇h(q) T , 0 ∩ ker(M).

Fig. 4 A system satisfying Proposition 4(ii); S = {q ∈ R n : h(q) ≥ 0} Thus if m = 1 (there is only one active constraint at q), one has to resort to condition (ii) of Proposition 4. This is illustrated in Fig. 4. In this case the tangent cone is the dashed half space, while the normal cone is a half line generated by -∇h(q). Proposition 5 (MLCP uniqueness of solutions) (a) Let q1 (t), q2 (t), and λ 1 (t), λ 2 (t) be solutions of the MLCP [START_REF] Brogliato | Kinetic quasi-velocities in unilaterally constrained Lagrangian mechanics with impacts and friction[END_REF].

Then ⎧ ⎪ ⎨ ⎪ ⎩ q1 (t) -q2 (t) ∈ ker(M(q)), λ 1 (t) -λ 2 (t) ∈ ker(∇h(q)), F (q, q, t) T ( q1 (t) -q2 (t)) = 0. (32)
(b) Suppose that ∇h(q) T z = 0, for all z ∈ ker(M(q)), F (q, q, t) T e = 0 for all e ∈ ker(M(q)), e = 0. (

) 33
Then the MLCP in [START_REF] Brogliato | Kinetic quasi-velocities in unilaterally constrained Lagrangian mechanics with impacts and friction[END_REF] has at most one solution q(t). Consequently, if solvability holds, ∇h(q)λ(t) is unique.

Proof The inclusions in (32) are a direct consequence of (c) in Proposition 2. Our first assumption in (b) implies that for all x ∈ Φ, z ∈ ker(M(q)):

∇h(q) T (x + z) + H (q, q) = ∇h(q) T x + H (q, q) ∈ K, see [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]. One has ϕ(x + z) = Ψ Φ (x + z) and ϕ(x) = Ψ Φ (x) = 0. Thus,

ϕ(x + z) = Ψ Φ (x + z) = ϕ(x).
It results that the first condition in (d) of Proposition 2 is satisfied. The second condition in (d) of Proposition 2 is a direct consequence of our second assumption in (b). The uniqueness result comes thus directly from Proposition 2.

The first assumption in (b) is equivalent to ∀i ∈ {1, . . . , m} : ∇h i (q) ∈ ker M(q) ⊥ .

Moreover, the matrix M(q) is symmetric and thus ker(M(q)) ⊥ = Im(M(q)). It results that the first assumption in (b) is equivalent to ∀i ∈ {1, . . . , m} : ∇h i (q) ∈ Im M(q) .

The second condition in ( 32) is equivalent to the uniqueness of ∇h(q)λ, the first one to the uniqueness of M(q) q(t), and the third one to the uniqueness of F (q, q, t) T q(t).

Remark 4

The first condition of Proposition 5(b) is equivalently rewritten as ker(M(q)) ⊆ ker(∇h(q) T ) ⇔ Im(M(q)) ⊇ Im(∇h(q)). Obviously, if M(q) has rank n, this is always satisfied. The condition requires that the rank of ∇h(q) be less or equal to that of M(q). In a system with two independent constraints like that of Fig. 3, the uniqueness condition of Proposition 5 thus implies that rank(M(q)) ≥ 2. If the system has more than 2 degrees of freedom, it is then allowed that M(q) be singular. Let us go a step further. The above condition is equivalent to z / ∈ ker(∇h(q) T ) ⇒ z / ∈ ker(M(q)). Moreau's viability Lemma (see [14, Propositions 5.1 and 5.2], [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF]) states that the right velocity q(t + ), which is in our case simply equal to q(t) since we do not consider impacts, has to satisfy q(t) ∈ T S (q(t)) for q(t) ∈ S to be satisfied (such velocities may be named viable). From ( 29) it follows that ∇h(q) T q(t) ≥ 0. Thus we deduce that q(t) / ∈ ker(∇h(q(t)) T ) ⇔ ∇h(q(t)) T q(t) > 0 ⇒ q(t) / ∈ ker(M(q(t))): velocities that point strictly inside the admissible set S and which satisfy the condition for uniqueness of accelerations have to produce positive kinetic energy. The solvability condition of Proposition 4(i) yields close conclusions. Indeed, it states that the nonzero viable velocities q(t) (i.e. the velocities that belong to the tangent cone) are not in ker(M(q)), consequently they produce positive kinetic energy. These are observations similar to those made in [START_REF] Garcia De Jalon | Multibody dynamics with redundant constraints and singular mass matrix: existence, uniqueness, and determination of solutions for accelerations and constraints forces[END_REF] for bilateral constraints (end of Sect. 2 in that paper). The other conditions of Propositions 4 and 5 impose some conditions on the generalized forces for the viable velocities which produce zero kinetic energy.

Mixed bilateral/unilateral frictionless constraints

It is often the case in practice that the system under study is subjected not only to m unilateral, but also to p bilateral constraints. For the sake of clarity of the presentation, let us denote the set of bilateral holonomic constraints as f (q) = 0, with f : R n → R p a differentiable function. The Lagrange multipliers associated with the bilateral constraints are denoted as μ ∈ R p . Therefore, the right-hand side of the dynamics in (1) becomes equal to ∇h(q)λ + ∇f (q)μ and the MLCP in [START_REF] Brogliato | Kinetic quasi-velocities in unilaterally constrained Lagrangian mechanics with impacts and friction[END_REF] becomes

⎧ ⎪ ⎨ ⎪ ⎩ M(q)x + F (q, q, t) -∇h(q)λ -∇f (q)μ = 0, 0 ≤ λ ⊥ ∇h(q) T x + H (q, q) ≥ 0, ∇f (q) T x + d dt (∇f (q) T ) q = 0 (34)
with x = q(t). It has been shown above that ∇h(q)λ ∈ -∂Ψ Φ ( q) with Φ in [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]. Quite similarly one may prove that ∇f (q)μ ∈ -∂Ψ Ω ( q), 2 with Ω = {x ∈ R n : ∇f (q) T x + d dt (∇f (q) T ) q = 0}. Then, ∇h(q)λ +∇f (q)μ ∈ -∂Ψ Φ ( q) -∂Ψ Ω ( q) = -∂Ψ Φ∩Ω ( q), where the last equality is obtained under the condition that Φ ∩ Ω = ∅ (which is an equivalent of Assumption 1), see [START_REF] Rockafellar | Convex Analysis[END_REF]Theorem 23.8]. Proceeding along the same steps as in 2 The minus sign is added here just to be coherent with the unilateral case.

Sect. 3, and replacing Φ by Φ ∩ Ω, one notes in particular that (dom(ϕ)) ∞ = (Φ ∩ Ω) ∞ = ker(∇f (q) T )∩T S (q), where we assumed as above that the m unilateral constraints are active and satisfy the MFCQ. We therefore obtain an extension of Proposition 4:

Proposition 6 (MLCP solvability) Suppose that all m unilateral constraints are active, i.e. h(q) = 0, that the MFCQ holds, and that the system is subjected to p bilateral holonomic constraints f (q) = 0. Let Assumption 1 hold, Φ ∩ Ω = ∅, and that one of the following conditions is satisfied:

(i) T S (q) ∩ ker(∇f (q) T ) ∩ ker(M(q)) = {0}, (ii) T S (q) ∩ ker(∇f (q) T ) ∩ ker(M(q)) = {0} and F (q, q, t) T v > 0 for all v ∈ T S (q) ∩ ker(M(q)) ∩ ker(∇f (q) T ), v = 0.

Then the MLCP in [START_REF] Rockafellar | Variational Analysis[END_REF] has at least one solution.

Another way to treat this extension is to rewrite the bilateral constraints as a set of unilateral constraints: f (q) ≥ 0 and -f (q) ≥ 0. This yields the same result. Uniqueness of solutions of the MLCP also follows as an extension of Proposition 5, where one rewrites the second line in [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] as λ 1 (t)λ 2 (t) ∈ ker(∇h(q)), μ 1 (t)μ 2 (t) ∈ ker(∇f (q)), and the first line in [START_REF] Pons | On Dirac's incomplete analysis of gauge transformations[END_REF] as ∇h(q) T z = 0 and ∇f (q) T z = 0, for all z ∈ ker(M(q)). This gives conditions such that the MLCP in [START_REF] Rockafellar | Variational Analysis[END_REF] has at most one solution x. Consequently, if solvability holds, ∇h(q)λ(t) + ∇f (q)μ(t) is unique.

Remark 5

The contact problem with both bilateral holonomic and unilateral constraints is analyzed in [START_REF] Brogliato | Inertial couplings between unilateral and bilateral holonomic constraints in frictionless Lagrangian systems[END_REF], with full rank M(q) and ∇f (q). Thus the bilateral multiplier μ can be calculated as a function of q, q and λ in order to obtain a modified contact LCP.

Transformation of the MLCP into an LCP

It is also of interest to investigate how the MLCP may be solved. A path is to transform it into an LCP, for which numerical solvers exist. Since M(q) is symmetric positive semidefinite of rank r, there exists a unitary matrix U(q) ∈ R n×n , U(q)U(q) T = I , such that U (q)M(q)U (q) T = Mr (q) 0 0 0 , where M r (q) is symmetric positive definite (it is even diagonal) [START_REF] Bernstein | Matrix, Mathematics. Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF]Corollary 5.4.4]. After the variable change y = U(q)x one may rewrite the MLCP [START_REF] Brogliato | Kinetic quasi-velocities in unilaterally constrained Lagrangian mechanics with impacts and friction[END_REF] as

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩
M r (q)y r + (U (q)F (q, q, t)) r = (U (q)∇h(q)) r λ, (U (q)F (q, q, t)) n-r = (U (q)∇h(q)) n-r λ, 0 ≤ λ ⊥ (U (q)∇h(q)) T r M r (q) -1 (U (q)∇h(q)) r λ -(U (q)∇h(q)) T r M r (q) -1 (U (q)F (q, q, t)) r + (U (q)∇h(q)) T n-r y n-r + H (q, q) ≥ 0.

(

) 35
The LCP matrix is the m × m matrix (U (q)∇h(q)) T r M r (q) -1 (U (q)∇h(q)) r . For a matrix A ∈ R n×m , the submatrix (A) r ∈ R r×m consists of its first r rows. A vector y ∈ R n is also decomposed as y = yr yn-r . If M(q) has full rank r = n then ( 35) becomes an LCP. This transformation may also be used to obtain the inclusion (similar to ( 21)): M r (q)y r + (U (q)F (q, q, t)) r (U (q)F (q, q, t)) n-r ∈ -∂Ψ Φ (y), [START_REF] Scheck | Mechanics. From Newton's Laws to Deterministic Chaos[END_REF] where Φ = {y ∈ R n : U(q) T y ∈ Φ} = U(q) -T (Φ), and the chain rule of convex analysis (see Theorem 1 in Appendix B) has been used to obtain [START_REF] Scheck | Mechanics. From Newton's Laws to Deterministic Chaos[END_REF]. This transformation does not simplify the expressions of the normal and tangent cones; however, it simplifies the expression of ker(M(q)). As we shall see in Sect. 6, the singularity of M(q) implies that y n-r = 0 in [START_REF] Rockafellar | Convex Analysis[END_REF].

The constrained Hamiltonian dynamics

It is of interest to investigate the Hamiltonian dynamics of systems with a singular mass matrix, in particular the transformation from Lagrangian to Hamiltonian dynamics, and vice versa. The theory of singular Lagrangian and Hamiltonian systems has a long history in Physics, where geometrical, coordinate-free analysis are performed [START_REF] Batlle | Equivalence between the Lagrangian and Hamiltonian formalism for constrained systems[END_REF][START_REF] Cariñena | Theory of singular Lagrangians[END_REF][START_REF] Kamimura | Singular Lagrangian and constrained Hamiltonian systems, generalized canonical formalism[END_REF][START_REF] Pons | On Dirac's incomplete analysis of gauge transformations[END_REF]. In particular, equivalence is shown in [START_REF] Batlle | Equivalence between the Lagrangian and Hamiltonian formalism for constrained systems[END_REF]. We adopt here a different path with generalized (local) coordinates q and quasi-Lagrange dynamics, where basic convex analysis allows us to extend the Legendre transformation to the singular case. It is commonly admitted that Lagrangian and Hamiltonian formalisms are equivalent, and the Fenchel (or Legendre-Fenchel) transformation (5) applied to the system's Lagrangian function allows one to pass from one formalism to the other, under the condition of strict convexity of the Lagrangian function with respect to q [6, §15.A] [36, §2.

. When M(q) has full rank n, the Legendre-Fenchel transform (5) of the Lagrangian function L(q, q, λ), defines the Hamiltonian function H (q, p, λ), which is the total mechanical energy of the system. One may write L(q, q, λ) = qT M(q) q -V (q)+ h(q) T λ and H (q, p, λ) = 1 2 p T M(q) -1 p + V (q)-h(q) T λ where V (q) is the potential energy.

Under the full rank condition r = n, one has L (q, q, λ) = H (q, p, λ) and by the Fenchel-Moreau theorem [13, Theorem I.10] L(q, q, λ) = H (q, p, λ), where both conjugacy calculations are made with q and p as variables, respectively. The complementarity Hamiltonian system is given by

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ q = ∂H ∂p = M(q) -1 p, ṗ = -∂H ∂q = -∂ ∂q ( 1 2 p T M(q) -1 p + V (q)) + ∇h(q)λ, 0 ≤ λ ⊥ h(q) ≥ 0. (37)
When r < n, the strict convexity is lost, so that the Lagrangian function is a convex degenerate quadratic function of q. Its Legendre-Fenchel transform L (q, q, λ) may be calculated following the material in [24, Chap. E, Example 1.1.4]. One obtains

H (q, p, λ) = 1 2 p T M(q) † p + V (q) -h(q) T λ, p ∈ Im(M(q)), (38)
where M(q) † is the Moore-Penrose generalized inverse of M(q). Since M(q) † is symmetric positive semi-definite, the Hamiltonian function H (q, p, λ) is convex in p, so its conjugate function is the original Lagrangian function, under the constraint q ∈ Im(M(q) † ). A constraint appears in [START_REF] Schutte | New approach to the modeling of complex multibody dynamical systems[END_REF], reflecting the loss of rank of M(q) that makes the Lagrangian function degenerate, i.e. ∂ 2 ∂ q2 L(q, q, λ) is not of full rank. Consequently, the Lagrange dynamics are an implicit system, see Remark 6. If this constraint is not satisfied then H (q, p, λ) = +∞. The basic idea of what comes next is to use the diagonalization of M(q) as in Sect. 5, so that the kinetic energy (more exactly the kinetic co-energy) expressed in the new variable (a quasi-velocity) has a simple form. To this end, let us introduce the quasi-velocity v = U(q) q = vr vn-r . Then the kinetic energy becomes

T (q, q) = 1 2 qT M(q) q = 1 2 qT U(q) T U (q)M(q)U (q) T U(q) q = 1 2 qT U(q) T M r (q) 0 0 0 U(q) q = 1 2 v T r M r (q)v r = T (q, v r ) (39)
One has M(q) † = U(q) T Mr (q) -1 0 0 0 U(q) [10, Fact 8.15.1 and §6.2]. Thus the kinetic energy is given by

T (q, p) = 1 2 p T M(q) † p = 1 2 p T U(q) T U (q)M(q) † U(q) T Up = 1 2 p T U(q) T M r (q) -1 0 0 0 U(q)p. ( 40
)
Let m = U(q)p. Doing the same partition as for the quasi velocity, one obtains

T (q, p) = 1 2 m T M r (q) -1 0 0 0 m = 1 2 m T r M r (q) -1 m r . (41)
Equaling the kinetic energy in [START_REF] Udwadia | Equations of motion for general constrained systems in Lagrangian mechanics[END_REF] and co-energy in [START_REF] Simeon | Computational Flexible Multibody Dynamics. A Differential-Algebraic Approach[END_REF], one finds that m r = M r (q)v r , which extends the relation between the generalized momentum and the generalized velocity for full rank M(q) [6, p. 65] to quasi momentum m and quasi-velocity v. If U(q) is a Jacobian then v is a generalized velocity and not merely a quasi-velocity. However, usually this is not the case. The quasi-Lagrange dynamics is obtained as follows, starting from (1):

M r (q) vr -M r (q)( U(q)U(q) T v) r + (U (q)F (q, U (q) T v, t)) r = (U (q)∇h(q)) r λ, (U (q)F (q, U (q) T v, t)) n-r = (U (q)∇h(q)) n-r λ.

(

) 42
It is clear that the loss of rank of M(q) imposes a constraint in the system, represented by the last nr lines in [START_REF] Wojtyra | Joint reaction in rigid body mechanisms with dependent constraints[END_REF]. In the unconstrained case, the nonlinear forces have to satisfy (U (q)F (q, U (q) T v, t)) n-r = 0. In the unilaterally constrained case, these forces have to comply with [START_REF] Rockafellar | Convex Analysis[END_REF], where the multiplier λ is a solution of the LCP, and (U (q)F (q, q, t)) n-r = (U (q)∇h(q)) n-r λ with λ ≥ 0.

Remark 6

The structure of the quasi-Lagrange dynamics in [START_REF] Wojtyra | Joint reaction in rigid body mechanisms with dependent constraints[END_REF] differs from other canonical dynamics obtained when the system has full rank inertia matrix, and is subject to holonomic equality constraints. Transformations applied to systems with holonomic constraints usually aim at decoupling the dynamics into a "tangent" (independent of the contact force) and "normal" (contact force dependent) dynamics, using various annihilators of ∇h(q) [START_REF] Brogliato | Inertial couplings between unilateral and bilateral holonomic constraints in frictionless Lagrangian systems[END_REF][START_REF] Brogliato | Kinetic quasi-velocities in unilaterally constrained Lagrangian mechanics with impacts and friction[END_REF][START_REF] Laulusa | Review of classical approaches for constraint enforcement in multibody systems[END_REF][START_REF] Mcclamroch | Feedback stabilization and tracking of constrained robots[END_REF]. Systems with a singular mass matrix are descriptor variable systems, or implicit systems in the quasi-linear form G( ẋ, x, t, λ) = E(x) ẋ + H (x, t, λ), where

E(x) = I 0 0 M(x 1 ) , x 1 = q, x 2 = q, x = x 1
x 2 , and with ∂G ∂ ẋ (x) = E(x) a singular Jacobian matrix (here E(x) has rank 2nr). Transformations aim at decoupling the dynamics according to the structure of E(x) [START_REF] Zimmer | On observing nonlinear descriptor systems[END_REF].

Let us now obtain the quasi-Hamilton dynamics. From the fact that q = ∂H ∂p = ∂p ∂m ∂H ∂m and p = U(q) T m, one obtains ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ q = U(q) T ∂H ∂m , ṁ = -U(q) ∂H ∂q -U(q) U(q)m, m ∈ U(q) Im(M(q)) = Im(U (q)M(q)) = Im(U (q)M(q)U (q) T ) ⇔ m = mr 0n-r . [START_REF] Yen | Constrained equations of motion in multibody dynamics as ODEs on manifolds[END_REF] Notice that if r = n then U(q) = I , and one recovers [START_REF] Shampine | Solving index-I DAEs in Matlab and Simulink[END_REF]. The system in ( 43) is called quasi-Hamiltonian because of the bias -U(q) U(q)m. The first equation in [START_REF] Yen | Constrained equations of motion in multibody dynamics as ODEs on manifolds[END_REF] rewrites v = ∂H ∂m , equivalently v r = M r (q) -1 m r , v n-r = 0. This reveals that r < n imposes the constraint v n-r = 0 in addition to m n-r = 0. One infers that q ∈ U(q) T Im(U (q)M(q)U (q) T ) = Im(M(q)) = Im(M(q) † ), and p ∈ U(q) T Im(U (q)M(q)U (q) T ) = Im(M(q)) = Im(M(q) † ) (as known from [START_REF] Schutte | New approach to the modeling of complex multibody dynamical systems[END_REF], and we used [10, Proposition 6.1.6(vii)]). Therefore, the admissible generalized velocities and momenta are those which correspond to motions that create kinetic co-energy and energy, respectively. Taking into account the velocity constraint allows one to conclude the reverse of (38): L(q, q, λ) = H (q, p, λ).

The generalized momentum constraint of ( 38) is transformed into an explicit quasimomentum constraint in [START_REF] Yen | Constrained equations of motion in multibody dynamics as ODEs on manifolds[END_REF]. Starting from [START_REF] Wojtyra | Joint reaction in rigid body mechanisms with dependent constraints[END_REF] the dynamics for m r is obtained as ṁr = Ṁr (q)M r (q) -1 m r + M r (q) U(q)U(q) T v -U(q)F q, U (q) T v, t r + U(q)∇h(q) r λ.

(

) 44
The last nr lines of the ṁ dynamics in [START_REF] Yen | Constrained equations of motion in multibody dynamics as ODEs on manifolds[END_REF] correspond to the same equality as in the second line of [START_REF] Wojtyra | Joint reaction in rigid body mechanisms with dependent constraints[END_REF], due to the constraint on m.

Remark 7 If r = 0 then M(q) = 0 and both q and p are identically zero, as well as their derivatives q and ṗ. The system is in an equilibrium state.

Let us recapitulate the developments of this section. The next dynamical equations are equivalent one to another, where U(q) is a unitary matrix such that U (q)M(q)U (q) T = diag(M r (q), 0), rank(M r (q)) = r < n:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ M(q) q + F (q, q, t) = ∇h(q)λ, 0 ≤ λ ⊥ h(q) ≥ 0, q ∈ Im(M(q)), L(q, q, λ) = 1 2 qT M(q) q -V (q) + h(q) T λ ⇔ ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ q = ∂H ∂p , ṗ = -∂H ∂q , 0 ≤ λ ⊥ h(q) ≥ 0, p ∈ Im(M(q)) H (q, p, λ) = 1 2 p T M(q) † p + V (q) -h(q) T λ (45) ⇔ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ m = U(q)p, m n-r = 0, ṁr = Ṁr (q)M r (q) -1 m r + M r (q)( U (q)U (q) T v) -(U (q)F (q, U (q) T v, t)) r
+ (U (q)∇h(q)) r λ, (U (q)F (q, U (q) T v, t)) n-r = (U (q)∇h(q)) n-r λ,

0 ≤ λ ⊥ h(q) ≥ 0 (46) ⇔ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ m = U(q)p, m n-r = 0, v= U(q) q, v n-r = 0, H (q, m, λ) = 1 2 m T r M r (q) -1 m r + V (q) -h(q) T λ, v = ∂H ∂m , ṁ = -U(q) ∂H ∂q -U(q) U(q)m, 0 ≤ λ ⊥ h(q) ≥ 0 (47) ⇔ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ v = U(q) q, v n-r = 0, L(q, v, λ) = 1 2 v T r M r (q)v r -V (q) + h(q) T λ, M r (q) vr -M r (q)( U(q)U(q) T v) r + (U (q)F (q, U (q) T v, t)) r = (U (q)∇h(q)) r λ, (U (q)F (q, U (q) T v, t)) n-r = (U (q)∇h(q)) n-r λ, 0 ≤ λ ⊥ h(q) ≥ 0. (48)
The LCP in [START_REF] Rockafellar | Convex Analysis[END_REF] may be obtained as the contact LCP of (48). The constraints p n-r = 0 are called the primary constraints in the literature [START_REF] Cariñena | Theory of singular Lagrangians[END_REF][START_REF] Pons | On Dirac's incomplete analysis of gauge transformations[END_REF]. Finally, the classical equality H (q, p, λ) = p T q -L(q, q, λ) remains valid, and equivalently H (q, m, λ) = m T v -L(q, v, λ).

Remark 8 There is nothing fundamentally new in the above developments. They are, however, tailored to multibody systems where the usefulness of coordinate-free, purely geometrical arguments as in [START_REF] Batlle | Equivalence between the Lagrangian and Hamiltonian formalism for constrained systems[END_REF][START_REF] Cariñena | Theory of singular Lagrangians[END_REF][START_REF] Kamimura | Singular Lagrangian and constrained Hamiltonian systems, generalized canonical formalism[END_REF][START_REF] Pons | On Dirac's incomplete analysis of gauge transformations[END_REF] stemming from Dirac's original results is unclear when time comes to study contact complementarity problems.

Examples

Many examples of systems with a singular mass matrix are given in [3-5, 25, 38-41], besides the well-known case when Euler angles are used and the mass matrix of a rigid body is obtained from the Olindes-Rodrigues matrix, rendering it singular at certain configurations [14, Eq. (4.3)], [25, §6]. To illustrate the foregoing developments let us present two simple examples.

Example 1

Let us consider the pendulum dynamics, taken from [1, § 3.2.3], as depicted in Fig. 5. Let the coordinates of G be x and y in the Galilean frame R 0 = (i, j), and those of O be x o and y o . The kinetic energy of the free pendulum is T ( θ, ẋ, ẏ) = 1 2 m( ẋ2 + ẏ2 ), as indeed the moment of inertia at G is zero, all the mass being concentrated at G. The potential energy is V (q) = Fig. 5 A simple pendulum mgx. The Lagrangian function is L(θ, θ, x, y, ẋ, ẏ) = T ( θ, ẋ, ẏ) -V (x). Since O is fixed, there are two bilateral constraints x O = xl cos θ = 0, y O = yl sin θ = 0, with which two Lagrange multipliers μ 1 and μ 2 are respectively associated. Let the unilateral constraint consist of a vertical wall passing through the point O, imposing the constraints h 1 (q) = y ≥ 0, h 2 (q) = θ ≥ 0, with which two Lagrange multipliers λ 1 and λ 2 are respectively associated. One has λ = (λ 1 , λ 2 ) T and μ = (μ 1 , μ 2 ) T . Thus

∇h(q) = ⎛ ⎝ 0 0 1 0 0 1 ⎞ ⎠
has rank 2, and the MFCQ holds. Let us choose the generalized coordinates of the free system as (x, y, θ) T . We can write the Lagrange equations for the pendulum evolving in the plane and subject to the bilateral and unilateral constraints as

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ m ẍ(t) = μ 1 (t) + mg + F 1 (t), m ÿ(t) = μ 2 (t) + λ 1 (t) + F 2 (t), 0. θ(t) = l sin θμ 1 (t) -l cos θμ 2 (t) + λ 2 (t), f 1 (q) = x -l cos θ = 0, f 2 (q) = y -l sin θ = 0, 0 ≤ λ 1 ⊥ h 1 (q) = y ≥ 0, 0 ≤ λ 2 ⊥ h 2 (q) = θ ≥ 0,
where F = (F 1 , F 2 ) T is the force applied at G, expressed in R 0 . Here M(q) = M = diag(m, m, 0), so One obtains T S (q) ∩ ker(M) ∩ ker(∇f (q) T ) = {0}, hence Proposition 6(i) applies for the MLCP in [START_REF] Rockafellar | Variational Analysis[END_REF] which is solvable. The uniqueness of λ follows from ker(∇h(q)) = {0}, and the uniqueness of μ follows from ker(∇f (q)) = {0}. One may also consider the unilateral constraint h 3 (q) = -θ + π ≥ 0. However, this does not change the analysis because the unilateral constraints 2 and 3 cannot be active at the same time. Thus the MFCQ still holds despite ∇h 2 (q) and ∇h 3 (q) are colinear.

Example 2

Let us now consider the system in Fig. 6, taken from [40, §5 (c)]. It consists of two masses m 1 > 0 and m 2 > 0 connected by linear springs. The unstretched lengths of the springs are l 10 and l 20 , and their stiffnesses are k 1 > 0 and k 2 > 0, respectively. According to the notations shown on the figure, one chooses q = ( x1 , q 1 , q2 ) T , with x1 = x 1l 10 , q2 = q 2l 20 . One analyzes the system as the interconnection of two subsystems, which are interconnected with a bilateral constraint f (q) = 0. The unilateral constraint is imposed on mass m 1 by a rigid barrier at a distance L from the wall. Then one obtains the constrained dynamics:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎛ ⎝ m 1 0 0 0 m 2 m 2 0 m 2 m 2 ⎞ ⎠ q(t) = ⎛ ⎝ -k 1 x1 (t) 0 -k 2 q2 (t) ⎞ ⎠ + ⎛ ⎝ -μ(t) -λ(t) μ(t) 0 ⎞ ⎠ , f (q) = q 1 -x1 -l 10 -d = 0, 0 ≤ λ ⊥ h(q) = L -x1 -l 10 -d ≥ 0. ( 49
)
One has ker(∇f (q) T ) = {z ∈ R 3 : z 1 = z 2 }, T S (q) = {z ∈ R 3 : -z 1 ≥ 0}, ker(M) = {z ∈ R 3 : z 1 = 0, z 2 = -z 3 } and obtains T S (q) ∩ ker(M) ∩ ker(∇f (q) T ) = {0}, hence Proposition 6(i) applies for the solvability of the MLCP in [START_REF] Rockafellar | Variational Analysis[END_REF]. The uniqueness of λ and μ follows from the fact that ker(∇h(q)) = {0} and ker(∇f (q)) = {0}.

Remark 9 It is noteworthy that both examples could be analyzed with a set of independent generalized coordinates, by eliminating the redundant coordinates due to the presence of the bilateral constraints. The singularity of the mass matrix comes from the generalized coordinates redundancy. In practice one often chooses to work with redundant coordinates, assembling subsystems by adding bilateral constraints. is also often hard, if not impossible, to eliminate redundant constraints.

Conclusions

The main objective of this article is to analyze the well-posedness of the complementarity problem which arises in unilaterally constrained systems, when the mass matrix and the constraints Jacobian matrix are not of full rank. This gives rise to a mixed linear complementarity problem. Conditions guaranteeing that the acceleration and the contact force Lagrange multiplier exist and are unique are stated. The case when bilateral, holonomic constraints are also present is also analyzed. Moreover, the relationships between the (singular) Lagrangian and the Hamiltonian formalisms are studied.

Fig. 1

 1 Fig. 1 Recession cones

⎞⎠

  , z 1 ∈ R, z 2 ∈ R ⎫ ⎬ ⎭ .

Fig. 6 -l sin θ l cos θ 1 ⎞⎠ z 3 , z 3

 6133 Fig.6Flexible two-degree-of-freedom system (from[START_REF] Udwadia | Explicit equations of motion for constrained mechanical systems with singular mass matrices and applications to multi-body dynamics[END_REF] Fig. 3])

The contact complementarity problemLet us consider[START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF] subjected to a set of m unilateral constraints h(q) ≥ 0, giving rise to complementarity conditions 0 ≤ λ ⊥ h(q) ≥ 0. During persistent contact phases with the m constraint surfaces, the complementarity can be written at the acceleration level as 0 ≤ λ(t) ⊥ d 2 dt 2 h(q(t)) ≥ 0 (if not all m constraints are active, then one constructs the LCP taking

Appendix A: The Mangasarian-Fromovitz Constraint Qualifications (MFCQ)

[19, pp. 17, 252]

, where h : R n → R m is continuously differentiable, and K is not necessarily convex. Suppose that there exists a vector v ∈ R n such that ∇h i (x) T v > 0 for all i ∈ I(x) = {i : h i (x) = 0}. Then the tangent cone to K at x, defined as the dual of the normal cone in ( 14) is equal to the linearization cone {z ∈ R n : z T ∇h i (x) ≥ 0, for all i ∈ I(x)}.

Appendix B: The chain rule of convex analysis

Theorem 1 [


Olivier Agullo

W.-C Müller

B Knaepen

D Carati

W.-C Mu

) B Knaepen

Large eddy

The numerical large eddy simulation ͑LES͒ technique is tested on decaying magnetohydrodynamic ͑MHD͒ turbulence. The LES approach allows for a strong reduction in computational cost compared to direct numerical simulations by modeling the effects of the smallest turbulent scales instead of computing them directly. Two small-scale models of eddy-viscosity type are presented for this purpose in combination with a procedure for the self-consistent calculation of the model parameters in the course of the simulation. The method is successfully tested by comparing the obtained results to a high-resolution direct numerical simulation of decaying three-dimensional MHD turbulence.

Many plasmas in turbulent motion observed in astrophysical systems, as well as nuclear fusion devices, can be described within the framework of incompressible magnetohydrodynamics ͑MHD͒. Since direct experimental access to fully developed MHD turbulence is often difficult or even impossible, one resorts to the numerical simulation of idealized systems by solving the incompressible MHD equations ץ t vϭϪٌ•͑vvϪbb͒ϩ⌬vϪٌ p, ͑1͒

ץ t bϭϪٌ•͑vbϪbv͒ϩ⌬b, ͑2͒ ٌ•vϭٌ•bϭ0, ͑3͒

written here with the fluid velocity v, the magnetic-field b expressed in Alfve ´n speed units ͑the constant mass density is rescaled to unity͒, the kinematic viscosity , and magnetic diffusivity . The dynamic pressure p is obtained by imposing the incompressibility of v.

Unfortunately, the numerical resolution necessary to render the discrete form of a turbulent MHD flow into a physically meaningful representation of the real system increases rapidly with the magnetic Reynolds number Rmϭv 0 l 0 /, v 0 and l 0 being a characteristic flow velocity and length scale, respectively. This limits direct numerical simulations ͑DNS͒ of three-dimensional MHD turbulence on today's supercomputers [1][2][3] to values of Rm that are still several orders of magnitude lower than the ones found in nature, e.g., in the solar convection zone. A possible alternative to DNS is the large-eddy simulation ͑LES͒ approach, where the Reynolds number gap between simulation and reality is narrowed by restricting direct numerical computation to the large spatial scales of turbulence in combination with a model to incorporate the effects of the small-scale motions. This approximation is based on the assumption that the large-scale dynamics reflects the characteristic physics of a specific turbulent system, e.g., the influence of geometry and boundary conditions, while the small-scale turbulent fluctuations act as a statistically quasi-uniform background.

The LES technique is considered a valuable tool for simulating rather large Reynolds number Navier-Stokes turbulence for more than thirty years now, 4 while work on adapting the method to MHD turbulence [5][6][7][8] has not received as much attention. This might in part be caused by the fact that LES small-scale models usually contain at least one free parameter, which has to be tuned to obtain optimal performance of the simulation, turning a LES into a trial and error process. However, not long ago a dynamic procedure was developed 9 that allows for a self-consistent calibration of the small-scale model during the simulation, leading to convincing results when used in Navier-Stokes turbulence. In this Brief Communication we present the generalization and application of this procedure to isotropic three-dimensional MHD turbulence, applying two generalized eddy-viscosity models to describe the small-scale dynamics. Scalar-level tests are performed by comparing the observed temporal and spectral behavior of the kinetic and magnetic energies to the results of a high-resolution DNS of decaying MHD turbulence. 3 The formal basis of LES is the application of a spatial gridfilter of width l ¯͑Ref. 10͒ to ͑1͒-͑3͒. where two unknown terms, usually referred to as subgridscale or filtered-scale stress tensors, have been introduced by the filtering operation and need to be modeled: v ϭ(vv Ϫv ¯v ¯)Ϫ(bbϪb ¯b ¯) and b ϭ(vbϪv ¯b ¯)Ϫ(bvϪb ¯v ¯). These tensors account for effects of the subgrid-scale motions on the flow at the resolved scales larger than l ¯. The filter chosen here is a sharp cubic Fourier cutoff, which sets all spatial Fourier modes with at least one wave-vector component larger than k c ϭ/ l ¯to zero. In order to assess the LES technique and the dynamic procedure in MHD we have performed simulations of decaying isotropic turbulence using the same pseudospectral method as the reference DNS, the same integration domain ͑a cube of linear extension 2 with periodic boundary con-ditions͒ and the same set of parameters (ϭϭ3•10 Ϫ4 ) at a significantly reduced numerical resolution. The initial condition has been generated from the reference DNS data at t ϭ1 when kinetic and magnetic energy dissipation, starting from smooth initial fields, have reached their maxima and the dissipative small-scale structures of turbulence are fully developed. The data set has been cutoff-filtered from 512 3 to 64 3 Fourier modes, removing about 99.8% of the originally available information. The influence of the filtered scales on the remaining large scales of motion, which still contain 90% of the total energy, has to be mimicked by the models for the filtered-scale stress tensors v and b .

The model expressions used here are based on the eddyviscosity assumption widely applied in LES of Navier-Stokes turbulence. 11 This class of models aims at reproducing the subgrid-scale kinetic-energy dissipation by assuming a linear relation between the deformation tensors at resolved and dissipative scales

v ϭϪ2 t S ¯, S ¯ϭٌ͑v ¯ϩٌ͓v ¯͔T ͒/2,
where the superscript ''T'' denotes the transposed matrix and t is the eddy-viscosity. When restricting the dependence of t to a characteristic mixing-length l ¯on the spatial scales under consideration and the subgrid kinetic-energy dissipation K , dimensional analysis suggests t ϳ l ¯4/3 ( K ) 1/3 . For the magnetic energy one can proceed analogously, taking into account that the property of the magnetic field, which gives rise to resistive dissipation, is its curl rather than its deformation. This yields the model prototype

b ϭϪ2 t J ¯, J ¯ϭٌ͑b ¯Ϫٌ͓b ¯͔T ͒/2,
with the eddy-resistivity t ϳ l ¯4/3 ( M ) 1/3 , M denoting the subgrid magnetic energy dissipation.

Depending on the way of estimating the subgrid dissipation, different eddy-viscosity models can be constructed. Here we consider two variants. If l ¯lies somewhere in the inertial scale range of fully developed turbulence and, as observed in the reference DNS, the nonlinear energy transfer between kinetic and magnetic energy is much smaller than the respective energy dissipation, K and M can be assumed as spatially constant. This yields the generalized Kolmogorov scaling subgrid model

M 1 : t ϭC 1 ͑ t ͒ l ¯Ϫ4/3 , t ϭD 1 ͑ t ͒ l ¯Ϫ4/3 ,
where C 1 and D 1 are unknown parameter functions of time.

The widely applied Smagorinsky model 12 approximates K as the product of the resolved local dissipation rate and the associated length scale l ¯, leading to K ϳ l ¯2(2S:S) 3/2 . Together with the corresponding extension to MHD 8 M ϳ l ¯2͉j ¯͉3 , j ¯ϭٌϫb ¯being the resolved electric current density, this gives

M 2 : t ϭC 2 ͑ t ͒ l ¯2͑2S:S͒ 1/2 , t ϭD 2 ͑ t ͒ l ¯2͉j ¯͉.
Both models contain parameters C i and D i (iϭ1,2), for which no direct relationship with the resolved large-scale quantities can be obtained without further assumptions. Therefore, one would be forced to prescribe these parameters in order to optimize the simulation results with respect to known reference characteristics. Apart from the arbitrariness introduced to the subgrid-modeling by such an approach, the lack of sufficiently large experimental MHD turbulence databases, like they exist for the Navier-Stokes case, prevents any a priori justification of freely prescribed subgrid-model parameters. Fortunately, this problem can be circumvented by exploiting the self-similarity of fully developed turbulence. A corresponding technique, that has proven to be successful in Navier-Stokes LES, is the dynamic procedure. [13][14][START_REF] Moin | AIAA 24th Fluid Dynamics Conference[END_REF][START_REF] Carati | [END_REF] The basic idea of this method is the application of a second filter-the test filter-of width l ˆϾ l ¯to the LES equations ͑4͒-͑6͒. The test filtering operation, here a cubic cutoff with l ˆϭ2 l ¯, will be indicated by the hat symbol, • ˆ. Introducing the filtered-scale stress tensors at test filter level as

T v ϭ(vv RϪv Rv RϪ(bb R Ϫb R b R ) and T b ϭ(vb R Ϫv Rb R )Ϫ(bv R Ϫb R v R) leads to the Germano identities 14 for MHD T v ϭ ˆvϩL v , T b ϭ ˆbϩL b , ͑7͒
where the Leonard tensors

L v ϭ(v ¯v ¯ Ϫv Rv R)Ϫ(bb R Ϫb R b R ) and L b ϭ(vb R Ϫv Rb R )Ϫ(bv R Ϫb R v R)
are known in terms of the resolved fields and, as such, do not require any modeling. In contrast, the filtered-scale stress tensors at grid scales ( v , b ) and test scales (T v ,T b ) have to be modeled. Since the turbulent fields are assumed to be self-similar on scales within the interval ͓ l ¯, l ˆ͔, the applied models and their parameters are independent of the test filter width. Hence, for the stress tensors T v and T b at test filter level the same model and coefficients should be used as for v and b , respectively. tered equations can then be solved using a coarser grid, since fluctuations on spatial scales smaller than ¯l -the subgrid scales-have been removed. In the following, ¯l -filtered quantities will be denoted by an overbar. The filtered equations for v and ¯b thus read Due to the modeling of the filtered-scale stress tensors the identities ͑7͒ are not exactly fulfilled. The errors can be quantified as

Q v ϭ ͵ V d 3 r͑T v Ϫ ˆvϪL v ͒ 2 , Q b ϭ ͵ V d 3 r͑T b Ϫ ˆbϪL b ͒ 2 ,
where the model parameters are assumed to be constant in space, which is justified when dealing with spatially homogeneous turbulence. Minimizing the deviations Q v and Q b by variation with respect to the model coefficients 13 yields an optimization ansatz for the filtered scale stress tensors. The implementation of the dynamic procedure is somewhat more difficult when there is no flow direction of homogeneity. However, several solutions have been considered in that case, like the local dynamic procedure 17 or the Lagrangian dynamic procedure. 18 The choice of the ratio ␣ϭ l ˆ/ l ¯ϭ2 as mentioned above has been guided by a series of test runs with model M 1 , confirming the self-similarity of the turbulent fields for ␣͓1.5,2.5͔.

A basic LES requirement is the reproduction of the temporal and spectral behavior of the filtered kinetic (E K ) and magnetic (E M ) energy. The time evolution of these quantities is shown in Figs. 1 and2, where t is given, like for the reference DNS, in units of the large eddy turnover time. Diamonds represent the filtered DNS results, the solid and the dashed lines give the energy development for LES with the models M 1 and M 2 , respectively. To evaluate the overall influence of the dynamic subgrid modeling the dotted lines show the energy curves for a LES with M 0 : v ϭ b ϭ0.

As expected, the M 0 -LES results strongly deviate from the filtered DNS data, since energy dissipation is mainly due to the high wave number modes, which have been cutoff by the gridfiltering operation. The lack of a well-developed dissipation range hampers energy decay and leads to an accumulation of kinetic and magnetic energy in the high wave number region of the LES system ͑see below͒.

A clear improvement is observed when the subgrid models M 1 and M 2 are applied. The evolution of E K is well reproduced in both cases, showing that the main influence of the small-scale velocity field fluctuations on the kinetic energy is dissipative. The temporal development of the magnetic energy in the LES with M 1 and M 2 also satisfactorily fits the reference data, though one observes an offset between the LES results and the DNS. Closer inspection reveals that this phenomenon is caused by an initial overshoot of the model resistivities estimated by the dynamic procedure. This effect is known from Navier-Stokes LES using dynamic models and can be reduced by allowing for a short model relaxation period at the beginning of the simulation. Due to the spatial homogeneity of the turbulent system, there is only a small systematic difference between the results obtained with the Kolmogorov scaling model M 1 , which mainly depends on the spatially averaged extrapolation process inherent to the dynamic procedure, and the Smagorinsky version M 2 , that additionally takes into account the local field gradients.

The application of the virtually parameter free dynamic subgrid models reproduces rather sensitive quantities like the angle-averaged energy spectra in good agreement with the DNS data. The kinetic-energy spectra E k K are shown in Fig. 3 at time tϭ6, when E K has decreased by a factor of about 6.5. Both models, M 1 and M 2 , lead to spectra that follow the filtered DNS data up to the high wave number range of the LES system. Nevertheless, the simple applied eddy-viscosity models provide too much dissipation in the high wave number range, leading to a slightly stronger fall-off of the LES spectra when compared to the reference data. One observes that the over-dissipation is reduced by the Smagorinsky term in model M 2 , yielding a better spectral performance.

The agreement of the LES results with the filtered DNS spectrum is evidently due to the subgrid models M 1 and M 2 , since the lack of such a model in the M 0 -LES causes a large accumulation of kinetic energy over nearly two-thirds of the LES wave number space, which is the consequence of the missing filtered-scale energy dissipation in combination with the direct spectral energy cascade.

The same trends are observed in Fig. 4 for the angleaveraged magnetic energy spectrum E k M . However, the LES show a spectral magnetic-energy distribution that is too low across a wide range of scales when compared to the DNS spectrum. Since the M 0 -LES without model exhibits the same behavior, the observation cannot be explained by the subgrid models being too dissipative. This implies that the model for the filtered-scale stress tensor b in Eq. ͑5͒ has to be more complex than the simple dissipative subgrid models M 1 and M 2 , e.g., allowing for a spectrally local inflow of magnetic energy. Nevertheless, the improvement in reproducing the magnetic energy reference spectrum using the presented subgrid-scale models compared to a M 0 -LES is still tremendous.

The results presented in this Brief Communication show that LES with dynamic subgrid modeling are able to reproduce the temporal and spectral energy evolution in isotropic MHD turbulence. Since decaying turbulence was considered, the LES subgrid models had to adapt to a nonstationary environment, where a continuous self-calibration of the model parameters as realized in the dynamic procedure is of particular importance. The achieved LES performance, especially when the highly reduced numerical resolution and the simplicity of the purely dissipative filtered-scale models are taken into account, encourages further work on dynamic subgrid modeling of MHD turbulence.

FIG. 1 .

 1 FIG. 1. Decay of the kinetic-energy E K . Diamonds correspond to the reference DNS filtered to 64 3 , the dotted line shows the result of the M 0 -LES without model. The solid and dashed lines represent the predictions of the LES with, respectively, model M 1 and model M 2 .

ACKNOWLEDGMENTS D.C. and B.K. are, respectively, ''Chercheur Qualifie ´'' and ''Charge ´de recherche'' of the ''Fonds National de la Recherche Scientifique'' ͑Belgium͒.