2.1 Aerodynamics model

Due to the high degree of rarefaction at high altitudes, the kinetic particle-based DSMC method is highly suitable for the calculation of the aerodynamics forces and torques on an orbiting spacecraft. This method proposed by Bird [Reference Bird39] provides a physics-based solution to the Boltzmann equation by mimicking the behaviour of molecules but without directly solving the Boltzmann equation. The DSMC method uses simulated molecules, each of which represents a cluster of actual molecules. In the DSMC approach, the behaviour of the molecules is modelled through their movement and intermolecular collisions. Their movement is modelled in a deterministic manner, whereas intermolecular collisions are modelled probabilistically. The primary assumption of the DSMC approach is that the movement and intermolecular collisions can be decoupled over a small time interval. DSMC simulation results have also been proven to converge towards solutions of the Boltzmann equation in the limit of an infinite number of particles, for a cell size smaller than the mean free path and a time step smaller than the mean collision time [Reference Bird39]. A proof of the convergence of the DSMC method and the Gas-Kinetic Unified Algorithm (GKUA) for the Boltzmann equation was investigated in the work of Li et al. [Reference Li, Fang, Jiang and Wu40]. The DSMC method can be used in all flow regimes in the dilute gas framework, but is primarily used in the rarefied and free molecular regime (Kn $ > 0.1$ ) due to its computational efficiency. Moreover, the development and application of coupled Navier–Stokes/DSMC solvers to simulate hypersonic flows in the transitional flow regime and multi-component plume flows from satellite attitude-control engines have been extensively reviewed in the work of Li et al. [Reference Li, Li, Wu and Peng41]. Furthermore, Li et al. [Reference Li, Cai and Li42] studied low-speed microchannel rarefied flows using the direct simulation Bhatnagar–Gross–Krook (DSBGK) model equation, also discussing the limitations associated with the use of the DSMC method to study such problems. A detailed explanation of the DSMC method can be found in Ref. [Reference Bird39].

For this work, the parallel, in-house, multi-species DSMC solver named Non-equilibrium Flow Solver (NFS) [Reference Kumar and Chinnappan43], which is based on Bird’s DSMC algorithm, is modified to develop a computationally efficient DSMC-based solver for free molecular flow [Reference Chinnappan, Kumar, Arghode and Myong44]. The modifications made to the DSMC algorithm are consistent with the assumption of free molecular flow. First, intermolecular collisions are not modelled in the modified DSMC solver, since the gas flow is assumed to be collisionless in the free molecular flow regime. Second, after the gas–surface interaction model calculations are performed for the interaction of a gas particle with a surface, its velocity is then assigned the same value as before the collision and it is relocated randomly inside the domain. This is to ensure that the incident gas remains unaffected by the presence of body, which is consistent with the assumptions made in free molecular flow theory [Reference Chambre and Schaaf45]. Thus, the core of the modified DSMC solver consists of initialisation, particle movement, gas–surface interaction and sampling of surface properties. A detailed explanation related to the algorithm of the solver can be found in Chinnappan et al. [Reference Chinnappan, Kumar, Arghode and Myong44].

2.2 Calculation of drag and torque

Aerodynamic forces are the main non-gravitational forces acting on a satellite due to the presence of the rarefied atmosphere in LEO. These forces can be resolved as follows (Fig. 1):

The drag force is a non-conservative force [Reference Vallado25] acting in the direction opposite to the satellite’s velocity vector, thereby resulting in energy dissipation from the orbit of the satellite. The dissipation of energy slows down the spacecraft, which then leads to orbital decay. The drag and lift forces are calculated in the DSMC framework by integrating the pressure and shear stress distribution around the body. The lift force is another part of the aerodynamic force, acting in a direction normal to the relative velocity vector. It is worth noting that there is no resultant lift acting on a satellite exhibiting uncontrolled tumbling motion, which is a valid assumption since the effect of the lift force is averaged out over its lifetime [Reference Cook46]. Thus, in the present work on the orbital decay analysis of a tumbling body, only the perturbation due to the drag force and the effect of the tumbling motion are considered. Moreover, it is shown below that the magnitude of the lift force is much smaller as compared with the drag force.

