# Analytical and Semi-Analytical Treatment of the Satellite Motion in a Resisting Medium

Authors: S. E. Abd El-Bar; F. A. Abd El-Salam
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

The orbital dynamics of an artificial satellite in the Earth's atmosphere is considered. An analytic first-order atmospheric drag theory is developed using Lagrange's planetary equations. The short periodic perturbations due to the geopotential of all orbital elements are evaluated. And to construct a second-order analytical theory, the equations of motion become very complicated to be integrated analytically; thus we are forced to integrate them numerically using the method of Runge-Kutta of fourth order. The validity of the theory is checked on the already decayed Indian satellite ROHINI where its data are available.


---

## Body

## 1. Introduction

The problem of the orbital motion of celestial bodies in a resisting medium can be traced back to Newton's Principle [1] and Laplace's analysis of a possibly finite velocity of gravitation [2]. However, with the advent of the space age there arose the need of precisely predicting the orbital behavior of the newly launched bodies.When the artificial satellites move in Low Earth Orbits (in brief LEO) their orbital motion is perturbed by the resisting force. Therefore, an adequate analytical description of the drag force is required to predict precisely the orbital perturbations due to atmospheric drag. These perturbations due to the atmospheric drag are more complicated and difficult to deal with than orbits perturbed by gravity because of the nonconservative property of the drag force. The adequate allowance for atmospheric drag turned out to be a problem for unforeseen complexity [3–8]. In order to be able to solve the equations of motion including the drag force analytically, one is forced to do many approximations. On the other hand, from the computational point of view, having an analytical solution at hand, we can jump from the initial conditions to the new state, which makes the computation extremely fast.The order of magnitude of the atmospheric drag acting on a satellite depends on the altitude of the satellite. The orbital behavior of an LEO satellite could be strongly influenced by atmospheric drag [4, 9, 10], and therefore affects the quality of the remote sensing of the satellite [11, 12]. Traditionally, the problem of orbit disturbance by atmospheric drag is solved by numerical integration [13–18]. An analytical solution gives the theoretical integrals and shows the physical effects with very clear spectral properties [19–24]. This may give a direct insight into the physical phenomenon of the disturbance [25–27]. Bezdeka and Vokrouhlicky [28] developed a semianalytic theory of motion for close-Earth spherical satellites including drag and gravitational perturbations. Xu et al. [29] derived an analytical solution of a satellite orbit disturbed by atmospheric drag. They first transformed the disturbance force vector and rotated it to the orbital frame so that it can be used in the simplified Gaussian equations of satellite motion. The disturbances are separated into three parts: shortperiodic terms with triangular functions of the mean anomalyM, longperiodic terms with triangular functions of the argument of perigee and inclination (ω,I), and secular terms [nonperiodic functions of semimajor axis and eccentricity (a,e)].The main effects of the force are to bring out large secular decreases in the elementsa, the semimajor axis, and e, the eccentricity, that cause the satellite plunges toward Earth. This perturbing influence is one of the most important perturbing forces in the altitude regime from 150 to 600 km which is called the drag force. The drag force,F→D, per unit mass of the satellite can be represented by(1.1)F→D=-12CD(Am)ρ‖V→‖V→,
where CD is a dimensionless drag coefficient, empirically determined to have a value close to 2; it can be taken as 2.2 with an error (standard deviation) which should not exceed 5% unless the assumption about the accommodation coefficient is greatly in error. But in fact the drag coefficient depends on the mean free path of the atmospheric molecules with respect to the linear dimensions of the satellite; therefore, the orbital theory should be applied with caution for large satellites in their last revolutions [4]. After that Sehnal [30] derived a theoretical expression for the dimensionless drag coefficient (see (3.1) there). V→ is the velocity of the satellite relative to the atmosphere (called the ambient velocity), A represents the effective cross-sectional area of the satellite, which is to be found by averaging all possible projected areas of the satellite onto a plane perpendicular to V→, m is the mass of the satellite, and ρ is the density of the atmosphere. The density ρ depends on position and time in a very complicated manner. Its variation with height depends on the form of the atmosphere, since it decreases rapidly with slight growing in the altitude; the effect of the flattening of the atmosphere is rather significant.


