1 Introduction

Interest in CubeSat development has been continuously growing, since the standard configuration was first introduced in the late 90s. A CubeSat has a cubic shape of 10 cm by side, and up to 1.3 kg of mass. These units (U) can be assembled together to form satellites of 2U, 3U, and beyond. This type of satellite was conceived for didactic purposes, however, as technology evolved, and component miniaturization continued, it was soon recognized that they had a great potential for applications with scientific and technological value. Some of these applications include: Earth observation [1], telecommunications [2], astronomy [3], and technology demonstrations. According to available records, around 700 CubeSats had been deployed into orbit as of July 2017 [4]. If those were not enough, some mega-constellations have been conceived. Such is the case of Starlink, an initiative of SpaceX, which intends to deploy 12,000 mini-satellites to provide satellite Internet, and OneWeb, which aims to deploy 2700 mini-satellites also for telecommunications purposes [5]. While the last two constellations are not composed by CubeSats, it is clear that the miniaturization trend came to stay, and that the number of mini/nano-satellites in orbit is going to grow exponentially.

This progress has come at a cost, as most nanosatellites are not equipped with any deorbiting system. Depending on the orbital altitude at which it is deployed, it can take a CubeSat from a few months to even centuries before re-entering Earth’s atmosphere. Therefore, there exists the concern that the CubeSat boom can have a negative impact in the space junk problem, which is already an issue recognized by the scientific community. A total of 24,000 objects the size of a baseball or larger were estimated to be in Earth orbit as of 2011 [6], most of them without any remaining means of orbital control, which increases the chances of collisions in space, as have already occurred [7]. Furthermore, the Kessler syndrome predicts that a critical density of objects in Low Earth Orbit could be reached, and that collisions between them could cause a cascade of impacts, that would eventually create enough debris to render some orbital ranges unusable within this century [8]. In response to this, international guidelines dictate that any artificial satellite must be deorbited within 25 years of the end of their mission [9]. These considerations make clear the necessity to develop technologies that enable CubeSats to be deorbited in a timely manner.

Some deorbiting systems aimed for CubeSats have been conceived already. Most notable of these approaches include: deorbit sails, inflatables, and electric tethers. Deorbit sails rely on aerodynamic drag force augmentation, such force can be expressed as:

$$\begin{aligned} \mathbf {F_{D}}=-\frac{1}{2}\rho \Vert {\mathbf {v}}_\text {rel}\Vert {\mathbf {v}}_\text {rel} C_{D} A, \end{aligned}$$
(1)

where \(\mathbf {F_{D}}\) is the drag force in Newtons (N), \(\rho \) is the mass density of the fluid in kg/m\(^3\), \({\mathbf {v}}_\text {rel}\) is the velocity vector relative to the atmosphere in m/s, \(\Vert \cdot \Vert \) is the 2-norm of the vector, \(C_{D}\) is the drag coefficient of the object, and A is the cross-sectional area of the object in m2.

In the case of inflatable structures, the aim is to decrease the ballistic coefficient of the spacecraft. The ballistic coefficient, BC, of an object is given by the following formula:

$$\begin{aligned} BC = \frac{M}{C_{D}A}, \end{aligned}$$
(2)

where BC is in kg/m\(^2\) and M is the mass of the object in kg.

Tethers, on the other hand, rely in the electrodynamic force created as a result of the interaction between a current carrying cable with Earth’s magnetic field. Such force can be expressed as:

$$\begin{aligned} \mathbf {F_{e}} =-I_\text {ave}L{\mathbf {b}}\times {\varvec{\ell }}, \end{aligned}$$
(3)

where \(\mathbf {F_{e}}\) is the electrodynamic force in N, \(I_\text {ave}\) is the average current across the electrodynamic tether in A, L is the length of the tether in meters, \({\mathbf {b}}\) is the local geomagnetic field vector in Tesla (T), and \({\varvec{\ell }}\) is the unit vector along the electrodynamic tether.

Some of the most notable missions and concepts related to each of these technologies are as follows: In the field of sails, NanoSail-D and NanoSail-D2 were missions developed by NASA [10]. NanoSail-D launch was unsuccessful. Its ground spare, NanoSail-D2 made it into orbit in November 2010. The experiment was successful, since the satellite was deorbited from an initial altitude of 650 km in just 240 days. This process would have taken around 25 years if left to natural orbital decay. Another example of this technology is University of Surrey’s DeorbitSail satellite [11], launched in 2015. However, this experiment was unsuccessful, since the sail failed to deploy. Finally, the University of Glasgow developed the Aerodynamic End of Life Deorbit System for CubeSats (AEOLDOS) concept [12]. To the best knowledge of the authors, such system has yet to be demonstrated in space.

