A recursive multibody formalism for systems with small mass and inertia terms

Abstract. Complex multibody system models that contain bodies with small mass or nearly singular inertia tensor may suffer from high frequency solution components that deteriorate the solver efficiency in time integration. Singular perturbation theory suggests to neglect these small mass and inertia terms to allow a more efficient computation of the smooth solution components. In the present paper, a recursive multibody formalism is developed to evaluate the equations of motion for a tree structured N body system with O(N) complexity even if isolated bodies have a rank-deficient body mass matrix. The approach is illustrated by some academic test problems in 2-D.


Introduction
Classical time integration methods in technical simulation are tailored to problems with smooth solution.Small system parameters in a mathematical model may introduce rapidly oscillating or strongly damped solution components that cause problems in time integration.Singular perturbation theory gives much insight in the analytical background of these phenomena and allows furthermore an efficient approximation of smooth solutions neglecting all terms that contain small parameters, see, e.g., Hairer and Wanner (1996).
The application of these classical results to multibody dynamics is non-trivial since the numerical algorithms for evaluating the equations of motion efficiently (multibody formalisms) are based on regularity assumptions that may be violated if small mass and inertia terms are neglected.A modified multibody formalism for chain structured multibody systems with an isolated "zero mass" body was developed in Arnold et al. (2010).
The present paper is a revised and extended version of the author's contribution to this conference paper.A recursive multibody formalism is developed for tree structured systems with bodies that suffer from a rank-deficient body mass matrix and may be considered as limit case of systems with bodies of (very) small mass or nearly singular inertia tensor.This research is guided by well known results from general singu-lar perturbation theory, see Hairer and Wanner (1996), and its extensions to singularly perturbed problems in the context of multibody dynamics by Lubich (1993) and Stumpp (2008).
Related problems are, e.g., the modelling of serial springdamper elements using an auxiliary zero mass body between spring and damper (Eich-Soellner and Führer, 1998, Section 1.3.4), the modification of inertia forces for highfrequency eigenmodes of flexible bodies in multibody system models for the analysis of elastohydrodynamic bearing coupling in Schönen (2003) and recently proposed methods from FE contact mechanics in Hager and Wohlmuth (2009).
For real-time applications in multibody dynamics, the neglection of inertia forces for small mass bodies was studied by Eichberger and Rulka (2004).A more detailed analysis shows, that this neglection of inertia forces is straightforward if all small mass bodies of the multibody system are leaf bodies in the kinematic tree.In numerical experiments for the model of a walking mobile robot (mobot) with stiff contact forces between lightweight legs and ground floor, the numerical effort was reduced by a factor of 4, see Weber et al. (2012).
The singular perturbation analysis is technically more complicated for multibody system models with small mass bodies having in the kinematic tree a successor of substantially larger mass since classical multibody formalisms and topological solvers are not applicable to the limit case of zero Published by Copernicus Publications.
mass bodies in a kinematic chain.In Burgermeister et al. (2011), a smoothed velocity approximation for small mass bodies was proposed as a work-around.
In the present paper, we discuss an alternative approach that extends the recursive multibody formalism directly to kinematic trees containing bodies of vanishing mass or rankdeficient inertia tensor.The classical set of second order equations of motion in the joint coordinates q(t) is substituted by a suitable combination of second and first order ODEs describing the system dynamics in the limit case of a zero mass body.These results extend the previous analysis for chain structured systems in Arnold et al. (2010) and may be considered as a next step to extend advanced multibody formalisms for flexible multibody systems to models with bodies of (very) small mass.
The remaining part of the paper is organized as follows: Basic results of singular perturbation theory and its application to multibody system dynamics are recalled in Sect. 2. A recursive multibody formalism and the resulting mixed set of first and second order equations of motion for a body of rankdeficient body mass matrix are derived in Sect.3. Finally, two simple test problems are discussed in Sect. 4. Some basic information on Moore-Penrose pseudo-inverses is provided in Appendix A.