## 2. Rotation of the Atmosphere

This effect is so small, so a simple approximation is justified. The differential rotation of the atmosphere complicates the numerical calculation and it gives no sense in this very small effect. So we will assume that the atmosphere rotates with the same angular velocity of the Earthσ→, then(2.1)V→=v→-σ→×r,→
where v→ is the orbital velocity of the satellite, and r→ its radius vector. Thus(2.2)V2=v2-2v→⋅(σ→×r→)+(σ→×r→)2.If(î,ĵ,k̂) are the base vectors of the rectangular geocentric equatorial frame associated with the spherical coordinates (r,λ,δ), then σ→=σk̂ with σ=7.27205217×10-5rad/s, and(2.3)V2=v2-2σhcos⁡I+σ2r2cos⁡2δ,
where I is the orbital inclination,δ is the satellite's declination, and h is the integral of areas. Now the ambient velocity can be written in the form:(2.4)V2=v2kR,where(2.5)kR=1-2σhcos⁡Iv2+σ2r2cos⁡2δv2.The third term in kR may usually be neglected, being less than about 1/250, and the second term which is of order of 1/15 may be evaluated at perigee, where the density is much greater than elsewhere on the orbit; thus we may effectively regard kR as a constant.


## 3. The Atmospheric Drag Force Model

We have assumed that the atmospheric density is approximated as an exponential function which is given by(3.1)ρ=ρpoexp⁡[β(rpo-r)],
where ρpo is the atmospheric density at the initial perigee point rpo, β is taken as a constant, 1/β is about 50 kilometers for most satellite orbits, and r is given by relation: (3.2)r=a(1-ecos⁡E),
where a, e are the orbital semimajor axis and eccentricity and E is the eccentric anomaly of the satellite. Now(3.3)ρ=ρpoexp⁡[β(ao-a-xo)+βxcos⁡E],wherex=ae,xo=aoeo,V→ is the velocity of satellite relative to ambient atmosphere and is generally through(3.4)V→=v→-ω→a×r,→
where v→ is the velocity of satellite with respect to the Earth's center, ω→a is the west-to-east angular velocity of the atmosphere and r→is the position vector of the satellite.The radialVs, transverse Vt, and orthogonal Vw components of the velocity are given by(3.5)Vs=μa(1-e2)esin⁡f,Vt=μa(1-e2)(1+ecos⁡f)-rωacos⁡I,Vw=ωarsin⁡Icos⁡(f+ω),whereμ=G(M⊕+m)≈GM⊕, the Earth's gravitational parameter, is the product of the universal constant of gravity G and the total mass of the Earth M⊕, the mass of the satellite m, and f,ω are the true anomaly and the argument of the perigee of the satellite. Then, the radial SD, transverse TD, and orthogonal WD components of the drag force are given by(3.6)SD=-12CDAmρ‖V→‖μa(1-e2)esin⁡f,TD=-12CDAmρ‖V→‖[μa(1-e2)(1+ecos⁡f)-rωacos⁡I],WD=-12CDAmρ‖V→‖rωasin⁡Icos⁡(f+ω),where(3.7)‖V→‖=[μa(1-e2)μa(1-e2)(e2+2ecos⁡f+1)+r2ωa2(cos⁡2⁡I+sin2Ico⁡s2(f+ω))-2rωacos⁡Iμa(1-e2)(1+ecos⁡f)]1/2.Using the relations(3.8)cos⁡f=cos⁡E-e1-ecos⁡E,(3.9)sin⁡f=1-e2sin⁡E1-ecos⁡E.Equations (3.6), can be written in the form(3.10)SD=-bn1-cos⁡Eaeρ‖V→‖sin⁡E,TD=-bn(1-e2)1-cos⁡Eaρ‖V→‖[1-d(1-edcos⁡E)21-e2],WD=-bωaaρ‖V→‖(1-ecos⁡E)sin⁡Icos⁡(f+ω),where(3.11)‖V→‖=μa1+ecos⁡E1-ecos⁡E{1-d1+ecos⁡E1-ecos⁡E+ωa22n2[1+ecos⁡E1-ecos⁡E]2+[(1-e2cos⁡2⁡E)sin2Icos⁡2(f+ω)+e2cos⁡2Isin2I]1-d1+ecos⁡E1-ecos⁡E+ωa22n2[1+ecos⁡E1-ecos⁡E]2[1+ecos⁡E1-ecos⁡E]}+⋯,d=ωan1-e2cos⁡I,wheren=μ/a3 is the mean motion.For a close Earth’s satellite the ratioωa/n is about 1/15; the third term in the expansion is thus only about 1/500 of the leading term and will be ignored. For satellites that are not close, the ratio will be larger but the atmospheric density that will multiply this expression will then be so small that the product will still be unimportant. we thus take(3.12)‖V→‖=μa1+ecos⁡E1-ecos⁡E[1-d1+ecos⁡E1-ecos⁡E],b=12CD(Am).


