No description has been provided for this image No description has been provided for this image




In [2]:
from IPython.core.display import HTML
style=open('notebooks.css', "r").read()
HTML(style)
Out[2]:
-- Mécanique Spatiale --

TN n°4

No description has been provided for this image

La rentrée Atmosphérique

- Utilisation et documentation -

Trajectoires et paramètres importants

simon.marie@lecnam.net
In [3]:
import numpy as np
import scipy

import matplotlib
import matplotlib.pyplot as plt

#Option pour afficher les figures dans le notebook et eviter le plt.show():
%matplotlib inline  
plt.rc('xtick',labelsize=14)
plt.rc('ytick',labelsize=14)
In [4]:
# Constante célestes:
G=6.67408e-11 # Constante de gravitation - 
Rt=6380000 # Rayon terrestre - m
Ms=1.9884e30 # Masse du soleil - kg
Mt=5.9722e24 # Masse de la terre - kg

Présentation

L'objectif de cette séance est d'étudier la trajectoire de rentrée atmosphérique d'un corps issue d'une orbite initiale connue.

La rentrée atmosphérique peut être annalysée suivant trois phases distinctes:

  1. La manoeuvre de désorbitation visant à quitter l'orbite initiale.
  2. La partie de la trajectoire visant à rejoindre les premières couches de l'atmosphère ($z\sim120$km). (Arc orbitale)
  3. La partie purement atmosphérique souvent dimensionnante pour la mission (Arc atmosphérique). Nous distinguerons l'entrée sans portance dite balistique de l'entrée avec portance dite pilotée.

Lors de cette séance, nous ferons de nombreuses hypothèses simplificatrices pour faciliter les calculs:

  • Orbites coplanaires (2D)
  • Vitesse de rotation des corps céléstes non prise en compte.
  • Coefficients de trainée et de portance indépendant du Mach (Valide en hypersonique mais faux en supersonique).

Ces hypothèses sont celles faites dans la théorie de la rentrée atmosphérique des travaux de Allen et Eggers et permettent néanmoins de conserver les phénomènes importants intervenant dans les trajectoires de rentrée.

Enfin, il ne sera pas question dans cette séance d'étudier en détail les écoulements hypersoniques mis en jeux lors d'une rentrée atmosphérique et faisant intervenir des effets importants de dissociation moléculaire et de ionisation dues aux très fortes températures atteintes. Nous considérerons ici le flux thermique de façon simplifié comme présenté dans la théorie de Allen et Eggers:

$$ \varphi=C_q\sqrt{\rho}u^3 $$

Ou $\varphi$ est en $W.m^{-2}$ et $C_q$ est un paramètre dépendant du type d'atmosphère et du véhicule de rentrée. Pour la terre, on prendra $C_q=1.705~ 10^{-4}$ USI.

On se focalisera dans cette séance sur le calcul des trajectoires possibles en fonction des paramètres de construction des capsules de rentrée (ou atterisseurs). Les trois principales grandeurs que l'on regardera pour dimensionner une mission de rentrée seront:

  • La décélération subit par la capsule (particulièrment pour les rentrées habitées [<10g])
  • Le flux thermique maximal pour le dimensionnement de la protection thermique (masse supplémentaire à prévoir)
  • La pression dynamique maximal pour le dimensionnement des efforts subit par le lanceur pendant la rentrée.
Notes sur l'utilisation du module SOLSTICE

Afin de pouvoir experimenter les différents scénari de rentrée et de mieux comprendre l'importance de certains paramètres, vous disposez d'un "simulateur" de rentrée en important la classe Capsule du module SOLSTICE. Dans ce module, les équations du mouvement sont intégrées dans un repère polaire géocentrique. Lors de votre première utilisation il vous faut télécharger la librairie SOLSTICE:

In [ ]:
#! wget https://hpp.education/Lessons/MecaSpace/Solstice/SOLSTICE.pyc

On peut ensuite l'importer:

In [5]:
# Importation de la classe capsule
from SOLSTICE import Capsule
SOLStiCE  v4.0

Comme toutes les classes en Python, on peut consulter la documentation. En python, la fonction help(Calss) renvoi la documentation de la class ainsi que la documentation de chaque méthode de classe:

In [6]:
help(Capsule)
Help on class Capsule in module SOLSTICE:

class Capsule(Spacedevice)
 |  Capsule(beta=0.005, u=7910, z0=100000.0, gi=0.0, fin=0.0, body='Terre', name='Rentree')
 |  
 |  **Class Capsule - CNAM IAS**
 |  
 |  *Simon Marie - nov.2018*
 |  
 |  La Class Capsule permet la creation d un objet Python modelisant
 |  une capsule de rentree.  Differentes characteristiques peuvent etre saisies a l initialisation
 |  
 |  * beta : parametre balistique de la capsule (default= 5e-3) 
 |  * z0   : altitude initiale (default=1e5 m)
 |  * u    : vitesse initiale (default=7910 m/s)
 |  * gi   : angle gamma initial (default=0)
 |  * fin  : finesse de la capsule (default=0)
 |  * body : Corp principal 
 |  * moon : Corp secondaire
 |  * name : nom du fichier resultats (default=Rentree)
 |  
 |  
 |  L instanciation se fait de la facon suivante:
 |  
 |  .. code-block:: python
 |  
 |     exemple=Capsule()
 |  
 |  On pourra modifier les parametres par defaut impose a l initialisation en passant des arguments lors de la creation:
 |  
 |  .. code-block:: python
 |  
 |     exemple=Capsule(u=12000,z0=400e3)
 |  
 |  Method resolution order:
 |      Capsule
 |      Spacedevice
 |      builtins.object
 |  
 |  Methods defined here:
 |  
 |  __init__(self, beta=0.005, u=7910, z0=100000.0, gi=0.0, fin=0.0, body='Terre', name='Rentree')
 |      Initialize self.  See help(type(self)) for accurate signature.
 |  
 |  __repr__(self)
 |      caps=Capsule()
 |      Permet d afficher l instance en faisant print(caps)
 |      Peut etre desactiver avec self.display=False
 |  
 |  mvr(self, dv, omega, ref=0)
 |      Effectue une manoeuvre en imposant:
 |      
 |      * dv: differentielle de vitesse
 |      * omega: angle du delta V
 |      * ref=0 : Manoeuvre par rapport au corps principal
 |      * ref=1 : Manoeuvre par rapport au corps secondaire
 |  
 |  ----------------------------------------------------------------------
 |  Methods inherited from Spacedevice:
 |  
 |  P2S(self)
 |      Converti les coordonees par rapport au corps principal (CP)
 |      en coordonnees par rapport au corps secondaire (CS)
 |      Retours:
 |      * r2      : Position / CS
 |      * theta2  : Azimute / CS
 |  
 |  RK2(self)
 |      Integration d'ordre 2 avec le schéma Runge-Kutta
 |      met a jour tp,rp,z et theta
 |  
 |  RK4(self)
 |      Integration d'ordre 1 avec le schéma d'euler
 |      met a jour tp,rp,z et theta
 |  
 |  RK4p(self)
 |      Integration d'ordre 4 avec le schéma RK4 vect
 |      met a jour tp,rp,z et theta
 |  
 |  S2P(self, r2, theta2)
 |      Converti les coordonnées par rapport au corps secondaire (CS)
 |      en coordonnees par rapport au corps principale (CP)
 |      Retours:
 |      * r      : Position / CP
 |      * theta  : Azimute / CP
 |  
 |  Sec2days(self, seconds)
 |  
 |  anim(self, kind='polar', skip=500, ref=0)
 |      Anim la trajectoire.
 |  
 |  atmosphere(self)
 |      Trace l'evolution de la masse volumique de 
 |      l'atmosphere en fonction de l'altitude.
 |  
 |  cart2pol(self, ux, uy, theta)
 |      Coordonnées cart vers corrd polaires
 |  
 |  e2m(self)
 |      Converti les vitesses par rapport au corps principal (CP)
 |      en vitesse  par rapport au corps secondaire (CS)
 |      Retours:
 |      * u_r     : Vitesse radiale   / CS
 |      * u_theta : Vitesse angulaire / CS
 |  
 |  euler(self)
 |      Integration d'ordre 1 avec le schéma d'euler
 |      met a jour tp,rp,z et theta
 |  
 |  get_planet(self, name='Terre')
 |      Base des charactéristiques physique et orbitale des corps celestes:
 |      
 |      Soleil
 |      Mercure
 |      Venus
 |      Terre Lune
 |      Mars Phobos
 |      Ceres
 |      Jupiter Ganymede
 |      Saturne Titan
 |      Uranus Titania
 |      Neptune Triton
 |      Pluton Charon
 |      Eris Dysnomie
 |  
 |  getgravity(self)
 |      Calcul des zones d influences des differents astres locaux et met a jour la force de gravite en fonction.
 |  
 |  linkal(self)
 |  
 |  load(self)
 |      Chargement des resultats:
 |  
 |  m2e(self, ur, ut)
 |      Converti les vitesses par rapport au corps secondaire (CS)
 |      en vitesse  par rapport au corps principal (CP)
 |      Retours:
 |      * u_r     : Vitesse radiale   / CP
 |      * u_theta : Vitesse angulaire / CP
 |  
 |  monit(self, fig)
 |      Represente l evolution des donnees trajectoire.
 |      Cette methode necessite un objet matplolib (fig) en input.
 |      * 11 :U(t)    12 :Z(U)
 |      * 21 :G(Z)    22 :Phi(Z)
 |      * 31 :P(t)    32 :gamma(Z)
 |  
 |  planet(self, name='Terre')
 |      Change les caractéristiques du corps principal:
 |          self.Mt
 |          self.Rt
 |          self.rho0
 |          self.Z0
 |          self.Ml
 |          self.Rl
 |          self.al
 |          self.at
 |          self.TL
 |          self.earth_lim1
 |          self.moon_inf
 |      
 |      Corps possibles:
 |      Soleil
 |      Mercure
 |      Venus
 |      Terre
 |      Mars
 |      Jupiter
 |      Ceres
 |      Saturne
 |      Uranus
 |      Neptune
 |      Pluton
 |      Eris
 |  
 |  plotraj(self, kind='rect', impref=0)
 |      Represente la trajectoire a partir des resultat de l integration 
 |      de la trajectoire.
 |      * kind='rect' (default) => Dans le plan (x,z)
 |      * kind='polar'          => Dans le plan (r,theta)
 |  
 |  pol2cart(self, ur, ut, theta)
 |      Coordonnées polaires vers corrd cart.
 |  
 |  rho(self)
 |      Modele d atmosphere
 |  
 |  rhs(self, state, t)
 |      Terme de droite de la conservation qdm
 |  
 |  rhs_old(self, u, t)
 |      Terme de droite de la conservation qdm
 |  
 |  rhsp(self, state, statep, t)
 |      Terme de droite de la conservation qdm
 |  
 |  run(self, ntmax=1000000.0, dt=0.01)
 |      Avance la trajectoire de ntmax iterations ou jusqu a z=0
 |  
 |  setmoon(self, name='Terre')
 |      Change les caractéristiques du corps secondaire:
 |          self.Ml
 |          self.Rl
 |          self.al
 |          self.at
 |          self.TL
 |          self.earth_lim1
 |          self.moon_inf
 |      
 |      Corps possibles:
 |      Mercure
 |      Venus
 |      Terre
 |      Mars
 |      Jupiter
 |      Saturne
 |      Uranus
 |      Neptune
 |  
 |  trajint(self)
 |      Permet d integrer les equations du mouvement dans un 
 |      repere geocentrique polaire 2D (er,etheta)
 |      Le schéma d'intégration est donné par self.integ
 |      Schemas dispo:
 |      Euler
 |      RK2
 |      RK4
 |      RK4 simplectique
 |  
 |  update_udyn(self)
 |      Mise a jour de la vitesse relative de l'air
 |      Cela permet de prendre en compte la vitesse de rotation de la terre
 |      En considerant que l'atmosphere tourne avec.
 |  
 |  update_wind(self)
 |      *** Ajoute du vent et modifie l'incidence ***
 |  
 |  write(self)
 |      Ecriture de la trajectoire
 |      
 |      "0-time | 1-r | 2-theta | 3-U | 4-Gamma | 5-Pdyn | 6-phi | 7-Acc | 8-Fl
 |      | 9-Alpha | 10-ref | 11-thetaL | 12-Ft | 13-Engine | 14-r2 | 15-theta2 |
 |      16-Moon "
 |  
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from Spacedevice:
 |  
 |  __dict__
 |      dictionary for instance variables (if defined)
 |  
 |  __weakref__
 |      list of weak references to the object (if defined)
 |  
 |  ----------------------------------------------------------------------
 |  Data and other attributes inherited from Spacedevice:
 |  
 |  G = 6.67408e-11
 |  
 |  Ms = 1.9884e+30
 |  
 |  thetaL = 1.5707963267948966