The efficiency of different shapes of inflatable structures has also been the subject of study. In their work, Maessen et al. [13] conclude that a pyramidal structure is optimal. Nakasuka et al. [14], on the other hand, investigate the effectiveness of a spherical structure, and propose a system that is claimed to be effective at initial altitudes of up to 800 km. Lokcu and Ash [15] experimented with three different geometries, spherical, pyramidal, and pillow shapes, claiming their system would be effective to deorbit CubeSats from an initial altitude of 900 km in a period of 30 years. Andrews et al. [16] propose a cone-shaped system that, besides deorbiting the nanosatellite, would also protect it during atmospheric re-entry, allowing payloads to be recovered. A combination between inflatables and sails has also been proposed by Viquerat et al. [17], for its InflatSail concept.

In the field of electrodynamic tethers, the study by Voronka et al. [18] claims that deorbit times of 25 years for CubeSats at an initial altitude of 1000 km can be attained deploying a 1 km-long tether. Zhu and Zhong have proposed different control approaches to the deorbit problem with electrodynamic tethers. One of these approaches is an On-Off current scheme [19], with which is claimed that a CubeSat would lose altitude at a rate of 100 km per 60 days. The second approach uses a finite receding horizon control [20]. In the latter case, the CubeSat tethered system would lose 100 km of altitude in 25 days. There are also approaches that focus in the tether interaction with space plasma rather than the Earth’s magnetic field, such is the case of the work by Janhunen and Khurshid et al. [21]. This approach is known as plasma brake, and will be tested in the Aalto-1 CubeSat mission [22].

All these three approaches require the deployment of actuators in space, such as sails, inflatables or tethers making the system prone to failures. They also require significant volumes for storage of their respective actuators. These are the major disadvantages of said approaches, reason why new alternatives are being investigated. Specifically, one of the technologies with more potential in this area is electric propulsion. Electric engines do not require the deployment of any actuator once in orbit, they are compact modules that only require the application of energy to generate thrust. Also, with a few exceptions, such as the micro Cathode Arc Thruster and micro Pulsed Plasma Thruster, which include springs in their design, they do not have movable parts. Important efforts have been made to develop electric thrusters that are small and efficient enough to be integrated in nanosatellites. Some of the most promising approaches of this technology are discussed in the next section.

One of the main challenges to use thrusters for deorbiting purposes, is that CubeSats often lack attitude determination and control capabilities [23]. In such cases, an alternative method is required to ensure that the thrusters are pointed in such a way that the thrust vector opposes in some degree to the velocity vector, causing the loss of orbital energy, eventually deorbiting the spacecraft. This paper presents a new approach that tackles this problem in two stages: first, to spin-stabilize the satellite, which inertially fixes the direction of the spinning axis; second, an orbit sampling phase, with the aim of identifying the portions of the orbit when the thruster is pointed in the required direction. A gyroless spin stabilization scheme is proposed, using only magnetometers and magnetorquers as sensors and actuators, respectively. As for the thrusting phase of the algorithm, only Global Positioning System (GPS) receiver data are needed as an input. In addition to the minimal requirements in terms of sensors and actuators, both algorithms are computationally cheap, which makes them ideal for their application in CubeSats.

The contents of this paper are organized as follows: Sect. 2 offers a summary of the hardware required by the proposed algorithm, as well as the basic satellite attitude kinematics and dynamics. The spin-stabilization algorithm is described in Sect. 3, whereas the deorbiting phase algorithm is presented in Sect. 4. Sections 3 and 4 present validation results of the proposed algorithms through numerical simulations. Finally, conclusions are presented in Sect. 5.

2 Hardware and Environment

Table 1 Electric engine types [24]
Full size table

2.1 Magnetic Torquers and Magnetometers