Singular perturbations in multibody system models
After a short introduction to singularly perturbed ODEs we consider in the present section singular perturbations in multibody system models that are caused by stiff potential forces, see Sect.2.1.Furthermore, some typical problems in the dynamical simulation of multibody system models with small mass bodies are illustrated by the analysis of two coupled oscillators in Sect.2.2.

Time integration of singularly perturbed problems in multibody dynamics
The generic form of singularly perturbed ODEs are partitioned systems with a small perturbation parameter ε > 0 that are considered at a finite time interval [0, t e ], see Hairer and Wanner (1996, Chapter VI) and the references therein.
The general solution of the singularly perturbed problem (1) has the form with smooth functions η ε (t/ε), ζ ε (t/ε) that decay like exp(−βt/ε) for some positive constant β ∈ (0, β), see Hairer and Wanner (1996, Theorem VI.3.2).For many stiff integrators, the numerical solution of (1) may be decomposed as well in a smooth part and a rapidly decaying part reflecting the transient behaviour for initial values y ε (0), z ε (0) that do not belong to a smooth solution.If there is no particular interest in this transient phase, an approximate numerical solution may be obtained much more efficiently solving for given initial values y 0 (0) := y ε (0) the DAE (2) by appropriate time integration methods (Hairer and Wanner, 1996, Chapter VI).Note, that the initial values z 0 (0) in DAE (2) are not free but have to satisfy the consistency condition γ(y 0 (0), z 0 (0)) = 0. Lubich (1993) extended these classical results to a class of singular singularly perturbed problems M(q) q = ψ(q, q) − 1 ε 2 ∇U(q) (4) with q(t) denoting the position coordinates of a multibody system.Matrix M(q) is the symmetric positive definite mass matrix and ψ(q, q) denotes a vector of forces and momenta.The crucial term in (4) are the (very) stiff potential forces −ε −2 ∇U(q) that depend on the perturbation parameter ε with 0 < ε 1.If the potential U(q) attains a local minimum on a manifold U and is strongly convex along all directions that are non-tangential to U, then the smooth solution of (4) may be approximated up to O(ε 2 ) by the solution q 0 (t) of a constrained (differential-algebraic) system with q 0 (t) ∈ U, ( t ≥ 0 ).In general, this constrained system may be solved much more efficiently than the original singularly perturbed problem (4), see Hairer and Wanner (1996).
In Stumpp (2008), these results were extended to mechanical systems with strong damping forces −ε −1 D(q) q in the right hand side of (4).In the limit case ε → 0, the solution of the singularly perturbed problem is again approximated by the solution of a DAE problem that can be obtained in a robust and efficient way by standard DAE time integration methods, see Hairer and Wanner (1996).
Example problem: two coupled oscillators with a fast oscillating small mass m 1 , see Burgermeister et al. (2011).

