Trajectory Approximation of a Coulomb Drag-Based Deorbiting

Abstract

The presence of a number of space debris in low Earth orbits poses a serious threat for current spacecraft operations and future space missions. To mitigate this critical problem, international guidelines suggest that an artificial satellite should decay (or be transferred to a graveyard orbit) within a time interval of 25 years after the end of its operative life. To that end, in recent years deorbiting technologies are acquiring an increasing importance both in terms of academic research and industrial efforts. In this context, the plasma brake concept may represent a promising and fascinating innovation. The plasma brake is a propellantless device, whose working principle consists of generating an electrostatic Coulomb drag between the planet’s ionosphere ions and a charged tether deployed from a satellite in a low Earth orbit. This paper discusses an analytical method to approximate the deorbiting trajectory of a small satellite equipped with a plasma brake device. In particular, the proposed approach allows the deorbiting time to be estimated through an analytical equation as a function of the design characteristics of the plasma brake and of the satellite initial orbital elements.

1. Introduction

The number of near-Earth space debris constituted of out-of-order satellites, launcher upper stages, fragments produced by collisions, and other man-made objects, has been dramatically increasing during the last decades [1]. Recent estimates, based on statistical models, suggest that the number of debris objects currently orbiting the Earth is above 131 million [1], of which about 1 million have a characteristic dimension in the range between $1\mathrm{cm}$ and $10\mathrm{cm}$.

This issue could jeopardize future launches and space operations [2,3] and is highly problematic for heavily populated orbital ranges [4,5], such as the Low Earth Orbit (LEO) region. The threat posed by space debris is further worsened by the Kessler syndrome scenario, where collisions between existing debris generate a dramatically increasing cascade effect [6]. Therefore, an end-of-life deorbiting strategy must be carefully considered during the mission design phase [7] in order to guarantee either a passive or an active decay within 25 years after the end of the operative life, in accordance with the Inter-Agency Space Debris Coordination Committee (IADC) guidelines [8,9].

In this context, one promising and innovative deorbiting propellantless concept is the electrostatic plasma brake (PB), which is theoretically capable of deorbiting a spacecraft from a LEO by providing a decelerating (drag) force and without the need for any propellant consumption. The PB concept derives from the working principle of the Electric Solar Wind Sail (E-sail) [10,11], an interplanetary propellantless propulsion system proposed in 2004 by Dr. Pekka Janhunen, which consists of a spinning grid of conducting tethers kept at a high potential by a voltage source and immersed in the solar wind flux. The incoming ions exchange momentum with the charged grid due to Coulomb collisions, thus generating a continuous propulsive acceleration [12]. The peculiarity of the E-sail enables heliocentric missions that could be difficult (or impossible) to achieve with conventional thrusters, including non-Keplerian orbits [13,14], asteroid deflection [15] and outer solar system exploration [16,17,18,19].

Unlike the grid of tethers used in an E-sail arrangement, the typical configuration of a PB device [20] consists of a single electrostatically charged tether, which is unreeled in the highly-ionized upper stages of the Earth’s atmosphere, as illustrated in Figure 1.

The interaction between the conducting tether and the upper Earth’s atmosphere ions produces a dissipative force, usually referred to as Coulomb drag, which reduces the spacecraft orbital energy and lowers its perigee altitude until the increasing effect of the atmospheric drag becomes sufficient to complete the orbital decay. Preliminary numerical simulations give encouraging results on the potentialities of the PB-enabled deorbiting concept [21,22], which will hopefully be validated by experimental evidence.