## 4. Lagrange's Planetary Equations

The orbit is uniquely determined by 6 elements, namely, the orbital elements. The Lagrange planetary equations are a set of 6 differential equations which measure the rate change of these orbital elementsa,e,i,ω,Ω and any time element, for example, the mean anomaly M(4.1)ȧ=2esin⁡fn1-e2S+2a1-e2nrT,ė=1-e2sin⁡fnaS+1-e2nae[a(1-e2)r-ra]T,I⋅=rcos⁡(f+ω)na21-e2W,Ω̇=rsin⁡(f+ω)na21-e2sin⁡IW,ω̇=-Ω̇cos⁡I-1-e2nae[Scos⁡f-T(1+ra(1-e2))sin⁡f],Ṁ=n-2na2Sr-ω̇1-e2-Ω̇1-e2cos⁡I.There are many forms of the Lagrange planetary equations. In many cases, some of orbital elements are undefined, and hence, we have to replace this set by a new set of elements valid to define the orbit, and the planetary equations will take a new form in terms of the new element. For the equatorial orbitsI=0 or I=180ο or even when I, e have small values we have near equatorial near circular orbits and the periapsis is undefined and hence the time of pericentric passage,T and ω are undefined or at least meaningless and should be replaced by other elements as the longitude of periapsis, ϖ, defined as ϖ=ω+Ω.


## 5. The Earth's Gravity Force Model

The gravitational field of the Earth is quite irregular. However, it can be represented by an infinite series of spherical harmonics, choosing the principal axes with origin of coordinates at centre of mass in addition assuming also that the gravitational field of the Earth is axially symmetric; then the form adopted is, [31](5.1)U⊕=-μr[1-∑n=2∞(R⊕r)nJnPn(sin⁡δ)],
where R⊕ is the equatorial radius of the Earth, μ is the Earth’s gravitational parameter, (r,δ) are the appeared geocentric coordinates of the satellite, Jnare zonal harmonic coefficients, and Pn(sinδ) the Legendre Polynomials are given by(5.2)Pn(sin⁡δ)=12n∑j=0N(-1)j(2n-2)!(sin⁡δ)n-2jj!(n-j)!(n-2j)!,(n=0,1,2,3,…),
where
(5.3)N={n2neven,n-12nodd.After some lengthy straightforward algebraic calculation we can evaluate the force components in the radialSG, transverse TG, and normal WG components as (5.4)SG=-34J2μR⊕2r4[2-3sin2I(1-cos⁡2(f+ω))]+J3R⊕3μ2r5sin⁡I[(15sin2I-12)sin⁡(f+ω)-5sin2Isin⁡3(f+ω)]+J45R⊕4μ64r6[24-120sin2I+105sin4I+20(6sin2I-7sin4I)×cos⁡2(f+ω)+35sin4Icos⁡4(f+ω)],TG=-3J2μR⊕2r4sin⁡Icos⁡Isin⁡F1,1cos⁡F1,1-J33R⊕3μ4r5sin⁡Icos⁡F1,1[-2+5sin2I-sin2Icos⁡F2,2]-J45R⊕4μ8r6sin2Icos⁡F1,1[(-12+21sin2I)sin⁡F1,1-7sin2Isin⁡F3,3],WG=-3J2μR⊕2r4sin⁡Icos⁡Isin⁡(f+ω)-J33R⊕3μ4r5cos⁡I[-2+5sin2I-sin2Icos⁡2(f+ω)]-J45R⊕4μ8r6sin⁡Icos⁡I[(-12+21sin2I)sin⁡(f+ω)-7sin2Isin⁡3(f+ω)].


## 6. The Analytical Solution

In this section, we proceed to calculate the first-order approximation of the orbital elements caused by the drag and the gravity forces during revolution by integrating Lagrange's planetary equations (4.1) from E=-π to E=π, and retain only terms to order O(e2), after substitution of {S=SD+SG,T=TD+TG,W=WD+WG} into the Lagrange planetary equations we have the perturbations in the elements as(6.1)Δσ=ΔσD+ΔσG,
where σ≡(a,e,I,ω,Ω,M) denotes any of the orbital elements, ΔσD the perturbations due to the drag force, and ΔσG the perturbations due to the gravity force is given by(6.2)ΔaD=-4πa2Ho{[(1-2ωancos⁡I)+32(1+ωancos⁡I)e2]Io+2eI1+12(32+ωancos⁡I)I2e2},ΔaG=J2∑i=15∑j=-22αij′′cos⁡Fi,j+J3∑k=17∑l=-33βkl′′sin⁡Fk,l+J4∑m=19∑n=-44γmn′′cos⁡Fm,n,ΔeD=-4πabH0{e2(1+5ωa2ncos⁡I)I0+[(1-2ωancos⁡I)-e2(58-2ωancos⁡I)]I1+e2(1+3ωa2ncos⁡I)I2+e28I3},ΔeG=J2∑i=15∑j=-22δ̅ijcos⁡Fi,j+J3{X1f+∑m=17∑n=-33σ̅mnsin⁡Fm,n}+J4{X2f+∑k=19∑l=-44Q̅klcos⁡Fk,l},ΔID=-πabωaHosin⁡In{[(94-3ωancos⁡I+178cos⁡2ω)e22(1-ωancos⁡I)+e22(94-3ωancos⁡I+178cos⁡2ω)+55ωa8ncos⁡Icos⁡2ωe22(94-3ωancos⁡I+178cos⁡2ω)]Io+e[-1+3ωancos⁡I+12(5ωancos⁡I-9)cos⁡2ω]I1+[e22(1-ωancos⁡I)+e22(5-17ωancos⁡I)+(1-21ωancos⁡I)cos⁡2ω(1-ωancos⁡I)+e22(5-17ωancos⁡I)]I2+e2(ωancos⁡I-1)cos⁡2ωI3+e216(9-11ωancos⁡I)cos⁡2ωI4},ΔIG=J2∑i=13Zicos⁡Fi,2+J3{X3f+∑k=15∑l=-13C̅̅klsin⁡Fk,l}+J4{X4f+∑p=17∑q=-24S̅̅pqcos⁡Fp,q},ΔΩD=-abωaH0ωan{([3e4(1-ωancos⁡I)+e28(41-64ωancos⁡I)]sin⁡2ω+e22(1-ωancos⁡I)cos⁡2ω)Io+[e2(9-7ωancos⁡I)-e22(6-11ωancos⁡I)]sin⁡2ωI1+([(1-ωancos⁡I)+e(1-ωancos⁡I)+e24(9-33ωancos⁡I)sin⁡2ω+e22(1-ωancos⁡I)cos⁡2ω])I2+[e2(3ωancos⁡I-1)-e24(4-9ωancos⁡I)]×sin⁡2ωI3+e8[(2-2ωancos⁡I)+(5-13ωancos⁡I)]sin⁡2ωI4+e24ωancos⁡Isin⁡2ωI5},ΔΩG=J2{X8f+∑b=13∑c=02ρ̅bcsin⁡Fb,c}+J3{X9f+∑α=15∑β=-13q̅αβcos⁡Fα,β}+J4{X10f+∑i=17∑j=-24p̅ijsin⁡Fi,j},ΔωD=ΔΩDcos⁡I,ΔωG=J2{X5f+∑p=15∑q=-22U̅pqsin⁡Fp,q}+J3{X6f+∑k=17∑l=-33V̅klcos⁡Fk,l}+J4{X7f+∑m=19∑n=-44Z̅mnsin⁡Fm,n},ΔM=-2πn,where(6.3)I0(x)=1π∫0πcos⁡nEexp⁡(xcos⁡E)dEis the first kind modified Bessel function of ordern, and by using the values of I0(c) and I1(c) one finds that(6.4)I2=I0-2cI1,I3=-4cI0+(1+8c2)I1,I4=(1+24c2)I0-(1-6c2)I1,I5=-12c(1+16c2)I0+(1+72c2+384c4)I1,whereHo is the scale height. It is given by(6.5)Ho=ρpoexe[β(ao-a-xo)].Finally the nonvanishing coefficients included in the gravity perturbations are previously computed by the authors [32].


## 7. The Numerical Solution

It is well known that the numerical solution methods provide accurate ephemeris of a satellite with respect to a force of a complex mathematical formulation. So in this section we will solve numerically the equations of motion after substituting the drag force componentsSD, TD, WD, given by (3.10) together with the Earth's gravity force components SG, TG, WG, given by (5.4) (we will use Earth's oblate model and retaining terms up to fourth zonal harmonics J4).There are many methods to get the numerical solutions of any first-order differential equations. We will use one of the best methods of the numerical integration schemes, namely, Runge-Kutta fourth-order method.If the Earth was assumed spherical, had no atmosphere, and the satellite is isolated from the other bodies in the solar system, the orbit of a satellite would be an ellipse of constant size and shape in a plane whose direction remains fixed relative to the stars. The orbital elementsa,e,i,ω,Ω would remain unchanged. This simple thought rather insipid situation is upset by the effects of the considered perturbing forces, the Earth's oblateness, and the atmospheric drag.The Earth's oblateness leads to two major gravitational perturbations on the satellite orbits, the first is the rotation of the orbital plane about the Earth's axis in the direction opposite to the satellite's motion so that(ifI<90o) the angle steadily decreases while I remains constant. The second effect is the rotation of the major axis in the orbital plane so that the argument of perigee ω increases. The chief effects of the atmospheric drag on a satellite orbit are fortunately quite different from those of the gravitational field. Since the atmospheric density decreases rapidly as height above the Earth increases, a satellite in an orbit of appreciable eccentricity is affected mostly by drag within a small section of the orbit where it is closest to the Earth. To first-order approximation, the effect of atmospheric drag is to retard the satellite as it passes perigee, with the result that it does not swing out so far from the Earth at the subsequent apogee. The apogee height is reduced while the perigee height remains almost constant. The orbit contracts and becomes more nearly circular; the orbital elements (a,e) decrease steadily [4].Now, a computer program of fourth-order Runge-Kutta method is applied for the solution of this system using the initial values and a fixed step size. Also, we will choose one of the numerical examples to show the decay in two orbital elements, namely,(a,e). Also we will give a numerical example to the analytical solution given by (6.1) in order to compare the results of the numerical and analytical solutions.


