Ph.D. Thesis:
The Dynamics of Tethers and Spacewebs
David J. McKenzie, BEng
Submitted in fulfilment of the requirements for the Degree of Doctor of Philosophy
Departments of Mechanical and Aeronautical Engineering,
Faculty of Engineering,
University of Glasgow, G12 8QQ, Scotland, UK
January 2010
.
David McKenzie, 20032010
Contents
Abstract 15
Acknowledgements 16
1 Introduction 19
1.1 Overviewof thesis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2 Literature review 23
2.1 Earlypioneers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.2 Tetherfundamentals ........................... 25
2.2.1 Fundamental tetherconcepts .................. 26
2.2.2 Nonrotating tethers ....................... 26
2.2.3 Momentumexchange ....................... 27
2.2.4 Motorized momentum exchange tethers . . . . . . . . . . . . . 29
2.2.5 Staging with momentumexchange . . . . . . . . . . . . . . . 30
2.2.6 Deploying and recovery of tethers . . . . . . . . . . . . . . . . 31
2.2.7 Captureand rendezvousof tethers . . . . . . . . . . . . . . . 33
2.2.8 WSB trajectories ......................... 33
2.3 Tetherderived structures. . . . . . . . . . . . . . . . . . . . . . . . . 35
2
2.3.1 Deployed spacestructures .................... 35
2.3.2 Spacewebs ............................ 36
3 Tether modelling 38
3.1 Constructing the equations of motion of a system . . . . . . . . . . . 38
3.1.1 Validating codemodules ..................... 39
3.2 Centreof massmodelling . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.2.1 Masspoints&body selection .................. 40
3.2.2 Centreof masscalculations ................... 41
3.2.3 Positionsof masspoints ..................... 43
3.3 Constructing Lagrange’sequations . . . . . . . . . . . . . . . . . . . 48
3.3.1 Energy modelling ......................... 49
3.3.2 Choosing generalised coordinate systems . . . . . . . . . . . . 50
3.4 Nonconservativeforces. . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.5 Rotationsystems ............................. 52
3.5.1 Rotationsummary ........................ 53
3.5.2 YZRotations ........................... 55
3.5.3 ZYRotations ........................... 56
3.5.4 Comparing the equations – position equations . . . . . . . . . 59
3.5.5 Comparing the equations – kinetic energies . . . . . . . . . . . 59
3.5.6 Comparing the equations – potential energies . . . . . . . . . 60
3.5.7 YZ rotationininertial axes ................... 61
3.5.8 Singularities in the rotational systems . . . . . . . . . . . . . . 61
3.5.9 Practical uses of the different equations . . . . . . . . . . . . . 62
3
3.6 Numericalintegrationtechniques .................... 62
4 Dynamics of tethers with inclination
64
4.1 Additionofinclinationtotethermodel . . . . . . . . . . . . . . . . . 65
4.1.1 5degreeoffreedommodel .................... 66
4.2 Constructing theequationsof motion . . . . . . . . . . . . . . . . . . 67
4.2.1 Rotations in the local coordinate system . . . . . . . . . . . . 68
4.2.2 Lagrange’sequations ....................... 71
4.2.3 Lagrange’sequations – local system . . . . . . . . . . . . . . . 71
4.2.4 Lagrange’s equations – global system . . . . . . . . . . . . . . 75
4.2.5 Nonconservativeforces. . . . . . . . . . . . . . . . . . . . . . 77
4.3 Analysisof tetheronaninclined orbit . . . . . . . . . . . . . . . . . 80
4.3.0.1 Behaviour of R, .
and .
with time .......... 81
4.3.0.2 Behaviour of a
with time ............... 81
4.3.0.3 Behaviour of .a
with time ............... 82
4.3.0.4
Behaviour of a
with .
. . . . . . . . . . . . . . . . . 82
.¨
4.3.0.5 Behaviour of ,
,
.
. . . . . . . . . . . . . . . . . 83
4.3.0.6 Behaviour of ,
,
.¨. . . . . . . . . . . . . . . . . . 83
4.3.0.7 Behaviour of V
with time .............. 84
4.3.0.8 Behaviour of tension with time . . . . . . . . . . . . 84
4.3.0.9 Behaviourof stresswith time . . . . . . . . . . . . . 84
4.3.1 Analysisof tensionand stressintether . . . . . . . . . . . . . 85
4.4 Tetherreleaseaiming
. . . . . . . . . . . . . . . . . . . . . . . . . . . 86
4.4.1 Tethererroranalysis ....................... 88
4
4.4.2 Lunarcapture. . . . . . . . . . . . . . . . . . . . . . . . . . . 90
4.5 A demonstration Lunar mission with a tether launch . . . . . . . . . 91
4.5.1 Spacecraftsizing ......................... 91
4.5.2 Analysis .............................. 93
4.6 Conclusions ................................ 94
4.7 Inclinationplots. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
5 Dynamics of tethers with length deployment 108
5.1 Deploymentand recovery withtheMMET . . . . . . . . . . . . . . .108
5.2 Addition of length as a generalised coordinate . . . . . . . . . . . . . 109
5.3 Constructing theequationsof motion . . . . . . . . . . . . . . . . . .110
5.4 Nonconservativeforces. . . . . . . . . . . . . . . . . . . . . . . . . .115
5.4.1 Tetherbraking ..........................116
5.5 Analysisoflengthdeploymentof tether . . . . . . . . . . . . . . . . .117
5.5.1 Discussion of results for deployment . . . . . . . . . . . . . . . 118
5.6 Refining theequationsof motion ....................120
5.7 Recovery of tethers ............................123
5.7.1 Recoveryprofile . . . . . . . . . . . . . . . . . . . . . . . . . .124
5.7.2 Theequationsof motionforrecovery . . . . . . . . . . . . . .125
5.7.3 Solving the equations of motion for recovery . . . . . . . . . . 126
5.8 Deployment and recovery on an elliptical orbit . . . . . . . . . . . . . 128
5.9 Centreof massmovement ........................130
5.10Conclusions ................................131
5.11Figures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .132
5
5.11.1 Linearkineticenergy .......................143
6 Dynamics of spacewebs 144
6.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .144
6.2 Modelling methodology . . . . . . . . . . . . . . . . . . . . . . . . . .146
6.3 Web meshing ...............................146
6.3.1 Web structure. . . . . . . . . . . . . . . . . . . . . . . . . . .146
6.3.2 Dividing theweb .........................147
6.3.3 Webdivisions ...........................148
6.4 Equationsof motion ...........................149
6.4.1 Centreof massmodelling . . . . . . . . . . . . . . . . . . . . .149
6.4.2 Rotations .............................150
6.4.3 Positionequations ........................152
6.5 Energy modelling .............................153
6.5.1 Virtual workterms ........................154
6.5.2 SimplifyingAssumptions .....................155
6.5.3 RobotPositionModelling ....................156
6.6 Spacewebdynamics ...........................157
6.6.1 Dynamical simulation. . . . . . . . . . . . . . . . . . . . . . .157
6.7 Investigating the stability of the spaceweb . . . . . . . . . . . . . . . 158
6.7.1 Massof theweb . . . . . . . . . . . . . . . . . . . . . . . . . .158
6.7.2 Robotmass ............................159
6.7.3 Web asymmetry . . . . . . . . . . . . . . . . . . . . . . . . . .160
6.7.4 Robotcrawl velocity .......................162
6
6.8 Statistical investigation into stability . . . . . . . . . . . . . . . . . . 162
6.8.1 Massesof theweb and satellite . . . . . . . . . . . . . . . . .163
6.8.2 Massesof theweb and robot. . . . . . . . . . . . . . . . . . .164
6.8.3 Massof theweb and angularvelocity . . . . . . . . . . . . . .165
6.8.4 Massof therobotand angularvelocity . . . . . . . . . . . . .167
6.9 Conclusionsand recommendations. . . . . . . . . . . . . . . . . . . .168
6.9.1 Conclusions ............................168
6.9.2 Recommendations. . . . . . . . . . . . . . . . . . . . . . . . .168
7 Conclusions 169
7.1 Tetherrotations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .169
7.2 Tetherswithinclination .........................169
7.3 Deploymentand recovery of tethers . . . . . . . . . . . . . . . . . . .170
7.4 Spacewebs ................................170
7.5 Furtherstudy ...............................171
8 Bibliography 172
Glossary 183
A Equations of motion with inclination 184
B Spacewebs – additional graphics 188
B.1 Case 1 – CoM plots of stability while increasing web mass . . . . . . 188
B.2 Case 1 – CoM plots of stability while increasing robot mass . . . . . . 190
B.3 CoMplotsof stability whilechanging robotpaths . . . . . . . . . . .191
7
B.4 Case 1 – CoM plots of stability while decreasing robot velocity . . . . 192
C Spacewebs – additional data 195
D Material strengths 200
E Motor properties 202
8
List of Figures
1.1
Diagramof theMMET. ......................... 20
2.1
Tether length increment L, from [Cosmo and Lorenzini, 1997, p175] 27
3.1
Symmetricaldumbbelltether ...................... 41
3.2
Symmetricaldumbbelltetheronacircularorbit(solid). Afterpayload
release,thefacility systempath atapogee(dash) andpayload
atperigee(dots) ............................. 42
3.3
Asymmetricaldumbbellininertial axes . . . . . . . . . . . . . . . . . 44
3.4
Local coordinateconstruction ...................... 45
3.5
Local coordinate construction after rotation through angle .
about
thefacility mass. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.6
Local coordinate construction after translation from the facility to
theCoM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.7
Local coordinate construction after translation from the CoM to the
inertial coordinatesystem ........................ 46
3.8
Local coordinate construction after rotation through .
about the inertial
coordinatesystem ......................... 47
9
3.9
Rotation sequence for YZ rotation: Rotating a
about Yaxis then
.
about Zaxis. a
is shown as a negative rotation here for ease of
visualisation. ............................... 57
3.10 Rotation sequence for ZY rotation.. Rotating .
about Zaxis then
a
about Yaxis. a
is shown as a negative rotation here for ease of
visualisation. ............................... 58
4.1
Ephemeris of Lunar orbital inclination [JPL, 2008] . . . . . . . . . . . 65
4.2
Rotation sequence for YZ rotation with inclination. a
is shown as a
negative rotation here for ease of visualisation. . . . . . . . . . . . . . 69
4.3
RotationsequenceforYZrotationwithinclination. a
and .
are shown
as negative rotations here for ease of visualisation. . . . . . . . . . . . 70
4.4
Eclipsecheckingfunction(nottoscale) . . . . . . . . . . . . . . . . . 79
4.5
Sample plot of the Eclipse function in Mathematica . . . . . . . . . . 79
4.6
Diagram showing an example of the WSB trajectory, from ESA Bulletin103:[
Biesbroek andJanin,2000]. . . . . . . . . . . . . . . . . . 87
4.7
Change in apogee of the payload after release from the tether system
after one half of an orbit compared to small changes of system
parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
4.8
Change in apogee and inclination of the payload after release from
the tether system after onehalf of an orbit compared to small changes
of systemparameters ........................... 96
4.9
The orbital debris environment from LEO to GEO, from [Wikipedia,
2008] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.10 Plots of R, .
and .
vs. timeforthefivecases. . . . . . . . . . . . . . 99
4.11 Plots of a
vs. timeforthefivecases. ..................100
10
4.12 Plots of .a
vs. timeforthefivecases. ..................101
4.13 Plots of a
vs. .
forthefivecases. ....................102
4.14 Phase plots of .
forthefivecases.. . . . . . . . . . . . . . . . . . . .103
4.15 Phase plots of a
forthefivecases.. . . . . . . . . . . . . . . . . . . .104
4.16 Plots of V
vs. timeforthefivecases. .................105
4.17PlotsofTensionvs. timeforthefivecases. . . . . . . . . . . . . . . .106
4.18PlotsofStressvs. timeforthefivecases. . . . . . . . . . . . . . . . .107
5.1 Discretised tether diagram with n
=5masspoints. . . . . . . . . . .112
5.2 InitialMMETSpinup ..........................132
5.3 InitialMMETSpinup ..........................133
5.4 InitialMMETSpinup ..........................133
5.5 Initial MMET Spinup, discretized tether mass model . . . . . . . . . 134
5.6 Initial MMET Spinup, discretized tether mass model . . . . . . . . . 135
5.7 Initial MMET Spinup, discretized tether mass model . . . . . . . . . 136
5.8 MMET reelin,firstphase ........................136
5.9 MMET reelin,firstphase ........................137
5.10MMET reelin,firstphase ........................138
5.11MMET reelin,secondphase. . . . . . . . . . . . . . . . . . . . . . .138
5.12MMET reelin,secondphase. . . . . . . . . . . . . . . . . . . . . . .139
5.13MMET reelin,secondphase. . . . . . . . . . . . . . . . . . . . . . .140
5.14MMETdeploymentwith elliptical orbit . . . . . . . . . . . . . . . . .140
5.15MMET recovery with elliptical orbit . . . . . . . . . . . . . . . . . .141
5.16 MMET deployment, centre of mass movements . . . . . . . . . . . . . 141
11
5.17 MMET recovery, centre of mass movements . . . . . . . . . . . . . . 142
5.18 MMET recovery, position and acceleration of centre of mass . . . . . 142
6.1 Concept rendering of the spaceweb with spiderbots on the surface. . 145
6.2 0divisions .................................148
6.3 1division .................................148
6.4 2divisions .................................148
6.5 3subspans – 0divs ...........................148
6.6 3subspans – 1div ............................148
6.7 3subspans – 2divs ...........................148
6.8 4subspans – 1div ............................148
6.9 5subspans – 2divs ...........................148
6.106 subspans – 3divs ...........................148
6.11 Simplified spaceweb layout with s
subspans. . . . . . . . . . . . . .149
6.12Midpointlocation. . . . . . . . . . . . . . . . . . . . . . . . . . . . .150
6.13 Spaceweb diagram with 3 subspans shown in inertial space. Note:
thelocal axesareintheYZplane. . . . . . . . . . . . . . . . . . . .152
6.14 Contours of CoM movement while varying web and satellite mass . . 164
6.15 Contours of CoM movement while varying web and robot mass . . . . 165
6.16 Contours of CoM movement while varying web mass and @.
. . . . .166
@t
6.17 Contours of CoM movement while varying robot mass and @.
. . . .167
@t
B.1 mW
eb
=27kg
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .188
B.2 mW
eb
=33.75kg
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .188
B.3 mW
eb
=37.8kg
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .189
12
B.4 mW
eb
=40.5kg
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .189
B.5 mW
eb
=44.82kg
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .189
B.6 mW
eb
=47.25kg
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .189
B.7 mW
eb
=54kg
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .189
B.8 mW
eb
=67.5kg
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .189
B.9 Mrobot
=1kg
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .190
B.10 Mrobot
=10kg
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .190
B.11 Mrobot
=15kg
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .190
B.12 Mrobot
=16kg
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .190
B.13 Mrobot
=20kg
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .190
B.14 Mrobot
=100kg
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .190
B.15Case1 ...................................191
B.16Case2 ...................................191
B.17Case3 ...................................191
B.18Case4 ...................................191
B.19Case5 ...................................191
B.20Case6 ...................................191
B.21 Vrobot
=10.0m/s
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .192
B.22 Vrobot
=2.0m/s
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .192
B.23 Vrobot
=1.0m/s
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .192
B.24 Vrobot
=0.2m/s
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .192
B.25 Vrobot
=0.1m/s
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .192
B.26 Case 1 – 3 robots moving symmetrically round the web perimeter . . 193
B.27 Case 6 – 1 robot moving asymmetrically along first subspan . . . . . 194
13
List of Tables
3.1 Orderof rotationsin3D and2D . . . . . . . . . . . . . . . . . . . . . 55
4.1 Reference table with figure and corresponding rotation . . . . . . . . 68
4.2 Initial conditions for the inclination numerical solutions. . . . . . . . 80
4.3 Initial conditions for the inclination numerical solutions. . . . . . . . 98
C.1
Influenceof web massonmaximumCoM . . . . . . . . . . . . . . . .196
C.2
Influence of symmetrical robot movement on maximum CoM . . . . . 196
C.3
Influence of simulation time(effectively robot speed across web) on
maximumCoM ..............................196
C.4
Influenceof robotmassonmaximumCoM . . . . . . . . . . . . . . .197
C.5
Statistical investigation into stability . . . . . . . . . . . . . . . . . . 198
C.6
Results of the statistical investigation into stability . . . . . . . . . . 199
D.1
Propertiesof selected materials . . . . . . . . . . . . . . . . . . . . .201
E.1
PropertiesofGEMotor5BC49JB1115 . . . . . . . . . . . . . . . . .202
14
Abstract
The thesis ‘The dynamics of tethers and spacewebs’ investigates the motion of the
Motorized Momentum Exchange Tether (MMET) on an inclined orbit, and while
deploying and retracting symmetricpayloads. TheMMET system isused asabasis
for examiningthe stability of spacewebs using atriangular structure oftethers while
rotating. The motion of small robots is introduced as they move on the spaceweb,
and their motions are found to influence the behaviour of the structure. Several
methodsof limiting thedestabilisinginfluencesof therobotsareconsidered,and are
found to stabilise the web in most circumstances.
Astructured methodfordescribing therotationsof atethersystem isoutlined,and
different rotational systems are compared. This lays the foundation for the further
chapters examining MMET dynamics on an inclined orbit and while deploying and
recovering the payloads. Lagrange’s equations are generated for the three cases,
and are solved using standard numerical integration techniques. To emphasise the
practical uses of the MMET system, several missions are analysed by using the
system as a reusable launcher for microsatellite payloads.
15
Acknowledgements
I would like to thank the EPSRC for providing a doctoral training grant, enabling
the research for this thesis.
I would like to thank the staff at the University of Glasgow for their thoughts and
discussions. Inparticular,DrMaxVasileforhissharp insights,DrDavidForehand
for discussions on many interesting problems and Dr Donald Ballance and Dr Ian
Watson for an introduction into teaching. Needless to say, I am also thankful for
thehard work of alltheadministrativestaff, inparticularMarilynDunlop andJane
Livingston.
I would like to thank Prof. Matthew Cartmell for the chance to study such an
interesting problem, for his kindness and patience, for his supervision, and for enthusing
me with thejoys of research.
I would like to thank Colin McInnes for encouraging me to study space systems
engineering. The love of space research, fostered under his guidance, played a big
part in my decision to undertake a Ph.D..
I would like to thank Chris Draper for the long and interesting conversations we
had while sharing an office, and for his enthusiasm in general.
Iwould liketothank thestaffatESAACTfortheopportunity tostudy spaceweb
dynamics.
I would like to thank my friends at pausegaming.com, for the camaraderie.
I would like to thank my families for their support. My dad and mum for encour
16
aging my talents inmathematics, engineering andLEGOfromavery earlyageright
through to the final chapter of this thesis.
I would like to thank my very good friends: Raj, Aaron, Ailsa, Dougie, George,
Rhys, and of course,Pob. Oursmallparty seemstohavegained afewmembers and
lost a few, but we seem to have weathered the Storm: good friends are as valuable
as gold.
I would like to thank Joanna, for her kindness and encouragement, for her love
and support, and for a thousand things too big and small to mention.
17
Declaration of Authorship
The author hereby declares that this thesis is a result of the work carried out in
the Department of Mechanical Engineering at the University of Glasgow during
the period between October 2003 and January 2010. This thesis is an original
contribution by the author unless otherwise indicated.
David McKenzie
January 2010
18
Chapter 1
Introduction
The field of space research is perhaps unique, in that many novel, different and
interesting concepts are commonly developed from a seed of an idea to a satellite
in a relatively short space of time. Additionally, mankind has been active in space
foravery shorttime,buthasmadeastoundingprogresses inunderstanding inboth
technical and social fields as a direct result of this.
The cost of these endeavours has been equally astounding, with many nations
contributing significant fractions of their GDP towards developing and maintaining
a space presence. Equally, all nations with a space presence continually reassess
theircapability –politicianswantdevelopmentsfaster,accountantswanteverything
donecheaper,scientistswanttodobetter–asadirect result,novelideascanflourish.
TheMotorizedMomentumExchangeTether(MMET) isoneofthese ideas,with
the potential to provide cheaper onorbit launch costs and a lighter mass when
compared to conventional rocket technologies.
A diagram of the MMET is shown in Figure 1.1. A payload mass is connected
to a large central facility mass with a high strength tether. A symmetrical mass
is attached to the opposite side of the facility and these two tethers and payloads
19
formthepayloadarm. Afurthertwotethersarejoinedtothefacility to increase
the inertia of facility, these are called the stator arm, and are similar in design to
the payload arm.
payload 2
facility
stator 1
stator 2
payload 1
Figure 1.1: Diagram of the MMET.
As with all alternative ideas, the MMET is best suited to a particular mission
profile. With the MMET, its strongest advantage is that it is reusable, and may be
employed to launch satellites and payloads as long as the tether is intact.
Unfortunately, the MMET does have drawbacks, the four tethers are susceptible
to severing from orbital debris, a symmetrical release strategy is required for each
payload release and the risk of interference between the payload and stator arms is
high. Additionally,themaximumorbital velocity thepayload canachieveis limited
by the strength of the tethers.
With these in mind, a large body of research into the performance of the MMET,
momentum exchange tethers, dumbbell tethers, and tethers in general has revealed
20
that there are definite advantages in employing these as a tool to enhance or enable
many satellite missions. As with all tools, spaceor Earthbased alike, it is
imperativeto usethe righttoolforthe rightjob.
The research into momentum exchange tethers has been focussed on constructing
mathematical models and analysing these models to determine the performance
and stability, to avoid any potential weaknesses and to gain a better insight into
the fundamental properties of the tether systems. Hardware testing is employed
whereverpractical(andfunds allow!) to validatethese models.
Several techniques are available to build these mathematical models. These include
Newtonian based techniques, using forces and accelerations, and energy based
techniques using Lagrangian formulations.
This thesis will develop a mathematical model of the MMET using a Lagrangian
based approach to further analyse the performance of the system in two key areas:
on an inclined orbit, and while deploying and recovering the tether.
A new structure – the spaceweb – is proposed and analysed, composed of a net
overlaid on a backbone of a rotating tether system. This can be employed as a base
on which to build large space structures such as solar sails and interferometers.
1.1 Overview of thesis
Chapter2 criticallyreviewsthecurrent literatureandgivesanoverviewof thewide
range of tether systems.
Chapter 3 introduces the rotational methods that underpin the rotational tether
system and rotational systems are examined and compared. A method of constructing
Lagrange’s equations for a rotating tether is outlined and the tether system is
mathematically constructed from its individual components.
Chapter 4 contains an analysis of the MMET dynamics on an inclined orbit, with
21
inplane and outofplane motions. A demonstration of the tether’s capabilities is
shown with an exemplar mission using aWeakStabilityBoundary(WSB)trajectory
from Medium Earth Orbit(MEO) to the Moon.
Chapter 5 studies the MMET while deploying and recovering the payload and
tethers. A mathematical model of the MMET system undergoing deployment is
derived, and the stability of the system is assessed using the movement of the centre
of mass.
Chapter 6 uses the fundamental tether models to construct a spaceweb – a structure
in space that may be used as a platform to assemble and deploy large space
structures.
22
Chapter 2
Literature review
The field of tether research is broad, and spans many branches of engineering,
physics and mathematics. Tethers may be broadly split into three categories: momentum
exchange tethers, electrodynamic tethers and structuresbuilt using tethers.
The notable events in tether history will be critically evaluated.
While this is not an exhaustive review of tethers, as provided in [Cartmell and
McKenzie, 2008] or [Cosmo and Lorenzini, 1997], it is intended to cover critically
the pertinent topics relevant to this thesis.
2.1 Early pioneers
In 1895, the father of astronautics, Konstantin Eduardovich Tsiolkovsky visited the
World’sFair inParis. SeeingtheEiffelTower,he imaginedavast structurestretching
from the tower to a celestial castle in geostationary orbit above the Earth. From
this ideaof atower,fundamentallyacompressivestructuredescribed in[Tsiolkovski,
1959],the ideaofaspaceelevatorwasborn. YuriArtsutanov,Tsiolkovsky’sstudent,
discussed lowering a structure from orbit to theground in[Artsutanov,1960].
23
Unfortunately, the good work in the USSR was not immediately available in the
West, and it was a number of years before the idea of largescale tethers and space
elevators was published in English. A brief, but informative, article published by
[Isaacs et al., 1966] appeared in Science,andprovided the impetusforacademics in
the UK and the Americas to study this new idea.
In one of the first articles in Acta Astronautica, [Pearson, 1975]identifies the main
problems(materials limits,vibrationalmodesand stability),and solvesagreat many
of them from first principles1, introducing tapered tethers and taper ratios. Much
of the terminology still used derives from that source, for example, characteristic
height of materials2 .
Hans Moravec metaphorically severs the link between the space elevator and the
Earth, postulating the tether as a momentum exchange device in [Moravec, 1977],
calling it the Skyhook3 . This giant spoke rotates, touching down once and interfacing
with the surface periodically. He identifies Mars as an ideal candidate for
the Skyhook, and provides a good grounding for the specific materials and optimal
layouts that are expanded in further research.
Several works of fiction have been inspired by these early pioneers, including
Clarke’s ‘The Fountains of Paradise’. [Clarke, 1979] brings a social aspect to the
space elevator research that is often overlooked, and is worth mentioning as the
source in which many young researchers first discover the field. Indeed, Clarke remarked
in his final interview, see [Das, 2008], that the space elevator will be built
‘about 10 years after everyone stops laughing’!
1Onthestateoftheartanaloguecomputer in1975,outpoweredby several ordersof magnitude
by today’s mobile phones!
2Pearson defines the characteristic height as ‘the height to which a constantdiameter tower
of the material could be built in a uniform oneg field without exceeding the stress limit of the
material at the base’.
3Moravec acknowledges in that paper that the Skyhook idea originated with John McCarthy
of Stanford.
24
2.2 Tether fundamentals
The foundation textbook, ‘Dynamics of Space Tether Systems’ by Beletsky and
Levin [Beletsky and Levin, 1993], rigorously introduces the dynamics of tethers, in
a progressive and pragmatic manner.
The book discusses the fundamental topics in tether research, such as material
density, strength, and orbital location, at a level accessible for a newcomer to tethers.
The book then develops equations of motion for a massless and massive flexible
tether with variations, covering perturbations and other environmental effects. The
dynamics are examined in terms of stability and oscillatory behaviour, based on a
Newtonianderivation. ElectroDynamic(ED)tethers, librationandrotation,deployment
and retrieval, andLunaranchored and satellitering systems servetocomplete
the coverage of the book. A very useful set of references is also provided, up to the
publication year of 1993.
Achapterof thetextbookbyMcInnes andCartmell[McInnes andCartmell,2006]
on the dynamics of propellantless propulsion systems covers both solar sails and
tethers, each sharing the common goal of overcoming the constraints of using a
reaction mass. This reviews missions involving tethers to date, then moves on to
summarisetheperformanceexpectationsofhanging, librating,and spinning tethers,
setting them in the context of results extant in the literature.
Carroll’s paper (see [Carroll, 1986]) gives an excellent background to the many
differenttypesoftethermissionproposedandflown, includingtheGemini12 mission
thatpioneered thetechnique ofgravitygradient stabilisation with theAgena vehicle.
This paper outlines the various roles of tethers, from the space elevator, librating
tethers, rotating tethers and momentum exchange tethers.
25
2.2.1 Fundamental tether concepts
Thereareafewfundamental conceptsthat limitthetether, itsperformanceand the
missions it may usefully provide.
Rotating tethers(whencomparedtogravitygradientstabilisedtethers) area logical
progression, and extend the usefulness of the tether system. Tether oscillations
and other nonlinear dynamics cause concern, and the everpresent possibility
of severing thetetherduetospacedebris isonlygoing toworsenwithbetteraccess
to space. Spinning the tether up to rotational speed is a notable problem; in the
vacuum of space, there is nothing to provide a base to rotate against.
A tether differs from a rocket engine, a massbased propulsion technique, in that
mass isnot ejected. Thisenablestetherstoseverthe linkbetweenpayload andfuel.
However tethers suffer from a few notable drawbacks: the tether line is vulnerable
todebrisdamage,tethersare limitedby theirfundamental materialproperties,and
rotating tethers must maximise their velocity increment.
Tethers can broadly be split into two categories: the rotating and nonrotating
tether. Therotatingtether(seeSection2.2.3)isprimarilyconcerned with imparting
a large velocity increment. The nonrotating tether has a variety of uses, with large
sections of the literature devoted to control and stabilisation of the tether, as is
covered in the next section.
2.2.2 Nonrotating tethers
Tethers mayprovide stabilisation, control ordamping to a structure or constellation
in space. Active stabilisation of the tether is usefully demonstrated by ED tethers
in a highly readable report, [Hoyt, 2002]. Pendulum librations, transverse wave
oscillations, and skiprope modes are all investigated and control laws proposed in
the form of feedback algorithms.
26
As Figure 2.1 (taken from [Cosmo and Lorenzini, 1997, p175]) illustrates, the
rotating tether provides a large advantage over the nonrotating tether in terms of
the velocity increment to the payload. A frequently cited result is that the apogee
gains 7L
in height from a hanging tether of length L, less than 14L
from a swinging
tether and more than 14L
from a rotating tether.
Figure 2.1: Tether length increment L, from [Cosmo and Lorenzini, 1997, p175]
2.2.3 Momentum exchange
Momentum exchange is, in general, the main method in which a tether can impart
a useful velocity increment to the payload.
Dumbbell tethers are essentially a coupled pendulum system, and have been well
coveredinthe literatureinthepast,(e.g.withHugen’sclocks,see[Bennettetal.,
2002]). More recently, [Xu et al., 2005] finds that the harmonically excited parametric
pendulum can be simplified to the Mathieu equation, opening up this area
to a new branch of analysis.
Cartmell investigates parametric excitation of the MMET system, finding that
the motor torque may be minimised while guaranteeing a monotonic spinup of the
system. This is achieved by a harmonic modulation of a portion of the payload
masses, and by manipulation of the radius of gyration of the additional lumped
27
masses. Alternatively, given a constant motor torque, the response of the tether
may be actively chosen by selecting a bifurcation path.
MisraandModi, in[MisraandModi,1992],provideanexcellentLagrangianformulationfor
the threedimensional motion of nbody tethered systems using a multiple
pendulum modelasthebasisforthe equations. The equations ofmotion are validfor
large motion and variable lengths, however, the tethers are massless and straight.
This links the nbody pendulum equations and tethered bodies, and provides an
excellent template for rigid body motion with many linked tethers.
Moravecpublished the ideaof themomentumexchangedevice in[Moravec,1977],
defining the Lunavator4 as a momentum exchange device that touches the surface
of the moon and picks up, or drops, a payload.
Thedeploymentcharacteristicsof atethersystemare investigated in[Modi and
Misra, 1979], with a thorough and indepth discussion on the derivation of Lagrange’s
equations of motion for a rotating tether with vibrations. The equations of
motionfor athreemass system(two end masses and a massive5 tether) are identified
as coupled and nonlinear. They study the tether close to the atmosphere, and
conclude that the outofplane motion decays provided the tether is outside Earth’s
atmosphere.
Modi andMisra continue their academicpartnership, and report on controlling the
tether in [Modi et al., 1982], control of tethers used in Space Shuttlebased tethers
in[Modi etal.,1992]and tethered elevators in[Modi etal.,1993].
Perhaps with the Shuttle retiring in 2010, the Shuttle will be proposed less frequentlyasanend
massforthetether,as in[KyroudisandConway,1988]and[Pascal
et al., 1999], instead replaced by the upper stage of an Ares rocket. It is interesting
to note that Pascal mentions ‘crawlers’ – moving payloads along the tether as
4The Lunavator works best on planetary bodies with no atmosphere, where there is no atmosphere
to drag against the tether, and on Mars, where the combination of favourable rotation and
relatively weak gravity suits the system. However, the name Martivator is clumsy, and has not
caught on.
5Massive in the sense of a body with a nonzero mass, not massive as in very large.
28
an alternative to full deployment. A variation of this moving mass system used to
control tethers may also be used to control spacewebs, as will be discussed later in
this thesis.
2.2.4 Motorized momentum exchange tethers
Torotatethetetherto launch velocity, it isnecessary toprovideameansof rotation.
On Earth, this is provided by the Earth, however on orbit, one must bring the
counterrotation device along with the launch device – in a similar manner to a
helicopter’s tail boom.
[Eiden and Cartmell, 2003] briefly summarises the possible role of a European
roadmap for nonconductive tethers, nominally based on momentum exchange, and
also for conductive tethers in which electrodynamics, alongside momentum exchange,
provide propulsion. In the case of the former class small and large payload
deorbitareseenasneartermgoals,withfree flying tetheredplatformsand artificial
gravity systems in the midterm, followed eventually by spinning tethers providing
interplanetary propulsion. Gravity gradient stabilisation is an important underpinning
phenomenon when considering spacecraft stability, and this is particularly the
casefor long momentumexchangetethers. In[Cartmellet al.,2003],dumbbellmodels
are considered for momentum exchange tethers, and offshoots and developments
of this work have shown conclusively that hanging, librating, and spinning tether
motions are intimately connected to this fundamental phenomenon, as reported in
[Ziegler and Cartmell, 2001]. Ziegler shows in [Ziegler and Cartmell, 2001] that
the rotating MMET outperforms a librating tether by two orders of magnitude in
boosting the apogee of a payload.
An indepth treatment of the rigid body dynamics of tethers in space is given
by [Ziegler, 2003]. In this work the dumbbell tether is modelled at various levels of
accuracy, and approximate analytical solutions are obtainedby means ofthe method
29
of multiple scales6 forperiodic solutions. Comprehensivedynamicalsystems analyses
are summarisedfordifferent configurations and models, andglobal stability criteria
for a rigid body dumbbell tether, in both passive and motorised forms, are defined
and investigated.
The method of multiple scales is described further in [Cartmell et al., 2003] in
a review article with several examples. With respect to tethers, the method of
derivation is covered in great depth. An analytical solution to the MMET system
is derived using the multiple scales method in [Cartmell et al., 2006], where the
equations describe a forced and parametrically excited Duffing oscillator.
Sorensen has authored two interesting papers ([Sorensen, 2003] and [Sorensen,
2001]) detailing a method of using Momentum eXchange/Electrodynamic Reboost
(MXER) tethers powered by ED tethers to provide interplanetary travel through
hyperbolic escape from Earth’s gravity. The second details a boost station that
sequentially raisesthepayload’s orbitthrough repeated rendezvous and momentum
exchanges. Of additional note are the appendices detailing analysis techniques and
mathematical derivations that can be used in tether facility design.
With any rotating system,thepossibility of rotational instability arises. [Valverde
and vanderHeijden,2008]find that amagneticbuckling of awelded rod isfound to
be described by a surprisingly degenerate bifurcation. This has serious implications
for using ED tethers as well as static (conducting and nonconducting) deployed
tethers for rotational launch systems.
[van der Heijden, 1995] covers the effects of bearing clearance driving a nonlinearly
coupled driven oscillator in rotordynamics, with close parallels to the MMET
system. Partial modelocking then leadstoquasiperiodicmotionof therotors,and
perturbations of the system lead to more complicated motions.
6Themultiplescalesmethod waspopularised in[Nayfeh andMook,1979].
30
2.2.5 Staging with momentum exchange
As the fundamental material limits of tethers limit the useful V
that may be
imparted, novel techniques such as staging tethers are explored.
A short and digestible article by Forward and Hoyt in Scientific American (see
[Forward and Hoyt, 1999]) first showed the ingenious concept behind the use of
multiple, staged tethers to increase velocity without the necessity of extreme design,
for EarthMoon payload transfer. [Forward, 1991] continues this work by pairing
theLunavator with aLowEarthOrbit(LEO)tether.
Theconceptof staging isexpandedfurther in[McInnes andCartmell,2006] and
[Cartmell et al., 2004].
Following upfromthegroundtestsof amomentum exchangetetherperformed on
ice, [Cartmell et al., 2003] takes this a stage further, and rigorously examines the
release and despin of payloads.
Largescaleelectrodynamicandmomentumtransferbetweenplanetswasstudiedin
theNASAMXERprogram,by Forward,Hoyt,Sorensen and others. Initial concept
design work in [Sorensen, 2001] outlined a massive tether system to send 20  80
tonnes to Lunar or Mars orbit, with a tether length of up to 100km. Sorensen
finds that hyperbolic injection is possible in [Sorensen, 2003], and this opens up
the MXER system to send payloads to Mars. It must be stated that the necessary
velocity to launch a payload to Mars is not attainable with current materials7, to
do so would require a significant increase in the specific strength of the tether.
2.2.6 Deploying and recovery of tethers
Further work on using tetherbased transfers was reportedby[Lorenzini et al.,2000]
in their landmark paper in which staged tethers in resonant orbits are proposed
7TheV
forLEOtoLMO(LowMarsOrbit)isapproximately6.1km/s
andthe maximum characteristic
velocity of a single strand tether is approximately 2.7km/s. See Appendix D, Table D.1
for current materialproperties.
31
as being more mass efficient than single tether systems8 . Lorenzini et al. briefly
refertotether orbit raising results citedby[Carroll,1986] for radial separation as a
function of tether length, and conclude that spinning staged tethers could provide
an ideal transferrateof fivetransfersperyear. Thetransferrateof astaged system
is determined by the periodic realignment of the apsidal lines of the two stages,
whereas in the case of a single tether it is dependent on the time required for re
boosting the stage.
Continuing with the theme of propulsion of a small payload tethered to a large
mass intheformof aspacestationor largeShuttle,[Pascal et al.,1999]investigated
the laws of deployment and retrieval by means of a three dimensional rigid body
model of adumbbelltether inboth circularand elliptical orbits,expanded in[Pascal
et al., 2001]. Several laws are proposed and analytical solutions for small planar
and nonplanar motions of the tether are given, showing that equilibrium tension
can be stated as a function of instantaneous tether length and corresponding axial
acceleration,forwhich controllawscanbestipulated. It isshownthatdeployment is
generallystablewhereasretrieval isnot. Various lawsareexaminedfordeployments
and retrievals, and also for crawler configurations in which the end payload moves
out along a predeployed tether and how this can mitigate the inherent instability
of retrieval.
The next conceptual step to take when considering deployment is to include some
form of flexibility within the tether, and an interesting study of this was published
by[Danilinetal.,1999], inwhichtheelastictethermodel of[NoandCochranJr.,
1995]isusedbut withdifferent variablesandderivation. Theobjectiveof thispaper
was to consider deployment of a completely flexible tether from a rigid rotating
space vehicle under the influence of a central gravitational field. The tether is
modelled using a Newtonian derivation as a series of discrete masses interconnected
by masslesselementsand with internal viscousdamping. Theauthorsmakethevery
8Theypropose using afacility topayload mass ratio of1 :3 using current materials.
32
important point that tether element forces cannot be compressive, so conditions
within the numerical solution algorithm have to be set up to accommodate the
consequential folding effects. Two numerical examples are summarised; one for a
swinging terrestrial cable with an end mass, which starts from a horizontal initial
condition, mainly as a verification of the model in those conditions; and the other
for plane motion of a space vehicle deploying a relatively short 3km
tether, with
elemental spacing of 100m. The deployment is linear and conditions are set up to
apply smooth braking of the tether to a halt at the end of the deployment.
2.2.7 Capture and rendezvous of tethers
An implicit assumption with tether staging is that the payload will, at some point,
have to rendezvous with a second system. [Lorenzini, 2004] provides an indepth
treatment of a spinning tether loop with an extended time opportunity for errortolerantpayload
capture withinhighV
propulsiontoGTOandEarth escape orbit.
Theconfiguration issuch thattheendsof the loop arefurthest awayfromthecentre
of mass, where the loop is at its narrowest. The concept is simple in principle. The
satellite contains a harpoon which shoots through the loop and hooks onto the loop
to capture the satellite. This makes it tolerant of large longitudinal position errors
and reasonable lateral errors as well as some outofplane error. The capture is soft,
and so caters for some velocity mismatch.
[Williams et al., 2003] outlines a ‘momentumenhanced gravityassist’ technique
to capture a payload on a previously hyperbolic orbit into planetary orbit. This is
feasible in theory, however requires accurate control and a large tether deployment
rate, which may limit the practical usefulness of the technique.
A subsequentpaper,[Williams et al.,2005], systematicallydescribestheproblems
in payload capture (spatial and temporal matching, postcapture dynamics and
failure concerns), and designs a suitable rendezvous and control system involving a
33
crawler to successfully capture a payload.
2.2.8 WSB trajectories
Perhapsoneof themostinteresting concepts intrajectory analysis,atleast asfaras
thetether isconcerned, istheWeakStabilityBoundary(WSB) trajectoryproposed
by [Belbruno, 1987]. This technique allows the majority of the V
for an Earth
Moon mission to be concentrated in the first burn, with a relatively minor capture
burn at the end of trajectory.
If this is compared to the V
requirement of the Hohmann twoburn technique9 ,
it becomes clear that a tetherbased launcher system could provide a more efficient
route to the Moon. A thoroughly good overview of the WSB method is given in
[Ross,2004,p112],where itisstated thatthefuel required isloweredby about20%.
As is common with large research groups, there are many valuable papers on WSB
transfers and the Interplanetary Transport Network from the CalTech researchers
that are listed in Ross’ thesis.
Apractical exampleof theWSB method isexplained clearly in[Koonet al.,1999],
giving an example of a trajectory similar to the Hiten Lunar probe. The WSB
trajectory10 has since been used to successfully launch NASA’s Genesis mission to
study the solar wind(see[Lo andRoss,2001]).
TheWSBtrajectoriesareanumbrellaforamultitudeofpossibleroutes, incontrast
to the simplicity of the Hohmann or bielliptic trajectories. A useful geometric
method of finding such a WSB trajectory is outlined in [G´omez et al., 1993], taking
thedestinationastheEarthSunL1point. Oncethispoint isreached,thepathways
are open to a great many places in the Solar system.
9Hohmann’sseminalbookon ‘TheAttainability ofHeavenlyBodies’maybefoundin itsoriginal
form[Hohmann,1925](inGerman) ortranslated toEnglish in[Hohmann,1960]
The WSB trajectory is hard to visualise and perhaps 2D plots of the spacecraft’s path in
inertial and rotating space are difficult to extrapolate mentally to 3D space. Thankfully Lo
has made videos of the spacecraft path on his website: http://www.gg.caltech.edu/~mwl/
communications/communications2.htm
34
A method of using a MMET with a WSB trajectory is outlined in [McKenzie and
Cartmell, 2004]. A small satellite MMET system, under 1000kg, is able to launch
a 10kg
microsatellite from MEO to Lunar capture orbit in around 100 days. This
system may be recovered and relaunched every month until the useful life of the
satellite is reached. This relies on a very small safety factor of 1.3 in the main
tethers, and does not account for the additional stator arm in the MMET, however
this shows that the current state of technology is almost ready to sustain an orbital
tether launch system.
2.3 Tether derived structures
In the same way that orbital tethers have evolved from the space elevator, space
structures may be constructed from tethers.
2.3.1 Deployed space structures
In1968,PeterGlaser suggested that a solar collector11 couldbeplacedingeostationary
orbit [Glaser, 1968]. Operating in high Earth orbit, this would use microwave
power transmission to beam solar power to a very large antenna on Earth where it
can be used in place of conventional power sources. The advantages to placing the
solar collectors in space are the unobstructed view of the Sun, the fact that it is
then unaffected by the day/night cycle, weather, or seasons. For efficient operation,
the satellite antenna must be between 1 and 1.5km
in diameter and the ground
rectenna around 14km
by 10km, generating between 5 and 10GW
of power. One
option of constructing this huge collector is manufacturing the cells on Earth, and
constructing them in space with the aid of robots.
An alternative approachtothe onorbit assembly ofsuch a massive structure would
betheuseofwebs,as in[Kayaet al.,2004a]and[Kayaet al.,2004b].Alargegeneric
11That is, a satellite for collecting solar power.
35
netlikestructurecould firstlybedeployed inorbit. Oncethenet isstabilised,robots
could be used to crawl over the net towards specified locations and move solar cells
into the desired positions. This enables the construction of such a satellite using
robot assistance to be faster and cheaper12 .
Other largespacestructureshavebeenproposed that wouldmakeuseof ageneric
reconfigurable platform. In the case of the orbiting stellar interferometer, [DeCou,
1989] showed that planar deformation of a spinning system comprising three collimating
telescopes at the corners of an equilateral triangle made up from three interconnecting
tethers would be inevitable due to the inertia of the tethers. Clearly
inertialess tethers will not deform centripetally and will, instead, merely stretch
intostraight linesduetothetensioncreated along their lengthby thecornermasses
as the whole system rotates. Perhaps the addition of a control system would make
this particular structure a more realisable option.
This directly led to the successful deployment of the Inflatable Antenna Experiment(
IAE),andled atrancheof successfullydeployed structuresinspace, including
the ESA Cluster mission as detailed in [Andi´on and Pascual, 2001].
Studies have previously shown that robots may be deployed in this way to reconfigure
the net structure [Kaya et al., 2004a], [Nakano et al., 2005]. A international
collaboration from the University of Glasgow, KTH Royal Institute of Technology
in Sweden and the ESA plans to test the concept on a sounding rocket in 2010, as
reported in a [University of Glasgow Press Release, 2009].
2.3.2 Spacewebs
If the generic satellite is to be used as an element of infrastructure, it is desirable
to start studying the dynamics and deployment of such a structure. A type of
12To a degree, the question of whether a method of power generation is ‘better’ is largely a
political decision.
36
satellite, called ‘Furoshiki’13 is outlined in [Nakasuka et al., 2001], with an overview
of the concept, a simulation of a membrane satellite and an evaluation of the control
strategy needed to stabilise such a satellite.
Further work on the deployment of a square membrane is examined in [Tibert
and G¨ardsback, 2006], where a membrane is attached to four daughter satellites
and deployed while the satellite is both free and rotating. It is shown that the free
deployment without damping is not possible, however a torque can be applied to
the central hub to stabilise the deployment.
Once the membrane is deployed, robots such as reported in [Kaya et al., 2004a]
may be tasked to assemble or reconfigure the structure.
A multitethered structure containing n
masses is investigated in [PizarroChong
andMisra,2008],for an openhub, closedhub(wheel) anddoublepyramidformations.
A system of Lagrange’s equations of motion is constructed and the stability
of the three configurations were tested in a variety of configurations – parallel and
perpendicular to the orbital plane and with and without an initial spinrate. The
dualpyramid formation was found to be the stable in most tests, with all configurations
benefiting from a rotating platform.
A similar analysis in [Wong and Misra, 2008], involving a multitethered structure
near the SunEarth L2 point, examines the planar dynamics of the satellite. Lagrange’s
equations are used to model the system, with a linear feedback controller to
control thetether lengthsand angulardisplacementsof thesystem. Theend masses
were successfully controlled, wielding a spiral pattern which allows the satellite to
be used as an optical interferometer.
13After the traditional Japanese wrapping cloth.
37
Chapter 3
Tether modelling
The steps to construct a mathematical model of a tether orbiting a planet are
outlined in this chapter. Starting with a simple model of a point mass orbiting a
massivebody,each stagebuildsonthe lastandprogressivelyadetailed(and more
importantly, validated) model will be outlined.
3.1
Constructingtheequationsof motionof asystem
Amethod todescribeasystem intermsofOrdinary DifferentialEquations(ODE)
is outlined below. For each mass point of the system:
.
the position of each mass point is described in local coordinate space
.
theCentreofMass(CoM) of thesystem isfound
.
the position is translated and rotated1 into inertial coordinate space, using
matrix and vector operations
1The centre of mass effectively links the two coordinate systems by providing a common point
of reference between the two.
38
then, for the system as a unit:
.
the Lagrangian is found by expressing the total energy of the system
.
generalised coordinates are chosen to suit the system.
.
velocities of each mass point are derived
.
Lagrange’s equations are constructed
.
Lagrange’s equations may be numerically integrated using suitable initial conditions
The remainder of this chapter will fully outline the steps needed to construct and
validate ageneral system of equations,focussing specifically on constructing aplanar
dumbbell tether on a circular equatorial Earthcentred orbit.
3.1.1 Validating code modules
It isessentialaspart ofthescientificmethodtoset up anexperiment,or inthiscasea
numerically solved mathematical model, that is both repeatable and verifiable. To
this end, a method of setting up equations describing each component’s position
using repeatable and testable modules is introduced here.
It is of fundamental importance, when building a mathematical model, to ensure
that the model is verified and validated insofar as practically possible. As such,
a method of constructing a system of validated modules is introduced to verify
that, for example, a module to rotate a vector about an axis does this accurately
and repeatably. Therefore, for each component of the system, a testable module is
produced and the system as a whole is much easier to validate as the sum of its
parts.
In Mathematica.
, modules of code may be configured and tested such that a
39
oneline testable routine may be evaluated. Assigning:
..
10 0
..
..
..
XRotation[]= 0 cos() sin() (3.1)
..
..
0 sin() cos()
it isthenpossibletoseparateand logicallyevaluatetheright and lefthand sidesof
Equation 3.1 to see if they are identically equal. If the right hand side is identically
equal to the left hand side, the code is validated, if not, the code fails validation.
In this way, it is possible to establish a repeatable and valid technique for constructing
the equations of motion of the system.
3.2 Centre of mass modelling
3.2.1 Mass points & body selection
It is essential when using the Lagrangian method that the movement of the bodies
concerned are described accurately. The modelling tasks can be broadly separated
into expressions concerning:
.
bodies containing mass which contribute to the kinetic and potential energies
.
external or internal conservative forces that may contribute to the energy sum
(e.g. springenergies)
.
external or internal forces that may contribute to the nonconservative right
handside terms in Lagrange’s equations
As such, the principles of system modelling will be demonstrated using a simple
exampleofthecoplanardumbbelltether incircularEarth orbit,shown inFigure3.1.
40
Figure 3.1: Symmetrical dumbbell tether
Central to the system, the facility mass houses the entire structure of the system;
the power generation, communications equipment, etc.. This is expected to form
the majority of the mass of the system.
Thetwopayloads arethe raison d’ˆetre of the momentum exchange tether system.
They are assumed to have mass and are rigidly connected to the facility by the
tethers and are usually a significantly smaller fraction of the system mass.
The two tethers are assumed to be massless. This is usually a significant assumption
to make because the fractional mass of a 100km
tether is nontrivial. In this
case, itisareasonableassumptiontodemonstratetheprinciplesof tethermodelling
without making the mathematics overly arduous.
Thetethersholdthepayloadsrigidlyand symmetrically inplaneaboutthefacility
mass, ensuring that the tethers are colinear about one common axis through the
centre.
Theotherbody of note is,of course,theEarth which acts indirectlyonthissystem
providing the gravitational potential energy.
3.2.2 Centre of mass calculations
The centre of mass can have a considerable effect on the motion of this system.
Consider an asymmetrical payload release, for example, when only one payload is
released instead of both, as normally occurs. The centre of mass would undergo an
instantaneous stepchange and, consequently, the orbit of the system would change.
Assuming the dumbbell is in a circular orbit aligned along the vertical (gravity
gradient) direction, if the uppermost payload is released the excess kinetic energy
41
of the orbit will cause the payload to describe an ellipse with the separation point
as the perigee. The remainder of the system2 will tend to deorbit as the potential
energy of the system is lower than needed to sustain a circular orbit. This separation
point becomes the apogee of the new system’s orbit, as shown in Figure 3.2.
Earth
Figure3.2: Symmetricaldumbbelltether on a circular orbit(solid). Afterpayload
release,thefacility systempath at apogee(dash) andpayload atperigee(dots)
Inasymmetrical system,thecentreof mass is located coincidentonthefacility. If
thesystem isasymmetrical,thenthecentreof masswillneed tobecalculated. This
is highly important to the accuracy of the model and in most cases necessary to
integrate the equations of motion of the system. Most of the problems encountered
whennumerically integratingtheequationsof motion inMathematicaweresolvedby
inserting a mathematical expression that accuratelydescribes the(variable)position
of the centre of mass in the positional equations which then feed into Lagrange’s
equations.
The Centre of Mass (CoM) calculations give the location of the CoM from an
arbitrarily chosen point. In this case it is convenient to chose the facility mass as
theoriginpoint(and thecentreof a local coordinatesystem) and refertoallother
2That is, the lower section of the system comprising the facility, two tethers and one payload.
42
points in the tether system relative to the facility location.
This may be alternatively expressed as constructing a local coordinate system
{Xlocal,Ylocal,Zlocal}centred on the facility mass, in this case, aligning the Xaxis
with the radius vector.
The position of the centre of mass about the facility for n
masses each with positions
{Xi,Yi,Zi}, in terms of the local {X,
Y,
Z}coordinate system centred on the
facility mass is:
.
.
n
n
n
.
.
.
.
.
MiXi
MiYi
MiZi
.
.
.
.
i=1 i=1 i=1 Pfacility.CoM = ,
,
(3.2) n
n
n
.
.
.
.
.
.
Mi
Mi
Mi
.
.
.
i=1 i=1 i=1
Once the CoM relative to the facility mass has been found, a vector describing the
position of any mass point in the inertial frame of reference3 may be constructed.
This is used in later equations by taking the position from the centre of mass to
the facility mass:
PCoM.facility = Pfacility.CoM (3.3)
3.2.3 Positions of mass points
A technique to obtain the location of each mass point relative to the facility is
outlined below, using a 2D dumbbell system as an illustrative example.
To assemble the equations of motion, a series of rotations and translations are
performed on the initial vectors to transform the local space coordinates into inertial
space coordinates. The initial vectors are chosen in synchronisation with the
generalised coordinates and the rotation angles.
Figure3.3showstheorbital motion,constrained tothe inertialXYplaneand the
local XY plane, of an asymmetric dumbbell tether. Note that the local Xaxis is
3In the Earth centred system.
43
Earth
LYlocalXinertialXlocalCoMRadius=P0!CoM
YinertialjPCoM!facilityEarth
LYlocalXinertialXlocalCoMRadius=P0!CoM
YinertialjPCoM!facility
Figure 3.3: Asymmetrical dumbbell in inertial axes
parallel to the radius vector.
The inertialpositionvectorPinertial of each masspoint, j, isassembledfromthevector
sum of the component vectors, taking into account the translation and rotation
of each component.
Pjinertial = R,Z · P0.CoM Pfacility.CoM + R j
,Z · Lj
(3.4)
Where P0.CoM = {R,
0,
0}, and represents the vector of the centre of mass from
Earth4 . The rotation matrices and matrix notation R j
,X are explained in Section
3.5.
The five stages in constructing the position of the mass point Pj
are shown in
Figures 3.4, 3.5, 3.6, 3.7 and 3.8, and the position vectors for an arbitrary mass
point are given in Equations 3.5, 3.6, 3.7, 3.8, 3.9 and 3.10.
Figure 3.4 shows the mass point aligned along the Xaxis of the local coordinate
system, in this case using a payload mass with length L
as an illustrative example.
4i.e. the Radius vector.
44
Ylocal
LPjXlocalfacility
Figure 3.4: Local coordinate construction
Ylocal
L
PjXlocalfacility
Figure 3.5: Local coordinate construction after rotation through angle .
about the
facility mass
45
YlocalXlocalPj
LPCoM.facilityfacilityCoMYlocalXlocalPj
LPCoM.facilityfacilityCoM
Figure 3.6: Local coordinate construction after translation from the facility to the
CoM
Earth
YlocalXlocalPj
LPCoM!facilityfacilityCoMYinertialXinertialR
Figure 3.7: Local coordinate construction after translation from the CoM to the
inertial coordinate system
46
Earth
YinertialYlocalXlocalPj
LPCoM!facilityfacilityRCoM
XinertialEarth
YinertialYlocalXlocalPj
LPCoM!facilityfacilityRCoM
Xinertial
Figure3.8: Local coordinate construction after rotation through .
about the inertial
coordinate system
Rotating this through an angle of .
about the facility as shown in Figure 3.5. This
is translated to align the origin with the centre of mass as shown in Figure 3.6.
A further translation is performed through the inertial Xaxis about a distance R,
giving the origin as the inertial coordinate system, as shown in Figure 3.7. Finally,
the inertial coordinatesystem isrotated abouttheoriginaround anangle .
as shown
in Figure 3.8.
47
The vectors of each are shown below:
..
L
..
..
..
P1 = 0 (3.5)
..
..
0
..
L
cos( )
..
..
..
P2 = R ,Z · P1 = L
sin( ) (3.6)
..
..
0
..
L
cos( ) CoMX
..
..
..
P3 = P2 PCoM.facility = L
sin( ) CoMY (3.7)
..
..
CoMZ
..
L
cos( ) CoMX +R
..
..
..
P4 = P3 + R
= L
sin( ) CoMY (3.8)
..
..
CoMZ
P5 = R,Z · P4 =
..
cos.
(L
cos( )CoMX + R)sin.
(L
sin( )CoMY)
..
..
..
sin.
(L
cos( )CoMX + R)+cos .
(L
sin( )CoMY) (3.9)
..
..
CoMZ
P5 = R,Z · P0.CoM Pfacility.CoM +(R ,Z · P1) (3.10)
The process is repeated for each mass point under consideration.
3.3 Constructing Lagrange’s equations
Lagrangian mechanics is primarily concerned with the trajectory of an object, derivedbyfinding
thepath which minimisestheaction,aquantity which isthe integral
48
of the Lagrangian over time. As this is an energy based approach, the process of
evaluating the kinetic and potential energy of the system will be outlined and the
choice of generalised coordinates will be explained.
3.3.1 Energy modelling
For every mass point in the system with position in the inertial frame, Pj, the
following steps are undertaken in order to find the Lagrangian energy expression L:
1. the respective velocities, Vj
are found:
.
Vj
= Pj
(3.11)
@t
2. thekinetic energies(linear Tlin = 21 mv2 and rotational Trot = 21I!2) are obtained
and summed:
n
1
Tlin = mj
Vj
· Vj
(3.12)
2
j=1
n
1
Trot = mj
Pilocal · Pjlocal !j
· !j
(3.13)
2
j=1
where the vector of angular velocity, !j
is the first derivative of the angular
position of each point, {j,j, j
}
3. the potential energies, U, are obtained and summed:
nn
µmj
µmj
U
=
.
=
.
(3.14)
Pj
Pj
· Pj
j=1 j=1
where µ
is the standard gravitational parameter. For the Earth, the value is
µ
=3.986* 1014 m3/s2 .
4. and the Lagrangian is found:
L
= Trot + Tlin U
(3.15)
49
For each mass point, the Lagrangian energy expression is constructed by considering
the total kinetic and potential energies of the system.
d
@T
@T
@U
 += Qj
