#### Abstract

The growing risk of falling debris from outer space as well as the atmospheric interaction effect makes the orbit decay prediction of large spacecraft in very low earth orbit (VLEO) increasingly significant. Focusing on the aerodynamic perturbation effects under multiscale and nonequilibrium states on the orbit decay of the large spacecraft in VLEO at the end of its lifetime, we developed a novel perturbation prediction model covering the entire altitude range before reentry to perform long-term and short-term predictions of the large-scale spacecraft. A unified local rapid engineering algorithm for aerodynamic force and moment coefficients covering all flow regimes is proposed. The orbit perturbation models, combining the components of aerodynamics solved by the engineering algorithm, are built for the large-scale spacecraft. For altitudes ranging from 350 km to 250 km, which we defined as the slow descending stage (SDS), the two-line orbital elements (TLEs) and simplified general perturbation 4 (SGP4) model were used for long-term prediction, whereas for altitudes from 250 km to 120 km, which we defined as the rapid falling stage (RFS), Kepler two-body motion dynamics with acceleration perturbations were applied. All the relevant orbital elements were analytically solved and numerically simulated by the Runge–Kutta integration method. Thus, the decay orbit for large spacecraft from 350 km to 120 km altitudes can be evaluated by the platform we built. All the predicted results were broadly consistent with the measurement data. The findings in this paper can be further applied into the orbit determination of noncooperative spacecraft.

#### 1. Introduction

The reentry of large-scale spacecraft into the earth atmosphere at the end of their lifetime can be divided into uncontrolled and controlled reentry. In uncontrolled reentry, which is mostly caused by communication problems, the spacecraft naturally deorbits and falls into the dense atmosphere [1]. In space history, there are many similar cases in which failed or end-of-life large-scale spacecraft reentered back to the earth, including the Skylab space station of the United States, which ablated on 12 July 1979 during the reentry process, although some pieces fell to the earth because of the large volume and mass of Skylab [2, 3]. The first Chinese space laboratory Tiangong-1 spacecraft, which was launched on 21 : 16 : 03 29 September 2011 (UTC +8) from the Jiuquan Satellite Launch Centre (JSLC) [4]. Its data services were terminated on 16 March 2016. Subsequently, the orbital altitude of spacecraft gradually decayed, and it finally reentered the atmosphere on 2 April 2018. The reentry area was located at the centre zone of the Pacific, and most of its components were ablated and destroyed during the reentry process. Other cases include the Salyut-7 space station of the former USSR and the Envisat environmental satellite of ESA [5, 6]. The flight attitude, trajectory, and dropping zone of these spacecraft are all out of control, dictated by the perturbations of aerodynamics and earth gravity, and only accurate forecasts can significantly reduce the threat of disintegration parts falling on densely populated areas [7]. Although hypersonic aerodynamic heating destroys most parts of large-scale spacecraft [8], and the probability of the practical damage of debris on the ground is extremely low, the orbit decay prediction accuracy improvement of large-scale spacecraft in very low earth orbit (VLEO) is still theoretically and practically significant and beneficial for orbit determination as well as the collision warning of noncooperative targets in VLEO [9, 10].

Generally, the natural reentry process of a spacecraft is related with three key stages: the orbit decay process before reentry, disintegration process, and the spread of wreckage pieces after reentry. The deorbit point of reentry for spacecraft can be treated as the separation between the first two stages. A more reliable position and velocity prediction of this point close to the true deorbit point of the spacecraft can be captured, and a higher accuracy reentry prediction can be realised. Thus, it is necessary to develop a proper orbit propagation model with high precision for the prediction of the position and velocity state vector of spacecraft and provide deorbit information (reentry date and time, velocity, and position) for the inputs of the calculation of disintegration and wreckage spread range. The running orbit of spacecraft before decay was approximately 350 km high; as the atmospheric density is very rarefied, and the perturbation magnitude of aerodynamic forces is extremely low, the decay rate of spacecraft in this altitude range is extremely slow, usually lasting for months or years. At this slow descending stage (SDS), only the secular perturbation factors need to be considered. A complex numerical propagation model has been proposed to predict the orbital elements of a target spacecraft. However, reliable predictions are required from the osculating orbital elements and the step-wise integration of the trajectory [11]. This disadvantage can be overcome using analytical orbit models. The general method for this long-term orbit determination is the use of the mean orbit elements model, in which the two-line element (TLEs) combined with the simplified general perturbation 4 (SGP4) model is a very efficient and feasible way to solve such problems. SGP4, the first North American Aerospace Defense Command (NORAD) orbit model [12], was developed in 1966 based on the analytical theory of Kozai [13]. This model is an analytical orbit model for low earth orbit (LEO) satellites. The SGP4 inputs consist of the TLEs, although the generating method has not been made public by NORAD. The main errors of the TLEs + SGP4 orbital prediction method come from two parts: one is the error of TLEs, and the other is that only second-order oblateness perturbation is considered, and the exponent atmospheric density is used [14]. Owing to the unknown drag term in TLE during propagation, it is easy to cause cumulative errors in propagation, especially for large-scale spacecraft, which was tumbling, and where the mass–area ratio changed significantly. Thus, the traditional fitting drag coefficient method using historical data hardly satisfies the prediction accuracy [15–17].