The small mass oscillator as singularly perturbed problem
The eigenfrequency of a harmonic oscillator is given by ω = √ c/m with mass m and spring constant c. High frequency oscillations in a mechanical system may not only be introduced by (very) stiff potential forces but also by potential forces of moderate size that act on a body with (very) small mass.
In Burgermeister et al. (2011), this phenomenon was studied for the simple model problem in Fig. 1.In two coupled oscillators, a small mass m 1 is connected to a large mass m 2 and the reference system by stiff springs with constants c 1 , c 2 and damping with damping parameters d 1 , d 2 .Both bodies can only move along the x-axis.Additional forces F(t) are only acting on m 1 .In absolute coordinates p(t) = (p 1 (t), p 2 (t)) , the equations of motion are given by The small mass can oscillate very fast depending on the ratio of the masses and spring parameters.If the perturbation parameter ε := m 1 gets smaller, the frequency of the oscillations increases and a time integration method with stepsize control would choose very small stepsizes to resolve these oscillations and meet the integration tolerances.The number of time steps and the computing time increase significantly, see Fig. 2. In the limit case ε = m 1 = 0, the inertia forces of the first mass point are neglected and (5a) results in an implicit first order differential equation The fast oscillations of the small mass disappear and the integrator can use large stepsizes, see Fig. 2.
In a system description by relative coordinates q(t) = (q 1 (t), q 2 (t)) with p 1 (t) =: q 1 (t) and p 2 (t) − p 1 (t) =: q 2 (t), the limit process ε → 0 causes substantial problems since the equations of motion (5b) of the large mass depend on p2 = q1 + q2 but q1 = p1 does not appear in (6).Furthermore, the differentiation of (6) w.r.t.time t shows that q1 (t) depends on the time derivative of F(t) in the limit case ε = m 1 = 0. 3 Mixed coordinate formulation of the equations of motion (Very) small masses or (nearly) singular inertia tensors in a multibody system model may be interpreted as small perturbations.The analysis of Sect.2.1 suggests to study a reduced system that neglects these perturbations.Therefore, we consider now the limit case of multibody systems that have one or more bodies with zero mass or singular inertia tensor, i.e., with a rank-deficient body mass matrix.
In that case, classical recursive multibody formalisms ("O(N)-formalisms") like the ones by Brandl et al. (1988), Lubich et al. (1992) and Eichberger (1994) fail since they use the inverse of projected body mass matrices.Udwadia and Phohomsiri (2006) consider multibody systems with singular mass matrix but do not exploit the system's topology to evaluate the equations of motion efficiently by a sequence of forward and backward recursions in the kinematic tree.
In this section, we describe the multibody system by the absolute position and orientation of all N bodies relative to the inertial frame and use joint coordinates as generalized coordinates for a tree structured system ( "mixed coordinates"), see also Schiehlen (1997).As in Lubich et al. (1992) and Eich-Soellner and Führer (1998), the recursive multibody formalism is interpreted as a block Gaussian elimination for a sparse, (very) large block-structured system of linear equations with the block structure being determined by the topology of the multibody system.
Using generalized inverses, the recursive formalism of Lubich et al. (1992) is modified to skip bodies with rankdeficient body mass matrix in the second forward recursion, see Sects.3.2 and 3.3.The resulting equations of motion are not longer explicit but form a (linearly) implicit set of first and second order differential equations (Sect.3.4).Finally, we show in Sect.3.5 by some analytical transformations that the equations of motion form a first order index-1 DAE if a certain regularity assumption is satisfied.

Tree structured multibody systems: kinematics
Recursive multibody formalisms are tailored to tree structured systems.Here, the term tree structure corresponds to the structure of the labelled graph being associated to the multibody system model.In this graph, each (rigid or flexible) body of the system is represented by a vertex.Two vertices of the graph are connected by an edge if and only if the corresponding bodies in the multibody system model are connected by a joint restricting their relative motion.
The graph of a tree structured multibody system is acyclic, i.e., it is free of loops.Then, the multibody system has a root body (•) (0) that corresponds to the root vertex of the tree structured graph and is supposed to be inertially fixed.All other bodies (•) (i) have a uniquely defined predecessor (•) (π i ) in the kinematic tree.Each body (•) (i) may have successors (•) ( j) being characterized by π j = i or, equivalently, by j ∈ I i := { k : π k = i } with an index set I i that represents the set of all successors of a given body (•) (i) in the multibody system model.Bodies without successors (I i = ∅) correspond to leafs of the kinematic tree and are therefore called "leaf bodies".
We suppose that position and orientation of body (•) (i) may be characterized by (absolute) position coordinates p i (t) ∈ R d with d = 6 for 3-D models and d = 3 for 2-D models (the position of point masses may be characterized by d = 3 absolute coordinates in 3-D and by d = 2 absolute coordinates in 2-D, see Sect. 4 below).The relative position and orientation of body (•) (i) w.r.t.its predecessor (•) (π i ) is characterized by joint coordinates q i (t) ∈ R n i representing the n i degrees of freedom of the joint connecting (•) (i) with (•) (π i ) : Here and in the following we suppose that ( 7) is locally uniquely solvable w.r.t.p i and that the Jacobian K i = ∂k i /∂p i is non-singular along the solution.In its most simple form, Eq. ( 7) defines p i explicitly by The kinematic relations (7) at the level of position coordinates imply relations at the level of velocity and acceleration coordinates that may formally be obtained by (total) differentiation of (7) w.r.t.time t: with It is supposed that the joint coordinates q i (t) are defined such that all Jacobians J i have full column rank: rank Functions k (I)  i := ∂k i /∂t and k (II)   i summarize partial time derivatives and all lower order terms in the first and second time derivative of (7), respectively.They may depend on the (absolute) coordinates p 0 of the root body, on the absolute coordinates p := ( p 1 , . . ., p N ) of the remaining N bodies in the system, on the corresponding joint coordinates q := (q 1 , . . ., q N ) and on ṗ0 , ṗ and q.
Note, that ( 9) is simplified substantially for all bodies (•) (i) that follow directly the root body (π i = 0) since the root body is inertially fixed ( pπ i = 0): In recursive multibody formalisms, it is supposed that position and velocity of the root body (p 0 (t), ṗ0 (t)) and all joint coordinates q i (t), qi (t), ( i = 1, . . ., N ), at a current time t are known.Starting from the root body, the absolute position and velocity coordinates p i (t), ṗi (t) of all N bodies (•) (i) , ( i = 1, . . ., N ), may then be computed recursively using ( 7) and ( 8), respectively, (forward recursion).

