Correction - Orbite d'un satellite dans le champs gravitationnel terrestre

Il est clair qu'en première approximation la trajectoire d'un satellite artificiel est une ellipse keplerienne dont l'un des foyers est occupé par la Terre. Lorsque l'on y réfléchit bien de nombreux phénomènes peuvent néanmoins entrer en ligne de compte pour venir perturber cette trajectoire ''naturelle'' : La Terre n'est ni sphérique ni formée d'une succession de couches concentriques, elle est dotée d'une atmosphère relativement dense qui peut contribuer à freiner les satellites surtout si ces derniers ont un périgée à faible altitude, enfin la pression de radiation, le vent solaire, le flux de particules cosmiques et les éventuels chocs avec tout un tas d'objets plus ou moins naturels doivent être pris en compte lors de la coûteuse mise en station de satellite artificiel. De tous ces effets la non sphéricité de la Terre et la déformation du champs de gravitation induite est le principal. Nous allons donc voir comment faire entrer ce dernier dans le modèle de la théorie planétaire de Lagrange.
 

Le véritable potentiel gravitationnel terrestre

Soit $Oxyz$, un référentiel d'origine $O$ le centre de la terre, et $\mathbf{r}$ le vecteur indiquant la position du satellite (supposé ponctuel) de masse $m$ dont nous étudions le mouvement. Soit$\mathbf{r}^{\prime}$ le vecteur indiquant la position d'un élément quelconque de masse $dm^{\prime}$ de la Terre de volume $V$ et de masse $$m_{\oplus}=%%{\displaystyle\int_{\mathbf{r}^{\prime}\in V}}dm^{\prime}\quad,$$ et $\delta$ la distance entre le satellite et $dm^{\prime}$. (voir figure 1)

L'élément de potentiel d'interaction s'écrit

$$dU=-G\frac{mdm^{\prime}}{\delta}=-G\frac{mdm^{\prime}}{\left\vert \mathbf{r}%%-\mathbf{r}^{\prime}\right\vert }$$ si $V$ est le volume occupé par la Terre, le potentiel de gravitation duquel dérive la force qu'exerce notre planète sur le satellite est donc

