Precise Analytical Computation of Frozen-Eccentricity, Low Earth Orbits in a Tesseral Potential

Abstract

Classical procedures for designing Earth’s mapping missions rely on a preliminary frozen-eccentricity orbit analysis. This initial exploration is based on the use of zonal gravitational models, which are frequently reduced to a simple 𝐽₂-𝐽₃ analysis. However, the 𝐽₂-𝐽₃ model may not be accurate enough for some applications. Furthermore, lower order truncations of the geopotential are known to fail in describing the behavior of elliptic frozen orbits properly. Inclusion of a higher degree geopotential, which also takes into account the short-period effects of tesseral harmonics, allows for the precise computation of frozen-eccentricity, low Earth orbits that show smaller long-period effects in long-term propagations than those obtained when using the zonal model design.

1. Introduction

The procedure of designing and maintaining mission orbits for artificial satellites is made of different phases that go from simple to complex. First of all it is essential to have a deep knowledge of the natural dynamics involved, which requires determining the origin and magnitude of the forces that prevent exact Keplerian motion to happen. The main forces are then used to create a model that, while being simple enough for providing analytical insight, must retain the bulk of the real dynamics. The simplified model allows for the selection of the range of trajectories useful for the mission, from which the nominal trajectory will be chosen and refined in a more sophisticated model. Finally, propagations in a real, ephemeris model are required for maneuver planning.

The selection of the simplified model is eased by the common practice of separating perturbations into secular, long-period, and short-period effects [1]. Then, long-term behavior is efficiently studied in mean elements, whose evolution equations are obtained by averaging the disturbing function over the fast rotating angles so as to remove short-period effects. Ideally, the long-term model would be solved analytically. However, as the analytical solution is rarely attained, one must be satisfied with taking advantage of the fair performances of semianalytical integration [2, 3].

However, the perturbation problem may admit Hamiltonian formulation in some cases, basically when the perturbative forces are of gravitational origin. Then, in analogy to the Hamilton-Jacobi method, which allows integrating a problem by finding the canonical transformation to action-angle variables [4], one can resort to classical perturbation theories where canonical transformations are used to construct a set of formal integrals, such that the new Hamiltonian only depends on actions—the momenta conjugate to the angles—up to the required truncation order [5].

Traditionally, perturbation approaches that proceed by averaging all the angles are used to obtain an explicit analytical time solution [6, 7]. However, averaging is only allowed when the involved angles rotate or “circulate,” as is the case of the mean anomaly. Hence, a full averaging process that removes both short- and long-period effects to obtain the so-called “proper” elements [8] may miss important pendular-type solutions related to the oscillation or “libration” of some angles, thus constraining the validity of the analytical solution to certain regions in the phase space. On the contrary, the averaging may only be applied to the circulating angles in order to reduce the Hamiltonian to one degree of freedom. Then, the Hamiltonian flow will be free of the constraints that apply to the full averaged problem. Even in the cases in which an explicit closed form solution to the reduced Hamiltonian is not known, the search for particular solutions of the reduced problem, and notably the stationary solutions, may provide important results.

A conspicuous example is the so-called Kozai mechanism [9], by means of which asteroids in high-inclination orbits may experience large orbital eccentricity and inclination variations under third-body perturbations. Kozai’s eccentricity-varying orbits stem from an unstable stationary solution of a reduced Hamiltonian, involving only the variation of the eccentricity and the argument of the periapsis (see also [10]). The Kozai mechanism is not only of interest for astronomers [11, 12], but also attracts the attention of aerospace engineers in view of its direct application to the case of science missions about natural satellites [13–18].

Similarly, orbits with stationary perigee and eccentricity, on average, are found to exist under the influence of the geopotential alone. These types of specialized orbits reduce stationkeeping requirements by using the natural dynamics. Furthermore, because almost constant eccentricity and fixed-perigee orbits minimize altitude variations, they are found quite useful as nominal orbits for Earth mapping missions [19]. In the usual aerospace engineering lingo, stationary solutions of these kinds are generally termed as frozen orbits and are loosely defined as orbits that, on average, maintain constant certain orbital design parameters, such as the inclination for sun-synchronous orbits, the node rotation rate as in the case of repeat-groundtrack orbits, or the eccentricity and argument of the perigee. The latter, which is the topic of this paper, are termed frozen-eccentricity orbits [20]. Historical details on the frozen orbit concept can be found in [21].

In the preliminary design of Earth mapping missions, the frozen-eccentricity orbit is usually computed in a 𝐽₂-𝐽₃ model [19], although the influence of higher order harmonics may be required in some applications [22, 23]. Specifically, it must be noted that second-order effects of 𝐽₂ and other zonal harmonics may introduce radical changes in the general frozen orbits picture when compared to the 𝐽₂-𝐽₃ case [21].

Because frozen orbits are determined in the mean elements space, whereas nominal orbits are realized by their osculating elements, the computation of precise initial conditions needs the readmission of the short-period effects that were removed in the averaging. This readmission is done by using the transformation from mean to osculating elements. In general, the frozen-eccentricity orbit design is efficiently carried out in a zonal model. However, in some cases the influence of tesserals may worsen the frozen-eccentricity condition considerably, thus making the effort of recovering short-period effects related only to zonals naive [24]. In the present paper we extend the definition of frozen-eccentricity orbits to include the short-period effects due to tesserals in the computation of initial conditions for the nominal orbits, in addition to the long- and short-period effects due to zonal harmonics.