Tree structured multibody systems: equilibrium conditions
The equations of motion of a multibody system with N bodies may be obtained from the equilibrium conditions for forces and momenta for each individual body that are formulated in absolute coordinates p i : The body mass matrix M i ∈ R d×d contains mass and inertia tensor of body (•) (i) and is supposed to be symmetric, positive semi-definite.The equilibrium conditions contain the reaction forces of the joints connecting body (•) (i) with its predecessor (K i µ i ) and with its successors in the kinematic tree (H j µ j , j ∈ I i ).All remaining forces and momenta acting on body (•) (i) are summarized in the force vector The specific structure of the joint reaction forces with Lagrange multipliers µ i (t) ∈ R d that satisfy results from the joint equations ( 7) and from d'Alembert's principle since the virtual work of constraint forces vanishes for all (virtual) displacements being compatible with (7).In (12), matrix J i denotes the Jacobian of the constraint function k i w.r.t.joint coordinates q i ∈ R n i , see Sect.3.1 above.
For leaf bodies (•) (i) , the equilibrium conditions (11) get a simpler form since I i = { j : π j = i } = ∅.We obtain One of the basic components of recursive multibody formalisms are algorithms to transform the equilibrium conditions (11) recursively for all bodies (•) (i) to the simpler form (13) with suitable Mi and f i .With the common assumption that all body mass matrices M i are non-singular, µ i may be expressed in terms of pπ i , Mi , f i and k (II)  i by block Gaussian elimination applied to (9), ( 12) and ( 13), see, e.g., Lubich et al. (1992).
It is an important observation that this backward recursion may be generalized to multibody systems with rankdeficient body mass matrices M i as long as all body mass matrices M i are symmetric, positive semi-definite.In the following, this will be shown by mathematical induction: let us suppose that the equilibrium conditions of all successors (•) ( j) of body (•) (i) are given in form (13), i.e., Since µ j belongs to the null space of J j , see ( 12), we get with the Moore-Penrose pseudo-inverse (J j M j J j ) + of the projected body mass matrix J j M j J j ∈ R n j ×n j , see Remark 1 in Appendix A. Because of π j = i, the term K j pj is given by see ( 9), and we obtain finally with C := M1/2 j J j , see Lemma 2a in Appendix A. Multiplying (11) from the left by K − i and inserting µ j from ( 16) for all j ∈ I i , the equilibrium conditions for body (•) (i) get the more compact form (13) with Lemma 1 in Appendix A shows that matrix Mi in (17a) is symmetric, positive semi-definite since it is a finite sum of symmetric, positive semi-definite matrices.Starting from the leaf bodies and following all branches of the kinematic tree to the root, the compact form (13) of the equilibrium conditions may be obtained recursively for all N bodies (•) (i) of the multibody system.The condensed mass matrices Mi ∈ R d×d .are symmetric, positive semi-definite.
The backward recursion with Mi , f i being defined by ( 17) is well-defined whenever all body mass matrices M i , ( i = 1, . . ., N ), are symmetric, positive semi-definite and matrices K i are non-singular.Assuming additionally that the body mass matrices M i ∈ R d×d are positive definite, matrices Mi ∈ R d×d and J j M j J j ∈ R n i ×n i in ( 17) are non-singular and (J j M j J j ) + may be substituted by (J j M j J j ) −1 .In that special case, the recursive definitions ( 17) are well known from classical multibody formalisms, see, e.g., Lubich et al. (1992).

