On the Motion of an Artificial Satellite in the Noncentral Gravitational Field of the Earth in the Presence of Atmospheric Resistance

G. P. Taratynova

The motion of an artificial satellite of the Earth is very complex, and an analytical investigation of motion of this type, taking into account atmospheric resistance and the noncentrality of the Earth’s gravitational field, presents considerable difficulties. Until recently, most works have been limited to special cases, such as the motion of satellites in circular orbits, analytical investigations of perturbations of satellite orbits caused only by atmospheric resistance, without taking into account the flattening of the Earth, etc. The solution of the problem in the general case, without the use of high-speed electronic machines, is extremely difficult. In view of this, it seemed advisable to create a method for computing, on an electronic machine, the orbit of an artificial satellite of the Earth with simultaneous allowance for the principal perturbing forces acting on it.

In the present article a method is set forth for computing, on a high-speed digital machine, the motion of an artificial satellite of the Earth with simultaneous allowance for atmospheric resistance, under the condition that the atmosphere moves together with the Earth, and for the deviation of the gravitational field from a central one. Since the perturbing action of the Sun and the Moon on the motion of the satellite relative to the Earth turns out to be very small, we shall neglect these perturbations.

§ 1. Equation of the Perturbed Motion of an Artificial Satellite of the Earth. Method of Integration

Let us consider the motion of an artificial satellite of the Earth in the noncentral field of terrestrial gravitation in the presence of atmospheric resistance. To describe the motion we shall use the differential equations in osculating elements1:

\[ \left. \begin{aligned} \frac{dp}{dt} &= 2r\widetilde{T},\\[4pt] \frac{de}{dt} &= \widetilde{S}\sin \vartheta +\left(1+\frac{r}{p}\right)\widetilde{T}\cos \vartheta +e\frac{r}{p}\widetilde{T},\\[4pt] \frac{d\omega}{dt} &= -\frac{1}{e}\widetilde{S}\cos \vartheta +\frac{1}{e}\left(1+\frac{r}{p}\right)\widetilde{T}\sin \vartheta -\frac{r}{p}\operatorname{ctg} i\,\widetilde{W}\sin u,\\[4pt] \frac{d\Omega}{dt} &= \frac{r}{p}\frac{\sin u}{\sin i}\widetilde{W},\\[4pt] \frac{di}{dt} &= \frac{r}{p}\widetilde{W}\cos u, \end{aligned} \right\} \tag{1} \]

where (Fig. 1)

\[ u=\omega+\vartheta,\quad r=\frac{p}{1+e\cos\vartheta},\quad \widetilde S=\frac{\sqrt p}{\sqrt{fM}}S,\quad \widetilde T=\frac{\sqrt p}{\sqrt{fM}}T,\quad \widetilde W=\frac{\sqrt p}{\sqrt{fM}}W, \]

\(f\) is the gravitational constant, \(M\) is the mass of the Earth, \(\vartheta\) is the true anomaly, \(p\) is the orbital parameter, \(e\) is the eccentricity, \(\omega\) is the angular distance of the perigee from the ascending node, \(\Omega\) is the longitude of the ascending node, \(i\) is the inclination of the orbit, \(S\), \(T\), and \(W\) are the projections of the perturbing acceleration onto the radius vector, onto the perpendicular to it in the plane of the orbit, and onto the normal to the plane of the orbit, respectively.

Fig. 1.

In order that the system of equations (1), describing the motion of the satellite, be closed, an additional differential equation has been introduced, determining the dependence of the true anomaly \(\vartheta\) on time:

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

where the derivatives \(\dfrac{d\omega}{dt}\) and \(\dfrac{d\Omega}{dt}\) are determined by the corresponding equations of system (1)\(^4\).

Let us transform the system of differential equations (1) to the independent variable \(\vartheta\), the true anomaly. On the basis of relation (2) we have

\[ \frac{dt}{d\vartheta} = \frac{1}{ \dfrac{\sqrt{pfM}}{r^2} -\dfrac{d\omega}{dt} -\cos i\,\dfrac{d\Omega}{dt} }. \tag{3} \]

Finally we obtain