## 8. Numerical Example

In this paragraph, we give a numerical example concerning the artificial satellite ROHINI (1980 62A) [33]. The data are as follows.The satellite mass and its cross-sectional areas are, respectively,(m)=35.443kg, (A)=0.319019×10-6. The initial orbital elements (semimajor axis, eccentricity, and the inclination) of the satellite are a0=6989.2057km, e0=0.04367712, I0=44.67198oω0=239.3378okm, Ω0=174.1602o, and M0=25.63974o. And one orbital revolution was elapsed in 1.58835208 hours. The required physical and dynamical constants are(8.1)RÅ=6378.135km,μ=398600.8km3/sec-2,J2=1.08263×10-3,J3=-2.53648×10-6,J4=-1.6233×10-6,ωa=7.27205217×10-6rad/sec,ρ0(300km)=0.0221kgkm-3,ρ0(250km)=0.0681kgkm-3,ρ0(200km)=0.271kgkm-3,CD=2.2.The results of the numerical example of the analytical solution as well as the numerical solution are represented on Figures1(a), 2(a), 3(a), 4(a), 5(a), 6(a), 7(a), 8(a), and 9(a) to show the decay of the semimajor axis, and Figures 1(b)– 9(b) to show the change in the eccentricity. To read the figures properly we used the letters B, C, and D to denote to the numerical example of the analytical solutions while we used the letters E, F, and G to denote for the numerical solutions. The description of these letters is as follows. B: with drag and rotation,C: with drag without rotation,D: without drag,E: with drag and rotation,F: with drag without rotation, andG: without drag.(a) For one hundred revolutions at 300 km height above the Earth's surface. (b) For one hundred revolutions at 300 km height above the Earth's surface.
(a)(b)(a) For five hundred revolutions at 300 km height above the Earth's surface. (b) For five hundred revolutions at 300 km height above the Earth's surface.
(a)(b)(a) For one thousand revolutions at 300 km height above the Earth's surface. (b) For one thousand revolutions at 300 km height above the Earth's surface.
(a)(b)(a) For one hundred revolutions at 250 km height above the Earth's surface. (b) For one hundred revolutions at 250 km height above the Earth's surface.
(a)(b)(a) For five hundred revolutions at 250 km height above the Earth's surface. (b) For five hundred revolutions at 250 km height above the Earth's surface.
(a)(b)(a) For one thousand revolutions at 250 km height above the Earth's surface. (b) For one thousand revolutions at 250 km height above the Earth's surface.
(a)(b)(a) For one hundred revolutions at 200 km height above the Earth's surface. (b) For one hundred revolutions at 200 km height above the Earth's surface.
(a)(b)(a) For five hundred revolutions at 200 km height above the Earth's surface. (b) For five hundred revolutions at 200 km height above the Earth's surface.
(a)(b)(a) For five thousand revolutions at 200 km height above the Earth's surface. (b) For five thousand revolutions at 200 km height above the Earth's surface.
(a)(b)All figures show the decrease in the semimajor axis and in the eccentricity for the rotating and nonrotating atmosphere at different altitudes.