$$U\left( \mathbf{r}\right) =-G{\displaystyle\int_{\mathbf{r}^{\pri......\frac{mdm^{\prime}}{\left\vert \mathbf{r}-\mathbf{r}^{\prime}\right\vert }%%$$ (1)

Soit $\psi$ l'angle entre $\mathbf{r\,}$et $\mathbf{r}^{\prime}$, on sait que$\delta^{2}=r^{2}+r^{\prime2}-2rr^{\prime}\cos\psi$ autrement dit

$$\frac{1}{\delta}=\frac{1}{\left\vert \mathbf{r}-\mathbf{r}^{\prim......}}{r}\right) \left(\left( \dfrac{r^{\prime}}{r}\right) -2\cos\psi\right) }}$$ en supposant que le satellite évolue toujours suffisamment loin de la Terre, c'est à dire $r\gg r^{\prime}$, on peut utiliser un développement limité pour l'expression de l'inverse de $\delta$, on trouve $$%%\begin{array}[c]{cl}%%\dfrac{1}{\delta}= & \dfrac{1}{r}......e}}{r}\right) -2\cos\psi\right)^{3}+...\right)\end{array}$$ $$%%\begin{array}[c]{cl}%%\Leftrightarrow\quad\dfrac{1}{\de......\dfrac{r^{\prime}}{r}\right) ^{5}\right) \right)\end{array}$$ une méthode moins expérimentale permet alors de voir que les coefficients des puissances de $r^{\prime}/r$ dans le développement ci-dessus ne sont autres que des polynômes de Legendre de première espèce $$%%\begin{array}[c]{lll}%%P_{0}\left( \cos\psi\right) =1 &......\frac{d^{n}}{dx^{n}}\left( x^{2}-1\right) ^{n}%%\end{array}$$ on montre ainsi que $$\frac{1}{\delta}=\frac{1}{\left\vert \mathbf{r}-\mathbf{r}^{\prim......^{\infty}\left( \dfrac{r^{\prime}}{r}\right)^{n}P_{n}\left( \cos\psi\right)$$ le potentiel terrestre $\left( \ref{potterre}\right)$ devient donc $$U\left( \mathbf{r}\right) =\sum_{n=0}^{\infty}U_{n}\left( \mathbf{r}%%\right)$$   avec $$U_{n}\left( \mathbf{r}\right) =-\frac{Gm}%%{r}\left\{ \frac{1}{......^{\prime}\in V}}r^{\prime n}P_{n}\left( \cos\psi\right) dm^{\prime}\right\}$$ Etudions à présent la nature des premiers termes de ce développement. Si $n=0$ on voit immédiatement que $$U_{0}\left( \mathbf{r}\right) =-\frac{Gm}{r}%%{\displaystyle\int_{r^{\prime}\in V}}dm^{\prime}=-\frac{Gmm_{\oplus}}{r}$$ est le potentiel de gravitation du problème non perturbé, lorsque la Terre est assimilée à un point de masse $m_{\oplus}$ ou bien constituée d'une succession de couches sphériques homogènes concentriques. Lorsque $n=1$, en utilisant la relation $$\cos\psi=\frac{\mathbf{r}.\mathbf{r}^{\prime}}{r\,r^{\prime}}=\frac{xx^{\prime}+yy^{\prime}+zz^{\prime}}{r\,r^{\prime}}$$ nous obtenons $$U_{1}\left( \mathbf{r}\right) =-\frac{Gm}{r^{3}}\left\{ x{\displa......y^{\prime}dm^{\prime}+z{\displaystyle\int_{V}}z^{\prime}dm^{\prime}\right\}$$ soit $$U_{1}\left( \mathbf{r}\right) =-\frac{Gmm_{\oplus}}{r^{3}}\left\{xx_{G}^{\prime}+yy_{G}^{\prime}+zz_{G}^{\prime}\right\}$$ où $\left( x_{G}^{\prime},y_{G}^{\prime},z_{G}^{\prime}\right)$ sont les coordonnées du centre de gravité de la Terre. Si l'on choisit ce point (supposé fixe) comme origine de notre référentiel, nous avons par conséquent $U_{1}\left( \mathbf{r}\right) =0$. Le cas $n=2$ est quant à lui plus complexe. Nous avons en effet $$U_{2}\left( \mathbf{r}\right) =-\frac{Gm}{r^{3}}\left\{{\displa......}\right) ^{2}}{\left( rr^{\prime}\right) ^{2}%%}\right) dm^{\prime}\right\}$$ en développant le carré nous avons donc (en posant $x_{\mu}%%(\mu=1,2,3)=x,y,z$ et $x_{\mu}^{\prime}(\mu=1,2,3)=x^{\prime},y^{\prime},z^{\prime}$) $$U_{2}\left( \mathbf{r}\right) =-\frac{Gm}{r^{3}}\left\{ -\frac{1}......{\displaystyle\int_{V}}x_{\mu}^{\prime}x_{\nu}^{\prime}dm^{\prime}\right\}$$ en introduisant les moments $I_{\mu\mu}$ et produits $I_{\mu\nu}(\mu\neq\nu)$ d'inertie de la Terre par rapport aux axes $Oxyz$ $$I_{\mu\mu}=%%{\displaystyle\int_{V}}\left( r^{\prime2}-x_{\mu}^{\prime2}\right) dm^{\prime}$$   et $$I_{\mu\nu}=%%{\displaystyle\int_{V}}x_{\mu}^{\prime}x_{\nu}^{\prime}\,dm^{\prime}$$ le troisième terme du développement du potentiel terrestre devient $$U_{2}\left( \mathbf{r}\right) =-\frac{Gm}{r^{3}}\left\{ -\tfrac{1......mu\mu}\right)+\sum_{\nu\neq\mu}^{3}x_{\mu}x_{\nu}I_{\mu\nu}\right] \right\}$$ On impose alors aux axes $Ox$, $Oy$ et $Oz$ de notre référentiel d'être confondus avec les axes principaux d'inertie de la Terre (supposés fixes), en pratique et en première approximation on confond donc $Oxy$ avec le plan équatorial terrestre et $Oz$ avec l'axe polaire. En conséquence, tout les produits d'inertie s'annulent et en posant $$A=%%{\displaystyle\int_{V}}\left( r^{\prime2}-x_{1}^{\prime2}......%{\displaystyle\int_{V}}\left( y^{\prime2}+z^{\prime2}\right) dm^{\prime}$$ $$B=%%{\displaystyle\int_{V}}\left( r^{\prime2}-x_{2}^{\prime2}......%{\displaystyle\int_{V}}\left( x^{\prime2}+z^{\prime2}\right) dm^{\prime}$$ $$C=%%{\displaystyle\int_{V}}\left( r^{\prime2}-x_{3}^{\prime2}......%{\displaystyle\int_{V}}\left( x^{\prime2}+y^{\prime2}\right) dm^{\prime}$$ il vient