The paper is organized as follows. First of all, several classical facts on perturbations of the orbital elements due to noncentralities of the geopotential are recalled in Section 2. Then, in Section 3 we describe the main lines in the computation of the perturbation theory, which consists of the long-period Hamiltonian and the transformation equations from mean to osculating elements, by means of the Lie transforms method [25–27]. The explicit time dependence introduced in the model by tesseral terms can be avoided by moving to a rotating frame. This change of the reference system requires the concomitant introduction of the Coriolis term into the Hamiltonian. In general, in order to make possible the elimination of short-period angles related to longitude-dependent terms from the Hamiltonian, a preliminary expansion of the true anomaly as a trigonometric series in the mean anomaly, whose coefficients are functions of the eccentricity, is required [28]. These kinds of expansions constrain the use of a perturbation theory to the case of small or moderate eccentricities [6]. However, in the case of close-Earth satellites, it is well known that the Coriolis term is of a higher order than the Keplerian [29]. This fact allows the short-period effects of the geopotential affecting low Earth orbits (LEO), including those related to longitude-dependent terms, to be averaged in closed form of the eccentricity, with a notable reduction in the number of literal terms required by the perturbation theory [7, 30–32]. Although this averaging is only valid away from tesseral resonances, deep tesseral resonances are not expected for the low-Earth orbits which are the topic of the present study.

The closed-form averaging is based on the computation of three canonical transformations of the Lie type. The standard elimination of the parallax [33] simplifies the original Hamiltonian by reducing the parallactic terms to the form 𝑟⁻², and by removing short-period terms related to the argument of the latitude. In the case of LEO orbits, a side effect of the standard elimination of the parallax is the removal of longitude-dependent terms, except for those related to the so-called 𝑚-dailies [34]. After this initial simplification, short-period terms depending on the mean anomaly are removed from the Hamiltonian via Delaunay normalization [35]. It follows a simple averaging of the remaining terms that depend on the argument of the node in the rotating frame. As expected from the known fact that the effects of tesseral average out in the long term, except for tesseral resonances, we find that the reduced Hamiltonian is exactly the same as the zonal problem [21].

Numerical experiments related to the computation of frozen-eccentricity orbits are presented in Section 4. Because the Lie transforms are carried out only up to second-order effects of 𝐽₂, the computed frozen orbits are affected by residual long-period effects derived from the truncation order of the perturbation theory. In the case of elliptic orbits, we show that short-period terms due to tesserals may worsen the frozen-eccentricity solution to some extent, thus magnifying the residual long-period effects. On the contrary, for a variety of cases tested, the short periods introduced by the tesserals do not seem to play any relevant role for almost circular frozen-eccentricity orbits. Therefore, the classical frozen orbits theory based on a zonal model would be enough for computing accurate initial osculating elements of LEO, frozen-eccentricity, almost circular orbits.

Finally, it must be mentioned that we limit the model to a fifth-degree and order truncation of the geopotential, whereas a ninth-degree model has been suggested to minimize variations of model predictions when compared with long-term observations of real satellites [21]. However, the fifth-degree and order model is simple enough to show all the relevant aspects of the frozen orbits picture and hence is considered sufficient for illustrating the theoretical facts discussed in the present research. On the contrary, higher order harmonics should be taken into account in the construction of a frozen-eccentricity orbit theory intended for practical applications.

2. Noncentralities of the Geopotential: Classical Facts

We only focus on the effects of the geopotential and assume that the Earth is in uniform rotation about the 𝑧-axis, without suffering from precessional or any other effects. Then, the Keplerian motion typical of mass-point attraction is disturbed by the noncentralities of the Earth’s gravitational field, which in an Earth-centered, Earth-fixed coordinate system, are derived from the disturbing potential [1] 𝜇𝑉=𝑟_𝑗≥2𝛼𝑟^𝑗𝐽_𝑗𝑃_𝑗,0−sin𝜑_𝑗𝑘=1𝐶_𝑗,𝑘cos𝑘𝜆+𝑆_𝑗,𝑘sin𝑘𝜆𝑃_𝑗,𝑘sin𝜑, (1) where 𝜇 is the Earth’s gravitational parameter, 𝛼 is the equatorial radius of the Earth, (𝑟,𝜑,𝜆) are spherical coordinates, with 𝑟 standing for geocentric radius, 𝜑 for latitude, and 𝜆 for longitude, 𝐽_𝑗=−𝐶_𝑗,0, 𝐶_𝑗,𝑘, and 𝑆_𝑗,𝑘 are harmonic coefficients, and 𝑃_𝑗,𝑘 are associated Legendre polynomials of degree 𝑗 and order 𝑘.

For orbital applications, it is customary to formulate the disturbing potential (1) in orbital elements (𝑎,𝑒,𝑖,𝜔,Ω, and 𝑀), respectively, standing for semimajor axis, eccentricity, inclination, argument of the perigee, longitude of the ascending node, and mean anomaly. The conversion from spherical coordinates to orbital elements is straightforward, based on the transformation from spherical to Cartesian coordinates in the Earth-centered Earth-inertial frame: 𝑥=𝑟cos𝜑cos(𝜆−𝜗),𝑦=𝑟cos𝜑sin(𝜆−𝜗),𝑧=𝑟sin𝜑, (2) where (𝑥,𝑦,𝑧) are usual Cartesian coordinates and 𝜗 is the Greenwich sidereal time, and the fact that the orbital frame, defined by the position vector 𝐫 and its derivative with respect to time 𝐫̇, and the Earth-centered, Earth-inertial frame are related by simple rotations: 𝑥𝑦𝑧=𝑅₃(−Ω)𝑅₁(−𝑖)𝑅₃()𝑟00−𝜃, (3) where 𝜃=𝑓+𝜔 is the argument of the latitude, 𝑓 is the true anomaly, 𝑟=‖𝐫‖, and 𝑅₁ and 𝑅₃, respectively, stand for standard rotation matrices about the 𝑥- and 𝑧- axes: 𝑅₁𝛽=,𝑅1000cos𝛽sin𝛽0−sin𝛽cos𝛽₃𝛽=.cos𝛽sin𝛽0−sin𝛽cos𝛽0001 (4)

