Dynamics of Satellite Formation Utilizing the Perturbed Restricted Three-Body Problem

Nabawia S. Khalifa

National Research Institute of Astronomy and Geophysics (NRIAG), Helwan, Egypt.

DOI: 10.4236/jamp.2022.101016                

Abstract

The dynamics of satellites formation is of great interest for the space mission. This work discusses a more efficient model of the relative motion dynamics of satellites formation. The model is based on employing the concepts restricted three-body problem (R3BP) and for more accuracy, it considers the effects of both oblateness and radiation pressure on deputy relative motion w.r.t the chief satellite. A model of deputy relative motion w.r.t the chief satellite is derived in the local-vertical local-horizontal system and simplified assuming the concept of the circular restricted three-body problem (CR3BP). The deputy equations of motion were rewritten in the form of recurrence relations and solved numerically using the Lie series approach. Assuming that the formation is revolving around the Moon in the Earth-Moon system, the effects of both oblateness and radiation pressure on the deputy satellite orbit were assessed through a particular example of satellites formation. A comparison between the perturbed and unperturbed R3BP shows a significant difference in the deputy relative position that has to be considered for the formation dynamics.

Keywords

Relative Motion, Restricted Three-Body Problem, Lie Series Integration

1. Introduction

Parallel to the start of the space programs, the study of spacecraft relative dynamics became one of the most important aspects in designing and analyzing space missions. Several authors were interested in this study to analyze rendezvous and docking of two spacecrafts in addition to maintenance of spacecraft formation. The most famous models of relative motion, Clohessy-Wiltshire and Tschauner-Hempel, have been used to analyze relative guidance, navigation and control systems. However, these two models assumed that the spacecraft relative distance is significantly small compared with its position w.r.t the centre of mass of the primary. Moreover, they assume pure Kepler motion (e.g. two-body problem without including any perturbations) [1] [2]. Afterwards, different models have been constructed to overcome the limitations of these two models [3] - [8].

In 2019, Giovanni Franzini and Mario Innocenti studied the relative motion dynamics using the classical restricted three-body problem (unperturbed problem). They found that a three-body scenario is more suitable than a two-body scenario to describe the dynamics of satellite relative motion and more efficient for studying the relative guidance and navigation system [9]. Recently, many authors are interested in modelling the formation dynamics employing the restricted three-body problem and studied the relative motion around the libration points [10] [11] [12] [13] [14].

The current study aims to get a more accurate formulation of the relative motion employing the perturbed restricted three-body problem. The dynamical model assumes that the primaries are radiating and oblate spheroids considering only zonal harmonics. Based on constructing recurrence formulas and applying the Lie series approach, the model is solved numerically. Finally, the numerical application is performed assuming a circular three-dimensional problem.

2. Formulation of the Problem

2.1. The Perturbed Restricted Three-Body Problem

The restricted three-body problem is one of the most famous dynamical modellings of celestial systems. Extensive studies were performed using a variety of methods assuming that the primaries are spheres and only their gravitational attraction is considered [15] - [20]. For a more accurate model, additional perturbing forces such as oblateness and radiation pressure are considered. The gravitational potential “ ${\varphi }_k$ ” of a massive body is given by [21] [22] [23]:

${\varphi }_k=-\frac{G{m}_k}{r_k}\left[1-{\displaystyle {\sum }_{n=2}^{\infty }{j}_n^k}{\left(\frac{R_k}{r_k}\right)}^n{P}_n^k\left(\sin \delta \right)\right]$ (1)

G is the universal gravitational constant, $R_k$ is the mean radius of each body, $\delta$ is the latitude of the infinitesimal body, $r_k$ is the separation between every two bodies and $j_n^k$ is the dimensionless coefficient which represents the non-spherical components of the potential. The Legendre polynomials “ $P_n^k\left(\sin \delta \right)$ ” of degree n is given by:

$P_n^k\left(\sin \delta \right)=\frac{1}{2^{n}n!}\frac{{\text{d}}^n}{\text{d}{\delta }^n}{\left({\sin }^2\delta -1\right)}^n$

The gravitational force “ ${\bar{F}}_{gr}$ ” exerted by the body is:

${\bar{F}}_{gr}=-\nabla {\varphi }_k$ (2)

where $\nabla$ denotes the vector differential operator. Apart from the gravitational potential, if the body is radiating, then it exerts a radiation force “ ${\bar{F}}_{rad}$ ” in the opposite direction of its gravitational force. Consequently, the total force “ ${\bar{F}}_k$ ” is:

${\bar{F}}_k={\bar{F}}_{gr}-{\bar{F}}_{rad}=q_k{\bar{F}}_{gr}$ (3)

where $q_k=\left(1-F_{rad}/F_{gr}\right)$ is the radiation factor $\in \left(0,1\right)$. Substitute (1) and (2) into (3):

${\bar{F}}_k=-q_{k}G{m}_k{\displaystyle {\sum }_{n=2}^{\infty }\left[\frac{1}{r_k^3}-A_n^{\left(k\right)}{\left(\frac{1}{r_k}\right)}^{n+3}\right]{\bar{r}}_k}$ (4)

where

$A_n^{\left(k\right)}=\left(n+1\right)j_n^{\left(k\right)}{R}_k^n{P}_n^{\left(k\right)}\left(\sin \delta \right)$

Consider a system of three bodies that have masses $m_1$, $m_2$ and $m_3$ such that the masses of primaries $m_1>m_2$ and $m_3$ is the mass of the chief satellite in a formation. Let $\left(I:\hat{i},\hat{j},\hat{k}\right)$ is a sidereal (inertial) coordinate system with the origin lies in the centre of mass of the two primaries, the motion of each body is given by:

${\left[{\ddot{\stackrel{¯}{r}}}_k\right]}_I=-G{\displaystyle {\sum }_{i=1}^3{\displaystyle {\sum }_{n=2}^{\infty }{q}_i{m}_i\left[\frac{1}{{\Vert {\bar{r}}_{ik}\Vert }^3}-A_n^{\left(i\right)}{\left(\frac{1}{\Vert {\bar{r}}_{ik}\Vert }\right)}^{n+3}\right]{\bar{r}}_{ik}}}$

where ${\bar{r}}_{ik}={\bar{r}}_k-{\bar{r}}_i$, $i=1,2,3$ and $k=1,2,3$ and $i\ne k$. For the restricted case of three bodies, the primaries are not affected by the gravitational influence of the satellite. Assuming that the formation is revolving around the smaller primary, the equation of motion of the chief satellite is (see Figure 1):

$\begin{array}{c}{\left[{\ddot{\stackrel{¯}{r}}}_{23}\right]}_I={\left[{\ddot{\stackrel{¯}{r}}}_3\right]}_I-{\left[{\ddot{\stackrel{¯}{r}}}_2\right]}_I\\ =-q_{1}G{m}_1{\displaystyle {\sum }_{n=2}^{\infty }\left[\frac{1}{{\Vert {\bar{r}}_{13}\Vert }^3}-A_n^{\left(1\right)}{\left(\frac{1}{\Vert {\bar{r}}_{13}\Vert }\right)}^{n+3}\right]{\bar{r}}_{13}}\\ -q_{2}G{m}_2{\displaystyle {\sum }_{n=2}^{\infty }\left[\frac{1}{{\Vert {\bar{r}}_{23}\Vert }^3}-A_n^{\left(2\right)}{\left(\frac{1}{\Vert {\bar{r}}_{23}\Vert }\right)}^{n+3}\right]{\bar{r}}_{23}}\\ -q_{1}G{m}_1{\displaystyle {\sum }_{n=2}^{\infty }\left[\frac{1}{{\Vert {\bar{r}}_{12}\Vert }^3}-A_n^{\left(1\right)}{\left(\frac{1}{\Vert {\bar{r}}_{12}\Vert }\right)}^{n+3}\right]{\bar{r}}_{12}}\end{array}$

