Quasi-satellite dynamics in formation flight

3.7 MB·

by  Seppo Mikkola

Mon.Not.R.Astron.Soc.000,1–8() Printed5January2016 (MNLATEXstylefilev2.2)
Quasi-satellite dynamics in formation flight
Seppo Mikkola1⋆ & Claudiu-Lucian Prioroc1,2†
6 1Department of Physics and Astronomy, University of Turku,
1 Tuorla Observatory, V¨aisa¨l¨antie 20, Piikki¨o , FI 21500, Finland
0 2GMV, Calea Victoriei 145, Bucuresti 010072, Romania
2

n
5January2016
a
J

1
ABSTRACT

Thequasi-satellite(QS)phenomenonmakestwocelestialbodiestoflyneareachother
]
P (Mikkola et al.2006)andthateffectcanbeusedalsotomakeartificialsatellitesmove
E in tandem. We consider formation flight of two or three satellites in low eccentricity
. near Earth orbits. With the help of weak ion thrusters it is possible to accomplish
h
tandem flight. With ion thrusters it is also possible to mimic many kinds of mutual
p
force laws between the satellites. We found that both a constant repulsive force or an
-
o attractive force that decreaseswith the distance are able to preservethe formationin
r which the eccentricities cause the actual relative motion and the weak thrusters keep
t
s the mean longitude difference small. Initial values are important for the formation
a
flight but very exact adjustment of orbital elements is not important. Simplicity is
[
oneofourgoalsinthisstudyandthisresultisachievedatleastinthewaythat,when

1 constant force thrusters are used, the satellites only need to detect the directions of
v the other ones to fly in tandem. A repulsive acceleration of the order of 10−6 times
6 the Earthattraction,is enoughto effectively eliminate the disruptive effects of allthe
6 perturbations at least for a timescale of years.
0
0 Key words: celestial mechanics – planets and satellites
0
.
1
0
6
1 INTRODUCTION ation. In the case of two satellites the force between the
1
satellites can been mimicked using ion thrusters. Both an
: In an attempt to understand better the quasi-satellite phe-
v attractive force 1/∆2 and repulsive constant force pro-
nomenon (Mikkola et al. 2006) we tested what would hap- ∝
i ducesimilar effects,which may look somewhat unexpected.
X peniftheforcebetweentwoco-orbitalbodiesweredifferent
Above ∆ means the distance of the satellites which we as-
r from the normal Newtonian gravity. We found that many
sumeto bequitesmall compared with thesize of theorbit.
a forces can produce similar relative motion of the bodies.
Werestrictourconsiderationtoloweccentricityorbits.High
Evenrepulsiveforcesresultintandemmotionofthebodies.
eccentricityorbitsaregenerallymorecomplicatedforforma-
At that point it became clear that constant weak repulsive
tion flying(Roscoe et al. 2013).
force could be used to keep artificial satellites flying near
The simplest of all methods to make satellites to fly in
eachother.Thispaperstudiesthatphenomenonmainlynu-
tandem is to put them into precisely same orbit to some
mericallybutwealsopresentsomesimpleanalyticalconsid-
distancefrom eachother.Anexampleofsucharealsystem
erations.
isNASA’sGrail mission in which twospacecrafts are flying
In quasi satellite motion an asteroid moves around the
in tandem orbits around the Moon to measure its gravity
Sun co-orbitally with a planet and remains near the planet
field in detail.
suchthatintherotatingcoordinatesystemitlookslikemov-
ing around the planet in a retrograde orbit. In that system This works perfectly as long as one can consider the
the mutual force of the bodies is attractive. On the other gravitational field to be rotationally symmetrical (like the
handitiswellknownthataspacecraft chasingasatellitein Earth’s field modeled without tesseral harmonics) in which
the same orbit must brake to catch it up. Thus, somewhat case the distance between the satellites varies only due to
counter intuitively, tandem flight of two satellites moving eccentricity caused speed differences. The simplest way to
around a central body in its gravitational field seems to be achieve tandem flight in that situation is possibly making
possible both with attractive or repulsive mutual acceler- ionthrusterstopointtowards,orawayfrom,theothersatel-
lite(s) and thisseems to work according tooursimulations.
One can keep the satellites to move in the same fashion
⋆ E-mail:seppo.mikkola@utu.fi asthenaturalquasi-satellites.Evenneartriangleconfigura-
† E-mail:[email protected] tion for three satellites is possible without any complicated
2 Mikkola and Prioroc
control, just constant acceleration or an acceleration that the mean anomaly, respectively. For the position vector we
dependson thedistanceof thesatellites. Nootherinforma- use r and for velocity v. In the used units the velocity and
tion is necessary but the relative directions (and distances momentum are the same i.e. v = p. The semi-major-axes
if the acceleration needs to depend on the distance). Long of the satellites are only very slightly different (and equal
lasting thrustersare necessary, butsuch exist as shown e.g. on the average) and therefore we can in most consideration
byESA’sSMART-1mission from anearEarth orbit tothe takethemtobepractically samealthoughthedifferenceδa
Moon. has an important role in some analytical considerations.
The very basic method to consider satellite relative Weneedalsothedifferenceofthemeanlongitudesλ=
motions near each other and/or rendezvous is naturally M + ω + Ω. We consider (almost) coplanar orbits which
the Clohessy-Wiltshire equations (Clohessy and Wiltshire means that Ω = Ω to high precision and therefore the
1 2
1960). To get satellites flying in tandem has been consid- difference of λ’s is essentially the same as the difference of
ered in numerous publications e.g. precise orbital elements theanglesM+ω.Thisquantityisoneofthemostimportant
adjustments have been discussed in e.g. (Alfriend et al. ones in our analytical consideration of the stability of the
2001; Beichman et al. 2004; Schaub and Alfriend 2001; tandem flight. Weuse for it thenotation
Gim and Alfried 2005).
θ=λ λ =M M +ω ω , (1)
Kristiansen et al.(2010)consideredtheformationflight 2− 1 2− 1 2− 1
usingdifferentialapproximationsofthepureKeplermotion, where M˙k =nk =1/ (a3k/m⊕) are the mean motions. An
alsoPrioroc and Mikkola(2015)discussthecomputationof otherimportant quanptityis theeccentricity vector
relative motion using simple methods. Roscoe et al. (2013)
discuss formation flight in highly eccentric orbits using dif- e=v×(r×v)/m⊕−r/r, (2)
ferential mean orbital elements as design variables and the and especially their difference
Gim and Alfried (2005) state transition matrix for relative
e=e e , (3)
motionpropagation.Useofacoordinatetransformationwas 2− 1
consider by Vallado and Alfano (2014). ewhich is a kindof relative eccentricity vector.
We first consider the problem in terms of a simple ap-
proximateanalyticaltheory.Theassumptionsinthattheory
are based on simulation results since a complete analytical 2.2 C-W/Hill theory
solutionseemslaborious.Aftertheanalyticalconsiderations
Clohessy and Wiltshire (1960) presented a theory of satel-
we present results from numerical simulation. The gravity
lite rendezvous (C-W hereafter). This theory used the vari-
field that we used in the tests was the spherical function
ational equations of thetwo-body problem
expansion, following Mikkola, Palmer and Hashida (2002),
plus Luni-Solar terms, but we also studied the simple J2 x˙ =w, w˙ =−m⊕(x/r3−3r·xr/r5), (4)
field model. The ion thruster effects we model using differ-
where x is the variation of the position and w = x˙ is the
ent model potentials could be engineered with adjustable
variation of velocity. For a circular orbit in therotating co-
thrusts. However, the constant force model is the simplest
ordinatesystem, wherer isconstant, theequations become
one and it seems to work.
linearandeasily solvable.Muchearlier Hillhadstudiedthe
Finallyweremarkthatthisisnotanengineeringpaper,
Lunarorbitinasomewhatsimilarway(Hill 1878).Thedif-
we just study thedynamics.
ferenceinthesetreatmentsisthatinthelatteronethereare
additionalnonlineartermsintheequations(4)ofmotion.As
showne.g.byPrioroc and Mikkola(2015)theC-W-solution
2 ANALYTICAL CONSIDERATIONS andinfactasomewhatmoregeneralonecanbeobtainedby
differentiating the two-body equations with respect to the
2.1 Units and orbital elements orbital elements. Let q k =,1,2..,6 be a set of elements
k
e.g.q=(M ,i,Ω,ω,a,e),whicharethemeananomalyM at
Wechoose theunits such that 0
epoch, inclination, ascending node, argument of pericentre,
G=1, R⊕ =1, m⊕ =1 semi-major-axis and eccentricity, respectively. The solution
for thevariational equation takes theform
where G=the gravitational constant, R⊕=radius of the
Earth (actually 6378.137 km) and m⊕= the mass of the x= ∂r δq and w= ∂vδq . (5)
Earth. The unit of time becomes such that the formal X∂qk k X∂qk k
period of a satellite with semi-major-axis a = 1 would The constants δq represent the element variations. A par-
k
be = 2π. In this system one time unit in seconds is = ticular solution of equation (4) is
806.811064922699 sec and thelength of a day is =107.088.
x =2r 3tv, (6)
Despite these choices we write in some formulae the radius a
−
of the Earth and its mass explicitly visible. In illustration, whichisrelatedtothevariationofthesemi-major-axis. Be-
tomakethemeasier tounderstand,weusealsootherunits, cause the time t appears as a factor of v it is clear that
e.g. kilometers. For orbital elements we use the standard without some regulating additional force the distance be-
notations tween the bodies would increase continuously unless the
semi-major-axis difference is precisely zero.
a, e, i, ω,Ω, M,
Additional information about the stability of the Hill-
which are the semi-major axis, eccentricity, inclination, ar- type and distant retrograde orbits can also be found in
gument of perihelion, longitude of the ascending node and Jackson (1913),Scheeres, (1998) and Mikkola et al. (2006).
Quasi-satellite dynamics in formation flight 3
2.3 Secular dynamics 150
2.3.1 Assumptions based on numerical results
100
Here we present some numerically obtained facts that can
beusedasassumptionsinouranalyticalconsiderations.The 50
force model we used consisted Earth potential model with Y
∆
the J terms and a weak ion thruster acceleration that was
2 0
either a repulsive constant acceleration or, for comparison,
an attraction that had the1/∆2 distance dependence.Also
-50
other models were brieflyconsider (see below).
The following was found tobe truein our simulations:
-300 -200 -100 0 100 200 300
-The orbital elements vary periodically but their mean ∆X
valuesremainthesameandtheJ2(andhigherorder)effects Figure1.Therelativemotionofthesatellitesforthecaseswhen
have minor influencein the relative motion of the satellites θ=±2e,θ=±eandθ=0inthecoordinatesystemthatrotates
(see also Prioroc and Mikkola (2015)). withoneofthesatellites(=the central dot=satellite 1).Herethe
-InclinationsiandlongitudesofascendingnodeΩremain eccentricity ofsatellite2ise=0.01,whilethesatellite1has ec-
veryclosetoeachotheriftheyaresoinitially.Thisisdueto centricity0.Thedistanceswereconvertedtokilometersassuming
similar precession rates. Thus the satellites remain almost thata=7000km.
co-planar. Also theω’s precess at thesame rate.
-The difference of the eccentricity vectors e remains al- ∆2 =a2((θ+2esin(M))2+e2cos2(M))+O(ε2). (8)
most constant in length. This is important in the theory
e Thegeneralization ofthistothecaseoftwoeccentricorbits
andinaccordancewiththeresultsobtainedinMikkola et al.
is quitesimple (see below).
(2006). See Fig.6 in thispaper.
Forstabletandemflightsitisclearlynecessarythatthe
-When the initial semi-major-axes are nearly the same
mean longitude difference θ indeed remains small. This is
theyremainsoaslongasthesatellitesflyintandem.Dueto
what we studynow.
the thruster effect the values fluctuate such that the mean
If both satellites are supposed to use the ion thrusters
angular speeds remain the same. Thus in an approximate
with equal accelerations, it is possible to model the system
theory we can take the (mean) semi-major-axes to be the
mathematically in terms of a Hamiltonian
same, only their difference varies periodically.
H =K +K R R R , (9)
Those observations can be used as assumptions in the 1 2− 1− 2− 12
next section where the analytical treatment is discussed. where
K =p2/2 m⊕, (10)
2.3.2 Analytical approximations ν ν − rν
andR ’saretheperturbingforce functionswithhigh order
k
Considerfirstasimplecaseinwhichweassumethemotions sphericalharmonics,includingtesseralonesandtherotation
tobepuretwo-bodymotionsandofloweccentricity.Further of the Earth, were included. The force function for both
we assume the motions to be co-planar. satellites was
Asafirstsimpleexampleweconsiderthecaseinwhich
Nx n
ocennetorifctithye.Ionrbtihtesrisotcairtcinuglacroaonrdditnhaeteostyhsetreomneinhwashiachsmthaellceicr-- R=XX Pnmr(nco+s1(θ))(Cnmcos(mψ)+Snmsin(mψ))+Rls.(11)
n=1m=0
cular orbit satellites coordinates are simply = a(0,1) when
Hereψ=φ n t η ,φisthelongitude,readeastwards,in
thedistanceisconsidertobethey-coordinate,theeccentric − E − 0
orbit has approximately thecoordinates 1 theEarthfixedcoordinates,tistime,nE isEarthsrotational
frequency,η initialphase,whileN definestheorderupto
0 x
a(2esin(M)+θ,1 ecos(M)), which terms are taken into account (we used N =36, like
− x
Mikkola, Palmer and Hashida (2002)).Thenumericalcoef-
when eccentricity and θ are consider tobequantities of the
ficients C and S define the gravitational potential. Pm’s
same order of magnitude (we call this order O(ε)) and so n n n
aretheassociated Legendrepolynomials.ThetermR nat-
small that the linearized approximation is valid. ls
urally signifies the effect due to the Sun and Moon (ls =
The differences of the coordinates of the satellites are
Luni-Solar).
thus
If we include only theleading term, theJ term, then
2
(∆x,∆y)=a(2esin(M)+θ, ecos(M))+O~(ε2), (7)
− R = m⊕(J2R⊕2 (1 3(zν)2), (12)
fromwhichweseethattherelativemotionismainlycaused ν rν 2rν2 − rν
bythe eccentricity as long as themean longitude difference
are the J terms (J = 0.0010826299890519 in this nota-
2 2
θ remainssmall.Theeffectofθ isillustrated inFig.1.From
tion). The ‘perturbing’ thruster acceleration is in the term
the above coordinate differences we have for the distance
R =R (∆) (which thusdependsonly on ∆), where
12 12
between thesatellites the form
∆= r r . (13)
2 1
| − |
1 ThecorrespondingformulaeinMikkolaetal.(2006)areiner- The secular (averaged over one period) form of the Kν
rorbutnotused. Hamiltonians in theJ2 approximation are
4 Mikkola and Prioroc
hKν−Rνi=−m2a⊕ν − J2R⊕24ma3ν⊕(1(cid:0)2−e2ν3)s3i/n22.(iν)(cid:1) (14) θ¨≈−a32(∂∂MR2 − ∂∂MR1)=−a62∂∂Rθ, (23)
− in which the last equality follows from thefact that ∂R =
TheperturbingfunctionR12tellshowthethrustersareused. ∂M1
For example if constant acceleration is used, then ∂R = ∂R.Notethattheinclusionofthe ∂R -termisdue
−∂M2 ∂θ ∂M1
to our assumption that both satellites have a thruster. At
R (∆)=ǫ∆, (15)
12 thispointoneseesthattheJ termhasnoeffectinthemo-
2
wheretheconstantǫisthesatellitesacceleration duetothe tion of θ. Howeverit has asmall indirect effect because the
ion thruster. We assume this to be constant (unless other- orbitsarenot Keplerian and thesecular theoryis notaccu-
wise indicated). rate enough. These facts are illustrated in Fig. 2 in which
FollowingMikkola et al.(2006)forco-planarorbitsone onecan seetheeffect of J2 ascompared with thepuretwo-
can derivetheapproximation body motion. In these experiments the thrusters perturbed
thesatellite motions in the same way.
∆2 a2(e2cos(w)2+(θ+2esin(w))2) (16) In addition we really have the term ω¨ ω¨ , but here
≈ 2− 1
where e e the terms nearly cancel each other and they are very small
=O(J2),sothattheequation(23)isagoodapproximation.
2
e= e2 e1 , w=M +constant. (17) Thesecular perturbation duetothethrustereffect can
| − |
now beobtained by averaging ∆ overthe angle w as
eThusthismore general approximation differs from thesim-
ple equation (8) only in the way that the eccentricity is re- 1 2π
placedbythe‘relativeeccentricity’eandtheanglewdiffers hR12i= 2π Z R12(∆(w))dw (24)
from M by a constant (of two-body motion). 0
e
As said above, for constant mutual force R =ǫ∆. Due to
12
thestructureof theexpression for ∆2 it is clear that R is
2.4 Equations of motion h i
an even function of θ, i.e.
Thesecularperturbingfunctionforthetwo-satellitemotion R =R(0)+R(2)θ2/2+.., (25)
is h 12i 12 12
2 J2m⊕R⊕2 2 3sin2(iν) Averaging over the angle w in Eq.(16) and expanding in
R= − +ǫ<∆>. (18) powers of θ one gets the secular approximation
Xν=1 4a3ν((cid:0)1−e2ν)3/2 (cid:1) R ǫa(1.542e+0.1426θ2/e+O(θ4)), (26)
12
The usual Lagrange equations for the orbital elements h i ≈
a, e, i, ̟, Ω and mean anomaly M can be written in the These equations areenaturally valied only to second order in
form θ since that is the case for Eq.(16). For the equation of θ-
motion we have
2 ∂R
a˙ = ,
na∂M θ¨= 6 ∂R12 −1.711ǫθ, (27)
√1 e2∂R 1 e2 ∂R −a2 ∂θ ≈ ae
e˙ = − − + −
na2e ∂ω na2e ∂M whichshowsthattheθ-motionisharmonicoscillationinthis
e
1 ∂R cos(i) ∂R approximation. In Fig. 4 the periods are compared for the
i˙ = − +
na2√1 e2sin(i)∂Ω na2√1 e2sin(i)∂ω differentmodels(usingEq.(11),theJ2onlymodelandpure
− − two-body motion) and the results are quite close to each
1 ∂R
Ω˙ = other, with only a few percent differences. Thus one may
na2√1 e2sin(i) ∂i
− concludethattheR12 =ǫ∆perturbationmainlydetermines
√1 e2∂R cos(i) ∂R themotion of themean longitude difference θ.
ω˙ = −
na2e ∂e − na2√1 e2sin(i) ∂i In general it is interesting to consider the effect of dif-
−
2 ∂R 1 e2∂R ferent kinds of mutual interactions. For a general power n
M˙ = n − (19) of thedistance onemay write the average of R=ǫ ∆n as
− na ∂a − na2e ∂e n
(20) 1 2π
R =ǫ ∆n(w)dw. (28)
Since the semi-major-axes a and a are equal on the aver- h i n2π Z
1 2 0
age, we write a for this average, and n=1/(a a/m⊕) for The sign of the second derivative of R , with respect
thecorresponding mean motion. p h i
toθ, tells if themean longitude difference θ can stay small.
The derivativeof θ reads
In fact is turnedout, somewhat unexpectedly,that thesec-
θ˙=λ˙ λ˙ =M˙ +ω˙ M˙ ω˙ . (21) ondderivativeisapositivenumberforallnonzeropowersn,
2 1 2 2 1 1
− − −
positiveornegative.ThisresultisillustratedinFig.7.Thus
Due to the smallness (O(J2)) and similarity (ω˙ ω˙ ) of
the perturbations the main term in the equation2 f≈or θ˙1can the tandem flight seems possible with any mutual accelera-
tionofthesatellitesiftheforceisderivablefromapotential
betaken as
of the form ǫ ∆n for (n = 0) where ǫ > 0. In passing we
n n
3 6
θ˙ n n nα (22) mentionthattheperturbingfunctioncaningeneralbeany-
≈ 2− 1 ≈−2 thingoftheformR= +∞ c ∆n,whichmaybeafinite,
n=−∞ n
whereα=(a2 a1)/aandbecauseoftheexpression(above) oratleastconvergent,Pexpressionstillgivingpositivesecond
fora˙ andthef−act thatM ispresentonlyinθ=M M + derivative [R(2) in (25)] at θ = 0. This happens at least if
2− 1 12
ω ω ,we have thecoefficients c are non negative.
2 1 n
−
Quasi-satellite dynamics in formation flight 5
2 2
1.5 1.5
1
1
θ (theta in degrees) 0 .05 θ (theta in degrees) - 00.. 550
-0.5
-1
-1 -1.5
-1.5 -2
95000 95500 96000 96500 97000 97500 98000 98500 99000 99500 100000 0 2000 4000 6000 8000 10000
Time Time
Figure 2. Oscillation of the angle θ in two models: with the Figure 3. Oscillation of the angle θ in two models: with the
entire Earth potential and with J2-term only. After about 2.5 entire Earth potential and with pure two-body model. Here the
years(whichcorrespondsnearthemiddleofthisfigure,inwhich time-span is about three months and the inclination and semi-
thetotal time-spanisabout1.5months)theoscillationsarestill major-axisdifferencearecorrespondinglyi=30o,δa0=1km.
almost the same within plotting precision. Here i = 30o, δa0 =
1km.
2
True
2B
anal
3 NUMERICAL EXPERIMENTS 1.5
1
In this section we first discuss the results with two-satellite
experiments and later this is extended to three satellites 0.5
moviOnguresnsuenmtiearlilcyalinexap(eflraimtteenntesd)cotrnifiarnmgleccleoanrfilygutrhaetiomna.in θ degrees 0
-0.5
point of the theory part: the mean longitude difference θ
fluctuates around zero and remains small. The amplitude -1
depends on the difference in the initial values of the semi- -1.5
major-axis, but the effect of the thruster makes this differ-
-2
0 200 400 600 800 1000 1200 1400
ence fluctuate around zero according to the equation (27) Time
and so the mean values are same. Due to the fact that the
Figure4.Inthedifferentapproximationstheθperiodsaresome-
twoorbitsareverysimilar,alltheperturbationsinthemare
what different, but no more than one could expect. The ’True’
similar and do not affect noticeably the relative motion of curveis fromour mostaccurate modelling(all perturbations in-
thesatellites. cluded), the ’2B’ means two-body motion model without any of
The initial values in the numerical experiments were theEarthharmonicsand’anal’pointstotheresultcomputedus-
produced using two-body formulae. We set initially the in- ingtheaboveequation(27).Thedifferencesintheperiodlength
clinations i = 10◦,30◦,50◦,70◦,90◦ and longitudes of as- areonlyafewpercent.
cending nodes Ω to same value for the satellites. For the
semi-major-axes we used typically one kilometer difference
intheinitialaxisvalue.Theinitialeccentricitiesweresetto InFig.(2)wecomparetheoscillationsoftheθangleina
the values of e = 0.01. The differences in the initial po- system with all theEarth potential term [Eq.(11)]with the
sitions of the satellites were obtained by adding 180◦ to modelthatincludesonlytheJ2-term.Thetimespanshown
the valueof ω and subtracting that amount from themean inthefigurecorrespondsapproximatelyto1.5monthsinthe
anomaly M. For the experiments with three satellites we end of a total interval of nearly 2.6 years. We see that the
used (k 1) 120◦ (k=1, 2, 3) ω differencesand removed θ motion is still very similar in both models after this long
− ∗
anequalamountfrom M.Theseoperationsmadethemean
longitudes equal in thebeginning.
300
250
i=30 ; δa=1 km
3T.h1e aTppwrooxsimataetleliatensalytical theory suggest that to study wo thrusters 112050000 red=J36, blue=J0
ttahhneedrtaehlseamotriavyle,litmthoriutsisoetnenrooupfgothwheetro.saSstceeallellicitntegssfoaomnrdeottshhmeeralcvloarelrucececestnnitsersiecsaistoyyf m, J, J, t360 5 00
Y k -50
using the theory. In this section we present some results of ∆
-100
numerical experiments and display them graphically.
-150
We used for eccentricity the value e = .01 and initial
-600 -400 -200 0 200 400 600
semi-major-axis difference of 1km. The size of the radial ∆X km
oscillation is about 2ea as one derives from the two-body
Figure 5. Comparisonof relative orbits at extreme θ-values for
motionandthealong-orbit motion istwicethisinthecoor-
the case of full perturbations (red) and in the two-body model.
dinate system in which one of the satellites is in the origin Herei=30o,δa0=1km.
and the coordinate system rotates with that satellite.
6 Mikkola and Prioroc
0.022
0.02
6
0.018
e_1, e_2, |e_1-e_2| 000...000111246 4
2
0.01
0.008 m Orbital motion
k 0
Y
0.006
0 10000 20000 30000 40000 50000 60000 70000 80000 90000 100000
Time Earth
-2
Figure 6.Evolution ofindividual eccentricities andthe relative
eccentricity. Whilethe individualeccentricities ek fluctuate con-
siderably,therelativeonee=|e1−e2|remainsalmostconstant. -4
In this example the initial semi-major axis difference was one
kilometer,ek =.01andω-edifferencewassettobe=180o.
-6
10
-15 -10 -5 0 5 10 15
X km
1
Figure 8.Threesatelliteskeepanalmosttriangleconfiguration
0.1
fromthestarttotheendofthissimulationafterabout2.6years.
22θdR/ 0.01 Iisnotnhliys6sys1t%emoftthheemsiozteioonftpheerpsyesntdeimcu.lTarhteoctohloermssehaonworthbeitaidlepnltaintye
0.001 ofthesatellitesandsizeofdotsgofromsmalltobiggeraccording
totime.
0.0001
1e-05 Thisresultisthesamewhatwasobservedlongagoinquasi-
-3 -2 -1 0 1 2 3
n satellite orbits (Mikkola et al. 2006). This numerical result
isquitecomplicatedtoproveanalyticallyandweignoresuch
Figure 7. Second derivative ofe2−nh∆niwithrespect toθ (at
considerations.
θ = 0) as function of the power n. One sees that it is positive
Wecomputedthesecondderivativeof ∆n withrespect
for all non-zero powers n, thus predicting possibility for stable h i
to θ for different powers n. The result is that this quantity
motionneareachother.
is positive for any value of n, positive or negative. Thus
one may conclude that any mutualforce between the satel-
integration. Thus the effect of the potential terms higher lites that are derivable from a power-potential could result
than J2 seems minor at least in the θ angle. intostablequasi-satellitelikemotion.Infacttherearemany
Figure (3) we compare the θ motion in the full Earth othersuchforces since at least any sumof power-potentials
potentialandpuretwo-bodymotionforaboutthreemoths. with positive coefficients would do here. In practice, how-
Thereiscleardifferenceinamplitudeandalsoaphasediffer- ever, the simple constant force is probably the most useful
ence.Duetothesimilarityofthemotionsonemayconclude for real satellites. The only interesting point here is that in
that these differences are mainly due to the J2 term that thecaseof realworld quasi-satellite motion theforce isr−2
makes the energy values of the satellites different and the (plustheindirect term whichhasonly asmalleffect).Thus
initial derivative values, as computed using the two-body thedynamicswith mutualforces, derivablefrom any power
formulae, somewhat erroneous. In any case the differences potential, has similar features.
are not very large and thus the two-body motion based es- In Figures (9) – (12) the relative motions of the satel-
timates remain qualitatively correct. lites are illustrated by plotting it at every maximum of θ
The Fig.(4) compares θ behavior for a shorter time in foroneperiodaroundtheEarth.Theintegrationlastedov|er|
thetwomodelsasabovewiththeanalyticalresultsfollowing the time interval up to 105 time units in our system. This
from the Eq.(27). One can see that the analytical approx- correspondstoabout2.6yearsfornearEarthsatellites.The
imation predicts the period with a few percent accuracy. resultslooklikejustoneplottedcurve,butareactuallymore
Oneconcludesthattheanalyticalestimates,althoughbased that 1400 curves. Thus one sees the stability of the sys-
only on the Keplerian motion formulae, are useful approxi- tem.Theresultsareverysimilar inthefullpotentialmodel
mations. (markedJ ,whichincludesalsoLuni-Solarterms)and the
36
The difference in the relative motion curves in the J modelswithoneortwothrusters.(Onethrustermeaning
2
full perturbation and the two-body model are compared in thatonlyoneofthesatelliteshasathruster).IntheFig.(12).
Fig.(5). Again thedifferences remain acceptably small. the motion with respect to the osculating orbital plane is
In the theory the relative eccentricity e = e1 e2 is illustrated for the full-potential two-thruster computation.
| − |
an important quantity.In Fig.(6) we demonstrate that this Here clearly (like in fact in all the cases) this xz-motion
e
quantity is almost constant all the time although the indi- motion is very small compared with thexy-motion.
vidualeccentricitiesinthesatelliteorbitsvaryconsiderably.
Quasi-satellite dynamics in formation flight 7
300 300 300
250 i=10,30,50,70,90; δa=0 250 i=10,30,50,70,90; δa=0.5 km 250 i=10,30,50,70,90; δa=1.0 km
∆Y km, J, two thrusters36 112- 55050 000000 ∆Y km, J, two thrusters36 112- 50505 000000 ∆Y km, J, two thrusters36 112- 50505 000000
-100 -100 -100
-150 -150 -150
-600 -400 -200 0 200 400 600 -600 -400 -200 0 200 400 600 -600 -400 -200 0 200 400 600
∆X km ∆X km ∆X km
Figure 9.Potential uptoJ36,xy-motion,twothrusters.Thenumberofcurvesislargerthan1400.
300 300 300
250 i=10,30,50,70,90; δa=0 250 i=10,30,50,70,90; δa=0.5 km 250 i=10,30,50,70,90; δa=1.0 km
∆Y km, J, two thrusters2 112- 55050 000000 ∆Y km, J, two thrusters2 112- 50505 000000 ∆Y km, J, two thrusters2 112- 55050 000000
-100 -100 -100
-150 -150 -150
-600 -400 -200 0 200 400 600 -600 -400 -200 0 200 400 600 -600 -400 -200 0 200 400 600
∆X km ∆X km ∆X km
Figure 10.Potential onlytoJ2,xy-motion,twothrusters.Thenumberofcurvesislargerthan1400.
300 300 300
250 i=10,30,50,70,90; δa=0 250 i=10,30,50,70,90; δa=0.5 km 250 i=10,30,50,70,90; δa=1.0 km
∆Y km, J, one thruster36 112- 55050 000000 ∆Y km, J, one thruster36 112- 50505 000000 ∆Y km, J, one thruster36 112- 50505 000000
-100 -100 -100
-150 -150 -150
-600 -400 -200 0 200 400 600 -600 -400 -200 0 200 400 600 -600 -400 -200 0 200 400 600
∆X km ∆X km ∆X km
Figure 11.Potential uptoJ36,xy-motion,onethruster.Thenumberofcurvesislargerthan1400.
300 300 300
250 i=10,30,50,70,90; δa=0 250 i=10,30,50,70,90; δa=0.5 km 250 i=10,30,50,70,90; δa=1.0 km
∆Z km, J, two thrusters36 112- 55050 000000 ∆Z km, J, two thrusters36 112- 50505 000000 ∆Z km, J, two thrusters36 112- 50505 000000
-100 -100 -100
-150 -150 -150
-600 -400 -200 0 200 400 600 -600 -400 -200 0 200 400 600 -600 -400 -200 0 200 400 600
∆X km ∆X km ∆X km
Figure 12.Potential uptoJ36,xz-motion,twothrusters.Thenumberofcurvesislargerthan1400.
8 Mikkola and Prioroc
3.2 More satellites aresurprisinglysmall.Inanycase,therequiredacceleration
is easy to find by numerical simulations. The experiments
For three satellite experiments we used the initial values as
withthehighorderexpansions,whichcontainstherelevant
explained in the beginning of this section. The thruster ac-
perturbationsfor thesatellite motion, leads to same results
celerations for all the satellites were directed according to
astheoneswithonlytheJ term,atleast withinthepreci-
the mean direction of the two other satellites. In some ex- 2
sionofourillustrations.Thisresultsuggeststhattherelative
periments we randomized these by modifying the direction
motionundertheactionofthethrusterforceisnotsensitive
cosines randomly by 10% (and then we renormalized the
to perturbations. As mentioned in section 3.2, even some-
directionvector).Theseoperationsdidnothaveanynotice-
what randomized directions of the thruster force does not
able effect.
havean effect in therelative motion of thesatellites.
The result was that the satellites move in a kind of
Oneadvantageofthismethodisthatthesatellitesneed
triangle configuration. Actually the orbits with respect to
only information about thedirections of other satellite(s).
thecenterofthetripleareapproximatelyellipses,butifthe
Acknowledgments: Claudiu Prioroc’s research has
scaleofthevariationofthedistancefromtheEarthisscaled
been funded by the European Commission through the As-
by a factor of 2 theconfiguration looks like a circle. This is
trodynamics Network AstroNet-II, under Marie Curie con-
illustrated in Fig. 8 where the smaller dots represent the
tract PITN-GA-2011-289240.
situation in the beginning and the larger ones show it after
about 2.6 years.
Verysimilarresultswereobtainedalsowithmoresatel-
REFERENCES
lites, for example systems of six satellites seems to behave
in a similar way. Thus it seems possible to keep a cluster Alfriend, K. T., Vadali, S. R., Schaub, H. 2001. Celestial
of satellites close to each other just by using constant re- Mechanics and Dynamical Astronomy 81, 57-62.
pulsiveforces. Thisis probably simplerand otherthinkable Beichman,C.,G´omez,G.,Lo,M.,Masdemont,J.,Romans,
ways for formation flight of satellites. However this method L.2004. Advancesin Space Research 34, 637-644.
justprovidesawaytokeepsatellitesinakindofellipticlike Clohessy, W. H. and Wiltshire, R. S., Journal of the
configuration. AerospaceSciences,Vol.27,No.5,1960,pp.653658, 674.
HillG.W.,AmericanJournalofMathematics.1,1878,5-26.
Gim, Dong-Woo and Alfriend, K.T. 2005 Celestial Me-
chanics and Dynamical Astronomy 92, 295336. ghEcc
4 DISCUSSION AND CONCLUSION Jackson, J. 1913, MNRAS,74, 62-82.
Kristiansen,K.U.,Palmer,P.L.,Roberts,M.2010.Celes-
The motions of the satellites in low Earth orbits are essen-
tial Mechanics and Dynamical Astronomy 106, 371-390.
tially just elliptic motions so that the relative motions are
Prioroc,C-L.,Mikkola,S.2015 NewAstronomy34,41-46.
due to differences in two-body motion. Our numerical ex-
Mikkola, S., Innanen, K., Wiegert, P., Connors, M.,
perimentsdemonstratethatthetandemflighthastheeffect
Brasser, R. 2006. Monthly Notices of the Royal Astro-
thattheperturbationsfromtheEarthharmonics,whichare
nomical Society 369, 15-24.
orders of magnitude larger than the needed thruster force,
Mikkola, S., Palmer, P., Hashida, Y. 2002. Celestial Me-
affect the motions of the satellites in the same way and do
chanics and Dynamical Astronomy 82, 391-411.
notmakemuchdifferenceintheirrelativemotion.Thiswas
Roscoe, C.W.T., Vadali, S.R., Alfriend, K.T., Desai, U.P.
alsoconfirmedtobethecasewhenhighorderexpansionfor
2013 ActaAstronautica 82, 1624.
theEarthgravitational potential(equation(11))isused.In
Schaub, H., Alfriend, K. T. 2001. Celestial Mechanics and
otherwordstherelativemotionofsatellitesismainlydueto
Dynamical Astronomy 79, 77-95.
eccentricities so that in the orbital coordinate system they
Scheeres, D.J. 1998. Celestial Mechanics and Dynamical
seem to move approximately along a small ellipse in which
Astronomy 70: 7598.
theaxis ratio is 1:2. The only thingthe thrustersmust do
Vallado, D. A., Alfano, S. 2014. Celestial Mechanics and
istokeepthemeanlongitudesnearlysame.Forthispurpose
Dynamical Astronomy 118, 253-271.
very small thrusts are enough when the orbits of the satel-
lites are close to each other. Little differences in the orbits
are, however, not important as we found with simulations.
Themethodwoksaswellfortwoorthreeorevenmore
satellites.
Surprisingly even near triangle configuration (in case
of three satellites) is stable in the sense that the relative
motionofthesatelliteskeepthesameconfigurationforvery
long times.
In our experiments, with the constant acceleration
model, the value of relative acceleration of order ǫ =
10−6g [g=the standard gravity value ( 10m/sec2)] is
∼
enough. Writing this in physical units we can say that
ǫ 10−5m/sec2 =.01mm/sec2 is enough, but even smaller
∼
values may do. Only if the initial difference in semi-major-
axis is more than a kilometer or so, a somewhat stronger
force is required, but even in such cases the forces needed