Frozen orbit

Classical theory

The classical theory of frozen orbits is essentially based on the analytical perturbation analysis for artificial satellites of Dirk Brouwer made under contract with NASA and published in 1959.

This analysis can be carried out as follows:

In the article orbital perturbation analysis, the secular perturbation of the orbital pole ${\displaystyle \mathrm{\Delta }\hat{z}}$ from the ${\displaystyle J_2}$ term of the geopotential model is shown to be

which can be expressed in terms of orbital elements thus:

Making a similar analysis for the ${\displaystyle J_3}$ term (corresponding to the fact that the earth is slightly pear shaped), one gets

which can be expressed in terms of orbital elements as

In the same article the secular perturbation of the components of the eccentricity vector caused by the ${\displaystyle J_2}$ is shown to be:

where:

• The first term is the in-plane perturbation of the eccentricity vector caused by the in-plane component of the perturbing force
• The second term is the effect of the new position of the ascending node in the new orbital plane, the orbital plane being perturbed by the out-of-plane force component

Making the analysis for the ${\displaystyle J_3}$ term one gets for the first term, i.e. for the perturbation of the eccentricity vector from the in-plane force component

For inclinations in the range 97.8–99.0 deg, the ${\displaystyle \mathrm{\Delta }\mathrm{\Omega }}$ value given by (6) is much smaller than the value given by (3) and can be ignored. Similarly the quadratic terms of the eccentricity vector components in (8) can be ignored for almost circular orbits, i.e. (8) can be approximated with

Adding the ${\displaystyle J_3}$ contribution ${\displaystyle 2\pi \frac{J_3}{\mu p^3}\frac{3}{2}\sin i\left(\frac{5}{4}{\sin }^{2}i-1\right)(1,0)}$

to (7) one gets

Now the difference equation shows that the eccentricity vector will describe a circle centered at the point ${\displaystyle \left(0,-\frac{J_3\sin i}{J_{2}2p}\right)}$; the polar argument of the eccentricity vector increases with ${\displaystyle -2\pi \frac{J_2}{\mu p^2}3\left(\frac{5}{4}{\sin }^{2}i-1\right)}$ radians between consecutive orbits.

As

${\displaystyle \mu =398600.440{\text{ km}}^3/s^2}$

${\displaystyle J_2=1.7555{10}^{10}{\text{ km}}^5/s^2}$

${\displaystyle J_3=-2.619{10}^{11}{\text{ km}}^6/s^2}$

one gets for a polar orbit (${\displaystyle i={90}^{\circ }}$) with ${\displaystyle p=7200\text{ km}}$ that the centre of the circle is at ${\displaystyle (0,0.001036)}$ and the change of polar argument is 0.00400 radians per orbit.

The latter figure means that the eccentricity vector will have described a full circle in 1569 orbits. Selecting the initial mean eccentricity vector as ${\displaystyle (0,0.001036)}$ the mean eccentricity vector will stay constant for successive orbits, i.e. the orbit is frozen because the secular perturbations of the ${\displaystyle J_2}$ term given by (7) and of the ${\displaystyle J_3}$ term given by (9) cancel out.

In terms of classical orbital elements, this means that a frozen orbit should have the following mean elements:

${\displaystyle e=-\frac{J_3\sin i}{J_{2}2p}}$

${\displaystyle \omega ={90}^{\circ }}$

Modern theory

The modern theory of frozen orbits is based on the algorithm given in a 1989 article by Mats Rosengren.

For this the analytical expression (7) is used to iteratively update the initial (mean) eccentricity vector to obtain that the (mean) eccentricity vector several orbits later computed by the precise numerical propagation takes precisely the same value. In this way the secular perturbation of the eccentricity vector caused by the ${\displaystyle J_2}$ term is used to counteract all secular perturbations, not only those (dominating) caused by the ${\displaystyle J_3}$ term. One such additional secular perturbation that in this way can be compensated for is the one caused by the solar radiation pressure, this perturbation is discussed in the article "Orbital perturbation analysis (spacecraft)".