I - Désorbitation

La première phase de la rentrée consiste à sortir de l'orbite de "stationnement" pour bifurquer sur une orbite elliptique sécante à la surface d'"atterissage".

Pour cela, on procède à une manoeuvre orbitale en injectant de la quantité de mouvement $\Delta U$ suivant un angle choisit $\omega$.

No description has been provided for this image

Avec les notations suivantes:

$$ \delta=\dfrac{\Delta U}{U_s} ~~~ \xi=\dfrac{r_0}{r_e} ~~~ k_0=\dfrac{U_0}{U_s} ~~~ k_e=\dfrac{U_e}{U_s} $$

On montre, en utilisant les relations (Al-Kashi) dans le triangle $(U_0,U_s,\Delta U)$ que:

$$ U_0^2=U_s^2+\Delta U^2-2U_s\Delta U \cos(\pi-\omega) $$

soit

$$ k_0=\sqrt{1+\delta^2+2\delta\cos(\omega)}\\ \sin(\gamma_0)=\dfrac{\delta \sin(\omega)}{k_0} $$

No description has been provided for this image

En pratique, on considère une orbite initiale circulaire avec $U_s=\sqrt{\dfrac{GM_t}{r_s}}=\sqrt{\dfrac{GM_t}{r_0}}$.

Ainsi, la conservation de l'énergie sur l'arc orbital donne:

$$ \dfrac{1}{2}U_e^2-\dfrac{GM}{r_e}=\dfrac{1}{2}U_0^2-\dfrac{GM}{r_0} $$

soit:

$$ U_e=U_s\sqrt{2\xi-1+\delta^2+2\delta\cos(\omega)} $$

De la même façon, en utilisant la conservation du moment cinétique $r\dot{\theta}$ sur l'arc orbital soit:

$$ r_e U_e\cos(\gamma_e)=r_0 U_0\cos(\gamma_0) $$

on peut exprimer l'angle de rentrée $\gamma_e$ en fonction des paramètres de manoeuvre $\delta$ et $\omega$ et du paramètre initial $\xi$:

$$ \cos(\gamma_e)=\xi\dfrac{1+\delta\cos(\omega)}{\sqrt{2\xi-1+\delta^2+2\delta\cos(\omega)}} $$

A partir des résultats précédents, créer une fonction permattant de trouver $U_e$ et $\gamma_e$ en fonction de $\Delta U$, $\omega$ et $r_0$

In [12]:
def rentree(r0,du,w):
    us=np.sqrt(G*Mt/r0)
    delta=du/us
    xi=r0/(Rt+120e3)
    rad=2*xi-1+delta**2+2*delta*np.cos(w*np.pi/180)
    ue=us*np.sqrt(rad)
    ge=np.arccos(xi*(1+delta*np.cos(w*np.pi/180))/np.sqrt(rad))
    return ue,ge*180/np.pi

Quel est l'angle et la vitesse de rentrée d'une capsule de rentrée quittant une orbite circulaire de 400km d'altitude avec une manoeuvre du type $(\Delta U, \omega)=(200,180°)$ ?

In [13]:
zs=400e3
rentree(Rt+zs,200,180) 
Out[13]:
(7799.162572274184, 2.9150199313235436)

Ainsi, il est possible de déterminer les paramètres d'entrée $(U_e,\gamma_e)$ à partir d'une condition initiale ($U_s$,$r_0$) et des paramètres de manoeuvre ($\Delta U$,$\omega$).

A l'aide du module SOLSTICE fourni, placer une capsule sur une orbite circulaire à 400km d'altitude

In [19]:
zs=400e3
us=np.sqrt(G*Mt/(Rt+zs))
orb=Capsule(z0=zs,u=us)
print(orb)
----   05, Nov 2024 | 09:44   ----
time= 0.0
----Position--(/Terre/)-------
u=7667.398378831312
z=400000.0
theta=0.0
Assiette=0.0 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.005
finesse=0.0

Calculer la position de la capsule au bout d'1 heure

In [20]:
dt=1.
orb.run(dt=dt,ntmax=60*60/dt)
----   05, Nov 2024 | 10:44   ----
time= 3601.0
----Position--(/Terre/)-------
u=7667.292369075144
z=400073.4224566519
theta=233.32301942439943
Assiette=0.0009743813070319577 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.005
finesse=0.0

Effectuer une manoeuvre de $\Delta u=200$m/s avec $\omega=180°$ puis calculer la nouvelle position de la capsule, 20 minutes après la manoeuvre:

In [21]:
orb.mvr(200,180)
orb.run(dt=1,ntmax=20*60)
----   05, Nov 2024 | 10:44   ----
time= 3601.0
----Position--(/Terre/)-------
u=7467.292369075144
z=400073.4224566519
theta=233.32301942439943
Assiette=0.0009743813070321456 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.005
finesse=0.0

----   05, Nov 2024 | 11:04   ----
time= 4802.0
----Position--(/Terre/)-------
u=7799.260728414549
z=119832.14040117431
theta=311.4120761660617
Assiette=2.916449858001694 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.005
finesse=0.0

In [22]:
orb.plotraj('polar')
plt.ylim(0,8000)
0-time | 1-r | 2-theta | 3-U | 4-Gamma | 5-Pdyn | 6-phi | 7-Acc | 8-Fl | 9-Alpha | 10-ref | 11-thetaL | 12-Ft | 13-Engine | 14-r2 | 15-theta2 | 16-Moon 
Out[22]:
(0.0, 8000.0)
No description has been provided for this image

Comparer les résultats à ceux calculés par la fonction rentree()

Nous allons à présent étudier l'importance des paramètres d'entrée $(U_e,\gamma_e)$ sur la trajectoire et le dimensionnement mécanique de notre capsule de rentrée.

II - Une question de freinage

Considérons tout d'abord le cas le plus simple, en partant d'un corps en chute libre laché d'une altitude de 120km avec une vitesse initiale $U_s$ sous un angle de $90°$ avec l'horizontal local. Ce corp est soumis à son poid $P$ et à une force de trainée $F_d$.

No description has been provided for this image

Le bilan s'écrit:

$$ m\dfrac{dv}{dt}=P-F_d=mg_0-\dfrac{1}{2}\rho(z)v(z)^2*S C_d $$

en définissant le coefficient balistique $\beta=\dfrac{S.C_d}{m}$, le bilan devient:

$$ \dfrac{dv}{dt}=\dfrac{mg_0-F_d}{m}=g_0-\dfrac{1}{2}\beta \rho(z)v(z)^2 $$

Ainsi la décélération subi par le corps est directement lié à la force de trainée:

$$ \Gamma=g_0-\dfrac{1}{2}\rho(z)v(z)^2\beta $$

De la même façon que pour l'étape du décollage, on se donnera un modèle d'atmosphère du type:

$$ \rho(z)=\rho_0 e^{-z/z_0} $$

Pour la terre on prendra les valeurs proposées par la théorie de Allen et Eggers $Z_0=6700$m et $\rho_0=1.735$ $kg.m^{-3}$.

Remarque: Il est important de remarquer que les modèles d'atmosphère pour la rentrée et pour le départ, ne sont pas nécessairement les mêmes. En effet, pour la rentrée atmosphérique, on a généralement besoin de décrire l'atmosphère "par le haut" avec des vitesses d'approches très importantes.

En utilisant les hypothèses précédentes on peut facilement simplifier les équations du mouvement et montrer que la décélération maximale vaut:

$$ \Gamma_{max}=\dfrac{U_0^2 \sin(\gamma_0)}{2eZ_0} $$

Elle ne dépend donc pas de la forme de la capsule !! En revanche cette décélération maximale est atteinte pour:

$$ Z_g=Z_0\ln\left(\dfrac{\rho_0 Z_0 \beta }{\sin(\gamma_0)}\right) $$

Qui dépend donc de la forme de la capsule à travers $\beta$.

Remarque: Le coefficient balistique est parfois noté $\zeta=\dfrac{1}{\beta}=\dfrac{m}{S C_d}$ et s'exprime en $kg/m^2$. Ici nous utiliserons $\beta$.

Nous allons essayer de retrouver ce résultat par des simulations.

Le tableau suivant donne quelques exemples de coefficients $\beta$ pour des attérisseurs célèbres:

CapsuleDate$C_d$$S_{ref}$ (m2)$m$ (kg)$\beta$
No description has been provided for this imageGemini1964-19661.52654.15471982$3.2$ $10^{-3}$
No description has been provided for this imageApollo1967-19741.284711.94595806$2.6$ $10^{-3}$
No description has been provided for this imageSoyouz1967-today1.32543.69832850$1.7$ $10^{-3}$
No description has been provided for this imageHuygens20041.14625.7255318$2.1$ $10^{-2}$
No description has been provided for this imageMars Insight20181.23585.4739358$1.8$ $10^{-2}$
No description has been provided for this imageCrew Dragon20201.641810.75217500$2.35$ $10^{-3}$
No description has been provided for this imageCST-10020221.513719.79238500$3.5$ $10^{-3}$

A l'aide du module SOLSTICE placer une capsule sur une orbite initiale de 120km d'altitude avec un angle de 90° (à la vertical)

In [46]:
gem=Capsule(gi=90,z0=120e3)
print(gem)
----   05, Nov 2024 | 10:54   ----
time= 0.0
----Position--(/Terre/)-------
u=7910
z=120000.0
theta=0.0
Assiette=90.0 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.005
finesse=0.0

Par défaut la création de l'instance gem affiche les coordonées de la capsule: $z_0,\theta_0$, son orientation $\gamma$ et sa vitesse initiale $u$.

On peut alors calculer la trajectoire de la capsule en partant de ces conditions initiales. C'est la méthode run() qui permet de faire cela. Cette méthode réalise l'intégration des équations du mouvement.

Lançons la capsule dans l'atmosphère:

In [47]:
gem.run()
----   05, Nov 2024 | 10:57   ----
time= 196.74999999996604
----Position--(/Terre/)-------
u=47.93468608469881
z=-0.3048667050898075
theta=4.037800905044168e-07
Assiette=89.99980547411715 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.005
finesse=0.0

En pratique l'intégration s'arrête lorsque la capsule a atteint le sol (z=0) ou que le nombre d'itérations a atteint une limite fixée (par défaut $10^6$).

Quelle est la vitesse de la capsule lors de l'impact au sol ? Quelle est la durée de la rentrée ?

On peut représenter la trajectoire avec la méthode plotraj():

In [48]:
gem.plotraj()
0-time | 1-r | 2-theta | 3-U | 4-Gamma | 5-Pdyn | 6-phi | 7-Acc | 8-Fl | 9-Alpha | 10-ref | 11-thetaL | 12-Ft | 13-Engine | 14-r2 | 15-theta2 | 16-Moon 
No description has been provided for this image

Rien de sensationnel en 1D....

Pour charger les données enregistrées pendant la simulation on utilise la méthode load():

In [49]:
data=gem.load() # On charge les données enregistrtées pendant le "vol"
0-time | 1-r | 2-theta | 3-U | 4-Gamma | 5-Pdyn | 6-phi | 7-Acc | 8-Fl | 9-Alpha | 10-ref | 11-thetaL | 12-Ft | 13-Engine | 14-r2 | 15-theta2 | 16-Moon 