As for the rapid falling stage (RFS), the decay of large-scale spacecraft in the outer atmosphere experienced a transition from near–free-molecule to near-continuum flow, and the density increases exponentially as the altitude decreases. This process lasts only several weeks or even several days. When there is a lack of communication information from the satellite, only the orbit dynamics model can be relied on to make the propagation as accurate as possible. Currently, nearly all high-precision ephemeris forecasts use numerical prediction, which can take full account of various perturbation factors, especially for the aerodynamic forces and moments, and are applicable for short-term orbit determination [18, 19].

Overall, the aerodynamic perturbation effects are extremely significant for uncontrolled spacecraft in VLEO. Considering a drag error of approximately 5%, the maximum position error of the VLEO target for 24 h propagation can reach 64 km. Preliminary research by Li and Zhang [20] showed that for orbit decay simulation of VLEO uncontrolled targets, the aerodynamic effects covering all flow regimes on orbit and attitude determinations must be considered. Thus, the orbit dynamic model combined with the aerodynamic force and moment coefficients should be a good way to solve such problems for large-scale spacecraft. Hart et al. [21] performed fast orbit propagation and sensitivity analysis of uncertain parameters of low-orbit targets using the analytical free-molecule aerodynamic model, as well as the orbital lifetime prediction. Sommer et al. [22] analysed the attitude motion and cross-section of the large-scale spacecraft reentry. Nazarenko utilised TLE data to predict the reentry time of large-scale spacecraft, considering the disturbances of a random atmospheric environment.

We focus on the perturbation effects of aerodynamic components on the orbit decay and reentry prediction of the uncontrolled large-scale spacecraft. A novel local rapid algorithm of aerodynamics covering all fluid regimes from free-molecule to near-continuum flow using the numerical simulation method for large-scale spacecraft is built and combined with the orbit dynamic models creatively. For the SDS orbital modelling, the drag term of the TLEs for SGP4 is calculated by performing the aerodynamic model at each time step. Then, for the RFS trajectory simulation, the aerodynamic perturbation and the attitude change influence increase significantly, and a six degree-of-freedom (DOF) orbit-attitude dynamics model of large-scale spacecraft is built and directly integrated to realise high-accuracy prediction. The three force and moment coefficient components acting on the target are also calculated and updated at each time step by applying the aerodynamic model. The entire decay altitude orbital prediction of an uncontrolled spacecraft can be predicted using the proposed model. The results show that our methodology can greatly improve the accuracy, which validates the feasibility of combining unified aerodynamics computation with the orbital decay and reentry prediction of near orbital targets.

#### 2. Methods

This section presents the modelling method, including the mean orbital elements model of SGP4, osculating orbital elemental methods, and an aerodynamics model covering all flow regimes for large-scale spacecraft. First, some hypotheses should be stated to simplify the present complex problem. As the uncontrolled spacecraft is not in a three-axis stable flight attitude, random tumbling occurs in the perturbation action, and the attitude dynamic mechanism of uncontrolled large-scale spacecraft is too complex to be fully solved. For the SDS modelling, we neglected the attitude change of the spacecraft and assumed that the drag area of large-scale spacecraft is constant in the long-term prediction. For the RFS modelling, the attitude initial information can be obtained from the filter parameters of the flight data, and 6-DOF dynamic equations are solved.

The sketch of the reference frames is shown in Figure 1; the orbital coordinate frame and body frame , as well as the J2000 earth-centred inertial (ECI) frame are used to illustrate the perturbation and attitude evolution of large-scale spacecraft before the reentry.

##### 2.1. Mean Orbital Elements’ Model