Applying this algorithm for the case discussed above, i.e. a polar orbit (${\displaystyle i={90}^{\circ }}$) with ${\displaystyle p=7200\text{ km}}$ ignoring all perturbing forces other than the ${\displaystyle J_2}$ and the ${\displaystyle J_3}$ forces for the numerical propagation one gets exactly the same optimal average eccentricity vector as with the "classical theory", i.e. ${\displaystyle (0,0.001036)}$.

When we also include the forces due to the higher zonal terms the optimal value changes to ${\displaystyle (0,0.001285)}$.

Assuming in addition a reasonable solar pressure (a "cross-sectional-area" of 0.05 m²/kg, the direction to the sun in the direction towards the ascending node) the optimal value for the average eccentricity vector becomes ${\displaystyle (0.000069,0.001285)}$ which corresponds to :${\displaystyle \omega ={87}^{\circ }}$, i.e. the optimal value is not ${\displaystyle \omega ={90}^{\circ }}$ anymore.

This algorithm is implemented in the orbit control software used for the Earth observation satellites ERS-1, ERS-2 and Envisat.

Derivation of the closed form expressions for the J3 perturbation

The main perturbing force to be counteracted in order to have a frozen orbit is the "${\displaystyle J_3}$ force", i.e. the gravitational force caused by an imperfect symmetry north–south of the Earth, and the "classical theory" is based on the closed form expression for this "${\displaystyle J_3}$ perturbation". With the "modern theory" this explicit closed form expression is not directly used but it is certainly still worthwhile to derive it. The derivation of this expression can be done as follows:

The potential from a zonal term is rotational symmetric around the polar axis of the Earth and corresponding force is entirely in a longitudinal plane with one component ${\displaystyle F_r\hat{r}}$ in the radial direction and one component ${\displaystyle F_{\lambda }\hat{\lambda }}$ with the unit vector ${\displaystyle \hat{\lambda }}$ orthogonal to the radial direction towards north. These directions ${\displaystyle \hat{r}}$ and ${\displaystyle \hat{\lambda }}$ are illustrated in Figure 1.

In the article Geopotential model it is shown that these force components caused by the ${\displaystyle J_3}$ term are

To be able to apply relations derived in the article Orbital perturbation analysis (spacecraft) the force component ${\displaystyle F_{\lambda }\hat{\lambda }}$ must be split into two orthogonal components ${\displaystyle F_t\hat{t}}$ and ${\displaystyle F_z\hat{z}}$ as illustrated in figure 2

Let ${\displaystyle \hat{a},\hat{b},\hat{n}}$ make up a rectangular coordinate system with origin in the center of the Earth (in the center of the Reference ellipsoid) such that ${\displaystyle \hat{n}}$ points in the direction north and such that ${\displaystyle \hat{a},\hat{b}}$ are in the equatorial plane of the Earth with ${\displaystyle \hat{a}}$ pointing towards the ascending node, i.e. towards the blue point of Figure 2.

The components of the unit vectors

${\displaystyle \hat{r},\hat{t},\hat{z}}$

making up the local coordinate system (of which ${\displaystyle \hat{t},\hat{z},}$ are illustrated in figure 2), and expressing their relation with ${\displaystyle \hat{a},\hat{b},\hat{n}}$, are as follows:

${\displaystyle r_a=\cos u}$

${\displaystyle r_b=\cos i\sin u}$

${\displaystyle r_n=\sin i\sin u}$

${\displaystyle t_a=-\sin u}$

${\displaystyle t_b=\cos i\cos u}$

${\displaystyle t_n=\sin i\cos u}$

${\displaystyle z_a=0}$

${\displaystyle z_b=-\sin i}$

${\displaystyle z_n=\cos i}$

where ${\displaystyle u}$ is the polar argument of ${\displaystyle \hat{r}}$ relative the orthogonal unit vectors ${\displaystyle \hat{g}=\hat{a}}$ and ${\displaystyle \hat{h}=\cos i\hat{b}+\sin i\hat{n}}$ in the orbital plane

Firstly

${\displaystyle \sin \lambda =r_n=\sin i\sin u}$

where ${\displaystyle \lambda }$ is the angle between the equator plane and ${\displaystyle \hat{r}}$ (between the green points of figure 2) and from equation (12) of the article Geopotential model one therefore obtains