A spacecraft undergoes tumbling motion throughout its orbit [Reference De Pontieu47]. In the present work, the influence of the aerodynamic forces and torques on the tumbling motion is examined. The in-house modified DSMC solver, as explained above, is used to calculate the aerodynamic properties of the spacecraft. The aerodynamic forces (F) are obtained by integrating the pressure and shear stress over the surface of the spacecraft. These aerodynamic forces act around the Centre of Mass (CM), producing a torque. The torque caused by the aerodynamic forces can bes calculated using the analytical expression $\vec{T}= \delta \vec{d}\times \vec{F}$ , where $\delta\vec{d}$ is the position vector of the centre of pressure (CP) from the CM, and $\vec{F}$ is the resultant aerodynamic force acting on the body, as shown in Fig. 2. The lift (L) and drag (D) forces are the x and y components of the resultant force vector, respectively. To simplify this problem, we consider only one-dimensional rotation of a spacecraft about its pitching axis (z-axis) to model the tumbling motion, as illustrated in Fig. 2.

Figure 1. The aerodynamic forces on a satellite [Reference De-Lafontaine and Garg15].

Figure 2. Schematic of a body experiencing torque.

A parameter called the pitching torque coefficient ( $C_T$ ) is defined to normalise the torque acting on the spacecraft:

(1) \begin{equation}C_T=\frac{T_z}{\rho U^{2}SL},\end{equation}

where $T_z$ is the pitching torque about the z-axis, S is the reference area of the spacecraft and L is its characteristic length.

2.3 Orbital dynamics model

Using Newton’s second law and the universal law of gravitation, the equations for the relative motion of two bodies can be derived by assuming them to be point masses:

(2) \begin{align}\ddot{\vec{r}}+\frac{\mu}{r^{3}}\vec{r}=0,\end{align}

where $\vec{r}$ is the position vector of the spacecraft measured from the centre of the attracting body and $\mu$ is the standard gravitational parameter. Considering Earth as the attracting body, one has $\mu = 3.99\times10^{14}\textit{m}^{3}/\mathrm{s}^{2}$ [Reference Curtis48].

Equation (2) is the basis of the unperturbed motion of a spacecraft, considering only the gravitational force. In real life, a spacecraft in LEO encounters various forces apart from the gravitational force. These various forces cause perturbations that lead to a deviation of the spacecraft from its ideal orbit. The equation of relative motion of two bodies under such orbital perturbation can be expressed as

(3) \begin{equation}\ddot{\vec{r}}+\frac{\mu}{r^{3}} \vec{r}=\vec{a},\end{equation}

where $\vec{a}$ is the sum of all the perturbing accelerations [Reference Curtis48].

In the current study, only the drag perturbation is considered, because for spacecraft in LEO, this perturbation is many times larger than other orbital perturbations [Reference King-Hele49]. The mathematical expression used to model the deceleration due to drag is given by Ref. [Reference Vallado25] as follows:

(4) \begin{equation}\vec{a}_{d}=-\frac{1}{2}\rho \frac{C_{D} A}{m} v_{R}^{2} {\hat{v}_{R }},\end{equation}

where $C_D$ is the drag coefficient, $\rho$ is the density of the atmosphere, A is the cross-sectional area and m is the mass of the satellite, $v_R$ is the magnitude of the satellite’s velocity relative to the atmosphere and $\hat{v}_{R }$ is the unit vector along the satellite’s relative velocity. The term $\frac{m}{C_D A}$ is known as the ballistic coefficient $\beta$ , being a measure of the spacecraft’s susceptibility to drag effects and indicating how fast the orbit of the spacecraft will decay. Spacecraft with low $\beta$ values respond quickly to the atmosphere, and their orbit will decay faster than those with high values of $\beta$ [Reference Wertz, Everett and Puschell50].

The orbital decay program is employed to estimate the orbital decay trajectory of a tumbling body. The basis of the orbital decay is Equations (3) and (5) (as given below). Equation (3) describes the relative motion of two bodies under orbital perturbations. In the present work, the drag force is the only orbital perturbation considered in the orbital decay analysis. The drag perturbation is modelled using Equation (4). Meanwhile, Equation (5) describes the rotational dynamics of the body and is used to model its tumbling motion, considering only rotation about its z-axis. In Equation (5), $\omega$ is the angular momentum and $I_z$ is the moment of inertia of the spacecraft about the z-axis.