The conversion from spherical coordinates to orbital elements shapes the disturbing potential in (1) as a trigonometric series whose arguments are linear, integer combinations of the angles 𝑓, 𝜔, Ω, and 𝜗, and whose coefficients depend on inclination and eccentricity. If, in addition, the true anomaly is expanded as a function of the mean anomaly, one gets [1] 𝜇𝑉=−𝑎^∞_𝑙=2𝛼𝑟^𝑗𝑗𝑗_𝑚=0𝑝=0𝐹_{𝑙,𝑚,𝑝}(𝑖)_{∞𝑞=−∞}𝐺_{𝑙,𝑝,𝑞}(𝑒)𝒬_{𝑙,𝑚,𝑝,𝑞}(𝑀,Ω−𝜗,𝜔), (5) where 𝑙−𝑚even⇒𝒬_{𝑙,𝑚,𝑝,𝑞}=𝐶_𝑙,𝑚cos𝜓_{𝑙,𝑚,𝑝,𝑞}+𝑆_𝑙,𝑚sin𝜓_{𝑙,𝑚,𝑝,𝑞},𝑙−𝑚odd⇒𝒬_{𝑙,𝑚,𝑝,𝑞}=−𝑆_𝑙,𝑚cos𝜓_{𝑙,𝑚,𝑝,𝑞}+𝐶_𝑙,𝑚sin𝜓_{𝑙,𝑚,𝑝,𝑞},𝜓_{𝑙,𝑚,𝑝,𝑞}=𝑙−2𝑝+𝑞𝑀+𝑚(Ω−𝜗)+𝑙−2𝑝𝜔. (6)(7) The formulation of the disturbing potential in orbital elements is efficiently carried out using recursion formulas for the computation of the inclination and eccentricity functions, respectively, 𝐹 and 𝐺 in (5), as is explained in full detail, for instance, in the popular book by Kaula [1]. Note that, in general, Kaula’s 𝐺 functions cannot be expressed in closed form of the eccentricity and, therefore, the use of (5) is limited to the case of small or moderate eccentricities.

In view of (6) and (7), the effects of the disturbing potential when expressed in orbital elements are generally classified in short-periodic effects, due to terms with 𝑙−2𝑝+𝑞≠0 which are related to the period of the “fast angle” 𝑀, long-periodic effects, due to terms with 𝑙−2𝑝+𝑞=𝑚=0 related to the period of the perigee rotation, and secular effects produced by terms with 𝑙=2𝑝 and 𝑞=𝑚=0. Specifically, even zonal harmonics contribute secular terms 𝜓_{𝑙,𝑚,𝑝,𝑞}=0, as shown in (6), while odd zonal harmonics introduce long-period terms.

The period of terms depending on Ω−𝜗 is close to one day and may be comparable to the orbital period in the case of high-altitude orbits. On the contrary, for LEO orbits their period is one order of magnitude higher than the orbital period and the inclusion of higher-order 𝑚-daily terms in the disturbing potential is occasionally recommended [34, 36]. Terms with 𝑙−2𝑝+𝑞=0 and 𝑚≠0 are sometimes classified as medium-period terms.

The above classification of the disturbing function in secular, long-period, and short-period terms is quite useful in studying orbit evolution in the long term. The amplitude of the short-periodic effects is normally small and orbit evolution is conveniently studied by retaining only the secular and long-periodic terms of the disturbing function (5). Except for the case of tesseral resonances 𝑗𝑀̇≈𝑘𝜗̇ with 𝑗 and 𝑘 integers, where such combinations of the indices that 𝑘(𝑙−2𝑝+𝑞)≈𝑗𝑚 makes 𝜓_{𝑙,𝑚,𝑝,𝑞} evolve slowly, therefore contributing long-period effects to the disturbing function, the long-term disturbing function is represented by the zonal harmonics (𝑚=0) contribution to (5), which is further “averaged” over the mean anomaly (𝑙−2𝑝+𝑞=0).

Lagrange planetary equations provide the rate of variations for the orbital elements, which require the computation of the partial derivatives of the disturbing function (5) with respect to the orbital elements themselves (see [5], for instance). It is easy to check that the partial derivatives of 𝑉 with respect to 𝑀, Ω, and 𝜔—appearing in the Lagrange planetary equations for 𝑎, 𝑒, and 𝑖—change sines into cosines in (7), and vice versa, thus preventing the appearance of secular terms in the semimajor axis, eccentricity, and inclination, as derived from (6). On the other hand, partial derivatives of 𝑉 with respect to 𝑎, 𝑒, and 𝑖—appearing in the Lagrange planetary equations for the mean anomaly, the longitude of the node, and the argument of the perigee—do not alter (7) and hence contribute secular terms to the rate of variation of 𝑀, Ω, and 𝜔, as derived from (6). When dealing only with long-period terms of the disturbing function, Lagrange planetary equations provide the evolution of a different set of elements, called mean orbital elements, which are sometimes written with primes to enhance the differences with respect to the instantaneous or osculating elements.

While the flow in mean elements is representative enough for long-term orbit evolution, the computation of initial conditions corresponding to particular solutions of the mean elements flow requires the knowledge of the transformation from osculating to mean elements: (𝑎,𝑒,𝑖,𝜔,Ω,𝑀)−−^𝒯→𝑎^′,𝑒^′,𝑖^′,𝜔^′,Ω^′,𝑀^′. (8) Besides, the long-term orbital evolution of some specific osculating elements can be efficiently explored by propagating corresponding mean elements, which in turn requires the knowledge of inverse transformation 𝒯⁻¹. As a closed form transformation is not available, one has to be satisfied with knowing some terms in the series expansion of the above transformation, which these days are customarily computed by the Lie transforms method [25–27].

