Prediction of the space debris spatial distribution on the basis of the evolution equations

Acta Astronautica 100 (2014) 47–56

Contents lists available at ScienceDirect

Acta Astronautica journal homepage: www.elsevier.com/locate/actaastro

Prediction of the space debris spatial distribution on the basis of the evolution equations A.I. Nazarenko n Scientific Technological Center KOSMONIT, Roscosmos, Profsoyuznaya Street 84/32, 117997 Moscow, Russia

a r t i c l e i n f o

abstract

Article history: Received 24 January 2014 Accepted 25 February 2014 Available online 19 March 2014

The numerical–analytical technique for long-term prediction of the space debris (SD) spatial distribution, based on derivation and solution of new evolution equations, is developed. These equations are represented in two forms — difference and differential. In the latter case the problem is reduced to integration of the system of two ordinary differential equations. A high efficiency of the proposed technique, as compared to the traditional approach, is demonstrated. & 2014 IAA. Published by Elsevier Ltd. All rights reserved.

Keywords: Space debris Orbital elements Statistical distribution Evolution equations

1. Introduction It is assumed that two factors have effect on the evolution of the space debris time–spatial distribution during the prediction of space debris contamination in the region of low orbits (LEO) [1–8]: the increment of a quantity of new objects as a result of launches, technological operations, explosions, accidents, etc., and the atmospheric drag that results in lowering the perigee altitude of space objects (SOs) and their reentry. Two approaches are used for solving the considered problem. The traditional (deterministic) approach is widely applied by specialists. For example, NASA's contemporary “EVOLVE” model simulates the after-effect of all known satellites' launching and destruction events, as well as similar events possible in the future. For each object (or a group of objects) the vector of initial conditions is formed. Prediction is performed using the traditional motion models. To estimate the danger of collision of a pair of satellites D. Kessler's methods are used. The results have been summarized for a great number of objects. The distinction in modeling over the future time interval lies

n

Tel.: þ 7 8 499 793 59 00. E-mail address: [email protected]

http://dx.doi.org/10.1016/j.actaastro.2014.02.023 0094-5765/& 2014 IAA. Published by Elsevier Ltd. All rights reserved.

in the fact that all fragmentation events are formed according to the Monte Carlo technique. In this case several prediction operations are performed. Obviously, this approach is very labor-consuming; it can be implemented on rather powerful computers only. However, even on the contemporary large computers the traditional approach does not correctly predict the distribution of small-sized space debris fragments. The labor consumption in solving the problem in question can be judged by quotation from paper [8] given below: “As one of the more time consuming operations of our model deals with the orbital propagation of the sixth orbital elements for each objects of the population, the code of MEDEE has been designed to take advantage of massively parallel, computer system available at CNES. This means that the orbital propagation module has been parallelized, in order to propagate the population at each time-step over all available cores. The computer system in which MEDEE is executed is formed by 360 cores summing a total RAM of 24 Go and an overall computing power of 4 Tflops/second.” In spite of its labor-consuming character, this approach does not guarantee the model's adequacy. The accuracy of destruction after-effect modeling is unknown. The model

48

A.I. Nazarenko / Acta Astronautica 100 (2014) 47–56

tuning according to available measurement data turns out to be a rather difficult task. Another (statistical) approach to formation of NearEarth Space (NES) contamination sources and to the situation prediction is applied to the Russian Space Debris Prediction and Analysis (SDPA) model. Its main distinction lies in the fact that, instead of the data on particular launches and on SOs destruction events, the model uses the averaged data on SOs spatial distribution and on the number of yearly formed objects of various sizes. This particular approach is substantiated by the following reasons:  The SD number varies insignificantly (by few percents only) during a year. Therefore, the more detailed (in time) modeling of SD sources is excessive; it highly complicates the model, virtually not influencing the accuracy. Note. This statement does not exclude the possibility and expediency of detailed modeling of breakup consequences for short time intervals, when the SD “cloud” remains rather compact. However, the cloud scattering process is known to proceed rather quickly, as a rule. The duration of this process is about 1 month.  The reasons and conditions of satellites' fragmentations, which resulted in forming the majority of smallsized SD fragments, are extremely diverse. Therefore, one can hardly expect that results of modeling the consequences of all known fragmentations (the number of particles, the fly-away velocity) are accurate enough. The level of errors of such modeling is unknown. Hence, the approach, based on the averaged data, does not seem to be worse, but looks even more preferable.  The previous statement is even more valid for the future time instants, where reasons and circumstances of fragmentations are unknown. The use of the averaged data on the intensity of new SD formation seems to be optimal for performing the predictions.  One more reason in favor of applying the averaged data on new SD formation intensity is based on the fact that the averaged data contain a smaller number of parameters in comparison with the detailed data. It is known that the minimum number of any prediction model's parameters to be updated on the experimental data basis is desirable under conditions of limited measurement information. Therefore, the model having a smaller number of parameters to be tuned will provide better prediction accuracy as compared to models having a larger number of such parameters.