Forward recursion: Absolute coordinates
The second time derivative (9) of the kinematic relations (7) defines the acceleration pi of body (•) (i) in terms of the acceleration pπ i of its predecessor (•) (π i ) and in terms of the corresponding joint coordinates qi .Eliminating the joint coordinates, all (absolute) accelerations pi , ( i = 1, . . ., N ), may be computed recursively starting at the root body since p0 ≡ 0 (forward recursion).
Left multiplication of ( 9) by (J i Mi J i ) + J i Mi results in 12) and ( 13).If Mi ∈ R d×d is symmetric, positive definite then J i Mi J i ∈ R n i ×n i is symmetric, positive definite as well and (18) defines an explicit expression for qi because of (J i Mi J i ) + (J i Mi J i ) = I n i in that case.In general, however, Eq. ( 18) determines only r i := rank (J i Mi J i ) ≤ n i components of qi ∈ R n i , see also Remark 1b in Appendix A.
Substituting (18) in ( 9), we get the expression that proves to be useful in the forward recursion if J i Mi J i is rank-deficient.The main difference between the full rank and the rank-deficient case is the additional non-zero term It is an important (and non-trivial) observation that this additional term does not affect the successors of body (•) (i)  in the kinematic tree (if there are any).Suppose I i ∅ and consider a (direct) successor (•) ( j) of body (•) (i) , π j = i.In the condensed equilibrium conditions ( 14), the reaction force between bodies (•) ( j) and (•) (i) is represented by µ j ∈ R d that depends on pi , f j and k (II)  j , see ( 16): From ( 19), we see that the right hand side of (21a) contains the product of matrices that is rewritten as . Lemma 3 from Appendix A proves J j = 0 d×n i resulting in see ( 21a) and ( 19).
The two alternative representations of µ j in (21a,b) provide two different ways to evaluate pj by forward recursion, see ( 14): ) To keep the presentation compact, we assume in the following that the bodies with rank-deficient body mass matrix M i are isolated in the kinematic tree, i.e., the predecessor of a body (•) (i) with rank-deficient M i is either the root body (π i = 0) or a body with non-singular body mass matrix: For π i 0, this assumption implies that Mπ i is non-singular as well since a symmetric, positive semi-definite matrix M π i with rank M π i = d is positive definite and Mπ i is defined by a sum of K − π i M π i K −1 π i and a finite number of symmetric, positive semi-definite matrices, see (17a).
The technical assumption ( 23) allows to evaluate recursively pj for all bodies with non-singular M j starting at the root body (•) (0) that was supposed to be inertially fixed ( p0 = 0).Let a body (•) ( j) be given with rank M j = d and denote its predecessor by i := π j .If π j = 0 or i = π j 0 and the condensed mass matrix Mi of the predecessor is non-singular then pi has been computed before and pj may be obtained from (22a) since K j was supposed to be non-singular: Otherwise, i = π j 0 and Mi is rank-deficient and the corresponding body mass matrix M i has to be rank-deficient as well.Then, body (•) (i) may be skipped in the forward recursion since assumption ( 23) guarantees that pπ i has been computed before and pj may be obtained from (22b): The proposed forward recursion algorithm evaluates pj for all bodies (•) ( j) with non-singular M j provided that the technical assumption ( 23) is satisfied and the body mass matrices M i are symmetric, positive semi-definite for all N bodies of the multibody system model, ( i = 1, . . ., N ).It provides furthermore the algorithmic basis to evaluate the equations of motion as mixed second and first order system for the joint coordinates q = (q 1 , . . ., q N ).