Attitude estimation algorithms such as TRIAD and QUEST [25] require a minimum of two vector measurements provided by separate sensors. Since CubeSats have tight constraints in terms of mass, volume, and energy, it is not always possible to equip them with enough of these sensors. Therefore, in many cases, they lack attitude determination and control capabilities. This fact makes any orbital maneuver, including deorbiting, a rather difficult task. There are, however, some components that are standard in many CubeSat missions, such as magnetorquers and magnetometers, which can provide a certain degree of attitude control.

Magnetorquers are a common choice for attitude control in CubeSats. Their design is simple, consisting of an electromagnetic coil, they are small and light, with no movable parts, and consume no fuel. They do consume electrical power, which can be generated from solar panels present in virtually any CubeSat mission. All of these features make magnetorquers ideal for their integration in nanosatellites.

Magnetorquers work by generating a magnetic dipole, which in turn interacts with Earth’s magnetic field, producing a net torque in the spacecraft which can be expressed as:

$$\begin{aligned} \varvec{\tau } = {\mathbf {m}}\times {\mathbf {b}}, \end{aligned}$$
(4)

where \(\varvec{\tau }\) is the \(3 \times 1\) torque vector in N\(\cdot \)m, \({\mathbf {m}}\) is the \(3 \times 1\) magnetic dipole vector produced by the magnetorquers in Ampere(A)\(\cdot \)m2, and \({\mathbf {b}}\) is the \(3 \times 1\) Earth’s magnetic field vector in T.

Magnetic torquers work in conjunction with magnetometers, which are sensors in charge of providing the geomagnetic field vector measurement, \({\mathbf {b}}\). It is worth noting that, in general, magnetic torquers and magnetometers cannot be activated at the same time, since the magnetic torquers dipole would interfere with the magnetometers readings. Therefore, a duty cycle approach is commonly implemented, where actuators and sensors are activated alternately.

2.2 Electric Propulsion

Up until now, most CubeSats have lacked any means of propulsion, and as a result, they do not have orbital maneuvering capabilities. Among others, this has prevented the development of efficient deorbiting systems for this kind of satellites, and usually, this process is left to natural orbital decay. While this method can be effective for CubeSats at an orbital altitude of 650 km or lower, for higher orbits, it would take more than 25 years for the deorbiting to be completed.

In a response to this problem, much research has been performed in the field of electric propulsion. While there are different technologies for electric thrusters, the following ones are being developed for its application in CubeSats: Electrospray [26], Micro Pulsed Plasma Thruster [27], Hall Effect Thruster [28], CubeSat Ambipolar Thruster [29], and Micro Cathode Arc Thruster [30]. Working principles and specifics of each of these technologies are beyond the scope of this paper, and they can be found in the cited references.

Table 1 [24] summarizes some of the main features for each particular type of electric engine, which provide an insight of the potential that this technology has for applications in the field of CubeSats.

Each of these technologies is at different stages of development and provides a variety of features and efficiency. In this paper, Electrospray is considered for the numerical evaluation of the proposed algorithm, because of their thrust levels, simplicity, and low volume. However, in principle, any kind of thruster can be used in conjunction with the algorithm to be described.

2.3 Coordinates Systems and Geomagnetic Field

Two coordinate systems are used in this work. The Earth Centered Inertial (ECI) frame has its origin at the center of the Earth, the positive x-axis points towards the point of vernal equinox, the positive z-axis is aligned with the Earth rotation axis, and the positive y-axis completes the right-handed frame. Likewise, the Body (B) frame has its origin at the center of mass of the satellite, and in this paper, its axes are assumed to be aligned with the principal axes of the satellite.

As mentioned above, magnetorquers dipoles interact with the Earth’s magnetic field to generate torques that act on the satellite. In this paper, the Earth’s magnetic field is simulated using the International Geomagnetic Reference Field (IGRF-12) model provided by the International Association of Geomagnetism and Aeronomy (IAGA) [31].

2.4 Attitude Kinematics and Attitude Dynamics

Satellite attitude kinematics are described by the following equation:

$$\begin{aligned} \dot{\mathbf {q}}=\frac{1}{2} \begin{bmatrix} -[{{\varvec{\omega }}}\times ] &{} \varvec{\omega } \\ -\varvec{\omega ^{T}} &{} 0 \end{bmatrix} {\mathbf {q}}, \end{aligned}$$
(5)