## 9. Analysis of the Solutions

The following dynamical effects are very clear. The decay of the semimajor axis is shown on Figures1(a) to  9(a). We depicted the analytical and numerical solutions on one figure at different altitudes. Similar situation is done for the eccentricity on Figures 1(b)–9(b).(1)
The decrease in the semimajor and in the eccentricity is significant (much bigger) at the high altitudes at the beginning of the motion but it oscillates afterwards, that is, the slope of the rate of change is sometimes high, sometimes moderate, and sometimes low. Therefore, we can conclude that the perturbation on the long range is nonlinear (see Figures10(a) and 10(b)).(2)
The effect of the atmospheric rotation is so small.
(a)
The atmospheric rotation appeared clearly on the figures of numerical solution because we have used the well-known fourth-order Runge-Kutta method.(b)
The atmospheric rotation disappeared on the figures of the analytical solution since the approximation made in the analytical solution is much bigger than this effect.(3)
Disregarding the drag force makes the rate of change of the semimajor axis and the rate of change of the eccentricity very little as seen clearly from the slopes of the figures of the solutions with drag and those without drag. This justifies the fact that at Low Earth Orbits (in brief LEO) the dominant secular perturbation is due to the drag force.(4)
The difference between the solutions with drag and those without drag is sometimes big, sometimes low, and sometimes moderate. This signifies that the perturbation is nonlinear.(a) Long range perturbation in the semimajor axis. (b) Long range perturbation in the eccentricity.
(a)(b)


## 10. Conclusion

In this paper, we studied the effect of the atmospheric drag and the asphericity of the Earth on the orbital dynamics of an artificial satellite in a low Earth orbit. We approximated the density of the atmosphere to an exponential function including the rotation of the atmosphere. We obtained a first-order analytic atmospheric drag perturbations. We evaluated the short periodic perturbations of all orbital elements due to the geopotential. To obtain a more precise treatment, we integrated the equations of motion numerically using the Runge-Kutta method of fourth order. The validity of the theory is checked using the already decayed Indian satellite ROHINI where its data are available. A good agreement between the observed line elements of the satellite and our predicted solutions encourage us in future to apply our theory on any similar future missions.



---

*Source: 489849-2012-02-12.xml*