The first test mission of the PB technology was planned to be the Estonian satellite EstCube-1 [23], but a failure occurred in the conducting tether unreel mechanism [24]. More recently, a similar problem in the unreeling motor caused the failure of a PB test [25] in the Finnish 3U-CubeSat Aalto-1 [26,27]. However, since such failures were related to technological issues and not to intrinsic problems of the PB concept, other in-orbit tests are planned to be repeated in the near future. In this regard, the Finnish Centre of Excellence in Research of Sustainable Space (FORESAIL) developed the FORESAIL-1 satellite, a 3U-CubeSat launched in May 2022 [28]. One of the FORESAIL-1 mission objectives consists of the deployment of a $40\mathrm{m}$-tether and in the implementation of a PB test at the end of the satellite operative life [29]. Moreover, the EstCube project, which was not stopped by the failure of EstCube-1 mission, has scheduled the launch of a new satellite (called EstCube-2) in 2023, with the aim of deploying a $300\mathrm{m}$-tether and performing a Coulomb drag-based deorbiting from $700\mathrm{km}$ to $500\mathrm{km}$ of altitude [28,30]. The growing interest in the Coulomb drag-based technologies is confirmed by the design of the private satellite AuroraSat-1, a 1.5U-CubeSat launched on May 2022, which is equipped with a PB tether to perform a deorbiting at the end of its operative phase. Figure 2 summarizes past, current, and future missions aimed at testing the PB technology.

The aim of this paper is to provide an approximate model capable of simulating the geocentric trajectory of a PB-based spacecraft orbiting in a LEO, by using a set of non-singular modified orbital parameters. The PB-induced force used in the simulations is estimated with an analytical mathematical model [31], which introduces some simplifying assumptions. In fact, the small magnitude of the Coulomb drag allows a perturbative approach and an asymptotic series expansion to be employed [32,33,34] in the spacecraft trajectory analysis. The resulting differential equations that describe the evolution of the osculating orbit are written in terms of a set of non-singular modified orbital parameters. These equations can be integrated, and the results can be used to analyze the spacecraft geocentric trajectory during the deorbiting phase. The accuracy of the method can be improved by periodically updating the initial conditions with a rectification procedure. In this case, the magnitude of the Coulomb drag can be adjusted at each rectification point and maintained constant between two consecutive rectifications.

This paper is structured as follows. The trajectory equations of a spacecraft with a small continuous drag acceleration, based on a perturbative approach, are first specialized to the PB case. The PB-induced drag is then estimated by means of the simplified mathematical model proposed in Ref. [31], which is here briefly summarized. Finally, the approximate model is validated by comparison with an orbital propagator for some deorbiting profiles of small satellites initially placed on a LEO, and the PB device performance is evaluated in a test scenario.

2. Trajectory Approximation

Consider a spacecraft of mass m, equipped with a PB device, which is initially tracing a LEO with orbital eccentricity $e_0$. At a given time instant $t_0\triangleq 0$, when the spacecraft altitude is $h\left(t_0\right)\triangleq h_0$ (i.e., the orbital radius is $r\left(t_0\right)\triangleq r_0=h_0+R_{\oplus }$, $R_{\oplus }$ being the Earth’s mean radius), and the true anomaly is $\nu \left(t_0\right)\triangleq {\nu }_0$, the PB conducting tether is unreeled and charged by a suitable electric voltage source. Note that if the tether polarity is negative, the spacecraft itself can act as an electric power supply, since it acquires a negative charge due to the high thermal mobility of the electrons, thus removing the necessity of a power source [21,35]. In the rest of the work, a negative voltage is accordingly assumed.

The interaction between the ionosphere and the tether generates a Coulomb drag that acts in the opposite direction with respect to the relative velocity between the spacecraft and the ions. If the thermal oscillations are neglected, the ion-spacecraft relative velocity coincides with the actual vehicle orbital velocity, and the direction of the PB-induced drag acceleration substantially coincides with the tangent to the spacecraft osculating orbit [20].

In analogy with Refs. [32,33,34], the spacecraft motion under the Earth’s gravity and the PB-induced (Coulomb) drag can be described with a set of non-singular orbital parameters $\{q_1,q_2,q_3\}$ that depends on the spacecraft osculating orbit eccentricity e and the specific angular momentum magnitude H, viz.

$$q_1\triangleq {\displaystyle \frac{e}{\tilde{H}}}cos\omega ,q_2\triangleq {\displaystyle \frac{e}{\tilde{H}}}sin\omega ,q_3\triangleq {\displaystyle \frac{1}{\tilde{H}}}$$

