Model reduction of the tippedisk: a path to the full analysis

The tippedisk is a mechanical-mathematical archetype for friction-induced instability phenomena that exhibits an interesting inversion phenomenon when spun rapidly. The inversion phenomenon of the tippedisk can be modeled by a rigid eccentric disk in permanent contact with a flat support, and the dynamics of the system can therefore be formulated as a set of ordinary differential equations. The qualitative behavior of the nonlinear system can be analyzed, leading to slow–fast dynamics. Since even a freely rotating rigid body with six degrees of freedom already leads to highly nonlinear system equations, a general analysis for the full system equations is not feasible. In a first step the full system equations are linearized around the inverted spinning solution with the aim to obtain a local stability analysis. However, it turns out that the linear dynamics of the full system cannot properly describe the qualitative behavior of the tippedisk. Therefore, we simplify the equations of motion of the tippedisk in such a way that the qualitative dynamics are preserved in order to obtain a reduced model that will serve as the basis for a following nonlinear stability analysis. The reduced equations are presented here in full detail and are compared numerically with the full model. Furthermore, using the reduced equations we give approximate closed form results for the critical spinning speed of the tippedisk.


Introduction
Various gyroscopic systems which are interacting with a horizontal frictional support, such as Euler's disk [19,22,26], the rattleback [5,13] and the tippetop [6,10,25,29], form a scientific playground for frictioninduced instability phenomena. Even though the early research on the tippetop dates back to the 1950s, it is still a topic of current scientific interest [4,18,20]. In [30] we introduced a new mechanical-mathematical archetype, called the tippedisk, to this scientific playground and provided a suitable model. The tippedisk is essentially an eccentric disk, for which the center of gravity (COG) does not coincide with the geometric center of the disk. Neglecting dissipation due to spinning friction (i.e., pivoting friction), two stationary motions can be distinguished. For 'noninverted spinning', the COG is located below the geometric center and the disk is spinning with a constant velocity about the in-plane axis through the COG and the geometric center. The second stationary motion is referred to as 'inverted spinning', being similar to 'noninverted spinning', but with the COG located above the geometric center of the disk. If the noninverted tippedisk is spun fast around an in-plane axis, one observes that the COG rises until the disk remains in an inverted configuration, similar to the inversion of the tippetop shown in Fig. 1. For dissipative systems, people tend to intuitively assume trajectories to end in equilibria or stationary solutions with lower potential energy. Since the potential energy for the noninverted configuration is lower than in the case of inverted solutions, this energetic intuition contradicts our experiments from Fig. 2. The experimental observations qualitatively indicate that for a fastly spinning tippedisk the noninverted spinning is unstable and a stable inverted spinning motion attracts nearly all trajectories.
In this paper we aim to conduct an in-depth stability analysis for the tippedisk. We therefore have to reduce the complexity of the model, to pave the way for a closed form stability analysis. This reduction has to be understood as simplification of the full model and will be validated through various numerical models and the results will be compared to experimental data in future research. In this paper, the main focus is on linear stability analysis and a physically motivated reduction of the system.
In Sect. 2, we introduce the kinematics and briefly recapitulate the model developed in [30]. Moreover, we discuss the dimensions of our specimen and illustrate the two stationary solutions whose stability is studied in the following. Section 3 attempts to explain the local stability behavior of the tippedisk using Lyapunov's indirect method. Physical constraints are numerically and analytically introduced in Sect. 4, to obtain a reduced system which qualitatively describes the inversion phenomenon. The corresponding simulation results are compared in Sect. 5. In Sect. 6, the reduced system is locally approximated to obtain the linear behavior. In addition, a closed form approxima-tion of the critical spinning speed is derived in this section. The behavior of the linear system is subsequently discussed together with the applied approximations in Sect. 7.

Model of the tippedisk
We introduce an orthonormal inertial frame I = (O, e I x , e I y , e I z ) attached to the origin O, where e I z is perpendicular to the flat support. The body fixed B-frame, B = (G, e B x , e B y , e B z ), is attached to the geometric center G of the disk, such that e B z is normal to the surface of the disk. The axis e B x is defined as the normalized vector of r G S , which points from the geometric center G to the center of gravity S. The inertia tensor with respect to G expressed in the body fixed B-frame is given as To describe the interaction between disk and support, the point with minimal height is introduced as C. The contact distance between the contact point C and the corresponding point D is measured as g N (Fig. 3).

