#### Abstract

The paper is devoted to the problem of space debris mitigation. Contactless method of the space debris deorbiting is considered. It is assumed that ion thrusters on the active spacecraft create the ion flow, which blows the debris and slows it down. The objectives of this work are the development of mathematical models and the research of space debris motion under the action of the ion flow. It is supposed that the space debris is a rigid body of a cylindrical shape. Calculation of ion beam force and torque was performed for a self-similar model of plasma plume expansion using the hypothesis of ion fully diffused reflection from a surface. A mathematical model describing plane motions of the cylindrical space debris under the influence of gravity gradient torque and the ion flux was constructed. It was shown that motion of the space debris around its center of mass has a significant effect on its removal time. Phase portraits, describing the motion of the space debris relative to its center of mass, were constructed. Comparison of the descent times in different motion modes was carried out. The results can be used to create new effective systems of large space debris removal.

#### 1. Introduction

During recent years in the scientific literature, considerable attention has been paid to the problem of transportation of nonfunctioning satellites and space debris removal. The majority of works were devoted to systems that imply a stage of docking or capturing a transported object by harpoons [1], a net [2], a tether [3], or a robotic manipulator [4]. A detailed overview of the capture tools and methods is given in the work [5]. Docking with an unmanaged object is a complex technical task. Failure at this stage is highly probable and can lead to the formation of new debris. An alternative is the use of noncontact transportation methods: based on the Coulomb interaction [6] and the ion beam created by the electric-reactive engine [7].

The use of the ion beam involves the placement of electric-reactive engines on an active spacecraft. These engines are not something exotic and are widely used in modern aerospace [8]. The engines “blow” on the transported object and thus change the parameters of its motion. To date, “Ion Beam Shepherd” is the best-designed project in this field. Existing studies show that the considered method allows removing an object with a mass of several tons from the orbit with a height of the order of 1000 km for several months [7]. Works [9, 10] are devoted to the estimation of the influence of the ion beam created by an electric-reactive engine, on objects of various forms. There are studies in which the optimal control laws of the active spacecraft for the space debris removal from orbit are developed [11, 12]. An interesting concept for detumbling spinning debris objects using the interaction between the thruster exhaust gases from the chaser and the debris object was considered in [13]. The authors have developed control function for detumble Envisat in a closed-loop simulation. Analysis of the literature shows that in existing studies due attention is not paid to the motion of the transported object relative to its center of mass.

The aim of this work is the development of mathematical models and research space debris motion under the action of the ion flow. It is supposed that the space debris is a rigid body of a cylindrical shape. The geometric parameters of the cylinder, which was taken as an example in this paper, correspond to the Cosmos 3M rocket stage. Currently, about 300 such stages are in orbit in some of the most crowded orbital regions [14]. The removal of these objects from the orbits is of great importance for the space debris mitigation issue.

#### 2. Materials and Methods

Mathematical models of a space debris plane motion under the action of the gravitational, ion flux forces, and torques will be developed in this section. Method of the ion beam resultant force and torque calculation will also be described.

##### 2.1. Equations of Plane Motion

The plane motion of space debris is considered. It is assumed that the space debris is a rigid body; the Earth does not rotate, and it has a spherical shape; the active spacecraft, which creates ion flow, is maintained in a fixed position relative to the center of mass of the space debris by its control system. It is supposed that only gravitational and ion forces and torques act on space debris.

Let us introduce the inertial coordinate system . Origin* O* is the center of the Earth. The axis passes through the pericenter of initial space debris orbit. The origin of the orbital frame is located at the center of mass of the space debris (Figure 1). The axis lies along the radius vector of the space debris center of mass. The axis is directed towards the orbital flight. The body frame is fixed relative to the space debris. The active spacecraft is a material point (Figure 1). Its radius vector has the following coordinates: in frame. The direction of the ion beam axis is defined by the angle . The axis is directed along ion beam axis to flight direction, and the axis completes the right-handed set.

The Lagrange formalism is used to obtain the equations of the space debris motionwhere is the Lagrangian of the system, is the kinetic energy, is the potential energy, is the nonpotential generalized forces, is a component of the generalized coordinates vector , is the distance between center of Earth and the space debris center of mass, is true anomaly angle, and is deflection of the space debris axis from local vertical (Figure 1). The kinetic energy of the space debris is the sum of the kinetic energy of its center of masses and the kinetic energy of the motion concerning its center of mass The potential energy is [16]where is the gravitational constant of the Earth, is longitudinal moment of inertia of the space debris, and and are transversal moments of inertia.