The knowledge of the flow in mean elements alone may be enough in some engineering applications. But even in these cases, the simple removal of short-period terms 𝑗−2𝑙+𝑚≠0 from the geopotential [22, 23] misses important second-order effects of the second zonal harmonic, such that the averaging of short periods must be properly accomplished with the inclusion of second-order effects of 𝐶_2,0=−𝐽₂.

3. Long-Term Disturbing Function

As the perturbation problem admits Hamiltonian formulation, we perform the averaging using perturbation theory by Lie transforms [25–27], which, in addition, provides the explicit transformation from mean to osculating elements, which is required in the computation of the nominal orbit. Besides, in view of 𝜗 that always appears as a subtraction of the argument of the node, the explicit time dependence introduced by longitude-dependent terms, that is, tesseral and sectorial harmonics, is avoided by formulating the disturbing function in orbital elements in the Earth-centered Earth-fixed frame, which is a rotating frame attached to the Earth. In consequence, the Hamiltonian must be supplemented with the Coriolis term.

Thus, the perturbation problem is materialized by the initial Hamiltonian: 𝜇ℋ=−𝜗̇2𝑎1+2𝑛𝜂𝑐+𝑉, (9) where 𝜂=√1−𝑒², 𝑛=𝜇^1/2𝑎^−3/2 is the orbit mean motion, and 𝑐=cos𝑖, and orbital elements must be understood as functions of certain sets of canonical variables, such as polar-nodal or Delaunay variables. Bear in mind that the polar-nodal set (𝑟,𝜃,𝜈,𝑅,Θ,𝑁) is formed by the radial distance 𝑟, the argument of the latitude 𝜃, the argument of the node 𝜈, the radial velocity 𝑅, the modulus of the angular momentum vector Θ, and the projection of the angular momentum vector along the Earth’s rotation axis 𝑁. The set of Delaunay variables (ℓ,𝑔,ℎ,𝐿,𝐺,𝐻) is formed by the mean anomaly ℓ, the argument of the perigee 𝑔, the argument of the node ℎ, the Delaunay action 𝐿=√𝜇𝑎, the modulus of the angular momentum vector 𝐺=𝐿𝜂, and the component of the angular momentum vector along the Earth’s rotation axis 𝐻=𝐺𝑐. Remember also that the dynamical relation 0≤|𝐻|≤𝐺≤𝐿 holds. Besides, since we are using the rotating frame 𝜈=ℎ=Ω−𝜗.

We limit our model to a fifth-degree and order truncation of the geopotential 𝑉, which incorporates the fundamental qualitative dynamics of the real problem. This model is sufficient to illustrate the importance of taking short-period effects due to tesseral harmonics into account in the transformation equations from mean to osculating elements. The consideration of higher degrees and orders does not introduce relevant variations in our approach, except for the increase of the computational burden derived from the rapid growth in the size of the series used.

Given that we limit the theory only to the case of LEO orbits, deep tesseral resonances are not expected to happen and, therefore, the short-period effects to be removed by averaging are those related to the mean anomaly as well as to the argument of the node in the rotating frame. In addition, the ratio 𝜗̇/𝑛∼𝐽₂=𝜖 for LEO orbits, such that it is convenient to arrange the Hamiltonian in the following perturbative order: ℋ=⁴_𝑚=01ℋ𝑚!^𝑚, (10) where ℋ₀=12𝑅²+Θ²𝑟²−𝜇𝑟 (11) is the Keplerian term, ℋ₁=−𝑁𝜗̇ (12) is the Coriolis term, and ℋ_2,0 carries the contribution of 𝐶_2,0: ℋ₂𝜇=2!𝑟𝛼²𝑟²𝐶_2,0231−2𝑠²+32𝑠²cos2𝜃, (13) where 𝑠=√1−𝑐²=sin𝑖. Because the contribution of other harmonics of the geopotential is of the second order of 𝐶_2,0, which is 𝒪(𝜖⁴), there are no terms of the order of 𝜖³ in the original Hamiltonian; hence, we choose ℋ₃ℋ=0,₄1=4!𝑉−ℋ2!², (14) where, instead of using (5), the disturbing potential 𝑉=𝑉(𝑟,𝜃,𝜈,−,Θ,𝑁) is expressed in polar-nodal variables using (2) and (3), and taking into account that 𝜈=Ω−𝜗 and cos𝑖=𝑁/Θ.

Averaging is then made by Lie transforms up to the fourth order, based on the so-called “Lie triangle” inductive algorithm; a full account of the method can be consulted in the original reference [26]. For the sake of making the analytical manipulation more efficient, the perturbation theory is split into three Lie transforms, each of which is of a different nature, as sketched out in the following.

3.1. Elimination of the Parallax

First of all, we carry out a preparatory simplification (𝑟,𝜃,𝜈,𝑅,Θ,𝑁)−−→^𝒯1𝑟^′,𝜃^′,𝜈^′,𝑅^′,Θ^′,𝑁^′, (15) called the elimination of the parallax [33]. Originally designed for the simplification of zonal models, this transformation only removes part of the short-period terms, those related to the argument of the latitude, while retaining the contribution of parallactic terms of the form 1/𝑟². This transformation was very soon extended to consider also tesseral effects, although the efficiency of the “extended parallax transformation” [37] may be questioned, due to the scarce simplification of tesseral terms and the large size of the generating function. However, in those cases in which the Coriolis term may be considered as a higher order than the Keplerian, the homological equation is released from its dependence on the argument of the node; therefore, the standard elimination of the parallax can be applied, providing a notable simplification of tesseral terms [7, 30–32].

