Determination of the Lifetime of an Artificial Earth Satellite and Study of the Secular Perturbations of Its Orbit

D. E. Okhotsimsky, T. M. Eneev, G. P. Taratynova

One of the important questions connected with the problem of creating an artificial Earth satellite is the sufficiently reliable determination of the time it will remain in orbit. Owing to atmospheric resistance, the satellite’s energy will be dissipated and it will gradually descend.

When the motion takes place at great altitudes in the rarefied layers of the atmosphere, the resistance is small and the satellite’s time of motion may prove to be very considerable. When the motion takes place at comparatively low altitudes, of the order of 100–150 km, the satellite’s lifetime is short, and with small transverse loads the satellite may fail to complete even one full revolution.

At present there is a large number of works devoted to the problem of determining the lifetime of an artificial satellite. In this connection, a sufficiently rigorous solution is given only for a circular orbit. To estimate the time of motion of a satellite in an elliptic orbit, various approximate methods are used, employing energy considerations and based on the fact that energy losses occur mainly in the region of perigee, when the satellite approaches the Earth most closely.

The application of these methods does not give a complete solution of the problem in the general case. Moreover, as analysis shows, the use of approximate methods for determining the lifetime may in a number of cases lead to substantial errors.

This situation has made it necessary to develop a method making it possible to determine, fairly rapidly and reliably, the lifetime of a satellite for the general case of its motion. The investigation revealed the existence of universal relationships between the principal parameters of the osculating ellipse, such as the heights of perigee and apogee or the parameter and eccentricity. These relationships are valid for any satellites and depend only on the law of distribution of air density with altitude. This circumstance made it possible to reduce the complete solution of the problem of the lifetime of artificial satellites to the construction of a one-parameter family of integral curves of a first-order equation. The graph and tables included in the article make it possible to determine the lifetime rapidly by multiplying the result taken from the table or graph by a certain number depending simply on the basic parameters of the satellite. The results presented make it possible not only to compute the lifetime, but also to determine the law

changes in the orbital parameters over time for any prescribed satellite parameters and for a sufficiently wide range of initial orbital parameters.

The integration of the equations was carried out on the high-speed electronic computing machine (BESM) of the Academy of Sciences of the USSR. The integration method used made it possible to reduce the volume of computational work as much as possible while retaining the required accuracy of the results. The method is sufficiently general and can apparently be successfully applied in many cases where the problem reduces to the integration of equations whose solution is a function close to periodic, with slowly varying parameters. Thus, a similar method was developed and successfully used to investigate the perturbed motion of a satellite in a noncentral force field[^5].

Below, in § 1, the initial data adopted for the computations are presented. In § 2 the equations of motion of a satellite in osculating elements are transformed to a new independent variable convenient for the investigation being carried out. In § 3 the method for integrating the differential equations for the secular perturbations of the parameter and eccentricity of the orbit due to the action of air resistance is described. In § 4 the results obtained are presented and discussed. In § 5 the main simplifying assumptions adopted in the computation are justified, and an estimate is given of the methodological error resulting from these assumptions. An approximate investigation is also given of the secular perturbations of the orbit due to the oblateness of the Earth and the rotation of the atmosphere.

We note that the numerical results given in the paper for the lifetime of the satellite are based on the use of certain assumptions about the structure of the upper layers of the atmosphere. The absence of reliable information on the parameters of the upper atmosphere makes the numerical results suitable only for rough estimates. It should, however, be pointed out that analysis of the motion of artificial satellites will make it possible substantially to supplement the data on the upper atmosphere and will make it possible to carry out, by the method described here, the computations necessary for further refinement.

§ 1. DEPENDENCE OF ATMOSPHERIC DENSITY ON ALTITUDE

To compute the magnitude of aerodynamic resistance it is necessary to know the density of the atmosphere at high altitudes. At present there is not sufficiently precise information concerning the physical parameters of the upper atmosphere. The temperature and composition of the atmosphere are not known with a sufficient degree of accuracy. In various works, values of the density have been obtained that differ by 1.5–2.5 orders of magnitude.

In the present work, data on the physical parameters of the upper atmosphere given in[^1] were used. The dependence of density on altitude \(\rho(y)\) was approximated by formulas of the form

\[ \rho = \rho_1 \Delta, \qquad \Delta = \frac{\chi}{\left(1+\dfrac{y-y_0}{\alpha}\right)^k}, \tag{1.1} \]

where \(\rho_1\) is the air density for some fixed altitude \(y_1\). For the constants \(\chi\), \(\alpha\), \(y_0\), and \(k\), the values indicated in Table I were adopted.

Table 1

The deviations of the density values computed by formula (1.1) from the values calculated from the data\(^1\) do not exceed 10–15% of the magnitude of the density. Such an error is quite admissible, taking into account that the atmospheric data themselves are very approximate.

§ 2. EQUATIONS OF MOTION

We shall study the motion of the satellite using the osculating elements of the orbit. The equations of motion in osculating elements have the form\(^2\):