Figure 1. Restricted three-body problem in an inertial frame.

As clear in Figure 1, ${\bar{r}}_{23}={\bar{r}}_c$ is the relative position of the chief w.r.t. the smaller primary and ${\bar{r}}_{13}={\bar{r}}_{12}+{\bar{r}}_c$, then

$\begin{array}{c}{\left[{\ddot{\stackrel{¯}{r}}}_c\right]}_I=-q_{1}G{m}_1\left(\frac{\left({\bar{r}}_{12}+{\bar{r}}_c\right)}{{\Vert {\bar{r}}_{12}+{\bar{r}}_c\Vert }^3}+\frac{{\bar{r}}_{12}}{{\Vert {\bar{r}}_{12}\Vert }^3}-{\displaystyle {\sum }_{n=2}^{\infty }{A}_n^{\left(1\right)}}\left[\frac{\left({\bar{r}}_{12}+{\bar{r}}_c\right)}{{\Vert {\bar{r}}_{12}+{\bar{r}}_c\Vert }^{n+3}}+\frac{{\bar{r}}_{12}}{{\Vert {\bar{r}}_{12}\Vert }^{n+3}}\right]\right)\\ -q_{2}G{m}_2{\displaystyle {\sum }_{n=2}^{\infty }\left[\frac{1}{{\Vert {\bar{r}}_c\Vert }^3}-A_n^{\left(2\right)}{\left(\frac{1}{\Vert {\bar{r}}_c\Vert }\right)}^{n+3}\right]{\bar{r}}_c}\end{array}$ (5)

Equation (5) can be normalized by assuming that the gravitational constant $G=1$, the distance between the two primaries is unity and the sum of their masses is unity. Let the normalized mass of the small primary parameter be $\mu =\frac{m_2}{m_1+m_2}$, then that of the big primary is “ $1-\mu$ ”. Consequently, the chief equation of motion will be:

$\begin{array}{c}{\left[{\ddot{\stackrel{¯}{r}}}_c\right]}_I=-q_1\left(1-\mu \right)\left(\frac{\left({\bar{r}}_{12}+{\bar{r}}_c\right)}{{\Vert {\bar{r}}_{12}+{\bar{r}}_c\Vert }^3}+\frac{{\bar{r}}_{12}}{{\Vert {\bar{r}}_{12}\Vert }^3}-{\displaystyle {\sum }_{n=2}^{\infty }{A}_n^{\left(1\right)}}\left[\frac{\left({\bar{r}}_{12}+{\bar{r}}_c\right)}{{\Vert {\bar{r}}_{12}+{\bar{r}}_c\Vert }^{n+3}}+\frac{{\bar{r}}_{12}}{{\Vert {\bar{r}}_{12}\Vert }^{n+3}}\right]\right)\\ -q_2\mu {\displaystyle {\sum }_{n=2}^{\infty }\left[\frac{1}{{\Vert {\bar{r}}_c\Vert }^3}-A_n^{\left(2\right)}{\left(\frac{1}{\Vert {\bar{r}}_c\Vert }\right)}^{n+3}\right]{\bar{r}}_c}\end{array}$ (6)

Similarly, the motion of a deputy satellite of the formation will be:

$\begin{array}{c}{\left[{\ddot{\stackrel{¯}{r}}}_d\right]}_I=-q_1\left(1-\mu \right)\left(\frac{\left({\bar{r}}_{12}+{\bar{r}}_d\right)}{{\Vert {\bar{r}}_{12}+{\bar{r}}_d\Vert }^3}+\frac{{\bar{r}}_{12}}{{\Vert {\bar{r}}_{12}\Vert }^3}-{\displaystyle {\sum }_{n=2}^{\infty }{A}_n^{\left(1\right)}}\left[\frac{\left({\bar{r}}_{12}+{\bar{r}}_d\right)}{{\Vert {\bar{r}}_{12}+{\bar{r}}_d\Vert }^{n+3}}+\frac{{\bar{r}}_{12}}{{\Vert {\bar{r}}_{12}\Vert }^{n+3}}\right]\right)\\ -q_2\mu {\displaystyle {\sum }_{n=2}^{\infty }\left[\frac{1}{{\Vert {\bar{r}}_d\Vert }^3}-A_n^{\left(2\right)}{\left(\frac{1}{\Vert {\bar{r}}_d\Vert }\right)}^{n+3}\right]{\bar{r}}_d}\end{array}$ (7)

Let $\left(\mathcal{S}:{\hat{i}}_2,{\hat{j}}_2,{\hat{k}}_2\right)$ is a synodic coordinate system with an origin that lies in the centre of mass of the smaller primary which is defined as follows:

${\hat{i}}_2=-\frac{{\bar{r}}_{12}}{r_{12}}$, ${\hat{j}}_2={\hat{k}}_2\times {\hat{i}}_2$, ${\hat{k}}_2=-\frac{{\bar{h}}_2}{h_2}$

where ${\bar{h}}_2={\bar{r}}_{12}\times {\stackrel{\dot{}}{\bar{r}}}_{12}$ is the angular momentum of the smaller primary w.r.t. the bigprimary. Let this frame rotates with angular velocity ${\bar{\omega }}_2={\omega }_2{\hat{k}}_2$ w.r.t the inertial frame $\left(I:X,Y,Z\right)$. Then, the equation of motion of the chief w.r.t. smaller primary is:

$\begin{array}{c}{\left[{\ddot{\stackrel{¯}{r}}}_c\right]}_{\mathcal{S}}={\left[{\ddot{\stackrel{¯}{r}}}_c\right]}_I-2{\bar{\omega }}_2\times {\left[{\stackrel{\dot{}}{\bar{r}}}_c\right]}_{\mathcal{S}}-{\left[{\stackrel{\dot{}}{\bar{\omega }}}_2\right]}_{\mathcal{S}}\times {\bar{r}}_c-{\bar{\omega }}_2\times \left({\bar{\omega }}_2\times {\bar{r}}_c\right)\\ =-q_1\left(1-\mu \right)\left(\frac{\left({\bar{r}}_{12}+{\bar{r}}_c\right)}{{\Vert {\bar{r}}_{12}+{\bar{r}}_c\Vert }^3}+\frac{{\bar{r}}_{12}}{{\Vert {\bar{r}}_{12}\Vert }^3}-{\displaystyle {\sum }_{n=2}^{\infty }{A}_n^{\left(1\right)}}\left[\frac{\left({\bar{r}}_{12}+{\bar{r}}_c\right)}{{\Vert {\bar{r}}_{12}+{\bar{r}}_c\Vert }^{n+3}}+\frac{{\bar{r}}_{12}}{{\Vert {\bar{r}}_{12}\Vert }^{n+3}}\right]\right)\\ -q_2\mu {\displaystyle {\sum }_{n=2}^{\infty }\left[\frac{1}{{\Vert {\bar{r}}_c\Vert }^3}-A_n^{\left(2\right)}{\left(\frac{1}{\Vert {\bar{r}}_c\Vert }\right)}^{n+3}\right]{\bar{r}}_c}-2{\bar{\omega }}_2\times {\left[{\stackrel{\dot{}}{\bar{r}}}_c\right]}_{\mathcal{S}}-{\left[{\stackrel{\dot{}}{\bar{\omega }}}_2\right]}_{\mathcal{S}}\times {\bar{r}}_c\\ -{\bar{\omega }}_2\times \left({\bar{\omega }}_2\times {\bar{r}}_c\right)\end{array}$ (8)

Similarly, the motion of the deputy satellite will be:

$\begin{array}{c}{\left[{\ddot{\stackrel{¯}{r}}}_d\right]}_{\mathcal{S}}=-q_1\left(1-\mu \right)\left(\frac{\left({\bar{r}}_{12}+{\bar{r}}_d\right)}{{\Vert {\bar{r}}_{12}+{\bar{r}}_d\Vert }^3}+\frac{{\bar{r}}_{12}}{{\Vert {\bar{r}}_{12}\Vert }^3}-{\displaystyle {\sum }_{n=2}^{\infty }{A}_n^{\left(1\right)}}\left[\frac{\left({\bar{r}}_{12}+{\bar{r}}_d\right)}{{\Vert {\bar{r}}_{12}+{\bar{r}}_d\Vert }^{n+3}}+\frac{{\bar{r}}_{12}}{{\Vert {\bar{r}}_{12}\Vert }^{n+3}}\right]\right)\\ -q_2\mu {\displaystyle {\sum }_{n=2}^{\infty }\left[\frac{1}{{\Vert {\bar{r}}_d\Vert }^3}-A_n^{\left(2\right)}{\left(\frac{1}{\Vert {\bar{r}}_d\Vert }\right)}^{n+3}\right]{\bar{r}}_d}-2{\bar{\omega }}_2\times {\left[{\stackrel{\dot{}}{\bar{r}}}_d\right]}_{\mathcal{S}}-{\left[{\stackrel{\dot{}}{\bar{\omega }}}_2\right]}_{\mathcal{S}}\times {\bar{r}}_d\\ -{\bar{\omega }}_2\times \left({\bar{\omega }}_2\times {\bar{r}}_d\right)\end{array}$ (9)

where

${\left[{\stackrel{\dot{}}{\bar{r}}}_c\right]}_{\mathcal{S}}={\left[{\stackrel{\dot{}}{\bar{r}}}_c\right]}_I-{\bar{\omega }}_2\times {\bar{r}}_c$

2.2. The Relative Motion in Perturbed Restricted Three Body Problem

To describe the relative motion of the deputy w.r.t. the chief, let $\left(\mathcal{L}:{\hat{i}}_3,{\hat{j}}_3,{\hat{k}}_3\right)$ is a Local Vertical Local Horizontal frame “LVLH” with an origin that lies in the chief centre of mass and is defined as (As is clear in Figure 2):

${\hat{i}}_3={\hat{j}}_3\times {\hat{k}}_3$, ${\hat{j}}_3=-\frac{{\bar{h}}_c}{h_c}$, ${\hat{k}}_3=-\frac{{\bar{r}}_c}{r_c}$ (10)

where ${\bar{h}}_c={\bar{r}}_c\times {\left[{\stackrel{\dot{}}{\bar{r}}}_c\right]}_{\mathcal{S}}$ is the chief angular momentum w.r.t. the second primary. The position of the deputy w.r.t. the smaller primary is given by (see Figure 2):

${\bar{r}}_d={\bar{r}}_c+\bar{\rho }$ (11)

Let $\left(\mathcal{L}:{\hat{i}}_3,{\hat{j}}_3,{\hat{k}}_3\right)$ rotates with angular velocity ${\bar{\omega }}_3$ w.r.t. the inertial frame $\left(I:X,Y,Z\right)$, then

${\left[{\stackrel{\dot{}}{\bar{r}}}_d\right]}_I={\left[{\stackrel{\dot{}}{\bar{r}}}_c\right]}_I+{\left[\stackrel{\dot{}}{\bar{\rho }}\right]}_I={\left[{\stackrel{\dot{}}{\bar{r}}}_d\right]}_I+{\left[\stackrel{\dot{}}{\bar{\rho }}\right]}_{\mathcal{L}}+{\bar{\omega }}_3\times \bar{\rho }$

${\left[{\ddot{\stackrel{¯}{r}}}_d\right]}_I={\left[{\ddot{\stackrel{¯}{r}}}_c\right]}_I+{\left[\ddot{\stackrel{¯}{\rho }}\right]}_{\mathcal{L}}+2{\bar{\omega }}_3\times {\left[\stackrel{\dot{}}{\bar{\rho }}\right]}_{\mathcal{L}}+{\left[{\stackrel{\dot{}}{\bar{\omega }}}_3\right]}_I\times \bar{\rho }+{\bar{\omega }}_3\times \left({\bar{\omega }}_3\times \bar{\rho }\right)$ (12)

Introducing (6) and (7) into (12), then

Figure 2. Local vertical local horizontal frame w.r.t. synodic frame.

$\begin{array}{c}{\left[\ddot{\stackrel{¯}{\rho }}\right]}_{\mathcal{L}}=-q_1\left(1-\mu \right)\left[\frac{1}{{\Vert {\bar{r}}_{12}+{\bar{r}}_d\Vert }^3}-{\displaystyle {\sum }_{n=2}^{\infty }{A}_n^{\left(1\right)}}\left[\frac{1}{{\Vert {\bar{r}}_{12}+{\bar{r}}_d\Vert }^{n+3}}\right]\right]\left({\bar{r}}_{12}+{\bar{r}}_d\right)\\ -q_2\mu {\displaystyle {\sum }_{n=2}^{\infty }\left[\frac{1}{{\Vert {\bar{r}}_d\Vert }^3}-A_n^{\left(2\right)}{\left(\frac{1}{\Vert {\bar{r}}_d\Vert }\right)}^{n+3}\right]{\bar{r}}_d}\\ +q_1\left(1-\mu \right)\left[\frac{1}{{\Vert {\bar{r}}_{12}+{\bar{r}}_c\Vert }^3}-{\displaystyle {\sum }_{n=2}^{\infty }{A}_n^{\left(1\right)}}\left[\frac{1}{{\Vert {\bar{r}}_{12}+{\bar{r}}_c\Vert }^{n+3}}\right]\right]\left({\bar{r}}_{12}+{\bar{r}}_c\right)\\ +q_2\mu {\displaystyle {\sum }_{n=2}^{\infty }\left[\frac{1}{{\Vert {\bar{r}}_c\Vert }^3}-A_n^{\left(2\right)}{\left(\frac{1}{\Vert {\bar{r}}_c\Vert }\right)}^{n+3}\right]{\bar{r}}_c}-W\left(x,y,z\right)\end{array}$ (13)

where

$\begin{array}{l}W\left(x,y,z\right)=2{\bar{\omega }}_3\times {\left[\stackrel{\dot{}}{\bar{\rho }}\right]}_{\mathcal{L}}+{\left[{\stackrel{\dot{}}{\bar{\omega }}}_3\right]}_I\times \bar{\rho }+{\bar{\omega }}_3\times \left({\bar{\omega }}_3\times \bar{\rho }\right)\\ =\left({\omega }_3^x{\omega }_3^{y}y+{\omega }_3^x{\omega }_3^{z}z-{\left({\omega }_3^y\right)}^{2}x-{\left({\omega }_3^z\right)}^{2}x+2\stackrel{\dot{}}{z}{\omega }_3^y-2\stackrel{\dot{}}{y}{\omega }_3^z+z{\stackrel{\dot{}}{\omega }}_3^y-y{\stackrel{\dot{}}{\omega }}_3^z\right){\hat{i}}_2\\ +\left({\omega }_3^x{\omega }_3^{y}x+{\omega }_3^y{\omega }_3^{z}z-{\left({\omega }_3^x\right)}^{2}y-{\left({\omega }_3^z\right)}^{2}y-2\stackrel{\dot{}}{z}{\omega }_3^x+2\stackrel{\dot{}}{x}{\omega }_3^z-z{\stackrel{\dot{}}{\omega }}_3^x+x{\stackrel{\dot{}}{\omega }}_3^z\right){\hat{j}}_2\\ +\left({\omega }_3^x{\omega }_3^{z}x+{\omega }_3^y{\omega }_3^{z}y-{\left({\omega }_3^x\right)}^{2}z-{\left({\omega }_3^y\right)}^{2}z+2\stackrel{\dot{}}{y}{\omega }_3^x-2\stackrel{\dot{}}{x}{\omega }_3^y+y{\stackrel{\dot{}}{\omega }}_3^x-x{\stackrel{\dot{}}{\omega }}_3^y\right){\hat{k}}_2\end{array}$

