University of Southern California Dissertations and Theses

Periodic motion analysis in spacecraft attitude dynamics

(USC Thesis Other)

Periodic motion analysis in spacecraft attitude dynamics

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Periodic Motion Analysis in Spacecraft Attitude Dynamics by Dayung Koh A thesis submitted to the Faculty of the Graduate School of the University of Southern California in partial fulfillment of the requirements of the degree of Doctor of Philosophy Department of Astronautical Engineering May 2016 Dean Prof. Henryk Flashner Committee Dr. Dan Erwin Dr. Mike Gruntman Dr. Rodney L. Anderson Dr. Robert Sacker Dedication I dedicate this dissertation to my parents 고희경 and 박 효경, to my sister 고다 혜, and to the rest of my family. Acknowledgments I am here with all the support of many people over the years. I really thank to everyone that I know. I would like to begin by acknowledging my immense gratitude to my advisor, Professor Henryk Flashner. He has always been there to guide me, and to encourage me throughout my graduate studies and life. I cannot think of my graduate years without him. I also wish to thank Rodney Anderson. He was a good listener whenever I have questions, and he was the best supporter to pursue my dream. His time and advice have been invaluable. I have been fortunate to have many great people help me along the way. I would like to thank Ryan Park for his valuable advice and his willingness to support me. He has taught me many things especially about the world out of the classroom. I had a great experience to work with David Barnhart for two projects. I would like to thank him for his trust, all his supports, and his friendship. I also would like to thank Bongsub Song who inspired me as I was studying in Dynamics and Control, Mechanical engineering. He helped me to find my interest. The support of the Astronautical engineering at USC for my study is gratefully acknowledged. I cannot thank Dr. Erwin and Dr. Gruntman enough for their trust and guidance. A number of people served on my committee and supported me in and out. Thank you Dr. Kunc, Dr. Sacker, Dr. Yang, Dell Cuason, Norma Orduna, Linda iv Ly, Ana Olivares, and Marrietta Penoliar. All of my friends here in California and back in Korea, thank you all for your friendship. I thank my current roommate, Kunhwa Kim for her kind concern and friendship. I would especially like to thank Jules Lee for sharing her happiness and sorrow. You are the best. For the last, I would like to thank my lovely family Heekyung Koh, Hyokyung Park, and Dahye Koh for their endless support of me through this process. I will see you soon. v Contents Acknowledgments iv 1 Introduction 1 1.1 Gravity Gradient Satellite . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Spinning Satellite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2 Model of Satellite in Orbit 8 2.1 Orbit definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2 Reference Frames . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 Transformation between Coordinate Frames . . . . . . . . . . . . . . . 12 2.3.1 Euler angle sequence . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3.2 Transformation from Inertial to Orbital Inertial Frame . . . . . 13 2.3.3 Transformation from Orbital Inertial to Local Orbital Frame . . 14 2.3.4 Transformation from Local Orbital to Reference Frame . . . . . 14 2.3.5 Transformation from Reference to Body Frame . . . . . . . . . . 15 2.4 Gravity Gradient Torque . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.5 Euler’s Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.5.1 Global Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 vi Contents 3 Analysis Methodology 23 3.1 Point Mapping Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.1.1 Periodic Solutions and Local Stability Properties . . . . . . . . 24 3.1.2 Bifurcation Conditions . . . . . . . . . . . . . . . . . . . . . . . 27 3.1.3 Truncated Point Mapping Method . . . . . . . . . . . . . . . . . 28 3.2 Cell Mapping Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.2.1 Periodic Motions and Domains of Attraction . . . . . . . . . . . 32 3.2.2 Global Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.3 Extended Cell Mapping Method . . . . . . . . . . . . . . . . . . . . . . 35 3.3.1 Computation of P-1 Solutions . . . . . . . . . . . . . . . . . . . 35 3.3.2 Computation of P-K Solutions . . . . . . . . . . . . . . . . . . . 37 4 Analysis of Pitch Motion of a Gravity-Gradient Satellite Stability 42 4.1 Equation of Motion for Pitch Motion . . . . . . . . . . . . . . . . . . . 42 4.2 Periodic Solutions of Pitch Oscillations . . . . . . . . . . . . . . . . . . 44 4.2.1 Periodic Solutions using Cell Mapping . . . . . . . . . . . . . . 44 4.3 Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.4 Bifurcations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.4.1 P-1 to P-1 bifurcation . . . . . . . . . . . . . . . . . . . . . . . 59 4.4.2 Bifurcation P-1 solutions to P-2 solutions . . . . . . . . . . . . . 63 4.4.3 P-1 to P-3 Bifurcation and P-3 to P-9 Bifurcation . . . . . . . . 66 4.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5 Analysis of Spinning Satellite 75 5.1 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.2 Equation of motion for spinning satellite . . . . . . . . . . . . . . . . . 77 vii Contents 5.3 Periodic Motion Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.3.1 Periodic solutions using Cell Mapping . . . . . . . . . . . . . . . 78 5.3.2 Periodic solutions verification procedure . . . . . . . . . . . . . 78 5.4 Effect of spin direction e = 0.1, k 1 = 2 . . . . . . . . . . . . . . . . . . 79 5.4.1 Periodic solutions for positive spin direction (σ≈ 0.5) . . . . . . 79 5.4.2 Bifurcations for 0.3<σ< 0.5 . . . . . . . . . . . . . . . . . . . 84 5.4.3 Periodic solutions for negative spin direction (σ≈−0.5) . . . . 85 5.4.4 Bifurcations for−0.6<σ<−0.5 for small motion . . . . . . . . 90 5.4.5 P - 1 to P - 1 Bifurcation for−0.6<σ<−0.5 for small motion 91 5.4.6 P− 1 to P− 2 Bifurcation for−0.6<σ<−0.5 for small motion 95 5.4.7 P− 1 to P− 3 Bifurcation for−0.6<σ<−0.5 for small motion 99 5.5 Effect of high spin rate . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.6 Effect of orbital eccentricity . . . . . . . . . . . . . . . . . . . . . . . . 108 5.7 Effect of satellite shape . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 5.7.1 k 1 < 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 5.7.2 k 1 > 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 5.7.3 Special case for k 1 = 1.25 . . . . . . . . . . . . . . . . . . . . . . 117 5.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 6 Conclusion and Future Research 121 6.1 Summary and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . 121 6.2 Future Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 6.2.1 Coupled Orbit-Attitude Dynamics Model . . . . . . . . . . . . . 122 6.2.2 Asteroid Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 Bibliography 124 viii Abstract ix Dayung Koh (Ph.D., Astronautical Engineering) Periodic Motion Analysis in Spacecraft Attitude Dynamics Thesis directed by Professor Henryk Flashner (Chair) In this dissertation, the complete dynamics of a spacecraft in a general orbit is studied. For this study, a new computational approach which combines analytical and numerical analysis is introduced to investigate the global behavior of the full nonlinear dynamical models of satellites. This approach consists of two methods: point mapping and cell mapping. The point mapping algorithm is employed to study local stability characteristics and bifurcation conditions, and the cell mapping method is utilized to analyze periodic solutions and their global regions of attraction. The proposed approach does not depend on the assumption of a small parameter and, therefore, no limitations are imposed on the magnitudes of parameters or amplitude of motion. This approach is applicable in a wide range of applications. Two classes of satellites are investigated for their global attitude dynamics: gravity gradient satellites and spin- ning satellites in the elliptic orbits. First, a gravity gradient satellite as an example of a one degree of freedom problem is analyzed. Then it is extended to a two degree of freedom problem in which coupling between spinning motion and pitch motion occurs. As a result, families of periodic solutions and rich dynamic phenomena, which have not been found before, are investigated. Global regions of attraction and evolution of peri- odic solutions are demonstrated using invariant surfaces, and stability characteristics of the solutions are studied. A comprehensive investigation of bifurcation diagrams was also performed for studying transitions to various types of bifurcations and chaos with varying parameters. x List of Figures 2.1 Elliptic Keplerian orbit characters: a,b, and p . . . . . . . . . . . . . . 8 2.2 Orbitalelementstodefinepositionoftheorbitinspacew.r.t. theinertial frame: i, Ω,ω,and ν . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3 Relation between the orbital inertial frame and the local orbital frame . 11 2.4 Relation between the local orbital frame and the reference frame . . . . 12 2.5 Spacecraft B in a gravitational field . . . . . . . . . . . . . . . . . . . . 17 3.1 Schematicphaseplaneforanequilibriumpointx ∗ andperiodicsolutions X ∗(i) ,i = 1,··· ,K for K = 3 . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Construction of a cell state space for N = 2 . . . . . . . . . . . . . . . 31 3.3 Schematic cell state space showing an equilibrium cell z ∗ and P−K cells z ∗(i) ,i = 1,··· ,K for K = 3 . . . . . . . . . . . . . . . . . . . . . 33 3.4 P− 1 solutionz ∗ of cell j with small perturbation z . . . . . . . . . . . 35 3.5 Neighborhood of cell j might have a P− 1 Solutionz ∗ with small per- turbation z . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.6 Cell j i are P−K cells with perturbation z i where i = 1,··· ,K . . . . 38 4.1 Geometry of gravity gradient satellite motion . . . . . . . . . . . . . . 43 4.2 Cell state space with e = 0.01, k 2 = 0.1 . . . . . . . . . . . . . . . . . 45 xi List of Figures 4.3 Global behavior in a cell state space for e = 0.01,k 2 = 0.1 . . . . . . . 46 4.4 Cell state space shownP−1 solutions fore = 0.01,k 2 = 0.1. A group of P− 1 cells from the cell mapping converge toA = (−1.95e− 7, 0.0283370) 48 4.5 Cell state space shown stable P− 2 solutions for e = 0.01, k 2 = 0.1. A group of P− 2 cells from cell mapping converge to B i = (−0.6131613, 0.0239830), (0.6131484, 0.0239292) . . . . . . . . . . . . . 50 4.6 Cell state space shown unstable P− 2 solutions for e = 0.01, k 2 = 0.1. A group of P− 2 cells from cell mapping converge to C i = (0.0001175, 0.3147778), (0.0008436,−0.2638848) . . . . . . . . . . . . . 51 4.7 Cell state space shown P− 3 solutions and trajectories in phase plane for e = 0.01, k 2 = 0.1 (S i : stable solutions, U i : unstable solutions) . . 52 4.8 Cell state space shown stable P− 3 solutions for e = 0.01, k 2 = 0.1 S i : (−0.0006549, 0.5332902), (1.1431803,−0.1327799), (−1.1426959,−0.1331036) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.9 Cell state space shown unstable P− 3 solutions for e = 0.01, k 2 = 0.1 U i : (1.1323459, 0.1684754), (−1.1339993, 0.1675323), (0.0009308,−0.4855702) 54 4.10 Fast Fourier analysis results fore = 0.01,k 2 = 0.1 at three fundamental periodic solutions and at a neighborhood point (0.2, 0) . . . . . . . . . 55 4.11 P− 1 toP− 1 bifurcation for three different ‘e’s by changingk 2 . (Solid line: stable solution, Dashed line: unstable solution) . . . . . . . . . . . 60 4.12 Numerical P− 1 solutions for different e and k 2 (Does not show the stability of solutions) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.13 Invariant surface of phase plane representation fore = 0.01. One stable P−1 solution bifurcates into two stableP−1 solutions and an unstable P− 1 solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 xii List of Figures 4.14 Bifurcation of unstableP−1 to stableP−1 and unstableP−2 solutions for k 2 = 0.1,e = 0.01∼ 0.7 at bifurcation point bfc1 (for bfc2 : see Fig.4.15). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.15 Bifurcation of stable P− 1 loses stability and stable P− 2 solutions emerge for k 2 = 0.1,e = 0.01∼ 0.7 at bifurcation point bfc2 (for bfc1: see Fig.4.14) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.16 Bifurcations for k 2 = 0.1: combined Fig.4.14 and Fig.4.15 . . . . . . . 65 4.17 Invariant surface of phase plane representation for k 2 = 0.1 . . . . . . . 66 4.18 P− 1 to P− 3 bifurcation and P− 3 to P− 9 bifurcation for e = 0.1 in 3-D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.19 P− 1 to P− 3 bifurcation and P− 3 to P− 9 bifurcation for e = 0.1 from the view of x 1 -axis . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.20 P− 1 to P− 3 bifurcation and P− 3 to P− 9 bifurcation for e = 0.1 from the view of x 2 -axis . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.21 Continuous trajectory for stableP−3 solutionsg i ,i = 1, 2, 3 ate = 0.1, k 2 = 0.19 and P− 1 solution ’+’ . . . . . . . . . . . . . . . . . . . . . . 69 4.22 True anomaly vs. angular responses and angular rate responses for g 1 showing the period of three motion . . . . . . . . . . . . . . . . . . . . 69 4.23 Continuous trajectory for unstable P− 3 solutions h i , i = 1, 2, 3 at e = 0.1, k 2 = 0.19 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.24 Continuous trajectory for P− 9 solutions at e = 0.1, k 2 = 0.21 . . . . . 71 4.25 True anomaly vs. angular response for a P− 9 solution at e = 0.1, k 2 = 0.21 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.26 True anomaly vs. angular rate response for a P− 9 solution ate = 0.1, k 2 = 0.21 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 xiii List of Figures 4.27 Fast Fourier analysis results for a P− 9 solution at e = 0.1, k 2 = 0.21 . 72 4.28 Invariant surface of phase plane representation, e = 0.1, k 2 = 0.21 . . . 73 5.1 Spinning satellite description . . . . . . . . . . . . . . . . . . . . . . . . 75 5.2 β−β 0 (left), γ−γ 0 (right) cell mapping stroboscopic phase plane for P− 1 solutions (k 1 = 2, e = 0.1) . . . . . . . . . . . . . . . . . . . . . . 80 5.3 Eigenvalues of matrixH in the complex domain for twoP−1 solutions, σ = 0.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.4 Trajectory inβ−β 0 phase plane (left) andβ response vs. true anomaly (right) of P− 1 solution for 1·T, σ = 0.4 . . . . . . . . . . . . . . . . 82 5.5 Trajectory in β−β 0 plane (left) and β vs. true anomaly (right) of the P− 1 solutions for 14·T, σ = 0.4 . . . . . . . . . . . . . . . . . . . . 83 5.6 Fast Fourier analysis for P− 1 solution for σ = 0.4 . . . . . . . . . . . 83 5.7 Bifurcation diagram for the origin with positive spin direction 0.3 < σ< 0.5 (k 1 = 2, e = 0.1) . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.8 Eigenvalues ofH at the origin showing P− 1 to P− 1 bifurcation for 0.43<σ< 0.45 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.9 Periodic solutions in cell mapping cross-section phase plane for different σ, e = 0.1, k 1 = 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.10 Eigenvalues of matrixH for two P− 1 solutions, σ =−0.5 . . . . . . . 88 5.11 Trajectoryinβ−β 0 phaseplaneandβ vs. trueanomalyofP−1solution for 1·T . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.12 Trajectoryinβ−β 0 phaseplaneandβ vs. trueanomalyofP−1solution for 30·T . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.13 Bifurcation map around the origin for negative spin (e = 0.1, k 1 = 2) . 91 5.14 Bifurcation sequence P− 1 to P− 1 around σ =−0.58 . . . . . . . . . 92 xiv List of Figures 5.15 Eigenvalues ofH evaluated at the origin, σ from−0.6 to−0.53 . . . . . 93 5.16 Bifurcated P− 1 solution, σ =−0.58 . . . . . . . . . . . . . . . . . . . 94 5.17 Fast Fourier analysis for stable P− 1 solution, σ =−0.58 . . . . . . . . 94 5.18 P− 1 to P− 2 bifurcation diagram as function of σ where σ≈−0.5 . . 95 5.19 Eigenvalues ofH with varying spinning parameter−0.52<σ<−0.48 96 5.20 Trajectory inβ−β 0 phase plane (left) andβ response vs. true anomaly (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.21 Fast Fourier analysis on claimed P− 2 solutions for four different σs . . 99 5.22 P− 1 to P− 3 bifurcation sequence at σ≈−0.54 . . . . . . . . . . . . 100 5.23 Eigenvalues ofH for spinning parameter−0.55<σ<−0.52 . . . . . 100 5.24 P− 3 solution for σ =−0.54 . . . . . . . . . . . . . . . . . . . . . . . . 102 5.25 P− 3 solution for σ =−0.55 . . . . . . . . . . . . . . . . . . . . . . . . 103 5.26 Fast Fourier analysis on P− 3 solutions . . . . . . . . . . . . . . . . . 103 5.27 P− 1 solutions in cell mapping stroboscopic phase plane for different σ (e = 0.1, k 1 = 2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.28 Trajectory inβ−β 0 phase plane (left) andβ response vs. true anomaly (right) of P− 1 solution . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.29 Fast Fourier analysis on P− 1 solutions . . . . . . . . . . . . . . . . . 107 5.30 P− 1 solutions in stroboscopic phase plane for three different orbital eccentricities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 5.31 Phase plane trajectory andβ response vs. true anomaly ofP− 1 solution109 5.32 Fast Fourier analysis on P− 1 solutions . . . . . . . . . . . . . . . . . 110 5.33 Stability analysis at the origin, 0<k 1 < 2 with σ = 0.5, 0,−0.5 . . . . 111 5.34 P− 1 solutions in stroboscopic phase plane (e = 0.1,σ = 0.5) . . . . . . 112 5.35 Stability analysis of equilibrium point, k 1 < 1 (unstable), σ = 0.5 . . . 113 xv List of Figures 5.36 Trajectory in β−β 0 phase plane and β response vs. true anomaly of P− 1 solution for 1·T, k 1 = 0.5 (e = 0.1, σ = 0.5) . . . . . . . . . . . 114 5.37 Trajectory in β−β 0 phase plane and β response vs. true anomaly of P− 1 solution for 1·T, k 1 = 0.75 (e = 0.1, σ = 0.5) . . . . . . . . . . . 115 5.38 P− 1 solutions in stroboscopic phase plane for k 1 > 1 (e = 0.1, σ = 0.5 )116 5.39 Trajectory in β−β 0 phase plane and β response vs. true anomaly of P− 1 solution for 1·T, k 1 > 1 (e = 0.1, σ = 0.5) . . . . . . . . . . . . 117 5.40 Eigenvalues evaluated at the origin 1.25<k 1 < 1.27 (left: magnified on the crowed region near λ =−1, right: full scale of unit circle) . . . . . 118 5.41 Trajectory in β−β 0 phase plane and β response vs. true anomaly of P− 2 solution for k 1 = 1.25 (e = 0.1, σ = 0.5) . . . . . . . . . . . . . . 119 5.42 Fast Fourier analysis on P− 2 solutions, k 1 = 1.25 (e = 0.1, σ = 0.5) . 119 6.1 Geometry of coupled orbit-attitude dynamics model in the three body problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 xvi List of Tables 4.1 Computed coefficients of the point mapping G for e = 0.01 and k 2 = 0.1 57 4.2 Periodic solutions’ stability check by computing eigenvalues of matrix H 59 5.1 Samples of P− 1 solution for σ = 0.4 . . . . . . . . . . . . . . . . . . . 82 5.2 Unstable P− 1 solution in Fig.5.9 . . . . . . . . . . . . . . . . . . . . 87 5.3 Eigenvalue pairs λ i (i = 1,...4) for the origin with varying σ showing P− 1→ P− 1 bifurcation . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.4 Bifurcated P− 1 solutions for σ =−0.58 . . . . . . . . . . . . . . . . . 93 5.5 Eigenvalue pairs λ i (i = 1,...4) for the origin with varying σ showing P− 1→ P− 2 bifurcation . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.6 P− 2 solutions for σ =−0.51 to−0.53 . . . . . . . . . . . . . . . . . 97 5.7 Eigenvaluesλ 1,...4 for the origin with varying σ showingP− 1→P− 3 bifurcation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.8 P− 3 solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.9 Examples of P− 1 solutions for larger σ (e = 0.1, k 1 = 2) . . . . . . . . 105 5.10 Examples of P− 1 solutions for larger e, σ = 0.5 . . . . . . . . . . . . . 109 5.11 Examples of P− 1 solutions for different k 1 , σ = 0.5 . . . . . . . . . . . 113 5.12 Examples of P− 1 solutions for different k 1 , σ = 0.5 . . . . . . . . . . . 118 xvii 1 Introduction Spacecraft attitude control to orient as desired has been accomplished via various ways including spin stabilization and three-axis stabilization (Hughes [36]). We are interested in the attitude control that exploits the characteristics of system dynamics for its benefits such as reduced power consumption or reletively simple control system hardware. There are two common methods for this: gravity-gradient stabilization and single-axis spin stabilization. The gravity-gradient satellites are stabilized about the local verical line to point to the Earth. By having this orientation of the spacecraft, pointing to a primary planet, communication ability beween the two is significantly improved. Spin-stabilization is an effective way (Gerlach [35]) to maintain the one axis (spin-axis) pointing in a required direction and control spacecraft’s orientation in space. Spin-stabilized satellites resist perturbations and require small control effort to maintain the orientation. Moreover, the spacecraft provide a continuous sweeping motion which is a requirement in many missions. Gravity-gradient attitude stabilization was used in-orbit since 1963. The first to successfully be gravity-gradient stablized was 1963-22A (Fischell [39]). The mission requirement for the satellite 1963-22A to have the antennas directed to the Earth within 20 degrees of the local vetical has been accomplished. The DODGE (Depart- ment of Defense Gravity Experiment) satellite assessed the potential of three axis 1 Introduction stabilization in geostationary orbit. Moreover, Salyut 6 had spent some time in a gravity gradient stabilized mode. More are described in Huges [36]. However, there still are pointing accuracy issues and three-axis stabilization has not functioned as intended. Consequently knowledge of gravity gradient torque and its variation in time are important in attitude control system designs for satellites. Spin-stabilization is the most commonly used method to maintain the spin axis in the desired direction, for example NASA’s Pioneer 10 and 11 (Fimmel et al. [42]), and Galileo Jupiter orbiter spacecraft and its atmospheric probe (Fischer [11]). There has been many applications with combination of spin-stabilization scheme and the others as well, including Garcia el al [17], Frost [18], Roldugin et al. [15], and Rarquhar [45]. To perform effective attitude control of both gravity-gradient and spin stabilized satellites, itiscrucialtounderstandthedynamicbehaviorofthesystem. Thedynamics of this class of satellites exhibits complex behavior. Attitude dynamics of the satellite couples with orbital dynamics resulting the dynamic behavior which includes stable and unstable periodic solutions with different periods that vary and bifurcate with changing system parameters. This complex and parameter dependent dynamics has a great influence on the design of the spacecraft and its attitude control system. Past studies contain several analytical approaches to the problem of analyzing the attitude dynamics in a slightly elliptical orbit. In most works, linearization was per- formed to find closed form solutions assuming small angular deviations in order to apply some form of perturbation analysis. In some cases, the orbital eccentricity was set to be small, and the method of harmonic balance was applied. Accordingly the analysiswasapplicableonlytolocalmotionabouttheequilibrium. Computingfamilies of periodic solutions is a straight forward numerical computation, such as integration using boundary conditions and Poincare map analysis. However, the randomly cho- 2 Introduction sen initial conditions did not exhaust all motion possibilities and therefore could not produce a reliable global motion representation. Existing methods cannot efficiently analyze important nonlinear modes of behavior, such as regions of attraction and oc- curing bifurcations. Since the exact equations of motion for attitude dynamics of a satellite and astrody- namics models are coupled, nonlinear, and time varying, analytic solutions are rarely possible. In this research we propose an approach that will enable global analysis of the attitude dynamics models. The proposed approach utilizes cell mapping and point mapping (Poincar` e map) techniques that form a combination of numerical and analyt- ical methods. A point mapping method is applied using a discrete time representation of system dynamics. Computing equilibrium points of the mapping corresponds to finding a periodic solution of the continuous time representation (see Flashner [23]). Stability of the solution can then be established analytically. The cell mapping method is used to analyze the global behaviour of the system (see Hsu [9, 10]). In the cell map- ping approach, the state variables are thought of as a collection of very small intervals. The cell state space that we are interested in is constructed by dividing each state variable component into uniform sizes, respectively. In the cell state space, a cell-to- cell mapping dynamics of the system is created. Solutions of any period of interest are foundwhenthecellrepeatsafterapplyingthemapK-times. Theunravellingalgorithm is used to find the periodic cells and a modification of cell mapping algorithm is used to compute the periodic solution and its stability accurately. Fast Fourier analysis is used to further characterize the motion of invariant surfaces surrounding the fundamental periodic solutions. The proposed approach has the following features: • Small parameter assumption is not needed. • The analysis is global. 3 1.1 Gravity Gradient Satellite • Regions of attraction can be computed. • All three axes can be analyzed simultaneously. To demonstrate the approach, two different classes of problems are considered: 1) pitch motion of a gravity gradient satellite in an elliptic orbit, 2) attitude dynamics of a spinning satellite in an elliptic orbit. This proposal is organized as follows. In the rest of Chapter 1, background of gravity gradient satellite and spinning satellite is introduced. In Chapter 2, the background for modeling attitude dynamics is derived. Defining orbitalelementsandreferenceframesareengagedforthemotionofanorbitingsatellite. Chapter3developstheanalysismethodsincludingpointmapping, cellmapping, and extended cell mapping. Solving for periodic solutions and properties of the solutions are explicitly explained. Chapters4and5discusstheresultsofanalysisofpitchmotionandspinningsatellite, respectively, using the methods from Chapter.3. Some of results are repeated from the literature to demonstrate our methods. New periodic solutions and invariant surfaces are revealed and bifurcation diagrams are computed. Approached methods and results are summarized and future work are discussed in Chapter 6. 1.1 Gravity Gradient Satellite ThestudyoftheeffectofgravitationaltorqueonaspacecraftwasinitiatedbyRoberson [40] using potential methods. The model of the torque was expanded to the first- order expressions in the ratio of satellite dimension to its distance from the center of attraction. Roberson [40], Nidey [38], and Hultquist [37] expressed the torque on a 4 1.1 Gravity Gradient Satellite rigid body in an inverse square field. Analysis of gravity gradient stabilized satellites was performed by assuming circular orbits. The circular orbit assumption lead to the time-invariant dynamic model (Kaplan [32], Hughes [36], and Sidi [33]). In addition, the motion of gravity gradient stabilized satellite models were developed for a small angle motion that together with the circular orbit assumption leads to a linear time- invariant dynamic representation (Kaplan [32], Sidi [33])). The representation allowed for classical analysis using Laplace and frequency domains. Moreover, in this form the pitch axis was decoupled from the roll and yaw axes allowing two separate control designs (see Wertz [27]). A considerable research effort was done to analyze pitch motion of a satellite in an elliptic orbit, extensively in the 1960’s. Pitch motion is most important in the design of a central-body pointing satellite. A stabilized satellite tumbles once an orbital pe- riod, while pointing to the central-body. This problem was first considered by Baker [41]. He pointed that the linearized pitch motion could be represented by a Mathieu equation with an added forcing term. Moran [26] employed a perturbation technique to study the planar motion of an idealized dumbbell satellite in a circular orbit. He showed that resonance conditions in the radial equation of motion can arise for certain values of the initial librational frequency. He also showed that the differential equa- tion for the first radial correction term gives rise to a divergent oscillatory solution as time increases. Schechter [20, 21] discussed the planar motion of dumbbell shaped satellite in an elliptic orbit, presenting the solution as a power series in eccentricity. Zlatoustov et al. [50] studied families of periodic solutions of the pitch oscillations by the method of numerical integration and investigated stability in linear approximation. Beletskii[55]derivedpitchmotionofagravitygradientsatelliteandshowedfundamen- tal analytic solutions for small eccentricities and small amplitude motion. Branches of 5 1.1 Gravity Gradient Satellite periodic solutions that depend on eccentricity and the body’s mass distribution were also discussed. Modi et al. [53] studied periodic solutions in a circular and elliptic orbits by employing the method of harmonic balance and boundary condition prob- lems. Following this study, stability of the periodic solutions was analyzed using linear perturbation analysis in Modi et al. [54]. Zlatoustov et al. [49] investigated stability of the period of 2π periodic solutions in nonlinear mode. Sarychev et al. [48] investigated the branching of periodic solutions (period of 2π, 4π, and 6π solutions) of the motion of the satellite. Simple boundary value problem was applied for slightly elliptical or- bits and numerical computions was used for analysis of large eccentricities. Zevin [1] obtained uniqueness conditions for 2π-periodic solution for the pitch osillations of the satellite qualitatively. Recent years, studies of pitch motion of gravity satellite utilized nonlinear analysis using Poincar 0 e maps and Melnikov’s analytic method to study the routes to chaos, presenting phase planes, bifurcation diagrams, and chaos plots (see Tong et al. [57] and Karasopoulos [19]). The global behavior of the pitch motion was successfully investigated in Koh [13]. Three-axis motion was also studied to analyze the stability of the pitch motion with out-of-plane perturbations. In his dissertation, Fleig [3] investigated the pitch motion as well as three dimensonal attitude motion of a gravity gradient spacecraft. The pitch motion was analyzed using linearization to form a Mathieu equation for motion with small eccentricity and applying asymptotic expansion method for motion with large eccentricities. In this work, three-axis motion was also showed by utilizing the averaging method. Kane [46] applied Floquet theory to three dimensional analysis and discovered instabilities in the nominally stable range of spacecraft parameters. Wallace et al. [16] showed the same phenomena with nonlinear equations using Liapunov direct method. Breakwell et al. [28] explained this phenomenon in terms of a nonlinear 6 1.2 Spinning Satellite resonance between the pitch and yaw frequencies. Beletsky et al. [56] and Celletti[2] proved numerically and analytically that the pitch periodic motions are unstable with respect to out-of-plane perturbations. 1.2 Spinning Satellite The studies of attitude dynamics of spinning satellites have mostly shown how stability is related to the rate of spin and to the inertia properties of a satellite. A linearized model of a spinning satellite in an elliptic orbit was first formulated by Kane et al. [47], and stability conditions were studied. Wallace et al. [16] developed a nonlinear model that investigated the effect of additional nonlinear terms on the stability analysis of the linearized model in Kane et al. [47]. Modi et al. [51] utilized cross-section analysis to find periodic solutions using randomly chosen initial conditions. This method was limited to circular orbits and by the fact that the initial conditions were randomly chosen, it did not exhaust all possible modes of behavior. McGill et al. [12] considered unsymmetrical rigid body to examine the classical stability problem. Modi et al. [52] included solar radiation pressure torques to analyze attitude motion of spinning satellite. The resonance analysis suggested the existence of critical combinations of system parameters leading to large amplitude oscillations. Moreover, Hitzl et al. [14] studiedtheattitudestabilityofaspinningsatelliteintherestrictedthreebodyproblem using the averaging method and generated stability charts for the trivial equilibrium point. Cole et al. [29] employed a linearized model to develop a modal control scheme to stabilize a spinning satellite. This study showed that nonlinear effects dominate the behavior of the closed-loop system, and therefore linear control design approach for stabilizing a spinning satellite is unsatisfactory. 7 2 Model of Satellite in Orbit 2.1 Orbit definition Consider a satellite in an elliptic orbit. The position of the satellite in the orbit is defined by a distance vector ¯ R c from the Earth center to the body center of mass and true anomaly v as shown in Fig.2.1. The elliptic orbit is characterized as follows. Figure 2.1: Elliptic Keplerian orbit characters: a,b, and p The semi-major and semi-minor axis are denoted by a and b respectively as shown in Fig.2.1. The orbital eccentricity denoted by e characterizes the oblateness of the orbit, and it is defined as e = v u u t 1− b a ! 2 . 8 2.1 Orbit definition (Eccentricity 0<e< 1 for eccentric orbit, e = 0 for circular orbit, e = 1 for parabolic trajectory, and e> 1 for hyperbolic trajectory) In addition, p is the notation for semi-latus rectum as shown in Fig.2.1. It has the relation with a and e as p =a(1−e 2 ). It is generally accepted to use a and e as the orbital elements defining the shape and size of the orbit (see for example Bate et al. [43] and Kaplan [32]). Figure2.2: Orbital elements to define position of the orbit in space w.r.t. the inertial frame: i, Ω,ω,and ν The orientation of the orbit with respect to the inertial frame is defined by the following four parameters (see also Bate et al. [43]) as shown in Fig.2.2: • Inclination(i) defines the orientation of the orbit with respect to the Earth’s equator. 9 2.2 Reference Frames • Longitude of the ascending node(Ω) defines the location of the ascending and descending orbit locations with respect to the Earth’s equatorial plane. • Argumentofperiapsis(ω)defineswherethelowpoint, perigee, oftheorbitiswith respect to the earth’s surface as an angle measured from the ascending node to the periapsis. • True anomaly(ν) defines the position of the satellite in orbit at each time with respect to periapsis direction. 2.2 Reference Frames We shall use the notation [¯ e F ] to denote the frame F and e F i where i = 1, 2, 3 to refer to the i th -axis unit vector of a Cartesian coordinate of the frame F. The motion of a satellite can be expressed in different frames listed in the following. Inertial Frame The inertial frame [¯ e I ] is known as ECI, a non-rotating Earth centered inertial frame. As shown in Fig.2.2, the unit vector e I 1 is permanently fixed in a direction relative to the celestial sphere in the Earth’s equatorial plane, e I 3 points the celestial north pole, and e I 2 completes the coordinate system in the Earth’s equatorial plane. 10 2.2 Reference Frames Orbital Inertial Frame and Local Orbital Frame Figure 2.3: Relation between the orbital inertial frame and the local orbital frame The orbital inertial frame [¯ e i ] and the local orbital frame [¯ e O ] are shown in Fig.2.3. The orbital inertial frame [¯ e i ] is located at the Earth’s center and non-rotating. The unit vector e i 1 is fixed in the direction of periapsis, e i 3 is normal to the orbital plane, the direction of the orbit’s north pole, ande i 2 completes a right handed coordinate set. The local orbital frame [¯ e O ] is at the Earth’s center with the unit vector e O 1 always in the direction of a satellite. Note that the direction of the unit vector e O 1 is called the local vertical. The normal vector e O 3 is same as e i 3 , and e O 2 completes a right handed coordinate set. 11 2.3 Transformation between Coordinate Frames Reference Frame and Body Frame Figure 2.4: Relation between the local orbital frame and the reference frame The reference frame [¯ e R ] is used to describe satellite attitude with respect to the orbital frame. Itsoriginisatthesatellitemasscenter, andtheunitvectore R 3 isinthedirection oflocalverticalpointingtoEarth,e R 1 pointsinthedirectionoftravel, ande R 2 completes the right handed orthonormal triad. The body frame [¯ e B ] is fixed to the body mass center. 2.3 Transformation between Coordinate Frames The various frames are related to each other by rotational transformations. In the following, these transformations are defined using Euler rotation sequences. 12 2.3 Transformation between Coordinate Frames 2.3.1 Euler angle sequence An Euler rotation sequence is defined using principal rotations. The corresponding rotation matrices are C 1 (ψ) =         1 0 0 0 cosψ sinψ 0 − sinψ cosψ         , C 2 (θ) =         cosθ 0 − sinθ 0 1 0 sinθ 0 cosθ         , C 3 (φ) =         cosφ sinφ 0 − sinφ cosφ 0 0 0 1         . C 1 , C 2 , and C 3 denote rotations about 1−, 2−, and 3−axis, respectively. Thesequenceoftherotationsrepresentstheorientationofanytwosuccessiveframes. For example, the transformation from the reference to the body frame is denoted by C(θ 1 ,θ 2 ,θ 3 ) , and given by: B C R =C 1 (θ 1 )C 2 (θ 2 )C 3 (θ 3 ) (2.1) =         cosθ 2 cosθ 3 cosθ 2 sinθ 3 − sinθ 2 sinθ 1 sinθ 2 cosθ 3 − cosθ 1 sinθ 3 sinθ 1 sinθ 2 sinθ 3 + cosθ 1 cosθ 3 sinθ 1 cosθ 2 cosθ 1 cosθ 2 cosθ 3 + sinθ 1 sinθ 3 cosθ 1 cosθ 2 sinθ 3 − sinθ 1 cosθ 3 cosθ 1 cosθ 2         2.3.2 Transformation from Inertial to Orbital Inertial Frame This transformation is defined by longitude of ascending node Ω, argument of periapsis ω, and inclination i. The sequence of the rotations is denoted by arrows the axis and angle of the rotation as follows, 13 2.3 Transformation between Coordinate Frames [¯ e I ] e I 3 ,Ω −−→ [¯ e a ] e a 1 ,i −−→ [¯ e b ] e b 3 ,ω −−→ [¯ e i ] where [¯ e a ] and [¯ e b ] are auxiliary frames for intermediate steps through the transfor- mation. Then the corresponding rotation matrix from inertial frame [¯ e I ] to orbital inertial frame [¯ e i ] is i C I =C 3 (ω)C 1 (i)C 3 (Ω). 2.3.3 Transformation from Orbital Inertial to Local Orbital Frame The local orbital frame is defined with respect to the orbital inertial frame by true anomaly ν about e i 3 . The unit axis e O 1 of the local orbital frame is pointing to the satellite as shown in Fig.2.3. Consequently, the rotation matrix from the orbital inertial frame to the local orbital frame is O C i =C 3 (ν). Orbital angular velocity with respect to the inertial frame, I ¯ ω O , expressed in the orbital frame is given by I ¯ ω O = ¯ e O 1 ¯ e O 2 ¯ e O 3         0 0 | I ¯ ω O |         . (2.2) 2.3.4 Transformation from Local Orbital to Reference Frame The reference frame [¯ e R ] is introduced for determining satellite’s attitude with respect to the local orbital frame. The rotation matrix from the local orbital frame to the 14 2.3 Transformation between Coordinate Frames reference frame (see Fig.2.4) is R C O =         0 1 0 0 0 −1 −1 0 0         . (2.3) 2.3.5 Transformation from Reference to Body Frame The orientation of the body with respect to the reference frame is determined by three Euler’s angles (θ 1 ,θ 2 ,θ 3 ). By convention, the angles are called (roll,pitch,yaw), respectively. Euler angle sequence for the position of the body axes is [¯ e R ] e R 3 ,θ 3 −−−→ [¯ e c ] e c 2 ,θ 2 −−−→ [¯ e d ] e d 1 ,θ 1 −−−→ [¯ e B ], (2.4) and the corresponding rotation matrix B C R is given in Eq. (2.1). The transformation from the local orbital frame to the body frame is given by multiplication, B C R · R C O , 15 2.3 Transformation between Coordinate Frames B C O = B C R · R C O (2.5) =C 1 (θ 1 )·C 2 (θ 2 )·C 3 (θ 3 )·         0 1 0 0 0 −1 −1 0 0         (2.6) =         sinθ 2 cosθ 2 cosθ 3 − cosθ 2 sinθ 3 − sinθ 1 cosθ 2 sinθ 1 sinθ 2 cosθ 3 − cosθ 1 sinθ 3 − sinθ 1 sinθ 2 sinθ 3 − cosθ 1 cosθ 3 − cosθ 1 cosθ 2 cosθ 1 cosθ 2 cosθ 3 + sinθ 1 sinθ 3 − cosθ 1 cosθ 2 sinθ 3 − sinθ 1 cosθ 3         . (2.7) To define body’s angular velocity, we use the notation R ¯ ω B B for the body’s angular velocity with respect to the reference frame [e R ] expressed in the body frame [e B ]. The total angular velocity can be obtained by a sum of angular velocities of each rotation in Eq. (2.4) as R ¯ ω B B = R ¯ ω c c + c ¯ ω d d + d ¯ ω B B , where the notation vector α ¯ ω β β represents the angular velocity of the frame β with respect to the frameα expressed in the frameβ. Each angular velocity is expressed in the body frame so that we have: R ¯ ω B B = ¯ e B 1 ¯ e B 2 ¯ e B 3         − ˙ θ 3 sinθ 2 +0 + ˙ θ 1 ˙ θ 3 sinθ 1 cosθ 2 + ˙ θ 2 cosθ 1 +0 ˙ θ 3 cosθ 1 cosθ 2 − ˙ θ 2 sinθ 1 +0         (2.8) Note that when the angular velocity is expressed in the body frame and it is obvious, 16 2.4 Gravity Gradient Torque the subscript B is omitted. 2.4 Gravity Gradient Torque Consider a spacecraft B is moving in a gravitational field as shown in Fig.2.5, see also NASA [34] and Hughes [36]. E a r t h R B d m Gravitational Field r R c Figure 2.5: Spacecraft B in a gravitational field The following four assumptions simplify the gravity-gradient torque expressions (see Hughes [36]) : 1. The gravitational field is uniform. 2. There is only one spherically symmetrical primary body, specifically the Earth in our case. 17 2.4 Gravity Gradient Torque 3. The spacecraft is small compared to its distance from Earth. 4. The spacecraft consists of a single body. According to the assumptions, gravitational torque about⊗ (see Fig.2.5) is g c =−μ B ¯ r× ¯ R R 3 dm (2.9) ¯ R = ¯ R c + ¯ r (2.10) where ¯ R c is the position vector of the mass center of spacecraft B with respect to the center of Earth, ¯ r is the location of the mass element dm with respect to⊗, and μ =G·m E where m E = 5.97e24kg μ = 3.986e5N·km 2 /kg. By substituting Eq. (2.10) into Eq. (2.9) and using binomial expansion, we get: R −3 =R −3 c " 1− 3(¯ r· ¯ R c ) R 2 c +O r 2 R 2 c !# . To this degree of accuracy, Eq. (2.9) becomes (see Hughes [36]) g c =− 3μ R 5 c ! ¯ R c × B ¯ r ¯ ·rdm· ¯ R c (2.11) Since ¯ R c is in the direction of (−e R 3 ) in the reference frame, the unit vector ˆ e R 3 = 18 2.5 Euler’s Equation − ¯ R c /| ¯ R c |. Therefore Eq. (2.11) can be written as g c = 3μ R 3 c ! ˆ e R 3 ×I· ˆ e R 3 , (2.12) where I is the [3x3] moment of inertia matrix. Note that gravitational torque is expressed in principal axes, and it depends on the direction of the local vertical, which varies throughout the orbit. 2.5 Euler’s Equation 2.5.1 Global Motion Euler’s equations of motion for a single rigid body expressed in principal axes when the external torque is present is (see Hughes [36]) I 1 I ˙ ω 1 B + (I 3 −I 2 ) I ω BI 3 ω B 2 =g c1 I 2 I ˙ ω 2 B + (I 1 −I 3 ) I ω BI 1 ω B 3 =g c2 (2.13) I 3 I ˙ ω 3 B + (I 2 −I 1 ) I ω BI 1 ω B 2 =g c3 . The gravitational torqueg c is expressed in the body frame by a relatione B 3 = B C R 3 ·e R 3 from Eq. (2.12) as g c = 3 μ R 3 c !         (I 3 −I 2 ) B C R 23 B C R 33 (I 1 −I 3 ) B C R 13 B C R 33 (I 2 −I 1 ) B C R 13 B C R 23         . (2.14) 19 2.5 Euler’s Equation Expressing gravity gradient torque in terms of Euler angles, we get: I 1 I ˙ ω 1 B + (I 3 −I 2 ) I ω BI 3 ω B 2 = 3 μ R 3 c (I 3 −I 2 ) sinθ 1 cosθ 1 cos 2 θ 2 I 2 I ˙ ω 2 B + (I 1 −I 3 ) I ω BI 1 ω B 3 = 3 μ R 3 c (I 1 −I 3 ) cosθ 1 sinθ 2 cosθ 2 (2.15) I 3 I ˙ ω 3 B + (I 2 −I 1 ) I ω BI 1 ω B 2 = 3 μ R 3 c (I 2 −I 1 ) sinθ 1 sinθ 2 cosθ 2 . Thevector I ¯ ω B ofEq. (2.15)denotesangularvelocityofthebodyframewithrespectto theinertialframeexpressedinthebodyframe. Whilethelocalverticalisrotatingalong the orbital path, the body’s orientation is changing. Therefore, dynamical analysis of the governing equation is required. To explicitly express the body’s angular velocity I ¯ ω B , the angular velocity is de- composed into two parts: orbiting motion and rotating motion, i.e. I ¯ ω B B = I ¯ ω O O + R ¯ ω B B . Here, I ¯ ω O O denotes the orbital angular velocity with respect to the inertial frame ex- pressed in the local orbital frame [¯ e O ] as Eq. (2.2), while R ¯ ω B B is the body’s angular velocity with respect to the reference frame expressed in the body frame [¯ e B ] as given in Eq. (2.8). The orbital angular velocity I ¯ ω O O is in the normal direction to the orbital plane which can be written as| I ¯ ω O |· ˆ k where the magnitude of the vector is (see Hughes [36]) | I ¯ ω O | = r μ a 3 (1 +e cosν(t)) 2 (1−e 2 ) 3 2 ! . The angular velocity I ¯ ω O expressed in the body frame denoted by I ¯ ω O B is given as: 20 2.5 Euler’s Equation I ¯ ω O B = [¯ e B ] B C R · R C O         0 0 | I ¯ ω O |         (2.16) = [¯ e B ]         s 2 c 2 c 3 −c 2 s 3 −s 1 c 2 s 1 s 2 c 3 −c 1 s 3 −s 1 s 2 s 3 −c 1 c 3 −c 1 c 2 c 1 c 2 c 3 +s 1 s 3 −c 1 c 2 s 3 +s 1 c 3                 0 0 | I ¯ ω O |         where s and c denote sine and cosine respectively and subscript 1, 2, and 3 are ab- breviations for θ 1 ,θ 2 , and θ 3 . The body’s rotational angular velocity R ¯ ω B given in equation (2.8) is already expressed in the body frame. The body’s angular velocity, finally, is I ¯ ω B = I ¯ ω O B + P ¯ ω B B =         −c 2 s 3 −s 1 s 2 s 3 −c 1 c 3 −c 1 c 2 s 3 +s 1 c 3         | I ¯ ω O | +         1 0 −s 2 0 c 1 s 1 c 2 0 −s 1 c 1 c 2                 ˙ θ 1 ˙ θ 2 ˙ θ 3         , (2.17) or using the kinematic relations Eq. (2.17) is written as 21 2.5 Euler’s Equation         ˙ θ 1 ˙ θ 2 ˙ θ 3         =         1 0 −s 2 0 c 1 s 1 c 2 0 −s 1 c 1 c 2         −1 I ¯ ω B −         1 0 −s 2 0 c 1 s 1 c 2 0 −s 1 c 1 c 2         −1         −c 2 s 3 −s 1 s 2 s 3 −c 1 c 3 −c 1 c 2 s 3 +s 1 c 3         | I ¯ ω O | =         1 s 1 s 2 /c 2 s 2 c 1 /c 2 0 c 1 −s 1 0 s 1 /c 2 c 1 /c 2         I ¯ ω B −         −(s 3 /c 2 )(c 2 2 +s 2 1 s 2 2 +c 2 1 s 2 c 2 ) −c 3 −(s 3 /c 2 )(s 2 1 s 2 +c 2 1 c 2 )         | I ¯ ω O |. 22 3 Analysis Methodology A computational approach that combines point mapping and cell mapping methods, and extended version of them will be customized to investigate the global behavior of problems. The point mapping method is utilized to analyze local stability and bifurcation conditions, and the cell mapping method is utilized for global behavior analysis of various applications. In this Section, implementation of the approach to the general dynamic problems and the computational procedure are presented. 3.1 Point Mapping Analysis Point mapping is a method to solve non-autonomous, non-linear periodic dynamical systemsinamannerofdiscretesequenceoftimeincludingfindingequilibriumsolutions and stability properties (see Flashner et al. [22, 24, 44]). Point mapping idea (Poincar´ e map) was introduced by Poincar´ e. Using a mapping transformation, the dynamics of the system is given in terms of a system of difference equations. Consider a dynamical system described by a set of N ordinary differential equations, ˙ x(t) =f(t,x(t),s) (3.1) wherex∈R N is a state vector, t∈R + denotes time,s∈R L is a parameter vector, 23 3.1 Point Mapping Analysis and f is a real valued vector function periodic in t with a period T. The state of the system is observed at discrete instances every period T unit of time and it sets difference equations of Eq. (3.1). Then the dynamic relationship between the states at t =nT and t = (n + 1)T is defined as point mapping and it can be written as x(n + 1) =G(x(n),s), n = 1, 2,··· (3.2) where x(n + 1) and x(n) are the states of the system at t = (n + 1)T and t = nT, respectively. The point mapping formulation defines the state of the system only at integervaluesof theindependentvariablen. Thecontinuoustimehistory ofthesystem fornT <t< (n + 1)T can be recovered at any stage of the analysis by integrating Eq. (3.1). Based on this concept, in this Section, computing periodic solutions by means of point mapping is addressed. 3.1.1 Periodic Solutions and Local Stability Properties The dynamics of the point mapping consists of equilibrium points and periodic solu- tions. An equilibrium pointx ∗ is given by x ∗ =G(x ∗ ). A periodic solution of period K or in short aP−K solution of the point mapping Eq. (3.2) for somes =s ∗ consists of K distinctP−K pointsx ∗ (j), j = 1,...,K such that x ∗ (j + 1) =G(x ∗ (j),s ∗ ), j = 1,...,K− 1 x ∗ (1) =G K (x ∗ (1),s ∗ ) (3.3) 24 3.1 Point Mapping Analysis (see Hsu [6]). From this equality, P−K solutions Eq. (3.3) for point mapping satisfy *(1) χ Χ 2 Χ 1 χ * *(2) χ *(3) χ Figure 3.1: Schematic phase plane for an equilibrium pointx ∗ and periodic solutions X ∗(i) ,i = 1,··· ,K for K = 3 x ∗ −G K (x ∗ ) = 0. (3.4) By solving the algebraic system Eq. (3.4), P−K solutionsx ∗ are determined. Define z(t) = x(t)−x ∗ (t;s) as the perturbation of the state x about a periodic solution x ∗ . In order to analyze periodic solutions, the system Eq. (3.1) can be expressed as ˙ z(t) =A(t,s)z(t) + ∞ X k=2 r k (t,z(t),s;x ∗ ), (3.5) where the matrix A(t)∈R N×N is given by 25 3.1 Point Mapping Analysis A(t,s) = " ∂f(t,x,s) ∂x # x=x ∗ , and r k (t,z(t),s;x ∗ ) is a vector of all polynomials of degree k in the components of z(t). Then the point mapping formulation of equation (3.5) is z m+1 =H(s)z m + ∞ X i=2 q i (z m ,s;x ∗ ), m = 1, 2,··· where H(s)∈R N×N is given by H(s) =H K (s)H K−1 (s)···H 1 (s), and H j (s) = " ∂G(x,s) ∂x # x=x ∗(j) . (3.6) Then, the local stability of a P−K solution is determined by the eigenvalues of matrix H(s) (see Flashner [22]). 1. A P−K solution of Eq. (3.2) is locally asymptotically stable if and only if all the eigenvalues of H(s) have moduli less than unity. 2. A P−K solution of Eq. (3.2) is unstable if at least one of the eigenvalues of H(s) has modulus larger than unity. 3. If some of the eigenvalues of H(s) have moduli equal to unity and all other eigenvalues havemoduli less than one, thelinear stability analysisis inconclusive. 26 3.1 Point Mapping Analysis 3.1.2 Bifurcation Conditions The bifurcation conditions of a given solution are discussed. If there exists an integer M such that det(I−H M ) = 0, then a bifurcation of aP−K solution to aP−MK solution may occur. (see Flashner et al. [24]) 1. Bifurcation of P−K to P−K: det(I−H) = 0. 2. Bifurcation of P−K to P− 2K: det(I +H) = 0. 3. Bifurcation of P−K to P−MK: det(I +H M/2 ) = 0, M :even, M/2 :even; det(I−H +H 2 −H 3 ··· +H M/2−1 = 0, M :even.M/2 :odd; det(I +H 2 +H 3 ··· +H M−2 = 0, M :odd. (3.7) 27 3.1 Point Mapping Analysis 3.1.3 Truncated Point Mapping Method Since the functionf(t,x(t),s) can be expressed by a series of the form as (see Flashner et al. [24]) f(t,x,s) =p(t,x,s) = M X in=1 b i 1 ···i N (t,s)x i 1 1 ···x i N N wherep(t,x,s) denotes a vector polynomial of degree M in the state variable x i , i = 1, 2,··· ,N. Note that a polynomial of degree M contains terms with sum of powers of x i (for i = 1, 2,··· ,N) less or equal to M. To abbreviate notation denote this set of indices i 1 ,··· ,i N by E r . The algorithm for point mapping is based on multinomial functions by truncating at some degree k, truncatedpointmapping. Letp(t,x,s) be a polynomial of degree greater than k. Define the truncation operation at degree k by p k (t,x,s) =p(t,x,s) k = k X in=1 b Er (t,s)x i 1 1 ···x iN N to be the polynomial formed by taking only those terms whose degree is k or less. The dependence on parameters of the coefficients can also be obtained as p(t,x,s) k,l = k X in=1 l X jn=1 σ Er,Eq (t)x i 1 1 ···x i N N s j 1 1 ···s j L L . (3.8) The algorithm devised separately keeps track of the powers of the components of the parameter vectors a polynomial degree up tol, in addition to the tracking of the state vectorx. 28 3.1 Point Mapping Analysis The numerical algorithm is based on the fact that the operations of truncation and telescoping on polynomials can be interchanged. Moreover, the Runge-Kutta method of integration can be expressed as a sequence of polynomial telescoping and truncation operations which finally results in a truncated polynomial expression for the point mappingG(x,s) of equation (3.2). To derive the algorithm, the period T is divided into n 0 sub-intervals of length h = T/n 0 . Then, the numerical algorithm starts with the formulation of the Runge-Kutta method (RK4) as x(t p +h) =x(t p ) +h M X m=1 d m k m (t,x,s) (3.9) whereM is the order of the Runge-Kutta method (here,M = 4),h is a fixed time step, and d m are certain constants determined by the Runge-Kutta method. The vectors k m for the differential equations are given by k m (t,x,s) =f(t p +ha m ,x(t p ) +hc m k m−1 ,s), m = 1, 2,··· ,M The constants a m and c m are given by the particular Runge-Kutta method. Letx(t p ) =x(nT + (i− 1)h) =x with 1≤ i≤ n 0 . Through the iteration of the integration, weareinterestedincomputingtruncatedpointmappingsothatk m (t,x,s) takes the form k m =p m (t,x,s) k,l =p k.l m 29 3.2 Cell Mapping Method as Eq. (3.8). Thus, consider Eq. (3.9) and observe that x(t p +ih) k,l =x + M X m=1 hd m p k,l m (i) = k,l X β Er,Eq (i)x i 1 1 ···x i N N s j 1 1 ···s j L L =Q k,l (i;x,s) Computation of the coefficients of the truncated point mapping over the period T is done by repeating (n 0 − 1) times to get x((n + 1)T ) =Q k,l (n 0 ; (Q k,l (n 0 − 1; (···Q k,l (1; (Q k,l (0;x,s))··· ))) k,l =G k,l (x,s) where Q k,l (0;x,s) =x. The corresponding truncated point mapping is exact up to thek th order for the state vector andl th order for the parameters of multinomial terms. 3.2 Cell Mapping Method A method for global analysis of nonlinear dynamical systems based on the idea of point mapping is introduced by Hsu [7, 8]. Because of accuracy limitations of numerical evaluations of point mapping, here the state variable is thought of as a collection of very small intervals. Consider a dynamical system governed by ˙ x(t) =f(t,x(t)) (3.10) wherex is a real-valued N-vector andf is a real-valued vector function periodic in t with a periodT. The cell state space is constructed by setting a bounded region from 30 3.2 Cell Mapping Method x 1 x 2 (L) x 1 (U) x 1 (L) x 2 (U) x 2 h 1 h 2 … N C-1 N C 1 2 3 … 1 2 3 N c1 N c2 1 2 3 … … (sink cell) (sink cell) Figure 3.2: Construction of a cell state space for N = 2 x L i to x U i , where L and U denote lower and upper bound, respectively. Then, divide each state variable component x i , i = 1, 2,··· ,N into a large number of intervals with uniform size h i (see Fig.3.2 an example of a two dimensional cell state space.) The interval z i of the x i -axis is defined to contain the state variables x i such that z i − 1 2 h i ≤x i < z i + 1 2 h i where z i is an integer. Each cell is considered as a separate entity identified by a cell number, and the state space S is regarded as a collection of all cells. If N c i denotes the number of cells along the x i -axis, then the cell space S contains a total of N c =N c 1 ×N c 2 ×···×N c N cells as shown in Fig.3.2. 31 3.2 Cell Mapping Method In the cell state space S, one can define a cell-to-cell mapping C of a dynamical system in the form z(n + 1) =C(z(n)),C :S→S,z(n) = 1, 2,··· ,N c . (3.11) It describes the evolution of a cell dynamical system inN-dimensional cell state space S. For a detailed treatment of cell mapping, an algorithm which allows us to determine the equilibrium states, periodic motions and their domains of attraction when they are asymptotically stable is used. 3.2.1 Periodic Motions and Domains of Attraction The dynamics of the cell mapping is characterized by singular cells consisting of equi- librium cells and periodic cells. An equilibrium cellz ∗ is given by z ∗ =C(z ∗ ). To define periodic cells, let C m denote the cell mapping C applied m times with C 0 understood to be the identity mapping. A sequence of K distinct cellsz ∗ (j), j = 1, 2,··· ,K which satisfy z ∗ (m + 1) =C m (z ∗ (1)), m = 1, 2,··· ,K− 1, z ∗ (1) =C K (z ∗ (1)) is a P-K motion. The domain of attraction of a P−K motion is determined using a r-step domain. If a cellz is mapped in r-steps into one of the P−K cells of theP−K motion, then 32 3.2 Cell Mapping Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x 1 x 2 (L) x 1 (U) x 1 (L) x 2 (U) x 2 (sink cell) z* *(1) z *(2) z *(3) z Figure 3.3: Schematic cell state space showing an equilibrium cellz ∗ andP−K cells z ∗(i) ,i = 1,··· ,K for K = 3 any further mapping will be locked in this P−K motion. The set of all cells which are r-steps or less removed cells from P−K motion is called the r-step domain of attraction. If the total domain of attraction is with r→∞ then it is associated with the asymptotically stable equilibrium states or periodic solutions. 3.2.2 Global Analysis In order to obtain the cell mapping associated with the periodic dynamical system Eq. (3.10), the following procedure is implemented. First, construct a cell space structure in the x i -direction with cell size h i . Let x i (t p ) denote the center point of the cell z p so that x i (t P ) =h i z i (t p ), i = 1, 2,··· ,N. Integrate the state equation (3.10) over one 33 3.2 Cell Mapping Method period of time with initial conditionsx(t p ) ofz p and find which cell it maps to. This process of finding image cells is repeated for each cell in the cell state spaceS, Fig.3.2. Once the mapping process is done, the properties of the cells are determined by an unraveling algorithm (see Hsu [8]). To delineate the global properties four terms are used, (1) Sink cell, (2) Group number, (3) Periodicity number, and (4) Step number: • Sink cell : In application of cell-to-cell mapping, one is interested in a bounded region of the state space which contains a finite number of cells. Once the state variable exceeds the region, it is no longer considered for further evolution of the system. It is called the sink cell. By definition the sink cell is a P− 1 cell. • Group number : The group numbers are positive integers which are assigned sequentially as the periodic motions are discovered. Each group has an invariant set in the form of a periodic motion which shares the same period. • Periodicity number : The periodicity number indicates the period of the cells in the group. • Step number : The step number shows how many steps to map each cell into a periodic cell as discussed in the domain of attraction Subsection 3.2.1. If the step number ofz is 0, it happens to be a periodic cell. Thus the global analysis is to determine those properties for every cell when the cell mappingC(z) is given. 34 3.3 Extended Cell Mapping Method 3.3 Extended Cell Mapping Method 3.3.1 Computation of P-1 Solutions Assume that a cellj is aP− 1 cell. Then the center of the cell is mapped back to the same cell to location φ ∗ j (T ), i.e. φ ∗ j (T ) =G(φ ∗ j (0)) as shown in Fig.3.4. However, the center of the cellφ ∗ j (0) is not an exactP− 1 point. We want to refine the center point using point mapping with a parameter z to find an exact P− 1 solutionz ∗ . The idea of extended cell mapping method is based on Golat et al. [31]. th j cell φ * ( 0 ) j φ * ( 1 ) j z * th j cell z Figure 3.4: P− 1 solutionz ∗ of cell j with small perturbation z Consider a case shown in Fig.3.5. The center point of cellj,φ ∗ j (0) is mapped to the neighbor cell, and it is again mapped to the cell next to itself as shown in Fig.3.5. We conjecture that there is a P− 1 solution close to φ ∗ j (0), and we are going to solve for it using point mapping method with a parameter z. Mathematically it means that we want to find location z such that 35 3.3 Extended Cell Mapping Method (Sink Cell) th j cell φ *(1) j φ *(0) j+2 φ *(0) j+1 φ *(2) j th j cell φ *(0) j z* z Figure 3.5: Neighborhood of cell j might have aP− 1 Solutionz ∗ with small pertur- bation z φ ∗ j (0) +z =G(φ ∗ j (0) +z) (3.12) Expansion of the function G in Eq. (3.12) in Taylor series about the solution φ ∗ j (0) yields φ ∗ j (0) +z =C(φ ∗ j (0)) + ∂G ∂z | φ ∗ j (0) ·z + ∂ 2 G ∂z 2 | φ ∗ j (0) ·z 2 +··· =φ ∗ j (T ) +H·z +J·z 2 +··· =φ ∗ j (T ) +G j(m) (z). To find a P− 1 solution we need to solve the equation 36 3.3 Extended Cell Mapping Method z−G j(m) (z) =φ ∗ j (T )−φ ∗ j (0). (3.13) Using the linear part of Eq. (3.13), then it can be written as (I−H)z =φ ∗ j (T )−φ ∗ j (0), and z = (I−H) −1 · (φ ∗ j (T )−φ ∗ j (0)) (3.14) This iteration process can continue in the neighborhood of cellj. Note that Eq. (3.14) will have a singularity problem if we are at a bifurcation point. 3.3.2 Computation of P-K Solutions Assume that cells j 1 ,··· ,j K constitute a P−K solution, i.e. φ ∗ j 1 G(φ ∗ j 1 ) −−−−→φ ∗ j 2 G(φ ∗ j 2 ) −−−−→··· φ ∗ j K G(φ ∗ j K ) −−−−→φ ∗ j 1 where φ ∗ j i is the center of the cell j i , i = 1,··· ,K. As shown in Fig.3.6, the cell j i ,i = 1,··· ,K is mapped back to itself in period K·T. We now look at the perturbationz 1 in cellj 1 that will generate perturbationz i , i = 2,··· ,K in the other cell that will give a precise P−K solution. By definition of a 37 3.3 Extended Cell Mapping Method j 2 φ * ( 0 ) j 2 G ( φ * ( 0 )+z ) j1 j 1 1 j 1 φ * ( 0 ) j 1 j K φ * ( 0 ) j K z 1 z 2 z K G ( φ * ( 0 )+z ) j2 j 2 2 j 1 j 2 j K (K) z* (1) z* (2) z* G ( φ * ( 0 )+z ) jk j k k Figure 3.6: Cell j i are P−K cells with perturbation z i where i = 1,··· ,K P−K solution we have G j ! (φ ∗ j 1 (0) +z 1 ) =φ ∗ j 2 (0) +z 2 G j 2 (φ ∗ j 2 (0) +z 2 ) =φ ∗ j 3 (0) +z 3 . . . G j K−1 (φ ∗ j K−1 (0) +z K−1 ) =φ ∗ j K (0) +z K G j K (φ ∗ j K (0) +z K ) =φ ∗ j 1 (0) +z 1 38 3.3 Extended Cell Mapping Method or in a compact notation; G j i (φ ∗ j i (0) +z i ) =φ ∗ j i+1 (0) +z i+1 , i = 1,··· ,K− 1 G j K (φ ∗ j K (0) +z K ) =φ ∗ j 1 (0) +z 1 (3.15) where we are going to solve forz 1 ,··· ,z K . Hence we haveN·K equations withN·K unknowns where N is the number of states. Using Taylor series we have G j i (φ j i (0) +z i ) =G j i (φ ∗ j i (0)) + " ∂G ∂x # φ ∗ j i z i +··· (h.o.t.) =φ ∗ j i (T ) +H j i z i +··· (h.o.t.) (3.16) =φ ∗ j i (T ) +G j i (m) (z), i = 1,··· ,K and substitute Eq. (3.16) into Eq. (3.15) yields: z i+1 =φ ∗ j i (T )−φ ∗ j i+1 (0) +G j i (m) (z), i = 1,··· ,K− 1 z 1 =φ ∗ j K (T )−φ ∗ j 1 (0) +G j K (m) (z). (3.17) We assume m = 1, i.e. linear approximation, then Eq. (3.17 ) takes the form z i+1 =4 i +H i ·z i z 1 =4 1 +H K ·z K where4 i =φ ∗ j i (T )−φ ∗ j i+1 (0), i = 1,··· ,K− 1 and4 1 =φ ∗ j K (T )−φ ∗ j 1 (0). 39 3.3 Extended Cell Mapping Method Explicitly we have z 2 =4 1 +H 1 z 1 z 3 =4 2 +H 2 z 2 . . . z K =4 K−1 +H K−1 z k−1 z 1 =4 K +H K z K , and substitution yields: z 2 =4 1 +H 1 z 1 z 3 =4 2 +H 2 (4 1 +H 1 z 1 ) =4 2 +H 2 4 1 +H 2 H 1 z 1 z 4 =4 3 +H 3 (4 2 +H 2 4 1 +H 2 H 1 z 1 ) =4 3 +H 3 4 2 +H 3 H 2 4 1 +H 3 H 2 H 1 z 1 z 5 =4 4 +H 4 (4 3 +H 3 4 2 +H 3 H 2 4 1 +H 3 H 2 H 1 z 1 ) =4 4 +H 4 4 3 +H 4 H 3 4 2 +H 4 H 3 H 2 4 1 +H 4 H 3 H 2 H 1 z 1 ... z K =4 K−1 +H K−1 4 K−2 +H K−1 H K−2 4 K−3 +··· +H K−1 ···H 2 4 1 +H K−1 H K−2 ···H 1 z 1 z 1 =4 K +H K z K 40 3.3 Extended Cell Mapping Method Hence z 1 =4 K +H K 4 K−1 +H K H K−1 4 K−2 +H K H K−1 H K−2 4 K−3 +··· +H K H K−1 ···H 2 4 1 +H K H K−1 ···H 1 z 1 or (I−H)z 1 =b =⇒z 1 = (I−H) −1 b (3.18) where H = Q K i=1 H i , b =4 K + P K−1 i=1 ( Q K K=i+1 H K )·4 i . As in the case of a P− 1 solution (see Eq. (3.14)) if there is no bifurcation Eq. (3.18) has a unique solution. 41 4 Analysis of Pitch Motion of a Gravity-Gradient Satellite Stability 4.1 Equation of Motion for Pitch Motion Consider the governing equations of attitude motion for a satellite given in Eq. (2.15). If the roll and yaw (θ 1 and θ 3 ) are initially quiescent, the pitch motion is decoupled from the other angles, yielding I 2 d 2 θ 2 dt 2 + 3 μ R c (t) 3 ! (I 1 −I 3 ) sinθ 2 cosθ 2 =I 2 I ˙ ω O . (4.1) Beyond this point, the subscript is omitted for the pitch angle. By Kepler’s laws of planetary motion, an elliptical orbit is described in polar coor- dinates using the variables (R c ,ν) as (see Fig.2.1) R c (t) = p 1 +e cosν(t) . (4.2) The dependence of ν(t) on time satisfies the law of motion of satellites, and the rate 42 4.1 Equation of Motion for Pitch Motion of change of ν in time is (see Bate et al. [43]) dν dt = √ μp p 2 (1 +e cosν) 2 . (4.3) Figure 4.1: Geometry of gravity gradient satellite motion It is convenient to replace the time t by true anomaly ν as an independent variable so that we investigate periodic solutions in terms of orbital period. Substituting Eq. (4.2) to Eq. (4.1) and using the chain rule, it is explicitly written in terms of ν, as follows: (1 +e cosν) d 2 θ dν 2 − 2e sinν dθ dν + 3k 2 sinθ cosθ =−2e sinν, (4.4) where (θ,ν) are respectively pitch angle and true anomaly angle and k 2 is the ratio of principal moments of inertia defined as k 2 , I 1 −I 3 I 2 . 43 4.2 Periodic Solutions of Pitch Oscillations 4.2 Periodic Solutions of Pitch Oscillations 4.2.1 Periodic Solutions using Cell Mapping In order to analyze the global dynamic behavior of the system in Eq. (4.4), cell mapping method is utilized. We define the states x 1 = θ, x 2 = θ 0 , and the system is expressed as x 0 1 =x 2 x 0 2 = 2e sinν 1 +e cosν x 2 − 3k 2 1 +e cosν sinx 1 cosx 1 − 2e sinν 1 +e cosν . (4.5) Note that the prime denotes derivative with respect toν, and (e,k 2 ) are system param- eters. Then, the cell state space discretization and an unraveling algorithm were imple- mented using MATLAB. The region of interest is a box within± π 2 rad and± π 2 rad/rad for θ and θ 0 , respectively, for considering practical motions. An uniform cell size is chosen for each value of the parameters (e,k 2 ) to adjust the accuracy of the solutions. Eq. (4.5) was integrated over one orbit period T = 2π with initial conditions at the center point of z p to find their image cells. The global properties of the cells were determined using the unraveling algorithm (see Hsu [8]). In previous studies, Beletskii [55], Modi and Brereton [53], and Zlatoustov and Markeev [49] found that there exist three families of periodic solutions depending on e and k 2 . For small eccentricity e = 0.01, k 2 = 1.12 3 is a resonance inertia ratio: if k 2 < 1.12 3 , there is only one stable periodic solution, else if k 2 > 1.12 3 , there are three solutions of which two of are stable and one is unstable. 44 4.2 Periodic Solutions of Pitch Oscillations Analysis for e=0.01 k 2 =0.1 We first study the case of e = 0.01 and k 2 = 0.1 to verify the families of periodic solutions computed in [53]. According to Modi and Brereton [53], one stable P− 1 solution is expected for k 2 < 1.12. ... … Nc= 175,00 0 ... ... 351 352 … 1 2 3 … 350 x 1 x 2 (L) x =-1.2 1 (U) x =1.2 1 (L) x =-1.5 2 (U) x =1.5 2 h =0.0069 1 h =0.0060 2 sink cell sink cell Figure 4.2: Cell state space with e = 0.01, k 2 = 0.1 A rectangular area was considered for the cell state space. As shown in Fig.4.2 in thex 1 -direction 350 cells are used coveringθ from−1.2 rad to 1.2 rad with an interval size h 1 = 0.0069. In the x 2 -direction 500 cells are used covering θ 0 from−1.5 to 1.5 rad/radwithanintervalsizeh 2 = 0.0060. ThusthetotalnumberofcellsN c is 175, 000. The result of global analysis of the cell mapping is presented in Fig.4.3. The white region in Fig.4.3 shows the cells at that are mapped eventually outside of the rectangular region, i.e. sink cells. The sign (+) at the center of circles, (A), is 45 4.2 Periodic Solutions of Pitch Oscillations ν (rad) dθ dν s₁ s₂ s₃ u₁ u₂ u₃ A B₁, B₂ C₁, C₂ Figure 4.3: Global behavior in a cell state space for e = 0.01,k 2 = 0.1 46 4.2 Periodic Solutions of Pitch Oscillations a P− 1 cell which agrees with the results of Beletskii [55], Modi and Brereton [53], and Zlatoustove and Markeev [49]. In addition to theP− 1 solution,P− 2 andP− 3 solutions are discovered: 1) stable P− 2 cells, B 1 and B 2 2) unstable P− 2 cells, C 1 and C 2 3) stable P− 3 cells S 1 ,S 2 , and S 3 4) unstable P− 3 cells U 1 ,U 2 , and U 3 (see Fig.4.3). Moreover, closed regions in Fig.4.3 surrounding the stable fundamental periodic cells show the motion around those cells is bounded. 47 4.2 Periodic Solutions of Pitch Oscillations -0.0615 -0.0339 -0.0063 0.0213 0.0489 0.0765 -0.09 -0.06 -0.03 0 0.03 0.06 0.09 ν (rad) dθ dν -1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5 P-1 Soluon ν (rad) dθ dν A (Close-up) Figure 4.4: Cell state space shownP− 1 solutions fore = 0.01,k 2 = 0.1. A group of P− 1 cells from the cell mapping converge to A = (−1.95e− 7, 0.0283370) As can be observed in Fig.4.3, P− 1 solution computed using the cell mapping consists of a group of cells as shown in Fig.4.4. The enlarged region of the group of cells A is shown in the lower part of the figure. The computation using cell mapping yields accuracy of the order of 2×10 −2 . In order to overcome the accuracy limitations, 48 4.2 Periodic Solutions of Pitch Oscillations we modified the cell mapping technique to refine the periodic solutions and find exact solutions. This technique is referred to as Extended cell mapping (see Section 3.3). Utilizing the extended cell mapping algorithm, P− 1 solution converged to A : (x ∗ 1 ,x ∗ 2 ) = (−1.95e− 07, 0.0283370) denoted red cross in Fig.4.4. The groups of P− 2 cells are shown in Fig.4.5 and Fig.4.6. One of the solutions is stable and the other one is unstable. As before, we applied the extended cell mapping algorithm and found exact P− 2 solutions as shown in the figures. 49 4.2 Periodic Solutions of Pitch Oscillations -0.638 -0.603 -0.569 -0.5342 -0.09 -0.06 -0.03 0 0.03 0.06 -0.09 -0.06 -0.03 0.03 0.06 0.508 0.5353 0.5698 0.6043 0.6388 0 P-2 Soluon dθ dν ν (rad) - 1 . 5 - 1 - 0 . 5 0 0 . 5 1 1 . 5 - 1 . 5 - 1 - 0 . 5 0 0 . 5 1 1 . 5 ν (rad) ν (rad) C₁ C₂ B₁ B₂ B₁ B₂ (Close-up) (Close-up) Figure 4.5: Cell state space shown stable P − 2 solutions for e = 0.01, k 2 = 0.1. A group of P − 2 cells from cell mapping converge to B i = (−0.6131613, 0.0239830), (0.6131484, 0.0239292) 50 4.2 Periodic Solutions of Pitch Oscillations -0.0615 -0.027 0.0144 0.0558 0.21 0.24 0.27 0.3 0.33 0.36 -0.0615 -0.027 0.0144 0.0558 -0.36 -0.33 -0.3 -0.27 -0.24 -0.21 1.5 -1.5 -1 -0.5 0 0.5 1 -1.5 -1 -0.5 0 0.5 1 1.5 P-2 Soluon d θ d ν ν (rad) ν (rad) d θ d ν d θ d ν C₁ C₂ C₁ C₂ B₁ B₂ (Close-up) (Close-up) (Close-up) Figure 4.6: Cell state space shown unstable P − 2 solutions for e = 0.01, k 2 = 0.1. A group of P − 2 cells from cell mapping converge to C i = (0.0001175, 0.3147778), (0.0008436,−0.2638848) TheglobalanalysisofcellmappingallowsustodiscovertwogroupsofP−3solutions 51 4.2 Periodic Solutions of Pitch Oscillations that were not found before. They are shown in Fig.4.7. The stable P− 3 group is denoted by S i and the unstable P− 3 group is denoted by u i where i = 1, 2, 3. These solutions are mapped back: S 1 →S 2 →S 3 →S 1 and U 1 →U 2 →U 3 →U 1 as shown on their trajectories in Fig.4.7. ν (rad) dθ dν - 1 . 5 - 1 - 0 . 5 0 0 . 5 1 1 . 5 - 0 . 8 - 0 . 6 - 0 . 4 - 0 . 2 0 0 . 2 0 . 4 0 . 6 S ₁ S ₂ S ₃ u ₁ u ₃ u ₂ Figure 4.7: Cell state space shownP− 3 solutions and trajectories in phase plane for e = 0.01, k 2 = 0.1 (S i : stable solutions, U i : unstable solutions) Again, we applied the extended cell mapping algorithm to find the exact solutions. The enlarged areas aroundP− 3 solutions in the cell state space are shown in Fig.4.8 and Fig.4.9. 52 4.2 Periodic Solutions of Pitch Oscillations -0.2 -0.1581 -0.117 -0.0753 -0.0339 0.0075 0.0489 0.3 0.33 0.36 0.39 0.42 0.45 0.48 0.51 0.54 0.57 0.6 0.63 0.66 P-3 Solution -1.1724 -1.131 -1.0896 -1.0482 -1.0068 -0.9654 -0.924 -0.27 -0.24 -0.21 -0.18 -0.15 -0.12 -0.09 -0.06 -0.03 0 0.03 0.06 P-3 Solution 0.9252 0.9666 1.008 1.409 1.0908 1.1322 1.1736 -0.27 -0.24 -0.21 -0.18 -0.15 -0.12 -0.09 -0.06 -0.03 0 0.03 0.06 P-3 Solution ν (rad) ν (rad) ν (rad) dθ dν dθ dν S ₁ S ₃ S ₂ Figure 4.8: Cell state space shown stable P− 3 solutions for e = 0.01, k 2 = 0.1 S i : (−0.0006549, 0.5332902), (1.1431803,−0.1327799), (−1.1426959,−0.1331036) 53 4.2 Periodic Solutions of Pitch Oscillations -0.57 -0.54 -0.51 -0.48 -0.45 -0.42 P-3 Solution -0.02 0.0075 0.042 0.076 0.111 -1.2 -1.1724 -1.131 -1.0896 -1.055 0.06 0.09 0.12 0.15 0.18 0.21 P-3 Solution 1.056 1.0908 1.1322 1.1736 1.2 0.06 0.09 0.12 0.15 0.18 0.21 P-3 Solution ν (rad) dθ dν ν (rad) ν (rad) dθ dν u 1 u ₂ u ₃ Figure 4.9: Cell state space shown unstable P− 3 solutions for e = 0.01, k 2 = 0.1 U i : (1.1323459, 0.1684754), (−1.1339993, 0.1675323), (0.0009308,−0.4855702) For the further study of the neighborhood of the fundamental periodic solutions, fast Fourier analysis was applied. As shown in Fig.4.10, the fundamental periodic solutions (A), (B), and (C) have one large peak atP−1,P−2, andP−3, respectively. However, (D) in Fig.4.10 which is tested at (0.2, 0) located in an invariant surface, close toP−1 solution but not a point of fundamental periodic solution shows a combination of three irrational frequencies. 54 4.3 Stability Analysis X: 0.3329 Y: 1.219 X: 1 Y: 0.08696 (C) P-3 Soluon (Period 1/0.333) (B) P-2 Soluon (Period 1/0.5) (A) P-1 Soluon (Period 1 ) Frequency * (2 π ) (D) Neighborhood of fundamental periodic soluon (0.7,0) is bounded and shows combinaon of irraonal frequencies 0 0 . 5 1 1 . 5 2 2 . 5 3 3 . 5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Y( ) ν Frequency * (2 π ) 0 0 . 5 1 1 . 5 2 2 . 5 3 3 . 5 X: 0.5031 Y: 0.3692 X: 0.4725 Y: 0.2402 X: 0.5277 Y: 0.09292 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Y( ) ν Frequency * (2 π ) 0 0 . 5 1 1 . 5 2 2 . 5 3 3 . 5 0 0.005 0.01 0.015 0.02 0.025 Y( ) ν Frequency * (2 π ) 0 0 . 5 1 1 . 5 2 2 . 5 3 3 . 5 0 0.1 0.2 0.3 0.4 Y( ) ν 0.5 0.6 0.7 Figure 4.10: Fast Fourier analysis results fore = 0.01,k 2 = 0.1 at three fundamental periodic solutions and at a neighborhood point (0.2, 0) 4.3 Stability Analysis To determine the stability of the periodic solutions, matrix H is computed using the point mapping algorithm described in Eq. (3.6). Let δ = 2θ to use the double angle formula, then Eq. (4.4) can be written as 55 4.3 Stability Analysis (1 +e cosν)δ 00 − 2e sinνδ 0 + 3k 2 sinδ =−4e sinν (4.6) where the prime denotes derivative with respect to true anomalyν. Defining the states x 1 =δ, x 2 =δ 0 , we get x 0 1 =x 2 x 0 2 = 1 (1+e cosν) (2e sinνx 2 − 3k 2 sinx 1 − 4e sinν) . (4.7) The discrete counterpart to Eq. (4.7) is then formed by the point mapping algorithm in Section 3.1, yielding x 1 (n + 1) =G 1 (x(n);e,k 2 ) x 2 (n + 1) =G 2 (x(n);e,k 2 ) . (4.8) To compute the approximate point mapping, Eq. (4.7) was integrated from τ = 0 to τ = 2π using time step size h = 2π/100 employing the coefficients of fourth order Runge-Kutta method. The expansion in Eq. 4.9 was carried out to include terms up to fifth power in the states. The mapping is written as x 1 (n + 1) = 5 X i,j=1 φ ij x i 1 (n)x j 2 (n) (4.9) x 2 (n + 1) = 5 X i,j=1 ψ ij x i 1 (n)x j 2 (n) and computed coefficients of the point mappingG for the case ofe = 0.01 andk 2 = 0.1 are shown in Table 4.1. 56 4.3 Stability Analysis i j φ ij ψ ij i j φ ij ψ ij 0 0 0.02808648 0.11079739 4 0 0.00453579 -0.00089991 1 0 -0.95984591 0.17149894 3 1 0.00308941 0.01809886 0 1 -0.48822485 -0.95460108 2 2 -0.06515537 -0.03487541 2 0 -0.02964038 -0.01137257 1 3 0.01423299 0.04633806 1 1 0.08001880 0.03309032 0 4 -0.17798911 -0.04779005 0 2 -0.19857046 0.00410416 5 0 0.02602411 0.00728993 3 0 -0.06651000 -0.13025951 4 1 0.03746501 0.04611338 2 1 0.36846353 -0.06631159 3 2 0.09570022 -0.06629808 1 2 -0.23487829 -0.36317193 2 3 0.4456393 0.17871561 0 3 1.22059444 -0.19700666 1 4 0.09751830 -0.18759918 0 5 0.72837480 0.14298076 Table 4.1: Computed coefficients of the point mapping G for e = 0.01 and k 2 = 0.1 In order to obtain the P− 1 solutions, x ∗ is determined which satisfies x ∗ −G(x ∗ ;e,k 2 ) = 0 (4.10) asx ∗ = (6.92e−10, 0.028337) which matches with the point (A) from the cell mapping. Furthermore,P− 2 andP− 3 solutions can be obtained byx ∗ −G 2 (x ∗ ;e,k 2 ) = 0 and x ∗ −G 3 (x ∗ ;e,k 2 ) = 0, respectively. Consider again the Eq. (4.7) for stability analysis of x ∗ . Stability analysis is based on the consideration of the spectral properties of Jacobian matrix H (see Eq. (3.6)) of the point mapping G. Let us define ξ(ν) =x(ν)−x ∗ as the perturbation of the state x about the solution x ∗ . The motion of ξ close to x ∗ is governed by 57 4.3 Stability Analysis ˙ ξ =f(ν,x ∗ (ν) +ξ(ν),s 0 )−f(ν,x ∗ ,s 0 ) = ∂f ∂x x ∗ ξ(ν) +o(ξ 2 ). Thus, the linear part of Eq. (4.7) about the P− 1 solution is ˙ ξ 1 =ξ 2 ˙ ξ 2 = −3k 2 cosx ∗ 1 (ν) 1 +e sinν ξ 1 + 2e sinν 1 +e cosν ξ 2 , (4.11) and the point mapping formulation of the Eq. (4.11) is ξ m+1 =H(s)ξ m where H(s)∈R N×N is given by H(s) =H K (s)H K−1 (s)···H 1 (s), and H j = " ∂G(x,s) ∂x # x=x ∗(j) . Since the system is Hamiltonian in this case, the periodic solution is linearly stable when the eigenvalues ofH lie on the unit circle in the complex plane and are pairwise different. Otherwise, the solution is unstable. 58 4.4 Bifurcations x ∗ H det H eig H P-1 (stable) A: [-1.95e-07, 0.0283] −0.9561 −0.4991 0.1722 −0.9561 1.0000 −0.9561+0.2932i −0.9561−0.2932i P-2 (stable) B 1 : [-0.6132, 0.0239] 0.9077 −0.2354 0.5735 0.9530 1.0000 0.9303+0.3667i 0.9303−0.3667i P-2 (unstable) c 1 : [0.0001, 0.3148] 1.0573 −2.6049 −0.0455 1.0578 1.0000 0.7134 1.4017 P-3 (stable) S 1 : [-0.0006, 0.5333] −0.4073 −20.7667 0.0374 −0.5485 1.0000 −0.4779+0.8784i −0.4779−0.8784i P-3 (unstable) u 1 : [1.1323, 0.1685] −4.2665 −8.1721 5.0553 9.4486 1.0000 0.2007 4.9814 Table 4.2: Periodic solutions’ stability check by computing eigenvalues of matrix H Table 4.2 shows the stability properties of the periodic solutions computed using the point mapping. 4.4 Bifurcations 4.4.1 P-1 to P-1 bifurcation Period 2π solutions computed using the cell mapping method for three different ec- centricities are shown in Fig.4.11. The results shown in the figure verify the results of Beletskii [55], Zlatoustov et al. [49] for the three families of P− 1 solutions. 59 4.4 Bifurcations 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 0 . 7 0 . 8 - 1 . 5 - 1 - 0 . 5 0 0 . 5 1 1 . 5 k 2 A m p l i t u d e e = 0 . 1 e = 0 . 0 5 e = 0 . 0 1 e = 0 . 0 1 e = 0 . 0 5 e = 0 . 1 Figure 4.11: P− 1 toP− 1 bifurcation for three different ‘e’s by changingk 2 . (Solid line: stable solution, Dashed line: unstable solution) A bifurcation diagram Fig.4.11 is presented with bifurcation parameter k 2 vs. the amplitudeofthepitchmotionwiththreedifferentorbitaleccentricities. Atsomepoints of k 2 a stable P− 1 solution bifurcates into three P− 1 solutions. Two of them are stable and one is unstable. At the largere, the bifurcation occurs at the largerk 2 . The bifurcation is emerging at k 2 ≈ 1.12 3 for e = 0.01, k 2 ≈ 1.36 3 for e = 0.05, and k 2 ≈ 1.61 3 for e = 0.1. The original branch of stable P− 1 solutions are the curves of the left upper side, and bifurcated stableP− 1 and unstableP− 1 solutions are the curves of the right bottom side. Solid lines present stable P− 1 solutions and dashed lines, the counterparts of the parabolas present unstable P− 1 solutions.. The numerical solutions that were investigated by Modi et al. [53] and Beletskii [55] are repeated here to verify the accuracy of our approach. They demonstrated the 60 4.4 Bifurcations solutions using boundary condition problem setting as θ(0) = 0, θ(2π) = 0 and the repeated results are shown in Fig.4.12 for the cases k 2 ≥ 0.4. It clearly shows where a P− 1 solution bifurcates into three P− 1 solutions. Note that the x−axis is k 2 in Fig.4.11 and e in Fig.4.12, but the solutions match with each other. 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 0 . 7 0 . 8 0 . 9 - 2 - 1 . 5 - 1 - 0 . 5 0 0 . 5 1 1 . 5 2 e k = 1 2 k = 0 . 8 2 k = 0 . 6 2 k = 0 . 5 2 k = 0 . 4 2 A m p l i t u d e Figure 4.12: Numerical P− 1 solutions for different e and k 2 (Does not show the stability of solutions) Analysis for e = 0.01 To further understand the dynamics around the bifurcation point, the invariant surface was studied. The 2π-period solutions and invariant surfaces surrounding the solutions for the case of e = 0.01 are shown in Fig.4.13. 61 4.4 Bifurcations - 1 . 5 - 1 - 0 . 5 0 0 . 5 1 1 . 5 - 1 . 5 - 1 - 0 . 5 0 0 . 5 1 1 . 5 - 1 . 5 - 1 - 0 . 5 0 0 . 5 1 1 . 5 - 1 . 5 - 1 - 0 . 5 0 0 . 5 1 1 . 5 - 1 . 5 - 1 - 0 . 5 0 0 . 5 1 1 . 5 - 1 . 5 - 1 - 0 . 5 0 0 . 5 1 1 . 5 2 2 2 -1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5 e=0.01, k=0.3 x ₁ x ₂ e=0.01, k=0.35 e=0.01, k=0.37 -1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5 e=0.01, k=0.38 2 e=0.01, k=0.4 2 e=0.01, k=0.6 2 stable P-1 solution unstable P-1 solution x ₁ x ₂ x ₁ x ₂ x ₁ x ₂ x ₁ x ₁ a ) b ) c ) d ) e ) f ) Figure 4.13: Invariant surface of phase plane representation for e = 0.01. One stable P− 1 solution bifurcates into two stable P− 1 solutions and an unstable P− 1 solution. In Fig.4.13a, for k 2 = 0.3, there exists an uniform invariant surface surrounding a stable P− 1 solution. However, for k 2 = 0.35, it can be observed that the invariant surface moves and deforms as shown in Fig.4.13b. The characteristics of the surface is 62 4.4 Bifurcations changing across the inner curve. While the inner surface is deforming, another stable P− 1 solution appeared under the deformed region in Fig.4.13c. The evolution of stable P− 1 solutions is shown in Fig.4.13c through f. As shown in Fig.4.13e and f, the invariant region around the new P− 1 solutions grows and pushes the surface surrounding the upper stable solution away as k 2 is getting larger. Another feature to note is that an unstableP− 1 solution is emerging atk 2 = 0.38, Fig.4.13d, marked as an empty circle. 4.4.2 Bifurcation P-1 solutions to P-2 solutions We have observed that different period solutions, P−K solutions, exist for different system parameters (e,k 2 ). It is important to understand when bifurcation occurs since these bifurcations may lead to complex dynamic behavior, e.g. chaotic behavior. Analysis for k 2 = 0.1 We consider the case for whichk 2 = 0.1 ande varies between 0.01 and 0.7. In the case of e = 0.01, a stable P− 1 solution was found as discussed in Fig.4.4. We also found stable and unstable P− 2 solutions (see Fig.4.5, Fig.4.6). Here, we are interested in the mechanism that leads to the creation of P− 2 solutions. Consider the unstable P− 2 branch in Fig.4.14 first. As e decreases from 0.7 to 0.01, θ 0 of P− 1 solution also decreases as shown in Fig.4.14. 63 4.4 Bifurcations 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 0 . 7 - 0 . 4 - 0 . 2 0 0 . 2 0 . 4 0 . 6 0 . 8 1 e b f c 1 b f c 2 x ₂ k = 0 . 1 2 c 1 c 2 Figure4.14: Bifurcation of unstableP−1 to stableP−1 and unstableP−2 solutions for k 2 = 0.1,e = 0.01∼ 0.7 at bifurcation point bfc1 (for bfc2 : see Fig.4.15). The decreasing unstableP− 1 solution becomes stable and unstableP− 2 solutions emerge at e = 0.1 which is bfc1 in Fig.4.14. Next, the stableP− 2 branch in Fig.4.15 is considered. We only consideredx 2 -axis so far, but the bifurcation for this case is interpreted in x 1 -axis. 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 0 . 7 - 0 . 8 - 0 . 6 - 0 . 4 - 0 . 2 0 0 . 2 0 . 4 0 . 6 0 . 8 e b f c 1 b f c 2 x ₁ k = 0 . 1 2 B 2 B 1 Figure 4.15: Bifurcation of stable P− 1 loses stability and stable P− 2 solutions emerge fork 2 = 0.1,e = 0.01∼ 0.7 at bifurcation pointbfc2 (forbfc1: see Fig.4.14) As e decreases from 0.7 to 0.01, the stable P− 1 solution loses stability and stable 64 4.4 Bifurcations P−2solutionsemergeate = 0.37(bfc2inFig.4.15). Thetwoindependentbifurcations are shown in 3−D plot as Fig.4.16. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 -1 -0.5 0 0.5 1 -0.5 0 0.5 1 e x ₁ x ₂ b f c1 b f c2 stable P-1 solution unstable P-1 solution stable P-2 solution unstable P-2 solution Figure 4.16: Bifurcations for k 2 = 0.1: combined Fig.4.14 and Fig.4.15 It is also important to see the invariant surface surrounding the periodic solutions. The phase plane, Fig.4.17 is obtained using cell mapping like Fig.4.2. Around the fundamental solutions, curves show invariant surfaces which contain quasi-periodic motions. For the case of k 2 = 0.1, as e gets smaller the bounded region increases. Where e > 0.1, there are no more bounded regions around the fundamental periodic solutions. 65 4.4 Bifurcations -1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5 x ₁ e=0.04, k =0.1 2 x ₂ -1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5 e=0.01, k =0.1 2 x ₁ x ₂ -1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5 e=0.07, k =0.1 2 x ₁ x ₂ -1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5 e=0.1, k =0.1 2 x ₁ x ₂ stable P-1 solution unstable P-1 solution stable P-2 solution unstable P-2 solution stable P-3 solution unstable P-3 solution Figure 4.17: Invariant surface of phase plane representation for k 2 = 0.1 4.4.3 P-1 to P-3 Bifurcation and P-3 to P-9 Bifurcation Analysis for e = 0.1 For a case of e = 0.1, it was found that P− 1 solutions as shown in Fig.4.11. One of stable P− 1 solutions bifurcates into two P− 3 solutions at k 2 ≈ 0.16: one is stable, g i , and other one is unstable, h i , where i = 1,..., 3. Bifurcation diagrams are given in Fig.4.18 showing a 3-Dimensional representation. Fig.4.19 and Fig.4.20 show its 66 4.4 Bifurcations projections on x 1 and x 2 vs. k 2 , respectively. -1 -0.5 0 0.5 1 0.16 0.18 0.2 0.22 -0.5 0 0.5 1 k 2 x 1 x 2 stable P-3 unstable P-3 stable P-1 g 2 g 1 g 3 h 2 h 1 h 3 k 1,2,3 k 4,5,6 k 7,8,9 unstable P-9 Figure 4.18: P− 1 to P− 3 bifurcation and P− 3 to P− 9 bifurcation for e = 0.1 in 3-D -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 0.15 0.16 0.17 0.18 0.19 0.2 0.21 0.22 k 2 x 1 stable P-3 unstable P-3 stable P-1 g 2 g 1 g 3 h 2 h 1 h 3 k 1,2,3 k 4,5,6 k 7,8,9 , unstable P-9 Figure 4.19: P− 1 to P− 3 bifurcation and P− 3 to P− 9 bifurcation for e = 0.1 from the view of x 1 -axis 67 4.4 Bifurcations 0.15 0.16 0.17 0.18 0.19 0.2 0.21 0.22 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 k 2 x 2 stable P-3 unstable P-3 stable P-1 h 3 k 1,2,3 k 4,5,6 k 7,8,9 g 1 g 2 g 3 h 1 h 2 unstable P-9 Figure 4.20: P− 1 to P− 3 bifurcation and P− 3 to P− 9 bifurcation for e = 0.1 from the view of x 2 -axis The continuous trajectories associated with theP− 3 solutions show the periodicity of the solution explicitly. As an example, a stable trajectory associated with g i , i = 1, 2, 3 is examined for the case ofk 2 = 0.19 (see Fig.4.21). Integration of one period of time for the initial point g 1 results g 2 , integration of one period of time for g 2 results g 3 , and the same step forg 3 tog 1 , and the mapping repeats as shown in Fig.4.21. For comparison of the trajectories,P−1 solution is also shown in Fig.4.21 with a ’+’ sign. The continuous response as true anomaly progresses is given in Fig.4.22. The initial point g 1 is integrated for t = 54π, and the response shows a period three motion. 68 4.4 Bifurcations g 1 g 2 g 3 - 0 . 8 - 0 . 6 - 0 . 4 - 0 . 2 0 0 . 2 0 . 4 0 . 6 dθ dν - 0 . 8 - 0 . 6 - 0 . 4 - 0 . 2 0 0 . 2 0 . 4 0 . 6 0 . 8 ν (rad) Figure4.21: Continuous trajectory for stableP−3 solutionsg i ,i = 1, 2, 3 ate = 0.1, k 2 = 0.19 and P− 1 solution ’+’ θ ν ( rev ) 0 5 1 0 1 5 2 0 2 5 3 0 - 0 . 8 - 0 . 6 - 0 . 4 - 0 . 2 0 0 . 2 0 . 4 0 . 6 0 . 8 - 0 . 8 - 0 . 6 - 0 . 4 - 0 . 2 0 5 1 0 1 5 2 0 2 5 3 0 0 0 . 2 0 . 4 0 . 6 ν ( rev ) dθ dν Figure 4.22: True anomaly vs. angular responses and angular rate responses for g 1 showing the period of three motion The unstableP− 3 solutions,h i ,i = 1, 2, 3 and a continuous trajectory atk 2 = 0.19 are shown in Fig.4.23. It is clear to see h 1 →h 2 →h 3 →h 1 in 3T. However, the true anomaly response diverges over time. 69 4.4 Bifurcations - 0 . 8 - 0 . 6 - 0 . 4 - 0 . 2 0 0 . 2 0 . 4 0 . 6 0 . 8 - 0 . 8 - 0 . 6 - 0 . 4 - 0 . 2 0 0 . 2 0 . 4 0 . 6 h 2 h 1 h 3 dθ dν ν (rad) Figure 4.23: Continuous trajectory for unstable P− 3 solutions h i , i = 1, 2, 3 at e = 0.1, k 2 = 0.19 The new stable P− 3 solution bifurcates into unstable P− 9 solutions (k i , i = 1,··· , 9) at k 2 ≈ 0.20 as shown in Fig.4.18. The continuous trajectories associated with the P− 9 solutions are presented in Fig.4.24. As can be observed, the initial point denoted as ‘0’ is integrated for t = 9T, and is mapped back to the initial state. The true anomaly response shows that angular and angular rate motions repeat with 9T, but start to be perturbed and diverging for its second orbit as shown in Fig.4.25 and Fig.4.26. 70 4.4 Bifurcations - 1 . 5 - 1 - 0 . 5 0 0 . 5 1 1 . 5 - 1 - 0 . 8 - 0 . 6 - 0 . 4 - 0 . 2 0 0 . 2 0 . 4 0 . 6 0 . 8 2 7 4 5 6 3 1 8 9 0 dθ dν ν (rad) Figure 4.24: Continuous trajectory for P− 9 solutions at e = 0.1, k 2 = 0.21 0 3 6 9 12 15 18 21 24 27 30 -1.5 -1 -0.5 0 0.5 1 1.5 True anomaly (rev) θ Figure 4.25: True anomaly vs. angular response for a P− 9 solution at e = 0.1, k 2 = 0.21 71 4.4 Bifurcations 0 3 6 9 12 15 18 21 24 27 30 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 True anomaly (rev) d θ d ν Figure 4.26: True anomaly vs. angular rate response for a P− 9 solution ate = 0.1, k 2 = 0.21 The fast Fourier analysis response verifies the period of the solution. As shown in Fig.4.27, a combination of fractions of period 9T is found for one of theP−9 solutions. Figure 4.27: Fast Fourier analysis results for a P− 9 solution at e = 0.1, k 2 = 0.21 The variations of invariant surfaces according to the bifurcations are shown in Fig.4.28. 72 4.4 Bifurcations Figure 4.28: Invariant surface of phase plane representation, e = 0.1, k 2 = 0.21 When new P− 3 solutions emerged at k 2 = 0.175, the associated deformation of invariant surfaces can be seen. Following the track of stable P− 3 solutions as k 2 is gradually increasing, they move apart and the size of surrounding invariant surfaces is decreasing. Finally, in between k 2 = 0.2 and k 2 = 0.21, the invariant surfaces around the stableP− 3 solutions disappear and thenP− 9 solutions are born. The invariant surface of the stableP− 1 solution, however, is growing in size ask 2 is increasing from k 2 = 0.175 to k 2 = 0.21. 73 4.5 Summary 4.5 Summary Gravity gradient torque stabilized satellites have been studied for their practical val- ues. The motion is described in three coupled nonlinear non-autonomous second order differential equations which are difficult to solve analytically. The motion of the pitch axis could be decoupled from the other axes motion which simplifies the dynamics system. Moreover, a study of the pitch motion is important for an Earth pointing stabilized satellite. In this study, the pitch motion of a gravity gradient satellite was investigated using a new approach which combines numerical and analytical analysis. The study was performed globally, and it was able to demonstrate the previous results and find new families of periodic solutions. The stability of the solutions was examined analytically using a point mapping technique and fast Fourier analysis. The invariant surfaces which give ideas on how the periodic solutions emerge and behave were found by using cell mapping analysis. Moreover, P− 1 to P− 2, P− 1 to P− 3, and P− 3 to P− 9 bifurcations were introduced. Finding bifurcation diagrams is important for study of complex dynamic behavior such as chaos. 74 5 Analysis of Spinning Satellite 5.1 Problem Formulation Consider a spinning satellite with the principal body coordinate frame, ¯ e B , e B 1 coin- ciding with the axis of symmetry as shown in Fig.5.1. Note that the coordinates are different from the gravity gradient pitch motion analysis, presented in Chapter 4. Let ¯ e I represent the inertial frame with e I 1 is normal to the orbital plane, e I 2 −e I 3 . The orbital reference frame, ¯ e O , is moving with true anomalyν whilee O 2 is aligned with the radius vector, e O 3 is the velocity direction, and e O 1 completes the orthogonal set. The attitude of the satellite is defined by Euler angles, γ, β, and α (3-2-1) relative to the orbital reference frame. Figure 5.1: Spinning satellite description 75 5.1 Problem Formulation Transformation from Orbital Reference Frame to Body Frame The first rotation,γ, aboute O 3 axis is referred to as roll, the second rotation, β, about the second axis is referred to as yaw, and the last rotation α about the symmetry axis is the spin. The sequence of the rotations from the orbital reference frame to the body frame is written using arrows with the reference axis rotating about and rotation angles as follows: [¯ e O ] e O 3 ,θ 3 −−−→ [¯ e a ] e a 2 ,θ 2 −−−→ [¯ e b ] e b 1 ,θ 1 −−−→ [¯ e B ]. The corresponding rotation matrices are B C b =         1 0 0 0 c 1 s 1 0 −s 1 c 1         , b C a =         c 2 0 −s 2 0 1 0 s 2 0 c 2         , a C O =         c 3 s 3 0 −s 3 c 3 0 0 0 1         Then the angular velocity of the system with respect to inertial frame is given by I ¯ ω B = I ¯ ω O + O ¯ ω B where O ¯ ω B = [¯ e a ]         0 0 ˙ θ 3         + [¯ e b ]         0 ˙ θ 2 0         = [¯ e B ]         −s 2 ˙ θ 3 ˙ θ 2 c 2 ˙ θ 3         76 5.2 Equation of motion for spinning satellite 5.2 Equation of motion for spinning satellite The independent variable is converted from time to true anomaly by the relation dν dt = √ μp p 2 (1 +e cosν) 2 d 2 ν dt 2 = −2e sinν 1 +e cosν dν dt ! 2 . Then the equations of motion for an axi-symmetric rigid spinning satellite with gravity gradient torque are: β 00 −γ 0 cosγ +k 1 (σ + 1) 1 +e 1 +e cosν 2 (γ 0 cosβ + sinβ cosγ) −(γ 0 cosβ + sinβ cosγ)(cosβ cosγ−γ 0 sinβ) − 3(k 1 − 1) 1 +e cosν ! sinβ cosβ sin 2 γ− 2e sinν 1 +e cosν (β 0 − sinγ) = 0, (5.1) γ 00 cosβ− 2β 0 (γ 0 sinβ− cosβ cosγ)−k 1 (σ + 1) 1 +e 1 +e cosν 2 (β 0 − sinγ) + 3(k 1 − 1) 1 +e cosν − 1 ! cosβ sinγ cosγ− 2e sinν 1 +e cosν (γ 0 cosβ + sinβ cosγ) = 0, where the superscript “ 0 ” denotes derivative with respect toν, andβ andγ are rotation angles about e 2 , e 3 axis, respectively. The three parameters that characterize the dynamic behavior of the satellite are: 1) orbital eccentricity e, 2) the satellite inertia parameter (k 1 = I 1 I 2 ) whereI is a body’s principal moment of inertia diagonal matrix, and 3) the dimensionless spin parameter σ, σ = ˙ α ˙ ν β=γ=0 77 5.3 Periodic Motion Analysis whereα represents spin angle aboute B 1 axis and the superscript “·” denotes derivative with respect to t. 5.3 Periodic Motion Analysis 5.3.1 Periodic solutions using Cell Mapping In order to analyze the system, the system dynamics described in Eq. (5.1) is written in state space form by defining the states asx 1 =β,x 2 =γ,x 3 =β 0 , andx 4 =γ 0 . The region of interest in the state space was chosen as±π/2 (rad) for angular positions (β and γ) and±2.51 (rad/rad) for angular velocities (β 0 and γ 0 ). The cell mapping is obtained for the four-dimensional dynamical system by integrating over one orbital period T = 2π with initial conditions as the center point of cells, z 0 , to find z T , for given system parameters, satellite’s spin parameter σ, its inertia parameter k 1 , and orbital eccentricitye. This process of searching image cells is repeated for every single cell through the cell state spaceS completing the mapC. Then the effect ofσ,k 1 , and e on the attitude dynamics of the satellite was investigated. The following parameter effects were discussed: • Effect of spin direction (sign(σ)), for|σ|≈ 0.5, k 1 = 2, e = 0.1 • Effect of high spin rate (σ), for 0.5≤σ≤ 2, k 1 = 2, e = 0.1 • Effect of orbital eccentricity (e), for σ = 0.5, k 1 = 2, 0.1≤e≤ 0.2 • Effect of satellite shape (k 1 ), for σ = 0.5, 0<k 1 < 2, e = 0.1 5.3.2 Periodic solutions verification procedure We presented the newly found solutions with the procedure as follows: 78 5.4 Effect of spin direction e = 0.1, k 1 = 2 1. Presenting the P−K solutions in β−β 0 and γ−γ 0 planes. 2. Demonstratingcontinuousbehaviorofasolutionrepresentingafamilyofperiodic solutions by selecting it as an initial condition for numerical integration. 3. Plotting angular position plane trajectories and β response versus true anomaly for the base period showing the periodicity clearly. 4. If necessary, repeating step 3 for longer time responses showing the solution remains confined or diverges. 5. Performing fast Fourier analysis to confirm frequency content of the given solu- tion. 5.4 Effect of spin direction e = 0.1, k 1 = 2 The effect of the satellite’s spin direction on the attitude dynamics is studied for fixed inertia parameter and orbital eccentricity set to bek 1 = 2 ande = 0.1 when|σ|≈ 0.5. 5.4.1 Periodic solutions for positive spin direction (σ≈ 0.5) The analysis was performed for σ within 0.4≤ σ≤ 0.5, i.e., about a half spin about body’s symmetry axis while one orbital revolution is performed. The resulting P− 1 solutions, i.e. period of T = 2π, are shown in Fig.5.2. 79 5.4 Effect of spin direction e = 0.1, k 1 = 2 (a) σ = 0.5 (b) σ = 0.46 (c) σ = 0.44 (d) σ = 0.4 Figure 5.2: β−β 0 (left), γ−γ 0 (right) cell mapping stroboscopic phase plane for P− 1 solutions (k 1 = 2, e = 0.1) The origin is a P− 1 solution in all four cases as it is shown as a point in the cross-section phase planes, since it is an equilibrium point for all values of system 80 5.4 Effect of spin direction e = 0.1, k 1 = 2 parameters. Note that the equilibrium point can be interpreted as a P−K solution for K = 1, 2,.... Two different amplitude solution families with period of 2π exist forming closed curves in the phase plane. A large amplitude solution family (outer curve) is within±0.8rad (±45.8366 o ), and a small amplitude solution family (inner curve) is within±0.3rad (±17.1887 o ), for both β and γ. Amplitudes of both families grow with a decreasing spin parameter from σ = 0.5 to σ = 0.4. To demonstrate the stability characteristics of the solutions, the eigenvalues of the monodromy matrix H, computed from Eq. (3.6), were studied for both families. A representative analysis for each family for σ = 0.4 results that all the eigenvalues are complex conjugates that reside on the unit circle as shown in Fig.5.3, i.e. the solutions are stable. Figure 5.3: Eigenvalues of matrixH in the complex domain for twoP− 1 solutions, σ = 0.4 To verify this analysis result, the continuous dynamics in Eq. (5.1) was simulated with the fixed point solution as an initial condition. Corresponding trajectories in β−β 0 phase plane and β response vs. true anomaly are shown for the case of σ = 0.4, as an example, in Fig.5.4. One of each family of solutions is used as the initial condition for the simulations and shown in Table 5.1. 81 5.4 Effect of spin direction e = 0.1, k 1 = 2 Table 5.1: Samples of P− 1 solution for σ = 0.4 σ K θ 1 (rad) θ 2 (rad) θ 0 1 θ 0 2 index 0.4 1 -0.359278382 -0.675066173 1.752896404 -1.098311492 Large 0.4 1 -0.131996734 -0.239868133 0.377344878 -0.449253523 Small (a) σ = 0.4, Large amplitude solution for 1·T (b) σ = 0.4, Small amplitude solution for 1·T Figure 5.4: Trajectory in β−β 0 phase plane (left) and β response vs. true anomaly (right) of P− 1 solution for 1·T, σ = 0.4 Trajectories from ν = 0 to ν =T are shown in Fig.5.4, and trajectories from ν = 0 toν = 14·T are shown in Fig.5.5 for both families. A small dot in the plots on the left hand side (phase planes) indicates that the initial and final states are the same. The solutions are repeating every period. In addition it can be observed that there exist higher frequency oscillations within the period T for both large and small amplitude solutions. 82 5.4 Effect of spin direction e = 0.1, k 1 = 2 (a) Large amplitude solution for 14·T (b) Small amplitude solution for 14·T Figure 5.5: Trajectory in β−β 0 plane (left) and β vs. true anomaly (right) of the P− 1 solutions for 14·T, σ = 0.4 To further investigate the high frequency content, i.e. frequencies higher than ω = 1/T, FFT analysis of the two solutions was performed. The results are shown in Fig.5.6. (a) Large amplitude solution (b) Small amplitude solution Figure 5.6: Fast Fourier analysis for P− 1 solution for σ = 0.4 83 5.4 Effect of spin direction e = 0.1, k 1 = 2 As can be observed, the largest magnitude is observed at ω = 4/T for the large amplitude solution family, andω = 3/T for the small amplitude solution family. These imply that the dominant periods for the two families are T/4 and T/3, respectively. It can also be observed that the two families contain frequencies with ω = n/T, n = 1, 2,..., 6, i.e. solution period contains T/n. It should be noted that all these solutions appear asP−1 solutions in discrete-time representation of system dynamics, and therefore, are found by the cell mapping approach. 5.4.2 Bifurcations for 0.3<σ< 0.5 As depicted in Fig.5.2, a small amplitude solution appears in the range 0.44 < σ < 0.46. Hence, we investigated the bifurcation of the small amplitude P− 1 solution from the origin, for spin rate range 0.4<σ< 0.5. A bifurcation diagram of the origin as function ofσ is shown in Fig.5.7. Fig.5.7a presents the solutions in theβ−β 0 plane vs. σ and Fig.5.7b shows the evolution of the maximum magnitude of β as function of σ. Figure 5.7: Bifurcation diagram for the origin with positive spin direction 0.3<σ< 0.5 (k 1 = 2, e = 0.1) 84 5.4 Effect of spin direction e = 0.1, k 1 = 2 As indicated in Eq. (3.7), a pair of eigenvalues ofH evaluated at the originalP− 1 solution satisfies λ 1 =λ 2 = 1 when P− 1 to P− 1 bifurcation occurs. Fig.5.8 shows that indeed one pair of eigenvalues at the origin moves along the unit circle, as σ decreases, and reaches Re(λ 1,2 ) = 1 for σ = 0.446. By continuing to decrease σ, the pair of eigenvalues is following the unit circle implying that no changes in stability characteristics by the bifurcation. Figure 5.8: Eigenvalues of H at the origin showing P− 1 to P− 1 bifurcation for 0.43<σ< 0.45 5.4.3 Periodic solutions for negative spin direction (σ≈−0.5) Periodic solutions for a negative spinning parameter in the range of−0.5≤σ≤−0.55 are investigated. The negative sign implies that a satellite spins in a direction opposite tothedirectionoforbitalrevolution. Periodicsolutionsofperiod 2π arefoundasshown in Fig.5.9. The empty circles represent unstable P− 1 solutions, and the disks are stable equilibrium points, the origin. It can be observed that the origin is a P− 1 solution in all cases, as it is an equilibrium point. A paired large amplitude P− 1 solutions is found, and their amplitude increases as the magnitude ofσ increases. The cross signs for σ = −0.5 and σ = −0.52 are P− 2 solutions and triangles when σ =−0.55 are P− 3 solutions. 85 5.4 Effect of spin direction e = 0.1, k 1 = 2 (a) σ =−0.5 (b) σ =−0.52 (c) σ =−0.55 Figure 5.9: Periodic solutions in cell mapping cross-section phase plane for different σ, e = 0.1, k 1 = 2 86 5.4 Effect of spin direction e = 0.1, k 1 = 2 Comparing these results to the results for positive spin direction, it is clear that for σ≈−0.5 there are no families of stable periodic solutions besides the origin. In addition, double-, triple-period of periodic solutions are shown around the origin as σ increases. All unstable P− 1 solutions are specified in Table 5.2 and two of the cases for σ =−0.5 and σ =−0.55 are studied in detail. Table 5.2: Unstable P− 1 solution in Fig.5.9 σ K β(rad) γ(rad) β 0 γ 0 stability -0.5 1 0 -0.533606179 -1.144287256 0 u -0.5 1 0 0.533606179 1.144287256 0 u -0.52 1 0 -0.552144603 -1.198526906 0 u -0.52 1 0 0.552144603 1.198526906 0 u -0.55 1 0 -0.582558403 1.283616124 0 u -0.55 1 0 0.582558403 1.283616124 0 u A stability analysis result for σ =−0.5 is shown in Fig.5.10. All the eigenvalues ofH for the large amplitude solution are on the real axis implying it is unstable, see Fig.5.10a. The eigenvalues for the equilibrium point are on the unit circle as shown in Fig.5.10b, implying that the origin is stable. 87 5.4 Effect of spin direction e = 0.1, k 1 = 2 Figure 5.10: Eigenvalues of matrixH for two P− 1 solutions, σ =−0.5 To verify the analysis result, the continuous dynamics given in Eq. 5.1 was simulated with initial conditions as the P− 1 solution for σ =−0.5. Corresponding trajectories in β−β 0 phase plane and β vs. true anomaly for 1·T are shown in Fig.5.11. The trajectories for ν = T come back to the same point after one revolution. However, long term behavior, for ν = 30·T, indicates that the solutions are not periodic any longer as shown in Fig.5.12. 88 5.4 Effect of spin direction e = 0.1, k 1 = 2 (a) σ =−0.5, Unstable P−1 solution for 1·T (b) σ =−0.55, Unstable P−1 solution for 1·T Figure5.11: Trajectory inβ−β 0 phase plane andβ vs. true anomaly ofP−1 solution for 1·T 89 5.4 Effect of spin direction e = 0.1, k 1 = 2 (a) σ =−0.5, Unstable P−1 solution for 30·T (b) σ =−0.55, Unstable P−1 solution for 30·T Figure5.12: Trajectory inβ−β 0 phase plane andβ vs. true anomaly ofP−1 solution for 30·T 5.4.4 Bifurcations for−0.6<σ<−0.5 for small motion Modi et al. [51] showed that small motions for e = 0.1 and k 1 = 2 are unstable for−1.9<σ <−0.45. This study shows that a rich dynamics exists in this range of spinningparameterasshowninFig.5.13: theoriginisstablefor−0.58<σ<−0.5and different types of bifurcations from the origin with negative direction of spin around σ =−0.5. P− 1 toP− 1,P− 1 toP− 2, andP− 1 toP− 3 bifurcations are found. A detailed study of these bifurcations follows. 90 5.4 Effect of spin direction e = 0.1, k 1 = 2 (a) σ vs. β (b) σ vs. γ Figure 5.13: Bifurcation map around the origin for negative spin (e = 0.1, k 1 = 2) 5.4.5 P - 1 to P - 1 Bifurcation for−0.6<σ<−0.5 for small motion A P− 1 to P− 1 bifurcation at the origin as a function of σ is presented in Fig.5.14. The solid disks indicate stable P− 1 solutions and empty circles indicate unstable P− 1 solutions. As the spinning rate increases fromσ =−0.52, the stable equilibrium point becomes unstable and new stable P− 1 solutions appear at σ≈−0.57. Note that in this study P− 1 solution is defined as |X(T )−X(0)|≈ 10 −5 where X is a state vector, [β,γ,β 0 ,γ 0 ] T . 91 5.4 Effect of spin direction e = 0.1, k 1 = 2 Figure 5.14: Bifurcation sequence P− 1 to P− 1 around σ =−0.58 To verify the P− 1 to P− 1 bifurcation, eigenvalues of H were evaluated at the origin (see Eq. (3.7)) as shown in Fig.5.15. As the magnitude of σ increases, a pair of eigenvalues follows the unit circle in the direction of Re(λ) = 1. The eigenvalues collide at λ = 1, then leave the unit circle along the real axis at σ≈−0.57. Since the system is area preserving, H is symplectic. One eigenvalue moves to the inside and the other one moves to the outside of the unit circle, satisfying λ 1 = 1/λ 2 . The equilibrium at the origin is no longer stable and two new stableP−1 solutions appear. 92 5.4 Effect of spin direction e = 0.1, k 1 = 2 Figure 5.15: Eigenvalues ofH evaluated at the origin, σ from−0.6 to−0.53 Table5.3: Eigenvalue pairsλ i (i = 1,...4) for the origin with varyingσ showingP−1→ P− 1 bifurcation σ λ 1,2 λ 3,4 Stability -0.53 0.995250567105360± 0.0973227184292254i -0.643958708863879± 0.765060102710648i s -0.54 0.996896949617243± 0.0787042524113157i -0.471778970330879± 0.881716012510921i s -0.55 0.998268684854318± 0.0587804010131666i -0.267324884429348± 0.963606804968395i s -0.56 0.999297774856389± 0.0374118004508307i -0.0291045113798667± 0.999576087592256i s -0.57 0.999888269726258± 0.0147689925690726i 0.243311117868802± 0.969948290721545i BFC -0.58 0.997050586268362, 1.00295273855304 0.550222846464345± 0.835017844304221i u -0.59 1.02563395170468, 0.975001497704819 0.891307458139955± 0.453399277183050i u -0.6 2.05448490467412, 0.486741186361070 0.997160070415447± 0.0754476044891167i u Bifurcated stable P− 1 solutions for σ = −0.58 are specified in Table 5.4 and the lower one, as an example, is demonstrated for ν = 1·T and ν = 30·T. The corresponding trajectories in β−β 0 phase plane and β response vs. ν are shown in Fig.5.16. Table 5.4: Bifurcated P− 1 solutions for σ =−0.58 σ P−K β(rad) γ(rad) β 0 γ 0 stability -0.58 1 4.73E-05 0.071540931 0.084438327 -8.60E-05 s -0.58 1 -4.73E-05 -0.071540931 -0.084438327 8.60E-05 s 93 5.4 Effect of spin direction e = 0.1, k 1 = 2 Figure 5.16: Bifurcated P− 1 solution, σ =−0.58 The motion remains periodic with period T = 2π for 30·T. It can also be noticed that the motion possesses additional frequency within ν = 1·T. Fast Fourier analysis was performed for the solution to further study the characteristics and the result is shown in Fig.5.17. Two spikes at ω = 1/T and ω = 2/T are observed as expected. The largest FFT amplitude at ω = 2/T explains the additional frequency observed in Fig.5.16. Figure 5.17: Fast Fourier analysis for stable P− 1 solution, σ =−0.58 94 5.4 Effect of spin direction e = 0.1, k 1 = 2 5.4.6 P− 1 to P− 2 Bifurcation for−0.6<σ<−0.5 for small motion As we recall Fig.5.9, a P− 1 to P− 2 bifurcation at the origin was observed as the magnitude of spinning rate increases from σ =−0.5. For the detailed study, a P− 2 bifurcation diagram as function of σ is shown in Fig.5.18. As the magnitude of σ increases from−0.48, the unstable origin, P− 1 solution, becomes stable and a P− 2 solution evolves at σ≈−0.498. Figure 5.18: P− 1 to P− 2 bifurcation diagram as function of σ where σ≈−0.5 This period-doubling bifurcation is verified by examining the eigenvalues ofH. The eigenvalues about the origin for the above cases are shown in Fig.5.19 and Table 5.5. One of the pairs of the eigenvalues on the unit circle moves toward Re(λ) =−1 as the magnitude of σ decreases. Then they merge at λ =−1 and are separated along the real axis as the bifurcation parameter continues decreasing in speed. 95 5.4 Effect of spin direction e = 0.1, k 1 = 2 Figure 5.19: Eigenvalues ofH with varying spinning parameter−0.52<σ<−0.48 Table5.5: Eigenvalue pairsλ i (i = 1,...4) for the origin with varyingσ showingP−1→ P− 2 bifurcation σ λ 1,2 λ 3,4 Stability -0.48 -1.44706193722838 , -0.691047742406332 0.985800593093026±0.167983007065857i u -0.49 -1.30901000186359 , -0.763936077310019 0.987622463716411±0.156834199806479i u -0.4947 -1.17402187742606 , -0.851772863550970 0.988505383737146±0.151169311862836i u -0.497 -0.999331872937309±0.0365483920998347i 0.988945481282108±0.148262876560312i BFC -0.5 -0.979589941681421±0.201006483660305i 0.989524264357135±0.144349271433463i s -0.513 -0.866009042664189±0.500028870532206i 0.992065971580157±0.125695710934521i s -0.52 -0.785026253748758±0.619462010715720i 0.993417774078554± 0.114529873658189i s BifurcatedP− 2 solutions for−0.53≤σ≤−0.5 are studied in detail. The periodic solutions and the determinant of its linear part of the monodromy matrix H are specified in Table 5.6. As it can be seen|H−I|< 10 −5 indicates the accuracy of the calculation. 96 5.4 Effect of spin direction e = 0.1, k 1 = 2 Table 5.6: P− 2 solutions for σ =−0.51 to−0.53 σ P−K β(rad) γ(rad) β 0 (rad/rad) γ 0 (rad/rad) −0.51 P−2 0.000523289 -0.058733093 0.093107597 -8.95E-05 P−2 -0.000689127 0.058756484 -0.093342598 -0.000263693 −0.52 P−2 0.003187074 -0.078719814 0.118331905 0.000671146 P−2 -0.003187074 0.078719814 -0.118331905 -0.000671146 −0.53 P−2 0.000532972 -0.099808332 0.137467735 0.000380074 P−2 -0.000532972 0.099808332 -0.137467735 -0.000380074 The trajectory in β−β 0 phase plane and β versus ν for the P− 2 solutions are studied forν = 2·T andν = 80·T as shown in Fig.5.20. It shows the solution’s 2·T periodicity but the long term result for ν = 80·T has distinctive characteristics. The long term response is not periodic with 2·T and not diverging, but rather showing quasi-periodic behavior. When σ is between−0.5 and−0.52, the response exhibits beat phenomenon, and there is no obvious pattern for the response for σ =−0.53. 97 5.4 Effect of spin direction e = 0.1, k 1 = 2 (a) 2·T, σ =−0.5 (b) 80·T, σ =−0.5 (c) 2·T, σ =−0.51 (d) 80·T, σ =−0.51 (e) 2·T, σ =−0.52 (f) 80·T, σ =−0.52 (g) 2·T, σ =−0.53 (h) 80·T, σ =−0.53 Figure 5.20: Trajectory inβ−β 0 phase plane (left) andβ response vs. true anomaly (right) To further characterize the solutions, fast Fourier analysis was performed. As shown in Fig.5.21, for σ =−0.5,−0.51, and−0.52, two close magnitude peaks are found close toω = 0.5/T which explains the beat phenomena observed in the simulation. For σ =−0.53, the FFT content is spread over large range of frequencies (see Fig.5.21d, 98 5.4 Effect of spin direction e = 0.1, k 1 = 2 and the motion appeared chaotic as shown in Fig.5.20h. Figure 5.21: Fast Fourier analysis on claimed P− 2 solutions for four different σs 5.4.7 P− 1 to P− 3 Bifurcation for−0.6<σ<−0.5 for small motion P− 1 toP− 3 bifurcation is studied about the origin for−0.56<σ<−0.5 as shown in Fig.5.22. Solid dots are the origin, equilibrium point, and black crosses are P− 3 solutions. As the satellite’s spinning rate increases from σ =−0.5 to σ =−0.56, the P− 3 solution evolves at σ =−0.54. 99 5.4 Effect of spin direction e = 0.1, k 1 = 2 Figure 5.22: P− 1 to P− 3 bifurcation sequence at σ≈−0.54 Three-period bifurcation possesses eigenvalues at 3 rd roots of unity (λ 3 = 1) and this is verified in Fig.5.23. The dotted line in the unit circle indicates the 3 rd roots (1,−0.5± √ 3/2i). It can be seen the eigenvalues are crossing the line at σ≈−0.538. Table 5.7 shows the eigenvalues explicitly. Figure 5.23: Eigenvalues ofH for spinning parameter−0.55<σ<−0.52 100 5.4 Effect of spin direction e = 0.1, k 1 = 2 Table 5.7: Eigenvalues λ 1,...4 for the origin with varying σ showing P− 1→ P− 3 bifurcation σ λ 1,2 λ 3,4 Stability -0.52 -0.784997822610847±0.619498508579087i 0.993420244358505±0.114503993919956i s -0.53 -0.643958708863879±0.765060102710648i 0.995250567105360± 0.0973227184292254i s -0.537 -0.526683235314628±0.850061619045628i 0.996431995768645±0.0843690732249635i s -0.538 -0.508663378944382±0.860965477579103i 0.996591165416195±0.0824676088354798i BFC -0.539 -0.490321635988220±0.871541553837268i 0.996747798713487±0.0805523308488839i s -0.54 -0.471678259291847±0.881770672483893i 0.996900762079159±0.0786378911618127i s -0.55 -0.267002255466185±0.963695897746117i 0.998278938546619±0.0585999429630683i s Table 5.8 gives the exact P− 3 solutions from the cell mapping. These bifurcated P− 3 solutions are verified with their β responses vs. true anomaly and trajectories in phase plane as the other cases. Table 5.8: P− 3 solutions σ K β(rad) γ(rad) β 0 (rad/rad) γ 0 (rad/rad) −0.53 1 -1.48E-09 -1.21E-08 1.21E-08 -3.48E-09 −0.54 3 0.008371525 -0.001023843 0.00156352 0.000512866 3 -0.000464049 0.002148555 -0.003493245 -6.46E-06 3 -0.007894679 -0.001213309 0.001865596 -0.00044302 −0.55 3 0.154342434 -0.006796917 0.046364095 0.009558272 3 -0.003473505 0.054803969 -0.047252179 0.000241366 3 -0.143096487 -0.00936239 0.050269285 -0.00803785 As shown in Fig.5.24 and Fig.5.25, the solution comes back to the same point after 3·T. However, they are not exactP−3 solutions as shown in the longer time response. 101 5.4 Effect of spin direction e = 0.1, k 1 = 2 (a) Trajectory ofP−3 solution in phase plane (left) andβ response vs. true anomaly for 3·T (b) Trajectory ofP−3 solution in phase plane (left) andβ response vs. true anomaly (right) for 100·T Figure 5.24: P− 3 solution for σ =−0.54 102 5.4 Effect of spin direction e = 0.1, k 1 = 2 (a) Trajectory ofP−3 solution in phase plane (left) andβ response vs. true anomaly (right) for 3·T (b) Trajectory ofP−3 solution in phase plane (left) andβ response vs. true anomaly (right) for 100·T Figure 5.25: P− 3 solution for σ =−0.55 The results of fast Fourier analysis have the dominant frequency at around ω = 0.33/T which means period of three. (a) σ =−0.54 (b) σ =−0.55 Figure 5.26: Fast Fourier analysis on P− 3 solutions 103 5.5 Effect of high spin rate 5.5 Effect of high spin rate The attitude dynamics of a satellite with three different spinning rates are demon- strated on a slightly eccentric orbit (e = 0.1) and a disk shape (k 1 = 2) body. Periodic solutions are found and compared in the order of increasing spinning velocity. The corresponding P− 1 solutions are shown in β−γ and β−β 0 point mapping section planes in Fig.5.27. (a) β−γ (b) β−β 0 Figure 5.27: P− 1 solutions in cell mapping stroboscopic phase plane for different σ (e = 0.1, k 1 = 2) As noted earlier, each point on the oval shaped curve represents a periodic solution in the cross section plane, and the origin is an equilibrium point for all the cases. For all three cases, the oval shaped curve belongs to a stable P− 1 solution family. As σ increases, the maximum amplitude of the solutions decreases. Spinning rate σ = 0.5 shown in Fig.5.2 can also be compared with this result. 104 5.5 Effect of high spin rate In order to study the characteristics of the solution family, a representative solution, specified in Table 5.9 is chosen from each family. Table 5.9: Examples of P− 1 solutions for larger σ (e = 0.1, k 1 = 2) σ P−K β(rad) γ(rad) β 0 (rad/rad) γ 0 (rad/rad) 1 P−1 -0.093123155 -0.555035636 1.868299865 -0.362434439 1.5 P−1 0.146141786 -0.421333233 1.808427853 0.698041729 2 P−1 -0.148915114 0.306770863 -1.594307473 -0.845948474 As can be observed in Fig.5.28, each example of an initial condition is integrated for 1·T and 14·T and showed the solutions have a periodT = 2π. The initial and the final points are marked with a small circle in the left β−β 0 plane with trajectories. β response versus true anomaly is presented on the right side. As the spinning parameter increases, the amplitude of motion decreases. In addition, it can be seen that higher frequency oscillations are present in the response. 105 5.5 Effect of high spin rate (a) σ = 1, P−1 solution for 1·T (b) σ = 1, P−1 solution for 14·T (c) σ = 1.5, P−1 solution for 1·T (d) σ = 1.5, P−1 solution for 14·T (e) σ = 2, P−1 solution for 1·T (f) σ = 2, P−1 solution for 14·T Figure 5.28: Trajectory inβ−β 0 phase plane (left) andβ response vs. true anomaly (right) of P− 1 solution These quasi-periodic motions are also found in the following fast Fourier analysis. Fourier analysis shown in Fig.5.29 confirms the simulation observations where su- perharmonic solutions with frequencies that are integer multiples of the fundamental frequency 1/T appear. The peak magnitude is shifting to the right which means the satellite oscillates faster. 106 5.5 Effect of high spin rate (a) σ = 1 (b) σ = 1.5 (c) σ = 2 Figure 5.29: Fast Fourier analysis on P− 1 solutions 107 5.6 Effect of orbital eccentricity 5.6 Effect of orbital eccentricity The effect of the orbital eccentricity on the satellite’s behavior is studied for three different orbital eccentricities (0.1 < e < 0.2) with σ = 0.5, k 1 = 2. As shown in Fig.5.30, two families of stable P− 1 solutions are found for each eccentricity value. (a) x 1 −x 2 phase plane (b) x 1 −x 3 phase plane Figure 5.30: P− 1 solutions in stroboscopic phase plane for three different orbital eccentricities The larger eccentricity we considered, the smaller maximum angular amplitude for both families. A representative solution for each eccentricity was studied in detail. One point of each case, as an example, is specified in Table 5.10 and studied in detail for its periodicity and stability analysis. 108 5.6 Effect of orbital eccentricity Table 5.10: Examples of P− 1 solutions for larger e, σ = 0.5 e P−K β(rad) γ(rad) β 0 (rad/rad) γ 0 (rad/rad) 0.1 P−1 0.608247123 0.108148065 -0.307642323 2.506594604 0.15 P−1 -0.44285226 0.107903366 -0.255600905 -1.560618597 0.2 P−1 -0.140297319 -0.027008624 0.051148881 -0.43206838 The corresponding trajectories inβ−β 0 phase plane and theβ response as function ofν forT and 14·T are shown in Fig.5.31. The solutions are periodic with period T, and as the previously studied solutions, higher frequencies thanω = 1/T are observed. (a) e = 0.15, P−1 solution for 1·T (b) e = 0.15, P−1 solution for 14·T (c) e = 0.2, P−1 solution for 1·T (d) e = 0.2, P−1 solution for 14·T Figure5.31: Phase plane trajectory andβ response vs. true anomaly ofP−1 solution The fast Fourier analysis results demonstrate this quasi-periodic motion. As the FFT analysis shown in ??, a combination of frequencies ofω =n/T wheren = 1, 2,.. is detected, i.e. periods of T/n are present in the solution. Larger eccentricity leads to higher frequency contents. 109 5.6 Effect of orbital eccentricity (a) e = 1.5 (b) e = 2 Figure 5.32: Fast Fourier analysis on P− 1 solutions 110 5.7 Effect of satellite shape 5.7 Effect of satellite shape The effect of the inertia parameter k 1 , that characterizes satellite’s shape, on the satellite motion was carried out. The satellite is a sphere shaped body when k 1 = 1 with no effect of gravity gradient torque. The case of k 1 = 2 was studied so far which tends to be more stable than where k 1 < 1. The relationship between this sort of instability and the inertia parameter k 1 is examined. In Fig.5.33, the stability property of the equilibrium point for three different spin rates where 0<k 1 < 2 is shown. Figure 5.33: Stability analysis at the origin, 0<k 1 < 2 with σ = 0.5, 0,−0.5 All three of spinning parameter σ, the equilibrium point is mostly unstable when k 1 < 1, pencil shape body. For σ =−0.5, the equilibrium point is mostly unstable even for k 1 > 1. For σ = 0, there is a small region in k 1 > 1 where the equilibrium point is unstable, while it is unstable for the entire region in k 1 < 1. For σ = 0.5, the origin is stable when k 1 > 1 but unstable at k 1 = 1.25. Here, in order to investigate the effect of satellite’s shape, we consider only the case σ = 0.5. Since there are significant stability differences between larger and smaller than k 1 = 1, the study is grouped accordingly. 111 5.7 Effect of satellite shape 5.7.1 k 1 < 1 As the first part, the results for k 1 < 1 are described. Two pairs of large amplitude unstable P− 1 solutions (marked as a and b) are found as shown in Fig.5.34. For clearer visualization of the solutions, they are shown in both β−γ and β−β 0 cell mapping cross section planes. The solutions are symmetric about β = 0 and γ = 0 axes, respectively, and the amplitude of motion increases with k 1 . (a) β−γ angular position plane (b) β−β 0 phase plane Figure 5.34: P− 1 solutions in stroboscopic phase plane (e = 0.1,σ = 0.5) All the solutions, include the equilibrium point, are unstable. The eigenvalues ofH evaluated along these solutions are shown in Fig.5.35. 112 5.7 Effect of satellite shape (a) k 1 = 0.5 (b) k 1 = 0.75 Figure 5.35: Stability analysis of equilibrium point, k 1 < 1 (unstable), σ = 0.5 It is clear that for all these cases, the solutions are unstable since one of the eigen- values is outside the unit circle. To be more specific, one of each paired solutions is specified in Table 5.11 and integrated for one period. Table 5.11: Examples of P− 1 solutions for different k 1 , σ = 0.5 k 1 P−K index β(rad) γ(rad) β 0 (rad/rad) γ 0 (rad/rad) 0.5 P−1 0.5a -2.16E-05 0.37844006 0.979020441 1.02E-05 0.5 P−1 0.5b -0.331544855 6.30E-05 7.43E-05 1.255389156 0.75 P−1 0.75a 3.74E-11 0.468174869 0.198074111 3.50E-09 0.75 P−1 0.75b -0.554438091 0.000116346 0.000134571 1.675079027 113 5.7 Effect of satellite shape (a) Initial condition with ’0.5a’ in Table 5.11 (b) Initial condition with ’0.5b’ in Table 5.11 Figure 5.36: Trajectory in β−β 0 phase plane and β response vs. true anomaly of P− 1 solution for 1·T, k 1 = 0.5 (e = 0.1, σ = 0.5) 114 5.7 Effect of satellite shape (a) Initial condition with ’0.75a’ in Table 5.11 (b) Initial condition with ’0.75b’ in Table 5.11 Figure 5.37: Trajectory in β−β 0 phase plane and β response vs. true anomaly of P− 1 solution for 1·T, k 1 = 0.75 (e = 0.1, σ = 0.5) As shown in Fig.5.36 and Fig.5.37, the solutions are periodic with period ofT. The maximum β and γ positions for ‘0.5a’ are 2 rad and 1.1 rad, respectively, whereas for ‘0.75a’ they are 0.71 rad and 0.25 rad, respectively. Note that, since these P− 1 solutions are unstable, they diverge after some revolutions. 5.7.2 k 1 > 1 For the next part, the results for k 1 > 1 are described. Again, P− 1 solutions are found using Cell mapping and they are shown in Fig.5.38. When k 1 > 1, the large amplitudeP− 1 solutions are stable forming a family of periodic solutions. Again, the points in Fig.5.38 presentP− 1 solutions inβ−γ plane (upper one) andβ−β 0 plane 115 5.7 Effect of satellite shape (lower one), i.e. each point in the figure repeats itself every 1·T. Ask 1 increases from 1.25 to 2, angular position β and γ tend to be symmetric in the position plane. (a) β−γ angular position plane (b) β−β 0 phase plane Figure 5.38: P− 1 solutions in stroboscopic phase plane fork 1 > 1 (e = 0.1,σ = 0.5 ) The solutions periodicity can be shown with their response by orbital resolutions. One point of each family is selected as an initial condition and its trajectory for 1·T is demonstrated. 116 5.7 Effect of satellite shape (a) k 1 = 1.25 (b) k 1 = 1.5 Figure 5.39: Trajectory in β−β 0 phase plane and β response vs. true anomaly of P− 1 solution for 1·T, k 1 > 1 (e = 0.1, σ = 0.5) 5.7.3 Special case for k 1 = 1.25 Considering the stability chart for the equilibrium point (see Fig.5.33) again, it is observed an unstable region for σ = 0.5, at k 1 = 1.25. The detailed dynamics is studied for understanding its stability transition. AP− 1 toP− 2 bifurcation is newly found for the region. Comparing Fig.5.33 and Fig.5.40, the stable equilibrium point is bifurcated to unstable equilibrium point and P− 2 solution at k 1 = 1.25 as eigenvalues ofH on the unit circle move along to the real axis through Re(λ) =−1. 117 5.7 Effect of satellite shape Figure 5.40: Eigenvalues evaluated at the origin 1.25<k 1 < 1.27 (left: magnified on the crowed region near λ =−1, right: full scale of unit circle) The bifurcatedP− 2 solutions are specified in Table 5.12 to verify their periodicity and trajectory. Since these two solutions are paired and symmetric, the one with positive β, as an example, is demonstrated. Table 5.12: Examples of P− 1 solutions for different k 1 , σ = 0.5 k 1 P−K β(rad) γ(rad) β 0 (rad/rad) γ 0 (rad/rad) 1.25 P−2 0.010926162 -3.05E-06 1.57E-06 0.020992411 1.25 P−2 -0.010926162 3.05E-06 -1.57E-06 -0.020992411 With this solution as an initial condition, Fig.5.41 shows its 2·T and 14·T responses in β−γ angular position plane and β−β 0 phase plane. It is clearly shown that the solution repeats after 2·T and there exists superposed oscillations. 118 5.7 Effect of satellite shape (a) Response trajectory for 2·T (b) Response trajectory for 14·T Figure 5.41: Trajectory in β−β 0 phase plane and β response vs. true anomaly of P− 2 solution for k 1 = 1.25 (e = 0.1, σ = 0.5) Fast Fourier analysis shown in Fig.5.42 agrees with the numerical analysis. Figure 5.42: Fast Fourier analysis on P− 2 solutions, k 1 = 1.25 (e = 0.1, σ = 0.5) 119 5.8 Summary Thepeakmagnitudesatω = (0.5·n)/T acombinationwheren = 1, 3, 5,aredetected. This presents P− 2 motion with a combination of sub-motions in the solution. 5.8 Summary The periodic motion analysis for spinning satellite was analyzed by a new method that combines numerical and analytical approach. Global analysis was performed, and new types of periodic solutions were found. Stable and unstable families ofP− 1 solutions were analyzed for various combinations of the system parameters. It was found that the direction of satellite’s spin as well as its speed affects the stability of the satellite. P− 1 to P− 1 and P− 1 to P− 2 bifurcations as a function of spinning rate and inertia parameter were conducted. Eigenvalues of the monodromy matrix obtained by point mapping method were used to verify the stability of solutions and bifurcation conditions. Further characterization of the solutions was performed using fast Fourier analysis showing inherent motion frequencies. 120 6 Conclusion and Future Research 6.1 Summary and Conclusion Unveiling the dynamics of spacecraft under the influence of gravity is important to design the spacecraft and its control system in the environment. In this dissertation, two classes of problems, gravity gradient torque stabilized satellites and dynamics of spinning satellites, were discussed. For the study, an approach with point mapping and cell mapping that combines numerical and analytical analysis was introduced. The analysis was performed globally, and it was able to demonstrate the previous results and find new families of periodic solutions. The stability of the solutions was examined analytically using point mapping technique and fast Fourier analysis. The invariant surfaces which give ideas on how the periodic solutions are emerging and behave were found by using cell mapping analysis. Moreover, the bifurcation diagram which is important to study a route to chaos was computed. 6.2 Future Research The results indicate that the method can be applied to various systems for global motion analysis including periodic solutions, double- and triple- (or more) periods of periodic solutions, and global regions of attraction. The study of such problems as dy- 121 6.2 Future Research namics of asteroids and spacecraft attitude models in which coupling between orbital dynamics and attitude occurs has great value for understanding complex astrodynam- ics models. Periodic solutions in a wide range of dynamics models will be useful for spacecraft design requirements as well. 6.2.1 Coupled Orbit-Attitude Dynamics Model We start with the attitude motion of a satellite within the context of the three body problem. ThespacecraftinthemutualgravityfielddepictedinFig.6.1willbeanalyzed using the proposed methods. Ashenberg [25] studied pitch dynamics at the Lagrange points in the elliptic three body problem. Some of periodic solutions were discussed but it was confined to circular problem and at the libration points. Knutson et al. [4] analyzed the stability of spacecraft dynamics in the Lyapunov orbits. Figure 6.1: Geometry of coupled orbit-attitude dynamics model in the three body problem The satellite motion in the three body problem is described using Euler’s equation of motion with kinematic relations as 122 6.2 Future Research I 1 I ˙ ω 1 B + (I 3 −I 2 ) I ω BI 3 ω B 2 =g 1 I 2 I ˙ ω 2 B + (I 1 −I 3 ) I ω BI 1 ω B 3 =g 2 I 3 I ˙ ω 3 B + (I 2 −I 1 ) I ω BI 1 ω B 2 =g 3 where the gravitational torque ¯ g by the primary and the secondary bodies ¯ g = 3μ 1 r 3 1         (I 3 −I 2 )ˆ r 12 ˆ r 13 (I 1 −I 3 )ˆ r 13 ˆ r 11 (I 2 −I 1 )ˆ r 11 ˆ r 12         + 3μ 2 r 3 2         (I 3 −I 2 )ˆ r 22 ˆ r 23 (I 1 −I 3 )ˆ r 23 ˆ r 21 (I 2 −I 1 )ˆ r 21 ˆ r 22         for ˆ r i = h e R i         x p ±D i y p z p         · 1 |¯ r i | . 6.2.2 Asteroid Model The recent successful missions to asteroids such as the NEAR mission(1996), Hayabusa 1(2003), andHayabusa2(2014, ongoing)inspiredustostudythesmallbodydynamics in more detail to expand mission possibilities. Some of the periodic solutions for the small body dynamics were studied by Giancotti et al. [30] and Katerine et al. [58] and Villac et el. [5]. However, a relatively small number of studies for the small body dynamics were performed due to the limitations of analysis methods. This study will contribute to find useful periodic solutions around asteroids as appropriate for new mission design requirements. 123 Bibliography [1] A.A.Zevin. Oscillations of a satellite in the plane of an elliptical orbit. Cosmic research, 19:455–460, 1981. [2] A.Celletti and V.Sidorenko. Some properties of the dumbbell satellite attitude dynamincs. Celest Mech Dyn Astr, 101:105–126, 2008. [3] A.J.Fleig. On the libration of a gravity gradient stabilized spacecraft in an ec- centric orbit. Technical report NASA TR R-350, National Aeronautics and space administration, Nov. 1970. [4] A.Knutson, D.Guzzetti, K.C.Howell, and M.Lavagna. Attitude responses in cou- pled orbit-attitude dynamical model in earth-moon lyapunov orbits. Journal of Guidance, Control, and Dynamics, 38(7):1264–1273, 2015. [5] B.F.Villac, R.L.Anderson, and A.J.Pini. Organizing ballistic orbit classes around smallbodies. InAAS/AIAAAstrodynamicsSpecialistConference, Vail, Colorado, August 9-13, 2015, number AAS 15-619, 2015. [6] C.S.Hsu. On non-linear parametric excitation problems. Advanced applied me- chanics, 17:245, 1977. 124 Bibliography [7] C.S.Hsu. A theory of cell-to-cell mapping dynamical systems. Journal of Applied Mechanics, 47:931–939, December 1980. [8] C.S.Hsu. An unravelling algorithm for global analysis of dynamical systems: An application of cell-to-cell mappings. Journal of Applied Mechanics, 47:940–948, December 1980. [9] C.S.Hsu. A theory of cell-to-cell mapping dynamical systems. Journal of Applied Mechanics, 47(4):931–939, December 1980. [10] C.S.Hsu. An unravelling algorithm for global analysis of dynamical systems: An application of cell-to-cell mappings. Journal of Applied Mechanics, 47(4):940–948, December 1980. [11] D.Fischer. The spectacular journey of the Galilo spacecraft. Springer Science & Businees Media, 2013. [12] D.J.McGill and L.S.Long. Stability regions for an unsymmetrical rigid body spin- ning about an axis through a fixed point. Acta Mechanica, 22:91–112, 1975. [13] D.Koh and H.Flashner. Global analysis of gravity gradient satellite’s pitch motion in an elliptic orbit. Journal of Computational and Nonlinear Dynamics, 10(6):1– 11, November 2015. [14] D.L.Hitzl and D.A.Levinson. Attitude stability of a spinning symmetric satellite in a planar periodic orbit. Celest Mech Dyn Astr, 20:179–205, 1979. [15] D.S.Roldugin and P.Testani. Spin-stabilized satellite magnetic attitude control scheme without initial detumbling. Acta Astronautica, 94(1):446–454, 2014. 125 Bibliography [16] F.B.Wallace and L.Meirovitch. Attitude instability regions of a spinning symmet- ric satellite in an elliptic orbit. AIAA Jounal, 5(9):1642–1650, September 1967. [17] R.V. Garcia, M.C. Zanardi, and H.K.Kuga. Spin-stabilized spacecraft: Ana- lytical attitude propagation using magnetic torques. Mathematical Problems in Engineering, 2009:1–18, 2010. [18] G.Frost. Attitude orientation control for a spinning satellite. Technical Report RAND/N-3250-AF, 1991. [19] H.A.Karasopoulos. Nonlinear Dynamics of the Planar Pitch Attitude Motion for a Gravity-Gradient Satellite. PhD thesis, University of Cincinnati, 1994. [20] H.B.Schechter. Satellite librations on an elliptic orbit. Rand Corp. RM-3632-PR, May 1963. [21] H.B.Schechter. Dumbbell librations in elliptic orbits,. AIAA Jounal, 2:1000–1003, 1964. [22] H.Flashner and C.S.Hsu. A study of nonlinear periodic systems via the point mapping method. International Journal for Numerical Methods in Engineering, 19:185–215, 1983. [23] H.Flashner and C.S.Hsu. A study of nonlinear periodic systems via the point mapping method. International Journal for Numerical Methods in Engineering, 19(2):185–215, 1983. [24] H.Flashner and R.S.Guttalu. Analysis of non-linear non-autonomous systems by truncated point mappings. International Journal of Non-Linear Mechanics, 24 (4):327–344, 1989. 126 Bibliography [25] J.Ashenberg. Satellite pitch dynamics in the elliptic problem of three bodies. Journal of guidance, control, and dynamics, 19(1):68–74, January-February 1996. [26] J.P.Moran. Effect of plane librations on the ororbit motion of a dumbbel satellite. ARS Journal, pages 1089–1096, August 1961. [27] J.R.Wertz. Spacecraft attitude determination and control (Astrophysics and space science library). D.Reidel, 1978. [28] J.V.Breakwell and R.Pringle. Nonlinear resonance affecting gravity-gradient sta- bility. Proc. XVIth International Astronautical Congress, Athenes,1965/ Parix: Gauthier-Villars, pages 305–325, 1966. [29] J.W.Cole and R.A.Calico. Nonlinear oscillations of a controlled periodic system. Journal of Guidance, Control, and Dynamics, 15:627–633, 1992. [30] M.Giancotti, S.Campagnola, Y.Tsuda, and J.Kawaguchi. Families of periodic orbits in hill’s problem with solar radiation pressure: application to hayabusa 2. Celest. Mech. Dyn. Astr., 120:269–286, 2014. [31] M.Golat and H.Flashner. A new methodology for the analysis of periodic systems. Nonlinear Dynamics, 28:29–51, 2002. [32] M.H.Kaplan. Modern Spacecraft Dynamics and Control. John Wiley & Sons, Inc.,New York, 1976. [33] M.J.Sidi. Spacecraft Dynamics and Control: A Practical Engineering Approach. Cambridge University Press, 2000. [34] NASA. Spacecraft gravitational torques. Technical Report SP-8024, May 1969. 127 Bibliography [35] O.H.Gerlach. Attitude stabilization and control of earth satellites. Space science reviews, 4(4):541–582, 1965. [36] P.C.Hughes. Spacecraft Attitude Dynamics. John Wiley & Sons, Inc.,New York, 1986. [37] P.F.Hultquist. Gravitationaltorqueimpulseonastabilizedsatellite. ARSJournal, 31:1506–1509, Nov. 1961. [38] R.A.Nidey. Gravitational torque on a satellite of arbitrary shape. ARS, 30:203– 204, 1960. [39] R.E.Fischell. Torques and attitude sensing in earth satellite, chapter Passive gravity-gradient stabilization for earth satellites, pages 12–29. Academic Press, New York and London, 1964. [40] R.E.Roberson. Gravitational torque on a satellite vehicle. Journal of the Franklin Institute, 265:13–22, 1958. [41] R.M.L.Baker. Librations on a slightly eccentric orbit. ARS Journal, 30:124–126, Jan. 1960. [42] R.O.Fimmel, W.Swindell, and E.Burgess. Nasa publication sp-349/396. Technical report, NASA, Washington, 1977. [43] R.R.Bate, D.Mueller, and J.E.White. Fundamentals of astrodynamics. Dover, 1971. [44] R.S.Guttalu and H.Flashner. Stability analysis of periodic systems by truncated point mapping. Journal of Sound and Vibration, 189:33–54, 1996. 128 Bibliography [45] R.W.Rarquhar. The flight of isee-3/ice: origins, mission history, and a legacy. Journal of the Astronautical Sciences, 49(1):23–74, 2001. [46] T.R.Kane. Attitude stability of earth-pointing satellites. AIAA Jounal, 3(3): 726–731, APRIL 1965. [47] T.R.Kane and P.M.Barba. Attitude stability of a spinning satellite in an elliptic orbit. Journal of Applied Mechanics, pages 402–405, June 1966. [48] V.A.Sarychev, V.V.Sazonov, and V.A.Zlatoustov. Periodic oscillations of a satel- lite in the plane of an elliptic orbit. Cosmic research, 15(6):698–719, 1977. [49] V.A.Zlatoustov and A.P.Markeev. Stability of panar oscillations of a satellite in an elliptic orbit. Celestial Mechanics, 7:31–45, 1973. [50] V.A.Zlatoustov, D.E.Okhotsimsky, V.A.Sarychev, and A.P.Torzhevsky. Investiga- tion of a satellite oscillations in the plane of an elliptic orbit. Applied Mechanics; proceedings of the eleventh international congress of applied mechanics, Munich (Germany), pages 436–439, 1964. [51] V.J.Modi and J.E.Neilson. On the periodic solutions of slowly spinning gravity gradient system. Celestial Mechanics, 5:126–143, 1972. [52] V.J.Modi and K.J.Pande. On periodic solutions and resonance of spinning satel- lites in near-circular orbits. Celestial Mechanics, 11:195–212, 1975. [53] V.J.Modi and R.C.Brereton. Periodic solutions associated with the gravity- gradient-oriented system: Part 1. analytical and numerical determination. AIAA Jounal, 7(7):1217–1225, July 1969. 129 Bibliography [54] V.J.Modi and R.C.Brereton. Periodic solutions associated with the gravity- gradient-oriented system: Part 2. stability analysis. AIAA Jounal, 7(8):1465– 1468, 1969. [55] V.V.Beletskii. The Satellite Motion about Center of Mass. Publishing House Science, Moscow, Russia, 1965. [56] V.V.Beletsky and E.M.Levin. Dynamics of space tether Systems, volume Sci 83 of Adv. Astronaut. AAS Publication, 1993. [57] X.Tong and F.Rimrott. Numerical studis on chaotic planar motion of satellite in an elliptic orbit. Chaos, Solitions and Fractals, 1(2):179–186, 1991. [58] Y.Katherine and B.F.Villac. Periodic orbits families in the hill’s three-body prob- lem with solar radiation pressure. Advances in the Astronautical Sciences Series, 136, 2010. 130

Abstract (if available)

Abstract In this dissertation, the complete dynamics of a spacecraft in a general orbit is studied. For this study, a new computational approach which combines analytical and numerical analysis is introduced to investigate the global behavior of the full nonlinear dynamical models of satellites. This approach consists of two methods: point mapping and cell mapping. The point mapping algorithm is employed to study local stability characteristics and bifurcation conditions, and the cell mapping method is utilized to analyze periodic solutions and their global regions of attraction. The proposed approach does not depend on the assumption of a small parameter and, therefore, no limitations are imposed on the magnitudes of parameters or amplitude of motion. ❧ This approach is applicable to a wide range of applications. Two classes of satellites are investigated for their global attitude dynamics: gravity gradient satellites and spinning satellites in the elliptic orbits. First, a gravity gradient satellite as an example of a one degree of freedom problem is analyzed. Then it is extended to a two degree of freedom problem in which coupling between spinning motion and pitch motion occurs. As a result, families of periodic solutions and rich dynamic phenomena, which have not been found before, are investigated. Global regions of attraction and evolution of periodic solutions are demonstrated using invariant surfaces, and stability characteristics of the solutions are studied. A comprehensive investigation of bifurcation diagrams was also performed for studying transitions to various types of bifurcations and chaos with varying parameters.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Techniques for analysis and design of temporary capture and resonant motion in astrodynamics

PDF

Relative-motion trajectory generation and maintenance for multi-spacecraft swarms

PDF

Spacecraft trajectory optimization: multiple-impulse to time-optimal finite-burn conversion

PDF

Cooperative localization of a compact spacecraft group using computer vision

PDF

New approaches to satellite formation-keeping and the inverse problem of the calculus of variations

PDF

Application of the fundamental equation to celestial mechanics and astrodynamics

PDF

Designing an optimal software intensive system acquisition: a game theoretic approach

PDF

Optimal guidance trajectories for proximity maneuvering and close approach with a tumbling resident space object under high fidelity J₂ and quadratic drag perturbation model

PDF

Trajectory mission design and navigation for a space weather forecast

PDF

Experimental and numerical investigations of charging interactions of a dusty surface in space plasma

PDF

Controlled and uncontrolled motion in the circular, restricted three-body problem: dynamically-natural spacecraft formations

PDF

Interactions of planetary surfaces with space environments and their effects on volatile formation and transport: atomic scale simulations

PDF

Control of spacecraft with flexible structures using pulse-modulated thrusters

PDF

Particle-in-cell simulations of kinetic-scale turbulence in the solar wind

PDF

New approaches in modeling and control of dynamical systems

PDF

Kinetic studies of collisionless mesothermal plasma flow dynamics

PDF

Motions and conformations of nucleic acids studied using site-directed spin labeling

PDF

Numerical and experimental investigations of dust-plasma-asteroid interactions

PDF

Numerical and experimental investigations of ionic electrospray thruster plume

PDF

Laboratory investigations of the near surface plasma field and charging at the lunar terminator

Asset Metadata

Creator Koh, Dayung (author)

Core Title Periodic motion analysis in spacecraft attitude dynamics

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Astronautical Engineering

Publication Date 03/17/2016

Defense Date 01/28/2016

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag Attitude,bifurcation,cell mapping,dynamics,gravity gradient satellite,OAI-PMH Harvest,periodic solution,point mapping,spin-stabilized satellite

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Flashner, Henryk (committee chair), Anderson, Rodney L. (committee member), Erwin, Dan (committee member), Gruntman, Mike (committee member), Sacker, Robert (committee member)

Creator Email bbjupiter@gmail.com,dayungko@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c40-222196

Unique identifier UC11276448

Identifier etd-KohDayung-4210.pdf (filename),usctheses-c40-222196 (legacy record id)

Legacy Identifier etd-KohDayung-4210.pdf

Dmrecord 222196

Document Type Dissertation

Format application/pdf (imt)

Rights Koh, Dayung

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA