Perturbed Orbital Dynamics

C H A P T E R 5 Perturbed Orbital Dynamics 5.1 OBJECTIVES In Chapter 3, we have studied the motion of a small body treated as a point mass and subjec...

Download PDF

2MB Sizes 0 Downloads 59 Views

Report

PDF Reader
Full Text

C H A P T E R

5 Perturbed Orbital Dynamics 5.1 OBJECTIVES In Chapter 3, we have studied the motion of a small body treated as a point mass and subject to the gravitational force generated by a larger-mass body. The basic assumption of the study was that the interaction between the two bodies only consisted of the central gravitational force of the larger body, which implied that the small-mass body moved along ideal orbits called the Keplerian orbits. In the present chapter, we will study the body motion in situations where the body is not only subjected to central gravitational forces but also to other forces that we call perturbing forces. For example, the actuating force delivered by the thrusters of a spacecraft represents a perturbation with respect to a Keplerian motion. Other examples, already discussed in Section 4.2, are the perturbations due to a nonperfect spherical shape of the larger body (the gravitational anomalies), or the perturbations caused by the gravitational attraction of third bodies. Further examples include perturbations caused by nonconservative forces such as the atmospheric drag and the pressure of the solar radiation. The atmospheric drag is relevant for low-Earth orbits, since it causes a progressive reduction of the orbit semimajor axis and, consequently, earlier re-entry and burning in the Earth’s atmosphere. The motivation for using the term perturbation is that these kinds of forces are signiﬁcantly weaker than the central gravitational force of the larger body. Orbits subjected to perturbations are known as non-Keplerian orbits. In Section 5.2, we will recall the restricted two-body problem in its general form and discuss how the relevant state equations can be used to predict and simulate the behavior of a body immersed in the ﬁeld of a central gravitational force and subjected to perturbing forces. In Section 5.3, we will introduce the Gauss and Lagrange equations, which express the perturbed motion in terms of the orbital elements. In Section 5.5, we will again derive the HilleClohessye Wiltshire (HCW) equation (see Section 3.6), which describes the relative motion of a deputy spacecraft in a perturbed orbit with respect to a given nominal/chief orbit. In the same section, a feedback control strategy is derived, which is capable of bounding the magnitude of the relative motion and, if necessary, of progressively reducing the motion close to zero. These maneuvers can be performed by a deputy spacecraft which either makes a rendezvous with the chief satellite or tracks a reference orbit. In Section 5.6, the classical restricted three-body problem will be addressed, where a spacecraft is subject to the gravitational attraction of two

Spacecraft Dynamics and Control https://doi.org/10.1016/B978-0-08-100700-6.00005-2

187

Copyright © 2018 Elsevier Ltd. All rights reserved.

188

5. PERTURBED ORBITAL DYNAMICS

signiﬁcantly larger bodies, such as the Earth and the Moon. The classical Lagrangian equilibrium points (collinear and Trojan) will be derived, their stability discussed and two classes of orbits around the collinear points, Lissajous and Halo, will be brieﬂy outlined.

5.2 PERTURBED ORBITS Following the formulation of Section 3.2, let us consider two points P0 and P1 of constant masses m0 and m1, with their positions denoted by ! r 0 and ! r 1 . As discussed in Section 3.2, the center of mass (CoM) of two-body system is the natural origin O of the trajectory inertial frame. In the case of a satellite (located at P1) orbiting around a planet (located at P0), O can be assumed to coincide with the planet CoM, which implies ! r 0 ¼ 0. The motion of the body with mass m1 is described by the restricted two-body Eq. (3.5) of Section 3.2.1, which holds: ! F ðtÞ ! €r ðtÞ ¼ m ! ; r ðtÞ þ 3 m rðtÞ

! r ð0Þ ¼ ! r 0;

! r_ ð0Þ ¼ ! r_ 0

(5.1)

! where F is the resultant of all the perturbations acting on the body, m ¼ m0 ¼ Gm0, and we have used the notations ! r ¼ ! r 1 and m ¼ m1.

5.2.1 Cowell’s Method In the case of generic perturbations, a simple and general method for determining the trajectory of a body, whose motion is described by Eq. (5.1), is the numerical integration of the perturbed acceleration in Eq. (5.1). The method is known as the Cowell’s method.

5.2.2 Encke’s Method A variant of the Cowell’s method, which is more efﬁcient from a computational point of view, is the Encke’s method. The basic idea is to replace the numerical integration of the € €r with the integration of the deviation ! perturbed acceleration ! d of the perturbed accelera-

€ €r from the acceleration ! tion ! R of a reference Keplerian orbit, known as the osculating orbit (from Latin osculari, to kiss, from osculum, kiss, literally “little mouth”). Deﬁnition 1

r , which is Consider a body with mass m ¼ m1, traveling on a perturbed orbit of radius ! ! the solution of Eq. (5.1). At any given time t, the osculating orbit is deﬁned as the solution R in the time interval [t,N) of the restricted two-body equation m ! ! € R ð sÞ ¼ R ðsÞ; 3 R ð sÞ

! R ðs ¼ tÞ ¼ ! r ðtÞ;

! _ R ðs ¼ tÞ ¼ ! r_ ðtÞ; s˛½t; NÞ:

(5.2)

According to this deﬁnition, at any given time t, the osculating orbit is the Keplerian orbit deﬁned by the state vector x(t) ¼ [r(t),v(t)] of the perturbed orbit at time t. Let us observe

5.2 PERTURBED ORBITS

189

that the osculating orbit changes in time but, at any time t, perturbed and osculating orbits are in contact. We say that they osculate or, in other words, they “kiss each other.” , The deﬁnition of the osculating orbit is rather delicate. In fact, let us consider at time t two orbits which are denoted by the subscripts 1 and 2 and are deﬁned by their state vectors x1(t) ¼ [r1,v1] and x2(t) ¼ [r2,v2], that is, by the inertial coordinates of their position and velocity. The orbits osculate each other if and only if x1(t) ¼ x2(t). So, which is the difference between the two orbits? They only differ in their accelerations, since a1(t) s a2(t). More speciﬁcally, we may assume that the perturbed orbit denoted by subscript 2 has a nonKeplerian acceleration, that is a2 ðtÞs mr2 =jr2 j3 ¼ mr1 =jr1 j3 , whereas the osculating orbit, denoted by subscript 1, has a Keplerian acceleration, that is a1 ðtÞ ¼ mr1 =jr1 j3 . Moreover, two Keplerian orbits with different state at the same time, i.e., x1(t) s x2(t), are not osculating also if a1(t) ¼ a2(t). This typically occurs if r1(t) ¼ r2(t), but v1(t) s v2(t), as in the orbits of the Newton’s cannonball (see Exercise 7 of Section 3.5.1), where only the magnitude of the initial velocity is different and so is the orbit eccentricity. Fig. 5.1 sketches an osculating orbit that at time t is perturbed by the drag acceleration D! a 2 ðtÞf ! v ðtÞ, thus causing the orbit to decay. By deﬁning the perturbed orbit deviation as ! ! d ¼ ! r R;

(5.3)

from Eqs. (5.1) and (5.2), we can write the equation of the relative motion as follows: ! m ! ! m ! F ! € d ¼ 3 r d 3 r þ R r m

(5.4)

or, equivalently, as m ! € d ¼ 3 R

FIGURE 5.1

! R3 F ! ! 1 r þ d þ . m r3

Sketch of an osculating orbit and of the perturbed orbit due to drag.

(5.5)

190

5. PERTURBED ORBITAL DYNAMICS

Since R3/r3 1 is the difference between two almost equal numbers, the integration of Eq. (5.5) may become dominated by numerical errors. The problem can be circumvented by deﬁning the scalar quantity ! ! d $ d 2! r q ¼ ; (5.6) r2 where now the subtraction is between quantities of different magnitude. It can be easily seen that q ¼ R2/r2 1, which identity allows us to write R3/r3 1 as R3 3 þ 3q þ q2 1 ¼ q . r3 1 þ ð1 þ qÞ3=2

(5.7)

In this way, R3/r31 can be computed through Eqs. (5.6) and (5.7), without subtracting close quantities. The Encke’s method can be summarized as follows. Consider the time interval [t,t þ Dt], with Dt sufﬁciently small. In this interval, the following operations are carried out: ! 1. R ðsÞ with s˛½t; t þ Dt is computed by the numerical integration of Eq. (5.2). 2. ! r ðsÞ ¼ ! r ðtÞ with s˛½t; t þ Dt is assumed constant. ! 3. d ðsÞ with s˛½t; t þ Dt is computed by the numerical integration of Eq. (5.5), with the help ! ! _ of Eqs. (5.6) and (5.7), and starting from the initial conditions d ðtÞ ¼ d ðtÞ ¼ 0. 4. The new values of the perturbed position and velocity are set equal to ! r ðt þ DtÞ ¼ ! ! ! _ ! _ R ðt þ DtÞ þ d ðt þ DtÞ and to ! r_ ðt þ DtÞ ¼ R ðt þ DtÞ þ d ðt þ DtÞ. 5. Starting from the time interval [0,Dt], these operations are iterated along the sequence of time intervals {[Dt,2Dt],...,[iDt,(i þ 1)Dt],...,[(N 1)Dt,NDt]}. The sequence of computations provides ! r ðtÞ and ! r_ ðtÞ for t ˛ [0,NDt], where NDt is the ﬁnal simulation time.

5.3 DYNAMICS OF THE ORBITAL ELEMENTS Cowell’s and Encke’s methods suggest a numerical solution of the perturbed two-body equation (5.1), where the perturbed orbits are described in terms of the position ! r ðtÞ and ! _ velocity v ðtÞ. A ﬁrst issue is that these methods do not give a direct indication on how

the perturbations affect the six orbital elements of p ¼ [U,i,u,a,e,q], which have been introduced in Section 3.4. Another issue is that they provide a numerical description of the motion, without any analytical insight except in a few particular cases. Gauss and Lagrange planetary equations, to be developed in the following sections, allow these issues to be solved, since their integration provides the time evolution of the orbital elements of the perturbed orbits.

5.3.1 Gauss Planetary Equations Consider a two-body system described by Eq. (5.1), where E ¼ P0 denotes the Earth point mass and P1 the spacecraft point mass. We will use two frames of reference deﬁned n !already ! !o in Sections 2.5 and 3.3.4: the Earth centered inertial (ECI) frame JE ¼ P0 ; j 1 ; j 2 ; j 3 and

5.3 DYNAMICS OF THE ORBITAL ELEMENTS

F3h3

h3 i

– F2 h2 – F h

1 1

F

h1 = r / r

j3 P

g

h2

Equator

j1

FIGURE 5.2

191

j2

E Ω

n Line of nodes

A Keplerian orbit and a perturbing force.

n ! ! ! ! o ! the Hill’s frame H ¼ P1 ; h 1 ¼ ! r =r; h 2 ; h 3 ¼ h =h . The ﬁrst two Hill-frame axes, h 1 , the ! radial axis, and h 2 , the longitudinal or along-track axis, are unit vectors in the plane of the ! osculating orbit deﬁned in Section 5.2, and the third axis h 3 , the cross-track axis, is the unit vector orthogonal to this plane and parallel to the angular momentum of the osculating orbit. The goal is to derive differential equations that describe how the orbital elements of ! p ¼ [U,i,u,a,e,q] change in time, in the presence of a nonnull perturbing force F in Eq. (5.1). The situation is sketched in Fig. 5.2. In the previous vector p, the argument of perigee u has been replaced by the argument of ! ! latitude u ¼ u þ q. F will be replaced by the perturbing acceleration ! a ¼ F =m and by the

Hill-frame coordinates ah ¼ [ah1,ah2,ah3]. The current parameter vector p(t) is the osculating orbit vector, and the perturbation dp(t), to be computed, is deﬁned as the increment dp(t) ¼ p(t þ dt) p(t) over the inﬁnitesimal time interval dt. In the case of the true anomaly, the increment dq(t) ¼ dqosc(t) þ dqpert(t) is written as the sum of two terms: the increment dqosc refers to the free response qosc of the osculating orbit and the increment dqpert refers to the perturbed ! anomaly qpert due to F . The force-free anomaly qosc(t) is the solution of the differential Eq. (3.43) of Section 3.3.5, which is here repeated: rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ dqosc ðtÞ h m pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃð1 þ e cos qðtÞÞ2 . ¼ 2 ¼ (5.8) 3 dt r a 1 e2 In the following, the distinction between force-free and perturbed anomaly applies also to the derivatives q_ osc and q_ pert , since the current anomaly q(t), i.e., the osculating orbit anomaly, may be thought of as the result of the past combination of free and forced responses. The same distinction applies to the argument of latitude u, which is written as u ¼ u_ pert þ q_ osc . The six perturbation equations to be derived in the next paragraphs can be arranged in the following normalized form: _ PpðtÞ ¼ AðpðtÞÞah ðtÞ; pð0Þ ¼ p0

_ di=dt; u_ pert ; a; _ e;_ q_ pert . _ pðtÞ ¼ U;

(5.9)

192

5. PERTURBED ORBITAL DYNAMICS

The matrices in Eq. (5.9) are as follows: 2 0 6 6 0 6 6 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 6 0 1 e2 6 6 AðpðtÞÞ ¼ 6 auo ð1 þ ecÞ 6 2esð1 þ ecÞ 6 6 6 6 sð1 þ ecÞ 4

sin u

0 0 0 2ð1 þ ecÞ2 cð2 þ ecÞ þ e

3

7 cos u sin i 7 7 7 7 sin u cos i 7 7 7; 7 0 7 7 7 7 0 5

cð1 þ ecÞ sð2 þ ecÞ P ¼ diag sin i; sin i; sin i; 1 e2 a; 1; e

(5.10)

0

where c ¼ cos q and s ¼ sin q. The eccentricity e has been treated as a factor of q_ pert to eliminate the ratio 1/e from A(p(t)). A similar scaling applies to {U,i,u}, since sin i is a common factor of the angle derivatives. Let us ﬁnally observe that the derivatives of the orbit shape and trajectory parameters in {a,e,q} are independent of the orbit orientation parameters in {U,i,u}. The following argumentation is a synthesis of those in Refs. [4], [8], and [11]. As necessary tools, we recall the following equations from Eqs. (3.22e24), (3.27), (3.38) and (3.40) of Chapter 3: a 1 e2 h2 m p ¼ r ¼ ¼ 1 þ e cos q 1 þ e cos q 1 þ e cos q qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ uo ¼ m=a3 ; h ¼ ma 1 e2 q_ ¼ h r2 h pe sin q he sin q he sin q : ¼ ¼ q_ ¼ 2 2 p r ð1 þ e cos qÞ a 1 e2 ð1 þ e cos qÞ 2sﬃﬃﬃﬃﬃﬃﬃﬃ 3 rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 3 m m e sin q 4 ¼ m=s5 ¼ a 1 e2 ms2

r_ ¼

rq_ osc ¼

pe sin q

2

h ¼ r

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ ma 1 e2 r

¼

(5.11)

rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ m ð1 þ e cos qÞ a 1 e2 j

From nSection 3.3.4, owe recall the transformation Rh from the Hill’s to ECI frame ! ! ! JE ¼ P0 ; j 1 ; j 2 ; j 3 , which holds: 2! 3 j1 6 ! 7 ! ! !

j Rh ¼ 4 j 2 5 h 1 h 2 h 3 ¼ ZðUÞXðiÞZðu þ qÞ ! j3 ; 2 3 cU cu sU ci su cU su sU ci cu sU si 6 7 ¼ 4 sU cu þ cU ci su sU su þ cU ci cu cU si 5 (5.12) si su si c u ci

193

5.3 DYNAMICS OF THE ORBITAL ELEMENTS

where the vectrix notation of Section 2.2.1 and the abbreviations cx ¼ cosx and sx ¼ sinx with ! x ¼ U,i,u have been employed. The notation j k of a generic inertial ECI axis, and the relevant subscript/superscript j, avoid confusion with the inclination i. The mean orbital angular rate pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ is uo ¼ m=a3 , and the orbital period is To ¼ 2p/uo. Semimajor Axis Perturbation We start from Eq. (3.87) in Section 3.5.1, which expresses the semimajor axis of the osculating orbit as a ¼ m/2E, where E is the total energy per unit mass associated with the orbit, not to be confused with the eccentric anomaly E of the Kepler’s equation. Hereafter, no mention of the eccentric anomaly will be done. It follows that a_ ¼

m E_ 2a2 _ ¼ E. 2 2E m

(5.13)

From classical mechanics, the power per unit mass holds E_ ¼ ! a $! v;

(5.14)

! where ! v ¼ ! r_ is the spacecraft velocity and ! a is the perturbing acceleration. Since ! v can be expressed in the Hill’s frame as r ¼ r h 1 , the velocity vector ! ! ! ! ! ! ! ! ! v ¼ r_ h 1 þ rd h 1 dt ¼ r_ h 1 þ rq_ osc h 3 h 1 ¼ r_ h 1 þ rq_ osc h 2 ; (5.15) and Eq. (5.13) becomes a_ ¼

2a2 ah1 r_ þ ah2 rq_ osc ; m

(5.16)

where q_ osc has been deﬁned in Eq. (5.8). By replacing rq_ osc and r_ from Eq. (5.11), the Gauss planetary equation of a is found: rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ a_ m a ¼ 2 ðe sin qah1 þ ð1 þ e cos qÞah2 Þ a að1 e2 Þ m . (5.17) pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 1 e2 ðe sin qah1 þ ð1 þ e cos qÞah2 Þ ¼ auo 1 e2 Eq. (5.17) shows that in the case of a small eccentricity, i.e., for e << 1, the major contribution ! to the semimajor axis perturbation comes from the along-track component a of ! a ¼ F =m. h2

Exercise 1 Assume in Eq. (5.17) a constant tangential perturbation ah2(t) ¼ ah0 ¼ 100 mm/s2 due to the aerodynamic drag, assume a null radial perturbation, i.e., ah1 ¼ 0, and approximate the orbital anomaly as q ¼ uot, with uo ¼ 1.17 mrad/s. Compute the initial osculating semimajor axis a0 from uo and, by assuming e ¼ e0 ¼ 0.005, solve Eq. (5.17) in terms of the perturbation da ¼ a a0. ,

194

5. PERTURBED ORBITAL DYNAMICS

Eccentricity Perturbation From Eq. (3.25) in Section 3.3.3, the magnitude of the angular momentum per unit mass of pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ the osculating orbit is given by h ¼ mað1 e2 Þ. The derivative holds: pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃ ma e e2 Þ _h ¼ mð1p ﬃﬃﬃ a_ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ e._ (5.18) 2 a 1 e2 In Section 7.2.3, the Newton’s equation of a point mass P1 rotating with respect to the origin P0 of an inertial framedhere the Earth’s CoMdtells us that the time derivative of ! _ ! the angular momentum m h is equal to the moment of force ! r F with respect to P0, which allows us to write ! _ ! ! h ¼ ! r ! a ¼ rah3 h 2 þ rah2 h 3 : (5.19) ! Since the Hill-frame expression of h is ! ! ! ! ! ! h ¼ ! r ! v ¼ r h 1 v1 h 1 þ v2 h 2 ¼ rv2 h 3 ¼ h h 3 ;

(5.20)

v in the Hill’s frame, we differentiate Eq. (5.20) to where vj, j ¼ 1,2,3 are the components of ! obtain the following alternative identity to Eq. (5.19): ! ! ! _ _ h ¼ h_ h 3 þ h h 3 .

(5.21)

By comparing Eqs. (5.19) and Eq. (5.21), we ﬁnd that h_ ¼ rah2 ;

! ! _ rah3 ! rah3 ! h2 ¼ h 1 h 3. h3 ¼ h h

The two expressions of h_ in Eqs. (5.18) and (5.22) can be combined to obtain pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃ ma e mð1 e2 Þ pﬃﬃﬃ a_ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ e_ ¼ rah2 : 2 a 1 e2

(5.22)

(5.23)

By making explicit Eq. (5.23) with respect to e,_ one ﬁnds with the help of Eq. (5.17) that sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 e2 að1 e2 Þ ah2 r 1 e2 e_ ¼ a_ pﬃﬃﬃﬃﬃ rah2 ¼ ah1 sin q þ . (5.24) 1 þ e cos q m a mae 2ae e Finally, Eq. (5.11) and a few manipulations provide the Gauss planetary equation of e: pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ cos q þ e 1 e2 _e ¼ ah2 . sin qah1 þ cos q þ (5.25) 1 þ e cos q auo The eccentricity equation is very similar to the semimajor axis Eq. (5.17), since only the ! in-plane coordinates of ! a ¼ F =m, that is the tangential and radial components, contribute

to the perturbation. The main difference is that, also for a small eccentricity e << 1, both tangential and radial components must be accounted for.

5.3 DYNAMICS OF THE ORBITAL ELEMENTS

195

Exercise 2 Assume in Eq. (5.25) a constant tangential perturbation ah2(t) ¼ ah0 ¼ 100 mm/s2 due to the aerodynamic drag, assume a null radial perturbation, i.e., ah1 ¼ 0, and approximate the orbital anomaly as q ¼ uot, uo ¼ 1.17 mrad/s. Compute the initial osculating semimajor axis a0 from uo and by assuming a ¼ a0 and e ¼ e0 ¼ 0, solve Eq. (5.25) in terms of the perturbation de ¼ e. Plot de during one orbital period and compute emax ¼ maxjdeðtÞj. , Inclination Perturbation ! ! j the Hill’s to ECI frame transformation Rh in Eq. (5.12), we ﬁnd that h 3 $ j 3 ¼ ! From ! h $ j 3 =h ¼ cos i. Time differentiation of this identity yields ! ! ! ! _ ! _ ! _ h $ j 3 h h $ j 3 h_ h $ j 3 h hhcos i di sin i ¼ ¼ ; (5.26) dt h2 h2 ! where j 3 has been treated as an inertial direction. Eqs. (5.19), (5.22) and (5.26) imply that (5.27) ! ! where the terms proportional to ah2 cancel each other and the identity h 2 $ j 3 ¼ sin i cos u in Eq. (5.12) has been used. Finally, replacement in Eq. (5.27) of the expression of r from Eq. (5.11) provides the Gauss equation of di/dt: sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ di að1 e2 Þ cos u cos u 1 e2 ¼ ah3 ¼ ah3 . (5.28) dt m 1 þ e cos q auo 1 þ e cos q ! Eq. (5.28) tells us that the inclination is only perturbed by the coordinate of ! a ¼ F =m which is orthogonal to the osculating plane. Let us assume small eccentricity, zero argument of perigee and constant semimajor axis, i.e., e << 1, u ¼ 0, and a ¼ a0. The third assumption requires that either the longitudinal perturbations are negligible or the semimajor axis is kept constant by some orbit control. Under these assumptions Eq. (5.28) simpliﬁes to di ah3 ðtÞ y cos qðtÞ. dt a 0 uo

(5.29)

Under a constant perturbation ah3 ¼ ah0 and by assuming q(t) ¼ uot, we ﬁnd the time response ah0 diðtÞ y i i0 ¼ sin uo t; (5.30) a0 u2o which proves that the inclination perturbation is close to be periodic and zero mean. Perturbation of the Right Ascension of the Ascending Node Let us consider the line of nodes ! n , which, from Fig. 5.2, can be written ! ! ! as n ¼ cosU j 1 þ sinU j 2 . With the help of Eq. (5.12) (speciﬁcally the third column of

196

5. PERTURBED ORBITAL DYNAMICS

! ! ! ! sinU j 1 cosU j 2 sin i þ cos i j 3 and ! n $ h 3 ¼ 0. In other words, the ! ! line of nodes is orthogonal to the orbit normal h 3 . Therefore, any perturbation of h 3 which destroys the orthogonality, must be compensated by a deviation of the line of nodes as the following differential equation proves:

! Rih ), we ﬁnd h 3 ¼

! d h 3 $! n ! _ ! _ n þ h 3 $! n ¼ 0: ¼ h 3 $! dt

(5.31)

An explicit expression of the two terms in Eq. (5.31) can be obtained with the help of ! _ Eq. (5.22) that provides h 3 , with the help of Eq. (5.12) which yields the scalar products between the axes of the ECI and Hill’s frames as the entries of Rih , and with the help of the ! _ ! derivative ! n_ ¼ sinU j þ cosU j U. At the end, the following identities are found: 1

2

! ! ! ! ! ! _ i h 3 $ n_ ¼ sinU h 3 $ j 1 þ cosU h 3 $ j 2 U_ ¼ Usin ! _ rah3 sin u n ¼ h 3 $! h

.

(5.32)

Replacement in Eq. (5.31) of the right-hand-side expressions in the ﬁrst and second row of Eq. (5.32), together with the expressions of r and h from Eq. (5.11), provide the Gauss planetary equation of U: sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 r sin u að1 e Þ sin u að1 e2 Þ sin u _ ¼ ah3 ¼ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ah3 ¼ ah3 U h sin i m ð1 þ e cos qÞsin i 2 ð1 þ e cos qÞsin i mað1 e Þ . pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ sin u 1 e2 ah3 ¼ uo að1 þ e cos qÞ sin i (5.33) Perturbation of the Argument of Latitude We look for the derivative u_ pert of the argument of latitude instead of the derivative u_ of the argument of perigee. Before proceeding, we decompose u_ pert in two different ways, in agreement with early considerations of this section:

! The components of j 3 in the Hill’s frame from Eq. (5.12) and the formula of the scalar triple product allow us to write 2 3 1 si su 0 ! ! ! 6 7 h 1 $ j 3 h 3 ¼ det4 0 si cu 0 5 ¼ si cu . (5.34) 0 ci 1

197

5.3 DYNAMICS OF THE ORBITAL ELEMENTS

Time differentiation of Eq. (5. 34) and di/dt from Eq. (5.27) provide the identities ! ! ! ! _ ! _ ! di h 1 $ j 3 h 3 þ h 1 $ j 3 h 3 ¼ cos i cos u sin i sin uu_ dt cos i cos2 u ¼ r ah3 sin i sin uu_ h

.

(5.35)

An alternative expression of the derivative of the right-hand-side term in Eq. (5.34) is ob! _ ! ! ! tained from the second identity in Eq. (5.22), from h 1 ¼ q_ osc h 3 h 1 ¼ q_ osc h 2 in Eq. (5.15), and from Eq. (5.12) as follows: ! _ ! ! ! ! _ ! h 1$ j 3 h 3 þ h 1$ j 3 h 3

! ! ! ! rah3 ! ! h 2 $ h 1 þ q_ osc j 3 h 3 $ h 2 ¼ j3 h ; rah3 _ cos i qosc sin i sin u ¼ h

(5.36)

where q_ osc has been deﬁned in Eq. (5.8). We now write the equality of Eq. (5.35) with Eq. (5.36). The two previous expressions of u_ pert and some manipulations lead to the intermediate equation: u_ pert ¼ u_ þ q_ pert ¼ u_ q_ osc ¼

r cos i sin u ah3 : h sin i

(5.37)

Substitution of r and h from Eq. (5.11) converts Eq. (5.37) into the Gauss planetary equation of the argument of latitude: pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ sin u cos i 1 e2 u_ pert ¼ ah3 . (5.38) auo ð1 þ e cos qÞsin i The equation of the argument of perigee can be obtained from Eq. (5.37). Perturbation of the True Anomaly We start from the second identity of the radius equation in Eq. (5.11), and by differentiating both sides with respect to time, we ﬁnd the equation: _ _ þ e cos qÞ ¼ 2hh=m. r e_ cos q eq_ sin q þ rð1

(5.39)

Eq. (5.39) can be solved for q_ as a function of r,_ e,_ and h_ as follows: q_ ¼

1 2h _ r_ _ ð1 þ e cos qÞ . h þ cos qe r e sin q mr

(5.40)

Eq. (5.11), the ﬁrst identity of Eq. (5.22), and Eq. (5.25) provide the derivatives of Eq. (5.40), that is

198

5. PERTURBED ORBITAL DYNAMICS

r_ ¼

he sin q að1 e2 Þ

. h_ ¼ rah2 h cos q þ e e_ ¼ ah1 sin q þ ah2 cos q þ m 1 þ e cos q

(5.41)

Replacement of the three identities in Eq. (5.40) and a few manipulations yield the intermediate equation: 2 2 h 2ah2 h _q ¼ h ah1 cos q þ ah2 cos qð1 þ e cos qÞ þ cos q þ e cos q þ me 1 þ e cos q me sin q r2 sin q ; h ah2 ðcos2 q 1Þð1 þ e cos qÞ þ ðcos2 q 1Þ ¼ ah1 cos q þ þ q_ osc me 1 þ e cos q sin q

(5.42)

where q_ osc ¼ h r2 is reported in Eq. (5.11). The identity q_ ¼ q_ pert þ q_ osc and further manipulations provide the Gauss planetary equation of the anomaly perturbation: eq_ pert ¼

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 þ e cos q 1 e2 sin q ; ah1 cos q ah2 1 þ e cos q auo

(5.43)

where e has been treated as a scale factor.

5.3.2 Lagrange Planetary Equations Generic Equations In the Lagrange planetary equations, unlike the Gauss equations of the previous section, ! the perturbing acceleration ! a ¼ F ðtÞ m in Eq. (5.1) is restricted to the gravity acceleration of third bodies and to the gravity anomalies of the planet in P0. Both kinds of perturbation have been treated in Sections 4.2.1 to 4.2.4. Let us decompose the total gravity potential V ! r affecting the small body in P1, which is orbiting around P0, into the central potential V0(r) ¼ m/r of P0, with r ¼ ! r , and into the perturbing potential U ! r . We prefer this notation [8] to the usual R [11], since R is mainly reserved to rotation matrices. The perturbed potential V and the gravity acceleration ! g are written as follows: m rÞ r Þ ¼ þ Uð! Vð! r Þ ¼ V0 ðrÞ þ Uð! r ; m ! g ð! r Þ ¼ VT Vð! r Þ ¼ 3! r VT Uð! rÞ r

(5.44)

where the sign conventions of Sections 3.2.2 and 4.2.1 are maintained. The inertial coordinates centered in P0, i.e., the ECI coordinates in the case of the Earth, allow the two-body perturbed equation in Eq. (5.1) to be rewritten as

199

5.3 DYNAMICS OF THE ORBITAL ELEMENTS

r_ ðtÞ ¼ v

2

vUðrÞ=vr1

3

7 6 m 7 6 _ vðtÞ ¼ 3 r 6 vUðrÞ=vr2 7; 5 4 r vUðrÞ=vr3

(5.45)

where r ¼ [r1,r2,r3] and v ¼ [v1,v2,v3] are the position and velocity coordinates of P1. Eq. (5.45) is an autonomous state equation. Deﬁnition 1 in Section 5.2 tells us that at any time t, the pair {r(t),v(t)} that satisﬁes Eq. (5.45) deﬁnes an osculating orbit, which can be univocally materialized by the vector p(r,v) ¼ [U,i,u,a,e,M] of the orbital elements. In the next derivations, the true anomaly q(t) will be replaced by the mean anomaly M(t) whose expression, from Eq. (3.49) of Section 3.3.5 is as follows: pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ MðtÞ ¼ uo t ¼ EðtÞ e sin EðtÞ; uo ¼ m=a3 ; (5.46) where E is the eccentric anomaly, and uo (denoted by n in other textbooks) is the mean orbital rate. Lagrange equations are concerned with long-period variations, which are typical of the orbital parameters that are constant for Keplerian orbits, or whose rate, always for Keplerian _ orbits, is constant like MðtÞ ¼ uo . We can write the pair {r(t),v(t)} and their time derivatives as functions of p(t) as follows: rðtÞ ¼ fðpðtÞÞ;

vðtÞ ¼ gðpðtÞÞ

r_ ðtÞ ¼

6 vfðpðtÞÞ X vfðpÞ dpk þ vt vpk dt . k¼1

_ vðtÞ ¼

6 vgðpðtÞÞ X vgðpÞ dpk þ vt vpk dt k¼1

(5.47)

The meaning of Eq. (5.47) is that the perturbed orbit represented by the pair {r(t),v(t)} is the envelope of a continuous sequence of osculating orbits represented by the vector p(t). Taking the equality of the position and velocity derivatives in Eqs. (5.45) and (5.47) provides the following identities: vfðpðtÞÞ ¼ vðtÞ vt 6 6 X X vfðpÞ vr vr _ pðtÞ ¼ 0 p_k ðtÞ ¼ p_k ðtÞ ¼ vp vp vp k k k¼1 k¼1 vgðpðtÞÞ m ¼ 3r vt r 6 6 X X vgðpÞ v_r v_r _ pðtÞ ¼ VT UðrÞ p_k ðtÞ ¼ p_k ðtÞ ¼ vp vp vp k k k¼1 k¼1

;

(5.48)

200

5. PERTURBED ORBITAL DYNAMICS

where the matrices

2

vr1 6 6 vp1 6 6 vr vr 6 2 ¼ 6 6 vp1 vp 6 6 vr 4 3 vp1

. . .

3 vr1 7 vp6 7 7 vr2 7 7 7; vp6 7 7 vr3 7 5 vp6

2

vv1 6 6 vp1 6 6 vv v_r vv 6 2 ¼ ¼ 6 6 vp1 vp vp 6 6 vv 4 3 vp1

. . .

3 vv1 7 vp6 7 7 vv2 7 7 7 vp6 7 7 vv3 7 5 vp6

(5.49)

are the 3 6 Jacobian matrices of the vectorial functions r ¼ f(p) and v ¼ g(p). Since the perturbing potential U(r) in Eq. (5.45) depends on r(p), the row vector gradient with respect to p is obtained by premultiplying the 3 1 column vector VTU(r) by the 6 3 transposed Jacobian matrix (vr(p)/vp)T as follows: 3 2 vU ðrðpÞÞ 6 vU 7 T T 7 6 vU ðrðpÞÞ vrðpÞ 7 6 T « (5.50) ¼ V U ðrÞ ¼ 6 7: 7 6 vp vp 4 vU ðrðpÞÞ 5 vM The implicit version of the Lagrange planetary equations is obtained in two steps. First, the last identity of the fourth row in Eq. (5.48) is premultiplied by (vr(p)/vp)T as in Eq. (5.50). _ Then, the zero identity ðv_r=dpÞT ðvr=vpÞpðtÞ ¼ 0 from the second row of Eq. (5.48) is subtracted, which yields T ! T T vr vv vv vr vUðrðpÞÞ _ _ pðtÞ ¼ ; (5.51) LpðtÞ ¼ vp vp vp vp vp where L is a 6 6 matrix to be explicitly computed and inverted. Lagrangian Brackets The generic element Lij of the matrix L in Eq. (5.51) has the following form, known as the Lagrangian bracket, and often denoted by [pi,pj]: ! ! T T vr vv vv vr Lij ¼ ½pi ; pj ¼ . (5.52) vpi vpj vpi vpj Some properties of the Lagrangian brackets, which are easily checked, are the following: 1. Lii ¼ [pi,pi] ¼ 0, 2. Lij ¼ [pi,pj] ¼ Lji ¼ [pj,pi], 3. Lij ¼ [pi,pj] does not explicitly depend since r and v only depend on p as on time t expressed by Eq. (5.47), that is vLij vt ¼ v pi ; pj vt ¼ 0.

201

5.3 DYNAMICS OF THE ORBITAL ELEMENTS

Exercise 3

Prove the previous third property by explicitly computing v pi ; pj vt. , As a consequence of these properties, L is a skew symmetric matrix (the main diagonal is zero). Thus, we just need to compute the N ¼ 6 5/2 ¼ 15 upper diagonal brackets Lij ¼ [pi,pj], i ¼ 1,., 6; j ¼ i þ 1,.,6. They can be calculated by knowing r(p) and v(p) at any point of the osculating orbit, for instance in the periapsis neighborhood. In the sequel, we decompose p into angular and trajectory parameters, namely as p ¼ [a,b], with a ¼ [U,i,u] (angular) and b ¼ [a,e,M] (trajectory), and we show that r and v can be j written as the product between the perifocal-to-ECI frame transformation Rp ðaÞ, function of the angular parameters of a, and the perifocal coordinates, functions of the trajectory parameters of the vector b. Since we work in the neighborhood of the periapsis, the inertial vector r is obtained from the periapsis coordinates in Eq. (3.45) of Section 3.3.5, as follows 3 2 að 1 e Þ þ o M 2 2 3 aðcos E eÞ 6 rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 7 6 6 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ2 7 3 7 i i 1þe 7 rðpÞ ¼ Rp ðaÞrp ðbÞ ¼ Rp ðaÞ4 a 1 e sin E 5 ¼ Rip ðaÞ6 (5.53) M þ o M 7: 6a 4 5 1e 0 0 The velocity vector v is obtained in the same way and holds: 3 2 rﬃﬃﬃﬃ 3 m M 3 2 6 a ð1 eÞ2 þ o M 7 7 6 aE_ sin E 7 6 sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 7 6 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 7 6 i i i 7 vðpÞ ¼ Rp ðaÞvp ðbÞ ¼ Rp ðaÞ6 2 7. (5.54) m ð 1 þ e Þ 4 a 1 e2 E_ cos E 5 ¼ Rp ðaÞ6 7 6 þo M 7 6 að1 eÞ 5 4 0 0 In Eqs. (5.53) and (5.54) the periapsis approximations sin E ¼ E þ o E3 ; cos E ¼ 1 þ o E2 ; M ¼ Eð1 eÞ þ o E3 rﬃﬃﬃﬃ uo m 1 _ þ o E2 ¼ E ¼ a 1e 1 e cos E j

j

have been used, and Rp ðaÞ is the same as Rh in Eq. (5.12) upon replacement of u with u, as follows: 2

c U c u sU c i su

6 Rjp ðaÞ ¼ ZðUÞXðiÞZðuÞ ¼ 4 sU cu þ cU ci su si su

cU su sU ci cu sU su þ cU ci cu si cu

with the abbreviations cx ¼ cos x and sx ¼ sin x for x ¼ U,i,u.

sU si

3

7 cU si 5; ci

(5.55)

202

5. PERTURBED ORBITAL DYNAMICS

The expressions in Eqs. (5.53e5.55) suggest that the Lagrange brackets can be subdivided into three groups: 1. three angular brackets [aj,ak], j,k ¼ 1,2,3: [U,i], [U,u], [i,u], 2. three trajectory brackets [bj,bk], j,k ¼ 1,2,3: [a,e], [a,M], [e,M], 3. nine mixed brackets [aj,bk], j,k ¼ 1,2,3: [U,a], [U,e], [U,M], [i,a], [i,e], [i,M], and [u,a], [u,e], [u,M]. Angular, trajectory, and mixed brackets have the following generic expressions: 0 1 T T vRip vRip vRip vRip C TB ½ai ; aj M¼0 ¼ rp @ A vp vai vaj vaj vai M¼0 T T vrp vvp vrp vvp ; ½bi ; bj M¼0 ¼ vbi vbj vbj vbi M¼0 T T i i vrp vvp p vRp p vRp ½bi ; aj M¼0 ¼ Ri vp Ri rp vbi vaj vbi vaj

(5.56)

M¼0

which must be computed at the orbital periapsis, that is for M ¼ 0. In the following, the notation h ¼ 1e2 will be used. Exercise 4 ANGULAR BRACKETS

Compute from Eq. (5.55) the derivatives vRip =vU, vRip =vi, and vRip =vu, and, using the ﬁrst row in Eq. (5.56), prove the following identities: pﬃﬃﬃﬃﬃﬃﬃﬃ L12 ¼ ½U; i ¼ hma sin i; L13 ¼ ½U; u ¼ 0; L23 ¼ ½i; u ¼ 0: , (5.57) Exercise 5 TRAJECTORY BRACKETS

Compute from Eqs. (5.53) and (5.54) the Jacobian matrices vrp/vb and vvp/vb, and using the second row in Eq. (5.56), prove the following identities: 1 pﬃﬃﬃﬃﬃ ma; L56 ¼ ½e; M ¼ 0: , L45 ¼ ½a; e ¼ 0; L46 ¼ ½a; M ¼ (5.58) 2a Exercise 6 MIXED BRACKETS

Using the third row in Eq. (5.56) and the derivative matrices computed in Exercise 4 and 5, prove the following identities: pﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃ hma hma cos i; L24 ¼ ½i; a ¼ 0; L34 ¼ ½u; a ¼ L14 ¼ ½U; a ¼ 2a 2a rﬃﬃﬃﬃﬃ rﬃﬃﬃﬃﬃ ma ma ., (5.59) cos i; L25 ¼ ½i; e ¼ 0; L35 ¼ ½u; e ¼ e L15 ¼ ½U; e ¼ e h h L16 ¼ ½U; M ¼ 0;

L26 ¼ ½i; M ¼ 0;

L36 ¼ ½u; M ¼ 0

203

5.3 DYNAMICS OF THE ORBITAL ELEMENTS

The brackets [pi,pj] in Eqs. (5.57e5.59) and their reciprocal [pj,pi] ¼ [pi,pj] are the components of the skew-symmetric matrix L in Eq. (5.51), which takes the form: 2

3

6 0 1 0 6 6 6 6 6 1 0 0 6 6 6 6 0 0 6 0 6 6 pﬃﬃﬃﬃﬃﬃﬃﬃ6 L ¼ E hma6 1 6 cot i 0 6 6 2a 2a 6 6 6 1 6 cot i 6 0 6 h h 6 6 6 6 4 0 0 0 E ¼ diagðsin i; 1; 1; 1; e; 1Þ;

cot i 2a

cot i h

0

0

1 2a

1 h

0

0

0

0

1 pﬃﬃﬃ 2a h

0

7 7 7 7 7 0 7 7 7 7 7 0 7 7 7 7 1 7 7E pﬃﬃﬃ 7 ; 2a h 7 7 7 7 7 0 7 7 7 7 7 7 0 5 0

(5.60)

h ¼ 1 e2

where the scale factor matrix E includes the eccentricity e and sini, as they may tend to zero. pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃ We recall from Eq. (3.25) of Section 3.3.3 that h ¼ mað1 e2 Þ ¼ hma is the orbit angular momentum per unit mass. The inversion of the matrix L in Eq. (5.60) provides the explicit version of the Lagrange planetary equations [11]: rﬃﬃﬃﬃ T vU m 1 _ EpðtÞ ¼ AðpÞE þB 3 vp a 2 0 1 0 6 6 1 0 cot i 6 6 6 0 cot i 0 1 6 6 AðpÞ ¼ pﬃﬃﬃﬃﬃﬃﬃﬃ 6 ham 6 0 0 0 6 6 6 6 0 0 h 4 0 p ¼ ½U; i; u; a; e; M;

0

0

0

0

0

0

0

h

0

0

0 0 pﬃﬃﬃ pﬃﬃﬃ 2a h h h

E ¼ diagðsin i; 1; 1; 1; e; 1Þ

0

3

7 7 7 7 7 0 7 7 ; pﬃﬃﬃ 7 2a h 7 7 7 pﬃﬃﬃ 7 h h7 5 0

0

2 3 0 6 7 607 6 7 6 7 6 7 607. 6 7 B ¼ 6 7 607 6 7 6 7 6 7 607 4 5 1

(5.61)

204

5. PERTURBED ORBITAL DYNAMICS

The inverse L1 ¼ E1AE1 is related to the matrix P of the Poisson brackets of the Hamiltonian dynamics, that is PT ¼ L1 [3]. The scale factor (ma)1/2 can be expressed in terms of the mean orbital rate uo as follows: 1 1 h s i . pﬃﬃﬃﬃﬃ ¼ ma uo a 2 m 2

(5.62)

Average Perturbations Due to the Earth’s Flatness Comparison of Eq. (5.44) with Eq. (4.14) of Section 4.2.2 allows us to write the perturbing potential U due to the Earth’s ﬂatness, that is the J2 zonal term, as follows: 2 m 3J2 Re 1 r3 2 UðrðpÞÞ ¼ UðpÞ ¼ r 2 3 r r ; 2 m 3J2 Re a 3 1 2 1 cos 2ðu þ qÞ sin i ¼ a 2 r 3 2 a

(5.63)

1 where r ¼ jrj, r ¼ ½r1 ; r2 ; r3 , r3 ¼ r sin i sin u, u ¼ u þ q, and a=r ¼ 1 e2 ð1 þ e cos qÞ. As already pointed out in Section 4.2.2, the zonal potential in Eq. (5.63) includes shortperiod and long-period variable parameters, where short-period means periodic variations during the orbital period such as those due to the true anomaly q, whereas long-period refers to variations of orbital parameters such as a, e, i and u that are constant in a Keplerian orbit. As explained in Ref. [7], long-period variations only refer to u, whereas the variations of a, e, and i are known as secular variations. Since we are concerned with long-period variations, we take the average of the short-period terms. Let us consider ﬁrst pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ the average of (a/r)3, which exploits the differential identity dq ¼ ða=rÞ2 1 e2 dM in Eq. (3.50) of Section 3.3.5: Z 2p 3 Z a3 uo 2p=uo a3 1 a ¼ dt ¼ dM r r 2p 0 r 2p 0 ; Z 2p 1 2 3=2 ð1 þ e cos qÞdq ¼ 1 e ¼ 3=2 0 2p 1 e2

(5.64)

where the overbar denotes average (see Refs. [7] and [11]). With a similar development, we ﬁnd the identities a3 r

cos 2 q ¼

a3 r

sin 2 q ¼ 0;

(5.65)

which together with Eq. (5.64) gives the following expression of the average potential: UðpÞ ¼

m að1 e2 Þ3=2

2 3J2 Re 1 sin2 i . 3 2 2 a

(5.66)

205

5.3 DYNAMICS OF THE ORBITAL ELEMENTS

Since UðpÞ only depends on the orbital elements in {i,a,e}, and the elements of the average

triad i; a; e are constant, the average of the Lagrange planetary equation in Eq. (5.61) becomes: 2 3 2

:

6 U: 6 i 6 : 6 _pðtÞ ¼ 6 u: 6 a 6 : 6 4 e: M

cosi 6 7 6 7 6 7 0 6 7 6 7 7 5 2 6 7 7 6 7 7 2 2 sin i 6 7 7 3J R u 2 2 e o 6 7ðtÞ ¼ 7; 7 7 2 26 2 a 0 7 7 1e 6 6 7 7 6 5 7 0 6 7 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 7 6 6 7 3 4 1 sin2 i 1 e2 5 2 3

pð0Þ ¼ p0 ;

(5.67)

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ where uo ¼ m=a3 . The ﬁrst-row equation of (5.67) has been already used in Section 3.5.2 to ﬁnd the inclination of Sun-synchronous orbits. The third-row equation has been used in Section 3.5.2 to ﬁnd the inclination of the highly eccentric Molnyia orbits, as they need a constant argument of perigee. Exercise 7 Derive Eq. (5.67) through the negative gradient of UðpÞ in Eq. (5.66). ,

5.3.3 Frozen Orbits The concept of frozen orbits dates back to the early times of artiﬁcial satellites. Deﬁnition 2 A frozen orbit is such that one or more mean orbital elements ideally remain constant under the inﬂuence of the central body (the Earth in this context) gravity anomalies. The concept can be extended to other perturbations. , We are interested in keeping the mean orbital eccentricity constant. We have seen from Eq. (5.67) that e_ ¼ 0 under the action of the Earth’s ﬂatness. We have also mentioned in Section 4.2.3 that the weak pear-shaped ﬁgure of the Earth was conﬁrmed by the eccentricity time variation of the second US satellite Vanguard I. These two results suggest of studying the mean Lagrange equations under the J2 and J3 perturbations. To this end, we complete Eq. (5.63) with the J3 potential term from Section 4.2.3, using the notation u ¼ u þ q: 2 3 m 3J2 Re 1 r3 2 m 5J3 Re r3 3 r3 2 UðrÞ ¼ r 2 3 r 2 r r r r 5 r 2 m 3J2 Re 1 1 cos 2 u a3 . (5.68) sin2 i ¼ a 2 3 2 r a 3 m 5J3 Re 3 3 2 1 a 4 sin i sin u þ sin2 i sin 3 u sin i a 2 5 4 4 r a

206

5. PERTURBED ORBITAL DYNAMICS

The average of the J2 potential can be found in Eq. (5.66). The average of the J3 potential needs the following averages [12]: a4 5=2 cos q ¼ e 1 e2 r ; (5.69) a4 a4 a4 sin q ¼ cos 3 q ¼ sin 3 q ¼ 0 r r r which have been computed by the same method of Eq. (5.64). The resulting average potential is found to be: ! 2 m Re J2 3J3 e Re 2 2 UðpÞy 2 3 sin i þ sin i 4 5 sin i sin u . 3=2 5=2 a a a 8ð1 e2 Þ 4ð1 e2 Þ (5.70) The key equations in Eq. (5.61), which must be converted into average equations, are the third one, that of the argument of perigee u, and the ﬁfth one, that of the eccentricity e. They are repeated below with the approximation h ¼ 1 e2 y 1, since we anticipate that for frozen near-polar orbitsdwe are interested in these orbitsdthe mean constant eccentricity must be very small, close to e ¼ 0.001, which implies that in the following derivation we can assume h y 1. The average equations are as follows: de 1 vU vU y pﬃﬃﬃﬃﬃ dt e am vu vM (5.71) : du 1 cos i vU 1 vU y pﬃﬃﬃﬃﬃ dt e ve am sin i vi The eccentricity derivative only involves the J3 component in Eq. (5.70). By setting h y 1 and marking the averages with an overbar, the derivative of the average eccentricity is given by [11]: rﬃﬃﬃﬃ 3 de m Re 3J3 y sini 4 5 sin2 i cosu. (5.72) 3 dt a 8 a Eq. (5.72) tells us that e can be made constant in three modes: (1) by making the orbit equapﬃﬃﬃﬃﬃﬃﬃﬃ torial, that is by i ¼ 0; (2) by imposing the average inclination i ¼ sin1 4=5, which is the same condition for freezing u in the case of the Molnyia orbits in Section 3.5.2; (3) by imposing cosu ¼ 0. We are only interested in the last condition as it applies to a generic average inclination i, the other two conditions deﬁning singular orbits. The generic condition is summarized as follows: de p ¼ 05cosu ¼ 00u ¼ ; dt 2

(5.73)

where the sign of the solution is still to be found. The derivative of the argument of perigee in Eq. (5.71) involves both J2 and J3 terms, and, by neglecting e2 with respect to the unit just after the derivation of UðpÞ, the derivative holds:

5.4 FROM N-BODY SYSTEM TO THREE-BODY SYSTEM

du y dt

rﬃﬃﬃﬃ 2 m Re 3J2 5 cos2 i 1 3 4 a a . rﬃﬃﬃﬃ 3 1 m 3J3 Re 2 2 2 2 2 i 5 cos i 1 e cos i 15 cos i 11 sinu þ sin a esini a3 8

207

(5.74)

If we restrict our choice to a near-polar orbit deﬁned by i ¼ p 2 þ Di and Di << 1 rad, and if we neglect second-order terms in Di and eDi, the identity du=dt ¼ 0 in Eq. (5.74) provides the stationary condition J3 R e ey (5.75) sini sinu. 2J2 a ! Since for near polar orbits we have sin i > 0 and, from Section 4.2.1, we have J3 =J2 > 0, Eq. (5.75) implies e > 0 and rules out u ¼ p=2 from Eq. (5.73). In conclusion, the mean eccentricity stationary conditions are: J3 R e u ¼ p=2; ey (5.76) siniz0:001: 2J2 a More accurate conditions should involve higher-order harmonics. Exercise 8 Compute from Eq. (5.75) the stationary e for the Sun-synchronous satellite GOCE at an orbit altitude h ¼ 260 km. , Exercise 9 Let us denote the mean perturbations of e and u, about the mean equilibrium deﬁned by Eq. (5.76), with De and Du, respectively. Derive the linear time invariant, second-order perturbation state equation of (5.72) and (5.74) for a polar orbit deﬁned by i ¼ p 2. Prove that the imaginary eigenvalues of the state equation are given by rﬃﬃﬃﬃ 2 m Re 3J2 :, (5.77) l1;2 ¼ j 3 4 a a

5.4 FROM N-BODY SYSTEM TO THREE-BODY SYSTEM Following the formulation of Section 4.2.5, let us consider n point masses Pi, i ¼ 0,...,n 1 of constant mass mi, with inertial position ! r i . The Newton’s equation of each point mass is ! n1 Fi ! €r ¼ X Gmj ! ! r þ r ; rji ¼ ! rj! r i ; i j i 3 mi rji jsi

(5.78)

! where F i is the total nongravitational force acting on the ith point mass. Eq. (5.78) describes the motion of n bodies with mass mi, which are interacting through central gravitational forces, and is known as the n-body system equation. The n-body CoM is deﬁned by

208

5. PERTURBED ORBITAL DYNAMICS n1 X mi ! ! r i; rc ¼ m i¼0

m ¼

n1 X

mi .

(5.79)

i¼0

As already discussed in Section 4.2.5, the n-body CoM is the natural origin O of the inertial frame associated to the n-body system. In the case of one or more satellites located at Pi, which are orbiting around the Earth centered in E ¼ P0, the origin O can be assumed to coincide with the Earth’s CoM, which implies ! r 0 ¼ 0. For n ¼ 3, the Newton’s equation of each point mass is given by ! Gmk ! Fi ! €r ¼ Gmj ! ! ! rj ri þ 3 rk ri þ ; i mi r3ji rki

(5.80)

with i,j,k ˛ {0,1,2}, i s j s k. This equation describes the motion of any of the three bodies, which interacts with the other two bodies through central gravitational forces. Eq. (5.80) is known as the three-body system equation. In this chapter, we will study two important cases covered by Eq. (5.80): 1. The ﬁrst case concerns the system made by a large body and two small bodies, for example, by the Earth and two spacecraft. The aim is to describe the relative motion between the two small bodies. The equation system to be written and studied is the HCW equation, which was already found in Section 3.6 and will be again derived and studied in Section 5.5. 2. The second case concerns the system made by a small body subject to the gravitational force of two large bodies, such as a satellite subject to the gravitational force of the Earth and the Moon. This leads to the restricted three-body equation to be studied in Section 5.6. The classical three-body system {Po ¼ S, P1 ¼ E, P2 ¼ M} is the SuneEartheMoon system, already studied by Newton and subsequently by Euler and Lagrange, until H. Poincaré in 1887 proved that no analytic solution of the three-body problem exists, which is given by algebraic expressions and integrals. These difﬁculties well explain the longstanding efforts since the Baylonian and Greek astronomers in developing a simple yet accurate mathematical model (the Lunar Theory [7]) of the complex Moon motion under Sun and Earth gravitation. The orbital motion of the third body, such as the Moon, can be described either in an inertial frame, the ecliptic frame, here denoted with the inertial frame notation n like ! ! !o I ¼ O; i 1 ; i 2 ; i 3 , or in the synodic or corotating frame (from ancient Greek synodikos,

of the SuneEarth pair. The origin O of related to the assembly) [11] O ¼ O; ! o ;! o ;! o 1

2

3

both frames is the center of mass of the pair (Sun, Earth); the ecliptic and synodic poles are coincident, i.e., i3 ¼ o3, and are orthogonal to the revolution plane of the Earth. The ﬁrst ! ecliptic axis i 1 is aligned with the mean Earth’s vernal equinox, . whereas the ﬁrst synodic axis, ! ! ! which is not inertial, is the Sun-to-Earth direction o 1 ¼ P0 P1 P0 P1 . Since the third body, the Moon, moves around the Earth, it is useful to deﬁne two main orbit periods. Let us denote ! the Earth-to-Moon direction with ! m ¼ ! r ! r ¼ P P . The Moon sidereal r , where ! 1

2

2

2

1 2

period Tm ¼ 2p/um y 27.3 days is the time interval between two successive alignments of

5.5 HILLeCLOHESSYeWILTSHIRE EQUATION

209

! ! m 1 with i 1 , whereas the Moon synodic period Tmo ¼ 2p/umo y 29.5 days is the time interval beo 1 . Clearly the Moon inertial angular rate tween two successive alignments of ! m 1 with ! um ¼ us þ umo is the sum of the Earth revolution rate us (see Eq. (3.92) in Section 3.5.2 and Eq. (4.40) in Section 4.2.5) and of the Moon synodic rate umo, which implies the period relation 1 1 . (5.81) Pmo ¼ 1=Pm ð2pÞ us The association of the Moon synodic period with the period of the Moon phases (fullefull Moon) is immediate.

5.5 HILLeCLOHESSYeWILTSHIRE EQUATION 5.5.1 State Equations and Stability Consider a three-body system {P0,P1,P2} with m0 [ mj, j ¼ 1,2, and the relevant Eq. (5.80). For this of system, as pointed out above, we have ! r 0 ¼ 0, since a generic inertial frame n kind ! ! !o J ¼ O; i 1 ; i 2 ; i 3 has the origin in O ¼ P0. Thus, for the point mass in P0, the second term on the right-hand side of Eq. (5.80) is negligible. In other words, we can assume that each mass mj, j ¼ 1,2, is subjected to the central gravitational ﬁeld generated by the point P0 with mass m0, while the mutual gravity interaction between P1 and P2 is negligible. From Eq. (5.80), the points P1 and P2 obey the restricted two-body equations: ! F1 ! €r ¼ m ! r þ 1 1 m1 r31 ; (5.82) ! m F 2 ! €r ¼ ! r2þ 2 m2 r31 where m ¼ Gm0 and rk ¼ ! r k ; k ¼ 1; 2. The second row of Eq. (5.82) describes the motion of the point mass m2 in the gravitational ﬁeld generated by m0, regardless of the motion of the other point P1. However, in several space applications, it is of interest to describe the relative motion of P2 with respect to the other point P1. A typical example occurs when a deputy spacecraft, represented by P2, has to rendezvous with a chief/reference spacecraft, represented by P1, and both the spacecraft are orbiting around the Earth’s point mass at P0 as in Fig. 5.3. Here, rendezvous (or more generically formation acquisition) means that the deputy satellite approaches the chief spacecraft up to a target distance dt and afterward other operations such as either formation keeping or docking can start. Nonlinear Differential Equations of Relative Motion To properly describe the relative motion between P2 and P1, we deﬁne their relative posir 1 , with r ¼ ! r , and we assume that the orbit of the chief satellite is tion ! r ¼ ! r2! ! ! ! Keplerian, which implies that F 1 ¼ 0 and F ¼ F 2 . Thus, we rewrite Eq. (5.82) as ! ! F m F ! €r ¼ ! €r þ ! €r ¼ m ð! ! ! ! r 1þ rÞþ ¼ ! ð r 1þ r Þþ . (5.83) 2 1 3 ! m2 m r32 2 jr1þ rj

210

5. PERTURBED ORBITAL DYNAMICS

FIGURE 5.3 Reference and perturbed orbits.

This equation describes the motion of P2 relative to P1 in an inertial frame with origin at

p 1; ! p 2; ! p 3 of the O ¼ P0. In the following, we will choose the perifocal frame P ¼ P0 ; ! chief satellite as the inertial frame. n ! ! !o A description of the relative motion in the noninertial trajectory frame T ¼ P1 ; t 1 ; t 2 ; t 3 deﬁned by the orbit of P1 is usually more convenient. To work out such a description, we €r can be written as observe, with the help of Section 6.2.2, that the inertial acceleration ! ! €r ¼

! €r þ 2! u1 ! r_ þ ! u1 ! u1 ! r þ! u_ 1 ! r; t

(5.84)

t

_ € where ! r t and ! r t are the velocity and acceleration vectors in the trajectory frame T , ! respectively, and u 1 and ! u_ 1 are the angular velocity and acceleration of the frame T in the inertial frame J, respectively. If we choose as a trajectory frame the Hill’s frame n ! ! ! o r 1 ; h 2 ; h 3 ¼ ! u 1 of the chief orbit, the subscript t is replaced H ¼ P1 ; h 1 ¼ ! r 1 = ! u 1 = ! by h, and we write the coordinate equations ! ! ! ! ! ! r 1 ¼ r1 h 1 ; ! r ¼ x h 1 þ y h 2 þ z h 3 ¼ Hr ! ! ! ! ! ! r_ h ¼ x_ h 1 þ y_ h 2 þ z_ h 3 ¼ H r_ ¼ H v ! ! ! ! ! €r ¼ x€ h 1 þ y€ h 2 þ z€ h 3 ¼ H r€ h ! ! ! ! ! F ¼ F1 h 1 þ F2 h 2 þ F3 h 3 ¼ H F

; (5.85)

whereh r ¼ [x,y,z] is the i coordinate vector of P2 in the Hill’s frame deﬁned by the vectrix ! ! ! ! H ¼ h 1 h 2 h 3 . The three coordinates are usually known as radial (x), longitudinal

211

5.5 HILLeCLOHESSYeWILTSHIRE EQUATION

or along-track (y), and cross-track or out-of-plane (z). The pair (radial, longitudinal) is known €r holds: as the in-plane pair. With the help Eq. (5.82) and of Eq. (5.83), the relative acceleration ! m! ! €r ¼ r€ ! : 1 1 h1 ¼ 2 h1 r1 ! ! ! m m ! €r ¼ m ! r2 ¼ 3 ! r1þ! r ¼ 3 ðr1 þ xÞ h 1 þ y h 2 þ z h 3 2 3 r2 r2 r2

(5.86)

! In addition, since the chief orbit is Keplerian, by recalling that the orbit pole h 3 is inertial and by denoting with h the magnitude of the chief-orbit angular momentum per unit mass (see Eq. (3.40) in Section 3.3.4) , we can write the differential equations: ! ! h! ! u 1 ¼ q_ 1 h 3 ¼ 2 h 3 ¼ u1 h 3 r1 ! : r_1 ! h ! ! u_ 2 ¼ 2 3 r_1 h 3 ¼ 2u1 h 3 ¼ u_ 1 h 3 r1 rq1

(5.87)

h ¼ r21 q_1 By replacing Eqs. (5.85e5.87) into Eq. (5.84) and after doing some simpliﬁcations, the nonlinear differential equations of the relative motion are found to be: m m F1 x€ ¼ 2u1 y_ þ u_ 1 y þ u21 x þ 3 r1 3 ðr1 þ xÞ þ m2 r1 r2 y€ ¼ 2u1 x_ u_ 1 x þ u21 y

m F2 yþ m2 r32

:

(5.88)

m F3 z€ ¼ 3 z þ m2 r2 They can only be solved by numerical integration. Linear State Equations of the Relative Motion Eq. (5.88) can be linearized for a small separation between the chief and deputy satellite, i.e., for r ¼ ! r2! r 1 << r1 . Using the binomial theorem, the term m r32 can be expanded and approximated in terms of r1 as follows: ! ! 3=2 m m r2 2! r 1 $! r m 3 r 1$ r r2 m x ¼ 3 1þ 2þ y 3 1 þ 2 2 y 3 13 . r32 r1 r1 2 r1 r1 r1 r1 r21 r21

(5.89)

The above approximation only retains a term of the order oðjxj=r1 Þ, where jxj is the radial component of ! r and neglects terms of the order o r2 r21 and higher. By assuming an orbit altitude of h hmin ¼ 300 km,a relative motion up to rmax 0.1 hmin may be of interest, which corresponds to 3 r2 r21 2 30 106 and justiﬁes the approximation.

212

5. PERTURBED ORBITAL DYNAMICS

Replacement of Eq. (5.89) in Eq. (5.88) and elimination of the second-order terms in the coordinates of ! r simpliﬁes the nonlinear equation into the linear time-varying and periodic form: x€ ¼ 2u1 y_ þ u_ 1 y þ u21 þ 2m r31 x þ F1 =m2 (5.90) y€ ¼ 2u1 x_ u_ 1 x þ u21 m r31 y þ F2 =m2 ; 3 z€ ¼ m r1 z þ F3 =m2 where the coefﬁcients u1, u_ 1 , and r1 are periodic with the chief-orbit period Po ¼ 2p/uo, qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ deﬁned by uo ¼ m a31 and by the semimajor axis a1 . By changing the above notations 2 2 into e ¼ e1, q ¼ q1, q_ ¼ u1 , u_ 1 ¼ 2q_ sin qð1 þ e cos qÞ1 and m r31 ¼ q_ ð1 þ e cos qÞ1 , where the last two identities should be proved by the reader, the state equation (5.90) becomes r dr=dt 0 I3 1 0 FðtÞ; (5.91) ðtÞ þ ðtÞ ¼ m2 B2 A21 ðqðtÞÞ A22 ðqðtÞÞ v dv=dt with the following vectors and matrices r ¼ ½x; y; z; v ¼ ½v1 ; v2 ; v3 ; F ¼ ½F1 ; F2 ; F3 ; 2 3 þ e cos q 2e sin q 2 6 q_ 6 e cos q A21 ðtÞ ¼ 6 2e sin q 1 þ e cos q 4 0 0

B 2 ¼ I3 3 2 0 0 7 6 7 6 0 7; A22 ðtÞ ¼ 2q_ 6 1 5 4 1 0

1 0 0

0

3.

7 7 07 5 0

(5.92)

The free response of Eq. (5.91) can be found as in Ref. [9] by changing the independent variable from t to q; in other words, by changing dt with dq in Eq. (5.91) and the matrices in Eq. (5.92) with the following ones: 2 3 3 þ e cos q 2e sin q 0 6 7 1 6 7 A21 ðqÞ ¼ e cos q 0 7 6 2e sin q 5 1 þ e cos q 4 0 0 1 . (5.93) 2 3 e sin q 1 þ e cos q 0 7 6 2 1 6 7 A22 ðqÞ ¼ e sin q 0 7; B2 ¼ 2 I3 6 ð1 þ e cos qÞ 5 1 þ e cos q 4 q_ 0 0 e sin q As a second step, the following transformation is applied to the state vector [r,dr/dq]: r ð1 þ e1 cos qðtÞÞI3 0 rs r ðqÞ ¼ SðqÞ : (5.94) ¼ e1 sin qðtÞI3 vs ð1 þ e1 cos qðtÞÞI3 dr=dq dr=dq The transformation in Eq. (5.94) converts the modiﬁed version of Eq. (5.91), with the new state vector [r,dr/dq] and the matrices of Eq. (5.93), into the TschaunereHempel equations [9].

5.5 HILLeCLOHESSYeWILTSHIRE EQUATION

213

Exercise 10 By taking into account that S(q) in Eq. (5.94) is a function of the independent variable q, prove from Eqs. (5.91), (5.93), and (5.94) the TschaunereHempel equations drs =dq rs 0 0 I3 1 ðqÞ þ ¼ FðqÞ; (5.95) m2 G2 ðqÞ dvs =dq F21 ðqÞ F22 vs with the following matrices 2 1 3ð1 þ e cos qÞ 6 F21 ðqÞ ¼ 4 0 0

0 0 0

3 2 0 2 0 7 6 0 5; F22 ¼ 4 2 0 0 0 1

0

3

7 0 5; G2 ðqÞ ¼ ð1 þ e cos qÞB2 : , 0 (5.96)

HilleClohessyeWiltshire Equation The normalized HCW equation assumes that the chief orbit is circular, that is e1 ¼ 0, and follows from Eq. (5.91) through the following identities and the state transformation S0: u1 ðtÞ ¼ uo ; m r31 ðtÞ ¼ m a31 ¼ u2o ; u_ 1 ðtÞ ¼ 0 ; 2 3 2 3 2 3 I3 0 r r (5.97) 5 4 5 ¼ S0 4 5; S0 ¼ 4 0 u1 v Dr 0 I3 where S0 replaces the relative velocity v with the relative increment Dr ¼ v/uo. Exercise 11 Prove from Eqs. (5.91), (5.92), and (5.97) that the normalized HCW equation can be written in the form: 0 I3 0 1 r_ r ðtÞ ¼ uo ðtÞ þ FðtÞ; (5.98) D_r D_ r m2 uo I3 A21 A22 with the following matrices and vectors 2 3 2 3 0 0 0 6 7 6 A21 ¼ 4 0 0 0 5; A22 ðtÞ ¼ 4 2 0 0

1

0

2 0 0

2 3 3 2 3 vx 0 Dx 1 6 7 7 6 7 0 5; Dr ¼ 4 Dy 5 ¼ 4 vy 5: , uo 0 Dz vz

(5.99)

Since Eq. (5.98) shows that the out-of-plane relative motion is decoupled from the in-plane pair, it is convenient to split Eq. (5.98) into two sets of equations, upon a reordering of the state variables. The in-plane equation holds x_ xy ðtÞ ¼ Axy xxy ðtÞ þ Bxy F12 ðtÞ;

(5.100)

214

5. PERTURBED ORBITAL DYNAMICS

where matrices and vectors are as follows: 3 2 2 0 1 0 0 0 63 0 0 27 6 1 61 7 6 Axy ¼ uo 6 7; Bxy ¼ 6 40 0 0 15 u o m2 4 0 0

2 0

0

0

0

3

2

x

3

6 Dx 7 07 F1 7 6 7 . 7; xxy ¼ 6 7; F12 ¼ 4 y 5 05 F2 Dy 1

The out-of-plane equation is the following undamped oscillator equation: 2 3 " # z z_ ðtÞ ¼ Az 4 5 þ Bz F3 Dz_ Dz 2 3 2 3: 0 1 0 5; Bz ¼ 1 4 5 A z ¼ uo 4 uo m 2 1 0 1

(5.101)

(5.102)

The eigenvalues of the out-plane equation are given by the imaginary pair Lz ¼ L(Az) ¼ {juo}. The eigenvalues of Axy may be obtained through the state transformation 3 2 0 2 0 2 62 0 0 0 7 7 6 (5.103) xxy ¼ Sxz zxy ¼ 6 7zxy ; 44 0 3 3 5 0

4 0

3

which converts Axy into a block diagonal matrix as follows: z_ xy ¼ Fxy xxy þ Gxy F12

Fxy ¼ S1 xy Axy Sxy

2

0

6 61 6 ¼ uo 6 6 60 4 0

1 0 0 0

0 0

3

2

1=2

0

3

7 7 6 7 6 0 0 07 1 7 7 6 7; Gxy ¼ S1 Bxy ¼ 1 6 7. xy 7 6 uo m2 6 2=3 1 7 7 0 17 5 5 4 0 0

0

(5.104)

1

Eq. (5.104) corresponds to the parallel of a double integrator and an oscillator, whose eigenvalues are

Lxy ¼ LðAxy Þ ¼ LðFxy Þ ¼ f0; 0; juo g.

(5.105)

This means that the in-plane relative motion is a combination of unstable and oscillatory components. The relevant motion trajectories become explicit by ﬁnding the free response. The normalized HCW equation in Eq. (5.100) and Eq. (5.102) is the same equation found in Section 3.6 by perturbing a circular orbit. The only difference with respect to Section 3.6 is the

5.5 HILLeCLOHESSYeWILTSHIRE EQUATION

215

! addition of the perturbing acceleration F m2 : The state-variable correspondence with Eq. (3.124), Section 3.6, is the following: x ¼ ½x; Dx; y; Dy; z; Dz5dx ¼ ½dr; dvr ; dy; dvy ; dz; dvz .

(5.106)

By assuming the initial conditions xð0Þ ¼ x0 ¼ ½x0 ; Dx0 ; y0 ; Dy0 ; z0 ; Dz0 ; ! the free response of Eqs. (5.100) and (5.102), under F 2 3 2 3 2 2 x x0 4 3c s 7 6 7 6 6 6 Dx 7 6 6 Dx0 7 3s c 7 6 7 6 6 7 ¼ eAxy t 6 6 7 ¼ 6 6 y 7 6 6ðs uo tÞ 2ð1 cÞ 6 y0 7 5 4 5 4 4 Dy

Dy0

6ð1 cÞ

zðtÞ ¼ zmax sinðuo ðt t0 Þ þ j0 Þ; zmax

¼ 0, is given by 32 3 x0 0 2ð 1 c Þ 76 7 76 Dx0 7 0 2s 76 7 76 7 6 y0 7 ; 1 4s 3uo t 7 54 5

(5.107)

(5.108)

Dy0 2s 0 3 þ 4c qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ z20 þ Dz20 ; j0 ¼ a tan 2ðz0 ; Dz0 Þ

where c ¼ cosðuo tÞ and s ¼ sin(uot). The equation, which is the same as Eq. (3.127) of Section 3.6, shows that the unstable longitudinal motion y has a linear free response in agreement with the double integrator and only affects the longitudinal coordinate y, whereas the radial coordinate x is biased but bounded. T.A Lovell and S. Tragesser [5] employed the periodicity of the free response in Eq. (5.108) to derive six relative orbit elements (ROE) in analogy with the six orbit elements of Keplerian orbits. In fact from Eq. (5.108), in-plane and out-of-plane motions can be rewritten as xðt t0 Þ ¼ xc xmax sinðuo ðt t0 Þ þ qr Þ yðt t0 Þ ¼ yc þ 2xmax cosðuo ðt t0 Þ þ qr Þ ;

(5.109)

zðt t0 Þ ¼ zmax sinðuo ðt t0 Þ þ j0 Þ where the six coefﬁcients, the ROEs, which depend on the initial condition x(t0) ¼ x0 of x in Eq. (5.106), are given by: xc ¼ 2ð2x0 þ Dy0 Þ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 xmax ¼ ðxc x0 Þ þ Dx20 3 yc ¼ y0 2Dx0 xc uo ðt t0 Þ 2 : qr ¼ a tan 2ðDx0 ; xc x0 Þ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ zmax ¼ z20 þ Dz20 j0 ¼ a tan 2ðz0 ; Dz0 Þ

(5.110)

216

5. PERTURBED ORBITAL DYNAMICS

The relative motion geometry becomes clear by observing Eq. (5.109). The out-of-plane motion in the third row is a periodic oscillation of magnitude zmax, which is orthogonal to the chief orbital plane. The in-plane motion in the ﬁrst and second row describes an ellipse with the longitudinal semimajor axis equal to ymax ¼ 2xmax and with the radial semiminor axis equal to xmax. From Eq. (5.110), the center {xc,yc} of the ellipse is drifting in the longitudinal direction (see the third row) at the rate y_c ¼ 3xc uo 2. The drift disappears if xc ¼ 0, which is a well-known fact. In summary, the radial free response is always bounded, and the longitudinal free response becomes bounded by ﬁxing the initial conditions to xc ¼ 2(2x0 þ Dy0) ¼ 0. As already mentioned in Section 3.6, only the longitudinal free response y(t) is unbounded. All the other response terms are bounded, since they are decoupled from y(t), and their local Lyapunov stability around a circular orbit has been proved by Theorem 1 of Section 3.6, which means that the shape and orientation of the perturbed/deputy orbit remains close to the chief/reference circular orbit. Furthermore, if the longitudinal perturbation y is expressed in terms of the true anomaly as in Section 3.6, that is y(t) ¼ adq(t) ¼ a(q2(t) q1(t)), where q2 refers to the deputy orbit and q1 ¼ q10 þ uot refers to the chief orbit, we can write the inequality jyj=a ¼ jdqj p, which is valid if all the other variables remain bounded. This means that the relative position of the two satellites is bounded but drifting along the orbit. This type of local instability is the same instability as a pair of oscillators having the same initial conditions but different natural frequencies. In summary, since |dq| p notwithstanding the local instability of dq, all the state variables remain bounded, which condition is referred to as the stability in the sense of Lagrange, as already mentioned in Section 3.6. Exercise 12 Consider a reference geocentric circular orbit, materialized by a chief spacecraft P1 traveling at an altitude of h ¼ 300 km from the Earth’s equatorial radius. The orbit is referred to as the reference or chief orbit. The chief and deputy orbit parameters are reported in Table 5.1. The deputy parameters are obtained from the chief parameters as p2 ¼ p1(1 þ vp), where p1 refers to the chief satellite, p2 to the deputy, and vp is the fractional difference. Exceptions are the eccentricity e2 and the time t2 of the perigee passage. The initial along-track distance dnom and the time t2 of the perigee passage of the deputy satellite are related by rﬃﬃﬃﬃ dnom m t2 y ; u2 ¼ . (5.111) a 1 u2 a32 The goal is the integration of the nonlinear and HCW Eqs. (5.82) and (5.98), respectively, during 10 orbit periods. SOLUTION

! ! The perturbations F 1 and F 2 in Eq. (5.82) are assumed to include the aerodynamic forces that are signiﬁcant at 300 km altitude. We assume that these forces are accurately cancelled by a drag-free control as explained in Section 11.2. In this case, the residual perturbing acceleration is the combination of the drag residuals of the order of few micronewtons, and of gravitational anomalies such as those due to Earth’s ﬂatness. Fig. 5.4, left, shows the time proﬁle of the relative motion ! r ¼ ! r2! r 1 with the coordinate vector r ¼ [x,y,z] in the chief-

217

5.5 HILLeCLOHESSYeWILTSHIRE EQUATION

TABLE 5.1

Chief and Deputy Orbital Parameters

No

Parameter

Symbol

Unit

Chief (e.g., U1 )

Deputy Fractional Difference (e.g., vU)

1

Right ascension of ascending node

U1,vU

rad

p/2

0.0001

2

Inclination

i1,vi

rad

p/2 þ 0.113

0.0001

3

Argument of perigee

u1,vu

rad

0

0

4

Semimajor axis

a1,va

km

6678.14

0.00005

5

Satellite mass

mk,k ¼ 1,2

kg

1000

0

Chief

Deputy

0

0.0001

6

Eccentricity

e1,e2

7

Time of the perigee passage

s

0

0.518

8

Initial distance (along track)

t1 ¼ t0,t2 d ¼ ! r ðt0 Þ

km

NA

4.36 (dnom ¼ 4 km)

9

Orbital angular rate

uo

mrad/s

1,16

NA

10

Orbital period

To ¼ 2p/uo

s

5430

NA

11

Simulation time unit

Ts

s

1

NA

12

Maximum command force

Fmax

N

2

NA

NA, not applicable.

FIGURE 5.4 Hill-frame coordinates of the relative motion (deputy minus chief). Left: time proﬁle. Right: Threedimensional curve.

orbit Hill’s frame H ¼

n ! ! ! o! r 1 ; h 2 ; h 3 ¼ ! u 1 r . The numerical solution P1 ; h 1 ¼ ! r 1 = ! u 1 = !

of the nonlinear state equations in Eq. (5.82) and that of the HCW equation in Eq. (5.98) are shown in Fig. 5.4, right. Both solutions overlap without visible differences. In reality, they are progressively diverging as soon as the distance between the satellites increases.

218

5. PERTURBED ORBITAL DYNAMICS

In terms of the orbital elements, the drift of the longitudinal coordinate in Fig. 5.4 is due to the semimajor axis fractional difference va. The radial and longitudinal oscillations are mainly due to the eccentricity difference e2 e1. The oscillations of the cross-track coordinate are due to the fractional differences vU and vi of the right ascension of the ascending node and of the inclination, respectively. We remark that the simulated chief orbit is exactly circular, and therefore can be compared with the HCW solution. An elliptical orbit should be compared with the linear time-varying Eq. (5.90). A chief circular orbit may be also interpreted as a virtual reference orbit. ,

5.5.2 Feedback Stabilization Section 5.5.1 has shown that the relative motion of two bodies traveling along two slightly different orbits is locally unstable, in the sense that the positions of the two bodies tend to diverge from each other as time increases. We are now going to see that the two bodies can be kept close to each other by means of an appropriate feedback law, which, in the case of a chief circular orbit, can be designed on the basis of the HCW Eq. (5.98). To this end, a transformation of Eq. (5.98) that was suggested in [1] replaces Dy in Eq. (5.106) with xc/2, which is a factor of the longitudinal drift in Eq. (5.110). By deﬁning the new longitudinal increment dy ¼ xc =2 ¼ 2x þ Dy;

(5.112)

and by updating x in Eq. (5.106) as follows x ¼ ½xx ; xy ; xz ¼ ½x; Dx; y; dy; z; Dz;

(5.113)

the HCW Eq. (5.98) becomes x_ ðtÞ ¼ AxðtÞ þ Bu ðuðtÞ þ dðtÞÞ; uðtÞ ¼

Fu ðtÞ Fd ðtÞ ; d ðtÞ ¼ ; m2 u2o m2 u2o

where the input force has been split into the command Fu ¼ ½Fux ; Fuy ; Fuz bance Fd ¼ [Fdx, Fdy, Fdz], the vectors u and d are in length units, and the follows: 2 3 2 3 Axx Axy 0 Bxx 0 0 6 7 6 7 7 6 7 A ¼ uo 6 4 Ayx Ayy 0 5; Bu ¼ uo 4 0 Byy 0 5 0 0 Azz 0 0 Bzz " # " # " # " 0 1 0 1 0 0 2 Axx ¼ Azz ¼ ; Ayy ¼ ; Axy ¼ ; Ayx ¼ 1 0 0 0 0 2 0 " # 0 B ¼ Bxx ¼ Byy ¼ Bzz ¼ 1

(5.114) and the disturmatrices are as

0 0

# ;

(5.115)

219

5.5 HILLeCLOHESSYeWILTSHIRE EQUATION

The main property of Eq. (5.114) is that the state matrix A is neither block diagonal nor block triangular, but the eigenvalues are those of the three diagonal blocks. The natural question is whether the same property could be maintained under state feedback. The answer is positive under conditions to be found in the next subsection on the decoupled state feedback. Decoupled State Feedback Let us denote the feedback matrix K of the state feedback command u(t) ¼ Kx(t) as follows: 2 3 2 3 Pxx Kxy 0 px1 px2 kx3 kx4 0 0 6 6 7 7 0 5; K ¼ 4 Kyx Pyy 0 5 ¼ 4 ky1 ky2 py3 py4 0 (5.116) 0

0

Pyy

0

0

0

0

pz1

pz2

where the set of feedback gains {Pxx,Pyy,Pzz} is the design set which must be given an expression in terms of the closed-loop eigenvalues: qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ lx1;x2 ¼ zx 1 z2x uo ; zx > 0

ly1;y2 ¼ uo gy1 ; uo gy2 ; gy1 > 0; gy2 > 0 . (5.117) qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ lz1;z2 ¼ zz 1 z2z uo ; zz > 0 The second pair of eigenvalues has the goal of bounding the longitudinal relative motion, that is of zeroing the drift scale factor xc in Eq. (5.110), whereas the complex eigenvalues in the ﬁrst and third row have the goal of bringing to zero or to a desired value the in-plane and out-of-plane radiuses xmax and zmax in Eq. (5.110). To this end, the closed-loop natural frequencies are kept equal to the chief orbital rate uo. Exercise 13 Find the expression of the scalar gains in {Pxx,Pyy,Pzz} in terms of the eigenvalues in Eq. (5.117). , The cross-coupling gains {Kxy,Kyx} must guarantee that the closed-loop eigenvalues coincide with those of the diagonal blocks of A BuK. The following Lemma has been proved in Ref. [1]: Lemma The eigenvalues of A BuK are those of the diagonal blocks {Ak BuPk, k ¼ xx,yy,zz} if and only if 1 0 Ayx BKyx ¼ 0; B ¼ Axy BKxy lI Ayy þ BPyy : (5.118) 1 The ﬁrst identity of Eq. (5.118) corresponds to the following nonlinear equations in terms of the scalar entries of {Kxy,Kyx}: 2kx3 ky1 ð2 kx4 Þ ¼ 0 kx3 ðky1 þ 2py2 Þ þ 2py1 ð2 kx4 Þ ¼ 0 ., ky2 kx3 ¼ 0 ky2 ð2 kx4 Þ ¼ 0

(5.119)

220

5. PERTURBED ORBITAL DYNAMICS

Eq. (5.119) possesses two sets of solutions, but only one of them, referred to as the position decoupling solution in [1] provides cross-coupling gains which depend on the design gains and therefore can be optimized. Exercise 14 Prove that the two sets which solve Eq. (5.119) are fkx4 ¼ 2; kx3 ¼ 0; ky1 ; ky2 anyg rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ o (5.120) n . ky2 ¼ 0; ky1 ¼ kc ¼ py2 1 1 4py1 p2y2 ; kx3 ¼ ky1 ð1 kx4 =2Þ; kx4 any . Then prove that only the second set with kx4 ¼ 0 and kx3 ¼ ky1 ¼ kc provides nonzero gains with degrees of freedom that allow optimization. Finally prove that

minfgy1 ;gy2 gjkc j ¼ 2min gy1 ; gy2 : , (5.121) It is straightforward to ﬁnd the values of the four eigenvalue parameters {zx,gy1,gy2,zz} in Eq. (5.117), which trade-off between the fastest response and the least command authority. By augmenting the state feedback with the disturbance rejection and by ﬁxing {xr,yr,0} as the reference values of the relative position coordinates {x,y,z} of x in Eq. (5.113), the following control law results:

Fux ¼ m2 u2o 2zx Dx 2min gy1 ; gy2 ðy yr Þ þ 3xr Fdx

Fuy ¼ m2 u2o gy1 gy2 ðy yr Þ þ gy1 þ gy2 ðdy 2xr Þ 2min gy1 ; gy2 ðx xr Þ Fdy . (5.122) Fuz ¼ m2 u2o 2zz Dz Fdz Let us remark that, for xr s 0, the radial command component Fux remains different from zero. This corresponds to the fact that xrs 0 is not a natural equilibrium, but a forced equilibrium of Eq. (5.114), unlike yr which is a natural equilibrium whichever be its value. Exercise 15 Prove that the control law in Eq. (5.122) guarantees that the tracking error e xr ¼ x xr of the LTI Eq. (5.114) asymptotically converges to zero for zx > 0; gy1 > 0; gy2 > 0; zz > 0 . , Exercise 16 Prove by simulation that the control law in Eq. (5.122), with the gains in Table 5.2, stabilizes the relative motion in Fig. 5.4 and forces the deputy spacecraft to approach the chief orbit with the target relative position in Table 5.2. Assume Fdj ¼ 0, j ¼ x,z,y, in Eq. (5.122). Prove that the maximum command force in Hill-frame coordinates satisﬁes jFdj j 2 N. The longitudinal gains in Table 5.2 have been chosen by a simple trade-off between the maximum command force and the chief-orbit acquisition time. The control law in Eq. (5.122) must be accompanied by a state predictor capable of predicting Fdj and by a reference generator which ﬁxes the deputy satellite trajectory to be tracked. Since both simulated chief

221

5.5 HILLeCLOHESSYeWILTSHIRE EQUATION

TABLE 5.2

Target Position and Control Gains of the Deputy Spacecraft

No

Parameter

Symbol

Unit

Value

1

Target relative position (chief Hill’s frame)

{xr,yr,zr}

m

{20,20,0}

2

Longitudinal gains

{gy1,gy2}

{0.2,0.2}

3

Radial and cross-track damping ratio

{zx,zz}

{0.5,0.5}

and deputy orbits are close to be geodesic (only subject to gravitational forces) because of a drag-free control (see Section 11.2), the only perturbing forces are of gravitational nature and correspond to gravity anomalies dominated by the Earth’s ﬂatness. For low-Earth orbits, the state vector x in Eq. (5.113) can be directly reckoned from the differential range and rate of chief and deputy Global Navigation Satellite System receivers (see Chapter 8) and by the radio transmission of the chief spacecraft data to the deputy. Spacecraft rendezvous guided by optical sensors is not explicitly treated. SOLUTION

Fig. 5.5 shows the time proﬁle and three-dimensional (3D) curve of the relative motion ! r ¼ ! r2! r 1 , which converges to the reference relative position in about six orbit periods. The converging 3D trajectory is the solid line in the north-west corner of Fig. 5.5, right. The long spiral-like trajectory (dashed line) is the open-loop trajectory of Fig. 5.4, right, for comparison. Fig. 5.6, left, which is the enlargement of Fig. 5.5, left, shows the reference acquisition of the three relative coordinates. Fig. 5.6, right, shows the Hill-frame coordinates of the command force. The dominating coordinate is the radial coordinate Fux, since the second term of Fux in Eq. (5.122) is proportional to the longitudinal tracking error e yr ¼ y yr . The maximum force value is 2 N as expected. At the onset, Fux undergoes a short saturation interval which

FIGURE 5.5

Closed-loop relative motion in Hill-frame coordinates. Left: time proﬁle. Right: 3D curve of the openloop (dashed line) and closed-loop (solid line) motion.

222

5. PERTURBED ORBITAL DYNAMICS

FIGURE 5.6 Left: reference acquisition from Fig. 5.5, left. Right: Force command in Hill-frame coordinates.

cannot be appreciated in Fig. 5.6. Let us observe that as previously said the radial command component converges to a nonzero value because of the nonzero radial reference. During a target-approaching maneuver like that shown in Fig. 5.6, left, a collision avoidance strategy should be in operation. It must be designed on the basis of measurement and control errors. Fig. 5.7 shows that the collision avoidance threshold of 20 m is neither touched nor crossed. The threshold value of 20 m should be compared with the braking stroke during approach. By assuming an approaching speed less than vmax 20 mm/s, which is achieved at 30 ks, the data of Table 5.1 guarantee that the braking time is bounded by sbrake m2vmax/ Fmax ¼ 10 s, which corresponds to a braking stroke well less than 0.2 m. ,

FIGURE 5.7 Distance [m], relative speed magnitude [mm/s] and collision avoidance threshold during the target approaching phase.

5.6 RESTRICTED THREE-BODY PROBLEM

223

5.6 RESTRICTED THREE-BODY PROBLEM 5.6.1 State Equations Consider the three-body system {P0,P1,P2} in Fig. 5.8 with m1 mi, i ¼ 0,2, and assume that the larger-mass points P0 and P2 (the primary bodies) travel along circular orbits about their common CoM denoted by O. The angular speed of rotation is pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ uo ¼ ðm0 þ m2 Þ=d3 ; (5.123) ! where mi ¼ Gmi and d ¼ P0 P2 is the distance between the two bodies.

! ! ! Let O ¼ O; o 1 ; o 2 ; o 3 be the common orbital frame of reference, chosen to coincide with the synodic frame already introduced in Section 5.4. We recall that the frame is deﬁned by the orbital plane of the primary bodies, that ! o 1 and ! o 2 lie on the orbital plane, that . ! ! ! ! o 1 ¼ P0 O P0 O is the ﬁrst synodic axis, directed from P0 to O on the segment P0 P2 joining the primary bodies, that ! o 3 is the inertial orbital pole and ! o2 ¼ ! o3! o 1 is the second syn! odic axis. Because of the rotating segment P0 P2 , the frame O is not inertial, and the pair

! ! o 1 ; o 2 rotates in the orbital plane around the pole with the angular rate uo of Eq. (5.123). The frame angular rate vector can be denoted by ! u o ¼ uo ! o 3 , and the derivatives of the frame axes are the following: ! o_ 1 ¼ ! uo ! o 1 ¼ uo ! o 2; ! o_ 2 ¼ ! uo ! o 2 ¼ uo ! o 1; ! o_ 3 ¼ 0:

(5.124)

This section is a synthesis of the treatments in Refs. [2,4,11,13]. The coordinates in the common orbital frame O of the position ! r 1 , velocity ! r_ 1 , and accel! € eration r 1 of the small point mass P1 must account for the frame rotation and hold:

Synodic frame

o2

P1

r10 r1 P0 Earth

r12

P2

O d 0o1

Spacecraft

CoM

– d 2o1

Moon

d

FIGURE 5.8

Three-body geometry: Earth, Moon, and spacecraft.

224

5. PERTURBED ORBITAL DYNAMICS

! r 1 ¼ x! o 1 þ y! o 2 þ z! o3 ! _r ¼ x_! ! ! o 1 þ y_ o 2 þ z_ o 2 þ x! o_ 1 þ y! o_ 2 þ z! o_ 3 ¼ ðx_ yuo Þ! o 1 þ ðy_ þ xuo Þ! o 2 þ z_! o 3; 1 ! ! ! €r ¼ x€ 2u y_ u2 x ! 2 o 1 þ y€ þ 2uo x_ uo y o 2 þ z€ o 3 1 o o (5.125) o 1 and 2uo x_! o 2 are the Coriolis accelerations, and u2o x! o 1 and u2o y! o 2 are the where 2uo y_! centripetal accelerations. With the help of Fig. 5.8, the relative positions of P1 with respect to P0 and P2, and their magnitudes, are found to be: ! r 1! r 0 ¼ d0 ! o1þ! r1 ! ! ! r 1 r 2 ¼ d2 ! o1þ r1 qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 ! ! r10 ¼ j r 1 r 0 j ¼ ðx þ d0 Þ þ y2 þ z2 ; qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 r ! r j ¼ ðx d Þ þ y2 þ z2 r ¼ j! 12

1

2

(5.126)

2

! ! where d0 ¼ P0 O , d2 ¼ P2 O , and the distance d already used in Eq. (5.123) holds d ¼ d 0 þ d 2. Replacement into Eq. (5.80) of the relative positions in Eq. (5.126) and of the acceleration coordinates from Eq. (5.125) provides the well-known circular restricted three-body equations: x€ 2uo y_ u2o x ¼ m0 y€ þ 2uo x_ u2o y ¼ m0

x þ d0 x d2 F11 m2 3 þ m1 r310 r12

y y F12 m2 3 þ ; 3 r10 r12 m1

z€ ¼ m0

z z F13 m2 3 þ 3 r10 r12 m1

(5.127)

! where the components of F 1 in the common orbital frame O are collected in F1 ¼ [F11,F12,F13]. _ y_ and z. _ Eq. (5.127) must be The state variables in Eq. (5.127) are given by x, y and z and by x; completed with the orbit equations of P0 and P2, which are given by Eq. (5.80), with ! mi [ m1, i ¼ 0,2 and F i ¼ 0; i ¼ 0; 2:

5.6.2 Free Response: The Unique Known Constant of Motion By setting F1 ¼ 0 in Eq. (5.127), one can derive the free response of the nonlinear restricted three-body equations. In principle, the method might be the same as for the two-body problem in Chapter 3. There, ﬁve nonlinear combinations of the state vector [r,v]which remain constant in time were found, namely the angular momentum per unit mass h (three constants of motion) and the eccentricity vector e (only two constants of motion because of the orthogonality to h). Then, the time response of the remaining degree of freedom, the eccentric anomaly E, was computed using Kepler’s equation. Unfortunately, only a single constant of motion is known for the present problem.

5.6 RESTRICTED THREE-BODY PROBLEM

225

Let us deﬁne the scalar function Ug ðx; y; zÞ ¼

m0 m2 þ ; r310 r312

(5.128)

where the subscript g stands for gravity, and whose gradient holds vUg x þ d0 x d2 ¼ m0 3 m2 3 vx r10 r12 vUg y y ¼ m0 3 m2 3 vy r10 r12

.

(5.129)

vUg z z ¼ m0 3 m2 3 vz r10 r12 Replacement of Eq. (5.129) in Eq. (5.127) and the identity F1 ¼ 0 yield the following autonomous equations: v x€ 2up y_ u2p x ¼ Ug ðx; y; zÞ vx y€ þ 2up x_ u2p y ¼ z€ ¼

v Ug ðx; y; zÞ: vy

v Ug ðx; y; zÞ vz

(5.130)

_ 2y, _ and 2z, _ respectively, and addition of the Multiplication of the above equations by 2x, resulting equations yield the following identities: _z u2p xx_ u2p yy_ ¼ x_ x_x€ þ y_y€ þ z€

v v v d Ug þ y_ Ug þ z_ Ug ¼ Ug . vx vy vz dt

Integration by parts of Eq. (5.131) results into the constant of motion: x_2 þ y_2 þ z_2 u20 x2 þ y2 2Ug ¼ constant ¼ CJ ;

(5.131)

(5.132)

where CJ is a constant, known as the Jacobi integral or Jacobi constant. Eq. (5.132) can be rewritten as 1 1 2 x þ y2 þ z 2 ; U ðx; y; zÞ CJ ¼ 2 2

(5.133)

where the functional Uðx; y; zÞ ¼ Ug ðx; y; zÞ þ

u2o 2 x þ y2 2

(5.134)

is called the effective potential. The name indicates that U(x,y,z), and not Ug, is the effective gravitational potential at the point P1, since the coordinates [x,y,z] of P1 are given in the rotating orbital frame O. The rotating frame is accounted for by the additional term u2o x2 þ y2 2.

226

5. PERTURBED ORBITAL DYNAMICS

5.6.3 Free Response: Lagrangian Equilibrium Points and Stability As a second step in the free-response study, we derive the equilibrium points of the circular ! restricted three-body equation, corresponding to F 1 ¼ 0. The points can be found by setting all the derivatives to zero in Eq. (5.127), which leads to the algebraic equations: x þ d0 x d2 þ m2 3 u2o x ¼ 0 3 r10 r12 y y m0 3 þ m2 3 u2o y ¼ 0 . r10 r12 z z m0 3 þ m2 3 ¼ 0 r10 r12

m0

(5.135)

h i The solutions of these equations correspond to the equilibrium points xj ; y ; zj ; j ¼ 1; .; n, j where n is the solution cardinality. As a ﬁrst remark, the third equation is solved by zj ¼ 0 whichever be j, which fact implies that all the equilibrium points lie in the {x,y} plane. The second equation in Eq. (5.135) is solved either by y ¼ 0 or by j

g1 ¼

m0 m2 þ ¼ u2o . r310 r312

(5.136)

Let us suppose that y s0. Replacement of u2o from Eq. (5.136) in the ﬁrst equation of Eq. j (5.135) yields m0

d0 d2 m2 3 ¼ 0: 3 r10 r12

(5.137)

Since the origin O of the common orbital frame O is the CoM of the larger-mass points P0 and P2, the identity m0d0 ¼ m2d2 holds, implying that r10 ¼ r12 . Moreover, from Eqs. (5.126) and ! ! ! ! (5.136), it follows that r10 ¼ r12 ¼ d ¼ d0 þ d2 , i.e., that r 1 r 0 ¼ r 1 r 2 ¼ d. The last equation identiﬁes two equilibrium points, denoted by L4 and L 5, which are the in ! tersections of the two circles deﬁned by r 1 ! r 0 ¼ d and ! r 1! r 2 ¼ d and shown in Fig. 5.9. Let us consider now the case y ¼ 0. The ﬁrst equation in Eq. (5.135) can be rewritten as j

xj þ d0 xj d2 2 f xj ¼ m 0 3 þ m2 3 uo x j ¼ 0 xj þ d0 xj d2 ; m0 m2 2 2 0 3 þ 3 ¼ uo þ Duo xj þ d0 xj d2

(5.138)

227

5.6 RESTRICTED THREE-BODY PROBLEM

FIGURE 5.9 The Lagrangian equilibrium points L1 to L5 of the EartheMoon system.

with 0 Du2o ¼

1

C 1B B m2 d2 m0 d0 C. 3 3A @ xj xj d2 xj þ d0

The term Du2o can be shown to be positive for all the three solutions of Eq. (5.138) (the collinear points deﬁned in the next paragraph) as reported in Table 5.3. No analytic solution can be found of this equation, but it can be numerically and graphically solved. Fig. 5.10 shows that Eq. (5.138) has three solutions, which identify the other three equilibrium points, denoted by L1, L2, and L3.

TABLE 5.3

EartheMoon Lagrangian Points

No.

Point

Earth to Point Distance [km]

Moon to Point Distance [km]

Du2o u2o

1

CoM

4700 (þ)

384,400

NA

2

L4, L5

384,400

384,400

0

3

L1

326,400

58,000

4.1

4

L2

448,900

64,500

2.2

5

L3

381,700

766,100

0.01

(þ), inside the Earth body; NA, not applicable.

228

FIGURE 5.10

5. PERTURBED ORBITAL DYNAMICS

The equilibrium points L1, L2, and L3 together with the Earth and Moon represented by point

masses.

In summary, the free circular restricted three-body Eq. n (5.127) o has the ﬁve equilibrium points shown in Fig. 5.9. They are denoted by Lj ¼ xj ; y ; 0 and are known as the j

Lagrangian points or libration points. They can be subdivided into two groups: 1. Three collinear points denoted by L1, L2, and L3. These are unstable equilibrium points (see Section 5.6.4 for a proof), but they can be stabilized by a small control action. Space applications include the SuneEarth L1 point for the observation of the Sun and solar wind, the SuneEarth L2 point for astronomical observation, the EartheMoon L2 point for communications with the “far side” of the Moon. 2. Two Trojan points L4 and L5 . The Trojan attribute comes from the asteroids that share the orbit of Jupiter around the Sun and librates around the SuneJupiter L4 (the Greek camp) and L5 (the Trojan camp) Lagrangian points. By convention, they are named from mythological characters of the Trojan war, hence the attribute Trojan. The Trojan points are marginally stable for m0/m2 > 24.96 (see Section 5.6.4 for a proof). The EartheMoon Trojan points are marginally stable if one neglects the gravitational attraction of the Sun. Possible perturbations of the free circular restricted three-body system stem from the Sun’s gravity ﬁeld and radiation pressure, which may generate large orbital motions around the Trojan points, with nonzero orbital eccentricity. Exercise 17 THE EARTHeMOON LAGRANGIAN POINTS

With reference to Fig. 5.8 and to Table 3.1 of Section 3.2, the EartheMoon system is characterized by the following data: m0 ¼ 0.3986 1015 m3/s2, m2 ¼ 0.0123 m0 ¼ 4.9 1012 m3/s2, d ¼ d0 þ d2, d0 ¼ 4760 km, d2 ¼ 380,000 km, and uo ¼ 2.662 mrad/s. The EartheMoon Lagrangian points are located as shown in Table 5.3. Carry out a numerical simulation to study the behavior of a spacecraft (the point mass P1) in a neighborhood of the EartheMoon Lagrangian point L4. The following initial conditions _ _ are assumed: x(0) ¼ 188,000 km, y(0) ¼ 337,000 km, z(0) ¼ 0 m, xð0Þ ¼ 30 m=s, yð0Þ ¼ 0 m=s,

5.6 RESTRICTED THREE-BODY PROBLEM

FIGURE 5.11

229

Left: simulated spacecraft trajectory around the EartheMoon Lagrangian point L4. Right:

enlargement.

! _ and zð0Þ ¼ 0 m=s. Numerically integrate Eq. (5.127) during 3 107 s under F 1 ¼ 0. The spacecraft 2D trajectory coming out of the simulation is the faint complex curve in Fig. 5.11, left and right. The curve is shown together with the two circles deﬁning the Lagrangian points L4 and L5. We observe that, starting from an initial position close to L4, the motion remains close to L4, thus conﬁrming the marginal stability of this equilibrium point. ,

5.6.4 Linearized Equations of Motion and Stability Analysis In this section, we linearize the circular restricted three-body Eq. (5.127). This operation is important for analyzing the stability properties of the Lagrangian points and for studying important types of orbits known as quasiperiodic (or Lissajous) orbits and Halo orbits. Let ! r ¼ x! o 1 þ y! o 2 þ z! o3

(5.139)

be the position vector of a generic Lagrangian equilibrium point with z ¼ 0 and let d! r ¼ dx! o 1 þ dy! o 2 þ dz! o ¼ ! r1! r

(5.140)

be the position vector of P1 relative to this equilibrium point, with dz ¼ z. With reference to the notations of Eq. (5.127), we rewrite r10 and r12 as 3=2 2 3 2 ! ! ! ! ! 3 2 r10 ¼ j r 1 r 0 j ¼ j r þ d r r 0 j ¼ ðx þ dx þ d0 Þ þ y þ dy þ z . (5.141) 3=2 2 3 2 ! ! ! ! ! 3 2 r12 ¼ j r 1 r 2 j ¼ j r þ d r r 2 j ¼ ðx þ dx d2 Þ þ y þ dy þ z The binomial theorem allows us to expand these equations and to keep the ﬁrst two terms of the series: 2 2 3 2 r3 þ o jdxj ; jdyj ; jdxdyj 10 ¼ r10 1 3r10 ðx þd0 Þdx þ y dy : 2 2 3 3 2 r12 ¼ r12 1 3r12 ðx d2 Þdx þ y dy þ o jdxj ; jdyj ; jdxdyj (5.142)

230

5. PERTURBED ORBITAL DYNAMICS

From Eqs. (5.127), (5.135), and (5.142), we obtain the linearized three-body differential equations: 2 2 3 5 5 dx€ 2uo dy_ u2o m0 r3 m r þ 3m r ðx þ d Þ þ 3m r ðx d Þ dx 0 2 2 12 0 10 2 12 10 F11 5 3y m0 r5 10 ðx þ d0 Þ þ m2 r12 ðx d2 Þ dy ¼ m1 5 dy€ þ 2uo dx_ 3y m0 r10 ðx þ d0 Þ þ m2 r5 12 ðx d2 Þ dx

.

(5.143)

5 F12 3 2 5 dy ¼ u2o m0 r3 10 m2 r12 þ 3y m0 r10 þ m2 r12 m1 F13 3 dz€ þ m0 r3 10 þ m2 r12 dz ¼ m1 By deﬁning the following four parameters 3 g1 ¼ m0 r3 10 þ m2 r12 5 g2 ¼ 3y2 m0 r5 10 þ m2 r12 ; 5 g3 ¼ 3y m0 r5 10 ðx þ d0 Þ þ m2 r12 ðx d2 Þ 5 g4 ¼ 3m0 r5 10 ðx þ d0 Þ þ 3m2 r12 ðx d2 Þ 2

(5.144)

2

Eq. (5.143) is transformed into the LTI differential equations: F11 dx€ ¼ 2uo dy_ þ u2o g1 þ g4 dx þ g3 dy þ m1 F12 dy€ ¼ 2uo dx_ þ g3 dx þ u2o g1 þ g2 dy þ . m1 d€ z ¼ g1 dz þ

F13 m1

(5.145)

The third equation is decoupled from the other two equations, and is marginally stable, being the equation of an undamped oscillator. The ﬁrst two equations can be studied as follows. We ﬁrst consider the case where the equilibrium point is a collinear point, either L1, L2, or L3, as deﬁned in Section 5.6.3 and shown in Fig. 5.9. In this case, we have that y ¼ 0, which implies from Eq. (5.144) that g2 ¼ g3 ¼ 0, r10 ¼ jx þ d0 j, and r12 ¼ jx d2 j. The last two identities in turn provide the new equality g4 ¼ 3g1. Hence, the ﬁrst two equations in Eq. (5.145) can be rewritten as follows: F11 dx€ ¼ 2uo dy_ þ u2o þ 2g1 dx þ m1 . 2 F12 dy€ ¼ 2uo dx_ þ uo g1 dy þ (5.146) m1

231

5.6 RESTRICTED THREE-BODY PROBLEM

The characteristic polynomial P123(s) of Eq. (5.146) is found to be: P123 ðsÞ ¼ s4 þ 2u20 g1 s2 þ u20 g1 u20 þ 2g1 ¼ s4 þ a2 s2 þ a0 .

(5.147)

The roots {li,i ¼ 1,.,4} of the polynomial are pﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃ 1 w1 > 0; l2 ¼ w1 < 0; w1 ¼ g 2u20 þ d > 0 2 1 ; (5.148) pﬃﬃﬃﬃﬃﬃﬃﬃ 1 2 l3;4 ¼ j jw2 j ¼ juxy ; w2 ¼ 2u0 g1 þ d < 0 2 qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ where d ¼ g1 9g1 8u20 . It has been shown in Eq. (5.138) that g1 ¼ u2o þ Du2o holds, and l1 ¼

in Table 5.3 that the inequality Du20 > 0 holds for any collinear point. Thus, the constant coefﬁcient a0 in Eq. (5.147) is negative, whereas a2 may be positive or negative depending on the collinear point. This implies that a22 4a0 > 0 and therefore that w1 > 0 and l1 > 0 as shown in Eq. (5.148). Thus, the Eq. (5.146) is proved to be unstable. The eigenvalues are related to the coefﬁcients of P123(s) by a0 ¼ u20 g1 u20 þ 2g1 ¼ w1 jw2 j ¼ w1 u2xy . (5.149) a2 ¼ 2u20 g1 ¼ w1 þ jw2 j ¼ w1 þ u2xy We now consider the case when the equilibrium point is a Trojan point, either L4 or L5, as deﬁned in Section 5.6.3. In this case, Eq. (5.136) and the ﬁrst row of Eq. (5.144) imply the equality g1 ¼ u2o . The ﬁrst two equations in Eq. (5.145) can thus be rewritten as follows: F11 m1 F12 dy€ ¼ 2uo dx_ þ g3 dx þ g2 dy þ m1 dx€ ¼ 2uo dy_ þ g4 dx þ g3 dy þ

:

The characteristic polynomial P45(s) of Eq. (5.150) is found to be P45 ðsÞ ¼ s4 þ 4u2o g2 g4 s2 þ g2 g4 g23 ¼ s4 þ a2 s2 þ a0 .

(5.150)

(5.151)

We have already proved the identities r10 ¼ r12 ¼ d and m0d0 ¼ m2d2 for the Trojan points. Since d is the Earth-to-Moon distance, the three points P0, P2, and either of Li, i ¼ 4,5, constitute an equilateral triangle, which in turn implies that x þ d0 ¼ x d2 ¼ d=2, . pﬃﬃﬃ x=d ¼ 1=2 d0 =d, and that y d ¼ 3 2. The positive and negative sign of the last identity corresponds to L4 and L5, respectively. Since the expression of the four coefﬁcients in Eq. (5.144) simpliﬁes to: pﬃﬃﬃ 9 2 3 3 2d0 2 3 2 g2 ¼ uo ; g2 ¼ uo ; g3 ¼ (5.152) 1 uo ; g4 ¼ u2o ; 4 4 4 d

232

5. PERTURBED ORBITAL DYNAMICS

the coefﬁcients of P45(s) in Eq. (5.151) hold: 9 3 2 a2 ¼ 4 uo ¼ u2o > 0 4 4 : 2 ! 27 d0 27 d0 d0 27 4 4 4 1 12 bð1 bÞuo > 0 1 a0 ¼ uo ¼ u ¼ 16 4 d0 þ d2 4 d0 þ d2 d0 þ d2 o (5.153) The coefﬁcient b in the second row follows from the identities m0d0 ¼ m2d2 and mj ¼ Gmj, j ¼ 0,2, and holds b ¼ d0 =d ¼ m2 =ðm0 þ m2 Þ ¼ m2 =ðm0 þ m2 Þ < 1. Both coefﬁcients a0 and a2 in Eq. (5.153) are positive, but this fact is not sufﬁcient to prove stability, since the polynomial is fourth order, and we also need the inequality 1 4a0 a22 ¼ 27bð1 bÞ > 0. The roots li, i ¼ 1,.,4 of P45(s) are as follows: pﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃ l1;2 ¼ w1 ; l3;4 ¼ w2 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ u2 w1 ¼ o 1 þ 1 27bð1 bÞ (5.154) 2 2 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ u w2 ¼ o 1 1 27bð1 bÞ ; 2 and two cases can be distinguished: 1. Marginal stability. Under 1 27b(1b) > 0, the pair w1,2 is negative, and the four roots pﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃ li, i ¼ 1,.,4 are imaginary: l1;2 ¼ j jw1 j; l3;4 ¼ j jw2 j. Thus, Eq. (5.146) is marginally stable. 2. Instability. In the opposite case of 27b(1b) > 1, the pair w1,2 is complex conjugate with negative real part, namely w1;2 ¼ wexpð jðp=2 þ qÞÞ with 0 < q < p/2 and w ¼ w1;2 . pﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃ The square root holds w1;2 ¼ wexpð jðp=4 þ q=2ÞÞ and includes a complex pair with positive real part, which implies instability. The marginal stability condition of the Trojan points is b2 b þ 1 27 0, with 0 < b < 1, which in turn splits into the following two-side inequalities: rﬃﬃﬃﬃﬃ ! 1 23 0> 1, that is to b << 1. The left-hand-side inequality of the ﬁrst row in Eq. (5.155), rewritten in terms of the mass ratio m0/m2, becomes: rﬃﬃﬃﬃﬃ !, rﬃﬃﬃﬃﬃ ! m0 23 23 > 1þ 1 ¼ 24:96; (5.156) 27 27 m2 and is consistent with the Earth/Moon mass ratio mEarth =mMoon y81:3.

233

5.6 RESTRICTED THREE-BODY PROBLEM

5.6.5 Lissajous and Halo Orbits The motion in a neighbourhood of a collinear point is described by the two equations in Eq. (5.146), together with the third equation in Eq. (5.145). The three equations can be rewritten in the state-space form as d_rxy ¼ Axy drxy þ d€ z ¼ g1 dz þ

u m1

F13 m1

;

(5.157)

where g1 ¼ u2o þ Du2o with Du2o > 0 and the following vectors and matrices hold: 3 3 2 3 2 2 0 0 1 0 0 dx 6 7 6 0 7 6 0 0 0 1 7 7 7 6 dy 7 6 6 ; Axy ¼ 6 2 drxy ¼ 6 7; u ¼ 6 7 7. 4 dx_ 5 4 F11 5 4 uo þ 2g1 0 0 2uo 5 dy_ 0 u2o g1 2uo 0 F12

(5.158)

The four eigenvalues of Axy are in Eq. (5.148). The eigenvalues of the second row in Eq. (5.157) are pﬃﬃﬃﬃﬃ l5;6 ¼ j g1 ¼ juz : (5.159) The free response of Eq. (5.157) in the (dx,dy) plane is unstable due to the positive eigenvalue l1 in Eq. (5.148). However, for particular choices of the initial conditions, the free response of Eq. (5.157) becomes bounded. To this aim, consider the modal decomposition of Axy: Axy ¼ VLxy U T ;

2

Lxy ¼ diagðl1 ; :::; l4 Þ;

V ¼ ½ v1

.

v4 ;

uT1

3

6 7 7 UT ¼ 6 4 « 5; uT4

(5.160)

where vi and li are the eigenvectors and eigenvalues of Axy, respectively, and ui is a left eigenvector of Axy. From linear algebra (see Section 2.3.1 and Section 13.2), the matrix exponential of Axy entering the free response can be written as expðAxy tÞ ¼ V expðLxy tÞU T .

(5.161)

Thus, the free response of the ﬁrst equation in Eq. (5.157) is given by drxy ðtÞ ¼ expðAxy tÞdrxy ð0Þ ¼ V expðLxy tÞUT drxy ð0Þ ¼

4 X i¼1

xi vi expðli tÞ;

(5.162)

234

5. PERTURBED ORBITAL DYNAMICS

where xi ¼ uTi drxy ð0Þ is a modal coefﬁcient. By choosing the initial condition drxy(0) in Eq. (5.158) such that x1 ¼ x2 ¼ 0, the free response (Eq. 5.162) simpliﬁes to include only oscillatory modes, and Eq. (5.162) becomes drxy ðtÞ ¼ x3 v3 expðjuxy tÞ þ x4 v4 expðjuxy tÞ ¼ 2ðReðx3 v3 Þcosðuxy tÞ Imðx3 v3 Þsinðuxy tÞÞ; (5.163) where uxy has been deﬁned in the second row of Eq. (5.148). Exercise 18 Prove that the initial conditions in Eq. (5.158) which force x1 ¼ x2 ¼ 0 in Eq. (5.162) are _ dxð0Þ ¼ k1 ðuxy Þuxy dyð0Þ ; (5.164) _ dyð0Þ ¼ kðuxy Þuxy dxð0Þ where kðuxy Þ ¼ ð2uo uxy Þ1 u2o g1 þ u2xy and ky2:91 for the EartheMoon system. HINT

With the eigenvalue notations of Eq. (5.148), prove that the left eigenvectors ui ; i ¼ 1; 2; of Axy hold 2 3 uo u2o þ 2g1 6 7 3 3 2 2 7 6 2 u11 u11 6 uo g1 u2o þ 2g1 w1 7 6u 7 6 u 7 7 pﬃﬃﬃﬃﬃﬃ 16 6 6 12 7 6 12 7 7 2 w1 (5.165) u1 ¼ 6 7 ¼ 6 7. 7; u2 ¼ 6 4 u13 5 4 u13 5 7 d6 pﬃﬃﬃﬃﬃﬃ 6 7 u0 w1 6 7 u14 u14 4 5 2 uo þ 2g1 w1 2 Then, solve the equations uTi drxy ð0Þ ¼ 0; i ¼ 1; 2. Symbolic tools may be of help. Prove also, with the help of Eqs. (5.163) and (5.164), that the free response of dx and dy takes the bounded form: cosðuxy tÞ sinðuxy tÞ dx dxð0Þ ðtÞ ¼ :, (5.166) sinðuxy tÞ cosðuxy tÞ k1 dyð0Þ k1 dy The free response of the second equation in Eq. (5.157) is given by dzðtÞ ¼ dzð0Þcosðuz tÞ þ

_ dzð0Þ sinðuz tÞ. uz

(5.167)

Lissajous Orbits By combining Eqs. (5.166) and (5.167), and by choosing the initial conditions dx(0) ¼ dz(0) ¼ 0 and dz_ ¼ uz yð0Þ, we obtain the following free response: dxðtÞ ¼ sinðuxy tÞk1 dyð0Þ dyðtÞ ¼ cosðuxy tÞdyð0Þ dzðtÞ ¼ sinðuz tÞdyð0Þ

;

(5.168)

5.6 RESTRICTED THREE-BODY PROBLEM

FIGURE 5.12

235

A Lissajous orbit around the EartheMoon L2 Lagrangian point.

which corresponds to the quasiperiodic orbits known as Lissajous orbits; an example is shown in Fig. 5.12. If the frequency ratio uxy/uz is not a rational number, Lissajous orbits are not closed. Halo Orbits If the frequency ratio uxy/uz is a rational number, Lissajous orbits become closed and periodic, and are known as Halo orbits; an example is shown in Fig. 5.13. Since in practical cases, the solution of Eq. (5.157) is not periodic, a feedback control is needed to force a rational frequency ratio and to achieve a periodic closed orbit, that is a Halo orbit. This kind of control is commonly known as frequency control or period control.

FIGURE 5.13

EartheMoon L2 Lagrangian point: a Halo orbit.

236

5. PERTURBED ORBITAL DYNAMICS

Lissajous and Halo orbits are solutions of the linearized Eq. (5.157), whereas the actual motion of a spacecraft near a Lagrange collinear point is determined by the nonlinear Eq. (5.127) and does not track Lissajous and Halo orbits. Despite this limitation, these orbits are important as they can be used as reference trajectories for several space missions, such as two recent European missions. LISA Pathﬁnder launched on 2 December 2015, operates (the nominal mission life is 11 months extendable to 17 months) in a Halo orbit around the SuneEarth Lagrangian point L1 [6]. GAIA, launched on 19 December 2013, operates (the nominal mission life is 5 years) in a Lissajous orbit around the SuneEarth Lagrangian point L2 [10]. Exercise 19 Consider the EartheMoon L2 Lagrangian point, for which we have uo ¼ 2.66 mrad/s, uxy/ uo ¼ 1.86, uz/uo ¼ 1.79, and k¼ 2.91 (see Exercise 18). Assume the initial condition y(0) ¼ 3500 km. Plot, from Eq. (5.168) and for t ¼ [0.16 2p/uo), the corresponding Lissajous orbit, which is shown in Fig. 5.12. Exercise 20 Consider the EartheMoon L2 Lagrangian point, for which we have uo ¼ 2.66 mrad/s, uxy/uo ¼ 1.86, and k ¼ 2.91 (see Exercise 18). Assume the initial condition y(0) ¼ 3500 km. Suppose that a frequency control is used, so that the out-of-plane frequency is changed to uz ¼ uxy. Plot, from Eq. (5.168) and for t ¼ [0.16 2p/uo), the corresponding Halo orbit, which is shown in Fig. 5.13.

References [1] E. Canuto, A. Molano -Jimenez, C. Perez-Montenegro, L. Massotti, Long-distance, drag-free, low-thrust, LEO formation control foe Earth gravity monitoring, Acta Astronautica 69 (2011) 571e582. [2] R. Fitzpatrick, An Introduction to Celestial Mechanics, Cambridge University Press, 2012. [3] H. Goldstein, Classical Mechanics, second ed., Addison-Wesley Pu. Co., Reading, MA, 1980. [4] M.H. Kaplan, Modern Spacecraft Dynamics & Control, John Wiley & Sons, 1976. [5] T.A. Lovell, S. Tragesser, Guidance for relative motion of low Earth orbit spacecraft based on relative orbit elements, in: AAS/AIAA Astrodynamics Specialist Conf., Providence, RI, 2004. AIAA 2004e4988. [6] P. McNamara, G. Racca, Introduction to LISA Pathﬁnder, in: ESA Doc. LISA-lpf-rp-0002, March 30, 2009. [7] A.E. Roy, Orbital Motion, fourth ed., IOP publishing, Bristol, 2005. [8] M.J. Sidi, Spacecraft Dynamics and Control, Cambridge University Press, 1997. [9] A. Sinclair, R.E. Sherrill, T.A. Lovell, Review of the solutions to the Tschauner-Hempel equations for satellite relative motion, in: 22nd AAS/AIAA Space Flight Mechanics Meeting, Charleston, South Carolina, 2012. AAS 12e149. [10] The GAIA collaboration team, The GAIA mission, Astronomy and Astrophysics 595 (A1) (2016) 1e36. [11] D.A. Vallado, Fundamentals of Astrodynamics and Applications, second ed., Microcosm Press, El Segundo, CA, 2001 (Kluwer Academic Pu., Dordrecht). [12] K.F. Wakker, Fundamentals of Astrodynamics, Institutional Repository Library, Delft University of Technology, Delft, The Netherlands, 2015. [13] B. Wie, Space Vehicle Dynamics and Control, in: AIAA Education Series, 1998.

Perturbed Orbital Dynamics

Perturbed Orbital Dynamics

Recommend Documents