Evolution of space debris for spacecraft in the Sun-synchronous orbit

Abstract In this paper, the evolution of space debris for spacecraft in the Sun-Synchronous orbit has been investigated. The impact motion, the evolution of debris from the Sun-Synchronous orbit, as well as the evolution of debris clouds from the quasi-Sun-Synchronous orbit have been studied. The formulas to calculate the evolution of debris objects have been derived. The relative relationships of the velocity error and the rate of change of the right ascension of the ascending node have been presented. Three debris objects with different orbital parameters have been selected to investigate the evolution of space debris caused by the Sun-Synchronous orbit. The debris objects may stay in quasi-Sun-Synchronous orbits or non-Sun-Synchronous orbits, which depend on the initial velocity errors of these objects.


Introduction
The increasing of population of the space debris poses a hazard to human spacecraft (Liou and Johnson 2006;ESA Space Debris O ce 2017). Previous literature investigates the orbit determination, collision probability, and removal of the space debris, as well as the impact and evolution of the debris cloud. About the observation and orbit determination of space debris, Ananthasayanam et al. (2006) applied the equivalent fragments to model the space debris scenario, and used the constant Kalman gain method to handle the errors caused by the nite bin size and the perturbations. Fujita et al. (2016) presented a computational method to calculate the plane of trajectory orbit of a micronsized debris cloud, and used the simulation test to demonstrate the e ectiveness of the method. Vallduriola et al. (2018) investigated the Optical IN-SITU Monitor project to detect space debris; several algorithms of data reduction are used for the image ltering.
Space debris may collide with spacecraft; spacecraft may also collide with each other. Kamm and Willemson (2015) studied the secure multiparty computation (SMC) to calculate the collision probability between two satellites. Kurihara et al. (2015) investigated the impact frequency calculated by the debris environment model and the impact craters on the surface of the exposed instrument of the Tanpopo mission; they suggested that the impact energy is in proportion to the crater volume, and is not related to the projectile materials and relative speed. Letizia et al. (2016) modeled the debris clouds as a uid with the density changed under the force caused by the atmospheric drag, they derived the analytical expression for the density evolution of the uid, and the method can be applied for the simulation of several collision scenarios in a short time. Francesconi et al. (2015) used the honeycomb sandwich panels to model the structure of the spacecraft and presented an engineering model describing debris clouds. The debris clouds include larger fragments and several small fragmented particles. Dutt and Anilkumar (2017) investigated the orbit prediction of space debris with the semianalytic theory; the results can be applied to analyze the lifetime estimation of debris. Zhang et al. (2020) studied the characteristics structure of the debris cloud generated by the impact of disk projectile on a thin plate.
To reduce the risk of space debris, previous studies investigate the debris removal. Bombardelli and Peláez (2011) assessed the feasibility of the ion beam shepherd spacecraft (IBS) concept for contactless maneuvering of space debris without considering the attitude motion of the debris. Aslanov and Yudintsev (2013) discussed the removal of large debris by the space tug to the surface of the Earth. The attitude motion of the debris and the gravita-tional torque acting on the debris is considered. Rubenchik et al. (2014) investigated the space-debris cleaning with the ground-based laser system, they found the self-focusing of the powerful laser pulses in the atmosphere can obviously decrease the laser intensity on debris. Botta et al. (2016) applied the Vortex Dynamics to the simulation for the capture of space debris using the tether-nets. Qi et al. (2017) investigated the dynamics and control of a double-tethered space-tug system which is used for the debris removal.
To study the characteristic of the evolution of space debris is useful for the design of the trajectories for the debris removal. However, the characteristics of the evolution of space debris in di erent kinds of orbits are di erent. In this paper, we choose the Sun-Synchronous orbit to investigate the evolution of space debris. The Sun-Synchronous orbit is the most common orbits in the Low Earth Orbit . The evolution of space debris is related to the relative motion between space targets (Cao and Chen 2016;Dang and Zhang 2018). The evolution of the inclination and the right ascension of the ascending node is under the perturbation including the nonspherical perturbation of the Earth, the third body perturbation, the drag perturbation, and the solar radiation pressure. The impact of a debris object and a spacecraft which are in the sun-synchronous orbit has been calculated by the material point method. We have chosen three debris objects to calculate the evolution of debris. The Earth gravitation with 10×10 EGM96, the Solar gravitation, Lunar gravitation, the atmospheric model with NRLMSISE00 and USSA1976, as well as solar radiation pressure are considered when calculating the orbital evolution of debris. The relative position and relative distance between di erent debris objects have been presented. In addition, the relative relationship of the objects' velocity errors and the rate of change of the right ascension of the ascending node has been calculated, the results showed the deviation of the orbital characteristic of the sun-synchronous. Relative trajectories between di erent debris objects may be quite di erent, which depends on the velocity error.