where \({\mathbf {q}}\) is the 4x1 quaternion representing the relative attitude of the B-frame with respect to the ECI frame, \(\dot{\mathbf {q}}\) is the time derivative of the quaternion, \(\varvec{\omega }\) is the \(3 \times 1\) angular velocity vector of the body frame with respect to the inertial frame expressed in the body frame, \((\cdot )^T\) is the transpose, and \([\varvec{\omega }\times ]\) is the \(3 \times 3\) matrix defined, such that \([{\mathbf {a}}\times ]{\mathbf {b}} = {\mathbf {a}} \times {\mathbf {b}}\), \({\mathbf {a}}\) and \({\mathbf {b}}\) are \(3 \times 1\) vectors.

Attitude dynamics of the satellite are expressed through the Euler’s Rotation Equation as follows:

$$\begin{aligned} \varvec{\dot{\omega }} = \mathbf {J^{-1}}(-\varvec{\omega }\times {\mathbf {J}}\varvec{\omega }+\varvec{\tau }_\text {tot}), \end{aligned}$$
(6)

where \({\mathbf {J}}\) is the \(3 \times 3\) inertia tensor of the satellite, and \(\varvec{\tau }_\text {tot}\) is the \(3 \times 1\) vector representing the sum of all external torques, which include control torques as well as perturbation torques caused by the environment such as solar radiation pressure, gravity gradient, and atmospheric drag. Only control torques were considered for the simulations that are presented later in this paper.

2.5 Orbital Dynamics

To be able to develop any deorbiting algorithm, first, it is necessary to define the equations that govern the motion of a body orbiting the Earth. There are a number of different perturbation accelerations that affect the orbit of a body around the Earth. Examples of these accelerations are: lunar and solar gravitational pulls, solar pressure, accelerations caused by Earth’s oblateness, of which \(J_{2}\) term is the most significant one, and atmospheric drag. In the Low Earth Orbit environment, however, other bodies gravitational pulls and solar pressure can often be neglected, and the \(J_{2}\) term and atmospheric drag are the perturbations which dominate over the rest. When the \(J_{2}\) acceleration is plugged into the equations of motion, the following set of equations are obtained in the ECI frame, defining the motion of a satellite around the Earth [32]:

$$\begin{aligned} a_{x}= & {} {} \mu \left[ -\frac{x}{r^3}+A_{J2}\left( 15\frac{xz^2}{r^7}-3\frac{x}{r^5}\right) \right] , \end{aligned}$$
(7)
$$\begin{aligned} a_{y}= & {} {} \mu \left[ -\frac{y}{r^3}+A_{J2}\left( 15\frac{yz^2}{r^7}-3\frac{y}{r^5}\right) \right] , \end{aligned}$$
(8)
$$\begin{aligned} a_{z}= & {} \mu \left[ -\frac{z}{r^3}+A_{J2}\left( 15\frac{z^3}{r^7}-9\frac{z}{r^5}\right) \right] , \end{aligned}$$
(9)

where x, y, and z are the coordinates in the ECI frame, r is the distance between of the center of mass of the two bodies, \(a_{x}, a_{y}\), and \(a_{z}\) are the acceleration components in the ECI frame, \(A_{J2}=\frac{1}{2}J_{2}R_{e}^2\) in m2, the coefficient \(J_{2}\) being equal to 1.08263 \(\times \) 10-3, and \(R_{e}\) is the mean radius of the Earth, being equal to 6371.2 \(\times \) 103 meters.

3 Spin-Stabilization Controller

In this section, a new control law for CubeSat spin-stabilization is introduced, which is later applied to the CubeSat deorbiting problem. The proposed controller is inspired in the work of Avanzini and Giulietti [33], who develop their version of a B-dot detumbling control algorithm, originally proposed by Stickler and Alfriend [34]. The control law in [33] is given by:

$$\begin{aligned} {\mathbf {m}} = -\frac{k}{ \Vert {\mathbf {b}} \Vert ^2}({\mathbf {b}}\times \varvec{\omega }). \end{aligned}$$
(10)

The objective of a B-dot controller is to dissipate the angular velocities along all of the three main axes of a satellite, usually during the early stages of a mission. In our case, however, the objective is to dissipate the angular velocities along two of the main axes of the CubeSat, while bringing the angular velocity of the remaining axis to a desired value greater than zero, to spin-stabilize this axis in the inertial frame. We, therefore, introduce the variable \(\varvec{\omega }_\text {d} = [\omega _{dx}~0~ 0]^T\), which represents the desired angular velocity vector. It must be noted that the spinning axis must be the thruster carrying axis. Then, we can propose the following control law:

$$\begin{aligned} {\mathbf {m}} = -\frac{k}{ \Vert {\mathbf {b}} \Vert ^2}[{\mathbf {b}}\times (\varvec{\omega }-\varvec{\omega }_\text {d})]. \end{aligned}$$
(11)

Rearranging:

$$\begin{aligned} {\mathbf {m}} = -\frac{k}{ \Vert {\mathbf {b}} \Vert ^2}[{\mathbf {b}}\times \varvec{\omega }-{\mathbf {b}}\times \varvec{\omega }_\text {d}]. \end{aligned}$$
(12)

We have the following basic theory of kinematics [34]:

$$\begin{aligned} \frac{d^{ECI}}{dt}{\mathbf {b}} = \frac{d^{B}}{dt}{\mathbf {b}}+\varvec{\omega }\times {\mathbf {b}}. \end{aligned}$$
(13)

Therefore, we can use of the approximation \({\mathbf {b}}\times \varvec{\omega } \approx \dot{{\mathbf {b}}}\), and the following is derived:

$$\begin{aligned} {\mathbf {m}} = \frac{k}{ \Vert {\mathbf {b}} \Vert ^2}[({\mathbf {b}}\times \varvec{\omega }_\text{ d})-\dot{\mathbf {b}}], \end{aligned}$$
(14)

where \(\dot{\mathbf{b}}\) denotes the Earth’s magnetic field time derivative in the body frame. It can be seen that this control law does not require readings of the angular velocity; therefore, there is no need for gyroscopes. The magnetic field derivative can be obtained numerically as follows:

$$\begin{aligned} \dot{\mathbf {b}} \approx \frac{{\mathbf {b}}(t) - {\mathbf {b}}(t-\Delta t)}{\Delta t}. \end{aligned}$$
(15)

It is well known that derivative methods are susceptible to noise. As pointed out in [35], this fact is accounted for with the inclusion of a low pass filter, easily integrable into a CubeSat. Robustness of the B-dot controller in the presence of noise is also studied in [36, 37], and [38].

To test the effectiveness of the algorithm, the following scenario is designed: A 3U CubeSat, equipped with three orthogonal magnetorquers, a three-axis magnetometer, and a GPS receiver is considered. The CubeSat has a mass of 3.5 kg and the inertia tensor is equal to diag(0.01, 0.0506, 0.0506) kg\(\cdot \)m2, where diag(\(\cdot \),\(\cdot \),\(\cdot \)) is the diagonal matrix whose diagonal terms are given in the arguments. It is deployed at an initial orbital altitude of 500 km, in a near circular orbit, and with an inclination of 65\(^\circ \). The desired angular velocity is set to \(\varvec{\omega }_\text {d}=[10~ 0~ 0]\) \(^\circ \)/s, considering an electrospray thruster to be mounted in the \(+x\) face of the CubeSat body frame. A duty cycle scheme is implemented to allow the magnetometers and magnetorquers to work in conjunction, where the magnetorquer is active for 4 s, and the magnetometer is active for 1 s at a time. The magnetic dipole generated by the magnetorquers is between the values of \(\pm 0.2\) A \(\cdot \) m\(^2\), a reasonable range for current magnetorquers [39]. This is the scenario considered for all simulations within this paper, and they were performed using Matlab/Simulink with the default numerical integration algorithm and the relative tolerance of 0.001.

Figure 1 shows the evolution of the angular velocities for each axis. Observe that the desired angular velocities are acquired after around 67 min or less than one orbital period. It can also be noted that \(\omega _{x}\) attains a value slightly higher than the desired 10\(^\circ \)/s. This is a consequence of the duty cycle approach, which prevents the control dipole from being applied continuously. However, its value is stable. High level of accuracy for the spinning rate is not critical for the deorbit phase of the algorithm, so this is not considered as a problem. Figure 2 depicts the time history of the control magnetic dipoles applied during the execution of the algorithm. It is worth noticing that even when we are dealing with a scenario where no external disturbances are considered, continuous control dipoles are needed. This is also due to the use of a duty cycle approach, as well as the fact that the CubeSat is underactuated, since the magnetorquers cannot generate torques in arbitrary directions, and thus, an ideal torque cannot be achieved. However, once the angular velocity vector is sufficiently close to the desired one, the torquers can be turned off to save electric power and the deorbiting phase of the algorithm can be executed, which is presented in Sect. 4.

