ABSTRACT

The relative motion between multiple satellites is a developed technique with many applications. Formation-flying missions use the relative motion dynamics in their design. In this work, the motion in invariant relative orbits is considered under the effects of second-order zonal harmonics in an equatorial orbit. The Hamiltonian framework is used to formulate the problem. All the possible conditions of the invariant relative motion are obtained with different inclinations of the follower satellite orbits. These second-order conditions warrantee the drift rates keeping two, or more, neighboring orbits from drifting apart. The conditions have been modeled. All the possibilities of choosing mean elements of the leader satellite orbit and differences in momenta between leader and follower satellites’ orbits are presented.

Keywords:

Invariant Relative Orbits, Formation Flying Satellites, Relative Motion

1. Introduction

As the geostationary Earth orbits (GEO) belt becomes more crowded it is increasingly difficult to acquire slots for new satellites. Consequently, many organizations choose to collocate their spacecraft in the same slot. Also for missions which a single satellite cannot accomplish, as global position satellite system (GPS), the needed of formation flight began.

The formation flight concept is the use of several small satellites, which work together in a group (or constellation) to accomplish the objective of one larger, usually more expensive, satellite. This increases the likelihood of mission success in the event of a malfunction Hughes [1] . Formation flights have invariant relative orbits for their satellites to ensure that they will not separate over time.

The invariant relative orbits have been studied for a long time, as the earlier work of Clohessy and Wiltshire [2] in addition to the studies of Tschauner and Hempel [3] . These models introduced conditions on the initial relative position and velocity so that the relative orbits result to be periodic, which are closed orbits. Recently, Schaub and Alfriend [4] , Abd El-Salam et al. [5] passing through Li and Li [6] until Abd El-Salam and El-Saftawy [7] in which they discussed the invariant relative orbits due to the influence of the perturbative effects of the asphericity of the Earth, the relativistic corrections and the direct solar radiation pressure. Rahoma [8] also, discussed the J₂ invariant relative orbits with the effect of lunisolar attraction.

In this paper, we extend the works of Schaub and Alfriend [2] and Abd El-Salam et al. [5] model by introducing an atlas for the curves of invariant relative orbits’ conditions. This atlas will be presented using Mathematica program to calculate and plot graphics of the initial conditions of invariant relative orbits. Those graphics will be shown as curves in 2D; in the case of the orbit of the leader satellite is equatorial.

2. Hamiltonian System

There are several ways to derive the equations of motion for any such system. We emphasized on the Hamiltonian structure of this system. The Hamiltonian formulation allows additional conservative forces to add to the Hamiltonian, thus the addition of complexity to the model can be incorporated with ease. Non-conservative forces can add in the momenta equations of motion. The Hamiltonian equations of motion allow us to directly use control and simulation techniques.

After expressing the Hamiltonian, as a series in power of J₂ (The second geopotential zonal harmonic) up to the second order, and using Lie-Deprit-Kamel perturbation method Kamel [9] to eliminate, in successive, the short and long periodic terms, the transformed Hamiltonian, H^**, for different orders 0, 1, and 2, are obtained by El-Saftawy et al. [10] .

$H^{**}=H_0^{**}+J_2{H}_1^{**}+\frac{J_2^2}{2}{H}_2^{**}$ (1)

where,

$H_0^{**}=-\frac{{\mu }^2}{2}{\eta }_{2,0}$ (1.1)

$H_1^{**}=A_{11}{\eta }_{3,3}\left(3{\mathbb{S}}^2-2\right)$ (1.2)

$\begin{array}{c}{H}_2^{**}=\frac{1}{128}\frac{A_{11}^2}{{\mu }^2}[\left(-5517{\eta }_{1,9}+246{\eta }_{3,7}-3456{\eta }_{4,6}+135{\eta }_{5,5}\right){\mathbb{S}}^4\\ +\left(11520{\eta }_{1,9}-2976{\eta }_{3,7}+8064{\eta }_{4,6}+672{\eta }_{5,5}\right){\mathbb{S}}^2\\ +\left(-4704{\eta }_{1,9}+576{\eta }_{3,7}-1536{\eta }_{4,6}+288{\eta }_{5,5}\right)]\\ +\frac{3}{2}{A}_{22}\left(5{\eta }_{3,7}-3{\eta }_{5,5}\right)\left(35{\mathbb{S}}^{\$}-40{\mathbb{S}}^2+8\right)\end{array}$ (1.3)

With, ${\eta }_{i,j}=L_1^{-i}{L}_2^{-j}$ , $\mathbb{S}$ = SinI (I is the inclination of the orbit), μ is the gravitational parameter of the planet and zero order quantities defined as:

$\begin{array}{l}{A}_{11}=\frac{{\mu }^4{r}_e^2}{4},\\ A_{22}=\frac{{\mu }^6{r}_e^4{J}_4}{32J_2^2}.\end{array}$

where $l_i$ and $L_i$ are the Delaunay elements ( $L_1=\sqrt{\mu a}$ , $L_2=L_1\sqrt{\left(1-e^2\right)}$ , $L_3=L_{2}CosI$ ), $r_e$ is the equatorial radius of the Earth, and J₂, J₄ are the second and fourth geopotential zonal harmonic respectively.

The problem of designing invariant relative orbits for spacecraft flying formations is outlined as follows:

1) Compute the secular drift of the longitude of the ascending node and the sum of the argument of perigee and mean anomaly.

2) These secular drift rates are set equal between two neighboring orbits.

3) Having both orbits drift at equal angular rates on the average, they will not separate over time due to the influence of the perturbative effects of the asphericity of the Earth up to the desired order of magnitude (or the accuracy) to the equations of motion.