(3.16)
dt
@q.j
@qj
@qj
Lagrange’s equations are generated for all the generalised coordinates as specified
in Equation 3.16.
The forces are then used to calculate the right hand side of Lagrange’s equations,
Qj
, through consideration of the virtual work.
Primarily, the main nonconservative force acting on the tether system will be the
motortorque. Thisconcept iscoveredfurther inZiegler’sPh.D.thesis[Ziegler,2003,
p4345].
3.3.2 Choosing generalised coordinate systems
Thegeneralised coordinatesused inderivingLagrange’sequations inSection3.3are
chosen to be independent coordinates. These coordinates need not be orthogonal
nor Cartesian, in fact it is limiting to assume they are either, however they must
describe an independent degree of freedom of the system.
In the case of the orbital variables R
and , these are chosen in preference to
the inertial variables Xinertial and Yinertial to enableperiodic analysis overthe orbital
cycle.
Inthe localcoordinatesystem,thechoiceof ageneralised coordinate is lessstraightforward.
In a two generalised coordinate system with an inplane and outofplane
motion, there are many angles that could be chosen as the generalised coordinate,
andthisposes aproblem over choice of variables. The equations of motion(and
therefore the corresponding set of Lagrange’s equations) that describe one system
of rotations would be different from the others, but the motion would have to be
50
identical when comparing the many possible rotations systems. This implies that
for the motion to be identical5 there could be many possible systems of equations
that could, when integrated, provide the motion of the system, and that there may
be one or more systems of equations that provide a more elegant solution than the
rest. Two systems of rotations are investigated in Section 3.5.
The goal in selecting a “good” choice of generalised coordinates is that they are
easily interpreted in graphical form, and the equations generated are compact and
meaningful to interpret.
3.4 Nonconservative forces
Any nonconservative forces that appear in the system need to be included in the
equations of motion in terms of virtual work. These are contained within Ziegler’s
Ph.D. Thesis [Ziegler, 2003, p4345], and will be contextualised here.
The nonconservative force is represented by Qj
and can be evaluated by considering
the virtual work done by the motor torque summed over each mass:
@Pmass
Qj
= F
· (3.17)
@qj
mass
with j
representing each of the n
chosen generalised coordinates and Pmass as the
position vector of the mass under consideration. F
is the nonconservative force.
The3×3rotationmatrixforageneralised rotation inthetetherbody axes6 (as will
be explained in Section 3.5) is multiplied by the length vector to give the position
5If the systems were nonidentical then, generally speaking, there would either be flaw when
choosing the generalised coordinates or a mistake in the equations. The process of simulating the
motion of a system should be entirely independent of the means used to arrive at the end result!
6Rotations are only provided here for the Yand Zaxis.
51
of each mass in the tether body axes :
..
Lmass
.
.
.
.
Pmass = R,Y · R ,Z · .
0 .
(3.18)
.
.
.
.
0
The resultant 3 × 1 vector is differentiated with respect to each of the n
chosen
generalisedcoordinates togive an n×3matrix. Theforcevectorexpressed intether
body axes isthenfoundby multiplying thegeneralised rotationmatrixby theforce
vector:
..
0
..
..
..
Force = R
,X · R,Y · R ,Z · 0 (3.19)
..
..
/Lmass
where t
is the nonconservative force converted into the rotating frame.
Thedotproduct of the3×1 force vector and the n
×3 matrixthengives an n
×1
vector containing the Qj
terms for the n
chosen generalised coordinates to be used
in the right hand side of the equations of motion.
3.5 Rotation systems
A comparison will now be made between two rotation systems that specify the
position of an arbitrary mass point in the local axes system.
As discussed earlier7 there are multiple ways to specify the position of the mass.
No one system may be termed the definite solution,however, one may bepreferable
to choose if it generates a more easily interpreted system of equations.
7See section 3.3.2.
52
3.5.1 Rotation summary
Three rotation matrices are set up to describe the direction of a single 3D vector
about its origin using a series of rotations about the local X, Yand Zaxes.
As the three rotations can rotate an arbitrary vector around the origin, the initial
vector must be defined. For the rotation systems that follow, a vector aligned with
the Xaxis in the positive direction starting at the origin has been chosen. The
series of rotations are then performed sequentially to give 3 new coordinate axes
(two intermediate andone final) for 3 rotations.
The Rotation R,B is defined here as a rotation matrix acting about the Baxis
and rotating about an angle .
The new vector is found by rotating the starting vector coordinates by means of
three rotation matrices about the local origin. The rotation system will be one of
the six permutations from the list of three possible rotations: R,X, R,Y and R ,Z
given in Equations 3.20, 3.21 and 3.22.
Therotationsystemthatrotatesthepoint first around theXaxis,thentheYaxis,
then the Zaxis is termed the XYZ rotation system.
..
10 0
..
..
..
R,X = 0 cos() sin() (3.20)
..
..
0 sin() cos()
..
cos() 0 sin()
..
..
..
R,Y =
.
0 10
.
(3.21)
..
sin() 0 cos()
..
cos( ) sin( )0
..
..
..
R ,Z = sin( ) cos( ) 0 (3.22)
..
..
0 01
53
The order of matrix multiplication is important; for the XYZ rotation sequence,
the vector multiplication is evaluated sequentially from the right hand side:
R ,Z · R,Y · R,X · Porig.
= R ,Z · R,Y · R,X · Porig.
Notewell,thisdiffersfromtheZYX rotationsystem intheorderof rotations. The
ZYX rotation system first rotates the vector by the Zaxis, then the Yaxis, then
the Xaxis.
R,X · R,Y · R ,Z · Porig.
= R,X · R,Y · R ,Z · Porig.
When looking at the matrix rotation equation, this may seem counterintuitive to
label the rotational sequence from right to left8 . However, it reflects the sequence
of rotations, even if it is backwards from a point of view of reading the matrix
multiplication subscripts.
Using thethreerotationsgivenhere inthethreeaxes,thereare6possiblecombinations
of defining a rotational sequence, as shown in Table 3.1. These describe the
sequence of rotations when rotating a point around the origin in 3D and the three
possible rotations in 2D that may be derived from the 3D rotation.
The 2D rotations are not unique to a single 3D rotation. For example there are
three XY rotations that may be made in 2D that correspond to XYZ, XZY or
ZXY rotations in 3D space.
Two of these sequences will be investigated using a simple 2D rotation with the
8Indeed,this isaratherarbitrarydecision,but itcouldbearguedthatthisfollowsthearbitrary
direction of text inkeeping with the custom of languages of the Western hemisphere. In choosing
the nomenclature, it makes sense to make XYZ define the rotational sequence rather than the
order in which the matrices are multiplied purely because most nonspecialists would grasp the
concept of order of rotations before learning threedimensional matrix algebra!
54
Order of rotation in . . .
3D 2D
R,
R,
R ,
1 2 3
X Y Z XY XZ YZ
X Z Y XZ XY ZY
Y X Z YX YZ XZ
Y Z X YZ YX ZX
Z X Y ZX ZY XY
Z Y X ZY ZX YX
Table 3.1: Order of rotations in 3D and 2D
sequenceof rotationsabouttheYZandZYaxes, inSection3.5.2andSection3.5.3
respectively.
3.5.2 YZ Rotations
The YZ series of rotations9 rotates a point lying on the Xaxis around first the Y
axis, then the Zaxis. The Xaxis rotation is not included in this system for ease of
visualisation.
Rotating the original vector Porig. = {L,
0,
0}T
, shown in Figure 3.9a about the
Yaxis through an angle of a
gives the following, as shown in Figure 3.9b:
..
L
cos()
.
.
PinterYZ = R,Y · Porig. = .
.
.
0 .
.
.
(3.23)
.
.
L
sin()
Rotating this intermediatestagearound theZaxisthrough anangleof .
gives the
9Or to be specific the XYZ series of rotations in 3D space when the X rotation is zero.
55
following, as shown in Figures 3.9c and 3.9d:
.
.
L
cos()cos( )
.
.
.
.
PfinalYZ = R ,Z · PinterYZ = .
.
L
cos()sin( ) .
.
(3.24)
.
.
L
sin()
A diagram of the complete rotation sequence is are shown in Figure 3.9e.
Theoutofplaneanglegivenby theYZ seriesof rotations istheangle . However,
the inplane angle is not the angle . The inplane angle lies in the plane normal
to the ZYZ plane, and intersects the XZ plane.
The inplaneangle(givenby the XYXYZ plane)isprojectedontotheXYinertial
plane to give the angle .
3.5.3 ZY Rotations
TheZY systemof rotations isgeneratedby rotating thestarting vectoraligned with
the Xaxis around first the Zaxis, then the Yaxis. This can be thought of as the
ZYX series of rotations with a zero Xaxis rotation – therefore giving only a ZY
rotation.
Rotating the original vector Porig. = {L,
0,
0}T
, shown in Figure 3.10a about the
Zaxis through an angle of .
gives the following, as shown in Figure 3.10b:
..
L
cos( )
..
..
..
PinterZY = R ,Z · Porig. = L
sin( ) (3.25)
..
..
0
Rotating this intermediate stage around the Yaxis through an angle of a
gives
56
ZZ Z
Y
y
a
a
(a)Rotations YZ – start
a
Y
a
a
ZYXY
(b)Rotations YZ –
X
ZZ
a
Y
y
a
a
XZXYZXYZY(c)Rotations YZ –
XX
(d)Rotations YZ – a
then .
y
a
Y
y
a
a
YZZYZXZXYZXYZY
Z
a
Y
y
a
a
ZYZYZXYZZYXY(e)Rotations YZ – full rotation shown
Figure 3.9: Rotation sequence for YZ rotation: Rotating about Yaxis then
about Zaxis. is shown as a negative rotation here for ease of visualisation.
57
X
ZZ
Y
h
y
(a)Rotations ZY – start
Y
XZ
XX
h
y
(b)Rotations ZY –
a
Y
h
y
XZYZY
(c)Rotations ZY –
ZZ
XZ
XZ
XX
a
Y
a
ZYXZYXY
(d)Rotations ZY – then
Figure 3.10: Rotation sequence for ZY rotation.. Rotating .
about Zaxis then a
about Yaxis. a
is shown as a negative rotation here for ease of visualisation.
the following, as shown in Figures 3.10c and 3.10d:
..
L
cos( )cos()
.
.
.
.
PfinalZY = R,Y · PinterZY = .
.
L
sin( ) .
.
(3.26)
.
.
L
cos( )sin()
A diagram of the complete rotation sequence is shown in Figure 3.10d.
Conversely to the YZ series of rotations, the outofplane angle given by the ZY
series of rotations is not the angle . However, the inplane angle is now the angle
. The outofplane angle is given by the angle in the XZY XZ plane between the
fully rotated vector and the X0 Y0 inertial plane.
The outofplane angle(givenby the XZY XZ plane) is projected onto the XY
inertial plane to give the angle .
58
The outofplane angle, , may be found by simple trigonometry as follows:
ZPfinalZY
sin.
= (3.27)
L
.
= sin1 (cos( )sin()) (3.28)
3.5.4 Comparing the equations – position equations
In each of the rotation sets, the xcoordinate is identical. The sets differ in the y
and zcoordinate sets as follows:
YZ ycoord: L
cos()sin( ) from Equation 3.24
ZY ycoord: L
sin( ) from Equation 3.26
The ycoordinates differ by a factor of cos()between the YZ and ZY systems.
YZ zcoord: L
sin() from Equation 3.24
ZY zcoord: L
cos( )sin() from Equation 3.26
The zcoordinates differ by a factor of cos( )between the YZ and ZY systems.
The two sets of 3D rotations are different because of the commutative rules of
matrix multiplication, in general when multiplying two rotation matrices A and B
the matrix product AB will be different from BA.
3.5.5 Comparing the equations – kinetic energies
Comparing the differences in the kinetic energies of the spinning systems in a local
coordinate system:
m
YZ Ek
= 2 L2 .2 + cos()2 .
2
m
L2 .2 2
ZY Ek
= 2 2 + cos( ).
Clearly, the two are similar – the difference essentially being the substitution of
one angle for the other between the two sets of equations.
59
3.5.6 Comparing the equations – potential energies
Comparing the differences in the potential energies of the spinning systems in an
inertial coordinatesystem isastagemorecomplicated,asthe local coordinatesystem
mustbetranslatedfromthe local centretothe inertial centreto incorporatethefull
inertial terms.
The potential energy term for each mass point j
µmj
Ep
= v
pj
· pj
v
depends strongly on the magnitude of distance, given by pj
· pj. This evaluates to
the following for YZ and ZY rotations10
YZ: pj
= {R
+ L
cos()cos( ),L
cos()sin( ),
L
sin()}
ZY: pj
= {R
+ L
cos()cos( ),L
sin( ),
L
cos( )sin()}
Therefore
YZ: v pj
· pj
= L2 + R2 +2LR
cos()cos( )
ZY: v pj
· pj
= L2 + R2 +2LR
cos()cos( )
For the potential energy term:
YZ Ep
= v µmj
L2+R2+2LR
cos()cos( )
ZY Ep
= v µmj
L2+R2+2LR
cos()cos( )
That is, the potential energy expression of the two sets of equations YZ and ZY
are identical for an inertially planar system.
When the orbital rotation .
is introduced to the potential energy equations, the
result will be identical to the expression above, as the inertial rotation through R,Z
is a rotation about one axes only, unlike thedual system of rotations outlined above.
Theexpressionforpotential energy willchangesignificantly whenrotated through
10With a simple offset from the origin by R, representing the orbital radius.
60
around the orbit in 3D when the orbital parameter .
is included.
3.5.7 YZ rotation in inertial axes
The full rotation sequence in inertial axes, for a rotation of an initial vector through
a Yaxis rotation, then a Zaxis rotation, translations in CoM and radius vector,
then a final rotation about the Zaxis is as follows:
Pfinal = R,Z · (P0.CoM Pfacility.CoM +(R ,Z · R,Y · P1)) (3.29)
3.5.8 Singularities in the rotational systems
The main disadvantage in using the method of rotational matrices is they have
singularities or ‘gimbal lock’ points where there is a lack of a unique output when a
rotation is performed that aligns it with the rotational axis.
In the case of YZ rotation, an initial rotation of 90. about the Yaxis would align
the original vector with the Zaxis. A second rotation about the Zaxis would be
redundant as the vector would not change value.
In the case of ZY rotation, an initial rotation of 90. about the Zaxis would align
the original vector with the Yaxis. A second rotation about the Yaxis would be
redundant as the vector would not change value.
In general, any rotation through 90., 270. (or any multiples thereof) so that the
original vector aligns itself with the new axes of rotation will produce a singularity.
The equations of motion usually can not be numerically integrated with the initial
condition as a singularity; the lack of uniqueness of the solution violates the method
of using the Lagrangian to solve the equations of motion.
61
3.5.9 Practical uses of the different equations
The two sets of equations above may find different applications, depending on the
desired investigative outcome required.
If a system of equations is needed to solve a solution space near an axis, then
naturally, the system of equations should be chosen to avoid a singularity. As the
choice of system of equations is entirely arbitrary, the problem of singularities has
been avoided here.
Quaternions were considered as an alternative to the threeangle rotation system
above they confer several advantages such as the ability to multiply the rotations
inonepolynomial and they arenotpronetogimbal lock. However,oneof themain
disadvantages of quaternions is that it is very difficult to extract physical meaning
fromthe individual componentswithout asignificantdegreeofpostprocessing. For
this reason, it was decided to remain with the threeangle rotation method.
The YZ system of equations is slightly more compact and easier to interpret than
the ZY equations, although both give identical results when integrated. The difference
in time taken to solve the two sets of equations differs by approximately 5%
using a wide range of initial conditions.
3.6 Numerical integration techniques
The equations of motion, once assembled, are integrated with Mathematica11 using
the NDSolve subroutine. This is a Jackofalltrades solution that effectively deals
withtheproblems of choosingtimesteps, any stiffness that may arise and a multitude
of otherproblemsthatplaguetheresearcher. However, caremustbetakentoassure
that the output of the numerical integrator is representative of real life – validation
and sanitychecking at allstageswillnot ensurepoorresultsareentirely eliminated,
11Specifically, using Mathematica version 5.1
62
but they will reduce the risk of spurious results.
The generated equations are solved with Mathematica’s general purpose solver
NDSolve[. . .] with the following parameters:
.
AccuracyGoal >
Automatic, PrecisionGoal >
Automatic – controls
the number of digits of accuracy, effectively the absolute and relative errors
respectively.
.
With PrecisionGoal>p and AccuracyGoal>a, Mathematica attempts
to make the numerical error in a result of size x
be less than 10a
+ x10p.
.
WorkingPrecision >
20 – limits internal computations to 20digit precision.
From the online Mathematica NDSolve FAQ webpage – [Wolfram, 2007]:
“For initial valueproblems,NDSolve uses anAdamsPredictorCorrector
method for nonstiff differential equations and backward difference formulas(
Gearmethod)forstiffdifferential equations. It switchesbetween
the two methods using heuristics based on the adaptively selected step
size. It starts with the nonstiff method under essentially all conditions
and checksfor the advisability of switching methods every10 or20 steps.
The algorithms and the heuristics for switching between algorithms
are described in the following references: [Hindmarsh, 1983], [Petzold,
1983]”
More indepth information about the NDsolve functions is also available from the
Mathematica Journal: [Keiper, 1992].
63
Chapter 4
Dynamics of tethers with
inclination
The simple dumbbell tether model can be expanded to include an inclination
component to the orbit.
Theprocessofaddingthe inclinationtermwillbeexemplified with atetherpayload
deliverytotheMoonunderatimevaryinginclinationgivenbytheSaroscycle,which
describes the inclination variation of the Moon’s orbit around the Earth.
The Saros cycle describes the Moon’s orbital inclination cycle, lasting 18 years 11
daysand8hours,cutting theSun’seclipticplane(usedasthe inertial coordinate
plane). This cycle determines, in the short term at least, the Lunar and Earth’s
eclipse frequencies and locations. This also allows the inclination of the Moon to
be accurately predicted from ephemeral data. Figure 4.1 shows the inclination as a
function oftime over oneSaros cycle withdata takenfrom theJPLHorizons website
[JPL, 2008].
As shown in Figure 4.1, the Saros cycle varies from a maximum inclination with
respect to thegeocentricplane of approximately 28.8. to a minimum of18.0., which
64
Moon DEC_(ICRF/J2000.0)
30
20
10
0
10
20
30
Moon Declination (deg)
01/01/2024
01/01/202301/01/202201/01/202101/01/202001/01/201901/01/201801/01/201701/01/201601/01/201501/01/201401/01/201301/01/201201/01/201101/01/201001/01/200901/01/200801/01/200701/01/200601/01/200501/01/2004
Date
Figure 4.1: Ephemeris of Lunar orbital inclination [JPL, 2008]
poses a problem to the trajectory and mission specialists of a longterm reusable
launch system, like the MMET. The position of the Moon changes, which makes
it more difficult to launch a Lunar mission, and the cost of the mission will change
due to a potentially expensive1 plane change manoeuvre.
Comparing the Lunar inclination to Earth’s tilt of 23.5., the Moon’s inclination
with respect to Earth’s equatorial plane varies by approximately ±5. . The V
required to perform an inclination change will be small, but the dynamics of the
system will depend on the inclination with respect to the inertial plane. In 2008,
this is approximately 28. .
4.1 Addition of inclination to tether model
The expansion of the tether model to include an inclination term is relatively
straightforward, at least in terms of the model outlined in Section 3.5.7. The inclination
term is added with a rotation matrix2 R,Y before the final inertial rotation
1In terms of fuel and therefore cost.
2Rotating around the Yaxis through the inclination angle, .
65
sequence R,Z is performed.
Thefullrotationsequence in inertialaxes,forarotationof an initial vectorthrough
a local Yaxis rotation, then a local Zaxis rotation, translations with the CoM and
radiusvectors,thentworotationsaboutthe(recentlycreated) inertialYaxis,then
the inertial Zaxis is as follows:
Pfinal = R,Z · R,Y · (P0.CoM Pfacility.CoM +(R ,Z · R,Y · P1)) (4.1)
As inSection3.3,Lagrange’sequationswillnowbederivedforthesystem including
the orbital inclination.
4.1.1 5 degree of freedom model
Five variables of interest are chosen – {R,
,
,
,
}– thatfulfil all the requirements
of generalised coordinates. An analysis will be performed to investigate the effect
the inclination variable has on the tether system orbit.
The mathematical model will be constructed to include a symmetrical dumbbell
tether with a facility mass at the centre of the system. The important mass points
included in the model are the two identical payload masses, Mpayload, the facility
mass, Mfacility, and the tether masses, Mtether, which are concentrated at the centre
of each tether.
The orbital system {R,
,
}is conceptually separated from the local rotation system
{ ,
}. The equations of motion for the orbit are constructed as one concentrated
point mass, and the local rotational system is superimposed on this orbit.
This allows a significantly simpler representation of the equations, however, this
does not allow the orbital motion to be influenced by the local rotations. Overall,
this assumption is valid and should give an insight into the broad motions of the
MMET on an inclined orbit.
66
The stator is not included in this model because this would significantly increase
the complexity of the system.
4.2 Constructing the equations of motion
The sequence of rotations follows the YZ system of rotations as described earlier in
Section 3.5.2.
Manipulating the position vector from the original Equation 4.1, the rotation
matrices are collected through the commutative properties of matrices and terms
containing the Radius vector, P0.CoM, and the local vector3 , P1, are separated:
Pfinal = R,Z · R,Y · (P0.CoM Pfacility.CoM +(R ,Z · R,Y · P1)) (4.2)
Pfinal = R,Z · R,Y · (P0.CoM Pfacility.CoM)+R,Z · R,Y · R ,Z · R,Y · P1 (4.3)
Taking the assumption that the dumbbell is symmetric about the facility mass:
.
.
.
.
CoMX 0
.
.
.
.
.
.
.
.
.
...
Pfacility.CoM =0 . CoMY =0
.
...
.
...
CoMZ 0
gives the simplified positional equation Pfinal for a single mass point, P1 in the
inertial coordinate system:
Pfinal = R,Z · R,Y · (P0.CoM)+R,Z · R,Y · R ,Z · R,Y · P1 (4.4)
3The local vector is the position vector between the facility and payload 1.
67
When evaluated, this becomes:
...
.
cos .
cos.
cos(cos.
cos .
cos.
sin.
sin )cos.
sina
sin.
...
.
...
.
...
.
Pfinal = R
cos.
sin .
+ L
cos(cos.
cos.
sin.
+ cos .
sin )sina
sin .
sin .
.
.
sin.
.
.
.
.
cos .
sina
cosa
cos.
sin.
.
.
(4.5)
The purpose of separating the inertial and local systems in Equation 4.5 is that
Lagrange’s equations are easier to derive and subsequently, the equations of motion
are computationally more efficient to solve.
4.2.1 Rotations in the local coordinate system
Figures 4.2a, 4.2b, 4.3a and 4.3b show the sequence of rotations needed to assemble
Equation 4.4: R,Z · R,Y · R ,Z · R,Y · P1.
Figures 4.2a and 4.2b show the steps needed to assemble the local coordinate
system, with Figure 4.2a duplicating the steps followed to arrive at the YZ rotation
sequence followed in Section 3.5.2. This is followed by a translation through
P0.CoM = {R,
0,
0}to transfer the local coordinate system to the inertial coordinate
system. Figures 4.3a and 4.3b show the tether rotating around the inertial axes
system – first through a rotation of .
around the inertial Yaxis, then a rotation of
.
around the inertial Zaxis.
Figure Rotational sequence
Figure 4.2a R ,Z · R,Y · P1
Figure 4.2b R+ R ,Z · R,Y · P1
Figure 4.3a R,Y · (R+ R ,Z · R,Y · P1)
Figure 4.3b R,Z · R,Y · (R+ R ,Z · R,Y · P1)
Table 4.1: Reference table with figure and corresponding rotation
68
Z
y
a
Y
y
a
a
ZYZYZXZXYZXYZY(a)YZ local rotation – R,Ythen R ,Z
X
y
y
q
a
i
i
a
a
Y
X
Y
Z
Z
X
XYZYZZYZXZXYZY
(b)Translation through R along Xaxis
Figure 4.2: Rotation sequence for YZ rotation with inclination. a
is shown as a
negative rotation here for ease of visualisation.
69
y
Z
Y
X
i
a
y
a
y
y
q
a
i
i
a
a
Y
X
Y
Z
Z
X
ZYXYZYZYZXZXYZXYZYXYZYZZYZXZXYZY
y
Z
Y
X
i
a
y
a
y
y
q
a
i
i
a
a
Y
X
Y
Z
Z
X
ZYXYZYZYZXZXYZXYZYXYZYZZYZXZXYZY
(a)Y inertial rotation through angle .
– R,Y
X
Y
Z
q
a
a
y
ay
y
Z
Y
X
i
a
y
a
y
y
q
a
i
i
a
a
Y
X
Y
Z
Z
X
XYZZYZYZXZXYYZZYZYZXXZZYXYZYZYZXZXYZXYZYXYZYZZYZXZXYZY
(b)Z inertial rotation through angle .
– R,Z
Figure 4.3: Rotation sequence for YZ rotation with inclination. a
and .
are shown
as negative rotations here for ease of visualisation.
70
4.2.2 Lagrange’s equations
Oncetheposition of the masspointshavebeendefined,theprocess of constructing
the equations of motion will follow that outlined in Section 3.3.
As the equations form two separate groups, this will be examined separately: in
Section 4.2.3 for the local system and Section 4.2.4 for the orbital system. The
equations are combined and listed in Appendix A.
4.2.3 Lagrange’s equations – local system
The mass points, as indicated earlier, are:
.
the two payload masses, each a distance L
from the facility mass and each
with mass Mpayload
.
the two tethers with masspoints each adistance L/2,halfway alongthetether,
each with mass Mtether
.
the facility mass, Mfacility.
Thetethersystem issymmetrical;theinplaneangleswillbedefined as .
and  ,
while the outofplane angles will be defined as a
and a
respectively.
Aspreviously,thetethersareassumed toberigidand inextensibletosimplify the
equations. The stator will not be included to simplify the equations. The masses of
the payloads are assumed to be equal and the masses of the tethers are assumed to
be equal such that:
Mpayload 1 = Mpayload 2 = Mpayload and Mtether 1 = Mtether 2 = Mtether
71
The locations of each of the mass points in the inertial plane are therefore:
...
.
cos.
cos.
cos(cos.
cos.
cos.
sin.
sin )cos.
sina
sin.
...
.
...
.
...
.
Ppayload
1
= R
cos.
sin.
+ L
cos(cos.
cos.
sin.
+ cos.
sin )sina
sin.
sin.
(4.6)
...
.
...
.
sin.
cos.
sina
cosa
cos.
sin.
...
.
cos.
cos.
cos(cos.
cos.
cos.
sin.
sin )cos.
sina
sin.
...
.
...
.
...
.
Ppayload
2
= R
cos.
sin.
L
cos(cos.
cos.
sin.
+ cos.
sin )sina
sin.
sin.
(4.7)
...
.
...
.
sin.
cos.
sina
cosa
cos.
sin.
...
.
cos.
cos.
cos(cos.
cos.
cos.
sin.
sin )cos.
sina
sin.
...
.
..
L
..
...
.
Ptether
1
= R
cos.
sin.
+ cos(cos.
cos.
sin.
+ cos.
sin )sina
sin.
sin.
(4.8)
...
.
2
...
.
sin.
cos.
sina
cosa
cos.
sin.
...
.
cos.
cos.
cos(cos.
cos.
cos.
sin.
sin )cos.
sina
sin.
...
.
..
L
..
...
.
Ptether
2
= R
cos.
sin.
 cos(cos.
cos.
sin.
+ cos.
sin )sina
sin.
sin.
(4.9)
...
.
2
...
.
sin.
cos.
sina
cosa
cos.
sin.
..
cos.
cos.
..
..
..
Pfacility
= R
cos.
sin.
(4.10)
..
..
sin.
The velocities of each of the mass points are then calculated. These are not listed
here because of the length of the equations.
The fivemasspointsareused to find thetranslationalkineticenergy, Ttrans, of the
system:
Mpayload 1 Mpayload 2
Ttrans = Vpayload 1 · Vpayload 1 + Vpayload 2 · Vpayload 2
22
Mtether 1 Mtether 2
+ Vtether 1 · Vtether 1 + Vtether 2 · Vtether 2
22
1
+(Mpayload 1 + Mpayload 2 + Mtether 1 + Mtether 2 + Mfacility)VCoM · VCoM (4.11)
2
72
which evaluates to:
1
Ttrans = cos 2().2 +.2 R2 + R.
2 (MFacility +2(Mpayload + Mtether))+
2
1
...
L2 16.2 +32.
sin()sin( ).a
4 2cos(2 )cos2()+cos(2)3 .2+
64
16cos2() .
2 +16.
.
2cos2()cos()cos( )sin(2)sin()+
.2 8cos2()cos(2 )sin2()+2cos(2)2cos(2)+
3cos(2(.
))+3cos(2(a
+ ))8cos( )sin(2)sin(2)+10)+
.
.
.
16..sin()sin(2 )cos2()+2cos( ).a
+ cos().+ .
sin(2)sin( )
(4Mpayload + Mtether) (4.12)
The rotational kinetic energy, Trot, of the system is found in a similar way:
Trot =(Ipayload 1 + Ipayload 2 + Itether 1 + Itether 2)· (!local · !local) (4.13)
where:
..
..
0 100
..
..
..
..
..
..
!local = .I3 = 010
..
..
..
..
.
.
001
Ipayload 1 =L2Mpayload 1 · I3 Ipayload 2 =L2Mpayload 2 · I3
Itether 1 =
1
L2Mtether 1 · I3 Itether 2 =
1
L2Mtether 2 · I3
33
The evaluated sum of the component rotational kinetic energies are:
1
Trot = L2 (3Mpayload + Mtether) .2 + .
2 (4.14)
3
73
The sum of the potential energies, U, are:
nn
µmj
µmj
U
=
.
=
.
(4.15)
Pj
Pj
· Pj
j=1 j=1
where Pj
are the magnitude of the mass point locations, given in Equations 4.6 to
4.10.
The full rotational sequence is not needed here, merely the distance between the
two points, which leads to a simpler expression of the magnitude of the distance.
Taking Equation4.2asthestartpoint,thedistance(i.e.,themagnitude) is:
.
Ppayload 1
.
=
.
R,Z · R,Y · (P0.CoM + R ,Z · R,Y · P1)
.
(4.16)
.
Ppayload 1
.
=
.
R,Z
.
· .
R,Y
.
· .
P0.CoM + R ,Z · R,Y · P1
.
(4.17)
The magnitude of a purely rotational matrix is, by definition, equal to one:
.
R,Z
.
=1 (4.18)
.
R,Y
.
=1 (4.19)
This simplifies Equation 4.17 to:
.
Ppayload 1
.
=
.
P0.CoM + R ,Z · R,Y · P1
.
(4.20)
...
.
.
R
L
cosa
cos.
...
..
...
..
...
..
.
Ppayload 1
.
=0+ L
cosa
sin.
(4.21)
...
..
...
.
.
0 L
sin a
.
Ppayload 1 = L2 +2LR
cosa
cos.
+ R2 (4.22)
74
The evaluated expression for the potential energy is therefore:
U
=  .
µ
Mfacility .
µ
.
v
1
+ v
1 .
Mpayload
R
L2 Lmed + R2 L2 + Lmed + R2
11
2µ
v + v Mtether
L2 2Lmed +4R2 L2 +2Lmed +4R2
(4.23)
where:
Lmed =2LR
cosa
cos.
The Lagrangian is then given by L
= T
U:
2 2
L2 (3Mpayload + Mtether) .+.
L
= Ttrans +
3
µMfacility 11
++ µ
v + v Mpayload
R
L2 Lmed + R2 L2 + Lmed + R2
11
+2µ
v + v Mtether (4.24)
L2 2Lmed +4R2 L2 +2Lmed +4R2
where Ttrans is given by Equation 4.12.
Lagrange’s equations are then generated for two out of the five generalised co
ordinates, { ,
}as follows:
d
@T
@T
@U
 += Qj
