Close relative motion on distant retrograde orbits

Chihang YANG, Ming WANG, Hao ZHANG

a University of Chinese Academy of Sciences, Beijing 100049, China

b Key Laboratory of Space Utilization, Technology and Engineering Center for Space Utilization, Chinese Academy of Sciences, Beijing 100094, China

KEYWORDS Distant retrograde orbits;Ephemeris model;Floquet theory;Formation flight;Relative motion;Spacecraft;Three-body problem;Transfer singularity

Abstract Distant Retrograde Orbits (DROs) in the Earth-Moon system have great potential to support varieties of missions due to the favorable stability and orbital positions.Thus,the close relative motion on DROs should be analyzed to design formations to assist or extend the DRO missions.However,as the reference DROs are obtained through numerical methods,the close relative motions on DROs are non-analytical,which severely limits the design of relative trajectories.In this paper, a novel approach is proposed to construct the analytical solution of bounded close relative motion on DROs.The linear dynamics of relative motion on DRO is established at first.The preliminary forms of the general solutions are obtained based on the Floquet theory.And the general solutions are classified as different modes depending on their periodic components.A new parameterization is applied to each mode, which allows us to explore the geometries of quasi-periodic modes in detail.In each mode, the solutions are integrated as a uniform expression and their periodic components are expanded as truncated Fourier series.In this way,the analytical bounded relative motion on DRO is obtained.Based on the analytical expression,the characteristics of different modes are comprehensively analyzed.The natural periodic mode is always located on the single side of the target spacecraft on DRO and is appropriate to be the parking orbits of the rendezvous and docking.On the basis of quasi-periodic modes, quasi-elliptical fly-around relative trajectories are designed with the assistance of only two impulses per period.The fly-around formation can support observations to targets on DRO from multiple viewing angles.And the fly-around formation is validated in a more practical ephemeris model.

1.Introduction

Close relative motion is the dynamical foundation of close spacecraft formation flying.The close spacecraft formation flying is the essential technique of many space mission scenarios,including rendezvous and docking, fly-around tasks, and onorbit services.The close spacecraft formation flying originated from the research for rendezvous and docking on geocentric orbit by Clohessy and Wiltshire in 1960.1Since then, the rendezvous and docking have become the key stages of the onorbit assembly and supply, and have been widely used to manned and unmanned space missions.2Two typical examples are the astronaut and cargo transfers in the Apollo program3and the establishment and maintenance of the International Space Stations.4Compared with the rendezvous and docking scenarios, more varieties of relative trajectories can be used in the fly-around tasks for different applications, such as fault monitoring and function extension to the target spacecraft.In 2016, a small subsatellite was released from the Tiangong-2 space laboratory to fly around and take photos to the space laboratory.5And a 2 m optical telescope is planned to be launched in 2024 and fly around with the Chinese Space Station, which will facilitate the refueling and on-orbit upgrade of the telescope.6The on-orbit services are a type of integrated application based on rendezvous and docking, and fly-around techniques.In 2020,the success of the Mission Extension Vehicle (MEV) opened the new era of commercial on-orbit services,7which implies that the on-orbit services may be extensively applied to future space missions.As a result, the close spacecraft formation flying is not only playing an essential role for now but will benefit more missions in the future.

To serve spacecraft near the Earth, most close relative motion researches focus on the two-body dynamics.With the development of space technology, deep space, especially the Earth-Moon space, is showing more and more potentials.8–9There are plenty of periodic orbits in the Earth-Moon space,of which DRO are a type of valuable orbits and attract a lot of attention in recent years.10–11DROs are periodic and orbit the Moon with a far distance in the plane of the Moon’s orbit.DROs have excellent stability and can keep bounded for thousands of years in the ephemeris model without orbital maintenance.12Moreover, spacecraft on DRO will have a good coverage to the space near the Earth and the surface of the Moon.13With these advantages, DRO can support varieties of large-size missions,such as cislunar space station,13–14Mars sample return mission,15and asteroid redirect mission.16The close spacecraft formation flying is indispensable to the establishment,the maintenance,and the function extension of these large-size missions.Therefore, the close relative motion on DROs is investigated in this paper.

In general, two targets, the boundedness and the specific geometry configuration,are required by a close spacecraft formation flying.The boundedness is the essential requirement to prevent the drift and collision between spacecraft.And specific geometry configurations are further desired by different missions.These two targets require a comprehensive understanding to the close relative motion.For simplicity, linearized dynamics are the most commonly used models to study the close relative motion.For linearized relative motion near a reference orbit in a six-dimensional Hamiltonian system,the solution set contains six linearly independent solutions.Numerical approaches are infeasible for such a high-dimensional solution set, while an analytical solution can significantly simplify the design of varieties of relative trajectories.Therefore, it is desired to obtain an analytical form to explore the close relative motion comprehensively.

The linearized relative motions in the two-body problem have been thoroughly investigated with analytical solutions.The most well-known linearized model is the Hill-Clohessy-Wiltshire (HCW) equations proposed by Hill17and the Clohessy-Wiltshire team.1HCW equations describe the close relative motion near a circular orbit.Due to the concise form and constant coefficients of HCW equations, the analytical solution can be obtained by directly solving the ordinary differential equations.1The boundedness condition of HCW equations is a linear constraint of the initial relative states.18And the bounded solutions of HCW equations are periodic orbits around the reference orbit.The so-called general circular orbit is one of the bounded solutions of HCW equations and has been applied to design gravitational wave detection missions.19As for the close relative motion near elliptical orbits, the linearized model is derived independently by Lawden20and the Tschauner-Hempel team,21thus consequently called Lawden’s equations or Tschauner–Hempel (TH) equations.The boundedness condition of TH equations is also a linear but much more complex equation.22–23As each elliptical orbit has a time-varying angular velocity, the analytical solution of TH equations is more complicated than that of HCW equations and is obtained through changing variables and constructing integrals.24Similar to HCW equations, the bounded solution of TH equations is also periodic but usually twisted,thus having more complex geometries.25The geometries of the bounded solution of TH equations are investigated by several researchers.Jiang et al.discovered that each periodic solution of TH equations is on a particular quadric surface with proper parameterization.25And Jiang et al.also analyzed the spatial self-intersections of the bounded solution.25Besides,Bando and Ichikawa created several ellipses to understand the evolution of the in-plane projection of the solution26.

