Dans les rapport précedents l'itérêt était d'étudier la convergence et la stabilitée de différentes schémas d'intégration numérique, pour
en tester la qualité et les différentes caractéristiques. Les resultats précedents nous on permit de déduire par exemple,
que le schéma de Runge-Kutta est assez bon et peut être implementé pour obtenir un ordre de convergence assez élévé.
Comme en physique les problèmes plus intéressantes ne peuvent pas être étudiés de façons purement analytique, un schéma numérique
performat est strictement nécessaire pour l'étude de ces problèmes. Dans le rapport qui suit sera étudié un problème de gravitation dans lequel trois
objets massiques interagissent entre eux. Le problème ne peut pas être étudié analytiquement, et donc le Schéma de Runge-Kutta d'ordre 4 sera implementé pour l'utiliser dans l'étude de ce problème. Le rapport en soit illustre différentes façons d'exploiter un tel schéma numérique.
\section{Analyse du problème}
\subsection{Description du problème}
L'intérêt est porté sur le mouvement de N ($N \leq3$) corps soumis aux forces de gravitation.
Comme le nombre d'objets est maximum 3, ils sont tous coplanaires, et on simplifie le problème en supposant qu'il sont contraint à
bouger sur un plan bidimensionnel. Ainsi on s'intèresse donc seulement aux coordonnée planaires $xOy$ du système.
Les forces agissantes dans le système sont :
\begin{itemize}
\item La force gravitationnelle exércée par le corps j sur le corps i :
où T représente la Terre, et $j$ l'objet subissant la force de trainée, et où
\begin{equation*}
\rho(r)=\rho_0 \exp{-(r-R_T)/\lambda}
\end{equation*}
\end{itemize}
Dans ces rélation $G=6.674\cdot10^{-11}$\si{\metre^3 \kilogram^{-1}\second^{-2}} est la constante de gravitation universelle, $\rho_0$
la valeur de $\rho(r)$ au niveau de la mer, $R_T =6.3781\cdot10^6$ le rayon de la Terre, $\lambda$ une épaisseur caractéristique, $S$ la surface
de section de l'objet j, $C_x$ le coefficient de traînée.
\subsection{Horbite régulière d'un objet autour de la terre}\label{horbitenormale}
Ici l'intérêt est celui de connaitre quelles doivent être les position et vitesse initiales d'un objet qui veut garder une horbite constante autour de la terre à
une distance minimale $h =10000$\si{\metre} de la surface de la Terre, et donc à $h_{min} = R_T + h$ du centre de la Terre. Ceci en connaissant
le module de la vitesse initiale $v_0$ et la distance initiale du centre de la Terre $r_0$.
Ici on suppose que la masse de l'objet $m_r \ll m_T$ où $m_T$ est la masse de la terre, donc l'action de gravitation subit par la Terre de la part de l'objet est négligeable. Dans ce cas la force de traînées n'est pas considérée, donc par conservation de l'énérgie mécanique la rélation suivante est obtenue :
Une ultérieure supposition est que l'objet, bougeant sur une horbite elliptique, a une vitesse exactement perpendiculaire au rayon liant l'objet au centre de la Terre. Par la conservation du moment angulaire la rélation suivante est détérminé :
où $\theta$ est l'angle entre la vitesse initiale $\vec{v_0}$ et le vecteur de position initiale $\vec{r_0}$ de l'objet.
L'équation (\ref{momang}) implique donc
\begin{equation}\label{teta}
\theta = \arcsin{\frac{h_{min}v_{max}}{r_0v_0}}
\end{equation}
Par le resultat en équation (\ref{teta}) il peut être déduit que le vecteur de vitesse initiale, si celui de position initiale $\vec{r_0}=(0,r_0)$, est donné par :
\begin{equation}\label{v0}
\vec{v_0} = v_0
\begin{pmatrix}
\sin{\theta}\\
-\cos{\theta}
\end{pmatrix}
\end{equation}
Voir Figure (\ref{}).
\subsection{Système Terre-Lune}\label{terrelune}
Ici, il est intéressant de considèrer le système binaire Terre-Lune, telle qu'il reste en mouvement circulaire uniforme par rapport au centre de gravité G.
Dans le reférentiel inertiel qui voit son origine en G, par définition du centre de gravité, les position de la Terre et de la Lune $\vec{r_T}$ et $\vec{r_L}$ sont liés par la rélation :
Ensuite pour détérminer la vitesse angulaire $\Omega$ du système il suffit d'appliquer le deuxième principe de newton en considèrant l'accéleration angulaire
$a_{a,\alpha} =\omega_{\alpha}^2r_{\alpha}$ et que la vitesse angulaire $\omega_\alpha=\Omega$ pour les deux corps $\alpha=$ T,L.
De ce resultat les modules des vitesses $v_T$ et $v_L$ respectivement de la Terre et de la Lune, qui sont opposée pour que la distance $d = cste$, sont données par :
\begin{equation}\label{vitesses}
v_T = \Omega r_T \qquad v_L = \Omega r_L
\end{equation}
Pour finir, une période de révolution $T$ pour ce système est donné par la rélation $T =2\pi/\Omega$.
\subsection{Points de Lagrange}
\begin{minipage}{\textwidth}
\begin{minipage}{0.6\textwidth}
Comme montré en cours, dans un {\it Système à 3 corps limité}, c'est à dire dans un système
où la masse du troisième corps $m_3\ll m_1,\, m_2$, existent des point d'équilibre par rapports
aux deux premièrs corps, pourlequels si l'on y place le corps
de masse plus basse, ceci reste dans ce point si l'on suit le référentiel binaire des corps 1 et 2.
Deux de ces points sont les Point de lagrange, caractérisées par le fait que le système tertiaire lorsque
le troisème corps est placé en ce point, décrit un triangle équilateral.
Pour le système Terre Lune donc, en se plaçant dans le référentiel R' tournant avec vitesse angulaire $\Omega$ (ici on prend les constantes de la séction précedente) et centré en G, les deux points de Lagranges peuvent être déduit par une simple application du théorème de pitaghore.
\end{minipage}
\hspace{0.02\textwidth}
\begin{minipage}{0.48\textwidth}
\includegraphics[scale=0.8]{schema_lagrange.PNG}
\captionof{figure}{Schématisation géométrique pour l'étude du {\it Point de Lagrange}\cite{ref1}}
\label{schemalagrange}
\end{minipage}
\end{minipage}
Le premier point de Lagrange est donné, dans le referentiel R', par la rélation suivante :
\begin{equation}\label{lagrange}
\vec{P_0} =
\begin{pmatrix}
d(\frac{1}{2} - \frac{m_L}{m_T + m_L})\\
\\
d\frac{\sqrt{3}}{2}
\end{pmatrix}
\end{equation}
où $d$ est la distance Terre-Lune.
Pour lancer une simulation correctement il faut que le corps situé en $P_0$ tourne correctement avec le système binaire Terre-Lune.
Il faut donc que sa vitesse initiale soit $\vec{v_0} =\Omega r \vec{e_\theta}$, où $\Omega$ est la vitesse angulaire du système, et
$r$ la norme de $\vec{P_0}$ (voir Figure (\ref{schemalagrange})).
$r$ est simplement donné par la rélation suivante :
\begin{equation}\label{rlagrange}
r = \sqrt{d^2(\frac{1}{2} - \frac{m_L}{m_T + m_L})^2 +d^2\frac{3}{4}} = d \sqrt{(\frac{m_L}{m_T+m_L})^2 - \frac{m_L}{m_T+m_L} +1}
\end{equation}
Pour finir, il est necéssaire de connaitre l'angle $\theta$ entre l'orizontale et $\vec{P_0}$ pour en déduire la direction $\vec{e_\theta}$.
Où $y_1$ et $y_2$ sont respectivement les résultat de la simulation pour un pas de temps où deux demi-pas de temps $dt$ et donc $|y_1- y_2|$ est l'erreur sur entre les deux.
{\bf Schema(}$dt_{old}${\bf)} modifie $(y_1,y_2)$ pour que les deux soyent les résultat de la simulation pour un nouveu $dt_{old}$ et deux $dt_{old}$.
$n$ est l'ordre de convergence du schéma numérique utilisé pour la simulation. $dt_{old}$ et $dt_{new}$ sont respectivement le pas de temps avant et après modification. En fin $y$ est la valeur prise par la mesuration simulée pour ce point de la simulation.
\section{Simulations Numériques}
Le schéma numérique utilisé dans les simulations est celui de Runge-Kutta d'ordre 4, qui converge en effet à l'ordre 4.
\subsection{Trajectoire Elliptique}
En première étude il est interéssant de verifier par une simulation numérique sur les donnée trouvée en section (\ref{horbitenormale}), que les conditions
imposées soyents respéctées.
Pour commencer une simulation de la trajectoire de l'Apollo 13 est lancée avec le centre de la terre situé en origine.
La masse de l'Apollo est $m_A =5809$\si{\kilogram}, à un distance verticale $r_0=3.14159\cdot10^8$\si{\metre} et une vitesse initiale de module $v_0=1.2\cdot10^3$\si{\metre\per\second} dont le vecteur est donné par les relation (\ref{v0}) et (\ref{teta}). Le temps totale de simulation est
$t_{fin} =2$ jours et la masse de la terre est $m_T =5.972\cdot10^{22}$.
%TODO inserire grafico Orbita
Le graphe en Figure (\ref{}) montre que la trajectoire simulée est en effet elliptique, en tout respect des lois de Kepler.
De plus le graphe montre aussi que 2 jours ne sont pas assez pour parcurir une révolution complète, et de plus comme l'Apollo s'éloigne de la Terre,
sa vitesse dinimue, et il prendra plus de temps pour parcourir le reste de la trajectoire.
Ensuite un étude de convergence est porté sur l'erreur des valeurs simulés par
rapport à celles cherchés analityquement, de $h_{min}$ et $v_{max}$ donnée par la rélation (\ref{vmax}) en diminuant le pas de temps.
%TODO inserire grafico Convergenze
Les graphes en Figure (\ref{}) et (\ref{}) montrent que le schéma converge en $dt$ et de plus les erreus sur les valeurs convergent à 0.
L'ordre de convergence, qui est donné par la pente de la droite de convergence en graphe loglog, est en effet 4, qui l'ordre de convergence
effectif du Schéma utilisé.
Pour finir la même simulation est lancée avec un pas de temps adaptif, est un étude de convergence sur la valeur de $\epsilon$ est mené.
%TODO inserire grafico convergenza su epsilon
L'ordre de convergence de Runge-Kutta sur $\epsilon$ pour un temps adaptif est %TODO stimare ordine di convergenza
Il est donné par la pente de la droite de convergence.
Un dernier étude interessant est celui de comparer le nombre de pas de temps utilisé dans le pas de temps adaptif et ceux pour un pas de temps fixe.
De plus, il est intéréssant de voir comment varie le pas de temps en fonction de la distance que l'Apollo 13 a par rapport à la Terre.
%TODO confrontare numero di passi fatti per dt fisso e per adattativo e commentare
%TODO inserire grafico dt adattativo in fuzione della distanza
Comme c'est possible de voir sur le graphe en Figure (\ref{}) le $dt$ est très petit lorsque la distance Terre-Apollo est petite. Ceci est du au fait que c'est en ces points que l'Apollo subit le plus d'accèlèration, et que donc si l'on prend un pas de temps lègèrement plus grands le resultat change beaucoup.
Comme en chaque pas de temps l'objet simulé par vers la tangente de la trajectoire, c'est dans les point où cette dernière change le plus que le pas de temps doit rester petit pour que une certaine précision soit respéctée.
\subsection{Atterrissage}
Pour cette simulation la force de trianée de l'air est prise en considération et agit sur Apollo 13 qui va donc atterir sur Terre sous l'action de freinage de cette force. Le but étant de trouver une direction initiale de l'Apollo pour laquelle ce dernier puisse atterir mais en subissant une accèlèration maximale qui soit le moindre possible. En donné il y a: $\rho_0=1.2$\si{\kilogram\per\metre^3}, $\lambda=7238.2$\si{\metre}, le diamètre de l'Apollo, par lequel est mésurée sa surface, est $d =3.9$\si{\metre} et $C_x =0.3$. Ainsi la force de trainée de l'air sur l'Apollo est calculée par la rélation (\ref{friction}).
Dans un premier temps, on considère la même vitesse initiale que celle prise dans la section précédente et une simulation est lancée.
La trajectoire de l'Apollo est illustré sur le graphe en Figure (\ref{}), et il peut être observé que la force de trainée est assez ample pour faire tomber à pique
l'Apollo vers le centre de la terre. En effet l'on voit que la trajectoire converge vers l'origine, où se trouve le centre de la Terre.
%TODO aggiungere considerazione sull'accellerazione massimale e potenza dell'attrito
%TODO studio di convergenza
Dans un deuxième temps, plusieures simulations pour des angles $\theta$ différent injéctés en (\ref{v0}) sont lancées. De ces simulation
l'acceleration maximale est collectée et progetée sur un graphe en fonction de cet angle $\theta$. Le graphe suit illustre cette courbe :
%TODO inserire grafico accellerazione minimale
L'analyse du graphe en Figure (\ref{}) aide dans la recherche de l'angle optimale pour la minimisation de l'accélération maximale $a_{max}$ subite par l'Apollo, et donc par les astronaute. Cette courbe tends à 0 pour des angles plus amples que celui défini en (\ref{teta}), mais ceci parceque en effet pour ces angles
l'Apollo ne rentre jamais dans l'atmosphère terrestre. Le minimum de $a_{max}$ est à chercher donc à l'intérieur du domaine où $a_{max}$ n'approche pas 0. %TODO continuare commento...
\subsection{Terre-Lune}
Pour tester les condition initiales sur les corps Terre et Lune trouvées en section \ref{terrelune}, une simulation avec ces conditions est lancée.
La masse de la lune est $m_L =7.3477\cdot10^22$\si{\kilogram} et la distance Terre-Lune $d =3.84748\cdot10^8$\si{\metre}.
%TODO inserire grafico terraluna
Le graphe en Figure (\ref{}) montre la trjectoire des deux corps dans les référentiel inertiel ayant pour origine le centre de masse G du système.
L'orbite décrite par les deux corps est en effet circulaire, et les graphes qui suivent
montrent que les gradeur comme $d$ (voir Figure (\ref{})) ou l'énérgie mécanique du système $E_m$ (voir Figure (\ref{})), reste constant. En effet
ces valeurs oscillents, mais l'amplitude de l'oscillation est d'un ordre de grandeur assez petit pour que les variations sur $d$ et $E_m$ soyent
négligéable.
% TODO inserire grafici conservazione di d et E_m e quantità di moto
% TODO commentare la non conservazione della quantità di moto?
\subsection{L'Apollo entre la Terre et la Lune}
\subsection{Point de Lagrange}
\section{Conclusions}
\section{Annexes}
\subsection{Changement de Référentiel}
De suite est présenté un algorithme permettant de convertir des coordonnée cartésiennes dans un référentiel inertiel R à origine O
dans les cordonnées cartésiennes d'un référentiel non inertiel R' à origine O et tournant à vitesse angulaire constante $\Omega$.
L'algoritme s'applique lorsque le R' et R sont confondus à temps t = 0, en connaissant $\Omega, t$ et $(x,y)$ c'est à dire les coordonnées
que l'on veut convertir.
\begin{tabular}{|l|}
\hline
{\bf Algorithme}\\
\hline
Entrées : $x, y, t, \Omega$\\
Sortie : $x', y'$\\
\hline
$r =\sqrt{x^2+ y^2}$
\\
\\
$\;\;${\bf Si} : $y >0$\\
$\qquad\theta=\arccos{\frac{x}{r}}$
\\
$\;\;${\bf Si} : $y \leq0$\\
$\qquad\theta=2\pi-\arccos{\frac{x}{r}}$\\
\\
$\alpha=\theta-\Omega t$\\
\\
\hline
$x' = r \cos{\alpha}$
\\
$y' = r \sin{\alpha}$
\\
\hline
\end{tabular}
\begin{thebibliography}{99}
\bibitem{ref1}
ECOLE POLYTECNIQUE FEDERALE DE LAUSANNE, Semestre d'automne 2018, Physique Numérique I - Exercice 4.
\bibitem{ref2}
L. VILLARD, PHYSIQUE NUMERIQUE I - II ,Swiss Plasma Center, Faculté des Sciences de Base, Section de Physique, Ecole Polytecnique Fédérale de Lausanne.
\bibitem{ref3}
L. VILLARD, Cours de Physique Numérique pour physiciens, Faculté des Sciences de Base, Section de Physique, Ecole Polytecnique Fédérale de Lausanne.