(1)

where $\tilde{H}\triangleq H/\sqrt{{\mu }_{\oplus }{r}_0}$ is the osculating orbit dimensionless angular momentum, and $\omega$ is the rotation angle of the eccentricity vector $\mathit{e}$, measured counterclockwise with respect to the direction of the eccentricity vector $\mathit{e}\left(t_0\right)\triangleq {\mathit{e}}_0$ at the initial time $t_0$; see Figure 3.

Now introduce an angular coordinate $\theta$, defined as the angle (measured counterclockwise) from the direction of ${\mathit{e}}_0$ to the current spacecraft position vector $\mathit{r}$, that is

$$\theta =\nu +\omega \mathrm{with}\theta \left(t_0\right)\triangleq {\theta }_0={\nu }_0$$

(2)

The spacecraft orbital radius r is given by the simple equation

$$r={\displaystyle \frac{r_0}{q_1{q}_{3}cos\theta +q_2{q}_{3}sin\theta +q_3^2}}$$

(3)

Using $\theta$ as the independent variable, the evolution of the non-singular modified orbital parameters may be written as [32,36]

$$\frac{\mathrm{d}}{\mathrm{d}\theta }\left[\begin{array}{c}{q}_1\\ q_2\\ q_3\end{array}\right]=-\frac{ϵ}{q_3{s}^3\sqrt{e^2+2ecos\nu +1}}\left[\begin{array}{cc}ssin\theta & (s+q_3)cos\theta \\ -scos\theta & (s+q_3)sin\theta \\ 0& -q_3\end{array}\right]\left[\begin{array}{c}esin\nu \\ 1+ecos\nu \end{array}\right]$$

(4)

where $\nu$ is the spacecraft true anomaly measured on the osculating orbit, and s is an auxiliary dimensionless parameter defined as

$$s\triangleq q_{1}cos\theta +q_{2}sin\theta +q_3$$

(5)

while $ϵ$ is the ratio of the Coulomb drag magnitude $D_c$ to the Earth’s gravitational attraction at $t_0$, or

$$ϵ\triangleq {\displaystyle \frac{D_c}{m{\mu }_{\oplus }/r_0^2}}$$

(6)

The parameter $ϵ$ can be thought of as the dimensionless form of the PB-induced drag acceleration, being the spacecraft mass m a constant.

The initial conditions of Equation (4) are found by calculating the values of the non-singular modified orbital parameters at $t_0$, viz.

$$q_1\left({\nu }_0\right)={\displaystyle \frac{e_0}{{\tilde{H}}_0}},q_2\left({\nu }_0\right)=0,q_3\left({\nu }_0\right)={\displaystyle \frac{1}{{\tilde{H}}_0}}$$

(7)

where the subscript 0 denotes the value on the initial orbit. The initial magnitude of the specific angular momentum may be written as

$${\tilde{H}}_0\triangleq \frac{H_0}{\sqrt{{\mu }_{\oplus }{r}_0}}\equiv \sqrt{1+e_{0}cos{\nu }_0}$$

(8)

so that Equation (7) becomes

$$q_1\left({\nu }_0\right)={\displaystyle \frac{e_0}{\sqrt{1+e_{0}cos{\nu }_0}}},q_2\left({\nu }_0\right)=0,q_3\left({\nu }_0\right)={\displaystyle \frac{1}{\sqrt{1+e_{0}cos{\nu }_0}}}$$

(9)

2.1. Asymptotic Series Expansion Approach

The analysis of the spacecraft trajectory can be simplified by assuming that the Coulomb drag magnitude is significantly smaller than the local gravitational gravity [31], that is, $ϵ\ll 1$ in Equation (4). In other terms, the PB-induced drag behaves like a perturbation acting on a Keplerian motion, in analogy with the approach highlighted in Ref. [32].