Motion Equation
For the sun-synchronous orbit, the rate of node precession and the Earth's mean motion relative to the sun are the same. The mean nodal precession rate  satis esΩ Where Ω represents the right ascension of the ascending node, a represents the semi-major axis, e represents the orbital eccentricity, i represents the inclination, J 2 represents the second-order zonal harmonic, n represents the mean angular velocity, Re represents the reference radius of Earth, ns represents the Earth's mean motion around the sun.
Using the above equation, one can get Here δ represents the di erence operator. From the above equation, one can see that the variety of semi-major axis, eccentricity, and inclination cause the variety of the mean nodal precession rate. Using the following equation One can get δφ = 3 2 δa a n Where φ = ω+M, ω is the argument of perigee, and M is the mean anomaly. The above equation implies that the variety of semi-major axis cause the variety of phase relative to the perigee (Yang et al. 2016;Ren and Shan 2018). Considering the semi-major axis and the inclination have perturbations with the change of time, let Thus, we have In this height, the atmospheric drag perturbation can be neglected, and we thus only consider the third-Body perturbation e ects to the inclination. Thus the above equation

Impact Motion
In this section, we use the material point method (MPM) to calculate the impact motion of the debris and a spacecraft (Gong et al. 2012;Lian et al. 2015). The material domain Ω can be discretized using a set of particles which is covered by the Eulerian background grid. The shape functions are de ned to link the mapping from the mass and momentum of particles to grid nodes.
Here m I and p iI represent the mass and the momentum of the particle on the I-th grid node, respectively. I represents the grid node, means the number of particles, m j as well as v ij represent the mass and velocity of the particle j, respectively. N Ij means the shape function at the location j. The nodal velocity for the background grid can be given by The momentum at the grid node can be calculated by the leap-frog algorithm by Here the superscript k ± 1 2 and k represent the time t k± 1 2 and t k , respectively. f k,ext iI and f k,int iI are the external and internal nodal force, respectively.
Here b k ip is the body force acting on the unit mass, and σ k ijp is the Cauchy stress acting on the particle p. The positions and velocities of the particle p can be calculated by Here Then the values at the grid nodes are reset to zero such that a new regular background grid can be used for the next integral step. We give a simple simulation of a debris object impact to a spacecraft. The numerical simulation here used 404224 material points. In the impact process, the debris object is neglected if the diameter is smaller than 1.0 cm. We use the step-size 30ms. Considering a spherical debris collide with a spacecraft which has the size of 1.0×0.8×0.5 m. The relative velocity is 3.0 km·s −1 . The impact direction is along the z axis, and the impact direction points through the centre of mass of the spacecraft. Figure 1 shows the impact process of the debris and the spacecraft.

Evolution of Debris from the Sun-Synchronous Orbit
The perturbation of the inclination (Jiang et al. 2013) of the debris on the Sun-synchronous orbit due to the Solar gravity is di dt = 3n 2 s 2n (cos βs cos Ω + sin βs cos is sin Ω) (cos βs sin i sin Ω − sin βs cos is sin i cos Ω + sin βs sin is cos i) , where βs represents the longitude of the ascending node of the apparent motion of the Sun, is represents the inclination of the ecliptic plane, ns represents the orbital angular velocity of the Earth relative to the Sun. i represents the inclination of the debris, and Ω represents the longitude of the ascending node of the debris. The orbital elements are expressed in the equatorial inertial system of the Earth. For the long-term motion, we don't consider the periodic term of βs, thus we have di dt = 3n 2 s 8n (sin 2Ω sin i + sin 2is sin Ω cos i − sin 2Ωcos 2 is sin i )︁ .
Considering that the debris is in the Sun-synchronous orbit, the above equation reduces to di dt = − 3n 2 s 16n sin i(1 + cos is) 2 sin (2βs − 2Ω) .
Considering the motion of the di erent debris after the impact, one can use the above equation and the following equation to calculate the relative evolution of orbital parameters for debris objects.