The SGP4 model is an orbital analytical prediction model released by NORAD for LEO targets and should be propagated combined with the TLE data for long-term prediction. The main perturbation factors considered in SGP4 are aerodynamic drag and the J2 term oblateness of the earth [14].

The osculating elements transformed from the position and velocity vector of large-scale spacecraft in orbit at the epoch time are given in the J2000 earth-centred inertial (ECI) frame. Note that the computation inputs and outputs of SGP4 are required to be expressed in the true equator mean equinox (TEME) frame. Therefore, transformation between J2000 and TEME was first required [11]. Then, the state variables can be described in the form of NORAD-No. 1 TLE as where is the state vector at time , represents the SGP4 function, and represents the mean orbital elements in TLE at time , including the drag term , which can be written as

The are the six mean Keplerian orbital elements, which is eccentricity, inclination, right ascension of the ascending node, argument of perigee, mean anomaly, and mean motion separately, and the drag term is related to the ballistic coefficient by [23].

The ballistic coefficient is given by , where is the dimensionless aerodynamic drag coefficient, is the area-to-mass ratio of spacecraft, and is the atmospheric reference density, , where is the radius of the earth [24].

##### 2.2. Osculating Orbital Elements’ Model

For the short-term prediction of the RFS before reentry, the orbital dynamic equations can be stated by position and velocity vectors in the J2000 ECI frame. Quaternions are used for the attitude dynamic model to avoid singularity. Thus, the 6-DOF dynamic equation of the spacecraft can be expressed as [25]. where is the position vector, is the velocity vector, is the quaternion from the orbital frame to the body frame, is the attitude angular velocity, is the earth’s gravitational constant, is the inertia matrix of the spacecraft, is the -body gravity perturbed acceleration, and is the oblateness perturbation acceleration of the earth. The joint gravity model 3 (JGM-3) gravitational potential model is applied in this study, is the aerodynamic lift perturbation in the J2000 coordinate system, is the solar radiation perturbation, and is the aerodynamic drag perturbation in the J2000 coordinate system, which can be transformed by to body frame as where is the direction cosine matrix of quaternions, is the transformation matrix from the orbital coordinate system to the J2000 coordinate system, and can be expressed as three components of the aerodynamic coefficients in the body coordinate system. where , , and are the axial force, normal force, and side force coefficients, respectively. is the drag area of the target spacecraft, is the atmospheric density (here, the density model NRLMSISE00 is applied in osculating orbital element models), and is the relative velocity between the spacecraft and atmosphere, which can be expressed as where is the rotational angular velocity of the earth.

As for the disturbed moments of Equation (4), is the gravity gradient moment, is the solar pressure moment, is the geomagnetic moment, and is the aerodynamic moment, which can also be listed as three components in the body coordinate system by applying the rolling moment coefficient , pitching moment coefficient , yawing moment coefficient , drag area , and reference length :

##### 2.3. Aerodynamic Characteristics Covering All Flow Regimes

The aerodynamic characteristics of spacecraft in VLEO vary significantly along the orbital altitude and attitude. The flight of large-scale spacecraft in VLEO is a hypersonic flow problem from hundreds of Knudsen to dozens of Knudsen in the rarefied free-molecular flow regime for complex configurations under multiphysics fields [26]; the refine of the aerodynamics characteristics of spacecraft in rarefied transitional regime can be achieved by the numerical method like direct simulation Monte Carlo (DSMC) and the gas-kinetic unified algorithm solving the Boltzmann model equation (GKUA). However, the time cost of numerical solving is so high that cannot be accepted during the reentry prediction, and the calculation efficiency of aerodynamics covering all flow regimes according to the fast response of spacecraft reentry is also vital in the orbit prediction process. Here, a unified fast algorithm computational method for aerodynamic force driven by fine numerical simulation covering various flow regimes developed by Li et al. is applied in this paper. The specific construction method can be found in References [20, 26–29]; here, we detailly describe the construction and application of this aerodynamics method in our paper.

###### 2.3.1. Hypersonic Aerodynamics Local Algorithm