Equations of motion
In [30] we provide a suitable mechanical model of the tippedisk, which is able to describe the inversion phenomenon. Using the parameterization of the geometric center G together with Euler angles ϕ = α, β, γ T , in z-x-z convention, we obtain equations of motion in the form of in the coordinates q = [x, y, z, α, β, γ ] T ∈ R 6 . For a detailed explanation of Euler angles, we refer to [27]. As mentioned in [17], Euler angles can lead to singularities. Euler angle singularities are often mistakenly confused with gimbal lock [17]. In [30] we additionally present a quaternion-based approach and showed that the model in Euler angles is suitable to describe the inversion of the tippedisk, since the motion occurs far from singularities of the chosen parameterization. The vector h(q,q) contains all gyroscopic terms. On the right-hand side of Eq. (2) appear the gravitational force f G , the normal contact force λ N ∈ R, with corresponding normal force direction w N , as well as the tangential contact force λ T ∈ R 2 in the contact plane (e I x , e I y ), with associated matrix of tangential force directions W T . Equation (2) represents the equations of motion of a single rigid body in the sense of nonlinear mechanics, which has structurally the same form as a general rigid multibody system. Constraint forces are thereby represented through Lagrange multipliers resulting from constraint equations and set-valued force laws [1,9,15]. In the field of Nonsmooth Mechanics, there are a vari-ety of set-valued frictional force laws such as dry Coulomb friction, spinning friction and contour friction [15,21,24]. Often different friction laws are considered in isolation. In [11] it was shown that spatial Coulomb friction and spinning friction (i.e., pivoting friction) must be considered in a coupled fashion. This leads to more advanced and realistic friction laws such as the Coulomb-Contensou friction presented in [23]. Simplified, Coulomb-Contensou friction law results from a microscopic consideration in which a tangential force distribution in the sense of Coulomb it is assumed to act on a contact patch. The macroscopic contact force and contact torque are then obtained by integrating the force distribution over the contact area. If the tippedisk is spun around an in-plane axis, the contact patch is forced into the state of sliding. This sliding then leads to a smoothing of Coulombs friction law for sliding friction. In [30] we showed that the smooth Coulomb friction law with friction coefficient μ and smoothing parameter ε is sufficient to describe the inversion phenomenon of the disk. The kinematic quantity γ T describes the relative slip velocity 1 between the contact points C and D in the inertial frame. Moreover, the normal contact force λ N is given through Signorini's law where the expression a ⊥ b means ab = 0.
To obtain more compact expressions in the following, we introduce the notations sϑ = sin ϑ, s 2 ϑ = sin 2 ϑ, cϑ = cos ϑ and c 2 ϑ = cos 2 ϑ. The corresponding mass matrix M(q) and the gyroscopic h-vector, as well as the generalized normal and tangential force directions w N and W T , are given in "Appendix A".

Dimensions
Similar to [30], a stainless steel disk is considered, with dimensions depicted in Fig. 4. The nonlinearities induced by the rounding of the edge are not responsible for the inversion of the tippedisk. Neglecting these additional geometric effects, we model the tippedisk as an infinitely thin disk with mass m, radius r and eccentricity e. For the experimental investigations, a disk with a rounded edge has been constructed, as a sharp edge cuts grooves in the support and leaves marks on it, resulting in inhomogeneous frictional properties. In addition, a sharp edge is not as resistant as a rounded one, such that a flattening of the tread can occur. Such flattening of the running surface can lead to micro-jumps of the tippedisk. The dimensions and inertia properties for the considered specimen are given in Table 1. For a more detailed derivation of the principal moments of inertia, we refer to [30]. In the following numerical simulations, the smoothing parameter is set to ε = 0.1 m s and the friction coefficient is chosen as μ = 0.3. In accordance with the model presented in [30], the tippedisk has two stationary solutions for which dissipation is absent. The stationary solutions are characterized by a constant spinning speedα = Ω, a constant height z = r and a constant inclination angle β = π 2 . The COG is below G (γ = − π 2 ) for noninverted spinning and above G (γ = + π 2 ) for inverted spinning. The stationary solutions can be expressed in the coordinates q(t) as with generalized velocities u * n (t) = u * i (t) = [0 0 0 Ω 0 0] T . In the following, we aim to analyze the stability properties of both stationary spinning solutions, depicted in Fig. 5.

Linear stability analysis (6 DOF)
As a first step in the analysis of the dynamic behavior of the tippedisk, the local stability properties of inverted an noninverted spinning solutions are analyzed using Lyapunov's indirect method. The experiments and simulations from [30] show that the inverted spinning solution seems to be locally attractive for a large spinning speed Ω. We may try to apply Lyapunov's indirect method, i.e., use an eigenvalue analysis, to prove local attractivity. Hence, we linearize the system equations around the noninverted and inverted stationary spinning solutions. Both stationary spinning solutions Eq. (5) are characterized by the constant rotational velocityα = Ω and thus by a linear time dependence α(t) = Ωt + α 0 of the angle α. If we linearize the equations of motion around the noninverted and inverted spinning solutions q * (t), we obtain a linear system of the form with y(t) = q(t) − q * (t) and M(t) = M(q * (t)). The time dependent matrices B(t) and C(t) are extracted from h(q * , u * ), the normal contact force λ N (t) = λ N (q * , u * ) and the friction force λ T (t) = λ(q * , u * ).
Since the mass matrix M(t) = M(q * (t)) here depends on α and is thus time-dependent, the linear system is nonautonomous. As a consequence, Lyapunov's indirect method is not applicable due to time-dependent matrices and one would need to resort to Floquet theory which is not feasible in closed form [12]. In contrast, the physical interpretation suggests that the system does not explicitly depend on time t. For this reason, it is presumed that a non-autonomous description exists, so that this time dependency vanishes. In Sect. 3.1, it is shown that a non-autonomous description can be achieved through a coordinate transformation, so that subsequently Lyapunov's indirect method can be applied. The above linearization motivates the introduction of new minimal coordinates for the nonlinear system.