We will therefore adopt the same procedure discussed in Refs. [33,34] for either a solar sail or an E-sail spacecraft, by specializing it to the case of a tangential drag [37]. To that end, the modified orbital parameters (1) are written with a series expansion in the form

$$q_i=q_{i_0}+ϵq_{i_1}+\mathcal{O}\left({ϵ}^2\right)withi\in \{1,2,3\}$$

(10)

where the generic $q_{i_j}$ coincides with the j-th order perturbative term of $q_i$.

When Equation (10) is substituted into Equation (4) and the perturbative terms of the same order are equated, we obtain a set of differential equations that approximate the variation of the perturbative terms with respect to the angular coordinate $\theta$. In particular, the zeroth order terms are found to be

$$\frac{\mathrm{d}{q}_{i_0}}{\mathrm{d}\theta }=0withi\in \{1,2,3\}$$

(11)

which, as expected, simply state that the modified orbital parameters are constant when the trajectory is unperturbed (Keplerian). Enforcing the initial conditions, we obtain

$$q_{1_0}={\displaystyle \frac{e_0}{\sqrt{1+e_{0}cos{\nu }_0}}},q_{2_0}=0,q_{3_0}={\displaystyle \frac{1}{\sqrt{1+e_{0}cos{\nu }_0}}}$$

(12)

A more interesting result comes from the equations involving the first-order perturbative terms, which may be written as follows

$$\frac{\mathrm{d}}{\mathrm{d}\theta }\left[\begin{array}{c}{q}_{1_1}\\ q_{2_1}\\ q_{3_1}\end{array}\right]=\frac{{\tilde{H}}_0^3}{{(1+e_{0}cos\theta )}^2\sqrt{e_0^2+2e_{0}cos\theta +1}}\left[\begin{array}{c}{e}_0+2cos\theta \\ 2sin\theta \\ -1\end{array}\right]$$

(13)

where ${\tilde{H}}_0$ is given by Equation (8). Assuming a constant magnitude of the Coulomb drag, and so a constant value of $ϵ$, Equation (13) can be integrated in the angular coordinate $\theta$ to obtain the evolution of the generic $q_{i_1}$ term. More precisely, with the aid of Equation (1), the modified orbital parameters can be written as

$$\begin{array}{cc}\hfill q_1& =-2ϵ\left(sin\theta -sin{\nu }_0\right)\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill q_2& =2ϵ\left(cos\theta -cos{\nu }_0\right)\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill q_3& ={\displaystyle \frac{1}{{\tilde{H}}_0}}+ϵ\left(\theta -{\nu }_0\right)\hfill \end{array}$$

(16)

which are valid in case of initial circular orbit ($e_0=0$). If, instead, the initial orbit is elliptic ($e_0\ne 0$), then

$$\begin{array}{cc}\hfill q_1& ={\displaystyle \frac{e_0}{{\tilde{H}}_0}}-ϵ{\displaystyle \frac{{\tilde{H}}_0^3}{{(1-e_0^2)}^2}}\left[Q_1(e_0,E)-Q_1(e_0,E_0)\right]\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill q_2& =-ϵ{\displaystyle \frac{{\tilde{H}}_0^3}{{(1-e_0^2)}^{3/2}}}\left[Q_2(e_0,E)-Q_2(e_0,E_0)\right]\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill q_3& ={\displaystyle \frac{1}{{\tilde{H}}_0}}-ϵ{\displaystyle \frac{{\tilde{H}}_0^3}{{(1-e_0^2)}^2}}\left[Q_3(e_0,E)-Q_3(e_0,E_0)\right]\hfill \end{array}$$

(19)

Note that in Equations (17)–(19) the angle E depends on the angular coordinate $\theta$ through the equation

$$E=E(e_0,\theta )\triangleq 2arctan\left[\sqrt{{\displaystyle \frac{1-e_0}{1+e_0}}}tan\left({\displaystyle \frac{\theta }{2}}\right)\right]$$

(20)

with

$$E_0=2arctan\left[\sqrt{{\displaystyle \frac{1-e_0}{1+e_0}}}tan\left({\displaystyle \frac{{\nu }_0}{2}}\right)\right]$$