Fig. 1
figure 1

Angular velocities

Full size image
Fig. 2
figure 2

Magnetic dipoles

Full size image

3.1 Stability

The stability proof of the spinning algorithm is provided by means of Lasalle’s invariance principle [40, 41]. This principle establishes that a controller is stable if it fulfills the following three conditions: (1) a scalar function V(x) (known as Lyapunov function), can be defined, such that it is positive definite, (2) its time derivative \(\dot{V}\) is negative semi-definite, and (3) no trajectory can stay at points where \(\dot{V}(x) = 0\), except at the origin.

Bearing this in mind, we can express the control torque \({{\mathbf {M}}}\) as follows:

$$\begin{aligned} {\mathbf {M}} = (-k I -\hat{{\mathbf {b}}}\hat{{\mathbf {b}}}^T)(\varvec{\omega }-\varvec{\omega }_\text {d}), \end{aligned}$$
(16)

where I is a \(3 \times 3\) identity matrix, and \(\hat{{\mathbf {b}}} = {\mathbf {b}}/\Vert {\mathbf {b}}\Vert \) is the normalized magnetic field vector in the B frame.

Define \(\varvec{\omega }_\text {e} = {\varvec{\omega }}-{\varvec{\omega }}_\text {d}\), and, hence, \(\dot{\varvec{\omega }}_\text {e}\) is equal to \(\dot{\varvec{\omega }}-\dot{\varvec{\omega }}_\text {d}\), where \(\dot{\varvec{\omega }}_\text {d} = 0\). Then, we propose the following Lyapunov function:

$$\begin{aligned} V = \frac{1}{2} \varvec{\omega }_\text {e}^T {\mathbf {J}} \varvec{\omega }_\text {e}. \end{aligned}$$
(17)

Taking its time derivative, we have:

$$\begin{aligned} \dot{V}=\varvec{\omega }_\text {e}^T {\mathbf {J}} \dot{\varvec{\omega }}_\text {e} = \varvec{\omega }_\text {e}^T {\mathbf {J}} \dot{\varvec{\omega }}. \end{aligned}$$
(18)

Then:

$$\begin{aligned} \dot{V}=\varvec{\omega }_\text {e}^T {\mathbf {M}}. \end{aligned}$$
(19)

And, finally, we obtain:

$$\begin{aligned} \dot{V}= -k \varvec{\omega }_\text {e}^T (I-\hat{{\mathbf {b}}}\hat{{\mathbf {b}}}^T)\, \varvec{\omega }_\text {e}, \end{aligned}$$
(20)

where \(\dot{V}\) is negative semi-definite.

Finally, it can be observed that \(\dot{V}(x) = 0\) in two scenarios: 1) at the origin \(\varvec{\omega }_\text {e} = 0\), and 2) when the Earth’s magnetic field \({\mathbf {b}}\) is parallel to \(\varvec{\omega }_\text {e}\). It is clear, however, that the second scenario cannot be maintained in practice, since the magnetic field vector follows a very complex pattern, which would prevent \(\varvec{\omega }_\text {e}\) from tracking it as the satellite orbits the Earth.

4 Deorbiting Algorithm

For any propulsive deorbiting system to be effective, proper orientation of the satellite is required. When no attitude determination capability is available, however, the development of deorbiting algorithms becomes a challenge. In this section, we present a new algorithm which allows a spin-stabilized CubeSat to be deorbited, without the need for attitude determination capabilities, using only GPS readings.

In this phase of the mission, the CubeSat is assumed to be spin-stabilized around its thruster carrying axis. The controller presented in Sect. 3 can be used to accomplish this goal. Without full attitude determination, all that is known is the thruster must point in a suitable direction for deorbiting during half an orbit, as depicted in Fig. 3, where the blue arrows represent the velocity vector, whereas the red arrows represent the thrust vector. In this case, the upper half of the orbit will be the optimal region to activate the thrusters, because this is when the thrust vector opposes the velocity vector. The goal is, therefore, to find such portion of the orbit. We take an orbit sampling approach, and the steps of the proposed algorithm are listed next.

Fig. 3
figure 3