\[ \left. \begin{aligned} \frac{dp}{dt} &= \frac{2r\sqrt{p}}{\sqrt{fM}}\,T,\\[4pt] \frac{de}{dt} &= \frac{\sqrt{p}}{\sqrt{fM}}\sin\vartheta\cdot S+ \frac{\sqrt{p}}{\sqrt{fM}} \left[\left(1+\frac{r}{p}\right)\cos\vartheta+e\frac{r}{p}\right]T,\\[4pt] \frac{d\omega}{dt} &= -\frac{\sqrt{p}}{e\sqrt{fM}}\cos\vartheta\cdot S+ \frac{\sqrt{p}}{e\sqrt{fM}}\left(1+\frac{r}{p}\right)\sin\vartheta\cdot T -\frac{\sqrt{p}}{\sqrt{fM}}\ctg i\sin u\cdot W,\\[4pt] \frac{d\Omega}{dt} &= \frac{\sqrt{p}}{\sqrt{fM}}\frac{r}{p}\frac{\sin u}{\sin i}\,W,\\[4pt] \frac{di}{dt} &= \frac{\sqrt{p}}{\sqrt{fM}}\frac{r}{p}\cos u\cdot W,\\[4pt] \frac{d\tau}{dt} &= \frac{r^2}{efM}\left[-(\cos\vartheta-e\sin\vartheta\cdot N)S+\frac{p}{r}NT\right], \end{aligned} \right\} \tag{2.1} \]

where

\[ u=\omega+\vartheta,\qquad r=\frac{p}{1+e\cos\vartheta}, \]

\[ N=-2\frac{p^2}{r^2}\int_{0}^{\vartheta} \frac{\cos\vartheta\,d\vartheta}{(1+e\cos\vartheta)^3} \tag{2.1a} \]

and \(\vartheta\) is related to \(t\) by the equation

\[ t-\tau=\frac{p^{3/2}}{\sqrt{fM}} \int_{0}^{\vartheta}\frac{d\vartheta}{(1+e\cos\vartheta)^2}. \tag{2.2} \]

Here, \(p\) is the parameter of the osculating ellipse, \(e\) its eccentricity, \(\omega\) the angular distance of perigee from the node, \(\Omega\) the longitude of the ascending node,

\(i\) is the inclination of the orbit, \(\tau\) is the time of passage through the perigee of the osculating ellipse, \(\vartheta\) is the true anomaly, \(u\) is the argument of latitude; \(S, T, W\) are the projections of the perturbing acceleration on the radius vector, on the perpendicular to it in the plane of the osculating ellipse, and on the perpendicular to the plane of the osculating ellipse; \(r\) is the radius vector, \(f\) is the constant of gravitation, \(M\) is the mass of the Earth. We note that the values of \(e\) entering into the integrand in the right-hand sides of equations (2.2) and formula (2.1a) correspond to the instant of time \(t\).

In the case where the right-hand sides of system (2.1) do not depend explicitly on time, it is convenient to introduce a new independent variable—the argument of latitude \(u\). In order to pass to this variable, let us derive the differential relation connecting \(u\) with the osculating elements of the orbit and with the time \(t\). For this purpose we shall use the basic rule applied in courses of celestial mechanics in deriving the equations for osculating elements. This rule consists in the following: the kinematic elements of absolute motion can be expressed through the osculating elements of the orbit by means of formulas relating these same kinematic elements to the elements of unperturbed elliptic motion. Using this rule, one can derive the following relation:

\[ r^{2}\,d\sigma=\sqrt{fMp}\,dt, \tag{2.3} \]

where \(d\sigma\) is the total angular displacement of the radius vector during the time \(dt\). We note that for unperturbed motion relation (2.3) becomes the area integral, and in this case, obviously, \(d\sigma=d\vartheta\).

Fig. 1.

Let us find an expression for \(d\sigma\). From Fig. 1 it is not difficult to see that