Reparametrization
The parameterization of the reference point G with respect to the I -frame leads to a time-dependent mass matrix M and therefore to a non-autonomous system in coordinates Since this time dependence is an artifact of the chosen parameterization, we introduce a second parameterization in new minimal coordinates where the position of the geometric center G is expressed with respect to the co-rotating R-frame, which results from rotation of the I -frame with the angle α around the e I z -vector.
The corresponding transformation matrix A I R yields the relationship between coordinates q and z. The kinematic relation betweenq andż can be derived by taking the time derivative of Eq. (10) aṡ from which one extracts the kinematic matrix Together with Eq. (10) and Eq. (12), the equation of motion from Eq. (2) transforms to which can be written in short form as with S. Sailer, R. I. Leinē where the notation cγ = cos γ and c 2 γ = cos 2 γ has been used and with Equations (16)-(29) define the left-hand side of the equation of motion, see Eq. (15). The generalized gravitational force yields The generalized normal and tangential force directions are given with respect to the coordinates z as and With γ T =W T Tż the smooth Coulomb friction law Eq. (3) can be evaluated.

Linear stability analysis of the autonomous system
In Sect. 3.1, the reparameterization of the equations of motion to new minimal coordinates z has been introduced. The resulting equations of motion Eq. (15) can be linearized together with the smooth Coulomb friction law (3) and the condition g N = z − r sβ = 0 around the inverted and noninverted stationary spinning solutions of the tippedisk. Due to symmetry, the linearized equations for noninverted spinning solutions are similar to the linearized equations of motions for inverted spinning and differ only in the sign of the eccentricity e. For this reason the linearization is here only performed for the inverted spinning tippedisk. In the inverted orientation β = +π/2 and γ = +π/2 are valid. Introducing the shifted anglesβ = β − π/2 1 andγ = γ − π/2 1, we linearize the trigonometric expressions Instead of the linearization around an isolated state, we assume only to be small and define the vector Note that z and α are not small. We defineẑ := xȳβγ T of small quantities and use O(||ẑ||) to qualify the deviation from the stationary solution. In [30], it is shown that during the inversion process the disk makes always contact with the flat support, i.e., the contact between the points C and D is closed so that g N = 0 holds. Accordingly, the vertical height of the geometric center G and therefore the coordinate z can be expressed as z = r sin β = r cosβ. With the assumption of slightly perturbed inverted motions, the second time derivativë of z vanishes in the linear analysis. From the third row of the equation of motion Eq. (15), we then conclude for motions in the vicinity of the inverted and noninverted stationary solutions of the tippedisk. If the disk moves in a stationary solution, the relative velocity γ T between contact points C and D is equal to zero. For motion in the vicinity of the stationary solutions, the smooth Coulomb friction law from Eq. (3) can be approximated as linear dissipative force law with dissipation constant d = μλ N ε . Inserting the relative velocity γ T =W T Tż yields the generalized forcē Neglecting higher-order terms in the fourth row of Eq. (15) we obtainα = 0 + O(||ẑ||) 2 , such thatα = Ω and α(t) = α 0 + Ωt follows. The associated solution which can be equivalently expressed in z-coordinates is similar to the inverted stationary solution of the tippedisk. Without loss of generality, the initial angle α 0 can be set to zero since we assume an isotropic friction law. The stationary solution from Eqs. (43)-(46) is derived from the assumptions (37) neglecting higherorder terms. Consequently, the normal contact force λ N used in the tangential Coulomb friction law and the rotational velocityα = Ω in the linearization are no longer degrees of freedom but become parameters of the reduced linear system. Since for each Ω ∈ R, there exists a pair of inverted/noninverted stationary solutions, the parameter Ω defines a foliation in the state-space. Together with the matrix  15), we have to pre-multiply with the C matrix, which resembles the deletion of the third and fourth line. Applying the linearization around zero and usinġ α = Ω yields the linear system with constant system matrices and In Eq. (51) the abbreviation D = (A − B − C) is used. The linearization of the gyroscopic h-vector leads to a skew-symmetric gyroscopic matrix G h 4×4 and a symmetric stiffness matrix K h 4×4 . The gravitational and normal contact forces induce a symmetric K f 4×4 matrix. The linearized Coulomb friction includes linear terms inż, which are distributed by the symmetric matrix D C 4×4 on generalized coordinatesẑ. Moreover, Coulomb friction leads to linear terms inẑ, which occur symmetrically with B 4×4 sym and asymmetrically with B 4×4 skew . Defining sym and N 4×4 := B 4×4 skew , the linearized equations of motion can be written as and form a system of linear ordinary differential equations with constant coefficients. The matrices M 4×4 , D 4×4 and K 4×4 have a symmetric structure. The matrices G 4×4 and N 4×4 are skew symmetric. Since these matrices and hence the eigenvalues, depend on the spinning velocity Ω, the linear stability behavior of the system can be studied with respect to the bifurcation parameter Ω.
Due to the high dimension of the corresponding firstorder system, the eigenvalues are calculated numerically. Figure 6 shows the real and imaginary part of the eigenvalues λ 1 − λ 6 for Ω ∈ [0, 50 rad/s]. Both real eigenvalues λ 7 ≈ −9·10 2 1/s and λ 8 ≈ −1.5·10 3 1/s are strongly negative. As the real part does remain almost constant and does not change sign and therefore is not mainly responsible for the qualitative stability behavior, the corresponding two-dimensional subspace remains stable. The comparison of the upper left and right graph shows that the eigenvalues λ 1 (red) and λ 2 (red, dashed) are crossing the imaginary axis at the critical spinning velocity Ω crit 2 ≈ 30.2 rad/s, such that the magnitude of their real part is negative for Ω > Ω crit 2 . Since the eigenvalues are crossing the imaginary axis as complex conjugated pair, the bifurcation can be identified as Hopf bifurcation with corresponding stable two-dimensional eigenspace for super critical spinning velocities. The green pair of conjugated eigenvalues λ 3 and λ 4 are purely imaginary and their magnitude is directly related to the spinning velocity Ω. The eigenvalues λ 5 (blue) and λ 6 (blue, dashed) are crossing the imaginary axis simultaneously at Ω crit 1 ≈ 27.1 rad/s. The bifurcation at Ω crit 1 can also be identified as a Hopf bifurcation leading to an unstable subspace for supercritical spinning velocities Ω > Ω crit 1 . In the lower graph of Fig. 6, both critical velocities Ω crit 1 and Ω crit 2 are depicted in a magnified plot. The numerical eigenvalue analysis indicates that the equilibriumẑ 0 =ż 0 = 0 corresponding to a inverted spinning solution is always unstable, since there exists for each Ω an unstable subspace. This local stability analysis does not reflect the physical observation of an inverted spinning motion which, loosely speaking, seems to attract almost all trajectories. According to this contradiction, it is not possible to analyze the asymptotic dynamical behavior applying Lyapunov's indirect method to the six-dimensional system Eq. (15). Just because of the vanishing real parts of λ 3 and λ 4 , a stability statement about the nonlinear system by neglecting higher-order terms is not possible.
For the sake of completeness, the eigenvalues for the noninverted spinning tippedisk are shown in Fig. 11 of "Appendix B". Here we observe that noninverted stationary solutions are always unstable because of the positive real part of λ 1 . This stability behavior is consistent with experimental observations, since it is not possible to rotate the disk so that it remains in a noninverted configuration. In summary, the linearized stability analysis contradicts the physical observation for inverted spinning solutions. The reason for this contradiction could be the restrictive assumption of small |x| 1 and |ȳ| 1. This hypothesis is supported by the fact that in the real experiment we observe horizontal movements of the tippedisk during the inversion phenomenon. Although the linear analysis cannot properly describe the qualitative dynamics, it provides important results. On the one hand we observe that the inverted spinning of the disk, which is basically a stationary solution, manifests as equilibrium in generalized coordinatedẑ ∈ R 4 . On the other hand the linearization of the six-dimensional system motivates a constant spinning velocityα ≈ Ω and a constant normal contact force λ N ≈ mg. Furthermore, the large difference in the eigenvalues implies that the qualitative behavior can be decomposed into slow and fast dynam- ics, suggesting the application of singular perturbation theory and thus the theory of slow-fast systems.

