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.