Evolution of Debris Clouds from the Sun-Synchronous Orbit
We consider the relative evolution of debris objects. The orbit coordinate system A-xyz is de ned as: Ax axis points to the reverse direction of the Earth centre to the object A, Az pints to the normal direction of the object A's orbital plane, Ay, Az, and Ax are right hand orthogonal. Let r be the distance between the object A and the Earth, f be the true anomaly, u = ω + f , λ = ω + M. The subscripts A and B represent the orbital elements of the objects A and B, respectively. In the Sun-Synchronous orbit, if the relative errors of the orbital elements are small after the impact, the distance between the two objects is much less than the distance between the object and the Earth centre. In the coordinate frame A-xyz, the relative position of B relative to A can be expressed as If at least one of relative errors of the orbital elements is large, or the relative distance of these two objects is large (Vadali 2002), the relative position of B relative to A can be calculated by We consider three debris objects in the debris clouds caused by the impact in the Sun-Synchronous orbit. The  initial orbital elements of these debris objects are listed in Table 1. The initial positions and velocities of debris objects can be given in the inertia space, and are presented in Ta-ble 2. The orbit prediction considered 10×10 EGM96 of Earth gravitation, Solar gravitation, Lunar gravitation, the atmospheric model with NRLMSISE00 and USSA1976, as well as the solar radiation pressure (El-Salama 2018; Frouard and Yokoyama 2013). The evolution of orbital elements of the object A has been presented in Figure 2, including the semi-major axis, the eccentricity, the inclination, the right ascension of the ascending node, and the Arg. periapsis. Figure 3 shows the relative position of object B relative to the object A in the coordinate frame A-xyz. Figure 4 shows the relative distance between these two objects A and B. We integrated for 86400 s. From Figure 4, one can see that the relative distance increases very fast during the evolution of orbital elements. This is because the error of the semi-major axis is large, which is 409255.34 m. The maximum error of the relative position in z axis is about 159.6 km, which is because the relative error of the inclination and the right ascension of the ascending node is small. The maximum error of the relative position in x and y axes are about 14881.5 km and 7865.1 km, respectively.   Figure 5 indicates the relative position of object C relative to object A in the coordinate frame A-xyz. Figure 6 shows the relative distance between these two objects A and C. The integral time is also 86400s. From Figure 5(a) and 5(b), one can see that the major evolution is along the y axis, i.e. the motion direction of the debris object. The projection of the trajectory of C relative to A in xz plane looks like an oblate ellipsoid. The maximum evolution of C relative to A is in the y axis, which is 26.856 after 86400s, which is quite small than the evolution of Brelative to A. At the initial time, the norm of the velocity error of objects A and B is 314.561 m·s −1 , while the norm of the velocity error of objects A and C is 0.384979 m·s −1 . Thus the evolution of the object C relative to the object A is slowly than the evolution of the object B relative to the object A. From Figure 6, one can see that the relative distance between objects C and A are not monotone increasing.
From the above numerical examples, one can conclude that the relative evolution of debris objects may be quite different, including the evolution and the relative trajectories. Evolution velocity in di erent coordinate axis depends on the initial relative velocity of the coordinate axis and the initial orbital parameters.
The initial orbit is in the Sun-Synchronous orbits. The objects A and C are in the quasi-Sun-Synchronous orbits. However, the velocity error between objects A and B are two large; this makes the object B orbits in a non-Sun-Synchronous orbit. Using the formulas presented in Section 2.2, one can calculate the relative relationship of the velocity error and the rate of change of orbital parameters. Figure 7 indicates the relative relationship of the velocity error and the rate of change of the right ascension of the ascending node. One can see that the velocity errors in a di erent coordinate axes can cause di erent variety of the rate of change of the right ascension of the ascending node. Figure 8 presents the 2D plot of the relative relationship of the velocity error and the rate of change of the right ascension of the ascending node. From Figure 8, one can see that the rate of change of the right ascension of the ascending node when the velocity errors in two coordinate axes coexist. The ranges of velocity errors are in [−2.0km/s, 2.0 km/s].
From the above results, one can have the following conclusion: In the debris clouds which is caused by the impact of a debris object and a spacecraft in the Sun-Synchronous  7. The relative relationship of the velocity error and the rate of change of the right ascension of the ascending node. The red point means that the rate of change of the right ascension of the ascending node equals 0.9856 deg/day when the velocity error is zero. (a) the velocity error is along the x axis of the inertial system; (a) the velocity error is along the y axis of the inertial system; (a) the velocity error is along the z axis of the inertial system. orbit, the debris objects can stay in quasi-Sun-Synchronous orbits or non-Sun-Synchronous orbits, which depends on the velocity error of debris objects.

Conclusions
This paper investigates the evolution of space debris caused by the impact between debris object and a spacecraft in the Sun-Synchronous orbit. The material point method is used to calculate the impact process. During the orbit prediction of the debris objects, the Earth gravitation with 10×10 EGM96, the Solar gravitation, Lunar gravitation, the atmospheric model with NRLMSISE00 and USSA1976, as well as solar radiation pressure are calculated. Relative trajectories between three debris objects are calculated to investigate the evolution of objects in the debris clouds. The relative trajectories depend on the initial velocity errors among debris objects. Although debris objects had caused by the spacecraft in the Sun-Synchronous orbit, the debris objects may stay in quasi-Sun-Synchronous orbits or non-Sun-Synchronous orbits, which depends on the initial velocity errors of these objects. About the relation between the initial velocity errors and derbis, if the velocity errors in two coordinate axes coexist, the rate of change of the right ascension of the ascending node looks obviously.