(21)

Finally, the expressions of $Q_i$ in Equations (17)–(19) are [32]

$$\begin{array}{c}{Q}_1(e_0,E)={\displaystyle \frac{2C_1\left(e_0\right)(2-e_0^2)-4C_2\left(e_0\right)}{\pi e_0}}E+\left({\displaystyle \frac{15e_0^4}{32}}+{\displaystyle \frac{3e_0^2}{4}}+2\right)sin\left(E\right)+\hfill \\ \hfill -\left({\displaystyle \frac{e_0^3}{8}}+{\displaystyle \frac{e_0}{2}}\right)sin\left(2E\right)+\left({\displaystyle \frac{5e_0^4}{64}}+{\displaystyle \frac{e_0^2}{12}}\right)sin\left(3E\right)-\left({\displaystyle \frac{e_0^3}{32}}\right)sin\left(4E\right)+\left({\displaystyle \frac{3e_0^4}{320}}\right)sin\left(5E\right)\end{array}$$

(22)

$$\begin{array}{c}{Q}_2(e_0,E)=-\left({\displaystyle \frac{3e_0^4}{32}}+{\displaystyle \frac{e_0^2}{4}}+2\right)cos\left(E\right)+\left({\displaystyle \frac{e_0^3}{8}}+{\displaystyle \frac{e_0}{2}}\right)cos\left(2E\right)+\hfill \\ \hfill -\left({\displaystyle \frac{3e_0^4}{64}}+{\displaystyle \frac{e_0^2}{12}}\right)cos\left(3E\right)+\left({\displaystyle \frac{e_0^3}{32}}\right)cos\left(4E\right)-\left({\displaystyle \frac{3e_0^4}{320}}\right)cos\left(5E\right)\end{array}$$

(23)

$$\begin{array}{c}{Q}_3(e_0,E)={\displaystyle \frac{2C_1\left(e_0\right)-4C_2\left(e_0\right)}{\pi }}E+\left({\displaystyle \frac{3e_0^3}{32}}+2e_0\right)sin\left(E\right)+\hfill \\ \hfill -\left({\displaystyle \frac{3e_0^2}{8}}+{\displaystyle \frac{7e_0^4}{32}}\right)sin\left(2E\right)+\left({\displaystyle \frac{e_0^3}{12}}\right)sin\left(3E\right)+\left(-{\displaystyle \frac{7e_0^4}{256}}\right)sin\left(4E\right)\end{array}$$

(24)

where $C_1\left(x\right)$ and $C_2\left(x\right)$ denote the complete elliptic integrals of the first and second kind of the argument x, respectively.

2.2. Rectification Procedure

The error introduced by the asymptotic series approximation increases when either the “perturbative” parameter $ϵ$ or the angular coordinate $\theta$ increases. Here, the term error is intended as the difference between the approximate value and the numerical solution of the spacecraft equation of motion. Moreover, since the plasma properties in the Earth’s upper atmosphere vary with altitude, the assumption of a constant Coulomb drag throughout the whole deorbiting profile is not realistic [22,31]. To overcome these issues, a rectification procedure can be easily added to the mathematical model. The necessary operative steps are thoroughly discussed in Refs. [32,33] and are here briefly summarized.

First, the angular coordinate $\theta ={\theta }_r$ at which the rectification takes place is selected, and the corresponding osculating orbital parameters are calculated as

$$e_r={\displaystyle \frac{\sqrt{q_{1_r}^2+q_{2_r}^2}}{q_{3_r}}},{\tilde{H}}_r={\displaystyle \frac{1}{q_{3_r}}},{\omega }_r=arctan\left({\displaystyle \frac{q_{2_r}}{q_{1_r}}}\right),{\nu }_r={\theta }_r-{\omega }_r$$

(25)

where $q_{i_r}=q_i\left({\theta }_r\right)$. Bearing in mind the expression of r from Equation (3), the spacecraft altitude $h\left({\theta }_r\right)\triangleq h_r$ at the rectification point is calculated from