In the Russian SDPA model, during the situation prediction with allowance for atmospheric drag of SOs, various objects with perigee altitude o2000 km are considered. We shall suppose that among all variable SO parameters only the perigee altitude essentially influences the evolution of altitude distribution of SO number. The other orbital elements will be designated by Э. We subdivide the whole set of objects with different elements Э into some finite number of sub-sets (groups) with

elements Эl, l ¼1,2,…,lmax . Let p(t,h) be the density of the perigee altitude distribution for objects from the selected group at time t. Then we state the problem of studying the laws of density variation in time. Index l will be omitted hereafter in analyzing the distribution evolution for some particular SO group. In calculating the evolution of altitude distribution of SO number the following factors are taken into account: the atmospheric drag at altitudes of up to 2000 km; the sub-division of all SOs in parameters into the groups which differ in size d, eccentricity e and ballistic factor kb; the initial altitude distribution of SOs of various types; the expected intensity of formation of new SOs of various types as a result of launches and explosions: p(h,t)new is the increment of SO number at various altitudes per time unit; the non-stationary character of factors is taken into account; namely, the atmospheric density in connection with solar activity variation during the 11-year cycle. The perigee altitudes are of special importance among the listed orbital parameters taken into account. While rendering the basic influence on orbit lowering, it changes (the perigee altitude is finally lowered resulting in SO reentry). Therefore, the perigee altitude is one of the arguments of evolution equation and it does not enter the partition. The ballistic factors and SO size have absolutely different character of influence: they do not virtually change during the orbital evolution process. The eccentricity has intermediate character of influence: it changes, generally speaking, under the atmospheric drag effect, but this change does not play an essential part, because large parts of SOs have orbits with low eccentricities.

2. Derivation and solution of evolution equations Let us consider the technique for prediction of the distribution of SOs p(h,t) in the perigee altitude. We will derive the relationship for determining this density at various altitudes at the time instant tþΔt. In the discrete partition of argument h over some specified interval with the step Δh the initial distribution p(h,t) is specified on this grid of the values of arguments. Fig. 1 depicts the values of distribution p(h,t) for two values of argument h and hþΔh. The number of objects in

Fig. 1. Scheme of change of SO distribution in the perigee altitude. Scheme of change of SO distribution in the perigee altitude.

A.I. Nazarenko / Acta Astronautica 100 (2014) 47–56

this interval of perigee altitudes is equal to NðtÞh;h þ Δh ¼ pðh; tÞΔh

ð1Þ

The rate of decrease of the perigee altitude of objects with altitude h is denoted as V(h,t). Similarly, the rate of decrease of the perigee altitude of objects with altitude hþΔh is denoted as V(h þΔh,t). Let us consider what will happen with the distribution p(h,t) in some time, namely, at time (tþΔt). It is obvious that the quantity of objects (1) in the considered range of altitudes will change as a result of three circumstances, namely: 1. The part of objects in the vicinity of altitude h will decrease so that their perigee altitude will be less than the altitude hþΔh. The number of these objects is equal to Nðt; t þ ΔtÞð1Þ h þ Δh

¼ Vðh; tÞΔt  pðh; tÞ

ð2Þ

2. The parts of the objects in the vicinity of altitude 4hþΔh will decrease to the extent that their perigee altitude will be less than the altitude hþΔh. The number of these objects is equal to ¼ Vðh þΔh; tÞΔt  pðhþ Δh; tÞ Nðt; t þ ΔtÞð2Þ h þ Δh

¼ pðh; tÞnew Δh  Δt Nðt; t þ ΔtÞð3Þ h þ Δh

ð4Þ

NðtÞh;h þ Δh Nðt; t þ ΔtÞhð1Þþ Δh

þ Nðt; t þ ΔtÞð2Þ þ Nðt; t þ ΔtÞhð3Þþ Δh h þ Δh

ð5Þ

In its content, Eq. (5) is just the evolution equation in discrete form for prediction of SOs distribution in perigee altitude per step in time. The required density of distribution of perigee altitudes at time instant tþΔt is determined based on the estimate of Eq.(5) by a simple formula: pðh; t þ ΔtÞ ¼

Nðt þ ΔtÞh;h þ Δh : Δh

¼ Vðh; tÞ U dt U pðh; tÞ; Nðt; t þ dtÞð1Þ h þ Δh Nðt; t þ dtÞð2Þ ¼ Vðhþ dh; tÞdt  pðhþ dh; tÞ h þ dh
 ∂pðh; tÞ dh ¼ Vðh þdh; tÞdt pðh; tÞ þ ∂h Nðt; t þ ΔtÞhð3Þþ dh ¼ pðh; tÞnew dh  dt; Nðt þ dtÞh;h þ dh ¼ pðh; t þ dtÞ Udh ¼ pðh; tÞ U dh þ

