ON SOME ASPECTS OF DYNAMICS, MODELLING, ATTITUDE ANALYSIS OF A N D SATELLITES Said Rashed Marandi .S. (Mechanical Engineering and Mathematics), M.A. (Mathematics), California Institute of Technology, 1981 University of California, Berkeley, 198S A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T O F THE REQUIREMENTS FOR THE D E G R E E OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies Department of Mechanical Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH C O L U M B I A July 1988 © S a i d Rashed M a r a n d i , 1988 In presenting degree freely at this the available copying of department publication of in partial fulfilment University of British Columbia, for this or thesis reference and thesis by this for his thesis study. scholarly or for her financial of The University of British 1956 Main Mall Vancouver, Canada V6T 1Y3 Date DE-6(3/81) I further purposes Columbia JJ\pft<K gain the shall requirements agree that agree may representatives. permission. Department I of be It not that the Library permission granted is by understood be for allowed an advanced shall for the that without make it extensive head of my copying or my written ABSTRACT The thesis identifies several limitations in the modelling and attitude stability analysis of two classes of spacecraft: rigid and flexible satellites. Attractive methods are proposed which promise to have far reaching consequences in spacecraft dynamics. These alternatives, developed based on techniques of differential equations, classical mechanics, and differential topology, are indicated below. (a) A n Alternate Transition from the Lagrangian of a Satellite to Equations of Motion The classical procedure requires the Lagrangian to be expressed in terms of the corresponding generalized coordinates of the problem. This requirement significantly complicates the derivation of the equations of motion through an introduction of a set of librational generalized coordinates, which is strictly not a part of the dynamical system. Using the Lagrangian in the natural variables (angular velocity, direction cosines, and vibrational coordinates), one develops a procedure for derivation of equations of motion without an a priori choice of rotational generalized coordinates. For the case of a satellite with two flexible plate-type appendages, for example, the approach reduced the formulation time to one-third. (b) Synthesis and Depiction of Rotational Motion of Satellites and Robots The rotational coordinates in use for numerical prediction of orientation of a satellite are either singular or redundant. Furthermore, they lack a convenient visual interpretation. A new set of coordinates is proposed and an associated representation is developed which avoids these limitations. The procedure is applied to represent and integrate numerically the librational response of the flexible satellite mentioned in (a). (c) Resolution of Attitude Stability of Delp Satellites The development here tackles a long outstanding problem in the area of attitude stability of satellites. The resolution of this problem through normalization of the Hamiltonian leads to a better appreciation of stability associated with the class of gravity gradient structures such as the proposed Space Station. 11 CONTENTS ABSTRACT ii CONTENTS iii LIST OF T A B L E S vi LIST OF F I G U R E S vii NOMENCLATURE x ACKNOWLEDGEMENT xiv 1. INTRODUCTION 1 2. D E V E L O P M E N T OF T H E L A G R A N G I A N 7 2.1 A Rigid Satellite in a Circular Orbit 7 2.1.1 Introduction 7 2.1.2 Gravitational potential energy 8 2.1.3 The Lagrangian 8 2.2 2.3 A Hybrid Satellite in a Circular Orbit 9 2.2.1 Satellite configuration 9 2.2.2 Computation of the inertia matrix 11 2.2.3 Kinetic energy 14 2.2.4 Discretization procedure 15 2.2.5 Strain energy 18 2.2.6 The Lagrangian 18 References 23 APPENDIX 24 2-1 24 EXPANSION OF POTENTIAL E N E R G Y ••• ill 3. A N A L T E R N A T E TRANSITION F R O M T H E L A G R A N G I A N TO T H E EQUATIONS OF MOTION F O R A LIBRATING S Y S T E M 25 3.1 25 3.2 28 The Energy Function From the Lagrangian 48 3.4 Hamiltonian From the Energy Function 50 3.5 The Canonical 2-form and the Differential of the Hamiltonian . . . 50 3.6 The Expression of the Hamiltonian Vector Field 51 3.7 Summary of the Procedure 53 3.8 Application of the Procedure to a Hybrid Satellite 57 3.9 The Software R E D U C E in Derivation of the Equations of Motion References . 59 61 A NEW LIBRATIONAL GENERALIZED COORDINATE (a-PARAMETERS) 62 4.1 Introduction 62 4.2 Formulation 64 4.3 Interpretation, Properties, and Generalization of the er-Parameters . 67 4.4 Two Illustrative Examples 70 4.5 Representation of Librational Trajectories 72 4.6 References 90 APPENDICES 91 4-1 91 4-II 5. Summary of the Mathematical Tools 3.3 3.10 4. Introduction D I R E C T I O N COSINES A N D Q U A T E R N I O N S ANGULAR VELOCITY AND a-PARAMETERS 93 4-III CONSERVATION OF E N E R G Y 94 4-IV C O N F O R M A L I T Y OF c o p " 95 1 U S E O F T H E E N E R G Y F U N C T I O N IN D E T E R M I N I N G B O U N D S iv O F S T A B I L I T Y F O R RIGID S A T E L L L I E S 98 5.1 Introduction 98 5.2 Formulation 100 5.3 Disturbance Energy and Attitude Perturbation 108 5.4 References 128 APPENDICES 5-1 5- II 6. 129 T H E CRITICAL POINTS OF E 129 p T H E L O C A L EXPRESSIONS OF E p AT EQUILIBRIA . . . . 130 A T T I T U D E S T A B I L I T Y OF RIGID S A T E L L I T E S T H R O U G H N O R M A L I Z A T I O N OF T H E H A M I L T O N I A N 133 6.1 Introduction 133 6.2 The Hamiltonian and the Choice of Generalized Coordinates 6.3 The Hamiltonian up to the Fourth Degree 139 6.4 Normalization Procedure 141 6.4.1 Normalization up to the second degree 141 6.4.2 Elimination of the third degree terms (First method) 6.4.3 Normal form through the fourth degree terms 6.4.4 Normalization beyond the second degree terms—the Lie series . . . . 138 151 153 154 6.5 Stability of an Equilibrium—the Delp Case 158 6.6 References 164 APPENDIX 167 6- 1 H I G H E R T E R M S O F T H E H A M I L T O N I A N 167 EPILOGUE 168 V LIST OF TABLES 4- 1 Parameters defining the configuration and the initial condition of the hybrid satellite 5- 1 84 Equilibria in SO (3) given in quaternion coordinates, in cr-parameters, as an orthogonal matrix, and by the image of the body frame in the orbital frame 104 vi LIST OF FIGURES 2-1 A schematic diagram of a rigid satellite with two flexible panels. . 2-2 Positions of various reference frames on the hybrid satellite: (a) body 10 frame at c and the principal frame of the rigid section at o; (b) local frames at A\ and A2 12 2- 3 Plate modes: (1,1), (1,2); (2,1), (2,2) 19 3- 1 (a) Complement of B in A, A intersection B, and A union B\ (b) Cartesian product of two sets; (c) an open set in R ; (d) (0,1) x [0,1] 2 is a neighbourhood of (1/2,1/2) € R 2 3-2 but not (1/2,0) 6 R . 2 . . . 30 (a) Schematics of a function: domain, transformation, and range; (b) the graph of a function and its inverse; (c) a function continuous at x\ (d) a smooth map of X at x 3-3 33 (a) Diffeomorphic one and two dimensional spaces constrasted with non-diffeomorphic spaces; (b) a trajectory on a torus 3-4 39 (a) Formation of manifolds by pasting; (b) a covering of a sphere by eight local coordinates 3-5 41 (a) A tangent vector to a manifold; (b) tangent spaces to one and two dimensional manifolds; (c) the construction of the tangent space to a manifold at a given point 3- 6 42 (a) The tangent bundle of an open interval and a circle; (b) a vector field on a two dimensional manifold; (c) an integral curve of a vector field 4- 1 45 A diagrammatic transition from <r-parameters to unit quaternions, in R , to the space of rigid motions with a fixed point, in R . 4 4-2 9 . . . A diagrammatic correspondance between SO(3) and RP(3) defined as a 3-ball with antipodal identification along its boundary 4-3 65 Time evolution of a rigid object seen form the isometric vantage point. Vll 73 Tick marks on the orientation trajectory indicate intervals of one-tenth of an orbit 4-4 74 Time evolution of the orientation of a rigid object seen from the top vantage point 4-5 76 Time evolution of the orientation of a rigid object seen from the side vantage point 4-6 77 Time evolution of the orientation and angular velocity of a rigid object in a circular orbit 4-7 78 Successive orientations of a rigid object, every one-fifth of an orbit, as seen from the orbital frame by an observer along the local vertical looking in the direction of decreasing 03 axis (opposite orbit normal). 4-8 Time evolution of the orientation of a rigid object in terms of successive rotations about body axes bi,b ,b 2 4-9 79 81 3 Time evolution of a rigid object seen from the isometric vantage point represented in cr-parameters 82 4-10 Orientation of the hybrid satellite in the <r-parameters 85 4- 11 Upper and lower panel plate mode amplitudes: (a) mode (1,1); (b) mode (1,2); (c) mode (2,1); (d) mode (2,2) 5- 1 Graph of the nondimensionalized inertia parameter T 3 as a function of Ti and T 5-2 86 102 2 Index of the equilibrium en, equal to the muliplicity of hatching, as a function of T and T : (a) i = l ; (b) i=2; (c) i=3; (d) i=4; (e) i=5; (f) x 2 i=6 5-3 109 Maximum allowable disturbance energy (nondimensional) before tumbling at equilibrium e n as a function of T\ and T 2 5-4 Optimal configuration contours obtained by equialtitude contours of the surface in Figure 5-3 5-5 116 117 (a) Illustration of the threshold of tumbling for a pendulum as a func••• vm tion of disturbance energy; (b) comparison of the energy of equilibria on the energy scale for a given satellite configuration in order to find the minimum perturbation energy 5-6 118 Equipotential contours in the configuration space SO(3) for various intertia parameters: (a) T =• (—0.3,0.4), equilibrium potential = 0.21; (b) f = (0.3,0.3), equilibrium potential = 0.16; (c) f = (0.4,-0.3), equilibrium potential = 0.42; (d) T = (0.3, —0.4), equilibrium potential = 0.46; (e) f = ( - 0 . 3 , - 0 . 3 ) , equilibrium potential = 0.45; (f) f = (-0.4, 0.3), equilibrium potential = 0.28 119 5- 7 The location of the octant o i > 0,o > 0, o > 0 in RP(3). 6- 1 Two classes of satellites for which the equilibrium orientation is stable 2 3 . . . in the linearized analysis 6-2 137 Some possible spectra for linear Hamiltonian systems of (a) two and (b) three degrees of freedom 6-3 125 144 The spectra of the linearized equations of motion on a lattice of points in the parameter space. The spectra in the first quadrant of the complex plane appear to the upper right of each point of the lattice. 6-4 . . . 145 Parameters for which the equilibrium is stable (hatched region) for the linearized equations of attitude motion 146 6-5 Variation of u>i as a function of T\ and T 6-6 Variation of io 6-7 Variation of <JJ% as a function of T\ and T 6-8 The geometry of Thkai's theorem 160 6-9 Resonance regions shown as a collection of curves in the parameter space 162 Status of the equilibrium as a function of T\ and T 163 6-10 148 2 2 a s a function of T\ and T 149 2 150 2 2 ix NOMENCLATURE Operations and functions A fa B A diffeomorphic to B fog composition of functions f and g df derivative of function f at a a / : A —• B a mapping f of set A to B / : a (-> b f maps element a to b /| A function / restricted to the set A VI the canonical 2-form on T * C Abbreviations V for 3 there exists; for some 3! there exists a unique e element of; in => implies; if then A all 2 Ti T ^2 2 T 3 u (1 + m ) 2 —* CLij 2 the homogeneous polynomial of degree i in the formal power H i series expansion of H center of mass cm. Sets < R ,+,*,• > algebra of quaternions Q s unit quaternions c configuration space; C — 5 0 ( 3 ) for a rigid satellite and C = 3 5 0 ( 3 ) x JR for a rigid core with two flexible panels—four modes 8 at each panel TC tangent bundle of C (the phase space of the dynamical system) T*C cotangent bundle of C SO (3) simple orthogonal group of dimension three so(3) Lie algebra of 5 0 ( 3 ) or the tangent space to 5 0 ( 3 ) at identity RP(3) real projective space of dimension three Scalars E E energy potential energy P G gravitational constant m mass of earth e angular speed of the c m . of a satellite in orbit around Earth 1 E kinetic energy k centroidal principal inertias of the rigid portion of a satellite {h,h,h} {ri,r ,r } 2 3 Rodrigues parameters xi {51,52,53} coordinates conjugate to H Hamiltonian {ctxa.2,a.z} sequence of rotations about L Lagrangian {u>i,u>2, W3} {r\,r2,rz} 01,02, and 03 characteristic frequencies of H (first order invariants) {w,y} second order invariants of H m mass of the rigid portion of a satellite m panel (plate-type flexible appendage) mass / panel length, Figure 2-1 h the distance of the lower edge of the panel from the c m . , Figure r p 2-1 b half the panel width, Figure 2-1 D parameter defining the rigidity of the panels p mass per unit area of the panel E s Aijk strain energy i = 1,2 for upper and lower panels, j = 1,2 for longitudinal first and second modes, A; = 1,2 for transverse first and second modes, respectively (3 angular motion about the axis of rotation {o"i,o"2,03} cr-parameters; a preferred local coordinate system on 50(3) e offset of the satellite center of mass w.r.t. the center of mass of the rigid section, Figure 2-2a d m max\d \ Xll Vectors —* 2-2a R c m . of the satellite expressed i n the body frame, Figure d position vector to an arbitrary point of the satellite from its c m . expressed i n the body frame, Figure {oi, o , 03} 2 unit vectors along local normal, local horizontal, and orbit norrespectively (orbital frame), Figure 2-1 mal, —* —* 2-2a —* {61,62^3} the principal body frame of the rigid section at the c m . of the satellite, Figure fi* 2-2a inertial angular velocity of the body expressed in the body frame h,b ,b 2 fi 3 angular velocity of the body frame w.r.t the orbital frame expressed in the body frame T angular velocity of the satellite w.r.t the orbital frame expressed in the orbital frame —* {x, V }, {£, £ } conjugate pairs {z, tu v}conjugate pairs used at various stages of the canonical transformation A ( - ^ l l l j - ^ l l , Au2, A 1 2 , -4l21)-^221i Ai2 , B conjugate pair to A f (r r ) 2 l5 2 A222) T 2 Matrices (bold capital letters) E n identity matrix of order n, abbreviated as E when the order is obvious from the context Il2 111 Inertia operator expressed in body frame, | 7 i / 2 2 ^31 Xlll 2 -^32 Il3 I23 -^33 ACKNOWLEGEMENT I am grateful to Professor V . J . Modi for his unwavering attention, support, and assistance throughout this thesis. I would like to thank Professor D. Rolfsen for many insightful disscussions which directly contributed to the quality of this work. Professor M . Davies, kindly, made himself available to periodically review my progress. Many thanks to all friends at the Department who facilitated the progress of this work by sharing their technical know-how. I would like to express my appreciation to K i l l a m Foundation for their generous support through the pre-doctoral K i l l a m award during the past four years. xiv To my parents, with love X V 1 . INTRODUCTION The transfer of knowlege from one field of specialization to another has often proved inadequate, and this has been particularly troublesome in theflowbetween the various branches of engineering and the underlying branches of mathematics. Certainly I would encourage a more open and receptive attitude toward the movements of thought and people across the boundaries that separate the many disciplines of science and engineering. Peter Likins What was once a scientific curiosity in the fall of 1957 has become an indispensible tool in exploration of earth resources, weather forcasting, and communication. Satellites will be increasingly relied on as the stepping stone for man to wider opportunities. Today, design and development of artificial satellites draws on every aspect of human scientific prowess: from classical mechanics to metallurgy, electronics, computer science, and others. Dynamics and differential equations provide the conceptual tool for the study of the motion of satellites. There is an ever increasing demand for larger satellites with a wide spectrum of requirements. The proposed space platform is perhaps a prime example of an orbiting object with an array of stringent dynamic requirements. This has emphasized a need for more efficient modelling and solution procedures. Perhaps in the science of satellite motion more than in other engineering branches, a need for efficient and accurate analytical and numerical methods of solution is felt as prototype modelling of satellite motion under real life operating conditions is, at present, enormously 1 2 expensive. In fact, it is virtually impossible to simulate environmental effects on dynamics and control of spacecraft during the ground tests. This is particularly true for the next generation of large flexible spacecrafts. To identify the areas of weakness in the current techniques, a sample dynamical problem is analyzed: a satellite with a rigid central body and two symmetrically deployed flexible panels is considered. This problem includes all the salient features of hybrid satellites and is reminiscent of the Canadian Communication Technology Satellite (CTS, Hermes) launched in 1977. The procedures which lead to differential equations describing the motion of the satellite about its center of mass (cm.) are scrutinized to seek improvement. Consideration of motion of a satellite about its c m . is known as librational study or attitude dynamics. Attitude dynamics, in our work, is comprised of the following aspects. (a) A n a priori decision is made about significant features of a satellite for an intermediate simulation which unearths the qualitative librational behaviour of an orbiting object. The choice is made based on the cumulative experience accrued from the past successful simulations; the decisions are empirical in nature. The balance lies between the following two types of models. First, a model which meticulously reflects the details of an orbiting object; such a model is not conducive to mathematical analysis due to its very complex nature. Second, a model which captures only the essential features of a satellite and yet yields some qualitative information which can be used as a guiding rule at the developmental stage of the design. These set of decisions leadis process leads to an abstract dynamical system characterized by a choice of a configuration space and a function (Lagrangian) defined on the tangent bundle of the configuration space. 3 (b) Formulation is the transition from the Lagrangian to the equation of motion. The equations are used to simulate motion of the satellite and predict the vibration response of the system to various disturbances. The method of transition from the Lagragian to the equations of motion, the choice of coordinates which define the orientation of a satellite (librational generalized coordinates), and a method of representing the time evolution of this orientation are some of the salient issues in the formulation process. The method of transition affects the final form of the equations as well as the time and effort involved in the derivation of these equations. The choice of librational angles dictates the degree of faithfulness of the model to the actual system and determines the required C P U time of simulations. Finally, the link between the mathematical model and the designer is established via a judicious representation of results. A perspicuous representation of orientation enhances the ability of the dynamicist to assess the results of a simulation quickly. (c) Analysis of the motion close to an equilibrium yields information concerning stability and its relation to various satellite parameters. This sheds light on the susceptibility of satellite to various perturbations and on ways of reducing this sensitivity by manipulating the free parameters of the satellite. In this study, stability analysis of equilibria of differential equations is the predominant tool. In fact, the most persistent questions of dynamics since the time of Galileo resurfaces in a seemingly simple problem of attitude stability of a rigid object in a circular orbit. The questions on the stability of equilibria of a conservative system, which have been responsible for the creation of the qualitative theory of differential equations, are burning topics of research in the present century. In this study, the entry to the field of attitude dynamics has been through a 4 careful study of the attitude motion of a rigid satellite in a circular orbit. Surprisingly enough, some of the problems of simulation and analysis of stability are not only unresolved but have been, in other forms, the inspiring areas of research and progress in the past three centuries. The unifying theme in this thesis has been the use of qualitative mathematical techniques in the solution of problems of attitude dynamics of satellites. Transfer and translation of techniques from mathematics to engineering has borne many fruits unkown to the engineering discipline as illustrated in the following chapters. The chapters are written as independent units to the extent possible. The number of cross references among chapters is kept to a minimum. Numbering of formulas and references are local to a given chapter. Chapter 2 contains the modelling of two large classes of satellites. The L a grangian of the attitude motion of a rigid satellite is obtained assuming absence of coupling between trajectory motion and the attitude motion. A further set of assumptions on the nature of vibrations of the flexible panels of an hybrid satellite enables one to write a Lagrangian for the more elaborate system of a satellite with two flexible plate-type panels or appendages. In Chapter 3 a method of transition from the Lagrangian to the equations of motion is presented. The method starts with a Lagrangian expressed in direction cosines, angular velocities, and vibrational amplitudes. The resulting equations of motion is found in terms of the same variables without recourse to any librational generalized coordinates. This direct method is first illustrated for the problem of a rigid satellite in a circular orbit and then extended to the case of a satellite with two flexible deployed panels. For the numerical integration of the equations of motion a choice of librational 5 coordinates is necessary. Chapter 4 contains a short review of the known coordinates and the desirable characteristics expected from these coordinates. It is, then, shown how these criteria lead one to a new set of coordinates which for large deviations from the original orientation exhibit the best characteristics for the numerical simulations of the system. The required C P U time for the simulation via the new coordinates and Euler angles is compared both for a rigid and an hybrid satellite. A related topic is the method of presentation of the time evolution of the orientation of a satellite. This has received, virtually, no attention in the literature. In connection with the new coordinates, a method of depicting orientation trajectory is developed which preserves the essential geometric features of relative orientations with ease of visualization. A n analysis of the attitude motion of a rigid satellite is carried out with the help of the energy function in Chapter 5. Bounds of stability in the parameter space for a given level of disturbance is found by comparing the potential energy at various equilibria. This technique owes its inception to Palais-Smale's theorem described in Section 5.3. Use of energy function and an approach to representation of the configuration space yields the bounds of attitude motion for a given level of disturbance. A n application of Morse's lemma (Section 5.3) to the energy function gives a complete local picture of the satellite attitude motion near various equilibria. In particular, a convenient and illuminating proof of the stability of the equilibrium in the Lagrange case is found. Chapter 6 completes the stability analysis of the attitude motion of rigid satellites. In an unresolved case (Delp), the equilibrium of the system for the range of parameters under consideration is stable for the linearized differential equations. However, the Hamiltonian is no longer a Liapunov function. The conclusive proof 6 of stability for this case should be found by examining the higher order terms of the differential equation. Two methods are presented by which the Hamiltonian of the system is brought to a standard form, called the normal form. Then, a stability criterion, based on the coefficients of the normalized Hamiltonian, is applied. This is tantamount to a stability criterion which takes into account higher degree terms of the governing differential equations. The thesis ends with a few suggestions for further research. These suggestions range from topics which though closely linked with the issues covered here require a completely new set of tools for their resolution to those which can be handled through a routine extension of the present techniques. 2. D E V E L O P M E N T OF T H E LAGRANGIAN The most general formulation of the law governing the motion of mechanical systems is the principle of least action, according to which every mechanical system is characterized by a definite function L(q, q, t) and the motion of the system is such that a certain condition is satisfied. L.D. Landau In this chapter, the Lagrangian for two distinct classes of satellites is derived: (a) rigid spacecraft and (b) satellites with two flexible panels connected to a central rigid body. 2.1 A Rigid Satellite in a Circular O r b i t In this section, the potential energy of a rigid satellite in a circular orbit is computed. The Lagrangian is formed as the difference between the kinetic and potential energy. 2.1.1 Any Introduction arbitrary motion of a satellite may be looked upon as a combination of the trajectory motion of its center of mass (cm.) c and the librational or the attitude motion, defined by relative rotations of the principal body frame respect to the right handed orbital frame 01,02,03 (Figure &i,6 ,& 2-1). 2 6*1 3 at c with is along the local vertical, the line connecting the c m . of Earth to c, pointing away from Earth; 7 8 the vector 62 is along the local horizontal, orthogonal to B\ and the orbit normal 03. The sense of 03 agrees with that of the angular momentum vector of a point mass at c about Earth. The center of mass c is assumed to go through a circular trajectory about Earth with an angular rate of 7. The rigid satellite is dynamically characterized by its three principal moments of inertia. 2.1.2 Gravitational potential energy Gravitational potential energy of a rigid object in the gravitational field of a central body of mass m is given by e —* = --^y —* (l-r2o .-+(-) ) ?dm. 2 1 Using Taylor's formula with remainder, one obtains (Appendix 2-1) P E = provided 0 < d /R m Gm m e r ^ 7 ~J < 1. Here, d 3 e 1 T 0 1 ^ ' 1 , 1 „, m^ ~R >' d f R is the maximum distance of a point on the body m from the c m . , 7 = y/Gm /R 3 2 2 is the angular speed of the satellite in a circular orbit as measured by the motion of its c m . , and I is the centroidal inertia operator. 2.1.3 The Lagrangian Dropping the constant terms from the expression of potential energy and preserving only the lowest terms in d /R, one obtains m 37 2 E = - y O i ' Ioip If fi* = (fi£,Cl2,^t) is the inertial angular velocity of the satellite expressed in 9 principal body coordinates (61,62,63), then the Lagrangian of the satellite becomes E -E = -J2l nf- - ^I al , L = 3 1 k p 1 i t=i l i t = i where o = (011,012,013) represents the unit local vertical in body coordinates, - > n* 61,62,63- Dividing the above equation by and renaming as Cl*, one obtains x L= - tl (nf-3al ). 1 y i i i—l Here, |f2*| represents the nondimensional inertial angular speed of the satellite. For example, | f i * | = 1 indicates an angular velocity equal in magnitude to that of the c m . of the satellite. 2.2 A H y b r i d Satellite in aC i r c u l a r O r b i t This section presents the steps leading to the Lagrangian of a satellite with two symmetrically deployed flexible panels. 2.2.1 Satellite configuration A schematic diagram of the satellite under consideration is given in Figure —* —* 2-1. —* In this diagram the body frame (61,62,63) coincides with the orbital frame (0*1,02,03). The satellite is comprised of an arbitrary rigid central body character- ized by its three principal moments of inertia. Two flexible panels of dimensions 26 x / are symmetrically deployed from a station at a distance h from the c m . o of the rigid section. The plates are of flexural rigidity D and of uniformly distributed mass m . The spine (center line) of the plates pass through o. The panel deflections p are assumed to be small; this assumption is more fully explained in Section 2.2.4. 10 F i g u r e 2 - 1 A schematic diagram of a rigid satellite with two flexible panels. 11 The center of mass c of the satellite is taken to be distinct from o. The orbital frame o o , o is placed at c. The body frame bi,b ,63 l 5 2 3 2 —# —* to the principal centroidal frame (X,Y,Z) 2.2.2 of the rigid section of the satellite. Computation of the inertia matrix As in the rigid case, potential energy of the satellite up to the second degree terms in d m is given by ^ P = ^ - O where located at c is defined parallel —* I - I O 1 1 , (1) is the inertia matrix of the satellite expressed in the body frame 61,62563. The following paragraphs elaborate on the derivation of the entries of I . First, the inertia matrix of the rigid body with respect to As shown in Figure 2-2a, frame X,Y,Z 61,62,63 at c is derived. is placed at o along the principal directions. Two parallel frames are placed at the base of the plates P\ and P with 2 their origins at Ai and A , respectively (Figure 2-2b). Let the deflection of plates 2 Pi and P with respect to 77 A if and r)A c coordinate systems be given by 2 2 /i(»7.?) = X ^ i » y ^ ) < M f ) > (2) { J where Auj and A ij 2 are the amplitudes of vibration of the upper and lower plates corresponding to the characteristic function il>i<f>j. respect to the coordinate system X,Y,Z e= 2m The location of point c with is given by — I / d dm p + m J (/ r 2m p -r m + r m Jp i p JPi&cPo r 1 — 2m i p r JB d dm, f +m d dm + / d dm) 2 Jp 1&lP2 12 F i g u r e 2-2 Positions of various reference frames on the hybrid satellite: (a) body frame at c and the principal frame of a rigid section at o; (b) local frames at A i and A 2 . 13 where P i and P refer to appendages and B refers to the rigid core. The c m . offset 2 is given by -^{Mii e = f + A) 2ij 2m„ + m r . . f i/>i{ri)dri ^y(f)df 0 JO J-b Let e be the first coordinate of e (the only nonzero coordinate). Then, the inertia matrix for the rigid core at c expressed i n the body coordinate system 61,62,63 is given by I fl Ii 0 0 \ = I 0 h 0 / o +e m 2 r o 0 1 o 0 l o o i 0 0 / 3 / Next, the moments of inertia of the panels at c expressed in the frame 61,62,63 are computed: J P = 2 K ( A + i ) + ^((26) + / )]; 2 1 1 I12 = ~P Iiz = -P /' I f Jo J-b (A I" a + i f i 2 f2)ridvds; / )(f + - 2 ^2 = P / ' I* [(h - e) + (S + h) ]dr,d< 2 P 2 2 Jo J-b + pf f l(f2-e) Jo J-b fl fb b I / 2 3 = -P 2 + (! + h) }dr,ds; 2 fl rb r){h + $)dnd$-p Izz = pf f l h Jo J-b / Jo J-b ^0 J-b [(fi - ) e 2 + n ]dndi -pf 3 n { - h - $)drjd$ = 0 ; f Jo J-b [(/a - 6 ) + r, ]dr,dt. 2 2 14 O n substitution of mode amplitudes these entries simplify to the following: A i = P 7 7 f + n 2 P P 2 2 / 3 P 2 P 2(/ +^) m ; + = -pJZi uj A 2 i 3 = J ^ 2 i p +A ) / - A j) / 2ij -pJZiAuj 2i / <t>j{$)d$; tpi{ri)dri[h = 2[m(A + - ) + ^/ ] + pj 2 J 2 <j>j{$)dt; + f^y(?)<*f]; [fl + ffidndc - 2(m + mp)e2; r = 0 _ 2m r p b +p f 2 Jo I h J-b (fl + fi)drid$ - 2(m + mp)e2. r The inertia matrix of the satellite expressed in the coordinate system 6\, o , 0 3 2 becomes I = I where I F and I s B + I , P are the inertia matrices of the panels and the rigid body, respec- tively. 2.2.3 Kinetic energy —* —* —t Velocity of any point on the satellite with respect to the body frame (61, b ,63) 2 is given by (—e + / , 0,0), where { fl f 0 2 on the upper panel, on the lower panel, on the main body. Absolute velocity of a generic point expressed in the body frame is v + (n* x d) + v, c 15 where • f-i + f 0 V o V = d= The expression of the kinetic energy takes the form i J{V + n* x d + V) dm = 2 c \ Jft? The + (ft* x d) + V + 2V 2 2 C • (ft* x J) + 2V • V + 2(& C x d) • V]dm. orbit is assumed circular. Having dropped the constant terms from the above relation, one obtains the following expression for the energy, in* 2 • I f2* - -(2m + m )e + i / f 2 2 Jp & p p dm + fi* • [ dx d dm. J 2 2 c r 1 c 2 Here / <f x <Fdm — I dx (—e) dm + I J - dx f dm JP icP l / Jp &p 1 2 dx f dm 2 Jo J-b \ -r,(f 1 + f) 2 J as expressed in the body frame (61,62,63). Therefore, fi* • j dx d = pVl* YtiAiij -pntJZiAiij ~ Mij) f +A ) 2ij / J-b 2.2.4 ^i{v)dri j 4>j{c){h + c)d$ ^i(n)ndn / 4>j{$)d$. JO Discretization procedure Under certain assumptions [2, pp. 254-258], the deformation / of a thin plate is governed by A •! 1 r q 16 Consider a rectangular plate, of dimensions 2b x I lying on rj, f plane, clamped along the n axis as shown in Figure 2-2(a). The associated boundary conditions are %{V,0) = f{v, 0 ) = 0 , -b< <b; M^v,l) -b<n<b; V = i? (77,0 = 0, f M„(±o,<r) =R {±b,s) = 0, TI 0 <<;</. Reaction, i2, and the bending moment, M, are respectively given by the following expressions: ,d f df 2 2 s Now, two denumerable sets of eigenvectors (tpi and ^y) for the beam motion are introduced to represent the plate dynamics. Any deformation of the cantilever plate is approximated by a linear combination of ipi<j>j. The beam differential equation ^ = P f, f'"(Tb) 4 = f'(Tb) = 0 : [—6,6] —>• R corresponding to eigenvalues has a denumerable set of eigenvectors /3f 6 R. The first two modes are given below: * = ^ j , ft = o ; The family (ipi) is orthonormal and total in the prehilbert space C[[—b, b],R] with norm ||V>|| = \Jji ^ dri [3; 4, pp. 295-97]. Furthermore, 2 b f ^/"dr, J —b b = - I" rpU/dT, J-b = -/?y % 2 I i/j/'ip/'dri = %l)iipj""dr) = -/?y <5,y. J-b J-b 4 17 In what follows, ^,'s will govern the transverse variation of the plate deflection. The beam differential equation /(0) = / ' ( 0 ) = 0, ^ T = T7, /'(/)=/'«(/)= 0 has a denumerable set of eigenvectors <j>i : [0,/] —• R corresponding to eigenvalues 1 E R. The first two eigenvectors are given below: A <£i(f) = £ > i ( c o s h ( 7 ! £ ) _ cosC/yxf) + A; (sinh(7if) - s i n ( 7 i f ) ) , 1 1 D = -j=, , k = -0.9825022, t 2.3650204 j ; 71 = t 0 2 (Y) = •£ (cosh(72C) - cos(7 f) + fc (sinh(7 f) - sin(7 f)), , 2 2 2 2 , k = -0.9999664, 1 D = 2 2 5.497039 72 = . 2 ; VI 1 The family (<^,) is a total orthonormal system with following properties: / Jo / Jo 4>i"<t>j"d$ = / Jo ifSij, <f>i'<J>j'dc = - / (t>i"4>jd^ = Jo It follows that ipi{i])(j>j{$) Jo <t>i(f)"d<;. is orthonormal and total in the prehilbert space C[[-b,b] x [0,1], R] with L 2 ipi^j's 4>i4>""d$ = norm | | / | | = yjj l Q $ _ pdnde b h [4, p. 56]. Furthermore, satisfy the geometric boundary conditions as well as the condition of zero bending moment at the free boundaries of the cantilever plate. If the condition of zero reaction at the free boundaries is satisfied then the expansion of a given deflection / in the form Yl o-ij^i^j would converge uniformly and absolutely [4, p. 363], i n addition to the present L 2 convergence. In the absence of any closed form sequence of eigenvectors for the cantilever plate, the characteristic functions ipi(j>j are suitable for discretization of the system [5, pp. 344-345; 6]. That is, deflections 18 /,-, i = 1,2 may be expressed in terms of characteristic functions (Equation (2)). and V 2^ 2- Figure 2-3 show the first four mode shapes tpi(f>i,ipifo,^2^1, 2.2.5 , ) Strain energy of bending and twisting The strain energy of a plate is given by [7, p. 46] which simplifies to D f f l b df d f. 2 f 2 2 j J once is expressed terms,/ one is led to as Ylk.ij kijfpi<f>j- Retaining only Ami, Aku, Ak2i, and Ak22 A D ^,,6.0880682 31.285245 2.6194167 *" ~ J 2J( £4 + 74 + —£272 ) k21 2 A k= l 6.0880682 + ( f4 31.285245 "75 2.2.6 913.60187 + 2 ^fcll ^ 4.0366263 J4 + ) k22 913.60187 Ji 2 A 1.8141853 2 A kl2 ^2j2 k21 k22\- A A The Lagrangian The Lagrangian corresponding to the attitude motion of the satellite may be written as L = E/c — E — E . p (3) s The kinetic energy takes the form E = i n * • i n * + K + n* K + n * x , k x 2 2 3 m o d e (1,1) F i g u r e 2-3 Plate modes: (1,1) and (1,2). F i g u r e 2-3 Plate modes: (2,1) and (2,2). 21 where 2 1 Ki = -^{2m + m )e p + ^ { A 2 r + A 2 + i 2 k l l kl2 2 2 1 + i 2 2 2 ); k=l K = 2 >/pm '[{A p + ( i Kz = U 2 - i )(0.830862/1 + 0.5749084/) 111 211 - A )(0.36377/1 + 0.0191105/)]; 212 - /pm^6[0.3367362(ii + i v 21 ) + 0.1474306(Ai + 2 2 1 22 e = 2m + m, •[0.830862(Aiu + ^ 2 1 1 ) + 0.36377(Au + 2 i 2 )]; 2 2 A )]. 212 p The gravitational potential energy is given by equation ( l ) , where /n=/i + ^(46 Ji 2 hz + / ) + 2(/ +i) m ; 2 2 l p = -& /p7rT -[0.3367362(Ai i + ^ 2 2 1 ) + 0.1474306(,4i + v p 2 = -y/pmp'iiAin + (A 112 hz 2 A )]; 22 222 ^2ii)(0.830862/i + 0.5749084/) - - A )(0.36377h + 0.0191105/)]; 212 = 0; I22 = h + 2m (h + -) 1 p 2 + ml - (2m + A\ + 2 p p + m )e 2 r 2 + P + A \ 1 2 2X A ); 2 k22 k=l hz = h + (2m + p m )e 2 r 2 Ak22)k=l The expressions Ki, i = 1,2,3, are homogeneous polynomials i n A ij. k be used to write L i n the form 1 /n*\ /n*\ 3 2 7 _ This may 22 where the 11 x 11 matrix (hi hi n= hi h3 0 0 V I12 I22 dK aX f 3\T aK aA 3 aK 2 dX 1 —IT Kn 0 —K21 + p 0 K n K 0 0 K 0 0 Kz\ 0 0 0 0 0 K31 0 Kz K32 12 J 0 0 Kn K12 Kn 0 —K22 22 + p 0 K12 Ki 2 K Kiz l2 K12 + P K13 K\2 K13 K13 + P 0 0 0 0 0 0 0 0 0 p 0 0 0 p 0 0 0 0 0 0 0 0 p 0 0 0 0 0 0 0 0 p P m ; p + m + r m ; r = -0.1323286—^—; 2m p + m r = v^mp"(0.830862/i + 0.5749084/); K22 = v^mp"(0.36377/i + 0.0191105/); = -0.3367362^/^/^6; = -0.1474306^^6. \ K32 0 p 3 2 0 0 2m ir 0 0 p 3i 0 0 2m K K31 - K 3 2 0 = -0.3022427—^ 21 0 K31 0 K K 0 0 0 = -0.603317 13 0 0 0 Kn K 0 0 0 0 Here 12 0 0 12 2 0 K 2 0 0 0 dK I33 0 —K i 21 I33 21 0 f 2\T aA ' 0 K 0 K 22 0 0 /13 0 22 h2 dA K —K 0 hz hi o (hi h2 J 23 References [l] Markushevich, A.I., Theory of Functions of a Complex Variable, Chelsea, New York, Vol.1, 1977. [2] Housner, G . W . and Vreeland, T., The Analysis of Stress and Deformation, 3 rd edition, 1975. [3] Taylor, A . E . and Lay, D . C , Introduction to Functional Analysis, Wiley, New York, 1980. [4] Courant, R. and Hilbert, D . , Methods of Mathematical Physics, Interscience Publishers, New York, 1953. [5] Timoshenko, S. and Woinowsky-Krieger, S., Theory of Plates and Shells, McGrawH i l l , New York, 1959. [6] Warburton, G . B . , "The vibration of rectangular plates," Proceedings of the Institute of Mechanical Engineers, Vol. 168, No. 12, 1954, pp. 371-84. [7] Landau, L . D . and Lifshitz, E . M . , Theory of Elasticity, Pergamon Press, Oxford, 2 n d edition, 1970. 24 A P P E N D I X 2-1: EXPANSION OF POTENTIAL ENERGY dm, where x = 2o\ • j| + |jf-| . Application of Taylor's Define A — J ( l + x)~? 2 formula to the inegrand (1 + x) 2" yields - A = J(1 - -x+ ^x + ^(1 + r))-lx ) dm, 2 X 3 where either 0 < rj(x) < x or x < rj(x) < 0. Integrating through second degree terms in the integrand, one obtains A = m + ^ j"(Ml - H°i • <*T) dm + B 2 r = m+ r (rr(I) - 3oi • K i ) + B, where < for J(^) [^(4 3 + If) + ^ a constant K independent of - I f provided + (^ )~^ 2 2 + < 1. Therefore, „ Gm m Gm , ,. . „ 1 „,,d„ E, = — ^ + ^ ( - 1 > ( D + 3ol • lol) + - 0 ( ( - | e r e m 3. A N A L T E R N A T E T R A N S I T I O N F R O M T H E L A G R A N G I A N OF ASATELLITE TO EQUATIONS OF M O T I O N We admit that the present state of the world only depends on the immediate past. Thanks to this postulate, instead of studying directly the whole succession of phemomena, we may confine ourselves to writing down its differential equation. Henri Poincare In this chapter, a procedure is outlined for a conservative system which makes it possible to go from a Lagrangian of a librating system to the corresponding equations of motion in the Eulerian form. The transition does not require a choice of rotational coordinates and makes use of angular velocities and direction cosines directly. The procedure thus synthesizes attractive features of two classical approaches and has far reaching consequences; it is particularly useful in formulating equations of motion for complex flexible systems of contemporary interest. The basic mathematical concepts are briefly touched upon in the beginning which help explain the subsequent development. 3.1 Introduction A Lagrangian L of a rotating conservative system, such as a satellite in orbit or —* a robot arm, has its natural expression in terms of angular velocity 17 and orientation of its body frame expressed in terms of direction cosines (a y)- Equations of motion, t H = /(fit, {aij)), expresssed in terms of angular velocity and direction cosines are desirable for their compactness. The passage from the Lagrangian to the equations 25 26 of motion in the desired form, however, involves an arbitrary choice of coordinates, such as Euler angles, on the space of rotations 5 0 ( 3 ) , i.e. a collection of all 3 x 3 orthogonal matrices of determinant one. This has been true of methods in use at present even though the initial and final products are devoid of such an a priori choice. The choice has been necessary even if one were to retain angular velocities as quasi-coordinates, as in the method developed by Boltzmann and Hamel [l]. Here, a direct transition from the Lagrangian to the equations of motion is achieved by a procedure which owes its inception to the natural geometry of the configuration and phase space [2] as well as techniques of calculus on manifolds [2], [3] and [4], or [5]. The method obviates a priori use of rotational coordinates and provides an efficient transit to the equations of motion thereby reducing the required set of manipulations. The approach makes use of direction cosines throughout the transition. As a result the arbitrary choice of coordinates on 50(3) is eliminated. In a test case, the modelling time of a satellite with rigid central body and two flexible panels has been reduced to one-third of that required by the present methods. First, the method is applied to a simple but nontrivial case of a rigid satellite in a circular orbit. Even at this level, there is no known prior attempt to go from the Lagrangian to the equations of motion using inherent coordinates (i.e., f2, (a»y)). The method consists of the following steps: • the energy function is computed from the Lagrangian; • the Hamiltonian is found; • the vector field on the co-phase space is computed via the canonical 2-form and the differential of the Hamiltonian; • the vector field is expressed as M=f(M,{ )), aij (1) 27 where M is determined from the relation M = 1(6 + 703). To numerically integrate the above equations, a choice of kinematical equations and a relation connecting the angular momenta to angular velocities is adjoined: 'q = ~K(q)fi, n = I- M-~i6 . 1 3 The kinematical equations are not affected by the system under study. Next, the method is extended to the problem of a rigid satellite with two flexible panels. The Lagrangian is of the form L((a,y), f l * , A, A), where A G R 8 refers to the vibrational amplitude of the first four admissible functions associated with each panel treated as a thin plate. A n analogous method to that used for the rigid satellite, yields the following equations of motion: M = f(( ),M,A,B), B = g((a.ij),M,A,B), aij where M is determined from —* —t Here B, the momenta associated with the vibrational amplitudes A, is conjugate to A. These equations adjoined with *-*«>«. ® - n - ' < 4 ( 5 ) - ( f l yield a complete set of first order differential equations. In modelling any conservative system by this method, the Hamiltonian is computed as a by-product. The Hamiltonian remains constant on a given trajectory 28 and so may be used as a check on the numerical integration routine. In the second illustrative example, the accuracy of results is ascertained via the Hamiltonian and also checked against the results obtained for a rigid satellite with identical inertia parameters. 3.2 S u m m a r y of the Mathematical Tools This section summarizes the mathematical tools which appear in the following sections. The reader who is familiar with the concepts may choose to briefly review the notation or skip this section and come back to it as the need arises. A brief verbal explanation in coordination with a visual representation elucidates each concept. For a more formal exposition the reader may refer to any of the following sources: [3] and [4], [2], or [5]. In mathematics a statement is a sentence which is true or false as it stands. Thus 1 < 3 and 4 + 3 = 5 are, respectively, true and false mathematical statements. Many sentences occurring in mathematics contain variables and are therefore not true or false as they stand, but become statements when the variables are given values. Simple examples are x < 5, x + y 2 2 = 10. Such sentences are called statement frames. If p(x) is a frame containing the variable x, then p(4) is the statement obtained by replacing x in p(x) by the numeral 4. One way to obtain a statement from the frame p(x) is to assert that p(x) is always true. This is done by prefixing the phrase 'for every x ' denoted by Vx. (Vx) x < 4 and (Vx, y) (x + y) 2 Then = x + y + 2xy are, respectively, false and 2 2 true. The universal quantifier Vx is also read 'for all x\ Another way to obtain a statement from a frame is to say it is sometimes true. This is done by prefixing the frame with 'for some x', symbolically, 3x. Then, (3x) x < 4 becomes true. The 29 existential quantifier 3x is also read 'there exists an x\ The symbol 3!x reads 'there exists a unique x ' and indicates that there is only one value for x. For example, the statement (3!x) Then, (3!x) x = 1 is false as both +1 and -1 satisfy the equation x = 1. 2 2 x > 0 and x = 1 is true. 2 One writes a G A to mean 'a is an element of A\ or 'a belongs to A\ does not belong to the set A, one writes a ^ A. If a Sets are sometimes designated by displaying the elements in braces; for example, the set of positive even integers less than 10 is denoted by {2,4,6,8}. If A is the collection of all a which satisfy a property p, one indicated by writing A = {a : a satisfies p}. From a given set one can form new sets, called subsets of the given set. For example, the closed interval [0,1] defined as {x : 0 < x < 1} is a subset of real numbers R. As another example, the open interval (0,1) defined as {x : 0 < x < 1} is a subset of [0,1]. One says A is a subset of B and writes A C B if Vsei xe B. Two sets A and B are equal if A C B and B C A. A \ B = {x : x G A and AUB = {x : x G A x G B} defines A union B. AnB or Given two sets A and B, x (fz B} defines the set 'complement of B in A\ = {x : x & A and x G B} defines A intersection B (Figure 3-la). The ordered pair a and b is denoted by (a, b), to be distinguished from the open interval (a, b) from the context. Given two sets A and B, the set of ordered pairs {(a,b) : a G A,b G B} is called the cartesian product of A and B, and is denoted by A x B. For example, i 2 = R x iE is the set of points on a Euclidean plane 2 or is the set of all complex numbers. If S S 1 XR C R 3 1 = {(x,y) G R 2 : x + y 2 2 = 1} then — R x R x i? is the set of all the points on an infinite cylinder of unit 30 1 (d) hWWWWVWWWWVl _ 1 Figure 3-1 (a) Complement of B in A, A intersection B, and A union B; (b) Cartesian product of two sets; (c) an open set in R \ 2 neighbourhood of (1/2,1/2) e R 2 but not (1/2,0) <E R . 2 (d) (0,1) x [0,1] is a 31 radius (Figure 3-lb). (ai, 61) x (<J2,62) is an open rectangle in R . is an open rectangle in R . z VxGX Similarly, (a\,b\) x 2 3{a b ) u («3 63) 5 A set X is open in R if n x • • • x (a , b ) x x (02,62) n x G (fli,&i) x • • • x ( a , 6 ) C X. n n n In other words, X is opeii if for every x in X there exists an open rectangle in R , n ( a i , 61) X • • • x ( a , 6 ) , such that x is in the rectangle and the rectangle is entirely n n contained in X (Figure 3-lc). One may verify that an n-dimensional open rectangle and the n-dimensional open unit ball {x G R w CR n n : \x\ < 1} are open in R . n A set is open in X if 30 open in R n W = O n X. For example, (0, l] = {x G R : 0 < x < 1} is open in [0,1] but is not open in R. A set B C X is a neighbourhood of x in space X if 3 0 open in X x £ O C. B. For example, (0, l ) x [0,1] is a neighbourhood of (1/2,1/2) G R neighbourhood of (1/2,0) G R 2 2 but is not a (Figure 3-ld). A function F is a set of ordered pairs (x, y), no two of which have the same first member. That is, if (x, y) G F and (x, z) G then y — z. It is customary to call y the value of F at x and to write y = .F(x) instead of (x, y) G F. The set of all elements x that occur as first members of pairs (x, y) G F is called the domain of JP, denoted by dora(.F). The set of second members y is called the range of F , and is denoted by ran(F). Let S C dom(F) and ran(F) C T . F ( 5 ) = { F ( x ) : x G 5} is called the image of 5. F is referred to as the mapping from S to T. This is often 32 denoted by writing F : S -> T, s where 5 (->• F(s), i-> shows the action of F on an element s € S. If = T , the mapping is said to be onto T. Given a set B C T, the preimage of B is denoted and defined as f- {B) = {xeS:f{x)eB}. 1 Consider the function of a complex variable defined by the equation F(z) = z . Let 2 a < —. This function maps every sector S = {z € R 2 sector F(S) = {z e R 2 : 0 < arg(z) < a} onto a : 0 < argr(z) < 2a} (Figure 3-2a). The function may be written as F : S -> R z i-> 2 z . 2 If Vx,y (F(x) = F(y) x = y) then F is said to be one to one ( l - l ) . One may readily check that if F is one to one, then {(a, b) : (b, a) £ F} is a function. This function is called the inverse of F. The inverse function is denoted by F . - 1 For example, the inverse of exp : R —> R, x \—y e x is (Figure 3-2b) In : (0,00) —• R, y 1—• /ray. As another example, consider the 1-1 function i i-> a,- which maps {1,2, • • •, n} onto {a,- : i = 1,2, - • • ,n}. («0?=i. This function is denoted by (at)?=i a n d i called a family s F i g u r e 3-2 (a) Schematics of a function: domain, transformation, and range; (b) the graph of a function and its inverse; (c) a function continuous at x; (d) a smooth map of X at x. 34 If two functions F and H satisfy H C F, H is said to be a restriction of F or that F is an extension of H . In particular, if S C dom(F) and if H is defined by the equation Vx € 5 H(x) = F(x) then fl" is called the restriction of F to 5, denoted by F\§. Let G be a set of elements. Assume that each pair of elements x,y 6 G can be combined to yield another element denoted by x • y. The set G with this operation, < G , • >, is called a group if the following axioms are satisfied: • x-(yz) = (x-y)-z; • 3e £ G Vx G G e • x = x • e = x (call e the identity element); • Vx G G 3y G G x • y — y • x = e (call y the inverse of x and denote by x _ 1 ). If in addition Vx, y G G x • y = y • x, then the group is called commutative or abelian. Given two groups < G i , - > and < G , • > a map / : G i — • G 2 2 is a group homomorphism if Vx,yGG! /(x-y) = /(x)./(y) and /(e) is the identity element of G , where e is the identity element of G\. The 2 set of real numbers R with the binary operation + is a commutative group with identity element 0. The set of positive real numbers (0,00) under multiplication is another commutative group with identity element 1. The set of eight elements {±l,±i,±j,±k} form a finite group, with identity element 1, under the following 35 multiplication table: ,* = f — k — -1; ij = k, ji = -k; jk = i, kj = -i; ki = j, ik = 2 -j. This group is not commutative. The set of quaternions Q = {a + a i + a j + a k : a - G R, 0 x 2 3 t i = 0 , . . . , 3} minus the zero element form another non commutative group. The multiplication of letters i, j, k follow the same rules as i n the previous example. Similar to complex numbers, one may associate a magnitude to each quaternion: \a + iai + ja-2 + ka \ = \J'a + a\ + a\ + o | . 2 0 3 The function I|:Q\{0}^(0,oo), Q \q\ is a homomorphism onto multiplicative group of positive real numbers. The set {? € Q : \q\ = 1}, known as the set of unit quaternions, is another noncommutative group. Let V be a set of elements, hereafter sometimes called points or vectors, and denoted by letters u, v, — Assume that each pair of elements u, v can be combined by a process called addition to yield another element w denoted by iw = u + v. Assume also that each real number a and each element u can be combined by a process called scalar multiplication to yield another element v denoted by v = au = a • u. The set V with these two operations, < V,-,+ >, is called a vector space or a linear space if the following axioms are satisfied: • u + v = v + u; 36 • u + (v + w) = (u + v) + w; u + 0 = u; • 30eV VueV • Vu G V 3v G V u + v = 0; (v is denoted — u.) • a(u + v) = au + ctv; • (a + (3)u = au + /3u; • a{fiu) = (a(3)u; • 1 • u = u; • 0 • u = 0. Consider V = R = { ( x . . . , x ) : n l 5 re x - G R} with operations t ( x . . . , x ) + ( y i , . . . , y ) = (x + y i , . . . , x + y ) , l 5 n n x n n a ( x i , : . . , x ) = ( a x i , . . . , ax ). n n W i t h these two operations R becomes a vector space. Given two vector spaces V\ n and V a map <f> : V\ —»• Vj is a n 2 Va,/?Gi2 UjivG^i homomorphism if ^(av + = a<£(u) + f3<j)(w). A homomorphism is an isomorphism if it is 1-1 and onto. vectorspace V is a subspace if A subset W of the is a vector space. For example, the 1-dimensional vector spaces {0} x R and R x {0} are two subspaces of R , respectively, known as 2 the x and y axes. A finite set of vectors v\,..., v G V is linearly independent if n re awi = 0 =>• V i a,- — 0. t=i A finite set of vectors v\,..., v G V is a 6as*s for V if they are linearly independent n and Vv G V 3(a,-)?=i n u = " ' ^ i = i 37 For example V\ — ( l , 0,..., 0), v basis for R . n = (0,1,0,..., 0),..., v 2 = (0, ...,0,1) form a n Sometimes, to distinguish vectors from scalars an arrow is placed on the vectors. A set A with three binary operations -,+,o is an algebra if < A,-,+ > is a vector space and Va G R Va, b E A Va, b,c (£ A a • (a o b) = (a • a) o b = a o (ab), co(a + b)=coa + cob, (a + b)oc — aoc + boc. For example the set of all n x n matrices is an algebra under the operations of scalar multiplication, addition, and the product of two matrices. One can define a binary operation on the space of matrices called the bracket as [A,B] = A B - B A . The vector space of matrices with this product is an algebra with the following properties: [A, [B, C]] + [C, [A, B]] + [B, [C, A]] = 0, [A,B] = - [ B , A ] . A n algebra with these properties is called a Lie algebra. For example, the set of 3 x 3 antisymmetric matrices (i.e.,A = — A ) , denoted by so(3), is a Lie algebra T under the above operations defined for n x n matrices. In fact, / o o o \ Vx = 0 \0 0 -1 1 , 0/ V 2 = /0 0 \1 0 0 0 - l \ 0 , V 0 J / s = o -1 \ 0 form a basis for so(3). Another common Lie algebra is < R ,-,+, 3 i o \ 0 0 0 OJ x >, where • is the multiplication by a scalar, + is the addition of vectors, and x is the cross product of vectors. 38 A mapping of X C R n to R is continuous at x G X if (Figure 3-2c) m ViV neighbourhoods of /(x) A mapping / : X / 1 ( i V ) is a neighbourhood of x in X. Y is continuous on X if for all x G X it is continuous at x. For example, the function sin(l/x), 0, i^O; x= 0 is not continuous at x = 0 but is continuous on (0, oo). A mapping / of an open set U C R to R n m is called smooth if it has continuous partial derivatives of all orders. A mapping / : X —• R m of an arbitrary subset X C i E is called smooth n if each point x £ X is contained in an open set U C R map F : U ^ R n and there exists a smooth n such that F equals f on U D X (Figure 3-2d). A smooth map / : X —• Y of subsets of two Euclidean spaces is a diffeomorphism if it is one to one and onto, and if the inverse / _ 1 : Y —• X is also smooth. X and Y are diffeomorphic - if such a map exists (Figure 3-3a). Let 50(3) = {A : A a 3 x 3 matrix, A A where E 3 T = E , det(A) = 1}, 3 is a 3 x 3 identity matrix. Let O G 5 0 ( 3 ) . Then 50(3) -» SO{3) A H-> O A is a diffeomorphism, mapping £ 0 ( 3 ) C R 9 onto itself. The map is called the left translation by O. The geometric model for the set of all idealized states is the state space. Changes in the actual state of the system are observed, and are represented as a curve in the state space. Each label records the time of the observation. This 39 Protype Diffeomorphic Copies Non-diffeomorphic Spaces F i g u r e 3-3 (a) Diffeomorphic one and two dimensional spaces contrasted with non-diffeomorphic spaces; (b) a trajectory on a torus. 40 is called the trajectory of the model (Figure 3-3b). Many phenomena require geometric models which are simply coordinate spaces. In dynamical systems theory, the geometric models used are manifolds. A n infinite cylinder and a torus are two examples of mainfolds. Manifolds are made of pieces of flat spaces, bent, and glued together (Figure 3-4a). More precisely, a subset M C R is called a smooth n manifold of dimension m if each x G M is contained in W n M diffeomorphic to an open set U of R , m where W is an open set in R . The diffeomorphism n <f> = (4>\,..., (f> ) : W fl M —• U is called a local coordinate system with $ -'s as local n t coordinates. The inverse of </>, (j>~ : U —• W fl M , is called a parametrization of the 1 region W fl M of M. The unit sphere S = {(x,y,z) G R 2 3 : x + y + z = 1} is a 2 2 2 smooth manifold of dimension two. In fact the diffeomorphism {{x,y)eR :x 2 2 + y < 1} -> 2 R, 3 (x, y) H-> (x, y, parametrizes the region z > 0 of S . 2 y/l-x -y ) 2 2 B y interchanging the roles of x,y,z, and changing the sign of the variables, one obtains similar parametrizations of regions x > 0, y > 0, x < 0, y < 0, and z < 0. Since these cover S , it follows that S 2 2 is a smooth manifold (Figure 3-4b). A Lie group is a group and a manifold on which the group operations a H-> a - 1 and (a,b) H-> a • 6 are smooth. A circle, a torus, the set of unit quaternions, and £ 0 ( 3 ) are some examples of Lie groups. Trajectories determine velocity vectors, by differentiation process of calculus (Figure 3-5a). Tangent vectors to curves passing through the same point make up the tangent space at that point (Figure 3-5b). More precisely, the tangent space to a manifold M at x G M, TM , X is defined as follows. Choose a parametrization g:U^MCR k F i g u r e 3-4 (a) Formation of manifolds by pasting; (b) a covering of a sphere by six local coordinates. 42 F i g u r e 3-5 (a) A tangent vector to a manifold; (b) tangent spaces to one and two dimensional manifolds; (c) the construction of the tangent space at a given point to a manifold. 43 of a neighborhood g(U) of x in M, with g(u) — x. Here U is an open subset of R. m Think of g as a mapping from U to R , so that the derivative k dg : R m u as a linear map of R to R , m is defined. Here ( | § r ) u is the k x m matrix of first k partials, evaluated at u. Set T M j equal to the image dg (R ) m u One can show that TM X is an m-dimensional subspace of R k the chosen parametrization. Each element of TM X of dg (Figure 3-5c). u and is independent of is a tangent vector to M at z. Corresponding to the coordinate system <j> = (<j>\,... ,<f> ) = <7 _1 m defined as dg {ei), where (ei)^L is the usual basis for R . m u 1 form a basis for TM . X The family let (d/d(f>i) be x {{d/d(j)i) )^_ x 1 Consider a curve in 5 0 ( 3 ) <f> : R -> 5 0 ( 3 ) , fl 0 !-»• 0 \ 0 0 cos(0) sin(0) 0 - sin(^) cos{6) passing through the identity matrix E £ 5 0 ( 3 ) C R at 6 = 0. One may see that /0 0 0 \ ^/(0) = I 0 0 —1 I £ i 2 is a tangent vector at E . Similarly, one may show that \0 1 o) the curves which represent rotations about y and z axes determine, respectively, 0 0 l\ /o -1 o\ and I 1 0 0 j as the corresponding tangent vectors. Therefore 0 0 0 \ 0 0 0 / - 1 0 0 9 9 form a basis for T 5 0 ( 3 ) E - The tangent space T 5 0 ( 3 ) E is also denoted by so(3). The union of the tangent spaces to a manifold M at each point, (JigA* TM , X has a natural smooth manifold structure, the dimension of which is twice the dimension 44 of M. This manifold is called the tangent bundle of M and is denoted by TM. The tangent space at x is called the fiber of the tangent bundle over x. The map 7T : TM —• M which maps the vector u G TM to the point x at which the vector is tangent (u G TAf^) is called the natural projection. Consider the open interval (0,1) as a manifold. Each tangent space at x G (0,1) is a one dimensional vector space isomorphic to R , denoted by T (0,l). x The tangent bundle T(0,1) is the infinite stripe (0,1) x R with each 'vertical' line at x composing the fiber over x (Figure 3-6a). Note how the union of fibers comprise the tangent bundle. Another example is the tangent bundle of a circle S . It is an infinite cylinder, that is, TS l 1 = S x R 1 (Figure 3-6a). It is possible to extend the definition of a derivative to a smooth map / : M —• AT, with f(x) = y, between two manifolds M C R df : TM x and N C R . The derivative k l TN X y is defined as follows. Since / is smooth there exists an open set W containing x and a smooth map F : W -> R l that coincides with / on WDM. Define df (v) to be equal to dF (v) for all v G TM . x One may prove that dF (v) belongs to TN x y x X and it does not depend on the particular choice of F; this justifies the definition of df . x The prescription of a velocity vector at each point i n the state space is called a velocity vectorfield.A vectorfieldis a set of vectors, one defined at each and every point of the state space (Figure 3-6b); more precisely, it is a smooth function X : M —• TM such that 7r O X = id. A local expression of a vector field is a system of differential equations. A curve in the state space is a trajectory, or integral curve, F i g u r e 3-6 (a) The tangent bundle of an open interval and a circle; (b) a vector field on a two dimensional manifold; (c) an integral curve of a vector field. 46 of the dynamical system if its velocity vector agrees with the vector field at each point along the curve (Figure 3-6c). The local expression of this curve becomes a solution to the system of differential equations corresponding to the vector field. Given a vector space V, one may form the dual space V* = {4>: <f> a linear functional on V } which itself is a vector space with (<t>i + M )) = <t>i{v) + 4> {v), v 2 (o^i)(w) = a<f>i[v), veV; aeR. Here is an example of a basis for so(3)*. O n so(3) define three linear functionals / i > / 2 > / 3 as follows: h : so{3) aiVi + a V 2 + a V 2 3 3 R, i-> a i ; f : so(3) -+ R, 2 aiVi + a V 2 2 / aiVi + a V 2 The linear functionals ( / i ) f (V ) 3 t = 1 m m € M 3 3 3 : ao(3) + a V 3 3 H-> a ; 2 iE, •-»• a . 3 form a basis for so(3)*. This basis is the dual basis to . The dual of a tangent space to a manifold M at m is called the cotangent space T M*. (J =1 2 + a V The cotangent bundle T*M comprises the union of cotangent spaces, T M * . The cotangent bundle of a circle is an infinite cylinder. m A 1-form on a vector space V is a linear functional on V. A 2-form on a vector space V is a function of a pair of vectors <f> : V x V —+ R, which is bilinear and skew 47 symmetric: VA A2 G R Vai,a ,a G V l5 2 3 ^(Axax + A a , a ) = Ax^(ai,a ) + A ^>(a ,a ), 2 2 3 3 2 2 3 <j>(a a ) = -<^(a ,ai). u 2 2 A differential 1-form on a manifold M is an assignment of a 1-form on the tangent space to the manifold at each point. The simplest example of a differential 1-form is the differential of a function. Let / :R -» R, 2 (x,y) !-»• x + y. 2 Then df = 2xdx + dy is a differential 1-form. That is 2xdx + dy : T( )R -* R, 2 x>y (vi, v ) !-)• 2xvi + v . 2 2 There is a way of obtaining a 2-form from two 1-forms: by wedge product. The wedge product of two 1-forms / and g is a 2-form defined as / Ag:R n xR n -* R, {v,w)~det({\% g \%). A differential 2-form on a manifold M is an assignment of a 2-form on the tangent space to M at each point m G M. For example, define a 2-form at (x, y) G R 2 by W( y) = X| T( )R -> R, 2 XtV (v,w)~(x 2 v Then u = (x +y)dxAdy 2 ' v + y)det( Vl ' \Wi is a differential 2-form on R . 2 V 2 ). w J 2 This concludes the summary of the mathematical tools which appear in the following sections. 48 3.3 The Energy Function From the Lagrangian In this section, the notions of fiber derivative and action associated with a given Lagrangian are defined and are used to compute the energy. The procedure is illustrated by the example of a rigid satellite in a circular orbit. Consider the Lagrangian L : TC —• R as a function on the tangent bundle TC of the configuration space C. Classically, the tangent bundle of the configuration space is referred to as the phase space. w h-» dL {w ) € T*C is called the fiber D e f i n i t i o n . The map FL : TC —• TC*, c c c derivative [2, p. 209] of L. L denotes the restriction of L to the fiber over c G C. c T*C is the contangent bundle of C. For example, for the Lagrangian of a rigid satellite in a circular orbit L : SO(3) x so(3) -»• (a,y), n R, 1 H+ - 3 ^ h (^» + 2 2n -a <: + t 3 a\i - 3a ) 2 u i=i the fiber derivative is found by differentiating with respect to the velocity terms, in this case the angular velocity terms: FL : SO{3) x so{3) -»• ,50(3) x so(3)*, ( a y ) , f 2 »-> t Here so(3) « R (aij),(Ii(ni + a \),h{^2 3 + 032)^3(03 + o 3 3 )). is the Lie algebra of the Lie group 5 0 ( 3 ) , i.e., the tangent space 3 to the group at identity provided with the commutator bracket operation ([A, B] = AB — BA). The cotangent space so(3)* is the dual algebra. Let / o o o A Vi= 0 \0 0 1 -1 , 0 J / 0 0 l \ V 2 = 0 \-l 0 0 0 Land OJ /O V 3 = 1 \0 -1 o\ 0 0 0 . OJ 49 Then, { V i , V 2 , V 3 } is a basis for so(3). The map so(3) -»• R 3 , niVi + n v + n v ^ n 2 2 3 3 is an isomorphism and is taken as a local coordinate on so(3). Left translations of this basis yields a frame field on 50(3) which is denoted with the same symbols [4, p. 319]. This frame field is also called the left-invariant frame field generated by ( V - ) ? . Denote the set of three 1-forms dual to { V , V , V } by { / i , / , / } . The t =1 x family ( / ; ) 3 = 1 2 3 2 3 is called the set of left-invariant differential 1-forms on 50(3) dual to ( V i ) ? . These 1-forms at E € SO(3) b ecome a basis for so(3)*. Henceforth = 1 elements of so(3) and so(3)* are denoted by their components with respect to the above bases. The components are referred to as body coordinates. The action A [2, p. 213] is defined as : T C A —• w c R, »-»• FL(w ) • w. c c For our example, A : S O (3) x so(3) R 3 3 (a -y),n i-> X)( *" t'=l n t + °sO f = n + f03i)n t=l Finally, the energy is given by E = A — L [2, p. 213]. For a rigid satellite in a circular orbit, the energy is found to be E : SO{3) x so{3) -+ R, 1 (o ), » vy 3 A-(n? - al + 3a ). 2 u 50 3.4 Hamiltonian From the Energy Function Hamiltonian is defined as the completing member of the following commuting triangle: TC The Hamiltonian corresponding to the Lagrangian L, then, becomes H = Eo FL' : SO[3) x so(3)* -> R, 1 (a y), M t Therefore, H=\ 3.5 "\Y, 7 »'[(^ - 3'") " °3i + a 2 3 a iil- £ = i ( ^ - 2M;a - + 3J -o? .). 3 3t t t T h e Canonical 2-Form and the Differential of the Hamiltonian in the B o d y Coordinates There is a canonical (closed nondegenerate) differential 2-form f2 on T*C for any manifold C [4, p. 202]. The expression of this 2-form on T*SO(3) » SO(3) x —* 5o(3)* in body coordinates at ((a,y),M) G 50(3) x so(3)* is given by [2, p. 315] 3 / Mi n((5,A),(^,0)) = ^ ( 0 . A - A A ) + def ^ i=i V ©i where {6, A ) , ( 0 , 0 ) G T (( M 0 2 02 2 M ^ 3 «5 ©3 J 3 , ^(50(3)xso(3)*) « r . S O ( 3 ) x s o ( 3 ) * and <5, <5, <5 (a y) are the components of 6 G T( ^.)50(3) with respect to frame field { V i , a X V2, 2 V3}. 3 51 Next, the differential 1-form dH is expressed in terms of the basis differential 1-forms on T*SO(3). For 1 1 A/f2 3 2iW 3 d-ff = - y^(—Y^-dMj 2 ,=i — 2azidMi — 2M,rfa ,- + 6 i , a i , d a i , ) , 3 J i or more explicitly d-ff + M a 3 2 + 3(J - / ) a a ) / i + (-M a i + Mia 3 3 + 3(/ =(-M a 2 3 3 3 3 3 +(-Mia 3 2 - -j- M a i + 3(Ji 2 3 — a )rfM 2 - 3 a )rfM 3 3 0 / Here, the relation d(a y) = (a-y) t - / 3 V ~h terms of independent 1-forms dMi,dM ,dM , 2 t 3.6 1 2 1 3 - ii)anai )/ 2 - 7 )anai )/ 3 3 2 2 3 / at each point ( ( a y ) , M ) <E 3 3 a i)dMi 32 t 2 s / 3 2 0 — fi fi 0 fi,fii I is used to express dH in fz spanning T* -T*SO(3) T*SO(3). T h e Expression of the Hamiltonian Vector Field Given a function H : T*C —• R, there is an associated vector field XH on T*C. The local expression of XH in canonical coordinates is classically referred to as the Hamiltonian differential equation. In this section, this vector field is expressed in body coordinates. 52 The vector field XJJ is defined by the following condition: V(0,6) 3!(«5,A) X H n{{6,A),{6,G)) = dH{6,e), : T*50(3) —> TT*SO{3), (a -y), M i — • ( a ) , M , (5, A ) . t i y There exists a unique nonsingular matrix f2 such that n((6,A),(e,Q))= ( | ) - n ( ^ - E ^ 3 J , where M = v - M 0 Mi 0 f M Here, f2 = I j , M -M 3 2 M — M i | and E 3 2 3 is an iden- 0 tity matrix of order three. Then, the above condition implies f2 (jj) = dH, where s , , ^ , dj dH = £ = i ( C / / i + G^AMi) is considered as a column vector ). Therefore, 3 : f t M' i) = n- dH. x As only A is of interest, A = (-B ,M)dH, 3 MM A i = ( —— 2 i.e., 3 3ai ai )(J - J ), 2 3 2 3 J -«3 2 A 2 A 3 = (^^--3ai au)(/ -/i), 3 .M M = ( - ^ iii X 2 3 - 3 a u a i ) ( / i - 7 ). 2 2 2 ^ :R — + M is an integral curve of a vector field X : M —• T M on a manifold M if = Xo(J>. Let g be a local coordinate for the manifold M . From the following commutative diagram TM — • R 2n TqoXoq- 1 53 one finds {<t>u • • •, <t>n, 0i, • • •, 4>n) = Tqo <}>' = TqoXo<}> = Tq o X o q ~ o qo <f> l = ((f>l,...,<f> ,Xi,...,Xn), n where qo cf> — ( c ^ i , . . . , <j> ) and (Xi,..., X ) are the local expressions of the integral n n curve and the vector field. Therefore, <t>i = Xi{4>i,...,<f> ), i=l,...,n n if and only if <j> is an integral curve of the vector field X. W i t h these preliminaries, the local expression for vector field XJI is expressed as • , M M 2 Mi = ( — — M = 2 M 3.7 3 2 (MlMl. _ ,MiM 3 3ai ai )(h = ( — — lil 2 n 3 - h), 3ai au)(/3 - ^l), 3 W r 3anoi2)(/i - _. h). 2 S u m m a r y of the Procedure In this section, the essential features of the transition from Lagrangian to the equations of motion are outlined. The method is generally applicable to Lagrangians of the form L((a,y), fl, A, A), where the vector A denotes the other generalized coordinates such as, for example, amplitudes of oscillation for elastic modes of the system. The procedure is illustrated by an outline of the computational steps required in derivation of equations of attitude motion for a rigid satellite. The sequence of transitions, illustrated i n the last section by the example of a rigid satellite i n a circular orbit, is applicable to any Lagrangian with dependence 54 on direction cosines and the associated angular velocity. Computationally, all the steps are trivial save the formation of dH from H. The fiber derivative FL of a Lagrangian is defined as the gradient of the L a grangian L with respect to the generalized velocities. For example, for the L a grangian of a rigid satellite in a circular orbit L{{ )),n) aij = Ii{n + 2JW- + a\i - Za ) 2 2 u t=i the fiber derivative of L, the gradient of L with respect to fi, becomes J\L((a -y),n) = ( ), {hiQi + a ),I (n t aij 31 2 + a ),I {Q 32 2 3 3 + a )). 33 The second element on the right hand side is referred to as the angular momentum —* M. Therefore I(n + o M = 3 ). The action A is formed by finding the dot product of velocity terms with the fiber derivative of the Lagrangian. In the above example A((a ),n) = FL{( tJ 3 aij ),n) • n = ^(n? + r w ) . The energy is given by E = A — L. For a rigid satellite i n a circular orbit the energy becomes » = i The Hamiltonian corresponding to the Lagrangian L is found as a composition of inverse of FL and E: H = Eo FL-\{ ),M) aij = \JZU[(^Z i=i U - a) 2 3i - af,- + 3a?,.]. 55 Therefore, 1 H/f 2 3 -r ~ = oE( H t=i z + 2Miasi A 3Iia u)- Next, dif is computed. The terms da-y may be expressed in terms of three linearly t independent differential 1-forms (/») . The differential 1-forms (/») 3 =1 3 =1 are related to daij by the following relationship / d{ ) = (o -y) -f 0 h 0 h -h z / V-/2 t aij 3 \ o ; . For the above example dH = - y^(——-dMj 2 i=i — 2a idMi — 2M,da - + 6audau), Z 3l I { or more explicitly dH = (-M a 2 33 +M a 3 + (-M a i -j- M i a 3 3 32 Mj + (-; H M + (-7 a )dM J 2 + (-= H 2 3 3 3 2 3 + 3(J - i i ) a u a i ) / 3 3 2 + M a i + 3(Ji - J )anai )/3 + (-Mia 2 + 3 ( / - 7 )ai ai )/i 3 2 2 3 2 2 a3i)rfMi 32 2 a )dM . zz 3 Then, it can be shown that (Section 3.6) M= if dH = YJ =i 3 + CM-dMi) M = {-E ,M)dH 3 is considered as a column vector ( ^ ) and 0 -M | M -M 0 M -Mi Mi 0 3 2 3 2 The sequence of transitions is indicated below: FL Y E V H dH expressed i n body coordinates Equations of Motion M- (~E ,M)dH 3 57 3.8 A p p l i c a t i o n of the P r o c e d u r e to aSatellite W i t h Flexible A p p e n d a g e s In this section, the equations of motion for a rigid satellite with two flexible panels are derived following the procedure outlined in the previous section. The starting point is the Lagrangian of the librating system which is derived in 2.2.5 (equation 3). The fiber derivative becomes FL : 50(3) x R x so(3) x R -> 50(3) x R x so(3)* x R *, 8 8 8 8 ( a y ) , A , n , A ^ (a ),A,in* + t K tJ 2 V^ y dA 3 The action takes the form R, A : 50(3) x R x so(3) x R 8 s [aij),A,n,A ^ FL(n,A).(n,A) ^ H-> rl / 0 • in* + n. K \+2K 2 l + n* K 2 2 where dA dK -^-•A - 2 = K, 2 dA dA The energy £ = A-L in • in - y o • io + 3 3 o + n • #|2 K 2 \ +E P + E, S #3 may be written in the form ( 0 ^2 \} + E + E . P A + n*K , 3 58 The Hamiltonian, then, becomes where (§) = n( d H Hence =U) • U ) 2 U) • n + - rr *rn ^ *( r 5 0 ( ) 3 x R 8 ) (1) 2 1 in terms of the 1-forms /i,/2, fz,dA,dM, T m 3 + - 7 o . d I o + dE - -7003 • M + S^doi o f >U) - ^ • rf(n a t e a c h 1 8 dB on T*(50(3) x i Z ) , which form a basis 8 P° i n t e T*(50(3) x ((0.7), fl ). 8 That is, 3 8 dH = Y, fifi C i=i 3 + E *i i=i C d A i + 8 Mi i t'=i C + Y, C dBi. dM t=i Bi (2) Analogous to the case of the rigid satellite / M 0 E V0 3 0 0 0 - E 8 - E 0 0 0 3 0 \ E 0 ' 0J 8 These adjoined with a choice of kinematical differential equations q = K(q)Q and the equation GH-'<)-(7> relating the angular velocities and the time rate of amplitudes to momenta, yield a complete set of differential equations (22 in all). The state of the system at each instant is found by a numerical integration of this set. 59 3-9 T h e Software R E D U C E in Derivation of the Equations of Motion To illustrate the procedure involved in derivation of the equations of motion the leading steps to the value of Bi are presented here. The use of a symbolic manipulator language (e.g., R E D U C E ) seems indispensible. From equation (2), in the previous section, B\ — CA - To find C x , one sets X 6 = A = 0 £ R ,X 3 x = ex = (1,0,0,0,0,0,0,0), and A = 0 £ R in dH ^^{6, 8 { From equation ( l ) , one obtains C On ^ = ' ^ (f) h * rf(E_1)l + 2 1 • + d E ^ = h - the other hand, dhz din dL 22 5/33 where, from Section 3-8, d I\ 2 = K idA\ \ 3 + K \dA \ 2 3 dlis = —K \dAm 2 + K idA ii 2 2 '11 + P Ku K K 0 0 0 0 12 dL 33 dL 22 2dA K dA , + -K^32^122 + 22 12 — K dAn 22 2 Ku Ku + P K\ K\ 0 0 0 0 2 2 32 + 222 K dA , 22 K i K\ K +p K 0 0 0 0 K K\ Ki i2 2 2 13 13 212 2 3 K 13 +p 0 0 0 0 0 0 0 0 P 0 0 0 0 0 0 0 0 P 0 0 0 0 0 0 0 0 P 0 0\ 0 0 0 0 0 0 pj A, A, A). 60 As a result dh2{ex) = 0, dli3{ei) - -K i, dl {e ) fK +p\ Ku = 2 K V K 2 u 22 x 12 12 On simplification 31.28245 Z + 3 2 7 [ K\ ^211 -DA 4 111 [a\ + a ^) - K i a i i O i ] , 2 2 2 3 2 #12 y v where -i> l A = e i di 13 1- —•ft 2 1 — — di 33 22 Now a symbolic manipulation code proves useful • to invert TI symbolically; • to find the partials of I I - 1 with respect I\z,I , 22 and 7 . 3 3 In fact, it is more expedient to work with the adjoint of II, |ITjIT , as the -1 entries of an adjoint are devoid of a denominator. Then a ( n - i ) _ ^ | n | - a r f nami ^ |IX| in which operations 81 12 and ai 12 calculation yields the entries of 2 are performed symbolically. A further numerical ^ • The matrix II and hence I I - 1 and adTl are symmetric. It suffices, then, to compute only the upper triangular entries of these matrices. 61 References [1] Whittaker, E . T . , A Treatise on the Analytical Dynamics of Particles and Rigid Bodies, Cambridge University Press, London, 1964 pp. 41-44. [2] Abraham, R . and Marsden, J . E . , Foundation of Mechanics, The Benjamin/ Cummings Publishing Company, Reading, Massachusetts, 2 nd Edition, 1980, pp. 208-23, 311-16. [3] Arnold, V.I., Ordinary Differential Equations, The M I T Press, Cambrige, Massachusetts, Third Printing, 1981. [4] Arnold, V . I . , Mathematical Methods of Classical Mechanics, Springer-Verlag, New York, 1980. [5] Loomis, L . H . and Sternberg, S., Advanced Calculus, Addison-Wesley, Reading, Massachusetts, 1968. [6] Marandi, S.R. and Modi, V . J . , " A preferred coordinate system and the associated representation in attitude synamics," Acta Astronautica, Vol. 15, No. 11, 1987, pp. 833-43. 4. N E W LIBRATIONAL GENERALIZED (^--PARAMETERS) COORDINATES It is not enough for a theory not to affirm false relations; it must not conceal true relations. Henri Poincare Librational motion when expressed in terms of local coordinates as discussed in this chapter yields a set of rational differential equations devoid of planes of singularity. For a specified level of accuracy in numerical integration, a rational set requires less C P U time than an equivalent transcendental set. To interpret the results of the integration, the time evolution of orientation is represented as a curve in the three dimensional topological space i?P(3). The new local coordinates represent a globally defined nonsingular rational parametrization of space of rotations, suitable for problems of dynamics involving general rotations. 4.1 Introduction Librational dynamics is influenced by a myriad of phenomena, each in need of classification according to its influence on the attitude of an orbiting object. However, independent of the complexity of a satellite and the degree of approximation in modelling such a system, there is almost always a platform whose orientation is of essence. To quantify orientation, consider a frame attached to the platform. Let the components of the unit vectors along the axes of this (body) frame in terms of the orbital frame be written as columns of a 3 x 3 orthogonal matrix (a ) with detertJ minant one. The matrix (a,y) is also referred to as the direction cosine matrix. The 62 63 collection of such matrices 50(3) is precisely the space of all orientation preserving isometries of R with a fixed point. 3 One of the objectives of most librational investigations is to study the time evolution of (a y) as a function of initial conditions and system parameters. The t time evolution of (a-y) is summarized as a second order differential equation [l] on t 50(3). A n immediate issue of interest to theoretical analysts, numerical investigators as well as design engineers is the choice of local coordinates [ l , p. 12], also referred to as generalized coordinates on 5 0 ( 3 ) , which yield a convenient and nonsingular set of differential equations. The Euler angles and their variants [2], twenty four i n all, are widely used in practice. They suffer from the following shortcomings. (a) They involve planes of singularity in their domain of definition. This poses an obstacle to the study of the motion in the proximity of the singularities. The singularities are not an inherent feature of the problem, e.g., of the vector field on the phase space T 5 0 ( 3 ) . (b) They introduce transcendental functions of coordinates in the equations of motion, an aspect which adds significantly to the computer simulation time. (c) The depiction of the time history of orientation becomes unwieldy: without a model, it is often difficult to recognize two "close" orientations by an inspection of the triple of the angles. For example, in the body: 3-1-3 convention (the sequence of rotations are about 63,61,63) the sequence (90°, 6°, 270°) is equivalent to a rotation of 6° about 62. Similarly the sequences (11°, 0°, 355°) and (2°,0°,4°) both define a rotation of 6° about 63. A n alternative to the Euler-like angles is to express a ' s i n terms of some set tJ 64 of coordinates appended with a number of relations [3]. The aiy's and quaternions belong to this category. The alternative, though it partially resolves the limitations mentioned in (a) and (b), introduces a higher cumulative error during a numerical integration compared to representations using a minimal number of variables. This chapter presents a set of coordinates for SO (3) To as an answer to (a) and (b). resolve (c), a method of depicting the time history of orientation is developed using the implicit interconnection of Euler's theorem and the above coordinates. The 4.2 details of the above interconnection are elucidated. Formulation Given a unit quaternion [4] g, the corresponding element of SO (3) 2-1 Lie group covering (i.e, a continuous Lie at 1 6 5 3 is an isomorphism) c: S 3 2 —• SO (3) is given by 2 3 2(gig3-9o? ) 2 The group homomorphism whose derivative 2(gig -g g ) go - °l +l\ ~ Qz 2(g g + g gi) ?o + 9 i - ° 2 - 9 3 2(gig + g go) 2 under the 3 3 2(gig + g g ) \ 2(g g - g gi) . g - q\ - ql + 9 J 0 2 3 0 3 0 2 (1) 2 3 0 columns of the above orthogonal matrix are components of the unit vectors along the axes of the body frame in the orbital frame. Details of the derivation appear in Appendix 4-1. The function c establishes a 2-1 correspondence between the points of S and the space of orientations 5 0 ( 3 ) . 3 The 3-sphere S 3 parametrized by the inverse of the stereographic projection p-'-.R 3 - 5 \{1}, p - ( a ) = (|a 3 J | - l,2a 2o ,2<7 )/(l + \a | ), 2 2 u 2 3 where a = (ai,o ,cr ) and \o | = o\ + o\ + o\. Therefore, 2 2 3 cop- -.R l 3 -* 5 0 ( 3 ) , eop- [a) 1 = {a j) i is now 65 q = ( q , fl, fl fl ) 2 3 2-1 F i g u r e 4-1 A diagrammatic transition from <7-parameters, in R , 3 to unit quaternions, in R , to the space of rigid motions with a fixed point, in R . 4 9 66 is an analytic local diffeomorphism (i.e., diffeomorphic on an open neighbourhood of each point), 1-1 on {x G R : \x\ < 1} (Figure 4-1). Here 3 u = (|a| + l ) ; 2 a n a 1 2 2 = l + 8(o- - |a| )/u; 2 a l3 = 4(2aiO" - a (\a | - l ) ) / u ; 2 3 2 = 4(2<7io- + o-(|a | - l ) ) / u ; 2 3 a i = 4(2<TIO2 a 2 2 + 2 o- (\a 3 2 = l + 8(o- -|a| )/u; 2 2 2 | - l))/u; (2) 2 023 = 4(2a a - o-i(|a | - l))/u; 2 2 3 a i = 4(2o"io- - o-2(|or | - l ) ) / « ; 2 3 3 a 32 = 4(2o- o + oi{\8 | - l ) ) / u ; a 3 3 = l + 8 ( a | - |^| )/a. 2 2 3 2 Note the absence of transcendental functions in the entries of (a y). t Given G 5 ' 0 ( 3 ) \ { E } there are precisely two elements (corresponding to ) antipodal elements of S ) which transform to (a y). The origin 0 G R 3 3 t is the only element which is transformed to E . Using (2), one obtains (Appendix 4-II) = B \a ) 2 3 3 2 j, (3) \n J 3 where — 4 B in terms of s is I - al - o\ + 1 a 2(o"io- - o ) 2(o-ia + cr ) ( 3 2(o-icr + 03) -o\ + o\ - o\ + 1 2(-ai + o-o-) 2 2 3 2(o-io- - cr ) 2(o-i + o-o-) -o- - a| + o-| + 1 3 2 2 3 2 3 (B) = > 1. Hence, the proposed parametrization has no singularity (a point 67 at which the linearization of the parametrization is not of full rank), a fact which also follows from local diffeomorphism of c o p - 1 . Let x denote an n-vector corresponding to the other generalized coordinates such as degrees of freedom for an appendage or amplitudes corresponding to modes of vibration. The following set arid relation (3) yield a complete set of equations for attitude motion of an orbiting object: x = v, v € R; n (j^=h(x,v,n,(a ),t), (4) heR . n+3 ij —* Here h is an explicit function of time if some of the appendages describe a predetermined motion with respect to the satellite or if the environmental forces are modeled in a time dependent fashion. 4.3 Interpretation, Properties, and Generalization of the Let q E S 3 ^-Parameters be written as (cos(/?/2), n i sin(/?/2), n sin(/?/2), n sin(/?/2)). The 2 3 orientation q is obtained by a rotation of /? about the unit vector n = ( n i , rc ,ra ). 2 3 The application of stereographic projection, p:S \{l}->R , 3 3 q) P{ = ^^i, to the components of q yields: = *isin(|)/(l= n sin(|)/(l 2 0-3 = n sin(|)/(l3 cos (2-)) = n x c o t ( - ) ; = n cot(^); 4 2 — n3 cotl—). Hence, the orientation given by a is obtained through a rotation of 4 cot \B\ about 68 a. This is equivalent to a rotation of 4 tan \a | about — a . Given two successive 1 rotations E and a, the composition is given by |a+E| 2 ((1 - lE| )a+ (1 - \o\ )t -23 2 or equivalently 2 xt) (|a| |£| -2£.a+l) 2 2 Note that a and —a/\a\ represent the same orientations as they are mapped to 2 antipodal elements of unit quaternions which in turn transform to a single element of 5 0 ( 3 ) . This confines the attitude calculations to \B \ < 1 and avoids large coordinates for a-parameters. This is effected by replacing any a-parameter outside the unit sphere \a\ — 1 with its equivalent element in the interior. There is a six dimensional family of parametrizations equivalent (in properties) to ^-parameters. These are obtained by a projection from an arbitrary point of S . 3 Let M be an orthogonal matrix of order 4. Then poMop- :^ 1 3 50(3), poMop- ^) 1 is another oMike parametrization. For example, if / - l M V 0 0 0 0 0 0\ 1 0 0 0 1 0 0 0 lj = (a,y) 69 one obtains the following: = l-8(|a| -a )/«; 2 2 «12 = 4(2o a + < 7 ( | a | - l ) ) / u ; ^13 = 4(2cT o- -o (|tT| -l))/u; a 2 1 = 4(2(7^2 - <7 (|a | - l ) ) / u ; a 2 2 = 2 1 2 3 2 1 3 2 2 «23 3 l-8(|a| -a )/u; 2 2 = 4(2a a + o - ( | a | - l ) ) / u ; 2 2 3 1 = 4 ( a ( | a | - l ) + 2 - - )/u; 2 2 a-32 0 = 4(^(1-|a | ) + 0-33 = 2 l0 3 2a a )/u; 2 3 l-8(|o| -a )/«. 2 2 Though all these parametrizations are nondegenerate, algebraic, and minimal in the number of components only two, corresponding to /l 0 0 0\ 1 0 0 , projection from the North pole, 0 1 0 Vo 0 0 lj (-\ 0 0 0 0 1 0 0 0 0 1 0 0 0 M and M Vo o o i yield a convenient geometric interpretation. This follows from the computation of p o M ( g ) for an arbitrary unit quaternion a = { {^)^ i {^)^ 2sin{-),n3sin{-)). cos n sin n In the latter parametrization, o represents a rotation of 4 tan 1 \a\ about a. This parametrization was also discovered, independently, by Milenkovic [5]. 70 4.4 T w o Illustrative Examples To assess relative merit of a-parameters, the Euler angles (cti, a ,ct ; 2 z Equation 5b), and quaternions (qo, qi, g , qz\ Equation 5c) a typical example of a rigid satellite 2 in a circular orbit is considered. The governing equations of attitude motion are as follows: fi* = (fi*n*: _ "& \t \3)T\\ a a ft = (ft^ft* °i3aii)r; (5) _3 2 0 3 = 2 (ft*ft2 _ 3 a n i2)l3; a (5a) ° with fi = fi* — I a a 31 ) , where B was given in Equation (3). The nondimensional 3 2 33. intertia parameters T i , T , 2 T z are denned as 3 yJ , ^J^ > repectively. l3 ^J^ , 1 I Equa- tion (5a) relates the time rate of a-parameters with the components of angular velocity fi. Moving to the Euler angles, the corresponding kinematical differential equations have the form co /co si -co si /co z —si /co co si si /co 2 z z z where co - = cos(o; ) and t 2 3 2 = sin(ai), t 2 z 2 2 o\ 0 j 1J fi, i = 1,2,3. The direction cosine matrix (aiy) is expressed i n terms of the Euler angles as [2], co co siisi co + si coi —coxsi co + si sii 2 2 z z 2 z z z (56) —co si —siisi si + si co\ coisi si + co si\ 2 2 2 z z z z z si \ si co I . co\co J 2 x 2 2 71 Kinematical differential equations for quaternions have the following form [3] f -qi /9o\ 9i 92 V93/ -q -q \ 2 3 9o -03 92 03 9o — 9i V-92 fl. (5c) 9o / 9i Ideally, |g| = 1 must always be satisfied. During the numerical integration an error correction routine with |g| = 1 as normalizing condition is used at every step. The energy function, 1 3 E = -J2 t'=i - 2n*a + 3a,-), 2 3i was utilized as an absolute measure of error in comparing the above three parametrizations. As shown in Appendix 4-III, E = 0 and hence E is a constant on trajectories of Equation (5). W i t h initial conditions of 6° tilt about the local horizontal and T = (0,0,-0.05), after five orbits the C P U time corresponding to the three cases as given by Equations (5a), (5b), and (5c) are of ratios 1 : 1.4 : 7.5, for comparable accuracies in numerical integration. The comparable accuracies were achieved by adjusting the local error limit and hence the step size to maintain E within a specified threshold. It appears that the use of cr-parameters promises to reduce computational effort significantly. Applying the same procedure to the example of a hybrid satellite (i.e., a satell i t e with a central rigid body and two plate-like flexible appendages) studied in section 3.7, one finds that use of cr-parameters and Euler angles do not result in significantly different C P U times for comparable accuracies in numerical integration. A significant portion of computational effort in the case of a hybrid satellite goes into panel vibrations, which remains independent of the choice of librational coordinates. It must further be pointed out that the required C P U time for a given run is not only a function of the complexity of the right hand side of a given differential 72 equation but also depends upon the step-size of integration at each time-step. As a result, the time advantage of a given coordinate system may vary both with a given model and the initial condition of the particular run. The absence of singularities and minimality of the number of coordinates for a-parameters are, however, factors which contribute to the ease of operation. 4.5 Representation of Librational Trajectories A concise and effective method of visualization, which clearly portrays the essential features of a given phenomenon, is naturally conducive to further understanding of the underlying principles. The following theorem by Euler is the mainspring of the proposed representation: any orientation in space can be attained by a rotation about some axis. From the discussion on the properties of a-parameters one may deduce that the body frame is realized by rotating the orbital frame about r\t) through \r(t)\ radians, where -4atan0, 1 (|a|)/(|a|7r), a 0; a = 0. In this fashion, the orientation of a satellite at any given instant is represented by a point whose distance from the origin is the amount of rotation in fractions of ir and whose direction is the axis of rotation about which the orbital frame should be rotated to obtain the desired orientation (Figure 4-2). In Figure 4-3 this procedure is applied to the time history of attitude motion of a rigid satellite in a circular orbit subject to initial conditions as used in the illustrative example. Time intervals are indicated on the curve f(t) by tick marks every 1/10 of an orbit. Antipodal points on the boundary of the unit closed ball {x G R 3 : \x\ — 1} represent the same element of 5 0 ( 3 ) ; the space thus obtained, homeomorphic 73 Any rigid motion is a r o t a t i o n about about a fixed some point axis. axis of rotation "amount of" in f r a c t i o n s rotation of TT F i g u r e 4 - 2 A diagramatic correspondence between SO(3) and RP3 defined as a 3-ball with antipodal identification along its boundary. 74 Longitude= 45. f=(-0.1 , 0.3 ) Latitude= 35.27° c={ 0. Position Scale= 10. f=( 0. -0.026, 0. , 0. ) -0.05) 4°3 F i g u r e 4 - 3 Time evolution of a rigid object seen from the isometric vantage point. Tick marks on the trajectory indicate intervals of one-tenth of an orbit. 75 to 5 0 ( 3 ) , is known as RP3 (Real Projective 3-space [6]). Isometric view of the o i , o , o -frame presented here is that of an observer stationed at longitude=45° 2 3 and Latitude=35.27°. Point A with coordinates (0,1,0) denotes a rotation of 7r/10 about the o 2 axis. The position scale is the number by which the lengths are magnified before plotting; for instance, at position scale of 1, point A would denote a rotation of TT about o axis. Similarly, point B with coordinates (0,0.5,0) denotes 2 a rotation of 7r/20 about o 2 axis. Point C , two thirds of the way between the origin and point B , which denotes the initial orientation of the satellite, represents a rotation of 7r/30 radians or 6° about o axis. 2 To enhance the three dimensional aspect of the band r\t), the hidden segments are omitted; the band is linearly thinned the farther it recedes from the observer. Though the representation is primarily for a qualitative assessment of the trajectory, one may directly read the data from the curves as seen from any two distinct views (e.g., top and side view, Figures 4-4 and 4-5, respectively). For instance, point D representing the orientation of the satellite at 4/5 of the first orbit has coordinates (—0.065,0.022,0.02) obtained from Figures 4-4 and 4-5. A s a result, orientation D is obtained by a rotation of 13° about the axis given by (—0.65,0.22,0.2). To monitor T, one may adjoin its components, as a triad, to the trajectory f[t). In Figure 4-6, angular velocity triads appear every 4/5 of an orbit. orientation C (the initial orientation) the angular velocity has only a — o 3 At com- ponent of length 1/5. Taking the velocity scale into account it is apparent that r = (0.00,0.00, —0.05) as it appears on the top of the diagram. Similarly, one may note that f = (0.05, -0.04,0.04) at orientation D . In Figure 4-7, successive orientations of a homogeneous rectangular object, representing the satellite in the illustrative example, is drawn every 1 /5 of an orbit 76 Longitude= 0? f=(-0.1 , 0.3 ) Latitude= 90? *=( 0. Position Scale= 10. f=( 0. -0.026, 0. , 0. ) -0.05) F i g u r e 4 - 4 Time evolution of the orientation of a rigid object seen from the top vantage point. 77 Longiiude= 90. f=(-0.1 , 0.3 ) Latitude= 0? o-{ Position Scale= 10. f =( 0. o. 40 ,-0.026, 0. , 0. -0.05) 3 0- F i g u r e 4-5 Time evolution of the orientation of a rigid object seen from the side vantage point. 78 Longitude= 45? f=(-0.1 , 0.3 ) Latitude= 35.27° *=( 0. Position Scale= 10. f =( 0. -0.026, 0. , 0. ) -0.05) Velocity Scale= 4. Figure 4-6 Time evolution of the orientation and angular velocity of a rigid object. 79 &={ 0. f-( 0. -0.026 , 0. , 0. <°i f =(-0.1 , 0.3 ) ) -0.05) o 2 © 2D01 ©o(3) \ 2DD1 ©o(3) © \ 2001 ©o(3) 1 © 2DD1 ©o(3) © 2DD1 ©o(3) © 2001 ©o(3) © 2001 ©o(3) © 2001 ©o(3) © 2001 ©o(3) © 2D01 ©o(3) © i 20D1 ©o(3) F i g u r e 4 - 7 Successive orientations of a rigid object, every one-fifth of an orbit, as seen from the orbital frame by an observer along the local vertical looking in the direction of decreasing o axis (opposite orbit normal). 3 80 as seen by an observer along the local vertical looking downwards, i.e., along the decreasing 03 axis. The sequence is from left to right and from top to bottom. The first subplot illustrates the dimensions of a generic homogeneous satellite with inertia parameters T\, T . The orientations C and D are shown in the second subplot 2 of the first and second row. For comparison, the same motion as in Figures 4-3 to 4-7 is depicted via the triple of angles 0.1,0.2,0.3 (Figure 4-8). Note the difficulty of visualizing the orientation of a satellite as depicted by a sequence of consecutive Euler rotations. A n alternative to the above representation would be to plot the time evolution of the orientation directly in terms of tr-parameters. That is, corresponding to each orientation (a y), we plot the point — a i n the closed unit ball xG R t 3 : \x\ < 1 ( Figure 4-9). The difference between the representations r and a lies in the way the amount of rotation is related to |rj and \o\. In the case of a, the represented orientation is obtained by a rotation of 4 t a n - 1 (|<J|) about — a, whereas for r, \r\ directly gives the amount of rotation i n fractions of TT about r. The advantage of this representation lies i n its cbnformality. In other words, the angle between two intersecting trajectories in this representation is equal to the angle between the angular velocities of the corresponding motions at the time of the common orientation determined by the intersection point of the trajectories (Appendix 4IV). A typical response of the hybrid satellite to an initial disturbance in the form of the flexible panel deflection is represented in Figures 4-10 to 4-14. The parameter defining the configuration of the satellite and the exact initial conditions appear in Table 4-1. Figure 4-10 depicts the orientation of the rigid section of the satellite about the c m . with respect to the orbital frame. A rise i n the deflection of the a(0)=( 0 . ° , 6 . ° , 0 . ° o 61 ) f=(-.l,.3) f =(0.,0.,-05) 1 F i g u r e 4-8 Time evolution of the orientation of a rigid object in terms of successive rotations about the axes 61,62,63. 82 Longitude= 45° f=(-0.1 , 0.3 ) Latitude= 35.27° ° = { 0. Position Scale= 10. f=( 0. •0 -0.026, 0. , 0. ) -0.05) 3 F i g u r e 4 - 9 Time evolution of a rigid object seen from the isometric vantage point represented in the cr-parameters. 83 upper plate from the undeformed position to its deflected form (1,1) (Figure 2-3) at the rate of 1 mm/s, as an initial condition, initiates a libration of the satellite about the oi axis which evolves to about 0.3° (180°/595) tilt of the satellite in various directions (Figure 4-10). This librational motion of the satellite in turn excites the other modes of the upper and lower plates (Figure 4-ll(a) to 4 - l l ( d ) ) . Note, for example, how the (1,1) mode of the lower plate (Figure 4-ll(a)) is excited out of phase with that of the upper plate. Furthermore, a rapid fall in the amplitudes of vibration for modes in the order (1,1), (1,2), (2,1), and (2,2) is observed; this is an indication of a successively lower contribution of the "higher modes" to the deflection of the plates. 84 Satellite Parameters Initial Amplitudes (m) m= A 400 r Kg n = i 0 m = 10 Kg 1= 7.3 m h= 0.7 m b= 0.6 m D= 0.003 N m Ij= 1071.3 Kg m I = 4761.6 Kg m I = 5500 Kg m Orb. Period- 120 min A A A A A A A Initial Conditions of Rigid Core Initial Time Rate of Amplitudes (m/s) a=0 A = 0.001 A = 0 2 2 2 2 n n 2 212 = = = = = = = 121 221 122 2 2 2 0 0 0 0 0 0 0 2 3 f=0 i n 211 A .= 0 I 8 ^221~ ^ 2 2 2 ^ =0 A ^ 1 2 2 U = 0 T a b l e 4-1 Parameters defining the configuration and initial condition of the hybrid satellite. Longitude= 45.° T=(-0.17, 0.76) Latitude^ 35.° Orbits= 2. Position Scale= 595. F i g u r e 4-10 Orientation of the hybrid satellite in the a-parameters Figure 4-11 ( ) a 87 F i g u r e 4 - 1 1 ( b ) Upper and lower panel plate mode amplitudes: mode (1,2). 88 F i g u r e 4 - 1 1 (c) Upper and lower panel plate mode amplitudes: mode (2,1). 90 References [l] Lang, S., Differential Manifolds, Springer-Verlag, New York, 1985, p. 90. [2] Kane, T., Likins, P., and Levinson, D . , Spacecraft Dynamics, McGraw-Hill, New York, 1983, Appendices I and II. [3] Williams, C . J . H . , "Dynamics modelling and formulation techniques for nonrigid spacecraft," Proceedings of the Symposium on the Dynamics and Control of Non-Rigid Spacecraft, E S A SP 117, Frascati, Italy, 24-26 May 1976, pp. 53-70. [4] Brocker, T . and Dieck, T., Representations of Compact Groups, SpringerVerlag, New York, 1985, p. 5. [5] Milenkovic, V . , "Coordinates suitable for angular motion synthesis in robots," Internal Report, Ford Motor Company, Detroit, Michigan, 1984. [6] Arnold, V . I . , Ordinary Differential Equations, M I T , Cambridge, 1981, p. 233. [7] Warner, F . , Differential Manifolds and Lie Groups, Scott-foresman, Glencoe, Illinois, 1972, p. 100. [8] Arnold, V . I . , Mathematical Methods of Classical Mechanics, Springer-Verlag, New York, 1980, Appendix 2. [9] Spivak, M . , A Comprehensive Introduction to Differential or Perish, Inc., Berkeley, 1979, Vol. 2, p. 84. Geometry, Publish 91 A P P E N D I X 4-1: DIRECTION COSINES A N D QUATERNIONS This appendix starts with a construction of the (universal) covering c : S 3 —• £ 0 ( 3 ) . The last paragraph contains a proof of the Euler's theorem established via c. The proof underlines the explicit connection between the entries of an orthogonal matrix and the corresponding axis of rotation. Let Q = {(go, 9 i , qi, 93) : 9; G R, i = 0,1,2,3} be the division ring of quater- nions [4]. Identify R with {q £ Q : qo = 0}. Given q £ 5 , 3 3 C : Q —> Q, C (a) = qaq q q is an orthogonal transformation, as | C ( a ) | = \qaq\ = \q\\a\\q\ = \a\. Here q = g C leaves R C Q invariant. Hence, it leaves the orthogonal (q , — <7i, — qi, — 93). 0 complement R 3 q invariant. Therefore, C restricted to R , c , is an orthogonal 3 q q transformation isotopic to identity (i.e, 3H : [0,1] x R —• R smooth 3 H(0, •) = 3 c , H(l, •) — id, Vi £ [0,1] H(t, •) is a diffeomorphism) as S is pathwise connected 3 g (i.e., Va,6 £ S 3 3<f> : [0,1] — • S smooth 3 <f>(0) = a,<f>(l) = b). Hence c is q orientation preserving and has determinant 1. Alternatively, one may check that CqC* = Id. Define i = (0,1,0,0), j = (0,0,1,0), and k = (0,0,0,1). In the basis i,j,k, the matrix of c appears as q ( 0.1 + 9.1 {aij) = ~ ?2 ~ ?3 2 qsQo) - qq) 2(?!0 + 2 V 2(9^3 0 2(oiff + q qi) \ ( o g - g 9i) ql - q\ - q\ + q% J ( ? i 9 2 - 9o?3) 9o - 9 i + ?2 - 3 ol 2 2{q q + q q{) 2 2 3 It is obtained by computing c (i),c (j),c (k). q q q 0 q 3 50(3), 0 q,p£ Q p Hence c :S 3 Now c {a) = qpa{qp) = qpapq = c c (a), qp 0 2 c(q) = c q • (I ~ 1) 92 is a Lie group homomorphism. Furthermore, from Equation (1-1) it follows that: 2go = +V {an + 0-12 + 033 + 1); 4gi = («32 - a )/g ; 23 0 a )/q ; 4q = (ai3 2 31 0 4q = (a i - a )/qo3 l2 2 Therefore, 6((a,y)) = (q ,9i, b : 50(3) -» g, 0 92,9s) defined with the positive choice of q is smooth (i.e., infinitely differentiate) about 0 Id. Futhermore, cob and bo c are identity functions about Id. In particular dc\ is an isomorphism. Therefore, c is a 2-1 Lie cover of 50(3) [7]. From the definition of c ,c (q — q ) = q(q — qo)q~ = 9 — 9o i-e., g — go G ii! is 3 q q 5 0 an eigenvector of c €E 50(3) with eigenvalue 1. Suppose go 7^ l(c<j 7^ Id). Define q / = (g — go)/y/(l — q ). Complete J to a right-handed orthonormal basis J, J, K of 2 R. z Then one calculates: jJ = (2g - 1)J + 2 cK q = -2q y/{l-q*)J Therefore, in basis /, J, if, by V ( l - * ) = cos(/V2). 0 0 c reads g 2q yj{l-ql)K; Q + (2g* - l)X. (\ 0 0 cos j3 y 0 sin (3 0 — sin /? J , where /? is defined cos /? 93 ANGULAR VELOCITY A N D CT-PARAMETERS A P P E N D I X 4-II: Consider a motion m : R —• 50(3) in the configuration space. 3 (Vi), the components of m in the left-invariant frame field [9; 8, p. 323] a s (^t') =i Interpret i = 1,2,3, generated by 0 1 0 -1 0 0 0 0 |G 0 T SO{3). Id To express the natural basis (d/dcr,) in terms of ( V i ) , (a y) t9(a,y)/c9cT , k = T t fc 1,2,3, are evaluated as follows: 2(aitT 0 2(-cT a X + 02) 3 a) -2(t7iO- 2 tr - 2 0 2(ai -2(ai + tT a ) + a| - al + 1 2 3 2 o\ 3 0 2(-o- a + tr ) 1 2 2(-CTI + a c r ) 2 o\ 3 + a\ - 1 2 3 0 2(tT! - tr o- ) 2 0 -2(trxcT 3 1 2(t7itr — t r ) - a \ - a \ + al + 1 0 - 3 o\ + o\-o\-\ 03) trf+ 2 0 + cr cr ) 2 + 2 tr - 2 ~o\ + t r | + al - 1 2(CTICT + 0-3) -a - 3 0 3 2(tTitT + 3 + tr ) a) 2 0 2 Components of these matrices in the basis (v ) yield the columns of the following t matrix 4 B _ 1 = (-al + ol + al-l u -2(tTiO- + tr ) 2 I 2{-o o x z 3 + a) 2 erf - 2(-o-!a + 2 a\ + t r f - 1 - 2 ( 0 - ! + cr cr ) 2 3 cr ) a 2 -2{a a 3 2(ai x - + CT| - 3 CT CT ) 2 3 a\ - 1 + a 2 ) 94 A P P E N D I X 4-III: CONSERVATION OF j ENERGY 3 2 . t=l 1 3 ^EMnf-^ + 2 . Sa^-ai,.] t=i 33 1 -2 ^ / t ( n 2 t=i + 3a -a3i ). 2 l t Therefore, E = U(tiiZli + 3 a a - - a a ). lt lt 3i (Ill 3i - 1) t'=i o -n Using the identity [8, p. 323] (a,y) = (a,y) I Q 0 3 n 2 from Equation (5) in (III-l), one obtains E = 0. Hi 3 n 2 — H i I , and substituting 0 95 A P P E N D I X 4-IV: C O N F O R M A L I T Y O F cop - First the conformality of p _ 1 l is shown, i.e., the angle between any two intersect- ing curves is equal to the angle between the image of such curves under the function p at the point of intersection. Tangent vectors along the increasing directions of -1 01 02 j o"3 are given as follows: 5 + al + al\ (l-a\ —2o a\ -2a a 2 2 u 3 V x J 2<7i -2a o x 3 2 do u 2 1+ \ 3 aj - a\ + a 2 —2a a 3 2a 2 J 2 \ -2o\a —2a a 3 2 da? u 3 l + a + al-al 2 J 2a 3 Defining g,y = < d/dai,d jdai >, where < , > is the usual dot product, one obtains { ) =-E. gij Therefore, p _ 1 u is conformal [9, p. 334]. following argument outlines. Let Z\,Z ,Z 2 The function c is also conformal as the be the left invariant frame field on S 3 3 generated by orthonormal tangent vectors i,j,k at 1 6 S , where 3 o\ i = 1 0 oJ , j = 0 1 k= Consider c as i2 \{0} — > R , where the space of 3 x 3 matrices is viewed as a 4 9 subset of R . The linearization of c at q is 9 dc: R —• R , 4 9 dc (v) = L v , q 96 where 90 9l -92 -93 92 9i -9o 92 93 9o 9i 93 92 9i 9o 9o -9l 92 -93 -9i -9o 93 92 -92 93 -9o 91 9i 9o 93 92 9o -9i -92 I Observing that Zi, -93^ 93 J i = 1,2,3, at q £ S are expressed as 3 -93 9o V J 91 respectively, one finds LZi 0 0 2(9o92 + 9i93) 2 ( - 9 o 9 i + 9293) 2(9o93 - "9o + 9? - 9i92) -2(9092 + 9i93) 9o + 9? - 92 - 9l 2 ( - 9 o 9 3 + 9i92) LZo = 2 "9o - 9293) 0 0 9o ~ 9 i + 9 2 - 0 -9o + 9 i + 9 2 ql 2v ; 2 2 ( - g 9 2 + 9i93) 0 93 - 2(9o93 + 9 i 9 2 ) 0 = 2t/i; -2(g 9i + 9293) 0 2(g 93 + 9 i 9 2 ) ql + 92 + 93 2 2 92 + 9s 2(?o9i - 0 LZ, 9 - 9 i + ql 9o 2 ( g 9 i + 9293) 0 2(go92 - 0 9i93) = 2v . 3 0 Hence, c is conformal with respect to metrics in which (Zi) and ( V , ) are orthonormal families [9, p. 334]. i — 1,2, with mi(0) = Consider two intersecting trajectories m,- : R —»• R , 3 m2(0). Define y - = c • p • m ; . From conformality of p _ 1 t Ang(rhi(0), 7722(0)) = Ang(y (0),y (0)). x - 1 and c it follows that Here, Ang( , ) is the measure of the angle 2 between the vectors appearing inside the parentheses. Express the two vectors yi (0) and y (0) in terms of the their components with respect to the frame field ( V , ) 2 (Appendix 4-II) as yi(0) = 91V1 + 02V2 + Then, g = (91,92,93) and G = (Gi,G2,G ) 3 O3V3 and y (0) = G 1 V 1 + G 2 V 2 2 + G3V3. denote the angular velocity vectors of 97 the two motions defined by yi and y% at the instant 0 . Hence, Ang(yi(0),2/2(0)) Ang{g, G) due to the isometry ffiVi Thus, A ? i o ( m i ( 0 ) , m 2 ( 0 ) ) = + 0 2 V 2 + gzV Ang{g,G). 3 (t7i,g ,g ). 2 3 = 5. USE OF T H E ENERGY FUNCTION IN DETERMINING BOUNDS OF STABILITY During the motion of a mechanical system the 2s quantities qi, tji i — 1,2,..., s which specify the state of the system vary with time. There exist, however, functions of these quantities whose values remain constant during the motion, and depend only on the initial L.D. In conditions. Landau this chapter, the energy function of a satellite in a circular orbit is used to shed light on various aspects of the satellite attitude motion. A theorem is formulated and proved on the stability of rigid objects in circular orbits: any nonsymmetric rigid object in a circular orbit has a stable orientation with respect to the orbital frame. Energy is used as a Morse function to demarcate the regions of stability in the parameter space for a given level of disturbance. This provides an effective method of finding the range of possible configurations which remain untumbled while librating about a given equilibrium orientation under perturbations below a prescribed threshold. The possible regions of attitude deviation are obtained by equipotential contour surfaces in the configuration space SO(3); this method focuses on the amount and the direction of possible deviations from equilibrium as a function of disturbance energy for a prescribed configuration. 5.1 Introduction An abstract approach [l] to attitude dynamics of an orbiting satellite, has led us to focus on the underlying properties of the configuration space SO(3). It is shown 98 99 that the dynamical system of a rigid body without axis of symmetry in circular orbit possesses 24 equilibria in 5 0 ( 3 ) ; these are classified into six dynamically distinct groups each containing four elements. A nondimensional parameter space for the dynamical system is defined by — 1 < T\ = / 2 ~ f 3 < 1,-1 < T = 2 f 3 ^ * f < 1. Dynamic symmetries about bi,b , and 63 are respectively expressed by T\ = 0, T 2 2 = 0, and T\ + T 2 = 0. Almost all configurations have no axis of symmetry, i.e., the subset of symmetric configurations is of measure zero. T h e o r e m . One of the six dynamically distinct equilibria is stable for any asymmetric rigid object in a circular orbit. In particular, almost all satellites have a stable orientation. Three proofs are presented. They all rely on a Liapunov function [2] for establishing the stability of a given equilibrium state. The proofs differ in the method of globalizing this local result. The above theorem may be derived from results obtained by Delp [2,3]. Our proof of the theorem underlines the global nature of the result. A n energy function E is introduced with E = 0 on trajectories of the differential equation. B y Liapunov's theorem [4], an equilibrium is stable if E has an isolated local minimum. Hence, the theorem is proved if for any set of parameters corresponding to asymmetric satellites one of the 24 equilibria is a strict local minimum of E. The proof (i) is the most rudimentary. The proof (ii) depends on Morse theory and homology of the configuration space SO(S). The computations leading to the proof (iii) give detailed information on the nature of E at equilibria. This information sheds light on the behaviour of the satellite close to the equilibrium and yields results on the maxmium allowed perturbation (given in terms of the amount 100 of imparted energy) to the satellite before it swings away from the equilibrium. The results are further corroborated by equipotential contour surfaces on the configuration space. For a discussion of notions of manifolds, local coordinates, and other mathematical tools which appear in this chapter the reader may refer to Section 3.2. The proof (i) is the most rudimentary. The computations leading to the proof (iii) give detailed information on the nature of E at equilibria. The results are further corroborated by equipotential contour surfaces on the configuration space. 5.2 Formulation Euler's equations yield the following nondimensional attitude equations of mo- tion for a rigid object in a circular orbit: fi* = (HjOg - 3a ai3)r 1 2 l 5 n^ = ( n ^ n * - 3 a a ) r , n 1 3 (i) 2 fig = (n £} — 3 a i i a i ) T 3 , x 2 2 f aA ( o 3 where fi = Cl* — a . Using (d y) = (a,j) I 3 2 \a 33 fl t J -n 3 Hi V~ 2 n n \ 2 3 0 — f)i 0 J [5], one may write equations (1) in terms of Cl: n f2 x = (f2 n + a 2 3 3 2 a 3 3 - 3 a a ) r + a n ( T i + 1) + a f2 (T' - 1), 12 13 x 3 3 2 = ( n n ! + a a i - 3 a a ) r + a n {T 3 2 3 3 f2 = (n n + a 3 x 2 3 1 3 a 3 2 1 3 n - 3a a )r n 1 2 2 3 31 3 32 1 + 1) + a n ! ( r - 1), (2a) 33 2 + a n i ( r + 1) + 32 3 3 2 031^2(^3 - ! ) • The following complementary kinematical equations are derived [l] from a particular choice of generalized coordinates—the choice, however, is incidental to the subject 101 at hand: ^fflj=B^naj, (26) where — 4 B in terms of a is given by o\-o\-ol +\ - 0 2 2(o-io- + a ) 3 2(0-103 - o- ) A 2(o-i + 0-2+0-3) . -al - a\ + o\ + 1J 1(0102 + 0-3) + 0 - 0J + 1 2(-<7i + a a ) 2(o-xO- - 0-3) 2 2 2 2 2 3 The above choice is derived from the following coordinates on SO (3). Here a = (01 o"2 03)5 5 5 \o 1 = al + a\ + o | . The parametrization 2 i2 — • 50(3), 3 a 1 — • (o -y) t is an analytic local diffeomorphism, 1-1 on {x G R : \x\ < 1} with 3 u = (l*| 2 + l) , 2 an = l + 8 ( a - | a | ) / « , 2 a i 2 = 4(2oia 2 - a (\a 2 | z - !))/«, 2 = 4(2o-ia + a ( | a | - !))/«, 2 3 2 a i = 4(2oia /u + 0-3), 2 2 «22 = l + 8(0 -|o| )/u, 2 2 2 ^23 = 4(20203 - 0"i(|0 | <*31 = 4(20103 — o- (|0 | —!))/«, 032 2 a 3 = l + 8(0 3 3 2 - !))/«, 2 2 = 4(2a az/u 2 + 01), -|0| )/u. 2 The system of equations (2) is the local expression of a 2-parameter vector field (on TSO{3)). The inequality \I -h\ 2 < h implies - 1 < T < 1. Similarly, - 1 < T, < 1, x i= 102 F i g u r e 5 - 1 Graph of the nondimensionalized inertia parameter Tz as a function of T\ and T . 2 103 2,3. The nondimensional inertia parameters (T ) satisfy the relation T\ + T + T + t T T Tz x 2 3 = 0; this relation defines a surface shown in Figure 5-1. The parameter 2 t = 1,2. Given T\ and T space is defined by —1 < T, < 1, eters Ii,I ,l3 2 the inertia param- are determined up to multiplication by a constant. Equations (2) 2 yield the following relations, satisfied by the equilibria: 0 = 0, 032^33 033031 a ia 3 provided 50(3). ( »j')i=i, e TiT T 2 3 2 — — 3ai2ai3 = (3a) 0, 3 a a n = 0, (36) 1 3 — 3 o a i 2 = 0, n ^ 0. The equation (3b) has exactly 24 distinct solutions in 3 These are partitioned in to six "dynamically distinct" classes of four y=r In Table 5-1, the body axes lying on the 01,02-plane are drawn next to each equilibrium element (a y) G 5 0 ( 3 ) . To indicate the equilibria in the RP(3) t repre- sentation of 50(3) [1], one of the two corresponding quaternions to a given (a,y) is included. The orbital frame drawn at en indicates the orientation of the observer in the orbital frame. The details of each subtable are given below: • matrix of direction cosines, ». • orientation expressed in quaternions, m orientation expressed in a-parameters,. body axes . • T h e o r e m . In the dynamical system given in local coordinates by equation (2), the equilibria e y are stable for some t provided T\T T3 ^ 0. t 2 104 / 1 0 0\ ( 0 1 0 ) V 0 /l 0 0 \ ( 0 - 1 0 ) Vo 0 -l) i/ 0 / - l 0 0\ ( 0 - 1 0 ) V o »>i| (o.o.o) '°' b \ ] 0 0 \ 0 -1 ] V o -l o ; f b, . ( v / 2 - 1,0,0) J 3 1 Voi / 0 0 1\ ( 0 1 0 ) /O 1\ 3 2 0 f1 0 - 1 0 \1 0 0 / 0 0 - 1 2 (0,-^,f) o y b lb] b V-ioo; 2 2 (0,1,0) /-l ( 0 \o-io) 2 1—- b2 (0,0,1) j 0 ) /1 0 0 \ / l 0 0\ ( 0 0 1 ) i) o k ( 1 /-l 0 0 \ ( 0 1 0 ) V 0 0-1) V b (0,l-v/2,0) /0 1 0 \ (10 0 ) V o 2 * (0,0,^-i)J o) 5L 0 \ 0 ) /O - 1 0\ ( 1 0 0) o - i ; f-1,1,-1) 1 i) Vo o 2 J (0,0,1-v^jbi / 0 -1 0\ ( 0 0 1 ) V-i I(l + i-J+fc) (1,-1,1) f 1 b V o i |t> 3 Li (1,-1,-1) o o; /O 0 - 1 \ ( 1 0 0 ) o; 1 I »(l-,-+J + fc) 3 (-1,1,1) bT"] / 0 0 -1 \ -1 0 0 ) V 0- 1 0 J ^ 2 /O - 1 0 \ ( 0 0 1 ) [\ (1,1,1) J B I / 0 1 0 \ ( 0 0 - 1 ) Voio; b *b2 V-io 2 b 1 Viooy 2 0 ) (0,x/2-l,0) 2 / 0 -1 ( - 1 0 2 /0 1 0 \ ( 0 0 1 ) \10 o; (f ,0,-f )^ V 0 01 ) Voo-i; b wo /O 0 - 1 \ ( 0 1 0 ) - &k 1 2 b!- / 0 0 -1\ ( 0 - 1 0 ) (1-V^,0,0)j ^ 3i 3 (0,-^,-^3 Vo-i o; i ( i - , - - y + ft) b, (-1,-1,1) lb ^ — i b3* 3 Table 5-1 The list of equilibria in 5 0 ( 3 ) given in quaternion coordinates, in cr-parameters, as an orthogonal matrix, and by the image of the body frame in the orbital frame. 105 ^ ( n , (a,ij)) is a first integral for equation (2), i.e., E = 0 on trajectories Proof: of (2) [1]. B y Liapunov's stability theorem, it suffices to show that E has an isolated (equivalently, strict) local minimum at e,y for some i. As X)f=i h^i l s a positive definite quadratic, focus on 2£ ((a -y)) for local minima. The critical points of E p t p are precisely the equilibria (e,y), provided T1T2T3 ^ 0 (Appendix 5-1). E :SO(3) p —* R remains invariant under the action of the subgroup Hence, the behavior of E {eij) . 4 =1 about the equilibria within a given row, in Table 5-1, p is identical (up to a diffeomorphism). Therefore, it suffices to examine the family ( *'i)f=i f ° local minima (strictness is a consequence of (e,y) forming a discrete e r subspace of SO(3)). The proof is concluded via three distinct methods. These methods give progressively more detailed information about the local behavior of E (i) The analytic (in particular, continuous) function E p p at equilibria. has a minimum on the compact (i.e., has an open complement and is bounded as a subset of R ) 9 set SO(3) [7, p. 18]. (ii) The equilibria are nondegenerate (nonsingular second derivative for E at the equilibria) provided T1T2T3 ^ 0 (Appendix 5-II). Hence, E is a Morse function p (i.e., all the equilibria are nonsingular). Morse inequalities [6] applied to this problem yield V{ > rii, —v + v > -n 0 vo - vi + -VQ Here n ; = dimHi(SO(3), x V2 0 Z ) = 1, 2 3 + rai; > n — n\ + n ; 0 + v\ - v + v = -n 2 i = 0,1,2,3; (4) 2 Q + nx - n + n . 2 3 * = 0,1,2,3 [8]. The variable denotes the 106 number of critical points with index i. The index of a nondegenerate critical point is the number of negative eigenvalues of the second derivative of E p at the critical point. This number is independent of the local coordinates [6]. The n,'s reflect some of the topological properties of 5 0 ( 3 ) . Equations (4) simplify to the following: (a) Vi>l, (b) - v + (c) (d) 0 i = 0,l,2,3; V l > 0; v - vi + v > 1; 0 (5) 2 - v + vi - v + v = 0. 0 2 z Recall that u,-'s are multiples of 4. Equation 5(a) implies v,- > 4, i = 0,1,2,3. Furthermore, a critical point of index zero is a minimum. (iii) B y examining the second order terms of E about the equilibria (e,-y) (Ap- p pendix 5-II), one determines the regions of T\, T2-space at which a given equilibrium has indices 0,1,2,3 (Figures 5-2a to 5-2f). The index of an equilibrium is denoted by the type of the hatching. For example, absence of hatching in the figure demarcates those (Ti,T2)'s at which the four equilibria (eiy)y =1 have index zero. The union of the regions of index zero, taken over the six dynamically distinct equilibria, is the whole of the parameter space minus the set described by T\T (T\ + T ) = 0. 2 2 Index of a function at an equilibrium gives information on the behaviour of the function close to the equilibrium. For example, at the equilibrium e n (Figure 2a) the energy function of a satellite with inertia parameters in the upper half of the second quadrant increases as one moves away from the equilibrium in the phase space. In other words, any perturbation in orientation or angular velocity away from the equilibrium will increase the corresponding energy function. When the index of the equilibrium is one (e.g., when the inertia parameters of the satellite 107 are in the first quadrant or the bottom of the second quadrant) the local behaviour of the energy function is different: there is one direction in the phase space along which the energy falls in magnitude as one recedes from the equilibrium. B y the same token, an energy function with index two or three at an equilibrium indicates the existance of two or three coordinate directions in the phase space along which the energy falls as one recedes form the equilibrium. More precisely, E has index n at a nondegenerate equilibrium if there exists a local coordinate system j / i , . • • ,j/6 about the equilibrium in which the expression of the energy function is n 6 i=l i=n+l where E is the energy at the equilibrium. 0 The variation of the index of the energy function at the equilibrium e i as 2 a function of the inertia parameters is given in Figure 2b. The significance of indices 0,1,2,3 are as described above save that now they refer to the behaviour of the energy function near the equilibrium e i . Moving clockwise about the point 2 T\ = T = 0 in the parameter space, one sees that the sequence of transitions is the 2 same as those of the equilibrium e n : the index varies in increments of one. A t the boundries of these six regions in the parameter space the index is undefined as the equilibrium is degenerate [E has a singular second derivative). The variation of the index is shown for one equilibrium from each of the six classes (Figures 2a-f); this will be identical for the equilibria within each class (rows of the Table 1). Though a function on a manifold of dimension six at a nondegenerate equilibrium can have an index varying between zero and six, the index of E on the phase space 50(3) x R 3 at the equilibria does not exceed three. This is because the energy function is a positive definite quadratic form in angular velocities; any 108 deviation from zero in angular velocities will result in a rise in the energy. Remarks. • K,t>i,i>2,v ) = (4,8,8,4) for all [T T ) 3 U with T i ^ T i + T ) ^ 0. 2 2 • The local maxima and minima are in fact the maxima and minima. • The distribution of indices for equilibria (e,i)®_ may be derived from that 2 of the equilibrium e n by a permutation of axes and a transformation. The procedure is detailed in Appendix 5-II. 5.3 Disturbance Energy and Attitude Perturbation B y Morse's lemma [6], about a nondegenerate equilibrium defined by n = 0, e 6 {e y}, t E is of the form n Eo 6 + Y,Vi- J2 i=l t'=n+l Vi> E 0 = E {e) p = E(0,e), 3 < n < 6 for some local coordinate (y ), where 6 — n is the index of the equilibrium e. t Let e — en,m — min <i<Q\E (en) manifolds 2 ? ( / i ) , E B y Palais-Smale theorem [7], 1 0 0 0 is diffeomorphic to i £ ( / i ) x (Eo,E - 1 0 1 p < h < E + m, are diffeomorphic and E~ ((E ,E -1 E~ ([Eo, — E (e)\. p 2 0 + m)) + m). Consider a connected component B of EQ + m)) containing e. If e is a stable equilibrium, then there exists a diffeomorphism B which maps E\g (E {x e R + h) onto {x £ R 1 0 6 6 : \x\ < m} (6) : \x\ = h} for any 0 < h < m. In Figure 5-3, m* = m/(Ii + I + l3) is plotted versus the pair ( T i , T ). Figure 5-4 depicts the 2 2 equialtitude contour lines of the above surface. The contour lines demarcate those F i g u r e 5-2(a) Index of the equilibrium e n , equal to the multiplicity of hatching, as function of T\ and T%. 110 F i g u r e 5 - 2 ( b ) Index of the equilibrium e i, 2 ing, as function of 7 \ and T . 2 equal to the multiplicity of hatch- Ill F i g u r e 5-2(c) Index of the equilibrium e i , equal to the multiplicity of hatch3 ing, as function of Ti and T%. 112 F i g u r e 5 - 2 ( d ) Index of the equilibrium e+i, equal to the multiplicity of hatching, as function of T\ and T. 2 113 F i g u r e 5-2(e) Index of the equilibrium esi, equal to the multiplicity of hatching, as function of T\ and T . 2 114 Figure 5-2 (f) Index of the equilibrium e&i, equal to the multiplicity of hatching, as function of T\ and T. 2 115 configurations which do not "tumble" about e n under a given perturbation energy h. Absence of tumbling about an equilibrium under a perturbation energy is to be interpreted as the presence of spherical ( 5 ) constant energy levels shrinking to the s equilibrium as h > 0 approaches zero (equation 6). A lower dimensional analogy, is the family of concentric topological circles (S ) 1 about the equilibrium point in the phase plane defined for a simple pendulum (Figure 5-5). More precisely, in the n-dimensional phase space of a Hamiltonian vector field a state 5 is tumbled about a stable equilibrium e if a deleted neighbourhood of e is diffeomorphic to 5 such that for all t G (0,1] the image of S n _ 1 n _ 1 x (0,1] x t under the diffeomorphism is a connected component of a level set of a first integral (e.g., energy function) of the Hamiltonian system and s belongs to the image of S ' n _ 1 x 1. This definition imitates the notion of a point belonging to the basin of a point attractor [9]. The hatched region in Figure 5-4 distinguishes a class of rigid objects which remain untumbled for perturbations of up to 0.02 in energy. This class shrinks as the perturbation increases. The intersection of energy levels in B with the configuration space, represented as RP(3) [1], is given in Figures 5-6a to 5-6f, for a selection of parameters from each of the six regions in the T i , jT -space. What is shown, is the cross-section of 2 equipotential surfaces of E /(Ii p + 7 + 2 I3) with the planes o*i,o ; 0 2 , 0 3 ; 2 seen from the isometric vantage point with respect to the octant of RP(3) 01,03 as defined by 01 > 0,02 > 0, o > 0. Compare with Figure 5-7 for the location of this octant 3 in RP(3). A l l equilibria appearing in the octant are marked in Figure 5-5a save e6i which coincides with e n from the assumed vantage point. Recall that any point Q = (°i) 92? 93) m this representation determines a unique element (a -y) G 5 0 ( 3 ) , t obtained by a rotation of 2 arcsin \q\ about q. For example, a rotation of ir/2 about 62 axis, e , has coordinates (0, V2/2,0) in this representation. 31 116 F i g u r e 5-3 Maximum allowable disturbance energy (nondimensional) before tumbling at equilibrium e n as a function of T\ and T . 2 117 F i g u r e 5-4 Optimal configuration contours obtained by equialtitude contours of the surface in 5-3. 118 F i g u r e 5-5 (a) Illustration of the threshold of tumbling for a pendulum as a function of disturbance energy; (b) comparison of the energy of equilibria for a given satellite configuration to find the minimum perturbation energy. 119 F i g u r e 5-6(a) Equipotential contours in the configuration space 50(3) for inertia parameter T = (—0.3,0.4); equilibrium potential=0.21. 120 F i g u r e 5-6(b) Equipotential contours in the configuration space SO(3) for inertia parameter T — (0.3,0.3); equilibrium potential=0.16. 121 F i g u r e 5-6(c) Equipotential contours in the configuration space SO (3) for inertia T = (0.4, —0.3); equilibrium potential=0.42. 122 F i g u r e 5-6(d) Equipotential contours in the configuration space 50(3) for inertia T = (0.3, —0.4); equilibrium potential=0.46. 123 Figure 5-G(e) Equipotential contours in the configuration space 50(3) for inertia f = (—0.3,-0.3); equilibrium potential=0.45. 124 Figure 5-6(f) Equipotential contours in the configuration space 50(3) for inertia T = (—0.4,0.3); equilibrium potential=0.28. 125 F i g u r e 5-7 The location of the octant oi > 0,o > 0,o > 0 in 2 3 RP(3). 126 These diagrams depict the possible regions of SO(3) within which the attitude motion can occur for a given level of energy. Once a particular shape is selected for which e n is stable (Figure 5-5a), a design engineer can use diagrams of this type to determine maximum amplitudes and directions of the attitude deviation for a given level of disturbance. For example, at an energy level of 0.244, possible orientations are confined to interior (containing en) of the 2-sphere marked 0.244 (Figure 5-6a). This level corresponds to a perturbation of 0.244 — 0.21 = 0.034 in energy. The equipotential contours confirm the information about the index of the energy function. For example, for a satellite with inertia parameters T\ = —0.3 and T 2 — 0.4 the equipotential contours close to the equilibrium e n (Figure 6a) are a sequence of concentric spheres with decreasing energy levels the closer the spheres are to the equilibrium. This confirms the predictions concerning the behaviour of the energy function about the equilibrium when the index of the energy function is zero. For satellite configuration T\ = 0.3 and T = 0.3 (Figure 6b) the energy falls 2 as the satellite rotates about the o\ axis away from the equilibrium orientation—a case of index one. Whereas the energy rises away from the equilibrium in the other five directions in the phase space: rotations about o , 0*3, and deviations in the 2 angular velocity components ( f i i , II2, $^3). When the energy function has index two at the equilibrium (e.g., T\ = 0.4 and T 2 = —0.3, Figure 6c) the energy falls in two directions (rotations about o*i and 02). Finally, when the index of the energy function is three (e.g., T\ = 0.3 and T = —0.4, Figure 6d) the energy falls in all 2 three directions in the configuration space; the equipotential contours close to the equilibrium are again concentric spheres with the difference that the energy falls the farther one recedes from the equilibrium. Continuing clockwise in the parameter space, one reaches the second region of index two. For a configuration in this region (e.g., T\ = —0.3 and T = —0.3, Figure 2e) two directions from the equilibrium lead 2 127 to lower energy levels (rotations about o and 03); note, however, the difference 2 between this situation and that of the other index two configuration (Figure 2c). A n analogous situation exists for the equipotential contours corresponding to configurations from the second index one region (Figure 2f); contrast this with Figure 6b. As a further crosscheck between the information given in Figures 2 and 6, note how the presence of a stable equilibrium affects the contours in the vicinity of en—see Figure 6a for the postion of equilibria which appear on the walls of the first octant. For example, for the configuration T\ = 0.3 and T = 0.3 Figure 2 2b indicates that the energy function has a minimum at equilibria e 2 i , e 2 2 , ^23, «24 (second row in Table 1) and hence the potential energy can only fall as one moves away from e n close to e22 (Figure 6b); a result which confirms the predictions made based on the index of E given in the previous paragraph. 128 References [l] Marandi, S.R. and Modi, V . J . , " A preferred coordinate system and the associated orientation representation in attitude dynamics," Acta Astronautica, Vol. 15, No. 11, 1987, pp. 833-843. [2] Debra, D . B . and Delp, R . H . , "Rigid Body Attitude Stability and Natural Frequencies in a Circular Orbit," Journal of the Astronautical Sciences, Vol. 8, N o . l , Spring 1960, pp. 14-17. [3] Kane, T., Likins, P., and Levinson, D . , Spacecraft Dynamics, McGraw-Hill, 1983, pp. 199-210. [4] L a Salle, J . and Lefschetz, S., Stability by Liapunov's Direct Method, Academic Press, 1961, pp. 33-37. [5] Arnold, V.I., Mathematical Methods of Classical Mechanics, Springer-Verlag, 1980, Appendix 2, p. 174. [6] Hirsch, M . W . , Differential Topology, Springer-Verlag, 1976, Chapter 6. [7] Choquet-Bruhat, Y . and Dewitt-Morette, C , Analysis, Manifolds, and Physics, North-Holland Publishing Company, 1982, pp. 568-570. [8] Vick, J . W . , Homology Theory, Academic Press, 1973, p. 85. [9] Abraham, R . H . and Shaw, C D . , Dynamics—The Geometry of Behaviour, Aerial Press, Inc., Santa Cruz, California, 1981, pp. 44-5. [10] Loomis, L . H . and Sternberg, S., Advanced Calculus, Addison-Wesley, Reading, Massachusetts, 1968, pp. 3-4. A P P E N D I X 5-1: T H E C R I T I C A L P O I N T S O FE t Direct calculation with differential forms [5] yields h 0 -fi fi 0 / 0 d(aij) = (dij) f \-h -f 3 3 where (/,) are the left invariant differential 1-forms on 50(3) dual at E G 5 0 ( 3 ) . Now, 3 dEp — 5^ ii{3aiidau — a ida i) 3 3 i=i = {h - ^3)(3ai ai3 2 + (I 3 - 0 203 )/l 3 3 - I i ) ( 3 a a i 3 - 031*133)h u + {h ~ ^ 2 ) ( 3 a n a i 2 - a i a ) / 3 . 3 Therefore, dE = 0 and T T T p X 2 3 ^ 0 imply 3 d i a i 3 = tX3 a33, 2 2 3 a n a i 3 = 03i033> 3anai 2 = 031032. 3 2 130 T H E EXPRESSION OF E AT A P P E N D I X 5-LT: EQUILIBRIA p Expansion of E at critical points e,i up to the second order in local coordinates p a = (o~i,o~2,o~z) yields the following. A t en, i(3/! (16(-7 + 7 ) 0 V 0 2 where A = - 7 ) + \t A.t, 0 64(-7i + 7 ) 0 3 t = a; T 3 0 0 48(-/ + / ) 3 1 2 At e , 21 2 T /2(/ -I ) 0 V 0 2 where A = t = a-{V2- - h) + ^(1 + V2) t At, |(3/i 0 - 7 7 + 47 + 37 3 2 2 -7x + 47 - 37 2 -I 3 1,0,0); 0 + 47 - 37 2 x 3 - 7/X + 472 + 3/3 3 Ate , 3 1 ^{-h + 3/ ) + ^(1 + V2)H At, / h + 27 0 \ h - 47 + 2 where A = t = 3 - (0,1 - y/2,0); T 3 2 37 h ~ 0 8(7! - 7 ) 0 3 37 3 4/ + 3/ 2 0 7 + 27 - 3 7 3 X 2 3 3 At e i, 4 -(37 where A = At / - 7 ) + -t At, 3 - 7 7 i + 27 + 5 7 5(7!- 27 + 7 ) 2 3 - 7 7 i + 27 + 5 7 2 5(7i - 27 + 7 ) 3 2 V(-7i r= a+(—,—,0); T 2 2 + 4 7 - 37 )\/2 2 3 3 (-7i+ 4 7 - 37 )\/2 2 3 3 (-7 + 4 7 - 37 )\/2 \ T 2 /4(7x-7 ) 0 \ 0 2 where A = 0 3(-7 + 7 ) 0 2 0 0 7i-7 3 2{-I t = a + (1,1,1); 3 A t e6i, U-h + 3 7 ) + \t At, T 3 3 2 esi, + 3 7 ) + \t At, 2 X {-h + 47 - 3I )V2 t = a - (1,1,1); x 3 - 4 7 + 57 ) J 2 3 131 '3(Ji-7 ) 0 0 0 - I + 7 0 3 where A = | Index of E p x 0 0 4(7 - 7 ) 2 2 3 at e n is zero if -7 + 7 > 0 & 2 3 h + 7 > 0 & - 7i + 7 > 0; 3 (1) 2 equivalently, Ti < 0 & T > 0 & T < 0. 2 3 Index—1 if (Ti > o & r > o & r < o) 2 3 or (Ti < 0 & T < 0 &: T < 0) or (Ti < 0 & T > 0 & Ts > 0). 2 3 2 Index=2 if (Ti > 0 & T < 0 & T < 0) 2 3 or ( i i > o & r > o & r > 0) or (Ti 2 < 0 3 & r < o & r > 0). 2 3 Index=3 if Ti > o & r < o & r > o. 2 3 Here " & " and "or" are the logical "and" and "or" [10]. A n analogous procedure applies to the other equilibria after diagonalization of the corresponding matrices. A n alternative to the use of expansion of E v at the equilibria (e,i)f_ 2 is the use of symmetries. For example, the regions of various indices in the parameter space for the equilibrium e±\ are obtained by analogy from (1) and use of Table 5-1. Substitute 7 and I\ for 7 and 7 in (1) to find 2 X -7 2 + 7 > 0& 3 2 - h + 7 > 0& - 7 + 3 2 h > 0 132 in order for E p to have index zero at e^\- Equivalently Ti < 0 & T > 0 & T > 0. 2 3 This leads to a partitioning of the parameter space for e4i which is the mirror image of that of e n with respect to diagonal Ti + T = 0 (Figures 5-2a and 5-2d). 2 6. ATTITUDE STABILITY OF RIGID SATELLLITES VIA A NORMALIZED HAMILTONIAN In studying the behaviour of solutions to Hamilton's equation near an equilibrium position it is insufficient to look only at the linearized equation. V.I. Arnold A breakthrough is made in the long outstanding problem of attitude stability of a rigid satellite in a circular orbit. The focus is on an equilibrium at which the linearized techniques are inconclusive and no Liapunov function is known. In this chapter, the Hamiltonian for the attitude motion of the satellite is expanded about the equilibrium up to the terms of degree four. B y a sequence of canonical transformations the expansion is brought to a normal form. A region in the parameter space for which the equilibrium orientation is stable is demarcated by an application of a recent extension of K A M theory to the normal form of Hamiltonians. This ap- proach resolves the question and opens the door for the stability analysis of higher order conservative systems, such as flexible satellites, not amenable to traditional tools. 6.1 Introduction The question of librational stability of satellites has a long history. It was known that the moon always presented its same face. But it was not till the advent of the telescope when the empirical laws governing the libration of the moon were dicovered. Libration refers to the oscillations of the moon about its c m . which during a period of each month causes sides of the moon, as seen from the earth, to 133 134 appear and disappear alternately. Galileo was the first to observe this librational phenomenon; however, it was left to Dominique Cassini to present a general and complete explanation of it. The explanation which was published by Jacque Cassini, the son of Dominique, in 1721 in the Memoirs of Academy of Sciences of Paris supposed that (a) the moon was spinning with a uniform angular rate about an axis inclined at 87.5° to the ecliptic and 82.5° to the orbital plane of the moon, and (b) one revolution about the spin axis was completed in a period of 27 days and 5 hours which was equal to the orbital period of the moon. The surprise of the exact equality of spin period with the orbital period was compounded by the discovery of two other moons in the solar system with the same properties: the first moon of Saturn and the fourth moon of Jupiter. In parallel with the above empirical data, attempts to give a sufficiently detailed model for the motion of the moon to fit the observed phenomena was on the way. Newton thought that the moon was elongated towards the earth and that fact was reponsible for the stable libration of the moon [l]. However, the actual treatment of the motion of a rigid object about its c m . had to wait the immense contributions of Euler. Finally, the librational stability of the moon during its orbital motion became the topic of a prize essay of the Paris Academy of Sciences in 1763. The prize was awarded to Louis Lagrange who with his mathematical model confirmed Newton's conjecture of the nonsphericity of the moon [2]. The work was further corroborated by d'Alembert and Euler. Lagrange's definitive work [3] on the subject appeared in 1780. In "the theory of the libration of the moon and the other phenomena which depend on the non- 135 spherical shape of the moon" he derived the four inequalities ( I - h ~ h) 2 3 + (/ 3 + 4 ( / - h)h - h)h 3 (/ -/ )(/ -/i) 3 [(/ - I ~ h) 2 3 2 + (/s - h)h 2 3 + 4 ( / - h)h\ 2 3 (/ 2 - h) > 0, >0, > 16/2/!(J - / ) ( / - / i ) , 3 2 3 > 0. These inequalities were to be satisfied in order for the moon to assume stable librational motion. The result was obtained after linearization of the equations of librational motion. In 1885, Henri Poincare [4] observed that the linearized analysis of Hamiltonian systems did not yield a conclusive proof of the stability of the original system. In the meantime, a number of researchers approached the problem of the motion of the moon about its centre of gravity, notably, Laplace [5] and Tisserand [6]. A summary of these contributions has been given by Routh [7]. Almost two centuries after Lagrange's work, Beletskii [8] and Delp [9] rederived some of the same results in connection with the rising interest in artificial earth satellites. The linearized analysis indicates that any rigid body in a circular orbit about a homogeneous spherical mass is stable in an orientation with the principal axis of the least inertia along the local vertical and the principal axis of maximum inertia along the orbit normal; a satellite in this stable orientation was said to have the Lagrange configuration [10J. Another class of satellite shapes (Delp satellites), with the principal axis of least inertia normal to the orbital plane and the axis of maximum inertia along the local horizontal, are also stable with respect to the linearized equations (Figure 6-1). Subsequent to the publication of Delp's paper the inadequacy of the linearized analysis as a conclusive proof of stabililty was widely acknowledged [10]. The use of the Hamiltonian as a Liapunov function establishes the stability of the Lagrange satellites. For the Delp satellites, however, the question 136 of stability or instability remained unanswered [10, 11]. A parallel problem which has stimulated much growth and controversy in dynamics since the early days of Newton is the stability of the solar system. In particular, the first intractable special case is the stability of the three body problem. This question was analyzed by Laplace, Lagrange, Poisson, Dirichlet and many others, all of whom claimed to have proved that solar system was stable. As Dirichlet died before writing down his proof, K i n g Oscar of Sweden offered a prize for its discovery, which was given to Poincare in 1889. Although Poincare had not solved the given problem, his prize essay contained an abundance of original ideas which are recorded in his New Methods of Celestial Mechanics [12]. Birkhoff [13], Carl Siegel, and Andrei Kolmogorov further extended these ideas in search of a criterion for stability of equilibria in Hamiltonian systems. Arnold [14] in the early 1960's proved a criterion, for systems with two degrees of freedom, the statement of which was announced by Kolmogorov a few years earlier. Arnold's work and its further refinement by Moser [15] threw some light on the behavior of the Hamiltonian systems close to integrable ones. In 1970, Arnold proved the algebraic unsolvability of the problem of Lyapunov stability in general [16]. This amounts to the nonexistence of an algorithm which leads in a finite number of steps to a decision on stability of an equilibrium in every dynamical system. This opened the possibility that the question of stability of the Delp satellites or the three body problem may have no complete answer. Nonetheless, an algorithmic sufficient condition was not ruled out and was found by Tkhai [17], in 1985, for arbitrary Hamiltonian systems in the normal form. The normalization techniques were implicit in the works of Lindestad and Poincare, and were explicitly developed by Birkhoff [13]. A complete classifica- 137 S A T E L L I T E F i g u r e 6-1 Two classes of satellites for which the equilibrium orientation is stable in the linearized analysis. 138 tion of the normal forms of linear dynamical systems is given by Williamson [18]; however, his work did not yield an algorithmic procedure. Laub and Meyer [19], and Burgoyne and Cushman [20] published an algorithmic solution to the normal form of linear Hamiltonian systems in 1974. For normalization of higher order terms of a Hamiltonian two distinct techniques emerged; one was a variation of Birkhoff's method as presented by Gustavsoh [21] and the other was based on the Lie structure of the space of polynomials. The latter appeared in the work of Deprit and Henrard [22] and is summarized by Churchill, Kummer, and R o d [23]. In this chapter, two distinct procedures are applied to the Hamiltonian of a rigid satellite about an equilibrium, in the Delp case, to normalize the Hamiltonian up to the terms of degree four. The sufficient stability condition of Tkhai is applied to obtain a set of parameters in the Delp region for which the equilibrium under consideration is stable. Markeev and Sokolskii [24] approached the Delp case with a similar technique—their choice of normalization procedure is not clear from the paper—but their research did not lead to a stability condition as Tkhai's sufficient criterion was not discovered till ten years later. 6.2 T h e Hamiltonian and the Choice of Generalized Coordinates The Rodrigues parameters r = ( r i , r , r ) provide a local coordinate system 2 3 about E G 5 0 ( 3 ) . The corresponding conjugate variables s = dL/dr&re found in terms of the angular velocity of the satellite in the following paragraph. The parametrization of 50(3) by Rodrigues parameters can be written as ( *j) = a / 1 + r\ - r\ - r f ( i 2 +- s ) * U + Mjy 2[r r -r ) h . M\P\2\ x ' 2 r r r 1 3 2 2{ r - r ) l-rl „2+ rl „2 -4Jl ri 2 I 2 ( r r + n) 2 3 3 2(nr + r ) 2(r r -n) 1 - r\ - r\ + r\ 3 2 2 3 where (a y) denotes the orientation matrix of the body frame with respect to the t 139 orbital frame. Using (a y) (d ) r t tJ one finds the time rate of Rodrigues parameters in terms of U (analogous to the method shown in Appendix 4-II): r T2 I = - I n r V h) r * 3 r 2 _ + r f 2 r l r — r r l 3 + 2 r r r r 2 l+ + 2T3 r Next, the conjugate variables s = dL/dr L(r, l 2 3 1 + ri 3 - r i 10. rl are computed, where L is thought of as r ) after the appropriate substitutions for U and (a,y): ^1 + ^31 n + 1 + |r 2 ^3 + 33 a Therefore, n = 6.3 1 + r ?|2 (1) ri/ 2 T h e Hamiltonian up to the Fourth Degree Equation (1) implies that the equilibrium ((a y),n) = (E,0) G SO(3) x so(3) t corresponds to (0,0,0,0,0,2/ ) in (r, s ) coordinates. To bring the equilibrium to 3 the origin the following canonical transformation is defined: The substitution of 0 and (a,-y) in terms of r and 5 followed by a further substitution in terms of f and s in the expression of the Hamiltonian (Chapter 3, Section 3.7) 140 yields the expression H = H + H + H + H + H + ---, 0 1 2 3 4 where #1 = 0, Here, I 2 0 h 0 _1+ 1 V ^ 0 l -12(/x 0 0 0 *3 1 i2_ 2J 0 0 o z h 0 1 /) 0 2 0 1 0 0 0 0 0 1 4/ 0 0 2 0 0 0 0 \ 2 0 0 12(-h + 0 2/ + 2 0 ^3 J The homogeneous polynomials of degree three H and four H4 are given in Appendix 3 6-1. The constant term Ho plays no role in what follows and is omitted. The absence of the first order terms (Hi = 0) confirms that (f, 5) = (0,0) is an equilibrium. The new variables (f, I) are renamed (f, s) in the following sections as the original variables are no longer needed. Unlike H, which is defined on the phase space, the expressions H , H ,... 2 3 are not invariantly defined. They are functions on the coordinate space Rf?g) whose points are denoted by (r, s); in particular, the expressions for H 2 and H 3 change depending on the chosen coordinate system. That is why it is important to recognize the variable with respect to which H is expanded. 141 6.4 Normalization Procedure In this section, a sequence of canonical transformations bring the successive homogeneous polynomials of a given degree of the Hamiltonian to a standard form called the normal form. 6.4.1 Normalization up to the second degree Second degree terms H W by l 2 can be brought to the form ( _ L _ ) + W 2 ( _ a linear canonical change of variables (r, s) ) + W g( _ ) (x, y). The procedure outlined in this chapter shows how a symplectic basis leads to this linear change of variables. The dynamical system associated with H is denned by (H-° o)(l> E Hence, the linear part of this system comes from H . As a result, 2 is the linearization of the equation (2). If the eigenvalues of L are distinct and pure imaginary, then by a linear canonical transformation /o, equation (3) is transformed to 142 where / x 0 L / o is of the form ( 0 0 0 0 0 0 0 0 0 0 0 0 0 W Vo 2 0 0 - W l 0 0 0 0 0 U>3 —(JJ 2 0 0 0 0 0 \ 0 -w 0 0 0 J 3 with w,-'s real. The sign of w,-'s affect the stability of the equilibrium for the complete system as will be shown in Section 6.5. The conjugate variable r 3 and s 3 are decoupled from the rest. The equation (2) may be written in the following form: _1 0 •1 + 0 0 3 V 0 0 0 0 0 0 0 0 0 0 0 12(A - / ) 1 0 4/ _ J3_ A (ri\ 2 0 h r 4 0 0 7*7 0 12(Ji - J ) 3 0 0 0 0 13_ 2/ J3_ •1 + 2/i 0 0 2 1 a oJ fri\ T2 \s J 2 =L T2 S\ \s J 2 It may be shown that the eigenvalues of a linearized Hamiltonian always appear in groups of A, A, — A, —A [25] (Figure 6-2); hence, in depicting the spectrum (location of eigenvalues in the complex plane) of a linear Hamiltonian system it suffices to focus on the first quadrant of the complex plane. This program is carried out for the system defined by equation (2): the first quadrant of the spectrum of equation (2) is plotted corresponding to each pair of parameters in the Ti,T2-space, at the chosen intervals, on the right upper corner of the point (Ti,T ) 2 S2 r \s J 3 3 Consider a subset of the above equations, •Si S\ (Figure 6-3). To avoid the confusion of adjacent spectra they are marked in alternating colors. The eigenvalues 143 of L are shown as crosses. The scale of the complex plane for the spectra is chosen so that the farthest eigenvalue from the origin, for all selections of parameters, is contained entirely within each subplot. The equation Det(L—A) has four pure imaginary roots —1X2,-1X1,1X1,1X2 -T1T2 = C = -AT1T2 2 - A C = + B X 4 2 + C = 0 provided 2 < 0, r T - 6Tir + 2 B X 3T + l > 0, + B = 2 2 2 9T + 6T + 1 > 0. 2 1 4 T i T 2 2 + 2 2 The region in the parameter space defined by Ti + T 2 > 0 and the above inequalities is hatched in Figure 6-4. This region characterizes all those Hamiltonian differential equations (2) whose linearization L at the considered equilibrium has only pure imaginary eigenvalues. From now on, the attention is confined to this region. Let Then 0 < A i < A . Let vi be the eigenvector corresponding to the eigenvalue t'Ai. 2 Define e i N ? = I m v i • I E 0 ) R e v i , where kA = 1 and iVi > 0. Define bi = — R e v i , I\i I m v \ . e i N i Then Lbi = - ^ - L R e v i = - ^ - I m v N i L64 — — — l i l r n v i e i N i = -Axei^"^ x N i = — 1 £ i N i — R e i N i e v 1 = Aiei&i. = - X i € i b 4 , 144 (a) x—*• x—x-x- -x-x -x—x X (b) F i g u r e 6-2 Some possible spectra for linear Hamiltonian systems of (a) two and (b) three degrees of freedom. F i g u r e 6-3 The spectra of the linearized equations of motion on a lattice of points in the parameter space. The spectra in the first quadrant of the complex plane appear to the upper right of each point of the lattice. F i g u r e 6-4 Parameters for which the equilibrium is stable (hatched region) for the linearized equations of attitude motion. 147 Similarly, define 6 = -rz-Revi, iv 2 2 1 Then, i n the basis -Imv-i ^ reads {61,62,64,65}, 0 0 eiAi 0 0 0 0 c A -«iAi 0 0 0 0 — 62^2 0 / V \ 2 2 J 0 Define 0 0 0 0 1 0 2V/I3A3" 0 0 0 0 where A3 = y/STs. In the basis / V ^ * * s {61,62,63,64,65,66}, 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 U>3 0 0 0 0 0 0 0 0 0 —0J 0 0 2 -w 3 0 n c a n o n i c a l form \ 0 J where u>\ = e i A i , u — e A , and W3 = A . W i t h a canonical change of variables 2 2 2 3 (61 the expression H 2 6 63 2 64 65 6 ) 6 ( .j), becomes A + yf i = i in x, y variables. The characteristic frequencies w 's are functions of T and T , t x 2 defined on the hatched region of Ti,T -parameter space. Variations of u>\, u , w 2 2 3 with respect to the parameters T\ and T are presented i n Figures 6-5 to 6-7. 2 Contour plots cover those regions of T\, T -space for which the corresponding w- is 2 t denned. Note the change in the sign of oj\ during transition from T > 0 to T < 0 . 2 2 148 F i g u r e 6-5 Variation of u>i as a function of T\ and T . 2 Ti -1.0 • f,, as a function of T and T . F i g u r e 6-6 Variation of « as a t V a 2 F i g u r e 6-7 variation of W3 as a function of T\ and T. 2 151 6.4.2 In cal Elimination of the third degree terms this subsection, a generating function K is found which yields a canoni- transformation leading to the elimination of terms of degree three in the new Hamiltonian. D e f i n i t i o n . The characteristic frequencies wi, u , u> 2 satisfy a resonance relation of 3 order M if there exist integers m,- not all equal to zero such that m i w j + m oj + m w = 0, 2 2 3 | m i | + | m | + | m | = M. 3 2 3 D e f i n i t i o n . A Birkhoff normal form of degree s for a Hamiltonian is a polynomial of degree 5 in the canonical coordinates (u -,v ) which is actually a polynomial (of t t degree [s/2]) in the variables r» = ( u + v ) / 2 . 2 2 It will be shown that if the characteristic frequencies w,- do not satisfy any resonance relation of order three or smaller, then there is a canonical coordinate system in a neighbourhood of the equilibrium position in which the expression of the Hamiltonian in the new variables is reduced to a Birkhoff normal form of degree three. First, make a complex linear canonical transformation: X z _ = —h tVJ, 2 y - i- + w. it W i t h this H 2 becomes ihJ\Z\VJ\ Note that ZJVJJ = + iu) Z VJ 2 2 + 2 ibJ ZzWz. 3 —i(x - + J/y)/2. B y Jacobi's theorem [26], the implicit relations 2 _ w dK = az define a canonical transformation (z, w) dK w u = ov (u,v) about (ZQ,WO), provided the gen- 152 erating function K(v, z) satisfies the condition dK det—^—^y ,z ) ^ 0, ovoz 2 0 with vj = dK/dz(v ,zo). Q 0 0 Let K = v • z -\- K {y,z), 3 where K 3 is a homogeneous polynomial of degree three. Then dK det——(0,0) ovoz 2 W = = det~E = 1, V + u = z+ --—. ov Since K is analytic, equation (4) determines z as a unique analytic function of u and v which maps (0,0) 0 [27] . L e m m a . Expansion of w and z in variables u and v up to terms of degree two is found by substituting u whenever z appears on the right hand side of the following expressions: z = u- dK 3 , Proof: a,Q = 0 in the expansion of 2 = a*o + A i v + A 2 w H about (u,v) = (0,0), as z(0,0) = 0. To obtain the first degree terms in the expansion of z it suffices to ignore dKs/dv i n the expression z — u — dK /dv(v, 3 z) as dK /dv 3 after substitution for z starts with the second degree terms. Hence z = u + 0(2), where 0(2) represents a polynomial of degree two or larger. To obtain the expansion of z to terms of degree two it suffices to substitute terms of degree one or less for z in dK /dv(z,v) 3 higher degree terms will lead to third or higher degree expansions of z. Hence, z=u-^-(u,v) ov + 0(3). as 153 Similarly, rf=v+^{u,v) + ou 0{3). Consider the expansion of H in terms of the new conjugate variables (u, v) up to the terms of degree three: H(z, w) = = H(u,v) X > y ( u y - ^ 7 ) ( » y + ^ 3= 1 v — dK j=l Here, u m v SfhH n m 3 * ( * > * ) + O f t dK 1 2 is the coefficient of the z v m n 3 n 3 v , with |m | + \n | = 3 which is an SmS^r ^ ^ ^ "™ ^ ) 2 H 3 Any term in K may be written as s^z 1 + 3 ^ 3 3= 1 abbreviation for ) 3 with mi+m + rn3 + ni+n2+n3 2 — 3. monomial. Hence, the coefficients of in the third degree terms of H reads as Smn[iui{ni where hfha l s ~ mi) + iu (n 2 2 - m) 2 the coefficient of the monomial »'w (ra - m )] + + 3 u v m n 3 in fc^g, 3 H (u,v). 3 (5) If there are no resonances of order smaller than or equal to three among w,'s, then may be chosen so that expresssion (5) becomes zero. Therefore, an appropriate choice of K eliminates the third degree terms of H. 6.4.3 N o r m a l form through the fourth degree terms In this subsection, a generating function is found which normalizes H through terms of degree four. previous subsection. The procedure is analogous to the one presented in the 154 Expand H through terms of degree four in variables (u, v) as follows: 3 H(u, v) = ibjjUjVj + H^u, v) + 0(5). ;'=i Let K (u, f ) be a homogeneous polynomial of degree four. The implicit relations 4 _ 8K ou ? 4 ^ dK 0$ define an analytic canonical transformation (tT, v) (£*, $) close to identity. As before u=i-^-(u,v) + 0(4), v = $+—^{u,v) + 0{4). ou Then 3 = Y iu} jtj$j E + { ms[»'wi(ni - mi) + iu (n s 2 2 - m) 2 |m | + |n 1=4 J= l + » w ( n - m )] + 3 3 + 0(5). 3 Unlike the third degree case, even in the absence of resonances of order four or less, it is not possible to eliminate all the fourth degree terms. For m,- = n,-, the monomials /imnf C* m n remain. Therefore K4 i = 1,2,3 may be chosen ( the choice is not unique) such that 3 B(t,$) = Y^iujtjtj + J2 j=l 6.4.4 hrnnrr +0(5). 1 m=n N o r m a l i z a t i o n b e y o n d t h e second degree t e r m s — L i e series In this section, a procedure is outlined which yields the normalized form without recourse to a generating function. First the general procedure is described. 155 Let P denote the set of homogeneous polynomials of degree r. The direct r sum of vector spaces P , P = ©£2_ P , is a positively graded real (or complex) Lie r 2 algebra [28] with r Y— — -—— \F,G] = d X J d y l dy 3 dX 3 W i t h Zj — Xj + iyj and Zj = Xj — iyj, Subscript elements H E P - Let K E P . Define r r a ad :P -> P, K [H, K] as a graded homomorphism of degree s — 2. Denote ad {H ) = j K r ad {ad -\(H )) j K 4 For 5 > 3 and fixed n, the equality r + j(s — 2) = n has only a finite number of solutions in the positive integers r and j. A s a consequence, we can define a linear operator exp(adx) • P —* P by ex (ad )(H) P = K f : ^ l j=o J D e f i n i t i o n . A n element H = ©£2_ .fir € P is i n the normal form through terms of 2 r order m > 2 if adn [H ) — 0 for 2 < r < m . Here, it is assumed that 2 r This definition agrees with the one used before. Proposition. Let H = © £ L i I order (m - 1) > 2. Let H 2 m r = F m E P be in the normal form through terms of + G, m where F m E Ker(ad \pm) H2 a n d m G G 156 Then exp(adK ){H) ran(adH \Pm)- E P is i n the normal form through terms of m 2 order m , agrees with H through terms of order m — 1, and has F as m m th term. Proof: exp(adK )(H) - H + H -\ m 2 hH 3 + adK m [H ) + (terms in Py with 2 m j > m + 1}. Therefore, H + ad (H ) m Km + G - ad {Km) = F 2 m m = H2 F. m A routine computation shows that ad {z*Z*) = -i(m H2 - ti).uzTZ . n As a result, if G =E - " ^ ^ 6 m e r a n i a d H \ p 2 ) , m then Km = iJ2K m - n)^\- b z' Z . 1 h ii mH The above procedure is applied to the Hamiltonian H = H + H3 2 Let H ^Y.^fhnZ Z . fh Then, K K z exp{K ){H) 3 + if + • • • • 4 = » E [ ( n » - n j . w ] - ^ . Therefore 1 3 = H + H + H + ad {H ) 2 3 On the other hand adK (H ) 3 terms i n exp(adK ){H) 3 2 4 Ks + ad {H ) 2 Ks = — adff (K ) 2 + \ad\ {H ) 3 z 3 Jf4 = H + ad (H ) = H + 4 — H< + Ks ^ad K + 0(5). — —H . Therefore, the fourth degree 3 become A 2 - 3 3 H 3 ^\H ,K \. 3 3 ^ad H Ks 3 157 Once H4 is computed and expressed in the complex form c %az Z m n r the fourth degree terms in the normal form expressed in complex variables are of the form Y^rn=n ™nZ Z . c rn n In summary, to compute the fourth degree terms in the normal form the following three steps are taken. i . The expression ii. H 4 [#3,^3] = H4+ ^[H ,K ] iii. The coefficients 3 3 Crnn is computed in complex variables. is found in complex variables Yl rhH^ Z . e m n with rh = n are selected. These yield w,y's as follows "11 = 4c 0020Cb 2 = W12 = w 1 3 = 2cuono, = 2c iioij 10 W 2 = 4c 200205 2 W32 - 0 W 3 2 ^33 = 2conoii, = 4C002002- There is a third method for deriving the normal form of Hamiltonians, which is described by Gustavson [20, pp. 674-75]; this method is a variant of Birkhoff's approach described in subsections 6.4.2 and 6.4.3 with the difference that the generating function is real and leads to a real canonical transformation . In Birkhoff's method the canconical transformation bringing a Hamiltonian into the normal form need not be real. However, all these methods yield the same normal form as the normal form is unique in the nonresonance case [20, p. 675]. 158 6.5 Stability of an E q u i l i b r i u m — t h e Delp Case In this section, a sufficient stability condition is used to identify a subset of the parameter space for which the equilibrium under consideration is stable. The method depends on the normal form of the Hamiltonian up to the terms of degree four. The normalized Hamiltonian of a satellite in a circular orbit in the canonical variables £, r? has the form t'=l i,j where w,-'s and w -y's are determined as functions of T\ and T . t 2 The following is the stability condition found by Thkai [5]: = (0,0) is stable if 3 Vr -> 0 t 3 u; r -= 0 and fir^O t w.yr.ry = 0. t »=i (6) i,j=l To verify this condition the following approach is adopted. Define <f> = r (k;,y)f. The minimum and maximum values of <j> subject to the conditions 3 r >0, t i= 1,2,3; 3 ^w r i=i t t = 0; J^ * = 7 1 t=i are found. It is required that either the minimum value be positive or the maximum value be negative. This is shown to be equivalent to condition (6). To find the extrema of <f> on the intersection P of planes P i : r i + r + r = 1 and P : Yli=i 2 0, form the function (the method of Lagrange multiplier) 3 V» = <f> + A i £ J2 3 n- l) + A 2 3 2 w i t = r 159 Denote the unique solution to /Wll <90 3(f, A i , A ) ->12 W13 1 1 ( T l \ U> W21 W22 <~23 W31 W32 w 33 1 1 1 1 0 0 W2 OJ3 0 0 2 V OJl 0 T2 2 + J T3 Al u / 2 + 0 - 1 v 0 y by ( f b , A i , A ) . 0 2 o The line P intersects the three planes r,- = 0 , 0 "3 w -u> -^2 3 t = 1,2,3 (Figure 6-8) at points and Da = Do = 2 CJ -WI --1 w -W! 0 2 2 w -u;i 3 If ct),-'s change sign, then the following argument shows that two of the points D,- satisfy r,- > 0, t = 1,2,3. Two of the u>i/w > W 1 / W 3 , W 2 / W 3 are negative when 2 w i , W2» and w are not all of the same sign. The reciprocal of the nonzero components 3 of Dt, i = 1,2,3 are respectively 1_ 1 fl 0J3 _ 1 > U>2 fl D3 0J3 1 - - * fl 1 _ U>2 Therefore, two of the D , ' s are in the first octant, as their nonzero components are positive. The above argument holds true in all but exceptional cases (of measure zero in the parameter space-T\, T ) for which w 's of the same sign have equal 2 t magnitudes. Assume, with out loss of generality, that D and D are in the first octant. Let 2 3 4>i — 4>(Di), i — 2,3. In the following chart, if the sequence of verifications does T3 D N 2 D1 F i g u r e 6-8 The geometry of Thkai's theorem. 161 not lead to "otherwise" the relation (6) is satisfied: (D 2 T) 0 •(£3 - ( \ ( I T) 0 mm(c/> ,02,^3) > 0 or max(<f> , <£ , <j>z) < 0, otherwise, min(<f)2,<f>3) > 0 or max(<l> ,4>z) < 0» otherwise. o 0 2 2 The chart is based on the following argument: for any f ^ 0 in the first octant satisfying <f>(f) = 0, there exists a A G R such that Ari + A r + Ar3 = 1 for some A G R. 2 of <f> on the segment D\D To ^ D1D2 and fo G P\DiD 2 2 <£(Ar) = \ <J>(T) 2 = 0 and The maximum M and minimum m value is detemined in each of the two cases corresponding to The condition (6) is satisfied when 0 G" [ m , M ] . The curves in the Delp region corresponding to which the characteristic frequencies satisfy a resonance relation are shown in Figure 6-9. Outside these resonance curves the normalization procedure up to the fourth degree is possible and the relation (6) demarcates a subset of the Delp region for which the considered equilibrium is stable (Figure 6-10). F i g u r e 6-9 Resonance regions shown as a collection of curves in the parameter space. F i g u r e 6-10 Status of the equilibrium as a function of T\ and T . 2 164 References [1] Newton, I., Principia, Vol.2: The system of the world, Trans. Motte, A . , Revised by Cajori, F . , University of California Press, Berkeley, 1962, pp. 48485. [2] Lagrange, J . L . , "Recharches sur la libration de la lune," Oeuvres VI, Editor: M . Serret, Gauthier-Villars, Paris, 1870, pp. 5-61. [3] Lagrange, J . L . , "Theorie de la libration de la lune," Oeuvres V, Editor: M . Serret, Gauthier-Villars, Paris, 1870, pp. 5-122. [4] Poincare, H . , "Sur l'equilibre d'une mass fluide," Acta Math. VII, 1885, pp. 293-99. [5] Laplace, P.S., Celestial Mechanics, Vol. 2, trans. Bowditch, N . , Chelsa Publishing Company, New York, 1966, pp. 933-67. [6] Tisserand, F . , Traite de Mechanique Celeste, tome 2, Gauthier-Villars et fils, Paris, 1891, pp. 444-75. [7] Routh, E . J . , Advanced Dynamics of a System of Rigid Bodies, Dover Publications, New York, 1955, pp. 369-84. [8] Beletiskii, V . V . , "The librations of satellite," Artificial Earth Satellites, Vol. 3, Editor: Kurnosova, L . V . , plenum press, New York, pp. 18-45. [9] Delp, R . H . , and DeBra, D . B . , "Rigid body attitude stability and natural frequencies in a circular orbit," The Journal of the Astronautical Sciences, V o l . 8, 1961, pp. 14-17. [10] Likins, P . W . , "The influence on dynamics and control theory of the spacecraft attitude control problem," Proceedings of the Sixth Canadian Congress of Applied Mechanics, Vancouver, B . C . , Canada, Editor: V . J . Modi, 1977, pp. 165 321-35. [11] Kane, T . R . , Likins, P . W . , and Levinson, D . A . , Spacecraft Dynamics, McGrawH i l l , New York, 1983, pp. 199-210. [12] Poincare, H . , Les Methodes Nouvelles de la Mecanique Celeste, V o l . 1, 2, 3, Gauthier-Villars, Paris, 1892, 1893, 1899. [13] Birkhoff, G . D . , Dynamical Systems, American Mathematical society, Providence, Rhode Island, Revised Edition 1966, pp. 59-85. [14] Arnold, V . I . , "Proof of A . N . Kolmogorov's theorem on the conservation of conditionally periodic motions with a small variation in the Hamiltonians," Russian Math, survey 18, 5, 1963, pp. 9-36. [15] Siegel, C . L . and Moser, J . K . , Lectures on Celestial Mechanics, Springer-Verlag, New York, 1971. [16] Arnold, V . I . , "Algebriac unsolvability of the problem of Lyapunov stability and the problem of topological classification of singular points of an analytic system of differential equations," Functional Analysis and Rs Applications, Vol.4, N o . l , 1970, pp. 173-80. [17] Tkhai, V . N . , "The stability of multidimensional Hamiltonian systems," PMM U.S.S.R., Vol. 49, No. 3, 1985, pp. 273-81. [18] Williamson, J . , "On the algebriac problem concerning the normal forms of linear dynamical systems," Amer. J. Math. 58, 1936, pp. 141-63. [19] Laub, A . J . and Meyer, K . , "Canonical forms for symplectic and Hamiltonian matrices," Celestial Mechanics, 8, 1974, pp. 213-38. [20] Burgoyne, N . and Cushman, R., "Normal forms for real linear Hamiltonian systems with purely imaginary eigenvalues," Celestial Mechanics, 8, 1974, pp. 166 435-43. [21] Gustavson, F . G . , "On constructing formal integrals of a Hamiltonian system near an equilibrium point," The Astronomical Journal, Vol. 71, No. 8, 1966, pp. 670-86. [22] Deprit, A . , Henrard, J . , Price, J . F . , and Rom, A.,"Birkhoff's normalization," Celestial Mechanics, 1, 1969, pp. 225-51. [23] Churchill, R . C . , Kummer, M . , and R o d , D . L . , "On averaging, reduction, and symmetry in Hamiltonian systems," Journal of Differential Equations, 49,1983, pp. 359-414. [24] Markeev, A . P . and Sokolskii, A . G . , "On the stability of relative equilibrium of a satellite in a circular orbit," translated from Kosmicheskie Issledovaniya, Vol. 13, No.2, 1975, pp. 139-46. [25] Abraham, R . and Marsden, J . E . , Foundations of Mechanics, The Benjamin/ Cummings Publishing Company, Reading, Massachusetts, 2 nd Edition, 1980, p. 169. [26] Arnold, V . I . , Mathematical Methods of Celestial Mechanics, Springer-Verlag, New York, 1980, p. 267. [27] Dieudonne, J . , Foundations of Modern Analysis, Academic Press, New York, 1969, p. 272. [28] Bourbaki, N . , Elements of Mathematics, Algebra, Part I, Addison-Wesley, Reading, Massachusettes, 1974, p. 457. 167 A P P E N D I X 6-1: H = 3 HIGHER TERMS O F T H E HAMILTONIAN —r sis 4iii2 A T -h T 3 + 2 + 2 s i s H 3 — ri~ si — 2iii 3 2 2i; r T 3 3 ^ h ( h - h ) H—77-7—r!S s 4i i , h 2 ,^3 + T 4Iil + h 2 A H 3 3 r i 5 3 + — 2 27r — r r s 2 3 2 , - 2 r 2 S 3 + 2 r 3 5 3 -12/x/f + 127x72/3 + / i / — - / / rxr r 2 2 3 2 3 2 3 ixJ2 HA h h , /s(/i-/ ) + h h - h h 7TT7 = 2 r 2~ s 3 + 3 2 rxr|5 3 4/1/2/3 2 l 2 + 2 5 3 2 3 3 T7- 3 3 + r T 7 s r l 2T S r 3 3 Jx-12 J . f or 1 f U , , / o j c r U . 2 1 -24/i/ + (-oh + i J r x + 2(dl2 - 51 )r r H 4i 2 4/x/2 —/l + 7 , - 2 2 , - 2 2 + -27-—-i-sra +^ r + — r , 1 , / ( / l ~ h) 3 2 3 3 3 l + 16/x/ + / | 2 3 — 2 + for 11 7 V + 3 6 / i / - 2 4 / | - 2 0 / / + / + (9ix - H i 3 j r + — 4 2 2 2 3 2 3 2 rxr 2 3 0 2 , 3(6/x - 8 / + / ) r r H r 2 2 3 3 Epilogue A science produced is impossible; nated; if we cleave immediate and truths there result, will with a view single are fruitful to those the only only connecting be no longer a to its if they of which links applications are concate- we expect will be an lacking, chain. Henri Poincare Suggested Further W o r k The following topics of research constitute a natural sequel to the subjects covered in this thesis. • In derivation of the potential energy of a satellite as a function of its orientation with respect to a central field, one ignores terms of order three or higher in d /R m (Section 2.1.2). These are the terms which contribute to subtle differ- ences such as the disparity between the gravitational potential field of a cube and a ball with equal moments of inertia. Moments of inertia is a tensor of rank two representable as a matrix (7,-y) with nine entries. In analogy, it might be possible to define a tensor of rank three representable as a triply indexed family (Iijk) of twenty-seven entries which reflects the finer characteristics of the potential field of an object. To give (Iijk) the property of a tensor—obeying the appropriate transformation laws—it might prove necessary to define Iijk by an expression more elaborate than — / 168 XiXjXkdm for repeated indices. Recall 169 that for the inertia tensor hi = / Xk ~ i j) - x x dm k=l J Here, the diagonal terms are defined by a different formula than the off-diagonal terms. A n adequate definition of (Iijk) a n d a study of its invariants and prin- cipal directions, parallel to the study of inertia tensor, could yield some new results relevant to the motion of rigid objects under potential field. • The method developed in Chapter 3 assumes a time independent Lagrangian. In practice, however, there are instances for which the dynamics of the system is best modelled by a time dependent Lagrangian: the motion of a satellite during deployment and retrieval of appendages or when the environmental effects are significant is often expressed in terms of a time dependent Lagrangian. To encompass these cases, it is desirable to extend the presented method. Some of the relevant material on integral curves of a time dependent Hamiltonian is developed in [1, pp. 233-45] and [2, pp. 265-67]. The material might give leads as to the method of deriving the equations of motion from a Lagrangian in the form L((a.ij) t U, A , A , t) without use of generalized coordinates for the space of rotations. • In Chapter 4 a set of rational generalized coordinates was introduced corresponding to which the evaluation of the right hand side of differential equations required a shorter C P U time than those of transcendental rotational coordinates such as the Euler angles. It was further verified that, in some instances, this translated to a significantly lower C P U time for a given integration. The required C P U time for numerical integration of x = / ( ^ in a time interval between ti and t , however, is not only a function of f(x) but also dependent 2 on the time step of the integration. The time step in a numerical integration 170 affects the local error and consequently determines the degree of accuracy of the final result. A study which estimates the time step of integration for a fixed level of error tolerance will shed further light on the effective utility of a given coordinate system. Such a study will further indicate how the speed of integration varies with the initial condition. • In connection with the material in Chapter 5, it is of interest to know if there are any other integrals of motion save the energy function for the attitude motion of a rigid satellite. It is likely that there are no further integrals of motion. Proof of the nonexistance of further integrals of motion or the development of other integrals, each has its own distinct technique described i n literature [3, p. 178] and [4-6]. The techniques for the former alternative are non algorithmic and complex in nature. • To carry out an investigation parallel to that of Chapter 5 for a satellite with two symmetrically deployed flexible panels one determines the equilibria by setting dH = 0, where dH is given in equation (2), Chapter 3. After a classification of equilibria (analogous to Table 5-1) a local analysis of the Hamiltonian at the equilibria identifies those regions in the appropriate parameter space for which a given equilibrium is stable (analogous to Figures 5-2a to 5-2f). Then through a comparison of energy levels at equilibria, as in section 5.3, one stratifies the parameter space according td the maximum allowable perturbation energy (analogous to Figure 5-4). • A straight forward extension of the investigation reported in Chapter 6 is to normalize the Hamiltonian through terms of degree 2n, n>3, and apply the generalized Tkhai's theorem [7, p. 277] to establish a more comprehensive set of parameters in the Delp region at which the equilibrium orientation is stable. 171 One may approach the problem from a different direction. A n adaptation of Arnold's proof of nonexistance of a finite step algorithm for establishment of stability of the equilibrium for his differential equation [8] to the stability of equilibrium orientation in the Delp case might show the impossibility of a complete answer to the stability of the equilibrium orientation in rigid satellites. This may render any further attempt for a complete answer futile. Concluding Remarks This paragraph outlines the original contributions of the thesis. • In the area of formulation, an alternative transition from a Lagrangian L((aij),U,A,A) to the corresponding equations of motion is devised based on geometrical properties of T*SO{3). • In the area of simulation of attitude motion, a nonsingular global coordinate system on the space of rigid motions with a fixed point (50(3)) is found which yields a rational expression for the entries of the rotation matrix. The parametrization is obtained as the composition of the inverse of stereograhic projection of 5 onto SO(3). 3 and a homomorphism mapping the universal cover of 50(3) For visualization of attitude motion, a homeomorphism between 50(3) and RP(3) is utilized to represent every element of 50(3) as a vector in the unit closed closed ball; the axis and the amount of rotation (in fractions of 7r) are determined by the direction and the magnitude of the vector, respectively. • In the area of analysis of attitude motion of rigid satellites in a circular orbit the following findings are notable. (a) There are 24 equilibria (eiy)f =1 f =x which may be classified into six group 172 of four based on the dynamical behaviour. (b) The existence of a stable equilibrium orientation is a result of the compactness of 5 0 ( 3 ) . (c) The notion of a state remaining untumbled about a stable equilibrium is introduced. The energy levels of equilibria are compared to find the maximum allowable perturbation energy before the satellite tumbles away from a given stable equilibrium orientation. This method is applicable to any Hamiltonian vector field. (c) The bounds of attitude motion about an equilibrium is found by considering equipotential surfaces in the configuration space. (d) The long outstanding problem of stability of the equilibrium e n is resolved by normalization of the Hamiltonian about e n through terms of degree four and application of Tkhai's theorem. 173 References [l] Arnold, V . I . , Mathematical Methods of Classical Mechanics, Springer-Verlag, New York, 1980. [2] Choquet-Bruhat, Y . and Dewitt-Morette, C , Analysis, Manifolds, and Physics, North-Holland, Amsterdam, 1982. [3] Arnold, V.I., "Small denominators and problems of stability of motion in classical and celestial mechanics," Russian Math. Surveys, 18, 1963, pp. 85-192. [4] Sumbatov, A . S . , "On the problem of determining the existance of ignorable coordinates in conservative dynamic systems, " PMM U.S.S.R., Vol. 42, No. 1, 1978, pp. 43-51. [5] Sumbatov, A . S . , "On ignorable coordinates of conservative and natural systems with three degrees of freedom," PMM U.S.S.R., Vol. 45, No. 5, 1982, pp. 58897. [6] Goldman, L . , "Integrals of multinomial systems of ordinary differential equations," Journal of Pure and Applied Algebra, 45, 1987, pp. 225-240. [7] Tkhai, V . N . , "The stability of multidimensional Hamiltonian systems," PMM U.S.S.R., Vol. 49, No. 3, 1985, pp. 273-281. [8] Arnold, V.I., "Algebraic unsolvability of the problem of Lyapunov stability and the problem of topological classification of singular points of an analytic system of differential equations," Functional Analysis and Rs Applications, Vol. 4, No. 1, 1970, pp. 173-80.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- On some aspects of dynamics, modelling, and attitude...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
On some aspects of dynamics, modelling, and attitude analysis of satellites Marandi, Said Rashed 1988
pdf
Page Metadata
Item Metadata
Title | On some aspects of dynamics, modelling, and attitude analysis of satellites |
Creator |
Marandi, Said Rashed |
Publisher | University of British Columbia |
Date Issued | 1988 |
Description | The thesis identifies several limitations in the modelling and attitude stability analysis of two classes of spacecraft: rigid and flexible satellites. Attractive methods are proposed which promise to have far reaching consequences in spacecraft dynamics. These alternatives, developed based on techniques of differential equations, classical mechanics, and differential topology, are indicated below. (a) An Alternate Transition from the Lagrangian of a Satellite to Equations of Motion The classical procedure requires the Lagrangian to be expressed in terms of the corresponding generalized coordinates of the problem. This requirement significantly complicates the derivation of the equations of motion through an introduction of a set of librational generalized coordinates, which is strictly not a part of the dynamical system. Using the Lagrangian in the natural variables (angular velocity, direction cosines, and vibrational coordinates), one develops a procedure for derivation of equations of motion without an a priori choice of rotational generalized coordinates. For the case of a satellite with two flexible plate-type appendages, for example, the approach reduced the formulation time to one-third. (b) Synthesis and Depiction of Rotational Motion of Satellites and Robots The rotational coordinates in use for numerical prediction of orientation of a satellite are either singular or redundant. Furthermore, they lack a convenient visual interpretation. A new set of coordinates is proposed and an associated representation is developed which avoids these limitations. The procedure is applied to represent and integrate numerically the librational response of the flexible satellite mentioned in (a). (c) Resolution of Attitude Stability of Delp Satellites The development here tackles a long outstanding problem in the area of attitude stability of satellites. The resolution of this problem through normalization of the Hamiltonian leads to a better appreciation of stability associated with the class of gravity gradient structures such as the proposed Space Station. |
Subject |
Satellites Artificial satellites -- Attitude control systems |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-10-07 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0098148 |
URI | http://hdl.handle.net/2429/29019 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1988_A1 M36_5.pdf [ 7.81MB ]
- Metadata
- JSON: 831-1.0098148.json
- JSON-LD: 831-1.0098148-ld.json
- RDF/XML (Pretty): 831-1.0098148-rdf.xml
- RDF/JSON: 831-1.0098148-rdf.json
- Turtle: 831-1.0098148-turtle.txt
- N-Triples: 831-1.0098148-rdf-ntriples.txt
- Original Record: 831-1.0098148-source.json
- Full Text
- 831-1.0098148-fulltext.txt
- Citation
- 831-1.0098148.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0098148/manifest