\[ d\sigma=\widehat{BC'}-\widehat{BC}=\widehat{BC'}-\widehat{OC} =\widehat{O'C'}-\widehat{OC}+\widehat{BO'}, \tag{2.4} \]

where \(\widehat{OC}\) and \(\widehat{BC'}\) are projections of arcs of osculating ellipses onto the unit sphere for the instants of time \(t\) and \(t+dt\). On the other hand, we have

\[ \widehat{OO'}=d\Omega,\qquad \widehat{BO'}=\cos i\,d\Omega, \]

\[ \widehat{OC}=u,\qquad \widehat{O'C'}=u+du. \tag{2.5} \]

Substituting (2.5) and (2.4), we obtain:

\[ r^{2}(du+\cos i\,d\Omega)=\sqrt{fMp}\,dt \]

or

\[ r^{2}\left(\frac{du}{dt}+\cos i\,\frac{d\Omega}{dt}\right) =\sqrt{fM}\cdot\sqrt{p}. \tag{2.6} \]

Finally, taking into account the fourth equation of system (2.1), we obtain the following final differential relation for \(u\):

\[ \frac{du}{dt} = \sqrt{fM}\,\frac{\sqrt{p}}{r^{2}} \left( 1-\frac{r^{3}}{fMp}\,\operatorname{ctg} i\,\sin u\cdot W \right). \tag{2.7} \]

Using formula (2.7), we transform system (2.1) to the form

\[ \left. \begin{aligned} \frac{dp}{du} &=-\frac{2\gamma}{fM}\,r^3 T,\\ \frac{de}{du} &=-\frac{r^2\gamma}{fM}\left[\sin\vartheta\cdot S+\cos\vartheta\left(1+\frac{r}{p}\right)T+e\frac{r}{p}T\right],\\ \frac{d\omega}{du} &=-\frac{r^2\gamma}{fMe}\left[\cos\vartheta\cdot S+\sin\vartheta\left(1+\frac{r}{p}\right)T-e\frac{r}{p}\operatorname{ctg} i\sin u\cdot W\right],\\ \frac{d\Omega}{du} &=-\frac{r^3\gamma}{fMp}\frac{\sin u}{\sin i}\,W,\\ \frac{di}{du} &=-\frac{r^3\gamma}{fMp}\cos u\cdot W, \end{aligned} \right\} \tag{2.8} \]

where

\[ \vartheta=u-\omega,\qquad \gamma=\frac{1}{1-\dfrac{r^3}{fMp}\operatorname{ctg} i\sin u\cdot W}. \tag{2.9} \]

System (2.8) is closed, and the equations of this system must in the general case be integrated jointly. The time \(t\) can be determined, after system (2.8) has been integrated, by means of the equation

\[ \frac{dt}{du}=\frac{r^2\gamma}{\sqrt{fMp}}. \tag{2.10} \]

The time of passage through perigee \((\tau)\) can also be determined after integration of system (2.8) by means of the corresponding equation. This equation is not given here, since it will not be needed later.

Instead of the argument of latitude \(u\), some other angular parameter may be taken as the independent variable, for example the true anomaly \(\vartheta\). In this case system (2.1), after transformations, will also be reduced to the form (2.8), with the only difference that in the left-hand sides of the equations there will stand derivatives not with respect to \(u\), but with respect to \(\vartheta\). In addition, the quantity \(\gamma\) with independent variable \(\vartheta\) will be determined not by formula (2.9), but by the formula

\[ \gamma=\frac{1}{1+\dfrac{r^3}{fMe}\cos\vartheta\cdot S-\dfrac{r^3}{fMe}\left(1+\frac{r}{p}\right)\sin\vartheta\cdot T}. \tag{2.11} \]

The equations thus obtained will differ from the equations with the true anomaly as independent variable given in \({}^{2}\). The equations indicated are obtained in \({}^{2}\) from equations (2.1) with the aid of the relation

\[ r^2\,d\vartheta=\sqrt{fMp}\,dt, \tag{2.12} \]

which is valid for osculating motion, but erroneous for perturbed motion, for which, instead of (2.12), relation (2.3) must be used. Consequently, the equations with respect to the true anomaly given in \({}^{2}\) are erroneous.

Let us note that an indication of the possibility of passing to the true anomaly as the independent variable by means of relation (2.12) is also contained in \({}^{3}\).

We also note that, in studying perturbations in the first approximation, the use of relation (2.12) is admissible, since it is satisfied for unperturbed elliptic motion. We shall make use of this remark in § 5 when determining secular perturbations.

It should be pointed out that the equations in the variable \(u\) are more convenient for computations than the equations in the variable \(\vartheta\), since for small values of the eccentricity, owing to the significant changes of \(r\) that occur, the derivative \(\dfrac{d\vartheta}{dt}\) may be very small and in some cases may vanish. Such a situation occurs, for example, in the motion of a satellite in a circular orbit in the equatorial plane under the action of the gravitational force of an oblate terrestrial spheroid.

§ 3. METHOD FOR DETERMINING THE LIFETIME OF AN ARTIFICIAL SATELLITE

Consider the motion of a satellite in the terrestrial atmosphere under the condition that the Earth’s gravitational field is central. We shall also neglect the rotation of the atmosphere together with the Earth in its diurnal motion. For estimating the lifetime of the satellite, such assumptions are quite acceptable. In this case, evidently, \(\widetilde W=0\), and system (2.8) takes the form:

\[ \left. \begin{aligned} \frac{dp}{du} &= \frac{2r^{3}}{fM}\,T,\\ \frac{de}{du} &= \frac{r^{3}}{fM}\left[\sin\vartheta\cdot S+\cos\vartheta\left(1+\frac{r}{p}\right)T+e\frac{r}{p}T\right],\\ \frac{d\omega}{du} &= \frac{r^{2}}{fMe}\left[\cos\vartheta\cdot S+\sin\vartheta\left(1+\frac{r}{p}\right)T\right],\\ \frac{d\Omega}{du} &= 0,\\ \frac{di}{du} &= 0,\\ \vartheta &= u-\omega. \end{aligned} \right\} \tag{3.1} \]

From (3.1) we see that

\[ \Omega=\Omega_{0}=\mathrm{const}, \qquad i=i_{0}=\mathrm{const}, \]

i.e., the resistance of the atmosphere does not cause secular perturbations of the longitude of the node or of the inclination of the orbit. An estimate of the effect of the rotation of the atmosphere will be given below, in § 5.

For the acceleration due to the resistance force we take:

\[ R_{x}=\frac{c_{x}F}{m}\,\frac{\rho v^{2}}{2}, \]

where \(m\) is the mass of the satellite, \(v\) is the velocity of the satellite relative to the air, \(c_x\) is the aerodynamic drag coefficient, and \(F\) is the area to which the aerodynamic coefficient \(C_x\) is referred. We have:

\[ \left. \begin{aligned} S&=-\frac{c_x F}{m}\,\frac{\rho v}{2}\,v_r,\\ T&=-\frac{c_x F}{m}\,\frac{\rho v}{2}\,v_n. \end{aligned} \right\} \tag{3.2} \]

Here \(v_r\) and \(v_n\) are the radial and transverse components of the velocity:

\[ \left. \begin{aligned} v_r &= \sqrt{\frac{fM}{p}}\, e \sin \vartheta,\\ v_n &= \sqrt{\frac{fM}{p}}\,(1+e\cos\vartheta), \end{aligned} \right\} \tag{3,3} \]