After applying the standard elimination of the parallax to the perturbed Hamiltonian (10), we find 𝒫=_𝑚≥01𝒫𝑚!^𝑚, (16) where, neglecting terms above the fourth order, we get 𝒫₀𝜇=−,𝒫2𝑎₁𝜇=−2𝜗̇2𝑎𝑛𝒫𝜂𝑐,₂𝜇=−𝛼2𝑎𝑟²²𝐶_2,01𝜂²1−3𝑐²,𝒫₃𝒫=0,₄=𝜇𝑎2𝑎𝑟²²𝛼⁴𝑎⁴𝐶²_2,02𝜂⁶𝑞^′_0,0+𝑞^′_2,0𝜂²+𝑞^′_2,2𝑒²+𝜇cos2𝜔𝑎2𝑎𝑟²²⁵_𝑘=2𝛼^𝑘𝑎^𝑘2𝜂^12𝑘−2_𝑚=0𝑌^𝑘_{𝑘,𝑚𝑗=−𝑘}𝑆_{𝑘,𝑚,𝑗}, (17)(18)(19)(20)(21) where 𝑆_{𝑘,𝑚,𝑗}=𝑄_𝑘,𝑗𝑞_{𝑘,𝑚,𝑗}cos(𝑘+2𝑚−4)𝜔+𝑗ℎ+𝑃_𝑘,𝑗𝑝_{𝑘,𝑚,𝑗}sin(𝑘+2𝑚−4)𝜔+𝑗ℎ (22) with 𝑘+𝑗even⇒𝑄_𝑘,𝑗=𝐶_𝑖,|𝑗|,𝑃_𝑘,𝑗=𝑆_𝑖,|𝑗|,𝑘+𝑗odd⇒𝑄_𝑘,𝑗=𝑆_𝑖,|𝑗|,𝑃_𝑘,𝑗=𝐶_𝑖,|𝑗|, (23) and the eccentricity polynomials 𝑌_𝑘,𝑚 are 𝑌_0,2=1,𝑌_0,3=𝑒,𝑌_0,4=2+3𝑒²,𝑌_0,5=4𝑒+3𝑒³,𝑌_1,𝑘=𝑒^𝑘−2. (24) Finally, the inclination polynomials 𝑞_𝑚,𝑗, 𝑝_{𝑘,𝑚,𝑗}, and 𝑞_{𝑘,𝑚,𝑗} are given in Tables 1, 2, and 3, respectively. It is worth noting that the orbital elements are now functions of the prime variables.

Table 1 

Inclination polynomials 𝑞^′_𝑚,𝑗 in (21) and 𝑑^′_𝑚,𝑗 in (28).

Table 2 

Inclination polynomials 𝑝_{𝑘,𝑚,𝑗} in (22). We abbreviate 𝛾≡1±𝑐.

Table 3 

Inclination polynomials 𝑞_{𝑘,𝑚,𝑗} in (22). We abbreviate 𝛾≡1±𝑐.

The Lie transforms method also provides the generating function of the transformation 𝒯₁, from which the explicit transformation equations 𝜉=𝜉𝑟^′,𝜃^′,𝜈^′,𝑅^′,Θ^′,𝑁^′,𝜉∈(𝑟,𝜃,𝜈,𝑅,Θ,𝑁) (25) are obtained by means of the straightforward application of the Lie triangle (see [26], for details).

3.2. Delaunay Normalization

First, the simplified Hamiltonian (16) is expressed in Delaunay elements. Then, the elimination of remaining short-period terms related to the mean anomaly is carried out by means of a new transformation: ℓ^′,𝑔^′,ℎ^′,𝐿^′,𝐺^′,𝐻−−→^′𝒯2ℓ^′′,𝑔^′′,ℎ^′′,𝐿^′′,𝐺^′′,𝐻^′′ (26) which, up to the truncation order, normalizes the Hamiltonian by introducing the Delaunay action as a formal integral of the perturbation problem [35]. The Delaunay normalization is carried out in closed form of the eccentricity based on the standard differential relations between the true and mean anomaly 𝑎²𝜂dℓ=𝑟²d𝑓, derived from the preservation of the angular momentum in the Kepler problem.

After the normalization, we obtain a new Hamiltonian: 𝒟=_𝑚≥01𝒟𝑚!^𝑚, (27) where, neglecting terms of orders higher than 4, we get 𝒟₁=𝒫₁,  𝒟₂=𝒫₂,  𝒟₃=𝒫₃, and 𝒟₄=𝜇𝛼2𝑎𝑎⁴⁴𝐶²_2,02𝜂⁷𝑑^′_0,0+𝑑^′_1,0𝜂+𝑑^′_2,0𝜂²+𝑑^′_2,2𝑒²+𝜇cos2𝜔2𝑎⁵_𝑘=2𝛼^𝑘𝑎^𝑘2𝜂^12𝑘−1_𝑚=0𝑌^𝑘_{𝑘,𝑚𝑗=−𝑘}𝑆_{𝑘,𝑚,𝑗}, (28) where 𝑌_𝑘,𝑚 and 𝑆_{𝑘,𝑚,𝑗} are the same as before, and the inclination coefficients 𝑑^′_𝑚,𝑗 are given in the right column of Table 1. Note that now the orbital elements and the different functions of them are expressed in the double-prime Delaunay variables.

Again, the Lie transforms method provides the generating function of the transformation 𝒯₂, from which the transformation 𝒯₂ is itself obtained by means of explicit equations.

3.3. Elimination of the Node