La méthode nous renvoi la structure du tableau de données "data". Ainsi data[:,0] contient le temps, data[:,1] l'altitude...

Tracer l'évolution de la décélération en fonction de l'altitude et comparer la valeur max à celle estimer plus haut.

In [50]:
fig=plt.figure(figsize=(10,6),dpi=120)
# On calcul la deceleration maximal et son altitude
Z0,u0=6700,7910
rho0=1.735
Gm=u0**2*np.sin(90*np.pi/180)/(2*np.exp(1)*gem.Z0) # Gamma max
Zg=gem.Z0*np.log(gem.rho0*gem.Z0*gem.beta/np.sin(90*np.pi/180)) # Zg
plt.plot(Gm/9.81,Zg/1000,'sr')

# Comparons avec notre simulation:
plt.plot(data[:,7]/9.81,(data[:,1]-Rt)/1000);
plt.ylabel('Alt (km)');plt.xlabel('Deceleration (g)');
No description has been provided for this image

Qu'observez-vous ?

Regardons alors l'importance de l'angle $\gamma_e$ dans une trajectoire plus réaliste en 2D.

III - Rentrée balistique (sans portance)

Ici nous considérons que la rentrée se fait sous un angle initial $\gamma_e$. La capsule est soumise aux deux mêmes forces que précédemment mais avec une angle $\gamma$ quelconque:

No description has been provided for this image

Les équations du mouvement dans le plan polaire $(r,\theta)$ s'écrivent:

$$ m(\ddot{r}-r\dot{\theta}^2) = -P+F_d\sin(\gamma_e)\\ m(r\ddot{\theta}+2\dot{r}\dot{\theta}) = -F_d\cos(\gamma_e) $$

soit en divisant par la masse et en utilisant les notations introduites précédemment:

$$ \ddot{r}-r\dot{\theta}^2 = -\dfrac{\mathcal{G}M_t}{(R_t+r)^2}+\dfrac{1}{2}\rho(r)U^2\beta\sin(\gamma_e)\\ r\ddot{\theta}+2\dot{r}\dot{\theta} = -\dfrac{1}{2}\rho(r)U^2\beta\cos(\gamma_e) $$

Ces équations sont intégrées dans la classe Capsule.

A l'aide du module SOLSTICE, placer une capsule sur la même orbite que précédemment en changeant uniquement l'angle de rentré à $\gamma_e=6$, puis réaliser le calcul de la trajectoire:

In [55]:
gem=Capsule(gi=3,z0=120e3)
print(gem)
----   05, Nov 2024 | 11:02   ----
time= 0.0
----Position--(/Terre/)-------
u=7910
z=120000.0
theta=0.0
Assiette=3.0 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.005
finesse=0.0

In [56]:
gem.run()
----   05, Nov 2024 | 11:11   ----
time= 494.85999999969494
----Position--(/Terre/)-------
u=47.93411861604016
z=-0.4555568331852555
theta=13.948180265153185
Assiette=89.99980886471765 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.005
finesse=0.0

Quels est la vitesse de la capsule lors de l'impact au sol ? Quelle est la durée de la rentrée ? Tracer la trajectoire à l'aide de la méthode plotraj()

In [57]:
fig=plt.figure(figsize=(15,5),dpi=120)
gem.plotraj()
0-time | 1-r | 2-theta | 3-U | 4-Gamma | 5-Pdyn | 6-phi | 7-Acc | 8-Fl | 9-Alpha | 10-ref | 11-thetaL | 12-Ft | 13-Engine | 14-r2 | 15-theta2 | 16-Moon 
No description has been provided for this image

Que remarquez-vous ?

A l'aide de la méthode monit(fig) commenter l'évolution de la décélération et du flux thermique absorbé par la capsule

In [58]:
fig=plt.figure(figsize=(12,14),dpi=90)
gem.monit(fig)
0-time | 1-r | 2-theta | 3-U | 4-Gamma | 5-Pdyn | 6-phi | 7-Acc | 8-Fl | 9-Alpha | 10-ref | 11-thetaL | 12-Ft | 13-Engine | 14-r2 | 15-theta2 | 16-Moon 
No description has been provided for this image

Conclusion ?

IV - Rentrée pilotée (avec portance)

Les équations du mouvement a intégrer sont sensiblement les mêmes que dans le cas précédent, en ajoutant la force de portance. Le coefficient de portance est supposé constant et proportionel au coefficient de trainée:

$$ finesse=\dfrac{C_L}{C_D} $$

No description has been provided for this image

La finesse d'un vehicule de rentrée est plutôt faible ($\sim0.3$) et pouvait atteindre l'unité pour les (feu-)navettes américaines.

En suivant la même procédure que précédemment, étudier la rentrée de la même capsule ayant une finesse de 0.3