Substituting (1,1), (3,2), and (3,3) into (3,1) and using the formulas of elliptic theory:

\[ r=\frac{p}{1+e\cos\vartheta}, \]

\[ v=\sqrt{\frac{fM}{p}}\sqrt{1+2e\cos\vartheta+e^2}, \]

we obtain:

\[ \left. \begin{aligned} \frac{dp}{du} &= -c\rho_1 \varphi(p,e,\omega,u),\\ \frac{de}{du} &= -c\rho_1 \psi(p,e,\omega,u),\\ \frac{d\omega}{du} &= -c\rho_1 \chi(p,e,\omega,u), \end{aligned} \right\} \tag{3,4} \]

where

\[ \left. \begin{aligned} \varphi &= \frac{p^3\Delta\sqrt{1+2e\cos\vartheta+e^2}} {(1+e\cos\vartheta)^2},\\ \psi &= \frac{p\Delta\sqrt{1+2e\cos\vartheta+e^2}} {(1+e\cos\vartheta)^2},\\ \chi &= \frac{p\Delta\sin\vartheta\sqrt{1+2e\cos\vartheta+e^2}} {2e(1+e\cos\vartheta)}; \end{aligned} \right\} \tag{3,5} \]

\(c\) is a constant, \(c=\dfrac{c_x F}{m}\), and \(\vartheta=u-\omega\).

Let us note that direct integration of equations (3,4) is impractical, since, owing to the large interval of integration and the limited step in \(u\), the number of steps would be very large, which would lead to a large amount of time required for integrating the system, and could also lead to a substantial accumulation of error in the sought functions. To avoid this, we proceed as follows.

Let us integrate equations (3,5) with respect to \(u\) from 0 to \(2\pi\). We obtain:

\[ \Delta p=-c\rho_1\int_0^{2\pi}\varphi(p,e,\omega,u)\,du, \tag{3,6} \]

\[ \Delta e=-c\rho_1\int_0^{2\pi}\psi(p,e,\omega,u)\,du, \tag{3,7} \]

\[ \Delta\omega=-c\rho_1\int_0^{2\pi}\chi(p,e,\omega,u)\,du, \tag{3,8} \]

where \(\Delta p\), \(\Delta e\), \(\Delta\omega\) are the changes of the parameter, the eccentricity, and the distance of perigee from the node over one revolution. Let us note that in computing the indicated

of the integrals, owing to the small change of the parameter \(p\), the eccentricity \(e\), and the distance of the perigee from the node \(\omega\) over the course of one revolution, these elements may be taken as constant. In this case it is easy to show that the integral appearing on the right-hand side of (3.8) vanishes:

\[ \Delta \omega = 0. \tag{3.9} \]

Let us note that for very small eccentricities \((e < 0.0001)\) the assumption adopted may prove to be quite crude for \(\omega\); however, in this case, for determining the lifetime of the satellite, only the change of the parameter \(p\) will be significant, which, as can be shown with the aid of the first equation of the system (3.1), will practically not depend on \(\omega\).

Since over one revolution the elements \(p\) and \(e\) change very little, it may be assumed with high accuracy that the indicated changes of these quantities over one revolution are equal to the derivatives of these elements with respect to the number of revolutions of the satellite \(N\); \(N = \dfrac{u}{2\pi}\). Taking (3.9) into account and putting

\[ \omega = \omega_0 = \mathrm{const}, \]

we obtain

\[ \left. \begin{aligned} \frac{dp}{dN} &= -c \rho_1 \int_{0}^{2\pi} \varphi(p,e,u)\,du, \\[6pt] \frac{de}{dN} &= -c \rho_1 \int_{0}^{2\pi} \psi(p,e,u)\,du. \end{aligned} \right\} \tag{3.10} \]

Putting \(\nu = Nc\), from equations (3.10) we obtain:

\[ \left. \begin{aligned} \frac{de}{dp} &= \frac{\displaystyle \int_{0}^{2\pi} \psi(p,e,u)\,du} {\displaystyle \int_{0}^{2\pi} \varphi(p,e,u)\,du}, \\[10pt] \frac{d\nu}{dp} &= -\frac{1} {\displaystyle \rho_1 \int_{0}^{2\pi} \varphi(p,e,u)\,du}. \end{aligned} \right\} \tag{3.11} \]

Thus, the problem of determining the lifetime of a satellite has been reduced to integrating a system of two differential equations (3.11), whose right-hand sides are expressed through definite integrals of the quantities \(\varphi\) and \(\psi\), which, according to formulas (3.5), are known functions of \(p\), \(e\), and \(u\). We also note that the right-hand sides of equations (3.11) do not depend on the aerodynamic coefficient \(c_x\), nor on the design parameters of the satellite: the weight \(G\) and the midsection area \(F\). Therefore, for a given atmosphere it is sufficient to integrate equations (3.11) once, and then, by a simple transition from \(\nu\) to \(N\) \(\left(N=\dfrac{\nu}{c},\ \text{where } c=\dfrac{c_x F}{G}\,g\right)\), obtain the number of revolutions of the satellite for specified values of the quantities \(c_x\), \(G\), and \(F\). Hence, assuming that one revolution around the Earth is completed by the satellite in approximately 90–100 minutes, one can obtain a sufficiently accurate estimate of its lifetime.

When integrating system (3,11), for convenience in solving it and to economize on machine time, the integrals, the functions \(\cos\), and the square root in the right-hand sides of the equations were introduced by means of differential equations. Finally, the basic (3,11) and auxiliary (3,13) systems of equations, programmed for solution on the machine, had the form

\[ \frac{de}{dp}=\frac{\zeta_k}{\eta_k},\qquad \frac{dv}{dp}=\frac{1}{\rho_1\eta_k}, \tag{3,12} \]

where

\[ \left. \begin{aligned} \frac{d\eta}{du} &=-\,\frac{\xi p^3\Delta}{(1+e\lambda_2)^3},& \frac{d\zeta}{du} &=-\,\frac{\xi p\Delta(e+\lambda_2)}{(1+e\lambda_2)^2},\\ \frac{d\xi}{du} &=-\,\frac{e\lambda_1}{\xi},& \frac{d\lambda_1}{du}&=\lambda_2,& \frac{d\lambda_2}{du}&=-\lambda_1 . \end{aligned} \right\} \tag{3,13} \]

In equations (3,12) and (3,13) the following notation is adopted:

\[ \lambda_2=\cos \vartheta,\qquad \lambda_1=\sin \vartheta,\qquad \xi=\sqrt{1+e^2+2e\lambda_2}, \]

\[ \eta=-\int_0^u \frac{\xi p^3\Delta}{(1+e\lambda_2)^3}\,du, \]

\[ \zeta=-\int_0^u \frac{\xi p\Delta(e+\lambda_2)}{(1+e\lambda_2)^2}\,du,\qquad \zeta_k=\zeta_{u=2\pi},\qquad \eta_k=\eta_{u=2\pi}. \]

\(\Delta\) is determined by formula (1,1). For \(\rho_1\) the density of the atmosphere at an altitude of \(100\) km was taken. In addition to the quantities \(p\), \(e\), \(v\), at each computational point of the basic system (3,12) the velocity at perigee \(v_\pi\), the apogee height \(h_a\), and the perigee height \(h_\pi\) were additionally computed from the formulas

\[ v_\pi=\sqrt{fM}\,\frac{1+e}{\sqrt{p}},\qquad h_a=\frac{p}{1-e}-R,\qquad h_\pi=\frac{p}{1+e}-R, \]

where \(R\) is the radius of the Earth.

The integration was carried out as follows. For prescribed values of \(p\) and \(e\), the auxiliary system (3,13) was integrated with respect to \(u\) over the interval from \(0\) to \(2\pi\) by the Runge—Kutta method with a constant step. Since, for the purposes of subsequent interpolation, it was necessary as a result of the calculations to obtain a sufficiently dense grid of values of the quantities \(p\) and \(e\), and since the required accuracy of the calculation was comparatively low, it proved expedient for integrating the basic system (3,12) to use Euler’s method with a constant step. For integrating the auxiliary system (3,13) the step \(\Delta\vartheta=6^\circ\) was chosen. In integrating the basic system (3,12) the step \(\Delta p=5\) km was chosen for values of the apogee height \(h_a\leq 700\) km, and \(\Delta p=10\) km for values \(h_a>700\) km. For the adopted step values, the methodological error in determining the time of motion of the satellite due to approximate integration did not exceed \(2\text{–}5\%\).

As the initial data for system (3,12), it proved more convenient to prescribe the initial values of the apogee height \(h_{a_0}\) and the perigee height \(h_{\pi_0}\), with which the initial values of the eccentricity and of the parameter are related by

\[ e_0=\frac{h_{a_0}-h_{\pi_0}}{h_{a_0}+h_{\pi_0}+2R},\qquad p_0=(h_{a_0}+R)(1-e_0). \]

Since each point on the integral curve can be regarded as the initial point for some other motion, in choosing the initial

D. E. OKHOTSIMSKII, T. M. ENEEV, G. P. TARATYNOVA

Fig. 2.

values of the parameters \(h_a\) and \(h_\pi\), there was no need to vary both parameters. It was sufficient to obtain, for example, a series of integral curves in \(h_\pi\) for the maximum initial value of \(h_a\). The computations were carried out for an initial apogee height \(h_{a_0}=1600\) km and initial perigee heights from the range

\[ 160\ \text{km} \leq h_{\pi_0} \leq 500\ \text{km}. \]

In view of the fact that, for heights greater than 900 km, there are no data concerning the density of the atmosphere, it was assumed that the law of variation of the density at such heights is the same as for heights

\[ 250\ \text{km} \leq y \leq 900\ \text{km}. \]

Integration of the system of equations (3,12) was carried out up to the moment when the satellite reached a height of 100 km. Consideration of the satellite’s motion after this was meaningless, since, owing to strong atmospheric braking, it could exist for only a very short time.

§ 4. RESULTS OF THE COMPUTATIONS AND THEIR DISCUSSION

The results of the computations are given in Table II and in Fig. 2. Table II gives the values of the quantity \(\nu\) \([\text{m}^3/\text{kg sec}^2]\) as a function of the initial values of the apogee height \(h_a\) and the perigee height \(h_\pi\). The same table gives the values of the velocity at perigee \(v_\pi\) [in m/sec] at the beginning of the satellite’s motion. With the aid of the results given in Table II, and using the relation

\[ N=\nu\,\frac{G}{F}\,\frac{1}{g c_x}, \tag{4,1} \]