Reduction
In the previous section, it has been shown that the qualitative dynamic behavior of the tippedisk cannot be identified by an eigenvalue analysis of the full system. For this reason, a system reduction is sought which preserves the qualitative stability. Before introducing new constraints, the bilateral constraint of the contact point C is applied, such that an ordinary differential equation is obtained. This bilateral constraint was motivated in [30] and numerically validated. After this first reduction step, new physically motivated constraints are introduced which are validated numerically using simulations with fixed initial conditions. For the sake of clarity, Table 2 introduces the following model names with the assumptions made.
As several constraints will be introduced in this section, the reduction procedure is explained in a general way. The starting point of any reduction step is an 'unconstrained' mechanical system of the forṁ − h(q, u) = f(q, u).
The addition of a generic mechanical constraint equation c(q, u, λ C ) = 0 ∈ R m induces the generalized constraint force f C = W C λ C into the right-hand side of the equation of motion (56), where W C is the matrix of generalized force directions and λ C ∈ R m is the associated constraint force (i.e., Lagrange multiplier). This addition yields the constrained systeṁ which forms a differential algebraic equation (DAE) [16]. In DAE-theory it is possible to obtain the underlying ordinary differential equation (ODE) by successive differentiation of the constraint equation. The (differentiation) index of a DAE counts the differentiations needed to obtain an ODE for all states and Lagrange multipliers [14,16]. In mechanical systems, we restrict us to constraints on position, velocity or acceleration level. Constraints on position level g C (q) = 0 can be differentiated once to obtain the corresponding constraint on velocity level Bilateral constraint g N = 0 Hor. fixed COG g S = 0 No tangential slip γ x = 0 Constant spinning speed γ α = 0 The meaning of the constraints will be explained in following sections Subsequently, the differentiation of a constraint on velocity level yields the constraint on acceleration level Since the mass matrix M is symmetric, positive definite and thus invertible, we obtain the generalized acceleratioṅ For perfect constraints it holds that the matrix W C of generalized force directions follows from the kinematics: Substitution of Eq. (60) allows to express the relative acceleration as a function of q, u and λ Ċ The square matrix D = W T C M −1 W C , defined in Eq. (62), is called the Delassus matrix [8]. For linearly independent constraintsγ , the matrix W C of generalized force directions is of full column rank. The Delassus matrix D is then symmetric, positive definite and thus invertible, such that an explicit equation for λ C can be found. Direct substitution of the explicit equation λ C (q, u) into the equations of motion yields the underlying ODE. However, the DAE (57) with constraint Eq. (62) is a DAE of index 1, as the constraint equation needs to be differentiated once to obtain a differential equation in (q, u, λ C ). In summary, a constraint on acceleration level describes an index 1 problem, while constraints on velocity or position level lead to index 2 and index 3 DAEs, respectively. Instead of solving the related ODE, we use a numerical scheme, which solves the DAE of index 1 directly. Considering the constraint on acceleration level Eq. (62) in system (57), we obtain the DAE of index 1 which can be rewritten in matrix form as ⎡ The square matrix A is of full rank, as it can be transformed into the upper triangular matrix with a submatrix A * 33 = W T C M −1 W C , which is equal to the Delassus matrix D an thus of full rank. Equation (64) can be solved for each time step and therefore integrated with any standard integrator for ordinary differential equations (e.g., the stiff integrator ode15s from MATLAB). To avoid numerical difficulties, the initial conditions must be admissible, i.e., they must be compatible with the applied constraints. The used initial conditions are therefore listed in Table 3 and are chosen such that the following constraints are not violated.