Using the canonical equations of motion,

${\stackrel{\dot{}}{l}}_i=\frac{\partial H^{**}}{\partial L_i},{\stackrel{\dot{}}{L}}_i=-\frac{\partial H^{**}}{\partial l},i=1,2,3$ (2)

Since the argument of mean latitude $\theta$ is the sum of the mean anomaly and the argument of perigee (l₁ + l₂). Evaluating the derivatives yields the sum of the argument of perigee and the mean anomaly rate of changes. Follows, the rate of change of mean latitude, $\stackrel{\dot{}}{\theta }$ , and the secular drift rates of the longitude of the ascending node, ${\stackrel{\dot{}}{l}}_3$ can be calculated, i.e. $\stackrel{\dot{}}{\theta }={\stackrel{\dot{}}{l}}_1+{\stackrel{\dot{}}{l}}_2$ .

So, using Equation (1) in Equation (2), the result can written in the form:

$\stackrel{\dot{}}{\theta }={\displaystyle \underset{n=0}{\overset{2}{\sum }}\frac{J_2^n}{n!}{D}_n^{\theta }}$ and ${\stackrel{\dot{}}{l}}_3={\displaystyle \underset{n=1}{\overset{2}{\sum }}\frac{J_2^n}{n!}{D}_n^{l_3}}$ (3)

With, $D_n^{\theta }$ and $D_n^{l_3}$ are published by Abd El-salam et al. [5] and given by:

$D_0^{\theta }={\mu }^2{\bar{K}}_0$

$D_1^{l_3}=A_{11}{Z}_1$

$D_1^{\theta }=A_{11}{\displaystyle \underset{i=1}{\overset{2}{\sum }}{\bar{K}}_0}$

$D_2^{\theta }=\frac{3}{128}\frac{A_{11}^2}{{\mu }^2}{\displaystyle \underset{i=2}{\overset{5}{\sum }}{Z}_i}+\frac{3}{2}{A}_{22}{\displaystyle \underset{i=6}{\overset{7}{\sum }}{Z}_i}$

$D_2^{\theta }=\frac{3}{128}\frac{A_{11}^2}{{\mu }^2}{\displaystyle \underset{i=3}{\overset{8}{\sum }}{\bar{K}}_i}+\frac{3}{2}{A}_{22}{\displaystyle \underset{i=9}{\overset{12}{\sum }}{\bar{K}}_i}$

where ${\bar{K}}_i$ and $Z_i$ are function of the action variable and introduced in Appendix.

To prevent the satellites from drifting apart over time, the average secular growth needs to be equal. So, it would be desirable to match all three rates ${\stackrel{\dot{}}{l}}_1,{\stackrel{\dot{}}{l}}_2$ and ${\stackrel{\dot{}}{l}}_3$ between the satellites in each formation. So $\stackrel{\dot{}}{\theta }$ and ${\stackrel{\dot{}}{l}}_3$ of all satellites in the formation should be equal.

$\begin{array}{l}{\stackrel{\dot{}}{\theta }}_i={\stackrel{\dot{}}{l}}_{1,i}+{\stackrel{\dot{}}{l}}_{2,i}={\stackrel{\dot{}}{l}}_{1,j}+{\stackrel{\dot{}}{l}}_{2,j}={\stackrel{\dot{}}{\theta }}_j\forall i\ne j\\ {\stackrel{\dot{}}{l}}_{3,i}={\stackrel{\dot{}}{l}}_{3,j}\forall i\ne j\end{array}$ (4)

Denoting the reference means orbit elements with the subscript “0”. Using Taylor expansion for the drift rate ${\stackrel{\dot{}}{\theta }}_i$ and ${\stackrel{\dot{}}{l}}_{3,i}$ of a neighboring orbit “i” about the reference orbital elements, retaining the terms up the second-order derivatives, can be simplified as:

$\begin{array}{c}\delta {\stackrel{\dot{}}{\theta }}_i=\left[{\frac{\partial {\stackrel{\dot{}}{\theta }}_i}{\partial L_1}|}_{x=x_0}+{\eta }_{-1,1}\left({\frac{\partial {\stackrel{\dot{}}{\theta }}_i}{\partial L_2}|}_{x=x_0}+ℂ{\frac{\partial {\stackrel{\dot{}}{\theta }}_i}{\partial L_3}|}_{x=x_0}\right)\right]\delta L_1\\ +{\eta }_{-1,0}\left[{\frac{\partial {\stackrel{\dot{}}{\theta }}_i}{\partial L_2}|}_{x=x_0}+ℂ{\frac{\partial {\stackrel{\dot{}}{\theta }}_i}{\partial L_3}|}_{x=x_0}\right]\delta {\eta }_{-1,1}-\left[{\eta }_{-2,1}\mathbb{S}{\frac{\partial {\stackrel{\dot{}}{\theta }}_i}{\partial L_3}|}_{x=x_0}\right]\delta I\\ +\frac{1}{2}[{\frac{{\partial }^2{\stackrel{\dot{}}{\theta }}_i}{\partial L_1^2}|}_{x=x_0}+{\eta }_{-2,2}\left({\frac{{\partial }^2{\stackrel{\dot{}}{\theta }}_i}{\partial L_2^2}|}_{x=x_0}+{ℂ}^2{\frac{{\partial }^2{\stackrel{\dot{}}{\theta }}_i}{\partial L_3^2}|}_{x=x_0}+2ℂ{\frac{{\partial }^2{\stackrel{\dot{}}{\theta }}_i}{\partial L_2\partial L_3}|}_{x=x_0}\right)\\ +2{\eta }_{-1,1}\left({\frac{{\partial }^2{\stackrel{\dot{}}{\theta }}_i}{\partial L_1\partial L_2}|}_{x=x_0}+ℂ{\frac{{\partial }^2{\stackrel{\dot{}}{\theta }}_i}{\partial L_1\partial L_3}|}_{x=x_0}\right)]{\left(\delta L_1\right)}^2+\frac{1}{2}{\eta }_{-2,0}[{\frac{{\partial }^2{\stackrel{\dot{}}{\theta }}_i}{\partial L_2^2}|}_{x=x_0}\\ +2ℂ{\frac{{\partial }^2{\stackrel{\dot{}}{\theta }}_i}{\partial L_2\partial L_3}|}_{x=x_0}+{ℂ}^2{\frac{{\partial }^2{\stackrel{\dot{}}{\theta }}_i}{\partial L_3^2}|}_{x=x_0}]{\left(\delta {\eta }_{-1,1}\right)}^2+\left[\frac{1}{2}{\eta }_{-4,2}{\mathbb{S}}^2{\frac{{\partial }^2{\stackrel{\dot{}}{\theta }}_i}{\partial L_3^2}|}_{x=x_0}\right]{\left(\delta I\right)}^2\end{array}$

$\begin{array}{l}+[{\eta }_{-2,1}({\frac{{\partial }^2{\dot{\theta }}_i}{\partial L_2^2}|}_{x=x_0}+{ℂ}^2{\frac{{\partial }^2{\dot{\theta }}_i}{\partial L_3^2}|}_{x=x_0}+2ℂ{\frac{{\partial }^2{\dot{\theta }}_i}{\partial L_2\partial L_3}|}_{x=x_0})\\ +{\eta }_{-1,0}\left({\frac{{\partial }^2{\dot{\theta }}_i}{\partial L_1\partial L_2}|}_{x=x_0}+ℂ{\frac{{\partial }^2{\dot{\theta }}_i}{\partial L_1\partial L_3}|}_{x=x_0}\right)]\left(\delta L_1\right)\left(\delta {\eta }_{-1,1}\right)\\ -{\eta }_{-3,1}\mathbb{S}\left[ℂ{\frac{{\partial }^2{\dot{\theta }}_i}{\partial L_3^2}|}_{x=x_0}-{\frac{{\partial }^2{\dot{\theta }}_i}{\partial L_2\partial L_3}|}_{x=x_0}\right]\left(\delta I\right)\left(\delta {\eta }_{-1,1}\right)\\ -\mathbb{S}[{\eta }_{-3,2}\left(ℂ{\frac{{\partial }^2{\dot{\theta }}_i}{\partial L_3^2}|}_{x=x_0}+{\frac{{\partial }^2{\dot{\theta }}_i}{\partial L_2\partial L_3}|}_{x=x_0}\right)+{\eta }_{-2,1}{\frac{{\partial }^2{\dot{\theta }}_i}{\partial L_1\partial L_3}|}_{x=x_0}]\left(\delta L_1\right)\left(\delta I\right)\end{array}$ (5)

With $ℂ$ = CosI, similarly:

$\begin{array}{c}\delta {\stackrel{\dot{}}{l}}_{3,i}=\left[{\frac{\partial {\stackrel{\dot{}}{l}}_{3,i}}{\partial L_1}|}_{x=x_0}+{\eta }_{-1,1}\left({\frac{\partial {\stackrel{\dot{}}{l}}_{3,i}}{\partial L_2}|}_{x=x_0}+ℂ{\frac{\partial {\stackrel{\dot{}}{l}}_{3,i}}{\partial L_3}|}_{x=x_0}\right)\right]\delta L_1\\ +{\eta }_{-1,0}\left[{\frac{\partial {\stackrel{\dot{}}{l}}_{3,i}}{\partial L_2}|}_{x=x_0}+ℂ{\frac{\partial {\stackrel{\dot{}}{l}}_{3,i}}{\partial L_3}|}_{x=x_0}\right]\delta {\eta }_{-1,1}-\left[{\eta }_{-2,1}\mathbb{S}{\frac{\partial {\stackrel{\dot{}}{l}}_{3,i}}{\partial L_3}|}_{x=x_0}\right]\delta I\\ +[\frac{1}{2}{\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_1^2}|}_{x=x_0}+\frac{1}{2}{\eta }_{-2,2}\left({\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_2^2}|}_{x=x_0}+{ℂ}^2{\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_3^2}|}_{x=x_0}2ℂ{\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_2\partial L_3}|}_{x=x_0}\right)\end{array}$

$\begin{array}{c}+{\eta }_{-1,1}\left({\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_1\partial L_2}|}_{x=x_0}+ℂ{\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_1\partial L_3}|}_{x=x_0}\right)]{\left(\delta L_1\right)}^2+\frac{1}{2}{\eta }_{-2,0}[{\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_2^2}|}_{x=x_0}\\ +2ℂ{\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_2\partial L_3}|}_{x=x_0}+{ℂ}^2{\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_3^2}|}_{x=x_0}]{\left(\delta {\eta }_{-1,1}\right)}^2+\left[\frac{1}{2}{\eta }_{-4,2}{\mathbb{S}}^2{\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_3^2}|}_{x=x_0}\right]{\left(\delta I\right)}^2\\ +[{\eta }_{-2,1}({\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_2^2}|}_{x=x_0}+{ℂ}^2{\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_3^2}|}_{x=x_0}+2{\eta }_{-2,1}ℂ{\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_2\partial L_3}|}_{x=x_0})\end{array}$