(4.25)
dt
@q.j
@qj
@qj
and are provided in Appendix A along with the nonconservative forcing terms.
4.2.4 Lagrange’s equations – global system
The global equations used to construct the motion of the facility around the Earth
are now outlined. The mass for the system is summed and used as the total mass
75
located at the CoM:
Mtotal =2Mpayload +2Mtether + Mfacility (4.26)
The locationsof thesinglerepresentativemasspoint inthe inertialplane istherefore:
..
R
cos.
cos .
..
..
..
Porbit = R
cos.
sin.
(4.27)
..
..
R
sin .
The velocities ofeach of the masspoints are then calculated4 . Equation4.28 shows
the velocity of the system CoM:
..
.
R
cos.
cos.
R.
cos.
sin.
+..
cos.
sin.
..
..
..
Vorbit = R.cos.
cos.
+ sinR.
cos.
R.sin.
(4.28)
..
..
R.cos.
R.
sin .
Thesinglemasspoint isused to find thetranslationalkineticenergy, Ttrans, of the
system:
Mtotal
Ttrans =(Vorbit · Vorbit) (4.29)
2
Mtotal
Ttrans = .2 cos 2()+.2 R2 + R.
2 (4.30)
2
The rotational kinetic energy is assumed to be zero:
Trot =0 (4.31)
4Note: as the inclination term is given the symbol iota, , the first derivative of the inclination
term, , is., and not the ninth letter of the Latin alphabet, i.
76
The expression for potential energy, U, is simply:
µMtotal µMtotal
U
= = (4.32)
R.R R
as the distance to the CoM is defined as R.
The Lagrangian is then given by L
= T
U:
Mtotal µMtotal
L
= .2 cos 2()+.2 R2 + R.
2  (4.33)
2 R
Lagrange’s equations are then generated for three out of the the five generalised
coordinates, {R,
,
}as follows:
.
.
d
@T
@T
@U
dt
@q.j
 @qj