Bilateral constraint
During the inversion process, one observes that the unilateral constraint g N ≥ 0 expressing the impenetrability of the contact remains closed, i.e., g N = 0 see [30].  The unilateral constraint can therefore be replaced with a bilateral one which forces the contact point to be in touch with the flat support and prevents penetration. Of course, this bilateral restriction does not apply to general motion of the tippedisk, but it is indeed fulfilled during the inversion process and therefore not an approximation. The equations of motion from (15) together with the constraint equation (66) can be written in a differential algebraic form The tangential generalized contact force λ T , described by the friction law (3), depends linearly on the normal contact force λ N . The index of the DAE (67) can be reduced by differentiating the constraint equation, i.e., replacing it with the velocity constraint, yielding the DAE of index 2 Subsequently, we reduce the index of system (73) by formulating the constraint on acceleration level to arrive at the DAE of index 1 As the first equation of (75) describes the equation of motion, which is basically a differential equation of second order, we introduce the kinematic equation such that the system (75) can be rewritten in first-order form aṡ Inserting the force law (3) of smooth Coulomb friction into the right-hand side of the kinetic equation, we obtain as right-hand side (RHS) depending on the kinematic quantitiesq,u and the normal force λ N . The λ N -dependence in the Coulomb friction law renders the matrix A 0 asymmetric and is therefore not guaranteed to be invertible. Using the results from the previous section, where the system was linearized with six degrees of freedom, the assumption λ N ≈ mg for calculation of the friction force is motivated. Using this approximation, Eq. (79) is rewritten as ⎡ with i.e., λ T = λ T (q, u). The assumption of a constant normal contact force mg in the Coulomb friction law has led to a symmetric matrix A 1 , having the same structure as the matrix A in Eq. (64), guaranteeing the solvability of Eq. (80). In Fig. 7, the angles β and γ are shown during the process of inversion. The γ -graph (red, dashed) starts from γ 0 = −π/2+ 0.1 rad and converges in the interval t ∈ [0 s, 0.5 s] to γ = π/2. As the angle β (blue, dotted) only changes slightly, we identify the increase in γ with the inversion of the tippedisk. Both angles are superimposed with small oscillations occurring from t ≈ 0.5 s. The colored angles were obtained by solving the DAE of index 1 from Eq. (79). The angles shown in black are calculated assuming a constant normal force λ N ≈ mg, i.e., by combining Eqs. (80) and (81). Since the normal force λ N is nearly constant, it is convenient to neglect the λ N -dependence in the tangential friction law and to assume a constant scaling. This has already been motivated by the linear stability analysis. The difference between the solutions with constant and varying normal force is nearly zero. To be more specific, the absolute difference | β| and | γ | lies in the range of 10 −2 rad. For the sake of simplicity, we will use only the approach of a constant normal force in the following.

