from IPython.core.display import HTML
style=open('notebooks.css', "r").read()
HTML(style)
TN n°4
La rentrée Atmosphérique
- Utilisation et documentation -
Trajectoires et paramètres importants
simon.marie@lecnam.net
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)
# 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:
- La manoeuvre de désorbitation visant à quitter l'orbite initiale.
- La partie de la trajectoire visant à rejoindre les premières couches de l'atmosphère ($z\sim120$km). (Arc orbitale)
- 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:
#! wget https://hpp.education/Lessons/MecaSpace/Solstice/SOLSTICE.pyc
On peut ensuite l'importer:
# 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:
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$.
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} $$

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$
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°)$ ?
zs=400e3
rentree(Rt+zs,200,180)
(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
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
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:
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
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
(0.0, 8000.0)
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$.

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:
| Capsule | Date | $C_d$ | $S_{ref}$ (m2) | $m$ (kg) | $\beta$ | |
|---|---|---|---|---|---|---|
![]() | Gemini | 1964-1966 | 1.5265 | 4.1547 | 1982 | $3.2$ $10^{-3}$ |
| Apollo | 1967-1974 | 1.2847 | 11.9459 | 5806 | $2.6$ $10^{-3}$ | |
![]() | Soyouz | 1967-today | 1.3254 | 3.6983 | 2850 | $1.7$ $10^{-3}$ |
![]() | Huygens | 2004 | 1.1462 | 5.7255 | 318 | $2.1$ $10^{-2}$ |
![]() | Mars Insight | 2018 | 1.2358 | 5.4739 | 358 | $1.8$ $10^{-2}$ |
| Crew Dragon | 2020 | 1.6418 | 10.7521 | 7500 | $2.35$ $10^{-3}$ | |
![]() | CST-100 | 2022 | 1.5137 | 19.7923 | 8500 | $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)
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:
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():
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
Rien de sensationnel en 1D....
Pour charger les données enregistrées pendant la simulation on utilise la méthode load():
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.
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)');
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:

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:
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
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()
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
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
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
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} $$

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
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
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
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
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
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
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
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
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
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$
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
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
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 silicium | 1900 | 0.89 | 1100 |
| niobium | 1300 | 0.89 | 310 |
| Fibre de zircone | 1400 | 0.84 | 380 |
| Tantale | 1500 | 0.64 | 360 |
En se basant sur vos résultat, estimer les températures de paroi des différentes rentrées que vous avez effectuées
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 FM5055 | 2000 | 1450 | 60000 |
| Super Lightweight Ablator SLA561V | 1500 | 260 | 49000 |
| ULD100 MDac | 1073 | 250 | 45000 |
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
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.
from IPython.display import YouTubeVideo
YouTubeVideo('AcAgnQ9K7UY',width=900,height=400)