+
@qj
= Qj
(4.34)
and are provided in Appendix A along with the nonconservative forcing terms.
4.2.5 Nonconservative forces
TheMMET system requires that aforcebe transferredfrom the motor to the tether
arms to ‘spinup’ the system. The mechanics of spinup of such a system is nontrivial,
especially when flexural5 terms and inertial terms are taken into account.
The spinup dynamics of length deployment is covered in Chapter 5.
Thetetherarmsareassumed toberigidand symmetrical at alltimes, limiting the
payload tethers to a planar disc. The motor torque is assumed to act in the same
plane as this disc. The stator arms are neglected.
The motor torque will exert a large power drain, therefore all the power for the
motor willbederivedfrom the solar arraysharnessing afree,plentiful and renewable
source of energy. The Forcing term in Lagrange’s equations must reflect the binary
5Flexural terms are not included in the analysis due to the complexities they introduce.
77
nature of the availability of the power supply. The Eclipse function modulates the
motor torque to take into account the eclipsing of the Sun by the Earth, reducing
the torque to t
= 0 when in darkness and taking the default t
= max when in
sunlight.
Eclipse is the binary function:
cos.
cos>
0
Eclipse =
.
.
.
.
.
.
.
.
.
.
.
.
1
.
.
.
.
.
.
.
R
cos .
2
.R2 = R2
.
Earth
.
.
RSun
.
.
.
.
.
.
.
.
1cos2()cos2() 1+
AND
(4.35)
.
.
.
.
.
.
.
.
.
.
.
.
0 otherwise
The radius of the Earth’s shadow is scaled from REarth by the factor RSun+R
cos .
RSun
as shown in Figure 4.4. The hypotenuse of the Yand Zcomponent of the CoM at
a time, t, in the orbit path is calculated and checked. If it is less than the shadow
distance, then the CoM is in shadow.
RSun + R
cos.
RShadow = REarth
RSun
.
(4.36)
R
cos.
RShadow = REarth 1+
RSun
Y
2 22
=(R
cos .
sin)+(R
sin )
CoM inertial + Z2
CoM inertial
(4.37)
= R2 (1cos2()cos2())
Hence, the eclipse checking function is calculated as:
2
...
R
cos()
R2 1cos2()cos2() = R2 1+ (4.38)
Earth
RSun
which gives two eclipses an orbit, therefore an additional constraint is added to
78
Y
Z
Sun
CoM Orbital Path
Cone of Eclipse
EarthR
Earth
P
X......................
....
..................
..............i q..............
..
....
RSun
Figure4.4:Eclipse checkingfunction(notto scale)
ensure that the CoM is eclipsed only once per orbit:
cos.
cos >
0 (4.39)
A diagram of the distances onorbit is shown in Figure 4.4 and the Mathematica
output for a 29. orbit is shown in Figure 4.5. The Sun is assumed to be a point
source and the Earth is assumed to describe a circular orbit around the sun.
Eclipse
EclipseEclipse1
11
0
00
Orbits
OrbitsOrbits
1
113
33
0
001
112
22
..
.....
...
....
.....
...
..
2
222
22
Figure 4.5: Sample plot of the Eclipse function in Mathematica
79
4.3 Analysis of tether on an inclined orbit
The generated equations are solved with Mathematica’s general purpose solver
NDSolve[. . .] with the following parameters:
AccuracyGoal >
Automatic
PrecisionGoal >
Automatic
WorkingPrecision >
20
To fully understand the interaction between the local and inertial systems, Lagrange’s
equations are solved with different initial conditions or ‘cases’. These are
given in Table 4.2.
.
case . ..
a
eccent t
I1 0 0 0 0 0 0 0 1000
I2 0 0 0 29/180 0 0 0 1000
I3 0 1 0 29/180 0 0 0 1000
I4 0 0 0 29/180 0 29/180 0 1000
I5 0 1 0 29/180 0 29/180 0 1000
Table 4.2: Initial conditions for the inclination numerical solutions.
Case I1 gives a baseline in the inertial plane: no inclination, no outofplane spin
and no initial inplane spin. Case I2 introduces inclination, while case I3 has an
inclined orbitwith asignificant initial inplanespin. CasesI4 andI5 arecomparable
tocasesI2 andI3,with thedifferenceof anoutofplanespinequal tothe inclination
(i.e. so the spinplane is coincidentalwith the inertial plane).
All five cases have an equal torque level to enable a fair comparison.
The initial conditions6 for the other parameters held constant are:
L
=1000m,
Mpayload =10kg,
Mtether =100kg,
Mfacility =500kg,
e
=0,
R.
[0] =0,.[0] =0.00113905rad/s,
[0] =0,R[0] =7.378* 106 m
6Note: the inclination term is represented by the Greek letter iota – .
– therefore the first
derivative with respect to time is iota dot – . – and the second derivative iota dot dot – ¨. This
.
was defined in order to remove the confusion between i dot and iota dot – i,
..
80
Figures4.11 to4.18 show the solution to the equations of motionfor10 orbits over
63069.5 seconds where:
.
Figure 4.11 shows plots of a
vs. time
.
Figure 4.12 shows plots of .a
vs. time
.
Figure 4.13 shows plots of a
vs. .
.
Figure 4.14 shows 3D phase plots of .
.
Figure 4.15 shows 3D phase plots of a
.
Figure 4.16 shows plots of V
vs. time
.
Figure 4.17 shows plots of tension vs. time
.
Figure 4.18 shows plots of stress vs. time
4.3.0.1 Behaviour of R, .
and .
with time
The orbital parameters R, .
and .
are shown in Figure 4.10. They are deliberately
chosensothattheorbitalparametersareindependentfromthe localparameters,as
such the MMET does not influence any of these generalised coordinates.
Together, the three generalised coordinates R, .
and .
describe a circular orbit in
the inertial plane in case I1 and a circular inclined orbit in cases I2 to I5.
4.3.0.2 Behaviour of a
with time
Figures 4.11a to 4.11e demonstrate the a
vs. time behaviour of the system. Clearly,
the baseline case I1 shows there is no interaction between the inclination and out
ofplane coordinates when both are initially zero.
Comparing CasesI2 andI3 inFigures4.11b and4.11cgivesan ideaof therelative
.
stabilising influence of the inplane spin velocity, . It is the gyroscopic action of
the inplane spin that limits the outofplane coordinate, .
The opposite appears to be true when Figures 4.11d and 4.11e – cases I4 and I5
81
– are compared. The outofplane spin angle is set such that the local system is
aligned with the inertialplane. Whenthis isacceleratedfromzerovelocity initially,
the outofplane angle returns to be coplanar with the inclination orbit7, with a
small nutation of(0.1 rad
˜ 6.). When started from a large initial inplane velocity
in the inertial plane, the opposite is true – it remains spinning in the inertial plane
with a small adjustment towards the orbital plane.
4.3.0.3 Behaviour of .with time
Figures 4.12a to 4.12e demonstrate the .vs. time behaviour of the system. As
previously, the baseline case I1shows there is no interaction between the inclination
and outofplane coordinates when both are initially zero.
Case I2 and I4 show the angular rate increasing from zero – this is larger than
desired for operating the MMET, especially in the early startup regime. However,
case I3 shows that the greater the inplane spin velocity, the lower the outofplane
nutation.
4.3.0.4 Behaviour of a
with .
Figures 4.13a to 4.13e show the .a
vs. .
interaction of the system.
The graphs here are inconclusive. There is a large interplay between a
and .
.
in cases I3 and I5, where the inplane velocity .
is large, and a correspondingly
.
smaller relationship between the two in cases I2 and I4 where the .
term is zero.
This, however, does not point to any concrete notions of interdependence.
CaseI4 clearly showstheprogression of a
froma largeoutofplanemotion(starting
in the bottom right of the graph) to a smaller angle as the inplane velocity
increases.
7The inclination orbit is the orbital plane, offset 29/180rad
from the inertial plane.
82
.¨
4.3.0.5 Behaviour of ,
,
.
The eclipse functionality is clearly shown in the .
phase plots in Figires 4.14a to
4.14e and the plot of V
in Figure 4.16a to 4.16e.
The .
phase plane plots illustrate the shadow periods of the graph where the in
plane angular acceleration ¨, shown in the horizontal axis, is approximately zero
compared to the illuminated periods where the tether is accelerated by the motor.
Theplots ofV
highlight the end result of the work done by the motor – the V
thesystemcan imparttothepayload. Inthiscase,aV
of approximately 400m/s
is attained for a modest acceleration period of 17.5 hours.
In cases I1 and I3, the outofplane coordinate a
is zero and the rate of increase
.
in .
is highest. The greater the outofplane spin angle, the less energy is available
totransfer into increasing the(useful) inplanespinrate.
Cases I4 and I5 show the more disordered plots when the outofplane coordinate
is nonzero. Clearly, having an clean acceleration profile is advantageous from a
stressminimisation point of view – these profiles have a spinup phase that feature
a large .a
term, which is difficult to aim when the payload must be released.
Similarly, case I2 does not have a clean spinup phase when compared to case I1.
This ispurelyduetothe largeoutofplanemotionthat interfereswith the inplane
spin.
4.3.0.6 Behaviour of ,
,
.¨
The a
phase plane behaviour of the system is shown in Figures 4.15a to 4.15e.
In case I3 – high .
and the outofplane initial angle, , set to zero – the out of
plane velocity is minimised. Marginally larger outofplane angular velocities and
accelerations are shown by case I4: zero initial .
and the outofplane initial angle
set to coincide with the inertial plane.
83
4.3.0.7 Behaviour of V
with time
The V
vs. time behaviour of the system is shown in Figures 4.16a to 4.16e.
The calculation ofV
isstraightforward, it istheproduct of thetether length and
the angular velocity. The maximum useful V
is obtained when the outofplane
angle is zero, and is analysed in Section 4.4. In this case, the maximum useful V
is simply the product of length and angular speed:
.
V
= L.
(4.40)
As thegraphs of V
show the magnitude of the resultant velocity, cases I1 I2 and
I4 should be almost identical, and cases I3 and I5 should also be identical. Whereas
the former is true, case I5 has a useful velocity that is lower than expected. This
is caused by rapid changes in the outofplane velocity ., which detracts from the
.
inplane spin velocity, .
4.3.0.8 Behaviour of tension with time
The tension profile of the tether is shown in Figures 4.17a to 4.17e.
Thecaseswith zeroinitialinplanespin – casesI1,I2 andI4 –exhibitsimilarlevels
of tension in the tether, approximately 10 kN. The tension appears to be slightly
lower in the baseline case I1, however, this may be due to the larger outofplane
movements in cases I2 and I4.
The remaining cases I3 and I5 start with a higher tension due to their large initial
inplane spin rate. Similarly to the V
plot in Figure 4.16e, the tension in case I5
rapidly fluctuates due to the large outofplane angular motions.
4.3.0.9 Behaviour of stress with time
The stress profile of the tether is shown in Figures 4.18a to 4.18e.
84
Thegraphsherearevery similartothetensiongraphs,asstress isequal totension
divided by the crosssectional area, 0.000064 m2 (equal to an 8mm
square section
of tether).
4.3.1 Analysis of tension and stress in tether
The tension in the tether ignores any outofplane motion of the tether, as it was
shown that the outofplane motion is not significant for small deviations from the
spinplane.
The tension in the tether is calculated as:
Lpayload
. 2
Tpayload = Mpayload 2Lpayload + A
.(4.41)
2
with the firsttermasthetensionduetothepayload,thesecond termasthetension
duetomassof thetether itself. Thesecond termdominates inthiscase,asthemass
of the tether is an order of magnitude higher than the mass of the payload.
The stress in the tether is calculated as:
Tpayload
s
= (4.42)
A
Themaximumallowablestress inthetether iscalculatedusingtheUltimateTensile
Stress(forZylon,theUTS is5.9* 109 N/m2)and aSafety Factor(SF=1.5) as:
s
5.9* 109
max
== =3.933GP
a
(4.43)
SF 1.5
The variables are defined as:
Mpayload =10kg, L
=1000m, .
=1570kg/m3 , A
=0.000064m2
using values ofthetether8 composed ofZylon.
fibres withpropertiesfromTableD.1,
8The area is equivalent to a rectangular extrusion of 0.008m.
85
Appendix D.
Themaximumallowablestressgivesanabsoluteceilingtothestresses inthetether.
This maximum stress, working backwards through the tension, is dependant on the
inplane velocity assuming the material and geometric properties of the system do
not change when it is inorbit. This interplay between .
and Lpayload
drives the
maximum V
and therefore the entire feasibility of the system.
4.4 Tether release aiming
The criteria to maximise the V
when using the motorised tether with a Hohmann
twoburn trajectory are as follows.
The tether will gain the most velocity from the release if the velocity vector is at
rightangles to the orbital plane9 . This is generally satisfied when the tether arm is
aligned with the radius vector, a criteria of .
.
0 and a
.
0
An excellent treatment of the sensitivity to errors in transfer orbits for Hohmann
transfers is given in [Kamel and Soliman, 2006]. This may be directly applied to
tethers for the same type of manoeuvres. For example, a 0.1% error in V
may
give up to a 5% error in radius and semimajor axis. A thorough treatment of the
errors due to an early or late release is needed for the tether.
Intermsofreleasestrategyforan inclined orbit,thiscriteriaforrelease isgenerally
only satisfied at one point in the cycle when the inclined orbital plane cuts the
inertialplane,thedirection isdeterminedfromtheclassicalHohmanntrajectory,as
in [McLaughlin, 2000].
An alternative to the Hohmann trajectory is the WSB method, which has been
describedinthe literature:[Belbruno,1987],[MillerandBelbruno,1991] and[Koon
9Thisassumesthatthemission istoprovidethemaximumV
;many alternativetether missions
exist, including a inclinationchangeusing atetherreleasing thepayload atanangletotheorbital
plane.
86
et al., 2001]. To summarise, the WSB utilises the gravitational interaction of the
EarthSunMoon system to lower the V
required to capture a satellite in Lunar
orbit. The WSB can be thought of as two Circular Restricted 3 Body Problems
(CR3BP)patched together with a small correction burn of less than 25m/s
at the
boundary between the SunEarth system and the EarthMoon system.
A patch point is located by integrating forward the dynamics of the spacecraft
in the EarthMoon system from the patch point until the spacecraft is ballistically
captured in Lunar orbit, and simultaneously integrating backward the dynamics of
the spacecraft in the SunEarth system until the spacecraft links up with a LEO
orbit around the Earth. Thus a patch point is found along with two halves of the
trajectory from LEO to ballistic capture in Lunar orbit via the patch point.
Figure4.6: Diagram showing an example of the WSB trajectory, fromESABulletin
103: [Biesbroek and Janin, 2000]
Timing the release of the tether payload to launch on a WSB transfer is a difficult
computational puzzle, but one which may be solved by optimising a numerical
solution to the equations of motion. A trajectory is planned for a future time and
the tether mission must provide a V
precisely with a very small margin of error
in release angle.
The useful tip velocity(V
)willbeapproximatelyequal tothe localbodycentred
87
velocity of the tip payload. The local velocity is given in Equation 4.44 below:
..
.
L
sin()cos( ).L
cos()sin( ).
..
..
Vpayload 1 = .L
sin()sin( ).+L
cos()cos( ) .
(4.44)
.
.
.
.
L
cos().a
This reduces to the simple onedimensional angular velocity when considering the
optimal release case of a perfectly aligned tether of .
. 0, a
. 0 and .a
. 0:
..
0
..
..
..
Vpayload 1 . L .
(4.45)
..
..
0
The major longterm perturbations such as J2, Lunar and Solar gravity etc. are
not covered here as they do not affect the shortterm dynamics of tether spinup.
[Parsons, 2006] covers some of these alongside the WSB trajectories.
4.4.1 Tether error analysis
Small changes in the tether variables could have a large effect on the performance
of the tether, therefore it is important to quantify these for the purposes of risk
mitigation. Ananalysiswasperformed toquantify theeffectsthat asmallchange in
thesystemparametershave onthepayload apogeeand inclination, whencalculated
at the orbital location directly opposite the separation point10 .
The Figures 4.7 and 4.8 are calculated by solution of the equations of motion for
MMET tether system with two payloads, with the default initial conditions:
Lpay
=1000m,
Mpay
=10kg,
Mtether
=100kg,
Mhub
=500kg,
e
=0,
..
.[0] =0rad/s,
[0] =0.1rad/s,
.[0] =0rad/s,
R[0] =0m/s,
10Calculatedas p
radians
roundthe orbitfrom the system centre of masspoint when thepayload
is separated from the tether.
88
.[0] =0.00113905rad/s,
[0] =29/180rad,
[0] =0rad,
[0] =0rad,
R[0] =0m,
[0] =0rad,
Torque =1000Nm,
NominalOrbitalPeriod =6306.95s.
The initial conditions above equate to an equivalent V
of 100m/s.
Thepayload isreleased whenthetether isaligned withtheradiusvector(.
=0
and a
= 0) and follows a ballistic trajectory thereafter. No allowance is made here
for angular momentum effects between the payload and tether system on release.
Errors are simulated by the addition of small offsets in the initial conditions at
the release point. The effect that the errors have on that variable is analysed by
comparing the apogee and the inclination to the unperturbed values measured at the
apogee point. The starting radius11 is 7,
378,
000m. With the default initial conditions
above, the apogee is 7,
747,
377m, a difference of 369,
377m
or approximately
5% of the radius vector. In terms of altitude, the starting altitude is 100,
000m
and
theapogeewhenthepayloadisreleased is469,
377m. The inclination of the apogee
point is equal to the magnitude of the perigee inclination of 29. .
Figure4.7ashowstheeffectthatasmallchange inradiushasonthetetherapogee.
A change in radius of 1000m
at the tether release point corresponds to a 7460m
increase intheapogeeof thepayload. Thedifferenceoftheapogee ishigherthanthe
expected first order approximation of 7L
as suggested12 in [Ziegler and Cartmell,
2001]by 460m, the higher calculated number is a result of the V
supplied by the
rotating tether.
The change in apogee due to small changes in inclination is shown in Figure 4.7b.
It follows from orbital mechanics that a small change in the inclination to flatten
the orbit will increase the apogee.
A small change in the inplane rotation angle, , causes minor changes to the apogee
altitude, as shown in Figure 4.7c. The relationship is predominantly parabolic
11That is, the perigee point.
12For a hanging tether.
89
based due to the increase in energy supplied to the orbit 12mV
2, which is itself
dependent on .
.
4.4.2 Lunar capture
Once the spacecraft is ballistically captured around the Moon, a small thrust is
required to ensure that the orbital energy is reduced below the threshold required
for the spacecraft to escape. A V
of less than 50m/s
will fix the spacecraft in a
highly elliptical(e
=0.9) but stable Lunar orbit.
A midpoint burn of V
= 25m/s
with a capture burn of V
= 50m/s
corresponds
to a fuel mass fraction of 9.6% assuming an Isp
rating of 75s
for a Nitrogen
coldgas thruster. Aconservative estimate of1kg
massfor the engine and associated
structure leavesthe10kg
original satellite mass with8.04kg
of usefulpayload mass.
When this is compared to the Hohmann transfer V
cost, [Koon et al., 2001,p71]
calculates that the WSB transfer requires 20% less fuel to execute. This necessarily
transfers into a lower mass of the payload and subsequent cost savings.
Once the payload is captured, it is possible to use tethers for the descent to the
Lunarsurface. Staged tethersas in[Cartmellet al.,2004]orLunarelevatorsmaybe
used to transfer the payload to the surface, and Lunar minerals to LEO. [Williams
et al., 2005, p 784] claims that the payload can be delivered within 1cm
for a
controlled handover of a staged tether system, and can be held in a 10m
window
for approximately 135 seconds, which is ideal for the multistaged tether system
provided a control system can be devised that will enable an automatic remote
handover of payloads.
90
4.5
A demonstration Lunar mission with a tether
launch
A demonstration mission to the Moon is outlined below, involving a tether launch
fromLEO toLunarparking orbit. It isexamined using hardwareavailable in2008,
as a lowcost microsatellite mission that could be launched in the short term, with
today’s technology.
4.5.1
Spacecraft sizing
The MMET system is sized to launch a 10kg
payload to Lunar capture orbit. This
is intended as a demonstration level mission, although the MMET system may be
scaled to launch larger payloads, as discussed previously.
The facility is proposed to be entirely housed within a cylinder containing the
motor andgearbox. The solar arrays will consist of35000 standard sized2cmx4cm
solar cells, whichprovide thepower requiredbythe motor, rechargethebatteries and
providepower to the spacecraftbus. The solar cells willbe aInGaP/GaAs/Ge
type
capable of 25% efficiency [Green et al., 2003], with a degradation rate of 1%/year
[GriffinandFrench,1991].Themission length of5years isbased ontheprobability
of a tether cut due to orbital debris when using multiple redundant tether strands,
coupled with the useful life of the solar arrays.
Thepowertothemotorwillonlybeavailableduring thetimethefacility isnot in
eclipse and this is factored into the equations using the function Eclipse described
in Section 4.2.5.
The totalusefulpower availablefor a2mx3m
cylinder augmentedby a1.5mx6m
fold out panel will be 4.3kW
when launched and 3.05kW
atthe end of a5year
mission. The end of life power allows the motor characteristics to be defined, in
this case the motor used will be a 3kW
rated motor capable of driving 1000Nm
91
of torque through a reducing gearbox. The gearbox will reduce the rotations from
the motor armature rotation of 430rad/s
to the required angular velocity of the
rotor arm of 20rad/s. This motor is based around GE Industrial DC motor model
5BC49JB1115, rated at 4hp
and 4150RP
M, which is designed to drive an electric
vehicle.
The rotor is comprised of two 1km
tether sections with each linking the facility
to the payloads. The stator is identical, with the payloads replaced by counterweights.
The tether itself is Zylon.
fibre with a tensile strength of 5.8GP
a, density
of 1570kg/m3 and cross sectional area of 64mm2 . This leads to a tether mass of
100kg
for each of the four tether lines, totalling 400kg. A safety factor of 1.3 is
applied to the calculation of the maximum breaking stress, however additional precautions
must be taken by careful design of the tether structure to ensure the tether
will remain intact if impacted by micrometeorites or debris [Forward and Hoyt,
July 1995].
The tether mass outweighs the payload mass in order to keep the breaking stress
to a manageable level. This is advantageous as the tethers may be adjusted to
compensate for any tether failure mode or asymmetric payload release.
The MMET is designed to be highly reusable. The 5 day spin up cycle must be
matched by an equal 5 day spin down cycle to return the rotor to a stationary
positionaligned with thegravity vector. Whenthis isachieved,thepayload willbe
much easier to dock than when in motion. For the WSB transfer, payloads can be
transferred to the Moon every Lunar month, leaving a window of 20 days to dock
the payloads to the tethers.
92
4.5.2 Analysis
The square of the characteristic velocity, Uvel, of the tether is proportional to the
tensile strength, T, and inversely proportional to the density, :
2T
Uvel = (4.46)
.
The characteristic velocitydeterminesthe maximum velocityattainablefor atether
supporting its own mass while rotating. Adding a payload limits the velocity a
tether can impart by increasing the tension and stress in the tether, therefore by
keeping the mass of the payload to a fraction of the tether’s mass, the decrease in
characteristic velocity is minimised.
A tether manufactured from Zylon.
with T
=5.9GP
a, .
= 1570kg/m3 with a
safety factor of 1.3 gives a characteristic velocity of 2.238km/s. With a payload
mass equivalent to 10% of the tether mass attatched to the end of the tether, the
characteristic velocity will drop by approximately 10% to 2.10km/s. To achieve
Lunar orbit, thepayloadwill need3.187km/s
[Sweetster,1991]fromLEO at167km
altitude. Therefore, the tether will have to occupy a higher earth orbit of 1000km
to successfully launch to the moon, or a staged system will be required [Cartmell
et al., 2004].
The necessary orbital altitude of 1000km
is unfortunate in terms of the orbital
debris that the tether will encounter. Figure 4.9 shows the debris environment from
LEO to Geostationary Earth Orbit (GEO) – the debris environment is relatively
dense until 3000km
altitude, meaning the risk of a severed tether is higher in this
region of space. This is compounded by the long time to clear the debris from the
higher orbits, meaning the orbitaldebris will continuetoposeproblemsfordecades
to come.
93
4.6 Conclusions
A demonstration mission launching a microsatellite from MEO to Lunar orbit is
achievable using current technology. The safety margins, however, are extremely
small and the MMET launcher would be located in an orbit with large amounts of
debris.
The inclination term does not significantly alter the dynamics of the MMET system.
However,anyoutofplanecomponentinthelocal axes –whetheronaninclined
orbitornot –interactswith theorbitalparameterstocreateanunsteady oscillation
in the tether. This is unacceptable for payload aiming and release, therefore the
local outofplane angle should be minimised insofar as possible.
94
Apogee change @m.
Apogee change @%.
1000 500 500 1000
6000
4000
2000
2000
4000
6000
Rad
0.01 0.005 0.005 0.01
0.01
0.005
0.005
0.01
Inc
change
change
@%D
@m.
(a) Payload apogee change (m) vs. change (b) Payload apogee change (%) vs. change
in radius (m). Solid line – calculated in inclination (%). Solid line – calcuresults,
dotted line – linear relationship lated results, dotted line – linear relationship
Apogee
=7* R
Apogee
= 1* .
1.5 1 0.5 0.5 1 1.5
175
150
125
100
75
50
25
Apogee change @m.
Apogee change @m.
Psi
@deg.
0.01 0.005 0.005 0.01
30000
20000
10000
10000
20000
30000
PsiDot
@rads.
(c) Payload apogee change (m) vs. change (d) Payload apogee change (m) vs. change
.
in .
(deg). Solid line – calculated resin .
(rad/s). Solid line – calculated
ults, dotted line – trigonometric relationship results, dotted line – linear relationship
Apogee
=(( ))2
*60 Apogee
=3500000* .
Figure 4.7: Change in apogee of the payload after release from the tether system
after one half of an orbit compared to small changes of system parameters
95
Apogee change @m.
Apogee change @m.
1 0.5 0.5 1
3000
2000
1000
1000
2000
3000
10 5 5 10
4000
2000
2000
4000
Length
Alpha
Change
@Deg.
@m.
(a) Payload apogee change (m) vs. change (b) Payload apogee change (m) vs. change
in a
(deg). Solid line – calculated in length (m). Solid line – calculated
results, dotted line – linear relationship results, dotted line – linear relationship
Apogee
= 3000* a
Apogee
=400L
Inclination change @%.
1 0.5 0.5 1
Alpha
@degD
0.0004
0.0002
0.0002
0.0004
(c) Payload inclination change (%) vs.
change in a
(deg). Solid line – calculated
results, dotted line – linear relationship
.
=
 1