\[ \begin{aligned} \frac{dp}{d\vartheta} & = F\cdot 2r\widetilde T,\\ \frac{de}{d\vartheta} & = F\left[\widetilde S\sin\vartheta+ \left(1+\frac{r}{p}\right)\widetilde T\cos\vartheta +e\,\frac{r}{p}\widetilde T\right],\\ \frac{d\omega}{d\vartheta} & = F\left[-\frac{1}{e}\widetilde S\cos\vartheta +\frac{1}{e}\left(1+\frac{r}{p}\right)\widetilde T\sin\vartheta -\frac{r}{p}\operatorname{ctg} i\,\widetilde W\sin u\right],\\ \frac{d\Omega}{d\vartheta} & = F\cdot\frac{r}{p}\,\frac{\sin u}{\sin i}\,\widetilde W,\\ \frac{di}{d\vartheta} & = F\cdot\frac{r}{p}\,\widetilde W\cos u, \end{aligned} \tag{4} \]

where

\[ F=\frac{dt}{d\vartheta}. \]

The system of equations (3), (4) completely determines the motion of the satellite.

Projections of the perturbing acceleration shall be represented in the form

\[ \begin{aligned} S &= S_1 + S_2,\\ T &= T_1 + T_2,\\ W &= W_1 + W_2, \end{aligned} \tag{5} \]

where \(S_1, T_1, W_1\) are the projections of the perturbing acceleration due to the noncentrality of the Earth’s gravitation; \(S_2, T_2, W_2\) are the projections of the perturbing acceleration produced by the force of atmospheric resistance (taking into account the rotation of the atmosphere together with the Earth).

To an accuracy up to terms of first order of smallness relative to the flattening, for the potential of the Earth’s gravitation we have³

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

where

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

\(\psi\) is the geocentric latitude of the satellite, \(a\) is the equatorial radius of the Earth, \(\alpha\) is the flattening, \(\alpha = \dfrac{a-b}{a}\), \(\Omega\) is the angular velocity of the Earth’s rotation, and \(g_a\) is the acceleration of the Earth’s gravitational force at the equator.

Differentiating expression (6) for the potential along the radius and along the meridian in the direction toward the north, we obtain the formulas for the projections of the acceleration of the Earth’s gravitational force on the radius vector and on the tangent to the meridian, respectively,

\[ \begin{aligned} g_r &= -\frac{fM}{r^3} + \frac{\varepsilon}{r^4}\left(3\sin^2\psi - 1\right),\\ g_m &= -\frac{\varepsilon}{r^4}\sin 2\psi. \end{aligned} \tag{7} \]

Hence, for the projections of the perturbing acceleration \(S\), \(T\), and \(W\) due to the deviation of the Earth’s gravitational field from a central one, taking into account terms of first order relative to \(\alpha\), we obtain

\[ \begin{aligned} S_1 &= \frac{\varepsilon}{r^4}\left(3\sin^2 i \sin^2 u - 1\right),\\ T_1 &= -\frac{\varepsilon}{r^4}\sin^2 i \sin 2u,\\ W_1 &= -\frac{\varepsilon}{r^4}\sin 2i \sin u. \end{aligned} \tag{8} \]

Let us derive formulas for the projections \(S_2\), \(T_2\), and \(W_2\). For the acceleration caused by the force of atmospheric resistance \(R\), we take

\[ R = \frac{c_x s}{m}\,\frac{\rho v_{\mathrm{rel}}^2}{2}, \tag{9} \]

where \(m\) is the mass of the satellite, \(v_{\mathrm{rel}}\) is the velocity of the satellite relative to the air, \(c_x\) is the coefficient of aerodynamic resistance, \(s\) is the area to which the aerodynamic coefficient is referred, and \(\rho\) is the density of the atmosphere.

For the vector of relative velocity \(\mathbf{v}_{\mathrm{rel}}\) we have

\[ \mathbf{v}_{\mathrm{rel}} = \mathbf{v} - \mathbf{v}_1, \tag{10} \]

where

\[ v_1 = \Omega r \cos\psi, \]

where \(\mathbf v_1\) is directed from west to east (Fig. 2); \(\mathbf v\) is the velocity of the satellite relative to fixed space.

Projecting \(\mathbf v_{\mathrm{rel}}\) onto the axes \(S\), \(T\), and \(W\), we obtain

\[ \left. \begin{aligned} v_{\mathrm{rel}\,S} &= v_r,\\ v_{\mathrm{rel}\,T} &= v_n - \Omega r\cos\psi\sin a,\\ v_{\mathrm{rel}\,W} &= \Omega r\cos\psi\cos a, \end{aligned} \right\} \tag{11} \]

where \(a\) is the azimuth of the direction of the absolute velocity \(v\), and \(v_r\) and \(v_n\) are the radial and transverse components of the absolute 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{12} \]