The total number of objects in the considered range of perigee altitudes at time tþΔt will be equal to Nðt þ ΔtÞh;h þ Δh ¼

the situation prediction possess a rather high uncertainty, which makes senseless excessive details. The general tendency of the p(h,t) distribution evolution in the region of low orbits is to lower the perigee altitude of satellites under an effect of atmospheric disturbances. This fact characterizes the process of space selfpurification from space debris as a result of objects drag in the atmosphere. The stronger this process proceeds, the lower the perigee altitude of satellites (the higher the perigee altitude lowering rate). Of particular interest is the transformation of relations (5) and (6) into the form of the differential equation, assuming the discrete values of Δh and Δt to be infinitely small quantities (dh and dt). In this case expressions (1)–(4) and (6) take the form NðtÞh;h þ dh ¼ pðh; tÞdh,

ð3Þ

3. As a result of launches and explosions, the new SOs will be added. The number of these objects is equal to

ð6Þ

Thus, relations (5) and (6) allow calculating the changes in the perigee altitude distribution of SOs number in the situation prediction for one step in time. Successive application of these relations in a cycle in altitude and in time ensures the solution of the problem under consideration. For solving Eq. (5) a specially developed numerical– analytical technique is applied. In choosing the number of SO set partitions into groups a compromise is required between the detailed character of the analysis and feasibility of the algorithm: an extremely detailed partition requires a lot of memory and increases the computation time. Also, one should keep in mind that the initial data for

49

∂pðh; tÞ U dt U dh: ∂t

Here Vðhþ dh; tÞ ¼ Vðh; tÞ þ ð∂Vðh; tÞ=∂hÞ U dh: Substituting these expressions into (5), after eliminating identical terms in the right- and left-hand parts, as well as the multiplier dt  dh, we obtain the following equation in partial derivatives, which describes the evolution of the distribution of SO number in the perigee altitude: ∂pðh; tÞ ∂pðh; tÞ ∂Vðh; tÞ ¼ Vðh; tÞ U þ pðh; tÞ U þ pðh; tÞnew ∂t ∂h ∂h

ð7Þ

Consider now the full differential of the p(h,t) distribution at the point (h,t): dpðh:tÞ ¼

∂pðh; tÞ ∂pðh; tÞ U dh þ Udt: ∂h ∂t

ð8Þ

With regard to determination of the perigee altitude lowering rate dh=dt ¼ Vðh; tÞ

ð9Þ

substitution of expression (7) into (8) and the use of designation ∂Vðh; tÞ=∂h ¼ Aðh; tÞ

ð10Þ

results in the equation dpðh; tÞ ∂Vðh; tÞ ¼ U pðh; tÞ þpðh; tÞnew ¼ Aðh; tÞ U pðh; tÞ þ pðh; tÞnew dt ∂h

ð11Þ

Here, in accordance with Eq. (9), altitude h is the wellknown function of time hðt Þ. In its content, expression (11) represents the evolution equation in the differential form. Namely, this is the linear inhomogeneous differential equation of the first order. As known, the general solution of Eq. (11) is the sum of the general solution of the

50

A.I. Nazarenko / Acta Astronautica 100 (2014) 47–56

ð12Þ

will occur to be in a wider range of altitudes. On the basis of difference equation, one can easily compose the differential equation for the range Δh(t). Applying the same approach in deriving Eq. (7), we obtain

ð13Þ

dΔhðtÞ dVðh; tÞ ¼ U ΔhðtÞ ¼  Aðh; tÞ UΔhðtÞ: dt dh

homogeneous equation dx=dt ¼ AðtÞ U x which has the form
Z t  xðtÞ ¼ x0 U exp AðξÞ Udξ ¼ x0 U uðt; t 0 Þ; ξ ¼ t0

and the partial solution of the inhomogeneous equation. By using the function uðt; t 0 Þ, the solution of Eq. (11) assumes the form Z t pðh:tÞ ¼ uðt; t 0 Þ Upðh0 ; t 0 Þ þ uðt; ξÞ Upðh; ξÞnew Udξ: ð14Þ ξ ¼ t0

An important result of the performed analysis is the fact that the solution of a considered problem was reduced to solution (14) of a rather simple ordinary differential Eq. (11). Let us consider in more detail the derivative ∂Vðh; tÞ=∂h, which in (11) is designated asAðh; t Þ. We will take into account the circumstance, that for the given orbital elements of SO and its known ballistic coefficient the rate of lowering the perigee altitude Vðh; tÞ is proportional to the atmospheric density at perigee ρðh; tÞ. As the altitude increases, the density of atmosphere decreases. For the exponential altitude dependence of atmospheric density the derivative dρðh; tÞ=dhhas the form dρðh; tÞ ρðh; tÞ ¼ dh Hðh; tÞ

ð15Þ