The continuum flows as well as the free-molecular flow regimens are firstly considered, for high rarefied flow, most of which are free-molecular fluid states; the modified Nocilla wall reflection model corresponding to different materials can be used to calculate the pressure and friction coefficients [30]. Under the approximation of the Maxwell equilibrium state gas molecular velocity distribution function, the pressure and friction coefficients on each surface element can be calculated and validated. For the continuum flow regime, the Knudsen number tends to be exceedingly small, and the modified Newton inviscid flow theory can be used. The pressure coefficients can be calculated by the modified Newton inviscid flow theory in a continuum flow regime as [21] where , , and are the pressure, density, and velocity of free inlet flow, respectively; is the surface angle; and is the pressure coefficient under the maximum surface pressure , which is also the stagnation point after a normal shock wave. can be calculated as where is the inlet flow Mach number, and is the specific heat ratio. In the approximation of the Maxwell equilibrium gas distribution, the pressure and friction coefficients for each surface element are where and are the surface and inlet temperatures, respectively; and are the normal and tangential momentum adaptation coefficients, respectively; is the error function; and is the speed ratio, which is defined as where is the universal gas constant. The semi-empirical bridge function can be applied in the transition region to smoothly connect the coefficients of continuum and free-molecular flow regimes, which can both tend to be stream functions. For a given surface angle , the pressure and friction coefficients on the local surface element are [31]. where and are the pressure and friction bridge functions, respectively, which independently rely on the parameters , , and .

As for the transition zones connecting the adjacent regions of rarefied and continuum flow regimes, the modified Boettcher/Legge theory was developed to describe the nonsymmetric bridge correction function of pressure and friction coefficients, of which the nonsymmetric pressure coefficients can be expressed as follow: where , , and are the adjustable related parameters corresponding with the configuration and surrounding flow characteristics of spacecraft.

As for the nonsymmetric friction coefficients, it can be expressed as follow: where , , and are the adjustable related parameters.

###### 2.3.2. Fast Algorithm Covering All Flow Regimes for Target Spacecraft

The six local correlation parameters , , , , , and should be adjusted and determined by the refined aerodynamics coefficients data from wind tunnel test or numerical method. According to the flow characteristics of rarefied transitional regime, refined calculation of aerodynamics by DSMC or GKUA can be launched on high-performance computers [27]. Considering the complex configuration and surrounding flow characteristics of spacecraft in this paper, the GKUA solving the Boltzmann model equation in References [20, 26, 28] was applied in this paper to achieve the refinement aerodynamics coefficients in the flow of 350 km~120 km altitudes, by which the six local correlation parameters can be adjusted and determined to best fit the numerical results.

Thus, a unified fast calculation technique for the aerodynamics of a VLEO spacecraft covering all flow regimes is developed. For the large-scale spacecraft in this paper, a triangular unstructured grid can be used to represent their irregular configurations, and when the surface element is sufficiently small, the errors of the integrated aerodynamics coefficients on the normal and tangential forces of all surface elements along the entire geometrical shape of the target spacecraft would also be small. Finally, the local fast algorithm of aerodynamics for the large-scale spacecraft specified configuration can be modified, and the aerodynamic coefficients of the target spacecraft on the corresponding orbital altitude and attitude can be obtained.

Aiming at the aerodynamic modelling for uncontrolled flight of large-scale spacecraft, which experiences various attitudes and multiphysics as well as multiscale nonequilibrium reentry state during the long orbital decay process, the three axial aerodynamic force coefficients, , and as well as the three aerodynamic moment coefficients , , and in body frame can be directly integrated and calculated by the pressure and friction coefficients on the local surface element of large-scale spacecraft. The aerodynamic frame of the large-scale spacecraft in this paper is shown in Figure 2, from which the transformation relationship between the aerodynamics in body frame and velocity frame can be clearly seen. Thus, for the drag coefficients and lift coefficients used in SDS model, it can be transformed by the axial aerodynamic force coefficients and the attitude angles including the angle of attack and angle of side slip , which are

To further validate the drag coefficient data over a range of altitudes is for the large-scale spacecraft in this paper, the numerical methods of GKUA as well as the unified local fast algorithm are applied to calculate the drag coefficient data at some typical state of height from 350 km~120 km with ; the comparison results are shown in Table 1. It can be found from the data that the local fast algorithms are very close to the numerical methods from 350 km~120 km, which could be a validation for our aerodynamics rapid algorithm.

##### 2.4. Integration Platform