2000
Figure 4.8: Change in apogee and inclination of the payload after release from
the tether system after one half of an orbit compared to small changes of system
parameters
96
Figure 4.9: The orbital debris environment from LEO to GEO, from [Wikipedia,
2008]
97
4.7 Inclination plots
Table 4.3 shows the varying cases with an explanation of their initial conditions.
This is an exact duplicate of Table 4.2, reprinted for reference.
.
case . ..
a
eccent t
I1 0 0 0 0 0 0 0 1000
I2 0 0 0 29/180 0 0 0 1000
I3 0 1 0 29/180 0 0 0 1000
I4 0 0 0 29/180 0 29/180 0 1000
I5 0 1 0 29/180 0 29/180 0 1000
Table 4.3: Initial conditions for the inclination numerical solutions.
The initial conditions for the other parameters held constant are:
L
=1000m,
Mpayload =10kg,
Mtether =100kg,
Mfacility =500kg,
e
=0,
R.
[0] =0,.
[0] =0.00113905rad/s,
[0] =0,R[0] =7.378* 106 m
98
Orbit Radius
Orbit Radius
7
1.4·10
6
7
7.378·10
1.2·10
7
1·10
6
7.378·10
6
8·10
6
6
7.378·10
6·10
6
4·10
6
7.378·10
6
2·10
t
t
21600 43200
21600 43200
(a) Case 1: R
(m)vs.time(s) (b) Cases 25: R
(m)vs.time(s)
TT
60
60
50
50
40
40
30
30
20
20
10
10
t
t
21600 43200
21600 43200
(c) Case 1: .
(rad)vs.time(s) (d) Cases 25: .
(rad)vs.time(s)
inclination inclination
1
0.5
t
21600
43200
0.5
1
21600 43200
t
0.4
0.2
0.2
0.4
(e) Case 1: .
(rad)vs.time(s) (f) Cases 25: .
(rad)vs.time(s)
Figure 4.10: Plots of R, .
and .
vs. time for the five cases.
99
21600 43200
t
1
0.5
0.5
1
.
(a) Case 1: (rad) vs. time (s)
21600 43200
t
0.2
0.1
0.1
0.2
.
(b) Case 2: (rad) vs. time (s)
21600 43200
t
0.0015
0.001
0.0005
0.0005
0.001
0.0015
.
(c) Case 3: (rad) vs. time (s)
21600 43200
t
0.5
0.4
0.3
0.2
0.1
0.1
.
(d) Case 4: (rad) vs. time (s)
21600 43200
t
0.4
0.2
0.2
0.4
.
(e) Case 5: (rad) vs. time (s)
Figure 4.11: Plots of vs. time for the five cases.
100
21600 43200
t
1
0.5
0.5
1
D@.,tD
(a) Case 1: . (rad/s) vs. time (s)
21600 43200
t
0.01
0.005
0.005
0.01
D@.,tD
(b) Case 2: . (rad/s) vs. time (s)
21600 43200
t
0.0015
0.001
0.0005
0.0005
0.001
0.0015
D@.,tD
(c) Case 3: . (rad/s) vs. time (s)
21600 43200
t
0.004
0.002
0.002
0.004
D@.,tD
(d) Case 4: . (rad/s) vs. time (s)
21600 43200
t
0.3
0.2
0.1
0.1
0.2
0.3
D@.,tD
(e) Case 5: . (rad/s) vs. time (s)
Figure 4.12: Plots of . vs. time for the five cases.
101
1 0.5 0.5 1
i
1
0.5
0.5
1
.
(a) Case 1: (rad) vs. (rad)
0.4 0.2 0.2 0.4
i
0.2
0.1
0.1
0.2
.
(b) Case 2: (rad) vs. (rad)
0.4 0.2 0.2 0.4
i
0.0015
0.001
0.0005
0.0005
0.001
0.0015
.
(c) Case 3: (rad) vs. (rad)
0.4 0.2 0.2 0.4
i
0.5
0.4
0.3
0.2
0.1
0.1
.
(d) Case 4: (rad) vs. (rad)
0.4 0.2 0.2 0.4
i
0.4
0.2
0.2
0.4
.
(e) Case 5: (rad) vs. (rad)
Figure 4.13: Plots of vs. for the five cases.
102
0.00001
5·106
0
5·106
0.00001
.''
0
0.1
0.2
0.3
0.4
.'
0
5000
10000
.
(a) Case 1: , .
, ¨ phase plot
0.0001
0.00005
0
0.00005
0.0001
.''
0
0.1
0.2
0.3
0.4
.'
0
5000
10000
.
(b) Case 2: , .
, ¨ phase plot
0
5·106
0.00001
.''
1
1.1
1.2
1.3
1.4
.'
0
20000
40000
60000
.
(c) Case 3: , .
, ¨ phase plot
0.00002
0
.'' 0.00002
0
0.1
0.2
0.3
0.4
.'
0
5000
10000
.
(d) Case 4: , .
, ¨ phase plot
0.05
0
.'' 0.05
1
1.2
1.4
.'
0
20000
40000
60000
.
(e) Case 5: , .
, ¨ phase plot
Figure 4.14: Phase plots of for the five cases.
103
1
0.5
0
0.5
1
.''
1
0.5
0
0.5
1
.'
1
0.5
0
0.5
1
.
(a) Case 1: , . , ¨ phase plot
0.002
0
.'' 0.002
0.01
0
0.01
.'
0.2
0
0.2
.
(b) Case 2: , . , ¨ phase plot
0.001
0
.'' 0.001
0.001
0
.' 0.001
0.001
0
0.001
.
(c) Case 3: , . , ¨ phase plot
0.001
0
.'' 0.001
0.005
0.0025
0
0.0025
0.005
.'
0.4
0.2
0
.
(d) Case 4: , . , ¨ phase plot
0.2
0
.'' 0.2
0.2
0
0.2 .'
0.5
0.25
0
0.25
0.5
.
(e) Case 5: , . , ¨ phase plot
Figure 4.15: Phase plots of for the five cases.
104
21600 43200
t
100
200
300
400
DeltaV
(a) Case 1: V (m/s) vs. time (s)
21600 43200
t
100
200
300
400
DeltaV
(b) Case 2: V (m/s) vs. time (s)
21600 43200
t
1100
1200
1300
1400
DeltaV
(c) Case 3: V (m/s) vs. time (s)
21600 43200
t
100
200
300
400
DeltaV
(d) Case 4: V (m/s) vs. time (s)
21600 43200
t
1000
1100
1200
1300
1400
DeltaV
(e) Case 5: V (m/s) vs. time (s)
Figure 4.16: Plots of V vs. time for the five cases.
105
21600 43200
t
2000
4000
6000
8000
10000
TetherTension N
(a) Case 1: Tension (N) vs. time (s)
21600 43200
t
2000
4000
6000
8000
10000
TetherTension N
(b) Case 2: Tension (N) vs. time (s)
21600 43200
t
70000
80000
90000
100000
110000
120000
TetherTension N
(c) Case 3: Tension (N) vs. time (s)
21600 43200
t
2000
4000
6000
8000
10000
TetherTension N
(d) Case 4: Tension (N) vs. time (s)
21600 43200
t
60000
70000
80000
90000
100000
110000
TetherTension N
(e) Case 5: Tension (N) vs. time (s)
Figure 4.17: Plots of Tension vs. time for the five cases.
106
21600 43200
t
2.5·107
5·107
7.5·107
1·108
1.25·108
1.5·108
Tether Stress HPa L
Tether Intact
(a) Case 1: Stress (N/m2) vs. time (s)
21600 43200
t
2.5·107
5·107
7.5·107
1·108
1.25·108
1.5·108
Tether Stress HPa L
Tether Intact
(b) Case 2: Stress (N/m2) vs. time (s)
21600 43200
t
1.2·109
1.4·109
1.6·109
1.8·109
Tether Stress HPa L
Tether Intact
(c) Case 3: Stress (N/m2) vs. time (s)
21600 43200
t
2.5·107
5·107
7.5·107
1·108
1.25·108
1.5·108
Tether Stress HPa L
Tether Intact
(d) Case 4: Stress (N/m2) vs. time (s)
21600 43200
t
1·109
1.2·109
1.4·109
1.6·109
1.8·109
Tether Stress HPa L
Tether Intact
(e) Case 5: Stress (N/m2) vs. time (s)
Figure 4.18: Plots of Stress vs. time for the five cases.
107
Chapter 5
Dynamics of tethers with length
deployment
Ashasbeenshown inChapter4, it ispossibletodescribethemotionoftheMMET
around complex orbital paths. Now the deployment characteristics of the MMET
will be explored.
5.1 Deployment and recovery with the MMET
Thetwopayload armsof theMMET maybothbegainfullyemployed to launch two
payloads,andthispresentsaset of uniquechallengesforthe launcher. Bothpayload
tethers must be able to demonstrate all the necessary properties of a single system
tether, with the additional demands of a duallaunch system overlaid. The tethers
must be strengthened for a potential asymmetric release, and the oscillations that
this imposes on the system and tethers.
Thegeneralised coordinateshavebeenchosen inordertoascertainthe impactthat
the length deployment has on the MMET without overcomplicating the dynamics
108
of the system. The global coordinates R
and .
have been included to model the
orbital dynamics. The primary rotational coordinate, , has been expanded to
include both payload arms in the form of 1 and 2. Similarly, each payload arm
has been assigned a length coordinate, L1 and L2, to enable independent, rigid
motion of the extendable tether arm.
To limitthenumberofgeneralised coordinates,theinclinationvariable, , and the
outofplane coordinate, , have not been included as generalised coordinates. A
six degree of freedom model will be adequate to examine the dynamics of length
deployment,howeveraninedegreeoffreedommodelhasthepotential to introduce
orbital resonance to the deployment. At this early stage of research, nine degrees
of freedom are too many to successfully isolate and examine the effects of length
deployment on the system.
5.2 Additionoflength asageneralised coordinate
The six variables of interest – {R,
,
1, 2,L1,L2}– fulfil all the requirements of
generalised coordinates. Ananalysiswillbeperformed to investigatetheeffectthat
deployment and recovery have on the tether system orbit.
A system of equations will be constructed to include a symmetrical dumbbell
tether with a facility mass at the centre of the system. The important mass points
included in the model are the two identical payload masses, Mpayload, the facility
mass, Mfacility, and the n
discretized tether masses, Mtether, which are evenly distributed
along each tether.
As before, the stator is not included in this model because this would significantly
increase the complexity of the system.
As previously, Lagrange’s equations will now be derived for the system.
109
5.3 Constructing the equations of motion
A system of Lagrange’s equations is constructed using the rotational method described
in Chapter 3, similarly employed for the inclination equations in Chapter 4.
The localpositionof themasseswith respecttothefacility massarederived. This
allows the centre of mass with respect to the facility mass to be calculated. The
orbital positions are then assembled from these equations.
Simplifying assumptions have been made to lower the computational cost of assembling
and solving these equations.
.
Theorbitallength used inthecalculationofpotential energy isassumed tobe
the orbital radius vector, R. This is acceptable as L
« 1.
.
R
The orbital rotations have been decoupled from the local rotations. This is
sufficient to analyse the fast rotations of the dualpayload system, but should
be treated with caution when analysing the stationary characteristics of the
MMET.
.
The tethers are assumed to be rigid and inelastic.
.
Thecrosssectional area(CSA) of thetether isassumedtobeconstantalong
thetether length. TheCSA willalmost certainlybeafunctionoflength when
thetether isbuilt(i.e.tapered),however,thisassumptionbuilds ina level of
conservatism into the analysis which can be refined at a later stage.
These assumptions aredesigned tomaximisethebenefit of thesolutiontotheequations,
while minimising the computational cost and minimising the risk of losing
that information while studying a multiparameter system.
As outlined in the above assumptions, the orbital and local rotations have been
decoupled. The positional equation of the first payload therefore is defined as:
Pfinal1 = R,Z · (P0.CoM)(Pfacility.CoM)+(R 1,Z · P1) (5.1)
110
Ascanbeseen inEquation5.1,therotationsarepurelybased onaseriesofZaxis
rotations.
Using the rotation matrices defined in Section 3.5, along with the following definitions1
for the payload position:
..
L1
..
..
..
P1 = 0 (5.2)
..
..
0
leads to the expanded version of the positional equation for payload 1: When
evaluated, this becomes:
.....
.
R
cos .
CoMX
L1 cos( 1)
.....
.
.....
.
.....
.
Pfinal1 = R
sin.
 CoMY
+ L1 sin( 1) (5.3)
.
.
.
.
.
.
.
.
.
.
.
.
0 CoMZ
0
Thepositionof thecentre of massfromthefacility, P0.CoM, is calculated with the
system mass points including the payload, facility and tether mass points.
Thetether isdiscretized with each tethercomprising n
equal masspoints, as shown
in Equation 5.4 and Figure 5.1. The aim in discretizing the tether is to enhance
the potential and rotational energy expressions. The n
discrete mass elements were
chosen sothecentre of mass of thetetherfallsonthegeometric centre of thetether,
necessitating n
+1 spaces between the elements.
..
n
n
(M1+Mtether 1)L1 cos( 1)+(M2+Mtether 2)L2 cos( 2)
22
..
M1+M2+Mfacility+nMtether 1+nMtether 2
n
n
(M1+2 Mtether 1)L1 sin( 1)+(M2+2 Mtether 2)L2 sin( 2)
.
M1+M2+Mfacility+n
Mtether 1+n
Mtether 2 .
.
.
0
.
.
P0.CoM
(5.4)
.
=
.
To create Lagrange’s equations, the position, velocity and energy of each point
1Where L1
is the length in absolute terms between the facility and the payload.
111
Figure 5.1: Discretised tether diagram with n= 5 mass points Figure 5.1: Discretised tether diagram with n= 5 mass points
must be determined.
The positions of all the mass points are given as:
.....
.
R
cos.
CoMX
L1 cos( 1)
.....
.
.....
.
.....
.
Ppayload 1 = R
sin .
 CoMY
+ L1 sin( 1) (5.5)
.....
.
.....
.
0 CoMZ
0
.....
.
R
cos.
CoMX
L2 cos( 2)
.....
.
.....
.
.....
.
Ppayload 2 = R
sin .
 CoMY
+ L2 sin( 2) (5.6)
.....
.
.....
.
0 CoMZ
0
.....
.
R
cos.
CoMX
L1 cos( 1)
.....
.
.....
.
j
.....
.
Ptether 1,j = R
sin .
 CoMY
+ L1 sin( 1) (5.7)
.....
.
n
+1
.....
.
0 CoMZ
0
.....
.
R
cos.
CoMX
L2 cos( 2)
.....
.
.....
.
j
.....
.
Ptether 2,j = R
sin .
 CoMY
+ L2 sin( 2) (5.8)
.....
.
n
+1
.....
.
0 CoMZ
0
.
..
.
R
cos.
CoMX
.
..
.
.
..
.
.
..
.
Pfacility = R
sin .
 CoMY
(5.9)
.
..
.
.
..
.
0 CoMZ
where each j
is the number of the n
tether mass points.
The velocity of each point is determined. The payload 1 is used as an illustrative
112
example: in local axes in Equation 5.10 and in inertial axes in Equation 5.11
..
L. 1 cos( 1)L1 sin( 1) . 1
..
..
..
Vpayload1,local = L. 1 sin( 1)+L1 cos( 1) . 1 (5.10)
..
..
0
Vpayload1,inertial =
...
.
R.
cos()R. sin() L1 . 1 sin( 1)+L. 1 cos( 1)
...
.
...
.
.
.
..
.
R
sin()+R. cos()+ L1 . 1 cos( 1)+L. 1 sin( 1)
...
.
...
.
00
1
+ ·
M1+M2+Mfacility+2nMtether
..
.
n
(5.11)
 M1 + Mtether L1 . 1 sin( 1)L. 1 cos( 1)
..
2
.
..
.
..
n
.
..
 M1 + 2 Mtether L1 .
1 cos( 1)+ L. 1 sin( 1)
.
..
.
0
.
..
n
M2 + 2 Mtether L2 . 2 sin( 2)L. 2 cos( 2)
.
..
.
..
.
n
..
+
.
M2 + 2 Mtether L2 . 2 cos( 2)L. 2 sin( 2)
..
.
..
0
The expressions for energy are then calculated from the velocity and position
information.
The translational kinetic energy, Ttrans, of the system is calculated for the MMET
113
tether with n
discrete masses per tether2:
M1 M2
Ttrans = Vpayload 1 · Vpayload 1 + Vpayload 2 · Vpayload 2
2
2
n
Mtether
+ Vtether 1j · Vtether 1j
2
j=1
n
Mtether
+ Vtether 2j · Vtether 2j
2
j=1
1
+(M1 + M2 + nMtether + nMtether + Mfacility)(VCoM · VCoM) (5.12)
2
The rotational kinetic energy, Trot, of the system is found in a similar way:
Trot =(Ipayload 1 + Ipayload 2)· (!local · !local)
n
n
X.
(5.13)
+
Itether 1j · (!localj · !localj)+ Itether 2j · (!localj · !localj)
j=1 j=1
where:
..
..
0
100
..
..
..
..
..
..
!local =0
I3 = 010
..
..
..
..
.
.
001
=L2
=L2
Ipayload 1 1Mpayload 1 · I3
Ipayload 2 2Mpayload 2 · I3
Itether 1j =L21j
Mtether 1j · I3
Itether 2j =L22j
Mtether 2j · I3
The evaluated sum of the component rotational kinetic energies for n
=10 tether
masses are:
22
L12 M1 . 1 L22 M2 . 2 35Mtether 22 2
Trot =++ L1 . 1 + L22 . 2 (5.14)
2 222
2Where Mtether
is the discrete mass point, this sums to the total tether mass M
= nMtether.
114
The sum of the potential energies, U, are:
n
µmj
U
=
.
(5.15)
j=1 Pj
where Pj
are the mass point locations given in Equations 5.5 to 5.9.
The potential energy term is simplified by using the assumption that the mass
points are sufficiently close to the centre of mass of the system, that is PR
j
« 1, that
the potential energy of each mass point of the system may be represented as if it
were at a distance R
fromthe majorgravitationalbody. The evaluatedform of one
tether arm is therefore:
U
=
n
.
µmj
R
=
n
Mtether
R
(5.16)
j=1
Summing over all the mass points3 gives the full form for the potential energy of
the MMET:
µM1 µM2 20µMtether µMfacility
U
=   (5.17)
RRR
R
The Lagrangian is then given by L
= T
U.
Lagrange’s equations are then generated for the five generalised coordinates,
{R,
,
1,
2,
L1,
L2}as follows:
.
.
d
@T
@T
@U
dt
@q.j
 @qj
+
@qj
= Qj
(5.18)
5.4 Nonconservative forces
The motor on the MMET spins up the payload arms, while the countertorque
rotates the stator in the opposite direction. This is a nonconservative force in
action, as the energy is transferred from the solar cells4 to the rotational motion of
3There are 10 discrete mass points on each tether.
4Or the equivalent driving power source.
115
the payloads, via the motor. The nonconservative forces are derived according to
the principles of virtual work outlined in Section 3.4:
..
0
..
..
..
0
..
..
..
.
t
cos(
)
.
..
Q
= (5.19)
..
.
t
cos(
)
.
..
..
2 .
..
.
 Ldamp
L. 1
.
..
.
2 .
 Ldamp
L. 2
where t
is the torque magnitude, .
is the angle that the torque acts relative to
the inplane spin and Ldamp
is a braking term modelling the action of the spoolout
mechanism. The form of the tether friction braking was taken after the velocity
based barberpole braking system analysed in [Lennert and Cartmell, 2006].
5.4.1 Tether braking
It is proposed to use a barber pole braking system for this tether analysis. This
is a passive system employed to minimise the need for an active control system
or complicated reelout mechanism. The kinetic energy of the tether and payload
spoolsoutthetetherwithfrictionforcesproviding resistancetoensurethetether is
deployed in a controlled manner.
The braking force on the tether will be proportional to the square of the speed of
deployment:
L.
2
F
=  Ldamp
(5.20)
where Ldamp
is aproportional constant modellingthefrictionproperties andphysical
characteristics of the braking system. This allows a flexibility to either define the
damping coefficient to reflect a physical braking system or to specify a damping
coefficient and match a braking system to the numerical modelling.
116
5.5 Analysis of length deployment of tether
The MMET is different to other tether systems by virtue of the dualpayload motorized
arrangement of its tethers. This confers several advantages through the
symmetric arrangement of the payloads, however the dynamics of the spinup are
untested in an orbital environment5 .
A goal is set to study the dynamics of the MMET system: it must be capable of
imparting aV
=1000m/s
tothepayloads. This isamoderatevelocity increment
that will demonstrate the usefulness of the MMET system, allowing inclusion on a
small launcher or as a secondary payload.
The system of equations is solved in Mathematica.
, using NDSolve, for a time
of t
= 208129s
(approx 33 orbits). Figures 5.2 and 5.3 show the results for the
following initial conditions:
100
M1 = M2 =10kg,
Mtether = 10 =10kg,
Mfacility =500kg,
R.
[0] =0,.
[0] =0.00113905rad/s,
. 1[0] =0.1rad/s,
. 2[0] =0.1rad/s,
L. 1[0] =0,L. 2[0] =0,R[0] =7.378* 106 m,
[0] =0, 1[0] =0.01., 2[0] =180.1. ,
L1[0] =100m,
L2[0] =100m,
Ldamp
=5* 108 kg/m,
max =1000Nm
A further explanation of the starting conditions is below:
.
the starting radius R[0]isequivalentto1000km
abovethe surface oftheEarth.
.
the motor torque is proportional to time, and increases from t
=0 at t
=0s
to t
=1000Nm
at t
=208129s
˜ 33 orbits.
.
the initial conditions for the tether positions, 1 and 2 are chosen to be
slightly asymmetrical.
.
the tether mass, Mtether, is the discretized mass of a section of the tether, the
total tether mass in this case is given by nMtether =100kg.
5Fordetails ofground testing and analysis, see[Cartmell et al.,2003].
117
An asymmetryin the initial tether angles was chosen to counter numerical integrationproblems
–thenumerical solverhasgreatdifficultiesinintegratingtheequations
of motion when the tether arms are exactly 180. apart. The small asymmetry is
acceptable, as the tether arms will never be exactly opposite.
The solution to the equations of motion for the above initial conditions are given
in Figures 5.2 to 5.4.
5.5.1 Discussion of results for deployment
As can be seen in Figures 5.2a, the radius vector R
is constant – this is due to the
assumption intheexpressionforpotential energy thatthemassesactatdistance R
from the inertial origin.
ThemotortorqueshowninFigures5.2bisincreased linearlytoamaximumtorque
to avoid any large startup transients. This is key to ensure a smooth deployment
phase. A slow increase in the torque results in a gradual increase in the tether
rotational velocity, and a gradual increase in the stress in the tether.
The length of the payload 1 tether smoothly reels out from the hub, as shown
in Figures 5.2c, first slowly, then speeds up to a moderate and predictable linear
increase due to the tether braking6 .
Ifthetorquewassuddenly increasedfromzerotomaximum,the initialfewminutes
ofdeploymentwouldbeuncontrolled,and thisis likelytocausesignificantproblems
in several areas – including tether asymmetry and nonlinear flexing in the tethers.
This would invalidate the previous assumptions, perhaps leading to payload/stator
interaction, which is to be avoided at all costs.
The linear scaling in torque undoubtedly delays the overall time taken to deploy
the tether, compared to a flat maximum torque profile, and could be optimised
6Recall that the tether braking force is proportional to the square of the tether velocity:
L.
2
F
.Ldamp
118
further. This isaconservativeapproach toensurethesafedeploymentof thetether,
and inevitably leads to a longer deployment time than is absolutely necessary.
The torque drives the tether rotational speed, as shown in Figures 5.2d: it first
increases, then trends back to an asymptote at around 1.35rad/s.
The profile of angular velocity varying with length is shown in Figure 5.3b. This
shows that the rate of increase of . 1 is within acceptable limits and will not place
onerous demands on the strength of the tether.
The potential V
the payload may achieve is defined as the dependence between
therotationandlength,V
= .
1L1,and isshown inFigure5.3a. TheV
increases
approximately linearly,incontrasttothenonlinearbehaviouroftheangularrotation
and length generalised coordinates. At the end of the numerical integration, V
˜
1000m/s.
One issue with the MMET system is that there is a possibility of mutual interference
with the tethers of the payload a stator arms. As Figures 5.3c and 5.3d
show, in this case, the difference between the angle of the tethers and the lengths
of the tethers are small. The angular difference is nominally zero after a startup
transient, with an initial maximum of ˜ 0.5., which is acceptable for the spinup
criteria7 . The difference in lengths of the tethers are minimal – this is likely a result
of rounding errors encountered using a numerical solver.
Analysing the tension and stress in the tether shows that the tether is intact and
withinmaterialbounds. Thetethertension,shown inFigure5.4a,ataround100kN
is high, but within acceptable limits. The stress, shown in Figure 5.4b, assumes a
CSA of 64mm2, equivalent to a single circular section8 of diameter 9mm.
The stress level in the tether is below failure point. The tensile strength of Zylon
7This may be unacceptable for the final ‘launch phase’ of the payloads where the required
accuracy of the tether angle is likely to be more stringent.
8This is a simplifying assumption – a tether launch satellite would be likely to employ a combination
of a tapered tether to reduce the mass along the span of the tether and a multistrand
tether to mitigate against singlestrand failure.
119
fibre is5.9GP
a;with a safetyfactor of1.5, the maximum acceptable tensile strength
is3.9GP
a.Themaximumcalculated stress inthetetheris1.6GP
a, whichprovides
an additional level of redundancy in the tether. The total safety factor is approximately
15.
.69 =3.7.
The stress in the tether has been designed to be linear, avoiding stepincreases in
thestressand allowing the load todistributeevenlyalong the fibre inaprogressive
way.
The tether tension is calculated as the sum of the tension due to the tether mass
and the payload:
1 22 2
T
=
AL1 . 1 + L1 M1 . 1 (5.21)
2
where A
isthe crosssectional area(CSA) of thetether, and .
is the tether density.
The tether stress is defined conventionally as:
T
s
= (5.22)
A
5.6 Refining the equations of motion
A refinement of the tether mathematical model is proposed. Whereas the previous
model describes a discretized tether of constant mass, Mtether, this can be more
accurately modelled as a mass proportional to the tether length:
AL
Mtether = (5.23)
n
where .
is the tether density, and Ln
is the effective distance between mass points.
As shown previously in Figure 5.1, the tether mass model allows the dynamics of
a varying tether mass to be investigated. This is more accurate, better reflecting
the realities of the mass distribution when the tether is deployed.
120
The previous set of equations modelled the mass points as fixed fractions of the
tether which moved proportional to each other and did not vary in mass. This is
unsatisfactory because the entire tether mass would be distributed over the length
of the tether, whether the tether length is a few metres or at full deployment.
Obviously, this is not representative of the tether mass therefore a new model is
proposed to better describe the tether mass in Equation 5.23. This refines the
tether mass model to increase the tether mass in proportion to the length of the
tether, giving a more representative model of the tether mass.
Aswith theprevious investigation intothedeployment of thetether,theequations
of motionaresolvedforthefollowinginitial conditionsforatimeperiod9 of69376s:
M1 = M2 =10kg,
Mtether = AL/n,
Mfacility =500kg,
.
=1570kg/m3,A
=(0.008)2 =0.000064m2,n
=10
R.
[0] =0,.
[0] =0.00113905rad/s,
.
1[0] =0.1rad/s,
.
2[0] =0.1rad/s,
L. 1[0] =0,L. 2[0] =0,R[0] =7.378* 106 m,
[0] =0, 1[0] =0.01., 2[0] =180.1. ,
L1[0] =100m,
L2[0] =100m,
max =1000Nm,
Ldamp =5* 108 kg/m
Toavoidtheproblemofinitiallyoveraccelerating thepayloads,anonlineartorque
is applied by the tether motor, as shown in Figure 5.5b. This effectively limits
the initial spinup torque, limiting the angular acceleration and providing a stable
basis for controlled deployment. The maximum value of angular acceleration is
0.025rad/s2, and the maximum angular speed is 10.2rad/s. Figure 5.5d shows a
timehistory of .
with a moderate angular acceleration in the first few seconds of
spinup,followedby aprogressive and controlleddeployment ofthetether thereafter.
The torque is controlled in such a way to avoid an initial spike of acceleration, to
progress rapidly to a hightorque plateaux and furthermore, to provide a smooth
9This corresponds to exactly 11 orbits.
121
transition between the two.
7500
t
= max
· exp +50 (5.24)
t
+1
where the maximum torque provided by the motor, max
˜ 1000Nm
and t
is the
elapsed time in seconds.
A significant advantage of the motor torque control is that the tether is deployed
linearly,with theintentionoflimiting thejerk oraccelerationtothetether,asshown
in Figure 5.5c. The linear deployment is very successful in limiting the tension and
stress in the tether, see Figure 5.7a and Figure 5.7b respectively.
The expression for the tether tension at the facility (i.e. the point of greatest
tension) has changed compared to the previous tether mass model described in
Equation 5.21. The tension is not calculated with respect to the centre of mass of
the tether in the centre of the tether, it is instead assumed that the tether mass acts
in its entirety at the root of the tether. The tether tension in the discretized tether
is calculated as:
T
= AL12 .
12
+ L1 M1 .
12
(5.25)
where A
is the crosssectional area of the tether, and .
is the tether density.
When examining stresses in the tether, the maximum stress levels for the exponential
torque case appears higher at 2.1GP
a, compared to 1.3GP
a
for the linear
torquemodel.This isduetothediffering tension(and therefore,stress) expressions
in the two cases, and can not be directly compared.
However, if the stress is qualitatively examined, it is interesting to compare the
stress profiles. The time at which the maximum stress is applied is very different;
with the linear torque case, the stress is maximum at the end of the spinupjust
before payload release; with the exponential torque case, the stress quickly rises to
maximumearly inthespinup,and iseased slightlybeforepayload releaseachieving
122
a stress before release of 1.75GP
a. The maximum allowable stress, for comparison,
is 3.9GP
a
The time taken to spin the payload to release velocity is significantly less in the
exponential torque case – the time taken to spinup is one third of the linear torque
case. This is partly due to the difference in mass modelling, and partly due to the
torque profile used. This represents a marked improvement, and will undoubtedly
enhance the economic strengths of the MMET as a renewable launch platform.
In terms of tether synchronisation, Figures 5.6c reveals that the tethers are well
separated by approximately 180., and furthermore, they tend to settle at a point
diametrically opposed to each other. This is an important result, especially when
the tethers are initially separated by a few fractions of degrees as would occur in a
real spinup scenario. The solution to the equations of motion were based on the
initial conditions 1 =0.01. and 2 =180.1. .
Likewise, inFigure5.6d,thedifferenceintetherlengthsareminimal.Thismaynot
be true in the real world, as material properties and coefficients of friction are very
difficult to exactly match between the two tether arms. Clearly, the asymmetrical
spinup scenario requires further study.
As it is clear that deploying the tethers with a MMET system is achievable for a
modest payload, the recovery of the tether model shall be examined.
5.7 Recovery of tethers
OneofthestrengthsoftheMMETis in itsroleasareusable launcherofpayloads. In
order to reuse the tethers, they must be despun and then recovered to the facility.
Thehardwarerequired toperformthedespinisalargebraking system –thedisc
brake system on a car axle is similar in operation to what may be required. There
are important differences, however, as the disc brake must be cooled. On Earth’s
123
roads,thecoolingisprovidedby conductiveheattransferof air; onorbit,conduction
and radiation are necessary to channel the excess heat away from the brake to the
vacuum of space.
The required braking energy is significant – approximately equal to the energy
supplied toacceleratethetethertoreleasevelocity. Thepowersuppliedby thesolar
cells to the motor will be comparable to the power required to brake the assembly,
assuming equal time spent accelerating and decelerating the payloads. Given that
thesolarcellsarecurrently inthe kW
range, thebraking system andheatdissipation
systems should be sized accordingly and should not burdensome to source.
5.7.1 Recovery profile
Recovery requires a different approach to deployment – the two are conceptually
opposite, but can not be treated as mirror images when analysing the technicalities
as there are important differences between the two cases.
Deployment simply requires a spinning motor and a braking system, the tether
arms held rigid by the rotation.
When recovering the arms, there are requirements that must be fulfilled:
.
the tethers must be slowed gradually to avoid backlash and tether tangling
.
there must be a minimum spin rate to keep the tether rigid
.
the tethers must be fully recovered to enable a subsequent deployment
The initial slowdown phase is similar to the spinup phase – energy is dissipated
during the slowdown to minimise tension in the tether, making recovery easier. The
motorrequiredtoreelthetether inforstoragemay require lesspower(doing less
work against the tether tension) and therefore will be lighter and cheaper.
124
5.7.2 The equations of motion for recovery
The equations of motion are the same as thedeployment equations of motionfor the
orbital generalised coordinates, R
and , and the spin generalised coordinates, 1
and 2.The lengthgeneralised coordinates, L1 and L2 shown inEquations5.26 and
5.27,differ inthefactthatthey arenotbraked tostopfastdeployment,butthey are
retracted undermotorpower. Similarly,the initial conditionsand nonconservative
forcing terms10 have been altered to reflect the differences between decelerating and
accelerating the tether arms.
The aim of the new length equations of motion is to provide a consistent length
.
recovery rate for the tethers: L
= v. This models a constant speed motor reeling
in the tethers. A small acceleration term, L¨
, isadded tomakeboth length equations
into two ODEs, such that the equations of motion may be solved numerically.
¨
0.000001 L1 + L. 1 = Lrec
(5.26)
¨
0.000001 L2 + L. 2 = Lrec
(5.27)
The nonconservative forces in the right hand side of Largange’s equations are:
..
0
..
..
..
0
..
..
..
.
 damp