The angular velocity of “ $\mathcal{L}$ ” w.r.t “I” is given by:

${\bar{\omega }}_3={\bar{\omega }}_{\mathcal{L}/\mathcal{S}}+{\bar{\omega }}_2$

${\left[{\stackrel{\dot{}}{\bar{\omega }}}_3\right]}_{\mathcal{L}}={\left[{\stackrel{\dot{}}{\bar{\omega }}}_{\mathcal{L}/\mathcal{S}}\right]}_{\mathcal{L}}+{\left[{\stackrel{\dot{}}{\bar{\omega }}}_2\right]}_{\mathcal{L}}={\left[{\stackrel{\dot{}}{\bar{\omega }}}_{\mathcal{L}/\mathcal{S}}\right]}_{\mathcal{L}}+{\left[{\stackrel{\dot{}}{\bar{\omega }}}_2\right]}_{\mathcal{S}}-{\bar{\omega }}_3\times {\bar{\omega }}_2$ (14)

where ${\bar{\omega }}_{\mathcal{L}/\mathcal{S}}$ and ${\bar{\omega }}_2$ are angular velocities of $\mathcal{L}$ w.r.t. $\mathcal{S}$ and $\mathcal{S}$ w.r.t. I respectively. To determine ${\bar{\omega }}_{\mathcal{L}/\mathcal{S}}$ and ${\stackrel{\dot{}}{\bar{\omega }}}_{\mathcal{L}/\mathcal{S}}$, a simple scheme based on the time derivatives of LVLH w.r.t. the synodic frame “ $\mathcal{S}$ ” [9] [24]:

${\left[{\widehat{\stackrel{\dot{}}{i}}}_3\right]}_{\mathcal{S}}={\bar{\omega }}_{\mathcal{L}/\mathcal{S}}\times {\hat{i}}_3$, ${\left[{\widehat{\stackrel{\dot{}}{j}}}_3\right]}_{\mathcal{S}}={\bar{\omega }}_{\mathcal{L}/\mathcal{S}}\times {\hat{j}}_3$, ${\left[{\widehat{\stackrel{\dot{}}{k}}}_3\right]}_{\mathcal{S}}={\bar{\omega }}_{\mathcal{L}/\mathcal{S}}\times {\hat{k}}_3$ (15)

Multiply (15) by the relative unit vector as follows:

${\hat{i}}_3\times {\left[{\widehat{\stackrel{\dot{}}{i}}}_3\right]}_{\mathcal{S}}={\hat{i}}_3\times \left({\bar{\omega }}_{\mathcal{L}/\mathcal{S}}\times {\hat{i}}_3\right)={\bar{\omega }}_{\mathcal{L}/\mathcal{S}}-\left({\bar{\omega }}_{\mathcal{L}/\mathcal{S}}\cdot {\hat{i}}_3\right){\hat{i}}_3$

${\hat{j}}_3\times {\left[{\widehat{\stackrel{\dot{}}{j}}}_3\right]}_{\mathcal{S}}={\hat{j}}_3\times \left({\bar{\omega }}_{\mathcal{L}/\mathcal{S}}\times {\hat{j}}_3\right)={\bar{\omega }}_{\mathcal{L}/\mathcal{S}}-\left({\bar{\omega }}_{\mathcal{L}/\mathcal{S}}\cdot {\hat{j}}_3\right){\hat{j}}_3$

${\hat{k}}_3\times {\left[{\widehat{\stackrel{\dot{}}{k}}}_3\right]}_{\mathcal{S}}={\hat{k}}_3\times \left({\bar{\omega }}_{\mathcal{L}/\mathcal{S}}\times {\hat{k}}_3\right)={\bar{\omega }}_{\mathcal{L}/\mathcal{S}}-\left({\bar{\omega }}_{\mathcal{L}/\mathcal{S}}\cdot {\hat{k}}_3\right){\hat{k}}_3$

By summing up the previous equations we get

${\hat{i}}_3\times {\left[{\widehat{\stackrel{\dot{}}{i}}}_3\right]}_{\mathcal{S}}+{\hat{j}}_3\times {\left[{\widehat{\stackrel{\dot{}}{j}}}_3\right]}_{\mathcal{S}}+{\hat{k}}_3\times {\left[{\widehat{\stackrel{\dot{}}{k}}}_3\right]}_{\mathcal{S}}=2{\bar{\omega }}_{\mathcal{L}/\mathcal{S}}$ (16)

Considering (10), the time derivative of the unit vectors of the LVLH frame is:

${\left[{\widehat{\stackrel{\dot{}}{k}}}_3\right]}_{\mathcal{S}}=-\frac{1}{r_c}\left({\left[{\stackrel{\dot{}}{\bar{r}}}_c\right]}_{\mathcal{S}}\cdot {\hat{i}}_3\right){\hat{i}}_3$

${\left[{\widehat{\stackrel{\dot{}}{j}}}_3\right]}_{\mathcal{S}}=-\frac{{\bar{r}}_c}{h_c}\left({\left[{\ddot{\stackrel{¯}{r}}}_c\right]}_{\mathcal{S}}\cdot {\hat{j}}_3\right){\hat{i}}_3$

${\left[{\widehat{\stackrel{\dot{}}{i}}}_3\right]}_{\mathcal{S}}=\frac{{\bar{r}}_c}{h_c}\left({\left[{\ddot{\stackrel{¯}{r}}}_c\right]}_{\mathcal{S}}\cdot {\hat{j}}_3\right){\hat{j}}_3+\frac{1}{r_c}\left({\left[{\stackrel{\dot{}}{\bar{r}}}_c\right]}_{\mathcal{S}}\cdot {\hat{i}}_3\right){\hat{k}}_3$ (17)

Substitute from (17), into (16), then:

${\bar{\omega }}_{\mathcal{L}/\mathcal{S}}=\left[\begin{array}{c}0\\ -\frac{h_c}{r_c^2}\\ -\frac{r_c}{h_c^2}\left({\left[{\ddot{\stackrel{¯}{r}}}_c\right]}_{\mathcal{S}}\cdot {\bar{h}}_c\right)\end{array}\right]$ and ${\stackrel{\dot{}}{\bar{\omega }}}_{\mathcal{L}/\mathcal{S}}=\left[\begin{array}{c}0\\ -\left(\frac{{\stackrel{\dot{}}{h}}_c}{r_c^2}-\frac{2{\stackrel{\dot{}}{r}}_c{h}_c}{r_c^3}\right)\\ -\left(\frac{{\stackrel{\dot{}}{r}}_c}{h_c^2}-\frac{2{\stackrel{\dot{}}{h}}_c{r}_c}{h_c^3}\right){\bar{h}}_c\cdot {\left[{\ddot{\stackrel{¯}{r}}}_c\right]}_{\mathcal{S}}-\frac{r_c}{h_c^2}{\bar{h}}_c\cdot {\left[{\stackrel{⃛}{\stackrel{¯}{r}}}_c\right]}_{\mathcal{S}}\end{array}\right]$ (18)

where ${\left[{\ddot{\stackrel{¯}{r}}}_c\right]}_{\mathcal{S}}$ and the jerk ${\left[{\stackrel{⃛}{\stackrel{¯}{r}}}_c\right]}_{\mathcal{S}}$ can be obtained by Equation (8) and its differentiation. Equation (13) along with (14) and (18) represent the equation of relative motion of the deputy w.r.t. the chief in the frame of the perturbed restricted three-body problem. It noted that the equation of relative motion is a nonlinear 2^nd order differential equation with time-varying parameters which can be simplified assuming the circular case of the restricted three-body problem.