As for the three-body problems, most existing researches focused on close relative motion near Lagrange points or periodic orbits.Since the position of each Lagrange point is fixed in the Circular Restricted Three-Body Problem (CRTBP), the linearized relative model is formed as homogeneous linear ordinary differential equations with constant coefficients,which can be solved analytically.27The close relative motions near collinear Lagrange points are unstable, and the bounded relative motion can be obtained by removing the unstable exponential terms from the initial states.Collange and Leitner found that some relative orbits near the Sun-Earth L2point are planar,elliptical,and periodic.28And the close relative motion near triangular Lagrange points is stable, naturally bounded,and composed of a long-period component and a shortperiod component.By isolating the long-period and shortperiod components, Catlin and McLaughlin investigated the close relative motion of two spacecraft near the Earth-Moon triangular Lagrange points.29The existence of analytical solutions simplifies the analyses of close relative motion and further formation design near Lagrange points.However, the fixed positions of Lagrange points limit the applications of close formations near them.Compared with Lagrange points,the periodic orbits in CRTBP are much more abundant and thus can support various formation missions.30Unlike the periodic orbits in the two-body problem, the periodic orbits in CRTBP are obtained through numerical iterations.Hence,there do not exist analytical solutions for close relative motion on periodic orbits in CRTBP,which leads to difficulties for relative motion analysis.To address this matter,the Floquet theory is applied by Howell and Marchand to explore the solutions for close relative motion on halo orbits.31And the basic geometries of the solutions are further analyzed by Simanjuntak32and Zhou33et al.The Floquet theory helps to understand the relative motion by classifying the solutions as different modes,34while the obtained solution set is only semianalytical but not fully analytical.Simanjuntak et al.treated the solutions based on the Floquet theory as piecewise functions and numerically fitted the solutions with a modified Fourier series.32However, the fitted results by Simanjuntak et al.have a complex form and a relatively low accuracy, but need tens of thousands of samples, and thus are inconvenient to use.32Therefore, it is still required to construct a concise,convenient, and accurate analytical solution for close relative motion on periodic orbits in CRTBP.

Moreover, the dynamical structures near periodic orbits in CRTBP are much more complex than those in the two-body problem.In the two-body problem, each periodic orbit lies in a five-parameter family.A commonly used group of parameters to describe the family is the classical orbital elements, of which the mean/true anomaly indicates the phase of the spacecraft on the orbit,and the other five elements determine a periodic orbit in the family.If we treat the periodic orbit after a phase shift as a new orbit, six new periodic orbits can be obtained by slightly perturbing each orbital element of a particular periodic orbit in the two-body problem.It is noticed that only the semi-major axis in the orbital elements can adjust the orbit period.Hence, the perturbation of the semi-major axis leads to a new periodic orbit with a different period.In contrast,the periods of the other five new periodic orbits keep the same.As a result, under linearization, the disturbance of the semi-major axis leads to a divergent solution,and each left bounded relative motion solution is periodic.This is the reason why the bounded solutions of HCW equations and TH equations are both periodic.Different with the two-body problem,the families of periodic orbits in CRTBP are one-dimensional and nearby orbits usually have different periods.Therefore,for the linearized relative motion near a periodic orbit in CRTBP,the nearby periodic orbits lead to a divergent solution.And a slight phase shift results in a periodic solution.Notice that the solution set of the linearized relative motion contains six independent solutions.The other four solutions arise from the manifolds of the periodic orbit.Stable, unstable, and center manifolds will lead to convergent, divergent, quasi-periodic solutions, respectively.35For DROs, the left four solutions are both quasi-periodic.Each quasi-periodic solution is composed of more than one periodic component, does not repeat itself, and consequently has a much more complicated geometry than that of periodic solutions.

To handle these issues, we reconstruct the general solution based on the Floquet theory as different modes and expand the non-analytical parts in each mode as truncated Fourier series.The specific works are shown below.The linearized model of relative motion on DRO is established at first.Based on the Floquet theory,we obtained a set of six independent solutions.By removing the divergent solution, the bounded solution set with five independent solutions is obtained.In the five bounded solutions, the quasi-periodic solutions with the same periodic components are unified as a quasi-periodic mode through proper parameterization.In this way, the general bounded solution is expressed as the linear combinations of a periodic mode and two quasi-periodic modes.And the non-analytical parts in each mode are expanded with truncated Fourier series.Consequently, the analytical bounded solution of linearized relative motion on DRO is derived.Based on the analytical solutions,the geometry of each mode is comprehensively analyzed, including the geometry shape of the periodic mode and bounded areas of the quasi-periodic modes.It is found that periodic mode can be the parking orbits of rendezvous and docking.Finally,a fly-around formation on DRO is designed to show the capability of the obtained general bounded solution.And the fly-around formation is transformed into a more realistic dynamics to obtain practical results.