. 1
.
..
Q
=
.
(5.28)
..
. .
damp
.
2
..
..
..
.
(Lrec)
.
..
(Lrec)
Thenonconservativeforcesforthe length equationsof motionprovideaconstant
reelin speed, modelled by the rate Lrec. Similarly, the spin generalised coordinates
The right hand side terms of Lagrange’s equations, shown in Equation 5.28.
125
dictate a linear braking proportional to the rotational speed of the tethers.
5.7.3 Solving the equations of motion for recovery
The recovery stage necessitates a twostage approach to recovering the tether. The
first stagerecoversthetethersclose intothefacility,at whichpoint aslowerrecovery
stage is used to safely guide the remaining length of tether into the facility.
The tether arms must be kept rigid11, and this requires a minimum rotational
velocity of thetethersystem. Thetether isdecelerated toasmallrotational velocity
in the first phase with a constant braking force. Once the tether arm is rotating
slowly, close to the facility, a much lower braking force and marginally lower tether
retractionspeed areemployed toensurethesaferecovery ofthe last sectionoftether.
The initial conditions for the first recovery phase, numerically solved for a time
period of 50000s, are:
M1 = M2 =10kg,
Mtether = AL/n,
Mfacility =500kg,
.
=1570kg/m3,A
=(0.008)2 =0.000064m2,n
=10,
R.
[0] =0,.
[0] =0.00113905rad/s,
.
1[0] =1.0rad/s,
.
2[0] =1.0rad/s,
L. 1[0] =0,L. 2[0] =0,R[0] =7.378* 106 m,
[0] =0, 1[0] =0.01., 2[0] =180.1. ,
L1[0] =1000m,
L2[0] =1000m,
Torquemax
=0Nm,
Lrec =0.015kg/m,
damp =1500Ns
Again, an asymmetry in the tether angles was chosen to reflect the unknowns that
affect the payload release dynamics. Additionally, the numerical solver has great
difficulties in integrating the equations of motion when the tether arms are exactly
180. apart.
The second phase of recovery has the initial conditions similar to the first phase,
and is numerically solved for a time period of 25000s. Initial conditions that differ
11This is an explicit assumption made when assembling the equations of motion.
126
are listed below:
. 1[0] =0.026707rad/s,
. 2[0] =0.026707rad/s,
L. 1[0] =0,L. 2[0] =0,
L1[0] =250m,
L2[0] =250m,
Lrec =0.010kg/m,
damp =100Ns
The goal for recovery of the tethers are to provide a stable and controlled reelin,
avoiding scenariosthat would lead toabreak of thetethersuch astetherasymmetry
or uncontrolled resonance in the tether.
Foracircularorbitand a lineartetherrecoveryprofile,thepayload tethersremain
separated by approximately 180. . The maximum range of tether movement from
thesymmetriccentre lineisapproximately0.43.,asshown inFigures5.9cand5.12c.
Furthermore, there are no significant deviations or other areas of concern in either
the spin angle or the tension or stress profiles.
Thetensionandstress inthetetheraregradually lowered,whilekeeping thetether
taut, as seen in Figures 5.10a, 5.13a, 5.10b and 5.13b. This will not guarantee a
successful recovery, but will reduce the risk of a transient phenomenon that may
cause the tether to momentarily rise above the ultimate strength of the material,
thus breaking the tether.
The rotational speed of the arms initially increase due to conservation of momentum,
asthe length(andtherefore inertia) ofthetethersarereduced. Angular
momentum isconserved,and asFigure5.9b clearlyshows,thetetherpayloadsrotate
at a faster rate while reeled in. This does not lead to an increase in the tension or
V
due to the careful choice of the rotationalbraking system(andbraking constant
damp)and the reelin speed(Lrec).
The second phase of tether recovery is shown in Figures 5.11 to 5.13. This occurs
with 250m
of tether still outside the facility, i.e. one quarter remaining to recover.
The tether arms are rotating slowly enough to reduce the tension in the tether to
127
minimal levels, while keeping the tether rigid enough to ensure the equations of
motion remain valid in describing the motion.
The tether arms are decelerated at a lower rate than the previous recovery phase,
with damp =100Ns. The tether recovery rate is slowed slightly to L.
=0.1m/s.
Similartothe first recoveryphase,theremaining tether isslowlyretracted intothe
facility with no adverse dynamic effects. The angle between the tethers remains at
approximately 180. ±0.5. . However, at such a low tension, it is unknown whether
the tether can remain rigid without coiling or exhibiting other nonlinear behaviour.
An indepth analysisshouldbeperformed toascertainthetether’sbehaviour insuch
an environment, as a study that involves nonrigid tether behaviour is beyond the
scope of this thesis.
Thetimetakentodeploythenrecoverthetetheris69376 +(50000 +25000) =
144376s, ˜ 40 hours, or a little under 22 orbits12, not taking into account the time
taken to release the payloads and the time to stabilise the tethers between the two
recovery stages. This short timespan allows a small swarm of 18 nanosatellites to
be launched from the MMET within a month.
As one of the key strengths of the MMET is the reusability of the platform, this
is an encouraging finding.
5.8
Deployment and recovery on an elliptical orbit
The MMET describing an elliptical orbit may interact with the tether while deploying
or retracting the tethers. In order to study this, a sensitivity analysis is carried
out to ascertain if this hypothesis is correct.
The deployment and recovery of the tether are carried out under identical initial
12
One orbit takes approx. 6307s
at an altitude of 1000 km.
128
conditions to the equations describing the discretized mass model with the exponential
tether tension model.
The nondimensional ellipse eccentricity parameter, e
=0.25, is supplied to the
initial conditions of Lagrange’s equations of motion through the orbital angle initial
condition: . =0.00111382rad/s.
This is calculated using the semilatus rectum and orbital parameters, as in [Szebehely
and Mark, 1998], to calculate the orbital angular velocity thus:
The semilatus rectum(SLR) is calculated:
SLR
= R(1+e
cos()) (5.29)
where e
is the eccentricity.
This isusedtocalculatetheorbital angularvelocity as in[Ziegler,2003]:
v 2
µ
(1+e
cos)
.= 3 (5.30)
SLR
While deploying the tether, the orbital parameter .
does not adversely affect the
dynamics of the tether deployment. As shown in Figure 5.14b, the key deployment
metrics – the smooth deployment of the entire tether length and a progressive increase
in stress over time – are realised. When the parameters for the circular and
elliptical are compared, the parameters share much in common. When comparing
the tether length with angular velocity, for example, Figures 5.14b and 5.6b are
almost identical.
Similarly, when the recovery of the tether undergoing a elliptical orbit is analysed,
the solutions to the equations of motion present a similar result. The length and
angular velocity are very similarly matched over the two cases, as is shown when
comparing Figures 5.15b and 5.9b.
129
2
5.9 Centre of mass movement
Additional stresses may be exerted on the tether by centre of mass movement. This
movementcausestheroot of thetetherto libratealong theorbitalpath,subjecting
the tether to additional loading.
The movement of the centre of mass is shown in Figure 5.16 for the deployment
cases and Figure 5.17 for the recovery case. The movement of the centre of mass
with respect to the facility mass is given by Equation 5.4.
The maximum movement of the CoM, for the deployment cases, are 0.07m
and
0.0017m
– the latter is the discretized mass model using exponential torque. The
effort expended to model the tether mass in more detail has stripped back a layer
of conservatism, and shown that the centre of mass movement is minimal. This
reinforcestheearlier findingsof therotationandlengthgeneralised coordinates; the
symmetrical angleof thetethersand theequal lengthsof thetethersboth contribute
to a very symmetrical system.
Themaximummovement oftheCoMfortherecovery case isanorderof magnitude
greater, at 0.7m
for the first phase and 0.08m
for the second phase. This is due to
the starting assumption of asymmetrical tether angles: 1[0] =0.01. and 2[0] =
180.1. .
Whereas the difference in the tether angle in Figure 5.6c at full deployment tends
to zero, the initial conditions to recover the tether are asymmetrical. This is an
allowanceforthepayload release,whichisanunknownquantityinMMETdynamics.
Therefore, an initial asymmetry of 0.1. and 0.01. was chosen to investigate the
recovery dynamics. To remedy this unknown, a separation study is recommended,
to fully investigate the effects that a symmetrical (or asymmetrical) dual release
have on the tethers, however this is outside the scope of this document.
Onfurther examination,thegraphs of the magnitude ofposition and acceleration
were plotted against time, and are show in Figures 5.18a and 5.18b.
130
The maximum acceleration, due to the moving CoM, encountered in the recovery
phase is approximately 0.8m/s2, therefore the force exerted by the moving centre
of mass,givenatotal systemmassof10+10+100+100+500 =720kg
is:
F
= ma
=0.8* 720 =576N
(5.31)
Giventhis isseveral ordersof magnitude lessthanthetethertensionof110,
000N
inFigure5.10a, itdoesnotpresent asignificant risk tothetether. Indeed, itexplains
the slight high frequency tension superimposed on the overall smooth downward
trend shown in that figure.
Onpayload release,the likelihood of asmalltetherasymmetry may besignificant.
Removing thisasymmetry, insofaraspossible,couldbeachievedby repeated cycling
of short recovery anddeploymentphases, oranactivecontrol systemonthepayload
release mechanism.
5.10 Conclusions
Deploying and recovering the tethers on the MMET system is not a trivial matter;
however, they can be achieved in a well controlled and structured way.
A method for deployment and recovery in the orbital plane have been outlined,
givingconfidencethattheMMETiscapableand suitableforuseasareusable launch
platform.
131
5.11 Figures
R’
10
1.5·10
10
1·10
11
5·10
11
5·10
10
1·10
10
1.5·10
t
86400 172800
(c) Tether 1 lengthL1
(m)vs. time t
(s) (d) Tether 1rotational velocity . 1
(rad/s)vs.
time t
(s)
Figure 5.2: Initial MMET Spinup
Motor Torque
1000
800
600
R
66667 77
2·104·106·108·101·101.2·101.4·10
400
200
t
86400 172800
(a) Phase plot: radius velocity .R
(m/s) vs. (b) Applied motor torque t
(Nm) vs. time
radius of CoM R
(m) t
(s)
L1 D@.1,t.
t
86400 172800
0.5
1.5
2
132
DeltaV L1
86400 172800
200
400
600
800
1000
1200
800
600
400
t
.1'
0.8
1.2 1.4 1.6 1.8 2
(a) Payload V
(m/s)vs. time t
(s) (b) Angular velocity . 1
(rad/s) vs. length
L1
(m)
.1.2
L1L2
200
50000
t
0.0075
0.005
5·10
0.0025
0.00001
t
100000 150000 200000
0.000015
0.0025
0.00002
0.005
0.0075
0.000025
50000 100000 150000 200000
6
(c) Difference in tether angles 1
 2
p
(rad) (d) Difference in tether lengths L1
L2
(m)
vs. time t
(s) vs. time t
(s)
Figure 5.3: Initial MMET Spinup
Tether 1 Stress HGPAL
TetherTension N
86400 172800
Tether 1 Intact
100000
9
1.5·10
80000
9
1.25·10
9
60000
1·10
8
7.5·10
40000
8
5·10
20000
8
2.5·10
t
t
86400 172800
(a) Tether tensionT
(N)vs. time t
(s) (b) Tether stresss
(N/m2)vs. time t
(s)
Figure 5.4: Initial MMET Spinup
133
R’
11
3·10
11
2·10
11
1·10
11
1·10
11
2·10
11
3·10
11
4·10
Torque
R
66667 77
2·104·106·108·101·101.2·101.4·10
800
600
400
200
t
10000 20000 30000 40000 50000 60000 70000
(a) Phase plot: radius velocity R.
(m/s) vs. (b) Applied motor torque t
(Nm) vs. time
radius of CoM R
(m) t
(s)
L1 D@.1,t.
1000
800
600
400
200
4
2
t
t
21600 43200 64800
21600 43200 64800
10
8
6
(c) Tether 1 lengthL1
(m)vs. time t
(s) (d) Tether 1rotational velocity . 1
(rad/s)vs.
time t
(s)
Figure 5.5: Initial MMET Spinup, discretized tether mass model
134
DeltaV
L1
1000
800
800
600
600
400
400
200
200
.1'
t
2 4 6 810
21600 43200 64800
(a) Payload V
(m/s)vs. time t
(s) (b) Angular velocity . 1
(rad/s) vs. length
L1
(m)
.1.2
0.008
0.006
0.004
0.002
0.002
0.004
L1L2
t
10000 20000 30000 40000 50000 60000 70000
0.00002
0.00004
t
10000 20000 30000 40000 50000 60000 70000 0.00006
0.00008
(c) Difference in tether angles 1
 2
p
(rad) (d) Difference in tether lengths L1
L2
(m)
vs. time t
(s) vs. time t
(s)
Figure 5.6: Initial MMET Spinup, discretized tether mass model
135
Tether 1 Stress HGPAL
TetherTension N
Tether 1 Intact
9
120000
2·10
100000
9
1.5·10
80000
60000
9
1·10
40000
8
5·10
20000
t
t
21600 43200 64800
21600 43200 64800
(a) Tether tensionT
(N)vs. time t
(s) (b) Tether stresss
(N/m2)vs. time t
(s)
Figure 5.7: Initial MMET Spinup, discretized tether mass model
D@.1,t.
L1
1000
1
900
800
0.8
700
0.6
600
0.4
500
400
0.2
t
43200
t
21600 43200
21600
(a) Tether 1 lengthL1
(m)vs. time t
(s) (b) Tether 1rotational velocity . 1
(rad/s)vs.
time t
(s)
Figure 5.8: MMET reel in, first phase
136
DeltaV
L1
21600 43200
200
400
600
800
1000
1000
900
800
700
.1'
1.08 1.09
1.11 1.12 1.13 1.14
t
(a) Payload V
(m/s)vs. time t
(s) (b) Angular velocity . 1
(rad/s) vs. length
L1
(m)
L1L2
.1.2
10000 20000 30000 40000 50000
0.0075
0.005
0.0025
0.0025
0.005
0.0075
1
0.5
t
t
10000
20000
30000
40000
50000
0.5
1
(c) Difference in tether angles 1
 2
p
(rad) (d) Difference in tether lengths L1
L2
(m)
vs. time t
(s) vs. time t
(s)
Figure 5.9: MMET reel in, first phase
137
21600 43200
t
20000
40000
60000
80000
100000
TetherTension N
(a) Tether tension T (N) vs. time t (s)
21600 43200
t
2.5·108
5·108
7.5·108
1·109
1.25·109
1.5·109
1.75·109
Tether 1 Stress HGPAL
Tether 1 Intact
(b) Tether stress (N/m2) vs. time t (s)
Figure 5.10: MMET reel in, first phase
3600 7200 10800 14400 18000 21600
t
50
100
150
200
250
L1
(a) Tether 1 length L1 (m) vs. time t (s)
3600 7200 10800 14400 18000 21600
t
0.005
0.01
0.015
0.02
0.025
0.03
D@.1,tD
(b) Tether 1 rotational velocity .1 (rad) vs.
time t (s)
Figure 5.11: MMET reel in, second phase
138
DeltaV
L1
t
3600 7200 10800 14400 18000 21600
(a) Payload V
(m/s)vs. time t
(s) (b) Angular velocity . 1
(rad/s) vs. length
L1
(m)
L1L2
.1.2
1
0.0075
0.005
0.5
0.0025
t
t
25000
5000
10000
15000
20000
25000
0.0025
0.5
0.005
0.0075
1
(c) Difference in tether angles 1
 2
p
(rad) (d) Difference in tether lengths L1
L2
(m)
vs. time t
(s) vs. time t
(s)
Figure 5.12: MMET reel in, second phase
0.005 0.01 0.015 0.02 0.025 0.03
.1'
50
100
150
200
250
5
4
3
2
1
5000 10000 15000 20000
139
TetherTension N
6
5
4
3
2
1
t
3600 7200 10800 14400 18000 21600
(a) Tether tensionT
(N)vs. time t
(s)
Tether 1 Stress HGPA.
Tether 1 Intact
100000
80000
60000
40000
20000
t
3600 7200 10800 14400 18000 21600
(b) Tether stresss
(N/m2)vs. time t
(s)
Figure 5.13: MMET reel in, second phase
R
7
1.2·10
7
1.1·10
7
1·10
6
9·10
6.28319 12.5664 18.8496 25.1327 31.4159 37.6991 43.9823
(a)
L1
1000
800
600
400
200
T
.1'
2 4 6 810
Orbital radius R
(m) vs. orbital angle (b) Angular velocity . 1
(rad/s) vs. length
.
(rad) L1
(m)
Figure 5.14: MMET deployment with elliptical orbit
140
6.283185301721.75965683674076161894.3854991575259952321955.31837257944132032818.741185394256950385897932385
T
9·106
1·107
1.1·107
1.2·107
R
(a) Orbital radius R(m) vs. orbital angle
(rad)
1.08 1.09 1.11 1.12 1.13 1.14
.1'
700
800
900
1000
L1
(b) Angular velocity .1 (rad/s) vs. length
L1 (m)
Figure 5.15: MMET recovery with elliptical orbit
0.1 0.05 0.05 0.1
CoM_x
0.1
0.05
0.05
0.1
CoM_y
(a) Centre of mass position, linear mass model
0.0015 0.001 0.0005 0.0005 0.001
CoM_x
0.0015
0.001
0.0005
0.0005
0.001
0.0015
CoM_y
(b) Centre of mass position, discretized mass
model
Figure 5.16: MMET deployment, centre of mass movements
141
1 0.75 0.5 0.25 0.25 0.5 0.75 1
CoM_x
1
0.75
0.5
0.25
0.25
0.5
0.75
1
CoM_y
(a) Centre of mass position, discretized mass
model, first stage
0.1 0.05 0.05 0.1
CoM_x
0.1
0.05
0.05
0.1
CoM_y
(b) Centre of mass position, discretized mass
model, second phase
Figure 5.17: MMET recovery, centre of mass movements
10000 20000 30000 40000 50000
t
0.1
0.2
0.3
0.4
0.5
0.6
0.7
AbsHCoML
(a) Absolute position of CoM, discretized
mass model, first phase
10000 20000 30000 40000 50000
t
0.2
0.4
0.6
0.8
AbsHCoM_aL
(b) Absolute acceleration of CoM, discretized
mass model, first phase
Figure 5.18: MMET recovery, position and acceleration of centre of mass
142
143
5.11.1 Linear kinetic energy
The linear kinetic energy of the MMET with length deployment and ten discrete tether masses, expanded from the analytical version
in Equation 5.12, is as follows:
2
Ttrans = Mtether(11M1 (M2 + Mfacility)+145M1 Mtether +5Mtether (7(M2 + Mfacility)+85Mtether))L. 1 +(11M2 (M1 + Mfacility)
2)L. 22 2 ..2 ..R.
2
+5(7M1 +29M2 +7Mfacility)Mtether +425Mtether +11M1 R2 +22M1 M2 R2 +11M2 R2 +22M1 Mfacility R2 +22M2 Mfacility
2 ....2 .
+11Mfacility R2 +440M1 Mtether R2 +440M2 Mtether R2 +440Mfacility Mtether R2 +4400Mtether R2 +11M12 R2 .2 +22M1 M2 R2 .2
+11M22 R2 .2 +22M1 Mfacility R2 .2 +22M2 Mfacility R2 .
2 +11Mfacility
2 R2 .2 +440M1 Mtether R2 .
2 +440M2 Mtether R2 .2
+440Mfacility Mtether R2 .
2 +4400Mtether
2 R2 .2 +22(M1 +5Mtether)(M2 +5Mtether)L1 sin( 1  2)L. 2 .
1
2222 22
2222 22
+11M1 M2 L1 . 1 +11M1 Mfacility L1 . 1 +145M1 Mtether L1 . 1 +35M2 Mtether L1 . 1 +35Mfacility Mtether L1 . 1 +425Mtether
2 L1 . 1
22(M1 +5Mtether)(M2 +5Mtether)cos( 1  2)L1 L2 . 1 . 2 +(11M2 (M1 + Mfacility)+5(7M1 +29M2 +7Mfacility)Mtether
2
2
+425Mtether
2)L2 . 2 22(M1+5Mtether)(M2+5Mtether)L. 1 (cos( 1 2)L. 2+L2 sin( 1 2) . 2)22(M1+M2+Mfacility+20Mtether)
(5.32)
Chapter 6
Dynamics of spacewebs
This chapter was funded as part of the ESA Ariadna study into spacewebs; the
final report waspresented as[McKenzie et al.,2006].The authorwishestoacknowledge
the hard work and contributions to this report from M.P Cartmell and M.
Vasile of the University of Glasgow and D. Izzo of the ESA ACT.
6.1 Introduction
As was outlined in Section 2.3, a generic satellite that may be constructed or re
configured asneeded isauseful concept. Studieshavepreviouslyshownthat robots
may be deployed in this way to reconfigure the net structure [Kaya et al., 2004a],
[Nakano et al., 2005].
This chapter will consider the fundamental conceptual design of an appropriate
and generic thin membrane, its orbital mechanics and control.
A model of a net in orbit is presented here, with robots moving along the surface
of the net, simulating reconfiguration of the system. Useful guidelines to the initial
conditions of the net and components arefound andgive a startingpointforfurther
144
studies into the field.
This new class of structures with robots moving over the surface will be defined
as ‘Spacewebs’, with the robots correspondingly named ‘Spiderbots’. A concept
illustration of the idea is shown in Figure 6.1.
Figure 6.1: Concept rendering of the spaceweb with spiderbots on the surface.
Spacewebs have one definite advantage over single use spacecraft: the ability to
reconfiguretheweb tosuitthemissionorenvironment. Duediligencemusttherefore
be observed to ensure that the act of reconfiguring the web does not endanger the
web or satellite. Investigating the movement of the robots while moving over the
web is essential: the movement of their mass and momentum may significantly affect
145
the web dynamics.
The stability of the system is investigated while the robots crawl along the web
surface in two predefined motions:
.
Three robots simultaneously moving along the outer catenaries of the web
.
Three robots simultaneously moving from the centre to the outside of the web
6.2 Modelling methodology
This study will aim to develop an analytical model of the fully deployed spaceweb
system that may be numerically integrated.
When constructing a numerical simulation of an existent system, it is important to
consider that the simulation can only be said to represent a real life scenario when
the model has been validated and the initial conditions used in the calculation are
representative of the simulation scenario.
This may seem like an obvious statement, but these facts are often overlooked for
convenience or speed.
In modelling the spaceweb, it is important to make the modelling process transparent
and verifiable. In this spirit, the steps followed to assemble the spaceweb
model are outlined paying particular attention to the steps followed to produce a
model in Mathematica and the process of validating the model against previous
models and the real life case where possible.
146
6.3 Web meshing
6.3.1 Web structure
Togain an accurate estimate ofthe energy of the web, the massdistribution over the
web is discretized [Ziegler, 2003]. An algorithm has been designed in Mathematica
to take the triangle bounded by the ith
subspan, the (i
+ 1)th
subspan and the
facility mass and divide this into equally sized smaller triangles.
The spaceweb is modelled as s
subspans, clustered around the central facility
mass. Each subspan is considered to be rigid and massless, held rigid by the
centripetal force of rotation around the centre of mass. An idealized point mass mi
isplaced attheend of thesubspan,at length li, with theposition of each endpoint
mass given by Pilocal in Equations 6.5, 6.6. In this form, this is very similar to
thetetherspreviously modelled,thedifferences aremultiple subspans and massless
tethers.
The web is stretched between the subspans i
and i
+1 and replicated to give s
triangular webs. Each triangular web section is idealised as an elastic plane containing
the mass of that section, however the internal elastic forces within the web
are not modelled at this stage.
6.3.2 Dividing the web
Each web section is divided into n
equal sections(referred tohereas ‘divisions’), in
thedirectiongivenbythe linejoining themidpointofthesubspanendsandthe
facility, as represented in Figures 6.2 to 6.4.
For each of the n
web sections;
147
triangles = (n
+ 1)2
triangles in top row = 2n
+ 1
row configuration = {2n
+ 1,
2(n
1)+ 1,
2(n
2)+ 1,
.
.
.
,
3,
1}
number of rows = n
+ 1
nodes = 1
2(n
+ 2)(n
+ 3)
midpoints = (n
+ 1)2
Figure 6.2: 0 divisions Figure 6.3: 1 division Figure 6.4: 2 divisions
6.3.3 Web divisions
A point mass is placed at each midpoint, therefore the total number of masses for
the webs are s(n
+1)2, where s
is the number of subspans. That is to say, the
number of masses in the web increases as the square of the number of divisions.
Thiscan lead toavery largenumberof masspointstoconsiderforafinegrainweb.
Assembling the spaceweb withthree subspans(n
= 3) and 0 divisions gives the
layout shown inFigure6.5,withFigure6.6showing thesame layout with1division,
and finally Figure 6.7 shows the same layout with 2 divisions. Each red dot is a
masspoint on the web.
Figure 6.5: 3 subspans
– 0 divs
Figure 6.6: 3 subspans
– 1 div
Figure 6.7: 3 subspans
– 2 divs
148
Figure 6.9: 5 subspans
Figure 6.8: 4 subspans
Figure 6.10: 6 sub
– 2 divs
– 1 div
spans – 3 divs
The spaceweb layout may be expanded to analyse different web configurations. A
square layout with1divisionpersection isshown inFigure6.8; apentagonallayout
with 2 divisions per section is shown in Figure 6.9; and a hexagonal layout with 3
divisions per section is shown in Figure 6.10 with total of 103 masses1 .
6.4 Equations of motion
6.4.1 Centre of mass modelling
The mass of the web is likely to be unevenly distributed. Modelling this requires
anexpressionforthecentreof mass(CoM) position, inthiscaseusing thefacility
mass as the origin.
For n
masses with positions {Xi,Yi,Zi}, the position of the centre of mass about
the facility, in terms of the local {X,
Y,
Z}coordinate system centred on the facility
mass, is:
.
.
n
n
n
.
.
.
.
.
MiXi
MiYi
MiZi
.
.
.
.
i=1 i=1 i=1 Pfacility.CoM
n
,
n
,
n
(6.1)
.
.
.
.
.
.
Mi
Mi
Mi
.
.
.
i=1 i=1 i=1
This isused whentaking theposition from the centre of mass to the facility mass,
as PCoM.facility
= Pfacility.CoM
1103 masses = 96 web midpoints + 6 subspan masses + 1 central facility mass.
149
P_1
CoM
P_s1
P_s
Facility
P_2
P_3
Figure 6.11: Simplified spaceweb layout with s
subspans
6.4.2 Rotations
Rotational matrices are used to rotate the starting vectors to their positions on the
local and inertial planes. The local position vector in Equation 6.5 is the matrix
product of a series of rotations. A diagram of the spaceweb configuration is shown
in Figure 6.13 with the web removed for clarity.
..
10 0
..
..
..
R i,X = 0 cos( i) sin( i) (6.2)
..
..
0 sin( i) cos( i)
..
cos(i) 0 sin(i)
..
..
..
Ri,Y = 0 10 (6.3)
..
..
sin(i) 0 cos(i)
150
(i,j+1)
(i+1,j)
(i,j)
L
2
L
2
h
2h
3
3
(i,j)
L
2
L
2
h
2h
3
3
Figure 6.12: Midpoint location
..
cos(i) sin(i)0
.
.
.
.
Ri,Z = .
.
sin(i) cos(i) 0 .
.
(6.4)
.
.
0 0 1
PiLocal = R i,X · Ri,Y · Ri,Z · Li
(6.5)
where Li
= {0,
0,li}T
, aligned along the local Zaxis
Piinertial = R,Z (P0.CoM Pfacility.CoM +(R i,X · Ri,Y · Ri,Z · Li)) (6.6)
Rotations may be performed using the rotation order Z
then Y
then X:
The subspan is rotated around the facility by an angle .
with the axis of rotation
intheXaxisonly. Thisvector isthenrotated around thefacilityby anangle a
with
the axis of rotation in the Yaxis only. Finally, this vector is then rotated around
the facility by an angle f
with the axis of rotation in the Zaxis only.
151
y2
y1
q
y3
localZ
Ylocal
Earth
inertialY
inertialX
....
....
Mass3
Mass2
CoM
....
..
Facility
Radius of CoM
Orbital Path of
SpaceWeb CoM
Mass1
....
....
Mass3
Mass2
CoM
....
..
Facility
Radius of CoM
Orbital Path of
SpaceWeb CoM
Mass1
Figure 6.13: Spaceweb diagram with 3 subspans shown in inertial space. Note:
the local axes are in the YZ plane.
The angles are projected about their respective axes i.e., the Zrotation, , is
contained in the X
 Y
