-- Cours/Projet - AERO-I3 --

Mécanique Spatiale

- Module Orbiter -

Utilisation et documentation

simon.marie@lecnam.net
In [2]:
# Importation de la classe capsule
from Orbiter import Capsule,Stage

%matplotlib inline  
plt.rc('xtick',labelsize=14)
plt.rc('ytick',labelsize=14)
In [3]:
# Constante célestes:                        
G=6.67408e-11 # Constante de gravitation -   
Rt=6380000 # Rayon terrestre - m             
Rl=1736000 # Rayon luniare - m               
Ms=1.9884e30 # Masse du soleil - kg          
Mt=5.9722e24 # Masse de la terre - kg        
Ml=7.3477e22  # Masse Lune - kg              
al=384e6      # Distance Terre-Lune - m      
at=149.598262e9  # Distance Terre-Soleil - m 

Le module Orbiter

Le module Orbiter permet la résolution simplifiée des équations du mouvement d'un objet spatial soumis à une force de gravité et à une force de poussée. Il tient également compte des effets atmosphériques au décollage et des flux thermique lors d'une rentrée atmosphérique. Le module contient 3 classes d'utilisation différentes:

  1. La class *capsule* : Elle permet l'integration du mouvement d'un corps celeste sans propulsion dans le référentiel géocentrique 2D ($r,\theta)$. Les maneuvres sont gérée via une méthode $mvr()$ permettant d'imposer un $\Delta v$.
  2. La class *Stage* : Elle permet l'intégration du mouvement d'un élément de lanceur spatial ayant des caractéristiques propulsives données ($ISP$, $qm$, $k_s$). Les $\Delta V$ doivent ici être transmis via des temps d'allumage moteurs.
  3. La class $Lanceur$ : Elle permet de construire un lanceur spatial complet constitué de plusieurs étages ($Stage$).

Vous pouvez trouver en ligne la documentation complète du module Orbiter. Ce notebook présente l'utilisation des deux première class. L'utilisation de la class Lanceur est présenté dans le notebook consacré au lanceur SaturnV.

La class Capsule

L'utilisation de la class capsule permet de se familiariser avec les trajectoires spatiales et les manoeuvres. Ici on se propose de placer une capsule sur une orbite initiale circulaire à 400km d'altitude. La class prend comme argument initiaux, l'altitude et la vitesse $r\dot{\theta}$ initiale:

In [4]:
zs=5000e3
us=np.sqrt(G*Mt/(Rt+zs))
cnap=Capsule(z0=zs,u=us) # cnap est le nom de notre capsule
print(cnap)
----   29, Oct 2019 | 17:32   ----
time= 0.0
----Position--(/earth/)-------
u=5918.22871188
z=5000000.0
theta=0.0
Assiette=0.0 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.005
finesse=0.0

Si on veut regarder l'évolution de la capsule pendant une periode entière de révolution, on intégre la trajectoire sur une periode avec la méthode run():

In [5]:
# Periode de l'orbite à zs:
T1=2*np.pi*(Rt+zs)/us
cnap.run(dt=1,ntmax=T1)
----   29, Oct 2019 | 20:53   ----
time= 12082.0
----Position--(/earth/)-------
u=5918.74863611
z=4998767.71903
theta=360.00617721
Assiette=0.00154797756078 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.005
finesse=0.0

La capsule a donc parcouru un tour complet ($\theta=360$), on peut tracer la trajectoire à l'aide de la méthode plotraj():

In [6]:
cnap.plotraj('polar')
plt.ylim(0,20000)
0-time | 1-r | 2-theta | 3-U | 4-Gamma | 5-Pdyn | 6-phi | 7-acc | 8-Fl | 9-Alpha | 10-ref
Out[6]:
(0, 20000)

On va maintenant changer d'orbite: Calculons le ΔV nécessaire pour bifurquer vers une orbite elliptique de demi-grand axe a2=20000km ayant la terre comme foyer et calculons la nouvelle période T2:

In [7]:
a1=Rt+cnap.z
a2=20000e3
dv=cnap.u*(np.sqrt(1+(G*Mt/cnap.u**2)*((1/a1)-(1/a2)))-1)
T2=2*np.pi*np.sqrt(a2**3/(G*Mt))
print(dv)
1161.6001793

On notera ici que l'on peut accéder à la vitesse de notre capsule en faisant cnap.u.

On va maintenant réaliser la manoeuvre à l'aide de la méthode mvr() qui prend comme argument le $\Delta V$ et l'angle de la manoeuvre (ici 0):

In [8]:
cnap.mvr(dv,0) # Manoeuvre 
----   29, Oct 2019 | 20:53   ----
time= 12082.0
----Position--(/earth/)-------
u=7080.34881541
z=4998767.71903
theta=360.00617721
Assiette=0.00154797756078 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.005
finesse=0.0

Notez que la vitesse à changée. On intègre maintenant pendant une periode $T_2$ complète et on trace la nouvelle trajectoire:

In [9]:
cnap.run(dt=1,ntmax=T2) # Intégration pendant la periode orbitale T2

# On trace la trajectoire:
cnap.plotraj('polar')
plt.ylim(0,30000)
----   30, Oct 2019 | 04:42   ----
time= 40231.0
----Position--(/earth/)-------
u=7091.73695951
z=4976688.97254
theta=718.783781865
Assiette=0.372027181802 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.005
finesse=0.0