Fig. 2.

On the basis of formulas (9) and (11), using the relations (Fig. 1)

\[ \cos\psi\sin a=\cos i, \]

\[ \cos\psi\cos a=\sin i\cos u, \]

for the projections of the perturbing acceleration of the resistance force we obtain the expressions:

\[ \left. \begin{aligned} S_2 &= -\frac{c_x s}{2m}\rho v_{\mathrm{rel}}v_r,\\ T_2 &= -\frac{c_x s}{2m}\rho v_{\mathrm{rel}}(v_n-\Omega r\cos i),\\ W_2 &= \frac{c_x s}{2m}\rho v_{\mathrm{rel}}\Omega r\sin i\cos u, \end{aligned} \right\} \tag{13} \]

where

\[ v_{\mathrm{rel}}= \sqrt{ v_r^2+(v_n-\Omega r\cos i)^2+\Omega^2 r^2\sin^2 i\cos^2 u }, \]

and \(v_r\) and \(v_n\) are determined by formulas (12).

As the curve of the distribution of atmospheric density with height, we shall adopt the dependence obtained on the basis of the data cited in \(^{2}\). We approximate this dependence by the formulas

\[ \left. \begin{aligned} \rho &= \rho_0\Delta,\\ \Delta &= \frac{x}{\left(1+\dfrac{y-y_0}{\xi}\right)^k};\qquad y=\frac{p}{1+e\cos\vartheta}-R, \end{aligned} \right\} \tag{14} \]

where \(\rho\) is the air density for some fixed height; the constants \(x\), \(\xi\), \(y_0\), \(k\) take the following values:

\[ \begin{aligned} &\text{for } 100\ \text{km}\le y\le 150\ \text{km}\quad x=1, \qquad \xi=55,\quad y_0=100,\quad k=8,\\ &\text{for } 150\ \text{km}\le y\le 250\ \text{km}\quad x=0.005667,\qquad \xi=100,\quad y_0=150,\quad k=7,\\ &\text{for } 250\ \text{km}\le y\quad x=0.00004428,\quad \xi=215,\quad y_0=250,\quad k=6. \end{aligned} \]

The system of equations (3), (4) is a system of nonlinear differential equations with respect to the osculating elements of the orbit \(p\), \(e\), \(\omega\), \(\Omega\), \(i\) and the time \(t\). For this system there is no analytic solution.

In direct integration of this system by numerical methods, a large error may accumulate in the determination of the initial parameters, since over the entire lifetime an artificial satellite may make thousands of revolutions around the Earth, while the step of integration in \(\vartheta\) is limited. Moreover, the time required in this case for solving the system (3), (4) on a machine proves to be unjustifiably large.

In order to avoid the indicated shortcomings, we transform the system of equations (3), (4) as follows. Introduce the functions \(z_1,\alpha_2,\ldots,\alpha_6\):

\[ \begin{gathered} \alpha_1=\int_{\vartheta_0}^{\vartheta}\frac{dp}{d\vartheta}\,d\vartheta;\quad \alpha_2=\int_{\vartheta_0}^{\vartheta}\frac{de}{d\vartheta}\,d\vartheta;\quad \alpha_3=\int_{\vartheta_0}^{\vartheta}\frac{d\omega}{d\vartheta}\,d\vartheta;\\ \alpha_4=\int_{\vartheta_0}^{\vartheta}\frac{d\Omega}{d\vartheta}\,d\vartheta;\quad \alpha_5=\int_{\vartheta_0}^{\vartheta}\frac{di}{d\vartheta}\,d\vartheta;\quad \alpha_6=\int_{\vartheta_0}^{\vartheta}\frac{dt}{d\vartheta}\,d\vartheta, \end{gathered} \tag{15} \]

where the derivatives \(\dfrac{dp}{d\vartheta}, \dfrac{de}{d\vartheta},\ldots\) are determined by equations (3), (4). Obviously, when \(\vartheta=\vartheta_0+2\pi\), the functions \(\alpha_1,\alpha_2,\ldots,\alpha_5\) are equal, respectively, to the values of the perturbations of the osculating elements \(p,e,\ldots,i\) over one revolution, while \(\alpha_6\) is equal to the corresponding time of motion.

