Dynamics of a Rolling and Sliding Disk in a Plane.
Asymptotic Solutions, Stability and Numerical Simulations
Maria Przybylska1* and Stefan Rauch-Wojciechowski2**
1Institute of Physics, University of Zielona G´ora,
Licealna 9, PL-65–417 Zielona G´ora, Poland
2
Department of Mathematics, Link¨oping University, 581 83 Link¨oping, Sweden
Received October 27, 2015; accepted February 15, 2016
Abstract—We present a qualitative analysis of the dynamics of a rolling and sliding disk in a horizontal plane. It is based on using three classes of asymptotic solutions: straight-line rolling, spinning about a vertical diameter and tumbling solutions. Their linear stability analysis is given and it is complemented with computer simulations of solutions starting in the vicinity of the asymptotic solutions. The results on asymptotic solutions and their linear stability apply also to an annulus and to a hoop.
MSC2010 numbers: 37J60, 37J25, 70G45 DOI: 10.1134/S1560354716020052
Keywords: rigid body, nonholonomic mechanics, rolling disk, sliding disk
1. INTRODUCTION
The problem of a disk rolling (without sliding) in a plane is integrable [1, 4, 9, 15]. It is described in the Euler angles by four dynamical equations and this system admits three additional integrals of motion. It can be in principle reduced to a single quadrature. The difficulty is, however, that only the energy integral is given explicitly, while two other integrals, of angular momentum type, are defined implicitly as constants of integrations for the Legendre equation for the dependence of ω3 on
cos θ [22] and [9,§8]. References to original articles treating a rolling disk up to 1992 are contained in the bibliography of the book [18] and more recent results can be found, e.g., in [2, 19, 22].
Equations for a rolling and sliding disk have two more dynamical variables to describe the motion of the center of mass. These equations are dissipative, have monotonically decreasing energy and do not admit any known integrals of motion. A variety of friction models are discussed in the literature. The laws of friction that are used most often in mechanics are the dry Coulomb friction law when the friction force is proportional to the reaction force gn(t) and the law of viscous friction when the
friction force is proportional to the sliding velocity of the rigid body. Coulomb friction considers two cases: rolling/sticking and sliding. The first one occurs when the sliding velocity is zero in the tangent plane of the contact frame and sliding friction occurs when the velocity is nonzero. The value of the sliding friction force is equal to Ff = μgn and for a rolling or sticking contact friction
the Coulomb law is given by the inequality Ff μgn, where μ is the constant friction coefficient.
The dry Coulomb friction law has been used to describe the dynamics of a disk on a plane, e.g., in [16, 17]. The viscous friction has been used for describing the motion of a disk in a plane in [20] and in Chapter 4, §4 of the book [18], see also [12] where a body of revolution was considered. We notice that [18, Ch. 4 §4] contains more general results concerning asymptotic motions of a solid body with an arbitrary shape moving on the horizontal plane and include results of Moshchuk. The value μ = 0 corresponds to the case of an absolutely smooth surface and an appropriate limit
μ→ ∞ corresponds to the case of an absolutely rough surface, see, e.g., [21].
*E-mail: M.Przybylska@if.uz.zgora.pl **
In this paper we study the rolling and sliding motion of a disk under the assumption of viscous friction force Ff =−μgnvA acting against the sliding velocity vAand proportional to the value of
the reaction force gn(t). This law of friction interpolates between the dry friction law Ff =−μgn|vvAA|
and the viscous friction law Ff =−μvA and it is convenient for numerical simulation as the value
of gn(t) is a dynamically determined variable, which, as it is easy to check, stays positive during
numerical simulation. The viscous friction of this form has been used in several papers on the tippe top [8, 24]. The qualitative behavior of the tippe top obtained numerically using the viscous friction law [24] and that obtained from calculations with the Coulomb friction law [5] are very similar. The possibility of detachment of the body from the plane and aspects of the friction law were discussed in [3, 11]. There are also models of friction when the contact between the rigid body and the plane is no longer a point but a small region and the friction force is obtained by integration over the whole contact region, see, e.g., [6, 14, 26].
The purpose of this work is to study behavior of nonbouncing solutions of a sliding disk for the concrete model of friction force Ff =−μgnvA and for initial conditions that correspond to typical
situations of how a disk may be launched into motion. The main tool to understand dynamics is the study of asymptotic (stationary) solutions and their stability and the study of the behavior of solutions in some vicinity of the asymptotic solutions. Stability properties of asymptotic solutions provide a framework for interpreting results of numerical simulations.
For fixing notation we rederive the Euler equations of motion in Section 2. In Section 3 we solve stationary equations with friction and obtain three types of stationary solutions (rolling, spinning and tumbling) with vertical reaction force that are well known for the problem of a rolling disk [1]. They are in agreement with the results of a general study by N. K. Moschuk [20] of asymptotic motion for a rigid body moving in a plane under the action of a viscous friction force as presented in the book by A. P. Markeev [18]. The linear stability analysis of the three types of asymptotic solutions is given in Section 4. It provides only necessary conditions for stability as Jacobian matrices for linearization admit zero eigenvalues. In Section 5 we present results of numerical simulations. The starting point are two types of initial conditions (IC), called here reference IC, which differ from the rolling and spinning solution by a small tilt angle of the symmetry axis. We study them for the range of initial angular velocities, which guarantee positivity of the vertical reaction force
gn> 0 to avoid bouncing from the plane. The results of numerical study of stability for tested
solutions are in agreement with the necessary conditions for linear stability, indicating that these conditions may also turn out to be sufficient.
In the asymptotic behavior the tumbling solutions play a central role as they are attractors for solutions starting far away from the stable rolling and spinning solutions and for solutions with low initial energy. The relation between initial conditions and the final asymptotic solutions is difficult to predict and it shows great sensitivity to small changes of initial conditions for certain ranges of initial velocities.
2. EQUATIONS OF MOTION OF A DISK
In this paper we consider a thin disk of radius R and mass m having a single point of contact
A with a horizontal plane. We assume that the mass distribution in the disk is homogeneous and
that the center of mass CM coincides with its geometric center.
To describe the motion of the disk, we use three reference frames, see Fig. 1. The first one is a fixed inertial Cartesian frame K0 = ( ˆX, ˆY, ˆZ) with basis vectors ˆX, ˆY in the horizontal plane
and ˆZ being perpendicular to this plane. The second reference frame ˜K = (ˆ1, ˆ2, ˆ3), with origin at
CM , has the basis vector ˆ1 is directed to the contact point A of the disk with the plane, the basis
vector ˆ3 is perpendicular to the surface of the disk and is a symmetry axis of the disk. Therefore, ˆ
2 = ˆ3× ˆ1 rests in the plane of the disk and is always parallel to the plane of support. We set a
third coordinate frame K = (ˆx, ˆy, ˆz) with origin at the contact point A so that ˆz = ˆZ, ˆy = ˆ2 and ˆ
x = ˆy× ˆz.
Although the frame ˜K is not a fixed principal reference frame, the symmetry axis ˆ3 is a main
body axis and by symmetry of the disk the inertia tensor is diagonal, I = diag(I1, I1, I3), where I1 := I11= I22= kmR2, I3 = I33= 2kmR2. (2.1)
Fig. 1. Geometry of a sliding and rolling disk.
The value k = 1/4 corresponds to a uniform distribution of the mass of the disk and k = 1/2 is taken for a disk with mass concentrated at the edge.
We denote by s the position of CM with respect to the origin of the frame ( ˆX, ˆY, ˆZ) and by a = Rˆ1 the vector connecting CM to the point A of contact with the horizontal plane.
Two external forces act on the disk: the gravity force −mgˆz applied to the center of mass and the friction-reaction force FR+ Ff applied at the point of contact A. It consists of the reaction
force FR = gnˆz, where we shall assume gn= gn(L, ˙s, s, ˆ3) 0 and the friction force Ff =−μgnvA
acting against the sliding velocity
vA= ˙s + ω× a, (2.2)
with the friction coefficient μ = μ(L, ˙s, s, ˆ3, t) 0, which we later shall consider to be constant.
Thus, in our considerations only the sliding friction is taken into account, while the rolling and the spinning friction is neglected.
The motion of the disk is described by a Newton equation for the linear momentum m˙s, an equation for the angular momentum L =Iω and a kinematic equation describing the motion of the symmetry axis ˆ3 with respect to the fixed inertial frame K0. In vector form the equations are
m¨s = FR+ Ff − mgˆz, L = a˙ × (FR+ Ff) , ˙ˆ3 = ω × ˆ3 =
1
I1
L× ˆ3. (2.3) The energy of the system is the sum of the kinetic energy of translational motion of the mass center Etrans, the rotational energy of the rigid body Erot and of the potential energy of the mass
center Epot
E = Etrans+ Erot+ Epot=
1 2m˙s
2+ 1
2ω· Iω + mgs · ˆz. (2.4)
The change of the energy in time is given by the derivative ˙ E = d dt 1 2m˙s 2+1 2ω· L + mgs · ˆz = ˙s· (FR+ Ff − mgˆz) + ω · [a × (FR+ Ff)] + mg ˙s· ˆz = (vA− ω × a) · (FR+ Ff)− mg˙s · ˆz + ω · [a × (FR+ Ff)] + mg ˙s· ˆz = vA· Ff =−μgnv2A,
where we have used Eqs. (2.3), as well as the equality ˙ω· L = ˙ω · Iω = ω · ˙L, valid due to the
axial symmetry of the disk, and vA· ˆz = 0. In this calculation, vA· Ff =−μgnv2A is the rate of
frictional loss of energy and mg ˙s· ˆz is the rate of energy transfer from the translational part into potential energy Epot. The term ω· [a × (FR+ Ff)] is the work performed in unit time by the
torque a× (FR+ Ff) and it transfers energy between the translational and rotational components.
• pure rolling condition when the sliding velocity vanishes identically in time
vA(t) = ˙s(t) + ω(t)× a(t) = 0, (2.5)
i.e., the contact point is not sliding. This vector equation gives three scalar conditions. Because (2.5) is an identity in time, all derivatives of vA must also vanish. From (2.5) we
obtain ˙s = a× ω and by differentiating this expression we get
m¨s = m(a× ω)·= FR+ Ff − mgˆz, (2.6)
where we have used the first Eq. (2.3). From this equality the friction-reaction force is determined dynamically as FR+ Ff = m(a× ω)·+ mgˆz meaning that at each instant of
time it is exactly such that it sustains pure rolling. Then the second and the third equation in (2.3) become a closed system of equations
˙ L = a× m(a× ω)·+ mgˆz , ˙ˆ3 = ω × ˆ3 = 1 I1 L× ˆ3, (2.7)
which in the Euler angles describes the classical problem [1, 4, 15] of a purely rolling disk having four degrees of freedom.
• the rolling and sliding condition
(s(t) + a(t))· ˆz = 0 (2.8)
holds all time and so all its time derivatives vanish. The first derivative 0 = (˙s + ˙a)· ˆz = (˙s + ω× a) · ˆz = vA· ˆz says that the sliding velocity, indeed, stays in the horizontal plane of
support. From the second derivative multiplied by mass m we obtain
m(˙s + ˙a)·· ˆz = m[¨s + (ω × a)·]· ˆz = (FR+ Ff − mgˆz + m(ω × a)·)· ˆz = gn(ˆz− μvA)− mgˆz + m(ω × a)· · ˆz = 0,
thus the value of the reaction force is equal to
gn= mg + m(a× ω)·· ˆz. (2.9)
From the constraint (2.8) we can calculate the vertical component z of motion of the center of mass, sz = s· ˆz = −a · ˆz and ˙sz=− ˙a · ˆz = − (ω × a) · ˆz. We introduce a vector r that is a projection
of s on the horizontal plane, i.e., r = s− szˆz, and then system (2.3) reduces to
m¨r = Ff, L = a˙ × (FR+ Ff) , ˙ˆ3 = ω × ˆ3 =
1
I1
L× ˆ3. (2.10)
The orientation of the disk with respect to K0 is described by the Euler angles (θ, ϕ, ψ). As
usual, the angle θ is the inclination of the symmetry axis ˆ3 with respect to ˆz, see Fig. 1. We note
that for θ =±π/2 the disk is vertical and for θ = 0, π it is horizontal. Since ˆy = ˆ2, the reference frame K is obtained from K by rotating about ˆy by an angle θ, we have
ˆ
x = cos θˆ1 + sin θˆ3, ˆ1 = cos θˆx− sin θˆz, ˆ
y = ˆ2, ˆ2 = ˆy,
ˆ
z =− sin θˆ1 + cos θˆ3, ˆ3 = sin θˆx + cos θˆz.
(2.11)
We set ϕ as the rotation angle around the ˆz-axis and ψ is the rotation around the the symmetry
The reference frame K rotates with angular velocity ˙ϕˆz with respect to the fixed frame K0
and K is obtained from K by rotation about ˆy by an angle θ. Thus, the angular velocity of the
frame ˜K with respect to K0 is equal to
ωref= ˙θˆ2 + ˙ϕˆz =− ˙ϕ sin θˆ1 + ˙θˆ2 + ˙ϕ cos θˆ3.
The total angular velocity of the disk is obtained by adding to ωrefthe angular velocity of rotation
about the symmetry axis ˆ3 by angle ψ:
ω = ωref+ ˙ψˆ3 =− ˙ϕ sin θˆ1 + ˙θˆ2 + ( ˙ψ + ˙ϕ cos θ)ˆ3.
As usual, we denote ω3= ˙ψ + ˙ϕ cos θ. The angular momentum of the disk is
L =Iω = −I1ϕ sin θˆ˙ 1 + I1θˆ˙2 + I3( ˙ψ + ˙ϕ cos θ)ˆ3.
The kinematic equations describing the rotation of the axes (ˆ1, ˆ2, ˆ3) are the following:
˙ˆ1 = ωref× ˆ1 = ˙ϕ cos θˆ2 − ˙θˆ3,
˙ˆ2 = ωref× ˆ2 = − ˙ϕ cos θˆ1 − ˙ϕ sin θˆ3,
˙ˆ3 = ωref× ˆ3 = ω × ˆ3 = ˙θˆ1 + ˙ϕ sin θˆ2.
The velocity of the point of contact with the horizontal plane (due to (2.11)) is
vA= vxˆx + vyy = vˆ xcos θˆ1 + vyˆ2 + vxsin θˆ3. (2.12)
Substitution of ω, vA and ˆz expressed in terms of the Euler angles into the equations of
motion (2.10) gives us a system of equations that, when solved with respect to (¨θ, ¨ϕ, ˙ω3, ˙vx, ˙vy),
takes the form ¨ θ = sin θ I1 I1ϕ˙2cos θ− I3ω3ϕ˙ +μgnvxR sin θ I1 − Rgncos θ I1 , (2.13) ¨ ϕ = 1 I1sin θ(I3ω3 ˙ θ− 2I1θ ˙˙ϕ cos θ), (2.14) ˙ ω3 =− μRgnvy I3 , (2.15) ˙vx=−( ˙θ2+ ˙ϕ2sin2θ)R cos θ− Rω3ϕ˙ I1 (I1− I3sin 2θ)−μgnvx mI1 (I1+ mR 2sin2θ) + R 2g nsin θ I1 cos θ + ˙ϕvy, (2.16) ˙vy =−R ˙ϕ ˙θ sin θ − (I3+ mR2)μgnvy mI3 − ˙ϕvx . (2.17)
It effectively is a system of sixth order w.r.t. θ, ˙θ, ˙ϕ, ω3, vx, vy. Here the value of the reaction force
is
gn=
mgI1− mRI1sin θ ˙θ2+ mRI1cos2θ ˙ϕ2sin θ− mRI3ω3ϕ sin θ cos θ˙ I1+ mR2cos2θ− μmR2vxsin θ cos θ
. (2.18)
Notice that this system has an invariant manifoldN ={(θ, ˙θ, ˙ϕ, ω3, vx, vy)∈R6| ˙ϕ =vy = ω3= 0}.
The energy of the system in the Euler angles takes the form
E = Etrans+ Erot+ Epot= m 2 v2xcos2θ + (vy− Rω3)2+ (vxsin θ + R ˙θ)2 +1 2 I3ω23+ I1θ˙2+ I1ϕ˙2sin2θ + mgR sin θ. (2.19)
3. ASYMPTOTIC SOLUTIONS
We define asymptotic solutions of Eqs. (2.13)–(2.17) as nondissipative solutions of constant energy ˙E =−μgnvA2 = 0, and so satisfying vA= 0. As a consequence, for asymptotic solutions also
the friction force vanishes Ff =−μgnvA= 0.
Proposition 1. Asymptotic solutions vA= 0 of the rolling and sliding disk are of the following
three types:
• straight rolling solutions with θ = π/2, ˙ϕ = 0 and ω3 an arbitrary constant, • spinning solutions with θ = π/2, ω3= 0 and ˙ϕ an arbitrary constant, • tumbling solutions with θ ∈ (0, π) and
ω3= 0, ϕ =˙ ±
gmR I1sin θ
. (3.1)
They are also stationary solutions of Eqs. (2.13)–(2.17). The reaction force for all these solutions is gn(t) = mg. For spinning and tumbling solutions the center of mass CM is fixed and for the rolling
solutions it is moving along a straight line with speed |˙s| = ω3R. Proof. The asymptotic solutions are determined by the system
˙ ω3 = 0, (3.2) ¨ ϕ = 1 I1sin θ(I3ω3 ˙ θ− 2I1θ ˙˙ϕ cos θ), (3.3) 0 = ˙vy =−R ˙ϕ ˙θ sin θ, (3.4) 0 = ˙vx=−( ˙θ2+ ˙ϕ2sin2θ)R cos θ− Rω3ϕ˙ I1 (I1− I3sin2θ) + R2gnsin θ I1 cos θ, (3.5) ¨ θ = sin θ I1 I1ϕ˙2cos θ− I3ω3ϕ˙ −Rgncos θ I1 , (3.6) where gn=
mgI1− mRI1sin θ ˙θ2− mR cos θ sin θ(I3ω3ϕ˙− I1cos θ ˙ϕ2)
I1+ mR2cos2θ
. (3.7)
From (3.2) ω3= const and then one can integrate Eq. (3.3) to get
(I1ϕ sin˙ 2θ + I3ω3cos θ)·= 0 ⇒ I1ϕ sin˙ 2θ + I3ω3cos θ = C = const . (3.8)
Thus, all asymptotic motions have a first integral C.
Let us note that sin θ = 0 does not give asymptotic solutions because then vanishing of ¨θ
determined in (3.6) requires that gn= 0, while in this case gn= mgI1/(I1+ mR2)= 0. Thus, we
can assume that sin θ = 0, which enables expressing ˙ϕ as a function of cos θ ˙
ϕ = C− I3ω3cos θ
I1sin2θ . (3.9)
Equation (3.4) is now satisfied when either ˙ϕ = 0 or ˙θ = 0. Then Eq. (3.9) implies that
θ = const for every asymptotic solution and ˙ϕ = const. Then the condition of vanishing of the
linear combination R sin θ(¨θ) + ( ˙vx) (see (3.6) and (3.5)) after dividing by −R gives
ω3ϕ = 0.˙ (3.10)
If we substitute gngiven by (3.7) and the condition ˙θ = 0 into (3.6), then we obtain the equality
mgR cos θ + sin θI3ω3ϕ˙− I1cos θ ˙ϕ2
= 0, (3.11)
which substituted into (3.7) gives gn= mg for every asymptotic solution. Since ω3ϕ = 0, condi-˙
tion (3.6) says that for every solution
cos θ[I1ϕ˙2sin θ− gmR] = 0. (3.12)
Thus, every asymptotic solution satisfies ˙
θ = 0, ω3ϕ = 0˙ and cos θ[I1ϕ˙2sin θ− gmR] = 0.
This gives either cos θ = 0 with solution θ =±π/2 or I1ϕ˙2sin θ− gmR = 0 with solution ˙ϕ =
± gmR
I1sin θ and then ω3 = 0. Solutions θ =±π/2 require either ω3 = 0 with arbitrary ˙ϕ or ˙ϕ = 0 with arbitrary ω3. Thus we obtain all solutions mentioned in Proposition 1.
Every asymptotic solution satisfies 0 = vA= ˙s + ω× Rˆ1 and either ˙s = 0 or for rolling solutions
the velocity of the center of mass is
˙s =−ω × Rˆ1 = −ω3ˆ3× Rˆ1 = −ω3Rˆ2
and CM moves with constant velocity directed parallel to the plane of support since here ˙ϕ = 0,
ω3 ∈ R. For spinning and tumbling solutions ˙s = 0, CM is fixed and ω ˆ1. For spinning solutions
cos θ = 0 and ˙ϕ∈ R may be arbitrary. For tumbling solutions cos θ = 0, ˙ϕ = ±
gmR I1sin θ
and
ω3 = cos θ ˙ϕ + ˙ψ = 0.
The asymptotic solutions of Proposition 1 turn the r.h.s. of Eqs. (2.13)–(2.17), together with the equation ˙θ = 0, into zero. So they are also stationary solutions. To show that they are all stationary
solutions, we need to solve the stationary equations of (2.13)–(2.17). From (2.15) ˙ω3 =−μRgnvy
I3 and for gn= 0 there is vy = 0. Then from the Eq. (2.17) 0 = ˙vy =
− ˙ϕvx either vx = 0 or ˙ϕ = 0.
Fig. 2. Bifurcation diagram.
The linear combination
0 = (¨θ)R sin θ + ( ˙vx) =−
μ
mgnvx− Rω3ϕ + ˙˙ ϕvy
Thus, stationary solutions have vx = vy = 0, ˙θ = 0 and they are described by Eqs. (3.2)–(3.7)
with ˙θ = 0. Obviously they have the same solutions.
If it were gn= 0, then by (2.17) there is vxϕ = 0 and either v˙ x = 0 or ˙ϕ = 0. If ˙ϕ = 0,
then all stationary equations are fulfilled but the formula (2.18) gives gn= 0 and we have a
contradiction. If vx= 0, then from (2.13) either sin θ = 0, this describes a disk laying on a plane,
or (I1ϕ˙2cos θ− I3ω3ϕ) = 0. By substituting (I1˙ ϕ˙2cos θ− I3ω3ϕ) = 0 into expression (2.18) we see˙
that again gn= 0 and the assumption gn= 0 leads to a contradiction.
Using the values of integral C of (3.8), we graphically represent asymptotic solutions of the rolling disk in Fig. 2. The bold lines indicate (possibly) stable solutions and the dotted lines unstable ones.
4. LINEAR STABILITY OF ASYMPTOTIC SOLUTIONS
In this section we analyze stability of the asymptotic solutions found in the previous section. In order to reduce the number of parameters, we rescale the variables of system (2.13)–(2.18) in the following way: ( ˙θ, ˙ϕ, ω3) = 1 T(θ , ϕ, Ω 3), (vx, vy) = L T(Vx, Vy), τ = T t, (4.1)
where L and T are the respective scales of length and time. The variables Ω3, Vx, Vy denote
now nondimensional quantities. The derivative denotes differentiation with respect to new nondimensional time τ . We choose normalization in such a way that T = R/g and L = R/g/μ.
In these variables the equations of motion take the form
θ = cos(θ) sin(θ) θ2+ kϕ2 − 1+ sin(θ) Vx− 2ϕkΩ3− θ2Vxsin(θ) k + cos2(θ)− V xsin(θ) cos(θ) , (4.2) ϕ= 2θ sin θ(Ω3− ϕ cos θ), (4.3) Ω3 = Vy sin(θ) θ2− φ2cos2(θ) + ϕΩ3sin(2θ)− 1 2 (cos2(θ) + k− V xsin(θ) cos(θ)) , (4.4) Vx =
4α sin θ cos θ− 4αθ2(k + 1) cos θ− 4αφ2k sin2θ cos θ− 2αφΩ3((2k + 1) cos(2θ)
+ 1) + 2αVx
2θ2(k + 1) sin θ− 2φ2k sin θ cos2θ + (2k + 1)(−1 + φΩ3sin(2θ)) (4.5)
+ cos(2θ)) + 4φVy
cos2θ + k− 4φVxVysin θ cos θ
4cos2θ + k− Vxsin θ cos θ
, Vy = − 8αθϕsin θcos2θ + k+ V x
8αθϕsin2θ cos θ− 8ϕcos2θ + k
+ Vy
4αθ2(2k + 1) sin(θ)− 4αϕ2(2k + 1) sin θ cos2θ + 4α(2k + 1)(ϕΩ3sin(2θ)− 1)
+ 8ϕVx2sin θ cos θ
8cos2θ + k− Vxsin θ cos θ
, (4.6)
where the parameter α is defined by
α = μ gR.
In these equations we have already substituted the expression for gn given by (2.18) and k is either
1/4 for a uniform disk or 1/2 for a hoop. After this rescaling of variables the description of three types of asymptotic solutions takes the form
• straight rolling solutions
• spinning solutions θ = π/2, θ= Ω3 = Vx= Vy = 0, ϕ = ϕ0, • tumbling solutions θ = θ0, Ω3 = Vx = Vy = 0, ϕ =± 1 √ k sin θ0.
Recall that all these asymptotic solutions are equilibria of the vector field v = [θ, θ, Ω3, Vx, Vy]T defined by system (4.2)–(4.6). Equations of the linearized vector field around an equilibrium x0
have the form
Δx = LΔx, (4.7)
with the Jacobian matrix
L = ∂v
∂x(x0), (4.8)
where Δx = [Δθ, Δθ, Δϕ, ΔΩ3, ΔVx, ΔVy]T denote variations of variables x = [θ, θ, ϕ, Ω3, Vx, Vy]T,
respectively, and v in (4.8) is a vector field defined by the r.h.s. of Eqs. (4.2)–(4.6).
The Jacobian matrices of linearization around our asymptotic solutions have the following forms: for straight rolling solutions
Lstraight= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 1 0 0 0 0 1 k 0 −2Ω30 0 1 k 0 0 2Ω30 0 0 0 0 0 0 0 0 0 −2k1 −α k 0 αΩ30 0 − (k+1)α k 0 0 0 0 0 0 −2kα+α2k ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (4.9)
for spinning solutions
Lspin= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 1 0 0 0 0 1 k− ϕ0 2 0 0 −2ϕ0 1k 0 0 0 0 0 0 0 0 0 0 0 0 −2k1 ϕ02−1k α 0 0 ϕ0α −(k+1)αk ϕ0 0 −ϕ0α 0 0 −ϕ0 −2kα+α2k ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (4.10)
for tumbling solutions
Ltumb = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 1 0 0 0 0 1−κ2 −κ3+kκ+κ 0 ± 2√k(κ−κ3) −κ2+k+1 ∓ 2 √ kκ −κ2+k+1 −κ2+k+1κ 0 0 ∓2 1−κ2 kκ3 0 0 0 0 0 0 0 0 0 −2k1 α(κ2−1) −κ2+k+1 0 ± 2ακ3/2√k−kκ2 κ2−k−1 α(κ2+k(2κ2−1)−1) √ kκ(−κ2+k+1) − (k+1)α −κ2+k+1 1 √ kκ 0 ∓α κk 0 0 ∓√1 kκ − 2kα+α 2k ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (4.11)
where we have introduced the notation κ = sin θ0 and double signs in some entries correspond to
two choices of the sign of ϕ for tumbling solutions.
As is well known, a sufficient condition for stability of an equilibrium point is that real parts of all eigenvalues of the Jacobian matrix are negative. For the eigenvalues of matrices (4.9), (4.10) and (4.11) we apply the Routh – Hurwitz conditions to its characteristic polynomials.
Recall that for a polynomial of degree n with real coefficients
f (z) = a0zn+ a1zn−1+· · · + an, (4.12)
with a0> 0, the Routh – Hurwitz conditions are necessary and sufficient for all roots of f (z) to
have negative real parts. The Routh – Hurwitz matrix related to polynomial (4.12) is
H = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ a1 a3 a5 a7 · · · 0 0 a0 a2 a4 a6 · · · 0 0 0 a1 a3 a5 · · · 0 0 0 a0 a2 a4 · · · 0 0 .. . ... ... ... · · · ... ... 0 0 0 0 · · · an−1 0 0 0 0 0 · · · an−1 an ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (4.13)
with entries ak= 0 for k > n. The Routh – Hurwitz conditions for all real parts of the eigenvalues
to be negative is a a sequence of inequalities
Δ1 > 0, Δ2> 0, . . . Δn> 0 (4.14)
for all principal minors Δi of the Routh – Hurwitz matrix H, where i = 1, . . . , n, see [10, Vol II,
Chap. XV, §13].
The characteristic polynomial for the Jacobian matrix related to the linearization along rolling solutions has the form
p1(λ) = λ2 k λ + (2k + 1)α 2k q1(λ), q1(λ) = kλ3+ (k + 1)αλ2+ (4kΩ230− 1)λ + α(2(2k + 1)Ω230− 1). (4.15)
We have here a double zero eigenvalue and a negative eigenvalue −(2k + 1)α/(2k). The principal minors of the Routh – Hurwitz matrix for the third-order polynomial q1(λ) are
Δ1 = α(k + 1), Δ2= α
2kΩ230− 1, Δ3 = α2
2kΩ230− 1 2(2k + 1)Ω230− 1. (4.16) The requirement that all these determinants be greater than zero gives Ω230> 1/(2k). Thus, for a
uniform disk with k = 1/4 asymptotic rolling solutions may be linearly stable only when Ω230> 2
and for a ring (k = 1/2) when Ω230> 1. Notice that this condition is independent of the value of
the parameter α.
For spinning solutions the characteristic polynomial of the Jacobian matrix has the form
p2(λ) = λ 2k2q2(λ), q2(λ) = 2k2λ5+ k(4k + 3)αλ4+ [α2+ k 4ϕ02k + α2(2k + 3)− 2 ]λ3+ α 4k ϕ02(k + 1) −1) − 1) λ2+ α2(2k + 1)ϕ 0 2 (k + 1)− 1 + 2ϕ02k ϕ02k− 1 λ + αϕ02 ϕ02k− 1 . (4.17)
We have a simple zero eigenvalue and the Routh – Hurwitz criterion applied to the fifth-order polynomial q2(λ) gives the conditions
Δ1 = αk(4k + 3) > 0, Δ2 = αk 4k ϕ02k(2k + 1)− 1 + α2(k + 1)(2k + 1)(4k + 3) > 0, Δ3 = α2k α2(2k + 1)(4k + 3) k ϕ02(k + 1)− 2 − 1+ 4k k ϕ04k(2k + 1)− 2ϕ02+ 4 + 1 > 0, Δ4= α2kC > 0, Δ5 = α3ϕ0 2 k ϕ02k− 1 C > 0, where C = 32ϕ02k3 ϕ02k− 1 − α4(2k + 1)2(4k + 3)2kϕ 0 2 (k + 1)− 1 − 1 − 4α2k(2k + 1)k2ϕ 0 4 k(2k + 1)− ϕ02(8k + 9) + 4 + 1 . (4.18) The regions, in a plane α, (ϕ0)2, corresponding to linearly stable asymptotic spinning solutions for a uniform disk and for a hoop are shown in Figs. 3a and 3b, respectively.
(a) for uniform disk (b) for hoop
Fig. 3. Regions of linear stability for spinning solutions.
A characteristic polynomial of the Jacobian matrix for tumbling solutions is the same for both choices of signs in the matrix Ltumb entries and it also has a factor λ
p3(λ) = λq3(λ) 2κ2k2(−κ2+ k + 1)2, q3(λ) = 2κ2k2 −κ2+ k + 12λ5+ ακ2k−κ2+ k + 1 4k2− κ2(2k + 1) + 5k + 1λ4 + κk−κ2+ k + 1 α2κ(k + 1)(2k + 1) + 2κ2(5k− 1) − 8k + 2λ3+ ακ − 2κ2 + κ4(1− k(4k + 7)) + κ2k(k(20k + 31) + 10)− k(4k(5k + 6) + 3) + 1 λ2 + 10kκ2(k + 2)− k − 1+ α2κ56k2+ k− 1+ α2κ3(2k + 1)(k(5k + 3) + 1) − 5α2κk(k + 1)(2k + 1)− 10κ4k− ακ2− 1 κ2(16k + 5)− 5(k + 1). (4.19)
For tumbling solutions the Routh – Hurwitz conditions applied to the fifth-order polynomial
q3(λ) give the following inequalities: Δ1 = ακ2k −κ2+ k + 1 4k2− κ2(2k + 1) + 5k + 1> 0, Δ2 = ακ3k2 −κ2+ k + 12α2κ(k + 1)28k2+ 6k + 1− κ3(k + 1)(α + 2αk)2 + 4κ4(2− 3k)k − 4κ2k(k + 4) + 8k(k + 1)2 > 0, Δ3 = α2κ4 −k2 −κ2+ k + 12 × α2κ5(k + 1)(2k + 1)kk24k2+ 58k + 7− 21− 6− α2κ(k + 1)48k2+ 6k + 1 + 4κ6kkk60k2+ 51k− 55− 52+ 8− 4κ4kkk80k2+ 6k− 119− 36+ 12 + α2κ9(2k + 1)3(3k− 1) − α2κ3(k + 1)2(2k + 1)(k(2k(8k− 3) − 17) − 4) − κ7(k(2k(9k + 11)− 3) − 4)(α + 2αk)2− 4κ8k(k(k(12k + 7)− 20) + 2) + 4κ2k(k + 1)(k(k(16k− 41) − 4) + 8) − 8k(k + 1)4 > 0, Δ4=−α2k2κ6D > 0, Δ5 = α3(κ− 1)κ6(κ + 1)k2 κ2(16k + 5)− 5(k + 1)D > 0, (4.20) where D =−κ2+ k + 12−6k2+ κ4(3k− 1) + κ2(k(5k + 4) + 2)− 7k − 1 κ432k3(k(15k + 19)− 39) + α4(k + 1)(2k + 1)2(k(4k + 3)(k(6k + 13)− 2) − 3) + 32k3(k + 1)(2k− 13) − 4α2κ7k(2k + 1)k12k2+ k− 23+ 2− 4α2κk(k + 1)(2k + 1) ×k4k2+ 66k + 15− 2+ κ632(13− 15k)k3− 3α4(2k + 1)36k3+ 8k2− 1 + κ232k3((7− 17k)k + 39) − α4(k + 1)2(2k + 1)2(4k + 1)6k2+ 2k− 1+ 4α2κ5k × (2k + 1)(k(3k(k(20k + 7) − 31) − 59) + 6) + 4α2κ3k(k(k(k(2(39− 56k)k + 417) + 273) + 37)− 6) + κ8(3k− 1)(α + 2αk)4 .
(a) for uniform disk (b) for hoop
Fig. 4. Regions of linear stability for tumbling solutions (κ = sin θ0).
The regions, in a plane (α, κ = sin θ0), corresponding to (possibly) linearly stable asymptotic
tumbling solutions for a uniform disk and for a hoop are shown in Figs. 4a and 4b, respectively. Since these regions are qualitatively different, we considered an annulus with external radius R and internal radius r. Its principal moments of inertia are equal to
I1 = I2 = m 4(R 2+ r2) = kmR2, I 3 = 2kmR2, k = 1 4 1 + r R 2 ,
where the internal radius r∈ [0, R], which gives k ∈14,12. The limiting cases r = 0 correspond to a solid disk with minimal value k = 1/4 and r = R corresponds to a hoop with maximal value of
k = 1/2. We notice that already very small perturbations of the asymptotic case of the maximal
value of the principal moments of inertia change qualitatively the region of linear stability, see Fig. 5a, where k = 12 −10001 corresponding to r/R = 0.997998 was chosen. In general the region of linear stability as a function of parameters (α, κ, k) with k∈14,12 is given in Fig. 5b. On the upper surface k = 12 we recognize the region 4b for a hoop and on the bottom surface k = 14 we observe the region 4a for a solid disk.
(a) for a very thin annulus with
k = 12−10001
(b) general case
Fig. 5. Regions of linear stability for tumbling solutions for a) a very thin annulus with k = 1 2−
1
1000 and b)
general case with k∈ [1 4,
1
2]. Here κ = sin θ0. Taking into account that
˙ ϕ = g Rϕ , ω3 = g RΩ3,
the necessary conditions of linear stability for rolling and spinning solutions are
ω32 > 1 2k g R, ϕ˙ 2> 1 k g R, (4.21)
respectively, and, moreover, for spinning solutions α must belong to the region given in Fig. 3. It should be made clear that the linear analysis has determined only possible regions of stability for all three types of asymptotic solutions. The criterion is not conclusive because the characteristic polynomial for each Jacobian matrix (4.9), (4.10) or (4.11) admits one or two zero eigenvalues. However, it determines regions of possible stability which we can relate to the behavior of solutions studied numerically in the next section.
5. NUMERICAL SIMULATIONS
The knowledge of asymptotic solutions and of their possible stability regions provides a framework for understanding and interpreting numerical simulations. In this section we simulate the dynamics of the sliding and rolling disk with parameters m = 0.02 kg, R = 0.02 m, μ = 0.3, I3 =
1
2mR
2, I
1= 14mR2, and g = 9.82 m/s2. They correspond to the value α = 3 √
491/500≈ 0.133 m/s. Taking into account that
g R =
√
491 s−1 and inequalities (4.21), the necessary conditions for linear stability of asymptotic solutions for our disk are the following:
• for rolling solutions: |ω30| > ω3tresh=√2√491≈ 31.337 rad/s,
• for spinning solutions: |ϕ0| > ϕ0tresh= 2√491≈ 44.31704 rad/s,
• for tumbling solutions: if α 1 the only condition is κ = sin θ0> 0, for α > 1 the threshold
curve is given by a solution of
3α4κ29κ2− 2κ2+ 5+ α2κ8κ4− 30κ2+ 34+ 8κ2κ2− 3+ 16 = 0, with κ > 0 for α > 0, see the dashed curve in Fig. 4a.
To learn about the behavior of solutions, we take reference initial conditions in a neighborhood of rolling or spinning asymptotic solutions and then we study what happens under variation of certain initial variables. In order to prevent the possibility of bouncing in our simulations, we bound the maximal values of the rolling and spinning angular velocities ω3(0) 130 rad/s, ˙ϕ(0) 80 rad/s so
that the value gn(t) > 0 of the reaction force is positive all time and goes in oscillatory way to the
weight mg of the disk (see Fig. 8b).
As a first reference initial condition we choose the values of variables from a neighborhood of an asymptotic rolling solution
ICroll(0) = θ(0) = π 2 − 0.1 rad, ˙θ(0) = ˙ϕ(0) = 0, ω3(0) = 90 rad/s, vx(0) = vy(0) = 0 m/s , (5.1) where the disk is tilted from the vertical by 0.1 rad and the only nonzero variable is the rolling angular velocity ω3. This choice of initial condition is motivated by the following observations:
a) all asymptotic solutions have ˙θ(∞) = vx(∞) = vy(∞) = 0,
b) usually a sliding and rolling disk is launched close to θ≈ π/2 ≈ 1.57 rad,
c) time evolution of ICroll(0) having energy about 8 times larger than the threshold energy Etresh= 12I1ϕ˙2tresh with ˙ϕtresh=
mgR
I1 illustrates how energy is dissipated by the disk approaching an asymptotic solution.
Figures 6 and 7 show that the time evolution consists of two distinct stages: dissipative rolling and sliding up to the critical time Tcr= 15 s when ω3 decreases from 90 to 0 rad/s. It ends up
with a drastic jump of the inclination angle θ from the initial value θ(0) = π/s− 0.1 rad to the neighborhood of the asymptotic tumbling angle θ(∞) ≈ 0.61 rad, see Fig. 6a. Also, the angular velocity ˙ϕ jumps from ˙ϕ(0) = 0 rad/s to the neighborhood of the asymptotic value
˙
ϕ(∞) = −
gmR
I1sin θ(∞) ≈ −58.55 rad/s,
see Fig. 6b and formula (3.1). An oscillatory approach to the asymptotic solution takes place after
Tcr= 15 s, while the amplitude of oscillations attenuates to zero. Figures 6 and 7 indicate that
after 25 s the disk gets close to the asymptotic tumbling state
Stumb =
θ(∞) ≈ 0.61 rad, ˙ϕ(0) = −58.55 rad/s, ˙θ(0) = ω3(0) = 0 rad/s, vx(0) = vy(0) = 0 m/s
.
(5.2) The initial value of the angular velocity ω3(0) = 90 rad/s decreases almost monotonically to
the asymptotic value ω3(∞) = 0 rad/s characterizing all tumbling asymptotic solutions. It is
accomplished at t = Tcr= 15 s, see Fig. 7a. Notice that for t < Tcr the angular velocity ˙ϕ remains
almost constant, while ˙ψ is decreasing from the initial value ˙ψ(0) = 90 rad/s to about ˙ψ = 20 rad/s.
For t > Tcrthe oscillations of ˙ψ and ˙ϕ match each other to keep ω3(t) = ˙ψ(t) + ˙ϕ cos θ(t) = 0 during
the final stabilization phase. Asymptotically ω3(∞) = 0, thus ˙ψ(∞) = − cos θ(∞) ˙ϕ(∞) = 48.12
(a) evolution of θ(t)
(b) evolution of ˙ϕ(t) and ˙ψ(t)
Fig. 6. Time evolutions of θ, ˙ϕ and ˙ψ. The plots were obtained from numerical integration of equations of
motion with ICroll(0).
Components of velocity that start with initial values vx(0) = vy(0) = 0 m/s acquire the positive
maximal value vx(t) = vy(t) = 0.4 m/s for times t < Tcr. At critical time Tcr they jump to a
neighborhood of its asymptotic values vx(∞) = vy(∞) = 0 m/s and oscillate with amplitudes going
to zero. The transversal velocity vx for most of the time interval [0, Tcr] takes the value vx(t)≈ 0.2
m/s, while vy(t)≈ 0.02 m/s is an order of magnitude smaller. Thus, it is the friction in the direction
(a) evolution of ω3(t)
(b) evolution of vx(t) and vy(t)
Fig. 7. Time evolutions of ω3, vxand vy. The plots were obtained from numerical integration of equations of
motion with ICroll(0).
In Fig. 8a the total energy decreases monotonically in the interval [0, Tcr]. After Tcr all energies
remain close to its asymptotic values: E(∞) = 0.00346 J for the total energy, Etrans(∞) = 0 for the energy of translational motion, Epot(∞) = 0.002306 J for potential energy, and Erot(∞) =
0.001153 J for the energy of rotational motion. But due to the oscillatory approach to the asymptotic state we see in Fig. 8a an oscillatory exchange of energy between the rotational and potential modes as the center of mass is wobbling up and down.
(a) evolution of total energy E(t) and its components
(b) evolution of the value of the normal force gn z
Fig. 8. Time evolution of a) total energy E(t) and its components Etrans(t), Erot(t), Epot(t) and b) the normal
force gn z which is positive at all times. Left y axis and bottom x axis concern ICspin(0) and ICmixed(0). The
right y axis and the top x axis concern ICroll(0).
Figures 6, 7 and 8a show that when the disk starts its time evolution with initial conditions
ICroll(0), then at the beginning up to Tcr its evolution takes place in a neighborhood of an
asymptotic rolling solution θ(∞) = π/2. The inclination angle stays close to π/2, the angular velocity is close to ˙ϕ(t) = 0 for t∈ [0, Tcr) and the value of ω3 decreases, but the solution has
still a rolling character. It is a stage when it dissipates most of excess energy, which initially has been about 15 times higher than the energy of the asymptotic solution Stumb given by (5.2).
At t = Tcrthe solution switches to a different regime of dynamics close to an asymptotic tumbling
solution and during final oscillations its energy approaches the energy of Stumb. Similar behavior
can be observed for many other values of the initial velocity ω3(0), see Fig. 9. For some values of ω3(0) the critical time is large and θ(t) may strongly oscillate.
Fig. 9. Curves of θ(t) escaping from a neighborhood of the asymptotic rolling solution θ(∞) = π/2 to certain
asymptotic tumbling solutions for chosen initial ω3(0), with remaining initial conditions as in ICroll(0). The
threshold value is ω3tresh=
√
2√491≈ 31.337 rad/s.
Interesting is the dependence of the asymptotic inclination angle θ(∞) on the initial value ω3(0) when keeping remaining initial conditions as in ICroll(0). The graph presented in Fig. 10 is quite
complicated:
• for 0 < ω3(0) < 19 rad/s the asymptotic inclination angle is growing up to θ(∞) → 0.67 rad. After reaching a maximum, in the interval 19 ω3(0) < 23 rad/s the graph is decreasing to θ(∞) = 0.47 rad,
• for 23 ω3(0) < 45 rad/s the asymptotic angle θ(∞) belongs to almost the whole range
θ(∞) ∈ [0, π]. The dependence θ(∞) (ω3(0)) shows great sensitivity to small changes of ω3(0),
as in Fig. 10 and in the magnification in Fig. 11(a).
Time evolutions of θ(t) with very close initial values of ω3(0) = 40.07731423 rad/s and ω3(0) =
40.077314235 rad/s that lie on opposite sides of one jump of the function θ(∞) = θ(∞) (ω3(0)) are presented in Fig. 11b. A very small change from ω3(0) = 40.07731423 rad/s to ω3(0) =
40.077314235 rad/s causes the inclination angles to evolve for about 3 seconds in the same way, but later both curves separate and reach different asymptotic values: θ(∞) = 2.65 rad and
θ(∞) = 0.32 rad, respectively.
This sensitivity shows that despite the fact that the disk is initially tilted to one side θ(0) =
π/2− 0.1 rad, a small change of initial ω3(0) may cause tipping of the sliding rolling disk onto the other side. This raises a question whether this behavior is not an artefact of a sharp edge of the disk.
• surprisingly starting from ω3(0) = 45 rad/s, the asymptotic angle is almost constant θ(∞) = 0.62 rad for the whole tested range 45 < ω3(0) < 130 rad/s. The time of rolling, before the
Fig. 10. Asymptotic angle θ(∞) as a function of initial condition ω3(0) with the remaining initial conditions
as in ICroll(0).
initial kinetic energy of a rolling and sliding disk has to be dissipated and it takes longer time as the average and maximal value of the transversal velocity vx(t) varies little for ω3(0)
in the tested interval 45 < ω3(0) < 130 rad/s.
Regions of large sensitivity of the asymptotic inclination angle θ(∞) to small changes of initial conditions are also visible in Fig. 12, where we keep ω3(0) = 90 rad/s but vary θ(0) and ˙ϕ(0). In all
these figures asymptotic values of θ(∞) are obtained by numerical integration until ω3(t) reaches the asymptotic value ω3(∞) = 0, with accuracy 10−6 in two successive steps.
In Fig. 13 we present the time evolution of an initially vertical disk with θ(0) = π2 rad, ˙
θ(0) = 10 rad/s for two different values of the rolling velocity ω3(0) = 90 > ω3tresh= √
2√491≈ 31.337 rad/s > ω3(0) = 20; one above and the other below the stability threshold value ω3tresh.
The solution with ω3(0) = 90 rad/s remains vertical, see in Fig. 13a within a simulation time of
20 s in the picture, while the other one with ω3(0) = 20 rad/s falls to the tumbling solution with ω3(∞) = 0 rad/s and θ(∞) = 0.371 rad/s, which is related to ˙ϕ(∞) by relation (3.1). In Fig. 13b
we see that ˙ϕ(t)≈ 0 for the solution ω3(0) = 90 and falls to the value ˙ϕ(∞) = −73.516 rad/s
for the solution ω3(0) = 20. It is also visible in Fig. 13b that ω3(t) = 90 rad/s remains constant
for the first solution. Also, sufficiently small perturbations of vy(0)= 0 do not cause escaping of
trajectories from the neighborhoods of rolling solutions in agreement with the necessary linear stability condition for rolling asymptotic solutions.
We do not plot graphs for initial conditions in a neighborhood of certain tumbling solutions because they qualitatively behave similarly as solutions with the initial condition ICroll(0) after t > Tcr. For our disk with α≈ 0.133 1 the conditions for linear stability of tumbling solutions are satisfied for all values of the asymptotic inclination angle θ(∞) and numerical experiments show that when we choose some θ(0) and calculate ˙ϕ according to (3.1) and disturb this value,
then these solutions immediately converge to a stable asymptotic tumbling solution. However, for a disk with α > 1 for sufficiently large θ(0) tumbling solutions are unstable and escape to other tumbling solutions contained in the linear stability region presented in Fig. 4a. This is visible in Fig. 14 presenting simulations for a disk with μ = 4.513 that corresponds to α = 2. A threshold value dividing stable and unstable regions is θtresh= 0.8553 rad. Calculations in Fig. 14 are made
(a) sensitivity interval
(b) time evolutions of θ(t)
Fig. 11. Magnification of the interval of great sensitivity to small changes of ω3(0) and time evolution of θ(t)
with two very close initial ω3(0) = 40.07731423 rad/s (upper curve) and ω3(0) = 40.077314235 rad/s (lower
curve).
(1, 48.216), (1.20, 45.8043), (1.45, 44.3794), (1.57, 44.217), (1.75, 44.5762). They show that tumbling solutions with θ(0) θtresh keep the angle constant all time.
A second natural class of initial conditions to study are spinning solutions with θ(0) = π/2− 0.1 being close to the vertical spinning solution and with exchanged initial angular velocities so that
˙
ϕ(0) = 80 rad/s and ω3(0) = 0 rad/s. As a spinning reference solution we take ICspin(0) =
θ(0) = π/2− 0.1 rad, ˙ϕ(0)=80 rad/s, ˙θ(0)=ω3(0) = 0 rad/s, vx(0) = vy(0) = 0 m/s
.
(a) dependence of θ(∞) on θ(0)
(b) dependence of θ(∞) on ˙ϕ(0)
Fig. 12. Asymptotic angle θ(∞) as a function of initial conditions θ(0) and ˙ϕ(0) with remaining initial
conditions as in ICroll(0).
This solution with θ(0) = π/2− 0.1 = 1.47 rises to a vertically spinning solution with θ(∞) =
π/2 as shown by the θ curve in Fig. 15. The angular velocity ˙ϕ remains almost constant very close
to ˙ϕ(0) = 80 rad/s as illustrated by the upper ˙ϕ curve plotted in Fig. 15 with respect to the right y axis. The other variables ˙ψ(t), ω3(t), vx(t), vy(t) take very small values that quickly go to zero.
This behavior is a feature of all solutions with θ(0) = π/2− 0.1 rad and ˙ϕ(0) > ˙ϕ0tresh=
(a) evolution of θ(t)
(b) evolution of ˙ϕ(t) and ω3(t)
Fig. 13. Time evolution of the inclination angle θ(t) and of angular velocities ˙ϕ(t) and ω3(t) for initial
conditions θ(0) = π
2 rad, ˙θ(0) = 10 rad/s and two different values of ω3(0) = 90 and 20 rad/s with zero
remaining initial values. In figure (b) the left y axis provides a scale for ˙ϕ and the right y axis provides
a scale for ω3.
and ˙ϕ(∞), which are (again) related by
˙ ϕ(∞) = mgR I1sin θ(∞) .
Figure 16a presents the time evolution for several initial values of ˙ϕ(0) = 80, 45, 2√491, 40 and 30 rad/s, while the remaining initial conditions are the same as in ICspin(0). As expected from
Fig. 14. Time evolutions of tumbling solutions for a disk with α = 2.
Fig. 15. Time evolutions of θ(t), ˙ϕ(t) and ˙ψ(t) from numerical integration of equations of motion with ICspin(0). The left y axis provides a scale for θ and the right y axis provides a scale for angular velocities ˙ϕ
and ˙ψ.
the stability analysis, the first two solutions approach the asymptotic spinning solution and the remaining three with ˙ϕ(0) ˙ϕ0tresh = 2
√
491 rad/s approach certain tumbling solutions with decreasing asymptotic inclination angle θ(∞) as ˙ϕ(0) is decreasing. The dependence of θ(∞) as a function of the initial angular momentum in ICspin(0) for ˙ϕ(0)∈ [−80, 80] rad/s is shown in
Fig. 16b. Regions of asymptotic spinning solutions and of tumbling solutions separated by the threshold value| ˙ϕ0tresh| = 2√491 rad/s are clearly visible in the graph.
(a) graphs of θ(t) for chosen values of ˙ϕ(0)
(b) dependence of θ(∞) on ˙ϕ(0)
Fig. 16. Time evolution of θ(t) for a sequence of values ˙ϕ(0) = 80, 45, 2√491, 40 and 30 rad/s and dependence of the asymptotic θ(∞) on ˙ϕ(0) ∈ [−80, 80] rad/s, with remaining initial conditions as in ICspin(0).
Perturbations of the initial value ˙θ(0) = 0, vx(0) = 0 and vy(0) = 0 give similar behavior of
dynamical variables as with the initial conditions ICspin(0), but the time of approaching an
asymptotic solution is longer since there is more energy to be dissipated.
What about solutions with both ˙ϕ(0)= 0 and ω3(0)= 0? Here the most interesting question is which type of stable asymptotic solution is approached depending on the values of ˙ϕ(0) and ω3(0). It seems that there is no any clear answer. As an example we give an interesting solution that
behaves in a nonoscillatory way
ICmixed(0) =
θ(0) = π/2− 0.35 rad, ˙ϕ(0)=50 rad/s, ˙θ(0)=ω3(0) = 0 rad/s, vx(0) = vy(0) = 0 m/s
,
(5.4) which is going to θ(∞) = 1.5 rad, see Figs. 17 and 18.
Fig. 17. Time evolutions of θ(t), ˙ϕ(t) and ˙ψ(t) from numerical integration of equations of motion with ICmixed(0). The left vertical y axis provides a scale for angle θ and the right vertical y axis provides a scale
for angular velocities ˙ϕ and ˙ψ.
Fig. 18. Time evolutions of ω3(t), vx(t) and vy(t) from numerical integration of equations of motion with
initial conditions ICmixed(0). The left vertical y axis provides scale for angular velocity ω3and the right vertical
Equations of motion were integrated numerically using a double precision version of the Bulrisch – Stoer extrapolation method with adaptive step-size control [23]. The local relative precision of integration was fixed to 10−12for most of the tests. It appeared that for some values of the parameters and initial conditions the system is stiff. However, comparison of results generated by the Bulrisch – Stoer integrator with numerical methods dedicated to integration of stiff problems shows differences within precision limits.
An analogous linear stability analysis, as in Section 4, made for the straight rolling and spinning motions of a purely rolling disk (see condition (2.5)) gives the well-known classical values of stability [25, §§242–244], see also [19, 22]. The results are summarized in Table 1 to see that the stability thresholds for the model with friction are higher than those for pure rolling. This is probably due to difference in the reaction force in both cases.
Table 1. Comparison of necessary stability conditions. Solution rolling & sliding pure rolling spinning ϕ˙2> 1 k g R ϕ˙ 2> 1 k + 1 g R rolling ω23> 1 2k g R ω 2 3 > 1 2(2k + 1) g R
Finally, let us check how the friction influences the stability of asymptotic solutions in our model. Notice that all results in Section 4 have been obtained under the assumption μ= 0 used in our rescaling, see Eq. (4.1) with T = R/g and L = R/g/μ. The case μ = 0 corresponds to the
motion of a disk on an absolutely smooth plane. Repeating linear stability analysis, as in Section 4 for Eqs. (2.13)–(2.17) with μ = 0, we have obtained the following results:
1. The spinning solution is linearly stable if ˙ϕ2 > 1kRg.
2. The rolling solution is linearly unstable. The matrix of the linearized system has the multiple eigenvalue λ = 0 with a two dimensional Jordan block.
3. The tumbling solution is linearly stable for all θ0 ∈ (0, π).
We notice that threshold for spinning solution agrees with the value for our model and results in [22]. For rolling solutions the condition that real parts of nonvanishing eigenvalues of linearization matrix are negative gives inequality ω32 > 4k1 Rg that agrees with the condition in [22]. But perturbations in the variables (θ, ˙θ) produce a two-dimensional Jordan block corresponding to
the multiple eigenvalue λ = 0. The presence of a Jordan block of dimension greater than one for an eigenvalue with a vanishing real part implies linear instability, see, e.g., Theorem 1 on page 86 in [7]. Nonlinear stability analysis for the motion of a disk on an absolutely smooth plane is given in [18, Chap. 2].
In the simulations we have tested angular momentum velocities in the range of 0 < ω3(0) <
130 rad/s or 0 < ˙ϕ < 80 rad/s to keep energy sufficiently low for having gn(t) > 0 all time and for
avoiding the possibility of bouncing solutions. In Fig. 8b we show the evolution of gn(t) for our
three basic reference initial conditions. Note that two y axes are used: the left one for solutions with initial conditions ICspin(0) and ICmixed(0) and the right one for solutions with ICroll(0). For
all these solutions the asymptotic value of gn(t) is the same and equal to mg. For this reason we
observe no loss of contact of the disk with the supporting plane such as may occur in the case of the Euler disk [3, 11].
In order to keep the results of numerical simulations understandable, we have chosen reference initial conditions ICroll(0), ICspin(0) close to linearly stable vertical rolling and spinning asymptotic
solutions by tilting the initial inclination angle θ(0) by 0.1 rad and having only either ω3(0) =
to other variables under the action of the friction force. As expected, the solutions acquire nonzero sliding velocity with the transversal component vx(t) being an order of magnitude larger than
the longitudinal vy(t) speed. The reference solution ICroll(0) tends asymptotically to a tumbling
solution with the inclination angle θ(∞) = 0.62 rad being (surprisingly) almost constant over the whole range 45 rad/s < ω3(0) < 130 rad/s. Quite unexpectedly the value of θ(∞) shows great
sensitivity to the value of the initial angular velocity ω3(0) in the range 25 rad/s < ω3(0) < 45 rad/s
and is given by a regular, increasing continuous function of ω3(0) in the range 0 < ω3(0) < 25 rad/s.
It would be impossible to deduce this picture using analytical methods. The asymptotic rolling solutions with ω3(0) > ωtresh=
√
2√491 rad/s behave stable in numerical simulations only for small perturbations. Otherwise they fall into tumbling solutions.
Greater stability regions have the spinning solutions as the chosen reference spinning solution
ICspin(0) starting with the same initial inclination angle θ(0) = π/2− 0.1 = 1.47 rad as for ICroll(0) goes asymptotically to a stable spinning solution with almost the same ˙ϕ(∞) ≈ 80 rad/s. The curve
16(b) for an asymptotic inclination angle θ(∞) agrees with the threshold value ˙ϕ(∞) = 2√491 rad/s for the stability of spinning solutions, but the shape of the dependence of θ(∞) on ˙ϕ(0) cannot be found analytically. The tumbling solutions appear to be stable in simulations for all values of
θ(∞) ∈]0, π/2] and small changes of initial conditions. This is in agreement with the linear stability
analysis. The tumbling solutions seem to play a central role as they are attractors for solutions starting far away from the stable rolling and spinning solutions and for solutions with low initial energy starting even close to rolling or spinning solutions.
ACKNOWLEDGMENTS
M.P. has been supported by grant No. DEC-2013/09/B/ST1/04130 of the National Science Center of Poland, and M.P. and S.R. gracefully acknowledge support from the Department of Mathematics of Link¨oping University.
The authors express their gratitude to Andrzej J. Maciejewski for his valuable comments and for his help with numerical simulations.
REFERENCES
1. Appell, P., Sur l’int´egration des ´equations du mouvement d’un corps pesant de r´evolution roulant par une arˆete circulaire sur un plan horizontal: cas particulier du cerceau, Rend. Circ. Mat. Palermo, 1900, vol. 14, no. 1, pp. 1–6.
2. Borisov, A. V., Mamaev, I. S., and Kilin, A. A., Dynamics of Rolling Disk, Regul. Chaotic Dyn., 2003, vol. 8, no. 2, pp. 201–212.
3. Borisov, A. V., Mamaev, I. S., and Karavaev, Yu. L., On the Loss of Contact of the Euler Disk, Nonlinear
Dynam., 2015, vol. 79, no. 4, pp. 2287–2294.
4. Chaplygin, S. A., On a Motion of a Heavy Body of Revolution on a Horizontal Plane, Regul. Chaotic
Dyn., 2002, vol. 7, no. 2, pp. 119–130; see also: Collected Works: Vol. 1, Moscow: Gostekhizdat, 1948,
pp. 51–57.
5. Cohen, C. M., The Tippe Top Revisited, Am. J. Phys., 1977, vol. 45, no. 1, pp. 12–17.
6. Contensou, P., Couplage entre frottement de pivotement et frottement de pivotement dans la th´eorie de latoupie, in Kreiselprobleme Gyrodynamics: IUTAM Symp. Celerina, Berlin: Springer, 1963, pp. 201– 216.
7. Demidovich, B. P., Lectures on the Mathematical Stability Theory, Moscow: Nauka, 1967 (Russian). 8. Ebenfeld, S. and Scheck, F., A New Analysis of the Tippe Top: Asymptotic States and Liapunov Stability,
Ann. Physics, 1995, vol. 243, no. 2, pp. 195–217.
9. Gallop, E. G., On the Rise of a Spinning Top, Proc. Camb. Phylos. Soc., 1904, vol. 19, no. 3, pp. 356–373. 10. Gantmacher, F. R., The Theory of Matrices: Vol. 2, New York: Chelsea, 1959.
11. Ivanov, A. P., On Detachment Conditions in the Problem on the Motion of a Rigid Body on a Rough Plane, Regul. Chaotic Dyn., 2008, vol. 13, no. 4, pp. 355–368.
12. Karapetian, A. V., On the Regular Precession of a Body of Revolution on a Horizontal Plane with Friction, J. Appl. Math. Mech., 1982, vol. 46, no. 4, pp. 450–453; see also: Prikl. Mat. Mekh., 1982, vol. 46, no. 4, pp. 568–572.
13. Karapetian, A. V., Global Qualitative Analysis of Tippe Top Dynamics, Mech. Solids, 2008, vol. 43, no. 3, pp. 342–348; see also: Izv. Akad. Nauk SSSR. Mekh. Tverd. Tela, 2008, no. 3, pp. 33–41.
14. Karapetian, A. V., A Two-Parameter Friction Model, J. Appl. Math. Mech., 2009, vol. 73, no. 4, pp. 367– 370; see also: Prikl. Mat. Mekh., 2009, vol. 73, no. 4, pp. 515–519.
15. Korteweg, D. J., Extrait d’une lettre `a M. Appell, Rend. Circ. Mat. Palermo, 1900, vol. 14, no. 1, pp. 7–8. 16. Leine, R. I., Le Saux, C., and Glocker, Ch., Friction Models for the Rolling Disk, in Proc. of the ENOC
Conf. (Eindhoven, Netherlands, 7–12 August, 2005), 10 pp.
17. Le Saux, C., Leine, R. I., and Glocker, Ch., Dynamics of a Rolling Disk in the Presence of Dry Friction,
J. Nonlinear Sci., 2005, vol. 15, no. 1, pp. 27–61.
18. Markeev, A. P., Dynamics of a Body Being Contiguous to a Rigid Surface, Moscow: Nauka, 1992 (Russian).
19. McDonald, A. J. and McDonald, K. T., The Rolling Motion of a Disk on a Horizontal Plane, arXiv:physics/0008227 (2000), 21 pp.
20. Moshchuk, N. K., Motion of a Heavy Rigid Body on a Horizontal Plane with Viscous Friction, J. Appl.
Math. Mech., 1985, vol. 49, no. 1, pp. 49–53; see also: Prikl. Mat. Mekh., 1985, vol. 49, no. 1, pp. 66–71.
21. Munitsyna, M. A., The Motions of a Spheroid on a Horizontal Plane with Viscous Friction, J. Appl. Math.
Mech., 2012, vol. 76, no. 2, pp. 154–161; see also: Prikl. Mat. Mekh., 2012, vol. 76, no. 2, pp. 214–223.
22. O’Reilly, O. M., The Dynamics of Rolling Disks and Sliding Disks, Nonlinear Dynam., 1996, vol. 10, no. 3, pp. 287–305.
23. Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P., Numerical Recipes: The Art of
Scientific Computing, 3rd ed., Cambridge: Cambridge Univ. Press, 2007.
24. Rauch-Wojciechowski, S. and Rutstam, N., Dynamics of the Tippe Top: Properties of Numerical Solutions Versus the Dynamical Equations, Regul. Chaotic Dyn., 2013, vol. 18, no. 4, pp. 453–467. 25. Routh, E. J., The Advanced Part of a Treatise on the Dynamics of a System of Rigid Bodies: Being
Part II of a Treatise on the Whole Subject, 6th ed., New York: Dover, 1955.
26. Zhuravlev, V. F., On a Model of Dry Friction in the Problem of the Rolling of Rigid Bodies, J. Appl.
Math. Mech., 1998, vol. 62, no. 5, pp. 705–710; see also: Prikl. Mat. Mekh., 1998, vol. 62, no. 5, pp. 762–