This paper is organized as follows.Section 2 describes the establishment of the linearized relative dynamics on DRO.Section 3 demonstrates the general solution based on Floquet theory and the analytical bounded solution with different modes.Section 4 introduces the design of a fly-around formation based on the analytical bounded solution.And the flyaround formation is validated in a more practical ephemeris model in Section 5.

2.Dynamics

In this section, the absolute and relative motion dynamics are introduced successively.As the classical absolute dynamics of the three-body problem, the CRTBP is presented at first.And to describe the close relative motion on DRO, the linearized relative motion dynamics is established.

2.1.CRTBP

CRTBP is a dynamics model to describe the motion of a massless spacecraft attracted by two massive bodies.Two massive bodies orbit each other following circular orbits.In the Earth-Moon system, two massive bodies are the Earth and the Moon,respectively.Rotating frames can simplify the analysis of motion in CRTBP.Since DROs orbit the Moon, a Moon-centered rotating frame M is adopted and shown in Fig.1: the origin OMis located in the center of the Moon;xM-axis points from the Earth to the Moon; zM-axis is along with the angular momentum of the Earth-Moon system; yMaxis follows the right-handed orthogonal system.As also shown in Fig.1, the inertial frame I is defined nearly in the same way as the rotating frameM, except that it is based on the state of the Earth-Moon system in a specific time and thus has fixed axis directions.

To simplify the calculation and analysis, all variables are expressed with normalized units.The normalized mass, position, and time units are shown as follows:

Fig.1 Moon-centered rotating frame M and inertial frame I.

2.2.Linear relative dynamics

Consider two spacecraft, including a chief spacecraft and a deputy spacecraft.The chief spacecraft is located on a DRO,and the deputy spacecraft is located nearby.The motion of the deputy spacecraft relative to the chief spacecraft is described in the Local Vertical Local Horizontal(LVLH)frameL:the origin OLlies in the center of the chief spacecraft;xL-axis follows the position vector of the chief spacecraft in the frameM; zL-axis is along the orbit normal; yM-axis completes the right-handed orthogonal system.

Based on Eq.(2), the absolute motion of two spacecraft in the rotating frame M are shown as follows:

where the subscripts‘‘c”and‘‘d”denote the chief and deputy spacecraft, respectively, and the notations‘‘[˙?]i,i=(M,I,L)”and‘‘[¨?]i,i=(M,I,L)”stand for the first/second-order derivative operations in framei, respectively.

Table 1 Parameters of Earth-Moon system.

It should be noticed that the vectors in Eqs.(4) and(8) are expressed in frame M and frame L, respectively.Thus the transformation matrices between these two frames are given here.

The relative motion on DRO can be analyzed using the linear dynamics.

3.Relative motion solution

The linearized relative dynamics in Eq.(15) has a periodic coefficient matrix.Therefore, the general solution is obtained using the Floquet theory.The bounded terms in the general solution are further classified as different modes.

3.1.General solution based on Floquet theory

where P( t )S=Φ(t,0 )Se-Jtis the so-called Floquet modal matrix and P(0 ) is the identity matrix.According to the Floquet theory,H=eJTis the Jordan canonical form of the monodromy matrix Φ(T,0),where J is the Jordan matrix and T is the period of the reference DRO.37And S is a nonsingular constant matrix satisfying

Since S is nonsingular and constant, Φ(t,0 )S=P(t )SeJtis also a fundamental matrix solution of Eq.(15).Based on the Floquet theory, it can lead to a semi-analytical general solution37.

To analyze the solution of Eq.(15), a reference DRO is chosen with a period of 13.64 days,which is 2:1 resonant with the Moon’s orbital period.The 2:1 DRO is shown in Fig.2.As can be seen,the initial position of the chief spacecraft is chosen at the perigee.

The monodromy matrix of the 2:1 DRO has two pairs of conjugated eigenvalues with unity magnitude and a pair of real unity eigenvalues.

Fig.2 2:1 reference DRO.

where p1,p3and p5are the corresponding eigenvectors of λ1,λ3,and λ4,respectively,and p4is the generalized eigenvector ofλ3.Besides,operators Re(?)and Im(?)take the real and imaginary parts of the vector, respectively.To simplify the following analysis,the eigenvectors are normalized as shown in Eq.(21).

Since H=eJT, there is

According to Eqs.(21) and(23), the Floquet modal matrix can be derived as

Since^ek(t ) (j=1,2,...,6)are periodic with periodT,x3t( )is periodic, x4(t ) is divergent, and xj(t ) (j=1,2,5,6) are quasi-periodic.As analyzed in Section 1, these solutions originate from different nonlinear dynamical structures in the vicinity of the reference DRO.35The periodic solution x3(t )corresponds to the case that the deputy spacecraft is located on the reference DRO with a phase shift.As for the divergent solutionx4(t ), the deputy spacecraft lies in another nearby DRO with a different period.And the quasi-periodic solutions originate from the center manifold of the reference DRO.In the frameL, the specific values of ^ej(0 ) (j=1,2,...,6) are shown as follows:

As can be seen,^ej(0 ) (j=1,2,3,4)are in the orbital plane,and ^ej(0 ) (j=5,6) are in the normal direction.Analysis of A(t ) shows that the in-plane and out-of-plane relative motion are independent of each other.Therefore, xj(t ) (j=1,2,3,4)are in-plane solutions and xj(t ) (j=5,6) are out-of-plane solutions.

In fact, monodromy matrices of stable DROs have the same eigenvalue structure, of which all lie in the unit circle,two pairs are conjugated complex and one pair is real.Only the corresponding eigenvectors of one pair of conjugated complex eigenvalues are in the normal direction, while the other eigenvectors are on the orbital plane.Thus,the relative motion on each stable DRO has the same form as Eq.(27).