(5) \begin{align}T_z=I_z\ddot{\alpha}=I_z\dot{\omega}.\end{align}

The computer program for the orbital decay simulation solves Equations (3) and (5) simultaneously. The Earth-Centred Inertial (ECI) coordinate system is employed in the orbital decay simulation. A circular orbit is assumed while developing the computer code. Note that this assumption is made to simplify the analysis to some extent. Figure 3 shows the algorithm of the orbital decay program.

Figure 3. Orbital decay program algorithm.

In Fig. 3, h and $h_{\rm Re}$ are the altitude and re-entry altitude of the spacecraft, respectively, from the Earth’s surface. Reentry is assumed to occur when the spacecraft descends to an altitude of 180km [Reference Kennewell51]. After reaching this altitude, the spacecraft will impact within a few hours. The look-up table for the aerodynamics coefficients obtained using the free molecular flow simulation and the US Standard Atmosphere 1976 [52] is used as input to this computer program. Finally, the orbital decay trajectory of the spacecraft is obtained by plotting the altitude of the spacecraft versus time.

3.1 Aerodynamics model validation

The modified DSMC solver is already validated for many flow problems, such as flow over a sphere at low and high speed ratios, and a flat plate at different angles of attack [Reference Chinnappan, Kumar, Arghode and Myong44]. Additionally, we present herein its validation for free molecular flow over a sphere by comparing the simulation results against theory. The validation case study involves flow over a spherical body at different speed ratios. The drag coefficient ( $C_D$ ) is compared with the theoretical expression for the drag experienced by a sphere (with diffuse reflection) in the free molecular flow regime given by Equation (6). Figure 4 shows the computational domain and the mesh employed for the validation case study.

(6) \begin{align}C_{D,{\rm diffuse}}&=\frac{e^{-S^{2}}}{\sqrt{\pi} S^{3}}\left(1+2 S^{2}\right)+\frac{4 S^{4}+4 S^{2}-1}{2 S^{4}} \operatorname{erf}(S)+\frac{2 \sqrt{\pi}}{3 S_{{w}}}\end{align}

Figure 4. (a) Computational domain. (b) Mesh for a sphere experiencing free molecular flow.

A cubic computational cell is considered with a spherical body at the centre. Periodic boundary conditions are employed on all six sides of the domain. Each simulation is performed with 20,000 simulated particles and 50,000 samples. The time step for the simulation is $1\times10^{-10}$ s, chosen such that it is less than the mean collision time of the gas. The Cercignani–Lampis–Lord (CLL) model is used for the gas–surface interaction with entirely diffuse reflection.

The freestream conditions used for the simulation are presented in Table 1. The wall temperature is kept the same as the freestream temperature.

Table 1. Freestream conditions for free molecular flow simulations over a sphere

The variation of the drag coefficient with the speed ratio is obtained from the in-house modified DSMC solver, as described in Section A. The results are plotted and compared with theory in Fig. 5. It is clearly observed that the $C_D$ profile is in excellent agreement with the analytical results over the considered range of speed ratio.

Figure 5. Comparison of DSMC and analytical results for spherical body.

3.2 Orbital decay validation

The orbital decay simulation is validated against a period reduction model given by the Australian IPS Radio and Space Services [Reference Kennewell51]. The governing equation for the period reduction model is expressed as

(7) \begin{align}\frac{dp}{dt}=-3\pi a\rho \left(\frac{A_e}{m}\right),\end{align}

where P is the orbital time period, m is the mass of the satellite, a is the semi-major axis of its orbit, $A_e$ is the effective cross-sectional area, defined as $A_e= A C_D$ with A being the cross-sectional area of the satellite and $C_D$ is the coefficient of drag experienced. In the period reduction model, drag is the only aerodynamics force taken into account. Equation (7) calculates the reduction in the orbital period in the time step dt by using the current altitude, the atmospheric density and the effective cross-sectional area of the satellite divided by its mass. The period of the orbit is updated, and the new semi-major axis (and hence altitude) of the orbit is calculated as $a={\left(\frac{p^{2}\mu}{4\pi^{2}}\right)}^{1/3} $ , then the process is repeated until the spacecraft reaches the re-entry altitude of 180km.