where \(g=9.81\ \text{m}/\text{sec}^2\), one can determine the number of revolutions of the satellite and, consequently, the time of its existence for any values of the aerodynamic coefficient \(c_x\) and the transverse loading \(\dfrac{G}{F}\) [in kg/m\(^2\)].

In Fig. 2, two families of curves are shown in the plane \((h_a, h_\pi)\). The first family corresponds to the dependences

\[ h_\pi=f(h_a), \tag{4,2} \]

which occur in the motion of the satellite around the Earth. Each point of the indicated curves corresponds to its own value of \(\nu\). The lines of the second family satisfy the relation

\[ \nu=\text{const} \tag{4,3} \]

and connect, on the curves of the first family, points corresponding to the same values of \(\nu\). Joint consideration of both families of curves makes it possible not only to estimate the total lifetime of the satellite for certain initial values of the apogee and perigee heights, but also to estimate the time during which the apogee and perigee heights vary within certain limits.

On the basis of the computational results presented, a number of conclusions may be drawn concerning the character of the change in the orbital parameters during the satellite’s motion. It is seen that the apogee and perigee heights decrease monotonically, and for all elliptical orbits the rate of decrease of the apogee height is greater than the rate of decrease of the perigee height. For highly elongated orbits this difference may be quite significant. Thus, for an orbit with perigee height 300 km and apogee height 700 km, a lowering of the apogee by 100 km corresponds to a lowering of the perigee by approximately 6 km. For larger

values of the apogee height this difference will be still more substantial. Thus, when there is a large difference between the apogee and perigee heights, the change in the orbital parameters will, over a long period of time, amount practically only to a decrease in the apogee height at an almost constant perigee height.

With such a change in the shape of the orbit, the eccentricity of the orbit will decrease all the time and tend to zero. The satellite’s orbit will

Table II

Values of the quantity \(\gamma\) and of the velocity at perigee \(v_\pi\) as functions of the perigee height \(h_\pi\) and the apogee height \(h_a\)

LIFETIME OF A SATELLITE AND SECULAR PERTURBATIONS OF ITS ORBIT

Table 11 (continued)

Table II (continued)

tend toward circular. In Fig. 2 this is manifested in the fact that the curves of the first family tend, as to an asymptote, to the straight line \(h_a = h_\pi\), corresponding to a circular orbit. Examination of Fig. 2 shows that, as the eccentricity decreases, the difference between the rates of decrease of apogee and perigee diminishes.

The results of the calculations make it possible to estimate the change in the satellite’s time of motion when the initial parameters \(h_{a0}\) and \(h_{\pi0}\) are changed, and also to indicate the values of the initial parameters that make it possible to ensure a prescribed time of motion in the simplest way.

The character of the curves of the second family shows that, as was to be expected, the lifetime of the satellite increases more strongly with an increase in the initial height of perigee and more weakly with an increase in the initial height of apogee. Thus, for example, for an orbit with perigee 360 km and apogee 1500 km, a change of the perigee height by 20 km causes a change in lifetime of about 40%, while the same change of apogee causes about 2%, i.e., 20 times less.

In view of the fact that the realization of satellite orbits with a large initial perigee height may encounter a number of difficulties, it is essential to note that a considerable increase in the duration of a satellite’s existence can also be achieved with an unchanged perigee height by increasing the initial height of apogee; moreover, this requires a comparatively small increase of the velocity at perigee. Thus, for example, for an orbit with parameters \(h_\pi = 360\) km and \(h_a = 700\) km, increasing the apogee height to 1000 km leads to an increase in the lifetime by a factor of 2.2; in this case an increase of the velocity at perigee by only 78 m/sec is required. The result given indicates the advisability of using elongated orbits, which may make it possible to obtain a substantial increase in the duration of existence of an artificial Earth satellite in a comparatively simple way.

As has already been indicated, the cited values of \(\nu\) were obtained for a definite scheme of the distribution of air density with height. If the actual values deviate from those adopted, the duration of the satellite’s existence will be different. Refinement of data on the density of the upper layers of the atmosphere will make it possible, by the same method, to carry out a refined calculation of the values of \(\nu\) for issuing more accurate predictions of the satellite’s time of existence.

Since the rate of change of the orbital parameters is proportional to the density of the atmosphere, which decreases rapidly with height, the change of the orbital parameters with time will at first be considerably slower than subsequently, when descending into denser layers. Therefore the satellite will spend the greater part of its time in the high layers of the atmosphere.

The rapid decrease of density with height and the slow change of the perigee height indicate that the principal significance for the lifetime of a satellite will be the value of the air density in the region of the initial perigee. This conclusion makes it possible to estimate the magnitude of the change in the calculated lifetime when the atmospheric-density data are changed. The duration of a satellite’s existence will be approximately inversely proportional to the air density in the region of the initial perigee.

This circumstance indicates the possibility of estimating the actual air density in the region of the initial perigee from the initial values of apogee and perigee and from the actual time of existence of the satellite in orbit. Processing the results for a series of launches with different initial perigee heights will make it possible to estimate the actual distribution of air with height.

To conclude, let us give several examples. Consider a spherical satellite with diameter \(d=0.5\ \text{m}\) and weight \(G=10\ \text{kg}\). We shall take the value of the drag coefficient to be \(c_x=2\). In this case the value of the factor by which the quantity \(\nu\) must be multiplied in order to obtain the satellite’s lifetime in revolutions will be \([\text{in } \text{kg}\,\text{s}\,\text{m}^{-3}]\) equal to