COG constraint
The horizontal position of the tippedisk does not affect the dynamical behavior, since the equations of motion (w.r.t the inertial I -frame) do not explicitly depend on the horizontal position, see "Appendix A". The only force acting in the horizontal direction is caused by friction and is therefore dependent on the relative sliding velocity at the contact point. The relative sliding velocity can be decomposed into a part caused by horizontal translation and a part induced by gyroscopic effects. As the gyroscopic effects dominate, the part induced by horizontal motion is negligible, motivating the assumption of a horizontally fixed center of gravity. In experiments, we observe that the disk only drifts slowly in the horizontal plane, supporting the hypothesis of a horizontally immobile COG. The COG constraint can be expressed on position level with

No tangential slip
The relative velocity between the moving contact point C and the inertial fixed support is equal to the horizontal projection of the velocity v C of the contact point C and can be decomposed into two orthogonal parts. For a more physical interpretation, we use the co-rotating R-frame with and define the projected relative velocity components γ x := e R x · v C and γ y := e R y · v C as tangential and lateral slipping velocity. Motivated through experiments we assume that the tippedisk inverts without any slip in e R x -direction, which justifies the velocity constraint At his point we mention that the assumption of pure rolling (v C = 0) does not lead to the inversion of the tippedisk and therefore does not describe the phenomenon under investigation. Figure 8 shows the simulation results of Model 3, using initial conditions of Table 3 under the assumption of a bilateral constraint, a horizontally fixed center of gravity and zero slip condition in e R x -direction, introduced in (87). Qualitatively, the simulation results of Model 3 are similar to Model 2.

Constant spinning velocity
The linearization of the system (15) reveals that the spinning velocityα is constant for motion in the vicinity of an inverted spinning solution, such thatα = Ω holds. Furthermore, simulation results of the full model show that the angle α seems to increase almost linearly during the whole inversion process of the tippedisk.
The assumption of a constant spinning velocityα = Ω can be written with the generalized force direction w α = 0 0 0 1 0 0 T (88) and χ α = Ω = const as a constraint on velocity level Since Eq. (89) describes a holonomic constraint, we integrate the constraint on velocity level to obtain the rheonomic constraint on position level [7]. The equations of motion do not depend on α explicitly. For this reason α 0 can be set to zero without loss of generality.

Reduced system
The application of a bilateral constraint, a horizontal fixed center of gravity, zero tangential slip in e R xdirection and a constant spinning velocity yields the constrained systeṁ of index 1, which can be written in linear matrix form ⎡ In (92) the constraints are applied on acceleration level so that the system can be solved numerically using Fig. 8 Comparison of simulation results of the various models explicit integration schemes. However, this formulation on acceleration level can lead to drift problems and therefore to excessive constraint violation (e.g., penetration). To prevent the constraints from drifting, one may choose a new set of minimal coordinates such that q = q(z) fulfills all constraints on position level intrinsically. In addition, we choose new minimal velocities v :=β and introduce the kinematic differential equation with such that the velocity constraint is automatically fulfilled for all z and for all v. As the velocity constraint from Eq. (98) is non-integrable, it is identified as a nonholonomic constraint [2,3]. The constraints from (91) are said to be perfect in the sense of d'Alembert, as we obtain the reduced systeṁ + β(z, t)) This reduced system in minimal coordinates z and minimal velocities v from Eq. (101) corresponds to a firstorder ordinary differential equation of the forṁ with mass matrix (here scalar) M = A cos 2 γ +B sin 2 γ + m(r + e sin γ ) 2 cos 2 β, generalized gravitational forcē and generalized friction forcew y λ T y with corresponding force direction lateral sliding velocitȳ and friction force In Eqs. (103) and (104) the abbreviationB = B − me 2 is used. The kinematic part of (102) yieldṡ with Equations (92) and (102) describe the bilaterally constrained tippedisk with horizontally fixed center of gravity, zero slip condition in e R x -direction and a constant spinning velocity. Equation (92) describes a DAE of index 1, which can be solved numerically since all constraints are formulated on acceleration level. In contrast, Eq. (102) corresponds to a formulation in minimal coordinates q and minimal velocities v satisfying the introduced constraints, so that the behavior of the system is described by a first-order ordinary differential equation. Since both approaches consider the same constraints, they lead (in the absence of numerical errors such as drift) to the same simulation results, which are shown in Fig. 8.