Optimal portion of the orbit to thrust (upper half)

Full size image
  1. 1.

    Define the values \(\nu _{\mathrm{on}}\) and \(\nu _{\mathrm{off}}\), where \(\nu _{\mathrm{off}}\) = \(\nu _{\mathrm{on}}\) + \(\pi \). These variables represent the true anomalies at which the thrusters will be turned on and off, respectively, effectively thrusting during half an orbit. At the beginning of the deorbiting operation, these values can be chosen arbitrarily.

  2. 2.

    Once \(\nu _{\mathrm{on}}\) is reached, turn on the thrusters. Keep track of the semi-major axis at the beginning of the thrusting phase, \(a_{i}\).

  3. 3.

    Once \(\nu _{\mathrm{off}}\) is reached, turn off the thrusters, and keep them in that state for the remainder of the orbit. Keep track of the semi-major axis at the end of the orbit, \(a_{f}\).

  4. 4.

    At the end of the orbit, and once \(a_{i}\) and \(a_{f}\) are determined, compare their values.

  5. 5.

    If \(a_{f}\) presents a decrease with respect to \(a_{i}\), this implies that, for the most part, we are thrusting in the correct portion of the orbit, and the orbit is descending. This scenario is depicted in Fig. 4. Retain values of \(\nu _{\mathrm{on}}\) and \(\nu _{\mathrm{off}}\) and go back to step 2.

  6. 6.

    If, on the other hand, \(a_{f}\) presents an increase with respect to \(a_{i}\), this implies that, for the most part, we are thrusting in the wrong portion of the orbit, and the orbit is rising. This scenario is depicted in Fig. 5. In this case, \(\nu _{\mathrm{on}}\) and \(\nu _{\mathrm{off}}\) must be updated in the following manner: \(\nu _{\mathrm{on}\_new}\) = \(\nu _{\mathrm{on}}\) + \(\pi \)/2 and \(\nu _{\mathrm{off}\_new}\) = \(\nu _{\mathrm{off}}\) + \(\pi \)/2. The reason to add \(\pi \)/2 to the new parameters is to sample the orbit in four quadrants. Go back to step 2.

The semi-major axis can be computed with the following formula:

$$\begin{aligned} a = \left( \dfrac{2}{r}-\dfrac{v^2}{\mu }\right) ^{-1}, \end{aligned}$$
(21)

where a is the semi-major axis in meters, r is the orbital radius magnitude in meters, v is the velocity vector magnitude in m/s, and \(\mu \) is Earth’s gravitational constant. Both r and v vectors are available from GPS readings. Notice how one of the effects that the \(J_{2}\) term introduces is an oscillation in the semi-major axis of the orbit. This is clearly seen in Figs. 4 and 5.

Fig. 4
figure 4

Sampling of the orbit. Orbit descending case, where \(a_{f} < a_{i}\). Note the sinusoidal component in the semi-major axis caused by the \(J_{2}\) term

Full size image
Fig. 5
figure 5

Sampling of the orbit. Orbit rising case, where \(a_{f} > a_{i}\). Note the sinusoidal component in the semi-major axis caused by the \(J_{2}\) term

Full size image

It is noted that, even once condition required by step 5 is met, as the deorbiting process continues, the orbital parameters will change over time, and condition in step 6 may arise again. However, it has been observed in simulations that the algorithm is robust to update the values of \(\nu _{\mathrm{on}}\) and \(\nu _{\mathrm{off}}\), such that the deorbiting maneuver is always successful.

One special case must be noted, when the spinning axis is perpendicular to the orbital plane, the deorbiting operation could not be executed with the algorithm described. It is, however, improbable to encounter this scenario in practice, and the solution of this special case is left as future work.

This algorithm was applied to the CubeSat system described in Sect. 3. The results of these simulations are shown in Figs. 6 and 7, which depict the evolution of the satellite orbital altitude and semi-major axis, respectively. Figure 6 shows how the perigee of the orbit is decreasing in a sustained way. There are a few points to be discussed in Fig. 7. First of all, the effects of the \(J_{2}\) perturbation are easily seen as they add a sinusoidal component to the semi-major axis; also, a key feature of the deorbit algorithm can be observed. At around three and a half days into the deorbiting operation, it is observed that the semi-major axis starts to increase, however, the algorithm is able to correct this condition, and bring the satellite back to a steady altitude decrease. In this case, the deorbit operation is achieved within 8 days. It is worth to stress the fact that the deorbiting time may vary from one scenario to another, as there is no control over the final orientation of the spinning axis. However, the algorithm is effective in deorbiting the CubeSat every time, which is demonstrated later when the results of the Monte Carlo simulations are shown.