\[ \frac{4G}{\pi d^2}\frac{1}{g c_x}\simeq 2.6. \]

We obtain an estimate of the lifetime in days by dividing the resulting number of revolutions by 16, which corresponds to the number of revolutions per day for a one-and-a-half-hour orbit \(\left(n\simeq \frac{N}{16}\right)\).

Let the perigee and apogee heights be \(h_\pi=360\ \text{km}\) and \(h_\alpha=800\ \text{km}\). Then from Table II we find \(\nu=2760\); the number of revolutions \(N\) will be \(N=2760\times 2.6\simeq 7200\) revolutions, and the number of days \(n\simeq 450\) days, which corresponds to approximately 1 year and 3 months.

For perigee and apogee heights equal to \(h_\pi=500\ \text{km}\) and \(h_\alpha=1500\ \text{km}\), we similarly have \(\nu=66\,600\), \(N=174\,000\) revolutions and \(n\simeq 11\,000\) days, which corresponds to a lifetime of the order of 30 years.

Thus, when sufficiently large initial perigee and apogee heights are chosen, the lifetime of an artificial Earth satellite may prove to be quite considerable.

Taking the perigee and apogee heights equal to \(h_\pi=200\ \text{km}\) and \(h_\alpha=400\ \text{km}\), we obtain \(\nu=16.8\), \(N\simeq 44\), and \(n=2.7\) days; i.e., for relatively small values of the perigee and apogee heights the duration of the satellite’s existence in orbit turns out to be small.

§ 5. On Secular Perturbations of the Orbital Parameters of an Artificial Satellite

The investigation carried out above of the question of the lifetime of an artificial satellite was based on the study of secular perturbations of the orbital elements due to atmospheric drag. In the present paragraph an estimate will be given of the secular perturbations of the orbital elements due to the influence of other perturbing factors. Periodic perturbations of the orbital elements are not considered. It is clear that periodic perturbations cannot have any appreciable influence on the lifetime of an artificial satellite.

We shall proceed from the equations of motion in the parameters of the osculating ellipse (2,8). In taking into account the influence of the deviation of the gravitational field from a central one, we shall start from the expression for the potential

\[ V=\frac{fM}{r}-\frac{\varepsilon}{3r^3}\left(3\sin^2\psi-1\right), \tag{5,1} \]

where

\[ \varepsilon=fMa^2\left(\alpha-\frac{m}{2}\right),\qquad m=\frac{\Omega^2 a}{g_a}, \tag{5,2} \]

\(\psi\) is the latitude, \(a\) is the equatorial radius of the Earth, \(\alpha\) is the flattening of the Earth, \(\Omega\) is the angular velocity of the Earth’s diurnal rotation, and \(g_a\) is the acceleration of terrestrial gravity at the equator. Formula (5,1) leads to the following expressions for the projections of the perturbing acceleration:

\[ \left. \begin{aligned} S&=-\frac{\varepsilon}{r^4}\left[3\sin^2 i\sin^2 u-1\right],\\ T&=-\frac{\varepsilon}{r^4}\sin^2 i\sin 2u,\qquad W=-\frac{\varepsilon}{r^4}\sin 2i\sin u. \end{aligned} \right\} \tag{5,3} \]

In estimating the additional accelerations caused by the effect of wind due to the rotation of the atmosphere together with the Earth, we shall proceed from

approximate formulas obtained from more exact ones by discarding terms containing the parameter \(\dfrac{\Omega r}{v}\) in powers higher than the first\(^5\):

\[ \begin{aligned} S&=\frac{c_rF}{G}-\frac{g}{2}\,\frac{v_n}{v}\,v\rho\Omega r\cos i, \qquad T=-\frac{c_xF}{G}-\frac{g}{2}\,\frac{v^2+v_n^2}{v}\,\rho\Omega r\cos i,\\ W&=-\frac{c_rF}{G}-\frac{g}{2}\,v\rho r\Omega\sin i\cos u. \end{aligned} \tag{5,4} \]

To estimate the secular perturbations, we shall integrate the right-hand sides of equations (2,8) over one revolution in \(u\) from \(u=0\) to \(u=2\pi\), considering the parameters of the ellipse constant and neglecting periodic perturbations of the parameters. Let us note that, in obtaining the indicated estimates in the first approximation, one may take \(\gamma=1\).

Using this method, we obtain that the oblateness of the Earth and the associated deviation of the potential cause secular perturbations of the longitude of the ascending node \(\Omega\) and of the angular distance of the perigee from the node \(\omega\). The other orbital parameters, including the parameter and the eccentricity, do not experience secular perturbations from the influence of the oblateness. This means that neglecting the oblateness in calculating the lifetime of the satellite is quite legitimate. The regression of the node shows that, under the influence of the oblateness, the plane of the orbit rotates about the axis of rotation of the Earth, maintaining a constant angle with this axis.

To calculate the magnitude of the regression of the ascending node and the magnitude of the displacement of the distance of the perigee from the ascending node during one revolution, we obtain the formulas:

\[ \frac{d\Omega}{dN}=-\frac{2\pi\varepsilon}{p^2 fM}\cos i, \tag{5,5} \]

\[ \frac{d\omega}{dN}=\frac{\pi\varepsilon}{p^2 fM}\left(5\cos^2 i-1\right). \tag{5,6} \]