3. The Relative Motion in the Circular Restricted Three-Body Problem

Assuming that the two primaries revolve in a circular orbit around their common centre of mass, then the following simplifications will be considered [10]:

${\bar{r}}_{12}=-{\hat{i}}_2$, ${\left[{\stackrel{\dot{}}{r}}_{12}\right]}_{\mathcal{S}}=0$, ${\bar{\omega }}_2={\hat{k}}_2$, ${\stackrel{\dot{}}{\bar{\omega }}}_2=0$, ${\ddot{\stackrel{¯}{\omega }}}_2=0$

Consequently, the angular velocity and acceleration of the LVLH frame w.r.t. the inertial frame is simplified as follows:

${\bar{\omega }}_3={\bar{\omega }}_{\mathcal{L}/\mathcal{S}}+{\hat{k}}_2$

${\left[{\stackrel{\dot{}}{\bar{\omega }}}_3\right]}_{\mathcal{L}}={\left[{\stackrel{\dot{}}{\bar{\omega }}}_{\mathcal{L}/\mathcal{S}}\right]}_{\mathcal{L}}-{\bar{\omega }}_{\mathcal{L}/\mathcal{S}}\times {\hat{k}}_2$

For more simplifications, assume that both primaries are radiating and only the second zonal harmonic is considered. Then

$A_2^{\left(k\right)}=3j_2^{\left(k\right)}{R}_k^2{P}_2^{\left(k\right)}\left(\sin \delta \right)$

Let

$r_c=\Vert {\bar{r}}_c\Vert$

${\bar{r}}_c=-r_c{\hat{k}}_2$

$\bar{\rho }=x{\hat{i}}_2+y{\hat{j}}_2+z{\hat{k}}_2$

${\bar{r}}_d={\bar{r}}_c+\bar{\rho }=x{\hat{i}}_2+y{\hat{j}}_2+\left(z-r_c\right){\hat{k}}_2$

$G_1={\Vert {\bar{r}}_d\Vert }^{-3}={\Vert {\bar{r}}_c+\bar{\rho }\Vert }^{-1}={\left(x^2+y^2+{\left(z-r_c\right)}^2\right)}^{-\frac{1}{2}}$

$G_2={\Vert {\bar{r}}_{12}+{\bar{r}}_c+\bar{\rho }\Vert }^{-1}={\left({\left(x-1\right)}^2+y^2+{\left(z-r_c\right)}^2\right)}^{-\frac{1}{2}}$

$G_3={\Vert {\bar{r}}_{12}+{\bar{r}}_c\Vert }^{-1}={\left(1+r_c{}^2\right)}^{-\frac{1}{2}}$

Under these assumptions, the equation of relative motion (13) will be reduced to:

$\begin{array}{c}{\left[\ddot{\stackrel{¯}{\rho }}\right]}_{\mathcal{L}}=[-\left(1-\mu \right)q_1\left[G_2^3-A_2^{\left(1\right)}{G}_2^5\right]\left(x-1\right)-\mu q_2\left[G_1^3-A_2^{\left(2\right)}{G}_1^5\right]x\\ -\left(1-\mu \right)q_1\left[G_3^3-A_2^{\left(1\right)}{G}_3^5\right]-{\omega }_3^x{\omega }_3^{y}y-{\omega }_3^x{\omega }_3^{z}z+{\left({\omega }_3^y\right)}^{2}x\\ +{\left({\omega }_3^z\right)}^{2}x-2\stackrel{\dot{}}{z}{\omega }_3^y+2\stackrel{\dot{}}{y}{\omega }_3^z-z{\stackrel{\dot{}}{\omega }}_3^y+y{\stackrel{\dot{}}{\omega }}_3^z]{\hat{i}}_2\\ +[-\left(1-\mu \right)q_1\left[G_2^3-A_2^{\left(1\right)}{G}_2^5\right]y-\mu q_2\left[G_1^3-A_2^{\left(2\right)}{G}_1^5\right]y-{\omega }_3^x{\omega }_3^{y}x\\ -{\omega }_3^y{\omega }_3^{z}z+{\left({\omega }_3^x\right)}^{2}y+{\left({\omega }_3^z\right)}^{2}y+2\stackrel{\dot{}}{z}{\omega }_3^x-2\stackrel{\dot{}}{x}{\omega }_3^z+z{\stackrel{\dot{}}{\omega }}_3^x-x{\stackrel{\dot{}}{\omega }}_3^z]{\hat{j}}_2\end{array}$

$\begin{array}{c}+[-\left(1-\mu \right)q_1\left[G_2^3-A_2^{\left(1\right)}{G}_2^5\right]\left(z-r_c\right)-\mu q_2\left[G_1^3-A_2^{\left(2\right)}{G}_1^5\right]\left(z-r_c\right)\\ -\left(1-\mu \right)q_1\left[G_3^3-A_2^{\left(1\right)}{G}_3^5\right]r_c+\mu q_2\left[\frac{1}{r_c{}^2}-A_2^{\left(2\right)}{\left(\frac{1}{r_c}\right)}^4\right]-{\omega }_3^x{\omega }_3^{z}x\\ -{\omega }_3^y{\omega }_3^{z}y+{\left({\omega }_3^x\right)}^{2}z+{\left({\omega }_3^y\right)}^{2}z-2\stackrel{\dot{}}{y}{\omega }_3^x+2\stackrel{\dot{}}{x}{\omega }_3^y-y{\stackrel{\dot{}}{\omega }}_3^x+x{\stackrel{\dot{}}{\omega }}_3^y]{\hat{k}}_2\end{array}$ (19)

where ${\omega }_3^x,{\omega }_3^y,{\omega }_3^z$ and ${\stackrel{\dot{}}{\omega }}_3^x,{\stackrel{\dot{}}{\omega }}_3^y,{\stackrel{\dot{}}{\omega }}_3^z$ are the components of the angular velocity ${\bar{\omega }}_3$ and ${\stackrel{\dot{}}{\bar{\omega }}}_3$.

4. Solution Algorithm

Power series approaches are widely used to solve different celestial mechanics problems. Many authors depend on that algorithm to find an approximate solution for their problems [25] [26]. The Lie-integration method is one of the most famous power series algorithms that can be applied to find both displacement and velocity components of the deputy satellite. The method is outlined in the following three steps [16] [17] [26]:

Step I: Construction of the Lie operator

Let $\begin{array}{cc}x=g_1& \stackrel{\dot{}}{x}={\stackrel{\dot{}}{g}}_1=g_2\\ y=g_3& \stackrel{\dot{}}{y}={\stackrel{\dot{}}{g}}_3=g_4\\ z=g_5& \stackrel{\dot{}}{z}={\stackrel{\dot{}}{g}}_5=g_6\end{array}$

Then

${\stackrel{\dot{}}{g}}_1=g_2$

$\begin{array}{c}{\stackrel{\dot{}}{g}}_2=-\left(1-\mu \right)q_1\left[G_2^3-A_2^{\left(1\right)}{G}_2^5\right]\left(g_1-1\right)-\mu q_2\left[G_1^3-A_2^{\left(2\right)}{G}_1^5\right]g_1\\ -\left(1-\mu \right)q_1\left[G_3^3-A_2^{\left(1\right)}{G}_3^5\right]-{\omega }_3^{g_1}{\omega }_3^{g_3}{g}_3-{\omega }_3^{g_1}{\omega }_3^{g_5}{g}_5+{\left({\omega }_3^{g_3}\right)}^2{g}_1\\ +{\left({\omega }_3^{g_5}\right)}^2{g}_1-2g_6{\omega }_3^{g_3}+2g_4{\omega }_3^{g_5}-g_5{\stackrel{\dot{}}{\omega }}_3^{g_3}+g_3{\stackrel{\dot{}}{\omega }}_3^{g_5}\end{array}$