plane. The angles themselves can be compared to the
standard aerospace rotations, with the subspan direction from the facility to the
edge taken as the tailnose direction. In this case: .
is the yaw angle, in the plane
of the spaceweb; a
is the pitch angle, out of the spaceweb plane; and f
is the roll
angle, the twist in the subspan.
Considering rotations about the .
direction only – that is constraining the motion
of the spaceweb to be inplane – gives the following equation
Piinertial = R,Z (P0.CoM Pfacility.CoM +(R i,X · Li)) (6.7)
Where P0.CoM = {R,
0,
0}, and representstheposition of the centre of massfrom
Earth.
152
6.4.3 Position equations
The positions of the masses in the local coordinate system must be translated and
rotated into the inertial (Earth centred) coordinate system. The local position
vectors are addedtotheposition vectorfromtheCoMtothefacility andtheposition
vectorfromtheEarthtotheCoM.Theresultant vector isthenrotatedintotheEarth
centred Inertial axis system.
6.5 Energy modelling
For every mass point in the system with position in the inertial frame, Pi, the
following steps are undertaken in order to find the Lagrangian energy expression L:
the respective velocities, Vi
are found:
.
Vi
= Pi
(6.8)
@t
thekinetic energies(linear Tlin and rotational Trot)are obtained and summed:
n
1
Tlin = miVi
· Vi
(6.9)
2
i=1
!2
n
1 . 1 + . 2 + . 3
· (6.10)
Trot = mi
Pilocal Pilocal
23
i=1
the potential energies, U, are obtained and summed:
n
µmi
U
= v (6.11)
Pi
· Pi
i=1
and the Lagrangian is found:
L
= Trot + Tlin U
(6.12)
153
Themomentofinertiaof each masspoint onthewebiscalculated usingtheparallel
axis theory about the central facility, then multiplied by the square of the average
angular velocities of the three subspans to give the rotational kinetic energy.
Theaverageangularvelocityofthethreesubspans,asshown inEquation6.13,was
chosen in preference to the actual calculated angular velocities of each mass point
tosimplify theLagrangianenergy expressionand to lightenthecomputationalload.
Instead of27individual angularvelocities,every masspoint ontheweb wasassumed
tohaveonecommonangularvelocity. There isanerror inthisassumption,butthis
will be acceptably small because the actual velocity is approximately equal to the
average velocity.
@.
1 @ 1 @ 2 @ 3
= + + (6.13)
@t
3 @t
@t
@t
ave
For each mass point, the Lagrangian energy expression is constructed by considering
the total kinetic and potential energies of the system.
d
@T
@T
@U
 += Qj
(6.14)
dt
@q.j
@qj
@qj
Lagrange’s equations are generated for all the generalised coordinates as specified
in Equation 6.14.
6.5.1 Virtual work terms
The elastic force the web imparts on the subspan ends are modelled as a simple
elastic tether [Miyazaki and Iwai, 2004], obeying Hooke’s Law: F
= Kx:
2p
F
= KH
(Pi
Pi+1) (6.15)
s
where H[.. .] is the Heaviside function.
The forces are then used to calculate the right hand side of Lagrange’s equations,
154
Qi, through consideration of the virtual work.
@i
Qi
= Ksci
; (6.16)
qi
where i
isthedifferencebetweenthesubspansdefinedby the inplane( ) coordinates
as:
1 =0 (6.17)
2 =0 (6.18)
2p
2p
3 = 2  1  H 2  1 
ss
2p
2p
 1  3  +2p
H 1  3  +2p
(6.19)
ss
2p
2p
4 = 3  2  H 3  2 
ss
2p
2p
 2  1  H 2  1  (6.20)
ss
2p
2p
5 = 1  3  +2p
H 1  3  +2p
ss
2p
2p
 3  2  H 3  2  (6.21)
ss
Note: iftheanglebetweensubspans1 and3( 1  3) is examined, this would
be usually be large and negative – e.g. for an equally spaced spaceweb. this angle
would be 240. instead of 120. . Therefore to counteract this, an offset of 2p
has
been added to Equation 6.19 to ensure that the angle is positive and similar in size
tothe other subspan angles( 2  1 and 3  2).
155
6.5.2 Simplifying Assumptions
TosolvethefullsetofLagrange’sequationsforthisspaceweb system inareasonable
time on a desktop PC requires some simplifying assumptions to be made.
Firstly, five generalised coordinates are chosen defining the space web rotating in
theplane normal to the radius vector: (R,
,
1, 2, 3). Thediagram inFigure6.13
shows the layout of the spaceweb, where R
is the orbital radius, .
is the true
anomaly, n
is the anglebetween the nth
subspan and the Xlocal
coordinate system.
Iftheoutofplanemotionofthespaceweb weretobe included,thiswould increase
thenumber ofgeneralised coordinatestoeight, aseach subspan musthavetheoutof
plane motion defined.
If thespaceweb orbitsexclusively inplane,asimplifying assumptionmay beperformed
in terms of the potential energy expression:
n
µmi
U
=
R
(6.22)
i=1
6.5.3 Robot Position Modelling
Robots may move across the spaceweb in order to reconfigure the web or other tasks
as discussed previously. To model this as an kinetic energy based term while retaining
the flexibility of defining the path without hardcoding this into the equations,
the robot position vector needs to be kept in a very general form.
Therobotpositionvector inEquation6.23 isdefined similarlytoEquation6.6: a
vector addition of the R vector, the CoM vector and the vector defining the local
position of the robot where
{P
RobotX
,
P
RobotY
,
P
RobotZ
}local
arethe local vectorpositionsof therobot inthe
156
X,
Y,
Z
axes with the origin coincident on the facility mass.
P
Robotinertial
= R,Z(P0.CoM
Pfacility.CoM
+
(6.23)
{P
RobotX
,
P
RobotY
,
P
RobotZ}local)
Therobot localpositionvectormaybeafunctionof time,positionof thesubspan
masses ( 1, 2, 3), an arbitrary path function or any other smooth generalised
function.
The kinetic energies of the robot (both linear and rotational) are derived from
the positions of the robot as before. Both the Lagrangian energy expression and
Lagrange’sequationscontaintermsforthe localposition,velocity,and acceleration
of the robots.
Keeping the robot position functions in the most general form allows for rapid
reconfiguration of the robot paths and may lead to, for example, robot control
studies in the future.
6.6 Spaceweb dynamics
6.6.1 Dynamical simulation
The spaceweb dynamics are heavily governed by the centreofmass movement of
the system. The spaceweb system is very rarely symmetrical(both in reality and
in simulations) and this asymmetry can lead to unstable and even chaotic motion
in certain configurations of the spaceweb.
Several different parameters were considered to be possible influences on the stability
of the spaceweb system, including:
.
variation of the masses of the web
.
the masses of the robots, the subspan and facility masses
157
.
the position and velocity of the robots
.
the general angular momentum of the robots
.
the angular velocity of the web
.
and the starting configuration of the web or ‘web symmetry’
Lagrange’s equations for the spaceweb system are solved using the Mathematica
numerical integrator NDSolve[. . .]. The package LiveGraphics3D [Kraus, 2007] is
then used to display any animated graphics in a webbrowser.
Additionally,thearrowsusedto indicatedirection intheCoMplotsweregenerated
with the CurvesGraphics6 package [Gorni, 2008] for Mathematica 5.1.
6.7 Investigating the stability of the spaceweb
To test the stability of the spaceweb with robot movement over the web, several
different numerical experiments have been considered, each investigating the impact
the following have on the centre of mass movement of the spaceweb:
.
web mass
.
robot mass
.
web asymmetry
.
robot crawl velocity
6.7.1 Mass of the web
The mass of the web was found to have a negative impact on the stability of the
system. The higher the mass of the web, the more likely the triangle is to deform
from the perfect equilateral shape, and the more likely the system is to exhibit
158
chaotic motion. This is likely to have been due to the increase in the ratio of the
web mass to the other masses, causing the web mass to dominate the other mass
terms. Figures B.1, B.2, B.3, B.4, B.5, B.6, B.7 and B.8 demonstrate this effect,
which were simulated with Case 12 in which only the mass of the web is varied.
The initial conditions used to numerically integrate the equations are given in
Table C.1.
The mass of the web reaches a limit of around 40 kg, after which, the web is no
longer able to keep its shape in this configuration.
Runs E01–E04 show the spaceweb can be stabilised while the robots move across
the surface, the mass of the web increasing to a critical value of approximately 25%
of the total system mass. The maximum CoM displacement is 0.85m
for these four
runs. In runs E05–E06, the mass of the web is increased to 44.82kg
and 48.25kg
respectively. This causes significant instability in the web for the 48.25kg
web,
showing that for certain configurations of the spaceweb system, the mass of the
web is a critical parameter. The marginal increase in web mass gives a noticeable
degradation in the behaviour of the spaceweb past that critical point.
6.7.2 Robot mass
The momentum of the robot was thought to be a significant parameter in the stability
of the spaceweb, therefore the robot mass was investigated as a possible
parameter. Runs 1821 tested values of Mrobot3 from 1kg
to 100kg
and the CoM
positions are shown in Figures B.9, B.10, B.11, B.12, B.13 and B.14.
The initial conditions used to numerically integrate the equations are given in
Table C.4.
Counterintuitively, increasing the mass of the robot caused the spaceweb to be
2With symmetrical robot paths.
3From 3% to 17% of the system mass.
159
comemorestable!This isduetothecentripetal effect of moremassontheoutsideof
the web causing the subspans to become more rigid and effectively act as a damper
to the system.
There is a sharp divide between the stable behaviour of the robot masses set
at 16kg
shown in Figure B.12, and the unstable behaviour with the masses at
15kg
shown in Figure B.11. This critical dividing line between the acceptable and
unacceptable system movement is clear as it is unexpected, and raises significant
issues for the construction of the spaceweb, especially with a lighter robot.
6.7.3 Web asymmetry
Thesinglebiggestdriverof instability onthespaceweb systemwasfound tobethe
asymmetry in the web configuration and mass distribution.
Very small asymmetries in the initial conditions of the spaceweb (compared to
a perfectly symmetrical equilateral triangle) are amplified in certain configurations
of the spaceweb, and may create large instabilities in the system, especially for
energetic configurations such as high web mass and/or large angular velocities. A
smalldifference inorientation(ontheorderof1.)of the subspans maylead to large
asymmetriesbetweenthesubspans inarelativelyshorttime(˜ 100s).
The cases chosen to investigate the effect of asymmetry on the web are as follows:
.
Case 1: Robots 1,2,3 walk along the edges of the web from subspan to next
subspan – (symmetrical)
.
Case2:Robots2,3walkalongtheedgesof theweb whileRobot1 isstationary
onsubspan1 –(unsym.)
.
Case 3: Robot 3 walks along the edges of the web while Robots 1,2 are stationary
onsubspans1,2 – (unsym.)
.
Case4:Robots1,2,3walkalong fromthecentretotheirsubspans – (sym.)
160
.
Case 5: Robot 1 spirals outwards from the centre to the edges, Robots 2,3 are
fixed on subspans – (unsym.)
.
Case 6: Robot 1 walks along from the centre to subspan1, Robots 2,3 remain
in centre –(unsym.)
The initial conditions used to numerically integrate the equations are given in
Table C.2.
Cases 1 and 4 are symmetrical, shown in Figures B.15 and B.18. Case 1 has a
stable CoM position throughout the simulation, with a maximum CoM travel of
only 0.46m. Case 4 is symmetrical, however has a maximum CoM travel of 10.96m
and is on the limits of what is defined as acceptable movement for the CoM. In
contrast, Cases 2, 3, 5 and 6 are asymmetrical, shown in Figures B.16, B.17, B.19,
and B.20. The asymmetric cases show a large CoM movement of between 18m
and
38m, and are unstable.
Compounding theproblemof unevenmass,the initial conditionsof thethreesubspans
were found to influence the stability of the system. The equations could not
be solved with .
= {0. ,
120. ,
240.}or spacing the web subspans by exactly 120. ,
for reasons not known at this time. Initial conditions for .
were implemented as a
triplet; the three subspan values of {0+ ,
120. ,
240.  }offer a solution to this
problem.
Moregenerally, it isexpected thatformostconfigurationsof thetriangularspaceweb,
holding the three subspans permanently rigid at exactly 120. separation will
be impossible as small perturbations to the web in the space environment will be
constantly experienced.
Adding the effect of many asymmetrical robot masses exacerbates the problem of
uneven mass distribution. To remedy this, the spaceweb must be configured to
occupy as low an energy state as the mission will allow: low web mass; light, slow
moving robots; low angular velocity of the web. In this state, the asymmetry of the
161
web still exists, but it is kept to a manageable level.
6.7.4 Robot crawl velocity
Thefastertherobotstravelled along theweb,thegreaterthe likelihood of unstable
behaviouroftheweb,astheCoMdisplacementplotsshow inFiguresB.21 andB.25.
Both are simulated with Case 1 – symmetrical robot paths.
The initial conditions used to numerically integrate the equations are given in
Table C.3.
When robot crawl velocity was increased from 0.1m/s
to 10m/s, by lowering the
time fixedfrom1000s
to10s,themaximumCoMdisplacement increased marginally
–
from 0.3m
to 1.5m.
These two cases of mass dependant stability and velocity dependant stability are
clearly energy related. The higher the energy contained in the system, the more
likely the system is to exhibit unstable behaviour. This boundary has not yet been
clearly defined. Generally speaking, however, the slower the robot movement across
the web, the more likely the web is to remain stable given a perfectly symmetrical
web. Given the levels of CoM displacement, the robot velocity is not a significant
destabilising factor in this scenario.
6.8 Statistical investigation into stability
A series of solutions to the spaceweb equations of motion were examined to examine
the effect of five parameters on the spaceweb stability. The mass of the web,
Mweb, the mass of the three daughter satellites, Msat, the robot mass, Mrobot,
the subspan angular coordinate, .
and the average subspan angular velocity @.
@t
were thought to have an effect on the maximum movement of the Centre of Mass
(CoM)of the spaceweb system. Therefore, the solutions to 36 different sets of ini
162
tial conditions were found, with 25 =32 runs complemented by 4 centre point runs.
The 36 runs were performed in Mathematica, and the results of the maximum CoM
displacement were entered into a statistical package, Design Expert 7.03. Design
Expert thenperforms an analysis of variance(ANOVA) calculation onthefactorial
data provided. ANOVA is a collection of statistical models and their associated
procedures which compare means by splitting the overall observed variance into differentparts.
A5degreeoffreedommodel(orless, ifrequested) isthenextrapolated
from the data and analysed for statistical significance.
The five most strongly correlated model parameters were:
Mweb Msat+ Mweb
* @.
+ Msat
* Mrobot
* @.
 Mrobot
@t
@t
+ and @.
 were statistically relevant in the model, but to a lesser degree than
@t
the other parameters. As they have a stronger effect when combined with other
variables,the(un)stabilising influencetheyhavemaybeovertakenby theseother
factors.
Parameters areshown with a(+)indicating astabilising effect ora()indicating
a destabilising effect, with the full table of values located in Appendix C.
As Mweb has the strongest influence over the model, it would be most beneficial
to the stability to reduce the mass of the web.
Thereverse istrueofthedaughtersatellitemass,Msat+, increasingthisparameter
will tend to stabilise the system.
For combinations of main parameters such as Mweb
* @.
+, it would be most
@t
beneficial to increase the product to stabilise the system. In this case, the angular
momentum of the web will tend to rigidize the system and enhance the stability.
6.8.1 Masses of the web and satellite
The influence the masses of the web and the daughter satellites have on the CoM
positionareshown inthecontourgraph,Figure6.14. Thecontoursareshadedfrom
163
blue(stable) throughgreen andyellow to red(unstable) andfollow the maximum
CoM displacement from the central facility.
DesignExpert® Software CoM Dist
10.00
CoM Dist
Design Points
55.7692
0.00324405 8.00
X1 =A:Mweb
10
15
20
25
30
35
40
45
4
X2 =B:Msat
Actual Factors
C: Mrobot = 5.50
D: psi = 0.06
E: psidot = 0.55
Msat B:
6.00
4.00
2.00
27.00 87.75 148.50 209.25 270.00
A: Mweb
Figure 6.14: Contours of CoM movement while varying web and satellite mass
For the small sample space of initial conditions, the most stable point would be a
high satellite mass(>
8kg)and low web mass(<
50kg). A CoM displacement of
above 10m
is unstable, and would be very difficult for a robot to manoeuvre over
the surface.
Ingeneralterms, this clearly shows the stabilisingeffect ofthe mass ofthe satellites
and the destabilising effect of the mass of the web.
6.8.2 Masses of the web and robot
The influence the masses of the web and the robot have on the CoM position are
shown inthecontourgraph,Figure6.15. Thecontoursareshadedfromblue(stable)
164
throughgreen andyellowto red(unstable) andfollowthe maximumCoMdisplace
ment from the central facility.
DesignExpert® Software CoM Dist
10.00
CoM Dist
55.7692
0.00324405
7.75
X1 =A:Mweb
X2 = C: Mrobot
27.00 87.75 148.50 209.25 270.00
10
15 20
25 30
35
40 45
50
Actual Factors
B: Msat = 6.00
D: psi = 0.05
Mrobot C:
5.50
E: psidot = 0.10
3.25
1.00
A: Mweb
Figure 6.15: Contours of CoM movement while varying web and robot mass
The relative strengths of the robot and web masses are shown, with the mass of
the web exerting a much stronger influence over the behaviour of the web. The
stability boundary isshown inthebottomlefthand cornerof thegraph. Astability
margin of 5m
is preferable, here the web mass must be kept below 50 kg
and the
massof therobot,although less influential, wouldbebest suited atbelow5kg. This
shows that if a low @.
is required, the web can be made more stable with careful
@t
consideration of variables.
6.8.3 Mass of the web and angular velocity
The influence the masses of the web and the angular velocity, @.
, have on the CoM
@t
positionareshown inthecontourgraph,Figure6.16. Thecontoursareshadedfrom
165
blue(stable) throughgreen andyellow to red(unstable) andfollow the maximum
CoM displacement from the central facility.
DesignExpert® Software CoM Dist
1.00
CoM Dist
55.7692
0.00324405
0.78
X1 =A:Mweb
5
10
15
20
25
30
35
40
45
50
X2 = E: psidot
Actual Factors
B: Msat = 10.00
C: Mrobot = 1.00
D: psi = 0.06
psidot E:
0.55
0.33
0.10
27.00 87.75 148.50 209.25 270.00
A: Mweb
Figure 6.16: Contours of CoM movement while varying web mass and @.
@t
In a highstability scenario, the options to configure the spaceweb are plentiful.
Carefully choosing the previous parameters that lead to a high probability of stability,
namely ahigh satellite mass(10kg)and a low robot mass(1kg), affords a
larger range of possible values for the mass of the web and the angular velocity of
the web.
If the largest acceptable CoM movement is limited to 5m, then there are three
choices available,depending on theprimary requirement of the spaceweb system. If
a low web mass is required, then a low @.
must match, and vice versa. Alternatively,
@t
@.
if a high mass is required, a high @t
must be specified, and vice versa. The final
option is for a low web mass and a high @.
, affording a large safety margin that
@t
may be advantageous in, for example, a pilot study mission.
166
6.8.4 Mass of the robot and angular velocity
The influence the mass of the robot and the angular velocity, @.
, have on the CoM
@t
positionareshown inthecontourgraph,Figure6.17. Thecontoursareshadedfrom
blue(stable) throughgreen andyellow to red(unstable) andfollow the maximum
CoM displacement from the central facility.
DesignExpert® Software CoM Dist
1.00
CoM Dist
55.7692
0.00324405
0.78
X1 = C: Mrobot
X2 = E: psidot
1.00 3.25 5.50 7.75 10.00
5
10
15
20
25
30
35
40
45
50
Actual Factors
A: Mweb = 270.00
B: Msat = 10.00
D: psi = 0.05
psidot E:
0.55
0.33
0.10
C: Mrobot
Figure 6.17: Contours of CoM movement while varying robot mass and @.
@t
The final comparison is for a high mass specification – a 270kg
web mass and
10kg
satellite mass. This combines the two strongest influences on the system, the
former destabilising and the later stabilising. For a stable system, it is essential to
ensure the robot mass is small and the angular rotation rate is large.
167
6.9 Conclusions and recommendations
6.9.1 Conclusions
The achievement of a stable, reconfigurable web orbitingin space with robots moving
along its surface is a realistic goal. There must, however, be limits to the
behaviour of the robots and the configuration of the web. Both the robots and the
web must be as light as possible and the robots must be as slow moving as the
mission allows, given that this limits the destabilising effects of these parameters.
The daughter satellites must be as heavy as possible, and the angular rotation rate
must be as large as possible to maximise these stabilising effects. In all cases, the
web configuration must be as symmetrical as practicable.
6.9.2 Recommendations
To capitalise on the knowledge gained in this research, there are a few areas that
may be expanded on for potential future exploration.
Constructing the web while on orbit would mitigate the deployment phase of the
mission, eliminating one of the major failure modes of the system. Using the robots
to weave the web could be inspired by spiders on Earth, for example, just as the
space web has been inspired by the Furoshiki cloth.
Using the knowledge gained through moving masses over the web, there are many
applications that would benefit from simulating webbased constructions with moving
masses. Solarsails and antennae could be assembled or reconfigured by robots
moving across their surface.
168
Chapter 7
Conclusions
7.1 Tether rotations
.
Rotations using a different system of equations are equally able to describe
the same system.
.
Therefore,thechoiceof rotational systemcanbetailored insuch away tosuit
the analysis.
7.2 Tethers with inclination
.
A demonstration mission launching a microsatellite from MEO to Lunar orbit
is achievable using current technology. The safety margins, however, are
extremely small and the MMET launcher would be located in an orbit with
large amounts of debris.
.
The inclination term does not significantly alter the dynamics of the MMET
system. However, any outofplane component in the local axes – whether on
169
an inclined orbit or not – interacts with the orbital parameters to create an
unsteady oscillation in the tether. This is unacceptable for payload aiming
and release, therefore the local outofplane angle should be minimised insofar
as possible.
7.3 Deployment and recovery of tethers
.
Deploying and recovering the tethers on the MMET system is not a trivial
matter; however, they can be achieved in a well controlled and structured
way.
.
Amethodfordeployment and recoveryintheorbitalplanehavebeenoutlined,
givingconfidencethattheMMET iscapableand suitableforuseasareusable
launch platform.
7.4 Spacewebs
.
The achievement of a stable, reconfigurable web orbitingin space with robots
moving along its surface is a realistic goal.
.
There are limits to the behaviour and configuration of the robots and web:
the robots must be light and as slow moving, the web must be as light, the
daughter satellites must be heavy and the angular rotation rate must be as
large as possible to maximise these stabilising effects. In all cases, the web
configuration must be as symmetrical as practicable.
170
7.5 Further study
It is my belief that tethers, and in particular the MMET, have a significant role to
play in mankind’s future in space. By addressing the following issues uncovered in
this thesis, the case for using the MMET will be strengthened.
The stator arm will require to be modelled and analysed if the concept of the
MMET is ever to be successful. The possibility of interference between the tether
and stator arms will logically follow on from this analysis.
Likewise, the concept of a spaceweb is both novel and interesting, and doubtless
will be studied further. It is essential to perform further ground tests or even zero
gravity studies of theMMET and spacewebstoprovide validation of the equations
of motion.
171
Bibliography
J.A. Andi´on and C. Pascual. Useful experiences in a series of deployable booms for
CLUSTER satellites. In R. A. Harris, editor, ESA SP480: 9th European Space
Mechanisms and Tribology Symposium, pages 113–120, September 2001. URL
http://adsabs.harvard.edu/abs/2001smt..conf..113A.
Yuri Artsutanov. V kosmos na electrovoze (to the cosmos by electric train).
Newspaper report, Komsomolskaya Pravda, July 31 1960. URL http://www.
robotstore.com/download/Artsutanov_Pravda_SE.pdf. Translated by Joan
Barth Urban, Roger G. Gilbertson and M. Molloy.
E.A. Belbruno. Lunar capture orbits, a method of constructing Earth Moon trajectories
andthe Lunar GAS mission. DGLR, and JSASS, International Electric
Propulsion Conference, 19th, Colorado Springs, CO, USA, AIAA1987(1054):1–
10, 1987.
Vladimir V. Beletsky and Evgenii M. Levin. Dynamics of Space Tether Systems
(Advances in the Astronautical Sciences). American Astronautical Society, 1993.
URL http://eleanor.lib.gla.ac.uk/record=b2210829.
Matthew Bennett, Michael F. Schatz, Heidi Rockwood, and Kurt Wiesenfeld. Huygens’s
clocks. Proceedings of the Royal Society of London. Series A: Mathematical,
Physical and Engineering Sciences, 458(2019):563–579, 2002. doi:
10.1098/rspa.2001.0888.
172
R. Biesbroek and G. Janin. ESA bulletin 103 “Ways to the Moon?”, August 2000.
URL http://www.esa.int/esapub/bulletin/bullet103/biesbroek103.pdf.
A. S. Brown. Spreading spectrum of reinforcing fibers,. Aerospace America, 27(1):
16, January 1989.
Joseph A. Carroll. Tether applications in space transportation. Acta Astronautica,
13(4):165–174, April 1986. doi: 10.1016/00945765(86)900615.
M.P. Cartmell and D.J. McKenzie. A review of space tether research. Progress in
Aerospace Sciences, 44(1):1–21, January 2008. doi: 10.1016/j.paerosci.2007.08.
002.
M.P. Cartmell, S.W. Ziegler, R. Khanin, and D.I.M. Forehand. Multiple scales
analyses of thedynamics of weakly nonlinear mechanical systems. Transactions of
ASME, AppliedMechanicsReviews,56(5):455–492,2003. doi: 10.1115/1.1581884.
M.P. Cartmell, S.W. Ziegler, and D.S. Neill. On the performance prediction and
scale modelling of a Motorised Momentum Exchange Propulsion Tether. In M. S.
ElGenk, editor, AIP Conf. Proc. 654: Space Technology and Applications International
Forum STAIF 2003, pages 571–579, January 2003. STAIF2003571.
M.P.Cartmell,C.R.McInnes, andD.J.McKenzie. Proposalsfor a continuousEarth
Moon cargo exchange mission using the Motorised Momentum Exchange Tether
concept. Proc. XXXII International Summer School and Conference on Advanced
Problems in Mechanics, APM2004, pages 72–81, June 24 July 1 2004.
M.P.Cartmell,D.I.M.Forehand,M.C.D’Arrigo,D.J.McKenzie,Y.Wang, andA.V.
Metrikine. Approximate analytical solution for a motorised momentum exchange
tether on a circular earth orbit. In Proceedings of theXXXIVth SummerSchool on
Advanced Problems in Mechanics, pages 119–130. Russian Academy of Sciences,
Repino, St.Petersburg, Russia, 2006.
173
Arthur C. Clarke. The Fountains of Paradise. Victor Gollancz, 1979. ISBN 057502520
4.
M.L. Cosmo and E.C. Lorenzini. TethersinSpaceHandbook(3rd ed). Number 3.
NASA Marshall Space Flight Center, Cambridge, MA, December 1997. URL
http://www.tethers.com/papers/TethersInSpace.pdf.
A.N. Danilin, T.V. Grishanina, F.N. Shklyarchuk, and D.V. Buzlaev. Dynamics of
a space vehicle with elastic deploying tether. Computers & Structures, 72(13):
141 – 147, 1999. ISSN 00457949. doi: 10.1016/S00457949(99)000395.
Saswato R. Das. Final thoughts from Sir Arthur C. Clarke. Online, March 2008.
URL http://spectrum.ieee.org/mar08/6075.
A.B. DeCou. Tether static shape for rotating multimass, multitether, spacecraft for
‘triangle’Michelson interferometer. Journal ofGuidance,Control, andDynamics,
12(2):273–275, 1989. URL http://www.aiaa.org/content.cfm?pageid=406&
gID=20401.
M. Eiden and M.P. Cartmell. Overcoming the challenges: Tether systems roadmap
for space transportation applications. AIAA/ICAS International Air and Space
Symposium and Exposition, Ohio, USA, 1417 July 2003.
R.L.Forward.TethertransportfromLEOtotheLunar surface. InProceedings of the
27thAIAA/ASME/SAE/ASEE JointPropulsionConference andExhibit.AIAA.,
1991. URL http://www.tethers.com/papers/LEO2Lunar’92.pdf. 912322.
R.L. Forward and R.P. Hoyt. Space tethers. Scientific American, February 1999.
R.L. Forward and R.P. Hoyt. Failsafe multiline Hoytether lifetimes. 1st AIAA/SAE/
ASME/ASEE Joint Propulsion Conference, San Diego, CA, AIAA Paper
9528903, July 1995.
174
P.E. Glaser. Power from the sun, its future. Science, 162(3856):857–861,1968. doi:
10.1126/science.162.3856.857.
G. G´omez, A. Jorba, J. Masdemont, and C. Sim´o. Study of the transfer from the
Earth to a halo orbit around the equilibrium point L1. Celestial Mechanics and
Dynamical Astronomy, 56(4):541–562, August 1993. doi: 10.1007/BF00696185.
Gianluca Gorni. Mathematica plots with directional arrows the curvesgraphics6package,
9thAugust2008. URL http://users.dimi.uniud.it/~gianluca.
gorni/Mma/Mma.html.
M.A. Green, K. Emery, D.L. King, S. Igari, and W. Warta. Solar cell efficiency
tables (version 22). Progress in Photovoltaics: Research and Applications, 11:
347–352, 2003.
M.D. Griffin and J.R. French. Space Vehicles Design and Construction. AIAA
Education Series, 1991.
A.C. Hindmarsh. ODEPACK, a systematized collection of ODE solvers. In
R.S.Stepleman etal., editor, ScientificComputing(vol.1 ofIMACSTransactions
on Scientific Computation), volume 1, pages 55–64, NorthHolland, Amsterdam,
1983.
Walter Hohmann. Die Erreichbarkeit der Himmelsk¨orper Untersuchungen ¨uber
das Raumfarhtproblem. Oldenbourg, R., M¨unchen und Berlin, 1925. ISBN 348623106
5.
Walter Hohmann. The Attainability Of Heavenly Bodies. NASA(CASI), November
1960. URLhttp://ntrs.nasa.gov/details.jsp?R=191962. ReportNASATTF
44, Translated from original German.
R.P. Hoyt. Stabilization of electrodynamic space tethers. In M.S. ElGenk, editor,
AIP Conf. Proc. 608: Space Technology and Applications International Forum
175
STAIF 2002, pages 570–578, January 2002. URL http://www.tethers.com/
papers/ED_Stabilization.pdf.
John D. Isaacs, Allyn C. Vine, Hugh Bradner, and George E. Bachus. Satellite
elongation into a true “skyhook”. Science, 151(3711):151, Feb 1966. doi: 10.
1126/science.151.3711.682.
Solar System Dynamics Group JPL. JPL’s HORIZONS system. Website, August
2008. URL http://ssd.jpl.nasa.gov/horizons.cgi.
Osman M. Kamel and Adel S. Soliman. Sensitivity of transfer orbits to small errors
in position and velocity at cutoff. Acta Astronautica, 58(5):5, March 2006. doi:
10.1016/j.actaastro.2005.09.006.
N. Kaya, M. Iwashita, S. Nakasuka, L. Summerer, and J. Mankins. Crawling robots
onlargeweb inrocket experimentonFuroshikideployment.In 55th International
AstronauticalCongress2004Vancouver,Canada,2004a.URL http://www.esa.
int/gsp/ACT/doc/POW/ACTRPRNRG2004IACSPSCrawling_robots.pdf.
N. Kaya, S. Nakasuka, L. Summerer, and J. Mankins. International rocket experiment
on a huge phased array antenna constructed by furoshiki satellite with
robots. In 24th International Symposium on Space Technology and Science,
Miyazaki, Japan, 2004b.
Jerry B. Keiper. Mathematica developer conference, Boston, MA, USA, June 1992.
URL http://library.wolfram.com/infocenter/Conferences/4687/.
Wang Sang Koon, Martin W. Lo, Jerrold E. Marsden, and Shane D. Ross. Dynamical
systems, the threebodyproblem and space mission design. InB.Fiedler,
K.Groger, andJ.Sprekels, editors,InternationalConference onDifferentialEquations,
Berlin, 1999, pages 1167–1181, Berlin, 1999. World Scientific.
W.S. Koon, M.W. Lo, J.E. Marsden, and S.D. Ross. Low energy transfer to the
176
Moon. Celestial Mechanics and Dynamical Astronomy, 81(12):63–73,2001. URL
www.gg.caltech.edu/~mwl/publications/papers/lowEnergy.pdf.
M. Kraus. LiveGraphics 3D webpage, 2007. URL http://wwwvis.informatik.
unistuttgart.de/~kraus/LiveGraphics3D/index.html.
G.A. Kyroudis and B.A. Conway. Advantages of tether release of satellites from
elliptic orbits. Journal of Guidance, Control, and Dynamics,11(5):441–448,1988.
SimoneLennertandMatthewP.Cartmell. Analysis anddesign of africtionbrakefor
momentum exchange propulsion tethers. Acta Astronautica, 59(811):923 – 930,
2006. ISSN00945765.doi: 10.1016/j.actaastro.2005.07.042.SelectedProceedings
of the Fifth IAA International Conference on Low Cost Planetary Missions.
Martin W. Lo and Shane D. Ross. The Lunar L1 gateway: Portal to the stars and
beyond. In AIAA Space 2001 Conference, Albuquerque, New Mexico, August
2830 2001.
E.C. Lorenzini. Errortolerant technique for catching a spacecraft with a spinning
tether. Journal of Vibration and Control, 10(10):1473 – 1491, 2004. ISSN 10775463.
doi: 10.1177/1077546304042062.
E.C. Lorenzini, M.L. Cosmo, M. Kaiser, M.E. Bingham, D.J. Vonderwell, and
L. Johnson. Mission analysis of spinning systems for transfers from low orbits
to geostationary. Journal of Spacecraft and Rockets, 37(2):165–172, 2000. URL
http://www.aiaa.org/content.cfm?pageid=406&gID=3562.
C.R. McInnes and M. Cartmell. Dynamics of propellantless propulsion systems.
Elsevier, 2006. URL http://eprints.cdlr.strath.ac.uk/6507/. ISBN 012373562
9.
D. McKenzie, M. Cartmell, G. Radice, and M. Vasile. Space webs final report
054109a. Ariadna study 054109a, University of Glasgow and ESA
177
ACT, 2006. URL http://www.esa.int/gsp/ACT/doc/ARI/ARIStudyReport/
ACTRPTMADARI054109bSpaceWebsGlasgow.pdf.
D.J.McKenzie andM.P.Cartmell.Ontheperformance of a motorized tether using
a ballistic launch method. 55th International Astronautical Congress, Vancouver,
Canada,(IAC04IAA3.8.2.10), Oct. 48 2004.
W.I. McLaughlin. Walter Hohmann’s roads in space. Space Mission Architecture,
2000. URL http://www.lpl.arizona.edu/undergrad/classes/
fall2007/Hubbard_195a10/Hohmann_bio.pdf.
J.K. Miller and E.A. Belbruno. A method for the construction of a Lunar transfer
trajectory using ballistic capture. Spaceflight Mechanics: Advances in the Astronautical
Sciences, AAS 91100, 75(1):97–109, 1991.
A.K. Misra and V.J. Modi. Threedimensional dynamics and control of tether
connected nbody systems. Acta Astronautica, 26(2):77 – 84, 1992. ISSN 00945765.
doi: 10.1016/00945765(92)90048N.
Y. Miyazaki and Y. Iwai. Dynamics model of solar sail membrane. In ISAS 14th
Workshop on Astrodynamics and Flight Mechanics, 2004. URL http://www.
hayabusa.isas.jaxa.jp/kawalab/astro/2004contents_e.html.
V.J. Modi and A.K. Misra. On the deployment dynamics of tether connected two
body systems. Acta Astronautica, 6(9):1183–1197,September 1979. doi: 10.1016/
00945765(79)90064X.
V.J. Modi, Geng ChangFu, A.K. Misra, and Da Ming Xu. On the control of the
space shuttle based tethered systems. Acta Astronautica, 9(67):437–443, 1982.
doi: 10.1016/00945765(82)900741.
V.J. Modi, P.K. Lakshmanan, and A.K. Misra. On the control of tethered satellite
systems. Acta Astronautica, 26(6):411–423, June 1992. doi: 10.1016/
00945765(92)90070Y.
178
V.J. Modi, S. Bachmann, and A.K. Misra. Dynamics and control of a space station
based tethered elevator system. Acta Astronautica,29(6):429–449,June1993. doi:
10.1016/00945765(93)90036V.
H. Moravec. A nonsynchronous orbital skyhook. The Journal of the Astronautical
Sciences, 25(4):307–322, 1977. URL http://www.frc.ri.cmu.edu/~hpm/
project.archive/1976.skyhook/papers/scasci.txt.
T. Nakano, O. Mori, and J. Kawaguchi. Stability of spinning solar sailcraft containing
a huge membrane. AIAA Guidance, Navigation, and Control Conference
and Exhibit, 2005. AIAA20056182.
S. Nakasuka, T. Aoki, I. Ikeda, Y. Tsuda, and Y. Kawakatsu. “furoshiki satellite”
– a large membrane structure as a novel space system. Acta Astronautica, 48:
461–468, March 2001. doi: 10.1016/S00945765(01)00056X.
A.H. Nayfeh and D.T. Mook. Nonlinear Oscillations. John Wiley and Sons, New
York, 1979.
T.
S. No and J. E. Cochran Jr. Dynamics and control of a tethered
flight vehicle. Journal of Guidance, Control, and Dynamics, 18(1):66–
72, 1995. URL http://www.aiaa.org/content.cfm?pageid=318&volume=18&
issue=1&pubid=23&paperid=56658.
Zoe Parsons. Lunar perturbations of a supersynchronous GEO transfer
orbit in the early orbit phase. Master’s thesis, Cranfield University,
2006. URL http://dspace.lib.cranfield.ac.uk:8080/bitstream/1826/
1767/1/ZParsons_ThesisMSc.pdf.
M. Pascal, A. Djebli, and L. El Bakkali. Laws of deployment/retrieval in tether
connected satellites systems. Acta Astronautica, 45(2):61–73, July 1999. doi:
10.1016/S00945765(99)001150.
179
M. Pascal, A. Djebli, and L. El Bakkali. A new deployment/retrieval scheme for a
tethered satellite system, intermediate between the conventional scheme and the
crawler scheme. Journal of Applied Mathematics and Mechanics, 65(4):689–696,
2001. doi: 10.1016/S00218928(01)000739.
Jerome Pearson. The orbital tower: A spacecraft launcher using the Earth’s rotational
energy. Acta Astronautica, 2(910):785 – 799, 1975. ISSN 00945765. doi:
10.1016/00945765(75)900211.
LindaPetzold.Automaticselection of methodsforsolving stiff and nonstiff systems
of ordinary differential equations. SIAM Journal on Scientific and Statistical
Computing, 4(1):136–148, 1983. doi: 10.1137/0904010.
Ary PizarroChong and Arun K. Misra. Dynamics of multitethered satellite formations
containing aparentbody. Acta Astronautica, 63(1112):1188 – 1202, 2008.
ISSN 00945765. doi: 10.1016/j.actaastro.2008.06.021.
Dominic V. Rosato. Rosato’s Plastics Encyclopedia and Dictionary. Hanser Publishers,
Munich, January 1993.
Shane David Ross. Cylindrical manifolds and tube dynamics in the restricted three
body problem. Ph.D. thesis, California Institute of Technology, 7th April 2004.
URL http://etd.caltech.edu/etd/available/etd05182004154045/.
K. Sorensen. Conceptual design and analysis of an MXER tether boost station. In
37th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, 2001.
URL www.tethers.com/papers/AIAA_20013915_MXER.pdf.
K.Sorensen.Hyperbolic injection issuesforMXER tethers.In 39th AIAA/ASME/SAE/
ASEE JointPropulsionConference,2003. URL http://www.tethers.com/
papers/AIAA_20035221.pdf.
T.H.Sweetster. An estimate of theglobal minimum V
needed for theEarthMoon
transfer. In Spaceflightmechanics1991;Proceedings of the1stAAS/AIAAAnnual
180
Spaceflight MechanicsMeeting, Houston, TX, Feb.1113, 1991. Pt. 1(A9317901
0513), p. 111120., pages 111–120, 1991. URL http://adsabs.harvard.edu/
abs/1991sfm..proc..111S. AAS/AIAA Paper No. 91101.
Victor G. Szebehely and Hans Mark. Adventures in Celestial Mechanics. Wiley
VCH, 2nd edition edition, 2 1998. ISBN 9780471133179. URL http://amazon.
co.uk/o/ASIN/0471133175/.
G. Tibert and M. G¨ardsback. Spacewebs. Technical report, KTH and ESA
ACT, 2006. URL http://www.esa.int/gsp/ACT/doc/ARI/ARIStudyReport/
ACTRPTMADARI054109aSpaceWebsKTH.pdf.
K.E. Tsiolkovski. Grezi o zemle i nene (daydreams of Heaven and Earth).
In U.S.S.R. Academy of Sciences edition, page 35, 1959. URL http://www.
daviddarling.info/encyclopedia/S/space_elevator.html.
University of Glasgow Press Release. Students in a spin after ESA commission
spaceweb. Online, 07 Apr 2009. URL http://www.gla.ac.uk/news/headline_
115222_en.html.
J. Valverde and G.H.M. van der Heijden. Magneticallyinduced buckling of a whirling
conducting rod with applications to electrodynamic space tethers, 2008. URL
http://www.citebase.org/abstract?id=oai:arXiv.org:0810.3668.
G.H.M. van der Heijden. Modelocking in nonlinear rotordynamics. Journal of
Nonlinear Science, 5(3):257–283, May 1995. doi: 10.1007/BF01212957.
Robert C. Weast. CRC Handbook of Chemistry and Physics. CRC Press, 66th
edition, June 1985. ISBN 0849304652,page E43.
Wikipedia. Space debris — Wikipedia, the free encyclopedia, September
2008. URL http://en.wikipedia.org/w/index.php?title=Space_debris&
oldid=239761439. [Online; accessed 28September2008].
181
Paul Williams, Chris Blanksby, and Pavel Trivailo. Tethered planetary capture:
Controlled maneuvers. Acta Astronautica, 53(410):681–708, 2003. doi: 10.1016/
S00945765(03)800292.
Paul Williams, Chris Blanksby, Pavel Trivailo, and Hironori A. Fujii. Inplane
payload capture using tethers. Acta Astronautica, In Press, Corrected Proof:–,
2005. doi: 10.1016/j.actaastro.2005.03.069. AA57772.
Wolfram. Mathematica NDSolve FAQ, 01 May 2007. URL http://support.
wolfram.com/mathematica/mathematics/numerics/ndsolvereferences.en.
html. Accessed on 20th June 2007.
Brian Wong and Arun Misra. Planar dynamics of variable length multitethered
spacecraft near collinear Lagrangian points. Acta Astronautica, 63(1112):1178 –
1187, 2008. ISSN 00945765. doi: 10.1016/j.actaastro.2008.06.022.
Xu Xu, M. Wiercigroch, and M.P. Cartmell. Rotating orbits of a parametrically
excited pendulum. Chaos, Solitons & Fractals, 23(5):1537 – 1548, 2005. ISSN
09600779. doi: 10.1016/j.chaos.2004.06.053.
S.W. Ziegler. The Rigid Body Dynamics of Tethers in Space. Ph.D. thesis, Department
of Mechanical Engineering, University of Glasgow, 2003.
S.W.Ziegler andM.P.Cartmell.Usingmotorizedtethersforpayload orbitaltransfer.
Journal of Spacecraft and Rockets, 38(6):904–913, 2001.
182
Glossary
Notation Description
LEO Low Earth Orbit
MEO Medium Earth Orbit
MMET Motorized Momentum Exchange Tether
MXER Momentum eXchange/Electrodynamic Reboost
WSB Weak Stability Boundary
183
184
Appendix A
Equations of motion with inclination
Theequationsof motiontheMMET(statorarmnot included) forthegeneralised coordinates R,
are:
R
equation:
µMtotal
+¨ 2
R2 RMtotal 
R
cos 2().2 +.Mtotal =0 (A.1)
.
equation:
.
R
cos() 2cos()R.+
R
cos()
.
¨ 2.
.sin() MFacility =0 (A.2)
185
.
equation:
.
2 2.
R2 cos()sin()MFacility+
R
R.+ R¨MFacility =0 (A.3)
.
equation:
.
.
1
..
2cos2()sin(2 ).2 + 4.
a
sin( )2 cos( ) .
sin(2)+.2cos(2 )sin()cos2()+cos()cos( )sin(2) .
8
...
4cos( ).
a
sin()+2.
sin(2)sin( )sin()+2cos2()sin(2 )sin2()+sin(2)sin(2)sin( )
1
¨
(4Mpayload + Mtether)L2 +
192
128
.
(3Mpayload + Mtether)+
48 2cos()
.
¨cos 2()+2
.
¨cos 2()cos( )
.
¨sin(2)sin()2. .
sin(2)+.
(cos()sin(2)+cos(2)cos( )sin()) +
¨
.
sin(2)sin( )+. .
sin(2)sin()sin( )+.
.
cos( ) .
sin(2).
2sin()cos2()+cos()cos( )sin(2) +2cos(2).
sin( )
(4Mpayload + Mtether))L2 + R
µ
cos()sin( )
1