$\begin{array}{c}+{\eta }_{-1,0}\left({\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_1\partial L_2}|}_{x=x_0}+ℂ{\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_1\partial L_3}|}_{x=x_0}\right)]\left(\delta L_1\right)\left(\delta {\eta }_{-1,1}\right)\\ -{\eta }_{-3,1}\mathbb{S}\left[ℂ{\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_3^2}|}_{x=x_0}+{\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_2\partial L_3}|}_{x=x_0}\right]\left(\delta I\right)\left(\delta {\eta }_{-1,1}\right)\\ -\mathbb{S}[{\eta }_{-3,2}\left(ℂ{\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_3^2}|}_{x=x_0}+{\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_2\partial L_3}|}_{x=x_0}\right)+{\eta }_{-2,1}{\frac{{\partial }^2{\stackrel{\dot{}}{l}}_{3,i}}{\partial L_1\partial L_3}|}_{x=x_0}]\left(\delta L_1\right)\left(\delta I\right)\end{array}$ (6)

where $\delta {\stackrel{\dot{}}{\theta }}_i={\stackrel{\dot{}}{\theta }}_i-{\stackrel{\dot{}}{\theta }}_0$ is the difference between the drift rates of the argument of mean latitude of the reference orbit and one of the neighboring orbits, and ${{\stackrel{\dot{}}{\theta }}_i|}_{X=X_0}={\stackrel{\dot{}}{\theta }}_0$ .

And $\delta {\stackrel{\dot{}}{l}}_{3,i}={\stackrel{\dot{}}{l}}_{3,i}-{\stackrel{\dot{}}{l}}_{3,0}$ is the difference between the drift rates of the ascending node of the reference orbit and one of the neighboring orbits, and ${{\stackrel{\dot{}}{l}}_{3,i}|}_{X=X_0}={\stackrel{\dot{}}{l}}_{3,0}$ .

Now, the conditions satisfying the invariance property for the relative orbits are:

$\delta {\stackrel{\dot{}}{\theta }}_i={\stackrel{\dot{}}{\theta }}_i-{\stackrel{\dot{}}{\theta }}_0=0$ (7)

$\delta {\stackrel{\dot{}}{l}}_{3,i}={\stackrel{\dot{}}{l}}_{3,i}-{\stackrel{\dot{}}{l}}_{3,0}=0$ (8)

Substituting the included derivatives of the last two equations into Equations (5) and (6) and after the needed mathematical manipulations we will get:

$\begin{array}{l}{A}_{L_1{L}_1}^{\theta }{\left(\delta L_1\right)}^2+\left[A_I^{\theta }\delta I+A_{II}^{\theta }{\left(\delta I\right)}^2\right]+A_{\eta \eta }^{\theta }{\left(\delta {\eta }_{1,-1}\right)}^2+A_{L_1\eta }^{\theta }\left(\delta L_1\right)\left(\delta {\eta }_{1,-1}\right)\\ +\left[A_{\eta }^{\theta }+A_{I\eta }^{\theta }\left(\delta I\right)\right]\left(\delta {\eta }_{1,-1}\right)+\left[A_{L_1}^{\theta }+A_{L_{1}I}^{\theta }\left(\delta I\right)\right]\left(\delta L_1\right)=0\end{array}$ (7')

$\begin{array}{l}{A}_{L_1{L}_1}^{l_3}{\left(\delta L_1\right)}^2+\left[A_I^{l_3}\delta I+A_{II}^{l_3}{\left(\delta I\right)}^2\right]+A_{\eta \eta }^{l_3}{\left(\delta {\eta }_{1,-1}\right)}^2+A_{L_1\eta }^{l_3}\left(\delta L_1\right)\left(\delta {\eta }_{1,-1}\right)\\ +\left[A_{\eta }^{l_3}+A_{I\eta }^{l_3}\left(\delta I\right)\right]\left(\delta {\eta }_{1.-1}\right)+\left[A_{L_1}^{l_3}+A_{L_{1}I}^{l_3}\left(\delta I\right)\right]\left(\delta L_1\right)=0\end{array}$ (8')

Multiplying Equation (7') by $A_{L_1{L}_1}^{l_3}$ and Equation (8') by $A_{L_1{L}_1}^{\theta }$ and subtracting we will get:

$a_1\delta L_1\delta {\eta }_{1,-1}+a_2{\left(\delta {\eta }_{1,-1}\right)}^2+a_3\delta L_1+a_4\delta {\eta }_{1,-1}+a_5=0$ (9)

with the coefficients $a_i$ ’s are:

$a_1={\mathcal{A}}_{L_1\eta }^{\theta }{\mathcal{A}}_{L_1{L}_1}^{l_3}-{\mathcal{A}}_{L_1\eta }^{l_3}{\mathcal{A}}_{L_1{L}_1}^{\theta }$

$a_2={\mathcal{A}}_{\eta \eta }^{\theta }{\mathcal{A}}_{L_1{L}_1}^{l_3}-{\mathcal{A}}_{\eta \eta }^{l_3}{\mathcal{A}}_{L_1{L}_1}^{\theta }$

$a_3=\left[{\mathcal{A}}_{L_1}^{\theta }+{\mathcal{A}}_{I{L}_1}^{\theta }\left(\delta I\right)\right]{\mathcal{A}}_{L_1{L}_1}^{l_3}-{\mathcal{A}}_{L_1{L}_1}^{\theta }\left[{\mathcal{A}}_{L_1}^{l_3}+{\mathcal{A}}_{I{L}_1}^{l_3}\left(\delta I\right)\right]$

$a_4=\left[{\mathcal{A}}_{\eta }^{\theta }+{\mathcal{A}}_{I\eta }^{\theta }\left(\delta I\right)\right]{\mathcal{A}}_{L_1{L}_1}^{l_3}-{\mathcal{A}}_{L_1{L}_1}^{\theta }\left[{\mathcal{A}}_{\eta }^{l_3}+{\mathcal{A}}_{I\eta }^{l_3}\left(\delta I\right)\right]$

$a_5=\left[{\mathcal{A}}_I^{\theta }\left(\delta I\right)+{\mathcal{A}}_{II}^{\theta }{\left(\delta I\right)}^2\right]{\mathcal{A}}_{L_1{L}_1}^{l_3}-{\mathcal{A}}_{L_1{L}_1}^{\theta }\left[{\mathcal{A}}_I^{l_3}\left(\delta I\right)+{\mathcal{A}}_{II}^{l_3}{\left(\delta I\right)}^2\right]$

and the derivatives of ${\mathcal{A}}^{\theta }$ ’s and ${\mathcal{A}}^{l_3}$ ’s are in the appendix. Solving Equation (9), for $\delta L_1$ , we get:

$\delta L_1=-\frac{a_2{\left(\delta {\eta }_{1,-1}\right)}^2+a_4\delta {\eta }_{1,-1}+a_5}{a_1\delta {\eta }_{1,-1}+a_3}$

Substituting the last results in Equation (7') we get the quartic equation:

$b_1{\left(\delta {\eta }_{1,-1}\right)}^4+b_2{\left(\delta {\eta }_{1,-1}\right)}^3+b_3{\left(\delta {\eta }_{1,-1}\right)}^2+b_4\left(\delta {\eta }_{1,-1}\right)+b_5=0$ ,

where the coefficients $b_i$ ’s are functions of $L_i$ ’s and given by:

$b_1=a_2^2{\mathcal{A}}_{L_1{L}_1}^{\theta }-a_2{a}_1{\mathcal{A}}_{L_1\eta }^{\theta }+a_1^2{\mathcal{A}}_{\eta \eta }^{\theta }$

$\begin{array}{l}{b}_2=2a_2{a}_4{\mathcal{A}}_{L_1{L}_1}^{\theta }-\left(a_2{a}_3+a_1{a}_4\right){\mathcal{A}}_{L_1\eta }^{\theta }-a_2{a}_1\left[{\mathcal{A}}_{L_1}^{\theta }+{\mathcal{A}}_{I{L}_1}^{\theta }\left(\delta I\right)\right]\\ +a_1^2\left[{\mathcal{A}}_{\eta }^{\theta }+{\mathcal{A}}_{I\eta }^{\theta }\left(\delta I\right)\right]+2a_1{a}_3{\mathcal{A}}_{\eta \eta }^{\theta }\end{array}$

$\begin{array}{c}{b}_3=\left(a_4^2+2a_2{a}_5\right){\mathcal{A}}_{L_1{L}_1}^{\theta }-\left(a_4{a}_3+a_5{a}_1\right){\mathcal{A}}_{L_1\eta }^{\theta }+a_3^2{\mathcal{A}}_{\eta \eta }^{\theta }\\ -\left(a_2{a}_3+a_1{a}_4\right)\left[{\mathcal{A}}_{L_1}^{\theta }+{\mathcal{A}}_{I{L}_1}^{\theta }\left(\delta I\right)\right]+2a_1{a}_3\left[{\mathcal{A}}_{\eta }^{\theta }+{\mathcal{A}}_{I\eta }^{\theta }\left(\delta I\right)\right]\\ +a_1^2\left[{\mathcal{A}}_I^{\theta }\left(\delta I\right)+{\mathcal{A}}_{II}^{\theta }{\left(\delta I\right)}^2\right]\end{array}$

$\begin{array}{c}{b}_4=2a_4{a}_5{\mathcal{A}}_{L_1{L}_1}^{\theta }-a_5{a}_3{\mathcal{A}}_{L_1\eta }^{\theta }-\left(a_4{a}_3+a_5{a}_1\right){\mathcal{A}}_{L_1}^{\theta }+{\mathcal{A}}_{I{L}_1}^{\theta }\left(\delta I\right)\\ +a_3^2\left[{\mathcal{A}}_{\eta }^{\theta }+{\mathcal{A}}_{I\eta }^{\theta }\left(\delta I\right)\right]+2a_3{a}_1\left[{\mathcal{A}}_I^{\theta }\left(\delta I\right)+{\mathcal{A}}_{II}^{\theta }{\left(\delta I\right)}^2\right]\end{array}$

$b_5=a_5^2{\mathcal{A}}_{L_1{L}_1}^{\theta }-a_5{a}_3\left[{\mathcal{A}}_{L_1}^{\theta }+{\mathcal{A}}_{I{L}_1}^{\theta }\left(\delta I\right)\right]+a_3^2\left[{\mathcal{A}}_I^{\theta }\left(\delta I\right)+{\mathcal{A}}_{II}^{\theta }{\left(\delta I\right)}^2\right]$

Solving the resulting equation we will get four roots of $\delta {\eta }_{1,-1}$ :

$\begin{array}{l}{\left(\delta {\eta }_{1,-1}\right)}_{1,2}=A_1-A_2\mp A_3,\\ {\left(\delta {\eta }_{1,-1}\right)}_{3,4}=A_1+A_2\mp A_4.\end{array}$ (10)

and four roots of $\delta L_1$ :

$\begin{array}{l}{\left(\delta L_1\right)}_1=-\frac{a_2\left(A_5-2A_1{A}_2+2A_2{A}_3-2A_1{A}_3\right)+a_4\left(A_1-A_2-A_3\right)+a_5}{a_1\left(A_1-A_2-A_3\right)+a_3},\\ {\left(\delta L_1\right)}_2=-\frac{a_2\left(A_5-2A_1{A}_2-2A_2{A}_3+2A_1{A}_3\right)+a_4\left(A_1-A_2-A_3\right)+a_5}{a_1\left(A_1-A_2-A_3\right)+a_3},\\ {\left(\delta L_1\right)}_3=-\frac{a_2\left(A_6+2A_1{A}_2-2A_2{A}_3-2A_1{A}_3\right)+a_4\left(A_1+A_2-A_3\right)+a_5}{a_1\left(A_1+A_2-A_3\right)+a_3},\\ {\left(\delta L_1\right)}_4=-\frac{a_2\left(A_6+2A_1{A}_2+2A_2{A}_3+2A_1{A}_3\right)+a_4\left(A_1+A_2-A_3\right)+a_5}{a_1\left(A_1+A_2-A_3\right)+a_3},\end{array}$ (10')

where,

$A_1=-\frac{b_2}{4b_1}$ , $A_2=\frac{1}{2}\sqrt{B_1+B_2+B_3}$ , $A_3=\frac{1}{2}\sqrt{2B_1-B_2-B_3-B_4}$ ,

$A_4=\frac{1}{2}\sqrt{2B_1-B_2-B_3+B_4}$ , $A_5=A_1^2+A_2^2+A_3^2$ , $A_6=A_1^2+A_2^2+A_4^2$ ,

$B_1=\frac{b_2^2}{4b_1^2}-\frac{2b_3}{3b_1}$ , $B_2=\frac{2^{1/3}{C}_1}{3b_1{C}_4}$ , $B_3=\frac{C_4}{3\times 2^{1/3}{b}_1}$ , $B_4=\frac{C_2}{8A_2}$ ,

$C_1=b_3^2-3b_2{b}_4+12b_1{b}_5$ , $C_2=-\frac{b_2^3}{b_1^3}+\frac{4b_2{b}_3}{b_1^2}-\frac{8b_4}{b_1}$ ,

$C_3=2b_3^3-9b_2{b}_3{b}_4+27b_1{b}_4^2+27b_2^2{b}_5-72b_1{b}_3{b}_5$ ,

$C_4={\left(C_3+\sqrt{-4{\left(C_1\right)}^3+C_3^2}\right)}^{\frac{1}{3}}$ .

3. Modeling of Invariant Relative Orbit Conditions for near Equatorial-Circular Case

When we choose the leader orbit to be circular equatorial then, we will obtain four solutions for $\delta {\eta }_{1,-1}$ (which we redefined as $\delta \eta$ for simplicity), as in Equation (10), and four solutions for $\delta L_1$ as in Equation (10'). Here, we will present the plots of these solutions, which we obtained in the last section (Equations (10) and (10')) to compare the effect of J₄ on the conditions of invariance for the formation. Before we introduce the graphs, it is important to mention that in all figures, we provided a set of curves in each condition for the inclination of each member of the formation with respect to the leader orbit $\left(-2^{\circ }\le \delta I\le 2^{\circ }\right)$ as example. By means that the follower satellites’ orbits, in the formation, will be inclined by a range of 2˚ with respect to the leader satellite orbit, of course we can extend this range. Also, we will introduce a comparison between the effects of $J_4$ in the formation.

The variation of the formation relative to the leader orbit in $\delta \eta$ is related to the variation in eccentricity for the followers through:

$\delta \eta =-\frac{e}{\sqrt{1-e^2}}\delta e=MAG\left(e\right)\delta e$ ,

i.e. the variation in the follower’s orbit in ${\eta }_{1,-1}$ ( $\delta \eta$ ) is scaled by $MAG\left(e\right)$ in their variations in eccentricities $\delta e$ . It is important to note that the scale function $MAG\left(e\right)$ is always negative for circular and elliptical orbits. In our case the eccentricity of the leader orbit, $e_0$ , is zero, then $\delta e=e-e_0$ is always positive.

While the variation in the $\delta L_1$ , for the followers, is scaled by $mag\left(a\right)$ for the variation of the semi-major axis through the relation $\delta L_1=\frac{1}{2\sqrt{\mu a}}\delta a=mag\left(a\right)\delta a$ . It is important to note that the scale function $mag\left(a\right)$ is always positive, but $\delta a$ maybe positive, negative or zero.

3.1. The First Solution (Circular Formation)

The first solution of Equations (10) for $\delta e$ and (10') for $\delta a$ for the orbits of the followers with respect to the leader orbit, for inclination range [−2˚, +2˚] in case of J₂ and the net effects of J₂ and J₄ can represent in the following curves.

Figure 1(a) and Figure 1(b) show the first choice of the invariant relative conditions, in the equatorial-circular case. In this choice, $\delta \eta$ and $\delta L_1$ get Zero for all $\delta a,\delta e$ and $\delta I$ values. That is mean, the formation for the followers orbit will be in the same orbital plane of the leader satellite (In Plane Circular Formation) with their eccentricities and semi-major axis was scaled by $MAG\left(e\right)$ and $mag\left(a\right)$ of the leader satellite. The scale function $MAG\left(e\right)$ , for this solution, is equal zero whatever the choosing the value of $\delta e$ . Also the scale function $mag\left(a\right)$ never equal zero, then $\delta a$ must equal zero. That can me conclude that the follower satellite’s must be in the same orbit of the leader one.

In this case, the effect of J₄ has no significant variation in the formation.

3.2. The Second Solution

Figure 2(a) and Figure 2(b) show the second choice of the invariant relative formation, in the equatorial-circular case. In this solution, the choice of eccentricites and semi-major axis for the followers orbit is not affect the formation in case of J₂ effect whatever choosing the inclination for folowers orbits. But the J₄ effects change the choosing of eccentricities slightly while the choosing the semi-major axis still unaffected whatever chosing the inclinations for the followers orbits. As we see in the vertical axis. The semi-magor axis and eccentricites of the followers orbits must be greater than those for the leader one.

The effects of J₄ are changing slightly while the choosing of eccentricies is not for choosing the semi-major axis.

3.3. The Third Solution

Figure 3(a) and Figure 3(b) show the third choice of the invariant relative formation, in the equatorial-circular case. In this solution, the value of the function $MAG\left(e\right)\delta e<-1.3$ which gives the limits for the eccentricities through the inqualety $e^4+1.69e^2-1.69\prec 0$ .

The solution for this inqualety has only one positive value 0.839936. With respect to choosing the semi-major axis, it is not depend on choosing the inclination, by mean that the formation will be in plane of the leader one.

In the case of the formation, under the effects of J₂ and J₄, the function $MAG\left(e\right)\delta e$ is modified the limit of the eccentricities to be 0.774904.

3.4. The Fourth Solution

Figure 4(a) and Figure 4(b) show the fourth choice of the invariant relative formation, in the equatorial-circular case. In this choice, the formation for defferent inclination is distributed about the leader orbit (with $\delta I=0$ ) and the function $MAG\left(e\right)\delta e$ is increasing for positive $\delta I$ while it is decreasing for negative $\delta I$ . For choosing the semi-major axis, it increases by increasing the semi-major axis of the leader orbit whatever choosing $\delta I$ .

In this formation choics, the effect of J₄ is significant for choosing the eccentricities for the followers as it clear from the first of Figure 4(b). In the first of

Figure 4(b), the function $MAG\left(e\right)\delta e$ is positive and increasing while $MAG\left(e\right)$ is always negative. And $\delta e$ , in our case, must be positive or zero. For that this solution is not mathematically accepted.

Also, the second of Figure 4(b) shows that this choice is not continuous for orbits with semi-major axis greater than 2.5 earth radii. And the function $mag\left(a\right)\delta a$ is decreasing by increasing the semi-major axis.

4. Conclusions

The problem formulated using the oblate Earth model, truncating its potential series at J₄ to the equations of motion, and then the canonical equations of motion and the Hamiltonian formed. In order to keep the relative motion invariable, eight-second order conditions between the differences in the semi-major axis a and the inclination I are obtained. These conditions guarantee that the drift rates of neighboring orbits are equal on the average. The resulting orbits require less control and maintenance fuel. Then we studied the curves of these conditions in equatorial-circular case. The plots of these cases are presented as relations between $\delta \eta$ or $\delta L_1$ and the semi-major axis of the leader satellite orbit, at different $\delta I$ .

In the first choice, we can conclude that the follower satellite’s must be in the same orbit of the leader one (on orbit formation).

In the second choice, whatever the choosing the inclination for the followers it is not affect the choosing the semi-major axis (In plane with different eccentrisities).

The third choice, under the effects of J₂ and J₄, the value of the function $MAG\left(e\right)\delta e<-1.3$ which gives the limits for the eccentricities through the inequality $e^4+1.69e^2-1.69\prec 0$ . The solution for this inqualety has only one positive value 0.839936. With respect to choosing the semi-major axis, it does not depend on choosing the inclination, by mean that the formation will be in plane of the leader one.

The fourth choice, the function $MAG\left(e\right)\delta e$ is positive and increasing while $MAG\left(e\right)$ is always negative. And $\delta e$ , in our case, must be positive or zero. For that, the eccentricities of the followers must be negative and that is not mathematically accepted.

Also, the second of Figures 4(b) shows that this choice is not continuous for orbits with semi-major axis greater than 2.5 earth radii. And the function $mag\left(a\right)\delta a$ is decreasing by increasing the semi-major axis.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Altalhi, S.A. and El Saftawy, M.I. (2019) Analytical Solution for Formation Flying Problem near Equatorial-Circular Reference Orbit. International Journal of Astronomy and Astrophysics, 9, 217-230. https://doi.org/10.4236/ijaa.2019.93016

References

Appendix

${\mathcal{A}}_{L_1}^{\theta }=\underset{n=0}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{J}}_n^{\theta }\right)$ , ${\mathcal{A}}_{\eta }^{\theta }=\underset{n=1}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{K}}_{\theta }^n\right)$ , ${\mathcal{A}}_I^{\theta }=\underset{n=1}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{L}}_n^{\theta }\right)$ ,