${\stackrel{\dot{}}{g}}_3=g_4$

$\begin{array}{c}{\stackrel{\dot{}}{g}}_4=-\left(1-\mu \right)q_1\left[G_2^3-A_2^{\left(1\right)}{G}_2^5\right]g_3-\mu q_2\left[G_1^3-A_2^{\left(2\right)}{G}_1^5\right]g_3-{\omega }_3^{g_1}{\omega }_3^{g_3}{g}_1\\ -{\omega }_3^{g_3}{\omega }_3^{g_5}{g}_5+{\left({\omega }_3^{g_1}\right)}^2{g}_3+{\left({\omega }_3^{g_5}\right)}^2{g}_3+2g_6{\omega }_3^{g_1}-2g_2{\omega }_3^{g_5}+g_5{\stackrel{\dot{}}{\omega }}_3^{g_1}-g_1{\stackrel{\dot{}}{\omega }}_3^{g_5}\end{array}$

${\stackrel{\dot{}}{g}}_5=g_6$

$\begin{array}{c}{\stackrel{\dot{}}{g}}_6=-\left(1-\mu \right)q_1\left[G_2^3-A_2^{\left(1\right)}{G}_2^5\right]\left(g_5-r_c\right)-\mu q_2\left[G_1^3-A_2^{\left(2\right)}{G}_1^5\right]\left(g_5-r_c\right)\\ -\left(1-\mu \right)q_1\left[G_3^3-A_2^{\left(1\right)}{G}_3^5\right]r_c+\mu q_2\left[\frac{1}{r_c{}^2}-A_2^{\left(2\right)}{\left(\frac{1}{r_c}\right)}^4\right]-{\omega }_3^{g_1}{\omega }_3^{g_5}{g}_1\\ -{\omega }_3^{g_3}{\omega }_3^{g_5}{g}_3+{\left({\omega }_3^{g_1}\right)}^2{g}_5+{\left({\omega }_3^{g_3}\right)}^2{g}_5-2g_4{\omega }_3^{g_1}+2g_2{\omega }_3^{g_3}-g_3{\stackrel{\dot{}}{\omega }}_3^{g_1}+g_1{\stackrel{\dot{}}{\omega }}_3^{g_3}\end{array}$

In general, the components of the equation of relative motion can be rewritten in a matrix form as:

$\stackrel{\dot{}}{g}=P+W$

where $\stackrel{\dot{}}{g}$, $P$ and $W$ are 6 × 1 matrices defined as follows:

$\stackrel{\dot{}}{g}={\left(\begin{array}{cccccc}{\stackrel{\dot{}}{g}}_1& {\stackrel{\dot{}}{g}}_2& {\stackrel{\dot{}}{g}}_3& {\stackrel{\dot{}}{g}}_4& {\stackrel{\dot{}}{g}}_5& {\stackrel{\dot{}}{g}}_6\end{array}\right)}^{\text{T}}$

and

$W={\left(\begin{array}{cccccc}{W}_1& W_2& W_3& W_4& W_5& W_6\end{array}\right)}^{\text{T}}$

With $W_1=W_3=W_5=0$. Based on (18), ${\omega }_3^{g_1}=0$, then

$W_2=-\left(1-\mu \right)q_1{G}_3^3+{\left({\omega }_3^{g_3}\right)}^2{g}_1+{\left({\omega }_3^{g_5}\right)}^2{g}_1-2g_6{\omega }_3^{g_3}+2g_4{\omega }_3^{g_5}-g_5{\stackrel{\dot{}}{\omega }}_3^{g_3}+g_3{\stackrel{\dot{}}{\omega }}_3^{g_5}$

$W_4=-{\omega }_3^{g_3}{\omega }_3^{g_5}{g}_5+{\left({\omega }_3^{g_5}\right)}^2{g}_3-2g_2{\omega }_3^{g_5}+g_5{\stackrel{\dot{}}{\omega }}_3^{g_1}-g_1{\stackrel{\dot{}}{\omega }}_3^{g_5}$

$\begin{array}{c}{W}_6=-{\omega }_3^{g_3}{\omega }_3^{g_5}{g}_3+{\left({\omega }_3^{g_3}\right)}^2{g}_5+2g_2{\omega }_3^{g_3}-g_3{\stackrel{\dot{}}{\omega }}_3^{g_1}+g_1{\stackrel{\dot{}}{\omega }}_3^{g_3}-\left(1-\mu \right)q_1{G}_3^3\\ +\mu q_2\left[\frac{1}{r_c{}^2}-A_2^{\left(2\right)}{\left(\frac{1}{r_c}\right)}^4\right]\end{array}$

However, $P$ is defined as:

$P=UV$

where is 6 × 6 matrix given by

$U=\left(\begin{array}{cccccc}{g}_2& 0& 0& 0& 0& 0\\ 0& \left(g_1-1\right)& 0& g_1& 0& 1\\ 0& 0& g_4& 0& 0& 0\\ 0& g_3& 0& g_3& 0& 0\\ 0& 0& 0& 0& g_6& 0\\ 0& \left(g_5-r_c\right)& 0& \left(g_5-r_c\right)& 0& r_c\end{array}\right)$

$V={\left(\begin{array}{cccccc}1& -\left(1-\mu \right)q_1\left[G_2^3-A_2^{\left(1\right)}{G}_2^5\right]& 1& -\mu q_2\left[G_1^3-A_2^{\left(2\right)}{G}_1^5\right]& 1& -\left(1-\mu \right)q_1\left[G_3^3-A_2^{\left(1\right)}{G}_3^5\right]\end{array}\right)}^{\text{T}}$

The Lie-Operator is defined as [16]:

$D=\frac{\text{d}}{\text{d}t}={\displaystyle {\sum }_{i=1}^6\frac{\partial }{\partial g_i}\frac{\text{d}{g}_i}{\text{d}t}}+\frac{\partial }{\partial t}$

For explicit time variables, then

$D={\displaystyle {\sum }_{i=1}^6{\stackrel{\dot{}}{g}}_i\frac{\partial }{\partial g_i}}$

For $j=1$, assuming that $P_{ij}=P_i$, $W_{ij}=W_i$ and ${\stackrel{\dot{}}{g}}_{ij}={\stackrel{\dot{}}{g}}_i$ then ${\stackrel{\dot{}}{g}}_i=\left(P_i+W_i\right)$ and

$\begin{array}{c}D={\displaystyle {\sum }_{i=1}^6\left(P_i+W_i\right)\frac{\partial }{\partial g_i}}\\ =\left(P_1+W_1\right)\frac{\partial }{\partial g_1}+\left(P_2+W_2\right)\frac{\partial }{\partial g_2}+\left(P_3+W_3\right)\frac{\partial }{\partial g_3}\\ +\left(P_4+W_4\right)\frac{\partial }{\partial g_4}+\left(P_5+W_5\right)\frac{\partial }{\partial g_5}+\left(P_6+W_6\right)\frac{\partial }{\partial g_6}\end{array}$ (20)

Step II: Construction of the recurrence relations for each variable

Applying the Lie operator Equation (20) on $g_i$, we obtain that

$D^n{g}_i=D^{n-1}\left(P_i+W_i\right)$

1) The recurrence formulas for $D{P}_i$

$P_i={\displaystyle {\sum }_{j=1}^6{U}_{ij}{V}_{ji}}$

$D{P}_i={\displaystyle {\sum }_{j=1}^6\left[U_{ij}D{V}_{ji}+V_{ji}D{U}_{ij}\right]}$

Generally,

$D^n{P}_i={\displaystyle {\sum }_{m=0}^n{\displaystyle {\sum }_{j=1}^6\left[D^{n-m}{V}_{ji}\right]\left[D^m{U}_{ij}\right]}}$ (21)

