Effects of Rotational Motion on the Ballistic Coefficient of Space Debris by Thomas Benjamin Walsh A thesis submitted to the Graduate Faculty of Auburn University in partial fulfillment of the requirements for the Degree of Master of Science Auburn, Alabama August 4, 2012 Approved by Dr. John Cochran, Jr., Chair, Professor of Aerospace Engineering Dr. David Cicci, Professor of Aerospace Engineering Dr. Andrew Sinclair, Associate Professor of Aerospace Engineering ii Abstract Accurate predication of the future states of operational spacecraft and space debris are necessary for conjunctional analysis. Prediction of the states of debris for a few orbital periods is possible with improved models of the upper atmosphere if the objects have spherical geometries. However, spacecraft and other objects with complex geometries that tumble throughout their orbit present a difficult problem because the rotational motion affects their orbital motion and vice versa. In this thesis, the coupled translational and rotational motions of objects, specifically space debris, are studied using a digital simulation based on a six-degree-of-freedom, rigid body model. In particular, focus is on variations due to rotation of an object's ballistic coefficient, which is the product of a coefficient of drag and a reference area divided by the object's mass. If the density is well modeled, the ballistic coefficient is the principal unknown in the drag force. The results of numerous simulations show a predictable relationship between an object's ballistic coefficient and its nodal regression. The simulation is also used to produce data for an orbit determination process in which the ballistic coefficient is estimated. Results are presented that show continuous estimation of the ballistic coefficient, which is fundamental to accurately predicting future states of debris. iii Acknowledgements The author wishes to express his most sincere appreciation for his advisor, Dr. John Cochran, Jr., for his valuable guidance throughout the development of this thesis. Dr. Cochran's unparalleled expertise in education and the subject matter covered in this thesis have given the author a rewarding experience. He would also like to thank the members of his committee, Dr. David Cicci and Dr. Andrew Sinclair, for their helpful insight throughout the research and review processes. The author would like to thank his family and friends in both Arkansas and Alabama for their patience and continuous support. iv Table of Contents Abstract ...............................................................................................................................ii Acknowledgments .............................................................................................................iii List of Symbols .................................................................................................................. vi List of Figures ...................................................................................................................vii Chapter 1: Introduction .................................................................................................. 1 1.1 Background .......................................................................................................... 1 1.2 Focus of this Research ......................................................................................... 2 Chapter 2: Review of Literature and Low Earth Orbit .................................................. 3 2.1 Space Debris .......................................................................................................... 3 2.2 The Space Environment ......................................................................................... 4 2.2.1 Forces........................................................................................................... 4 2.2.2 Torques......................................................................................................... 6 Chapter 3: Development of the Analysis Program ........................................................ 8 3.1 Introduction.......................................................................................................... 8 3.2 Equations of Motion............................................................................................. 8 3.2.1 Translational Equations................................................................................ 9 3.2.2 Rotational Equations .................................................................................. 18 3.3 Effective Area .................................................................................................... 25 3.4 Numerical Integration ........................................................................................ 27 v 3.4.1 Integration Performance............................................................................. 27 Chapter 4: Analysis Theory and Results...................................................................... 30 4.1 Analysis Approach............................................................................................. 30 4.1.1 Average Effective Area.............................................................................. 31 4.1.2 Existing Ballistic Coefficient Data............................................................. 34 4.2 Relationship of Effective Area Variations with Inclination............................... 37 4.3 Impact Analysis.................................................................................................. 52 Chapter 5: Orbit Determination ................................................................................... 58 5.1 Methodology ...................................................................................................... 58 5.1.1 Data Generation Program........................................................................... 58 5.1.2 Power-Density Matrix and Process Noise.................................................. 59 5.2 Orbit Determination Results .............................................................................. 67 Chapter 6: Conclusions ............................................................................................... 65 Bibliography...................................................................................................................... 67 Appendix A: ECI Coordinate System and Alternate Relative Motion Equations Derivation.......................................................................................................................... 69 Appendix B: Attitude Dynamics................................................................................... 71 vi List of Symbols X ECI X-position Y ECI Y-position Z ECI Z-position i? x-direction unit vector j? y-direction unit vector k? z-direction unit vector r Geocentric range rr Geocentric position vector r? Geocentric position unit vector Vr ECI Velocity relV r Relative velocity between object and atmosphere vX, Y, Z Velocity components ar ECI Acceleration a X, Y, Z Acceleration components draga r Acceleration due to drag 2Ja r Acceleration due to zonal harmonic J 2 ? Standard gravitational parameter G Universal gravitational constant m Mass of object m1 Mass of primary source of gravitational attraction m2 Mass of object undergoing gravitational attraction Re Radius of Earth U Gravitational potential U2 Gravitational potential due to zonal harmonic J2 Pnm Legendre polynomial Jn Zonal harmonic coefficient Fr Force on object Fr Force due to drag CD Drag coefficient A Drag area, Effective Area Atmospheric density Torque on object Torque due to gravity-gradient vii List of Figures Figure 1: Earth-centered inertial coordinate system with orbital elements....................... 10 Figure 2: Orbital plane with semi-major axis, a, occupied focus, F, and perigee, P....... 10 Figure 3: ECI coordinate system with two bodies ............................................................ 12 Figure 4: ECI Coordinate System with Spherical Coordinates......................................... 15 Figure 5: Body centered coordinate system...................................................................... 18 Figure 6: Rotation of body centered coordinate system through Euler angles ................. 19 Figure 7: Coordinate System for Gravity-Gradient Torque.............................................. 23 Figure 8: Performance Indicators vs. Integration Time .................................................... 28 Figure 9: Spacecraft Energies vs. Time ............................................................................ 29 Figure 10: Effective Area vs. Time Over 1 Orbit for Cylinder......................................... 32 Figure 11: Average Effective Area Overlaid on Instantaneous Effective Area for Cylinder............................................................................................................................. 33 Figure 12: Ballistic Coefficient vs. Time for Satellite 00461 ........................................... 35 Figure 13: Ballistic Coefficient vs. Time for Satellite 12184 ........................................... 35 Figure 14: Ballistic Coefficient vs. Time for Satellite 10230 ........................................... 36 Figure 15: Ballistic Coefficient vs. Time for Satellite 20857 ........................................... 36 Figure 16: Average Area & Sine of R.A.A.N. vs. Time for 10? Inclination for Cylinder..................................................................................................... 40 Figure 17: Average Area & Sine of R.A.A.N. vs. Time for 30? Inclination for Cylinder..................................................................................................... 40 viii Figure 18: Average Area & Sine of R.A.A.N. vs. Time for 45? Inclination for Cylinder.................................................................................................... 41 Figure 19: Average Area & Sine of R.A.A.N. vs. Time for 60? Inclination for Cylinder..................................................................................................... 41 Figure 20: Average Area & Sine of R.A.A.N. vs. Time for 80? Inclination for Cylinder..................................................................................................... 42 Figure 21: Average Area & Sine of R.A.A.N. vs. Time for 10? Inclination for Cone ..... 43 Figure 21: Average Area & Sine of R.A.A.N. vs. Time for 10? Inclination for Cone ..... 43 Figure 23: Average Area & Sine of R.A.A.N. vs. Time for 45? Inclination for Cone ..... 44 Figure 24: Average Area & Sine of R.A.A.N. vs. Time for 60? Inclination for Cone ..... 44 Figure 25: Average Area & Sine of R.A.A.N. vs. Time for 80? Inclination for Cone ..... 45 Figure 26: Average Area & Sine of R.A.A.N. vs. Time for 10? Inclination for Flat Plate ................................................................................................... 46 Figure 27: Average Area & Sine of R.A.A.N. vs. Time for 30? Inclination for Flat Plate ................................................................................................... 46 Figure 28: Average Area & Sine of R.A.A.N. vs. Time for 45? Inclination for Flat Plate ................................................................................................... 47 Figure 29: Average Area & Sine of R.A.A.N. vs. Time for 60? Inclination for Flat Plate ................................................................................................... 47 Figure 30: Average Area & Sine of R.A.A.N. vs. Time for 80? Inclination for Flat Plate ................................................................................................... 48 Figure 31: Average Area & Sine of R.A.A.N. vs. Time for 45? Inclination for Cylinder with Random Initial Attitude Conditions for 250 Cases................................................... 49 Figure 32: Average Area & Sine of R.A.A.N. vs. Time for 45? Inclination for Cylinder with Random True Anomaly for 250 Cases...................................................................... 50 Figure 33: Average Area & Sine of R.A.A.N. vs. Time for 45? Inclination for Cone with Random True Anomaly and Initial Attitude Conditions for 250 Cases............................ 51 Figure 34: Average Area & Sine of R.A.A.N. vs. Time for 45? Inclination for Flat Plate with Random True Anomaly and Initial Attitude Conditions for 250 Cases.................... 52 ix Figure 35: Impacted Cylinder Case 1 vs. Time................................................................. 54 Figure 36: Impacted Cylinder Case 2 vs. Time................................................................. 55 Figure 37: Impacted Cylinder Case 3 vs. Time................................................................. 55 Figure 38: Impacted Cylinder Case 4 vs. Time................................................................. 56 Figure 39: Impacted Cylinder Case 5 vs. Time................................................................. 56 Figure 40: MEKF Results for 10 Days with Synthetic Ballistic Coefficient.................... 62 Figure 41: MEKF Results for One Half Year with Convergent Ballistic Coefficient ...... 63 Figure 42: MEKF Results for One Half Year with Dynamic Ballistic Coefficient .......... 63 Figure B.1: Euler Angles vs. Time ................................................................................... 71 Figure B.2: Rotational Rates vs. Time.............................................................................. 72 Figure B.3: Angular Momentum vs. Time........................................................................ 73 Figure B.4: Instantaneous Area vs. Time.......................................................................... 73 1 Chapter 1 Introduction 1.1 Background Since mankind started launching spacecraft into geocentric orbits in the 1950's, the threat of orbital collisions has been of increasing prevalence and importance. Increased understanding in the fields of geopotential and atmospheric density have led to highly accurate orbital propagation schemes. Although these schemes take into account a great many perturbations, they typically do not account for the variability of the area used in the calculation of atmospheric drag. This area is the projection of the surface area of the object onto a plane normal to the relative velocity of the object with respect to the atmosphere. Herein, it is called the "projected area" or "effective area." Unless a satellite is spherical or tidally locked, this area can play a major role in changing the magnitude of the drag force. An example of such change occurred in late September, 2011, when NASA's Upper Atmospheric Research Satellite plummeted back to Earth over the pacific ocean. Scientists and engineers tracking the satellite had difficulty predicting re-entry locations because the satellite began tumbling, so that the area used in determining the magnitude of drag was highly variable. This problem extends beyond re-entry events to objects in the upper reaches of Earth's atmosphere. Anytime the projected area is not known, accurate drag calculations cannot be performed. Historically, the calculation of the variable projected area has been left out of propagation schemes due to the complex time varying nature of the variables necessary to compute it. With modern computational 2 power, it is feasible to use high fidelity dynamic models to compute variables that were once deemed constant. 1.2 Focus of this Research To create a better understanding of the projected area, a digital simulation was developed that includes coupling of the rotational or "tumbling" motion of a satellite with the orbital motion of its center of mass. The focus of this work is the variability of the projected area and its effects on the ballistic coefficient. Since the other parameters in the calculation of the ballistic coefficient are constant, the difference between the two parameters is simply a ratio. Characteristic trends of the ballistic coefficient are shown that demonstrate a relationship between the ballistic coefficient's time evolving behavior and the orientation of the orbit it follows. Analyses are also conducted that provide insight into the emergence of certain characteristics, such as frequency and amplitude changes, that are observed in actual ballistic coefficient data. An orbit determination program is employed to show that with a known atmospheric density, estimates of the average ballistic coefficient with an averaged projected area can be obtained. The analysis and results in this thesis should prove useful in not only determining the ballistic coefficient of satellites with unknown attitude conditions, but also in predicting what the ballistic coefficient should be in the near future. 3 Chapter 2 Review of Literature and Low Earth Orbit 2.1 Space Debris The term space debris is restricted to manmade objects that are typically divided into three size categories: "large" objects greater than 10 cm diameter, "risk" objects between 1 and 10 cm diameter, and "small" objects less than 1 cm diameter. While all three sizes of debris can cause catastrophic damage to active orbital payloads, it is only feasible to catalog large and some risk sized objects. Objects in the risk category are typically the most dangerous, as they are too small to track but large enough to cause substantial damage. It is estimated that over 19,000 large objects, 500,000 risk objects, and upwards of tens of millions of small objects are currently in orbit about the Earth [1]. Sources for debris can include derelict spacecraft, depleted upper-stages of launch vehicles, and ejecta from orbital collisions. It is clear that without mitigation, the number of debris objects, and the threat to operating spacecraft will continue to increase. Smaller objects may be harmless, but the relative speed in a collision of a small object and a spacecraft in orbit can be nearly 15 km/s and lead to enormous releases of energy. The largest concentrations of orbital debris lie in low Earth orbit, around 200 to 2,000 km altitude, and in geostationary orbit, around 35,000 km altitude. Nearly all debris can be attributed to Russia, China, and the United States [2]. Due to the threat of debris colliding with operating spacecraft and the manned International Space Station, it is necessary that debris objects are tracked accurately and their positions catalogued. Several U.S. governmental and military agencies track more than 20,000 of the larger objects. They also perform analyses and predict close misses 4 and collisions between objects. Of course, active payloads must be monitored closely to avoid collisions with objects in the debris field, but it is also helpful to know of collisions between two debris objects as well. Simulations can be used to predict where the new debris from these collisions will reside as well as the number of objects created. 2.2 The Space Environment The motion of any object in geocentric orbit may be described using well known physical laws. If the object's physical characteristics and those of its environment are known, the study of its motion can be broken down into linear and angular displacements caused by forces and moments, respectively. It is necessary to understand the phenomena that causes these actions in order to model them mathematically with an acceptable degree of accuracy. 2.2.1 Forces Generally, the largest magnitude environmental force exerted on an object in orbit is gravity. Gravity is an attractive force between any two objects in the universe that have mass. As such, all objects in the universe that have mass attract the object under inspection at any given time. For orbital velocities, Newton's law of gravitation for two point masses can be used to model spacecraft motion. This is an acceptable way to achieve accurate results if the masses, velocities, and energies of the objects in the system are sufficiently small. 5 The principal aerodynamic force experienced by objects in low Earth orbit is drag. Strictly dissipative, drag is responsible for circularizing such an object's orbit and eventually causing it to reenter the Earth's atmosphere. For decades, perhaps the greatest challenge to be overcome in determining an accurate model for orbital motion has been that of accurately modeling the upper atmosphere. Atmospheric densities vary not only with altitude, but with the complex solar cycle. Solar output, the intensity of the radiation emitted by the sun, varies in an approximate 11 year cycle of the number and locations of sunspots. Many atmospheric density models, such as those in the Jacchia family (Jacchia 70, Jacchia 77) include models of the solar cycle. However, a substantial amount of error in density models for periods of time several days in length can still be experienced. The drag coefficient multiplied by the projected area and divided by the mass of the object form the ballistic coefficient. Along with density models, the ballistic coefficient has been a value long sought after by investigators to increase the accuracy of atmospheric drag models. Historically, accurate estimates of the drag coefficient have been obtained in laboratory tests for atmospheric vehicles. However, accurate recreation of the environment of space in a laboratory is much more difficult and expensive, leaving the drag coefficient to be estimated in other ways. Values of CD are typically between 2 and 3 when the projected area of an object is used as the reference area. The most common value for the drag coefficient used in orbital analyses is 2.2. This value was recommended by Cook as acceptable for satellites of unknown geometries [4]. Assuming that the mass of the object is known, the remaining parameter in the calculation of drag is the projected area, A. For most satellites, space stations, and payloads the geometries are well known. 6 Hence, for any given orientation of such an object with respect to the velocity, an A can be calculated. For an operational spacecraft, the challenge in determining the projected area in these cases is determining the orientation. However, due to the uncertain histories of many objects in the debris field, accurate geometry and orbital orientation are generally unknown, making the estimation of A very difficult. The last force that is considered of importance in this investigation is that due to the oblateness of the Earth. Due to the Earth not being entirely spherical and having a rough surface with large amounts of mass displacements, i.e. mountains, oceans, continents; its gravitational field is not uniform. The conventional approach to modeling the Earth's gravitational field is to use Legendre polynomials. The J2 zonal harmonic, which is symmetric and invariant about the Earth's axis of rotation, is used to model the Earth's oblateness, or bulging about the equator. The force due to this bulging results in periodic changes in the eccentricities and semi-major axes of the orbits spacecraft and debris, as well as secular changes in the right ascension of the ascending node, known as nodal precession or regression. Depending on an orbit's inclination and shape, the ascending node can precess or regress at rates greater than 6 degrees per day. 2.2.2 Torques There are several environmental factors that can cause changes in the orientation of an object. The two that have the greatest effect on an object in low Earth orbit are the gravity-gradient and aerodynamic torques. Gravity-gradient torque, like gravitational perturbations, is due to the non-uniformity of the Earth's gravitational field and the finite size of an orbiting object. The gradient is due primarily to the decreasing magnitude of 7 Earth's gravitational field with distance, which results in a greater force on portions of an orbiting object that are closer to the center of the Earth. The generation of a torque due to this effect of course requires that the object be non-symmetrical about at least one axis. The direction in which the torque acts on an object depends entirely on its current orientation. In Chapter 3 it will be shown that the magnitude of the gravity-gradient torque is inversely proportional to the cube of the magnitude of the object's geometric position vector [3]. For orbits below about 400 km altitude, aerodynamic torque becomes the primary perturbation affecting rotational motion [5]. Atmospheric densities at these altitudes allow for interactions to be handled by using a free-molecule flow model. Free- molecular flow, derived from the kinetic theory of gases, is based on the fact that in near- Earth space, the molecular mean free path is much larger than any particular object's physical dimensions. This treatment allows for incoming and outgoing particles to be handled separately. As such, two limiting cases arise in the study of the molecular reflections: specular and diffuse. Specular reflections dictate that the molecule bounces off a surface without any energy loss, and that it leaves the surface in the same plane it arrived in with equal speed. The momentum exchanged is in the surface-normal direction, orthogonal to the surface plane. In the case of diffuse reflections, molecules adhere to the surface upon collision. Then, they depart the surface at a later time with an energy that is determined by the surfaces temperature. Few molecules exhibit either model exactly, but fall somewhere in between [6]. The calculation of aerodynamic torques on orbiting objects is difficult for complex object geometries, but formulas are available for simple geometries like cylinders, cones, and flat plates. 8 Chapter 3 Development of the Analysis Program 3.1 Introduction The numerical analysis of low Earth orbit debris presented in this thesis was performed by a program written in Fortran 77. This program is a six degree-of-freedom (6DoF) digital simulation that numerically integrates a set of ordinary differential equations of motion. At user specified time increments, information such as position, attitude, and angular momentum can be extracted from the simulation and recorded. 3.2 Equations of Motion Herein, orbiting objects are modeled as rigid bodies. Hence, the equations of motion that determine the linear and angular motion experienced by an orbiting object as functions of time are well known. The equations describe the translational motion of the object's center of mass and the rotation of its principal axis system. The integration scheme is designed to integrate first-order, ordinary differential equations. Thus, the equations are written to describe changes in velocity (acceleration) and changes in position (velocity) for both translational and rotational motion. The phrase "six-degrees- of-freedom" stems from three possible linear movements and three possible rotational movements. The equations that describe these movements are derived about two sets of orthogonal, ordered triplet axes. Six first-order differential equations are required to describe the point mass motion of the object in one set of axes. Typically, an additional six equations are required to describe the rotational motion of the object in the other set 9 of axes. Conventionally Euler angles are used to define orientation of the rigid bodies. To overcome singularities in the rotational motion equations, the three equations that describe the rotational velocity are replaced by four equations for Euler parameters. All tolled, there are thirteen equations of motion: six for linear motion and seven for rotational motion. 3.2.1 Translational Equations The application of Newton's laws of motion requires the use of an inertial coordinate system. Although an Earth-centered coordinate system is not truly inertial, accelerations due to forces neglected in considering the Earth and an orbiting object as a system are small in magnitude compared to those due to the interaction between the Earth and the object. Thus, to study the relative motion of a spacecraft or debris object, an Earth-centered coordinate system is appropriate. The Earth-centered system can also be fixed in orientation. So, the six translational equations are derived using an Earth- centered inertial (ECI) coordinate system. An explanation of why the ECI system is appropriate to use is given in Appendix A. 10 Figure 1: Earth-centered inertial coordinate system with orbital elements Figure 2: Orbital plane with semi-major axis, a, occupied focus, F, and perigee, P Figures 1 and 2 show five orbital elements and the Cartesian coordinate system OXYZ used as the reference system for the equations of motion. Elements shown are the inclination, i, the right-ascension of the ascending node, ? or R.A.A.N., argument of the perigee, ?, the true anomaly, ?, and the semi-major axis, a. The sixth orbital element, 11 eccentricity, or simply e, falls between 0 and 1 and is a measure of how much the orbit deviates from a circle. Together they form the shape of an orbit and the position of the object along the orbit. The position is also given by the coordinates X, Y, and Z or the position vector r(X,Y,Z). The first three equations for linear motion are rather trivial. The time rate of change of a position component is equal to the velocity corresponding to that component, vx, y, z. The three equations which govern the acceleration of the center of mass of the object are broken down into four parts: two-body gravitational force due to the Earth, gravitational perturbations due to the Earth, aerodynamic drag, and forces due to other sources, such as the Sun and moon. otherdragoblategrav FFFFam rrrrr +++= 3.4 The forces acting on the object are labeled according to their source. The forces are caused by gravity, grav, the oblateness of the Earth, oblate, atmospheric drag, drag, and third-body objects, other. The sum of the forces is equal to the mass, m, times the acceleration of the object, ar , according to Newton's second law. The equations of motion for the acceleration caused by two-body gravitational force are derived from Newton's law of universal gravitation. 3.5 Z Y X vZ vY vX = = = & & & 3.1 3.2 3.3 12 The force gravFr is a function of the universal gravitational constant, G, the primary and secondary masses, m1 and m2, respectively, and the relative position vector of the two masses, rr . Consider an inertial coordinate system that contains two point masses. Figure 3: Inertial coordinate system with two bodies The origin, I, is a fixed point. 1R and 2R represent the inertial position vectors of the two point masses, while cmR is the inertial position vector of the system's center of mass. Taking into account mutual gravitational attraction only, the inertial equations of motion for the position vectors of the masses are given by: rrmmGRm ?2 2111 =&& 3.6 3.7 13 Here, the subscripts 1 and 2 refer to the two point masses, m1 and m2. The vector pointing from m1 to m2 is rv , and r? is the magnitude of that vector. Through vector addition, it is clear that: 21 RrR =+ v 12 RRr ?= v 3.8 3.9 1R r and 2R r are the inertial position vectors of the two point masses. Taking the second derivative of equation 3.9 yields the acceleration of the relative position vector of m2 with respect to m1. ( ) r r mmGr r mGr r mGr RRr ??? 2 212221 12 +?=??= ?= &&r &&&&&&r 3.10 3.11 We now assume that m1 is the mass of the Earth and m2 is the mass of the satellite. Therefore the mass of m2 is negligible in comparison to m1 and the following approximation is made: ( ) ?=??+ 121 mGmmG 3.12 This product will be referred to as the standard gravitational parameter, or simply ?. The standard gravitational parameter used in all simulations has a value of 398600.4415 km3/s2. We can now state the equation for the relative motion of the satellite when it and the Earth are both considered to attract like point masses: 3.13 To obtain a more general result, it is convenient to assume that the center of mass of the Earth-object system is O in the ECI system and is "inertial." It is also assumed that the 14 Earth is a rigid body with a gravitational potential U. For two-body gravity, part of U is U2B so that: BUr 2?= &&r where [ ] 21222 222 2 ??? ZYXr Z U Y U X UU BBB B ++= ? ?+ ? ?+ ? ?=? r KJI and rU B ?= 2 3.14 3.15 3.16 For the case of a more general gravitational field, we consider U to be a more general function of the position vector r. Here, the del symbol represents the mathematical gradient operator. In Cartesian coordinates, the acceleration components are the components of BU2? . [ ] ( ) 3 2 3 2 3 232222 2 2 1 r ZvZ Z U r YvY Y U r XXZYXvX X U Z B Y B X B ? ? ?? ?=== ? ? ?=== ? ? ?=++?=== ? ? ? &&& &&& &&& 3.17 3.18 3.19 Since the Earth is non-spherical, terms must be added to the right-hand side of the above potential to accurately model the gravitational field. It is convenient to express these terms in spherical coordinates (r, ?, ?), due to the nature of planets being nearly spherical. U for the Earth can now be expressed in the form [7]: 15 ( )??? ,,rBrU += where B is: 3.20 ( ) ( ) ( ) ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? +? ? ?? ? ? ?? ? ?? ? ? ?= ? ? ? = = 2 1 0 sinsincos sin ),,( n n m nmnmnm n e nn n e PmSmCrR PJrR rrB ??? ?? ?? 3.21 Figure 4: ECI coordinate system with spherical coordinates This full form version of the approximation of Earth's gravitational field consists of Legendre polynomials, Pnm, and the radius of Earth, Re. Zonal harmonic coefficients are Jn, Cnm, and Snm are tesseral harmonic coefficients for n ? m and sectoral harmonic coefficients for n = m. It should be noted that the cases for n = 0 is simply the two-body gravitation and n = 1 is not present because the origin coincides with the center of Earth's mass. Tesseral and sectoral harmonics are generally smaller than zonal harmonics and not nearly as important for accurate simulations of low Earth orbits. As such, the potential function can be approximated as: 3.22 16 From Fig. 4, it can be seen that rZ=?sin . The second zonal harmonic is of concern for this program, so for n = 2 the Legendre polynomial is: ( ) 212321sin23sin 2 2 0,2 ??? ?? ? ???= r ZP ?? 3.23 It follows that the second term in the above approximation for the gravitational potential, U2, is: ??? ? ??? ? ??= ??? ? ??? ? ? ?????????????? 35 22 2 2 2 2 2 13 22 1 2 3 rr ZRJ r ZJ r R rU ee ?? 3.24 Taking the gradient gives the component accelerations of the J2 zonal harmonic to be added to the two-body gravitational accelerations. ?? ? ?? ? ?= ?? ? ?? ? ?= ?? ? ?? ? ?= 57 32 2 ,2 57 22 2 ,2 57 22 2 ,2 915 2 315 2 315 2 r Z r ZRJa r Y r YZRJa r X r XZRJa e zJ e yJ e xJ ? ? ? 3.25 3.26 3.27 This completes the second phase of summing the forces that affect low Earth orbits in the analysis program. The next accelerations to be added to the object's equations of motion are those caused by aerodynamic drag. Using a similar approach to two-body gravitation, the acceleration experienced by objects in the upper atmosphere is given by: 3.28 The atmospheric density is found using an exponential model that reads in an altitude and outputs a density. This model is formed by a combination of the U.S. 17 Standard Atmosphere (1976), CIRA-72, and CIRA-72 with exospheric temperature with the temperature set at 1000 Kelvin [7]. The velocity term, Vrel, represents the relative velocity between the object and the atmosphere. This is calculated using the Earth's rotational velocity, e?r . k?ee erel rVV ?? ? = ??= r rrrr 3.29 3.30 Another modification to the drag equations is the use of the previously defined ballistic coefficient, ? . Taking the partial derivatives of the total accelerations gives the component accelerations to be added into the linear equations of motion. ( ) ( )[ ] ( ) ( ) ( )[ ] ( ) ( ) ( )[ ] ( )ZZeYeXZdrag eYZeYeXYdrag eXZeYeXXdrag vvXvYva XvvXvYva YvvXvYva ?+?++?= ??+?++?= +?+?++?= 21222 , 21222 , 21222 , 2 1 2 1 2 1 ???? ????? ????? 3.31 3.32 3.33 Adding these accelerations according to their components results in the final three equations for linear motion. These first-order ordinary differential equations are to be numerically integrated to propagate the location of the center of mass of objects in low Earth orbit. YdragYJY XdragXJX aarYv aar Xv ,,23 ,,23 ++?= ++?= ? ? & & 3.34 3.35 3.36 18 3.2.2 Rotational Equations As mentioned previously, there are seven equations of motion that are used to describe the rotational motion of an object in this investigation. For convenience, these equations are derived using the center of mass of the object and assuming that the body- fixed axis system CxByBzB is a principal coordinate system of the object. Figure 5: Body centered coordinate system This coordinate system coincides with the Earth-centered inertial system when the attitude angles are all zero. In a manner similar to the derivation of the translational equations of motion, the first equations considered are those that define the angular position variables. First, Euler angles are considered. These angles are functions of themselves and the angular velocity of the object. Euler angles are used to define the initial angular orientation of the object. The set of three angles (?, ?, ?) that are used define a 3-2-1 transformation from the ECI system to the principal CxByBzB system. 19 Figure 6: Rotation of body centered coordinate system through Euler angles That is, the z-axis rotation takes place first through an angle of ?, followed by a rotation about the new y-axis through ?, and finally a rotation about the new x-axis through ?. These angles are converted into a quaternion, the elements of which are Euler parameters, for integration purposes. The advantage of using quaternion algebra in attitude dynamics is that there are no singularities in the differential equation for the Euler parameters. Quaternion algebra was developed in 1843 by Sir William Rowan Hamilton [8]. A quaternion consists of four terms, i.e. a vector and scalar part. The vector part can be directly correlated to the axis that a given rotation is about, while the scalar part is the cosine of one-half the angle swept through in the rotation. The quaternion used in this analysis is a unit quaternion whose values range from -1 to 1. The initial quaternion for integration is created through the use of a direction cosine matrix, or DCM. A DCM is a matrix that houses the cosines of the angles between a vector and the coordinate system. For example, it can transform the components of a vector in the inertial frame to a body- fixed frame. The transpose of a DCM transforms that vector back to the original 20 coordinate frame. The following DCM is created using the three Euler angles that describe the attitude of the principle axis system. As denoted by the subscript, it transforms vectors in the inertial frame to the body frame. For brevity, cosine has been abbreviated as "c" and sine as "s" in the following expression. ? ? ? ? ? ? ? ? ? ? ?+ +? ? = )()()()()()()()()()()()( )()()()()()()()()()()()( )()()()()( ???????????? ???????????? ????? ccscsscsccss scsssccscssc scscc BIC 3.37 The components of this DCM are used to create the components of the unit quaternion that is to be integrated. The components of the unit quaternion are Euler parameters as defined in equations 3.38 through 3.41. The subscripts of Cij denote the location of the element in the DCM, where i is the row and j is the column. 4 2112 3 4 1331 2 4 3223 1 3322114 4 4 4 121 q CCq q CCq q CCq CCCq ?= ?= ?= +++= 3.38 3.39 3.40 3.41 The differential equation for a unit quaternion is a function of angular velocity and the current quaternion. This is where most of the advantages can be realized. Using the tilde notation, a matrix of angular velocities is formed as such: ? 3.42 The differential equation is given by [8]: 21 2 1=q& ?q 3.43 where ? ? ? ? ?? ? ? ?= 0 ~ T? ??v r 3.44 Equation 3.43 provides four scalar differential equations that propagate the components of a quaternion, q. ( ) ( ) ( ) ( )3322114 3412213 2413312 1423321 2 1 2 1 2 1 2 1 ??? ??? ??? ??? qqqq qqqq qqqq qqqq ???= +?= ++?= +?= & & & & 3.45 3.46 3.47 3.48 After integration, the Euler parameters are transformed into Euler angles and recorded as an output. Transformation formulas may be obtained by equating components of the DCM in terms of quaternion components with that of the DCM in terms of the Euler angles. The quaternion form of the DCM is: ( ) ( ) ( ) ( ) ( ) ( ) ?? ? ? ? ? ? ? ? ? ??+?+ +??+? ?+??+ = 2 2 2 1 2 3 2 441324231 4132 2 3 2 1 2 2 2 44321 42314321 2 3 2 2 2 1 2 4 22 22 22 qqqqqqqqqqqq qqqqqqqqqqqq qqqqqqqqqqqq BIC 3.49 The Euler angles are found by relating selected elements of each DCM. ( ) ( ) ??? ? ??? ? ??+ +=? ??+ +== ? 2222 43211 2222 432112 2tan2)tan( qqqq qqqq qqqq qqqq C C ?? 3.50 3.51 3.52 22 The last three equations that describe rotational motion are for the propagation of the rotational angular momentum vector. This angular momentum is the product of the object's inertia dyadic, I, with its angular velocity. ?rr I=H 3.53 It is a vector quantity that is conserved if there are no external torques on the object. Since there are torques acting the objects in this simulation, the angular momentum can vary greatly in magnitude and direction depending on the case. Since the principal axis system is rotating with respect to the inertial frame, the time rate of change of angular momentum is equal to the torque about the center of mass [8]. THH rr&r +?= ?~ 3.54 In component form, Eq. 3.54 gives the final three differential equations for the rotational motion of the object, zyxxyz yxzzxy xzyyzx THHH THHH THHH +?= +?= +?= ?? ?? ?? & & & 3.55 3.56 3.57 where Tx, Ty, and Tz are body-fixed components of Tr . After propagating the angular momentum, the new angular velocities are formed by multiplying the inverse of the inertia matrix, I rr , times the new angular momentum vector. ? ? ? ? =? ? ? ? ? xx H ? ? 3.58 The final step in completing the equations of motion is to determine expressions for the torques acting on an object. The gravity-gradient torque stems from the fact that the force due to gravity is non-uniform over a body. The derivation begins by using 23 Newton's law of universal gravitation in mass element form to write the element forces, iFd v , as a function of mass elements dm i. For the case of an arbitrary origin of the body axes, Fig. 7 is used in the derivation. i i i i dmr rFd 3 rv ? ?= 3.59 The torque on an object is due to this force acting on an element of mass at a distance di from the geometric center. The following equation takes advantage of ?r , the vector from the geometric center to the center of mass, and idr? , the vector from the center of mass to the mass element. iiiii FddFddTd rrrrrr ??+=?= )( ? 3.60 \ Figure 7: Coordinate System for Gravity-Gradient Torque The differential torque in Eq. 3.60 may be integrated over the entire object to obtain: 3.61 The position of the mass element in the inertial frame can be expressed in terms of the position of the center of mass. This is important because torques act about the center of mass. 24 ' iii drdrr ++=+= ? rrrrr 3.62 For cases involving space debris, the dimensions of the object are much smaller than the distances between the gravitational sources, i.e. 'ii dr +>> ?rr . This allows the use of an approximation of the distance between the mass element and the inertial origin. ( ) ( ) ( ) ( ) ??? ? ??? ? +??? ?? ??? ?? ??? ? ? ? ? ? ? ? ? + ++?+=?= ? ? ?? 2 ' 3 23 2 2' 2 ' 2233 3121 r drr r d r drrrrr iii iii rrrrrrrr rr ??? 3.63 Substituting the above equation into the mass element integral yields the gravity-gradient torque for the coordinate system in Fig. 7. ( ) ( )( ) iiiGG dmbrbrrbrmT ??3? 32 ??+?= ? rrrr ??? 3.64 To further simplify the expression, the first term can be eliminated by choosing the geometric center of objects to coincide with their center of mass. The unit vector b? is the orientation of the body axes with respect to the orbiting frame. This vector can be determined by taking the first column of the product of two direction cosine matrices. T OIBIBO CCC = 3.65 The integral in Eq. 3.64 may be evaluated and expressed in terms of the moment of inertia tensor to yield the gravity-gradient torque [7]. ( )[ ]bbrTGG ??3 3 ??= I?r 3.66 The Cartesian components of the above torque are the expressions for gravity-gradient torque that are part of the components Tx, Ty, and Tz for Eqs. 3.55 through 3.57. 25 ( ) ( ) ( ) ?? ? ? ? ? ? ? ? ? ? ? ? = 21 31 32 3 3 bbII bbII bbII rT XXYY ZZXX YYZZ GG ?r 3.67 The aerodynamic torque model adopted for this effort takes advantage of the aerodynamic force per unit mass calculated in the equations of motion. Multiplying the accelerations by the mass gives the aerodynamic forces in the inertial coordinate system. The forces must then be transferred into the body centered coordinate system by use of a direction cosine matrix. If we let CPrr be the vector pointing from the center of mass to the user-defined center of pressure, the cross-product of it and the aerodynamic forces results in the aerodynamic torques. [ ] [ ] ? ? ? ? ? ? ? ? ? ? ??? ??? ??? ?=?= ?= = xdragCPyydragCPx zdragCPxxdragCPz ydragCPzzdragCPy dragCPaero zdragydragxdragdrag CPzCPyCPxCP arar arar arar mFrT aaamF rrrr ,, ,, ,, ,,, rr rr rr rrr rrrr r 3.70 The resultant torque added to the rate of change of the angular momentum is now the summation of the gravity gradient and aerodynamic torques. zaerozGGz yaeroyGGy xaeroxGGx TTT TTT TTT ,, ,, ,, += += += 3.71 3.72 3.73 3.3 Effective Area To determine the effective area of a tumbling object in orbit, it is convenient to take advantage of the simple shapes that make up the greater object. In the case of a 26 right-circular cylinder, it becomes a circle when viewed from the x-direction, and a rectangle when viewed from the y or z-directions. For a right-circular cone, these become a circle and triangle, respectively. The cosine of the angle between the normal direction of these shapes and the relative wind determine the magnitude of the contribution they make to the total area. Since the normal direction from the circular base of the cylinder or cone are in the x-direction, the dot product between an x-direction unit vector and the relative wind give the cosine of the angle between the two vectors. [ ] V V TTBI rr 001cos ??= C? where [ ]zyx VVVV =r V CVCVCV zyx r 131211cos ++=? 3.74 3.75 The cosine of the angle, ?cos , is simply a function of the inertial velocity and three direction cosines. It can now be applied to the individual areas to determine the entire total effective area. 2 cos cos1 ?? ??+?= sidebase AAA 3.76 If the angle between the normal direction of the base of the object reaches 90 degrees, the cosine of that angle will be zero and the resultant total area will be the area of the side of the object. Likewise, any deviation of the angle from zero degrees will add to the total area via the side without regard to the orientation. This allows the use of the Pythagorean identity as a means of finding the multiplier for the contribution of the side, Aside. 27 3.4 Numerical Integration Propagation of the equations of motion through time is achieved by numerical integration. The algorithm used is a modified version of the Runge-Kutta-Fehlberg fifth order solution method. It is a six-step, variable time-step algorithm that products a fourth and fifth order solution. The difference of these solutions is compared to a tolerance which determines if the time-step will increase or decrease. When the comparison is below the desired accuracy tolerance, the fifth-order solution is reported and the algorithm continues, with an increase in the time step. This algorithm proves to be time efficient while maintaining accuracy [8]. Depending on the initial conditions, mainly the rotational velocities, one low Earth orbit simulation takes a fraction of a second to run. Modifications include logic to set a maximum number of iterations and a maximum time step. 3.4.1 Integration Performance To determine the rate at which integration error contaminates the results, stability tests must be performed on the analysis program. To have constant orbital elements, angular momentum, and energy, a number of subroutines are withheld from the program's normal operation. In effect, it reduces the program to a two-body propagation scheme with an object tumbling at a constant rate. An arbitrary integration time of two years was chosen because it is much longer than any proposed ballistic coefficient or position forecasting window. The orbit propagated is a near-circular, 700 km altitude orbit with the initial right ascension of the ascending node at 0 degrees and inclination of 45 degrees. Three attributes, which would be constant in an ideal scenario, were chosen 28 as indicators of integration accuracy decay: angular momentum magnitude, H, specific orbital energy, ?, and rotational kinetic energy, Trot. All three indicators increased secularly throughout the integration, but at very slow rates. Translational kinetic energy, Ttrans, and potential energy, Utrans, varied periodically as expected for an elliptic orbit. The specific orbital energy is a combination of the translational and rotational energies. rottranstrans TUT ++=? 3.77 The angular momentum magnitude increased by 4.3E-3 J-s, an increase of 0.0135 percent. Rotational kinetic energy increased by 1.5E-4 J, an increase of 0.0257 percent. The specific orbital energy increased in accordance with the increase in rotational kinetic energy. The following plot show these performance parameters plotted against time. Figure 8: Performance Indicators vs. Integration Time The next plot shows the energy components versus time. It is worth noting that the translational kinetic and potential energies counter-balance each other and account for a constant total translational energy throughout the integration. The orbital elements, except for the true anomaly, remained nearly constant as in an ideal scenario. 29 Figure 9: Spacecraft Energies vs. Time 30 Chapter 4 Analysis Theory and Results 4.1 Analysis Approach The principal objective of the this thesis is to obtain a greater understanding of the factors that cause time variations of the ballistic coefficients of debris objects. Of particular interest are the causes of another projected area variations. Another objective is to investigate the feasibility of using orbit determination techniques to detect rapid changes in the ballistic coefficient. A possible reason for sudden changes in a ballistic coefficient a collision of the tracked object of debris with another that results in a change of the angular momentum of the tracked object. Another reason is that there are un- modeled density variations that change the drag on the object, but the changes are attributed to the ballistic coefficient. The projected area of a tracked object is an important variable in the drag model, but often it is overlooked and assumed to be constant. A better understanding of variables such as the projected area will lead to more accurate force and prediction models that can be used to predict the state of an object for use in evasive maneuvers. Sibert, et al. concluded that accurate position calculations as well as an increase in state vector update frequency can result in a nearly complete reduction of near collision events [9]. Due to the number of pieces of debris and the limited information available about each individual piece, certain assumptions must be made to obtain typical results. Here, it is assumed that all tracked objects in the simulation program have a mass of 100 31 kilograms. Three geometries are used to represent the objects: a cylinder, a right circular cone, and a flat plate. The cylinder is representative of a satellite bus or spent launch vehicle stage. Right circular cones are also common in rocketry, especially for re-entry vehicles. A flat plat geometry is representative of a separated solar panel or a portion of a spacecraft. In all simulations, the three geometries have lengths extending in the body-x direction of 10 meters. The cylinder and cone have base diameters of 1 meter, while the flat-plate has a width of 1 meter. It is assumed that the flat-plate is infinitely thin. The three geometries are meant to give diversity to the simulations, that is, of course still less than the diverse set of geometries of actual space debris. 4.1.1 Average Effective Area Since the simulation program integrates non-linear equations of motion with high accuracy, a very small time step is generally required. It is not uncommon for this time- step to be less than one second, meaning there are thousands of opportunities to record the dynamical states and effective area per orbit. When collecting data on the dynamical states, this can be very advantageous and result in high fidelity area data as shown in Fig. 10. The orbit propagated to generate data for Fig. 10 is nearly circular with the object at an altitude of 700 km. The orbit is inclined at 45 degrees with initial right ascension of ascending node, true anomaly, and argument of perigee of zero degrees. 32 Figure 10: Effective Area vs. Time Over 1 Orbit for Cylinder However, the effective area may vary so much over an extended period of time that without manipulation, finding trends can be difficult. This problem is alleviated by calculating the average area over a time period, and reporting that value less frequently. The average effective area, A , is calculated by summing the area at each time, A(t), times the time-step, t? , and dividing the sum by the difference of the total observational time as indicated in Eq. 4.1. ( ) initialfinal t tt tt ttA A final initial ? ?? = ? = 4.1 33 Figure 11: Average Effective Area Overlaid on Instantaneous Effective Area for Cylinder The observational time used to obtain the average effective area in Fig. 11 is 8 hours. Clearly, an average effective area is a much more representative factor that scales the amount of drag on an object is observed. After propagating the motion for ten days using both the actual and average effectives areas, the final positions are compared. Table 1: Difference in Positions Determined Using Instantaneous Area vs. Average Area x y z rrelative Act. Area (km) 5922.5495 -1246.1790 3658.2549 - Ave. Area (km) 5922.5400 -1246.1057 3658.2954 - Difference (m) 9.1268699998 -71.233669999 -39.30757999 81.869538554 34 Linear interpolation was used between average area data points to obtain the instantaneous area for continuous drag calculations. Since the magnitude of the position error is only 81.87 meters, using and studying the average effective area can be considered reasonable for some purposes. For collision avoidance the actual area is probably better if characterized well. 4.1.2 Existing Ballistic Coefficient Data The Air Force Space Command Space Analysis Center is one of the governmental agencies tasked with tracking debris and performing conjunctional analyses on the catalog of objects including active payloads. Researchers at the center have developed a method to extract the ballistic coefficient estimates for the objects they track. Although the Air Force's software and methods are not public domain, they have provided a sample of ballistic coefficient data sets for 17 pieces of debris. The samples are for the period January, 2005 to January, 2007 and are presented as percent deviations from the mean for the two years. The following figures are presented without edit. The title contains the NORAD catalog number and name of the object. Also stated in the title is the term "DCA Values," where DCA stands for dynamic calibration atmosphere. This refers to the Air Force's method of calibrating the temperature values for their synthetic atmosphere by observing the forces on satellites whose characteristics are well established. 35 00461 THOR DEBRIS DB/Bave 2005-6 DCA Values -100 -80 -60 -40 -20 0 20 40 60 80 100 105.0 105.5 106.0 106.5 107.0 YEAR DB /B Av e % Figure 12: Ballistic Coefficient vs. Time for Satellite 00461 12184 DELTA DEBRIS DB/Bave 2005-6 DCA Values -60 -40 -20 0 20 40 60 80 100 DB /B Av e % Figure 13: Ballistic Coefficient vs. Time for Satellite 12184 36 10230 DELTA DEBRIS DB/Bave 2005-6 DCA Values -100 -80 -60 -40 -20 0 20 40 60 80 100 105.0 105.5 106.0 106.5 107.0 YEAR DB /B Av e % Figure 14: Ballistic Coefficient vs. Time for Satellite 10230 20857 CZ-4 DEBRIS DB/Bave 2005-6 DCA Values -60 -40 -20 0 20 40 60 80 100 DB /B Av e % Figure 15: Ballistic Coefficient vs. Time for Satellite 20857 37 These four examples were chosen because they easily exhibit the extreme variability of the ballistic coefficient as well as some insight into the nature. The data in Fig. 12 shows well-behaved periodicity for the first year, followed by a large increase in the mean value of the ballistic coefficient in November of the first year. That increase is followed by high frequency variation at what appears to be a different mean value. The data in Fig. 13 exhibits an almost periodic relationship that looks like it should be predictable. The data in Fig. 14 shows a much more complex time variation that looks to be somewhat predictable. However, the data in Fig. 15 shows extreme jumps in value with noisy data throughout. Some of the ballistic coefficient data exhibits almost period variations with periods on the order of months. Interest was found in these variations because the period of the nodal regression of an object orbiting the Earth varies from 1 to 8 months depending on the orbit's inclination. Thus, an investigation into the variation of the effective area of an object and its relation to inclination and the regression of the ascending node was conducted. Of interest in the following section are the inclinations of the orbits of the objects that were used to calculate ballistic coefficient plots. Table 2: Inclinations of Example Satellites [10] Satellite 00461 Satellite 12184 Satellite 10230 Satellite 20857 Inclination 67.5086? 98.7720? 29.0228? 98.6656? 4.2 Relationship of Effective Area Variations with Inclination The right ascension of the ascending node of a prograde orbit about the Earth varies periodically and secularly. The secular change is a rotation of the orbital plane 38 about the Z-axis of the ECI coordinate system and is known as nodal regression (or precession for retrograde orbits). The cause of the regression is the non-sphericity of the Earth and in turn its gravitational field. The rate at which an orbit regresses depends principally on the semi-major axis, eccentricity of the orbit, and its inclination. The average precession rate of the right ascension of the ascending node is given by Eq. 4.2 [8]. ( )i ae rJ e cos )1(2 3 2722 2 2 ? ? ? ? ? ? ? ? ? ?=? ?& 4.2 Since the orbits simulated in this thesis are nearly circular, inclination is the main factor determining the period of the motion of the right ascension of the ascending node in the results presented here. The "nodal rate" is generally between 1.5 and 8 degrees per day, resulting in periods ranging from 240 to 45 days, respectively. Figures 16 through 30 contain plots of time histories of the average effective area and the sine of the right ascension of the ascending node. In each case, the initial orbit of the object is nearly circular with the object initially placed at 700 kilometers altitude and at the ascending node. For calculating the semi-major axis in all simulations, the radius of the Earth used is 6378.137 km. The initial true anomaly, right ascension of ascending node, and argument of perigee are all zero. The initial attitude with respect to the inertial frame and angular rates are zero. After the propagation begins, environmental torque cause the object to tumble in a "random" manner. All simulations were propagated for 365 days with various initial inclinations. The total propagation time was chosen to exhibit the periodicity of the effective area and its relationship to the nodal regression. 39 During a year the right ascension of the ascending completes at least one period for all inclinations. The average of the effective area was taken every 8 hours. The sine of the right ascension of the ascending node is presented with a negative phase shift of 45 degrees. This phase shift aligns the effective area and the sine of the ascending node to show their similar periodicities. Figures 16 through 34 contain the phase shifted sine of R.A.A.N. 40 Figure 16: Average Area & Sine of R.A.A.N. vs. Time for 10? Inclination for Cylinder Figure 17: Average Area & Sine of R.A.A.N. vs. Time for 30? Inclination for Cylinder 41 Figure 18: Average Area & Sine of R.A.A.N. vs. Time for 45? Inclination for Cylinder Figure 19: Average Area & Sine of R.A.A.N. vs. Time for 60? Inclination for Cylinder 42 Figure 20: Average Area & Sine of R.A.A.N. vs. Time for 80? Inclination for Cylinder As can be seen, the period of the average effective area is clearly related to the period of the right ascension of the ascending node. Also, the amplitude grows with the increases in inclination. These factors should prove useful in determining a bandwidth in which the ballistic coefficient should reside. Of course, the patterns may change with varying initial conditions. The following figures contain area data for a right circular cone with the same initial conditions. Since the center of pressure of the cone is not at the same location as its center of mass, aerodynamic torques play a role in causing rotational motion of the cone. 43 Figure 21: Average Area & Sine of R.A.A.N. vs. Time for 10? Inclination for Cone Figure 22: Average Area & Sine of R.A.A.N. vs. Time for 30? Inclination for Cone 44 Figure 23: Average Area & Sine of R.A.A.N. vs. Time for 45? Inclination for Cone Figure 24: Average Area & Sine of R.A.A.N. vs. Time for 60? Inclination for Cone 45 Figure 25: Average Area & Sine of R.A.A.N. vs. Time for 80? Inclination for Cone The behavior of the variable effective area of the right circular cone is markedly similar to that of the cylinder. The same symmetries in moments of inertia are present in both the cone and the cylinder, but only half of the side area is present in the cone. Again, the amplitudes of the effective areas increased with the inclination. The aerodynamic torque seems to have played a minor role in the rotational motion at this altitude. The next set of plots is generated from the same simulations using a flat plate as the object. A flat plate has the same moment of inertia symmetries as the other geometries, so motion is expected to be similar. The main difference occurs in the area calculation, which is now a function of only one surface. 46 Figure 26: Average Area & Sine of R.A.A.N. vs. Time for 10? Inclination for Flat Plate Figure 27: Average Area & Sine of R.A.A.N. vs. Time for 30? Inclination for Flat Plate 47 Figure 28: Average Area & Sine of R.A.A.N. vs. Time for 45? Inclination for Flat Plate Figure 29: Average Area & Sine of R.A.A.N. vs. Time for 60? Inclination for Flat Plate 48 Figure 30: Average Area & Sine of R.A.A.N. vs. Time for 80? Inclination for Flat Plate The plots also exhibit periodic trends that follow the phase shifted sine of the ascending node. Although the trends are similar to those of the cylinder and right circular cone, the peaks corresponding to the minimum of the sine of R.A.A.N. are far less amplified. Since the geometry is so different from a cylinder or cone, it is expected to see different characteristics. However, seeing the same qualitative trends corresponding to the ascending node is surprising and encouraging to anyone seeking the ballistic coefficient. This also bodes well for the application of this work to the LEO debris field since the geometries are unknown and almost certainly quite varied. Trends emerged between the average effective area and the ascending node for a number of different inclinations with the initial attitude conditions set to zero for all cases. To determine if these trends holds for other conditions, further testing is required. In the following simulation, a uniform pseudo-random number generator was employed 49 to randomize the initial conditions of the attitude of the object before integration began. As before, the object was placed at the ascending node. The year-long simulation was repeated 250 times, each with unique initial attitude conditions. Figure 31: Average Area & Sine of R.A.A.N. vs. Time for 45? Inclination for Cylinder with Random Initial Attitude Conditions for 250 Cases Obviously, the same trend of following the phase shifted sine of the ascending node occurs. The area falls within a bandwidth that is largest when the phase shifted sine of R.A.A.N. reaches its lowest value, corresponding to a -45 degree right ascension of the ascending node. A similar simulation was conducted to determine if the initial position of the satellite plays a role in the behavior of the average effective area. This time, the true anomaly was randomized. 50 Figure 32: Average Area & Sine of R.A.A.N. vs. Time for 45? Inclination for Cylinder with Random True Anomaly for 250 Cases As can be seen in Fig. 32, the previously stated trend holds. The characteristics of the plot are very similar to Fig. 31, which strengthens the argument that the effective average area, and in turn the ballistic coefficient, has the same frequency as the sine of R.A.A.N. In another benchmark of validity, agreements with the actual ballistic coefficient data can be seen. For a year's worth of data from satellite 12184 in Fig. 13, the ballistic coefficient peaked three times. Satellite 12184's orbit is inclined 98.772?, or nearly 10? from polar. This is the same case as an 80? orbit, which for all three geometries the average area peaked three times as well. The same conclusion can be drawn between satellite 00461 in Fig. 12 and the characteristics of the 60? effective average area plots. 51 For Figs. 33 and 34, simulations were run for the right circular cone and the flat plate. The propagation time was one half year, repeated 250 times. The true anomaly and initial attitude conditions were selected from a pseudo-random uniform distribution ranging from -180? to 180?. Figure 33: Average Area & Sine of R.A.A.N. vs. Time for 45? Inclination for Cone with Random True Anomaly and Initial Attitude Conditions for 250 Cases 52 Figure 34: Average Area & Sine of R.A.A.N. vs. Time for 45? Inclination for Flat Plate with Random True Anomaly and Initial Attitude Conditions for 250 Cases It appears that no matter the geometry, there is a substantial dependence on the sine of the phase shifted right ascension of the ascending node. Like the cylinder, the average area falls within a bandwidth for altering initial conditions of the cone and flat plate. 4.3 Impact Analysis One of the possible reasons for a sudden change in the ballistic coefficient seen in Figs. 12 and 15 is an impact with another object. Relative speeds for objects in low Earth orbit can approach 15 km/s, resulting in large amounts of kinetic energy in the smallest of objects. Of interest to this investigation is the change in angular momentum due to the impact. Another way in which a change in the ballistic coefficient could occur is the 53 addition, or reduction, of mass by an inelastic collision. Not only would this change the ballistic coefficient inherently, but could quite possibly affect the rotational behavior by altering the moments of inertia. The impact as simulated herein is based on the assumption that a small orbiting mass is absorbed by the object under inspection. The total mass change of the system is negligible, and the impact takes place a certain distance from the center of mass of the object, causing a change in angular momentum. The object used for this analysis is the cylinder. The angular momentum added to the system is formed by taking the cross product of the impact vector and the linear momentum exchanged between the impacting particle and the cylinder. relpimpactimpact VmrH rrr ?= 4.3 In the above equation, relVr is the relative velocity between the cylinder and the impacting particle, impactrr is the vector pointing from the center of mass to the position of impact, mp is the mass of the impacting particle, and impactHr is the angular momentum added to the object. The mass of the impacting particle was chosen to be 1 gram, and the impact location was 4 meters from the center of mass. These parameters were held constant while the impacting particle's orbit was chosen using the pseudo-random number generator. Impacting particle velocities ranged from 0 to 100 percent of the object's velocity since they have the same eccentricity and altitude. The propagation for each impacted cylinder case lasted one half year. The resulting average area is plotted over the average area for the same case without an impact, the control case. The cylinder's orbit was inclined 45? for each case with a nearly circular shape and an altitude of 700 km. 54 The angular orbital elements are initialized at 0 degrees, and the impact occurred after one day of integration. The cases shown are for different completions of the algorithm in which the relative velocity is unique in each case. These cases demonstrate the similarities of the effective area's time-varying behavior between unique impacts, as well as the differences that occur. Figure 35: Impacted Cylinder Case 1 vs. Time 55 Figure 36: Impacted Cylinder Case 2 vs. Time Figure 37: Impacted Cylinder Case 3 vs. Time 56 Figure 38: Impacted Cylinder Case 4 vs. Time Figure 39: Impacted Cylinder Case 5 vs. Time 57 It is clear that a substantial reduction in the average effective area occurs after impact for each case. No cases arose where the orbits of the object and impacting particle were the same, resulting in zero relative velocity. The periodic trend occurs at the same frequency, but the amplitude range is much larger. The trend in Fig. 37 is of utmost importance, as it shows that in impact is a plausible cause for the erratic behavior of the ballistic coefficient exhibited by satellite 00461 in the second year of Fig. 12. Also, the behavior of the area in Fig. 36 is very similar to that of satellite 10320 in Fig. 14. This indicates that satellite 10320 was tumbling with substantial angular momentum during the two-year period since there are no qualitative changes in its behavior. 58 Chapter 5 Orbit Determination 5.1 Methodology In this chapter, orbit determination techniques are shown to prove useful in determining the average effective area, which further validates the use of an average area for tracking purposes. Also, simulations are produced to test the orbit determination program's ability to detect large characteristic changes in the ballistic coefficient. The process outlined here makes use of a modified extended Kalman filter (MEKF) written in Fortran 77. The program tracks and updates seven states: the inertial position vector, the inertial velocity vector, and the ballistic coefficient. The MEKF is a well-known tool used for orbit and trajectory determination. The particular algorithm used was taken from Tapley [11]. Since the MEKF estimates and update the position to the most accurate location, it does not rely on a variable ballistic coefficient to calculate the correct drag accelerations to propagate the orbital position accurately. In other words, the MEKF is a tool used in this work to estimate values of the ballistic coefficient, not estimate accurate positions of the object. The ballistic coefficient that the MEKF determines is a best fit to the data it is given. It will be shown later that this "best fit" is very close to the ballistic coefficient calculated using the average effective area. 5.1.1 Data Generation Program The MEKF essentially operates by propagating the state vector, reading an observation, and applying updates to the state after every observation. This event driven 59 process continues, for simulation purposes, until there are no more observations to read. Thus, observations must be generated for the program to read. The data generation program uses the 6DoF simulation to provide the motion of the cylinder for a given time period. Every 15 seconds, a range measurement is taken from the nearest radar station on Earth. This only occurs if the cylinder is in the line of sight to the radar station, which is restricted to 60? downward from zenith in any direction. The location of these radar stations are generated randomly and stored before the process begins. The 10 radar stations modeled are comparable in number to the stations in the space surveillance network used by the U.S. government [12]. The average effective area is also recorded during this process for later comparison with the estimated ballistic coefficient. This area is simply multiplied by the ratio of the drag coefficient to the cylinder mass to yield the average ballistic coefficient. The orbit used for data generation was nearly circular with an initial altitude of 700 km. The orbit was inclined 45 degrees with all other angular orbital elements initialized at zero degrees. 5.1.2 Power-Density Matrix and Process Noise One drawback that plagues the MEKF is a tendency to diverge, or produce incorrect state updates, after many observations. This occurs because the program is based on a variable covariance matrix. After time has passed and many observations have been processed, this covariance matrix becomes small. The result is the MEKF having more confidence in the current state than the observations. To avoid this problem, a diagonal matrix the same size as the covariance matrix is added to the differential equation governing the covariance matrix's rate of change. 60 QPAAPP ++= T& 5.1 In Eq. 5.1, P is the covariance matrix, A is the state relation matrix, and Q is the power- density matrix. This diagonal matrix is known as the power-density matrix, Q, described by Vallado as the second moment (or covariance) of the process noise [13].The process noise represents un-modeled accelerations acting on the object. It is the fact that process noise is present in the model that allows the MEKF to update the ballistic coefficient. The ballistic coefficient does not change through propagation, but its value is estimated using the observed range to the range calculated by the model used to propagate the six translational states. The addition of the power-density matrix to the covariance's differential equation has a large effect when the covariance becomes small, which is when it is most necessary to take action. The effect is an increase in the elements of the covariance matrix. This process is repeated every time the covariance matrix becomes small, resulting in non-divergence of the MEKF. The values that occupy the power-density matrix were determined by a binary- encoded genetic algorithm made available at Auburn University. The optimization program's goal is to minimize the root-mean-square error of the difference between the range measurements and the ranges calculated in the filter. The design variables that the genetic algorithm varied to meet this goal are the 7 diagonal elements of the power- density matrix. The reference orbit, which was nearly circular with an initial object altitude at 700 km, was propagated for 50 days. The algorithm evaluated 200 members for 20 generations. Therefore, the values for the power-density matrix found by the genetic algorithm are the best found for the 4000 cases evaluated. The last diagonal position, which is coupled to the ballistic coefficient, can be left at zero. This allows the 61 estimated ballistic coefficient to converge to a value, and not vary every time the covariance matrix was affected by the power-density matrix. These two options will be referred to as a convergent ballistic coefficient option and a dynamic ballistic coefficient option. The diagonals of the power-density matrix determined by the genetic algorithm are shown in Table 5.1. Table 5.1: Diagonal Values of Power-Density Matrix (Q) Q11 Q22 Q33 Q44 0.99139344E-02 0.20386984E-03 0.20386984E-03 0.67142769E-09 Q55 Q66 Q77 0.10986999E-08 0.14038943E-07 0.97215438E-05 5.2 Orbit Determination Results The first simulation used a synthetic ballistic coefficient that was not calculated by the 6DoF. This was done to show that the MEKF does an exemplary job of finding the mean of the ballistic coefficient. The orbit used to generate the observations was inclined at 45 degrees. The power-density option of a convergent ballistic coefficient was used. The output of the orbit determination program consists of two plots: the average and estimated ballistic coefficients versus time and the measurement residuals (O-C) versus time. The O-C plot shows the difference between the range measurement calculated by the MEKF and the range measurement given by the observation at that time. If the filter is working correctly, these values are typically less than 10 meters and only influenced by the random noise introduced in the data generation program. 62 Figure 40: MEKF Results for 10 Days with Synthetic Ballistic Coefficient The synthetic ballistic coefficient consisted of a mean value of 0.2 m2/kg plus a sine function with a period of 10 days and amplitude of 0.04 m2/kg. Figure 40 shows that the MEKF converged on the mean value of 0.2 m2/kg and that the filter did not diverge as shown by the O-C plot. The MEKF converged to a near perfect answer within 2 days. A similar simulation was repeated twice with actual ballistic coefficient data. The simulations lasted one-half year, with the only difference being the power-density matrix option. 63 Figure 41: MEKF Results for One Half Year with Convergent Ballistic Coefficient Figure 42: MEKF Results for One Half Year with Dynamic Ballistic Coefficient 64 As can be seen from Figs. 41 and 42, the MEKF does an excellent job of estimating the ballistic coefficient with either power-density matrix option. In Fig. 41, the estimated ballistic coefficient converged to the mean value of the sinusoidal average ballistic coefficient. Figure 42 exhibits the variability introduced by using the dynamic ballistic coefficient option. This option is not ideal for instantaneous retrieval of the actual ballistic coefficient, but should prove useful in determining time varying trends that the actual ballistic coefficient exhibits. 65 Chapter 6 Conclusions An investigation of the effects of gravitational and aerodynamic forces and torques on space debris was conducted to gain a better understanding of the time variations of the ballistic coefficient of pieces of debris. A six-degree-of-freedom, rigid body model digital simulation was developed as well as an orbit determination program that uses a modified extended Kalman filter. The simulation program is capable of calculating the effective drag area at every time step of integration and provides accurate drag calculations. These calculations are important for the determination of conjunctional events and re-entry of objects into the Earth's atmosphere. Results obtained show that with knowledge of the atmospheric density, the ballistic coefficient can be calculated for a number of different geometries. The attitude of a tumbling object, and thus the behavior of the ballistic coefficient, exhibit certain observed characteristics that are reproducible in the digital simulation. These characteristics that are dependent on the orbital inclination and the regression (or precession) of the right ascension of the ascending node should prove useful in predicting the ballistic coefficient not only for a 10-day predication window, but perhaps for much longer. It has also been shown that using an average effective drag area is an accurate substitute for the instantaneous effective area. The instantaneous effective area of a tumbling object yields high amplitude and frequency that requires more computation. The orbit determination program that uses a modified extended Kalman filter was employed to show that determining and tracking the ballistic coefficient of a tumbling object is feasible. Although this process is too computationally intensive to apply to the 66 entire catalog of trackable objects, it could be applied to objects that are "at risk" of entering a conjunctional event. The orbit determination program yielded both the average ballistic coefficient and data that represented the trend of the coefficient. There are number of ideas that could increase the fidelity of the analysis program and this work if further research on the subject is pursued. The effective area can be calculated numerically, instead of analytically, to increase the accuracy and allow for more complex geometries. This can be coupled with "shadowing" algorithms that take into account portions of the object that are blocked by other portions such as solar panels. These ideas can also be applied to determine accurate aerodynamic torques for complex geometries. A continuance of the current effort should be to determine the underlying cause of the periodicity of the effective area and its relationship to its orbital inclination. Another objective of such work should be to determine the frequency relationship between the right ascension of the ascending node and the average effective area. 67 Bibliography [1] NASA Orbital Debris Program Office, "Orbital Debris Frequently Asked Questions", July 2009, http://orbitaldebris.jsc.nasa.gov/faqs.html [retrieved 2 March 2012]. [2] NASA Orbital Debris Program Office, "Orbital Debris Quarterly News", Volume 12, Issue 1, January 2008, p. 3, http://orbitaldebris.jsc.nasa.gov/newsletter/pdfs/ODQNv12i1.pdf [retrieved 23 January 2012]. [3] Zanardi, M., Real, F., "Environmental Torques Acting on a Low Earth Orbiter Cylindrical Spacecraft," Advances in Space Research, Vol. 31, No. 8, pp. 1981- 1986, 2003. [4] Cook, G. E., ?Drag coefficients of Spherical Satellites,? Ann. Geophys., Vol. 22, p. 53, 1966. [5] Beletskii, V. V., "Motion of a Artificial Satellite About Its Center of Mass," NASA TT F-429, 1965. [6] Hughes, P. C., Spacecraft Attitude Dynamics, John Wiley & Sons, Inc., New York, NY, 1986. [7] Wertz, J. R., Spacecraft Attitude Determination and Control, Vol. 73, D. Reidel Publishing Company, Dordrecht, Holland, 1978. [8] Curtis, H. D., Orbital Mechanics for Engineering Students, 2nd ed., Elsevier, Burlington, MA, 2010. [9] Sibert, D., Borgeson, Maj. D., Peterson, G., Jenkin, A., Sorge, M., "Operational Impact of Improved Space Tracking on Collision Avoidance in the Future of LEO Space Debris Environment," Advanced Maui Optical and Space Surveillance Technologies Conference, Maui, HI, Sept., 2010. [10] Satellite Tracking, Prediction and Informations about Objects in the Sky, http://www.infosatellites.com/ [retrieved 29 March 2012]. [11] Tapley, B. D., Schutz, B. E., Born, G. H., Statistical Orbit Determination, Elsevier, Burlington, MA, 2004. [12] United States Strategic Command Factsheet, http://www.stratcom.mil/factsheets/USSTRATCOM_Space_Control_and_Space_ Surveillance/ [retrieved 7 April 2012]. 68 [13] Vallado, D. A., Fundamentals of Astrodynamics and Applications, 3rd ed., Springer, New York, NY, 2007. [14] Etkin, B., Reid, L., Dynamics of Flight: Stability and Control, 3rd ed., John Wiley & Sons, Inc., Hoboken, NJ, 1996. 69 Appendix A ECI Coordinate System and Alternate Relative Motion Equations Derivation The ECI coordinate system is used to write the equations of motion for the spacecraft or object in Chapter 2. The ECI system has its origin at the center of the Earth, which is also assumed to be the geometric center of the Earth. Despite its name, the origin of the ECI system is not inertial. Its axes are oriented so that the Z-axis points towards the celestial north pole and its X-axis is pointed toward the vernal equinox. Let 1rr and 2rr be position vectors of the center of mass of the Earth, O in Fig. 1, and the center of mass of the object, r(X,Y,Z) in Fig. 1, with respect to a true inertial point. Whether the Earth and object are modeled as point masses or rigid bodies, their total accelerations are modeled as other other FFrm FFrm ,21/222 ,12/111 rr &&r rr&&r += += A.1 A.2 where m1 and m2 are the masses of the Earth and object, respectively. Also, 2/1Fr is the force on the Earth due to the object, 1/2Fr is the force on the object due to the Earth, otherF ,1 r is the remaining force on the Earth and otherF ,2 r has the same meaning for the object. The forces otherF ,1r and otherF ,2r are due primarily to the Sun, Moon, and Jupiter. The relative acceleration ra &&rr = of the object with respect to the center of the Earth is A.3 By Newton's Third Law of Motion, the force on the Earth due to the object is equal and opposite to . Hence, 70 otherother FmFmFmma ,1 1 ,2 2 1/2 12 1111 rrrr ?+? ? ? ?? ? += A.4 The force 1/2Fr includes the force on the object due to the Earth's gravitational field and the force on the object due to the Earth's atmosphere, i.e., the atmospheric drag on the object. The forces otherF ,1r and otherF ,2r , divided by the masses of the Earth and the object, respectively, are essentially the same because of the close proximity of the object to the Earth. Thus, for this analysis, the relative acceleration of the object with respect to the Earth is approximated by: 1/2 12 11 F mma rr ?? ? ?? ? += A.5 Since the force 1/2Fr is the sum of the gravitational force and the atmospheric drag, ?? ? ?? ? ??? ? ? ?? ? += relrel D grav VV ACF mma rrr ? 2 11 12 A.6 where gravFr is the gravitational force, ? is the atmospheric density, A is the "flat-plate" or effective area for calculating the drag, CD is the drag coefficient, relVr is the velocity of the center of mass of the object with respect to the atmosphere and Vrel is the magnitude of relVr . Since the mass of the Earth is approximately 5.974 ? 1024 kilograms [8], 111 mmm ??? ? ?? ? + A.7 Finally, with m = m2, A.8 71 Appendix B Attitude Dynamics The one purpose of the analysis program developed for this investigation and described in Chapter 3 is to integrate the equations of motion for the attitude of an object in orbit. A closer look at the variables involved in this integration shows that, as expected, Euler angles and rotational rates are continuous. The angular momentum is determined about the body principal axis system and is related to the rotational rates by the principal moments of inertia. The orbit for the following simulation is inclined at 45 degrees and the body is placed at the ascending node with zero initial attitude conditions. Figure B.1: Euler Angles vs. Time The Euler angles shown in Fig. B.1 were obtained from the DCM constructed using the attitude quaternion. Since this transformation uses inverse trigonometric identities, the 72 range for the Euler angles is between -180 and 180 degrees, but the angle ? continuously increases. Figure B.2: Rotational Rates vs. Time As shown in Fig. B.2, the angular velocity components ?y and ?z vary with time due to the gravity-gradient torque. However, ?x remains zero because the torque about the body x-axis, the symmetry axis of the cylinder, is zero. The Euler angles vary according to the following kinematical relationships [14]. ( ) )sin()cos( )tan()cos()sin( ????? ??????? zyx ?= ++= & & B.1 B.2 B.3 Clearly, the object is tumbling slowly, as the rates do not exceed 0.15 degrees per second. 73 Figure B.3: Angular Momentum vs. Time Figure B.4: Instantaneous Area vs. Time 74 Continuity in the area algorithm can be seen by running the simulation without torques. Figure B.4 shows the reference area for the cylinder with no torques, propagated over one period. Without torques, the object's attitude is fixed in inertial space. As the object progresses through it's orbit, the instantaneous projected area changes smoothly and predictably. Areas corresponding to the end of the cylinder, less than 1 m2, and the side of the cylinder, 10 m2, are both perpendicular to the oncoming wind at different times throughout the orbit.