$$U_{2}\left( \mathbf{r}\right) =-\frac{Gm}{r^{3}}$$

$$\left\{ -\tfrac{\left( A+B+C\right) }{4}+\tfrac{3\left[ x^{2}\lef......+x^{2}\left( \left( \frac{A+B+C}{2}\right) -C\right) \right] }{2r^{2}}\right\}$$

soit

$$U_{2}\left( \mathbf{r}\right) =-\frac{Gm}{r^{3}}\left\{ \frac{\le......C\right) }{2}-\frac{3\left( Ax^{2}+By^{2}+Cx^{2}\right) }{2r^{2}%%}\right\}$$ on introduit alors les coordonnées équatoriales $\left(r,\lambda,\phi\right)$, où $\lambda$ est appellé longitude et $\phi$ latitude, (voir figure 2)

le potentiel s'écrit

$$%%\begin{array}[c]{ll}%%U_{2}\left( \mathbf{r}\right) & =......A-B\right)\cos^{2}\phi\cos2\lambda}{4}\right\}\end{array}$$ Si finalement on admet que la Terre présente une symétrie de révolution autour de son axe polaire, nous avons $A=B$ et le premier terme non nul de perturbation du potentiel ne dépend que du rayon et de la latitude et s'écrit

$$U_{2}\left( \mathbf{r}\right) =U_{2}\left( r,\phi\right) =-\frac{......>120 {r^{3}}\left( C-A\right) \left( \frac{1}{2}-\frac{3\sin^{2}\phi}{2}\right)$$ (2)

en introduisant la constante

$$J_{2}=\frac{C-A}{m_{\oplus}r_{e}^{2}}$$ où $m_{\oplus}=5,976\,\,10^{24}$Kg et $r_{e}=6378$ Km sont respectivement la masse et le rayon équatorial de la Terre, on obtient (le signe $-$ étant absorbé par le polynôme de Legendre) $$U_{2}\left( r,\phi\right) =\frac{Gmm_{\oplus}}{r}\left( \frac{r_{e}}%%{r}\right) ^{2}J_{2}P_{2}\left( \sin\phi\right)$$ dans ce même cadre d'approximation, on peut montrer (on pourra consulter par exemple l'article de Bretagnon et al.  pour avoir un apercu global de ce genre de problème) que ce résultat n'est pas fortuit et qu'il est le cas particulier $n=2$ d'une relation générale pour $n\geq2$ de la forme $$U_{n}\left( r,\phi\right) =\frac{Gmm_{\oplus}}{r}\left( \frac{r_{e}}%%{r}\right) ^{n}J_{n}P_{n}\left( \sin\phi\right)$$ le potentiel de gravitation de notre modélisation de la Terre s'écrit donc $$U\left( r,\phi\right) =-\frac{Gmm_{\oplus}}{r}\left\{ 1-\sum_{n=2......}\left( \frac{r_{e}}{r}\right) ^{n}J_{n}P_{n}\left( \sin\phi\right)\right\}$$ les coefficients $J_{n}$ peuvent être déduits d'une théorie plus ou moins approximative ou directement issus de l'observation du mouvement de satellites artificiels, on trouve (dans [#!bel!#] par exemple) $$%%\begin{array}[c]{c}%%J_{2}=0,0010826\\J_{3}=-2,5\,\,10^{-6}\\J_{4}=-1,6\,\,10^{-6}%%\end{array}$$ La force de gravitation exercée par la Terre sur le satellite est donc $$\overrightarrow{F}=-\mathrm{grad}_{\mathbf{r}}U\left( r,\phi\righ......}\left( \frac{r_{e}}{r}\right) ^{n}J_{n}P_{n}\left( \sin\phi\right)\right]$$ le principe fondamental de la dynamique donne alors

$$\dfrac{d^{2}\mathbf{r}}{dt^{2}}=-Gm_{\oplus}\dfrac{\mathbf{r}}{\l......ert \mathbf{r}\right\vert ^{3}}+\mathrm{grad}_{\mathbf{r}}\left( R\right)%%$$ (3)

où l'on voit sourdre la fonction perturbatrice du problème

$$R=R\left( r,\phi\right) =-\frac{Gm_{\oplus}}{r}\sum_{n=2}^{\infty}\left(\frac{r_{e}}{r}\right) ^{n}J_{n}P_{n}\left( \sin\phi\right)$$

Expression des éléments osculateurs

Comme le montre la relation $\left( \ref{probsat}\right)$, si $R$ est petit devant le potentiel à deux corps, le problème du mouvement d'un satellite dans le champ de gravitation de la Terre non ponctuelle peut se traiter par la théorie planétaire de Lagrange. Il ne reste plus qu'à exprimer la fonction perturbatrice en fonction des éléments elliptiques.

Comme on peut le voir sur la figure $\left( \ref{movsat}\right)$, si le mouvement non perturbé est repéré par rapport au plan équatorial, l'inclinaison $i$ de l'orbite du satellite par rapport au plan équatorial terrestre, l'argument $\omega$ de son périgée et son anomalie vraie $f$ vérifient la relation 

$$\sin\phi=\sin i\sin\left( \omega+f\right)$$ En ne considérant que le premier terme ($n=2$)de la fonction perturbatrice, nous avons donc  $$R=\frac{GJ_{2}m_{\oplus}r_{e}^{2}}{2}\left( \frac{1}{r}\right) ^{3}\left(1-3\sin^{2}i\sin^{2}\left( \omega+f\right) \right)$$ Comme dans la théorie de la Lune nous allons décomposer la fonction perturbatrice $R$ en une partie séculaire $\overline{R}$ et une partie périodique $\widetilde{R}$ telles que  $$\widetilde{R}=R-\overline{R}$$   avec $$\overline{R}=\frac{2\pi}%%{T}\frac{1}{2\pi}%%{\displaystyle\i......{0}^{2\pi}}Rdt=\frac{n}{2\pi}%%{\displaystyle\int\limits_{0}^{2\pi}}Rdt$$ en écrivant la forme différentielle de la loi des aires (cf relation $\left( \ref{loiaire}\right)$) nous avons  $$dt=\frac{r^{2}}{C}df=\frac{r^{2}}{\sqrt{p\mu}}df=\frac{1}{\sqrt{\mu}}%%\frac{r^{2}}{\sqrt{a\left( 1-e^{2}\right) }}df$$ ainsi  $$\overline{R}=\frac{nGJ_{2}m_{\oplus}r_{e}^{2}}{2\sqrt{\mu}}\frac{......pi}}\frac{1}{r}\left( 1-3\sin^{2}i\sin^{2}\left( \omega+f\right) \right) df$$ et comme (trajectoire elliptique)  $$\frac{1}{r}=\frac{1+e\cos\left( f\right) }{p}=\frac{1+e\cos\left( f\right)}{a\left( 1-e^{2}\right) }$$ nous avons finalement  $$\overline{R}=\frac{nGJ_{2}m_{\oplus}r_{e}^{2}}{2\sqrt{\mu}}\frac{......\right) \right] \left[ 1-3\sin^{2}i\sin^{2}\left(\omega+f\right) \right] df$$ c'est à dire après intégration  $$\overline{R}=\frac{nGm_{\oplus}}{\sqrt{\mu}}\frac{J_{2}r_{e}^{2}}......c{\left( 3\cos^{2}i-1\right) }{\left[ a\left(1-e^{2}\right) \right] ^{3/2}}$$ en posant  $$\gamma=n\frac{Gm_{\oplus}}{\sqrt{\mu}}\frac{J_{2}}{4}r_{e}^{2}$$ nous remarquons que $\overline{R}$ ne dépend ni de $\Omega$ ni de $\omega$ ni de $\tau$, on en déduit automatiquement  $$\left( \frac{d\overline{a}}{dt}\right) =\left( \frac{d\overline{e}}%%{dt}\right) =\left( \frac{d\overline{i}}{dt}\right) =0$$ le demi grand axe, l'exentricité et l'inclinaison ne présentent donc pas de variation séculaire (dans notre cadre d'approximation). Le mouvement séculaire du périgée est décrit par l'équation différentielle  $$%%\begin{array}[c]{lll}%%\dfrac{d\overline{\omega}}{dt} &......=\dfrac{2\pi}%%{T}=\sqrt{\dfrac{\mu}{a^{3}}}%%\end{array}$$ on a donc  $$\overline{\omega}=\alpha t+\beta$$ avec $$\alpha=\dfrac{\gamma}{\sqrt{\mu}}\dfrac{3\left( 5\cos^{2}i-1\righ......t) }\right] ^{2}\frac{m_{\oplus}}{m+m_{\oplus}%%}\left( 5\cos^{2}i-1\right)$$ De même, le mouvement séculaire du noeud ascendant est décrit par l'équation  $$%%\begin{array}[c]{lll}%%\dfrac{d\overline{\Omega}}{dt} &......rac{6\cos i}{a^{2}\left( 1-e^{2}\right)^{2}}%%\end{array}$$ on a donc $$\overline{\Omega}=\lambda t+\delta$$ avec  $$\lambda=-\dfrac{\gamma}{\sqrt{\mu}}\dfrac{6\cos i}{a^{2}\left( 1-......}}{a\left( 1-e^{2}\right)}\right] ^{2}\frac{m_{\oplus}}{m+m_{\oplus}}\cos i$$ Les perturbations de l'orbite sont dans notre cadre d'approximation, de même nature que celle trouvées pour notre théorie de la Lune. Pour un satellite de masse $m=1000Kg$, sur une orbite inclinée de $i=45{{}^{\circ}}$ d'excentricité $e=0.01$ et de demi-grand axe $a=2,5$$r_{e}$ on trouve une période de  $$T=2\pi\frac{a^{3/2}}{\sqrt{G\left( m+m_{\oplus}\right) }}=20034s=5$$h$$%%33^{\prime}54^{\prime\prime}$$ correspondant à un moyen mouvement de  $$n=\frac{360{{}^{\circ}}}{T}=1552{{}^{\circ}}33^{\prime}38.32^{\prime\prime}\text{\textsf{par jour soit 4,3 tours journaliers.}}$$ ainsi la précession du périgée est de $\alpha=0{{}^{\circ}%%}18^{\prime}8,5^{\prime\prime}$ par jour et la rétrogradation du noeud ascendant de $\lambda=-0{{}^{\circ}}17^{\prime}6,45^{\prime\prime}$ par jour. Ce type d'analyse a été rendu nécessaire par la mise en orbite des premiers satellites : Les premiers articles sur ce sujet sont russes et datent de la fin des années 50. Après l'applatissement de la Terre, la deuxième source de perturbation à prendre en compte est l'effet de freinage dû à son l'atmosphère. Sur de petites échelles de temps (quelques dizaines de révolutions), les résultats que nous avons obtenu sont malgré tout en assez bon accord avec les mesures effectuées sur les satellites. 
Jérôme Perez 26/09/04

Références
G. Pascoli, Eléments de mécanique céleste, Armand Colin, 1993
P.Bretagnon, P. Rocher, J.L.\Simon, Theory of the rotation of the rigid earth, Astronomy \&Astrophysics
V. Beletski, 1977, Essai sur le mouvement des corps cosmiques, Editions Mir, Moscou