Mathematical Methods Applied to the Celestial Mechanics of Artificial Satellites 2013
View this Special IssueResearch Article  Open Access
Martin Lara, Juan F. SanJuan, Luis M. LópezOchoa, "Precise Analytical Computation of FrozenEccentricity, Low Earth Orbits in a Tesseral Potential", Mathematical Problems in Engineering, vol. 2013, Article ID 191384, 13 pages, 2013. https://doi.org/10.1155/2013/191384
Precise Analytical Computation of FrozenEccentricity, Low Earth Orbits in a Tesseral Potential
Abstract
Classical procedures for designing Earth’s mapping missions rely on a preliminary frozeneccentricity 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 shortperiod effects of tesseral harmonics, allows for the precise computation of frozeneccentricity, low Earth orbits that show smaller longperiod effects in longterm 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, longperiod, and shortperiod effects [1]. Then, longterm 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 shortperiod effects. Ideally, the longterm 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 HamiltonJacobi method, which allows integrating a problem by finding the canonical transformation to actionangle 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 longperiod effects to obtain the socalled “proper” elements [8] may miss important pendulartype 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 socalled Kozai mechanism [9], by means of which asteroids in highinclination orbits may experience large orbital eccentricity and inclination variations under thirdbody perturbations. Kozai’s eccentricityvarying 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 fixedperigee 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 sunsynchronous orbits, the node rotation rate as in the case of repeatgroundtrack orbits, or the eccentricity and argument of the perigee. The latter, which is the topic of this paper, are termed frozeneccentricity orbits [20]. Historical details on the frozen orbit concept can be found in [21].
In the preliminary design of Earth mapping missions, the frozeneccentricity 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 secondorder 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 shortperiod effects that were removed in the averaging. This readmission is done by using the transformation from mean to osculating elements. In general, the frozeneccentricity orbit design is efficiently carried out in a zonal model. However, in some cases the influence of tesserals may worsen the frozeneccentricity condition considerably, thus making the effort of recovering shortperiod effects related only to zonals naive [24]. In the present paper we extend the definition of frozeneccentricity orbits to include the shortperiod effects due to tesserals in the computation of initial conditions for the nominal orbits, in addition to the long and shortperiod 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 longperiod 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 shortperiod angles related to longitudedependent 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 closeEarth satellites, it is well known that the Coriolis term is of a higher order than the Keplerian [29]. This fact allows the shortperiod effects of the geopotential affecting low Earth orbits (LEO), including those related to longitudedependent 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 lowEarth orbits which are the topic of the present study.
The closedform 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 shortperiod 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 longitudedependent terms, except for those related to the socalled dailies [34]. After this initial simplification, shortperiod 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 frozeneccentricity orbits are presented in Section 4. Because the Lie transforms are carried out only up to secondorder effects of , the computed frozen orbits are affected by residual longperiod effects derived from the truncation order of the perturbation theory. In the case of elliptic orbits, we show that shortperiod terms due to tesserals may worsen the frozeneccentricity solution to some extent, thus magnifying the residual longperiod 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 frozeneccentricity orbits. Therefore, the classical frozen orbits theory based on a zonal model would be enough for computing accurate initial osculating elements of LEO, frozeneccentricity, almost circular orbits.
Finally, it must be mentioned that we limit the model to a fifthdegree and order truncation of the geopotential, whereas a ninthdegree model has been suggested to minimize variations of model predictions when compared with longterm observations of real satellites [21]. However, the fifthdegree 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 frozeneccentricity 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 masspoint attraction is disturbed by the noncentralities of the Earth’s gravitational field, which in an Earthcentered, Earthfixed coordinate system, are derived from the disturbing potential [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, , , 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 Earthcentered Earthinertial frame: 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 Earthcentered, Earthinertial frame are related by simple rotations: where is the argument of the latitude, is the true anomaly, , and and , respectively, stand for standard rotation matrices about the  and  axes:
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] where 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 shortperiodic effects, due to terms with which are related to the period of the “fast angle” , longperiodic effects, due to terms with related to the period of the perigee rotation, and secular effects produced by terms with and . Specifically, even zonal harmonics contribute secular terms , as shown in (6), while odd zonal harmonics introduce longperiod terms.
The period of terms depending on is close to one day and may be comparable to the orbital period in the case of highaltitude orbits. On the contrary, for LEO orbits their period is one order of magnitude higher than the orbital period and the inclusion of higherorder daily terms in the disturbing potential is occasionally recommended [34, 36]. Terms with and are sometimes classified as mediumperiod terms.
The above classification of the disturbing function in secular, longperiod, and shortperiod terms is quite useful in studying orbit evolution in the long term. The amplitude of the shortperiodic effects is normally small and orbit evolution is conveniently studied by retaining only the secular and longperiodic terms of the disturbing function (5). Except for the case of tesseral resonances with and integers, where such combinations of the indices that makes evolve slowly, therefore contributing longperiod effects to the disturbing function, the longterm disturbing function is represented by the zonal harmonics () contribution to (5), which is further “averaged” over the mean anomaly ().
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 longperiod 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 longterm 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: Besides, the longterm 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 shortperiod terms from the geopotential [22, 23] misses important secondorder effects of the second zonal harmonic, such that the averaging of short periods must be properly accomplished with the inclusion of secondorder effects of .
3. LongTerm 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 longitudedependent terms, that is, tesseral and sectorial harmonics, is avoided by formulating the disturbing function in orbital elements in the Earthcentered Earthfixed 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: where , is the orbit mean motion, and , and orbital elements must be understood as functions of certain sets of canonical variables, such as polarnodal or Delaunay variables. Bear in mind that the polarnodal 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 holds. Besides, since we are using the rotating frame .
We limit our model to a fifthdegree 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 shortperiod 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 shortperiod 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: where is the Keplerian term, is the Coriolis term, and carries the contribution of : where . Because the contribution of other harmonics of the geopotential is of the second order of , which is , there are no terms of the order of in the original Hamiltonian; hence, we choose where, instead of using (5), the disturbing potential is expressed in polarnodal variables using (2) and (3), and taking into account that and .
Averaging is then made by Lie transforms up to the fourth order, based on the socalled “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 called the elimination of the parallax [33]. Originally designed for the simplification of zonal models, this transformation only removes part of the shortperiod terms, those related to the argument of the latitude, while retaining the contribution of parallactic terms of the form . 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 where, neglecting terms above the fourth order, we get where with and the eccentricity polynomials are 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.
The Lie transforms method also provides the generating function of the transformation , from which the explicit transformation equations 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 shortperiod terms related to the mean anomaly is carried out by means of a new transformation: 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 , derived from the preservation of the angular momentum in the Kepler problem.
After the normalization, we obtain a new Hamiltonian: where, neglecting terms of orders higher than 4, we get , , , and 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 doubleprime 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 shortperiod 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 in order to remove from the Hamiltonian, and thus introducing as a new formal integral of the doubleaveraged problem.
The final Hamiltonian is where, up to the fourth order, for , and where , , 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 tripleprime Delaunay variables.
Except for tesseral resonances, the effects of tesseral perturbations average out in the longterm. 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 longterm flow is derived from Hamilton equations of the final Hamiltonian (30). Because of the averaging, (30) is of the form 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 for the initial conditions , to get Then, the evolution of and is obtained by the following quadratures:
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 longterm 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 frozeneccentricity 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 and . The latter is given by where stands for the inclination of circular orbits, a case in which .
4. Frozen Orbits Computation: Numerical Experiments
Because the reduced Hamiltonian is of just onedegreeoffreedom 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 inclinationeccentricity diagrams and local eccentricity vector diagrams , 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 frozeneccentricity orbits which appear in the “meridian” , diagrams of frozeneccentricity orbits can be obtained by the simple, inexpensive evaluation of the contour level of (33) in the variation range , [40].
A sample inclinationeccentricity diagram of frozeneccentricity, direct inclination orbits is provided in Figure 1 for a mean semimajor axis km, where each point of the respective curves represents a frozeneccentricity orbit for different values of the (averaged) energy and (or ). An almost symmetric plot can be obtained analogously for the case of retrograde, frozeneccentricity 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)
Thus, in Figure 1 we observe two different lines of frozeneccentricity orbits. The first starts from an equatorial almost circular orbit and is comprised of loweccentricity orbits with the perigee fixed, on average, at degrees, which exist for increasing inclination. Close to the critical inclination (due to the influence of higher order harmonics, the value of the critical inclination is slightly different from the value , which corresponds to the problem), the eccentricity of the frozeneccentricity orbits increases with almost constant inclination. A new line starts from a frozeneccentricity orbit close to the critical inclination, where , deg., and deg. (the large red point in the bottom plot of Figure 1). One branch of this line is comprised of frozeneccentricity orbits whose eccentricity increases with almost constant inclination, whereas the other branch is comprised of loweccentricity orbits, which exist for increasing inclinations. For the latter, the eccentricity of the orbits reduces with inclination until finding an exact circular orbit at deg. Then, for increasing inclinations, eccentricity becomes greater and the argument of the perigee changes to deg. We do not attempt to find other frozeneccentricity orbits that could exist with very high eccentricities.
Numerical computations based on a truncation of the GRACE Gravity Model GGM02C (http://www.csr.utexas.edu/grace/gravity/), which uses the values km, km^{3}/s^{2}, and the unnormalized values of the harmonic coefficients, are presented in Table 4. Besides, for the Earth rotation rate, we used the value deg/day [41].