In [59]:
cap=Capsule(gi=3,z0=120e3,fin=0.3,name='Rentree_port.dat')
print(cap)
----   05, Nov 2024 | 11:12   ----
time= 0.0
----Position--(/Terre/)-------
u=7910
z=120000.0
theta=0.0
Assiette=3.0 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.005
finesse=0.3

In [60]:
cap.run(dt=0.1)
----   05, Nov 2024 | 11:25   ----
time= 806.3000000001148
----Position--(/Terre/)-------
u=46.85338003589934
z=-4.086200397461653
theta=23.957759124820384
Assiette=73.02402167940645 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.005
finesse=0.3

In [61]:
fig=plt.figure(figsize=(15,5),dpi=90)
cap.plotraj()
0-time | 1-r | 2-theta | 3-U | 4-Gamma | 5-Pdyn | 6-phi | 7-Acc | 8-Fl | 9-Alpha | 10-ref | 11-thetaL | 12-Ft | 13-Engine | 14-r2 | 15-theta2 | 16-Moon 
No description has been provided for this image
In [62]:
fig=plt.figure(figsize=(12,14),dpi=90)
cap.monit(fig)
gem.monit(fig)
0-time | 1-r | 2-theta | 3-U | 4-Gamma | 5-Pdyn | 6-phi | 7-Acc | 8-Fl | 9-Alpha | 10-ref | 11-thetaL | 12-Ft | 13-Engine | 14-r2 | 15-theta2 | 16-Moon 
0-time | 1-r | 2-theta | 3-U | 4-Gamma | 5-Pdyn | 6-phi | 7-Acc | 8-Fl | 9-Alpha | 10-ref | 11-thetaL | 12-Ft | 13-Engine | 14-r2 | 15-theta2 | 16-Moon 
No description has been provided for this image

Comparer alors les deux trajectoires. Que pouvez-vous dire ?

Maintenant, réaliser la même simulation en choisissant un angle de rentrée très faible de l'ordre du demi-degrès

In [63]:
cap2=Capsule(gi=0.5,z0=120e3,fin=0.3,name='Rentree_port2.dat')
cap2.run(dt=0.1)
----   05, Nov 2024 | 17:10   ----
time= 20480.59999998246
----Position--(/Terre/)-------
u=46.86339528015038
z=-1.291457887738943
theta=1359.9208022550126
Assiette=73.02410181325249 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.005
finesse=0.3

In [64]:
fig=plt.figure(figsize=(12,4),dpi=90)
cap2.plotraj()
0-time | 1-r | 2-theta | 3-U | 4-Gamma | 5-Pdyn | 6-phi | 7-Acc | 8-Fl | 9-Alpha | 10-ref | 11-thetaL | 12-Ft | 13-Engine | 14-r2 | 15-theta2 | 16-Moon 
No description has been provided for this image
In [65]:
fig=plt.figure(figsize=(12,18),dpi=90)
cap2.monit(fig)
0-time | 1-r | 2-theta | 3-U | 4-Gamma | 5-Pdyn | 6-phi | 7-Acc | 8-Fl | 9-Alpha | 10-ref | 11-thetaL | 12-Ft | 13-Engine | 14-r2 | 15-theta2 | 16-Moon 
No description has been provided for this image

Qu'observez-vous ?

V - Mission Complète

Maintenant vous devez pouvoir revenir sur terre en un point précis à partir d'une orbite donnée:

En partant d'une capsule de type "Orion" de finesse $0.26$, initialement sur une orbite circulaire de 350 km d'altitude, vous devez atterir au point $\theta = 272°$ avec les contraintes suivantes: $\Gamma<10g$ et $\varphi<1000 kW/m^2$, précision$<100m$

In [83]:
cap3=Capsule(gi=0,z0=350e3,beta=2.35e-3,fin=0.26)
#cap3.run(dt=0.1,ntmax=60*60*5.65)
cap3.run(dt=0.1,ntmax=1)
cap3.mvr(200,210)
cap3.run(dt=0.1)
cap3.plotraj()
print("precision= "+str(Rt*(((180*cap3.theta/np.pi)-272)*np.pi/180))+" m")
----   05, Nov 2024 | 11:43   ----
time= 0.2
----Position--(/Terre/)-------
u=7909.999995597803
z=350000.00993262324
theta=0.013468339247935265
Assiette=-0.0007194652798927814 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.00235
finesse=0.26

----   05, Nov 2024 | 11:43   ----
time= 0.2
----Position--(/Terre/)-------
u=7737.441150296928
z=350000.00993262324
theta=0.013468339247935265
Assiette=-0.7412404408671102 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.00235
finesse=0.26

----   06, Nov 2024 | 15:30   ----
time= 100000.3000013329
----Position--(/Terre/)-------
u=7781.655443726914
z=309487.95656148996
theta=6447.161983373231
Assiette=-0.2844004272965388 deg 
Incidence=0.0deg 
--------Parametres-----------
beta=0.00235
finesse=0.26

0-time | 1-r | 2-theta | 3-U | 4-Gamma | 5-Pdyn | 6-phi | 7-Acc | 8-Fl | 9-Alpha | 10-ref | 11-thetaL | 12-Ft | 13-Engine | 14-r2 | 15-theta2 | 16-Moon 
precision= 687616675.935539 m
No description has been provided for this image
In [84]:
fig=plt.figure(figsize=(16,16))
cap3.monit(fig)
0-time | 1-r | 2-theta | 3-U | 4-Gamma | 5-Pdyn | 6-phi | 7-Acc | 8-Fl | 9-Alpha | 10-ref | 11-thetaL | 12-Ft | 13-Engine | 14-r2 | 15-theta2 | 16-Moon 
No description has been provided for this image

Y'a t'il unicité de la solution ? Si non comment pourrait-on calculer la manoeuvre optimale nécéssitant le moins de carburant ?

VI - Quelques notions sur les dispositifs de freinages

Comme nous l'avons vu, le flux thermique auquel est soumis le véhicule de rentrée peut rester important même avec des angles de rentrée relativement faibles ($500-800$ $kW/m^2$). Ainsi la protection thermique devient indispensable pour assurer l'intégrité du véhicule.

Deux solutions sont possibles:

  • les protections rayonnantes pour les rentrées plannées à faibles $\gamma_e$
  • les protections ablatives pour les rentrées balistiques à $\gamma_e$ plus élevé engeandrant des flux thermique plus importants.

1 - Protections thermiques Rayonnantes

Pour ce type de protection, la paroi externe du bouclier emet un flux thermique par rayonement qui est en équilibre thermique avec le flux convectif incident. D'après la loi de Stephan-Boltzmann, on a donc en notant $T_p$ la température de paroi:

$$ \varphi_{i}=\varphi_{rayonement}=\epsilon\sigma T_p^4 $$

soit:

$$ T_p=\left(\dfrac{\varphi_i}{\epsilon\sigma}\right)^{(1/4)} $$

$\sigma=5.670367~ 10^{-8}$ $W.m^{-2}.K^{-4}$ $\epsilon$: émissivité du materiaux

Voici quelques candidats pour la construction de ce type de protection

Materiaux$T_{max}$ °C$\epsilon$$\sim \varphi_{max}$ ($kW/m^2$)
Carbure de silicium19000.891100
niobium13000.89310
Fibre de zircone14000.84380
Tantale15000.64360

En se basant sur vos résultat, estimer les températures de paroi des différentes rentrées que vous avez effectuées

In [ ]:
traj_data=cap.load()
phi=traj_data[:,6]
Tp=...

2 - Protections thermiques Ablatives

Ce deuxième type de protection est basée sur le principe de l'absorbtion de la chaleur par des réactions endothermiques de type changement d'état (fusion, sublimation, vaporisation...).

Dans ce cas, la surface subissant le flux thermique reste à température constante $T_a$ appellée température d'ablation. On peut écrire un modèle très simplifié d'évolution du front d'ablation $s$:

$$ \varphi(t)\sim\rho_a H_a \dfrac{\partial s}{\partial t} $$

Materiaux$T_{a}$ (°C)$\rho$ (kg/$m^3$)$H_a$ (kJ/kg)
Carbon Phenolic FM50552000145060000
Super Lightweight Ablator SLA561V150026049000
ULD100 MDac107325045000

On peut donc estimer globalement:

$$ \mathcal{e}\sim\int_{t_0}^{t_1}\dfrac{\varphi(t)}{\rho_a H_a}dt $$

A partir des résultats de la dernière simulation, estimer la taille et la masse du bouclier thermique a utiliser pour absorber le flux thermique incident

In [ ]:
e=...
print('ep= ',e)

Conculsion

3 - Parachute supersonique

Lorsque le nombre de Mach local a atteint des niveaux supersoniques raisonable ($\sim 1.5$), on peut utiliser un système de freinage basé sur le déploiement d'un parachute supersonique. En pratique la décélération de basse atmosphère dépend énormément de la mission et du type d'atmosphère. Sur ce sujet, on pourra consulter les références données en conclusion.

In [1]:
from IPython.display import YouTubeVideo
YouTubeVideo('AcAgnQ9K7UY',width=900,height=400)
Out[1]:

Conclusion¶

Rédigez une synthèse de toutes les notions que vous avez parcourues dans ce projet.

Références¶

Voici quelques références qui ont servi a mettre au point ce notebook. Les lecteurs intéressés pourront approfondir de nombreux sujets en les parcourant.

  • Allen & Eggers - 1953
  • Systèmes Spatiaux - D.Marty - 1994
  • Coming Home (NASA)
  • Undocking from ISS (ESA)
  • Atmospheric rentry (WIKI)
  • LDSD (WIKI)
No description has been provided for this image No description has been provided for this image
Retour en haut de la page