Fig. 6
figure 6

CubeSat altitude. Perigee of the CubeSat is continuously decreasing until reentering the atmosphere at an altitude of 100 km

Full size image
Fig. 7
figure 7

Semi-major axis. It can be seen that three and a half days into the deorbiting process, the algorithm has to perform a new sampling of the orbit to continue with the deorbiting process

Full size image

An interesting parameter to look at is the efficiency per orbit that this algorithm can achieve, this is, how much of the time the thrusters actual state (on/off) coincides with the desired state. This metric is dependent on initial conditions, and, therefore, will be slightly different each time; however, Fig. 8 depicts the evolution of this value for the current scenario. In this case, initial efficiency is high. However, as the process continues, and orbital parameters evolve, this efficiency eventually drops to a value of around 50%. At this point, the algorithm detects that the values of \(\nu _{\mathrm{on}}\) and \(\nu _{\mathrm{off}}\) need to be updated, and starts sampling the orbit again. After four orbits, it converges to an efficiency of above 90%. Once this is achieved, the high efficiency remains for the rest of the deorbiting phase. Notice how the semi-major axis evolution shown in Fig. 7 is consistent with the efficiency depicted in Fig. 8.

Fig. 8
figure 8

Deorbiting algorithm efficiency per orbit

Full size image
Table 2 Uncertainties vector [24]
Full size table

4.1 Robustness Analysis

Monte Carlo simulations were executed to prove the robustness of the algorithm. An uncertainties vector is defined as per Table 2 [24], which includes reasonable levels of uncertainties for real life missions. The results of the Monte Carlo simulations are depicted in Fig. 9. Deorbiting times between 6 and 35 days are attained. Such wide range of deorbiting times comes from the fact that the algorithm is sensitive to the orientation of the spinning axis and we have no control over the latter. The more this spinning axis tends to lie in the orbital plane, less time it will take to deorbit the CubeSat, since thrust is being used more efficiently. Also, notice how these times are considerable shorter than those achievable with the alternative technologies mentioned in the Introduction of this paper. The most important fact to highlight from these simulations, however, is that the algorithm is successful in deorbiting the CubeSat every time, even in the presence of uncertainties.

Fig. 9
figure 9

Robustness analysis through Monte Carlo simulations

Full size image

4.2 Deorbiting Times Comparison

For the sake of performance comparison of the different deorbiting approaches, Table 3 lists common deorbiting times for each of the approaches discussed in this paper. It is clearly seen that the electric propulsion not only overcomes the practical disadvantages of the other technologies, but also outperforms them in terms of deorbiting times.

Table 3 Deorbiting times of different technologies
Full size table

5 Conclusion

Both a spin-stabilization controller and a deorbiting algorithm for CubeSats have been presented. They can be used in conjunction to achieve the deorbiting of a satellite when full attitude information is not available. Both are ideal for application in CubeSats, as they take into account hardware and software limitations present in this type of satellites. The spin-stabilization algorithm has the advantage of being gyroless, and only requires the magnetic field measurement as an input, whereas the deorbiting algorithm only requires GPS information readings. The effectiveness of both algorithms is proven through numerical simulations. CubeSat deorbiting has been suggested in the past as one of the potential applications for novel electric propulsion technology; however, this is one of the first works that explores practical methods for its implementation. While electrospray features were used for the simulations presented in this paper, effectiveness of the deorbiting algorithm is independent from the type of electric engine being used. Stability of the spinning algorithm is provided by the means of Lasalle’s invariance principle. Robustness against model uncertainties is proven through Monte Carlo simulations. Deorbit times between 6 and 35 days are obtained, which represent a substantial improvement respect to the performance of other approaches such as sails, inflatables, and electric tethers. Nonetheless, there exists the potential for the spin-stabilization controller to also be employed in conjunction with either sails or inflatables, to improve the performance of such approaches. This could be the subject of future work, as well as looking for ways to improve the efficiency of the orbit sampling algorithm. These algorithms can be applied as part of the efforts to tackle the space junk problem.