Based on the orbital perturbation and rarefied aerodynamic models build earlier, the simulation platform for the orbit decay prediction of large-scale spacecraft before the reentry can be integrated, as shown in Figure 3. Firstly, the initial ephemeris and attitude of angle (AOA) of spacecraft are input to calculate the orbital elements for the integration, then the AOA including and ; the velocity and orbital height at time are transferred into the rarefied aerodynamic model to calculate the aerodynamics coefficients that the orbital perturbation models used. Then, the orbital perturbation model combined with the rarefied aerodynamics model begins to iterate and integrate. The represents the average orbital height of the spacecraft, the SDS model is used when the , or the RFS model is used below the 250 km orbital height. And the integration is finished when the decayed to 120 km, which we marked as the reentry height of the large-scale spacecraft.

#### 3. Input Data Filtering

To further validate the accuracy of the integrated orbit decay prediction model we developed, the external measurements TLE data are applied for the comparison of our integrated model. Before the discrete time series TLE data was used, the orbit perturbation model for LEO and UKF filtering method is used to correct the measurement errors of the original TLE data [32]. Table 2 illustrates the initial values of TLEs with drag term for simulation before and after filtering.

#### 4. Results and Discussion

Based on the integration platform we developed, both the SDS and RFS models were utilised to propose the orbit decay prediction of the large-scale spacecraft from approximately 350 km to 120 km height to verify the dynamic models and methodology we employed.

##### 4.1. Aerodynamic Results

Figure 4(a) illustrates the variation in the aerodynamic coefficients as well as the moment coefficients along with the variation in the attack angle. Figure 4(b) plots the drag coefficients along with the altitude from 350 to 120 km. Figure 4(c) plots the lift coefficients along with the altitude from 350 to 120 km. It can be seen from Figure 4(b) that the drag coefficients rapidly increased from 350 to 300 km, and then gradually increased as the orbital altitude decreased below 150 km. As for the lift coefficients, which is very small and may not play a significant role in this work to influence the altitude/attitude dynamics of the spacecraft during its orbital decay, mostly for the reason that the rotating velocity of the spacecraft is too high to make its shape like a sphere, and the orbit drifts driven by the lift are very small and mostly offset during the period of near circular orbit in VLEO.

**(a)**

**(b)**

**(c)**

##### 4.2. Orbit Decay Prediction

For the orbit decay process above 250 km orbital altitude, the SDS model was solved to predict the secular perturbation effects on the orbit of large-scale spacecraft, especially for the aerodynamic drag and lift perturbations. To reduce the error of the initial values of the formatted TLE data set for our modified SGP4 model, the UKF method was utilised to minimise the errors of mean orbital elements after iteration combined with the SDS model, and a series of orbit decay prediction data propagated by the TLE + SGP4 model filtered by TLE data sets and UKF are applied for the compassion from height of 350-120 km. The initial ephemeris and AOA data of large-scale spacecraft are carefully chosen for the orbit decay propagation.

As the altitude of large-scale spacecraft decreased, the influence of aerodynamic forces increased dramatically. We first predicted the orbit of large-scale spacecraft above 250 km in the SDS, for which the altitude range of the effects of aerodynamics is relatively low. Figures 5 and 6 illustrate the comparison data of orbital information between the propagation and observation data from TLE sets of NORAD, including the orbital average altitude, perigee altitude, and apogee altitude, as well as the orbital elements including semi-major axis, eccentricity, and inclination. The results show that the simulation data are in exceptionally good agreement with the observation data, which proves the effectiveness of the SDS propagation model we built.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

The influence of aerodynamic forces plays a more significant role below 250 km. The prediction results of the RFS model from approximately 250 km to 120 km, including the orbital altitude and element information, are illustrated in Figures 7 and 8. It can be seen from the results that the propagation data trends are also in good agreement with the observation data transferred from the TLE sets. This supports a good validation of the reentry date and location predictions of large-scale spacecraft using this method.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

Based on the aerodynamic results, the 6-DOF orbit-attitude coupling model can be solved to analyse the attitude motion characteristics of large-scale spacecraft, although the precise attitude estimation for uncontrolled spacecraft is extremely difficult. We found a significant rolling phenomenon of uncontrolled large-scale spacecraft’s attitude, which was also validated by Lin et al. [33]. The analysis dates ranged from 12 August 2017 to 4 April 2018, and the statistical interval in prophase was one week and then decreased to one day in anaphase close to reentry time. Figure 9 shows a comparison of the prediction results and observation data, which shows that there is a significant change in large-scale spacecraft’s attitude from prophase to anaphase, and the overall trend of prediction coincides with the observation data. The rolling velocity of large-scale spacecraft slowly increased in prophase, from about 33°/min to 102°/min, and then secularly changes; the trend of rolling velocity demonstrates a relatively stationary periodic rolling. Close to the reentry time, the maximum rolling velocity surge increased sharply, approximating a linear increase from 121°/min to approximately 215°. This type of accelerated rolling phenomenon increases the difficulty of reentry point prediction.

