# Attitude Motion of Cylindrical Space Debris during Its Removal by Ion Beam

Authors: Vladimir S. Aslanov; Alexander S. Ledkov
Journal: Mathematical Problems in Engineering
 (2017)

---

## 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.


---

## Body

## 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 systemOxpyp. OriginO is the center of the Earth. The axis xp passes through the pericenter of initial space debris orbit. The origin of the orbital frame CxOyO is located at the center of mass of the space debris (Figure 1). The axis xO lies along the radius vector r of the space debris center of mass. The axis yO is directed towards the orbital flight. The body frame Cxbyb is fixed relative to the space debris. The active spacecraft is a material point B (Figure 1). Its radius vector has the following coordinates: ρB=CB→=0,ρBT in CxOyO frame. The direction of the ion beam axis is defined by the angle α. The axis xa is directed along ion beam axis to flight direction, and the axis ya completes the right-handed set.Figure 1
Plane motion of the space debris under the influence of ion beam forces and torques.The Lagrange formalism is used to obtain the equations of the space debris motion(1)ddt∂L∂q˙j-∂L∂qj=Qj,where L=T-U is the Lagrangian of the system, T is the kinetic energy, U is the potential energy, Qj is the nonpotential generalized forces, qj is a component of the generalized coordinates vector q=[r,ν,φ]T, r 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 xb from local vertical (Figure 1). The kinetic energy of the space debris is the sum of the kinetic energy of its center of masses TC and the kinetic energy of the motion concerning its center of mass Tr(2)TC=mr2ν˙2+r˙22,Tr=Jzν˙-φ˙22.The potential energy is [16](3)U=-μmr+3μIx-Iycos2⁡φ2r3,where μ is the gravitational constant of the Earth, Ix is longitudinal moment of inertia of the space debris, and Iy and Iz are transversal moments of inertia.To calculate the nonpotential generalized forces, the virtual work of nonpotential ion force and torque should be written(4)δW=Fi·δr+Li·δψ=∑jQjδqj.Here Fi=Fix,FiyT is the ion beam force, which is specified in the orbital reference frame, Li is the ion beam torque, ψ=ν-φ is the angle of the axis xb deflection from the axis yp, and δr is variation of r. It follows from (4) that(5)Qr=Fix,Qν=Fiyr+Liz,Qφ=-Liz.Substitute expressions (2) and (3) into (1) and the expression of second derivatives gives(6)r¨=rν˙2-μr2+9μIx-Iycos2⁡φ2mr4+Fixm,(7)ν¨=-2r˙ν˙r+3μIx-Iycos⁡φsin⁡φmr5+Fiymr,(8)φ¨=-2r˙ν˙r+3μIx-Iycos⁡φsin⁡φIz+mr2Izr5m-LizIz+Fiymr.Equations (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 as(9)Fj=-njm0Vj·NjVj,where Vj=[-u0,u0yaj/-xaj,u0zaj/-xaj] is the velocity of the ion flux at the point Pj and Pj is the barycenter of jth triangle. The velocity components are calculated according to the self-similar model [15, 19]. The barycenter point has coordinates -xaj,yaj,zaj in the frame Bxayaza; nj is the plasma density at the point Pj. The plasma density can be calculated within the self-similar model as [19](10)nj=n0R02xaj2tan2⁡α0exp⁡-3yaj2+zaj2xaj2tan2⁡α0,where n0 is the plasma density at the beginning of the far region, m0 is the mass of particle, R0 is the radius of the beam at the beginning of the far region, u0 is the axial component of the ion flus velocity, α0 is the divergence angle of the beam, αsj is the angle between the velocity Vj and the normal unit vector Nj of the jth triangle (Figure 2), and Aj is the area of jth triangle.Figure 2
Force acting onjth triangle.The ion flux torque vector, relative to the space debris center of mass, can be calculated as(11)Li=∑j∈JCP→j×Fj,where J is an index set indicating the subset of faces in the polyhedron that is inside the ion beam flow. The ion beam force is (12)Fi=∑j∈JFj.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 CxOyOzO. In the case of plane motion, generalized forces (5) contain projections Fix, Fiy, and Liz. If the considered space debris has a symmetrical shape and its axis of symmetry lies in the plane of orbital flight, then Fiz=0, Lix=Liy=0, 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.


## 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 systemOxpyp. OriginO is the center of the Earth. The axis xp passes through the pericenter of initial space debris orbit. The origin of the orbital frame CxOyO is located at the center of mass of the space debris (Figure 1). The axis xO lies along the radius vector r of the space debris center of mass. The axis yO is directed towards the orbital flight. The body frame Cxbyb is fixed relative to the space debris. The active spacecraft is a material point B (Figure 1). Its radius vector has the following coordinates: ρB=CB→=0,ρBT in CxOyO frame. The direction of the ion beam axis is defined by the angle α. The axis xa is directed along ion beam axis to flight direction, and the axis ya completes the right-handed set.Figure 1
Plane motion of the space debris under the influence of ion beam forces and torques.The Lagrange formalism is used to obtain the equations of the space debris motion(1)ddt∂L∂q˙j-∂L∂qj=Qj,where L=T-U is the Lagrangian of the system, T is the kinetic energy, U is the potential energy, Qj is the nonpotential generalized forces, qj is a component of the generalized coordinates vector q=[r,ν,φ]T, r 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 xb from local vertical (Figure 1). The kinetic energy of the space debris is the sum of the kinetic energy of its center of masses TC and the kinetic energy of the motion concerning its center of mass Tr(2)TC=mr2ν˙2+r˙22,Tr=Jzν˙-φ˙22.The potential energy is [16](3)U=-μmr+3μIx-Iycos2⁡φ2r3,where μ is the gravitational constant of the Earth, Ix is longitudinal moment of inertia of the space debris, and Iy and Iz are transversal moments of inertia.To calculate the nonpotential generalized forces, the virtual work of nonpotential ion force and torque should be written(4)δW=Fi·δr+Li·δψ=∑jQjδqj.Here Fi=Fix,FiyT is the ion beam force, which is specified in the orbital reference frame, Li is the ion beam torque, ψ=ν-φ is the angle of the axis xb deflection from the axis yp, and δr is variation of r. It follows from (4) that(5)Qr=Fix,Qν=Fiyr+Liz,Qφ=-Liz.Substitute expressions (2) and (3) into (1) and the expression of second derivatives gives(6)r¨=rν˙2-μr2+9μIx-Iycos2⁡φ2mr4+Fixm,(7)ν¨=-2r˙ν˙r+3μIx-Iycos⁡φsin⁡φmr5+Fiymr,(8)φ¨=-2r˙ν˙r+3μIx-Iycos⁡φsin⁡φIz+mr2Izr5m-LizIz+Fiymr.Equations (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 as(9)Fj=-njm0Vj·NjVj,where Vj=[-u0,u0yaj/-xaj,u0zaj/-xaj] is the velocity of the ion flux at the point Pj and Pj is the barycenter of jth triangle. The velocity components are calculated according to the self-similar model [15, 19]. The barycenter point has coordinates -xaj,yaj,zaj in the frame Bxayaza; nj is the plasma density at the point Pj. The plasma density can be calculated within the self-similar model as [19](10)nj=n0R02xaj2tan2⁡α0exp⁡-3yaj2+zaj2xaj2tan2⁡α0,where n0 is the plasma density at the beginning of the far region, m0 is the mass of particle, R0 is the radius of the beam at the beginning of the far region, u0 is the axial component of the ion flus velocity, α0 is the divergence angle of the beam, αsj is the angle between the velocity Vj and the normal unit vector Nj of the jth triangle (Figure 2), and Aj is the area of jth triangle.Figure 2
Force acting onjth triangle.The ion flux torque vector, relative to the space debris center of mass, can be calculated as(11)Li=∑j∈JCP→j×Fj,where J is an index set indicating the subset of faces in the polyhedron that is inside the ion beam flow. The ion beam force is (12)Fi=∑j∈JFj.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 CxOyOzO. In the case of plane motion, generalized forces (5) contain projections Fix, Fiy, and Liz. If the considered space debris has a symmetrical shape and its axis of symmetry lies in the plane of orbital flight, then Fiz=0, Lix=Liy=0, 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.Table 1
Space debris (Cosmos 3M) parameters.
Parameter
Value
Massm
1400 kg
Longitudinal moment of inertiaIx
1300 kg m2
Transversal moment of inertiaIy
6800 kg m2
Transversal moment of inertiaIz
6800 kg m2
Products of inertiaIxy, Ixz, Iyz
0 kg m2
Length
6.5 m
Radius
1.2 mTable 2
Parameters for ion beam force and torque calculation.
Parameter
Value
Plasma densityn0
2.6
·
1
0
16 m−3
Mass of particle (xenon)m0
2.18
·
1
0
-
25
Radius of the beam at the beginning of the far regionR0
0.1 m
Axial component of the ion flus velocityu0
38000 m/s
Divergence angle of the beam
15
°
Distance between thruster exit and debris centreρB
15 mTable 3
Parameters of existing thrusters [15].
Parameter
SPT100
T100
Thrust
84.9 mN
82.4 mN
Discharge voltage
300 V
300 V
Discharge current
4.5 A
4.5 A
Power
1350 W
1390 W
Anode flow rate
4.99 mg/s
4.93 mg/s
Cathode flow rate
0.38 mg/s
0.41 mg/s
Specific impulse
1610 s
1570 s
Thrust efficiency
49.7%
47.1%Figure 3
Sample of mesh model for ion flux force and torque calculation (856 triangles).Let us consider plane motion of the system. Figures4–6 show projections of (11) and (12) on the axis of orbital frame CxOyO 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 Cxbyb, it suffices to consider only the positive values of the angle. The values for the negative angle can be obtained from the relations(13)Fix-α=-Fixπ-α,Fiy-α=Fiyπ-α,Liz-α=-Lizπ-α.Figures 4 and 5 show that increasing the angle α leads to a shift Fix downward and a decrease in the modules Fiy. 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 α=1.87° 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 α=0, the maximum Fi force is observed. Let us consider points Ak in Figure 5. Points A1, A2, and A4 are local maxima and minima of the curve Fiy(φ) (for α=0). Figure 7 demonstrates orientations of the Cosmos 3M stage corresponding to these points in the orbital frame.Figure 4
Dependence of ion beam force projectionFix on deflection angle φ.Figure 5
Dependence of ion beam force projectionFiy on deflection angle φ.Figure 6
Dependence of ion flux torque projectionLiz on deflection angle φ.Figure 7
The positions of the Cosmos 3M stage.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 φ=0. In this case, the time of descent to an altitude of 100 km is 85 days. At the least effective case φ=86°, the descent takes 125 days.Table 4
Time of the debris removal from 500 km to 100 km altitude.
Case
φ
0
φ
˙
0
F
i
x
F
i
y
Time
Fixed orientationA1
0
0
0
−0.0444 N
84 days
Fixed orientationA2
86°
0
−0.00046 N
−0.0301 N
125 days
Fixed orientationA3
90°
0
0
−0.0306 N
123 days
Fixed orientationA4
94°
0
0.00046 N
−0.0301 N
125 days
Oscillations inZ1near separatrix (Figure 8)
45°
0
0
−0.0439 N
86 days
Oscillations inZ2near separatrix (Figure 8)
55°
0
2.027
·
10
-
5 N
−0.0429 N
88 days
Oscillations inZ1near stable pointc1(Figure 8)
1°
0
0
−0.0444 N
85 days
Oscillations inZ2near stable pointc2(Figure 8)
89°
0
5.591
·
10
-
6 N
−0.0306 N
123 days
Rotation inZ4(Figure 8)
0
0.015 s−1
1.059
·
10
-
7 N
−0.0407 N
88 daysTable4 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 α=0. For the oscillations and rotation cases Fix and Fiy columns contain average values on a period or turnover at the beginning of the space debris removing.It should be noted that forceFix takes values of the same order with Fiy (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 ρB is invariant in the orbital coordinate system CxOyOzO.In order to analyze the motion of the cylinder relative to its center of mass, let us consider Figure6. Two types of graphs are observed: when function Liz(φ) has three roots (for example, α=0 on Figure 6) when it has two roots (for example, α=5° 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 s1, s2 and three centers c1, c2, and c3 (Figure 8). This case can exist only for α=0°. The separatrices connect the saddle points and divide the phase portrait into five areas: the phase portrait contains three oscillation areas Z1, Z2, and Z3 and two rotation areas Z4 and Z5. If initial orientation and angular velocity of the space debris correspond to area Z1, Z2, or Z3, then the space debris will oscillate with relatively small amplitude around the center c1, c2, or c3, respectively. Areas Z4 and Z5 correspond to the space debris rotations in opposite directions.Figure 8
Phase portrait for case 1 (α=0).If resultant moment has four roots, there are two saddle pointss3, s4 and two centers c4, c5 (Figure 9). The phase portrait contains two oscillation areas. If resultant moment has two roots, there are one saddle point s5 and one center c6 (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.Figure 9
Phase portrait for case 2 (α=0.5°).Figure 10
Phase portrait for case 3 (α=12°).It should be noted that the centersc1 and c3 in the phase portraits are located close to the minimum points on the curve Fiy, and the center c2 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 c2 (in Z2 zone) is significantly less effective than cases of oscillations around points c1 and c3 (Table 4). It makes sense to develop control laws for angle α, which provide stabilization of the cylindrical space debris near the center c1 or c3. 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 pointc2 is significantly less effective than oscillations around stable equilibrium point c1 or c3. The results can be used to create new efficient systems of large space debris removal.



---

*Source: 1986374-2017-11-26.xml*