Programme de calcul de vol d'une fusée

Le problème exposé dans ce sujet a été résolu.

Bonjour,

Nous avons eu en classe un TP dans lequel nous devions réaliser le calcul de vol d’une fusée en phase de montée en Python. On nous a donnée le système différentiel suivant :

{dvdt=Dumgk0exp(za0)×v2mdmdt=Ddzdt=v\begin{cases} \frac{dv}{dt} &= \frac{Du}{m} - g - k_0 \exp\left(-\frac{z}{a_0}\right) \times \frac{v^2}{m}\\ \frac{dm}{dt} &= -D\\ \frac{dz}{dt} &= v \end{cases}

Avec les données suivantes :

  • D=4 kg.s1D = 4\ kg.s^{-1}, le débit massique des gaz
  • a0=8.103 ma_0 = 8.10^3\ m
  • g=9.81 m.s2g = 9.81\ m.s^{-2}
  • k0=0.1 N.m2s2k_0 = 0.1\ N.m^{-2}s^2,
  • u=2.10³ m.s1u = 2.10³\ m.s^{-1}
  • vv, la vitesse, initialement à 0 m.s10\ m.s^{-1}
  • mm la masse initialement à 400 kg400\ kg et à 80 kg80\ kg à la fin de la combustion
  • zz l’altitude, initialement à 0 m0\ m

Je souhaiterais également réaliser la phase de descente de la fusée, il faut donc modifier notre système différentiel (notamment supprimer la poussée). J’ai réalisé le code suivant :

from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math


D = 4
g = 9.81
a0 = 8e3
k0 = 0.1
u = 2e3

# Y[0] => vitesse : v
# Y[1] => masse : m
# Y[2] => altitude : z


def subdivise(a, b, N):
   return [a + k * ((b - a) / (N - 1)) for k in range(N)]


def compute_dv(Y, t):
   if Y[1] > 80:  # Si la masse est supérieure à 80
       dv = ((D * u) / Y[1]) - g - k0 * math.exp(-(Y[2] / a0)) * (Y[0] ** 2 / Y[1])
       dm = -D
       dz = Y[0]
   else:
       dv = -g - k0 * math.exp(-(Y[2] / a0)) * (Y[0] ** 2 / Y[1])
       dm = 0
       dz = Y[0]
   
   if Y[2] < 0:  # Si l'altitude est inférieure à 0
       dv = 0
       dz = 0
       dm = 0

   return dv, dm, dz


if __name__ == "__main__":
   times = subdivise(0, 300, 250)
   Y0 = (0, 400, 0)

   Ye = odeint(compute_dv, Y0, times)

   velocity = [y[0] for y in Ye]
   weight = [y[1] for y in Ye]
   altitude = [y[2] for y in Ye]
   
   plt.subplot(2, 2, 1)
   plt.plot(times, velocity, label="vitesse")
   plt.legend()

   plt.subplot(2, 2, 2)
   plt.plot(times, weight, label="masse")
   plt.legend()

   plt.subplot(2, 1, 2)
   plt.plot(times, altitude, label="altitude")
   plt.legend()

   plt.show()

Ce qui me donne les courbes suivantes :

Vol d'une fusée
Vol d'une fusée

On peut voir que la vitesse donne n’importe quoi et je n’arrive pas à savoir pourquoi. Auriez-vous une idée ?

Merci pour votre aide !

Salut,

Le solveur te donne un indice, il donne ce warning :

lsoda--  warning..internal t (=r1) and h (=r2) are
      such that in the machine, t + h = t on the next step           
      (h = step size). solver will continue anyway
     in above,  r1 =  0.2657747190447D+03   r2 =  0.2339161631143D-13

Il dit que la machine n’est pas capable de différencier tt de t+δtt+\delta t parce que δt\delta t est relativement trop petit. Ça vient, je pense, de ta condition qui passe toutes les dérivées à 00 lorsque tu arrives sur le sol. Remarque que ton code fonctionne très bien si tu t’arrêtes plus tôt (j’ai pris times = subdivise(0, 260, 100)) :

trajectoire avant de toucher le sol
trajectoire avant de toucher le sol

EDIT : en fait, cela vient de ton accélération qui tend vers l’infini, probablement à cause du v2v^2 dans le terme de frottements. Il y a quelque chose qui cloche dans tes équations… On voit sur le plot de la vitesse que l’accélération augmente furieusement sur la fin.

+2 -0

Salut,

J’ai aussi des soupçons sur cette ligne-là :

dv = -g - k0 * math.exp(-(Y[2] / a0)) * (Y[0] ** 2 / Y[1])

La force de frottement représentée par le deuxième terme doit être opposée à la vitesse. Or, l’équation donnée dans ton énoncé fait implicitement l’hypothèse que la vitesse est positive. L’équation donnée est donc fausse pour la phase de vitesse négative qui arrive quand on a perdu tout l’élan et qu’on tombe.

Il est assez facile de prendre en compte le signe de la vitesse. Il y a peut-être d’autre subtilités, mais on voit déjà que c’est plus normal. Le résultat :

Avec correction de l'équation pour la rendre valide aussi pour les vitesses négatives.
Avec correction de l'équation pour la rendre valide aussi pour les vitesses négatives.

Si on zoome, on voit à quelle vitesse on touche le sol. Ouille !

Avec ma modification, on n’a pas d’avertissement du solveur.

Pour ceux qui comme moi se sont interrogé sur la diminution de la vitesse (en valeur absolue) à la fin de la chute, je pense que ça simule la rentrée dans l’atmosphère.

Une fois les moteurs coupés, la vitesse de la fusée va diminuer, passer par 0 et atteindre une vitesse limite due aux frottements si ceux-ci ne dépendent que de la vitesse, donc jusqu’au choc, la vitesse devrait augmenter (en valeur absolue) jusqu’à un plateau. Là, on voit que la vitesse augmente bien dans un premier temps (en absolue), mais diminue ensuite ! Si je suis bien, ça vient du terme en zz dans l’équation : les frottements ne dépendent pas que de la vitesse ici. Les frottement diminuent à haute altitude, on simule ici la rentrée dans l’atmosphère (moins d’air, donc moins de frottement à haute altitude). On retrouve d’ailleurs ce a0a_0 à partir de la formule de la pression en fonction de l’altitude. gg, lui, ne change pas, car 60 km, c’est peu devant les 6 000 km du rayon de la Terre.

Je n’ai souvenir d’avoir fait ça ; c’est amusant, en tout cas.

+1 -0

Merci pour votre aide ! En effet, les frottements sont toujours opposés à la vitesse, je n’hésiterais pas à le noter dans mon rapport !

fusee.png
fusee.png

Je me posais la même question que @Gabbro et j’en suis venu aux mêmes conclusions.. le point jaune correspond bien à une rentrée dans une atmosphère plus épaisse ?

Le point jaune correspond bien sûr à la rentrée dans l’atmosphère (je ne comprends pas ce que tu veux dire par plus épaisse ?). Ce qui se passe, c’est que lorsque l’altitude diminue jusqu’à être de l’ordre de a0a_0, le terme exponentiel devient important et force la vitesse à se caler sur une valeur plus faible (évidemment, si l’atmosphère était plus épaisse, ie a0a_0 plus grande, cette vitesse terminale pourrait potentiellement être plus grande que la vitesse d’entrée).

L’expression exponentielle vient du fait que la densité atmosphérique est grosso-modo elle-même une fonction exponentielle de l’altitude.

J’ai dit plus épaisse car notre fusée ne dépasse pas les 100 km et est très loin des 800 km. Par conséquent elle est toujours dans l’atmosphère, du coup je trouve le terme rentrée dans l’atmosphère injustifié. ;)