Equations of motion: mixed second and first order system of differential equations
Let us consider again a body (•) ( j) with non-singular M j and denote i = π j .Multiplying (15) from the left by J j M j , we get since J j M j K j pj = J j f j − J j µ j = J j f j , see ( 14) and ( 12).As before, we get an alternative expression from substituting in the right hand side of (25a) the term pi according to (19): The simpler expression (25a) can be used whenever π j = 0 or i = π j 0 and Mi is non-singular since pi has been evaluated by forward recursion in that case.Eq. ( 25b) shows that the situation is substantially more complicated if i = π j 0 and Mi is rank-deficient since in that case body (•) (i) and the corresponding acceleration term pi have been skipped in the forward recursion and pπ i has to be used instead.
In the latter case, the technical assumption (23) guarantees that pπ i is really available from the forward recursion since either π i = 0 or π i 0 and Mπ i is non-singular.Therefore, the right hand side of (25a) with ( j, i) being substituted by (i, π i ) may be evaluated straightforwardly: In ( 26), the coefficient J i Mi J i of qi is symmetric, positive semi-definite and may therefore be diagonalized, see Remark 1 in Appendix A: Here, Λi ∈ R r i ×r i with r i = rank (J i Mi J i ) ≤ n i is a positive diagonal matrix containing the non-zero eigenvalues of 25b) may be expressed as Multiplying ( 26) from the left by X i , we end up with a decoupled system of r i linearly implicit second order differential equations and n i − r i additional equations that do not contain qi and may further be simplified to The equations of motion for the multibody system model are given by (25b) for all bodies (•) ( j) with i := π j 0 and rank Mi < d and by (25a) for the remaining bodies, ( j = 1, . . ., N ).They are composed of i r i linearly implicit second order differential equations (25b), (28a) and i (n i − r i ) additional equations (28b) to define ( 0 (n i −r i )×r i I n i −r i )X i qi , ( i = 1, . . ., N ), see also the detailed discussion in Sect.3.5 below.Similar to a classical residual formalism, see Eichberger (1994), the residuals in (25) may be used to integrate the equations of motion by general purpose DAE solvers like D, see Brenan et al. (1996).