Figures 10 and 11 plot the attitude angle as well as the angular velocity variations of large-scale spacecraft in this paper by the 6-DOF orbit-attitude coupling model. As for the attitude numerical integration algorithm, the time-step needs to be small enough to avoid singularity phenomenon during the attitude iteration. Here, we only display the results of attitude angles during about 4.5 days before the reentry. It can be seen from Figures 10 and 11 that the variations of both attitude angles as well as the attitude angular velocity have a certain degree of randomness, because the limitations of attitude measuring and observations techniques for the noncooperative large-scale spacecraft are still not precise enough currently worldwide, only the preliminary investigations of attitude information are presented in this paper. And both the stochastic theory as well as the statistical methods still needs to be applied into the attitude angles investigation in the future. All in all, the attitude integration results are converged and reasonable.

##### 4.3. Reentry Prediction and Error Comparison

In this study, we further analysed the reentry time prediction error using our methodology comparing with the other representative algorithm like conventional SGP4 perturbation method. Figure 12 shows the evolution of the relative error in the estimation of the residual lifetime for the large-scale spacecraft’s reentry prediction, where we compare two methods of orbital prediction: one is to make the reentry prediction solely using the rapid algorithm of perturbation model built by us in this paper from 350 km to 120 km, and the other is to predict the reentry using the representative SGP4 model from 350 km to 120 km. The comparison results show that the error when using only the SGP4 model is much larger than that when applying our model below 250 km altitude. In addition to the inherent error in the short-term prediction of mean orbital perturbation models, the aerodynamics play a more significant role as the altitude decreases, and our rapid model combined with the three aerodynamic forces as well as the three aerodynamic torques are more suitable for orbital prediction below 250 km owing to the closer integration with aerodynamics.

##### 4.4. Simulation Rapidness of Integration Platform

The rapidness of our integration platform based on the rapid algorithm of rarefied aerodynamics is crucial for the reentry prediction, especially for the last few days before the reentry of uncontrolled large-scale spacecraft. Because the orbit decay speediness of the large-scale spacecraft is very high and the prediction time must be faster than the real altitude decay time of the spacecraft before the reentry to the atmosphere as much as possible, which could earn more time for the risk assessment of the fall and disintegration of large-scale spacecraft after the reentry. Besides on, the computation efficiency is also a representative factor of the coupling integration process like the accuracy. The most time cost part of the platform is the simulation of 3D surrounding flow of complex configuration of spacecraft, comparing with the numerical method of GKUA; the calculation of the unified local fast algorithm could improve as much as 70%. As for the overall integration platform, the calculation for 7 days’ orbit decay of real physical time would cost 1 day integration time, which could be fast enough for the reentry prediction of large-scale spacecraft.

#### 5. Conclusion

In this paper, a unified modelling method, including orbit and aerodynamics, was proposed to propagate the large-scale spacecraft’s decay orbit from 340 km to 120 km. The UKF method was also used to filter the discrete measurement TLEs for the simulation results comparison and verification of our proposed integrated orbit perturbation model combined with the rarefied aerodynamic algorithm, and the results showed good agreement. This supports the fact that the method and model proposed in this paper can be applied in the orbital prediction for large-scale uncontrolled spacecraft after the service is terminated, which is the foundation of point prediction of remaining debris after the disintegration of ablation as well as the falling risk assessment. The methodology proposed in this paper can be also further developed for orbit determination and collision alerts of spacecraft and debris in VLEO, even for noncooperative space targets. The study findings in this paper also could highly promote the combination of elaborate rarefied aerodynamics models and results with the orbit in aerospace field, which could be a strong support for the high accuracy of large-scale spacecraft flight control, as well as the aerodynamic design of novel reusable spacecraft manoeuvring on orbit.

#### Data Availability

The simulation data used to support the findings of this study were supplied by Zhihui Li under license and so cannot be made freely available. Requests for access to these data should be made to Zhihui Li ([email protected]).

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The authors are grateful for the financial support from the China Manned Space Engineering Office, the National Natural Science Foundation of China (11902331, 91530319, and 52106215, 11325212), and the National Key Basic Research and Development Program (2014CB744100). They would also like to thank the Beijing Aerospace Control Center for providing support.