Let us note that, since the true anomaly is reckoned from the direction to perigee, and the direction to perigee itself changes, in absolute space, when \(\vartheta\) changes from \(\vartheta_0\) to \(\vartheta_0+2\pi\), the satellite will not make exactly one complete revolution around the Earth. Therefore the expression given above is not entirely rigorous. However, in what follows, for convenience, we shall use this term, which plays a purely auxiliary role. (The time of motion of an artificial satellite along the orbit is determined quite accurately with the aid of the corresponding differential equation.)

The perturbations of the osculating elements of the satellite’s orbit over one revolution, for motion at altitudes where the influence of atmospheric resistance is smaller than or comparable with the influence of the noncentrality of the Earth’s gravitational field, are very small. Therefore it may be assumed that throughout the entire motion (except, perhaps, for the last several tens of revolutions) the indicated perturbations are, with great accuracy, equal to the derivatives of the orbital elements with respect to the number of revolutions \(N\). Then, on the basis of the equalities (15), we shall have:

\[ \begin{gathered} \frac{dp}{dN}=\alpha_{1k},\quad \frac{de}{dN}=\alpha_{2k},\quad \frac{d\omega}{dN}=\alpha_{3k},\\ \frac{d\Omega}{dN}=\alpha_{4k},\quad \frac{di}{dN}=\alpha_{5k},\quad \frac{dt}{dN}=\alpha_{6k}, \end{gathered} \tag{16} \]

where

\[ \alpha_{ik}=[\alpha_i]_{\vartheta=\vartheta_0+2\pi}\quad (i=1,2,\ldots,6). \]

The obtained system of differential equations (16), whose right-hand sides are expressed in terms of the values of the combined system of integrals (15), where \(\vartheta=\vartheta_0+2\pi\), describes the change of the satellite’s orbit with time. Let us note that the solution of the system of equations (16) is a discrete sequence of values of the osculating orbital elements for integer values of \(N\), or for values of the true anomaly \(\vartheta=\vartheta_0+2\pi k\) \((k=1,2,\ldots)\).

The transition from the solution of the system of differential equations (3), (4) to the solution of the system of differential equations (16) is indicated geometrically in Fig. 3. In Fig. 3 the solid line shows the curve, show-

producing periodic measurements of the element \(p\), described by the first equation of system (4), over the first four revolutions. The dotted line shows the curve of variation of the element \(p\) according to the first equation (16).

The values of the integrals in the right-hand sides of equations (16) should be computed by integrating the differential equations (3), (4) with respect to \(\vartheta\) over the interval from \(\vartheta_0\) to \(\vartheta_0+2\pi\).

Fig. 3.

In this case, the solution of the system of differential equations (16) reduces to double integration: the outer one with respect to the argument \(N\), and the inner one with respect to the argument \(\vartheta\); moreover, the latter must be performed at each point of the outer integration when computing the right-hand sides of the equations.

For integrating the system of equations (16) with respect to the number of revolutions \(N\), any method of numerical integration of ordinary differential equations with a variable integration step may be used. At high altitudes this step may be considerable, of the order of tens and hundreds of revolutions (when integrating by the Runge—Kutta method), since here the perturbations of the elements over one revolution will be very small and will differ little from one another in passing from revolution to revolution. In passing to low altitudes, the step must be decreased.

From the fact that, when integrating the system of equations (16) by the Runge—Kutta method, the right-hand sides must be computed at the points \(N=N_0+\dfrac{\Delta N}{2}\), where \(\Delta N\) is the integration step, and \(N\) must be an integer, it follows that the step \(\Delta N\) must be an even number. Therefore, when solving the system of equations (16) on an electronic machine with automatic variation of the step depending on the methodological error of the quantities sought (decreasing or increasing it by a factor of 2), it is advisable to choose the initial step \(\Delta N\) equal to \(2^n\).

We note that, as the original differential equations describing the motion of an artificial satellite, one could use the equations for the osculating elements of the orbit with respect to the parameter \(u\), and integrate this system from \(u_0\) to \(u_0+2\pi\) at each corresponding point of the outer system. Since the parameter \(u\) is measured from the equatorial plane, which is fixed in absolute space, and the inclination of the orbit changes very insignificantly with time, when the parameter \(u\) changes by \(2\pi k\), where \(k\) is an integer, the satellite will be located approximately at one and the same latitude.

Therefore the argument of latitude \(u\), as the independent variable, may prove more convenient than the true anomaly \(\vartheta\).

§ 2. CALCULATION OF THE ORBIT OF AN ARTIFICIAL SATELLITE