The orbital decay simulation is performed on the same prototypical satellite analysed using the period reduction model in Ref. [Reference Kennewell51]. The prototypical satellite has an effective cross-section $A_e = 1\mathrm{m}^{2}$ and mass m = 100kg. The US Standard Atmosphere 1976 [52] is applied to acquire the air density at a given altitude. A starting altitude of 200km is chosen for both simulations.

Figure 6 compares the results obtained from our orbital decay simulation with those of the period reduction model employed on the prototypical satellite. Evidently, the orbital decay trajectory obtained from the orbital decay simulation is in good agreement with the results of the periodic reduction model.

Figure 6. Comparison between orbital decay simulation and periodic reduction model.

Figure 7. Apollo space capsule geometry [Reference Griffith and Boylan53].

4.1 Orbital decay analysis

Three test cases are considered in the orbital decay analysis. In the first test case, the orbital decay simulation is performed on a sphere and the tumbling Apollo space capsule to analyse their orbital trajectory. In the second test case, simulations on both bodies are performed, assuming a constant value of $C_D=2.2$ . This value of $C_D$ is chosen because it is commonly used for orbital decay analysis, based on the studies carried out by Cook [Reference Cook54]. In the third test case, the simulation is performed on both bodies but assuming the same reference area and weight for the sphere as for the Apollo space capsule, to reveal the effect of shape on the orbital decay.

The input parameters for the simulations are presented in Table 2. These input parameters are the same for all the test cases, except for $\alpha$ and $\omega$ , which are not required for the second test case since it does not consider any tumbling motion. The initial position of the spacecraft is 200km from the Earth’s surface. It is assumed that the spacecraft is initially at a zero angle-of-attack, which is not a binding constraint but rather simply an initial condition that is assumed in the present analysis. The initial angular acceleration is chosen as $24^{\circ}/\mathrm{s}$ based on observations of orbiting spacecraft [Reference De Pontieu47].

Table 2. Input parameters for orbital decay simulation

First, a comparison of the orbital decay trajectory of the sphere and Apollo space capsule is discussed to investigate the effect of the geometrical parameters on the trajectory. Next, the effect of assuming a constant $C_D$ on the trajectories of both bodies is explored, followed by a comparison of the trajectory of the Apollo space capsule with that of a sphere having the same mass and reference area to analyse the effect of shape on the trajectory. These three case studies are discussed below.

Case study I: comparison of trajectories of sphere and Apollo space capsule

A comparison of the orbital decay trajectories of the two bodies (the sphere and the Apollo space capsule) is shown in Fig. 11 as the variation of altitude with time. Figure 12 shows the drag experienced by the Apollo space capsule and the sphere respectively during their orbital decay. The drag force experienced by the Apollo space capsule is much larger than that experienced by the sphere due to its significantly larger size. On the other hand, the magnitude of the drag deceleration experienced by the Apollo space capsule is comparatively smaller than that experienced by the sphere, as shown in Fig. 13, since the mass of the Apollo space capsule is much larger than that of the sphere. Consequently, the Apollo space capsule has a larger value of the ballistic coefficient ( $\beta$ ) compared with the sphere. Therefore, de-orbiting occurs much faster in the case of the sphere.

Figure 11. Orbital decay trajectory of Apollo space capsule and sphere.

Figure 12. Drag over time experience by Apollo space capsule and sphere. (a) Apollo space capsule. (b) Sphere.

Figure 13. Drag deceleration over time experience by Apollo space capsule and sphere. (a) Apollo space capsule. (b) Sphere.

The fluctuations in the drag and the drag deceleration for the Apollo space capsule, as shown in Figs 12(a) and 13(a), can be explained by the rapid rotation of the body due to the tumbling motion. The variation of the angle-of-attack, $\alpha$ , as a result of the tumbling motion is shown in Fig. 14, as the orbit decays with time. The zoomed-in view shown alongside provides further confirmation of the same. The tumbling motion is not applicable for the case of the sphere.