$$h_r=r_r-R_{\oplus }\equiv {\displaystyle \frac{r_0}{q_{1_r}{q}_{3_r}cos{\theta }_r+q_{2_r}{q}_{3_r}sin{\theta }_r+q_{3_r}^2}}-R_{\oplus }$$

(26)

in such a way that a new set of auxiliary variables may be defined for $\theta \ge {\theta }_r$, that is

$$\overline{\theta }\triangleq \theta -{\omega }_r,\overline{\omega }\triangleq \omega -{\omega }_r,{\overline{q}}_1\triangleq {\displaystyle \frac{e}{\tilde{H}}}cos\overline{\omega },{\overline{q}}_2\triangleq {\displaystyle \frac{e}{\tilde{H}}}sin\overline{\omega },{\overline{q}}_3\triangleq {\displaystyle \frac{1}{\tilde{H}}}$$

(27)

whose initial conditions at the rectification point are

$${\overline{\theta }}_0={\nu }_r,{\overline{\omega }}_0=0,{\overline{q}}_{1_0}={\displaystyle \frac{e_r}{{\tilde{H}}_r}},{\overline{q}}_{2_0}=0,{\overline{q}}_{3_0}={\displaystyle \frac{1}{{\tilde{H}}_r}}$$

(28)

with $e_r$, ${\omega }_r$, and ${\tilde{H}}_r$ given by Equation (25).

The definitions of Equation (27) generate a new set of initial conditions for the auxiliary parameters, given by Equation (28), that are analogous to the ones expressed by Equation (12). Moreover, when the rectification procedure is applied, a new value of the dimensionless propulsive parameter $\overline{ϵ}$ is calculated, as is discussed in the next section. The new value of the dimensionless propulsive parameter is used for the trajectory determination just after the rectification point, that is, for $\theta \ge {\theta }_r$. In other terms, the geocentric trajectory of the PB-based spacecraft is obtained by considering a sort of piecewise-constant [38] Coulomb drag acceleration.

Given all the previous considerations, the spacecraft geocentric trajectory for $\theta \ge {\theta }_r$ can be estimated from Equations (14)–(16) if $e_r=0$, or from Equations (17)–(19) if $e_r\ne 0$, provided that $e_0$, ${\tilde{H}}_0$, ${\nu }_0$, $ϵ$, and $\theta$ are substituted with $e_r$, ${\tilde{H}}_r$, ${\nu }_r$, $\overline{ϵ}$, and $\overline{\theta }$, respectively. Finally, the last step is to apply a rotation matrix to the auxiliary variables in order to obtain the original modified orbital parameters for $\theta \ge {\theta }_r$, viz.

$$\left[\begin{array}{c}{q}_1\left(\theta \right)\\ q_2\left(\theta \right)\\ q_3\left(\theta \right)\end{array}\right]=\left[\begin{array}{ccc}cos{\omega }_r& -sin{\omega }_r& 0\\ sin{\omega }_r& cos{\omega }_r& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{\overline{q}}_1\left(\overline{\theta }\right)\\ {\overline{q}}_2\left(\overline{\theta }\right)\\ {\overline{q}}_3\left(\overline{\theta }\right)\end{array}\right]for\theta \ge {\theta }_r$$

(29)

where ${\omega }_r$ is given by the third of Equation (25).

The rectification procedure described here is illustrated in the block scheme of Figure 4. The same procedure can be repeated during the numerical simulation to obtain a more accurate spacecraft (deorbiting) trajectory with a more realistic estimation of the PB-induced drag, as described in the next section. Even if a large number of rectifications tends to increase the computation time, the latter remains two orders of magnitude smaller than that required by a numerical integration of the nonlinear equations of motion.

2.3. PB-Induced Drag Model

The mathematical method used for the trajectory analysis must be completed with a thrust model for the PB-induced drag. According to Ref. [35], the magnitude of the Coulomb drag $D_c$ generated by a PB conducting tether can be written as