To calculate the nonpotential generalized forces, the virtual work of nonpotential ion force and torque should be writtenHere is the ion beam force, which is specified in the orbital reference frame, is the ion beam torque, is the angle of the axis deflection from the axis , and is variation of . It follows from (4) thatSubstitute expressions (2) and (3) into (1) and the expression of second derivatives givesEquations (6)–(8) describe plane motion of the space debris.

##### 2.2. Method of the Ion Beam Force and Torque Calculation

Plasma plume expansion into vacuum is a complex physical phenomenon, which is studied quite extensively in [15, 17]. An overview of simplified plume exploration and debris interaction models is given in [11]. Following the self-similar model, at a distance exceeding several diameters of the thruster’s nozzle, the plume has Gaussian density profile, a constant axial velocity, and a linearly increasing radial velocity [15].

Ion beam produces high velocity rarefied flow. Two simplified models of the interaction of ions with the surface are considered in the literature. The first one supposes ions spectral reflection when the resulting force is normal to the surface. The second one assumes fully diffused reflection when the force acting on the surface is directed along the velocity of the incoming particles. The experimental results reported in [18] indicate that the diffuse reflection model is in better agreement with the experimental data. The sputtering of the target material and the escaping ions from the space debris surface are neglected in the framework of this study. The surface of the body is divided into triangles. For the case of fully diffused reflection, the force acting on each triangle can be calculated aswhere is the velocity of the ion flux at the point and is the barycenter of th triangle. The velocity components are calculated according to the self-similar model [15, 19]. The barycenter point has coordinates in the frame ; is the plasma density at the point . The plasma density can be calculated within the self-similar model as [19]where is the plasma density at the beginning of the far region, is the mass of particle, is the radius of the beam at the beginning of the far region, is the axial component of the ion flus velocity, is the divergence angle of the beam, is the angle between the velocity and the normal unit vector of the th triangle (Figure 2), and is the area of th triangle.

The ion flux torque vector, relative to the space debris center of mass, can be calculated aswhere is an index set indicating the subset of faces in the polyhedron that is inside the ion beam flow. The ion beam force is

The authors developed a program in MATLAB that performs the calculation of (11) and (12) depending on the orientation of the body frame relative to the orbital frame and angle of the ion beam axis deflection . The ion flux torque and force are calculated in projections on the axis of the orbital frame . In the case of plane motion, generalized forces (5) contain projections , , and . If the considered space debris has a symmetrical shape and its axis of symmetry lies in the plane of orbital flight, then , , and ion flux torque and force do not tend to lead the space debris body out of the plane of flight. The developed program of ion flux torque and force calculation can be used in general three-dimensional cases.

#### 3. Results and Discussion

Let us consider the deorbiting of Cosmos 3M stage from the circular orbit of 500 km height. Equations (6)–(8) were integrated numerically in MATLAB. For the ion beam force and torque calculation, the program described in Section 2.2 was used. The stage was considered as a cylinder whose center of mass lies on the axis and plane of symmetry. Parameters of the stage are given in Table 1. Table 2 contains the parameters of some hypothetical electric-reactive engine, which is close in parameters to the existing engines [15]. Table 3 contains the nominal values of the main characteristics of existing thrusters, which can provide an ion beam with the parameters given in Table 2. When generating a grid of triangles, the size of the triangles of the cylindrical part was selected first. Based on this size, the mesh for the end faces is generated. An example of a grid is shown in Figure 3. The size of the triangles was chosen in such a way that the standard deviation between the ion forces, which were obtained on the current grid and the grid with twice the number of triangles, did not exceed 10^{−6} N. For calculation of ion beam force and torque the cylinder was divided into 54836 triangles. Such a large number of triangles cannot be shown in the figure, so Figure 3 is generated for 856 triangles.

Let us consider plane motion of the system. Figures 4–6 show projections of (11) and (12) on the axis of orbital frame as functions of angle for various values of , which defines ion beam axis direction. Let us investigate the influence of on the force. Since the considered body is symmetric with respect to the plane , it suffices to consider only the positive values of the angle. The values for the negative angle can be obtained from the relationsFigures 4 and 5 show that increasing the angle leads to a shift downward and a decrease in the modules . In other words, as the angle increases, the force tending to withdraw the space debris from the ion beam axis will grow, and the braking force will decrease. It should be noted that even at an angle there is a partial blowing of the stage when part of the cylinder is outside the ion flow. The kinks on the curves in the figures are caused by the transition of the end and lateral surfaces of the cylinder into the shadow. At the angle , the maximum force is observed. Let us consider points in Figure 5. Points , , and are local maxima and minima of the curve (for ). Figure 7 demonstrates orientations of the Cosmos 3M stage corresponding to these points in the orbital frame.