3.2.Transformation to elliptic-orbit SDR coordinates

Different solutions of the relative motion are given in Eq.(27)based on the Floquet theory.However,only the general forms of the solution are given instead of specific characteristics.In this subsection, the bounded solution is analyzed in detail.In accordance with Eq.(27), the bounded solution is composed of three different modes

where the subscripts‘‘pe”,‘‘p”, and‘‘n”denote the planar periodic mode, the planar quasi-periodic mode, and the normal quasi-periodic mode, respectively, and kj(j=0,1,2) are the amplitudes of three modes.Three different modes are shown in

where c21+c22=1 andc25+c26=1.Here, cj(j=1,2,5,6) are designed to combine planar quasi-periodic solutions and normal quasi-periodic solutions into two uniform expressions in the following subsections, respectively.According to Eq.(27), the non-analytical parts in three modes, that is^ej(t ) (j=1,2,3,5,6), are all T-periodic.The truncated Fourier series is a powerful tool to approximate periodic functions with arbitrary accuracies.Therefore, the analytical form can be obtained by expanding^ej(t )into truncated Fourier series.In this way, the geometry characteristics of three different modes are analyzed successively in the following subsections.

3.2.1.Planar periodic mode

The planar periodic mode has the potential to support formation missions due to the periodic property.To apply the Fourier series expansion, an angular variable θ0=2πt/T is introduced as the independent variable instead of the time.Thus, the planar periodic mode becomes

where ajand bjcan be obtained through discrete Fourier transform,HOT represents the Higher Order Terms.In general,the nth truncated order requires 2n+1 samples to the orbit.And higher truncated order leads to smaller truncated error.When the 25-th truncated order is applied, the relative error is below10-9.In this way, the absolute error is sub-millimeter for a formation with a scale of hundreds of kilometers, which is enough for close relative motion analysis and formation design.Therefore,only dozens of samples are needed to obtain an accurate analytical form of the planar periodic mode.In fact,the truncated order should be set based on the nonlinearity of the reference orbit.Higher nonlinearity requires larger truncated orders to achieve the same relative error.And for most planar DROs,a 50th truncated order can lead to the relative error below 10-9.

The analytical form with truncated Fourier series gives us the freedom to compute the relative states of arbitrary time points without numerical integration and dramatically reduces the computational cost.Thus, the analytical form facilitates the analysis of the periodic mode.The specific geometry of the planar periodic mode is shown in Fig.3.As can be seen,the planar periodic mode is a closed curve with an intersection,and thus can be divided into an inside curve and an outside curve.The intersection point on the periodic mode corresponds to two points on DRO.These two points separate the reference DRO into two parts, of which the left one and the right one correspond to the inside curve and the outside curve, respectively.Analysis shows that passing through the inside curve and outside curve takes 7.09 days and 6.55 days,respectively.Since the outside curve is longer with shorter passing time,the deputy spacecraft moves faster in the outside curve.The intersection point not only helps to understand the geometry of the planar periodic mode,but also may be the key point of formation missions.It divides the planar periodic mode into two parts, and meanwhile it is the nearest position to the chief spacecraft.Considering the distance,positions near the intersection point are better choices to observe or communicate with the chief spacecraft.

Fig.3 Planar periodic mode.

The planar periodic modes with different amplitudes are illustrated in Fig.4.It is shown that the planar periodic mode always lies near the y-axis of the frame L.By choosing different amplitudes, the periodic mode has different scales but the same view angle as seen from the chief spacecraft.And this characteristic is reasonable since the planar relative motion corresponds to the formations that the deputy spacecraft is located on the reference DRO with small phase shifts.It makes the planar periodic mode appropriate to be the parking orbit of rendezvous and docking missions on DRO.Since the intersection point lies on the y-axis of the frame L, it has potential to be the starting point of the final approaching phase for docking to the chief spacecraft from the positive or negative y-axis.

3.2.2.Planar quasi-periodic mode

Fig.4 Planar periodic mode with different amplitudes.

Unlike the periodic mode,the quasi-periodic modes have multiple periodic components and thus are much more complicated.In Eq.(27), both x1(t ) and x2(t ) are composed of two periodic components.^e1(t ) and ^e2(t ) have a period of T, and cos (α1t/T) and sin (α1t/T) have a period of 2πT/α1.Combining x1(t ) and x2(t ), the planar quasi-periodic mode can be obtained as

Fig.6 Invariant curves of θ1 of planar quasi-periodic mode.

Fig.5 Invariant curves of θ0 of planar quasi-periodic mode.

Fig.7 Schematic process to obtain bounded area of planar quasi-periodic mode.

Fig.8 Bounded area and evolution of planar quasi-periodic mode with θ1,0 = π/3.

Fig.6 illustrates the invariant curves of θ1with different constant θ0.Different from the invariant curves of θ0, each invariant curve of θ1is approximately elliptical with similar scales and orbits the chief spacecraft anticlockwise.With the increase of θ0, invariant curves of θ1rotate clockwise.there needs dozens of invariant curves and hundreds of sample points on each invariant curve to obtain the bounded area with enough high accuracy, which is colored with gray in Fig.8.And Fig.8 also demonstrates the evolution of the planar quasi-periodic mode during 7T.It can be seen that the trajectory is always inside the bounded area.Under the joint influence of two periodic components, the planar quasi-periodic mode orbits the chief spacecraft and oscillates between the inside and outside boundaries simultaneously, which leads to

3.2.3.Normal quasi-periodic mode

Fig.9 Planar quasi-periodic mode during a DRO period with different θ1,0.

