O N S O M E A S P E C T S O F D Y N A M I C S , M O D E L L I N G , A N D A T T I T U D E A N A L Y S I S O F S A T E L L I T E S S a i d R a s h e d M a r a n d i .S. (Mechanical Engineering and Mathematics), California Institute of Technology, 1981 M.A. (Mathematics), University of California, Berkeley, 198S A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E DEGREE OF D O C T O R OF P H I L O S O P H Y in The Faculty of Graduate Studies Department of Mechanical Engineering We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y OF BRITISH C O L U M B I A July 1988 © Said Rashed Marandi, 1988 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of The University of British Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date JJ\pft<K DE-6(3/81) A B S T R A C T The thesis identifies several limitations in the modelling and attitude stability analysis of two classes of spacecraft: rigid and flexible satellites. Attractive meth-ods are proposed which promise to have far reaching consequences in spacecraft dynamics. These alternatives, developed based on techniques of differential equa-tions, 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 signif-icantly complicates the derivation of the equations of motion through an intro-duction 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 C O N T E N T S A B S T R A C T ii C O N T E N T S ii i LIST OF T A B L E S vi LIST OF F I G U R E S vii N O M E N C L A T U R E x A C K N O W L E D G E M E N T xiv 1. I N T R O D U C T I O N 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 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 2.3 References 23 A P P E N D I X 24 2-1 E X P A N S I O N OF P O T E N T I A L E N E R G Y 24 • • • i l l 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 T O T H E E Q U A T I O N S OF M O T I O N F O R A L I B R A T I N G S Y S T E M 25 3.1 Introduction 25 3.2 Summary of the Mathematical Tools 28 3.3 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 . 59 3.10 References 61 4. A N E W L I B R A T I O N A L G E N E R A L I Z E D C O O R D I N A T E ( a - P A R A M E T E R S ) 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 A P P E N D I C E S 91 4-1 D I R E C T I O N COSINES A N D Q U A T E R N I O N S 91 4-II A N G U L A R V E L O C I T Y A N D a - P A R A M E T E R S 93 4-III C O N S E R V A T I O N 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 " 1 95 5. U S E OF 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 OF 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 A P P E N D I C E S 129 5-1 T H E C R I T I C A L POINTS OF Ep 129 5- II T H E L O C A L E X P R E S S I O N S OF Ep A T E Q U I L I B R I A . . . . 130 6. 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 . . 138 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) . . 151 6.4.3 Normal form through the fourth degree terms 153 6.4.4 Normalization beyond the second degree terms—the Lie series 154 6.5 Stability of an Equilibrium—the Delp Case 158 6.6 References 164 A P P E N D I X 167 6- 1 H I G H E R T E R M S OF T H E H A M I L T O N I A N 167 E P I L O G U E 168 V LIST OF TABLES 4- 1 Parameters defining the configuration and the initial condition of the hybrid satellite 84 5- 1 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. . 10 2-2 Positions of various reference frames on the hybrid satellite: (a) body 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 R2; (d) (0,1) x [0,1] is a neighbourhood of (1/2,1/2) € R2 but not (1/2,0) 6 R2. . . . 30 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 33 3-3 (a) Diffeomorphic one and two dimensional spaces constrasted with non-diffeomorphic spaces; (b) a trajectory on a torus 39 3-4 (a) Formation of manifolds by pasting; (b) a covering of a sphere by eight local coordinates 41 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 to a manifold at a given point 42 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 45 4- 1 A diagrammatic transition from <r-parameters to unit quaternions, in R4, to the space of rigid motions with a fixed point, in R9. . . . 65 4-2 A diagrammatic correspondance between SO(3) and RP(3) defined as a 3-ball with antipodal identification along its boundary 73 4-3 Time evolution of a rigid object seen form the isometric vantage point. V l l Tick marks on the orientation trajectory indicate intervals of one-tenth of an orbit 74 4-4 Time evolution of the orientation of a rigid object seen from the top vantage point 76 4-5 Time evolution of the orientation of a rigid object seen from the side vantage point 77 4-6 Time evolution of the orientation and angular velocity of a rigid object in a circular orbit 78 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 03 axis (opposite orbit normal). 79 4-8 Time evolution of the orientation of a rigid object in terms of successive rotations about body axes bi,b2,b3 81 4-9 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) 86 5- 1 Graph of the nondimensionalized inertia parameter T 3 as a function of Ti and T2 102 5-2 Index of the equilibrium en, equal to the muliplicity of hatching, as a function of T x and T2: (a) i = l ; (b) i=2; (c) i=3; (d) i=4; (e) i=5; (f) i=6 109 5-3 Maximum allowable disturbance energy (nondimensional) before tum-bling at equilibrium en as a function of T\ and T2 116 5-4 Optimal configuration contours obtained by equialtitude contours of the surface in Figure 5-3 117 5-5 (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 118 5-6 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 oi > 0,o2 > 0, o 3 > 0 in RP(3). . . . 125 6- 1 Two classes of satellites for which the equilibrium orientation is stable in the linearized analysis 137 6-2 Some possible spectra for linear Hamiltonian systems of (a) two and (b) three degrees of freedom 144 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. . . . 145 6-4 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 T2 148 6-6 Variation of io2 a s a function of T\ and T2 149 6-7 Variation of <JJ% as a function of T\ and T2 150 6-8 The geometry of Thkai's theorem 160 6-9 Resonance regions shown as a collection of curves in the parameter space 162 6-10 Status of the equilibrium as a function of T\ and T2 163 ix N O M E N C L A T U R E O p e r a t i o n s a n d f u n c t i o n s A fa B A diffeomorphic to B fog composition of functions f and g dfa derivative of function f at 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 A b b r e v i a t i o n s V for all 3 there exists; for some 3! there exists a unique e element of; in => implies; if then A 2 T i T 2 2^ T3 u (1 + m 2 ) 2 CLij —* H i c m . the homogeneous polynomial of degree i in the formal power series expansion of H center of mass S e t s Q s3 c TC T*C SO (3) so(3) RP(3) S c a l a r s E EP G me 1 Ek {h,h,h} { r i , r 2 , r 3 } < R ,+,*,• > algebra of quaternions unit quaternions configuration space; C — 5 0 ( 3 ) for a rigid satellite and C = 5 0 ( 3 ) x JR 8 for a rigid core with two flexible panels—four modes at each panel tangent bundle of C (the phase space of the dynamical system) cotangent bundle of C simple orthogonal group of dimension three Lie algebra of 5 0 ( 3 ) or the tangent space to 5 0 ( 3 ) at identity real projective space of dimension three energy potential energy gravitational constant mass of earth angular speed of the c m . of a satellite in orbit around Earth kinetic energy centroidal principal inertias of the rigid portion of a satellite Rodrigues parameters xi { 5 1 , 5 2 , 5 3 } coordinates conjugate to {r\,r2,rz} H Hamiltonian {ctxa.2,a.z} sequence of rotations about 0 1 , 0 2 , and 03 L Lagrangian {u>i,u>2, W3} characteristic frequencies of H (first order invariants) { w , y } second order invariants of H mr mass of the rigid portion of a satellite mp 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 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 Es strain energy Aijk 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 dm max\d \ X l l V e c t o r s —* R c m . of the satellite expressed in the body frame, Figure 2 - 2 a d position vector to an arbitrary point of the satellite from its c m . expressed in the body frame, Figure 2 - 2 a { o i , o 2 , 03} unit vectors along local normal, local horizontal, and orbit nor-mal, respectively (orbital frame), Figure 2 - 1 —* —* —* {61 ,62^3} the principal body frame of the rigid section at the c m . of the satellite, Figure 2 - 2 a fi* inertial angular velocity of the body expressed in the body frame h,b2,b3 fi angular velocity of the body frame w.r.t the orbital frame ex-pressed 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 transfor-mation A ( - ^ l l l j - ^ l l , Au2, A 2 1 2 , -4l21)-^221i Ai22, A222)T B conjugate pair to A f (rl5r2) M a t r i c e s ( b o l d c a p i t a l l e t t e r s ) E n identity matrix of order n, abbreviated as E when the order is obvious from the context 111 Il2 Il3 Inertia operator expressed in body frame, | 7 2 i / 2 2 I23 ^31 -^ 32 -^ 33 X l l l ACKNOWLEGEMENT I am grateful to Professor V . J . Modi for his unwavering attention, sup-port, 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 Ki l lam Foundation for their gen-erous support through the pre-doctoral Ki l lam award during the past four years. xiv To my parents, with love X V 1 . INTRODUCTION The transfer of knowlege from one field of specializa-tion to another has often proved inadequate, and this has been particularly troublesome in the flow between 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 indispen-sible tool in exploration of earth resources, weather forcasting, and communication. Satellites will be increasingly relied on as the stepping stone for man to wider op-portunities. Today, design and development of artificial satellites draws on every aspect of human scientific prowess: from classical mechanics to metallurgy, elec-tronics, 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 orbit-ing 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 vibra-tion 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 deriva-tion 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 per-spicuous 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. Sur-prisingly 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 math-ematics 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 La-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 coordi-nates 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 devia-tions 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 orienta-tion 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 satel-lites. 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 O F T H E L A G R A N G I A N 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 R i g i d S a t e l l i t e i n a C i r c u l a r 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 I n t r o d u c t i o n Any 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 & i , 6 2 , & 3 at c with respect to the right handed orbital frame 0 1 , 0 2 , 0 3 (Figure 2 - 1 ) . 6*1 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 G r a v i t a t i o n a l p o t e n t i a l e n e r g y Gravitational potential energy of a rigid object in the gravitational field of a central body of mass me is given by —* —* = --^y ( l - r 2 o 1 . - + ( - ) 2 ) ?dm. Using Taylor's formula with remainder, one obtains (Appendix 2-1) Gmemr 7 2 3 T^ , 1 „,fdm^ EP = ^ ~J 21 0 1 ' 1 R ~R >' provided 0 < dm/R < 1. Here, dm is the maximum distance of a point on the body from the c m . , 7 = y/Gme/R3 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 T h e L a g r a n g i a n Dropping the constant terms from the expression of potential energy and pre-serving only the lowest terms in dm/R, one obtains 3 7 2 Ep = - y O i ' Ioi-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 L = Ek-Ep=1-J2linf-3-1^Ilali, t = i t = i where o x = (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 L=1-ytli(nf-3ali). 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 S a t e l l i t e i n a C 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 S a t e l l i t e c o n f i g u r a t i o n 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 mp. The spine (center line) of the plates pass through o. The panel deflections are assumed to be small; this assumption is more fully explained in Section 2.2.4. 10 Figure 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 l 5 o 2 , o 3 is placed at c. The body frame bi,b2,63 located at c is defined parallel —# —* —* to the principal centroidal frame (X,Y,Z) of the rigid section of the satellite. 2 . 2 . 2 C o m p u t a t i o n o f t h e i n e r t i a m a t r i x As in the rigid case, potential energy of the satellite up to the second degree terms in d m is given by ^ P = ^ - O 1 - I O 1 , (1) where I 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 61,62,63 at c is derived. As shown in Figure 2-2a, frame X,Y,Z is placed at o along the principal directions. Two parallel frames are placed at the base of the plates P\ and P2 with their origins at Ai and A2, respectively (Figure 2-2b). Let the deflection of plates Pi and P2 with respect to 77 A if and r)A2c coordinate systems be given by / i (»7.? ) = X ^ i »y ^ ) < M f ) > { J (2) where Auj and A2ij are the amplitudes of vibration of the upper and lower plates corresponding to the characteristic function il>i<f>j. The location of point c with respect to the coordinate system X,Y,Z is given by e = / d dm 2m p — I + mr J (/ d dm + / d dm) + m r JPi&cPo JB 2mp -r mr Jp1iip2 — f 2mp + mr Jp1&lP2 d dm, 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 2 refer to appendages and B refers to the rigid core. The c m . offset is given by e = 2m„ + m -^{Mii + A2ij) f i/>i{ri)dri f ^y(f)df r . . J-b JO 0 Let e be the first coordinate of e (the only nonzero coordinate). Then, the inertia matrix for the rigid core at c expressed in the body coordinate system 61,62,63 is given by Ii 0 0 \ / o o o I f l = I 0 h 0 + e 2m r 0 1 0 0 0 / 3 / l o o i Next, the moments of inertia of the panels at c expressed in the frame 61,62,63 are computed: J 1 P 1 = 2K (A+ i ) 2 + ^ ((26)a + /2)]; I12 = ~P I f (A + f2)ridvds; Jo J-b Iiz = -P /' I" i f i - / 2)(f + ^2P2 = P / ' I* [(h - e)2 + (S + h)2]dr,d< Jo J-b + pf f b l ( f 2 - e ) 2 + (! + h)2}dr,ds; Jo J-b fl fb fl rb I 2 3 = -P / r){h + $)dnd$-p / n { - h - $)drjd$ = 0 ; ^0 J-b Jo J-b Izz = pflfh [(fi - e ) 2 + n3]dndi -pf f [(/a - 6) 2 + r,2]dr,dt. Jo J-b Jo J-b 14 On substitution of mode amplitudes these entries simplify to the following: A P i = ^ 2 + n + 2 ( / i + ^ ) 2 m p ; 7 f 2 = -pJZiAuj + A2ij) / / <t>j{$)d$; 7 i P 3 = -pJZiAuj - A2ij) / tpi{ri)dri[h <j>j{$)dt; + f ^ y (? )< * f ] ; J 2 P 2 = 2[m(A + -) 2 + ^ /2] + pj J [fl + ffidndc - 2(mr + mp)e2; / 2 P 3 = 0 P _ 2 m p r h b2 + p f I (fl + fi)drid$ - 2(mr + mp)e2. Jo J-b The inertia matrix of the satellite expressed in the coordinate system 6\ , o 2 , 0 3 becomes I = I B + I P , where I F and Is are the inertia matrices of the panels and the rigid body, respec-tively. 2 .2 .3 K i n e t i c e n e r g y —* —* —t Velocity of any point on the satellite with respect to the body frame (61, b2,63) is given by (—e + / , 0 ,0) , where { fl on the upper panel, f2 on the lower panel, 0 on the main body. Absolute velocity of a generic point expressed in the body frame is vc + (n* x d) + v, 15 where • f-i + f V = d= 0 V o The expression of the kinetic energy takes the form i J{Vc + n* x d + V)2dm = \ J ft? + (ft* x d)2 + V2 + 2VC • (ft* x J) + 2VC • V + 2(& x d) • V]dm. The orbit is assumed circular. Having dropped the constant terms from the above relation, one obtains the following expression for the energy, i n * • Icf2* - -(2mp + mr)e2 + i / f2 dm + fi* • [ dx d dm. 2 2 2 Jp1&cp2 J Here /<f x <Fdm — I dx (—e) dm + I dx f dm J JPlicP2 - / dx f dm Jp1&p2 Jo J-b \ -r,(f1 + f2) J as expressed in the body frame (61,62,63). Therefore, fi* • j dx d = pVl* YtiAiij ~ Mij) f ^i{v)dri j 4>j{c){h + c)d$ -pntJZiAiij + A2ij) / ^i(n)ndn / 4>j{$)d$. J-b JO 2.2 .4 D i s c r e t i z a t i o n p r o c e d u r e Under certain assumptions [2, pp. 254-258], the deformation / of a thin plate is governed by A •! r 1 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<V<b; M^v,l) = i?f(77,0 = 0, -b<n<b; M„(±o,<r) =RTI{±b,s) = 0, 0 <<;</. Reaction, i2, and the bending moment, M, are respectively given by the following expressions: ,d2f d2fs 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 ^ = P4f, f'"(Tb) = f'(Tb) = 0 has a denumerable set of eigenvectors : [—6,6] —>• R corresponding to eigenvalues /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>|| = \Jjib^2dri [3; 4, pp. 295-97]. Furthermore, fb ^/"dr, = - I" rpU/dT, = - / ? y 2 % J —b J-b I i/j/'ip/'dri = - %l)iipj""dr) = -/?y4<5,y. J-b J-b 17 In what follows, ^,'s will govern the transverse variation of the plate deflection. The beam differential equation ^ T = T 7 , /(0) = / ' ( 0 ) = 0, / ' ( / ) = / ' « ( / ) = 0 has a denumerable set of eigenvectors <j>i : [0,/] —• R corresponding to eigenvalues 1A E R. The first two eigenvectors are given below: <£i(f) = £>i(cosh (7!£) _ cosC/yxf) + A; 1(sinh(7if) - s in (7 i f ) ) , 1 , 2.3650204 Dt = -j=, kt = -0.9825022, 71 = j ; 0 2 ( Y ) = •£ , 2 (cosh(72C) - cos(7 2f) + fc2(sinh(72f) - sin(7 2 f)) , 1 , 5.497039 D2 = k2 = -0.9999664, 72 = ; . VI 1 The family (<^ ,) is a total orthonormal system with following properties: / 4>i"<t>j"d$ = / 4>i4>""d$ = ifSij, Jo Jo / <f>i'<J>j'dc = - / (t>i"4>jd^ = - <t>i(f)"d<;. Jo Jo Jo It follows that ipi{i])(j>j{$) is orthonormal and total in the prehilbert space C[[-b,b] x [0,1], R] with L2 norm | | / | | = yjjlQ $b_h pdnde [4, p. 56]. Furthermore, ipi^j's 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], in addition to the present L2 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)). Figure 2-3 show the first four mode shapes tpi(f>i,ipifo,^2^1, and V,2^)2-2 . 2 . 5 S t r a i n e n e r g y o f b e n d i n g a n d t w i s t i n g The strain energy of a plate is given by [7, p. 46] which simplifies to D f l f b fd2f d 2f . 2 j J once / is expressed as Ylk.ij Akijfpi<f>j- Retaining only Ami, Aku, Ak2i, and Ak22 terms, one is led to D ^ , ,6 .0880682 31.285245 2.6194167 2 *" ~ J 2J( £4 + 74 + —£272 )Ak21 k = l 6.0880682 913.60187 4.0366263 2 + ( f4 + J4 + )Ak22 31.285245 2 913.60187 2 1.8141853 "75 ^ f c l l ^ Ji Akl2 ^2j2 Ak21Ak22\-2 . 2 . 6 T h e L a g r a n g i a n The Lagrangian corresponding to the attitude motion of the satellite may be written as L = E/c — Ep — Es. (3) The kinetic energy takes the form Ek = i n * • i n * + Kx + n*2K2 + n*x 3 , m o d e ( 1 , 1 ) Figure 2-3 Plate modes: (1,1) and (1,2). Figure 2-3 Plate modes: (2,1) and (2,2). 21 where 1 2 Ki = -^{2mp + mr)e2 + ^ { A 2 k l l + A2kl2 + i 2 2 1 + i 2 2 2 ) ; k=l K2 = >/pmp'[{A111 - i 2 1 1 )(0.830862 /1 + 0.5749084/) + ( i U 2 - A212)(0.36377/1 + 0.0191105/)]; Kz = - v /pm^6[0.3367362(i i 2 1 + i 2 2 1 ) + 0.1474306(Ai22 + i 2 2 2 ) ] ; e = 2m p + m, •[0.830862(Aiu + ^ 2 1 1 ) + 0.36377(Au 2 + A212)]. The gravitational potential energy is given by equation ( l ) , where / n = / i + ^ ( 4 6 2 + / 2) + 2 ( / l + i ) 2 m p ; J i 2 = -&v/p7rTp-[0.3367362(Ai2i + ^ 2 2 1 ) + 0.1474306(,4i2 2 + A222)]; hz = -y/pmp'iiAin - ^2ii)(0.830862/i + 0.5749084/) + (A112 - A212)(0.36377h + 0.0191105/)]; hz = 0; I22 = h + 2mp(h +1-)2 + mpl2 - (2mp + mr)e2 2 + P + A \ 1 2 + A\2X + A2k22); k=l hz = h + (2m p + mr)e2 2 Ak22)-k=l The expressions Ki, i = 1,2,3, are homogeneous polynomials in Akij. This may be used to write L in the form 1 /n*\ /n*\ 3 7 2 _ 22 where the 11 x 11 matrix n = (hi hi hi o h2 h2 0 aX hz 0 I33 dK3 dX fdK2\T a A ' faK3\T a A a2K1 —IT dA J (hi I12 / 1 3 0 0 0 0 0 0 0 0 \ hi I22 0 K21 —K2i K22 —K22 0 0 0 0 h3 0 I33 0 0 0 0 K31 K31 - K 3 2 K32 0 K21 0 Kn + p Kn K12 K12 0 0 0 0 0 —K21 0 Kn Kn + p Ki2 K\2 0 0 0 0 0 K22 0 K12 Kl2 Kiz + P K13 0 0 0 0 0 —K22 0 K12 K12 K13 K13 + P 0 0 0 0 0 0 Kz\ 0 0 0 0 p 0 0 0 0 0 K31 0 0 0 0 0 p 0 0 0 0 Kz2 0 0 0 0 0 0 p 0 V 0 0 K32 0 0 0 0 0 0 0 p J Here Kn = -0.603317 P m p ; 2mp + mr K12 = -0.3022427—^ ; 2mp + mr K13 = -0.1323286—^—; 2m p + mr K21 = v^ mp"(0.830862/i + 0.5749084/); K22 = v^ mp"(0.36377/i + 0.0191105/); K3i = -0.3367362^/^/^6; i r 3 2 = -0.1474306^^6. 2 3 R e f e r e n c e s [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, 3rd 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, McGraw-Hi 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. 2 4 A P P E N D I X 2 - 1 : E X P A N S I O N O F P O T E N T I A L E N E R G Y Define A — J ( l + x)~? dm, where x = 2o\ • j| + |jf-|2. Application of Taylor's formula to the inegrand (1 + x) -2" yields A = J(1 - X-x + ^x 2 + ^(1 + r))-lx3) dm, where either 0 < rj(x) < x or x < rj(x) < 0. Integrating through second degree terms in the integrand, one obtains A = mr + ^ j"(Ml2 - H°i • <*T) dm + B = mr + (rr(I) - 3 o i • K i ) + B, where < J(^)3[^(4 + I f ) + ^ - I f + ( ^ 2 ) ~ ^ 2 + for a constant K independent of provided < 1. Therefore, „ Gmemr Gme, m , . . „ 1 „,,d„ E, = — ^ + ^ ( - 1 > ( D + 3ol • lol) + - 0 ( ( - | 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 O F A S A T E L L I T E T O E Q U A T I O N S O F M O T I O N We admit that the present state of the world only de-pends on the immediate past. Thanks to this postulate, instead of studying directly the whole succession of phe-momena, 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 equa-tions 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 ap-proaches 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 I n t r o d u c t i o n 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 (aty)- Equations of motion, 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 2 5 2 6 of motion in the desired form, however, involves an arbitrary choice of coordinates, such as Euler angles, on the space of rotations 50(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) 2 7 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-1M-~i63. 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 R8 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((aij),M,A,B), B = g((a.ij),M,A,B), 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 com-puted as a by-product. The Hamiltonian remains constant on a given trajectory 2 8 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 o f t h e M a t h e m a t i c a l T o o l s 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 2 + y2 = 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. Then (Vx) x < 4 and (Vx, y) (x + y)2 = x 2 + y2 + 2xy are, respectively, false and 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) x 2 = 1 is false as both +1 and -1 satisfy the equation x2 = 1. Then, (3!x) x > 0 and x2 = 1 is true. One writes a G A to mean 'a is an element of A\ or 'a belongs to A\ If a does not belong to the set A, one writes a ^ 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 V s e i xe B. Two sets A and B are equal if A C B and B C A. Given two sets A and B, A \ B = {x : x G A and x (fz B} defines the set 'complement of B in A\ AUB = {x : x G A or x G B} defines A union B. AnB = {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 2 = R x iE is the set of points on a Euclidean plane or is the set of all complex numbers. If S1 = {(x,y) G R2 : x 2 + y 2 = 1} then S1 X R C R3 — R x R x i? is the set of all the points on an infinite cylinder of unit 3 0 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 R2\ (d) (0,1) x [0,1] is a neighbourhood of (1/2,1/2) e R2 but not (1/2,0) <E R2. 3 1 radius (Figure 3-lb). (ai , 61) x (<J2,62) is an open rectangle in R2. Similarly, (a\,b\) x (02,62) x («3 563) is an open rectangle in Rz. A set X is open in Rn if V x G X 3{aubx) x • • • x (an, bn) x G (fli,&i) x • • • x ( a n , 6 n ) C X. In other words, X is opeii if for every x in X there exists an open rectangle in Rn, (ai , 61) X • • • x ( a n , 6 n ) , such that x is in the rectangle and the rectangle is entirely contained in X (Figure 3-lc). One may verify that an n-dimensional open rectangle and the n-dimensional open unit ball {x G Rn : \x\ < 1} are open in Rn. A set w C Rn is open in X if 30 open in Rn 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 R2 but is not a neighbourhood of (1/2,0) G R2 (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 i-> F(s), where 5 (->• 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-1{B) = {xeS:f{x)eB}. Consider the function of a complex variable defined by the equation F(z) = z2. Let a < —. This function maps every sector S = {z € R2 : 0 < arg(z) < a} onto a sector F(S) = {z e R2 : 0 < argr(z) < 2a} (Figure 3-2a). The function may be written as F : S -> R2 z i-> 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 ex 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}. This function is denoted by (at)?=i a n d i s called a family («0?=i. 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. 3 4 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 2 , • > a map / : G i — • G 2 is a group homomorphism if V x , y G G ! / ( x - y ) = / ( x ) . / ( y ) and /(e) is the identity element of G 2 , where e is the identity element of G\. The 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 3 5 multiplication table: ,* = f — k2 — -1; ij = k, ji = -k; jk = i, kj = -i; ki = j, ik = -j. This group is not commutative. The set of quaternions Q = {a0 + axi + a2j + a3k : at- G R, i = 0 , . . . , 3} minus the zero element form another non commutative group. The multiplication of letters i, j, k follow the same rules as in the previous example. Similar to complex numbers, one may associate a magnitude to each quaternion: \a0 + iai + ja-2 + ka3\ = \J'a2 + a\ + a\ + o | . The function I | : Q \ { 0 } ^ ( 0 , o o ) , 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; 3 6 • u + (v + w) = (u + v) + w; • 30eV VueV u + 0 = u; • 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 = Rn = { ( x l 5 . . . , x r e) : xt- G R} with operations ( x l 5 . . . , x n ) + ( y i , . . . , y n ) = (x x + y i , . . . , x n + y n ) , a ( x i , : . . , x n ) = ( a x i , . . . , axn). With these two operations Rn becomes a vector space. Given two vector spaces V\ and V2 a map <f> : V\ —»• Vj is a n homomorphism if V a , / ? G i 2 U j i v G ^ i ^(av + = a<£(u) + f3<j)(w). A homomorphism is an isomorphism if it is 1-1 and onto. A subset W of the vectorspace V is a subspace if is a vector space. For example, the 1-dimensional vector spaces {0} x R and R x {0} are two subspaces of R2, respectively, known as the x and y axes. A finite set of vectors v\,..., vn G V is linearly independent if re awi = 0 =>• V i a,- — 0. t = i A finite set of vectors v\,..., vn G V is a 6as*s for V if they are linearly independent and n Vv G V 3(a,-)?=i u = " ' ^ i = i 3 7 For example V\ — ( l , 0,..., 0), v2 = (0,1,0,..., 0), . . . , v n = (0, ...,0,1) form a basis for Rn. 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 a • (a o b) = (a • a) o b = a o (ab), Va, b,c (£ A 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 T ) , denoted by so(3), is a Lie algebra under the above operations defined for n x n matrices. In fact, / o o o \ /0 0 - l \ / o i o \ V x = 0 0 1 , V 2 = 0 0 0 , V s = -1 0 0 \0 -1 0/ \1 0 0 J \ 0 0 OJ form a basis for so(3). Another common Lie algebra is < R3,-,+, 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 Rn to Rm is continuous at x G X if (Figure 3-2c) ViV neighbourhoods of /(x) / 1 ( iV) is a neighbourhood of x in X. A mapping / : X Y is continuous on X if for all x G X it is continuous at x. For example, the function s in( l /x) , i ^ O ; 0, x = 0 is not continuous at x = 0 but is continuous on (0, oo). A mapping / of an open set U C Rn to Rm is called smooth if it has continuous partial derivatives of all orders. A mapping / : X —• Rm of an arbitrary subset X C i E n is called smooth if each point x £ X is contained in an open set U C Rn and there exists a smooth map F : U ^ Rn 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 is a diffeomorphism, mapping £ 0 ( 3 ) C R9 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 50(3) = {A : A a 3 x 3 matrix, A A T = E 3 , det(A) = 1}, where E 3 is a 3 x 3 identity matrix. Let O G 50(3) . Then 50(3) -» SO{3) A H-> O A 3 9 Non-diffeomorphic Protype Diffeomorphic Copies 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. 4 0 is called the trajectory of the model (Figure 3-3b). Many phenomena require geo-metric 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 Rn is called a smooth manifold of dimension m if each x G M is contained in W n M diffeomorphic to an open set U of Rm, where W is an open set in Rn. The diffeomorphism <f> = (4>\,..., (f>n) : W fl M —• U is called a local coordinate system with $t-'s as local coordinates. The inverse of </>, (j>~1 : U —• W fl M , is called a parametrization of the region W fl M of M. The unit sphere S2 = {(x,y,z) G R3 : x2 + y2 + z2 = 1} is a smooth manifold of dimension two. In fact the diffeomorphism {{x,y)eR2 :x2 + y2 < 1} -> R3, (x, y) H-> (x, y, y/l-x2-y2) parametrizes the region z > 0 of S2. By 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 S2, it follows that S2 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, TMX, is defined as follows. Choose a parametrization g:U^MCRk F i g u r e 3-4 (a) Formation of manifolds by pasting; (b) a covering of a sphere by six local coordinates. 4 2 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 Rm. Think of g as a mapping from U to Rk, so that the derivative dgu : Rm as a linear map of Rm to Rk, is defined. Here ( | § r)u is the k x m matrix of first partials, evaluated at u. Set T M j equal to the image dgu(Rm) of dgu (Figure 3-5c). One can show that TMX is an m-dimensional subspace of Rk and is independent of the chosen parametrization. Each element of TMX is a tangent vector to M at z. Corresponding to the coordinate system <j> = (<j>\,... ,<f>m) = <7_1 let (d/d(f>i)x be defined as dgu{ei), where (ei)^L1 is the usual basis for Rm. The family {{d/d(j)i)x)^_1 form a basis for TMX. Consider a curve in 50(3) <f> : R -> 50(3) , fl 0 0 0 !-»• 0 cos(0) - sin(^) \ 0 sin(0) cos{6) passing through the identity matrix E £ 50(3) C R9 at 6 = 0. One may see that / 0 0 0 \ ^/(0) = I 0 0 —1 I £ i 2 9 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 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, ( J i g A * TMX, 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 Tx(0,l). 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 Sl. It is an infinite cylinder, that is, TS1 = S1 x R (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 Rk and N C Rl. The derivative dfx : TMX TNy is defined as follows. Since / is smooth there exists an open set W containing x and a smooth map F : W -> Rl that coincides with / on WDM. Define dfx(v) to be equal to dFx(v) for all v G TMX. One may prove that dFx(v) belongs to TNy and it does not depend on the particular choice of F; this justifies the definition of dfx. The prescription of a velocity vector at each point in the state space is called a velocity vector field. A vector field is 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, Figure 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 + Mv)) = <t>i{v) + 4>2{v), veV; (o^i)(w) = a<f>i[v), aeR. Here is an example of a basis for so(3)*. On so(3) define three linear functionals / i > / 2 > / 3 as follows: h : so{3) R, a i V i + a 2 V 2 + a 3 V 3 i-> a i ; f2 : so(3) -+ R, a i V i + a 2 V 2 + a 3 V 3 H-> a 2 ; / 3 : ao(3) iE, a i V i + a 2 V 2 + a 3 V 3 •-»• a 3 . The linear functionals ( / i ) f = 1 form a basis for so(3)*. This basis is the dual basis to ( V t ) 3 = 1 . The dual of a tangent space to a manifold M at m is called the cotangent space TmM*. The cotangent bundle T*M comprises the union of cotangent spaces, ( J m € M T m M * . The cotangent bundle of a circle is an infinite cylinder. 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 4 7 symmetric: VA l 5 A2 G R Vai ,a 2 , a 3 G V (^Axax + A 2 a 2 ,a 3 ) = Ax^(ai,a3) + A2^>(a2,a3), <j>(aua2) = -<^(a2,ai). 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 / : R2 -» R, (x,y) !-»• x2 + y. Then df = 2xdx + dy is a differential 1-form. That is 2xdx + dy : T(x>y)R2 -* R, (vi, v2) !-)• 2xvi + v2. 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 / A g : Rn x Rn -* 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 R2 by W( X |y) = T(XtV)R2 -> R, (v,w)~(x2 + y)det(Vl V 2 ) . v ' v ' \Wi w2 J Then u = (x2+y)dxAdy is a differential 2-form on R2. This concludes the summary of the mathematical tools which appear in the following sections. 48 3.3 T h e E n e r g y F u n c t i o n F r o m t h e L a g r a n g i a n 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. Def in i t ion . The map FL : TC —• TC*, wc h-» dLc{wc) € T*C is called the fiber derivative [2, p. 209] of L. Lc denotes the restriction of L to the fiber over c G 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) -»• R, 1 3 (a,y), n H + - ^ h (^ »2 + 2n t -a 3 <: + a\i - 3a2u) 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 t y),f2 »-> (aij),(Ii(ni + a3\),h{^2 + 0 3 2 ) ^ 3 ( 0 3 + o 3 3 ) ) . Here so(3) « R3 is the Lie algebra of the Lie group 50(3) , i.e., the tangent space 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 / 0 0 l \ /O -1 o\ V i = 0 0 -1 , V 2 = 0 0 0 L a n d V 3 = 1 0 0 . \0 1 0 J \ - l 0 OJ \0 0 OJ 49 Then, { V i , V 2 , V 3 } is a basis for so(3). The map so(3) -»• R 3 , n iVi + n 2v 2 + n 3v 3 ^ n 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 t - ) ? = 1 . Denote the set of three 1-forms dual to { V x , V 2 , V 3 } by { / i , / 2 , / 3 } . The family ( / ; ) 3 = 1 is called the set of left-invariant differential 1-forms on 50(3) dual to ( V i ) ? = 1 . These 1-forms at E € SO(3) b ecome a basis for so(3)*. Henceforth 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 A : T C —• R , wc »-»• FL(wc) • wc. For our example, A : S O (3) x so(3) R 3 3 (at-y),n i-> X)(n*" + °sO nf = + nf03i)-t'=l 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 3 (ovy), » - A-(n? - al + 3a2u). 5 0 3 . 4 H a m i l t o n i a n F r o m t h e E n e r g y F u n c t i o n 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'1 : SO[3) x so(3)* -> R, (a ty), M "\Y, 7»'[(^ - a3'")2 " °3i + 3 a i i l -Therefore, H=\ £ 3 = i ( ^ - 2M;a 3 t - + 3Jt-o?t.). 3 . 5 T h e C a n o n i c a l 2 - F o r m a n d t h e D i f f e r e n t i a l o f t h e H a m i l t o n i a n i n t h e B o d y C o o r d i n a t e s 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 / M i M 2 M 3 ^ n((5,A),(^,0)) = ^ ( 0 . A - A A ) + def ^ 02 «53 , i=i V ©i 02 ©3 J where {6, A ) , ( 0 , 0 ) G T(( ^(50(3)xso(3)*) « r ( a . y ) SO(3)xso (3 )* and <5X, <52, <53 are the components of 6 G T( a^.)50(3) with respect to frame field { V i , V 2 , V 3 } . 5 1 Next, the differential 1-form dH is expressed in terms of the basis differential 1-forms on T*SO(3). For 1 3 A/f2 1 3 2iW d-ff = - y^ (—Y -^dMj — 2azidMi — 2M,rfa3,- + 6i ,a i ,dai , ) , 2 ,=i J i or more explicitly d-ff = ( - M 2 a 3 3 + M 3 a 3 2 + 3 ( J 2 - / 3 ) a 1 2 a 1 3 ) / i + ( - M 3 a 3 i + M i a 3 3 + 3 ( / 3 - i i ) a n a i 3 ) / 2 + ( - M i a 3 2 -j- M 2 a 3 i + 3(Ji - 7 2 ) a n a i 2 ) / 3 - a3i)dMi — a 3 2 ) r f M 2 - a 3 3 ) r f M 3 Here, the relation d(a ty) = (at-y) / 0 - / 3 / 2 / 3 0 — fi I is used to express dH in V ~h fi 0 terms of independent 1-forms dMi,dM2,dMs, fi,fii fz spanning T* -T*SO(3) at each point ((a t y),M) <E T*SO(3). 3.6 T h e E x p r e s s i o n o f t h e H a m i l t o n i a n V e c t o r F i e l d 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: 3!(«5,A) V(0,6) n{{6,A),{6,G)) = dH{6,e), XH : T*50(3) —> TT*SO{3), (at-y), M i — • ( a i y ) , M , (5, A ) . There exists a unique nonsingular matrix f2 such that n((6,A),(e,Q))= ( | ) - n ( ^ M - E 3 Here, f2 = I j , ^ J , where M = f 0 - M 3 M 2 M 3 0 — M i | and E 3 is an iden-v - M 2 M i 0 tity matrix of order three. Then, the above condition implies f2 (jj) = dH, where s , , ^ : , fdj M' dH = £ 3 = i ( C / t / i + G^AMi) is considered as a column vector ). Therefore, As only A is of interest, i) = n-xdH. A = (-B3,M)dH, M2M3 i.e., A i = ( —— 3 a i 2 a i 3 ) ( J 2 - J 3 ) , J 2-«3 A 2 = ( ^ ^ - - 3 a i 3 a u ) ( / 3 - / i ) , . M X M 2 A 3 = ( - ^ - 3 a u a i 2 ) ( / i - 7 2). i i i 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 — • R2n TqoXoq-1 5 3 one finds {<t>u • • •, <t>n, 0i, • • •, 4>n) = Tqo <}>' = TqoXo<}> = Tq o X o q ~l o qo <f> = ((f>l,...,<f>n,Xi,...,Xn), where qo cf> — ( c^ i , . . . , <j>n) and (Xi,..., Xn) are the local expressions of the integral curve and the vector field. Therefore, <t>i = Xi{4>i,...,<f>n), i=l,...,n if and only if <j> is an integral curve of the vector field X. With these preliminaries, the local expression for vector field XJI is expressed as • , M 2 M 3 Mi = ( — — 3ai2ai3)(h - h), M2 = (MlMl. _ 3 a i 3 a u)(/3 - ^ l ) , ,MiM2 n W r _ . M 3 = ( — — 3 a n o i 2 ) ( / i - h). lil2 3 . 7 S u m m a r y o f t h e P r o c e d u r e 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), f l , 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 in the last section by the example of a rigid satellite in a circular orbit, is applicable to any Lagrangian with dependence 5 4 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 La-grangian L with respect to the generalized velocities. For example, for the La-grangian of a rigid satellite in a circular orbit L{{aij)),n) = Ii{n2 + 2JW- + a\i - Za2u) t = i the fiber derivative of L, the gradient of L with respect to fi, becomes J\L((at-y),n) = (aij), {hiQi + a31),I2(n2 + a32),I3{Q3 + a33)). The second element on the right hand side is referred to as the angular momentum —* M. Therefore M = I(n + o 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 3 A((a t J),n) = FL{(aij),n) • n = ^ (n? + rw) . The energy is given by E = A — L. For a rigid satellite in 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-\{aij),M) = \JZU[(^- - a3i)2 - af,- + 3a?,.]. Z i=i U 5 5 Therefore, 1 3 H/f 2 H = o E ( -r ~ 2Miasi + 3Iiau)-z t=i A Next, dif is computed. The terms dat-y may be expressed in terms of three linearly independent differential 1-forms (/»)3=1. The differential 1-forms (/»)3=1 are related to daij by the following relationship / 0 -fz h \ d{aij) = (ot-y) / 3 0 - h . V-/2 h o ; For the above example dH = - y^ (——-dMj — 2aZidMi — 2M,da3l- + 6audau), 2 i=i I { or more explicitly dH = (-M2a33 + M 3a 3 2 + 3 ( / 2 - 7 3)ai 2ai 3)/i + (-M 3a 3i -j- M i a 3 3 + 3(J 3 - ii)a u a i 3 ) / 2 + (-Mia 3 2 + M 2a 3i + 3(Ji - J 2)anai 2)/3 Mj + (-; a3i)rfMi H M 2 + ( -7 a 3 2)dM 2 J 2 + (-= azz)dM3. H Then, it can be shown that (Section 3.6) M= {-E3,M)dH if dH = YJ3=i + CM-dMi) is considered as a column vector ( ^ ) and 0 -M 3 M 2 M = | M 3 0 -Mi -M 2 Mi 0 The sequence of transitions is indicated below: FL Y E V H dH expressed in body coordinates Equations of Motion M- (~E3,M)dH 5 7 3 . 8 A p p l i c a t i o n o f t h e P r o c e d u r e t o a S a t e l l i t e W i t h F l e x i b l e 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 R8 x so(3) x R8 -> 50(3) x R8 x so(3)* x R8*, ( a t y ) , A,n , A ^ (a t J),A,in* + K 2 V ^ 3 y dA The action takes the form A : 50(3) x R8 x so(3) x Rs R, [aij),A,n,A ^ FL(n,A).(n,A) ^ / 0 H-> rl • in* + n. K2 \ + 2 K l + n*2K2 + n*K3, where dA dK2 --^-•A = K2, dA dA The energy £ = A-L in • in - y o 3 • io 3 + + n • | K2 \ +EP + ES, o # 2 # 3 may be written in the form ( 0 ^ 2 \} + EP + EA. 5 8 The Hamiltonian, then, becomes where (§) = n( Hence d H =U) • n U ) + 2 U) • r f ( n >U) - ^ • m - 3 - -7003 • M + S^doi + - 7 2 o 1 . d I o 1 + dE8 (1) in terms of the 1-forms /i,/2, fz,dA,dM, dB on T*(50(3) x i Z 8 ) , which form a basis o f Trr * r n ^ r * ( 5 0 ( 3 ) x R 8 ) a t e a c h P ° i n t ( (0 .7) , e T*(50(3) x fl8). That is, 3 8 3 8 dH = Y, Cfifi + E C * i d A i + CMidMi + Y, CBidBi. (2) i=i i=i t'=i t=i Analogous to the case of the rigid satellite / M 0 - E 3 0 \ 0 0 0 E 8 E 3 0 0 0 ' V 0 - E 8 0 0 J These adjoined with a choice of kinematical differential equations q = K(q)Q and the equation G H - ' < ) - ( 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. 5 9 3 - 9 T h e S o f t w a r e R E D U C E i n D e r i v a t i o n o f t h e E q u a t i o n s o f M o t i o n 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\ — CAX- To find C x x , one sets 6 = A = 0 £ R3,X = ex = (1,0,0,0,0,0,0,0), and A = 0£ R8 in dH{ ^ ^ { 6 , A, A , A). From equation ( l ) , one obtains C^ = 'rf(E_1)l^ (f)+ h2*1 • + d E ^ = h -On the other hand, din dhz dL 22 5/33 where, from Section 3-8, d I\2 = K3idA\2\ + K3\dA22\ + -K^32^122 + K32dA222, dlis = —K2\dAm + K2idA2ii — K22dAn2 + K22dA212, dL 33 dL 22 2dA '11 + P Ku K i 2 Ki2 0 0 0 0\ Ku Ku + P K\2 K\2 0 0 0 0 K12 K\2 K13 + p Ki3 0 0 0 0 K12 K\2 K13 K13 + 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 6 0 As a result dh2{ex) = 0, dli3{ei) - -K2i, dl22{ex) = 2 fKu+p\ Ku K12 V K12 On simplification 31.28245 Z4 - D A 111 + 3 72 [ K\2 v #12 y ^211 [a\2 + a2^) - K 2iaiiOi 3], where -i> l A=e i —•ft 21—— 1-di 13 di 22 33 Now a symbolic manipulation code proves useful • to invert TI symbolically; • to find the partials of I I - 1 with respect I\z,I22, and 7 3 3 . In fact, it is more expedient to work with the adjoint of II, |ITjIT - 1, as the entries of an adjoint are devoid of a denominator. Then a ( n - i ) _ ^ | n | - a r f n ^ ami in which operations 81 12 and ai 12 |IX|2 are performed symbolically. A further numerical calculation yields the entries of ^ • 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. 6 1 R e f e r e n c e s [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, 2nd Edition, 1980, pp. 208-23, 311-16. [3] Arnold, V. I . , Ordinary Differential Equations, The M I T Press, Cambrige, Mas-sachusetts, 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 associ-ated representation in attitude synamics," Acta Astronautica, Vol . 15, No. 11, 1987, pp. 833-43. 4 . N E W L I B R A T I O N A L G E N E R A L I Z E D C O O R D I N A T E S ( ^ - - P A R A M E T E R S ) 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 rep-resent a globally defined nonsingular rational parametrization of space of rotations, suitable for problems of dynamics involving general rotations. 4 . 1 I n t r o d u c t i o n 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. How-ever, 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 t J ) with deter-minant one. The matrix (a,y) is also referred to as the direction cosine matrix. The 6 2 6 3 collection of such matrices 50(3) is precisely the space of all orientation preserving isometries of R3 with a fixed point. One of the objectives of most librational investigations is to study the time evolution of (aty) as a function of initial conditions and system parameters. The time evolution of (at-y) is summarized as a second order differential equation [l] on 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 50(3) , which yield a convenient and nonsingular set of differential equations. The Euler angles and their variants [2], twenty four in 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 T50(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 inspec-tion 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 t J ' s in terms of some set 6 4 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) as an answer to (a) and (b). To 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 details of the above interconnection are elucidated. 4 . 2 F o r m u l a t i o n Given a unit quaternion [4] g, the corresponding element of SO (3) under the 2-1 Lie group covering (i.e, a continuous Lie group homomorphism whose derivative at 1 6 5 3 is an isomorphism) c: S3 —• SO (3) is given by ?o + 9 i - ° 2 - 9 3 2(gig 2-g 3g 0) 2(gig3 + g0g2) \ 2(gig2 + g3go) go - °l +l\ ~ Qz 2(g2g3 - g0gi) . (1) 2 ( g i g 3 - 9 o ? 2 ) 2(g2g3 + g0gi) g 2 - q\ - ql + 9 3 J The 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 S3 and the space of orientations 50(3) . The 3-sphere S3 is now parametrized by the inverse of the stereographic projection p-'-.R3 - 5 3 \ {1} , p - J (a ) = (|a |2 - l,2au2o2,2<73)/(l + \a |2), where a = (ai,o 2 ,cr 3) and \o |2 = o\ + o\ + o\. Therefore, cop-l-.R3 -* 50(3) , eop-1[a) = {aij) 65 q = ( q , fl, fl2 fl3) 2-1 F i g u r e 4-1 A diagrammatic transition from <7-parameters, in R3, to unit quaternions, in R4, to the space of rigid motions with a fixed point, in R9. 66 is an analytic local diffeomorphism (i.e., diffeomorphic on an open neighbourhood of each point), 1-1 on {x G R3: \x\ < 1} (Figure 4-1). Here u = ( | a | 2 + l ) 2 ; a n = l + 8(o-2 - |a|2)/u; a 1 2 = 4(2aiO"2 - a3(\a | 2 - l ) ) / u ; al3 = 4(2<7io-3 + o-2(|a | 2 - l ) ) /u ; a 2i = 4(2<TIO-2 + o-3(\a | 2 - l ) ) /u ; a 2 2 = l + 8(o-2-|a|2)/u; (2) 023 = 4(2a 2a 3 - o-i(|a | 2 - l))/u; a 3i = 4(2o"io-3 - o-2(|or | 2 - l ) ) / « ; a 3 2 = 4(2o-2o3 + oi{\8 | 2 - l ) ) /u ; a 3 3 = l + 8 ( a | - | ^ | 2 ) / a . Note the absence of transcendental functions in the entries of (a ty). Given G 5 ' 0 ( 3 ) \ { E } ) there are precisely two elements (corresponding to antipodal elements of S3) which transform to (a ty). The origin 0 G R3 is the only element which is transformed to E . Using (2), one obtains (Appendix 4-II) = B ( j , ( 3 ) \a3) \ n 3 J where — 4 B in terms of s is a I - al - o\ + 1 2(o-icr2 + 03) 2(o-io-3 - cr2) 2(o"io-2 - o3) -o\ + o\ - o\ + 1 2(o-i + o-2o-3) 2(o-ia3 + cr2) 2(-ai + o-2o-3) -o-2 - a| + o-| + 1 3 (B) = > 1. Hence, the proposed parametrization has no singularity (a point 6 7 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 € Rn; (j^=h(x,v,n,(aij),t), heRn+3. (4) —* Here h is an explicit function of time if some of the appendages describe a pre-determined motion with respect to the satellite or if the environmental forces are modeled in a time dependent fashion. 4.3 I n t e r p r e t a t i o n , P r o p e r t i e s , a n d G e n e r a l i z a t i o n o f t h e ^ - P a r a m e t e r s Let q E S3 be written as (cos(/?/2), n i sin(/?/2), n2 sin(/?/2), n 3 sin(/?/2)). The orientation q is obtained by a rotation of /? about the unit vector n = (ni , rc2,ra3). The application of stereographic projection, p:S3\{l}->R3, P{q) = ^^i, to the components of q yields: = * i s i n ( | ) / ( l - cos(2-)) = nxcot ( - ) ; = n 2 s i n ( | ) / ( l - = n 2 cot(^) ; 4 0-3 = n 3 s i n ( | ) / ( l - — 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 1 \a | about — a . Given two successive rotations E and a, the composition is given by or equivalently | a + E | 2 ((1 - lE| 2)a+ (1 - \o\2)t -23 xt) ( | a | 2 | £ | 2 - 2 £ . a + l ) Note that a and —a/\a\2 represent the same orientations as they are mapped to antipodal elements of unit quaternions which in turn transform to a single element of 50(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 S3. Let M be an orthogonal matrix of order 4. Then p o M o p - 1 : ^ 3 50(3) , poMop-1^) = (a,y) is another oMike parametrization. For example, if M / - l 0 0 0\ 0 1 0 0 0 0 1 0 V 0 0 0 lj 6 9 one obtains the following: = l - 8 ( | a | 2 - a 2 ) / « ; «12 = 4(2o 1 a 2 + < 7 3 ( | a | 2 - l ) ) / u ; ^13 = 4 ( 2 c T 1 o - 3 - o 2 (| tT |2 - l ) ) / u ; a 2 1 = 4(2(7^2 - <73(|a | 2 - l ) ) / u ; a 2 2 = l - 8 ( | a |2 - a 2 ) / u ; «23 = 4(2a 2 a 3 + o - 1 ( | a | 2 - l ) ) / u ; = 4 ( a 2 ( | a | 2 - l ) + 2 0 - l 0 - 3 ) /u; a-32 = 4 ( ^ ( 1 - | a | 2 ) + 2a2a3)/u; 0-33 = l - 8 ( | o | 2 - a 2 ) / « . Though all these parametrizations are nondegenerate, algebraic, and minimal in the number of components only two, corresponding to M / l 0 0 0\ 0 1 0 0 0 0 1 0 Vo 0 0 lj (-\ 0 0 0 , projection from the North pole, and M 0 1 0 0 0 0 1 0 V o o o i yield a convenient geometric interpretation. This follows from the computation of p o M(g) for an arbitrary unit quaternion a = { cos{^)^ ni sin{^)^ n2sin{-),n3sin{-)). In the latter parametrization, o represents a rotation of 4 tan 1 \a\ about a. This parametrization was also discovered, independently, by Milenkovic [5]. 7 0 4 . 4 T w o I l l u s t r a t i v e E x a m p l e s To assess relative merit of a-parameters, the Euler angles (cti, a2,ctz; Equation 5b), and quaternions (qo, qi, g 2 , qz\ Equation 5c) a typical example of a rigid satellite in a circular orbit is considered. The governing equations of attitude motion are as follows: fi* = (fi*n*: _ "&a\ta\3)T\\ ft2 = (ft^ft* _ 3°i3aii)r2; (5) 03 = (ft*ft2 _ 3 a n a i 2 ) l 3 ; (5a) ° 3 1 with fi = fi* — I a 3 2 ) , where B was given in Equation (3). The nondimensional a 3 3 . intertia parameters T i , T 2 , Tz are denned as ^J^3, l3yJ1, I^J^ > repectively. 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 coz/co2 —siz/co2 o\ siz co 3 0 j fi, (56) -cozsi2/co2 sizsi2/co2 1J where cot- = cos(o;t) and = sin(ai), i = 1,2,3. The direction cosine matrix (aiy) is expressed in terms of the Euler angles as [2], co2coz —co2siz si2 \ siisi2coz + sizcoi —siisi2siz + sizco\ sixco2 I . —coxsi2coz + sizsii coisi2siz + cozsi\ co\co2 J 7 1 Kinematical differential equations for quaternions have the following form [3] /9o\ 9i 92 f -qi -q2 - q 3 \ 9o - 0 3 92 03 9o — 9i fl. (5c) V93/ V-92 9i 9o / 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 - 2n*a3i + 3a2,-), t'=i was utilized as an absolute measure of error in comparing the above three parametriza-tions. As shown in Appendix 4-III, E = 0 and hence E is a constant on trajectories of Equation (5). Wi th 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 compara-ble 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 satel-l 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 sig-nificantly 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 coor-dinates. 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 7 2 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 R e p r e s e n t a t i o n o f L i b r a t i o n a l T r a j e c t o r i e s A concise and effective method of visualization, which clearly portrays the essential features of a given phenomenon, is naturally conducive to further under-standing of the underlying principles. The following theorem by Euler is the main-spring 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 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 R3 : \x\ — 1} represent the same element of 50(3); the space thus obtained, homeomorphic -4atan-1 (|a|)/(|a|7r), a 0; 0, a = 0. 7 3 A n y r i g i d mot ion a b o u t a f i x e d po in t is a r o t a t i o n a b o u t some axis. axis o f r o t a t i o n "amount o f " r o t a t i o n in f r ac t i ons 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. 7 4 Longitude= 45. f=(-0.1 , 0.3 ) Latitude= 35.27° c={ 0. -0.026, 0. ) Position Scale= 10. f=( 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. 7 5 to 50(3) , is known as RP3 (Real Projective 3-space [6]). Isometric view of the o i , o 2 , o3-frame presented here is that of an observer stationed at longitude=45° and Lati tude=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 2 axis. Similarly, point B with coordinates (0,0.5,0) denotes 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 2 axis. 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. As 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. At orientation C (the initial orientation) the angular velocity has only a — o3 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 7 6 Longitude= 0? f=(-0.1 , 0.3 ) Latitude= 90? *=( 0. -0.026, 0. ) Position Scale= 10. f=( 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. 7 7 Longiiude= 90. Latitude= 0? Position Scale= 10. 0-f=(-0.1 , 0.3 ) o-{ o. ,-0.026, 0. f =( 0. , 0. -0.05) 4 0 3 F i g u r e 4-5 Time evolution of the orientation of a rigid object seen from the side vantage point. 7 8 Longitude= 45? Latitude= 35.27° Position Scale= 10. Velocity Scale= 4. f=(-0.1 , 0.3 ) *=( 0. -0.026, 0. ) f =( 0. , 0. -0.05) Figure 4-6 Time evolution of the orientation and angular velocity of a rigid object. &={ 0. -0.026 , 0. ) f-( 0. , 0. -0.05) o2 7 9 <°i f =(-0.1 , 0.3 ) 2D01 ©o(3) \ © 2001 \ ©o(3) 1 © 2DD1 ©o(3) © 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 3 axis (opposite orbit normal). 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\, T2. The orientations C and D are shown in the second subplot 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 ty), we plot the point — a in the closed unit ball xG R 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 in fractions of TT about r. The advantage of this representation lies in 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 4-IV) . 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 in the deflection of the a(0)=( 0 . ° , 6 . ° , 0 . ° ) f=(-.l,.3) f = ( 0 . , 0 . , - 0 5 ) o 6 1 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. -0.026, 0. ) Position Scale= 10. f=( 0. , 0. -0.05) • 0 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. 8 3 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 r = 400 Kg A n i = 0 m = 10 Kg A 2 n = 0 1= 7.3 m A n 2 = 0 h= 0.7 m A 2 1 2 = 0 b= 0.6 m A 1 2 1 = 0 D= 0.003 N m A 2 2 1 = 0 Ij= 1071.3 Kg m 2 A 1 2 2 = 0 I2= 4761.6 Kg m 2 A 2 2 2 = 0 I3= 5500 Kg m 2 Orb. Period- 120 min Initial Conditions Initial Time Rate of Rigid Core of Amplitudes (m/s) a = 0 A i n = 0.001 f=0 A 2 1 1 = 0 A I 8 . = 0 ^ 2 2 1 ~ ^ A = 0 ^ 1 2 2 U ^ 2 2 2 = 0 Table 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. Figure 4-10 Orientation of the hybrid satellite in the a-parameters Figure 4-11 ( a) 8 7 F i g u r e 4 - 1 1 ( b ) Upper and lower panel plate mode amplitudes: mode (1,2). 88 Figure 4 - 1 1 (c) Upper and lower panel plate mode amplitudes: mode (2,1). 9 0 R e f e r e n c e s [l] Lang, S., Differential Manifolds, Springer-Verlag, New York, 1985, p. 90. [2] Kane, T., Likins, P., and Levinson, D. , Spacecraft Dynamics, McGraw-Hil l , New York, 1983, Appendices I and II. [3] Williams, C . J .H . , "Dynamics modelling and formulation techniques for non-rigid 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, Springer-Verlag, 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 Geometry, Publish or Perish, Inc., Berkeley, 1979, Vol . 2, p. 84. 9 1 A P P E N D I X 4 - 1 : D I R E C T I O N C O S I N E S A N D Q U A T E R N I O N S This appendix starts with a construction of the (universal) covering c : S3 —• £ 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 R3 with {q £ Q : qo = 0}. Given q £ 5 3 , Cq : Q —> Q, Cq(a) = qaq is an orthogonal transformation, as |C g (a) | = \qaq\ = \q\\a\\q\ = \a\. Here q = (q0, — <7i, — qi, — 93). Cq leaves R C Q invariant. Hence, it leaves the orthogonal complement R3 invariant. Therefore, Cq restricted to R3, cq, is an orthogonal transformation isotopic to identity (i.e, 3H : [0,1] x R3 —• R3 smooth H(0, •) = c g , H(l, •) — id, Vi £ [0,1] H(t, •) is a diffeomorphism) as S3 is pathwise connected (i.e., Va,6 £ S3 3<f> : [0,1] — • S3 smooth <f>(0) = a,<f>(l) = b). Hence cq is 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 cq appears as ( 0.1 + 9.1 ~ ?2 ~ ?3 2 ( ? i 9 2 - 9o?3) 2(oiff3 + q0qi) \ {aij) = 2 ( ? ! 0 2 + qsQo) 9o - 9 i + ?2 - ol 2 ( o 2 g 3 - g 0 9 i ) • (I ~ 1) V 2 ( 9 ^ 3 - q0q2) 2{q2q3 + q0q{) ql - q\ - q\ + q% J It is obtained by computing cq(i),cq(j),cq(k). Now cqp{a) = qpa{qp) = qpapq = cqcp(a), q,p£ Q Hence c : S3 50(3) , c(q) = cq 9 2 is a Lie group homomorphism. Furthermore, from Equation (1-1) it follows that: 2go = +V {an + 0-12 + 033 + 1); 4gi = («32 - a 2 3)/g 0; 4q2 = (ai3 - a31)/q0; 4q3 = (a 2i - al2)/qo-Therefore, b : 50(3) -» g, 6((a,y)) = (q0,9i, 92,9s) defined with the positive choice of q0 is smooth (i.e., infinitely differentiate) about 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 cq,cq(q — q0) = q(q — qo)q~ = 9 — 9o 5 i-e., g — go G ii! 3 is an eigenvector of cq €E 50(3) with eigenvalue 1. Suppose go 7^ l(c<j 7^ Id). Define / = (g — go)/y/(l — q2). Complete J to a right-handed orthonormal basis J, J, K of Rz. Then one calculates: jJ = (2g2 - 1)J + 2qQyj{l-ql)K; cqK = -2q0y/{l-q*)J + (2g* - l)X. (\ 0 0 Therefore, in basis /, J, if, c g reads 0 cos j3 — sin /? J , where /? is defined y 0 sin (3 cos /? by V ( l - * 0 ) = cos(/V2). 9 3 A P P E N D I X 4-II: A N G U L A R V E L O C I T Y A N D C T - P A R A M E T E R S Consider a motion m : R —• 50(3) in the configuration space. Interpret ( ^ t ' ) 3 = i a s the components of m in the left-invariant frame field [9; 8, p. 323] ( V i ) , i = 1,2,3, generated by 0 - 1 0 1 0 0 | G TIdSO{3). 0 0 0 To express the natural basis (d/dcr,) in terms of ( V i ) , (a ty)Tt9(a,y)/c9cT f c, k = 1,2,3, are evaluated as follows: 0 2 ( a i t T 3 - a2) - 2 ( t 7 i O - 2 + 0 3 ) 2 ( - c T X a 3 + 0 2 ) 0 t r 2 - t r 2 - t r f + 1 2(CTICT2 + 0-3) ~o\ + t r | + al - 1 0 0 2(ai + c r 2 c r 3 ) o\ - o\ + a\ - 1 - 2 ( a i + t T 2 a 3 ) 0 2 ( t7 i tr 2 — t r 3 ) - a 2 + a | - a l + 1 2(-o- 1 a 2 + t r 3 ) 0 0 - a \ - a \ + al + 1 2(tT! - tr 2 o- 3 ) o\ + o\-o\-\ 0 2(tTitT 3 + a2) 2(-CTI + a 2 c r 3 ) - 2 ( t r x c T 3 + t r 2 ) 0 Components of these matrices in the basis (v t) yield the columns of the following matrix 4 (-al + ol + al-l 2(-o-!a 2 + c r 3 ) -2{axa3 + a 2 ) B _ 1 = - - 2 ( t T i O - 2 + t r 3 ) erf - a\ + trf - 1 2(ai - CT2CT3) u I 2{-oxoz + a 2 ) - 2 ( 0 - ! + c r 2 c r 3 ) a 2 + CT| - a\ - 1 9 4 A P P E N D I X 4 - I I I : C O N S E R V A T I O N O F E N E R G Y j 3 2 . t=l 1 3 ^ E M n f - ^ + Sa^-ai,.] 2 . t=i 3 2 1 3 - ^ / t ( n 2 + 3 a l t - a 3 i 2 ) . t=i Therefore, E = U(tiiZli + 3a l t a l t - - a3ia3i). (Ill - 1) t'=i o - n 3 n 2 Using the identity [8, p. 323] (a,y) = (a,y) I Q3 0 —H i I , and substituting n 2 H i 0 from Equation (5) in (III-l), one obtains E = 0. 9 5 A P P E N D I X 4 - I V : C O N F O R M A L I T Y O F cop - l First the conformality of p _ 1 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 - 1 at the point of intersection. Tangent vectors along the increasing directions of 01 5 02 j o"3 are given as follows: 3 do2 da? 2 u 2 u u (l-a\ + al + al\ —2o2a\ -2a3ax V 2<7i J -2axo3 \ 1 + aj - a\ + a2 —2a3a2 2a2 J -2o\a3 \ —2a2a3 l + a2 + al-al 2a3 J Defining g,y =< d/dai,d jdai >, where < , > is the usual dot product, one obtains { g i j ) = - E . u Therefore, p _ 1 is conformal [9, p. 334]. The function c is also conformal as the following argument outlines. Let Z\,Z2,Z3 be the left invariant frame field on S3 generated by orthonormal tangent vectors i,j,k at 1 6 S 3 , where i = o \ 1 0 0 , j = 1 oJ k = Consider c as i24\{0} —> R9, where the space of 3 x 3 matrices is viewed as a subset of R9. The linearization of c at q is dc: R4 —• R9, dcq(v) = Lv , 96 where 90 9 l - 9 2 - 9 3 ^ - 9 3 92 9 i - 9 o 92 93 9o 9 i 93 92 9 i 9o 9o - 9 l 92 - 9 3 - 9 i - 9 o 93 92 - 9 2 93 - 9 o 91 9 i 9o 93 92 I 9o - 9 i - 9 2 93 J Observing that Zi, i = 1,2,3, at q £ S3 are expressed as - 9 3 9o V 91 J respectively, one finds 0 0 LZi LZ, LZo = 2 2 ( - 9 o 9 i + 9293) "9o + 9? - 92 + 9s 9o2 9 2 - 9 i + ql -2(g 09i + 9293) -9o + 9 i + 9 2 ql 2(9o92 + 9i93) 2(9o93 - 9i92) -2 ( 9 0 9 2 + 9i93) 2 ( ? o 9 i - 9293) 0 0 0 9o + 9? - 92 - 9 l 2(g 093 + 9i92) 2 ( - g 0 9 2 + 9i93) 2 ( - 9 o 9 3 + 9i92) 9o ~ 9 i + 9 2 - 93 2 ( g 0 9 i + 9293) "9o - ql + 92 + 93 - 2(9o93 + 9i92) 2(go92 - 9i93) = 2 t / i ; 2 v 2 ; = 2v3. 0 0 0 Hence, c is conformal with respect to metrics in which (Zi) and (V,) are orthonormal families [9, p. 334]. Consider two intersecting trajectories m,- : R —»• R3, i — 1,2, with mi(0) = m2(0). Define yt- = c • p _ 1 • m;. From conformality of p - 1 and c it follows that Ang(rhi(0), 7722(0)) = Ang(yx(0),y2(0)). Here, Ang( , ) is the measure of the angle between the vectors appearing inside the parentheses. Express the two vectors yi (0) and y 2 (0) in terms of the their components with respect to the frame field (V,) (Appendix 4-II) as yi(0) = 91V1 + 02V2 + O 3 V 3 and y 2 (0) = G 1 V 1 + G 2 V 2 + G 3 V 3 . Then, g = (91,92,93) and G = (Gi,G2,G3) denote the angular velocity vectors of 9 7 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 + 0 2 V 2 + gzV3 (t7i,g2,g3). Thus, A ? i o ( m i ( 0 ) , m 2 ( 0 ) ) = Ang{g,G). 5 . U S E O F T H E E N E R G Y F U N C T I O N I N D E T E R M I N I N G B O U N D S O F S T A B I L I T Y During the motion of a mechanical system the 2s quan-tities 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 conditions. L . D . Landau In 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 non-symmetric 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 perturba-tions 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 equi-librium as a function of disturbance energy for a prescribed configuration. 5 . 1 I n t r o d u c t i o n A n 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 9 8 9 9 that the dynamical system of a rigid body without axis of symmetry in circular orbit possesses 24 equilibria in 50(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,b2, and 63 are respectively expressed by T\ = 0, T2 = 0, and T\ + T2 = 0. Almost all configurations have no axis of symmetry, i.e., the subset of symmetric configurations is of measure zero. Theo rem. One of the six dynamically distinct equilibria is stable for any asym-metric 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 glob-alizing 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 differ-ential 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 parame-ters 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 infor-mation 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 1 0 0 of imparted energy) to the satellite before it swings away from the equilibrium. The results are further corroborated by equipotential contour surfaces on the configu-ration 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 F o r m u l a t i o n Euler's equations yield the following nondimensional attitude equations of mo-tion for a rigid object in a circular orbit: fi* = (HjOg - 3 a 1 2 a i 3 ) r l 5 n^ = ( n^n * - 3 a 1 3 a n ) r 2 , (i) fig = (n x£} 2 — 3a i i a i 2 )T3 , f a3A ( o - n 3 n2 \ where fi = Cl* — a 3 2 . Using (dty) = (a,j) I f l 3 0 — f)i [5], one may \a33 J V ~ n 2 Hi 0 J write equations (1) in terms of Cl: n x = (f2 2n 3 + a 3 2 a 3 3 - 3a 1 2 a 1 3 )r x + a 3 3 n 2 ( T i + 1) + a 3 2 f2 3 (T' 1 - 1), f22 = (n 3n! + a 3 3 a 3 i - 3 a 1 3 a n ) r 2 + a31n3{T2 + 1) + a 3 3 n!( r 2 - 1), (2a) f23 = (n xn 2 + a 3 1 a 3 2 - 3 a n a 1 2 ) r 3 + a 3 2ni(r 3 + 1) + 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 1 0 1 at hand: ^fflj=B^naj, (26) where — 4 B in terms of a is given by o \ - o \ - o l + \ 1(0102 + 0-3) 2(0-103 - o-2) A 2(o-xO-2 - 0-3) - 0 2 + 0 2 - 0J + 1 2(o-i + 0-2+0-3) . 2(o-io-3 + a2) 2(-<7i + a2a3) -al - a\ + o\ + 1J The above choice is derived from the following coordinates on SO (3). Here a = (01 5 o"25 03)5 \o 12 = al + a\ + o | . The parametrization i23 —• 50(3), a 1 — • (o t -y) is an analytic local diffeomorphism, 1-1 on {x G R3: \x\ < 1} with u = ( l * | 2 + l ) 2 , an = l + 8 ( a 2 - | a | 2 ) / « , a i 2 = 4(2oia2 - az(\a | 2 - !))/«, = 4(2o - ia 3 + a 2 ( | a | 2 - !))/«, a 2 i = 4(2oia2/u + 0-3), « 2 2 = l + 8 ( 0 2 2 - | o | 2 ) / u , ^23 = 4(20203 - 0"i(|0 | 2 - !))/«, <*31 = 4(20103 — o-2(|0 | 2 — !))/«, 032 = 4(2a2az/u + 01), a 33 = l + 8(0 3 2 - | 0 | 2 ) / u . The system of equations (2) is the local expression of a 2-parameter vector field (on TSO{3)). The inequality \I2-h\ < h implies - 1 < Tx < 1. Similarly, - 1 < T, < 1, i = 102 Figure 5 -1 Graph of the nondimensionalized inertia parameter Tz as a func-tion of T\ and T 2 . 1 0 3 2,3. The nondimensional inertia parameters (T t) satisfy the relation T\ + T2 + T 3 + TxT2Tz = 0; this relation defines a surface shown in Figure 5-1. The parameter space is defined by —1 < T, < 1, t = 1,2. Given T\ and T2 the inertia param-eters Ii,I2,l3 are determined up to multiplication by a constant. Equations (2) yield the following relations, satisfied by the equilibria: 0 = 0, (3a) 032^33 — 3 a i 2 a i 3 = 0, 033031 — 3 a 1 3 a n = 0, (36) a 3 i a 3 2 — 3 o n a i 2 = 0, provided TiT2T3 ^ 0. The equation (3b) has exactly 24 distinct solutions in 50(3) . These are partitioned in to six "dynamically distinct" classes of four ( e»j ') i=i, y=r In Table 5-1, the body axes lying on the 01,02-plane are drawn next to each equilibrium element (aty) G 50(3) . To indicate the equilibria in the RP(3) 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 ety are stable for some t provided T\T2T3 ^ 0. 104 / 1 0 0\ ( 0 1 0 ) V 0 0 i / »>i| (o.o.o) / l 0 0 \ ( 0 - 1 0 ) Vo 0 -l) ( 1 ' ° ' 0 ) j b2 / - l 0 0 \ ( 0 - 1 0 ) V o o i) k 1—-(0,0,1) b 2 / - l 0 0 \ ( 0 1 0 ) V 0 0-1) (0,1,0) \ ] / l 0 0 \ ( 0 0 1 ) \o-io) 2 2 b, . (v / 2 - 1,0,0) J b 3 1 / 1 0 0 \ f 0 0 - 1 V o i o y 3 2 f 1 (1-V ,^0,0)j ^ / - l 0 0 \ ( 0 0 - 1 ] V o - l o ; ( 0 , - ^ , f ) b3 lb] ( 0 , - ^ , - ^ 3 / 0 0 1\ ( 0 1 0 ) V - i o o ; 2 V b 3 i ( 0,l-v/ 2 , 0 ) 2 1 /O 0 1\ 0 - 1 0 \ 1 0 0 / / 0 0 - 1 \ ( 0 - 1 0 ) wo o ; - &k ( f , 0 , - f ) ^ b 2 /O 0 - 1 \ ( 0 1 0 ) \10 0 ) 2 2 J (0,x/2-l ,0) b 2 / 0 1 0 \ ( 1 0 0 ) V o o - i ; b2 b!- 1 V 0 0 1 ) ( 0 , 0 , ^ - i ) J B I / 0 - 1 0 \ ( - 1 0 0 ) V o o - i ; 2 * 2 J *b2 /O - 1 0\ ( 1 0 0 ) Vo o i) (0 ,0 ,1 -v^ jb i I / 0 1 0 \ ( 0 0 1 ) V i o o y / 0 1 0 \ ( 0 0 - 1 ) V - i o o ) [\ /O - 1 0 \ ( 0 0 1 ) I ( l + i - J + f c ) (1,-1,1) 1 b 3 / 0 - 1 0 \ ( 0 0 1 ) V - i o o; »(l-,-+J + fc) (-1,1,1) bT"] V o i o ; (1,1,1) ^ 5L 1 V 0 - 1 0 J f-1,1,-1) |t>3 Li / 0 0 - 1 \ f - 1 0 0 ) V o i o ; (1,-1,-1) 1 b, lb3 /O 0 - 1 \ ( 1 0 0 ) V o - i o ; i ( i - , - -y + ft) (-1,-1,1) ^ — i b3* Table 5-1 The list of equilibria in 50(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 Proof : ^ (n , (a,ij)) is a first integral for equation (2), i.e., E = 0 on trajectories of (2) [1]. By 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£p((a t-y)) for local minima. The critical points of Ep are precisely the equilibria (e,y), provided T1T2T3 ^ 0 (Appendix 5-1). Ep:SO(3) —* R remains invariant under the action of the subgroup {eij)4=1. Hence, the behavior of Ep about the equilibria within a given row, in Table 5-1, is identical (up to a diffeomorphism). Therefore, it suffices to examine the family (e*'i)f=i f ° r local minima (strictness is a consequence of (e,y) forming a discrete subspace of SO(3)). The proof is concluded via three distinct methods. These methods give pro-gressively more detailed information about the local behavior of Ep at equilibria. (i) The analytic (in particular, continuous) function Ep has a minimum on the compact (i.e., has an open complement and is bounded as a subset of R9) 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, Ep is a Morse function (i.e., all the equilibria are nonsingular). Morse inequalities [6] applied to this prob-lem yield V{ > rii, i = 0,1,2,3; —v0 + vx > -n0 + rai; vo - vi + V2 > n0 — n\ + n 2 ; (4) -VQ + v\ - v2 + v3 = -nQ + nx - n 2 + n3. Here n ; = dimHi(SO(3), Z2) = 1, * = 0,1,2,3 [8]. The variable denotes the 1 0 6 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 Ep at the critical point. This number is independent of the local coordinates [6]. The n,'s reflect some of the topological properties of 50(3) . Equations (4) simplify to the following: (a) Vi>l, i = 0,l,2,3; (b) - v0 + V l > 0; (c) v0 - vi + v2 > 1; (5) (d) - v0 + vi - v2 + vz = 0. 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 Ep about the equilibria (e,-y) (Ap-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\T2(T\ + T2) = 0. Index of a function at an equilibrium gives information on the behaviour of the function close to the equilibrium. For example, at the equilibrium en (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 1 0 7 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. By 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 where E0 is the energy at the equilibrium. The variation of the index of the energy function at the equilibrium e 2 i as 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 2 i . Moving clockwise about the point T\ = T 2 = 0 in the parameter space, one sees that the sequence of transitions is the same as those of the equilibrium en : 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 equi-librium can have an index varying between zero and six, the index of E on the phase space 50(3) x R3 at the equilibria does not exceed three. This is because the energy function is a positive definite quadratic form in angular velocities; any n 6 i=l i=n+l 1 0 8 deviation from zero in angular velocities will result in a rise in the energy. R e m a r k s . • K,t>i,i>2 ,v 3) = (4,8,8,4) for all [TUT2) with T i ^ T i + T 2 ) ^ 0. • The local maxima and minima are in fact the maxima and minima. • The distribution of indices for equilibria (e,i)®_ 2 may be derived from that of the equilibrium en by a permutation of axes and a transformation. The procedure is detailed in Appendix 5-II. 5 . 3 D i s t u r b a n c e E n e r g y a n d A t t i t u d e P e r t u r b a t i o n By Morse's lemma [6], about a nondegenerate equilibrium defined by n = 0, e 6 {e ty}, E is of the form n 6 Eo + Y,Vi- J2 Vi> E0 = Ep{e) = E(0,e), 3 < n < 6 i=l t'=n+l for some local coordinate (y t), where 6 — n is the index of the equilibrium e. Let e — en,m — min2<i<Q\Ep(en) — Ep(e)\. By Palais-Smale theorem [7], manifolds 2? - 1 ( / i ) , E0 < h < E0 + m, are diffeomorphic and E~1((E0,E0 + m)) is diffeomorphic to i £ - 1 ( / i ) x (Eo,E0 + m). Consider a connected component B of E~1([Eo, EQ + m)) containing e. If e is a stable equilibrium, then there exists a diffeomorphism B {x e R6 : \x\ < m} (6) which maps E\g1(E0 + h) onto {x £ R6 : \x\ = h} for any 0 < h < m. In Figure 5-3, m* = m/(Ii + I2 + l3) is plotted versus the pair (T i , T2). Figure 5-4 depicts the 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 hatch-ing, as function of T\ and T%. 1 1 0 F i g u r e 5 -2 (b ) Index of the equilibrium e2i, equal to the multiplicity of hatch-ing, as function of 7 \ and T2. Ill F i g u r e 5-2(c) Index of the equilibrium e 3 i , equal to the multiplicity of hatch-ing, as function of Ti and T%. 1 1 2 F i g u r e 5 -2 (d ) Index of the equilibrium e+i, equal to the multiplicity of hatch-ing, as function of T\ and T2. 1 1 3 F i g u r e 5-2(e) Index of the equilibrium esi, equal to the multiplicity of hatch-ing, as function of T\ and T2. 1 1 4 Figure 5-2 (f) Index of the equilibrium e&i, equal to the multiplicity of hatch-ing, as function of T\ and T2. 1 1 5 configurations which do not "tumble" about en 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 s ) constant energy levels shrinking to the equilibrium as h > 0 approaches zero (equation 6). A lower dimensional analogy, is the family of concentric topological circles (S1) 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 n _ 1 x (0,1] such that for all t G (0,1] the image of S n _ 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 , jT2-space. What is shown, is the cross-section of equipotential surfaces of Ep/(Ii + 7 2 + I3) with the planes o*i,o2; 0 2 , 0 3 ; 0 1 , 0 3 as seen from the isometric vantage point with respect to the octant of RP(3) defined by 01 > 0,02 > 0, o 3 > 0. Compare with Figure 5-7 for the location of this octant in RP(3). A l l equilibria appearing in the octant are marked in Figure 5-5a save e6i which coincides with en from the assumed vantage point. Recall that any point Q = (°i) 92? 93) m this representation determines a unique element (at-y) G 50(3), obtained by a rotation of 2 arcsin \q\ about q. For example, a rotation of ir/2 about 62 axis, e 3 1 , has coordinates (0, V2/2,0) in this representation. 1 1 6 F i g u r e 5-3 Maximum allowable disturbance energy (nondimensional) before tumbling at equilibrium en as a function of T\ and T2. 1 1 7 F i g u r e 5-4 Optimal configuration contours obtained by equialtitude contours of the surface in 5-3. 1 1 8 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. 1 1 9 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. 1 2 0 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. 1 2 1 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. 1 2 2 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. 1 2 3 Figure 5-G(e) Equipotential contours in the configuration space 50(3) for inertia f = (—0.3,-0.3); equilibrium potential=0.45. 1 2 4 F i g u r e 5 - 6 ( f ) Equipotential contours in the configuration space 50(3) for inertia T = (—0.4,0.3); equilibrium potential=0.28. 1 2 5 F i g u r e 5-7 The location of the octant oi > 0,o2 > 0,o3 > 0 in RP(3). 1 2 6 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 en 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 en-ergy function. For example, for a satellite with inertia parameters T\ = —0.3 and T2 — 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 T2 = 0.3 (Figure 6b) the energy falls 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 o2, 0*3, and deviations in the angular velocity components ( f i i , II2, $^3). When the energy function has index two at the equilibrium (e.g., T\ = 0.4 and T2 = —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 T2 = —0.4, Figure 6d) the energy falls in all 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 T2 = —0.3, Figure 2e) two directions from the equilibrium lead 1 2 7 to lower energy levels (rotations about o 2 and 03); note, however, the difference between this situation and that of the other index two configuration (Figure 2c). A n analogous situation exists for the equipotential contours corresponding to con-figurations 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 T2 = 0.3 Figure 2b indicates that the energy function has a minimum at equilibria e 2 i , e 2 2 , ^23, « 2 4 (second row in Table 1) and hence the potential energy can only fall as one moves away from en close to e22 (Figure 6b); a result which confirms the predictions made based on the index of E given in the previous paragraph. 1 2 8 R e f e r e n c e s [l] Marandi, S.R. and Modi , V . J . , " A preferred coordinate system and the asso-ciated 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 Fre-quencies 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-Hil l , 1983, pp. 199-210. [4] La 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 F Et Direct calculation with differential forms [5] yields / 0 -f3 h d(aij) = (dij) f3 0 -fi \ - h fi 0 where (/,) are the left invariant differential 1-forms on 50(3) dual at E G 50(3) . Now, 3 dEp — 5^ ii{3aiidau — a3ida3i) i = i = {h - ^3)(3ai 2 ai3 - 0 3 2 0 3 3 ) / l + (I3 - I i ) ( 3 a u a i 3 - 031*133)h + {h ~ ^ 2 ) ( 3 a n a i 2 - a 3ia 3 2 ) / 3 . Therefore, dEp = 0 and TXT2T3 ^ 0 imply 3 d i 2 a i 3 = tX32a33, 3 a n a i 3 = 03 i033> 3 a n a i 2 = 031032. 1 3 0 A P P E N D I X 5 - L T : T H E E X P R E S S I O N O F Ep A T E Q U I L I B R I A Expansion of Ep at critical points e,i up to the second order in local coordinates a = (o~i,o~2,o~z) yields the following. A t en, i ( 3 / ! - 73) + \tTA.t, t = a; (16(-7 2 + 73) 0 0 where A = 0 64(-7i + 73) 0 V 0 0 48(-/1 + /2) A t e 2 1 , |(3/i - h) + ^ (1 + V2)2tTAt, t = a-{V2- 1,0,0); /2(/ 2-I 3) 0 0 where A = 0 - 7 7 2 + 4 7 2 + 3 7 3 -Ix + 4 7 2 - 3 7 3 V 0 - 7 x + 4 7 2 - 3 7 3 - 7/X + 472 + 3/3 A t e 3 1 , ^{-h + 3 / 3 ) + ^ (1 + V2)HTAt, t = 3 - (0,1 - y/2,0); / h + 272 - 3 7 3 0 h ~ 4 / 2 + 3 / 3 where A = 0 8(7! - 73) 0 \ h - 4 7 2 + 3 7 3 0 7X + 272 - 3 7 3 At e 4i, - ( 3 7 2 - 73) + -t TAt, r = a + ( — , — , 0 ) ; / 5(7!- 272 + 73) - 7 7 i + 272 + 5 7 3 (-7X + 4 7 2 - 373)\/2 \ where A = - 7 7 i + 272 + 5 7 3 5(7i - 272 + 73) {-h + 4 7 2 - 3I3)V2 V ( - 7 i + 4 7 2 - 373)\/2 (-7i+ 4 7 2 - 373)\/2 2{-Ix - 4 7 2 + 57 3 ) J A t esi, + 3 7 2 ) + \tTAt, t = a + (1,1,1); /4 ( 7 x-7 2) 0 0 where A = 0 3(-7 2 + 73) 0 \ 0 0 7 i - 7 3 A t e6i, U-h + 3 7 3 ) + \tTAt, t = a - (1,1,1); 1 3 1 ' 3 ( J i - 7 3 ) 0 0 where A = | 0 - Ix + 7 2 0 0 0 4(72 - 7 3 ) Index of Ep at e n is zero if -7 2 + 7 3 > 0 & - h + 7 3 > 0 & - 7i + 7 2 > 0; (1) equivalently, Index—1 if Index=2 if Index=3 if Ti < 0 & T 2 > 0 & T 3 < 0. (Ti > o & r 2 > o & r 3 < o) or (Ti < 0 & T 2 < 0 &: T 3 < 0) or (Ti < 0 & T 2 > 0 & Ts > 0). (Ti > 0 & T 2 < 0 & T 3 < 0) or ( i i > o & r 2 > o & r 3 > 0) or (Ti < 0 & r 2 < o & r 3 > 0). T i > o & r 2 < o & r 3 > o. 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 Ev 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 2 and I\ for 7X and 7 2 in (1) to find - 7 2 + 7 3 > 0 & - h + 7 3 > 0 & - 7 2 + h > 0 1 3 2 in order for Ep to have index zero at e^\- Equivalently Ti < 0 & T 2 > 0 & T 3 > 0. This leads to a partitioning of the parameter space for e4i which is the mirror image of that of en with respect to diagonal Ti + T 2 = 0 (Figures 5-2a and 5-2d). 6 . A T T I T U D E S T A B I L I T Y O F R I G I D S A T E L L L I T E S V I A A N O R M A L I Z E D H A M I L T O N I A N In studying the behaviour of solutions to Hamilton's equa-tion 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. By a sequence of canonical trans-formations 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 I n t r o d u c t i o n 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 1 3 3 1 3 4 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-1 3 5 spherical shape of the moon" he derived the four inequalities ( I 3 - h ~ h)2 + (/3 - h)h + 4 ( / 3 - h)h > 0, ( / 3 - / 2 ) ( / 3 - / i ) >0 , [(/3 - I2 ~ h)2 + (/s - h)h + 4 ( / 3 - h)h\2 > 16/2/!(J3 - / 2 ) ( / 3 - / i ) , (/2 - h) > 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 1 3 6 of stability or instability remained unanswered [10, 11]. A parallel problem which has stimulated much growth and controversy in dy-namics 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 prob-lem. 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, King 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-1 3 7 S A T E L L I T E Figure 6-1 Two classes of satellites for which the equilibrium orientation is stable in the linearized analysis. 1 3 8 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 struc-ture of the space of polynomials. The latter appeared in the work of Deprit and Henrard [22] and is summarized by Churchill, Kummer, and Rod [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 H a m i l t o n i a n a n d t h e C h o i c e o f G e n e r a l i z e d C o o r d i n a t e s The Rodrigues parameters r = ( r i , r 2 , r 3 ) provide a local coordinate system about E G 50(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 . / 1 + r\ - r\ - r f 2{rir2 - r 3 ) 2 ( n r 3 + r 2 ) x ' - * - „2 I „2 Jl (a*j) = h M\P\2\ 2 ( r i r 2 + r s) l-rl + rl-4 2 ( r 2 r 3 - n ) U + M j y 2[r1r3-r2) 2 ( r 2 r 3 + n) 1 - r\ - r\ + r\ where (a ty) denotes the orientation matrix of the body frame with respect to the 1 3 9 orbital frame. Using (a t y) r (d t J ) one finds the time rate of Rodrigues parameters in terms of U (analogous to the method shown in Appendix 4-II): r l r 2 — r3 r l r 3 + r2 r 2 r 3 - r i 10. 1 + r i T2 I = - I n r 2 + r h ) V r* r3 _ f2 r l + r2T3 l + r l Next, the conjugate variables s = dL/dr are computed, where L is thought of as L(r , r ) after the appropriate substitutions for U and (a,y): 1 + |r ^1 + ^31 n 2 + ^3 + a33 Therefore, n = 1 + r ?|2 r i / 2 (1) 6 . 3 T h e H a m i l t o n i a n u p t o t h e F o u r t h D e g r e e Equation (1) implies that the equilibrium ((aty),n) = (E,0) G SO(3) x so(3) corresponds to (0,0,0,0,0,2/ 3 ) in (r, s ) coordinates. To bring the equilibrium to 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) 1 4 0 yields the expression H = H0 + H1 + H2 + H3 + H4 + ---, where #1 = 0, Here, I2 h 0 o 0 0 1 + 2 / 2 0 \ 0 -12(/x lz h 0 0 0 0 0 1 2 ( -h + /2) 0 0 0 0 1 *3 0 1 0 0 _ 1 + i 2 _ 1 ^ 2 J 2 0 0 0 1 4 / 2 0 V 0 0 0 0 0 ^3 J The homogeneous polynomials of degree three H3 and four H4 are given in Appendix 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 H2, H3,... 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 H2 and H3 change depending on the chosen coordinate system. That is why it is important to recognize the variable with respect to which H is expanded. 1 4 1 6 . 4 N o r m a l i z a t i o n P r o c e d u r e 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 N o r m a l i z a t i o n u p t o t h e s e c o n d d e g r e e Second degree terms H2 can be brought to the form W l ( _ L _ ) + W 2 ( _ ) + W g ( _ ) by a linear canonical change of variables (r, s) (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 - ° E o ) ( l> Hence, the linear part of this system comes from H2. As a result, is the linearization of the equation (2). If the eigenvalues of L are distinct and pure imaginary, then by a linear canon-ical transformation /o, equation (3) is transformed to 1 4 2 where / 0 x L / o is of the form (0 0 0 - W l 0 0 \ 0 0 0 0 —(JJ2 0 0 0 0 0 0 - w 3 0 0 0 0 0 0 W 2 0 0 0 0 V o 0 U>3 0 0 0 J 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: h r 3 V 0 •1 + 4 0 0 0 0 0 12(Ji - J 3 ) -0 0 7*7 0 _1 0 0 13_ •1 + 0 0 J3_ 2/i 1 4/2 _ J3_ 2/ 2 0 0 0 0 0 0 0 0 12(A - / a ) 0 A 0 0 0 1 o J (ri\ S\ S2 r 3 \s3J Consider a subset of the above equations, fri\ T2 = L T2 •Si S\ \s2J \s2J 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,T2) (Figure 6-3). To avoid the confusion of adjacent spectra they are marked in alternating colors. The eigenvalues 1 4 3 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) = X 4 + B X 2 + C = 0 has four pure imaginary roots —1X2,-1X1,1X1,1X2 provided B = -T1T2 + 3T 2 + l > 0, C = -AT1T2 < 0, B 2 - A C = r 2 T 2 2 - 6Tir22 + 1 4 T i T 2 + 9T 2 2 + 6T 2 + 1 > 0. 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 < Ai < A 2 . Let vi be the eigenvector corresponding to the eigenvalue t'Ai. Define e i N ? = I m v i • I 0 E ) 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 x = - A x e i ^ " ^ 1 = - X i € i b 4 , N i N i £ i N i L64 — — — l i l r n v i = — — R e v 1 = Aie i&i . e i N i e i N i 1 4 4 (a) x—*• -x—x ( b ) x—x-x- -x-x X F i g u r e 6-2 Some possible spectra for linear Hamiltonian systems of (a) two and (b) three degrees of freedom. Figure 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. 1 4 7 Similarly, define 6 2 = -rz-Revi, i v 2 1 -Imv-i Then, in the basis {61,62,64,65}, ^ reads / 0 0 e iAi 0 \ 0 0 0 c 2 A 2 - « i A i 0 0 0 V 0 — 62^2 0 0 J Define 0 0 1 2V/I3A3" 0 0 0 0 0 0 0 where A3 = y/STs. In the basis {61,62,63,64,65,66}, ^ * s * n c a n o n i c a l form / 0 0 0 0 0 \ 0 0 0 0 0 0 0 0 0 0 U>3 0 0 0 0 0 0 —0J2 0 0 0 0 V 0 0 - w 3 0 0 0 J where u>\ = e i A i , u2 — e 2 A 2 , and W3 = A 3 . Wi th a canonical change of variables (61 6 2 63 64 65 6 6 ) ( .j), the expression H2 becomes i = i A + yf in x, y variables. The characteristic frequencies w t 's are functions of Tx and T2, defined on the hatched region of Ti,T 2-parameter space. Variations of u>\, u2, w 3 with respect to the parameters T\ and T2 are presented in Figures 6-5 to 6-7. Contour plots cover those regions of T\, T2-space for which the corresponding wt- is denned. Note the change in the sign of oj\ during transition from T2 > 0 to T2 < 0 . 1 4 8 F i g u r e 6-5 Variation of u>i as a function of T\ and T2. -1.0 T i • f,, as a function of T V and T 2 . F i g u r e 6-6 Variation of « a as a t F i g u r e 6-7 variation of W3 as a function of T\ and T2. 1 5 1 6 . 4 . 2 E l i m i n a t i o n o f t h e t h i r d d e g r e e t e r m s In this subsection, a generating function K is found which yields a canoni-cal transformation leading to the elimination of terms of degree three in the new Hamiltonian. Def in i t ion . The characteristic frequencies wi , u2, u>3 satisfy a resonance relation of order M if there exist integers m,- not all equal to zero such that miwj + m2oj2 + m 3 w 3 = 0, | m i | + | m 2 | + | m 3 | = M. Def in i t ion . A Birkhoff normal form of degree s for a Hamiltonian is a polynomial of degree 5 in the canonical coordinates (u t-,v t) which is actually a polynomial (of degree [s/2]) in the variables r» = (u 2 + v 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: z _ X = —h tVJ, 2 y - i- + w. it With this H2 becomes ihJ\Z\VJ\ + iu)2Z2VJ2 + ibJ3ZzWz. Note that ZJVJJ = —i(x2- + J/y)/2. By Jacobi's theorem [26], the implicit relations _ dK w dK w = u = az ov define a canonical transformation (z, w) (u,v) about (ZQ,WO), provided the gen-1 5 2 erating function K(v, z) satisfies the condition d2K det—^—^y0,z0) ^ 0, ovoz with vjQ = dK/dz(v0,zo). Let K = v • z -\- K3{y,z), where K3 is a homogeneous polynomial of degree three. Then d2K det——(0,0) = det~E = 1, ovoz W = 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: - dK3 z = u , Proof: a,Q = 0 in the expansion of 2 = a*o + A i v + A2wH 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 in the expression z — u — dK3/dv(v, z) as dK3/dv 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 dK3/dv(z,v) as higher degree terms will lead to third or higher degree expansions of z. Hence, z=u-^-(u,v) + 0(3). ov 1 5 3 Similarly, rf=v+^{u,v) + 0{3). ou 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 + ^ ) + H * ( * > * ) + O f t 3 = 1 3 3 v — 3 ^ dK dK 3 = 1 j=l 3 3 Any term in K may be written as s^z m v n, with |m | + \n | = 3 which is an abbreviation for S m S ^ r 1 ^ 2 ^ 3 ^ 1 " ™ 2 ^ 3 ) with mi+m2 + rn3 + ni+n2+n3 — 3. Here, SfhH is the coefficient of the z mv n monomial. Hence, the coefficients of u m v n in the third degree terms of H reads as Smn[iui{ni ~ mi) + iu2(n2 - m2) + »'w3(ra3 - m 3)] + fc^g, (5) where hfha l s the coefficient of the monomial u mv n in H3(u,v). 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 f o r m t h r o u g h t h e f o u r t h d e g r e e t e r m s In this subsection, a generating function is found which normalizes H through terms of degree four. The procedure is analogous to the one presented in the previous subsection. 1 5 4 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 K4(u, f ) be a homogeneous polynomial of degree four. The implicit relations _ 8K4 ? ^ dK ou 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 { sms[»'wi(ni - mi) + iu2(n2 - m2) J = l |m | + |n 1=4 + »w 3 (n 3 - m 3)] + + 0(5). 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,-, i = 1,2,3 the monomials / i m n f m C * n remain. Therefore K4 may be chosen ( the choice is not unique) such that 3 B(t,$) = Y^iujtjtj + J2 hrnnrr1+0(5). j = l m=n 6 . 4 . 4 N o r m a l i z a t i o n b e y o n d t h e s e c o n d d e g r e e t e r m s — L i e s e r i e s In this section, a procedure is outlined which yields the normalized form with-out recourse to a generating function. First the general procedure is described. 1 5 5 Let Pr denote the set of homogeneous polynomials of degree r. The direct sum of vector spaces Pr, P = ©£2_ 2P r, is a positively graded real (or complex) Lie algebra [28] with \F,G] = Y— — -—— d X J d y l dy3 dX3 With Zj — Xj + iyj and Zj = Xj — iyj, Subscript elements Hr E Pr- Let K E Pa. Define adK :P -> P, [H, K] as a graded homomorphism of degree s — 2. Denote adjK{Hr) = adK{adj-\(H4)) 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. As a consequence, we can define a linear operator exp(adx) • P —* P by exP(adK)(H) = f : ^ l j=o J Def in i t ion . A n element H = ©£2_2.firr € P is in the normal form through terms of order m > 2 if adn2 [Hr) — 0 for 2 < r < m. Here, it is assumed that This definition agrees with the one used before. P r o p o s i t i o n . Let H = © £ L 2 i I r E P be in the normal form through terms of order (m - 1) > 2. Let Hm = Fm + Gm, where Fm E Ker(adH2\pm) a n d Gm G 1 5 6 ran(adH2\Pm)- Then exp(adKm){H) E P is in the normal form through terms of order m, agrees with H through terms of order m — 1, and has Fm as m th term. Proof: exp(adKm)(H) - H2 + H3-\ h H m + adKm [H2) + (terms in Py with j > m + 1}. Therefore, Hm + adKm(H2) = Fm + G m - adH2{Km) = Fm. A routine computation shows that adH2{z*Z*) = -i(m - ti).uzTZn. As a result, if Gm = E 6 - " ^ ^ e r a n i a d H 2 \ p m ) , then Km = iJ2Km - n)^\-1bmHz'hZii. The above procedure is applied to the Hamiltonian H = H2 + H3 + i f 4 + • • • • Let Hz ^Y.^fhnZfhZK. Then, K3 = » E [ ( n » - n j . w ] - 1 ^ . Therefore exp{K3){H) = H2 + H3 + H4 + adKs{H2) + adKs{H3) + \ad\z{H2) + 0(5). On the other hand adK3(H2) = — adff2(K3) — —H3. Therefore, the fourth degree terms in exp(adK3){H) become Jf4 = HA + adKs(H3) - ^adKsH3 = H4 + ^adK 3H3 — H< + ^\H3,K3\. 1 5 7 Once H4 is computed and expressed in the complex form cr%azmZn the fourth degree terms in the normal form expressed in complex variables are of the form Y^rn=nc™nZrnZn. In summary, to compute the fourth degree terms in the normal form the following three steps are taken. i . The expression [ # 3 , ^ 3 ] is computed in complex variables. i i . H4 = H4+ ^[H3,K3] is found in complex variables YlerhH^mZn. i i i . The coefficients Crnn with rh = n are selected. These yield w,y's as follows "11 = 4c20020Cb = W12 = 2cuono, = w 1 3 = 2 c 1 0 i i o i j W 2 2 = 4c 0200205 W32 - W 2 3 = 2conoii , ^33 = 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 gen-erating 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]. 1 5 8 6 . 5 S t a b i l i t y o f a n E q u i l i b r i u m — t h e D e l p C a s e 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 wt-y's are determined as functions of T\ and T 2 . The following is the stability condition found by Thkai [5]: = (0,0) is stable if 3 3 Vr t -> 0 fir^O u; tr t-= 0 and w.yr.ry = 0. (6) »=i 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 3 r t > 0 , i= 1,2,3; ^ w t r t = 0; J^ 7 * = 1 i=i 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 2 + r 3 = 1 and P 2 : Yli=i w i r t = 0, form the function (the method of Lagrange multiplier) 3 3 V» = <f> + A i £ n - l ) + A 2 J2 1 5 9 Denote the unique solution to <90 3(f, A i , A 2 ) / W l l ->12 W13 1 ( T l \ W21 W22 <~23 1 U>2 T2 0 W31 W32 w 3 3 1 + T3 + 0 1 1 1 0 0 A l - 1 V OJl W2 OJ3 0 0 J u 2 / v 0 y by ( f b , A i 0 , A 2 o ) . The line P intersects the three planes r,- = 0 , t = 1,2,3 (Figure 6-8) at points 0 "3 w 3 - u > 2 - ^ 2 Do = and Da = w 3 - u ; i CJ 2-WI - - 1 w 2 - W ! 0 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/w2> W 1 / W 3 , W 2 / W 3 are negative when wi , W2» and w 3 are not all of the same sign. The reciprocal of the nonzero components of Dt, i = 1,2,3 are respectively D3 1 _ f l 0J3 f l 0J3 f l U>2 1 _ 1 > U>2 1 - - * 1 _ 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 2 ) for which w t 's of the same sign have equal magnitudes. Assume, with out loss of generality, that D 2 and D 3 are in the first octant. Let 4>i — 4>(Di), i — 2,3. In the following chart, if the sequence of verifications does T3 D2N D1 F i g u r e 6-8 The geometry of Thkai's theorem. 1 6 1 not lead to "otherwise" the relation (6) is satisfied: ( mm(c/>o,02,^3) > 0 or max(<f>0, <£2, <j>z) < 0, \ otherwise, ( min(<f)2,<f>3) > 0 or max(<l>2,4>z) < 0» I otherwise. 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 <£(Ar) = \2<J>(T) = 0 and Ari + Ar 2 + Ar3 = 1 for some A G R. The maximum M and minimum m value of <f> on the segment D\D2 is detemined in each of the two cases corresponding to To ^ D1D2 and fo G P\DiD2- The condition (6) is satisfied when 0 G" [m,M] . The curves in the Delp region corresponding to which the characteristic fre-quencies satisfy a resonance relation are shown in Figure 6-9. Outside these reso-nance 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). (D2 - T0) • ( £ 3 - T0) 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 . 1 6 4 R e f e r e n c e s [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. 484-85. [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 Pub-lishing 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 Publica-tions, 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 fre-quencies in a circular orbit," The Journal of the Astronautical Sciences, Vol . 8, 1961, pp. 14-17. [10] Likins, P .W., "The influence on dynamics and control theory of the space-craft attitude control problem," Proceedings of the Sixth Canadian Congress of Applied Mechanics, Vancouver, B . C . , Canada, Editor: V . J . Modi , 1977, pp. 1 6 5 321-35. [11] Kane, T.R. , Likins, P .W., and Levinson, D .A . , Spacecraft Dynamics, McGraw-Hi l l , New York, 1983, pp. 199-210. [12] Poincare, H . , Les Methodes Nouvelles de la Mecanique Celeste, Vol . 1, 2, 3, Gauthier-Villars, Paris, 1892, 1893, 1899. [13] Birkhoff, G.D. , Dynamical Systems, American Mathematical society, Provi-dence, 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 Rod, 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, 2nd 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, Read-ing, Massachusettes, 1974, p. 457. 167 A P P E N D I X 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 H3 = A T T — r 3 s i s 2 + A T T r 2 s i s 3 H — — r i ~ 3 s i 4 i i i2 4Iil3 2iii 3 -h + h ^ h ( h - h ) H — 7 7 - 7 — r ! S 2 s 3 H — — r 2 r 3 s 2 4 i 2 i 3 , h 2 , ^ 3 2 , - 2 + 2 i ; r i 5 3 + 2 7 r r 2 S 3 + 2 r 3 5 3 -12/x/f + 127x72/3 + /i/ 3 2 - / 2/ 3 2 + — r x r 2 r 3 ixJ2 HA h h + h h - h h , / s ( / i - / 2 ) 2 = 7TT7 r 2 ~ 3 s 2 3 3 + r x r | 5 2 4 / 1 / 2 / 3 4/x/2 — / l + 73 2 , - 2 2 , - 2 2 + - 2 7 - — - i - s r a + ^ r l 5 3 + — r 2 , 3 1 2 2 , / 3 ( / l ~ h) + T 7 - r 3 s 3 + T 7 r l r 2 T 3 S 3 4 i 3 Jx-12 J. f or 1 f U , , / o j c r U . 2 1 - 2 4 / i / 2 + 16/x/3 +/| 2 2 + (-oh + i 3 J r x + 2(dl2 - 513)rlr2 H — r x r 3 + for 11 7 V 4 + 36 / i / 2 - 24/ | - 20 / 2 / 3 + / 3 2 0 2 , 3(6/x - 8/ 2 + / 3 ) + (9ix - H i 3 j r 2 + — r 2 r 3 H - r E p i l o g u e A science produced with a view single to its applications is impossible; truths are fruitful only if they are concate-nated; if we cleave to those only of which we expect an immediate result, the connecting links will be lacking, and there will be no longer a chain. Henri Poincare S u g g e s t e d F u r t h e r 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 dm/R (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 — / XiXjXkdm for repeated indices. Recall 1 6 8 1 6 9 that for the inertia tensor hi = / Xk ~ xixj)dm-J k=l 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 corre-sponding to which the evaluation of the right hand side of differential equations required a shorter C P U time than those of transcendental rotational coordi-nates 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 2 , however, is not only a function of f(x) but also dependent on the time step of the integration. The time step in a numerical integration 1 7 0 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 in 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 classifi-cation 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 strat-ifies 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. 1 7 1 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. C o n c l u d i n g R e m a r k s 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 3 and a homomorphism mapping the universal cover of 50(3) onto SO(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 frac-tions 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 1 7 2 of four based on the dynamical behaviour. (b) The existence of a stable equilibrium orientation is a result of the com-pactness of 50(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 consid-ering equipotential surfaces in the configuration space. (d) The long outstanding problem of stability of the equilibrium en is resolved by normalization of the Hamiltonian about en through terms of degree four and application of Tkhai's theorem. 1 7 3 R e f e r e n c e s [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 clas-sical 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. 588-97. [6] Goldman, L . , "Integrals of multinomial systems of ordinary differential equa-tions," 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 |
AggregatedSourceRepository | 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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0098148/manifest