Figure 14. Time variation of the angle-of-attack due to tumbling motion as the orbit decays for Apollo space capsule. A zoomed-in view is shown alongside. (a) Apollo space capsule. (b) Sphere.

Figure 15. Comparison of trajectories for Apollo space capsule and sphere with a constant drag coefficient of 2.2. (a) Apollo space capsule. (b) Sphere.

It is observed from Fig. 12 that the drag force of both bodies increases over time, which leads to an increase in the rate of loss of altitude as the spacecraft descends. This accounts for the increased rate of de-orbiting of both the Apollo space capsule and the sphere as they descend, as is evident from Fig. 11.

The orbital decay of the sphere from a circular orbit at 200km to an altitude of 180km will occur within a day, whereas for a larger spacecraft such as the Apollo space capsule, this will take about twice as long. The chances of surviving re-entry for smaller spacecraft, such as the sphere, are very slim due to aerodynamic heating. As a result, ground impacts of such spacecraft pose a lower risk. Larger spacecraft, such as the Apollo space capsule, have much higher chances of surviving re-entry, thereby posing a risk of ground impact in populated areas. Therefore, controlled re-entry of such spacecraft is necessary to avoid the potential for damage and hazard.

Case study II: effect of drag coefficient on trajectory

Figure 15 shows the orbital decay trajectory of the sphere and the Apollo space capsule for a constant $C_D$ of 2.2 and compares the results with those obtained when $C_{D}$ varies. The variation of $C_D$ with the speed ratio is considered for the sphere, whereas the variation of $C_D$ with the speed ratio and the angle-of-attack is considered for the Apollo space capsule. It is observed that, for both objects, the calculated orbital decay times for constant $C_D$ deviate from those estimated in the previous section with varying $C_D$ . Figure 16 shows that, in the case of the sphere, the deviation is comparatively smaller than that for the Apollo space capsule. This can be attributed to the fact that the $C_D$ value obtained for the sphere is close to 2.2, the value assumed for the case of constant $C_{D}$ . However, when the $C_D$ values obtained for the Apollo space capsule are analysed, they are found to deviate widely from the 2.2 level. Moreover, the $C_D$ value of the sphere changes only with the speed ratio but not with the angle-of-attack, while the $C_D$ of the Apollo space capsule changes with both the speed ratio and the angle-of-attack. This is clearly evident from Fig. 16. The deviation in the time required for orbital decay indicates the importance of using accurate estimates of $C_D$ when calculating decay trajectories.

Figure 16. Variation of $C_D$ with angle-of-attack as the body tumbles at a speed ratio of 10. The constant value of $C_D=2.2$ is shown for reference.

Figure 17. Comparison of trajectories for sphere and Apollo space capsule having same mass and reference area.

Case study III: effect of geometrical shape on trajectory

This case study is undertaken to examine the effect of geometrical shape on the trajectory of tumbling bodies. For this purpose, the reference area and mass of the sphere are kept equal to those of the Apollo space capsule. This allows us to compare the decay trajectories for both bodies having different shapes.

Figure 17 shows the trajectory of the sphere and Apollo space capsule. The current value of $C_D$ for both bodies is the same as that in case study I discussed above. It is found that the decay time of the sphere in this case is greatly increased as compared with a regular sphere as discussed in case study I, because of its increased mass and size, which are kept the same as those of the Apollo space capsule. The difference in the orbital decay time between the sphere and Apollo space capsule with the same size and mass is shown in Fig. 17. This occurs because of the difference in the value of $C_D$ , which is a consequence of the shape of the objects. As seen in case study I, the sphere has a larger $C_D$ value as compared with the Apollo space capsule, thus causing a larger drag perturbation to the sphere, as illustrated in Fig. 18. Therefore, the orbital decay rate is faster for the sphere of the same size and mass as the Apollo space capsule.

Figure 18. Drag deceleration over time experienced by Apollo space capsule and sphere having same mass and reference area. (a) Apollo space capsule. (b) Sphere.