Equations of motion: formal analysis
For a formal analysis of equations of motion (25), we introduce velocity coordinates v 0 := ṗ0 , ( j = 1, . . ., N ), with Because of (29), the joint coordinates q j may also be expressed in terms of v := (v 1 , . . ., v N ): since rank Mi < d and the technical assumption (23) imply qi = X i v i .For all bodies (•) (i) with rank Mi < d, vector v i is split according to with In the full rank case (rank Mi = d), we set η i := v i ∈ R n i and leave ζ i "empty" since r i = rank (J i Mi J i ) = n i , i.e., n i − r i = 0.In the rank-deficient case (rank Mi < d ⇒ r i < n i ), the technical assumption (23) guarantees that M j is nonsingular resulting in η j = v j .With (29) and the diagonalized projector in ( 27), we see that q j may be written as a linear combination of η j and ζ i that is independent of η i : The time derivative of v j in (29) depends on time derivatives of X j = X j (p 0 (t), p(t), q(t), t) , T j = T j ( p 0 (t), p(t), q(t), t) with p = p( p 0 , q, t), see (7).Let Ẋ j = Ẋ j ( p 0 , v 0 , q, q, t) and Ṫ j = Ṫ j ( p 0 , v 0 , q, q, t) be defined such that For bodies (•) ( j) with i := π j 0 and rank Mi < d, we get from ( 29), ( 30) and from the product rule v j = X j q j − T j qi + Ẋ j q j − Ṫ j qi = X j q j − T j qi + Ẋ j X j (v j + T j X i v i ) − Ṫ j X i v i .
Multiplying the equations of motion (25b) from the left by Λ −1 j X j , we observe Λ −1 j X j (J j M j J j ) = Λ −1 j X j X j Λ j X j = X j and end up with Since pπ i has been evaluated by forward recursion and the joint coordinates q are given in terms of q, v, p 0 (t), v 0 (t) and t, see (30), the system of n j first order differential equations (32) may be written as η j = ϕ [2]  j (q, v, t) if i := π j 0 and rank Mi < d .In the same way, (28a) is seen to imply a system of r i ≤ n i first order differential equations ηi = ϕ [2]  i (q, v, t) Finally the n i − r i equations (28b) are written as algebraic equations With these transformations, the equations of motion are re-formulated as differential-algebraic system (2) with i (n i + r i ) differential variables y 0 = (q 1 , . . ., q N , η 1 , . . ., η N ) satisfying (2a) with right hand sides ϕ [1]   i of dimension n i and right hand sides ϕ [2]  i of dimension r i , ( i = 1, . . ., N ), and i (n i − r i ) algebraic variables z 0 = ζ := (ζ 1 , . . ., ζ N ) satisfying (2b) with functions γ i , ( i = 1, . . ., N ).
The algebraic equations (2b) define implicitly the "algebraic" velocity components z 0 = ζ, if the Jacobian ∂γ/∂ζ is non-singular along the solution.In practical applications, this regularity assumption will typically be satisfied if ∂γ i /∂ζ i is non-singular for all bodies (•) (i) with r i = rank Mi < n i , ( i = 1, . . ., N ), which may be achieved by appropriate damping terms in the force vector f i that should depend on the velocity coordinates A more detailed analysis of the regularity of Jacobian ∂γ/∂ζ is subject of further research.