0-time | 1-r | 2-theta | 3-U | 4-Gamma | 5-Pdyn | 6-phi | 7-acc | 8-Fl | 9-Alpha | 10-ref
Out[9]:
(0, 30000)

La class Stage

On peut réaliser les simulations précédentes en utilisant un niveau supérieure de modélisation. La class Stage permet de spécifier les caractéristiques propulsives comme l'ISP, le débit massique de carburant $q_m$ et la paramètre constructif $k_s$ représentant le ratio entre la masse des structures et la masse d'ergols.

On va donc initialiser notre véhicule spatial sur la même orbite que précedemment. Pour les paramètres propulsifs, nous allons d'abord estimer la masse $\delta m$ de carburant nécessaire pour réaliser la manoeuvre précédente de $\Delta V=1161$m/s. Pour cela on utilise la formule de Tsiolkovsky vue dans le notebook précédent:

$$ \Delta V =g_0 ISP \ln{\left(\dfrac{m_i}{m_f}\right)}=g_0 ISP \ln{\left(\dfrac{m}{m-\delta m}\right)} $$

soit

$$ \dfrac{\delta m}{m_e} =(k_s+1)\left(1-e^{\frac{-\Delta V}{g_0 ISP}}\right) $$

Ainsi pour $k_s=0.24$ et $ISP=445$ (typiquement 3ème étage d'Ariane5) on a:

In [10]:
ks=0.24
ISP=445
dv=1161
(ks+1)*(1-np.exp(-dv/(9.81*ISP)))
Out[10]:
0.28956982786072571

Ainsi, en prenant 1 tonne d'ergols au départ (avant manoeuvre) on consommera donc:

In [11]:
1000*0.28956982786072571
Out[11]:
289.5698278607257

soit 289.5 kg. En prenant comme référence le débit du moteur vinci: $qm=280$kg/s, la manoeuvre durera donc:

In [12]:
tm=289.5698278607257/280
print(tm)
1.034177956645449

On peut à présent initialiser notre objet Stage et réaliser la simulation:

In [13]:
cnap=Stage(ISP=445,IC=0.24,qm=280,mc=1000,z0=zs,u0=us,gamma0=0)
print(cnap)
----   29, Oct 2019 | 17:32   ----
time= 0.0
----Position--(/earth/)-------
u=5918.22871188
z=5000000.0
theta=0.0
Assiette=0.0 deg 
Incidence=0.0deg 
---Propulsion Etage en cours---
Indice propulsif (Ft/mg)=100.48387096774194
Duree Propulsion (s): 3.5714285714285716
Moteur: Off
ISP: 445

On remarque que par défaut, le moteur est éteint. On peut donc réaliser une orbite sans rien faire comme dans la première partie:

In [14]:
cnap.run(dt=1,ntmax=T1-1)
cnap.plotraj('polar')
plt.ylim(0,20000)
----   29, Oct 2019 | 17:32   ----
time= 12081.0
----Position--(/earth/)-------
u=5918.74858582
z=4998767.87894
theta=359.976374394
Assiette=0.00154996603555 deg 
Incidence=0.0deg 
---Propulsion Etage en cours---
Indice propulsif (Ft/mg)=100.48387096774194
Duree Propulsion (s): 3.5714285714285716
Moteur: Off
ISP: 445

0-time | 1-r | 2-theta | 3-U | 4-Gamma | 5-Pdyn | 6-phi | 7-acc | 8-Fl | 9-Alpha | 10-ref
Out[14]:
(0, 20000)

Maintenant nous allons réaliser la manoeuvre en allumant le moteur pendant le temps $t_m$ calculé précedemment. Attention, la précision d'intégration est très importante ici !! Il faudra donc adapter le pas de temps pour atteindre la précision voulue:

In [15]:
cnap.engine_start()
cnap.run(dt=0.001,ntmax=tm/0.001)
cnap.engine_stop()
----   29, Oct 2019 | 17:32   ----
time= 12082.03500000021
----Position--(/earth/)-------
u=7080.95622967
z=4998767.9187
theta=360.010117473
Assiette=-0.00391028702886 deg 
Incidence=0.0deg 
---Propulsion Etage en cours---
Indice propulsif (Ft/mg)=131.13028836034127
Duree Propulsion (s): 2.5364285714286723
Moteur: On
ISP: 445

Puis on réalise une periode complète sur la nouvelle orbite:

In [16]:
cnap.run(dt=1,ntmax=T2)
----   29, Oct 2019 | 17:32   ----
time= 40231.03500000021
----Position--(/earth/)-------
u=7091.87892346
z=4977695.03198
theta=718.133895093
Assiette=0.564314510691 deg 
Incidence=0.0deg 
---Propulsion Etage en cours---
Indice propulsif (Ft/mg)=131.13028836034127
Duree Propulsion (s): 2.5364285714286723
Moteur: Off
ISP: 445

In [17]:
cnap.plotraj('polar')
plt.ylim(0,30000)
0-time | 1-r | 2-theta | 3-U | 4-Gamma | 5-Pdyn | 6-phi | 7-acc | 8-Fl | 9-Alpha | 10-ref
Out[17]:
(0, 30000)

On pourra refaire tous les calculs en changeant les caractéristiques propulsives (changer le moteur Vinci).

In [1]:
from IPython.core.display import HTML
style=open('notebooks.css', "r").read()
HTML(style)
Out[1]: