Effect of atmospheric drag on artificial satellite orbits

An artificial satellite that is in a low-altitude orbit around the Earth interacts with the upper atmosphere. Although the air in the upper atmosphere has a very low mass density, when integrated over a sufficiently long period of time, this interaction can significantly modify the satellite orbit, causing it to decay (i.e., to decrease in altitude). Let us investigate this effect. We shall work in the geocentric frame of reference introduced in the previous section.

The force per unit mass that the atmosphere exerts on an orbiting satellite can be modeled as (Batchelor 2000)

$${\bf F} = -\frac{1}{2}\,\frac{\rho\,C_D\,A}{m}\,\vert\dot{\bf r}\vert\,\dot{\bf r},$$ (10.130)

where $\rho$ is the mass density of the atmosphere (at the satellite's position), $C_D$ the drag coefficient, $A$ the cross-sectional area of the satellite perpendicular to its direction of motion, $m$ the satellite mass, and $\dot{\bf r}$ the satellite velocity. It can be seen that the force is proportional to the square of the satellite speed, and is oppositely directed to its instantaneous direction of motion (i.e., the force constitutes a drag). The drag coefficient is an order unity dimensionless constant that depends primarily on the satellite's shape. For the case of a spherical satellite, $C_D\simeq 2.2$ (Cook 1965).

The density distribution of the terrestrial atmosphere is conveniently modeled as (de Pater and Lissauer 2010)

$$\rho(r) = \rho_0\exp\left[-\frac{(r-R)}{H}\right],$$ (10.131)

where $r$ measures radial distance from the Earth's center, $R=6.4\times 10^6\,{\rm m}$ is the terrestrial radius, $\rho_0\simeq 1.3\,{\rm kg}\,{\rm m}^{-3}$ the atmospheric density at ground level (i.e., $r=R$), and $H\simeq 8.5\times 10^3\,{\rm m}$ the atmospheric scale height. Obviously, the previous formula is only valid when $r\geq R$.

Assuming that the atmospheric drag force, (10.130), is small compared to the force of gravitational attraction between the Earth and the satellite—and can, thus, be treated as a perturbation—the satellite's orbit can be modeled as Keplerian ellipse whose six elements evolve slowly in time under the influence of the drag. Let $r$,  $\theta$,  $z$ be cylindrical coordinates in a reference frame that is aligned with the satellite's instantaneous orbital plane, as described in Section I.1. Here, $\theta$ is the satellite's true anomaly. (See Section 4.11.) Making use of the analysis of Appendix I, we can write

$$\dot{\bf r} = \skew{5}\dot{r}\,{\bf e}_r + r\,\skew{5}\dot{\theta}\,{\bf e}_\theta,$$ (10.132)

$${\bf F} = F_r\,{\bf e}_r+ F_\theta\,{\bf e}_\theta+ F_z\,{\bf e}_z,$$ (10.133)

$$r =a\,(1-e\,\cos E),$$ (10.134)

$$\skew{5}\dot{r} = \frac{h\,e}{a\,(1-e^{\,2})^{1/2}}\,\frac{\sin E}{1-e\,\cos E},$$ (10.135)

$$h \equiv r^{\,2}\,\skew{5}\dot{\theta} = (1-e^{\,2})^{1/2}\,(\mu\,a)^{1/2}.$$ (10.136)

Here, $a$, $e$, $h$, and $E$ are the satellite's orbital major radius, eccentricity, angular momentum per unit mass, and eccentric anomaly, respectively. (See Section 4.11.) Moreover, $\mu= G\,(M+m)$, where $M$ is the terrestrial mass. It follows from the previous equations that

$$F_r = -\frac{1}{2}\,\frac{C_D\,A\,\rho(a)}{m}\,\frac{\mu}{a}\,{\rm e}^{\,\alpha\,\cos E}\,\frac{(1+e\,\cos E)^{1/2}}{(1-e\,\cos E)^{3/2}}\,e\,\sin E,$$ (10.137)

$$F_\theta = -\frac{1}{2}\,\frac{C_D\,A\,\rho(a)}{m}\,\frac{\mu}{a}\,\,{\rm ... ...a\,\cos E}\,\frac{(1+e\,\cos E)^{1/2}}{(1-e\,\cos E)^{3/2}}\,(1-e^{\,2})^{1/2},$$ (10.138)

$$F_z = 0,$$ (10.139)

where $\alpha = a\,e/H$.

As has already been mentioned, the satellite's instantaneous orbit is a Keplerian ellipse characterized by six orbital elements, which we choose to be the major radius, $a$, the mean anomaly at epoch, ${\cal M}_0$, the eccentricity, $e$, the argument of the perigee, $\omega$, the inclination (to the Earth's equatorial plane), $I$, and the longitude of the ascending node (measured with respect to the vernal equinox), ${\mit \Omega }$. (See Section 4.12.) These elements evolve slowly in time, under the influence of the atmospheric drag force. For the case in hand, this time evolution is most conveniently specified in terms of the Gauss planetary equations, which take the form (see Appendix I)