Secondly the projection of direction north, ${\displaystyle \hat{n}}$, on the plane spanned by ${\displaystyle \hat{t},\hat{z},}$ is

${\displaystyle \sin i\cos u\hat{t}+\cos i\hat{z}}$

and this projection is

${\displaystyle \cos \lambda \hat{\lambda }}$

where ${\displaystyle \hat{\lambda }}$ is the unit vector ${\displaystyle \hat{\lambda }}$ orthogonal to the radial direction towards north illustrated in figure 1.

From equation (11) we see that

${\displaystyle F_{\lambda }\hat{\lambda }=-J_3\frac{1}{r^5}\frac{3}{2}\left(5{\sin }^2\lambda -1\right)\cos \lambda \hat{\lambda }=-J_3\frac{1}{r^5}\frac{3}{2}\left(5{\sin }^2\lambda -1\right)(\sin i\cos u\hat{t}+\cos i\hat{z})}$

and therefore:

In the article Orbital perturbation analysis (spacecraft) it is further shown that the secular perturbation of the orbital pole ${\displaystyle \hat{z}}$ is

Introducing the expression for ${\displaystyle F_z}$ of (14) in (15) one gets

The fraction ${\displaystyle \frac{p}{r}}$ is

${\displaystyle \frac{p}{r}=1+e\cdot \cos \theta =1+e_g\cdot \cos u+e_h\cdot \sin u}$

where

${\displaystyle e_g=e\cos \omega }$

${\displaystyle e_h=e\sin \omega }$

are the components of the eccentricity vector in the ${\displaystyle \hat{g},\hat{h}}$ coordinate system.

As all integrals of type

${\displaystyle \underset{0}{\overset{2\pi }{\int }}{\cos }^{m}u{\sin }^{n}udu}$

are zero if not both ${\displaystyle n}$ and ${\displaystyle m}$ are even, we see that

and

It follows that

where

${\displaystyle \hat{g}}$ and ${\displaystyle \hat{h}}$ are the base vectors of the rectangular coordinate system in the plane of the reference Kepler orbit with ${\displaystyle \hat{g}}$ in the equatorial plane towards the ascending node and ${\displaystyle u}$ is the polar argument relative this equatorial coordinate system

${\displaystyle f_z}$ is the force component (per unit mass) in the direction of the orbit pole ${\displaystyle \hat{z}}$

In the article Orbital perturbation analysis (spacecraft) it is shown that the secular perturbation of the eccentricity vector is

where

• ${\displaystyle \hat{r},\hat{t}}$ is the usual local coordinate system with unit vector ${\displaystyle \hat{r}}$ directed away from the Earth
• ${\displaystyle V_r=\sqrt{\frac{\mu }{p}}\cdot e\cdot \sin \theta }$ - the velocity component in direction ${\displaystyle \hat{r}}$
• ${\displaystyle V_t=\sqrt{\frac{\mu }{p}}\cdot (1+e\cdot \cos \theta )}$ - the velocity component in direction ${\displaystyle \hat{t}}$

Introducing the expression for ${\displaystyle F_r,F_t}$ of (12) and (13) in (20) one gets

Using that

${\displaystyle \frac{V_r}{V_t}=\frac{e_g\cdot \sin u-e_h\cdot \cos u}{\frac{p}{r}}}$

the integral above can be split in 8 terms:

Given that

${\displaystyle \hat{r}=\cos u\hat{g}+\sin u\hat{h}}$

${\displaystyle \hat{t}=-\sin u\hat{g}+\cos u\hat{h}}$

we obtain

${\displaystyle \frac{p}{r}=1+e\cdot \cos \theta =1+e_g\cdot \cos u+e_h\cdot \sin u}$

and that all integrals of type

${\displaystyle \underset{0}{\overset{2\pi }{\int }}{\cos }^{m}u{\sin }^{n}udu}$

are zero if not both ${\displaystyle n}$ and ${\displaystyle m}$ are even:

Term 1

Term 2

Term 3

Term 4

Term 5

Term 6

Term 7

Term 8

As

It follows that