To estimate the time of the debris removal from the orbit, let us simulate the motion of the system using (6)–(8) at a fixed angle . Table 4 contains results of simulations. The most effective is orientation of the stage over the angle . In this case, the time of descent to an altitude of 100 km is 85 days. At the least effective case , the descent takes 125 days.

Table 4 contains three types of cases: fixed orientation, oscillations, and rotation. Fixed orientation cases assume a calculation of fixed values of angle. In the oscillation cases, oscillations along the angle near a certain equilibrium position are observed during the whole time of descent. This equilibrium is slowly displaced during the descent process. In the case of rotation, the stage rotates relative to its center of mass during the entire time of descent. All results were obtained for . For the oscillations and rotation cases and columns contain average values on a period or turnover at the beginning of the space debris removing.

It should be noted that force takes values of the same order with (Figures 4 and 5). By analogy with the case of motion in the atmosphere, we can talk about the presence of a lift force that will move the considered cylinder from the ion beam axis. In this regard, additional efforts should be taken to control the active spacecraft position and orientation. This question remains behind the scopes of this study since, in the developed models, the position of the active spacecraft is invariant in the orbital coordinate system .

In order to analyze the motion of the cylinder relative to its center of mass, let us consider Figure 6. Two types of graphs are observed: when function has three roots (for example, on Figure 6) when it has two roots (for example, on Figure 6). As the gravity gradient torque acts in addition to the ion beam torque, their ratio determines the form of the phase portrait of (8). A series of numerical calculations showed that one of three types of phase portrait could be realized. If resultant moment has three roots, there are two saddle points , and three centers , , and (Figure 8). This case can exist only for . The separatrices connect the saddle points and divide the phase portrait into five areas: the phase portrait contains three oscillation areas , , and and two rotation areas and . If initial orientation and angular velocity of the space debris correspond to area , , or , then the space debris will oscillate with relatively small amplitude around the center , , or , respectively. Areas and correspond to the space debris rotations in opposite directions.

If resultant moment has four roots, there are two saddle points , and two centers , (Figure 9). The phase portrait contains two oscillation areas. If resultant moment has two roots, there are one saddle point and one center (Figure 10). The phase portrait contains one oscillation area. Figures 9 and 10 allow us to conduct that if at the initial time the space debris rotates in the direction corresponding to the upper part of the phase portrait, then as a result of the action of the ion flow the rotation will first slow down and then will spin up in the opposite direction.

It should be noted that the centers and in the phase portraits are located close to the minimum points on the curve , and the center is located close to maximum points on this curve. In this regard, it can be concluded that the removal of space debris in the regime of oscillations near the center (in zone) is significantly less effective than cases of oscillations around points and (Table 4). It makes sense to develop control laws for angle , which provide stabilization of the cylindrical space debris near the center or . This topic will be the subject of further research.

These phase portraits were plotted on a relatively short time interval (about 3 hours). On this interval the change in the altitude of the orbit is insignificant. During the process of the space debris deorbiting, the magnitude of the gravity gradient torque will slowly increase, and the phase portraits will deform.

In this study, the self-similar model of plasma plume expansion and the hypothesis of ion fully diffused reflection from a surface were used to calculate ion beam force and torque acting on the cylindrical space debris. In subsequent studies, more attention should be paid to influence of the ions interaction with the surface of the space debris on its behavior.

The study focuses on the dynamics of a solid under the influence of ion flow and gravitational momentum. However, in low orbits, the atmosphere begins to exert a significant influence. The study of the effect of the atmosphere along with the ion flux will be the topic of our next study.

#### 4. Conclusions

The noncontact method of the space debris deorbiting by an active spacecraft with ion thrusters was considered in this paper. This is a relatively new method of space debris removal, whose main advantage means there is no need for a complex and dangerous docking stage. The mathematical models developed within the framework of the article allow studying the attitude motion of space debris under the influence of the ion beam on the orbit. Calculation of ion beam force and torque was performed for a self-similar model of plasma plume expansion using and the hypothesis of ion fully diffused reflection from a surface. It was shown that motion of the space debris around its center of mass has a significant effect on its removal time. For the case of plane motion, phase portraits describing the motion of the space debris relative to its center of mass were constructed. Comparison of the descent times in different motion modes was carried out. It was shown that, in terms of descent time, the space debris oscillations in a small neighborhood of stable equilibrium point is significantly less effective than oscillations around stable equilibrium point or . The results can be used to create new efficient systems of large space debris removal.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest regarding this research article.

#### Acknowledgments

This study was supported by the Russian Foundation for Basic Research (RFBR 17-41-630274 r-a).