Comparison
In the previous sections, we have presented several models by applying additional constraints on the system (15). The introduced constraints are used to reduce the number of degrees of freedom, step by step, and thus to obtain a minimal model which qualitatively describes the inversion behavior of the tippedisk. In this section, the solutions of all reduction stages are compared, starting from the initial conditions given in Table 3. The evolution of the angle γ and the inclination angle β is shown in Fig. 8 for Models 1-4.
Qualitatively, the trajectories of the inclination angle β remain close to π/2 for all introduced models. The qualitative comparison of the derived models shows that the behavior is similar. However, it can be noted that the additional constraints lead to less overshooting and faster convergence to the inverted spinning solution. In Fig. 9, the evolution of the spinning velocityα is depicted for the presented models.
The red graph describes the spinning velocityα under the assumption of a bilateral constraint of the contact point, the green under assumption of a bilateral constraint and a horizontally fixed center of gravity. The blueα-graph is obtained from the model with an additional constraint on the tangential sliding velocity. The black graph is given using the model from Eq. (102), which assumes a constant spinning velocity and therefore implies directly a linear time dependency of α. According to Fig. 9, the spinning velocityα(t) drops slightly during the inversion phenomenon froṁ α 0 = 40 rad/s toα ≈ 35 rad/s. This drop inα stems from the fact that the kinetic energy has to decrease when the potential energy rises during the inversion process. The exact decrease here depends on the respective model, but is in a similar order of magnitude. After the disk inverts its orientation, the rotational velocityα is kept nearly constant.
6 Linear stability analysis of the 3 DOF system Similar to Sect. 3, the shifted anglesβ = β − π/2 andγ = γ − π/2 are introduced. The linearization of Eq. (102) around the equilibrium 'inverted spinning' yields the linear homogeneous systeṁ The eigenvalues λ n i from Eq. (111) are computed numerically and shown in color in Fig. 10 for varying Ω. The left graph shows the real part R(λ i ) of the eigenvalues λ i for Ω ∈ [0, 50 rad/s]. The right graph of Fig. 10 represents the corresponding imaginary part I(λ i ). Since the real part of the third eigenvalue λ n 3 yields R(λ n 3 ) < −100 1/s ∀ Ω ∈ [0, 50 rad/s], it is not shown in the left diagram. For Ω = 0 the eigenvalues λ n 1 and λ n 2 are purely real and R(λ n 1 ) > 0, R(λ n 2 ) = 0 holds. If the rotation speed Ω is increased, these two eigenvalues λ 1 and λ 2 with positive real part converge so that from Ω 1 ≈ 12.5 rad/s on they form a complex conjugate pair with equal positive real part. If the spinning velocity Ω is further increased, the complex conjugate pair of eigenvalues λ n 1 and λ n 2 is crossing the imaginary axis, changing the sign of the real part at the numerically determined critical spinning velocity Ω n crit ≈ 31 rad/s. For Ω > Ω n crit all eigenvalues have negative real part, such that the inverted spinning equilibrium is locally asymptotically attractive. As a complex conjugate pair of eigenvalues (i.e., the pair of λ n 1 and λ n 2 ) changes the stability behavior of the equilibrium, we identify the corresponding bifurcation as a Hopf bifurcation.
Due to the dimension of the system (111), a linear stability analysis in closed form is not feasible. However, we can use the smoothing parameter ε as perturbation parameter with lim ε ↓ 0 to approximate the eigenvalues. The characteristic polynomial p(λ) can be calculated using the Laplace expansion as Since the linear first-order system (111) has dimension three, there is at least one pure real eigenvalue such that the characteristic polynomial p(λ) can be written as polynomial decomposition The comparison of coefficients from Eqs. (113) and (114) yields with Eq. (112) from which we derive the mathematical orders From Eqs. (118), (117) and (112) we obtain the coefficients According to the polynomial decomposition (114), the purely real eigenvalue is identified as Solving the quadratic polynomial factor λ 2 +bλ+c = 0 results in In Eqs. (122) and (123), the eigenvalues λ 1,2,3 are approximately given in a closed form under the assumption of a small smoothing parameter ε. These has to vanish. For the limit of ε ↓ 0 ⇒ b → 0, the following expression has to hold. The critical spinning velocity Ω crit is thus given by The numerically determined critical spinning velocity Ω n crit ≈ 31 rad s is close to Ω crit from Eq. (126), such that Ω crit ≈ Ω n crit is a valid approximation. The corresponding complex eigenvalues are given as which also agrees with the numerically calculated eigenvalues.
With prior knowledge of a Hopf bifurcation, λ = iω can be substituted directly into the characteristic polynomial (114), which yields Setting both real and imaginary part to zero, the two equations are obtained. From Eq. (129) it follows which can be inserted into Eq. (130) and solved for Ω = Ω crit , yielding Hence, the critical spinning velocity Ω crit given in Eq. (126) is exact and not an approximation for small values of the smoothing parameter ε. Moreover, it can be seen from Eq. (132) that the critical spinning velocity does not depend on the friction coefficient μ and the smoothing parameter ε.