By the method presented above, an example of the orbit of an artificial satellite was calculated on the high-speed electronic computer of the Academy of Sciences of the USSR. The orbit was calculated for a spherical satellite weighing 10 kg and having a diameter of 0.5 m. The aerodynamic drag coefficient \(C_x\) was taken equal to 2. The initial values of the orbital elements were taken to be:

\[ h_a = 1285\ \text{km},\quad h_{\pi_0}=320\ \text{km},\quad i_0=45^\circ,\quad \omega_0=90^\circ,\quad \Omega_0=129^\circ \]

(\(h_a\) and \(h_{\pi_0}\) are the initial values of the heights of apogee and perigee, respectively).

The results of the calculations are given in Fig. 4, where the curves of variation of the parameter \(p\) (in km), eccentricity \(e\), angular distance of perigee from the node \(\omega\) (in degrees), and longitude of the ascending node \(\Omega\) (in degrees) are shown over a time interval equal to 700 days.

Fig. 4.

A characteristic feature of the curves for the parameter \(p\) and eccentricity \(e\) is that they are oscillatory in character. The period of the oscillations is \(\sim 36\) days and coincides with the period during which—

the osculating element \(\omega\), which determines the position of the perigee of the osculating ellipse with respect to the ascending node, changes by \(\pi\). In contrast to the short-period oscillations of the osculating elements \(p\) and \(e\), which occur when the true anomaly \(\vartheta\) changes over the interval from \(\vartheta_0+2\pi(k-1)\) to \(\vartheta_0+\pi k\), we shall call the oscillations indicated above long-period oscillations.

The presence of long-period oscillations is explained by the fact that, when the true anomaly \(\vartheta\) changes by \(2\pi\), both the osculating ellipse of the orbit and its orientation in absolute space change; at the same time the osculating element \(\omega\), which determines the position of the perigee of the osculating ellipse with respect to the equatorial plane, changes. Therefore, when \(\vartheta\) changes by \(2\pi k\) \((k=1,2\ldots)\), the satellite will occupy different positions relative to the equatorial plane, and the action of the force due to the deviation of the Earth’s gravitational field from a central one, which depends on latitude, will be different.

If only the perturbing force due to the noncentrality of the Earth’s gravitational field is present, the corresponding curves of variation of the osculating elements \(p\) and \(e\) at the points \(\vartheta=\vartheta_0+2\pi k\) must have the form of oscillatory curves with a horizontal mean line. This follows from the fact that the force due to the noncentrality of the Earth’s gravitational field is conservative. The perturbing action of atmospheric resistance leads to dissipation of the satellite’s energy and causes a monotonic change of the osculating elements (at the points \(\vartheta=\vartheta_0+2\pi k\)). The combined action of the two principal perturbing forces indicated leads to the result that the curves of variation of the osculating elements of the orbit of the artificial satellite at points equidistant from the perigees of the osculating ellipses have the form shown in Fig. 4.

The curves given in Fig. 4 make it possible to judge what the secular perturbations of the osculating elements of the orbit of an artificial satellite are over a certain interval of time. Indeed, it is not difficult to see that over a period of time equal to 700 days the secular perturbations of the orbital elements are:

\[ \Delta p=-414\ \text{km},\quad \Delta e=-0.0564,\quad \Delta\omega=-3860^\circ,\quad \Delta\Omega=-3529^\circ . \]

It should be noted that, in the motion of an artificial satellite in the noncentral field of the Earth’s gravitation, for the orbit under consideration the perigee of the osculating ellipse changes its position with respect to the plane of the equator over time. Over a period of 700 days, the perigee of the osculating ellipses will make about 11 revolutions around the Earth.

The ascending node for the orbit under consideration recedes in the direction opposite to the motion of the Earth, with a rate of \(\sim 5^\circ\) per day.

In conclusion, we note that in computing the indicated orbit of the satellite on a machine by the method described above, over the time of its motion during 700 days, i.e., of the order of two years, about four hours of machine time were expended.

References

1. G. N. Duboshin, Introduction to Celestial Mechanics, Moscow–Leningrad, ONTI, 1938.
2. S. K. Mitra, The Upper Atmosphere, Moscow, IL, 1955.
3. N. I. Idelson, Potential Theory with Applications to the Theory of the Figure of the Earth and Geophysics, Moscow, ONTI, 1936.
4. D. E. Okhotsimskii, T. M. Eneev, G. P. Taratynova — in this issue of Uspekhi Fizicheskikh Nauk.