Fig.10 Invariant curves of θ0 and θ2 of normal quasi-periodic mode.

Fig.11 Evolution of normal quasi-periodic mode with θ2,0=π/6.

The analytical expression of xn(θ0,θ2) can be obtained by expanding ^e5and ^e6as 25-th order truncated Fourier series.Same as the analysis of planar quasi-periodic mode,the analytical form is also essential to calculate invariant curves of the normal quasi-periodic mode.

The normal quasi-periodic mode only has zL-component.The amplitudes of invariant curves of θ0and θ2are illustrated in Fig.10.It is shown that the invariant curves of θ0are symmetrical with zL=0,since they are simply sinusoidal functions according to Eq.(35).As illustrated in Fig.10, the lower and upper bounds of the normal quasi-periodic mode are -0.665 and 0.665, respectively.The evolution of the normal quasiperiodic mode during 20T is shown in Fig.11.The values of zLare always between the lower and upper bounds and have a relatively stable fluctuation.

3.2.4.Reparameterized bounded solution

Combining Eqs.(31), (33), and (35), Eq.(29) can be transformed to a new reparameterized form as follows:

which contains five design parameters, including three amplitudes kj(j=0,1,2) and two initial angles θ1,0andθ2,0.

Eq.(37)is a complete analytical form of the bounded solution of the linearized relative motion on a DRO.It replaces and surpasses the function of complex numerical integration and gives us a new perspective to understand two quasiperiodic modes.By choosing the five design parameters, we can obtain arbitrary bounded relative trajectories to satisfy the requirements of different formation missions.

The proposed procedure can simplify the design of close relative trajectories on DRO.In the practical design, some alterations may be needed on the basis of Table 2.The firstone is that some steps may be ignored if the corresponding modes do not benefit the desired formation.And some parameters may be designed at the same time due to the coupled influence on the formation.In addition,impulses may be introduced to this procedure to complete, adjust, or maintain the relative geometries.To show the capability of the reparameterized bounded solution in Eq.(37), a fly-around formation is designed in the next section following this procedure.

Table 2 A general procedure to determine design parameters.

4.Formation design

4.1.Planar formation

Fig.12 Minimum distance between two spacecraft and distance between initial and final positions of natural trajectories based on planar quasi-periodic mode (k1 = 1).

Fig.13 Best natural trajectories based on planar quasi-periodic mode (k1 = 1).

Fig.14 min‖ρ(t )‖ and ‖ρ(0 )-ρ(T )‖ of natural trajectories based on planar periodic mode and planar quasi-periodic mode(k1 = 1 and θ1,0 = 0.625π).

where Δv1and Δv2are the maneuver impulses in the initial and final positions of the transfer trajectory, respectively.For the planar fly-around mission,ρ and ρvcan be calculated through Eq.(37) by setting k0/k1=-0.197k2=0 and θ1,0=0.625π-nα1.And Eq.(40) can be solved with root-finding algorithms.Although different parameters are applied to calculate the natural trajectory in each lap, the specific trajectory keeps repeated.Therefore,the two impulses in each lap remain the same.The total impulse Δv=‖Δv1‖+‖Δv2‖in each lap with respect to different Δt are demonstrated in Fig.16,where k1is set as 8.624×10-6to obtain a 1 km scale formation.The impulse decreases rapidly at first with the increase of Δt.When Δt is larger than 0.081T,the total impulse is below 1 cm/s.And the minimum impulse appears to be 1.57 mm/s when Δt is equal to 0.490T.And the fly-around trajectories with respect to different Δt are illustrated in Fig.17.It is reasonable that longer transfer time leads to a longer transfer trajectory.Besides,longer transfer time will shorten the minimum distance between two spacecraft.Therefore, considering the specific trajectory and fuel consumption, transfer time between 0.1T and 0.3T is practical choices.In these cases, the complete planar flyaround trajectories are quasi-elliptical with centers near the origin.As for the case that θ1,0=1.625π and k0=0.197, centrosymmetric trajectories can be obtained due to the linear property of the dynamics.

In conclusion,the planar fly-around formations are accomplished to have quasi-elliptical trajectories with two impulses per DRO period.And for a 1-km scale trajectory,the fuel consumption per period is less than 1 cm/s.

4.2.3D formation

Fig.15 Evolution of natural trajectory with the decrease of k0 (k1 = 1 and θ1,0 = 0.625π).

Fig.16 Total impulses per period of planar fly-around formation (k1 = 8.624 × 10-6, k0/ k1 = -0.197, k2 = 0, and θ1,0 = 0.625π－nα1).

On the basis of planar fly-around trajectories, the normal quasi-periodic mode is added to design 3D fly-around trajectories.And the normal quasi-periodic mode is analyzed separately because the three modes in Eq.(37) are independent of each other.Similar to the planar quasi-periodic mode, the normal quasi-periodic mode has two design parameters,including the amplitude k2and the initial angle θ2,0.As the amplitude k2does not affect the geometry, it is set as 1 and only the initial angle θ2,0is analyzed.To preserve the geometry characteristics of planar fly-around trajectories, the distance between the initial and final positions of the normal quasiperiodic mode is desired to be minimum, that is,