$$D_c=3.864L_t{m}_{i}nv\sqrt{{\displaystyle \frac{{ϵ}_0{V}_a}{eln}}}exp\left(-{\displaystyle \frac{m_i{v}^2}{2el{V}_a}}\right)$$

(30)

where $L_t$ is the tether length, $m_i$ is the ions mass, n is the plasma bulk number density, v is the relative velocity of the ions relative to the spacecraft, $el$ is the elementary charge, and ${ϵ}_0$ is the vacuum permittivity. In Equation (30), the auxiliary voltage $V_a$ is defined as

$$V_a\triangleq {\displaystyle \frac{2|V_t|}{ln\left({\displaystyle \frac{{ϵ}_0\left|V_t\right|}{el{n}_0{b}_t{r}_w}}\right)}}$$

(31)

where $V_t$ is the tether voltage, which is assumed negative according to Refs. [21,35].

Note that Equation (30) is based on the following simplifying assumptions: (i) the PB-induced drag force per unit of tether length is constant, (ii) the effects of the magnetic field are neglected [22], and (iii) the PB conducting tether has a Heytether structure [39]. As sketched in Figure 5, a Heytether consists of a principal line with radius $r_w$ and multiple interconnections introduced to increase the tether resistance against possible micrometeoroid impacts. The whole arrangement makes the total conducting tether width equal to $b_t$.

In analogy with Ref. [31], the Coulomb drag magnitude $D_c$ at a generic spacecraft altitude can be written as a function of the value at a reference altitude as

$$D_c\simeq D_{c_0}\sqrt{{\displaystyle \frac{n}{n_0}}}{\left({\displaystyle \frac{v}{v_0}}\right)}^2{f}_1(v,v_0)f_2(n,n_0,v,v_0)f_3(n,n_0)$$

(32)

where the subscript 0 indicates that the initial altitude is taken as a reference value. The explicit expressions of the auxiliary functions $f_i$ are reported in Ref. [31] and are here omitted for the sake of conciseness. Equation (32) introduces a further assumption; that is, the ions’ mass is constant and equal to the atomic mass of the dominating species in LEO, the atomic oxygen; that is, $m_i\simeq 16\mathrm{u}$.

The work by Orsini et al. [31] also makes an analysis of the order of magnitudes of the terms involved in Equation (32), from which it is shown that a conservative approximation in the LEO range is to assume $f_1{f}_2{f}_3\simeq 1$. When the relative velocity v is approximated with the local circular velocity (thus neglecting the thermal oscillations of ions and the orbital eccentricity), it is possible to assume that the variation of the velocity modulus in the LEO range is small; that is, ${(v/v_0)}^2\simeq 1$. Enforcing these approximations into Equation (32) yields

$$D_c\simeq D_{c_0}\sqrt{{\displaystyle \frac{n}{n_0}}}$$

(33)

Assuming that the plasma temperature T is constant and estimating the plasma density with the geopotential ionosphere model, Equation (33) can be rewritten as

$$D_c\simeq D_{c_0}exp\{-\frac{m_i{\mu }_{\oplus }}{4k_{B}T}\left[\frac{h}{{(R_{\oplus }+h)}^2}-\frac{h_0}{{(R_{\oplus }+h_0)}^2}\right]\}$$

(34)

where h is the spacecraft altitude, and $k_B$ is the Boltzmann constant. The latter equation gives the variation of the PB-induced force with respect to the reference value $D_{c_0}$ as a function of the spacecraft altitude, which can be calculated at every rectification point by means of Equation (26). Recall that $D_{c_0}$ may be calculated with Equation (30). Figure 6 shows the variation of $D_c$ as a function of the spacecraft altitude under the assumption of a mean solar activity. The value of $D_c$ calculated at each rectification point can be substituted into Equation (6) to update the “perturbative” parameter $ϵ$ and simulate a piecewise constant Coulomb drag acceleration. The procedure is illustrated in Figure 4.

3. Case Study