where H(h,t) is the so-called scale height (the homogeneous atmosphere altitude). Taking into account these data the derivative can be written as Aðh; tÞ ¼

dVðh; tÞ Vðh; tÞ ¼ dh Hðh; tÞ

ð16Þ

In the particular case when coefficient A does not depend on time, the solution of (14) is simplified: pðh:tÞ ¼ exp½AU ðt  t 0 Þ Upðh0 ; t 0 Þ Z t exp½AU ðt  ξÞU pðh; ξÞnew U dξ þ

ð17Þ

ξ ¼ t0

Because the value of coefficient A is negative, the value of the function exp½Aðt  t 0 Þ decreases in time. So, two methods of solving the problem in question are substantiated. The first one is based on solution of the evolution equation (5) in a discrete form. The second method consists in applying the analytical solution (14). To conclude this section, we consider the features of the evolution equations. It is obvious from Fig. 1 that all objects, which at time instant t were located in the perigee altitude region (h(t), h(t)þΔh(t)), move with time into the region of lower altitudes. Namely, at time t þΔt they will be in the region (h(tþΔt), h(tþΔt)þΔh(tþΔt)). The width of this region of altitudes Δh(t þΔt) depends on the difference of velocities V(h,t) and V(h þΔh,t)shown as follows: Δhðt þ ΔtÞ ¼ ΔhðtÞ þ ½Vðh; tÞ  Vðhþ ΔhðtÞ; tÞ UΔt:

ð18Þ

It is obvious that, because the value of the difference in square brackets is positive, this quantity increases in time. Thus, the objects, which at the initial time instant were located in some altitude range (h, hþΔh), in some time

ð19Þ

This equation differs from Eq. (12) only in sign. Its general solution has the form ΔhðtÞ ¼ Uðt; t 0 Þ U Δhðt 0 Þ; where


Z Uðt; t 0 Þ ¼ exp 

t ξ ¼ t0

 AðξÞ Udξ :

ð20Þ

ð21Þ

One can easily show that general solutions of homogeneous Eqs. (12) and (19) are related by the equation Uðt; t 0 Þ ¼ uðt; t 0 Þ  1 ¼ uðt 0 ; tÞ

ð22Þ

Thus, we have constructed the dependence for determining the range Δh(t) at any time instant. It was noted above that all objects which at the initial time instant were located in the altitude range (h0, h0 þΔh (t0)) in some time move into other range of altitudes (h(t), h(t)þΔh(t)). We will make sure that solution (14) agrees with this statement. We multiply the right- and the lefthand parts of solution (14) by (20), assuming pðh; tÞnew ¼ 0 and using formula (22). We get pðh; tÞ U ΔhðtÞ ¼ uðt; t 0 Þ U pðh0 ; t 0 Þ U Uðt; t 0 Þ UΔhðt 0 Þ ¼ pðh0 ; t 0 Þ UΔhðt 0 Þ:

ð23Þ

3. Determination of the perigee altitude lowering rate The motion of many SOs occurs in rarefied layers of the upper atmosphere, where the aerodynamic forces are small as compared to their values in the lower atmosphere. However, the durable flight of SOs and the dissipative character of the effect of air drag forces result in essential impact of atmospheric drag on SOs' orbits evolution at altitudes up to 600–1000 km. The action of the drag force on SO of mass m causes acceleration: F ¼ kb ρV 2rel

ð24Þ

where the quantity kb ¼

CxS ðm2 =kgÞ 2m

ð25Þ

is called ballistic coefficient. Here C x is a dimensionless aerodynamic drag coefficient, S is the characteristic area of SO, ρ is the density of the atmosphere, and V rel is the incident gas flow velocity that is equal to velocity of SO flight relative to air. Coefficient C x is a function of many quantities: the geometric shape and orientation of SO, the properties of its surface's material, the composition of the atmosphere and its parameters. In the majority of cases, for the upper layers of the atmosphere, this coefficient varies within the limits 2.0–2.5. Area S is the area of maximum SO's cross section normal to the velocity vector V rel . For oriented spacecraft (SC) it is strictly determined by SC geometry; for non-oriented SC the area S occurs to be

A.I. Nazarenko / Acta Astronautica 100 (2014) 47–56

variable, which results in the necessity of applying its averaged values. The atmospheric density ρ is the spatial–temporal function ρ ¼ f ðh; α; δ; F 10:7 ; ap ; C i ; i ¼ 1; 2::Þ:

ð26Þ

The particular form of this function is specified by the upper atmosphere models [9–11]. Basic arguments of upper atmosphere's dynamic models are the following quantities: h is the altitude of a point above the Earth's surface; α; δ are spherical coordinates of a point in the geocentric inertial coordinate system (GICS); F 10:7 is solar activity index that is equal to the solar radiation intensity emission (1 Solar Flux Unit (SFU)¼ 10  22 W/(m2 Hz) at the wave of 10.7 cm); ap (or K p ) is the index characterizing the geomagnetic activity; t is the time that is used in Eq.(26) in calculating a half-year effect; C i are parameters of the model. The altitude h depends on the radius-vector (r) and latitude (φ) of a point: h ¼ r RE ð1 ε sin 2 φÞ:

ð27Þ

The most essential argument in the model (Eq. (26)) is the altitude. For example, if it changes in the range of 200– 600 km, the density varies 700–3100 times. Within a relatively small range of altitudes the altitude dependence is rather well approximated by the expression   h  h0 ρðhÞ ¼ ρðh0 Þ U exp  : ð28Þ H HereHis the so-called homogeneous atmosphere altitude (the scale height). The effects of α and δ coordinates on the atmospheric density are associated with the so-called diurnal effect. The origin of this term is explained by various degrees of the upper atmosphere warming during day and night. The thermal inertia effect leads to a situation where the maximum of diurnal variation of the density falls on 14–15 h of local time. The diurnal effect amplitude depends on point's altitude and on the solar activity level. The maximum ratio of the day-time density to the nighttime one is achieved at altitudes of 500–600 km and equals 2–3. Consider now, how the perigee altitude changes under the atmospheric drag effect. We will make use of formulas for disturbances of a semi-axis (δa) and eccentricity (δe) per revolution, published in the papers [12–14].

51

of these errors has the order of 10%. The calculated values of atmospheric density have the same order of errors. Thus, in accordance with the definition (9) and expression (29), one applies the following formula for determining the perigee altitude lowering rate: Vðh; tÞ ¼ δh=T

ð30Þ

where T is the satellite revolution period. The basic errors of applying formula (29) for calculating the space debris evolution are associated with a great scattering of possible values of ballistic coefficients. This scattering reaches 4 orders of magnitude or more. 4. Accounting for scattering of possible values of ballistic coefficients Consider the data on scattering of estimates of a drag of cataloged objects. Fig. 2 presents the estimates of period variation under the atmosphere effect per revolution for all cataloged objects with perigee altitudes from 100 to 1000 km. These estimates were determined according to catalog's data for October 2011 [16]. Decrease in the average values of this parameter with increasing altitude is explained by lowering the density of the atmosphere. At all altitudes the scattering of estimates equals several orders of magnitude. At low altitudes the scattering is smaller, than at high altitudes, due to two reasons. First, objects with large ballistic coefficients do not exist for a long time at all altitudes up to 300–400 km. Second, at altitudes higher than 600–700 km, where the atmospheric drag is relatively weak, the errors of determining parameter ΔT may reach 100% and more. Therefore, the data on drag scattering at the altitude of 400–500 km most objectively characterize the scattering of possible values of ballistic coefficients of cataloged objects. In the majority of cases this scattering does not exceed two orders of magnitude. In modeling the evolution of the distribution of SO number in the perigee altitude one should take into account the statistical distribution pðd; kb Þ of possible values of ballistic coefficients for objects of various sizes. This distribution is presented in Table 1. Possible kb values are separated

a δh ¼ ð1  eÞ Uδa  a U δe ¼  4π Uðkb U ρ UpÞ1 þ e Uexpð zÞ U

U fI 0 ðzÞ  I 1 ðzÞ þe U½I 1 ðzÞ 0:5 UI 0 ðzÞ  0:5 U I 2 ðzÞ þ :::g ð29Þ Here a and p are semi-axes and parameters of the orbit; ρ is the density of atmosphere at the perigee; z ¼ a Ue=H, I j ðzÞ are Bessel's functions of the imaginary argument of the order j, and factor ðkb ρpÞ is dimensionless. It characterizes the level of atmospheric disturbances. The product expð zÞ U f:::g takes into account the influence of the shape of an orbit. For circular orbits (e¼0) it equals unity. With increasing eccentricity its value decreases. Formula (29) is approximate [15]. In the figure brackets there are no terms proportional to e2, e3, etc. In addition, formula (29) ignores the effect of “swelling” of the atmosphere and the orbit distinction from ellipse. The magnitude

Fig. 2. Estimates of parameter del T for SOs with vatious perigee altitudes.

52

A.I. Nazarenko / Acta Astronautica 100 (2014) 47–56

Table 1 Statistical distribution pðd; kb Þ of SOs of various sizes SO size (cm)

0.1–0.25 0.25–0.5 0.5–1.0 1.0–2.5 2.5–5.0 5.0–10 10–20 420

Ballistic coefficient values (m2/kg) 0.005

0.015

0.05

0.15

0.5

1.5

0 0 0 0 0 0.059 0.157 0.05

0 0 0 0.077 0.202 0.235 0.210 0.35

0 0.08 0.272 0.308 0.267 0.235 0.210 0.40

0.43 0.35 0.364 0.308 0.267 0.235 0.210 0.15