The Hamiltonian is now free from short-period effects related to the mean anomaly; however, it still depends on the argument of the node in the rotating frame. Therefore, we carry out one more Lie transform −,𝑔^′′,ℎ^′′,𝐿^′′,𝐺^′′,𝐻−−→^′′𝒯3−,𝑔^′′′,ℎ^′′′,𝐿^′′′,𝐺^′′′,𝐻^′′′ (29) in order to remove ℎ^′′′ from the Hamiltonian, and thus introducing 𝐻^′′′ as a new formal integral of the double-averaged problem.

The final Hamiltonian is 𝒦=_𝑚≥01𝒦𝑚!^𝑚, (30) where, up to the fourth order, 𝒦_𝑚=𝒟_𝑚 for 𝑚=0,…,3, and 𝒦₄=𝜇𝛼2𝑎𝑎⁴⁴𝐶²_2,02𝜂⁷𝑑^′_0,0+𝑑^′_1,0𝜂+𝑑^′_2,0𝜂²+𝑑^′_2,2𝑒²+𝜇cos2𝜔2𝑎⁵_𝑘=2𝛼^𝑘𝑎^𝑘2𝜂^12𝑘−1_𝑚=0𝑌_𝑘,𝑚𝑆_𝑘,𝑚,0, (31) where 𝑌_𝑘,𝑚, 𝑆_𝑘,𝑚,0, and the inclination polynomials 𝑑^′_𝑚,𝑗 are the same as before, and the orbital elements, as well as the different functions of them, are expressed now in the triple-prime Delaunay variables.

Except for tesseral resonances, the effects of tesseral perturbations average out in the long-term. Therefore, it is no surprise to find that the long term Hamiltonian (30) is exactly the same as the one provided in the appendix of [21] for a zonal model, up to the truncation degree of the geopotential used in this paper. This agreement is an additional, important check on the validity of the transformations carried out in this paper.

The long-term flow is derived from Hamilton equations of the final Hamiltonian (30). Because of the averaging, (30) is of the form 𝒦=𝒦−,𝑔^′′′,−,𝐿^′′′,𝐺^′′′,𝐻^′′′, (32) where 𝐿^′′′=𝐿^′′ and 𝐻^′′′ are constant. Therefore, the time evolution of the reduced phase space (𝑔^′′′,𝐺^′′′) decouples from the motion of ℓ^′′′ and ℎ^′′′. Thus, we first solve the one degree of freedom system 𝑔̇^′′′=d𝑔^′′′=d𝑡𝜕𝒦𝜕𝐺^′′′,𝐺̇^′′′=d𝐺^′′′d𝑡=−𝜕𝒦𝜕𝑔^′′′, (33) for the initial conditions 𝑔₀^′′′, 𝐺₀^′′′ to get 𝑔^′′′=𝑔^′′′𝑡;𝑔₀^′′′,𝐺₀^′′′;𝐿^′′′,𝐻^′′′,𝐺^′′′=𝐺^′′′𝑡;𝑔₀^′′′,𝐺₀^′′′;𝐿^′′′,𝐻^′′′. (34) Then, the evolution of ℓ^′′′ and ℎ^′′′ is obtained by the following quadratures: ℓ^′′′=ℓ₀^′′′+𝜕𝒦𝑔^′′′(𝑡),𝐺^′′′(𝑡);𝐿^′′′,𝐻^′′′𝜕𝐿^′′′ℎd𝑡,_′′′=ℎ₀^′′′+𝜕𝒦𝑔^′′′(𝑡),𝐺^′′′(𝑡);𝐿^′′′,𝐻^′′′𝜕𝐻^′′′d𝑡. (35)

Although the analytical solution (34) may be computed in some cases, it would require the use of special functions, which do not provide the desired insight into the long-term evolution. On the other hand, particular solutions of the reduced flow may be quite useful in specific satellite missions. This is exactly the case of the stationary solutions of (33), corresponding to the frozen-eccentricity orbits which are the topic of this research.

Finally, we note that the reduced Hamiltonian (30) can be written in terms of the mean orbital elements as 𝒦=𝒦(𝑔^′′′,𝐺^′′′;𝐿^′′′,𝐻^′′′)=𝒦(𝑒,𝜔;𝑎,𝜅), where 𝑒=𝑒(𝐺^′′′,𝐿^′′′), 𝜔=𝑔^′′′, and the two formal integrals of the reduced flow 𝐿^′′′ and 𝐻^′′′ have been replaced by 𝑎=𝐿^′′′2/𝜇 and 𝜅=𝜅(𝑒,𝑖). The latter is given by 𝐻𝜅=^′′′𝐿^′′′=1−𝑒²cos𝑖=cos𝑖_𝑐, (36) where 𝑖_𝑐=arccos(𝐻^′′′/𝐿^′′′) stands for the inclination of circular orbits, a case in which 𝐿=𝐺.

4. Frozen Orbits Computation: Numerical Experiments

Because the reduced Hamiltonian is of just one-degree-of-freedom in the modulus of the angular momentum and the argument of the perigee, then for each pair of the formal integrals (𝐿^′′′,𝐻^′′′) the reduced flow is comprised of closed curves and equilibria. Note that the cylindrical map (𝑔^′′′,𝐺^′′′) introduces singularities in the phase space representation, which should be properly discussed in the invariants of the problem [38]. For engineering applications, where rectilinear trajectories are not of practical interest, the reduced flow can also be studied in a sphere [39].

Alternatively, insightful information can be obtained directly in orbital variables by means of global inclination-eccentricity diagrams (𝑖,𝑒) and local eccentricity vector diagrams (𝑒cos𝜔,𝑒sin𝜔), which give a thorough picture of the solution space for a given semimajor axis 𝑎₀. Indeed, after the straightforward reformulation of the relevant equations in orbital elements and in view of the frozen-eccentricity orbits which appear in the “meridian” 𝑔=±𝜋/2, (𝑖,𝑒)-diagrams of frozen-eccentricity orbits can be obtained by the simple, inexpensive evaluation of the contour level 0=𝑔(̇𝑒,𝑖;𝜔=±𝜋/2;𝑎=𝑎₀) of (33) in the variation range 0≤𝑒<1, 0≤𝑖<𝜋 [40].

A sample inclination-eccentricity diagram of frozen-eccentricity, direct inclination orbits is provided in Figure 1 for a mean semimajor axis 𝑎₀=8000km, where each point of the respective curves represents a frozen-eccentricity orbit for different values of the (averaged) energy and 𝐻^′′′ (or 𝜅). An almost symmetric plot can be obtained analogously for the case of retrograde, frozen-eccentricity orbits. Note the agreement of Figure 1 with results in [21] except for quantitative variations due to the different truncation of the model and the different value selected for the mean semimajor axis.

(a)

(b)

Figure 1 

(a) Inclination-eccentricity diagram of frozen-eccentricity orbits of the 5×5 tesseral problem, for a mean semimajor axis of 8000km. (b) Detail in the vicinity of the critical inclination. Full lines are used for 𝜔=90deg., whereas dashed lines indicate 𝜔=270deg.

Thus, in Figure 1 we observe two different lines of frozen-eccentricity orbits. The first starts from an equatorial almost circular orbit and is comprised of low-eccentricity orbits with the perigee fixed, on average, at 90 degrees, which exist for increasing inclination. Close to the critical inclination 𝑖=arccos1/5 (due to the influence of higher order harmonics, the value of the critical inclination is slightly different from the value 𝑖=arccos1/5, which corresponds to the 𝐽₂ problem), the eccentricity of the frozen-eccentricity orbits increases with almost constant inclination. A new line starts from a frozen-eccentricity orbit close to the critical inclination, where 𝑒=0.023, 𝑖=63.45deg., and 𝜔=270deg. (the large red point in the bottom plot of Figure 1). One branch of this line is comprised of frozen-eccentricity orbits whose eccentricity increases with almost constant inclination, whereas the other branch is comprised of low-eccentricity orbits, which exist for increasing inclinations. For the latter, the eccentricity of the orbits reduces with inclination until finding an exact circular orbit at 𝑖=64.3533 deg. Then, for increasing inclinations, eccentricity becomes greater and the argument of the perigee changes to 90 deg. We do not attempt to find other frozen-eccentricity orbits that could exist with very high eccentricities.