${\mathcal{A}}_{L_1{L}_1}^{\theta }=\underset{n=0}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{M}}_n^{\theta }\right)$ , ${\mathcal{A}}_{\eta \eta }^{\theta }=\underset{n=1}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{N}}_n^{\theta }\right)$ , ${\mathcal{A}}_{II}^{\theta }=\underset{n=1}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{P}}_n^{\theta }\right)$ ,

${\mathcal{A}}_{L_1\eta }^{\theta }=\underset{n=1}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{Q}}_n^{\theta }\right)$ , ${\mathcal{A}}_{I\eta }^{\theta }=\underset{n=1}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{R}}_n^{\theta }\right)$ , ${\mathcal{A}}_{I{L}_1}^{\theta }=\underset{n=1}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{S}}_n^{\theta }\right)$ .

${\mathcal{A}}_{L_1}^{l_3}=\underset{n=1}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{J}}_n^{l_3}\right)$ , ${\mathcal{A}}_{\eta }^{l_3}=\underset{n=1}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{K}}_n^{l_3}\right)$ , ${\mathcal{A}}_I^{l_3}=\underset{n=1}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{L}}_n^{l_3}\right)$ ,

${\mathcal{A}}_{L_1{L}_1}^{l_3}=\underset{n=1}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{M}}_n^{l_3}\right)$ , ${\mathcal{A}}_{\eta \eta }^{l_3}=\underset{n=1}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{N}}_n^{\theta }\right)$ , ${\mathcal{A}}_{II}^{l_3}=\underset{n=1}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{P}}_n^{l_3}\right)$ ,

