- Email: [email protected]
THE DRIVEN ROLLING TORUS
THE DRIVEN ROLLING TORUS John Hauser ∗ Alessandro Saccon ∗∗ ∗ Electrical and Computer Engineering University of Colorado, Boulder, CO 80309-0425 [email protected] ∗∗ Dipartimento di Ingegneria dell’Informazione Universit` a di Padova, 35131 Padova, Italy [email protected]
Abstract: We develop the equations of motion for a torus with a general contact profile that rolls without slipping and that is actuated using drive and steer torques. The form of the dynamics plus nonholonomic constaints allows an easy reduction to a dynamic control system of dimension 4 including the three angular rates and the roll (or lean) angle. Using a standard torus with geometric and inertial parameters similar to those of a motorcycle tire and wheel, we find that the rolling torus exhibits, with increasing forward speed, an unexpected transition c from stability to instability and back to stability! Copyright °2007 IFAC Keywords: Nonlinear control systems, Lagrangian dynamics, reduction, integrability, motorcycle tires, instability.
We develop the equations of motion for a torus with a general contact profile that rolls without slipping and that is actuated using drive and steer torques that are somehow magically applied. The purpose of this study is to work toward a better understand of the dynamics and control of motorcycles operating in high traction/low slip conditions. Proceeding carefully to encourage the discovery of key kinematic and dynamic features, we find that, for a general class of contact surface (or tire) profiles, the form of the dynamics plus nonholonomic constaints allows an easy reduction to a dynamic control system of dimension 4 including the three angular rates and the roll (or lean) angle. We show, by elementary means, that the undriven system admits to classical integrability so that the pitch and yaw rates can be given as (instantaneous) functions of the roll angle. Using this, together with conservation of energy, we find that the stability of steady motions
is easy to study and, furthermore, that almost every (undriven) motion will be periodic. For each roll angle, steady turning is only possible when the forward velocity (or wheel pitch rate) is at or above a critical value determined by a simple equation, quadratic in the yaw rate, required for equilibrium. Above that value, there are two equilibria. Using a standard torus with geometric and inertial parameters similar to those of a motorcycle tire and wheel, we demonstrate that dynamics can be strange and intuition lacking. Indeed, for the parameters chosen and a range of steady roll angles, the rolling torus exhibits, with increasing forward speed, a transition from stability to instability and back to stability! The study of rigid bodies subject to constraints, especially nonholonomic constraints, has a long and exciting history and continues to this day, see, e.g., (Bloch, 2003), (Bullo and Lewis, 2004), and (Arnold et al., 2006), which include many
classical and modern references. Indeed, (classical) integrability and the ubiquity of periodic solutions for rolling disks and bodies of revolution were established in the late 1800’s and early 1900’s, with special emphasis on expressing these solutions using, e.g., hypergeometric or Legendre functions. Though the dynamics of a rolling disk are quite similar, the bifurcation analysis in (O’Reilly, 1996) does not indicate the presence of the stability-instability-stability transition that we observe in our numerical example. Generalized coordinates for the rolling torus For the purpose of deriving the equations of motion, we consider the system as a rigid body subject to the constraint of rolling without slipping. We will see that, since the region of operation is restricted lean angles less than 90 degrees, there is a convenient set of generalized coordinates q ∈ R6 parametrizing the configuration of the rolling torus, enabling the (global) determination of the rolling torus dynamics using the Lagrangian method. In this approach, the rolling contact constraint is enforced by contact point forces acting as Lagrange multipliers is the system dynamic equations. Our goal, then, is to establish this set of generalized coordinates, express the kinetic and potential energies, T and V , in terms of the generalized coordinates and velocities, and determine the (mostly nonholonomic) constraints with respect to same. Let x ∈ R3 and R ∈ SO(3) denote the position and orientation of a frame fixed at the wheel center of mass with respect to a fixed inertial (or spatial ) frame with x-y-z axes oriented in a north-east-down (or, locally, forward-right-down) fashion. Within this scheme, the y-axis of the body frame is located at the virtual wheel axle. R maps vectors in the body frame to vectors in the spatial frame so that, for instance, the spatial angular velocity ω and the body angular velocity vector ωb are related by ω = Rωb and ωb = RT ω. Here the adjectives spatial and body refer to the frame (or basis) in which the vector is expressed. Similarly, p = x + R pb gives the spatial coordinates of a point on the body with body coordinates pb ∈ R3 . Abusing notation, we will sometimes write x = (x, y, z)T with context indicating whether x is the position vector or the scalar first component of the position vector. Within the expected operating region, the orientation R can be parametrized as R = Rz (ψ) Rx (ϕ) Ry (θ) cθ cψ − sϕ sθ sψ −cϕ sψ sθ cψ + sϕ cθ sψ = cθ sψ + sϕ sθ cψ cϕ cψ sθ sψ − sϕ cθ cψ −cϕ sθ sϕ cϕ cθ where cϕ = cos ϕ, etc. The order of the rotations expresses the fact that rotation of the wheel about
its axle (the body y-axis) occurs first, setting the wheel pitch angle θ, followed with a lean of the wheel to set the roll angle ϕ and a yaw to set the heading angle ψ. The mapping (ϕ, θ, ψ) 7→ R is one-to-one and onto in the region |ϕ| < π/2 which includes the desired operating region. Note that θ decreases with forward motion, ϕ is positive when leaned right, and ψ increases in a right hand turn. Abusing notation, we write φ = (ϕ, θ, ψ)T for this vector of Euler angles. Although this notation leads to the somewhat amusing ϕ = φ1 , the two forms of the Greek letter phi (scalar ϕ and vector φ) should be easy enough to keep straight by form and context. Note that the expressions used for principal axis rotations (Rx (·), etc.) can be obtained by setting two of the three angles to zero in the above expression for R. These rotations follow the standard right hand rule. We will use ek , k = 1, 2, . . ., to denote the standard basis vectors with 1 in the kth position and zero otherwise so that, e.g., φ1 = eT1 φ with e1 = (1, 0, 0)T . The vector q = (x, y, z, ϕ, θ, ψ)T provides a valid set of generalized coordinates for dynamics calculations on the region {q ∈ R6 : |q4 | < π/2}. In addition to the spatial frame and the body frame (with orientation R and origin located at x), we find it useful to define some additional frames. The frame with origin x and orientation Rz (ψ) is called the wheel heading frame. We call the frame with origin x and orientation Rz (ψ)Rx (ϕ) the rolled heading frame. Vectors that are represented in frames other than spatial are given an appropriate subscript. Thus, for instance, the angular velocity satisfies ω = Rz (ψ) ωh = Rz (ψ)Rx (ϕ) ωrh = Rz (ψ)Rx (ϕ)Ry (θ) ωb in the wheel heading (h), rolled heading (rh), and body (b) frames, respectively. We consider a generalized torus, geometrically and dynamically symmetric about its primary axis of revolution (the axle). For the sake of concreteness, it may be helpful to keep in mind a standard torus, representing a wheel for which the tire has a constant profile radius of ρ and the wheel plane radius is r + ρ. In this and the general case, the center of mass (with mass m) is located at the origin of the body frame and the (body frame) inertia tensor Ib is diagonal with components Ix , Iy , and Iz with Iz = Ix providing rotational symmetry. This dynamic symmetry implies that Ry (θ) Ib Ry (θ)T = Ib . Geometric (contact point) symmetric with respect to rotation about the body y-axis will be discussed below. Roughly speaking, the generalized torus needs to be a surface of revolution.
Kinetic energy
Contact point geometry
We begin by developing the formula for the kinetic energy with respect to the generalized coordinates and velocities. Along the way, we will develop some experience with some rigid body coordinate frames of interest. The kinetic energy of a rigid body system is the sum of the the translational and rotational kinetic energies
Given a horizontal road surface, the contact point between the tire (with quite general profile) and the road occurs at the bottom of the tire so that the uniquely defined contact point and the axle (the body y-axis) lie within a common vertical plane. This plane is, in fact, the x-z-plane of the wheel heading frame. In particular, for a general class of torus (or tire) profiles, there is a roll dependent mapping ch (ϕ) with eT1 ch (ϕ) ≡ 0 specifying the contact point with respect to the wheel heading frame.
T = mx˙ T x/2 ˙ + ωbT Ib ωb /2 ωb is the angular velocity vector expressed in the body frame. Thus, we need to express the rotational kinematic energy in terms of the generalized coordinates and velocities.
yh
To this end, let ω ∈ R3 denote the angular velocity of the wheel expressed in the spatial frame so that R˙ = ω × R = ω bR and ω b = R˙ RT
where
Note that Jω (φ) can be written as Jω (φ) = Rz (ψ) Jωh (ϕ) where 1 0 0 Jωh (ϕ) = 0 cϕ 0 0 sϕ 1 maps φ˙ to ωh , the angular velocity expressed in the wheel heading frame. Similarly, using ωb = RT ω, the Jacobian mapping φ˙ to ωb is given by Jωb (φ) = Ry (θ)T Rx (ϕ)T Rz (ψ)T Rz (ψ)Jωh (ϕ) = Ry (θ)T Rx (ϕ)T Jωh (ϕ) . Exploiting the rotational symmetry of Ib , we see that the rotational kinetic energy is given by ˙ ωbT Ib ωb /2 = φ˙ T I(ϕ)φ/2 where the generalized coordinates inertia tensor I(ϕ) is given by I(ϕ) = Jωh (ϕ)T Rx (ϕ) Ib Rx (ϕ)T Jωh (ϕ) or, equivalently, I(ϕ) = Jωrh (ϕ)T Ib Jωrh (ϕ) 10 0 Jωrh (ϕ) = 0 1 sϕ 0 0 cϕ
zh
ch ()
0 −ω3 ω2 ω b = ω3 0 −ω1 −ω2 ω1 0 is the (skew-symmetric) matrix representation of the linear mapping v 7→ ω × v = ω b v. Direct calculations show that the Jacobian Jω (·) mapping φ˙ to ω satisfies ϕ˙ cψ −sψ cϕ 0 ω = Jω (φ) φ˙ = sψ cψ cϕ 0 θ˙ . 0 sϕ 1 ψ˙
where
Fig. 1. The contact point vector ch (ϕ) in the wheel heading frame. Now, the generalized torus is a surface of revolution. How general is the curve specifying the contacting surface? Clearly the surface should be strictly convex to ensure that the contact point is unique and the surface should be sufficiently smooth to ensure that the dynamics are well defined and perhaps even differentiable. To this end, we will suppose that the contact surface curve is such that the instantaneous curvature is well defined at each contact point giving a radius of curvature function ρ(ϕ), |ϕ| ≤ π/2. We will suppose that ρ(ϕ) is continuously differentiable and bounded according to 0 < ρ0 ≤ ρ(ϕ) ≤ ρ1 < ∞ on its domain. Integrating the instantaneous radius function, the contact profile is given, in a local frame, by ] [ [ ] ∫ ϕ cos ζ w(ϕ) ρ(ζ) dζ = sin ζ h(ϕ) 0 so that
w0 (ϕ) = ρ(ϕ) cos ϕ h0 (ϕ) = ρ(ϕ) sin ϕ and the slope of the tangent line at the contact point (w(ϕ), h(ϕ)) is h0 (ϕ)/w0 (ϕ) = tan ϕ as expected. h
maps φ˙ to ωrh , the angular velocity in the rolled heading frame. Explicitly, Ix 0 0 . Iy sϕ I(ϕ) = 0 Iy 2 2 0 Iy sϕ Ix cϕ + Iy sϕ
w
Fig. 2. The torus tire profile, (w(ϕ), h(ϕ)). Now, letting (w0 , h0 ) denote the (fixed) location of the center of mass, the contact point mapping ch (ϕ) is given by
0 ch (ϕ) = Rx (ϕ) w(ϕ) − w0 −(h(ϕ) − h0 ) where the minus sign on the h portion reflects the required flip when going from local w-h coordinates with h up to local y-z coordinates with z down. Differentiating ch (ϕ), we obtain c0h (ϕ) = eb1 ch (ϕ) + ρ(ϕ) e2 , c00h (ϕ) = −ch (ϕ) + ρ0 (ϕ) e2 + ρ(ϕ) e3 .
(1) (2)
For the sake of brevity, we will refer to the components of ch (ϕ) according to ch (ϕ) = (0, c2 (ϕ), c3 (ϕ))T rather than, e.g., ch2 (ϕ), and, at times, suppress the ϕ dependence, writing c2 , c3 , etc. With this notation, we see that hch (ϕ), c0h (ϕ)i = ch (ϕ)T (b e1 ch (ϕ) + ρ(ϕ)e2 ) = ρ(ϕ) c2 (ϕ) = ρ c2 and, from (1), c03 = c2 and c02 = −c3 + ρ. Define, for future use, c4 (ϕ) = cϕ c3 (ϕ)−sϕ c2 (ϕ) and note that c4 (ϕ) = eT3 crh (ϕ) is the z component of the contact point vector in the rolled heading frame. Direct computation shows that c04 = −sϕ ρ. Contact profile information enters the dynamic description of the rolling torus exclusively through the mapping ch (·). In general, the generalized torus need not be geometrically symmetric with respect to the body x-z-plane. Even in the case that the torus has a standard symmetry shape, the body x-z-plane (which passes through the center of mass) need not be a geometric symmetry plane. In any case, there should a roll angle within the valid region for which the moment due to gravity is zero, that is, there is a ϕ0 such that c2 (ϕ0 ) = 0. In this paper, we consider the general symmetric case where the curvature mapping ρ(·) is a well defined even function of ϕ with domain that includes [−π/2, π/2] and the mass is geometrically centered, w0 = 0. Furthermore, we require that the contact point always occur on the same side of the axle (the body y-axis) so that c4 (ϕ) > 0 for all ϕ ∈ (−π/2, π/2). For this class of symmetric profiles, c2 (·) is an odd function of ϕ, and c3 (·) and c4 (·) are even functions of ϕ. With the mapping ch (·) in hand, we are ready to determine the nature of the rolling constraint. The velocity at the contact point is vc = x˙ + ω × Rz (ψ) ch (ϕ) = x˙ + Rz (ψ) (ωh × ch (ϕ)) = x˙ − Rz (ψ) b ch (ϕ) ωh = x˙ − Rz (ψ) b ch (ϕ) Jωh (ϕ) φ˙ so that vc = Jvc (q) q˙ where [ ] ch (ϕ) Jωh (ϕ) Jvc (q) = I −Rz (ψ) b is the contact point velocity Jacobian, which we write as [ ] Jvc (q) = I Rz (ψ) Jcωh (ϕ) .
where Jcωh (ϕ) = −b ch (ϕ) Jωh (ϕ). The subscript cωh does not follow the subscript convention above, but rather reminds us that Jcωh (ϕ) φ˙ is a velocity vector in the wheel heading frame orientation that is involved in the contact point velocity calculation. The constraint for rolling without slipping requires that the velocity at the contact point on the wheel be the same as the velocity at the contact point on the road, vc = 0, leading to the (mostly) nonholonomic constraint Jvc (q) q˙ = 0 .
(3)
We say mostly since the portion normal to the road surface is actually the holonomic constraint z = eT3 ch (ϕ). We will also make use of the differential constraint obtained by differentiating (3) with respect to time to obtain Jvc (q) q¨ + (DJvc (q) · q) ˙ q˙ = 0 , which we write as ¨ ˙ ˙ ˜ x ¨ +Rz (ψ) Jcωh (ϕ)φ+R z (ψ)Jcωh (ϕ, φ)φ = 0 , (4) where ˙ = ψ˙ eb3 Jcω (ϕ) + ϕ˙ J 0 (ϕ) . J˜cωh (ϕ, φ) cωh h Lagrangian dynamics We now determine the dynamics of the constrained system. The Lagrangian is L = T − V ˙ where T = m x˙ T x/2 ˙ + φ˙ T I(ϕ)φ/2 is the kinetic energy of the rolling wheel and and V = −mg eT3 x is the potential energy where the ‘−’ sign is due to the fact that z = eT3 x is positive down. Using L and Jvc , the equations of motion satisfy d ∂L T ∂L T − = (0, τ )T − JvTc f dt ∂ q˙ ∂q
(5)
where −f is the contact point force (in spatial frame) needed to enforce the rolling constraint (3) and τ ∈ R3 is a (to be specified) torque applied at the center of mass. The sign on the force has been chosen so that the tire pushes on the road with force f . Equations (5) and (4) together result in the partitioned system x ¨ mI 0 I 0 T T ¨ I(ϕ) J (ϕ) R (ψ) φ cωh z I Rz (ψ)Jcωh (ϕ) +
0 −mg e3 ˙ Nφ (ϕ, φ)
˙ φ˙ Rz (ψ)J˜cωh (ϕ, φ)
f 0 τ = ,
0
˙ This ˙ = (ϕI where Nφ (ϕ, φ) ˙ − e1 φ˙ T /2) I0 (ϕ)φ. system of equations can be solved for the force f ¨ Moreover, the structure and the accelerations x ¨, φ. of these equations allows us to explicitly eliminate the constraints and translation variables from
dynamic consideration. To this end, note that the center of mass acceleration satisfies ˙ φ] ˙ x ¨ = −Rz (ψ)[Jcω (ϕ)φ¨ + J˜cω (ϕ, φ) h
h
so that the contact force is given by ˙ φ] ˙ . f = mg e3 + mRz (ψ)[Jcωh (ϕ)φ¨ + J˜cωh (ϕ, φ) Applying the contact force to the rotational dynamics, we obtain [ ] ˙ I(ϕ) + mJcωh (ϕ)T Jcωh (ϕ) φ¨ + Nφ (ϕ, φ) T ˜ ˙ φ˙ + mg ch2 (ϕ) e1 = τ + mJcω (ϕ) Jcω (ϕ, φ) h
h
where ch2 (ϕ) = eT2 ch (ϕ) that Jcωh (ϕ)T Rz (ψ)T e3
and we have used the fact = eT2 ch (ϕ) e1 . For future reference, note that the contact force expressed in the wheel heading frame orientation is given by ˙ φ] ˙ . fh = mg e3 + m [Jcω (ϕ)φ¨ + J˜cω (ϕ, φ) h
h
To this point, we’ve considered the the externally applied torque τ ∈ R3 to be general. In fact, the driven torus is acted upon by two exogenous inputs: a drive torque τdr that acts around the wheel axle and a steer torque τst that acts about a ficticious steering axis set at a castor angle θst . In the spatial frame orientation, the drive torque is applied about the axis ξdr = Rz (ψ)Rx (ϕ) e2 . T ˙ we find Equating power hξdr τdr , ωi = hJdr τdr , φi, that
Jdr (ϕ)T = Jωrh (ϕ)T e2 = (0, 1, sϕ )T .
so that T
Jst (ϕ) = Jωrh (ϕ) Ry (θst ) e3 = (sθst , 0, cθst cϕ ) . The applied torque is thus ] [ 0 sθst τ ¯ 0 dr = B(ϕ) u. τ = 1 τst sϕ cθst cϕ ˙ ψ) ˙ T, Thus, the reduced dynamics, with state (ϕ, ϕ, ˙ θ, is given by ˙ + G(ϕ) ¯I(ϕ)φ¨ + N ¯ (ϕ, φ) ¯ ¯ = B(ϕ) u. (6) where
˙ T = m x˙ T x/2 ˙ + φ˙ T I(ϕ)φ/2 ˙ = φ˙ T (mJcω (ϕ)T Jcω (ϕ) + I(ϕ))φ/2 h
h
˙ = φ˙ T ¯I(ϕ)φ/2 so that the same generalized inertia tensor ¯I(ϕ) appears both in the dynamics and the kinetic energy as one might expect. Now, the potential energy is given by V = −mg eT3 x = mg eT3 ch (ϕ) where we use the fact that the z component of the contact point is constant at z = 0. Setting E = T + V and u = 0, we compute E˙ = φ˙ T [¯I0 (ϕ)ϕ] ˙ φ˙ /2 + φ˙ T ¯I(ϕ)φ¨ + mg ϕ˙ eT c0 (ϕ) 3 h
0 = φ˙ T [I0 (ϕ)ϕ˙ + 2mJcωh (ϕ)T Jcω (ϕ)ϕ] ˙ φ˙ /2 h
˙ φ˙ ˙ − mφ˙ T Jcω (ϕ)T J˜cω (ϕ, φ) − φ˙ T Nφ (ϕ, φ) h h + mg ϕ˙ (eT3 c0h (ϕ) − eT2 ch (ϕ)) 0 = φ˙ T [I0 (ϕ)ϕ˙ + 2mJcωh (ϕ)T Jcω (ϕ)ϕ] ˙ φ˙ /2 h
0 − mφ˙ T Jcωh (ϕ)T [ψ˙ eb3 Jcωh (ϕ) + ϕ˙ Jcω (ϕ)] φ˙ h
˙ T eb3 (Jcω (ϕ) φ) ˙ = 0 = −mψ˙ (Jcωh (ϕ) φ) h
ξst = Rz (ψ)Rx (ϕ)Ry (θst ) e3 T
Since the contact point force f that ensures that the torus rolls without slipping does no work, we expect the unforced system (u = 0) to conserve energy. Although verification of energy conservation is, in principle, a straighforward calculation, there are some interesting features. First, using x˙ from (7), we find that
− φ˙ T [ϕI ˙ − e1 φ˙ T /2] I0 (ϕ)φ˙
The steering torque is applied, in spatial reference frame, about the virtual steering axis
T
Conservation of energy
[ ] ¯I(ϕ) = I(ϕ) + mJcω (ϕ)T Jcω (ϕ) , h h
˙ = Nφ (ϕ, φ) ˙ + mJcω (ϕ)T J˜cω (ϕ, φ) ˙ φ˙ , ¯ (ϕ, φ) N h h and ¯ G(ϕ) = mg c2 (ϕ) e1 . The translational motion is obtained by integrating the rolling constraint x˙ = −Rz (ψ)Jcω (ϕ) φ˙ (7) h
from which we see that, kinematically speaking, it is necessary to keep track of the heading angle ψ as well. It is not necessary to keep track of θ unless one desires, for example, to trace out the contact point path on the surface of the torus.
where the last equality follows from the fact that eb3 is skew symmetric. We have also made use of the fact that (1) implies that eT3 c0h (ϕ) = eT2 ch (ϕ) to eliminate the terms involving gravity. When the externally applied torque is not zero, the power supplied to the system is given by ¯ E˙ = φ˙ T τ = φ˙ T B(ϕ) u. Structure and detail of the dynamic equations The reduced dynamics (6) has a structure that is interesting and exploitable. Since Ix +mkch k2 0 0 ¯I(ϕ) = 0 Iy +mc24 Iy sϕ −mc2 c4 0 ∗ Ix c2ϕ +Iy s2ϕ +mc22 we see that its inversion does not mix up the apparent torques too much: the expressions for θ¨ and ψ¨ do not depend on expressions from the first ˙ are ¯φ (ϕ, φ) row (ϕ) ¨ of (6). The components of N given by ( ) ˙ ¯φ (ϕ, φ) N = m ρ c2 ϕ˙ 2 − (Iy cϕ + m c3 c4 ) θ˙ψ˙ 1 + ((Ix − Iy )cϕ sϕ + m c2 c3 ) ψ˙ 2
( ) ˙ ¯φ (ϕ, φ) N
[ = ϕ˙ C(ϕ) 23
θ˙ ψ˙
]
where C(ϕ) is the 2 × 2 matrix ] [ −mρsϕ c4 Iy cϕ +m(2c3 −ρ)c4 . Iy cϕ +mρsϕ c2 2(Iy −Ix )cϕ sϕ −m(2c3 −ρ)c2
θ˙2 ≥ 4mgc2 ((Ix −Iy )cϕ sϕ +mc2 c3 )/(Iy cϕ +mc3 c4 )2 . When the inequality is strict, there will be two ˙ for each roll angle. distinct turn rates, ψ,
Integrability of the undriven rolling torus It turns out that (almost) every solution of (6) with u(t) ≡ 0 that satisfies |ϕ(t)| < π/2, t ≥ 0, is, in fact, periodic. The key is to note that (6) with u = 0 has the form ϕ¨ = f (ϕ, ϕ, ˙ w) w˙ = ϕ˙ A(ϕ) w
every −θ˙ > 0. In contrast, when turning a certain level of speed is needed to keep the wheel at the chosen lean angle. Indeed, for nonzero roll ϕ, the requirement that ψ˙ in (10) be real requires that
(8)
where w ∈ Rp and A(ϕ) is a bounded p×p matrix, ˙ ψ) ˙ T in this case. with p = 2 and w = (θ,
Motorcycle wheel example Consider now a standard torus with constant profile radius ρ and center of mass located at (w0 , h0 ) = (0, r + ρ) so that the overall wheel radius is r + ρ and ch (ϕ) = (0, −rsϕ , ρ + rcϕ )T . Numerically, let r = 0.2, ρ = 0.1, Iy = 0.64, Ix = 0.38, m = 14.7, and g = 9.81. The values are 10
Let W (·) satisfy W 0 (ϕ) = A(ϕ) W (ϕ),
W (0) = I,
5
on the ϕ region of interest and define w(ϕ; ¯ ϕ0 , w0 ) = W (ϕ) W (ϕ0 )−1 w0 . It follows easily that every trajectory of (6) is determined by a trajectory of the (parametrized) two dimensional system
0
ϕ¨ = f (ϕ, ϕ, ˙ w(ϕ; ¯ ϕ0 , w0 )) .
−5
(9)
Furthermore, since (9) is conservative, the energy provides a first integral making stability analysis trivial and periodic solutions transparent. In particular, an equilirium point of (9) is unstable if the linearization has real eigenvalues and Lyapunov stable otherwise. Furthermore, these eigenvalues are the nonzeros eigenvalues of the linearization of (6). Equilibria We now consider the set of constant operating conditions (xe , ue ) for the reduced dynamics (6). Each constant state has the form xe = ˙ ψ) ˙ T . Since the second and third compo(ϕ, 0, θ, ˙ are zero by ϕ˙ = 0 and rows ¯ (ϕ, φ) nents of N ¯ 2 and 3 of B(ϕ) are linearly independent on (−π/2, π/2), we find that ue = 0 for every constant operating condition of the dynamic control ˙ and ψ˙ system (6). Thus, for equilibrium, ϕ, θ, must satisfy ((Ix − Iy )cϕ sϕ + m c2 c3 ) ψ˙ 2 − (Iy cϕ + m c3 c4 ) θ˙ψ˙ + mg c2 = 0
(10)
where we recall that c2 , c3 , and c4 are roll dependent components of the contact point vector: ch (ϕ) = (0, c2 , c2 )T and c4 = eT3 crh (ϕ). We are interested in the equilibrium motions with nonzero forward velocity, −θ˙ > 0. For the symmetric wheel case considered, straight running motion occurs with ϕ = 0 and ψ˙ = 0 for
0
50
100
150
200
250
300
350
Fig. 3. Equilibrium stability: linearized natural ˙ for roll angles ϕ from zero frequency versus |θ| to 40 degrees in increments of 2 degrees. representative of those for a rear motorcycle tire and wheel (Sharp et al., 2004). Figure 3 presents quantitative equilibrium stability results for the rolling wheel with these parameters. The natural frequency of an imaginary pair is plotted positive while the frequency of a hyperbolic pair is plotted negative, implying unstable. Blue indicates the ˙ condition, and more rapidly turning (higher ψ) green indicates the usual turning rate. Thus, at ϕ = 18deg, we see that the green branch is stable ˙ ≈ 107deg/sec to 118deg/sec, unstable from |θ| from 118 to 152deg/sec, and stable for higher speeds. Every blue branch is always stable. References Arnold, V. I., V. V. Kozlov and A. I. Neishtadt (2006). Mathematical Aspects of Classical and Celestial Mechanics. Springer-Verlag. Bloch, A. M. (2003). Nonholonomic Mechanics and Control. Springer-Verlag. Bullo, F. and A. D. Lewis (2004). Geometric Control of Mechanical Systems. Springer-Verlag. O’Reilly, O.M. (1996). The dynamics of rolling disks and sliding disks. Nonlinear Dynamics 10, 287–305. Sharp, R.S., S. Evangelou and D.J.N. Limebeer (2004). Advances in the modelling of motorcycle dynamics. Multibody System Dynamics 12(3), 251–283.
Abstract: We develop the equations of motion for a torus with a general contact profile that rolls without slipping and that is actuated using drive and steer torques. The form of the dynamics plus nonholonomic constaints allows an easy reduction to a dynamic control system of dimension 4 including the three angular rates and the roll (or lean) angle. Using a standard torus with geometric and inertial parameters similar to those of a motorcycle tire and wheel, we find that the rolling torus exhibits, with increasing forward speed, an unexpected transition c from stability to instability and back to stability! Copyright °2007 IFAC Keywords: Nonlinear control systems, Lagrangian dynamics, reduction, integrability, motorcycle tires, instability.
We develop the equations of motion for a torus with a general contact profile that rolls without slipping and that is actuated using drive and steer torques that are somehow magically applied. The purpose of this study is to work toward a better understand of the dynamics and control of motorcycles operating in high traction/low slip conditions. Proceeding carefully to encourage the discovery of key kinematic and dynamic features, we find that, for a general class of contact surface (or tire) profiles, the form of the dynamics plus nonholonomic constaints allows an easy reduction to a dynamic control system of dimension 4 including the three angular rates and the roll (or lean) angle. We show, by elementary means, that the undriven system admits to classical integrability so that the pitch and yaw rates can be given as (instantaneous) functions of the roll angle. Using this, together with conservation of energy, we find that the stability of steady motions
is easy to study and, furthermore, that almost every (undriven) motion will be periodic. For each roll angle, steady turning is only possible when the forward velocity (or wheel pitch rate) is at or above a critical value determined by a simple equation, quadratic in the yaw rate, required for equilibrium. Above that value, there are two equilibria. Using a standard torus with geometric and inertial parameters similar to those of a motorcycle tire and wheel, we demonstrate that dynamics can be strange and intuition lacking. Indeed, for the parameters chosen and a range of steady roll angles, the rolling torus exhibits, with increasing forward speed, a transition from stability to instability and back to stability! The study of rigid bodies subject to constraints, especially nonholonomic constraints, has a long and exciting history and continues to this day, see, e.g., (Bloch, 2003), (Bullo and Lewis, 2004), and (Arnold et al., 2006), which include many
classical and modern references. Indeed, (classical) integrability and the ubiquity of periodic solutions for rolling disks and bodies of revolution were established in the late 1800’s and early 1900’s, with special emphasis on expressing these solutions using, e.g., hypergeometric or Legendre functions. Though the dynamics of a rolling disk are quite similar, the bifurcation analysis in (O’Reilly, 1996) does not indicate the presence of the stability-instability-stability transition that we observe in our numerical example. Generalized coordinates for the rolling torus For the purpose of deriving the equations of motion, we consider the system as a rigid body subject to the constraint of rolling without slipping. We will see that, since the region of operation is restricted lean angles less than 90 degrees, there is a convenient set of generalized coordinates q ∈ R6 parametrizing the configuration of the rolling torus, enabling the (global) determination of the rolling torus dynamics using the Lagrangian method. In this approach, the rolling contact constraint is enforced by contact point forces acting as Lagrange multipliers is the system dynamic equations. Our goal, then, is to establish this set of generalized coordinates, express the kinetic and potential energies, T and V , in terms of the generalized coordinates and velocities, and determine the (mostly nonholonomic) constraints with respect to same. Let x ∈ R3 and R ∈ SO(3) denote the position and orientation of a frame fixed at the wheel center of mass with respect to a fixed inertial (or spatial ) frame with x-y-z axes oriented in a north-east-down (or, locally, forward-right-down) fashion. Within this scheme, the y-axis of the body frame is located at the virtual wheel axle. R maps vectors in the body frame to vectors in the spatial frame so that, for instance, the spatial angular velocity ω and the body angular velocity vector ωb are related by ω = Rωb and ωb = RT ω. Here the adjectives spatial and body refer to the frame (or basis) in which the vector is expressed. Similarly, p = x + R pb gives the spatial coordinates of a point on the body with body coordinates pb ∈ R3 . Abusing notation, we will sometimes write x = (x, y, z)T with context indicating whether x is the position vector or the scalar first component of the position vector. Within the expected operating region, the orientation R can be parametrized as R = Rz (ψ) Rx (ϕ) Ry (θ) cθ cψ − sϕ sθ sψ −cϕ sψ sθ cψ + sϕ cθ sψ = cθ sψ + sϕ sθ cψ cϕ cψ sθ sψ − sϕ cθ cψ −cϕ sθ sϕ cϕ cθ where cϕ = cos ϕ, etc. The order of the rotations expresses the fact that rotation of the wheel about
its axle (the body y-axis) occurs first, setting the wheel pitch angle θ, followed with a lean of the wheel to set the roll angle ϕ and a yaw to set the heading angle ψ. The mapping (ϕ, θ, ψ) 7→ R is one-to-one and onto in the region |ϕ| < π/2 which includes the desired operating region. Note that θ decreases with forward motion, ϕ is positive when leaned right, and ψ increases in a right hand turn. Abusing notation, we write φ = (ϕ, θ, ψ)T for this vector of Euler angles. Although this notation leads to the somewhat amusing ϕ = φ1 , the two forms of the Greek letter phi (scalar ϕ and vector φ) should be easy enough to keep straight by form and context. Note that the expressions used for principal axis rotations (Rx (·), etc.) can be obtained by setting two of the three angles to zero in the above expression for R. These rotations follow the standard right hand rule. We will use ek , k = 1, 2, . . ., to denote the standard basis vectors with 1 in the kth position and zero otherwise so that, e.g., φ1 = eT1 φ with e1 = (1, 0, 0)T . The vector q = (x, y, z, ϕ, θ, ψ)T provides a valid set of generalized coordinates for dynamics calculations on the region {q ∈ R6 : |q4 | < π/2}. In addition to the spatial frame and the body frame (with orientation R and origin located at x), we find it useful to define some additional frames. The frame with origin x and orientation Rz (ψ) is called the wheel heading frame. We call the frame with origin x and orientation Rz (ψ)Rx (ϕ) the rolled heading frame. Vectors that are represented in frames other than spatial are given an appropriate subscript. Thus, for instance, the angular velocity satisfies ω = Rz (ψ) ωh = Rz (ψ)Rx (ϕ) ωrh = Rz (ψ)Rx (ϕ)Ry (θ) ωb in the wheel heading (h), rolled heading (rh), and body (b) frames, respectively. We consider a generalized torus, geometrically and dynamically symmetric about its primary axis of revolution (the axle). For the sake of concreteness, it may be helpful to keep in mind a standard torus, representing a wheel for which the tire has a constant profile radius of ρ and the wheel plane radius is r + ρ. In this and the general case, the center of mass (with mass m) is located at the origin of the body frame and the (body frame) inertia tensor Ib is diagonal with components Ix , Iy , and Iz with Iz = Ix providing rotational symmetry. This dynamic symmetry implies that Ry (θ) Ib Ry (θ)T = Ib . Geometric (contact point) symmetric with respect to rotation about the body y-axis will be discussed below. Roughly speaking, the generalized torus needs to be a surface of revolution.
Kinetic energy
Contact point geometry
We begin by developing the formula for the kinetic energy with respect to the generalized coordinates and velocities. Along the way, we will develop some experience with some rigid body coordinate frames of interest. The kinetic energy of a rigid body system is the sum of the the translational and rotational kinetic energies
Given a horizontal road surface, the contact point between the tire (with quite general profile) and the road occurs at the bottom of the tire so that the uniquely defined contact point and the axle (the body y-axis) lie within a common vertical plane. This plane is, in fact, the x-z-plane of the wheel heading frame. In particular, for a general class of torus (or tire) profiles, there is a roll dependent mapping ch (ϕ) with eT1 ch (ϕ) ≡ 0 specifying the contact point with respect to the wheel heading frame.
T = mx˙ T x/2 ˙ + ωbT Ib ωb /2 ωb is the angular velocity vector expressed in the body frame. Thus, we need to express the rotational kinematic energy in terms of the generalized coordinates and velocities.
yh
To this end, let ω ∈ R3 denote the angular velocity of the wheel expressed in the spatial frame so that R˙ = ω × R = ω bR and ω b = R˙ RT
where
Note that Jω (φ) can be written as Jω (φ) = Rz (ψ) Jωh (ϕ) where 1 0 0 Jωh (ϕ) = 0 cϕ 0 0 sϕ 1 maps φ˙ to ωh , the angular velocity expressed in the wheel heading frame. Similarly, using ωb = RT ω, the Jacobian mapping φ˙ to ωb is given by Jωb (φ) = Ry (θ)T Rx (ϕ)T Rz (ψ)T Rz (ψ)Jωh (ϕ) = Ry (θ)T Rx (ϕ)T Jωh (ϕ) . Exploiting the rotational symmetry of Ib , we see that the rotational kinetic energy is given by ˙ ωbT Ib ωb /2 = φ˙ T I(ϕ)φ/2 where the generalized coordinates inertia tensor I(ϕ) is given by I(ϕ) = Jωh (ϕ)T Rx (ϕ) Ib Rx (ϕ)T Jωh (ϕ) or, equivalently, I(ϕ) = Jωrh (ϕ)T Ib Jωrh (ϕ) 10 0 Jωrh (ϕ) = 0 1 sϕ 0 0 cϕ
zh
ch ()
0 −ω3 ω2 ω b = ω3 0 −ω1 −ω2 ω1 0 is the (skew-symmetric) matrix representation of the linear mapping v 7→ ω × v = ω b v. Direct calculations show that the Jacobian Jω (·) mapping φ˙ to ω satisfies ϕ˙ cψ −sψ cϕ 0 ω = Jω (φ) φ˙ = sψ cψ cϕ 0 θ˙ . 0 sϕ 1 ψ˙
where
Fig. 1. The contact point vector ch (ϕ) in the wheel heading frame. Now, the generalized torus is a surface of revolution. How general is the curve specifying the contacting surface? Clearly the surface should be strictly convex to ensure that the contact point is unique and the surface should be sufficiently smooth to ensure that the dynamics are well defined and perhaps even differentiable. To this end, we will suppose that the contact surface curve is such that the instantaneous curvature is well defined at each contact point giving a radius of curvature function ρ(ϕ), |ϕ| ≤ π/2. We will suppose that ρ(ϕ) is continuously differentiable and bounded according to 0 < ρ0 ≤ ρ(ϕ) ≤ ρ1 < ∞ on its domain. Integrating the instantaneous radius function, the contact profile is given, in a local frame, by ] [ [ ] ∫ ϕ cos ζ w(ϕ) ρ(ζ) dζ = sin ζ h(ϕ) 0 so that
w0 (ϕ) = ρ(ϕ) cos ϕ h0 (ϕ) = ρ(ϕ) sin ϕ and the slope of the tangent line at the contact point (w(ϕ), h(ϕ)) is h0 (ϕ)/w0 (ϕ) = tan ϕ as expected. h
maps φ˙ to ωrh , the angular velocity in the rolled heading frame. Explicitly, Ix 0 0 . Iy sϕ I(ϕ) = 0 Iy 2 2 0 Iy sϕ Ix cϕ + Iy sϕ
w
Fig. 2. The torus tire profile, (w(ϕ), h(ϕ)). Now, letting (w0 , h0 ) denote the (fixed) location of the center of mass, the contact point mapping ch (ϕ) is given by
0 ch (ϕ) = Rx (ϕ) w(ϕ) − w0 −(h(ϕ) − h0 ) where the minus sign on the h portion reflects the required flip when going from local w-h coordinates with h up to local y-z coordinates with z down. Differentiating ch (ϕ), we obtain c0h (ϕ) = eb1 ch (ϕ) + ρ(ϕ) e2 , c00h (ϕ) = −ch (ϕ) + ρ0 (ϕ) e2 + ρ(ϕ) e3 .
(1) (2)
For the sake of brevity, we will refer to the components of ch (ϕ) according to ch (ϕ) = (0, c2 (ϕ), c3 (ϕ))T rather than, e.g., ch2 (ϕ), and, at times, suppress the ϕ dependence, writing c2 , c3 , etc. With this notation, we see that hch (ϕ), c0h (ϕ)i = ch (ϕ)T (b e1 ch (ϕ) + ρ(ϕ)e2 ) = ρ(ϕ) c2 (ϕ) = ρ c2 and, from (1), c03 = c2 and c02 = −c3 + ρ. Define, for future use, c4 (ϕ) = cϕ c3 (ϕ)−sϕ c2 (ϕ) and note that c4 (ϕ) = eT3 crh (ϕ) is the z component of the contact point vector in the rolled heading frame. Direct computation shows that c04 = −sϕ ρ. Contact profile information enters the dynamic description of the rolling torus exclusively through the mapping ch (·). In general, the generalized torus need not be geometrically symmetric with respect to the body x-z-plane. Even in the case that the torus has a standard symmetry shape, the body x-z-plane (which passes through the center of mass) need not be a geometric symmetry plane. In any case, there should a roll angle within the valid region for which the moment due to gravity is zero, that is, there is a ϕ0 such that c2 (ϕ0 ) = 0. In this paper, we consider the general symmetric case where the curvature mapping ρ(·) is a well defined even function of ϕ with domain that includes [−π/2, π/2] and the mass is geometrically centered, w0 = 0. Furthermore, we require that the contact point always occur on the same side of the axle (the body y-axis) so that c4 (ϕ) > 0 for all ϕ ∈ (−π/2, π/2). For this class of symmetric profiles, c2 (·) is an odd function of ϕ, and c3 (·) and c4 (·) are even functions of ϕ. With the mapping ch (·) in hand, we are ready to determine the nature of the rolling constraint. The velocity at the contact point is vc = x˙ + ω × Rz (ψ) ch (ϕ) = x˙ + Rz (ψ) (ωh × ch (ϕ)) = x˙ − Rz (ψ) b ch (ϕ) ωh = x˙ − Rz (ψ) b ch (ϕ) Jωh (ϕ) φ˙ so that vc = Jvc (q) q˙ where [ ] ch (ϕ) Jωh (ϕ) Jvc (q) = I −Rz (ψ) b is the contact point velocity Jacobian, which we write as [ ] Jvc (q) = I Rz (ψ) Jcωh (ϕ) .
where Jcωh (ϕ) = −b ch (ϕ) Jωh (ϕ). The subscript cωh does not follow the subscript convention above, but rather reminds us that Jcωh (ϕ) φ˙ is a velocity vector in the wheel heading frame orientation that is involved in the contact point velocity calculation. The constraint for rolling without slipping requires that the velocity at the contact point on the wheel be the same as the velocity at the contact point on the road, vc = 0, leading to the (mostly) nonholonomic constraint Jvc (q) q˙ = 0 .
(3)
We say mostly since the portion normal to the road surface is actually the holonomic constraint z = eT3 ch (ϕ). We will also make use of the differential constraint obtained by differentiating (3) with respect to time to obtain Jvc (q) q¨ + (DJvc (q) · q) ˙ q˙ = 0 , which we write as ¨ ˙ ˙ ˜ x ¨ +Rz (ψ) Jcωh (ϕ)φ+R z (ψ)Jcωh (ϕ, φ)φ = 0 , (4) where ˙ = ψ˙ eb3 Jcω (ϕ) + ϕ˙ J 0 (ϕ) . J˜cωh (ϕ, φ) cωh h Lagrangian dynamics We now determine the dynamics of the constrained system. The Lagrangian is L = T − V ˙ where T = m x˙ T x/2 ˙ + φ˙ T I(ϕ)φ/2 is the kinetic energy of the rolling wheel and and V = −mg eT3 x is the potential energy where the ‘−’ sign is due to the fact that z = eT3 x is positive down. Using L and Jvc , the equations of motion satisfy d ∂L T ∂L T − = (0, τ )T − JvTc f dt ∂ q˙ ∂q
(5)
where −f is the contact point force (in spatial frame) needed to enforce the rolling constraint (3) and τ ∈ R3 is a (to be specified) torque applied at the center of mass. The sign on the force has been chosen so that the tire pushes on the road with force f . Equations (5) and (4) together result in the partitioned system x ¨ mI 0 I 0 T T ¨ I(ϕ) J (ϕ) R (ψ) φ cωh z I Rz (ψ)Jcωh (ϕ) +
0 −mg e3 ˙ Nφ (ϕ, φ)
˙ φ˙ Rz (ψ)J˜cωh (ϕ, φ)
f 0 τ = ,
0
˙ This ˙ = (ϕI where Nφ (ϕ, φ) ˙ − e1 φ˙ T /2) I0 (ϕ)φ. system of equations can be solved for the force f ¨ Moreover, the structure and the accelerations x ¨, φ. of these equations allows us to explicitly eliminate the constraints and translation variables from
dynamic consideration. To this end, note that the center of mass acceleration satisfies ˙ φ] ˙ x ¨ = −Rz (ψ)[Jcω (ϕ)φ¨ + J˜cω (ϕ, φ) h
h
so that the contact force is given by ˙ φ] ˙ . f = mg e3 + mRz (ψ)[Jcωh (ϕ)φ¨ + J˜cωh (ϕ, φ) Applying the contact force to the rotational dynamics, we obtain [ ] ˙ I(ϕ) + mJcωh (ϕ)T Jcωh (ϕ) φ¨ + Nφ (ϕ, φ) T ˜ ˙ φ˙ + mg ch2 (ϕ) e1 = τ + mJcω (ϕ) Jcω (ϕ, φ) h
h
where ch2 (ϕ) = eT2 ch (ϕ) that Jcωh (ϕ)T Rz (ψ)T e3
and we have used the fact = eT2 ch (ϕ) e1 . For future reference, note that the contact force expressed in the wheel heading frame orientation is given by ˙ φ] ˙ . fh = mg e3 + m [Jcω (ϕ)φ¨ + J˜cω (ϕ, φ) h
h
To this point, we’ve considered the the externally applied torque τ ∈ R3 to be general. In fact, the driven torus is acted upon by two exogenous inputs: a drive torque τdr that acts around the wheel axle and a steer torque τst that acts about a ficticious steering axis set at a castor angle θst . In the spatial frame orientation, the drive torque is applied about the axis ξdr = Rz (ψ)Rx (ϕ) e2 . T ˙ we find Equating power hξdr τdr , ωi = hJdr τdr , φi, that
Jdr (ϕ)T = Jωrh (ϕ)T e2 = (0, 1, sϕ )T .
so that T
Jst (ϕ) = Jωrh (ϕ) Ry (θst ) e3 = (sθst , 0, cθst cϕ ) . The applied torque is thus ] [ 0 sθst τ ¯ 0 dr = B(ϕ) u. τ = 1 τst sϕ cθst cϕ ˙ ψ) ˙ T, Thus, the reduced dynamics, with state (ϕ, ϕ, ˙ θ, is given by ˙ + G(ϕ) ¯I(ϕ)φ¨ + N ¯ (ϕ, φ) ¯ ¯ = B(ϕ) u. (6) where
˙ T = m x˙ T x/2 ˙ + φ˙ T I(ϕ)φ/2 ˙ = φ˙ T (mJcω (ϕ)T Jcω (ϕ) + I(ϕ))φ/2 h
h
˙ = φ˙ T ¯I(ϕ)φ/2 so that the same generalized inertia tensor ¯I(ϕ) appears both in the dynamics and the kinetic energy as one might expect. Now, the potential energy is given by V = −mg eT3 x = mg eT3 ch (ϕ) where we use the fact that the z component of the contact point is constant at z = 0. Setting E = T + V and u = 0, we compute E˙ = φ˙ T [¯I0 (ϕ)ϕ] ˙ φ˙ /2 + φ˙ T ¯I(ϕ)φ¨ + mg ϕ˙ eT c0 (ϕ) 3 h
0 = φ˙ T [I0 (ϕ)ϕ˙ + 2mJcωh (ϕ)T Jcω (ϕ)ϕ] ˙ φ˙ /2 h
˙ φ˙ ˙ − mφ˙ T Jcω (ϕ)T J˜cω (ϕ, φ) − φ˙ T Nφ (ϕ, φ) h h + mg ϕ˙ (eT3 c0h (ϕ) − eT2 ch (ϕ)) 0 = φ˙ T [I0 (ϕ)ϕ˙ + 2mJcωh (ϕ)T Jcω (ϕ)ϕ] ˙ φ˙ /2 h
0 − mφ˙ T Jcωh (ϕ)T [ψ˙ eb3 Jcωh (ϕ) + ϕ˙ Jcω (ϕ)] φ˙ h
˙ T eb3 (Jcω (ϕ) φ) ˙ = 0 = −mψ˙ (Jcωh (ϕ) φ) h
ξst = Rz (ψ)Rx (ϕ)Ry (θst ) e3 T
Since the contact point force f that ensures that the torus rolls without slipping does no work, we expect the unforced system (u = 0) to conserve energy. Although verification of energy conservation is, in principle, a straighforward calculation, there are some interesting features. First, using x˙ from (7), we find that
− φ˙ T [ϕI ˙ − e1 φ˙ T /2] I0 (ϕ)φ˙
The steering torque is applied, in spatial reference frame, about the virtual steering axis
T
Conservation of energy
[ ] ¯I(ϕ) = I(ϕ) + mJcω (ϕ)T Jcω (ϕ) , h h
˙ = Nφ (ϕ, φ) ˙ + mJcω (ϕ)T J˜cω (ϕ, φ) ˙ φ˙ , ¯ (ϕ, φ) N h h and ¯ G(ϕ) = mg c2 (ϕ) e1 . The translational motion is obtained by integrating the rolling constraint x˙ = −Rz (ψ)Jcω (ϕ) φ˙ (7) h
from which we see that, kinematically speaking, it is necessary to keep track of the heading angle ψ as well. It is not necessary to keep track of θ unless one desires, for example, to trace out the contact point path on the surface of the torus.
where the last equality follows from the fact that eb3 is skew symmetric. We have also made use of the fact that (1) implies that eT3 c0h (ϕ) = eT2 ch (ϕ) to eliminate the terms involving gravity. When the externally applied torque is not zero, the power supplied to the system is given by ¯ E˙ = φ˙ T τ = φ˙ T B(ϕ) u. Structure and detail of the dynamic equations The reduced dynamics (6) has a structure that is interesting and exploitable. Since Ix +mkch k2 0 0 ¯I(ϕ) = 0 Iy +mc24 Iy sϕ −mc2 c4 0 ∗ Ix c2ϕ +Iy s2ϕ +mc22 we see that its inversion does not mix up the apparent torques too much: the expressions for θ¨ and ψ¨ do not depend on expressions from the first ˙ are ¯φ (ϕ, φ) row (ϕ) ¨ of (6). The components of N given by ( ) ˙ ¯φ (ϕ, φ) N = m ρ c2 ϕ˙ 2 − (Iy cϕ + m c3 c4 ) θ˙ψ˙ 1 + ((Ix − Iy )cϕ sϕ + m c2 c3 ) ψ˙ 2
( ) ˙ ¯φ (ϕ, φ) N
[ = ϕ˙ C(ϕ) 23
θ˙ ψ˙
]
where C(ϕ) is the 2 × 2 matrix ] [ −mρsϕ c4 Iy cϕ +m(2c3 −ρ)c4 . Iy cϕ +mρsϕ c2 2(Iy −Ix )cϕ sϕ −m(2c3 −ρ)c2
θ˙2 ≥ 4mgc2 ((Ix −Iy )cϕ sϕ +mc2 c3 )/(Iy cϕ +mc3 c4 )2 . When the inequality is strict, there will be two ˙ for each roll angle. distinct turn rates, ψ,
Integrability of the undriven rolling torus It turns out that (almost) every solution of (6) with u(t) ≡ 0 that satisfies |ϕ(t)| < π/2, t ≥ 0, is, in fact, periodic. The key is to note that (6) with u = 0 has the form ϕ¨ = f (ϕ, ϕ, ˙ w) w˙ = ϕ˙ A(ϕ) w
every −θ˙ > 0. In contrast, when turning a certain level of speed is needed to keep the wheel at the chosen lean angle. Indeed, for nonzero roll ϕ, the requirement that ψ˙ in (10) be real requires that
(8)
where w ∈ Rp and A(ϕ) is a bounded p×p matrix, ˙ ψ) ˙ T in this case. with p = 2 and w = (θ,
Motorcycle wheel example Consider now a standard torus with constant profile radius ρ and center of mass located at (w0 , h0 ) = (0, r + ρ) so that the overall wheel radius is r + ρ and ch (ϕ) = (0, −rsϕ , ρ + rcϕ )T . Numerically, let r = 0.2, ρ = 0.1, Iy = 0.64, Ix = 0.38, m = 14.7, and g = 9.81. The values are 10
Let W (·) satisfy W 0 (ϕ) = A(ϕ) W (ϕ),
W (0) = I,
5
on the ϕ region of interest and define w(ϕ; ¯ ϕ0 , w0 ) = W (ϕ) W (ϕ0 )−1 w0 . It follows easily that every trajectory of (6) is determined by a trajectory of the (parametrized) two dimensional system
0
ϕ¨ = f (ϕ, ϕ, ˙ w(ϕ; ¯ ϕ0 , w0 )) .
−5
(9)
Furthermore, since (9) is conservative, the energy provides a first integral making stability analysis trivial and periodic solutions transparent. In particular, an equilirium point of (9) is unstable if the linearization has real eigenvalues and Lyapunov stable otherwise. Furthermore, these eigenvalues are the nonzeros eigenvalues of the linearization of (6). Equilibria We now consider the set of constant operating conditions (xe , ue ) for the reduced dynamics (6). Each constant state has the form xe = ˙ ψ) ˙ T . Since the second and third compo(ϕ, 0, θ, ˙ are zero by ϕ˙ = 0 and rows ¯ (ϕ, φ) nents of N ¯ 2 and 3 of B(ϕ) are linearly independent on (−π/2, π/2), we find that ue = 0 for every constant operating condition of the dynamic control ˙ and ψ˙ system (6). Thus, for equilibrium, ϕ, θ, must satisfy ((Ix − Iy )cϕ sϕ + m c2 c3 ) ψ˙ 2 − (Iy cϕ + m c3 c4 ) θ˙ψ˙ + mg c2 = 0
(10)
where we recall that c2 , c3 , and c4 are roll dependent components of the contact point vector: ch (ϕ) = (0, c2 , c2 )T and c4 = eT3 crh (ϕ). We are interested in the equilibrium motions with nonzero forward velocity, −θ˙ > 0. For the symmetric wheel case considered, straight running motion occurs with ϕ = 0 and ψ˙ = 0 for
0
50
100
150
200
250
300
350
Fig. 3. Equilibrium stability: linearized natural ˙ for roll angles ϕ from zero frequency versus |θ| to 40 degrees in increments of 2 degrees. representative of those for a rear motorcycle tire and wheel (Sharp et al., 2004). Figure 3 presents quantitative equilibrium stability results for the rolling wheel with these parameters. The natural frequency of an imaginary pair is plotted positive while the frequency of a hyperbolic pair is plotted negative, implying unstable. Blue indicates the ˙ condition, and more rapidly turning (higher ψ) green indicates the usual turning rate. Thus, at ϕ = 18deg, we see that the green branch is stable ˙ ≈ 107deg/sec to 118deg/sec, unstable from |θ| from 118 to 152deg/sec, and stable for higher speeds. Every blue branch is always stable. References Arnold, V. I., V. V. Kozlov and A. I. Neishtadt (2006). Mathematical Aspects of Classical and Celestial Mechanics. Springer-Verlag. Bloch, A. M. (2003). Nonholonomic Mechanics and Control. Springer-Verlag. Bullo, F. and A. D. Lewis (2004). Geometric Control of Mechanical Systems. Springer-Verlag. O’Reilly, O.M. (1996). The dynamics of rolling disks and sliding disks. Nonlinear Dynamics 10, 287–305. Sharp, R.S., S. Evangelou and D.J.N. Limebeer (2004). Advances in the modelling of motorcycle dynamics. Multibody System Dynamics 12(3), 251–283.