Neglecting inertia forces in multibody systems: two examples
The theoretical analysis of Sect. 3 generalizes the results of Arnold et al. (2010) from chain structured systems to general tree structured systems.In this section, we recall two academic test problems from Arnold et al. (2010) to illustrate the basic steps of these investigations.We consider a chain of two mass points (•) (i) , (•) ( j) in 2-D with i = π j and π i = 0. I.e., body (•) (i) follows in the kinematic chain directly the inertial system ("root") and is the predecessor of body (•) ( j) .In Sect.3, there are no specific physical assumptions on the joints between bodies (•) (0) and (•) (i) and between bodies (•) (i) and (•) ( j) , respectively.Therefore, the resulting set of r i explicit second order differential equations (28a) and n i − r i implicit first order differential equations (28b) for body (•) (i) describes arbitrary joint configurations and has a substantially more complex mathematical structure than the corresponding equations of motion in a classical multibody formalism.
In the test problems, bodies (•) (i) and (•) ( j) are represented by point masses with d = 2 degrees of freedom, see Fig. 3.The root body (•) (0) is inertially fixed resulting in p π i (t) ≡ 0. The absolute coordinates of bodies (•) (i) and (•) ( j)  are denoted by p i = (p i,x , p i,y ) , p j = (p j,x , p j,y ) ∈ R 2 .In this simplified setting, the diagonal mass matrices M i , M j have format 2 × 2 and the joint Jacobians satisfy J i ∈ R 2×n i and J j ∈ R 2×n j .
Bodies (•) (0) and (•) (i) are connected by two linear springdamper elements acting parallel to the x-axis and y-axis, see Fig. 3.The n i = 2 degrees of freedom in this joint are represented by joint coordinates q i (t) = p i (t) ∈ R 2 such that the functions k (I)  i , k (II)   i in (8), ( 9) and (10) vanish identically, K i = I 2 , H i = J i = −I 2 .The free motion of body (•) ( j) in y-direction is represented by the joint coordinate q j,y (t) := p j,y (t) − p i,y (t).
For the configuration in the upper plot of Fig. 3, bodies (•) (i) and (•) ( j) are connected by another linear spring-damper element acting parallel to the x-axis.
The joint has n j = 2 degrees of freedom with joint coordinates q j = (q j,x , q j,y ) ∈ R 2 and q j,x (t) := p j,x (t) − p i,x (t).Functions k (I)  j , k (II) j in ( 8) and ( 9) vanish identically, K j = I 2 , H j = J j = −I 2 .
In absolute coordinates, the equations of motion are given by The equations of motion for coordinates q i (t) are composed of r i = 0 second order differential equations (28a) and the n i − r i = 2 implicit first order differential equations see ( 28b).This result is in perfect agreement with (33a,b) in the limit case m i = 0.The joint coordinates q j = p j − q i are not defined explicitly if Mi = 0.The equations of motion (25b) yield m j ( q j,x + qi,x ) = −d j,x q j,x − c j,x q j,x , (35a) m j ( q j,y + qi,y ) = 0 . (35b) In the lower plot of Fig. 3, the relative motion of body (•) ( j) w.r.t.body (•) (i) is restricted in x-direction by the scalar constraint p j,x (t) = p i,x (t) + l i .The joint has only n j = 1 degree of freedom q j (t) = q j,y (t) with a joint Jacobian J j = −( 0 , 1 ) , K j = I 2 , H j = −I 2 .The update formula (17a) with M i = m i I 2 , M j = m j I 2 results in Mi = M i + M j − M j 0 1 ( 0 1 ) M j 0 1 + ( 0 1 ) M j = m i + m j 0 0 m i .
In (17b), we have f i = f i since k (I) j = k (II) j = 0 and f j = 0.The equations of motion in absolute coordinates are m i pi,x − µ j,x = −d i,x ṗi,x − c i,x p i,x , (36a) m i pi,y = −d i,y ṗi,y − c i,y p i,y , (36b) m j p j,x + µ j,x = 0 , (36c) m j p j,y = 0 . (36d) Because of p j,x (t) = pi,x (t), we obtain µ j,x (t) = −m j p j,x (t) = −m j pi,x (t) and the equations of motion (36a,b) get the form (m i + m j ) pi,x = −d i,x ṗi,x − c i,x p i,x , m i pi,y = −d i,y ṗi,y − c i,y p i,y .
In the limit case m i = 0, a combined set of r i = 1 second order differential equation for q i,x (t) and n i − r i = 1 first order differential equation for q i,y (t) is obtained, see also (28a) and (28b) with J i = −I 2 , X i = I 2 and Λi = m j : m j qi,x = −d i,x qi,x − c i,x q i,x , 0 = −d i,y qi,y − c i,y q i,y .
For this second test problem, we have a scalar joint coordinate q j = q j,y ∈ R that is again not explicitly defined but has to satisfy the linearly implicit second order differential equation m j ( q j,y + qi,y ) = 0 , see ( 25b) and (35b).

Conclusions
Motivated by results from singular perturbation theory, multibody system models with bodies of small mass or nearly singular inertia terms are analysed considering the limit case of systems with rank-deficient body mass matrices.Replacing in a classical recursive multibody formalism the inverse of condensed body mass matrices by their Moore-Penrose pseudo-inverse, the backward recursion phase may be adapted to the rank-deficient case.The crucial point in the analysis is the evaluation of accelerations for successors of bodies with rank-deficient body mass matrix in the forward recursion phase.It was shown that bodies with rank-deficient body mass matrix may simply be skipped in forward recursion.The acceleration coordinates of joints leaving such bodies to one of its successors are not given in explicit form but satisfy a linearly implicit equation that may be handled conveniently by common general purpose DAE solvers.
For each body with rank-deficient body mass matrix, a mixed system of first and second order differential equations is obtained resulting in a first order DAE that has index 1 for multibody system models with appropriate damping terms in the force elements acting at the "zero mass" body.Further investigations will be necessary to analyse practical aspects of this index-1 assumption in more detail.
In future research, the basic framework that has been developed in the present paper for tree structured rigid multibody system models will be extended to flexible systems and to multibody system models with (holonomic) constraints.These additional results will provide the algorithmic basis for a reference implementation in industrial multibody system simulation software.