${\mathcal{A}}_{L_1\eta }^{l_3}=\underset{n=1}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{Q}}_n^{l_3}\right)$ , ${\mathcal{A}}_{I\eta }^{l_3}=\underset{n=1}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{R}}_n^{l_3}\right)$ , ${\mathcal{A}}_{I{L}_1}^{l_3}=\underset{n=1}{\overset{2}{{\displaystyle \sum }}}\left(\frac{J_2^n}{n!}{\mathcal{S}}_n^{l_3}\right)$ .

With,

${\mathcal{J}}_0^{\theta }=\frac{\partial D_0^{\theta }}{\partial L_1}$ ,

${\mathcal{J}}_n^{\theta }=\frac{\partial D_n^{\theta }}{\partial L_1}+{\eta }_{1,-1}\frac{\partial D_n^{\theta }}{\partial L_2}+{\eta }_{1,-1}ℂ\frac{\partial D_n^{\theta }}{\partial L_3}$ ,

${\mathcal{K}}_n^{\theta }=L_1\frac{\partial D_n^{\theta }}{\partial L_2}+L_1ℂ\frac{\partial D_n^{\theta }}{\partial L_3}$ ,

${\mathcal{L}}_n^{\theta }=-L_2\mathbb{S}\frac{\partial D_n^{\theta }}{\partial L_3}$ ,