0.43 0.43 0.272 0.230 0.200 0.176 0.157 0.05

0.14 0.14 0.091 0.077 0.066 0.060 0.056 0

Fig. 4. Normalized altitude distributions.

magnitude. The presented data on possible values of ballistic coefficients of SOs of various size are used in constructing the model of their spatial distribution. In conclusion, we note that accounting for scattering of possible values of ballistic coefficients (Table 1) is the basic factor that determines the accuracy of prediction of a spatial distribution of space debris. Naturally, particular values of the distribution pðd; kb Þ may be updated in due time. 5. Examples of prediction of a spatial distribution of space debris Fig. 3. Calculated values of parameter del T for cataloged SOs.

into 6 ranges. The average kb values lie in the range from 0.005 m2/kg to 1.5 m2/kg, i.e. the maximum possible values are 300 times the minimum ones. Here the sum of density values pðd; kb Þ in each line is equal to 1. The last line presents the data for cataloged SOs. Their application for calculating the drag characteristics (ΔT) of SOs at various altitudes lead to the results presented in Fig. 3. It is seen that these calculated data are in good agreement with the real estimates of a drag parameter at various altitudes (Fig. 2). For smaller-size (non-cataloged) SOs the scattering of possible values of ballistic coefficients was determined from statistical modeling results. Namely, 3 types of objects have been considered:  a solid sphere of diameter d,  a spherical shell of thickness δ and diameter d,  a thin panel of diameter d with shell's thickness δ.

The thickness of a spherical shell and plate was assumed to be within the limits from 0.01 cm to 0.1d. The results of calculations are presented in Table 1. All these data show that, when the size of objects decreases from 10 to 0.1 cm, the range of possible values of their ballistic coefficients also decreases by two orders of

5.1. Modeling of cataloged objects on a previous interval of time In the modeling process the altitude distributions of an annual number of SOs have been updated in such a manner that the consent of modeled and real altitude distributions of a number of cataloged SOs in various years is ensured. The change of a number of SOs in time was taken into account by using the formula pðh; t i Þnew ¼ pðhÞ0 Ukðt i Þ

ð31Þ

As a result, the following estimates were constructed: (a) the nominal altitude distribution of the annual growth p (h)0, and b) the estimates of coefficients k(ti), which were used in calculating the annual growth distributions in various years. Fig. 4 presents the normalized altitude distributions of: a) the number of SOs in the catalog at various altitudes (at the end of 2012), and b) the nominal annual growth of a number of SOs. The estimate of this growth constituted 413 objects. For each curve the sum of presented estimates is equal to 1. It is seen from the data of this figure that at altitudes up to 1000 km the character of distributions is different. Here the fraction of annual growth is essentially greater than SOs fraction at the same altitudes in the catalog. This effect is explained by the influence of SO drag in the atmosphere. At low altitudes many SOs have dropped so deeply that it terminated their lifetime.

A.I. Nazarenko / Acta Astronautica 100 (2014) 47–56

53

Table 2 Values of coefficient k(ti) Data for various years Years k(ti)

1960–1990 1.0

1990–2006 0.8

2007 5.0

2008 1.5

2009 3.0

2010 2.0

2011 0.8

2012 1.5

readers, publications [18–23], where these methodological issues are discussed in detail. Fig. 7 presents the obtained altitude distribution of a number of SOs of various areas, formed as a result of one collision of SO larger than 20 cm in size. This altitude distribution of collision fragments is taken into account in the right-hand parts of evolution equations during their integration as an additional source of space debris (SD) formation. Such technology does not require repeated predictions. In addition, it allows taking into account the small-sized fragments. Two scenarios were applied here for performing longterm predictions on the basis of SDPA model application which are as follows:

Fig. 5. Modeling of a number of catalogued objects.

The values of coefficient k(ti) are presented in Table 2. High values of coefficient k(ti) in 2007 and 2009 are explained by unique events of SO fragmentation during these years. Fig. 5 presents the results of modeling the number of cataloged SOs over the time interval since 1960 to the end 2012. Periodic variations of SOs number in prediction results are explained by the solar activity effect on their atmospheric drag. The growth of a number of objects in 2007 and 2009 is caused by unique events of fragmentation. The results presented here were obtained by using the evolution equations considered above. A 2.40 GHz personal computer was applied here. Computation time for implementing the prediction was about 1 s. Fig. 6 presents the corresponding NASA's data [16]. The data of this figure agree to an acceptable degree with the prediction results. Here the number of SOs in the catalog is slightly greater, because it includes the objects at various altitudes – up to the GEO region.

5.2. Prediction of space debris of various sizes The basic problem that arises in modeling the evolution of small-size (non-classified) space debris fragments is the formation of the initial distribution for objects of various sizes. The difficulties in solving this problem are associated with the lack of experimental data, with the impossibility of predicting emergencies, with the necessity of estimating the probability of collisions and with application of any fragmentation models. The detailed description of the technique for overcoming these difficulties in the SDPA model is beyond the scope of this paper. We refer to