a) The recurrence formulas for $D{V}_{ji}$

$V_{21}=-\left(1-\mu \right)q_1\left[G_2^3-A_2^{\left(1\right)}{G}_2^5\right]$

$D{V}_{21}=T_1$

$\begin{array}{c}{T}_1=\left(1-\mu \right)q_1\left[3G_2^5-5A_2^{\left(1\right)}{G}_2^7\right]\left[\left(g_1-1\right)D{g}_1+g_{3}D{g}_3+\left(g_5-r_c\right)D{g}_5\right]\\ =\left(1-\mu \right)q_1\left[3G_2^5-5A_2^{\left(1\right)}{G}_2^7\right]\left(g_1{g}_2+g_3{g}_4+\left(g_5-r_c\right)g_6\right)\end{array}$

$V_{41}=-\mu q_2\left[G_1^3-A_2^{\left(2\right)}{G}_1^5\right]$

$D{V}_{41}=T_2$

$\begin{array}{c}{T}_2=\mu q_2\left[3G_1^5-5A_2^{\left(2\right)}{G}_1^7\right]\left(g_{1}D{g}_1+g_{3}D{g}_3+\left(g_5-r_c\right)D{g}_5\right)\\ =\mu q_2\left[3G_1^5-5A_2^{\left(2\right)}{G}_1^7\right]\left(g_1{g}_2+g_3{g}_4+\left(g_5-r_c\right)g_6\right)\end{array}$

The higher powers of $D{V}_{21}$ and $D{V}_{41}$ are computed as

$D^n{V}_{21}=D^{n-1}{T}_1$ and $D^n{V}_{41}=D^{n-1}{T}_2$ (22)

The recurrence relation of the rest of the elements of $V$ is zero.

b) The recurrence formulas for $D{U}_{ij}$

By definition of $g_i$, its noted that $D^n{g}_i=D^{n-1}{g}_{i+1}$

$U_{11}=g_2$ then $D^n{U}_{11}=D^n{g}_2$

$U_{22}=g_1-1$ and $U_{24}=g_1$ then $D^n{U}_{22}=D^n{U}_{24}=D^n{g}_1=D^{n-1}{g}_2$

$U_{33}=g_4$ then $D^n{U}_{33}=D^n{g}_4$

$U_{42}=U_{44}=g_3$

then

$D^n{U}_{42}=D^n{g}_3=D^{n-1}{g}_4$ (23)

$U_{55}=g_6$

Then

$D^n{U}_{55}=D^n{g}_6$ (24)

$U_{62}=U_{64}=g_5-r_c$

Then

$D^n{U}_{62}=D^n{U}_{64}=D^n{g}_5=D^{n-1}{g}_6$ (25)

The recurrence relation of the rest of the elements of $U$ is zero.

2) The recurrence formulas for $D{W}_{ij}$

$D^n{W}_1=D^n{W}_2=D^n{W}_3=0$ (26)

$D{W}_2={\beta }_2$

$\begin{array}{c}{\beta }_2=\left[{\left({\omega }_3^{g_3}\right)}^2+{\left({\omega }_3^{g_5}\right)}^2\right]D{g}_1-2{\omega }_3^{g_3}D{g}_6+2{\omega }_3^{g_5}D{g}_4-{\stackrel{\dot{}}{\omega }}_3^{g_3}D{g}_5+{\stackrel{\dot{}}{\omega }}_3^{g_5}D{g}_3\\ =\left[{\left({\omega }_3^{g_3}\right)}^2+{\left({\omega }_3^{g_5}\right)}^2\right]g_2-2{\omega }_3^{g_3}D{g}_6+2{\omega }_3^{g_5}D{g}_4-{\stackrel{\dot{}}{\omega }}_3^{g_3}{g}_6+{\stackrel{\dot{}}{\omega }}_3^{g_5}{g}_4\end{array}$

$D{W}_4={\beta }_4$

$\begin{array}{c}{\beta }_4=-{\omega }_3^{g_3}{\omega }_3^{g_5}D{g}_5+{\left({\omega }_3^{g_5}\right)}^{2}D{g}_3-2{\omega }_3^{g_5}D{g}_2-{\stackrel{\dot{}}{\omega }}_3^{g_5}D{g}_1\\ =-{\omega }_3^{g_3}{\omega }_3^{g_5}{g}_6+{\left({\omega }_3^{g_5}\right)}^2{g}_4-2{\omega }_3^{g_5}D{g}_2-{\stackrel{\dot{}}{\omega }}_3^{g_5}{g}_2\end{array}$

$D{W}_6={\beta }_6$

${\beta }_6=-{\omega }_3^{g_3}{\omega }_3^{g_5}D{g}_3+{\left({\omega }_3^{g_3}\right)}^{2}D{g}_5+2{\omega }_3^{g_3}D{g}_2+{\stackrel{\dot{}}{\omega }}_3^{g_3}D{g}_1$

The higher powers of $D{W}_2,D{W}_4$ and $D{W}_6$ are computed as

$D^n{W}_2=D^{n-1}{\beta }_2$, $D^n{W}_4=D^{n-1}{\beta }_4$ and $D^n{W}_6=D^{n-1}{\beta }_6$ (27)

The recurrence relation of the rest of the elements of $W$ is zero.

Step III: Find the Lie-series solution

The solution is given by:

$\bar{X}\left(g_1,g_2,g_3,g_4,g_5,g_6\right)={\left[\left\{\exp \left[\left(t-t_o\right)D\right]\right\}X\right]}_{\bar{X}={\bar{X}}_o}={\displaystyle {\sum }_{k=0}^{\infty }{\left[D^k\bar{X}\right]}_{{\bar{X}}_o}\frac{{\tau }^k}{k!}}$

Then

$\begin{array}{c}{g}_i={\displaystyle {\sum }_{k=0}^{\infty }\frac{{\tau }^k}{k!}{\left[D^k{g}_i\right]}_{g_i=g_o}}\\ ={\left[g_i\right]}_{g_o}+{\displaystyle {\sum }_{k=1}^{\infty }\frac{{\tau }^k}{k!}{\left[D^{k-1}\left(P_i+W_i\right)\right]}_{g_i=g_o}}\\ ={\left[g_i\right]}_{g_i=g_o}+\tau {\left[P_i\right]}_{g_i=g_o}+\tau {\left[W_i\right]}_{g_i=g_o}+{\displaystyle {\sum }_{k=2}^{\infty }\frac{{\tau }^k}{k!}{\left[D^{k-1}{W}_i\right]}_{g_i=g_o}}\\ +{\displaystyle {\sum }_{k=2}^{\infty }{\displaystyle {\sum }_{j=1}^6{\displaystyle {\sum }_{m=0}^n\frac{{\tau }^k}{k!}\left[D^{k-m-1}{G}_{ji}\right]{\left[D^m{U}_{ij}\right]}_{g_i=g_o}}}}\end{array}$ (28)

5. Numerical Application

Consider a circular restricted three-body problem in the Earth-Moon system where the period of the Moon about the Earth is 27.23 days. The parameters of the primaries, the properties of the chief orbit and the initial condition of the deputy in the LVLH frame are tabulated in Table 1 and Table 2.

The calculations have been made using 10⁻² step size and five terms calculations over 24 hours. To assess the motion of the deputy employing the concepts of the restricted three-body problem, the solution algorithm is applied for both the unperturbed (classical) and the perturbed cases. A set of curves represent the relative motion of the deputy satellite w.r.t. the chief in both cases.

As is clear in Figures 3-5, in all panes of motion there is a significant difference in the components of the deputy position between the classical and perturbed cases. Consequently, the magnitude of its relative position vector is changed as is clear in Figure 6 where the solid curve represents the perturbed motion and the dotted curve represents the unperturbed motion and their difference is represented in Figure 7.