$$\frac{\skew{3}\dot{a}}{a} = \frac{2\,h}{\mu\,(1-e^{\,2})}\left[e\,\sin\theta\,F_r + (1+e\,\cos\theta)\,F_\theta\right],$$ (10.140)

$$\skew{5}\dot{\cal M}_0 = \frac{h}{\mu}\,\frac{(1-e^{\,2})^{1/2}}{e}\left(\left[\cos\thet... ... -\left[1+\frac{1}{(1-e^{\,2})}\,\frac{r}{a}\right]\sin\theta\,F_\theta\right),$$ (10.141)

$$\skew{3}\dot{e} = \frac{h}{\mu}\left[\sin\theta\,F_r+(\cos\theta+\cos E)\,F_\theta \right],$$ (10.142)

$$\skew{3}\dot{\omega} = -\frac{h}{\mu}\,\frac{1}{e}\left[\cos\theta\,F_r-\left(\frac{2+... ...theta\,F_\theta \right] -\frac{\cos I\,\sin(\omega+\theta)\,r\,F_z}{h\,\sin I},$$ (10.143)

$$\dot{I} = \frac{\cos(\omega+\theta)\,r\,F_z}{h},$$ (10.144)

$$\skew{5}\dot{\mit\Omega} =\frac{\sin(\omega+\theta)\,r\,F_z}{h\,\sin I}.$$ (10.145)

Given that (see Section 4.11 and Exercise 16)

$$\cos \theta = \frac{\cos E-e}{1-e\,\cos E},$$ (10.146)

$$\sin \theta = \frac{(1-e^{\,2})^{1/2}\,\sin E}{1- e\,\cos E},$$ (10.147)

the Gauss planetary equations can be combined with Equations (10.137)–(10.139) to give

$$\frac{\skew{3}\dot{a}}{a} = -\frac{C_D\,A\,a\,\rho(a)}{m}\,n\,{\rm e}^{\,\alpha\,\cos E}\,\frac{(1+e\,\cos E)^{3/2}}{(1-e\, \cos E)^{3/2}},$$ (10.148)

$$\skew{5}\dot{\cal M}_0 = \frac{C_D\,A\,a\,\rho(a)}{m}\,n\,e^{\,-1}\,{\rm e}^{\,\alpha\,\... ...\frac{(1+e\,\cos E)^{1/2}}{(1-e\, \cos E)^{3/2}} \,(1-e^{\,3}\,\cos E)\,\sin E,$$ (10.149)

$$\skew{3}\dot{e} =-\frac{C_D\,A\,a\,\rho(a)}{m}\,n\,(1-e^{\,2})\,{\rm e}^{\,\alpha\,\cos E}\,\frac{(1+e\,\cos E)^{1/2}}{(1-e\, \cos E)^{3/2}}\,\cos E,$$ (10.150)

$$\skew{3}\dot{\omega} = -\frac{C_D\,A\,a\,\rho(a)}{m}\,n\,\frac{(1-e^{\,2})^{1/2}}{e}\,... ...}^{\,\alpha\,\cos E}\,\frac{(1+e\,\cos E)^{1/2}}{(1-e\, \cos E)^{3/2}}\,\sin E,$$ (10.151)

$$\dot{I} = 0,$$ (10.152)

$$\skew{5}\dot{\mit\Omega} = 0.$$ (10.153)

Here, $n=(\mu/a^{\,3})^{1/2}$ is the unperturbed mean orbital angular velocity. It is immediately apparent that the atmospheric drag force does not give rise to any change in the inclination of the satellite's orbital plane (because $I$ and ${\mit \Omega }$ are both constant). This is not surprising, because the drag force has no component normal to this plane (i.e., $F_z=0$).

As explained previously, we are working in the limit in which the atmospheric drag force is perturbative. This limit requires the relative changes in the satellite's orbital elements induced by the drag in an orbital period, $T=2\pi/n$, to all be small. It is apparent from Equations (10.148)–(10.151) that this is the case provided $2\pi\,a\,A\,\rho(a)\ll m$; in other words, as long as the mass of air encountered by the satellite in a single orbit, which is $2\pi\,a\,A\,\rho(a)$, is much less than the satellite's mass. Assuming that we are in the perturbative limit (i.e., $2\pi\,a\,A\,\rho(a)\ll m$), the evolution of the satellite's orbital elements takes place on a timescale that is much longer than its orbital period. We can concentrate on this evolution, and filter out any relatively short-term oscillations in the elements, by averaging expressions (10.148)–(10.151) over an orbital period. A suitable orbit-average operator is

$$\langle\cdots\rangle \equiv \frac{1}{T}\int_0^T (\cdots)\,dt = \oint (\cdots)\,(1-e\,\cos E)\,\frac{dE}{2\pi}.$$ (10.154)

Here, we have made use of the fact that $E-e\,\sin E = {\cal M}$, and $d{\cal M}/dt = n$, for a Keplerian orbit. (See Section 4.11.) Thus, we deduce that

$$\left\langle\frac{\skew{3}\dot{a}}{a}\right\rangle = -\frac{C_D\,A\,a\,\rho(a)}{m}\,n\oint {\rm e}^{\,\alpha\,\cos E}\,\frac{(1+e\,\cos E)^{3/2}}{(1-e\, \cos E)^{1/2}}\,\frac{dE}{2\pi},$$ (10.155)

$$\langle \skew{3}\dot{e}\rangle =-\frac{C_D\,A\,a\,\rho(a)}{m}\,n\,(1-e^{\,2})\oint{\rm e}^{\,\al... ... E}\,\frac{(1+e\,\cos E)^{1/2}}{(1-e\, \cos E)^{1/2}}\,\cos E\,\frac{dE}{2\pi},$$ (10.156)

$$\langle\skew{5}\dot{\cal M}_0 \rangle = 0,$$ (10.157)

$$\langle \skew{3}\dot{\omega} \rangle = 0.$$ (10.158)

It follows that the atmospheric drag force causes the major radius and eccentricity of the satellite's orbit to both decay monotonically in time [because the right-hand sides of Equations (10.155) and (10.156) are both negative]. On the other hand, the force does not produce any precession in the perigee of the orbit (because $\langle \skew{3}\dot{\omega} \rangle =0$). Furthermore, the fact that $\langle\skew{5}\dot{\cal M}_0 \rangle=0$ implies that the drag force does not modify the Keplerian result that the mean orbital angular velocity is $(\mu/a^{\,3})^{1/2}$. Now, it is easily demonstrated (see Exercise 4) that the satellite's orbit-averaged kinetic energy is

$$\langle K\rangle =\frac{\mu}{2\,a}.$$ (10.159)

Hence, we deduce that, despite the negative work done on the satellite by the atmospheric drag force, its kinetic energy actually increases as its orbit decays (i.e., as $a$ decreases). Evidently, the negative work done on the satellite by the drag force is more than offset by the positive work done by gravity as its altitude decreases.