For the 3D transfer trajectory, the impulses can be calculated based on Eq.(40) as well, where the parameters in ρ and ρvare set as k0/k1=-0.197k2/k1=0.5,θ1,0=0.625π-nα1and θ2,0=0.266π-nα2.Here, k2is chosen to be 0.5k1to balance the amplitudes in planar and normal components.In fact,different values are free to be chosen for k2/k1to satisfy varieties of mission requirements.And k1is still set as 8.624×10-6to obtain a 1-km scale formation.The total impulses are shown in Fig.20.Since the relative motion inside and outside the orbital plane are independent of each other,the planar impulse components calculated in the last subsection are still feasible,and only the impulses in the zL-axis need to be added.Therefore, for the same transfer time, the total impulse of the 3D fly-around trajectory is definitely larger than that of the planar fly-around trajectory.Specifically, with the increase of transfer time, the impulse in the zL-axis pulls up fuel consumption in Fig.20 compared with that in Fig.16.And a singularity exists when Δt=0.565T.It is shown in Fig.21 that the directions of transfer trajectories overturn at the singularity.And the amplitudes in zL-axis are extremely high near the singularity,and thus the corresponding impulses are extremely large.When Δt<0.565T, the amplitude in the zL-axis increases with the rise of the transfer time, leading to an increasing impulse in zL-axis.Therefore, considering the total impulses and the specific geometries, the transfer time between 0.1T and 0.2T corresponds to practical 3D flyaround trajectories.In these cases, the total impulses are less than 1 cm/s and the quasi-elliptical shapes of the planar flyaround trajectories are preserved.Although the fluctuation in the zL-axis is non-symmetry with a larger amplitude in the negative direction, the safety of spacecraft can still be guaranteed by the planar components.And for the cases that θ2,0=1.266π-nα2, similar formations can be designed with trajectories symmetry with respect to the plane zL=0.

Consequently, the planar fly-around formations are extended to 3D cases with the introduction of normal quasiperiodic mode.And for a 1-km scale 3D fly-around trajectory,the fuel consumptions per period are still less than 1 cm/s.

4.3.Singularity of normal transfers

In the frameL,the transfer between two points may be singular for some specific initial and transfer time due to the normal relative motion,as can be seen in Fig.20.In the frame L,the linearized relative motions in the xy-plane and normal direction are fully decoupled.Therefore, only the normal relative motion needs to be considered to investigate the normal transfer singularity.

Fig.17 Planar fly-around trajectories (k1 = 8.624 × 10-6, k0/ k1 = -0.197, k2 = 0, and θ1,0 = 0.625π-nα1).

Fig.18 Distance between initial and final positions of normal quasi-periodic mode (k2 = 1).

The normal transfer problem is actually to find the initial velocity ˙zL(t0) given the initial time t0, the transfer time Δt,the initial normal position zL(t0), and the final normal positions zL(tf).Here, tf=t0+Δt.Using the state transition matrix to propagate the linearized relative motion,the transfer problem is expressed as39

Here, Φz(tf,t0) is the state transition matrix of the normal component, and is actually a submatrix of Φ(tf,t0), which is shown as

Fig.19 Trajectories with minimum distance between initial and final positions of normal quasi-periodic mode (k2 = 1).

Fig.20 Total impulses per period of 3D fly-around formation.

According to Eq.(44), if and only if φ36is equal to 0 and the right-hand side term of Eq.(44) is not equal to 0, there does not exist ˙zL(t0) to satisfy Eq.(44), that is, the normal transfer problem is singular.When φ36is equal to 0, Eq.(44)becomes

which means that the initial and final normal positions are constrained to have a constant ratio φ33(tf,t0).Therefore, for most initial and final positions, the normal transfer problem is singular when φ36is equal to 0.

Fig.22(a)shows the variation of φ36over initial time t0and transfer time Δt.The zero points of φ36exist in the crossing lines of the curved surface φ36(tf,t0) and the plane φ36=0,and are illustrated in Fig.22(b).As can be seen, there exist multiple zero points for arbitraryt0.And the transfer time in zero points varies with the initial time.The fluctuation of φ36over Δt makes the zero points keep emerging.On average,approximately-one additional zero point appears for every 0.65T increase in transfer time.For the normal transfers in the fly-around formations, there is a linear constraint between the initial time and the transfer time,that is t0+Δt/2=T.The linear constraint intersects with the zero-point curves of φ36at Δt=0.565T,which is consistent with the position of the singularity in Fig.22.

The value of φ33at zero points of φ36are shown in Fig.23.When Δt is equal to 0, zL(tf)=zL(t0).Therefore, φ33and φ36are constant 1 and 0, respectively.When Δt is not equal to 0,φ33and φ36vary over t0and Δt.

Now, the singularities of normal transfer problem can be obtained through Fig.22 and Fig.23.Since φ36varies continuously over t0and Δt,the absolute values of φ36are very small near the zero points, which may lead to large impulses as shown in Fig.20.To investigate this phenomenon, impulses around zero points are further analyzed.Near a zero point,the initial velocity ˙zL(t0) should satisfy

Fig.21 3D fly-around trajectories.

Fig.22 Variation and zero points of φ36.

Fig.23 Value of φ33 at zero points of φ36.

5.Fly-around formation in ephemeris model

In this section, the fly-around formation in CRTBP model is transformed into a more practical ephemeris model.The ephemeris of the Earth, the Moon, and the Sun are considered in the ephemeris model.The three celestial bodies are seen as point mass due to the far distance to spacecraft on DRO.In the Moon-centered J2000 frame, the dynamics of the spacecraft on DRO is

where μEμSand μMare the gravitational constant of the Earth,the Sun, and the Moon, respectively, and rEand rSare the position vectors of the Earth and the Sun, respectively.The positions of the three celestial bodies are computed from the planetary and lunar ephemerides DE430 of JPL.Two frame transformations are needed to obtain relative states in the frame L.The first one is to transform states of two spacecraft from the Moon-centered J2000 frame to the frame M.The second one is to transform the differential states between two spacecraft into the frame L as shown in

Fig.24 Value of ?φ36/?tf at zero points of φ36.

5.1.Design of fly-around formations

