- Email: [email protected]
A linear model for relative motion in an elliptical orbit based on a spherical coordinate system
Accepted Manuscript A linear model for relative motion in an elliptical orbit based on a spherical coordinate system Chao Han, Huan Chen, Gustavo Alonso, Yinrui Rao, Javier Cubas, Jianfeng Yin, Xiaohui Wang PII:
S0094-5765(18)30910-X
DOI:
https://doi.org/10.1016/j.actaastro.2019.01.016
Reference:
AA 7282
To appear in:
Acta Astronautica
Received Date: 1 June 2018 Revised Date:
21 November 2018
Accepted Date: 10 January 2019
Please cite this article as: C. Han, H. Chen, G. Alonso, Y. Rao, J. Cubas, J. Yin, X. Wang, A linear model for relative motion in an elliptical orbit based on a spherical coordinate system, Acta Astronautica (2019), doi: https://doi.org/10.1016/j.actaastro.2019.01.016. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
RI PT
A Linear Model for Relative Motion in an Elliptical Orbit Based On a Spherical Coordinate System Chao Hana , Huan Chena , Gustavo Alonsob , Yinrui Raoc , Javier Cubasb , Jianfeng Yind , Xiaohui Wanga,∗ a School
M AN U
SC
of Astronautics, Beihang University, No.37, Xueyuan Road, Haidian District, Beijing, 100191, China b Faculty of Aerospace Engineering, Technical University of Madrid, Madrid, 28040, Spain c Institute of Systems Engineering, China Academy of Engineering Physics, Mianyang, Sichuan Province, 621999, China d Beijing Institute of Spacecraft System Engineering, China Academy of Space Technology, Beijing, 100094, China
Abstract
When compared with Cartesian coordinates, curvilinear coordinates have shown significant superior accuracy when applied to a linear model for relative motion in a circular reference orbit. Experts assume the same superiority for an elliptical reference orbit. However, the curvilinear model in an elliptical reference
TE D
orbit has not been established, and the superiority has not been proved. After a strictly theoretical derivation this paper establishes a linear model with spherical coordinates for relative motion in an elliptical reference orbit. The differential equation has the same form as the Lawden or Tschauner-Hempel
EP
equation and can be easily solved with the solutions from the literature. In addition both the theoretical derivation and the numerical examples prove the model’s superior accuracy in leader-follower formation and formations with a
AC C
large along-track distance. The proposed model reduces the approximate limits of close relative motion, and it is expected to be applied in the propagation, guidance, and control for longer term relative motion with large along-track distances.
∗ Corresponding
author Email addresses: [email protected] (Chao Han), [email protected] (Huan Chen), [email protected] (Gustavo Alonso), [email protected] (Yinrui Rao), [email protected] (Javier Cubas), [email protected] (Jianfeng Yin), [email protected] (Xiaohui Wang)
Preprint submitted to Journal of LATEX Templates
January 12, 2019
ACCEPTED MANUSCRIPT
Keywords: formation flying, relative motion, linear model, spherical
RI PT
coordinate 2010 MSC: 00-01, 99-00
1. Introduction
Relative motion dynamics is a key space mission subject area, highly im-
SC
portant to efforts related to formation flying, spacecraft rendezvous, terminal
guidance, space debris evolution and collision avoidance. For the relative motion 5
in a central gravity field Joshua Sullivan [1] provides a solid review of several
M AN U
representative models, including the Hill [2] or Clohessy-Wiltshire[3] (HCW) equations for the circular reference orbit and Lawden [4] or Tschauner-Hempel [5] (LTH) equations for the elliptical reference orbit. Both the HCW and LTH equations describe the relative motion in Cartesian coordinates and they are lin10
earized based on the assumption of close relative motion. For HCW equations in the circular case Alfriend [6] analyzed the effects of linearizing the initial condi-
TE D
tion mapping between the inertial state and relative state, and he found that the linearized mapping would cause considerable radial error when the along-track distance is large, which will limit the models’s application. 15
Many researchers have devoted effort to solving the above problem. A feasible approach is to establish the relative motion model in curvilinear coordinates,
EP
such as cylindrical coordinates or spherical coordinates. So far as we know, Berreen, T. and Crisp, J [7] were the first to set up a first order relative mo-
AC C
tion model based on cylindrical coordinates for a circular reference orbit. They 20
applied the model to the mission of releasing a detector from the space station, which is on a circular orbit. Geller, D. K. [8] proved that using cylindrical coordinates can remove the nonlinear terms related to the along-track distance and cross-track velocity. The resulting linear model in the HCW form is valid for arbitrarily large along-track distances and cross-track velocities. Alfriend [6] de-
25
rived a first order model, which has the same form as the HCW equations, based on the spherical coordinate system for a circular reference orbit, . De Bruijin,
2
ACCEPTED MANUSCRIPT
F. [9] compared the accuracy of the first order relative motion models based on
RI PT
both curvilinear coordinates and Cartesian coordinates that were in the form of the HCW equation. He found that for a formation with a large along-track 30
distance, the model based on curvilinear coordinates showed superior accuracy and reduced the fuel cost of formation maintenance maneuvers. However, for other formations with small along-track distances, the model based on curvilin-
SC
ear coordinates had no superiority or performed even worse. Further research, also based on curvilinear coordinates, achieved better accuracy by obtaining the 35
analytical solution to the second order model [10] and the approximate analyt-
M AN U
ical solution to a highly nonlinear model [11], or by taking into account the J2 perturbation [12]. When researching relative motion in elliptical reference orbits, Yamanaka [13] claimed that using cylindrical coordinates can significantly lessen the proximity requirement, although this work actually has not been con40
ducted or examined. Because formation missions have recently been extended to elliptical orbits recently [14, 15, 16], it is of great importance to investigate
reference orbit.
TE D
the application of curvilinear coordinates to the relative motion in an elliptical
This paper has two major innovations. First, the linearized equation for 45
relative motion in an elliptical reference orbit is obtained with a strictly dynamic derivation based on spherical coordinates. Second, the theoretical derivation
EP
and numerical examples show that the model based on spherical coordinates has superior accuracy in leader-follower formation when compared with Cartesian coordinates. The manuscript is organized as follows. The first section is the introduction. The second section is the theoretical part. First, we present the
AC C
50
definition of spherical coordinates, and we present the mappings between the coordinates and the typical inertial orbital parameters, using both the forms of position-velocity vectors and the classical orbital elements. Then the differential equations for relative motion in an elliptical reference orbit are derived, whose
55
analytical solutions can be found in the literature [13, 17, 18, 19, 20] because they are similar in form to the LTH equation. In addition the solutions’ accuracy are investigated, and the theoretical derivation shows its superior accuracy in 3
ACCEPTED MANUSCRIPT
leader-follower formation. In the third section, numerical examples examine the
60
RI PT
model’s accuracy in comparison with the Cartesian model. The last section comes to a conclusion.
2. Relative Motion Model 2.1. Definition of Spherical Coordinates
SC
The reference satellite A is denoted as the chief, and another satellite B is
the deputy. In the following sections the subscript ∗ (∗ = a, b) indicates the chief or the deputy. In Figure 1 two satellites are currently located at points A and
M AN U
65
B, respectively. The ascending nodes of the two orbits are Ma and Mb . The orbit planes of A and B intersect at point J, and the interior angle between the orbit planes is ∆j. The angle separation between point J and the ascending node M∗ is γ∗ . The angle separation between point J and the satellite is β∗ . 70
Suppose the projection of the radius vector of satellite B onto the orbit plane
TE D
of satellite A intersects with satellite A’s orbit at point B 0 . Orbit Plane of Satellite B
B
zi
rˆb
ha
EP
hb
AC C
O
r
b a
A
r
rˆbi
J
rˆa
a
Orbit Plane of Satellite A
B
j
b
ib
ia
Equator
Mb
Ma
Figure 1: Definition of Spherica Coordinate
The spherical coordinates are defined similarly to those in reference [6]. The parameter ρ denotes the radius difference between the chief and deputy. Denote 4
ACCEPTED MANUSCRIPT
the angle difference between OA and OB 0 as θr , and the angle difference between
RI PT
OB and OB 0 as φr . The coordinate transformation between the Earth center orbital (ECO) coordinates of the chief and the deputy can be handled as follows: Let the chief’s ECO frame Soa rotates by θr around the z axis, and then rotate by −φr around the y axis.
SC
Soa −−→ [Lz (θr )] −−→ [Ly (-φr )] Sob
In this way, ρ, θr , and φr indicate the deputy’s relative motion in the radial direction, along-track, and normal direction of the chief’s orbit, respectively.
75
M AN U
The deputy’s relative state parameters with respect to the chief are ρ, θr , φr , ρ, ˙ θ˙r and φ˙ r . Here the dot on the top indicates the derivative with respect to time.
zo (z'o)
y'o
O z r r xy
TE D
A
yo
B
xo
B' x'o
Figure 2: Geometric Description of ECO Coordinate
EP
The corresponding Cartesian coordinates in the local vertical local horizontal
AC C
(LVLH) coordinate system can be expressed by Equation (1): x = (ra + ρ) cos φr cos θr − ra y = (ra + ρ) cos φr sin θr z = (r + ρ) sin φ a r
(1)
And the relative velocity can be obtained from Equation (1) as given by Equa-
5
ACCEPTED MANUSCRIPT
RI PT
tion (2): ˙ r sin φr cos θr + θ˙r cos φr sin θr − r˙a x ˙ = ( r ˙ + ρ) ˙ cos φ cos θ − (r + ρ) φ a r r a y˙ = (r˙a + ρ) ˙ cos φr sin θr − (ra + ρ) φ˙ r sin φr sin θr − θ˙r cos φr cos θr z˙ = (r˙a + ρ) ˙ sin φr + (ra + ρ) φ˙ r cos φr (2) Equations (1) and (2) will be used for accuracy comparison of the spherical
SC
and Cartesian models.
2.2. Mapping between the Inertial Radius-Velocity Vectors and Spherical Coordinates
M AN U
80
Here denote a vector variable with bold type as x = xˆ x, where x = kxk is ˆ with a hat is the unit vector along x. The inertial radius its square norm and x and velocity vectors are denoted as r∗ and v∗ (∗ = a, b).
2.2.1. Mapping from the inertial radius-velocity vectors to spherical coordinates 85
When the chief and deputy’s inertial radius and velocity vectors, which are
TE D
ra , va , rb and vb , are known, the deputy’s relative state parameters with respect to the chief, i.e. ρ, θr , φr , ρ, ˙ θ˙r and φ˙ r , remain to be calculated. The velocity projection in the radial direction is given by Equation (3): v∗r = v∗ · ˆr∗
(3)
EP
And the radial gradient is given by Equation (4): r˙∗ = v∗r
(4)
AC C
Then ρ and ρ˙ can be obtained from Equations (5) and (6): ρ = rb − ra
(5)
ρ˙ = vbr − var
(6)
The angular momentum vector of chief satellite is given by Equation (7): ha = ra × va 6
(7)
ACCEPTED MANUSCRIPT
ˆ a of chief’s The projections of rb and vb onto the perpendicular direction h
RI PT
orbit plane are given by Equations (8) and (9): ˆa h ˆa rbh = rb · h
(8)
ˆa h ˆa vbh = vb · h
(9)
(10) and (11):
M AN U
rbi = rb − rbh
SC
The projections of rb and vb in the chief’s orbit plane are given by Equations
vbi = vb − vbh
(10)
(11)
Then θr can be obtained from Equation (12):
θr = sign (θr ) · arccos (ˆra · ˆrbi )
(12)
where sign (θr ) is the sign of θr , which is positive when ra rotates clockwise to
TE D
ˆ a , and is determined by Equation (13): rbi around h sign (θr ) = sign
ˆa (ra × rbi ) · h
(13)
The angular velocities of vectors ra and rbi are given by Equations (14) and
AC C
EP
(15):
kra × va k θ˙a = ra2
(14)
krbi × vbi k θ˙bi = 2 rbi
(15)
Then θ˙r can be obtained from Equation (16) θ˙r = θ˙bi − θ˙a
(16)
rbh rb
(17)
From Equation (17) sin φr =
7
ACCEPTED MANUSCRIPT
and its derivative, Equation (18),
vbh − vbr sin φr φ˙ r = rb cos φr
SC
we obtain φr and φ˙ r from Equations (19) and (20) rbh φr = arcsin rb
(18)
RI PT
vbh rb − rbh vbr φ˙ r cos φr = rb2
(19)
(20)
M AN U
2.2.2. Mapping from spherical coordinates to the inertial radius-velocity vectors When the chief’s inertial radius and velocity vectors and the deputy’s relative 90
state, given by ra , va , ρ, θr , φr , ρ, ˙ θ˙r and φ˙ r , are known, the deputy’s inertial radius and velocity vector, i.e. rb and vb , remain to be calculated. First, from Equations (5) and (6), we have Equations (21) and (22):
TE D
rb = ra + ρ
vbr = var + ρ˙
(21)
(22)
Then from Equations (19) and (20), we have Equations (23) and (24): (23)
vbh = φ˙ r rb cos φr + vbr sin φr
(24)
AC C
EP
rbh = rb sinφr
Therefore, we produce Equation (25): rbi =
q
2 rb2 − rbh
(25)
And the unit vector along rbi can be obtained from Equation (26): ˆrbi = Lz (θr ) ˆra
8
(26)
ACCEPTED MANUSCRIPT
where Lz (θr ) is the coordinate transition matrix for rotation along the z axis
cos θr Lz (θr ) = sin θr 0
0 0 1
− sin θr cos θr 0
ˆ a + rbi ˆrbi rb = rbh h
SC
Consequently, we have Equation (28):
RI PT
as shown by Equation (27):
(27)
(28)
M AN U
And using Equation (16), we get Equations (29) and (30): θ˙bi = θ˙a + θ˙r
(29)
vbt = rbi θ˙bi
(30)
Finally, we obtain Equation (31):
(31)
TE D
vb = [vbr , vbt , vbh ]T
2.3. Mapping between the Classical Orbital Elements and Spherical Coordinates 2.3.1. Mapping from the classical orbital elements to spherical coordinates When the chief and deputy’s orbital elements, which are the semi-major axis a∗ , eccentricity e∗ , inclination i∗ , right ascension of ascending node Ω∗ ,
EP
95
argument of perigee ω∗ and true anomaly θ∗ , are known, the deputy’s relative
AC C
state parameters with respect to the chief, i.e. ρ, θr , φr , ρ, ˙ θ˙r and φ˙ r , remain to be calculated.
For the Keplerian motion, we use Equation (32): r∗ =
p∗ 1 + e∗ cos θ∗
(32)
where p is the semi-latus rectum given by Equation (33): p∗ = a∗ 1 − e2∗
9
(33)
ACCEPTED MANUSCRIPT
(34)
RI PT
The radial velocity v∗r is given by Equation (34): r µ v∗r = r˙∗ = e∗ sin θ∗ p∗
Then ρ and ρ˙ can be obtained from Equations (5) and (6).
In addition the transverse velocity v∗u is given by Equation (35): h∗ v∗u = r∗ θ˙∗ = r∗
SC
(35)
where h∗ is the module of the angular momentum as given by Equation (36): √
µp∗
M AN U
h∗ =
(36)
Consequently, θ˙∗ can be obtained from Equation (37): h∗ θ˙∗ = 2 r∗
(37)
Spherical geometry gives Equations (38) and (39):
TE D
cos ∆j = − cos ia cos (π − ib ) + sin ia sin (π − ib ) cos ∆Ω
sin γa sin γb sin ∆Ω = = sin ∆j sin (π − ib ) sin ia
(38)
(39)
The parameters ∆j, γa and γb can be obtained from Equations (40), (41),
EP
and (42):
AC C
∆j = sign (∆i) · arccos (cos ia cos ib + sin ia sin ib cos ∆Ω)
100
γa = arcsin γb = arcsin
sin ∆Ω sin ib sin ∆j
sin ∆Ω sin ia sin ∆j
(40)
(41)
(42)
For Equations (41) and (42), when ∆j = 0, that is when ∆i = ∆Ω = 0, −−→ −−→ vectors OB and OB 0 coincide, and thus γa = γb = 0.
10
It’s obvious that Equation (43) holds: βa = ωa + θa − γa
(43)
βb = ωb + θb − γb
Spherical geometry gives Equations (44) and (45): tan βb0 = tan βb cos ∆j
(44)
SC
sin φr = sin βb sin ∆j
RI PT
ACCEPTED MANUSCRIPT
(45)
M AN U
The parameters βb0 , θr and φr can be obtained from Equations (46), (47), and (48):
(46)
θr = βb0 − βa
(47)
φr = arcsin (sin βb sin ∆j)
(48)
TE D
βb0 = arctan (tan βb cos ∆j)
Using Equations (46), (47), and (48), we obtain the angles’ derivative as given by Equations (49), (50), and (51): (49)
θ˙r = β˙ b0 − θ˙a
(50)
cos βb φ˙ r = θ˙b sin ∆j cos φr
(51)
AC C
EP
2 cos2 βb0 ˙b sin βb0 1 cos ∆j = θ β˙ b0 = θ˙b cos2 βb sin2 βb cos ∆j
Here β˙ b0 is presented in two forms to avoid the singularity when βb =
N.
11
n 2 π,
n∈
ACCEPTED MANUSCRIPT
2.3.2. Mapping from spherical coordinates to the classical orbital elements When the chief’s orbital elements, i.e. aa , ea , ia , Ωa , ωa and θa , and the
RI PT
105
deputy’s relative state parameters, i.e. ρ, θr , φr , ρ, ˙ θ˙r and φ˙ r , are known, the deputy’s orbital elements, i.e. ab , eb , ib , Ωb , ωb and θb , are to be calculated. Firstly from Equation (50), we get Equation (52): β˙ b0 = θ˙a + θ˙r Spherical geometry gives Equation (53):
M AN U
tan φr = sin βb0 tan ∆j
SC
(52)
(53)
and its derivative, given by Equation (54):
φ˙ r = β˙ b0 cos βb0 tan ∆j cos2 φr
TE D
we obtain βb0 and ∆j with Equations (55) and (56), ! β˙ b0 sin 2φr βb0 = arctan 2φ˙ r v u u ∆j = arctan ttan2 φr +
φ˙ r β˙ b0 cos2 φr
(54)
(55)
!2 (56)
EP
Spherical geometry gives Equation (57): sin βb =
sin φr sin ∆j
AC C
Then βb can be obtained as Equation (58): sin φr βb = arcsin sin ∆j
(57)
(58)
Second, from Equations (47) and (43), we obtain Equations (59) and (60): βa = βb0 − θr
(59)
γa = ωa + θa − βa
(60)
12
ACCEPTED MANUSCRIPT
Then spherical geometry gives Equations (61) and (62): (61)
sin ∆Ω sin γa sin γb = = sin ∆j sin ib sin ia
(62)
RI PT
cos (π − ib ) = − cos ia cos ∆j + sin ia sin ∆j cos γa
We obtain ib , ∆Ω and γb using Equations (63), (64), and (65):
sin γa sin ∆j sin ib
γb = arcsin
sin γa sin ia sin ib
(63)
M AN U
∆Ω = arcsin
SC
ib = arccos (cos ia cos ∆j − sin ia sin ∆j cos γa )
(64)
(65)
Therefore, we get Equation (66):
ub = ωb + θb = γb + βb
(66)
TE D
Taking the derivative of Equations (66) and (57) provides Equation (67): θ˙b = β˙ b =
φ˙ r cos φr cos βb sin ∆j
(67)
Given Equation (56), when φr = 0 and φ˙ r = 0, ∆j = 0, and we find
EP
110
Equations (58) and (67) are singular. In this case, ∆i = ∆Ω = 0, and the −−→ −−→ vector OB and OB 0 coincide. Then Equations (58) and (67) can be replaced with βb = ωa + θa − θr and θ˙b = β˙ b = θ˙a − θ˙r .
AC C
From Equations (5) and (6), we get Equations (68) and (69): rb = ra + ρ
(68)
r˙b = r˙a + ρ˙
(69)
Substitute Equations (67) and (68) into Equations (37) and (36) to get
Equation (70): pb =
rb4 θ˙b2 µ
13
(70)
ACCEPTED MANUSCRIPT
Then substitute Equations (68), (69), and (70) into Equations (32) and (34)
s eb =
2
(pb − rb ) pb + r˙b2 2 rb µ
tan θb =
rb r˙b pb − rb
r
pb µ
RI PT
to get Equations (71) and (72): (71)
(72)
SC
Finally, ab and ωb can be obtained with Equations (33) and (66). 2.4. Relative Motion Equations
M AN U
ˆo, y ˆo The unit vectors of the chief’s ECO coordinate axes are denoted as x ˆo . Then the chief’s unit radial vector, angular velocity and acceleration and z
TE D
vectors can be expressed by Equation (73): ˆra = x ˆo ˆo wa = θ˙a z w ˙ a = θ¨a z ˆo
(73)
Using the coordinate transformation, the deputy’s unit radial vector can be
(74)
EP
expressed in the chief’s ECO coordinates as Equation (74): 1 cos φr cos θr ˆrb = Lz (−θr ) Ly (φr ) 0 = cos φr sin θr 0 sin φr
AC C
The deputy’s angular velocity is given by Equation (75): −φ˙ r sin θr wb = wa + wb/a = φ˙ r cos θr θ˙a + θ˙r
where
wb/a
0 0 −φ˙ r sin θr = 0 + Lz (−θr ) φ˙ r = φ˙ r cos θr θ˙r 0 θ˙r 14
(75)
(76)
ACCEPTED MANUSCRIPT
Then the deputy’s angular acceleration is given by Equation (77):
RI PT
˙b=w ˙a+w ˙ b/a = w ˙ a|a + wa × wa + w ˙ b/a|a + wa × wb/a w ˆo = − φ¨r sin θr + θ˙a + θ˙r φ˙ r cos θr x ˆ o + θ¨a + θ¨r z ˆo + φ¨r cos θr − θ˙a + θ˙r φ˙ r sin θr y
(77)
The chief’s radial vector and its first and second-order derivatives are given
rb = rb ˆrb
M AN U
r˙ b = r˙b ˆrb + rb wb × ˆrb
SC
by Equations (78), (79), and (80):
¨rb = r¨b ˆrb + 2r˙b wb × ˆrb + rb w ˙ b × ˆrb + rb wb × (wb × ˆrb )
(78)
(79)
(80)
Under the condition of close relative motion, Equation (80) can be linearized
TE D
to form Equation (81): r¨a + ρ¨ − 2r˙a θ˙a θr − ra θ˙a2 + 2θ˙a θ˙r + θ¨a θr − θ˙a2 ρ ¨rb = r¨a θr + 2 r˙a θ˙a + r˙a θ˙r + ρ˙ θ˙a + ra −θ˙a2 θr + θ¨a + θ¨r + θ¨a ρ r¨a φr + 2r˙a φ˙ r + ra φ¨r
(81)
The forces on the deputy include the central gravity gb and the control force
EP
fb as given by Equation (82):
¨rb = gb + fb
(82)
AC C
where the central gravity is given by Equation (83): cos φr cos θr µ µ gb = − 2 ˆrb = − cos φ sin θ r r 2 rb (ra + ρ) sin φr and gb can be linearized as Equation (84): 1 µ µ gb = − 3 (ra − 2ρ) θr = − 3 ra ra φr 15
ra − 2ρ ra θr ra φr
(83)
(84)
ACCEPTED MANUSCRIPT
For the chief with only central gravity we have Equation (85): ¨ra − ga = 0
RI PT
(85)
That is spelled out as Equation (86): µ r¨a − ra θ˙a2 + 2 = 0 ra 2r˙ θ˙ + r θ¨ = 0 a a
SC
a a
(86)
Because
d dt
d = θ˙ dθ ,
d2 dt2
M AN U
Subtracting Equation (82) from Equation (86) we get Equation (87): ρ¨ − 2θ˙a r˙a θr + ra θ˙r − θ˙a2 ρ − θ¨a ra θr −2ρ µ ¨ ra θr + 2r˙a θ˙r + r¨a θr + 2θ˙a ρ˙ − θ˙a2 ra θr + θ¨a ρ + 3 ra θr = fbc ra ra φr ra φ¨r + 2r˙a φ˙ r + r¨a φr d2 ˙ dθ˙ d = θ˙2 dθ 2 +θ dθ dθ , we can denote
00
d(∗) dθ
0
= (∗) ,
d2 (∗) dθ 2
=
(88)
TE D
(∗) , to give Equation (88): µ2 θ˙ = 3 (1 + e cos θ)2 h 2 0 θ˙ = − 2µ e (1 + e cos θ) sin θ 3 h
(87)
EP
Then Equation (87) can be rearranged to form Equation (89): 0 ((1 + ea cos θa ) ρ00 − 2ea sin θa ρ0 ) − (3 + ea cos θa ) ρ − 2 (1 + ea cos θa ) (ra θr ) + 2ea sin θa (ra θr ) µ 00 0 0 (1 + ea cos θa ) (ra θr ) − 2ea sin θa (ra θr ) − ea cos θa (ra θr ) + 2 (1 + ea cos θa ) ρ − 2ea sin θa ρ = fbc ra3 00 0 (1 + ea cos θa ) (ra φr ) − 2ea sin θa (ra φr ) + (ra φr ) (89)
AC C
In addition, we know that Equation (90) holds: ((1 + e cos θ) x)0 = (1 + e cos θ) x0 − e sin θx
00
(90)
((1 + e cos θ) x) = (1 + e cos θ) x00 − 2e sin θx0 − e cos θx
Substitute Equation (90) into (89) to get Equation (91): 00 0 ((1 + ea cos θa ) ρ) − 3ρ − 2((1 + ea cos θa ) ra θr ) µ 00 0 = fbc ((1 + e cos θ ) r θ ) + 2((1 + e cos θ ) ρ) a a a r a a 3 ra 00 ((1 + ea cos θa ) ra φr ) + (1 + ea cos θa ) ra φr 16
(91)
ACCEPTED MANUSCRIPT
Rearrange the above equation to get Equation (92):
µ h2a ra3 µ
ρ ra
00
3ρ ra (1+ea cos θa ) 0 θr00 + 2 rρa
−
φ00r + φr
− 2θr0
ρ ra
00
= ra θ˙a2
3ρ ra (1+ea cos θa ) 0 θr00 + 2 rρa
−
− 2θr0
RI PT
= fbc
φ00r + φr
(92)
In the following sections, the subscript ’a’ referring to the chief satellite is
M AN U
SC
omitted for simplicity. Denote ρ¯ = ρr and k = 1 + e cos θ to give Equation (93): 00 0 ρ¯ − 3¯ ρ/k − 2θr 2 ˙ 00 (93) rθ = fbc θr + 2¯ ρ0 φ00r + φr For the relative motion without the control force, i.e. fbc = 0, we have Equation (94):
θr00 + 2¯ ρ0 = 0
(94)
φ00r + φr = 0
TE D
ρ¯00 − 3¯ ρ/k − 2θr0 = 0
For the circular reference orbit, e = 0, θ˙ = n, θ¨ = 0, rb = a + ρ,
µ r3
= n2 ,
EP
where n is the mean motion. Equations (94) can be simplified to Equation (95): ρ¨ − 3n2 ρ − 2naθ˙r = 0 θ¨r + 2nρ/a ˙ =0 φ¨ + n2 φ = 0 r
(95)
AC C
Equations (94) and (95) have the same form as the LTH and HCW equations, 115
respectively. And the result of the circular case is the same as Equation (5.88) in reference [6], while ignoring the typo that the sign of 2n0 ρ˙ in the second equation of (5.88) should be positive. 2.5. Solutions to the Relative Motion The relative motion equations have the same form as the LTH equations, and many solutions have been proposed in previous research [13, 17, 18, 19, 20].
17
ACCEPTED MANUSCRIPT
is adopted in this paper. ρ¯ = c1 sθ + c2 cθ + c3 (2 − 3esθ L) θr = c1 qcθ − c2 qsθ − 2c3 k 2 L + c4 φr = c5 cos θ + c6 sin θ ρ¯0 = c1 s0θ + c2 c0θ − 3c3 e s0θ L + sθ /k 2
(96)
SC
RI PT
Yamanaka’s solution [13] is shown as Equation (96). His solution is concise and
θr0 = −2c1 sθ − c2 (2cθ − e) + 3c3 (2esθ L − 1) φ0r = −c5 sin θ + c6 cos θ
M AN U
where ci are the constant parameters determined by the initial state, and for the rest we have Equations (97) and (98): L=
(97)
q = 1 + 1/k sθ = k sin θ
TE D
µ2 (t − t0 ) h3
cθ = k cos θ
(98)
s0θ = cos θ + e cos 2θ
c0θ = − (sin θ + e sin 2θ)
EP
¯ = [¯ Denote the state parameters as x ρ, θr , φr , ρ¯0 , θr0 , φ0r ]T and the constant parameters in the solution as C = [c1 , c2 , c3 , c4 , c5 , c6 ]T , then the relative motion
AC C
solutions can be expressed as Equation (99) [13, 19, 20]: x ¯ = Φθ C
(99)
The state transition equation is given by Equation (100): ¯ 0 = Φθθ0 x x ¯ = Φθ Φ−1 ¯0 θ0 x
(100)
where the parameters Φθ , Φθ0 , and Φ−1 θ0 are given by Equations (101), (102),
18
ACCEPTED MANUSCRIPT
and (103): 2 − 3eLsθ
0
0
−qsθ
−3k 2 L
1
0
0
cos θ
0
0
0 c0θ
+ sθ /k
3 (2eLsθ − 1)
0
0
0
0
0
− sin θ
s θ qcθ 0 = 0 sθ −2sθ 0
cθ
2
0
−qsθ
0
1
0
0
c0θ
−3esθ /k
e − 2cθ
−3
0
0
2
(101)
θ
0
0
0
0
cos θ
0
0
0
0
0 − sin θ
0 sin θ 0 0 cos θ
0
0
cθ − 2e
−qsθ
0
0
−sθ
−e − qcθ
0
0
esθ
k2
w
0
ecθ − 2
−eqsθ
0
w cos θ
0
0
0
w sin θ
0
0
TE D
−3 1 + e2 /k sin θ −3 (e + cos θ) 3k + e2 − 1 1 = w −3eq sin θ 0 0
0 sin θ 0 0 cos θ
(102)
θ0
0
0 0 −w sin θ w cos θ 0
EP
Φ−1 θ0
−3e
2
e − 2cθ
Φθ0
0 Ls0θ
0
RI PT
qcθ 0 Φθ = 0 sθ −2sθ 0
cθ
SC
sθ
M AN U
θ0
(103)
AC C
where w = 1 − e2 . 120
2.6. Analysis of Solution Accuracy A closed relative orbit with no long-term relative drift is called a periodic
solution. In order to evaluate the accuracy of the proposed model, this section investigates the real relative drift of the periodic solution to the model. Following the same procedure as in reference [21], the relative motion drift
125
per orbit can be obtained with Equations (96) throuth (103) to give Equations
19
ACCEPTED MANUSCRIPT
(104), (105), and (106):
θr|drift = −
RI PT
6πc3 esθ w5/2 0
(104)
6πc3 2 k w5/2 0
φr|drift = 0 where Equation (107) gives c3 :
(105)
SC
ρ¯drift = −
M AN U
0 c3 = (3k0 − w) ρ¯0 + esθ0 ρ¯0 + k02 θr0
(106)
(107)
The periodic condition, i.e. ρ¯drift = θr|drift = φr|drift = 0, results in c3 = 0 [21] , and leads to Equations (108) and (109): 0 θr0 =−
1 ((3k0 − w) ρ¯0 + esθ0 ρ¯00 ) k02
(109)
TE D
1 θ˙r0 = − 2 (3k0 − w) θ˙0 ρ¯0 + esθ0 ρ¯˙ 0 k0 ! esθ0 1 (3k0 − w) r0 θ˙0 − esθ0 r˙0 ρ0 + =− 2 ρ˙ 0 k0 r02 r0
(108)
Then the periodic solutions to the LTH equations can be obtained using Equation (110):
ρ = γx sin (θ + αx )
EP
AC C
where
θr =
1 + k0 γx cos (θ + αx ) + γy p
(110)
φr = γz sin (θ + αz ) q ρ20 + ρ00 2 γ = x 1 + k0 0 γy = θr0 − ρ0 p q γz = φ2r0 + φ0 2r0 ρ0 αx = arctan ρ0 − θ0 0 φ r0 αz = arctan 0 − θ0 φr0 20
(111)
ACCEPTED MANUSCRIPT
into Equation (112) [22]: ∂2a h i ∆r 1 2 ∂a ∂a ∆a = ∂r ∂v + [∆r, ∆v] ∂r2 ∂ a 2 ∆v ∂r∂v
∂2a ∂v 2
∆v
=
2a2 r2
∂a ∂v
4a3 v 2 µr 3
∂2a ∂r∂v
=
2a2 v µ 8a3 v µr 2
=
∂2a ∂v 2
=
+ ···
(112)
(113)
2a2 (8a−3r) µr
SC
∂2a ∂r 2
=
∂2a ∂r∂v ∆r
where ∂a ∂r
RI PT
The semi-major axis between the chief and deputy satellite can be expanded
Denote the n order expansion of variable ξ as ξ (n) , i.e. ξ = ξ (1) + ξ (2) + · · · ,
M AN U
then we have Equation (114):
∆a = ∆a(1) + ∆a(2) + · · ·
(114)
i ∆r(1) 2 ∂a = 2a µ∆r(1) + vr2 ∆v (1) ∂v µr2 ∆v (1)
(115)
where ∆a(1) =
∂a ∂r
∂a ∂r
i ∆r(2) h ∂a + 1 ∆r(1) ∂v 2 ∆v (2)
TE D
∆a(2) =
h
h
EP
2a2 2a2 v a2 = 2 ∆r(2) + ∆v (2) + 3 r µ µr
i
∆v (1)
2
2av (∆r
∂2a ∂r 2 ∂2a ∂r∂v
(1) 2
∂2a (1) ∂r∂v ∆r ∂2a ∂v 2
) + 8arv∆r
∆v (1)
(1)
∆v
(1)
2
+ r (8a − 3r) ∆v (116)
The precise periodic condition should satisfy the energy match condition as
AC C
defined by Equation (117) [6]: ∆a = 0
(117)
The first order approximate periodic condition is expressed as Equation
130
(118):
∆a(1) = 0
(118)
which gives Equation (119): ∆v (1) = −
µ ∆r(1) vr2
21
(119)
(1)
2
ACCEPTED MANUSCRIPT
In spherical coordinates, the deputy’s radius and velocity can be expressed
RI PT
by Equations (120) and (121): rb = ra + ρ
2
θ˙a + θ˙r
2
Consequently, we have Equation (122):
M AN U
∆r = ∆r(1) = ρ
cos2 φr + φ˙ 2r
SC
2
vb 2 = (r˙a + ρ) ˙ + (ra + ρ)
(120)
(121)
(122)
Denote the satellite velocity as Equation (123):
2 v 2 = r˙ 2 + r2 θ˙in cos2 φr + φ˙ 2r
(123)
where θin is the true anomaly’s projection on the chief’s orbital plane, and
TE D
φr is the projection vertical to the plane. And the state parameters are x = T r, θin , φr , r, ˙ θ˙in , φ˙ r . Then ∆v can be expanded as Equation (124): ∆v = ∆v (1) + ∆v (2) + · · · =
The first and second order partial derivatives
∂v ∂x
and
∂2v ∂2x
are attached in
EP
the appendix.
(124)
2 1 ∂v T ∂ v ∆x + (∆x) · 2 · ∆x + · · · ∂x 2 ∂ x
For the deputy’s relative motion with respect to the chief, we have Equation
AC C
(125):
T ˙ x = r, θ, 0, r, ˙ θ, 0 a T ∆x = ρ, θr , φr , ρ, ˙ θ˙r , φ˙ r x = x + ∆x b a
(125)
Substitute the above condition into the partial derivatives in the appendix,
then
∂v ∂x
and
∂2v ∂2x
can be simplified to Equations (126) and (127): ∂v r r˙ r2 ˙ = [ θ˙2 , 0, 0, , θ, 0] ∂x v v v 22
(126)
ACCEPTED MANUSCRIPT
0 0 ∂2v 1 = 3 ∂2x v −rr˙ θ˙2 2 v + r˙ 2 rθ˙ 0
0
0
−rr˙ θ˙2
v 2 + r˙ 2 rθ˙
0
0
0
0
0
−2r2 v 2 θ˙2
0
0
2 ˙2
0
0
r θ
−r2 r˙ θ˙
0
0
−r2 r˙ θ˙
r2 r˙ 2
0
0
0
0
Therefore, we have Equations (128) and (129):
2 1 T ∂ v (∆x) · 2 · ∆x 2 ∂ x 2 2rθ˙ r2 ˙2 2 1 ˙ ρθ˙r + −2θ φr + φ˙ 2r = 3 r˙ θρ − rθ˙ρ˙ + rr˙ θ˙r + v v v
∆v (2) =
0 0 (127) 0 0 r2 v2
∂v 1 ˙2 ∆x = rθ ρ + r˙ ρ˙ + r2 θ˙θ˙r ∂x v
M AN U
∆v (1) =
0
RI PT
r˙ 2 θ˙2
SC
(128)
(129)
Substituting Equations (122) and (128) into Equation (119) will also result
TE D
in Equations (108) and (109). In fact, the periodic condition c3 = 0 for the linearized LTH equation equals the first order approximate periodic condition for the full nonlinear relative model, i.e. Equation (119), and therefore results in a non-periodic solution and causes a small relative drift [22]. Under the condition of c3 = 0 or Equation (119), the real semi-major axis difference is
EP
estimated with Equation (130): (2)
=a
AC C
∆a ≈ ∆a
2 ! 2a ∆r(2) 4a − 2r ∆v (2) a (3a + 2r) ∆r(1) + + (130) r r r v r (2a − r) r
and the real relative drift rate can be estimated with Equations (104), (104), and (106), where c3 can be estimated with Equation (131) [21]: 2a c3 ≈ ∆a w
(131)
For the spherical model, ∆r(2) = 0, and there is no along-track term θr in 135
∆v (2) and ∆a(2) . However, for the Cartesian model, there is an along-track term y in ∆r(2) , ∆v (2) and ∆a(2) ( see Equations (57), (59), and (62) in reference 23
ACCEPTED MANUSCRIPT
[21]). Therefore, for the formation with large along-track distances, i.e. y x,
RI PT
y z or equivalently θr ρ/r, θr φr , ∆a(2) in equation (62) of reference [22] has a higher order than ∆a(2) in Equation (130) in this paper. 140
Consider a special case when the coplanar leader-follower formation with only a different argument of perigee, i.e. aa = ab , ea = eb , ia = ib , Ωa = Ωb , ωa 6= ωb , and θa = θb . It fulfills the energy match condition in Equation (117).
SC
And its spherical relative state, i.e. ρ = 0, θr 6= 0, φr = 0, ρ˙ = 0, θ˙r = φ˙ r = 0,
fulfills the first order periodic condition in Equation (118). This means that 145
the coplanar leader-follower formation does not drift if designed with spherical
3. Numerical Examples
M AN U
coordinates.
This section evaluates the proposed solution’s accuracy for relative motion in both circular and elliptical reference orbit, including the formations with a 150
difference in a single orbital element, and with differences in combined orbital
TE D
elements. To compare the accuracy, the spherical solution is first transformed into LVLH Cartesian coordinates using Equations (1) and (2). Then the solutions’ error is evaluated with respect to the accurate LVLH Cartesian coordinates, which are transformed from the chief and deputy’s inertial positions that 155
are obtained with Keplerian propagation. The parameters of the circular and
EP
elliptical reference orbits are listed in Table 1. Table 1: orbital elements of the Chief Satellite Orbit
a/m
e
i/◦
Ω/◦
ω/◦
θ/◦
Circular
1.0 × 107
0
40
10
0
0
Elliptical
1.0 × 107
0.5
40
10
0
0
AC C
Orbital elements
24
ACCEPTED MANUSCRIPT
3.1. Example 1. Error analysis for formations with differences in orbital ele-
RI PT
ments This example evaluates the solution’s accuracy for relative motion of forma160
tions with a difference in a single orbital element. Table 2 presents the differ-
ences of the orbital elements of different cases. Table 3 and Table 4 present the long term relative state errors for the formation in a circular reference orbit and
SC
an elliptical reference orbit, respectively. The first column gives the formation
cases, and the other rows show the corresponding components’ maximum errors 165
for solutions based on spherical coordinates and Cartesian coordinates.
M AN U
In Table 3 for the circular reference orbit, the first row is a formation with ∆a = -100 m. The spherical solution has a periodically varying error with a maximum of 0.0092 m in the radial (x) direction and a periodically increasing error of 0.0471 m per orbit period in the along-track (y) direction. The Cartesian 170
solution has an accelerated increase in error in the radial direction (greater than 0.0471m for the first period, as shown in Fig. 3), and a periodically increasing error of 0.0093 m per orbit period in the along-track direction. The long-term
TE D
accuracy of the spherical solution is much better than the Cartesian solution because of significant priority of the spherical solution in the radial direction. 175
For the formation with differences in other elements, the Cartesian solution has a periodically varying error in the radial direction and periodically increasing
EP
error in the along-track direction. The spherical solution has slightly worse accuracy when ∆e 6= 0, the same accuracy when ∆i 6= 0, and a little better accuracy when ∆Ω 6= 0. For the leader-follower formation when ∆ω 6= 0 or ∆θ 6= 0, the spherical solution is accurate in the radial and the along-track
AC C
180
direction. However, for the Cartesian solution, the radial and along-track errors are nearly quadratically related to the phase difference. In Table 4 for the elliptical reference orbit, when ∆a 6= 0, the spherical
solution performs worse than the Cartesian solution in the along-track direction,
185
but the 3-dimensional accuracy has the same magnitude. Similar to the cases in the circular reference orbit, the spherical solution has worse accuracy when ∆e 6= 0, same accuracy when ∆i 6= 0, and a little better accuracy when ∆Ω 6= 0. 25
ACCEPTED MANUSCRIPT
When ∆ω 6= 0 or ∆θ 6= 0, the errors of the Cartesian solution in the radial
190
RI PT
and along-track directions are more significant than the cases in the circular reference orbit because of the eccentricity of the reference orbit. When ∆ω 6= 0, the spherical solution is accurate. When ∆θ 6= 0, the spherical solution’s errors
periodically increase in the radial and the along-track directions, but are still one order lower than the Cartesian solution.
195
SC
In all the simulation cases, the spherical and Cartesian solutions are both accurate in the normal direction.
From this example, it is shown that the spherical solution has a significant
M AN U
accuracy advantage in the leader-follower formation (∆ω 6= 0 or ∆θ 6= 0). And in the cases when ∆Ω 6= 0, the accuracies of the spherical solutions are slightly better.
Table 2: Orbital elements’ differences of the Formations
∆a/m
formation with ∆a
−100
formation with ∆e
0.00
∆e
∆i/◦
∆Ω/◦
∆ω/◦
∆θ/◦
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0000
0.0000
0.0000
0.0000
TE D
Cases
0.0000
0.0057
0.0000
0.0000
0.0000
formation with ∆Ω
0.00
0.0000
0.0000
0.0057
0.0000
0.0000
formation with ∆ω
0.00
0.0000
0.0000
0.0000
0.0057
0.0000
formation with ∆θ
0.00
0.0000
0.0000
0.0000
0.0000
0.0057
EP
0.00
3.2. Example 2. Analysis for leader-follower formations
AC C
200
formation with ∆i
Example 1 shows the accuracy advantage of spherical solutions in the leader-
follower formation. Example 2 evaluates the solutions’ accuracy in the leaderfollower formation (∆ω 6= 0 or ∆θ 6= 0) with differences in other orbital elements. The circular and elliptical reference orbital parameters are the same as
205
example 1. Table 5 lists the differences for the orbital elements in two cases: the in-plane drift case with ∆a and ∆ω, and the 3-D periodic case with ∆θ, ∆i,
∆Ω and ∆Ω. The relative trajectories of the two cases in circular and elliptical
26
RI PT
ACCEPTED MANUSCRIPT
Table 3: Maximum errors of relative state components for formations in the circular orbit
Spherical
Cartesian
Along-track
Normal
Radial
Along-track
Normal
formation with ∆a
0.0092m
+0.0471m/p
0m
See Fig. 3
+0.0093m/p
0m
formation with ∆e
0.9997m
+4.7125m/p
0m
0.2248m
+0.9425m/p
0m
formation with ∆i
0.1980m
+0.9330m/p
0m
0.1980m
+0.9330m/p
0m
formation with ∆Ω
0.0823m
+0.3868m/p
0m
0.2147m
+1.4800m/p
0m
formation with ∆ω
0m
0m
0m
0.2969m
+1.8660m/p
0m
formation with ∆θ
0m
0m
0m
0.2969m
+1.8660m/p
0m
TE D
M AN U
Radial
SC
Cases
Table 4: Maximum errors of relative state components for formations in the elliptical orbit
Spherical
Orbit element difference
Cartesian
Along-track
Normal
Radial
Along-track
Normal
formation with ∆a
See Fig. 4
See Fig. 4
0m
See Fig. 4
+0.0389m/p
0m
formation with ∆e
+27.66m/p
+82.71m/p
0m
+5.842m/p
+17.42m/p
0m
AC C
EP
Radial
formation with ∆i
+1.6190m/p
+4.8490m/p
0m
+1.6190m/p
+4.8490m/p
0m
formation with ∆Ω
0.6698m
+2.0030m/p
0m
+2.0320m/p
+6.0770m/p
0m
formation with ∆ω
0m
0m
0m
+2.7000m/p
+8.0780m/p
0m
formation with ∆θ
+0.1802m/p
+0.5385m/p
0m
+1.8020m/p
+5.3850m/p
0m
27
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
Figure 3: Relative state errors for the deputy with a semi-major axis difference with respect
AC C
EP
TE D
to the circular reference orbit
Figure 4: Relative state errors for the deputy with a semi-major axis difference with respect to the elliptical reference orbit
28
ACCEPTED MANUSCRIPT
reference orbits are presented in Fig. 5 through 8, and the maximum errors of
210
RI PT
the components are listed in Table 6. It can be seen that the spherical solution’s accuracy are one to two orders higher than the Cartesian solution.
Table 5: Orbital elements difference of the leader-follower formations
∆a/m
∆e
∆i/◦
∆Ω/◦
∆ω/◦
∆θ/◦
Coplanar Drift Case
−100
0.0000
0.0000
0.0000
0.0573
0.0000
Non-coplanar Periodical Case
0.0000
0.0001
0.0057
SC
Orbit element difference
0.0057
0.0000
0.0573
M AN U
Table 6: Maximum errors of relative state components for the leader-follower formation
Spherical
Orbit element difference Radial 0.0083m
Circular - Non-coplanar Periodical
+0.8874m/p
Elliptical - Coplanar Drift
0.9259m
Elliptical - Non-coplanar Periodical
+8.395m/p
4. Conclusion
Along-track
Normal
Radial
Along-track
Normal
+0.0471m/p
0m
13.87m
+188.5m/p
0m
+1.2360m/p
0m
+271.0m/p
+816.3m/p
0m
+4.163m/p
0.3637m
34.790m
+218.7m/p
0.2728m
+25.27m/p
1.7270m
+213.10m/p
+637.2m/p
1.6900m
TE D
Circular - Coplanar Drift
Cartesian
EP
This paper establishes a linear model using spherical coordinates for relative motion in an elliptical orbit. First, we present the definition of spherical coordinates and their mapping between typical inertial orbital parameters. Then a linear model is derived through a strictly dynamical approach. The linear
AC C
215
equation has the same form as the Lawden or Tschauner-Hempel equation and can be easily solved with solutions from earlier literature. Then the model’s relative drift rate, periodic solution, and accuracy are investigated. Both the theoretical derivation and numerical examples show the model’s superior accu-
220
racy. For leader-follower formation with only a difference in the argument of perigee, the linear model is accurate. For formation with a large along-track distance, both the theory and the examples show the orders’ superior accuracy 29
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
AC C
EP
TE D
Figure 5: Relative trajectory for coplanar drift leader-follower formation in a circular orbit
Figure 6: Relative trajectory for coplanar drift leader-follower formation in an elliptical orbit
30
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
AC C
EP
TE D
(a) x-z plane
(b) y-z plane
Figure 7: Relative trajectory for con-coplanar periodic leader-follower formation in a circular orbit
31
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
AC C
EP
TE D
(a) x-z plane
(b) y-z plane
Figure 8: Relative trajectory for con-coplanar periodic leader-follower formation in an elliptical orbit
32
ACCEPTED MANUSCRIPT
compared with the linear model based on Cartesian coordinates. The proposed
225
RI PT
model relaxes the approximate limits of close relative motion, and shows advantages in the long term relative motion propagation. It can be applied in the evolution of a debris cloud, and can improve the accuracy of terminal guidance and fuel efficiency of interception or rendezvous.
SC
Acknowledgments
This research is supported by National Natural Science Foundation of China (61703380).
AC C
EP
TE D
M AN U
230
33
ACCEPTED MANUSCRIPT
Appendix: First and Second Order Derivatives of the Velocity
RI PT
Listed here are the derivatives of the velocity equations mentioned above.
AC C
EP
TE D
M AN U
SC
∂v r ˙2 θin cos2 φr + φ˙ 2r = ∂r v ∂v =0 ∂θin ∂v r2 2 = − θ˙in sin 2φr ∂φr 2v r˙ ∂v = ∂ r˙ v ∂v r2 = θ˙in cos2 φr v ∂ θ˙in ∂v r2 ˙ = φr v ∂ φ˙ r 2 2 ∂ v 1 2 ˙2 2 2 2 ˙2 2 2 ˙ ˙ = v θ cos φ + φ − r θ cos φ + φ r r in r in r ∂2r v3 ∂2v =0 ∂r∂θin r3 2 ∂2v r 2 2 sin 2φr + 3 θ˙in = − θ˙in cos2 φr + φ˙ 2r θ˙in sin 2φr ∂r∂φr v 2v rr˙ 2 ∂2v = − 3 θ˙in cos2 φr + φ˙ 2r ∂r∂ r˙ v ∂2v r r3 2 = 2 θ˙in cos2 φr − 3 θ˙in cos2 φr + φ˙ 2r θ˙in cos2 φr v v ∂r∂ θ˙in 3 2 r ∂ v r 2 cos2 φr + φ˙ 2r φ˙ r = 2 φ˙ r − 3 θ˙in v v ∂r∂ φ˙ r 2 2 ∂ v ∂ v ∂2v ∂2v ∂2v ∂2v = 2 = = =0 = = ∂θin ∂r ∂ θin ∂θin ∂φr ∂θin ∂ r˙ ∂θin ∂ θ˙in ∂θin ∂ φ˙ r ∂2v r 2 r3 2 2 = − θ˙in sin 2φr + 3 θ˙in cos2 φr + φ˙ 2r θ˙in sin 2φr ∂φr ∂r v 2v ∂2v =0 ∂φr ∂θin
(132)
(133)
(134)
(135) (136) (137) (138) (139) (140) (141) (142) (143) (144) (145) (146) (147)
34
ACCEPTED MANUSCRIPT
(148)
EP
TE D
M AN U
SC
RI PT
∂2v r2 ˙2 r4 ˙4 = −2 θ cos 2φ − θ sin2 2φr r ∂ 2 φr v in 4v 3 in ∂2v r2 r˙ 2 = 3 θ˙in sin 2φr ∂φr ∂ r˙ 2v r2 r4 3 ∂2v sin 2φr cos2 φr = − θ˙in sin 2φr + 3 θ˙in v 2v ∂φr ∂ θ˙in ∂2v r4 2 ˙ φr sin 2φr = 3 θ˙in 2v ∂φr ∂ φ˙ r ∂2v rr˙ 2 cos2 φr + φ˙ 2r = − 3 θ˙in ∂ r∂r ˙ v ∂2v =0 ∂ r∂θ ˙ in ∂2v r2 r˙ 2 sin 2φr = 3 θ˙in ∂ r∂φ ˙ r 2v ∂2v 1 r˙ 2 = − 3 2 ∂ r˙ v v ∂2v r2 r˙ = − 3 θ˙in cos2 φr v ∂ r∂ ˙ θ˙in 2 ∂ v r2 r˙ = − 3 φ˙ r v ∂ r∂ ˙ φ˙ r ∂2v r˙ r3 2 = 2 θin cos2 φr − 3 θ˙in cos2 φr + φ˙ 2r θ˙in cos2 φr v v ∂ θ˙in ∂r 2 ∂ v =0 ∂ θ˙in ∂θin r4 3 r2 ∂2v sin 2φr cos2 φr = − θ˙in sin 2φr + 3 θ˙in v 2v ∂ θ˙in ∂φr
AC C
∂2v r2 r˙ = − 3 θ˙in cos2 φr v ∂ θ˙in ∂ r˙ 2 2 r ∂ v r4 2 = cos2 φr − 3 θ˙in cos4 φr v v ∂ 2 θ˙in r4 ∂2v = − 3 θ˙in φ˙ r cos2 φr v ∂ θ˙in ∂ φ˙ r ∂2v r ˙ r3 2 = 2 φr − 3 θ˙in cos2 φr + φ˙ 2r φ˙ r v v ∂ φ˙ r ∂r
(149) (150) (151)
(152) (153) (154) (155) (156) (157) (158) (159) (160) (161) (162) (163) (164)
∂2v =0 ∂ φ˙ r ∂θin
(165)
∂2v r4 2 ˙ = 3 θ˙in φr sin 2φr 2v ∂ φ˙ r ∂φr
(166) (167) 35
ACCEPTED MANUSCRIPT
∂2v r2 = − 3 r˙ φ˙ r v ∂ φ˙ r ∂ r˙ r4 ∂2v = − 3 θ˙in φ˙ r cos2 φr v ∂ φ˙ r ∂ θ˙in ∂2v r2 r4 = − 3 φ˙ 2r v v ∂ 2 φ˙ r
(169) (170)
SC
References
RI PT
(168)
[1] J. Sullivan, S. Grimberg, S. DAmico, Comprehensive survey and assessment
235
M AN U
of spacecraft relative motion dynamics models, Journal of Guidance Control & Dynamics 40 (8) (2017) 1–23. doi:10.2514/1.G002309. [2] G. W. Hill, Researches in the lunar theory, American Journal of Mathematics 1 (2) (1878) 5–26. doi:10.2307/2369430.
[3] W. H. Clohessy, Terminal guidance system for satellite rendezvous,
240
TE D
Aerospace Sci 29 (9) (1960) 79–86. doi:10.2514/8.8704. [4] D. F. Lawden, Optimal Trajectories for Space Navigation, Butterworths, 1963.
[5] J. Tschauner, P. Hempel, Rendezvous with a target in an elliptical orbit,
EP
Acta Astronautica 11 (2) (1965) 104–109. [6] K. T. Alfriend, S. R. Vadali, P. Gurfil, J. How, L. S. Breger, Spacecraft formation flying: Dynamics, control and navigation, Butterworths, 2010.
AC C
245
[7] T. F. Berreen, J. D. C. Crisp, An exact and a new first-order solution for the relative trajectories of a probe ejected from a space station, Celestial Mechanics 13 (1) (1976) 75–88. doi:10.1007/BF01228535.
[8] D. K. Geller, T. A. Lovell, Relative orbital motion and angles-only relative
250
state observability in cylindrical coordinates, in: Proceedings of the AASAIAA Space Flight Mechanics Meeting, Santa Fe, New Mexico, 2014.
36
ACCEPTED MANUSCRIPT
[9] F. D. Bruijn, E. Gill, J. P. How, Comparative analysis of cartesian and
RI PT
curvilinear clohessy-wiltshire equations, Journal of Aerospace Engineering, Sciences and Applications 3 (2) (2011) 1–15. doi:10.7446/jaesa.0302. 255
01.
[10] C. D. Karlgaard, F. H. Lutze, Second-order relative motion equations, Journal of Guidance Control & Dynamics 26 (1) (2001) 41–49.
SC
10.2514/2.5013.
doi:
[11] C. Bombardelli, J. L. Gonzalo, J. Roa, Compact solution of circular orbit relative motion in curvilinear coordinates, in: Astrodynamics Specialist Conference, 2016.
M AN U
260
[12] A. C. Perez, Applications of relative motion models using curvilinear coordinate frames, Ph.D. thesis (2017).
[13] K. Yamanaka, F. Ankersen, New state transition matrix for relative motion 265
on an arbitrary elliptical orbit, Journal of Guidance Control & Dynamics
TE D
25 (1) (2002) 60–66. doi:10.2514/2.4875.
[14] J. S. Ardaens, S. D’Amico, A. Cropp, Gps-based relative navigation for the proba-3 formation flying mission , Acta Astronautica 91 (2013) 341–355. doi:10.1016/j.actaastro.2013.06.025. [15] J. L. Burch, T. E. Moore, R. B. Torbert, B. L. Giles, Magnetospheric
EP
270
multiscale overview and science objectives, Space Science Reviews 199 (1-
AC C
4) (2016) 5–21. doi:10.1007/s11214-015-0164-9. [16] A. W. Koenig, T. Guffanti, S. DAmico, New state transition matrices for spacecraft relative motion in perturbed orbits, Journal of Guidance Control
275
& Dynamics 40 (7) (2017) 1749–1768. doi:10.2514/1.G002409.
[17] T. Carter, M. Humi, Fuel-optimal rendezvous near a point in general keplerian orbit, Journal of Guidance Control & Dynamics 10 (6) (1987) 567–573. doi:10.2514/3.20257.
37
ACCEPTED MANUSCRIPT
[18] T. E. Carter, New form for the optimal rendezvous equations near a keplerian orbit, Journal of Guidance Control & Dynamics 13 (1) (1990) 183–186. doi:10.2514/3.20533.
RI PT
280
[19] T. E. Carter, State transition matrices for terminal rendezvous studies:
Brief survey and new example, Journal of Guidance Control & Dynamics
285
SC
21 (1) (2015) 148–155. doi:10.2514/2.4211.
[20] Z. Dang, Solutions of tschaunerhempel equations, Journal of Guidance Control & Dynamics 40 (11) (2017) 2953–2957. doi:10.2514/1.G002774.
M AN U
[21] P. Sengupta, S. R. Vadali, Relative motion and the geometry of formations in keplerian elliptic orbits with arbitrary eccentricity, Journal of Guidance Control & Dynamics 30 (4) (2007) 953–964. doi:10.2514/1.25941. 290
[22] F. Jiang, J. Li, H. Baoyin, Approximate analysis for relative motion of satellite formation flying in elliptical orbits, Celestial Mechanics & Dynamical
AC C
EP
TE D
Astronomy 98 (1) (2007) 31–66. doi:10.1007/s10569-007-9067-8.
38
S0094-5765(18)30910-X
DOI:
https://doi.org/10.1016/j.actaastro.2019.01.016
Reference:
AA 7282
To appear in:
Acta Astronautica
Received Date: 1 June 2018 Revised Date:
21 November 2018
Accepted Date: 10 January 2019
Please cite this article as: C. Han, H. Chen, G. Alonso, Y. Rao, J. Cubas, J. Yin, X. Wang, A linear model for relative motion in an elliptical orbit based on a spherical coordinate system, Acta Astronautica (2019), doi: https://doi.org/10.1016/j.actaastro.2019.01.016. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
RI PT
A Linear Model for Relative Motion in an Elliptical Orbit Based On a Spherical Coordinate System Chao Hana , Huan Chena , Gustavo Alonsob , Yinrui Raoc , Javier Cubasb , Jianfeng Yind , Xiaohui Wanga,∗ a School
M AN U
SC
of Astronautics, Beihang University, No.37, Xueyuan Road, Haidian District, Beijing, 100191, China b Faculty of Aerospace Engineering, Technical University of Madrid, Madrid, 28040, Spain c Institute of Systems Engineering, China Academy of Engineering Physics, Mianyang, Sichuan Province, 621999, China d Beijing Institute of Spacecraft System Engineering, China Academy of Space Technology, Beijing, 100094, China
Abstract
When compared with Cartesian coordinates, curvilinear coordinates have shown significant superior accuracy when applied to a linear model for relative motion in a circular reference orbit. Experts assume the same superiority for an elliptical reference orbit. However, the curvilinear model in an elliptical reference
TE D
orbit has not been established, and the superiority has not been proved. After a strictly theoretical derivation this paper establishes a linear model with spherical coordinates for relative motion in an elliptical reference orbit. The differential equation has the same form as the Lawden or Tschauner-Hempel
EP
equation and can be easily solved with the solutions from the literature. In addition both the theoretical derivation and the numerical examples prove the model’s superior accuracy in leader-follower formation and formations with a
AC C
large along-track distance. The proposed model reduces the approximate limits of close relative motion, and it is expected to be applied in the propagation, guidance, and control for longer term relative motion with large along-track distances.
∗ Corresponding
author Email addresses: [email protected] (Chao Han), [email protected] (Huan Chen), [email protected] (Gustavo Alonso), [email protected] (Yinrui Rao), [email protected] (Javier Cubas), [email protected] (Jianfeng Yin), [email protected] (Xiaohui Wang)
Preprint submitted to Journal of LATEX Templates
January 12, 2019
ACCEPTED MANUSCRIPT
Keywords: formation flying, relative motion, linear model, spherical
RI PT
coordinate 2010 MSC: 00-01, 99-00
1. Introduction
Relative motion dynamics is a key space mission subject area, highly im-
SC
portant to efforts related to formation flying, spacecraft rendezvous, terminal
guidance, space debris evolution and collision avoidance. For the relative motion 5
in a central gravity field Joshua Sullivan [1] provides a solid review of several
M AN U
representative models, including the Hill [2] or Clohessy-Wiltshire[3] (HCW) equations for the circular reference orbit and Lawden [4] or Tschauner-Hempel [5] (LTH) equations for the elliptical reference orbit. Both the HCW and LTH equations describe the relative motion in Cartesian coordinates and they are lin10
earized based on the assumption of close relative motion. For HCW equations in the circular case Alfriend [6] analyzed the effects of linearizing the initial condi-
TE D
tion mapping between the inertial state and relative state, and he found that the linearized mapping would cause considerable radial error when the along-track distance is large, which will limit the models’s application. 15
Many researchers have devoted effort to solving the above problem. A feasible approach is to establish the relative motion model in curvilinear coordinates,
EP
such as cylindrical coordinates or spherical coordinates. So far as we know, Berreen, T. and Crisp, J [7] were the first to set up a first order relative mo-
AC C
tion model based on cylindrical coordinates for a circular reference orbit. They 20
applied the model to the mission of releasing a detector from the space station, which is on a circular orbit. Geller, D. K. [8] proved that using cylindrical coordinates can remove the nonlinear terms related to the along-track distance and cross-track velocity. The resulting linear model in the HCW form is valid for arbitrarily large along-track distances and cross-track velocities. Alfriend [6] de-
25
rived a first order model, which has the same form as the HCW equations, based on the spherical coordinate system for a circular reference orbit, . De Bruijin,
2
ACCEPTED MANUSCRIPT
F. [9] compared the accuracy of the first order relative motion models based on
RI PT
both curvilinear coordinates and Cartesian coordinates that were in the form of the HCW equation. He found that for a formation with a large along-track 30
distance, the model based on curvilinear coordinates showed superior accuracy and reduced the fuel cost of formation maintenance maneuvers. However, for other formations with small along-track distances, the model based on curvilin-
SC
ear coordinates had no superiority or performed even worse. Further research, also based on curvilinear coordinates, achieved better accuracy by obtaining the 35
analytical solution to the second order model [10] and the approximate analyt-
M AN U
ical solution to a highly nonlinear model [11], or by taking into account the J2 perturbation [12]. When researching relative motion in elliptical reference orbits, Yamanaka [13] claimed that using cylindrical coordinates can significantly lessen the proximity requirement, although this work actually has not been con40
ducted or examined. Because formation missions have recently been extended to elliptical orbits recently [14, 15, 16], it is of great importance to investigate
reference orbit.
TE D
the application of curvilinear coordinates to the relative motion in an elliptical
This paper has two major innovations. First, the linearized equation for 45
relative motion in an elliptical reference orbit is obtained with a strictly dynamic derivation based on spherical coordinates. Second, the theoretical derivation
EP
and numerical examples show that the model based on spherical coordinates has superior accuracy in leader-follower formation when compared with Cartesian coordinates. The manuscript is organized as follows. The first section is the introduction. The second section is the theoretical part. First, we present the
AC C
50
definition of spherical coordinates, and we present the mappings between the coordinates and the typical inertial orbital parameters, using both the forms of position-velocity vectors and the classical orbital elements. Then the differential equations for relative motion in an elliptical reference orbit are derived, whose
55
analytical solutions can be found in the literature [13, 17, 18, 19, 20] because they are similar in form to the LTH equation. In addition the solutions’ accuracy are investigated, and the theoretical derivation shows its superior accuracy in 3
ACCEPTED MANUSCRIPT
leader-follower formation. In the third section, numerical examples examine the
60
RI PT
model’s accuracy in comparison with the Cartesian model. The last section comes to a conclusion.
2. Relative Motion Model 2.1. Definition of Spherical Coordinates
SC
The reference satellite A is denoted as the chief, and another satellite B is
the deputy. In the following sections the subscript ∗ (∗ = a, b) indicates the chief or the deputy. In Figure 1 two satellites are currently located at points A and
M AN U
65
B, respectively. The ascending nodes of the two orbits are Ma and Mb . The orbit planes of A and B intersect at point J, and the interior angle between the orbit planes is ∆j. The angle separation between point J and the ascending node M∗ is γ∗ . The angle separation between point J and the satellite is β∗ . 70
Suppose the projection of the radius vector of satellite B onto the orbit plane
TE D
of satellite A intersects with satellite A’s orbit at point B 0 . Orbit Plane of Satellite B
B
zi
rˆb
ha
EP
hb
AC C
O
r
b a
A
r
rˆbi
J
rˆa
a
Orbit Plane of Satellite A
B
j
b
ib
ia
Equator
Mb
Ma
Figure 1: Definition of Spherica Coordinate
The spherical coordinates are defined similarly to those in reference [6]. The parameter ρ denotes the radius difference between the chief and deputy. Denote 4
ACCEPTED MANUSCRIPT
the angle difference between OA and OB 0 as θr , and the angle difference between
RI PT
OB and OB 0 as φr . The coordinate transformation between the Earth center orbital (ECO) coordinates of the chief and the deputy can be handled as follows: Let the chief’s ECO frame Soa rotates by θr around the z axis, and then rotate by −φr around the y axis.
SC
Soa −−→ [Lz (θr )] −−→ [Ly (-φr )] Sob
In this way, ρ, θr , and φr indicate the deputy’s relative motion in the radial direction, along-track, and normal direction of the chief’s orbit, respectively.
75
M AN U
The deputy’s relative state parameters with respect to the chief are ρ, θr , φr , ρ, ˙ θ˙r and φ˙ r . Here the dot on the top indicates the derivative with respect to time.
zo (z'o)
y'o
O z r r xy
TE D
A
yo
B
xo
B' x'o
Figure 2: Geometric Description of ECO Coordinate
EP
The corresponding Cartesian coordinates in the local vertical local horizontal
AC C
(LVLH) coordinate system can be expressed by Equation (1): x = (ra + ρ) cos φr cos θr − ra y = (ra + ρ) cos φr sin θr z = (r + ρ) sin φ a r
(1)
And the relative velocity can be obtained from Equation (1) as given by Equa-
5
ACCEPTED MANUSCRIPT
RI PT
tion (2): ˙ r sin φr cos θr + θ˙r cos φr sin θr − r˙a x ˙ = ( r ˙ + ρ) ˙ cos φ cos θ − (r + ρ) φ a r r a y˙ = (r˙a + ρ) ˙ cos φr sin θr − (ra + ρ) φ˙ r sin φr sin θr − θ˙r cos φr cos θr z˙ = (r˙a + ρ) ˙ sin φr + (ra + ρ) φ˙ r cos φr (2) Equations (1) and (2) will be used for accuracy comparison of the spherical
SC
and Cartesian models.
2.2. Mapping between the Inertial Radius-Velocity Vectors and Spherical Coordinates
M AN U
80
Here denote a vector variable with bold type as x = xˆ x, where x = kxk is ˆ with a hat is the unit vector along x. The inertial radius its square norm and x and velocity vectors are denoted as r∗ and v∗ (∗ = a, b).
2.2.1. Mapping from the inertial radius-velocity vectors to spherical coordinates 85
When the chief and deputy’s inertial radius and velocity vectors, which are
TE D
ra , va , rb and vb , are known, the deputy’s relative state parameters with respect to the chief, i.e. ρ, θr , φr , ρ, ˙ θ˙r and φ˙ r , remain to be calculated. The velocity projection in the radial direction is given by Equation (3): v∗r = v∗ · ˆr∗
(3)
EP
And the radial gradient is given by Equation (4): r˙∗ = v∗r
(4)
AC C
Then ρ and ρ˙ can be obtained from Equations (5) and (6): ρ = rb − ra
(5)
ρ˙ = vbr − var
(6)
The angular momentum vector of chief satellite is given by Equation (7): ha = ra × va 6
(7)
ACCEPTED MANUSCRIPT
ˆ a of chief’s The projections of rb and vb onto the perpendicular direction h
RI PT
orbit plane are given by Equations (8) and (9): ˆa h ˆa rbh = rb · h
(8)
ˆa h ˆa vbh = vb · h
(9)
(10) and (11):
M AN U
rbi = rb − rbh
SC
The projections of rb and vb in the chief’s orbit plane are given by Equations
vbi = vb − vbh
(10)
(11)
Then θr can be obtained from Equation (12):
θr = sign (θr ) · arccos (ˆra · ˆrbi )
(12)
where sign (θr ) is the sign of θr , which is positive when ra rotates clockwise to
TE D
ˆ a , and is determined by Equation (13): rbi around h sign (θr ) = sign
ˆa (ra × rbi ) · h
(13)
The angular velocities of vectors ra and rbi are given by Equations (14) and
AC C
EP
(15):
kra × va k θ˙a = ra2
(14)
krbi × vbi k θ˙bi = 2 rbi
(15)
Then θ˙r can be obtained from Equation (16) θ˙r = θ˙bi − θ˙a
(16)
rbh rb
(17)
From Equation (17) sin φr =
7
ACCEPTED MANUSCRIPT
and its derivative, Equation (18),
vbh − vbr sin φr φ˙ r = rb cos φr
SC
we obtain φr and φ˙ r from Equations (19) and (20) rbh φr = arcsin rb
(18)
RI PT
vbh rb − rbh vbr φ˙ r cos φr = rb2
(19)
(20)
M AN U
2.2.2. Mapping from spherical coordinates to the inertial radius-velocity vectors When the chief’s inertial radius and velocity vectors and the deputy’s relative 90
state, given by ra , va , ρ, θr , φr , ρ, ˙ θ˙r and φ˙ r , are known, the deputy’s inertial radius and velocity vector, i.e. rb and vb , remain to be calculated. First, from Equations (5) and (6), we have Equations (21) and (22):
TE D
rb = ra + ρ
vbr = var + ρ˙
(21)
(22)
Then from Equations (19) and (20), we have Equations (23) and (24): (23)
vbh = φ˙ r rb cos φr + vbr sin φr
(24)
AC C
EP
rbh = rb sinφr
Therefore, we produce Equation (25): rbi =
q
2 rb2 − rbh
(25)
And the unit vector along rbi can be obtained from Equation (26): ˆrbi = Lz (θr ) ˆra
8
(26)
ACCEPTED MANUSCRIPT
where Lz (θr ) is the coordinate transition matrix for rotation along the z axis
cos θr Lz (θr ) = sin θr 0
0 0 1
− sin θr cos θr 0
ˆ a + rbi ˆrbi rb = rbh h
SC
Consequently, we have Equation (28):
RI PT
as shown by Equation (27):
(27)
(28)
M AN U
And using Equation (16), we get Equations (29) and (30): θ˙bi = θ˙a + θ˙r
(29)
vbt = rbi θ˙bi
(30)
Finally, we obtain Equation (31):
(31)
TE D
vb = [vbr , vbt , vbh ]T
2.3. Mapping between the Classical Orbital Elements and Spherical Coordinates 2.3.1. Mapping from the classical orbital elements to spherical coordinates When the chief and deputy’s orbital elements, which are the semi-major axis a∗ , eccentricity e∗ , inclination i∗ , right ascension of ascending node Ω∗ ,
EP
95
argument of perigee ω∗ and true anomaly θ∗ , are known, the deputy’s relative
AC C
state parameters with respect to the chief, i.e. ρ, θr , φr , ρ, ˙ θ˙r and φ˙ r , remain to be calculated.
For the Keplerian motion, we use Equation (32): r∗ =
p∗ 1 + e∗ cos θ∗
(32)
where p is the semi-latus rectum given by Equation (33): p∗ = a∗ 1 − e2∗
9
(33)
ACCEPTED MANUSCRIPT
(34)
RI PT
The radial velocity v∗r is given by Equation (34): r µ v∗r = r˙∗ = e∗ sin θ∗ p∗
Then ρ and ρ˙ can be obtained from Equations (5) and (6).
In addition the transverse velocity v∗u is given by Equation (35): h∗ v∗u = r∗ θ˙∗ = r∗
SC
(35)
where h∗ is the module of the angular momentum as given by Equation (36): √
µp∗
M AN U
h∗ =
(36)
Consequently, θ˙∗ can be obtained from Equation (37): h∗ θ˙∗ = 2 r∗
(37)
Spherical geometry gives Equations (38) and (39):
TE D
cos ∆j = − cos ia cos (π − ib ) + sin ia sin (π − ib ) cos ∆Ω
sin γa sin γb sin ∆Ω = = sin ∆j sin (π − ib ) sin ia
(38)
(39)
The parameters ∆j, γa and γb can be obtained from Equations (40), (41),
EP
and (42):
AC C
∆j = sign (∆i) · arccos (cos ia cos ib + sin ia sin ib cos ∆Ω)
100
γa = arcsin γb = arcsin
sin ∆Ω sin ib sin ∆j
sin ∆Ω sin ia sin ∆j
(40)
(41)
(42)
For Equations (41) and (42), when ∆j = 0, that is when ∆i = ∆Ω = 0, −−→ −−→ vectors OB and OB 0 coincide, and thus γa = γb = 0.
10
It’s obvious that Equation (43) holds: βa = ωa + θa − γa
(43)
βb = ωb + θb − γb
Spherical geometry gives Equations (44) and (45): tan βb0 = tan βb cos ∆j
(44)
SC
sin φr = sin βb sin ∆j
RI PT
ACCEPTED MANUSCRIPT
(45)
M AN U
The parameters βb0 , θr and φr can be obtained from Equations (46), (47), and (48):
(46)
θr = βb0 − βa
(47)
φr = arcsin (sin βb sin ∆j)
(48)
TE D
βb0 = arctan (tan βb cos ∆j)
Using Equations (46), (47), and (48), we obtain the angles’ derivative as given by Equations (49), (50), and (51): (49)
θ˙r = β˙ b0 − θ˙a
(50)
cos βb φ˙ r = θ˙b sin ∆j cos φr
(51)
AC C
EP
2 cos2 βb0 ˙b sin βb0 1 cos ∆j = θ β˙ b0 = θ˙b cos2 βb sin2 βb cos ∆j
Here β˙ b0 is presented in two forms to avoid the singularity when βb =
N.
11
n 2 π,
n∈
ACCEPTED MANUSCRIPT
2.3.2. Mapping from spherical coordinates to the classical orbital elements When the chief’s orbital elements, i.e. aa , ea , ia , Ωa , ωa and θa , and the
RI PT
105
deputy’s relative state parameters, i.e. ρ, θr , φr , ρ, ˙ θ˙r and φ˙ r , are known, the deputy’s orbital elements, i.e. ab , eb , ib , Ωb , ωb and θb , are to be calculated. Firstly from Equation (50), we get Equation (52): β˙ b0 = θ˙a + θ˙r Spherical geometry gives Equation (53):
M AN U
tan φr = sin βb0 tan ∆j
SC
(52)
(53)
and its derivative, given by Equation (54):
φ˙ r = β˙ b0 cos βb0 tan ∆j cos2 φr
TE D
we obtain βb0 and ∆j with Equations (55) and (56), ! β˙ b0 sin 2φr βb0 = arctan 2φ˙ r v u u ∆j = arctan ttan2 φr +
φ˙ r β˙ b0 cos2 φr
(54)
(55)
!2 (56)
EP
Spherical geometry gives Equation (57): sin βb =
sin φr sin ∆j
AC C
Then βb can be obtained as Equation (58): sin φr βb = arcsin sin ∆j
(57)
(58)
Second, from Equations (47) and (43), we obtain Equations (59) and (60): βa = βb0 − θr
(59)
γa = ωa + θa − βa
(60)
12
ACCEPTED MANUSCRIPT
Then spherical geometry gives Equations (61) and (62): (61)
sin ∆Ω sin γa sin γb = = sin ∆j sin ib sin ia
(62)
RI PT
cos (π − ib ) = − cos ia cos ∆j + sin ia sin ∆j cos γa
We obtain ib , ∆Ω and γb using Equations (63), (64), and (65):
sin γa sin ∆j sin ib
γb = arcsin
sin γa sin ia sin ib
(63)
M AN U
∆Ω = arcsin
SC
ib = arccos (cos ia cos ∆j − sin ia sin ∆j cos γa )
(64)
(65)
Therefore, we get Equation (66):
ub = ωb + θb = γb + βb
(66)
TE D
Taking the derivative of Equations (66) and (57) provides Equation (67): θ˙b = β˙ b =
φ˙ r cos φr cos βb sin ∆j
(67)
Given Equation (56), when φr = 0 and φ˙ r = 0, ∆j = 0, and we find
EP
110
Equations (58) and (67) are singular. In this case, ∆i = ∆Ω = 0, and the −−→ −−→ vector OB and OB 0 coincide. Then Equations (58) and (67) can be replaced with βb = ωa + θa − θr and θ˙b = β˙ b = θ˙a − θ˙r .
AC C
From Equations (5) and (6), we get Equations (68) and (69): rb = ra + ρ
(68)
r˙b = r˙a + ρ˙
(69)
Substitute Equations (67) and (68) into Equations (37) and (36) to get
Equation (70): pb =
rb4 θ˙b2 µ
13
(70)
ACCEPTED MANUSCRIPT
Then substitute Equations (68), (69), and (70) into Equations (32) and (34)
s eb =
2
(pb − rb ) pb + r˙b2 2 rb µ
tan θb =
rb r˙b pb − rb
r
pb µ
RI PT
to get Equations (71) and (72): (71)
(72)
SC
Finally, ab and ωb can be obtained with Equations (33) and (66). 2.4. Relative Motion Equations
M AN U
ˆo, y ˆo The unit vectors of the chief’s ECO coordinate axes are denoted as x ˆo . Then the chief’s unit radial vector, angular velocity and acceleration and z
TE D
vectors can be expressed by Equation (73): ˆra = x ˆo ˆo wa = θ˙a z w ˙ a = θ¨a z ˆo
(73)
Using the coordinate transformation, the deputy’s unit radial vector can be
(74)
EP
expressed in the chief’s ECO coordinates as Equation (74): 1 cos φr cos θr ˆrb = Lz (−θr ) Ly (φr ) 0 = cos φr sin θr 0 sin φr
AC C
The deputy’s angular velocity is given by Equation (75): −φ˙ r sin θr wb = wa + wb/a = φ˙ r cos θr θ˙a + θ˙r
where
wb/a
0 0 −φ˙ r sin θr = 0 + Lz (−θr ) φ˙ r = φ˙ r cos θr θ˙r 0 θ˙r 14
(75)
(76)
ACCEPTED MANUSCRIPT
Then the deputy’s angular acceleration is given by Equation (77):
RI PT
˙b=w ˙a+w ˙ b/a = w ˙ a|a + wa × wa + w ˙ b/a|a + wa × wb/a w ˆo = − φ¨r sin θr + θ˙a + θ˙r φ˙ r cos θr x ˆ o + θ¨a + θ¨r z ˆo + φ¨r cos θr − θ˙a + θ˙r φ˙ r sin θr y
(77)
The chief’s radial vector and its first and second-order derivatives are given
rb = rb ˆrb
M AN U
r˙ b = r˙b ˆrb + rb wb × ˆrb
SC
by Equations (78), (79), and (80):
¨rb = r¨b ˆrb + 2r˙b wb × ˆrb + rb w ˙ b × ˆrb + rb wb × (wb × ˆrb )
(78)
(79)
(80)
Under the condition of close relative motion, Equation (80) can be linearized
TE D
to form Equation (81): r¨a + ρ¨ − 2r˙a θ˙a θr − ra θ˙a2 + 2θ˙a θ˙r + θ¨a θr − θ˙a2 ρ ¨rb = r¨a θr + 2 r˙a θ˙a + r˙a θ˙r + ρ˙ θ˙a + ra −θ˙a2 θr + θ¨a + θ¨r + θ¨a ρ r¨a φr + 2r˙a φ˙ r + ra φ¨r
(81)
The forces on the deputy include the central gravity gb and the control force
EP
fb as given by Equation (82):
¨rb = gb + fb
(82)
AC C
where the central gravity is given by Equation (83): cos φr cos θr µ µ gb = − 2 ˆrb = − cos φ sin θ r r 2 rb (ra + ρ) sin φr and gb can be linearized as Equation (84): 1 µ µ gb = − 3 (ra − 2ρ) θr = − 3 ra ra φr 15
ra − 2ρ ra θr ra φr
(83)
(84)
ACCEPTED MANUSCRIPT
For the chief with only central gravity we have Equation (85): ¨ra − ga = 0
RI PT
(85)
That is spelled out as Equation (86): µ r¨a − ra θ˙a2 + 2 = 0 ra 2r˙ θ˙ + r θ¨ = 0 a a
SC
a a
(86)
Because
d dt
d = θ˙ dθ ,
d2 dt2
M AN U
Subtracting Equation (82) from Equation (86) we get Equation (87): ρ¨ − 2θ˙a r˙a θr + ra θ˙r − θ˙a2 ρ − θ¨a ra θr −2ρ µ ¨ ra θr + 2r˙a θ˙r + r¨a θr + 2θ˙a ρ˙ − θ˙a2 ra θr + θ¨a ρ + 3 ra θr = fbc ra ra φr ra φ¨r + 2r˙a φ˙ r + r¨a φr d2 ˙ dθ˙ d = θ˙2 dθ 2 +θ dθ dθ , we can denote
00
d(∗) dθ
0
= (∗) ,
d2 (∗) dθ 2
=
(88)
TE D
(∗) , to give Equation (88): µ2 θ˙ = 3 (1 + e cos θ)2 h 2 0 θ˙ = − 2µ e (1 + e cos θ) sin θ 3 h
(87)
EP
Then Equation (87) can be rearranged to form Equation (89): 0 ((1 + ea cos θa ) ρ00 − 2ea sin θa ρ0 ) − (3 + ea cos θa ) ρ − 2 (1 + ea cos θa ) (ra θr ) + 2ea sin θa (ra θr ) µ 00 0 0 (1 + ea cos θa ) (ra θr ) − 2ea sin θa (ra θr ) − ea cos θa (ra θr ) + 2 (1 + ea cos θa ) ρ − 2ea sin θa ρ = fbc ra3 00 0 (1 + ea cos θa ) (ra φr ) − 2ea sin θa (ra φr ) + (ra φr ) (89)
AC C
In addition, we know that Equation (90) holds: ((1 + e cos θ) x)0 = (1 + e cos θ) x0 − e sin θx
00
(90)
((1 + e cos θ) x) = (1 + e cos θ) x00 − 2e sin θx0 − e cos θx
Substitute Equation (90) into (89) to get Equation (91): 00 0 ((1 + ea cos θa ) ρ) − 3ρ − 2((1 + ea cos θa ) ra θr ) µ 00 0 = fbc ((1 + e cos θ ) r θ ) + 2((1 + e cos θ ) ρ) a a a r a a 3 ra 00 ((1 + ea cos θa ) ra φr ) + (1 + ea cos θa ) ra φr 16
(91)
ACCEPTED MANUSCRIPT
Rearrange the above equation to get Equation (92):
µ h2a ra3 µ
ρ ra
00
3ρ ra (1+ea cos θa ) 0 θr00 + 2 rρa
−
φ00r + φr
− 2θr0
ρ ra
00
= ra θ˙a2
3ρ ra (1+ea cos θa ) 0 θr00 + 2 rρa
−
− 2θr0
RI PT
= fbc
φ00r + φr
(92)
In the following sections, the subscript ’a’ referring to the chief satellite is
M AN U
SC
omitted for simplicity. Denote ρ¯ = ρr and k = 1 + e cos θ to give Equation (93): 00 0 ρ¯ − 3¯ ρ/k − 2θr 2 ˙ 00 (93) rθ = fbc θr + 2¯ ρ0 φ00r + φr For the relative motion without the control force, i.e. fbc = 0, we have Equation (94):
θr00 + 2¯ ρ0 = 0
(94)
φ00r + φr = 0
TE D
ρ¯00 − 3¯ ρ/k − 2θr0 = 0
For the circular reference orbit, e = 0, θ˙ = n, θ¨ = 0, rb = a + ρ,
µ r3
= n2 ,
EP
where n is the mean motion. Equations (94) can be simplified to Equation (95): ρ¨ − 3n2 ρ − 2naθ˙r = 0 θ¨r + 2nρ/a ˙ =0 φ¨ + n2 φ = 0 r
(95)
AC C
Equations (94) and (95) have the same form as the LTH and HCW equations, 115
respectively. And the result of the circular case is the same as Equation (5.88) in reference [6], while ignoring the typo that the sign of 2n0 ρ˙ in the second equation of (5.88) should be positive. 2.5. Solutions to the Relative Motion The relative motion equations have the same form as the LTH equations, and many solutions have been proposed in previous research [13, 17, 18, 19, 20].
17
ACCEPTED MANUSCRIPT
is adopted in this paper. ρ¯ = c1 sθ + c2 cθ + c3 (2 − 3esθ L) θr = c1 qcθ − c2 qsθ − 2c3 k 2 L + c4 φr = c5 cos θ + c6 sin θ ρ¯0 = c1 s0θ + c2 c0θ − 3c3 e s0θ L + sθ /k 2
(96)
SC
RI PT
Yamanaka’s solution [13] is shown as Equation (96). His solution is concise and
θr0 = −2c1 sθ − c2 (2cθ − e) + 3c3 (2esθ L − 1) φ0r = −c5 sin θ + c6 cos θ
M AN U
where ci are the constant parameters determined by the initial state, and for the rest we have Equations (97) and (98): L=
(97)
q = 1 + 1/k sθ = k sin θ
TE D
µ2 (t − t0 ) h3
cθ = k cos θ
(98)
s0θ = cos θ + e cos 2θ
c0θ = − (sin θ + e sin 2θ)
EP
¯ = [¯ Denote the state parameters as x ρ, θr , φr , ρ¯0 , θr0 , φ0r ]T and the constant parameters in the solution as C = [c1 , c2 , c3 , c4 , c5 , c6 ]T , then the relative motion
AC C
solutions can be expressed as Equation (99) [13, 19, 20]: x ¯ = Φθ C
(99)
The state transition equation is given by Equation (100): ¯ 0 = Φθθ0 x x ¯ = Φθ Φ−1 ¯0 θ0 x
(100)
where the parameters Φθ , Φθ0 , and Φ−1 θ0 are given by Equations (101), (102),
18
ACCEPTED MANUSCRIPT
and (103): 2 − 3eLsθ
0
0
−qsθ
−3k 2 L
1
0
0
cos θ
0
0
0 c0θ
+ sθ /k
3 (2eLsθ − 1)
0
0
0
0
0
− sin θ
s θ qcθ 0 = 0 sθ −2sθ 0
cθ
2
0
−qsθ
0
1
0
0
c0θ
−3esθ /k
e − 2cθ
−3
0
0
2
(101)
θ
0
0
0
0
cos θ
0
0
0
0
0 − sin θ
0 sin θ 0 0 cos θ
0
0
cθ − 2e
−qsθ
0
0
−sθ
−e − qcθ
0
0
esθ
k2
w
0
ecθ − 2
−eqsθ
0
w cos θ
0
0
0
w sin θ
0
0
TE D
−3 1 + e2 /k sin θ −3 (e + cos θ) 3k + e2 − 1 1 = w −3eq sin θ 0 0
0 sin θ 0 0 cos θ
(102)
θ0
0
0 0 −w sin θ w cos θ 0
EP
Φ−1 θ0
−3e
2
e − 2cθ
Φθ0
0 Ls0θ
0
RI PT
qcθ 0 Φθ = 0 sθ −2sθ 0
cθ
SC
sθ
M AN U
θ0
(103)
AC C
where w = 1 − e2 . 120
2.6. Analysis of Solution Accuracy A closed relative orbit with no long-term relative drift is called a periodic
solution. In order to evaluate the accuracy of the proposed model, this section investigates the real relative drift of the periodic solution to the model. Following the same procedure as in reference [21], the relative motion drift
125
per orbit can be obtained with Equations (96) throuth (103) to give Equations
19
ACCEPTED MANUSCRIPT
(104), (105), and (106):
θr|drift = −
RI PT
6πc3 esθ w5/2 0
(104)
6πc3 2 k w5/2 0
φr|drift = 0 where Equation (107) gives c3 :
(105)
SC
ρ¯drift = −
M AN U
0 c3 = (3k0 − w) ρ¯0 + esθ0 ρ¯0 + k02 θr0
(106)
(107)
The periodic condition, i.e. ρ¯drift = θr|drift = φr|drift = 0, results in c3 = 0 [21] , and leads to Equations (108) and (109): 0 θr0 =−
1 ((3k0 − w) ρ¯0 + esθ0 ρ¯00 ) k02
(109)
TE D
1 θ˙r0 = − 2 (3k0 − w) θ˙0 ρ¯0 + esθ0 ρ¯˙ 0 k0 ! esθ0 1 (3k0 − w) r0 θ˙0 − esθ0 r˙0 ρ0 + =− 2 ρ˙ 0 k0 r02 r0
(108)
Then the periodic solutions to the LTH equations can be obtained using Equation (110):
ρ = γx sin (θ + αx )
EP
AC C
where
θr =
1 + k0 γx cos (θ + αx ) + γy p
(110)
φr = γz sin (θ + αz ) q ρ20 + ρ00 2 γ = x 1 + k0 0 γy = θr0 − ρ0 p q γz = φ2r0 + φ0 2r0 ρ0 αx = arctan ρ0 − θ0 0 φ r0 αz = arctan 0 − θ0 φr0 20
(111)
ACCEPTED MANUSCRIPT
into Equation (112) [22]: ∂2a h i ∆r 1 2 ∂a ∂a ∆a = ∂r ∂v + [∆r, ∆v] ∂r2 ∂ a 2 ∆v ∂r∂v
∂2a ∂v 2
∆v
=
2a2 r2
∂a ∂v
4a3 v 2 µr 3
∂2a ∂r∂v
=
2a2 v µ 8a3 v µr 2
=
∂2a ∂v 2
=
+ ···
(112)
(113)
2a2 (8a−3r) µr
SC
∂2a ∂r 2
=
∂2a ∂r∂v ∆r
where ∂a ∂r
RI PT
The semi-major axis between the chief and deputy satellite can be expanded
Denote the n order expansion of variable ξ as ξ (n) , i.e. ξ = ξ (1) + ξ (2) + · · · ,
M AN U
then we have Equation (114):
∆a = ∆a(1) + ∆a(2) + · · ·
(114)
i ∆r(1) 2 ∂a = 2a µ∆r(1) + vr2 ∆v (1) ∂v µr2 ∆v (1)
(115)
where ∆a(1) =
∂a ∂r
∂a ∂r
i ∆r(2) h ∂a + 1 ∆r(1) ∂v 2 ∆v (2)
TE D
∆a(2) =
h
h
EP
2a2 2a2 v a2 = 2 ∆r(2) + ∆v (2) + 3 r µ µr
i
∆v (1)
2
2av (∆r
∂2a ∂r 2 ∂2a ∂r∂v
(1) 2
∂2a (1) ∂r∂v ∆r ∂2a ∂v 2
) + 8arv∆r
∆v (1)
(1)
∆v
(1)
2
+ r (8a − 3r) ∆v (116)
The precise periodic condition should satisfy the energy match condition as
AC C
defined by Equation (117) [6]: ∆a = 0
(117)
The first order approximate periodic condition is expressed as Equation
130
(118):
∆a(1) = 0
(118)
which gives Equation (119): ∆v (1) = −
µ ∆r(1) vr2
21
(119)
(1)
2
ACCEPTED MANUSCRIPT
In spherical coordinates, the deputy’s radius and velocity can be expressed
RI PT
by Equations (120) and (121): rb = ra + ρ
2
θ˙a + θ˙r
2
Consequently, we have Equation (122):
M AN U
∆r = ∆r(1) = ρ
cos2 φr + φ˙ 2r
SC
2
vb 2 = (r˙a + ρ) ˙ + (ra + ρ)
(120)
(121)
(122)
Denote the satellite velocity as Equation (123):
2 v 2 = r˙ 2 + r2 θ˙in cos2 φr + φ˙ 2r
(123)
where θin is the true anomaly’s projection on the chief’s orbital plane, and
TE D
φr is the projection vertical to the plane. And the state parameters are x = T r, θin , φr , r, ˙ θ˙in , φ˙ r . Then ∆v can be expanded as Equation (124): ∆v = ∆v (1) + ∆v (2) + · · · =
The first and second order partial derivatives
∂v ∂x
and
∂2v ∂2x
are attached in
EP
the appendix.
(124)
2 1 ∂v T ∂ v ∆x + (∆x) · 2 · ∆x + · · · ∂x 2 ∂ x
For the deputy’s relative motion with respect to the chief, we have Equation
AC C
(125):
T ˙ x = r, θ, 0, r, ˙ θ, 0 a T ∆x = ρ, θr , φr , ρ, ˙ θ˙r , φ˙ r x = x + ∆x b a
(125)
Substitute the above condition into the partial derivatives in the appendix,
then
∂v ∂x
and
∂2v ∂2x
can be simplified to Equations (126) and (127): ∂v r r˙ r2 ˙ = [ θ˙2 , 0, 0, , θ, 0] ∂x v v v 22
(126)
ACCEPTED MANUSCRIPT
0 0 ∂2v 1 = 3 ∂2x v −rr˙ θ˙2 2 v + r˙ 2 rθ˙ 0
0
0
−rr˙ θ˙2
v 2 + r˙ 2 rθ˙
0
0
0
0
0
−2r2 v 2 θ˙2
0
0
2 ˙2
0
0
r θ
−r2 r˙ θ˙
0
0
−r2 r˙ θ˙
r2 r˙ 2
0
0
0
0
Therefore, we have Equations (128) and (129):
2 1 T ∂ v (∆x) · 2 · ∆x 2 ∂ x 2 2rθ˙ r2 ˙2 2 1 ˙ ρθ˙r + −2θ φr + φ˙ 2r = 3 r˙ θρ − rθ˙ρ˙ + rr˙ θ˙r + v v v
∆v (2) =
0 0 (127) 0 0 r2 v2
∂v 1 ˙2 ∆x = rθ ρ + r˙ ρ˙ + r2 θ˙θ˙r ∂x v
M AN U
∆v (1) =
0
RI PT
r˙ 2 θ˙2
SC
(128)
(129)
Substituting Equations (122) and (128) into Equation (119) will also result
TE D
in Equations (108) and (109). In fact, the periodic condition c3 = 0 for the linearized LTH equation equals the first order approximate periodic condition for the full nonlinear relative model, i.e. Equation (119), and therefore results in a non-periodic solution and causes a small relative drift [22]. Under the condition of c3 = 0 or Equation (119), the real semi-major axis difference is
EP
estimated with Equation (130): (2)
=a
AC C
∆a ≈ ∆a
2 ! 2a ∆r(2) 4a − 2r ∆v (2) a (3a + 2r) ∆r(1) + + (130) r r r v r (2a − r) r
and the real relative drift rate can be estimated with Equations (104), (104), and (106), where c3 can be estimated with Equation (131) [21]: 2a c3 ≈ ∆a w
(131)
For the spherical model, ∆r(2) = 0, and there is no along-track term θr in 135
∆v (2) and ∆a(2) . However, for the Cartesian model, there is an along-track term y in ∆r(2) , ∆v (2) and ∆a(2) ( see Equations (57), (59), and (62) in reference 23
ACCEPTED MANUSCRIPT
[21]). Therefore, for the formation with large along-track distances, i.e. y x,
RI PT
y z or equivalently θr ρ/r, θr φr , ∆a(2) in equation (62) of reference [22] has a higher order than ∆a(2) in Equation (130) in this paper. 140
Consider a special case when the coplanar leader-follower formation with only a different argument of perigee, i.e. aa = ab , ea = eb , ia = ib , Ωa = Ωb , ωa 6= ωb , and θa = θb . It fulfills the energy match condition in Equation (117).
SC
And its spherical relative state, i.e. ρ = 0, θr 6= 0, φr = 0, ρ˙ = 0, θ˙r = φ˙ r = 0,
fulfills the first order periodic condition in Equation (118). This means that 145
the coplanar leader-follower formation does not drift if designed with spherical
3. Numerical Examples
M AN U
coordinates.
This section evaluates the proposed solution’s accuracy for relative motion in both circular and elliptical reference orbit, including the formations with a 150
difference in a single orbital element, and with differences in combined orbital
TE D
elements. To compare the accuracy, the spherical solution is first transformed into LVLH Cartesian coordinates using Equations (1) and (2). Then the solutions’ error is evaluated with respect to the accurate LVLH Cartesian coordinates, which are transformed from the chief and deputy’s inertial positions that 155
are obtained with Keplerian propagation. The parameters of the circular and
EP
elliptical reference orbits are listed in Table 1. Table 1: orbital elements of the Chief Satellite Orbit
a/m
e
i/◦
Ω/◦
ω/◦
θ/◦
Circular
1.0 × 107
0
40
10
0
0
Elliptical
1.0 × 107
0.5
40
10
0
0
AC C
Orbital elements
24
ACCEPTED MANUSCRIPT
3.1. Example 1. Error analysis for formations with differences in orbital ele-
RI PT
ments This example evaluates the solution’s accuracy for relative motion of forma160
tions with a difference in a single orbital element. Table 2 presents the differ-
ences of the orbital elements of different cases. Table 3 and Table 4 present the long term relative state errors for the formation in a circular reference orbit and
SC
an elliptical reference orbit, respectively. The first column gives the formation
cases, and the other rows show the corresponding components’ maximum errors 165
for solutions based on spherical coordinates and Cartesian coordinates.
M AN U
In Table 3 for the circular reference orbit, the first row is a formation with ∆a = -100 m. The spherical solution has a periodically varying error with a maximum of 0.0092 m in the radial (x) direction and a periodically increasing error of 0.0471 m per orbit period in the along-track (y) direction. The Cartesian 170
solution has an accelerated increase in error in the radial direction (greater than 0.0471m for the first period, as shown in Fig. 3), and a periodically increasing error of 0.0093 m per orbit period in the along-track direction. The long-term
TE D
accuracy of the spherical solution is much better than the Cartesian solution because of significant priority of the spherical solution in the radial direction. 175
For the formation with differences in other elements, the Cartesian solution has a periodically varying error in the radial direction and periodically increasing
EP
error in the along-track direction. The spherical solution has slightly worse accuracy when ∆e 6= 0, the same accuracy when ∆i 6= 0, and a little better accuracy when ∆Ω 6= 0. For the leader-follower formation when ∆ω 6= 0 or ∆θ 6= 0, the spherical solution is accurate in the radial and the along-track
AC C
180
direction. However, for the Cartesian solution, the radial and along-track errors are nearly quadratically related to the phase difference. In Table 4 for the elliptical reference orbit, when ∆a 6= 0, the spherical
solution performs worse than the Cartesian solution in the along-track direction,
185
but the 3-dimensional accuracy has the same magnitude. Similar to the cases in the circular reference orbit, the spherical solution has worse accuracy when ∆e 6= 0, same accuracy when ∆i 6= 0, and a little better accuracy when ∆Ω 6= 0. 25
ACCEPTED MANUSCRIPT
When ∆ω 6= 0 or ∆θ 6= 0, the errors of the Cartesian solution in the radial
190
RI PT
and along-track directions are more significant than the cases in the circular reference orbit because of the eccentricity of the reference orbit. When ∆ω 6= 0, the spherical solution is accurate. When ∆θ 6= 0, the spherical solution’s errors
periodically increase in the radial and the along-track directions, but are still one order lower than the Cartesian solution.
195
SC
In all the simulation cases, the spherical and Cartesian solutions are both accurate in the normal direction.
From this example, it is shown that the spherical solution has a significant
M AN U
accuracy advantage in the leader-follower formation (∆ω 6= 0 or ∆θ 6= 0). And in the cases when ∆Ω 6= 0, the accuracies of the spherical solutions are slightly better.
Table 2: Orbital elements’ differences of the Formations
∆a/m
formation with ∆a
−100
formation with ∆e
0.00
∆e
∆i/◦
∆Ω/◦
∆ω/◦
∆θ/◦
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0000
0.0000
0.0000
0.0000
TE D
Cases
0.0000
0.0057
0.0000
0.0000
0.0000
formation with ∆Ω
0.00
0.0000
0.0000
0.0057
0.0000
0.0000
formation with ∆ω
0.00
0.0000
0.0000
0.0000
0.0057
0.0000
formation with ∆θ
0.00
0.0000
0.0000
0.0000
0.0000
0.0057
EP
0.00
3.2. Example 2. Analysis for leader-follower formations
AC C
200
formation with ∆i
Example 1 shows the accuracy advantage of spherical solutions in the leader-
follower formation. Example 2 evaluates the solutions’ accuracy in the leaderfollower formation (∆ω 6= 0 or ∆θ 6= 0) with differences in other orbital elements. The circular and elliptical reference orbital parameters are the same as
205
example 1. Table 5 lists the differences for the orbital elements in two cases: the in-plane drift case with ∆a and ∆ω, and the 3-D periodic case with ∆θ, ∆i,
∆Ω and ∆Ω. The relative trajectories of the two cases in circular and elliptical
26
RI PT
ACCEPTED MANUSCRIPT
Table 3: Maximum errors of relative state components for formations in the circular orbit
Spherical
Cartesian
Along-track
Normal
Radial
Along-track
Normal
formation with ∆a
0.0092m
+0.0471m/p
0m
See Fig. 3
+0.0093m/p
0m
formation with ∆e
0.9997m
+4.7125m/p
0m
0.2248m
+0.9425m/p
0m
formation with ∆i
0.1980m
+0.9330m/p
0m
0.1980m
+0.9330m/p
0m
formation with ∆Ω
0.0823m
+0.3868m/p
0m
0.2147m
+1.4800m/p
0m
formation with ∆ω
0m
0m
0m
0.2969m
+1.8660m/p
0m
formation with ∆θ
0m
0m
0m
0.2969m
+1.8660m/p
0m
TE D
M AN U
Radial
SC
Cases
Table 4: Maximum errors of relative state components for formations in the elliptical orbit
Spherical
Orbit element difference
Cartesian
Along-track
Normal
Radial
Along-track
Normal
formation with ∆a
See Fig. 4
See Fig. 4
0m
See Fig. 4
+0.0389m/p
0m
formation with ∆e
+27.66m/p
+82.71m/p
0m
+5.842m/p
+17.42m/p
0m
AC C
EP
Radial
formation with ∆i
+1.6190m/p
+4.8490m/p
0m
+1.6190m/p
+4.8490m/p
0m
formation with ∆Ω
0.6698m
+2.0030m/p
0m
+2.0320m/p
+6.0770m/p
0m
formation with ∆ω
0m
0m
0m
+2.7000m/p
+8.0780m/p
0m
formation with ∆θ
+0.1802m/p
+0.5385m/p
0m
+1.8020m/p
+5.3850m/p
0m
27
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
Figure 3: Relative state errors for the deputy with a semi-major axis difference with respect
AC C
EP
TE D
to the circular reference orbit
Figure 4: Relative state errors for the deputy with a semi-major axis difference with respect to the elliptical reference orbit
28
ACCEPTED MANUSCRIPT
reference orbits are presented in Fig. 5 through 8, and the maximum errors of
210
RI PT
the components are listed in Table 6. It can be seen that the spherical solution’s accuracy are one to two orders higher than the Cartesian solution.
Table 5: Orbital elements difference of the leader-follower formations
∆a/m
∆e
∆i/◦
∆Ω/◦
∆ω/◦
∆θ/◦
Coplanar Drift Case
−100
0.0000
0.0000
0.0000
0.0573
0.0000
Non-coplanar Periodical Case
0.0000
0.0001
0.0057
SC
Orbit element difference
0.0057
0.0000
0.0573
M AN U
Table 6: Maximum errors of relative state components for the leader-follower formation
Spherical
Orbit element difference Radial 0.0083m
Circular - Non-coplanar Periodical
+0.8874m/p
Elliptical - Coplanar Drift
0.9259m
Elliptical - Non-coplanar Periodical
+8.395m/p
4. Conclusion
Along-track
Normal
Radial
Along-track
Normal
+0.0471m/p
0m
13.87m
+188.5m/p
0m
+1.2360m/p
0m
+271.0m/p
+816.3m/p
0m
+4.163m/p
0.3637m
34.790m
+218.7m/p
0.2728m
+25.27m/p
1.7270m
+213.10m/p
+637.2m/p
1.6900m
TE D
Circular - Coplanar Drift
Cartesian
EP
This paper establishes a linear model using spherical coordinates for relative motion in an elliptical orbit. First, we present the definition of spherical coordinates and their mapping between typical inertial orbital parameters. Then a linear model is derived through a strictly dynamical approach. The linear
AC C
215
equation has the same form as the Lawden or Tschauner-Hempel equation and can be easily solved with solutions from earlier literature. Then the model’s relative drift rate, periodic solution, and accuracy are investigated. Both the theoretical derivation and numerical examples show the model’s superior accu-
220
racy. For leader-follower formation with only a difference in the argument of perigee, the linear model is accurate. For formation with a large along-track distance, both the theory and the examples show the orders’ superior accuracy 29
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
AC C
EP
TE D
Figure 5: Relative trajectory for coplanar drift leader-follower formation in a circular orbit
Figure 6: Relative trajectory for coplanar drift leader-follower formation in an elliptical orbit
30
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
AC C
EP
TE D
(a) x-z plane
(b) y-z plane
Figure 7: Relative trajectory for con-coplanar periodic leader-follower formation in a circular orbit
31
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
AC C
EP
TE D
(a) x-z plane
(b) y-z plane
Figure 8: Relative trajectory for con-coplanar periodic leader-follower formation in an elliptical orbit
32
ACCEPTED MANUSCRIPT
compared with the linear model based on Cartesian coordinates. The proposed
225
RI PT
model relaxes the approximate limits of close relative motion, and shows advantages in the long term relative motion propagation. It can be applied in the evolution of a debris cloud, and can improve the accuracy of terminal guidance and fuel efficiency of interception or rendezvous.
SC
Acknowledgments
This research is supported by National Natural Science Foundation of China (61703380).
AC C
EP
TE D
M AN U
230
33
ACCEPTED MANUSCRIPT
Appendix: First and Second Order Derivatives of the Velocity
RI PT
Listed here are the derivatives of the velocity equations mentioned above.
AC C
EP
TE D
M AN U
SC
∂v r ˙2 θin cos2 φr + φ˙ 2r = ∂r v ∂v =0 ∂θin ∂v r2 2 = − θ˙in sin 2φr ∂φr 2v r˙ ∂v = ∂ r˙ v ∂v r2 = θ˙in cos2 φr v ∂ θ˙in ∂v r2 ˙ = φr v ∂ φ˙ r 2 2 ∂ v 1 2 ˙2 2 2 2 ˙2 2 2 ˙ ˙ = v θ cos φ + φ − r θ cos φ + φ r r in r in r ∂2r v3 ∂2v =0 ∂r∂θin r3 2 ∂2v r 2 2 sin 2φr + 3 θ˙in = − θ˙in cos2 φr + φ˙ 2r θ˙in sin 2φr ∂r∂φr v 2v rr˙ 2 ∂2v = − 3 θ˙in cos2 φr + φ˙ 2r ∂r∂ r˙ v ∂2v r r3 2 = 2 θ˙in cos2 φr − 3 θ˙in cos2 φr + φ˙ 2r θ˙in cos2 φr v v ∂r∂ θ˙in 3 2 r ∂ v r 2 cos2 φr + φ˙ 2r φ˙ r = 2 φ˙ r − 3 θ˙in v v ∂r∂ φ˙ r 2 2 ∂ v ∂ v ∂2v ∂2v ∂2v ∂2v = 2 = = =0 = = ∂θin ∂r ∂ θin ∂θin ∂φr ∂θin ∂ r˙ ∂θin ∂ θ˙in ∂θin ∂ φ˙ r ∂2v r 2 r3 2 2 = − θ˙in sin 2φr + 3 θ˙in cos2 φr + φ˙ 2r θ˙in sin 2φr ∂φr ∂r v 2v ∂2v =0 ∂φr ∂θin
(132)
(133)
(134)
(135) (136) (137) (138) (139) (140) (141) (142) (143) (144) (145) (146) (147)
34
ACCEPTED MANUSCRIPT
(148)
EP
TE D
M AN U
SC
RI PT
∂2v r2 ˙2 r4 ˙4 = −2 θ cos 2φ − θ sin2 2φr r ∂ 2 φr v in 4v 3 in ∂2v r2 r˙ 2 = 3 θ˙in sin 2φr ∂φr ∂ r˙ 2v r2 r4 3 ∂2v sin 2φr cos2 φr = − θ˙in sin 2φr + 3 θ˙in v 2v ∂φr ∂ θ˙in ∂2v r4 2 ˙ φr sin 2φr = 3 θ˙in 2v ∂φr ∂ φ˙ r ∂2v rr˙ 2 cos2 φr + φ˙ 2r = − 3 θ˙in ∂ r∂r ˙ v ∂2v =0 ∂ r∂θ ˙ in ∂2v r2 r˙ 2 sin 2φr = 3 θ˙in ∂ r∂φ ˙ r 2v ∂2v 1 r˙ 2 = − 3 2 ∂ r˙ v v ∂2v r2 r˙ = − 3 θ˙in cos2 φr v ∂ r∂ ˙ θ˙in 2 ∂ v r2 r˙ = − 3 φ˙ r v ∂ r∂ ˙ φ˙ r ∂2v r˙ r3 2 = 2 θin cos2 φr − 3 θ˙in cos2 φr + φ˙ 2r θ˙in cos2 φr v v ∂ θ˙in ∂r 2 ∂ v =0 ∂ θ˙in ∂θin r4 3 r2 ∂2v sin 2φr cos2 φr = − θ˙in sin 2φr + 3 θ˙in v 2v ∂ θ˙in ∂φr
AC C
∂2v r2 r˙ = − 3 θ˙in cos2 φr v ∂ θ˙in ∂ r˙ 2 2 r ∂ v r4 2 = cos2 φr − 3 θ˙in cos4 φr v v ∂ 2 θ˙in r4 ∂2v = − 3 θ˙in φ˙ r cos2 φr v ∂ θ˙in ∂ φ˙ r ∂2v r ˙ r3 2 = 2 φr − 3 θ˙in cos2 φr + φ˙ 2r φ˙ r v v ∂ φ˙ r ∂r
(149) (150) (151)
(152) (153) (154) (155) (156) (157) (158) (159) (160) (161) (162) (163) (164)
∂2v =0 ∂ φ˙ r ∂θin
(165)
∂2v r4 2 ˙ = 3 θ˙in φr sin 2φr 2v ∂ φ˙ r ∂φr
(166) (167) 35
ACCEPTED MANUSCRIPT
∂2v r2 = − 3 r˙ φ˙ r v ∂ φ˙ r ∂ r˙ r4 ∂2v = − 3 θ˙in φ˙ r cos2 φr v ∂ φ˙ r ∂ θ˙in ∂2v r2 r4 = − 3 φ˙ 2r v v ∂ 2 φ˙ r
(169) (170)
SC
References
RI PT
(168)
[1] J. Sullivan, S. Grimberg, S. DAmico, Comprehensive survey and assessment
235
M AN U
of spacecraft relative motion dynamics models, Journal of Guidance Control & Dynamics 40 (8) (2017) 1–23. doi:10.2514/1.G002309. [2] G. W. Hill, Researches in the lunar theory, American Journal of Mathematics 1 (2) (1878) 5–26. doi:10.2307/2369430.
[3] W. H. Clohessy, Terminal guidance system for satellite rendezvous,
240
TE D
Aerospace Sci 29 (9) (1960) 79–86. doi:10.2514/8.8704. [4] D. F. Lawden, Optimal Trajectories for Space Navigation, Butterworths, 1963.
[5] J. Tschauner, P. Hempel, Rendezvous with a target in an elliptical orbit,
EP
Acta Astronautica 11 (2) (1965) 104–109. [6] K. T. Alfriend, S. R. Vadali, P. Gurfil, J. How, L. S. Breger, Spacecraft formation flying: Dynamics, control and navigation, Butterworths, 2010.
AC C
245
[7] T. F. Berreen, J. D. C. Crisp, An exact and a new first-order solution for the relative trajectories of a probe ejected from a space station, Celestial Mechanics 13 (1) (1976) 75–88. doi:10.1007/BF01228535.
[8] D. K. Geller, T. A. Lovell, Relative orbital motion and angles-only relative
250
state observability in cylindrical coordinates, in: Proceedings of the AASAIAA Space Flight Mechanics Meeting, Santa Fe, New Mexico, 2014.
36
ACCEPTED MANUSCRIPT
[9] F. D. Bruijn, E. Gill, J. P. How, Comparative analysis of cartesian and
RI PT
curvilinear clohessy-wiltshire equations, Journal of Aerospace Engineering, Sciences and Applications 3 (2) (2011) 1–15. doi:10.7446/jaesa.0302. 255
01.
[10] C. D. Karlgaard, F. H. Lutze, Second-order relative motion equations, Journal of Guidance Control & Dynamics 26 (1) (2001) 41–49.
SC
10.2514/2.5013.
doi:
[11] C. Bombardelli, J. L. Gonzalo, J. Roa, Compact solution of circular orbit relative motion in curvilinear coordinates, in: Astrodynamics Specialist Conference, 2016.
M AN U
260
[12] A. C. Perez, Applications of relative motion models using curvilinear coordinate frames, Ph.D. thesis (2017).
[13] K. Yamanaka, F. Ankersen, New state transition matrix for relative motion 265
on an arbitrary elliptical orbit, Journal of Guidance Control & Dynamics
TE D
25 (1) (2002) 60–66. doi:10.2514/2.4875.
[14] J. S. Ardaens, S. D’Amico, A. Cropp, Gps-based relative navigation for the proba-3 formation flying mission , Acta Astronautica 91 (2013) 341–355. doi:10.1016/j.actaastro.2013.06.025. [15] J. L. Burch, T. E. Moore, R. B. Torbert, B. L. Giles, Magnetospheric
EP
270
multiscale overview and science objectives, Space Science Reviews 199 (1-
AC C
4) (2016) 5–21. doi:10.1007/s11214-015-0164-9. [16] A. W. Koenig, T. Guffanti, S. DAmico, New state transition matrices for spacecraft relative motion in perturbed orbits, Journal of Guidance Control
275
& Dynamics 40 (7) (2017) 1749–1768. doi:10.2514/1.G002409.
[17] T. Carter, M. Humi, Fuel-optimal rendezvous near a point in general keplerian orbit, Journal of Guidance Control & Dynamics 10 (6) (1987) 567–573. doi:10.2514/3.20257.
37
ACCEPTED MANUSCRIPT
[18] T. E. Carter, New form for the optimal rendezvous equations near a keplerian orbit, Journal of Guidance Control & Dynamics 13 (1) (1990) 183–186. doi:10.2514/3.20533.
RI PT
280
[19] T. E. Carter, State transition matrices for terminal rendezvous studies:
Brief survey and new example, Journal of Guidance Control & Dynamics
285
SC
21 (1) (2015) 148–155. doi:10.2514/2.4211.
[20] Z. Dang, Solutions of tschaunerhempel equations, Journal of Guidance Control & Dynamics 40 (11) (2017) 2953–2957. doi:10.2514/1.G002774.
M AN U
[21] P. Sengupta, S. R. Vadali, Relative motion and the geometry of formations in keplerian elliptic orbits with arbitrary eccentricity, Journal of Guidance Control & Dynamics 30 (4) (2007) 953–964. doi:10.2514/1.25941. 290
[22] F. Jiang, J. Li, H. Baoyin, Approximate analysis for relative motion of satellite formation flying in elliptical orbits, Celestial Mechanics & Dynamical
AC C
EP
TE D
Astronomy 98 (1) (2007) 31–66. doi:10.1007/s10569-007-9067-8.
38