Scenario 1. Total termination of all launches over the prediction interval, excluding the possibility of explosions of spacecraft and launch vehicles (the “ideal” scenario). Scenario 2. Continuation of launches at a medium rate, excluding the possibility of explosions of spacecraft and launch vehicles and without applying deflections of spacecraft and launch vehicles (the pessimistic scenario). These scenarios take into account, to a maximum extent, such measures on mitigating human-produced space contamination, as preventing space debris formation during a regular operation of a space system (SS), and preventing possible SS destruction as a result of explosions. In all scenarios the consequences of mutual collisions of cataloged SOs (larger than E20 cm in size) were taken into account. For the 1st scenario the basic prediction results are presented in Fig. 8. For 4 ranges of sizes this figure presents the estimates of a number of SOs, which exist in the space over the prediction interval, as well as the estimates of a number of collisions of SOs larger than E20 cm (cataloged objects) and larger than 10 cm in size. It is seen from these data that on the prediction interval the number of SOs larger than E20 cm in size will decrease 2-fold. This effect is explained by the dissipative effect of SO drag in the atmosphere, as a result of which the altitude of SOs drops so low that they cease orbital flight. This effect results in gradual decreasing of an annual increment of a number of collisions. Nevertheless, the number of collisions will grow monotonously, and in 200 years it will reach the value of 14. These collisions will result in growing the number of SOs of smaller size. Whereas for SOs sizing from 10 to 20 cm this growth will be insignificant (by 24%), the growth of a number of smaller-size fragments will be essential: their quantity will increase 4–6 times! The result stated here, namely, the significant growth of a number of small-size SOs over the prediction interval, as

54

A.I. Nazarenko / Acta Astronautica 100 (2014) 47–56

Fig. 6. NASA’s data on the growth of a number of SOs in the catalog.

Fig. 7. Altitude distribution of a number of fragments of single collision.

Fig. 8. Change of a number of SOs in scenario 1.

Fig. 9. Change of a number of SOs in scenario 2.

a result of collisions, seems to be quite important. Over the considered time interval it is inevitable. Even complete termination of launches of new SOs and exclusion of other sources of space debris formation cannot stop, in the foreseeable future, this process of growing a number of small-size SOs as a result of collisions. Nevertheless, the results stated above demonstrate a possibility of preventing the exponential growth of a number of small-size fragments – the tendency of stabilizing a number of collisions' fragments is observed. This is achieved on the basis of complete termination of launches of spacecraft and launch vehicles. The basic prediction results for the 2nd scenario are presented in Fig. 9. As one should expect, in this scenario the growth of a number of cataloged SOs (larger than E20 cm in size) will continue. In 200 years their quantity will increase about 1.5 times. As a consequence, the annual increment of a number of collisions will grow, and

A.I. Nazarenko / Acta Astronautica 100 (2014) 47–56

55

Table 3 Comparison of a number of catastrophic collisions and growth of a number of SOs of size 410 cm in prediction for 200 years, calculated by different models Model

SDPA

Scenario Number of collisions Enhancement of debris population

1 19 0.77

in 200 years the total number of collisions of cataloged SOs will reach the value of 48, that is 3.4 times greater than the corresponding estimate for Scenario 1. The number of SOs of 10–20 cm in size will increase 3.2 times, and the increase of a number of smaller-size fragments will be quite essential: their quantity will grow 13–20 times! A prominent feature of the data of Fig. 9 is the exponential growth of both a number of collisions and a number of small-size space debris fragments over the prediction interval. This result confirms once again the necessity of taking cardinal measures on preventing the formation of new space debris. Thus, if the intensity of growth of a number of cataloged SOs is kept at the existing level, the number of small-size, non-cataloged objects will grow exponentially. This implies that the avalanche growth of human-produced space contamination has turned from a hypothesis (the Kessler syndrome) to reality: this process has already begun. The possibility of preventing the exponential growth of a number of small-size fragments can be achieved by zeroing the increment of a number of new SOs of large size. However, under real conditions this measure may occur to be insufficient, since the technique of estimating the consequences of collisions, applied at the given stage of studies, does not take into account the “contribution” of collisions of objects smaller than 20 cm in size. The results presented here were obtained using the evolution equations considered above. The environment prediction for 100 years takes not more than 1–2 min. Let us compare the above results with the materials of report [24]. – Our scenarios (1) and (2) reflect the extreme situations that include 90% PMD scenario of the mentioned IADC report. – The data about a number of collisions and enhancement of debris population in case of fragment's size larger than 10 cm at the end of the 200-year period are presented in Table 3 below. The data were calculated using the SDPA model and the models listed in the table.