Discussion
Experiments show that at high spinning velocities the tippedisk inverts so that the center of gravity lies above the geometric center of the disk, see Fig. 2. For lower spinning velocities the disk does not invert. These observations indicate that the inverted spinning solution is unstable for low spinning speeds and becomes attractive when the disk is spun fast enough. In Sect. 3.2, the local stability of the inverted spinning solution of the full system was investigated using a linear eigenvalue analysis. Numerical computations for Ω ∈ [0, 50] rad s showed at low spinning velocities two eigenvalues with a pronounced positive real part. For larger values of Ω, the real part of these eigenvalues becomes negative. Furthermore, the real part of a second complex conjugate pair of eigenvalues increases from a noticeable negative value at lower speeds to a very small positive value at higher rotational speeds. Consequently, the inverted spinning solution of the full system is always unstable, since for all Ω ∈ [0, 50] rad s there exists at least one eigenvalue with positive real part, indicating an unstable subspace.
Numerical simulations of the full nonlinear system show the inversion during a first phase of motion. However, after the solutions reach the inverted spinning solution, they are eventually repelled. The eigenvalues of the full system indeed indicate that there are three different timescales involved. For supercritical spinning velocities, two strongly negative eigenvalues λ 7,8 with R(λ 7,8 )| 1 corresponding to fast asymptotic dynamics are obtained. Moreover, two complex conjugate eigenvalues λ 1,2 with negative real part describe a stable subspace with intermediate spinning velocity (i.e., the asymptotic behavior is neither very slow nor extremely fast). The third and thus slow timescale characterizes the long-term behavior, as the slightly positive real parts of λ 5,6 define an unstable subspace, such that the inverted tippedisk is not stable on longer timescales (because of the eigenvalues with slightly positive real part and also because the spinning speedα decreases as the system dissipates energy), see [30]. Since we are only interested in the inversion process during the first stage of motion, the focus is on the fast and intermediate timescales on which the inversion process occurs. Therefore, we need to somehow 'block' the slow unstable motion related to the eigenvalues with small positive real part. This is exactly what we do by adding physically motivated constraints, as done in Sect. 4. Clearly, these assumptions affect somewhat the quantitative behavior, as the comparison shows that the more constrained models converge faster to the inverted spinning solution. However, the qualitative behavior is preserved, i.e., the inversion process occurs only for supercritical initial spinning velocities and the evolution of the inclination angle β and the angle γ is not strongly affected by the four constraints which we introduced. Two constraints, namely the horizontally fixed center of gravity and the constant spinning speedα = Ω, require further discussion. As the inversion behavior does not strongly depend on the horizontal motion, the constraint of a horizontally fixed center of gravity seems justified, although actually not present for the real system. By introducing the horizontal constraint, the state space is reduced such that the inversion behavior is preserved and superfluous horizontal dynamics is neglected. According to the graphs in Fig. 9, the spinning velocityα remains nearly constant (on an intermediate timescale) after the disk spins in an inverted orientation. Since a raise of the center of gravity implies an increase of the potential energy, the kinetic energy and thus the rotational speed must naturally decrease as a consequence of energetic consistency. The assumption of a constant spinning velocityα correctly describes the qualitative dynamics, but leads to inconsistencies in energetic considerations. In [30] it is shown that the inversion phenomenon cannot be explained from an energetic viewpoint. Therefore, we argue that the violation of the energetic consistency is acceptable.

Conclusions and merits
The tippedisk shows an interesting inversion phenomenon that cannot be explained intuitively and must be studied in the context of nonlinear dynamics. From the stability analysis conducted in this paper, which serves as stepping stone to a more complete global analysis, the following conclusions are drawn: • The asymptotic behavior of the tippedisk strongly depends on the considered timescale. On a 'fast' and 'intermediate' timescale, the inversion can be observed. On a 'slow' timescale, the motion reaches a static equilibrium state, for which the disk lies flat on the table. • Comparison of the models at different reduction levels shows that the qualitative dynamical behavior is well described by the three-dimensional reduced system (102). • Based on the reduced model, the linearization around the inverted spinning motion yields the local stability behavior. For Ω = Ω crit we identify an eigenvalue λ 3 with negative real part and a complex conjugate pair λ 1,2 crossing the imaginary axis, implying a Hopf bifurcation. The classification as sub-or supercritical Hopf bifurcation will be studied in a follow-up paper. • A closed form expression of the critical spinning velocity Ω crit = (r + e) 2 r mḡ B = 30.92 rad s was derived. The critical spinning velocity Ω crit does not depend on the friction coefficient μ and smoothing parameter ε. However, these parameters affect the real part of the eigenvalue and thus characterize the asymptotic behavior of solutions, i.e., how fast solutions converge to the inverted spinning solution.
The reduced system will set the starting point for a further nonlinear dynamic analysis in follow-up papers.