In CRTBP, there exist fixed mappings between the periodic DRO and analytical bounded relative motion, since they both have a T-periodic component.Given the natural trajectory and the transfer time,the absolute and relative maneuver positions can be easily obtained.However, in the ephemeris model, the reference orbit is no longer periodic, and the relative motion is non-analytical with slightly fluctuated dynamics.Therefore,the natural trajectory and the transfer trajectory of the flyaround formation in CRTBP will both become transfer trajectories in the ephemeris model.And the fly-around trajectories in the ephemeris model should be obtained by successively solving the transfer problem.As a result, the key to design the fly-around formation is to find proper maneuver positions in the ephemeris model with respect to those positions in the CRTBP.

To find reasonable absolute maneuver positions,there lacks a bond to connect the reference orbits in two dynamical models.Considering that the reference orbits in two dynamical models both orbit the Moon with similar geometries, a new geometric phase f in the frame M is defined as the angle of clockwise rotation from the positive x-axis to the position vector projected on the xy-plane as shown in

For fly-around formations in CRTBP, the two absolute maneuver positions rm1/rm2and their geometric phases fm1/fm2can be easily obtained given the transfer timeΔt.Assume that the absolute maneuver positions in the ephemeris model have the same geometric phases fm1/fm2, and the corresponding absolute maneuver positions r1,r2,r3,...,rnand maneuver time t1,t2,t3,...,tncan be computed, as seen in Fig.25(a).In the ephemeris model, the reference orbit successively passes through fm1and fm2in each circle.Therefore, the absolute maneuver positions will be alternatively fm1and fm2.The specific order is decided by the initial state of the reference orbit.

To make the fly-around trajectories in two dynamical models close enough, the relative maneuver positions are set as anchor points,and thus keep the same in two dynamical models.The two relative maneuver positions are shown in Fig.25(b), and are denoted as ρm1and ρm2, respectively.

With the absolute and relative maneuver positions and the maneuver time, the natural trajectory and the transfer trajectory of the fly-around formation in CRTBP can be transformed as the transfer trajectories from ρm1to ρm2and from ρm2to ρm1in the ephemeris model,respectively.In the transfer problem from the maneuver time tkto the maneuver time tk+1,the impulse Δvkat the time tkshould satisfy

With the above designs, the fly-around formation in the ephemeris model can be obtained.The planar and 3D flyaround formations are successively designed in the following two subsections.

5.2.Planar formation

To design the fly-around formation in the ephemeris model, a reference orbit is needed.Although periodic DROs do not exist in the ephemeris model, there exist orbits close to DROs in CRTBP while keeping stable for a long time due to the favorable stability of DRO.A reference orbit close to the 2:1 DRO in CRTBP is obtained through the optimization method in Ref.40The initial epoch and optimized initial state are shown in Table 3.The evolution of the optimized results during 150 days is shown in the left subfigure of Fig.25.As can be seen, the reference orbit keeps close to the 2:1 DRO in CRTBP.

Based on the analysis in Section 4.1,the minimum distance between two spacecraft achieves maximum for the planar flyaround formation when the coefficients are set as k0/k1=-0.197k2=0 andθ1,0=0.625π-nα1.For linearized relative motion, the relative trajectories are strictly linear with the scale.For this reason, only the 1-km scale fly-around formation is shown in Section 4.1.However,the strict linear relationship does not exist in the ephemeris model.Therefore, to show the feasibility of the fly-around formation in the ephemeris model, different scales from 1 km to 1000 km are analyzed.Besides, as the reference orbit is not periodic in the ephemeris model,the fuel consumptions and fly-around trajectories in different circles do not keep the same anymore.To compare formations in two dynamical models, ten circles of fly-around trajectories during around 138 days are considered,and the average,maximum,and minimum fuel consumption in all circles are computed.Fig.26 shows the total impulses per circle of the planar fly-around formation in the ephemeris model with different transfer time and scales (Eph represents ephemeris).The upper and lower bounds of each yellow area are the maximum and minimum impulses per circle, respectively.The horizontal axis of Fig.26 is still the transfer time Δt in CRTBP,as the actual transfer time fluctuates in different circles in the ephemeris model.In the ephemeris model,the two types of transfers,that is,from ρm2to ρm1and from ρm1to ρm2,have transfer time around Δt and T-Δt, respectively.In Fig.26, impulses are almost linear with the scale.And it can be seen in the local zoomed subfigure that the average impulses are nearly the same with the impulses in CRTBP,and the fluctuation of impulses are less than 10% in most cases.

Fig.25 Absolute and relative maneuver positions.

Table 3 Epoch time and initial state of reference DRO in frame M and in ephemeris model.

Fig.26 Total impulses per circle of planar fly-around formation in ephemeris model.

The planar fly-around trajectories with different transfer time and the same scale are illustrated in Fig.27.As can be seen, the trajectories in the ephemeris model are always close to the trajectories in CRTBP.And the planar fly-around trajectories with the same transfer time and different scales are shown in Fig.28.The trajectories with different scales nearly keep the same geometry.It means that the fly-around formations in the ephemeris model remain the favorable linear features, which is not strictly but approximate enough to simplify the analysis of fuel consumptions and trajectories.

Fig.27 Planar fly-around trajectories with different transfer time the same scale in ephemeris model.

Fig.28 Planar fly-around trajectories with same transfer time and different scales in ephemeris model.

Fig.29 Two planar fly-around trajectories in spikes of time-impulse curve in ephemeris model.

5.3.3D formation

Fig.30 Total impulses per circle of 3D fly-around formation in ephemeris model.