Let us consider, as an example, an orbit with a mean altitude of the order of 500 km. We obtain for the regression of the node and the perigee

\[ \frac{d\Omega}{dN}\simeq -0^\circ,54\cos i,\qquad \frac{d\omega}{dN}\simeq 0^\circ,27\left(5\cos^2 i-1\right). \]

For an inclination of \(45^\circ\) we obtain, per revolution, \(\Delta\Omega_{\mathrm{rev}}\simeq -0^\circ,38\), \(\Delta\omega_{\mathrm{rev}}\simeq 0^\circ,4\), and per day \(\Delta\Omega_{\mathrm{day}}\simeq -6^\circ,1\), \(\Delta\omega_{\mathrm{day}}\simeq 6^\circ,5\).

Formula (5,5) shows that the rate of regression of the ascending node depends substantially on latitude and will be greatest for orbits close to equatorial ones, and equal to zero for an orbit passing through the poles.

Let us now consider the influence on the secular perturbations of the satellite’s orbital parameters of the rotation of the atmosphere together with the Earth in its diurnal motion. Using formulas (5,4), substituting \(\rho\) from formula (1,1) and integrating with respect to \(u\) from \(\omega\) to \(\omega+2\pi\), we obtain the following estimate formulas:

\[ \left|\frac{dp}{dN}\right|< \frac{2cp^3\sqrt{p}}{V\sqrt{fM}}\,\rho_1\Omega I\cos i, \tag{5,7} \]

\[ \left|\frac{de}{dN}\right|< \frac{cp^2\sqrt{p}}{2V\sqrt{fM}}\,\rho_1\Omega I\cos i, \tag{5,8} \]

\[ \left|\frac{d\Omega}{dN}\right|< \frac{cp^3\sqrt{p}}{2V\sqrt{fM}}\,\rho_1\Omega I\sin 2\omega, \tag{5,9} \]

\[ \left|\frac{di}{dN}\right|< \frac{cp^3\sqrt{p}}{2V\sqrt{fM}}\,\rho_1\Omega I\sin i, \tag{5,10} \]

\[ \left|\frac{d\omega}{dN}\right|< \frac{cp^3\sqrt{p}}{2V\sqrt{fM}}\,\rho_1\Omega I\cos i\sin 2\omega, \tag{5,11} \]

where

\[ c=\frac{C_xF}{G}\,g. \tag{5,12} \]

As the calculation shows, the perturbations of the parameter and eccentricity due to the rotation of the atmosphere do not exceed 10–12% of the corresponding perturbations for a stationary atmosphere. This is also understandable, since the change in dynamic pressure in passing from absolute motion to relative motion should not exceed this amount. The effect will be greatest for an equatorial orbit; moreover, when the satellite moves eastward the drag force will decrease, and when it moves westward it will increase in comparison with the drag for a stationary Earth. For orbits passing near the poles, the influence of the Earth’s rotation on the parameter and eccentricity is insignificant.

The result obtained means that, for the eastward motion of a satellite, the satellite’s lifetime, for the same initial values of apogee and perigee, will be greater than for a satellite launched westward. The difference in lifetime from the case in which the rotation of the atmosphere is neglected, greatest for orbits close to equatorial ones, will not exceed 10–12%. For other orbits the difference will be smaller. For a polar orbit the difference will be zero. At present, given the very low accuracy with which the density of the upper layers of the atmosphere is known, the above-mentioned magnitudes of the errors in determining the lifetime should be regarded as negligibly small.

The magnitudes of the secular perturbations of the longitude of the node, the inclination, and the distance of perigee from the node due to the rotation of the atmosphere turn out to be very small. To estimate these perturbations it is convenient to estimate the magnitude of the total displacement of the indicated parameters over the entire lifetime of the satellite. Such an estimate is possible because both the indicated perturbations and the secular perturbations of the orbit that directly affect the lifetime are proportional to the same quantities—the atmospheric density, the reciprocal of the transverse loading, and the coefficient of aerodynamic drag. In view of this, the estimate obtained will be valid for any satellites with different design parameters and different initial values of the orbital parameters.

On the basis of the calculations we obtain the following estimates:

\[ |\Delta\Omega|<0.2^\circ,\quad |\Delta i|<0.1^\circ,\quad |\Delta\omega|<0.2^\circ . \tag{5.13} \]

The results presented show that the deviations obtained are very small. The influence of the rotation of the atmosphere on the secular displacements of the longitude of the node, the inclination, and the angular distance from the node is insignificant.

We note that, in carrying out the computations, values of the integral \(I\) obtained from an analysis of the results of machine calculations were used.

The results presented above show that the method used in the present work for calculating the lifetime is sufficiently well justified; this was also confirmed by calculations carried out using the exact equations of motion in osculating elements\(^5\). This means that the accuracy of the method set forth is quite sufficient for obtaining reliable forecasts of the lifetime of artificial Earth satellites.

References

1. S. K. Mitra, The Upper Atmosphere, Moscow, IL, 1955.
2. G. N. Duboshin, Introduction to Celestial Mechanics, Moscow–Leningrad, ONTI, 1938.
3. M. F. Subbotin, A Course in Celestial Mechanics, Moscow–Leningrad, ONTI, 1937.
4. N. I. Idelson, Theory of the Potential with Applications to the Theory of the Figure of the Earth and Geophysics, Moscow, ONTI, 1936.
5. G. P. Taratynova, Article in the present issue of UFN.