1. Introduction
1.1. Importance of Relative Motion Satellites
Relative motion is considered one of the fundamental problems in orbital mechanics, it has different applications that arise in the literature of astrodynamics for rendezvous and formation of flying satellites. The first idea concerning this problem was explained by the author of [
1] in the late 19th century, who proposed analyzing the Moon’s motion. His objective was to establish a more mathematically sound method for developing tables of lunar motion, this was at the time based on “practical astronomy rather than of mathematics” in his talk.
The first practical space of relative motion was in the field of intercept and space rendezvous throughout the late 1950s and extends up to the current time. In the intercept model, a chase satellite is compelled to move in such a way that its orbit intersects the trajectory of the chief satellite at a certain time. In the rendezvous problem the relative velocities of both vehicles must be pushed to zero at the intersection instant; thereby, a docking or berthing procedure or other activities may be carried out.
The authors in [
2] analyzed the relative motion within the frame of developing a guidance algorithm for a rendezvous task with an assumption that the chief satellite moves in a circular orbit, and acts as a control center. It sends messages with relative location and velocity data to slave satellites. Based on that data, the satellites perform a docking maneuver and rendezvous by using an on-board propulsion system.
The comprehension of relative motion satellites is substantial for two significant applications of orbital mechanics:
The first problem investigates the designing of centralization and decentralization algorithms to dock a satellite successfully, for instance, the space shuttle with the International Space Station. These algorithms demand the modeling of relative motion satellites, without needing permanent communication with ground-based stations. Rendezvous orbital dynamics and control is one of the key elements of space flight technology for operating space rendezvous and docking missions. Rendezvous demands a specific congruent of the orbital velocities and position vectors of the two vehicles, permitting them to stay at a constant separation during orbital station-keeping. It may or may not be followed by docking or berthing, procedures which bring the spacecraft into physical connection and generate a link between them. The same rendezvous method could be used for a spacecraft “landing” on a natural celestial body with a weak gravitational field, for example, landing on one of the Martian moons will demand the same conveniency of orbital velocities, followed by a “descent” that shares some similarities with docking.
The second problem is related to the formation of a satellite system, where the dynamics of the system are exploited to place a number of satellites in space. Formation of satellites is a novel concept that distributes a large spacecraft function among several small satellites, which are less expensive and cooperative. Most applications to relative motion are in the formation satellites, which have a considerable significance, because the use of a large number of satellites with low impact satellites running, in a concrete style, could introduce a better outcome than a single high implementation satellite. Furthermore, these formations do not require a big cost, which means a great opportunity for space mission success with more resiliency. Satellite formations have good benefits for Earth observation of space missions, where the group distribution of low resolution devices, operating in conjunction with each other, may provide a higher overall information quality than a single high resolution device.
1.2. Mathematical Models of Relative Motion Satellites
The relative motion between two satellites is described by several mathematical models [
3,
4,
5,
6]. The “Hill–Clohessy–Wiltshire model (HCW)” is as a classical approach described by a set of linearized, time-invariant ordinary differential equations. However, this model assumes a circular chief orbit and it is valid only for small initial separations between the satellites.
An assumption of an elliptical shape of the chief satellite orbit results in time-varying differential equations. In this work such a model is represented by the Tschauner–Hempel (TH) equations. However, since the applicability of this model is also limited to small initial separations between the satellites, this investigation deploys it only for comparative purposes.
In cases where a chief satellite (target or passive vehicle) is moving in a circular orbit, the orbital relative motion can be described using HCW equations [
2]. A principal assumption for the application of HCW equations is a distance of separation between chief and deputy (chase or active vehicle) satellites of less than 1 km.
The authors of [
7] were the first to find a closed-form solution for linearized relative motion within the frame of elliptic orbits. The authors of [
8] found a singularity in Lawden’s solution and then constructed a state-transition matrix that depends explicitly on the true and the eccentric anomaly, overcoming the singularity. We find a closed-form solution for a linearized relative motion within the frame of elliptic orbits in [
7] for the first time. The author of [
8] found a singularity in Lawden’s solution. Afterwards, they proposed a new state-transition matrix without a singularity by explicitly using the true and the eccentric anomaly data. They derived the state transition matrix (STM) of the (TH) equations for orbits with small eccentricity, based mainly on the solutions of the HCW equations using a perturbation theory. In the framework of improving the obtained solution, approximation solutions of the second–order relative motion equations set in spherical coordinates are constructed [
9].
Several versions of the relative motion dynamics have been developed to involve the perturbations effect [
10,
11,
12,
13]. The inhomogeneous gravitational field generated by the presence of the non-sphericity of objects, such as the Earth’s oblateness, is based on an assumption that the gravitational field is a sum of zonal harmonics coefficients; the authors of [
10] derived a set of relative dynamics time-variant equations under the effect of the zonal harmonic parameter. An analytical solution of relative motion subject to Lorentz-force perturbations is constructed in [
14]. Within the frame of developing the analysis of relative motion using state transition, matrices for the elliptical case are found, by the authors in [
15,
16], to be under the perturbation effect.
Furthermore, an outstanding piece of work studied the reformulation of relative motion based on non-linear system dynamics; see [
17]. Another excellent piece of work proposed to analyze relative motion dynamics, control techniques and specific applications for formation flying, deployment, station keeping, and reconfiguration; see [
18,
19,
20,
21,
22].
The aim of the proposed paper is to find the Periodic solutions of nonlinear relative motion using the perturbation technique of the Lindstedt–Poincaré method, which enable us to remove the secular terms, which will be growth with increasing time and will go to infinity with increasing time to infinity. However, there are some numerical methods, which can be used to estimate such solutions; for example: The Verrlet method, the step size control algorithm and the explicit Runge–Kutta method; for details, comparisons and effective methods, see [
23,
24]. However, we emphasize that the numerical methods are not used to verify the analytical solutions in most cases, and vice versa. In addition, there is no full warranty for convergence of numerical solutions in most cases. Thus, we are interested to find the periodic solution of non-linear relative motion, which is more realistic in practical applications by using the perturbation technique of Lindstedt–Poincaré.
2. Formulation of Motion Model
We consider a set of two satellites model, the first is called target satellite (or chief or passive vehicle) and the second is called chaser (or deputy or active vehicle). In general, this model is implemented to accomplish some space mission such as space rendezvous. The mission between two space crafts or satellites is successful when the two satellites have the same velocity and position vectors at a certain time. In the initiates of rendezvous time the two parts in a rendezvous sequence are:
Phasing of maneuvers to accomplish a rendezvous in the timing sequence, which will bring the two satellites into close proximity.
Maneuvers of finalizing rendezvous that involve the relative motion between the two satellites for rendezvous and docking.
Now we will assume that the inertial references frame with origin, which coincides with the center of the main body, say Earth’s center; the orthogonal unit vector is denoted by
. The vector
in the direction of the equinoxes and the vectors
is perpendicular to it, while
passes through the North Pole. To describe a relative motion satellite in Local-Vertical-Local-Horizontal (LVLH) coordinates, we proposed that the orientation of these frame by the unit vector triad
, where the vector
lies in the chief radial direction,
lies in the direction of the chief angular momentum, and
completes the right-handed orthogonal triad, such that
; see
Figure 1.
According to
Figure 1, the deputy satellite relative position vector in LVLH coordinates can be written as:
where
,
, and
are the components of the relative position vector along the radial, long track, and out of plane, respectively.
Since
, then
where
represents the current radius of the chief’s satellite orbit. The LVLH frame rotation with angular velocity
is given by
where
and
, while the attached frame rotates with the angular velocity
. We have to note that
in which
and
f are the argument of perigee and the true anomaly, respectively, within the frame of two-body motion. The angular momentum magnitude of the chief satellite is expressed by
where
r is the trajectory of the chief (target) satellite, which may be circular or elliptical. On the other hand, the expression of the chief’s true anomaly rate
can be written in terms of a semi-latus rectum as
wherein the elliptical case
,
a are the semi-major axes and
e is the eccentricity,
is the gravitational constant and
. The angular velocity and the acceleration of the the LVLH reference frame are given by
Now, the general vectorial equations that model the relative motion in the LVLH reference frame are
where
is the sum of total external forces acting on the deputy and chief, respectively;
and
denote the first and second derivative with respect to time.
Substituting Equations (
1) and (
2) into Equation (
3), we obtain
Now we can write the right-hand side of Equation (
4) in the form
where
is the external forces acting on the deputy and
is the external forces acting on the chief satellite, which can be written in the case of the gravitational field being generated by a spherical object, such as
Substitution Equation (
5) into Equation (
4); then, the vector form or relative motion equations will be controlled by
The vector Equation (
6) can be degraded to
Equation (
7) represents the general equations of the satellite relative motion in LVLH coordinates.
2.1. General Circular Relative Motion Satellites
Now, we denote
and
as the potential functions for gravitational forces
and
, acting on the chief and deputy satellite, respectively; hence,
Note that
; Equation (
8) is a generating function for Legendre polynomials with argument
; consequently,
where
is the Legendre polynomial with degree
n, where
Then,
where
is the higher-order of a Legendre polynomial given by
Then the total external forces can be rewritten as
where
Since
is the gravitational acceleration acting on the deputy and
is the gravitational acceleration acting on the chief’s satellite, while
is the total relative gravitational acceleration, then
and
Utilizing Equations (
9) and (
10) with Equation (
4), the vector equation of relative motion can be rewritten as
Hence, Equation (
11) can be written as three differential equations of second order
Equation (
12) represents the general dynamical system of the relative motion satellite. However, we have to observe that the quantities in the right-hand sides,
, represent the components of perturbing gravitational accelerations, which are considered the main contribution for higher order of nonlinearity to a dynamical system of relative motion.
2.2. Linear Circular Relative Motion Satellites
In the case where the chief’s orbit is moving in circular orbits, then
and
are constants. Thereby, the rotation frame with these properties is referred to as Hill’s frame and rotates with angular velocity
, and the resulting equations are referred to as HCW equations, as in [
2]. In addition, if the perturbing acceleration or the higher order on nonlinearity terms (
) are neglected. Then the linearized relative motion equations will be governed by
Let
,
,
; thereby, the above equations can be written as
The above dynamical system is famous and is called the Hill–Clohessy–Wiltshire equations, and represents the linearity of relative satellite motion. Furthermore, the dynamical system of relative motion (HCW) is invariant under the symmetry .
3. Solutions of HCW Equations
In order to construct the solutions of the dynamical system which represent the linear relative motion satellite, we will assume that the deputy satellite starts its motion under the following set of the initial conditions
Thereby, after integration of the second equation in System (
13), we obtain
where
is a constant which can be determined from the initial conditions, so one obtains
and Equation (
15) will take the following form
Substituting
into the first equation in System (
13), we get
Let
, then
and
and, the Equation (
17) becomes
Hence, Equation (
18) has the solution
Then the final form solution is
Substituting Equation (
19) into Equation (
16) and integrating it, we get
where
the constant of the integration.
To complete our steps, we have to integrate the third equation in System (
13), which represents the out-plane motion; hence, we get
Using the obtained solutions of Equations (
19)–(
21), the solutions of System (
13) are summarized in the following formulas
where
and
are arbitrary constants that must match the initial conditions.
Here, we observe that the second formula in Solutions (
22) contains an unbounded term, which will grow to infinity with increasing time. To get bounded solutions, in particular, periodic solutions, the coefficient of the unbounded term must equal zero; these will give an additional condition to get periodic and bounded solutions; thereby, we will obtain the new conditions on the initial conditions as
So we rewrite the periodic solutions of System (
13) in the following formulas
Applying the first derivative to Equations (
24), we obtain the relative velocities in the form
If we apply the initial conditions on Equation (
24) and the condition in Equation (
23), then the constants of integration can be determined by
and
Equations (
24) and (
25) represent the relative positions and velocities of the chaser satellite for space rendezvous within the framework of linearized motion.
With the help of Equations (
24)–(
27), the initial amplitude of the chaser satellite for space rendezvous, i.e.,
is taken as
, with different values of phase angle. The projection of the periodic solution in the xy-plane is shown in
Figure 2, whereas the orbit in three dimensions is shown in
Figure 3 at different values of the phase angle.
It is clear from
Figure 2a that the analytical solution of linear relative motion in two dimensions is periodic, and the obtained solution by a numerical integration of the explicit Euler method using step size
at a set of initial condition
as
is also periodic; this solution is shown in
Figure 2b. In addition, the two solutions are symmetric due to the y-axis
.
The linear motion in three dimensions is also periodic due two both solutions. We remark that both analytical and numerical solutions of linear relative motion in three dimensions are also periodic. The analytical solution is given in
Figure 3a, while the solution by a numerical integration of the explicit Euler method using step size
at a set of initial condition
as
is given in
Figure 3b.
Nonlinear Circular Relative Motion Satellites
In the case of considering the perturbed potential
, and keeping the nonlinearity up to second order terms, the relative motion satellite equations of System (
12) are approximated to the following equations
The above dynamical system consists of the three second-order differential equations with second degree, which represent the relative motion of the deputy satellite.
Now we intend to find the periodic solutions of these equations to find the relative position and velocity of the deputy satellite to accomplish space rendezvous. The solutions of these equations will be generalized or extended for the linear solutions; then, we have to find the solutions under the previous conditions of linear relative motion. Unfortunately, there are no exact analytical solutions for these equations, and the exact solutions are not possible. Thereby, we must introduce techniques or perturbation methods, which can be used to evaluate approximate or semi-analytical solutions; in particular, we must construct periodic solutions. So, the next section will be devoted to analysis of a perturbation technique, which enables us to find an appropriate approximated analytical solution to perform a space rendezvous mission.
5. Periodic Solutions by the Lindstedt–Poincar é Technique
Now we will apply the Lindstedt–Poincaré technique on the relative satellite motion to find the periodic solution of the equations of motion of the second degree, say System (
30). For simplicity, we suppose that
; then, System (
30) can be written in the form
To account for the nonlinear dependence of the frequency, we will introduce the stretch variable,
; we denote
and
; then, System (
31) can be written in the form
Now we look for approximated solutions in the form of a power series as
and
The initial conditions of position components can be rewritten in the form
while the initial conditions of velocity components can also be rewritten in the form
Substituting Equation (
33) into Equation (
32) with the help of Equation (
34) and keeping the terms up to first order in
, we obtain
where the frequency is
The set of initial conditions in Equations (
35) and (
36) may be rewritten as
and
Equating each of the coefficients of
that have the same power in System (
37), using Conditions (
38) and (
39), we get the unperturbed equations system in the form
and the set of initial conditions is
while the equations which represent the first corrections are
with the initial conditions
To accomplish our procedure and find the approximated periodic solutions, we have to find solutions for both the dynamical systems in Equations (
40) and (
42)
Zero order solutions
According to the solutions of linear motion, the general periodic solutions of the unperturbed relative motion will be controlled by
With a help of Conditions (
41), the constants of integration are given by
and
First order solutions
Substituting Equation (
44) into Equation (
42), we obtain
where
,
,
and
are constants, which are calculated by the following relations
Equation (
45) consists of a dynamical system of second-order non-homogeneous differential equations. After integrating this system with the help of Conditions (
43), we observe that the terms with coefficients
,
, and
represent secular terms, as well as the extra term with coefficient
. Therefore, if these coefficients are nonzero, then the correction solutions will involve secular terms in the scale
. However, that is exactly what we need to avoid. Hence we will remove these terms by equating their coefficients by zeros. These procedures ensure that the solutions will not include secular terms and at least no secularities appear in the first two terms in the perturbation series. Furthermore, the terms with coefficients
,
, and
represent solutions for the associated homogeneity of Equation (
45); hence, we set the coefficients to zeros before integration.
In this context, the removal of the secular terms yields
Then the first correction will be controlled by
Therefore, the periodic solutions of the nonlinear relative motion satellite can be written in the form
where
.
Since we use the method of place keeping parameters, we have to take
; then, Solutions (
47) will be rewritten as
Substituting Equations (
44) and (
46) into Equation (
48), the periodic solutions of the nonlinear relative motion satellite, with the first corrections, can be estimated approximately by using the Lindstedt–Poincaré technique.
With the help of Equation (
48), the initial amplitude of the chaser satellite for space rendezvous, i.e.,
, works similar to in the previous section
with different values for the phase angle. The projection of the periodic solution in the xy-plane is shown in
Figure 4, whereas the orbit in three dimensions is shown in
Figure 5 with different values of phase angle.
A comparison between the nonlinear relative motion solutions of the Lindstedt–Poincaré technique obtained by the numerical integration of the explicit Euler method are shown in
Figure 4, in two dimensions, and in
Figure 5, in three-dimensional motion. While numerical solutions for linear motion are periodic, the numericality of nonlinear motion is not periodic; this can be observed in
Figure 4b and
Figure 5b.
The projections in the xy-plane of the periodic solution of linear and nonlinear motion are shown in
Figure 6 at different phases, whereas the differences between orbits in three dimensions of the linear and nonlinear motion at the same different phases angle with
are shown in
Figure 7.
The graphical investigations in both
Figure 6 and
Figure 7 show that the analytical solutions of linear and nonlinear motions are periodic; however, there are differences between the style of the two motions. The periodicity of the nonlinear relative motion solutions gives a considerable verification due to the solution obtained with the Lindstedt–Poincaré technique.
6. Conclusions
In this paper, we studied the relative motion of an outline of the space rendezvous problem, assuming that the chief satellite was in circular symmetric orbits. It was found that the equations that describe the motion of the deputy satellite with respect to the references fixed at the center of the chief satellite is nonlinear in the general case. Then the periodic solutions first for linear relative motion satellite were evaluated. Here a key role is the symmetry of the periodic solution joined with the fact that the system HCW is invariant under negative symmetry.
The Lindstedt–Poincaré technique was used to get the approximated periodic solutions for nonlinear differential equations of the second degree. This method is referred to by many sources as a method of successive approximations. In our study, we use it to find the first corrections of a relative motion satellite. In this context, the obtained solutions are more general compared to the linear solutions. The analysis of a linear relative motion satellite and its periodic solutions are constructed. The approximate periodic solutions of the nonlinear relative motion satellite are constructed using the Lindstedt–Poincaré technique.
Comparisons among the analytical solution of linear motion and the obtained solution by the numerical integration of the explicit Euler method are investigated, and we show that bth solutions are periodic in two and three dimensions. Moreover, comparisons among the obtained solutions of nonlinear relative motion using the Lindstedt–Poincaré technique and the numerical solution by integration of the explicit Euler method were studied. However, the numerical solutions are not periodic in both motions in two and three dimensions. Thus, the Lindstedt–Poincaré technique is recommended for designing the periodic solutions.