Following the analyzed results in Section 4.2, the coefficients of the 3D fly-around formation are chosen as k0/k1=-0.197k2=0.5θ1,0=0.625π-nα1and θ2,0=0.266π-nα2, so that the two maneuver positions have the same normal components.The total impulses per circle with different scales and transfer time are shown in Fig.30.Same as planar formations, the impulses are approximately linear with the scale.And in most cases, the average impulses in the ephemeris model are quite close to the results in CRTBP.The fluctuation of impulses are below 10%in most cases when Δt is far away from the two singularities.

However, the difference with planar formations is that the singularities have much larger impacts on the 3D formations as shown in Fig.30.Compared with the impulses nearΔt=0.306T,the impulses near Δt=0.565T fluctuate with larger amplitude and larger transfer time span.In CRTBP,the transfer problem from ρm1to ρm2satisfies Eq.(45),and thus is not singular.Therefore, the reason for the impulse fluctuation near Δt=0.306T is the limited zMcomponent of the reference orbit in the ephemeris model.However, the transfer problem from ρm2to ρm1is already singular in CRTBP, which leads to more dramatic fluctuations in the ephemeris model.The fly-around trajectories in two fluctuation segments are shown in the first row and second row in Fig.31,respectively.Overall,the closer it is to the singularity,the more sensitive the normal transfer problem is, and the more dramatic the effect of the dynamical differences between different circles of the reference orbit on the normal transfer is.Since Δt=0.565T is the singularity of the transfer from ρm2to ρm1in CRTBP, the normal components have opposite directions before and after the singularity based on Eq.(50), as shown in the subfigures with Δt equal to 0.53 T and 0.595 T of Fig.31.When close to a singularity enough,the different circles of the reference orbit lead to normal components with different directions due to the high sensitivity of the normal transfer to the dynamical differences,as illustrated in the subfigures of Fig.31 with Δt equal to 0.56 T and 0.31 T.Because the transfer problem from ρm1to ρm2is not singular when Δt=0.306T in CRTBP, the normal components before and after 0.306 T are still with the same direction, as shown in the subfigures of Fig.31 with Δt equal to 0.3 T and 0.32 T.Besides,it is found that the transfer problem may not converge near the singularities due to the high sensitivity.

Fig.31 3D fly-around trajectories in two fluctuation segments of time-impulse curve in ephemeris model.

Fig.32 3D fly-around trajectories with different transfer time and the same scale in ephemeris model.

Fig.33 3D fly-around trajectories with the same transfer time and different scales in ephemeris model.

The 3D fly-around trajectories far from singularities are illustrated in Fig.32 and Fig.33.Fig.32 shows 3D flyaround trajectories with different transfer time and the same scale in the ephemeris model.As can be seen,when the transfer time is far from singularities, the trajectories in the ephemeris model are always located near the trajectories in CRTBP.And Fig.33 demonstrates 3D fly-around trajectories with the same transfer time and different scales in the ephemeris model.Same as the planar formations, the 3D fly-around trajectories with different scales keep the favorable linear property.

In this section,the planar and 3D fly-around formations are transformed into a more realistic ephemeris model.The trajectories and impulses are consistent in two dynamical models,except the cases near singularities.And the results near singularities in the ephemeris model are explained through the singularity analysis in CRTBP and the difference between reference orbits in two dynamical models.

6.Conclusions

This paper presents a novel approach to construct the analytical bounded relative motion on DROs.The bounded relative motion is composed of three modes with different periodic components through a new parameterization.The parameterization allows us to thoroughly explore the geometries of quasi-periodic modes through the invariant curves.To construct the analytical form, the periodic components in each mode are expanded as truncated Fourier series.With the analytical form, each mode is analyzed in detail.The planar periodic mode lies near the y-axis of the frame L and has the same view angle as seen from the chief spacecraft.These characteristics make the planar periodic mode appropriate to be the parking orbits of the rendezvous and docking missions.The planar quasi-periodic mode is within a bounded hallow area and each invariant curve of θ0orbits the chief spacecraft clockwise.The normal quasi-periodic mode has relatively stable amplitudes and hence can provide normal component for 3D formations.

Based on the three analytical modes, a fly-around formation is further designed based on the analytical bounded relative motion and short transfer trajectories.The fly-around formations support the observations to targets on DROs from different view angles with very few fuel consumptions.Both planar and 3D fly-around formations are designed.Normal transfer singularities are found in 3D fly-around formations.Analysis shows that singularities of normal transfers generally exist for arbitrary initial time and keep coming up with the increase of the transfer time.And the directions of the normal transfer trajectories are opposite before and after each singularity.Moreover, the fly-around formations are transformed into the ephemeris model.The trajectories and impulses in two dynamical models are consistent with each other.Singularities also exist in the ephemeris model, and affect not only the normal transfers but also the planar transfers,and yet they are still located near the singularities in CRTBP.Therefore,the singularity analysis in CRTBP provides essential and quantitative basis for singularity avoidance in practical formation designs to prevent dramatically increasing impulses, nonconvergent transfers, and highly uncertain transfer trajectories.The design results demonstrate that the basic dynamical features of the linearized relative motion near DRO in CRTBP still approximately exist in the ephemeris model, including the favorable linear property.As a result, the linearized relative dynamics in CRTBP is capable of supporting preliminary analysis for close relative motion on DRO.

In conclusion, the proposed approach significantly boosts and extends the relative motion investigation and formation design on DROs.And this approach can also be applied to other non-analytical periodic orbits with center manifolds,such as halo orbits, tadpole orbits, and orbits surrounding asteroids.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgement

This study was supported by the Strategic Priority Research Program of Chinese Academy of Sciences (No.XDA30010200).