1
Mpayload+
3/23/2
(L2 2R
cos()cos( )
L
+ R2)(L2 +2R
cos()cos( )
+ R2)
11
4  Mtether
3/23/2
(L2 4R
cos()cos( )
L
+4R2)(L2 +4R
cos()cos( )
L
+4R2)
= cosa
Forcepayload L(1Eclipse) (A.4)
186
a
equation:
1
((6cos(2)+cos(2(
.
 ))2cos(2 )+cos(2(
.
+ ))+2)sin(2)+8cos(2)cos( )sin(2)).2+
32
16 .
(cos()sin(2)+cos(2)cos( )sin()).8.2 sin(2)sin2( )+8 .
2 sin(2)
.
8.
.
2cos(2)
.
sin( )+.(2cos(2)cos()sin( )sin(2)sin()sin(2 )) (4Mpayload + Mtether)L2+
.
.
128¨(3Mpayload + Mtether)+96 ¨+ cos( )¨
.
+ cos( ).
.
sin()+.
.
cos(). .
sin( )+
.
¨sin()sin( ) (4Mpayload + Mtether) L2+
192
11
R
µ
cos( )sin()  Mpayload+
3/23/2
(L2 2R
cos()cos( )
L
+ R2)(L2 +2R
cos()cos( )
L
+ R2)
11
4  Mtether L
3/23/2
(L2 4R
cos()cos( )
L
+4R2)(L2 +4R
cos()cos( )
L
+4R2)
= Forcepayload L(1Eclipse) (A.5)
The equations use the following substitutions:
Mtotal = MFacility +2Mpayload +2Mtether (A.6)
187
and Eclipse is the binary function:
Eclipse =
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1
0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
R2
.
1cos2()cos2()
.
= R2
Earth
.
1+
R
cos()
RSun 2
AND
cos()cos()
>
0
.
otherwise
(A.7)
Note: the inclination generalised coordinate is
.
theninth letterof theGreek alphabet(iota). The firstderivative is @
= . and the
@
second derivative is @2.
=¨.
@t2
Appendix B
Spacewebs – additional graphics
B.1
Case1 –CoMplotsof stability whileincreasing
web mass
0.3
0.2
0.1
0
0.1
0.2
0.3
ZCoM Position
Z
CoM Position
0.4
0.2
Y
0
0.2
0.4
Y
0.3 0.2 0.1 0 0.1 0.2 0.3
0.4 0.2 0 0.2 0.4
Figure B.1: mW
eb
=27kg
Figure B.2: mW
eb
=33.75kg
188
Z
ZCoM Position
CoM Position
0.75 0.5 0.25 0 0.25 0.5 0.75
0.75
0.5
0.25
0
0.25
0.5
0.75
Y
0.4
0.2
0
0.2
0.4
Y
0.4 0.2 0 0.2 0.4
Figure B.3: mW
eb
=37.8kg
Figure B.4: mW
eb
=40.5kg
Z
CoM Position Z
CoM Position
3
2
1
0
1
2
3
15
10
5
Y
0
5
10
1510 5 0 5
3210 1 2 3
Y
10 15 20
Figure B.5: mW
eb
=44.82kg
Figure B.6: mW
eb
=47.25kg
ZZ
CoM Position CoM Position
40
30
20
10
Y
30
20
10
0
10
20
Y
10
20
30
40 30 20 10 0
0
30
10 20 30
302010 0 10 20 30
Figure B.7: mW
eb
=54kg
Figure B.8: mW
eb
=67.5kg
189
B.2 Case1 –CoMplotsof stability whileincreasing
robot mass
40
20
0
20
40
ZCoM Position
Z
CoM Position
20
0
20
40
Y
Y
40 20 0 20 40 4020 0 20 40
Figure B.9: Mrobot
=1kg
Figure B.10: Mrobot
=10kg
Z
CoM Position
Z
CoM Position
40
1.5
1
20
0.5
Y
Y0
0.5
20
1
1.5
40
1.5 1 0.5 0 0.5 1 1.5
4020 0 20 40
Figure B.11: Mrobot
=15kg
Figure B.12: Mrobot
=16kg
ZZ
CoM Position CoM Position
0.8
0.2
0.6
0.6
0.4
0.4
Y
0.2
0
Y0
0.2
0.2
0.4
0.4
0.6
0.6
0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.6 0.4 0.2 0 0.2 0.4 0.6
Figure B.13: Mrobot
=20kg
Figure B.14: Mrobot
=100kg
190
B.3 CoM plots of stability while changing robot
paths
Z
ZCoM Position
CoM Position
20
10
Y
0
10
20
0.4 0.2 0 0.2 0.4
2010 0 10 20
0.4
0.2
0
0.2
0.4
Y
Figure B.15: Case 1
Z
CoM Position
15
10
5
0
5
10
15
15105 0 5 10 15
Figure B.17: Case 3
Z
CoM Position
30
20
10
0
10
20
30
40302010 0 10 20 30
Figure B.16: Case 2
Z
CoM Position
10
5
Y
0
Y
5
10
1050 510
Figure B.18: Case 4
Z
CoM Position
20
10
Y
0
Y
10
20
30
3020 0 20
10 10
Figure B.19: Case 5 Figure B.20: Case 6
191
B.4 Case 1 – CoMplots of stability while decreasing
robot velocity
0.3 0.2 0.1 0 0.1 0.2 0.3
0.3
0.2
0.1
0
0.1
0.2
0.3
Y
CoM PoZsition
Figure B.21: Vrobot = 10.0m/s
0.3 0.2 0.1 0 0.1 0.2 0.3
0.3
0.2
0.1
0
0.1
0.2
0.3
Y
CoM PoZsition
Figure B.22: Vrobot = 2.0m/s
0.4 0.2 0 0.2 0.4
0.4
0.2
0
0.2
0.4
Y
CoM PoZsition
Figure B.23: Vrobot = 1.0m/s
1 0.5 0 0.5 1
1
0.5
0
0.5
1
Y
CoM PoZsition
Figure B.24: Vrobot = 0.2m/s
1.5 1 0.5 0 0.5 1 1.5
1.5
1
0.5
0
0.5
1
1.5
Y
CoM PoZsition
Figure B.25: Vrobot = 0.1m/s
192
X
100
0
50
100
100
50
50
50
100
Z
0
00
50
5050
100
100100
50
0
50
Y
100
Figure B.26: Case 1 – 3 robots moving symmetrically round the web perimeter
193
X
100
0
50
100
100
50
50
50
100
Z
0
00
50
5050
100
100100
50
0
50
Y
100
Figure B.27: Case 6 – 1 robot moving asymmetrically along first subspan
194
Appendix C
Spacewebs – additional data
Data generated during the first run phase – initially examining the dynamics of
the spaceweb.
The units for the variables are as follows:
Mweb Msat MRobot psi psidot k CoM disp
kg
kg
kg
degrees
radians/s
N/m
m
The initial conditions of all other parameters are:
L
=100.0m
; Mfacility =100.0kg
; eccent
=0.0;
tend =100.0s
; R
=6.578* 106 m
; K
=122.69
K
is the spring stiffness, given in Equation 6.15. This is derived from the Young’s
Modulus, E
=3* 1010 and the CSA of each strand of 4 * 106 .
The webs are configured with3 web sections(i.e. 2divisions).
The angulardisplacements aregiven with one angle only: . This may be converted
to give the angles of the three subspans such that: {0+ ,
120. ,
240.  }.
The three angular rates are equal, and are quoted once in the following tables..
195
Datageneratedduring thestatistical investigation intostability(seeSection6.8).
run A:Mweb. B:Msat C:Mrobot D:psi E:psidot Max CoM
E01 27. 10 2 1 0.1 0.35
E02 33.75 10 2 1 0.1 0.41
E03 37.8 10 2 1 0.1 0.85
E04 40.5 10 2 1 0.1 0.39
E05 44.82 10 2 1 0.1 3.74
E06 47.25 10 2 1 0.1 20.44
E07 54. 10 2 1 0.1 34.82
E08 67.5 10 2 1 0.1 39.40
E09 135. 10 2 1 0.1 46.49
E10 270. 10 2 1 0.1 51.03
Table C.1: Influence of web mass on maximum CoM
run symmetrical A:Mweb B:Msat C:Mrobot D:psi E:psidot Max CoM
Case 1 yes 40.5 10 10 1 0.1 0.46
Case 2 no 40.5 10 10 1 0.1 22.61
Case 3 no 40.5 10 10 1 0.1 18.52
Case 4 yes 40.5 10 10 1 0.1 10.95
Case 5 no 40.5 10 10 1 0.1 38.24
Case 6 no 40.5 10 10 1 0.1 30.99
Table C.2: Influence of symmetrical robot movement on maximum CoM
run A:Mweb B:Msat C:Mrobot D:psi E:psidot time Robot speed Max CoM
FG1 40.5 10 2 1 0.1 10 10.0 0.36
FG2 40.5 10 2 1 0.1 50 2.0 0.36
FG3 40.5 10 2 1 0.1 100 1.0 0.40
FG4 40.5 10 2 1 0.1 500 0.2 1.40
FG5 40.5 10 2 1 0.1 1000 0.1 1.70
Table C.3: Influence of simulation time (effectively robot speed across web) on
maximum CoM
196
run A:Mweb B:Msat C:Mrobot D:psi E:psidot Max CoM
Frobot01 135 10 0.1 1 0.1 45.09
Frobot02 135 10 1 1 0.1 45.89
Frobot03 135 10 5 1 0.1 46.66
Frobot04 135 10 10 1 0.1 46.08
Frobot05 135 10 12 1 0.1 50.53
Frobot06 135 10 15 1 0.1 48.44
Frobot07 135 10 16 1 0.1 1.75
Frobot08 135 10 17 1 0.1 1.08
Frobot09 135 10 18 1 0.1 0.86
Frobot10 135 10 19 1 0.1 0.82
Frobot11 135 10 20 1 0.1 0.77
Frobot12 135 10 21 1 0.1 0.66
Frobot13 135 10 50 1 0.1 0.67
Frobot14 135 10 100 1 0.1 0.75
Table C.4: Influence of robot mass on maximum CoM
197
run A:Mweb B:Msat C:Mrobot D:psi E:psidot F: CoM posn predicted Diff %
1 27 2 1 0.01 0.1 11.49 12.29 107%
2 270 2 1 0.01 0.1 49.41 46.71 95%
3 27 10 1 0.01 0.1 0.00 2.70 83118%
4 270 10 1 0.01 0.1 52.76 53.56 102%
27 2 10 0.01 0.1 28.09 25.39 90%
6 270 2 10 0.01 0.1 53.20 54.00 102%
7 27 10 10 0.01 0.1 0.00 0.80 19141%
8 270 10 10 0.01 0.1 55.75 53.05 95%
9 27 2 1 0.1 0.1 11.28 10.16 90%
270 2 1 0.1 0.1 49.83 47.37 95%
11 27 10 1 0.1 0.1 0.03 0.12 369%
12 270 10 1 0.1 0.1 52.79 51.67 98%
13 27 2 10 0.1 0.1 27.49 26.72 97%
14 270 2 10 0.1 0.1 52.99 51.65 97%
27 10 10 0.1 0.1 0.04 0.05 111%
16 270 10 10 0.1 0.1 55.61 54.83 99%
17 148.5 6 5.5 0.055 0.55 47.15 31.10 66%
18 148.5 6 5.5 0.055 0.55 47.15 31.10 66%
19 148.5 6 5.5 0.055 0.55 47.15 31.10 66%
148.5 6 5.5 0.055 0.55 47.15 31.10 66%
21 27 2 1 0.01 1 49.41 46.71 95%
22 270 2 1 0.01 1 49.86 50.66 102%
23 27 10 1 0.01 1 0.00 0.80 24649%
24 270 10 1 0.01 1 0.00 2.70 64127%
27 2 10 0.01 1 33.10 33.90 102%
26 270 2 10 0.01 1 49.41 46.71 95%
27 27 10 10 0.01 1 49.83 47.13 95%
28 270 10 10 0.01 1 55.77 56.57 101%
29 27 2 1 0.1 1 49.41 48.63 98%
270 2 1 0.1 1 49.82 48.69 98%
31 27 10 1 0.1 1 0.03 1.09 3346%
32 270 10 1 0.1 1 0.03 0.74 2287%
33 27 2 10 0.1 1 33.11 31.98 97%
34 270 2 10 0.1 1 49.41 48.63 98%
27 10 10 0.1 1 0.00 0.77 18329%
36 270 10 10 0.1 1 55.76 54.63 98%
Table C.5: Statistical investigation into stability
198
Term Effect Stdized Effect SumSqr % Contribtn
AMweb 27.21 7182.9 37.72
BMsat + 16.57 2447.8 12.85
AE + 15.32 1869.6 9.82
BCE 14.34 1656.8 8.70
CMrobot 10.90 872.3 4.58
BC 9.93 666.4 3.50
BE + 8.69 686.3 3.60
AB 7.21 435.0 2.28
ACE 6.98 346.3 1.82
CE 5.04 232.9 1.22
AC 4.66 193.2 1.01
ABE + 3.76 118.3 0.62
ABC 3.75 118.6 0.62
BDE + 3.37 93.4 0.49
ADE 3.29 117.6 0.62
ABCDE 3.19 85.1 0.45
ACD 3.17 99.5 0.52
BCD + 3.15 102.7 0.54
CDE + 3.12 92.4 0.49
ABDE 3.12 81.5 0.43
CD + 3.09 110.1 0.58
Dpsi + 3.08 79.6 0.42
BCDE + 3.08 83.7 0.44
ACDE 3.04 89.0 0.47
ABCD 3.02 61.5 0.32
DE + 3.00 75.6 0.40
ABD 2.97 74.0 0.39
AD 2.93 98.7 0.52
BD + 2.84 34.2 0.18
Epsidot 1.48 18.3 0.10
Lack of fit 9.6 0.05
Pure error 0.0 0.00
Table C.6: Results of the statistical investigation into stability
The ‘Effect’ term incolumn2 isareflectiononthestabilising(+) ordestabilising
()effect of the term.
199
Appendix D
Material strengths
The information in Table D.1 is provided courtesy of the Island One Society at
http://www.islandone.org/LEOBiblio/SPBI1MA.HTM, and is compiled from data
from the following sources: [Brown, 1989, pp 1418], [Rosato, 1993, p 638], [Weast,
1985, p E43].
1aramid fiber
2gelspun polyethylene
3gelspun polyethylene
4plastic fiber, Zylon is brand name of PBO
5theoretical data
200
Material Property
tensile Young’s density () characteristic speed of
strength
(P)
modulus
(Y)
velocity
.
2P/.
sound in thin
rods .
Y/.
[GPa] [GPa] [kg/m3] [m/s] [m/s]
steel 15 200 7900 5031125 5032
aluminum alloys 0.10.7 72 2700 272720 5270
titanium alloys 0.61.3 110 5000 490721 4690
berylium fiber 3.3 310 1870 1879 12870
boron fiber 3.5 400 2450 1690 12778
fused silica 73 2200 5760
pyrex glass 62 2320 5170
Eglass fiber 2.4 72.4 2540 1375 5339
Sglass fiber 4.5 85.5 2490 1901 5860
Kevlar 491 3.6 130 1440 2236 9502
Spectra 1000 fiber2 3.0 170 970 2487 13239
Spectra 2000 fiber3 3.51 970 2690 
PBO (Zylon)4 5.8 280365 15601580 27102727 1339715199
carbon fiber 16.5 250830 1850 10402651 1160021200
buckytube cable5 150 630 1300 15191 22014
Table D.1: Properties of selected materials
201
Appendix E
Motor properties
Specification
List Price: $910
Horsepower 4hp
RPM 4150
mass 23kg
Armature Volts 48V
Torque 6.86Nm
Full Load Current 72.0A
Table E.1: Properties of GE Motor 5BC49JB1115
The information in Table E.1 is provided courtesy of the GE industrial products
guide at:
www.geindustrial.com/cwc/marketing/Motors/catalog/CrossReference.pdf
and is compiled from data from the following additional source:
www.emotorstore.com/productdetail.asp_Q_brandID_E_1_A_catID_E_23_A_
subCatID_E_354_A_productID_E_583_A_skuID_E_30698
The motor is a GE Industrial Motor, Product Number D379, Model Number
5BC49JB1115. ThisisaDC motor, used topower electric vehicles(commonlygolfcarts).
The specifications are used in Section 4.5.1 to gain an approximate mass, torque and
rotational speedforan equivalent spacebased motor. Thetorqueisquoted as81.0ozft,
and converted to Nm.
202