The approximate method is used to estimate the deorbiting profile of a set of small satellites, initially placed in a typical LEO. The initial (reference) altitude is set equal to $h_0=1000\mathrm{km}$, the middle of the LEO orbital range [1], with a corresponding local circular velocity $v_0=\sqrt{{\mu }_{\oplus }/(R_{\oplus }+h_0)}$. A mean solar activity is considered, which corresponds to a plasma temperature of $T=1011.5\mathrm{K}$ in the upper stages of the Earth’s atmosphere. The ion number density at the reference altitude is estimated as $n_0=3\times {10}^{10}{\mathrm{m}}^{-3}$, which is in good accordance with the international reference ionosphere (IRI) [40,41] for a mean solar activity. The spacecraft parameters used in the simulations have been selected in analogy with the actual values of the PB-equipped missions and in accordance with the current or near-term technology level for nanosatellites, as reported in

Table 1. Characteristics of the PB-based small satellites in the deorbiting simulations.

Configuration	m [kg]	${\mathit{L}}_{\mathit{t}}$ [m]	$|{\mathit{V}}_{\mathit{t}}|$ [V]
spacecraft ①	1	25	500
spacecraft ②	4	100	1000
spacecraft ③	10	300	1000

.

The variations of the orbital mean altitude are presented in Figure 7, where the outputs of the approximate method (with 100 rectifications per year) are compared with those from an orbital propagator that numerically integrates the equations of motion and updates the PB-induced drag acceleration magnitude at each time step. In this comparison, the spacecraft nonlinear equations of motion have been numerically integrated with a variable order Adams–Bashforth–Moulton PECE solver [42,43], with absolute and relative errors equal to ${10}^{-10}$.

The results from the two methods are so similar that the different lines are nearly coincident, although the computational cost of the approximation is about two orders of magnitude smaller than the numerical integration. Note that the decay rate is slower at higher altitudes when the plasma density is smaller and is steeper in the final part. Since the atmospheric drag is neglected in this analysis and its effects are significant for $h\le 500\mathrm{km}$, it is reasonable to conclude that this effect should be even more evident in real cases.

Finally,

Table 2. Decay time from $1000\mathrm{km}$ to $300\mathrm{km}$ estimated with an orbital propagator and the approximate method. A mean solar activity is assumed.

Configuration	Numerical	Approximate	Percentage Error
spacecraft ①	$1317\mathrm{days}$	$1320\mathrm{days}$	$0.26\%$
spacecraft ②	$924\mathrm{days}$	$928\mathrm{days}$	$0.38\%$
spacecraft ③	$770\mathrm{days}$	$774\mathrm{days}$	$0.45\%$

shows the times required for the spacecraft to lower its altitude below $300\mathrm{km}$. Note that the decay time is estimated by the approximate method with high accuracy when compared with the numerical results. It is evident that the PB is capable of deorbiting a nanosatellite in a LEO range and to comply with the international guidelines, even if the contribution of the atmospheric drag is (conservatively) neglected. Unluckily, the procedure discussed here is not capable of providing a quick estimate of the required PB characteristics, given the desired decay time. However, the interested reader may refer to Ref. [44], where an estimate of the deorbiting time of a PB-equipped spacecraft is obtained with a different approach.

4. Conclusions

An approximate method for the trajectory analysis of a spacecraft deorbiting from a LEO is presented. The analytical approximation uses a perturbative approach to estimate the evolution of the spacecraft osculating orbit parameters based on the small magnitude of the Coulomb drag acceleration. The latter is evaluated by means of a simplified mathematical model that relates the drag force magnitude to the spacecraft altitude. A number of rectifications are included in the proposed model, and the value of the PB-induced drag acceleration is updated at each rectification point and maintained constant between two successive rectifications.

The proposed procedure allows a very accurate estimation of the deorbiting profile to be obtained along with the decay time of a small spacecraft orbiting in LEO. A preliminary estimation of the PB performance is also possible in an efficient way, as the computation time required by the proposed approach is about two orders of magnitude smaller than those associated to a numerical integration of the equations of motion with an orbital propagator.