Numerical computations based on a 5×5 truncation of the GRACE Gravity Model GGM02C (http://www.csr.utexas.edu/grace/gravity/), which uses the values 𝛼=6378.1363km, 𝜇=398600.4415 km³/s², and the unnormalized values of the harmonic coefficients, are presented in Table 4. Besides, for the Earth rotation rate, we used the value 𝜗̇=360.9856235 deg/day [41].

Table 4 

Unnormalized harmonic coefficients from GGM02.

On the other hand, local eccentricity vector diagrams can be constructed for each pair of the formal integrals (𝐿^′′′,𝐻^′′′), or, equivalently, (𝑎,𝜅), by means of simple contour plots of the reduced Hamiltonian (30). Eccentricity vector diagrams cannot show a clear representation of the reduced flow for orbits with high eccentricities. In spite of that, these kinds of diagrams succeed in presenting valuable information for usual aerospace engineering missions requiring low or moderate eccentricities, and hence their use is quite common in mission design and analysis procedures (see [42, 43], for instance). Applications in this paper are limited to 𝑒<0.2 and, therefore, corresponding eccentricity vector diagrams clearly show the frozen orbits’ behavior.

Thus, Figure 2 shows a collection of eccentricity vector diagrams that illustrate the typical bifurcation sequence of frozen-eccentricity orbits in the vicinity of the critical inclination. The mean semimajor axis, a dynamical parameter of the reduced flow, is fixed at 8000km, corresponding to circular orbits’ altitude of about 1 600 km over the Earth’s surface—about 400 km below the upper limit for LEO. The variation of the other dynamical parameter of the reduced flow 𝜅 produces the well-known sequence of bifurcations of frozen orbits in the vicinity of the critical inclination. Further details on the influence of different truncations of a zonal geopotential in the frozen-eccentricity orbits picture can be consulted in [21].

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2 

Bifurcation sequence of frozen-eccentricity orbits of the 5×5 tesseral problem, for a mean semimajor axis 𝑎=8000km. From left to right and from top to bottom, the inclination of the circular (nonfrozen) orbits is 𝑖=63.33, 63.43, 63.456, 63.46, 63.50, and 63.61 deg.

While local plots are useful in illustrating qualitative dynamical aspects of the problem and for estimating mean values of frozen-eccentricity orbits of interest, the computation of nominal orbits requires precise knowledge of corresponding mean elements. This is done by computing stationary solutions of the Hamilton equations of the reduced long-term Hamiltonian 𝒦=𝒦(𝑔,𝐺;𝐿,𝐻) in (33). Therefore, we must solve the system 𝐺̇^′′′=−𝜕𝒦𝜕𝑔^′′′=0,𝑔̇^′′′=𝜕𝒦𝜕𝐺^′′′=0. (37) Stationary solutions are conveniently found by means of standard root-finding procedures that can be fed by approximate values of fixed points on eccentricity vector diagrams, which are directly estimated from the figures.

Thus, for instance, in reference to the right-bottom plot of Figure 2, we find three fixed points. Two of them are of the elliptic type with 𝑒≈0.12 and 𝜔=90deg, 𝑒≈0.004 and 𝜔=270deg. The other, which happens in 𝑒≈0.11 and 𝜔=270deg, is of the hyperbolic type. Exact values fulfilling (37) give 𝑒=0.120130,𝜔=90^∘,𝑖=63^∘.4024,𝑒=0.00342451,𝜔=270^∘,𝑖=63^∘.6098, (38) for the elliptic type solutions, corresponding to stable frozen-eccentricity orbits, and 𝑒=0.113231,𝜔=270^∘,𝑖=63^∘.4258, (39) for the hyperbolic-type solution, corresponding to an unstable frozen-eccentricity orbit.

Even after the improvements, stationary solutions of (37) are just mean values. Therefore, direct propagation of these values in the nonaveraged problem will provide orbits that are only approximately frozen. This is illustrated in Figure 3 for the stable, elliptic frozen-eccentricity orbit, where we see the long-period oscillations of the argument of the perigee with an amplitude of about 30deg and a period of about one hundred years. Correspondingly, the osculating eccentricity is affected by nonnegligible oscillations between the nominal value and 𝑒≈0.128.

(a)

(b)

Figure 3 

Long-term evolution of the stable, elliptic frozen-eccentricity orbit in a 5×5 tesseral model when the osculating elements are taken the same as mean elements.

Classically, the nominal orbit is computed by recovering the short-period effects due to zonal harmonics, which are the only relevant perturbations in the long-term, assuming that there are no resonances. The readmission of zonal harmonics’ short-period effects in the computation of osculating elements, as derived from the generating function of each Lie transform, indeed improves the frozen orbit condition, as shown in Figure 4, where we note that now the amplitude of the long-term oscillations in the eccentricity is almost negligible when compared with the amplitude of the inevitable short-period oscillations. The amplitude of the long-term oscillations of the perigee is also notably reduced to less than 2 deg.

(a)

(b)

Figure 4 

Long-term evolution of the stable, elliptic frozen-eccentricity orbit in a 5×5 tesseral model when short-periodic effects from zonal harmonics are included in the computation of osculating elements.

As expected, the amplitude of the long-period oscillations can be further reduced when the computation of the osculating elements of the nominal orbit includes the readmission of short-period effects caused by tesseral harmonics. Thus, Figure 5 shows that now the amplitude of the oscillations in the argument of the perigee is further reduced to about one degree, while the long-period oscillations of the eccentricity are masked by the short periods, whereas the remaining long-period effects should be attributed to the truncation order of the perturbation theory, as illustrated, for instance, in [44].

(a)

(b)

Figure 5 

Long-term evolution of the stable, elliptic frozen-eccentricity orbit in a 5×5 tesseral model when short-periodic effects from zonal and tesseral harmonics are included in the computation of osculating elements.

The improvements in the frozen orbit computation are further illustrated in the eccentricity vector diagram of Figure 6, where we note that the readmission of short-period effects due to tesseral harmonics centers the osculating elements propagation over the predicted mean value.

Figure 6 

Eccentricity vector evolution of the orbits in Figure 5 (red points) and Figure 4 (blue points). Contour levels are mean elements.

The benefits of recovering short-periodic effects due to tesseral perturbations are not so evident in the case of the other two stationary solutions. Thus, as shown in the left plot of Figure 7, either taking or neglecting tesseral effects in the computation of the osculating orbit provides very good results, with slight improvements when recovering the tesseral effects (red plot) that provide an orbit better centered over the mean frozen value and with more regular shape. In the case of the unstable stationary solution (right plot of Figure 7), the osculating orbit derails over the unstable manifold in both cases, such that the corrections added by the tesserals simply change the unstable branch that the orbit follows.

(a)

(b)

Figure 7 

Eccentricity vector evolution of the frozen-eccentricity orbits when incorporating short-period effects from zonals (blue points) and also from tesserals (red points) superimposed to the eccentricity vector diagram in mean elements. Left: almost circular frozen-eccentricity orbit. Right: elliptic, unstable frozen-eccentricity orbit.

5. Conclusions

The classical engineering definition of frozen-eccentricity orbits as those orbits where secular effects due to even zonal harmonics are canceled by the long-period contribution of odd zonal harmonics, which has been more precisely defined in the past as stationary solutions of an averaged zonal model, is now amended to include the short-period effects from tesseral harmonics. Stationary solutions of the averaged model always carry residual long-period effects derived from the truncation order of the perturbation theory. But these undesired long-period effects can be minimized by recovering the short-period effects due to the tesseral harmonics in the computation of the osculating elements of the nominal frozen-eccentricity orbit. This theoretical result has been demonstrated for a five-degree and order truncation geopotential, which shares the main characteristics of the full model. To this effect, we have constructed a perturbation theory that considers up to second-order 𝐽₂ effects. In the case of LEO orbits, the perturbation theory can be constructed in closed form of the eccentricity, thus notably reducing the number of terms required in the formal series, and, therefore, making their evaluation more efficient. The perturbation theory, of course, is not limited to the computation of frozen-eccentricity orbits and can be applied to the fast and efficient semianalytical integration of LEO orbits. Thus, long-term evolution could be used to explain some observable drift in the eccentricity/argument of perigee of debris orbiting in special conditions.

We approached only the mathematical problem of the computation of frozen-eccentricity orbits under the noncentralities of the geopotential. The real, engineering problem is more complex, because second-order effects due to the non-fixed spin axis, lunisolar attraction, solar radiation pressure, or drag may destroy the frozen orbit condition, thus making it necessary to solve an optimization problem. Nevertheless, the proposed solution of this research will be an excellent starting solution in the optimization process.

Acknowledgment

Part of this research has been supported by the Government of Spain (Projects AYA 2009-11896, AYA 2010-18796, and Grant Gobierno de La Rioja Fomenta 2010/16).