Table 1. The parameters of the primaries.

Table 2. Initial Conditions of the chief and deputy satellites.

Figure 3. Motion of the deputy satellite w.r.t. the chief in x-y plane for both the unperturbed and perturbed cases of the three-body problem.

Figure 4. Motion of the deputy satellite w.r.t. the chief in x-z plane for both the unperturbed and perturbed cases of the three-body problem.

Figure 5. Motion of the deputy satellite w.r.t. the chief in z-y plane for both the unperturbed and perturbed cases of the three-body problem.

Figure 6. Comparison between the deputy relative motion for the unperturbed and perturbed cases of the three-body problem.

Figure 7. Relative distance difference between the unperturbed and perturbed cases.

6. Conclusion

The motion of the deputy satellite w.r.t. the chief of the formation is modelled in frame of the perturbed restricted three-body problem using the LVLH frame. The model is a more accurate formulation compared with the previous work where it considers the effects of radiation of both primaries in addition to their oblateness second zonal harmonic. The system is simplified assuming the circular problem and solved numerically using the Lie series approach. The solution is tested using suitable initial conditions and applied for both the classical and perturbed restricted three-body problems. By comparing the two cases, the results show that there is a significant difference in the deputy relative distance. Consequently, the model will be suitable for a more accurate study of the different space mission issues (e.g. satellite rendezvous and control).

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Clohessy, W.H. and Wiltshire, R.S. (1960) Terminal Guidance System for Satellite Rendezvous. Journal of the Aerospace Sciences, 27, 653-674. https://doi.org/10.2514/8.8704
[2]	Tschauner, J. and Hempel, P. (1965) Rendezvous zu einem in elliptischer Bahn umlaufenden Ziel [Rendezvous with a Target in an Elliptical Orbit]. Astronautica Acta, 11, 104-109.
[3]	Di Mauro, G., Lawn, M. and Bevilacqua, R. (2018) Survey on Guidance Navigation and Control Requirements for Spacecraft Formation-Flying Missions. Journal of Guidance, Control, and Dynamics, 41, 583-602. https://doi.org/10.2514/1.G002868
[4]	Sullivan, J., Grimberg, S. and D’Amico, S. (2017) Comprehensive Survey and Assessment of Spacecraft Relative Motion Dynamics Models. Journal of Guidance, Control, and Dynamic, 40, 1837-1860. https://doi.org/10.2514/1.G002309
[5]	Djojodihardjo, H. (2014) Formation Flying Dynamics of Micro-Satellites near Equatorial Low Orbits under the Influence of J2 and J3. 2nd IAA Conference on Dynamics and Control of Space Systems (DYCOSS), Roma, 24-26 March 2014, 1-15. https://doi.org/10.4304/jait.5.1.5-14
[6]	Schaub, H. and Junkins, J.L. (2009) Analytical Mechanics of Space Systems. 2nd Edition, American Institute of Aeronautics and Astronautics, Reston. https://doi.org/10.2514/4.867231
[7]	Roscoe, C.W.T., Vadali, S.R. and Alfriend, K.T. (2013) Third-Body Perturbation Effects on Satellite Formations. The Journal of the Astronautical Sciences, 60, 408-433. https://doi.org/10.1007/s40295-015-0057-x
[8]	Fehse, W. (2003) Automated Rendezvous and Docking of Spacecraft. Cambridge University Press, New York. https://doi.org/10.1017/CBO9780511543388
[9]	Franzini, G. and Innocenti, M. (2019) Relative Motion Dynamics in the Restricted Three-Body Problem. Journal of Spacecraft and Rockets, 56, 1-16. https://doi.org/10.2514/1.A34390
[10]	Qian, Y.J., Yang, X.D., Jing, W.X. and Zhang, W. (2018) An Improved Numerical Method for Constructing Halo/Lissajous Orbits in a Full Solar System Model. Chinese Journal of Aeronautics, 31, 1362-1374. https://doi.org/10.1016/j.cja.2018.03.006
[11]	Ferrari, F. and Lavagna, M. (2014) Formation Flying and Relative Dynamics under the Circular Restricted Three-Body Problem Formulation. Conference: 24th AAS/AIAA Space Flight Mechanics Meeting, Santa Fe, 26-30 January 2014, 185-204.
[12]	Heritier A. and Howell K.C. (2012) Regions near the Libration Points Suitable to Maintain Multiple Space. AIAA/AAS Astrodynamics Specialist Conference, Minneapolis, 13-16 August 2012, 1-20. https://doi.org/10.2514/6.2012-4666
[13]	Gomez, G., Marcote, M., Masdemont, J.J. and Mondelo, J.M. (2006) Natural Configurations and Controlled Motions Suitable for Formation Flying. Advances in the Astronautical Sciences, 123, 513-1530.
[14]	Barden B.T. and Howell, K. (1998) Formation Flying in the Vicinity of Libration Point Orbits. Advances in the Astronautical Sciences, 99, 969-988.
[15]	Szebeheley, V. (1976) Theory of Orbits. Academic Press, New York.
[16]	Magda, D. (1984) Integration of the Elliptic Restricted Three Body Problem with Lie Series. Journal of Celestial Mechanics, 34, 145-154. https://doi.org/10.1007/BF01235797
[17]	Abouelmagd, E.I., Guirao, J.L. and Mostafa, A. (2014) Numerical Integration of the Restricted Three-Body Problem with Lie Series. Astrophysics and Space Science, 354, 369-378. https://doi.org/10.1007/s10509-014-2107-4
[18]	Sanam Suraj, M., Amit, M., Krishan, P. and Deepak, M. (2020) Stability of the Artificial Equilibrium Points in the Low-Thrust Restricted Three-Body Problem When the Smaller Primary Is an Oblate Spheroid. Nonlinear Dynamics and Systems Theory, 20, 439-450.
[19]	Arredondo, J.A. and Jimenez-Cardenas, J. (2021) About the Restricted Three-Body Problem with the Schwarszchild-de Sitter Potential. Nonlinear Dynamics and Systems Theory, 21, 126-137.
[20]	Tang, H. (2019) Koopman Reduced Order Control for Three Body Problem. Modern Mechanical Engineering, 9, 20-29. https://doi.org/10.4236/mme.2019.91003
[21]	Abouelmagd, E.I. and Guirao, J.L. (2016) On the Perturbed Restricted Three-Body Problem. Applied Mathematics and Nonlinear Sciences, 1, 123-144. https://doi.org/10.21042/AMNS.2016.1.00010
[22]	Abouelmagd, E.I., Alhothuali, M., Guirao, J.L. and Malaikah, H. (2015) The Effect of Zonal Harmonic Coefficients in the Framework of the Restricted Three-Body Problem. Advances in Space Research, 55, 1660-1672. https://doi.org/10.1016/j.asr.2014.12.030
[23]	Montenbruck, O. and Gill, E. (2001) Satellite Orbits. Models, Methods, Applications. Corrected Second Printing. Springer-Verlag, Berlin. https://doi.org/10.1007/978-3-642-58351-3
[24]	Casotto, S. (2016) The Equations of Relative Motion in the Orbital Reference Frame. Celestial Mechanics and Dynamical Astronomy, 124, 215-234. https://doi.org/10.1007/s10569-015-9660-1
[25]	Alghamdi, M.H. and Alshaery, A.A. (2020), A Series Solution Approach to the Circular Restricted Gravitational Three-Body Dynamical Problem. Journal of Applied Mathematics and Physics, 8, 2703-2712. https://doi.org/10.4236/jamp.2020.812200
[26]	Rajib, M. (2020) Lie-Series Solution of Restricted Three-Body Problem: Application to Binary Stellar Systems. The Journal of the Astronautical Sciences, 67, 59-76. https://doi.org/10.1007/s40295-019-00172-5