It is seen from table's data that all results of calculations by the NASA, ESA and BNCS models lie in the interval between the results of calculations according to the optimistic and pessimistic scenarios of prediction by the SDPA model. This testifies that the results of calculation by different models (with respect to estimation of a number of SOs larger than 10 cm in size) are in acceptable agreement.

2 66 2.1

NASA

ESA

BNCS

90% PMD scenario 36 1.35

50 1.40

32–40 1.0–1.35

6. Conclusion 1. When using the traditional approach the implementation of long-term predictions of space debris environment is a quite a labor-consuming computation task. Even the application of contemporary supercomputers does not allow taking into account, to a sufficient extent, the essential influencing factors, such as the consequences of mutual collisions of objects of various sizes. 2. The application of evolution equations, considered in this paper, allows several orders of magnitude lower the computer time expenses and takes into account the main influencing factors in the prediction of space debris environment.

References [1] A.I. Nazarenko Prediction and Analysis of Orbital Debris Environment Evolution. First European Conference on Space Debris, Darmstadt, Germany, April 1993. [2] G.M. Chernyavskiy and A.I. Nazarenko. Modeling of near-Earth space contamination. Collection Collisions in the near-Earth space. Space debris. Institute of astronomy, RAS, Kosmoinform,1995. [3] B. Reynolds. Documentation of Program EVOLVE: A Numerical Model to Compute Projections of the Man-Made Orbital Debris Environment, System Planning Corp. Technical Report OD91-002-U-CSP, February 1991. [4] P. Eichler, B. Reynolds. Mid- and Long-Term Debris Environment Projection using the EVOLVE and CHAIN Models. 46th International Astronautical Congress, October 2-6, 1995, Oslo, Norway. [5] N.N. Smirnov, and A.I. Nazarenko, A.B. Kiselev. LEO technogeneous contaminants evolution modeling with account of satellites collisions. IAA-00-IAA.6.4.10. [6] A.I. Nazarenko., The solution of applied problems using the space debris prediction and analysis model, in: Nickolay N. Smirnov (Ed.), Space Debris: Hazard Evaluation and Mitigation, Edited by, Taylor & Francis Inc., 2002. (Chapter 4.). [7] A.I. Nazarenko. Application of average contamination sources for the prediction of space debris environment. AAS/AIAA Space Flight Mechanics Meeting. Monterey, CA, February 1998. AAS 98-161. [8] Dolado-Perez JC et al. Introducing MEDEE – a new orbital debris evolution model. Proc. ‘6th European Conference on Space Debris’ Darmstadt, Germany, April 2013. [9] L.G. Jacchia, New Static Models of the Thermosphere and Exosphere with Empirical Temperature ProfilesSmithonian. Astrophys. Report, 313, . [10] Earth's Upper Atmosphere. Density Model for Ballistic Support of Flights of Artificial Earth Satellites GOST 25645.115-84, Publishing House for the Standards, Moscow, 1985. [11] J.M. Picone, A.E. Hedin, D.P. Drob, A.C. Aikin, NRLMSISE-00 empirical model of the atmosphere: statistical comparisons and scientific issues, J. Geophys. Res. 107 (A12) (2002) 1468, http://dx.doi.org/ 10.1029/2002JA009430. [12] D.G. King-Hele, The descent of an Earth-satellite through the upper atmosphere,, J. Brit. Interpl. Soc. 15 (1956) 314. [13] C.U Okhotsymsky, T.M. Yeneyev, G.P. Taratynova, Determination of the lifetime of the artificial Earth satellite and the study of secular disturbances, Uspekhi fiz. nauk 63 (1A) (1957).

56

A.I. Nazarenko / Acta Astronautica 100 (2014) 47–56

[14] Elyasberg P.E. Dependence of secular changes of orbital elements of the of the air drag. In: «Artificial Earth satellites», 1958, no. 1, p. 21. [15] A.I. Nazarenko, B.S. Skrebushevsky, Evolution and stability of satellite systems, MASHINOSTROYENIE (1981). [16] 〈http://www.space-track.org〉. [18] A.I. Nazarenko, “Space debris status for 200 years taking into account collisions and estimating ahead & the Kessler effect», 29th IADC 2011. [19] A.I. Nazarenko, «Estimation of the contribution of the effect of collisions of objects larger than 1 cm in area», 30-th IADC 2012.

[20] A.I. Nazarenko, and I.V. Usovik. Space debris evolution modeling with allowance for mutual collisions of objects larger than 1 cm in area. Proceedings of the sixth European conference on space debris. ESA/ESOC, Darmstadt 2013. [21] A.I. Nazarenko. The prediction of near-Earth space contamination for 200 years and the Kessler syndrome. [22] 〈http://www.satmotion.ru/en/basic〉. [23] A.I. Nazarenko, Space debris modeling, M.: IKI RAN (2013) 216. [24] Report of IADC Study, Stability of the Future Leo Environment, 50th Session of the Scientific and Technical Subcommittee – UN COPUOS. 11–22 February 2013.