Analytical study of a hovering magnetic system (Levitron)

i ** Equation Chapter 1 Section 1 End of Degree Project Degree in Industrial Technologies Engineering Author: Alfonso García-Agúndez Blanco Tutors: Antonio González Fernández Manuel Toscano Jiménez Tutor: Javier Payán Somet Applied Physics III Department Higher Technical School of Engineering University of Seville Seville, 2017 Analytical study of a hovering magnetic system (Levitron) ii List of figures iii iii End of Degree Project Degree in Industrial Technologies Engineering Analytical study of a hovering magnetic system (Levitron) Author: Alfonso García-Agúndez Blanco Tutors: Prof. Dr. Antonio González Fernández Prof. Dr. Manuel Toscano Jiménez Applied Physics III Department Higher Technical School of Engineering University of Seville Seville, 2017 iv List of figures v v End of Degree Project: Analytical study of a hovering magnetic system (Levitron) Author: Alfonso García-Agúndez Blanco Tutors: Antonio González Fernández Manuel Toscano Jiménez El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: vi List of figures vii vii Agradecimientos Quiero dar las gracias, en primer lugar, a mis padres, que dan todo por mi bienestar y buscan por encima de cualquier cosa la felicidad de sus hijos. No hay palabras para agradecer todo vuestro esfuerzo diario. También agradecer a mis hermanos por compartir su día a día conmigo apoyándome y dándome cariño. A mis amigos, tanto a los de siempre como a los que he tenido la oportunidad de conocer en la Escuela, por ser mi segunda familia estos años. Especialmente, gracias a aquellos que han apostado por mi siempre y a los que se han interesado en el día a día sobre la marcha del trabajo, preguntándome e interesándose por el progreso. Gracias a Antonio por su total implicación en el trabajo desde el día que me ofreció el tema, por las reuniones de todos los jueves, por transmitirme conocimientos de diversas áreas (desde electromagnetismo hasta cuaterniones) y por el continuo interés mostrado. Gracias a Manuel, por transmitirme su entusiasmo también desde el primer día, por orientarme y aconsejarme en distintos puntos del trabajo y por confiar en mi, dándome tranquilidad y seguridad. Por último, quiero mencionar a Paula, que junto a mis padres, es la persona que ha seguido conmigo la evolución de este trabajo desde el inicio. Agradecerle sus consejos, sobre todo de programación, y el apoyo diario. viii List of figures ix ix Abstract (English) The scope of this project is to perform an analytical study of the hovering magnetic system, known as Levitron. Firstly, we will start with an introduction where we explain how the toy works, the typical elements of the kit and its history. In chapter 2, the analytical study starts: we begin describing the free motion, in absence of the external magnetic field and solely submitted to the action of gravity, using the Classic Euler angles. Nevertheless, the use of these angles leads to a singularity that we need to avoid, so in chapter 3 we present the Tait-Bryan angles to solve this problem. In chapter 4 we obtain the expressions for the magnetic field generated by the base, and in chapter 5 we outline Earnshaw’s Theorem and study the static stability by means of the potential energy of the system. In chapter 6, we derive the equations of motion using vector mechanics, and in chapter 7 the same is done using the Lagrangian and Hamiltonian formulation of the analytical mechanics. Once the equations of the system are obtained, we define nondimensional variables in chapter 8 to perform the numerical simulations. In chapter 9 we study the linear stability of the system, obtaining the stability region for which stable hovering is possible and the normal modes. In chapter 10 we analyze the linear and nonlinear coupling, and numerically simulate the trajectory of the spinning top considering different initial conditions that reproduce real situations when one plays with the toy. In chapter 11, we present the constants of the motion of the system, and in chapter 12 a simple model considering air friction is shown. Finally, chapter 13 contains some instructions to master the toy and we attach a flowchart where we describe the usual situations that a player has to face to achieve levitation. Lastly, we include three addendums: the first contains the experimental adjustment of the magnetic field generated by the base in the 0𝑍1 axis; the second is dedicated to the quaternions, containing the main properties of its algebra and the set of equations of the system in terms of the Euler parameters; and finally, in the third we study the dynamics of three systems that share some of the characteristics of the hovering magnetic device. xvi LIST OF FIGURES Figure 1-1: Levitron toy. 1 Figure 1-2: Versions of the toy with changing electromagnetic fields. 1 Figure 1-3: Contents of a Levitron kit. 2 Figure 1-4: View of the base magnet and legs. 2 Figure 1-5: Levitator system with rectangular base. 3 Figure 2-1: Classical Euler angles. 6 Figure 2-2: Important magnitudes in the free-motion analysis. 10 Figure 2-3: Angles 𝛼 and 𝛽. 10 Figure 2-4: Body and space cones. 11 Figure 3-1: Yaw-pitch-roll sequence used in navigation. 13 Figure 3-2: Rotation 21. 13 Figure 3-3: Rotation 32. 14 Figure 3-4: Rotation 43. 15 Figure 3-5: Successive rotations and different frames. 16 Figure 3-6: Fixed frame ‘1’ and body frame ‘4’ in the real toy. 17 Figure 4-1: Analogy between permanent magnet and coil. 21 Figure 4-2: Magnetic field created by a current loop in the 𝑂𝑍1-axis. 21 Figure 4-3: Magnetic field created by two concentric loops in the 𝑂𝑍1-axis. 22 Figure 4-4: Section of the magnetic base of the system. 22 Figure 4-5: Ilustration of the magnetic field created by the base and the top. 26 Figure 5-1: Potential energy for different values of 𝛽. 30 Figure 6-1: Free body diagram of the top. 34 Figure 6-2: Magnetic torque acting on the spinning top. 35 Figure 9-1: Equilibrium points, given by 𝑧=ze-axis 49 Figure 9-2: Interpretation of the complex variables. 57 Figure 9-3: Stability region. 62 Figure 9-4: Stability region in terms of the mass and rotational speed. 64 Figure 9-5: Campbell diagram for ze=3.1 65 Figure 9-6: Centered point in the stability region. 66 Figure 9-7: Natural frequencies for a centered value in the stability region. 66 Figure 9-8: Diagram Campbell for ze=3.19. 67 Figure 9-9: Evolution of 𝜆1, 𝜆2 and 𝜆3 for ze=3.19 68 Figure 9-10: Noncentered point in the stability region. 68 Figure 9-11: Natural frequencies for a noncentered point in the stability region. 69 Figure 9-12: Natural frequencies for a noncentered point in the stability region. 69 Figure 9-13: Identification of the motions using real and imaginary parts 71 List of figures xvii xvii Figure 9-15: Circular trajectories of the normal modes. 72 Figure 9-16: Circumferences of the four normal modes. 72 Figure 9-17: Circumference corresponding with the fast mode. 73 Figure 9-18: Vertical mode 73 Figure 9-19: Time evolution of 𝑧(𝑡) and 𝑝𝑧(𝑡). 74 Figure 9-20: Three-dimensional trajectories of the normal modes, adding a vertical perturbation. 74 Figure 9-21: Projection of the trajectories on the XZ plane: Lissajous curves. 75 Figure 9-22: Time evolution of 𝑥, 𝑦, 𝑝𝑥, 𝑝𝑦. 75 Figure 9-23: Time evolution of 𝜙,𝜃. 76 Figure 9-24: Variation of the normal modes radius with the spin velocity. 77 Figure 9-25: Horizontal perturbation of the top. 77 Figure 9-26: Trajectories of the linear system for horizontal perturbations. 78 Figure 9-27: Pitch perturbation of the top. 78 Figure 9-28: Trajectory of the linear system with pitch perturbation. 79 Figure 9-29: Trajectories for 𝜔<𝜔𝑚𝑖𝑛 and 𝜔=𝜔𝑚𝑖𝑛. 80 Figure 9-30: Trajectories within the stable rotational speed range. 80 Figure 9-31: Trajectories within the stable rotational speed range. 81 Figure 9-32: Trajectories for 𝜔=𝜔𝑚𝑎𝑥 and 𝜔>𝜔𝑚𝑎𝑥. 81 Figure 10-1: Circumferences in the nonlinear system. 83 Figure 10-2: Linear coupling: time evoluton of 𝑥, 𝑦. 84 Figure 10-3: Nonlinear coupling 84 Figure 10-4: Horizontal perturbation. 85 Figure 10-5: Trajectories of the nonlinear system with horizontal perturbation. 86 Figure 10-6: Pitch perturbation. 86 Figure 10-7: Trajectory of the nonlinear system with pitch perturbation. 86 Figure 10-8: Flower-shape trajectories. 87 Figure 10-9: Nonlinear trajectories near the equilibrium height. 88 Figure 10-10: Nonlinear trajectories near the equilibrium height. 88 Figure 10-11: Flower shape trajectories. 89 Figure 10-12: Time evolution of 𝑥 for 𝑧0=3 and𝑧0=3.08 89 Figure 10-13: Time evolution of 𝑦 for 𝑧0=3 and 𝑧0=3.08 90 Figure 10-14: Time evolution of 𝑧 for 𝑧0=3 and 𝑧0=3.08 90 Figure 10-15: Time evolution of 𝜙 for 𝑧0=3 and𝑧0=3.08 90 Figure 10-16: Time evolution of 𝜃 for 𝑧0=3 and 𝑧0=3.08 91 Figure 10-17: Time evolution of 𝑝𝑥 for 𝑧0=3 and 𝑧0=3.08 91 Figure 10-18: Time evolution of 𝑝𝑦 for 𝑧0=3 and𝑧0=3.08 91 Figure 10-19: Lower limit and stable trajectory of the nonlinear system. 92 Figure 10-20: Stable trajectories of the nonlinear model. Flower shape trajectories. 93 Figure 10-21: Stable trajectory and upper limit of the nonlinear model. 93 xviii Figure 11-1: Mechanical energy of the system. 95 Figure 12-1: Magnus Force on a rotating sphere. 100 Figure 12-2: Trajectory of the linear system. 105 Figure 12-3: Trajectory of the nonlinear system. 105 Figure 12-4: Trajectory considering air friction. 105 Figure 12-5: Trajectory and fall of the spinning top due to air friction. 106 Figure 13-1: Bubble levelling of the magnetic base. 107 Figure 13-2: Flowchart to achieve levitation. 109 Figure AdI-1: Experimental adjustment. 115 Figure AdII-1: Euler’s rotation theorem. 118 Figure AdII-2: Successive rotations in terms of the quaternions. 119 Figure AdIII-1: Martin Ekman: precession-nutation and polar motion. 122 Figure AdIII-2: Classical spinning top. 123 Figure AdIII-0-3: Intepretation of variable u. 125 Figure AdIII-4: Plot of 𝑓(𝑢) 126 Figure AdIII-5: Types of precession. 126 Figure AdIII-6: Earth motions in term of the Classic Euler angles. 127 Figure AdIII-7: Sun-Earth system. 128 Figure AdIII-8: Sun-Earth-Moon system. 131 Figure AdIII-9: Sun-Earth-Moon system. 131 Figure AdIII-10: Scheme of a Tokamak reactor. 133 Figure AdIII-11: Scheme with constant vertical and azimuthal magnetic fields. 134 Figure AdIII-12: Helical trajectory considering vertical constant magnetic field. 135 Figure AdIII-13: Trajectory considering vertical and azimuthal magnetic fields. 136 1 1 INTRODUCTION 1.1 What is a Levitron? The hovering magnetic device consists on a permanent magnetic top and its corresponding magnetic base plate, with a ring-shape or others alternate geometric configurations. In specific conditions and under top rotation the magnetic, gyroscopic and gravitational forces are balanced, and the top remains in levitation. The stabilizing rotational speed of the hovering top decreases due to air friction, becoming gradually slower, so that the hovering phenomenon fails within two minutes unless external power is supplied to sustain rotation. Figure 1-1: Levitron toy. This method, with moving permanent magnets, is quite distinct from other versions which use changing electromagnetic fields, levitating various items in a permanent way. This kind of device isn’t in the scope of this study. Figure 1-2: Versions of the toy with changing electromagnetic fields. The most known brand of spin-stabilized magnetic levitation device is Levitron. Levitron® is a brand of levitating toys and gifts in science and educational markets marketed by Creative Gifts Inc. and Fascination Toys & Gifts. The typical contents of a Levitron kit are: Introduction 2 2 Figure 1-3: Contents of a Levitron kit. 1.2 How does it work? In order to achieve the hovering, a plastic lifter plate is placed on top of the magnetic base, and the top is spun on the plate with a rotational speed between 25-50 rotations per second (1500-3000 rpm). If too slow, the top falls over and flies off sideways; if too fast it does not orient itself to follow the magnetic field lines as it moves, and slides off. Since it can be difficult to spin the top fast enough by hand, Creative Gifts makes a batterypowered, hand held device to spin the top with an electric motor. Next, the plate is lifted by hand until, if conditions are right, the top rises above it to an equilibrium point. The top must also be weighted with washers of various sizes supplied in the kit. If too heavy it will not rise above the plate; if too light it flies off. After a few minutes, the top falls when air friction slows it below a lower critical speed. Figure 1-4: View of the base magnet and legs. The top is a rotationally symmetric rigid body with mass m, polar and transversal moments of inertia 𝐼𝑍 and 𝐼𝑋, angular momentum with respect to the center of mass 𝐿 󰇍 𝐺 and angular speed 𝜔 󰇍 󰇍 whose centre of mass G is located at 𝑟 ={𝑥,𝑦,𝑧}. It can be regarded as a magnetic dipole with magnetic moment 𝜇 centered in 𝑟 and directed along its axis of symmetry. The gradients of the magnetic field 𝐵(𝑟) compensate for the gravitational force by generating a repulsive force that acts on 𝜇 in the presence of 𝐿 󰇍 𝐺 (whose gyroscopic effect prevents the top from overturning and falling) and must provide the mechanism for the top to levitate in a stable way above the base. 1.3 Levitron’s history According to [1], the first spin-stabilized permanent magnet levitation device was invented (1976) and patented (1983) by inventor Roy Harrigan, of Vermont. In the mid-1990s, Seattle entrepreneur Bill Hones, who was himself exploring the possibility of permanent magnet levitation, discovered Harrigan’s patent. Hones subsequently contacted Harrigan and later met him. Upon Hones’ request, Harrigan permitted him to borrow his Introduction 3 3 prototype with the understanding that they were entering into a business arrangement. Hones, with the help of his father, a Physicist and employee at Los Alamos National Laboratory, analyzed the physics of the prototype, and then filed for an “improvement patent”. In 1984, independent of Roy Harrigan, inventor Joseph Chieffo, of Pennsylvania, also discovered spin-stabilized magnetic levitation. Chieffo then developed his own spinstabilized magnetic levitation device and attempted to obtain a patent thereon. Employing an attorney to conduct a U.S. Patent and Trademark records search, Chieffo was informed of the existence of the Harrigan patent; he thus concluded his efforts to secure a patent. In a final assessment, Chieffo’s attorney noted that his device, although apparently unpatentable, could be marketed without infringing upon the Harrigan patent. In 1988, Chieffo marketed his device in kit form. Contrasting with the dished supporting magnet of the Harrigan invention, the base magnet of this latter device was rectangular and planar of upper surface, not unlike the base magnet of Hones’ later-patented device, the now-popular Physics toy known as the Levitron. Figure 1-5: Levitator system with rectangular base. In the 1990s, Michael and Karen Sherlock formed the company they named “UFO” in New Mexico to market the Levitron under an oral agreement in partnership with Hones’ company, Creative Gifts, Inc. Efforts to formalize the agreement in writing fell apart and grew acrimonious after UFO’s principals learned about the device’s earlier invention by Harrigan, and redesigned their website to incorporate the exposé-style article “THE HIDDEN HISTORY OF THE LEVITRON!” [2], which accused Hones of stealing the invention from Harrigan. Creative Gifts, in turn, filed a trademark infringement suit in United States District Court of New Mexico against UFO and its owners. At trial and on appeal to the Tenth Circuit Creative Gifts’ trademark claims were upheld, and all of UFO’s counterclaims were rejected after UFO, which had been representing itself as a “pro se” defendant, was sanctioned by the court for abuse of discovery. The appeals court, noting that UFO had submitted a one-page opening brief with no citations to the record or discussion of the relevant law, commented in its ruling, “they have shot themselves in the foot.” Note that in this moment (2017) it’s not possible to buy a Levitron device. This magnet-mechanical problem has been studied by several authors. M. V. Berry, in [3], describes the system using adiabatic approximations and provides the mass interval in which stable hovering is achieved. Genta et al. [4] provides an analytical expression to compute the lower and upper stable rotational speed bounds and Dullin et al. [5] and Gov et al. [6] find the stability range where hovering is possible. Moreover, A. T. Pérez et al. [7] analyze the misalignment between the magnetic dipole moment and the mechanical axis of the top, and achieve stable levitation for hours and even days using an alternating magnetic field. Introduction 4 4 2 Free motion. Description with Euler angles 5 5 2 FREE MOTION. DESCRIPTION WITH EULER ANGLES We begin our dynamical analysis describing the motion of the spinning top with the classical Euler angles, in absence of the external magnetic field and solely submitted to the action of gravity. Euler’s equations are particularised for the case of inertial symmetry in the body (like our spinning top), where two of the principal moments of inertia are the same. We then analyse the free-motion, the constants of the motion and give a geometric visualization by means of the axoids (body and space cones). 2.1. Euler angles To describe the orientation of the body, there are twelve choices in which no two adjacent rotation indices are the same. A sequence is given by three numbers, which represent each of the rotations (1 denotes a rotation about a X axis, 2 about a Y axis and 3 is a transformation about a Z axis). These twelve sets are called Euler angles sequences and they can be classified into two main categories. Depending on the application, one select the most suitable set, such that the sequence used gives a better physical visualization or leads to fewer singularities. - The first category group the transformations where the first and third indices are the same. - The second group consists of rotations where the first and the third indices are different. These are usually known as Tait-Bryan angles. The most common used Euler angle sequence is the 3-1-3 sequence, which is beforehand the most suitable choice to visualize the orientation of the body with respect the fixed coordinate system. In the 3-1-3 transformation, the rotation angles 𝜙, 𝜃 and 𝜓 are known as the precession, nutation and spin angles, respectively. 2.1.1. Rotation {21}: Precession angle 𝝓 Starting in the fixed frame, our first rotation is counterclockwise (if we look down the axis) through an angle 𝜙 about the 𝑍1-axis. The new frame is numbered as ‘2’ and its unit vectors are 𝑖 2, 𝑗 2 and 𝑘 󰇍 2. The relation between frames ‘1’ and ‘2’ is given by 𝑖 2=cos(𝜙)𝑖 1+sin(𝜙)𝑗 1 𝑗 2=−sin(𝜙)𝑖 1+cos(𝜙)𝑗 1 𝑘 󰇍 1=𝑘 󰇍 2 (2.1) 2.1.2. Rotation {32}: Tilt angle 𝜽 The second rotation is counterclockwise (if we look down the axis) through an angle θ about the 𝑋2-axis, changing from frame ‘2’ to ‘3’. The unit vectors are 𝑖 3, 𝑗 3 and 𝑘 󰇍 3. The transformation of coordinates can be represented as follows: 𝑖 3=𝑖 2 𝑗 3=cos(𝜃)𝑗 2+sin(𝜃)𝑘 󰇍 2 𝑘 󰇍 3=−sin(𝜃)𝑗 2+cos(𝜃)𝑘 󰇍 2 (2.2) 2 Free motion. Description with Euler angles 6 6 2.1.3. Rotation {43}: Roll angle 𝝍 The third rotation is counterclockwise (if we look down the axis) through an angle ψ about the 𝑍3-axis. The new frame is the body frame ‘4’ and the transformation is given by 𝑖 4=cos(𝜓)𝑖 3+sin(𝜓)𝑗 3 𝑗 4=−sin(𝜓)𝑖 3+cos(𝜓)𝑗 3 𝑘 󰇍 4=𝑘 󰇍 3 (2.3) In figure 2-1, we can see the successive rotations Figure 2-1: Classical Euler angles. The line of nodes is given by the instersection of the planes 𝝅𝒙𝟏𝒚𝟏(which contains the axis 𝑋1 and 𝑌1 of the fixed frame ‘1’) and 𝝅𝒙𝟒𝒚𝟒. To obtain the combined transformation from the frame ‘1’ to the local frame ‘4’, we combine the three previous rotations, {𝑋4}=𝑹𝑻{𝑋1} (2.4) with 𝑹𝑻 being 2 Free motion. Description with Euler angles 7 7 𝑹𝑻=𝑹𝝍𝑹𝜽𝑹𝝓 =(cos𝜓cos𝜙−cos𝜃sin𝜙sin𝜓 cos𝜓sin𝜙+cos𝜃cos𝜙sin𝜓sin𝜓sin𝜃 −sin𝜓cos𝜙−cos𝜃sin𝜙cos𝜓 −sin𝜓sin𝜙+cos𝜙cos𝜃cos𝜓 cos𝜓sin𝜃 sin𝜃sin𝜙 −sin𝜃cos𝜙 cos𝜃 ) (2.5) 2.2. Angular velocity The angular velocity of the top with respect the fixed frame ‘1’ is given by the composition of the previous rotations: 𝜔 󰇍 󰇍 =𝜔 󰇍 󰇍 41=𝜔 󰇍 󰇍 43+𝜔 󰇍 󰇍 32+𝜔 󰇍 󰇍 21=𝜙󰇗𝑘 󰇍 1+𝜃󰇗𝑗 2+𝜓󰇗𝑘 󰇍 3 (2.6) with its components being 𝜔𝑋=𝜃󰇗cos(𝜓)−𝜙󰇗 sin(𝜃)sin(𝜓) 𝜔𝑌=−𝜃󰇗sin(𝜓)−𝜙󰇗 sin(𝜃)cos(𝜓) 𝜔𝑍=𝜓󰇗+𝜙󰇗cos(𝜃) (2.7) Inverse relations give 𝜙󰇗=−𝜔𝑋sin(𝜓)+𝜔𝑌cos(𝜓) sin(𝜃) 𝜃󰇗=𝜔𝑋cos(𝜓)−𝜔𝑌sin(𝜓) 𝜓󰇗=𝜔𝑍+cos(𝜃) sin(𝜃)(𝜔𝑋sin(𝜓)+𝜔𝑌cos(𝜓)) (2.8) Note that the presence of the sine functions in the denominators of equations (2.8) leads to a singularity when 𝜃=0°. This value of 𝜃 corresponds with the spinning top in vertical position, with its local axis aligned with the 𝑂𝑍1-axis. This fact obscures the description of the problem with these angles, as we are not able to perform numerical integrations with 𝜃 around this value. Consequently, in chapter 3 we derive the matrix of rotation for the known ‘Gimbal set’, that corresponds with a yaw-pitch-rolling sequence. 2.3. Free motion with inertial symmetry The motion of the hovering spinning top, in absence of external magnetic field and in presence of gravity, is going to be analyzed. The translational equations are straightforward, as we are only considering the gravity force. Newton’s second law imply the conservation of the horizontal components of linear momentum (𝐶 ℎ𝑜𝑟=𝐶𝑇𝐸 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 ), and the zcomponent of the center of mass of the body describes a parabolic path. 3 Tait-bryan angles 14 14 The relation between the two systems due to this rotation is: 𝑖 2=𝑖 1 𝑗 2=cos(𝜙)𝑗 1+sin(𝜙)𝑘 󰇍 1 𝑘 󰇍 2=−sin(𝜙)𝑗 1+cos(𝜙)𝑘 󰇍 1 (3.1) with the inverse relation 𝑖 1=𝑖 2 𝑗 1=cos(𝜙)𝑗 2−sin(𝜙)𝑘 󰇍 2 𝑘 󰇍 1=sin(𝜙)𝑗 2+cos(𝜙)𝑘 󰇍 2 (3.2) Angular velocity of this rotation 𝜔 󰇍 󰇍 21 is 𝜔 󰇍 󰇍 21=𝜙󰇗 𝑖 1=𝜙󰇗𝑖 2 (3.3) 3.1.2. Rotation {32}: 𝜽 angle We arrive to frame ‘3’ making a counterclockwise rotation through an angle 𝜃 about 𝑂𝑌2=𝑂𝑌3 axis (line of nodes). Figure 3-3: Rotation {32}. The relation between both bases in this rotation is 𝑖 3=cos(𝜃)𝑖 2−sin(𝜃)𝑘 󰇍 2 𝑗 3=𝑗 2 𝑘 󰇍 3=sin(𝜃)𝑖 2+cos(𝜃)𝑘 󰇍 2 (3.4) with the inverse relation 3 Tait-bryan angles 15 15 𝑖 2=cos(𝜃)𝑖 3+sin(𝜃)𝑘 󰇍 3 𝑗 2=𝑗 3 𝑘 󰇍 2=−sin(𝜃)𝑖 3+cos(𝜃)𝑘 󰇍 3 (3.5) Angular velocity of this rotation, 𝜔 󰇍 󰇍 32 is given by 𝜔 󰇍 󰇍 32=𝜃󰇗𝑗 2=𝜃󰇗𝑗 3 (3.6) 3.1.3. Rotation {43}: 𝝍 angle The last rotation leads to the body-frame ‘4’ by means of a counterclockwise rotation 𝜓 about 𝑂𝑍3=𝑂𝑍4 axis. Figure 3-4: Rotation {43}. The relation between both bases is given by 𝑖 4=cos(𝜓)𝑖 3+sin(𝜓)𝑗 3 𝑗 4=−sin(𝜓)𝑖 3+cos(𝜓)𝑗 3 𝑘 󰇍 4=𝑘 󰇍 3 (3.7) and 𝑖 3=cos(𝜓)𝑖 4−sin(𝜓)𝑗 4 𝑗 3=sin(𝜓)𝑖 4+cos(𝜓)𝑗 4 𝑘 󰇍 3=𝑘 󰇍 4 (3.8) 3 Tait-bryan angles 16 16 The angular velocity corresponding with the last rotation is given by 𝜔 󰇍 󰇍 43 𝜔 󰇍 󰇍 43=𝜓󰇗𝑘 󰇍 3=𝜓󰇗𝑘 󰇍 4 (3.9) As we did with the 3-1-3 transformation, the objective is to obtain the combined rotation matrix 𝑹𝑻 which completely describes the relative orientation between the fixed and local frames: {𝑋1}=𝑹𝑻{𝑋4} (3.10) with 𝑹𝑻 𝑹𝑻=( cos𝜃cos𝜓 −cos𝜃sin𝜓 sin𝜃 sin𝜙sin𝜃cos𝜓+cos𝜙sin𝜓 −sin𝜙sin𝜃sin𝜓+cos𝜙cos𝜓 −sin𝜙cos𝜃 −cos𝜙sin𝜃cos𝜓+sin𝜙sin𝜓 cos𝜙sin𝜃sin𝜓+sin𝜙cos𝜓 cos𝜙cos𝜃) (3.11) Because of the rotational symmetry of the spinning top, we are going to use frame ‘3’ in the next chapters. Therefore, we give the expression of matrix 𝑹𝟑, which relates frames ’1’ and ‘3’: {𝑋1}=𝑹𝟑{𝑋3} (3.12) 𝑹𝟑=( cos𝜃 0 sin𝜃 sin𝜙sin𝜃 cos𝜙 −sin𝜙cos𝜃 −cos𝜙sin𝜃 sin𝜙 cos𝜙cos𝜃) (3.13) Figure 3-5: Successive rotations and different frames. 3 Tait-bryan angles 17 17 Figure 3-6: Fixed frame ‘1’ and body frame ‘4’ in the real toy. Despite the numerical advantage provided by the yaw-pitch-roll sequence, there are disadvantages. These angles provide a worse visualization of the problem than the 3-1-3 transformation. In the description with the ‘Classical Euler angles’, the angle between 𝑘 󰇍 4 and 𝑘 󰇍 1 is θ and it represents the nutation motion. In this case, visualization is not so clear and the nutation motion is given by a combination of 𝜙 and 𝜃. This can be seen computing the dot product 𝑘 󰇍 4·𝑘 󰇍 1: 𝑘 󰇍 4·𝑘 󰇍 1=cos𝜙cos𝜃=cos𝛿 (3.14) Therefore, the angle between both vectors is given by 𝛿=arccos (cos𝜙cos𝜃) (3.15) 3.2. Angular velocity 𝝎 󰇍 󰇍 󰇍 The angular velocity of the top with respect to the fixed frame ‘1’ is given by the composition of the previous rotations: 𝜔 󰇍 󰇍 =𝜔 󰇍 󰇍 41=𝜔 󰇍 󰇍 43+𝜔 󰇍 󰇍 32+𝜔 󰇍 󰇍 21=𝜓󰇗𝑘 󰇍 3+𝜃󰇗𝑗 3+𝜙󰇗𝑖 2 (3.16) Expressed in the frame ‘3’, we have 𝜔 󰇍 󰇍 |3=𝜔 󰇍 󰇍 41|3=𝜓󰇗𝑘 󰇍 3+𝜃󰇗𝑗 3+𝜙󰇗(cos(𝜃)𝑖 3+sin(𝜃)𝑘 󰇍 3)=𝜙󰇗cos(𝜃)𝑖 3+𝜃󰇗𝑗 3+(𝜓󰇗+𝜙󰇗sin(𝜃))𝑘 󰇍 3 (3.17) 3 Tait-bryan angles 18 18 or, in components 𝜔𝑋=𝜙󰇗cos(𝜃) 𝜔𝑌=𝜃󰇗 𝜔𝑍=𝜓󰇗+𝜙󰇗sin(𝜃) (3.18) The inverse relation provides 𝜃󰇗=𝜔𝑌 𝜙󰇗=𝜔𝑋 cos(𝜃) 𝜓󰇗=𝜔𝑍−𝜔𝑋tan(𝜃) (3.19) Note that equations (3.19) are regular for 𝜃=0°. Arrived to this point, it is important to make a clarification between 𝜔 󰇍 󰇍 31 and 𝜔 󰇍 󰇍 |3, because the nomenclature might be confusing: - On the one hand, 𝜔 󰇍 󰇍 |3=𝜔 󰇍 󰇍 41|3 is the angular velocity of the top with respect to the frame ‘1’, expressed in the frame ‘3’, as we have just seen. - On the other hand, 𝜔 󰇍 󰇍 31is the angular velocity of the frame ‘3’ with respect to the fixed frame ‘1’. Vectors at base 3 move with respect to the fixed frame with an angular velocity which is different to that of the top, being the difference the angular velocity of the top around its own axis. Its expression will be necessary in chapter 6 and it is given by 𝜔 󰇍 󰇍 31=𝜔 󰇍 󰇍 32+𝜔 󰇍 󰇍 21=𝜃󰇗𝑗 3+𝜙󰇗𝑖 2 (3.20) To take advantage of the rotational symmetry of the top, 𝜔 󰇍 󰇍 31 is going to be expressed precisely at base ‘3’: 𝜔 󰇍 󰇍 31|3=Ω 󰇍 󰇍 =𝜙󰇗cos(𝜃)𝑖 3+𝜃󰇗𝑗 3+𝜙󰇗sin(𝜃)𝑘 󰇍 3 (3.21) The components of 𝜔 󰇍 󰇍 31can be written in terms of the components of the angular velocity 𝜔𝑋, 𝜔𝑌 of the top as follows: Ω𝑋=𝜙󰇗cos(𝜃)=𝜔𝑋 Ω𝑌=𝜃󰇗=𝜔𝑌 Ω𝑍=𝜙󰇗sin(𝜃)=𝜔𝑋tan(𝜃) (3.22) 3.3. Kinetic moment 𝑳 󰇍 󰇍 𝑮 Assuming rotational symmetry of the spinning top about its axis (which coincides with the direction of the magnetic pole), the inertia moment tensor of the dipole can be written, in a system of coordinates that moves with the top, as: 3 Tait-bryan angles 19 19 𝐼=(𝐼𝑋0 0 0 𝐼𝑋0 0 0 𝐼𝑍) Due to the rotational symmetry of the top, this expression is also valid in the frame ‘3’. (3.23) Assimilating the top to a disk of mass m and radius r, the inertia tensor can be written as 𝐼=14𝑚𝑟2(1 0 0 0 1 0 0 0 2) (3.24) Thus, the kinetic moment is related with the angular velocity through the inertia tensor. Its expression, in frame ‘3’, is 𝐿 󰇍 𝐺|3=𝐼|3·𝜔 󰇍 󰇍 |3 (3.25) and we get 𝐿 󰇍 𝐺|3=(𝐼𝑋0 0 0 𝐼𝑋0 0 0 𝐼𝑍)( 𝜙󰇗cos(𝜃) 𝜃󰇗 𝜓󰇗+𝜙󰇗sen(𝜃))=( 𝐼𝑋𝜙󰇗cos(𝜃) 𝐼𝑋𝜃󰇗 𝐼𝑍(𝜓󰇗+𝜙󰇗sin(𝜃))) (3.26) 3.4. Kinetic energy 𝑻 To compute the total kinetic energy, we have to consider the translational and rotational contributions 𝑇=𝑇𝑇+𝑇𝑅 (3.27) The translational kinetic energy 𝑇𝑇 is given by 𝑇𝑇=12𝑚𝑣𝐺2=12𝑚(𝑥󰇗2+𝑦󰇗2+𝑧󰇗2) (3.28) On the other hand, the rotational kinetic energy 𝑇𝑅 is 𝑇𝑅=12𝜔 󰇍 󰇍 𝑇·𝐼·𝜔 󰇍 󰇍 =12 𝜔 󰇍 󰇍 ·𝐿 󰇍 𝐺=12(𝜙󰇗cos(𝜃) 𝜃󰇗 𝜓󰇗+𝜙󰇗sin(𝜃))·( 𝐼𝑋𝜙󰇗cos(𝜃) 𝐼𝑋𝜃󰇗 𝐼𝑧(𝜓󰇗+𝜙󰇗sin(𝜃))) (3.29) obtaining 𝑇𝑅=12𝐼𝑋(𝜃󰇗2+𝜙󰇗2cos2(𝜃))+12𝐼𝑍(𝜓󰇗+𝜙󰇗sin(𝜃))2 (3.30) Therefore, we have 𝑇=12𝑚(𝑥󰇗2+𝑦󰇗2+𝑧󰇗2)+12𝐼𝑋(𝜃󰇗2+𝜙󰇗2cos2(𝜃))+12𝐼𝑍(𝜓󰇗+𝜙󰇗sin(𝜃))2 (3.31) 3 Tait-bryan angles 20 20 4 Magnetic field 21 21 4 MAGNETIC FIELD Before obtaining the expression of the magnetic field created by the base, some of the key aspects of the theory of the magnetic fields are going to be remembered: - The magnetic field created by a bar magnet form lines without extremes. The direction of the magnetic field is taken to be outward from the north pole and in to the south pole of the magnet. - The magnetic field produced by a circular bar magnet is similar to that of a cylindrical coil with the same diameter. Therefore, if the thickness of the bar magnet tends to zero and becomes a magnetized disk, its magnetic field is equivalent to the field to that of a current loop. Figure 4-1: Analogy between permanent magnet and coil. 4.1. Magnetic field at Z-axis In this way, if the base were a magnetized disk we could consider it as a closed current-carrying segment. Since it has a circular bore can be modelled as a superposition of two concentrical circular loops. The expression of the magnetic field produced by a line current is given by the Biot-Savart law: 𝐵 󰇍 (𝑟 )=𝜇0 4𝜋∮𝐼𝑑𝑟 ×(𝑟 −𝑟 ′) |𝑟 −𝑟 ′|3 (4.1) If we apply the law on the centerline of a current loop and integrate the z-component, we get: 𝐵 󰇍 0(𝑧)=𝜇0𝐼 2𝑅2 (𝑅2+𝑧2)3/2 𝑘 󰇍 1 (4.2) Figure 4-2: Magnetic field created by a current loop in the Z-axis. 4 Magnetic field 22 22 In the case of the Levitron, as we have seen in chapter 1, the base magnet consists of an annulus shaped magnetized block. In this way, it can be modeled as two concentric current loops of radius a and b, being a the inner radius and b the outer radius of the annulus. To take into account the central hole of the annulus (which does not contribute to the magnetic field along the axis), the outer ring is going to be modeled as a loop that transports a counterclockwise current intensity 𝐼, whereas the inner ring is going to be considered as a loop along which a clockwise current intensity 𝐼 flows. Figure 4-3: Magnetic field created by two concentric loops in the Z-axis. Considering a generic point P of the Z-axis, the use of the right-hand rule allows us to see that the contributions of the two current loops to the magnetic field at this point have opposite senses, so that the inner loop weakens the magnetic field of the outer ring (in the same way the central hole is a desmagnetized area in the base magnet). Thus, the analytical expression for the magnetic field along the Z-axis is given by: 𝐵 󰇍 0(𝑧)=𝜇0𝐼 2(𝑏2 (𝑏2+𝑧2)3/2 −𝑎2 (𝑎2+𝑧2)3/2 )𝑘 󰇍 1 (4.3) Figure 4-4: Section of the magnetic base of the system. 4 Magnetic field 23 23 In order to study the behavior of the system, it is necessary to experimentally measure the magnetic field along the Z-axis and try to make a regression adjustment between the experimental data and the analytical expression given by the theory. We have to take into account some details to do this: - First of all, the base is not a perfect annulus. The central hole can be considered as a circle of radius 𝑎=25 mm, but we cannot measure the radius of the outer loop b directly on the base. - The second thing to consider is that our analytical model consists of two concentric rings placed on a plane, whereas the base magnet is a block with height. If we want to make an adjustment between the experimental data and the analytical model, it will be necessary a variable change working with 𝑧−𝑧0. The magnetic field produced by the base was measured in the work ‘Efectos de pequeñas corrientes de aire sobre un Levitron’ [8] using a teslameter. A table with the experimental data and the detailed attainment of the parameters for the adjustment are attached at addendum I. The expression of 𝐵 󰇍 0(𝑧) is therefore 𝐵 󰇍 0(𝑧)=𝐴(𝑏2 (𝑏2+(𝑧−𝑧0)2)3/2 −𝑎2 (𝑎2+(𝑧−𝑧0)2)3/2 )𝑘 󰇍 1 with 𝑏=67.6864 mm, 𝑧0=−11.4067 mm, 𝐴=1.9608· mT·m (4.4) 4.2. Magnetic field near the Z-axis The magnetic field caused by the magnetized base acting on the dipole fulfils Maxwell’s equations for magnetostatics: ∇·𝐵 󰇍 =0 ∇×𝐵 󰇍 =0 󰇍 (4.5) Moreover, we assume cylindrical symmetry of the magnetic field about the 𝑂𝑍1 axis. Using cylindrical coordinates, it can be written as: 𝐵 󰇍 =𝐵𝜌(𝜌,𝑧)𝑢 󰇍 𝜌+𝐵𝑧(𝜌,𝑧)𝑘 󰇍 1 (4.6) Where 𝜌=√𝑥2+𝑦2 𝑢 󰇍 𝜌=𝑥𝜌𝑖 1+𝑦𝜌𝑗 1 (4.7) The magnetostatic equations lead to the following set: 1𝜌𝜕 𝜕𝜌(𝜌𝐵𝜌)+𝜕𝐵𝑧 𝜕𝑧=0 𝜕𝐵𝜌 𝜕𝑧−𝜕𝐵𝜌 𝜕𝑧=0 (4.8) As we know the magnetic field along the Z-axis of the system, we can integrate these equations in order to determine the field at points near the axis of the system. The following procedure is going to be performed: 1. We start with the expression of the magnetic field created by the base along the Z-axis, given by eq. (4.4). Because of the symmetry of the system, this field can only be vertical 𝐵 󰇍 (0,𝑧)=𝐵0(𝑧)𝑘 󰇍 (4.9) 2. Considering points near the axis, i.e, with 𝜌≪1, at these points 𝐵𝑧 approximately equals to the field at the points of the Z-axis, so we can do 5 Study of the static stability 30 30 Figure 5-1: Potential energy for different values of 𝛽. We can look for the minimum value of 𝛽 for which a stable static equilibrium position can exist. Computing the first and second derivatives of 𝑈∗ with respect to 𝑧∗, this value corresponds to the curve where the minimum and the maximum in the potential energy coalesce, becoming this stationary point also an inflexion point. 𝑑𝑈∗ 𝑑𝑧∗=1+𝛽𝐵0∗′(𝑧∗)=1+3𝛽( 1 (1+𝑧∗2)52−𝛾2 (𝛾2+𝑧∗2)52) (5.13) 𝑑2𝑈∗ 𝑑𝑧∗2=𝛽𝐵0∗′′(𝑧∗)=3𝛽( 1 (1+𝑧∗2)52−𝛾2 (𝛾2+𝑧∗2)52+5𝛾2𝑧∗2 (𝛾2+𝑧∗2)72−5𝑧2 (1+𝑧∗2)72) (5.14) and solving the set of equations { 𝑑𝑈∗ 𝑑𝑧∗=0 𝑑2𝑈∗ 𝑑𝑧∗2=0 (5.15) we get the value of 𝑧∗ for which the curve 𝑈∗(𝛽,𝑧∗) has an inflexion point. Parameter 𝛽 verifies the relation 5 Study of the static stability 31 31 𝛽=− 1 𝐵0∗′(𝑧∗) (5.16) The solutions are 𝑧∗=2.9235 and 𝛽𝑚𝑖𝑛=30.4311. In figure 5-2 we can see the potential energy for 𝛽𝑚𝑖𝑛. This minimum value of z corresponds to 7.31 cm. Below this height levitation is not possible. Figure 5-2: Potential energy for 𝛽𝑚𝑖𝑛 From the definition of the parameter 𝛽, we can see that 𝛽𝑚𝑖𝑛 provides the minimum ratio 𝜇 𝑚 in order to have a stable equilibrium. Moreover, if we assume a constant value of the magnetic moment 𝜇, one can see that the only variable in the definition of 𝛽 is the mass of the spinning top 𝑚. This mass can be changed by the player adding the washers described in the Introduction, so that a correct election is of great importance to achieve the hovering of the top. 5 Study of the static stability 32 32 6 Analysis with vector mechanics 33 33 6 ANALYSIS WITH VECTOR MECHANICS The motion of the hovering top is governed by the equation of the dynamics, that is, Newton’s second law and Euler’s second law. 𝑚𝑎 𝐺=𝐹 d𝐿 󰇍 𝐺 d𝑡=𝑀 󰇍 󰇍 𝐺 (6.1) where G is the center of mass of the dipole. 6.1. Translational motion: Newton’s second law As the top hovers over the magnetized plate, it is submitted to the gravity and magnetic force: 𝐹 =𝑚𝑔 +∇(𝜇 ·𝐵 󰇍 ) (6.2) The application of Newton’s second law provides the following vectorial equation: 𝐹 =𝑚𝑔 +∇(𝜇 ·𝐵 󰇍 )=𝑚𝑎 𝐺 (6.3) To compute the magnetic force, we need the dot product 𝜇 ·𝐵 󰇍 . As we saw in the previous chapters, 𝜇 has the same direction as the unit vector 𝑘 󰇍 4, pointing downwards. This unit vector, expressed in frame ‘1’, is given by: 𝑘 󰇍 4=( sin(𝜃) −sin(𝜙)cos(𝜃) cos(𝜙)cos(𝜃)) (6.4) The dot product 𝜇 ·𝐵 󰇍 is therefore 𝜇 ·𝐵 󰇍 =−𝜇(𝐵𝑥sin(𝜃)−𝐵𝑦sin(𝜙)cos(𝜃)+𝐵𝑧cos(𝜃)cos(𝜙)) (6.5) 𝐹 =𝑚𝑔 +∇(𝜇 ·𝐵 󰇍 )=( 00 −mg)−𝜇 ( 𝜕𝐵𝑥 𝜕𝑥sin(𝜃)−𝜕𝐵𝑦 𝜕𝑥sin(𝜙)cos(𝜃)+𝜕𝐵𝑧 𝜕𝑥cos(𝜙)cos(𝜃) 𝜕𝐵𝑥 𝜕𝑦sin(𝜃)−𝜕𝐵𝑦 𝜕𝑦sin(𝜙)cos(𝜃)+𝜕𝐵𝑧 𝜕𝑦cos(𝜙)cos(𝜃) 𝜕𝐵𝑥 𝜕𝑧sin(𝜃)−𝜕𝐵𝑦 𝜕𝑧sin(𝜙)cos(𝜙)+𝜕𝐵𝑧 𝜕𝑦cos(𝜙)cos(𝜃) ) (6.6) Thus, the equilibrium of forces applied to the top leads to the following three scalar equations: 𝑚𝑥󰇘=−𝜇(𝜕𝐵𝑥 𝜕𝑥sin(𝜃)−𝜕𝐵𝑦 𝜕𝑥sin(𝜙)cos(𝜃)+𝜕𝐵𝑧 𝜕𝑥cos(𝜙)cos(𝜃)) (6.7) 𝑚𝑦󰇘=−𝜇(𝜕𝐵𝑥 𝜕𝑦sin(𝜃)−𝜕𝐵𝑦 𝜕𝑦sin(𝜙)cos(𝜃)+𝜕𝐵𝑧 𝜕𝑦cos(𝜙)cos(𝜃)) (6.8) 6 Analysis with vector mechanics 34 34 𝑚𝑧󰇘=−𝑚𝑔−𝜇(𝜕𝐵𝑥 𝜕𝑧sin(𝜃)−𝜕𝐵𝑦 𝜕𝑧sin(𝜙)cos(𝜙)+𝜕𝐵𝑧 𝜕𝑦cos(𝜙)cos(𝜃)) (6.9) Figure 6-1: Free body diagram of the top The previous expressions can be rewritten as a first order set using the three components of the center of mass velocity: {𝑥󰇗=𝑣𝑥 𝑦󰇗=𝑣𝑦 𝑧󰇗=𝑣𝑧 (6.10) { 𝑣󰇗𝑥=1 𝑚𝜕 𝜕𝑥(𝐵𝑥𝜇𝑥+𝐵𝑦𝜇𝑦+𝐵𝑧𝜇𝑧) 𝑣󰇗𝑦=1 𝑚𝜕 𝜕𝑦(𝐵𝑥𝜇𝑥+𝐵𝑦𝜇𝑦+𝐵𝑧𝜇𝑧) 𝑣󰇗𝑧=−𝑔+1 𝑚𝜕 𝜕𝑧(𝐵𝑥𝜇𝑥+𝐵𝑦𝜇𝑦+𝐵𝑧𝜇𝑧) (6.11) 6 Analysis with vector mechanics 35 35 6.2. Rotational motion: Modified Euler’s equations The resultant torque acting on the center of mass is only caused by the magnetic field, because the weight of the top is applied at G. The magnetic torque can be related to the characteristics of the dipole by magnetic dipole moment define above. Watching the figure, we can see that the torque exerted by the magnetic field 𝑀 󰇍 󰇍 𝐺 tends to align the magnetic dipole moment 𝜇 with the magnetic field 𝐵 󰇍 , so this represents its lowest energy configuration. Figure 6-2: Magnetic torque acting on the spinning top. Considering the torque as a vector quantity, it can be written as the vector product: 𝑀 󰇍 󰇍 𝐺=𝜇 ×𝐵 󰇍 (6.12) Applying the angular momentum theorem, we have: d𝐿 󰇍 𝐺 d𝑡=𝑀 󰇍 󰇍 𝐺 (6.13) Because of the rotational symmetry of the spinning top, it is highly recommendable the expression of the vector equation using frame ‘3’, which allows us to work with an easier set of equations. When we take advantage of this rotational symmetry, these equations are known as the Modified Euler’s equations d𝐿 󰇍 𝐺 d𝑡|1=d𝐿 󰇍 𝐺 d𝑡|3+𝜔 󰇍 󰇍 31×𝐿 󰇍 𝐺=𝜇 |3×𝐵 󰇍 |3 (6.14) To compute the absolute time derivative of the angular momentum of the top, we need to know the time derivative of 𝐿 󰇍 𝐺 in the frame ‘3’ and the angular velocity of this frame with respect to the fixed frame ‘1’. Both vectors were computed in chapters 3.3 and 3.1, respectively. d𝐿 󰇍 𝐺 d𝑡|3=(𝐼𝑋(𝜙󰇘cos(𝜃)−𝜙󰇗𝜃󰇗sin(𝜃)) 𝐼𝑋𝜃󰇗 𝐼𝑍(𝜓󰇘+𝜙󰇘sin(𝜃)+𝜙󰇗𝜃󰇗cos(𝜃))) (6.15) 6 Analysis with vector mechanics 36 36 The cross product between both vectors is given by 𝜔 󰇍 󰇍 31|3×𝐿 󰇍 𝐺|3=| 𝑖 3𝑗 3𝑘 󰇍 3 𝜙󰇗cos(𝜃)𝜃󰇗𝜙󰇗sin(𝜃) 𝐼𝑋𝜙󰇗cos(𝜃)𝐼𝑋𝜃󰇗𝐼𝑍(𝜓󰇗+𝜙󰇗sin(𝜃))|=( 𝐼𝑍𝜃󰇗(𝜓󰇗+𝜙󰇗sin(𝜃))−𝐼𝑋𝜙󰇗𝜃󰇗sin(𝜃) 𝐼𝑋𝜙󰇗2sin(𝜃)cos(𝜃)−𝐼𝑍𝜙󰇗(𝜓󰇗+𝜙󰇗sin(𝜃))cos(𝜃) 0) (6.16) On the other hand, the magnetic torque is 𝑀 󰇍 󰇍 𝐺=𝜇 |3×𝐵 󰇍 |3=−𝜇𝑘 󰇍 3×𝐵 󰇍 |3=−𝜇𝐵𝑌𝑖 3−𝜇𝐵𝑋𝑗 3 (6.17) where 𝐵𝑋 and 𝐵𝑌 are the components in the frame ‘3’. Projecting in this frame, we get 𝐵𝑋=𝐵 󰇍 ⋅𝑖 3=𝐵 󰇍 ·(cos(𝜃)𝑖 1−sin(𝜃)(−sin(𝜙)𝑗 1+cos(𝜙)𝑘 󰇍 1)) =𝐵𝑥cos(𝜃)+𝐵𝑦sin(𝜃)sin(𝜙)−𝐵𝑧sin(𝜃)cos(𝜙) (6.18) 𝐵𝑌=𝐵 󰇍 ⋅𝑗 3=𝐵 󰇍 ·(cos(𝜙)𝑗 1+sin(𝜙)𝑘 󰇍 1)=𝐵𝑦cos(𝜙)+𝐵𝑧sin(𝜙) (6.19) where 𝐵𝑥, 𝐵𝑦 and 𝐵𝑧 are the components of the field derived in chapter 4.2. The torque is therefore 𝑀 󰇍 󰇍 𝐺=( 𝜇(𝐵𝑦cos(𝜙)+𝐵𝑧sin(𝜙)) 𝜇(𝐵𝑧cos(𝜙)sin(𝜃)−𝐵𝑥cos(𝜃)−𝐵𝑦sin(𝜙)sin(𝜃)) 0) (6.20) Substituting these expressions in eq. (6.14), we get the following second order differential equation set 𝐼𝑋𝜙󰇘cos(𝜃)+(𝐼𝑍−2𝐼𝑋)𝜙󰇗𝜃󰇗sin(𝜃)+𝐼𝑍𝜃󰇗𝜓󰇗−𝜇(𝐵𝑦cos(𝜙)+𝐵𝑧sin(𝜙))=0 (6.21) 𝐼𝑋𝜃󰇘+(𝐼𝑋−𝐼𝑍)𝜙󰇗2sin(𝜃)cos(𝜃)−𝐼𝑍𝜙󰇗𝜓󰇗cos(𝜃) +𝜇(𝐵𝑥cos(𝜃)+𝐵𝑦sin(𝜙)sin(𝜃)−𝐵𝑧cos(𝜙)sin(𝜃))=0 (6.22) 𝑑 𝑑𝑡(𝐼𝑧(𝜓󰇗+𝜙󰇗sin(𝜃)))=0 (6.23) As it was done with the translational motion, we can rewrite the second order differential equations set as a first order set using the components of the angular velocity 𝜔𝑋, 𝜔𝑌, 𝜔𝑍. d𝐿 󰇍 𝐺 d𝑡|1=d𝐿 󰇍 𝐺 d𝑡|3+𝜔 󰇍 󰇍 31×𝐿 󰇍 𝐺=𝐼𝑋(𝜔󰇗𝑋𝑖 3+𝜔󰇗𝑌𝑗 3)+𝐼𝑍ω󰇗𝑍𝑘 󰇍 3+|𝑖 3𝑗 3𝑘 󰇍 3 𝜔𝑋𝜔𝑌𝜔𝑋tan(𝜃) 𝐼𝑋𝜔𝑋𝐼𝑋𝜔𝑌𝐼𝑍𝜔𝑍| (6.24) Solving and separating in components, we have 6 Analysis with vector mechanics 37 37 {𝐼𝑋 𝜔󰇗𝑋+𝜔𝑌(𝐼𝑍𝜔𝑍−𝐼𝑋𝜔𝑋tan(𝜃)) = 𝑀𝑋 𝐼𝑋 𝜔󰇗𝑌−𝜔𝑋(𝐼𝑍𝜔𝑍−𝐼𝑋𝜔𝑋tan(𝜃)) = 𝑀𝑌 𝐼𝑍ω󰇗𝑍= 𝑀𝑍 (6.25) The first-order set is therefore { 𝜙󰇗=𝜔𝑋 cos(𝜃) 𝜃󰇗=𝜔𝑌 𝜓󰇗=𝜔𝑍−𝜔𝑋tan(𝜃) (6.26) {𝐼𝑋 𝜔󰇗𝑋+𝜔𝑌(𝐼𝑍𝜔𝑍−𝐼𝑋𝜔𝑋tan(𝜃)) = 𝜇(𝐵𝑦cos(𝜙)+𝐵𝑧sin(𝜙)) 𝐼𝑋 𝜔󰇗𝑌−𝜔𝑋(𝐼𝑍𝜔𝑍−𝐼𝑋𝜔𝑋tan(𝜃)) = −𝜇(𝐵𝑥cos(𝜃)+𝐵𝑦sin(𝜃)sin(𝜙)−𝐵𝑧sin(𝜃)cos(𝜙)) 𝐼𝑍ω󰇗𝑍= 0 (6.27) Taking advantage of the fact that the torque equation has been used to describe the motion of the spinning dipole about its center of mass, we are going to obtain briefly the spin precession. In the limit of large spin angular velocity ω, the angular momentum 𝐿 󰇍 𝐺 may be approximated by 𝐿 󰇍 𝐺=𝐼𝑍𝜔 󰇍 󰇍 =𝐼𝑍𝜔𝜇 𝜇 (6.28) Then, the torque equation d𝐿 󰇍 𝐺 d𝑡=𝜇 ×𝐵 󰇍 (6.29) can be rewritten as d𝜇 d𝑡=−𝜇𝐵 󰇍 𝐼𝑍𝜔×𝜇 (6.30) Thus, one can see that the motion of the dipole relative to its center of mass consists of a precession about the local direction of the magnetic field with an angular velocity 𝛺𝑝𝑟𝑒𝑐=−𝜇𝐵 𝐼𝑍𝜔 (6.31) 6 Analysis with vector mechanics 38 38 7 Analysis with analytical mechanics 39 39 7 ANALYSIS WITH ANALYTICAL MECHANICS In the previous chapters a description of the system using vector quantities of motion was made using Newton and Euler laws. By contrast, there is an alternative formulation of classical mechanics that uses scalar properties of motion (such as its kinetic and potential energy) to describe the system as a whole: analytical mechanics. We are going to derive the equations of motion of the Levitron using the two dominant branches of analytical mechanics: the Lagrangian formulation, using generalized coordinates and the corresponding generalized velocities, and the Hamiltonian mechanics, which uses coordinates and corresponding conjugate momenta. 7.1. Lagrange’s equations Lagrange equations for a system without friction are: 𝑑 𝑑𝑡(𝜕ℒ 𝜕𝑞󰇗𝑘)−𝜕ℒ 𝜕𝑞𝑘=0 (7.1) where 𝑞𝑘 denotes each of the generalized coordinates, i.e 𝑞𝑘={𝑥,𝑦,𝑧,𝜙,𝜃,𝜓 } and ℒ is the Lagrangian of the system, which is given by ℒ=T−U (7.2) Substituting the expressions of T and U obtained in the previous chapters, we have ℒ=12𝑚(𝑥󰇗2+𝑦󰇗2+𝑧󰇗2)+12𝐼𝑋(𝜃󰇗2+𝜙󰇗2cos2(𝜃))+12𝐼𝑍(𝜓󰇗+𝜙󰇗sin(𝜃))2−𝑚𝑔𝑧− 𝜇 (𝐵𝑥sin(𝜃)−𝐵𝑦sin(𝜙)cos(𝜃)+𝐵𝑧cos(𝜙)cos(𝜃)) (7.3) Computing the partial derivatives for each of the coordinates, 𝒒=𝒙 𝜕ℒ 𝜕𝑥󰇗=𝑚𝑥󰇗→𝑑 𝑑𝑡(𝜕ℒ 𝜕𝑦󰇗)= 𝑚𝑥󰇘 (7.4) 𝜕ℒ 𝜕𝑥=−𝜇(𝜕𝐵𝑥 𝜕𝑥sin(𝜃)−𝜕𝐵𝑦 𝜕𝑥cos(𝜙)sin(𝜃)+𝜕𝐵𝑧 𝜕𝑥cos(𝜙)cos(𝜃)) (7.5) 𝒒=𝒚 𝜕ℒ 𝜕𝑦󰇗=𝑚𝑦󰇗→𝑑 𝑑𝑡(𝜕ℒ 𝜕𝑦󰇗)=𝑚𝑦󰇘 (7.6) 𝜕ℒ 𝜕𝑦=−𝜇 (𝜕𝐵𝑥 𝜕𝑦sin(𝜃)−𝜕𝐵𝑦 𝜕𝑦cos(𝜙)sin(𝜃)+𝜕𝐵𝑧 𝜕𝑦cos(𝜙)cos(𝜃)) (7.7) 8 Nondimensionalization and summary of equations 46 46 have a nondimensional expression for the potential energy: 𝛾=𝑏𝑎 𝐵0=𝐴𝑎 𝛽=𝜇𝐵0 𝑚𝑔𝑎 𝜎=𝐼𝑧 𝐼𝑥=2 𝑅=𝑚𝑎2 𝐼𝑥=𝑚𝑎2 14𝑚𝑟2=16 (8.3) To obtain the value of 𝑅, it has been used the relation 𝑎𝑟=2, obtained by simple inspection of the toy. 8.2. Summary of equations 8.2.1. Vectorial mechanics x󰇗∗=𝑣𝑥∗ (8.4) y󰇗∗=𝑣𝑦∗ (8.5) z󰇗∗=𝑣𝑧∗ (8.6) 𝜙󰇗∗=𝜔𝑥∗ cos(𝜃∗) (8.7) 𝜃󰇗∗=𝜔𝑦∗ (8.8) 𝜓󰇗∗=𝜔𝑧∗−𝜔𝑥∗tan(𝜃∗) (8.9) 𝑣𝑥∗󰇗=𝛽(12𝑑𝐵0∗ 𝑑𝑧∗sin(𝜃∗)+𝑥∗ 2 𝑑2𝐵0∗ 𝑑𝑧∗2cos(𝜙∗)cos(𝜃∗)) (8.10) 𝑣𝑦∗󰇗=−𝛽(12𝑑𝐵0∗ 𝑑𝑧∗sin(𝜙∗)cos(𝜃∗)−𝑦∗ 2 𝑑2𝐵0∗ 𝑑𝑧∗2cos(𝜙∗)cos(𝜃∗)) (8.11) 𝑣𝑧∗󰇗=−1−𝛽(−12𝑥∗ 𝑑2𝐵0∗ 𝑑𝑧∗2sin(𝜃∗) +12𝑦∗ 𝑑2𝐵0∗ 𝑑𝑧∗2sin(𝜙∗)cos(𝜃∗) +𝑑𝐵0∗ 𝑑𝑧∗cos(𝜙∗)cos(𝜃∗)−𝑥∗2+𝑦∗2 4 𝑑3𝐵0∗ 𝑑𝑧∗3cos(𝜙∗)cos(𝜃∗)) (8.12) 8 Nondimensionalization and summary of equations 47 47 𝜔𝑥∗󰇗=𝛽𝑅(−12𝑦∗𝑑𝐵0∗ 𝑑𝑧∗cos(𝜙∗)+𝐵0∗sin(𝜙∗)−𝑥∗2+𝑦∗2 4 𝑑2𝐵0∗ 𝑑𝑧∗2sin(𝜙∗)) −𝜔𝑦∗(𝜎𝜔𝑧∗−𝜔𝑥∗tan(𝜃∗)) (8.13) 𝜔𝑦∗󰇗=−𝛽𝑅(−12𝑥∗𝑑𝐵0∗ 𝑑𝑧∗cos(𝜃∗)−12𝑦∗𝑑𝐵0∗ 𝑑𝑧∗sin(𝜃∗)sin(𝜙∗)−𝐵0∗sin(𝜃∗)cos(𝜙∗) +𝑥∗2+𝑦∗2 4 𝑑2𝐵0∗ 𝑑𝑧∗2sin(𝜃∗)cos(𝜙∗))+𝜔𝑥∗(𝜎𝜔𝑧∗−𝜔𝑥∗tan(𝜃∗)) (8.14) 𝜔𝑧∗󰇗=0 (8.15) 8.2.2. Hamilton’s equations x󰇗∗=𝑝𝑥∗ (8.16) y󰇗∗=𝑝𝑦∗ (8.17) z󰇗∗=𝑝𝑧∗ (8.18) 𝜙󰇗∗=𝑅𝑝𝜙 ∗−𝑝𝜓 ∗sin(𝜃∗) cos2(𝜃∗) (8.19) 𝜃󰇗∗=𝑅𝑝𝜃∗ (8.20) 𝜓󰇗∗=𝑅𝜎𝑝𝜓 ∗−𝑅sin(𝜃∗) cos2(𝜃∗)(𝑝𝜙 ∗−𝑝𝜓 ∗sin(𝜃∗)) (8.21) 𝑝𝑥∗󰇗=−𝛽(−12𝑑𝐵0∗ 𝑑𝑧∗sin(𝜃∗)−𝑥∗ 2 𝑑2𝐵0∗ 𝑑𝑧∗2cos(𝜙∗)cos(𝜃∗)) (8.22) 𝑝𝑦∗󰇗=−𝛽(12𝑑𝐵0∗ 𝑑𝑧∗sin(𝜙∗)cos(𝜃∗)−𝑦∗ 2 𝑑2𝐵0∗ 𝑑𝑧∗2cos(𝜙∗)cos(𝜃∗)) (8.23) 𝑝𝑧∗󰇗=−1−𝛽(−12𝑥∗ 𝑑2𝐵0∗ 𝑑𝑧∗2sin(𝜃∗) +12𝑦∗ 𝑑2𝐵0∗ 𝑑𝑧∗2sin(𝜙∗)cos(𝜃∗) +𝑑𝐵0∗ 𝑑𝑧∗cos(𝜙∗)cos(𝜃∗)−𝑥∗2+𝑦∗2 4 𝑑3𝐵0∗ 𝑑𝑧∗3cos(𝜙∗)cos(𝜃∗)) (8.24) 8 Nondimensionalization and summary of equations 48 48 𝑝𝜙 ∗󰇗=−𝛽(12𝑦∗𝑑𝐵0∗ 𝑑𝑧∗cos(𝜙∗)cos(𝜃∗)−𝐵0∗sin(𝜙∗)cos(𝜃∗) +𝑥∗2+𝑦∗2 4 𝑑2𝐵0∗ 𝑑𝑧∗2sin(𝜙∗)cos(𝜃∗)) (8.25) 𝑝𝜃∗󰇗=𝑅 cos3(𝜃∗)(𝑝𝜓 ∗sin(𝜃∗)−𝑝𝜙 ∗)(𝑝𝜙 ∗sin(𝜃∗)−𝑝𝜓 ∗) −𝛽(−12𝑥∗𝑑𝐵0∗ 𝑑𝑧∗cos(𝜃∗)−12𝑦∗𝑑𝐵0∗ 𝑑𝑧∗sin(𝜙∗)sin(𝜃∗)−𝐵0∗cos(𝜙∗)sin(𝜃∗) +𝑥∗2+𝑦∗2 4 𝑑2𝐵0∗ 𝑑𝑧∗2cos(𝜙∗)sin(𝜃∗)) (8.26) 𝑝𝜓 ∗󰇗=0 (8.27) In both cases, the nondimensional expression of 𝐵0∗ and its derivatives with respect to 𝑧∗are 𝐵0∗=(𝛾2 (𝛾2+𝑧∗2)3/2 −1 (1+𝑧∗2)3/2 ) (8.28) d𝐵0∗ d𝑧∗=3( 𝑧∗ (1+𝑧∗2)52−𝛾2𝑧∗ (𝛾2+𝑧∗2)52) (8.29) d2𝐵0∗ d𝑧∗2=3( 1 (1+𝑧∗2)52−𝛾2 (𝛾2+𝑧∗2)52+5𝛾2𝑧∗2 (𝛾2+𝑧∗2)72−5𝑧∗2 (1+𝑧∗2)72) (8.30) d3𝐵0∗ d𝑧∗3=3( 15𝛾2𝑧∗ (𝛾2+𝑧∗2)72−15𝑧∗ (1+𝑧∗2)72+35𝑧∗3 (1+𝑧∗2)92−35𝛾2𝑧∗3 (𝛾2+𝑧∗2)92) (8.31) From this point forward, the asterisk * is supressed in order to ease the notation. The objective of the next chapters is study the Hamiltonian set of equations, ignoring the trivial equations for ψ and using 𝑝𝜓 as a parameter. Study of the linear system 49 49 9 STUDY OF THE LINEAR SYSTEM 9.1. Position of equilibrium The positions of equilibrium of the spinning top, in the Hamiltoninan formulation, are given by the set 𝛪={𝑥=𝑦=0; 𝑧=𝑧𝑒; 𝜙=𝜃=0; 𝑝𝑥=𝑝𝑦=0; 𝑝𝑧=0; 𝑝𝜙=𝑝𝜃=0} (9.1) We must check that the time derivatives of each of the variables used to define 𝛪 are equal to zero at each point in the set 𝛪. For example, 𝑥󰇗=𝑝𝑥, and since 𝑝𝑥=0 on 𝛪, it follows that 𝑥󰇗=0. The same can be done for the rest of variables. As ψ is a cyclic coordinate, its corresponding conjugate momenta is a constant of the motion that can be treated as a parameter. The value of ψ, if it was necessary, could be obtained later. At this position, the spinning top has its center of mass at rest at the z-axis, and it is spinning with its local axis 𝑍4 fixed along the 𝑂𝑍1-axis. Figure 9-1: Equilibrium points, given by 𝑧=𝑧𝑒 9.2. Linearization about the equilibrium position We now define the local variables 𝑞𝑖=𝑞𝑖−𝑞𝑖𝑒, 𝑝𝑖=𝑝𝑖−𝑝𝑖𝑒 and Taylor expand each of the equations of motion about the equilibrium position given by 𝛪, retaining the linear terms. This leads to the linear system Study of the linear system 50 50 { 𝑥󰇗=𝑝𝑥    𝑚 𝑦󰇗=𝑝𝑦    𝑚 𝑧󰇗=𝑝𝑧    𝑚 𝜙󰇗=𝑝𝜙     𝐼𝑋−𝑝𝜓 𝐼𝑋𝜃 𝜃󰇗=𝑝𝜃  𝐼𝑋 𝑝𝑥    󰇗=𝜇𝐵0′(𝑧𝑒) 2𝜃+𝜇𝐵0′′(𝑧𝑒) 2𝑥 𝑝𝑦    󰇗=−𝜇𝐵0′(𝑧𝑒) 2𝜙+𝜇𝐵0′′(𝑧𝑒) 2𝑦 𝑝𝑧    󰇗=−𝜇𝐵0′′(𝑧𝑒)𝑧 𝑝𝑧    󰇗=−𝜇𝐵0′′(𝑧𝑒)𝑧 𝑝𝜙     󰇗=−𝜇𝐵0′(𝑧𝑒) 2𝑦+𝜇𝐵0(𝑧𝑒)𝜙 𝑝𝜃    󰇗=𝜇𝐵0′(𝑧𝑒) 2𝑥+(𝜇𝐵0(𝑧𝑒)−𝑝𝜓 2 𝐼𝑋)𝜃+𝑝𝜓 𝐼𝑋𝑝𝜙     (9.2) In matrix format, the set of equations can be written as 𝒒󰇗=𝑳 𝒒 (93) where L is the matrix that gathers all the coefficients of the linearization ( 𝑥󰇗𝑦󰇗𝑧󰇗𝜙󰇗𝜃󰇗 𝑝𝑥    󰇗 𝑝𝑦    󰇗𝑝𝑧    󰇗 𝑝𝜙     󰇗 𝑝𝜃    󰇗 ) = ( 0 0 0 0 0 1 𝑚0 0 0 0 0 0 0 0 0 0 1 𝑚0 0 0 0 0 0 0 0 0 0 1 𝑚0 0 0 0 0 0 −𝑝𝜓 𝐼𝑋0001 𝐼𝑋0 0 0 0 0 0 0 0 0 0 1 𝐼𝑋 𝜇𝐵0′′(𝑧𝑒) 20 0 0 𝜇𝐵0′(𝑧𝑒) 20 0 0 0 0 0𝜇𝐵0′′(𝑧𝑒) 20−𝜇𝐵0′(𝑧𝑒) 20 0 0 0 0 0 0 0 −𝜇𝐵0′′(𝑧𝑒)0 0 0 0 0 0 0 0−𝜇𝐵0′(𝑧𝑒) 20 𝜇𝐵0(𝑧𝑒) 0 0 0 0 0 0 𝜇𝐵0′(𝑧𝑒) 20 0 0 (𝜇𝐵0(𝑧𝑒)−𝑝𝜓 2 𝐼𝑋) 0 0 0 𝑝𝜓 𝐼𝑋0 ) ( 𝑥𝑦𝑧𝜙𝜃 𝑝𝑥    𝑝𝑦    𝑝𝑧    𝑝𝜙     𝑝𝜃    ) (9.4) In terms of the nondimensional variables, we have: Study of the linear system 51 51 { 𝑥󰇗=𝑝𝑥    𝑦󰇗=𝑝𝑦    𝑧󰇗=𝑝𝑧    𝜙󰇗=𝑅𝑝𝜙     −𝑅𝑝𝜓𝜃 𝜃󰇗=𝑅𝑝𝜃  𝑝𝑥    󰇗=𝛽𝐵0′(𝑧𝑒) 2𝜃+𝛽𝐵0′′(𝑧𝑒) 2𝑥 𝑝𝑦    󰇗=−𝛽𝐵0′(𝑧𝑒) 2𝜙+𝛽𝐵0′′(𝑧𝑒) 2𝑦 𝑝𝑧    󰇗=−𝛽𝐵0′′(𝑧𝑒)𝑧 𝑝𝜙     󰇗=−𝛽𝐵0′(𝑧𝑒) 2𝑦+𝛽𝐵0(𝑧𝑒)𝜙 𝑝𝜃    󰇗=𝛽𝐵0′(𝑧𝑒) 2𝑥+(𝛽𝐵0(𝑧𝑒)−𝑅𝑝𝜓 2)𝜃+𝑅𝑝𝜓𝑝𝜙     (9.5) Using that 𝛽𝐵0′(𝑧𝑒𝑞)=−1, the system of equations can be simplified to: ( 𝑥󰇗𝑦󰇗𝑧󰇗𝜙󰇗𝜃󰇗 𝑝𝑥    󰇗 𝑝𝑦    󰇗𝑝𝑧    󰇗 𝑝𝜙     󰇗 𝑝𝜃    󰇗 ) = ( 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 −𝑅𝑝𝜓0 0 0 𝑅 0 0 0 0 0 0 0 0 0 0 𝑅 𝛽𝐵0′′(𝑧𝑒) 20 0 0 −120 0 0 0 0 0𝛽𝐵0′′(𝑧𝑒) 20120 0 0 0 0 0 0 0 −𝛽𝐵0′′(𝑧𝑒)0 0 0 0 0 0 0 0120 𝛽𝐵0(𝑧𝑒) 0 0 0 0 0 0 −120 0 0 (𝛽𝐵0(𝑧𝑒)−𝑅𝑝𝜓 2) 0 0 0 𝑅𝑝𝜓0 ) ( 𝑥𝑦𝑧𝜙𝜃 𝑝𝑥    𝑝𝑦    𝑝𝑧    𝑝𝜙     𝑝𝜃    ) (9.6) The solutions of the characteristic polynomial of the matrix 𝑳 (its eigenvalues) provide the ten natural frequencies of the system: |𝑳−𝜆𝑰|=0 (9.7) Moreover, the solutions of the equation 𝑳𝜑 󰇍 𝑖 =𝜆𝑖𝜑 󰇍 𝑖 for each of the natural frequencies 𝜆𝑖 give the ten eigenvectors 𝝋𝒊, which represent the normal modes of the linear system. The analysis of these natural modes will be done at chapter 9.8. From now on, the constants 𝑆 and 𝐵 are introduced to simplify the expressions: 𝑆=𝛽d2𝐵0 d𝑧2(𝑧𝑒)=−𝐵0′′ 𝐵0′(𝑧𝑒) 𝐵=𝛽𝐵0(𝑧𝑒)=−𝐵0 𝐵0′(𝑧𝑒) (9.8) As all the coefficients of the system of equations are real, given a real initial condition the solution of the set will be real at any time. From this point forward, linear coordinates 𝑞𝑖 and conjugate momenta 𝑝𝑖 are going to be written without the upper line in  in order to ease notation. Study of the linear system 52 52 9.3. Solution of the system: expression with real trigonometric functions We can try to find a basis of real solutions which enables to express any other as a combination of them. As we can see from equations (9.5) , vertical motion is uncoupled from the horizontal one, which involves the coordinates x, y and its respective conjugate momenta 𝑝𝑥, 𝑝𝑦; and the whirl one, given by the degrees of freedom 𝜙, 𝜃 and its conjugate momenta 𝑝𝜙, 𝑝𝜃 . Therefore, on one hand we have 𝑧󰇗=𝑝𝑧 𝑝󰇗𝑧=−𝑆𝑧 (9.9) that results in a simple harmonic motion in vertical direction with natural frequency 𝑆=𝛽d2𝐵0 d𝑧2(𝑧𝑒) (9.10) and on the other hand we have the rest of equations of the system. As we are looking for periodic orbits, we start with a solution with the following form: 𝑥(𝑡)=𝑥0cos(𝜆𝑡) 𝑦(𝑡)=𝑦0sin(𝜆𝑡) (in principle, 𝑦(𝑡) might also have a cosine component, but we will see that it does not appear). Here, if λ and 𝑦 have the same sign, the orbit is traversed in the positive sense (counterclockwise) around the 𝑂𝑍1axis, whereas opposite sign shows that the trajectory goes in the negative sense (clockwise). Differentiating once, we obtain 𝑝𝑥=−𝜆𝑥0sin(𝜆𝑡) 𝑝𝑦=𝜆𝑦0cos(𝜆𝑡) and differentiating again and solving for 𝜃(𝑡), 𝜙(𝑡): 𝜃(𝑡)=(𝑆+2𝜆2)𝑥0cos(𝜆𝑡) 𝜙(𝑡)=−(𝑆+2𝜆2)𝑦0sin(𝜆𝑡) (9.11) Note that, for each mode, the vector 𝜃𝑖 1−𝜙𝑗 1 is proportional to 𝑥𝑖 1+𝑦𝑗 1. This implies that the dipole precesses while it is spinning, moving on the surface of an imaginary cone. Furthermore, as the proportionality constant is positive (𝜆2 is positive and S is also positive in the region of interest) the unit vector 𝑘 󰇍 3 always points outwards. Conjugates momenta 𝑝𝜃(𝑡) and 𝑝𝜙(𝑡) are obtained differentiating the expressions of 𝜃(𝑡) and 𝜙(𝑡) and solving: 𝑝𝜃(𝑡)=−𝜆𝑅(𝑆+2𝜆2)𝑥0sin(𝜆𝑡) (9.12) 𝑝𝜙(𝑡)=(𝑝𝜓𝑥0−𝜆𝑅𝑦0)(𝑆+2𝜆2)cos(𝜆𝑡) (9.13) Study of the linear system 53 53 and leading these results to the last two equations of the set, we have: −𝜆(𝑝𝜓𝑥0−𝜆𝑅𝑦0)(𝑆+2𝜆2)=𝑦0(12−𝐵(𝑆+2𝜆2)) (9.14) −𝜆2 𝑅(𝑆+2𝜆2)𝑥0=𝑥0((𝐵−𝑅𝑝𝜓2)(𝑆+2𝜆2)−12)+𝑅𝑝𝜓(𝑝𝜓𝑥0−𝜆𝑅𝑦0)(𝑆+2𝜆2) (9.15) Simplifying and rewriting (9.15) in order to have a similar expression as (9.14), −𝜆(𝑝𝜓𝑥0−𝜆𝑅𝑦0)(𝑆+2𝜆2)=𝑦0(12−𝐵(𝑆+2𝜆2)) (9.16) 𝜆(𝑝𝜓𝑦0−𝜆𝑅𝑥0)(𝑆+2𝜆2)=𝑥0(𝐵(𝑆+2𝜆2)−12) (9.17) and dividing both expressions, we finally get 𝑝𝜓𝑦0−𝜆𝑅𝑥0 𝑝𝜓𝑥0−𝜆𝑅𝑦0=𝑥0 𝑦0→𝑝𝜓𝑦02−𝜆𝑅𝑥0𝑦0=𝑝𝜓𝑥02−𝜆𝑅𝑥0𝑦0 𝑦0=±𝑥0 (9.18) We can see that oscillations along x and y have the same amplitude, and therefore the normal modes correspond with circular orbits. The same result can be obtained, as it could not be otherwise, considering the linear set of equations derived from the equations of vector mechanics. For the vertical motion, we have 𝑧󰇗=𝑣𝑧 𝑣󰇗𝑧=−𝑆𝑧 (9.19) and the horizontal and whirl motions are described by the following equations { 𝑥󰇗=𝑣𝑥 𝑦󰇗=𝑣𝑦 𝑣󰇗𝑥=−12𝜃+𝑆2𝑥 𝑣󰇗𝑦=12𝜙+𝑆2𝑦 𝜙󰇗=𝜔𝑋 𝜃󰇗=𝜔𝑌 𝜔󰇗𝑋=−𝜎ω𝜔𝑌+𝑅2𝑦+𝑅𝐵𝜙 𝜔󰇗𝑌=𝜎ω𝜔𝑋−𝑅2𝑥+𝑅𝐵𝜃 (9.20) Study of the linear system 54 54 By simple inspection, we can see that 𝑣𝑥=𝑝𝑥, 𝑣𝑦=𝑝𝑦 and 𝑣𝑧=𝑝𝑧, so 𝑣𝑥=−𝜆𝑥0sin(𝜆𝑡) 𝑣𝑦=𝜆𝑦0cos(𝜆𝑡) (9.21) Nevertheless, one can see that the equations for the whirl motion in the vector mechanics set are simpler than in the Hamiltonian formulation , and the expressions of 𝜔𝑋 and ωY are different to 𝑝𝜙 and 𝑝𝜃: ωY=−𝜆(2𝜆2+𝑆)𝑥0sin(𝜆𝑡) 𝜔𝑋=−𝜆(2𝜆2+𝑆)𝑦0cos(𝜆𝑡) (9.22) Leading these results to the last two equations, ((𝜆2+𝑅𝐵)(2𝜆2+𝑆)−𝑅2)𝑦0=𝜎ω𝜆(2𝜆2+𝑆)𝑥0 ((𝜆2+𝑅𝐵)(2𝜆2+𝑆)−𝑅2)𝑥0=𝜎ω𝜆(2𝜆2+𝑆)𝑦0 (9.23) And dividing both expressions, we get again the same relation between the amplitudes 𝑦0 𝑥0=𝑥0 𝑦0 ⇒ 𝑦0=±𝑥0 (9.24) One can see that the yardstick to determine whether the orbit is counterclockwise or clockwise is given by the factor Δ=((𝜆2+𝑅𝐵)(2𝜆2+𝑆)−𝑅2) (9.25) If Δ>0, 𝑦0 and 𝜆 have the same sign and the orbit is counterclockwise, whereas a negative sign means clockwise sense. Note that the sign of 𝜆 also plays a crucial role. We are assuming that 𝜆 is positive, so that when we refer to positive spin we actually mean that “the top spins in the same sense to the proper rotation” and negative means that “it spins in the opposite sense to the proper rotation”. If we took the proper rotation in the opposite sense, signs of −𝜆 would change. Considering 𝑦0=+𝑥0, we have the characteristic equation (𝜆2+𝑅𝐵)(2𝜆2+𝑆)−𝑅2=𝜎ω𝜆(2𝜆2+𝑆) (9.26) that provides the four eigenvalues. To complete the basis of solutions, for each solution of the form 𝑥=𝑥0cos(𝜆𝑡) 𝑦=𝑥0sin(𝜆𝑡) there is another Study of the linear system 55 55 𝑥=−𝑥0sin(𝜆𝑡) 𝑦=𝑥0cos(𝜆𝑡) that fulfills the same equations. The complete solution is a linear combination of the 8 functions 𝑥=∑(𝑎𝑖cos(𝜆𝑖𝑡)−𝑏𝑖sin(𝜆𝑖𝑡)) 4 𝑖=1 (9.27) and the corresponding for the rest of variables. 9.4. Expression with complex exponentials The previous results can also be written in term of complex exponential functions. Using Euler’s formula cos(𝜆𝑡)=12ej𝜆𝑡+12e−j𝜆𝑡 sin(𝜆𝑡)=−j2ej𝜆𝑡+j2e−j𝜆𝑡 (9.28) the expression for x is 𝑥=∑(𝑎𝑖+j𝑏𝑖 2ej𝜆𝑖𝑡+𝑎𝑖−j𝑏𝑖 2e−j𝜆𝑖𝑡) 4 𝑖=1 (9.29) Introducing 𝑐𝑖 𝑐𝑖=𝑎𝑖+j𝑏𝑖 2 𝑐𝑖∗=𝑎𝑖−j𝑏𝑖 2 we have 𝑥=∑(𝑐𝑖ej𝜆𝑖𝑡+𝑐𝑖∗e−j𝜆𝑖𝑡) 4 𝑖=1 (9.30) The same can be done for the rest of variables: 𝑦=∑(−j𝑐𝑖ej𝜆𝑖𝑡+j𝑐𝑖∗e−j𝜆𝑖𝑡) 4 𝑖=1 (9.31) 𝜃=(2𝜆2+𝑆)∑(𝑐𝑖ej𝜆𝑖𝑡+𝑐𝑖∗e−j𝜆𝑖𝑡) 4 𝑖=1 (9.32) etc. If we consider the solutions for λ and −𝜆 as independent modes, the following eigenvector corresponds Study of the linear system 62 62 9.7. Stability region In eq. (9.54) the values of the parameters 𝜎 and 𝑅 were given in chapter 8.1. Moreover, as 𝛽 is related with the height of levitation 𝑧𝑒 by means of eq. (5.15), we can look for the range 𝜔𝑚𝑖𝑛≤𝜔≤𝜔𝑚𝑎𝑥 that provides a stable solution (four real roots 𝜆𝑖) for each 𝑧𝑒. Also, in chapter 5.3 we determined the minimum value of 𝑧𝑒 for which levitation is possible (remember that below 𝑧𝑒=2.9235 the spinning top cannot hover), so the region starts in this value. The physical reason for the existence of the lower limit 𝜔𝑚𝑖𝑛 is obvious. Since the effect on stability due to the gyroscopic effect is proportional to the spin velocity, if the top is spun too slowly then it becomes unstable against rotation and tips over. On the other hand, if the top is spun too fast, its axis becomes too stiff and cannot respond quickly enough to the changing direction of the magnetic field. Then the magnetic dipole moment, which has the same direction of the axis of the top, can be considered as fixed in space and, according to Earnshaw’s theorem, becomes unstable against translations. This explains the existence of the upper limit 𝜔𝑚𝑎𝑥. The stability region has the following form In terms of the 𝛽 parameter, linear stability is possible within the range 30.43≤𝛽≤31.68 This corresponds with a stable range of 𝑧𝑒 Figure 9-3: Stability region. Study of the linear system 63 63 2.9235≤𝑧𝑒≤3.22 or, in dimensional variables 7.31 𝑐𝑚≤𝑧𝑒≤8.05 𝑐𝑚 Genta et al. [4] provide an analytical approximation of the expressions for the lower and upper stability limits. The lower bound takes its origin from the coupling between the two rotational modes. To analytically compute the value, it is sufficient to study the stability of the rotational motion (described by the 𝜙 and 𝜃 degrees of freedom), considered as uncoupled with the vibrational one. Using the complex variable 𝜃, the equation of the whirl motion becomes 𝐼𝑋𝜃󰇘−𝑗𝜔𝐼𝑍𝜃󰇗−𝜇𝐵0𝜃=0 (9.55) A solution 𝜃=𝜃0𝑒𝑗𝜆𝑡 provides 𝐼𝑋𝜆2−𝐼𝑍𝜔𝜆+𝜇𝐵0(𝑧𝑒)=0→𝜆=𝐼𝑍𝜔±√(𝐼𝑍𝜔)2−4𝐼𝑋𝜇𝐵0(𝑧𝑒) 2𝐼𝑋 (9.56) For the roots 𝜆𝑖 to be real, the discriminant of the quadratic equation must be positive. This condition provides the lower limit of 𝜔 to have a stable solution: 𝜔𝑚𝑖𝑛=2√𝜇𝐵0(𝑧𝑒)𝐼𝑋 𝐼𝑍2 (9.57) To obtain the upper limit, one can use the fast top condition. Doing 𝐼𝑋=0, equations become: {𝑚𝑟󰇘−𝜇𝐵0′′(𝑧𝑒) 2𝑟−𝜇𝐵0′(𝑧𝑒) 2 𝜃=0 𝐼𝑋𝜃󰇘−𝐼𝑍𝜔𝑗𝜃󰇗−𝜇𝐵0(𝑧𝑒)𝜃−𝜇𝐵0′(𝑧𝑒) 2𝑟=0 (9.58) The characteristic equation of the previous set is the following cubic function: 4𝑚𝐼𝑧𝜔𝜆3−4𝑚𝜇𝐵0(𝑧𝑒)𝜆2+2𝜇𝐵0′′(𝑧𝑒)𝐼𝑧𝜔𝜆+𝑚2𝑔2−2𝜇2𝐵0′′(𝑧𝑒)𝐵0(𝑧𝑒)=0 (9.59) Using that the discriminant of a cubic equation with the form 𝑎3𝑧3+𝑎2𝑧2+𝑎1𝑧+𝑎0=0 is given by 𝐷3=𝑎12𝑎22−4𝑎0𝑎23−4𝑎13𝑎3+18𝑎0𝑎1𝑎2𝑎3−27𝑎02𝑎32 and introducing the variables 𝜁= 𝐵0′(𝑧𝑒)2 2𝐵0(𝑧𝑒)𝐵0′′(𝑧𝑒)−1 𝑓(𝜁)=1−18𝜁−27𝜁2 Study of the linear system 64 64 𝛾𝑚𝑎𝑥=√1+𝜁 2[𝑓(𝜁)+√𝑓(𝜁)2+64𝜁] One gets the expression for the maximum rotational speed is 𝜔𝑚𝑎𝑥=𝛾𝑚𝑎𝑥 𝑔𝐼𝑍√𝜇3𝐵0(𝑧𝑒)3 𝑚 (9.60) Based on figure 9-3 and the definition of 𝛽 in (5.11), we can define the dimensionless variable 𝑚0 𝑚0=𝜇𝐵0 𝑔𝑎 and plot the nondimensional variable 𝑚 𝑚0 against the rotational speed. Figure 9-4: Stability region in terms of the mass and rotational speed. Therefore, figure 9-4 provides the mass tolerance of the toy Δ𝑚 𝑚=100𝑚𝑚𝑎𝑥−𝑚𝑚𝑖𝑛 𝑚𝑚𝑎𝑥 ≃3.8 % Following Gov et al. [6], experimentally the tolerance is only about 1%. Note, however, that the lower mass region is difficult to access. Study of the linear system 65 65 Two numerical examples are given to clarify all the ideas presented: Example 9.7.1 Evolution of the natural frequencies for a centered value of 𝑧𝑒 in the stability region: 𝑧𝑒=3.1 The centered value 𝑧𝑒=3.1 (that equals to a dimensional height of 7.75 cm) corresponds with 𝛽=30.9078 For this 𝑧𝑒, the range of 𝜔 that provides linear stability is given by 𝜔𝑚𝑖𝑛=6.13≡1159 𝑟𝑝𝑚 𝜔𝑚𝑎𝑥=26.67≡5045 𝑟𝑝𝑚 The evolution of the four natural frequencies of the system for that height of levitation as a function of ω (with ω moving in the previous range) is depicted in the following Campbell diagram: Figure 9-5: Campbell diagram for 𝑧𝑒=3.1 We can compute these natural frequencies for a particular value of 𝜔 within the stable range. Therefore, for 𝜔=20 (for which we are situated in the middle of the stability region) Study of the linear system 66 66 Figure 9-6: Centered point in the stability region. the natural frequencies 𝜆𝑖 are Figure 9-7: Natural frequencies for a centered value in the stability region. 𝜆1=−0.1122 𝜆2=0.2481 𝜆3=0.8294 𝜆4=39.0347 Note that the natural frequency 𝜆1 corresponds with the slowest and retrograde mode (because of the minus sign). Study of the linear system 67 67 Example 9.7.2 Evolution of the natural frequencies for a noncentered value of 𝑧𝑒 in the stability region: 𝑧𝑒=3.19 The noncentered value 𝑧𝑒=3.19 (that equals to a dimensional height of 7.98 cm) corresponds with 𝛽=31.4581 For this 𝑧𝑒, the range of 𝜔 that provides linear stability is given by 𝜔𝑚𝑖𝑛=20.80≡3935 𝑟𝑝𝑚 𝜔𝑚𝑎𝑥=28.78≡5444 𝑟𝑝𝑚 The evolution of the four natural frequencies of the system for that height of levitation as a function of ω (with ω moving in the previous range) is given in the following Campbell diagram: Figure 9-8: Diagram Campbell for 𝑧𝑒=3.19. As 𝜆4 (the natural frequency corresponding with the fast mode) is much larger than the others, we show the evolution of 𝜆1, 𝜆2, 𝜆3 separately from the previous one. Study of the linear system 68 68 Figure 9-9: Evolution of 𝜆1, 𝜆2 and 𝜆3 for 𝑧𝑒=3.19 Computing these natural frequencies for 𝜔=24.79 (in the middle of 𝜔𝑚𝑖𝑛 and 𝜔𝑚𝑎𝑥), we are situated in Figure 9-10: Noncentered point in the stability region. and the natural frequencies are Study of the linear system 69 69 Figure 9-11: Natural frequencies for a noncentered point in the stability region. Figure 9-12: Natural frequencies for a noncentered point in the stability region. 𝜆1=0.0444 𝜆2=0.1242 𝜆3=0.6481 𝜆4=45.1834 Study of the linear system 70 70 9.8. Analysis of the normal modes As we said in 9.2, solving the characteristic equation |𝑳−𝜆𝑰|=0, where L is the 10x10 linear matrix in the stability region previously bounded, we obtain 5 pairs of conjugate imaginary frequencies. Moreover, each of the frequencies 𝜆𝑖 have its corresponding eigenvector 𝜑 󰇍 𝑖, and these eigenvectors correspond with the normal modes. Nevertheless, as we have complex eigenvectors, an immediate identification with real motions is not possible. Taking into account the linearity of the equations, we can combine the corresponding +𝜆 and −𝜆 solutions so that we still get oscillating solutions with the same frequencies. In practice, this is the same as analyzing real and imaginary parts separately and study their associated motions. These ideas are clarified in the following example: Example 9.8.2 Identification of the motions given two complex conjugate eigenvectors We can numerically compute the eigenvalues and eigenvectors of matrix L. To continue with example 9.7.1, that corresponds with a centered value of the stability region, we choose 𝑧𝑒=3.1≡7.75 cm 𝜔=20≡3783 rpm and particularize matrix L for this values. The identification of the motions is going to be explained considering two of the ten complex cojugate eigenvectors. Naming them 𝜑1and 𝜑2, we have 𝜑1= ( 0 −0.0001 0 0.3545 0 0.0045 0 0 0 0.8649 ) +𝑗 ( −0.0001 000 −0.3545 0 −0.0045 0 −0.0214 0 ) 𝜑2= ( 0 −0.0001 0 0.3545 0 0.0045 000 0.8649 ) +𝑗 ( 0.0001 000 0.3545 0 0.0045 0 0.0214 0 ) As it was said before, we take the real and imaginary parts separately of one of the complex eigenvectors and study their associate motions. Therefore, taking into accout that the components of the real or imaginary parts correspond with 𝜑 󰇍 𝑖= ( 𝑥0 𝑦0 𝑧0 𝜙0 𝜃0 𝑣𝑥0 𝑣𝑦0 𝑣𝑧0 𝜔𝑋0 𝜔𝑌0 ) we consider 𝜑1 (it would be the same for 𝜑2). Its real part provides Study of the linear system 71 71 𝑦0=−0.0001 and 𝑣𝑥0=0.0045 so we start in the negative part of the Y-axis (as we have 𝑦0<0) with a positive velocity in the X-direction (as 𝑣𝑥0>0). Therefore, the trajectory corresponds with a counterclockwise circumference (see figure 9-13, left). Similarly, the imaginary part provides 𝑥0=−0.0001 and 𝑣𝑦0=−0.0045 and we start in the negative part of the X-axis with a negative velocity in the Y-direction (figure 9-13, right). Note that, as in the case mentioned above, we have a counterclockwise circumference and the interpretation of the real and imaginary parts results in the same trajectory. Figure 9-13: Identification of the motions using real and imaginary parts Therefore, the five independent normal modes of the Hamiltonian system for 𝑧𝑒=3.1≡7.75 cm 𝜔=20≡3783 rpm scaled to 𝜃0=1°=𝜋 180 rad, are 𝝋 󰇍 󰇍 𝟏= ( 0.0920 000 0.0175 0 −0.0103 0 0.0438 0 ) 𝝋 󰇍 󰇍 𝟐= ( 0.0607 000 0.0175 0 0.0151 0 0.0434 0 ) 𝝋 󰇍 󰇍 𝟑= ( 0.0113 000 0.0175 0 0.0094 0 0.0427 0 ) 𝝋 󰇍 󰇍 𝟒= ( 5.7· 10−6 000 0.0175 0 2.2· 10−4 0 0.0011 0 ) 𝝋 󰇍 󰇍 𝟓= ( 00 𝑧0 0000000 ) where we have continued with the color pattern used to compute the evolution of the natural frequencies in examples 9.7.1 and 9.7.2. Thus, 𝜑 󰇍 1 corresponds with the slow, retrograde mode and 𝜑 󰇍 4 is the fast mode. Also note that 𝜑 󰇍 5 only has as nonzero component 𝑧0, as it corresponds with the vertical motion (remember that it is uncoupled with the horizontal and whirl motions). Arrived to this point, we numerically solve the linear Hamiltonian system introducing the previous eigenvectors as initial conditions. Following eq. (9.18), the trajectories correspond with four circumferences: Study of the linear system 78 78 In figure 9-26 we can see the trajectories of the spinning top for three different values of rotational speed within the range of 9.7.1: Figure 9-25: Trajectories of the linear system for horizontal perturbations 9.9.2. Trajectory of the top considering a pitch perturbation We next consider a pitch perturbation. This corresponds with the initial condition 𝒒𝟎 and figure 9-27. 𝑞0= ( 𝑥0 𝑦0 𝑧0 𝜙0 𝜃0 𝑝𝑥0 𝑝𝑦0 𝑝𝑧0 𝑝𝜙0 𝑝𝜃0 ) = ( 00 3.1 0 0.035 0 0 0 0 0 ) Figure 9-26: Pitch perturbation of the top Study of the linear system 79 79 The trajectory is shown in figure 9-28. We note that, despite the top starts in the 𝑂𝑍1 axis, the perturbation in 𝜃0 causes small displacements, moving the top away from the origin (remember that, in the linear system, translational and whirl motions continue being coupled). Figure 9-27: Trajectory of the linear system with pitch perturbation It should be mentioned that the representation of the horizontal and pitch perturbations shown in figures 9-25 and 9-27 has been exagerated in order to provide a better understanding. In the real toy, these two situations can be matched with the influence of weak air streams or a bad use of the starter, placing it off-centered or not vertically in the lifter plate. 9.9.3. Validation of the rotational speed range For 𝑧𝑒=3.1, we saw that 𝜔𝑚𝑖𝑛=6.13≡1159 𝑟𝑝𝑚 𝜔𝑚𝑎𝑥=26.67≡5045 𝑟𝑝𝑚 are the lower and upper rotational speed bounds in order to have stability. Considering a more general initial condition 𝒒𝟎, adding vertical, horizontal and rotational perturbations and also nonzero initial values of the conjugate momenta, we are going to check this range, showing the shape of the trajectories for different values of 𝜔. 𝒒𝟎= ( 𝑥0 𝑦0 𝑧0 𝜙0 𝜃0 𝑝𝑥0 𝑝𝑦0 𝑝𝑧0 𝑝𝜙0 𝑝𝜃0 ) = ( 0.02 03 0.01 0.02 −0.01 0.04 0.01 0.03 −0.01 ) In figure 9-29, the trajectories for a rotational speed under 𝜔𝑚𝑖𝑛 and for this limit are shown. In the first situation, stable levitation is not possible and the top moves away from the 𝑂𝑍1 axis, whereas for 𝜔=𝜔𝑚𝑖𝑛 we see that hovering becomes stable. Study of the linear system 80 80 Figure 9-28: Trajectories for 𝜔<𝜔𝑚𝑖𝑛 and 𝜔=𝜔𝑚𝑖𝑛 In figures 9-30 and 9-31, we show the stable trajectories for values of 𝜔 within the stability range: Figure 9-29: Trajectories within the stable rotational speed range Study of the linear system 81 81 Figure 9-30: Trajectories within the stable rotational speed range Finally, in figure 9-32 we arrive to the upper bound 𝜔𝑚𝑎𝑥. In the first situation, which corresponds with this limit, we still have a stable motion, but for larger rotational speed we see that hovering is not possible. Figure 9-31: Trajectories for 𝜔=𝜔𝑚𝑎𝑥 and 𝜔>𝜔𝑚𝑎𝑥 Study of the linear system 82 82 10 Study of the nonlinear system 83 83 10 STUDY OF THE NONLINEAR SYSTEM At chapter 9.8, we gave the expression of the five independent normal modes. The analytical and numerical study showed that these modes correspond with circular orbits. In this chapter, we numerically solve the nonlinear system with initial conditions proportional to these eigenvectors and study the structural stability of the system. In this way, we are going to start with the study of the linear coupling, continue with the nonlinear coupling and finally see the behavior of the system in interesting situations. 10.1. Study of the linear coupling To study the linear coupling, we numerically solve the nonlinear system choosing as initial condition the expression for the normal modes given in chapter 9.8 and see if the expected circular orbits take place. In figure 10-1 we can see that these circular orbits are well reproduced. In particular, the smaller the circumference is, the better the circumference is replayed. Figure 10-1: Circumferences in the nonlinear system. As happened in 9.8, the circumference corresponding with the fast mode cannot be seen with the others, so we plot it separetly. 10 Study of the nonlinear system 84 84 We also show the time evolution of the 𝑥 and 𝑦 degrees of freedom for the different modes. In figure 10-2 we can see the oscillations of these modes, obtaining similar results as those of figure 9-22. Figure 10-2: Linear coupling: time evoluton of 𝑥, 𝑦. 10.2. Study of the nonlinear coupling After the study of the linear coupling, we can verify if the qualitative behavior of the trajectories is unaffected by small perturbations, in order to verify its structural stability. Therefore, we are going to raise the amplitudes of the previous initial conditions and see when the circular orbits stop being reproduced. As we scaled the normal modes making 𝜃0=1°, this is equivalent to find the maximum pitch angle 𝜀𝑖𝜃0=𝜀𝑖° for which we still have circumferences. For these values 𝜀𝑖, the linear approximation sin(𝜃)≃𝜃 is not true, as the other terms of the Taylor expansion become important. We therefore obtain, for the different modes 𝜀1=2.6° 𝜀2=1.9° 𝜀3=9.4° 𝜀4=5° Figure 10-3: Nonlinear coupling 10 Study of the nonlinear system 85 85 10.3. Behaviour of the system in interesting situations 10.3.1. Base not levelled with gravity When one plays with the device, in all likelihood the spinning top will fly off to a side at some time. If this happens multiple times, it means that the base is not leveled with gravity. Thus, one of the things that we can study is the maximum percentaje of misalignment between the base and gravity that allows stable hovering. It can be done adding a constant 𝜀 in the equations of motion for eather x or y degrees of freedom and testing different values in order to find the maximum 𝑝𝑥∗󰇗=−𝛽(−12𝑑𝐵0∗ 𝑑𝑧sin(𝜃∗)−𝑥∗ 2 𝑑2𝐵0∗ 𝑑𝑧2cos(𝜙∗)cos(𝜃∗))+𝜀 (10.1) 𝑝𝑦∗󰇗=−𝛽(12𝑑𝐵0∗ 𝑑𝑧sin(𝜙∗)cos(𝜃∗)−𝑦∗ 2 𝑑2𝐵0∗ 𝑑𝑧2cos(𝜙∗)cos(𝜃∗)) (10.2) The numerical simulations provide a maximum 𝜀=0.0018. This corresponds with a percentage of misalignment of 0.18 % , which is equivalent to an inclination along this axis of 0.12 °. Therefore, we have numerically verified the great importance of having a levelled base, as it is usually the first problem that one faces with when plays with the toy. 10.3.2. Trajectory of the top considering a horizontal perturbation In the study of the linear system, we obtained the trajectories simulating the influence of weak air streams or a bad use of the starter, placing it off-centered in the lifter plate, by means of a horizontal perturbatioin. The same can be done with the nonlinear system. The initial condition is given by 𝒒𝟎 and figure 10-4 shows the initial position of the top. 𝒒𝟎= ( 𝑥0 𝑦0 𝑧0 𝜙0 𝜃0 𝑝𝑥0 𝑝𝑦0 𝑝𝑧0 𝑝𝜙0 𝑝𝜃0 ) = ( 0 0.02 3.1 000 0 0 0 0 ) Figure 10-4: Horizontal perturbation. 10 Study of the nonlinear system 86 86 The results, for two different rotational speeds, are shown in figure 10-5. The trajectories are similar to those obtained for the linear system (figure 9-26). Figure 10-5: Trajectories of the nonlinear system with horizontal perturbation. 10.3.3. Trajectory of the top considering a pitch perturbation We can also see how the nonlinear system responds to a pitch perturbation. Considering the initial condition given by 𝒒𝟎, the trajectory can be seen in figure 10-7. 𝑞0= ( 𝑥0 𝑦0 𝑧0 𝜙0 𝜃0 𝑝𝑥0 𝑝𝑦0 𝑝𝑧0 𝑝𝜙0 𝑝𝜃0 ) = ( 00 3.1 0 0.035 0 0 0 0 0 ) Figure 10-7: Trajectory of the nonlinear system with pitch perturbation Figure 10-6: Pitch perturbation 10 Study of the nonlinear system 87 87 As happened in 9.9.2, we can see that 𝜃0 causes small displacements, moving the top away from the origin. 10.3.4. Influence of the initial height in the trajectory Another aspect that can be analyzed is how the trajectories change depending on the initial condition 𝑧0. In the nonlinear model, vertical motion is coupled with horizontal and whirl motion, so that different values of 𝑧0 provide different trajectories. The change of 𝑧0 corresponds in the real toy with placing the lifter plate at different heights. We consider again an equilibrium height 𝑧𝑒=3≡7.75 𝑐𝑚. In figure 10-8, we can see that values of 𝑧0=3 and 𝑧0=3.02 provide a flower-shape trajectory in the 𝑋1𝑌1 plane. Figure 10-8: Flower-shape trajectories. Nevertheless, as 𝑧0 tends to 𝑧𝑒, one can see that the three-dimensional trajectories and their projections on the 𝑋1𝑌1 change, losing the flower shape and approaching to the characteristic form of the linear trajectories (see figures 9-31 and 9-31). These changes can be seen in figures 10-9 and 10-10. 10 Study of the nonlinear system 94 94 11 Integrals of the motion 95 95 11 INTEGRALS OF THE MOTION In this chapter we show the constants of motion (quantities that are conserved throughout the trajectory of the dipole) and demonstrate their constancy. 11.1. Conservation of mechanical energy E The constancy of the mechanical energy of the system E is justified by the fact that the gravitational force is a conservative force and the magnetic force does not work. Therefore, E is a constant of motion. Moreover, the Levitron toy is an example of system submitted to ideal links, independent coordinates, conservative active forces and the Lagrangian is explicitly time-independent and quadratic velocity dependant. Therefore, the Hamiltonian is a conserved quantity and it coincides with the total energy of the system: ℋ=𝐸=𝑇+𝑈 (11.1) This constancy can be numerically verified for the example considered in previous chapters. Figure 11-1: Mechanical energy of the system. 11.2. Conservation of 𝑳𝟎𝒁𝟏 We next demonstrate that the Z-component of the angular momentum with respect to the origin of coordinates 𝐿0𝑍1is a constant of the motion, provided that we have cylindrical symmetry. The latter condition implies that 11 Integrals of the motion 96 96 𝐵𝑥 𝑥=𝐵𝑦 𝑦 (11.2) We start from the angular momentum theorem d𝐿 󰇍 0 d𝑡|1=𝑀 󰇍 󰇍 𝐺+𝑟 ×𝐹 (11.3) Using Newton’s Second law 𝐹 =𝑚𝑔 +∇(𝜇 ·𝐵 󰇍 )=𝑚𝑎 󰇍 𝐺 (11.4) Expressing 𝜇 in frame ‘1’ 𝜇 =−𝜇𝑘 󰇍 3=−𝜇( sin(𝜃) −sin𝜙cos𝜃 cos𝜙cos𝜃) (11.5) and computing the cross products 𝜇 ×𝐵 󰇍 and 𝑟 ×𝑚𝑎 󰇍 𝐺 in the fixed frame, we have 𝜇 ×𝐵 󰇍 |1=−𝜇| 𝑖 1𝑗 1𝑘 󰇍 1 sin(𝜃)−sin𝜙cos𝜃 cos𝜙cos𝜃 𝐵𝑥𝐵𝑦𝐵𝑧|=(𝜇(𝐵𝑦cos(𝜙)cos(𝜃)+𝐵𝑧sin(𝜙)cos(𝜃)) 𝜇(𝐵𝑧sin(𝜃)−𝐵𝑥cos(𝜙)cos(𝜃)) −𝜇(𝐵𝑦sin(𝜃)+𝐵𝑥sin(𝜙)cos(𝜃)) ) (11.6) 𝑟 ×(𝑚𝑔 +∇(𝜇 ·𝐵 󰇍 ))|1=𝑟 ×𝑚𝑎 󰇍 𝐺|1=𝑚|𝑖 1𝑗 1𝑘 󰇍 1 𝑥 𝑦 z 𝑥󰇘 𝑦󰇘 𝑧󰇘|=𝑚(𝑦𝑧󰇘−𝑦󰇘𝑧 𝑥󰇘𝑧−𝑥𝑧󰇘 𝑥𝑦󰇘−𝑥󰇘𝑦) (11.7) Substituting in the expression of the angular momentum theorem d𝐿 󰇍 0 d𝑡|1=(𝜇(𝐵𝑦cos(𝜙)cos(𝜃)+𝐵𝑧sin(𝜙)cos(𝜃))+𝑚(𝑦𝑧󰇘−𝑦󰇘𝑧) 𝜇(𝐵𝑧sin(𝜃)−𝐵𝑥cos(𝜙)cos(𝜃))+𝑚(𝑥󰇘𝑧−𝑥𝑧󰇘) −𝜇(𝐵𝑦sin(𝜃)+𝐵𝑥sin(𝜙)cos(𝜃))+𝑚(𝑥𝑦󰇘−𝑥󰇘𝑦)) (11.8) To show the constancy of 𝐿0𝑍1, we substitute the following relations 𝑚𝑥󰇘=−𝜇(𝜕𝐵𝑥 𝜕𝑥sin(𝜃)+𝜕𝐵𝑧 𝜕𝑥cos(𝜙)cos(𝜃)) 𝑚𝑦󰇘=−𝜇(−𝜕𝐵𝑦 𝜕𝑦sin(𝜙)cos(𝜃)+𝜕𝐵𝑧 𝜕𝑦cos(𝜙)cos(𝜃)) (11.9) 11 Integrals of the motion 97 97 in the third component of d𝐿 󰇍 0 d𝑡|1 d𝐿0𝑍1 d𝑡 =−𝜇(𝐵𝑦sin(𝜃)+𝐵𝑥sin(𝜙)cos(𝜃))+𝜇𝑥𝜕𝐵𝑦 𝜕𝑦sin(𝜙)cos(𝜃)−𝜇𝑥𝜕𝐵𝑧 𝜕𝑦cos(𝜙)cos(𝜃) +𝜇𝑦𝜕𝐵𝑥 𝜕𝑥sin(𝜃)+𝜇𝑦𝜕𝐵𝑧 𝜕𝑥cos(𝜙)cos(𝜃) (11.10) The components of the magnetic field, differentiating eq. with respect x and y, verify: 𝑦𝜕𝐵𝑥 𝜕𝑥=𝐵𝑦+𝑥𝜕𝐵𝑦 𝜕𝑥 𝑥𝜕𝐵𝑦 𝜕𝑦=𝐵𝑥+𝑦𝜕𝐵𝑥 𝜕𝑦 Therefore d𝐿0𝑍1 d𝑡 =−𝜇(𝐵𝑦sin(𝜃)+𝐵𝑥sin(𝜙)cos(𝜃))+𝜇𝐵𝑥sin(𝜙)cos(𝜃)+𝜇𝑦𝜕𝐵𝑥 𝜕𝑦sin(𝜙)cos(𝜃) −𝜇𝑥𝜕𝐵𝑧 𝜕𝑦cos(𝜙)cos(𝜃)+𝜇𝐵𝑦sin(𝜃)+𝜇𝑥𝜕𝐵𝑦 𝜕𝑥sin(𝜃) +𝜇𝑦𝜕𝐵𝑧 𝜕𝑥cos(𝜙)cos(𝜃) (11.11) Differentiating with respect to z the relation provided by the cylindrical symmetry in eq. (11.2) and taking into account the properties of the partial derivatives, 𝑥𝜕𝐵𝑦 𝜕𝑧=𝑦𝜕𝐵𝑥 𝜕𝑧 𝜕𝐵𝑥 𝜕𝑦=𝜕𝐵𝑦 𝜕𝑥 𝜕𝐵𝑥 𝜕𝑧=𝜕𝐵𝑧 𝜕𝑥 𝜕𝐵𝑧 𝜕𝑦=𝜕𝐵𝑦 𝜕𝑧 (11.12) we arrive to d𝐿0𝑍1 d𝑡 =0 QED (11.13) Note that this constant of motion had not been obtained using the equations of motion in any of the previous formulations studied. 11.3. Conservation of 𝑳𝑮𝒁𝟑 The constancy of 𝐿𝐺𝑍3 can be easily deduced seeing that the gravitational force is parallel to 𝐺𝑍3=𝐺𝑍4, and the magnetic force cuts this axis. Moreover, as 𝜇 =−𝜇𝑘 󰇍 3, the magnetic torque 𝑀 󰇍 󰇍 𝐺=𝜇 ×𝐵 󰇍 is perpendicular 11 Integrals of the motion 98 98 to the axis. In contrast to 𝐿0𝑍1, this constant of motion has been obtained in the vectorial mechanics and hamiltonian formulations: 𝑑 𝑑𝑡(𝐼𝑧(𝜓󰇗+𝜙󰇗sin(𝜃)))=0⇒𝐿𝐺𝑍3=𝐼𝑧(𝜓󰇗+𝜙󰇗sin(𝜃))=𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡 (11.14) 12 Analysis including friction 99 99 12 ANALYSIS INCLUDING FRICTION Once the top is hovering over the base, air friction causes the loss of energy of the spinning top. If spin velocity decreases under the lower bound of stability (remember the stability region in figure 9-3), hovering ends and the top falls. We could not finish this study without creating a model taking into account the effects of the air. These effects are going to be considered computing a drag force, a drag torque and Magnus force and introducing them in the Hamiltonian system as a generalized force. 𝑑𝑞𝑖 𝑑𝑡=𝜕ℋ 𝜕𝑝𝑖 𝑑𝑝𝑖 𝑑𝑡=−𝜕ℋ 𝜕𝑞𝑖+𝑄𝑖 (12.1) where 𝑄𝑖 is the generalized force associated with the generalized coordinate 𝑞𝑖. The expressions of these forces were stabilished by Rubinow et al. in [10] for a sphere moving in a viscous fluid. In short, the results are the Stokes force, which is proportional to velocity 𝐹 𝐷=−𝑘𝑣 (12.2) a force which depends on the cross product, known as Magnus force 𝐹 𝑀=𝐴 𝜔 󰇍 󰇍 ×𝑣 (12.3) and a torque against the spinning 𝑀 󰇍 󰇍 𝐷=−𝐵𝜔 󰇍 󰇍 (12.4) Note that there is no torque caused by 𝑣 . In these results, constants k, A and B have known values for a sphere and it has been assummed low Reynolds number. To study the effect of friction on our spinning top, we have to consider the motion in air, so that these expressions can be modelled as quadratic functions. For this case, Luckerchenko et al. in [11] provide the expressions of drag force, Magnus force and drag torque. Generalized forces can be obtained from the computation of the virtual work δW of the applied forces. 12.1. Drag force The expression of this force is 𝐹 𝐷=−34𝜌𝐶𝑑 𝑑|𝑣 |𝑣 =−𝑘|𝑣 |𝑣 (12.5) 12 Analysis including friction 100 100 where 𝜌 is the air density, 𝑑 is the diameter of the sphere and 𝐶𝑑 is a drag coefficient, dependant of the traslational Reynolds number. All these constants are going to be grouped in the constant 𝐶𝑑𝑟𝑎𝑔. Therefore, this force can be written as 𝐹 𝐷=−𝐶𝑑𝑟𝑎𝑔√𝑝𝑥2+𝑝𝑦2+𝑝𝑧2(𝑝𝑥 𝑝𝑦 𝑝𝑧) (12.6) Computing its virtual work 𝛿𝑊=𝐹 𝐷·𝛿𝑟 =𝑄𝐷·𝛿𝑞 (12.7) we get 𝑄𝐷=(𝑄𝐷𝑋 𝑄𝐷𝑌 𝑄𝐷𝑍)=−𝐶𝑑𝑟𝑎𝑔√𝑝𝑥2+𝑝𝑦2+𝑝𝑧2(𝑝𝑥 𝑝𝑦 𝑝𝑧) (12.8) 12.2. Magnus force The usual explanation of this phenomenon is that the rotating body produces a swirl of fluid around it, modifying the uniform flow of fluid. As a result, the pressure is diminished on one side of the body and increased on the other side, giving a side thrust. Figure 12-1: Magnus Force on a rotating sphere. Following Lukerchenko et al. [11], the lateral Magnus force per unit volume is 𝐹 𝑀=𝐶𝑚𝜌𝜔 󰇍 󰇍 ×𝑣 (12.9) where 𝐶𝑚 is a Magnus coefficient, whose expression can be found in [15]. Grouping again all the constants in 𝐶𝑚𝑎𝑔, we have 12 Analysis including friction 101 101 𝐹 𝑀=𝐶𝑚𝑎𝑔𝜔 󰇍 󰇍 ×𝑣 =𝐶𝑚𝑎𝑔|𝑖 1𝑗 1𝑘 󰇍 1 𝜔𝑋𝜔𝑌𝜔𝑍 𝑣𝑋𝑣𝑌𝑣𝑍|=𝐶𝑚𝑎𝑔(𝜔𝑌𝑣𝑍−𝜔𝑍𝑣𝑌 𝜔𝑍𝑣𝑋−𝜔𝑋𝑣𝑍 𝜔𝑋𝑣𝑌−𝜔𝑌𝑣𝑋) (12.10) To compute the expression of the force, we express the angular velocity in frame ‘1’ 𝜔 󰇍 󰇍 |1=( 𝜙󰇗+𝜓󰇗sin(𝜃) 𝜃󰇗cos(𝜙)−𝜓󰇗sin(𝜙)cos(𝜃) 𝜃󰇗sin(𝜙)+𝜓󰇗cos(𝜙)cos(𝜃)) (12.11) and rewriting 𝜙󰇗, 𝜃󰇗 and 𝜓󰇗 in terms of the conjugate momenta { 𝜙󰇗=𝑝𝜙−𝑝𝜓sin(𝜃) 𝐼𝑥cos2(𝜃) 𝜃󰇗=𝑝𝜃 𝐼𝑋 𝜓󰇗=𝑝𝜓 𝐼𝑧−𝑝𝜙−𝑝𝜓sin(𝜃) 𝐼𝑥cos2(𝜃)sin(𝜃) (12.12) the components of 𝜔 󰇍 󰇍 |1 are { 𝜔𝑋=𝑝𝜙 𝐼𝑥+(1 𝐼𝑍−1 𝐼𝑋)𝑝𝜓sin(𝜃) 𝜔𝑌=𝑝𝜃 𝐼𝑋cos(𝜙)+𝑝𝜙 𝐼𝑥sin(𝜙)tan(𝜃)−(cos(𝜃) 𝐼𝑍+sin(𝜃)tan(𝜃) 𝐼𝑋)𝑝𝜓sin(𝜙) 𝜔𝑍=𝑝𝜃 𝐼𝑋sin(𝜙)−𝑝𝜙 𝐼𝑥cos(𝜙)tan(𝜃)+(cos(𝜃) 𝐼𝑍+sin(𝜃)tan(𝜃) 𝐼𝑋)𝑝𝜓cos(𝜙) (12.13) Computing its virtual work 𝛿𝑊=𝐹 𝑀·𝛿𝑟 =𝑄𝑚·𝛿𝑞 (12.14) The generalized force corresponding with the Magnus Force is therefore 𝑄𝑚=𝐶𝑚𝑎𝑔(𝜔𝑌𝑝𝑍−𝜔𝑍𝑝𝑌 𝜔𝑍𝑝𝑋−𝜔𝑋𝑝𝑍 𝜔𝑋𝑝𝑌−𝜔𝑌𝑝𝑋) (12.15) 12.3. Drag Torque The expression provided again by [11] for the drag torque is 12 Analysis including friction 102 102 𝑀 󰇍 󰇍 𝐷=−𝐶𝜔𝜌2𝑟5|𝜔 󰇍 󰇍 |𝜔 󰇍 󰇍 =−𝐶𝑡𝑜𝑟𝑞𝑢𝑒|𝜔 󰇍 󰇍 |𝜔 󰇍 󰇍 (12.16) where we have grouped in 𝐶𝑡𝑜𝑟𝑞𝑢𝑒 the constants 𝐶𝜔 (dependant of the rotational Reynolds number), 𝜌 and the radius of the sphere r. The virtual work of the torque is given by 𝛿𝑊=𝑀 󰇍 󰇍 𝐷·𝛿𝜑 󰇍 =𝑀 󰇍 󰇍 𝐷|3𝜔 󰇍 󰇍 |3 ·𝑑𝑡=𝑄𝑀·𝛿𝑞 (12.17) where 𝜔 󰇍 󰇍 |3 ·𝑑𝑡 is 𝜔 󰇍 󰇍 |3 𝑑𝑡=(𝜙󰇗cos(𝜃)𝑖 3+𝜃󰇗𝑗 3+(𝜓󰇗+𝜙󰇗sin(𝜃))𝑘 󰇍 3)𝑑𝑡 (12.18) Note that in the computation of 𝛿𝑊vectors have been expressed in frame ‘3’ to take advantage, once again, of the rotational symmetry of the top. Using the relations {𝛿𝜙=𝜙󰇗𝑑𝑡; 𝛿𝜃=𝜃󰇗𝑑𝑡; 𝛿𝜓=𝜓󰇗𝑑𝑡} we have 𝜔 󰇍 󰇍 |3 𝑑𝑡=𝛿𝜙cos(𝜃)𝑖 3+𝛿𝜃𝑗 3+(𝛿𝜓+𝛿𝜙sin(𝜃))𝑘 󰇍 3 (12.19) Therefore, virtual work is 𝛿𝑊=−𝐶𝑡𝑜𝑟𝑞𝑢𝑒|𝜔 󰇍 󰇍 |((𝜙󰇗+𝜓󰇗sin(𝜃))𝛿𝜙+𝜃󰇗𝛿𝜃+(𝜓󰇗+𝜙󰇗sin(𝜃))𝛿𝜓) =𝑄𝑀𝜙𝛿𝜙+𝑄𝑀𝜃𝛿𝜃+𝑄𝑀𝜓𝛿𝜓 (12.20) and the components of the generalized torque are {𝑄𝑀𝜙=−𝐶𝑡𝑜𝑟𝑞𝑢𝑒|𝜔 󰇍 󰇍 |(𝜙󰇗+𝜓󰇗sin(𝜃)) 𝑄𝑀𝜃=−𝐶𝑡𝑜𝑟𝑞𝑢𝑒|𝜔 󰇍 󰇍 |𝜃󰇗 𝑄𝑀𝜓=−𝐶𝑡𝑜𝑟𝑞𝑢𝑒|𝜔 󰇍 󰇍 |(𝜓󰇗+𝜙󰇗sin(𝜃)) (12.21) In terms of the conjugate momenta, the generalized torque results in: { 𝑄𝑀𝜙=−𝐶𝑡𝑜𝑟𝑞𝑢𝑒|𝜔 󰇍 󰇍 |(𝑝𝜙 𝐼𝑥+(1 𝐼𝑍−1 𝐼𝑋)𝑝𝜓sin(𝜃)) 𝑄𝑀𝜃=−𝐶𝑡𝑜𝑟𝑞𝑢𝑒|𝜔 󰇍 󰇍 |𝑝𝜃 𝐼𝑥 𝑄𝑀𝜓=−𝐶𝑡𝑜𝑟𝑞𝑢𝑒|𝜔 󰇍 󰇍 |𝑝𝜓 𝐼𝑥 (12.22) 12 Analysis including friction 103 103 12.4. Set of equations considering air friction The Hamiltonian system, including the generalized forces computed in the previous sections, is x󰇗∗=𝑝𝑥∗ y󰇗∗=𝑝𝑦∗ z󰇗∗=𝑝𝑧∗ 𝜙󰇗∗=𝑅𝑝𝜙 ∗−𝑝𝜓 ∗sin(𝜃∗) cos2(𝜃∗) 𝜃󰇗∗=𝑅𝑝𝜃∗ 𝜓󰇗∗=𝑅𝜎𝑝𝜓 ∗−𝑅sin(𝜃∗) cos2(𝜃∗)(𝑝𝜙 ∗−𝑝𝜓 ∗sin(𝜃∗)) 𝑝𝑥∗󰇗=−𝛽(−12𝑑𝐵0∗ 𝑑𝑧∗sin(𝜃∗)−𝑥∗ 2 𝑑2𝐵0∗ 𝑑𝑧∗2cos(𝜙∗)cos(𝜃∗))−𝐶𝑑𝑟𝑎𝑔𝑝𝑥∗ √𝑝𝑥∗2+𝑝𝑦∗2+𝑝𝑧∗2 +𝐶𝑚𝑎𝑔(𝜔𝑌∗𝑝𝑧∗−𝜔𝑍∗𝑝𝑦∗) 𝑝𝑦∗󰇗=−𝛽(12𝑑𝐵0∗ 𝑑𝑧∗sin(𝜙∗)cos(𝜃∗)−𝑦∗ 2 𝑑2𝐵0∗ 𝑑𝑧∗2cos(𝜙∗)cos(𝜃∗))−𝐶𝑑𝑟𝑎𝑔𝑝𝑦∗ √𝑝𝑥∗2+𝑝𝑦∗2+𝑝𝑧∗2 +𝐶𝑚𝑎𝑔(𝜔𝑍∗𝑝𝑥∗−𝜔𝑋∗𝑝𝑧∗) 𝑝𝑦∗󰇗=−𝛽(12𝑑𝐵0∗ 𝑑𝑧∗sin(𝜙∗)cos(𝜃∗)−𝑦∗ 2 𝑑2𝐵0∗ 𝑑𝑧∗2cos(𝜙∗)cos(𝜃∗))−𝐶𝑑𝑟𝑎𝑔𝑝𝑦∗ √𝑝𝑥∗2+𝑝𝑦∗2+𝑝𝑧∗2 +𝐶𝑚𝑎𝑔(𝜔𝑍∗𝑝𝑥∗−𝜔𝑋∗𝑝𝑧∗) 𝑝𝑧∗󰇗=−1−𝛽(−12𝑥∗ 𝑑2𝐵0∗ 𝑑𝑧∗2sin(𝜃∗) +12𝑦∗ 𝑑2𝐵0∗ 𝑑𝑧∗2sin(𝜙∗)cos(𝜃∗) +𝑑𝐵0∗ 𝑑𝑧∗cos(𝜙∗)cos(𝜃∗)−𝑥∗2+𝑦∗2 4 𝑑3𝐵0∗ 𝑑𝑧∗3cos(𝜙∗)cos(𝜃∗)) −𝐶𝑑𝑟𝑎𝑔𝑝𝑧∗ √𝑝𝑥∗2+𝑝𝑦∗2+𝑝𝑧∗2+𝐶𝑚𝑎𝑔(𝜔𝑋∗𝑝𝑦∗−𝜔𝑌∗𝑝𝑥∗) 𝑝𝜙 ∗󰇗=−𝛽(12𝑦∗𝑑𝐵0∗ 𝑑𝑧∗cos(𝜙∗)cos(𝜃∗)−𝐵0∗sin(𝜙∗)cos(𝜃∗) +𝑥∗2+𝑦∗2 4 𝑑2𝐵0∗ 𝑑𝑧∗2sin(𝜙∗)cos(𝜃∗)) −𝐶𝑡𝑜𝑟𝑞𝑢𝑒|𝜔 󰇍 󰇍 |(𝑅𝑝𝜙 ∗+𝑅(1𝜎−1)𝑝𝜓 ∗sin(𝜃∗)) 𝑝𝜃∗󰇗=𝑅 cos3(𝜃∗)(𝑝𝜓 ∗sin(𝜃∗)−𝑝𝜙 ∗)(𝑝𝜙 ∗sin(𝜃∗)−𝑝𝜓 ∗) −𝛽(−12𝑥∗𝑑𝐵0∗ 𝑑𝑧∗cos(𝜃∗)−12𝑦∗𝑑𝐵0∗ 𝑑𝑧∗sin(𝜙∗)sin(𝜃∗)−𝐵0∗cos(𝜙∗)sin(𝜃∗) +𝑥∗2+𝑦∗2 4 𝑑2𝐵0∗ 𝑑𝑧∗2cos(𝜙∗)sin(𝜃∗))−𝐶𝑡𝑜𝑟𝑞𝑢𝑒|𝜔 󰇍 󰇍 |𝑅𝑝𝜃∗ 𝑝𝜓 ∗󰇗=−𝐶𝑡𝑜𝑟𝑞𝑢𝑒|𝜔 󰇍 󰇍 |𝑅𝑝𝜓 ∗ (12.23) 13 Let’s play! 110 110 14 Conclusions and future steps 111 111 14 CONCLUSIONS AND FUTURE STEPS During this project, an analytical study of the hovering magnetic system (known as Levitron) has been developed: - In chapter 1, we describe the typical contents of a Levitron kit, how the toy works and briefly expose the history of the device. - In chapter 2, an analysis of the free-motion, in absence of the external magnetic field and solely submitted to the action of gravity is performed using the Classical Euler angles. We derive the constants of the free-motion and describe the motion in term of the axoids (also known as body and space cones). We also see that Classical Euler angles provide a nice visualization of the problem, but the set of equations is singular for 𝜃=0°. - In chapter 3, we present the Tait-Bryan angles, describing the system with a yaw-pitch-roll sequence in order to solve the singularity previously mentioned. We obtained the relation between the fixed frame ‘1’ and the local body-frame ‘4’, and computed the angular velocity 𝜔 󰇍 󰇍 , kinetic moment 𝐿 󰇍 𝐺 and kinetic energy 𝑇 using this set of angles. - In chapter 4, we find the expression of the magnetic field at 0𝑍1 axis, derive the equations of the magnetic field 𝐵 󰇍 for points near the axis and first approximate the spinning top as a magnetic dipole of magnetic moment 𝜇 . - In chapter 5, we analyze the static stability of the system. Firstly, the potential energy is computed and we see that its expression validates Earnshaw’s theorem. Nevertheless, we see that the Laplacian of the potential energy considering the alignment between 𝜇 and 𝐵 󰇍 because of gyroscopic effect is positive, so that the previous theorem is not verified. Finally, we compute the potential energy 0𝑍1 axis and found the minimum height for which hovering is posible. - In chapter 6, the equations of motion of the system are derived using vector mechanics. Newton’s and Euler’s Second law provide a second-order differential set of equations, that is also rewritten as a firstorder set. - In chapter 7, we get the same set of equations using the Lagrangian formulation of the analytical mechanics and Hamilton’s equations are derived. The Hamiltonian framework allow us to work with a first-order set of ten equations, as we can ignore the equation for ψ. - In chapter 8, nondimensional variables are defined in order to perform numerical integrations and we give the expressions for the vector mechanics and Hamilton’s set of equations in term of these variables. - In chapter 9, the analysis of the linear stability is performed. We linearize the first-order set of equations about the position of equilibrium in the 0𝑍1 axis, solve the linear system with real variables, complex exponentials and complex variables and find a basis of solutions to express any other as a combination of them. Moreover, we derive mass, stiffness and gyroscopic matrixes to describe the system and its characteristic equation, that allow us to find the stability region in which stable hovering is possible. We also numerically solve two examples for different height of levitation, obtaining Campbell’s diagrams and the natural frequencies for a certain rotational speed. 14 Conclusions and future steps 112 112 Computing the eigenvectors of the linear matrix, we explain the interpretation of the complex conjugate eigenvectors and find the five normal modes of the system. An analysis of these modes is also done, showing their trajectories and the time evolution. Finally, we numerically show the trajectories for two given initial conditions and numerically validate the stable rotational speed for a height of levitation centered in the stability region. - In chapter 10, we study the nonlinear system, verifying the linear and the nonlinear coupling. Moreover, we numerically solve the set of equations with initial conditions that try to reproduce the real situations that one faces when plays with the toy. - In chapter 11, constants of the motion are presented. We numerically check the constancy of the mechanical energy and demonstrate the constancy of 𝐿0𝑍1 and 𝐿𝐺𝑍3. - In chapter 12, a simple model is presented to consider the effect of air friction. We compute the drag force, Magnus force and drag torque and their corresponding generalized forces are included in the Hamiltonian system. The values of three coefficients, 𝐶𝑑𝑟𝑎𝑔, 𝐶𝑚𝑎𝑔 and 𝐶𝑡𝑜𝑟𝑞𝑢𝑒 are found in order to provide a time of levitation around one minute. - Finally, in chapter 13 we give some instructions in order to use the toy properly and we present a flowchart that summarizes all the situations that can be found when playing with the toy, indicating the corresponding chapter where we numerically study the situation. Some recommendations for a future research should be given: - Our analytical study is valid near the 0𝑍1 axis, as the expressions of the magnetic field can only be used for 𝜌≪1. Therefore, a study using the exact expressions of the magnetic field should be done. - We have considered the spinning top as a magnetic dipole, neglecting its real size. This is another effect that should be taken into account in next studies. - A more detailed model of the air friction is another possibility for future works. 14 Conclusions and future steps 113 113 Addendum I: Experimental adjustment 114 114 ADDENDUM I: EXPERIMENTAL ADJUSTMENT To make the experimental adjustment of the magnetic field, the experimental data taken in the job ‘Efectos de pequeñas corrientes de aire sobre un Levitron’ [8], which was obtained using a teslameter, is going to be used: z(mm) B(mT) z(mm) B(mT) z(mm) B(mT) z(mm) B(mT) 10 8,40 31 7,36 52 7,29 73 5,01 11 -7,21 32 7,54 53 7,17 74 4,91 12 -5,61 33 7,70 54 7,05 75 4,83 13 -4,55 34 7,82 55 6,95 76 4,73 14 -3,14 35 7,90 56 6,83 77 4,63 15 -1,84 36 7,98 57 6,74 78 4,53 16 -1,01 37 8,03 58 6,64 79 4,46 17 0,10 38 8,06 59 6,52 80 4,37 18 0,91 39 8,07 60 6,40 81 4,27 19 1,73 40 8,07 61 6,29 82 4,18 20 2,52 41 8,06 62 6,16 83 4,10 21 3,25 42 8,01 63 6,05 84 4,03 22 3,90 43 7,98 64 5,95 85 3,95 23 4,46 44 7,93 65 5,85 86 3,86 24 5,05 45 7,87 66 5,74 87 3,79 25 5,54 46 7,78 67 5,64 88 3,70 26 5,87 47 7,72 68 5,53 89 3,63 27 6,25 48 7,65 69 5,44 90 3,58 28 6,52 49 7,55 70 5,32 91 3,49 29 6,91 50 7,47 71 5,20 30 7,16 51 7,39 72 5,11 One has to find the values of the parameters 𝐴, 𝑏 and 𝑧0 in the following expression: Addendum I: Experimental adjustment 115 115 𝐵 󰇍 0(𝑧)=𝐴(𝑏2 (𝑏2+(𝑧−𝑧0)2)3/2 −𝑎2 (𝑎2+(𝑧−𝑧0)2)3/2 )𝑘 󰇍 1 1- The inner radius of the platform 𝑎 is approximately equal to 25 mm. 2- To calculate 𝑏 and 𝑧0, we solve the set of equations given by (AI.1) and (AI.2) 𝐵0=0 (AI.1) 𝜕𝐵0 𝜕𝑧=0 (AI.2) obtaining 𝑏=67.6864 mm and 𝑧0=−11.4067 mm 3- Parameter A can be obtained evaluating those points where 𝐵0 reaches its maximum value. Hence, 𝐴=1.9608·103 mm·mT In figure, we can see a comparison between the distribution of points experimentally obtained and the adjustment given by eq. Figure AdI-1: Experimental adjustment. Addendum I: Experimental adjustment 116 116 Addendum II: Analysis with quaternions and euler parameters 117 117 ADDENDUM II: ANALYSIS WITH QUATERNIONS AND EULER PARAMETERS In the previous chapters, we studied the Euler angles to relate the angular velocities to the rotation angles 𝜙, 𝜃, 𝜓 and their rates. These equations have singularities, which makes it difficult to do any numerical work with them. Even with the advantage of the yaw-pitch-roll sequence, which eliminated the problem in the denominators by shifting the singularity to 𝜃=𝜋2, the equations are highly nonlinear. To alleviate these difficulties, we are going to express the motion of system with another set, called Euler parameters. These parameters increase the number of variables from three to four, but in many cases they eliminate the nonlinearities and many of the numerical problems. Because a set of three variables (𝜙, 𝜃,𝜓) is being expressed in terms of four variables, the use of Euler parameters introduces a redundancy. This implies that there is no unique way of expressing the Euler parameters, existing other commonly used quantities which are, essentially, different forms of these parameters. The vector formulation, known as quaternions, was developed by Heaviside. Nevertheless, from the inspection of the Hamiltonian system of equations derived in chapter 7.2, one can see that there are some equations with cosine functions in the denominators. This fact allays the significant advantage provided by the use of quaternions. The most important ideas and expressions of quaternions used to derive the equations of motion are - Generalization of complex numbers to 4D q=𝑞0+𝑞 =𝑞0+𝑞𝑥𝑖 +𝑞𝑦𝑗 +𝑞𝑧𝑘 󰇍 - Conjugate quaternion: q∗=𝑞0−𝑞 - Scalar and vector parts of a quaternion q: 𝑞0=𝑆(q)=q+q∗ 2 𝑞 =𝑉(q)=q−q∗ 2 - Product of two quaternions q and p: qp =(𝑞0𝑝0−𝑞 ·𝑝 )+(𝑞0𝑝 +𝑝0𝑞 +𝑞 ×𝑝 ) - Modulus of a quaternion |q|=√qq∗=√𝑞02+|𝑞 |2 Addendum II: Analysis with quaternions and euler parameters 118 118 - Inverse q−1=q∗ |q|2 - The direction cosine matrix 𝑅, in terms of the parameters {𝑞0,𝑞𝑥,𝑞𝑦,𝑞𝑧} has quadratic expressions 𝑅=(𝑞02+𝑞𝑥2−𝑞𝑦2−𝑞𝑧22(𝑞𝑥𝑞𝑦−𝑞0𝑞𝑧) 2(𝑞𝑥𝑞𝑧+𝑞0𝑞𝑦) 2(𝑞𝑥𝑞𝑦+𝑞0𝑞𝑧) 𝑞02−𝑞𝑥2+𝑞𝑦2−𝑞𝑧22(𝑞𝑥𝑞𝑦−𝑞0𝑞𝑧) 2(𝑞𝑥𝑞𝑧−𝑞0𝑞𝑦) 2(𝑞𝑦𝑞𝑧+𝑞0𝑞𝑥) 𝑞02−𝑞𝑥2−𝑞𝑦2+𝑞𝑧2) - The consideration of unit quaternions leads to |q|=1⇒𝑞02+𝑞𝑥2+𝑞𝑦2+𝑞𝑧2=1 Following Euler’s rotation theorem, a rotation through an angle 𝜃 about an axis defined by the unit vector 𝑢 󰇍 Figure AdII-1: Euler’s rotation theorem. can be written as q𝜃=cos(𝜃2)+sin(𝜃2)𝑢 󰇍 Denoting as 𝑟 0 and 𝑟 the initial and final positions, we have the following relations 𝑟 =q𝑟 0q∗ 𝑟 0=q∗𝑟 q As we saw in chapter 3, we use the yaw-pitch-roll sequence to describe the motion of the spinning top, so that we arrive from frame ‘1’ to the body-frame ‘4’ by means of the successive rotations: Addendum II: Analysis with quaternions and euler parameters 119 119 Figure AdII-2: Successive rotations in terms of the quaternions. In terms of the quaternions, we can compound these three rotations by multiplying the quaterions: q=q𝜙q𝜃q𝜓 Differentiating with respect to time 𝑟 =q𝑟 0q∗ we have 𝑟 󰇗=q󰇗𝑟 0q∗+q𝑟 0q󰇗∗=q󰇗q∗𝑟 +𝑟 qq󰇗∗=𝑉(2q󰇗q∗𝑟 ) The angular velocity can be written as 𝜔 󰇍 󰇍 =2q󰇗q∗ Therefore, we arrive to the following relation between the time derivatives of Euler parameters and the angular velocity components. q󰇗=12𝜔 󰇍 󰇍 q As we saw in chapter 6.1, the traslational motion of the system is given by Newton’s Second law 𝑚𝑎 =𝑚𝑔 +𝛻(𝜇 ·𝐵 󰇍 ) and the rotational motion, studied in 6.2 d𝐿 󰇍 𝐺 d𝑡=𝜇 ×𝐵 󰇍 The components of the magnetic field 𝐵 󰇍 and the magnetic moment 𝜇 verify 𝐵 󰇍 =𝐵𝑥𝑖 1+𝐵𝑦𝑗 1+𝐵𝑧𝑘 󰇍 1 𝜇 =−𝜇𝑘 󰇍 4=−𝜇𝑅−1𝑘 󰇍 1 Addendum III: Description of other systems 126 126 Figure AdIII-5:Plot of 𝑓(𝑢) The precession rate, from eq. is 𝜙󰇗=𝑎−𝑏𝑢 1−𝑢2 (AdIII.16) so that this rate depends on the value 𝑎−𝑏𝑢. Defining the quantity 𝑢0, based on initial conditions as 𝑢0=𝑝𝜙 𝑝𝜓=𝑎𝑏=𝐼𝑋sin2(𝜃0)𝜙󰇗0 𝐼𝑍(𝜙󰇗0cos(𝜃0)+𝜓󰇗0) (AdIII.17) It turns out that characteristics of the motion of the top depend on 𝑢0, 𝑢1 and 𝑢2. We have the following different types of motion based on the initial conditions: Unidirectional precession Looping precession Cuspidal motion 𝑢0>𝑢2 𝑜𝑟 𝑢0<𝑢1 𝑢1<𝑢0<𝑢2 𝑢0=𝑢2 Figure AdIII-6: Types of precession. Addendum III: Description of other systems 127 127 AdIII.2 Forced precession and nutation of Earth In this section, we follow [22] in order to obtain a model to describe precession and nutation motions and compute the precession and nutation periods of our planet. To do this, we use again the proper or classic Euler angles (chapter 2.1). Figure AdIII-7: Earth motions in term of the Classic Euler angles. We start considering the Sun-Earth system and then we increase the accuracy of the model adding the moon influence. The potential energy of the Earth-Sun system, using McCullough’s formula, is 𝑈=−𝐺𝑀𝑆𝑀 𝑎𝑆+𝐺𝑀𝑆(𝐼𝑍−𝐼𝑋) 𝑎𝑆3𝑃2(cos𝛾𝑆) (AdIII.18) where 𝑀𝑆 is the mass of the Sun, 𝑀 the mass of the Earth, 𝐼𝑍 the Earth’s moment of inertia about its axis of rotation and 𝐼𝑋 the Earth’s moment of inertia about an axis lying on its equatorial plane. Function 𝑃2(𝑥) is given by 𝑃2(𝑥)=12(3𝑥2−1) (AdIII.19) and 𝛾𝑆 is the angle subtended between 𝜔 󰇍 󰇍 and the position vector of the Sun relative to the Earth. Being 𝑎𝑆 the radius of the apparent orbit and 𝜆𝑆 the ecliptic longitude , 𝑟 𝑆 is given by 𝑟 𝑆=𝑎𝑆(cos𝜆𝑆,sin𝜆𝑆,0) (AdIII.20) Computing the dot product to obtain cos𝛾𝑆, we have 𝜔 󰇍 󰇍 ⋅𝑟 𝑆 =|𝜔 󰇍 󰇍 ||𝑟 𝑆|cos𝛾𝑆→cos𝛾𝑆=𝜔 󰇍 󰇍 ⋅𝑟 𝑆 |𝜔 󰇍 󰇍 ||𝑟 𝑆|=sin(𝜃)sin(𝛾𝑆) (AdIII.21) Addendum III: Description of other systems 128 128 In figure AdIII-8 we can see all the magnitudes previously defined: Figure AdIII-9: Sun-Earth system. Potential energy is therefore 𝑈=−𝐺𝑀𝑆𝑀 𝑎𝑆+𝐺𝑀𝑆(𝐼𝑍−𝐼𝑋) 2𝑎𝑆3(3sin2(𝜃)sin2(𝜆𝑆)−1) (AdIII.22) As we are primarily interested in the motion of the Earth’s axis of rotation on timescales that are much longer than a year, we can average the expression over the Sun’s orbit, using that the average of sin2(𝜆𝑆) over a year is ½ 𝑈=−𝐺𝑀𝑆𝑀 𝑎𝑆+𝐺𝑀𝑆(𝐼𝑍−𝐼𝑋) 2𝑎𝑆3(32sin2(𝜃)−1) (AdIII.23) After some mathematical manipulation, we get 𝑈=−𝐺𝑀𝑆𝑀 𝑎𝑆−3𝐺𝑀𝑆𝐼𝑍(𝐼𝑍−𝐼𝑋) 8𝑎𝑆3𝐼𝑍(13+cos(2𝜃)) (AdIII.24) The preceding expression can be rewritten as 𝑈=𝑈0−𝜖𝛼𝑆cos(2𝜃) (AdIII.25) 𝑈0 is a constant that will disappear when we derive the equations of motion. Moreover, the following constants have been introduced: 𝛼𝑆=38𝐼𝑍ω𝑆2 where ω𝑆 is the Sun’s apparent orbital angular velocity, whose expression can be obtained by equating the gravitational force to the centripetal force of the Earth Addendum III: Description of other systems 129 129 𝐺𝑀𝑆𝑀 𝑎𝑆2=𝑀𝜔𝑆2𝑎𝑆→𝜔𝑆=𝑑𝜆𝑆 𝑑𝑡=(𝐺𝑀𝑆 𝑎𝑆3)12 Finally, 𝜖 is the Earth’s dynamical ellipticity and its value is determined from the Earth’s observed flattening 𝜖=𝐼𝑍−𝐼𝑋 𝐼𝑍 The rotational kinetic energy is 𝑇=12(𝐼𝑋ω𝑋 2+𝐼𝑌ω𝑦2+𝐼𝑍ω𝑍2) (AdIII.26) and in terms of the Classic Euler angles 𝑇=12[𝐼𝑋(𝜃󰇗2+𝜙󰇗2sin2(𝜃))+𝐼𝑍(𝜓󰇗+𝜙󰇗cos(𝜃))2] (AdIII.27) Hence, neglecting any constant terms, the Earth’s Lagrangian is ℒ=12[𝐼𝑋(𝜃󰇗2+𝜙󰇗2sin2(𝜃))+𝐼𝑍(𝜓󰇗+𝜙󰇗cos(𝜃))2]+𝜖𝛼𝑆cos(2𝜃) (AdIII.28) As ℒ does not depend on the angular coordinate 𝜓, it follows that the corresponding conjugate momentum is a constant of the motion, implying that ω𝑍 is also a constant of the motion. Note that ω𝑍 is the Earth’s angular velocity of rotation about its axis 𝑝𝜓=𝜕ℒ 𝜕𝜓󰇗=𝐼𝑍ω𝑍 Using Lagrange equation for 𝜃 𝑑 𝑑𝑡(𝜕ℒ 𝜕𝜃󰇗)−𝜕ℒ 𝜕𝜃=0 we have 𝐼𝑋𝜃󰇘−𝜕ℒ 𝜕𝜃=0 (AdIII.29) To compute the precession period, we consider steady precession of the rotation axis of the Earth about the normal of the ecliptic plane. In terms of the Proper Euler angles, this means 𝜙󰇗=𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡→𝜙󰇘=0 𝜃=𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡→𝜃󰇗=𝜃󰇘=0 𝜓󰇗=𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡→𝜓󰇘=0 and therefore Addendum III: Description of other systems 130 130 𝜕ℒ 𝜕𝜃=0→12𝐼𝑋𝜙󰇗2sin(2𝜃)−𝐼𝑍ω𝑍𝜙󰇗sin(𝜃)−2𝜖𝛼𝑆sin(2𝜃)=0 (AdIII.30) Using that ω𝑍 is much larger than precession speed, i.e |𝜙󰇗|≪𝜔, the first term of the previous equation can be neglected and obtain 𝜙󰇗 as a function of the rest of the variables: 𝜙󰇗=−4𝜖𝛼𝑆cos(𝜃) 𝐼𝑋𝜔=−3𝜖ω𝑆2 2ω𝑍cos(𝜃)=−𝛺𝜙 (AdIII.31) with 𝛺𝜙=3𝜖ω𝑆2 2ω𝑍cos(𝜃) (AdIII.32) Thus, Earth’s axis of rotation describes a steady precession about the normal of the ecliptic plane with speed −𝛺𝜙. Minus sign means that the precession sense is opposite to the Sun’s apparent orbit and the Earth’s diurnal rotation: 𝑇𝜙=ω𝑆 𝛺𝜙=ω𝑆 3𝜖ω𝑆2 2ω𝑍cos(𝜃)=2ω𝑍/ω𝑆 3𝜖cos(𝜃) (AdIII.33) with ω𝑍 ω𝑆=366.26 𝑟𝑒𝑣/𝑦𝑒𝑎𝑟 1 𝑟𝑒𝑣/𝑦𝑒𝑎𝑟 =366.26 For the rest of the parameters 𝜖=0.00335 𝜃=23.44° (mean value) Therefore, we get 𝑇𝜙≃79400 years Nevertheless, the observed precession period of the Earth’s axis of rotation about the normal to the ecliptic plane is approximately 25.800 years, so we need to refine the model considering the influence of the Moon. The potential energy of the Earth-Moon system is 𝑈=−𝐺𝑀𝑚𝑀 𝑎𝑚+𝐺𝑀𝑚(𝐼𝑍−𝐼𝑋) 𝑎𝑚 3𝑃2(cos𝛾𝑚) (AdIII.34) Addendum III: Description of other systems 131 131 where 𝑀𝑚 is the lunar mass, 𝑎𝑚the radius of the Moon’s orbit (assuming that is approximately circular) and 𝛾𝑚 is the angle subtended between the Earth’s angular velocity 𝜔 󰇍 󰇍 and the position vector of the Moon relative to the Earth 𝑟 𝑚. The Moon’s orbital plane is slightly inclined to the ecliptic plane, being the angle of inclination 𝐼𝑚=5.16° 𝑟 𝑚=𝑎𝑚(cos𝜆𝑚,sin𝜆𝑚,sin (𝐼𝑚)sin(𝜆𝑚−𝛼𝑛)) Approximating to first order in 𝐼𝑚, we can write: 𝑟 𝑚≃𝑎𝑚(cos𝜆𝑚,sin𝜆𝑚,𝐼𝑚sin(𝜆𝑚−𝛼𝑛)) where 𝜆𝑚 is the Moon’s ecliptic longitude and 𝛼𝑛 is the ecliptic longitude of the lunar ascending node, which is defined as the point of the lunar orbit where the Moon crosses the ecliptic plane from south to north. All these variables can be seen in figure AdIII-10 and AdIII-11: Figure AdIII-12: Sun-Earth-Moon system. Figure AdIII-13: Sun-Earth-Moon system. Computing Moon’s apparent angular velocity 𝜔𝑚 and cos𝛾𝑚 Addendum III: Description of other systems 132 132 𝜔𝑚=𝑑𝜆𝑚 𝑑𝑡 =(𝐺𝑀 𝑎𝑚 3)12 𝜔 󰇍 󰇍 ⋅𝑟 𝑚=|𝜔 󰇍 󰇍 ||𝑟 𝑚|cos𝛾𝑚→cos𝛾𝑚=𝜔 󰇍 󰇍 ⋅𝑟 𝑚 |𝜔 󰇍 󰇍 ||𝑟 𝑚| =sin(𝜃)sin(𝜆𝑚−𝜙)+𝐼𝑚cos(𝜃)sin(𝜆𝑚−𝛼𝑛) (AdIII.35) The potential energy of Moon’s-Earth system is 𝑈=−𝐺𝑀𝑚𝑀 𝑎𝑚+𝐺𝑀𝑚(𝐼𝑍−𝐼𝑋) 2𝑎𝑚 3[3sin2(𝜃)sin2(𝜆𝑚−𝜙) +3𝐼𝑚sin(2𝜃)sin(𝜆𝑚−𝜙)sin(𝜆𝑚−𝛼𝑛)−1] (AdIII.36) We are interested in the motion of the Earth’s axis of rotation on timescales that are much longer than a month 𝑈=−𝐺𝑀𝑚𝑀 𝑎𝑚+𝐺𝑀𝑚(𝐼𝑍−𝐼𝑋) 2𝑎𝑚 3[3sin2(𝜃)sin2(𝜆𝑚−𝜙)+3𝐼𝑚sin(2𝜃)sin(𝜆𝑚−𝜙)−1] (AdIII.37) where it has been used that the average of sin2(𝜆𝑚−𝜙)=1/2 sin(𝜆𝑚−𝜙)sin(𝜆𝑚−𝛼𝑛)=12cos(𝛼𝑛−𝜙) Moreover, defining 𝛼𝑚=38𝐼𝑍𝛾𝑚𝜔𝑚 2 𝛽𝑚=34𝐼𝑍𝐼𝑚𝛾𝑚𝜔𝑚 2 Therefore, Moon’s-Earth potential energy can be shortened to 𝑈=𝑈0′−𝜖𝛼𝑚cos(2𝜃)+𝜖𝛽𝑚sin(2𝜃)cos(𝛼𝑛−𝜙) (AdIII.38) To compute the total potential energy of the Earth we use that gravity is a superposable force, so we can add the potential energy of the Earth-Sun and Eath-Moon systems: 𝑈=𝑈0′′−𝜖𝛼𝑆cos(2𝜃)−𝜖𝛼𝑚cos(2𝜃)+𝜖𝛽𝑚sin(2𝜃)cos(𝛼𝑛−𝜙) (AdIII.39) where 𝑈0′′ is another constant Therefore, the Lagrangian of the Earth, considering the influence of Sun and Moon, is ℒ=12[𝐼𝑋(𝜃󰇗2+𝜙󰇗2sin2(𝜃))+𝐼𝑍(𝜓󰇗+𝜙󰇗cos(𝜃))2] +𝜖𝛼cos(2𝜃)−𝜖𝛽𝑚sin(2𝜃)cos(𝛼𝑛−𝜙) (AdIII.40) Addendum III: Description of other systems 133 133 Following in an analogous way to [13], Lagrange’s equations 𝑑 𝑑𝑡(𝜕ℒ 𝜕𝜃󰇗)−𝜕ℒ 𝜕𝜃󰇗=0 𝑑 𝑑𝑡(𝜕ℒ 𝜕𝜙󰇗)−𝜕ℒ 𝜕𝜙=0 provide 𝑇𝜙≃24800 years 𝑇𝜃≃18.6 years Both stimations are quite close to the mark. The inaccuracies are due to small eccentricities of the Earth’s orbit around the Sun and the Moon’s orbit around the Earth, that have been neglected, and the consideration of the Earth as a homogeneous body. AdIII.3 Trajectories of particles in a heli-toroidal magnetic field Nuclear fusion is a magnet-mechanical problem with cylindrical symmetry (as the Levitron case) of great cientific repercussion. We are going to construct a simplified model of a Tokamak reactor in order to study the motion of a single particle. For that purpose, a constant magnetic field acting in the Z-direction 𝐵𝑍 and an azimuthal component 𝐵𝜃 are going to be considered (the real problem is, of course, so much complex). A basic scheme of the reactor can be seen in figures AdIII-14 and Figure AdIII-0-15 Figure AdIII-16: Scheme of a Tokamak reactor. Addendum III: Description of other systems 134 134 Figure AdIII-17: Scheme with constant vertical and azimuthal magnetic fields. We assume that both fields are constants 𝐵𝑍=𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡 1 𝐵𝜃=𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡 2 In a more detailed model, the variable 𝑟 in figure AdIII-18 would represent small variations around the stationary solution (circumference about 𝑍1-axis). To take into account the magnetic field acting in the azimuthal direction, we can use the Ampere’s circuit law, which relates the integrated magnetic field around a closed loop (our simplified Tokamak) to the electric current passing through the loop: ∮𝐵 󰇍 𝑑𝑙 =𝜇0∑𝐼𝑖𝑖𝑛𝑡 𝑖→𝐵𝜃∙2𝜋𝜌=𝜇0𝐼 (AdIII.41) Therefore, 𝐵𝜃=𝑘𝜌 with 𝑘=𝜇0𝐼/2𝜋 Using cylindrical coordinates, 𝐵 󰇍 can be expressed as 𝐵 󰇍 =𝐵𝜃𝑢 󰇍 𝜃+𝐵𝑍𝑘 󰇍 The acceleration of the particle is 𝑎 =(𝜌󰇘−𝜌𝜃󰇗2 𝜌𝜃󰇘+2𝜌󰇗𝜃󰇗 𝑧󰇘 ) (AdIII.42) Addendum III: Description of other systems 135 135 To compute the force acting on the particle of charge 𝑞 due to the magnetic field 𝐵 󰇍 we use the expression of the Lorentz force: 𝐹 =𝑞𝑣 ×𝐵 󰇍 =𝑞|𝑢 󰇍 𝜌𝑢 󰇍 𝜃𝑘 󰇍 𝜌󰇗 𝜌𝜃󰇗𝑧󰇗 0 𝐵𝜃𝐵𝑍|=𝑞(𝜌𝜃󰇗𝐵𝑍−𝑧󰇗𝐵𝜃 −𝜌󰇗𝐵𝑍 𝜌󰇗𝐵𝜃) (AdIII.43) This leads to the second-order differential equations {𝑚(𝜌󰇘−𝜌𝜃󰇗2)=𝑞(𝜌𝜃󰇗𝐵𝑍−𝑧󰇗𝐵𝜃) 𝑚(𝜌𝜃󰇘+2𝜌󰇗𝜃󰇗)=−𝑞𝜌󰇗𝐵𝑍 𝑚𝑧󰇘=𝑞𝜌󰇗𝐵𝜃 (AdIII.44) In order to perform numerical integrations, we rewrite it as a first-order set of equations { 𝜌󰇗=𝑣𝜌 𝜃󰇗=𝑣𝜃 𝑧󰇗=𝑣𝑧 𝑣󰇗𝜌=𝜌𝑣𝜃2 𝑚+𝑞 𝑚(𝜌𝑣𝜃𝐵𝑍−𝑣𝑧𝐵𝜃) 𝑣󰇗𝜃=−2𝑣𝜌𝑣𝜃 𝜌−𝑞𝑣𝜌𝐵𝑍 𝜌 𝑣󰇗𝑧=𝑞𝑣𝜌𝐵𝜃 (AdIII.45) In the following figures, we can see that the motion consists on the superposition of the helical trajectory caused by 𝐵𝑍 and the contribution due to the azimuthal field 𝐵𝜃 Figure AdIII-19: Helical trajectory considering vertical constant magnetic field.