${\mathcal{M}}_0^{\theta }=\frac{1}{2}\frac{{\partial }^2{D}_0^{\theta }}{\partial L_1\partial L_1}$ ,

$\begin{array}{l}{\mathcal{M}}_n^{\theta }=\frac{1}{2}\frac{{\partial }^2{D}_n^{\theta }}{\partial L_1\partial L_1}+\frac{1}{2}{\eta }_{2,-2}\left[\frac{{\partial }^2{D}_n^{\theta }}{\partial L_2\partial L_2}+{ℂ}^2\frac{{\partial }^2{D}_n^{\theta }}{\partial L_3\partial L_3}+2ℂ\frac{{\partial }^2{D}_n^{\theta }}{\partial L_2\partial L_3}\right]\\ +{\eta }_{1,-1}\left[ℂ\frac{{\partial }^2{D}_n^{\theta }}{\partial L_1\partial L_3}+\frac{{\partial }^2{D}_n^{\theta }}{\partial L_1\partial L_2}\right]\end{array}$ ,

${\mathcal{N}}_n^{\theta }=\frac{1}{2}{L}_1^2\left[\frac{{\partial }^2{D}_n^{\theta }}{\partial L_2\partial L_2}+{ℂ}^2\frac{{\partial }^2{D}_n^{\theta }}{\partial L_3\partial L_3}+2ℂ\frac{{\partial }^2{D}_n^{\theta }}{\partial L_2\partial L_3}\right]$ ,

${\mathcal{P}}_n^{\theta }=\frac{1}{2}{\eta }_{0,-2}{\mathbb{S}}^2\frac{{\partial }^2{D}_n^{\theta }}{\partial L_3\partial L_3}$ ,

${\mathcal{Q}}_n^{\theta }=L_2\left[\frac{{\partial }^2{D}_n^{\theta }}{\partial L_2\partial L_2}+{ℂ}^2\frac{{\partial }^2{D}_n^{\theta }}{\partial L_3\partial L_3}+2ℂ\frac{{\partial }^2{D}_n^{\theta }}{\partial L_2\partial L_3}\right]+L_1\left[\frac{{\partial }^2{D}_n^{\theta }}{\partial L_1\partial L_2}+ℂ\frac{{\partial }^2{D}_n^{\theta }}{\partial L_1\partial L_3}\right]$ ,

${\mathcal{R}}_n^{\theta }=-{\eta }_{-1,-1}\left[ℂ\frac{{\partial }^2{D}_n^{\theta }}{\partial L_3\partial L_3}+\mathbb{S}\frac{{\partial }^2{D}_n^{\theta }}{\partial L_2\partial L_3}\right]$ ,

${\mathcal{S}}_n^{\theta }=-{\eta }_{1,-2}\mathbb{S}\left[ℂ\frac{{\partial }^2{D}_n^{\theta }}{\partial L_3\partial L_3}+\frac{{\partial }^2{D}_n^{\theta }}{\partial L_2\partial L_3}+{\eta }_{-1,1}\frac{{\partial }^2{D}_n^{\theta }}{\partial L_1\partial L_3}\right]$ ,

${\mathcal{J}}_n^{l_3}=\frac{\partial D_n^{l_3}}{\partial L_1}+{\eta }_{1,-1}\left[\frac{\partial D_n^{l_3}}{\partial L_2}+ℂ\frac{\partial D_n^{l_3}}{\partial L_3}\right]$ ,

${\mathcal{K}}_n^{l_3}=L_1\left[\frac{\partial D_n^{l_3}}{\partial L_2}+ℂ\frac{\partial D_n^{l_3}}{\partial L_3}\right]$ ,