Merci pour tes explications en tout cas !

Ah, tu veux dire plus dense, alors, pas plus épaisse. En fait, la notion de rentrée de l’atmosphère a du sens quand même, précisément parce qu’elle est associée, dans ce contexte, à une forte augmentation des frottements. La notion d’épaisseur d’une atmosphère est nécessairement arbitraire, l’atmosphère d’une planète n’a pas de borne supérieure absolue (celle de la Terre va allègrement taquiner la Lune par exemple). Donc ici, ce qu’on prend comme épaisseur de l’atmosphère est de l’ordre de a0a_0 (par exemple la position de ton trait jaune, à environ 4a04a_0), l’épaisseur caractéristique au-delà de laquelle les frottements sont faibles devant les autres forces en jeu.

Une question intéressante peut être : est-ce que lors d’une rentrée atmosphérique, le point bas d’inversion de la courbe des vitesses correspond-il à max Q ? En tous cas, tu peux au moins indiquer que tes pointillés jaunes correspondent au point où l’augmentation des frottements (dûs à l’augmentation de la vitesse et de la densité atmosphérique) compensent l’accélération dûe à la gravité, et donc où ta fusée ralentit au lieu de continuer à accélérer. Ce sera sans doute plus clair si tu traces en plus la courbe d’accélération.

PS : on entend parler souvent d’atmosphère qui s’arrête à 100 km ou d’espace qui commence à 100 km (avec des variantes à 80 km ou 60 miles). C’est en fait à cause de la ligne de Kármán (et sa variante américaine) qui définit la limite officielle internationale (et aux USA) de l’espace. C’est pas un choix au hasard. Plus on va haut, moins l’atmosphère est dense, et plus on doit aller vite pour qu’un aéronef soit pilotable à l’aide de surface aérodynamiques. Arrive une altitude à laquelle cette vitesse est à peu près la même que celle de satellisation à la même altitude : c’est un point à partir duquel on ne peut plus piloter d’aéronef via des surfaces aérodynamiques aussi grandes soient-elles. On a arrondit tout ça à la louche pour donner un nombre rond facile à retenir, et voilà.

Il y a un projet de recalcul de cette altitude pour la définition internationale, qui viendrait probablement l’abaisser.

PPS : l’ISS comme presque tous les satellites en orbite basse ne sont même pas dans la couche la plus haute de notre atmosphère.

PPPS : en atmosphère terrestre, les freinages de rentrées atmosphériques depuis l’orbite commencent vers 80 km d’altitude et atteignent leurs maximums d’intensité vers 40 km. C’est différent avec ta fusée parce que n’étant pas orbitale, elle est beaucoup plus lente.

Connectez-vous pour pouvoir poster un message.
Connexion

Pas encore membre ?

Créez un compte en une minute pour profiter pleinement de toutes les fonctionnalités de Zeste de Savoir. Ici, tout est gratuit et sans publicité.
Créer un compte